魏若巖，阮曉鋼，肖 堯，朱曉慶，3，黃 靜，4

（1.北京工業(yè)大學(xué)電子信息與控制工程學(xué)院，北京100124；2.河北經(jīng)貿(mào)大學(xué)信息技術(shù)學(xué)院，河北石家莊050061；3.麥吉爾大學(xué)機(jī)器人機(jī)械系統(tǒng)實(shí)驗(yàn)室，魁北克蒙特利爾H3A0C3；4.北京工業(yè)大學(xué)實(shí)驗(yàn)學(xué)院，北京101101）

基于單幅圖像且避封閉環(huán)境的星體表面著陸區(qū)選取方法

為了避免著陸器在著陸過(guò)程中落于地勢(shì)平坦的封閉環(huán)境，本文提出了一種基于單幅圖像的星表著陸區(qū)選擇方法。首先進(jìn)行障礙檢測(cè)，并且提取出圖像中非障礙區(qū)的骨架，然后以骨架三線交點(diǎn)為圓心得到多個(gè)具有允許著陸范圍的著陸圓，最后利用計(jì)分加權(quán)法從著陸圓的半徑、連通性、紋理，以及著陸圓心與障礙物的集群關(guān)系四個(gè)方面得到最佳著陸圓。通過(guò)與其他方法的對(duì)比實(shí)驗(yàn)發(fā)現(xiàn)，本文方法能選擇出障礙物稀少的區(qū)域，能夠避免著陸器落入地勢(shì)封閉的環(huán)境。

著陸區(qū)；星體表面；避障；封閉環(huán)境

魏若巖1，2，阮曉鋼1，肖 堯1，朱曉慶1，3，黃 靜1，4

（1.北京工業(yè)大學(xué)電子信息與控制工程學(xué)院，北京100124；2.河北經(jīng)貿(mào)大學(xué)信息技術(shù)學(xué)院，河北石家莊050061；3.麥吉爾大學(xué)機(jī)器人機(jī)械系統(tǒng)實(shí)驗(yàn)室，魁北克蒙特利爾H3A0C3；4.北京工業(yè)大學(xué)實(shí)驗(yàn)學(xué)院，北京101101）

0 引 言

在探索外天體的過(guò)程中，具有科學(xué)研究?jī)r(jià)值的區(qū)域往往位于一些特殊地形的周圍，這些地形包括隕石坑、山脊和丘陵等，并且復(fù)雜的地形會(huì)給著陸器帶來(lái)風(fēng)險(xiǎn)，再加上所探索的天體距離地球較遠(yuǎn)，地面控制中心無(wú)法對(duì)著陸器進(jìn)行實(shí)時(shí)控制，所以如何保證其自主選擇安全的著陸區(qū)是當(dāng)今深空探測(cè)中的一個(gè)重要課題。目前，關(guān)于星體表面安全軟著陸區(qū)的選擇方法主要有以下兩種［1-3］：主動(dòng)視覺(jué)和被動(dòng)視覺(jué)。

主動(dòng)視覺(jué)是通過(guò)安裝在著陸器的無(wú)線電雷達(dá)和激光雷達(dá)進(jìn)行三維地形數(shù)據(jù)的獲取，并且從中得到理想的著陸區(qū)域［1，49］，該方法已被應(yīng)用到了2012年成功著陸火星表面的“好奇號(hào)”，在此次任務(wù)中，火星著陸器在動(dòng)力下降段便開(kāi)始了無(wú)線電雷達(dá)對(duì)著陸地形的獲取，并且在著陸段實(shí)施避障。

被動(dòng)視覺(jué)是通過(guò)相機(jī)成像的方法對(duì)著陸區(qū)的地形進(jìn)行分析［1022］，基于被動(dòng)視覺(jué)的方法可以分為以下7個(gè)步驟［3］：

步驟1障礙檢測(cè)，該步驟主要分析障礙物的分布；

步驟2地形分析，該步驟包括著陸區(qū)域的紋理分析和坡度估計(jì)等；

步驟3在當(dāng)前探測(cè)器的燃料和動(dòng)力約束下對(duì)到達(dá)著陸區(qū)各個(gè)位置的約束性進(jìn)行分析；

步驟4在當(dāng)前探測(cè)器的位置，速度和姿下對(duì)到達(dá)著陸區(qū)域各個(gè)位置的精確度進(jìn)行分析；

步驟5基于以上4步的內(nèi)容進(jìn)行信息融合，得到整個(gè)著陸區(qū)域的著陸風(fēng)險(xiǎn)分布圖；

步驟6在著陸風(fēng)險(xiǎn)分布圖中搜索可能的著陸點(diǎn)；

步驟7選取最優(yōu)的著陸點(diǎn)。

1 問(wèn)題提出

