GWR+克里金联合建模：解决空气质量空间非平稳预测

简介：本资源是一份面向地理信息科学、环境科学及空间统计方向本科生与初阶研究者的课程论文实践资料，聚焦于利用地理加权回归（GWR）模型与克里金插值法协同预测空气质量指数（AQI），解决传统全局模型难以刻画空间异质性的问题。压缩包为1个1.08MB的docx文档，完整包含引言、研究方法（含GWR数学推导与带宽优化策略、克里金变异函数建模）、GWR与克里金双案例分析（含系数显著性解读、预测结果对比）、总结讨论，以及附录中可直接运行的Python源代码（GWR建模、克里金插值、地图可视化三部分）。内容结构严谨，公式推导清晰，代码注释详实，特别适合学习空间回归建模逻辑、理解局部参数估计思想、掌握GIS空间预测全流程的读者。目前已有1258人学习下载，是理论结合编程实践的典型教学范例。

1. 为什么用GWR+克里金做AQI预测，而不是直接上LSTM或XGBoost？

你手头有一批带坐标的空气质量监测点数据——比如北京16个区各一个站点，每小时记录AQI、温度、气压、湿度、风速、风向、海拔。你想回答的不是“明天全城平均AQI多少”，而是：“如果朝阳区某新建工地开始施工，它周边500米内居民暴露的AQI会升高多少？这个效应在通州和密云是否一样强？”

这时候扔进一个全局模型（比如OLS或XGBoost），它会告诉你“湿度每升1%，AQI平均降0.3”——但实际中，这个关系在北京南部可能是-0.8，在延庆却变成+0.2（地形抬升导致湿空气滞留）。空间非平稳性才是AQI建模真正的拦路虎。GWR不强行统一参数，而是给每个站点配一套“本地化系数”：朝阳站的湿度系数是-0.78，密云站却是+0.15，模型自己学出来。而克里金不关心“为什么”，只解决“没测的地方值多少”——它用已知点的空间自相关结构（变异函数），把离散监测点“织”成连续污染分布图，让决策者一眼看到污染热区在哪、扩散方向朝哪。

这套组合拳专治三类典型课设/毕设痛点：

• 数据少还带坐标：不像深度学习要上万条样本，GWR用20+站点+7个变量就能跑出可解释结果；
• 导师强调空间思维：回归系数能画成地图（比如temp系数图显示“北部降温对AQI改善更敏感”），答辩时比一串数字有说服力；
• 代码必须可复现：附录里的Python脚本不调黑盒API，核心是mgwr库的GWR类和pykrige的UniversalKriging，所有参数（核函数、带宽搜索策略、变异函数模型）都明文写死，你换自己数据改3行路径就能跑通。

这不是炫技，是用地理信息专业的“基本功”解决真实问题：当环保部门需要定位治理优先级时，他们要的不是RMSE降低0.5，而是“昌平回龙观片区风速系数显著为正，建议排查该区域夜间渣土车通行路线”。下面我们就从需求落地的角度，拆解这套方案怎么搭、怎么调、怎么避坑。

2. 总体架构设计：三层数据流与四个核心模块的职责划分

整个预测流程不是“先GWR再克里金”的简单串联，而是按数据生命周期分层处理：原始观测层 → 空间建模层 → 插值增强层 → 可视化输出层。每一层都有明确输入输出和容错机制，避免学生常犯的“代码跑通但结果无法解释”问题。

2.1 原始观测层：坐标标准化与变量清洗的硬性门槛

空气质量监测数据最常翻车的不是模型，而是数据本身。以论文中提到的北京16个站点为例，原始CSV可能包含：

• 坐标格式混乱：有的用WGS84经纬度（如116.4074,39.9042），有的用GCJ02偏移坐标，甚至混入UTM东距/北距；
• 变量缺失模式复杂：temp字段在冬季某站点连续72小时为空，但hpa有值；mm（降水量）全为0——这说明该站点无雨量计，必须剔除而非插补。

关键操作：

