一種改進的病態問題截斷奇異值方法

謝 建,李明樾

(1.湖南科技大學 煤炭資源清潔利用與礦山環境保護湖南省重點實驗室，湖南 湘潭 411201;2.廣西交通職業技術學院 路橋工程系, 廣西 南寧 530023)

病態問題在GNSS快速定位[1]，回歸分析[2]，衛星重力向下延拓[3-6]，PolInSAR植被參數反演[7]等測繪領域廣泛存在，嚴重影響參數估計的精度和可靠性。病態性產生的原因可能是參數之間存在相關關系，或者觀測值采樣不足，或者觀測結構不合理[8]。處理病態性的方法主要包括嶺估計、正則化方法、截斷奇異值或修正奇異值方法、附加等式或不等式約束、將參數視為隨機變量的Bayes方法或最小二乘配置方法[3]。Bayes方法和最小二乘配置需要假設參數的先驗統計信息，在實際應用中一般不容易預先得到。當正則化矩陣為單位陣時，嶺估計和單參數正則化解具有相同的形式。正則化方法是在殘差加權平方和的最小二乘準則上，加入參數線性組合的最小范數，然后確定合適的因子(正則化參數)來調節數據的擬合量和平滑量[9-10]。線性組合方法一般根據參數的特點，選擇一階或者二階差分操作算子[11]。當設計陣的零空間和正則化矩陣的零空間不相交時，能得到唯一的正則化解。截斷奇異值方法是將系數矩陣進行奇異值分解后，將嚴重擴大參數方差的小奇異值截去，從而得到可靠的解。而修正奇異值方法是對小的奇異值進行修正改善系數的病態性[6]。截斷奇異值方法簡單明了，僅需要考慮小奇異值對方差的影響，從而在大地測量領域得到廣泛關注。截斷奇異值方法的關鍵問題在于截斷參數的選取。常用的有L曲線法、廣義交叉核實(GCV)法、F假設檢驗法和均方誤差極小化方法[4]。林東方提出了顧及截斷偏差影響的截斷奇異值截斷參數確定方法，并與上述常用方法進行比較[7]。目前的研究大部分聚焦于截斷參數的選取本身，而對病態問題的奇異值解的最優化準則缺乏研究。

在正則化方法中，對于某些實際問題，附加合理的正則化矩陣能得到比單位正則化矩陣更好的結果[11]。文中通過對截斷奇異值理論分析發現，對設計矩陣進行奇異值分解后，刪除較小的奇異值及對應的奇異值向量，組成新的秩虧系數矩陣后，截斷奇異值解就是秩虧觀測系統下的最小范數最小二乘解。為了改善截斷奇異值估計的精度，借鑒正則化矩陣的選取方法，將最小范數條件‖x‖2最小替換為參數線性組合的最小范數條件‖Lx‖2最小，并在此條件下給出改進的截斷奇異值分解方法的參數估計表達式，提出用L曲線法確定截斷參數。數值模擬實驗表明，修正的截斷奇異值分解算法能得到比截斷奇異值方法在均方誤差意義下更好的解。

1 截斷奇異值分解及其等價形式

1.1 最小二乘解的奇異值分解表達式

測量數據處理中，一般是將觀測值表示成待定參數的線性或線性化函數，即經典的Gauss-Markov模型：

y=Ax+e,

(1)

(2)

A=UΛVT.

(3)

在最小二乘準則eTPe=min下，Gauss-Markov模型的最小二乘參數估值為：

(4)

將奇異值分解表達式(3)代入最小二乘估計式(4)，可以得到最小二乘解的奇異值表達式及相應的方差矩陣為[4]：

(5)

(6)

1.2 病態問題的截斷奇異值方法

minxTx，

(7)

s.t.φ(x)=(Akx-y)T(Akx-y)=min.

(8)

由于rank(Ak)=k

(9)

(10)

具體證明見文獻[12]和文獻[13]。另外，由奇異值分解定理有[12]：

(11)

將式(11)代入式(10)，可以得到最小范數最小二乘解的奇異值分解形式為：

(12)

從式(12)可以看出，將系數陣A中的小奇異值刪除構造近似矩陣Ak，將Ak代替式(1)中的設計矩陣，得到的最小范數最小二乘解與式(5)比較， 相當于截去對應于最小二乘解特征值較小部分的解分量，從而達到穩定解的目的。與正則化估計和嶺估計類似，截斷奇異值解比最小二乘解具有更小的均方誤差[4]，因此成為解病態問題的一種重要方法。令V2=[vk+1vk+2…vn]，U2=[uk+1uk+2…un]，Λ2=diag[λk+1λk+2…λn]，對式(12)求期望，有：

(13)

可見，TSVD解不再是無偏估計量，其偏差為：

(14)

1.3 截斷參數的確定

