基于重力卫星几何轨道线性化的地球重力场反演方法

陈秋杰，沈云中＊，张兴福

1 同济大学测绘与地理信息学院，上海 200092

2 同济大学空间信息与可持续发展应用中心，上海 200092

3 广东工业大学测绘工程系，广州 510006

1 引 言

随着低轨重力卫星CHAMP、GRACE和GOCE卫星陆续发射及相关数据的释放，国际相关研究机构已使用该类数据反演了多个高精度的地球重力场模型.其主要代表是GFZ（Geo Forschungs Zentrum Potsdam）、CSR（Center for Space Research）、伯尔尼大学、波恩大学所分别反演的 EIGEN［1－2］、GGM［3－4］、AIUB［5］和 ITG－GRACE［6－8］系 列 的 重 力场模型以及多家机构联合（慕尼黑工业大学、波恩大学、格拉茨工业大学、奥地利科学院和伯尔尼大学）所反演的GOCO［9］系列重力场模型，其中，EIGEN和GGM模型均是采用动力学法，ITG－GRACE模型采用短弧长积分法，AIUB模型采用天体力学方法，GOCO模型是联合GRACE和GOCE以及其他重力数据所反演的地球重力场模型.动力学法基于卫星初始状态向量、先验地球重力场模型及其他摄动力模型计算出参考轨道后，以参考轨道为初值对非线性卫星运动方程进行线性化［10－11］，涉及变分方程积分计算，运算量较大.天体力学法同样要用先验地球重力场模型计算参考轨道，它与动力学法的区别在于采用了不同的参数化方法，其卫星状态采用开普勒轨道根数，并增加了经验加速度参数和其他动力学参数［5］.短弧长积分法将弧段的轨道向量表示成边界轨道向量、重力场待估参数和其他力模型参数相关的积分方程［6－8，12］，其线性化的参考轨道由边界轨道参数和已知力模型计算，为减小线性化误差，还可直接用作了梯度改正的几何观测轨道［7］.Ditmar利用三点数值微分公式，根据3个相邻点的卫星轨道观测值计算卫星加速度，并将求得的加速度表示成3点轨道弧段上卫星加速度的加权平均值，建立待估重力场参数与平均加速度之间的关系进行求解［13］，但其解算模型并没有以几何轨道作为初值进行线性化.Xu（2011）给出的模型直接以卫星观测轨道进行线性化［14］，但该模型将线性化轨道参数的改正量重新表示成初始状态参数和重力场参数的函数，这需要进行4重积分.由于GNSS（Global Navigation Satellite System）确定的低轨重力卫星几何轨道精度可达1～2cm，因此直接以几何轨道进行线性化时，其二阶以上的改正量可以忽略.考虑到线性化的一阶改正量完全与定轨精度有关，本文将该改正量作为间接观测值的误差直接由最小二乘求解，除了求解的重力场模型参数外，只需要考虑非保守力加速度计观测值的尺度和偏差参数.与天体力学法和动力学法不同，本文方法不需要卫星的初始状态参数、初始重力场模型和解算变分方程，因此解算模型非常简单、便于程序实现.Mayer－Gürr等（2005）采用短弧长积分法解算ITG－CHAMP01模型时，直接用卫星轨道计算重力场参数的系数矩阵［6］，其本质是忽略轨道误差对系数矩阵影响的一种线性化方案，系数矩阵误差将被观测误差吸收，影响模型的解算精度.因此，Mayer－Gürr（2006）对轨道进行梯度改正［7］，减小其误差对系数矩阵影响后，才用于GRACE高精度的低低卫卫跟踪数据的解算.本文所提出方法本质上是对短弧长积分法的改进，不需要梯度改正就能用GRACE低低卫卫跟踪数据解算地球重力场模型.利用2008年全年的2颗GRACE卫星轨道数据解算了90阶次的地球重力场模型TJGRACE01S，并与其他模型进行比较，结果显示，该模型精度优于同阶次的EIGENCHAMP03S和EIGEN－CHAMP05S模型，其前27阶位系数整体精度甚至优于EIGEN－GRACE01S，其前15阶位系数整体精度与EIGEN－GRACE02S模型精度大致相当.

2 观测与解算模型

对于由m＋1个离散点组成的卫星轨道弧段，其相邻点坐标与加速度间的关系可表示为［15］

其中，r为卫星坐标，ti为第i历元时刻，Δt为积分步长，等于相邻离散点的时间间隔，a（tj）为卫星加速度，等于单位质量作用力，kj为Cowell积分器系数，n为积分器长度，J（i）定义如下：

卫星的加速度a（tj）与卫星位置r（tj）、重力位系数u和加速度计尺度和偏差参数w有关.考虑到

