徐衛(wèi)紅，彭家真，高紅旗

(浙江華東測(cè)繪與工程安全技術(shù)有限公司，310014，杭州)

0 引言

當(dāng)前遙感點(diǎn)云數(shù)據(jù)獲取最廣泛的手段包括激光掃描和攝影測(cè)量，這2種方式都能獲得高精度的點(diǎn)云數(shù)據(jù)，生成的點(diǎn)云已廣泛用于各種生產(chǎn)應(yīng)用[1]。LiDAR使用主動(dòng)式的探測(cè)方式，基于極坐標(biāo)幾何點(diǎn)位原理快速獲取空間點(diǎn)云信息，存在成本高、點(diǎn)云密度較低、缺乏地物光譜信息等缺點(diǎn)[2]。無(wú)人機(jī)攝影測(cè)量則是一種被動(dòng)測(cè)量目標(biāo)物的手段，其通過(guò)影像密集匹配可以得到比LiDAR點(diǎn)云密度更高、地物幾何特征更加明顯的空間點(diǎn)云。影像中豐富的光譜和紋理特征可以作為L(zhǎng)iDAR點(diǎn)云數(shù)據(jù)的補(bǔ)充，有效地解決了單一數(shù)據(jù)源信息不足問(wèn)題。通過(guò)LiDAR點(diǎn)云與影像融合，能產(chǎn)生更豐富的地物空間信息和語(yǔ)義信息，而2種點(diǎn)云數(shù)據(jù)融合的前提條件是進(jìn)行點(diǎn)云配準(zhǔn)。因此，對(duì)2種不同的點(diǎn)云數(shù)據(jù)進(jìn)行配準(zhǔn)是一個(gè)至關(guān)重要的環(huán)節(jié)。

點(diǎn)云配準(zhǔn)是指通過(guò)某種配準(zhǔn)算法計(jì)算得到2組不同點(diǎn)云的剛體變換參數(shù)，然后再利用變換參數(shù)將2組點(diǎn)云統(tǒng)一到同一坐標(biāo)系下。國(guó)內(nèi)外研究人員提出了許多點(diǎn)云配準(zhǔn)算法，其中ICP算法具有簡(jiǎn)單、復(fù)雜度低等優(yōu)點(diǎn)，是點(diǎn)云匹配最為常用的算法之一。然而，ICP算法存在計(jì)算效率低、魯棒性差、無(wú)初值時(shí)收斂慢等問(wèn)題。因此，針對(duì)ICP算法效率不高的問(wèn)題，許多研究人員提出了大量的改進(jìn)算法[3]。Li等通過(guò)動(dòng)態(tài)調(diào)整剛體變換參數(shù)的因子減少ICP算法迭代次數(shù)，實(shí)驗(yàn)結(jié)果表明改進(jìn)的算法提高了對(duì)應(yīng)點(diǎn)搜索效率，點(diǎn)云配準(zhǔn)時(shí)間相比于傳統(tǒng)方法要更短，且精度更高[4]。孫煒等針對(duì)ICP配準(zhǔn)時(shí)間長(zhǎng)問(wèn)題，提出了融合輪廓特征快速點(diǎn)云配準(zhǔn)方法[5]。李運(yùn)川等利用重心重合法提高點(diǎn)云重疊度，利用RANSAC算法得到點(diǎn)云配準(zhǔn)較好的變換參數(shù)初始值，最后改進(jìn)后的ICP算法實(shí)現(xiàn)點(diǎn)云高精度配準(zhǔn)，且配準(zhǔn)效率高于傳統(tǒng)的ICP算法[6]。這些改進(jìn)的算法主要是通過(guò)減少在2組點(diǎn)云中對(duì)應(yīng)點(diǎn)的搜索時(shí)間來(lái)提高配準(zhǔn)效率。

