王 彬 高井祥 劉 超 王 堅(jiān) 周 鋒

(1)中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，徐州 221116 2)江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室，徐州221116)

L曲線法在等價(jià)權(quán)抗差嶺估計(jì)模型中的應(yīng)用*

王 彬1，2)高井祥1)劉 超1，2)王 堅(jiān)1)周 鋒1，2)

(1)中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，徐州 221116 2)江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室，徐州221116)

借鑒相關(guān)抗差估計(jì)的思想，導(dǎo)出一種基于相關(guān)觀測(cè)的等價(jià)權(quán)抗差嶺估計(jì)模型，將L曲線法引入該模型中用于嶺參數(shù)的求解。針對(duì)抗差迭代過(guò)程中等價(jià)權(quán)陣可能不正定而無(wú)法使用經(jīng)典Cholesky分解進(jìn)行單位化的問(wèn)題，根據(jù)雙因子模型的特點(diǎn)，提出一種等價(jià)權(quán)陣的分解算法對(duì)等價(jià)權(quán)陣進(jìn)行單位化，使得L曲線法在迭代過(guò)程中得以成功運(yùn)用。算例結(jié)果表明:該模型能夠有效克服系數(shù)陣病態(tài)和粗差的影響。

抗差嶺估計(jì);L曲線法;嶺參數(shù);雙因子模型;等價(jià)權(quán)陣分解

1 引言

在控制網(wǎng)平差、大地測(cè)量反演、GPS快速定位等數(shù)據(jù)處理中存在大量的病態(tài)問(wèn)題，當(dāng)系數(shù)矩陣病態(tài)時(shí)，經(jīng)典最小二乘得到的解是極不穩(wěn)定的，需對(duì)系數(shù)矩陣的病態(tài)性進(jìn)行處理，一般采用有偏估計(jì)的方法，如嶺估計(jì)，廣義嶺估計(jì)等。另外，由于受到周圍環(huán)境等因素的影響，觀測(cè)值可能會(huì)受到粗差的污染，即使使用有偏估計(jì)，結(jié)果可能也不理想，需要引入抗差有偏估計(jì)方法。

當(dāng)前關(guān)于抗差有偏估計(jì)的研究主要有抗差嶺估計(jì)［1］、抗差泛嶺估計(jì)［2］、正則化抗差估計(jì)［3］以及針對(duì)GPS快速定位的抗差部分嶺估計(jì)［4］等，這些模型還存在著一些問(wèn)題，主要是在抗差迭代過(guò)程中嶺參數(shù)或正則化參數(shù)的選取，大多對(duì)其未作詳細(xì)探討或者仍采用最小二乘嶺估計(jì)所采用的嶺參數(shù)，少有的一些迭代過(guò)程中嶺參數(shù)選取方法精度也不高且較為繁瑣。嶺參數(shù)的選取對(duì)估計(jì)結(jié)果是至關(guān)重要的，由于在抗差迭代過(guò)程中，等價(jià)權(quán)陣可能是不斷變化的，所以在理論上需要不斷求取相對(duì)于每次迭代最合適的嶺參數(shù)。而如何在算法上實(shí)現(xiàn)嶺參數(shù)的較高精度的自動(dòng)選取，是抗差嶺估計(jì)等方法實(shí)際運(yùn)用的關(guān)鍵。

Hansen［5，6］針對(duì)不適定問(wèn)題提出的L曲線法是求解嶺參數(shù)較好的方法，在最小二乘嶺估計(jì)中已得到了成功的運(yùn)用。針對(duì)相關(guān)觀測(cè)，本文將L曲線法引入相關(guān)抗差嶺估計(jì)模型中，根據(jù)雙因子模型的特點(diǎn)，提出一種針對(duì)迭代過(guò)程中等價(jià)權(quán)陣可能不正定而無(wú)法進(jìn)行經(jīng)典喬列斯基(Cholesky)分解的等價(jià)權(quán)陣分解算法，使得等價(jià)權(quán)陣得以順利單位化，從而使L曲線法成功應(yīng)用于抗差迭代過(guò)程中嶺參數(shù)的求解。最后通過(guò)算例，驗(yàn)證了該模型是可靠有效的。

2 等價(jià)權(quán)抗差嶺估計(jì)模型

設(shè)平差的函數(shù)模型為:

式中，ε為相關(guān)觀測(cè)隨機(jī)誤差向量，且E(ε)=0，先驗(yàn)權(quán)陣為P，A為n×t階系數(shù)矩陣，X為t維參數(shù)向量，對(duì)應(yīng)的誤差方程為:

以式(3)作為估計(jì)準(zhǔn)則得到參數(shù)估計(jì)，抵抗粗差的能力很弱。借鑒相關(guān)抗差估計(jì)的思想，用可變函數(shù)ρ(vi，vj)代替vivj，這樣采取的估計(jì)準(zhǔn)則為:

其中ρ為任意選取的函數(shù)，Vi=ai-li，式(5)對(duì)求導(dǎo)并令導(dǎo)數(shù)為零，則:

式中:

由于抗差函數(shù)具有對(duì)稱性［7］，式(6)變?yōu)?

若令Ψj(vi，vj)/vj=wij，則可令=pijwij= pijΨj(vi，vj)/vj

則式(7)可變?yōu)?

將式(2)代入式(8)可得:

采用楊元喜教授的雙因子等價(jià)權(quán)模型［8］，取:

式中wii和wjj為等價(jià)權(quán)因子，采用IGGⅢ方案計(jì)算:

式中，Ngg=A(ATPA+αI)-1AT。

σ0是由中位數(shù)計(jì)算的單位權(quán)中誤差，其計(jì)算公式為［9］:

式中，k0和k1為閾值，實(shí)際計(jì)算中，k0取1.5～2.5，k1取2.5～4.0。

3 L曲線法的基本原理及等價(jià)權(quán)陣分解算法

3.1 L曲線法確定嶺參數(shù)的基本原理

兩邊取對(duì)數(shù)，得:

L曲線上點(diǎn)的曲率的計(jì)算公式為［10］:

該過(guò)程是在權(quán)陣為單位陣的前提下推導(dǎo)的，若權(quán)陣不是單位陣，需先對(duì)其進(jìn)行單位化［11］，一般采用Cholesky分解方法，如對(duì)于式(2)的誤差方程，其單位化過(guò)程為:

對(duì)先驗(yàn)權(quán)陣P進(jìn)行Cholesky分解，使P=CTC，C為上三角陣，這樣式(4)就變?yōu)?

令B=CA，W=Cl，即實(shí)現(xiàn)了權(quán)陣的單位化，再使用L曲線法即可求解嶺參數(shù)。

3.2 等價(jià)權(quán)陣的分解算法

對(duì)于最小二乘嶺估計(jì)，由于先驗(yàn)權(quán)陣P一般為對(duì)稱正定矩陣，可直接使用Cholesky分解進(jìn)行單位化。在抗差嶺估計(jì)的迭代中，先驗(yàn)權(quán)陣P變?yōu)榈葍r(jià)權(quán)陣且是變化的，因此仍采用最小二乘嶺估計(jì)所采用的嶺參數(shù)是不合理的，而需在每次迭代中都使用L曲線法確定合適的嶺參數(shù)。但是如果觀測(cè)值中含有明顯的粗差，權(quán)因子的計(jì)算出現(xiàn)了淘汰段，必然會(huì)出現(xiàn)全零行和列，為對(duì)稱不正定矩陣，無(wú)法使用經(jīng)典的Cholesky分解對(duì)進(jìn)行單位化，從而導(dǎo)致L曲線法不能順利進(jìn)行。

針對(duì)以上問(wèn)題，本文提出一種等價(jià)權(quán)陣的分解算法。如果第i個(gè)觀測(cè)值中含有明顯的粗差，則在權(quán)因子wii的計(jì)算過(guò)程中，必然會(huì)處于淘汰段，使得wii=0，對(duì)于雙因子模型，則的第i行和第i列必然都為零。根據(jù)等價(jià)權(quán)陣的這種特點(diǎn)，按以下步驟對(duì)進(jìn)行分解:

1)將所有數(shù)值為零的權(quán)因子在觀測(cè)值中的序號(hào)存入數(shù)組ZeroW。設(shè)n個(gè)觀測(cè)值中共有m個(gè)觀測(cè)值所對(duì)應(yīng)的權(quán)因子為零，則將其序號(hào)i1，i2，…，im存入ZeroW。

3)對(duì)矩陣G進(jìn)行Cholesky分解得到上三角矩陣M，使得G=MTM，M∈R(n-m)×(n-m)。

4)在矩陣M中插入m行全零行和m列全零列，且使插入后得到的矩陣的全零行和列的序號(hào)對(duì)應(yīng)且與數(shù)組ZeroW中的各元素相同。設(shè)得到的上三角矩陣為C，C∈Rn×n，則C的第i1，i2，…，im行和列為全零行和全零列。

4 算例分析

本算例是對(duì)文獻(xiàn)［10］中的算例5.1進(jìn)行改化而得的，病態(tài)設(shè)計(jì)矩陣為:

其中A為系數(shù)矩陣，L為觀測(cè)值。未知參數(shù)有三個(gè)，X=(X1X2X3)T，其真值為X=(10.0 15.0 6.0)T。

將觀測(cè)值的權(quán)陣改化為對(duì)角線元素為50，非對(duì)角線元素為10的對(duì)稱方陣，設(shè)參數(shù)的近似值為X0=(10.5 14.4 5.7)T。求得法方程的條件數(shù)cond(ATPA)=6.233 5×104，呈病態(tài)。

