曲率約束的激光點(diǎn)云全局優(yōu)化配準(zhǔn)算法

馬偉麗，2，王健，孫文瀟

(1.山東科技大學(xué)，山東 青島 266590；2.山東省水利勘測(cè)設(shè)計(jì)院，濟(jì)南 250101)

0 引言

在物體表面三維重構(gòu)[1]中，點(diǎn)云數(shù)據(jù)配準(zhǔn)是最為重要的步驟。使用三維激光掃描技術(shù)[2-3]，可以得到被測(cè)物體表面點(diǎn)云數(shù)據(jù)，通過(guò)在多個(gè)視角下測(cè)量的點(diǎn)云數(shù)據(jù)，然后對(duì)多測(cè)站數(shù)據(jù)進(jìn)行配準(zhǔn)，最終獲得完整的三維物體模型。點(diǎn)云數(shù)據(jù)配準(zhǔn)精度，對(duì)三維物體模型重構(gòu)和應(yīng)用有著直接的影響。

目前，國(guó)內(nèi)外提出了很多數(shù)據(jù)配準(zhǔn)的算法，其中使用最為廣泛的算法就是1992年Besl等[4]提出的迭代最近點(diǎn)算法(ite-rative closest point，ICP)以及它的改進(jìn)算法[5-6]。為了使數(shù)據(jù)配準(zhǔn)的結(jié)果更加準(zhǔn)確，ICP算法使用時(shí)必須要滿足2個(gè)條件：點(diǎn)云之間存在包含關(guān)系、點(diǎn)云之間初始位置不可以相差過(guò)大。數(shù)據(jù)配準(zhǔn)常用方法有：特征點(diǎn)提取法[7-10]，根據(jù)曲率特征、法向量等信息提取特征點(diǎn)，由得到的特征點(diǎn)對(duì)完成配準(zhǔn)；質(zhì)心重合法[11]，計(jì)算點(diǎn)集的質(zhì)心，并將2個(gè)點(diǎn)集作相對(duì)于質(zhì)心的平移；標(biāo)志點(diǎn)法[12-13]，人工粘貼標(biāo)志點(diǎn)，通過(guò)標(biāo)志點(diǎn)配準(zhǔn)。上述配準(zhǔn)方法，對(duì)點(diǎn)云間位姿和噪聲程度有嚴(yán)格要求，同時(shí)容易陷入局部最優(yōu)。

針對(duì)上述點(diǎn)云數(shù)據(jù)在配準(zhǔn)中存在的問(wèn)題，本文采用散亂點(diǎn)云建立的曲率信息研究點(diǎn)云的配準(zhǔn)方法，在基于曲率特征配準(zhǔn)的基礎(chǔ)上，引入現(xiàn)在比較成熟的模擬退火法(simulated annealing，SAA)對(duì)點(diǎn)云數(shù)據(jù)拼接誤差的參數(shù)進(jìn)行尋優(yōu)，避免局部最優(yōu)。因此提出基于隨機(jī)抽樣一致性算法(random sampling consensus，RANSAC)提取曲率特征點(diǎn)，然后用距離最近原則尋找對(duì)應(yīng)點(diǎn)對(duì)，最后使用模擬退火法防止陷入局部最優(yōu)完成最終配準(zhǔn)。

1 經(jīng)典ICP算法

經(jīng)典ICP算法[14]利用最小二乘最優(yōu)匹配原理，通過(guò)迭代進(jìn)行點(diǎn)云之間的剛體變換。首先已知兩塊待配準(zhǔn)的重疊點(diǎn)云M和N，其中M為待配準(zhǔn)點(diǎn)云，N為參考點(diǎn)云。M點(diǎn)集中的每一點(diǎn)Mi以最小距離原則在N中尋找對(duì)應(yīng)點(diǎn)，對(duì)應(yīng)點(diǎn)構(gòu)成新的點(diǎn)集。所有對(duì)應(yīng)點(diǎn)對(duì)的拼接誤差[15]

(1)

式中：Mi為目標(biāo)點(diǎn)集；Ni為待配準(zhǔn)點(diǎn)集；n為匹配點(diǎn)對(duì)總數(shù)；λ、R、T分別表示點(diǎn)集M與點(diǎn)集N之間剛體變換的縮放參數(shù)、旋轉(zhuǎn)參數(shù)、平移參數(shù)。反復(fù)迭代求出滿足式(1)且收斂于給定閾值的轉(zhuǎn)換參數(shù)，即為ICP算法配準(zhǔn)的最終目的。ICP算法過(guò)程描述如下：

①尋找對(duì)應(yīng)點(diǎn)對(duì)M’、N’；

