基于ICP算法的非常規(guī)航空影像與機(jī)載LiDAR數(shù)據(jù)配準(zhǔn)

杜全葉，麥曉明，褚福俠，彭向陽(yáng)，王銳

(1.中國(guó)測(cè)繪科學(xué)研究院，北京 100830；2.廣東電網(wǎng)公司電力科學(xué)研究院，廣州 510080；3.濟(jì)寧光大房地產(chǎn)服務(wù)有限責(zé)任公司，濟(jì)寧 272000)

1 引 言

隨著機(jī)載激光掃描(LiDAR)技術(shù)推廣使用，短時(shí)間內(nèi)可以獲取大范圍地面離散點(diǎn)的三維坐標(biāo)。航空影像在地表紋理表達(dá)上有著不言而喻的明顯優(yōu)勢(shì)。航空影像與機(jī)載LiDAR數(shù)據(jù)配準(zhǔn)，有利于發(fā)揮兩種數(shù)據(jù)的優(yōu)勢(shì)，為建模和紋理映射等應(yīng)用，打下良好的基礎(chǔ)。

如何配準(zhǔn)激光掃描數(shù)據(jù)和對(duì)應(yīng)的光學(xué)影像，很多學(xué)者都做過(guò)一些研究。Habbib使用LiDAR和影像數(shù)據(jù)上人字形房屋頂交線作為配準(zhǔn)的基元，將這些直線特征的點(diǎn)坐標(biāo)代入變換函數(shù)，利用最小二乘原理解算出配準(zhǔn)參數(shù)[1]。吳波利用影像的平面特征和LiDAR數(shù)據(jù)中的離散點(diǎn)作為配準(zhǔn)的基元[2]，配準(zhǔn)航空影像和機(jī)載LiDAR數(shù)據(jù)。Mwafag Ghanma在其博士論文中，討論了基于直線段和面片的航空影像與LiDAR數(shù)據(jù)的配準(zhǔn)[3]。這些配準(zhǔn)方法局限于有直線特征或平面特征的區(qū)域。

航空攝影測(cè)量一般需按航空攝影規(guī)范進(jìn)行航線的設(shè)計(jì)，然后按預(yù)設(shè)的航線進(jìn)行飛行和影像獲取。但是在發(fā)生地震、滑坡等自然災(zāi)害的應(yīng)急響應(yīng)中，為了第一時(shí)間獲取通往災(zāi)區(qū)的主要道路、橋梁等交通設(shè)施的毀壞情況，及時(shí)地為救援行動(dòng)提供決策參考，就必須進(jìn)行“非常規(guī)”的航空攝影，即沿著主要的道路進(jìn)行飛行[4]。使用機(jī)載激光掃描設(shè)備和數(shù)碼成像設(shè)備同時(shí)進(jìn)行常規(guī)作業(yè)時(shí)，也可能會(huì)因相機(jī)像幅不大，視場(chǎng)角相對(duì)激光掃描儀的視場(chǎng)角小，而造成影像重疊度達(dá)不到53%。這些非常規(guī)航空影像無(wú)法進(jìn)行傳統(tǒng)的區(qū)域網(wǎng)空中三角測(cè)量，它們與機(jī)載LiDAR數(shù)據(jù)的配準(zhǔn)更加困難。

2 原理和方法

利用攝影測(cè)量中的共線方程[5]描述配準(zhǔn)后的LiDAR點(diǎn)和影像的對(duì)應(yīng)關(guān)系，即配準(zhǔn)后LiDAR腳點(diǎn)，影像上對(duì)應(yīng)的像點(diǎn)，以及攝影中心在同一條直線上。此時(shí)，配準(zhǔn)航空影像與LiDAR數(shù)據(jù)，即轉(zhuǎn)化為求解影像在LiDAR坐標(biāo)系下方位元素的問(wèn)題。首先，航空影像上提取密集特征點(diǎn)，在LiDAR點(diǎn)云和定位定向系統(tǒng)(Position and Orientation System，POS)數(shù)據(jù)輔助下，利用匹配和前方交會(huì)得到匹配點(diǎn)云，然后配準(zhǔn)匹配點(diǎn)云與LiDAR點(diǎn)云，最后，篩選良好的配準(zhǔn)點(diǎn)解算出影像方位元素，實(shí)現(xiàn)影像與LiDAR數(shù)據(jù)配準(zhǔn)，流程如圖1所示。

圖1 航空影像與LiDAR數(shù)據(jù)配準(zhǔn)流程

