ABSTRACT

The Gravity Recovery and Climate Experiment (GRACE) mission was launched in 2002. The 15-year GRACE time-variable gravity map tells us how mass transports on or in the Earth with the precision of up to a few centimeters of equivalent water height. The global gravity data have been released in the standard form of spherical harmonic Stokes coefficients as well as in surface Mascon (short for mass concentration) solutions. In this paper, we study the energy integral approach through the estimation of Mascon parameters from GRACE L1B data directly, aiming to improve the spatial and temporal resolutions of local gravity signals. Firstly, we determine the geopotential differences between the twin satellites. Secondly, we recover surface mass anomalies, namely Mascon parameters, from the geopotential difference anomalies in a downward continuation process. We present how to derive in situ geopotential differences and an iterative inversion regularization method to stabilize the Mascon solutions. We employ a closed-loop simulation study over South America to validate the efficiency of this Mascon approach. We demonstrate that by introducing a priori expectation and the associated covariance, the iterative inversion method can effectively reduce striping errors inherent in the measurement configuration and temporal aliasing. The ten-day interval time series of South America revealed a seasonal water cycle and had meaningful variations.

Keywords:GRACE, Mascon, Geopotential Difference, Spatial Constraints, Terrestrial Water Mass Redistributions

利用Mascon方法反演区域地表流体质量变化的研究

俞瑶

武汉大学测绘学院，湖北 武汉

收稿日期：2018年4月6日；录用日期：2018年4月20日；发布日期：2018年4月26日

摘 要

GRACE卫星重力计划提供了高精度的全球时变重力场信息，可应用于研究地表物质迁移。目前主要有两种研究方法，包括利用时变重力场模型反演地表质量变化和利用研究区域的卫星观测数据直接反演局部地表质量变化。本文选择后一种方法，基于能量位差法计算研究区域上方的引力位差，再将位差观测值与局部地表的Mascon参数联系起来并解算病态方程。本文介绍了引力位差计算方法及基于空间约束的病态问题解算方法，利用模拟观测数据，反演了南美洲的水储量变化并对病态问题解算方法进行分析。研究结果表明在解算过程中，引入先验模型的空间约束迭代计算法可以有效提高反演结果精度。模拟计算结果验证了利用能量位差法反演局部地表水储量变化的有效性，2014年南美洲水储量变化结果表明，南美洲水储量呈季节性变化，在10天的时间尺度上也能明显反映出水储量变化特征。

关键词 :GRACE，Mascon，引力位差，空间限制，水储量变化

Copyright © 2018 by author and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

GRACE (Gravity Recovery and Climate Experiment)卫星重力计划由美国宇航局和德国空间飞行中心联合实施，其科学目标是在400到40,000千米的空间尺度上，30天的时间尺度上提供全球重力场信息 [1] ，当前广泛应用于研究全球物理环境和气候变化，如研究全球水循环及区域水储量变化、极地与高山冰川消融、全球海平面变化、固体地球形变等 [2] 。

利用卫星观测数据(Level-1B)直接进行局部求解地表质量变化的优势包括：1) 仅采用研究区域上方的观测值，解算结果的空间分辨率得到了提高；2) 将局部地表质量变化作为待估参数，与球谐系数法相比，参数个数大大减小；3) 解算结果不需要进行滤波处理，高频信息得到保留 [3] ；4) 此外在解算过程中可以加入先验物理模型信息对参数进行约束，或者利用质量变化在空间时间上的相关性进行约束。

将研究区域按一定规则划分为若干格网，每个格网的地表质量变化视为均匀分布的水柱，等效水柱高即为该格网的mascon参数，反映特定格网范围特定时间段内的地表质量变化 [4] 。基于mascon参数利用卫星观测数据直接进行局部求解的方法广义上都可称为mascon方法。可以采用动力法 [2] [4] [5] 、能量位差法 [4] [6] [7] [8] 、在轨视线加速度法 [9] 等将GRACE L1B观测值与mascon参数联系起来，实现局部时变重力场的反演。

