青藏高原水资源变化是国家经济和社会发展最关心的问题之一，陆地水储量（terrestrial water storage，TWS）变化的准确估计对气候变化、农业生产、洪灾以及其他自然灾害的研究和预报具有重要意义^[1]。重力反演与气候实验(gravity recovery and climate experiment, GRACE)卫星计划的出现，为监测大范围陆地水储量变化提供了有效技术手段，极大程度上弥补了遥感卫星只能观测地表十几厘米厚度土壤水、观测台站空间分布不均匀等不足。水储量变化通常采用滤波球谐系数估算，以降低GRACE高阶噪声对反演结果的影响^[2]。

近年来，利用GRACE滤波时变重力场监测青藏高原水储量变化的研究工作陆续开展，相关研究主要集中在水储量变化特征分析等方面^[3-5]，而针对该地区滤波影响的定量研究则相对较少。滤波影响表现为研究区信号泄漏，振幅明显衰减^[6]。为减小泄漏误差，目前多数研究采用尺度因子法恢复滤波衰减信号^[7]。如Klees等^[8]提出了多时刻尺度因子法；冯伟等^[9]研究了GRACE反演亚马逊流域陆地水变化的单一尺度因子法。两种方法均较依赖水文模型。为减小水文模型不确定性影响，Chen等^[10]提出了Forward‐Modeling方法。

本文侧重Forward‐Modeling方法的实现，针对GRACE泄漏误差影响开展定量研究，并采用Forward‐Modeling方法改正，验证该方法在青藏高原地区的适用性。在顾及冰川均衡调整(glacial isostatic adjustment, GIA)效应情况下，估算水储量时空变化，将结果与全球水评估与预测模型（water-global assessment and prognosis hydrology model, WGHM）进行比较。分析基于短时傅里叶变换（short-time Fourier transformation, STFT）提取的时序信号时频谱特征，并结合全球降水气候中心（global precipitation climatology centre, GPCC）降水数据探讨水储量与降水的关系。

1 数据和反演方法 1.1 GRACE卫星反演方法

GRACE月重力场模型采用美国德克萨斯大学空间研究中心提供的2级产品第5版本的2003—2014年数据。该数据为正则化球谐系数，扣除了固体潮、海潮、固体极潮、海洋极潮、非潮汐大气和海洋影响，以及日月等行星引起的重力扰动；C₂₀项由卫星激光测距系数代替^[11]，同时改正一阶项。

为避免高阶次噪声及奇偶阶相关误差的影响，本文截断60阶位系数，同时采用高斯滤波与15阶以上位系数通过3次多项式拟合去相关滤波法进行去条带处理。经过处理后，位系数最终可转换为分辨率为1°×1°的全球分布等效水高变化：

$\begin{array}{*{20}{c}} {\Delta \sigma (\theta , \lambda ) = }\\ {\frac{{R{\rho _e}}}{3}\sum\limits_{n = 0}^N {{W_n}} \sum\limits_{m = 0}^n {\frac{{2n + 1}}{{1 + {k_n}}}} \left[ {\Delta {C_{mm}}\cos (m\lambda ) + } \right.}\\ {\left. {\Delta {S_{mn}}\sin (m\lambda )} \right] \cdot {{\bar P}_{mn}}(\cos \theta )} \end{array}$

(1)

式中，(θ, λ)为计算点地心余纬和经度; $\left({{\rm{\Delta }}{C_{nm}}, {\rm{\Delta }}{S_{nm}}} \right)$为经P3M15去相关滤波方法处理后的n阶m次归一化位系数变化；${\bar P_{nm}}\left(\cdot \right)$为规格化的n阶m次缔合勒让德函数；k_n为n阶负荷勒夫数；${\rho _e} \approx 5.5 \times {10^3}\; {\rm{kg}}/{{\rm{m}}^3}$为固体地球平均密度；R为地球平均半径；W_n为n阶高斯平滑系数。

1.2 水文模型

水文模型采用德国法拉克福大学自然地理研究所提供的2.2c版本的WGHM模型^[12]。该模型联合月时序参数以及气候和地理数据建立而成，数据时段为2003—2014年，空间分辨率为0.5°×0.5°，每月一值。

1.3 降水模型

降水模型采用GPCC提供的2015年发布的全球降水数据（GPCC full data reanalysis version 7）。该版本数据准确度较高，更适合应用于水文气象模型验证和水循环研究，数据更新至2013年12月^[13]。空间分辨率为1°×1°，每月一值。

1.4 GIA模型

