农业活动是导致地下水硝酸盐污染的主要因素：白洋淀流域中的机器学习与基于风险信息的监测优化方法

《Agricultural Water Management》：Agricultural legacies drive groundwater nitrate contamination: Machine learning and risk-informed monitoring optimization in the Baiyangdian Lake watershed

　　张志雄|王世勤|张子龙|周景雄|刘炳霞
中国科学院遗传与发育生物学研究所农业资源研究中心，河北省农业节水重点实验室，石家庄050022

**摘要**
地下水硝酸盐污染对集约化农业地区的水资源安全构成了严重威胁。由于对时空动态的理解不足，尤其是历史耕作方式的遗留影响以及非饱和带调节作用的理解有限，有效的管理措施受到阻碍。本研究以华北平原（NCP）典型的农业区白洋淀流域为研究对象，利用包含20年作物变化数据、详细非饱和带地层信息及278个地下水样本的独特数据集，进行了全面分析。采用可解释的机器学习框架（XGBoost-SHAP）来量化历史农业活动和非饱和带特性的相对重要性，从而开发出SHAP-WHO综合风险指数（SWRI），该指数将预测结果与特征贡献相结合，用于高分辨率风险制图。模型显示，21世纪初的灌溉和施肥率是当前硝酸盐水平的主要预测因子。研究发现存在数十年滞后现象：21世纪初的灌溉和施肥率对当前硝酸盐水平具有显著影响。非饱和带特性（尤其是全新世沉积物的厚度）和长期平均地下水位深度是控制硝酸盐迁移和衰减的关键因素。SWRI具有强大的预测能力，并通过基于风险加权的多目标空间聚类算法优化了地下水监测网络，从而在保持极高风险区域全覆盖的同时，减少了70.5%的监测站点数量。这项研究为历史硝酸盐污染的机制提供了理解，并为全球可持续地下水管理提供了一个可推广的、数据驱动的框架。

**1. 引言**
地下水在全球淡水安全中发挥着不可替代的作用，满足了超过25亿人的用水需求，尤其是在地表水资源稀缺的干旱和半干旱地区（Famiglietti, 2014; Niazi et al., 2024）。在这些水资源匮乏的地区，地下水往往是灌溉的生命线，直接关系到作物产量和粮食安全。然而，这种过度依赖加上不可持续的农业管理方式（如过量施用氮肥和低效灌溉）导致了普遍且紧迫的问题：地下水硝酸盐污染严重，其浓度经常超过世界卫生组织规定的50毫克/升的安全阈值（Wu et al., 2024; Zhou et al., 2024a）。

华北平原（NCP）正是这一全球性难题的典型代表。作为全球灌溉最集约化的农业区之一，NCP约70%的用水依赖地下水（Liu et al., 2011; Yang et al., 2021）。数十年的集约化耕作导致地下水硝酸盐污染严重，带来了深远后果：威胁饮用水安全（例如，高铁血红蛋白血症和某些癌症的风险（Chaudhary et al., 2025）），并加剧了相关地表水体的富营养化，包括生态脆弱的白洋淀湖（Hao et al., 2024; Zhu et al., 2019）。尽管人们对这一问题的认识日益增强，但由于对硝酸盐从源头到含水层迁移的时空动态了解有限（Xu et al., 2025），有效管理仍受到限制。传统方法通常假设当前土地利用与地下水质量之间存在直接且同时发生的关联（例如，使用单一年份的土地覆盖图），忽视了地表氮肥施用与其进入饱和带之间的时间滞后（Basu et al., 2022; Gao et al., 2021; Zhang et al., 2023）。这种滞后时间可能从几年到几十年不等，主要受非饱和带过程控制，包括渗透率、土壤质地和非饱和带厚度（Ascott et al., 2021; Guo and Wu, 2024; He et al., 2024）。在NCP，第四纪沉积物厚度可超过50米，由淤泥、粘土和沙子交替组成，非饱和带起到复杂的过滤和延迟作用，使地表活动与地下水响应脱钩（Liu et al., 2022）。因此，由于传统研究方法忽略了氮迁移的时间滞后效应，加上非饱和带的复杂结构和水文过程，硝酸盐及其在非饱和带中的迁移特性尚未得到充分阐明。

近年来机器学习（ML）的进步为建模这些非线性、多尺度相互作用提供了新的机会。随机森林（RF）、XGBoost和人工神经网络（ANNs）因能够处理高维、非线性数据且无需强分布假设而越来越被应用于地下水质量预测（Alam et al., 2025; Deng et al., 2023; Pham et al., 2022）。例如，ML模型已用于农业流域的地下水硝酸盐风险制图，其准确性优于传统的基于指数的方法（如DRASTIC，Baalousha, 2025; Rodriguez-Galiano et al., 2014）。然而，大多数ML应用仍依赖于静态或短期输入变量，往往忽视了历史农业管理的累积影响。此外，尽管ML模型提供的变量重要性指标（如RF中的Gini杂质减少）很有用，但它们可能偏向于相关或高基数特征（Saarela and Jauhiainen, 2021）。基于合作博弈论的SHAP值通过考虑所有可能的特征组合，提供了更稳健和可解释的特征贡献度量（Aas et al., 2021; Li et al., 2024），有助于深入理解机制。

