0 引 言

地球重力场反映了地球物质的空间分布、运动和变化.地球物质系统的质量及其分布是随时间和空间不断变化的，质量的重新分布将导致不同时间尺度的地球重力场变化.稳态的质量分布决定了平均地球重力场信息，非稳态的质量重新分布导致了地球重力场的时变部分，可由时变重力位模型描述(朱广彬，2007).随着卫星重力测量技术的发展，该技术为获取全球范围内的质量迁移和地球动态响应提供了一种有效途径，研究时变地球重力场成为现实.2002年3月GRACE(Gravity Recovery And Climate Experiment)卫星的成功发射开创了高精度全球重力场观测与气候变化试验的新纪元，该计划旨在恢复高精度和高分辨率的全球静态重力场，并估计地球重力场的时变特征，目前已广泛应用于全球及区域陆地水储量变化、地震同震、海洋环流等研究(王武星等，2010；邢乐林等，2012；李琼等，2013).当前，利用GRACE卫星数据反演地球质量变化的方法主要有两种：一种是利用发布的Level-2数据的时变重力场模型，该方法是最常用的方法，解算结果的时间分辨率在30天左右，空间分辨率800 km左右(黄振威，2008)；另一种是利用发布的Level-1B数据的星间距离变率(K-B and Range Rate，KBRR)数据推求地球质量变化，即Mascon(Mass Concentration)方法，采用该方法反演地球质量变化的结果时间分辨率约为10天，空间分辨率约为200 km(张坤，2011)，是当前时变重力场研究的热点之一.与利用时变重力场模型反演地球质量变化的方法相比，Mascon方法在反演精度上具有更高的时间和空间分辨率.国外学者针对Mascon方法的研究取得了许多成果：Rowlands等(2005)提出利用Mascon方法能够反演400 km、 10天分辨率的区域水储量变化；Klosko等(2009)将Mascon方法应用于密西西比流域，验证了和评估了Mascon方法，实验结果表明Mascon方法的反演结果和GLDAS水文模型结果相一致；Krogh等(2010)应用改进的Mascon方法对马拉维湖水储量变化进行了评估，改进的方法中在时空约束条件的条件上增加了区域约束条件，取得了更好的反演结果.由于Mascon方法算法实现上较为复杂，国内相关参考文献也较少.因此，本文针对Mascon方法的数学模型详述了推导过程，分析了Mascon方法的特点，并总结了Mascon方法的解算处理步骤.

1 Mascon方法基本原理

地球质量变化及其重新分布对全球重力场的影响可以通过相对于稳态背景地球重力场模型球谐系数的改变量来描述.将研究区域按一定大小划分为若干个网格，并视同一个网格区域内的质量变化是相同的.t时刻由网格j内质量变化引起的球谐系数的改变量{ΔC_jlm(t)，ΔS_jlm(t)}可以按下面公式计算(Rowlands et al., 2010)：

式中，k_l是l阶负荷LOVE数，σ_j(t)是网格区域j在t时刻单位面积上质量变化，单位是kg/m²，Ω是积分区域且dΩ=sinθdθdλ，l和m分别是地球重力场模型的阶数和次数，θ和λ分别是地心余纬和地心经度，M是地球质量，R是地球半径，P_lm(cosθ)是完全规格化缔合勒让德函数.

另一方面，σ_j(t)可写成σ_j(t)=10×H_j(t).比例因子10 的含义是1 m²区域内质量增加或减少10 kg等效为区域内水柱升高或减少1 cm，H_j(t)为t时刻网格j内质量变化的等效水柱高度，单位是cm，H_j(t)即Mascon参数，且H_j(t)是个阶跃函数，其公式为

由于GRACE卫星对地球表层区域质量变化的时间分辨率 是有限的，故认为在一定时间段内同一网格区域 j的Mascon参数H_j(t)是不变的.因此(1)式又可改写为如下形式：

从(3)式可以看出，对于网格j，该网格区域内质量变化引起的重力场球谐系数改变量与质量变化的等效水柱高度成正比.因此，取Mascon参数H_j(t)=1 cm，令相应质量变化引起的重力场模型球谐系数的改变量作为基础位系数{f^c_jlm，f^s_jlm}，那么Mascon参数H_j(t)可视为相对于这组基础位系数{f^C_jlm，f^S_jlm}的比例因子.

对研究区域内所有网格采取以上相同处理，可得到(5)式，式中N₁为Mascon参数个数.通过上述处理即可将不同区域质量变化对地球重力场的影响分离出来，从而提高反演结果的空间分辨率.

2 观测方程

GRACE卫星的K波段星间测距系统能够精确测定两颗卫星的星间距离变率，精度可达到μm/s.星间距离变率数据KBRR观测值对地球表层质量变化非常敏感，因此可在星间距离变率和Mascon参数之间建立函数关系式，从而解出Mascon参数.根据基于星间距离变率数据确定地球重力场的数学模型的推导(肖云，2006)，可建立星间距离变率观测值与Mascon参数间的关系式为

式中，x₀是两颗卫星的初始状态参数，α是两卫星的加速度计参数，H_i是Mascon参数，ρ · ₀为星间距离变率计算值.由(5)式可知，星间距离变率对Mascon参数H_i的偏导数是星间距离变率对重力场模型球谐系数{C_lm，S_lm}的偏导数的线性组合：

