一種基于協(xié)方差矩陣的快速抗旋轉(zhuǎn)景象匹配算法

符艷軍, 孫開鋒, 劉軍蘭

(1 西安外事學(xué)院工學(xué)院，西安 710077；2 西安精密機械研究所，西安 710077)

0 引言

無人機在飛行過程中，受地形地貌變化、平臺姿態(tài)變化以及慣性導(dǎo)航系統(tǒng)累加誤差等因素影響，所拍實測圖與基準圖相比可能存在遮擋、噪聲以及旋轉(zhuǎn)變化，為了能夠正確地進行導(dǎo)航定位，要求所采用的景象匹配算法能夠?qū)崪y圖的這些變化具有魯棒性。目前，基于歸一化積相關(guān)(normalized cross-correlation, NCC)的灰度相關(guān)匹配方法已經(jīng)廣泛應(yīng)用在目標跟蹤以及景象匹配中[1-5]，該方法對圖像發(fā)生的小變化(包括噪聲、灰度變化及2°內(nèi)的旋轉(zhuǎn))具有一定的魯棒性，但另一方面，NCC作為相似性測度，在真實匹配位置其值接近于1，但對于非匹配區(qū)域其取值往往也接近于1[6-7]，也就是說NCC作為相似性測度缺乏很好的區(qū)分度，這樣當基準圖上存在多個相似區(qū)域且實測圖存在局部變化時往往會導(dǎo)致錯誤的定位結(jié)果。協(xié)方差矩陣最小特征值是一種能夠?qū)D像局部變化做出顯著響應(yīng)的相似性度量，可以較好的區(qū)分出真實匹配位置和非匹配位置[8]，在目標識別、目標跟蹤中取得了較好的結(jié)果[9-10]，但作為一種灰度相關(guān)匹配算法仍面臨計算量大、對大角度旋轉(zhuǎn)敏感等問題。為了減少匹配耗時，文獻[11]利用歸一化協(xié)方差矩陣特性，以5個數(shù)學(xué)條件是否滿足作為是否匹配的判斷依據(jù)，避開了協(xié)方差矩陣特征值的計算，減少了計算量，但該匹配依據(jù)對旋轉(zhuǎn)非常敏感，不適合應(yīng)用于景象匹配場合。文獻[12]借助GPU通用計算技術(shù)，結(jié)合CUDA編程模型，通過并行計算實現(xiàn)圖像協(xié)方差矩陣高效運算，該方法對計算的硬件條件有要求，且主要適用于多通道彩色圖像，對景象匹配用的遙感灰度圖像不能發(fā)揮其優(yōu)越性。總之，已有的基于協(xié)方差矩陣的快速匹配算法或者通過并行計算或通過修改相似性測度本身達到減少匹配耗時的目的，這些改進方法不太適用于無人機景象匹配場合。

文中以協(xié)方差矩陣最小特征值為相似性度量，從實時性和抗旋轉(zhuǎn)兩方面著手進行改進，提出了一種快速抗旋轉(zhuǎn)匹配算法。

1 協(xié)方差矩陣及其特征值

圖像協(xié)方差矩陣的特征值可用于提取對應(yīng)圖中灰度分布的形狀信息。較大特征值代表兩匹配圖像灰度對應(yīng)圖中沿形狀長軸的方差，較小特征值代表沿形狀短軸的方差，基于最小特征值的相似性度量能夠顯著響應(yīng)兩幅被比較圖像之間的局部偏差[13]。

設(shè)實測圖尺寸為m×n，w(i,j)為實測圖上點(i,j)處的灰度值，f(i,j)為與實測圖尺寸相同的基準子圖上點(i,j)處的灰度值，則實測圖與基準子圖的協(xié)方差矩陣C可表示為：

(1)

其中：

式(1)所示的協(xié)方差矩陣C是一個正定的、對稱的矩陣，有兩個特征值，其中較大特征值表示沿著形狀長軸的方差，較小特征值表示沿著形狀短軸的方差。對于相互匹配的圖像，理想情況下，沿著形狀短軸方向的方差應(yīng)該為0，因此，較小特征值可以用作兩幅圖像的相似性度量。記較小特征值為λS，計算公式為：