航空影像需要有全球定位系統(tǒng)(Global Positioning System，GPS)和IMU測(cè)量的POS數(shù)據(jù)作為影像的初始外方位元素，若沒(méi)有POS數(shù)據(jù)，則需要人工選取影像與LiDAR數(shù)據(jù)的同名點(diǎn)，解算影像初始方位元素。

2.1 影像匹配

采用立體像對(duì)密集匹配方式獲得影像匹配點(diǎn)云。如圖2所示，航攝區(qū)域連續(xù)影像，前4張標(biāo)記為Photo1、Photo2、Photo3和Photo4，斜線區(qū)域?yàn)橄噜徲跋裰丿B區(qū)。

圖2 立體像對(duì)密集匹配示意圖

首先，在每張影像上提取密集Harris特征點(diǎn)。然后，將Photo1與Photo2在LiDAR點(diǎn)云生成的數(shù)字表面模型(Digital Surface Model，DSM)和POS數(shù)據(jù)輔助下進(jìn)行特征點(diǎn)匹配，得到兩張影像匹配點(diǎn)云——“Photo1-Photo2匹配點(diǎn)云”，將Photo2與Photo3進(jìn)行特征點(diǎn)匹配，得到“Photo2-Photo3匹配點(diǎn)云”，以此類推。

2.2 點(diǎn)云配準(zhǔn)

計(jì)算機(jī)視覺(jué)界已經(jīng)對(duì)三維離散點(diǎn)云的拼接做了大量的研究工作，最常用的點(diǎn)云拼接方法是1992年由Besl教授提出的迭代最臨近點(diǎn)算法(ICP，Iterative Closest Point)[6]。它是一個(gè)尋找兩個(gè)三維表面點(diǎn)集，最優(yōu)匹配幾何變換的迭代優(yōu)化過(guò)程。目的是確定相鄰兩個(gè)坐標(biāo)系之間的轉(zhuǎn)換參數(shù)F，包括旋轉(zhuǎn)矩陣R和平移向量T，以使得位于坐標(biāo)系1中的點(diǎn)P與位于坐標(biāo)系2中的點(diǎn)Q重合，即FP=Q，表達(dá)成R和T函數(shù)關(guān)系式即為RP+T=Q。P與Q不可能完全重合，因此，該問(wèn)題就轉(zhuǎn)化為尋找參數(shù)F，使得|FP-Q|*|FP-Q|的值最小。具體實(shí)現(xiàn)是通過(guò)迭代的方法找到均方差最小的變換。具體步驟如下：

①k為第k次迭代，旋轉(zhuǎn)參數(shù)Rk和平移參數(shù)Tk，Pk和Qk分別代表參與計(jì)算的兩個(gè)點(diǎn)集，那么兩個(gè)點(diǎn)集的函數(shù)關(guān)系式為RkPk+Tk=Qk，初始化k=0。

②尋找Pk中對(duì)應(yīng)的每個(gè)最近點(diǎn)Qk(k為第k次迭代)；

③尋找互換最鄰近點(diǎn)Pεk和Qεk(同時(shí)互為最近點(diǎn)且距離小于ε時(shí)才被標(biāo)注)；

④解算Pεk和Qεk間的均方距離dk(Pεk和Qεk是第k次迭代中互換的最鄰近點(diǎn))；

⑤解算Pε0和Qεk間最小二乘意義下的三維變換參數(shù)Rk、Tk；

⑥執(zhí)行變換：Qk+1=RkPεk+Tk；

目前，地面激光掃描點(diǎn)云拼接，機(jī)載LiDAR數(shù)據(jù)的拼接和星載LiDAR數(shù)據(jù)的拼接都應(yīng)用此算法。本文也將ICP算法應(yīng)用于密集匹配點(diǎn)云與LiDAR點(diǎn)云的配準(zhǔn)。ICP算法中所要求的初值，由影像初始方位元素確定。

2.3 單像空間后方交會(huì)

單像空間后方交會(huì)是利用控制點(diǎn)及其相應(yīng)的像點(diǎn)求解航片外方位元素。本文篩選點(diǎn)云配準(zhǔn)后距離小于平均LiDAR點(diǎn)間隔的配準(zhǔn)點(diǎn)作為控制，采用基于單位四元數(shù)的方法進(jìn)行單像空間后方交會(huì)[7]。基本思路為：先求出攝站點(diǎn)與控制點(diǎn)的距離；然后利用單位四元數(shù)表示旋轉(zhuǎn)，從而獲得旋轉(zhuǎn)矩陣，得到外方位元素；最后再求出攝站點(diǎn)的坐標(biāo)。這種方法不需要對(duì)共線條件方程進(jìn)行線性化，且無(wú)需提供外方位元素初值，計(jì)算簡(jiǎn)便，對(duì)于大傾角攝影的航片也能適用。

