城市表层土壤重金属污染分析（2011A）· 范文三（空间插值 / 扩散反演主线）

本范文为写作范式示范：由本站基于 2011A 赛题要点撰写，非真实参赛论文；核心数值取自本站合成练习集 cumcm2011a-points.csv（96 采样点，含坐标与 8 种重金属浓度），用于展示「地统计插值 + 扩散反演 + 源定位」类题型的建模与表述结构，正式参赛请以官方数据与真实方法为准。

一、摘要

针对 2011A「城市表层土壤重金属污染分析」，本文以 地统计空间插值 + 指数扩散反演 + 源定位 为主线，系统回答原题 4 个子问题：① 以变异函数与 IDW 插值刻画 Cd 的空间连续分布并评估各功能区污染程度；② 由扩散衰减模型与源反演锁定污染源位置、判定污染主因；③ 建立"局部源 + 指数衰减"的污染传播模型并量化扩散尺度；④ 给出模型优缺点与应补充信息。

计算表明：Cd 的实验变异函数以指数模型拟合最优（SSE=0.002），块金/基台比仅 0.0014，空间自相关极强，变程约 7.5 km；IDW 留一法交叉验证最优邻点数 K=2（RMSE=0.112）。污染源反演锁定主源（工业区）位于 (7.08, 2.54) km、Cd=1.701 mg/kg，次源（交通干线）位于 (7.72, 2.68) km、距主源仅 0.65 km；扩散衰减模型 $\text{Cd}(d)=1.701\,e^{-0.25d}+0.280$ 的残差标准差 0.279、最大绝对残差 0.70。源定位 bootstrap（2000 次）中位偏移 0.292 km、97.5% 分位 1.301 km；衰减系数 λ 的 95% CI 为 [0.220, 0.270]。结论与范文一（熵权—改进内梅罗）、范文二（PCA/APCS 源解析）三方互证"工业区—主干道复合源"的污染机制。

二、问题重述

赛事背景：全国大学生数学建模竞赛 2011 年 A 题。
本文视角：前两篇分别从"区域综合评价"与"因子源解析"切入；本文以空间连续场与扩散物理过程为视角，把"污染在哪、怎么传、源在哪"落到地理坐标与可积扩散模型上。
需回答的原题子任务：① 空间分布与各功能区污染程度；② 污染主要原因；③ 传播特征与污染源定位；④ 模型优缺点与应补充信息。

三、文献综述

空间统计与大气扩散是城市污染研究的另一支柱，本文建立于以下理论与文献：

变异函数（variogram）：Matheron（1963）提出半变异函数 $\gamma(h)=\frac12\mathrm{Var}[Z(x)-Z(x+h)]$ ，用以刻画空间自相关随距离的变化。其三个核心参数——块金（nugget，微观变异/测量误差）、基台（sill，总变异）、变程（range，自相关消失距离）——是地统计建模的基础。块金/基台比越小，空间结构性越强。
克里金（Kriging）与 IDW：Journel & Huijbregts（1978）系统化的克里金以变异函数为权重做最优无偏插值；Shepard（1968）的反距离加权（IDW）则以距离幂次反比赋权，实现简单、对局部极值敏感，是本文插值对比的基线方法。
高斯烟羽扩散：Pasquill–Gifford 框架下的指数/高斯衰减描述污染物由点源向外浓度衰减，其核心是衰减系数与扩散尺度；将其用于土壤累积浓度的空间格局反演，可在"源强—距离"间建立可积模型。
源反演：基于"观测浓度 = 源强 × 衰减核"的反问题，以最速下降或网格搜索估计源位置与源强（类比大气源反演，Okubo, 1980）。
方法手册衔接：本站点《地统计与变异函数》《空间插值 IDW/Kriging》手册提供配套代码与数据集。

四、模型假设与符号

假设 1：浓度场在空间上连续，可用变异函数描述其自相关结构。
假设 2：存在有限个局地源，观测浓度近似为各源指数衰减场的叠加（忽略平流各向异性时取径向对称）。
假设 3：以 Cd 为代表污染物（交通/工业源指纹明显、变异大），其空间格局可代表整体污染的空间形态。
主要符号： $\gamma(h)$ 半变异函数； $c_0$ 块金、 $c$ 基台、 $a$ 变程； $\hat Z(x_0)$ IDW 预测值； $\lambda$ 扩散衰减系数； $s$ 源坐标。

五、数据说明

采用 cumcm2011a-points.csv（96 点，含 $x\_km,y\_km$ 坐标及 8 金属）。以 Cd 为代表污染物（权重高、交通源特征明显）开展空间分析，浓度单位 mg/kg。

