卫星重力径向梯度数据的最小二乘配置调和分析

吴 星,张传定,刘晓刚

1.总装备部工程设计研究总院,北京100028;2.海洋测绘研究所,天津300061;3.信息工程大学测绘学院,河南郑州450052

卫星重力径向梯度数据的最小二乘配置调和分析

吴 星1,2,张传定3,刘晓刚3

1.总装备部工程设计研究总院,北京100028;2.海洋测绘研究所,天津300061;3.信息工程大学测绘学院,河南郑州450052

研究利用卫星重力梯度径向分量确定地球引力场位系数的最小二乘配置(LSC)调和分析方法。首先论述最小二乘配置法的原理,推导扰动引力梯度观测量与球谐系数之间的协方差和自协方差矩阵。在扰动引力梯度观测数据为等经差规则网格数据的情况下,引力位与扰动引力梯度之间的协方差矩阵具有分块Toeplitz循环阵的结构,有效地利用FFT变换技术将其降阶。研究利用截断奇异值分解法(TSVD)减弱协方差阵的病态性。最后得到引力梯度径向分量的最小二乘配置调和分析的完整计算公式。模拟试算结果表明,基于TSVD的最小二乘配置调和分析方法,能够以较高的精度还原全球重力场,验证本文算法的有效性和实用性。

地球重力场模型;最小二乘配置;卫星重力梯度;调和分析;协方差阵;Toeplitz矩阵;TSVD

1 引 言

国际上首颗重力梯度测量卫星 GOCE已经于2009年3月17日成功发射,其科学目标是在100 km的空间分辨率下,实现大地水准面精度达厘米级、重力异常精度达1～2 mGal的地球重力场(1 mGal=10-3cm/s2)。如何利用重力梯度观测数据恢复地球重力场早已得到了广泛研究,目前主要有两大类方法:时域法和空域法,两类方法各有优缺点[1-2]。空域法中的最小二乘配置法作为确定全球或区域地球重力场的主要方法之一备受关注,其理论与方法也得到了完善和补充[3-4]。罗志才提出的频域最小二乘配置法可以有效地提高计算速度[5],张传定提出的最小二乘复数配置法解决了重力梯度张量水平分量的最小二乘配置问题[6]。Tscherning和Sanso研究了能够处理各类重力径向分量的最小二乘配置法和快速配置法[7-8];汪海洪等针对不同分辨率的数据融合问题,提出多分辨率最小二乘配置法[9];Arabelos等和Balmino分别研究了误差协方差矩阵的计算问题[10-11];Kotsakis提出具有协方差匹配约束的最小二乘配置法[12];Reguzzoni等提出最优多步配置法[13]。这些研究工作都是利用了最小二乘配置法能综合各类观测数据确定地球重力场,并且在数据处理中能顾及观测量的误差和计算待估量的误差协方差矩阵的特点。但是,配置法要求有可靠的先验协方差模型,需要求解协方差矩阵的逆矩阵。在局部重力场逼近中,观测量协方差矩阵维数不高,可直接求解,因此最小二乘配置在局部重力场逼近得到了广泛应用。在全球重力场逼近中由于协方差矩阵维数较大,必须对其进行合理分解,才能给出稳定解。

为此,本文主要研究卫星重力径向梯度规则网格观测数据的最小二乘配置调和分析的算法,推导引力位系数与格网平均扰动引力梯度间的协方差矩阵,研究协方差矩阵降阶求逆方法及其数值稳定性,给出仿真试验结果。

2 解算重力位系数的最小二乘配置法

卫星重力梯度测量可以利用卫星上的重力梯度仪直接测定重力位梯度张量的所有分量或部分分量。按照物理大地测量的传统作法,首先选择一合适的参考重力场,将观测量转换为扰动重力位梯度张量,然后将观测数据归算到以地球质心为球心,卫星标称球带空间的平均球面上,并进行网格化处理,将得到的球面规则网格平均扰动引力梯度张量数据作为已测信号,利用最小二乘配置法求解扰动重力位系数。

2.1 最小二乘配置解法的原理

