李明磊，李廣云，王 力，李海波，范哲瑞

(1.信息工程大學(xué)導(dǎo)航與空天目標(biāo)工程學(xué)院，河南鄭州450052;2.68029部隊，甘肅蘭州730000;3.66240部隊，北京 100042)

一、引 言

近年來，激光掃描測量技術(shù)得到了飛速發(fā)展，測量速度從每秒幾千個點(diǎn)提升到百萬量級，測量精度在一定距離內(nèi)可達(dá)到毫米級甚至亞毫米級;其應(yīng)用范圍幾乎涵蓋所有傳統(tǒng)測繪的應(yīng)用領(lǐng)域，包括地形測繪［1］、礦山測量［2］、變形監(jiān)測［3］、文物保護(hù)［4］、虛擬現(xiàn)實［5］、數(shù)字城市［6］、逆向工程［7］等。激光掃描測量技術(shù)大幅減少了外業(yè)數(shù)據(jù)采集的時間和工作量，內(nèi)業(yè)數(shù)據(jù)處理成為項目的工作重心，也是難點(diǎn)所在。對激光點(diǎn)云數(shù)據(jù)進(jìn)行處理、挖掘點(diǎn)云數(shù)據(jù)信息是當(dāng)前激光掃描測量技術(shù)的研究熱點(diǎn)之一。特征提取方法作為二維圖像處理的重點(diǎn)研究對象，在三維激光點(diǎn)云數(shù)據(jù)處理中也存在類似的研究價值，其不僅可以用來進(jìn)行多站點(diǎn)云的快速自動拼接，還可以用于點(diǎn)云的簡化及后續(xù)的建模等任務(wù)。Paul Hough［8］于 1962 年提出 Hough transform(HT)，用于二維圖像中的直線檢測，目前已發(fā)展出多種形式，包括直接由二維HT導(dǎo)出的standard Hough transform(SHT)，以及改進(jìn)的 fast Hough transform(FHT)［9］、probabilistic Hough transform(PHT)［10］和randomized Hough transform(RHT)［11］等。二維 HT 可以容易地擴(kuò)展到三維空間，即三維Hough transform(3D-HT)。本文在分析三維激光點(diǎn)云特點(diǎn)的基礎(chǔ)上，以點(diǎn)云中平面特征的提取為例，研究將3D-HT用于點(diǎn)云特征提取所涉及的關(guān)鍵技術(shù)。

二、Hough transform原理

1.二維HT直線探測原理

在二維圖像空間，每一個像素點(diǎn)在笛卡爾坐標(biāo)系中都有一個明確的二維坐標(biāo)(x，y)，如果圖像中某一特征能用函數(shù)F進(jìn)行數(shù)學(xué)表示，則此特征上的所有像素點(diǎn)都應(yīng)滿足由函數(shù)F變換得到的等式。以二維直線特征為例，函數(shù)F用式(1)表示。式中，θ取任意角度圖像中每一個像素點(diǎn)(x，y)都能在二維Hough空間找到對應(yīng)的Hough變換后的坐標(biāo)(θ，ρ)。其中θ和ρ的幾何意義是:從原點(diǎn)O向過點(diǎn)P的直線F作垂線L，L與x軸之間的夾角即為θ，ρ為原點(diǎn)O到直線F的垂距，如圖1所示，圖中FON(foot of normal)表示L與F的交點(diǎn)。

之后利用二維累加器(θ方向按角度劃分為若干段，ρ方向按距離數(shù)值劃分為若干段)將二維Hough空間的點(diǎn)坐標(biāo)θ和ρ分別按照角度和距離進(jìn)行累加統(tǒng)計，即對于累加器中的任一角度θ1，利用式(1)求出對應(yīng)的 ρ1，累加器中(θ1，ρ1)處進(jìn)行一次累加，這個過程也被稱作投票。如果二維圖像中存在一條直線，其參數(shù)為(θ0，ρ0)，則累加器累加結(jié)果在(θ0，ρ0)處存在一個峰值。每個峰值代表一條直線，因此下一步任務(wù)就是對累加器的累加結(jié)果進(jìn)行局部峰值的探測。峰值探測最常用的方式是通過劃定一個閾值來確定各個峰值的分布，提取出各個票數(shù)高于閾值的峰值。

