企业官网建设流程全解析

简介：围绕地理加权回归与克里金插值在空气质量指数预测中的应用，这份大学课程论文以完整的 docx 文档呈现，面向地理信息系统、空间分析与环境统计方向的学习者和研究者，帮助理解空间异质性建模与空间插值的实现路径。全文从研究背景出发，依次梳理 GWR 模型的一般形式、bi-square 核函数与自适应带宽的 CV、AICc 优化过程，并结合案例展示回归结果、回归系数显著性的空间差异以及预测环节；随后转入克里金插值，讨论变异函数与最优参数选择下的污染分布制图。文末附有三段 Python 源代码，分别对应地理加权回归预测、克里格插值与地图可视化，便于读者对照论文复现实验、理解参数设置与结果解读。资源包仅含 1 个 docx 文件，体积约 1.08MB，结构紧凑、目录层级清晰。目前已有 1260 人学习，适合作为空间分析课程设计或相关论文写作的参考样本。

1. GWR 不是插值器：把空气质量指数预测拆成「趋势面 + 残差面」才走得通

第一次把 GWR 模型和克里金法放进同一条流程做空气质量指数预测的人，多半会卡在同一个岔口：这两个方法都能输出一张栅格图，到底交哪个。我见过两种典型做法，一种是只跑 GWR，出一张系数图就宣布完工；另一种更干脆，直接拿普通克里金对站点 AQI 做空间插值，把回归那部分省掉，结果谁也说服不了谁。问题出在对两层逻辑的混淆：GWR 解释的是污染物浓度随位置变化的趋势，克里金补的是趋势之外那部分还带空间相关性的残余。把两者串成回归克里金，这套方法在 AQI 预测里才站得住脚。这篇写给手里已有站点监测数据和边界矢量、准备出一张可用 AQI 分布图的人，从数据表一直写到验证曲线。

2. 建模前的数据地基：把站点监测值拼成 GWR 能吃的空间表

GWR 对数据的要求比普通回归苛刻得多：它要坐标、要空间权重、要局部样本量，任何一步没对齐，后面带宽搜出来的结果都是垃圾。我一般把这一章当作整个项目的分水岭，数据表搭错了，后面再调参数也救不回来。

2.1 因变量选 AQI 还是原始浓度，这是第一个需要确定的点

很多方案直接拿 AQI 当因变量，理由是预测目标就是它。这在空间回归里是个陷阱。AQI 是六项污染物分指数取最大值再分段映射出来的分级指标，本身带阶梯和截断，PM2.5 浓度从 35 涨到 36，AQI 可能跳 5 个单位，也可能纹丝不动。GWR 假设局部线性关系，把这种非线性映射塞进去，残差会呈现明显的系统性结构，而你误以为是空间效应没解释完，于是去调带宽，越调越偏。

我通常的做法是把 PM2.5 或 PM10 的原始浓度做因变量，必要时取对数压右偏，最后再把预测浓度按标准分段换算成 AQI。这样回归端是连续的，换算端是确定性的，两件事分清楚。如果研究区里臭氧在夏季主导，因变量要跟着换，不要全年混着 PM2.5 建一个模型——空气质量指数的首要污染物在季节上切换，混建等于让模型同时拟合两套机制。

辅助变量的挑选遵循「越贴近排放和扩散机制越好」的原则，常见组合是三类：气象类包括风速、相对湿度、温度、边界层高度，负责解释扩散条件；地表类包括 NDVI、不透水面比例、土地覆盖类型，负责解释下垫面；人为活动类包括人口密度、路网密度、工业源距离，负责解释排放强度。变量不是越多越好，超过八个还全都显著，通常说明共线性已经在作乱。

2.2 投影、坐标系和样本量，三个必须提前把关的地方

带宽的单位是距离，一旦数据还挂在经纬度上，带宽就变成了「度」，而一度纬度约 111 公里、一度经度随纬度收缩，核函数算出来的权重完全没有物理意义。这一步翻车的人特别多，因为代码不报错，只是结果莫名其妙。所以站点表读完就先投影到等距投影，坐标单位统一到米。

样本量是第二个坎。GWR 每个位置都要拟一个局部回归，样本少于三十个点，局部方程直接欠定；超过两百个点，带宽搜索的计算量又会明显上升。实操里五十到一百五十个站点是比较舒服的区间，站点太密可以按网格抽稀，太疏就得换思路，改用分区建模或者降级到全局回归加克里金。

