多尺度PCA-HOG遙感異源圖像匹配算法 *

韓松來，王鈺婕，王 星，羅世彬，董 晶

(1. 中南大學 航空航天學院， 湖南 長沙 410083； 2. 復雜系統控制與智能協同技術重點實驗室， 北京 100074；3. 國防科技大學 空天科學學院， 湖南 長沙 410073)

由于成像的物理機制不同，異源圖像可反映出地物的不同特性，并通過各類信息的互補獲得更加完整的圖像信息。遙感異源圖像匹配是各類遙感信息處理和視覺導航應用中的關鍵技術之一[1-5]。一方面，由于成像傳感器的差異，異源圖像之間存在復雜的非線性灰度畸變，這會嚴重降低相同對象在不同圖像間的相似性；另一方面，遙感圖像可能因傳感器或傳輸引入強噪聲干擾，對圖像信息產生嚴重破壞，以光學影像和合成孔徑雷達(Synthetic Aperture Radar，SAR)影像為例(如圖1所示)。因此，盡管圖像匹配經過了長期廣泛的研究，遙感異源圖像匹配仍然是一項極具挑戰的難題。

1 相關研究

(a) 光學影像(a) Visible image (b) SAR影像(b) SAR image圖1 噪聲干擾和非線性灰度畸變Fig.1 Noise corruption and nonlinear gray distortion

圖像匹配可以分為基于局部特征的匹配方法和基于模板的匹配方法。基于局部特征的匹配方法首先對圖像提取局部特征，然后生成局部特征描述，再進行特征匹配，最后根據特征對應關系計算圖像的全局變換得到匹配結果。由于有效地利用了特征的不變性，并能適應尺度和視角變化，基于局部特征的圖像匹配方法在計算視覺領域中得到了廣泛應用[6-8]。然而，遙感異源圖像之間存在強噪聲干擾和復雜的非線性灰度畸變，這會導致局部特征發生嚴重變化，從而無法在圖像間檢測到足夠多的重復局部特征，并且變化也使得正確匹配局部特征變得非常困難[9-11]。尺度比不變特征變換算法(Scale Invariant Feature Transform, SIFT)是一種經典的局部特征匹配算法。圖2使用SIFT對遙感SAR和可見光圖像進行匹配，可以看到幾乎所有的局部特征都匹配錯誤。

圖2 基于SIFT特征的遙感SAR(左)與可見光(右)圖像匹配Fig.2 Remote sensing SAR (left) and visible light (right) image matching based on SIFT feature

基于模板的匹配方法也可以稱為相關匹配方法，它將模板圖像與基準圖像上每個候選位置區域的圖像進行比較，并選擇具有最大相似性系數或最小差別系數的對應位置作為匹配結果。直接基于灰度的模板匹配很難有效適應復雜的非線性灰度畸變，但是利用結構特征描述圖像信息可以改善這個問題，因為結構的幾何信息能對灰度畸變保持不變。并且，基于模板的匹配方法可以利用密集描述(每個像素都提取)最大限度地提取圖像的所有結構特征信息，而基于局部特征的匹配方法往往是基于稀疏的特征對應，因此前者往往具有更強的可靠性。相關研究中的比較測試結果表明，基于模板的匹配方法比基于局部特征的匹配方法能更有效地適應強噪聲干擾和復雜的非線性灰度畸變，從而能更有效地匹配遙感異源圖像[9-11]。

為了適應異源圖像間的復雜非線性灰度畸變，一些模板匹配算法對相似性或差別度量準則進行了改進[12-13]。歸一化互相關(Normalized Cross Correlation, NCC)是一種圖像匹配常用的相似性度量方法，它能對圖像間的線性灰度畸變保持不變，并且因為圖像間的單調非線性灰度變化通常是局部線性的，對這種情況，其也能有效匹配圖像[14]。然而，當圖像間發生非單調非線性的灰度畸變時，NCC可能無法正確匹配圖像[15]。對此，Hel-Or等提出了一種基于色調映射(Matching by Tone Mapping, MTM)的圖像匹配方法[14]，該方法可以看作是NCC的非線性范化，而NCC是該方法的線性簡化版。MTM匹配方法的計算時間與NCC接近，但是能有效適應圖像間的非單調非線性灰度畸變。MTM仍然需要假設圖像間的灰度畸變滿足函數映射關系，但是異源圖像間的灰度畸變通常不滿足這種關系?；バ畔?Mutual Information, MI)是一種常用的異源圖像匹配相似度度量方法。MI通過計算兩幅圖像之間的統計相關性來確定相似性，它不需要圖像間的灰度畸變滿足函數映射關系[16]。根據相關研究，MI比NCC和MTM具有更好的異源圖像匹配性能[13]。但是在實際使用中，MI匹配方法往往存在匹配搜索計算量過大和參數敏感問題[12]。

