利用混合深度学习方法，将青藏高原地表水储量变化归因于气候因素和人为因素

《Journal of Hydrology》：Attributing terrestrial water storage changes on the Tibetan Plateau to climate and human drivers using a hybrid deep learning approach

　　杨志强|杨爱丽|周雄
厦门理工大学环境科学与工程学院，中国福建省厦门市

**摘要**
青藏高原（TP），被誉为亚洲的“水塔”，对气候变化和人为活动都极为敏感。然而，评估陆地水储量（TWS）的时空动态及其在TP上的变化仍然具有挑战性。这主要是由于卫星观测的不连续性以及气候驱动因素与潜在人为影响之间的复杂相互作用。为了解决这些限制，我们开发了一个基于变分模态分解（VMD）-卷积神经网络（CNN）-长短期记忆（LSTM）的TWS预测和层次归因（VCLHA）框架。VCLHA能够量化气候因素（降水、蒸散、径流和温度）以及可能与人类活动相关的残差成分对TWS变化的独立贡献。评估结果表明，该框架性能稳健，纳什-萨特克利夫效率（NSE）值为0.983。通过与独立的重力恢复与气候实验（GRACE）类似重建方法（如GRACE-REC、CSRM-Bayesian和BNML_TWSA）在GRACE/GRACE-FO观测间隙期（2017年7月至2018年5月）的比较进一步证明了其稳健性，显示出一致的高原平均变化趋势。预测结果显示出明显的空间异质性：北部TP的水量增加了1.87毫米/年，而南部则减少了8.17毫米/年。降水被确定为主要的气候驱动因素，解释了TWS变化的43.2%。尽管与人类活动相关的残差成分没有明显的长期趋势（0.05毫米/年），但它表现出强烈的季节性特征，春季为正效应，夏季和秋季为负效应（分别为-2.81毫米和-1.62毫米）。这种季节性变化加剧了夏季至秋季（7月至10月）的水资源短缺，尤其是在南部和流出盆地。这些发现为气候变异性和人为影响在区域水资源变异性中的相对作用提供了定量证据，从而支持在TP上制定适应性和特定区域的水资源管理策略。

**通俗语言总结**
青藏高原为亚洲许多主要河流提供淡水，支持着数亿下游人口的生活。然而，确定近期水资源变化与气候变化和人为影响的关联程度颇具挑战性。在这项研究中，我们开发了一个混合深度学习框架，用于重建卫星重力观测数据，并评估自然过程和可能与人类活动相关的残差成分在水资源变化中的相对作用。2003年至2023年间，与气候相关的过程导致水资源储量持续长期下降（每年-0.10毫米）。虽然与人类活动相关的残差成分没有明显的长期趋势，但它具有显著的季节性影响，在春季有助于水资源增加，但在夏季和秋季加剧了水资源损失（分别为-2.81毫米和-1.62毫米）。这些季节性短缺在南部下游盆地最为严重，而北部和内陆地区在春季有所恢复。降水被认为是影响水资源储存的最重要气候因素。这些见解有助于提高干旱应对能力、灌溉规划和水库管理，因为气候变化正在重塑区域水循环。

**1. 引言**
陆地水储量（TWS）是全球水循环和陆气相互作用过程的关键组成部分（Rodell等人，2018年）。青藏高原以其丰富的陆地水资源而著称（Immerzeel等人，2020年；Yao等人，2012a年），在亚洲水文系统中发挥着关键作用（Meng等人，2019年）。最近的研究表明，TP上的TWS存在显著变化（Li等人，2022年）。越来越多的证据表明，由于全球变暖的持续影响，TP上的TWS正在持续减少，从而危及这一亚洲“水塔”的长期水资源可持续性（Liu等人，2024年）。全球多模型分析进一步表明，持续变暖将加速TWS的枯竭并加剧干旱的严重程度，强调了区分气候驱动效应和人为影响的必要性（Zhou等人，2022年）。这种枯竭可能对下游地区造成水资源短缺的风险（Pritchard，2019年），并可能影响数亿人（Zhang等人，2023年）。此外，由于气候变化的加速，TP上TWS的稳定性面临更多挑战（Yao等人，2019年）。研究TWS的动态并确定其变化机制至关重要。

**2. VMD-CNN-LSTM TWS预测和层次归因（VCLHA）框架的开发**
本研究开发了VCLHA（基于VMD-CNN-LSTM的TWS预测和层次归因）框架，用于重建、分析和归因TP上陆地水储量异常（TWSAs）的动态。该框架整合了先进的信号处理、深度学习和统计方法，以解决数据不连续性问题，表征时空变异性，并量化驱动因素。VCLHA框架包括四个相互连接的模块，以协同方式全面理解TWSA动态：
（i）模块1（GRACE数据填补）利用VMD分解TWSA，然后通过CNN-LSTM进行填补，从而获得连续数据。
（ii）模块2（时空变化特征分析）使用Mann-Kendall、Sen斜率、STL和DWT分析时空变化。
（iii）模块3（量化气候变化对TWS的影响）通过偏相关性和层次划分评估气候对TWS的影响。
（iv）模块4（量化人为活动对TWS的影响）将GRACE观测数据与Noah陆地表面模型的仅气候模拟数据进行比较，以推断可能与人类活动相关的残差成分。

