GRACE (Gravity Recovery And Climate Experiment) 重力卫星探测任务由美国NASA和德国DLR于2002年发起 (Tapley et al., 2004), 迄今为止已经累计了大量的观测数据，三大官方数据处理机构美国德克萨斯大学空间中心CSR (Center for Space Research)、美国宇航局喷气推进实验室JPL (Jet Propulsion Laboratory) 和德国地学研究中心GFZ (GeoFoschungs Zentrum) 利用GRACE观测值和球谐基 (Spherical Harmonic) 函数反演的level 2全球时变重力场产品 (Bettadpur, 2012; Dahle et al., 2014; Watkins and Yuan, 2012) 已经广泛应用在海洋学 (Chambers et al., 2004)、水文学 (Richey et al., 2015)、地震学 (Han and Simons, 2008) 等各个学科 (Kusche et al., 2009).根据ICGEM (International Centre for Global Earth Models) 统计结果，自2000年以来共有超过150个静态或时变重力场模型发布，其中绝大部分是利用球谐基函数建立的模型，包括国内基于GRACE数据解算的时变重力场 (杨帆等, 2017；Chen et al., 2015a, 2015b；王长青等, 2015；冉将军等, 2014).除了常用的球谐函数这种作用域为全球的基函数外，还有另一类局部基函数也逐渐引起关注，因为其在重力场建模中有如下优势：仅需覆盖研究区域而非全球覆盖的数据支持，从而降低计算负荷 (Eicker, 2008; Naeimi, 2013)；可以有效地融合不同分辨率的观测数据，如稀疏的卫星观测值和稠密的地面测站，完成数据融合 (Younis, 2013)；根据已知的地球物理信息作为先验信息进行空间约束求解，从而揭示出更加符合真实情况的信号 (Luthcke et al., 2013; Watkins et al., 2015).因此，发展局部基函数来进一步挖掘GRACE观测数据是十分有意义的，也是当前国际的研究热点.

最早的关于局部重力场恢复的理论可以追溯到Weightman (1967)，其中讨论了利用点质量函数反演地球重力场的方法.除了点质量函数之外，比较重要的局部基函数还包括但不仅限于以下几种：Moritz (1980)用最小二乘配置法研究重力场反演，其中重力场协方差函数被用来作为局域基函数；Mascon局域基函数则最早用于月球重力场恢复，现在也被广泛用于GRACE地球时变重力场解算 (Rowlands et al., 2010; Watkins et al., 2015)，尤其以JPL最新发布的RL05M mascon时变重力场模型为代表，证明了全球无约束mascon解和无约束球谐基函数解等效，但是经过合理的先验信息约束后的mascon算法可以揭示更高分辨率的细节信号；除此之外，国际同行 (Chambodut et al., 2005; Panet et al., 2011; Schmidt et al., 2007) 对小波基函数进行了探讨，并基于GRACE卫星实测或模拟数据，对静态重力场恢复进行了实验；另外，Wieczorek和Simons (2005)提出了一种在局部时频域最优化的Slepian基函数，然而其更多是被应用在GRACE时变重力场的后处理过程如监测地震信号 (Han and Simons, 2008) 和格陵兰冰川消融信号 (Harig and Simons, 2012).从广义的概念上讲，上述局部基函数均可以认为是具有特殊解析表达形式的径向基函数RBF (radial basis function)，这种解析形式的RBF也称为球谐展开到无穷阶 (non-band limited) 的RBF.除此之外，还有另一种狭义的严格数学意义上的RBF，即球谐展开到一定阶次 (band limited) 的RBF，也是本文的研究对象.这种有限带宽的RBF数学性质灵活，可以自由定义各种相关参数，根据其观测数据特征和应用需求设计特定的RBF，而不是仅仅限于上面所提到的几种固定解析形式的RBF.然而，无论是有限带宽的RBF或者解析形式的RBF，因为其参数设置过于多样，利用RBF进行GRACE时变重力场解算还是处于探索阶段，并没有形成统一的理论框架，国际同行对RBF的参数设置进行了大量讨论：Klees等 (2008)基于地面测站数据研究了RBF的格网点深度对重力场恢复的影响；Eicker (2008)利用样条球谐形式的RBF，基于GRACE实测数据反演了2003—2008年的静态重力场，并且研究了格网类型对重力场恢复的影响；Wittwer (2009)基于GRACE从2003到2006年的观测数据，利用相对稀疏的RBF格网分布反演了等同于球谐函数30阶的全球时变重力场，并且研究了格网和深度对重力场恢复的影响；Naeimi (2013)通过模拟GRACE数据研究了观测数据范围和格网范围对局部重力场反演的影响；Bentel等 (2013)通过模拟GRACE数据讨论了具有不同核函数形式的RBF对重力场反演的影响；Weigelt等 (2012)和Antoni (2012)基于GRACE模拟数据，在重力场反演过程中将RBF的核函数形状因子和格网点位置作为待求参数同时反演以获得最优的RBF参数设置.

