面向山地區域光照變化下的魯棒景象匹配方法

王華夏，程詠梅，劉楠

1.西北工業大學 自動化學院，西安 710072 2.太原科技大學 電子信息學院，太原 030024

面向山地區域光照變化下的魯棒景象匹配方法

王華夏1,2,*，程詠梅1，劉楠1

1.西北工業大學 自動化學院，西安 710072 2.太原科技大學 電子信息學院，太原 030024

提出了一種面向山地區域光照變化下的魯棒景象匹配導航方法。該方法針對山區景象匹配導航中光照環境不同造成的基準圖與實時圖像不一致問題，采取在高程數據上使用光照模型生成光照明暗圖，利用數字高程圖(DEM)得到水流匯集數據并生成山谷顯著圖，將光照明暗圖與山谷顯著圖融合作為基準圖；針對基準圖與實時圖中細節邊緣的差異帶來的誤匹配問題，提出了基于形態學約束的Hausdorff距離邊緣匹配算法。采用LANDSAT圖像與ASTERDEM高程數據進行實驗分析，結果表明提出的方法匹配正確率高且魯棒性好。

景象匹配導航；光照模型；顯著邊緣；Hausdorff距離；圖像融合

圖像匹配導航是利用機載基準圖像與實時獲取的圖像進行匹配，確定出飛行器位置。該技術由于具有自主性和定位精度高的優點成為一種重要的導航方式，與慣性導航系統(INS)組合可實現長航時、高精度的自主導航。圖像匹配導航技術可以分為地形高程匹配導航和景象匹配導航兩大類，其中地形高程匹配導航是采用激光雷達、合成孔徑雷達干涉(INSAR)或序列圖像重建地形的方法[1-5]獲取實時地形，與機載基準數字高程圖(Digital Elevation Model,DEM)進行匹配，實現飛行器定位，激光雷達與INSAR獲取的地形數據精度高但設備成本昂貴，而可見光序列圖像對地形進行恢復要求有充足的圖像重合區域和稠密穩定的特征點，在山區中不易獲取[6-7]；景象匹配導航技術是利用合成孔徑雷達(SAR)或光學圖像傳感器實時獲取地面景物圖像與機載基準圖像進行匹配，實現飛行器定位，目前景象匹配主要針對具有特定結構以及人造目標的城鎮區域，在缺乏人造建筑、顏色特征單一的山區，難以提取穩定而獨特的結構與色彩特征[8-9]，并且由于成像時間不同、光照條件不同使得實時獲取的圖像與基準圖像在灰度以及邊緣信息上存在較大差異，若采用固定圖像作為基準圖，易發生誤匹配，故對于山區景象匹配的基準圖像制備是一個研究的熱點和難點。文獻[10]針對山地區域的景象匹配導航提出了一種基于GRRR(Geo-Registered Radar Returns)數據構造可見光仿真圖像的景象匹配方法，但構造基準圖時沒有充分考慮山區中反射系數差異較大的山谷區域。而利用DEM數據可獲取山谷水文信息[11]，結合DEM數據與光照模型以及水文信息構造一種抗光照變化、計算快速、可在線生成的基準圖是一種新穎的思路。

山區景象的邊緣特征豐富，通過模型構造的基準圖像與實時圖像在圖像的細節邊緣上存在差異，如何在這種差異存在的情況下實現正確的局部匹配是景象匹配的另一個問題。由于霍斯多夫距離(Hausdorff Distance，HD)具有對邊緣局部形變不敏感、不要求所有匹配點嚴格對應、計算快捷的優點，因此在圖像匹配中得到了廣泛的應用，但HD極易受到出格點的影響，當邊緣圖像受到噪聲污染或存在遮擋等情況時，易造成誤匹配。文獻[12]提出了PHD (Partial HD)，在距離集中選擇特定排位的距離值作為單向HD以減小出格點的影響；文獻[13]提出了MHD (Modified HD)，對距離集取均值以減小出格點的影響，但PHD與MHD對出格點的排除能力有限；文獻[14]提出了LTSHD (Least Trimmed Square HD)與M-HD (M-estimation HD)，這兩種算法采用魯棒回歸與野點檢測的思想進一步提升了HD排除出格點的能力；文獻[15]在M-HD的基礎上提出了AM-HD (Accurate M-HD)距離，引入邊緣的方向信息對M-HD進行權重調整，雖然以上方法在處理遮擋與退化圖像的匹配問題時已具有較好的性能，但當基準圖與實時圖的細節邊緣差異較大、實時圖存在噪聲的情況下，誤匹配率較高。