僅僅改進匹配的測量準則并直接對灰度圖像進行匹配，效果往往不太理想，這樣忽略了圖像的紋理結構信息。鑒于此，一些模板匹配算法先使用結構特征描述提取圖像中的紋理結構信息得到本征圖像，然后再對本征圖像進行匹配。

方向梯度直方圖(Histogram of Oriented Gradients，HOG)利用梯度的方向和歸一化梯度強度來描述圖像的結構特征[17]，它可以很好地適應光照和對比度變化，在圖像匹配中經常使用。但是，HOG描述基于圖像梯度，因此對圖像噪聲比較敏感。局部自相似(Local Self Similarity, LSS)描述[18]和稠密自適應自相關(Dense Adaptive Self-Correlation, DASC)描述[19]都利用自相關系數的變化來提取圖像的結構特征，可以有效適應非線性灰度畸變，但是因為相關系數對噪聲也非常敏感，這兩種方法都很難用于遙感異源圖像匹配。Kovesi[20]提出的相位一致性模型可以捕獲圖像的結構幅值，該模型對灰度畸變能較好地保持不變性，并且不容易受噪聲影響，但是結構幅值沒有方向性，因此它包含的圖像結構信息比較有限，這不利于匹配中區分重復模式。因此，Ye等擴展了相位一致性模型，提出了相位一致性直方圖[21](Histogram of Orientated Phase Consistency, HOPC)。他們利用Log-Gabor奇對稱小波計算相位一致性的方向，并根據相位一致性的方向和幅值構造描述子。他們的實驗表明：相比于MTM、MI和NCC，HOPC能更可靠地匹配遙感異源圖像。

基于模板的匹配方法還存在搜索過程計算量太大的問題。全搜索的模板匹配是一種逐像素的匹配，特別是對復雜的全局變換(如透視變換、仿射變換)，搜索時要考慮平移、縮放、旋轉等多個維度。當搜索空間維度加大，全搜索計算復雜性會呈指數級增長。而基于局部特征的匹配算法只要在特征集里搜索，特征描述的維數遠遠低于圖像維數，并且不需要在圖像變換空間搜索，所以計算速度優勢明顯。但是如果圖像在匹配前經過校正，消除了縮放、旋轉、仿射、透視等變形，則模板匹配可以只做平移搜索，從而能大大減少計算量，經過搜索優化的模板匹配算法的計算速度往往快于基于特征的圖像匹配算法[22]。但是對于一些實時應用，平移搜索的模板匹配計算量仍然需要進一步優化。

針對遙感異源圖像匹配存在的問題，本文提出了一種基于主成分分析(Principal Component Analysis，PCA)增強的HOG特征描述的快速圖像匹配算法(PCA-HOG)，并在實驗部分對該算法的匹配準確率與計算效率進行了驗證。

2 基于PCA-HOG描述的結構特征提取

2.1 傳統HOG算子描述

HOG[17]是目標檢測和圖像匹配中常用的特征描述子。HOG通過統計局部梯度強度和方向構成直方圖，并對圖像的結構特征進行描述。HOG特征與圖像的整體灰度變化相獨立，因此對于異源圖像存在的非線性輻射畸變，HOG能保持良好的不變性。HOG特征的提取主要分成以下幾個步驟：

1)使用Sobel或Laplacian梯度算子與圖形進行卷積，得梯度圖像。

2)按照圖3，將方向分成8個bin，并對梯度按照方向進行統計，得到8維的直方圖。

3)對直方圖進行歸一化處理，降低灰度畸變引起的梯度強度變化。

4)一個局部特征區域可以分成若干個block，每個block又可以分成若干個patch，每個patch可以按照步驟2統計得到一個8維直方圖，然后進行拼接得到最終的HOG描述。若一個block有16個patch，那么它對應的HOG描述就是128維。

可見HOG描述是基于局部梯度而形成的，而梯度方向對于噪聲干擾十分敏感[23]。由于HOG描述是按照梯度方向進行統計得到，它的抗噪聲能力并不強，對遙感異源圖像的匹配正確率不高。

圖3 HOG的方向劃分Fig.3 Orientation division of HOG

2.2 PCA增強的HOG特征描述