过去的研究对构建RBF算法的各种参数做了探讨，然而迄今为止，尚无基于GRACE实测观测值反演的等效或者优于目前level 2产品的RBF全球时变重力场模型发布，而国内对于局部基函数反演GRACE时变重力场的相关研究更加匮乏.因此本文对RBF反演时变重力场模型的数学理论进行了讨论，并且基于GRACE L1b实测数据反演了全球无约束RBF模型Hust-IGG03，利用Tikhonov正则化得到全球约束RBF模型Hust-IGG04，同时将以上两模型与90阶球谐函数解GFZ RL05a进行全面细致的对比.需要注意的是，由于构建RBF的相关参数比较复杂，这些参数对重力场恢复的影响已经超出本文的研究范畴，因此本文在GRACE L1b反演中使用了固定的RBF设置，另外，由于局部区域重力场的反演依赖于研究区域的信号质量和其观测数据覆盖范围，为了杜绝选择区域而造成的信号丢失，本文仅利用全球覆盖的GRACE L1b数据反演全球时变重力场模型而非限于某一局部，这样更利于与GFZ RL05a产品进行比较.本文工作将为局部重力场反演的数学理论提供参考，并且为我国将来的重力卫星观测数据和高分辨率地面测站的数据融合提供新的思路.

局部重力场建模要求选用的基函数能量集中在一定范围，RBF因其在空间域和频率域高度局部化的特性，因此常常被选作为局部基函数.下面阐述RBF的构建过程，单位球面上的任意函数F(在本文中F特指重力场)，可以用无穷阶球谐函数的形式展开 (Moritz, 1980)，但在实际应用中往往仅展开到一定阶次求其近似，如下：

(1)

其中，e代表单位球面上任意一点的位置向量，e=(sinθcosλ, sinθsinλ, cosθ); λ, θ分别表示该点的经度和余纬；Y_nm即是常用的球谐函数，f_nm为其相应的展开系数；n, m为球谐函数对应的阶和次，N_max代表在球谐展开实际应用中截断的最高阶次, 本文为了和GFZ RL05a时变重力场模型保持一致，因此N_max=90.球谐函数和完全规格化连带勒让德函数P_nm的数学关系表示如下：

(2)

除了式 (1) 中以球谐基函数表达的球面重力场外，还可以表示成以径向基函数RBF的数学展开形式，如下：

(3)

其中，i代表球面上的某个节点，在利用RBF进行重力场建模之前，必须先将球面按照一定的规律分割成格网，格网上每个节点即是RBF所在位置，I_max是格网点总数，常用的格网类型包括Regular格网，Gauss格网，Discoll-Heally格网，Icosahedron格网，Fibonacci格网和Reuter格网等等.我们列举以上6种格网于图 1, 可以发现，Regular、Gauss和Discoll-Heally这三种格网在高纬度地区的分布过于密集，而在低纬度附近的分布过于分散，这种不均匀性不利于RBF反演重力场，相比之下Icosahedron、Fibonacci、Reuter三种格网分布较均匀，更加适合作为本文RBF算法的几何架构.因此，本文采用了Icosahedron格网 (level 30) 作为RBF的节点分布，其中level直接决定了Icosahedron格网点总数和分布方式 (Eicker, 2008; Sadourny et al., 1968).Φ_i(e_i, e)代表节点i处的球面RBF函数，该函数能量聚集在节点i及其附近，且具有各向同性特征，在任一待估点e的函数值仅与该点和RBF所在节点的球面距离e·e_i有关；a_i是节点i处RBF函数的尺度因子，该因子表征了每个节点RBF对总体重力场贡献的权重.与球谐函数对应的Stokes系数一样，a_i将作为待求参数在GRACE L1b反演过程中求解.RBF的具体表达形式如引言中所述有很多种，本文仅讨论如下比较普遍的有限带宽RBF表达形式：

(4)

图 1

Fig. 1

图 1 构建RBF几何网络结构的格网类型 Fig. 1 Several types of network geometries of RBF global gridding (a) Regular gridding; (b) Gauss Gridding; (c) Discoll-Heally Gridding; (d) Icosahedron Gridding; (e) Reuter Gridding; (f) Fibonacci Gridding.

其中，φ_n作为形状因子表征了RBF的形状；P_n是勒让德多项式.根据球谐函数的加法原理 (Addition theorem)，式 (4) 可以展开为球谐函数形式如下：

(5)

其中，形状因子φ_n基本决定了RBF在频域和空域的特性.为了更加直观地理解形状因子的物理含义，这里引入另一种局部基函数，即球面重力场的协方差函数 (Moritz, 1980)：

(6)

