非線性不等式約束平差的一種新嶺估計(jì)算法

李文娜 宋迎春 鄧 偉 謝雪梅

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院,長(zhǎng)沙市麓山南路932號(hào)，410083 2 中南林業(yè)科技大學(xué)土木工程學(xué)院,長(zhǎng)沙市韶山南路498號(hào)，410004

處理測(cè)量數(shù)據(jù)病態(tài)問題時(shí)常用嶺估計(jì)[1]方法，但該方法僅從數(shù)學(xué)角度對(duì)平差模型的病態(tài)性進(jìn)行改善，缺乏對(duì)參數(shù)實(shí)際情況的考慮，難以保證顯著提高平差精度。解決病態(tài)問題的另一個(gè)思路是充分利用一些基于前期觀測(cè)和附加(或額外)的先驗(yàn)信息。這些附加的先驗(yàn)信息常表現(xiàn)為參數(shù)的約束條件。

近年來國(guó)內(nèi)學(xué)者對(duì)等式約束[2]、不等式約束[3]、區(qū)間約束[4]的平差理論與算法的研究取得了許多新成果。但這些約束基本上都是線性約束，用來描述先驗(yàn)信息存在明顯不足，對(duì)于解決數(shù)據(jù)處理中的問題也有一定的局限性。楊婷等[5]指出，當(dāng)系數(shù)矩陣具有近似復(fù)共線性(模型病態(tài))時(shí)，參數(shù)的最小二乘估計(jì)的長(zhǎng)度將在向量空間中的某些方向上偏大，此時(shí)一般的線性約束條件不再適合。 最近，學(xué)者們提出一種新的非線性約束形式的“二次約束”，如范數(shù)約束、球約束和橢球約束，其“二次”特性更容易與平差準(zhǔn)則融合，但直接將它們線性化有時(shí)會(huì)使約束條件失去意義。針對(duì)病態(tài)問題的非線性不等式約束平差模型的解算，左廷英等[6]提出一種帶有球約束的迭代算法，該算法只對(duì)參數(shù)進(jìn)行迭代。在實(shí)際應(yīng)用中，通常采用參數(shù)的最小二乘解作為迭代的初始值。但當(dāng)法方程系數(shù)矩陣嚴(yán)重病態(tài)時(shí)，參數(shù)的最小二乘解與參數(shù)真值相距甚遠(yuǎn)，此時(shí)不僅會(huì)延長(zhǎng)迭代的時(shí)間，也會(huì)降低計(jì)算結(jié)果的精度。肖兆兵等[7]給出參數(shù)帶橢球約束平差的具體算法，該算法獲取高精度解的前提是需要選擇一個(gè)合適的拉格朗日乘子迭代初始值。對(duì)于不同的觀測(cè)數(shù)據(jù)，應(yīng)該有不同的初始值，由于缺乏對(duì)平差模型的深入了解，該算法通常選取0作為迭代初始值，導(dǎo)致可能會(huì)出現(xiàn)迭代不收斂的現(xiàn)象。綜合以上問題，本文基于嶺參數(shù)與約束條件有關(guān)這一特性，提出一種新的嶺估計(jì)算法求解非線性不等式約束平差模型，不僅可以保證迭代收斂，而且還能提高計(jì)算效率。

1 非線性不等式約束平差模型及其算法

一般的平差模型可以表示為：

L=AX+e

(1)

式中，L為m×1的線性化后的常數(shù)向量，A為m×n的觀測(cè)值系數(shù)矩陣，設(shè)m≥n，rank(A)=n，X為n×1的未知參數(shù)向量，e為隨機(jī)誤差向量，期望為0(e～N(0,σ2I),Σ=σ2I為協(xié)因數(shù)矩陣)。

如果參數(shù)帶有不確定性，對(duì)X加以先驗(yàn)約束，得到：

式中，‖·‖ 表示二范數(shù)，c≥0且已知。根據(jù)最小二乘準(zhǔn)則，將式(2)轉(zhuǎn)化為非線性不等式約束最小二乘問題。當(dāng)X沒有約束時(shí)，直接利用最小二乘法可求得參數(shù)解為XLS。若‖XLS‖2≤c，則該最小二乘解就是平差問題(2)的解；當(dāng)系數(shù)矩陣病態(tài)時(shí)，往往有‖XLS‖2>c，式(2)中的非線性不等式約束問題可轉(zhuǎn)化為非線性等式約束問題：

min(L-AX)T(L-AX)

s.t.‖X‖2=c