LiDAR點(diǎn)云與影像匹配點(diǎn)云的配準(zhǔn)方法主要分為3種。第1種是將影像匹配得到的點(diǎn)云與LiDAR點(diǎn)云進(jìn)行配準(zhǔn)。劉燕萍等基于影像密集匹配得到三維點(diǎn)云，通過(guò)粗配準(zhǔn)建筑物屋頂?shù)奶卣鹘屈c(diǎn)，再根據(jù)提取屋檐特征線以實(shí)現(xiàn)點(diǎn)云的精配準(zhǔn)[7]。李彩林等通過(guò)影像生成三維稀疏點(diǎn)云，以影像三維點(diǎn)云和激光點(diǎn)云擬合的最鄰近曲面為約束條件，結(jié)合共線條件方程實(shí)現(xiàn)影像三維點(diǎn)云和激光點(diǎn)云的高精度配準(zhǔn)[8]。Stamos等利用運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)(Structure from Motion，SfM)算法恢復(fù)物方稀疏點(diǎn)云，最后將LiDAR點(diǎn)云和攝影測(cè)量點(diǎn)云進(jìn)行配準(zhǔn)[9]。第2種是將LiDAR點(diǎn)云轉(zhuǎn)化為二維圖像，再利用點(diǎn)云配準(zhǔn)方法進(jìn)行匹配。趙中陽(yáng)等先利用LiDAR點(diǎn)云生成的深度圖像，基于航空影像和深度圖像的面特征采用SIFT算子提取點(diǎn)特征，最后通過(guò)面特征到點(diǎn)特征的配準(zhǔn)策略實(shí)現(xiàn)LiDAR點(diǎn)云與航空影像的配準(zhǔn)[10]。徐景中等通過(guò)LiDAR點(diǎn)云生成深度圖像，提取深度圖像與對(duì)應(yīng)影像的結(jié)構(gòu)特征，采用粗匹配-精匹配的循環(huán)迭代匹配策略，來(lái)實(shí)現(xiàn)LiDAR點(diǎn)云與航空影像的自動(dòng)配準(zhǔn)[11]。第3種是建立影像與LiDAR點(diǎn)云之間的直接配準(zhǔn)關(guān)系。Zhu等提出一種穩(wěn)健有效的由粗到細(xì)的配準(zhǔn)方法，利用空間約束和Gabor結(jié)構(gòu)特征進(jìn)行配準(zhǔn)，結(jié)果表明該配準(zhǔn)方法對(duì)幾何失真和輻射變化具有較高的魯棒性[12]。Habib等利用直線特征建立遙感影像與激光點(diǎn)云之間的直接映射關(guān)系，這種方法對(duì)于變化檢測(cè)和系統(tǒng)校準(zhǔn)應(yīng)用很有用[13]。王育堅(jiān)等提出一種基于保局(Principal Component Analysis，PCA)的三維點(diǎn)云配準(zhǔn)方法，利用特征尋找點(diǎn)云之間的對(duì)應(yīng)關(guān)系來(lái)實(shí)現(xiàn)配準(zhǔn)[14]。然而，大多LiDAR點(diǎn)云與影像配準(zhǔn)方法存在數(shù)據(jù)冗余、自動(dòng)匹配精度不高或效率低下等問(wèn)題。

本文基于ICP算法，采用粗-精配準(zhǔn)匹配策略來(lái)實(shí)現(xiàn)LiDAR點(diǎn)云與正射影像的配準(zhǔn)，使得點(diǎn)云與影像匹配兼顧精度與效率。首先利用體素濾波器對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理，再利用PCA算法壓縮數(shù)據(jù)量和獲得主分量進(jìn)行點(diǎn)云粗配準(zhǔn)來(lái)獲得點(diǎn)云剛體變換的初始值，最后采用ICP算法對(duì)點(diǎn)云進(jìn)行高效地精配準(zhǔn)。

1 LiDAR點(diǎn)云與正射影像配準(zhǔn)

1.1 點(diǎn)云預(yù)處理

本文通過(guò)無(wú)人機(jī)采集得到影像，利用Pix4D軟件經(jīng)過(guò)密集匹配得到點(diǎn)云數(shù)據(jù)，研究區(qū)域的正射影像如圖1所示。同時(shí)在激光雷達(dá)設(shè)備上輸出得到同一研究區(qū)域的激光點(diǎn)云，激光點(diǎn)云高度渲染圖見(jiàn)圖2。由正射影像和LiDAR渲染圖可以看出，正射影像含有豐富的紋理信息，而LiDAR點(diǎn)云具有高精度的空間幾何信息。

圖1 正射影像

圖2 LiDAR點(diǎn)云渲染圖

