王華夏，程詠梅，劉楠

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

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

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

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

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

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

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

山區(qū)景象的邊緣特征豐富，通過模型構(gòu)造的基準(zhǔn)圖像與實時圖像在圖像的細節(jié)邊緣上存在差異，如何在這種差異存在的情況下實現(xiàn)正確的局部匹配是景象匹配的另一個問題。由于霍斯多夫距離(Hausdorff Distance，HD)具有對邊緣局部形變不敏感、不要求所有匹配點嚴(yán)格對應(yīng)、計算快捷的優(yōu)點，因此在圖像匹配中得到了廣泛的應(yīng)用，但HD極易受到出格點的影響，當(dāng)邊緣圖像受到噪聲污染或存在遮擋等情況時，易造成誤匹配。文獻[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的基礎(chǔ)上提出了AM-HD (Accurate M-HD)距離，引入邊緣的方向信息對M-HD進行權(quán)重調(diào)整，雖然以上方法在處理遮擋與退化圖像的匹配問題時已具有較好的性能，但當(dāng)基準(zhǔn)圖與實時圖的細節(jié)邊緣差異較大、實時圖存在噪聲的情況下，誤匹配率較高。

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

1 光照變換下的魯棒基準(zhǔn)圖制備方法

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

MatchmethodNormalizedcorrelationHDHDMapstomatchFig．1(a)，F(xiàn)ig．1(d)Fig．1(b)，F(xiàn)ig．1(e)Fig．1(c)，F(xiàn)ig．1(f)Accuracy/%43824

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

1.1 光照明暗圖

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

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

L=Ecosγ

(1)

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

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

(2)

式中：θp為地形坡度角，是地形切平面與水平面的夾角;φo為地表法線方位角，為ON在水平面投影向量與正北方向的夾角，θp、φo與DEM數(shù)據(jù)相關(guān)，通過DEM象元及其鄰域象元的關(guān)系求?。沪誥為太陽方位角;θz為太陽天頂角;θH為太陽高度角，其中θz+θH=π/2，θH、φa與太陽位置和當(dāng)?shù)匚恢糜嘘P(guān)，通過時間和當(dāng)?shù)亟?jīng)緯可推算得出。

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

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

(3)

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

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

1.2 山谷顯著圖

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

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

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

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

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

MD(i,j)=

(4)

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

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

flag(i,j,Sk)=

(5)

(6)

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

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

Map_valley(i,j)=

(7)

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

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

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

Map_fusion=

W·Map_vally+(1-W)·Map_hillshade

(8)

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

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

(9)

圖5 不同融合權(quán)重下融合圖像與真實圖像的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)不僅具有光照產(chǎn)生的明暗邊緣，還包含較為完善的山谷顯著邊緣信息，與可見光邊緣圖7(e)更具一致性。在317×317大小的圖7(e)中隨機截取60×60大小區(qū)域作為實時邊緣圖像，采用4種HD方法分別與圖7(f)、圖7(g)、圖7(h)進行局部匹配，設(shè)定定位誤差小于5個像素單位為正確匹配，對匹配正確率進行蒙特卡羅統(tǒng)計實驗，實驗結(jié)果如表2所示，以圖7(h)為基準(zhǔn)圖的局部匹配正確率最高，而以圖7(f)為基準(zhǔn)圖的匹配正確率高于以圖7(g)為基準(zhǔn)圖的匹配正確率，由此說明：光照明暗圖產(chǎn)生的邊緣信息在邊緣匹配中做了主要貢獻；由于山谷邊緣信息稀疏，使得山谷顯著邊緣作為基準(zhǔn)圖的匹配正確率很低；而采用融合圖像的邊緣圖像作為基準(zhǔn)圖得到了提升匹配正確率的效果；同時也側(cè)面驗證了光照明暗圖在融合圖像中起的作用較大，應(yīng)給予更高的權(quán)值。

圖6 不同融合權(quán)重對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邊緣圖像局部匹配結(jié)果Table 2 Local matching results of edge images in Fig.7

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

