一種三維地形特征提取和匹配方法

田 陽，李國慶,2，宋 新

(1. 哈爾濱工業(yè)大學深空探測基礎(chǔ)研究中心，哈爾濱 150080；2. 北京空間飛行器總體設(shè)計部，北京 100094；3. 哈爾濱工業(yè)大學電子與信息工程學院，哈爾濱 150080)

0 引 言

利用光學設(shè)備感知環(huán)境是航天器尤其是行星漫游器和著陸器等工作在陌生環(huán)境中的航天器必備能力之一[1]。例如NASA的Spirit，Opportunity和Curiosity火星漫游器，以及我國的玉兔月球漫游器[2]都裝備了雙目相機獲取周圍的三維地形信息，NASA最新的OSIRIS-REx小天體探測任務(wù)中，著陸器搭載了工作距離超過7 km的LiDAR敏感器，在著陸階段測量著陸器高度和天體表面地形[3]。除了利用三維地形數(shù)據(jù)感知巖石、隕石坑和陡坡等對航天器構(gòu)成危險的地形和障礙之外，三維數(shù)據(jù)還可以用來進行航天器導航定位。利用連續(xù)的三維地形數(shù)據(jù)序列能夠確定漫游器或是著陸器的幀間相對運動，這與Spirit和Opportunity漫游器采用的視覺里程計類似。另一種方式是將獲得的三維地形數(shù)據(jù)與導航地形數(shù)據(jù)庫對比，確定航天器在導航坐標系下的位置和姿態(tài)。不過這都建立在準確匹配三維地形數(shù)據(jù)的基礎(chǔ)上。

LiDAR或雙目相機獲取的三維地形數(shù)據(jù)一般是三維點云，匹配三維點云可以采用迭代最近點(Iterative Closest Point，ICP)方法[4]和該方法的多種變形[5-6]。ICP其實是一種非線性最小二乘方法，通過迭代旋轉(zhuǎn)和平移參數(shù)不斷縮小點云間的空間距離，在三維點云匹配中得到廣泛應(yīng)用，但是該方法需要猜測合適的迭代初值，效果較好的ICP變形方法將考察點到點距離改為點到平面距離，這需要對點云數(shù)據(jù)逐點計算平面法線方向，計算量較大。另一類三維點云匹配方法與二維圖像特征匹配方法類似，首先在點云數(shù)據(jù)中檢測形狀或曲率具有顯著性的點作為特征，利用特征鄰域內(nèi)的點云形成描述矢量(描述符)，通過計算描述矢量的相似性判斷特征匹配與否。

出于任務(wù)安全性的考慮，行星漫游器或著陸器工作環(huán)境都經(jīng)過精心設(shè)計，如行星著陸任務(wù)的著陸點一般都選擇地形平坦的區(qū)域，行星漫游器在路徑規(guī)劃過程中也盡量避免通過崎嶇的地形，這使得三維地形特征提取變得十分困難，更為重要求的是這些提取出的特征缺乏顯著性，容易出現(xiàn)誤匹配問題。

本文以行星漫游器和著陸器地形相關(guān)導航為應(yīng)用背景，提出一套三維特征提取和匹配方法，解決平坦地形條件下三維特征魯棒匹配問題。首先在第1節(jié)中介紹基于點云鄰域散布矩陣特征值的三維特征提取方法，接著將在第2節(jié)中構(gòu)建基于特征間有向距離的特征描述方法，最后通過仿真實驗校驗所提方法的有效性。本文提出的三維特征提取與匹配方法不但可以應(yīng)用在深空天體巡視采樣和著陸導航定位中，還在物體識別以及SLAM(Simultaneous Localization and Mapping)中有重要作用。

1 地形特征檢測