②計(jì)算配準(zhǔn)轉(zhuǎn)換參數(shù)λ、R、T；

③修正待配準(zhǔn)點(diǎn)云M得到新的點(diǎn)云M’；

④計(jì)算相鄰兩次迭代點(diǎn)集的距離差值di，點(diǎn)對(duì)匹配均方差收斂于閾值τ，即di-di+1<τ，則停止迭代，否則令i=i+1，并返回步驟①。

2 改進(jìn)ICP算法

2.1 曲率特征點(diǎn)集的提取

數(shù)據(jù)配準(zhǔn)時(shí)，對(duì)于待配準(zhǔn)點(diǎn)云M和基準(zhǔn)點(diǎn)云N，如果對(duì)所有點(diǎn)云進(jìn)行匹配點(diǎn)的搜索，會(huì)非常繁瑣且浪費(fèi)時(shí)間，同時(shí)因?yàn)樵肼暤拇嬖谝矔?huì)對(duì)特征點(diǎn)的選取產(chǎn)生影響。本文在進(jìn)行曲面擬合局部點(diǎn)云時(shí)，加入RANSAC算法，以消除噪點(diǎn)的影響，最后利用曲率極值法來(lái)提取曲率特征點(diǎn)。

在三維歐幾里得空間中，幾何體結(jié)構(gòu)的特征主要由曲率和法矢量來(lái)描述[16]。曲率是一種不隨平移、旋轉(zhuǎn)、縮放變化的特征，該特征可用高斯曲率KG、平均曲率KM、主曲率k進(jìn)行描述，首先引入曲率度量公式[17]

(2)

(3)

式中：k1代表最大曲率；k2表示最小曲率；α代表度量曲率相似度。通過(guò)公式(2)、公式(3)度量曲率相似度，提取出具有相似結(jié)構(gòu)的點(diǎn)云，然后配準(zhǔn)2片點(diǎn)云。特征點(diǎn)提取的具體步驟描述如下：

①選任意一點(diǎn)為種子點(diǎn)，選取K領(lǐng)域計(jì)算包含的點(diǎn)個(gè)數(shù)，計(jì)算迭代次數(shù)L[18]由式(4)得出：

(4)

(5)

式中：A為點(diǎn)云數(shù)；a為擬合曲面方程所需最小點(diǎn)個(gè)數(shù)，在本文中a=6，P為a個(gè)點(diǎn)均為有效點(diǎn)的概率；ε是數(shù)據(jù)的錯(cuò)誤率，表示局外點(diǎn)與總數(shù)據(jù)量之比。

②對(duì)點(diǎn)云任意一點(diǎn)，建立鄰域點(diǎn)集合，設(shè)曲面方程[19]為

(6)

在K鄰域內(nèi)隨機(jī)選取a個(gè)點(diǎn)，計(jì)算Qij的初始參數(shù)。

④重復(fù)步驟②～③，直到循環(huán)次數(shù)達(dá)L次，找到最大局內(nèi)點(diǎn)，由參數(shù)曲面方程解出最優(yōu)參數(shù)Qij。

(7)

(8)

⑥由求得的高斯曲率與平均曲率可直接求解主曲率

(9)

(10)

⑦判斷該點(diǎn)主曲率k1或k2在對(duì)應(yīng)的主方向上是否為極值，若為極值，則保留該點(diǎn)，進(jìn)行下一個(gè)種子點(diǎn)判斷。重復(fù)步驟⑦～⑧，直到遍歷完所有的數(shù)據(jù)。

⑧曲率極值越大特征越明顯。為曲率極值設(shè)定一個(gè)閾值ω，大于該設(shè)置閾值則保留曲率特征點(diǎn)，反之刪除。選取曲率特征點(diǎn)集流程圖如圖1所示。

圖1 選取曲率特征點(diǎn)集流程圖

2.2 獲取全局最優(yōu)轉(zhuǎn)換參數(shù)

配準(zhǔn)的過(guò)程即是求解待配準(zhǔn)點(diǎn)云到基準(zhǔn)點(diǎn)云的縮放參數(shù)λ、剛性旋轉(zhuǎn)R和平移T的過(guò)程。待配準(zhǔn)點(diǎn)云M中的A點(diǎn)和基準(zhǔn)點(diǎn)云N中的B點(diǎn)，表示的是地球表面物體的同一個(gè)點(diǎn)，可以通過(guò)縮放、旋轉(zhuǎn)和平移變換把A點(diǎn)和B點(diǎn)重合在一起。A點(diǎn)坐標(biāo)XA，YA，ZA，B點(diǎn)坐標(biāo)XB，YB，ZB，轉(zhuǎn)換的公式為：

