地球重力场是地球固有的物理特性之一，它反映了地球内部物质分布、运动及其变化的状态，并制约地球自身及其邻近空间中的一切物理事件，因此确定地球重力场历来是大地测量学、地球物理学、海洋学、空间科学以及地球动力学研究的重要课题.准确测定地球重力场及其变化，对于深入研究地球形状与内部构造、防震减灾、探查矿产资源、航天保障和战略武器的发射等诸多领域都具有非常重要的物理意义.

自2002年GRACE 重力卫星发射以来，已能够稳定地提供空间分辨率为几百公里时间分辨率为一个月的地球重力场模型，国内许多学者也进行了卫星重力相关理论和数值解算研究，但目前未见有连续的时变重力场资料公布^[1~6].卫星重力技术被认为是目前中长波长尺度上恢复高精度地球重力场及地球陆地水质量(含冰川)空间分布及变化探测研究中，最为高效、经济和具发展潜力的方法之一.过去，由于资料的周期较短和数据处理方法的限制，数据反演结果并没有达到预计的精度，人们主要利用重力卫星资料研究诸如大尺度的季节信号变化、陆地水储量变化的检测、极地冰川冰盖的消融以及其他地球物理变化等^[7~12].国内许多学者也在大尺度下研究了长江流域水系水储量变化的季节特征^{[13, 14]}.随着资料的积累和数据后处理方法的改进，人们已开始关注GRACE 卫星重力长期变化的意义与解释^[15~17].

本文拟利用现有近7年的GRACE 重力卫星资料，研究广西省红水河上游地区重力场变化，结合全球陆地资料同化模型GLDAS数据计算了土壤湿度对该地区重力场的影响，分析引起该地区重力场变化的可能原因.

本文采用的重力数据由美国德克萨斯大学空间中心CSR 提供，时间为2002年8月到2009年5月GRACE 所测得的共82 个近似月平均时变重力场模型.每个重力场模型对应的数据长度以及重力场之间的时间间隔并不完全一致.潮汐影响，包括海潮、固体潮和地球自转产生的极潮都已经在GRACE 数据处理过程中扣除.非潮汐的大气和海洋影响也在数据处理过程中扣除^[18].这意味着除了大气和海洋模型的误差以外，GRACE 时变重力场反映的是非大气、非海洋的质量变化.在季节性或更短的时间尺度上，主要是陆地水储量的变化.

土壤湿度变化对重力场的影响，本文采用的是2002年8月至2009年5月GLDAS全球陆地资料同化模型数据计算得到.全球陆地资料同化模型GLDAS同化了4个不同的水文模型，并采用了NASA新一代的地面和空间观测技术得到的数据来约束地球表面的状态，进而获得地球表面的近实时信息，是目前最好的全球水文模型之一^[19].为便于和GRACE资料比较，我们对GLDAS 数据进行处理得到每月土壤水的质量变化，并对估算的地表水质量变化做球谐展开，采用与处理GRACE 数据同样的方法进行截断和平滑处理.即在利用展开系数恢复土壤湿度变化对重力场的影响时，仅采用了2~60 阶的展开系数进行计算.其次，由于GLDAS数据计算得到的重力场变化中并不存在条带误差的影响，因此，在数据处理中并没有采用多项式拟合去相关滤波，仅采用了扇形滤波，以保持GLDAS 数据各阶次展开系数与GRACE 模型中各阶次系数相同的权重.

由时变重力场模型计算地球重力异常的公式可表示为

(1)

其中GM为地球引力常数，a为地球赤道平均半径.θ，Φ 分别为余纬和经度，ΔC_nm，ΔS_nm为相对于静止的或时间平均值的位系数的变化量，P_nm为正则化的勒让德函数.由于重力场模型不可能给出从0 到无穷所有阶次的位系数而必须截止于某一阶次，根据(1)式，重力异常的变化也只能计算到模型截止的阶次，这样必然会产生截止误差，以有别于重力场模型系数误差.考虑到模型系数的误差随着阶n的增大而迅速增加，且高阶项系数对重力异常的贡献不可忽略，由(1)式直接计算得到的重力异常的精度对于单点而言比较低，但是对地球物理学者而言，感兴趣的并不是某一点的重力变化，而是某区域内的平均重力变化.即

(2)

其中，区域特征函数u(θ，Φ)在积分区域内为1，在区域外为0，Ω 为积分区域的面积.由于u(θ，Φ)在积分边界不连续，Wahr等提出引入一个平滑函数来代替，即所谓的各向同性的高斯平滑(仅依赖于模型的阶n).之后，许多学者意识到GRACE 数据反演结果中条带形状具有非各向同性，而提出了多种非各向同性的高斯平滑方法，这些方法虽然在其文献中取得了很好的效果，但是由于这些方法的实现过程过于繁杂，不易与其他方法联合处理，因而没有得到广泛应用^{[20, 21]}.扇形滤波的实现相对比较简单，它分别对模型的阶n和次m各应用一次高斯平滑，以达到滤波的效果，特别是对于赤道区域的条带噪声效果显著，其本质上也是一种非各向同性的高斯滤波^[22].扇形滤波可表示为