圖1 二維空間直線的Hough變換

2.三維SHT平面識別原理

對于三維激光點(diǎn)云的特征識別，對二維HT增加一個參數(shù)φ就可以比較容易地擴(kuò)展到三維HT。如果三維空間中某一特征可以用函數(shù)F進(jìn)行數(shù)學(xué)表示，如平面、球、圓柱等，3D-HT就可以以某種適當(dāng)?shù)姆绞綄Υ颂卣鬟M(jìn)行識別。以平面特征為例，在三維笛卡爾坐標(biāo)系中，F(xiàn)一般表示為

為確定某一平面，必須利用這4個參數(shù)a、b、c、d，由于4個參數(shù)的取值范圍是不確定的，直接以這4個參數(shù)進(jìn)行投票時累加器無法確定范圍和投票間隔。因此，3D-HT采用式(3)表示平面特征，將平面參數(shù)轉(zhuǎn)換為角度信息。

式中，參數(shù)θ表示平面的法向n在xoy平面中的投影和x軸正向之間的夾角;φ表示n和xoy平面的夾角;ρ表示原點(diǎn)O距離平面的距離，如圖2所示。

以nφ表示在φ方向?qū)⑷SHough空間劃分的段數(shù)，以nθ表示在θ方向劃分的段數(shù)，以nρ表示在ρ方向劃分的段數(shù)，則三維Hough空間共被三維劃分為nθ×nφ×nρ塊，用作累加器的分區(qū)，以進(jìn)行下一步投票。事實上三維空間中過點(diǎn)P的平面有無窮多個，而由于對三維Hough空間的離散化劃分，平面?zhèn)€數(shù)被限制到 nθ×nφ個(θ和φ 確定后代入式(3)，ρ也可確定，故個數(shù)為 nθ×nφ)，這也就意味著對于每一個點(diǎn)P，需要找到這nθ×nφ個滿足式(3)的平面(按劃分間隔逐個尋找)，即每個點(diǎn)P在三維累加器中需要投票nθ×nφ次。假設(shè)點(diǎn)云中共有N個點(diǎn)，則累加器整個過程就需要N×nθ×nφ次。最終在累加器累加結(jié)果中，局部峰值點(diǎn)表示可能存在的平面。

圖2 三維空間平面的Hough變換

上述過程即為利用三維標(biāo)準(zhǔn)Hough變換(3DSHT)對平面特征進(jìn)行識別。從其過程的實現(xiàn)來看，此方法耗費(fèi)的時間和計算資源較大，因此許多改進(jìn)算法被用于改善標(biāo)準(zhǔn)Hough變換效率低、耗空間等問題。

三、3D-HT平面提取新方法

僅用三維點(diǎn)的笛卡爾坐標(biāo)進(jìn)行3D-SHT平面特征識別還可能會存在一些問題。如圖3所示，內(nèi)層平面和外層平面(就平面本身來說是三維無限的)雖然相交，但事實上卻沒有相交的點(diǎn)，但3D-SHT識別的結(jié)果是將外層平面上與內(nèi)層平面理論相交的點(diǎn)也識別為內(nèi)層平面上的點(diǎn)。為解決這個問題就需要對3D-SHT額外增加相應(yīng)的限制條件，這樣的限制條件并不容易找到，而且適用性不一定廣泛。從3D-SHT原理來看，此方法僅利用各個點(diǎn)的三維坐標(biāo)信息，并且將各個點(diǎn)單獨(dú)考慮，沒有利用各點(diǎn)之間的相互關(guān)系，而且每個點(diǎn)對360°范圍內(nèi)任意傾角φ和θ的平面都要投票(求出任意傾角φ和θ相應(yīng)的距離常數(shù)ρ)，存在變換效率低、耗空間等缺點(diǎn)。事實上對激光點(diǎn)云進(jìn)行處理時至關(guān)重要的一步就是確定各點(diǎn)之間的關(guān)系，以此為基礎(chǔ)才能對點(diǎn)云進(jìn)行下一步分析處理，本文嘗試以此種思路出發(fā)提出3D-HT的新方法。

