Python实战GTWR模型从数据预处理到参数优化的全流程解析时空地理加权回归GTWR作为地理加权回归GWR的时空扩展版本能够捕捉变量关系在空间和时间维度上的异质性。本文将带您从零开始实现一个完整的GTWR分析流程重点解决实际应用中常见的数据格式、参数调优和结果解读问题。1. 环境准备与数据生成在开始GTWR建模前需要确保Python环境中已安装必要的库。推荐使用conda创建虚拟环境conda create -n gtwr_env python3.8 conda activate gtwr_env pip install mgtwr geopandas numpy pandas matplotlib模拟数据生成是理解GTWR模型的基础步骤。我们创建一个包含时空属性的合成数据集import numpy as np import pandas as pd np.random.seed(42) n_samples 1000 # 生成空间坐标模拟经纬度 lng np.random.uniform(115, 117, n_samples) lat np.random.uniform(35, 37, n_samples) # 生成时间序列模拟天数为单位 time np.random.randint(1, 365, n_samples) # 生成解释变量 x1 np.random.normal(10, 2, n_samples) x2 np.random.poisson(5, n_samples) # 生成具有时空异质性的系数 beta0 5 0.1*lng 0.05*time/365 beta1 3 0.2*np.sin(lat/10) 0.1*np.cos(time/30) beta2 -2 0.15*(lng-116)**2 0.1*np.log(time1) # 生成因变量 epsilon np.random.normal(0, 1, n_samples) y beta0 beta1*x1 beta2*x2 epsilon # 构建DataFrame data pd.DataFrame({ lng: lng, lat: lat, time: time, x1: x1, x2: x2, y: y })注意实际应用中空间坐标通常需要转换为投影坐标系如UTM的欧式距离而非直接使用经纬度。这是GTWR建模中常见的第一个坑。2. 数据预处理关键步骤GTWR对输入数据格式有特定要求不当的数据准备会导致模型无法运行或结果不可靠。2.1 坐标系统转换经纬度坐标需要转换为平面坐标系如UTM的x,y坐标from pyproj import Proj, transform # WGS84经纬度转UTM 50N适用于中国大部分地区 wgs84 Proj(initepsg:4326) utm50n Proj(initepsg:32650) # 批量转换坐标 x, y transform(wgs84, utm50n, data[lng].values, data[lat].values) data[x_proj] x data[y_proj] y2.2 时空数据标准化由于空间和时间维度的量纲不同需要进行标准化处理from sklearn.preprocessing import MinMaxScaler scaler MinMaxScaler() data[[x_proj,y_proj]] scaler.fit_transform(data[[x_proj,y_proj]]) data[time_norm] scaler.fit_transform(data[[time]])2.3 输入数据结构准备GTWR模型要求的输入格式为NumPy数组# 准备模型输入 coords data[[x_proj, y_proj]].values t data[time_norm].values.reshape(-1,1) X data[[x1, x2]].values y data[y].values.reshape(-1,1)3. 带宽参数优化技术带宽参数(bw)和时间比例参数(tau)的选择直接影响模型性能。mgtwr包提供了Sel_bws类进行参数搜索。3.1 黄金分割搜索法from mgtwr.sel_bws import Sel_bws # 初始化参数搜索器 sel Sel_bws(coords, t, y, X, kernelgaussian, fixedTrue) # 执行参数搜索 bw, tau sel.search(bw_max20, tau_max5, verboseTrue) print(f最优带宽: {bw:.2f}, 最优时间比例: {tau:.2f})参数搜索过程中需要注意bw_max和tau_max的设置应基于数据特性对于大型数据集搜索过程可能耗时较长可以尝试不同的核函数gaussian或bisquare3.2 交叉验证方法除了内置的黄金分割法也可以实现k折交叉验证from sklearn.model_selection import KFold def cv_gtwr(coords, t, y, X, bw_candidates, tau_candidates, k5): kf KFold(n_splitsk) scores np.zeros((len(bw_candidates), len(tau_candidates))) for i, bw in enumerate(bw_candidates): for j, tau in enumerate(tau_candidates): r2_scores [] for train_idx, test_idx in kf.split(X): # 训练模型 model GTWR(coords[train_idx], t[train_idx], y[train_idx], X[train_idx], bw, tau, kernelgaussian).fit() # 测试集预测 pred model.predict(coords[test_idx], t[test_idx], X[test_idx]) # 计算R2 r2 1 - ((y[test_idx] - pred)**2).sum() / ((y[test_idx] - y[test_idx].mean())**2).sum() r2_scores.append(r2) scores[i,j] np.mean(r2_scores) best_idx np.unravel_index(scores.argmax(), scores.shape) return bw_candidates[best_idx[0]], tau_candidates[best_idx[1]]4. 模型拟合与结果分析获得最优参数后可以进行GTWR模型拟合from mgtwr.gtwr import GTWR # 使用最优参数拟合模型 gtwr_model GTWR(coords, t, y, X, bwbw, tautau, kernelgaussian, fixedTrue).fit() # 模型摘要 print(f模型R平方: {gtwr_model.R2:.4f}) print(f调整后R平方: {gtwr_model.adj_R2:.4f}) print(fAIC: {gtwr_model.AIC:.2f})4.1 回归系数分析GTWR模型的系数具有时空异质性可以通过空间可视化来解读import matplotlib.pyplot as plt # 提取第一个解释变量的系数 beta_x1 gtwr_model.betas[:,1] # 创建空间系数分布图 plt.figure(figsize(10,8)) sc plt.scatter(data[lng], data[lat], cbeta_x1, cmapcoolwarm, s20) plt.colorbar(sc, labelx1系数值) plt.xlabel(经度) plt.ylabel(纬度) plt.title(x1变量的空间系数分布) plt.show()4.2 结果可视化技巧时空分析结果可以通过动画展示时间维度变化from matplotlib.animation import FuncAnimation # 按时间切片可视化 fig, ax plt.subplots(figsize(10,8)) time_slices np.linspace(data[time].min(), data[time].max(), 12) def update(frame): ax.clear() mask (data[time] frame-15) (data[time] frame15) sc ax.scatter(data.loc[mask, lng], data.loc[mask, lat], cgtwr_model.betas[mask,1], cmapcoolwarm, s20) ax.set_title(f时间窗口: {frame-15}-{frame15}天) ax.set_xlabel(经度) ax.set_ylabel(纬度) return sc, ani FuncAnimation(fig, update, framestime_slices, interval500) plt.close()提示对于大型数据集建议将动画保存为文件而非实时显示以节省内存。5. 模型诊断与优化GTWR模型需要验证其假设是否成立并进行必要的优化。5.1 残差空间自相关检验使用Morans I检验残差的空间自相关from esda.moran import Moran from libpysal.weights import DistanceBand # 创建空间权重矩阵 w DistanceBand(coords, threshold0.1, binaryFalse) # 计算Morans I moran Moran(gtwr_model.resid, w) print(fMorans I: {moran.I:.3f}, p-value: {moran.p_norm:.4f})5.2 多重共线性诊断检查解释变量间的共线性问题from statsmodels.stats.outliers_influence import variance_inflation_factor # 计算VIF vif_data pd.DataFrame() vif_data[feature] [x1, x2] vif_data[VIF] [variance_inflation_factor(X, i) for i in range(X.shape[1])] print(vif_data)5.3 模型比较与普通线性回归(OLS)和GWR模型进行比较import statsmodels.api as sm from mgtwr.gtwr import GWR # OLS模型 X_ols sm.add_constant(X) ols sm.OLS(y, X_ols).fit() # GWR模型 gwr_model GWR(coords, y, X, bw0.1, kernelgaussian).fit() # 比较指标 results pd.DataFrame({ Model: [OLS, GWR, GTWR], R2: [ols.rsquared, gwr_model.R2, gtwr_model.R2], AIC: [ols.aic, gwr_model.AIC, gtwr_model.AIC] }) print(results)6. 实际应用案例以一个模拟的房价数据集展示GTWR的实际应用价值。假设我们收集了某城市2年间的房价数据包含因变量每平方米房价解释变量到市中心距离、周边学校数量、房龄时空坐标小区经纬度和交易月份# 案例数据准备 np.random.seed(123) n 1500 lng np.random.uniform(116.2, 116.6, n) lat np.random.uniform(39.8, 40.0, n) month np.random.randint(1, 25, n) # 转换为投影坐标 x, y transform(wgs84, utm50n, lng, lat) # 生成解释变量 dist_center np.sqrt((x - 434000)**2 (y - 4430000)**2)/1000 schools np.random.poisson(3 2*(lat - 39.9), n) age np.random.randint(1, 30, n) # 生成具有时空异质性的房价 price_base 50000 10000*(lat - 39.9) price_trend 1000*(month/12) price_dist -8000*dist_center*(1 0.2*np.sin(month/3)) price_school 3000*schools*(1 0.1*np.cos(dist_center)) price_age -500*age*(1 0.05*(month-12)) price price_base price_trend price_dist price_school price_age np.random.normal(0, 3000, n) # 准备GTWR输入 coords np.column_stack([x, y]) t month.reshape(-1,1) X np.column_stack([dist_center, schools, age]) y price.reshape(-1,1)通过GTWR分析我们可以发现房价影响因素的空间分布和时间变化规律为房地产投资决策提供更精细化的参考。