針對以上問題，本文提出一種新型的景象匹配導航方法。首先在DEM數據上采用光照模型生成明暗圖像；然后由DEM數據提取山谷信息，以表達反射特性與周圍環境差異較大的山谷河流區域；隨后融合光照明暗圖與山谷信息圖像生成融合圖像作為基準圖，對實測圖像與基準圖像采用MLSEC(Multi-local Level-Set Extrinsic Curvature)方法提取顯著的邊緣信息；在邊緣匹配過程中，采用了形態學方法對初始邊緣點集進行約束，利用約束前后的公共鄰近邊緣點數量對HD進行權值校正，進一步減少了出格點、提高了正確匹配位置的HD峰值特性。最后采用LANDSAT圖像與ASTERDEM數據驗證了本文方法的魯棒性與有效性。

1 光照變換下的魯棒基準圖制備方法

由于光照條件不同，三維物體的成像不同，如圖1所示，圖1(a)、圖1(d)為不同時間段拍攝的秦嶺山區可見光圖像，兩幅圖像成像的光照條件差異較大(圖1(a)拍攝時間為2014/12/30，當地太陽高度角為30°，太陽方位角為156°；圖1(d)拍攝時間為2016/7/28，當地太陽高度角為64°，太陽方位角為118°)；圖1(b)、圖1(e)為圖1(a)、圖1(d)對應的CANNY[17]邊緣特征；圖1(c)、圖1(f)為圖1(a)、圖1(d)對應的MLSEC邊緣特征。本文采用互相關方法與Hausdorff測度，分別對圖1中的3種圖像進行了計算分析：圖1(a)與圖1(d)全局互相關系數為0.47，圖1(b)與圖1(e)的全局HD距離為27.4，圖1(c)與圖1(f)的全局HD距離為6.8；隨機選取圖1(d)、圖1(e)、圖1(f)中的部分區域分別與圖1(a)、圖1(b)、圖1(c)進行局部匹配定位并進行蒙特卡羅仿真， 由表1可見，3種圖像的局部匹配正確率都非常低。

MatchmethodNormalizedcorrelationHDHDMapstomatchFig．1(a)，Fig．1(d)Fig．1(b)，Fig．1(e)Fig．1(c)，Fig．1(f)Accuracy/%43824

由此可知，光照條件差異較大造成圖像灰度差異大、全局相關性低、邊緣特征不一致、圖像不適配。若采用固定時刻拍攝的圖像作為基準圖，無法保障正確匹配以及后續的飛行器準確定位，因此需要尋找一種抗光照變化的基準圖制備方法。

1.1 光照明暗圖

由于三維物體所成的圖像取決于：① 物體的形狀；② 物體表面對光的反射特性；③ 光源的分布。地物所在的地形環境及相對太陽的位置不同，使得同一時刻不同地形條件下的地表所接收到的光輻射強度不同，而在不同時刻同一地形的輻射強度也會有差異[18]。

Lambert光照模型如圖2所示，設輻照強度為E的光源S照射在Lambert表面的象元O上，灰色圓形區域是光線入射點O所處地形區域的切平面，則該象元的輻射強度L為

L=Ecosγ

(1)

式中:γ為光線OS與物體表面法向量fON的夾角。當象元O處的切平面水平時，fON與天向軸fOZ重合，γ與太陽天頂角θz相同；當象元O處的切平面不水平時，以O為原點的單位球面與fOS、fON、fOZ相交于點IO、NO、ZO，構成球面三角形ΔIONOZO，根據球面三角形關于邊的余弦定理[19]，可知

cosγ=cosθpcosθz+sinθpsinθzcos(φa-φo)

(2)

式中：θp為地形坡度角，是地形切平面與水平面的夾角;φo為地表法線方位角，為ON在水平面投影向量與正北方向的夾角，θp、φo與DEM數據相關，通過DEM象元及其鄰域象元的關系求取；φa為太陽方位角;θz為太陽天頂角;θH為太陽高度角，其中θz+θH=π/2，θH、φa與太陽位置和當地位置有關，通過時間和當地經緯可推算得出。

在天氣晴好的情況下，基于平行光源假設的地表光照模型反映了物體表面所接收到的太陽輻射量，可通過為DEM柵格數據中的每個高程象元指定太陽高度、太陽方位角，來計算地表的模擬亮度值，設光照強度E=1，觀測方向為正下視，根據成像位置DEM數據以及成像時間與當地的經緯可生成光照明暗灰度圖像：

Map_hillshade(i,j)=255cosγi,j

(3)

式中：γi,j為DEM(i,j)處表面法向量與太陽光線的夾角。