还有一点容易忽略：辅助变量的空间分辨率要匹配。用 30 米 NDVI 和 10 公里再分析气象场拼在一张表上，等于告诉模型「近处的地表差异很重要，远处的气象差异不重要」，这个先验是人为塞进去的。常见做法是把所有栅格重采样到同一格网，或者统一到与带宽同量级的尺度。

2.3 用代码把站点表和栅格辅助变量拼成一张表

下面这段是整条流程的入口，作用是读站点、投投影、逐点取栅格值、清掉无效样本。

import geopandas as gpd import pandas as pd import numpy as np import rasterio # 1. 读站点表，至少包含 lon / lat / 监测时间 / 污染物浓度 st = pd.read_csv("station_hourly.csv", parse_dates=["time"]) st = st.dropna(subset=["lon", "lat", "pm25"]) # 时间做聚合，避免小时数据把样本量虚高成几万条 st = st.groupby(["station_id", "lon", "lat"], as_index=False)["pm25"].mean() # 2. 转点要素并投影到等距坐标系，带宽单位才能是米 pts = gpd.GeoDataFrame( st, geometry=gpd.points_from_xy(st["lon"], st["lat"]), crs="EPSG:4326" ).to_crs(gpd.GeoSeries(pts_geom).estimate_utm_crs()) # 按研究区自动估算投影带 # 3. 逐点从栅格采样辅助变量 covars = ["dem", "ndvi", "pop_den", "road_den", "blh", "wind", "rh"] for name in covars: with rasterio.open(f"raster/{name}.tif") as src: pts_ll = pts.to_crs(src.crs) vals = np.array([v[0] for v in src.sample( [(g.x, g.y) for g in pts_ll.geometry])], dtype="float32") pts[name] = vals # 4. 剔除 NoData（栅格边缘、水域、缺测站常见） pts = pts.replace(-9999.0, np.nan).replace(9999.0, np.nan) pts = pts.dropna(subset=covars + ["pm25"]) print(f"可用站点数: {len(pts)}, 投影: {pts.crs}")

代码的逻辑顺序不能乱：先聚合时间维度再做空间采样，否则一小时一条记录会把同一个站点重复算几百次，局部回归的权重全被自己吃掉。estimate_utm_crs()是按数据范围自动推投影带，比手写 EPSG 编码稳妥，如果研究区跨带，还是要手动指定。src.sample返回的是栅格原始坐标下的值，所以采样前必须把点转成src.crs，不然取到的是隔壁城市的数据而你不会收到任何报错。第四个参数dtype不写，采样结果可能是 object 数组，后面标准化时直接崩。

2.4 共线性诊断：VIF 超过 10 的变量必须处理

GWR 的局部样本本来就少，变量之间再一打架，局部系数会剧烈震荡，出现「同一个变量在相邻两个站点符号相反」的离谱结果。

from statsmodels.stats.outliers_influence import variance_inflation_factor X = pts[["dem", "ndvi", "pop_den", "road_den", "blh", "wind", "rh"]].copy() X = (X - X.mean()) / X.std() # 标准化，VIF 对量纲敏感 vif = pd.Series( [variance_inflation_factor(X.values, i) for i in range(X.shape[1])], index=X.columns, ).sort_values(ascending=False) print(vif.round(2))

标准化这一步不是可选项。高程动辄上千，湿度只有几十，不标准化的 VIF 反映的是量纲差异而不是相关性。诊断出 VIF 高的变量之后，处理顺序建议是：先看能不能合并（比如高程和地形起伏度），再考虑删掉机制上更外围的那个，最后才用主成分降维。主成分虽然能压共线性，但解释性会大幅下降，系数图就没法看了，做空气质量指数这类要出政策解读的题目，我一般尽量保留原始变量。

3. GWR 的带宽怎么定：核函数、AICc 搜索和局部共线性

带宽是 GWR 里唯一的超参数，也是唯一一个能让结果从「看起来很合理」变成「完全没法用」的旋钮。它的含义是：在拟合某个位置的回归时，多远范围内的站点还算「邻居」。带宽太宽退化成全局回归，太窄则每个位置只拟合自己，残差接近零但预测毫无泛化能力。

3.1 固定核和自适应核，选错会直接反映在城郊交界

固定核用统一的空间距离划窗，适合站点在空间上分布均匀的情形。自适应核用「最近 k 个站点」定义邻域，k 就是带宽，适合站点疏密不均——城市核心区站点密集，自适应核自动收缩；郊区站点稀疏，邻域自动放大。空气质量监测网几乎都存在城乡分布不均，所以自适应核是默认选择，只有在站点经过严格网格化布设时才考虑固定核。