為了使HOG特征描述更好地適應噪聲，可利用PCA對圖像的梯度方向進行增強，然后再按增強后的梯度方向統計得到HOG描述。本文提出的PCA-HOG描述的計算主要分成3個步驟：首先，使用 Sobel算子計算圖像梯度Gx與Gy，然后對梯度Gx、Gy進行多尺度PCA增強計算，最后按照梯度方向統計計算HOG描述。其中第一個步驟和最后一個步驟與傳統HOG算子的計算沒有區別。PCA增強計算可以分成基于奇異值分解(Singular Value Decomposition，SVD)和基于梯度求和的兩種方法。

2.2.1 基于SVD的 PCA計算

PCA常被用于計算給定數據集的優勢向量。因此，可以利用PCA對圖像局部梯度處理，計算出局部梯度主向量，其方向即是局部主方向。主方向比原來的梯度方向更加穩定，不容易受到噪聲干擾。圖4(a)給出了直接利用Sobel算子計算的梯度方向，圖4(b)給出了使用SVD計算的PCA提取的主方向，可見在噪聲干擾下，主方向能更加準確地提取圖像中結構紋理的方向。

(a) Sobel直接提取的梯度方向(a) Orientation maps estimated with Sobel

(b) SVD的PCA提取的梯度方向(b) Orientation maps estimated with PCA used SVD

(c) 多尺度梯度求和的PCA提取的梯度方向(c) Orientation maps estimated with PCA used method of multi-scale gradients summation圖4 提取圖像局部方向的對比Fig.4 Comparison of the image local orientation extraction methods

對于給定圖像I中的一個像素點(x,y)，在N×N的鄰域內，其垂直與水平方向的梯度可以構成一個大小為N2×2的局部梯度矩陣G：

(1)

對G矩陣進行SVD分解：

G=USVT

(2)

2.2.2 PCA計算的多尺度梯度求和法

(3)

式中，∑表示對i×i鄰域進行求和，尺度i決定了鄰域范圍大小。 最終的主方向計算是基于多個尺度平均梯度的加權融合計算得到的Gxx、Gyy、Gxy：

(4)

式中，權值wi的定義如下：

(5)

最后像素點的主方向φ計算方法如下:

(6)

(7)

圖4(c)給出了使用多尺度梯度求和計算的主方向，可見多尺度梯度求和計算的主方向與SVD計算的主方向(圖4(b))基本一致，都對噪聲干擾有良好適應性，但是梯度求和法的計算量要比SVD少，并且可以通過積分圖像優化進一步減少計算量。

3 匹配算法的實現與加速方法

3.1 基于多尺度PCA-HOG描述的模板匹配算法

本文提出的匹配算法是一種全搜索的模板匹配算法，算法計算模板與基準圖上候選窗口之間的相關性，并逐像素地移動候選窗口以尋找最佳匹配目標窗口。圖5所示為全搜索模板匹配的示意。這種全搜索方法是假設圖像已經進行了幾何校正，模板和基準圖之間只存在平移變換，這樣的全搜索匹配算法常常在視覺導航和制導領域使用。

圖5 全搜索模板匹配示意Fig.5 Diagram of full-search template matching

設A={a(x1,y1),a(x2,y2),…,a(xn,yn)}為模板窗口，大小為n×n；基準圖為B={b(α1,β1),b(α2,β2),…,b(αm,βm)}，大小為m×m。則所提出的模板匹配方法可以定義為如式(8)所示的最優化問題：

(8)

式中：P代表匹配輸出的位置結果；Ii為基準圖B上候選窗口，Uw為基準圖B上候選窗口的集合；T(·)表示將灰度圖像轉換為PCA-HOG特征向量；D(·)為兩特征向量之間的相似性度量函數，本文使用NCC[9]作為相似性度量函數。

圖6 算法流程Fig.6 Flowchart of the proposed algorithm

算法的整體流程如圖6所示，關鍵步驟包括：

1)使用Sobel算子提取x與y方向的圖像梯度Gx與Gy。

2)利用PCA計算每個像素的局部主方向。

3)利用主方向與梯度幅值計算PCA-HOG特征。

4)利用快速傅里葉變換(Fast Fourier Transform，FFT)加速的NCC計算模板與基準圖中候選窗口之間的相似性系數。

5)輸出相似性系數最大的候選窗口位置作為匹配結果。

3.2 算法加速

計算效率優化可以使用金字塔策略或者基于梯度下降的搜索優化策略，但是這樣的優化方法與全搜索方法不等價，可能錯失最優結果。本文對算法計算效率進行了兩方面的改進，能夠在與全搜索方法等價的條件下(確保不遺漏最優結果)大幅提升計算速度。

1)使用積分圖像對基于梯度求和的PCA計算進行加速。

