一種探測器自主安全著陸區確定及精確避障的策略

金 濤 章登義 蔡 波

1(武漢大學計算機學院 武漢 430072)2(武漢大學國家網絡安全學院 武漢 430072)(792042768@qq.com)

發射星際探測器并實現在目標星球表面的軟著陸是我國探索深空、發展航天技術的重點需求.自主安全著陸是未來著陸器進行安全著陸的重要能力，這種能力主要體現在著陸器能夠對傳感器采集的信息進行快速處理、判斷并自動形成新的著陸指令.為了達到預定的科學目標，地面預先選擇的軟著陸區并不是完全平坦的，它存在著一定分布的石塊或隕石坑.這些石塊和隕石坑的存在會給探測器安全著陸帶來危險，因此，著陸器需要具有發現和識別障礙并進行機動避障的能力.為了實現下落過程中對障礙物的躲避,需要使用探測器上攜帶的電荷耦合器件(charge coupled device,CCD)相機進行拍照，根據拍攝的月面圖像來判斷月表的地形.

至今為止,各國實施深空探測任務的著陸器還不具有對著陸區的潛在危險進行檢測和規避的閉環控制能力.為了使深空探測器具有自主軟著陸及避障的能力，噴氣推進實驗室(Jet Propulsion Laboratory,JPL)以未來火星和小天體著陸任務為背景，開發了基于光學敏感器的精確導航和障礙規避算法[1-3].歐洲航天局(European Space Agency,ESA)也開發了具有高精度著陸和自主避障能力的新一代著陸系統，以探索科學價值高且危險的區域[4-5].從當前的研究現狀可以看出，基于灰度圖像的避障技術和基于激光雷達(light detection and ranging,LIDAR)的避障技術已經被考慮用于探測器自主障礙中，并且得到了進一步發展[6].文獻[7]綜合使用多因子分析評估了月球探測的著陸區域和巡航路線.

1 相關工作

目標天體的地表障礙主要是巖石、隕石坑等.對于巖石的檢測，Ma等人[8]提出了邊緣光流算法，該算法是基于紋理的圖像分割，光照影響較大；Matthies等人[9]通過立體視覺三維重構來進行障礙檢測，該方法較為準確可靠，可以檢測的障礙類型較多，但是涉及的數據計算量大，導致檢測速度較慢；Huertas等人[10]提出基于模版匹配的自主障礙檢測方法，這種方法的有效性取決于模板建模的精確程度，模板建模越精確，但匹配所需要的時間越長.對于最大安全區域的計算及選擇并沒有相應的算法；而文獻[11]也只是對安全著陸圓進行了粗略的估計.為了達到預定的科學目標，初避障通過圖像匹配識別所選定的初選著陸區域會存在一定的誤差，導致探測器安全著陸存在一定的風險，因此，探測器需要有更高精度的自主障礙識別與規避能力.為了實現探測器的高精度避障，需要對初選著陸區域進行高精度三維地形重構，然后根據重構后的模型，進行螺旋搜索，根據探測器安全著陸閾值來選定目標著陸點.

在精避障階段，高精度的障礙識別算法還未能成熟的應用于航天探測及仿真項目中.南加州大學的Leroy等人[12-14]都對圖像匹配與識別的方法進行了研究與論述.在圖像匹配與識別方法之后，激光測距儀的思想被學者們提出，即通過儀器獲得探測器到目標點之間的距離，然后對采樣值進行處理以獲得目標星球表面信息，從而提取出星表特征進行障礙檢測.這一思想最初是由JPL實驗室的Andrew等人[15]提出的，之后JPL實驗室又對該思想進行了進一步的研究[16]，設定了探測器障礙檢測的相關參數和技術指標.之后美國 JPL實驗室、MIT Lincoln實驗室以及美國、英國和加拿大的多所大學開始對多傳感器協同工作避障的方法展開研究.