六、模型与方法（空间插值 / 扩散反演主线）

6.1 变异函数与空间结构

计算实验半变异函数：对所有点两两滞后距 $h$ ，取 $\gamma(h)=\frac12(h\text{ 仓内})(Z_i-Z_j)^2$ 的均值。以三类理论模型拟合：

球状（spherical）： $\gamma(h)=c_0+c\,(1.5h/a-0.5(h/a)^3),\ h\le a$ ；
指数（exponential）： $\gamma(h)=c_0+c\,(1-e^{-h/a})$ ；
高斯（gaussian）： $\gamma(h)=c_0+c\,(1-e^{-(h/a)^2})$ 。

以网格搜索最小化 SSE 定参，SSE 最小者胜出。

6.2 IDW 插值

对未知点 $x_0$ 取最近 $K$ 邻点，预测

$\hat Z(x_0)=\frac{\sum_{j=1}^K Z(x_j)/d_{0j}^p}{\sum_{j=1}^K 1/d_{0j}^p},\qquad p=2$

幂 $p=2$ 固定，邻点数 $K$ 由留一法交叉验证（LOO-CV）选优，RMSE $=\sqrt{\frac1n\sum_i(\hat Z_{-i}(x_i)-Z(x_i))^2}$ 。

6.3 扩散反演与源定位

设单主源模型 $\widehat{\text{Cd}}(d)=A\,e^{-\lambda d}+B$ ，其中 $d$ 为到源的距离， $A=Z_{\max}-Z_{\min}$ 、 $B=Z_{\min}$ 。以网格搜索 $\lambda$ 最小化拟合 SSE 锁定衰减系数；源坐标取使"观测浓度最高"与"全局拟合最优"一致的位置。进一步以 bootstrap 重采样（2000 次）取 Cd 最高点相对名义源的偏移，刻画定位不确定性。

七、模型求解

7.1 空间分布与区域污染程度（任务一）

Cd 的实验变异函数（图 1）以指数模型拟合最优（SSE=0.002，优于球状 0.006、高斯 0.005），拟合得块金 $c_0=0.0005$ 、基台 $c=0.350$ 、变程 $a=7.5\ \text{km}$ 。块金/基台比仅 $0.0005/0.350\approx 0.0014$ ，说明 Cd 空间自相关极强、微观随机变异极小，浓度场高度结构化——污染呈"成片连片"而非"随机散点"。变程 7.5 km 给出污染影响的特征空间尺度。

IDW 插值面（图 2）显示高值区位于工业区—主干道一带并向外围递降。留一法交叉验证（图 6）中，邻点数 K=2 时 RMSE 最小（0.112），K 增大后 RMSE 缓升至 0.124，说明少量近邻即可刻画 Cd 的空间细节、过多邻点会平滑掉局部峰值。由此刻画的各功能区污染程度与范文一、范文二一致（城区高、山区低）。

7.2 污染主要原因（任务二）

扩散反演锁定 主源（工业）位于 (7.08, 2.54) km，该点 Cd=1.701 mg/kg 为全域最高；次源（交通干线）位于 (7.72, 2.68) km，Cd=1.663 mg/kg，距主源仅 0.65 km（图 3）。两源紧邻且均落在工业区—主干道交叉带，与范文二 APCS 判定的"交通/燃煤源 + 工业源复合"完全对应：主源 Cd 极高指向工业企业排放，次源紧邻主干道指向机动车源。由此可见 Cd（及伴生 Zn、Pb）的高浓度主要由该复合源区产生，局部源的人为排放是污染主因。

7.3 传播特征与源定位（任务三）

扩散衰减模型（图 4）

$\widehat{\text{Cd}}(d)=1.701\,e^{-0.25d}+0.280$

以距离 $d$ 解释 Cd 浓度：残差标准差 0.279、最大绝对残差 0.70，拟合合理。衰减系数 $\lambda=0.25\ \text{km}^{-1}$ 对应的特征衰减距离 $1/\lambda=4.0\ \text{km}$ ，与变异函数变程 7.5 km 同量级，二者独立印证污染扩散的"影响半径"约 4–7 km——即污染物以源为中心、在数公里尺度内经大气沉降与地表径流向外衰减，超出该范围浓度回落至背景。源定位 bootstrap（图 7）中位偏移 0.292 km、97.5% 分位 1.301 km，说明源区定位在亚公里至 1–2 km 精度内稳定，非个别高值点驱动。

7.4 模型优缺点与拓展（任务四）

见第十一、十二节。

八、结果与分析

表 1 汇总三视角证据对"源区—扩散"的一致刻画：

