病態總體最小二乘問題的譜修正迭代法

于冬冬 王樂洋，2，3

1 東華理工大學測繪工程學院，南昌市廣蘭大道418號，330013

2 江西省數字國土重點實驗室，南昌市廣蘭大道418號，330013

3 流域生態與地理環境監測國家測繪地理信息局重點實驗室，南昌市廣蘭大道418號，330013

在測量數據處理中，病態問題會造成參數估值與真值相差較大且參數解極不穩定。當前病態總體最小二乘問題的處理方法主要有嶺估計法［1］、Tikhonov 正 則 化 方 法［2］、廣 義 正 則 化 方法［3］等。然而，嶺估計法及正則化方法破壞了方程原有的等量關系且解是有偏估計［4－5］。針對此不足，王新洲等［5］提出病態問題的最小二乘譜修正迭代法，不僅保持了方程的等量關系，且解是無偏估計。文獻［6－7］分別對譜修正迭代法和改進的譜修正迭代法的收斂性進行研究。目前譜修正迭代法及其改進算法已應用在GPS快速定位研究［8］、Bursa模型空間直角坐標轉換［9］、有理多項式參數求解［10］等問題中。本文在最小二乘譜修正迭代法基礎上，提出病態總體最小二乘問題的譜修正迭代法，保留了最小二乘譜修正迭代法解無偏和保持方程的等量關系等優點。

1 最小二乘譜修正迭代法

對于觀測方程：

式中，L為n×1 觀測向量；e為觀測向量的隨機誤差；B為n×m系數矩陣；X為m×1未知參數向量。

根據最小二乘原理，有法方程：

式中，N＝BTB，U＝BTL。將式（2）兩邊同時加上，整理得［5］：

式中，Im為m階單位陣。

LS譜修正迭代法公式為［5］：

或

2 病態總體最小二乘譜修正迭代法

總體最小二乘平差問題的函數模型是：

式中，L∈Rn×1為觀測向量；e為觀測向量的隨機誤差；B∈Rn×m為列滿秩系數矩陣；EB為系數矩陣的誤差；X為m×1 待估參數向量。假設系數矩陣和觀測向量相互獨立且為等精度觀測，則其隨機模型為：

式中，eB是將矩陣EB按列拉直得到的列向量；為觀測值和系數矩陣元素的驗前單位權方差；Im和In分別為m階和n階單位陣；vec（·）表示矩陣的拉直運算；?表示Kronecker積。

根據總體最小二乘平差準則：

得法方程為［11］：

將式（8）兩邊同時加上，整理得：

病態總體最小二乘譜修正迭代法步驟如下：

4）重復第2）、3）步，當

時，迭代終止。ε可選一較小正值。

（BTB＋Im）是秩為m的滿秩矩陣［5］。由于式（8）的解與奇異值分解（SVD）法等價［11］，此時的總體最小二乘解為無偏［12］，譜修正迭代法只是在式（8）的兩端同時加上，因此譜修正迭代法得到的解也是無偏的。與嶺估計及正則化法相比，病態總體最小二乘譜修正迭代法沒有改變方程的等量關系，也不存在正則化因子的確定問題。

3 總體最小二乘譜修正迭代改進算法

譜修正迭代算法雖然對法方程系數矩陣的病態性進行了一定的修正，但在有些情況下，修正后的法方程系數矩陣（BTB＋Im）的條件數仍然很大，病態性依然存在。為此，本文在最小二乘譜修正迭代的改進算法［6－7］基礎上提出總體最小二乘譜修正迭代改進算法。

總體最小二乘譜修正迭代改進算法的計算步驟如下：

4）重復第2）、3）步，當

時，迭代終止。ε可選一較小正值。

