预处理共轭梯度法解病态问题及在GPS中的应用

张同兵,谢 建,朱建军

(1.中铁二局机械筑路工程公司,四川 成都 610031;2.中南大学 信息物理工程学院,湖南 长沙 410083)

预处理共轭梯度法解病态问题及在GPS中的应用

张同兵1,谢 建2,朱建军2

(1.中铁二局机械筑路工程公司,四川 成都 610031;2.中南大学 信息物理工程学院,湖南 长沙 410083)

介绍测量数据处理中病态问题的岭估计思想及确定岭参数的方法。引入预处理共轭梯度法求解病态法方程组,说明算法的基本原理和它能有效减弱病态性的原因。用一个GPS差分动态定位的实例,通过与不同岭估计方法的对比分析,验证该方法是一个有效的数值迭代算法。

病态问题;岭估计;预处理共轭梯度法;GPS差分动态定位

测量平差中,当法方程的系数阵或者常数项存在微小的扰动时,会引起解的剧烈变化,这种问题称为病态问题。病态问题广泛存在于大地测量数据处理的各个领域,如控制网平差、GPS快速定位、摄影测量的自检校平差和大地测量反演等[1]。多年来,数学和大地测量工作者进行了大量的研究,并取得了丰富的成果。这些工作主要集中在两方面:一是病态系统的诊断,二是病态问题的求解[2]。病态方程用最小二乘估计会造成解的均方误差很大而严重失真。目前最常用的是损失解的无偏性而减小均方误差的有偏估计,包括岭估计、主成分估计、特征根估计等。近年来又发展了直接解算的(截断)奇异值分解法,矩阵正交化方法,遗传算法[3],谱迭代法[4],加权迭代法[5]等。而最常用的岭估计一直是研究的热点,围绕岭参数的选取[6-7],模型中含粗差时的解算[8],以及岭估计与主成分估计的结合[9]等方面取得了丰富的成果。这里将引入一种迭代的预处理共轭梯度算法,通过减弱法方程的病态性,来获得较好的解。通过GPS差分动态定位的实例解算,并与各种岭估计方法做比较,验证了该方法能够获得均方误差较小的有偏估计值。

1 岭估计原理及岭参数的确定

平差观测方程为

设计阵A呈病态使最小二乘解

均方误差过大而不可靠。岭估计是从减小均方误差的角度出发而提出的一种有偏估计[10]。它是通过在法矩阵的对角线加上一个岭参数,从而改善法方程系数的病态性,使岭估计的均方误差小于最小二乘估计的均方误差,定义为

可以证明,岭估计^X(k)是最小二乘估计^X的线性组合,是一种有偏估计。存在k>0,使得

解此问题的关键是确定岭参数。下面分别简单介绍3类具有代表性并经过实践证明有效的岭参数选取方法。

1.1 岭迹法

岭迹法是以^X(k)的分量^X(ki)(i=1,2,…,t,t为参数个数)取不同的k值时,以k为横坐标,将每个参数所对应的^X(k)作为纵坐标画出图形,并将t条岭迹画在一幅图上,选取它们大体稳定时的那个k值。这种方法缺少理论依据,选择比较主观随意。

1.2GCV法

运用GCV(广义交叉核实)法确定岭参数就是找到GCV函数的最小值。根据式(1)和式(2)可以得到GCV函数为[7]

式中:H(k)=A(ATPA+k E)-1ATP,n为观测值个数,tr表示矩阵的迹。根据式(3)求解的最小值就是GCV法所确定的岭参数。GCV法的优点是:当式(3)的最小值存在时,可以选择一个最优的岭参数,它的缺点是:GCV函数有时变化过于平缓,这时定位它的最小值很困难。

1.3 L曲线法

由式(2)可知,‖X(k)‖和‖A X(k)-L‖都是k的函数,选择不同的k值,以A X(k)-L为横坐标,以‖X(k)‖为纵坐标画出类似L形状的曲线,找到曲率最大的点作为岭参数的方法称为L曲线法[6]。

令η=‖X(k)ξ=‖A X(k)-L‖,取对数得到^η=logη,^ξ=logξ。则L曲线是由许多点(^ξ/2,^η/2)拟合而成,用^ξ′、^η′、^ξ″、^η″分别表示^ξ和^η的一、二阶导数,则L曲线上任一点的曲率为

对式(4)求最大值,找到对应的k值就是要求的岭参数。

2 预处理共轭梯度法