本文采用能量位差法反演局部地表质量变化。能量位差法的优点在于引力位与坐标系统选取(地心地固系或者地心惯性系)无关，并且计算过程不需要进行轨道积分 [8] 。利用能量位差法反演局部地表水储量变化主要分为两步，一是计算由水储量变化引起的GRACE两颗卫星间的引力位差，二是建立位差与mascon参数的关系，解算病态方程。本文采用星间距离变率观测数据对星间状态向量校正的方法来计算引力位差；基于牛顿引力定律建立位差与地表水储量变化之间的关系。本文考虑了质量变化的负荷效应影响，并推导了位差与mascon参数之间的线性关系式。由于利用轨道观测数据解算地表质量变化是病态的向下延拓过程，这一过程中，观测数据噪声会被放大，如果不进行约束，得到的结果将不稳定，出现条带现象 [6] [8] 。为获得稳定的结果，本文在解算过程中，引入先验模型进行空间约束以改善法方程求解的病态性。

2. 局部地表水储量变化反演方法

2.1. 基本原理

根据能量守恒原理，在惯性系中卫星引力位V可表示为 [10] ：

$V=E^{kin}+E^{rot}-E^{diss}-C_0=\frac{1}{2}{\left|\dot{r}\right|}^2-{\omega }_e\left(r_x{\dot{r}}_y-r_y{\dot{r}}_x\right)-{\displaystyle \underset{t_0}{\overset{t}{\int }}a}\cdot \dot{r}\text{d}t-C_0$ (1)

其中 ${\omega }_e$ 是地球自转角速度， $r={\left[r_x,r_y,r_z\right]}^{\text{T}}$ 表示卫星位置向量， $\dot{r}={\left[{\dot{r}}_x,{\dot{r}}_y,{\dot{r}}_z\right]}^{\text{T}}$ 表示卫星速度向量， $E^{kin}=\frac{1}{2}{\left|\dot{r}\right|}^2$ 表示动能， $E^{rot}=-{\omega }_e\left(r_x{\dot{r}}_y-r_y{\dot{r}}_x\right)$ 表示地球旋转位能， $E^{diss}={\displaystyle \underset{t_0}{\overset{t}{\int }}a}\cdot \dot{r}\text{d}t$ 表示从 $t_0$ 到t时刻非保守力引起的耗散能， $C_0$ 表示 $t_0$ 时刻的总能量。如果已知卫星运行的位置、速度、所受非保守力，由上式可以计算得到引力位。可将能量常数 $C_0$ 作为待求参数或者通过相邻历元间作双差消去。

引力位V中包含了多种物理信号的影响，可表示为 [11] ：

$V=V^{\text{Earth}}+V^{\text{N-body}}+V^{\text{Tides}}+V^{\text{Ocean}}+V^{\text{Atmosphere}}+V^{\text{Hydrology}}+V^{\text{Others}}+e_v$ (2)

等号右边分别表示参考重力场引力位、N体摄动、潮汐(固体潮、海潮、极潮)、海洋环流、大气、水文、及其他效应(广义相对论效应、冰后回弹、地震、冰雪消融等)对引力位的影响， $e_v$ 是误差项。

由于本文研究的是地表水储量变化，需要扣除引力位V中除 $V^{\text{Hydrology}}$ 外其余信号的影响，由水文效应引起两颗卫星间的引力位差 $V_{AB}^{\text{Hydrology}}$ 可表示为：

$V_{AB}^{\text{Hydrology}}=V_B^{\text{Hydrology}}-V_A^{\text{Hydrology}}$ (3)

将研究区域划分为n个格网，认为所有的格网均在以地球半长轴为半径的球面上，第i个格网面积为 $d{S}_i$ ，格网中点的地心余纬、经度为 ${\theta }_i,{\lambda }_i$ ，在 $\Delta t$ 时间段内该格网表面的水储量变化 $d{m}_i$ 可以用等效水高 $h_i$ 的形式表示为

$d{m}_i\left({\theta }_i,{\lambda }_i,\Delta t\right)={\rho }_{w}d{S}_i{h}_i\left({\theta }_i,{\lambda }_i,\Delta t\right)$ (4)

其中 ${\rho }_w$ 是水密度。根据牛顿引力定律，在地心地固系中，区域地表水储量变化直接引起的GRACE A、B卫星间的引力位差 ${\left\{d{V}_{AB}\right\}}_{\text{surfacemass}}$ 为：