(11)

在實(shí)際測(cè)量中，點(diǎn)云的采集不可避免帶入噪聲，雖然根據(jù)曲率特征選取特征點(diǎn)，以距離最近為判斷依據(jù)獲得相對(duì)應(yīng)的特征點(diǎn)對(duì)，但是點(diǎn)云間數(shù)據(jù)沒(méi)有絕對(duì)的一一對(duì)應(yīng)關(guān)系，因此引入模擬退火法，得到使目標(biāo)函數(shù)最小的最優(yōu)轉(zhuǎn)換參數(shù)λ、R、T，避免配準(zhǔn)時(shí)出現(xiàn)局部最優(yōu)解。目標(biāo)函數(shù)的形式多樣化，本文使用拼接誤差公式(1)作為目標(biāo)函數(shù)。

在1953年，Metropolis[20]等率先提出了模擬退火法。模擬退火法是一種非線性反演，能夠避免使反演陷入局部極大值，是一個(gè)不斷尋優(yōu)的過(guò)程。在此過(guò)程中，擬合度隨迭代次數(shù)的增加呈現(xiàn)一種跳躍起伏的現(xiàn)象，但總體的趨勢(shì)是變小或變大，然而正是由于模擬退火允許擬合誤差變大或擬合度變小，進(jìn)而使得模型從局部的最優(yōu)值中跳出來(lái)，最終得到全局的最優(yōu)化模型。因此引入模擬退火法，在點(diǎn)云拼接過(guò)程中以拼接誤差Eλ，R，T為目標(biāo)函數(shù)，進(jìn)行優(yōu)化求解。

在模擬退火法的迭代過(guò)程中，要求溫度隨著迭代次數(shù)的增加而緩慢降低，由冷卻進(jìn)度表控制溫度t的逐步下降，同時(shí)依概率接受惡化解。接受概率參數(shù)是關(guān)于溫度t的函數(shù)，隨著t不斷的降低，接受概率也不斷降低，直到最后不再接受任何惡化解。本文采用雙曲下降型[21]退火方案。模擬退火法尋找轉(zhuǎn)換參數(shù)的最優(yōu)具體解步驟：

①初始化目標(biāo)函數(shù)模型參數(shù)λ0、R0、T0，確定參數(shù)變化范圍，選擇初始溫度t0=1°。

②計(jì)算目標(biāo)函數(shù)值E0λ0，R0，T0。

③對(duì)當(dāng)前模型參數(shù)進(jìn)行擾動(dòng)，產(chǎn)生新的模型參數(shù)λ1、R1、T1。

④計(jì)算新模型參數(shù)下的目標(biāo)函數(shù)值E1(λ1，R1，T1)。

⑤計(jì)算目標(biāo)函數(shù)差ΔE，ΔE=E1-E0。

⑥判斷新模型參數(shù)是否接受，依Metropolis準(zhǔn)則[22]為判斷依據(jù)：若ΔE<0則接受；若ΔE>0，則新模型參數(shù)以概率P進(jìn)行接受。

P=exp-ΔE′/t

(12)

式中：t為溫度。

⑦接受新的模型參數(shù)時(shí)，設(shè)置λ0=λ1，R0=R1，T0=T1。

⑧在溫度t下重復(fù)一定次數(shù)的擾動(dòng)接受。

⑨降溫t

tG=t0αG1/a

(13)

式中：t0為初始溫度；G為迭代次數(shù)；α為給定常數(shù)，本文α=1。

⑩重復(fù)步驟③至⑨，直到收斂條件滿足。

2.3 改進(jìn)ICP算法流程

改進(jìn)配準(zhǔn)點(diǎn)云的ICP算法步驟如下：

①讀取需要配準(zhǔn)的2個(gè)點(diǎn)云：待配準(zhǔn)點(diǎn)云M，基準(zhǔn)點(diǎn)云N。

②利用RANSAC算法的二次曲面法對(duì)其K鄰域進(jìn)行曲面擬合。

③根據(jù)結(jié)合式(9)、公式(10)計(jì)算該點(diǎn)主曲率k1、k2，并判斷是否為曲率極值點(diǎn)，遍歷所有點(diǎn)云數(shù)據(jù)，并刪除小于閾值μ的曲率特征點(diǎn)，得到基準(zhǔn)點(diǎn)集M’。

④根據(jù)②、③計(jì)算待配準(zhǔn)點(diǎn)云數(shù)據(jù)的曲率特征點(diǎn)，得到待配準(zhǔn)點(diǎn)集N’。

⑤令為M’0、N’0迭代運(yùn)算初始點(diǎn)集