由于基于平行光源的Lambert模型未考慮陽光被大面積遮擋或入射光線為散射光時的情況，故陰天以及雨雪天氣不在本文研究范圍內。由此模型在不同時間生成的光照明暗圖如圖3所示，展現了由光照模型與地形數據生成的明暗信息。由于光照明暗圖的生成是基于物體表面反射系數一致的情況，而山區包括山谷、道路與植被覆蓋的山體，反射系數不完全一致，因此明暗圖不能完整的描述成像信息。

1.2 山谷顯著圖

山區可見光圖像的顯著邊緣除了由光照產生的明暗邊緣外，其余的顯著邊緣主要集中于山谷，其原因是山谷區域表面無植被覆蓋，其光學反射率與植被覆蓋的山地差異較大。山谷大多由水流沖刷產生，而DEM數據本身包含了地形的高低變化的規律[20]，因此可通過分析DEM數據的高度變化而獲取水文信息并得到山谷顯著圖。

圖2 Lambert光照模型Fig.2 Lambert illumination model

圖3 同一地形區域不同時間產生的光照明暗圖Fig.3 Hillshade map at different times in same terrain

山谷信息及其顯著圖像的獲取步驟如下：

1) 水流方向。首先依據DEM像素單元八鄰域的數據下降梯度來指定水流方向，并進行水流方向的編碼：

MD(i,j)=

(4)

式中:DEM(i,j)為該單元的高程值;S為元素DEM(i,j)八鄰域的像素位置集合，第k個鄰域的位置為Sk=[pk,qk](k=0,1,…,7);水流方向MD(i,j)記錄最大梯度下降方向的鄰域編碼位置。若中心像素位置對應的高程數據低于周圍八鄰域的高程數據，該中心的高程值設定為八鄰域中最低高程；若同時存在多個最大下降梯度方向，采用查表方式來進行水流方向的判斷；若所有鄰域的梯度相同，采用遞歸的方式進行水流方向的計算[21]。

2) 水流匯集數據。由水流方向數據可知像素單元(i,j)的八鄰域的水流方向，標記其中水流方向指向中心(i,j)的單元格為1，其他為0；按行或列依據標記數據矩陣flag對水流匯集數據FA進行遍歷累加，在遍歷過程中更新并得到水流匯集數據，水流匯集數據如式(5)和式(6)所示，其中FA(i,j)的初始值為0。

flag(i,j,Sk)=

(5)

(6)

3) 山谷顯著圖。水流匯集數據FA(i,j)的單元累加值越大，說明匯集到該像素中的水流越多，越可能屬于山谷地區，由于水流匯集數值并不均勻分布，若直接將水流匯集數據轉化為灰度圖像，將得不到理想可視化山谷信息(如圖4(a)所示)。

采用閾值th對水流匯集數據進行規范后作為山谷信息的顯著圖像，其中閾值th為DEM總像素數量的1.5%，即當水流匯集數據大于th時不再增大其在圖像中的亮度影響，減少灰度級別，增大主要山谷區域的對比度。山谷信息圖像生成算法如式(7)所示。圖4(b)為由DEM數據生成的山谷信息顯著圖像，該圖像以灰度圖的形式展現了DEM數據中的山谷。

Map_valley(i,j)=

(7)

圖4 山谷顯著圖像Fig.4 Images of significant valleys

1.3 光照明暗圖與山谷顯著圖的融合圖像

光照明暗圖像缺失了山谷信息所呈現的顯著邊緣，而山谷顯著圖無法表達由光照引起的明暗交界的顯著邊緣。故采用加權融合方法，將光照明暗圖與山谷顯著圖融合，生成一種包含兩種信息的融合圖像，融合方法如式(8)所示。

Map_fusion=

W·Map_vally+(1-W)·Map_hillshade

(8)

式中:W(0≤W≤1)為融合權重，當W趨近0時，融合圖像退化為光照明暗圖，當W趨近1時，融合圖像退化為山谷顯著圖。

從真實圖像與融合圖像相似性的角度出發，采用最小均方根誤差(Root Mean Square Error,RMSE)方法考察融合權值W對圖像質量的影響，圖5為不同W權重下，歸一化的融合圖像與真實圖像的RMSE，圖6展示了W權重對Hausdorff算法局部匹配正確率的影響，統計數據表明在W取值為0.3～0.4時，融合圖像與真實圖像的RMSE最小，同時對應的4種Hausdorff匹配算法的局部匹配正確率也最高。通過大量真實圖像與融合圖像的RMSE后驗計算可獲取W的經驗值為

(9)