另一个关键缺口在于监测网络设计。当前的地下水监测系统主要采用基于便利性或空间均匀性的采样策略（Feng et al., 2023b; Levintal et al., 2022），导致资源分配不合理，高风险区域（如农业热点和脆弱含水层补给区）的覆盖不足（Farshad et al., 2025; Song et al., 2025; Zhou et al., 2024b）。虽然最近有研究提出了利用K均值聚类（Teimoori et al., 2023）或基于熵的空间指标（Xiong et al., 2022）的优化框架，但这些方法大多静态且缺乏动态、基于过程的风险指数整合。关键在于，它们忽略了将源特征因素（如施肥强度、土地利用模式）与路径脆弱性（如沉积物厚度、水文路径）相结合，因此产生的监测网络缺乏对污染热点的预测能力。这种方法论限制阻碍了能够优先考虑污染潜在高风险区域的自适应监测系统的发展，最终影响了环境风险评估的及时性和准确性。

为解决这些关键问题，我们对白洋淀流域的地下水硝酸盐动态进行了全面研究，利用了涵盖20年作物类型变化、详细非饱和带地层信息和密集地下水采样（n=278个样本）的综合性数据集。我们的目标有三个：（1）使用XGBoost-SHAP框架量化历史农业活动和非饱和带特性在控制地下水硝酸盐浓度方面的相对重要性；（2）开发一个基于物理原理、数据驱动的硝酸盐风险指数，以捕捉遗留效应和地下异质性；（3）通过将风险指数整合到多目标空间聚类算法中，优化长期监测网络。通过结合历史背景数据、物理过程和监测数据，这项工作可能有助于深入理解农业系统中的地下水污染机制，并为可持续水资源管理提供可推广的框架。

**2. 材料与方法**
**2.1. 研究区域**
白洋淀流域位于华北平原（NCP）西北部（38.02°–39.86°N, 114.30°–116.49°E），面积约为16,100平方公里，代表了NCP典型的山前冲积平原。该地区属于温暖温带季风气候，年平均气温为13.3°C，降水量为504毫米，其中70%发生在夏季（6月至8月）。地形从西北向东南逐渐倾斜，形成了典型的冲积-坡积扇、洪泛平原和湖泊洼地序列（图1b）。白洋淀湖是华北最大的淡水湖，位于该系统的末端，具有重要的生态缓冲和储水功能。区域地表水网络由几条受严格管理的季节性河流组成，这些河流最终流入湖泊。由于农业和家庭用水的回流，这些地表水体经常表现出较高的氮含量，加剧了区域富营养化风险。

该地区覆盖着厚度从200到600米的第四纪沉积物。非饱和带主要由全新世（Q4）、上更新世（Q3）和中更新世（Q2）沉积物组成（图1b）。这些沉积物表现出强烈的垂直和横向异质性：低洼地区的地层以细淤泥和粘土为主，而古河道和扇形沉积物中则含有较粗的沙子和砾石。本研究的主要研究对象——非限制含水层位于地表下10–120米深处，主要由灌溉回流和降水补给。地下水从西南和西北方向流向中心区域，最终流入白洋淀湖，平均水力梯度为0.09%。值得注意的是，该地下水系统受到人类活动的严重影响。地下水位存在显著的空间变化（平均变化8.14米），在1980–2015年间因过度开采而急剧下降。最近观察到恢复趋势，2022年3月的平均地下水位为19.7米（图1d）。农业耕作是主要土地利用方式，占流域面积的60%-77%（图1c）。典型的种植制度包括冬季小麦（10月至6月）和夏季玉米（6月至9月）的双季轮作，在城郊地区逐渐转向高价值蔬菜生产（如温室番茄、叶类蔬菜）。自20世纪80年代以来，氮肥施用量大幅增加，在集约化区域平均为350–500公斤/公顷/年，远超过作物吸收需求（Feng et al., 2023a; Tan et al., 2025）。

**2.2. 数据收集与处理**
为了支持本研究，我们整合了来自多个来源的综合性数据集：

- **地下水化学**：2022年3月进行了空间均匀的采样，收集了278个地下水样本，使用标准离子色谱和电位法分析了主要离子指标（包括硝酸盐NO??、氯化物Cl?和氧化还原电位ORP）。
- **历史土地利用和作物类型**：基于2018至2023年在流域内的实地调查和地下水采样，记录了采样点的作物类型和轮作模式。利用Landsat 8影像（空间分辨率30米）收集了340个参考样本，并按6:4的比例分为训练集和验证集。利用MODIS NDVI数据的23个周期年度时间序列（捕捉生长季节的物候动态），在ENVI平台上使用决策树算法将2000年、2010年和2020年的农业土地利用类型分为九类（包括七种单一种植类型、两种双季轮作和城市区域）。2020年的分类与地面真实验证点的Kappa系数为0.818，显示出高准确性。2000年和2010年的分类同样使用Landsat 7影像和相同的划分协议进行。
- **地质和水文数据**：关键变量（如全新世沉积物厚度d_qh、晚更新世沉积物厚度d_q3以及全新世-中更新世古环境分区（图1b）来自区域调查报告（Zhang et al., 2009）。1980年、1990年、2000年、2006年、2015年和2021年的地下水位数据来自62至95个监测井（https://geocloud.cgs.gov.cn/），通过空间插值计算长期平均水位（WTD_avg）。年降水量（Avg_pre_y）和灌溉数据来自当地气象站和农业统计资料。
- **历史农业活动**：变量如Fertilization2000、Irrigation2000、Fertilization2010、Irrigation2010、Fertilization2020和Irrigation2020表示Thiessen多边形内的单位面积施肥和灌溉率。变量2000_7、2000_9、2010_7、2010_9、2020_7和2020_9表示同一Thiessen多边形内的种植区域，“7”代表蔬菜种植，“9”代表冬季小麦-夏季玉米轮作。表1总结了数据集的特征。数据分析中使用的数据来源和变量总结