基于被動(dòng)視覺(jué)方法已有多個(gè)文獻(xiàn)提出了自己的方案。針對(duì)被動(dòng)視覺(jué)步驟1中的障礙物檢測(cè)，文獻(xiàn)［4］提出了一種隕石坑和巖石檢測(cè)方法，該方法通過(guò)檢測(cè)圖像中的陰影區(qū)域預(yù)測(cè)隕石坑和巖石的位置以及尺寸，其中陰影區(qū)的檢測(cè)是基于最大熵［5］的方法求取陰影區(qū)域的灰度閾值，再利用K均值的方法將接近陰影區(qū)域像素值的像素點(diǎn)進(jìn)行聚類從而確定最終的陰影區(qū)域。該方法在文獻(xiàn)中的實(shí)驗(yàn)均取得了良好的效果，但是此方法沒(méi)有考慮針對(duì)強(qiáng)光照區(qū)域像素值的閾值求取，因?yàn)橹憛^(qū)域中的障礙物不但會(huì)形成陰影區(qū)，也會(huì)形成強(qiáng)光照區(qū)，雖然文獻(xiàn)［1］中通過(guò)陰影區(qū)域?qū)r石的位置進(jìn)行了預(yù)測(cè)，但是考慮到地形復(fù)雜的不確定性，此方法很難對(duì)障礙物進(jìn)行全面的檢測(cè)。文獻(xiàn)［6- 7］提出的障礙物（巖石）檢測(cè)方法也是基于陰影區(qū)域的檢測(cè)，該方法首先通過(guò)邊緣檢測(cè)檢測(cè)出陰影區(qū)域，然后根據(jù)陰影區(qū)域的形狀和光照方向?qū)r石的位置和形狀大小進(jìn)行估計(jì)，但是地形復(fù)雜的區(qū)域紋理較為復(fù)雜，對(duì)障礙物的提取會(huì)造成影響。文獻(xiàn)［8］提出了一種基于灰度變化標(biāo)準(zhǔn)差（intensity stander deviation，ISD）的障礙物檢測(cè)方法，該方法假設(shè)非障礙物區(qū)域具有高度的統(tǒng)一性，但是實(shí)際的著陸區(qū)域地形復(fù)雜，紋理較多，再加上圖像中存在一定的噪聲點(diǎn)使得障礙物檢測(cè)的準(zhǔn)確度下降，并且障礙物的提取需人工設(shè)定閾值，閾值的不同使得所提取到的障礙物的信息具有明顯的差異。針對(duì)步驟2中的地形分析主要有紋理分析和著陸區(qū)的坡度分析，文獻(xiàn)［9- 10］中利用ISD的紋理分析能夠有效的檢測(cè)出粗糙度較小的區(qū)域。針對(duì)坡度估計(jì)問(wèn)題，文獻(xiàn)［9］提出了一種基于序列圖像的著陸區(qū)坡度估計(jì)方法，并且文獻(xiàn)［10］在此基礎(chǔ)上進(jìn)行了改進(jìn)，從而減少了對(duì)探測(cè)器內(nèi)部的狀態(tài)測(cè)量?jī)x的依賴，但是基于序列圖像的三維地形重建往往需要較多的計(jì)算量，而且該算法對(duì)誤差較為敏感，所以很難實(shí)現(xiàn)地形的三維信息識(shí)別［1，11］。文獻(xiàn)［2，12］提出了一種基于陰影恢復(fù)形狀（shape from shade，SFS）的坡度估計(jì)方法，該方法通過(guò)地形的陰影和光照特點(diǎn)對(duì)地形進(jìn)行三維重建并且估計(jì)著陸區(qū)域的坡度，但是應(yīng)用該方法的前提是需要圖像中的景物具有光滑的表面［13］，大部分星體表面是粗糙的。文獻(xiàn)［14］提出了一種基于光流法的三維地形重建方法，但是該方法要求相機(jī)距離星體表面較近，當(dāng)距離較遠(yuǎn)時(shí)，三維重建的效果不明顯。針對(duì)步驟6中的著陸點(diǎn)搜索問(wèn)題文獻(xiàn)［15］提出了一種螺旋式的安全著陸區(qū)搜索方法，該方法已經(jīng)用在了2013年成功著陸在月球表面的“嫦娥3號(hào)”的粗避障階段［16］。針對(duì)步驟7中的最優(yōu)著陸點(diǎn)選取問(wèn)題，文獻(xiàn)［17］提出了一種基于形態(tài)學(xué)的最大著陸圓搜索方法。文獻(xiàn)［4，8］提出了一種基于模糊規(guī)則的最佳著陸區(qū)域選取方法，該方法將各個(gè)指標(biāo)分為好、中、差3個(gè)等級(jí)，但是指標(biāo)等級(jí)的劃分需要提前設(shè)定，這對(duì)于探索復(fù)雜的未知地形不具備靈活性。上述安全著陸點(diǎn)搜索方法均未考慮如何避免探測(cè)器落入地勢(shì)平坦而環(huán)境封閉的特殊地形，如圖1所示的白色圓圈區(qū)域，該圖像中隕石坑內(nèi)部的地勢(shì)較為平坦，無(wú)論從圖像紋理，障礙物分布以及面積大小都符合著陸器的著陸條件，但是該區(qū)域周邊被隕石坑的邊緣所包圍，若著陸器選擇該區(qū)域降落，巡視器（一種巡視車，在著陸器安全降落后釋放）將無(wú)法脫離該區(qū)域，所以巡視范圍將受到限制，嚴(yán)重的情況下將導(dǎo)致整個(gè)天體探測(cè)任務(wù)的失敗。

圖1 具有平坦地形的封閉環(huán)境（圖中白色圓圈區(qū)域）

本文針對(duì)被動(dòng)圖像方法所涉及的問(wèn)題在文獻(xiàn)［17］的基礎(chǔ)上提出了一種基于單幅圖像且避封閉環(huán)境的星體表面安全著陸區(qū)選擇方法，該方法在著陸器的粗壁障過(guò)程中進(jìn)行，不涉及三維地形信息，也不涉及探測(cè)器的燃料以及狀態(tài)的約束問(wèn)題。該方法主要有以下兩方面的創(chuàng)新：第一，針對(duì)被動(dòng)圖像方法第一步中的障礙檢測(cè)所涉及的問(wèn)題提出了一種雙閾值的陰影區(qū)和強(qiáng)光照區(qū)障礙物檢測(cè)算法，該算法根據(jù)圖像的直方圖對(duì)陰影區(qū)域閾值和強(qiáng)光照閾值進(jìn)行自適應(yīng)提?。坏诙?，針對(duì)被動(dòng)圖像方法中步驟6和步驟7中的避封閉區(qū)域問(wèn)題提出了一套新的方案，首先通過(guò)尋找非障礙物區(qū)域的骨架三線交點(diǎn)確定可能的著陸圓心，并且通過(guò)提出的外周長(zhǎng)擴(kuò)展法得到每個(gè)交點(diǎn)的最大著陸圓半徑。然后通過(guò)著陸圓半徑的大小，連通性，紋理特征以及每個(gè)著陸圓心與障礙物的集群關(guān)系確定最佳著陸圓。

本文在實(shí)驗(yàn)階段首先將提出的基于雙閾值二值化的障礙物檢測(cè)方法與文獻(xiàn)［1］和文獻(xiàn)［8］的相關(guān)方法進(jìn)行對(duì)比，從對(duì)比中發(fā)現(xiàn)本文提出的對(duì)于對(duì)大部分星體表面圖像能有效的提取出圖像中的障礙物，并且該算法不需任何參數(shù)的設(shè)置。然后將本文的安全著陸區(qū)的選取方法與文獻(xiàn)［1］和文獻(xiàn)［17］中的相關(guān)方法進(jìn)行對(duì)比，發(fā)現(xiàn)本文方法在避障階段不但能選出相對(duì)平坦和障礙物稀少的區(qū)域，并且能夠避免著陸器落入地勢(shì)封閉的區(qū)域。