圖5 不同融合權重下融合圖像與真實圖像的RMSEFig.5 RMSE of fusion and real images with different fusion weights

融合圖像與其邊緣圖像如圖7所示，其中圖7(a)為可見光圖像，圖7(b)為光照明暗圖，圖7(c)為山谷顯著圖像，圖7(d)為融合圖像，圖7(e)、圖7(f)、圖7(g)、圖7(h)分別是圖7(a)、圖7(b)、圖7(c)、圖7(d)提取的MLSEC邊緣圖像，其中對圖7(c)采用MLSEC邊緣提取是為了保持邊緣機理的一致性。由圖可見，圖7(h)不僅具有光照產生的明暗邊緣，還包含較為完善的山谷顯著邊緣信息，與可見光邊緣圖7(e)更具一致性。在317×317大小的圖7(e)中隨機截取60×60大小區域作為實時邊緣圖像，采用4種HD方法分別與圖7(f)、圖7(g)、圖7(h)進行局部匹配，設定定位誤差小于5個像素單位為正確匹配，對匹配正確率進行蒙特卡羅統計實驗，實驗結果如表2所示，以圖7(h)為基準圖的局部匹配正確率最高，而以圖7(f)為基準圖的匹配正確率高于以圖7(g)為基準圖的匹配正確率，由此說明：光照明暗圖產生的邊緣信息在邊緣匹配中做了主要貢獻；由于山谷邊緣信息稀疏，使得山谷顯著邊緣作為基準圖的匹配正確率很低；而采用融合圖像的邊緣圖像作為基準圖得到了提升匹配正確率的效果；同時也側面驗證了光照明暗圖在融合圖像中起的作用較大，應給予更高的權值。

圖6 不同融合權重對HD局部匹配正確率的影響Fig.6 Influence of different fusion weights on accuracy of HD local matching

圖7 可見光圖像、光照明暗圖、山谷顯著圖、融合圖像及其MLSEC邊緣圖像Fig.7 Vision maps,hillshade maps,valley maps,fusion maps and MLSEC edge images

表2 圖7邊緣圖像局部匹配結果Table 2 Local matching results of edge images in Fig.7

雖然融合圖像與真實可見光圖像在顯著邊緣上有較大的一致性，但由于不具備詳盡的光照模型與地面反射系數，融合圖像無法還原所有真實圖像的細節，且由于實時圖像噪聲的影響，使得融合圖像與真實圖像在邊緣檢測階段仍存在細節邊緣的差異，這種差異減少了相近邊緣點的數量和比例，并增大出格點的數量與比例，這兩種因素都會降低Hausdorff匹配的正確率。

2 EHD邊緣匹配算法

設P、Q為兩個邊緣圖像的邊緣點集為

傳統Hausdorff雙向距離的定義如式(10)所示。

H(P,Q)=max{h(P,Q),h(Q,P)}

(10)

式中:h(P,Q)為點集P到點集Q的單向HD用來衡量兩個點集的相似程度，雖然目前已有的改進Hausdorff算法具有一定的克服出格點與處理差異點集的能力，但未從全局的角度考慮誤匹配的問題。

形態學中的膨脹操作是將結構元素B平移后與原始點集E求交集，所獲取的非空點集Y為膨脹操作構成的點集。膨脹操作的定義為

Y=B⊕E={y:B(y)∩E≠?}

(11)

首先對基本邊緣點集P、Q進行膨脹操作得到PB、QB，其中PB=B⊕P，QB=B⊕Q，B為圓形結構元素，圓形半徑為r。

其次分別求取新的邊緣點集P′、Q′，其中P′=P∩QB，Q′=Q∩PB，求取新點集的目的是為了生成P、Q兩幅邊緣圖像的鄰近公共邊緣，由此得到4個邊緣點集各自的邊緣點數量：NP、NQ為點集P、Q的邊緣點數，NP′、NQ′為新點集P′、Q′的邊緣點數。NP′、NQ′將P、Q邊緣圖像的細節差異排除在有效點集數量之外。

記DQ(px,py)為點集P中的一點p與點集Q所有點的最小歐式距離，即

(12)

定義EHD(Expansion HD)的單向HD距離為

(13)

式中:ρ為對稱凸函數，使膨脹半徑外的出格點衰減為距離值r,其表達式為

(14)

最終的EHD距離為

HEHD(P,Q)=

max{h(P,Q),h(Q,P)}·min{RP,RQ}

(15)

式中:RP與RQ為校正系數,表達式分別為

(16)

采用HD距離的圖像匹配導航的問題描述為

(17)