变量类别 变量来源/描述
地下水化学 2021年12月测量的硝酸盐（mg/L）
作物类型（遥感） c2000, c2010, c2020；Landsat数据，1：城市；2：水域；3：森林；4：果园；5：苗圃；6：草地；7：蔬菜；8：蔬菜-玉米；9：小麦-玉米；10：玉米
作物面积（遥感） 2000_7, 2010_7, 2020_7；2000年、2010年和2020年的蔬菜种植面积（m2）
冬小麦-夏玉米种植面积（m2） 2000年、2010年和2020年的冬小麦-夏玉米种植面积（m2）

农业投入 施肥 2000年；灌溉 2000年；施肥 2010年；灌溉 2010年；施肥 2020年；灌溉 2020年
县统计数据，空间插值；2000年、2010年和2020年的施肥和灌溉情况

地貌与地质 全新世、上更新世、中更新世 1：剥蚀山地/梯田；2：黄土沉积物/相关梯田；3：冲积扇；4：河床；5：洪泛区；6：湖泊和洼地
全新世、上更新世和中更新世的厚度（m）

地下水位深度 WTD2021 2021年的平均地下水位深度（m）
WTD_avg 1980–2021年的平均地下水位深度（m）

气候 avg_pre_y 气象站，年平均降水量（mm）

氧化还原状态 ORP 田间测量，氧化还原电位（mV）

2.3. 机器学习建模框架
实施了一套严格的数据预处理和建模流程，以确保准确识别影响地下水硝酸盐含量的因素。该框架整合了全面的数据清洗、特征选择和先进的机器学习技术，以应对水文地质数据的复杂性。

2.3.1. 数据预处理和特征选择
在建模之前，对数据集进行了系统性的预处理以解决数据质量问题。关键水质参数（T、pH、EC、ORP）中的缺失值使用中值进行插补，而缺失率超过75%的变量被排除在分析之外。表示地质单元的分类变量使用独热编码进行编码，以保留其离散特征。

通过方差膨胀因子（VIF）分析严格评估了预测变量之间的多重共线性，并应用迭代特征消除方法移除了VIF大于10的变量。这一过程确保了保留的特征能够为硝酸盐预测提供独立信息。此外，通过方差阈值识别并移除了常数和接近常数的变量，以提高模型稳定性。

2.3.2. 机器学习方法
采用了三种不同的机器学习算法来预测地下水硝酸盐浓度：随机森林（RF）、极端梯度提升（XGBoost）和反向传播神经网络（BPNN）。选择这些方法是因为它们在捕捉环境系统特有的复杂非线性关系方面具有互补优势。

对每种算法进行了超参数优化，使用随机搜索和5折交叉验证，评估了50-100种参数组合，以高效探索超参数空间。RF配置考虑了估计器数量（50-500）、最大深度（3-50）和最小样本分割（2-20）等参数。XGBoost优化涵盖了学习率（0.01-0.3）和正则化项。BPNN架构搜索包括隐藏层配置（每层50-100个神经元，1-3层）、激活函数和L2正则化参数。

2.3.3. 模型验证和解释
所有模型都使用了空间分层的80/20训练-测试分割，并进行了五次重复实验，以确保统计稳健性。通过决定系数（R2）和均方根误差（RMSE）等指标评估模型性能。通过SHapley Additive Explanations（SHAP）分析增强了模型的可解释性，该分析通过计算所有可能特征组合的边际贡献来量化特征重要性。对于RF和XGBoost模型使用了TreeExplainer，而对于BPNN模型使用了KernelExplainer，从而实现了不同方法之间特征效应的一致解释。这种方法既提供了全局特征重要性排名，也提供了对个别预测的局部解释。

综合建模框架在保持可解释性的同时确保了最佳性能，从而能够准确识别研究区域内控制地下水硝酸盐污染的主要因素。严格的预处理、自动超参数优化和先进的解释方法的结合代表了环境机器学习应用的最新技术。

2.4. 硝酸盐风险指数开发
开发并评估了三种不同的方法来构建硝酸盐污染风险指数。第一种方法是WHO综合风险指数，它建立了预测硝酸盐浓度（[NO??]）和风险得分之间的直接非线性映射。该方法根据WHO标准（50 mg/L）定义风险等级，并使用分段函数强调超过该指导值的严重性。第二种方法是简单线性标准化，通过将模型预测的[NO??]从最小值线性缩放到0-100的风险范围内来提供基准：
(1) RIsimple = 100 × (NO??pred ? NO??min) / (NO??max ? NO??min)

第三种也是最复杂的方法是SHAP-WHO综合风险指数（SWRI），它将模型预测与可解释的人工智能相结合。该方法将基于WHO的风险成分与累积的正SHAP值结合起来，后者量化了每个特征对增加预测硝酸盐浓度的边际贡献。综合指数表示为加权和：
(2) SWRI = α × RIWHO + β × ?total，
其中最优权重α和β通过网格搜索确定，以最大化与实际[NO??]的相关性。这种混合方法不仅反映了绝对浓度，还结合了模型对驱动高风险条件的因素的解释性见解。

2.5. 监测网络优化
本研究开发了一个高级的多目标优化框架，用于通过整合硝酸盐污染风险和空间分布标准来重新设计地下水监测网络。该方法采用了加权K均值聚类算法，根据采样点的空间坐标（经度/纬度）进行聚类，并根据XGBoost-SHAP分析得出的硝酸盐风险指数对聚类过程进行加权。首先将空间坐标投影到Beijing_1954_GK_Zone_20N坐标系（EPSG: 21480）中，以确保等面积表示，从而进行准确的距离计算。

