简介：这份资源是一篇完整的课程论文，围绕地理加权回归（GWR）模型与克里金插值法在空气质量指数预测中的应用展开，面向地理信息科学、空间统计与遥感方向的学生及研究者，适合具备一定Python与GIS基础、希望掌握空间非平稳性建模与地统计插值的读者。压缩包内仅含1个docx文档，约1.08MB，正文涵盖引言、研究方法、GWR案例分析、回归系数显著性检验、预测评估、克里金案例及总结讨论，并附有GWR预测、克里格插值、地图可视化三部分Python源代码，便于读者复现完整流程。论文详细推导了GWR模型形式、核函数与带宽优化（CV、AICc）过程，对比了GWR与OLS在捕捉空间局部效应上的差异，并借助克里金插值生成连续污染分布图，涉及RMSE、R²等精度评价指标。目前已有1256人学习下载，可作为空间分析课程设计、论文写作与实操练习的参考范例。
1. 从一张空气质量分布图说起：这套 GWR + 克里金方案到底解决什么问题
空气质量指数预测这件事，最容易被低估的坑是：你拿普通线性回归跑一遍，R² 只有 0.4 出头，模型告诉你「气温每升高 1 度，AQI 平均降多少」，可这个「平均」在城区和山区根本不是一回事。空间异质性摆在那里，全局模型给出的系数是假的——它把不同地点的规律强行揉成了一个数。这套方案要解决的就是这个问题：用地理加权回归（GWR）把回归参数从全局拆到局部，让每个监测点都有自己的系数；再用克里金插值把离散站点的预测值铺成连续曲面，补上没建站区域的空白。适合正在做空间分析方向课程设计、毕设的同学，也适合需要出一份「能跑通、能解释、能画图」的空气质量分析方案的从业者。整份方案包含 GWR 建模、克里金插值、显著性检验和地图可视化四块，附录里给了完整的 Python 源代码，照着改数据和坐标就能复现。
2. 总体技术路线：GWR 与克里金的分工和衔接
2.1 为什么不是「二选一」而是「两段式」
很多人第一次接触这两个方法时会纠结：到底用 GWR 还是克里金？其实它们解决的是两个不同阶段的问题。GWR 负责建模——给定一组自变量（气温、气压、湿度、风速、风向、高度等）和因变量（AQI），在每一个空间位置上估计一组局部回归系数，回答「哪些因素在哪些区域对 AQI 影响更大」。克里金负责插值——当你手上只有有限个监测点的值时，用它估计未采样点的值，回答「两个站点之间的区域 AQI 大概是多少」。
这套方案的分工是这样的：先用 GWR 对每个观测点做回归预测，得到每个点的 AQI 预测值；再用克里金对这些预测值做空间插值，生成连续分布图。两段式的好处是，GWR 阶段保留了自变量对 AQI 的解释能力，克里金阶段补上了空间连续性。如果只用 GWR，你得到的是散点预测值，画不出连续曲面；如果只用克里金，你只是在插值已知的 AQI，没有利用气象因子做因果解释。
2.2 数据准备与变量筛选
原始数据里通常包含多个气象变量。这套方案里出现了 7 个自变量：temp（气温）、hpa（气压）、wet（湿度）、speed（风速）、dir（方向距离值）、mm、height（高度）。其中 mm 变量所有观测值都是 0，属于无效变量，直接剔除。这一步看起来简单，但实际项目中经常有人忘了检查常量列，导致回归矩阵奇异，报错信息还特别隐晦。
数据预处理的完整流程通常是：
import numpy as np import pandas as pd from scipy import stats # 读取原始数据，包含坐标、AQI 和气象变量 df = pd.read_csv('air_quality_data.csv') # 检查并剔除常数列（方差为 0 的变量） variances = df.var(numeric_only=True) constant_cols = variances[variances == 0].index.tolist() df = df.drop(columns=constant_cols) print(f"剔除常量列: {constant_cols}") # 提取坐标和变量 coords = df[['x', 'y']].values # 空间坐标 y = df['aqi'].values # 因变量 X = df[['temp', 'hpa', 'wet', 'speed', 'dir', 'height']].values # 自变量 # 标准化处理（GWR 对量纲敏感，建议标准化） X_mean, X_std = X.mean(axis=0), X.std(axis=0) X_norm = (X - X_mean) / X_std坐标列的选择取决于你的数据格式。如果是经纬度，直接用经纬度作为空间坐标；如果是投影坐标，用投影后的 x、y。标准化这一步不是必须的，但做了之后回归系数的量级更可比，解释起来方便。
2.3 GWR 模型校准：核函数、带宽与求解
GWR 的核心公式是：
$$\hat{\beta}(u_i, v_i) = (X^T W_i X)^{-1} X^T W_i y$$
其中 $W_i$ 是第 $i$ 个观测点的空间权重矩阵，对角线元素由核函数决定。这套方案选用 bi-square 核函数：
$$w_{ij} = \begin{cases} [1 - (d_{ij}/b)^2]^2 & d_{ij} < b \ 0 & d_{ij} \geq b \end{cases}$$
$d_{ij}$ 是点 $i$ 和点 $j$ 之间的距离，$b$ 是带宽。带宽的选择直接决定模型的平滑程度：带宽太大，退化成全局回归；带宽太小，每个点只用自己的数据拟合，过拟合。方案里用的是 CV（交叉验证）最小值法来搜索最优带宽，同时也可以用 AICc 作为替代指标。
from scipy.spatial.distance import cdist from scipy.optimize import minimize_scalar def gwr_bisquare(coords, X, y, bw): """GWR 核心求解，bi-square 核函数，自适应带宽""" n = len(y) betas = np.zeros((n, X.shape[1] + 1)) dist_matrix = cdist(coords, coords, metric='euclidean') for i in range(n): d_i = dist_matrix[i] # 自适应带宽：bw 为最近邻个数，取第 bw 个距离作为阈值 b = np.sort(d_i)[min(bw, n - 1)] w = np.where(d_i < b, (1 - (d_i / b) ** 2) ** 2, 0) W = np.diag(w) Xw = np.column_stack([np.ones(n), X]) # 加权最小二乘求解局部系数 XtWX = Xw.T @ W @ Xw XtWy = Xw.T @ W @ y try: betas[i] = np.linalg.solve(XtWX, XtWy) except np.linalg.LinAlgError: betas[i] = np.linalg.lstsq(XtWX, XtWy, rcond=None)[0] return betas def cv_score(bw, coords, X, y): """交叉验证：留一法计算 CV 值""" n = len(y) dist_matrix = cdist(coords, coords, metric='euclidean') residuals = np.zeros(n) for i in range(n): d_i = dist_matrix[i] b = np.sort(d_i)[min(int(bw), n - 1)] w = np.where(d_i < b, (1 - (d_i / b) ** 2) ** 2, 0) w[i] = 0 # 留一：排除自身 W = np.diag(w) Xw = np.column_stack([np.ones(n), X]) XtWX = Xw.T @ W @ Xw XtWy = Xw.T @ W @ y try: beta_i = np.linalg.solve(XtWX, XtWy) except np.linalg.LinAlgError: beta_i = np.linalg.lstsq(XtWX, XtWy, rcond=None)[0] residuals[i] = y[i] - Xw[i] @ beta_i return np.sqrt(np.mean(residuals ** 2)) # 搜索最优带宽（自适应带宽，范围为 10 到 n-1） result = minimize_scalar( lambda bw: cv_score(bw, coords, X_norm, y), bounds=(10, len(y) - 1), method='bounded' ) optimal_bw = int(result.x) print(f"最优带宽（最近邻数）: {optimal_bw}") # 用最优带宽拟合最终模型 betas = gwr_bisquare(coords, X_norm, y, optimal_bw)这段代码里有两个关键设计点。第一，自适应带宽的含义是「每个点取最近的 bw 个邻居作为有效范围」，而不是固定距离——这样在数据密集区带宽自动缩小，在稀疏区自动放大，比固定带宽更稳健。第二，CV 搜索用的是留一法，每次排除当前点自身，避免用自己预测自己导致 CV 值虚低。搜索范围从 10 开始，是因为邻居太少会导致局部矩阵接近奇异，求解不稳定。
2.4 克里金插值的实现要点
克里金插值的核心是变异函数。这套方案用的是高斯模型拟合半变异函数，拟合参数为：块金常数 3.88e+01，基台值 2.95e-02，变程 7.59。判断参数 Q1=0.63、Q2=1.12、cR=66.68，Q1 越接近 0 越好，Q2 越接近 1 越好，cR 越小越好——从结果看拟合效果一般，说明空间自相关性不是特别强，这在空气质量数据里其实常见，因为污染源分布本身就不均匀。
from pykrige.ok import OrdinaryKriging def kriging_predict(coords, values, target_coords): """普通克里金插值，高斯变异函数模型""" OK = OrdinaryKriging( coords[:, 0], coords[:, 1], values, variogram_model='gaussian', verbose=False, enable_plotting=False ) # 对目标点进行插值 z_pred, ss_pred = OK.execute('points', target_coords[:, 0], target_coords[:, 1]) return z_pred, ss_pred # 用 GWR 预测值作为克里金输入 gwr_pred = np.array([np.column_stack([np.ones(len(X_norm)), X_norm])[i] @ betas[i] for i in range(len(y))]) # 对第 23 期坐标点做插值 target_coords = df_time23[['x', 'y']].values z_pred, ss_pred = kriging_predict(coords, gwr_pred, target_coords)变异函数模型的选择有讲究：高斯模型适合平滑变化的变量，指数模型适合变化剧烈的变量，球状模型介于两者之间。空气质量数据通常用高斯或指数模型，具体选哪个看拟合参数。如果 Q1 偏大、Q2 偏离 1 较远，说明变异函数拟合不好，可以换模型再试。
3. 结果解读与显著性分析：怎么从系数图里读出真东西
3.1 OLS 与 GWR 的对比：R² 从 0.417 到 0.916
方案里给出的对比数据很直观：OLS 的 R² 是 0.417，GWR 的 R² 是 0.916；AICc 从 3067.762 降到 2567.705。R² 翻了一倍多，说明空间异质性确实显著——不同地点的 AQI 形成机制不一样，用一套全局系数去拟合，丢掉了大量局部信息。AICc 下降超过 500，在模型选择意义上也是压倒性优势。
但这里有个容易翻车的地方：GWR 的 R² 高不代表模型一定好。GWR 的有效参数个数（ENP）远大于 OLS，它用更多的自由度换来了更好的拟合。如果 ENP 接近样本量，那就是过拟合。通常的做法是同时报告 ENP 和 AICc，AICc 已经惩罚了复杂度，所以 AICc 低才是真的优。
3.2 回归系数的空间分布说明了什么
表 3 给出了各变量回归系数的描述统计。以 temp 为例，均值 -56.438，标准差 61.919，最小值 -224.324，最大值 6.342。标准差比均值还大，说明气温对 AQI 的影响在不同地点差异极大——有的地方气温升高显著降低 AQI，有的地方影响很弱甚至反向。这种信息在 OLS 里完全看不到，OLS 只会给你一个平均系数。
再看 hpa（气压），均值 -134.426，范围从 -417.69 到 137.722，跨越了正负。这意味着气压对 AQI 的影响在空间上存在方向性反转：南部区域气压升高降低 AQI，通州附近反而升高 AQI。如果做全局回归，正负效应互相抵消，系数可能接近 0，你会误以为气压不重要。
3.3 显著性检验与分区解读
对每个变量的局部系数做 t 检验，在 α=0.05 水平上标记显著区域，这是 GWR 分析的标准动作。方案里的结论可以归纳为：
| 变量 | 显著正向影响区域 | 显著负向影响区域 | 影响强度特征 |
|---|---|---|---|
| 截距 | 北部（怀柔、密云、平谷） | 通州周边 | 北部系数 300+，中南部接近 0 |
| temp | 无 | 除房山外大部分区域 | 南部影响大，北部影响小 |
| hpa | 通州 | 南部、海淀、朝阳、房山、大兴 | 北部负影响大于南部 |
| wet | 石景山、丰台、房山 | 北部（密云北部系数 -343） | 北部负影响最大 |
| speed | 东北部（密云南部、怀柔、平谷） | 石景山、延庆、朝阳、大兴 | 方向性分化明显 |
| dir | 延庆（较弱） | 朝阳、大兴、石景山、丰台、房山、海淀 | 南部负影响大 |
| height | 大兴、石景山、房山、海淀 | 密云及周边 | 南部负影响小于北部 |
这张表是整份方案的核心产出。它告诉你：不同区域治理空气污染要抓不同的变量。北部山区重点看湿度和高度，南部城区重点看气压和方向距离值。这种分区施策的依据，全局回归给不了。
4. 预测精度与避坑指南：均方误差 449 和 423 背后的门道
4.1 两套预测结果的对比
GWR 方案用第 23 期数据预测，均方误差 449.44；克里金方案用第 23 期坐标预测，均方误差 423.29。两个误差在同一量级，克里金略低。但注意，这两个误差不能直接比较优劣——GWR 的预测值来自回归模型外推，克里金的预测值来自空间插值，输入信息不同。GWR 用了自变量信息，克里金只用了 AQI 的空间自相关。如果自变量对 AQI 的解释力强，GWR 应该更好；如果空间自相关强，克里金更好。这里两者接近，说明两种信息都有贡献。
4.2 避坑与常见问题排查
现象一：GWR 求解时报「矩阵奇异」错误。原因通常是局部窗口内样本太少，或者自变量之间存在严重共线性。解决方法是增大带宽下限，或者先做 VIF 检验剔除共线性变量。我一般会在建模前跑一遍statsmodels的 VIF，超过 10 的变量直接拿掉。
现象二：CV 搜索带宽时结果在边界值震荡。原因可能是搜索范围设置不合理，或者数据本身空间结构不明显。解决方法是扩大搜索范围，同时检查坐标单位——如果经纬度没转投影坐标，距离计算会失真，带宽搜索自然不稳定。
现象三：克里金插值结果出现异常极值。原因通常是变异函数拟合失败，块金效应过大或变程设置不当。解决方法是换变异函数模型（高斯换指数或球状），或者手动指定变程参数。另外，输入数据里如果有离群点，先做异常值处理再插值。
现象四：GWR 的 R² 很高但预测误差很大。这是典型的过拟合信号。检查 ENP 是否接近样本量，如果是，说明带宽太小。通常的做法是限制带宽下限不低于样本量的 10%，或者改用 AICc 代替 CV 做带宽优化。
现象五：显著性检验结果与业务认知矛盾。比如某个公认重要的变量在所有区域都不显著。原因可能是变量量纲差异太大，标准化没做，导致 t 统计量计算失真。解决方法是统一标准化后再做检验，或者检查该变量是否存在大量缺失值被填充为 0。
5. 验收与答辩：怎么证明这套方案真的跑通了
验收这套方案，核心看三件事：模型能不能复现、结果能不能解释、图能不能画出来。答辩时老师最常问的三个问题是：「为什么用 GWR 不用 OLS」「带宽怎么选的」「克里金插值的精度怎么评估」。这三个问题对应三个验收点。
第一个验收点：模型对比。你需要准备一张 OLS 和 GWR 的对比表，包含 R²、AICc、ENP、残差平方和。R² 提升幅度和 AICc 下降幅度是最直观的证据。如果老师追问「R² 高是不是过拟合」，你就报 ENP 和 AICc，说明 AICc 已经惩罚了复杂度。
第二个验收点：带宽选择的合理性。准备一张 CV 值随带宽变化的曲线图，标出最优带宽的位置。如果曲线在最优值附近比较平坦，说明模型对带宽不敏感，结果稳健；如果曲线尖锐，说明带宽选择很关键，需要解释为什么选这个值。通常的做法是同时报告 CV 最优带宽和 AICc 最优带宽，如果两者接近，说服力更强。
第三个验收点：预测精度的评估。均方误差（MSE）或均方根误差（RMSE）是标配，但光报一个数不够。你需要准备预测值与真实值的散点图，看是否沿对角线分布；再准备残差的空间分布图，看残差是否有空间聚集——如果残差在某个区域系统性偏大，说明模型在该区域拟合不足，需要引入新的自变量或调整带宽。
答辩时还有一个高频问题：「你的结论有什么实际意义？」这时候不要泛泛说「对空气质量治理有参考价值」，要具体到区域和变量。比如「密云北部湿度对 AQI 的负影响系数达到 -343，说明该区域湿度升高会显著改善空气质量，但同一变量在石景山是正向影响，治理策略需要分区制定」。这种具体到系数和区域的回答，比任何套话都有说服力。
我自己的习惯是，每次交方案前强制走一遍「三图一表」检查：CV 曲线图、系数空间分布图、预测散点图，加上 OLS/GWR 对比表。这四样齐了，答辩基本不会被问倒。希望帮到你。
本文还有配套的精品资源，点击获取