岭估计是一种直接计算方法。近年来,一些文献开始引入数值迭代方法加权迭代改善法解病态线性方程组,并取得了较好的效果。但权因子选取没有定论,解决病态严重问题也无能为力。共轭梯度法是解决系数矩阵为大型对称正定线性方程组的有效方法,当矩阵条件数较大时,尽管收敛速度比较慢,但结果还是令人满意。如果引入预处理技术,降低矩阵条件数,再用共轭梯度法求解,有望得到好的结果。这里先介绍共轭梯度法和预处理的基本原理。

2.1 共轭梯度法

共轭梯度法最初由 Hesteness和 Stiefel于1952年为求解线性方程组而提出。解线性方程组A x=b(A对称正定)相当于求解二次函数

的极小化问题。

计算方法为:首先给定初始值x(1),把-▽f(x(1))作为初始下降方向d(1),若‖d(1)‖=0,则停止计算;否则令

式中:r(1)=b-A x(1)为第1次迭代的残量。然后沿方向d(1)搜索得到点x(2),构造d(2)继续搜索。一般地,若已知点x(k)和搜索方向d(k),构造

其中步长λk满足f(x(k+1))

的极小值。对式(8)求λk的导数,并令它为0,可以计算得到

代入式(7)可以得到x(k+1),计算 ‖r(k+1)‖=‖b-A x(k+1)‖,若 ‖r(k+1)‖=0,停止计算;若‖r(k+1)‖≠0,用‖r(k+1)‖和d(k)构造下一个搜索方向d(k+1),即

重复使用式(7)、式(9)、式(10)、式(11)可以得到线性方程组的解。

2.2 预处理方法及与共轭梯度法的结合

当A x=b中A的条件数很大时,直接用共轭梯度法收敛慢,难以得到好的结果。如果先对A的条件数进行改善,化为

是容易求解的方程。如果条件数 cond(-A)小于cond(A),那么用共轭梯度法求解式(12),将得到较好的解y,代入式(13)可以得到原方程的解。这就是预处理共轭梯度法的基本思想。构造预处理矩阵的方法如下:

其中Q称为预处理矩阵,它是对称正定的,则存在对称正定矩阵C,使得

用C-1左乘A x=b可得

卫管毕业生在就业观念上大致可分为两类。第一，就业期望过高。经济社会逐渐发展，社会资源分配出现不均，大城市有更受期待的生活条件，居住环境、医疗环境、子女教育环境、起点待遇都比一些地级市、县城要优越一些。因此，许多卫管专业的毕业生也更倾向于重大城市的、高层次的卫生行政机构与医药相关机构。第二，就业途径单一。在求职过程中遭遇挫折后，有些毕业生便失去继续求职的动力，转而期待校方推荐的岗位与行业。而大多数毕业生依赖于校园招聘、朋友介绍等方式，导致在就业问题上比较被动。

预处理矩阵Q的选择有非常多的方法,包括对角预处理法,不完全 Cholesky分解法,矩阵分裂法等。本文取应用比较广泛的对称超松弛法(SSOR)[11-12],它是取

a22,…,ann),0≤ω≤2,一般取ω=1。

利用式(17)、(18)构造出形如式(12)的方程组,利用共轭梯度法解出间接参数y,然后代入式(13)得到病态线性方程组的解。

3 预处理共轭梯度法在GPS中的应用

高精度的GPS定位均采用载波相位作为观测量,它具有定位速度快、精度高、灵活性强等特点,从而在测量和导航中得到广泛运用。而初始相位模糊度的正确确定是利用相位观测量进行精密 GPS相对定位的关键。传统的方法是首先进行最小二乘估计,得到模糊度的浮动解,然后结合各种解算模糊度的方法来确定模糊度,如FARA方法和LAMBDA方法等。准确和快速地解算整周模糊度有两个前提,一个是较准确的模糊度实数解,另一个是较好的模糊度搜索方法。以往的研究侧重于寻找确定模糊度的更好办法,而本文的研究重点在于减弱法方程的病态性来改善模糊度实数解。

由于双差GPS可以消去电离层和对流层延迟误差、接收机和卫星钟差,从而得到广泛应用。设某历元两台单频 GPS接收机可共视k+1颗卫星,则可组成k个线性化后的双差方程为

式中:i为历元,Ai为k×3维系数阵,B为k×k模糊度系数阵,Xi和N分别为基线向量改正数和整周模糊度。Li和Δi分别为k维双差观测值和误差向量。假设有s个历元共视了k+1颗卫星,则总的观测方程为

式(20)可以简化写为

式中:A为n×m维(n=k×s,m=k+3)系数阵,Y为m维参数,L和Δ分别为n维双差观测值和误差向量。