(3)

其中，w_n，w_m分别为阶n和次m方向上的高斯滤波系数.

在各向同性的高斯滤波引入之后，许多学者研究了不同尺度下(平滑半径从400~1000km)陆地水储量变化的季节信号，从而验证了不同尺度下GRACE 卫星重力的检测能力.为进一步压制条带噪声，提高信噪比，Swenson 等研究发现，GRACE反演结果中的条带误差与模型中某些系数之间存在相关误差有关，进而提出了滑动窗多项式去相关误差滤波方法，该方法在赤道两侧区域取得了显著效果^[23].Swenson等认为，当把重力场模型系数C_nm的次m固定时，偶数阶系数之间以及奇数阶系数之间存在相关误差，提出采用滑动窗多项式最小二乘拟合技术来消除这些系数间的相关误差，达到滤波的效果.但因其在文献中并没有公开具体的实现步骤，因而许多学者按照其思想并未实现文献中的理想效果，目前国际上根据其思想常用的做法是采用多项式拟合去相关技术^{[8, 24]}.

作者根据Swenson的思想，分析了滑动窗多项式去相关误差的特点，针对滑动窗滤波技术的不足，在应用滑动窗去相关误差技术时对数据处理方法做了改进，在赤道两侧区域取得了较好的去条带效果^[25].值得说明的是，虽然目前提出了各种各样的滤波器，但是对于目前GRACE 数据而言，仍没有单独的一种滤波方法能够达到令人满意的滤波效果，必须采用组合滤波的方式才能有效地压制条带噪声提高信噪比.去相关误差方法的优劣，直接决定了高斯平滑核半径选取的大小，如果平滑半径选取过大，则平滑结果中有效信号的位置和形状有可能产生偏差或偏移，其次，相对较弱的信号有可能会被平滑掉.本文数据处理中，滑动窗的宽度为7 点，多项式拟合的次数为3.对应于扇形滤波器阶n和次m方向的滤波半径分别为100km 和300km.