**2.1 GRACE数据填补**
为了减少单个GRACE（-FO）重力异常解的不确定性，我们采用了三角帽（TCH）不确定性估计方法，随后采用贝叶斯三角帽（BTCH）融合方案（Zhang等人，2024年）融合CSR、JPL和GSFC的陆地水储量异常（TWSA）产品。该方法通过比较三个产品的差异来估计每个产品的相对误差水平，无需外部真实值。

**2.2 VCLHA框架的实现**
VCLHA框架利用变分模态分解（VMD）处理复杂信号，卷积神经网络（CNN）捕捉空间模式，长短期记忆（LSTM）提取时间趋势，确保TWSA重建的准确性。具体目标如下：
（i）使用开发的VCLHA框架重建2003年至2023年的连续TWSA记录；
（ii）描述TP上TWSA变化的年际、季节性和月度模式；
（iii）通过层次划分量化降水、蒸散、径流和温度对TWS变化的独立贡献；
（iv）通过将GRACE观测数据与Noah陆地表面模型的仅气候模拟数据进行比较，评估自然因素和可能与人类活动相关的残差成分对TWS的相对影响。通过结合先进的机器学习技术和归因分析，本研究提出了一个全面的框架，用于监测和理解TP上的TWS动态。预期这些结果将有助于更好地理解这一气候敏感区域水资源变化的机制，并促进亚洲水资源安全的适应性和可持续管理策略的发展。张等人（2019年）比较了GRACE在不同时间点的不同解决方案，并确认了mascon方法比球谐解更能有效检测到TWSA的趋势（ResearchGate）。李等人（2022年）报告了基于JPL的趋势与宾夕法尼亚州立大学对各个时间点的水量平衡估计之间具有高度一致性。这些研究共同支持将JPL的mascon方法与增益校正结合使用，以实现大规模水文目标。为了解决GRACE和GRACE-FO任务之间11个月的数据空白期（2017年7月至2018年5月），本研究采用了一种混合方法，将VMD与CNN-LSTM模型相结合。所提出的方法旨在重建2002年至2023年期间整个时间点的连续可靠的TWSA时间序列。该方法建立在分解-整合框架的基础上（李等人，2022年；刘等人，2025年），其中非平稳信号首先被分解为更平稳的成分，以提高模型性能并减少预测不确定性。填补空白的过程分为三个主要阶段，如下所述：

（i）使用VMD进行数据预处理：变分模态分解（VMD）被用来将GRACE总水储量异常（TWSA）时间序列分解为一组有限的带限模态（Dragomiretskiy和Zosso，2014年）。通过将观测信号分解为特定模态的窄带成分，VMD有助于减轻非平稳性，并隔离与不同频率范围（即特征时间尺度）相关的变异性，从而为后续基于学习的重建和预测提供更易处理的表示。进行了一组敏感性实验来评估重建保真度与计算负担之间的权衡。根据这些测试，模态数量被固定为K=6，二次惩罚参数被设置为α=2000，并对所有网格点保持不变，以确保空间上一致的分解设置。图S2展示了TP的一个代表性VMD分解。提取的模态（IMF1-IMF6）按频率升序排列。在这个例子中，IMF1主要捕捉低频的多年变化（Phillips等人，2012年；Guo等人，2021年）。IMF2-IMF3对应于主要的季节性成分（年度和半年周期）。相比之下，IMF4-IMF6描述了更高频率的变化：IMF4反映了季节内的变化，IMF5代表高频波动，IMF6对应于残差变化。在随后的重建/预测阶段，每个模态都独立建模，最终TWSA估计是通过将预测的模态相加得到的，遵循标准的基于分解的、逐个成分的建模范式（例如，He等人，2019年；Zuo等人，2020年）。这种策略减少了成分之间的光谱干扰，并使学习模型能够捕捉特定模态的动态，这对于具有明显多尺度变异性的TWSA信号特别有益。

（ii）模型训练和预测：从VMD得到的固有模态函数（IMFs）作为独立的“子序列”。这种分离允许模型隔离特定频率的水文变异性，并在学习过程中防止尺度混合（Dragomiretskiy和Zosso，2014年；Zuo等人，2020年）。为每个单独的子序列训练一个专用的CNN-LSTM模型，使用12个月的滑动窗口。输入向量包括之前的12个月值（Xt-12,…,Xt-1），而目标定义为下一个月的值（Yt）。网络架构从一个序列输入层开始，接着是一个一维卷积层（32个滤波器，核大小=3），用于提取局部特征（LeCun等人，1998年）。在LSTM层之前应用ReLU激活函数（Krizhevsky等人，2017年），该层有100个隐藏单元，用于捕捉时间依赖性（Hochreiter和Schmidhuber，1997年）。最后，一个全连接层生成回归输出。历史GRACE观测数据覆盖了2003年至2017年6月，作为训练数据集。为了填补2017年7月至2018年5月的11个月空白期，本研究采用了一种递归预测策略（Taieb等人，2012年）。在这种方法中，预测值被迭代地反馈到模型中作为后续步骤的新输入。具体来说，第一个缺失的月份（2017年7月）是使用2016年7月至2017年6月的观测数据完全组成的窗口进行预测的。然后滑动窗口前进以预测2017年8月。在这个阶段，预测的7月值作为输入来替换缺失的观测值。这个递归循环一直持续到整个空白期被填满。补充图S3可视化了这一过程，其中蓝色立方体代表观测输入，橙色立方体代表预测输入。最后，通过汇总所有子系列的预测值来重建完整的连续TWSA序列。

