Tikhonov正则化方法在GOCE重力场求解中的模拟研究

徐新禹,李建成,王正涛,邹贤才

1.武汉大学测绘学院,湖北武汉430079;2.武汉大学地球空间环境与大地测量教育部重点实验室,湖北武汉430079

Tikhonov正则化方法在GOCE重力场求解中的模拟研究

徐新禹1,2,李建成1,2,王正涛1,2,邹贤才1,2

1.武汉大学测绘学院,湖北武汉430079;2.武汉大学地球空间环境与大地测量教育部重点实验室,湖北武汉430079

给出四类可用于重力场解算的正则化矩阵(零次、一次、二次和Kaula),以及用于确定正则化参数的L曲线法和GCV方法的数学模型。基于SA方法,利用模拟数据分析讨论了零次、一次以及Kaula正则化矩阵应用于GOCE全球重力场模型确定的有效性,并由Kaula正则化矩阵分析L曲线法和GCV方法确定正则化参数的可行性。数值结果表明三类正则化矩阵获得的最优解(以大地水准面MSE最小为准则确定)的精度水平相近,关键在于相应正则化参数的确定,数值结果同时说明了GCV方法和L曲线法可用于确定正则化参数,且前者较后者具有更好的稳定性。

Kaula;正则化;GOCE;GCV;L曲线

1 引 言

卫星大地测量领域的线性问题往往是不适定的(病态的)[1-3],对卫星重力梯度测量任务来说,导致病态的主要原因可归结为以下几个方面:第一,重力信号向下延拓通常不稳定,该过程不仅放大观测误差,同时也会放大模型误差(未能用模型描述的信号);第二,观测数据的不规则分布破坏了原有球谐函数的正交性,例如 GOCE(gravity field and steady-state ocean circulation explorer)任务的PG(polar gap)问题和观测值间断等问题都使得观测值不能全球均匀覆盖;第三,重力观测量本身不能完全包含所有频谱的重力场信息,例如重力梯度张量观测值的某个分量仅对重力场模型的某些频段比较敏感,利用其恢复全频谱的重力场模型往往导致不适定;第四,卫星重力梯度测量的有色噪声特性也会使法方程系数阵的条件数增大,例如重力梯度仪并不是等精度获得全频段的重力信息。GOCE卫星已于2009年3月17日成功发射 (http:∥www.esa.int/SPECIALS/ GOCE),该任务的目标是实现厘米级精度的大地水准面,相应重力异常的精度为 1～2 mGal (1 mGal=10-3cm/s2),同时空间分辨率达到100 km。为实现该目标,研究利用其观测值恢复重力场的相关理论与方法至关重要,存在许多问题有待于进一步研究和讨论,其中病态问题就是一个需要深入研究的课题。

在物理大地测量领域,许多学者对病态问题都曾进行了深入细致的研究[4-9],提出很多计算正则化解和确定最优正则化参数的方法,在利用卫星重力观测值(卫星跟踪卫星和卫星重力梯度)求解重力场方面,主要有最小二乘配置法、Tikhonov正则化方法以及有偏估计、基于奇异值分解的方法等,正则化参数的选择准则主要有均方误差最小、基本迹估计器、Morozov不一致原理、方差分量估计等。本文主要研究物理大地测量领域广泛采用的 Tikhonov正则化方法在GOCE全球重力场确定中的应用,利用模拟观测值基于半解析法(semi-analytical,SA)讨论不同正则化矩阵的优劣和最优正则化参数的确定。正则化矩阵主要讨论零次、一次、二次和 Kaula正则化矩阵;正则化参数的选择主要讨论L曲线法[10]、广义交叉检验(generalized cross-validation,GCV)[11]。

2 基本原理

2.1 Tikhonov正则化方法

Tikhonov正则化方法是在经典最小二乘准则的基础上,加入参数向量范数最小的条件,同时引入正则化矩阵和正则化参数,对于 Gauss-Markov线性模型有准则如下[12]:

其中,α是正则化参数;K为正则化矩阵,其为正定对称矩阵 K=UTU;x为待求参数,y为观测值; A为观测方程的设计矩阵;P=,其中为观测值误差协方差矩阵为单位权方差。由该准则可求得正则化解如下:

式中,若正则化矩阵 K=I(为单位阵),则Tikhonov正则化就变为有偏估计;如果正则化矩阵的对角线元素取重力场信号的Kaula阶方差模型的倒数,其他元素为零,就得到著名的Kaula正则化方法。

式(1)中右边第一项(残差平方和最小)要求模型与观测值有较好的一致性,第二项要求信号的范数最小,信号的范数最小可提高求解的稳定性,同时考虑两项则需要找到这两种标准的平衡点,不仅使模型与观测值有较好的一致性,同时信号的范数较小,因而需要选择合适的正则化矩阵和确定最优正则化参数,下面介绍几种物理大地测量领域常用的正则化矩阵和正则化参数的确定方法[9]。

2.1.1 零次正则化矩阵

假设式(1)中正则化矩阵 K为单位阵,则右边第二项球谐系数的范数可表示为引力位函数的平方在半径为 R的球面上的积分,进一步推导可得

其中,I为单位矩阵。需要指出,式中的常数因子4π(GM)2在本文的研究中是将其包含在正则化参数中,后面两个正则化矩阵也是将常数部分包括到正则化参数中。

2.1.2 一次正则化矩阵

根据零次正则化矩阵的推导过程,若将引力位的一阶导数作为目标函数,取引力位面梯度的平方函数在半径为R的球面上的积分,则有

