地球物理学报  2015, Vol. 58 Issue (9): 3061-3071   PDF    
顾及多方向观测值权比反演地球重力场的动力积分法
罗志才1,2,3, 周浩1, 钟波1,2, 李琼1    
1. 武汉大学测绘学院, 武汉 430079;
2. 武汉大学地球空间环境与大地测量教育部重点实验室, 武汉 430079;
3. 武汉大学测绘遥感信息工程国家重点实验室, 武汉 430079
摘要:考虑到不同坐标系下各个方向观测值对反演地球重力场的频谱贡献不同,建立了顾及多方向观测值权比的动力积分法,并利用该方法反演了高精度的GOCE HL-SST卫星重力场模型.首先,分析了不同坐标系下各个方向观测值与地球重力场信息的响应关系,其中惯性系(IRF)下X、Z方向的观测值分别对扇谐系数、带谐系数最为敏感,Z方向的解算精度在全频段均略高于X、Y方向;地固系(EFRF)下各个方向的独立解算精度均与能量守恒法的解算精度相当;局部指北坐标系(LNOF)下X、ZY三个方向的解算精度依次递减,且Y方向在47阶附近有明显"驼峰"现象.其次,比较了不同坐标系下顾及三个方向观测值权比的加权解算模型,其中加权联合解算模型精度在20至70阶次均明显优于等权解算模型,在带谐项和共振阶次精度提升明显,且LNOF下的加权联合解算精度要优于IRF和EFRF.最后,比较了GOCE和CHAMP卫星的模型解算精度,采用本文计算方法,仅利用2个月GOCE轨道观测值解算的模型精度优于包含更长观测时段信息的AIUB-CHAMP01S和EIGEN-CHAMP03S模型,且略优于ASU-GOCE-2months模型.
关键词动力积分法     加权     GOCE     地球重力场    
Dynamic integral approach based on weighted multi-direction observations for inversion of the earth's gravity field
LUO Zhi-Cai1,2,3, ZHOU Hao1, ZHONG Bo1,2, LI Qiong1    
1. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;
2. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China;
3. State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, Wuhan 430079, China
Abstract: The precise GPS high-low satellite-to-satellite tracking (HL-SST) data is the key supplement to recover the long-wavelength information of the earth's gravity field for the gravity field and steady-state ocean circulation explorer (GOCE) mission. The dynamic integral approach is one of the efficient techniques to determine the spherical harmonic coefficients (SHCs) from GOCE HL-SST observations. However, the traditional dynamic integral approach is based on the inertial reference frame (IRF), and the discrepancy of spectral contribution of multi-direction observations are not fully considered in the SHCs determination.
To analyze the contribution of observations in different directions, we build a new dynamic integral approach which takes the weighted ratios of multi-direction observations into consideration. Meanwhile, a dynamic integral approach at different reference frames is also established, which is used to analyze the response relationship between the observations at different frames. The high precision GOCE HL-SST data are used to evaluate the new dynamic integral approach. In the SHCs determination, the conservative perturbing forces are modeled by precise background models, and the non-conservative perturbing forces are observed by onboard instruments.
Using the new dynamic integral approach, the response relationship between multi-direction observations in different frames and earth gravity information is analyzed firstly. In the IRF, the observations from X and Z directions are most sensitive to sectorial and zonal coefficients, respectively, and the accuracy of the model recovered from Z direction is higher than those recovered from X and Y directions in the whole frequency band. In the earth-fixed reference frame (EFRF), all of the solutions recovered from individual components are of similar accuracy over all harmonic degrees, and approximately equal to the solution recovered from the energy balance approach (EBA). In the local north-oriented frame (LNOF), the accuracy of solutions decreases from X, Z to Y directions, and there is a hump peaking around degrees 47 in the solution recovered solely from Y direction. Secondly, in terms of the weighted solutions which take the different contribution of each component into consideration, their accuracies between 20 to 70 degrees are higher than the equal weighted ones, and the coefficients in the zonal and resonant area are improved obviously. At the north pole and south pole, the gravity anomalies derived from weighted solutions also present better performance than equal weighted solutions. Although the equal weighted solutions have a good consistency in different frames due to the rotation invariant features, the weighted solutions, which adequately take the relationships between all components, have separate accuracy curves, and they are decreased from IRF, EFRF to LNOF. Thirdly, we compare the gravity field solutions based on 2 months of GOCE HL-SST with those obtained from the challenge mini-satellite payload (CHAMP) mission. Because of considering the contribution of observations in different directions, the GOCE HL-SST solution determined with our dynamic integral approach is slightly better than ASU-GOCE-2months model. Although only using 2 months of GOCE orbits, our solutions are much better than AIUB-CHAMP01S and EIGEN-CHAMP03S models, which are determined from 1 year and 2.8 year of CHAMP data, respectively.
From the comparison with equal weighted dynamic integral approach, the numerical results indicate the new dynamic integral approach, which fully considers the contribution of multi-directions, is more suitable for the SHCs determination than the traditional approach. The numerical results also demonstrate that the new dynamic integral approach performance in LNOF is superior to that in IRF and EFRF. Therefore, it is more preferable to adopt the weighted dynamic integral approach in LNOF for SHCs determinations with HL-SST data.
Key words: Dynamic integral approach     Weighting     GOCE     Earth gravity field model    
1 引言