其中，σ_n代表预估重力场的阶方差，可从式 (1) 该预估重力场的球谐展开Stokes系数导出.可以发现，协方差函数C(e_i, e) 和本文所用形式的RBF之间可以通过φ_n=σ_n²/(2n+1) 互相转化，因此在这个意义上RBF的形状因子φ_n表征了待求重力场的阶方差信息.如果对待求重力场的阶方差信息有一定了解，则可以控制RBF的形状因子以得到更加真实的解，这种方式常用在球谐样条类型的RBF中 (Eicker, 2008).除此之外，形状因子的获取也可以通过某种特定的核函数获取，核函数直接确定了形状因子的生成，如Cubic-polynomial, Blackman, Shannon核函数等等 (Naeimi, 2013).在图 2a给出了以上核函数对应的形状因子, 可以发现除了Shannon函数外，其他核函数产生的形状因子对待求重力场的阶方差信息均有一定程度的低通滤波效果，使得其对应的RBF在空域中更加平滑.尽管Shannon核函数在空域的震荡比较大，参见图 2b，但由于其在限定的频域里不含滤波效果，对于本文RBF算法的初步研究更加有利，因此本文仍然选取Shannon核函数来构建RBF的形状因子，该RBF的空间特征和本文所使用格网分布见图 3.

图 2

Fig. 2

图 2 (a) 不同核函数 (Harmonic Spline, Cubic-polynomial, Blackman, Shannon等)2阶到90阶的形状因子, 纵坐标为经过单位正规化后的值；(b) 相邻RBF (由Shannon核函数构建) 在空间的表现，因为RBF各向同性，其值仅和待估点到RBF中心的球面距离 (横坐标) 有关，横坐标单位为度，纵坐标为经过单位正规化后的值 Fig. 2 (a) Several types of kernel functions (Harmonic Spline, Cubic-polynomial, Blackman, Shannon, et al.) up to d/o 90 to produce shape coefficients of RBF, and the Y-axis denotes the values after unit-normalization; (b) Spatial performance of adjacent RBFs in this study (constructed by Shannon kernel function), note that the value of RBF solely depends on the spherical distance between estimated point and central RBF point; unit of X-axis (spherical distance) is degree, and the Y-axis denotes the values after unit-normalization

图 3

Fig. 3

图 3 (a) 单个RBF在空间的表现；(b) 全球RBF的格网分布，由Icosahedron (level 30) 算法生成 Fig. 3 Left panel (a) is the spatial performance of single RBF and the right panel (b) demonstrates the global distribution of our particular RBFs

经过以上设计，我们可以在球面上建立单个RBF的函数模型，但是为了将其与GRACE在轨观测数据联系起来，我们还必须将球面RBF式 (5) 延拓到球外空间任意点r处，如下：

(7)

其中，GM为地球重力常数，R为地球半径.因此将延拓后的RBF式 (7) 代入式 (3) 后，可将球面重力场(e)也延拓为球外重力场(r)：

(8)

利用式 (8) 求出变分方程法所需的所有相关偏导之后，即可利用GRACE观测值建立如下观测方程，经过线性化后表示为

(9)

其中，(Δr, Δ)代表GRACE观测值残差，即双星的轨道和高精度星间距离变率KBRR (k-band range rate) 残差；a代表节点i=1, …, I_max上所有RBF的尺度因子；b代表除RBF尺度因子之外的所有待解参数，如轨道初始位置，加速度计偏差参数等等；v为线性化误差.需要注意的是，由于RBF方法的优势主要体现细节信号上，因此其需要合适的数据处理方案来首先确保总体解算重力场信号水平的正确.关于准确反演GRACE时变重力场的讨论参见作者计算的60阶球谐基函数无约束全球时变重力场Hust-IGG01模型 (杨帆等, 2017).因此，本文所有背景模型 (参见表 2)，待估参数策略，观测值选择等等均与Hust-IGG01模型保持一致，以确保数据处理过程足够准确.除此之外，关于本文RBF反演重力场 (Hust-IGG03, Hust-IGG04) 的策略归纳于表 2，并且与GFZ RL05a球谐函数解法进行比较.

表 1(Table 1)

表 1 背景模型介绍 Table 1 Background models 力模型 备注 描述 静态重力场模型 GIF48 使用2~160阶 (Ries et al., 2011) 固体潮 IERS2010 使用2、3、4阶 (Petit and Luzum, 2010) 极潮 IERS2010 海洋极潮使用2~30阶 (Petit and Luzum, 2010) 海潮 EOT11 使用2~120阶 (Rieser et al., 2012) 大气海洋去混频 AOD1B RL05 handbook 使用2~100阶 (Flechtner et al., 2013) 三体摄动 JPL DE405 仅计算太阳和月亮相关摄动 相对论力 IERS2010 仅考虑太阳和地球因素 (Petit and Luzum, 2010) 非保守力 ACC1B measurements 采样率5 s (Case et al., 2002)

