论文推荐 | 刘焕玲,文汉江,徐新禹,赵永奇,蔡剑青：GOCE实测数据反演高阶重力场模型的Torus方法

《测绘学报》

构建与学术的桥梁        拉近与权威的距离

复制链接，关注《测绘学报》抖音！

【测绘学报的个人主页】长按复制此条消息，长按复制打开抖音查看TA的更多作品##7NsBSynuc88##[抖音口令]

本文内容来源于《测绘学报》2020年第8期，审图号GS（2020）4062号。

GOCE实测数据反演高阶重力场模型的Torus方法

刘焕玲^1,2,3,文汉江^1,2,徐新禹⁴,赵永奇⁴,蔡剑青⁵                            

1.中国测绘科学研究院, 北京 100036;
2.自然资源部测绘科学与地球空间信息技术重点实验室, 北京 100036;
3.城市空间信息工程北京市重点实验室, 北京 100038;
4.武汉大学测绘学院, 湖北 武汉 430079;
5.德国斯图加特大学大地测量研究所, 斯图加特 70174

基金项目：高分辨率对地观测系统重大专项（42-Y20A09-9001-17/18）；国家自然科学基金（41574019；41774020）；DAAD战略合作伙伴与专题研究合作网络项目（57421148）；中国测绘科学研究院基本科研业务费（AR1922）；城市空间信息工程北京市重点实验室经费资助项目（2019211）；武汉大学地球空间环境与大地测量教育部重点实验室开放基金（18-01-05）

摘要：不同于当前广泛使用的空域法、时域法、直接解法，本文尝试采用Torus方法处理GOCE实测数据，利用71 d的GOCE卫星引力梯度数据反演了200阶次GOCE地球重力场模型，实现了对参考模型的精化。首先，采用Butterworth零相移滤波方法加移去—恢复技术，处理引力梯度观测值中的有色噪声，并利用泰勒级数展开和Kriging方法对GOCE卫星引力梯度数据进行归算和格网化，计算得到了名义轨道上格网点处的引力梯度数据。然后，利用2D-FFT技术和块对角最小二乘方法处理名义轨道上数据，获得了200阶次的GOCE地球重力场模型GOCE_Torus。利用中国和美国的GPS/水准数据进行外部检核结果说明，GOCE_Torus与ESA发布的同期模型的精度相当；GOCE_Torus模型与200阶次的EGM2008模型相比，在美国区域精度相当，但在中国区域精度提高了4.6 cm，这充分体现了GOCE卫星观测数据对地面重力稀疏区的贡献。Torus方法拥有快速高精度反演卫星重力场模型的优势，可以在重力梯度卫星的设计、误差分析及在轨快速评估等方面得到充分应用。

关键词：GOCE    地球重力场    Torus方法    精化

引文格式：刘焕玲, 文汉江, 徐新禹, 等. GOCE实测数据反演高阶重力场模型的Torus方法. 测绘学报，2020，49(8)：965-973. DOI: 10.11947/j.AGCS.2020.20200044.阅读全文：http://xb.sinomaps.com/article/2020/1001-1595/2020-8-965.htm
全文概述

2009年3月17日，欧空局(European Space Agency, ESA)成功发射了GOCE(gravity field and steady-state ocean circulation explorer)卫星。该卫星同时采用高-低卫星跟踪卫星(satellite-to-satellite tracking in high-low mode，SST-hl)和卫星重力梯度测量技术(satellite gravity gradiometry，SGG)，可实现对地球重力场的直接探测。