# 使用pyproj统一转为WGS84，强制校验坐标范围（北京应在115.5~117.5, 39.5~41.0） import pyproj wgs84 = pyproj.CRS("EPSG:4326") gcj02 = pyproj.CRS("EPSG:4490") # 注意：国内常用GCJ02对应EPSG:4490 transformer = pyproj.Transformer.from_crs(gcj02, wgs84, always_xy=True) # 对每个站点坐标执行转换，并添加校验 def validate_coords(df): valid_mask = ( (df['lon'] >= 115.5) & (df['lon'] <= 117.5) & (df['lat'] >= 39.5) & (df['lat'] <= 41.0) ) if not valid_mask.all(): print(f"警告：{(~valid_mask).sum()}个站点坐标越界，已标记为NaN") return df.assign(valid_coord=valid_mask) # mm列全零检测——直接删除该列，不参与任何建模 if df['mm'].nunique() == 1 and df['mm'].iloc[0] == 0: df = df.drop('mm', axis=1) print("已移除全零变量：mm")

提示：坐标转换必须用pyproj而非手动加减偏移量。国内GCJ02到WGS84的偏移是非线性的，手动公式在郊区误差可达500米，会导致GWR权重矩阵完全失真。

2.2 空间建模层：GWR参数选择的物理意义驱动逻辑

GWR不是调参游戏，每个参数背后都有环境科学依据。论文中选用bisquare核函数和CV最小化带宽，其设计逻辑如下：

实现时的关键约束：

from mgwr.gwr import GWR from mgwr.sel_bw import Sel_BW import numpy as np # 构造GWR输入：y为AQI，X为其他6个变量（已剔除mm），coords为(N,2)数组 y = df['aqi'].values.reshape(-1, 1) X = df[['temp','hpa','wet','speed','dir','height']].values coords = df[['lon','lat']].values # 强制使用自适应带宽（knn模式），避免固定距离在稀疏区失效 sel = Sel_BW(coords, y, X, kernel='bisquare') bw = sel.search(criterion='cv') # 使用CV而非AICc，因样本量小（N=16） # 实例化GWR模型，显式指定参数避免默认陷阱 gwr_model = GWR( coords, y, X, bw=bw, fixed=False, # 自适应带宽 kernel='bisquare', constant=True, # 包含截距项 spherical=False # 坐标已是平面投影，不启用球面距离 ) gwr_results = gwr_model.fit()

注意：spherical=False必须显式设置。mgwr默认开启球面距离计算，但当坐标是WGS84经纬度时，它会错误地将经纬度当作平面坐标算欧氏距离（导致北京东西向1度≈111km，南北向1度≈85km，严重扭曲权重）。正确做法是先用geopandas转为UTM投影，或直接设spherical=False并确保坐标已转为平面米制。

2.3 插值增强层：克里金法如何补足GWR的“空间盲区”

GWR只能预测已有站点位置的AQI，但环保规划需要知道“未布点区域”的污染水平。克里金在此不是替代GWR，而是用GWR残差的空间结构指导插值。论文中泛克里金（Universal Kriging）的妙处在于：它把GWR的预测值作为趋势面（trend surface），再对残差进行空间插值，最终结果 = GWR趋势 + 克里金残差修正。

为什么不用普通克里金？
普通克里金假设均值恒定，但AQI存在明显空间趋势（如城区普遍高于郊区）。泛克里金用多项式拟合趋势，再对去趋势后的残差建模，论文中pykrige的UniversalKriging正是实现此逻辑：

from pykrige.uk import UniversalKriging import numpy as np # 步骤1：用GWR预测所有站点的AQI（已得gwr_results） gwr_pred = gwr_results.predy.flatten() # 形状(N,) # 步骤2：计算残差 = 观测值 - GWR预测值 residuals = y.flatten() - gwr_pred # 步骤3：对残差进行泛克里金插值（趋势设为'linear'，即用坐标一次项拟合） UK = UniversalKriging( coords[:, 0], coords[:, 1], residuals, variogram_model='gaussian', # 论文中用高斯模型，因其平滑性适合AQI drift_terms=['regional_linear'], # 关键！用坐标线性项拟合残差趋势 verbose=False, enable_plotting=False ) # 步骤4：生成规则网格用于插值（例如0.01度间隔的北京网格） lon_grid, lat_grid = np.meshgrid( np.arange(115.5, 117.6, 0.01), np.arange(39.5, 41.1, 0.01) ) grid_points = np.column_stack([lon_grid.ravel(), lat_grid.ravel()]) # 步骤5：预测残差网格，叠加GWR趋势得到最终AQI场 k_residual, _ = UK.execute('points', grid_points[:, 0], grid_points[:, 1]) # 注意：此处需将GWR趋势外推到网格点——用IDW插值GWR系数场，再计算

设计说明：drift_terms=['regional_linear']是泛克里金的核心。它告诉模型：“残差不是随机噪声，而是随地理位置线性变化的”，这比假设残差均值为零更符合实际（例如GWR在山区预测偏高，残差系统性为负）。