${\left\{d{V}_{AB}\right\}}_{\text{surfacemass}}=G{\displaystyle \sum _{i=1}^n\left(\frac{1}{l_i^B}-\frac{1}{l_i^A}\right)d{m}_i}$ (5)

其中G是万有引力常数， $l_i^A$ 、 $l_i^B$ 分别表示第i个格网中点到A、B卫星之间的距离。可将其球谐展开为：

$\frac{\text{1}}{l_i^{A,B}}=\frac{1}{R}{\displaystyle \sum _{j=0}^{\infty }{\left(\frac{R}{r^{A,B}}\right)}^{j+1}{P}_j\left(\cos {\psi }_i^{A,B}\right)}$ (6)

$\cos {\psi }_i^{A,B}=\cos {\theta }_i\cos {\theta }_{A,B}+\sin {\theta }_i\sin {\theta }_{A,B}\cos \left({\lambda }_i-{\lambda }_{A,B}\right)$ (7)

其中R是地球半长轴， $P_j$ 表示j阶勒让德函数， ${\psi }_i^{A,B}$ 为A、B卫星与第i个格网中点的球面角距， $r^{A,B}$ 、 ${\theta }_{A,B}$ 、 ${\lambda }_{A,B}$ 分别表示A、B卫星的地心距、地心余纬、经度。考虑到地表质量变化产生的负荷效应影响，该效应引起的的引力位差 ${\left\{d{V}_{AB}\right\}}_{\text{solidE}}$ 可表示为 [3] [11] ：

${\left\{d{V}_{AB}\right\}}_{\text{solidE}}=\frac{G}{R}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=2}^{\infty }d{m}_i{k}_j}}\left({\left(\frac{R}{r^B}\right)}^{j+1}{P}_j\left(\cos {\psi }_i^B\right)-{\left(\frac{R}{r^A}\right)}^{j+1}{P}_j\left(\cos {\psi }_i^A\right)\right)$ (8)

其中 $k_j$ 表示j阶负荷勒夫数。那么由区域地表水储量变化引起的卫星间的引力位差 $V_{AB}^{Hydrology}$ 为：

$V_{AB}^{\text{Hydrology}}={\left\{d{V}_{AB}\right\}}_{\text{surfacemass}}+{\left\{d{V}_{AB}\right\}}_{\text{solidE}}$ (9)

令 $H={\left[h_1,h_2,\cdots ,h_n\right]}^{\text{T}}$ 表示 $\Delta t$ 时间段内研究区域的n个mascon参数， $Y={\left[y_1,y_2,\cdots ,y_m\right]}^{\text{T}}$ 为该时段m个历元的 $V_{AB}^{\text{Hydrology}}$ 观测值，将公式(5)、(8)代入公式(9)，得到：

$\Gamma H=Y$ (10)

G是m行n列的设计矩阵，其中：

$\begin{array}{c}{\Gamma }_{p,q}=\frac{G{\rho }_{w}d{S}_q}{R}({\displaystyle \sum _{j=0}^1\left({\left(\frac{R}{r^B}\right)}^{j+1}{P}_j\left(\cos {\psi }_q^B\right)-{\left(\frac{R}{r^A}\right)}^{j+1}{P}_j\left(\cos {\psi }_q^A\right)\right)}\\ +{\displaystyle \sum _{j=2}^{\infty }\left(1+k_j\right)}\left({\left(\frac{R}{r^B}\right)}^{j+1}{P}_j\left(\cos {\psi }_q^B\right)-{\left(\frac{R}{r^A}\right)}^{j+1}{P}_j\left(\cos {\psi }_q^A\right)\right))\end{array}$ (11)

公式(10)即为基于能量位差法的mascon方法反演局部地表水储量变化的基本方程。

2.2. 引力位差的解算

基于能量守恒原理，可以利用高精度的KBRR观测数据对星间状态向量进行校正以提高引力位差计算精度的方法 [11] 。该方法的基本原理如下：