为了简便,假设不含非随机参数向量,此时最小二乘配置的数学模型简化为[14]

式中,l为观测值向量;t是已测信号向量;s为待估信号向量;I为单位阵;v为随机的观测误差向量。

随机模型为

其中,Ctt为已测信号协方差阵;Cvv为观测误差协方差阵。

利用最小二乘配置准则,可得未测信号s的推估值和误差协方差阵分别为

其中,Cst为已测信号和待估信号的互协方差阵。

显然,对于由卫星重力梯度数据解算地球重力场重力位系数的情况下,已测信号为卫星重力梯度观测信息,待估信号为所求的重力位系数。由此可见,利用最小二乘配置法求解位系数的关键问题是协方矩阵与误差矩阵之和的求逆问题。特别是已测信号的数据量较大时,逆矩阵的求解方法及其稳定性将直接影响配置结果的可靠性。

2.2 二阶径向扰动重力梯度的级数展开式

地球外部任意一点 P(r,θ,λ)处扰动位 T(r, θ,λ)的球谐函数级数表达式为[6]

当且仅当是规则球面点值或均值时,协方差矩阵才具有可分解的结构。为此,采用习用球面划分[15-16],即把整个球面用经、纬线划分成经度差为Δ θ、纬度差为Δ λ的 N×Nl个网格,其中, N、Nl分别为纬度方向和经度方向的网格分化数,Δ θ=π/N,Δ λ=2π/Nl。同时假设卫星轨道面上的引力梯度数据经过必要的归算后得到全球覆盖(非极轨卫星,两极空白地区用模型填充)的网格平均引力梯度观测值¯fi,k(i、k代表该网格所在的位置,i=0,1,…,N-1;k=0,1,…,Nl-1)。则有

式中,函数sinc(x)的表达式为

并且可以通过递推公式计算[6,17]。

2.3 扰动引力位系数与扰动引力梯度的协方差

在研究地球重力场有关问题时,扰动位 T之间的协方差是最基本的,其他扰动场元之间的协方差函数都可以由它求得。设 T(P)和 T(Q)是空间两点 P(r,θ,λ)和Q(r′,θ′,λ′)的扰动位,ψ为P、Q两点球面角距,则扰动位空间协方差函数可表示为[14]

式中,Pn(cosψ)为勒让德多项式;符号“＊”表示复数共轭;Kn为扰动位阶方差,一般可采用参考引力位模型系数计算,计算式为

另外,可用重力异常阶方差模型求得,其关系式为

例如采用如下6参数重力异常阶方差模型—T/R模型[18]

其中,α1、α2、S1、S2为实常数;A、B为整常数。取α1=3.405,α2=140.03,S1=0.998 006,S2= 0.914 232,A=1,B=2。α1、α2的量纲是(mGal)2(1 mGal=10-3cm/s2),其他量是无量纲的数。将式(12)代入式(11)并代入式(10)后,完成级数求和,即得全球协方差的解析函数。

设Si为任一重力场元信号与扰动位的线性泛函算子,则信号SiT(P)和 SjT(Q)间的协方差函数可由扰动位的协方差函数表示[4],即有

根据球谐函数正交性,由式(5)可得

把式(9)代入式(16),并且令 dnm=Kn/((2n+1)· (2-δ0m)),则可得

这说明了扰动引力位系数的协方差是一对角阵。

由协方差函数传播公式(13)和式(17),可得扰动引力位系数与格网平均扰动引力梯度径向分量(Q)间互协方差为

格网平均扰动引力梯度径向分量间的自协方差为

从式(19)可以看出,对于等经差规则网格数据而言,其第i个纬度圈与第i′个纬度圈间的协方差子矩阵是一个T oeplitz循环阵,利用傅里叶变换阵可将其化算为对角阵,降低矩阵求逆的维数,提高计算结果的稳定性和可靠性。正是利用这一特性,才使全球协方差矩阵的稳定求逆成为可能。

2.4 最小二乘配置调和分析解

将规则扰动引力梯度径向分量观测数据作以下排列

其中,¯fi是第i(i=0,1,2,…,N-1)个纬度圈上所有观测量所组成行向量

