預(yù)處理共軛梯度法解病態(tài)問題及在GPS中的應(yīng)用

張同兵,謝 建,朱建軍

(1.中鐵二局機(jī)械筑路工程公司,四川 成都 610031;2.中南大學(xué) 信息物理工程學(xué)院,湖南 長沙 410083)

預(yù)處理共軛梯度法解病態(tài)問題及在GPS中的應(yīng)用

張同兵1,謝 建2,朱建軍2

(1.中鐵二局機(jī)械筑路工程公司,四川 成都 610031;2.中南大學(xué) 信息物理工程學(xué)院,湖南 長沙 410083)

介紹測量數(shù)據(jù)處理中病態(tài)問題的嶺估計(jì)思想及確定嶺參數(shù)的方法。引入預(yù)處理共軛梯度法求解病態(tài)法方程組,說明算法的基本原理和它能有效減弱病態(tài)性的原因。用一個(gè)GPS差分動(dòng)態(tài)定位的實(shí)例,通過與不同嶺估計(jì)方法的對(duì)比分析,驗(yàn)證該方法是一個(gè)有效的數(shù)值迭代算法。

病態(tài)問題;嶺估計(jì);預(yù)處理共軛梯度法;GPS差分動(dòng)態(tài)定位

測量平差中,當(dāng)法方程的系數(shù)陣或者常數(shù)項(xiàng)存在微小的擾動(dòng)時(shí),會(huì)引起解的劇烈變化,這種問題稱為病態(tài)問題。病態(tài)問題廣泛存在于大地測量數(shù)據(jù)處理的各個(gè)領(lǐng)域,如控制網(wǎng)平差、GPS快速定位、攝影測量的自檢校平差和大地測量反演等[1]。多年來,數(shù)學(xué)和大地測量工作者進(jìn)行了大量的研究,并取得了豐富的成果。這些工作主要集中在兩方面:一是病態(tài)系統(tǒng)的診斷,二是病態(tài)問題的求解[2]。病態(tài)方程用最小二乘估計(jì)會(huì)造成解的均方誤差很大而嚴(yán)重失真。目前最常用的是損失解的無偏性而減小均方誤差的有偏估計(jì),包括嶺估計(jì)、主成分估計(jì)、特征根估計(jì)等。近年來又發(fā)展了直接解算的(截?cái)?奇異值分解法,矩陣正交化方法,遺傳算法[3],譜迭代法[4],加權(quán)迭代法[5]等。而最常用的嶺估計(jì)一直是研究的熱點(diǎn),圍繞嶺參數(shù)的選取[6-7],模型中含粗差時(shí)的解算[8],以及嶺估計(jì)與主成分估計(jì)的結(jié)合[9]等方面取得了豐富的成果。這里將引入一種迭代的預(yù)處理共軛梯度算法,通過減弱法方程的病態(tài)性,來獲得較好的解。通過GPS差分動(dòng)態(tài)定位的實(shí)例解算,并與各種嶺估計(jì)方法做比較,驗(yàn)證了該方法能夠獲得均方誤差較小的有偏估計(jì)值。

1 嶺估計(jì)原理及嶺參數(shù)的確定

平差觀測方程為

設(shè)計(jì)陣A呈病態(tài)使最小二乘解

均方誤差過大而不可靠。嶺估計(jì)是從減小均方誤差的角度出發(fā)而提出的一種有偏估計(jì)[10]。它是通過在法矩陣的對(duì)角線加上一個(gè)嶺參數(shù),從而改善法方程系數(shù)的病態(tài)性,使嶺估計(jì)的均方誤差小于最小二乘估計(jì)的均方誤差,定義為

可以證明,嶺估計(jì)^X(k)是最小二乘估計(jì)^X的線性組合,是一種有偏估計(jì)。存在k>0,使得

解此問題的關(guān)鍵是確定嶺參數(shù)。下面分別簡單介紹3類具有代表性并經(jīng)過實(shí)踐證明有效的嶺參數(shù)選取方法。

1.1 嶺跡法

嶺跡法是以^X(k)的分量^X(ki)(i=1,2,…,t,t為參數(shù)個(gè)數(shù))取不同的k值時(shí),以k為橫坐標(biāo),將每個(gè)參數(shù)所對(duì)應(yīng)的^X(k)作為縱坐標(biāo)畫出圖形,并將t條嶺跡畫在一幅圖上,選取它們大體穩(wěn)定時(shí)的那個(gè)k值。這種方法缺少理論依據(jù),選擇比較主觀隨意。

