地球物理学报  2011, Vol. 54 Issue (4): 966-976   PDF    
基于非全张量卫星重力梯度数据的张量不变量法
吴星1, 王凯2, 冯炜3, 汪涛1     
1. 总装备部工程设计研究总院,北京 100028;
2. 信息工程大学测绘学院,郑州 450052;
3. 北京环球信息应用开发中心,北京 100094
摘要: 在非全张量卫星重力梯度观测数据的处理过程中,由于卫星姿态角误差、梯度观测数据误差和非全张量观测等原因,重力梯度值从卫星重力梯度仪系转换到地固系后,精度损失严重.本文研究了张量不变量法以解决上述问题.首先在重力梯度张量不变量线性化的基础上,建立了基于卫星轨道面的不变量观测模型,完整地推导了两类重力梯度张量不变量的球近似和顾及地球扁率影响的球面边值问题的求解公式.针对GOCE卫星任务非全张量观测数据类型,分析了张量不变量的计算误差;结果表明,重力梯度观测误差在不变量的计算中并没有被放大.最后运用广义轮胎调和分析方法进行了模拟试验,数值试验证明,在卫星姿态误差较大时,处理张量不变量比处理张量分量更具优势,并且张量不变量法能有效地解决非全张量观测的问题.
关键词: 卫星重力梯度      张量不变量      地球重力场      调和分析      GOCE     
Method of tensor invariant based on non-full tensor satellite gravity gradients
WU Xing1, WANG Kai2, FENG Wei3, WANG Tao1     
1. Center for Engineering Design and Research under the Headquarters of General Equipment, Beijing 100028, China;
2. Institute of Surveying and Mapping,Information Engineering University,Zhengzhou 450052,China;
3. Global Information Application and Development Center of Beijing, Beijing 100094,China
Abstract: Due to the errors of satellite attitude, gradient observation as well as non-full tensor observation, the precision of gravity gradient values decreased greatly after transformation from satellite gradiometer coordinate system to earth-fixed coordinate system in the processing of non-full tensor observational data of satellite gravity gradient. The method of tensor invariant is mainly studied to overcome the issue mentioned above. Firstly, with the linearization of gravity gradient tensor, invariant observational model is established on the satellite orbital plane. Formulae solutions of sphere boundary problem of two invariants both in sphere approximation and considering earth eccentricity are given completely. Then, computing errors of tensor invariant are analyzed aiming at non-full tensor observational data of GOCE mission. The result shows that observational errors are not magnified in the processing. Finally, general torus harmonic analysis is introduced in the simulation, which shows that the processing result of tensor invariant is better than that of tensor components with big satellite attitude errors, and the method of tensor invariant could deal with the non-full tensor observational data effectively.
Key words: Satellite gravity gradient      Tensor invariant      Earth's gravity field      Harmonic analysis      GOCE     
1 引言

国际上首颗重力梯度卫星GOCE (Gravityfieldandsteady-stateOceanCirculationExplorer)已经于2009年3月成功发射升空.它将以前所未有的精度和分辨率测定地球重力场,有望恢复200 阶地球重力场模型和1~2cm 精度的大地水准面.这将对大地测量学、地球物理学、海洋学等相关地学学科产生深远意义.

GOCE 计划数据分析工作主要包括预处理[1]、卫星的位置及速度计算[2]、地球重力场的计算以及相关数据的其他科学应用[3].其中,运用重力梯度数据反演地球重力场得到了国内外学者的广泛研究,目前主要有三大类方法:直接法[4]、时域法[5]和空域法[6].