式中：U為實時邊緣圖像；I為基準邊緣圖像；Vt為I中的遍歷子圖；t為子圖V在I中的偏移位置；PU為U的邊緣點集，QVi為Vi的邊緣點集；Rθ,s為U與Vi之間的旋轉尺度變化矩陣，表達式為

本文提出的融合圖像基準圖可在線制備，假設飛行器航向與高度數據較為精準，Rθ,s可近似為單位矩陣，式(17)簡化為求取平移因子t的問題:

(18)

通過遍歷計算每一個子圖Vt與實時圖U的EHD距離，尋找全局最小EHD的子圖，確定出位置偏移t。由于細節邊緣的差異對EHD距離產生干擾，可能出現多個子圖Vt在式(15)的第1項有相同或相近值而導致錯誤匹配的發生，此時則依據最大的公共邊緣點比率來進行權值調整，也即選取其中具有最小R值的子圖，目的是為了更大程度減小公共邊緣點比率高的HEHD值，因此式(15)的第2項是在RP與RQ中選取較小的一個。公共點集比例越大，給予的權值越小、公共點集比例越小，給予的權值越大，將具有相似HD匹配位置依據最大鄰近點集的數量進行權值調整與區分，排除了細節邊緣差異對HD的干擾，增強了正確匹配位置的區分度。

目前已有的HD方法只能給出經驗參數來設置參與計算的邊緣點比例，無法進行自適應調整，本文采用膨脹約束的方法給出動態調整HD的方法，以提高HD對細節邊緣差異的魯棒性。

3 魯棒景象匹配方法

由以上分析可知，利用DEM數據與太陽光照信息以及山谷顯著邊緣生成的融合圖像，可解決實時圖與基準圖在光照變化下顯著邊緣不一致的問題；利用形態學邊緣點集約束的方法調整HD距離的系數，可解決融合圖像與實時圖像細節邊緣差異帶來的誤匹配問題。

基于此，提出面向山地區域光照變化下的魯棒景象匹配方法，該方法表述如下：

1) 飛行器在山區飛行過程中，依據當前指示的位置以及定位誤差上限，劃定一個較大的區域作為基準圖范圍。根據飛行參數中的時間、經緯度，可推算所處區域的太陽方位角與太陽高度角，采用光照模型對基準圖范圍內的DEM數據生成光照明暗圖。

2) 對基準圖范圍的DEM數據計算其水流方向與水流匯集數據并生成山谷顯著圖像。

3) 融合光照明暗圖與山谷顯著圖生成融合圖像作為基準圖。

4) 對基準圖與實時獲取的實測圖像進行邊緣提取。

5) 采用EHD邊緣匹配方法對邊緣圖像進行匹配，實現飛行器定位并修正定位誤差。

基于光照模型與水文信息的抗光照變化的基準圖生成方法與圖像匹配流程圖如圖8所示。

雖然慣導系統指示的經緯位置是存在誤差的，但在漂移誤差范圍內對太陽高度角與太陽方位角計算所產生的影響非常小(北緯度33°～35°范圍內，5 km的緯度誤差在太陽方位角上產生的偏差僅有0.15°，在太陽高度角上產生的偏差最大僅有0.17°)，反映到光照模型上帶來的光照明暗誤差微乎其微。

選擇邊緣提取算法時考慮兩個方面：① 由于實時圖與基準圖來源不同，實時可見光圖像具備更多的細節紋理信息，而基準可見光圖在細節紋理信息上缺失，提取更多的細節邊緣對匹配效果不會有更大的提升，并且Hausdorff邊緣匹配的效率與點集數量相關;② 為了精確定位，盡量不采用大尺度的濾波器對原圖像進行高斯平滑。CANNY算法在小尺度上可以得到精準的邊緣信息，但帶來干擾的細節信息較多；若采用大尺度的CANNY算法提取圖像邊緣，由于平滑作用，原圖像的邊緣位置產生較大的位置偏移，給匹配定位精度帶來一定的損失。因此本文采用MLSEC算法提取圖像邊緣，該方法利用灰度圖梯度向量場的方法計算圖像中梯度方向匯集/相反的邊緣，對圖像中的顯著邊緣有良好的效果。

圖8 基準圖生成方法與圖像匹配流程圖Fig.8 Flowchart of reference image generation and image matching

4 景象匹配實驗