表 1 背景模型介绍 Table 1 Background models

表 2(Table 2)

表 2 本文RBF重力场模型建模方法 Table 2 The RBF modelling in this study and its descriptions 重力场模型 基函数描述 解算方案 后处理/正则化 GFZ RL05a 截止到90阶球谐函数；作用范围全球 Stokes系数作为待解参数，一共8366个；解算范围全球 J1项改正；SLR C20替换；高斯滤波 Hust-IGG03 截止到90阶径向基函数RBF；作用范围局部；Shannon核函数；Icosahedron格网，Level 30 尺度因子作为待解参数，一共9002个；解算范围全球 J1项改正；SLR C20替换；高斯滤波 Hust-IGG04 截止到90阶径向基函数RBF；作用范围局部；Shannon核函数；Icosahedron格网，Level 30 尺度因子作为待解参数，一共9002个；解算范围全球 J1项改正；SLR C20替换；无滤波；Tikhonov正则化

表 2 本文RBF重力场模型建模方法 Table 2 The RBF modelling in this study and its descriptions

通过式 (9) 在每个历元建立观测方程后，利用最小二乘构建法方程即可解算得到所有RBF的尺度因子.由于目前发布的主流重力场模型和已知的GRACE后处理技术都普遍以球谐基函数为基础，因此，为了使本文RBF模型更具普适性并且能直接与当前主流模型进行对比，我们需要将解算的RBF模型再转化成球谐函数表达形式.从RBF重力场建模出发，地球重力场式 (8) 可以改写成

(10)

用球谐基函数表达的球面地球重力场已经在式 (1) 中表达.同样地，我们需要将地球重力场从球面式 (1) 延拓至球外空间任意位置，表达如下:

(11)

因此，式 (10) 和式 (11) 从理论上是等效的.我们对比两式可以发现，球谐函数的系数f_nm与RBF的尺度因子a_i有如下转换关系:

(12)

通过上式，我们可以将任意以RBF表达的重力场再次转换成常用的球谐函数表达形式.需要注意的是，本文所计算的RBF模型Hust-IGG03和Hust-IGG04均是先经过式 (12) 转换成球谐函数表达之后，再进行接下来的对比和评价工作.

利用GRACE level 1b原始数据，基于作者独立开发的重力场反演平台Hawk (Hust orbit and gravity Analysis frame Wor K)，我们成功解算了2009—2010一共24个月的无约束RBF月解时变重力场模型，记为Hust-IGG03.该模型经过系数转换 (式 (12)) 后，等效于90阶球谐函数模型.需要注意的是，采用此“03”命名规则的原因是，作者基于Hawk软件平台，解算了60阶无约束球谐函数SH时变重力场模型Hust-IGG01和Hust-IGG02，并且验证了Hust-IGG01的解算精度与国际同行一致 (杨帆等, 2017).为了检验Hust-IGG03的准确性，我们首先选取任意两个月如2009年5月、2010年10月为例，计算其大地水准面阶方差，并且与GFZ RL05a (90阶) 进行对比.作为参照，本文也对比了Hust-IGG01和CSR RL05 60阶重力场模型.

从图 4a和图 4b可以看出，RBF模型Hust-IGG03与GFZ RL05a 90阶球谐模型基本一致，尤其是在低阶项d/o 20之前，两者基本完全重合；Hust-IGG03 60阶球谐模型与CSR RL05 60阶球谐模型亦有同样表现，在低阶项d/o 20之前基本完全重合.这表明了RBF模型成功探测到了地球质量变化的主要信号 (反映在四个模型的中低阶项里)，并且在细节信号 (反映在高阶项) 的揭示上与GFZ RL05也基本一致.为了更直观地观察RBF模型Hust-IGG03对于质量信号估计的准确程度，本文在图 5中分别展示了2009年5月、2010年10月Hust-IGG03和GFZ RL05a两模型经500km平滑滤波后的空间质量异常分布图.需要注意的是，文中所有空间质量异常分布图均经过如下后处理：(1) J₁项替换 (Swenson et al., 2008); (2) C₂₀项替换 (Cheng and Tapley, 2004); (3) 对于无约束模型进行高斯平滑滤波; (4) 转换成等效水高 (Wahr et al., 2004).

图 4

Fig. 4

图 4 CSR RL05 60阶模型，GFZ RL05a 90阶模型，球谐基函数算法Hust-IGG01 60阶模型，RBF基函数算法Hust-IGG03 90阶模型的大地水准面阶方差 (相对于静态重力场模型GIF48) 对比，单位:m (a) 2009年5月; (b) 2010年10月; (c) 2009年全年平均; (d) 2010年全年平均. Fig. 4 Comparison among monthly gravity products for CSR RL05 up to d/o 60, GFZ RL05a up to d/o 90, SH-based Hust-IGG01 up to d/o 60 and RBF-based Hust-IGG03 up to d/o 90 with respect to static gravity model GIF48, in terms of degree geoid height (m) (a) May 2009; (b) October 2010; (c) Mean of 2009; (d) Mean of 2010.