由于重力场位系数是在地固坐标系中解算的,因此必须将梯度仪参考系(GRF:GradimeterReferenceFrame)中观测得到的重力梯度张量转换到地固坐标系中,因此不但需要全张量梯度观测,而且需要确定卫星的姿态.然而姿态角误差在均方根超过3~5个角秒时,将成为地球重力场数据处理中的主要误差[7].GOCE 卫星设计的姿态误差控制在几个角秒内,因此对于地球重力场解而言,姿态信息误差的影响可能大于重力梯度数据误差的影响.另外,GOCE 任务的梯度观测并不是全张量观测,仅能有效给出4个分量,而其余2 个非对角分量观测误差扩大近千倍,并且其误差将在新坐标系下放大,并几乎与其本身信号具有相同的数量级[8].因此非全张量观测也会增加数据处理的难度.由此可见,由于姿态角误差、梯度观测数据误差和非全张量观测等原因,在坐标系之间的转换将大大增加重力梯度张量的误差.然而目前的大多数方法忽略了坐标系转换带来的误差.

因此,寻求与坐标系选择无关的数据处理方法就显得尤为重要,而其中最为有效的方法就是运用重力梯度张量不变量法.早在1986 年,Rummel就研究过重力梯度张量不变量[9],随后,Holota、Sacerdote和Sansò、Vermeer等[10~12]学者都曾讨论过重力梯度张量不变量,但由于当时计算机技术的限制,他们的研究仅停留在理论与概念阶段,没有得到实用可行的方法,也没有进行有效的数值试验;直到最近几年,才有了实质性发展,Baur首次运用摄动理论得到了重力梯度张量不变量的线性化模型,并采用了最小二乘迭代法进行了数值分析与模拟试算[1314],于锦海和赵东明[15]在球坐标系下的导出了球近似模型与顾及J 2 项影响的模型,进一步简化了不变量的线性模型.

但是,现有文献没有针对卫星重力梯度张量非全张量观测、不同卫星姿态误差等情况进行深入研究,本文针对上述问题,详细完整的给出了两类不变量边值问题的求解公式,运用EGM2008模型[16]模拟得到了全球分布的卫星重力梯度张量数据,并采用广义轮胎调和分析方法[17],比较分析了不同卫星姿态误差对重力梯度张量分量和不变量解的影响,进一步证实了当存在卫星姿态误差的情况下,重力梯度张量不变量方法更具有优势.

2 卫星姿态误差对重力梯度观测的影响分析

O-X1X2X3 表示直角形式的地固坐标系,O′-x1x2x3 表示梯度仪参考系(GRF),V表示地球引力位,VijVij(ij= 1,2,3)分别表示在O-X1X2X3O′-x1x2x3 坐标系下的重力梯度;VV为相应的重力梯度张量矩阵.若用S′表示GOCE卫星的轨道面,则GOCE 卫星计划的观测结果就给出了S′ 的形状与S′ 上Vij值.设两坐标系的欧拉旋转角表示为ε= {αβγ},则有

(1)

其中R为旋转矩阵:

(2)

一般ε比较小,因此有sinεε,cosε≈1,另外忽略非线性项εi·εj(ij=1,2,3),可得如下线性关系式:

(3)

根据误差传播定律,由(3)式可得重力梯度张量分量Vij的误差传播公式;以分量V33 为例,假设各个分量互不相关,则可得V33 的误差传播公式:

(4)

在平均值的情况下,V 可近似地表示为[8]

(5)

式中,E=10-9m/s2.由(4)式和(5)式可以估算,当姿态角误差分别为10、20、50角秒时,对V33可造成1、2、5 mE 的误差(mE=10-12m/s2),显然,这对精度达mE 级重力梯度观测而言,是不能忽视的.由此可见,Vij转换为Vij的计算精度不仅受到观测量Vij的误差影响,还受到卫星姿态控制误差的影响.然而,目前关于重力梯度的讨论多数是从Vij出发展开的,忽略了卫星姿态控制误差的影响.

因此本文将研究直接由Vij计算得到重力梯度的张量不变量作为计算量反演地球重力场的方法.

3 重力梯度张量的不变量

梯度张量的特征方程[13]

(6)