在每个历元中，设置7个待估参数a，分别是三个星间相对位置参数 $r_{AB}$ ，三个星间相对速度参数 ${\dot{r}}_{AB}$ ，以及一个新的参数 $P_{AB}$ ，包括引力位差 $V_{AB}$ ，卫星耗散能之差 $E_{AB}^{diss}$ 和能量常数之差 $C_{AB}$ 之和。将星间距离变率观测值 ${\dot{\rho }}_{obs}$ 、精密轨道数据观测值 $PO{D}_{obs}$ (包括位置差观测值 ${\left(r_{AB}\right)}_{obs}$ 和速度差观测值 ${\left({\dot{r}}_{AB}\right)}_{obs}$ )作为观测值，构建观测方程：

${\dot{\rho }}_{obs}+e_{\rho }=f\left(\alpha \right)={\dot{r}}_{AB}{e}_{AB}$ (12)

$PO{D}_{obs}+e_0=g\left(\alpha \right)=\left[\begin{array}{l}{\left(r_{AB}\right)}_{obs}\\ {\left({\dot{r}}_{AB}\right)}_{obs}\end{array}\right],e_0~N\left(0,C_x\right)$ (13)

其中 $e_{AB}=\frac{r_{AB}}{\Vert r_{AB}\Vert }$ 是卫星连线方向的单位矢量， $e_{\rho }$ 表示KBRR观测误差， $e_0$ 表示轨道观测值误差， $C_x$ 是其协方差阵。将能量守恒方程 $h\left(\alpha \right)=0$ 作为限制条件，代入待求参数初值 ${\alpha }_0$ ，构建误差方程及限制方程，可以迭代计算待求参数改正数 $\Delta \alpha$ 进一步得到星间引力位差。如果GRACE卫星的相对位置精度在2到3 cm，相对速度精度在20 mm/s，星间距离变率测量精度在0.25到0.35 mm/s，不考虑其他噪声的影响，那么计算得到的引力位差精度在2.0 − 3.0 × 10⁻³ m²/s²，与KBRR观测噪声对位差计算精度的影响相当 [11] 。

2.3. 病态问题解算

利用轨道观测数据反演地表质量变化是病态的向下延拓过程。解算过程中观测数据噪声会被放大。为获得稳定的结果，可以引入先验模型及其协方差阵对解算过程进行约束，迭代求解待求参数。先验模型观测方程可表示为：

$x_0=IH+e_0,e_0~\left(0,C_B\right)$ (14)

其中 $x_0$ 是mascon参数的先验值，为n维列向量， $C_B$ 是先验模型协方差阵。观测值方差 ${\sigma }_y^2$ 未知，因此先验模型权阵 $P_B$ 未知，需要迭代求解mascon参数与 ${\sigma }_y^2$ ：

$\widehat{H}={\left({\Gamma }^{\text{T}}\Gamma +{\sigma }_y^2{C}_B^{-1}\right)}^{-1}\left({\Gamma }^{\text{T}}Y+{\sigma }_y^2{C}_B^{-1}{x}_0\right)$ (15)

${\widehat{\sigma }}_y^2=\frac{Y^{\text{T}}Y-{\widehat{H}}^{\text{T}}\left(2{\Gamma }^{\text{T}}Y\right)+{\widehat{H}}^{\text{T}}\left({\Gamma }^{\text{T}}\Gamma \right)\widehat{H}+{\left(x_0-\widehat{H}\right)}^{\text{T}}{C}_B^{-1}\left(x_0-\widehat{H}\right){\sigma }_y^2}{m}$ (16)

利用轨道观测数据反演地表质量变化是病态的向下延拓过程。解算过程中观测数据噪声会被放大。为获得稳定的结果，可以引入先验模型及其协方差阵对解算过程进行约束，迭代求解待求参数。先验模型观测方程可表示为：给定 $x_0,C_x$ 及 ${\sigma }_y^2$ 初值，采用公式(15)和(16)迭代计算mascon参数估值 $\widehat{H}$ 和位差方差估值 ${\widehat{\sigma }}_y^2$ ，直到 ${\widehat{\sigma }}_y^2$ 与上一次计算值之差小于一定限差时停止迭代。如果不能获得先验模型 $x_0$ 及其协方差阵 $C_B$ ，可以将 $x_0$ 初值设为0， $C_B$ 初值设为单位阵，在迭代过程中，将计算得到的 $\widehat{H}$ 作为 $x_0$ ，并计算经验协方差函数 $C_B$ ，方法如下：

