整合水化学与可解释机器学习用于土耳其比斯米尔平原地下水质量评估

《Water》：Integrating Hydrochemistry and Explainable Machine Learning for Groundwater Quality Assessment in the Bismil Plain, Türkiye

　　

本研究利用2022年至2024年期间从26口井采集的208个样本（涵盖8个季节采样周期），评估了比斯米尔平原（土耳其东南部迪亚巴克尔省）的地下水质量。在每个样本中，分析了pH、电导率（EC）以及主要离子Ca²⁺、Mg²⁺、Na⁺、K⁺、Cl^?、SO₄^2?、HCO₃^?和NO₃^?，并基于世界卫生组织（WHO）标准计算了每个观测点的水质指数（WQI）。研究将经典水化学解释方法（包括描述性统计、层次相关分析、方差膨胀因子、Piper图和Gibbs图）与可解释机器学习框架相结合，该框架将基于SHAP的特征选择整合到随机森林（Random Forest）、XGBoost、支持向量回归（SVR）和堆叠集成（stacking ensemble）模型中。此外，利用Moran's I和普通克里金（ordinary kriging）评估了空间残差，通过孤立森林（Isolation Forest）和局部异常因子（Local Outlier Factor）算法识别了异常值，并使用bootstrap重采样量化了预测不确定性。WQI值范围为79.37至125.48（均值99.67），所有样本仅被归类为“良好”（49.5%）和“较差”（50.5%）两个质量等级，表明含水层已接近水质临界阈值。在空间上，最高（最差）WQI值在平原西南部和中部形成连贯区域，而中东部井的WQI值最低；四种模型均再现了这一空间格局。XGBoost和堆叠集成模型表现出可比的预测性能（R²分别为0.911和0.910；均方根误差RMSE分别为3.29和3.27），而SHAP分析将EC识别为最主要的控制因素，其次为NO₃^?、SO₄^2?、Ca²⁺、Mg²⁺和Cl^?（平均|SHAP|值分别为6.86、1.57、1.10、1.09、0.85和0.72个WQI单位）。对XGBoost和堆叠集成模型的残差场计算的Moran's I分别为-0.067（p=0.275）和-0.068（p=0.273），因此对这些残差进行普通克里金后几乎未产生校正，而SVR残差仍存在空间自相关（I=0.242；p=0.001），并通过地统计步骤得到了有意义的校正。本研究的创新之处在于，将可解释机器学习、地统计残差分析、异常检测和基于bootstrap的不确定性评估整合到一个统一框架中，用于多季节地下水数据集，同时利用基于Moran's I的方法评估了空间校正的有效性。

本研究以《Water》期刊发表的论文《Integrating Hydrochemistry and Explainable Machine Learning for Groundwater Quality Assessment in the Bismil Plain, Türkiye》为基础，系统评估了土耳其东南部比斯米尔平原（Bismil Plain）的地下水质量。该平原位于底格里斯河流域，是迪亚巴克尔省（Diyarbak?r）最集约化的农业区之一，约95%的土地用于农业，主要种植谷物和棉花，依托GAP灌溉基础设施。半干旱气候与灌溉回流水、化肥施用共同导致盐度和硝酸盐在含水层中积累。此前已有研究对该平原的地表水进行了表征，揭示了盐度胁迫和高SO₄^2?、Cl^?浓度，并关联了健康风险，但对该平原含水层系统的多季节地下水质量评估尚属空白。现有区域水质指数（WQI）研究存在两个主要局限性：一是许多研究仅对比机器学习模型精度，而未量化驱动WQI的关键输入变量；二是混合机器学习与残差克里金（residual kriging）的框架常缺乏对残差空间结构的检验（如Moran's I），导致地统计步骤可能无效。因此，本研究旨在填补这些空白，通过整合水化学解释、基于SHAP的可解释机器学习、地统计残差分析、异常检测和bootstrap不确定性评估，建立统一分析框架，用于评估多季节地下水数据集，并明确检验残差克里金是否真正提升预测精度。