2.4 可视化输出层：从系数图到污染热力图的四步渲染链

答辩时打动导师的不是代码，而是能讲出故事的地图。论文中图4.1残差图只是起点，完整输出应包含：

1. GWR系数空间分布图：每个变量（temp/hpa等）单独成图，用matplotlib.colors.LinearSegmentedColormap定制色带（如temp系数用蓝→白→红，表示负向→无影响→正向）；
2. 显著性掩膜图：叠加t检验p值<0.05的区域，用半透明黑色遮罩非显著区；
3. 克里金插值AQI场：用contourf绘制等值线，叠加监测点位置（红色三角形）；
4. 预测误差空间图：计算GWR+克里金联合预测值与实测值的绝对误差，标识误差>50的异常区（可能指示未纳入的污染源）。

关键代码片段（生成系数图）：

import matplotlib.pyplot as plt import matplotlib.colors as mcolors # 假设gwr_results.params为(N,7)数组，第0列是截距，第1列是temp系数 temp_coef = gwr_results.params[:, 1] # 创建自定义色带：深蓝(-200)→浅蓝(0)→浅红(100)→深红(300) colors = ["navy", "skyblue", "white", "lightcoral", "red"] n_bins = 256 cmap = mcolors.LinearSegmentedColormap.from_list("temp_coef", colors, N=n_bins) fig, ax = plt.subplots(figsize=(10, 8)) scatter = ax.scatter( coords[:, 0], coords[:, 1], c=temp_coef, cmap=cmap, s=120, edgecolors='black', linewidth=0.5, vmin=-200, vmax=300 ) ax.set_title("Temp系数空间分布（单位：℃→AQI）", fontsize=14) plt.colorbar(scatter, ax=ax, label="系数值", extend='both') # 添加显著性圈选（p<0.05的点用星号标注） t_val = gwr_results.t_stat[:, 1] # temp对应的t统计量 p_val = 2 * (1 - stats.t.cdf(np.abs(t_val), df=gwr_results.df_model)) sig_mask = p_val < 0.05 ax.scatter( coords[sig_mask, 0], coords[sig_mask, 1], marker='*', s=200, c='yellow', edgecolors='black', linewidth=1 )

3. GWR建模全流程：从数据加载到系数显著性检验的七步实操

GWR不是“fit一下就完事”，它的每一步输出都服务于后续解释。下面以论文中北京16站点数据为例，展开完整建模链路，重点标注哪些输出必须保留、哪些中间结果决定模型成败。

3.1 数据加载与结构校验：拒绝“能跑就行”的侥幸心理

学生最容易忽略的是数据结构校验。GWR要求：

• y必须是列向量（N×1），不能是pandas Series（会报ValueError: Expected 2D array）；
• X的列名不能含空格或特殊字符（mgwr内部用字符串匹配列名）；
• coords必须是numpy二维数组（N×2），不能是GeoDataFrame的geometry列。

安全加载模板：

import pandas as pd import numpy as np def safe_load_data(filepath): """严格校验数据结构，返回标准化的y, X, coords""" df = pd.read_csv(filepath) # 1. 强制转换y为列向量 if 'aqi' not in df.columns: raise ValueError("CSV中必须包含' aqi'列（注意空格）") y = df['aqi'].values.astype(float).reshape(-1, 1) # 2. 构建X：排除非数值列、全零列、含空值列 candidate_vars = ['temp','hpa','wet','speed','dir','height','mm'] X_cols = [col for col in candidate_vars if col in df.columns] X_df = df[X_cols].copy() # 移除全零列（如mm） zero_cols = X_df.columns[(X_df == 0).all()].tolist() if zero_cols: print(f"移除全零变量：{zero_cols}") X_df = X_df.drop(zero_cols, axis=1) # 检查缺失值：GWR不支持nan，必须提前处理 if X_df.isnull().values.any(): print("警告：X中存在缺失值，将用列均值填充") X_df = X_df.fillna(X_df.mean()) X = X_df.values.astype(float) # 3. 提取坐标并校验 if not all(col in df.columns for col in ['lon','lat']): raise ValueError("CSV中必须包含'lon'和'lat'列") coords = df[['lon','lat']].values.astype(float) # 4. 最终校验维度 n = len(y) assert X.shape[0] == n, f"X行数({X.shape[0]})≠y长度({n})" assert coords.shape[0] == n, f"coords行数({coords.shape[0]})≠y长度({n})" return y, X, coords, X_df.columns.tolist() # 调用 y, X, coords, var_names = safe_load_data("beijing_aqi.csv")

