趙邵杰 宋迎春 李文娜

(1. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院, 湖南 長(zhǎng)沙 410083; 2. 有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室(中南大學(xué)), 湖南 長(zhǎng)沙 410083)

0 引言

病態(tài)問(wèn)題是大地測(cè)量的數(shù)據(jù)處理中常出現(xiàn)的棘手的問(wèn)題,廣泛存在于控制網(wǎng)平差、精密軌道解算、重力場(chǎng)向下延拓、變形監(jiān)測(cè)、極化干涉合成孔徑雷達(dá)(Polarimetric Interferometric Synthetic Aperture Radar,PolInSAR)植被參數(shù)反演等領(lǐng)域。病態(tài)問(wèn)題法方程的條件數(shù)遠(yuǎn)大于設(shè)計(jì)矩陣的條件數(shù),誤差方程系數(shù)陣或常數(shù)向量的微小擾動(dòng),便會(huì)造成參數(shù)解劇烈變化,難以得到可靠的結(jié)果,主要原因是:系數(shù)陣出現(xiàn)了相對(duì)較小的奇異值,且最大與最小的奇異值相差了幾個(gè)甚至十幾個(gè)量級(jí),這將導(dǎo)致計(jì)算過(guò)程中,未知參數(shù)的方差被較小奇異值過(guò)度放大,顯著降低平差結(jié)果的精度[1]。在病態(tài)問(wèn)題的處理上,目前有許多成熟的方法。如:嶺跡法、GCV法(Generalized Cross-Vali-dation)、Tikhonov正則化法[2]、截?cái)嗥娈愔捣╗3]等。這些方法在一定的條件下可以降低平差模型的病態(tài)性,得到較可靠的平差結(jié)果,不足之處在于無(wú)法利用測(cè)繪工程中的先驗(yàn)信息。

根據(jù)未知參數(shù)的先驗(yàn)信息建立不等式約束,并同誤差方程聯(lián)合平差的方法,稱(chēng)為不等式約束平差法。這種方法補(bǔ)充了平差問(wèn)題的信息量,對(duì)未知參數(shù)形成有效約束,深受?chē)?guó)內(nèi)外專(zhuān)家學(xué)者的廣泛關(guān)注[4-15]。由于不等式約束的存在,它給平差問(wèn)題的解算帶來(lái)了一定的難度,在已有的文獻(xiàn)中,許多學(xué)者重點(diǎn)研究了不等式平差模型的解算方法,然而,他們很少或者沒(méi)有針對(duì)系數(shù)矩陣病態(tài)問(wèn)題展開(kāi)研究。如,在橢球約束算法[15]計(jì)算過(guò)程中,將不等式約束表示為橢球約束的形式,通過(guò)計(jì)算廣義型嶺估計(jì)得到平差結(jié)果,但由于轉(zhuǎn)化過(guò)程會(huì)弱化不等式約束的條件,可能導(dǎo)致平差結(jié)果不在不等式約束范圍內(nèi),從而嚴(yán)重影響參數(shù)解的精度,故這類(lèi)方法受主觀條件的影響較大。在有效約束算法[14]和規(guī)劃類(lèi)算法[11]計(jì)算過(guò)程中就有法矩陣求逆運(yùn)算,或法矩陣的子矩陣求逆運(yùn)算。當(dāng)系數(shù)矩陣病態(tài)時(shí),這些求逆運(yùn)算就會(huì)出現(xiàn)異常,使得解產(chǎn)生嚴(yán)重的偏離。

當(dāng)不等式約束作為一個(gè)先驗(yàn)信息納入平差運(yùn)算時(shí),在病態(tài)的平差模型中可以增加一些虛擬的觀測(cè)信息,這顯然可以對(duì)觀測(cè)信息進(jìn)行補(bǔ)充,從而降低平差模型的病態(tài)性,使得平差結(jié)果更加可靠。這需要在構(gòu)建不等式約束平差算法時(shí),盡量避免對(duì)病態(tài)系數(shù)矩陣進(jìn)行求逆運(yùn)算,否則不等式約束信息還沒(méi)有利用,其病態(tài)問(wèn)題的影響已經(jīng)在參數(shù)解中存在了。因此,避免對(duì)病態(tài)矩陣求逆是求解病態(tài)問(wèn)題的一個(gè)較理想的途徑。