2 EHD邊緣匹配算法

設(shè)P、Q為兩個邊緣圖像的邊緣點集為

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

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

(10)

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

形態(tài)學(xué)中的膨脹操作是將結(jié)構(gòu)元素B平移后與原始點集E求交集，所獲取的非空點集Y為膨脹操作構(gòu)成的點集。膨脹操作的定義為

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

(11)

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

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

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

(12)

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

(13)

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

(14)

最終的EHD距離為

HEHD(P,Q)=

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

(15)

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

(16)

采用HD距離的圖像匹配導(dǎo)航的問題描述為

(17)

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

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

(18)

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

目前已有的HD方法只能給出經(jīng)驗參數(shù)來設(shè)置參與計算的邊緣點比例，無法進行自適應(yīng)調(diào)整，本文采用膨脹約束的方法給出動態(tài)調(diào)整HD的方法，以提高HD對細節(jié)邊緣差異的魯棒性。

3 魯棒景象匹配方法

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

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

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

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

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

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

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

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

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

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

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

4 景象匹配實驗

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

為了驗證在光照差異情況下采用融合圖像作為基準(zhǔn)圖和采用固定時間圖像作為基準(zhǔn)圖的匹配性能差異，選取秦嶺同一地區(qū)不同時間的兩幅可見光圖像進行實驗， 如圖9所示， 圖9(a)拍攝時間為2015/12/17，當(dāng)?shù)靥柛叨冉菫?0°，太陽方位角158°，圖9(b)拍攝時間為2016/7/28，當(dāng)?shù)靥柛叨冉菫?4°，太陽方位角118°，圖9(c)為依據(jù)圖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)分別作為基準(zhǔn)圖，采用5種HD方法進行局部匹配，局部匹配結(jié)果如表3所示。由圖可見，圖9(d)、圖9(e)在明暗邊緣處的一致性較差，圖9(d)、圖9(f)在明暗邊緣與河谷邊緣處的一致性均較好；以圖9(f)作為基準(zhǔn)圖的5種HD方法的匹配正確率均高于以圖9(e)作為基準(zhǔn)圖的匹配正確率。驗證了本文提出的基于光照模型與山谷信息的基準(zhǔn)圖制備方法在光照變化情況下的邊緣特征與實時圖一致性更高，匹配性能更好。

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

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

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

(19)

(20)

(21)

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

表4 不同基準(zhǔn)圖像下幾種HD算法的匹配結(jié)果對比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(融合圖像)為基準(zhǔn)圖像的匹配正確率均高于以Map_hillshade(光照明暗圖)為基準(zhǔn)圖的匹配正確率，圖10中展示了兩組實時圖(未加噪/加噪)在Map_fusion中正確匹配而在Map_hillshade中錯誤匹配的現(xiàn)象，驗證了融合圖像有效的引入了山谷的顯著邊緣信息，與真實可見光圖像更具一致性。

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

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

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

圖11 實測圖像與仿真基準(zhǔn)圖邊緣匹配結(jié)果Fig.11 Matching results between vision image and reference map

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

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

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

5 結(jié) 論

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

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

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

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

[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] 崔平遠, 馮軍華, 朱圣英, 等. 基于三維地形匹配的月球軟著陸導(dǎo)航方法研究[J]. 宇航學(xué)報, 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] 李立春, 苑云, 李由, 等. 三維地形不變性特征描述及其在地形匹配中的應(yīng)用[J]. 航空學(xué)報, 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] 馬旭, 程詠梅, 郝帥, 等. 面向未知區(qū)域深度測量的序列圖像稠密點特征生成算法[J]. 航空學(xué)報, 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圖像與光學(xué)圖像匹配[J]. 宇航學(xué)報, 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]. 自動化學(xué)報, 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">網(wǎng)絡(luò)出版時間

時間：2017-03-201514

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

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

*

.E-mailwanghx_sam@163.com

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

(責(zé)任編輯：蘇磊)