简介：这份资源是一篇完整的大学课程论文，主题为应用地理加权回归（GWR）模型与克里金插值法对空气质量指数进行预测，面向地理信息科学、空间统计与遥感方向的学生及研究者，帮助其理解空间异质性建模与地统计插值的完整流程。压缩包内仅含1个docx文档，约1.08MB，正文涵盖引言、研究方法、GWR案例分析、回归系数显著性分析、预测评估、克里金案例及总结讨论，并附有GWR预测程序、克里格插值程序与地图可视化三部分Python源代码，便于读者复现实验。论文详细推导了GWR模型形式、核函数与带宽优化（CV、AICc）过程，对比了GWR与OLS在捕捉空间非线性与局部效应上的差异，并利用克里金插值生成连续污染分布图。目前已有1256人学习下载，适合作为空间分析课程设计、论文写作与Python空间建模的实操参考。
1. 从一份课程设计说起：GWR 与克里金如何联手预测空气质量指数
去年帮一个做环境数据分析的朋友看代码，他扔过来一份课程设计，标题是“应用 GWR 模型和克里金法对空气质量指数进行预测”。我第一反应是：这题目听着像 GIS 专业课的期末作业，能跑通吗？结果翻完附录里的 Python 源码，发现它把地理加权回归、泛克里金插值、地图可视化串成了一条完整的链路，数据预处理、带宽优化、显著性检验、残差分析全都有。这不是那种只贴公式的论文，而是一份能直接复现的空间分析项目。如果你手头正好有空气质量监测数据，想从“全局回归”升级到“局部空间建模”，或者想看看克里金插值怎么和 GWR 配合做预测，这份资源值得拆开细读。它适合两类人：一是做 GIS 或环境统计的学生，需要一份带源码的完整案例；二是想快速上手空间回归的开发者，不想从零推导公式，只想把模型跑起来、把图出出来。
2. 项目定位与技术选型：为什么是 GWR 加克里金
2.1 从 OLS 到 GWR：空间非平稳性的解法
普通线性回归（OLS）假设自变量对因变量的影响在整个研究区域内是恒定的，但空气质量指数（AQI）受气温、气压、湿度、风速、风向、高度等因素影响，这些因素在不同区域的作用强度明显不同。比如风速对北部郊区 AQI 的正面影响可能很强，对南部城区却是负面影响。GWR 的核心思路就是把地理坐标嵌入回归参数，让每个观测点都有自己的一套回归系数。项目里用 bisquare 核函数做距离衰减加权，通过 CV 值最小化搜索最优自适应带宽，最终 GWR 的 R² 达到 0.916，而 OLS 只有 0.417，差距非常直观。
2.2 克里金插值的角色：填补空间空白
GWR 输出的是每个监测点上的局部回归系数和预测值，但监测点之间还有大量未采样位置。克里金插值依据变异函数提供的空间变异性信息，为未采样点提供局部最优线性无偏估计。项目里用的是泛克里金，因为数据存在未知的漂移成分，需要先用多项式拟合趋势项，再对残差做克里金插值。高斯模型被选来拟合半变异函数，块金常数、基台值、变程三个参数决定了插值的平滑程度和影响范围。
2.3 技术栈与依赖
项目源码是 Python 写的，核心依赖包括：
| 库名 | 用途 |
|---|---|
| numpy | 矩阵运算、线性代数求解 |
| pandas | 数据读取与预处理 |
| scipy | 优化算法（Golden Section 搜索带宽） |
| matplotlib | 残差图、系数分布图 |
| pykrige | 克里金插值计算 |
| geopandas | 地图可视化与空间数据操作 |
| statsmodels | OLS 对比与显著性检验 |
提示：pykrige 和 geopandas 的版本兼容性容易出问题，建议用虚拟环境固定版本，避免依赖冲突。
3. 核心模块拆解：从数据预处理到 GWR 系数估计
3.1 数据预处理与变量筛选
项目以 AQI 为因变量，初始选取 7 个自变量：temp（气温）、hpa（气压）、wet（湿度）、speed（风速）、dir（方向距离值）、mm（某个变量）、height（高度）。检查后发现 mm 变量所有值都是 0，直接剔除，最终保留 6 个自变量。数据按时间期数组织，用第 22 期做模型拟合，第 23 期做预测验证。
import pandas as pd import numpy as np # 读取监测点数据，包含坐标、AQI 和自变量 data = pd.read_csv('air_quality.csv') # 剔除全为 0 的 mm 变量 data = data.drop(columns=['mm']) # 提取第 22 期数据用于模型拟合 train_data = data[data['time'] == 22].copy() test_data = data[data['time'] == 23].copy() # 标准化自变量，避免量纲差异影响带宽搜索 feature_cols = ['temp', 'hpa', 'wet', 'speed', 'dir', 'height'] train_X = train_data[feature_cols].values train_y = train_data['aqi'].values coords = train_data[['x', 'y']].values # 空间坐标这段代码做了三件事：剔除无效变量、按期数切分训练与测试集、提取坐标矩阵。标准化那步我一般会加上，因为 temp 和 hpa 的数值范围差好几个数量级，不标准化会让距离计算偏向大数值变量。
3.2 空间权重矩阵与带宽优化
GWR 的核心是权重矩阵 W，每个观测点对回归点 i 的权重由 bisquare 核函数决定。带宽 bs 控制衰减速度，bs 越大，权重衰减越慢，模型越平滑；bs 越小，局部特征越突出，但容易过拟合。项目用 CV 值最小化来搜索最优带宽，具体用 Golden Section 方法迭代求解。
from scipy.optimize import minimize_scalar def bisquare_weight(distance, bandwidth): """bisquare 核函数：距离超过带宽的权重为 0""" w = np.zeros_like(distance) mask = distance < bandwidth w[mask] = (1 - (distance[mask] / bandwidth) ** 2) ** 2 return w def cv_score(bandwidth, coords, X, y): """交叉验证得分：留一法预测误差平方和""" n = len(y) errors = [] for i in range(n): # 计算其他点到 i 的距离 dist = np.sqrt(np.sum((coords - coords[i]) ** 2, axis=1)) w = bisquare_weight(dist, bandwidth) w[i] = 0 # 留一法排除自身 # 加权最小二乘求解局部系数 W = np.diag(w) XtWX = X.T @ W @ X XtWy = X.T @ W @ y beta = np.linalg.solve(XtWX, XtWy) y_pred = X[i] @ beta errors.append((y[i] - y_pred) ** 2) return np.sqrt(np.mean(errors)) # Golden Section 搜索最优带宽 result = minimize_scalar(cv_score, bounds=(0.1, 10.0), method='bounded', args=(coords, train_X, train_y)) optimal_bandwidth = result.xbisquare_weight 函数里，距离超过带宽的观测点权重直接置零，这保证了局部性。cv_score 用留一法逐个预测，每次排除自身点，避免信息泄漏。minimize_scalar 在 0.1 到 10.0 之间搜索，这个范围需要根据坐标单位调整，如果坐标是经纬度，范围要小得多。
3.3 局部系数估计与显著性检验
得到最优带宽后，对每个回归点 i 构造权重矩阵 W_i，用加权最小二乘求解局部系数。项目里还计算了每个系数的标准误和 t 值，用来判断哪些自变量在哪些区域显著。
def gwr_estimate(coords, X, y, bandwidth): """对每个观测点估计局部回归系数和标准误""" n, k = X.shape betas = np.zeros((n, k)) se = np.zeros((n, k)) for i in range(n): dist = np.sqrt(np.sum((coords - coords[i]) ** 2, axis=1)) w = bisquare_weight(dist, bandwidth) W = np.diag(w) # 加权最小二乘 XtWX = X.T @ W @ X XtWy = X.T @ W @ y beta = np.linalg.solve(XtWX, XtWy) betas[i] = beta # 计算标准误：sigma^2 * diag((X'WX)^-1 X'W W X (X'WX)^-1) residual = y - X @ beta sigma2 = np.sum(w * residual ** 2) / (np.sum(w) - 2 * np.trace( X @ np.linalg.inv(XtWX) @ X.T @ W)) C = np.linalg.inv(XtWX) @ X.T @ W var_beta = sigma2 * np.diag(C @ C.T) se[i] = np.sqrt(var_beta) return betas, se betas, se = gwr_estimate(coords, train_X, train_y, optimal_bandwidth) t_values = betas / se # t 统计量这段代码是 GWR 最核心的部分。betas 矩阵每一行对应一个观测点的局部系数，se 是对应标准误。t_values 大于 1.96 就说明在 0.05 水平上显著。项目里发现截距项在北部区域显著为正，temp 对 AQI 普遍负显著，hpa 在南部负显著，这些结论都来自这个 t 值矩阵。
3.4 预测与误差评估
用第 22 期拟合的模型去预测第 23 期的 AQI，需要把第 23 期的坐标和自变量代入，但局部系数是基于第 22 期观测点估计的。常见做法是：对第 23 期每个预测点，找到第 22 期中距离最近的若干观测点，用它们的局部系数做距离加权平均，或者直接用第 22 期对应位置的系数。
# 用第 22 期系数预测第 23 期 test_X = test_data[feature_cols].values test_coords = test_data[['x', 'y']].values test_y = test_data['aqi'].values predictions = [] for i, coord in enumerate(test_coords): # 找到第 22 期中距离最近的观测点 dist = np.sqrt(np.sum((coords - coord) ** 2, axis=1)) nearest_idx = np.argmin(dist) # 用最近点的局部系数预测 pred = test_X[i] @ betas[nearest_idx] predictions.append(pred) mse = np.mean((test_y - np.array(predictions)) ** 2) print(f'均方误差: {mse:.4f}')项目里报告的均方误差是 449.44，这个值不算小，说明 AQI 的波动性很强，单靠这 6 个自变量和空间位置还不够。但对比 OLS 的预测误差，GWR 已经有明显改善。
4. 克里金插值模块：泛克里金与半变异函数拟合
4.1 泛克里金原理与漂移项处理
普通克里金假设数据没有全局趋势，但空气质量数据往往存在大尺度的漂移，比如从北到南 AQI 逐渐升高。泛克里金通过引入多项式漂移项来处理这种非平稳性。项目里用一次多项式拟合趋势，然后对残差做克里金插值。
from pykrige.ok import OrdinaryKriging from pykrige.uk import UniversalKriging # 提取第 22 期数据 kriging_data = train_data[['x', 'y', 'aqi']].values x = kriging_data[:, 0] y = kriging_data[:, 1] z = kriging_data[:, 2] # 泛克里金：漂移项用一次多项式 UK = UniversalKriging(x, y, z, variogram_model='gaussian', drift_terms=['regional_linear']) z_pred, ss = UK.execute('grid', grid_x, grid_y)UniversalKriging 的 drift_terms 参数设为 regional_linear，表示用坐标的一次函数拟合趋势。variogram_model 选 gaussian，因为项目里高斯模型的拟合参数是：块金常数 3.88e+01，基台值 2.95e-02，变程 7.59e+00。块金常数相对基台值偏大，说明数据存在一定的随机噪声。
4.2 半变异函数拟合与模型判断
半变异函数刻画了空间相关性随距离的变化。项目用 Q1、Q2、cR 三个指标判断拟合效果：Q1 越接近 0 越好，Q2 越接近 1 越好，cR 越小越好。高斯模型的 Q1=0.63，Q2=1.12，cR=66.68，拟合效果一般。这说明空气质量数据的空间结构比较复杂，单一高斯模型可能不够。
# 手动计算半变异函数并拟合高斯模型 from scipy.optimize import curve_fit def gaussian_variogram(h, nugget, sill, range_): return nugget + sill * (1 - np.exp(-h**2 / range_**2)) # 计算实验半变异函数 lag, semivariance = [], [] for i in range(len(x)): for j in range(i+1, len(x)): h = np.sqrt((x[i]-x[j])**2 + (y[i]-y[j])**2) lag.append(h) semivariance.append(0.5 * (z[i]-z[j])**2) lag = np.array(lag) semivariance = np.array(semivariance) # 拟合高斯模型 popt, _ = curve_fit(gaussian_variogram, lag, semivariance, p0=[1.0, 1.0, 1.0]) nugget, sill, range_ = popt这段代码展示了半变异函数的计算和拟合过程。实际项目中可以直接用 pykrige 的 Variogram 类，但手动实现一遍有助于理解参数含义。nugget 是块金常数，反映随机误差；sill 是基台值，反映总变异；range_ 是变程，超过这个距离空间相关性就消失了。
4.3 克里金预测与残差分析
用第 23 期的坐标点做预测，项目报告的均方误差是 423.29，比 GWR 的 449.44 略好。残差图显示预测值在 38 到 66 之间波动，实际值在 46 到 65 之间，整体趋势吻合，但个别点偏差较大。
# 对第 23 期坐标做克里金预测 test_coords_kriging = test_data[['x', 'y']].values z_pred_kriging, ss_kriging = UK.execute('points', test_coords_kriging[:, 0], test_coords_kriging[:, 1]) # 计算均方误差 mse_kriging = np.mean((test_y - z_pred_kriging) ** 2) print(f'克里金均方误差: {mse_kriging:.4f}') # 残差图 import matplotlib.pyplot as plt residuals = test_y - z_pred_kriging plt.scatter(z_pred_kriging, residuals) plt.axhline(y=0, color='r', linestyle='--') plt.xlabel('预测值') plt.ylabel('残差') plt.title('克里金残差图') plt.show()残差图是判断模型是否系统性偏差的重要工具。如果残差随机分布在 0 附近，说明模型没有明显偏倚；如果呈现 U 形或倒 U 形，说明模型可能漏掉了非线性项。
5. 避坑与常见问题：跑通这份代码要注意什么
5.1 坐标单位与带宽范围不匹配
现象：CV 搜索出来的带宽要么是边界值 0.1，要么是 10.0，模型效果很差。 原因：坐标如果是经纬度（范围 -180 到 180），带宽搜索范围 0.1 到 10.0 太大，导致所有点权重都接近 1，退化成全局回归。 解决：先看坐标的量纲。经纬度的话，带宽范围设成 0.01 到 1.0；如果是投影坐标（米），范围可以设成 1000 到 50000。我一般会先画个距离分布直方图，取中位数附近作为搜索中心。
5.2 pykrige 的 drift_terms 参数报错
现象：调用 UniversalKriging 时提示 drift_terms 不支持 regional_linear。 原因：pykrige 版本差异，老版本用 'regional_linear'，新版本可能改成 'regional_linear' 或需要传入自定义函数。 解决：先pip show pykrige看版本，然后查对应版本的文档。如果不想折腾，直接用 OrdinaryKriging 加手动去趋势：先用线性回归拟合坐标和 AQI 的关系，对残差做普通克里金，最后把趋势加回去。
5.3 显著性检验的 t 值全为 NaN
现象：计算 t_values 时出现大量 NaN 或 Inf。 原因：标准误计算中 sigma2 为负，通常是因为分母sum(w) - 2*trace(...)小于等于 0，说明有效样本数太少或带宽太小。 解决：检查带宽是否过小，导致权重矩阵中非零元素太少。可以设一个最小带宽下限，比如坐标范围的 5%。另外，确保 X 矩阵没有完全共线的列，否则 XtWX 不可逆。
5.4 预测时用错系数矩阵
现象：预测均方误差极大，甚至超过 OLS。 原因：用第 23 期的坐标去匹配第 22 期的系数时，如果直接取最近点，可能匹配到空间上很远但编号相邻的点。 解决：不要用索引匹配，一定要用坐标距离匹配。我一般会算第 23 期每个点到第 22 期所有点的距离，取最近的 3 到 5 个点做距离加权平均，而不是只取一个最近点。
5.5 地图可视化坐标轴颠倒
现象：geopandas 画出来的地图上下颠倒或左右镜像。 原因：坐标参考系（CRS）设置错误，或者 x 和 y 列在读取时被交换。 解决：先print(data.head())确认 x 对应经度、y 对应纬度。geopandas 的set_crs要跟数据实际坐标系一致，WGS84 用 EPSG:4326，投影坐标用对应的 EPSG 代码。画图前用gdf.plot()快速检查形状。
6. 二次开发与改造思路：从课程设计到可复用工具
这份代码的骨架已经搭好了，但离一个能复用的空间分析工具还有距离。我一般会从三个方向改造。
第一，把 GWR 和克里金封装成类，统一接口。现在代码是脚本式的，变量散落在全局。可以定义一个SpatialAQIModel类，初始化时传入数据，提供fit_gwr()、fit_kriging()、predict()三个方法。这样换一份数据只需要改初始化参数，不用动核心逻辑。
class SpatialAQIModel: def __init__(self, data, feature_cols, target_col='aqi'): self.data = data self.feature_cols = feature_cols self.target_col = target_col self.gwr_betas = None self.kriging_model = None def fit_gwr(self, bandwidth=None): # 自动搜索带宽或使用指定值 if bandwidth is None: bandwidth = self._optimize_bandwidth() self.gwr_betas, self.gwr_se = gwr_estimate( self.coords, self.X, self.y, bandwidth) return self def fit_kriging(self, variogram_model='gaussian'): from pykrige.uk import UniversalKriging self.kriging_model = UniversalKriging( self.coords[:, 0], self.coords[:, 1], self.y, variogram_model=variogram_model, drift_terms=['regional_linear']) return self def predict(self, new_coords, new_X, method='gwr'): if method == 'gwr': # 距离加权平均局部系数 pass elif method == 'kriging': z_pred, _ = self.kriging_model.execute( 'points', new_coords[:, 0], new_coords[:, 1]) return z_pred第二，加交叉验证和模型对比。现在只用了第 22 期训练、第 23 期测试，样本量太小。可以改成 k 折交叉验证，每次留出一部分监测点做验证，计算 RMSE 和 R² 的均值和标准差。这样得到的模型评估更可靠。我一般会写一个cross_validate(model, data, k=5)函数，返回每个折的误差指标。
第三，把可视化做成可配置的。现在地图代码是硬编码的，可以改成传入column参数，自动根据系数值或预测值着色。用matplotlib的Normalize和cmap做连续色阶，加个 colorbar 标注单位。如果要做交互式地图，可以接folium或plotly，但那是另一个话题了。
注意：改造时不要动原始数据的坐标列名，很多空间库对列名有隐式依赖，改完容易出玄学 bug。
从那以后我每次拿到空间数据，都强制先跑一遍坐标范围检查和距离分布直方图，确认带宽搜索范围合理再往下走。这份课程设计的代码不算长，但把 GWR 和克里金的核心链路都串起来了，适合作为二次开发的起点。希望帮到你。
本文还有配套的精品资源，点击获取
