月面著陸多尺度邊緣光流特征提取與跟蹤方法

李云天，穆榮軍，單永志，崔乃剛

(1. 哈爾濱工業(yè)大學(xué)航天工程系，哈爾濱150001；2. 中國兵器工業(yè)集團(tuán)航空彈藥研究院，哈爾濱150001)

0 引 言

月面著陸按降落時序可大致分為離軌再入、動力下降與最終著陸三個階段。其中，動力下降段飛行時間最長，空域最廣，該段的導(dǎo)航定位精度直接決定了著陸器的最終著陸精度[1]。由于月球尚未建成完備可靠的絕對定位系統(tǒng)，地面測控站與著陸器之間又存在較大的通信延遲，著陸器的導(dǎo)航精度極大程度依賴于慣性測量單元(Inertial measurement unit，IMU)、著陸相機(jī)(Landing camera，LCAM)、激光探測與測距(Light detection and ranging，LiDAR)等傳感器組成的自主導(dǎo)航系統(tǒng)[2]。其中，LCAM及基于LCAM的視覺導(dǎo)航系統(tǒng)具有造價低、重量輕、功耗小、可同時獲取著陸器位姿等優(yōu)勢[3]，和慣性導(dǎo)航系統(tǒng)并列為動力下降段的首選導(dǎo)航系統(tǒng)。

月面著陸視覺導(dǎo)航系統(tǒng)可分為視覺絕對導(dǎo)航系統(tǒng)和視覺相對導(dǎo)航系統(tǒng)兩種。視覺絕對導(dǎo)航系統(tǒng)在動力下降過程中使用正下視安裝的LCAM獲取著陸器當(dāng)前位置的月面圖像及地形特征(如隕石坑特征[4]、陰影特征[5]、混合特征[6]等)，隨后通過與現(xiàn)有月面地形數(shù)據(jù)庫匹配，獲取著陸器在月球固連坐標(biāo)系下的絕對位置。中國的嫦娥-4(CE- 4)探測器、美國NASA的“自動著陸與風(fēng)險規(guī)避技術(shù)”(Autonomous landing and hazard avoidance tech-nology，ALHAT)、德國DLR的“自動光學(xué)相對地形導(dǎo)航”(Autonomous terrain based optical navigation，ATON)等均采用了這一方式。雖然視覺絕對導(dǎo)航系統(tǒng)可獲得著陸器在月球固連坐標(biāo)系下的絕對位姿，但仍存在計算負(fù)擔(dān)重、依賴地形數(shù)據(jù)庫精度、相機(jī)斜視情況下匹配困難、平坦區(qū)域無法工作、需要耗費(fèi)大量數(shù)據(jù)存儲空間等缺點(diǎn)。

與之相比，基于序列圖像的視覺相對導(dǎo)航系統(tǒng)雖無法獲取月球固連坐標(biāo)系下的絕對位置，但由于具有不依賴月面數(shù)據(jù)庫、在平坦區(qū)域能獲取充足的特征點(diǎn)，對LCAM朝向無要求等優(yōu)勢，近年來逐漸受到了重視[6-7]。視覺相對導(dǎo)航系統(tǒng)通過提取和匹配序列圖像特征，進(jìn)而最小化重投影誤差或灰度誤差的方式得到著陸器在兩幀圖像之間的相對位姿估計。因此特征提取與跟蹤的精度和計算效率很大程度上決定了視覺相對導(dǎo)航系統(tǒng)的精度和實(shí)時性。

由于隕石坑和石塊是月面地貌的主要組成部分，視覺相對導(dǎo)航系統(tǒng)通常采用特征角點(diǎn)(Harris、FAST等)和特征氣泡點(diǎn)(SURF、KAZE等)作為跟蹤特征。其中，以FAST為代表的特征角點(diǎn)因具有較高的計算效率和良好的光照不變性[8]，被廣泛應(yīng)用于視覺相對導(dǎo)航系統(tǒng)中[9-10]。但傳統(tǒng)特征角點(diǎn)方法提取的特征點(diǎn)過于集中，存在特征聚集問題，且在尺度不變性方面仍有待進(jìn)一步提高。