目前，三維特征檢測方法主要有兩類，一類是固定尺度算法，這類算法在固定的尺度下通過基于顯著性的非極大值抑制(Saliency-based NMS)方法檢測三維特征，如LSP(Local Surface Patches)方法[7]，ISS(Intrinsic Shape Signatures)方法[8]，KPQ(KeyPoint Quality)方法[9]和Heat Kernel Signature(HKS)方法[10]等；另一類是自適應(yīng)尺度算法，與SIFT算法思路類似，通過生成尺度空間并選取不同的尺度自適應(yīng)檢測三維特征，包括Laplace-Beltrami Scale-space (LBSS)[11]，MeshDoG[12]，KeyPoint Quality-Adaptive-Scale (KPQ-AS)[13]，Salient Points (SP)[14]方法等。此外文獻[15]針對月球表面地形特征提取問題，提出采用高程數(shù)據(jù)方差檢測三維特征。從特征重復檢出率等性能指標來看，ISS方法是較為理想的檢測方法[16]。

1.1 ISS三維特征檢測方法

ISS特征檢測方法的核心思想是通過特征點鄰域內(nèi)的點云位置協(xié)方差刻畫特征點的顯著性。設(shè)待考察的三維點坐標為p，q為點p鄰域N(p)內(nèi)的三維點坐標，那么點p對應(yīng)的散布矩陣表示為

(1)

式中：M為鄰域內(nèi)點的個數(shù)。若用λ1，λ2和λ3分別表示散布矩陣從大到小的三個特征值，那么滿足如下條件的三維點將被選取為特征點

(2)

式中：τ12和τ23是設(shè)定的檢測閾值?？梢钥闯?，ISS方法的特征選取策略是從特征的顯著性角度出發(fā)的，要求特征鄰域內(nèi)的點分布在半長軸有明顯區(qū)別的散布橢球內(nèi)，這使不同特征對應(yīng)的描述矢量具有一定的差異，同時也使得以散布矩陣特征向量為坐標軸的特征局部坐標系更加穩(wěn)定。在結(jié)構(gòu)化的場景中，這種特征選取策略十分有效，但在火星、月球等深空天體表面，主要分布著隕石坑、巖石等地形特征，這些特征的主要特點是突出或凹進地表，那么根據(jù)式(2)給出的特征點選取標準，一些有效的特征地形將被剔除。因此，本文將針對深空天體表面三維地形的特點在下一節(jié)中改進ISS檢測算法。

1.2 基于散布矩陣的起伏地形特征檢測方法

考慮到深空天體表面地形的特點，將ISS特征選取的策略修改為

(3)

此時，只要保證特征鄰域內(nèi)的點具有一定的起伏，即可被選取為地形特征。同時，為了去除在一處起伏地形附近檢出過多的特征點，利用散布矩陣特征值的乘積作為非極大值抑制的手段，也就是在特征點鄰域內(nèi)，只選取特征值乘積最大的點作為地形特征。

激光雷達或是雙目相機采集三維地形，是在離散空間內(nèi)進行的，運動中多次測量可能造成特征點位置在采樣間隔內(nèi)發(fā)生改變，因此這里用鄰域內(nèi)三維點位置的均值μp表示特征位置，以提高觀測過程中特征位置的穩(wěn)定性。

事實上，修改后的特征選取策略放寬了對特征鄰域內(nèi)三維點分布的要求，可以檢測出更多的地形特征，但是檢出特征的顯著性也隨之降低，這對于基于特征鄰域內(nèi)三維點分布或法線信息形成的特征描述方法是不可接受的，本文將在下一節(jié)建立一種基于有向距離的特征描述方法，由于采用特征間相對位置關(guān)系描述特征，因此可以有效的區(qū)分相似特征。

2 地形特征描述

已有的三維特征描述方法大多以特征附近點云空間分布或直方圖等統(tǒng)計手段作為描述，如Johnson提出的Spin Image(SI)方法[17]，F(xiàn)rome等提出的3D Shape Context(3DSC)方法[18]，以及與3DSC相類似的SHOT(Signature of Histograms of Orientations)[19]，Radu提出的Point Feature Histogram(PFH)[20]方法等。這些方法在描述密集的三維點云時都有良好的效果，但在解決稀疏點云匹配問題時，正確匹配率不理想[21]。其原因是稀疏的三維點不足以穩(wěn)定的描述鄰域形狀，采樣位置發(fā)生輕微改變就會導致特征描述符較大變化，這一現(xiàn)象在漫游器工作的平坦地形環(huán)境中更為嚴重。