（iii）训练策略和计算环境：模型参数使用Adam优化器（Kingma和Ba，2014年）进行优化，初始学习率为1×10^-4。采用了分段学习率计划，每100个周期将学习率减少0.5倍。最大周期数为200，小批量大小为16。训练/验证分割遵循时间顺序策略（80%/20%），以避免时间序列学习中的信息泄露。为了减少过拟合，应用了提前停止策略，如果验证损失连续10次评估没有改善，则终止训练。

（iv）性能评估：使用三个标准指标来评估模型性能：Nash-Sutcliffe效率（NSE）、相关系数（CC）和均方根误差（RMSE）。

（v）不确定性评估：使用深度集成框架（Lakshminarayanan等人，2017年；Klotz等人，2022年）量化VMD-CNN-LSTM重建中的认知不确定性。设y?m(x)表示第m个独立初始化模型的预测值，其中m=1，…，M。集成平均值被采用为最终重建的TWSA：y?(x)=1/M∑m=1My?m(x)，集成方差（标准差）被用作原始认知不确定性估计：σraw(x)=1/M-1∑m=1M(y?m(x)-y?(x))2。为了提高可靠性，使用在校准子集上估计的后验全局缩放因子c进一步校准原始不确定性。因此，95%预测区间（PI）定义为：PI95(x)=[y?(x)-1.96cσraw(x),y?(x)+1.96cσraw(x)]，其中PI95(x)表示位置x处的校准后95%预测区间，y?(x)是方程（9）中定义的集成平均值，σraw(x)是方程（10）中定义的集成方差，c是在校准子集上估计的后验全局缩放因子。因此，cσraw(x)表示用于构建预测区间的校准不确定性。

2.2. 空间时间变化特征分析
TWSA的趋势分析使用了两种广泛采用的非参数方法：Theil-Sen斜率估计器和Mann-Kendall（MK）检验。这两种方法由于其对抗异常值和非正态分布的鲁棒性，在环境时间序列分析中得到了广泛应用（Chen等人，2022年）。统计显著性在95%置信水平（p<0.05）下确定。TWSA动态使用Loess（STL，Cleveland等人，1990年）和离散小波变换（DWT，Mallat，1989年）进行季节-趋势分解。STL将月度TWSA序列分解为趋势、季节性和残差成分，从而将长期变化与周期性波动分离。相比之下，DWT使用Daubechies小波函数（db4）在五个分解级别上实现。这产生了细节系数（D1-D5）来捕捉高频扰动，以及近似系数（A5）来表示低频趋势（Don和Johnstone，1994年）。然后通过细节系数和近似系数之间的能量比来量化高频扰动和低频趋势的相对贡献（Kumar和Foufoula-Georgiou，1997年；Labat等人，2005年）。

2.3. 气候变化对TWSA影响的量化
为了评估气候因素的影响，首先根据方程（12）（Nie等人，2016年）从基于GRACE的TWSA序列中推导出2003-2023年的月度和年度TWSC：TWSC≈ds/dt≈dTWSAdt≈TWSAt+1-TWSAt-1/Δt。应用偏相关分析来量化四个气候变量（包括降水（P）、蒸散（ET）、径流（R）和温度（T）的独立效应。计算Durbin-Watson统计量（Durbin和Watson，1950年）来检测P、ET、R、T和TWSC序列中的时间自相关性（Anderson，1942年）。然后计算偏相关系数，并使用t检验在95%置信水平（p<0.05）下评估其显著性，以提供每个驱动因素对TWSC的直接影响的稳健度量（Leguina，2015年）。为了进一步量化每个变量的相对贡献，进行了层次分割（HP）分析（Chevan和Sutherland，1991年）。HP是一种用于解决预测变量之间多重共线性的强大工具，通过确定每个解释变量对响应变量的独立贡献并将其与联合贡献分开（Mac Nally，2002年）。这种方法已成功应用于将陆地水储量变化归因于各种河流流域的气候驱动因素（Zhang等人，2022年）。HP的基本理论将所有N个解释变量解释的总方差分为独立成分和联合成分。这种关系数学表达为（Chevan和Sutherland，1991年；Zhang等人，2022年）：∑c=1/NRc=∑c=1NIc+∑c=1NJc，其中N是解释变量的数量（这里N=4）。Rc表示第c个变量在单变量回归中的拟合优度（决定系数）。Ic和Jc分别表示第c个变量的独立贡献和联合贡献。为了明确量化联合贡献（Jc），计算它作为变量在单变量回归中解释的总方差（Rc）与其独立贡献（Ic）之间的差异：Jc=Rc-Ic。这种计算确保独立贡献和联合效应之和等于该变量的总解释力。