的三个解,即为梯度张量的三个特征值λi(i=1,2,3),是不变量;特征方程系数Ii(i=1,2,3)也是不变量并且有

(7)

(8)

(9)

(7) ~(9)式是坐标系O-X1X2X3 下的表示形式,若用坐标系O-x1x2x3 来表示,则仅需将VijVij代替即可.它们的值与坐标系O-X1X2X3O-x1x2x3的选择无关.显然,I1 就是梯度的迹,其值必为零.因此,本文重点从I2I3 的值出发来讨论求解地球重力场的问题.

4 不变量观测方程的建立 4.1 线性化

由(8)式和(9)式可知,不变量I2I3 关于引力位V是非线性的,因此需要进行线性化.在此引入所谓参考引力位U,一般可以选择现有的地球重力场模型,也可取如下旋转椭球的重力场[18]:

(10)

式中,(rθλ)表示球坐标系,μ =GM为地球引力常数,R为地球平均半径,P2n(cosθ)为2n阶勒让德多项式;则扰动位为T=V-U,相应的扰动梯度张量为Tij=Vij-Uij.设参考重力梯度张量对应的不变量分别为I2 refI3 ref ,并将Vij=Uij+Tij分别代入(8)式和(9)式,略去O(T2)以上的高阶小量,顾及Laplace方程ΔT=T11 +T22 +T33 =0,可得

(11)

(12)

(11) 式和(12)式就是关于扰动梯度张量Tij在轨道面S′ 上线性化公式,精度与O(T2)同阶.

为了便于计算位系数,坐标系O-X1X2X3 一般采用对应的球坐标系(rθλ),顾及ΔT=0,ΔU=0,由梯度张量坐标转换式(详见文献[19],(2.40)式),可得

(13)

(14)

将(13)式和(14)式代入(11)式和(12)式将得到更为复杂的形式.因此下面将进一步研究球近似下的不变量观测方程.

4.2 球近似下的不变量观测方程

球近似就是在线性化的基础上,将U近似的取为$\tilde{U}$ =μ/r,则(14)式可简化为

(15)

将(13)式和(15)式代入(11)和(12)式,并顾及ΔT=0,可得球近似下的不变量:

(16)

(17)

根据(16)式和(17)式,卫星轨道面S′ 上点P的观测方程可表示为

(18)

(19)

式中,P表示重力梯度观测量的位置,可由星载GPS数据确定,而δI2(P)与δI3(P)则直接来自于重力梯度仪在P点的观测数据,v2(P)和v3(P)为误差改正项.

又因为扰动重力位可表示成[19]

(20)

式中,(CnmSnm)为扰动重力位系数,Pnm(cosθ)为nm次完全正常化Legendre函数.为了表达简洁,分别将Pnm(cosθ),$\frac{{d{{\bar P}_{mn}}\left( {\cos \theta } \right)}}{{d\theta }},\frac{{{d^2}{{\bar P}_{mn}}\left( {\cos \theta } \right)}}{{d{\theta ^2}}}$简写为Pnm(θ),Pnm(1)(θ),Pnm(2)(θ).

由(20)式对r二次求导可得

(21)

把(20)式分别代入(18)式和(19)式,可得

(22)

(23)

(22) 式和(23)式可进一步简写为

(24)

(25)

其中,L2L3 为重力梯度张量不变量向量,$\tilde{A}$2 和$\tilde{A}$3为设计矩阵,Y为重力场模型参数向量,共有(N+1)2 个,其中(N+1)(N+2)/2 个CnmN(N+1)/2个Snm(m>0).

球近似使原来线性化的精度由O(T2)量级降为O(J2·T)量级.为了进一步提高精度,在(11)式和(12)式中对于参考引力位的选择必须考虑J 2 项的影响.下面给出顾及J 2 项的不变量观测方程.

4.3 顾及J 2 项的不变量观测方程