ESA官方发布的GOCE卫星重力场模型主要采用了直接法、时域法、空域法^[1-3]。此外，还有1D-FFT的半解析法^[4]、张量不变量法^[5]和短弧法^[6]等。直接法采用直接求逆方法解算法方程，在获取精确模型的同时，可提供模型的误差方差-协方差矩阵；时域法将引力梯度的球函数由轨道根数表示，即看作时间t的函数，解算过程中不采用任何先验重力信息，可利用纯GOCE卫星梯度数据确定重力场模型；空域法将沿轨观测数据归算到球面格网上，利用球谐分析方法求解，采用蒙特卡洛(Monte Carlo)统计方法估计模型误差方差-协方差阵。直接法最小二乘解算过程中，在测量带宽内对梯度观测值采用自回归滑动平均(autoregressive moving average, ARMA)滤波；时域法在整个谱域内对观测值进行ARMA滤波；空域法则采用沿轨Wiener滤波和迭代削弱有色噪声的方法。此外，滤波方法还有AR模型滤波^[7]、FIR滤波、移去恢复法和零相位的数字滤波方法^[8]等。本文采用的Torus方法是一种半解析方法，结合二维快速傅里叶变换(2-dimensional fast Fourier transform，2D-FFT)和块对角最小二乘，与直接法和时域法相比，在保证模型精度的情况下，计算效率得到有效提高。而与基于1D-FFT的半解析法相比，本文方法不需要沿轨卫星观测数据的连续性和重复周期的要求^[9-11]。

Torus方法处理卫星引力梯度数据在理论上是严密的^[9-10]，且数值模拟表明该方法是一种有效的GOCE数据处理方法^[11-12]。但该方法要求观测值位于名义轨道，而GOCE卫星实际轨道并非名义轨道，且卫星重力梯度数据的噪声特性较为复杂，使得没有将该方法用于处理GOCE实测数据。现有研究结果已经表明，GOCE引力梯度观测值有助于提高地球重力场中短波信号(特别是80~160阶内系数)的反演精度^[3]，所以本文将基于Torus方法，采用GOCE引力梯度数据对地球重力场模型进行精化，在介绍Torus方法的理论和算法基础上，研究利用Torus方法处理GOCE实测引力梯度观测值的关键步骤处理策略，最终采用GOCE卫星前71 d(2009年11月1日—2010年1月10日)的实测卫星引力梯度数据对地球重力场模型进行了反演，获得了200阶次的高精度GOCE地球重力场模型。

1  Torus方法的基本原理

对于近圆轨道，局部轨道坐标系(local orbit reference frame，LORF)下卫星引力梯度分量与卫星轨道根数之间的关系可表示为^[9-10]

(1)

式中，V_ij表示6个卫星引力梯度分量，ij分别对应(xx、yy、zz、xy、xz、yz)；r为地心向径；I为轨道倾角；u为升交角距，u=ω+f；ω和f分别为近地点角距和真近点角；Λ为升交点经度，Λ=Ω-θ_G；Ω为升交点赤经；θ_G为格林尼治恒星时角。m为球谐函数的次，N为球谐展开的最高阶数，A_mk、B_mk为集总系数，是转换系数H_nmk^V_ij和球谐系数(α_nm和β_nm)的线性组合

(2)

式中，n_min=max(|k|, m)+δ，若k与max(|k|, m)奇偶性相同，则δ=0，否则δ=1。6个卫星引力梯度分量对应的转换系数见式(3)^[10]

(3)

式中，μ为地心引力常数，即万有引力常数与地球质量的乘积；R为地球平均半径；F_nmk^′(I)为倾角函数F_nmk(I)关于轨道倾角I的一阶导数。

式(2)中，球谐系数α_nm和β_nm随n-m的奇偶性取不同的值，具体关系如下

(4)

(5)

式中，C_nm、S_nm为完全规格化球谐系数。

式(1)是SA方法的基本表达式，可看作傅里叶级数，集总系数A_mk和B_mk为傅里叶系数。对名义轨道上的观测值(轨道为圆形、轨道倾角为常量、轨道受J₂项影响而长周期进动)，式(1)可简化为式(6)

(6)

此时卫星引力梯度分量是关于u和Λ的函数，u和Λ的取值范围均为[0, 2π)，这两个不同方向上的圆周可形成一个封闭的圆环，因此该方法称为Torus方法。若观测值均匀分布在Torus圆环面上，则可利用2D-FFT解算集总系数^[13-15]，所以Torus方法可看作2D-FFT对应的SA方法。

2  数据处理方法

2.1  主要流程