其中，ag（tj）为保守力加速度，af（tj）为非保守力加速度.若直接以卫星轨道观测量rg（tj）以及重力场位系数和加速度计尺度和偏差参数的先验值u0，w0为初值对ag（tj）和af（tj）进行线性化，则

其中，δμ，δw分别为重力位系数参数和加速度计尺度与偏差参数的改正数，vr（tj）为轨道观测值改正数，

顺便指出，由于利用GNSS定轨结果rg（tj）作为初始轨道，引力向量ag（tj）是引力位系数的线性函数，因此系数矩阵H（tj）与先验位系数u0无关，即使没有先验重力场模型（即u0＝0）时，展开式（3）依然成立.将（3）、（4）式代入（2）式，得

其中

将（5）代入（1），则有

对于有m＋1个观测值的轨道弧段，可建立m－n＋1个观测方程，将其表示为

相对于传统动力学法，观测模型（7）的主要特点如下：

（1）将位系数直接在重力卫星观测轨道处展开，而不需表示为卫星初始状态参数及变分方程解算的形式.

（2）由于H（tj）与先验重力场模型的位系数参数无关，因此（7）式是位系数参数的线性观测方程，引入先验重力场模型仅改变常数项，其解算结果不受影响.

如果共有K个弧段，每个弧段k（k＝1，2，…，K）都可以组成如（7）式所示的观测方程：

其中，下标k表示与该弧段相关的量，其他含义与（7）相同.若Qk为观测向量lk的协因数阵，wk为第k弧段加速度计尺度和偏差参数，则第k弧段法方程为

（8）式法方程可分块表示为

位系数改正量δu为全局参数，加速度计尺度和偏差参数改正量δwk为局部参数，Nukuk为位系数参数在该弧段的法方程块，Nwkwk为尺度和偏差参数的法方程块，Nukwk和Nwkuk为全局参数与局部参数间联系数的法方程块，luk和lwk为相应的常数项.考虑到法方程的可叠加性，消去每个弧段的局部参数后，得全局参数的法方程［16］

将每一弧段的全局参数法方程进行叠加，得到求解位系数的总法方程式：

由（11）式可求得位系数参数的改正量.由于本文方法直接以几何观测轨道为初值进行线性化，不需要初始状态向量参数，因此避免了初始状态向量参数与重力场模型参数相关性的影响.但本文的误差方程（7）式的残差向量具有系数矩阵B，导致形成（8）式的法矩阵和常数项时，均需的逆阵，该矩阵的阶数约为弧段历元数的3倍，因此考虑到高维逆矩阵求解的数值稳定性与解算效率，弧段长度不宜设置得过长.同时，采用Cholesky矩阵分解方法进行求逆，可以提高逆矩阵求解的稳定性，相应算法可参考相关文献［17－19］.

3 反演弧段长度

为了选择合适的弧段长度，采用官方JPL（Jet Propulsion Laboratory）所公布的2008年1月份5s采样率的RL02约化动力学轨道数据（RL02版本轨道计算采用GIF48为先验地球重力场模型）和1s采样率的加速度计数据，轨道精度约为2cm，加速度计数据采样率进一步压缩为5s.分别采用不同弧段长度，解算60阶次的地球重力场模型.为减小重力场高阶信号截断对位系数解算精度的影响，第61阶至150阶采用ITG－GRACE2010模型的重力场位系数，要反演的前60阶系数不用先验重力场模型.ITG－GRACE2010模型为波恩大学采用7.5年GRACE低低卫卫跟踪数据反演的静态地球重力场模型.卫星轨道、加速度及姿态观测数据先要进行预处理，包括粗差探测、数据内插、数据间断处理、坐标系转换等.考虑了日月引力、海潮、固体潮、极潮、大气与海洋的非潮汐变化影响等保守力和加速度计测定的非保守力（具体信息见表1），非保守力的尺度和偏差参数先验值采用GFZ公布值，每个弧段估计一组尺度和偏差参数.形成总法方程后采用Cholesky矩阵分解法进行重力位系数的解算，根据解算结果与EGM2008模型的差值，进行相关分析.不同弧段长度求得的位系数与EGM2008模型位系数差值的阶方差和相应的大地水准面累积误差分别如图1和图2所示，其结果表明，30min弧段长度的位系数差值的阶方差和大地水准面累积误差最小，因此实际数据处理时，采用30min弧段长度进行地球重力场反演.

表1 摄动力模型信息Table 1 The information of perturbation force models

4 反演结果与分析

