目前陆地水储量的观测方法主要是利用遥感卫星观测反演和以气象和水文观测资料为基础，结合相关物理规律的模式模拟(钟敏等，2009)，前者只能得到十几厘米厚度的土壤含水量，后者在观测资料稀疏的地区不确定性较大，导致大气与陆地水文模式输出的陆地水量变化的结果存在一定的差异； GRACE(Gravity Recovery and Climate Experiment)重力卫星的出现，极大程度上弥补了遥感卫星只能观测地表十几厘米厚度的土壤湿度和地表观测台站空间分布不均匀等不足，它在全球观测分布均匀，并且观测尺度统一，为定量研究大、中尺度陆地水储量的变化提供了可能.Swenson等(2003)研究发现GRACE时变重力场反演中长空间尺度地表水储量变化时，当流域区域大于4.0×10⁵ km²时，所估计的水储量周年变化精度可达到0.7 cm；当流域区域为3.9×10⁶ km²或更大时，其精度可达到0.5 cm，并可通过定量确定流量和水储量来完善水文模型.Swenson等(2006)在面积约为2.8×10⁵ km²的Illinois区域，首次将GRACE反演水储量结果与实测的土壤水和地下水数据进行了验证比较，两者的不符值均方根误差在2 cm左右.杨元德等(2009)利用GRACE的RL04数据分析了2003年1月—2007年12月全球27条流域和陆地水储量的季节性和年际变化，GRACE与GLDAS数据均表明2006年后陆地水储量年际变化明显增加.近年来，利用GRACE卫星时变重力观测研究中国区域陆地水和地下水储量变化的工作也陆续展开(Su et al.，2011；罗志才等，2012；李琼等，2013；尼胜楠等，2014).Hu等(2006)利用2002年4 月—2003年12 月的GRACE数据分析了长江流域水储量变化，与GLDAS和CPC水文模型结果差别小于1 cm等效水高.罗志才等(2012)利用2002年8 月—2011年6 月的时变重力场模型GRGS-EIGEN-RL04反演了黑河流域陆地水储量变化，给出了黑河流域地下水储量的时空变化，并利用张掖地区23口地下水测井数据对地下水反演结果进行了初步验证.尼胜楠等(2014)利用2002年4 月—2011年9 月的GRACE数据反演了长江、黄河流域水储量变化，与水文模型(GLDAS、WGHM)结果相关系数达到0.78以上.

黑河流域地处河西走廊和祁连山中段，是我国西北典型的内陆河流域，总面积约为2.71×10⁵ km²，主要受中部高纬度的西风带环流控制和极地冷气团影响，气候干旱少雨，生态环境脆弱，水资源短缺，属于典型的资源型缺水地区.因而，研究黑河流域陆地水储量变化，包含了水循环过程中降雨量、蒸散发量、地表径流、土壤水与地下水含量等各部分变化的总和，不但可以为水文模型提供数据，也可验证模型输出结果的精度；此外，测量水储量的变化还有助于对干旱洪涝等灾害进行监测及评估(曹艳萍和南卓铜，2011).

尽管杨元德等(2009)考虑到利用水文模型对信号泄漏进行改正，但时间跨度短，其他文献则均没有考虑.Chen等(2014)针对GRACE数据处理中，由于采用滤波方法造成的信号衰减和泄露效应，提出了恢复初始信号、增强被湮没信号以及减小泄漏信号误差的Forward-Modeling方法，本文主要侧重于对Forward-Modeling方法的实现，以及利用该方法对GRACE RL05数据反演得到的黑河流域陆地水储量长时间序列(2003年1月—2013年12 月)的变 化结果进行了恢复改正，并与全球水文模型GLDAS(Global L and Data Assimilation System)、CPC(Climate Prediction Center)结果作了对比，对其时空分布变化特征进行了精细分析，给出了新的物理解释.同时，联合GRACE数据和GLDAS数据反演了黑河 流域地下水储量变化，与全球降雨数据GPCC(Global Precipitation Climatology Centre)进行了比较分析，两者的相关性达到0.88以上.

GRACE提供的时变重力场反映了地球系统的质量分布及变化，在陆地区域，陆地水的变化是引起重力场变化的主要原因.根据重力场球谐系数反演以面密度变化表示的地表质量变化的模型为(Wahr et al.，1998; Swenson and Wahr，2002)：