在完成特征點(diǎn)提取后，視覺相對導(dǎo)航系統(tǒng)還需對所提取出的特征點(diǎn)進(jìn)行跟蹤以確定特征點(diǎn)之間的匹配關(guān)系，作為下一步位姿估計的基準(zhǔn)?；诿枋鲎拥膫鹘y(tǒng)特征跟蹤方法工作流程如圖1(a)所示。首先，基于描述子生成算法計算前后兩幀中所有特征點(diǎn)的描述子；為了保證描述子的準(zhǔn)確性和魯棒性，這一過程需要進(jìn)行大量的統(tǒng)計計算(例如BRIEF描述子需要統(tǒng)計特征點(diǎn)附近256個像素對的明暗分布情況[11])。隨后，對于參考幀中的每一個特征點(diǎn)，遍歷當(dāng)前幀中的所有特征描述子，篩選出距離(通常采用馬氏距離或漢明距離)最近的特征點(diǎn)作為匹配點(diǎn)對。上述過程的計算耗時通常占整個特征提取和跟蹤算法耗時的3/4以上，同時需要大量存儲空間記錄所有特征描述子，嚴(yán)重降低了相對視覺導(dǎo)航系統(tǒng)的實(shí)時性。

由于月球表面不存在大氣且動力下降段時間較短，可認(rèn)為LCAM捕獲的序列圖像滿足灰度不變假設(shè)，這為光流法的應(yīng)用帶來了可能，涌現(xiàn)出了若干以Kanada-Lucas-Tomasi(KLT)光流跟蹤為代表的，使用光流預(yù)測代替特征匹配進(jìn)行特征跟蹤的研究成果[12-13]。由于跳過了特征描述子的計算和匹配過程，基于光流預(yù)測的特征跟蹤算法計算效率相較基于描述子的特征跟蹤算法有了極大提升。但傳統(tǒng)光流算法基于小運(yùn)動假設(shè)，對大尺度運(yùn)動的魯棒性較差；且求解過程涉及圖像二次導(dǎo)數(shù)的計算[14]和二維圖像平面的匹配搜索[15]，對于系統(tǒng)的計算能力依舊有一定要求。

針對上述問題，本文提出一種多尺度邊緣光流(Multi-Scale Edge-Flow，M-EF)特征提取與跟蹤方法，工作流程如圖1(b)所示。首先采用帶方向描述的ORB特征替換傳統(tǒng)的FAST特征，并使用多級掩膜方法優(yōu)化特征點(diǎn)分布。隨后，基于圖像邊緣直方圖，采用灰度差平方和(Sum of square grayscale difference，SSGD)滑窗搜索算法分別計算特征點(diǎn)在圖像坐標(biāo)系兩個方向上的邊緣光流，將二維光流迭代計算問題簡化為一維直方圖匹配問題，以提高計算效率。同時，通過引入多尺度邊緣直方圖和迭代搜索算法，在有效抑制大尺度運(yùn)動下特征點(diǎn)跟蹤易丟失問題的同時，將光流計算精度由整像素級拓展至亞像素級。最終，基于嫦娥4號探測器動力下降段LCAM獲取的真實(shí)圖像序列對所提出算法的性能進(jìn)行了仿真校驗(yàn)。

圖1 傳統(tǒng)特征跟蹤算法與M-EF特征跟蹤算法流程對比Fig.1 Comparison of working flow between traditional feature tracking method and M-EF feature tracking method

1 多級掩膜ORB特征提取

