0 引言

地表物质的迁移和质量的重新分布会引起重力场的变化，而且同时由于这种地表负荷的变化会引起地表形变，地表质量、地表形变和地球重力场三者的变化可以建立联系，所以通过研究地球重力场的变化可以得到地表质量变化与地表形变信息.2002年3月17日GRACE重力卫星发射，为我们得到更高精度和更高时间分辨率的全球重力场模型提供了途径，使得人们对地球物理的研究更进一步，获得的信息更精细、准确和及时.现今，GRACE重力卫星数据已经广泛应用于研究海洋质量、海洋环流、极地冰川质量变化、土壤湿度、陆地水储量等大尺度地球物信息的研究，同时应用于地表负荷变化形变、同震形变和地壳构造等研究中.

由于GRACE卫星轨道和重力场解算方法等技术特点，我们得到的时变重力场高阶项误差偏大，全球空间精度分布不均等原因为反演造成种种限制，产生诸如条带现象、信号衰减等问题 (Wahr et al., 1998).在利用GRACE数据时根据卫星时变重力场的特点和研究的需要，研究学者提出了多种时变重力场位系数后处理方法.按照处理对象的不同大致分为两类：一类是在频域内对位系数进行处理，如经验性的高斯滤波 (各项同性 (Jekeli, 1981; Wahr et al., 1998)、非各向同性 (Han et al., 2005)、扇形滤波 (Zhang et al., 2009))，去相关滤波 (固定窗口滑动去相关 (Swenson and Wahr, 2006)、变窗口滑动去相关 (Duan et al., 2009)、P4M6(Chen et al., 2006a, b, 2007a, b, 2009)，根据位系数精度的滤波如DDK方法 (Kusche et al., 2009)，球面小波方法 (Schmidt et al., 2008)，维纳滤波 (Sasgen et al., 2007)，统计滤波 (Davis et al., 2008) 等；另一类是在空域内对计算结果根据其空间分布特征提高信噪比，如EOF方法 (De Viron et al., 2008; 周江存等, 2013)，平滑先验信息方法 (詹金刚等, 2015) 等.然而，在上述针对月重力场位模型系数处理过程中，通常会进行截断高阶项以及滤波处理等操作，这将会导致真实信号衰减，或是会在提高信噪比的同时使得真实信号都会有所改变，改变其空间分布与信号值 (高春春等, 2015; 郑秋月和陈石, 2015).这种现象在反演区域水文信号时会导致区域内信号泄漏到区域外，或区域外信号泄漏到研究域内的情况.虽然GRACE反演陆地质量变化通常缺乏精确质量变化量对其检核，但将其获取的反演径向形变与其他监测手段测量 (如GPS) 获取的U方向形变进行对比分析 (魏娜等, 2015)，则可获取GRACE位系数处理过程中信号衰减比例，从而可精化及恢复GRACE反演地表质量变化中的真实信号.

据此，本文选择美国加州中央山谷地区作为研究区域，针对GRACE数据处理中存在的截断高阶项和高斯滤波，以及去条带滤波导致的信号衰减问题，提出了一种附加GPS时序约束的GRACE反演陆地水文信息修正法，以此标定衰减比例恢复真实信号.

1 GRACE反演 1.1 GRACE反演地表质量迁移和负荷形变原理

利用GRACE重力场提供的球谐系数可以计算地表任一点的质量变化及其由此引起的地表位移.根据Wahr等 (1998)理论，GRACE时变重力场能够反演得到地表质量重新分布情况，则地球表面薄层上面密度变化可表示为

(1)

式中, θ和φ分别为地心余纬和东经，a=6378.137 km为地球平均半径，ρ_ave=5517 kg/m³为地球平均密度，l和m为阶和次，为正规化的连带勒让德函数，和分别为时变重力场模型的无量纲正规化球谐系数变化量，k_l为l阶勒夫数.

同时，地表的某一点质量负荷会造成地表弹性形变，地表质量负荷与地表位移间关系Farrell (1972)给出了如下公式为

(2)

(3)

式中u(θ) 和v(θ) 分别为地表垂直向和水平向位移，θ为计算点到质量负荷位置的球面距离，a为地球平均半径，m_e为负荷质量，h_n和l_n为n阶勒夫数，P_lm(cosθ) 为勒让德函数.利用GRACE重力场提供的球谐系数可以计算地表任一点的质量变化及其由此引起的地表位移.在假设地球为弹性体的前提下，质量变化的球谐系数与地表位移的数学关系 (Blewitt, 2003; Blewitt and Clarke, 2003; Wahr et al., 1998) 如下所示为 (本文研究只关注地表径向位移)：

(4)

式中，dr(θ，φ) 为地表径向位移，R为地球平均半径, 数值同公式 (1) 中a、h_l和k_l为l阶负荷勒夫数, 其他字母和函数与公式 (1) 中相同.

在利用上述 (1)、(4) 两公式计算地球表面任意点某时刻的质量变化和负荷形变时，需要利用到共同的和.由于时变卫星重力场模型球谐系数具有高阶精度差，同次系数奇、偶阶各自相关的特点，使得解算的重力场在空间上表现出严重的条带现象 (Swenson and Wahr, 2006).相应的重力场反演结果也受这些误差的干扰，可利用引言中所列诸多方法对其处理.然而，这些方法在处理过程中会削弱目标结果的真实信号强度，本文针对此问题拟加入GPS时序信息对GRACE反演结果进行约束.基本思想为：利用GPS监测地表形变结果 (U向) 与GRACE反演得到的地表负荷形变作比较，从而得到由于球谐系数预处理产生的信号缩减比例，再依照比例来修正GRACE反演的地表质量变化结果.

1.2 GRACE反演计算

本文重力卫星数据采用美国德克萨斯大学空间研究中心 (CSR) 提供的GSM_RL-05月重力场模型，时间跨度为2002年12月至2015年8月14年共139个月的月重力场模型.计算时位系数先进行C₂₀项替换，后经过与连续整年数据2004年1月至2010年12月的月重力场平均值做差，再将差值做P4M6去相关和半径300 km的扇形滤波，最后利用经前处理过的和求得地表质量变化和地表径向位移.在计算地表负荷形变时需要进过特殊处理，首先是计算时需要将GAC系数加入滤波处理后的系数内计算，再者卫星重力场测量坐标系是地球质心坐标系，计算时通过一阶项改正 (Swenson et al., 2008) 将形变转到GPS测量用的CTRF框架内做比较分析.

为验证上述滤波处理方式的合理性，我们利用上述方法求得全球地表质量变化.根据Wahr等 (1998)的理论可将此质量变化认为是地球表层水质量的变化.将质量除以水的密度得到等效水高，从而获得2002年至2015年全球水储量变化速率图 (如图 1).计算结果表明两极地区和格陵兰岛水储量高速率损失，揭示了与全球气候变暖相关的冰川融化现象.此外，美国加州和中国华北地区也表现出了较明显的水储量下降趋势.上述全球水储量变化空间分布与区域性地下水严重耗损等实际情况相符 (Chen et al., 2016)，也在一定程度上表明了本文采用滤波处理方式的合理性.

2 地表形变分析——“缩减”系数计算 2.1 GPS站点形变分析

本文选取了美国西部监测网络中位于加州的监测网PBO数据中解算好的时间序列 (WNAM_Clean_TrendNeuTimeSeries_sopac).该数据为去掉极值以及消除同震形变等影响的数据，原始数据为天解，通过低通滤波消除高频变化得到2002年至2015每月一值的时间序列.在加州地区中选取了800多个GPS站点 (图 2中绿色十字) 记录，并与同址GRACE负荷形变比较.考虑到监测时间与GRACE观测数据时间段相近和两种形变的相关性，选取符合相关系数大于0.6且GPS监测时间段在2002年到2014年这两个条件的58个点 (图 2中红色圆圈).在这些点中选择具有代表性站点SHIN (位置如图中所示) 结果分析如图 3、图 4所示.图 3中红色的曲线代表由GRACE数据计算得到的地表径向负荷形变，时间段为CSR所公布的月重力场始末，其中除了2004年至2010年为连续数据外其他各年份均有中断.图中黑色曲线代表GPS检测结果，由于考虑到所有站点，数据时段选择始末选在2002年至2015年，GPS数据也有中断，各站点情况有所不同，特别是去月值后中断情况明显减少，且整体连续性好于GRACE计算结果的情况.从图 3中可以看出SHIN站点形变两种监测结果变化趋势基本一致，整体趋势也相同.所选58个站点情况和SHIN类似，相位相差都不大且均好于其余的点，空间分布也较为均匀，基本可以反映加州地表垂直形变的季节性特征.

2.2 “缩减”系数的计算

计算用GRACE月重力场和GPS监测数据均有中断并非连续时间序列，在对计算结果进行频谱分析时，无论是对原始数据进行插值为等间隔数据分析，还是其他常规的方法都会影响分析结果精度 (安涛等, 2016).本文利用可对非等间隔观测序列进行频谱分析的归一化Lomb_Scargle频谱分析 (Lomb, 1976; Scargle, 1982; Press and Rybicki, 1989).图 4中是站点 (SHIN) L_S形变归一化功率谱分析，(a) 为GRACE结果，(b) 为GPS结果.根据Press等人提出的显著性检测方法绘制出显著性水平为95%所在的红色虚线，两图均有两个峰值在线上，标注出的是周年项，未标注的是半周年项.站点两种形变的功率比为1.105，将此数值认为是真实值与GRACE重力场反演结果的比值.对所选58个GPS站点进行L_S频谱分析，得到GPS观测序列与GRACE重力场负荷形变序列的功率比平均值P_gps/P_grace为1.1127，此即所谓的“缩减系数”.进而可依据此缩减系数恢复重力场位系数处理计算过程中损失的真实信号.

利用GRACE反演地表质量恢复真实信号的问题处理方法有利用GRACE实测数据估计泄露误差理想核函数法 (Wahr et al., 1998)，利用模式数据来估算泄漏误差，Forward-Modeling方法 (Chen et al., 2005) 等.本文为了验证计算的“缩减系数”且为简化计算过程，借助利用模式数据法估算泄漏误差的原理来计算“尺度因子”.由于泄露误差涉及区域内部真是信号泄漏到外部和外部信号泄漏到研究区域内两部分，这两种泄露误差难以辨识分离.在计算尺度因子时又需要考虑到水文分布情况等，所以本文在计算尺度因子时不做过多分类，而是将反演后结果值与真实值得比例记为尺度因子.本文采取的方法是，将计算得到的格网值数据GRID1展开成为球谐系数阶段截取到60阶，然后将此系数进行相同的位系数处理步骤，即P4M6去相关和半径300 km的扇形滤波得到重新计算的格网值GRID2.将区域内的两种格网值平均值的比作为尺度因子，其计算结果为1.2881与之前的缩减系数1.1127较为相近.

3 加州中央山谷陆地水储量变化特征

美国加州地区是全球典型的地下水、地面沉降研究区域之一，特别近些年来加州干旱少雨人们为了保障生活用水和农业灌溉用水而大量开采地下水，致使地区地面沉降加速.在此地区有着丰富的地面监测、空间遥感数据，可以更好地验证反演结果.本文选取与加州地下水资源研究 (Faunt, 2009)、Famiglietti等 (2011)等的研究中同一区域，对加州中央山谷区域水文变化信息进行研究.本文采用基于Noah陆地表面模型的GLDAS水文模式数据分辨率为1°×1°的月数据，取其中地表积雪和土壤湿度数据来代替地表水.基于GRACE反演陆地水储量变化为地表水变化与地下水变化总和的关系，即可得到加州中央山谷地下水动态变化.根据由上得到的缩减系数反算出消除由于位系数滤波处理造影响的陆地水储量变化 (TWS)(图 5蓝色曲线).再用此结果减去同化数据算得的地表水 (GLDAS)(图 5绿色曲线) 得到中央山谷的地下水 (GW) 动态变化 (图 5红色曲线).

从图 5可以看出加州中央山谷区域陆地水储量整体呈下降趋势，地表水变化长期呈现出均衡状态但略有下降，并且变化趋势与陆地水储量变化相近.地下水变化与陆地水变化趋势相似.但是它们的变化相位有所不同，陆地水储量的变化与地表水变化几近一致，只有个别时间段表现出相位滞后或超前 (蓝、绿曲线)，陆地水与地下水之间看不出明显的规律 (蓝、红曲线)，地表水与地下水的变化基本是反相 (绿、红曲线)，这和实际中人们抽取地下水补充生活、农业用水的事实一致.研究的侧重点是加州地下水变化，因此将地下水变化单独拿出来进行分析如图 6所示.

1)图 6显示出加州中央山谷地下水动态变化有几个转折点，2006、2010、2012年，各阶段地下水损失速率不同.从地下水的动态变化基本可以反映出加州近些年来的干旱情况 (Faunt, 2009; Famiglietti et al., 2011; Swain et al., 2014).其中包括2006年加州的大旱和2011年以来至今连续五年的严重旱情.2011年地下水储量有所回升，情况比较反常未能查找到原因.分析其可能有两种因素：一是2010年后加州政府的地下水利用政策和气候变化使得地下水开采得到遏制，造成地下水有所补充；二是GRACE卫星数据精度开始有所降低；三是研究区域水资源得到补充.整体变化情况与陈剑利 (2016)等对加州地下水的研究是一致的，但是在2012年前后也略有差别分析其原因可能也是上述个别因素造成.