利用Torus方法处理GOCE卫星引力梯度观测值，需要获得名义轨道面上LORF下卫星引力梯度格网数据。本文所用数据包括：EGG_NOM_2产品中的卫星引力梯度观测值(GRF坐标系)和IRF向GRF转换的姿态四元数(子产品EGG_IAQ_2)，采样间隔均为1 s；精密科学轨道数据产品SST_PSO_2中的简化动力学轨道(子产品SST_PRD_2)和EFRF向IRF转换的四元数(子产品SST_PRM_2)，采样间隔为10 s。基于Torus方法精化GOCE地球重力场模型的主要流程如图 1所示。

图 1                         Torus方法精化GOCE地球重力场模型的主要流程 Fig. 1     Flow chart of Torus approach to refine the GOCE earth gravity field model      

图选项

2.2  处理策略

2.2.1  卫星引力梯度观测值的滤波

以利用EGM2008模型(前250阶次的系数)模拟的卫星引力梯度值为准确信号(图 2)分析了2009年11月1日—2009年11月10日GRF下GOCE卫星引力梯度数据的精度，如图 3所示。低频部分(约5 mHz以内)，观测值误差为有色噪声，随着频率的降低，误差迅速增大，呈现1/f的变化特点；5~200 mHz测量带宽内，误差近似白噪声，V_xy和V_yz精度较差，噪声水平在几百mE的量级，远大于信号，另外4个分量的精度较高，V_xx和V_yy的噪声水平约10 mE，V_zz、V_xz的噪声水平约20 mE，大于GOCE卫星的设计指标5 mE，相比图 2中相应分量信号的大小，V_xz的信噪比较低；高频部分(> 200 mHz)，误差表现为有色噪声。

图 2 模拟卫星引力梯度值 Fig. 2     The simulated satellite gravity gradiometry data      

图选项

图 3                         GOCE卫星引力梯度观测值中噪声的功率谱密度 Fig. 3     The PSD of the noise in GOCE gravity gradiometry observations      

图选项

Torus方法要求使用LORF下的卫星引力梯度值，直接将GRF下的卫星引力梯度观测值转换至LORF下，精度较差的V_xy和V_yz分量会影响其他分量的精度，且低频有色噪声可能会泄露到测量带宽内的信号中^[16-17]，因此在坐标转换之前，需要对GRF下的引力梯度观测值进行处理。本文利用EGM2008模型前250阶次的系数计算GRF下的引力梯度值，替换观测精度较低的分量V_xy和V_yz；对于其他4个分量中有色噪声的滤波方法，根据GOCE卫星引力梯度观测值中噪声的特性，本文采用带通滤波器，并设计了无限长单位脉冲响应(infinite impulse response，IIR)带通滤波器和有限长单位脉冲响应(finite impulse response，FIR)带通滤波器。对IIR带通滤波器，分别比较了Butterworth、Chebyshev TypeⅠ、Chebyshev Type Ⅱ和Elliptic等不同滤波器原型的滤波效果；对FIR带通滤波器，分别比较了Kaiser、Hamming、Gaussian、Chebyshev、Blackman、Blackman-Harris等不同滤波器原型的滤波效果。为提高滤波的精度，采用移去—恢复策略，为消除相移的影响，对数据进行了零相移滤波。结果表明，采用Butterworth零相移滤波和移去—恢复技术对低频、高频部分的噪声都具有良好的抑制效果^[18]，GRF下V_xx、V_yy、V_zz和V_xz分量观测值残差的滤波结果见图 4。

图 4                         GRF下V_xx、V_yy、V_zz和V_xz分量残差的Butterworth零相移滤波效果 Fig. 4     The filter result of observation residual in GRF using the Butterworth filter with zero-phase combining with remove-restore      

图选项

GOCE卫星引力梯度观测值中，V_xy、V_yz精度较差(滤波时利用模型值代替)，V_xz的信噪比也比较低。而且在GOCE卫星重力场反演中，V_xx、V_yy、V_zz 3个分量的贡献最大，约占总信号的86.86%^[19]，所以本文计算卫星重力场模型时，仅使用了这3个分量。将GRF下的观测值转换至LORF^[19-20]，以EGM2008为参考模型，利用其前200阶次的系数计算模型观测值，计算观测值残差，并绘制功率谱密度图(图 5)，观测值残差中除包含模型误差外，主要反映的是观测值的噪声水平。测量带宽内V_xx、V_yy、V_zz残差的功率谱密度均方根分别在和，所以在模型计算时，3个分量的权重比设为4:4:1。