2.4. 人为因素对TWSA影响的量化
为了评估气候变异性和可能与人类活动相关的残差成分对陆地水储量变化（TWSC）的相对作用，采用了基于模型的比较方法。使用Noah陆地表面模型来模拟一个基线情景，在该情景中只考虑气候驱动因素，即降水（P）、蒸散（ET）和径流（R），从而分离出气候驱动的TWSC成分（Hua等人，2024年）。观测到的TWSC与模拟的气候驱动成分之间的残差计算如下：ΔTWSC=P-ET-R。ΔTWSH=ΔTWS-ΔTWSC，其中ΔTWSH是与人类活动可能相关的残差成分；ΔTWSC是气候驱动的陆地水储量变化；ΔTWS是GRACE观测到的陆地水储量的总变化。应该注意的是，ΔTWSH也可能反映与水预算不闭合、强迫误差、模型结构限制以及Noah和GRACE之间存储表示差异相关的不确定性。因此，ΔTWSH不应被解释为人类引起的水储量变化的直接或唯一衡量。

3. 研究区域和数据收集
青藏高原（TP），通常被称为“亚洲水塔”，是区域和全球气候系统的关键组成部分（Wu等人，2021年）。TP覆盖面积约为250万平方公里，位于北纬26°00′12″-39°46′50″和东经73°18′52″-104°46′59″之间，平均海拔超过4000米（Hu等人，2023年）。其地形特征是高山山脉、广阔的盆地、深谷和广泛的冰川系统，包括喜马拉雅山脉和昆仑山脉（Yao等人，2012a）。作为许多亚洲主要河流的源头地区，包括印度河、恒河、布拉马普特拉河、湄公河、萨尔温河、黄河和长江，TP是一个关键的淡水储备库（Huang等人，2023年）。其水储量，包括冰川、积雪和地下水，对于维持区域水文循环和确保下游水资源至关重要（Immerzeel等人，2020年）。从气候上看，TP表现出明显的空间异质性，这是由其复杂的地形和陡峭的海拔梯度驱动的。主要的东-西山脉将TP划分为几个不同的气候带。北部地区具有寒冷干旱的气候，年平均温度在-7.6°C至0.1°C之间，年降水量在229毫米至594毫米之间（Huang等人，2023年；Yao等人，2012年）。南部地区受温带至亚热带季风气候的影响，年平均温度较高（-2.1至8.1°C），年降水量较大（287毫米至1,346毫米）（Wang等人，2020年）。这种气候变异性使TP对环境变化非常敏感，并突显了其在维持亚洲水文过程中的重要性（Immerzeel等人，2020年；Yao等人，2022年）。为了监测TP上的TWSA，本研究利用了重力恢复和气候实验（GRACE）及其后续任务（GRACE-FO）的数据。本研究没有简单地平均数据，而是采用了贝叶斯三角帽（BTCH）方法（Zhang等人，2024年）来融合CSR、JPL和GSFC的TWSA产品。这种结合的mascon解决方案旨在最小化不确定性并提高数据的整体可靠性。这种集成方法利用了每个中心处理方法的独特优势，以增强数据的稳健性（Loomis等人，2020年；Scanlon等人，2018年）。分析得出的空间平均权重分别为：JPL约为0.472，GSFC约为0.318，CSR约为0.211，空间分布详见补充图S1（Zhang等人，2024年）。由此产生的数据集提供了0.5° × 0.5°网格上的月度TWSA估计值。需要注意的是，该网格代表的是输出分辨率，而非GRACE任务的固有空间分辨率，由于卫星重力测量的特性，其固有空间分辨率约为300公里（Tapley等人，2019年；Sakumura等人，2014年）。为了支持仅考虑气候因素的水分储存动态模拟并分离人为信号，采用了全球陆地数据同化系统（GLDAS）。GLDAS整合了卫星和地面观测数据，使用四种陆地表面模型来模拟陆地表面状态和通量，包括Noah陆地表面模型、流域陆地表面模型（CLSM）、社区陆地模型（CLM）和可变渗透能力（VIC）模型（Rodell等人，2004年）。在这四种模型中，本研究选择了GLDAS-Noah模型，因为它在之前的研究中表现出较高的性能，并且与GRACE推导出的TWSA有很强的相关性（r > 0.85，Scanlon等人，2018年）。具体来说，本研究中使用了GLDAS-Noah L4月度0.25° × 0.25° V2.1产品（GLDAS_NOAH025_M_2.1），数据覆盖时间为2003年至2023年。