2.1 三維特征間距離描述集合

為了解決稀疏點云不易匹配的問題，利用每次觀測檢出特征之間的距離以及特征相對位置矢量與特征主方向夾角(以下簡稱為方向夾角)作為基本元素，構(gòu)成特征描述集合，如圖1所示。這里之所以采用距離與方向夾角的集合描述特征，是由于地形遮擋以及特征檢測算法重復檢出率的客觀限制，不能保證地形特征每次都被檢出，特征數(shù)量在觀測過程中是不斷變化的，形成的特征描述矢量長度也不斷變化，因此這里采用相對距離與方向夾角的集合作為特征描述。

事實上，如果能夠保證較高的特征重復檢出率，前后兩次觀測中大量公共特征將被重復檢出，僅利用特征間的相對距離就可以很好的描述地形特征。而當重復檢出率較低時(根據(jù)文獻[16]中的測試結(jié)果，ISS特征檢測算法對激光雷達數(shù)據(jù)的重復檢出率僅有40%左右)，此時由距離構(gòu)成的描述集合中大部分元素將是相對新檢出特征的距離，這些距離對匹配是無效的，甚至會導致誤匹配率的增加。因此需要引入更多的信息加強描述集合中距離元素的“獨特性”，考慮到漫游器與著陸器均工作在平坦的地形環(huán)境中，地形特征近似分布在平面內(nèi)，本文選取相對位置矢量與特征主方向夾角作為輔助信息對相對距離加以區(qū)分。這里特征主方向的定義與ISS算法對主方向的定義一致，為散布矩陣最大特征值對應(yīng)的特征向量。

圖1 特征間有向距離描述Fig.1 Orientative distance description of features

設(shè)三維特征X與其余m個同時檢測出的特征之間的距離集合為{d1,d2,…,dm}，結(jié)合相對位置矢量與特征主方向夾角θi，其中i=1,2,…m，特征X的完整描述為

D={(d1,θ1),(d2,θ2),…,(dm,θm)}

(4)

2.2 三維地形特征匹配

以特征間距離和方向夾角為元素形成的特征描述集合與以往的特征描述不同，無法利用多維矢量距離評估集合的相似性，因此引入集合相似理論計算描述集合的相似性。

首先給出特征描述集合元素相似性的定義

(5)

(6)

式中：τ1和τ2分別為兩個設(shè)定的閾值，M為第m個特征描述集合中的元素數(shù)量。即要求相似度最高與次高的兩個特征描述集合，在相似度方面應(yīng)有一定的差距，同時相同元素個數(shù)應(yīng)在第m個特征描述集合中占有一定的比例。

基于有向距離描述的特征匹配算法流程如下：

1)輸入兩次觀測檢出特征的描述集合，元素相似性閾值τ，以及描述集合相似性閾值τ1和τ2；

3)跳轉(zhuǎn)至步驟2)直至遍歷第一次檢測出的全部特征點。

3 三維地形特征誤匹配剔除

利用上述方法匹配的三維地形特征仍存在一定的誤匹配現(xiàn)象，不能直接用來進行漫游器或著陸器導航，這里給出基于最小中值(Least Median Squares)的誤匹配剔除方法[22]。基于最小中值的誤匹配剔除算法基本流程如下：

(7)

4 仿真校驗

為了測試三維地形特征檢測與匹配算法，生成不同起伏和粗糙程度的模擬地形20幀，模擬地形尺寸為25.6 m×25.6 m，如下圖所示，并模擬LiDAR測量原理對三維地形進行采樣，LiDAR的分辨率設(shè)定為256×256，最大測距距離為50 m。

圖2 三維地形和模擬激光雷達測量數(shù)據(jù)Fig.2 3D terrain and simulated LiDAR measurement data

