用Python搞定时空数据分析：手把手教你用mgtwr包跑通GTWR模型（附完整代码）

阿丁的猫

296人浏览 · 2026-06-06 15:11:08

阿丁的猫 · 2026-06-06 15:11:08 发布

Python实战：用GTWR模型解锁时空数据分析的完整流程

时空数据分析正在成为城市规划、环境科学和商业决策中的关键工具。想象一下，你手头有一组包含地理位置和时间戳的房价数据，如何从中挖掘出空间位置和时间变化对价格的综合影响？这正是地理时空加权回归（GTWR）模型的用武之地。与传统的GWR模型相比，GTWR不仅考虑空间异质性，还引入了时间维度，能够更精确地捕捉变量关系的时空演变特征。

本文将带你从零开始，使用Python中的mgtwr包完成GTWR建模全流程。不同于理论讲解，我们聚焦于 实战操作 ，涵盖数据准备、参数优化、模型拟合到结果可视化的每个技术细节。无论你是地理信息专业的研究生，还是需要处理时空数据的分析师，这份指南都能让你快速获得可复用的代码解决方案。

1. 环境配置与数据准备

1.1 安装必要工具包

GTWR分析需要一组专门的地理统计工具链。推荐使用Anaconda创建专属Python环境：

conda create -n gtwr_env python=3.8
conda activate gtwr_env
pip install mgtwr geopandas pandas numpy matplotlib statsmodels

核心工具包功能说明：

mgtwr ：提供GTWR模型实现
geopandas ：处理地理空间数据
statsmodels ：辅助统计检验

1.2 数据格式要求

GTWR模型需要四种核心数据要素：

import numpy as np
import pandas as pd

# 模拟数据生成示例
np.random.seed(2023)
coordinates = np.random.uniform(0, 10, (100, 2))  # 空间坐标(x,y)
timestamps = np.arange(0, 10, 0.1).reshape(-1,1)  # 时间序列
predictors = np.random.normal(0, 1, (100, 3))     # 解释变量
target = np.random.poisson(10, (100, 1))          # 因变量

关键检查点 ：

坐标需转换为平面坐标系（如UTM），直接使用经纬度会导致距离计算偏差
时间变量建议归一化到相同量级（如0-1范围）
缺失值需提前处理，GTWR不支持自动插补

1.3 数据可视化探索

建模前应通过可视化检验时空分布特征：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(coordinates[:,0], coordinates[:,1], timestamps.flatten(), 
           c=target.flatten(), cmap='viridis')
ax.set_xlabel('X坐标')
ax.set_ylabel('Y坐标')
ax.set_zlabel('时间')
plt.colorbar(ax.collections[0], label='目标变量值')
plt.show()

此三维散点图可直观显示目标变量在时空立方体中的分布模式，帮助识别异常值和空间聚集趋势。

2. 模型参数优化实战

2.1 理解关键超参数

GTWR模型有两个核心参数需要优化：

参数	作用	典型范围	影响
带宽(bw)	控制空间邻域范围	0.1-50	值越小局部性越强
时间比例(τ)	平衡时空权重	0.1-10	决定时间衰减速度

专业提示 ：当τ=1时，时空距离被同等对待；τ>1表示时间影响更强；τ<1则空间因素占主导

2.2 自动参数搜索技术

mgtwr包提供了基于黄金分割搜索的自动化参数优化：

from mgtwr.sel_bws import Sel_bws

selector = Sel_bws(
    coords=coordinates,
    t=timestamps,
    y=target,
    X=predictors,
    kernel='gaussian',
    fixed=True
)

best_bw, best_tau = selector.search(
    bw_max=30,
    tau_max=5,
    verbose=True
)

print(f"最优参数：带宽={best_bw:.2f}, τ={best_tau:.2f}")

搜索过程会输出不同参数组合的AICc值（越小越好），最终返回最优参数。对于大型数据集，可以设置 tol=1e-3 提前终止搜索以节省计算时间。

2.3 核函数选择策略

权重函数决定了邻近点的衰减方式，常见选项对比：

高斯核 ：平滑衰减，适合连续变化现象
```
kernel='gaussian'  # 默认选项
```
Bi-square ：明确边界，适合离散分区
```
kernel='bisquare'  # 需指定cutoff距离
```
指数核 ：快速衰减，强调局部特征

可通过交叉验证比较不同核函数的模型表现：

kernels = ['gaussian', 'bisquare', 'exponential']
results = {}

for kern in kernels:
    model = GTWR(coordinates, timestamps, target, predictors,
                bw=best_bw, tau=best_tau, kernel=kern).fit()
    results[kern] = model.R2

pd.DataFrame.from_dict(results, orient='index', columns=['R-squared'])

3. 模型拟合与诊断

3.1 完整建模流程

获得最优参数后，即可拟合最终模型：

from mgtwr.gtwr import GTWR

final_model = GTWR(
    coords=coordinates,
    t=timestamps,
    y=target,
    X=predictors,
    bw=best_bw,
    tau=best_tau,
    kernel='gaussian',
    fixed=True,
    constant=True
).fit()

3.2 关键输出解读

模型对象包含丰富的信息：

print(f"模型R²: {final_model.R2:.3f}")
print(f"调整后R²: {final_model.adj_R2:.3f}")
print(f"AICc: {final_model.AICc:.1f}")

# 回归系数分布统计
betas = pd.DataFrame(final_model.betas)
print(betas.describe().T)

重点关注：

R² >0.7表示模型解释力良好
系数标准差 ：波动过大可能预示过拟合
残差空间自相关 ：Moran's I检验应不显著

3.3 残差时空分析

检验模型残差的时空分布模式：

residuals = final_model.resid_response

# 空间自相关检验
from esda.moran import Moran
moran = Moran(residuals, final_model.W)
print(f"Moran's I: {moran.I:.3f}, p-value: {moran.p_norm:.4f}")

# 时间自相关检验
from statsmodels.tsa.stattools import acf
acf_values = acf(residuals, nlags=5)
print("残差自相关系数:", acf_values)

健康模型的残差应满足：

空间随机性（Moran's I不显著）
无时间自相关（ACF接近0）
正态分布（QQ图检验）

4. 结果可视化与应用

4.1 时空变化模式制图

将回归系数映射到地理空间：

import geopandas as gpd
from shapely.geometry import Point

# 创建地理数据框
geometry = [Point(xy) for xy in coordinates]
gdf = gpd.GeoDataFrame(
    pd.DataFrame(final_model.betas),
    geometry=geometry
)

# 绘制系数空间分布
fig, ax = plt.subplots(1, 3, figsize=(18,6))
for i, col in enumerate(['截距', '变量1', '变量2']):
    gdf.plot(column=i, legend=True, ax=ax[i],
            cmap='coolwarm', markersize=50)
    ax[i].set_title(f'{col}系数空间分布')
plt.tight_layout()

4.2 时间动态可视化

展示系数随时间的变化趋势：

trend_data = pd.DataFrame({
    '时间': timestamps.flatten(),
    '变量1系数': final_model.betas[:,1],
    '变量2系数': final_model.betas[:,2]
})

fig, ax = plt.subplots(figsize=(10,6))
for col in ['变量1系数', '变量2系数']:
    ax.plot(trend_data['时间'], trend_data[col], label=col)
ax.set_xlabel('时间')
ax.set_ylabel('回归系数')
ax.legend()
plt.show()