(3)

式(3)也可通過凸優(yōu)化理論得到，因?yàn)?L-AX)T(L-AX)是一個(gè)凸函數(shù)，其最優(yōu)值將在邊界上取得。使用拉格朗日乘數(shù)法求解式(3)：

X(λ)=(ATA+λI)-1ATL

(4)

式(4)為通常的嶺估計(jì)，其中，λ的取值與觀測(cè)值相關(guān)，即與參數(shù)的先驗(yàn)信息有關(guān)。

為研究λ的確定方法，采用Householder變換將矩陣A進(jìn)行雙對(duì)角化[8]：

式中，P為m階的正交矩陣，Q為n階的正交矩陣，B為n×n的上雙對(duì)角矩陣，且B和A的奇異值相同。因此，如果B的奇異值分解為:

B=CSDT

(6)

那么，

X1(λ)=(BTB+λI)-1BTL1

(7)

式(7)等價(jià)于求解：

為減少計(jì)算量，利用Givens變換簡(jiǎn)化式(8)。設(shè)W為Givens矩陣的乘積，即為2n階的正交矩陣，則有：

式中，Bλ為n階的上雙對(duì)角矩陣。因此，X1(λ)可由式(10)計(jì)算：

BλX1(λ)=Z1

(10)

由X1(λ)=QTX(λ)得‖X(λ)‖2=‖X1(λ)‖2。

設(shè)

ω(λ)=‖X(λ)‖2=‖X1(λ)‖2=

將式(6)代入式(11)：

式中，σi為B的奇異值，=CTL1。

顯然， 對(duì)于λ>0，ω(λ)是單調(diào)遞減函數(shù)。由于ω(0)=‖XLS‖2>c，因此存在唯一的λ*>0，使ω(λ*)=c。對(duì)于非線性方程ω(λ)=c，很難求得其精確解。通常情況下，求出的λ只需要使‖X1(λ)‖2/c與1足夠接近即可。

根據(jù)Halley迭代公式[9]可以得到λ的迭代公式為：

具體迭代步驟為：1)利用先驗(yàn)信息確定約束值c；2)設(shè)定λ的初始值λ0和一個(gè)非常小的限值ε；3) 將λn代入式(9)求出Bλn，根據(jù)式(10)得到X1(λn)，從而求出ω(λn)，計(jì)算|ω(λn)-c|<ε是否成立，若成立取出X1(λn)，則X(λn)=QX1(λn)作為最終解，迭代終止，反之則進(jìn)行步驟4)；4)利用X1(λn)求出ω′(λn)、ω″(λn)，然后再利用式(13)求出λn+1，回到步驟3)。

在上述算法中，迭代初始值λ0的確定非常關(guān)鍵，只有選取適當(dāng)?shù)某跏贾?，才能保證其迭代的收斂性，同時(shí)提高計(jì)算效率。根據(jù)上文可知，該算法中嶺參數(shù)λ與參數(shù)的先驗(yàn)信息有關(guān)，基于此特性，下面推導(dǎo)嶺參數(shù)的估計(jì)式，將其作為迭代的初始值，并闡述其合理性。

2 參數(shù)的確定方法

假設(shè)矩陣B的奇異值從大到小排列為：

σ1≥…≥σn-1≥σn

對(duì)于測(cè)量平差數(shù)據(jù)處理中的一些病態(tài)問題，其系數(shù)矩陣中往往會(huì)存在一個(gè)接近于0的特征值，也就是說σn?σn-1，當(dāng)n≠0時(shí)，有：

即

由此可得到λ的估計(jì)值：

λ*

(14)

ω()的值不會(huì)遠(yuǎn)大于c，否則令作為迭代初始值相比于0作為迭代初始值的優(yōu)勢(shì)不明顯。

定義一個(gè)新函數(shù)：

φ(σ)=(σ+/σ)2

由此函數(shù)的性質(zhì)可知，當(dāng)σ>0時(shí)，φ(σ)有唯一的最小值，假設(shè)該最小值在σ=σ*處取得，則：

φ′(σ*)=2(σ*+/σ*)(1-/(σ*)2)=0

即

圖1 θ>1時(shí)φ(σ)的圖像Fig.1 Image of φ(σ) when θ>1

圖2 θ≤1時(shí)φ(σ)的圖像Fig.2 Image of φ(σ) when θ≤1

1)假設(shè)φ(σi)≥φ(σn),i

(σi-σn)(σi-σnθ)≥0

(15)