本文采用Landsat8衛星圖像作為實測圖像數據，Landsat數據也為計算光照明暗圖提供了拍攝時間、太陽高度角、太陽方位角等信息，DEM數據采用ASTER2DEM作為基準數據，并將DEM圖像與可見光圖像都規范至30 m的圖像分辨率，選取秦嶺山區典型山地地形作為仿真實驗區域。實驗平臺采用i5-3230m 2.6 GHz CPU，內存8 G，操作系統為Windows 7的筆記本電腦。在生成光照明暗圖時，由于DEM數據本身的不平滑造成在一些點位上的光照信息出現突變，會產生椒鹽噪聲點，為保證后續生成仿真圖像的平滑性，采用中值濾波濾去光照明暗圖的噪聲點。生成融合圖像時的融合權值W均設定為經驗值0.4。MLCES邊緣提取的濾波參數設置為4，在EHD方法中采用的膨脹結構為圓形，半徑r設置為4。

為了驗證在光照差異情況下采用融合圖像作為基準圖和采用固定時間圖像作為基準圖的匹配性能差異，選取秦嶺同一地區不同時間的兩幅可見光圖像進行實驗， 如圖9所示， 圖9(a)拍攝時間為2015/12/17，當地太陽高度角為30°，太陽方位角158°，圖9(b)拍攝時間為2016/7/28，當地太陽高度角為64°，太陽方位角118°，圖9(c)為依據圖9(a)的太陽高度角與太陽方位角采用本文方法生成的融合圖像，圖9(d)、圖9(e)、圖9(f)分別是由圖9(a)、圖9(b)、圖9(c)提取的MLSEC邊緣圖像。在317×317大小的圖9(d)中隨機截取100幅60×60大小邊緣圖像作為實時邊緣圖像，將圖9(e)、圖9(f)分別作為基準圖，采用5種HD方法進行局部匹配，局部匹配結果如表3所示。由圖可見，圖9(d)、圖9(e)在明暗邊緣處的一致性較差，圖9(d)、圖9(f)在明暗邊緣與河谷邊緣處的一致性均較好；以圖9(f)作為基準圖的5種HD方法的匹配正確率均高于以圖9(e)作為基準圖的匹配正確率。驗證了本文提出的基于光照模型與山谷信息的基準圖制備方法在光照變化情況下的邊緣特征與實時圖一致性更高，匹配性能更好。

圖9 同一區域不同時段可見光圖、融合圖像及其邊緣圖像Fig.9 Two vision maps at different time,a fusion map and their edge images

表3 圖9邊緣圖像局部匹配結果Table 3 Local matching results of edge image in Fig.9

為了驗證本文算法在實時圖與基準圖像在細節邊緣上差異較大時的魯棒性，選取10組秦嶺典型山地區域的可見光圖進行局部匹配實驗，在317×317大小的可見光圖中隨機選取60×60大小的圖像作為實時圖，對實時圖像施加不同的加性高斯噪聲(強度分別為1%～5%)，分別采用融合圖像與光照明暗圖作為基準圖，采用5種HD方法對基準圖與實時圖像進行局部匹配，對匹配正確率Accuracy與定位誤差PE兩個指標進行蒙特卡羅實驗統計。單次匹配結果的誤差Error表示匹配位置(x,y)與圖像實際位置(xT,yT)的歐式距離，定義匹配誤差超過5個像素為錯誤匹配，小于等于5個像素單位為正確匹配，ACC為正確匹配次數與所有匹配次數的比率，PE為所有正確匹配實驗的誤差平均值。

(19)

(20)

(21)

限于篇幅，本文給出部分圖像匹配的結果如表4、圖10和圖11所示，實驗的統計結果如圖12所示。由實驗數據可得出以下結論：

表4 不同基準圖像下幾種HD算法的匹配結果對比Table 4 Comparison of matching results of HD algorithms in different reference images

圖10 兩組典型地形的實驗圖像Fig.10 Two groups of experiments on typical terrain

1) 以Map_fusion(融合圖像)為基準圖像的匹配正確率均高于以Map_hillshade(光照明暗圖)為基準圖的匹配正確率，圖10中展示了兩組實時圖(未加噪/加噪)在Map_fusion中正確匹配而在Map_hillshade中錯誤匹配的現象，驗證了融合圖像有效的引入了山谷的顯著邊緣信息，與真實可見光圖像更具一致性。

2) 在不同的噪聲環境下，對于兩種不同的基準圖像，EHD方法的匹配正確率均高于其他4種HD方法。實驗結果表明EHD更有效的消除了出格點，提高了匹配面的峰值區分度，降低了誤匹配率。

3) 隨著噪聲的增大，實時圖與基準圖的邊緣差異增大，除EHD方法外的4種HD匹配方法的匹配正確率的指標下降顯著，而EHD保持了較高的匹配正確率，說明EHD具有更好的魯棒性。