2 障礙檢測(cè)

對(duì)障礙物背光與強(qiáng)光照區(qū)的提取完成障礙物檢測(cè)。首先對(duì)文獻(xiàn)［18］的陰影區(qū)域灰度閾值確定（determination of binary threshold，DBT）算法進(jìn)行改進(jìn)以適用陰影區(qū)域和強(qiáng)光照區(qū)域的提取，然后找出大于障礙物尺寸閾值的陰影區(qū)以及強(qiáng)光照區(qū)。

2.1 陰影區(qū)與亮區(qū)二值化閾值的自適應(yīng)提取

在DBT算法的基礎(chǔ)上，加入了強(qiáng)光照區(qū)域的提取，這里命名為陰影區(qū)和強(qiáng)光照區(qū)灰度閾值確定法（determination of binary threshold for shadow and strong illumination areas，DBTSSI），具體步驟如下。

步驟1計(jì)算平均灰度值：

式中，M和N分別為像素點(diǎn)在該圖像中的行列數(shù)目；f（x，y）是位于圖像坐標(biāo)（x，y）的像素值。

步驟2得到圖像中小于f－的像素值向量b1和不小于f－的像素值向量b2，并且對(duì)b1和b2進(jìn)行升序和降序的排列：

式中，i＝1，2，…，n1；j＝1，2，…，n2。n1和n2分別為向量b1和b2中像素值的數(shù)量。

步驟3計(jì)算b1和b2與f－的差值向量b′1和b′2：

式中，i＝1，2，…，n1，j＝1，2，…，n2。

步驟4計(jì)算b′1和b′2的變化率，于是得到兩組新的向量v1和v2，并且計(jì)算其平均值v－1和v－2：

式中，i＝1，2，…，n1－1；j＝1，2，…，n2－1。

式中，i＝1，2，…，n1－1；j＝1，2，…，n2－1。所以，L1和L2即為v1和v2的高變化區(qū)間長(zhǎng)度。

步驟6計(jì)算L1和L2在v1和v2中所占的長(zhǎng)度比例p1和p2：

步驟7在b1和b2中計(jì)算高變化區(qū)域的標(biāo)準(zhǔn)差f′1和f′2：

步驟8得到陰影閾值α1和強(qiáng)光照閾值α2：

圖2展示了星表圖像的二值化效果，從中看出利用DBTSSI算法可以很好地提取到圖像中陰影區(qū)域和強(qiáng)光照區(qū)域。

圖2 二值化效果圖

2.2 障礙物尺寸閾值的設(shè)定

考慮到可能存在的噪聲以及微小障礙物，需要設(shè)定一個(gè)障礙物尺寸的閾值，然后將小于該尺寸閾值的障礙物去除。尺寸閾值通常參照相機(jī)的焦距，著陸器相機(jī)距離地面的垂直距離，以及著陸器的尺寸進(jìn)行設(shè)定，閾值公式可參考式（18）：

式中，Tarea（pixel2）為尺寸閾值；f為光學(xué)相機(jī)的焦距，單位：m；d為著陸器相機(jī)距離地面的垂直距離，單位：km；若著陸器為一個(gè)長(zhǎng)方體，且L和S分別為著陸器的長(zhǎng)和寬，單位：m；t為像素點(diǎn)數(shù)量和米的轉(zhuǎn)換關(guān)系，即一米長(zhǎng)度有多少個(gè)像素點(diǎn)，單位：pixels／m；λ為著陸器的膨脹系數(shù)，即障礙物的大小尺寸不能大于λLS，單位：m2，在本文中λ符合下面的關(guān)系式：

如果著陸器相機(jī)距離地面的垂直距離d為1 km，相機(jī)焦距f為2 cm，著陸器的長(zhǎng)寬L和S均為3 m，λ為10，并且相機(jī)中像素點(diǎn)的大小為54μm，即t＝1 000×1 000／54＝18 519，根據(jù)式（18）得到Tarea約為124，則圖2（a）經(jīng)過(guò)尺寸閾值篩選后得到的障礙物二值圖像如圖3所示。

圖3 經(jīng)尺寸閾值篩選后得到的障礙物二值圖像

3 備選著陸圓的提取

3.1 非障礙區(qū)域的骨架提取

骨架是表示平面區(qū)域結(jié)構(gòu)形狀的一種重要的方法，并且該方法是通過(guò)一種細(xì)化算法實(shí)現(xiàn)，關(guān)于骨架的算法細(xì)節(jié)可參考文獻(xiàn)［19］。圖4展示了非障礙物區(qū)域的骨架。

圖4 圖2（a）的非障礙物區(qū)域骨架

3.2 骨架交點(diǎn)的提取

若將每個(gè)障礙物看做一個(gè)質(zhì)點(diǎn)，根據(jù)三角形的三邊中垂線交點(diǎn)為該三角形的外接圓心的定理確定具有較大半徑著陸圓的圓心一般位于骨架的三線交點(diǎn)處。為了提取到骨架的交點(diǎn)，采用卷積的形式，如果圖像用函數(shù)f（x，y）表示，即圖像在（x，y）的像素值為f（x，y），模板函數(shù)用δ（x，y）表示，則圖像與模板的卷積表示為

式中，m和n為模板函數(shù)的尺寸。

為了提取骨架三線交點(diǎn)，令δ為2×2的模板，如圖5（a）所示，當(dāng)δ位于三線交點(diǎn)時(shí)，有圖5（b）的4種形式，并且交點(diǎn)灰度值為3×256／4＝192，由于兩線交點(diǎn)的灰度值為2× 256／4＝128，所以灰度閾值在［128，192］范圍內(nèi)。如果設(shè)定閾值為190，圖6中的白色像素點(diǎn)為骨架三線交點(diǎn)。

圖5 卷積模板及其三線交點(diǎn)的4種形式

圖6 提取到的骨架三線交點(diǎn)

3.3 基于外周長(zhǎng)擴(kuò)展的安全半徑確定

本文提出了一種基于外周長(zhǎng)擴(kuò)展的最大半徑搜索方法（perimeter extend，PE），該方法類似螺旋式方法［15，20-21］，原理如圖7所示，步驟如下。

圖7 外周長(zhǎng)擴(kuò)展法示意圖