冰川均衡调整(GIA)是指黏弹地球对末次冰进期与冰退期地表冰和海水负荷改变的响应，是一种重要地球动力学现象。GIA过程主要体现在地幔物质运动、地壳运动及地球重力场变化。本文GIA效应采用Paulson等^[14]的地壳均衡模型扣除。

2 Forward-Modeling方法

为了定量估算滤波泄漏误差影响，验证Forward-Modeling方法在研究区的适用性，本文模拟了区域重力场信号经滤波及泄漏信号恢复后的结果。研究区位置分布如图 1所示。

由于去相关滤波影响难以评估，且方法本身具有对数据自适应性和非线性特点^[10]，因此本文着重对高斯滤波信号影响进行恢复。通过对模拟初始信号作球谐展开得到重力场60阶模拟球谐系数，进一步利用式（1）计算研究区的模拟水储量变化1°×1°格网分布，如图 2所示。图 2(a)为模拟初始信号，图 2(b)~图 2(d)分别为球谐展开及不同平滑半径滤波结果。由图 2可知，滤波后初始信号振幅明显被削弱，且信号向周边区域有所泄漏。500 km滤波对初始信号损失较严重。相关研究表明，滤波半径越大时，高阶球谐系数权重越小，从而降低了GRACE反演水储量变化的空间分辨率^[15]。因此本文选择了300 km滤波核函数，采用Forward‐Modeling方法对滤波泄漏误差进行改正，其迭代过程如下^[16]：（1）将图 2(c)的模拟结果作为研究区GRACE的估算结果，记作M，并将M作为Forward‐Modeling方法所要恢复的信号I的初始值。（2）将信号I球谐展开至60阶，令零阶与一阶项为零，经滤波后计算相同区域1°×1°格网分布，记作T。（3）将M与T按对应格网作差，将差值乘以1.2（1.2是经验值^[16]，以加快迭代速度）回加到步骤（1）的I，作为新初始值，进行新一次迭代。当差值小于某一特定阈值或迭代到一定次数，迭代停止。

基于以上步骤对模拟信号滤波泄漏误差进行改正，结果如图 3所示。图 3(a)为GRACE滤波结果M；图 3(b)为Forward-Modeling方法恢复的信号I；图 3(c)为信号I经球谐展开及滤波后得到的结果T；图 3(d)为图 3(a)与图 3(c)的结果差值。由图 3(d)可知，多次迭代滤波结果与GRACE解算结果差异很小。当差值很小时，可认为恢复信号I是可信的。由图 2和图 3可以发现，经Forward-Modeling方法恢复的信号其泄漏误差明显收敛，恢复效果较为显著。GRACE解算结果与图 2(a)真实信号的偏差均方根为0.1 cm/a，泄漏误差改正后结果均方根为0.07 cm/a，精度提高了近30%。上述结果说明Forward-Modeling方法可适用于青藏高原地区GRACE滤波泄漏误差的处理。本文进一步采用该方法定量估计研究区水储量变化。

3 结果分析 3.1 水储量空间变化

为了定量研究GRACE泄漏误差影响，本文比较分析了经泄漏误差改正前后的青藏高原水储量时空变化。在顾及GIA均衡效应情况下，利用滤波时变重力场按式（1）估算了研究区水储量变化的等效水高。进一步对不同时段水储量的空间变化拟合其变化速率。拟合的时段包括2003—2014年、2003—2004年、2005—2010年、2011—2014年。

为了与GRACE结果对比，同样对WGHM模型作球谐展开，其展开阶数与GRACE相同。各时段变化速率如图 4~图 7所示。图 4(a)、图 5(a)、图 6（a）、图 7（a）为GRACE滤波估算的水储量变化速率；图 4(b)、图 5(b)、图 6（b）、图 7（b）为改正泄漏误差后水储量同期变化速率；图 4(c)、图 5(c)、图 6（c）、图 7（c）为WGHM模型结果变化速率。由图 4~图 7可知，恢复泄漏信号后，研究区水储量各时段空间分布变化仍然较大，具体表现为：（1）2003—2014年期间，水储量分布格局整体上仍保持从南到北、从西至东依次增长的趋势；东南区域水储量的最大增长约为0.6 cm/a，东北部亏损加剧，正北区域盈余信号加强。（2）2003—2004年，东北部水储量快速增加；正南、西南区域的正信号幅值明显增大，最大达到6 cm/a。水储量亏损信号依然从西北延伸到东南部，亏损速率降低至-6 cm/a。（3）2005—2010年，西南至东北方向水储量亏缺严重，达到-2 cm/a；西北、东南区域水储量逐年增加，最大增长率达到2 cm/a。（4）2011—2014年，正北、东南区域水储量正增长信号显著增强，最大达到4 cm/a；亏损信号仍集中在整个区域的正南、西南区域，最大亏损速率为-4 cm/a。

