靶向譜修正迭代法

吳光明,魯鐵定,2,3,鄧小淵,4,吳建江,5

(1.東華理工大學 測繪工程學院,江西 南昌 330013;2.流域生態與地理環境監測國家測繪地理信息局重點實驗室,江西 南昌 330013;3.江西省數字國土重點實驗室,江西 南昌 330013;4.浙江省地理信息中心，浙江 杭州 310000;5.浙江省地球物理技術應用研究所,浙江 杭州 310014)

病態問題在測量數據處理過程例如在控制網平差[1]、GPS快速定位[2-3]、航空重力向下延拓[4-6]等領域時常碰見。處理病態問題常見方法有Tikhonov正則化法[7-8]、嶺估計[9-10]等，然而嶺估計和正則化法改變了方程的等量關系，使得估計結果有偏并且嶺參數和正則化參數確定困難[11-12]。

1 譜修正迭代法

對于經典平差模型

L+Δ=AX.

(1)

式中：L為n×1觀測值向量，A為m×n系數矩陣，X為n×1待定參數，Δ為觀測向量的噪聲，Δ～N(0,σ02I)。其最小二乘估計及估計的協方差為

(2)

(3)

當法矩陣AΤA出現病態(條件數一般大于103)，則法矩陣求逆將會表現得不穩定，導致求解出的參數估值不可靠[12]。

(4)

或者

(5)

(6)

2 靶向譜修正迭代法

譜修正迭代中的(ATA+αI)-1與嶺估計相似，嶺估計式為

(7)

(8)

整理后

(9)

式中：R是靶向矩陣，其構造方法是基于較小特征值對應的特征向量而構成，

(10)

式中：R1是正則化矩陣，Gi是法矩陣ATA的特征向量。小特征值的判定方法可用特征值標準差分量之和占標準差比重達到95%以上[19]，即

(11)

式中：Λi是法矩陣ATA的特征值。將R1作為靶向矩陣并代入式(8)則整理得到本文靶向譜修正迭代式

(12)

譜修正迭代結果是無偏估計，根據文獻[13]證明方法，本文靶向修正迭代結果也是無偏的。只要迭代矩陣α(ATA+αR)-1·R的譜半徑小于1，迭代就能收斂[20]。在對α(ATA+αR)-1·R進行譜分解，得到

α(ATA+αR)-1·R=

(13)

譜修正參數α選擇有多種方法[16-18],沒有統一結論。因此本文根據文獻[18]的方法，選擇一個較長區間、一定步長的參數，討論不同參數條件下，兩種方法迭代結果的比較。

3 算例及分析

3.1 算例1

采用文獻[4]中的模擬病態問題算例，法矩陣條件數是4.1847×105,嚴重病態，其中未知參數的真值為X=[111111111]T。為比較兩種方法的估計結果，取迭代初值為X=[0.80.80.80.80.80.80.80.80.8]T或其他初值均可，不能選擇最小二乘估值作為迭代初值，若選擇最小二乘估值將無法迭代。根據文獻[18]的方法，取譜修正參數α=1，10，20，30，分別用最小二乘估計、譜修正迭代法、靶向譜修正迭代法對這個問題進行解算，結果見表1。

表1 兩種算法的解算結果

圖1 兩種方法結果

從圖1可以看出，兩種方法的差值范數基本相等，結果基本一致；見圖2(a)圖是α=1～50，由于豎軸刻度較大，兩種方法迭代次數也近乎相等；由圖2(b)看出，在α>11時，本文方法迭代次數少于譜修正迭代法。根據這個算例，靶向譜修正迭代法與譜修正迭代法相比，在參數估計和差值范數上沒有較大變化；但在參數逐漸增大時，迭代次數降低，計算效率提高，體現出本文方法的優勢。

圖2 兩種方法迭代次數

3.2 算例2

表2 兩種算法的解算結果

從圖3可以看出，圖3(a)是α=1～17，兩種方法的差值范數基本相等；由圖3(b)看出，在α>17時，本文方法差值范數小于譜修正迭代法。見圖4，1<α<17時，兩種方法迭代次數基本相等；而α>17時，本文方法迭代次數遠小于譜修正迭代法。這個算例表明，參數逐漸增大，兩種方法在參數估計、差值范數、迭代次數上先是沒有較大變化，但參數繼續增大時，本文方法解算的差值范數和迭代次數均降低，進一步驗證靶向譜修正迭代法高效、快速解算。

圖3 兩種方法解算結果

圖4 兩種方法迭代次數

4 結束語

譜修正迭代法是修正法矩陣所有譜且結果是無偏估計，而一般病態問題是法矩陣的幾個譜奇異，因此存在譜多余修正問題。針對該問題，本文提出靶向譜修正迭代法，即將譜修正迭代矩陣中的單位陣變換成靶向矩陣，迭代過程中只修正法矩陣奇異的譜。通過模擬算例對比兩種方法解算的估計結果、偏差范數和迭代次數，得到以下結論：

1)譜修正參數相同條件下，兩種方法迭代矩陣的譜半徑相等，且估計結果均是無偏的；

2)伴隨著參數增大，本文方法的解算的參數估值將優于譜修正迭代，并且偏差范數更低、迭代次數相比也更少。

因此，在探討處理病態問題時，相比譜修正迭代法，本文方法在估計結果更優且快速迭代，計算效率高。本文遞增選取參數，參數選擇未深入分析，因而參數如何選取有待進一步研究。