⑥以距離最近為依據(jù)，查詢特征點(diǎn)對(duì)(M’0、N’0)。

⑦根據(jù)步驟⑥獲得的對(duì)應(yīng)點(diǎn)集使用四元數(shù)法計(jì)算，獲取縮放參數(shù)λ、旋轉(zhuǎn)參數(shù)R、平移參數(shù)T。

⑧使用步驟⑥獲得的λ、R、T對(duì)點(diǎn)云M利用公式(11)進(jìn)行變換，得到新的點(diǎn)云。

⑨依模擬退火法判斷是否為全局最優(yōu)，既公式(1)全局最小，是則保留，否則重復(fù)⑤至⑧直到拼接誤差為全局最優(yōu)，最終完成配準(zhǔn)，算法流程圖如圖2所示。

圖2 算法流程圖

3 實(shí)驗(yàn)驗(yàn)證與結(jié)果分析

為了驗(yàn)證本文算法的正確性、有效性和實(shí)用性，進(jìn)行以下2組實(shí)驗(yàn)。實(shí)驗(yàn)數(shù)據(jù)采用斯坦福大學(xué)的Horse點(diǎn)云模型和Trimble TX8掃描儀采集某教學(xué)樓點(diǎn)云數(shù)據(jù)進(jìn)行數(shù)據(jù)配準(zhǔn)實(shí)驗(yàn)。所選用的2組數(shù)據(jù)都有十多種點(diǎn)云視角數(shù)據(jù)，在每組數(shù)據(jù)中隨機(jī)選取2種視角進(jìn)行配準(zhǔn)，并與經(jīng)典ICP算法進(jìn)行比較。本文算法是在經(jīng)典ICP算法基礎(chǔ)上進(jìn)行的改進(jìn)，算法中加入了噪聲點(diǎn)的剔除和避免出現(xiàn)局部最優(yōu)，使其相對(duì)于經(jīng)典ICP算法針對(duì)部分重疊點(diǎn)云的配準(zhǔn)問(wèn)題得到較好的結(jié)果。實(shí)驗(yàn)在Matlab 2010b軟件環(huán)境下編程實(shí)現(xiàn)。

3.1 實(shí)驗(yàn)一

圖3 Horse點(diǎn)云配準(zhǔn)結(jié)果

為驗(yàn)證本文方法的正確性，本文對(duì)所選用斯坦福大學(xué)的Horse點(diǎn)云數(shù)據(jù)進(jìn)行配準(zhǔn)，并將本文方法與經(jīng)典ICP算法配準(zhǔn)結(jié)果進(jìn)行比較。圖3為Horse點(diǎn)云配準(zhǔn)結(jié)果圖，紅色點(diǎn)云為基準(zhǔn)數(shù)據(jù)，綠色點(diǎn)云為待配準(zhǔn)數(shù)據(jù)，其中圖3(b)為經(jīng)典ICP算法配準(zhǔn)結(jié)果，圖3(c)為本文所提改進(jìn)算法配準(zhǔn)結(jié)果。

從圖3可以直觀地看出采用經(jīng)典ICP算法配準(zhǔn)結(jié)果較差，兩站點(diǎn)云不能很好地貼合在一起，雖然兩站點(diǎn)云有一部分可以配準(zhǔn)在一起，但是仍有大量點(diǎn)云數(shù)據(jù)配準(zhǔn)失敗；而本文所提算法配準(zhǔn)結(jié)果較好，經(jīng)計(jì)算對(duì)應(yīng)特征點(diǎn)對(duì)距離中誤差為2.7×10-3mm。雖然本文所提算法會(huì)增加配準(zhǔn)時(shí)間，但配準(zhǔn)結(jié)果可以體現(xiàn)出本文所改進(jìn)的算法可以很好地克服經(jīng)典ICP算法容易陷入局部最優(yōu)解的情況。

3.2 實(shí)驗(yàn)二

為驗(yàn)證本文算法的通用性和實(shí)用性，選取Trimble TX8掃描儀采集某教學(xué)樓的點(diǎn)云數(shù)據(jù)為實(shí)例進(jìn)行驗(yàn)證。所選用的某教學(xué)樓點(diǎn)云數(shù)據(jù)有十多站點(diǎn)云視角，在實(shí)驗(yàn)數(shù)據(jù)中隨機(jī)選取兩站數(shù)據(jù)進(jìn)行配準(zhǔn)，并與經(jīng)典ICP算法進(jìn)行比較。圖4為兩站點(diǎn)云配準(zhǔn)前后對(duì)比圖，紅色點(diǎn)云為待配準(zhǔn)點(diǎn)云M，綠色點(diǎn)云為基準(zhǔn)點(diǎn)云N。未配準(zhǔn)的兩站數(shù)據(jù)如圖4(a)所示，用本文算法配準(zhǔn)后點(diǎn)云數(shù)據(jù)如圖4(b)所示。