图 5 滤波后LORF下观测值残差的功率谱密度 Fig. 5     The PSD of the filtered observation residuals in LORF      

图选项

2.2.2  名义轨道的确定

将地固系下卫星位置和速度转换至惯性系下^[22]，计算轨道根数，以平均轨道半径和倾角为名义轨道的半径和倾角，利用泰勒级数展开公式将实际轨道上的GOCE卫星引力梯度观测值归算至名义轨道^[12]，图 6(a)、图 6(b)分别显示了一天内卫星在地球和Torus圆环面上的轨迹，受轨道倾角影响，图 6(a)中存在极空白，纬度θ∈[－90°+θ₀, 90°－θ₀]，其中，θ₀为极空白范围，θ₀=90°－I；而在Torus圆环面上，卫星轨迹覆盖地球两次，Λ和u的取值范围均为0, 360°。对照图 6(a)和图 6(b)，可以看出GOCE卫星轨迹在地球表面及Torus圆环面上的对应关系，方框内均体现了数据异常。图 6中，红色代表上升弧，蓝色代表下降弧。

图 6 卫星轨迹 Fig. 6     Satellite orbits      

图选项

2.2.3  格网化

因在地球重力场反演中引力梯度分量V_zz的影响较大，故以该分量为例，利用无噪声的模拟卫星引力梯度数据分析了顾及移去—恢复技术的Kriging方法的格网化效果^[12]，结果表明该方法的格网化误差在0.110 mE以内，远小于GOCE实测引力梯度观测值中的噪声水平，所以本文利用顾及移去—恢复技术的Kriging方法(EGM2008前200阶次系数作为参考模型)进行格网化，得到名义轨道格网上的引力梯度残差。对引力梯度残差格网进行球谐分析，利用式(1)和式(2)计算球谐系数改正值，结合参考模型得到高精度GOCE地球重力场模型。

3  结果分析

GOCE任务的目标是要恢复100 km空间分辨率的全球重力场模型，因此本文将反演GOCE地球重力场模型的最高阶次设定为200。考虑到本文是为了验证Torus方法对实测数据的处理效果，所以采用与GOCE第一代模型GO_CONS_GCF_2_DIR_R1、GO_CONS_GCF_2_TIM_R1、GO_CONS_GCF_2_SPW_R1(文中分别简称为GO_DIR_R1、GO_TIM_R1、GO_SPW_R1)相同时间段内的观测值，即71 d的GOCE卫星引力梯度数据(2009年11月1日—2010年1月10日)。

地球重力场模型EGM2008由美国国家地理空间局(National Geospatial-intelligence Agency，NGA)研制发布，该模型完全阶次至2159(扩展阶次至2190)，空间分辨率约5′×5′，采用了GRACE卫星跟踪卫星数据、卫星测高数据以及地面重力数据。为了分析GOCE卫星引力梯度观测值对重力场反演的贡献，客观分析Torus方法的反演精度，这里选择没有使用GOCE卫星观测值的EGM2008模型(前200阶次)作为参考模型。

为了解决极空白引起的法方程病态问题，提高带谐项、近带谐项和高阶系数解算的稳定性，对低次项^[20-21]和高阶项(大于170)采用Kaula正则化^[22]，并利用L曲线方法确定正则化参数^[23-24]，最终所得模型简称GOCE_Torus。

以GOCE第5代时域解模型GO_TIM_R5为真实模型，计算模型系数误差，并绘制误差谱和阶误差RMS，分别见图 7(a)和图 7(b)。正则化后，低次项和高阶项的精度明显提高。对于80~160阶系数，由于使用了GOCE卫星引力梯度数据，所以GOCE_Torus模型与GO_TIM_R5更接近，阶误差小于EGM2008。

图 7                         Kaula正则化后Torus模型的误差谱和阶误差RMS Fig. 7     Error spectrum and degree error RMS of Torus model after Kaula's regularization      

图选项

