具有最小截斷二乘穩(wěn)健初值的點云平面擬合算法

程效軍，李 杰

（同濟大學 測繪與地理信息學院，上海200092）

三維激光掃描作為一種新型的測量技術(shù)，可以快速獲取對象表面的點云數(shù)據(jù)，為目標的三維重建提供了一種便捷方法.真實場景中含有大量平面特征，點云平面擬合是三維重建的基礎步驟，在此基礎上衍生了眾多應用.蔡來良等［1］利用點云平面擬合方法進行了建筑物變形監(jiān)測應用研究；陳磊［2］和葉珉?yún)危?］提出了基于點云平面擬合的點云濾波算法；程效軍［4］利用最小二乘平面擬合方法來檢測墻面平整度.因此，探討點云平面擬合的穩(wěn)健算法對三維激光掃描的應用具有重要意義.

對于含有粗差的點云平面擬合，目前常用的方法有粗差探測法、選權(quán)迭代法、穩(wěn)健特征值法和隨機采樣一致性算法等.粗差探測法在假定只有一個觀測粗差的前提下其理論上是嚴密的，但實際獲取的點云數(shù)據(jù)往往不止一個粗差.若逐次采用粗差探測法來剔除粗差，則存在以下問題：由于粗差對每個觀測值都有影響，尤其當存在多個粗差時，第一步數(shù)據(jù)探測中統(tǒng)計量最大的那個觀測值不一定含有粗差，如果將其剔除就會造成錯誤的判斷［5］，這就造成了數(shù)據(jù)探測法應用的局限性.選權(quán)迭代法是目前測量數(shù)據(jù)處理中剔除粗差點的穩(wěn)健估計方法，可以保證所估計的參數(shù)少受模型誤差尤其是粗差的影響.一般通過最小二乘估計來確定第一次平差后的殘差，由于極少的粗差也可能使最小二乘方法崩潰，因此，在一定程度上削弱了選權(quán)迭代法的抗差能力，而且隨著粗差含量的增加，選權(quán)迭代法得到的擬合平面將明顯偏離真實值.穩(wěn)健特征值法能抵抗一定含量的粗差影響，得到較為準確的平面.隨機采樣一致性算法的采樣次數(shù)視具體數(shù)據(jù)而定，沒有固定的參考值，而且恰當?shù)拈撝惦y以確定.

最小截斷二乘法（least trimmed squares，LTS）是一種穩(wěn)健估計方法，常用于統(tǒng)計學中處理含有局外點的線性回歸問題［6-9］.LTS對局外點的影響不敏感，可以達到50%的崩潰點［10］.與LTS方法相近的還有最小中位二乘法（least median squares，LMS），Rousseeuw［8］指出，LTS回歸算法比 LMS更有優(yōu)勢：LTS的目標函數(shù)更平滑，使其對局部影響較??；其次LTS估值近似正態(tài)，收斂速度比LMS更快，因而LTS的統(tǒng)計效果更好.Doornik［9］討論了樣本量與計算時間的關(guān)系，當樣本量為105時，計算時間為6 800S.因此當數(shù)據(jù)量較大時，LTS估計效率不高.在處理含有粗差的點云平面擬合問題上，本文嘗試在選權(quán)迭代法的基礎上結(jié)合LTS，以提高選權(quán)迭代法的穩(wěn)健性，同時保持較高的估計效率.

1 具有LTS穩(wěn)健初值的點云平面擬合算法

常用的選權(quán)迭代法能抵抗點云中一定含量的粗差，當粗差含量達到一定值后，選權(quán)迭代法便不能得到正確的結(jié)果.為了得到統(tǒng)計意義上更為穩(wěn)健的估計方法，同時保持較高的估計效率，考慮將LTS估值作為初始值，計算第一次平差后的殘差，選擇較穩(wěn)健 的 絕 對 偏 差 中 位 數(shù) MAD（the median of all absolute deviations）作為第一次計算的單位權(quán)中誤差，然后再利用此殘差進行選權(quán)迭代，這樣既保留了選權(quán)迭代法較高的估計效率，又繼承了LTS的高崩潰點，因而在一定程度上改善了選權(quán)迭代法的估計性能.