2)FFT對NCC相似性系數計算進行加速。

3.2.1 積分圖像加速PCA

積分圖像中每個像素的值是從圖像原點(左上角)到這一點所構成矩形內所有像素值的和[24]。因此，圖7中，原圖像矩形區域內的像素值之和sum可以按照下式求得：

sum=s3-s2-s4+s1

(9)

式中，s1、s2、s3、s4是矩形區域4個角位置對應在積分圖像上的像素值。因此無論矩形區域大小，使用積分圖像求和都只需要三次加法計算。對于大小為M×N的圖像，從原圖像得到積分圖的計算復雜性為O(M×N)。所以如果需要對圖像上很多矩形區域求和，并且矩形區域較大時，積分圖像求和方法可以節省大量計算。

梯度求和法的PCA計算的主要計算量為每個像素的鄰域梯度求和。本文利用積分圖像對式(3)中的局部梯度進行求和。假設鄰域矩形大小為5像素×5像素，則直接求和的計算量約為25×M×N，而使用積分圖像求和的計算量約為4×M×N，計算速度能提高5倍多。

(a) 原圖像(a) Original image (b) 積分圖像(b) Integral image圖7 基于積分圖像求和Fig.7 Summation based on integral image

3.2.2 FFT加速NCC相似性計算

卷積定理可以利用FFT對卷積過程進行加速[22]。假設卷積操作中的基準圖大小為m像素×m像素, 模板大小n像素×n像素，若搜索時確保模板包含在基準圖內，則在空間域直接計算卷積的復雜度為O((m-n+1)2n2)，而使用FFT加速后的卷積計算復雜度為O(2m2log2m)。比較可以發現，當模板較大時，使用FFT加速后的卷積可以大幅減少計算量。NCC的計算公式如下：

(10)

4 實驗結果及分析

4.1 數據集

使用20對遙感異源圖像對提出的算法和現有的模板匹配算法進行比較測試。測試圖像包括SAR、紅外、激光雷達和可見光圖像。圖像中的地物類型主要包括：城市建筑、道路、橋梁、河流和山脈等。圖8給出了一些測試匹配對的例子，可見這些遙感異源圖像之間存在嚴重的灰度畸變，并且SAR圖像還存在嚴重的噪聲畸變和細節缺失，這些都給圖像匹配帶來了很大的困難。

(a) 可見光圖像與紅外圖像對比(a) Comparison between visible image and infrared image

(b) 可見光圖像與激光雷達圖像對比(b) Comparison between visible image and laser radar image

(c) 可見光圖像與SAR圖像對比(c) Comparison between visible image and SAR image圖8 測試集圖例Fig.8 Samples from testing data

每一對異源圖像都包括實時圖像和基準圖像，從實時圖像中選取模板匹配到基準圖上。選取的模板分為4種不同的大小，分別為32像素×32像素、64像素×64像素、96像素×96像素、128像素×128像素。對于各尺度的模板，在實時圖內隨機選取25個模板。基準圖像的尺寸都是320像素×320像素。此外，對測試圖像進行歸一化處理后添加了三種不同方差(0.01,0.03,0.05)的高斯噪聲來測試算法的抗噪聲能力。圖9給出了原圖像和添加噪聲圖像的例子，其中圖9(a)為原圖像，圖9(b)～(d)添加了不同方差的高斯噪聲。

圖9 噪聲示例Fig.9 Noise samples

4.2 評價標準

使用匹配正確率(Correct Matching Rate，CMR)作為算法表現的評價標準。CMR=CM/R，其中R是總匹配對次數，CM是匹配正確的次數。在本文，當匹配位置與實際位置重疊面積比(Overlapping Area Ratio, OAR)達到90%以上時，匹配結果被認定為正確的。

(11)

式中：TW為模板的邊長；Δx、Δy為匹配位置和實際位置的誤差；f(x)為截斷函數，

(12)

4.3 算法比較結果與分析

本節通過與現有算法MI[16]、基于方向梯度通道特征(Channel Features of Orientated Gradients，CFOG)的模板匹配算法[12]、HOG-NCC、HOPC[21]的實驗結果相比較，來分析所提算法PCA-HOG在遙感異源圖像匹配中的優越性。下面簡要介紹以上算法：

1)CFOG：是一種基于圖像方向梯度的像素化特征表示，是對HOG特征描述符的擴展，基于CFOG的模板匹配具有優越的計算效率，但是因為CFOG直接基于梯度構建，存在對噪聲適應性不足的問題。

2)HOG-NCC：該匹配算法使用HOG作為特征描述，并使用NCC作為相似度度量方法。相比于本文算法，HOG-NCC并沒有使用PCA提取主方向，也沒有使用積分圖像和FFT加速方法。