21世纪以来,随着CHAMP、GRACE和GOCE卫星计划的成功实施,人类对地球重力场的认知达到了前所未有的高度.其中,各个重力卫星均搭载了高频GPS接收机,不仅为各个卫星计划的有效载荷提供了高精度的位置信息,又可应用于改善地球重力场中长波信息(Bock et al.,2011Jäggi et al.,2011aBezděk et al.,2014).

基于卫星轨道获取地球重力场模型一直是卫星重力学的研究热点,特别是GOCE卫星成功发射之后,国内外不仅致力于梯度数据确定地球重力场的研究(Pail et al.,2011郑伟等,2011万晓云等,2012),也掀起了基于GOCE高低卫-卫跟踪数据(HL-SST,High-Low Satellite-to-Satellite Tracking)反演地球重力场模型研究的新热潮,研究内容主要包括两个方面:一方面是基于不同方法的研究,另一方面是基于不同卫星任务的比较.Baur等(2012)详细讨论了采用GOCE几何轨道恢复地球重力场模型的加速度法,分析了GOCE卫星获取地球重力场长波信息的能力;Jäggi等(2011b)基于天体力学法反演GOCE HL-SST地球重力场模型,其结果表明几何轨道的协方差信息对恢复地球重力场模型的影响较小;Baur等(2014)基于GOCE几何轨道数据分析了天体力学法、短弧法、平均加速度法、瞬时加速度法和能量守恒法恢复HL-SST地球重力场信息的能力,结果表明除能量守恒法之外的其他方法解算精度均相当;Jäggi et al.,2011aBezděk et al.,2014等分别基于天体力学法、加速度法比较了CHAMP、GRACE和GOCE的反演精度,结果表明GOCE卫星的解算精度优于CHAMP和GRACE.特别地,Visser等(2014)Jäggi等(2015)探讨了GOCE卫星轨道数据恢复时变重力场模型的可行性,也是近年来探寻GRACE卫星任务和GRACE Follow-On卫星计划之间空白时期时变信号获取方法的重要研究内容.此外,钟波等(2013)游为等(2012)Su和Fan(2013)Huang和Fan(2012)黄强等(2013)也分别利用加速度法、短弧积分法和能量守恒法等方法反演了多组GOCE HL-SST模型.

动力积分法虽然计算较为复杂,但其理论严密,具有可同时获取高精度地球重力场模型和高精度动力学轨道等优点(邹贤才,2007).由于GPS卫星观测的固有特性,各个方向观测值的精度各不相同,但在现有研究中,尚未讨论动力积分法中不同坐标系下各个方向观测值与地球重力场信号的响应关系和贡献权比等关键问题.同时,Bruinsma等(2013)基于动力积分法,联合激光测卫、GRACE和GOCE等观测数据获取了GOCE卫星直接解模型,但未对GOCE HL-SST模式下的动力积分法进行详细讨论.除此之外,国内外尚无关于动力积分法获取GOCE HL-SST地球重力场模型的详细研究.鉴于此,本文将建立顾及多方向观测值权比的动力积分法,并以GOCE HL-SST为例,详细讨论不同坐标系下各个方向观测值与地球重力场信息的响应关系,并引入多方向观测值的联合解算公式,比较加权和等权的联合解算模型精度;最后与CHAMP卫星任务恢复的地球重力场模型以及国际上采用GOCE HL-SST数据恢复的最新重力场模型进行比较.

2 动力积分法反演地球重力场模型的基本原理

动力积分法是将卫星星历扰动作为地球重力场信息的泛函,通过建立轨道摄动与地球重力场位系数之间的关系,精密获取地球重力场模型的方法.动力积分轨道的精度不仅受到初始状态的影响,还会受到各种先验力模型参数不确定性的影响,因此在地球重力场模型解算过程中通常一并求解.

2.1 基于HL-SST技术的动力积分反演方法

根据动力积分法的基本原理可知,积分轨道误差主要源于初始状态误差ΔX0和先验力模型参数误差Δp,因此将轨道摄动ΔX(t)表示为二者的线性组合(周旭华,2005肖云,2006邹贤才,2007张兴福,2007游为,2011):

其中,X(t)为卫星的状态向量(x,);(t)为参考轨道,可以基于弧段初始状态X(t0)和先验力模型参数p(重力场位系数、加速度计校准参数等),采用数值积分器积分求得;X(t)X(t0)、X(t)p分别为状态转移矩阵和参数敏感矩阵,可以由数值积分器同时求得.

利用公式(1)建立的观测方程,可通过最小二乘法直接获得位系数修正值、初始状态误差和力模型误差等参数.考虑到力模型误差的累积影响,动力积分法通常分弧段进行:根据公式(2)将分弧段形成的法方程进行约化、叠加、求解(公式(2)和公式(3)中N均为法方程阵分块,W为伴随矩阵分块);利用求得的全局变量Xe修正地球重力场位系数,然后根据公式(3)更新初始状态向量和加速度计校准参数等局部变量Xl,代入下一次循环中,直至获得各类参数的最佳估值.

