GTWR模型实战：从理论到Python代码实现

1. GTWR模型基础理解时空加权回归GTWRGeographically and Temporally Weighted Regression是地理加权回归GWR在时空维度的扩展。我第一次接触这个模型是在分析城市房价时空变化时发现传统回归方法完全无法捕捉区域差异。GTWR的核心思想很简单不同地点、不同时间影响因素的作用强度会变化。比如地铁开通对房价的影响在市中心和郊区、在开通前后都会呈现不同规律。模型公式看起来复杂但其实拆解后很好理解y_i β_0(u_i,v_i,t_i) Σ[β_k(u_i,v_i,t_i) * x_ik] ε_i其中(u_i,v_i)是空间坐标t_i是时间点。与普通线性回归最大的区别在于每个系数β都带着时空坐标下标这意味着在北京市朝阳区2020年和海淀区2022年同一个影响因素比如学区的系数可以完全不同。时空权重的计算是GTWR的精华所在。我常用的是高斯核函数w_ij exp(-(d_ij^2/h^2))这里的d_ij包含空间距离和时间距离h是带宽参数。实际项目中我发现时空比例参数τ的选择会显著影响结果。有次分析空气质量数据时τ值过小导致时间维度被过度压缩最终模型把雾霾归因到了完全错误的因素上。2. 环境准备与数据预处理2.1 安装必要的Python库推荐使用conda创建虚拟环境避免包冲突conda create -n gtwr_env python3.8 conda activate gtwr_env pip install mgtwr geopandas pandas numpy matplotlib我测试过多个版本组合发现mgtwr 1.0.4与Python 3.8的兼容性最好。曾经在3.10环境下遇到C扩展编译失败的问题折腾半天才发现是Python版本太新。2.2 数据格式要求GTWR需要四种核心数据空间坐标经度纬度建议先投影为平面坐标时间戳统一为数值型如2020.5表示2020年6月自变量数值矩阵因变量一维数组常见坑点直接用经纬度会导致距离计算错误。我在某次项目中使用WGS84坐标直接计算结果模型完全失效。正确做法是先转为UTM等投影坐标import geopandas as gpd from pyproj import Transformer # WGS84转UTM transformer Transformer.from_crs(EPSG:4326, EPSG:32650) points gpd.points_from_xy(df[lng], df[lat]) coords np.array([transformer.transform(lng, lat) for lng, lat in zip(df[lng], df[lat])])时间数据也要标准化。处理COVID-19数据时我最初用日期字符串导致计算异常后来转换为相对天数df[days] (pd.to_datetime(df[date]) - pd.Timestamp(2020-01-01)).dt.days3. 参数优化实战技巧3.1 带宽选择黄金分割搜索法带宽bw控制着时空影响的衰减速度。mgtwr库提供了自动搜索方法from mgtwr.sel_bws import Sel_bws sel Sel_bws(coords, t, y, X, kernelgaussian, fixedTrue) bw, tau sel.search(bw_max40, tau_max5, verboseTrue)几个经验值城市尺度数据bw通常在0.5-5公里省级范围可能需10-50公里时间带宽按数据周期比例调整我曾用交叉验证对比过不同方法发现AICc准则在中小样本量时最稳定。当数据点超过1万时改用BIC准则计算效率更高。3.2 时空比例τ的调优τ值决定时空影响的相对权重。一个实用技巧是先做空间GWR和时间TWR比较效果# 纯空间模型 gwr_model GWR(coords, y, X, bw5).fit() # 纯时间模型 twr_model TWR(t, y, X, bw2).fit() # 比较R平方 print(fGWR R2: {gwr_model.R2:.3f}, TWR R2: {twr_model.R2:.3f})如果两者R平方差距大τ应偏向表现好的维度。我处理交通流量数据时发现空间维度解释力更强最终设置τ0.3效果最佳。4. 完整建模流程示例4.1 生成模拟数据创建包含时空趋势的测试数据np.random.seed(42) n 1000 coords np.random.uniform(0, 10, (n, 2)) # 10km范围 t np.random.uniform(0, 365, (n, 1)) # 1年时间 # 创建时空相关的系数 beta1 3 (coords[:,0] t[:,0])/20 beta2 2 * np.sin(coords[:,1]/5) * np.cos(t[:,0]/30) X np.column_stack([np.random.normal(5, 2, n), np.random.poisson(3, n)]) y 5 beta1*X[:,0] beta2*X[:,1] np.random.normal(0, 1, n)4.2 模型训练与评估完整训练流程# 参数搜索 sel Sel_bws(coords, t, y, X, kernelbisquare) bw, tau sel.search(bw_max20, tau_max3) # 模型拟合 model GTWR(coords, t, y, X, bwbw, tautau, kernelbisquare).fit() # 结果分析 print(fR-squared: {model.R2:.3f}) print(f参数分布:\n 截距: {model.betas[:,0].mean():.2f}±{model.betas[:,0].std():.2f}) print(f 变量1: {model.betas[:,1].mean():.2f}±{model.betas[:,1].std():.2f})可视化系数空间分布plt.scatter(coords[:,0], coords[:,1], cmodel.betas[:,1], cmapcoolwarm) plt.colorbar(label变量1系数) plt.title(系数空间分布)4.3 结果解读要点系数稳定性检查健康的模型各点系数应呈现连续变化。某次分析中我发现系数剧烈震荡检查发现是带宽过小导致过拟合。时空模式识别用3D散点图观察系数时空趋势from mpl_toolkits.mplot3d import Axes3D fig plt.figure() ax fig.add_subplot(111, projection3d) ax.scatter(coords[:,0], coords[:,1], t[:,0], cmodel.betas[:,1])残差空间自相关检验用Morans I指数检查是否遗漏空间特征from esda.moran import Moran moran Moran(model.resid, gwr_model.spatial_weights) print(fMorans I: {moran.I:.3f}, p-value: {moran.p_sim:.3f})5. 常见问题解决方案5.1 计算效率优化当数据量超过5000点时使用fixedFalse自适应带宽换用kernelbisquare减少远端点影响分块计算后合并结果我曾用Dask并行处理10万房地产数据import dask.array as da # 将数据转为Dask数组 coords_da da.from_array(coords, chunks(5000, 2)) t_da da.from_array(t, chunks5000) X_da da.from_array(X, chunks(5000, X.shape[1]))5.2 缺失值处理GTWR对缺失值敏感。推荐两种方法时空插值from sklearn.impute import KNNImputer imputer KNNImputer(n_neighbors5) X_filled imputer.fit_transform(X)构建时空缓冲区均值def spatial_buffer_avg(df, radius1000): gdf gpd.GeoDataFrame(df, geometrygpd.points_from_xy(df.lng, df.lat)) gdf[filled] gdf.geometry.apply( lambda p: df[df.geometry.distance(p) radius][value].mean()) return gdf5.3 多重共线性诊断时空变系数会放大共线性问题。解决步骤计算局部VIFfrom statsmodels.stats.outliers_influence import variance_inflation_factor vif [variance_inflation_factor(X[i:i100], j) for i in range(0, len(X), 100) for j in range(X.shape[1])]若VIF10考虑移除高相关变量使用主成分分析降维加入岭回归惩罚项6. 进阶应用案例6.1 房价时空预测以北京二手房为例的关键步骤特征工程# 计算到地铁站距离 df[subway_dist] df.apply(lambda x: get_distance(x[lng], x[lat], subway_lng, subway_lat), axis1) # 时间特征 df[month] (df[year] - 2015)*12 df[month]交互项处理# 学区与时间的交互 X[school*time] X[school] * (t - t.min())结果可视化# 绘制系数时间趋势 plt.plot(pd.to_datetime(df[date]), model.betas[:,1]) plt.fill_between(df[date], model.betas[:,1] - 1.96*model.std_err[:,1], model.betas[:,1] 1.96*model.std_err[:,1], alpha0.2)6.2 空气质量分析处理PM2.5数据时的特殊处理风速方向特征# 转换为矢量分量 X[wind_x] X[wind_speed] * np.cos(X[wind_dir]*np.pi/180) X[wind_y] X[wind_speed] * np.sin(X[wind_dir]*np.pi/180)空间自相关项# 添加邻近点均值 gdf gpd.GeoDataFrame(X, geometrygpd.points_from_xy(coords[:,0], coords[:,1])) X[neighbor_avg] gdf.geometry.apply( lambda p: gdf[gdf.geometry.distance(p) 5000][value].mean())时空滞后项from pysal.model import spreg # 创建时空权重矩阵 w spreg.util.lat2W(coords[:,0], coords[:,1]) X[lagged] spreg.lag_spatial(w, y)