一般采用等效水高(Equivalent-Water-Height)ΔH(θ，λ)来表示地球表面质量密度变化(单位cm)：

Chen等(2005)发现对于重力场球谐系数的截断和高斯平滑都会不同程度削弱原始信号，引起信号失真(泄漏).本文计算模拟了黑河区域重力场信 号(以等效水高变化率表示)的处理过程，如图 1所示.

图 1

Fig. 1

图 1 不同条件下信号数据处理结果 Fig. 1 The result from different condition

图 1a为模拟的初始信号，图 1b为球谐阶数截断到60阶的结果，图 1c为60阶内300 km半径滤波结果，图 1d为60阶内500 km半径的滤波结果.发现在上述处理过程中，信号明显削弱且泄漏到周围区域.针对这种泄漏误差的改正，Chen等设计了一种正演模型(Forward-Modeling)，通过以下4个步骤可以恢复到最初的信号信息(Chen et al.，2005，2009，2011，2013，2014)：

(1)对于1°×1°的每个格网点的真实质量变化值(或质量变化率)，可以令其初始值等于零或等于GRACE算的等效水高(或等效水高变化率)，结果记为Model00；

(2)将(1)中的初始值采用与GRACE数据同样的处理方法：对全球1°×1°格网点数据进行球谐展开，得到正则化的球谐系数，系数截断到60阶次，0阶和1阶项设为零，然后进行半径为500km的高斯滤波，结果记为Model01；

(3)对于每一个格网点的质量变化(或变化率)值，将GRACE解算的结果与Model01之间的差值加到Model00，得到新的初始值，然后重复步骤(2)；

(4)当GRACE与Model01之间的差值低于某一个特殊的阙值或者迭代一定次数，则停止迭代.

图 2为采用Forward-Modeling方法，对黑河区域模拟数据的500 km滤波结果进行改正前后的对比效果，图 2a为500 km半径滤波后的信号结果，图 2b为Forward-Modeling方法恢复后的结果，即Model00，图 2c为得到的真实信号500 km滤波处理后结果，即Model01，图 2d为GRACE的结果与Model01的差值，差值很小时，认为恢复的真实信号 是可信的.由图 2a、2b可以看出，经Forward-Modeling 方法后，泄漏明显收敛，恢复效果显著，观测信号采用Forward-Modeling方法前后与真实信号不符值 的均方根误差分别为0.14 cm·a^-1和0.09 cm·a^-1，再结合图 2d的误差值很小，说明图 2b恢复的信号可以明显减小泄漏效应带来的误差.

图 2

Fig. 2

图 2 信号经Forward-Modeling方法改正结果 Fig. 2 The correction based on Forward-Modeling

(1)GLDAS水文模型

GLDAS水文模型是由美国宇航局哥达航空中心(Goddard Space Flight Center，NASA)和美国国家环境预报中心(NCEP，National Centers of Environmental Prediction)共同建立的全球水文模式(Rodell et al.，2004)，通过近实时的地面和空间数据约束其模型，利用数据同化技术，输出陆地表面各项参数(如土壤湿度和温度、蒸发量等)，进而获得陆地表面变化的近实时信息，主要反映土壤水分和冰雪变化.空间分 辨率为1°×1°，每月一值，土壤水分为0～10 cm、10～40 cm、 40～100 cm和100～200 cm四层.

(2)CPC水文模型

CPC水文模式是由美国国家海洋和大气管理局(NOAA)气候预报中心(CPC)提供的陆地同化数据集(Fan and van den Dool，2004)，根据全球观测到的降水分布而建立的，采用每小时和每日的降水分析结果、太阳辐射分布、地表大气压、温度和水平风速等数据，提供的产品包括土壤温度和地表以下4层(0～1.6 m)的土壤水含量等，该模式输出结果的空间分辨率为0.5°×0.5°，每月一值.