本文針對(duì)未知參數(shù)具有不等式約束的情形,通過(guò)KKT(Karush-Kuhn-Tucker)條件將平差問(wèn)題轉(zhuǎn)化為線性互補(bǔ)問(wèn)題(Linear Complementary Problem，LCP)[16-18],轉(zhuǎn)換后的LCP的系數(shù)矩陣是非對(duì)稱(chēng)的、病態(tài)的,現(xiàn)有的平差算法無(wú)能為力。本文借助Lemke算法給出了針對(duì)非對(duì)稱(chēng)的、病態(tài)的LCP的一個(gè)算法,同時(shí),為利用Lemke算法解決不等式約束平差問(wèn)題提供了一個(gè)新的途徑。新的算法不同于已有的不等式約束算法,在算法構(gòu)建過(guò)程中完全避免了矩陣求逆運(yùn)算,使得病態(tài)性在計(jì)算過(guò)程中對(duì)迭代解的影響大大降低。論文通過(guò)幾個(gè)實(shí)例驗(yàn)證了算法在處理病態(tài)問(wèn)題的有效性,豐富了利用不等式約束解決病態(tài)問(wèn)題的算法。

1 病態(tài)性影響分析

平差模型

L=AX+e權(quán)陣P

(1)

其中,A是m×n維的列滿秩矩陣,rank(A)

min ‖L-AX‖2

(2)

式(2)中,‖·‖2表示2范數(shù)。式(2)稱(chēng)為病態(tài)問(wèn)題,在矩陣A無(wú)病態(tài)條件下,最小二乘解為

(3)

我們對(duì)病態(tài)問(wèn)題的系數(shù)陣A進(jìn)行奇異值分解

(4)

(5)

其中,UA=[u1,…,um]是m階的正交矩陣;GA=[g1,…,gn]是n階的正交矩陣;n×n維的對(duì)角陣Σ=diag(λ1,λ2,…,λn),對(duì)角線元素為A的奇異值,且λ1≥λ2≥…≥λn≥0。

由于系數(shù)陣A病態(tài),我們?cè)O(shè)λN1、λNn為法方程系數(shù)陣N=ATPA的最大奇異值、最小奇異值,則法方程系數(shù)陣的條件數(shù)為

(6)

統(tǒng)計(jì)應(yīng)用經(jīng)驗(yàn)表明,法方程條件數(shù)小于100時(shí)可以認(rèn)為沒(méi)有病態(tài)性;條件數(shù)位于100到1 000之間認(rèn)為存在中等程度的病態(tài)性;條件數(shù)大于1 000認(rèn)為存在嚴(yán)重的病態(tài)性[19]。若條件數(shù)cond(N)很大,即使N和L的擾動(dòng)很小,求解式(3)時(shí)也會(huì)引起X很大的偏差。因此,改善病態(tài)問(wèn)題的一個(gè)可行途徑,就是避免對(duì)法方程中的不穩(wěn)定(條件數(shù)很大)矩陣N求逆,增強(qiáng)解的抗干擾能力。

目前,病態(tài)問(wèn)題的研究主要在兩個(gè)方面,一是對(duì)平差模型進(jìn)行診斷,判斷是否病態(tài);另一方面是研究病態(tài)問(wèn)題的解算方法,例如:正則化方法,截?cái)?修正奇異值方法。附加約束方法等確定性方法和隨機(jī)方法。下面,我們將對(duì)求解病態(tài)問(wèn)題的新算法進(jìn)行研究。

2 不等式約束病態(tài)模型與Lemke算法

2.1 不等式約束病態(tài)模型

最小二乘平差準(zhǔn)則為

minf(X)=(L-AX)TP(L-AX)

(7)

因在目標(biāo)函數(shù)(L-AX)TP(L-AX)=XT(ATPA)X-2LTPAX+LTPL,LTPL是一個(gè)常量,故可以把式(6)轉(zhuǎn)化為一個(gè)二次規(guī)劃問(wèn)題

(8)