在GPS快速定位中,为了在短时间取得足够多的观测量而设定较高的采样率,因为卫星位置在短时间内几乎不变,相邻历元间的坐标变化很小。在误差方程系数矩阵中,对应于测站点坐标的元素是卫星的观测方向相对于x、y、z3个坐标轴的方向余弦,它们之间近似线性相关,即系数矩阵A出现复共线性关系,这将导致法方程系数阵呈病态,这时使用最小二乘估计求解得到的模糊度的实数解与模糊度真值相差较大,在此基础上使用搜索法很难得到正确的整周模糊度值,或是造成搜索空间过大,严重影响搜索效率。而依靠延长观测时间来减弱这种复共线程度,既不经济也不及时。所以有必要改善病态方程的解。

表1 预处理共轭梯度法和岭估计各种算法的结果对比

从表1中可以看出,在GPS快速动态定位模糊度解算中,采用预处理共轭梯度法能够得到一组模糊度均方误差较优的解。它比常用的岭估计法效果更好,并且不需要确定岭参数。

4 结论与展望

预处理共轭梯度法是处理大型病态方程组的有效数值迭代算法,对法矩阵进行预处理使条件数大大减小,从而减弱其病态性。而GPS动态定位中采样率高,要同时用多个历元来解算整周模糊度,病态性非常严重。引入预处理共轭梯度法有如下优点:

1)对于大型病态问题,采用迭代计算,不需要求矩阵的逆,存储量小。

2)不需要确定未知岭参数,算法简便,易于编程实现。

3)预处理矩阵改善了原法方程的病态性,在计算机上采用长字符数据运算减小舍入误差,可以得到更好的解。

当然,由于预处理共轭梯度法是一种数值迭代方法,不能够得到参数估值和观测之间的明显表达关系,难以评定解的统计性质。在今后的应用中,可以尝试用不同类型的预处理方法,选择不同的预处理因子,针对病态严重程度的不同采取不同方法计算,有望丰富测量平差的病态数据处理理论,得到更为准确的解。

[1]王振杰.大地测量中不适定问题的正则化研究[D].武汉:中国科学院测量与地球物理研究所,2003.

[2]郭建峰.测量平差系统病态性的诊断与处理[D].郑州:信息工程大学,2002.

[3]曾群意,欧吉坤.用遗传算法解病态方程[J].大地测量与地球动力学,2003(3):98-100.

[4]王新洲 ,刘丁酉.最小二乘估计中法方程的迭代解法[J].湖北民族学院学报:自然科学版,2003(3):1-4.

[5]张春梅,姚晟.用加权迭代改善法解病态方程[J].安庆师范学院学报:自然科学版,2004(1):78-79.

[6]HANSEN P.C.Analysis of discrete ill-posed p roblem s by means of the L-curve[J].SIAM Review,1992,34(4): 561-581.

[7]GOLUB G.H.,HEA TH M.,WAHBA G.Generalized cross-validation asamethod fo r choosing a good ridge parameter[J].Techno metrics,1979,21(2):215-223.

[8]隋立芬.抗差岭估计原理及其应用[J].测绘通报,1994 (1):9-12.

[9]刘雁雨.有偏估计理论研究及其在 GPS数据处理中的应用[D].郑州:信息工程大学,2005.

[10]崔希璋,於宗俦,陶本藻,等.广义测量平差[M].武汉:武汉大学出版社,2005.

[11]范啸涛,季光明.预优矩阵及其构造技术[J].成都理工大学学报:自然科学版,2003(4):432-435.

[12]孙永杰,孙秦.预处理矩阵及其构造方法[J].长春理工大学学报,2006(4):128-130.

Solving the ill-posed problem by preconditioned con jugate gradient and itsapplication in GPS

ZHANG Tong-bing1,XIE Jian2,ZHU Jian-jun2

(1.CREGC Machinery Engineering Co.,L td.,Chengdu 610031,China;2.Department of Info-physics Geomatics Engineering, Central South University,Changsha 410083,China)

The ridge estimation and determination of ridge parameter in ill-posed p roblem in surveying data p rocessing are briefly introduced.We bring the p reconditioned conjugate gradient to solve the ill-posed no rmal equation,and the basic p rincip le and its ability to dim inish the illnessof this algo rithm are illustrated first.Then it is app lied to a kinematic DGPSand it p roves to be an effective iterative algo rithm through comparison w ith different ridge estimations.

ill-posed p roblem;ridge estimation;p reconditioned conjugate gradient;kinematic DGPS

P207

A

1006-7949(2010)04-0060-04

2009-09-07

国家自然科学基金资助项目(40574003;40874005)

张同兵(1983-),男,工程师.

[责任编辑:刘文霞]