1.1 LTS估計

Rousseeuw［11］提 出 LTS 的 定 義 式 為 arg min（本文中帶括號的下標表示排 序后的 順序），將平方后的殘差r2按照從小到大的順序進行排列后得到，即r2（1）≤…≤r2（i）≤…≤r2（h）.其中，r為殘差，h=int［（n+p+1）／2］，n為點的個數(shù)，p為待估參數(shù)個數(shù).LTS的主要思想是使目標函數(shù)升序排列的前h個殘差平方值之和達到最小［12］.LTS估計當h=int［（n+p+1）／2］時其崩潰點在合理情況下可達到50%，是一種非常穩(wěn)健的方法.

遍歷n個觀測值中任意h個觀測值的所有可能組合，對每種組合進行最小二乘估計，取其中最小殘差平方和對應的解作為LTS估計的準確解［6］.從n個觀測值中任取h個觀測值一共有Chn種可能，當n和h較大時，計算量非常大，因而也導致了LTS的估計效率不高.對此，為提高計算效率，王彤等［7］提出了重復抽樣算法以得到LTS的近似解.Rousseeuw等［8］利用選擇迭代和數(shù)據(jù)分塊的策略提出了快速LTS方法來提高LTS的效率.由于點云數(shù)據(jù)量較大，因此采用王彤等［7］提出的方法得到LTS近似解以提高計算效率.步驟如下：

（1）確定重復抽樣次數(shù)SN；

（2）從n個觀測值中任取p+1個觀測值（p為待估參數(shù)個數(shù)，如平面方程為ax+by+cz=1，其中a，b，c為平面方程系數(shù)，則p=3）進行平面擬合，得到相應的平面方程系數(shù)；

（3）利用步驟（2）得到的平面方程計算所有點的殘差平方值r2i，i=1，…，n，并排序為；

（5）重復步驟（2）～（4）直至抽樣次數(shù)SN，在所得到的SN個殘差平方和中取最小值，此時對應的平面方程系數(shù)即作為擬合平面的LTS估計.

1.2 權(quán)函數(shù)的確定

權(quán)函數(shù)的選擇對選權(quán)迭代法的抗差性能十分關(guān)鍵.李德仁等［13］指出，關(guān)于權(quán)函數(shù)的選擇應滿足三點要求：①通過迭代，含粗差觀測值的權(quán)應逐步趨近于零；②不含粗差觀測值的權(quán)，在迭代終止時應等于該組觀測值的權(quán)；③權(quán)函數(shù)的選擇應保證迭代過程能以較快的速度收斂.常見的權(quán)函數(shù)有Huber權(quán)函數(shù)、Hampel權(quán)函數(shù)、丹麥權(quán)函數(shù)、IGG（institute of geodesy and geophysics）權(quán)函數(shù)等.其中 Huber權(quán)函數(shù)和丹麥權(quán)函數(shù)沒有淘汰段，因而抗差能力減弱.Hampel權(quán)函數(shù)將閾值區(qū)間分為四段，計算較為復雜.王任享等［14］指出殘差并不能使改正數(shù)在同一概率水平，標準化殘差屬于同方差母體的子樣，用標準化殘差代替殘差作為權(quán)函數(shù)的參數(shù)在統(tǒng)計判斷上更為嚴密.標準化殘差的計算式為

式中：ωi為標準化殘差；ri為殘差；gii為標準化因子，其值為殘差協(xié)因數(shù)陣Qvv的主對角線元素；σ0為單位權(quán)中誤差.因此采用改進的IGG權(quán)函數(shù).

式中：P（ω）為權(quán)函數(shù)；ω為標準化殘差；k為常數(shù).

第一次平差計算時單位權(quán)中誤差可以采用Rousseeuw［15］推薦的尺度估計S：

式中：1.482 6為調(diào)制因子；n為觀測值個數(shù)；p為待估參數(shù)個數(shù)；r為殘差；median為r2i（i=1，2，…，n）的中值.也可以采用較穩(wěn)健的絕對偏差中位數(shù)MAD作為初始的單位權(quán)中誤差M［16］.

