程玉芳，李景山

(1.中國科學院空天信息創(chuàng)新研究院，北京 100094；2.中國科學院大學，北京 100094)

0 引言

多源遙感圖像匹配是指將不同傳感器、不同視角、不同時相獲取的同一地域的兩幅或多幅圖像通過匹配算法實現(xiàn)自動獲取圖像間同名點的過程[1]。

常用的圖像匹配方法主要有基于特征的方法和基于區(qū)域的方法?；谔卣鞯钠ヅ渌惴ㄆ浠舅枷胧翘崛D像中顯著不變的幾何特征，利用描述子描述特征，然后建立相似性度量準則，根據特征之間的相似性并結合其他的約束條件進行匹配和篩選。其中，尺度不變特征變換[2](scale invariant feature transform，SIFT)的描述子具有旋轉不變、亮度魯棒、對噪聲和視角的微小變化穩(wěn)定等優(yōu)點，在遙感圖像匹配領域得到了廣泛應用。基于區(qū)域的圖像配準方法的主要思想是根據已知的模板圖像，按照特定的相似性準則計算模板圖像與搜索圖像區(qū)域子圖之間的相關程度，對比所有子圖的相似度確定最佳的子圖位置?；趨^(qū)域灰度的配準方法是早期圖像配準中最常用的方法，常用的相似性度量方法有互相關[3]、互信息[4]等方法。

傳統(tǒng)的利用圖像灰度梯度變化來描述幾何結構的方法很容易受到圖像的輻射照度和視角變化的影響。相較而言，在描述圖像的結構信息方面相位一致性與人類神經視覺具有基本一致的敏感性[5]，對圖像中的輻射亮度變化、對比度變化和噪聲干擾等均不敏感。Ye等[6]采用相位一致性模型的強度和方向信息構建了相位一致性方向直方圖(histogram of orientated phase congruency，HOPC)，利用該特征描述輸入圖像的幾何結構信息，通過幾何結構的相似性進行匹配，能有效應用于光學與SAR影像之間的配準。王新生等[7]通過相位一致性模型構造出最大索引圖，在最大索引圖上使用分布直方圖構造特征描述符來完成圖像匹配，取得了較好的實用效果。

在現(xiàn)有研究成果基礎上，本文針對多源光學遙感圖像間存在的非線性輻射畸變、幾何定位差異大等問題，提出了一種利用圖像相位一致性特征進行同名點預測及匹配的方法。首先，通過遙感圖像自帶的地理位置信息將圖像劃分為規(guī)則格網，通過圖像的紋理結構進行相位相關得到初始的相對偏移量，以此預測同名點的空間搜索范圍；然后，以特征點為中心，構建點鄰域內的方向相位一致性結構特征描述符，采用傅里葉變換方法進行模板窗口的特征匹配；最后，通過最小二乘法剔除誤差較大的點，得到最終的同名點對，實現(xiàn)圖像間的自動匹配。

1 多源遙感圖像匹配算法及實現(xiàn)

本文多源遙感圖像匹配的總體算法設計如下。

1)圖像分塊相位一致性特征提取。將輸入圖像劃分為W×W規(guī)則大小的格網塊，通過地理位置映射對應參考影像區(qū)域并裁剪同樣大小圖像，利用相位一致性模型計算二者邊緣圖像以及待配準圖像格網塊中的角點圖像用于角點檢測。

2)邊緣圖像相位相關。將分割好的待配準圖像塊與對應基準影像格網塊的邊緣圖像進行相位相關，將相關脈沖尖峰位置記為兩張格網圖像間的相對偏移量。

3)誤相關剔除。利用RANSAC隨機采樣一致性模型作為幾何約束條件，篩選出正確的相關峰，即有效的格網偏移量集合。

4)特征點提取?；?)提取的角點圖像檢測輸入圖像中顯著的特征角點，并利用非極大值抑制在輸入圖像上提取出均勻分布的特征點。

5)方向相位一致性特征提取。以特征點為中心，統(tǒng)計特征點鄰域內的離散方向相位一致性特征，作為該點的特征模板。

6)幾何約束及匹配。利用3)得到的格網偏移量，預測特征點在基準影像上的搜索位置，提取方向相位一致性結構特征，并將其變換到頻率空間，提高模板匹配的計算效率，快速定位同名點。

7)誤匹配剔除?；谧钚《斯烙嬈ヅ潼c殘差，迭代剔除，得到最終的同名點集。

匹配算法總體設計見圖1。

圖1 自動匹配算法實現(xiàn)流程

1.1 相位一致性特征提取

相位一致性(phase congruency，PC)[8]模型是通過計算局部能量函數來描述圖像的結構信息，是一個無量綱量，與信號的整體大小無關，這個特性保證了相位一致性特征對于光照和對比度具有很好的不變性[9]。