3.2 带宽搜索：CV值曲线比单个最优值更重要

论文中只给出最优带宽bw=4，但实际调试中，必须绘制CV值随带宽变化的曲线。原因：CV极小值可能出现在平坦谷底，此时带宽微小变动不影响结果；也可能在陡峭峰顶，此时bw=4和bw=5的系数差异巨大。

绘制CV曲线的代码：

import matplotlib.pyplot as plt # 搜索带宽范围：从2到min(10, N-1)，步长1 knn_range = list(range(2, min(11, len(y)))) cv_scores = [] for knn in knn_range: try: sel_temp = Sel_BW(coords, y, X, kernel='bisquare') bw_temp = sel_temp.search(bw_min=knn, bw_max=knn, criterion='cv') cv_scores.append(sel_temp.cv) except: cv_scores.append(np.nan) # 某些带宽下计算失败 # 绘图 plt.figure(figsize=(8, 5)) plt.plot(knn_range, cv_scores, 'bo-', label='CV Score') plt.axvline(x=knn_range[np.argmin(cv_scores)], color='r', linestyle='--', label=f'Optimal BW={knn_range[np.argmin(cv_scores)]}') plt.xlabel('Number of Nearest Neighbors (k)') plt.ylabel('CV Score') plt.title('CV vs Bandwidth Search') plt.legend() plt.grid(True) plt.show() # 输出诊断 opt_knn = knn_range[np.argmin(cv_scores)] print(f"CV最小值：{min(cv_scores):.4f}，对应k={opt_knn}") if len([x for x in cv_scores if abs(x - min(cv_scores)) < 0.01]) > 3: print("✓ 带宽选择稳健：CV谷底宽度>3个k值") else: print("⚠ 带宽敏感：建议用AICc二次验证")

血泪经验：当CV曲线在k=3和k=4处几乎持平（差值<0.001），选k=4更稳妥——它用更多邻居提升稳定性；若k=3的CV显著更低，但k=2和k=4都飙升，说明数据稀疏，需检查坐标是否聚集（如16个点中有12个在五环内）。

3.3 模型拟合与诊断：三个必须检查的GWR专属指标

GWR拟合后，gwr_results对象包含20+属性，但以下三项是答辩必问：

诊断代码：

from mgwr.gwr import GWR from mgwr.sel_bw import Sel_BW from spreg import OLS # 1. 计算OLS基准 ols = OLS(y, X) print(f"OLS AICc: {ols.aicc:.3f}, R²: {ols.r2:.3f}") # 2. GWR诊断 print(f"GWR ENP: {gwr_results.ENP:.2f} (OLS参数个数=7)") print(f"GWR AICc: {gwr_results.aicc:.3f} (比OLS低{ols.aicc-gwr_results.aicc:.1f})") print(f"GWR局部R²均值: {np.mean(gwr_results.localR2):.3f}") # 3. 绘制局部R²空间分布（识别薄弱区） plt.figure(figsize=(10, 4)) plt.subplot(1, 2, 1) plt.hist(gwr_results.localR2, bins=10, alpha=0.7) plt.xlabel('Local R²') plt.ylabel('Count') plt.title('Distribution of Local R²') plt.subplot(1, 2, 2) scatter = plt.scatter(coords[:, 0], coords[:, 1], c=gwr_results.localR2, cmap='viridis', s=80) plt.colorbar(scatter, label='Local R²') plt.title('Local R² Spatial Distribution') plt.show() # 4. 标识R²<0.2的站点（需重点分析） low_r2_mask = gwr_results.localR2 < 0.2 if low_r2_mask.any(): print(f"警告：{low_r2_mask.sum()}个站点局部R²<0.2，位置：") print(df[low_r2_mask][['lon','lat','aqi']].round(3))

3.4 系数显著性检验：t值地图比表格更有杀伤力

论文中表3只列出系数均值/标准差，但答辩时导师会问：“为什么朝阳的temp系数显著，而密云不显著？” 这需要t值空间可视化。

t值计算与绘图：