圖2中重疊區(qū)域的匹配點(diǎn)云分別與對(duì)應(yīng)的LiDAR數(shù)據(jù)配準(zhǔn)后，篩選良好的配準(zhǔn)點(diǎn)，利用單像空間后方交會(huì)解算出影像Photo2、Photo3等的方位元素，實(shí)現(xiàn)航空影像與LiDAR數(shù)據(jù)的配準(zhǔn)。Photo1影像上僅右邊有少部分配準(zhǔn)點(diǎn)，必要時(shí)也可后方交會(huì)解算其方位元素。

3 配準(zhǔn)實(shí)驗(yàn)一

3.1 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)數(shù)據(jù)為使用Leica ALS50-II獲取的某平原區(qū)LiDAR數(shù)據(jù)和隨帶RCD105數(shù)碼相機(jī)拍攝的航空影像，航高為500m，影像重疊度約為20%。影像像幅為7162×5389像素，像元大小為0.0068mm，影像地面分辨率為0.1m，平均LiDAR點(diǎn)間隔為0.4m。以連續(xù)4張影像和對(duì)應(yīng)的LiDAR數(shù)據(jù)為實(shí)驗(yàn)對(duì)象。

3.2 人工配準(zhǔn)

從影像和對(duì)應(yīng)的LiDAR點(diǎn)云數(shù)據(jù)中手工選取9個(gè)同名點(diǎn)，獲得同名點(diǎn)像坐標(biāo)和物方坐標(biāo)，經(jīng)過(guò)單像空間后方交會(huì)獲得影像外方位元素。第二張和第三張影像手工配準(zhǔn)過(guò)程中，后方交會(huì)中誤差(Root Mean Square Error，RMSE)如表1所示。

3.3 影像密集匹配

自動(dòng)配準(zhǔn)過(guò)程中，相鄰兩張影像密集匹配獲得的匹配點(diǎn)云如圖3所示。

表1 手工配準(zhǔn)后方交會(huì)中誤差

圖3 影像立體匹配點(diǎn)云

密集匹配點(diǎn)使用影像初始方位元素前方交會(huì)出點(diǎn)的物方三維坐標(biāo)。相鄰兩張影像重疊度不大，一個(gè)像對(duì)僅匹配出單張影像上20%左右的密集點(diǎn)云。

3.4 點(diǎn)云配準(zhǔn)結(jié)果

影像匹配點(diǎn)云與LiDAR點(diǎn)云配準(zhǔn)中誤差如表2所示。影像匹配點(diǎn)云與LiDAR點(diǎn)云配準(zhǔn)后，Photo2和Photo3單像空間后方交會(huì)中誤差如表3所示。

表2 點(diǎn)云配準(zhǔn)中誤差

表3 后方交會(huì)中誤差

3.5 配準(zhǔn)后點(diǎn)云投影結(jié)果

配準(zhǔn)之后，選擇LiDAR點(diǎn)云數(shù)據(jù)中的一些平頂房房頂點(diǎn)云，投影到對(duì)應(yīng)影像上檢查配準(zhǔn)效果。其中一個(gè)房頂點(diǎn)云數(shù)據(jù)如圖4中紅色區(qū)域所示。

圖4 LiDAR點(diǎn)云中的屋頂數(shù)據(jù)

第2張影像和第3張影像中，配準(zhǔn)前后其中一個(gè)屋頂數(shù)據(jù)的投影結(jié)果如圖5、圖6、圖7和圖8所示。

圖5 Photo2影像配準(zhǔn)前點(diǎn)云投影結(jié)果

圖6 Photo2影像配準(zhǔn)后點(diǎn)云投影結(jié)果

圖7 Photo3影像上配準(zhǔn)前點(diǎn)云投影結(jié)果

圖8 Photo3影像上配準(zhǔn)后點(diǎn)云投影結(jié)果

其中，紅色點(diǎn)為L(zhǎng)iDAR腳點(diǎn)投影到影像上的位置。從圖中可以明顯看出，在配準(zhǔn)后房屋的錯(cuò)位現(xiàn)象已經(jīng)消除，航空影像與LiDAR數(shù)據(jù)已經(jīng)配準(zhǔn)好。

4 配準(zhǔn)實(shí)驗(yàn)二

4.1 實(shí)驗(yàn)數(shù)據(jù)

如同實(shí)驗(yàn)一，只是相鄰影像重疊度增大，約為60%。

4.2 密集匹配結(jié)果

連續(xù)3張影像進(jìn)行立體像對(duì)密集匹配，經(jīng)過(guò)前方交會(huì)生成的點(diǎn)云如圖9所示。