在進(jìn)行點(diǎn)云配準(zhǔn)之前，為了去除冗余點(diǎn)云數(shù)據(jù)，減少對(duì)應(yīng)點(diǎn)對(duì)搜索時(shí)間，本文利用體素濾波方法對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理。體素濾波在點(diǎn)云下采樣過(guò)程中的優(yōu)勢(shì)在于去除點(diǎn)云噪聲點(diǎn)和離群點(diǎn)的同時(shí)，保留了點(diǎn)云本身的幾何結(jié)構(gòu)。體素濾波的原理是根據(jù)輸入的點(diǎn)云計(jì)算得到包圍整個(gè)點(diǎn)云的立方體，將立方體根據(jù)設(shè)定的分辨率分割成多個(gè)小立方體，最后將小立方體內(nèi)點(diǎn)云的質(zhì)心坐標(biāo)近似為小立方體內(nèi)的點(diǎn)云。見(jiàn)圖3，體素濾波的計(jì)算過(guò)程如下：1)計(jì)算三維點(diǎn)云數(shù)據(jù)坐標(biāo)軸上的最大值和最小值；2)設(shè)計(jì)體素小柵格的邊長(zhǎng)；3)根據(jù)得到的坐標(biāo)軸最值計(jì)算點(diǎn)云最小包圍盒的邊長(zhǎng)；4)計(jì)算體素網(wǎng)格的尺寸；5)計(jì)算每個(gè)點(diǎn)云在體素小柵格內(nèi)的索引；6)從大到小排列索引里的元素，并以各體素重心代替其柵格內(nèi)的所有點(diǎn)，若重心不存在，則用距重心最近的數(shù)據(jù)點(diǎn)代替柵格內(nèi)所有點(diǎn)。

1.2 迭代最近點(diǎn)算法

ICP算法是被提出用于三維點(diǎn)云配準(zhǔn)[15]，該算法通過(guò)迭代搜索2個(gè)點(diǎn)云距離最近的對(duì)應(yīng)點(diǎn)，以各對(duì)應(yīng)點(diǎn)對(duì)歐氏距離平方和最小為目標(biāo)函數(shù)，再利用最小二乘計(jì)算最優(yōu)的2組點(diǎn)云配準(zhǔn)剛體變換參數(shù)。

假定有待配準(zhǔn)點(diǎn)云集P和基準(zhǔn)點(diǎn)云集M，M和P分別包含Nm和Np個(gè)點(diǎn)云。以m=[mR|mT]來(lái)表示2個(gè)點(diǎn)云集間的剛體變換參數(shù)。利用矩陣矢量m在基準(zhǔn)點(diǎn)云集中搜索對(duì)應(yīng)于待匹配點(diǎn)云集中最近的點(diǎn)，在這個(gè)過(guò)程當(dāng)中會(huì)搜尋到多個(gè)點(diǎn)。為了得到最優(yōu)變換參數(shù)，歐氏距離平方和最小的目標(biāo)函數(shù)表達(dá)式為：

圖3 體素濾波過(guò)程

(1)

式中mi和pi表示點(diǎn)云集M和P中的各個(gè)點(diǎn)云。

計(jì)算待匹配點(diǎn)云集和基準(zhǔn)點(diǎn)云集的重心：

(2)

(3)

計(jì)算2個(gè)點(diǎn)云集間的互協(xié)方差陣，并利用四元數(shù)法計(jì)算得到一個(gè)4×4的對(duì)稱矩陣：

(4)

(5)

再利用公式(6)計(jì)算平移向量T：

mT=μm-R(mR)μp

(6)

對(duì)點(diǎn)集P中的每一個(gè)點(diǎn)pi應(yīng)用變換函數(shù)E得到新的點(diǎn)云集Pnew，并在基準(zhǔn)點(diǎn)云集M中搜索Pnew的對(duì)應(yīng)點(diǎn)得到Mnew。

(7)

經(jīng)過(guò)以上過(guò)程的多次迭代計(jì)算，直到Pnew和Mnew的誤差小于給定的閾值或者達(dá)到設(shè)定的最大迭代次數(shù)時(shí)終止迭代，輸出R和T。

1.3 基于PCA的粗配準(zhǔn)

