干旱区生态系统服务与地下水储量异常的非线性关联及候选统计阈值：基于CatBoost–SHAP的分析

《Water》：Nonlinear Associations and Candidate Statistical Thresholds Between Ecosystem Services and Groundwater Storage Anomalies: An Analysis Using CatBoost–SHAP

　　

： 干旱区生态修复可增强生态系统功能，但可能伴随地下水压力。研究人员以塔里木河流域九大源流为研究区，分析了2003–2025年基于GRACE/GLDAS的地下水储量异常（GWSA）代理指标，以及2000–2025年四种生态系统服务与综合指数（ESI）；联合分析时段为2003–2025年。通过相关分析、联合趋势分区及CatBoost–SHAP方法，刻画了空间错配与非线性关联。GWSA代理指标以?0.31 mm month?1的速率下降，2020年后负异常加剧，并呈现西北减少、东南增加的空间格局。与地下水水位异常的对比显示，阿克苏站一致性中等（r = 0.52），而策勒站一致性弱且不显著（r = ?0.08），表明该代理指标可靠性存在空间异质性。碳储量、生境质量、净初级生产力及ESI均呈增加趋势，而防风固沙量下降。ESI增加而GWSA下降的区域占研究区面积的42.5%。五个CatBoost模型的决定系数（R2）介于0.868–0.977之间。就ESI而言，SHAP识别出约?18.6 mm处的候选转变点，低于该值时GWSA与负向模型贡献相关联；人类足迹贡献随干扰增强而持续转负。上述统计关联而非因果效应表明，生态系统服务的改善并不必然意味着地下水的恢复。

论文解读

一、研究背景与问题提出

干旱区约占全球陆地面积的40%，是受水资源约束最为强烈、生态退化风险最高的陆地系统之一。植被生产、土壤稳定性、生物多样性与人类生计均高度依赖有限且波动显著的水资源供给。在全球变暖、干旱加剧与高强度土地利用的叠加影响下，干旱区植被覆盖、土壤肥力及生态系统多功能性可能跨越临界干旱阈值后发生阶段性、非线性退化。过去二十年间，中国大规模生态修复与土地管理推动了广泛植被绿化，但干旱区降水稀少、潜在蒸散强烈，植被恢复同时改变了降水在蒸散发、径流、土壤水与地下水补给间的分配，可能产生生态功能增益与水资源消耗之间的权衡。因此，干旱区生态管理的核心问题已不再仅仅是植被是否增加，而是多种生态系统服务能否在区域水资源承载力范围内得以维持。

地下水是连接山区径流、绿洲生产与荒漠生态系统的关键水文纽带。全球地下水依赖生态系统制图表明，约53%具有地下水趋势数据的地下水依赖生态系统位于地下水下降区域。植被对地下水变化的响应可能呈现临界转变、恢复滞后与时间延迟等特征，使同期线性相关难以充分表征其内在过程。与此同时，不同生态系统服务在形成机制与水分敏感性上存在差异，可能表现出协同、权衡或空间错配。GRACE与GRACE-FO卫星任务为观测大尺度陆地水储量变化提供了关键基础，减去GLDAS模拟的土壤水、雪水与冠层水储量异常后，可估算区域地下水储量异常，但该残差可能仍包含未去除的冰川与地表水储量变化，不能直接解释为地下水深度或纯地下水储量变化。

塔里木河流域具有高山冰雪融水、绿洲灌溉需水与强烈荒漠蒸散发共存的特有水循环结构，其九大源流水系在水源供给、绿洲开发强度、土地利用与植被类型上差异显著。已有研究记录了该流域地下水储量重建、生态输水、生态系统服务权衡与环境响应等方面的显著变化，但在GRACE系统原生空间覆盖范围内，对九大源流生态系统服务与地下水之间联合趋势、空间错配、滞后响应及非线性关联的综合分析仍较为有限。本研究正是基于上述背景展开，旨在提升对干旱区内陆河流域生态系统服务演变与大尺度地下水储量变率之间统计关系的理解，为水资源承载力约束下的绿洲保护与荒漠边缘管理提供支撑。

二、研究内容与主要结论

