GRACE和GRACE-FO卫星的观测数据为监测陆地水储量变化提供了重要的数据支撑^[1]。但GRACE和GRACE-FO卫星间存在约1 a的数据空白，不利于构建更长时间连续的陆地水储量变化序列。因此，有学者提出采用低轨卫星的观测数据填补该数据空白^[2]，也有学者尝试将深度神经网络应用于填补陆地水储量变化的空缺^[3]。

长江是中国第一长河，利用GRACE和GRACE-FO模型分析长江流域的陆地水储量变化情况，可为保障社会经济的可持续发展提供一定的参考资料，具有重要意义。本文综合GRACE模型、GRACE-FO模型、GLDAS模型、实测降水量和实测气温等数据，分析基于不同激励函数的BP和RBF神经网络算法用于填补GRACE与GRACE-FO间空缺的精度及可行性。在优选填补方案后，利用连续长时间序列的时变模型分析长江流域陆地水储量的时空变化情况。

1 原理与方法 1.1 利用GRACE与GRACE-FO卫星反演区域陆地水储量变化

利用时变重力场模型计算区域平均陆地水储量变化(等效水柱高)的表达式如下^[4-5]：

$ \begin{aligned} \Delta h=& \frac{a \rho_{\mathrm{ave}}}{3 \rho} \sum\limits_{l=0}^{\infty} \sum\limits_{m=0}^{l} \bar{P}_{l m}(\cos \theta) \frac{2 l+1}{1+k_{l}} \omega_{l} \cdot \\ &\left(\Delta \bar{C}_{l m} \cos m \lambda+\Delta \bar{S}_{l m} \sin m \lambda\right) \end{aligned} $

(1)

式中，l、m分别为重力场模型的阶与次，a为地球赤道处的半径，θ、λ分别为地心余纬与地心经度，ρ_ave为地球平均密度，ω_l为高斯平滑系数，ΔC_lm与ΔS_lm为球谐系数变化量，k_l是l阶勒夫数，P_lm为归一化缔合勒让德多项式，ρ=1 000 kg/m³为水的密度。

1.2 GRACE与GRACE-FO间空缺数据的填补方法

选用BP和RBF神经网络用于填充GRACE与GRACE-FO卫星间空缺数据，其中BP神经网络为全局逼近前馈神经网络，RBF神经网络为局部逼近前馈神经网络。以GLDAS模型中的地表水变化(其中包括土壤水含量和雪水含量)、降水异常量(每月值扣除平均值)、气温异常量(每月值扣除平均值)为输入因子，陆地水储量变化为输出因子，构建BP和RBF神经网络填补模型。BP与RBF神经网络选用的隐含层激励函数见表 1。

另外，在采用BP和RBF神经网络处理数据前，为了消除不同输入因子与目标因子量纲等的影响，对变量进行归一化处理。

2 GRACE与GRACE-FO卫星间的数据空缺填补方案的确定 2.1 数据源

本文所采用的数据有GRACE/GRACE-FO月时变重力场模型、GLDAS水文模型、降水和气温等，具体见表 2。

2.2 GRACE与GRACE-FO月时变重力场模型预处理

由于GRACE与GRACE-FO月时变重力场模型存在高频噪声和南北条带误差等，需要对模型进行以下处理：1)采用Sun等^[6]计算的地心改正项填补GRACE/GRACE-FO模型的一阶项；2)以SLR测定的C₂₀项替换GRACE和GRACE-FO模型的C₂₀项^[7]；3)对时变模型进行冰川均衡改正(GIA)，以减少冰后回弹的影响^[8]；4)采用Swenson去条带方法^[9]减小模型的条带误差；5)利用半径为300 km的高斯滤波降低模型的高频噪声^[10]。

2.3 基于BP和RBF神经网络算法的数据处理方案

陆地水储量变化受到地表水、降水、气温变化等因素的综合影响，图 1给出了在长江流域GRACE模型反演的陆地水储量、GLDAS模型计算的地表水、降水异常量和气温异常量的变化情况。由图 1可知，在长江流域，各变量均呈现出明显的季节性变化，在地表水及降水较多的时期，陆地水储量也较高；在气温较高的时期，陆地水储量较低；GRACE陆地水储量与地表水、降水、气温的相关系数分别为0.83、0.66和0.26。本文选择4种不同输入因子方案进行精度分析，见表 3。

为了验证不同隐含层激励函数的BP和RBF神经网络用于填补GRACE与GRACE-FO间空缺数据的精度及可行性，取3个与空缺数据同时段(2017-07~2018-05)的GRACE数据进行模拟验证分析(表 4)。分别以模拟空缺时期前后各1 a、前后各2 a及除空缺外的所有GRACE数据作为输入的目标数据，经分析，BP神经网络取空缺前后各2 a数据、RBF神经网络取空缺前后各1 a数据作为输入的目标数据所得到的效果最佳。