这个选择的实际影响在城郊过渡带最明显。用固定核时，郊区站点周边样本太少，局部方程勉强拟合，系数图会在城乡交界处出现一块突兀的斑块，看起来像发现了一个污染热点，其实是样本不足导致的估计噪声。

3.2 带宽搜索：AICc 最小不等于最稳

mgwr 提供的搜索准则是 AICc、AIC、BIC 和 CV，其中 AICc 在中小样本下最常用，它对有效参数数量做了惩罚，能抑制带宽过小。但 AICc 曲线经常是平坦的，相邻几个带宽的 AICc 差值不到 1，这时候选哪个对系数图的影响可能很大。我的做法是除了看最小值，还把整条 AICc 曲线画出来，如果最优带宽位于曲线的陡降段边缘，就往后挪几格取一个稍大的值，用一点偏差换稳定性。

3.3 用 mgwr 跑一次带自适应核的 GWR

import numpy as np from mgwr.gwr import GWR from mgwr.sel_bw import Sel_BW coords = np.column_stack([pts.geometry.x, pts.geometry.y]) y = pts[["pm25"]].values X = pts[["dem", "ndvi", "pop_den", "road_den", "blh", "wind", "rh"]].values X = (X - X.mean(0)) / X.std(0) # 自变量标准化，系数可比 # 自适应核：bandwidth 是邻居个数，fixed=False selector = Sel_BW(coords, y, X, kernel="bisquare", fixed=False) bw = selector.search(criterion="AICc", search_method="interval", interval=10) print("最优带宽（邻居数）:", bw) model = GWR(coords, y, X, bw=bw, kernel="bisquare", fixed=False).fit() print(f"R2={model.R2:.3f} adjR2={model.adj_R2:.3f} " f"AICc={model.aicc:.2f} ENP={model.enp:.2f}") resid = model.resid_response.reshape(-1)

kernel="bisquare"是高斯核之外最常用的选择，它在带宽之外权重直接截断为零，计算快，对异常站点的抗性也略好；如果你的站点分布很不均匀又想保留远端影响，可以换成gaussian，但要接受带宽搜索变慢。search_method="interval"配interval=10表示先粗扫再细扫，站点少的时候可以把 interval 调小到 5。enp（有效参数数）这个指标要盯一眼，它接近样本量说明模型已经退化成逐点插值，这时候 R2 再高也不能用。

3.4 系数图能说明什么，不能说明什么

GWR 输出的每个站点的局部系数，常被解读成「某因素在该位置的影响强度」。这种解读有边界：局部系数是带方差的估计，两个相邻站点的系数差异如果落在置信区间内，就不该说成「影响在此处发生了转变」。mgwr 的结果对象里带filter_tvals，可以直接筛出显著区域，不显著的地方在出图时留白比涂色更诚实。

另外，系数正负翻转是最需要警惕的信号。风速对 PM2.5 的系数在绝大多数位置应该是负的，如果某片区域翻正，先别急着写「局地特殊机制」，八成是那块区域样本太少或者共线性在局部爆发。

4. 残差克里金：半变异函数怎么拟合，各向异性怎么处理

到这一步，GWR 已经给出了趋势面，剩下的残差里还有一部分是空间自相关的——相邻站点的模型误差往往同向，这部分信息不该浪费。克里金就是干这个的：对残差建模、插值，再把残差面叠回趋势面。

4.1 为什么插残差而不是直接插 AQI

直接对 AQI 做普通克里金，等于假设「浓度只跟距离有关，跟地形、气象、排放都无关」，回归项提供的信息全部被丢弃。回归克里金把空间变异拆成两部分：趋势部分由辅助变量解释，残差部分由半变异函数解释。这样做的好处不只是精度，还包括可诊断——如果残差的半变异函数几乎是纯块金效应（变程趋近于零），说明辅助变量已经把空间结构吃干净了，剩下的就是随机噪声，这时候继续插残差意义不大；如果变程很长，说明还有大尺度的空间过程没被变量捕捉到，应该回去补变量。

4.2 块金、基台、变程：三个参数各自的含义

实验半变异函数是把所有站点两两配对，按距离分组算半方差，再拟合一条理论曲线。球状模型最常用，三个参数里块金值代表微观尺度变异加测量误差，基台是总变异水平，变程是空间自相关消失的距离。三个参数的关系比数值本身更重要：块金除以基台得到的块金效应比超过 0.75，说明空间相关性很弱，插值结果基本不可信；变程要跟研究区尺度对比，变程只有几公里而研究区几十公里，残差插值只会在极局部起作用，出一张全区域图就是在放大噪声。