通过综合评估四个关键指标来确定最佳监测点数量：轮廓分数（空间聚类质量）、Calinski-Harabasz指数（簇分离度）、风险解释方差（风险分布的保持）以及使用多种方法（包括核密度估计、最近邻距离、Delaunay三角剖分和网格覆盖分析）计算的新空间均匀性指标。综合优化分数在这些标准之间进行了平衡，权重偏向于风险表示（30%），同时保持空间完整性（聚类质量和均匀性各占25%，空间分布占20%）。

对于每个最优簇，通过自适应策略选择了代表性监测点：在高风险区域（风险指数>70），优先选择风险值最大的点以确保污染热点覆盖；在中等至低风险区域，选择距离簇中心最近的点以保持空间代表性。这种混合方法有效地将网络密度降低了40-60%，同时保持了超过90%的高风险区域覆盖，并保持了原始的风险分布特征，与传统仅基于空间的优化方法相比有显著改进。

优化后的网络在监测效率、风险检测能力和空间覆盖之间实现了最佳平衡，为可持续的地下水质量管理提供了科学上稳健的框架。

3. 结果
3.1. 地下水硝酸盐的空间分布
流域内的地下水硝酸盐浓度表现出显著的空间异质性，从西到东呈梯度分布，西部冲积扇区的浓度较高，而东部洪泛区的浓度较低（图2a）。硝酸盐浓度范围从0.0 mg/L到244.1 mg/L，平均值为33.3 mg/L。大约26.3%的采样井超过了50 mg/L的WHO阈值。值得注意的是，西南部冲积扇区的硝酸盐浓度最高，平均为50.3 mg/L，42.2%的样本超过了WHO限制。该区域记录的最高浓度为244.1 mg/L，这归因于二十年来密集的蔬菜种植和高土壤渗透性，促进了硝酸盐的快速淋溶。相比之下，中部和东部洪泛区的硝酸盐浓度较低（平均：9.0 mg/L），只有3.0%的样本超过了阈值，这可能是由于厚实的非饱和带有效减缓了污染物迁移。然而，在北部洪泛区，由于地下水位较浅（4-20 m）且长期进行蔬菜种植，硝酸盐浓度较高（平均：39.4 mg/L），表明快速渗透和有限的反硝化能力。湖泊洼地通常具有低渗透性和高还原条件，但硝酸盐浓度却较高（平均：33.3 mg/L），28.6%的样本超过了WHO限制。这种差异归因于从冲积扇区通过含水层向湖泊洼地的侧向地下水流动，携带了上游地区的较高硝酸盐负荷。单向ANOVA确认了不同地貌单元之间硝酸盐浓度的显著差异（F = 13.19，p < 0.001），事后测试表明冲积扇区的硝酸盐浓度显著高于洪泛区（p < 0.01）。

有趣的是，一些目前农业活动较低的地区仍然表现出较高的硝酸盐水平，表明历史土地利用的影响持续存在。例如，尽管某些冲积扇区内的地点改种了较不密集的作物，但仍保持较高的硝酸盐浓度，这表明了过去农业活动的遗留影响。这些发现强调了在地下水质量评估中考虑时间动态的重要性，并突出了需要长期监测策略来考虑历史土地管理决策。这些见解对于制定有效的缓解措施和优化监测网络以应对当前和残留污染风险至关重要。

3.2. 机器学习模型性能和SHAP分析
三种机器学习模型在独立测试数据集上的预测性能总结见表2。XGBoost和随机森林（RF）均表现出强大且相当的回归能力，测试R2值分别为0.739和0.740，相应的RMSE分别为17.81 mg/L和17.77 mg/L。这种高一致性突显了基于树的集成方法在捕捉控制地下水硝酸盐浓度的复杂非线性关系方面的稳健性。相比之下，尽管进行了严格的超参数调整并确保了训练收敛（详见补充表S1），反向传播神经网络（BPNN）的表现明显较差，测试R2仅为0.580，RMSE为22.59 mg/L。这种显著的性能差距表明，基于树的集成方法在捕捉这种表格水文地质数据的特定非线性和异质性方面本质上更为合适。

鉴于RF和XGBoost的准确性和几乎相当的准确性，我们使用SHAP（Shapley Additive Explanations）分析来探究模型可解释性并阐明特征重要性。虽然两种模型一致地将历史农业活动和非饱和带属性识别为主要预测因子，但在对关键水文变量的解释上出现了关键分歧，特别是在地下水位深度指标方面。XGBoost-SHAP分析（图3b）产生的排名与水文地质理论更为一致。它将WTD_avg（1980年至2021年的长期平均地下水位深度）赋予了第7高的重要性（平均|SHAP| = 1.07），适当地反映了其作为累积补给历史和非饱和带厚度的稳健代理的作用——这是控制污染物迁移时间和衰减的关键因素。相反，单一年份的测量WTD2021的排名较低（第12位），这与预期一致，即瞬时快照不如多年平均值能更好地反映长期水文过程。相比之下，RF-SHAP分析（图3a）显示这两个变量之间的区分不那么明显，WTD2021获得了较高的重要性得分。这种差异表明，XGBoost的固有正则化和顺序提升机制可能更有效地区分了持久的水文条件和瞬态噪声。因此，由于其更大的物理合理性和解释一致性，选择了XGBoost模型进行后续分析。

2.3.3. 模型解释
通过空间分层的80/20训练-测试分割和五次重复实验评估了所有模型，以确保统计稳健性。性能通过决定系数（R2）和均方根误差（RMSE）等指标进行评估。通过SHAP分析增强了模型的可解释性，该分析通过计算所有可能特征组合的边际贡献来量化特征重要性。对于RF和XGBoost模型使用了TreeExplainer，而对于BPNN模型使用了KernelExplainer，从而实现了不同方法之间特征效应的一致解释。这种方法提供了全局特征重要性排名和对个别预测的局部解释。