图 1为采用组合滤波方法计算得到的红水河上游区域2002年8月至2009年5月期间重力变化趋势空间分布图.从图中可以看出，在广西省红水河上游地区有一明显的重力上升趋势，范围波及云南省南盘江下游区域、贵州省的北盘江以及广西红水河上游大部分地区，形成以广西天鹅县、贵州黔南布依族苗族自治州罗甸县和北盘江流域为中心的上升区域，中心上升幅度平均每年可达0.6μGal.图 2 为GLDAS模型计算得到的该区域土壤湿度变化对重力场的影响.由于模型系数间不存在相关误差，因而计算中没有采用去相关误差技术.从图中可以看出，红水河上游区域的重力变化趋势并不明显，仅在以昆明和南盘江上游地区为中心有一小的重力上升趋势，上升的幅度平均每年不到0.2μGal，表明昆明和南盘江上游区域的降雨量相对充沛，土壤湿润.而在北盘江至红水河上游区域的重力上升趋势很小，平均每年上升约0.1μGal，说明该地区多年降水基本持平，土壤湿度并无显著变化.那么卫星重力捕捉到该地区的重力上升趋势是由什么原因造成的呢?为此，我们在南盘江下游、北盘江以及红水河上游地区(104.75°E~108.75°E，23.75°N~27.75°N)，将GRACE 资料和GLDAS资料的结果以面积为权重取平均，得到该区域卫星重力变化及地表水变化对重力影响的时间序列，结果如图 3 中实线和虚线所示，实线为卫星重力结果，虚线为土壤湿度变化对重力的影响.从图 3中可以看出，卫星重力变化具有明显的上升趋势，特别是在2006年底之后，有一明显的重力上升，而土壤湿度变化对重力的影响则没有这一现象.作为检验，我们估计了两种数据的年周期信号变化，如图 4 所示，其中实线为卫星重力结果，虚线为GLDAS模型结果.从图中可以看出，两者的年周期信号符合很好，计算结果显示两者振幅差仅为0.18μGal，相位差18.8天，上述结果表明土壤湿度变化对重力的影响更多地表现为该地区重力变化的季节特征，而对该地区重力上升趋势并没有太多的贡献. 为分析该地区重力上升的原因，我们将土壤湿度变化对重力的影响从卫星重力结果中扣除，得到该区域内剩余重力变化的时间序列，如图 5所示.从图中可以发现，剩余重力变化仍然具有上升趋势，其中在2003年下半年重力值出现上升，2004 年年中重力值出现回落，之后重力值维持在相对正常水平，直到2006年10月开始，剩余重力异常出现明显的连续抬升，之后虽有回落，但基本维持在高位.我们以2006年10月为节点，分别对该节点前后的重力变化作平均，结果如图 5中直线所示，计算结果显示2006年10月前后平均重力的变化达2.17μGal.为研究产生这一现象的原因，我们查阅了大量网上资料和相关报道，确认这一现象主要系我国在建第三大水电站———龙滩水电站首次下闸蓄水导致水库及红水河上游区域水储量的增加所致.据中国政府网报道，作为西部大开发标志性工程，龙滩水电站工程位于珠江干流红水河上游的广西省天峨县境内，是国内在建的仅次于长江三峡电站，以及位于金沙江上的溪洛渡电站的特大型水电工程.主体工程于2001年7月1日正式开工，2003年11月6 日实现大江截流，2006年9月30日正式开始下闸蓄水，从起水位215m 蓄至大坝底孔最低泄流高程290.5m(http://www.gov.cn/jrzg/2006-09/30/content_403003.htm).为分析龙滩水库蓄水对周围地区重力场的影响，依据资料，我们采用4km 长、2km 宽、70 m 深的水体模拟计算了水库蓄水对周围地区(104.75°E~108.75°E，23.75°N~27.75°N)重力场的影响.首先将质量异常作球谐展开，之后采用2~60阶的展开系数用来计算储水量变化引起的重力场变化.计算时采用与GRACE 数据处理同样的扇形滤波器，滤波器参数在阶n、次m方向上分别为100km 和300km，同时考虑了负荷形变的影响^[26].计算结果如图 6所示，从图中可以看出，水库储水量变化对周围地区重力的影响与GRACE 结果在空间上非常相似，水库蓄水可引起周围地区平均重力场5.17μGal的重力变化.上述资料和模拟结果表明，该地区的重力上升主要系龙滩水库下闸蓄水引起，而卫星重力能够在中长波长尺度上捕捉到龙滩水库蓄水导致该地区约为2.17μGal重力的变化.

图 1

Fig. 1

图 1 红水河上游地区2002年8月至 2009年5月卫星重力变化趋势 Fig. 1 The satellite gravity changes trend in upstream of Hongshui river area from August 2002 to May 2009

图 2

Fig. 2

图 2 红水河上游地区2002年8月至2009年5月 期间土壤湿度变化对重力变化的影响 图5扣除土壤湿度影响后的剩余重力变化 Fig. 2 The effect of soil moisture changes on gravity field in upstream of Hongshui river area from August 2002 to May 2009

图 3

Fig. 3

图 3 红水河地区重力场变化：GRACE结果（实线）， GLDAS结果（虚线） Fig. 3 Gravity changes in the upstream of Hongshui river area；GRACE data (solid line)，GLDAS data (dotted line)

图 4

Fig. 4

图 4 周年信号变化估计:GRACE结果(实线）， GLDAS结果（虚线） Fig. 4 Change of annual signal； GRACE data (solid line)， GLDAS data (dotted line)

图 5

Fig. 5

图 5 扣除土壤湿度影响后的剩余重力变化 Fig. 5 Residual gravity changes after deducting the effect of soil moisture

图 6

Fig. 6

图 6 水库储水量变化对周围地区重力场的影响 Fig. 6 The effett of water storage variation on gravity field of Longtan Reservoir

广西红水河上游地区的土壤湿度变化主要反映了该地区的季节变化特征，对该地区长期重力上升并没有明显的贡献.这是因为以雷达和微波辐射计资料为主的气候模型，主要反映了浅层地表土壤湿度及地表积雪厚度随时间的变化情况，体现了降雨和蒸发等季节信号特征.

重力场的变化因地因时而异.红水河上游地区土壤湿度变化反映了重力场变化的季节特征，龙滩水电站大坝下闸蓄水引起该地区储水量的增加是该地区重力异常增加的主要原因.这一结果说明GRACE 卫星不仅能够反映大尺度的季节变化信息，而且能够为监测大型水库水储量变化提供近实时的空间对地测量资料.

由于GRACE 卫星的飞行高度和轨道设计以及仪器的灵敏度等特点，决定了GRACE 卫星对重力场的中长波长信息敏感，反映了较大空间尺度上质量随时间的变化情况.随着科技的发展和更加先进的卫星重力计划的实施，我们期待卫星重力技术能够在更小空间和时间尺度上获得重力场变化信息，监测小区域内地球内部物质运移和交换过程.