研究人员在2022年至2024年间，从26口井采集了208个地下水样本，覆盖8个季节周期（每年冬、春、夏、秋）。井深估计为135–269米，平均约200米，抽取深层盆地含水层（包括碳酸盐和蒸发岩单元）。每份样本分析pH、电导率（EC）及主要离子（Ca²⁺、Mg²⁺、Na⁺、K⁺、Cl^?、SO₄^2?、HCO₃^?、NO₃^?），并计算基于WHO的加权算术WQI。研究结合经典水化学方法（描述性统计、层次相关分析、方差膨胀因子VIF、Piper图和Gibbs图）与可解释机器学习框架：首先用随机森林（Random Forest）计算SHAP值，基于平均绝对SHAP贡献超过中位数筛选出6个关键预测因子（EC、Ca²⁺、Mg²⁺、Cl^?、SO₄^2?、NO₃^?）；然后训练随机森林、XGBoost、支持向量回归（SVR）和堆叠集成（stacking ensemble）四种模型，采用空间块交叉验证（3×3网格）和80%/20%训练/测试划分；之后对残差进行Moran's I检验（基于K=5最近邻权重矩阵），对存在显著空间自相关的残差拟合球状变差函数并进行普通克里金；同时采用孤立森林（Isolation Forest）和局部异常因子（LOF）进行异常检测（取两者交集），并通过100次bootstrap重采样量化预测不确定性。所有计算在Python 3中完成，使用scikit-learn、XGBoost、SHAP、Optuna、libpysal、esda和PyKrige等库。

**研究结果**
**3.1 水化学特征**：描述性统计显示，平均离子总和为1486 mg/L，对应总溶解固体约1200–1300 mg/L，表明含水层为中等矿化。平均pH 7.61，EC均值1499 μS/cm，最大2174 μS/cm。NO₃^?变异系数最高（40.87%），说明氮输入空间异质性大；7个样本（3.4%）超过WHO 50 mg/L限值，108个样本（51.9%）超过25 mg/L的欧盟指导值。季节均值显示EC和Cl^?在夏季最高（约1599 μS/cm和315 mg/L），冬季最低（约1425 μS/cm和278 mg/L），而HCO₃^?在春季峰值（约472 mg/L）。年际差异可忽略（Y1和Y2平均WQI分别为99.82和99.52），表明季节信号主导。空间分布上，高EC、Cl^?、SO₄^2?聚集在西南部和中部，东部矿化度较低。

**3.2 相关结构与多重共线性**：层次聚类相关系数矩阵显示EC与SO₄^2?、Cl^?、Ca²⁺、Mg²⁺、HCO₃^?强相关（r=0.55–0.86），NO₃^?与WQI弱相关（r=0.19）。VIF值均低于4（EC最高3.95），未超过阈值5，因此所有变量可保留。

**3.3 水化学相与主导过程**：Piper图显示样本落在混合Ca²⁺–Mg²⁺–Cl^?–SO₄^2?区域，夏季向氯化物端元漂移。Gibbs图表明岩石-水相互作用为主，夏季向蒸发-结晶端元偏移。平均摩尔Ca²⁺/SO₄^2?比为1.36，指示石膏溶解叠加碳酸盐钙源。

**3.4 WQI分布与分类**：WQI范围79.37–125.48，均值99.67。49.5%样本为“良好”（50–100），50.5%为“较差”（100–200），无“优秀”或“极差”样本。夏季平均WQI 102.34，冬季96.83，季节性偏移使部分样本跨过100界限。

**3.5 机器学习模型性能**：测试集上，XGBoost和堆叠集成表现最佳（R²分别为0.911和0.910，RMSE 3.29和3.27，MAE 2.38和2.36），随机森林次之（R²=0.859，RMSE=4.16），SVR最差（R²=0.819，RMSE=4.70）。但模型间差异在bootstrap置信区间内，可视为可比。关键区别在于残差空间结构：XGBoost、堆叠集成和随机森林的Moran's I为负且不显著（p>0.17），而SVR残差Moran's I=0.242（p=0.001），仍存在空间自相关。

**3.6 空间残差分析**：对XGBoost和堆叠集成，残差克里金几乎无校正（残差场接近零值）；对SVR，球状变差函数拟合有效（基台值14.78，偏基台4.79，变程约0.48°），克里金后产生有意义的校正。所有模型的空间预测图均显示西南部和中部为高WQI（差水质）区域，中东部最低。

**3.7 基于SHAP的驱动因子归因**：全局SHAP分析显示，EC贡献最大（平均|SHAP|=6.86 WQI单位），其次为NO₃^?（1.57）、SO₄^2?（1.10）、Ca²⁺（1.09）、Mg²⁺（0.85）、Cl^?（0.72）。EC单独贡献约为其余五个因子之和的四倍。NO₃^?虽与WQI线性相关弱（r=0.19），但通过高浓度阈值效应在SHAP中排名第二。