3)HOPC：該算法先使用log-Gabor提取方向，然后按方向統計相位一致性，得到相位一致性直方圖特征，最后使用NCC作為相似性度量進行匹配。HOPC利用了圖像的相位信息，但是相位對噪聲干擾也比較敏感。

圖10顯示了對原圖像(沒有添加噪聲)，各種算法使用不同大小模板的匹配正確率。結果表明本文提出的匹配算法PCA-HOG整體高于其他幾種匹配算法，平均正確率高出第二名(CFOG)3.5%。另外可以看出，基于結構特征的匹配算法(CFOG、HOPC、PCA-HOG、HOG-NCC)的正確率都相比于直接基于圖像灰度的匹配算法(MI)高。這一結果表明，對于非線性灰度畸變非常明顯的遙感異源圖像，直接對圖像灰度進行相似性度量的匹配算法不是非常有效，應該使用特征描述先提取結構信息，然后再匹配。

圖10 原圖像匹配結果Fig.10 Matching result on original images

此外，對測試算法在不同程度噪聲干擾下的結果進行了比較，結果如圖11所示?？梢钥闯觯兴惴ǖ腃MR都隨著噪聲水平的增加而降低，而PCA-HOG相較于其他幾種算法，整體匹配正確率更高。另外，隨著噪聲強度的增加，相比于基于圖像灰度的算法(MI)，基于結構特征的算法(CFOG、HOPC、PCA-HOG、HOG-NCC)正確率下降得都更快。尤其是當高斯噪聲方差為0.05時，一些基于結構特征的匹配算法正確率要小于MI。這是因為這些算法的結構特征提取方法都是基于梯度或相位，梯度或相位對噪聲干擾非常敏感。PCA-HOG算法的匹配正確率明顯優于HOG-NCC算法，后者沒有使用PCA提取主方向，這驗證了使用PCA對HOG特征增強可以改善匹配性能。

圖12給出了部分PCA-HOG匹配結果圖，其中的SAR圖像帶有嚴重的噪聲，且和可見光圖像、紅外圖像間有著很大的灰度差異，使用PCA-HOG算法仍然可以正確匹配圖像，這說明該算法能有效匹配帶有灰度畸變和噪聲干擾的遙感異源圖像。

(a) 模板尺寸： 32像素×32像素(a) Template size: 32 pixel×32 pixel

(b) 模板尺寸： 64像素×64像素(b) Template size: 64 pixel×64 pixel

(c) 模板尺寸： 96像素×96像素(c) Template size: 96 pixel×96 pixel

(d) 模板尺寸： 128像素×128像素(d) Template size: 128 pixel×128 pixel圖11 添加噪聲情況下測試算法的匹配結果Fig.11 Matching results of tested algorithms on images with noise

(a) SAR-可見光(將SAR圖像疊加到可見光圖上)(a) SAR-visible (bottom layer is visible image; top layer is SAR image)

(b) SAR-紅外(將SAR圖像疊加到紅外圖像上)(b) SAR-infrared (bottom layer is infrared image; top layer is SAR image)圖12 PCA-HOG算法的匹配結果示例Fig.12 Samples of PCA-HOG matching results

對提出的算法加速方法也進行了測試。測試平臺為Intel CPU i7 單核2.10 GHz，對使用加速方法和不使用加速方法的PCA-HOG算法的運行時間進行了測試，結果如表1所示，其中基準圖尺寸為320像素×320像素。

表1 算法運行時間比較

從表1可以看出，在使用加速方法的情況下，算法的運行速度有1～2個數量級的提升，這與本文3.2節對算法復雜性的分析基本一致。

5 結論

針對遙感異源圖像匹配存在的問題，提出了一種基于PCA增強HOG特征描述的快速遙感異源圖像匹配算法。該算法利用PCA對局部主方向進行提取，改進了HOG描述對噪聲的適應性。實驗對算法在遙感異源圖像數據集上進行了評估，并且與現有的算法進行了比較。實驗結果表明，本文提出的匹配算法在正確率上明顯優于其他測試算法。本文算法還利用積分圖像對特征提取過程進行加速，并使用FFT對匹配搜索過程進行了加速。實驗部分，對比了加速前后的匹配算法，測試結果表明加速方法能將算法的計算效率提高1～2個數量級。本文提出的匹配算法和其他4種比較的算法都是按照整像素平移搜索實現的模板匹配算法，因此，不能實現亞像素匹配定位。后續將對基于本算法的亞像素精確匹配開展研究。