4) 由于噪聲對邊緣位置帶來影響，各種HD算法的定位精度都隨噪聲的增大而下降，以Map_fusion為基準圖的匹配精度要高于以Map_hillshade為基準圖的匹配精度，驗證了融合圖像引入了有效、一致的邊緣信息。

圖11 實測圖像與仿真基準圖邊緣匹配結果Fig.11 Matching results between vision image and reference map

圖12 不同方法匹配正確率與定位誤差統計圖Fig.12 Graphs of correct rate and positioning error of different methods

5) 5種Hausdorff方法的匹配運行時間如表5所示。由于AM-HD采用了梯度濾波器帶來了更大的計算量，最為耗時；在邊緣點數量較少的情況下PHD、LTS-HD效率優于EHD，但由于PHD與LTS-HD計算中采用了排序算法，在邊緣點數量增大的情形下，相較EHD方法不具效率優勢；而EHD在邊緣點集數量較大的情形下效率僅次于M-HD。

表5 不同HD方法的平均運行時間Table 5 Average running time of different HD methods

5 結 論

1) 采用光照模型與山谷信息生成的融合圖像解決了不同光照環境引起的基準圖與實時圖在顯著邊緣處不致性的問題。

2) EHD方法克服了實時圖與基準圖由細節邊緣差異引起的匹配正確率下降的問題，對噪聲也具有良好的魯棒性。并且EHD方法保持了良好的效率與精度，對解決異構圖像邊緣匹配的問題也具有通用性。

3) 為構造低成本高精度的無人機景象匹配自主導航方法提供了一種新穎的思路。

4) 本文算法主要針對能見度好，云量少的拍攝環境，當出現大面積云層遮擋以及雨雪天氣的時，由于光源以及地表物體的反射率變得復雜，造成模型條件與實際情況差異大，使得匹配結果不佳，故算法仍存在一定的局限性，后續將進一步針對復雜光照條件的問題展開研究。

[1] ZHU Z S, SU A, LIU H B, et al. Vision navigation for aircrafts based on 3D reconstruction from real-time image sequences[J]. Science China Technological Sciences, 2015, 58(7): 1196-1208.

[2] MANCINI F, DUBBINI M, GATTELLI M, et al. Using unmanned aerial vehicles (UAV) for high-resolution reconstruction of topography:The structure from motion approach on coastal environments[J]. Remote Sensing, 2013, 5(12): 6880-6898.

[3] 崔平遠, 馮軍華, 朱圣英, 等. 基于三維地形匹配的月球軟著陸導航方法研究[J]. 宇航學報, 2011, 32(3): 470-476.

CUI P Y, FENG J H, ZHU S Y, et al. 3D terrain feature matching based navigation for lunar soft landing[J]. Journal of Astronautics, 2011, 32(3): 470-476 (in Chinese).

[4] 李立春, 苑云, 李由, 等. 三維地形不變性特征描述及其在地形匹配中的應用[J]. 航空學報, 2009, 30(11): 2143-2148.