為分析比較相關(guān)抗差嶺估計(jì)模型的結(jié)果，設(shè)計(jì)以下方案:①不添加任何模擬粗差;②在第二個(gè)觀測(cè)值上加入模擬粗差，且使Δ2=+5.0;③在方案②的基礎(chǔ)上，在第十個(gè)觀測(cè)值上加入模擬粗差，且使Δ10=-3.5;④方案③的基礎(chǔ)上，在第三個(gè)觀測(cè)值上加入模擬粗差，且使Δ3=+4.5。

分別使用最小二乘估計(jì)，最小二乘嶺估計(jì)和本文基于L曲線法的抗差嶺估計(jì)模型對(duì)方案①～④求解參數(shù)，抗差嶺估計(jì)中取k0=2.0，k1=3.0，迭代收斂精度eps=0.000 1。得到的結(jié)果見(jiàn)表1，求得的與真值X差值的二次范數(shù)為為后驗(yàn)單位權(quán)中誤差，m為抗差迭代計(jì)算完成后對(duì)應(yīng)權(quán)因子為零的觀測(cè)值個(gè)數(shù)。

通過(guò)對(duì)表1計(jì)算結(jié)果的對(duì)比分析，可以得到:

1)方案①中觀測(cè)值不含粗差，由于受系數(shù)矩陣病態(tài)的影響，LS估計(jì)與真值偏差較為嚴(yán)重，達(dá)到3.219 7，方案②～④觀測(cè)值中粗差的數(shù)量從一個(gè)增加到三個(gè)，由于受到系數(shù)矩陣病態(tài)和粗差的雙重影響，LS估計(jì)所得的分別為 15.821 8、18.432 1、16.150 8，參數(shù)估值已經(jīng)完全失真，單位權(quán)中誤差也明顯增大。

2)對(duì)于不含粗差的方案①，LS嶺估計(jì)求得的ΔX僅為0.114 3，可見(jiàn)LS嶺估計(jì)有效克服了系數(shù)矩陣的病態(tài)性，而對(duì)于方案②～④，由于LS嶺估計(jì)抵抗粗差的能力很弱，求得的分別為0.963 8、 1.015 8與0.995 7，較LS估計(jì)的結(jié)果有了明顯的改善，但是單位權(quán)中誤差仍然很大，結(jié)果也不是十分理想。

3)對(duì)于方案①，由于觀測(cè)值中不含粗差，所有觀測(cè)值所對(duì)應(yīng)的權(quán)因子均為1，抗差嶺估計(jì)所采用的嶺參數(shù)和各待求參數(shù)的估值與LS嶺估計(jì)的結(jié)果完全相同。對(duì)于含有粗差的方案②～④，采用本文的抗差嶺估計(jì)模型得到的 ΔX分別為0.088 7、0.080 3和0.089 4，各參數(shù)估值與真值十分接近，單位權(quán)中誤差也有明顯改善，觀測(cè)值的精度得到了顯著提高。若在抗差迭代過(guò)程中始終使用LS嶺估計(jì)所采用的嶺參數(shù)則無(wú)法達(dá)到如此的計(jì)算效果，這是由于使用L曲線法求取適合每次迭代的嶺參數(shù)在理論上更加合理。

4)對(duì)于方案②～④，滿足收斂條件時(shí)抗差嶺估計(jì)的最終嶺參數(shù)與LS嶺估計(jì)的嶺參數(shù)相比有了較大的變化，這是在迭代過(guò)程中反復(fù)運(yùn)用等價(jià)權(quán)陣分解算法對(duì)不正定的等價(jià)權(quán)陣進(jìn)行單位化并使用L曲法后所得的最終結(jié)果，圖1是L曲線法求解這三種方案LS嶺估計(jì)的嶺參數(shù)與抗差嶺估計(jì)的最終嶺參數(shù)的示意圖比較。

5 結(jié)語(yǔ)

表1 各方案計(jì)算結(jié)果Tab.1 Computation result under each scheme

圖1 用L曲線法求解兩種估計(jì)方法的嶺參數(shù)示意圖Fig.1 Determining ridge parameters with two estimation methods by use of L-curve method

抗差迭代過(guò)程中嶺參數(shù)的選取是影響抗差嶺估計(jì)結(jié)果的關(guān)鍵因素，本文將L曲線法引入相關(guān)抗差嶺估計(jì)模型中用于嶺參數(shù)的求解，針對(duì)迭代過(guò)程中L曲線法的使用問(wèn)題，提出等價(jià)權(quán)陣分解算法對(duì)等價(jià)權(quán)陣進(jìn)行單位化。算例結(jié)果表明，當(dāng)觀測(cè)值中不含粗差時(shí)，抗差嶺估計(jì)的結(jié)果與LS嶺估計(jì)等價(jià)且明顯優(yōu)于LS估計(jì);當(dāng)觀測(cè)值中含有粗差時(shí)，LS估計(jì)與LS嶺估計(jì)的結(jié)果均失真，而抗差嶺估計(jì)模型有效克服了系數(shù)陣病態(tài)和粗差的雙重影響，得到的各參數(shù)的估值與真值十分接近，觀測(cè)值的精度也得到了顯著提高。