1.2GCV法

運(yùn)用GCV(廣義交叉核實(shí))法確定嶺參數(shù)就是找到GCV函數(shù)的最小值。根據(jù)式(1)和式(2)可以得到GCV函數(shù)為[7]

式中:H(k)=A(ATPA+k E)-1ATP,n為觀測值個(gè)數(shù),tr表示矩陣的跡。根據(jù)式(3)求解的最小值就是GCV法所確定的嶺參數(shù)。GCV法的優(yōu)點(diǎn)是:當(dāng)式(3)的最小值存在時(shí),可以選擇一個(gè)最優(yōu)的嶺參數(shù),它的缺點(diǎn)是:GCV函數(shù)有時(shí)變化過于平緩,這時(shí)定位它的最小值很困難。

1.3 L曲線法

由式(2)可知,‖X(k)‖和‖A X(k)-L‖都是k的函數(shù),選擇不同的k值,以A X(k)-L為橫坐標(biāo),以‖X(k)‖為縱坐標(biāo)畫出類似L形狀的曲線,找到曲率最大的點(diǎn)作為嶺參數(shù)的方法稱為L曲線法[6]。

令η=‖X(k)ξ=‖A X(k)-L‖,取對(duì)數(shù)得到^η=logη,^ξ=logξ。則L曲線是由許多點(diǎn)(^ξ/2,^η/2)擬合而成,用^ξ′、^η′、^ξ″、^η″分別表示^ξ和^η的一、二階導(dǎo)數(shù),則L曲線上任一點(diǎn)的曲率為

對(duì)式(4)求最大值,找到對(duì)應(yīng)的k值就是要求的嶺參數(shù)。

2 預(yù)處理共軛梯度法

嶺估計(jì)是一種直接計(jì)算方法。近年來,一些文獻(xiàn)開始引入數(shù)值迭代方法加權(quán)迭代改善法解病態(tài)線性方程組,并取得了較好的效果。但權(quán)因子選取沒有定論,解決病態(tài)嚴(yán)重問題也無能為力。共軛梯度法是解決系數(shù)矩陣為大型對(duì)稱正定線性方程組的有效方法,當(dāng)矩陣條件數(shù)較大時(shí),盡管收斂速度比較慢,但結(jié)果還是令人滿意。如果引入預(yù)處理技術(shù),降低矩陣條件數(shù),再用共軛梯度法求解,有望得到好的結(jié)果。這里先介紹共軛梯度法和預(yù)處理的基本原理。

2.1 共軛梯度法

共軛梯度法最初由 Hesteness和 Stiefel于1952年為求解線性方程組而提出。解線性方程組A x=b(A對(duì)稱正定)相當(dāng)于求解二次函數(shù)

的極小化問題。

計(jì)算方法為:首先給定初始值x(1),把-▽f(x(1))作為初始下降方向d(1),若‖d(1)‖=0,則停止計(jì)算;否則令

式中:r(1)=b-A x(1)為第1次迭代的殘量。然后沿方向d(1)搜索得到點(diǎn)x(2),構(gòu)造d(2)繼續(xù)搜索。一般地,若已知點(diǎn)x(k)和搜索方向d(k),構(gòu)造

其中步長λk滿足f(x(k+1))

的極小值。對(duì)式(8)求λk的導(dǎo)數(shù),并令它為0,可以計(jì)算得到

代入式(7)可以得到x(k+1),計(jì)算 ‖r(k+1)‖=‖b-A x(k+1)‖,若 ‖r(k+1)‖=0,停止計(jì)算;若‖r(k+1)‖≠0,用‖r(k+1)‖和d(k)構(gòu)造下一個(gè)搜索方向d(k+1),即

重復(fù)使用式(7)、式(9)、式(10)、式(11)可以得到線性方程組的解。

2.2 預(yù)處理方法及與共軛梯度法的結(jié)合

當(dāng)A x=b中A的條件數(shù)很大時(shí),直接用共軛梯度法收斂慢,難以得到好的結(jié)果。如果先對(duì)A的條件數(shù)進(jìn)行改善,化為

是容易求解的方程。如果條件數(shù) cond(-A)小于cond(A),那么用共軛梯度法求解式(12),將得到較好的解y,代入式(13)可以得到原方程的解。這就是預(yù)處理共軛梯度法的基本思想。構(gòu)造預(yù)處理矩陣的方法如下:

其中Q稱為預(yù)處理矩陣,它是對(duì)稱正定的,則存在對(duì)稱正定矩陣C,使得

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

衛(wèi)管畢業(yè)生在就業(yè)觀念上大致可分為兩類。第一，就業(yè)期望過高。經(jīng)濟(jì)社會(huì)逐漸發(fā)展，社會(huì)資源分配出現(xiàn)不均，大城市有更受期待的生活條件，居住環(huán)境、醫(yī)療環(huán)境、子女教育環(huán)境、起點(diǎn)待遇都比一些地級(jí)市、縣城要優(yōu)越一些。因此，許多衛(wèi)管專業(yè)的畢業(yè)生也更傾向于重大城市的、高層次的衛(wèi)生行政機(jī)構(gòu)與醫(yī)藥相關(guān)機(jī)構(gòu)。第二，就業(yè)途徑單一。在求職過程中遭遇挫折后，有些畢業(yè)生便失去繼續(xù)求職的動(dòng)力，轉(zhuǎn)而期待校方推薦的崗位與行業(yè)。而大多數(shù)畢業(yè)生依賴于校園招聘、朋友介紹等方式，導(dǎo)致在就業(yè)問題上比較被動(dòng)。

預(yù)處理矩陣Q的選擇有非常多的方法,包括對(duì)角預(yù)處理法,不完全 Cholesky分解法,矩陣分裂法等。本文取應(yīng)用比較廣泛的對(duì)稱超松弛法(SSOR)[11-12],它是取

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

利用式(17)、(18)構(gòu)造出形如式(12)的方程組,利用共軛梯度法解出間接參數(shù)y,然后代入式(13)得到病態(tài)線性方程組的解。

3 預(yù)處理共軛梯度法在GPS中的應(yīng)用

高精度的GPS定位均采用載波相位作為觀測量,它具有定位速度快、精度高、靈活性強(qiáng)等特點(diǎn),從而在測量和導(dǎo)航中得到廣泛運(yùn)用。而初始相位模糊度的正確確定是利用相位觀測量進(jìn)行精密 GPS相對(duì)定位的關(guān)鍵。傳統(tǒng)的方法是首先進(jìn)行最小二乘估計(jì),得到模糊度的浮動(dòng)解,然后結(jié)合各種解算模糊度的方法來確定模糊度,如FARA方法和LAMBDA方法等。準(zhǔn)確和快速地解算整周模糊度有兩個(gè)前提,一個(gè)是較準(zhǔn)確的模糊度實(shí)數(shù)解,另一個(gè)是較好的模糊度搜索方法。以往的研究側(cè)重于尋找確定模糊度的更好辦法,而本文的研究重點(diǎn)在于減弱法方程的病態(tài)性來改善模糊度實(shí)數(shù)解。

由于雙差GPS可以消去電離層和對(duì)流層延遲誤差、接收機(jī)和衛(wèi)星鐘差,從而得到廣泛應(yīng)用。設(shè)某歷元兩臺(tái)單頻 GPS接收機(jī)可共視k+1顆衛(wèi)星,則可組成k個(gè)線性化后的雙差方程為

式中:i為歷元,Ai為k×3維系數(shù)陣,B為k×k模糊度系數(shù)陣,Xi和N分別為基線向量改正數(shù)和整周模糊度。Li和Δi分別為k維雙差觀測值和誤差向量。假設(shè)有s個(gè)歷元共視了k+1顆衛(wèi)星,則總的觀測方程為

式(20)可以簡化寫為

式中:A為n×m維(n=k×s,m=k+3)系數(shù)陣,Y為m維參數(shù),L和Δ分別為n維雙差觀測值和誤差向量。