图 5

Fig. 5

图 5 各重力场模型500 km平滑滤波后的质量异常，以等效水高来表示 (mm) (a) 2009年5月RBF基函数算法Hust-IGG03; (b) 2009年5月球谐基函数算法GFZ RL05a; (c) 2010年10月RBF基函数算法Hust-IGG03; (d) 2010年10月球谐基函数算法GFZ RL05a. Fig. 5 Mass anomaly in terms of equivalent water height (mm) after 500 km smooth filtering (a) RBF-based solution Hust-IGG03 on May 2009; (b) SH-based solution GFZ RL05a on May 2009; (c) RBF-based solution Hust-IGG03 on October 2010; (d) SH-based solution GFZ RL05a on October 2010.

由图 5可以看出，在2009年5月和2010年10月这两个随机选取的月份，Hust-IGG03和GFZ RL05a显示出几乎一样的质量异常信号.这意味着通过RBF或者球谐函数建模得到的无约束时变重力场模型并没有实质型的差别，实际上，类似的现象也在Mascon时变重力场模型中出现 (Rowlands et al., 2010; Watkins et al., 2015).为了进一步验证这个结论，在图 4(c，d) 中导出了2009和2010年Hust-IGG03、GFZ RL05a全年的平均重力场系数的阶方差.同样地，作为参考组，我们也导出了Hust-IGG01、CSR RL05全年的平均重力场系数的阶方差.图 4(c，d) 显示，Hust-IGG03和GFZ RL05a年平均阶方差曲线相关系数在2009年和2010年分别高达0.969、0.967；而Hust-IGG01和CSR RL05阶方差曲线相关系数在这两年也分别高达0.993、0996.本文结果表明Hust-IGG03、GFZ RL05a保持了高度一致，说明基函数的改进或变化并不能直接对无约束重力场信号产生明显影响.根据前人的研究，局部基函数相对于全球基函数 (球谐函数) 的真正优势在于在约束 (正则化) 时变重力场模型中提取到更多细节信号.因此，本文进一步对约束时变重力场进行探讨.

一般情况下，由于向下延拓的不稳定性以及各种其他误差源，GRACE重力场反演是一个典型的病态问题，往往需要利用正则化来克服病态 (Tapley et al., 2004)：

(13)

这里，H是GRACE观测值 (GPS轨道伪观测值，K波段测距或测速观测值) 相对于待解参数(RBF尺度因子或者球谐函数Stokes系数) 的偏微分矩阵，W是不同类型观测值的权矩阵，是待解参数的先验信息，N为包含待解参数先验协方差信息的正则化矩阵，λ为正则化参数用于调节正则化的强度.针对(λ, N, )的定义可以提出各种正则化方案，其中应用最广泛最成功的是Tikhonov正则化 (Tikhonov and Arsenin, 1977)，即要求=0.在这里，我们进一步简化Tikhonov正则化过程，即令正则化矩阵N为单位矩阵.式 (13) 中，正则化参数可利用L-curve或者GCV等数学方法寻找到最优值 (Save, 2009)，但是需要额外进行大量的计算，因此预先为λ确定合适的搜索范围是有必要的，可以显著提高计算效率.为此，我们首先计算了2009年5月的RBF法方程左阵H^TWH，在图 6中将其部分奇异值序列从大到小排列，观察其病态性可以发现条件数大于1×10⁵.根据奇异值的数值范围和计算经验，建议λ的搜索范围为λ=10²⁰~10²²，以较快寻找数学意义上的最优正则化参数，即最优信噪比.本文分别采用了L-curve和GCV两种方法寻找该参数，由于L-curve方法在GRACE解算过程中被认为可能存在过平滑的问题 (Kusche and Klees, 2002)，而本文利用GCV方法亦存在最小极值点不明显的问题，因此采用两种方法联合的做法确定最优参数.

图 6

Fig. 6

图 6 2009年5月解算重力场的法方程左阵的部分奇异值排列 Fig. 6 Singular values of the left hand of normal equations obtained from RBF gravity inversion on May 2009