本文根據中國深空探測的現狀，提出了一種基于CCD著陸相機生成表面遙感影像并進行著陸區巖石檢測與安全著陸區確定的方法.這種方法不僅提出了遙感影像的生成方法，同時對著陸區巖石的檢測、識別進行了研究，并且提出了確定最大安全著陸區的算法.在精避障階段本文提出了一種基于三維重構的探測器精避障螺旋掃描方法.這種方法可以通過激光三維成像儀獲取的數據，構造高精度的表面初選著陸區域三維網格模型，并針對模型的分塊進行高精度螺旋掃描，結合探測器安全著陸閾值選定目標著陸點.其整個流程構成了一個完整的探測器軟著陸精避障過程實現的解決方法.

Fig.1 The technical roadmap of autonomous safe landing simulation algorithm圖1 自主安全著陸仿真算法技術路線

2 系統模型

軟著陸的最后垂直著陸段將進行障礙物的檢測與規避，由于對垂直著陸段的高度有一定的約束，因而希望能夠快速地生成遙感影像，有效地檢測和識別巖石障礙物，相對精確地確定安全著陸區，以滿足制導控制的需要.本文首先提出了基于CCD相機的遙感影像的生成；然后用障礙識別算法對著陸表面的巖石、彈坑等障礙進行提取，主要是通過基于Canny理論的自適應邊緣檢測方法提取巖石輪廓，并利用橢圓擬合方法確定巖石的大小位置，簡單有效；最后提出算法確定最大無障礙著陸圓，對著陸點的選擇提供參考.在確定巖石輪廓的基礎上，搜索最大的安全著陸區.初步確定了著陸區域后，在探測器精確避障階段.首先利用Delaunay三角剖分算法以及分塊插值策略，對初選著陸區域進行高精度三維建模，形成精細的星表DEM模型；然后針對高精度的模型，根據探測器當前的星下點位置，進行分塊螺旋搜索，根據設定的探測器安全著陸閾值來選定安全的目標著陸點.螺旋搜索獲得的目標著陸點是距離探測器當前位置最近的符合著陸閾值的著陸點.其技術方案如圖1所示.

2.1 遙感影像生成

探測器在軟著陸過程中，在下落到10～2 000 m之間時，其攜帶的CCD相機會對星球表面進行拍照.本節根據探測器所在位置，確定所攜帶的CCD相機拍攝的區域，基于光線跟蹤算法結合高程數據對區域內的遙感影像圖片進行仿真模擬.光線跟蹤主要是判斷光線與球面的相交情況，如圖2所示：

Fig.2 The illustration of ray tracing圖2 光線跟蹤示意圖

根據代數學的方法，設置光線方程：

P=P0+tV，

(1)

其中，t為待定參數，求解P的方程：

|P-O|2-r2=0.

(2)

消去O，可以得到：

|P0+tV-O|2-r2=0.

(3)

根據式(3)得到關于t的一元二次方程：

at2+bt+c=0.

(4)

求解方程組：

(5)

得到t，即可求得點P.

衛星相機拍照的效果主要受光照、氣候以及相機特性影響.這里只考慮光照和相機的影響.通過計算光照模型以及相機拍照范圍的計算，算法流程如圖3所示：

Fig.3 Flow chart of remote sensing image generation圖3 遙感影像生成流程圖

2.2 障礙檢測

本文提出的是基于邊緣檢測用橢圓擬合探測器目標表面的巖石及隕石坑等障礙進行檢測.邊緣是圖像最基本的特征.所謂邊緣是指圖像像素灰度有階躍變化或屋頂狀變化的像素的集合.本文采用的是基于Canny算子的自適應邊緣檢測方法[17].這種方法將整幅圖像分割為若干子圖像，再根據非極大值抑制后的結果自適應地設定各子圖像的高、低閾值.自適應Canny算法即二維高斯函數：

(6)

對圖像進行平滑，計算像素的梯度幅值和方向：

(7)

其中，Px(i,j)為x方向偏導數，Py(i,j)為y方向偏導數，P135°(i,j)為135°方向偏導數，P45°(i,j)為45°方向偏導數.

(8)

通過圖像的直方圖進行高低閾值的確定.獲得圖像的邊緣信息之后，然后利用橢圓擬合的方法來估測隕石坑的大小和位置.