圖4 整體點(diǎn)云配準(zhǔn)前后對(duì)比圖

本文為驗(yàn)證該方法的穩(wěn)定性和避免局部最優(yōu)，采用2種算法對(duì)某教學(xué)樓的2個(gè)不同視角的點(diǎn)云進(jìn)行配準(zhǔn)。圖5為2種算法局部配準(zhǔn)結(jié)果圖，其中圖5(a)為ICP經(jīng)典算法的精確配準(zhǔn)圖，圖5(b)為為本文算法的精確配準(zhǔn)圖。

圖5 2種方法局部配準(zhǔn)結(jié)果

從圖5可以直觀地看出利用經(jīng)典ICP算法配準(zhǔn)后，待配準(zhǔn)數(shù)據(jù)和基準(zhǔn)點(diǎn)云數(shù)據(jù)之間有傾斜、錯(cuò)層現(xiàn)象。本文所提算法使兩站點(diǎn)云數(shù)據(jù)相互重合優(yōu)于經(jīng)典ICP算法。因此，利用本文算法可以更為精準(zhǔn)地進(jìn)行配準(zhǔn)，同時(shí)穩(wěn)健性更強(qiáng)。

為定量描述本文方法的配準(zhǔn)效果，在配準(zhǔn)后的點(diǎn)云數(shù)據(jù)中選取房屋角點(diǎn)作為特征點(diǎn)，計(jì)算配準(zhǔn)后對(duì)應(yīng)角點(diǎn)的距離，繪制距離誤差圖，如圖6所示。在獲取到的點(diǎn)云數(shù)據(jù)中，很難找到2個(gè)對(duì)應(yīng)的點(diǎn)數(shù)據(jù)表示真實(shí)物體的角點(diǎn)，本文利用平面擬合獲取待評(píng)價(jià)角點(diǎn)周?chē)?個(gè)平面信息，進(jìn)而通過(guò)3個(gè)面相交得到角點(diǎn)坐標(biāo)。

圖6 配準(zhǔn)后對(duì)應(yīng)角點(diǎn)距離殘差

從圖6可以清晰看出利用經(jīng)典ICP算法配準(zhǔn)后，大部分檢核點(diǎn)的配準(zhǔn)誤差在2 mm左右，但有部分檢核點(diǎn)配準(zhǔn)誤差大于5 mm，這是由于通過(guò)經(jīng)典ICP算法進(jìn)行配準(zhǔn)，出現(xiàn)局部最優(yōu)情況而產(chǎn)生影響。本文算法在提取的特征點(diǎn)時(shí)加入結(jié)合RANSAC算法，從而去除噪聲的影響提取特征點(diǎn)，同時(shí)利用模擬退火算法得到最優(yōu)的轉(zhuǎn)換參數(shù)，使配準(zhǔn)誤差優(yōu)于2.5 mm，且誤差分布較均勻。實(shí)驗(yàn)配準(zhǔn)圖和實(shí)驗(yàn)獲得的對(duì)比數(shù)據(jù)，充分表明本文所提配準(zhǔn)方法比經(jīng)典ICP算法更穩(wěn)定、精度更高。

4 結(jié)束語(yǔ)

在多視角測(cè)得散亂點(diǎn)云數(shù)據(jù)的情況下，利用經(jīng)典的ICP算法得到的配準(zhǔn)結(jié)果穩(wěn)定性不高，容易陷入局部最優(yōu)。針對(duì)這一問(wèn)題，本文在經(jīng)典ICP算法的基礎(chǔ)上，引入了一種模擬退火算法尋找全局最優(yōu)的配準(zhǔn)方法。該方法的核心思想是結(jié)合RANSAC算法剔除粗差點(diǎn)從而達(dá)到穩(wěn)健的效果，采用模擬退火算法避免點(diǎn)云配準(zhǔn)陷入局部最優(yōu)使得拼接誤差最小，得到全局最優(yōu)完成數(shù)據(jù)精確的配準(zhǔn)。實(shí)驗(yàn)結(jié)果表明，本文數(shù)據(jù)配準(zhǔn)算法提高了配準(zhǔn)的穩(wěn)定性和配準(zhǔn)精度，能夠滿足多視角點(diǎn)云數(shù)據(jù)的配準(zhǔn)要求，驗(yàn)證了本文提出算法的可行性。