研究结果表明地球表层区域的质量变化对GRACE卫星KBRR观测数据的影响在空间上的分布是有限的(Lemoine et al., 2007)，因此Mascon方法在反演地球表层区域质量变化时仅采用卫星飞临研究区域上空时的KBRR观测数据.故(6)式中N为观测时刻卫星飞临区域相临近的网格个数，而不是研究区域的网格总数.(6)式将Mascon参数与星间距离变率观测值联系起来，为Mascon反演地球表层区域质量变化的观测方程.

为进一步提高Mascon方法反演地球表层区域质量变化的精度，(Rowlands et al., 2010)提出两步法的Mascon方法.具体步骤为：

第一步：处理KBRR和GRACE卫星精密星历数据，估计两颗卫星初始状态参数和加速度计参数，并积分变分方程得到状态转移矩阵和参数敏感矩阵，用于在第二步中组成法方程组.

式中，x₀是两颗卫星的初始状态参数，α是两卫星的加速度计参数，γ是重力场模型球谐系数.该步骤涉及到地球重力场模型球谐系数的求解，但主要目的是为了进一步提高轨道精度.

第二步：将第一步解算得到的两颗卫星的初始状态向量转化为基线参数，并固定第一步解算得到的加速度计参数.第二步中只处理GRACE卫星的星间距离变率KBRR观测值，待估参数为Mascon参数和基线参数P₈、P₁₀和P₁₁，这是因为在解算Mascon参数时同时解算这个三个状态基线参数可以提高Mascon参数的求解精度.

(9)式中，y₀为基线参数P₈、P₁₀和P₁₁，H_i是Mascon参数，N为观测时刻卫星飞临区域相临近的网格个数.

3 时间-空间约束方程

从实际角度考虑，在空间上邻近的网格区域，受相似的环境影响，相邻网格内质量变化应该在数值上是连续的，而不应存在阶跃性的变化.另一方面，对同一个网格区域，邻近时间段内，该网格区域内的质量变化数值上也是趋于一致的，不存在阶跃变化.因此，基于这样考虑，可引入时间-空间约束方程(Rowlands et al.， 2010)，以期待获得更真实的解算结果.时间-空间约束方程形式为

对任意两个Mascon参数H_n和H_k，其中n≠k，满足上述约束方程，可以认为两者在数值上趋于一致的“可能性”满足一个与它们之间的空间和时间“距离”相关的指数函数(张坤，2011).

式中，t_nk是两个Mascon参数在时间上的“距离”，d_nk是两Mascon参数所在网格的空间距离.T和D是表征对Mascon参数的解的约束程度的参数，与研究区域的网格划分有关.一般取T=10 day，D=250到400 km，对应于Mascon方法反演地球表层区域质量变化的时空分辨率(Rowlands et al., 2010).从(11)式可以看出，两个Mascon参数在时空上的“距离”越小，它们在数值上接近的可能性越大；反之，它们二者差异越大.

4 参数解算

将时间-空间约束方程(10)看作虚拟观测方程.在对Mascon参数进行平差解算时，存在两类观测方程：KBRR观测方程和虚拟观测方程.因此，需要对虚拟观测方程赋予合适的权因子P_nk，由第3节分析知，可设置为(Klosko et al., 2009):

式中，S是虚拟观测方程相对于KBRR观测方程取单位权时的比例因子，取值一般是经验性的，也可以采用方差分量估计.对于KBRR观测方程，权值一般设为观测精度平方的倒数，即10¹².

联立(9)式和(10)式，写成误差方程组的形式：

其中，δ X 是待估参数，包括Mascon参数 H和基线参数y₀，A和B分别是观测方程(9)和(10)线性化后的系数矩阵，L是观测值.设P₁和P ₂分别是KBRR观测方程组和虚拟观测方程组的权阵，对(13)式进行平差解算可得：

解线性方程组(14)，即可求得Mascon参数解，解算步骤如图 1.

综上所述，采用Mascon方法反演地球表层区域质量变化时，首先将所研究区域按一定大小划分为若干网格，视每个网格区域内具有相同的质量变化，在星间距离变率观测值和Mascon参数建立关系，组成KBRR观测方程组；并建立时间-空间约束方程，组成虚拟观测方程组；再通过赋予适当权因子，将两类法方程组进行叠加；最后，对叠加后的法方程组进行平差解算，求得Mascon参数的解.

5 结 语

本文详述了利用GRACE卫星的星间距离变率观测值反演地球表层区域质量变化的Mascon方法，并对该方法的数学模型进行了推导分析.Mascon方法的优点在于通过分离不同区域质量变化对重力场的影响，解算时仅采用卫星飞临研究区域上空的观测数据，使得观测数据尽可能少的受到研究区域外质量变化的影响，从而提高了反演结果的时间和空间分辨率.并且，通过增加时间-空间约束方程，可获得更加稳定、符合实际的反演结果.

致 谢 西安测绘研究所翟振和博士的对的研究提出了有益的建议，同时匿名审稿专家也对本文提出了宝贵意见，在此一并致谢！