巖石檢測方法則是基于陰影檢測，該檢測技術主要分為2個過程，即陰影檢測和巖石區域確定.陰影檢測過程就是圖像分割過程，圖像閾值分割是一種廣泛使用的圖像分割技術，它利用了圖像中要提取的目標物與其背景在灰度特性上的差異，把圖像視為具有不同灰度級的2類區域的組合，選取一個合適的閾值，以確定圖像中每一個像素點應該屬于目標還是背景區域，從而產生相應的位置圖像.這里我們用橢圓擬合的方式利用巖石陰影和太陽光角度來估計巖石區域的范圍.

巖石的大體形狀可以通過巖石陰影估計出來，巖石和陰影模型如圖4所示.最合適的擬合橢圓能夠有效的表示巖石的形狀和大小[18].

Fig.4 The rock and shadow model圖4 巖石和陰影模型

橢圓擬合技術是指找到一個參數集合使得給定的數據點到橢圓的距離最短，求解過程:

FK(X)=KX=ax2+bxy+cy2+dx+ey+f=0,

(9)

其中

(10)

這里定義FK(X)為點(x,y)到曲線FK(X)=0的代數距離.可以計算多個點到曲線代數距離平方和的最小值來擬合接近一般的圓錐曲線：

(11)

為了避免零解或者不同K值表示相同的圓錐曲線，對參數向量K應該有一定的限制.為了在擬合橢圓的同時保證解決線性最小平方問題的效率，這里對向量K有一定的限制從而使得一般二次曲線表示為一個橢圓.最適合的限制眾所周知，即b2-4ac≥0.盡管在一般情況下，這種不等式的約束條件實行是比較困難的，在這種情況下，我們可以自由任意地縮放參數從而簡單地引入一定的等式約束4ac-b2=1.這個二次的等式限制也可以表示為矩陣形式KTCK=1，即：

(12)

根據式(9)的定義，我們構造矩陣D：

(13)

可以重寫為

SK=λCK,

(14)

KTCK=1，

(15)

(16)

2.3 著陸區的確定

將整個圖像區域分為2部分：障礙區域和非障礙區域，目標是要求在非障礙區域尋找最大圓.關鍵一點是要確定最大區域的質心.由于整個圖像僅分成2部分，可以用黑白顏色值來表示，白色區域表示非障礙區域.而圖像形態學經常在通過閾值化獲得的二進制圖像中完成，因而本節將介紹利用形態學操作在白色區域中確定最大區域的質心.

基于數學形態學的圖像處理是用具有一定形態的結構元素去度量和提取圖像中的對應形狀以達到對圖像分析和識別的目的[20].在形態學中，結構元素的形態決定了運算所提取的圖像中形態信息.這里利用數學形態學的膨脹和腐蝕基本變化操作，分割出其中的白色區域，并確定白色區域中連續的最大區域.假定結構元素為g，圖像為f，則灰度形態膨脹[21]：

(17)

灰度形態腐蝕:

(18)

膨脹是指將一些圖像與核進行卷積.核可以是任何的形狀或大小，核可以視為模版或掩碼，膨脹是求局部最大值的操作.腐蝕是膨脹的反操作.通過膨脹黑色區域或是腐蝕白色區域，不斷逼近縮小白色區域從而確定最大的白色區域；然后利用橢圓擬合來確定不規則區域的質心作為最大圓的中心；最后從該中心擴散，尋找能夠僅僅覆蓋白色區域的最大圓.

算法1.基于數學形態學的著陸區確定算法.

① 設置搜索半徑radus的初值，一般設為圖片大小較短邊值的一半；

② 以radus為搜索半徑，找到質心點center(x,y)中心，在區域(x±radus,y±radus)內開始搜索；

③ 若碰到障礙區域的點Pointi(xi,yi)，設置radus=distance(pointi,center)-1，轉步驟②，否則繼續搜索；

④ 設radus_found=radus，radus_found≤0表示算法搜索失敗(但由于是在最大連續非障礙區域中尋找，不會出現radus_found≤0)；

⑤ 確認以center(x,y)為中心、radus_found為半徑的圓形區域均為非障礙區.