基于上述方法,即可根据HL-SST技术获取初始状态、加速度校准因子以及位系数等参数的最佳估值,实现加速度计的精密校准和动力学轨道的精密确定,并更新地球重力场模型.

2.2 不同坐标系下的动力积分法

CHAMP、GRACE和GOCE等重力卫星提供的是地固系EFRF(Earth-Fixed Reference Frame)下的卫星几何轨道,而现有大部分学者均是在惯性系IRF(Inertial Reference Frame)下实现的动力积分法.为了分析不同坐标系下各个观测值分量与地球重力场信息的响应关系,本文分别实现了地固系EFRF、惯性系IRF和局部指北坐标系LNOF(Local North-Oriented Frame)下的动力积分法,计算中需要用到的转换关系如公式(4)所示:

其中,Rba表示从a坐标系到b坐标系的转换矩阵,各坐标系之间转换矩阵的具体形式可参见文献(McCarthy and Petit,2004Gruber et al.,2010).若已知a坐标系下的观测值,而要在b坐标系下完成地球重力场反演,则可以基于公式(4)得到b坐标系下的观测轨道,并利用积分器获得b坐标系下的积分轨道和变分方程解(积分中的力模型及其偏导数也均在b坐标系下),然后代入公式(1),即可实现b坐标系下的动力积分法.

2.3 多方向观测值的联合反演

HL-SST技术给出了三个方向的观测值,在不同坐标系下均可依次表示为X、Y、Z三个正交方向.由于天线相位中心变化以及GPS观测的固有特点,低轨卫星的几何轨道精度在各个方向不一致(Bock et al.,2011Jäggi et al.,2009),因此需要单独考虑各个方向观测值对反演结果的影响.首先将公式(2)简化为:

其中,第一项N-1为法方程阵的逆,对应公式(2)右边第一个乘数;W为伴随矩阵,对应公式(2)右边第二个乘数;X为解算的位系数;第二项为间接平差的随机模型,未知数的协因数阵等于N-1D为方差阵,为单位权方差,可由第三项求得;第三项中V为验后轨道残差,n为观测值个数,s为待求解的位系数个数;P为权阵.

利用公式(5)可以获得各个方向观测值独立解算的地球重力场模型,也可以获取各个方向观测值等权联合解算模型.结合各个方向的独立解算模型Xi和验后单位权方差,可以获得加权后的联合解算模型Xc(Yi et al.,2013):

其中,Nc为联合解的法方程阵,是关于各个独立解算模型法方程阵Ni的加权平均;i表示参与联合解算的观测值方向,对应于不同坐标系下的X、Y、Z三个正交方向;Ri为分辨率矩阵(Jackson,1972).

特别地,由公式(6)可知,联合解算模型Xc可表示为关于独立计算模型Xi的加权平均,各个分量权值为Ri.Ri为对角线占优矩阵,可以利用Ri的对角线元素分析各个正交方向观测值对各个联合解算位系数的相对贡献.考虑到:

也可以直接使用Ri的对角线元素之和分析各个分量整体上对联合解算的平均贡献.

3 数值解算与分析

本文基于动力积分法编写了HL-SST模式下的地球重力场反演软件,所有解算过程均采用联合OpenMP和MPI编写的并行函数库(周浩等,20112015),轨道积分和变分方程解算均采用Gauss-Jackson数值积分器(罗志才等,2013),摄动力模型包括三体摄动力、固体潮、海潮、大气潮、极潮和相对论效应,非保守力由加速度计观测数据提供,具体描述见表 1.

表 1 各项摄动力模型及其主要参数 Table 1 Summary of perturbation models and their key elements

解算获取了l阶m次位系数后,均采用公式(8)计算相对于“真实重力场模型” >的阶误差RMS(Root Mean Square)el,以评估反演结果的解算精度(钟波,2010).考虑到ITG-Grace2010s模型为纯GRACE卫星重力场模型,不包含GOCE卫星观测信息且在中低阶次解算精度较高,在后文的精度评定中,均取该模型为“真实重力场模型”.

为分析三个方向观测值的贡献权比以及对联合解的影响,本文将以GOCE HL-SST数据为例进行研究.其中,基于GOCE HL-SST技术恢复地球重力场模型需要Level 2的几何轨道SST_PKI_2、简化动力学轨道SST_PRD_2以及Level 1b的加速度计数据EGG_NOM_1b(Frommknecht et al.,2011Gruber et al.,2010).精密轨道数据由瑞士伯尔尼大学天文研究所(AIUB)提供,精度为2 cm;三个加速度计数据可以由梯度仪3个坐标轴上对称分布的6个加速度计给出.为减小计算量,在不影响解算精度的前提下,将剔除粗差的几何法轨道和加速度计观测数据降采样至5 s.需要说明的是,虽然GOCE卫星引入无阻力控制技术,可以抵消大部分的大气阻力等非保守力的作用,但仍有部分残余,因此在表 1中依然引入了加速度计的校准因子.下文的数值解算与分析均采用上述预处理后的观测数据完成.