本文采用UTCSR(University of Texas at Austin，Center for Space Research)提供的2003年1 月—2013年12月共124个月(2003年6月，2011年1、6月，2012年5、10月以及2013年3、8和9月数据缺失)的GRACE RL05数据对黑河流域水储量变化进行了反演.首先采用SLR获得的C₂₀项替代GRACE数据的C₂₀项，求得球谐系数平均值，并将124个月的球谐系数分别减去平均值得到球谐系数变化量(ΔC_lm，ΔS_lm)，然后对时变重力场采用平均半径为500 km的高斯平滑和P4M8去相关滤波，得到等效水高变化量.最后采用Forward-Modeling 方法对数据结果进行了改正，并与水文模型GLDAS、CPC 的结果进行了比较，如图 3所示.计算中，缺失的月份进行插值得到.

图 3

Fig. 3

图 3 GRACE、Forward-Modeling、GLDAS和CPC得到的黑河流域水储量变化 Fig. 3 Results of the Heihe River storage variation from GRACE, Forward-Modeling, GLDAS and CPC

图 3a为等效水高表示的黑河流域水储量月变化，变化幅度在±5 cm内，图 3b为周年、半周年变化，振幅变化在±1.5 cm.从图 3a的时间序列上看出2003年初到2006年底呈现下降趋势，2007年初到2010年底呈现平稳变化，而从2011年初到2013年底则呈现上升趋势，而水文模型的周期性变化则比较明显.从图 3b的周期项可以看出每年的4月份水储量亏损达到最大，8月份盈余量最大，而相应的水文结果则要滞后一个月.这可能与模型本身的建立相关，对于时变信号的反映比较滞后.比较结果如表 1、表 2所示.

表 1(Table 1)

表 1 黑河流域水储量变化等效水高的周期项特征及不同数据结果间的相关系数 Table 1 The trend of TWS change in the Heihe River and the correlation coefficients of different methods 周期项及相关性 周年项 半周年项 GRACE-Observed Forward-Modeling GLDAS 振幅/cm 相位/(°) 振幅/cm 相位/(°) GRACE-Observed 1.2 -3.8 0.4 -1.3 1 Forward-Modeling 1.5 -3.5 0.8 -1.1 0.95 1 GLDAS 0.4 -1.9 0.2 -0.9 0.53 0.64 1 CPC 1.3 -2.2 0.3 -0.4 0.63 0.71 0.89

表 1 黑河流域水储量变化等效水高的周期项特征及不同数据结果间的相关系数 Table 1 The trend of TWS change in the Heihe River and the correlation coefficients of different methods

表 2(Table 2)

表 2 黑河流域水储量变化率 Table 2 The TWS rate in the Heihe River area 趋势项/(cm·a-1) 不同时间段 2003年1 月—2006年12 月 2007年1 月—2010年12 月 2003年1 月—2010年12 月 2011年1 月—2013年12 月 GRACE-Observed -0.13 0.14 -0.06 0.56 Forward-Modeling -0.86 0.07 -0.35 0.14

表 2 黑河流域水储量变化率 Table 2 The TWS rate in the Heihe River area

从表 1结果看，GRACE的结果与水文模型(GLDAS、CPC)的相关性分别为0.53和0.63，采用Forward-Modeling后相关性有所提高，达到0.63和0.71.从周年项上看，CPC的结果与GRACE和Forward-Modeling的更为接近，GLDAS的振幅要偏小一些，这估计与模型的建立有关，水文模型只包含了土壤水文和积雪，并不包含地表水和地下水.在半周年项上，Forward-Modeling的结果要高于其他三种结果，为GRACE结果的2倍，这也说明了采用这种方法对信号泄漏的恢复在半周年上更为明显.表 2中给出了不同时间段的水储量变化率，两者的相关系数为0.95，这说明对GRACE的结果恢复是可信的，可能正是由于半周年振幅的差别，导致相同时间内的变化率大小不一致，为此给出了不同时间段的变化率空间分布，如图 4所示.

图 4

Fig. 4

图 4 黑河流域水储量变化率在不同时间段内的空间分布 Fig. 4 The spatial distribution of TWS variation rate in the Heihe River from different periods