在確定初始的著陸區后，采用圖像金字塔繼續對選中的區域進行處理.圖像金字塔是以多分辨率來解釋圖像的一種有效的結構.一幅圖像的金字塔是一系列金字塔形狀排列的分辨率逐步降低的圖像集合.

3 三維重構算法

在著陸區域基本確定之后的精避障階段，探測器在精避障初始時刻，開啟激光三維成像儀，對初選目標著陸區域進行激光三維成像，測量目標范圍內散亂點的高程數據，然后建立初選目標著陸區域的TIN模型，再根據建立的TIN模型，利用插值算法，計算出規則網格模型的坐標及高程值，形成初選著陸區域的DEM規則網格模型，最后通過三維引擎對規則網格模型進行渲染，繪制出初選著陸區域的三維地形網格.其主要流程如圖5所示：

Fig.5 Three-dimensional reconstruction algorithm flow圖5 三維重構算法流程

3.1 三維TIN模型構建

算法2.TIN模型構建算法(如圖6所示).

Fig.6 The TIN model building algorithm圖6 TIN模型構建算法

① 初始化三角網格，尋找散點集最外圍的4個點構建矩形外包圍框，連接任意一條對角線，對劃分出的2個三角形及其邊編號.

② 判斷是否有未處理的散亂點，如果所有散亂點均處理完畢，轉到步驟⑥.

③ 插入一個新的散亂點，得到所有外接圓包含該點的三角形，記錄三角形的編號.

④ 遍歷步驟③ 中記錄的三角形，刪除公共邊，形成多邊形空腔.

⑤ 連接點P與多邊形空腔的各個頂點，形成新的TIN三角網，轉到步驟②.

⑥ 遍歷所有邊，刪除包含矩形外包圍框中點的邊，獲得最終的三維TIN模型.

3.2 DEM規則網格構建

形成不規則網格TIN模型后，需要根據實際仿真的需求，確立所需DEM規則網格模型的精度，從而確定待插值DEM規則網格點的坐標.然后，開始進行DEM規則網格的構建.主要分2個步驟：

1)計算待插點在不規則三角網格模型中的三角形編號；

2)根據記錄編號的三角形各個頂點的高程值，利用插值算法計算出待插點的高程值，形成DEM規則網格模型.

3.2.1 判斷點在TIN模型三角形中算法

通過插值算法計算DEM規則網格點的高程值，首先需要判斷規則網格點在TIN模型的哪個三角形中，針對N個散亂點所構建的TIN模型一共有2N-6個三角形，且DEM規則網格模型的精度要求較高，若每個待插值點都掃描一遍TIN模型中的所有三角形，那么重復冗余的計算量是相當大的.

為了減少搜索的次數，降低程序的冗余度，采用TIN模型分塊搜索的策略，即將TIN模型按照經緯度分為大小相等的正方形塊，然后只掃描遍歷一次各個三角形的頂點，將每塊中包含的三角形的編號記錄下來(只要三角形有一個頂點在該塊中則認為該塊包含該三角形)，之后判斷規則網格點時，先根據經緯度確定規則網格點在TIN模型的哪個分塊中，然后將點P所在的塊以及其鄰接塊所包含的三角形作為搜索的范圍.針對點所在的塊的位置：內部塊，需要搜索9個分塊；邊界塊，需要搜索6個分塊；頂角塊，需要搜索4個分塊.由此，可以大大減少TIN模型中三角形搜索量，提高DEM規則網格模型的構建效率，算法示意圖如圖7所示：

Fig.7 The judgment point algorithms in TIN model triangle圖7 判斷點在TIN模型三角形中算法

3.2.2 DEM內插建網算法

算法3.DEM內插建網算法(如圖8所示).

① 當AB與AC不平行于x軸時，過待插值規則網格點P作1條平行于x軸的直線與AB和AC(或其延長線)分別相交于點L和點R.當AB和AC其中有1條邊平行于x軸時，則過點P作平行于x軸的直線與另外2條邊分別相交于點L和點R；

② 計算點L和點R的坐標，然后計算點L和點R的高程值:

ZL=ZA+|ZB-ZA|×|XL-XA|/|XB-XA|,