ORB特征是Ethan Rublee等在FAST基礎(chǔ)上，引入方向描述和Harris角點(diǎn)評價提出的一種具有較好旋轉(zhuǎn)和尺度不變性的視覺特征[16]?，F(xiàn)有ORB算法在篩選特征點(diǎn)時完全依賴于其Harris評分，而圖像強(qiáng)特征附近的特征點(diǎn)通常會具有十分相近的Harris評分。這就帶來了特征點(diǎn)聚集問題：大量特征點(diǎn)集中于圖像中的一個或某幾個強(qiáng)特征區(qū)域，其他區(qū)域則只存在零散的特征點(diǎn)；分布缺乏均勻性，進(jìn)而影響后續(xù)特征跟蹤和位姿估計的精度。

針對這一問題，提出一種構(gòu)建圖像金字塔，并采用多級掩膜(Mask)，分級分批提取特征點(diǎn)的多尺度ORB(Multi-scale ORB，M-ORB)特征提取算法，流程概述如下：

1)以當(dāng)前幀為底層，逐層下采樣(Down Sample)構(gòu)建圖像金字塔；

2)在金字塔的頂層進(jìn)行ORB特征提取，并根據(jù)Harris評分篩選出評分靠前的若干特征點(diǎn)；

3)將所提取特征點(diǎn)的坐標(biāo)上采樣(Up Sample)至下一層；

4)在該層圖像上以上采樣坐標(biāo)為圓心繪制圓形掩膜，掩膜內(nèi)的區(qū)域在本層特征提取過程中將被忽略；

5)重復(fù)步驟2)至3)直至第0層；

6)從第0層未被掩膜覆蓋的圖像區(qū)域中一次性提取出剩余數(shù)量的特征點(diǎn)；

7)將每層提取出的特征點(diǎn)坐標(biāo)統(tǒng)一上采樣至第0層，即為所需的特征點(diǎn)集。

1.1 圖像金字塔構(gòu)建

記需要提取特征點(diǎn)的原始圖像為圖像金字塔的第0層I0(x,y)，則圖像金字塔可通過下述下采樣過程構(gòu)建：

1)擴(kuò)充第L-1層圖像邊緣，將圖像四周高斯卷積核半窗口大小內(nèi)的所有像素灰度填充為0；

2)對圖像中的所有像素點(diǎn)應(yīng)用式(1)所示的高斯卷積核G重新計算灰度；

(1)

3)按照式(2)計算第L層圖像對應(yīng)位置的像素灰度值，式中x和y分別為圖像點(diǎn)在第L層圖像坐標(biāo)系上的橫縱坐標(biāo)；

(2)

4)重復(fù)步驟1-3直到頂層，得到圖像金字塔集合{IL}L=0,…,k。

顯然，對于金字塔的第k層圖像來說，特征檢測和光流計算的范圍將分別增大4k和2k倍。雖然理論上金字塔層數(shù)越多，特征檢測和跟蹤的魯棒性越高，但相應(yīng)的計算負(fù)擔(dān)也會加重。因此在實(shí)際應(yīng)用中，k的取值范圍通常為2～4。

1.2 多級掩膜特征提取

所謂多級掩膜特征提取，即在完成圖像金字塔的構(gòu)建后，首先在尺度最小的頂層圖像上提取若干特征點(diǎn)；這部分特征點(diǎn)所處的區(qū)域在下一個尺度圖像的特征提取過程中將被設(shè)置為掩膜從而被忽略；循環(huán)上述過程直至遍歷所有尺度的圖像，獲得最終的特征點(diǎn)集合。

圖2 采用多級掩膜的多尺度ORB提取算法流程圖Fig.2 Flow chat of multi-scale ORB extraction method based on multi-level mask

多級掩膜特征提取中最重要的兩個參數(shù)為每層的特征點(diǎn)數(shù)目nL和掩膜半徑rm，由于每次下采樣后圖像尺寸變?yōu)樵瓐D像的1/4，這里不失一般性地假設(shè)每層的特征數(shù)目和掩膜半徑滿足：

nL-1=4nL, (rm)L-1=2(rm)L

(3)

