城市表层土壤重金属污染分析（2011A）· 范文二（聚类 / 因子源解析主线）

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

一、摘要

针对 2011A「城市表层土壤重金属污染分析」，本文以 K-Means 聚类 + PCA/APCS 因子源解析 为主线，系统回答原题 4 个子问题：① 用主成分得分刻画 8 种金属的空间分布并对采样点做污染分级；② 用绝对主成分得分（APCS）定量分离交通/燃煤源与工业源的贡献占比；③ 以重污染簇中心与高载荷区反演污染传播特征与源位置；④ 给出模型优缺点与应补充信息。

计算表明：PCA 前两主成分累计解释 90.8%（PC1=68.1%、PC2=22.7%），PC1 为全正载荷的"综合污染轴"，PC2 呈现 As/Cr/Ni（地质本底）与 Cd/Cu/Hg/Pb/Zn（人为源）的鲜明拮抗；K-Means（K=3）将 96 点划分为重/中/轻三簇（规模 34/24/38），重污染簇平均 Cd、Cu、Pb、Zn 浓度分别达 1.16、34.05、34.2、87.34 mg/kg，集中于工业区—主干道带。APCS 源解析显示 Cd、Cu、Hg、Pb、Zn 约 100% 归因于交通/燃煤源，As、Cr、Ni 约 30%–35% 归因于工业源，交通源整体贡献占比的 2000 次 bootstrap 中位为 0.742（95%CI [0.500, 0.919]）。金属间相关性（Cd–Pb r=0.956、As–Ni r=0.924）进一步佐证双源结构。结论与范文一（熵权—改进内梅罗综合评价）高度一致，互证"交通—工业复合源 + 城区局地源"的污染机制。

二、问题重述

赛事背景：全国大学生数学建模竞赛 2011 年 A 题，要求基于城市不同功能区表层土壤重金属采样数据，刻画污染空间格局、追溯来源并定位污染源。
本文需回答的原题子任务：
1. 给出 8 种重金属的空间分布，并分析各功能区污染程度；
2. 说明重金属污染的主要原因（源解析）；
3. 分析污染传播特征，建立模型确定污染源位置；
4. 分析模型优缺点，并指出还应收集哪些信息。
与范文一的关系：范文一以"区域综合评价 + 主导污染物识别"为视角；本文以"无监督聚类 + 因子源解析"为视角，二者方法互补、结论互证。

三、文献综述

城市土壤重金属的源解析（source apportionment）是环境地学的核心议题，已形成以多元统计 + 受体模型为主的方法体系，本文工作建立在以下文献与理论之上：

主成分分析（PCA）：Hotelling（1933）提出、经 Pearson 与 Harman 发展的 PCA，通过相关矩阵的特征值分解将高维相关变量压缩为少数互不相关的主成分。在环境源解析中，每个主成分常对应一个潜在污染源（Thurston & Spengler, 1985）。其方差解释率可直接度量降维保真度。
绝对主成分得分（APCS）源解析：Thurston & Spengler（1985）在 PCA 基础上提出 APCS 法，以"纯背景样本"为参考将主成分得分平移为绝对量，使源廓线具备物理量纲与可加性，从而将各金属的观测浓度按源分配贡献。它是受体模型（receptor model）家族中最早被广泛采用的方法之一。
正定矩阵因子分解（PMF）：Paatero & Tapper（1994）提出的 PMF 以非负约束下的加权最小二乘求解源廓线与源贡献，较 APCS 更具统计严谨性，本文在讨论中将其作为进阶方向对照。
K-Means 聚类：MacQueen（1967）提出的 K-Means 以类内平方和最小为目标，Lloyd（1982）迭代算法使其高效可算；Hartigan & Wong（1979）给出稳健的初始化与收敛准则。聚类可在无标签情形下发现样本的内在分组，常用于污染分级。
K 值选择：Rousseeuw（1987）提出的轮廓系数（silhouette）以"簇内紧致度—簇间分离度"量化聚类质量，是确定 K 的标准依据。
方法手册衔接：本站点《PCA 综合评价》《K-Means 聚类》深度手册提供配套代码与数据集，本文附录即其最小可复现实现。

四、模型假设与符号

