李文娜 宋迎春 鄧 偉 謝雪梅

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

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

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

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

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

L=AX+e

(1)

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

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

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

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

s.t.‖X‖2=c

(3)

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

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

(4)

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

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

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

B=CSDT

(6)

那么，

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

(7)

式(7)等價于求解：

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

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

BλX1(λ)=Z1

(10)

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

設(shè)

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

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

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

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

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

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

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

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

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

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

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

即

由此可得到λ的估計值：

λ*

(14)

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

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

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

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

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

即

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

圖2 θ≤1時φ(σ)的圖像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)簡化為：

即

‖X1()‖2/c≤η

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

因此

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

3 數(shù)值實驗

3.1 數(shù)值模擬實驗

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

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

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

3.2 病態(tài)測邊網(wǎng)實驗

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

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

表1 已知點的坐標和邊長觀測值[11]Tab.1 Coordinates of known points and side length observations

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

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

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

4 結(jié) 語

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