可以发现，泄漏误差影响主要体现对水储量信号振幅削弱，采用Forward-Modeling方法之后，泄漏误差明显减小，信号振幅绝对值增大近2倍，恢复效果显著；水储量变化具有更明显的空间差异性，个别时段空间变化速率最大值与最小值之差可达12 cm/a。通过与水文模型的比较，发现GRACE结果空间变化特征与WGHM模型基本吻合。但在局部长期变化速率上仍存在较大差异，可能由于构建WGHM模型时采用的全球气候强度数据在局部空间上具有不确定性，从而导致模型本身存在一定误差。

3.2 水储量时间变化

为进一步分析泄漏误差对水储量时间变化的影响，本文对改正泄漏误差前后的研究区水储量时间变化进行定量描述，时间序列可通过对GRACE结果计算空间平均获得。图 8中蓝色曲线为GRACE滤波估算的研究区水储量变化，红色曲线为恢复泄漏信号后水储量变化，绿色曲线为WGHM模型估算结果。

为了定量分析GRACE改正泄漏误差前后以及WGHM模型结果时间序列的周期项变化，本文拟合了周年项、半周年项、S₂潮汐项，这里只给出周年与半周年信号振幅、相位，如表 1所示。关于水储量时间变化特征，Ray等^[17]通过研究卫星重力潮汐模型，发现S₂、K₂、K₁潮汐混频问题分别会产生周期为161 d、3.7 a、7.4 a的影响。本文研究时间尺度为12 a，故这里不考虑K₁、K₂潮汐引起的7.4 a和3.7 a周期。研究区格网点水储量变化的等效水高可通过最小二乘算法拟合，公式为：

$\Delta h = c + Rt + \sum\limits_{k = 1}^3 {{A_k}} \sin \left( {{w_k}t} \right) + \tau $

(2)

式中，t为时间变量；c为常数；R为长期速率；A_k为周期项w_k的振幅，其中k= 1和k= 2时分别表示周年与半周年项，k= 3时表示与S₂半波相关的161 d周期项；S₂为拟合误差。

结合图 8和表 1可以发现，改正泄漏误差对于滤波信号振幅的增强具有显著效果。泄漏信号恢复后，水储量变化季节性波动明显，周年与半周年振幅分别达到约7.5 cm、2.7 cm等效水高；时序异常特征信号更为凸显。在2004年和2010年，水储量实现了12年中最大盈余和亏损。周年相位与WGHM结果十分接近，年内变化趋势基本一致。但水文模型结果振幅要偏小一些，这可能与模型的建立有关，WGHM模型虽已包含土壤水、积雪、地表水以及地下水等多种水文要素，但其构成含量可能相对较小。然而泄漏信号恢复结果其半周年相位存在明显滞后，原因有待研究。此外，研究区水储量在不同时段具备不同的变化趋势，如图 8中各趋势项所示。其中细虚线表示泄漏误差改正之前的结果，粗虚线表示泄漏误差改正之后的结果。各时段变化率如表 2所示。

由表 2可知，改正泄漏误差前后，研究区水储量在2003—2014年均呈现长期弱下降趋势，泄漏误差对整体变化趋势影响较小。结合图 3发现，这种小幅变化可能与部分区域被淹没的正增长信号得到恢复有关，从而导致总体变化趋势减缓。从分段趋势上看，2003—2004年期间，泄漏误差改正后，研究区水储量表现为快速增长。正是由于该时段泄漏的正信号被有效恢复，整体水储量增加速率显著提升。2005—2010年，陆地水储量逐年减小，泄漏信号恢复后，水储量逐年亏损趋势有所加剧。2011—2014年，陆地水储量变化相对平稳，改正后水储量呈现明显的正增长趋势，这可能与局部区域水储量盈余信号的大幅加强有关。

为了展现信号在时频域上的能量分布，本文引入STFT算法获取研究区陆地水储量时间变化的时频谱特征，如图 9所示。