模型性能评估使用了三个标准统计指标（NSE、CC和RMSE）。提出的VMD-CNN-LSTM模型与VMD-LSTM模型（不含CNN组件）和CNN-LSTM模型（不含VMD预处理）进行了比较，这种比较有助于评估空间特征提取和信号分解对高原地区TWSA预测准确性的单独贡献。图3a-c展示了所有网格单元的NSE、CC和RMSE的箱形图。VMD-CNN-LSTM模型展示了最高的中位数NSE（0.983）和CC（0.992），以及最低的中位数RMSE（7.12毫米），并且四分位数范围最窄，表明其在不同水文气候区域具有很强的空间一致性和鲁棒性。VMD-LSTM模型也取得了有竞争力的性能（中位数NSE = 0.967；CC = 0.984；RMSE = 9.77毫米），而CNN-LSTM模型的结果明显较差（NSE = 0.876；CC = 0.943；RMSE = 37.49毫米），且指标分散度较大，显示出在高原地区性能的显著差异。

为了量化VMD-CNN-LSTM重建中的认知不确定性，采用了深度集成框架（Lakshminarayanan等人，2017年），该框架在水文深度学习中得到广泛应用，通常比单一确定性模型提供更可靠的不确定性估计（Klotz等人，2022年）。训练了一组独立初始化的网络来表示随机优化引起的变化。集成平均值被用作重建的TWSA，而集成方差（成员间的标准差）被用作原始的认知不确定性估计。通过将RMSE图（图4a）与不确定性图（σ；图4b）进行比较来评估空间可靠性，观察到强烈的空间一致性（r = 0.871），表明重建误差较大的区域通常伴随着较高的不确定性。原始不确定性进一步通过在后验全球缩放方法上进行校准，并在独立测试集上进行评估。校准后，95%预测区间实现了可靠的样本外覆盖（PICP = 0.955；图5b）。此外，校准后的不确定性与测试集上的绝对误差呈正相关（Spearman ρ = 0.60；图5c），支持其作为预测置信度的指标。图5a展示了代表性的时间序列示例。

为了评估GRACE/GRACE-FO观测间隙期间（2017年7月至2018年5月）的重建可靠性，将VMD-CNN-LSTM产品与三个公开可用的类似GRACE的数据集进行了比较，这些数据集代表了不同的重建范式：GRACE-REC（Humphrey和Gudmundsson，2019年）、CSRM-Bayesian（Rateb等人，2022年）和BNML_TWSA（Mandal等人，2025年）。为了保持一致性，所有产品都被重新网格化为0.5°的分辨率，并使用相同的TP掩码和异常定义进行评估。图6a表明，重建的高原平均TWSA在整个记录期间再现了参考产品的季节性变化，并在间隙期间保持了连贯的季节性减少-恢复行为。图6b进一步显示，差异在零附近波动，没有持续的偏差，在间隙期间，重建的TWSA与参考产品高度一致，相关系数为r = 0.827–0.939。间隙期间的空间一致性在图7a–c中总结：内部高原的RMSE较低，而冰川主导边缘的RMSE较高，表明差异主要发生在具有强烈冰冻圈信号的区域内。

图8展示了2003–2023年期间高原TWSA趋势的时空分布。空间分析揭示了水分储存动态的显著异质性，表现为明显的南北梯度。北部地区（北纬33°以北）主要呈现正趋势，年平均增加量为1.87毫米/年。大约69.8%的网格显示出水分储存增加，增幅最大的区域位于中北部和东北部。相反，南部地区（北纬33°以南）普遍出现水分减少，平均趋势为-8.17毫米/年。在南部区域的84.1%范围内观察到相应的降水负趋势，特别是在布拉马普特拉山脉和西南部地区。图9显示了2003至2023年期间GRACE/GRACE-FO间隙期间现有类似GRACE的重建产品的空间分布。季节性TWSA的合成图揭示了年内的显著变化和空间不对称性。平均而言，高原在夏季增加水分（+5.09毫米），在冬季减少水分（-36.59毫米），春季（-16.19毫米）和秋季（-13.75毫米）作为过渡阶段。空间上，春季的亏缺在中部和南部地区最为明显，而西北部地区仅有轻微的夏季增加。秋季标志着广泛减少的开始，特别是在西部、中部和南部地区达到最低点。干旱的西北部地区季节性变化最小，而东南部地区的储存变化最为显著。

