張同兵,謝 建,朱建軍

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

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

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

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

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

病態(tài)問題;嶺估計;預處理共軛梯度法;GPS差分動態(tài)定位

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

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

平差觀測方程為

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

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

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

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

1.1 嶺跡法

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

1.2GCV法

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

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

1.3 L曲線法

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

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

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

2 預處理共軛梯度法

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

2.1 共軛梯度法

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

的極小化問題。

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

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

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

的極小值。對式(8)求λk的導數(shù),并令它為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)構(gòu)造下一個搜索方向d(k+1),即

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

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

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

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

其中Q稱為預處理矩陣,它是對稱正定的,則存在對稱正定矩陣C,使得

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

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

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

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

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

3 預處理共軛梯度法在GPS中的應用

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

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

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

式(20)可以簡化寫為

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

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

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

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

4 結(jié)論與展望

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

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

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

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

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

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

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

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

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

[5]張春梅,姚晟.用加權迭代改善法解病態(tài)方程[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數(shù)據(jù)處理中的應用[D].鄭州:信息工程大學,2005.

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

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

[12]孫永杰,孫秦.預處理矩陣及其構(gòu)造方法[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-),男,工程師.

[責任編輯:劉文霞]