本研究聚焦塔里木河流域九大源流水系，整合多源遥感、气象、水文、土地利用与人类活动数据，评估2000–2025年生境质量（HQ）、碳储量（CS）、防风固沙量（SR）、净初级生产力（NPP）及生态系统服务综合指数（ESI）；2003–2025年月尺度GRACE/GLDAS地下水储量异常（GWSA）代理指标用于时间分析，2003–2025年GWSA与生态系统服务联合年分析用于相关、趋势分区与CatBoost建模。研究得出以下主要结论：2003–2025年GWSA代理指标以?0.3095 mm month^?1线性下降，2020年后负异常加剧，下降集中于西北部；2000–2025年CS、HQ、NPP与ESI总体增加，而SR波动并总体下降。联合趋势分区显示，ESI增加而GWSA下降的区域占42.5%，两者同增区域占26.3%，同降区域占19.7%，GWSA增加而ESI下降区域占11.5%。五个CatBoost模型测试集R²介于0.868–0.977。SHAP依赖分析识别出ESI预测中GWSA约?18.6 mm、NDVI约0.13、人类足迹指数约4.2的候选转变点。上述统计关联而非因果效应表明，生态系统服务改善并不必然指示地下水恢复。

三、关键技术与方法

研究采用以下主要技术方法：基于GRACE/GRACE-FO五种解算方案（CSR、GFZ、JPL球谐解及CSR、JPL mascon解）的月尺度陆地水储量异常（TWSA），对2017年7月至2018年5月任务间隙采用确定性全条件设定（FCS）线性回归重建，减去GLDAS Noah 2.1版土壤水、雪水当量与冠层水储量分量后获得GWSA代理指标；以国家生态系统科学数据中心两口监测井（阿克苏、策勒）地下水水位异常进行时间一致性验证；采用InVEST模型评估碳储量与生境质量，基于修正风蚀方程评估防风固沙量，基于CASA模型估算净初级生产力，经min–max归一化等权集成ESI；采用Sen斜率估计趋势并划分四种联合趋势类型；采用CatBoost回归模型拟合ESI与12个预测因子（气候、地下水、植被、地形、土壤、人类活动）的非线性统计关系，以空间分块80%训练、20%独立测试，五折交叉验证与贝叶斯优化调参，并以随机森林与XGBoost进行算法敏感性检验；采用SHAP TreeExplainer计算特征贡献、依赖关系与交互值，识别候选转变点与交互模式。

四、研究结果

4.1 GWSA代理指标与生态系统服务的时空演变

2003–2025年，TWSA与GWSA代理指标均呈总体下降趋势，GWSA以?0.3095 mm month^?1速率下降，2020年后负异常加剧。月尺度上，TWSA与GWSA自冬至春增加，5月前后达年最大值，夏季下降，晚秋至初冬达最小值。空间上，2002年4月至2006年12月期间正异常局地位于北部，负异常主导南部与西南部；2021–2025年负异常在北部与西北部扩展加剧，东南部部分区域转正。趋势分析显示显著下降集中于北部与西部，局部速率超过25 mm year^?1，显著增加主要位于东南部，最大速率约14 mm year^?1。盆地尺度上，阿克苏、迪那、渭干–库车下降最强，且末与车尔臣上升。与监测井对比显示，阿克苏站足迹均值GWSA与地下水水位异常（GWLA）中度相关（Pearson r = 0.524，p < 0.001），策勒站同期相关性弱且不显著（r = ?0.079，p = 0.404），表明GWSA代理指标可靠性存在空间异质性。

2000–2025年，五种生态系统服务指标均呈现显著空间异质性，高值主要分布于山区、山前绿洲与沿河地带，低值集中于沙漠内部与荒漠边缘。时间上，CS呈最稳定上升趋势，HQ总体增加但2015年后阶段性下降、2021年后恢复，NPP波动明显但2010年后相对较高，ESI总体增加，而SR年际波动最强且无持续单调增加。

4.2 生态系统服务与GWSA代理指标的统计关系

2003–2025年同期年相关分析显示，GWSA与各生态系统服务的相关关系在空间上呈正负斑块交错格局，且服务间存在差异。GWSA与SR相关性总体较弱且空间破碎，北部局地负相关；HQ与CS模式相似，负相关集中于西部、北部及绿洲边缘；GWSA与NPP正相关覆盖范围较广，尤其北部与东部；ESI关联总体北部为负、南部与东部为正。联合趋势分区显示，ESI增加而GWSA下降的Type II区域占42.5%，两者同增的Type I区域占26.3%，同降的Type III区域占19.7%，GWSA增加而ESI下降的Type IV区域占11.5%。上述联合趋势仅表明空间共现，不证明生态系统服务变化导致了GWSA变化。