3.1 不同坐标系下各方向观测值的解算精度

以2009年12月的观测数据为例,采用不同方向的观测值解算了最大截断阶次为100的地球重力场模型,解算结果如图 1图 2所示.其中,图 1反映了解算模型相对于“真实重力场模型”的位系数误差谱,图 1a、1b和1c分别表示仅采用惯性系下X、Y、Z方向观测值计算的位系数误差谱,图 1d是三方向等权解算的位系数误差谱;图 2图 4给出了在不同坐标系下采用各个观测分量及其联合解的阶误差RMS.

图 1 IRF下各种解的位系数误差谱(a) X方向; (b) Y方向; (c)Z方向; (d)等权联合解.Fig. 1 Coefficient error spectrum in IRF(a)—(c) Solutions solely derived from X, Y, Z component respectively; (d) Solution derived from the combination of three components with equal weight

在人卫轨道研究中,惯性坐标系IRF的XY平面为地球赤道面,X轴指向J2000.0历元的春分点,Z轴指向地球平均旋转轴.在IRF下,卫星可看作是沿着一个椭圆面飞行的,而地球相对于这一轨道面旋转,旋转速度为地球自转平均角速度.因此,若轨道倾角为90°,X、Y方向的轨道摄动完全反映了东西向的重力场信息,而Z方向的轨道摄动完全反映了南北向的重力场信息.已有研究结果表明,GRACE飞行模式对带谐系数部分较为敏感,而SWARM飞行模式对扇谐系数部分较为敏感(Wang,2011Zhou et al.,2014).根据上述讨论,若轨道倾角为90°,理论上IRF下X、Y方向的观测值对扇谐系数敏感,而Z方向的观测值对带谐系数敏感.GOCE卫星轨道倾角为96.7°,相当于绕X轴作6.7°旋转.因此,X方向的观测值主要反映了东西向的重力场信息,对扇谐系数最为敏感,对带谐系数最不敏感,导致扇谐系数解算精度最高,带谐系数解算精度最低(图 1a),反映在图 2所示的阶误差RMS,表现为其低阶次误差最大;Y方向同时反映了东西向和南北向的重力场信息,25阶之前的带谐系数与Z方向的精度相当,60阶之前的扇谐系数与X方向的精度相当(图 1b);Z方向也同时反映了东西向和南北向的重力场信息,但以南北向为主,导致带谐系数精度最高,扇谐部分精度较低,特别是75阶之后误差增加较为明显(图 1c).特别地,Y方向50阶之后的田谐系数误差增加明显,Z方向在50至75阶的近带谐项误差增加明显,具体原因有待进一步研究.综合上述讨论,反映在图 2所示的阶误差RMS上,其结果是:在73阶之前,Y方向的解算精度优于X方向的解算精度;73阶之后,Y方向的解算精度不及X方向的解算精度;Z方向的解算精度在全频段均较高.

图 2 IRF下不同方向的独立解和联合解模型的精度Fig. 2 Accuracy of individual and combined solutions in IRF

GOCE卫星任务提供的轨道数据是地固系EFRF下的.在EFRF下,卫星的飞行轨迹是围着地球的一组均匀分布的包络曲线,其分布特性与轨道重复周期相关(Klokoník et al,2013).由于其均匀分布的特性,基于三个正交方向观测值反演的地球重力场模型精度在全阶次均相近(如图 3所示).此外,能量守恒法(EBA,Energy Balance Approach)是确定地球重力场的具有代表性的方法之一(Jekeli,1999王正涛,2005郑伟等,2006郑伟,2015),由于其解算原理是将轨道位置和速度转换成观测位置处的能量值,损失了部分有用的观测信息.图 3中EFRF(ew)表示地固系下等权联合解算模型,在其基础上乘以得到能量守恒法的理论解算精度EBA(Imaginary).对比计算结果可知,能量守恒法的理论解算精度与EFRF下三个方向的独立解算精度整体上相当,这也表明相比于联合三个正交方向观测值的动力积分法,由于多余观测数的减少,

图 3 EFRF下不同方向的独立解和联合解模型的精度Fig. 3 Accuracy of individual and combined solutions in EFRF

能量守恒法的理论解算精度是其解算精度的倍;同时,能量守恒法的理论解算精度略优于单方向观测值的解算精度,这是因为能量守恒法充分考虑了三个正交方向观测值物理层面的相关性,其解算精度在部分阶次的精度损失小于1/,这与Baur等(2014)的解算结果也是一致的.

局部指北坐标系LNOF的原点为卫星质心,X轴指向北,Z轴沿径向朝外,Y轴垂直于卫星质心的子午面,并与X、Z轴构成右手坐标系.由图 4可知,LNOF下各方向观测值解算的模型精度差异较为明显:X方向在全阶次均最优,Z方向次优,Y方向最差.特别地,Y方向在40阶至58阶和80阶至95阶处的误差异常较为明显,表现为明显的“驼峰”现象,该现象在47阶达到最大值,可认为是轨道共振现象的影响(Klokoník et al.,2013);由于Y方向解算误差的影响,等权模式下的联合解LNOF(ew)也在该部分阶次存在明显“驼峰”现象.