实际上，更加合理的正则化过程还需融合待解参数的先验协方差信息到正则化矩阵N中，然而球谐基函数的正则化矩阵的设计难度较大，因为其作为一种全球基，在实际工作中将局部先验信息融合到全球基函数的协方差信息中有较大难度.但是局部先验信息的引入在RBF下比较容易实现：RBF作为局域基函数可根据区域信号特征调节正则化矩阵，将RBF基函数网络分成不同的子区域，根据不同区域的先验质量变化信号幅度来近似RBF尺度因子的协方差信息，比如海洋质量变化普遍比陆地水质量变化幅度小而且平滑，故可以给海洋上RBF更强的正则化约束，同时也不会影响陆地区域的信号估计水平.如何结合各种实际水文模型等物理模型以及其他独立观测技术所获得的先验信息来为RBF法方程创建合适的正则化矩阵是本文RBF方法今后发展的重点，但是限于篇幅，本文仅讨论单位正则化矩阵约束RBF时变重力场的解算效果.因此，以无约束RBF时变模型Hust-IGG03为基础，我们结合标准Tikhonov正则化来克服法方程矩阵的病态性，成功获取了2009年至2010年正则化RBF模型Hust-IGG04，并且达到了压制空间条带误差和频域高阶信号误差的效果.我们在图 7中对比了2009年4月正则化RBF模型Hust-IGG04和不同高斯滤波半径处理后的GFZ RL05a结果，同样的在图 8中也对比了2010年9月的结果.

图 7

Fig. 7

图 7 正则化RBF时变模型Hust-IGG04(无滤波) 与球谐基函数时变模型GFZ RL05a (搭配不同半径高斯滤波) 的比较：2009年4月的质量异常，表示为等效水高 (mm) Fig. 7 Comparison between regularized RBF-based monthly gravity model Hust-IGG04 without filtering and SH-based model GFZ RL05a with various Gauss filtering: Mass anomaly in terms of equivalent water height (mm) on April 2009

图 8

Fig. 8

图 8 正则化RBF时变模型Hust-IGG04(无滤波) 与球谐基函数时变模型GFZ RL05a (搭配不同半径高斯滤波) 的比较：2010年9月的质量异常，表示为等效水高 (mm) Fig. 8 Comparison between regularized RBF-based monthly gravity model Hust-IGG04 without filtering and SH-based model GFZ RL05a with various Gauss filtering: Mass anomaly in terms of equivalent water height (mm) on September 2010

在图 7和图 8中，可以发现正则化RBF模型Hust-IGG04未经后处理滤波已经有效地压制了条带误差，并且准确显示出全球主要质量变化信号.与GFZ RL05a不同半径高斯滤波结果相比，本文RBF结果的噪音水平和400~500 km滤波后的GFZ RL05a结果更加接近，但Hust-IGG04更好地恢复了信号的强度.由于GRACE的近极轨设计，高纬度地区的观测数据更加密集，在格陵兰岛和南极洲等高纬度地区，基于正则化RBF模型Hust-IGG04的信号幅度比传统基于全球球谐函数的反演结果明显加强，并且沿着海岸线的信号泄漏现象也明显减弱，Hust-IGG04在这些方面的提高符合已知冰川消融的真实情况 (Velicogna and Wahr, 2013).

考虑到Hust-IGG04和400 km高斯滤波后的GFZ RL05a的噪音水平最为接近，本文在图 9中展示了Hust-IGG04和400 km高斯滤波后的GFZ RL05a在2009—2010两年的质量异常年际趋势图和周年振幅图.从图 9年际趋势图可以看出，Hust-IGG04在大部分的区域质量变化信号均和400 km高斯滤波GFZ RL05a结果一致，在海洋上的噪音也基本处于同一水平，然而Hust-IGG04更加准确地揭示了细节信号，在某些局部区域Hust-IGG04表现出了更加明显的信号幅度增强和分辨率增强，如格陵兰南部区域的质量流失，格陵兰中部地区的质量增加，南极洲半岛附近的质量流失，南美地区和亚洲地区的信号空间分辨率增强等等.另外，从图 9中周年振幅图中我们也可获得类似结论：首先两模型在大部分区域信号一致，但是Hust-IGG04揭示了更多信息，如格陵兰南部地区的周年信号增强，北美的西北海岸线的周年振幅明显增强，亚马逊、南非、孟加拉海岸线等地信号分辨率显著提高.本文从年际趋势图和周年振幅图所观察到的上述地区的信号改善和Watkins等 (2015)用JPL Mascon模型分析所得到的结论类似，如格陵兰南部区域的质量流失趋势和周年振幅的增强，格陵兰中部质量增加趋势的增强，北美的西北海岸线的周年振幅的增强，孟加拉海岸线区域周年振幅的增强等等.

图 9

Fig. 9