根据本次迭代过程中得到的mascon参数估值 $\widehat{H}$ ，计算经验协方差值 ${\widehat{C}}_B$ [12] ：

${\widehat{C}}_B\left({\phi }_{i,j}\right)=\frac{1}{N}{\displaystyle \sum _{i,j}^N{\widehat{h}}_i{\widehat{h}}_j},\phi -\frac{\Delta \phi }{2}\le {\phi }_{i,j}\le \phi +\frac{\Delta \phi }{2}$ (17)

其中j是相关距离， $\Delta \phi$ 是距离间距， ${\phi }_{i,j}$ 是格网 $i,j$ 中点的球面距离，N是符合这样距离范围的格网对数。对计算得到的经验协方差值进行拟合，得到经验协方差函数 $C_B$ [6] ：

$C_B\left(i,j\right)={\sigma }_x^2\exp \left(-\frac{{\phi }_{i,j}}{\beta }\right)$ (18)

其中 ${\sigma }_x^2$ 为经验协方差函数方差， $\beta$ 为相关距离。因此，如果不能获得先验模型协方差阵，那么在迭代过程中，需要先计算mascon参数 $\widehat{H}$ ，利用 $\widehat{H}$ 计算经验协方差函数 $C_B$ ，进而计算 ${\widehat{\sigma }}_y^2$ ，直至收敛。

3. 数值模拟与分析

3.1. 计算方案与流程

本文的计算流程如图1所示。主要步骤包括利用采用全球陆地资料同化系统发布的GLDAS (Global Land Data Assimilation System)全球水文模型数据模拟地表水储量变化，模拟卫星轨道及位差观测数据，选取研究区域，构建观测方程，解算mascon参数，并对结果进行比较分析。

3.2. 数据模拟

本研究采用GLDAS NOAH全球水文模型数据 [13] 模拟地表水储量变化。数据时间分辨率为3小时，空间分辨率为1度 × 1度，采用2014年的土壤水和等效积雪和等效积雪，扣除年均值，将每天8组数据取平均，得到日均水储量变化格网数据。在此基础上将其转换成60阶的球谐系数 [14] 。

将通过上述计算得到的365组日异常重力场模型与60阶的EIGEN-GL04C背景重力场模型叠加，得到模拟的真实重力场模型。不考虑卫星运行所受的非保守力及固体潮、N体摄动等保守力影响，不考虑岁差、章动、极移，假设地球匀速自转，采用14阶的Gauss-Jackson积分器模拟GRACE卫星的轨道。积分步长取5秒，弧长为1天，每天多计算一个历元，作为下一天的初始轨道参数。星间距离变率观测数据根据模拟的轨道数据进行计算。

模拟的卫星轨道和真实轨道相似，卫星轨道高度约为480千米，轨道倾角89度，偏心率0.1，运行

图1. 数据模拟及计算流程

速率约为7600米每秒，两颗卫星间距离约为220千米。模拟的10天轨道观测数据可以实现较为均匀且密集的全球覆盖。考虑到在真实数据处理中存在的各种误差，位差计算精度在0.001 m²/s²量级[7]，本次模拟在水储量变化作用下的位差中加入0.001 m²/s²的白噪声，得到2014年全年的位差观测数据。

3.3. 研究区域选取

本次模拟选择南美洲作为研究区域。南美洲包括世界上最大的流域，亚马逊河流域和拉普拉塔河-巴拉那河流域，流域季节性水储量变化大，且南美洲周围主要是海洋，质量变化相对于陆地来说可忽略不计，在反演过程中质量泄漏效应较小 [2] 。考虑到本次模拟中球谐系数均展开到60阶，对应空间分辨率为330千米，按照4度 × 4度将研究区域均匀划分为221个格网，如图2所示。由于位差观测数据对星下点附近质量变化敏感，选取的轨道范围与研究区域一致 [3] [8] 。由于10天的轨道可以较为均匀且密集的覆盖研究区域，而且位差观测值的采样率为5秒，研究区域上方有近1万个位差观测值，因此本文选择以10天的时间间隔反演南美洲的水储量变化。为验证反演过程的正确性及分析计算结果精度，将10天日均水储量变化格网数据取平均，转换为60阶的球谐系数，在此基础上计算南美洲区域4度格网中点的等效水高作为mascon参数的“真值”。