图 4a为GRACE的结果，图 4b为Forward-Modeling 的结果，从上至下分别为2003—2006、2007—2010和2011—2013三个时间段的水储量变化率空间分布.在2003—2006时间段内，水储量变化率呈现从区域的西南方向到东北方向依次递减的阶梯分布，区域范围内主要以下降为主，在东北区域达到最大-0.5 cm·a^-1，采用Forward-Modeling方法后区域内的这种亏损信号得到了增强，在保持东北区域-0.5 cm·a^-1下降趋势和西南区域0.5 cm·a^-1上升趋势的情况下，最大亏损集中出现在区域中部，达到了-1.5 cm·a^-1，可能正是这种初始信号向周围的泄漏，导致了之前出现的阶梯分布.在2007—2010年间，GRACE结果同样仍呈现西南到东北依次下降的趋势，如果把区域从中间分为东、西两部分，可以看出西部以盈余为主，东部以亏损为主，而采用Forward-Modeling方法后在西北区域出现了部分亏损，以及东部的亏损集中到了东南区域，整个东北部区域保持稳定，变化率基本为0.西北部亏损的出现可能由于整个西部两端的盈余泄漏引起的，把这种信号湮没了.从2011—2013年，可以看出整个分布格局出现了由西北到东南依次增长的趋势，采用Forward-Modeling后这种趋势更为显著，亏损由西北延伸到了区域中部，说明了该方法有助于恢复被湮没的信号.结合表 2给出的变化率，对于2003—2006年间正是由于区域中部亏损信号的恢复使得Forward-Modeling结果远大于GRACE结果，对于2007—2010和2011—2013年间，Forward-Modeling的值要小于GRACE的结果，或许也正是对部分湮没了的亏损信号恢复，使得整体变化减缓.因此可以认为表 2给出的结果是可信的.

黑河流域作为西北干旱地区较大的内陆河流域，其水资源特性决定了地下水通常是干旱地区最重要的水源和供水选择，准确判断地下水在干旱地区的空间分布及其变化程度是制定区域地下水合理利用规划的关键.由GRACE反演得到的陆地水储量变化包含了地下水、冰雪、土壤水、地表水(水库、河流等)和地表生物含水量的变化，但考虑到黑河流 域实际情况和现有资料，利用全球水文模型GLDAS 得到黑河流域土壤水和冰雪变化，从GRACE时变重力反演的陆地水储量变化中扣除土壤水和冰雪变化即可得到黑河流域的地下水储量变化.

根据水量平衡方程ΔS=P+Q－E－R－W，其中ΔS为流域陆地水储量变化，P为降水，Q为冰川融水，E为蒸散发，R为地表径流，W为人类活动用水(曹艳萍和南卓铜，2011).黑河流域作为内陆河流域，除部分冰川融水外，降水对水储量变化起主导作用，是地下水补给的主要来源.为对比黑河流域地下水储量变化与降水变化的相关性，本文采用全球 降水气候学中心(GPCC)提供的2003年1月—2013年12月的全球180°×360°月降水格网数据(与GRACE 水 储量数据具有相同的分辨率)(Schneider et al.，2014)，得到黑河流域降水变化(以降水量距平值表示，单位cm).

图 5a为黑河流域地下水储量2003—2013年的时间变化序列和扣除平均值后的降水变化序列，对于地下水的获取本文采用了两种方法：第一种是将水文模型GLDAS的格网数据转换为与GRACE重力场同阶次的球谐系数，采用GRACE数据的处理方法获得土壤水和冰雪变化量，然后从GRACE结果中扣除得到地下水变化量；另一种是将GLDAS的数值除以水密度转换为等效水高值，然后从Forward-Modeling结果中直接扣除得到地下水变化.图 5b为地下水变化和降水变化的周年特征曲线.可以看出，由GRACE月模型直接计算得到的黑河流域地下水变化结果、通过Forward-Modeling方法计算得到的地下水变化结果以及扣除平均值后的降水变化序列等三种结果在周年变化上比较接近.结合表 3，尽管GRACE的结果与GPCC的相关系数达到了0.96，而Forward-Modeling的结果只为0.88，但从振幅和相位的总体趋势分析可以看出，Forward-Modeling 结果与GPCC更为接近.采用Forward-Modeling后，在周年和半周年项上，振幅均有显著提高.造成这样的原因可能是某些月份经过泄漏改正后与降水差别加剧，降低了数据间的相关性，但降水只是地下水变化的主要因素，它受到冰雪融化、蒸发、径流等影响，这种经过恢复后差距的加大，也许正体现了这些信号的影响.在表 3中同时给出了趋势项，在2007年以后，两者结果比较相近，而对于2003—2006年间则存在明显差异.如图 6所示，给出了不同时间段内的地下水储量变化率的空间分布.