圖3 3D-SHT平面識別結(jié)果

首先利用數(shù)據(jù)結(jié)構(gòu)快速恢復(fù)點(diǎn)云之間的鄰近關(guān)系，進(jìn)而估計每個點(diǎn)的法矢(法向計算與重定向)。激光點(diǎn)云中每個點(diǎn)的法向是指該點(diǎn)及其臨近的若干點(diǎn)確定的三維表面在此點(diǎn)處的法向。激光點(diǎn)云的法向信息是激光點(diǎn)云數(shù)據(jù)處理中的一個極其重要的信息，對激光點(diǎn)云進(jìn)行法向計算也是其數(shù)據(jù)處理中很重要的一步，之后利用恢復(fù)出的每個點(diǎn)的法向信息進(jìn)行HT的投票。每個點(diǎn)只投票一次，降低投票的次數(shù);利用法向信息進(jìn)行投票，利用到了點(diǎn)與點(diǎn)之間的鄰接關(guān)系，每個點(diǎn)與其鄰近點(diǎn)的法向比較近似，一般不會發(fā)生突變，圖3所示的錯誤識別問題也可解決。

具體過程如下:

1)利用數(shù)據(jù)結(jié)構(gòu)對點(diǎn)云進(jìn)行分塊，從而可以迅速確定點(diǎn)之間的鄰接關(guān)系，加快對每個點(diǎn)鄰近點(diǎn)的搜索速度，本文中數(shù)據(jù)結(jié)構(gòu)采用的是八叉樹［12］。利用八叉樹對某模型的點(diǎn)云數(shù)據(jù)進(jìn)行分塊的效果如圖4所示(視圖的方向逆著y軸方向，故只能看到二維效果)。

圖4 八叉樹分割效果

2)利用數(shù)據(jù)結(jié)構(gòu)尋找點(diǎn)P0的k個鄰近點(diǎn)，用這些點(diǎn)進(jìn)行最小二乘平面擬合(按式(2)進(jìn)行)［13］，并依次以此方式計算各個點(diǎn)的法向。算得的法向方向具有二向性(可能與實際情況相反)，通過從某點(diǎn)出發(fā)進(jìn)行鄰近傳遞可對其進(jìn)行調(diào)整。而對于激光點(diǎn)云，每個點(diǎn)的法向n=(nx，ny，nz)T必定逆向激光的入射方向r=(x，y，z)T，因此只需保證式(4)中 d為負(fù)，即可調(diào)整法向到正確指向。

3)將各個點(diǎn)的法向轉(zhuǎn)換為式(5)所示的以角度和距離參數(shù)表示的形式，以利于投票過程的執(zhí)行。

4)確定累加器的劃分間隔，利用各點(diǎn)變換得到的(θ，φ，ρ)進(jìn)行投票。

5)局部峰值探測，從而找出點(diǎn)云中可能存在的各個平面。

6)利用歸類到各個平面上的點(diǎn)集重新按照式(2)進(jìn)行抗差最小二乘平面擬合［13］，以使提取的平面具有更好的平面度。

從執(zhí)行過程可以看出，該方法比較充分地利用了激光點(diǎn)云預(yù)處理得到的成果(鄰近關(guān)系的確立、法向的計算與重定向等)，在步驟4)的投票過程中假設(shè)點(diǎn)云共有N個點(diǎn)，此方法也就僅需要投票N次，遠(yuǎn)遠(yuǎn)少于3D-SHT的N×nθ×nφ次，大大降低了投票過程的復(fù)雜度。