为了描述2003至2023年期间TWSA的时间动态，时间序列被分解为长期趋势、季节性周期和残差分量。季节性周期的平均幅度约为58.2毫米，8月达到峰值，1月达到最低值。长期趋势显示出统计上显著的下降趋势，每年减少2.77毫米（R2 = 0.674，p < 0.001），到2023年底减少到约-57.5毫米。残差遵循近似正态分布（μ = 1.23毫米，σ = 11.25毫米），证实STL方法能够有效分离确定性信号。图10展示了使用STL方法对2003–2023年期间高原每月TWSA的分解。图11使用五层DWT进一步分析了频域特征。大约70.5%的总信号能量集中在长期近似分量（A5，对应于超过64个月的周期），证实了低频变异性占主导地位。细节分量（D1-D5）代表了不同时间尺度上的变化，其中中间周期（8–16个月）对应于季节性和年际变化，而高频分量（2–8个月）表示短期、亚年波动。这些结果表明，高原上的TWSA变化受季节性周期、年际变化和持续的长期下降趋势的共同影响。图14显示了2003年至2023年间TWSC与关键气候因素（降水量P、蒸散量ET、径流量R和温度T）之间的偏相关系数的空间分布模式。每个相关性都控制了其他影响因素。降水量（图14a）在78.2%的网格单元中显示出显著的正相关性（p < 0.05），这些区域的平均相关性为0.298。强正相关性主要集中在东部和南部地区，强调了降水量在增加水资源储存量中的关键作用。蒸散量（图14b）在16.8%的网格单元中显示出显著负相关性（p < 0.05），这些区域的平均相关性为-0.177。这种负相关性在西部地区尤为明显，表明蒸散量在区域水资源枯竭中起着重要作用。径流量（图14c）在31.54%的网格单元中显示出显著负相关性（p < 0.05），平均相关性为-1.85。这种模式主要集中在东南部边缘，表明径流量的增加可能会由于水分输出或渗透减少而减少水资源储存量。温度（图14d）的相关性较弱且空间变化较大，在20.19%的网格单元中显示出显著相关性（p < 0.05），平均相关性为-0.333。这种异质的空间变异性表明温度通过复杂机制影响TWSC，例如加速雪/冰川融化或增加蒸散量。

图14. 在控制其他三个变量的情况下，TWSC与降水量（a）、蒸散量（b）、径流量（c）和温度（d）之间的偏相关性。十字标记表示p < 0.05的区域。

层次分割分析用于确定2003年至2023年间每个气候变量对TWSC变异性的相对贡献。该方法将TWSC的方差分解为每个因素的独立贡献，并考虑了它们之间的相互依赖性。分析表明降水量是主要驱动因素，解释了TWSC变异性的43.2%（图15a）。其影响在南部和中部地区尤为显著（>60%），这两个地区占研究区域的13.6%。温度是第二重要的因素，平均解释了24.9%的方差（图15d）。其贡献在西北部地区最为明显，在4%的网格单元中超过70%。这可能反映了温度驱动过程（如雪和冰川融化）在干旱的高海拔环境中的重要性。蒸散量平均贡献了17.4%的TWSC变异性，在中部地区影响稍大（图15b）。径流量是影响最小的因素，仅解释了14.5%的方差，在整个区域有轻微的空间变化（图15c）。

为了进一步研究区域差异，评估了12个主要河流流域中气候驱动因素的相对贡献（图16）。结果证实，在大多数流域中，降水量是最主要的因素，特别是在受季风影响的湿润的东部和南部流域，如布拉马普特拉河、恒河、长江和萨尔温江。降水量也在中部和北部流域（如河西内陆河流域、柴达木河、黄河和内陆流域）中起着重要作用。相比之下，在较干燥的西部流域（如阿姆河、印度河和塔里木河），温度成为主要因素，由于其对融雪和蒸散量的控制，贡献了41.9-52.2%的TWSC变性。降水量在东部和中部流域的贡献较低，但在长江和黄河等地区仍占10-25%。

图17和表1展示了12个主要流域中气候驱动因素对TWSC变异性的季节性贡献。在受季风影响的湿润流域（如恒河、布拉马普特拉河和湄公河），降水量仍然是主要的季节性驱动因素，占方差的52.9-67.8%。相比之下，阿姆河、印度河和塔里木河等干旱的高海拔流域的季节性TWSC幅度有限（<50毫米），这突显了内陆排水系统的水文刚性。在这些地区，温度成为主要驱动因素，约占方差的44.0-56.4%。同时，在其他干旱地区，如河西和柴达木河，降水量和蒸散量共同调节水资源的季节性变化。具体来说，降水量贡献了约50%的方差，而蒸散量约占23%。

为了评估气候变异性和可能与人类活动相关的残余成分的相对作用，我们将基于GRACE的TWSC数据与Noah地表模型的仅气候模拟进行了比较。图18a显示了2003年至2023年间气候驱动（ΔTWSC）成分和可能与人类活动相关的残余成分（ΔTWSH）的长期趋势。气候相关贡献显示出水储存量的一致下降趋势，每年下降0.10毫米（p < 0.01）。相比之下，可能与人类活动相关的残余成分表现出显著的年际变化，从2017年的最大值+5.0毫米到2004年的最小值-2.64毫米。季节性分解（图18b和表2）显示，气候因素对水储存量的影响最大，夏季平均增加+7.40毫米，秋季平均减少-12.12毫米。与此一致的是，可能与人类活动相关的残余成分在春季有正面贡献（+2.09毫米），但在夏季和秋季有负面影响（-2.81毫米和-1.62毫米）。

图18c和表3揭示了2003年至2023年间12个流域中可能与人类活动相关的残余变异性（ΔTWSH）的显著空间差异。年际幅度在南部和流出流域最大（阿姆河60.07毫米；恒河58.46毫米；印度河45.20毫米；布拉马普特拉河38.14毫米），在东南部流域适中（湄公河39.37毫米；萨尔温江34.86毫米），而在北部和内陆-东部流域较弱，长江为39.49毫米，黄河为31.02毫米，柴达木河为16.64毫米，塔里木河为19.61毫米，河西为8.48毫米。月度分析（图18d）表明北部和内陆流域在春季（3月至5月）有增加，而在夏季和秋季（7月至10月）有持续损失。