综合建模框架在保持可解释性的同时确保了最佳性能，从而能够准确识别控制研究区域内地下水硝酸盐污染的主要因素。严格的预处理、自动超参数优化和先进的解释方法的结合代表了环境机器学习应用的最新技术。

2.4. 硝酸盐风险指数开发
开发并评估了三种不同的方法来构建硝酸盐污染风险指数。第一种方法是WHO综合风险指数，它建立了预测硝酸盐浓度（[NO??]）和风险得分之间的直接非线性映射。该方法根据WHO标准（50 mg/L）定义风险等级，并使用分段函数强调超过指导值的严重性。第二种方法是简单线性标准化，通过将模型预测的[NO??]从最小值线性缩放到0-100的风险范围内来提供基准：
(1) RIsimple = 100 × (NO??pred ? NO??min) / (NO??max ? NO??min)

第三种也是最复杂的方法是SHAP-WHO综合风险指数（SWRI），它将模型预测与可解释的人工智能相结合。该方法将基于WHO的风险成分与累积的正SHAP值结合起来，后者量化了每个特征对增加预测硝酸盐浓度的边际贡献。综合指数表示为加权和：
(2) SWRI = α × RIWHO + β × ?total，
其中最优权重α和β通过网格搜索确定，以最大化与实际[NO??]的相关性。这种混合方法不仅反映了绝对浓度，还结合了模型对驱动高风险条件的因素的解释性见解。

2.5. 监测网络优化
本研究开发了一个高级的多目标优化框架，用于通过整合硝酸盐污染风险和空间分布标准来重新设计地下水监测网络。该方法采用了加权K均值聚类算法，根据采样点的空间坐标（经度/纬度）进行聚类，并根据XGBoost-SHAP分析得出的硝酸盐风险指数对聚类过程进行加权。首先将空间坐标投影到Beijing_1954_GK_Zone_20N坐标系（EPSG: 21480）中，以确保等面积表示，从而进行准确的距离计算。

通过综合评估四个关键指标来确定最佳监测点数量：轮廓分数（空间聚类质量）、Calinski-Harabasz指数（簇分离度）、风险解释方差（风险分布的保持）以及使用多种方法计算的新空间均匀性指标（包括核密度估计、最近邻距离、Delaunay三角剖分和网格覆盖分析）。综合优化分数在这些标准之间进行了平衡，权重偏向于风险表示（30%），同时保持空间完整性（聚类质量和均匀性各占25%，空间分布占20%）。

对于每个最优簇，通过自适应策略选择了代表性监测点：在高风险区域（风险指数>70），优先选择风险值最大的点以确保污染热点覆盖；在中等至低风险区域，选择距离簇中心最近的点以保持空间代表性。这种混合方法有效地将网络密度降低了40-60%，同时保持了超过90%的高风险区域覆盖，并保持了原始的风险分布特征，显示出比传统的仅基于空间的优化方法有显著改进。

优化后的网络在监测效率、风险检测能力和空间覆盖之间实现了最佳平衡，为可持续的地下水质量管理提供了科学上稳健的框架。

3. 结果
3.1. 地下水硝酸盐的空间分布
流域内的地下水硝酸盐浓度表现出显著的空间异质性，从西到东呈梯度分布，西部冲积扇区的浓度较高，而东部洪泛区的浓度较低（图2a）。硝酸盐浓度范围从0.0 mg/L到244.1 mg/L，平均值为33.3 mg/L。大约26.3%的采样井超过了50 mg/L的WHO阈值。值得注意的是，西南部冲积扇区的硝酸盐浓度最高，平均为50.3 mg/L，42.2%的样本超过了WHO限制。该区域记录的最高浓度为244.1 mg/L，这归因于二十年来密集的蔬菜种植和高土壤渗透性，促进了硝酸盐的快速淋溶。相比之下，中部和东部洪泛区的硝酸盐浓度较低（平均：9.0 mg/L），只有3.0%的样本超过了阈值，这可能是由于厚实的非饱和带有效减缓了污染物迁移。然而，在北部洪泛区，由于地下水位较浅（4-20 m）且长期进行蔬菜种植，硝酸盐浓度较高（平均：39.4 mg/L），表明快速渗透和有限的反硝化能力。湖泊洼地通常具有低渗透性和高还原条件，但硝酸盐浓度却较高（平均：33.3 mg/L），28.6%的样本超过了WHO限制。这种差异归因于从冲积扇区通过含水层到湖泊洼地的侧向地下水流动，携带了上游地区的较高硝酸盐负荷。单向ANOVA确认了不同地貌单元之间硝酸盐浓度的显著差异（F = 13.19，p < 0.001），事后测试表明冲积扇区的硝酸盐浓度显著高于洪泛区（p < 0.01）。

有趣的是，一些目前农业活动较低的地区仍然表现出较高的硝酸盐水平，这表明历史土地利用的影响持续存在。例如，尽管某些冲积扇区内的地点改种了较不密集的作物，但仍保持较高的硝酸盐浓度，表明了过去农业活动的遗留影响。这些发现强调了在地下水质量评估中考虑时间动态的重要性，并突出了需要长期监测策略来考虑历史土地管理决策。这些见解对于制定有效的缓解措施和优化监测网络以应对当前和残留污染风险至关重要。