实际数据处理采用 JPL（Jet Propulsion Laboratory）官方公布的2008年全年的RL02版本的GRACE约化动力学轨道和加速度计数据.根据上一节结果，我们采用30min的弧段长度反演90阶次地球重力场模型.其数据预处理方法、保守力与非保守力计算模型如前所述，为减小重力场高阶信号截断对重力场解算精度的影响，第91阶至150阶的重力场位系数采用高精度的ITG－GRACE2010模型系数，前90阶重力位系数的初值设置为零，迭代解算了90阶次的地球重力场模型TJGRACE01S，其大地水准面累积误差为17.6cm.为了分析比较TJGRACE01S模型的精度，本文给出TJGRACE01S、EIGEN－CHAMP03S、EIGEN－CHAMP05S、EIGENGRACE01S和EIGEN－GRACE02S模型与EGM2008模型差值的阶方差如图3所示，相应的大地水准面累积误差如图4所示.

其中，EIGEN－CHAMP03S、EIGEN－CHAMP05S分别为GFZ采用3年、6年的CHAMP卫星轨道解算的地球重力场模型，EIGEN－GRACE01S、EIGENGRACE02S分别为GFZ采用39天、110天的低低卫星跟踪数据解算的地球重力场模型.图3和图4表明，采用一年的GRACE轨道所解算出的重力场模型TJGRACE01S的精度比使用多年CHAMP轨道数据所反演的EIGEN－CHAMP03S和EIGENCHAMP05S模型都要高.从大地水准面累积误差看，前27阶的TJGRACE01S模型整体精度比EIGEN－GRACE01S模型的高；从位系数差值阶方差看，前15阶的TJGRACE01S模型精度与EIGENGRACE02S模型大致相当.显然TJGRACE01S模型的低阶位系数精度比较高，其高阶位系数精度不如 EIGEN－GRACE01S和 EIGEN－GRACE02S，主要原因为GRACE轨道数据适合反演重力场的长波长位系数，而中阶位系数主要依靠GRACE KBR观测数据.

为了进一步分析TJGRACE01S模型的精度与可靠程度，选择TJGRACE01S、EIGEN－CHAMP03S、EIGEN－CHAMP05S、EIGEN－GRACE01S、EIGENGRACE02S和EGM2008模型，分别取不同阶次计算区域：纬度－90°～90°，经度－180°～180°，间隔为1°的网格点的大地水准面高.以EGM2008模型计算得到的大地水准面高为基准值，其他模型求得的结果与基准值之差的统计结果如表2所示.由表2可见，对于不同地球重力场模型，截断至相同阶次（10、30、50、70、90），从所计算的大地水准面高差值的标准差来看，TJGRACE01S模型的精度都比相同阶次的EIGEN－CHAMP03S和EIGEN－CHAMP05S模型的高.

为了检验模型的外符合精度，选择美国8221个A级GPS水准点的实测高程异常作为基准，将EIGEN－CHAMP03S、 EIGEN－CHAMP05S、 EIGENGRACE01S、EIGEN－GRACE02S 和 TJGRACE01S模型所解算的这些点的高程异常与基准值进行比较.为避免各模型的谱泄露，采用谱组合技术，即被分析模型所选择阶次以外的位系数均用EGM2008模型扩充到2160阶次，如被分析模型选择阶次为N，则N＋1阶至2160阶用EGM2008位系数补充，得到各模型的外符合精度如表3所示.表3结果表明，对于不同地球重力场模型，截断至相同阶次（10、30、50、70、90）并与 EGM2008模型进行谱组合，从所计算的模型高程异常与实测高程异常差的标准差来看，除了阶次截断为70外，TJGRACE01S模型所计算的标准差均小于EIGEN－CHAMP03S和EIGEN－CHAMP05S模型，因此可以说TJGRACE01S模型的精度总体高于EIGEN－CHAMP03S和EIGEN－CHAMP05S模型，验证了 TJGRACE01S模型的可靠性以及本文算法的有效性.

表2 不同模型大地水准面高度差统计结果Table 2 The statistical results of geoid height difference between different models

5 结 论

本文提出的基于重力卫星几何轨道线性化的地球重力场反演方法不需要卫星初始轨道参数及先验地球重力场模型和变分方程解算，方法简单，解算效率高.利用2008年全年的GRACE双星卫星轨道数据，采用本文方法解算了90阶次的地球重力场模型TJGRACE01S，结果表明，TJGRACE01S模型在90阶次的大地水准面累积误差为17.6cm，其精度优于同阶次的EIGEN－CHAMP03S和EIGEN－CHAMP05S模型；从大地水准面累积误差看，前27阶位系数精度优于EIGEN－GRACE01S；从位系数差值阶方差看，前15阶位系数精度与EIGEN－GRACE02S模型精度大致相当.利用EGM2008模型和美国8221个GPS水准数据进一步验证了本文方法的有效性.