from scipy import stats # t_stat形状为(N, k+1)，k为变量数，+1为截距 t_stats = gwr_results.t_stat # (16, 7) # 对每个变量绘制t值地图（以temp为例，索引为1） temp_t = t_stats[:, 1] p_values = 2 * (1 - stats.t.cdf(np.abs(temp_t), df=gwr_results.df_model)) # 创建显著性掩膜 sig_mask = p_values < 0.05 # 绘图：t值大小 + 显著性标记 fig, ax = plt.subplots(figsize=(10, 8)) scatter = ax.scatter( coords[:, 0], coords[:, 1], c=temp_t, cmap='coolwarm', s=100, vmin=-4, vmax=4, alpha=0.7 ) ax.scatter( coords[sig_mask, 0], coords[sig_mask, 1], marker='o', s=200, facecolors='none', edgecolors='black', linewidth=2 ) plt.colorbar(scatter, ax=ax, label='t-statistic') ax.set_title("Temp系数t值空间分布（*表示p<0.05）", fontsize=14) plt.show() # 输出显著性统计 print(f"Temp系数：{sig_mask.sum()}/{len(sig_mask)}个站点显著（p<0.05）") print(f"其中t>2的站点：{np.sum(temp_t > 2)}/{len(temp_t)}")

玄学提示：t值分布常呈“中心高、边缘低”。若北京东南部（通州、大兴）t值普遍>3，而西北部（延庆、怀柔）t值接近0，暗示GWR在城区拟合更好——这很正常，因为城区监测点密度高，带宽内邻居多，估计更稳定。

3.5 预测与误差分析：用GWR残差指导克里金参数选择

GWR预测不是终点，而是克里金的起点。论文中表4的MSE=449看似合理，但需检查残差是否满足克里金前提：

1. 残差空间自相关：用Moran's I检验，I>0.3才适合克里金；
2. 残差分布正态性：Shapiro-Wilk检验p>0.05；
3. 残差无趋势：残差vs坐标散点图应无明显斜线。

残差诊断代码：

from libpysal.weights import Queen from esda.moran import Moran from scipy.stats import shapiro # 计算残差 residuals = y.flatten() - gwr_results.predy.flatten() # 1. 空间自相关检验 w = Queen.from_dataframe(df) # 基于邻接关系构建空间权重 moran = Moran(residuals, w) print(f"Moran's I: {moran.I:.3f} (p={moran.p_sim:.3f})") if moran.I < 0.2: print("⚠ 残差空间自相关弱，克里金增益有限") # 2. 正态性检验 shap_p = shapiro(residuals).pvalue print(f"Shapiro-Wilk p-value: {shap_p:.3f}") if shap_p < 0.05: print("⚠ 残差非正态，考虑对AQI取log变换后重跑GWR") # 3. 绘制残差空间分布 plt.figure(figsize=(12, 4)) plt.subplot(1, 3, 1) plt.hist(residuals, bins=10, alpha=0.7) plt.title('Residuals Histogram') plt.subplot(1, 3, 2) plt.scatter(coords[:, 0], residuals, alpha=0.6) plt.xlabel('Longitude') plt.ylabel('Residual') plt.title('Residual vs Longitude') plt.subplot(1, 3, 3) plt.scatter(coords[:, 1], residuals, alpha=0.6) plt.xlabel('Latitude') plt.ylabel('Residual') plt.title('Residual vs Latitude') plt.show()

4. 克里金插值实战：从变异函数拟合到网格预测的五步闭环

克里金不是“选个模型点确定”，其核心是变异函数（Variogram）的物理合理性。论文中表5的高斯模型参数（块金=38.8，基台=0.029，变程=7.59）看似枯燥，实则暗含北京地形线索：变程7.59（单位？）若为经纬度，则约840km，远超北京范围——这说明论文中坐标单位是度，但变异函数拟合时应转为米。下面还原真实操作链。

4.1 坐标投影转换：为什么必须用UTM而非WGS84

变异函数的变程（range）单位必须是米，否则：

• WGS84下1度≈111km，变程7.59度=842km，意味着任意两点间相关性都接近1，克里金退化为均值插值；
• UTM下坐标单位为米，变程7590米=7.6km，符合城市尺度污染扩散认知。

正确投影代码：

import geopandas as gpd from shapely.geometry import Point # 创建GeoDataFrame geometry = [Point(xy) for xy in zip(df['lon'], df['lat'])] gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") # 投影到北京UTM带（WGS84 / UTM zone 50N, EPSG:32650） gdf_utm = gdf.to_crs(epsg=32650) coords_utm = np.column_stack([gdf_utm.geometry.x, gdf_utm.geometry.y]) print(f"投影后坐标单位：米，范围 x:{gdf_utm.geometry.x.min():.0f}-{gdf_utm.geometry.x.max():.0f}") # 用UTM坐标计算距离矩阵（供变异函数拟合） from sklearn.metrics import pairwise_distances dist_matrix = pairwise_distances(coords_utm, metric='euclidean')