图 8显示了不同模型相对于模型GO_TIM_R5的阶误差RMS。图中所有模型的阶误差RMS均小于模型GO_TIM_R5的信号。由于GOCE_Torus模型受先验模型EGM2008的影响，GO_SPW_R1模型受先验模型EGM2008(仅20阶以内的长波部分)、EIGEN5C和ITG_GRACE2010的影响^{[3, 24]}，GO_DIR_R1模型受先验模型EIGEN5C的影响，3个模型小于80阶系数的精度高于纯GOCE卫星重力场模型GO_TIM_R1。在GOCE卫星引力梯度数据的贡献下，80~160阶内，GOCE_Torus模型精度高于参考模型，但受参考模型的影响，GOCE_Torus模型、GO_DIR_R1模型的精度略低于其他同期模型。而高阶部分，特别是170阶以上的部分，主要体现了正则化的影响。GO_DIR_R1模型利用了球冠正则化方法，由此引入了先验重力场模型Eigen-51c的影响^[3]，GOCE_Torus模型在170阶以上的高阶项进行了Kaula正则化，所以二者在该高频部分的精度较高。而GO_TIM_R1模型的正则化策略与GOCE_Torus模型相似，但精度略低，这是因为GO_TIM_R1模型无先验模型，因此这只反映了GOCE卫星引力梯度数据的反演误差，而且这也说明在170阶以上采用正则化约束的情况下，模型反演的结果受先验模型的影响较大。GO_SPW_R1模型是在整个频域内依据建模信号和噪声的协方差确定是否需要正则化，并确定正则化的权重，精度略低于GO_TIM_R1模型。

图 8 不同模型的阶误差RMS Fig. 8     Degree error RMS of different models      

图选项

利用美国6169个以及中国649个GPS/水准点，对重力场模型进行外部检核，结果分别见表 1和表 2。从结果中可以看出，融入了GOCE观测数据的模型GOCE_Torus、GO_TIM_R1、GO_SPW_R1精度相当，差别仅在1 cm以内，这间接说明Torus方法反演GOCE地球重力场的可靠性。由于GO_DIR_R1模型在解算时采用了GRACE数据，所以其精度略优于其他同期模型的精度。

表 1 地球重力场模型(截断至200阶次)确定的大地水准面高与美国GPS/水准的差值统计Tab. 1 Statistics of difference between geoid height from earth gravity models up to d/o 200 and GPS-leveling data in USAm

Model	Mean	Max	Min	RMS	STD
GO_TIM_R5	-0.567	2.243	-3.056	0.765	0.513
EGM2008	-0.567	2.277	-3.046	0.766	0.516
GO_DIR_R1	-0.570	2.227	-3.018	0.768	0.515
GO_TIM_R1	-0.571	2.274	-3.029	0.775	0.524
GO_SPW_R1	-0.570	2.216	-2.960	0.777	0.528
GOCE_Torus	-0.569	2.221	-3.027	0.769	0.517

表选项

表 2 地球重力场模型(截断至200阶次)确定的似大地水准面高与中国GPS/水准的差值统计Tab. 2 Statistics of difference between geoid height from earth gravity models up to d/o 200 and GPS-leveling data in Chinam

Model	Mean	Max	Min	RMS	STD
GO_TIM_R5	0.047	3.232	-3.007	0.570	0.569
EGM2008	0.047	3.831	-2.882	0.603	0.602
GO_DIR_R1	0.053	3.405	-3.134	0.577	0.575
GO_TIM_R1	0.048	3.345	-3.106	0.578	0.576
GO_SPW_R1	0.044	3.328	-3.062	0.576	0.575
GOCE_Torus	0.048	3.305	-2.750	0.579	0.577

表选项

因中国地形比较复杂，所以相比美国地区，各模型在中国的精度略低，特别是EGM2008模型在反演时未使用GOCE卫星观测数据，且缺乏中国地面重力数据，其误差达0.602 m。相比EGM2008模型，本文模型GOCE_Torus的误差由0.602 m降低至0.577 m，精度提高了2.5 cm，这说明本文基于Torus方法利用GOCE卫星实测引力梯度数据实现了对EGM2008模型的精化；同时，说明在地面重力数据稀少或空白的区域，GOCE卫星数据对提升现有地球重力场模型有明显的贡献。