其中，1.483為改正因子，使MAD在正態(tài)分布情況下無偏，r為殘差，median表示取中值.本文采用MAD作為初始單位權(quán)中誤差.綜上所述，具有LTS穩(wěn)健初值的點云平面擬合算法的計算步驟如下：

（1）計算點云平面方程系數(shù)的LTS初值，迭代次數(shù)i=1，得到殘差r（1），并計算初始單位權(quán)中誤差M；

（2）計算標準化殘差ω，根據(jù)IGG權(quán)函數(shù)確定權(quán)陣P；

（3）迭代次數(shù)i=i+1，計算平面方程系數(shù)的加權(quán)最小二乘估計p（i）p，并計算殘差r（i）；

（6）輸出平面方程系數(shù)的最佳估值pp.

2 算法驗證

為了驗證本文算法的可靠性和穩(wěn)健性，利用地面三維激光掃描儀掃描一平面，截取其中的一部分作為實驗數(shù)據(jù)，所有算法采用Matlab編程實現(xiàn).為檢驗擬合平面的準確性，在點云處理軟件Geomagic Studio中進行去噪，得到凈化的點云數(shù)據(jù)（圖1），再根據(jù)最小二乘方法計算得到標準平面方程為：

對粗差的模擬，利用Matlab向凈化的點云數(shù)據(jù)中隨機添加粗差，粗差含量從5%逐次添加到40%（如圖2所示為隨機添加5%粗差后的點云）.首先采用選權(quán)迭代法和特征值迭代法對不同粗差含量的點云數(shù)據(jù)進行擬合，結(jié)果見表1和表2，表中，a，b，c為平面方程系數(shù)，σ為中誤差.

圖1 不含粗差的點云平面圖Fig.1 Point clouds without outliers

圖2 含5%粗差的點云平面圖Fig.2 Point clouds with 5%outliers

表1 選權(quán)迭代法的擬合結(jié)果Tab.1 Fitting result of selecting weight iteration

根據(jù)圖表可以得出：當粗差含量在15%以內(nèi)時，選權(quán)迭代法和特征值迭代法均穩(wěn)健，當粗差含量達到20%時，選權(quán)迭代法和特征值迭代法的計算結(jié)果與標準系數(shù)的偏差分別為21.99%和7.04%，迭代崩潰.為定出這兩種方法的迭代崩潰點，對粗差含量為16%～19%的點云再次進行測試.試驗結(jié)果見表3和表4.若以偏差2%為限，則選權(quán)迭代法的崩潰點可定為18%，特征值迭代法的崩潰點可定為17%.將表1和表2中的數(shù)據(jù)繪制成圖，如圖3所示.

表2 特征值迭代法的擬合結(jié)果Tab.2 Fitting result of eigenvalue iteration

圖3 與標準系數(shù)的偏差及迭代次數(shù)的變化Fig.3 Deviation from standard plane and iteration times variation

表3 選權(quán)迭代法的擬合結(jié)果（16%～19%粗差含量）Tab.3 Fitting result of selecting weight iteration（With 16%～19%outliers）

表4 特征值迭代法的擬合結(jié)果（16%～19%粗差含量）Tab.4 Fitting result of eigenvalue iteration（With 16%～19%outliers）

試驗結(jié)果表明，隨著粗差含量的增加，迭代計算并不能取得滿意的結(jié)果，當粗差含量達到崩潰點后，由于初始迭代模型偏差過大，導致迭代失??；兩種方法解算出的系數(shù)與標準系數(shù)的偏差總體趨勢是增加的，但由于粗差的隨機性，在某些點的偏差有所減小；選權(quán)迭代法的迭代次數(shù)隨著粗差含量的增加先增加后減小，而特征值迭代法的迭代次數(shù)相對較穩(wěn)定.因此，當粗差含量在崩潰點以內(nèi)時，選權(quán)迭代法和特征值迭代法都可以得到準確的結(jié)果.在實際應用中，由于點云的粗差含量事先無法預知，為了驗證本文算法的有效性，對算法的測試從粗差含量為5%開始.