方法	关键参数	源区/尺度	推论
变异函数（图 1）	块金/基台 0.0014、变程 7.5 km	连片结构化	污染成片、影响半径 ~7 km
IDW 插值（图 2/6）	K=2 RMSE 0.112	工业区—主干道高值	城区高、外围低
扩散反演（图 3/4）	主源(7.08,2.54)、λ=0.25	双源紧邻（0.65 km）	局部源指数衰减、半径 ~4 km

三者共同支撑"工业区—主干道复合源 + 数公里尺度指数衰减"的污染传播机制。

九、结果可视化

图1 Cd 实验变异函数与球状拟合

图2 Cd 浓度 IDW 插值面

图3 污染源反演定位

图4 扩散衰减与95%置信带

图5 变异函数三种模型拟合 SSE 对比

图6 IDW 留一法交叉验证 RMSE（邻点数 K）

图7 源定位 bootstrap 稳定性

图8 扩散衰减系数 λ 蒙特卡洛分布

十、方法横评（灵敏度与对比）

变异函数模型选择：球状/指数/高斯三模型 SSE 分别为 0.006/0.002/0.005，指数模型最优。三者的变程估计均在 4–7.5 km 区间，结论（空间自相关强、影响半径数公里）对模型选择稳健。
IDW 邻点数 K 的稳健性：K=2→12 的 RMSE 在 0.112–0.124 窄带内波动，最小值稳定落在 K=2，说明插值结论不依赖 K 的精细取值；若改用克里金，预计因充分利用变异函数结构而略有改善，但不会改变高值区位置。
衰减模型 vs 纯空间相关：单源指数衰减以仅一个参数 λ 即解释 Cd 的径向格局（残差 0.279），与"局部源主导"的物理直觉一致；相比把浓度直接当作坐标的任意函数，衰减模型可外推至无采样点区域，实用性更强。

十一、模型检验

源定位 bootstrap（图 7）：2000 次有放回重采样取 Cd 最高点，相对名义源的偏移中位 0.292 km、95% 区间 [0.000, 1.301] km，表明源区定位在亚公里精度内稳定，不会因个别高值点漂移而改变。
衰减系数 λ 蒙特卡洛（图 8）：2000 次重拟合得 λ 中位 0.246、95% CI [0.220, 0.270]，区间窄且包含名义值 0.25，说明扩散尺度估计稳健。
IDW 留一法（图 6）：LOO-CV RMSE 全 K 均 < 0.13，交叉验证误差小，插值面可信；残差无系统性空间聚集（最大 0.70 出现在次源附近的高梯度区），模型偏差可控。

十二、模型评价（优缺点）

优点：① 把污染落到地理坐标与可积扩散模型，结论直观、可外推；② 变异函数量化空间自相关结构（块金/基台/变程），信息丰富；③ 留一法、bootstrap、蒙特卡洛三重检验闭环，结论稳健；④ 与范文一、范文二跨方法互证。
不足：① 以 Cd 单指标代表整体，未对 8 金属分别建场（各金属扩散尺度可能不同）；② 单源衰减模型忽略多源叠加与风向各向异性（真实城市扩散非径向对称）；③ IDW 为局部插值、在采样空白区可能失真，理想应上克里金；④ 衰减系数 λ 的物理含义（沉降速率/风速）未与现实气象参数挂钩。

十三、还应收集的信息

建议补充：① 高密网格或规则栅格采样，提升插值分辨率；② 气象（风向、风速、降水）以建模各向异性扩散；③ 点源（工厂坐标、交通流量）以做多源叠加反演；④ 不同深度土壤剖面以区分表层累积与下层背景；⑤ 时间序列以估计源强年际变化。

十四、结论

以地统计插值 + 指数扩散反演 + 源定位为主线，本文完整回答 2011A 的 4 个子问题：Cd 空间自相关极强（变异函数块金/基台 0.0014、变程 7.5 km，指数模型 SSE=0.002 最优），IDW 最佳 K=2（RMSE 0.112）；污染源反演锁定主源（工业，1.701 mg/kg）与次源（交通，距主源 0.65 km）于工业区—主干道带；扩散衰减 $\text{Cd}(d)=1.701e^{-0.25d}+0.280$ 、特征半径 ~4 km；源定位 bootstrap 中位偏移 0.292 km、λ 的 95%CI [0.220, 0.270]。方法可复现、结论与另两篇范文三方互证。

十五、参考文献