K中的元素kij满足kij=(δijn(i)(n(i)+1),δij为克罗内克符号;n(i)表示第i行对应的阶数。

2.1.3 二次正则化矩阵

类似于一次正则化矩阵,将面Laplace算子作用于引力位函数,并在半径为 R的球面上积分可得二次正则化矩阵表达式如下:

K中的元素kij满足kij=δijn2(i)(n(i)+1)2。

2.1.4 K aula正则化

正则化矩阵同样可基于阶方差模型,既可是经验性的重力信号阶方差,也可基于现有的卫星重力模型,根据 Kaula规则[13],重力场信号阶次方差满足

将式(6)取倒数即得 Kaula正则化矩阵,将常数因子1010包括到正则化参数α中,则可得到正则化矩阵元素的形式为 kij=δijn4(i),式中符号的含义同上。从kij的形式可看出其与二次正则化矩阵比较接近,其高阶趋近相同。另外,从数值上分析,正则化矩阵中相应低阶次的元素乘以正则化参数的量级远小于矩阵ATPA中的元素,对求解过程的影响基本相同,因而可认为 Kaula正则化矩阵与二次正则化矩阵是一致的。

2.2 最优正则化参数的确定方法

2.2.1 L曲线法

L曲线法的基本思路是将参数的信号范数η(α)=‖xα‖K与平差残差的范数ρ(α)=‖Axα -y‖P分别取常用对数^η(α)=logη(α),^ρ(α)= logρ(α),并将其绘制成二维曲线图,经证明对于不连续的病态问题,该曲线通常成L形,该方法定义最优正则化参数对应曲线的最大曲率点(即L曲线的顶点),曲线的曲率可表示为

2.2.2 GCV方法

GCV方法选择正则化参数的基本思想是:首先将某一个观测值yk从观测值序列中除去,并利用剩下的观测值求解某个正则化参数对应的解x[k]α,再用该解估计被剔除的观测值,求其与真实观测值之差Δyk=(Ax[k]α)k-yk,如果某个正则化参数能够使所有的Δyk都较小,则认为该正则化参数是最优的,最优选择准则如下[7]:

其中,Qα=A(ATPA+αK)-1ATP为影响矩阵,满足Axα=Qαy;“tr”表示求矩阵的迹;n表示观测值个数。考虑到Qα的计算中需要对法方程矩阵乘积和求逆,且Qα为n×n维矩阵,因此式(8)不适合大量观测数据的情况,为了便于计算,将式(8)进一步变换,由等式tr(AB)=tr(BA)可得

其中,Tα=tr(U(ATPA+αK)-1UT),K=UTU,Tα为u×u的矩阵;u为未知参数个数。

3 数值结果与分析

为分析比较不同正则化矩阵和正则化参数确定方法在 Tikhonov正则化过程中的作用和效果,本文模拟了29 d的 GOCE观测数据,数据采样间隔为5 s,选择EGM96为轨道模拟的参考模型。为使轨道模拟更具真实性,积分时模型的最高阶取至250,轨道倾角约为96.7°,是一个近似重复轨道;重力梯度张量观测值模拟时选择EIGEN_CG03C模型,取至180阶。本文下面的数值分析仅以重力梯度的径向分量为例,在观测值中加入根据GOCE任务解析形式误差功率谱密度基于自回归(auto regressive,AR)模型模拟的有色噪声时间序列[14]。

首先分析比较几类Tikhonov正则化矩阵对求解结果的影响。第二节中已经指明 Kaula正则化矩阵与二次正则化矩阵基本相同,因此这里仅比较零次、一次和 Kaula正则化矩阵。由模拟观测数据,利用零次、一次和Kaula正则化矩阵,分别选择不同的正则化参数(在某个区间内变化),基于半解析法(semi-analytical,SA)方法分别恢复180阶重力场模型。为了评价解算模型的精度,由其计算纬度-80°×80°范围内0.5°×0.5°大地水准面格网点值,并与“真”实重力场(这里指EIGEN_CG03C)对应值求差,计算大地水准面误差RMS。图1分别给出了零次、一次和Kaula正则化矩阵对应大地水准面误差RMS随正则化参数变化的曲线。

图1 零次、一次和Kaula正则化矩阵对应的大地水准面误差RMS随正则化参数变化的曲线Fig.1 The curve of the geoid error RMS with the regularization parameters based on zero-order, first-order and Kaula regularization matrix

图1中横坐标表示Kaula正则化矩阵对应的正则化参数α,零次和一次正则化矩阵对应的参数α需要分别乘以106和103,其中三角形符号表示大地水准面误差RMS最小值对应的点,Kaula正则化对应最小大地水准面误差RMS的正则化参数为α=1010.7,圆点表示采用GCV方法确定的最优正则化参数所对应的点。从图中可看出,基于这三类正则化矩阵求得的模型大地水准面误差RMS随正则化参数变化的曲线具有相同的趋势,均随着正则化参数的增大,先减小后增大;从不同曲线对应大地水准面误差最小值的量级来看,三类正则化矩阵均能较好达到稳定求解目的,特别是一次和Kaula正则化矩阵,若能选择最优的正则化参数,可使大地面误差从几十厘米降到10 cm左右。从图中还可看出,采用GCV方法确定最优参数时,一次正则化矩阵获得的参数值最接近大地水准面误差最小的点(对应的正则化参数可看作是最优正则化参数),Kaula正则化方法确定的正则化参数对应大地水准面误差RMS与一次正则化矩阵的最优值基本一致。此外,当正则化参数小于最优正则化参数时,Kaula正则化矩阵对应解的大地水准面误差变化较慢,而正则化参数较大时,Kaula正则化矩阵对应解的大地水准面误差变化较快,这说明在SA方法中,Kaula正则化矩阵对较大正则化参数的选择较零次和一次正则化矩阵更为敏感。

图1中误差曲线随正则化参数的变化趋势说明,只有选择合适的正则化参数才能使大地水准面误差RMS最小。下面以 Kaula正则化矩阵分析不同正则化参数对求解结果的影响。从理论分析可知,正则化过程可看作是一个滤波过程,小正则化参数很难使求解结果趋于稳定,因而导致求解结果误差较大,而大正则化参数会使求解过程过于稳定,即求解模型将会受到正则化带来的过强“平滑效应”约束,从而使解过于平滑,尤其是高频信号。图2中给出了正则化参数分别选择α=108、α=1010.7、α=1012时的大地水准面的误差。从图中可以看出,α=108时大地水准面误差在近两极地区明显比α=1010.7时大,其他地区则差异不大,说明小的正则化参数对解的稳定性提高有限,仍然暴露出PG问题对求解结果的影响;当α=1012时,大地水准面在重力场高频信号较明显的地区(例如高山地区)的误差较大,反映出高频信号被过度平滑,同时大的正则化参数给接近两极的地区也带来较大的误差;当α=1010.7时,求解的结果在比较的范围内大地水准面误差均较小。

图2 不同正则化参数下求得模型的大地水准面误差RMSFig.2 The geoid error RMS of the models estimated with differentα

除了从大地水准面误差图中可看出不同正则化参数对重力场模型不同频段的影响外,从模型的误差谱图以及模型的阶误差RMS也可看出相同的特点。如图3和图4所示,这里误差谱图是将系数误差的绝对值取常用对数得到的值,n、m分别表示阶次。图中α=108对应解的低次部分系数误差明显较大,α=1012时,图3中(右图)对应解的高频部分的误差普遍比α=1010.7偏大,图4中模型的阶误差RMS在高于142阶时比α=108对应模型的阶误差RMS大,这些都说明大的正则化参数对重力场模型的求解进行了强平滑滤波,从而损失了高频信号,这与图2(c)中的结果一致。

图3 正则化参数不同时求得模型的误差谱Fig.3 The error spectrum of the models estimated with differentα

图4 α=108、α=1010.7、α=1012时求得模型的阶误差RMSFig.4 The degree error RMS of the models estimated with differentα

前面分析的结果均说明选择合适的正则化参数对获得较好求解结果的重要性,在数值模拟分析时可用大地水准面MSE最小准则采用测试法确定最优的正则化参数,例如前文讨论不同正则化矩阵时的比较分析方法,但该准则无法直接应用于实测数据处理,因为无法获知真实重力场模型,因此必须找到确定最优正则化参数的方法。目前大地测量领域讨论较多正则化参数选择方法的是L曲线法和 GCV方法,本文分别对这两种方法在利用梯度观测数据恢复重力场中应用的可行性和有效性进行数值分析讨论。采用的模拟数据同前,选择 Kaula正则化矩阵,分别采用L曲线法和GCV方法确定最优的正则化参数。L曲线如图5所示,图中五角星表示L曲线的最大曲率点,对应正则化参数为1011.3,从图中可看出L曲线的顶点(最大曲率点)并不明显。GCV函数曲线如图6所示,图中五角星表示函数最小值对应的点,正则化参数α为1010.3。

图5 Kaula正则化矩阵对应的L曲线Fig.5 L-curve according to Kaula matrix

图6 Kaula正则化矩阵对应的GCV函数曲线Fig.6 The curve of GCV function according to Kaula matrix

基于这两种方法确定的正则化参数求解模型的大地水准面误差的统计结果如表1所示。表中RMS的下标±90°和±80°分别表示在计算范围为纬度 -90°～90°和-80°～80°,经度取0°～360°,表中还给出了由大地水准面误差的RMS最小为准则确定的最优正则化参数对应模型的大地水准面误差的统计结果。

表1 GCV方法和L曲线法求解模型的大地水准面误差RMS统计Tab.1 The statistics of the geoid error RMS of model estimated with GCVand L-curve

从表1中可看出,L曲线方法确定的正则化参数α相对最优值1010.7偏大,GCV方法确定的偏小,两者均不是理论上的最优。虽然 GCV方法确定模型的精度与最优值更为接近,但两种方法获得解的精度与最优解基本在同一量级,这说明了L曲线法和GCV方法在用Kaula正则化矩阵时确定最优正则化参数的可行性和有效性。从数值的稳定性来讲,GCV方法相对更优,很容易找到GCV函数曲线的最小值,而L曲线法有些时候会由于曲线过于平滑而难于准确找到顶点,缺乏较好的稳定性。

4 结 论

本文在介绍 Tikhonov正则化方法、正则化矩阵的选择以及正则化参数确定方法的基础上,利用GOCE模拟数据基于SA方法分析比较了Tikhonov正则化方法中选择零次、一次和 Kaula正则化矩阵对求解结果的影响,结果说明若能获得最优正则化参数,三类矩阵均能获得较好的结果;基于Kaula正则化矩阵分析比较了确定正则化参数的GCV方法和L曲线法,数值结果说明了两种方法的可行性。需要指出,本文采用的SA方法是一种快速近似方法,本身具有模型误差,这对求解结果会有一定影响,因此文中的分析更多的是定性分析不同正则化矩阵以及正则化参数选择方法的优劣,但总的来说,模拟结果已经说明在利用GOCE观测数据求解高阶全球重力场时采用 Kaula正则化矩阵的有效性,以及采用GCV方法确定最优正则化参数的可行性。

[1] RUMMEL R.Determination of Short-wavelength Components of the Gravity Field from Satellite-to-Satellite Tracking or Satellite Gradiometry:An Attempt to an Identification of Problem Areas[J].Manuscripta Geodaetica,1979,4(2): 107-148.

[2] BOUMAN J,KOOP R.Quality Differences between Tikhonov Regularization and Generalized Biased Estimation in Gradiometric Analysis[J].DEOS Progress Letters,1997,97(1): 42-48.

[3] WANG Zhenjie.Research on the Regularization Solutions of Ill-posed Problems in Geodesy[D].Wuhan:Institute of Geodesy and Geophysics,Chinese Academy of Sciences, 2003.(王振杰.大地测量中不适定问题的正则化解法研究[D].武汉:中国科学院测量与地球物理研究所,2003.)

[4] XU P L.Determination of Surface Gravity Anomalies Using Gradiometric Observables[J].Geophys J Int,1992,110(2): 321-332.

[5] IL K K H.Regularization for High Resolution Gravity Field Recovery by Future Satellite Techniques[C]∥ANGER G,et al.Inverse Problems:Principles and Applications in Geophysics, Technology and Medicine:74.Berlin:Akademie Verlag, 1993:189-214.

[6] KUSCHE J,KLEES R.On the Regularization Problem in Gravity Field Determination from SatelliteGradiometric Data[C]∥ÁDÁM J,SCHWARZ K.Vistas for Geodesy in the New Millennium:IGA 2001 Scientific Assembly.Budapest:Springer,2001:175-180.

[7] KUSCHEJ,KLEES R.Regularization of the Gravity Field Estimation from Satellite Gravity Gradients[J].Journal of Geodesy,2002,76(6-7):359-368.

[8] YANG Y X,SONG L,XU T H.Robust Estimator for Correlated Observations Based on Bifactor Equivalent Weights[J].Journal of Geodesy,2002,76 (6-7): 353-358.

[9] DITMAR P,KUSCHE J,KLEES R.Computation of Spherical Harmonic Coefficients from Gravity Gradiometry Data to be Acquired by the GOCE Satellite:Regularization Issues[J].Journal of Geodesy,2003,77(7-8):465-477.

[10] HANSEN P C.The L-curve and Its Use in the Numerical Treatment of Inverse Problems[C]∥JOHNSTON P. ComputationalInverse Problems in Electrocardiology. Southampton:WIT Press,2001:119-142.

[11] GOLUB G H,HEATH M,WAHBA G.Generalized Cross-validation as a Method for Choosing a Good Ridge Parameter[J].Technometrics,1979,21(2):215-223.

[12] TIKHONOV A N.Regularization of Ill-posed Problem [J].Dokl Akad Nauk SSSR,1963,151(1):49-52.

[13] KAULA W M.Theory ofSatelliteGeodesy[M]. Waltham:Blaisdell Publishing Company,1966:98.

[14] XU Xinyu,WANG Zhengtao,ZOU Xiancai,et al.The Research on theGOCE SatelliteGravity Gradiometry Error Analysis and Simulation[J].Journal of Geodesy and Geodynamics,2010,30(1):1-5.(徐新禹,王正涛,邹贤才,等.GOCE卫星重力梯度测量误差分析及其模拟研究[J].大地测量与地球动力学,2010,30(1):1-5.)

(责任编辑:丛树平)

The Simulation Research on the Tikhonov Regularization Applied in Gravity Field Determination of GOCE Satellite Mission

XU Xinyu1,2,LI Jiancheng1,2,WANG Zhengtao1,2,ZOU Xiancai1,2
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

The Tikhonov regularization is widely applied in the geodesy,the principle of which is discussed in this paper,including the mathematical models of four types of regularization matrices(zero-order,first-order,secondorder and Kaula)and the regularization parameter selection methods:L-curve and GCV.The validation of zeroorder,first-order and Kaula regularization matrices applied in the gravity field determination with GOCE simulated data is analyzed based on the SA method.And the applicability of L-curve and GCV is also discussed using the simulated data.The results show that the accuracies of the optimized solutions(selected by minimizing geoid MSE) with the three types of regularization matrices are at the same level.The key point is the selection of the corresponding regularization parameter.The results also show that GCV and L-curve can be applied in the regularization parameter estimation,and the former method is more stable than the latter one.

Kaula;regularization;GOCE;GCV;L-curve

XU Xinyu(1978—),male,PhD,lecturer, majors in satellite geodesy.

E-mail:xyxu@sgg.whu.edu.cn

1001-1595(2010)05-0465-06

P223

A

国家自然科学基金(40904003,40637034,40704004);国家863计划(2006AA12Z309)

2009-11-02

2010-03-15

徐新禹(1978—),男,博士,讲师,研究方向为卫星大地测量。