為了獲得穩定且質量可靠的參數估值，關鍵在于確定合適的截斷參數。若k值過大，則部分小奇異值未刪除，解不穩定。反之，則丟棄有用的觀測信息從而損失解的精度。這與Tikhonov正則化參數的選取原理相同，過小的正則化參數達不到穩定解的作用，而過大的參數造成過度平滑。L曲線法是確定正則化參數的有效方法，它能合理地平衡病態模型的擬合部分和平滑部分[14-16]。對于正則化解：

(15)

(16)

(17)

2 修正的截斷奇異值分解

由1.2節易知，截斷奇異值方法減弱病態性的方法與Tikhonov正則化方法類似，都是通過附加參數的范數最小條件，取觀測殘差的范數達到最小，從而求得最小范數最小二乘解。Tikhonov正則化解式(15)，實際上是條件極值問題的解[14]。

(18)

式中:ε是一個很小的正數。根據約束極值問題的Kuhn-Tucker條件，式(18)等價于求下列函數的極小值：

(19)

式中:α是Lagrange乘子(即正則化參數)，應滿足非負條件。正則化矩陣L可以取不同的形式。當L=I(表示單位矩陣)時，表示附加參數的范數最小條件。當已知參數間的某些特點時，常用一階或二階差分算子作為正則化矩陣，能得到比單位矩陣作為正則化矩陣更好的估計結果[11]。在重力場確定中，一般用Kaula準則作為正則化矩陣，它是以階方差模型或先驗重力場模型確定的信號方差作為對角元素。類似地，將式(8)的目標函數進行如下修正，從而建立起改進的奇異值算法的目標函數為：

(20)

(21)

與1.2節相同，首先求得式(21)的通解為式(9)，將其代入式(20)有：

(22)

可以證明，當N(Ak)∩N(L)=0，即Ak和L的零空間線性無關時，目標函數式(22)極小條件下的t值為[13]：

t=-(LP)+Lg.

(23)

將式(23)代入式(9)，可以得到x的極小值為：

(24)

(25)

(26)

MTSVD解的偏差為：

(27)

3 算例分析與比較

第一類Fredholm積分方程是典型的病態問題。在大地測量中，由衛星重力數據恢復局部重力場實質上就是對Fredholm積分方程的解算。該方程的一般形式為：

(28)

其中，z(y)表示含有誤差的觀測值;K(x,y)是已知的核函數;f(x)是待定的函數，由觀測值和核函數通過積分方程來估計。本文算例采用文獻[5]提供的核函數和精確已知的函數f(x)來構造病態模型。其中，

(29)

K(xi+1,yj)f(xi+1))Δx.

(30)

式中:j=1,2,…,401。式(30)可以進一步化為：

(31)

首先繪制f(x)的真實曲線f(x)-real，見圖1。可以看到該曲線呈二次拋物線形狀，可以采用二階離散差分算子來表達其先驗信息，它的具體形式為：

(32)

圖1 一次實驗中3種方法的計算結果比較

在大地測量實際應用中，常能預先知道參數間的某些先驗信息，比如重力場中采用Kaula準則等，可以根據問題的實際情況選擇合適的正則化矩陣。分別采用Tikhonov正則化方法(取正則化矩陣為單位矩陣)，截斷奇異值分解方法和改進的截斷奇異值分解方法進行計算。選取其中一次計算的結果展示于圖1。可以看出，在一次實驗中，3種方法都能減弱病態性的影響，得到與真值相近的估計。Tikhonov正則化方法和TSVD方法都是以參數的二次范數最小為平差準則，未顧及參數間的先驗信息，而MTSVD算法由于顧及參數間的先驗信息，其結果與真值更為吻合。

表1 不同方法估值的均方根誤差比較

圖2 重復實驗的誤差估計情況

4 結 論

測量平差模型的病態性主要是由于觀測方程系數矩陣中含有較小的奇異值，這部分較小的奇異值嚴重放大是觀測誤差的影響造成的，常用的截斷奇異值分解方法是通過刪除小奇異值對應的解分量，從而提高解的精度和可靠性。為了充分利用參數中可能存在的可靠先驗信息，進一步提高解的精度，本文提出改進的截斷奇異值分解方法，主要解決以下問題：

1)將系數矩陣進行奇異值分解，刪除小奇異值及對應的奇異向量，得到與系數矩陣近似的秩虧矩陣，證明該秩虧平差模型的最小范數最小二乘解就是常用的截斷奇異值解；

2)借鑒Tikhonov正則化算法的思想，對于能預先已知參數間先驗信息的病態問題，附加合適的正則化矩陣比單位正則化矩陣能獲得更優的解。本文將常規截斷奇異值問題中的最小范數條件擴展為參數間線性函數的最小范數條件，提出改進的截斷奇異值算法的準則，并得到MP廣義逆形式的解；

3)通過對典型病態問題的解算，比較Tikhonov正則化方法、TSVD方法和本文提出的MTSVD算法的性能，從均方根誤差的角度，證明新方法優于前兩種方法，能有效改善常規TSVD解的精度和可靠性。