3.2. 机器学习模型性能和SHAP分析
三种机器学习模型在独立测试数据集上的预测性能总结在表2中。XGBoost和随机森林（RF）均表现出强大且相当的回归能力，测试R2值分别为如图4b所示，与土壤相关的特征贡献从山前近缘区域的强烈正相关（距离山体前沿约50公里范围内）转变为远端洪泛区的负相关，这反映了从粗粒、高渗透性的扇形沉积物向细粒、低渗透性的洪泛区沉积物的转变。在50公里以外，正相关性的急剧下降突显了有利渗滤条件的空间局限性。图4c展示了Cl?的SHAP依赖性，显示出随着氯化物浓度的增加而呈现明显的正趋势。这种共变表明NO??和Cl?具有共同的人为来源（例如，农业回流水、生活污水），尽管中等程度的相关性表明它们在进入地下水位带后的衰减机制可能不同，这可能归因于不同地下水位带环境中的反硝化潜力差异。图4d展示了d_qh的影响：较厚的全新世沉积物（>15米）显示出负的SHAP值，表明硝酸盐迁移受到显著延迟，这可能是由于在细粒沉积物中通过吸附或反硝化等过程增强了衰减作用。图4f和4g展示了农业驱动因素的阈值效应。两组比较（高于阈值与低于阈值）证实，对于年灌溉量超过300毫米的情况，灌溉的正面贡献在统计上显著（p<0.05）。同样，对于氮负荷，当年负荷量超过265公斤/公顷时，其正面影响更为显著（p<0.05），并且当负荷量超过350公斤/公顷时，污染生成的加速效应更加明显，这突显了非线性和临界点的存在。最后，图4h通过组间比较显示，只有当地下水位深度超过16米时，地下水位带（WTD_avg）的保护功能才在统计上变得显著（p<0.05），这强调了足够的地下水位带厚度在减轻污染物流入含水层方面的关键作用。

下载：下载高分辨率图像（599KB）
下载：下载全尺寸图像

图4. SHAP总结（a）以及经度（b）、Cl?（c）、d_qh（d）、2010_9（e）、Irrigation_2000（f）、Fertilization_2000（g）和WTD_avg（h）的SHAP依赖性图。注意：每个特征变量的含义可以在表1.3.3中找到。

硝酸盐风险指数和空间风险模式
使用Combined SHAP-WHO方法（α = 0.5）构建的风险指数是一种非常有效的工具，用于分层与硝酸盐污染相关的健康风险。该指数与实际硝酸盐浓度表现出强烈的正相关性（R = 0.8822，图5a），验证了其在反映暴露严重程度方面的准确性。随后，根据预定的阈值将指数值在研究区域内的分布分为五个风险等级（图5c），与WHO标准对齐：非常低（0?10）、低（10?20）、中等（20?40）、高（40?70）和非常高（>70）。

下载：下载高分辨率图像（381KB）
下载：下载全尺寸图像

图5. SHAP-WHO综合风险指数（SWRI）的全面评估。(a) SWRI与观测到的地下水硝酸盐浓度之间的相关性；(b) 五个级别的硝酸盐污染风险类别的分布；(c) 预测的硝酸盐浓度与SWRI之间的关系；(d) 实际硝酸盐浓度与WHO硝酸盐阈值（50毫克/升）的超出情况。
实际值与预测的WHO标准超出值之间的高度一致性证明了该模型在识别危险区域方面的可靠性，而风险指数统计数据显示平均风险水平为中等（21.60），最高值达到105.50，表明存在显著的污染热点。这种分类揭示了风险的关键空间分布。大约62.7%的区域（包括非常低和低风险类别）处于较低风险水平。中等风险类别最为普遍，占总面积的32.15%（5343.0平方公里）（图6）。

下载：下载高分辨率图像（319KB）
下载：下载全尺寸图像

图6. 白羊淀湖流域硝酸盐污染风险水平的地理空间分布。
至关重要的是，16.1%的区域（包括高和非常高风险类别，总计2675.5平方公里）被认定为面临较高风险，其中非常高风险区域（1.01%，168.1平方公里）在空间上较为集中，风险地图显示这些区域围绕特定的水文特征聚集。结果强调，虽然大部分区域的风险较低，但对识别出的高风险区域进行有针对性的监测和干预策略是必要的，因为这些区域对公共健康构成了重大威胁。该风险指数为政策制定者提供了一个科学依据、空间明确的框架，以便优先分配资源用于环境保护和公共卫生保障。

3.4. 优化监测网络
多目标优化框架的应用成功重新设计了硝酸盐污染的地下水监测网络，在监测效率和风险评估能力之间取得了更好的平衡。优化过程结合了空间分布和基于SHAP的硝酸盐风险指数，确定了最佳簇数（k）为82，使得监测点数量从原来的278个减少了70.5%。
k=82的选择是通过多项指标的严格评估确定的，详细信息见图7。综合优化得分在0.694时达到峰值，这个最佳点代表了帕累托最优解，有效平衡了空间覆盖范围和风险敏感性之间的权衡。具体来说，虽然空间聚类质量（轮廓分数）随着k的增加而自然下降，但风险覆盖范围在达到某个点后显示出特征性的增加，表明超过最佳点后收益递减。综合得分与风险覆盖范围之间的强正相关（r=0.93）强调了该框架在优先检测污染热点方面的成功。相反，与轮廓系数的显著负相关（r=-0.75）证实了在空间紧凑性和基于风险的代表性之间的预期折中。

下载：下载高分辨率图像（498KB）
下载：下载全尺寸图像

