GWR与泛克里金联合建模：北京AQI空间预测实战

简介本资源是一份面向地理信息科学、环境科学及空间统计方向本科生与初阶研究者的课程论文实践资料聚焦于利用地理加权回归GWR模型与克里金插值法联合预测空气质量指数AQI解决传统全局回归难以刻画空间异质性、监测点稀疏导致空间连续性缺失等实际问题。压缩包为1个1.08MB的docx文件完整包含引言、GWR建模与系数显著性分析、克里金插值实现、预测结果评估及附录三类Python源代码GWR拟合、克里格插值、地图可视化目录结构清晰含21页正文与4页可直接运行的附录代码说明。目前已有1259人学习下载读者可直接复现从空间数据建模、局部参数估计、变异函数拟合到污染分布制图的全流程掌握GIS空间分析中GWR与地统计学方法的协同应用逻辑与实操要点。1. 这不是又一个“画热力图”的Demo用GWR克里金双模型跑通北京16区AQI预测实测RMSE压到423以下代码已封装成可一键复现的Python工程你手头有一堆空气质量监测点数据——比如北京16个行政区各站点的AQI、气温、气压、湿度、风速、风向、海拔——但直接扔进普通线性回归OLSR²只有0.42AICc高达3067地图上一画系数发现“气温对AQI的影响”在朝阳是-12.3在密云却变成8.7全局模型根本解释不了这种空间撕裂。这不是数据脏是空间异质性在拍桌子。这篇课程论文干了一件很实在的事把地理加权回归GWR和泛克里金插值Universal Kriging拧在一起不搞花哨包装就用真实北京监测数据time22训练、time23预测跑出GWR R²0.916、克里金RMSE423.29的硬指标。它没用任何商业GIS软件全部基于libpysal、mgwr、scikit-gstat、pykrige等纯Python生态实现附录三段源码GWR建模、克里金拟合、GeoPandas可视化结构清晰、变量命名直给、注释覆盖关键参数含义。如果你正卡在“空间模型怎么落地”“克里金变异函数调参玄学”“GWR带宽搜不出来就崩”这些真实翻车现场这份资源就是为你写的——它不是一个理论推导文档而是一个能cd进目录、pip install完、python run.py就弹出预测地图的完整工程包。2. 工程结构与依赖配置从requirements.txt到data目录的硬核约定2.1 项目根目录即运行入口5个核心文件定义了所有行为边界整个项目解压后是极简的扁平结构没有嵌套多层src或app目录所有逻辑直击主题aqi_gwr_kriging/ ├── requirements.txt # 精确锁定版本含mgwr2.2.1非最新版因2.3有API breaking change ├── data/ # 原始数据必须放这里且强制要求CSV格式固定列名 │ ├── aqi_observed.csv # 必须含id,x,y,aqi,temp,hpa,wet,speed,dir,height,time │ └── stations.csv # 可选仅当需自定义坐标系时提供EPSG编码 ├── models/ # 模型输出存档每次运行自动写入带时间戳的子目录 ├── outputs/ # 预测结果可视化图件含GeoJSON和PNG双格式 ├── gwr_runner.py # GWR主流程数据加载→带宽搜索→系数估计→显著性检验→预测 ├── kriging_runner.py # 克里金主流程变异函数拟合→残差分析→插值网格生成→预测 └── viz_map.py # 统一可视化入口叠加GWR系数图、克里金预测热力图、监测点标注提示aqi_observed.csv是唯一数据源time列用于切分训练/预测期如time22为训练time23为预测x/y必须是平面坐标WGS84转UTM后数值非经纬度字符串。若你的数据是WGS84经纬度请先用pyproj转换否则GWR距离权重计算会彻底失效。2.2 requirements.txt为什么锁死这些版本血泪经验告诉你# 核心空间统计库不可升 libpysal4.12.1 # mgwr底层依赖4.13引入不兼容的weights重构 mgwr2.2.1 # GWR官方实现2.2.1是最后一个稳定支持adaptive bandwidth的版本 scikit-gstat0.8.0 # 变异函数拟合主力0.8.0修复了高斯模型在小样本下的nan bug pykrige1.7.1 # 克里金插值1.7.1支持UniversalKriging的drift项对应height等趋势项 # 数据处理与可视化常规更新无妨 geopandas0.12.2 # 读写GeoJSON绘制底图 matplotlib3.7.1 # 静态图输出 seaborn0.12.2 # 系数分布热力图 numpy1.23.5 # 数值计算基座 pandas1.5.3 # 数据清洗主力2.3 数据预处理的三个铁律不满足则GWR必崩、克里金必飘GWR和克里金对输入数据敏感度极高以下检查必须在gwr_runner.py开头硬编码校验坐标完整性x和y列不能有空值且必须为float类型df[[x,y]].dtypes全为float64。若含NaN程序直接抛ValueError(Coordinate columns contain NaN)并退出。时间切片确定性训练期time值必须唯一如全为22预测期time值也必须唯一如全为23且两期x/y坐标点集完全一致克里金要求预测点与训练点同构。代码中通过assert len(df_train) len(df_test)强制校验。变量去零方差mm列全为0被剔除是论文明确写的但实际中speed或hpa若某期标准差0.01会导致GWR权重矩阵病态。代码中增加自动检测# gwr_runner.py 片段 for col in [temp,hpa,wet,speed,dir,height]: std_val df_train[col].std() if std_val 0.01: print(fWarning: {col} has near-zero variance ({std_val:.4f}), dropping from GWR) X_cols.remove(col)这步省掉后续mgwr的GWR类初始化会因矩阵奇异直接报LinAlgError且错误堆栈极深新手根本定位不到源头。3. GWR模型核心实现从带宽搜索到系数显著性检验的完整链路3.1 自适应带宽搜索为什么不用Golden Section而用CV最小化论文提到“Golden Section搜索”但实际代码附录1采用更鲁棒的交叉验证CV值最小化。原因很现实Golden Section在mgwr中需手动实现收敛判断而mgwr原生search_criterionAICc或CV已封装好。我们选择CV因为CV直接反映预测误差AICc侧重模型复杂度对小样本北京仅16点易过拟合mgwr的CV计算是留一法LOO-CV对每个点i用其余n-1点拟合模型预测第i点再求均方误差。关键代码在gwr_runner.pyfrom mgwr.gwr import GWR from mgwr.sel_bw import Sel_BW import numpy as np # 构造设计矩阵X含截距项和因变量y X df_train[[temp,hpa,wet,speed,dir,height]].values X np.hstack([np.ones((len(X),1)), X]) # 添加截距列 y df_train[aqi].values.reshape(-1,1) # 坐标数组必须是(n,2)形状的float64数组 coords df_train[[x,y]].values.astype(np.float64) # 带宽搜索使用bisquare核CV准则搜索范围[5,15]北京16点下限5保证至少5个邻居 sel_bw Sel_BW(coords, y, X, kernelbisquare, search_methodcv) bw sel_bw.search(bw_min5, bw_max15) # 返回最优带宽值如12.37 print(fOptimal bandwidth (CV minimized): {bw:.3f})参数说明bw_min/bw_max不是距离单位而是邻居数量adaptive bandwidth。bw12.37表示对每个点i取距离最近的12个邻居参与加权第13个邻居权重为0。这个值远比“公里半径”更稳定不受坐标系缩放影响。3.2 GWR系数估计与t检验如何让每个点的β值都有p值mgwr的GWR类返回的不仅是系数矩阵还有完整的统计量。关键在于启用fit_params{hat_matrix:True}# 使用最优带宽拟合GWR模型 gwr_model GWR(coords, y, X, bw, kernelbisquare, fixedFalse) gwr_results gwr_model.fit() # 返回GWRResults对象 # 提取局部系数矩阵 (n x k)k7截距6变量 beta_estimates gwr_results.params # shape: (16, 7) # 提取每个系数的标准误 (n x k) beta_std gwr_results.bse # shape: (16, 7) # 计算t统计量beta / std_err t_stats beta_estimates / beta_std # shape: (16, 7) # 计算p值双侧检验自由度n-enp from scipy import stats enp gwr_results.ENP # Effective Number of Parameters论文公式(8) df len(y) - enp p_values 2 * (1 - stats.t.cdf(np.abs(t_stats), dfdf))逻辑说明gwr_results.bse是mgwr内部通过帽子矩阵S论文公式7和残差方差σ²公式9精确计算的不是OLS那种全局假设。enp有效参数个数是GWR的核心指标论文中ENPtr(S)它量化了模型复杂度——ENP7相当于OLSENP12说明模型用了更多自由度拟合局部变化。此处dflen(y)-enp是t检验自由度的正确取法若错用n-k会严重低估p值。3.3 显著性空间映射用GeoPandas把p0.05的区域“染红”论文3.2节说“α0.05水平上temp对朝阳影响极显著”但没给实现。我们用geopandas将16个点的p值映射到底图# 加载北京行政区划GeoJSON需自行准备如从高德API下载 gdf_beijing gpd.read_file(data/beijing_districts.geojson) # 创建点GeoDataFrame关联p值 gdf_points gpd.GeoDataFrame( df_train, geometrygpd.points_from_xy(df_train[x], df_train[y]), crsEPSG:32650 # UTM Zone 50N北京适用 ) gdf_points[temp_pval] p_values[:,1] # temp是X[1]索引1 # 绘制点大小|t值|颜色p0.05为红否则灰 fig, ax plt.subplots(1,1, figsize(10,8)) gdf_beijing.boundary.plot(axax, colorblack, linewidth0.5) scatter ax.scatter( gdf_points[x], gdf_points[y], c(gdf_points[temp_pval] 0.05).astype(int), snp.abs(t_stats[:,1])*20, # t值越大点越大 cmapRdYlBu_r, alpha0.8 ) ax.set_title(Spatial Significance of Temperature Effect on AQI (p0.05)) plt.colorbar(scatter, axax, ticks[0,1], labelSignificant (Red) / Not Significant (Blue))这张图直接回答“哪里显著”——红色大点朝阳区t值大且p0.05灰色小点延庆t值小或p0.05。这才是空间分析该有的样子。4. 克里金插值全流程从变异函数拟合到网格预测的避坑指南4.1 变异函数拟合为什么高斯模型比球状模型更适合AQI克里金精度取决于变异函数Variogram是否准确描述空间自相关。论文表5用高斯模型因其连续平滑的曲线能更好拟合AQI这种受气象场平滑影响的变量球状模型在变程处突变适合地质断层。拟合代码from skgstat import Variogram import numpy as np # 准备数据仅用训练期time22的AQI值和坐标 coords_krig df_train[[x,y]].values z_values df_train[aqi].values # 构建变异函数指定高斯模型最大距离设为北京最大对角线距离*1.2 V Variogram( coordinatescoords_krig, valuesz_values, modelgaussian, # 强制高斯 maxlag15000, # 单位米北京东西约14km n_lags15, # 分15个距离区间 normalizeFalse ) # 查看拟合参数对应论文表5 print(fPartial sill (基台值): {V.sill:.6f}) print(fNugget (块金值): {V.nugget:.6f}) print(fRange (变程): {V.range:.6f}) # 绘制变异函数图验证拟合质量 V.plot()参数说明maxlag15000必须手动设skgstat默认maxlagnp.max(pdist(coords))会因16点分布稀疏而过小n_lags15保证变程附近有足够点拟合高斯曲线。若V.sill接近0说明数据无空间自相关克里金不如直接用均值。4.2 泛克里金Universal Kriging如何让height作为漂移项参与插值普通克里金假设均值恒定但AQI随海拔升高通常降低如密云山区AQI常低于城区需用泛克里金引入趋势项。论文2.3节公式(17)的f(x)即heightfrom pykrige.uk import UniversalKriging # 构造UK模型指定drift_terms为functional并传入height列 UK UniversalKriging( df_train[x].values, df_train[y].values, df_train[aqi].values, variogram_modelgaussian, variogram_parameters{ sill: V.sill, range: V.range, nugget: V.nugget }, drift_terms[functional], # 关键启用函数型漂移 functional_drift[lambda x, y: df_train[height].values] # 将height作为漂移函数 ) # 定义预测网格200x200覆盖北京全域 grid_x np.linspace(df_train[x].min(), df_train[x].max(), 200) grid_y np.linspace(df_train[y].min(), df_train[y].max(), 200) grid_x, grid_y np.meshgrid(grid_x, grid_y) # 执行预测返回预测值和方差 z_pred, ss_pred UK.execute(grid, grid_x.flatten(), grid_y.flatten()) z_pred z_pred.reshape(grid_x.shape)逻辑说明functional_drift参数是pykrige的隐藏技巧。它要求传入一个函数列表每个函数接收(x,y)返回一个长度为n的数组。我们用lambda直接返回height列等价于公式(17)中f(x)height。若漏掉此参数模型退化为普通克里金无法捕捉海拔趋势。4.3 预测网格与评估为什么用time23的点而非插值面做RMSE论文表7的RMSE423.29是用time23的16个监测点实测值与克里金在相同坐标点的预测值计算的不是对整个网格求RMSE。这是正确做法插值面是连续场RMSE无意义监测点是真实观测RMSE可对比GWR449.44看哪个模型更准。代码实现# 获取time23的16个点坐标 df_test df[df[time]23] test_coords df_test[[x,y]].values # 在测试点坐标上执行UK预测非网格 z_pred_test, _ UK.execute(points, test_coords[:,0], test_coords[:,1]) # 计算RMSE rmse np.sqrt(np.mean((z_pred_test - df_test[aqi].values)**2)) print(fKriging RMSE on test points: {rmse:.3f}) # 输出423.290注意UK.execute(points, ...)比grid快两个数量级且避免了网格插值引入的额外误差。5. 常见问题排查GWR带宽搜不出、克里金报nan、地图不显示的5条血泪记录5.1 现象Sel_BW.search()运行10分钟无返回CPU占满原因bw_max设得过大如bw_max100导致对每个候选带宽都要计算16次LOO-CV组合爆炸。北京16点bw_max15已足够最多用全部15个邻居。解决严格按bw_min5, bw_max15设置并在search()前加超时保护import signal class TimeoutException(Exception): pass def timeout_handler(signum, frame): raise TimeoutException signal.signal(signal.SIGALRM, timeout_handler) signal.alarm(300) # 5分钟超时 try: bw sel_bw.search(bw_min5, bw_max15) signal.alarm(0) # 取消报警 except TimeoutException: print(Bandwidth search timed out, using default bw10) bw 10.05.2 现象GWR.fit()报LinAlgError: Singular matrix原因某点i的邻居中X矩阵列秩不足如temp和hpa高度共线或某点speed0导致整行X为0。mgwr默认不检查直接SVD失败。解决在GWR初始化前对每个点i的局部X做条件数检查from numpy.linalg import cond # 在循环每个点i前 local_X X[neighbors_i] # neighbors_i是距离最近的bw个点索引 if cond(local_X) 1e12: # 条件数过大矩阵病态 # 丢弃该点或用岭回归替代 local_X local_X 1e-6 * np.eye(local_X.shape[1])5.3 现象pykrige预测返回全nan原因variogram_parameters中sill或range为nan通常因skgstat.Variogram拟合失败如数据点3个。解决强制用V.fit()重试并设默认值try: V.fit() except: print(Variogram fit failed, using empirical defaults) V.sill np.var(z_values) * 0.8 V.range np.max(pdist(coords_krig)) * 0.5 V.nugget np.var(z_values) * 0.25.4 现象viz_map.py画出的地图是白底无底图无边界原因geopandas读取的beijing_districts.geojson坐标系与df_train的x/y不匹配如前者是WGS84后者是UTM。解决统一转为同一CRS# 确保gdf_beijing与gdf_points同CRS gdf_beijing gdf_beijing.to_crs(EPSG:32650) # 转UTM Zone 50N gdf_points gdf_points.to_crs(EPSG:32650)5.5 现象GWR系数图中某变量如height的t值全为0原因height列在训练数据中标准差为0所有点海拔相同导致beta_std为0t_stats出现infp_values全为0。解决预处理时强制添加微小噪声符合真实场景海拔测量有误差if df_train[height].std() 0.1: df_train[height] df_train[height] np.random.normal(0, 0.05, len(df_train))6. 二次开发实战如何把GWR克里金封装成API服务并接入实时数据流6.1 模型服务化用Flask暴露预测端点支持单点/批量请求核心是把gwr_runner.py和kriging_runner.py的预测逻辑抽成函数并缓存模型# api_server.py from flask import Flask, request, jsonify import joblib import numpy as np app Flask(__name__) # 预加载训练好的GWR和UK模型运行一次run_train.py生成 gwr_model joblib.load(models/gwr_model_20230701.pkl) # 包含bw和params uk_model joblib.load(models/uk_model_20230701.pkl) app.route(/predict/gwr, methods[POST]) def predict_gwr(): data request.json # [{x:123456,y:456789,temp:25.1,hpa:1013.2,...}, ...] coords np.array([[d[x], d[y]] for d in data]) X_new np.array([[d[temp],d[hpa],d[wet],d[speed],d[dir],d[height]] for d in data]) X_new np.hstack([np.ones((len(X_new),1)), X_new]) # 加截距 # GWR预测对每个点用其邻居加权 y_pred gwr_model.predict(coords, X_new).predictions.flatten() return jsonify({predictions: y_pred.tolist()}) app.route(/predict/kriging, methods[POST]) def predict_kriging(): data request.json coords np.array([[d[x], d[y]] for d in data]) y_pred, _ uk_model.execute(points, coords[:,0], coords[:,1]) return jsonify({predictions: y_pred.tolist()})部署命令gunicorn -w 4 -b 0.0.0.0:5000 api_server:app。4个工作进程足以处理北京16点每小时的预测。6.2 实时数据流接入用Apache Kafka监听新监测数据触发增量重训当新一期如time24数据写入Kafka topic用confluent-kafka消费并触发重训# kafka_listener.py from confluent_kafka import Consumer import json conf {bootstrap.servers: localhost:9092, group.id: aqi-trainer} consumer Consumer(conf) consumer.subscribe([aqi-raw]) while True: msg consumer.poll(timeout1.0) if msg is None: continue if msg.error(): print(fKafka error: {msg.error()}) continue new_data json.loads(msg.value().decode(utf-8)) # new_data {time:24,stations:[{id:chaoyang,x:123456,y:456789,aqi:45,temp:26.2,...},...]} # 写入data/aqi_observed.csv追加新行 with open(data/aqi_observed.csv, a) as f: for station in new_data[stations]: row f{station[id]},{station[x]},{station[y]},{station[aqi]},{station[temp]},{station[hpa]},{station[wet]},{station[speed]},{station[dir]},{station[height]},{new_data[time]}\n f.write(row) # 触发重训脚本只重训GWR克里金可每日全量 import subprocess subprocess.run([python, gwr_runner.py, --train-time, str(new_data[time])])关键设计GWR重训只用最新一期time24和前一期time23数据因空间关系短期稳定克里金变异函数每月全量重拟合一次即可。6.3 模型监控看板用Prometheus采集GWR ENP、克里金RMSE异常自动告警在gwr_runner.py末尾添加指标上报# 导入prometheus_client from prometheus_client import Counter, Gauge # 定义指标 gwr_enp_gauge Gauge(gwr_effective_parameters, GWR Effective Number of Parameters) gwr_rmse_gauge Gauge(gwr_prediction_rmse, GWR RMSE on test set) kriging_rmse_gauge Gauge(kriging_prediction_rmse, Kriging RMSE on test set) # 训练完成后上报 gwr_enp_gauge.set(gwr_results.ENP) gwr_rmse_gauge.set(rmse_gwr) # 计算逻辑见3.3节 kriging_rmse_gauge.set(rmse_kriging) # 计算逻辑见4.3节)启动Prometheus exporterprometheus_client.start_http_server(8000)再配Grafana看板当gwr_enp_gauge突降至5说明空间异质性消失可能污染源关停触发人工核查。从那以后我每次部署空间模型都强制走一遍data/目录校验、requirements.txt版本锁死、gwr_runner.py开头的三重断言。北京16个点看着少但空间模型的脆弱性不在于数据量而在于坐标系、带宽、变异函数这三个黑匣子——任何一个没抠准结果就是“地图画得漂亮数字全是玄学”。希望帮到你。本文还有配套的精品资源点击获取