選擇譜修正參數α的基本思想為：令α的取值區間為0至一個相對較大的數值，選擇一個相對較小的步長，讓α以該步長取區間內的所有值，對每一個確定的α值都可以得到相對應的修正法方程系數矩陣（BTB＋αIm）的條件數。以譜修正參數α為橫坐標，（BTB＋αIm）的條件數為縱坐標，畫出修正后的法矩陣的條件數隨加入的譜修正參數α的變化曲線圖，根據對法矩陣條件數的改善程度來確定α的最終取值。

4 算例及分析

4.1 算例1

模擬平差模型為病態的情況，病態設計矩陣和觀測量的真值為［13］：

未知 參 數X＝［x1x2x3x4x5］T，真 值Xtrue＝［1 1 1 1 1］T。加入隨機噪聲，對于觀測值L，其觀測噪聲。設計矩陣B中的元素與觀測值之間的元素相互獨立，且其誤差。隨機誤差由MATLAB隨機數發生器產生。加入隨機誤差后的設計矩陣和觀測向量為：

法方程系數陣N＝BTB的條件數為2.083 7×104，病態性嚴重，利用譜修正參數進行改正后，法方程系數矩陣（BTB＋αIm）的條件數隨譜修正參數的變化如圖1所示。在不同迭代初值和譜修正參數α下分別利用LS算法、LS譜修正迭代算法、TLS 算法及TLS 譜修正迭代算法進行求解，結果見表1。

圖1 法方程系數矩陣條件數隨譜修正參數α的變化Fig.1 The picture of the condition number of the corrected normal equation coefficient matrix byα

4.2 算例2

采用文獻［14］“Regularization Tools”中病態模擬數據“ilaplace”。其中，系數矩陣條件數的數量級為1033，病態性非常嚴重。分別對系數矩陣和觀測向量加入σ0＝0.001的隨機誤差。加入譜修正參數后的修正法矩陣（BTB＋αIm）的條件數與譜修正參數α的關系如圖2所示，最小二乘譜修正迭代法和總體最小二乘譜修正迭代法分別在不同迭代初值和譜修正參數下得到的結果如圖3（a）、（b）、（c）所示。圖中，IMCCV－LS（0.4）表示最小二乘譜修正迭代法在初值為0.4＊ones（60，1）下得到的結果、IMCCV－TLS（0.4）表示總體最小二乘譜修正迭代法在初值為0.4＊ones（60，1）下得到的結果，圖中其他注釋同上。各方法所得參數估值與真值的偏差范數及迭代次數見表2。

圖2 修正法方程系數矩陣條件數隨譜修正參數α的變化圖Fig.2 The picture of the condition number of thecorrected normal equation coefficient matrix byα

表1 不同方法得到的結果Tab.1 The results from different methods

圖3 解算結果對比圖Fig.3 Comparison of different methods

表2 偏差范數與迭代次數Tab.2 Bias norm and iteration steps

4.3 算例分析

2）TLS譜修正迭代法盡管能得到較好的效果，但受迭代初值的影響很大，不同的迭代初值得到的結果差異較大，在初值選擇不合理的情況下，TLS譜修正迭代方法甚至不收斂。

5 結 語

本文在最小二乘譜修正迭代算法的基礎上對其進行拓展研究，提出總體最小二乘的譜修正迭代法及其改進算法。該方法保持了最小二乘譜修正迭代方法中方程的等量關系，解是無偏估計。通過算例分析了TLS譜修正迭代法的優勢和缺點。文中的迭代初值為多次試驗下的經驗值，譜修正參數α選取的主觀性較強，如何更加合理有效地選取迭代初值及譜修正參數α還需進一步研究。

［1］王樂洋，許才軍，魯鐵定.病態加權總體最小二乘平差的嶺估計解法［J］.武漢大學學報：信息科學版，2010，35（11）：1 346－1 350（Wang Leyang，Xu Caijun，Lu Tieding.Ridge Estimation Method in Ill－posed Weighted Total Least Squares Adjustment［J］.Geomatics and Information Science of Wuhan University，2010，35（11）：1 346－1 350）