圖9 影像立體匹配點(diǎn)云

4.3 點(diǎn)云配準(zhǔn)結(jié)果

密集匹配點(diǎn)云與LiDAR點(diǎn)云配準(zhǔn)中誤差如表4所示。影像匹配點(diǎn)云與LiDAR點(diǎn)云配準(zhǔn)后，Photo2影像后方交會(huì)中誤差為0.005114mm。

表4 點(diǎn)云配準(zhǔn)中誤差

5 實(shí)驗(yàn)分析

從Photo2影像對(duì)應(yīng)的LiDAR數(shù)據(jù)中選取房屋屋頂數(shù)據(jù)，按照配準(zhǔn)前方位元素、手工配準(zhǔn)后方位元素和通過(guò)點(diǎn)云配準(zhǔn)后的方位元素，投影到影像上，如圖10、圖11和圖12所示。

通過(guò)對(duì)比可以看出，配準(zhǔn)前LiDAR點(diǎn)云投影到影像上有較大錯(cuò)位，人工配準(zhǔn)的方式，已經(jīng)消除了錯(cuò)位現(xiàn)象，LiDAR點(diǎn)投影到影像上后誤差在一個(gè)LiDAR點(diǎn)間距以內(nèi)，而通過(guò)匹配點(diǎn)云配準(zhǔn)后屋頂?shù)腖iDAR點(diǎn)投影到影像后誤差更小，航空影像與LiDAR數(shù)據(jù)配準(zhǔn)的精度更高。

圖10 配準(zhǔn)前屋頂LiDAR點(diǎn)投影結(jié)果

圖11 人工配準(zhǔn)后屋頂LiDAR點(diǎn)投影結(jié)果

圖12 自動(dòng)配準(zhǔn)后屋頂LiDAR點(diǎn)投影結(jié)果

通過(guò)表1手工配準(zhǔn)后方交會(huì)中誤差和表 3 后方交會(huì)中誤差(點(diǎn)云配準(zhǔn)方式)對(duì)比，可以看出，采用手工配準(zhǔn)的方式一般情況下可以達(dá)到影像的4個(gè)像素，也就是一個(gè)LiDAR點(diǎn)間隔的精度。而采用影像匹配點(diǎn)云配準(zhǔn)的方式即便只有20%的重疊度的情況下也能達(dá)到影像4個(gè)像素以內(nèi)，良好情況下達(dá)到2個(gè)像素的精度，即半個(gè)LiDAR點(diǎn)間隔的精度。如果重疊度達(dá)到60%，那么點(diǎn)云配準(zhǔn)的精度會(huì)更高，解算的影像方位元素精度可以到達(dá)1個(gè)像素以內(nèi)。

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

通過(guò)影像匹配點(diǎn)云和LiDAR點(diǎn)云配準(zhǔn)，然后利用配準(zhǔn)點(diǎn)后方交會(huì)解算影像方位元素的配準(zhǔn)方法，相對(duì)于手工配準(zhǔn)航空影像和機(jī)載LiDAR數(shù)據(jù)，自動(dòng)化程度大大提高，并且達(dá)到的精度也更高。這種方法對(duì)于解決非常規(guī)攝影測(cè)量的航空影像與LiDAR數(shù)據(jù)的自動(dòng)高精度配準(zhǔn)非常有效，也適用于重疊度更大的常規(guī)航空影像。

參考文獻(xiàn)：

[1] HABIB A F，SHIN S，KIM C，et al.Integration of Photogrammetric and LIDAR Data in a Multi-Primitive Triangulation Environment[J].Innovations in 3D Geo Information Systems，2006：29-45.

[2] 吳波.自適應(yīng)三角形約束下的立體影像可靠匹配方法[D].武漢：武漢大學(xué)博士學(xué)位論文，2006.

[3] GHANMA M.Integration of photogrammetry and LIDAR[D].Calgary：Universityof Calgary (Canada).Ph.D.dissertation，2006.

[4] 張祖勛，郭大海，柯濤，等.抗震救災(zāi)中航空攝影測(cè)量的應(yīng)急響應(yīng)[J].遙感學(xué)報(bào)，2008，12(6)：852-857.

[5] 王之卓.攝影測(cè)量原理[M].武漢：武漢大學(xué)出版社，2007.20-21.

[6] BESL P J，MCKAY N D.A method for registration of 3-D shapes[J]，IEEE Transactions on Pattern Analysis and Machine Intelligence，1992，14(2)：239-256.

[7] 官云蘭，程效軍，周世健，等.基于單位四元數(shù)的空間后方交會(huì)解算[J].測(cè)繪學(xué)報(bào)，2008，37(1)：30-35.