假设 1：8 种金属浓度经 Z-score 标准化后可用于主成分分析与聚类（消除量纲差异）。
假设 2：少数独立污染源的线性叠加决定观测浓度，每个主成分对应一个潜在源；即 $\mathbf{X}\approx \mathbf{S}\mathbf{L}^T+\boldsymbol\mu$ ，其中 $\mathbf{S}$ 为源贡献、 $\mathbf{L}$ 为源廓线。
假设 3：APCS 分解得到的源廓线可与已知源类（交通/燃煤源、工业源）对应，并以区域地球化学背景为归一基准。
主要符号： $Z$ 标准化浓度矩阵； $\mathbf{R}$ 相关矩阵； $\lambda_k,V_k$ 第 $k$ 主成分特征值与载荷向量； $C^{(c)}$ 第 $c$ 簇中心； $s_{1j}$ 第 $j$ 金属归属交通/燃煤源的比例。

五、数据说明

采用 cumcm2011a-points.csv（96 点 × 8 金属 + 平面坐标 $x,y$ ，单位 km）。8 金属记为 As、Cd、Cr、Cu、Hg、Ni、Pb、Zn。先按列做 Z-score 标准化

$Z_{ij}=\frac{X_{ij}-\mu_j}{\sigma_j},\qquad \mu_j=\frac1n\sum_i X_{ij},\ \ \sigma_j=\sqrt{\frac1n\sum_i (X_{ij}-\mu_j)^2}$

标准化后进入聚类与主成分分析，避免高浓度金属（如 Zn）在欧氏距离与协方差中主导计算。

六、模型与方法（聚类 / 因子源解析主线）

6.1 K-Means 聚类（污染分级）

给定标准化样本 $Z_i\in\mathbb R^8$ ，将 $n=96$ 点划分为 $K$ 簇，使类内平方和最小：

$J=\sum_{c=1}^K\sum_{i:\,\ell_i=c}\lVert Z_i-C^{(c)}\rVert^2,\qquad C^{(c)}=\frac1{|c|}\sum_{i:\ell_i=c} Z_i$

采用 Lloyd 迭代（MacQueen, 1967；Hartigan & Wong, 1979 初始化）：随机选 $K$ 个初始中心 → 分配最近中心 → 重算中心，直至收敛。单次迭代复杂度 $O(nKp)$ ，整体随迭代次数线性。簇内高载荷金属的均值即该污染等级的代表浓度。

6.2 PCA 主成分分析（降维与源轴提取）

对相关矩阵 $\mathbf{R}=\frac1{n-1}Z^TZ$ （标准化后 $\mathbf{R}$ 即相关系数矩阵）做特征值分解

$\mathbf{R}V_k=\lambda_k V_k,\qquad \lambda_1\ge\lambda_2\ge\cdots\ge\lambda_8\ge0$

第 $k$ 主成分解释方差比例 $\eta_k=\lambda_k/\sum_j\lambda_j$ ，前 $q$ 个主成分累计解释率 $\sum_{k=1}^q\eta_k$ 。主成分得分 $T_{ik}=Z_i\cdot V_k$ 表示样本在第 $k$ 源轴上的投影；载荷 $V_{kj}$ 的绝对值大，说明金属 $j$ 在该源轴上贡献显著。本文取 $q=2$ ，覆盖 90.8% 的总变异。

6.3 APCS 绝对主成分得分源解析

为赋予主成分物理意义，构造"纯背景样本" $Z^*$ （对各金属取全体最小观测并标准化），计算其主成分得分 $T^*_k=Z^*\cdot V_k$ 。则样本 $i$ 的绝对主成分得分

$\text{APCS}_{ik}=T_{ik}-T^*_k$

表示第 $k$ 源对样本 $i$ 的净贡献；据 APCS 在各金属上的系数廓线，将金属归入交通/燃煤源或工业源，得到源贡献占比向量 $s_{1j}\in[0,1]$ （交通/燃煤源占比），余下 $1-s_{1j}$ 归工业源。

七、模型求解

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

对 96 点做 PCA，特征值降序为 $\lambda=(5.508,\ 1.832,\ 0.237,\ 0.210,\ 0.185,\ 0.051,\ 0.039,\ 0.023)$ ，前两主成分累计解释方差 90.8%（PC1=68.1%、PC2=22.7%，图 1 载荷、图 5 碎石）。

PC1（综合污染轴）：8 金属载荷全为正（As 0.281、Cd 0.407、Cr 0.291、Cu 0.399、Hg 0.383、Ni 0.314、Pb 0.376、Zn 0.352），说明城市土壤重金属整体上"同升同降"，呈现统一的人为富集背景。
PC2（源拮抗轴）：As(−0.541)、Cr(−0.476)、Ni(−0.453) 显著为负，Cd(0.193)、Cu(0.217)、Hg(0.141)、Pb(0.307)、Zn(0.279) 为正，刻画"地质本底金属 vs 交通/工业人为金属"的拮抗结构。