**3.8 异常样本**：孤立森林（固定污染率0.10）标记21个全局异常，LOF（数据驱动阈值）标记4个局部异常。两者交集（同时标记）对应EC、Cl^?、SO₄^2?或NO₃^?接近上十分位的样本，主要位于夏季中部集约农业区，表明农业负荷与季节蒸发浓缩在时空上重合。

**3.9 预测不确定性**：bootstrap 95%置信区间平均宽度为6.97 WQI单位，约为领先模型RMSE的两倍。最宽区间对应高NO₃^?异常和夏季高EC样本，即非线性效应最强的样本，提示需对这些井进行针对性重采样。

**3.10 水化学-机器学习综合解释**：经典水化学与数据驱动分析一致表明，含水层受碳酸盐风化、石膏/硬石膏溶解和季节性蒸发浓缩共同作用，无单一主导过程。SHAP分析将EC列为系统性驱动因子，NO₃^?为局部高影响因子，为管理提供了可操作信号：盐度需网络级管理，硝酸盐需热点干预。

**讨论与结论**
讨论部分总结了水化学驱动因子与过程：Bismil含水层为中等矿化，混合碳酸盐-硫酸盐系统，石膏/硬石膏溶解提供SO₄^2?，Cl^?部分来自人为输入（肥料、灌溉回流水）。NO₃^?作为局部农业污染标志，其高变异性和阈值效应解释了SHAP排名第二的原因。可解释机器学习的重要性在于：四种模型精度相近，但SHAP提供了超越线性相关的驱动因子排序，EC主导，NO₃^?为高影响局部变量。残差克里金的有效性通过Moran's I检验判定：仅当残差存在空间自相关时（如SVR），地统计校正才有意义；对于XGBoost和堆叠集成，残差无结构，校正无效。因此推荐未来研究将Moran's I作为常规预检。与区域先前工作相比，本研究结果与地表水盐度压力模式一致，但地下水响应时间尺度更长。硝酸盐阈值效应与当地化肥施用局部化特征相符。研究局限性包括：网络规模有限（26井），缺少微量元素和稳定同位素，SHAP特征降维可能略压缩WQI信息，阴离子轻微过量（平均CBE=-3.74%）可能来自未分析物种。未来需延长记录、整合同位素和微量元素，并开展贝叶斯空间模拟。

**研究结论**（翻译自5. Conclusions）：
通过经典水化学解释与可解释机器学习、残差地统计、异常筛选和bootstrap不确定性的结合，对比斯米尔平原地下水质量进行了评估。2022至2024年间从26口井采集的208个样本（覆盖8个季节）的WQI值处于狭窄范围（79.37–125.48），所有样本仅属于“良好”或“较差”类别，表明含水层处于这两个类别之间的边界而非严重退化状态。WQI变异的主要系统性驱动因子为EC、SO₄^2?和Cl^?，而NO₃^?作为高影响但局部集中的控制因子：7/208个样本超过WHO 50 mg/L限值，108/208个样本超过25 mg/L的欧盟指导值，因此硝酸盐管理需采取热点靶向策略而非网络级策略。SHAP返回的综合排序（EC > NO₃^? > SO₄^2? > Ca²⁺ > Mg²⁺ > Cl^?）同时捕捉了盐度驱动的系统性信号和硝酸盐的阈值效应，而线性相关分析单独会遗漏后者。方法学上，本文为小样本、多季节地下水数据集构建了透明流水线，整合了可解释ML、残差克里金、异常筛选和bootstrap不确定性。关键方法学信息是：残差克里金仅当底层ML模型未解释空间结构时才有意义；对于XGBoost和堆叠集成基线，残差场实际上无结构，混合步骤产生近乎零的校正。作者直接报告了这一结果，而非将其视为统一的空间改进，并建议任何未来的ML加克里金研究将Moran's I作为常规预检。对比斯米尔平原地下水管理，提出三项操作性建议：第一，夏季应优先监测已接近WQI阈值100的井的EC、Cl^?和SO₄^2?，因为季节性蒸发浓缩足以将边界样本推入“较差”类别；第二，通过联合孤立森林–LOF筛选识别的NO₃^?热点应进行实地调查，而非通过网络级化肥政策解决；第三，应在下一个监测周期向现有参数集添加同位素和微量元素分析，以明确区分自然和人为负荷。更广泛而言，本文描述的工作流程可移植到其他小样本农业含水层，在这些情况下，决策需基于有限但结构良好的数据，且透明的驱动因子归因比排行榜精度更有用。