步驟1確定安全著陸區(qū)所要求的最小正方體區(qū)域。假設(shè)著陸器為一個(gè)長(zhǎng)方體，L和S分別為著陸器的長(zhǎng)和寬，并且L≥S，考慮到一系列誤差源的存在，如：導(dǎo)航誤差、通信誤差等，必須根據(jù)著陸器的尺寸設(shè)計(jì)一個(gè)可允許著陸的最小正方形區(qū)域，該區(qū)域要大于著陸器的尺寸，所以需設(shè)定一個(gè)膨脹參數(shù)，該參數(shù)能夠保證實(shí)際著陸區(qū)尺寸與著陸器尺寸之間存在一個(gè)比例關(guān)系，則該最小正方形著陸區(qū)域的邊長(zhǎng)為

式中，L′為最小正方形的邊長(zhǎng)，單位：pixel；f為相機(jī)的焦距單位：m；d為相機(jī)距離行星表面的垂直距離，單位：km；t為像素點(diǎn)數(shù)量和米的轉(zhuǎn)換關(guān)系，單位：pixels／m；λ為著陸器的膨脹系數(shù)，并且符合關(guān)系式（19）。

步驟2以骨架三線交點(diǎn)為圓心，每次將正方形邊長(zhǎng)向相反方向擴(kuò)展L′個(gè)長(zhǎng)度。假設(shè)O為障礙物像素點(diǎn)集合，P為每次擴(kuò)展后的外周長(zhǎng)像素點(diǎn)集合，如果P∩O＝?，外周長(zhǎng)繼續(xù)擴(kuò)展，否則停止擴(kuò)展。

步驟3計(jì)算骨架三線交點(diǎn)到外周長(zhǎng)的垂直距離，該距離長(zhǎng)度即為以該交點(diǎn)為圓心的安全著陸圓的半徑。

3.4 備選著陸圓的確定

考慮到在導(dǎo)航和控制過(guò)程中可能存在的誤差，著陸位置要距離障礙物較遠(yuǎn)，如果著陸圓心視為著陸器的著陸點(diǎn)，則備選著陸圓的選擇需考慮著陸圓的半徑大小，所以下面介紹一種備選著陸圓的選擇方法，該方法分為以下兩個(gè)步驟。

步驟1對(duì)所有圓按照半徑的大小進(jìn)行升序排列，并且選出M個(gè)著陸圓，該M個(gè)著陸圓須可由下面的算法得出。

步驟2將M個(gè)著陸圓進(jìn)行降序排列得到新的組合M′，令T為最終備選著陸圓數(shù)量，T符合式（22）的條件：

式中，ri為第i（i＝1，2，…，M）個(gè)著陸圓的半徑。

圖8（a）為經(jīng)過(guò)以上兩個(gè)步驟提取到的37個(gè)半徑由大到小的備選著陸圓，從中可以看出，這些著陸圓相比較障礙物區(qū)域具有較為平整的地形；圖8（b）為備選著陸圓半徑歸一化后的排列關(guān)系，圖中，縱坐標(biāo)ncr表示歸一化半徑。

圖8 備選著陸圓及半徑

4 最佳著陸圓的選取

最佳著陸圓的選取需要從4個(gè)方面進(jìn)行考慮：著陸圓的半徑、著陸圓的連通性、著陸圓的紋理以及著陸圓心與障礙物的集群關(guān)系。由于著陸圓的半徑已在上一節(jié)進(jìn)行了分析，下面將對(duì)另外3個(gè)指標(biāo)進(jìn)行討論。

4.1 著陸圓的連通性

考慮到著陸器著陸后釋放的巡視器要有較大的巡視范圍，所以著陸器著陸的著陸圓附近要有較少的障礙物，并且著陸圓之間要有較好的連通性。著陸圓的連通性在本文定義為：假設(shè)存在3個(gè)著陸圓A、B、C，并且O1、O2、O3分別為3個(gè)著陸圓心。如果A∩B≠?并且A∩C＝?和B∩C＝?，如圖9（a）所示，則A與B相連通，C與A和B均不具備連通性；如果A∩B≠?并且B∩C≠?，但是A∩C＝?，如圖9（b）所示，則A與C相連通。在圖9（a）中的AB以及圖9（b）中的ABC均稱為連通區(qū)域。單獨(dú)的圓可認(rèn)為與自身連通，如圖9（a）中的圓C。

圖9 著陸圓的連通關(guān)系示意圖

本文基于廣度優(yōu)先搜索遍歷算法提出了一種連通區(qū)域的查找方法，該方法可分為以下3個(gè)步驟。

步驟1假設(shè)有T個(gè)備選著陸圓，對(duì)該備選著陸圓以半徑大小進(jìn)行升序排列，于是得到備選著陸圓的半徑排列向量R和圓心位置排列向量P。

步驟2建立一個(gè)T×T矩陣W，該矩陣記錄了當(dāng)前備選著陸圓的連接關(guān)系，連接關(guān)系符合下面的式子：

式中，i，j＝1，2，…，T。

步驟3基于下面的算法對(duì)連通區(qū)域進(jìn)行查找。____

圖10展示了圖8（a）中連通區(qū)域查找結(jié)果，連通區(qū)域用黃色線段對(duì)其圓心進(jìn)行了連接。

圖10 連通區(qū)域查找結(jié)果

假設(shè)存在T個(gè)備選著陸圓，且有Q個(gè)連通區(qū)域，則每個(gè)連通區(qū)域的連通性表示為

式中，acj為第j個(gè)連通區(qū)域的連通性；Kj表示第j個(gè)連通區(qū)域有Kj個(gè)著陸圓；area為第j個(gè)連通區(qū)域中的第k個(gè)著陸圓。本文規(guī)定每個(gè)著陸圓的連通性等同于所在連通區(qū)域的連通性，即：cci＝acj且cci?acj，其中cci為第i個(gè)著陸圓的連通性，i＝1，2，…，M，j＝1，2，…，Q。若將著陸圓的連通性進(jìn)行歸一化，可以表示為

式中，ncci表示第i個(gè)著陸圓歸一化后的連通性。

圖11為著陸圓的歸一化連通性，從中看出在相同連通區(qū)域的著陸圓具有相同的連通性。

圖11 著陸圓的歸一化連通性

4.2 著陸圓的紋理

