線性回歸模型的穩(wěn)健總體最小二乘解算

汪奇生 楊德宏 楊騰飛

1 湖南軟件職業(yè)學(xué)院，湘潭市寶馬西路，411100

2 昆明理工大學(xué)國(guó)土資源工程學(xué)院，昆明市文昌路68號(hào)，650093

測(cè)量數(shù)據(jù)處理中，諸如線性擬合、點(diǎn)云平面擬合、多項(xiàng)式擬合等問(wèn)題都可以看成是一個(gè)線性回歸模型來(lái)進(jìn)行平差處理。線性回歸模型有一個(gè)顯著的特點(diǎn)，即其平差模型中的系數(shù)矩陣都含有一列常數(shù)列，長(zhǎng)期以來(lái)這類問(wèn)題的平差都是以最小二乘（LS）為基礎(chǔ)進(jìn)行的，即考慮線性回歸模型中的因變量含有隨機(jī)誤差，而不考慮系數(shù)矩陣中自變量的誤差。相對(duì)于最小二乘，總體最小二乘（TLS）［1］能顧及系數(shù)矩陣含有的隨機(jī)誤差，在理論上更為嚴(yán)密。所以，眾多學(xué)者對(duì)此進(jìn)行研究，并提出一些針對(duì)測(cè)量數(shù)據(jù)處理的總體最小二乘平差算法［2－6］。對(duì)于線性回歸模型的總體最小二乘算法也有相關(guān)文獻(xiàn)進(jìn)行了研究［7－9］。

上述研究都只考慮平差模型中的隨機(jī)誤差，而在總體最小二乘平差時(shí)，這些線性回歸模型的自變量和因變量都可看成是通過(guò)測(cè)量獲取的觀測(cè)值。由于測(cè)量?jī)x器和觀測(cè)等原因，觀測(cè)值中不免含有粗差。一般來(lái)講，對(duì)于線性回歸模型的粗差處理都是采用穩(wěn)健最小二乘方法（RLS）。也有文獻(xiàn)對(duì)同時(shí)考慮自變量誤差的穩(wěn)健總體最小二乘問(wèn)題進(jìn)行了研究。文獻(xiàn)［10－11］從幾何角度計(jì)算每次平差后點(diǎn)到擬合方程的距離，從而刪除那些異常點(diǎn)，然后重新進(jìn)行平差。但這種方法并沒有利用到測(cè)量平差的優(yōu)勢(shì)，且處理較為麻煩。文獻(xiàn)［12］提出的抗差總體最小二乘法只是針對(duì)觀測(cè)向量含有的粗差，而對(duì)系數(shù)矩陣含有的粗差沒有涉及。文獻(xiàn)［13］以非線性平差模型為基礎(chǔ)將總體最小二乘平差模型看成是非線性的，用選權(quán)迭代的方法來(lái)求解三維坐標(biāo)轉(zhuǎn)換的穩(wěn)健參數(shù)。但文獻(xiàn)中沒有將公式統(tǒng)一成矩陣形式，即每一次迭代都需重新構(gòu)造矩陣。本文針對(duì)線性回歸模型將其進(jìn)行等價(jià)轉(zhuǎn)換，根據(jù)文獻(xiàn)［13］的思想，將平差模型看成是非線性的，并通過(guò)選權(quán)迭代，提出線性回歸模型的穩(wěn)健總體最小二乘（LRRTLS）算法，并通過(guò)模擬算例驗(yàn)證了算法的有效性和可行性。

1 穩(wěn)健總體最小二乘算法原理

1.1 總體最小二乘算法

線性回歸數(shù)學(xué)模型為：

由于其平差模型的系數(shù)矩陣含有一列常數(shù)，故在進(jìn)行總體最小二乘平差時(shí)，需要顧及系數(shù)矩陣的常數(shù)列。為了更方便地進(jìn)行處理，將式（1）進(jìn)行等價(jià)轉(zhuǎn)換：

此時(shí)，將原系數(shù)矩陣的常數(shù)列分離出來(lái)，而線性回歸中的自變量和因變量都可看成是觀測(cè)值。當(dāng)有多組觀測(cè)值時(shí)，式（2）可表示為矩陣形式：

式（3）即為等價(jià)轉(zhuǎn)換后的線性回歸總體最小二乘平差模型。其中，A為由因變量和自變量組成的m×n觀測(cè)矩陣，E為A的誤差矩陣，W為由常數(shù)1組成的m×1 常數(shù)向量，X為n×1 待估參數(shù)。根據(jù)總體最小二乘原理，在加權(quán)情況下相應(yīng)的隨機(jī)模型為：