式(7)中,c=-(LTPA)T為n維的列向量,N=ATPA為n×n維的對(duì)稱(chēng)正定矩陣,故f(X)是嚴(yán)格凸二次函數(shù),即f(X)的局部最小值為全局最小值。我們將先驗(yàn)信息表示為不等式約束GX≤h,X≥0,其中,G為s×n維的行滿秩矩陣,h為s維的列向量,并聯(lián)合式(1),得到具有不等式約束的病態(tài)模型

L=AX+e權(quán)陣P

(9)

s.tGX≤h,X≥0

(10)

此處,沒(méi)有把式(10)中的兩個(gè)不等式合成一個(gè)不等式,是為了便于式(8)和(9)的轉(zhuǎn)換。由上面的分析可知式(8)和(9)可轉(zhuǎn)換為如下的二次規(guī)劃問(wèn)題(Quadratic Programming):

(11)

s.tGX≤h,X≥0

(12)

二次規(guī)劃問(wèn)題是非線性規(guī)劃中一種特殊的情形,典型的算法有:Lagrange方法、起作用集方法、路徑跟蹤法、Wolfe算法等[22-23]。前三種方法無(wú)法回避對(duì)病態(tài)矩陣求逆,Wolfe算法應(yīng)用廣泛,卻要求N半正定,故不適用于常規(guī)的病態(tài)問(wèn)題。因此,尋找一種不求逆的、大多數(shù)情況適用的解算方法是很有必要的[20],由于Lemke算法既可避免對(duì)病態(tài)矩陣求逆,又可處理N正定和半正定情形,故本文深入研究了這種方法在病態(tài)問(wèn)題中的應(yīng)用。

2.2 線性互補(bǔ)在平差問(wèn)題中的應(yīng)用

本文利用不等式約束求解病態(tài)問(wèn)題,需將已得到的二次規(guī)劃問(wèn)題轉(zhuǎn)化為L(zhǎng)CP,再結(jié)合Lemke算法求解。由于A列滿秩,其法方程系數(shù)陣N是一個(gè)正定矩陣,f(X)為嚴(yán)格凸二次函數(shù),故函數(shù)的局部最小值也是全局最小值,K-T點(diǎn)是式(10)和(11)的最優(yōu)解。根據(jù)Karush-Kuhn-Tucker條件

NX+GTY+c-u=0,v=h-GX≥0

(13)

vTY=0,uTX=0

(14)

u≥0,v≥0,X≥0,Y≥0

(15)

式(10)整理得

(16)

(17)

MZ+q≥0,Z≥0,ZT(MZ+q)=0

(18)

其中,M∈R(s+n)×(s+n)是非對(duì)稱(chēng)矩陣;q∈Rs+n;Z∈Rs+n。

式(9)到式(18)的推導(dǎo),目的是把不等式約束條件下的病態(tài)模型轉(zhuǎn)化為線性互補(bǔ)問(wèn)題求解。式(18)的等價(jià)形式為

w-MZ=q,w≥0,Z≥0,wTZ=0

(19)

由于M是非對(duì)稱(chēng)的、病態(tài)的矩陣,利用一般的線性互補(bǔ)問(wèn)題算法是不能求解式(18)的。對(duì)于式(19)需要研究特殊情形的線性互補(bǔ)問(wèn)題,如非對(duì)稱(chēng)LCP,秩虧LCP等,因此,這些算法直接移植到測(cè)量平差中來(lái)是非常復(fù)雜的,從式(16)和(17)可以看出,w和Z不僅要滿足m+n維未知變量的線性方程組,而且要保證wTZ的內(nèi)積為零。1962年,Lemke基于單純型算法的思想,提出了求解線性互補(bǔ)問(wèn)題式(19)的Lemke算法。這種算法建立在表1的基礎(chǔ)上,可用于處理非對(duì)稱(chēng)的LCP,輸出的結(jié)果完全符合式(19)解的性質(zhì)。由于表1中只含有一個(gè)線性方程組的系數(shù)矩陣,沒(méi)有對(duì)應(yīng)目標(biāo)函數(shù)的系數(shù)向量,故它在結(jié)構(gòu)上比線性單純形表的設(shè)計(jì)更為簡(jiǎn)單。下面,給出了該算法的具體步驟[24]。