为分析不同的神经网络隐含层激励函数与不同输入因子组合方案的精度，以获得填补GRACE与GRACE-FO之间空缺数据的最优方案，设计了基本流程(图 2)。

2.4 BP和RBF神经网络填补数据空缺的精度分析

以长江流域为例，根据4种输入因子方案，统计3个时期GRACE反演的陆地水储量变化与各神经网络拟合陆地水储量变化的相关系数(表 5)。由表 5可知，不同输入因子、不同时期各神经网络拟合与GRACE反演的陆地水储量变化均有较强的相关性，相关系数均大于0.85。进一步将各神经网络拟合所得3个时段的陆地水储量变化与相应的GRACE反演的陆地水储量变化作差，差值的标准差见表 6(单位cm)。由表 6可知，方案1、3和4的标准差较小，说明这3种方案拟合效果比较好。图 3给出方案4拟合结果与GRACE反演结果的对比情况。由图可见，各神经网络拟合的陆地水储量变化与GRACE反演的陆地水储量变化较为符合，二者振幅、相位相近，说明拟合效果较好。结合表 5、表 6数据可知，时段3的拟合效果较时段1、2差，这与GRACE任务结束前2~3 a数据质量较差有关。综合来看，方案4作为输入因子时，隐含层为ReLU的BP神经网络算法的结果较优，将作为后续GRACE和GRACE-FO空白数据的填充方案。另外，由于GRACE卫星传感器性能下降和供能不足等问题，2002-04~2017-06期间ITSG-Grace2018时变重力场模型共有21个月数据缺失，GRACE-FO卫星由于微波仪器出现异常导致2018-08、2018-09数据缺失，以上缺失数据也采用该神经网络进行填补。

利用上述最优的BP神经网络拟合的3个时段陆地水储量变化数据，分别替代该时段原有GRACE反演的陆地水储量变化数据，并求得替代前后2002-04~2017-06长江流域的陆地水储量变化率(表 7)。由表 7可见，替代GRACE反演的数据前后，陆地水储量变化率的差异较小，说明采用隐含层激励函数为ReLU的BP神经网络填补GRACE与GRACE-FO间的数据空缺后对长江流域陆地水储量变化进行连续长时间序列的时空分析具有可行性。

3 2002~2020年长江流域陆地水储量变化的时空分析 3.1 长江流域陆地水储量的空间变化

利用GRACE/GRACE-FO模型(GRACE与GRACE-FO间空缺利用隐含层激励函数为ReLU的BP神经网络进行填补)、GLDAS模型、降水和气温数据，分别计算长江流域陆地水、地表水、降水和气温变化趋势的空间分布，其中，将长江流域分为上、中、下游进行分析，结果见图 4。由图可见，从整体上看，近20 a来长江流域陆地水储量基本呈上升趋势，尤其是长江上游东部及中游陆地水储量明显增多，速率达到5 mm/a以上，但存在空间分配不均的情况。20世纪60年代以来，青藏高原呈现升温趋势^[11]，导致该地区冰川消融，如唐古拉山的冰川消融^[12]使长江上游中西部区域陆地水储量减少，速率约为-5 mm/a。大量的冰雪融水很好地补充了长江流域其他区域的陆地水储量，唐古拉山冰川消融对整个长江流域水储量贡献达到15.35%^[12]，可见长江发源地的冰川消融是长江部分区域水储量上升的重要原因之一。长江流域的上游东部、中游西部部分地区降水量增加，该区域陆地水储量也出现明显上升，说明降水也是这些区域陆地水储量上升的重要驱动因素之一。

3.2 长江流域陆地水储量的时间变化

图 5给出长江流域不同区域的陆地水、地表水、降水和气温异常量变化情况。由图 5可见，长江流域不同区域的陆地水、地表水、降水、气温均呈现明显的季节性变化，且4者的变化周期较为符合，具有一定相关性。经计算，2002-04~2020-06长江流域全流域、上游、中游、下游的陆地水储量变化速率分别为2.6 mm/a、1.8 mm/a、4.3 mm/a、0.9 mm/a。据中国水资源公报^[13]的统计数据显示，2002~2019年长江流域水资源总量以约62.19亿m³/a的速率增长，而GRACE/GRACE-FO卫星探测到该时期的长江流域陆地水储量亦处于增长趋势，其中长江中游增长趋势较大。