(19)

ZR=ZA+|ZC-ZA|×|XR-XA|/|XC-XA|.

(20)

③ 依照點P的坐標值，計算出點P的高程值:

ZP=ZL+|ZR-ZL|×|XP-XL|/|XR-XL|.

(21)

4 精避障掃描策略

由于初避障過后，探測器是位于初選著陸區域中心上方的，為了使精避障目標著陸區域搜索的時間效率最高，我們采用螺旋搜索的策略，即從探測器的星下點開始搜索計算，計算的最小單位是探測器所要求的控制精度的范圍作為邊長的一個小正方形分塊，然后逐層向外螺旋式擴展搜索，一旦搜索到符合需求的著陸點分塊，就停止搜索.這樣搜索到的目標著陸點離探測器的距離最近，且搜索時間最短，可以較優的實現探測器軟著陸精避障的目的.策略如圖9所示.

針對螺旋搜索算法流程中判斷分塊是否符合軟著陸需求的問題，我們需要根據探測器的安全著陸參數來采取相應的安全判定規則.一般情況下，主要通過最大最小高度差h1和h2、坡度δ和粗糙度α來衡量當前區域是否符合軟著陸的需求：

(22)

在判斷式(22)閾值的過程中，按照最大最小高度差、坡度和粗糙度的順序依次判定，當某一項不符合探測器安全著陸參數的要求時，立刻停止搜索當前塊，不再計算后面的參數.這樣可以縮短判定安全區域的時間，一旦當前塊不符合著陸要求，則立即計算下一塊的參數，從而盡可能地優化精避障安全著陸區域判定的性能.

4.1 最大最小高度差

最大高度差和最小高度差為當前分塊的DEM規則網格點中的最大高程值和最小高程值與當前分塊的DEM規則網格點的平均高程值的差.根據上述定義，首先求得當前分塊的DEM規則網格點的平均高程值.

(23)

其中，havr表示當前分塊的平均高程值，n表示當前分塊DEM規則網格點的行列號，hi,j表示當前分塊DEM規則網格模型中第i行、第j列的點的高程值.

然后，搜索當前分塊的DEM規則網格模型中所有點的高程值，查找出max(hi,j)和min(hi,j)，將其與havr做差值，求得最大最小高度差：

h1=max(hi,j)-havr，
h2=min(hi,j)-havr，
|h1|≤hmax，
|h2|≤hmax,

(24)

其中，max(hi,j)和min(hi,j)是當前分塊DEM規則網格模型中點的高程值中的最大值和最小值，hmax是探測器安全降落參數所規定的最高或最低障礙物閾值.

4.2 坡度值

坡度表示地表地形的傾斜程度，如圖10所示：

Fig.10 Slope value圖10 坡度值

一般采用度數法來表示相應地形分塊的坡度值，即月表當前分塊與水平面的夾角，等同于月表分塊的法向量與z軸的夾角.

首先計算其法矢量ni,j:

(25)

然后計算其坡度值，即法向量ni,j與z軸的夾角φ.

(26)

4.3 粗糙度

粗糙度是反映地表凹凸程度的指標，本文采用平面夾角法來計算分塊的粗糙度，即將分塊的4個頂點對角線相連，確定4個平面，計算斜邊相鄰的2個平面的夾角，作為粗糙度參數.

設2個平面的夾角為θ，且θ∈[0,π]，則定義粗糙度R=1-cosθ∈[0,2].當2個平面的夾角θ=0°時，粗糙度R=0；當2個平面的θ=180°時，粗糙度R=2.這里計算2組平面的夾角粗糙度，取2組結果中較大的一個作為最終結果，提高了單一方法計算粗糙度的精度.具體計算流程如下.

1)計算構成4個平面的向量

例如平面ABC的構成向量為

IAB=(Xi,j+1-Xi,j,Yi,j+1-Yi,j,Zi,j+1-Zi,j)，
IBC=(Xi+1,j+1-Xi,j+1,Yi+1,j+1-Yi,j+1,
Zi+1,j+1-Zi,j+1).

(27)