表1 線性互補(bǔ)問(wèn)題對(duì)應(yīng)的單純形表

2.3 Lemke算法步驟

第1步將線性互補(bǔ)問(wèn)題的系數(shù)陣M和向量q寫(xiě)入一個(gè)單純形表中。若單純形表的最右列q的分量非負(fù),即q≥0,停止計(jì)算。否則,添加人工變量z0,使之對(duì)應(yīng)的系數(shù)為-1,并將該列寫(xiě)于q的左邊一列。

第2步將q中對(duì)應(yīng)負(fù)數(shù)絕對(duì)值最大的元素所在的行作為主行,選取z0中主行上的元素作為

主元消去的元素,進(jìn)行主元消去,并記錄此時(shí)w出基變量的序號(hào)。

第3步設(shè)出基變量為ws(或zs),則要求下一次必須zs(或ws)將作為進(jìn)基變量。主元消去元素應(yīng)在該列中選取,其對(duì)應(yīng)的主行r應(yīng)滿足:

(20)

其中,dj為當(dāng)前表中即將進(jìn)基變量所在列對(duì)應(yīng)的第j個(gè)分量;pj為當(dāng)前單純形表最右列q對(duì)應(yīng)的第j個(gè)分量。

第4步進(jìn)行主元消去,若此時(shí)人工變量z0出基,轉(zhuǎn)第6步;否則,第5步。

第5步記錄主元消去的出基變量的序號(hào),轉(zhuǎn)第3步。

第6步輸出對(duì)應(yīng)的最優(yōu)的平差結(jié)果Z=(X,Y)T。

算法步驟分析,以下面方程組為例

(21)

算法示意圖分析:

圖1(a)→(b)為第1步至第2步的過(guò)程:建立Lemke算法表,首先根據(jù)q的負(fù)數(shù)絕對(duì)值最大的元素“-10”確定z0列的主行元素w2,接著w2出基z0進(jìn)基,并進(jìn)行主元消去。根據(jù)式(14)定位第3步將主元消去的位置,即(w1,z2)。

圖1(c)→(d)為第3步至第6步的過(guò)程:在對(duì)(w1,z2)進(jìn)行主元消去后,w1出基z2進(jìn)基,我們發(fā)現(xiàn)z0并未出基,故根據(jù)式(20)重新定位第3步主元消去的位置(w3,z1),進(jìn)行主元消去,w3出基z1進(jìn)基。接著重復(fù)第3步定位主元消去的位置,完成進(jìn)基、出基的過(guò)程,直至z0出基,圖1(d)展示了算法結(jié)束時(shí)的狀態(tài)。根據(jù)Lemke算法表,我們得到了LCP的互補(bǔ)可行解

(22)

圖1 Lemke算法示意圖

由于Lemke算法應(yīng)用不等式約束解病態(tài)問(wèn)題過(guò)程,并未涉及到對(duì)病態(tài)矩陣的求逆運(yùn)算,故平差結(jié)果不會(huì)受到病態(tài)矩陣求逆的影響,保證了未知參數(shù)X的穩(wěn)定性;在二次規(guī)劃式(13)中,法方程系數(shù)陣N是正定矩陣,二次規(guī)劃為嚴(yán)格的凸二次規(guī)劃,因此,保證了未知參數(shù)X的唯一性[20]。綜合未知參數(shù)具有穩(wěn)定性和唯一性的特點(diǎn)可知,用Lemke算法求解不等式約束病態(tài)問(wèn)題是一種理論清晰、目的明確的新方式。嶺跡法、GCV法、L-曲線法等是大地測(cè)量數(shù)據(jù)處理中,使用頻率很高的處理病態(tài)問(wèn)題的算法,但它們無(wú)法回避對(duì)病態(tài)矩陣的求逆,嶺估計(jì)方法雖然可以應(yīng)用到病態(tài)矩陣的求逆問(wèn)題,但嶺參數(shù)的確定方法是人為的方法,無(wú)法利用未知參數(shù)X的先驗(yàn)信息,因此,存在著局限性。Lemke已在測(cè)繪數(shù)據(jù)處理中得到應(yīng)用[11],這些應(yīng)用主要還是針對(duì)不等式約束的算法研究,而不是針對(duì)病態(tài)問(wèn)題進(jìn)行研究。下面,結(jié)合兩個(gè)模擬算例進(jìn)行實(shí)驗(yàn)分析,驗(yàn)證Lemke算法的有效性。