[1] Matheron G. Principles of geostatistics[J]. Economic Geology, 1963, 58(8): 1246-1266.
[2] Journel A G, Huijbregts C J. Mining Geostatistics[M]. Academic Press, 1978.
[3] Shepard D. A two-dimensional interpolation function for irregularly-spaced data[C]//Proc. 23rd ACM National Conference, 1968: 517-524.
[4] Pasquill F, Smith F B. Atmospheric Diffusion[M]. 3rd ed. Wiley, 1983.
[5] Okubo A. Diffusion and Ecological Problems: Mathematical Models[M]. Springer, 1980.
[6] 中国环境监测总站. 中国土壤元素背景值[M]. 中国环境科学出版社, 1990.
[7] 全国大学生数学建模竞赛组委会. 2011 年高教社杯全国大学生数学建模竞赛 A 题及优秀论文选编[C]. 2011.
[8] 本站点《地统计与变异函数》《空间插值 IDW/Kriging》算法深度手册（配套代码与数据集）.

附录：核心 Python 实现（变异函数 + IDW + 扩散反演 + 源定位）

# 纯标准库，读取 cumcm2011a-points.csv 复现正文数值
import csv, math, random

rows = []
with open("../data/cumcm2011a-points.csv") as f:
    for d in csv.DictReader(f):
        rows.append([float(d[m]) for m in ["As","Cd","Cr","Cu","Hg","Ni","Pb","Zn"]]
                    + [float(d["x_km"]), float(d["y_km"])])
cd = [r[1] for r in rows]; xs = [r[8] for r in rows]; ys = [r[9] for r in rows]
n = len(cd)
def mean(xs): return sum(xs)/len(xs)

# 实验变异函数（14 仓，仓宽 1km）
pairs = [(math.hypot(xs[i]-xs[j], ys[i]-ys[j]),
          0.5*(cd[i]-cd[j])**2) for i in range(n) for j in range(i+1,n)]
nb, w = 14, 1.0; bins = [[] for _ in range(nb)]
for h, g in pairs:
    bins[min(nb-1, int(h/w))].append(g)
hc = [(b+0.5)*w for b in range(nb) if bins[b]]
gc = [mean(bins[b]) for b in range(nb) if bins[b]]

def fit_sse(model, c0v, cv, av):
    sse = 0.0
    for h, g in zip(hc, gc):
        if model == "sph":
            pred = c0v + (cv*(1.5*h/av - 0.5*(h/av)**3) if h <= av else cv)
        elif model == "exp":
            pred = c0v + cv*(1 - math.exp(-h/av))
        else:
            pred = c0v + cv*(1 - math.exp(-(h/av)**2))
        sse += (pred - g)**2
    return sse

# 网格搜索（指数模型示例）
best, bp = 1e9, None
for c0v in [0.0001, 0.0002, 0.0005]:
    for cv in [i/100 for i in range(5, 40)]:
        for av in [i/2 for i in range(3, 16)]:
            s = fit_sse("exp", c0v, cv, av)
            if s < best: best, bp = s, (c0v, cv, av)
print("指数模型 SSE=%.4f 参数(块金,基台,变程)=%.4f,%.4f,%.2f" % (best, bp[0], bp[1], bp[2]))

# IDW 留一法（K=2）
def idw_loo(K):
    err = []
    for i in range(n):
        dsts = sorted((math.hypot(xs[i]-xs[j], ys[i]-ys[j]), j)
                      for j in range(n) if j != i)[:K]
        num = sum(cd[j]/(dd**2+1e-6) for dd, j in dsts)
        den = sum(1/(dd**2+1e-6) for dd, j in dsts)
        err.append((num/den - cd[i])**2)
    return math.sqrt(mean(err))
print("IDW K=2 RMSE=%.3f" % idw_loo(2))

# 主源与扩散反演
src = max(range(n), key=lambda k: cd[k]); sx, sy = xs[src], ys[src]
print("主源 id=%d (%.2f,%.2f) Cd=%.3f" % (src+1, sx, sy, cd[src]))
zmax, zmin = max(cd), min(cd)
dists = [math.hypot(xs[i]-sx, ys[i]-sy) for i in range(n)]
best, lam = 1e9, 0.1
for L in [k/100 for k in range(5, 120)]:
    sse = sum((cd[i]-(zmax*math.exp(-L*dists[i])+zmin))**2 for i in range(n))
    if sse < best: best, lam = sse, L
resid = [cd[i]-(zmax*math.exp(-lam*dists[i])+zmin) for i in range(n)]
print("λ=%.2f 残差标准差=%.4f 最大|残差|=%.4f" % (lam, math.sqrt(mean([r*r for r in resid])), max(abs(r) for r in resid)))

本范文为写作示范，数值为合成数据，仅用于展示建模与表述范式。