4.2 变异函数拟合：高斯模型的选择依据与Q1/Q2/cR解读

论文中表6的Q1=0.63、Q2=1.12、cR=6.68，其含义是：

拟合与诊断代码：

from pykrige.variogram_models import gaussian_variogram_model from scipy.optimize import curve_fit # 计算实验变异函数 def calc_experimental_variogram(coords, z, max_lag=5000, n_lags=10): """coords: (N,2)米制坐标，z: (N,)残差值""" from sklearn.metrics import pairwise_distances dist_matrix = pairwise_distances(coords, metric='euclidean') np.fill_diagonal(dist_matrix, np.inf) # 排除自身 lags = np.linspace(0, max_lag, n_lags+1)[1:] gamma_vals = [] for i in range(len(lags)-1): mask = (dist_matrix >= lags[i]) & (dist_matrix < lags[i+1]) if mask.sum() > 0: gamma = 0.5 * np.mean((z[:, None] - z[None, :])**2 * mask) gamma_vals.append(gamma) else: gamma_vals.append(np.nan) return lags[:-1], np.array(gamma_vals) # 计算实验变异函数 exp_lags, exp_gamma = calc_experimental_variogram(coords_utm, residuals) # 拟合高斯模型：gamma(h) = nugget + sill * (1 - exp(-(h/range)^2)) def gaussian_model(h, nugget, sill, range_val): return nugget + sill * (1 - np.exp(-(h/range_val)**2)) # 初始参数：块金≈残差方差的10%，基台≈残差方差，变程≈最大距离的1/3 p0 = [np.var(residuals)*0.1, np.var(residuals), np.max(exp_lags)/3] popt, _ = curve_fit(gaussian_model, exp_lags, exp_gamma, p0=p0) # 计算Q1/Q2/cR fitted_gamma = gaussian_model(exp_lags, *popt) Q1 = np.median(np.abs(exp_gamma - fitted_gamma)) Q2 = np.std(exp_gamma - fitted_gamma) cR = np.sum(((exp_gamma - fitted_gamma) / (0.01 * fitted_gamma + 0.001))**2) print(f"高斯模型拟合参数：块金={popt[0]:.2f}, 基台={popt[1]:.2f}, 变程={popt[2]:.0f}m") print(f"Q1={Q1:.2f}, Q2={Q2:.2f}, cR={cR:.2f}")

4.3 泛克里金执行：drift_terms的三种模式对比

泛克里金的drift_terms参数决定趋势面形式，论文中用'regional_linear'，但需根据数据选择：

执行泛克里金：

from pykrige.uk import UniversalKriging # 使用UTM坐标和拟合的变异函数参数 UK = UniversalKriging( coords_utm[:, 0], coords_utm[:, 1], residuals, variogram_model='gaussian', variogram_parameters={'sill': popt[1], 'range': popt[2], 'nugget': popt[0]}, drift_terms=['regional_linear'], verbose=False ) # 生成预测网格（UTM坐标） x_grid = np.linspace(coords_utm[:, 0].min(), coords_utm[:, 0].max(), 100) y_grid = np.linspace(coords_utm[:, 1].min(), coords_utm[:, 1].max(), 100) xx, yy = np.meshgrid(x_grid, y_grid) grid_points = np.column_stack([xx.ravel(), yy.ravel()]) # 执行插值 k_residual, ss = UK.execute('grid', x_grid, y_grid) # ss为标准差网格

4.4 结果融合：GWR趋势 + 克里金残差的物理叠加逻辑

最终AQI场 = GWR预测值（在监测点） + 克里金残差（在网格点）。但GWR只给出监测点预测，需将其空间插值到网格，再叠加克里金残差：

from scipy.interpolate import griddata # 步骤1：将GWR预测值插值到网格（用IDW，因GWR本身是空间加权） gwr_grid = griddata( coords_utm, gwr_results.predy.flatten(), (xx, yy), method='linear' # linear比cubic更稳定 ) # 步骤2：叠加克里金残差（注意：k_residual是(100,100)数组） aqi_final = gwr_grid + <p> <a href="https://download.csdn.net/download/qq_37041207/12898257" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>