表3 各模型确定的高程异常与美国实测GPS／水准结果的比较Table 3 The statistic results of the height anomaly differences between GPS／leveling of the USA and the height anomaly from all models

TJGRACE01S模型仅仅采用1年的GRACE双星卫星轨道数据，若采用更长时间的观测数据，并与GRACE KBR数据进行融合处理，TJGRACE01S模型的精度将会更进一步得到提高.

（References）

［1］Reigber C，Schmidt R，Flechtner F，et al.First GFZ GRACE gravity field model IGEN－GRACE01S.http：／／op.gfz－potdam.de／grace／results.

［2］Reigber C，Schmidt R，Flechtner F，et al.An Earth gravity field model complete to degree and order 150from GRACE EIGEN－GRACE02S.Journal of Geodynamics，2005，39（1）：1－10.

［3］Tapley B，Ries J，Bettadpur S，et al.GGM02－An improved Earth gravity field model from GRACE.Journal of Geodesy，2005，79（8）：467－478.

［4］Chambers D P.Evaluation of new GRACE time－variable gravity data over the ocean.Geophys.Res.Lett.，33，L17603，doi：10.1029／2006GL027296.

［5］Beutler G，Jaggi A， Mervart L，et al.The celestial mechanics approach：theoretical foundations.Journal of Geodesy，2010，84（10）：605－624.

［6］Mayer－Gürr T，Ilk KH，Eicker A，et al.ITG－CHAMP01：a CHAMP gravity field model from short kinematic arcs over a one－year observation period.Journal of Geodesy，2005，78（7－8）：462－480.

［7］Mayer－Gürr T.Gravitationsfeldbestimmung aus der Analyse kurzer Bahnboegen am Beispiel der Satellitenmissionen CHAMP und GRACE.Dissertation at the Institute of Theoretical Geodesy，University Bonn，2006.

［8］Mayer－Gürr T，Eicker A，Kurtenbach E，et al.ITGGRACE：Global Static and Temporal Gravity Field Models from GRACE Data.Earth and Environmental Science，2010，（Part 2）：159－168.

［9］Pail R，Bruinsma S，Migliaccio F，et al.First GOCE gravity field models derived by three different approaches.Journal of Geodesy，2011，85（11）：819－843.

［10］沈云中.应用CHAMP卫星星历精化地球重力场模型的研究.武汉：中国科学院测量与地球物理研究所，2000.Shen Y Z.Study of recovering gravitational potential model from the ephemerides of CHAMP（in Chinese）.Wuhan：The Institute of Geodesy and Geophysics，Chinese Academy of Sciences，2000.

［11］周旭华，吴斌，许厚泽等.数值模拟估算低低卫卫跟踪观测技术反演地球重力场的空间分辨率.地球物理学报，2005，48（2）：282－287.Zhou X H，Wu B，Xu H Z，et al.Resolution estimation of Earth gravity field recovery through the low－low satellite to satellite technology by numerical simulation.Chinese J.Geophys.（in Chinese），2005，48（2）：282－287.

［12］Ilk K H， Mayer－Gürr T，Feuchtinger M.Gravity field recovery by analysis of short arcs of CHAMP.Earth Observation with CHAMP，2005：127－132.

［13］Ditmar P，Kuznetsov V，van Eck van der Sluijs A A，et al.＇DEOS＿CHAMP－01C＿70＇：a model of the Earth＇s gravity field computed from accelerations of the CHAMP satellite.Journal of Geodesy，2006，79（10－11）：586－601.

［14］Xu P L. Position and velocity perturbations for the determination of geopotential from space Geodetic measurements.Celestial Mechanics and Dynamical Astronomy，2008，100（3）：231－249.

［15］王正涛.卫星跟踪卫星测量确定地球重力场的理论与方法.武汉：武汉大学，2005.Wang Z T.Theory and methodology of Earth gravity field recovery by satellite to satellite tracking data （in Chinese）.Wuhan：Wuhan University，2005.

［16］张兴福.应用低轨卫星跟踪数据反演地球重力场模型.上海：同济大学，2007.Zhang X F.The earth′s field model recovery on the basis of satellite－to－satellite tracking missions （in Chinese）.Shanghai：Tongji University，2007.

［17］Davis T A，Hager W W.Dynamic supernodes in sparse Cholesky update／downdate and triangular solves.Math Software，2009，35（4）：1－17.

［18］Chen Y，Davis T A，Hager W W，et al.Algorithm 887：CHOLMOD，supernodal sparse Cholesky factorization and update／downdate.Math Software，2009，35（3）：2－14.

［19］Davis T A，Hager W W.Row modifications of a sparse Cholesky factorization.SIAM Journal on Matrix Analysis and Applications，2005，26（3）：621－639.