在計算LTS初值時，首先需要確定重復抽樣次數(shù).對平面擬合來說，每次應從所有觀測值中隨機抽取4個觀測值進行計算.設粗差含量為ν，則抽樣一次抽到粗差的概率為ν，隨機抽取u次，則抽到不含粗差的樣本的概率為P=1-νu，當樣本的粗差含量為40%時，取u=40，計算可得P≈1，因此在測試數(shù)據(jù)中確定重復抽樣次數(shù)為40.LTS初值結(jié)果見表5.

表5 LTS初值Tab.5 LTS initial value

由于LTS初值與標準平面系數(shù)仍有偏差，因此在LTS初值基礎上再進行選權(quán)迭代，對模型進行精煉.最終獲得平面方程系數(shù)見表6.

根據(jù)表6可以得出：對于粗差含量較低的點云，具有LTS初值的選權(quán)迭代法可以得到準確的平面方程系數(shù)，其精度與選權(quán)迭代法和特征值迭代法相當.當粗差含量為20%～40%時，本文算法依然穩(wěn)健，且點云平面擬合系數(shù)與標準系數(shù)的偏差均小于2%，滿足精度要求.在整個計算過程中，最大迭代次數(shù)為34次.因此，具有LTS初值的選權(quán)迭代法較上述兩種迭代法具有更高的崩潰點，更為穩(wěn)健，且保證了計算效率.

表6 具有LTS初值的選權(quán)迭代法的擬合結(jié)果Tab.6 Fitting result of selecting weight iteration with LTS initial value

為探討20%和35%處出現(xiàn)的偏差波動現(xiàn)象是由粗差的隨機性引起，還是由系統(tǒng)誤差造成，選擇另一平面按照具有LTS初值的選權(quán)迭代法進行試驗.標準平面方程為：0.023 910 4x+0.133 495 1y-0.003 445z=1，試驗結(jié)果見表7.

表7 具有LTS初值的選權(quán)迭代法的擬合結(jié)果Tab.7 Fitting result of selecting weight iteration with LTS initial value

根據(jù)表7可知，點云平面擬合系數(shù)與標準系數(shù)的偏差均小于2%，滿足精度要求.隨著粗差含量的增加，偏差值并未呈現(xiàn)相應的增大，而是呈現(xiàn)不規(guī)律的波動，故推斷由于粗差的隨機性而導致偏差變化的不規(guī)律.

3 結(jié)語

對于粗差含量不高的點云平面擬合問題，采用選權(quán)迭代法、特征值迭代法和具有LTS初值的選權(quán)迭代法都可以獲得理想的結(jié)果.當粗差含量增加到一定值后，由于初始迭代模型偏差過大，導致選權(quán)迭代法和特征值迭代法失效.具有LTS初值的選權(quán)迭代法可以獲得穩(wěn)健的平面初值，在此基礎上再進行選權(quán)迭代，這樣既采納了LTS的高崩潰點，又利用了選權(quán)迭代法的高估計效率，從而獲得較為精確的平面方程系數(shù).在進行LTS初值計算時采用了隨機抽樣的方法求得近似解，以提高計算效率，這對于數(shù)據(jù)量龐大的點云來說十分重要.實驗表明，當粗差含量達到很高的比例（如40%）時，具有LTS初值的選權(quán)迭代法依然穩(wěn)健.由于本文算法采用的平面方程為通式ax+by+cz=1，因此該算法可用于水平或鉛垂等特殊點云平面擬合.如果將平面模型改為曲線曲面等模型，則具有LTS穩(wěn)健初值的選權(quán)迭代法還可以應用到其他工程數(shù)據(jù)擬合方面.

［1］ 蔡來良，吳侃，張舒.點云平面擬合在三維激光掃描儀變形監(jiān)測中的應用［J］.測繪科學，2010，35（5）：231.CAI Lailiang，WU Kan，ZHANG Shu.Application of point cloud plane fitting to deformation monitoring using 3D laser scanner［J］.Science of Surveying and Mapping，2010，35（5）：231.