为了与球近似下的引力位U 珦有所区别,设顾及J 2 项的引力位为$\hat U$,即

(26)

式中,J 2≈0.0 0 1 0 8 2 6 3.顾及关系式(1 4) ,可得

(27)

将(1 3)式和(2 7)式分别代入(1 1)式和(1 2)式,可得顾及J 2 项影响的重力梯度不变量为

(28)

(29)

其中,δ T r r,2δ T r r,3 分别表示为

(30)

(31)

由(20)式,分别可得

(32)

(33)

把(21)式、(32)式和(33)式分别代入(28)式和(29)式,可得

(34)

(35)

由于地球引力位V中除去$\hat U$ 外的各项系数大致都是O(J22)量级,因此,(34)式和(35)式的精度可以达到O(J22·T)量级,基本上与O(T2)等价.类似于(22)式和(23)式,可以得到形如(24)式和(25)式的观测方程:

(36)

(37)

设计矩阵$\hat A$2 和$\hat A$3 的精度较$\tilde{A}$2 和$\tilde{A}$3 更高,但计算也更为繁琐.

求解重力梯度张量不变量观测方程(24)、(25)、(36)和(37),与常见的重力梯度分量的观测方程相同,可以采用离散化的方法直接求解,比如最小二乘方法、最小二乘配置法等,所不同的是计算量是张量不变量,并且与参考坐标系的选择无关.然而,由于在每一个观测点均要建立方程,对于长期的海量观测数据以及求解高阶重力位模型参数,构建法方程耗时巨大,求解困难.

4.4 不变量球面边值问题的求解

由于GOCE 卫星轨道面非常接近于一个球面,因此可以将沿卫星轨道的重力梯度张量不变量延拓到轨道平均球面上,然后将其内插成格网点值.需要指出的是,由于GOCE 卫星设计倾角为96.5°的缘故,在地球两极处还有观测的球冠极区空白区域,目前较为简单的方法就是采用精度较高的重力场模型计算两极地区的重力梯度张量不变量来填补.

由(16)式可得,在球近似下重力梯度张量不变量δI2 的球面边值问题为

(38)

同样,由(1 7)式可得不变量δI3 的球面边值问题为

(39)

上述边值问题与重力位二阶径向导数的边值问题完全一致,可以采用已有的各类调和分析方法[20]解算得到重力场位系数 (CnmSSnmS)及其相应的球近似解TS.

下面进一步讨论不变量δI2δI3 顾及J2 项影响的边值问题的解算.

根据(28)式,在保持精度O(J22·T)的前提下,可以得到顾及J 2 项影响的重力梯度张量不变量δI2的球面边值问题:

(40)

将球近似解(CnmSSnmS)代入(40)式,可得边值条件

(41)

又因为,存在如下关系式[15]:

(42)

(43)

式中,

(44)

(45)

对(43)式求导可得

(46)

运用关系式(42)~(46),由(41)式可得

(47)

因此,设边值问题(40)式的解为(CnmI2SnmI2),则可表示成

(48)

其中,

(49)

对于不变量δI3 边值问题的解算方法与δI2 相同.根据(29)式,可以得到顾及J 2 项影响的重力梯度张量不变量δI3 的球面边值问题:

(50)

将球近似解(CnmSSnmS)代入(50)式,可得边值条件:

(51)

同样,根据(42)~(46)式,由(51)式可得

(52)

因此,顾及J 2 项重力梯度边值问题(50)式的解可表示为

(53)

其中,

(54)

至此,已经获得了顾及J 2 项影响的引力张量不变量I2I3 的球面边值问题的完整解.具体解算过程为:首先采用调和分析方法得到球近似解,然后运用(49)式或(54)式求得顾及J 2 项的系数改正项,最后根据(48)式或(53)式得到边值问题(40)式或(50)式的解.

4.5 GOCE任务中不变量的计算及误差分析