Log-Gabor小波可以真實自然地描述圖像的頻率響應[10]，它能夠利用不同方向、不同中心頻率的濾波器來獲取圖像的局部紋理特征，因此本文選擇Log-Gabor提取圖像的相位一致性特征。指定方向θ、不同尺度s的Log-Gabor濾波器計算相位一致性PC(s，θ)，參照Jain[11]的方法，濾波器方向間隔選擇30°，即將θ分為6個方向(0°，30°，60°，90°，120°，150°)，濾波器尺度s的最大值設置為4，共24個Log-Gabor濾波器對輸入圖像進行濾波，提取圖像的相位一致性特征。

將待匹配遙感影像劃分為W×W大小互不重疊的格網塊，并通過格網左上角點的地理位置映射對應參考影像點圖像坐標，裁剪同樣尺寸大小圖幅，計算兩張小圖的相位一致性特征。劃分成格網塊的目的一方面是適用于本文的匹配算法，另一方面也是為了加速圖像的傅里葉變換過程。

1.2 邊緣圖像相位相關

圖像信號的傅里葉變換包含幅度與相位兩個部分，根據傅里葉變換相關理論[12]可知，函數位移不會改變傅里葉變換的模(幅值)，但是會改變實部與虛部之間的能量分布，其結果是產生一個與位移量成正比的相移。

因此通過計算兩頻譜之間的互功率譜，并進行傅里葉逆變換，得到的相位相關脈沖圖像除了在它們相對偏移(u0，v0)處不為零外其余位置都基本接近零值，這樣就很容易得到兩幅圖像之間的相對平移參數。

由于多種外部因素的影響，異源遙感影像上同一地物往往存在很大的非線性輻射差異。盡管在一定程度上相位相關算法能夠容忍兩幅圖像間存在的非均勻輻照差異，但是隨著圖像間光照差異的逐漸增大，相位相關的顯著度也會下降[13]。

PC圖的最大矩M代表圖像的邊緣特征，而邊緣結構具有較好的抗輻射畸變能力，地物的幾何結構信息在異源影像上也基本保持不變，因此可以通過影像中的幾何結構進行相關，初步獲取地物間幾何位置的偏移量。即可以通過計算兩幅圖像最大矩的互功率譜得到相位相關圖像，此相位相關脈沖圖像的脈沖峰位置即為兩幅影像之間的偏移量(圖2)(注：圖2(e)中在(25，443)處有一個尖銳的峰值，兩幅圖像大小為500像素×500像素，左上起點坐標為(1，1)，它們的相對偏移量為(24，-58))。

圖2 相位相關算法示例

為確保分塊相位相關方法的有效性，還需要為算法結果設置閾值進行篩選。

1)相關強度閾值Tmag=0.03。相位相關圖中最大響應峰值大小表征圖像相關性的強弱程度，若響應峰值小于給定的強度閾值Tmag，則表明它們之間相關程度較弱，相關失敗。

2)相關比例閾值Tratio=0.75。若相位相關脈沖圖像中次響應峰值與最大響應峰值之比大于Tratio，則認為兩幅圖像間的相關性不夠顯著，判定為相關失敗。

對W×W大小的格網塊提取出的邊緣結構特征計算相位相關圖像，記錄脈沖尖峰位置和其平面坐標，即為兩圖像之間的相對平移量。

1.3 誤相關剔除

由于遙感圖像的復雜性，通過幾何結構信息進行相位相關也可能存在誤相關情況，如大面積水域、荒漠等特征不明顯區(qū)域。因此仍需對獲取的初始脈沖集合進一步篩選，得到準確有效的偏移量集合。

RANSAC(random sample consensus)算法[14]通過迭代的方式從一組包含離群的被觀測數據中估計出一個認為符合此樣本的數學模型，從而將符合數學模型的點標記為內點(正常數據)，偏差大的點標記為外點(異常數據)。

分塊相位相關得到的初始脈沖峰集合實際上表征的是輸入圖像各格網塊與參考圖像各格網塊之間的相對偏移量，而整幅圖像各格網的偏移量應該基本符合某個變換模型，所以本文采用RANSAC算法提純分塊相位相關結果，避免對之后的匹配過程提供錯誤的幾何約束。

1.4 特征點提取

大多數特征檢測都是通過尋找局部極值作為特征點，這會導致特征點分布受圖像對比度影響在高對比度區(qū)域分布密集、低對比度區(qū)域分布稀疏的情況。

PC圖的最小矩m為圖像的角點特征，通過閾值T提取特征點，即：對于任意位置(x，y)，如果m(x，y)>T，則該位置將被標識為候選特征點。但直接通過設定閾值提取特征點同樣會造成關鍵點在圖像上非均勻分布。因此本文通過非極大值抑制方法，鄰域極值的查詢半徑為r，抑制區(qū)域大小為(2r+1)×(2r+1)，剔除容易受到噪聲干擾的低對比度的點，篩選出圖像中高對比度的點。然后，在每d×d的格網范圍內，取局部極大值中的最大值點，記錄其位置信息作為最終的特征點，使檢測到的關鍵點能夠均勻分布在整幅遙感圖像上。