LI L C, YUAN Y, LI Y, et al. Invariant feature vector description for 3D terrain and its application to terrain contour matching[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(11): 2143-2148 (in Chinese).

[5] LI L C, YU Q F, SHANG Y, et al. A new navigation approach of terrain contour matching based on 3-D terrain reconstruction from onboard image sequence[J]. Science China Technological Sciences, 2010, 53(5): 1176-1183.

[6] 馬旭, 程詠梅, 郝帥, 等. 面向未知區域深度測量的序列圖像稠密點特征生成算法[J]. 航空學報, 2015, 36(2): 596-604.

MA X, CHENG Y M, HAO S, et al. Dense point feature generation algorithm based on monocular sequence images for depth measurement of unknown zone[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(2): 596-604 (in Chinese).

[7] SIM D G, PARK R H, KIM R C, et al. Integrated position estimation using aerial image sequences[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002, 24(1): 1-18.

[8] 于秋則, 程輝, 柳健, 等. 基于改進Hausdorff 測度和遺傳算法的SAR圖像與光學圖像匹配[J]. 宇航學報, 2006, 27(1): 130-134.

YU Q Z, CHENG H, LIU J, et al. Matching SAR image to optical image using modified Hausdorff distance and genetic algorithms[J]. Journal of Astronautics, 2006, 27(1): 130-134 (in Chinese).

[9] 陳天澤， 李燕. 一種高性能SAR圖像邊緣點特征匹配方法[J]. 自動化學報, 2013, 39(12): 2051-2063.

CHEN T Z, LI Y. A high performance edge point feature match method of SAR images[J]. Acta Automatica Sinica, 2013, 39(12): 2051-2063 (in Chinese).

[10] WANG T, CELIK K, SOMANI A K. Characterization of mountain drainage patterns for GPS-denied UAS navigation augmentation[J]. Machine Vision and Applications, 2016, 27(1): 87-101.

[11] WANG L, LIU H. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling[J]. International Journal of Geographical Information Science, 2006, 20(2): 193-213.

[12] HUTTENLOCHER D P, KLANDERMAN G A, RUCKLIDGE W J. Comparing images using the Hausdorff distance[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1993, 15(9): 850-863.

[13] DUBUISSON M P, JAIN A K. A modified Hausdorff distance for object matching[C]∥Proceedings of 12th International Conference on Pattern Recognition, 1994: 566-568.

[14] SIM D G, KWON O K, PARK R H. Object matching algorithms using robust Hausdorff distance measures[J]. IEEE Transactions on Image Processing, 1999, 8(3): 425-429.

[15] KWON O K, KIM D, SUH J W. Accurate M-Hausdorff distance similarity combining distance orientation for matching multi-modal sensor images[J]. Pattern Recognition Letters, 2011, 32(7): 903-909.

[16] LOPEZ A M, LLORET D, SERRAT J, et al. Multilocal creaseness based on the level-set extrinsic curvature[J]. Computer Vision and Image Understanding, 2000, 77(2): 111-144.

[17] CANNY J. A computational approach to edge detection[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1986(6): 679-698.

[18] CHEN H F, BELHUMEUR P N, JACOBS D W. In search of illumination invariants[C]∥IEEE Conference on Computer Vision and Pattern Recognition,2000: 254-261.

[19] HOM B. Robot vision[M]. Beijing: China Youth Press, 2014: 246-247.

[20] BONHAM-CARTER G F. Geographic information systems for geoscientists: Modelling with GIS[M]. Amsterdam: Elsevier, 2014: 295-297.

[21] JENSON S K, DOMINGUE J O. Extracting topographic structure from digital elevation data for geographic information system analysis[J]. Photogrammetric Engineering and Remote Sensing, 1988,54(11): 1593-1600.

Arobustscenematchingmethodformountainousregionswithilluminationvariation

WANGHuaxia1,2,*,CHENGYongmei1,LIUNan1

1.CollegeofAutomation,NorthwesternPolytechnicalUniversity,Xi’an710072,China2.CollegeofElectronicandInformationEngineering,TaiyuanUniversityofScienceandTechnology,Taiyuan030024,China

Thispaperpresentsarobustscenematchingnavigationmethodformountainousregionswithilluminationvariation.Inordertosolvetheproblemoftheinconsistencybetweenthereferencemapandtherealmapcausedbydifferentlightingconditionsinmountainousareas,thereferenceimageisobtainedbyfusingthehillshademap,generatedfromelevationdatausingilluminationmodel,andthevalleysaliencymap,generatedfromwaterflowdatacollectedbyDigitalElevationModel(DEM).Toavoidthemismatchcausedbythedifferenceofedgebetweenthereferenceandtherealmaps,aHausdorffdistanceedgematchingalgorithmbasedonmorphologicalconstraintsisproposed.TheexperimentalresultsofLANDSATimageandASTERDEMdatashowthattheproposedmethodhashighaccuracyandgoodrobustness.

scenematchingnavigation;illuminationmodel;significantedge;Hausdorffdistance;imagefusion

2017-01-04；Revised2017-02-06；Accepted2017-03-13;Publishedonline2017-03-201514

URL：http://hkxb.buaa.edu.cn/CN/html/20171019.html

s：NationalNaturalScienceFoundationofChina(61603364);Xi’anScienceandTechnologyProject(CXY1436(9)，CXY1350(2))

.E-mailwanghx_sam@163.com

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.321101

V249.32

A

1000-6893(2017)10-321101-13

2017-01-04；退修日期2017-02-06；錄用日期2017-03-13；< class="emphasis_bold">網絡出版時間

時間：2017-03-201514

http://hkxb.buaa.edu.cn/CN/html/20171019.html

國家自然科學基金(61603364)； 西安市科技計劃項目(CXY1436(9)，CXY1350(2)).

*

.E-mailwanghx_sam@163.com

王華夏，程詠梅，劉楠．面向山地區域光照變化下的魯棒景象匹配方法J． 航空學報,2017,38(10):321101.WANGHX,CHENGYM,LIUN.ArobustscenematchingmethodformountainousregionswithilluminationvariationJ.ActaAeronauticaetAstronauticaSinica,2017,38(10):321101.

(責任編輯：蘇磊)