不变量的计算需要地球重力场全张量观测,然而GOCE 任务只能有效给出4 个分量:V11V22V33V13,而V12V23 的观测误差扩大近千倍,因此认为是无效的观测量.下面具体分析解决办法.

由(5)式可知,在地球重力梯度张量(Vij)中,对角线元素是主项,V33 是最大项,V33≈-2740E ,而V11V22≈1370E;而非对角线梯度值量级约为对角线梯度值的10-3.由I2I3 的计算公式(8)和(9)可以看出,非对角项运算结果在I2I3 中只能占到约10-3量级,是小量.因此用已知高精度重力场模型计算V12V23 的值,再代入I2I3 的计算式中,基本可以满足精度要求,但问题是利用模型计算的是在地固坐标系(O-X1X2X3)下的V12,V23,还需要转换成GRF下的V12V23.由(3)式可知,由Vij转换为Vij的计算同样受到卫星姿态误差的影响,在只考虑卫星姿态误差影响的情况下,其误差传播公式为

(55)

(56)

然而,由于非对角线梯度值量级较小,因此大大降低了对卫星姿态的精度要求.但是,当存在更大的姿态误差时,是必须考虑的,因此不变量法还是需要依靠一定的姿态控制精度或者需要全张量梯度测量.另一种计算V12V23 值的有效方法就是迭代法[13],经过2~3次迭代后,即可得到满意的结果.

下面分析不变量δI2δI3 组合生成时的误差,主要包括观测量Vij的观测误差和线性化误差.前文已分析得到在球近似和顾及J 2 项时的线性化误差分别为O(J 2·T)量级和O(J22·T)量级.下面以观测方程(18)和(19)为例,分析由Vij的观测误差引起不变量的误差.

假设GRF中的观测量Vij的误差是vijv2v3分别为观测误差传播给不变量δI2δI3 的误差项,则由(11)式和(12)式,并略去O(vij2)和O(T·vij)等以上的小量,可得

(57)

(58)

顾及(15)式,可得

(59)

(60)

假设V11V22V33 具有相同的观测精度,即v11 =v22=v33,按照误差传播定律,v2v3v11v22v33 的加权组合传播给δI2δI3,并可得δI2δI3 的精度分别为$\frac{{2\mu }}{{{r^3}}}{v_{33}},\frac{{2{\mu ^3}}}{{{r^6}}}{v_{33}}$,由观测方程(18)和(19)可知,观测量Vij的误差并没有使得观测误差放大,后文的数值试验也证实了这一点.

5 模拟试算及结果分析 5.1 重力梯度张量分量和不变量的模拟

本文以δI2 为例进行数值试验.首先采用EGM2008模型(取前360阶次)作为真实的地球重力场V,并采用文献[21]中的球谐综合算法计算得到平均轨道半径为RS=6621.00877km 的全球分布30′×30′分辨率的格网均值VijA(共有259200个值),模拟精度可达10-4mE[21].为分析观测量Vij的误差对生成不变量I2 的影响,在VijA的基础上加入误差为1.0mE的零均值白噪声得到数据VijB,并由(8)式分别计算出不变量I2AI2B .

在模拟试算中采用了两种参考重力场,第一种取水准旋转椭球正常重力场(见(10)式),记为U′,第二种取EGM96模型[22](360阶次)所产生的参考重力场,记为U″,显然后者更加接近于真实的重力场.采用同样的方法模拟出参考不变量I12 refI22 ref .根据(11)式,分别计算出扰动不变量:δI2A1 =I2A -I2 ref1δI2A2 =IA2 -I2 ref2δI2B1 =IB2 -I2 ref1δI2B2 =IB2 -I2 ref2 ,其中δI2A1δI2A2 的等值线图如图 1图 2所示.

图 1 重力梯度张量不变量δI2A1 的等值线图 Fig. 1 Isoline of the gravity gradient tensor invariant
图 2 重力梯度张量不变量δI2A2 的等值线图 Fig. 2 Isoline of the gravity gradient tensor invariant §I?2