由于法向計算方式的影響，在墻角等特征突變點(diǎn)附近的法向計算結(jié)果一般是不正確的，因此墻角等處的點(diǎn)集也就不能通過HT進(jìn)行有效識別。本文方法在步驟6)重新進(jìn)行了抗差最小二乘平面擬合，得到的不僅有平面法向參數(shù)(a、b、c)，還有平面到原點(diǎn)的距離參數(shù)d，即激光三維點(diǎn)云笛卡爾坐標(biāo)系下完整的平面方程式，在平面上的點(diǎn)一定滿足相應(yīng)的平面方程。據(jù)此可對法向計算結(jié)果錯誤的特征突變點(diǎn)及其附近點(diǎn)等進(jìn)一步歸類，分離出墻角點(diǎn)等特征突變點(diǎn):

1)由于墻角點(diǎn)滿足步驟6)結(jié)果中兩個以上的平面方程(平面的交點(diǎn))，依此可以有效地將其分離出。

2)將墻角點(diǎn)附近點(diǎn)的坐標(biāo)分別代入步驟6)平面方程式，可以將其重新歸類到其應(yīng)在的平面，并且其法向計算結(jié)果也可以得到修正。

四、試驗驗證

首先是累加器劃分間隔的具體確定方法。本文利用Rigel公司的VZ-400掃描儀對某實驗室的普通墻面(一般建筑內(nèi)部的常用墻面具有較好的代表性)進(jìn)行精掃，如圖5所示。圖5(a)是對墻面拍攝的照片;圖5(b)中上圖為在掃描結(jié)果中選取墻面比較平整的點(diǎn)云塊，下圖為點(diǎn)云法向計算結(jié)果添加光照顯示的效果;圖5(c)為將對每個點(diǎn)算得的單位化后的法向(包含3個方向的分量)作為三維坐標(biāo)顯示(分布在半徑為1的球面上)的效果，由于墻面不平整、激光掃描儀誤差、法向計算誤差等因素的影響，對墻面點(diǎn)云算得的法向并不完全相同(若完全相同，圖5(c)所有點(diǎn)都重合為一個點(diǎn))。

圖5 實驗室墻面照片與點(diǎn)云

對墻面的法向計算結(jié)果進(jìn)行統(tǒng)計，結(jié)果見表1。由于y軸與墻面的夾角最大，接近垂直，故其方向誤差較小。由于 sin(1°)=0.017 5，對于兩個單位向量，坐標(biāo)差別0.017 5相當(dāng)于夾角差別接近1°。根據(jù)法向計算的中誤差將累加器中角度的劃分間隔設(shè)置為1°(法向夾角小于1°的兩個點(diǎn)當(dāng)作處于兩個平行平面上的點(diǎn)，根據(jù)距離參數(shù)確定兩平面是否重合)，根據(jù)實際平面特征分布，將距離劃分間隔設(shè)為5 cm(距離小于5 cm的兩平行平面當(dāng)作同一平面，間隔太小無意義且增加空間消耗，太大則不能有效識別平面特征)。

表1 墻面法向中誤差

為測試本文方法的有效性，首先利用兩層立方體模擬數(shù)據(jù)(共12個平面，303 612個點(diǎn))進(jìn)行試驗。圖6即為平面識別的效果，不同的平面之間用不同的顏色來區(qū)分。圖6(a)所示，由于一般的鄰近點(diǎn)平面擬合不能有效計算得到正確的法向，故在平面相交處的點(diǎn)(邊角點(diǎn))并未被有效識別。而通過本文特征突變點(diǎn)及其附近點(diǎn)重新歸類的思路(修正方案)，可以有效地對其重新歸類，將平面和墻角分離，效果如圖6(b)所示。表2為對模擬數(shù)據(jù)進(jìn)行平面提取的輸出結(jié)果(轉(zhuǎn)換成式(2)的4個參數(shù)，并記錄平面上識別到的點(diǎn)數(shù))，與圖6數(shù)據(jù)模擬時的設(shè)定值基本符合(內(nèi)層平面完全符合，外層平面由于內(nèi)層平面與其交線的影響有個別點(diǎn)仍未有效歸類)。

圖6 兩層立方體數(shù)據(jù)平面提取效果

表2 兩層立方體數(shù)據(jù)平面提取結(jié)果