图 5

Fig. 5

图 5 黑河流域地下水储量变化和降水量变化 Fig. 5 The TWS variation of the Heihe River and its GPCC result

表 3(Table 3)

表 3 黑河流域地下水储量变化的周期性特征 Table 3 The trend of TWS variation in the Heihe River 周期项及趋势 周年项 半周年项 相关系数GPCC 趋势项/(cm·a -1) 振幅/cm 相位/(°) 振幅/cm 相位/(°) 2003—2006 2007—2010 2011—2013 GRACE-Observed 0.7 -4.6 0.3 -1.8 0.96 0.22 0.14 0.31 Forward-Modeling 1.4 -4.1 0.6 -1.3 0.88 -0.12 0.16 0.32 GPCC 1.4 -4.8 0.5 -1.3 1 — — —

表 3 黑河流域地下水储量变化的周期性特征 Table 3 The trend of TWS variation in the Heihe River

图 6

Fig. 6

图 6 黑河流域水储量变化率在不同时间段内的空间分布(扣除GLDAS模型) Fig. 6The spatial distribution of TWS variation rate in the Heihe River from different periods (GLDAS processed)

图 6a为GRACE扣除GLDAS滤波处理后的地下水变化率分布，图 6b为Forward-Modeling的结果，从上至下依次为2003—2006、2007—2010和2011—2013三个时间段.从图 6a中可以看出，2003—2006年期间 地下水变化率整体保持正增长，中部区域基本保持了零变化，地下水变化率最大变化集中在东南和西南区域，达到0.4 cm·a^-1；而图 6b则分辨比较明显，中部区域被湮没的负增长信号得到恢复，最大可达到-1 cm·a^-1，从而使得区域呈现下降趋势.对于2007—2010和2011—2013时间段，从图 6b恢复结果可以看出，正是部分区域正增长信号的泄漏，使得周围亏损信号被湮没，采用Forward-Modeling方法可以有效改善这种结果，增强信号.另外，从3个时间段恢复后的结果来看，都基本保持了上游区域正增长的趋势，该地区位于祁连山脉，植被条件好、降雨多、气温低、蒸发量小，且受到冰雪消融的影响，从而保证了地下水的盈余.中部多为干旱、大风和沙漠化区域，降雨少、蒸发量大、天然植被少，使得地下水下降比较严重.

本文通过计算发现，利用GRACE数据反演陆地水储量变化的处理过程中进行球谐阶数截断、滤波都会不同程度的引起原始信号泄漏，从而会湮没周围区域的信号，造成信号失真.针对此，文中采用了Chen等提出的Forward-Modeling方法对黑河流域陆地水储量(包括地下水储量)进行了泄漏改正，增强被湮没的信号，使得原始信号得到有效恢复，改善反演结果.并将本文结果与GLDAS和CPC水文数据结果以及GPCC降水数据结果进行了比较分析，研究结果表明：

(1)泄漏误差改正的量级平均接近1 cm(约为8 mm)，是利用GRACE数据反演水储量变化需要考虑的重要因素；并且采用Forward-Modeling方法后可有效地恢复初始信号，增强被湮没的信号，提高水储量变化的空间分辨效果.

(2)通过对黑河流域2003—2013年长时间序列的水储量变化分析，发现其变化具有明显的阶段性特征，在2003—2006年水储量下降趋势严重，可达到-0.86 cm·a^-1，在2007—2010年间基本保持 不变，然而从2011—2013年则呈现缓慢上升趋势，为0.14 cm·a^-1并且在周年季节变化上与水文结果较为接近，相关系数达到0.64以上，在8、9月份达到最大盈余，4、5月份亏损最为严重.

(3)降雨作为黑河流域地下水补给的主要来源，将GPCC月降雨变化量与地下水储量变化比较发现，两者相关性可达到0.88以上，在周年季节变化上基本一致，Forward-Modeling结果与GPCC在振幅和相位更为接近，在秋冬季普遍下降，2月份亏损最大，春夏逐渐增加，7月份达到最大盈余.