在GPS快速定位中,為了在短時(shí)間取得足夠多的觀測量而設(shè)定較高的采樣率,因?yàn)樾l(wèi)星位置在短時(shí)間內(nèi)幾乎不變,相鄰歷元間的坐標(biāo)變化很小。在誤差方程系數(shù)矩陣中,對(duì)應(yīng)于測站點(diǎn)坐標(biāo)的元素是衛(wèi)星的觀測方向相對(duì)于x、y、z3個(gè)坐標(biāo)軸的方向余弦,它們之間近似線性相關(guān),即系數(shù)矩陣A出現(xiàn)復(fù)共線性關(guān)系,這將導(dǎo)致法方程系數(shù)陣呈病態(tài),這時(shí)使用最小二乘估計(jì)求解得到的模糊度的實(shí)數(shù)解與模糊度真值相差較大,在此基礎(chǔ)上使用搜索法很難得到正確的整周模糊度值,或是造成搜索空間過大,嚴(yán)重影響搜索效率。而依靠延長觀測時(shí)間來減弱這種復(fù)共線程度,既不經(jīng)濟(jì)也不及時(shí)。所以有必要改善病態(tài)方程的解。

表1 預(yù)處理共軛梯度法和嶺估計(jì)各種算法的結(jié)果對(duì)比

從表1中可以看出,在GPS快速動(dòng)態(tài)定位模糊度解算中,采用預(yù)處理共軛梯度法能夠得到一組模糊度均方誤差較優(yōu)的解。它比常用的嶺估計(jì)法效果更好,并且不需要確定嶺參數(shù)。

4 結(jié)論與展望

預(yù)處理共軛梯度法是處理大型病態(tài)方程組的有效數(shù)值迭代算法,對(duì)法矩陣進(jìn)行預(yù)處理使條件數(shù)大大減小,從而減弱其病態(tài)性。而GPS動(dòng)態(tài)定位中采樣率高,要同時(shí)用多個(gè)歷元來解算整周模糊度,病態(tài)性非常嚴(yán)重。引入預(yù)處理共軛梯度法有如下優(yōu)點(diǎn):

1)對(duì)于大型病態(tài)問題,采用迭代計(jì)算,不需要求矩陣的逆,存儲(chǔ)量小。

2)不需要確定未知嶺參數(shù),算法簡便,易于編程實(shí)現(xiàn)。

3)預(yù)處理矩陣改善了原法方程的病態(tài)性,在計(jì)算機(jī)上采用長字符數(shù)據(jù)運(yùn)算減小舍入誤差,可以得到更好的解。

當(dāng)然,由于預(yù)處理共軛梯度法是一種數(shù)值迭代方法,不能夠得到參數(shù)估值和觀測之間的明顯表達(dá)關(guān)系,難以評(píng)定解的統(tǒng)計(jì)性質(zhì)。在今后的應(yīng)用中,可以嘗試用不同類型的預(yù)處理方法,選擇不同的預(yù)處理因子,針對(duì)病態(tài)嚴(yán)重程度的不同采取不同方法計(jì)算,有望豐富測量平差的病態(tài)數(shù)據(jù)處理理論,得到更為準(zhǔn)確的解。

[1]王振杰.大地測量中不適定問題的正則化研究[D].武漢:中國科學(xué)院測量與地球物理研究所,2003.

[2]郭建峰.測量平差系統(tǒng)病態(tài)性的診斷與處理[D].鄭州:信息工程大學(xué),2002.

[3]曾群意,歐吉坤.用遺傳算法解病態(tài)方程[J].大地測量與地球動(dòng)力學(xué),2003(3):98-100.

[4]王新洲 ,劉丁酉.最小二乘估計(jì)中法方程的迭代解法[J].湖北民族學(xué)院學(xué)報(bào):自然科學(xué)版,2003(3):1-4.

[5]張春梅,姚晟.用加權(quán)迭代改善法解病態(tài)方程[J].安慶師范學(xué)院學(xué)報(bào):自然科學(xué)版,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ì)原理及其應(yīng)用[J].測繪通報(bào),1994 (1):9-12.

[9]劉雁雨.有偏估計(jì)理論研究及其在 GPS數(shù)據(jù)處理中的應(yīng)用[D].鄭州:信息工程大學(xué),2005.

[10]崔希璋,於宗儔,陶本藻,等.廣義測量平差[M].武漢:武漢大學(xué)出版社,2005.

[11]范嘯濤,季光明.預(yù)優(yōu)矩陣及其構(gòu)造技術(shù)[J].成都理工大學(xué)學(xué)報(bào):自然科學(xué)版,2003(4):432-435.

[12]孫永杰,孫秦.預(yù)處理矩陣及其構(gòu)造方法[J].長春理工大學(xué)學(xué)報(bào),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

國家自然科學(xué)基金資助項(xiàng)目(40574003;40874005)

張同兵(1983-),男,工程師.

[責(zé)任編輯:劉文霞]