(2)

在景象匹配過程中，式(1)中基于實測圖的c22只需計算一次，而基于基準子圖的c11以及與實測圖和基準子圖都相關(guān)的c12則需要在每個匹配位置進行計算，若采用傳統(tǒng)方法以在線方式計算每個匹配位置的c11和c12，則很難滿足無人機景象匹配導(dǎo)航對實時性的要求，減少每個匹配位置c11和c12的計算量是提高景象匹配導(dǎo)航系統(tǒng)實時性的有效途徑。

2 基準子圖協(xié)方差矩陣的快速算法

根據(jù)協(xié)方差矩陣的計算特點以及匹配過程中相鄰基準子圖協(xié)方差的關(guān)系，通過在離線方式下預(yù)先創(chuàng)建基準圖累加和矩陣，減少匹配過程中每個匹配位置協(xié)方差矩陣的在線計算量來提高匹配實時性。

2.1 創(chuàng)建基準圖累加和矩陣

設(shè)基準圖尺寸為M×N,實測圖尺寸為m×n，以基準圖左上頂點為(1,1)坐標點，建立如圖1所示坐標系，點(u,v)(u≥m,v≥n)表示待匹配點，以點(u-m+1,v-n+1)及點(u,v)為對角線的小矩形表示與實測圖尺寸相等的基準子圖。

圖1 基準圖及基準子圖

匹配前，首先在離線狀態(tài)下建立基準圖的2個M×N的累加和矩陣S，T，可以看出，只需在基準圖上分別遍歷一次即可求出累加和矩陣S，T中的所有元素，即

(3)

(4)

2.2 任一點匹配位置協(xié)方差矩陣的快速求解

匹配過程中，對于圖1中任一點匹配位置(u,v)處大小為m×n的基準子圖有：

(5)

(6)

根據(jù)式(5)和式(6)，對任一點匹配位置的協(xié)方差矩陣，c11,c12的計算可以簡化為：

(7)

(8)

從上述推導(dǎo)過程可以看出，在每一個匹配位置，基于預(yù)先建立的累加和矩陣S和T，c11的在線計算只需7次加減運算和3次乘法運算即可精確求出，且計算量不隨基準子圖尺寸的增大而增加，而傳統(tǒng)方法則需要(mn+3)次乘法和(2mn+1)次加法，計算量隨實測圖尺寸的增加而增加；c12(或c21)的計算中，基準子圖矩陣元素的和tlocal(u,v)也可以由式(4)直接算出，相比傳統(tǒng)方法加法運算次數(shù)減少了(mn-3)次。

3 基準子圖虛擬圓形窗口的快速創(chuàng)建

在景象匹配過程中，獲取的實測圖通常為矩形圖像。矩形圖像在旋轉(zhuǎn)前后會產(chǎn)生裁剪誤差，造成實測圖和基準子圖所包含的圖像內(nèi)容有所差別，從而影響匹配結(jié)果，這也是造成很多匹配算法對旋轉(zhuǎn)敏感的因素之一。圓形窗口可以減少因旋轉(zhuǎn)而產(chǎn)生的圖像內(nèi)容的變化，為此，文中采用虛擬圓形窗口進行匹配。為了能夠快速的在每個待匹配位置定位出虛擬圓窗口，采用基于索引模板的快速虛擬圓定位算法[14]，具體步驟為：

1) 根據(jù)所拍實測圖大小，從其左上角點截取最大方形圖像作為匹配用實測圖，計算實測圖中每個像素點到中心的距離以定位出實測圖的內(nèi)切圓；

2) 將實測圖內(nèi)切圓外的所有像素點的灰度值賦為0，將這些0值點在原實測圖中的坐標保存在索引表Index中；

