線性回歸模型的穩健總體最小二乘解算

汪奇生 楊德宏 楊騰飛

1 湖南軟件職業學院，湘潭市寶馬西路，411100

2 昆明理工大學國土資源工程學院，昆明市文昌路68號，650093

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

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

1 穩健總體最小二乘算法原理

1.1 總體最小二乘算法

線性回歸數學模型為：

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

此時，將原系數矩陣的常數列分離出來，而線性回歸中的自變量和因變量都可看成是觀測值。當有多組觀測值時，式（2）可表示為矩陣形式：

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

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

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

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

根據平差準則可構造目標函數：

根據拉格朗日求極值原理，要使上式求得最小值，則φ關于V和的偏導要等于零，即

將式（9）化簡整理可得：

根據式（10），并結合式（6）可得：

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

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

其單位權中誤差求解公式為：

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

式中，為改正后的觀測矩陣，即＝A＋E，則根據協因數傳播定律可得：

則參數精度評定公式為：

1.2 穩健總體最小二乘解算

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

式中，i為迭代次數為單位權中誤差，L（t）為第t個觀測值。則線性回歸穩健總體最小二乘解算的具體步驟如下。

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

2）設E0＝0，考慮自變量與因變量為獨立等精度，則協因數初值Q0＝Im×n。再根據回歸參數的初值X0構造式（12），由式（12）可以計算回歸參數的改正數和拉格朗日乘數的初值k，相應的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（·）的逆運算，即將mn×1的列向量重新構造成m×n的矩陣，再根據式（17）計算新的權陣。

4）輸出參數估值，按式（13）計算單位權中誤差，按式（16）計算參數精度。

式中

2 模擬算例及分析

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

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

表1 不同估計方法的參數均方誤差及單位權中誤差均值Tab.1 Mean square error of parameters with various estimators

表2 觀測值Tab.2 Observation values

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

表4 不同估計方法的改正數Tab.4 Correction with various estimators

3 結 語

本文將線性回歸模型進行等價變換，根據文獻［13］的思想，將其平差模型看成是非線性的，并以選權迭代為基礎，推導了線性回歸模型的穩健總體最小二乘算法。其算法推導過程簡單，迭代格式較為方便且容易理解。

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

［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］.武漢大學學報：信息科學版，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］.武漢大學學 報：信 息 科 學 版，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］許超鈐，姚宜斌，張豹，等.基于整體最小二乘的參數估計新方法及精度評定［J］.測繪通報，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］邱衛寧，齊公玉，田豐瑞.整體最小二乘求解線性模型的改進算法［J］.武漢大學學報：信息科學版，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］邱衛寧，陶本藻，姚宜斌，等.測量數據處理理論與方法［M］.武漢：武漢大學出版社，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］.遼寧工程技術大學學報：自然科學版，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］孫同賀，羅志才.數字化曲線的整體最小二乘平差［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］.大地測量與地球動力學，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］官云蘭，劉紹堂，周世健，等.基于整體最小二乘的穩健點云數據平面擬合［J］.大地測量與地球動力學，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］官云蘭，周世健，張立亭，等.穩健整體最小二乘直線擬合［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］.大地測量與地球動力學，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］陳義，陸玨.以三維坐標轉換為例解算穩健總體最小二乘方法［J］.測繪學報，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