图 9 上半部分 (a)(b) 分别表示2009—2010两年期间Hust-IGG04(无滤波) 和GFZ RL05a (400 km高斯滤波) 质量异常的长期趋势 (cm·a-1); 下半部分 (c)(d) 分别表示2009—2010两年期间Hust-IGG04(无滤波) 和GFZ RL05a (400 km高斯滤波) 质量异常的周年振幅 (cm); 红色圈内区域为待研究区域, 参见表 3、4 Fig. 9 The top panel (a)(b) demonstrates the 2009—2010 trend maps for unfiltered Hust-IGG04 and GFZ RL05a with 400 km Gauss filtering, expressed in equivalent water height (cm·a-1); The bottom panel (c)(d) demonstrates the 2009—2010 annual amplitude maps for Hust-IGG04 without filtering and GFZ RL05a with 400 km Gauss filtering, expressed in equivalent water height (cm). Four regions of interest are displayed within the red circles, which are further assessed in Table 3 as well as in Table 4

我们进一步对图 9中四个红线区域的局部信号进行分析，来观察Hust-IGG04模型相对于GFZ RL05a的提高效果.分析结果列在表 3和表 4，其中，区域1为格陵兰主要冰川消融区域 (格陵兰岛南部地区，纬度低于70°)，边界数据来源于Forootan等 (2014)，区域2为南极洲主要冰川消融区域 (Lon [-142°, -77°], Lat [-78°, -72°])，区域3为北美西北海岸区域 (Lon [-151°, -124° ], Lat [54°, 63°])，区域4为澳洲北部区域 (Lon [130°, 146°], Lat [-16°, -12°]).选择区域1、2、3是考虑到JPL mascon模型也在这几处区域有更好的表现，选择区域4是因为澳洲北部有较大的周年振幅，但是利用传统球谐函数反演的官方level 2模型并没有很好捕捉到这里的信号 (Tregoning et al., 2012)，RBF模型有望在这个区域的信号得到提升.根据表 3结果显示，Hust-IGG04相比GFZ RL05a在这四个区域的提高程度分别为24.2%、45.5%、17.1%、3.3%，这些显示RBF更加充分挖掘了GRACE观测值蕴藏的信号.在以上几个区域里，我们对格陵兰岛的冰川消融进行进一步的分析和讨论，在图 10中列出了格陵兰岛南部地区平均等效水高的时间序列，并对比了球谐函数模型和RBF模型在该区域的冰川消融趋势.大量的文献认为由于信号泄漏和相对较低的球谐展开阶数等原因，造成直接从GRACE level 2球谐函数时变重力场产品中所估计的格陵兰岛地区的冰川消融趋势远小于其真实值，因此需要通过后处理信号泄漏修正技术来得到更加准确的估值，但是本文为了实验结果的公平，所有模型均未经过信号泄漏修正.可以发现，GFZ RL05a经400 km高斯滤波和Hust-IGG04两模型在格陵兰南部地区的周年趋势估计值 (以平均等效水高表示) 分别为-17.8 cm·a^-1和-22.1 cm·a^-1.Hust-IGG04在此区域的冰川消融长期趋势较GFZ RL05a提高了24.2%，表明了在球谐基函数重力场产品z在格陵兰区域被低估的信号在RBF时变模型中得到一定程度恢复.同时结合图 9也可以看出，RBF模型中由于正则化而泄漏到海洋区域的质量信号明显低于球谐基函数模型中由于平滑滤波而产生的信号泄漏，从而改善了冰川消融趋势的低估现象.

表 3(Table 3)

表 3 Hust-IGG04和GFZ RL05a区域1和区域2的质量异常年际趋势信号对比 (区域1：格陵兰南部，区域2：南极洲北部) Table 3 Trend assessments on Hust-IGG01 and GFZ RL05a models (Region 1: Southern Greenland, Region 2: Antarctica peninsula) Region Hust-IGG04 (cm·a-1) GFZ RL05a (cm·a-1) Improvement 1 -22.1 -17.8 24.2% 2 -14.4 -9.9 45.5%

表 3 Hust-IGG04和GFZ RL05a区域1和区域2的质量异常年际趋势信号对比 (区域1：格陵兰南部，区域2：南极洲北部) Table 3 Trend assessments on Hust-IGG01 and GFZ RL05a models (Region 1: Southern Greenland, Region 2: Antarctica peninsula)

表 4(Table 4)

表 4 Hust-IGG04和GFZ RL05a区域3和区域4的质量异常周年振幅信号对比 (区域3：北美西北海岸区域 4：澳洲北部) Table 4 Annual amplitude assessments on Hust-IGG01 and GFZ RL05a models (Region 3: the northwestern coast of North America, Region 4: the Northern Australia) Region Hust-IGG04 (cm) GFZ RL05a (cm) Improvement 3 9.6 8.2 17.1% 4 9.4 9.1 3.3%

表 4 Hust-IGG04和GFZ RL05a区域3和区域4的质量异常周年振幅信号对比 (区域3：北美西北海岸区域 4：澳洲北部) Table 4 Annual amplitude assessments on Hust-IGG01 and GFZ RL05a models (Region 3: the northwestern coast of North America, Region 4: the Northern Australia)

图 10

Fig. 10