拟合的时候还有个容易被跳过的步骤：各向异性检查。污染物扩散受主导风向影响，残差的空间相关性常常是椭圆而不是圆。如果直接按各向同性拟合，变程会被拉成一个折中值，两个方向都不准。处理方式要么在拟合前把坐标旋转到主风向方向再压缩一个轴，要么直接用带各向异性支持的模型。

4.3 用 pykrige 拟合残差并插值

import numpy as np from pykrige.ok import OrdinaryKriging # resid 来自上一章的 GWR 结果 ok = OrdinaryKriging( coords[:, 0], coords[:, 1], resid, variogram_model="spherical", nlags=12, weight=True, # 按 lag 内点对数加权拟合，远距离 lag 更稳 anisotropy_scaling=1.6, # y 方向压缩比，来自主风向分析 anisotropy_angle=30, # 主方向角度，单位度 enable_plotting=False, verbose=False, ) print("模型参数 [psill, range, nugget]:", ok.variogram_model_parameters) # 在研究区格网上插值残差 xmin, xmax = pts.geometry.x.min(), pts.geometry.x.max() ymin, ymax = pts.geometry.y.min(), pts.geometry.y.max() gridx = np.arange(xmin, xmax, 500.0) # 500 米格网 gridy = np.arange(ymin, ymax, 500.0) z_res, ss = ok.execute("grid", gridx, gridy) # ss 是克里金方差

nlags控制距离分组的数量，站点少于五十个时设 8 到 10 就够，设太大会出现只有一两个点对的空 lag。weight=True用每个 lag 的点对数做权重，能让拟合不被远处稀疏的点对带偏。anisotropy_scaling和anisotropy_angle这两个参数不要凭感觉填，先用实验半变异函数在几个方向上分别看变程，再取比值。execute返回的第二个数组是克里金方差，这张图本身就是个诊断工具——方差大的区域就是外推区，出图时应该标注或者打码，而不是和插值结果一起当成事实发布。

4.4 趋势面叠残差面：最后一步的坐标必须一致

GWR 的predict可以在网格点上给出趋势预测，但要先把网格点构造成与建模时同样的标准化自变量矩阵。这里最常见的翻车是标准化用的均值和标准差取了训练集，预测集却重新算了一遍，两张图的量纲错位，叠加结果整体偏移几百微克。

# 预测网格也必须用训练集的 mean/std 标准化 mean_, std_ = X.mean(0), X.std(0) # X 是上一章标准化前的原始矩阵 grid_X = (grid_covars - mean_) / std_ trend = model.predict(grid_coords, grid_X)[0].reshape(len(gridy), len(gridx)) aqi_surface = trend + z_res # 趋势 + 残差 if log_transform: # 若因变量取过对数，要还原 aqi_surface = np.exp(aqi_surface)

model.predict的返回值是一个元组，第一项才是预测值，直接拿去 reshape 会报维度错。如果因变量做了对数变换，记得在叠加之后统一还原，先还原再叠加会让残差失去可加性。克里金方差不要丢，它可以和残差面一起构成不确定性图层，是这套流程比单一插值更有说服力的地方。

5. 排查清单：GWR-Kriging 流程里最容易失效的五件事

5.1 系数图一片斑驳，R2 却很高

现象是模型报告的 R2 到了 0.85 以上，但局部系数图看起来像噪声，相邻站点符号来回跳。原因通常是带宽被 AICc 搜到了最小值，模型退化成逐点插值，有效参数数接近样本量。解决办法是强制给带宽设下界，一般不小于三十个邻居，同时把enp打印出来和样本量对比，超过样本量的 70% 就说明模型已经没在「回归」了。这个坑我在早期项目里踩过，当时还以为是发现了局地机制，白写了两周分析。

5.2 残差半变异函数几乎是纯块金

现象是拟合出来的变程只有一两公里，块金效应比超过 0.8，插出来的残差面全是小斑点。原因要么是关键辅助变量缺失，空间结构还留在残差里但表现为短程噪声；要么是站点坐标存在系统性偏差，比如站点搬迁后坐标没更新。解决办法是先补变量（把排放源、道路密度这类强空间变量加进来），再核对站点坐标与影像是否对得上，最后才考虑是不是研究区本身空间相关性就弱。

5.3 插值结果在边界外大面积外推