2) 加州地下水2003年至2015年损年失率为16.43±1.98 mm/a，2003年至2006年为24.42±16.41 mm/a，2006年至2010年为47.90±7.76 mm/a，2010年至2012年为-50.66±19.98 mm/a，2012年至2015年为57.93±14.08 mm/a.由于第一、三、四时段过短拟合得到的地下水变化精度较差.第一、二时段整体拟合后得到2003年至2010年加州中央山谷地下水年损失率为29.58±4.06 mm/a.其中2006年到2010年加州中央山谷地下水损耗为约为7.3 km³/a，与Scanlon等 (2012)等研究得到的2006年4月至2009年9月地下水损耗率7.7 km³/a基本一致，量值略小的原因可能是本文拟合时段较长所致.

3) 本文结果与2011年Famiglitetti等发表的对加州中央山谷地下水研究成果较为相近.其成果显示2003年10月至2010年3月间地下水损失率为20.4±3.9 mm/a，2006年4月至2010年3月为38.9±9.5 mm/a.上图 6结果均比这个值大9 mm左右，详细的每月结果与其文中结果略有差别，主要表现在2006年和2010年的波谷间距小于其文中图 4.造成这些差别的原因可能由于GRACE反演方法不同，Famiglietti采用的固定窗口滑动去相关 (Swenson and Wahr 2006) 和最佳平均核函数 (Swenson and Wahr, 2002) 的处理方法，另外水文数据基本为实测数据较为详实精度高，本文所采用的陆面同化数据空间分辨率较低.