则协方差矩阵Ctt可分解成N2个 Nl×Nl维的子矩阵,用表示(i,i′=0,1,…,N-1),其元素是所有第i行和第i′行平均引力梯度之间协方差函数值。另外,对于近圆极轨道卫星,由于同一纬度圈上,每个子块面积相同,卫星重力梯度测量中一个完整覆盖周期内落入子块内的测量点数大致相同,子块平均值的误差方差大致相同[6]。第 i个纬度带的误差方差用表示,则误差方差阵 Cvv同样可分解成N2个 Nl×Nl维的子矩阵,用表示(i,i′=0,1,…,N-1),并有

因此由(3)式可知,最小二乘配置调和分析解为

其中,

如此,对于 N2个子对角阵Dii′,若追踪 m指标,则共有Nl个N×N维的满秩矩阵

对这 Nl个Rm求逆,再按照对应项完成矩阵的乘法运算后,得到位系数的估值为

相应地,由式(4)、式(22)与式(25),得位系数估值误差的方差为

2.5 矩阵TSVD正则化求解

尽管利用分块循环阵解决了协方差阵的降阶求逆,但由于两极地区子午线收敛问题,即极区网格之间的球面角距ψ变化范围较小,将导致矩阵Rm存在病态性问题。实际数据表明,对于50′× 50′或40′×40′分辨率网格,矩阵 Rm的奇异值数值非常小,并接近于阶梯型分布。为了解算该类病态方程,本文采用截断奇异值分解(TSVD)正则化方法[21]。该方法是把容易造成不稳定的较小的奇异值直接截去,使原来的不适定问题转化为一个适定问题来求解。

设矩阵Rm的奇异值分解式(SVD)为

其中,正交矩阵U=(u0,u1,…,uN-1),V=(v0,v1,…,vN-1),Σ=diag(σ0,…,σN-1),奇异值满足σ0≥…≥σN-1>0。

则式(32)可表示为

由式(35),可得最小二乘配置的 TSVD正则化解为

其中,整数k称为截断参数,也就是正则化参数。

正则化参数的选取通常有容差原理法、广义交叉验证法,还有L曲线法等[22]。容差原理法需要知道数据项中噪声大小,而实际问题中噪声水平是未知的;广义交叉验证法操作较为烦琐,且结果不稳定;L曲线法的关键是定位L曲线的最大曲率的那个点,选择所对应的k值作为截断参数。L曲线的拐角点所对应的解不但平衡了解范数和残差范数,而且趋于平衡正则化误差和扰动误差。本文采用L曲线法确定正则化参数,具体方法详见[21,23]。

3 仿真试算与分析

3.1 协方差函数的计算

全球协方差函数模型采用两种方法计算,第一种方法是采用当前最新、精度最高的地球重力场模型EGM2008计算[23]。EGM2008模型是美国国家地理空间情报局在充分利用包括地面重力、卫星测高、卫星重力等最新数据基础上研制出的新一代地球重力场模型。GPS水准点外部检核结果表明,EGM2008模型具有很高的精度,在全球范围内标准差达到13 cm[23],在我国高程异常的总体精度也高达20 cm[24]。第二种方法是采用T/R全球重力异常阶方差模型转换得到。扰动位的阶方差统计如图1所示。

图1 扰动位阶方差统计图Fig.1 Chart of potential degree-variance

3.2 扰动引力梯度径向分量数据的仿真计算

利用截至300阶次的 EGM2008模型,采用文献[25]中球谐综合方法分别模拟了全球分辨率分别为40′×40′和50′×50′的格网扰动引力径向梯度径向分量数据,卫星轨道高度设为250 km。如表1所示,共产生了4套数据,即数据1～4。数据1、2是由2～300阶模型生成的,而数据3、4模拟生成时,截掉了低阶部分(2-50阶)的影响。由文献[25]可知,数据 A1～A4的精度优于10-7mE(1 mE=10-12s-2),其误差可以忽略不计。为了使模拟观测值具有实测性,需要向模拟梯度观测值中加入观测误差。根据卫星重力梯度的观测精度、采样间隔、飞行时间等不同,可得不同精度的重力梯度格网平均值,计算公式为[6]