為測試本文方法的實用性，利用某實驗室內(nèi)部激光掃描實測數(shù)據(jù)進(jìn)行試驗。圖7即為平面識別的效果，圖7(a)為僅依據(jù)法向投票識別的效果，圖7(b)為利用修正方案識別的效果，圖7(c)為將分離出的墻角顯示的效果。最終平面提取結(jié)果為6個，與實際符合，而且墻面墻角特征可以有效分離。

圖7 激光實測數(shù)據(jù)平面提取效果

五、結(jié)束語

本文根據(jù)激光點(diǎn)云數(shù)據(jù)處理的特點(diǎn)，結(jié)合激光點(diǎn)云數(shù)據(jù)預(yù)處理中點(diǎn)與點(diǎn)之間鄰近關(guān)系的確立及法向的計算與重定向，提出了直接根據(jù)每個點(diǎn)法向進(jìn)行Hough變換的方法，極大地降低了3D HT的復(fù)雜度，并可以有效解決3D-HT存在的錯誤識別問題。對特征突變點(diǎn)及其附近點(diǎn)的重新歸類方法可以進(jìn)一步識別大部分沒有正確識別的特征突變點(diǎn)，并可以對其法向計算結(jié)果進(jìn)行修正。本文提出的方法僅受法向計算結(jié)果制約，具有良好的適用性。試驗中首先利用普通常見墻面的激光掃描儀精掃數(shù)據(jù)確定了三維累加器的劃分間隔，之后分別利用模擬數(shù)據(jù)和實測數(shù)據(jù)證明了該方法的有效性和實用性。

本文的工作主要集中于激光三維點(diǎn)云平面特征及相應(yīng)墻角輪廓的提取，下一步可以以此為基礎(chǔ)研究更多種類特征的提取方法。

［1］許映林.3維激光掃描技術(shù)在溫泉水電站大比例尺地形圖測量中的應(yīng)用［J］.測繪通報，2007(6):9-13.

［2］楊青，張亮，王振.基于KD樹和CDT的露天礦三維建模［J］.計算機(jī)技術(shù)與發(fā)展，2011，21(12):224-226.

［3］譚國銓.3D激光掃描儀在電廠冷卻塔變形監(jiān)測中的應(yīng)用探討［J］.勘測設(shè)計，2009(4):30-32.

［4］周俊召，鄭書民，胡松，等.地面三維激光掃描在石窟石刻文物保護(hù)測繪中的應(yīng)用［J］.測繪通報，2008(11):68-69.

［5］李暉，吳祿慎.三維激光掃描技術(shù)在虛擬現(xiàn)實中的應(yīng)用［J］.南昌大學(xué)學(xué)報:工科版，2007，29(3):239-242.

［6］趙海瑩，張正鵬.三維激光點(diǎn)云數(shù)據(jù)在城市建模中的應(yīng)用［J］.城市勘測，2009(1):69-72.

［7］馬有良，李占峰.三維激光掃描儀在曲面逆向工程中的應(yīng)用研究［J］.機(jī)械設(shè)計與制造，2009(10):85-87.

［8］HOUGH PV C.Method and Means for Recognizing Complex Patterns:US，Patent 3069654［P］.1962-12-18.

［9］OGUNDANA OO，COGGRAVE CR，BURGUETE R L，et al.Automated Detection of Planes in 3-D Point Clouds Using Fast Hough Transforms［J］.Optical Engineering，2011，50(5):053609.

［10］KIRYATI N，ELDAR Y，BRUCKSTEIN A M.A Probabilistic Hough Transform ［J］.Pattern Recognition，1991，24(4):303-316.

［11］XU L，OJA E，KULTANEN P.A New Curve Detection Method:Randomized Hough Transform(RHT)［J］.Pattern Recognition Letters，1990(11):331-338.

［12］邢坤.海量激光掃描測量數(shù)據(jù)的處理［D］.鄭州:信息工程大學(xué)，2009.

［13］李明磊.點(diǎn)云平面擬合新方法［J］.測繪通報，2012(S1):84-87.