图 4 LNOF下不同方向的独立解和联合解模型的精度Fig. 4 Accuracy of individual and combined solutions in LNOF
3.2 不同坐标系下顾及权比的联合解精度

基于3.1节的讨论可知,不同坐标系下各个方向观测值与地球重力场信息的频谱响应各不相同;同时,考虑到GPS观测的固有特性,本文计算了顾及各个方向观测权比的联合解.首先,基于公式(6)得到了各个方向观测值对应的分辨率矩阵Ri,三者之和与单位阵I之差均达到了10-12量级,可认为符合公式(7)的检核条件;其次,将各个分辨率矩阵的对角线元素相加,可以得到X、Y、Z三个方向观测值对联合解的贡献权比;最后,基于该权比计算了IRF、EFRF和LNOF下的加权联合解,其精度分别见图 2的IRF(optimal)、图 3的EFRF(optimal)和图 4的LNOF(optimal).显然地,采用等权处理方式的联合解精度均优于采用单方向观测值的解算精度,而加权联合解精度均优于等权联合解.

为了更加清晰地比较两种不同处理方式下的解算精度,图 5给出了LNOF下等权和加权联合解相对于“真实模型”的位系数误差谱,对比图 5a图 5b可知,顾及各个方向观测值的权比后,其联合解算模型的带谐项系数明显优于等权模式;由图 5c中二者的差值可知,顾及各个方向观测值的权比也能够提高共振阶次附近的部分位系数精度.此外,图 6给出了相对于极区“真实模型”的重力异常.由于GOCE卫星轨道倾角为96.7°,其星下点轨迹在两极均存在6.7°的极空白,所以利用解算模型计算的该区域重力异常的误差较大(图 6中黑色圆内),即使采用加权处理可提高该部分的精度,为了进一步提高极区重力异常精度,需要引入其他类型观测数据联合求解;而在极空白影响区域外,采用加权处理方式,也能够提高部分区域的重力异常精度.

图 5 LNOF下各种解的位系数误差谱(a)等权;(b)加权;(c)二者之差.Fig. 5 Coefficient error spectrum in LNOF(a)—(b) Solutions derived from the combination of three components with equal weight and different weight respectively; (c) Differences of these two solutions.

图 6 等权解和加权解计算的极区重力异常差异Fig. 6 Difference of polar gravity anomaly computed from equal and different weight

图 7给出了在不同坐标系下,采用等权处理方式得到的联合解算模型之间的精度差异.由于旋转矩阵的不变特性,三种坐标系下的解算精度差值均在10-16量级以下,可认为在等权情况下(各个方向观测值的贡献权比均为33.3%),三种坐标系下的解是等价的.

图 7 不同坐标系下等权解算结果的差异Fig. 7 Comparison of the equal weight solutions derived from different frame

将公式(6)中各个分辨率矩阵Ri的对角线元素相加,可以得到X、Y和Z三个方向观测值对联合解的贡献权比.由于各个坐标系下不同方向观测值与地球重力场的频谱响应不同,基于公式(6)计算的三个方向观测值对联合解的贡献比例各不相同,具体比值见表 2.其中,不论在何种坐标系下,X方向的权最大,均超过了40%;Y和Z方向的权比均小于33.3%,即在基于公式(6)的加权解算中,Y和Z方向的观测值均经过了降权处理.

表 2 不同坐标系下各个分量对联合解的贡献权比 Table 2 Contribution of each component for the combined solutions in different frame

为了便于比较不同坐标系下顾及各个方向观测值权比的联合解精度,图 8给出了三种坐标系下加权联合解的阶误差RMS.不同于等权联合解,加权联合解在不同坐标系下的精度各不相同:在前10阶和70阶次之后,三种坐标系下的联合解与等权解的精度相当;在10至70阶次,EFRF与IRF下的加权联合精度几乎一致,只有在极少数阶次略优;而在10至70阶次,LNOF下的加权联合解精度略优于EFRF和IRF下的联合解精度.

图 8 不同坐标系下加权联合解算结果的差异Fig. 8 Comparison of the weighted solutions derived from different frame

图 2可知,在IRF下仅使用X方向观测值得到的解算模型的阶误差RMS最大,根据误差传播定律,理论上应该减少该方向的权值;而由表 2可知,IRF系下的X方向权比为42.0%,即基于公式(6)的联合处理方式对X方向观测值做了升权处理.为检验本文加权的合理性,在解算中将表 2中IRF系下X方向和Z方向的权比互换,即对X方向做降权处理,对Z方向做升权处理,得到的解算模型精度如图 8中IRF(XZ inverse)所示.采用该联合解算方式得到的精度不但低于加权联合解,在高阶次还不及等权联合解.根据上述讨论,基于公式(6)的联合解,不仅考虑了各个方向独立解的精度,还考虑了各个方向观测值之间的相互关系.因此,在后续研究中建议均采用本文的加权联合解算方式;而在GOCE HL-SST反演中,则建议在LNOF下实现加权联合解算.