PCA是一種經(jīng)典的高維數(shù)據(jù)降維方法[16]。PCA算法在點(diǎn)云配準(zhǔn)過(guò)程中主要用于給精配準(zhǔn)階段提供剛體變換參數(shù)的初始值，提高點(diǎn)云配準(zhǔn)效率。其算法原理是通過(guò)計(jì)算點(diǎn)云的主方向，然后根據(jù)主方向計(jì)算得到點(diǎn)云配準(zhǔn)的初始值并對(duì)初始值進(jìn)行校正。首先計(jì)算2個(gè)點(diǎn)云集的協(xié)方差陣，通過(guò)協(xié)方差陣求出點(diǎn)云的主軸方向，然后利用主軸方向計(jì)算點(diǎn)云集的旋轉(zhuǎn)矩陣，并通過(guò)2個(gè)點(diǎn)云集的中心偏移量求出點(diǎn)云集的平移向量[17]?；赑CA方法進(jìn)行點(diǎn)云粗配準(zhǔn)的步驟如下：

計(jì)算2組點(diǎn)云的中心：

(8)

(9)

計(jì)算2個(gè)點(diǎn)云集間的協(xié)方差陣：

(10)

對(duì)CP和CM進(jìn)行奇異值分解，將協(xié)方差陣轉(zhuǎn)化為由2組點(diǎn)云的主方向UP和UM構(gòu)成的矩陣：

(11)

通過(guò)點(diǎn)云主方向計(jì)算旋轉(zhuǎn)矩陣和平移向量的初始值R0和T0：

(12)

最后，針對(duì)初始值R0和T0進(jìn)行校正，由于各個(gè)主軸方向存在反向的可能性，求出的初始R0和T0并不一定能用于點(diǎn)云精配準(zhǔn)。因此，還需要對(duì)R0和T0的初始值進(jìn)行校正，校正過(guò)程如下。

利用得到的R0和T0和公式(14)對(duì)點(diǎn)云進(jìn)行旋轉(zhuǎn)：

pnew=P*R0+T0

(13)

在基準(zhǔn)點(diǎn)云集M中搜索Pnew的最近距離對(duì)應(yīng)點(diǎn)，計(jì)算得到對(duì)應(yīng)點(diǎn)間的距離：

H=pnew-mi

(14)

計(jì)算初始配準(zhǔn)誤差：

erri=H*HT

(15)

(16)

式中erri為第i個(gè)對(duì)應(yīng)點(diǎn)間的誤差，error為2組點(diǎn)云間的平均均方誤差。

將UP中3個(gè)列向量(待配準(zhǔn)點(diǎn)云的主方向)與目標(biāo)點(diǎn)云各個(gè)軸的方向進(jìn)行對(duì)比，若方向相反，則把相反方向所對(duì)應(yīng)的UP的列向量取反，如公式(17)：

UP(i)=-UP(i)

(17)

1.4 ICP精配準(zhǔn)

得到預(yù)處理后的點(diǎn)云數(shù)后，利用迭代最近點(diǎn)算法對(duì)預(yù)處理后的點(diǎn)云進(jìn)行匹配，見(jiàn)圖4，配準(zhǔn)過(guò)程如下。

1)體素濾波精簡(jiǎn)點(diǎn)云。利用體素濾波對(duì)待匹配點(diǎn)云P和基準(zhǔn)點(diǎn)云M進(jìn)行采樣，計(jì)算獲得精簡(jiǎn)后的待匹配點(diǎn)云P′和基準(zhǔn)點(diǎn)云M′。

2)ICP迭代搜索最近對(duì)應(yīng)點(diǎn)。對(duì)于待匹配點(diǎn)云中的每一個(gè)點(diǎn)云pi，在基準(zhǔn)點(diǎn)云集中搜索其對(duì)應(yīng)點(diǎn)mi。

4)實(shí)施點(diǎn)云變換。對(duì)待匹配點(diǎn)云集中的點(diǎn)云pi進(jìn)行變換得到新的點(diǎn)云集P2。

5)計(jì)算變換后的點(diǎn)對(duì)距離。

6)迭代終止判斷。若通過(guò)變換參數(shù)計(jì)算得到的Pnew和Mnew的距離誤差小于給定的閾值或超過(guò)最大迭代次數(shù)時(shí)，終止迭代運(yùn)算。

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

為驗(yàn)證本文方法的有效性，本文實(shí)驗(yàn)選取具有代表性的四川某一地區(qū)的LiDAR點(diǎn)云與無(wú)人機(jī)數(shù)字正射影像進(jìn)行配準(zhǔn)實(shí)驗(yàn)。該測(cè)區(qū)包含豐富的地物信息如：耕地、林地、建筑物、道路、水塘等地表信息，覆蓋范圍約為1 km2。