圖3為SIFT算法與非極大值抑制角點檢測方法的對比圖，表明非極大值抑制角點檢測能夠從圖像中提取均勻分布的關鍵點。

圖3 特征檢測算法對比圖(影像大?。? 000像素×1 000像素)

1.5 方向相位一致性特征

相位一致性的方向代表著圖像結構特征沿此方向劇烈變化，類似于梯度方向，可采用Log-Gabor小波的奇對稱濾波器來計算，其卷積結果表示影像在某個方向的能量變化，各個方向上的相位一致性可以反映不同方向的邊緣強度信息。因此，本文利用離散方向集的相位一致性描述圖像中每個點的特征。

構造不同方向的Log-Gabor濾波器與輸入圖像進行卷積得到各方向上的相位一致性PC(θ)，θ共6個方向(0°，30°，60°，90°，120°，150°)，即可以得到一組離散方向集的相位一致性特征圖像(PCθi，i=1，2，3，4，5，6)。

圖4 方向相位一致性特征提取

在構造區(qū)域的特征時，為使其更加穩(wěn)定，需要充分考慮特征點鄰域的結構特征。因此，本文提出一種基于方向相位一致性特征的匹配模板，構造以關鍵點為中心的區(qū)域方向相位一致性特征，確定模板圖像半徑和搜索半徑大小，再對不同角度的邊緣強度進行組合，形成關鍵點鄰域的方向相位一致性特征描述子，使其能夠抵抗非均勻光照變化所引起的圖像間非線性輻射差異。

1.6 幾何約束及匹配

考慮傳統(tǒng)模板匹配方法逐點滑動計算，搜索范圍大，計算量大，本文為了提高同名點匹配效率，先將待匹配圖像與參考圖像劃分為規(guī)則格網塊，提取邊緣圖像進行相位相關，獲取各格網塊之間的偏移參數。在特征點匹配的過程中，從已得到的偏移量中取距該特征點最近的格網相對偏移量，將其認為是該特征點與其同名點之間的幾何誤差，從而預測特征點對應的同名點位置，縮小搜索范圍。

在左上角點經緯度及圖像分辨率完全相同的情況下，特征點坐標(x，y)與預測同名點位(x′，y′)的對應圖像坐標關系為式(1)。

(1)

式中：(dx，dy)表示點(x，y)的偏移量，即距(x，y)最近的有效相關峰位置。

由于特征點與其預測的同名點的區(qū)域模板的方向相位一致性特征是逐像素的結構特征描述符，數據量大，計算耗時，因此本文采用傅里葉變換將特征模板從空間域轉換至頻率域，加速模板匹配過程，提高匹配性能。

誤差平方和算法(sum of squared differences，SSD)是一種比較經典判斷模板圖像與搜索圖像之間的相似度度量算法，計算方法見式(2)。

SSD(x，y)=∑i∑j[S(i+x，j+y)-T(i，j)]2

(2)

式中：T為匹配的模板圖像；S是搜索圖像。上式展開，可得式(3)。

SSD=S*S-2S*T+T*T

(3)

式中：*為卷積運算符。第一項S*S和第三項T*T為固定值，因此可以繼續(xù)將相似性度量Sim簡化為式(4)。

Sim=S*T=FS×FT

(4)