如式(2)，分別計算出平面ABC，ADC，ABD，CBD的構成矢量.

2)計算相應平面的法向量

平面ABC法矢量n1和平面ADC的法矢量n2為

(28)

同理，計算出平面ABD和平面CBD的法矢量n3和n4.

3)計算法矢量的夾角余弦值

例如平面ABC的法矢量n1和平面ADC的法矢量n2的夾角的余弦值為

(29)

同理，計算出平面ABD和平面CBD的法矢量夾角的余弦值cosθ2.

4)計算2組地表粗糙度值

R1=1-cosθ1,
R2=1-cosθ2.

(30)

5)確定最終星球表面分塊粗糙度

根據2組粗糙度取較大值的原則，確定最終的月表分塊粗糙度.若R1

5 本文算法實驗結果

選取月面不同區域的地形測量數據作為輸入參數，對本文算法進行了實驗，驗證了本算法的正確性，結果如表1所示，在每組圖中從左至右依次為：遙感影像生成、自適應邊緣檢測圖、橢圓擬合障礙檢測圖、最大安全著陸區確定圖.

Table 1 Algorithm Verification Results表1 算法驗證結果

圖11為DEM規則網格三維重構數據，左側為散亂點經緯度及高程值，右側為算法執行后構建的DEM規則網格點的經緯度及高程值.

Fig.11 DEM 3D reconstructed data圖11 DEM三維重構數據

圖12為利用三維引擎所繪制出的三維重構后的高精度初選著陸區域DEM規則三角網格模型.

Fig.12 DEM 3D mesh model圖12 DEM三維網格模型

圖13為螺旋搜索安全判定策略所計算的各項參數，依次為當前分塊中心點經緯度、最大高度差、最小高度差、粗糙度和坡度值.

為了度量算法的準確性及算法效率，選取了月面不同區域的地形測量數據作為輸入參數，對算法進行了實驗，并與基于立體視覺和陰影分析的避障算法以及基于多因子著陸區分析評估方法進行了實驗分析對比.實驗結果如表2所示，結果表明本文算法在著陸區選擇準確性上高于基于多因子著陸區分析評估方法的精度，達到了基于立體視覺和陰影分析算法的實驗精度.在計算效率方面均高于基于多因子著陸區分析評估方法和基于立體視覺和陰影分析算法.

Fig.13 Spiral search security decision policy parameters圖13 螺旋搜索安全判定策略參數

Table 2 Efficiency and Accuracy Comparison of Algorithms表2 算法效率及準確性對比

6 結 論

探測器軟著陸障礙檢測與安全著陸是未來的深空探測的關鍵技術之一，本文提出了一種兩階段的探測器自主安全著陸仿真算法.首先給出了一種基于遙感影像生成的著陸區巖石障礙的檢測以及安全著陸區的確定方法;然后根據探測器軟著陸精避障的基本要求，提出了一種基于三維重構的探測器精避障螺旋掃描方法.對初選著陸區域進行高精度月面三維重構，螺旋搜索算法提高了目標著陸點搜索的時間效率25%以上，且尋找到最近的符合安全著陸閾值的目標著陸點，節約探測器能源.

本文的算法方案通過確定CCD相機拍攝的區域，根據高程數據對區域內的遙感影像圖片進行模擬，利用自適應的邊緣檢測方法和陰影檢測以及橢圓擬合技術對著陸區的障礙進行了檢測.提出了基于數學形態學的最大著陸圓算法，該算法快速有效，可用于實時障礙檢測，同時有效確定相對安全著陸區，為探測器的精確避障提供了良好的依據.該算法在探測器下降過程中，能夠對多幀圖像進行障礙檢測，從而保證著陸區選擇的有效性.基于三維重構的探測器精避障螺旋掃描方法.對初選著陸區域進行高精度月面三維重構提高了軟著陸的避障精度.本文算法準確性達到93%以上.經對比分析，本文算法在著陸區選擇準確性上高于基于多因子著陸區分析評估方法的精度，達到了基于立體視覺和陰影分析算法的實驗精度.在計算效率方面均高于基于多因子著陸區分析評估方法和基于立體視覺和陰影分析算法.