設(shè)需要提取的特征點(diǎn)總數(shù)為N，金字塔頂層為k層，則第L層需要提取的特征點(diǎn)數(shù)目nL為：

(4)

同理，第L層的掩膜半徑(rm)L為：

(5)

式中：w和h分別為圖像的寬度和高度。

2 多尺度邊緣光流特征跟蹤

隨著相機(jī)的移動，像素點(diǎn)在兩幀圖像之間的運(yùn)動就像水流一樣，故稱之為光流(Optical Flow)，單位為pixel/s。由于光流表征了像素點(diǎn)的運(yùn)動方向和大小，基于光流預(yù)測的特征跟蹤方法可有效縮小特征跟蹤范圍，在提高特征跟蹤準(zhǔn)確度的同時有效減輕計算負(fù)擔(dān)。

基于同一特征區(qū)域在前后兩幀圖像中應(yīng)具有相似的空間邊緣分布這一假設(shè)，文獻(xiàn)[17]提出了一種基于邊緣直方圖的光流計算方法——邊緣光流(Edge-Flow，EF)。EF算法在小尺度運(yùn)動下具有較好的跟蹤效果和較低的計算耗時，但在大尺度運(yùn)動下的計算效果并不理想，且只能獲得整像素級光流。本文在此基礎(chǔ)上，進(jìn)一步提出了多尺度邊緣光流算法(Multi-scale Edge-Flow，M-EF)，在提高大尺度運(yùn)動下特征跟蹤魯棒性的同時，將光流計算精度由整像素級拓展至亞像素級。

2.1 邊緣光流

EF光流計算的核心是圖像邊緣直方圖，而圖像的邊緣可由圖像一階導(dǎo)數(shù)描述。因此在獲取序列圖像后，分別對前后兩幀圖像應(yīng)用式(6)所示的Sharr算子計算圖像x和y方向上一階導(dǎo)數(shù)。隨后，分別將圖像中所有像素點(diǎn)的一階導(dǎo)數(shù)值按列、按行求和作為邊緣數(shù)值，即得到圖像在x和y方向上的邊緣直方圖。

(6)

顯然，相同圖像點(diǎn)在兩幀邊緣直方圖上應(yīng)具有相同的邊緣點(diǎn)數(shù)目。但由于直方圖中具有相同邊緣數(shù)值的坐標(biāo)位置并不唯一，具有相同邊緣數(shù)值僅是兩個坐標(biāo)位置是同一特征點(diǎn)的必要而非充分條件。因此，需要采用圖3所示的灰度差平方和(SSGD)滑窗搜索，綜合比較待不同坐標(biāo)位置周圍的灰度分布，減少誤匹配，提高光流計算精度。

記q為圖像點(diǎn)在參考幀邊緣直方圖中的坐標(biāo)，則其在當(dāng)前幀邊緣直方圖中的坐標(biāo)t應(yīng)滿足：

(7)

式中:p為滑動搜索區(qū)間的半長度，(q,t,s)為灰度差平方和(SSGD)函數(shù)：

(8)

式中:hr(·)為參考幀的邊緣直方圖，hc(·)為當(dāng)前幀的邊緣直方圖，s為搜索窗口的半長度。

圖3 邊緣光流SSGD滑窗搜索算法Fig.3 SSGD based sliding window search for EF

2.2 多尺度邊緣光流

原始EF算法存在兩個缺陷:1)為了跟蹤大尺度運(yùn)動，需要設(shè)置較大的搜索窗口和搜索區(qū)間，一定程度上增加了計算負(fù)擔(dān)；2)滑窗搜索以像素坐標(biāo)為中心，只能獲取整像素級光流。本文通過構(gòu)建多尺度邊緣直方圖金字塔，并逐層搜索計算邊緣光流，在一定程度上解決了這一問題，算法流程如圖4所示。