3.3 顾及权比的GOCE解精度评定

CHAMP、GRACE和GOCE卫星虽然任务各有不同,但均搭载了高精度GPS接收机,可同时实现卫星精密定轨和地球重力场低频信息的获取.其中,GRACE和CHAMP卫星轨道倾角分别为89.0°和87.3°,可以不考虑极空白的影响,二者均搭载了BlackJack接收机,采用了Chockering天线.考虑到上述相似性,且CHAMP卫星轨道高度(450 km)低于GRACE卫星(500 km),因此本文仅比较了CHAMP与GOCE卫星任务的HL-SST解算模型精度.

AIUB-CHAMP01S和EIGEN-CHAMP03S模型分别采用1年和2.8年的CHAMP观测数据反演,其中AIUB系列采用天体力学法计算,EIGEN系列采用动力积分法计算,且EIGEN-CHAMP03S从70阶次开始做正则化处理,各个模型的阶误差RMS如图 9所示.由于HL-SST模式恢复地球重力场模型的能力随着高度增加而衰减,GOCE卫星的轨道高度(260 km)远低于CHAMP卫星(450 km),因此,在顾及三个方向观测值权比的前提下,仅采用2个月观测值解算的WHU-GOCE-2months模型精度优于AIUB-CHAMP01S和EIGEN-CHAMP03S.

图 9 WHU-GOCE系列模型与其他模型的精度比较Fig. 9 Comparison between WHU-GOCE and other typical solutions

ASU-GOCE-2months模型是Bezděk等(2014)基于加速度法计算的重力场模型,计算中采用了2009年11月至12月共计两个月的GOCE HL-SST数据,截断阶次为75.对比采用同时期观测数据计算得到的WHU-GOCE-2months模型可知,二者在前60阶的阶误差RMS趋于一致;采用本文顾及多方向观测权比的动力积分法,60阶之后的WHU-GOCE-2months模型优于加速度法解算的ASU-GOCE-2months模型.此外,对比分别采用2个月和3个月观测数据得到的WHU-GOCE-2months和WHU-GOCE-3months可知,在观测精度一致的前提下,增加多余观测可有效提升解算精度.根据上述讨论,采用共计4年的GOCE HL-SST观测数据,解算精度有望实现全频段提升.

4 结论

本文在基于HL-SST技术恢复地球重力场基本原理的基础上,给出了不同坐标系下的动力积分法,并引入了多方向观测值的联合解算公式,通过计算多组GOCE HL-SST卫星重力场模型,得到如下结论:

(1)不同坐标系下各个方向观测值与地球重力场信号的响应关系各不相同,IRF下X方向主要反映了东西向的重力场信息,对扇谐系数最为敏感;Y方向同时反映了东西向和南北向的重力场信息;Z方向主要反映了南北向的重力场信息,对带谐系数最为敏感;Z方向的阶误差RMS在全频段优于X、Y方向;EFRF下各个方向的解算精度差别较小,且与能量守恒法的解算精度相当;LNOF系下X方向的阶误差RMS最小,Z方向次之,Y方向最差,且在共振阶次47阶出现最大的“驼峰”现象.

(2)顾及三个方向观测值权比的联合解算方法能够提高解算精度,且在带谐系数和共振阶次部分表现最为明显;在等权联合解算中,不同坐标系下的解算精度相同,而由于顾及了三方向观测值的相关性,LNOF系下的加权联合解精度最高.

(3)采用2个月 GOCE卫星HL-SST数据,基于本文方法计算的WHU-GOCE-2months重力场模型精度与基于加速度法计算的ASU-GOCE-2month模型精度大体相当,且在60阶次之后更优;WHU-GOCE-2months模型精度优于分别采用1年和2.8年CHAMP卫星观测值解算的AIUB-CHAMP01S和EIGEN-CHAMP03S模型.

总体而言,采用本文的顾及多方向观测值权比的动力积分法解算精度较高,建议在后续的研究中采用本文的处理方法.目前,对不同精度观测值的加权处理方式较多,且方法之间的数学原理各不相同,同时局限于科研人员对卫星重力数据的数学物理特性的认知,针对卫星重力数据的加权处理方式仍需要进一步展开.

致谢 感谢ESA提供解算所需的GOCE卫星观测数据,感谢Bezděk博士对重力场解算过程中相关问题的帮助,感谢两位匿名审稿专家提出的宝贵修改意见.