本文第1節(jié)介紹了障礙檢測(cè)，但是這并不意味著沒(méi)有障礙物的區(qū)域一定是適合著陸器著陸的，因?yàn)橐恍┎黄教沟膮^(qū)域無(wú)法通過(guò)障礙檢測(cè)的方法進(jìn)行確定，所以本文介紹一種通過(guò)檢測(cè)著陸圓紋理的方法確定其平坦性。

利用文獻(xiàn)［19］提出的一致性度量對(duì)著陸圓的紋理特征進(jìn)行統(tǒng)計(jì)，該方法可以用式（26）表示：

式中，u代表一致性度量；z為代表灰度級(jí)的變量；p（zi）為第i個(gè)灰度級(jí)zi的像素點(diǎn)數(shù)量在整個(gè)圖像直方圖中所占的比例，i＝1，2，…，L，L為可區(qū)分的灰度級(jí)數(shù)目。

假設(shè)存在T個(gè)備選著陸圓，若將著陸圓一致性度量歸一化，可表示為

式中，nct表示著陸圓歸一化后的一致性度量，ui為序號(hào)為i的著陸圓。圖12為圖8（a）中著陸圓歸一化一致性度量，從中看出，一致化度量高的著陸圓在圖像中紋理單一，地形平坦，而一致化度量較低的著陸圓在圖像中的紋理較為復(fù)雜。

圖12 著陸圓的歸一化一致性度量

4.3 著陸圓心的障礙物集群性

本節(jié)提出了一種表示障礙物集群性方法，該方法能表示圖像中的任意位置與障礙物的位置關(guān)系。假設(shè)計(jì)算圖像中某一點(diǎn)P與圖像障礙物的位置關(guān)系wp，如果圖像中存在N個(gè)障礙物體，并且點(diǎn)P到這N個(gè)障礙物的重心的歐式距離表示為：di（i＝1，2，…，N），當(dāng)di≤10個(gè)像素時(shí)，di＝10；障礙物的大小表示為：wi（i＝1，2，…，N），即為每個(gè)障礙物所占的像素點(diǎn)的數(shù)量，點(diǎn)P與圖像障礙物的位置關(guān)系表示為

圖13展示了圖2（a）中每一個(gè)坐標(biāo)點(diǎn)的wp值，其中圖像的橫縱坐標(biāo)分別表示圖像的橫縱尺寸，從中看出，藍(lán)色深的區(qū)域具有較小的wp值，并且障礙物分布越稀少，距離圖像中其他的障礙物較遠(yuǎn)，反之，淺藍(lán)色和非藍(lán)色區(qū)域有較大的wp值，障礙物分布較多，距離圖像障礙物較近。

圖13 wp值高程圖

下面分析備選著陸圓心的wp值，假設(shè)存在T個(gè)備選著陸圓，若將著陸圓圓心的wp值歸一化，可表示為

式中，ncwi為第i個(gè)著陸圓心歸一化的wp值。

圖14展示了圖8（a）中著陸圓心的歸一化wp值，結(jié)合圖13可以看出具有較大的ncw值的著陸圓心基本分布在藍(lán)色較淺的區(qū)域，而較小的ncw值的著陸圓心分布在藍(lán)色較深的區(qū)域。

圖14 著陸圓心的歸一化wp值

4.4 最佳著陸圓的提取

本文根據(jù)著陸圓的半徑大小，連通性，紋理以及著陸圓心的障礙物集群性提出了一種最佳著陸圓的選區(qū)方法：指標(biāo)計(jì)分加權(quán)法，步驟如下。

步驟2根據(jù)指標(biāo)的平均值，計(jì)算每個(gè)備選著陸圓的分?jǐn)?shù)。

所以備選著陸圓的分?jǐn)?shù)為指標(biāo)的分?jǐn)?shù)加和。

步驟3輸出最佳著陸圓。

若N′為分?jǐn)?shù)最高的備選著陸圓的數(shù)量，若N′＝1，輸出當(dāng)前著陸圓；若N′＞1，設(shè)定綜合指標(biāo)ind，其中，i＝1，2，…，N′，j＝1，2，…，n，n為指標(biāo)項(xiàng)數(shù)，令ind＝ncri，ind＝，然后計(jì)算綜合指標(biāo)的加和，即最后以具有最大綜合指標(biāo)和的著陸圓為最佳著陸圓。

圖15為最佳著陸圓選取示意圖，其中圖15（a）為步驟2中計(jì)分最高的著陸圓，著陸圓的序號(hào)分別為：2、4、7，根據(jù)圖10，2號(hào)和7號(hào)著陸圓在同一個(gè)連通區(qū)域，而4號(hào)著陸圓在另外的連通區(qū)域，并且2號(hào)和7號(hào)著陸圓所在的連通區(qū)域面積大于4號(hào)著陸圓所在的連通區(qū)域；觀察圖11，2號(hào)著陸圓的一致性度量在3個(gè)著陸圓中處于最高值；觀察圖14，4號(hào)著陸圓的wp值最高，2號(hào)和7號(hào)著陸圓的wp值基本持平，這說(shuō)明4號(hào)著陸圓的地面特征在三者中最為復(fù)雜；由于2號(hào)著陸圓的半徑在三者中最大，所以結(jié)合步驟3得出2號(hào)著陸圓為最佳著陸圓，如圖15（b）所示。

圖15 最佳著陸圓選取示圖

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

5.1 針對(duì)障礙檢測(cè)方法的對(duì)比

圖16給出了文獻(xiàn)［5］中3幅星體表面圖像，其中圖16（a）和圖16（b）兩幅圖像為星體表面的模擬圖像，從中可以看出由于光照角度的不同，障礙物（石塊）的陰影區(qū)的大小不同。圖16（c）［23］為433 Eros的星體表面圖像，從中看出石塊在光照下形成較為明顯的陰影區(qū)。下面實(shí)驗(yàn)用本文的障礙物檢測(cè)方法與文獻(xiàn)［1］和文獻(xiàn)［8］的相關(guān)方法進(jìn)行對(duì)比。

圖16 3幅星體表面地形圖像