图 10 (a) 2009—2010年Hust-IGG04(无滤波) 和GFZ RL05a (400 km高斯滤波) 在格陵兰南部地区的平均等效水高时间序列；(b) 两模型时间序列导出的冰川消融年际趋势统计值 Fig. 10 (a) The 2009—2010 time-series of average EWH (equivalent water height) (cm) over southern Greenland derived from Hust-IGG04 (no filtering) and GFZ RL05a (400 km Gauss filtering), respectively; (b) The corresponding statistical of ice-melting yearly trend

另外，我们发现RBF算法 (Hust-IGG04) 不仅在信号幅度上比球谐函数算法有所提高，并且在信号的空间分布上也有所改善.图 11提供了中国以及中国附近区域的质量异常年际趋势分布，我们发现两个模型在这个区域的主要的质量变化信号特征基本一致，但是相比400 km滤波后的GFZ RL05a模型，Hust-IGG04明显表现出整体细节信号的提高，减少了后处理平滑滤波而造成的信号丢失，如Hust-IGG04在中国东南海岸线附近信号泄漏的减少，以及新疆西部质量增加信号的分辨率提高等等.但是，关于信号空间分布提高的验证仍有待下一步利用其他独立观测技术进行具体的量化研究.

图 11

Fig. 11

图 11 2009—2010在中国以及附近地区Hust-IGG04(无滤波) 和GFZ RL05a (400 km高斯滤波) 的质量异常长期趋势，表示为等效水高 (cm·a-1) Fig. 11 The 2009—2010 trend maps over regions around China for unfiltered Hust-IGG04 and for GFZ RL05a with 400 km Gauss filtering, expressed in (cm·a-1) of equivalent water height

局域基函数如mascon和RBF在国际范围已被用于GRACE时变重力场的解算，例如JPL最新发布的RL05M mascon模型，显示出局域基函数在挖掘GRACE数据上比通常使用的球谐基函数更具潜力.然而迄今为止，国内探讨局域基函数反演GRACE时变重力场的文献较少，因此本文对一种RBF特殊局域基径向基函数进行了初步的研究.首先，基于构建RBF的一般理论，我们利用Shannon核函数控制RBF的阶方差信息 (2阶至90阶) 亦即待求重力场的协方差信息，并且利用Icosahedron (level 30) 格网算法生成了本文RBF的几何网络结构；其次，我们详细介绍了Hawk软件系统中GRACE level 1b的数据处理流程以及RBF与球谐函数之间的转换关系；最后，基于以上准备工作，本文利用GRACE观测值成功反演了全球RBF时变重力场模型，结果显示: (1) Hust-IGG03和GFZ RL05a年平均阶方差曲线相关系数在2009年和2010年分别高达0.969、0.967，表明本文无约束RBF模型完全等同于球谐基函数90阶时变重力场模型，与国际上无约束mascon模型的结论一致，即仅仅基函数的改变无法改善目前GRACE重力场的解算精度; (2) 结合标准Tikhonov正则化反演的约束RBF模型Hust-IGG04显示出较高的信噪比，已经无需去相关滤波和高斯平滑等后处理过程，减少了GRACE用户使用重力场产品的复杂程度; (3) 实验表明Hust-IGG04在局部信号的分辨率或信号幅度上明显增强，如在格陵兰南部区域冰川消融趋势估计上Hust-IGG04比GFZ RL05a提高了24.2%.

上述实验结果证明了本文解算的RBF模型相对于传统球谐基函数模型GFZ RL05a对当前GRACE反演是有意义的，但是值得注意的是，本文仅作为RBF解算重力场算法的初探仍不足以揭示RBF全部的潜力，作者认为RBF解算模型结果有望在诸如以下方面得到进一步提高: (1) 针对局部地区重力场解算设计算法，如中国以及附近区域而非全球，这样可以极大降低解算复杂度减少计算时间，并且可以结合当地其他种类的重力观测值精化地区重力场模型; (2) 设计新型RBF形状和格网设置，优化RBF算法; (3) 本文所使用的单位矩阵Tikhonov正则化相对理想，如能结合其他独立观测技术或者已知的物理模型合理添加先验信息于正则化矩阵，将有望进一步提高RBF信噪比.综上所述，我们期望RBF局域基函数算法在未来能成为独立于球谐基函数算法的有效手段，进一步挖掘GRACE或者其他类似观测技术所包含的重力场信息，并且期待结合中国地区特殊的质量变化情况精化中国区域重力场模型.

感谢JPL和GFZ提供的GRACE L1B数据, 以及德国格拉茨大学提供的运动学轨道数据.感谢德国波恩大学APMG提供的超算平台和国家留学基金委CSC对本人在德国所有研究工作的资助.特别感谢中国科学院测量与地球物理研究所冯伟博士和两位匿名评审专家对本文结构和内容提出的宝贵意见.