其中,NT是整个飞行期的采样点数;Nik为落入第(i,k)网格的点数。

为简单起见,本文假设得到扰动引力梯度张量模拟观测值是等精度观测的,而且全球均匀分布,假定得到的40′×40′格网平均值的精度为1.0 mE,则根据式(37)可知,和50′×50′格网平均值的精度为0.84 mE,故数据B是在数据A的基础上分别均加入了标准差为1.0 mE和0.84 mE的零均值白噪声而产生的。

表1 EGM2008模拟的卫星重力径向梯度数据Tab.1 Radialsatellite gravity gradientssimulated using EGM2008

3.3 调和仿真分析计算

以计算位系数误差的阶均方根(error degree RMS)

为统计量对最小二乘配置调和分析方法进行评估分析。式(38)中,(,)为位系数估值, (,)表示EGM2008模型位系数。

由数据1～4采用本文的最小二乘配置调和分析法(LSCHA),协方差函数利用截至360阶次的EGM2008模型计算,所得统计结果见图2。

图2 最小二乘配置调和分析结果Fig.2 Least-squares collocation harmonic analysis of radial gravity gradients

在相同数据的情况下,采用截断至360阶的T/R全球重力异常阶方差模型转换得全球协方差函数,运用LSCHA方法,所得结果与图1几乎完全一致。分析图2可知:

1.对于截掉低阶部分影响的数据3、4的计算结果要优于数据1、2的计算结果。产生这一差异的原因是在数据1、2中考虑了重力场的低频部分,由于低频部分与高频部分量级相差较大,这影响了高阶位系数的计算结果。重力梯度数据反映了重力场的中频特征,不能分辨重力场的低频部分。

2.分辨率为40′×40′全球重力梯度数据比50′×50′的数据的调和分析结果有显著的提高。主要原因是实际地球重力场是连续分布的,当采用一定分辨率的离散格网数据表示时,存在数据的离散化误差。因此采用更高分辨率的观测数据表示,可以减小离散化误差,提高调和分析精度。

3.即使是采用精度较高的数据 A1～A4,其最小二乘配置调和分析结果仍然均与EGM2008模型存在差异,即不能100%地一致恢复位系数,这主要是由下列因素引起的:①由于协方差矩阵存在病态性问题,其求逆势必存在一定的计算误差,即使采用正则化方法也不能完全避免;②需要采用更高分辨率的数据,以进一步减小离散化误差;③配置方法是重力场数据处理的重要手段,但它毕竟是统计方法,能给出信号的估值,但未必能给出其真值,存在误差是容许的。

在数据相同的情况下,采用了文献[27]中的轮胎调和分析方法(THA)进行了仿真试算,结果如图3所示。

图3 轮胎调和分析结果Fig.3 Torus harmonic analysis of radial gravity gradients

比较图2和图3可知,本文提出的基于 TS-VD的最小二乘配置调和分析结果略优于轮胎调和分析结果,尤其是在观测数据精度较高的情况下(数据A1～A4),更为明显,而当采用数据B1～B4时,两种方法所得结果几乎一致,但轮胎调和分析的计算速度更快。

4 结 论

通过本文的理论推导和数值试验分析可得如下结论:

1.在重力梯度数据满足一定条件(等经差网格分布,以及同纬度带具有相同误差)的情况下,其协方差子矩阵是 Toeplitz循环矩阵,并可利用傅里叶变换阵将其化算为对角阵,降低矩阵求逆的维数,从而提高计算结果的稳定性和可靠性。正是这种特性才使得超大型矩阵的求逆成为可能,从而使最小二乘配置法真正用于恢复全球重力场。

2.由EGM2008模型和 T/R模型计算得到的扰动位阶方差在低阶部分吻合较好。

3.本文基于 TSVD的最小二乘配置调和分析方法计算精度略优于轮胎调和分析方法,但计算速度比轮胎调和分析方法慢。