LiDAR點(diǎn)云數(shù)據(jù)獲取采用SKY-LARK輕小型無(wú)人機(jī)LiDAR測(cè)量系統(tǒng)，該系統(tǒng)采用橢圓軌跡掃描方式，最高點(diǎn)頻能獲得高達(dá)60萬(wàn)點(diǎn)/s，最大測(cè)量距離可達(dá)500 m。采用大疆無(wú)人機(jī)航拍影像，飛行2個(gè)架次(每架次飛行時(shí)間約為18 min)，規(guī)劃航線共8條，相對(duì)航高設(shè)為150 m，航向與旁向重疊度均設(shè)為80%，地面分辨率約為2 cm，整個(gè)研究區(qū)的航拍影像為870幅，在Pix4D完成立體像對(duì)匹配、相對(duì)定向、空中三角測(cè)量等處理，并生成數(shù)字正射影像圖。本文方法LiDAR點(diǎn)云與正射影像配準(zhǔn)后，點(diǎn)云三維渲染圖見(jiàn)圖5，紋理映射后結(jié)果圖見(jiàn)圖6。

圖5 LiDAR三維渲染圖

圖6 附帶紋理的三維圖

此外，為對(duì)比分析本文方法針對(duì)原始的ICP算法的改進(jìn)效果，在同等軟硬件環(huán)境下，采用相同的數(shù)據(jù)對(duì)ICP算法進(jìn)行LiDAR點(diǎn)云與無(wú)人機(jī)正射影像配準(zhǔn)，并采集了10個(gè)地面控制點(diǎn)(Ground Control Point, GCP)進(jìn)行精度檢驗(yàn)見(jiàn)圖7，表1為本文方法與原始ICP算法在迭代運(yùn)算次數(shù)、平均誤差、效率方面的對(duì)比。

圖7 LiDAR三維渲染圖

表1 ICP算法與本文方法對(duì)比

由對(duì)比結(jié)果可知，本文提出的方法配準(zhǔn)結(jié)果明顯優(yōu)于ICP算法，能實(shí)現(xiàn)LiDAR點(diǎn)云與正射影像高精度配準(zhǔn)。在配準(zhǔn)效率方面，采取基于粗-精的配準(zhǔn)策略，能去除點(diǎn)云噪聲，并加快配準(zhǔn)時(shí)間，本文采用體素濾波方法對(duì)原始點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理，并且利用PCA算法進(jìn)行粗配準(zhǔn)，以獲取精配準(zhǔn)過(guò)程中所需的剛體變換參數(shù)的初始值，相比于未經(jīng)過(guò)預(yù)處理和粗配準(zhǔn)的點(diǎn)云配準(zhǔn)，處理后的點(diǎn)云配準(zhǔn)速度有明顯的加快，這有利于進(jìn)行點(diǎn)云實(shí)時(shí)配準(zhǔn)。在精度方面，本文通過(guò)PCA粗配準(zhǔn)和ICP精配準(zhǔn)的點(diǎn)云配準(zhǔn)策略，實(shí)現(xiàn)了點(diǎn)云與正射影像高精度配準(zhǔn)。相比原始ICP算法精度和效率分別提高了約1.8倍和3倍。

3 結(jié)論

針對(duì)傳統(tǒng)的ICP算法難以滿足LiDAR點(diǎn)云與正射影像高精度、高效的配準(zhǔn)的需求，本文結(jié)合PCA粗配準(zhǔn)與ICP精配準(zhǔn)策略來(lái)實(shí)現(xiàn)LiDAR點(diǎn)云與無(wú)人機(jī)正射影像配準(zhǔn)。為去除LiDAR點(diǎn)云數(shù)據(jù)中的離群點(diǎn)和噪聲點(diǎn)，本文采用體素濾波方法對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理以消除冗余數(shù)據(jù)，提高點(diǎn)云與影像配準(zhǔn)效率，然后，利用PCA算法計(jì)算點(diǎn)云配準(zhǔn)剛體變換參數(shù)初值，并校正后作為ICP算法的輸入值，從而提高ICP算法配準(zhǔn)的準(zhǔn)確性和收斂效率。實(shí)驗(yàn)結(jié)果表明，本文的配準(zhǔn)策略在點(diǎn)云與正射影像配準(zhǔn)方面能獲得了較好的配準(zhǔn)結(jié)果，能滿足包含復(fù)雜地表地物分布的LiDAR點(diǎn)云與影像多源數(shù)據(jù)配準(zhǔn)需求。