于冬冬 王樂洋，2，3

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

2 江西省數(shù)字國土重點(diǎn)實(shí)驗(yàn)室，南昌市廣蘭大道418號，330013

3 流域生態(tài)與地理環(huán)境監(jiān)測國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室，南昌市廣蘭大道418號，330013

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

1 最小二乘譜修正迭代法

對于觀測方程：

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

根據(jù)最小二乘原理，有法方程：

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

式中，Im為m階單位陣。

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

或

2 病態(tài)總體最小二乘譜修正迭代法

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

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

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

根據(jù)總體最小二乘平差準(zhǔn)則：

得法方程為［11］：

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

病態(tài)總體最小二乘譜修正迭代法步驟如下：

4）重復(fù)第2）、3）步，當(dāng)

時(shí)，迭代終止。ε可選一較小正值。

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

3 總體最小二乘譜修正迭代改進(jìn)算法

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

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

4）重復(fù)第2）、3）步，當(dāng)

時(shí)，迭代終止。ε可選一較小正值。

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

4 算例及分析

4.1 算例1

模擬平差模型為病態(tài)的情況，病態(tài)設(shè)計(jì)矩陣和觀測量的真值為［13］：

未知 參 數(shù)X＝［x1x2x3x4x5］T，真 值Xtrue＝［1 1 1 1 1］T。加入隨機(jī)噪聲，對于觀測值L，其觀測噪聲。設(shè)計(jì)矩陣B中的元素與觀測值之間的元素相互獨(dú)立，且其誤差。隨機(jī)誤差由MATLAB隨機(jī)數(shù)發(fā)生器產(chǎn)生。加入隨機(jī)誤差后的設(shè)計(jì)矩陣和觀測向量為：

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

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

4.2 算例2

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

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

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

圖3 解算結(jié)果對比圖Fig.3 Comparison of different methods

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

4.3 算例分析

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

5 結(jié) 語

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

［1］王樂洋，許才軍，魯鐵定.病態(tài)加權(quán)總體最小二乘平差的嶺估計(jì)解法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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］葛旭明，伍吉倉.病態(tài)總體最小二乘問題的廣義正則化［J］.測繪學(xué)報(bào)，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］王新洲，劉丁酉，張前勇，等.譜修正迭代法及其在測量數(shù)據(jù)處理中的應(yīng)用［J］.黑龍江工程學(xué)院學(xué)報(bào)，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ìn)［J］.四川教育學(xué)院學(xué)報(bào)，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īng)儉，廉琦，等.基于法方程病態(tài)的譜修正迭代算法的探討［J］.測繪科學(xué)，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］劉斌，龔健雅，江萬壽，等.基于嶺參數(shù)的譜修正迭代法及其在有理多項(xiàng)式參數(shù)求解中的應(yīng)用［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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］魯鐵定，寧津生.總體最小二乘平差理論及其應(yīng)用［M］.北京：中國科學(xué)技術(shù)出版社，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］王樂洋，于冬冬.病態(tài)總體最小二乘問題的虛擬觀測解法［J］.測繪學(xué)報(bào)，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