4.由于两极地区子午线收敛问题,导致了全球协方差矩阵病态性问题,采用截断奇异值分解正则化方法能够有效地减弱法矩阵的病态性。

[1] KOOP R.Global Gravity Field Modelling Using Satellite Gravity Gradiometry[R]. Amsterdam: Netherlands Geodetic Commission,1992.

[2] RUMMEL R,VAN GELDEREN M,KOOP R,et al. Spherical Harmonic Analysis of Satellite Gradiometry[R]. Amsterdam:Netherlands Geodetic Commission,1993.

[3] KRARUP T.A Contribution to the Foundation of Physical Geodesy[R].Copenhagen:Danish Geodetic Institute,1969.

[4] MORTIZ H.Advanced Physical Geodesy[M].England: Abacus Press,1980.

[5] LUO Zhicai.Theories and Method of the Determination of the Earth’s Gravity Field Using Satellite Gravity Gradient Data[D].Wuhan:Wuhan Technical University of Surveying and Mapping,1996.(罗志才.利用卫星重力梯度数据确定地球重力场的理论与方法[D].武汉:武汉测绘科技大学,1996.)

[6] ZHANG Chuanding.Satellite-borne Gravity Measurement: Foundation,Modeling Methods,and Data Processing Algorithms[D].Zhengzhou:Information Engineering University,2000.(张传定.卫星重力测量:基础、模型化方法与数据处理算法[D].郑州:信息工程大学,2000.)

[7] TSCHERNING C C.Computation of Spherical Harmonic Coefficients and Their Error Estimate Using Least-squares Collocation[J].Journal of Geodesy,2001,75(1):12-18.

[8] SANSO F,TSCHERNING C C.Fast Spherical Collocation:Theory and Examples[J].Journal of Geodesy,2003, 77(1-2):101-112.

[9] WANG Haihong,LUO Zhicai,LUO Jia,et al.Preliminary Research on Multiresolution Least Square Collocation [J].Journal of Geodesy and Geodynamics,2006,26(1): 115-118.(汪海洪,罗志才,罗佳,等.多分辨率最小二乘配置法初探[J].大地测量与地球动力学,2006,26(1): 115-118.)

[10] ARABELOS D N,TSCHERNING C C.Error-covariances of the Estimates of Spherical Harmonic Coefficients Computed by LSC,Using Second-order Radial Derivative Functional Associated with Realistic GOCE Orbits[J]. Journal of Geodesy,2009,83(5):419-430.

[11] BALMINO G.Efficient Propagation of Error Covariance Matrices of Gravitational Models:Application to GRACE and GOCE[J].Journal ofGeodesy,2009,83(10): 989-995.

[12] KOTSAKIS C.Least-squares Collocation with Covariance-matching Constrains[J].Journal of Goedesy,2007, 81(10):661-677.

[13] REGUZZIONI M,TSELFES N.Optimal Multi-step Collocation:Application to the Space-wise Approach for GOCE Data Analysis[J].Journal of Goedesy,2009,83 (1):13-29.

[14] HEISKANEN WA,MORTIZ H.Physical Geodesy[M]. San Francisco:WH Freeman,1967.

[15] COLOMBO O L.Optimal Estimation from Data Regularly Sampled on a Sphere with Applications in Geodesy:Rep 291[R].Columbus:Ohio State University,1979.

[16] COLOMBO O L.Numerical Methods for Harmonic Analysison the Sphere[R]. Columbus:Ohio State University,1981.

[17] WU Xing.Research of Methods of Spherical Harmonic Analysis of the Earth’s Gravity Field[D].Zhengzhou: Information Engineering University,2005.(吴星.地球重力场调和分析方法研究[D].郑州:信息工程大学,2005.)

[18] TSCHERNING C C,RAPP R H.Closed Covariance Expressions for Gravity Anomalies,Geoid Undulations, and Deflections of the VerticalImplied by Anomaly Degree-variance Models:Rep 208[R].Columbus:Ohio State University,1974.