采用较差法分析不变量组合生成的误差影响,记

(61)

η′ 和η″ 表示由重力梯度观测量VijB生成δI2 时的误差,其值描述了计算量δI2 的精度.表 1给出了其统计结果.

表 1 η′ 和η″ 的统计结果(单位:mE) Table 1 Statistical result of y and 7 (Unit: mE)

表 1 可知,当重力梯度观测量中含有1 mE的误差时,引起η′和η″ 的标准差均为0.816 mE,由此可见,重力梯度的观测误差并没有被放大.

同样采用较差法分析球近似边界条件模型误差,记

(62)

ε′,ε″ 表示球近似边界条件的误差,其值描述了边界条件(38)式的精度.表 2给出了其统计结果.

表 2 ε′ 和ε″ 的统计结果(单位:mE) Table 2 Statistical resultof r and e" (Unit: mE)

表 2可知,参考重力场越是逼近实际重力场,则对应的球近似边界条件的精度越高,可以肯定,考虑J 2 项影响的边界条件的精度比球近似边界条件的精度更高.

5.2 不变量边值问题的广义轮胎调和分析解

为了有效地验证不同方法的优劣,采用重力场模型系数估值误差阶方差(ErrorDegreeRMS)进行精度评价,其计算公式为

(63)

式中,(Cnm(est)Snm(est))为计算所得的重力场模型系数的估值,(Cnm(EGM)Snm(EGM))为EGM2008 地球重力场模型.

首先不考虑卫星姿态误差,分别对不变量δI2A1δI2A2δI2B1δI2B2 运用广义轮胎调和分析方法[17]求解球近似边值问题(38)式,得到球近似解,然后根据(48)式和(49)式求解得到顾及J 2 项影响的边值问题(40)式的解,它们的误差阶方差统计见图 34.

图 3 张量不变量δI2A1δI2A2 广义轮胎调和分析解的误差阶方差 Fig. 3 Error degree RMS of generalized torus harmonic analysis solution of the tensor invariant δI2A1δI2A2
图 4 张量不变量δI2B1δI2B2 广义轮胎调和 分析解的误差阶方差 Fig. 4 Error degree RMS of generalized torus harmonic analysis solution of the tensor invariant δI2B1δI2B2

图 3、4可知,无论是球近似解还是顾及J 2 项影响的解,采用参考场U″ 比U′ 所得精度要高,这是因为采用更加接近真实重力场的参考模型,所得的边界条件的精度也越高的缘故.另外,顾及J 2 项影响的调和分析解比球近似解的精度要高,这也是由于顾及J 2 项的边界条件精度比球近似要高的缘故.

为了比较分析卫星姿态误差、非全张量观测的影响,本文模拟如下情景:重力梯度观测值V11V22V33V13 的误差为1mE;V12V23 的误差根据卫星姿态误差由(55)式和(56)式估算,当卫星姿态误差分别取10、20、50、100、200 角秒(arcsec)时,V12 的误差分别约为0.1、0.2、0.5、1、2E,V23 的误差分别约为0.14、0.28、0.7、1.4、2.8E;参考重力场取EGM96模型,即U″,针对不同的卫星姿态误差,运用广义轮胎调和分析方法求得顾及J 2 项影响的不变量边值问题的解,其结果如图 5所示.

图 5 重力梯度张量不变量I2 广义轮胎调和分析解的误差阶方差 Fig. 5 Error degree RMS of generalized torus harmonicanalysis solution of the gravity gradient tensor invariant I2

为了与采用重力梯度张量分量(或分量组合)方法的比较,本文以重力梯度张量分量V33 为例,在相同的情形下,运用相同的广义轮胎调和分析方法计算,结果如图 6所示.

图 6 重力梯度分量V33 广义轮胎调和分析解的误差阶方差 Fig. 6 Error degree RMS of generalized torus harmonic analysis solution of the gravity gradient component V33