圖17給出了3個(gè)算法的對(duì)比結(jié)果，它們的運(yùn)行時(shí)間分別為：文獻(xiàn)［1］：1.5 s，文獻(xiàn)［8］：2.3 s，本文算法：1.4 s。從效果上看，文獻(xiàn)［1］可以根據(jù)當(dāng)前陰影區(qū)的大小和形狀能估計(jì)出巖石的位置和尺寸（圖中的橢圓為估計(jì)出的巖石，橢圓的方向，大小和形狀體現(xiàn)了估計(jì)出的巖石的放置位置，形狀以及大?。?，但是障礙物尺寸的估計(jì)往往受光照角度的影響較大，圖16（a）和圖16（b）的光照角度存在明顯的不同，利用文獻(xiàn)［1］的方法估計(jì)出的石塊的大小和尺寸存在明顯的差異。文獻(xiàn)［8］對(duì)于障礙物的檢測(cè)具有較好的結(jié)果，但是檢測(cè)出的障礙物圖像為非二值化圖像，還需設(shè)定閾值對(duì)障礙物的信息進(jìn)行提取，但是障礙物的提取效果對(duì)于閾值的設(shè)定較為敏感。雖然本文算法檢測(cè)出的障礙物存在較多的噪聲，但是和文獻(xiàn)［1］以及文獻(xiàn)［8］的方法相比巖石的輪廓效果較好，在不同光照角度下檢測(cè)到的強(qiáng)光照區(qū)域具有相似的位置和尺寸，然而利用文獻(xiàn)［1］，巖石的位置和尺寸在不同的光照條件下的估計(jì)差異較大。

圖17 3種算法的對(duì)比結(jié)果

5.2 安全著陸區(qū)的選取對(duì)比

針對(duì)安全著陸區(qū)的選取本節(jié)進(jìn)行兩組對(duì)比實(shí)驗(yàn)，第一組的對(duì)比對(duì)象為文獻(xiàn)［1］的相關(guān)方法，由于文獻(xiàn)［1］對(duì)備選著陸區(qū)域的選取方法在細(xì)節(jié)上沒(méi)有詳細(xì)的介紹，所以本文只能針對(duì)文獻(xiàn)［1］中用到的圖像進(jìn)行實(shí)驗(yàn)，然后將實(shí)驗(yàn)結(jié)果與文獻(xiàn)［1］中的進(jìn)行對(duì)比。由于本文方法是基于文獻(xiàn)［17］的改進(jìn)，并且為了說(shuō)明避封閉區(qū)域的必要性，第二組實(shí)驗(yàn)與文獻(xiàn)［17］的相關(guān)方法進(jìn)行對(duì)比。

5.2.1 第一組對(duì)比實(shí)驗(yàn)

圖18給出了文獻(xiàn)［1］的一幅圖像，圖像中的巖石在一定的光照條件下文獻(xiàn)的實(shí)驗(yàn)部分用提出的方法對(duì)安全著陸區(qū)域進(jìn)行了提取，如圖19（a）所示，從中看出文獻(xiàn)［1］的方法選出了多個(gè)備選著陸區(qū)。圖19（b）為應(yīng)用本文方法的到的備選著陸圓，一共有38個(gè)半徑不同的著陸圓被選出，并且這些備選著陸圓包含了大部分文獻(xiàn)［1］得到的著陸區(qū)域。圖19（c）給出了利用本文方法得到的最終著陸圓，從中可以看出最終著陸圓相比其他著陸圓具有較大的半徑，內(nèi)部的紋理相對(duì)簡(jiǎn)單，著陸圓的圓心距離周圍的障礙物較遠(yuǎn)，并且該著陸圓與其他備選著陸圓具有較好的連通性，所以可以作為理想的著陸區(qū)域。

圖18 文獻(xiàn)［1］中的一幅圖像

圖19 本文方法與文獻(xiàn)［1］方法在著陸區(qū)域選擇上的對(duì)比

5.2.2 第二組對(duì)比實(shí)驗(yàn)

圖20展示了4種不同的星表地形，圖20（a）（http：∥nssdc.gsfc.nasa.gov／planetary／mission／near／descent_images／near_descent_157416593.html）為Eros433的地表圖像，該區(qū)域分布有一些大小不等的石塊，并且部分非障礙區(qū)域的表面較為粗糙。圖20（b）為火星某區(qū)域的地表圖，該區(qū)域在火星上的坐標(biāo)為：21°38′7.52″N，173°36′22.45″W，該區(qū)域隕石坑內(nèi)部的地形較為平坦，并且隕石坑的邊緣部分形成了明顯的環(huán)形山，隕石坑外部分布有地形平坦的區(qū)域。圖20（c）也為火星某區(qū)域地表圖，該區(qū)域在火星的坐標(biāo)為：24°12′6.30″N，176°35′18.90″W，在圖中有兩個(gè)明顯的隕石坑，并且隕石坑的特征明顯不同。左上方的隕石坑所形成的環(huán)形山較為明顯，并且隕石坑內(nèi)部地勢(shì)凹凸明顯，右下方的隕石坑所形成的環(huán)形山相比第一個(gè)較為模糊，并且隕石坑周邊形成了復(fù)雜多變的地形，但是隕石坑內(nèi)部卻有著非常平坦的地形。圖20（d）為月球某區(qū)域圖像（http：∥lroc.sese.asu.edu／），該區(qū)域較之前3幅，其地形最為復(fù)雜，并且分布有許多大小不等的隕石坑，隕石坑內(nèi)部和周邊均較少分布有地形平坦的區(qū)域?，F(xiàn)在假設(shè)在這四幅圖像中距離地面的高度均為2 km，相機(jī)的焦距為5 cm，并且著陸器的長(zhǎng)寬L和S均為3 m，λ為0.1，λ為著陸器的膨脹系數(shù)，在本實(shí)驗(yàn)中所有的著陸區(qū)域均為算法的自動(dòng)選取，未考慮特定的科學(xué)興趣區(qū)域。實(shí)驗(yàn)環(huán)境為Matlab2013a，Win7操作系統(tǒng)，1G內(nèi)存，Core i3處理器。

圖20 4種不同星體的地形圖像