式中：FS和FT分別為S、T傅里葉變換結果。當模板圖像大小為N×N，搜索半徑為M/2時，搜索圖像大小為(M+N)×(M+N)，空間域計算SSD的時間復雜度為O(M2N2)，而將其計算過程簡化到頻率域，時間消耗為O((M+N)2log(M+N)，計算效率可以得到顯著提升。

圖5為基于方向相位一致性特征的模板匹配示例，在待匹配圖像中檢測到的關鍵點A通過幾何約束條件預測其在參考圖像中的同名點B’位置，分別以A、B’為中心取(M+N)×(M+N)鄰域計算該區(qū)域的方向相位一致性特征，再利用傅里葉加速的誤差平方和算法即可求得兩模板特征間相似性，相似性最大值處即為關鍵點A實際對應的同名點位置B。

注：模板半徑為50×50，搜索半徑10×10，此時搜索圖像大小為120×120。左上起點坐標為(1，1)，關鍵點A(61，61)的實際對應的同名點位置為B(65，59)。圖5 基于方向相位一致性特征的模板匹配示例

1.7 誤匹配剔除

盡管基于相位一致性特征的匹配流程已經具有較高的正確匹配概率，但在實際匹配過程中發(fā)現(xiàn)仍然會有部分偏差，即通過相似性度量測度之后，仍然不可避免地會有少量的錯誤匹配。為了剔除這些錯誤匹配，需要再次使用最小二乘法通過最小化誤差對得到的初始匹配點集進行提純。在通過最小二乘解算出模型參數后，計算每個點的殘差，并剔除殘差較大的同名點，直至殘差小于閾值時退出迭代過程，得到最終的匹配點。

2 實驗分析

為了驗證本文方法在多源光學遙感影像匹配中的有效性，以中國遙感衛(wèi)星地面站通用遙感衛(wèi)星數據預處理系統(tǒng)生產的國產高分系列衛(wèi)星遙感影像作為研究對象。參考影像數據為91衛(wèi)圖下載的地圖數據，具有多種分辨率且影像地物清晰、定位精度高，可作為系統(tǒng)幾何校正產品精度檢驗的參考基準。3組實驗數據的待配準圖像分別來自于高分一號寬幅相機、高分二號多光譜相機、高分七號后視多光譜相機，實驗所用的輸入影像數據均已經過系統(tǒng)幾何校正，投影坐標系為UTM坐標系，參考影像也經重采樣到同一分辨率下，消除了影像間的旋轉差異和較大的幾何差異。實驗數據基本信息見表1。

表1 實驗數據基本信息

2.1 實驗參數

本文特征點提取設置的檢測范圍為250像素，即每250×250 網格內提取1個候選點。分塊格網大小設置為1 000像素×1 000像素，則每塊區(qū)域最多能提取16個特征點。模板匹配的模板圖像設為50像素×50像素，搜索半徑為10，即搜索圖像大小為120像素×120像素。借助于GDAL、FFTW等源碼庫，采用C++編程，設計并行程序(核心數為14)實現(xiàn)本文的匹配流程。為驗證本文多源光學遙感圖像匹配方法的有效性，實驗將所提出的匹配方法與傳統(tǒng)SIFT算法進行對比分析，SIFT實驗環(huán)境參數與本文方法一致。

實驗機群刀片服務器配置：Intel Xeon Gold 6 132(2.6 GHz/14 c)/2 666 MHz/10.4 GT×2，內存為16 GB DDR4 2 666 ECC REG×8。

2.2 實驗結果分析

3組實驗數據同時存在傳感器成像物理特性、時相和輻射差異。實驗一圖像數據包含耕地和山脈，實驗二的地物覆蓋類型為耕地和水域，實驗三影像中心為城市，包含大量耕地，且有部分山脈，地物信息較為豐富。3組實驗數據的匹配結果見圖6。表2為實驗結果對比分析。

表2顯示，第一組實驗得到的匹配點數量基本一致，但第二組、第三組實驗數據由于圖像間時間間隔長，存在較大的輻射差異，SIFT匹配得到的同名點數量很少，而本文方法并未受傳感器物理特性、成像角度、時相等因素的影響，仍然可以得到較多的同名點。從圖6可知，本文的匹配方法在三組實驗數據中得到的匹配點都能夠基本均勻分布在整幅遙感圖像上。

將參考影像作為輸入產品的幾何精度檢驗基準，可以看出本文方法和SIFT方法的幾何精度檢驗結果基本接近，差距在亞像素以內，并且RMSE最大為1.03像素，說明本文的匹配方法能夠達到較高的匹配精度。此外，三組實驗檢查到的最大幾何誤差為28.62像素，表明本文方法有較大的位移檢測范圍。

從匹配效率來講，3組實驗中本文方法耗時均大于SIFT，這主要是在相位相關及特征匹配過程中都需要計算圖像4個尺度、6個方向上的結構信息。但通過并行程序設計，二者仍處于同一數量級，差距并不明顯，可以滿足系統(tǒng)幾何校正產品的實時匹配需求。

圖6 實驗數據匹配結果圖

表2 實驗結果對比分析

3 結束語

對于多源遙感圖像匹配來講，相比于灰度信息，其幾何結構和形狀等屬性更加穩(wěn)定，即這些特性對傳感器成像差異、光照、時相等變化不敏感。一方面，本文利用圖像的邊緣結構進行分塊相位相關，不受圖像噪聲、照度變化的干擾，預測同名點位置，使得本文方法在較小的搜索半徑內仍然具有較大的位移檢測范圍。另一方面，離散方向集的相位一致性特征可以反映不同方向上的邊緣結構特征，這種特征幾乎不依賴于圖像灰度，因而具有很強的抗輻射畸變能力，從而能夠實現(xiàn)同名點的精確匹配。

實驗結果表明，本文方法可以有效應對多源遙感圖像的非均勻輻射畸變，并克服圖像間較大的幾何誤差，從而獲得均勻分布的同名點，實現(xiàn)遙感圖像間的自動匹配。