-
随着星载全球导航卫星系统(global navigation satellite system,GNSS)精密定轨技术的不断发展,挑战性小卫星有效载荷(challenging mini-satellite payload,CHAMP)、重力反演和气候实验(gravity recovery and climate experiment,GRACE)、地球重力场和海洋环流探测卫星(gravity field and steady-state ocean circulation explorer,GOCE)、GRACE Follow-On重力卫星任务相继成功实施,地球重力场的中长波探测精度有了2~3个量级的提高[1-2]。作为直接探测地球重力场的GOCE卫星任务,其径向梯度观测值精度仅能达到10~20
(1 E=1×10-9/s2)的精度水平,未达到重力梯度仪的设计指标。纯GOCE卫星重力场模型GO_TIM_R5在100 km空间分辨率(除极地外)上的全球大地水准面误差可达到2.4 cm[3],这与满足相关科学领域的需求、实现大地水准面精度100 km半波长上优于1 cm的大地测量总体目标尚存在一定差距。为进一步提升全球重力场的探测精度,服务于全球高程系统统一等应用需求,中国重力梯度测量卫星任务已提上发展日程,开展利用卫星重力梯度数据确定地球重力场的理论方法研究仍是需要继续探索的基础科学问题。 在利用卫星重力梯度数据确定地球重力场的经典理论方法中,有直接法、空域法和时域法等[1-4]。为了构建观测模型,梯度分量需要在模型参考框架和梯度观测框架间进行转换,此时参考框架的旋转量的确定至关重要。实际上,在参考框架旋转过程中,单个梯度分量是变化的,引力梯度张量是不变的。为了实现参考框架间的转换,需要获取额外的定向信息,也就是需要确定梯度仪参考框架相对于某一确定的参考框架(例如惯性参考框架)的定向参数[5]。星敏感器可以提供该类信息,但是其测量精度有限,飞行器本身的定向就是一个主要的误差源。
张量不变理论在连续介质力学、计算力学等诸多科学领域得到了广泛应用[6]。随着GOCE重力梯度测量任务的实施,其在精化地球重力场方面的应用研究逐渐增多。文献[5,7]较早地将张量不变理论应用于利用卫星重力梯度技术确定地球重力场,文献[8-12]对张量不变函数模型的线性化问题等进行了系统研究。文献[13-14]推导了球近似条件下的张量不变方法的线性化模型表达式,给出了相应的球谐级数解,并利用GOCE卫星数据解算得到了地球引力场模型;文献[15]针对GOCE卫星任务非全张量观测数据类型,分析了张量不变量的计算误差;文献[16]提出了一种利用半解析法解算GOCE引力梯度张量不变量观测值来恢复全球卫星重力模型的方法。
本文应用张量不变理论对利用卫星重力梯度数据确定地球重力场的方法和关键问题进行了探讨,包括张量不变观测方程的线性化处理、非全张量观测值的数据处理策略以及采用白噪声特性下的梯度观测值恢复地球重力场的数值分析等,相关结论可为卫星重力梯度反演地球重力场的研究提供参考。
-
根据张量不变理论,对称二阶引力梯度张量
具有3个独立的标量不变量,即一、二、三阶主子式之和,它们组成了一个完备的不变系统。该不变系统可表达为[9]: 式中,Vij(i,j=1,2,3)表示引力梯度张量的6个分量;I1、I2、I3为引力梯度张量的3个不变量。
在局部指北坐标系中,引力位可表达为:
$$ \begin{gathered} V(r, \theta, \lambda)=\frac{G M}{R} \sum_{n=0}^{\infty} \sum_{m=0}^{n}\left(\frac{R}{r}\right)^{n+1}\left(\bar{C}_{n m} \cos m \lambda+\right. \\ \left.\bar{S}_{n m} \sin m \lambda\right) \bar{P}_{n m}(\cos \theta) \end{gathered} $$ (2) 式中,
为地心引力常数; 为地球平均半径; 分别为地心向径、地心余纬和地心经度; 为完全规格化缔合勒让德函数,简记为 ; 、 为完全规格化的引力位球谐系数;n、m为球谐展开的阶和次。将引力位截断至一定的阶次N,代入式(1),可得到3个不变量的表达式[9]。 -
该不变量是引力梯度张量系数矩阵3个主对角元素的和,且引力位系数与不变量观测值间的函数模型为线性模型,即:
式中,
为余纬 和阶次 、 的函数: 式(4)即勒让德差分方程。
I1不变量为引力梯度张量的迹。由于
,因此无法利用其进行地球重力场的恢复。但是它可以作为最小二乘方法的先验约束条件,从而提高求解的精度。 -
该不变量为引力梯度张量6个分量的组合形式,是引力梯度分量平方的线性组合形式。但相对于球谐系数,I2为一次非线性模型。
球谐系数的乘积形式,例如
、 等,导致了不变量 的非线性特征。其中,n2、m2分别为球谐展开的阶和次,且: 式中,
表示完全规格化缔合勒让德函数 相对于地心余纬 的导数。 $$ \begin{gathered} K_{n_{1} m_{1} n_{2} m_{2}}^{2}(\theta)=-\frac{1}{2}\left\{\left[-\left(n_{1}+1\right) \bar{P}_{n_{1} m_{1}}-\right.\right. \\ \left.\frac{m_{1}^{2}}{\sin ^{2} \theta} \bar{P}_{n_{1} m_{1}}+\cot \theta \frac{\mathrm{d} \bar{P}_{n_{1} m_{1}}}{\mathrm{~d} \theta}\right] \times \end{gathered} $$ -
不变量I3亦为一非线性函数模型:
相对于不变量I2,I3的非线性特征由球谐系数的乘积形式
等引起,为二次非线性模型。其中, 为球谐展开的阶和次,且: 需要指出的是,梯度仪相对于惯性空间的旋转特性所带来的旋转效应会对张量不变量的计算带来比较大的影响。因此,在利用卫星重力梯度数据恢复地球重力场时,需要首先从观测张量中消去欧拉和离心张量[5,9,17-19]。为了消去这种旋转效应,在张量不变方法中,角速度和角加速度信息同样需要获取,但不需要精确的定向参数信息,这与经典的重力梯度反演方法有别。下文的研究中,假定已经对引力梯度张量和旋转张量进行了分离处理。
-
为了在最小二乘准则下利用张量不变量进行地球重力场的恢复研究,寻找有效的线性化方法,实现函数模型式(5)和式(8)的线性化处理至关重要。用于评价线性化效果的标准一般有3个:(1)线性化误差的量级;(2)收敛特性的优劣;(3)数值运算的效率。对于大型的最小二乘问题,第3个标准特别重要。Brute-force线性化方法,即利用泰勒级数展开实现模型的线性化,并不适用于大型的最小二乘问题,其带来的嵌套求和过程将会增加大量的计算量,即便对于中低阶的重力位模型的求解亦会存在这种问题。
针对张量不变量观测值恢复地球重力场中的实际问题,本文采用文献[9]中的线性化方法对张量不变量进行了线性化处理[9,20],最大限度地减少了海量数据带来的计算负担,即利用一先验的参考重力位模型计算张量不变量的扰动量,同时达到线性化处理的目的。从第二次迭代开始,可以利用上次迭代的求解结果来建立线性函数模型。张量不变量扰动量的表达式为[9]:
$$ \begin{gathered} \delta I_{2}=I_{2}-I_{2}^{\mathrm{ref}}=-U_{11} \delta V_{11}-U_{22} \delta V_{22}- \\ U_{33} \delta V_{33}-2\left(U_{12} \delta V_{12}+U_{13} \delta V_{13}+\right. \\ \left.U_{23} \delta V_{23}\right)+o^{2}\left(\delta V_{i j}\right) \end{gathered} $$ (11) $$ \begin{aligned} \delta I_{3}=& I_{3}-I_{3}^{\mathrm{ref}}=U_{11} U_{22} \delta V_{33}+U_{11} \delta V_{22} U_{33}+\\ & \delta V_{11} U_{22} U_{33}+2\left(U_{12} U_{13} \delta V_{23}+\right.\\ &\left.U_{12} \delta V_{13} U_{23}+\delta V_{12} U_{13} U_{23}\right)=-\delta V_{11} U_{23}^{2}-\\ & \delta V_{22} U_{13}^{2}-\delta V_{33} U_{12}^{2}-2\left(U_{11} U_{23} \delta V_{23}+\right.\\ &\left.U_{22} U_{13} \delta V_{13}+U_{33} U_{12} \delta V_{12}\right)+\\ & o^{2}\left(\delta V_{i j}\right)+o^{3}\left(\delta V_{i j}\right) \end{aligned} $$ (12) 式中,
分别为利用参考模型进行球谐综合得到的张量不变量;Uij(i,j=1,2,3)为参考模型的引力梯度张量分量,且 , 为对应于参考量Uij的改正量。 该线性化处理可以在任意的参考系下进行,即线性化结果与引力梯度张量的方向无关。当参考梯度值的频谱分辨率与参数估计的频谱分辨率相等时,即阶数Nref=N时,线性化效果最优。该线性化方法增加的计算量主要集中在每一次迭代前参考模型值Uij的球谐综合计算。
本文利用EIGEN-6C4的前300阶计算得到了30 d、5 s采样的卫星引力梯度观测值,卫星轨道仿真计算时的初始轨道根数为:轨道长半径6 628.000 km,偏心率0.001,轨道倾角96.6°,升交点赤经0.0°,近地点角距0.0°,平近点角0.0°。图 1给出了全张量梯度观测值下,利用先验重力场模型进行模型线性化后的求解结果。为有效测试解的收敛特性,初始线性化位模型选为EGM96模型的前200阶。张量不变量选为I2,解算阶数为200阶。从图 1中可以看出,在第2次迭代后,精度不再有明显提高,即第2次迭代已满足收敛条件。此外,第1次迭代后的模型阶误差明显小于单独利用V33分量的求解精度,说明了线性化方法的误差量级非常小,可有效实现张量不变模型的线性化处理。
图 2显示了分别利用EGM96模型和GGM05C模型作为初始线性化模型的第一次迭代后的阶误差。从图 2中可以看出,以GGM05C作为线性化模型第一次迭代后的精度较EGM96模型要好,其差异主要集中在长波部分,且GGM05C模型第一次迭代阶误差非常接近于线性化后的最终收敛值。这是由于GGM05C与EIGEN-6C4模型中均使用了GRACE观测数据参与解算,相对于EGM96模型,GGM05C初始线性化模型更为精确。初始线性化模型的优劣与初始迭代精度是密切相关的,但对最终收敛结果的影响非常小,不会导致收敛特性的改变。在实测数据处理时,为了提高收敛速度,保证收敛特性不变,应选取精度更高的重力场模型作为初始线性化模型。
从图 1和图 2中还可以看出,张量不变方法的重力场模型解算精度要优于V33分量的解算精度。一方面,本文梯度观测值仿真计算时采用了300阶模型,而解算只解算到了200阶,这主要是考虑到卫星在轨的实际情况,这种频谱混淆在实际重力场模型解算中也是存在的。另一方面,由于张量不变法采用了先验模型进行线性化,观测数据的数学统计特性同时得到了改善和优化,加之线性化误差较为微小,因此其精度结果优于V33分量。
-
函数模型(5)和(8)的计算需要全张量梯度仪的观测值。不同梯度张量分量的误差特性会统一反映在张量不变量中。因此,所有的梯度张量分量应具有同等量级的测量精度。对于GOCE重力任务,张量分量V12和V23的精度较之其他分量要低,需要寻找一定的方法对这种情况进行处理。本文中,在进行不变量的计算时,首先假定V12和V23观测值不可知,利用先验模型,结合球谐综合方法得到或直接将其赋零,然后进行重力场模型的解算。
图 3给出了将V12和V23置为零后,利用张量不变方法解算得到的重力场模型阶误差。从图 3中可以看出,第一次迭代的精度较之全张量梯度观测值的求解结果要低一个数量级,但第二次迭代即满足收敛条件,精度仅仅稍逊于全张量的求解结果。图 4给出了利用EGM96位模型的球谐综合值代替V12和V23分量后的张量不变方法的解算结果。从图 4中可以看出,第一次迭代结果较之赋零策略的结果精度要高,且第二次迭代收敛值与全张量求解结果吻合很好。这说明,对于非全张量梯度观测值,在利用张量不变方法求解时,借助先验位模型的球谐综合值代替低精度分量的数据处理策略是可行和有效的。
图 4 V12和V23分量球谐综合计算后的张量不变解
Figure 4. Tensor Invariant Solution After Spherical Harmonic Synthesis of V12 and V23
GOCE卫星搭载的梯度仪具有天底指向特性,使得引力梯度张量非对角元素的变化范围在几个Eötvös左右。与之相比,主对角线分量的量级要大3个量级左右。图 5给出了将低精度观测值V12和V23分别置为零以及利用EGM96位模型球谐综合值进行代替时的张量不变量的相对误差。从图 5中可以看出,将V12和V23置为零或者利用先验位模型的球谐综合值进行代替对张量不变量的量级影响很微小,这种良好的数学特性也是选择不变系统进行重力场研究的初衷之一。因此,在利用张量不变方法求解卫星引力梯度重力场时,建议借助先验位模型的球谐综合值代替低精度分量进行处理。
-
对梯度张量模拟值的V11、V22、V33和V13分量分别加入标准差为3.2 mE的白噪声,利用张量不变方法进行重力位模型的解算。其中,V12和V23分量采用EGM96模型的球谐综合值代替,正则化方法选用Kaula正则化,结果见图 6。从图 6中可以看出,白噪声情况下,张量不变方法的模型阶误差要优于单分量的解算结果,这与无噪声情况下的结果一致。但迭代次数的增加并不能改进模型阶误差的精度。这说明模型的线性化误差以及采用先验模型的球谐综合值代替低精度分量所引入的误差较之白噪声水平要低,利用先验重力场模型进行线性化,以及球谐综合值代替低精度分量的数据处理策略仍然是有效的。
-
本文对张量不变理论在卫星重力梯度数据恢复地球重力场中的应用进行了数值分析。针对观测模型的非线性特征以及GOCE任务的非全张量观测特征,对基于先验重力场模型的线性化方法以及球谐综合方法的数据处理策略进行了探讨,结果表明:
1)张量不变方法不依赖于梯度仪的空间定向信息的获取,较之传统的地球重力场解算方法具有特别的优越性。该方法实现了不同观测值的联合求解,降低了计算负担,特别是对于海量的卫星重力梯度观测值尤为可贵,且张量不变解的精度明显高于单梯度分量的求解结果。
2)基于先验重力场模型的线性化方法能够在保证计算效率的前提下,有效实现观测模型的线性化,解的收敛非常快,精度较高,线性化误差较小,且初始线性化模型的选择对解的收敛特性影响较小。
3)利用先验位模型的球谐综合值代替低精度梯度分量实现非全张量梯度观测值的处理是有效的,不会改变解的收敛速度和最终精度。这主要源于梯度仪的天底指向特性。
4)白噪声情况下,基于先验重力场模型的线性化方法以及采用先验模型的球谐综合值代替低精度分量的数据处理策略仍然是有效的,其引入的误差较之白噪声水平要低,从而迭代次数的增加并不能带来模型阶误差精度的显著提高。
Tensor Invariant Method for Determining the Earth Gravity Field from Satellite Gravitational Gradient Data
-
摘要: 应用张量不变理论对利用卫星重力梯度数据确定地球重力场的方法进行了研究,对张量不变观测方程的线性化处理、非全张量观测值的数据处理策略以及采用白噪声特性下的梯度观测值恢复地球重力场的精度等进行了数值分析。结果表明,张量不变解实现了不同观测值的联合求解,基于先验重力场模型的线性化方法在实现张量不变观测模型线性化处理的同时,提升了张量不变方法的模型解算精度,且初始参考模型的选择对解的收敛特性影响较小;在无噪声和白噪声条件下,利用先验位模型的球谐综合值代替低精度梯度分量实现非全张量梯度观测值的处理均是可行的,不会改变解的收敛速度和最终精度。Abstract:
Objectives The objective is to study the tensor invariant method for determining the Earth's gravity field with satellite gravitational gradient data. Methods The tensor invariant method does not depend on the acquisition of the spatial orientation information of the gradiometer, so it has a special advantage over the traditional methods in determining the earth gravity field. Numerical analyses are made for the linearization of tensor invariant observation equation, the processing of non-full tensor observation data, and the accuracy of gravity field recovery with the gradient data containing white noise. Results The results show that the tensor invariant solution realizes the combine processing of the different tensor components. The tailored linearization effectively achieves the linearization of the tensor invariant observation model while the accuracy of gravity field model determined using tensor invariants is significantly higher than that of the single tensor component solution. Also, the choice of initial reference model exerts minor influence on the total iterative process. Conclusions Due to the non-full tensor gradient observations, it is proposed to replace the low precision gradient component using the spherical harmonic synthesis of a priori gravity field model, which will not change the convergence and the accuracy of the solution even if the white noise exists in the observations because of the nadir-pointing characteristics of the gradiometer. -
Key words:
- tensor invariant /
- satellite gravitational gradient /
- linearization /
- gravity field model
-
-
[1] Bruinsma S L, Förste C, Abrikosov O, et al. ESA's Satellite-Only Gravity Field Model via the Direct Approach Based on all GOCE Data[J]. Geophysical Research Letters, 2014, 41(21): 7508-7514 doi: 10.1002/2014GL062045 [2] Pail R, Bruinsma S, Migliaccio F, et al. First GOCE Gravity Field Models Derived by Three Different Approaches[J]. Journal of Geodesy, 2011, 85(11): 819-843 doi: 10.1007/s00190-011-0467-x [3] Sneeuw N. Space-Wise, Time-Wise, Torus and Rosborough Representations in Gravity Field Modelling [J]. Space Science Reviews, 2003, 108(1/2): 37-46 doi: 10.1023/A:1026165612224 [4] 苏勇, 范东明, 蒲星钢, 等. 联合GOCE卫星数据和GRACE法方程确定SWJTU-GOGR01S全球重力场模型[J]. 武汉大学学报·信息科学版, 2018, 43(3): 457-463 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201803019.htm Su Yong, Fan Dongming, Pu Xinggang, et al. New Static Gravity Field Model SWJTU-GOGR01S Derived from GOCE Data and GRACE Normal Equation[J]. Geomatics and Information Science of Wuhan University, 2018, 43(3): 457-463 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201803019.htm [5] Rummel R. Satellite Gradiometry[J]. Lecture Notes in Earth Sciences, 1986, 7: 317-363 [6] 宋来忠. 张量分析及其应用[M]. 武汉: 武汉大学出版社, 2016: 233-238 Song Laizhong. Tensor Analysis and Application [M]. Wuhan: Wuhan University Press, 2016: 233-238 [7] Holota P. Boundary Value Problems and Invariants of the Gravitational Tensor in Satellite Gradiometry [M]//Sanso F, Rummel R. Tensor in Satellite Gradiometry. Berlin: Springer, 1989 [8] Baur O, Grafarend E W. High-Performance GOCE Gravity Field Recovery from Gravity Gradient Tensor Invariants and Kinematic Orbit Information[J]. Obversion of the Earth System from Space, 2006, 32(11): 239-253 doi: 10.1007/3-540-29522-4_17 [9] Baur O, Sneeuw N, Grafarend E W. Methodology and Use of Tensor Invariants for Satellite Gravity Gradiometry[J]. Journal of Geodesy, 2008, 82(4/5): 279-293 [10] Baur O, Seneeuw N, Cai J, et al. GOCE Data Analysis: Realization of the Invariants Approach in a High Performance Computing Environment[J]. Signal, 2010, 12(6): 79-86 http://www.gis.uni-stuttgart.de/res/papers/BAUR_2010f.pdf [11] Oruç B. Depth Estimation of Simple Causative Sources from Gravity Gradient Tensor Invariants and Vertical Component[J]. Pure and Applied Geophysics, 2010, 167(10): 1259-1272 doi: 10.1007/s00024-009-0021-4 [12] Cai J, Sneeuw N. Stochastic Modeling of GOCE Gravitational Tensor Invariant[C]//Observation of the System Earth from Space - CHAMP, GRACE, GOCE and Future Missions, Berlin, Germany, 2014 [13] 于锦海, 万晓云. 利用引力梯度不变量解算的GOCE引力场模型[J]. 中国科学: 地球科学, 2012, 42(9): 1450-1458 https://www.cnki.com.cn/Article/CJFDTOTAL-JDXK201209014.htm Yu Jinhai, Wan Xiaoyun. Recovery of the Gravity Field from GOCE Data by Using the Invariants of Gradient Tensor[J]. Scientia Sinica (Terrae), 2012, 42(9): 1450-1458 https://www.cnki.com.cn/Article/CJFDTOTAL-JDXK201209014.htm [14] 于锦海, 赵东明. 引力梯度不变量与相关边界条件[J]. 中国科学: 地球科学, 2010, 40(2): 178-187 doi: 10.3969/j.issn.0253-2778.2010.02.013 Yu Jinhai, Zhao Dongming. the Gravitational Gradient Tensor's Invariants and the Related Boundary Conditions[J]. Scientia Sinica(Terrae), 2010, 40(2): 178-187 doi: 10.3969/j.issn.0253-2778.2010.02.013 [15] 吴星, 王凯, 冯炜, 等. 基于非全张量卫星重力梯度数据的张量不变量法[J]. 地球物理学报, 2011, 54(4): 966-976 doi: 10.3969/j.issn.0001-5733.2011.04.011 Wu Xing, Wang Kai, Feng Wei, et al. Method of Tensor Invariant Based on Non-full Tensor Satellite Gravity Gradients[J]. Chinese Journal of Geophysics, 2011, 54(4): 966-976 doi: 10.3969/j.issn.0001-5733.2011.04.011 [16] 李建成, 徐新禹, 赵永奇, 等. 由GOCE引力梯度张量不变量确定卫星重力模型的半解析法[J]. 武汉大学学报·信息科学版, 2016, 41(1): 21-26 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201601004.htm Li Jiancheng, Xu Xinyu, Zhao Yongqi, et al. Approach for Determining Satellite Gravity Model from GOCE Gravitational Gradient Tensor Invariant Observations[J]. Geomatics and Information Science of Wuhan University, 2016, 41(1): 21-26 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201601004.htm [17] 罗志才. 利用卫星重力梯度数据确定地球重力场的理论和方法[D]. 武汉: 武汉大学, 1996 Luo Zhicai. Theory and Method for Determining the Earth's Gravity Field from Satellite Gravity Gradiometry[D]. Wuhan: Wuhan University, 1996 [18] 徐新禹. 卫星重力梯度及卫星跟踪卫星数据确定地球重力场的研究[D]. 武汉: 武汉大学, 2008 Xu Xinyu. Study of Determining the Earth' s Gravity Field from Satellite Gravity Gradient and Satelliteto-Satellite Tracking Data[D]. Wuhan: Wuhan University, 2008 [19] 朱广彬. 卫星重力梯度测量技术确定地球重力场的理论方法研究[D]. 武汉: 武汉大学, 2010 Zhu Guangbin. Research on the Theory and Methodology for the Earth' s Gravity Field Determination Using Satellite Gravity Gradiometry Measurement[D]. Wuhan: Wuhan University, 2010 [20] Kern M. An Analysis of the Combination and Downward Continuation of Satellite, Airborne and Terrestrial Gravity Data[D]. Calgary: University of Calgary, 2003 -