3) 在基準圖上截取與實測圖大小相同的方形子圖后，首先計算基準子圖中Index表對應(yīng)位置所有像素的灰度值的和Eind以及灰度值的平方和Dind,再將這些位置對應(yīng)像素的灰度值賦為0即得到虛擬的圓形子窗口。

圖2所示為采用虛擬圓索引模板前后圖像變化情況，可以看出，圖像發(fā)生旋轉(zhuǎn)后，如果直接使用矩形窗口進行匹配，則旋轉(zhuǎn)前后圖像內(nèi)容相差較大(如圖2(a)和圖2(b)所示)，而采用虛擬圓索引模板后(如圖2(c)和圖2(d)所示)，圖像旋轉(zhuǎn)前后圖像內(nèi)容變化不大。

圖2 虛擬圓形子圖

4 景象匹配算法流程

基于無人機景象匹配導(dǎo)航系統(tǒng)對匹配算法實時性和魯棒性的要求，文中算法流程如圖3所示。具體實現(xiàn)步驟為：

圖3 文中方法流程圖

1) 根據(jù)實測圖尺寸大小建立虛擬圓索引表Index并應(yīng)用于實測圖，再計算c22；

2) 根據(jù)實測圖尺寸大小，在離線方式下，依據(jù)s(u,v)，t(u,v)建立基于基準圖的累加和矩陣S，T；

3) 從基準圖上第1個匹配位置開始遍歷所有待匹配點，以在線方式分別計算每個匹配位置對應(yīng)的圓形匹配窗口的協(xié)方差矩陣的c12(或c21)，c11以及λS；

4) 找出3)中最小的λS所對應(yīng)的匹配點作為最終匹配位置。

從上述步驟可以看出，在每一個匹配位置，協(xié)方差矩陣中需要在線計算的元素主要集中在c12或c21上。設(shè)實測圖尺寸為m×n，則在匹配過程中，文中方法在每個搜索位置創(chuàng)建協(xié)方差矩陣C的在線計算量(包括虛擬圓索引處理計算量)約為(1.2m2+3)次乘法和(1.2m2+7)次加法，而傳統(tǒng)方法的在線計算量約為(2m2+2)次乘法和(3m2-1)次加法，兩種方法都是對協(xié)方差矩陣最小特征值λS的精確計算；另一方面，采用虛擬圓形窗口匹配，可以提高匹配方法的抗旋轉(zhuǎn)性能。

5 實驗仿真

為了驗證文中方法的有效性，以Matlab R2017為仿真平臺，分別在可見光圖像(如圖4(a)所示)和SAR圖像(如圖5(a)所示)上進行了兩組仿真實驗。為了方便比對，每組仿真中的不同方法采用相同的基準圖和實測圖。

圖4 可見光圖像匹配

圖5 SAR圖像匹配

第一組仿真實驗以協(xié)方差矩陣最小特征值λS為相似性測度對3種方法進行了仿真。方法1是傳統(tǒng)匹配方法，方法2是單純采用文中提出的累加和矩陣的匹配方法，方法3即文中第5節(jié)所述的綜合采用累加和矩陣和虛擬圓形窗口的匹配方法。其中對方法2進行仿真的目的主要是為了突出累加和矩陣策略在減少匹配耗時方面的顯著性。

第一組仿真實驗的結(jié)果如表1所示，可以看出，方法2采用累加和矩陣策略后的匹配耗時遠小于傳統(tǒng)方法1，且實測圖尺寸越大，耗時減少的愈加明顯，而匹配精度與方法1完全相同，但在實測圖存在5°～10°旋轉(zhuǎn)時，方法1和方法2的失配率都較高；方法3在實測圖存在5°～10°旋轉(zhuǎn)時其匹配誤差不超過5個像素，但相比方法2因多了虛擬圓處理的環(huán)節(jié)其匹配耗時較方法2要長，但相比傳統(tǒng)方法1在耗時方面仍有明顯優(yōu)勢。

表1 基于最小特征值λS的3種匹配方法比較