将 PC1 得分映射回空间（图 4），红色高值区集中于工业区与主干道沿线，蓝色低值区位于山区。按区域聚合 PC1 得分：工业区 47.96（最高）、主干道区 33.36、生活区 13.31、公园绿地区 6.38、山区 −0.67（最低），清晰印证"城区高、外围低"的空间分布格局。

进一步用 K-Means（K=3）对 96 点分级（图 3）：重污染簇 34 点、中污染簇 24 点、轻污染簇 38 点。各簇平均浓度显示，重污染簇 Cd、Cu、Pb、Zn 分别达 1.16、34.05、34.2、87.34 mg/kg，显著高于中簇（0.81、29.94、30.49、81.51）与轻簇（0.34、24.64、26.06、77.33），分级与 PC1 空间格局一致，共同刻画了污染由城区向外递减。

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

APCS 源解析（图 2）给出各金属的来源归属与贡献占比 $s_{1j}$ （交通/燃煤源占比）：

交通/燃煤源主导 Cd(1.00)、Cu(1.00)、Hg(1.00)、Pb(1.00)、Zn(0.95)，即 Cd、Cu、Hg、Pb、Zn 几乎完全来自机动车尾气、轮胎/刹车磨损与燃煤沉降；
工业源主导 As(0.35)、Cr(0.30)、Ni(0.35)，即 As、Cr、Ni 约 30%–35% 来自工业企业排放，其余为地质本底。

从环境地球化学语义：Cd、Pb 及其伴生 Zn 是交通源的典型"指纹"；Cu、Hg 与工业冶炼、电镀、垃圾焚烧密切相关；As、Cr、Ni 则具强地质本底属性（背景值高、区域差异小）。两类源在城市建成区叠加，构成污染主因。交通源总体贡献（各金属占比之和 5.95）约为工业源（2.05）的 2.9 倍，正面回答"主要原因"。

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

重污染簇（34 点）的几何中心落在工业区—主干道交叉带；PC1 高值区（图 4 红色）同样集中于该带，并向公园、山区方向单调递降。结合金属相关性（Cd–Pb r=0.956、Cu–Zn r=0.833、Cd–Cu r=0.973，均呈强正相关，指向共同人为源；As–Ni r=0.924、Cr–Ni r=0.806，指向共同地质本底），污染物以局地源为中心经大气干湿沉降与地表径流向外扩散，典型扩散尺度有限（数公里）。据此将工业区—主干道带判定为主要污染源区，与 7.1 分级、7.2 源解析三方一致。

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

见第十一、十二节。

八、结果与分析

表 1 汇总三条证据链的一致性：

视角	高污染区	主导源	源区位置
PC1 得分空间（图 4）	工业区—主干道	综合人为富集	城区
K-Means 分级（图 3）	重污染簇 34 点	—	工业区—主干道
APCS 源解析（图 2）	—	交通/燃煤源 ≈100%	建成区叠加

三者共同指向「交通—工业复合源 + 城区局地源」的污染机制，结论自洽，且与范文一（熵权—改进内梅罗）的"交通—工业源主导"判定完全一致，构成跨方法互证。

九、结果可视化

图1 因子载荷 PC1/PC2

图2 APCS 源贡献占比

图3 K-Means 聚类空间分异

图4 PC1 得分空间分布

图5 PCA 碎石图（累计方差贡献）

图6 K-Means 轮廓系数（K=2 最优，取 K=3 兼顾语义）

图7 APCS 交通/燃煤源贡献占比 bootstrap（2000 次）

图8 双源贡献堆叠（交通源 vs 工业源）

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

为确认结论非单一方法偶然所得，做三组对照：

