李海波， 曹云峰， 丁 萌， 莊麗葵(南京航空航天大學， 江蘇 南京 210016)

1 引 言

為了對火星進行更全面的科學考察，通常需要探測器著陸到火星表面。探測器在著陸過程中，應避免著陸到坡度較大的斜坡上，以免造成側滑或傾翻[1]。因此探測器在著陸階段，需要對候選著陸區(qū)域的坡度進行估計?；鹦潜砻娲嬖谳^大的巖石[2],巖石會使巖石區(qū)域的測量數(shù)據相對周圍區(qū)域測量數(shù)據產生較大的突兀值，影響坡度估計。目前,對坡度的估計方法主要有基于機器視覺的方法與基于激光雷達(light detection and ranging，LIDAR)的方法。機器視覺在計量領域有多種應用[3～5]，其設備具有體積小、重量輕的優(yōu)點，但獲得的數(shù)據受光照變化影響大;激光雷達可以應用到不同測量目標上[6,7]，其測量數(shù)據不受光照影響，是一種理想的坡度估計信息源;文獻[8]提出了一種利用激光雷達測量數(shù)據進行坡度估計的方法，該方法采用最小二乘法擬合平面，擬合平面的傾角即為該區(qū)域的坡度;文獻[9]也提出了一種基于激光雷達測量的高程數(shù)據進行坡度估計的方法，該方法采用最小平方中位數(shù)法求取局部地形的法向量，法向量與天體重力方向的夾角即認為是局部地形的坡度，但該方法需要預先設定一些先驗參數(shù);文獻[10]針對激光雷達測量的數(shù)據利用最小二乘法擬合一個平面，求取該平面的法向量，定義法向量與當?shù)刂亓铀俣确较蛑g所夾的銳角即為局部地形坡度，也就是擬合平面與基準面之間的夾角;此外其他相關文獻中關于坡度估計的方法大都基于上述方法，如文獻[11]等。但這些方法對坡度進行估計均存在估計誤差較大的問題，為此，本文研究了一種基于三維點云數(shù)據聚類與隨機搜索最優(yōu)擬合平面的坡度估計方法，可以顯著降低估計誤差，提高坡度估計精度。

2 坡度估計原理

2.1 測量數(shù)據點的稀疏表示

首先定義測量區(qū)域空間坐標系，其定義及數(shù)據點三維表示如圖1所示?？臻g直角坐標系在三維激光掃描測量中起著關鍵作用[12]，以測量區(qū)域的左下角為原點，x軸與水平邊重合指向正東；y軸與x軸垂直，指向正北；z軸與x軸和y軸垂直并按右手定則確定，見圖1(a)。通過空間坐標系，激光雷達的測量數(shù)據可用三維坐標表示，即Di=[xi,yi,zi]，其中，D表示數(shù)據點；i為測量點數(shù)據序號;x、y、z分別表示測量點數(shù)據在三軸上的坐標，見圖1(b)，圖中模擬的數(shù)據點可用三維坐標表示。

圖1 坐標系定義及數(shù)據點三維表示

X[x1…xN]=[X1…Xn]Γ

(1)

式中：Xl∈R3×Nl；Γ∈RN×N，是未知的置換矩陣。由于子空間的基是未知的，無法確定測量數(shù)據點屬于哪個子空間。為解決這一問題，首先將數(shù)據點進行稀疏化表示，使每一個數(shù)據點可以用空間中少數(shù)其他數(shù)據點進行表示，即

xi=Xci,cii=0

(2)

式中：ci[ci1ci2…ciN]T；cii=0為限制條件，防止得到平凡解。由于在子空間中數(shù)據點的數(shù)量大于數(shù)據點的維數(shù)，解的形式不唯一，即每個數(shù)據點由其他數(shù)據點表示的方式不唯一，具有多種方式。在這些解中，存在一種稀疏解ci，其解的非零部分與來自同一子空間的數(shù)據點xi相對應[13]， 即對于子空間Sl中的任一數(shù)據點xi，可以用來自同一子空間中的其他數(shù)據點組合表示。因此，離散數(shù)據點的稀疏表示能夠發(fā)現(xiàn)屬于同一子空間的數(shù)據點。由于式(2)中存在多個解，故采用最小化解的1范數(shù)方式獲得稀疏解，即

(3)

式(3)的解可通過凸規(guī)劃工具求得[14～16]。通過式(3)，將所有數(shù)據點的稀疏優(yōu)化寫成矩陣形式：

(4)

式中:C[c1c2…cN]∈RN×N，其列向量對應數(shù)據點的稀疏表示。式(4)的解即為所有數(shù)據點的稀疏表示，解中非零元素與空間中數(shù)據點相對應。

2.2 測量數(shù)據點分割

通過構造好的非負對稱相似度矩陣W構造對角矩陣A，形式如下：

(5)