表 1和表 2中各地球重力场模型均截断至200阶次，考虑到截断误差，利用EGM2008模型200阶以上的系数对各模型进行补充，与GPS/水准数据差异统计结果见表 3。可以看出差异明显减小，例如：EGM2008、GOCE_Torus模型在美国地区的误差分别由0.516、0.517 m减小到了0.284、0.283 m，在中国地区的误差分别由0.602、0.577 m减小到了0.240、0.201 m。各模型在美国地区的精度相当。在中国地区，GO_TIM_R5模型的标准差仅为0.161 m，明显小于其他模型。第一代GOCE模型中，直接解精度仍然最高，为0.179 m，而时域解和空域解的精度与GOCE_Torus模型相当，差别在毫米量级，GOCE_Torus模型的精度比EGM2008模型的高了3.9 cm，这充分说明GOCE数据提高了参考模型EGM2008在中国区域的精度。

表 3 地球重力场模型(200阶以上系数采用EGM2008模型)确定的大地水准面高与GPS/水准的差值统计Tab. 3 Statistics of difference between geoid height from earth gravity models (The omission errors were compensated using the EGM2008 coefficients up to d/o 2190) and GPS-leveling datam

Model	Mean (USA)	STD (USA)	Mean (China)	STD (China)
GO_TIM_R5	-0.511	0.281	0.239	0.161
EGM2008	-0.511	0.284	0.239	0.240
GO_DIR_R1	-0.514	0.284	0.245	0.179
GO_TIM_R1	-0.515	0.295	0.240	0.191
GO_SPW_R1	-0.514	0.298	0.236	0.195
GOCE_Torus	-0.513	0.283	0.241	0.201

表选项

4  结论与展望

本文介绍了Torus方法反演卫星重力场模型的基本原理，给出了利用Torus方法由GOCE卫星引力梯度观测值反演地球重力场模型的数据处理策略，利用GOCE卫星前71 d的引力梯度数据，基于Torus方法对地球重力场模型EGM2008进行了精化，得到了200阶次高精度GOCE地球重力场模型GOCE_Torus。主要结论有：

(1) 将Torus方法用于GOCE实测数据处理中，制定了有效的数据处理流程和策略，联合2D-FFT/IFFT技术和最小二乘方法，可以实现高阶次地球重力场模型的快速高精度反演和精化。

(2) 中国和美国的GPS/水准数据对GOCE_Torus模型精度检核说明，各个卫星重力场模型在美国地区模型精度相当，但在地面重力数据覆盖较少的中国，相比200阶次的EGM2008模型，GOCE_Torus模型精度提高明显，充分体现了GOCE卫星引力梯度观测值的贡献。

Torus方法可以快速高精度地反演卫星重力场，可用于重力梯度卫星的设计、误差分析、在轨快速评估以及卫星重力场模型恢复。本文仅采用了卫星重力梯度观测数据，考虑到卫星跟踪卫星数据对地球重力场的长波部分较为敏感，下一步工作将研究Torus方法在联合利用GRACE的卫星跟踪卫星数据、GOCE卫星引力梯度观测值确定高精度高阶卫星重力场模型中的应用。



作者简介

第一作者简介：刘焕玲（1986—），女，博士，助理研究员，主要研究方向为卫星大地测量、物理大地测量，Email: [e-mail]。

第二作者简介：文汉江（1966—），男，博士，研究员，主要研究方向为卫星测高、卫星大地测量、物理大地测量，Email: [e-mail]。

第三作者（通信作者）简介：徐新禹（1978—），男，博士，副教授，研究方向为地球重力场模型确定及其应用，Email：[e-mail]。

第四作者简介：赵永奇（1990—），男，博士，研究方向为卫星重力测量，Email：[e-mail]。

第五作者简介：蔡剑青（1963—），男，德国工学博士，主要研究方向为经典和现代测量数据处理理论与方法、地球参考坐标系统的建立和转换、卫星大地测量与导航、物理大地测量和地球动力学解释等，Email：[e-mail]。

权威 | 专业 | 学术 | 前沿

微信、抖音小视频投稿邮箱 | [e-mail]

微信公众号中搜索「测绘学报」，关注我们，长按上图二维码，关注学术前沿动态。

欢迎加入《测绘学报》作者QQ群： 751717395

进群请备注：姓名+单位+稿件编号