4. 计算结果与分析

4.1. 病态问题解算方法分析

图3表示2014年2月上旬南美洲水储量变化“真值”以及利用水文效应影响下的引力位差观测值，在不加入任何约束时，直接解得到的水储量变化结果。可以看到由于观测噪声和混频效应的影响，直接

图2. 南美洲格网划分及地形和主要河流示意图[2]

图3. 2014年2月1日至10日南美洲水储量变化模拟计算结果。左图：“真值”；右图：最小二乘直接解结果

解算的结果中有明显的南北方向的条带现象。

利用mascon参数先验值及先验协方差进行空间约束的方法反演得到的水储量变化结果如图4所示。根据计算结果与“真值”之差，可以看到不加先验模型得到的结果中多处区域仍有异常信号，加入先验模型信息后计算结果显著改善，和“真值”相比差别非常小，误差均方根仅有8毫米，这说明引入合适的先验水文模型可以有效提高反演结果精度，抑制观测噪声的影响，削弱条带现象。

4.2. 水储量变化结果分析

上述分析表明，利用能量位差法反演局部地表质量变化的解算过程中，采用引入合适的先验水文模型并迭代求解观测值方差及mascon参数的方法，可以有效的提高反演结果精度。本文基于该方法，采用模拟的2014年全年的水储量变化作用下的位差数据，以60阶的南美洲格网中点10天日均值水储量变化平均值为水文模型先验值，并计算相应的经验协方差函数作为先验协方差阵，以10天时间间隔反演南美洲的水储量变化。2014年南美洲水储量变化反演结果见图5，反演误差统计结果见图6。

南美洲的水储量变化幅度在正负250毫米之间，可以看到亚马逊河、奥里诺科河、托坎廷斯河、巴拉那河、圣弗朗西斯科河流域的水储量有明显的季节性变化。另外，10天的尺度上也能明显反映出水储量变化特征，相邻时间段的水储量变化有较为明显的差异，如12月上旬圣弗朗西斯科河流域无明显质量变化，中旬该区域有明显的质量增加趋势，下旬这一变化更为显著。如果时间分辨率取为一个月，则不能准确反映更短时间内的质量变化，因此提高反演局部地表质量变化的时间分辨率非常必要。

5. 结论

本文介绍了基于能量位差法的mascon方法反演局部地表水储量变化的基本原理，该方法主要分为两步，一是根据轨道观测值计算水文效应影响下的引力位差，二是由建立位差观测值与mascon参数之间的关系，解算病态方程。基于牛顿引力定律建立位差与mascon参数关系的过程中，将卫星与地面格网中点之间的距离进行球谐展开，以便于加入负荷效应改正。本文利用GLDAS水文模型数据模拟2014年全球陆地地表水储量变化，在此基础上模拟了GRACE卫星的轨道数据及水储量变化影响下的位差数据，反演了南美洲的水储量变化，验证了该方法的有效性，并对病态问题解算方法进行分析。在解算过程中采

图4. (a) 无先验信息解算结果；(b) 加入先验信息解算结果；(c) 无先验信息计算结果与“真值”之差；(d) 加入先验信息结果与“真值”之差

图5. 2014年南美洲水储量变化模拟计算结果

图6. 2014年南美洲水储量变化模拟计算结果

用空间约束迭代计算法，并引入合适的先验水文模型，可以有效提高反演结果精度。模拟计算得到的2014年全年南美洲水储量变化结果表明，南美洲水储量呈季节性变化，在10天的尺度上也能明显反映出水储量变化特征。先验水文模型的经验协方差阵参数与水储量变化信号强弱有直接关系。当位差观测值误差在0.001 m²/s²时，采用引入先验模型的空间约束迭代计算法，得到的水储量变化结果的相对误差不超过20% (见图6)。

基金项目

国家自然科学基金项目(41474019)。

文章引用

俞 瑶. 利用Mascon方法反演区域地表流体质量变化的研究
Recovery of Terrestrial Water Mass Redistributions Using Mascon Approach[J]. 地球科学前沿, 2018, 08(02): 331-341. https://doi.org/10.12677/AG.2018.82035

参考文献