图 9中时间轴表示144个月时序；频率轴表示周期项；幅值轴则表示时序在时频域上的能量强度。经比较发现，改正泄漏误差后，GRACE结果在低频域的信号能量显著增强，这表明GRACE估算区域水储量变化是以低频特征信号为主，季节性变化波动明显。同时陈超等^[18]也监测到该区域南部地表存在明显的季节性垂直形变。但在2004年8月，水储量的季节信号强度开始骤减，接近2006年4月，减小幅度有所缓解并开始逐渐回升，说明该时段内研究区水储量季节性变化相对薄弱。由此可见，泄漏误差对水储量季节信号强度影响较大。

3.3 降水影响

降水是研究区陆地水补给的主要来源，本文进一步分析了两者之间的关系，采用2003—2013年GPCC降水数据得到研究区降水变化。首先将泄漏误差改正后的结果移去各周期项信号，包括季节项、S₂潮汐项。其中，S₂潮汐项按照式（2）进行计算与扣除，再采用3个月滑动平均得到水储量年际变化，如图 10所示。

图 10显示了研究区每年6月—9月存在大幅持续降水，GRACE监测结果较降水数据时间上基本一致。陆地水储量年变化特征与降水的年变化规律比较吻合。特别是在2006年和2009年，该地区降水量较多年平均量减少约75 mm，同时GRACE监测的陆地水储量均发生了下降突变。2003—2004年期间，研究区降水明显增多而后增幅减缓，同期水储量也表现出快速回升之后波动稳定。2007—2008年降水偏多，水储量持续下降的趋势得到抑制，并在2008年末出现小幅回升。2010年降水大幅增多，同期陆地水储量回升显著。2011—2013年降水持续小幅增多，陆地水储量增加的趋势相对平稳。

综上所述，改正泄漏误差后，研究区水储量时空分布与WGHM及降水模型的变化特征整体上比较一致。可以说明，利用Forward-Modeling方法改正泄漏误差后的GRACE反演结果及数据在表征青藏高原水储量变化方面具有一定的可靠性。

恢复泄漏信号后的GRACE结果在时空分布上与WGHM模型仍然存在局部差异，这种差异可能由二者空间分辨率不同，以及GRACE数据后处理产生的过程误差和WGHM模型本身不确定性综合引起。另外，GRACE反演的研究区整体水储量年际变化与降水在部分时段内的相关关系比较薄弱，这可能与青藏高原陆地水、冰川、湖泊、地下水变化规律以及受降水的影响机理不同有关。Yi等^[19]分别讨论了冰川、湖泊和地下水的信号变化，其中，南部冰川主要受夏季印度季风产生的强降水影响，帕尔米高原冰川受冬季西风带降水影响，高原内部冰川主要受大陆气候条件影响。总体上，冰川和地下水变化与降水具有强相关性，相关系数分别可达0.59、0.89。此外，Yi等^[19]还发现湖泊水位上升是导致高原内部存在较大正信号（+30 Gt/a）的原因之一，并指出同期大幅降水可能是引起湖泊水位增加的重要因素，这种相关性在高原中部体现尤为明显。

4 结语

本文使用GRACE重力数据，针对其滤波泄漏误差影响，采用Forward‐Modeling方法予以改正，并验证了该方法在青藏高原地区的适用性，以此反演研究区水储量时空变化；对泄漏误差改正前后结果进行定量描述与时空特征分析，并将其结果与WGHM模型进行比较；同时采用STFT算法提取时序信号的时频谱特征，结合GPCC降水数据分析水储量和降水的关系。研究结果表明：（1）Forward‐Modeling方法可减小青藏高原水储量变化的滤波泄漏误差。改正泄漏误差后，重力场信号均方根误差降低近30%，可见泄漏误差是利用GRACE数据反演区域水储量变化应考虑的重要因素。（2）青藏高原水储量变化具有明显空间差异性，研究区多数时段以水储量正负变化信号占半，正北、东南部以盈余信号为主，正南、西南部则以亏损信号为主。其空间变化特征与WGHM模型具有较好的一致性。此外，水储量时间序列具有显著的季节性和阶段性变化规律，周年变化振幅达到7.5 cm。2003—2004年期间，水储量增加速率约为3.89 cm/a；2005—2010年，亏损速率为-0.31 cm/a；2011—2014年，上升速率为0.22 cm/a。（3）青藏高原水储量年变化与降水特征吻合较好，在降水明显偏少的2006年和2009年，水储量均发生下降突变；在降水偏多的2003年，水储量逐渐增加而后波动稳定；在降水大幅增多的2010年，水储量显著回升。

致谢: 感谢德国法拉克福大学自然地理研究所Döll教授提供的WGHM (v2.2c)全球陆地水文模型。