首先對本文改進的三維地形特征檢測算法進行測試，考察在觀測方向改變時檢測算法的重復檢出率。設(shè)置激光雷達繞三維地形運動，且激光雷達的光軸始終指向三維地形中心，以光軸垂直地形基準平面時的觀測為參考，記此時的觀測方向為0°度，觀測方向每改變10°進行一次測量，直到改變50°。選取特征點鄰域半徑ρ=0.4 m，特征檢測算法閾值τd取值由0.05至0.13。改進的三維特征提取算法平均重復檢出率如圖3所示。重復檢出率定義為重復檢出三維特征數(shù)量與 0°觀測方向條件下檢出特征數(shù)量的比值。特征重復檢出的標準定義為，檢出特征的位置經(jīng)位姿變換后，與0°觀測方向檢出特征位置距離小于ρ/2。利用相同的模擬地形庫比較了ISS算法的特征重復檢出率，如圖4所示，ISS算法特征點鄰域半徑同樣設(shè)定為0.4 m，特征檢測閾值設(shè)定為τ12=τ23=τISS，其中τISS的取值范圍為0.5至0.9。由圖3和圖4可以看出改進的特征檢測算法閾值取0.05時重復檢出率最高，ISS算法閾值取0.9時重復檢出率最高。表1給出了兩種算法平均重復檢出率的具體數(shù)值，改進的特征檢測算法重復檢出率在各觀測方向條件下均高于ISS算法。

圖3 改進的特征檢測算法重復檢出率Fig.3 Feature repeatability rate of the improved feature detection algorithm

進一步對改進的三維地形特征檢測算法在不同噪聲水平下的重復檢出率進行測試。利用20幀模擬地形在0°至50°觀測角度內(nèi)模擬生成觀測數(shù)據(jù)，并疊加不同水平的高斯噪聲，噪聲標準差包括ρ/150，2ρ/150，3ρ/150，4ρ/150和5ρ/150，特征檢測算法閾值設(shè)置為0.05，不同噪聲水平下的平均重復檢出率如圖5所示，與未加噪聲的仿真結(jié)果對比

圖4 ISS算法重復檢出率Fig.4 Feature repeatability rate of the ISS algorithm

觀測方向改變本文算法ISS算法10°0.50810.434120°0.39590.329030°0.31290.253540°0.22520.203350°0.19010.1803

可以看出，噪聲對重復檢出率影響不大。

圖5 不同噪聲條件下的重復檢出率Fig.5 Feature repeatability rate under noise conditions

接著測試基于特征間距離的特征描述方法，并與PFH(Point Feature Histogram)特征描述方法進行比較。由于這里構(gòu)建的特征描述方法是基于有向距離集合的，匹配的判據(jù)是相同元素的數(shù)量是否占優(yōu)，這與基于特征描述矢量距離的匹配方法不同，無法比較兩種方法的PRC (Precision-Recall Curve)。因此這里僅比較兩種描述方法的正確匹配特征數(shù)量，也就是基于特征間有向距離的特征描述方法在匹配時只要相同元素最多則記相應(yīng)特征點匹配，而PFH描述法只要描述矢量距離最近則記相應(yīng)特征點匹配，最后將根據(jù)匹配特征實際位置誤差是否小于鄰域半徑的一半確定是否正確匹配?；谟邢蚓嚯x描述集合的匹配方法元素相同判定閾值設(shè)為0.1，PFH方法計算點云法線的半徑設(shè)為0.1 m，兩種方法的鄰域半徑均為0.4 m。特征檢測方法采用本文改進的方法，檢測閾值設(shè)定為0.05，鄰域半徑0.4 m。圖6給出了PFH描述法和基于特征間距離的描述方法在觀測方向改變且疊加不同噪聲條件下的正確匹配特征平均數(shù)量。

圖6 正確匹配的三維地形特征數(shù)量Fig.6 Number ofthe correct matched features