［2］ 陳磊，趙書河.一種改進的基于平面擬合的機載LiDAR點云濾波方法［J］.遙感技術(shù)與應用，2011，26（1）：117.CHEN Lei，ZHAO Shuhe.An improved plane fitting based filtering algorithm for airborne LiDAR data［J］.Remote Sensing Technology and Application，2011，26（1）：117.

［3］ 葉珉?yún)?，花向紅，陳西江，等.基于正交整體最小二乘平面擬合的點云數(shù)據(jù)去噪方法研究［J］.測繪通報，2013（11）：37.YE Minlv，HUA Xianghong，CHEN Xijiangetal.Research on method for the denoising of point cloud based on orthogonal total least squares fitting［J］.Bulletin of Surveying and Mapping，2013（11）：37.

［4］ 程效軍，唐劍波.基于最小二乘擬合的墻面平整度檢測方法［J］.測繪信息與工程，2007，32（4）：19.CHENG Xiaojun，TANG Jianbo.Method for estimating metope smoothing grade based on least squares fitting［J］.Journal of Geomatics，2007，32（4）：19.

［5］ 張勤，張菊清，岳東杰.近代測量數(shù)據(jù)處理與應用［M］.北京：測繪出版社，2011.ZHANG Qin，ZHANG Juqing，YUE Dongjie.Advanced theory and application of surveying data［M］.Beijing：Surveying and Mapping Press，2011.

［6］ Giloni A，Padberg M.Least trimmed squares regression，least median squares regression，and mathematical programming［J］.Mathematical and Computer Modeling，2002，35（9）：1043.

［7］ 王彤，何大衛(wèi).LTS回歸與M估計穩(wěn)健性的比較［J］.中國衛(wèi)生統(tǒng)計，1999，16（2）：79.WANG Tong，HE Dawei.Comparison of M-estimator and LTS estimator in regression ［J］.Chinese Journal of Health Statistics，1999，16（2）：79.

［8］ Rousseeuw P J，Van Driessen K.Computing LTS regression for large data sets［J］.Data Mining and Knowledge Discovery，2006，12（1）：29.

［9］ Doornik J A.Robust estimation using least trimmed squares［EB／OL］.（2011-03-14）［2014-08-08］.http：／／econ.au.dk／fileadmin／site＿files／filer＿oekonomi／subsites／creates／Seminar＿Papers／2011／ELTS.pdf.

［10］ Mount D M，Netanyahu N S，Piatko C D，etal.On the least trimmed squares estimator［J］.Algorithmica，2014，69（1）：148.

［11］ Rousseeuw P J.Least median of squares regression［J］.Journal of the American Statistical Association，1984，79（388）：871.

［12］ 楊飚，張曾科，孫政順.基于LTS穩(wěn)健初值的選權(quán)迭代法［J］.科學技術(shù)與工程，2005，5（22）：1720.YANG Biao，ZHANG Zengke，SUN Zhengshun.Selecting weight iteration method with initial value by LTS［J］.Science Technology and Engineering，2005，5（22）：1720.

［13］ 李德仁，袁修孝.誤差處理與可靠性理論［M］.武漢：武漢大學出版社，2002.LI Deren，YUAN Xiuxiao.Error processing and reliability theory［M］.Wuhan：Wuhan University Press，2002.

［14］ 王任享.選權(quán)迭代定位粗差時權(quán)函數(shù)參數(shù)之功能［J］.武漢測繪科技大學學報，1988，13（4）：42.WANG Renxiang.Effects of parameters of weight function for blunder detection by the recursive weighted least squares method［J］.Journal of Wuhan Technical University of Surveying and Mapping，1988，13（4）：42.

［15］ Rousseeuw P J，Leroy A M.Robust regression and outlier detection［M］.New Jersey：John Wiley &Sons，2005.

［16］ Rousseeuw P J，Hubert M.Robust statistics for outlier detection［J］.Wiley Interdisciplinary Reviews：Data Mining and Knowledge Discovery，2011，1（1）：73.