图 5可知,针对非全张量卫星重力梯度观测,当卫星姿态误差在50角秒内,不变量调和分析解几乎不受姿态误差的影响,当达到100角秒时,调和分析解精度略有下降.而张量分量V33 的调和分析解的精度随着卫星姿态误差的增大,而逐渐降低.比较图 5图 6可知,不存在卫星姿态误差的情况下,张量分量与张量不变量的调和分析解具有相当的精度,但是当存在卫星姿态误差时,张量不变量的调和分析解的精度要明显优于张量分量.不变量δI3 的模拟试验结果几乎与δI2 相同,在此不再详述.

6 结论与展望

本文深入研究了利用重力梯度张量不变量恢复地球重力场的理论与方法.首先在已有的研究基础上,进一步建立了在卫星轨道面上的不变量观测模型,详细完整地推导了重力梯度张量不变量δI2δI3 的球近似和顾及J 2 项影响的球面边值问题的求解公式.针对GOCE 任务的非全张量观测特性,进行了误差分析,理论分析与数值结果表明,不变量δI2δI3的计算虽然会引入复杂的误差传播,但是重力梯度观测误差在δI2δI3 的计算中并没有得到放大,并且该方法能够有效克服非全张量观测的不足.最后运用广义轮胎调和分析方法进行了模拟试验,结果表明,线性化过程中,采用精度较高的参考重力场能够得到更高精度的边界条件,从而得到更高精度的解;另外,顾及J 2 项影响的调和分析解比球近似解的精度更高.不存在卫星姿态误差的情况下,张量分量与张量不变量的调和分析解具有相当的精度,但在卫星姿态控制不理想的情况下,采用重力梯度张量不变量恢复地球重力场的精度更高.

另外,球近似下的不变量观测量δI2δI3 等效于重力梯度径向分量Trr,因此也可尝试用不变量作为观测质量检核的一种手段.

致谢 感谢中国科学院研究生院于锦海教授给予的指导和帮助;感谢审稿专家的修改意见.