［2］Golub G H，Hansen P C，O’Leary D P.Tikhonov Regularization and Total Least Squares［J］.SIAM Journal on Matrix Analysis and Applications，1999，21（1）：185－194

［3］葛旭明，伍吉倉.病態總體最小二乘問題的廣義正則化［J］.測繪學報，2012，41（3）：372－377（Ge Xuming，Wu Jicang.Generalized Regularization to Ill－posed Total Least Squares Problem［J］.Acta Geodaetica et Cartographica Sinica，2012，41（3）：372－377）

［4］Shen Y Z，Xu P L，Li B F.Bias－Corrected Regularized Solution to Inverse Ill－posed Models［J］.Journal of Geodesy，2012，86（8）：597－608

［5］王新洲，劉丁酉，張前勇，等.譜修正迭代法及其在測量數據處理中的應用［J］.黑龍江工程學院學報，2001，15（2）：3－6（Wang Xinzhou，Liu Dingyou，Zhang Qianyong，et al.The Iteration by Correcting Characteristic Value and Its Application in Surveying Data Processing［J］.Journal of Heilongjiang Institute of Technology.2001，15（2）：3－6）

［6］潘朝毅.譜修正迭代法的收斂分析及其改進［J］.四川教育學院學報，2009，25（5）：112－113（Pan Chaoyi.An Improvement on the Iteration Method by Correcting Characteristic Value［J］.Journal of Sichuan College of Education，2009，25（5）：112－113）

［7］黃德倫，史經儉，廉琦，等.基于法方程病態的譜修正迭代算法的探討［J］.測繪科學，2013，38（2）：35－37（Huang Delun，Shi Jingjian，Lian Qi，et al.Iterative Algorithm Based on the Morbid Equation［J］.Science of Surveying and Mapping，2013，38（2）：35－37）

［8］劉立龍，唐詩華，文鴻雁.基于譜迭代修正法GPS快速定位研究［J］.工程勘察，2007（12）：54－56（Liu Lilong，Tang Shihua，Wen Hongyan.The Rapid Positioning of GPS Based on Iteration Method by Correcting Characteristic Value［J］.Geotechnical Investigation ＆Surveying，2007（12）：54－56）

［9］侯東陽，張書畢，萬亞豪，等.基于譜修正迭代的Bursa模型［J］.海洋測繪，2011，31（5）：15－17（Hou Dongyang，Zhang Shubi，Wan Yahao，et al.Bursa Model Based on the Spectrum Correction Iteration Method［J］.Hydrographic Surveying and Charting，2011，31（5）：15－17）

［10］劉斌，龔健雅，江萬壽，等.基于嶺參數的譜修正迭代法及其在有理多項式參數求解中的應用［J］.武漢大學學報：信息科學版，2012，37（4）：399－402（Liu Bin，Gong Jianya，Jiang Wanshou，et al.Improvement of the Iteration by Correcting Characteristic Value Based on Ridge Estimation and Its Application in RPC Calculating［J］.Geomatics and Information Science of Wuhan University，2012，37（4）：399－402）

［11］魯鐵定，寧津生.總體最小二乘平差理論及其應用［M］.北京：中國科學技術出版社，2011（Lu Tieding，Ning Jinsheng.Total Least Squares Adjustment Theory and Its Applications［M］.Beijing：China Science and Techndogy Press，2011）

［12］Van Huffel S，Vandewalle J.The Total Least Squares Problem：Computational Aspects and Analysis［M］.SIAM，Philadelphia，1991

［13］王樂洋，于冬冬.病態總體最小二乘問題的虛擬觀測解法［J］.測繪學報，2014，43（6）：575－581（Wang Leyang，Yu Dongdong.Virtual Observation Method to Ill－posed Total Least Squares Problem［J］.Acta Geodaetica et Cartographica Sinica，2014，43（6）：575－581）

［14］Hansen P C.Regularization Tools：A Matlab Package for Analysis and Solution Discrete Ill－posed Problems［J］.Numerical Algorithms，1994，6（1）：1－35