4.3 预测因子重要性与非线性统计关联

五个CatBoost模型测试集R²分别为ESI 0.973、HQ 0.977、CS 0.976、NPP 0.868、SR 0.870。算法敏感性检验显示，XGBoost对ESI、HQ、CS、SR、NPP的测试集R²分别为0.974、0.978、0.976、0.886、0.873，CatBoost分别为0.973、0.977、0.976、0.870、0.868，随机森林分别为0.973、0.976、0.974、0.836、0.868。CatBoost在五个目标上均呈现最小的训练–测试R²差距，过拟合程度较低，故保留用于SHAP解释。

SHAP依赖分析显示，距沙漠距离具有最大的SHAP依赖幅度，靠近沙漠位置与负向贡献相关，距离增加后转正并约在1 km以上趋于稳定；NDVI在0.13以下主要呈负贡献，以上转正；海拔呈倒U型关联，约1.4–4.3 km间为正贡献；潜在蒸散（PET）在约792、1200与1800 mm year^?1处出现转变；人类足迹指数在4.2以上贡献转负且随指数增大进一步下降；黏粒含量在19.3%以上由正转负；较陡坡度与负贡献相关；降水在中等值处为正贡献，约503 mm year^?1以上转负；距河流距离关联较弱；GWSA低于?18.6 mm时与显著负向ESI预测贡献相关，高于该候选转变点则贡献为正。

SHAP交互分析显示，残差依赖分数介于0.01–0.35之间。NDVI与黏粒含量交互最强（0.35），其次为海拔与沙粒含量（0.30）。NDVI还与土壤有机碳、海拔、坡度及沙粒含量存在交互；海拔与人类足迹指数及黏粒含量交互；距沙漠距离与人类足迹指数及土壤质地交互。最强交互模式涉及植被–土壤质地与地形–人类活动。

五、讨论与结论

讨论部分指出，2000–2025年ESI总体增加而SR下降，表明区域绿化并未同步增强所有生态系统服务。防风固沙能力更依赖于植被位置、群落结构与防护体系连续性，而非单纯的植被覆盖增加。GWSA与生态系统服务未同步变化，ESI增加而GWSA下降区域的存在与生态效益–水资源消耗权衡的可能性一致，但不证明修复导致了地下水损失。GWSA是大尺度残差储量代理指标，所识别的?18.6 mm转变点是模型依赖的候选阈值，并非直接的地下水位阈值，管理上应将其作为筛查信号，结合井观测、用水记录与局地水文评估后再调整修复强度。

研究局限性包括：GWSA为GRACE/GLDAS派生代理指标而非地下水位的直接测量，冰川、湖泊、水库等残差信号及GLDAS分量误差可能贡献于估计信号；仅阿克苏与策勒两站可用于独立对比，稀疏站点网络与较短重叠记录无法验证全部九个流域；GRACE有效空间支持远粗于局地灌溉与地下水开采过程，GDP、夜间灯光与人类足迹指数代表广义人类压力，未直接量化开采量、灌溉强度或用水效率；等权集成方案可能影响ESI量值；转变点与交互模式仅从CatBoost解释，可能因算法与空间划分而异；生态系统服务为年尺度数据，无法解析季节尺度地下水–植被滞后关系。

研究结论可归纳为：（1）2003–2025年GWSA代理指标线性下降约?0.31 mm month^?1，2020年后负异常加剧，下降集中于西北部；2000–2025年CS、HQ、NPP与ESI总体增加，SR波动并总体下降。（2）2003–2025年生态系统服务与GWSA代理指标呈现显著空间错配，ESI增加而GWSA下降区域占42.5%，同增、同降及GWSA增加–ESI下降区域分别占26.3%、19.7%与11.5%，这些数值表明空间共现而非生态改善对地下水储量的因果效应。（3）五个生态系统服务模型测试集R²介于0.868–0.977；ESI的SHAP依赖分析识别出GWSA约?18.6 mm、NDVI约0.13、人类足迹指数约4.2的模型特定候选转变点，这些是条件统计关联而非因果或普适生态阈值。总体而言，干旱区生态管理应协调生态效益与地下水安全，在区域尺度应用上述结论，并在制定局地管理决策前以更多地下水井、开采记录与过程水文分析加以验证。