如果θ≤1，則式(15)成立；如果θ>1，且σi>σnθ，說明B的奇異值分布范圍不在區(qū)間Y內(nèi)，即σi?Y，則式(15)仍成立。在這2種情況下有：

令η=‖‖2/，可以將式(16)簡(jiǎn)化為：

即

‖X1()‖2/c≤η

2)假設(shè)θ>1，且B的奇異值有一些分布在區(qū)間Y內(nèi)，則存在：

因此

綜上所述，在θ≤1(cLS≤4c)的情況下，‖X1()‖2/c≤η。在θ>1(cLS>4c)的情況下，當(dāng)σi?Y(i

3 數(shù)值實(shí)驗(yàn)

3.1 數(shù)值模擬實(shí)驗(yàn)

為驗(yàn)證本文新嶺估計(jì)算法在處理病態(tài)問題中的有效性，采用文獻(xiàn)[10]給出的第一類Fredholm積分方程的具體函數(shù)形式模擬數(shù)據(jù)進(jìn)行計(jì)算，并將所得結(jié)果與利用L-曲線法確定嶺參數(shù)的嶺估計(jì)方法計(jì)算的結(jié)果進(jìn)行比較，見圖3。

圖3 本文算法與嶺估計(jì)法的比較Fig.3 Comparison of our algorithm and ridge estimation algorithm

從圖3可以看出，用本文方法計(jì)算出的函數(shù)曲線與真值的函數(shù)曲線吻合程度更高，本文方法計(jì)算的結(jié)果與真值的誤差平方和為0.080 1，而利用嶺估計(jì)法計(jì)算的結(jié)果與真值的誤差平方和達(dá)到0.713 6。雖然兩者都能夠?qū)ο禂?shù)矩陣的病態(tài)性進(jìn)行改善，但本文算法考慮了參數(shù)的先驗(yàn)信息，所以精度更高。

3.2 病態(tài)測(cè)邊網(wǎng)實(shí)驗(yàn)

病態(tài)測(cè)邊網(wǎng)見圖4，P10、P11為未知點(diǎn)，假設(shè)其坐標(biāo)真值分別為(0 m，0 m，0 m)和(7 m， 10 m，-5 m)，其坐標(biāo)近似值分別取(0.01 m，-0.01 m，0.02 m)和(7.01 m，9.99 m，-5.01 m)。P1，P2，…，P9為已知點(diǎn)，其坐標(biāo)及到2個(gè)未知點(diǎn)的觀測(cè)距離見表1。P10與P11之間的觀測(cè)距離為d10,11=13.107 9 m，各觀測(cè)量精度相等，中誤差為0.001 m。

圖4 空間測(cè)邊網(wǎng)[11]Fig.4 Spatial geodesic netwo

表1 已知點(diǎn)的坐標(biāo)和邊長(zhǎng)觀測(cè)值[11]Tab.1 Coordinates of known points and side length observations

表2 各種方法對(duì)病態(tài)問題解算結(jié)果的比較Tab.2 Comparison of solution results among various methods for ill-posed problem

表3 初始值對(duì)迭代次數(shù)的影響Tab.3 Influence of initial values on the number of iterations

從表3可以看出，將嶺參數(shù)的估計(jì)式作為迭代的初始值是合理的，不僅能保證迭代收斂，同時(shí)還能提高計(jì)算效率。

4 結(jié) 語

大地測(cè)量中存在大量關(guān)于參數(shù)的先驗(yàn)信息，利用這些先驗(yàn)信息對(duì)參數(shù)加以約束，可以有效改善平差模型的病態(tài)性。本文針對(duì)參數(shù)帶有非線性不等式約束的最小二乘問題，建立相應(yīng)的平差模型，給出求解該平差模型的一種迭代算法。該算法在解算效率和解算精度上有較明顯的優(yōu)勢(shì)，但由于誤差的隨機(jī)性和參數(shù)真值的未知性，導(dǎo)致參數(shù)的約束范圍通常需要采用前期數(shù)據(jù)處理的結(jié)果來確定，所以此時(shí)解算結(jié)果的質(zhì)量會(huì)受到前期數(shù)據(jù)處理提供的先驗(yàn)信息的影響。通常在實(shí)際應(yīng)用中遇到的問題比函數(shù)模型所展示的問題要復(fù)雜得多，如何充分利用先驗(yàn)信息獲取可靠的參數(shù)信息，進(jìn)而提高算法的準(zhǔn)確性，是接下來研究工作的重點(diǎn)。