由于邊緣直方圖僅有一維，這里采用一維高斯加權(quán)矢量替代式(1)的高斯卷積核進(jìn)行邊緣直方圖金字塔的構(gòu)建，第L層邊緣直方圖x坐標(biāo)處的邊緣數(shù)值可用L-1層邊緣直方圖表示為：

6hL-1(2x)+4hL-1(2x+1)+

hL-1(2x+2)]

(9)

圖4 多尺度邊緣光流算法流程圖Fig.4 Flow chart of Multi-Scale Edge-Flow method

邊緣直方圖金字塔構(gòu)建完成后，首先在最上層(k層)進(jìn)行滑窗搜索，將得到的搜索結(jié)果變換至k-1層作為該層的搜索起點(diǎn)；隨后，逐層重復(fù)上述過程直至獲取第0層坐標(biāo)t0；此時獲取的光流v0=t0-q為與真實(shí)亞像素級光流最為接近的整像素級光流，直接對完整邊緣直方圖進(jìn)行上采樣會帶來許多不必要的計算。因此這里首先對以t0為中心，半徑為[-1,+1]區(qū)間內(nèi)所有坐標(biāo)位置處的邊緣數(shù)值進(jìn)行線性插值，隨后在該區(qū)間內(nèi)尋找滿足式(10)的亞像素修正量Δt，最終得到亞像素級光流vs=t0+Δt-q。

(10)

3 仿真校驗(yàn)

為驗(yàn)證所提出的多尺度邊緣光流特征提取與跟蹤方法的有效性和可行性，基于嫦娥4號探測器LCAM獲取的真實(shí)圖像序列進(jìn)行了仿真校驗(yàn)。仿真程序基于C++開發(fā)并部署在Linux環(huán)境的Robot Operation System中；仿真平臺硬件配置為Intel(R) Core(TM) i5-8300@2.3 GHz CPU+8 GB RAM+GTX 750M GPU；為保證算法耗時比較的準(zhǔn)確性，所有GPU加速功能在仿真中均未開啟，全部運(yùn)算均由CPU獨(dú)立完成。

CE-4動力下降時序如圖5所示，主要包括著陸準(zhǔn)備、主減速、快速調(diào)整、垂直接近、懸停避障和緩速下降六個階段。其中，著陸器在著陸準(zhǔn)備段、主減速段和垂直接近段的飛行時間最長；在快速調(diào)整段和垂直下降段分別存在由持續(xù)姿態(tài)機(jī)動和快速水平機(jī)動引起的大尺度運(yùn)動(大于10 pixel/s)[18](圖5中雙線框部分)。

圖5 CE-4動力下降段降落時序Fig.5 Landing sequence of CE-4 power descent

圖6 M-ORB與ORB特征點(diǎn)坐標(biāo)分布直方圖對比Fig.6 Histogram of feature distribution between M-ORB and ORB

CE-4動力下降過程中所使用的LCAM的視場角為45.3°，焦距為8.5 mm，有效像素數(shù)為1024×1024，圖像采集幀率為10 FPS。圖7給出了所提出的M-ORB算法和傳統(tǒng)ORB算法在不同特征點(diǎn)數(shù)目(依次為200、500和1000)下的兩軸坐標(biāo)分布直方圖?？梢钥闯?，三種情況下M-ORB算法所提取特征點(diǎn)分布的均值和方差均小于傳統(tǒng)ORB算法，特征點(diǎn)的空間分布更加均勻。

隨后，進(jìn)一步對M-ORB+BRIEF、M-ORB+KLT、M-ORB+EF及M-ORB+M-EF四種算法在真實(shí)場景下的性能進(jìn)行了仿真校驗(yàn)和比較分析(由于特征提取部分均為M-ORB算法，為書寫方便，下文以追蹤部分算法名稱代指)，各項(xiàng)算法仿真參數(shù)如表1所示。跟蹤精度評估方面，在獲得特征跟蹤結(jié)果后，采用如式(11)所示絕對差值和(Sum of absolute difference，SAD)算法對每種跟蹤算法的正確跟蹤特征點(diǎn)數(shù)目進(jìn)行重校驗(yàn)以保證評估指標(biāo)的一致性，式中(xr,yr)和(xc,yc)為分別為特征點(diǎn)在前后幀圖像中的像素坐標(biāo)，k為窗口尺寸，通常取奇數(shù)。計算耗時評估方面，則采用相對計算耗時(特征跟蹤耗時占總時間的比例)作為評估指標(biāo)以避免隨機(jī)性影響。