对陆地水储量变化与地表水储量变化的时间序列进行拟合^[14]，表 8给出2002~2020年长江流域不同区域陆地水储量变化与地表水储量变化的周年振幅与相位。由表可见，长江全流域和中游的陆地水储量变化的振幅大于地表水，这是因为，陆地水储量是该区域水文情况的综合体现，而地表水只包含了土壤水与雪水。长江上游和下游区域的地表水变化振幅大于陆地水储量变化，原因在于，上游区域受人为因素影响较小，使地下水等变化量较小，且气候变化使得该地区冰雪融化量较大；而长江下游存在南水北调工程，且三峡大坝的运行会使下游的陆地水与地表水减少^[15]，但较多的降水又补充了该地区的陆地水与地表水，综合各种水文变化过程，出现地表水变化振幅大于陆地水储量变化的情况。从长江全流域来看，陆地水储量与地表水储量变化的周年相位差只有约0.1°，两者的年际性变化符合；但在中、下游周年相位差大于30°，这是由于中、下游的地表水受降水影响较大，而区域陆地水储量变化的信号比降水变化信号滞后。

降水的时空分配不均是造成干旱、洪涝等自然灾害的因素之一，利用GRACE/GRACE-FO卫星可以监测干旱、洪涝灾害。经计算，在长江流域各个区域GRACE/GRACE-FO卫星观测的陆地水储量变化与降水变化的周年相位差约为30°~40°，这是因为，GRACE/GRACE-FO某月的观测量包含了该月之前的降水部分，导致陆地水储量变化相对于降雨变化有一定的滞后^[16]。长江流域水资源公报显示^[17]，长江流域在2006年降水比常年偏少10.3%，在6~9月出现严重旱情；在2016年降水比常年偏多10.9%，5~8月发生严重洪涝灾害；2017-06~2017-08降水明显偏多，使得长江中游发生大洪水，2017-09~2017-10三峡水库发生4次入库洪峰约30 000 m³/s的洪水，长江上游发生秋汛；2019年夏季长江中游陆地水储量增多，中下游在7月迎来洪水，但在同年秋、冬两季雨水偏少，长江中下游发生持续干旱。以上的洪涝、干旱灾害过程在GRACE/GRACE-FO卫星反演的长时间序列水文信号中均有体现(图 5)。长江上、中、下游的陆地水储量对于旱、涝灾害与降水变化的响应程度不同，上游对两灾与降水变化响应程度较小，年际变化相对平稳；中、下游对两灾与降水变化响应程度较大，年际变化相对剧烈。三峡大坝能够对蓄水量进行调控，且中、下游降水变化量较大，故造成了长江流域上、中、下游陆地水储量变化对两灾与降水异常变化响应程度的差别。长江流域上、中、下游对于降水异常量变化的响应不一致，但其陆地水储量变化与降水相关系数分别达到0.62、0.51、0.44，相关性达中度至强相关，可见降水仍是长江流域陆地水储量年际变化的重要影响因素之一。

气温可通过影响水的蒸发量、冰雪融化量等进而影响陆地水储量变化。由图 5可见，2006年夏季的高温天气对长江流域旱灾的发生有一定贡献，而在2011年冬季长江流域气温较低且9月长江上游北部降水异常偏多^[17]，导致2011年秋冬季的陆地水储量减少得较慢。在地势较高、冰川较多的长江上游区域，气温是其陆地水储量变化的主要因素，而其融化的冰雪水注入长江，使得其他区域陆地水储量增多。

4 结语

本文利用隐含层激励函数为ReLU的BP神经网络算法填补了ITSG-Grace2018和ITSG-Grace_operational模型间的数据空缺，基于填补后完整的数据反演了长江流域陆地水储量变化，并结合GLDAS水文模型、降水和温度等数据，综合分析近20 a来长江流域水储量的时空变化特征，得到如下结论：

1) 对不同隐含层激励函数的BP与RBF神经网络进行对比分析表明，以ReLU作为隐含层激励函数，GLDAS模型的地表水、降水和气温异常量作为输入因子的BP神经网络拟合精度较高，利用该神经网络方案与GRACE模型分别计算的陆地水储量变化的相关系数在长江流域不同时段均不低于0.96，且其差值的标准差均小于1.00 cm，说明利用该神经网络方案填补GRACE与GRACE-FO之间的空缺是可行的。

2) 在2002~2020年期间，长江流域陆地水储量整体呈上升趋势，速率为2.6 mm/a。长江流域上游中西部区域的水储量明显减小，上游东部与中游明显增大，下游地区基本保持平衡，说明长江流域的水储量在空间上分布不均匀。经分析，长江流域上游中西部区域的陆地水储量下降与长江发源地唐古拉山冰川消融有关，长江流域上游东部与中游陆地水储量上升与降水量增大、冰川融化、三峡大坝的运行等有关。长时间序列的GRACE/GRACE-FO时变模型可监测到长江流域2019年的干旱与2017年、2019年的洪涝等灾害。

致谢: 感谢ICGEM提供时变重力场模型数据，GSFC提供GLDAS水文模型数据，中国气象数据网提供降水与气温数据，冯伟研究员提供GRACE_MATLAB_Toolbox软件。