[19] XU Zhong,ZHANG Kaiyuan,LU Quan.Fast Arithmetic of Toeplitz Band Matrices[M].Xi’an:Northwestern Polytechnical University Press,2001:31-33.(徐仲,张凯院,陆全.TOEPLITZ矩阵类的快速算法[M].西安:西北工业大学出版社,2001:31-33.)

[20] HANSEN P C.The Truncated SVD as a Method for Regularization:NA-M-86-36[R].Stanford:Stanford University,1986,27:534-553.

[21] WANG Zhenjie.Regularization of Ill-posed Problems in Surveying[M].Beijing:Science Press,2006.(王振杰.测量中不适定问题的正则化解法[M].北京:科学出版社,2006.)

[22] HANSEN P C.Analysis of Discrete Ill-posed Problems by Means of L-curve[J].SIAM Review,1992,34(4): 561-580.

[23] PAVILIS N K,HOLMES S A,KENYON S C,et al.An Earth Gravitational Model to Degree 2160:EGM2008 [R].Vienna:2008 General Assembly of the European Geosciences Union,2008.

[24] ZHANG Chuanyin,GUO Chunxi,CHEN Junyong. EGM2008 and ItsApplication Analysis in Chinese Mainland[J].Acta Geodaetica et Cartographica Sinica, 2009,38(4):283-289.(章传银,郭春喜,陈俊勇. EGM2008地球重力场模型在中国大陆适用性分析[J].测绘学报,2009,38(4):283-289.)

[25] WU Xing,ZHANG Chuanding,YE Xiusong,et al. Simulation of Satellite Gravity Gradient Tensor Observations[J].Journal of Geomatics Science and Technology, 2008,25(6):391-395.(吴星,张传定,叶修松,等.卫星重力梯度数据的模拟研究[J].测绘科学技术学报, 2008,25(6):391-395.)

[26] WUXing,ZHANG Chuanding,ZHAO Dongming. Generalized T orus Harmonic Analysis of Satellite Gravity Gradients Component[J].Acta Geodaetica et Cartographica Sinica,2009,38(2):101-107.(吴星,张传定,赵东明.卫星重力梯度分量的广义轮胎调和分析方法[J].测绘学报,2009,38(2):101-107.)

(责任编辑:丛树平)

Least-squares Collocation Harmonic Analysis of the Radial Satellite Gravity Gradients

WU Xing1,2,ZHANG Chuanding3,LIU Xiaogang3
1.Center for Engineering Design and Research under the Headquarters of General Equipment,Beijing 100028,China;2.Institute of Hydrographic Surveying and Charting;Tianjin 300061,China;3.Institute of Surveying and Mapping,Information Engineering University,Zhengzhou 450052,China

The least-squares collocation harmonic analysis method,which is used to determine the earth geopotential coefficients from the radial satellite gravity gradient,is deeply studied.The principle of the least-squares collocation is firstly dissertated,and the covariance and self-covariance matrix between the disturbing gravity gradients and spherical harmonics are derived.When the disturbing gravity gradients are in the regularized equi-longitude grid,the covariance matrix between geopotential and disturbing gravity gradients has the configuration of blocked Toeplitz circulation matrix,and its degree can be lowered by using the fast Fourier transform technology effectively.Truncated singular value decomposition(TSVD)which is used to solve the ill-posed problem of covariance matrix is studied.The complete formula of least-squares collocation harmonic analysis of the radial gravity gradient is finally obtained.The simulation experiment results show that the least-squares collocation harmonic analysis based on TSVD can recover the global gravity field in a rather high precision and validity and practicability of algorithms in this paper are also testified.

Earth gravity field model;least-squares collocation;satellite gravity gradient;harmonic analysis; covariance matrix;Toeplitz matrix;TSVD

WU Xing(1979—),male,PhD,engineer, majors in Geodesy and Surveying Engineering.

E-mail:wuxing1979@163.com

1001-1595(2010)05-0471-07

P223

A

国家自然科学基金(40774031);全国优秀博士论文专项基金(200344);信息工程大学博士生创新基金(200707);中科院动力大地测量学重点实验室开放基金(L06-06)。

2009-07-16

2010-03-02

吴星(1979—),男,博士,工程师,研究方向为大地测量学与测量工程。