第二組仿真實驗通過將文中第5節(jié)所述方案通過與傳統(tǒng)NCC方法進行比較，一方面是驗證文中方法對景象匹配導(dǎo)航過程中實測圖可能發(fā)生的畸變(包括噪聲、遮擋、旋轉(zhuǎn)、比例變化)的適應(yīng)性能。第二組仿真實驗結(jié)果如表2及圖4、圖5所示。

表2 文中方法與傳統(tǒng)基于NCC匹配方法比較

圖4中：圖(a)為原始可見光基準圖，圖(b)為采用文中方法時實測圖1、實測圖2及實測圖3在基準圖(a)上的匹配結(jié)果，圖(c)、圖(d)、圖(e)分別為實測圖1、實測圖2、實測圖3及其虛擬圓處理后的結(jié)果圖，其中，圖(c)為存在10°旋轉(zhuǎn)實測圖，圖(d)為存在10°旋轉(zhuǎn)、遮擋面積5×5像素的實測圖，圖(e)為存在10°旋轉(zhuǎn)、比例變化為1.05倍、遮擋面積5×5像素、高斯噪聲0.01的實測圖。圖(f)、圖(g)、圖(h)分別為實測圖1、實測圖2以及實測圖3采用NCC方法時的匹配相關(guān)面。

圖5中：圖(a)為原始SAR基準圖，圖(b)為采用文中方法時實測圖1、實測圖2及實測圖3在基準圖(a)上的匹配結(jié)果，圖(c)、圖(d)、圖(e)分別為實測圖1、實測圖2、實測圖3及其虛擬圓處理后的結(jié)果圖，其中，圖(c)為發(fā)生10°旋轉(zhuǎn)、比例變化0.9倍實測圖，圖(d)為發(fā)生10°旋轉(zhuǎn)、比例變化1.1倍、遮擋面積為10×10個像素的實測圖，圖(e)為發(fā)生10°旋轉(zhuǎn)，比例變化為1.05倍、高斯噪聲0.01的實測圖，圖(f)、圖(g)、圖(h)分別為實測圖1、實測圖2以及實測圖3采用NCC方法時的匹配相關(guān)面。

從表2以及圖4、圖5的匹配結(jié)果可以看出，文中方法與傳統(tǒng)NCC方法對噪聲、遮擋及比例變化方面的適應(yīng)性能相當，都對噪聲(方差在0.01內(nèi))、小范圍遮擋(10像素×10像素內(nèi))及小比例變化(正負0.05倍內(nèi))有一定的魯棒性，但在抗旋轉(zhuǎn)方面，當實測圖存在5°～10°左右的旋轉(zhuǎn)時，即使此時的實測圖不存在噪聲、沒有比例變化及遮擋，NCC方法的失配率仍然很高，此時其匹配相關(guān)面上會出現(xiàn)多個高度相近的峰值，而文中方法在一定誤差范圍內(nèi)(小于5個像素)可以正確匹配。

6 結(jié)論

1) 基于協(xié)方差矩陣最小特征值的相似性度量，能夠?qū)D像局部變化做出顯著響應(yīng)，同時對待匹配圖像間存在的噪聲、小范圍遮擋、小比例變化等有一定的魯棒性，比較適合無人機景象匹配導(dǎo)航的應(yīng)用場景。

2) 針對協(xié)方差矩陣計算量大、實時性差的問題，利用協(xié)方差矩陣計算特點，提出基于累加和矩陣減少匹配在線計算量的方法來改善實時性，仿真結(jié)果表明，匹配耗時明顯減少，且隨著基準圖和實測圖尺寸的增大，效果更加顯著。

3) 針對協(xié)方差矩陣最小特征值對圖像旋轉(zhuǎn)較為敏感的問題，提出采用虛擬圓形窗口進行匹配。仿真結(jié)果表明，采用虛擬圓型窗口匹配，能夠抵抗實測圖5°～10°度左右的旋轉(zhuǎn)變化，可以滿足中等精度慣性導(dǎo)航系統(tǒng)對系統(tǒng)漂移的補償。