1 隋立芬.抗差嶺估計(jì)原理及其應(yīng)用［J］.測(cè)繪通報(bào)，1994，(1):6-9.(Sui Lifen.The principle of robust ridge estimation and its applications in geodetic adjustments［J］.Bulletin of Surveying and Mapping，1994，(1):6-9)

2 歸慶明，等.抗差泛嶺估計(jì)［J］.測(cè)繪學(xué)報(bào)，1998，(3):211 -216.(Gui Qinming，et al.Robust universal ridge estimation［J］.Acta Geodaetica et Catographica Sinica，1998，(3):211-216)

3 常志巧，等.GPS快速定位中病態(tài)問(wèn)題的正則化抗差解法［J］.大地測(cè)量與地球動(dòng)力學(xué)，2008，(3):83-86.(Chang Zhiqiao，et al.Regularization combined with robust estimation and its application for GPS rapid positioning［J］.Journal of Geodesy and Geodynamics，2008，(3):83-86)

4 歸慶明，等.抗差部分嶺估計(jì)及其在GPS快速定位中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，2006，(2):62-65.(Gui Qinming，et al.Robust partial ordinary ridge estimation and its applications in GPS positioning［J］.Journal of Geodesy and Geodynamics，2006，(2):62-65)

5 Hansen P C.Analysis of discrete ill-posed problems by means of the L-curve［J］.SIAM Review，1992，34(4):561 -580.

6 Hansen P C and O’Leary D P.The use of the L-curve in the regularization of discrete ill-problems［J］.SIAM Review，1993，14(6):1 487-1 503.

7 余學(xué)祥，等.GPS變形監(jiān)測(cè)數(shù)據(jù)處理自動(dòng)化-似單差法的理論與方法［M］.徐州:中國(guó)礦業(yè)大學(xué)出版社，2004.(Yu Xuexiang，et al.Automation of GPS deformation monitoring data processing，the theory and method of SSDM［M］.Xuzhou:China University of Mining and Technology Press， 2004)

8 Yang Y，Song L and Xu T.Robust estimation for correlated observations based on bifactor equivalent weights［J］.J Geodesy.，2002，76:353-358.

9 宋力杰.測(cè)量平差程序設(shè)計(jì)［M］.北京:國(guó)防工業(yè)出版社，2009.(Song Lijie.Programming of surveying adjustment［M］.Beijing:National Defence Industry Press，2009)

10 王振杰.測(cè)量中不適定問(wèn)題的正則化解法［M］.北京:科學(xué)出版社，2006.(Wang Zhenjie.Regularization solutions of ill-posed problem in geodesy［M］.Beijing:Science Press，2006)

11 周江文，等.測(cè)量誤差理論新探［M］.北京:地震出版社，1999.(Zhou Jiangwen，et al.The new research of errors theory［M］.Beijing:Seismological Press，1999)

APPLICATION OF L-CURVE METHOD TO EQUIVALENT WEIGHT ROBUST RIDGE ESTIMAION MODEL

Wang Bin1，2)，Gao Jingxiang1)，Liu Chao1，2)，Wang Jian1)and Zhou Feng1，2)

(1)China University of Mining and Technology，CESSI，Xuzhou 221116 2)Jiangsu Key Laboratory of Resources and Environmental Information Engineering，Xuzhou221116)

An equivalent weight robust ridge estimation model is derived according to the idea of robust estimation for correlated observations，and then the L-curve method is introduced to this model for the solution of ridge parameters.In the iterative process，equivalent weight matrix may be non positive definite so it can not be unitized through classical Cholesky decomposition.According to the characteristic of bifactor model，an algorithm of equivalent weight matrix decomposition is proposed to solve this problem，then the L-curve method can be successfully used in this process.The calculated example manifests that this model can effectively eliminate the influence of coefficient matrix’s morbidity and gross errors.

robust ridge estimation;L-curve method;ridge parameter;bifactor model;equivalent weight matrix decomposition

1671-5942(2012)03-0097-05

2012-01-08

國(guó)家自然科學(xué)基金(41074010);江蘇省高校研究生科研創(chuàng)新計(jì)劃項(xiàng)目(CXLX11_0321)

王彬，男，1988年生，碩士生，研究方向:測(cè)量數(shù)據(jù)處理理論.E-mail:rainkingwang881107@163.com

P207

A