圖21為文獻(xiàn)［17］和本文方法結(jié)果的對(duì)比，在文獻(xiàn)［17］中，所有的最終著陸圓為探測(cè)到的半徑為最大的圓，在20（a）中，文獻(xiàn)［17］的方法所得到的圓內(nèi)部的紋理相比本文方法要復(fù)雜，并且周邊分布有較多的障礙物石塊，雖然本文得到的著陸圓的右上方有一個(gè)大的石塊，但是所選著陸圓的紋理較為單一，地形較為平坦。所以本文算法結(jié)果針對(duì)圖20（a）得到了更為合適的著陸區(qū)。在圖20（b）中，文獻(xiàn)［17］得到的著陸圓位于隕石坑的內(nèi)部，雖然該區(qū)域范圍較大，并且區(qū)域的紋理特征單一，地形平坦，但是該區(qū)域在隕石坑的內(nèi)部，如果著陸器降落于此會(huì)導(dǎo)致釋放的巡視器無(wú)法走出隕石坑，從而造成巡視器的探測(cè)范圍受限，無(wú)法向其他區(qū)域進(jìn)行探索。本文算法得到的著陸圓則選擇了隕石坑的周邊區(qū)域，雖然該著陸圓的半徑與文獻(xiàn)［17］得到的最大著陸圓相比較小，但是可以使得巡視器有更加廣闊的巡視范圍，本文算法之所以未把著陸圓選擇在隕石坑內(nèi)是因?yàn)殡E石坑的邊界布滿了突起的環(huán)形山，這些環(huán)形山所形成的障礙物使得隕石坑內(nèi)部的wp值要明顯大于隕石坑的周邊地區(qū)，所以具有較高wp值的著陸圓心不易在最佳著陸圓選取的過(guò)程中勝出。在圖20（c）中，著陸圓選擇在了邊緣較為模糊的隕石坑內(nèi)，雖然在該區(qū)域中巡視器能通過(guò)搜尋路徑繞出隕石坑，但是本文算法得到的著陸圓區(qū)域則避免了搜尋過(guò)程。在圖20（d）中，本文方法和文獻(xiàn)［17］方法均選擇了相同的區(qū)域，即該區(qū)域地形較為平坦，并且均躲避了隕石坑。所以從該實(shí)驗(yàn)中可以看出，本文方法不但能選擇地勢(shì)平坦的區(qū)域，也考慮了接下來(lái)的巡視器的巡航范圍，使其能著陸在在障礙物較為稀少，并且避免降落在地勢(shì)平坦但是環(huán)境封閉的區(qū)域。

圖21 文獻(xiàn)［17］方法與本文方法的對(duì)比

針對(duì)兩種方法在運(yùn)行時(shí)間的上的比較如表1所示，從中可以發(fā)現(xiàn)本文方法所需的時(shí)間要明顯多于文獻(xiàn)［17］的方法，這是因?yàn)楸疚牡姆椒ㄊ窃谖墨I(xiàn)［17］的基礎(chǔ)上進(jìn)行提出的，更確切的說(shuō)是文獻(xiàn)［17］方法的一種再加工，多出了著陸圓之間的連通性分析，著陸圓的紋理分析以及最終著陸圓的選擇上的運(yùn)行時(shí)間，所以本文方法雖然在避封閉環(huán)境的功能上具有一定的優(yōu)勢(shì)，但是方法的時(shí)效性還有待于提高。

表1 本文方法和文獻(xiàn)［17］的運(yùn)行時(shí)間的比較 s

6 結(jié) 論

本文在文獻(xiàn)［17］的基礎(chǔ)上提出了一種新的安全著陸區(qū)的搜索方法，該方法首先在障礙物檢測(cè)中運(yùn)用了本文提出的DBTSSI算法，在著陸圓的選擇方面考慮了圓半徑，著陸圓的連通關(guān)系，著陸圓的紋理以及著陸圓心得障礙物集群性4方面的決策因素，并且對(duì)于最佳著陸圓的選擇提出了一種加分加權(quán)法。最后通過(guò)對(duì)比實(shí)驗(yàn)表明，本文算法不但能為著陸器選擇出地勢(shì)平坦和范圍較大的區(qū)域，也能使其落入障礙物較為稀少，并且避免其落入地勢(shì)封閉的環(huán)境，為接下來(lái)的巡視器的探索帶來(lái)更廣闊的巡航空間。但是，本文所提方法的時(shí)效性目前還不能滿足于實(shí)際的工程應(yīng)用中，所以如何提高方法的時(shí)效性將是未來(lái)研究工作的重點(diǎn)。

［1］Huertas A，Cheng Y，Madison R.Passive imaging based multicue hazard detection for spacecraft safe landing［C］∥Proc.of the IEEE Aerospace Conference，2006.

［2］Brady T，Robertson E，Epp C.Hazard detection methods for lunar landing［C］∥Proc.of the IEEE Aerospace Conference，2009.

［3］Neveu D，Mercier G，Hamel J F，et al.Passive versus active hazard detection and avoidance systems［C］∥Proc.of the 9th International ESA Conference on Guidance，Navigation&Control Systems，2014：1- 8.

［4］Serrano N，Seraji H.Landing site selection using fuzzy rulebased reasoning［C］∥Proc.of the IEEE International Conference on Robotics and Automation，2007：4899- 4904.

［5］Kapur J N，Sahoo P K，Wong A K C.A new method for gray-level picture thresholding using the entropy of the histogram［J］.Computer Vision，Graphics，and Image Processing，1985，29（3）：273- 285.

［6］Crane E S，Rock S M.Guidance augmentation for reducing uncertainty in vision-based hazard mapping during lunar landing［C］∥Proc.of the Aerospace Conference，2013：1- 12.

［7］Crane E S，Rock S M.Influence of trajectory on accuracy of hazard estimation during lunar landing［C］∥Proc.of the AIAA Guidance，Navigation，and Control Conference，2012：1- 15.

［8］Feng J，Cui P Y，Cui H T.Autonomous hazard detection and land-ing point selecting for planetary landing［C］∥Proc.of the Systems and Control in Aeronautics and Astronautics，2010：1292- 1296.［9］Cheng Y，Johnson A E，Mattheis L H.Passive imaging based hazard avoidance for spacecraft safe landing［C］∥Proc.of the Aerospace Conference，2006：1- 8.

［10］Ding M，Cao Y F，Wu Q X.Passive image-based safe landing site searching and slope estimation in probe landing［J］.Journal of Astronautics，2009，30（6）：2258- 2263.（丁萌，曹云峰，吳慶憲.月球著陸器著陸中基于被動(dòng)圖像的著陸場(chǎng)搜索及坡度估計(jì)［J］.宇航學(xué)報(bào)，2009，30（6）：2258- 2263.）

［11］Mahmood W，Shah S M A.Vision based hazard detection and obstacle Avoidance for planetary landing［C］∥Proc.of the Nonlinear Dynamics and Synchronization，2009：175- 181.