现象是研究区边缘外的栅格值出现极端高低值，或者克里金方差大得离谱。原因是克里金在站点凸包之外完全靠趋势外推，没有数据约束。解决办法是把插值范围限制在站点凸包加一个缓冲区，超出范围直接掩膜为 NoData；如果确实要覆盖行政区全境，就明确在图上标注外推区，并且用克里金方差做透明度。这是最容易被忽略的一步，出的图看着完整，其实一半是猜的。

5.4 换一年数据，带宽和变程全变了

现象是同一个研究区，去年最优带宽 60，今年变成 120，变程也翻倍。原因可能是站点增减、也可能是模型对样本扰动过于敏感，这在样本量不足一百时很常见。解决办法是把带宽固定在两年数据上的折中值，或者干脆用时间分块交叉验证来选带宽，而不是每年单独搜。参数年际跳变本身就是一个信号，说明模型稳健性不够，结论不该下得太硬。

5.5 交叉验证分数漂亮，实际图对不上

现象是十折交叉验证 R2 有 0.8，但把预测图和独立年份的实际监测值一比，误差集中在城郊。原因是随机划分折数时，相邻站点被分到了训练集和验证集两边，空间自相关让验证集「偷看」了答案。解决办法是换成空间分块交叉验证，按地理网格或聚类划块，保证验证块与训练块在空间上分离。这一步不做，前面所有指标都不可信。

6. 验证与进阶：留一站交叉验证、空间分块和可复用管线

模型跑通只是起点，能不能拿来出图、能不能支撑结论，全看验证怎么做。

6.1 三种验证方式的适用边界

留一站交叉验证在站点不多时最实在，因为每次只留一个点，训练集几乎完整，误差估计偏差小。它的代价是要拟合 n 次 GWR，n 是一百时大概几分钟到十几分钟，可以接受。空间分块更适合回答「这套模型能不能推到没站点的区域」这个问题，分块尺度建议从带宽的 1.5 倍开始试。

from sklearn.model_selection import KFold import numpy as np def loo_gwr_krig(coords, y, X, bw, kernel="bisquare"): n = len(y) pred = np.zeros(n) for i in range(n): tr = np.ones(n, dtype=bool); tr[i] = False m = GWR(coords[tr], y[tr], X[tr], bw=bw, kernel=kernel, fixed=False).fit() # 用训练好的局部模型预测留出的点，再叠加该点的残差克里金估计 trend_i = m.predict(coords[i:i+1], X[i:i+1])[0][0, 0] pred[i] = trend_i + m.resid_response.mean() # 简化处理，残差取均值 rmse = np.sqrt(np.mean((pred - y.ravel()) ** 2)) return pred, rmse pred, rmse = loo_gwr_krig(coords, y, X, bw) print(f"LOO RMSE = {rmse:.2f} μg/m³")

这段是简化版，残差部分用训练集残差均值代替完整的克里金预测，实际项目里应该在每次留出时对训练集残差重新拟合半变异函数再插值，代价是时间翻几倍。如果只是为了比较不同带宽的稳定性，简化版够用；如果要写进报告，就用完整版。

6.2 把流程封成可复用管线

真正耗时间的不是第一次跑通，而是参数微调后要重跑十几次。我一般把「数据准备 → 带宽搜索 → GWR 拟合 → 残差克里金 → 网格叠加 → 验证」封成一个函数，输入是站点表和栅格目录，输出是 GeoTIFF 加一份指标 JSON。封装时有两个细节值得注意：一是把标准化用的均值和标准差一起存下来，预测新数据时直接复用，不要重算；二是把克里金方差一并输出，作为不确定性图层。

6.3 带宽和变程的耦合，是最需要经验的地方

带宽变大，GWR 趋势面更平滑，残差里保留的空间结构更多，克里金的变程通常会跟着变长；带宽变小，趋势面贴合数据，残差接近白噪声，变程缩短。这两个参数不是独立的，调一个另一个就得跟着看。我的习惯是固定克里金模型类型，只扫带宽，画出「带宽 - LOO RMSE」曲线，取曲线平坦段的中点，然后再在这个带宽下重新拟合一次半变异函数。反过来先定变程再调带宽，很容易掉进局部最优。

说个实在的教训：早些年我总觉得参数搜到最优就完事了，后来发现每一次参数年际跳变、每一张边界外推图，都是模型在提醒我它对数据边界的依赖比我想的更强。现在我做这类空间预测，出图之前一定先看三张诊断图——残差半变异函数、克里金方差、局部系数显著性，三张里有一张不对劲就不往下走。这套流程不难，难的是忍住不跳过验证那一步。希望帮到你。

本文还有配套的精品资源，点击获取