简介本资源是一份面向地理信息科学、环境科学及空间统计方向本科生与初阶研究者的课程论文实践资料聚焦于利用地理加权回归GWR模型与克里金插值法协同预测空气质量指数AQI解决传统全局模型难以刻画空间异质性的问题。压缩包为1个1.08MB的docx文档完整包含引言、研究方法含GWR数学推导与带宽优化策略、克里金变异函数原理、GWR与克里金双案例分析含系数显著性解读、预测结果对比、总结讨论以及附录中可直接运行的Python源代码GWR建模、克里金插值、地图可视化三部分。内容结构严谨公式推导清晰代码注释详实特别适合学习空间回归建模逻辑、理解局部参数估计思想及掌握GIS空间预测全流程。目前已有1258人学习下载是理论结合编程、论文写作与实操复现一体化的优质教学参考材料。1. 这不是又一个“画热力图”的Demo它用GWR克里金双模型把北京16个监测点的AQI预测误差压到±20.6以内你手头有一堆空气质量监测站数据——位置、时间、AQI值还有气温、气压、湿度、风速这些气象变量。想做个预测直接扔进LSTM或XGBoost先别急。这篇课程论文干了一件更“地理”的事它没把空间当背景板而是让每个监测点都拥有自己的一套回归参数GWR再用克里金把残差里的空间结构挖出来补上泛克里金插值。最终对第23期的预测GWR均方误差449.4克里金423.3——换算成AQI单位就是±20.6以内。这不是理论推演是实打实跑在真实北京16个站点朝阳、海淀、密云、延庆…上的双模型闭环。它不依赖卫星遥感或高精度气象再分析数据只靠常规地面观测站坐标基础气象要素就能给出带空间解释性的预测结果。适合正在做环境GIS、城市健康评估、或者需要向环保部门交付“为什么东边预测准、西边不准”这类可解释结论的开发者。如果你的项目卡在“模型黑匣子输出一个数但领导问‘昌平区为什么比大兴区预测偏差大’时答不上来”那这个源码包就是你的后悔药。2. 项目结构与技术栈Python生态下的轻量级空间建模组合拳这个项目不是用ArcGIS Pro点几下出来的也不是调用某个云API封装的黑盒。它是一套可调试、可替换、可嵌入你现有数据流水线的Python工程。核心逻辑全部落在三个附录代码文件里没有隐藏依赖没有私有包全是PyPI上能pip install的标准库和地理计算主力——libpysal、mgwr、scikit-learn、numpy、pandas、geopandas、rasterio。我拆开整个资源包后确认它只有3个核心.py文件GWR建模、克里金插值、地图可视化1个CSV格式的原始数据集含time、lat、lon、aqi及7个气象变量外加1份说明性PDF。没有数据库、没有Web服务、没有Dockerfile——它就是一个“拿来即跑”的分析脚本集合但恰恰因为足够轻反而更容易被你塞进自己的Flask API、Airflow DAG或者Jupyter Notebook工作流里。2.1 目录结构与文件职责三文件驱动完整流程整个项目目录极简但分工明确/air_quality_gwr_kriging/ ├── data/ │ └── beijing_aqi_weather.csv # 核心数据16个站点 × 多期观测含lat/lon/aqi/temp/hpa/wet/speed/dir/height/mm ├── src/ │ ├── gwr_model.py # 主体GWR建模、带宽优化、系数估计、t检验、预测 │ ├── kriging_interpolation.py # 主体变异函数拟合高斯模型、泛克里金权重求解、残差插值 │ └── viz_map.py # 辅助生成回归系数空间分布图、预测值热力图、显著性mask图 ├── requirements.txt # 明确列出libpysal4.12.1, mgwr2.2.3, scikit-learn1.3.2, ... └── README.md # 一句话说明如何准备数据、如何顺序执行三个脚本提示beijing_aqi_weather.csv是项目的生命线。它不是模拟数据而是真实采样——每行是一个站点在某期time的观测lat/lon是WGS84坐标aqi是因变量其余7列是自变量mm全为0故被剔除。你替换成自己城市的监测站CSV只要字段名一致几乎不用改代码。2.2 技术选型深意为什么是mgwr而不是statsmodels为什么泛克里金不用gstools很多人看到“地理加权回归”第一反应是statsmodels但它只支持全局OLS。GWR必须用专为空间建模设计的库。这里选mgwrMulti-scale Geographically Weighted Regression而非更老的gwr包是因为它原生支持AICc带宽优化、ENP有效参数个数计算、并内置了bisquare核函数——这正是论文里明确使用的。而克里金部分没选gstools功能强但偏重地质建模也没用pykrige接口略重而是手写最小二乘求解泛克里金方程组见kriging_interpolation.py中solve_kriging_system()函数。原因很实在论文明确要求“泛克里金”即带漂移项trend的克里金而pykrige的UniversalKriging虽支持但其内部用scipy.optimize求解对16个点的小样本容易数值不稳定手写矩阵求解np.linalg.solve反而更可控、更透明方便你debug权重是否发散、变异函数参数是否合理。2.3 数据预处理逻辑坐标标准化与变量清洗是预测稳定的地基GWR对坐标尺度极度敏感。如果直接把lat/lon度和temp℃、hpa百帕扔进同一个模型距离权重会完全被数值大的变量主导。论文没明说但代码里藏着关键一步# src/gwr_model.py 片段 def prepare_coordinates(coords): 将经纬度转换为平面米制坐标近似 使用北京中心点(39.9, 116.3)作为原点简化投影 避免调用pyproj降低依赖 lat_center, lon_center 39.9, 116.3 # 近似1度纬度 ≈ 111km1度经度 ≈ 111km * cos(lat) km_per_lat 111.0 km_per_lon 111.0 * np.cos(np.radians(lat_center)) x (coords[:, 1] - lon_center) * km_per_lon * 1000 # 转为米 y (coords[:, 0] - lat_center) * km_per_lat * 1000 # 转为米 return np.column_stack([x, y])这段代码把lat/lon转成以北京为中心的局部平面坐标单位米让bisquare核函数计算的空间距离d有物理意义。如果你的数据不在北京只需改lat_center/lon_center两行。同样mm降水量列全为0代码里直接df df.drop(mm, axis1)——这种“硬规则”清洗比用df.dropna()更安全避免误删有效行。2.4 模型输入输出契约定义清楚才能无缝集成这三个脚本不是孤立的。它们通过约定好的DataFrame结构传递数据形成清晰的输入输出契约脚本输入输出关键字段gwr_model.pypd.DataFrame含lat,lon,aqi,temp,hpa, ...dict含coefficients(16×8),t_values(16×8),predictions(16,)coefficients是核心每行对应一个站点的8个回归系数截距7变量kriging_interpolation.pypd.DataFrame含lat,lon,residuals来自GWR预测残差np.array长度为16的插值后残差修正值residuals必须与GWR输出的predictions严格对齐索引顺序不能乱viz_map.py合并后的gwr_resultkriging_correctionPNG图片系数空间图、最终AQI预测图所有图都带crsepsg:4326确保地理坐标正确这个契约意味着你可以把gwr_model.py的输出存成gwr_output.pkl再喂给kriging_interpolation.py最后用viz_map.py画图——完全解耦。如果你想把GWR换成你自己的神经网络模型只要保证输出predictions和residuals字段后续流程一动不动。3. GWR模型核心实现从带宽搜索到系数t检验的逐行拆解GWR不是“调个包就完事”。它的灵魂在于带宽bandwidth的选择——太小过拟合太大退化成OLS。论文用CV交叉验证最小化代码里实现了完整的Golden Section Search黄金分割法。这不是教科书伪代码是真正在16个点上跑出最优bs12450.3米的生产级实现。下面带你逐段看透gwr_model.py最硬核的200行。3.1 带宽优化黄金分割法不是玄学是收敛保证mgwr库本身提供SearchGWRParameter类但论文代码选择手写原因很实际便于打印每一步的CV值便于你观察搜索轨迹也便于在带宽不收敛时加print调试。核心函数find_optimal_bandwidth()如下# src/gwr_model.py def find_optimal_bandwidth(X, y, coords, kernelbisquare, max_iter30): 黄金分割法搜索最优带宽最小化CV值 X: (n, k) 自变量矩阵已标准化 y: (n,) 因变量向量 coords: (n, 2) 平面坐标米 # 定义搜索区间基于点间最大距离的10%~90% dist_matrix squareform(pdist(coords)) max_dist np.max(dist_matrix) a, b 0.1 * max_dist, 0.9 * max_dist # 黄金分割常数 r (np.sqrt(5) - 1) / 2 x1 a (1 - r) * (b - a) x2 a r * (b - a) f1 _cv_score(x1, X, y, coords, kernel) f2 _cv_score(x2, X, y, coords, kernel) for i in range(max_iter): if f1 f2: b x2 x2 x1 f2 f1 x1 a (1 - r) * (b - a) f1 _cv_score(x1, X, y, coords, kernel) else: a x1 x1 x2 f1 f2 x2 a r * (b - a) f2 _cv_score(x2, X, y, coords, kernel) # 返回CV更小的带宽 optimal_bs x1 if f1 f2 else x2 return optimal_bs def _cv_score(bandwidth, X, y, coords, kernel): 计算给定带宽下的CV值 n len(y) cv_sum 0 for i in range(n): # 留一法排除第i个点 mask np.ones(n, dtypebool) mask[i] False X_train, y_train X[mask], y[mask] coords_train coords[mask] coords_test coords[[i]] # 计算第i点的权重仅对训练点 weights _compute_weights(coords_test, coords_train, bandwidth, kernel) # 加权最小二乘 W np.diag(weights) try: beta np.linalg.solve(X_train.T W X_train, X_train.T W y_train) y_pred X[[i]] beta except np.linalg.LinAlgError: y_pred np.mean(y_train) # 退化为均值 cv_sum (y[i] - y_pred) ** 2 return np.sqrt(cv_sum / n)逻辑说明与参数说明a, b不是随便设的而是基于数据点间最大欧氏距离动态确定避免在稀疏数据上搜过宽、在密集数据上搜过窄。_cv_score里np.linalg.solve可能失败矩阵奇异代码做了try/except兜底为均值预测——这是真实项目必须的鲁棒性不是教科书理想情况。weights计算调用_compute_weights它实现bisquare核w (1 - (d/b)^2)^2 if d b else 0其中d是测试点到训练点的距离。bandwidth b就是这里的bs单位是米所以前面坐标必须转成米3.2 系数估计与t检验ENP有效参数个数是GWR可信度的命门GWR的系数不是独立的它们通过空间权重耦合。判断一个系数是否真的显著不能用OLS的t检验公式。必须用帽子矩阵S的迹trace来校正自由度。论文公式8ENP tr(S)在这里落地# src/gwr_model.py def fit_gwr_model(X, y, coords, bandwidth, kernelbisquare): 执行GWR拟合返回所有站点的系数、t值、标准误 n, k X.shape coefficients np.zeros((n, k)) t_values np.zeros((n, k)) std_errors np.zeros((n, k)) # 预分配帽子矩阵Sn x n用于计算ENP S np.zeros((n, n)) for i in range(n): weights _compute_weights(coords[[i]], coords, bandwidth, kernel) W np.diag(weights) # 局部加权最小二乘 try: XtWX X.T W X XtWy X.T W y beta_i np.linalg.solve(XtWX, XtWy) coefficients[i] beta_i # 计算帽子矩阵第i行Si xi (XtWX)^-1 XtW inv_XtWX np.linalg.inv(XtWX) Si X[i:i1] inv_XtWX X.T W S[i] Si # 计算该点残差和标准误 y_pred_i X[i] beta_i residuals y - (X beta_i) # 注意这是用beta_i预测所有点非严格留一 sigma2 np.sum(residuals**2) / (n - np.trace(S)) # 用当前S的迹校正 var_beta_i np.diag(inv_XtWX) * sigma2 std_errors[i] np.sqrt(var_beta_i) t_values[i] beta_i / (std_errors[i] 1e-8) # 防0除 except np.linalg.LinAlgError: coefficients[i] np.nan t_values[i] np.nan std_errors[i] np.nan # 最终ENP trace(S)用于模型比较AICc enp np.trace(S) return { coefficients: coefficients, t_values: t_values, std_errors: std_errors, enp: enp, sigma2: sigma2 }关键点解析S[i] Si这一行把每个点的帽子矩阵行存下来最后np.trace(S)就是ENP。这是GWR区别于OLS的核心——OLS的ENP恒等于k变量数而GWR的ENP在k到n之间反映模型复杂度。residuals计算用了X beta_i即用第i点的系数预测所有点。这虽不是严格留一但对16个点的小样本计算量可接受且论文结果表明其t检验结论稳定。1e-8防0除是血泪经验当某个变量在局部区域完全不变化如某站点speed全为0var_beta_i可能为0t值爆炸加微小扰动保程序不死。3.3 预测与评估为什么论文用第23期数据因为它在时间序列上是“下一个”GWR本身是横截面模型single time slice但论文用time22拟合time23预测本质是假设空间关系在短期内稳定。代码里predict_aqi_for_time()函数明确体现这一逻辑# src/gwr_model.py def predict_aqi_for_time(df, time_train22, time_pred23): 用time_train期数据拟合GWR预测time_pred期的AQI train_mask df[time] time_train pred_mask df[time] time_pred df_train df[train_mask].copy() df_pred df[pred_mask].copy() # 准备训练数据 coords_train prepare_coordinates(df_train[[lat, lon]].values) X_train df_train[[temp, hpa, wet, speed, dir, height]].values y_train df_train[aqi].values # 添加截距项 X_train np.column_stack([np.ones(len(X_train)), X_train]) # 搜索最优带宽 optimal_bs find_optimal_bandwidth(X_train, y_train, coords_train) # 拟合模型 gwr_result fit_gwr_model(X_train, y_train, coords_train, optimal_bs) # 预测用训练期拟合的系数乘以预测期的自变量 X_pred df_pred[[temp, hpa, wet, speed, dir, height]].values X_pred np.column_stack([np.ones(len(X_pred)), X_pred]) predictions X_pred gwr_result[coefficients].T # 注意coefficients是16x8X_pred是16x8 # 计算MSE y_true df_pred[aqi].values mse np.mean((y_true - predictions)**2) return predictions, mse, gwr_result注意predictions X_pred gwr_result[coefficients].T这一行是关键。gwr_result[coefficients]是16×8矩阵每一行是对应站点的系数X_pred也是16×8每一行是该站点在time_pred期的自变量。矩阵乘法后得到16×16预测矩阵再取对角线np.diag()才是每个站点用“自己的系数”预测“自己的AQI”。论文表4的16个预测值就是这么来的。4. 克里金插值核心实现从变异函数拟合到泛克里金方程组求解GWR解决了“主趋势”的空间异质性但残差里还藏着未被捕捉的空间自相关——比如密云和怀柔的残差总是同号而通州和石景山总是反号。克里金就是来捕获这个的。论文用泛克里金Universal Kriging因为它允许残差存在一个全局漂移trend这比普通克里金Ordinary Kriging更符合空气质量残差的实际它往往随地形高度有缓慢变化。kriging_interpolation.py没有调用高级库而是用scipy.optimize.curve_fit拟合变异函数再手解线性方程组全程可控。4.1 变异函数拟合高斯模型参数怎么来的Q1/Q2/cR指标怎么看变异函数Semivariogram是克里金的地基。它描述“两点间差异的平方期望值”随距离的变化。论文选用高斯模型Gaussian model因其在短距离平滑在长距离渐近于基台值符合大气扩散的物理直觉。拟合代码如下# src/kriging_interpolation.py def fit_variogram_model(distances, semivariances, modelgaussian): 用curve_fit拟合变异函数模型 distances: 1D array, 所有点对距离 semivariances: 1D array, 对应点对的半方差 def gaussian_model(h, nugget, sill, range_param): 高斯模型γ(h) nugget sill * (1 - exp(-(h/range)^2)) return nugget sill * (1 - np.exp(-(h / range_param) ** 2)) # 初始参数块金semivariances[0], 基台semivariances.max(), 变程distances.max()/3 p0 [np.min(semivariances), np.max(semivariances), np.max(distances)/3] try: params, pcov curve_fit(gaussian_model, distances, semivariances, p0p0, maxfev5000) nugget, sill, range_param params except RuntimeError: # 拟合失败返回默认参数 nugget, sill, range_param 38.8, 0.0295, 7.59 # 论文表5的值 # 计算拟合优度指标Q1, Q2, cR y_pred gaussian_model(distances, nugget, sill, range_param) Q1 np.sum((semivariances - y_pred) ** 2) / np.sum((semivariances - np.mean(semivariances)) ** 2) Q2 1 - np.sum((semivariances - y_pred) ** 2) / np.sum((y_pred - np.mean(semivariances)) ** 2) cR np.sum(np.abs(semivariances - y_pred)) / np.sum(semivariances) return { model: gaussian, nugget: nugget, sill: sill, range: range_param, Q1: Q1, Q2: Q2, cR: cR, fitted_curve: (distances, y_pred) } # 调用示例 coords prepare_coordinates(df[[lat, lon]].values) # 同GWR用米制坐标 dist_matrix squareform(pdist(coords)) # 计算所有点对的半方差γ(h) 0.5*(z_i - z_j)^2 semivariances [] for i in range(len(coords)): for j in range(i1, len(coords)): h dist_matrix[i, j] gamma_h 0.5 * (residuals[i] - residuals[j]) ** 2 semivariances.append((h, gamma_h)) distances, semis zip(*semivariances) variogram_result fit_variogram_model(np.array(distances), np.array(semis))参数说明与避坑nugget块金值代表微观尺度变异或测量误差论文值38.8说明残差噪声不小。sill基台值是变异函数的上限论文0.0295很小意味着残差整体空间相关性弱——这解释了为什么克里金提升有限MSE从449.4降到423.3。range变程7.59单位是“千公里”不是论文里没写单位但结合distances是米制range_param应为7590米7.59km即在此距离内残差仍有相关性。Q1越小越好论文0.63Q2越接近1越好论文1.12略超1说明拟合稍过但可接受cR越小越好论文66.68较大提示拟合一般——这就是论文说“模型拟合效果一般”的依据。4.2 泛克里金方程组为什么必须解 (nq)×(nq) 的大矩阵泛克里金比普通克里金多q个漂移项trend terms。论文中q2因为漂移设为线性μ(s) β0 β1*x β2*yx,y是平面坐标。因此待求未知数是n个权重λi q个漂移系数βj共nq个。方程组如下[ γ(s_i, s_j) x_j y_j ] [ λ1 ] [ γ(s_0, s_1) ] [ x_i 0 0 ] [ λ2 ] [ x_0 ] [ y_i 0 0 ] * [ ...] [ y_0 ]其中γ(s_i, s_j)是点i,j间的变异函数值γ(s_0, s_i)是待估点s0到点i的变异函数值。代码实现# src/kriging_interpolation.py def universal_kriging(coords, residuals, variogram_params, coords_pred): 泛克里金插值 coords: (n, 2) 训练点坐标米 residuals: (n,) 训练点残差 coords_pred: (m, 2) 待预测点坐标米 n len(coords) m len(coords_pred) q 2 # 线性漂移x, y # 构建变异函数矩阵Gamma (n x n) Gamma np.zeros((n, n)) for i in range(n): for j in range(n): h np.linalg.norm(coords[i] - coords[j]) Gamma[i, j] variogram_params[nugget] \ variogram_params[sill] * (1 - np.exp(-(h / variogram_params[range])**2)) # 构建漂移矩阵F (n x q): [x, y] F coords.copy() # (n, 2) # 构建右端项向量gamma_0 (nq, 1) gamma_0 np.zeros(n q) for i in range(n): h np.linalg.norm(coords_pred[0] - coords[i]) # 只预测第一个点循环处理 gamma_0[i] variogram_params[nugget] \ variogram_params[sill] * (1 - np.exp(-(h / variogram_params[range])**2)) # 构建左端矩阵A (nq) x (nq) A np.zeros((n q, n q)) A[:n, :n] Gamma A[:n, n:nq] F A[n:nq, :n] F.T # 解方程组 A * [lambda; beta] gamma_0 try: solution np.linalg.solve(A, gamma_0) lambdas solution[:n] # 预测值 sum(lambda_i * residual_i) pred_residual np.sum(lambdas * residuals) except np.linalg.LinAlgError: # 退化为普通克里金忽略漂移 pred_residual np.mean(residuals) # 或用IDW return pred_residual # 主函数对所有预测点循环 def interpolate_residuals(coords_train, residuals, coords_pred, variogram_params): corrections [] for i in range(len(coords_pred)): corr universal_kriging(coords_train, residuals, variogram_params, coords_pred[i:i1]) corrections.append(corr) return np.array(corrections)为什么这样设计universal_kriging函数每次只预测一个点是为了避免构建巨大的(m*n) x (m*n)矩阵。对16个点A是18×18np.linalg.solve毫秒级。如果一次预测16个点A会变成(16*162) x (16*162)内存和计算量爆炸。这是小样本克里金的务实选择。5. 避坑指南我在复现时踩过的5个真实坑省你三天调试时间这个项目看着只有3个脚本但空间建模的坑又深又隐蔽。我按论文步骤在本地复现时前两次全军覆没。以下是5条血泪经验每一条都对应一个真实报错和解决方案5.1 坑1坐标系不统一导致GWR带宽搜索发散CV值一路飙升到1e6现象find_optimal_bandwidth()运行中f1和f2在迭代中越来越大最终返回bandwidth1e9后续_compute_weights全为0系数全nan。原因beijing_aqi_weather.csv里的lat/lon是WGS84EPSG:4326但pdist(coords)计算欧氏距离时把度当成了米。1度≈111km所以max_dist算出来是1e6量级a0.1*max_dist就到了1e5米远超北京实际范围约100km导致权重永远不衰减。解决必须用prepare_coordinates()把lat/lon转成局部平面坐标米。检查coords数组的值正常应在[-50000, 50000]范围内而不是[116, 117]或[39, 40]。加一句print(Coord range:, coords.min(), coords.max())即可定位。5.2 坑2mgwr版本冲突导致SearchGWRParameter找不到但论文代码没用它现象pip install mgwr后运行gwr_model.py报错ModuleNotFoundError: No module named mgwr.search。原因论文写于2020年用的是mgwr2.0而当前PyPI最新版是2.2.3API已重构。但论文附录代码是纯手写GWR根本没导入mgwr。你看到requirements.txt里有mgwr其实是作者混淆了——删掉它pip install libpysal4.12.1足矣libpysal提供pdist等空间工具。解决删除requirements.txt中的mgwr行pip uninstall mgwr。gwr_model.py里所有import都只用numpy、scipy、libpysal。5.3 坑3克里金变异函数拟合失败curve_fit抛RuntimeError现象fit_variogram_model()中curve_fit报Optimal parameters not found程序走except RuntimeError分支返回默认参数但后续universal_kriging中np.linalg.solve又报奇异矩阵。原因residuals向量里有inf或nan可能来自GWR的std_errors计算失败导致semivariances计算出错。解决在interpolate_residuals()开头加清洗residuals np.nan_to_num(residuals, nannp.nanmean(residuals), posinfnp.nanmean(residuals), neginfnp.nanmean(residuals))同时print(Residuals stats:, np.min(residuals), np.max(residuals), np.isnan(residuals).sum())监控。5.4 坑4预测时维度错位X_pred coefficients.T得到16×16矩阵而非16×1现象predictions是形状(16, 16)的二维数组取np.diag(predictions)才得到16个值但论文表4是列向量。原因coefficients是(16, 8)X_pred是(16, 8)运算符在numpy中是矩阵乘X_pred coefficients.T结果是(16, 16)。正确做法是逐点相乘predictions np.array([X_pred[i] coefficients[i] for i in range(16)])。解决修改predict_aqi_for_time()中预测部分predictions np.zeros(len(X_pred)) for i in range(len(X_pred)): predictions[i] X_pred[i] gwr_result[coefficients][i]5.5 坑5地图可视化报错ValueError: Invalid crsgeopandas找不到CRS现象viz_map.py运行到gdf gpd.GeoDataFrame(..., crsepsg:4326)时报错。原因geopandas新版要求CRS必须是pyproj.CRS对象字符串epsg:4326不再被接受。解决升级pyproj并改写import pyproj gdf gpd.GeoDataFrame(df, geometrygpd.points_from_xy(df.lon, df.lat), crspyproj.CRS(epsg:43本文还有配套的精品资源点击获取