重建的基于GRACE的TWS数据显示，2003年至2023年间高原地区TWSC长期下降了2.77毫米/年。这一趋势与之前的基于GRACE的评估结果一致（Li等人，2022年；Wang等人，2020年），这些评估也发现整个流域大部分地区存在持续的水资源损失。空间异质性表现为北部略有增加（1.87毫米/年），南部显著减少（-8.17毫米/年），这与Li等人（2022年）的发现一致。Immerzeel等人（2020年）记录了这些地区对气候强迫、冰川动态和地形控制的响应差异。本研究中的定量归因分析改进了Wang等人（2020年）之前的基于相关性的方法。本研究应用层次分割来区分气候驱动因素的独立效应。降水量成为年际TWSC变异性的主要控制因素，占方差的43.2%，而温度贡献了24.9%。这些结果扩展了Deng等人（2022年）的发现，表明温度驱动的蒸发需求在整个流域具有广泛的影响。值得注意的是，本研究还明确量化了可能与人类活动相关的残余成分，这是以往评估中较少考虑的方面（Li等人，2022年；Yao等人，2019年）。结果表明，尽管与气候强迫相比长期趋势较弱，但可能与人类活动相关的残余成分在整个流域表现出明显的季节性变化。该残余成分没有明显的十年趋势（0.05毫米/年），但在春季有正面效应，在夏季和秋季有负面效应（-2.81毫米和-1.62毫米）。这种季节性行为加剧了夏季-秋季（7月至10月）的水资源短缺，尤其是在南部和流出流域，而北部和内陆流域通常在春季有所恢复。

主导驱动因素的空间变异性反映了特定地区的水文气候特征。在受季风影响的南部和东部流域（如恒河、布拉马普特拉河和湄公河），TWS变化主要由降水量控制，占方差的52.9%至67.8%（Wang等人，2024年）。相比之下，干旱的西部流域（如阿姆河、印度河和塔里木河）对温度非常敏感（41.9%至52.2%），这突显了雪和冰川融化等冰冻圈过程的作用（Uddin等人，2025年）。在北部地区，降水量增加（Yao等人，2019年）和融水渗透（Zou等人，2022年）似乎有助于水资源储存的增加。在南部流域，水资源枯竭主要由冰川退缩和在变暖条件下增加的蒸散量共同驱动（Li等人，2022年；Pritchard，2019年）。可能与人类活动相关的残余成分在夏季-秋季（7月至10月）的季节性加剧可能反映了自然水文供应与灌溉需求、水库调节和其他人类相关过程之间的时间不匹配，从而加剧了南部和流出流域的气候驱动下降。这种不匹配在干旱时期增加了区域水资源的不稳定性。

基于残差的方法将气候驱动成分和可能与人类活动相关的残余成分在TWSC中的分离受到几个不确定性的影响。首先，当降水量、蒸散量、径流量和储存量变化来自不同的观测或模型产品时，水预算不闭合是常见的（Pan等人，2012年；Abolafia-Rosenzweig等人，2021年；Lehmann等人，2022年）。先前的研究表明，多源分析中的此类残差误差可能仍然很大，而且表面上良好的闭合可能是由于预算成分之间的误差抵消，而不是数据质量的真正改善（Lehmann等人，2022年）。其次，残差项不应仅解释为人类活动的效果。现有的基于GRACE的研究指出，区分气候变异性和人类活动对水资源储存的影响仍然很困难，基于模型的估计可能与GRACE观测结果有很大差异（Liu等人，2021年）。此外，由于大多数地表模型不包括地表水和地下水，基于模型的储存估计与GRACE推导的储存变化之间的差异也可能反映储存表示的差异。因此，这里从GRACE观测和气候驱动的Noah模拟之间的差异推断出的残差项主要从其相对空间和季节性对比来解释。它不应被视为人类引起的水资源变化的直接或绝对估计。通过与独立的类似GRACE的重建/产品（如GRACE-REC和其他可用数据集）进行比较，进一步检验了其稳健性。高原范围内的长期趋势在各种产品之间总体上是一致的，而在年际时间尺度上以及受边界/冰冻圈影响的区域则存在较大差异。这种差异是预期之中的，因为不同的产品采用了不同的约束条件、正则化和平滑策略，这些策略会改变年际振幅和长期成分（Humphrey和Gudmundsson，2019；Sun等人，2021；Rateb等人，2022）。增益因子处理对于复杂地形下的流域尺度解释也非常重要。增益因子旨在减少GRACE质量异常解决方案中的衰减/泄漏效应，这在TP（青藏高原）地区尤为重要（Wiese等人，2016）。此外，使用BTCH策略进行多产品集成可以避免假设任何单一产品为真实情况，而是利用产品之间的相对一致性来获得具有不确定性信息的最佳估计值。这提高了高起伏地区流域尺度分析的稳健性。通过将GRACE卫星观测数据与基于深度学习的重建模型和归因分析相结合，本研究加深了对青藏高原水资源动态多方面驱动因素的理解。研究结果强调了青藏高原对季风强度和变化时间以及持续变暖导致的冰冻圈退化的脆弱性（Zhang等人，2023）。同时，随着持续变暖，冰川退缩、积雪减少和蒸发蒸腾作用加剧，预计将进一步加剧陆地水分损失。这种趋势威胁着青藏高原作为亚洲“水塔”的作用，并可能减少超过十亿下游居民的淡水供应（Immerzeel等人，2020）。在重建和归因方面仍存在不确定性。GRACE的有效分辨率和潜在的泄漏效应可能会影响流域尺度信号，尤其是在边界附近（Landerer和Swenson，2012；Wiese等人，2016）。归因还取决于仅基于气候的Noah模拟和气象强迫的质量，这些因素可能会将误差传播到非气候成分中。此外，分层划分量化了独立解释能力，但本身并不能确定因果关系。评估可能与人类活动（如地下水抽取和水库运营）相关的残差成分是一个重要的进展，尽管仍存在不确定性。需要通过高分辨率建模和地面监测进一步研究，以更精确地区分气候驱动因素和人为因素的影响。未来的研究应优先开发整合基于卫星观测、物理模型和社会水文数据的综合评估框架。这样的系统对于在高气候压力和人类活动压力下支持青藏高原的适应性流域水资源管理至关重要。