其中，vec（E）是將矩陣E按列從左到右拉直得到的列向量化矩陣。V是平差模型mn×1 的誤差向量，V＝vec（E）。Q為mn×mn觀測(cè)值的協(xié)因數(shù)陣。根據(jù)文獻(xiàn)［14］，將總體最小二乘平差模型看成是非線性的，則可將式（3）在X＝X0＋處展開得：

顧及EX0＝（（X0）T?Im×n）vec（E）＝FV，AX0＝（（X0）T?Im×n）vec（A）＝FL，其中?為矩陣的克羅內(nèi)克積。可將上式進(jìn)一步表示為：

式（6）即為線性展開后的平差模型，其中L＝vec（A）為由因變量和自變量組成的mn×1的觀測(cè)向量。根據(jù)總體最小二乘原理，其平差準(zhǔn)則為：

根據(jù)平差準(zhǔn)則可構(gòu)造目標(biāo)函數(shù)：

根據(jù)拉格朗日求極值原理，要使上式求得最小值，則φ關(guān)于V和的偏導(dǎo)要等于零，即

將式（9）化簡(jiǎn)整理可得：

根據(jù)式（10），并結(jié)合式（6）可得：

根據(jù)式（11），式（10）的第二式可表示成：

通過(guò)逐次迭代的方法，得到總體最小二乘解。即根據(jù)給定的初值E和X0，可以求解式（12），得到待求參數(shù)和拉格朗日常數(shù)k1，以及誤差向量V。將更新得到的E和X0作為新的初值代入式（12）重新計(jì)算，直到滿足收斂條件為止。

其單位權(quán)中誤差求解公式為：

式中V＝Q（X0?Im×n）k。由式（12）根據(jù)方塊矩陣求逆，可得到迭代到最后一步時(shí)參數(shù)改正數(shù)的表達(dá)式為：

式中，為改正后的觀測(cè)矩陣，即＝A＋E，則根據(jù)協(xié)因數(shù)傳播定律可得：

則參數(shù)精度評(píng)定公式為：

1.2 穩(wěn)健總體最小二乘解算

在只考慮線性回歸模型中變量含有隨機(jī)誤差，而不考慮粗差的情況下，通過(guò)以上的迭代計(jì)算可以求解出線性回歸模型的總體最小二乘解。如果線性回歸模型中變量含有粗差，則可根據(jù)選權(quán)迭代的思想通過(guò)迭代算法來(lái)求得穩(wěn)健總體最小二乘解［13］。本文選取Huber權(quán)函數(shù)：

式中，i為迭代次數(shù)為單位權(quán)中誤差，L（t）為第t個(gè)觀測(cè)值。則線性回歸穩(wěn)健總體最小二乘解算的具體步驟如下。

1）由式（1）并根據(jù)最小二乘原理求得回歸參數(shù)估值a0、a1、…、an，再根據(jù)式（2）將其變換為b0、b1、…、bn，并組成回歸參數(shù)的初值X0＝［b0b1…bn］T。

2）設(shè)E0＝0，考慮自變量與因變量為獨(dú)立等精度，則協(xié)因數(shù)初值Q0＝Im×n。再根據(jù)回歸參數(shù)的初值X0構(gòu)造式（12），由式（12）可以計(jì)算回歸參數(shù)的改正數(shù)和拉格朗日乘數(shù)的初值k，相應(yīng)的X0（i＋1）＝X0（i）＋（i＋1），V（i＋1）＝Q0（X0?Im×n）k（i＋1），E（i＋1）＝vec－1（V（i＋1）），其中vec－1（·）表示vec（·）的逆運(yùn)算，即將mn×1的列向量重新構(gòu)造成m×n的矩陣，再根據(jù)式（17）計(jì)算新的權(quán)陣。

4）輸出參數(shù)估值，按式（13）計(jì)算單位權(quán)中誤差，按式（16）計(jì)算參數(shù)精度。

式中

2 模擬算例及分析