式中對角元素aii(i=1,…,N)為相似度矩陣W中列向量ci所有元素之和。利用對角矩陣A構建Laplacian矩陣，即L=I-A-1/2WA-1/2，I為單位矩陣。對構建的Laplacian矩陣L進行SVD(singular value decomposition)分解，對分解后得到的右正交矩陣V應用K-means算法進行聚類，從而完成對數(shù)據點的分割。

2.3 子空間測量點平面擬合

對分割后的子空間測量點單獨進行平面擬合，設某一子空間測量數(shù)據點預確定的平面的方程為：

ax+by+cz=d

(6)

式中：a、b、c為所設定平面單位法向量中的元素值，即單位法向量為(a,b,c)，則有a2+b2+c2=1；d表示原點到所設定平面的距離，d≥0。

在子空間中隨機選取一個數(shù)據點，計算與其歐式距離最近的4點，設這5點為{(xi,yi,zi),i=1,…,5}，由式(6)可知每個點到平面的距離為

hi=|axi+byi+czi-d|

(7)

式中h表示點到平面的距離。

(8)

當e的值最小時，平面擬合最好。為求取e的最小值，根據約束條件：a2+b2+c2=1，利用拉格朗日乘子法構造函數(shù)如下：

λ(a2+b2+c2-1)

(9)

式中λ為拉格朗日乘子。

將式(9)對d求偏導,并令其為0，即

(10)

整理得

(11)

將式(11)代入式(7)得

(12)

將式(11)代入式(9)得

(13)

將式(13)分別對a、b、c求偏導,并令其為0，

(14)

將式(14)進一步整理轉換得：

(15)

Bx=μx

(16)

B為3×3實對稱矩陣，由矩陣理論可知：

(17)

式中(，)表示兩個向量做點積運算。由于a2+b2+c2=1，所以(x,x)=1，則

(18)

當e取最小值時擬合的平面最好，所以μ應取最小值，其對應的特征向量為a、b、c的值。

利用正交三角分解可求出矩陣B的所有特征值，設最小特征值為μmin，I為單位矩陣。求其對應的特征向量，令

(B-μminI)x=0

(19)

通過式(19)求得的非零解即為對應的特征向量，從而獲得參數(shù)a、b、c的值，即平面的法向量元素值。將所求的a、b、c值代入式(11)中，可求出d的值。

設子空間中另一數(shù)據點為Dj={(xj,yj,zj),j=1,…,n,j≠i}，由式(7)可知，其與平面的距離為hj。當hj的值小于設定閾值時，則該點屬于統(tǒng)計的點； 當hj的值大于或等于設定閾值時，該點不屬于統(tǒng)計的點。

按此方法計算子空間中其余點與該平面的距離，統(tǒng)計與平面距離小于設定閾值的點，當統(tǒng)計的數(shù)目大于預先設定值時，便認為該平面為最優(yōu)擬合平面，其法向量為(a,b,c)；若統(tǒng)計的數(shù)目小于或等于預先設定值，便認為該點不是最優(yōu)擬合平面，需重新隨機選取一點，重復上述步驟，直至統(tǒng)計的數(shù)目大于預先設置值為止。同理，用該方法可以求得另一子空間數(shù)據點的最優(yōu)擬合平面及法向量。

2.4 坡度角估計

在獲得兩個最優(yōu)擬合平面后，可以通過其法向量夾角求得兩個平面的夾角，設兩平面的法向量分別為n1、n2，則其夾角可通過下式計算。

(20)

式中θ為兩法向量之間夾角。根據幾何關系，θ在數(shù)值上等于兩平面夾角，即坡度角。

3 仿真實驗及分析

隨機選取兩個斜坡角度進行算法測量實驗，模擬生成的待測量斜坡角度如圖2所示。

如圖2所示兩種斜坡角度對應的模擬激光雷達測量數(shù)據點如圖3所示。

圖2 模擬生成的待測量斜坡角度

圖3 斜坡角度測量數(shù)據點

對于圖3的測量數(shù)據點，首先利用2.1節(jié)所述方法進行稀疏表示；在此基礎上，利用2.2節(jié)所述方法進行數(shù)據分割，劃分子空間數(shù)據點；對子空間數(shù)據點使用2.3節(jié)所述方法進行平面擬合，求出平面法向量；根據所求法向量，利用公式(20)計算兩平面夾角，即坡度角，平面擬合及計算結果如圖4所示。

對斜坡角度測量結果進行統(tǒng)計分析見表1。

表1 測量結果分析

由表1可以看出，本文所述方法測量誤差極小，測量相對誤差較小。

考慮到在實際工作中坡度可能很不規(guī)則及火星地表存在較大巖石等情況，在34°斜坡上模擬添加存在上述情況產生的突出測量值，如圖5(a)中橢圓內的測量值。利用2.2、2.3節(jié)所述方法對數(shù)據點進行分割與尋找最優(yōu)平面擬合數(shù)據點，數(shù)據點尋找結果如圖5(b)所示，可以看出突出測量值數(shù)據點沒有被統(tǒng)計到擬合數(shù)據點范圍內，利用圖5(b)所示的數(shù)據點進行平面擬合，結果如圖5(c)所示。