图7. 通过多目标优化框架确定最佳地下水监测点数量。(a) 基于综合得分的最佳簇数（k）选择；(b) 随簇数（k）变化的风险覆盖范围；(c) 空间聚类质量（轮廓分数）与簇数（k）；(d) 空间分布均匀性与簇数（k）；(e) 不同簇数（k）的标准化优化指标；(f) 优化指标的相关矩阵。
优化网络的有效性在图8中得到了直观展示。原始网络（278个点）显示出密集且不规则的分布。相比之下，优化后的网络（82个点）在空间分布上更加均匀，同时准确保留了高风险区域的地理分布（风险指数>70）。重要的是，优化后的网络实现了对原始七个非常高风险点的100%覆盖，确保了对最关键区域的持续监控。最终网络保持了稳健的风险表示，平均风险指数为34.94（高于原始网络的21.60），高风险点占比8.5%，中等风险点占比47.6%，低风险点占比43.9%。这表明优化后的网络不仅在成本效益上更优，而且在战略上也更加集中，增强了主动管理地下水质量的能力。

下载：下载高分辨率图像（181KB）
下载：下载全尺寸图像

图8. 原始（a）和优化后的（b）地下水监测网络的空间比较。

4. 讨论
本研究利用了空间广泛的数据集和机器学习来解码白羊淀湖流域地下水硝酸盐污染的复杂、历史驱动的模式。结合多十年的历史数据，可以对其中的时空机制进行有力的推断。以下讨论综合了研究结果，探讨了潜在的过程，将其置于现有文献的背景下，并讨论了我们方法的影响和局限性。

4.1. 历史农业实践的遗留效应
本研究的一个关键发现是当前农业管理与当代地下水质量之间存在显著的时间断层。机器学习模型一致地将历史变量，特别是Irrigation2000和2010年的冬小麦-夏玉米轮作面积（2010_9），作为顶级预测因子，其表现显著优于2020年的对应变量。这一结果提供了强有力的数据驱动证据，表明2021年测量的硝酸盐浓度在很大程度上是十年前密集农业实践的遗留效应。为了有效排除土地利用惯性的混淆影响，对这三个十年间的农业投入进行了时间相关性分析（补充图S1）。虽然相邻十年在空间上具有很强的连续性，但2000年和2020年投入之间的相关性仅处于中等水平（例如，灌溉的相关性r=0.44，施肥的相关性r=0.52）。这种时间上的差异证实了历史变量保留了独特的空间变异。因此，模型对2000年代初预测因子的明显偏好代表了真实的物理时间滞后，而不是共线性的统计假象。这种滞后是该地区固有的缓慢地下传输过程的直接后果。高硝酸盐浓度与2000年代的历史蔬菜种植强度之间的显著正相关（r=0.82，p<0.001），通过空间自相关分析得到了直观的证实，进一步强调了这种遗留效应。这一发现与NCP中的地下水时间滞后研究结果一致，这些研究指出停留时间超过10年，通常跨越数十年（Liu等人，2025年；Zhang等人，2023年）。2000年代初的特点是漫灌和氮利用效率低，导致硝酸盐的大量渗漏。这种信号的持续存在凸显了含水层系统的水文惯性，并对管理提出了重大挑战：仅关注当前肥料使用的政策可能需要数十年才能在水质改善上见到效果。在SHAP依赖性图中观察到的灌溉和施肥的非线性、阈值样响应进一步表明存在临界点（灌溉为300毫米/年，施肥为265公斤/公顷），超过这些临界点后，土壤-地下水位带系统的自然衰减能力会被超出，从而形成长期的污染源。

4.2. 地下水位带的异质性和硝酸盐传输
地下水位带特征的重要性，特别是全新世（d_qh）和晚更新世（d_q3）沉积物的厚度，证实了它们对硝酸盐命运和传输的关键控制作用。SHAP分析显示，这些地质特征主要作为重要的物理缓冲层，引入了旅行时间的显著空间变异性。较厚的全新世沉积物（>15米）的负SHAP值表明了强烈的保护效应，这主要是由于低渗透性细粒沉积物内的物理延迟作用，显著限制了污染物的向下流动（Liu等人，2022年）。然而，这种储存能力也可能成为一种负担，因为这些区域可能会积累历史遗留的硝酸盐，并随着时间的推移逐渐释放，解释了为什么当前输入量适中的区域仍然可能出现高含水层污染。氯化物（Cl-）和硝酸盐的协同解释进一步揭示了机制上的见解。它们在SHAP依赖性图中的强正趋势表明它们具有共同的人为来源，例如农业回流水。然而，它们之间的中等相关性意味着地下路径的不同，硝酸盐的衰减程度取决于异质地下水位带中的氧化还原条件（ORP）和碳的可用性（Zhang等人，2024年）。这一点得到了ORP总体特征重要性的支持（补充表S2），表明生化去除在高度异质的地下水位带中起着次要但可观察的作用。这种复杂性在传统的脆弱性评估中常常被简化。例如，长期平均地下水位深度（WTD_avg）的显著保护作用仅在超过16米时才变得明显，这突显了像DRASTIC这样的静态指数的一个关键局限性，它可能会错误地将深层地下水位带分类为低风险，而没有考虑其污染负荷的历史（Baalousha，2025年）。