運(yùn)用MATLAB軟件模擬一個(gè)一元線性回歸方程。設(shè)回歸方程為y＝a＋bx，回歸參數(shù)的真值為1.52、2.10，取x為［－10，10］區(qū)間內(nèi)均勻分布的21個(gè)數(shù)，并通過(guò)以上方程計(jì)算y的值，組成21對(duì)點(diǎn)。在21對(duì)點(diǎn)中隨機(jī)選取兩個(gè)點(diǎn)，在其中一個(gè)點(diǎn)的x加入0.9的粗差，y加入－0.9的粗差；另一個(gè)點(diǎn)x加入－0.9的粗差，y加入0.9的粗差。然后在其余各點(diǎn)的x和y上均加入均值為0、方差為0.3 的隨機(jī)誤差。分別采用最小二乘法（LS）、穩(wěn)健最小二乘法（RLS）、總體最小二乘法（TLS）以及本文提出的線性回歸穩(wěn)健總體最小二乘法（LRRTLS）進(jìn)行參數(shù)估計(jì)，進(jìn)行10 000次模擬計(jì)算。然后對(duì)各種方法計(jì)算的結(jié)果進(jìn)行統(tǒng)計(jì)，計(jì)算其均方誤差，以及單位權(quán)中誤差均值結(jié)果見表1。表2為從10 000次模擬計(jì)算中隨機(jī)選擇的一組數(shù)據(jù)組成的觀測(cè)值，其中帶＃的點(diǎn)為加入過(guò)粗差的數(shù)據(jù)。表3為該組數(shù)據(jù)采用不同估計(jì)方法得到的參數(shù)估計(jì)結(jié)果及單位權(quán)中誤差。表4為采用不同估計(jì)方法得到的含粗差點(diǎn)位的改正數(shù)，其中LS和RLS的改正數(shù)為y方向vy，TLS 和LRRTLS的改正數(shù)為x和y方向。

從表1 可以看出，在10 000 次試驗(yàn)中，按LRRTLS法得到的回歸參數(shù)均方誤差和單位權(quán)中誤差均值最小。說(shuō)明在線性回歸中變量含有隨機(jī)誤差同時(shí)又有粗差的情況下，按LRRTLS法求得的回歸參數(shù)與真值相差最小。這是因?yàn)長(zhǎng)RRTLS法通過(guò)選權(quán)迭代對(duì)粗差進(jìn)行探測(cè)定位，能得到較好的穩(wěn)健參數(shù)。LS法只能考慮因變量的隨機(jī)誤差，而不能顧及自變量的誤差，也不能對(duì)含有粗差的數(shù)據(jù)進(jìn)行處理。RLS 法雖然能對(duì)因變量的粗差進(jìn)行處理，但不能顧及自變量的粗差。TLS法只能對(duì)自變量和因變量中的隨機(jī)誤差進(jìn)行處理，而不能對(duì)其含有的粗差進(jìn)行處理。這從表3中也可以看出。表3的平差結(jié)果是以表2的數(shù)據(jù)為基礎(chǔ)的，而表2的數(shù)據(jù)是從10 000次試驗(yàn)中隨機(jī)選出的一組數(shù)據(jù)，具有代表性。表2結(jié)果表明，按LRRLS法得到的回歸參數(shù)與真值相差最小，且其單位權(quán)中差誤差也最小，而另外3種方法得到的回歸參數(shù)與真值偏差較大，這也進(jìn)一步驗(yàn)證了表1的結(jié)論。在加入粗差時(shí)，對(duì)2號(hào)點(diǎn)加入偏離直線右下方距離1.27的誤差，對(duì)15號(hào)點(diǎn)加入偏離直線左上方距離1.27的誤差。使用不同方法對(duì)含粗差點(diǎn)位的改正數(shù)可以從表4 中看出，LRRTLS法得到的改正數(shù)與加入的粗差最接近，而其他3種方法則相差較大。綜上所述，本文提出的LRRTLS法對(duì)解決變量含粗差的線性回歸問(wèn)題是可行的。

表1 不同估計(jì)方法的參數(shù)均方誤差及單位權(quán)中誤差均值Tab.1 Mean square error of parameters with various estimators

表2 觀測(cè)值Tab.2 Observation values

表3 不同估計(jì)方法得到的參數(shù)值及其中誤差Tab.3 Parameters and variance with various estimators

表4 不同估計(jì)方法的改正數(shù)Tab.4 Correction with various estimators

3 結(jié) 語(yǔ)

本文將線性回歸模型進(jìn)行等價(jià)變換，根據(jù)文獻(xiàn)［13］的思想，將其平差模型看成是非線性的，并以選權(quán)迭代為基礎(chǔ)，推導(dǎo)了線性回歸模型的穩(wěn)健總體最小二乘算法。其算法推導(dǎo)過(guò)程簡(jiǎn)單，迭代格式較為方便且容易理解。

通過(guò)模擬算例進(jìn)行分析，對(duì)于因變量和自變量都可能含粗差的線性回歸模型，常規(guī)的最小二乘法（LS）、穩(wěn)健最小二乘法（RLS）、總體最小二乘法（TLS）都不能對(duì)其含有的粗差進(jìn)行合理的處理，得到的參數(shù)與真實(shí)值偏離較大。而本文提出的線性回歸模型的穩(wěn)健總體最小二乘法（LRRTLS），能對(duì)其進(jìn)行較好的處理，得到的參數(shù)較貼近真值，驗(yàn)證了本文方法的有效性和可行性。

［1］Golub G H，Loan C F.An Analysis of the Total Least Squares Problem［J］.SIAM Journal Numerical Analysis，1980，17（6）：883－893