從圖5(c)可以看出，坡度測量結果為33.970 7°，與圖4(b)的測量結果相同，這表明本方法對突出的測量值具有一定的魯棒性。

圖5 存在突出測量值的坡度估計

4 與常用方法比較

利用激光雷達測量數(shù)據進行坡度估計的常用方法是最小二乘法(least square method,LSM)，將文中提出的方法與最小二乘法進行比較。隨機選取一個坡度，模擬生成的斜坡角度及測量數(shù)據點如圖6所示。分別用最小二乘法和本方法對測量數(shù)據進行處理，計算結果如圖7所示。對測量結果進行比較分析，結果如表2所示。從表2可以看出，本文提出的方法測量相對誤差較??；與LSM方法比較，測量誤差得到了顯著降低。

圖6 44°斜坡及測量數(shù)據點

圖7 測量結果比較

表2 不同方法測量結果比較分析

5 結 論

為降低火星探測器著陸時對坡度估計的誤差，研究了一種基于三維點云數(shù)據聚類與隨機搜索最優(yōu)擬合平面的坡度估計方法。對于三維點云測量數(shù)據，首先進行數(shù)據點的稀疏表示；然后根據稀疏表示系數(shù)對數(shù)據點進行聚類與分割，劃分子空間；對子空間中的數(shù)據點進行平面擬合，隨機搜索最優(yōu)擬合平面；根據最優(yōu)擬合平面計算其法向量及向量間夾角，完成坡度估計。實驗表明：本文方法可以對坡度進行較為準確的估計,與常用的坡度估計方法相比較，可以顯著降低估計誤差，提高坡度估計精度，為探測器著陸時決策是否需要規(guī)避機動提供更加準確的信息。

[參考文獻]

[1] Lunghi P, Ciarambino M, Lavagna M. A multilayer perceptron hazard detector for vision-based autonomous planetary landing[J].AdvancesinSpaceResearch,2016,58(1):131-144.

[2] Craddock R A, Golombek M P. Characteristics of terrestrial basaltic rock populations:Implications for Mars lander and rover science and safety[J].Icarus,2016,274:50-72.

[3] 劉慶民,張蕾,吳立群,等.基于機器視覺的非均勻分布點圓度誤差評定[J].計量學報,2016,37(6):567-570.

[4] 劉源泂,孔建益,徐福軍,等.基于圖像金字塔的鋼板表面深度信息提取方法[J].計量學報,2015,36(4):356-359.

[5] 朱奇光,張興家,陳衛(wèi)東,等.基于顏色矩的改進尺度不變特征變換的移動機器人定位算法[J].計量學報,2016,37(2):118-122.

[6] 王濤,仲思東.彈丸三維輪廓激光掃描測量方法[J].計量技術,2015,(2):14-18.

[7] 魏凱,宋述古,劉子勇.基于三維激光掃描原理的球形罐容量計量方法研究[J].計量學報,2015,36(6):607-609.

[8] 梁棟,王鵬基,劉良棟.一種基于LIDAR的精確月球軟著陸目標點選定方法[J].空間控制技術與應用,2009,35(6):24-29.

[9] 姜肖楠,韓誠山,李祥之.多傳感器月面障礙模糊識別方法[J].計算機仿真,2013,30(4):97-102.

[10] 張澤旭,王衛(wèi)東,崔平遠,等.一種行星軟著陸地形風險評估方法[J].哈爾濱工業(yè)大學學報,2011,43(5):25-29.

[11] 于正湜,朱圣英,崔平遠.基于LIDAR的月球著陸區(qū)評估與選擇方法[C]// 中國宇航學會深空探測技術專業(yè)委員會第九屆學術年會.杭州,2012.

[12] 李明慈,黃桂平,李冰,等.坐標轉換在應用中的問題及解決方法[J].計量技術,2015,(8):25-28.

[13] Elhamifar E, Vidal R. Sparse subspace clustering:algorithm,theory,and applications[J].IEEETransactionsonPatternAnalysis&MachineIntelligence,2013,35(11):2765-2781.

[14] Boyd S, Vandenberghe L, Faybusovich L. Convex Optimization[J].IEEETransactionsonAutomaticControl,2006,51(11):1859-1859.

[15] Kim S J, Koh K, Lustig M,etal. An Interior-Point Method for Large-Scale l1-Regularized Least Squares[J].IEEEJournalofSelectedTopicsinSignalProcessing,2007,1(4):606-617.

[16] Boyd S, Parikh N, Chu E,etal. Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[J].FoundationsandTrendsinMachineLearning,2011,3(1):1-122.