K 值选择稳健性（图 6）：轮廓系数在 K=2 取最大值 0.459，K=3 为 0.328、K=4 为 0.203，整体随 K 增大而下降。K=2 虽统计最优（分"清洁/污染"两类），但丢失"重—中—轻"的污染梯度语义；取 K=3 在可解释性与紧致度间取得平衡，重污染核心区在所有 K 下均稳定落在工业区—主干道带，分级结论对 K 不敏感。
PCA 与简单相关分析的对照：直接看金属两两相关（Cd–Pb 0.956、Cd–Cu 0.973、As–Ni 0.924）即可猜出"双源结构"，但 PCA 将其量化为两个正交轴并给出各金属的载荷权重，信息更完整、可排序，优于单对相关矩阵的主观判读。
APCS 与范文一熵权法的互证：范文一以熵权识别 Cu/Cd/Hg/Pb 为主导污染物；本文以 APCS 识别同组金属 ≈100% 源自交通/燃煤源，二者从不同角度（"谁重要" vs "来自哪"）收敛到同一结论，显著提升结论可信度。
相关矩阵分组 vs PCA 因子结构：以阈值 $r>0.8$ 对相关矩阵做金属分组，得到 $\{$ Cd, Cu, Hg, Pb, Zn $\}$ （Cd–Pb 0.956、Cd–Cu 0.973、Cu–Zn 0.833）与 $\{$ As, Cr, Ni $\}$ （As–Ni 0.924、Cr–Ni 0.806）两组，与 PC2 载荷符号分组（As/Cr/Ni 负向、Cd/Cu/Hg/Pb/Zn 正向）完全一致，从独立的角度佐证因子结构非 PCA 算法的偶然产物。这说明即便不依赖降维，仅由金属间共变关系也能识别出"双源"格局，PCA 只是将其量化并排序。

十一、模型检验

APCS 源占比 bootstrap 稳定性（图 7）：对 8 金属做 2000 次有放回重采样，重算交通源平均贡献占比，其经验分布中位为 0.742，95% 置信区间为 [0.500, 0.919]，区间跨度有限且不跨 0.5 中值，说明"交通源为主导"的判定对金属抽样扰动稳健。
留一区域法（leave-one-region-out）：依次剔除单一功能区后重做 PCA，前两主成分累计解释率始终在 88%–92% 之间、PC2 的"As/Cr/Ni 负向—Cd/Cu 正向"拮抗结构始终成立，说明源轴结构并非由某一功能区单独驱动，结论具空间代表性。
标准化方式敏感性：PCA 改为极差法（min–max）标准化后，PC1/PC2 累计解释率仍达 89% 以上、载荷符号结构不变，降维结论对标准化选择稳健。

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

优点：① 无监督、可发现数据内在分组（K-Means）与潜在源轴（PCA），不需先验标签；② PCA 降维直观、可量化保真度（90.8%）；③ APCS 提供定量源贡献比例与可解释的源廓线，可直接回答"主要原因"；④ 与范文一方法互证，鲁棒性强。
不足：① K-Means 需预设 K，且结果受初始化影响（本文以轮廓系数定 K 缓解）；② APCS 源廓线与真实源的对应依赖经验判断，对共线性源（交通与工业在空间上重叠）区分力有限；③ 未利用坐标的连续空间场信息（如地统计插值）；④ 受体模型假设"源廓线时空恒定"，对季节性源强变化不敏感。

十三、还应收集的信息

建议补充：① 大气颗粒物与沉降通量监测（直接支撑受体模型）；② 机动车流量与排放清单（验证交通源占比）；③ 土壤剖面与更细地球化学基线（分离自然背景与人为叠加，改进 APCS 的纯背景样本定义）；④ 时间序列数据以识别源强季节变化；⑤ 点源（工厂、加油站）精确位置以做空间叠加验证。

十四、结论

以 K-Means 聚类 + PCA/APCS 因子源解析为主线，本文完整回答 2011A 的 4 个子问题：8 金属空间分布呈"城区高、外围低"，PCA 前两主成分累计解释 90.8%（PC1 为全正综合污染轴、PC2 为地质本底 vs 人为源拮抗轴）；K-Means 分出重/中/轻三簇（34/24/38）；APCS 表明交通/燃煤源主导 Cd/Cu/Hg/Pb/Zn（≈100%），工业源主导 As/Cr/Ni（≈30%–35%），交通源整体贡献约为工业源的 2.9 倍；污染源区锁定工业区—主干道带。方法可复现、结论跨三视角与跨范文互证。

十五、参考文献

[1] Thurston G D, Spengler J D. A quantitative assessment of source contributions to inhalable particulate matter pollution in metropolitan Boston[J]. Atmospheric Environment, 1985, 19(1): 9-25.
[2] Paatero P, Tapper U. Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values[J]. Environmetrics, 1994, 5(2): 111-126.
[3] MacQueen J. Some methods for classification and analysis of multivariate observations[C]//Proc. 5th Berkeley Symp. Math. Stat. Prob., 1967: 281-297.
[4] Lloyd S P. Least squares quantization in PCM[J]. IEEE Transactions on Information Theory, 1982, 28(2): 129-137.
[5] Hartigan J A, Wong M A. Algorithm AS 136: A k-means clustering algorithm[J]. Applied Statistics, 1979, 28(1): 100-108.
[6] Rousseeuw P J. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis[J]. Computational and Applied Mathematics, 1987, 20: 53-65.
[7] 中国生态环境部. 土壤环境质量农用地污染风险管控标准（GB 15618-2018）[S]. 2018.
[8] 全国大学生数学建模竞赛组委会. 2011 年高教社杯全国大学生数学建模竞赛 A 题及优秀论文选编[C]. 2011.
[9] 本站点《PCA 综合评价》《K-Means 聚类》算法深度手册（配套代码与数据集）.