［2］魯鐵定，周世健.總體最小二乘的迭代解法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2010，35（11）：1 351－1 354（Lu Tieding，Zhou Shijian.An Iteration for the Total Least Squares Estimation［J］.Geomatics and Information Science of Wuhan University，2010，35（11）：1 351－1 354）

［3］孔建，姚宜斌，吳寒.整體最小二乘的迭代解法［J］.武漢大學(xué)學(xué) 報(bào)：信 息 科 學(xué) 版，2010，35（6）：711－714（Kong Jian，Yao Yibin，Wu Han.Iterative Method for Total Least－Squares［J］.Geomatics and Information Science of Wuhan University，2010，35（6）：711－714）

［4］許超鈐，姚宜斌，張豹，等.基于整體最小二乘的參數(shù)估計(jì)新方法及精度評(píng)定［J］.測(cè)繪通報(bào)，2011（10）：1－4（Xu Chaoqian，Yao Yibin，Zhang Bao，et al.New Method of Parameters Estimation and Accuracy Evaluation Based on TLS［J］.Bulletin of Surveying and Mapping，2011（10）：1－4）

［5］邱衛(wèi)寧，齊公玉，田豐瑞.整體最小二乘求解線性模型的改進(jìn)算法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2010，35（6）：708－710（Qiu Weining，Qi Gongyu，Tian Fengrui.An Improved Algorithm of Total Least Squares for Linear Models［J］.Geomatics and Information Science of Wuhan University，2010，35（6）：708－710）

［6］邱衛(wèi)寧，陶本藻，姚宜斌，等.測(cè)量數(shù)據(jù)處理理論與方法［M］.武漢：武漢大學(xué)出版社，2008（Qiu Weining，Tao Benzao，Yao Yibin，et al.The Theory and Method of Surveying Data Processing［M］.Wuhan：Wuhan University Press，2008）

［7］丁克良，沈云中，歐吉坤.整體最小二乘法直線擬合［J］.遼寧工程技術(shù)大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，29（1）：44－47（Ding Keliang，Shen Yunzhong，Ou Jikun.Methods of Line－Fitting Based on Total Least－Squares［J］.Journal of Liaoning Technical University：Natural Science，2010，29（1）：44－47）

［8］孫同賀，羅志才.數(shù)字化曲線的整體最小二乘平差［J］.工程勘察，2013（8）：71－73（Sun Tonghe，Luo Zhicai.Total Least Squares Method for Digitized Curve Fitting［J］.Geotechnical Investigation ＆Surveying，2013（8）：71－73）

［9］汪奇生，楊德宏，楊建文.基于總體最小二乘的線性回歸迭代算法［J］.大地測(cè)量與地球動(dòng)力學(xué)，2013，33（6）：112－114（Wang Qisheng，Yang Dehong，Yang Jianwen.An Iteration Algorithm of Linear Regression Based on Total Least Squares［J］.Journal of Geodesy and Geodynamics，2013，33（6）：112－114）

［10］官云蘭，劉紹堂，周世健，等.基于整體最小二乘的穩(wěn)健點(diǎn)云數(shù)據(jù)平面擬合［J］.大地測(cè)量與地球動(dòng)力學(xué)，2011，31（5）：80－83（Guan Yunlan，Liu Shaotang，Zhou Shijian，et al.Robust Plane Fitting of Point Clouds Based on TLS［J］.Journal of Geodesy and Geodynamics，2011，31（5）：80－83）

［11］官云蘭，周世健，張立亭，等.穩(wěn)健整體最小二乘直線擬合［J］.工程勘察，2012，40（2）：60－62（Guan Yunlan，Zhou Shijian，Zhang Liting，et al.A Robust Method for Fitting a Line to Point Clouds Based on TLS［J］.Geotechnical Investigation ＆Surveying，2012，40（2）：60－62）

［12］陳瑋嫻，袁慶.抗差總體最小二乘方法［J］.大地測(cè)量與地球動(dòng)力學(xué)，2012，32（6）：111－113（Chen Weixian，Yuan Qing.A Robust Total Least－Squares Method［J］.Journal of Geodesy and Geodynamics，2012，32（6）：111－113）

［13］陳義，陸玨.以三維坐標(biāo)轉(zhuǎn)換為例解算穩(wěn)健總體最小二乘方法［J］.測(cè)繪學(xué)報(bào)，2012，41（5）：715－722（Chen Yi，Lu Jue.Performing 3DSimilarity Transformation by Robust Total Least Squares［J］.Acta Geodaetica et Cartographica Sinica，2012，41（5）：715－722）

［14］Shen Y，Li B，Chen Y.An Iterative Solution of Weighted Total Least－Squares Adjustment［J］.Journal of Geodesy，2011，85（4）：229－238