3 算例1

Hilbert矩陣是經(jīng)典的病態(tài)矩陣,定義為

(23)

(24)

圖3 L-曲線法示意圖

圖4 嶺跡法示意圖

從平差結(jié)果可以看出,由于系數(shù)矩陣的高度病態(tài)性,使得最小二乘解XLS嚴(yán)重失真,‖XLS-Xreal‖的值高達(dá)3.418 3×1011;方法1和方法2雖然利用了未知參數(shù)X的先驗(yàn)信息GX≤W,在一定程度上提高了解的精度,卻未能改變平差模型AX=L的病態(tài)性,因此,平差的結(jié)果也嚴(yán)重失真。

表2 幾種算法的平差結(jié)果比較

4 算例2

根據(jù)GX≤h,X≥0,將先驗(yàn)信息Δx≤2.3 m,Δy≤4.0 m,Δz≤6.8 m,表示為不等式

表3 系數(shù)矩陣和觀測(cè)向量 單位：m

表4 GPS快速定位的算法結(jié)果與比較 單位：m

(25)

從表4的平差結(jié)果,可以得出以下的結(jié)論:

(1)在本算例中,法方程的條件數(shù)為cond(N)=2.480 9×109,因此,法方程嚴(yán)重病態(tài)。此情況下,按最小二乘的準(zhǔn)則對(duì)法方程進(jìn)行求逆運(yùn)算,必然會(huì)導(dǎo)致未知參數(shù)的解不穩(wěn)定,實(shí)驗(yàn)結(jié)果也驗(yàn)證了最小二乘解嚴(yán)重偏離未知參數(shù)的真值。

(2)GCV法、L-曲線法、嶺跡法均未利用到未知參數(shù)的先驗(yàn)信息,故未得到理想的平差結(jié)果。嶺跡法相對(duì)于其他兩種方法,偏離真值的程度要大得多,一方面與嶺參數(shù)在選擇上受主觀意識(shí)影響較大有關(guān),另一方面也與數(shù)據(jù)結(jié)構(gòu)有關(guān)。對(duì)于不同的數(shù)據(jù),嶺估計(jì)對(duì)病態(tài)問(wèn)題的改善程度也是不同的。

(3)方法1和方法2屬于廣義型嶺估計(jì),具有一定的主觀性。從平差結(jié)果可以看出這兩種方法均未克服法方程的病態(tài)性,方法2在觀測(cè)數(shù)據(jù)和約束條件理想的條件下,也可以在一定程度上對(duì)參數(shù)解改善。

(4)本文提出的不等式約束病態(tài)平差模型的Lemke算法充分利用了先驗(yàn)信息,對(duì)未知參數(shù)形成了有效約束,平差結(jié)果同真值最為接近且效果最好,顯示出了這種算法的優(yōu)越性。

5 結(jié)束語(yǔ)

大地測(cè)量中,未知參數(shù)的先驗(yàn)信息反映了被測(cè)目標(biāo)的幾何或物理特征、內(nèi)在相關(guān)性、測(cè)量精度等,在數(shù)據(jù)處理中,結(jié)合有效的先驗(yàn)信息以等式或不等式約束的形式進(jìn)行平差,可以補(bǔ)充觀測(cè)信息的不足,對(duì)未知參數(shù)形成約束,有效提高參數(shù)估計(jì)的效率。本文針對(duì)病態(tài)問(wèn)題,提供了一種新的平差方法——Lemke算法。本算法沒(méi)有涉及病態(tài)矩陣N=ATPA的求逆,因此,A的病態(tài)性對(duì)Lemke算法影響不大,有效抑制系數(shù)矩陣的復(fù)共線性,降低解的不穩(wěn)定性。其次,Lemke算法在構(gòu)造上簡(jiǎn)明易懂,充分利用測(cè)繪工程中的先驗(yàn)信息解決病態(tài)問(wèn)題,降低了平差模型的復(fù)雜性,豐富了病態(tài)問(wèn)題的算法理論。