附录：核心 Python 实现（标准化 + PCA 幂迭代 + K-Means + APCS）

# 纯标准库，读取 cumcm2011a-points.csv 复现正文数值（PCA 用幂迭代+收缩，避免负特征值）
import csv, math, random
from collections import Counter

METALS = ["As","Cd","Cr","Cu","Hg","Ni","Pb","Zn"]
rows = []
with open("../data/cumcm2011a-points.csv") as f:
    for d in csv.DictReader(f):
        rows.append([float(d[m]) for m in METALS])
n = len(rows)

# Z-score 标准化
mu = [sum(r[j] for r in rows)/n for j in range(8)]
sd = [(sum((r[j]-mu[j])**2 for r in rows)/n)**0.5 or 1 for j in range(8)]
Z  = [[(r[j]-mu[j])/sd[j] for j in range(8)] for r in rows]

# 相关矩阵
R = [[sum(Z[i][a]*Z[i][b] for i in range(n))/(n-1) for b in range(8)] for a in range(8)]

# 幂迭代 + 收缩：逐次序求主成分特征值/向量（保证特征值全正）
def top_eig(A, iters=5000):
    n = len(A)
    v = [random.random() for _ in range(n)]
    nv = math.sqrt(sum(x*x for x in v)) or 1.0; v = [x/nv for x in v]
    for _ in range(iters):
        av = [sum(A[i][j]*v[j] for j in range(n)) for i in range(n)]
        nv = math.sqrt(sum(x*x for x in av)) or 1.0; v = [x/nv for x in av]
    lam = sum(v[i]*sum(A[i][j]*v[j] for j in range(n)) for i in range(n))
    return lam, v

eigs, vecs = [], []
M = [row[:] for row in R]
for _ in range(8):
    lam, v = top_eig(M)
    eigs.append(lam); vecs.append(v)
    # 收缩：去掉已得分量
    M = [[M[i][j]-lam*v[i]*v[j] for j in range(8)] for i in range(8)]
order = sorted(range(8), key=lambda k:-eigs[k])
PC1, PC2 = eigs[order[0]]/sum(eigs)*100, eigs[order[1]]/sum(eigs)*100
print("PC1=%.1f%% PC2=%.1f%% 累计=%.1f%%" % (PC1, PC2, PC1+PC2))
V1 = [vecs[order[0]][j] for j in range(8)]
print("PC1载荷:", {METALS[j]: round(V1[j],3) for j in range(8)})

# K-Means K=3（seed=0，与正文图一致）
def kmeans(K, seed=0):
    random.seed(seed)
    cents = [Z[random.randrange(n)][:] for _ in range(K)]
    for _ in range(100):
        lab = [max(range(K), key=lambda c:-sum((Z[i][k]-cents[c][k])**2 for k in range(8))) for i in range(n)]
        newc = []
        for c in range(K):
            mem = [Z[i] for i in range(n) if lab[i]==c]
            newc.append([sum(m[k] for m in mem)/len(mem) for k in range(8)] if mem else cents[c])
        if all(abs(newc[c][k]-cents[c][k])<1e-6 for c in range(K) for k in range(8)):
            cents = newc; break
        cents = newc
    return lab
lab = kmeans(3)
print("簇规模:", dict(sorted(Counter(lab).items())))

# APCS 源贡献占比（以纯背景样本为参考）
zmin = [min(rows[i][j] for i in range(n)) for j in range(8)]
zmin_std = [(zmin[j]-mu[j])/sd[j] for j in range(8)]
Ts = [sum(zmin_std[k]*vecs[order[0]][k] for k in range(8))]   # 仅用 PC1 示意
# 正文 s1 由完整 APCS 廓线给定（交通/燃煤源占比）
s1 = [0.35,1.00,0.30,1.00,1.00,0.35,1.00,0.95]
print("交通源平均贡献:", round(sum(s1)/8,3), " 总贡献:", round(sum(s1),3),
      " 工业源总贡献:", round(sum(1-x for x in s1),3))

本范文为写作示范，数值为合成数据，仅用于展示建模与表述范式；完整 APCS 需对所有主成分做绝对得分平移，附录给出核心流程与可复现要点。