8. 结论

本研究开发了VCLHA（基于VMD-CNN-LSTM的分层归因）框架，用于重建和分析2003年至2023年的TWSA时间序列。该框架结合了变分模态分解（VMD）、卷积神经网络（CNN）和长短期记忆（LSTM）网络，以解决GRACE观测数据中的缺失和复杂动态问题。VCLHA采用分层划分和基于模型的比较方法，量化了气候驱动因素和可能与人类活动相关的残差成分对TWSC的相对贡献。模型稳定性得到了强有力的验证结果支持（中位数NSE=0.983；CC=0.992；RMSE=7.12毫米），并且不同网格单元之间的指标分散度较小，表明在异质水文气候条件下具有稳健的空间一致性。不确定性估计也具有信息性且校准合理（PICP=0.955），不确定性场与重建误差之间的空间一致性很高（r=0.871），支持了流域尺度解释的信心筛选。此外，在GRACE/GRACE-FO数据间隔期间（2017年7月至2018年5月），与独立的类似GRACE的重建结果相比，进一步证明了重建结果的稳健性。

主要发现总结如下：
(i) 观察到TWSA长期下降2.77毫米/年，存在显著的空间异质性。北部地区水资源有所增加（+1.87毫米/年，占网格单元的69.8%），而南部地区则出现广泛的水资源损失（-8.17毫米/年，占网格单元的84.1%）。
(ii) 降水是主要的气候驱动因素，解释了年际TWSC变化的43.2%（在青藏高原南部和中部部分地区超过60%）。温度总体贡献了24.9%，在西北部地区这一比例超过70%。蒸发蒸腾作用和径流分别贡献了17.4%和14.2%。
(iii) 与人类活动相关的残差成分加剧了季节性水资源短缺，尤其是在温暖季节。该成分在春季为正值，但在夏季（-2.81毫米）和秋季（-1.62毫米）转为负值，其中南部和流出流域的影响最为明显。

引用文献：
Donoho和Johnstone，1994；Gao等人，2023；Gneiting和Raftery，2007；Li等人，2022a；Li等人，2022b；Tavella和Premoli，1994；Taylor，2001；Velicogna，2009；Yan等人，2025；Yao等人，2012b；Zhang等人，2020。

数据可用性声明：
本研究使用的所有数据均来自开源仓库。JPL-M产品由物理海洋学分布式活动档案中心提供（https://doi.org/10.5067/TEMSC-3JC62）。CSR-M GRACE/GRACE-FO质量异常产品由德克萨斯大学空间研究中心提供（https://www2.csr.utexas.edu/grace）。GSFC RL06v2数据来自NASA戈达德地球科学部门项目（https://earth.gsfc.nasa.gov/geo/data/grace-mascons）。GLDAS-2.1 Noah陆地表面模型数据集来自NASA GES DISC（https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/）。
此外，用于外部基准测试的类似GRACE的重建数据集也是公开可用的：GRACE-REC发布在Figshare上（Humphrey & Gudmundsson，2019；https://doi.org/10.6084/m9.figshare.7670849）；贝叶斯CSR质量异常填充重建（CSRM-Bayesian）由德克萨斯数据仓库托管，并包含不确定性信息（Rateb等人，2022；https://doi.org/10.18738/T8/5MPOJU）；BNML_TWSA也发布在Figshare上，包括网格化重建结果及其相关不确定性（Mandal等人，2025；https://doi.org/10.6084/m9.figshare.2537669）。

CRediT作者贡献声明：
杨志强：撰写——原始草稿、可视化、验证、软件开发、方法论研究、数据管理。
杨艾丽：撰写——审稿与编辑、资金筹集。
周雄：撰写——审稿与编辑、监督、资金筹集、概念构思。