程效軍，李 杰

（同濟(jì)大學(xué) 測(cè)繪與地理信息學(xué)院，上海200092）

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

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

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

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

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

1.1 LTS估計(jì)

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

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

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

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

（3）利用步驟（2）得到的平面方程計(jì)算所有點(diǎn)的殘差平方值r2i，i=1，…，n，并排序?yàn)椋?/p>

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

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

權(quán)函數(shù)的選擇對(duì)選權(quán)迭代法的抗差性能十分關(guān)鍵.李德仁等［13］指出，關(guān)于權(quán)函數(shù)的選擇應(yīng)滿(mǎn)足三點(diǎn)要求：①通過(guò)迭代，含粗差觀(guān)測(cè)值的權(quán)應(yīng)逐步趨近于零；②不含粗差觀(guān)測(cè)值的權(quán)，在迭代終止時(shí)應(yīng)等于該組觀(guān)測(cè)值的權(quán)；③權(quán)函數(shù)的選擇應(yīng)保證迭代過(guò)程能以較快的速度收斂.常見(jiàn)的權(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ù)沒(méi)有淘汰段，因而抗差能力減弱.Hampel權(quán)函數(shù)將閾值區(qū)間分為四段，計(jì)算較為復(fù)雜.王任享等［14］指出殘差并不能使改正數(shù)在同一概率水平，標(biāo)準(zhǔn)化殘差屬于同方差母體的子樣，用標(biāo)準(zhǔn)化殘差代替殘差作為權(quán)函數(shù)的參數(shù)在統(tǒng)計(jì)判斷上更為嚴(yán)密.標(biāo)準(zhǔn)化殘差的計(jì)算式為

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

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

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

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

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

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

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

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

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

2 算法驗(yàn)證

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

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

圖1 不含粗差的點(diǎn)云平面圖Fig.1 Point clouds without outliers

圖2 含5%粗差的點(diǎn)云平面圖Fig.2 Point clouds with 5%outliers

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

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

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

圖3 與標(biāo)準(zhǔn)系數(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）

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

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

表5 LTS初值Tab.5 LTS initial value

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

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

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

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

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

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

3 結(jié)語(yǔ)

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

［1］ 蔡來(lái)良，吳侃，張舒.點(diǎn)云平面擬合在三維激光掃描儀變形監(jiān)測(cè)中的應(yīng)用［J］.測(cè)繪科學(xué)，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］ 陳磊，趙書(shū)河.一種改進(jìn)的基于平面擬合的機(jī)載LiDAR點(diǎn)云濾波方法［J］.遙感技術(shù)與應(yīng)用，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)危ㄏ蚣t，陳西江，等.基于正交整體最小二乘平面擬合的點(diǎn)云數(shù)據(jù)去噪方法研究［J］.測(cè)繪通報(bào)，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］ 程效軍，唐劍波.基于最小二乘擬合的墻面平整度檢測(cè)方法［J］.測(cè)繪信息與工程，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］ 張勤，張菊清，岳東杰.近代測(cè)量數(shù)據(jù)處理與應(yīng)用［M］.北京：測(cè)繪出版社，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估計(jì)穩(wěn)健性的比較［J］.中國(guó)衛(wèi)生統(tǒng)計(jì)，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］ 楊飚，張?jiān)疲瑢O政順.基于LTS穩(wěn)健初值的選權(quán)迭代法［J］.科學(xué)技術(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］.武漢：武漢大學(xué)出版社，2002.LI Deren，YUAN Xiuxiao.Error processing and reliability theory［M］.Wuhan：Wuhan University Press，2002.

［14］ 王任享.選權(quán)迭代定位粗差時(shí)權(quán)函數(shù)參數(shù)之功能［J］.武漢測(cè)繪科技大學(xué)學(xué)報(bào)，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.