圖6中虛線為PFH描述方法的平均正確匹配數(shù)，只有在觀測方向改變小于20°，噪聲水平小于3ρ/150的情況下正確匹配數(shù)量超過3個，考慮到導航定位至少需要3個匹配特征，那么PFH描述方法在以激光雷達為敏感器的導航應(yīng)用中是不合適的。圖6中實線為基于特征間有向距離的描述方法在不同噪聲水平下的平均正確匹配數(shù)，當觀測方向改變達到40°時，各噪聲水平對應(yīng)的正確匹配數(shù)降到了10個以下，其余情況均有10個以上的正確匹配。不同噪聲水平下的正確匹配數(shù)量均隨觀測角度增大而快速下降，其原因在于以下兩方面:一是觀測角度的增加使得部分特征由于地形遮擋而無法檢出；二是傾斜的觀測方向使激光雷達測量數(shù)據(jù)不均勻，在傾斜方向的遠端獲得的數(shù)據(jù)更加稀疏，使得特征檢出率下降，且檢出特征位置也不夠精確。因此在導航任務(wù)中，尤其是以激光雷達獲取的稀疏三維點云為觀測數(shù)據(jù)的導航任務(wù)，盡量在觀測方向改變小于20°的條件下應(yīng)用基于有向距離的特征描述和匹配方法。

最后，應(yīng)用第3節(jié)中給出的最小中值法進行了漫游器導航仿真。利用最小中值法可以剔除誤匹配特征，同時也給出了兩次觀測間相對運動的粗略估計，在初始位置上累加相對運動即可得到漫游器位置。在圖2給出的模擬地形上，設(shè)置了20個路徑點，并在這些點上進行導航觀測，應(yīng)用三維地形特征提取和匹配方法在相鄰的觀測數(shù)據(jù)中確定匹配點對，借助最小中值法估計幀間運動。

圖7 漫游器定位結(jié)果Fig.7 Positioning result of the rover

圖7給出了基于幀間運動估計的漫游器定位結(jié)果，開始的幾次定位結(jié)果較為精確，但是由于每次運動估計的誤差不斷累積，在后面的路徑點上導航誤差不斷增大。這個簡單的導航仿真試驗僅校驗了在導航定位中應(yīng)用本文提出的地形特征提取與匹配算法的可行性，如果想獲得更加精確的定位結(jié)果，可以在若干次運動估計后利用全局地圖進行特征匹配，提高漫游器全局定位精度。

5 結(jié) 論

針對行星漫游器地形相關(guān)導航應(yīng)用，提出了一種在稀疏點云數(shù)據(jù)中提取和匹配地形特征的方法。該方法利用特征點鄰域散布矩陣最大和最小特征值比檢測具有凹陷或凸起的地形，提高了平坦地形中三維特征的重復檢出率，為了解決稀疏點云特征匹配困難，提出了基于特征間有向距離的特征描述方法，通過與PFH三維特征描述方法的對比可以看出，利用多個特征間的有向距離集合能夠有效描述地形特征，滿足以稀疏點云為觀測數(shù)據(jù)的導航應(yīng)用需求。

參 考 文 獻

[1] 崔平遠, 袁旭, 朱圣英, 等. 小天體自主附著技術(shù)研究進展[J]. 宇航學報, 2016, 37(7):759-767.[Cui Ping-yuan, Yuan Xu, Zhu Sheng-ying, et al. Research progress of small body autonomous landing techniques[J]. Journal of Astronautics，2016, 37(7):759-767.]

[2] 吳偉仁, 周建亮, 王保豐, 等. 嫦娥三號“玉兔號”巡視器遙操作中的關(guān)鍵技術(shù)[J]. 中國科學:信息科學, 2014, 44(04): 425-440. [Wu Wei-ren, Zhou Jian-liang, Wang Bao-feng, et al. Key technologies in the teleoperation of Chang’E-3 “Jade Rabbit” rover[J]. Scientia Sinica Informationis, 2014, 44(04): 425-440.]