参考文献
[1] Bouman J, Rispens S, Koop R. GOCE gravity gradients for use in Earth sciences. In: Proceedings of the 3rd International GOCE User Workshop, 6-8 November 2006, Frascati, Italy, (ESA SP-627, January 2007), 135~139
[2] Visser PNAM, van den IJssel J, van Helleputte T, et al. Rapid and precise science orbit determination for the GOCE satellite. In: Proceedings of the 3rd International GOCE User Workshop, 6–8 November 2006, Frascati, Italy (ESA SP-627, January 2007), 235~239
[3] Gruber T, Rummel R, Koop R. How to use GOCE level 2 products. In: Proceedings of the 3rd International GOCE User Workshop, 6-8 November 2006, Frascati, Italy (ESA SP-627,January 2007),205~211
[4] Bruinsma S, Marty J C, Balmino G. Numerical simulation of the gravity field recovery from GOCE mission data. In: Proceedings of the 2nd International GOCE User Workshop, 8-10 March 2004,Frascati, Italy, (ESA SP-569, June 2004)
[5] Pail R, Schuh W D,Wermuth M. GOCE gravity field processing. In: Jekeli C, Bastos L, Fernandes J (eds) International association of geodesy symposia, "Gravity, Geoid and Spacemissions", vol 129. Springer, Berlin, 2005.36~41
[6] Migliaccio F, Reguzzoni M, Sansò F. Space-wise approach to satellite gravity field determination in the presence of coloured noise. J Geod , 2004, 78: 304-313. DOI:10.1007/s00190-004-0396-z
[7] Pail R. A parametric study on the impact of satellite attitude errors. J Geod , 2005, 79: 231-241. DOI:10.1007/s00190-005-0464-z
[8] Müller. GOCE gradients in various reference frames and their accuracies. Adv Geosci , 2003, 1: 33-38. DOI:10.5194/adgeo-1-33-2003
[9] Rummel R. Satellite gradiometry. In: Sünkel H(ed) Mathematical and Numerical Techniques in Physical Geodesy. Lecture Notes Earth Sciences, Springer, Heidelberg, 1986, 7: 317~363
[10] Holota P. Boundary value problems and invariants of the gravitational tensor in satellite gradiometry. In: Sansò F, Rummel R(eds) Theory of satellite geodesy and gravity field determination, Lecture Notes Earth Sciences, Springer, Heidelberg, 1989, 25:447~457
[11] Sacerdote F F. Sansò. Some problems related to satellite gradiometry. Bull Géod , 1989, 63: 405-415. DOI:10.1007/BF02519638
[12] Vermeer M. Observable quantities in satellite gradiometry. Bull Géod , 1990, 64: 347-361.
[13] Baur O N, Sneeuw E W. Grafarend. Methodology and use of tensor invariants for satellite gravity gradiometry. J Geod , 2008, 82: 279-293. DOI:10.1007/s00190-007-0178-5
[14] Baur O E W. Grafarend. High-performance GOCE gravity field recovery from gravity gradient tensor invariants and kinematic orbit information. In: Flury J, Rummel R, Reigber C, Rothacher M, et al (eds) Observation of the Earth System from Space. Springer, Heidelberg, 2006.239~253
[15] 于锦海, 赵东明. 引力梯度不变量与相关边界条件. 中国科学: 地球科学 , 2010, 40(2): 178–187. Yu J H, Zhao D M. The gravitational gradient tensor's invariants and the related boundary conditions. Sci. in China, Earth Sci. (in Chinese) , 2010, 40(2): 178-187.
[16] Pavilis N K, Holmes S A, Kenyon S C. An Earth Gravitational Mode to Degree 2160: EGM2008, presented at the 2008 General Assembly of the European Geosciences Union, Vienna, Austria, April 13~18, 2008
[17] 吴星, 张传定, 赵东明. 卫星重力梯度分量的广义轮胎调和分析方法. 测绘学报 , 2009, 38(2): 101–107. Wu X, Zhang C D, Zhao D M. Generalized tours harmonic analysis of satellite gravity gradients component. Acta Geodaetica et Cartographica Sinica (in Chinese ) (in Chinese) , 2009, 38(2): 101-107.
[18] 陆仲连. 地球重力场理论与方法. 北京: 解放军出版社, 1996 . Lu Z L. Theories and Methods of the Earth's Gravity Field (in Chinese). Beijing: People's Liberation Army Publishing House, 1996 .
[19] 张传定.卫星重力测量——基础、模型化方法与数据处理算法.郑州:信息工程大学测绘学院,2000. Zhang C D. Satellite gravimetry: foundation, modeling methods, and data processing algorithms (in Chinese), Zhengzhou: Institute of Surveying and Mapping, 2000.
[20] 宁津生, 罗志才, 陈永奇. 卫星重力梯度数据用于精化地球重力场的研究. 中国工程科学 , 2002, 4(7): 23–28. Ning J S, Luo Z C, Chen Y Q. Application of satellite gravity gradiometry data to the refinement of the Earth's gravity field. Engineering Science (in Chinese) , 2002, 4(7): 23-28.
[21] 吴星, 张传定, 叶修松, 等. 卫星重力梯度数据的模拟研究. 测绘科学技术学报 , 2008, 25(6): 391–395. Wu X, Zhang C D, Ye X S, et al. Simulation of satellite gravity gradient tensor observations. Journal of Geomatics Science and Technology (in Chinese) , 2008, 25(6): 391-395.
[22] Lemoine F G, Kenyon S C, Factor J K, et al. The Development of the Joint NASA GSFC and NIMA Geopotential Model EGM96. NASA GSFC, Greenbelt, Maryland, July 1998