4 结论 4.1

本文提出的利用附加GPS时序约束的方法来恢复GRACE反演陆地水文信号，从一定的程度上解决了GRACE反演结果的不确定性问题，从实测的角度降低了重力场位系数处理过程中产生的误差，从而在量值上得到较为真实的目标结果.本文算例中利用附加GPS时序约束的GRACE陆地水储量反演方法就是在这种思想的基础上提出并得到验证，在忽略了部分外在信息的条件下得到理想状态水储量的真实值，也为GRACE反演地表质量变化提供了一种思路.

4.2

文中对GPS站点形变的分析忽略了其受季节性、大气压强、温度、潮汐等多种因素影响 (魏娜等, 2015)，另外GPS解算过程中采用相关模型扣除了这些因素，但是这也与GRACE重力场解算中使用模型近似并不完全一致.水文模型过于单一，缺少地表水体变化资料，空间分辨率也较低.水储量变化缺少实测数据和相关加州中央山谷地区水文资料验证.因此，可以考虑选取合适的GPS站点和适宜的处理策略以及和更为适宜的时间序列处理方式来提高精度.例如，选取合适全球IGS站点连续观测数据，解算策略中扣除模型与GRACE重力场解算使用模型相同；考虑温度、站点位置和水文环境等因素；GPS时间序列直接在原始解算时得到；频谱分析可使用小波分析、EMD、窗口Fourier变换等时频分析方法得到时时的频谱做对比.这样将得到的缩减比例应用于大区域和全球范围内的GRACE反演研究中会使得结果更为真实可靠.

致谢 感谢CSR提供的GRACE RL05时变重力场数据，PBO网络GPS数据，以及NASA提供的GLDAS水文模型数据.感谢三位匿名审稿专家及编辑老师对本文提出的宝贵修改意见和建议！