[3] Lauretta D S, Balram-Knutson S S, Beshore E, et al. OSIRIS-REx: Sample return from asteroid (101955) Bennu[J]. Space Science Reviews, 2017, 212(1-2): 925-984.

[4] 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.

[5] Kok-Lim L.Linear least-squares optimization for point-to-plane ICP surface registration[R]. Technical Report TR04-004, Department of Computer Science, University of North Carolina at Chapel Hill, February 2004.

[6] Pomerleau F, Colas F, Siegwart R, et al. Comparing ICP variants on real-world data sets[J]. Autonomous Robots, 2013, 34(3):133-148.

[7] Chen H, Bhanu B. 3D free-form object recognition in range images using local surface patches[J]. Pattern Recognition Letters, 2007, 28(10): 1252-1262.

[8] Yu Z. Intrinsic shape signatures: a shape descriptor for 3D object recognition[C]. International Conference on Computer Vision Workshops, Kyoto, Japan, September 27- October 4, 2009.

[9] Mian A, Bennamoun M, Owens R. On the repeatability and quality of keypoints for local feature-based 3D object retrieval from cluttered scenes[J]. International Journal of Computer Vision, 2010, 89(2-3), 348-361.

[10] Sun J, Ovsjanikov M, Guibas L. A concise and provably informative multi-scale signature based on heat diffusion[J]. Eurographics Symposium on Geometry Processing, 2009, 28(5): 1383-1392.

[11] Unnikrishnan R, Hebert M. Multi-scale interest regions from unorganized point clouds[C]. Computer Society Conference on Computer Vision and Pattern Recognition Workshops, Anchorage, AK, United States, June 23-28, 2008.

[12] Zaharescu A, Boyer E, Varanasi K, et al.Surface feature detection and description with applications to mesh matching[C]. Computer Society Conference on Computer Vision and Pattern Recognition Workshops, Miami, FL, United States, June 20-25, 2009.

[13] Mian A S, Bennamoun M, Owens R A. On the repeatability and quality of keypoints for local feature-based 3D object retrieval from cluttered scenes[J]. International Journal of Computer Vision, 2010, 89(2-3): 348-361.

[14] Boyer E, Bronstein A M, Bronstein M M, et al. SHREC 2011: Robust feature detection and description benchmark[C]. Eurographics Workshop on 3D Object Retrieval, Llandudno, United Kingdom, April 10, 2011.

[15] 崔平遠, 馮軍華, 朱圣英, 等. 基于三維地形匹配的月球軟著陸導航方法研究[J]. 宇航學報, 2011, 32(3):470-476.[Cui Ping-yuan, Feng Jun-hua, Zhu Sheng-ying, et al. 3D terrain feature matching based navigation for lunar soft landing[J]. Journal of Astronautics，2011, 32(3):470-476.]

[16] Tombari F, Salti S, Di L. Performance evaluation of 3D keypoint detectors[J]. International Journal of Computer Vision, 2013, 102(1-3): 198-220.

[17] Johnson A E, Hebert M. Surface matching for object recognition in complex three-dimensional scenes[J]. Image and Vision Computing, 1998, 16(9-10): 635-651.

[18] Frome A, Huber D, Kolluri R, et al. Recognizing objects in range data using regional point descriptors[C]. European Conference on Computer Vision, Prague, Czech Republic, May 11-14, 2004.

[19] Tombari F, Salti S, Di S. Unique signatures of histograms for local surface description[C]. European Conference on Computer Vision, Crete, Greece, September 5-11, 2010.

[20] Radu R B, Marton Z C, Blodow N, et al.Persistent point feature histograms for 3D point clouds[C]. 10th International Conference on Intelligent Autonomous Systems, Baden-Baden, Germany, July 23- 25, 2008.

[21] Guo Y L, Bennamoun M, Sohel F, et al. A comprehensive performance evaluation of 3D local feature descriptors[J]. International Journal of Computer Vision, 2016, 116(1): 66-89.

[22] 馬頌德, 張正友. 計算機視覺:計算理論與算法基礎(chǔ)[M]. 科學出版社, 1998: 160-161.