表1 各算法仿真參數(shù)設(shè)置Table 1 Simulation parameters for all tracking methods

圖7 CE-4動力下降各段不同跟蹤算法平均跟蹤特征點(diǎn)數(shù)目Fig.7 Average tracked feature numbers of different tracking methods during entire CE-4 power descent

圖8 CE-4動力下降各段不同跟蹤算法相對計算耗時Fig.8 Relative elapsed time of different tracking methods during entire CE-4 power descent

表2 CE-4動力下降各段平均跟蹤特征點(diǎn)數(shù)目Table 2 Summary of average tracked feature numbers of different tracking methods during entire CE-4 power descent

表3 CE-4動力下降各段不同跟蹤算法相對計算耗時Table 3 Summary of relative elapsed time of different tracking methods during entire CE-4 power descent

(11)

四種算法在動力下降各段的平均跟蹤特征點(diǎn)數(shù)目和相對計算耗時分別如圖7～8和表2～3所示。仿真結(jié)果表明，基于描述子的BRIEF算法在所有階段的平均跟蹤特征點(diǎn)數(shù)目變化幅度最小，對于大尺度運(yùn)動的魯棒性更好，但隨之帶來了較大的計算負(fù)擔(dān)，相對計算耗時在13.86%左右。

三種基于光流的跟蹤算法在較為平穩(wěn)的降落階段具有更高的平均跟蹤特征點(diǎn)數(shù)目，其中EF及改進(jìn)的M-EF算法略高于KLT算法。但當(dāng)存在較大尺度運(yùn)動的快速調(diào)整和垂直接近段，三種光流跟蹤算法的平均跟蹤特征點(diǎn)數(shù)目均出現(xiàn)不同程度下降，其中以EF算法最為嚴(yán)重。KLT及所提出的M-EF算法由于引入了多尺度策略，在一定程度上提高了大尺度運(yùn)動下的特征跟蹤穩(wěn)定性，相比BRIEF算法下降在15%左右。

在相對計算耗時方面，三種光流跟蹤算法均不超過BRIEF算法的50%，其中未采用多尺度策略的EF算法相對耗時僅為2.56%左右，采用多尺度策略的M-EF和KLT算法計算耗時則依次略有增加，分別在5.12%和6.83%左右。

4 結(jié) 論

基于描述子的傳統(tǒng)視覺特征提取和跟蹤方法存在耗時長、誤匹配率高等問題，不能滿足月面著陸器動力下降段自主導(dǎo)航的精度和魯棒性要求。本文結(jié)合圖像邊緣直方圖和圖像金字塔策略，提出了一種多尺度邊緣光流特征提取與跟蹤方法、該方法在有效實(shí)現(xiàn)光流計算降維，提高計算效率的同時，改善了大尺度運(yùn)動下特征跟蹤的魯棒性，將光流計算精度拓展至亞像素級。仿真結(jié)果表明，本方法計算耗時不超過傳統(tǒng)描述子方法的50%；在非大尺度運(yùn)動下具有最高的穩(wěn)定跟蹤特征點(diǎn)數(shù)目；大尺度運(yùn)動下相比描述子方法下降不超過15%，優(yōu)于原始EF算法，與同樣采用多尺度策略的KLT算法基本持平。本方法在跟蹤效率和穩(wěn)定性方面取得了較好的平衡，是對傳統(tǒng)方法的有效改進(jìn)，具有一定的工程應(yīng)用價值。