4.3. 方法论综合与应用
机器学习模型的比较分析显示，尽管随机森林（RF）和XGBoost的预测准确性几乎相同（测试R2分别为0.740和0.739），但XGBoost提供了更符合水文地质学的解释（Gholami等人，2025年）。关键的区别在于处理地下水位深度指标的方式：XGBoost正确地将长期平均值（WTD_avg）视为比单一年份快照（WTD2021）更为重要，这与持续的水文条件比瞬态状态更相关的概念一致。这突显了像SHAP这样的模型可解释性工具的价值，它超越了黑箱预测，生成了关于系统行为的可测试的物理假设（Wang和Peng，2024年）。这种机制理解直接转化为实用的管理工具。SHAP-WHO综合风险指数（SWRI）的开发是一个重要的进步。通过将WHO指南值与累积的正SHAP值（φ?total）相结合，SWRI不仅捕捉到了预测的浓度，还捕捉到了驱动风险的因素。这种混合方法产生了一个强大的工具，其与观测到的硝酸盐水平有很强的相关性（r = 0.882），有效地将流域分层为不同的风险类别。由此产生的风险地图显示，虽然62.7%的区域处于低风险，但令人担忧的16.1%（2675.5平方公里）处于高或非常高风险，为优先干预提供了精确的空间模板。监测网络的多目标优化展示了这种基于风险的方法的实际应用。通过将SWRI集成到加权K-means聚类算法中，网络规模减少了70.5%（从278个点减少到82个点），同时实现了对原始高风险区域的100%覆盖。关键的是，我们的多目标优化通过严格实施空间均匀性约束，减轻了静态风险指标的局限性。将优化后的网络（图8b）与区域流场（图1b）对齐，证实了代表性点在低风险的东部洪泛区得到了有效保留。这些点作为下游的哨兵井，确保了能够及早检测到来自上游污染区的平流羽流的迁移。这种方法解决了环境监测中的一个关键问题，从基于便利性或均匀分布的设计转向了一种战略性、成本效益高的框架，以最大化检测污染热点（Feng等人，2025年；Fu等人，2025年）。

4.4. 局限性、影响和未来方向
本研究的主要局限性在于它依赖于空间快照来推断时间过程。虽然使用历史土地利用数据可以强烈暗示遗留效应，但由于缺乏多时间尺度的地下水硝酸盐测量或同位素数据（例如δ1?N-NO??），无法确定污染源的分配和进行直接年龄测定。此外，尽管地质数据详细，但这些数据来自区域调查，由于地下异质性的存在而引入了不确定性。另一个关键局限性源于XGBoost框架本质上是基于点的。该模型假设地下水质量主要受垂直渗透控制，忽略了横向平流传输和上游溶质负荷的影响。因此，模型可能会显示出空间结构化的残差，可能低估了区域排放区（例如湖泊洼地）的实际硝酸盐水平，因为在这些地方横向地下水流动补充了局部垂直渗透。这些局限性为未来的研究指明了方向。首先，同位素分析对于验证推断的污染源和传输时间至关重要。其次，为了解决基于点的机器学习的拓扑限制，未来的研究应该将机器学习得出的风险指标与基于物理的分布式模型（例如用于垂直传输的HYDRUS-1D和用于横向对流通量的MODFLOW/MT3DMS）结合起来，以创建一个更全面的预测框架。最后，通过实施带有连续传感器的优化网络，可以将其转变为一个自适应的早期预警系统。

从政策角度来看，研究结果表明，有效的地下水保护需要从根本上转向管理遗留污染。这种转变应包括在我们风险地图（图6）中划定的高风险区域采取有针对性的干预措施，例如推广控释肥料或将高投入作物转换为低强度的轮作制度。这里开发的框架——用于理解遗留问题、评估动态风险和优化监测——为全球农业流域实现可持续的地下水管理提供了一个科学上可靠且可转移的蓝图。

5. 结论
本研究通过整合多十年的农业数据、非饱和带特征和可解释的机器学习方法，阐明了典型农业流域中地下水硝酸盐污染的时空机制。我们的研究发现表明，当前的地下水硝酸盐浓度是历史农业活动的遗留结果，并且存在显著的多十年时间滞后。选择XGBoost-SHAP框架是因为其比随机森林和BPNN具有更好的物理合理性，它识别出了历史农业活动和非饱和带特性作为主要控制因素。重要的是，21世纪初的变量（例如Irrigation2000和Fertilization2000）在预测2021年的硝酸盐浓度方面名列前茅，突显了明显的遗留效应。分析还确定了关键的非线性阈值：年灌溉量超过约300毫米/年和氮肥施用量超过265公斤/公顷/年会成为硝酸盐渗漏风险急剧增加的转折点。这表明，历史上的高强度农业活动造成了长期的污染负担，这种负担将在未来几十年内影响地下水质量。

非饱和带被证实是一个关键且异质性的调节器，而不仅仅是简单的通道。全新世沉积物的厚度（d_qh）和长期平均地下水位深度（WTD_avg）成为关键控制因素，其中地下水位深度超过16米在硝酸盐衰减方面具有统计显著性。开发的SHAP-WHO综合风险指数（SWRI）有效地将这些机制洞察转化为一个空间明确的风险地图——该地图与观测到的硝酸盐浓度高度相关，并准确划分了高风险区域，尤其是在东部洪泛区。最后，应用风险加权的多目标空间聚类算法优化了监测网络，减少了70.5%的监测站点数量，同时确保了对高风险区域的全面覆盖。这反映了从空间均匀监测向基于风险的自适应监测的战略性转变。

最终，这项研究表明，有效的地下水治理必须应对历史农业实践造成的长期污染负担，而不仅仅是关注当前的输入。我们整合的框架推动了一个关键的战略转变：从被动的、空间均匀的采样转向主动的、基于风险的监测。通过明确划分多十年的遗留热点并考虑非饱和带的缓冲效应，本研究为政策制定者提供了一个科学上可靠的蓝图。它不仅优化了区域监测资源，还指导了有针对性的、特定地点的干预措施（例如定制的作物轮作和局部肥料限制），为全球集约化农业流域的可持续地下水保护提供了可行的路径。

**作者贡献声明**
王世勤：写作——审稿与编辑、监督、项目管理、资金获取。
张志雄：写作——审稿与编辑、原始草稿撰写、可视化、验证、方法论、调查、数据管理、概念化。
周敬雄：写作——审稿与编辑。
张子龙：可视化、调查。
刘炳霞：项目管理、调查。