［12］Rogata P，Di S E，Camara F，et al.Design and performance assessment of hazard avoidance techniques for vision-based landing［J］.Acta Astronautica，2007，61（1）：63- 77.

［13］Horn B K P.Height and gradient from shading［J］.International Journal of Computer Vision，1990，5（1）：37- 75.

［14］Janschek K，Tchernykh V，Beck M.Performance analysis for visual planetary landing navigation using optical flow and dem matching［C］∥Proc.of the AIAA GNC Conference，2006：1- 17.

［15］Zhang Z X，Cui P Y.An algorithm of rock detection and avoidance for planetary soft landing based on CCD landing camera［J］.Journal of Astronautics，2008，29（6）：1510- 1516.（張澤旭，崔平遠(yuǎn).基于CCD著陸相機(jī)的行星軟著陸巖石檢測(cè)與規(guī)避方法［J］.宇航學(xué)報(bào)，2008，29（6）：1510- 1516.）

［16］Huang X Y，Zhang H H，Wang D Y.Autonomous navigation and guidance for Chang’e3 soft landing［J］.Journal of Deep Space Exploration，2014，1（1）：52- 59.（黃翔宇，張洪華，王大軼.“嫦娥三號(hào)”著陸器軟著陸自主導(dǎo)航與制導(dǎo)技術(shù)［J］.深空探測(cè)學(xué)報(bào)，2014，1（1）：52- 59.）

［17］Shao W，Cui P Y，Cui H T.An algotithm of site detection for spacecraft safe landing based on morphology［J］.Journal of Astronautics，2008，29（2）：659- 664.（邵巍，崔平遠(yuǎn)，崔祜濤.一種基于形態(tài)學(xué)的著陸器安全著陸點(diǎn)選擇方法［J］.宇航學(xué)報(bào)，2008，29（2）：659- 664.）

［18］Ruan X G，Wei R Y，Li J G.Method of tracking for soft landing on small body based on shadow areas［J］.Control and Decision，2014，29（9）：1581- 1586.（阮曉鋼，魏若巖，李建更.一種小天體軟著陸中基于地表陰影區(qū)的跟蹤算法［J］.控制與決策，2014，29（9）：1581- 1586.）

［19］Gonzalez R C，Woods R E.Digital image processing［M］.New Jersey：Pearson Education，2002.

［20］Geng D W，Song Z J，Jiang X W.Scanning arithmetic for softlanding based on time optimization［J］.Chinese Space Science and Technology，2006，26（4）：6- 10.（耿殿伍，宋政吉，姜興渭.月球軟著陸的時(shí)間最優(yōu)安全著陸區(qū)搜索算法［J］.中國(guó)空間科學(xué)技術(shù)，2006，26（4）：6- 10.）

［21］Wang D Y，Li J，Huang X Y.A point autonomous navigation and hazard avoidance method for lunar soft landing［J］.Journal of Deep Space Exploration，2014，1（1），51- 54.（王大軼，李驥，黃翔宇.月球軟著陸過(guò)程高精度自主導(dǎo)航避障方法［J］.深空探測(cè)學(xué)報(bào)，2014，1（1），51- 54.）

［22］Scott M.NEAR image of the day for 2001 Feb 12（E）［EB／OL］.［2014- 12- 31］.http：∥near.jhuapl.edu／iod／20010212e／index.html.

Planetary surface landing area selection based on single image and closed environment avoidance

WEI Ruo-yan1，RUAN Xiao-gang1，XIAO Yao1，ZHU Xiao-qing1，3，HUANG Jing1，4
（1.School of Electronic Information and Control Engineering，Beijing University of Technology，Beijing 100124，China；2.College of Information Technology，Hebei University of Economics and Business，Shijiazhuang 050061，China；3.Center for Intelligent Machine，McGill University，Montreal H3A0C3，Canada；4.The Pilot College，Beijing University of Technology，Beijing 101101，China）

For avoiding the lander landing on the flat areas in a closed environment in the process of soft landing on asteroid，a landing area selection method based on a single image is proposed.Firstly，the proposed method of binary is applied to extract hazards，and the skeleton of areas with no hazards is extracted.Secondly，take the intersection points with three lines of the skeleton as the center of landing circles.Finally，the best landing circle can be obtained by the scoring method with the help of radius，connection among circles，circle’s texture characters，and the cluster feature of circle’s center with all the hazards.Compared with other methods the proposed landing area selection method is highlighted with fewer obstacles，and can avoide landing on the place in a closed environment.

landing area；surface of asteroid；hazards avoidance；closed environment

TP 391

A

10.3969／j.issn.1001-506X.2015.12.21

魏若巖（198-4- ），男，講師，博士，主要研究方向?yàn)槟Ｊ阶R(shí)別與智能系統(tǒng)。

E-mail：weiruoyan1984＠163.com

阮曉鋼（1958- ），男，教授，博士研究生導(dǎo)師，博士，主要研究方向?yàn)榭刂瓶茖W(xué)、模式識(shí)別。

E-mail：adrxg＠bjut.edu.cn

肖 堯（198-8- ），男，博士研究生，主要研究方向?yàn)槟Ｊ阶R(shí)別與智能系統(tǒng)。

E-mail：siaoyao＠126.com

朱曉慶（1987- ），男，博士，主要研究方向?yàn)闄C(jī)器人、機(jī)器學(xué)習(xí)。

E-mail：553455117＠qq.com

黃 靜（1979- ），女，博士研究生，主要研究方向?yàn)闄C(jī)器學(xué)習(xí)。

E-mail：757488601＠qq.com

1001-506X（2015）12-2799-11

2014- 11- 27；

2015- 04- 29；網(wǎng)絡(luò)優(yōu)先出版日期：2015- 07- 06。

網(wǎng)絡(luò)優(yōu)先出版地址：http：／／www.cnki.net／kcms／detail／11.2422.TN.20150706.1606.010.html

國(guó)家自然科學(xué)基金（61075110，61375086）；國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（973計(jì)劃）（2012CB720000）；北京市自然科學(xué)基金項(xiàng)目／北京市教育委員會(huì)科技計(jì)劃重點(diǎn)項(xiàng)目（KZ201210005001）；高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金（20101103110007）；北京高等學(xué)校青年英才計(jì)劃（YETP1610）資助課題