参考文献
[1] Baur O, Reubelt T, Weigelt M, et al. 2012. GOCE orbit analysis: Long-wavelength gravity field determination using the acceleration approach. Advances in Space Research, 50(3): 385-396.
[2] Baur O, Bock H, Hck E, et al. 2014. Comparison of GOCE-GPS gravity fields derived by different approaches. Journal of Geodesy, 88(10): 959-973.
[3] Bezděk A, Sebera J, Klokoník J, et al. 2014. Gravity field models from kinematic orbits of CHAMP, GRACE and GOCE satellites. Advances in Space Research, 53(3): 412-429.
[4] Bock H, Jäggi A, Meyer U, et al. 2011. GPS-derived orbits for the GOCE satellite. Journal of Geodesy, 85(11): 807-818.
[5] Bruinsma S L, Frste C, Abrikosov O, et al. 2013. The new ESA satellite-only gravity field model via the direct approach. Geophysical Research Letters, 40(14): 3607-3612.
[6] Frommknecht B, Lamarre D, Meloni M, et al. 2011. GOCE level 1b data processing. Journal of Geodesy, 85(11): 759-775.
[7] Gruber T, Rummel R, Abrikosov O, et al. 2010. GOCE level 2 product data handbook. GO-MA-HPF-GS-0110.
[8] Huang Q, Fan D M, You W. 2013. Recovery of earth′s gravitational field model based on GOCE satellite orbits. Geomatics and Information Science of Wuhan University (in Chinese), 38(8): 907-910.
[9] Huang Q, Fan D M. 2012. The average acceleration approach applied to gravity coefficients recovery based on GOCE orbits. Geodesy and Geodynamics, 3(4): 18-22.
[10] Jackson D D. 1972. Interpretation of inaccurate, insufficient and inconsistent data. Geophysical Journal International, 28(2): 97-109.
[11] Jäggi A, Dach R, Montenbruck O, et al. 2009. Phase center modeling for LEO GPS receiver antennas and its impact on precise orbit determination. Journal of Geodesy, 83(12): 1145-1162.
[12] Jäggi A, Bock H, Prange L, et al. 2011a. GPS-only gravity field recovery with GOCE, CHAMP, and GRACE. Advances in Space Research, 47(6): 1020-1028.
[13] Jäggi A, Prange L, Hugentobler U. 2011b. Impact of covariance information of kinematic positions on orbit reconstruction and gravity field recovery. Advances in Space Research, 47(9): 1472-1479.
[14] Jäggi A, Bock H, Meyer U, et al. 2015. GOCE: assessment of GPS-only gravity field determination. Journal of Geodesy, 89(1): 33-48.
[15] Jekeli C. 1999. The determination of gravitational potential differences from satellite-to-satellite tracking. Celestial Mechanics and Dynamical Astronomy, 75(2): 85-101.
[16] Klokoník J, Gooding R H, Wagner C A, et al. 2013. The use of resonant orbits in satellite geodesy: a review. Surveys in Geophysics, 34(1): 43-72.
[17] Luo Z C, Zhou H, Zhong B, et al. 2013. Analysis and validation of Gauss-Jackson integral algorithm. Geomatics and Information Science of Wuhan University (in Chinese), 38(11): 1364-1368.
[18] McCarthy D D, Petit G. 2004. IERS conventions (2003). International Earth Rotation and Reference Systems Service, GERMANY.
[19] Pail R, Bruinsma S, Migliaccio F, et al. 2011. First GOCE gravity field models derived by three different approaches. Journal of Geodesy, 85(11): 819-843.
[20] Su Y, Fan D M. 2013. Gravity field recovery from GOCE orbits using the energy conservation approach. Geodesy and Geodynamics, 4(2): 40-46.
[21] Visser P N A M, van der Wal W, Schrama E J O, et al. 2014. Assessment of observing time-variable gravity from GOCE GPS and accelerometer observations. Journal of Geodesy, 88(11): 1029-1046.
[22] Wan X Y, Yu J H, Zeng Y Y, et al. 2012. Frequency analysis and filtering processing of gravity gradients data from GOCE. Chinese J. Geophys. (in Chinese), 55(9): 2909-2916, doi: 10.6038/j.issn.0001-5733.2012.09.010 .
[23] Wang X X. 2011. Time-variable gravity field determination from satellite constellations [Ph. D. thesis]. Technische Universitt München.
[24] Wang Z T. 2005. Theory and methodology of earth gravity field recovery by satellite-to-satellite tracking data [Ph. D. thesis] (in Chinese). Wuhan: Wuhan University.
[25] Xiao Y. 2006. Research on the earth gravity field recovery from satellite-to-satellite tracking data [Ph. D. thesis] (in Chinese). Zhengzhou: PLA Information Engineering University.
[26] Yi W Y, Rummel R, Gruber T. 2013. Gravity field contribution analysis of GOCE gravitational gradient components. Studia Geophysica et Geodaetica, 57(2): 174-202.
[27] You W. 2011. Theory and methodology of earth′s gravitational field model recovery by LEO data [Ph. D. thesis] (in Chinese). Chengdu: Southwest Jiaotong University.
[28] You W, Fan D M, He Q B. 2012. Recovering earth′s gravitational field model using GOCE satellite orbits. Geomatics and Information Science of Wuhan University (in Chinese), 37(3): 294-297.
[29] Zhang X F. 2007. The earth′s field model recovery on the basis of satellite-to-satellite tracking missions [Ph. D. thesis] (in Chinese). Shanghai: Tongji University.
[30] Zheng W, Shao C G, Luo J, et al. 2006. Numerical simulation of Earth′s gravitational field recovery from SST based on the energy conservation principle. Chinese J. Geophys.(in Chinese), 49(3): 712-717.
[31] Zheng W, Hsu H T, Zhong M, et al. 2011. Accurate and rapid determination of GOCE Earth′s gravitational field using time-space domain method associated with Kaula regularization. Chinese J. Geophys. (in Chinese), 54(1): 14-21, doi: 10.3969/j.issn.0001-5733.2011.01.003.
[32] Zheng W. 2015. Theory and Approach of Satellite Gravity Inversion on the Basis of Energy Balance Principle (in Chinese). Beijing: Science Press.
[33] Zhong B. 2010. Study on the determination of the earth′s gravity field from satellite gravimetry mission GOCE [Ph. D. thesis] (in Chinese). Wuhan: Wuhan University.
[34] Zhong B, Luo Z C, Zhou H. 2013. Gravity field recovery using satellite average accelerations derived from high-low satellite-to-satellite tracking. Geomatics and Information Science of Wuhan University (in Chinese), 38(8): 902-906.
[35] Zhou H, Zhong B, Luo Z C, et al. 2011. Application of parallel algorithms based on OpenMP to satellite gravity field recovery. Journal of Geodesy and Geodynamics (in Chinese), 31(5): 123-127.
[36] Zhou H, Luo Z C, Zhong B, et al. 2015. MPI parallel algorithm in satellite gravity field model inversion on the basis of least square method. Acta Geodaetica et Cartographica Sinica (in Chinese), 44(8): 833-839.
[37] Zhou H, Luo Z C, Zhong B, et al. 2014. A simulation study of Earth gravity field model inversion from SWARM mission. EGU General Assembly Conference Abstracts, 16: 5888.
[38] Zhou X H. 2005. Study on satellite gravity and its application [Ph. D. thesis] (in Chinese). Wuhan: Institute of Geodesy and Geophysics, Chinese Academy of Sciences.
[39] Zou X C. 2007. Theory of satellite orbit and earth gravity field determination [Ph. D. thesis] (in Chinese). Wuhan: Wuhan University.
[40] 黄强, 范东明, 游为. 2013. 利用GOCE卫星轨道数据反演地球重力场模型. 武汉大学学报·信息科学版, 38(8): 907-910.
[41] 罗志才, 周浩, 钟波等. 2013. Gauss-Jackson积分器算法分析与验证. 武汉大学学报·信息科学版, 38(11): 1364-1368.
[42] 万晓云, 于锦海, 曾艳艳. 2012. GOCE引力梯度的频谱分析及滤波. 地球物理学报, 55(9): 2909-2916, doi: 10.6038/j.issn.0001-5733.2012.09.010.
[43] 王正涛. 2005. 卫星跟踪卫星测量确定地球重力场的理论与方法[博士论文]. 武汉: 武汉大学.
[44] 肖云. 2006. 基于卫星跟踪卫星数据恢复地球重力场的研究[博士论文]. 郑州: 中国人民解放军信息工程大学.
[45] 游为. 2011. 应用低轨卫星数据反演地球重力场模型的理论和方法[博士论文]. 成都: 西南交通大学.
[46] 游为, 范东明, 贺全兵. 2012. 利用GOCE卫星轨道反演地球重力场模型. 武汉大学学报·信息科学版, 37(3): 294-297.
[47] 张兴福. 2007. 应用低轨卫星跟踪数据反演地球重力场模型[博士论文]. 上海: 同济大学.
[48] 郑伟, 邵成刚, 罗俊等. 2006. 基于卫-卫跟踪观测技术利用能量守恒法恢复地球重力场的数值模拟研究. 地球物理学报, 49(3): 712-717.
[49] 郑伟, 许厚泽, 钟敏等. 2011. 基于时空域混合法利用Kaula正则化精确和快速解算GOCE地球重力场. 地球物理学报, 54(1): 14-21, doi: 10.3969/j.issn.0001-5733.2011.01.003.
[50] 郑伟. 2015. 基于能量守恒原理的卫星重力反演理论与方法. 北京: 科学出版社.
[51] 钟波. 2010. 基于GOCE卫星重力测量技术确定地球重力场的研究[博士论文]. 武汉: 武汉大学.
[52] 钟波, 罗志才, 周浩. 2013. SST-HL模式下利用卫星均值加速度反演地球重力场. 武汉大学学报·信息科学版, 38(8): 902-906.
[53] 周浩, 钟波, 罗志才等. 2011. OpenMP并行算法在卫星重力场模型反演中的应用. 大地测量与地球动力学, 31(5): 123-127.
[54] 周浩, 罗志才, 钟波等. 2015. 利用最小二乘直接法反演卫星重力场模型的MPI并行算法. 测绘学报, 44(8): 833-839.
[55] 周旭华. 2005. 卫星重力及其应用研究[博士论文]. 武汉: 中国科学院测量与地球物理研究所.
[56] 邹贤才. 2007. 卫星轨道理论与地球重力场模型的确定[博士论文]. 武汉: 武汉大学.