一種基于相位一致性相關的多源遙感影像配準方法

范登科，潘 勵，葉沅鑫

(武漢大學遙感信息工程學院，武漢 430079)

1 概述

圖像配準作為一項基礎性的影像幾何處理技術，在醫學、計算機視覺等學科領域中發揮著重要作用。盡管它在遙感中的研究和應用起步較晚，但是由于遙感影像在采集環境、觀測角、地形起伏等多方面因素作用下所體現出的復雜性遠遠大于其他圖像，再加上其本身海量及多源異構的特性，對配準技術在該領域的研究和應用提出新的挑戰。目前，國內外眾多學者針對遙感影像匹配問題的不同方面，提出了很多具有可行性和研究價值的配準策略和方法，也因此產生了多種將配準方法進行分類的方式，包括基于配準影像數據源類型的分類、基于影像間幾何畸變類型的分類等。最常見的是基于匹配時所使用影像信息類型的分類方式［1]，該方式將配準分為特征匹配和區域匹配兩類。

特征匹配是當前配準技術研究的熱點方向，它通過計算影像的梯度、矩等強度或對象的幾何形態信息來描述局部范圍內的顯著性特征，然后選擇適用于該特征匹配的相似性測度建立同名點的對應關系，最終完成對幾何畸變的改正。常見的用于遙感影像配準的點特征有:Harris、Forsterner、Morovac，線形或邊緣特征有:Canny、Sobel，面狀或區域特征有Hu不變矩、形狀上下文。SIFT［2]是近幾年圖像配準研究領域興起的一種具有尺度和旋轉不變性描述的特征算子，它以其多尺度、穩定的點特征表達特性贏得許多研究學者的青睞，并由此衍生出許多基于SIFT的配準方法，如PCA-SIFT［3]、CSIFT［4]、SR-SIFT［5]等。然而，由于大多數特征算子都是以圖像梯度的強度和方向為基礎的，且很少顧及區域性信息，造成多源異構遙感影像配準的同名點匹配率和糾正精度不能滿足實際應用生產的需要，因此特征匹配方法具有應用局限性和問題針對性。

區域匹配作為一種早期發展成熟的、方便易用的圖像配準方法，被廣泛應用于各類商業化圖像處理軟件中。它通過統計運算的方式對特定區域內的圖像信息進行描述，并借助相關關系模型度量彼此間的相似程度，相對于特征匹配需要穩定且充足的特征作為配準前提條件來說，區域匹配具有很強的通用性。目前研究和應用最多的區域匹配方法有相關系數和互信息，它們分別是基于影像原始灰度和信息熵進行計算和匹配的。Gao使用相關系數作為匹配策略［6]，借助影像金字塔和二次多項式實現對Landsat GeoCover數據子像素精度的自動配準;Peter［7]同樣以相關系數為核心匹配手段，完成對HyMap、AIRSAR等24組同源和非同源影像對的幾何校正;Sahil［8]通過迭代計算TerraSAR-X和Ikonos影像間的互信息，研究分析了城鎮區域異構圖像配準的精度。由于非同源異構影像間在原始灰度上通常是非線性變化的，且信息熵的計算量比較大，導致區域匹配方法的推廣應用受到一定的制約。

本文綜合特征匹配和區域匹配的優勢，提出一種適用于非同源異構遙感影像的配準方法，該方法以不同濾波尺度下位置穩定不變的角點作為控制點，借助金字塔分層映射搜索策略，通過計算領域和搜索區域內的相位一致性，引入相關系數作為相似性測度實現同名點的匹配。并最終給出利用可見光、近紅外及SAR影像對該方法進行配準實驗的驗證和分析。

2 配準方法

2.1 控制點提取

為了保證提取足夠多且位置相對穩定的控制點進行匹配和模型參數估計，需要使用通用性強、可靠性高的角點特征提取算子在待配準的輸入影像上執行計算和檢測。當前，Harris［9]和 FAST［10]被證明具有較高的角點提取重復率，適用于各類光譜組成結構的影像。本文選擇Harris算子在輸入影像上提取特征控制點，Harris通過局部圖像信息的自相關方式，在各個方向都計算了灰度變化，其角點提取的原理是:首先使用標準偏離差為σ的高斯濾波模板與原始灰度影像I進行卷積運算，得到平滑后的影像ˉI，然后計算ˉI的一階梯度矩陣M，公式為

其中k為介于0.04～0.06的常數，如果2個特征值都比較大，說明在該點處沿任意方向運動會引起判別函數R的急劇變化，判斷為具有強響應的角點，否則為弱邊緣或平滑域點，在使用閾值T對該值進行局部非最大抑制后，可以獲得穩定的特征點。

由于使用不同標準偏差σ的高斯卷積模板會使平滑處理后的影像ˉI變模糊，在特定濾波器尺度下檢測的角點可能會在其他尺度的濾波影像中發生位置偏移甚至消失。為了提取不受高斯濾波尺度影響的穩定角點，需要使用臨近尺度的高斯模板重新提取角點進行檢測和判斷。假設角點(x0，y0)是使用尺度(標準偏差值)為σ0的高斯濾波模板所提取的Harris特征點，構造尺度為 σ1=σ0-δ×σ0和 σ2=σ0+δ×σ0的高斯濾波模板對原始影像I進行卷積運算，提取的角點分別為(x1，y1)和(x2，y2)，當提取的角點滿足條件

并且

時，認為(x0，y0)是一個不受卷積尺度影響的穩定角點，選為控制點，否則被排除。公式(3)、(4)中D代表歐式距離，即以角點(x1，y1)和(x2，y2)相對于(x0，y0)的位置偏移作為特征點的判定依據。圖1展示了分別使用原始Harris算子和多尺度高斯濾波后Harris提取的特征點分布情況，可以看出，在加入了公式(3)、(4)的判別后，位于道路邊緣和開闊地中心的點受平滑尺度變化的影響，位置變化顯著，被判斷為不穩定的特征點排除，保留的點為穩定的不受平滑濾波干擾的特征點。

2.2 同名點匹配

圖1 Harris算子提取特征點分布示意

由于Harris特征算子提取的控制點為梯度變化顯著的角點，因此在其領域內影像的紋理信息比較豐富，灰度變化顯著，相關系數在小范圍內容易形成變化劇烈的極值，可以選擇作為同名點的匹配測度。然而，非同源異構遙感影像間的光譜組成存在明顯差異，在灰度強度上呈現非線性變化，導致相關系數不能穩定描述相似度高的同名區域。針對這一問題，本文引入相位一致性對影像進行轉換，得到相關性程度高的特征圖像。

相位一致性是一種使用傅里葉諧波分量描述信號局部強度的特征，它利用特定方向θ不同尺度s的Log Gabor濾波器，將信號分解為頻率域下與θ和s相對應的傅里葉諧波分量，通過加權疊合不同尺度的分量，獲得具有原始信號強弱特性描述的響應。Morrone證明了相位一致性在描述圖像特征信息上與人類神經視覺具有相一致的敏感度［11]。Wong指出，相位一致性特征具有抵抗圖像對比度和亮度變化的能力，并利用其構造的最小和最大距特征，分別實現了遙感影像特征點的提取和同名點的匹配［12]。應用于圖像二維離散信號，相位一致性的計算公式為

式中，PC(x，y)為點(x，y)處的相位一致性，θ為濾波方向，在區間［0，π]之間以固定步長變化，n為濾波器尺度s的最大值，Wθ(x，y)為頻率擴展的權重系數，As，θ(x，y)代表使用沿 θ方向尺度為 s的 Log Gabor濾波器進行分解后諧波分量的振幅值，」表示僅取正值的運算，負值全部賦為0，T為噪聲閾值，ε是一個避免除零的常數，Δφs，θ(x，y)可由分解后諧波分量的相位角計算獲得，公式為

圖2 Landsat 5 TM影像及其相位一致性特征圖像

從圖2可以看出，邊緣和角點處的相位一致性明顯大于平坦開闊地區，通常在角點處取得極值。同時，相位一致性特征影像反映的并不是二值化的邊緣，而代表的是不同尺度和方向頻率域卷積運算下傅里葉諧波分量的累積值。因此，以特征點為中心的窗口發生位置改變時，不會引起相關系數的劇烈變化，且由于角點處相位一致性的顯著特征，使得極值比較容易獲得。經計算對應圖像的相關系數，圖2中(a)與(c)之間的相關系數為-0.045，而(b)與(d)達到0.605，可見采用相位一致性特征度量影像間的相關性明顯優于原始灰度。

遙感影像在采集過程中一般表現出大視角、寬視場的特性，成像的圖幅范圍大，使全局匹配的效率降低，計算復雜度增大，需要使用特定的搜索策略優化同名點匹配過程。基于金字塔結構的同名區域預測與分層索引映射是目前最受歡迎的、應用于遙感影像匹配中的搜索策略，方法是以一定比例逐級降低圖像空間分辨率，縮小圖幅范圍。首先，在輸入影像的頂層金字塔圖像上提取特征點，并搜索參考影像對應尺度的金字塔，找出相關系數大的窗口中心作為匹配點，使用簡單的二維平面轉換關系模擬初始的幾何變形，根據最小二乘準則刪除殘差大的誤匹配點對;而后，利用不同尺度金字塔圖像的映射關系，在下一層更高分辨率尺度的金字塔圖像上獲取同名區域，執行與上述一致的精確匹配，直到達到底層即原始分辨率影像為止。兩層金字塔結構的同名區域映射關系及幾何變換模型估計如圖3所示。

2.3 模型參數估計與重采樣

利用匹配獲取的同名點對，依據最小二乘原理估計幾何畸變模型的參數，剔除殘差大的同名點對以保證模型估計精度在可接受的范圍內，最終使用特定的插值方法實現對輸入影像的重采樣處理，輸出在地理位置上與參照影像配準的校正影像。在這一過程中，全局性幾何畸變模型的選擇顯得十分重要。

圖3 基于金字塔結構的遙感影像同名區域索引和搜索策略

由于傳感器姿態和側視角的差異，采集獲得的多源異構遙感影像在全局幾何畸變上表現出不同形式。中分辨率光學影像的成像方式為推掃式，其傳感器所處軌道高度及線陣掃描方式決定了影像幾乎不受側視角和地形起伏的影響，在全局影像上僅存在平移、旋轉和尺度畸變，因此可以使用仿射模型描述待配準光學影像間的幾何變形關系。仿射變換模型公式為

式中，(x，y)和(X，Y)分別為對應同名點在輸入影像和參考影像的圖像坐標，λ代表尺度縮放因子，θ表示影像的相對旋轉角，(c，r)是圖像在二維平面的相對平移量。因此僅需要兩對同名點即可以求解該模型參數。而同分辨率SAR影像的采集使用主動傳感方式，雖然軌道高度避免了地形起伏的影響，但是其成像圖幅范圍較小，并且存在一定的側視角，因此對于光學影像與SAR影像以及不同時相SAR影像間的幾何畸變，通常采用投影變形模型來描述，一次項投影模型公式為

式中包括a～h總計8個模型參數，需要4對同名點進行求解，分母部分描述了成像平面的傾斜變形。在最小二乘原理基礎上對模型參數進行估計后，使用參考影像中的控制點坐標反算輸入影像控制點的模型坐標，并與真實坐標求差后計算均方根誤差RMSE，RMSE值的大小代表了該同名點對偏離模型估計的程度，刪除RMSE值最大的多個點對，重新估計模型參數，迭代執行這一過程直到同名點對集合中最大的RMSE值小于預先設定的閾值，以獲得滿足糾正精度要求的模型參數。最終，雙線性插值方法被用于重采樣過程輸出配準結果。

2.4 總體配準策略

本文提出的基于相位一致性相關的多源遙感影像配準方法可概括為以下7個步驟:(1)根據影像大小，設定特定的分辨率比例，分別構建輸入和參考影像的多級金字塔圖像;(2)在輸入影像的頂層金字塔上使用多尺度Harris提取控制點，并計算其領域內窗口的相位一致性;(3)計算參照影像頂層金字塔的相位一致性特征圖像，使用相關系數作為匹配測度全局搜索最大相關的窗口，以其中心作為匹配同名點;(4)用簡單二維平面變換模型估計畸變改正關系，預測下一級金字塔圖像的同名區域;(5)在下級金字塔圖像的子區域內迭代執行步驟(2)～(4)，直到完成原始分辨率影像上同名點的匹配;(6)根據待配準影像的類型確定全局幾何變換模型，采用最小二乘原理迭代解算和修正模型精度;(7)重采樣輸出配準后的結果影像。圖4展示了配準策略的完整流程。

圖4 配準方法流程

3 實驗與分析

3.1 實驗背景和數據

實驗以Windows XP系統為平臺，在Visual C++2008集成開發環境中設計開發程序實現上述配準方法，程序運行的硬件環境為Inter Core2 CPU 2.66 GHz，2 GB內存。實驗數據選用3組遙感影像，每一組中輸入影像和參考影像分別采集自異構光譜范圍或不同傳感器。前2組數據均由光學遙感影像組成，采集的區域為湖北省大東湖和梁子湖水域，其中第1組數據來自于同一傳感器，為了體現空間位置上的差異性，手工對輸入影像旋轉10°，第2組數據來自不同傳感器，存在微小旋轉及10～20像素的平移變形。第3組數據為光學影像與SAR影像之間的配準，為了正確計算相關系數，對SAR影像進行重采樣，確保其與光學影像的分辨率一致。3組實驗數據的詳細信息如表1所示。

在進行影像配準運算處理過程中，需要對相關參數進行初始化，具體項目及其初始值為:金字塔各級比例1∶3，金字塔級數3，原始分辨率下搜索區域大小200像素，相關系數閾值0.75，相關系數計算窗口大小13像素，RMSE閾值1.5。

3.2 配準精度分析

使用以上3組數據對本文配準方法進行實驗驗證，并對該方法的適用性和精度進行評價。表2列出了經實驗獲得的配準評價結果。

由表2可以看出，在配準精度方面，第3組數據的效果最好，最大殘差和總體均方根誤差指標都控制在1個像素以內，然而同名點對的誤匹配率卻達到了95%以上，僅有5對正確匹配的同名點用于模型參數估計，且它們分布集中于圖5(e)、(f)中左下角的島嶼海岸線上，占影像中主體的陸地范圍內并沒有正確匹配的同名點對，這是因為大陸中地物類型豐富而海島內單一，SAR影像之于前者表現出更多的噪聲，而相位一致性的計算結果對噪聲十分敏感，造成陸地中相位一致性的相關程度并不是十分顯著;反之，海島在SAR影像中紋理結構簡單，海岸線輪廓特征明顯，有利于相位一致性穩定描述局部圖像信息。

表1 實驗數據信息概覽

表2 配準方法的實驗評價結果

圖5 3組實驗數據影像及同名點分布情況

前2組實驗數據的配準結果在精度和誤匹配率上達到一致水平，盡管第2組數據的總體均方根中誤差略微偏高，但其匹配得到的同名點對數目是第1組的3倍，是因為第2組數據的影像對均由近紅外波段構成，相對于第1組的可見光-近紅外組合，相關程度更高。然而，在匹配獲得大量同名點對的同時，發生誤匹配的點對數目也相應提升，且從2組數據對比情況來看，這一比例水平并不受光譜組成結構的影響，因此要進一步降低誤匹配率，仍需要對匹配方法和策略進行改進。另一方面，由于第2組影像采集自不同時間，地物變化也會造成同名點的誤匹配。從圖5(a)～(d)中可以看出，正確匹配的同名點均勻分布于整個圖幅范圍內，且在對比度高的局部區域內相對集中，符合相位一致性描述亮度和對比度差異的不變特征。

綜合分析3組實驗數據的配準結果，本文方法更適用于對非同源異構的兩景光學影像進行配準，對于光學與SAR影像的配準有一定局限性和區域選擇性，同時將相關系數作為匹配測度引入到該方法中，對于全局仿射和投影變形有一定的抵抗能力。

3.3 對比驗證

為了進一步驗證本文方法的配準效果，引入傳統基于灰度的相關系數(NCC)、尺度不變特征(SIFT)、改進尺度不變特征(SR-SIFT)3種方法分別對3組實驗數據進行配準實驗，統計各項實驗結果的誤匹配率和均方根中誤差(RMSE)，得到的對比情況如表3所示。

表3 3種常用方法實驗結果對比(“—”表示配準失敗)

實驗過程中設定RMSE閾值為1.0像素，以保證幾何校正模型參數估計的準確性和一致性。當3種方法達到相同精度水平時，結果顯示實驗2的誤匹配率最低，原因是兩景近紅外影像間灰度相關性程度最高，且與特征表述相關的梯度強度和方向變化較為穩定，與此同時，NCC和SR-SIFT方法的誤匹配率明顯低于本文方法。相反，3種方法對實驗1數據進行配準的誤匹配率高于或接近本文方法，可見，對于光譜組成差異越顯著的影像來說，本文方法的適用性更強，效果更好。3種方法不能完成對實驗3數據配準的事實更進一步驗證了這一結論。因此，與3種配準方法進行對比分析得出的結論是，基于相位一致性相關的配準方法可以應對傳統方法所不能解決的異構遙感影像間難以配準的問題。

4 結語

本文提出一種基于相位一致性相關的多源遙感影像配準方法，該方法以多級圖像金字塔結構為索引和搜索策略，首先使用多尺度高斯濾波判別和遴選穩定的Harris特征點，在控制點領域和相關系數極值搜索區域內，引入相位一致性特征進行相似性度量，在降低誤匹配概率的同時獲得更多精確匹配的同名點對。經實驗分析和對比驗證，該方法適用于中分辨率多源異構光學-光學、光學-SAR遙感影像間的配準，精度滿足實際生產過程中影像幾何校正處理的需要，且高于其他基于區域匹配的配準方法。下一步的研究工作將專注于提高運算效率，改進算法以降低誤匹配率，并適用于不同分辨率、地形起伏等復雜條件下的多源遙感影像配準。同時希望本文方法能為圖像配準領域的學者提供研究基礎和借鑒。

［1]Zitová B，Flusser J.Image registration methods:A survey［J].Image and Vision Computing，2003，21(11):977-1000.

［2]Lowe D G.Distinctive image features from scale-invariant keypoints［J].International Journal of Computer Vision，2004，60(2):91-110.

［3]Ke Y，Sukthankar R.PCA-SIFT:A more distinctive representation for local image descriptors［C].Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition，2004:506-513.

［4]Abdel-Hakim A E，Farag A A.CSIFT:A SIFT descriptor with color invariant characteristics［C].Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition，2006:1978-1983.

［5]Yi Z，Zhiguo C，Yang X.Multi-spectral remote image registration based on SIFT［J].Electronics Letters，2008，44(2):107-108.

［6]Gao F，Masek J，Wolfe R E.Automated registration and orthorectification package for Landsat and Landsat-like data processing［J].Journal of Applied Remote Sensing，2009，3(1).

［7]Bunting P，Labrosse F，Lucas R.A multi-resolution area-based technique for automatic multi-modal image registration［J].Image and Vision Computing，2010，28(8):1203-1219.

［8]Suri S，Reinartz P.Mutual-information-based registration of Terra-SAR-X and Ikonos imagery in urban area［J].IEEE Transactions on Geoscience and Remote Sensing，2010，48(2):939-949.

［9]Harris C，Stephens M K.A combined corner and edge detector［C].Proc.Alvey Vision Conf.，1988:147-152.

［10]Rosten E，Drummond T.Machine learning for high-speed corner detection［C].Proceedings of the European Conference on Computer Vision，2006:430-443.

［11]Morrone M C，Owens R A.Feature detection from local energy［J].Pattern Recognition Letters，1987，6(5):303-313.

［12]Wong A，Clausi D A.AISIR:Automated inter-sensor/inter-band satellite image registration using robust complex wavelet feature representations［J].Pattern Recognition Letters，2010，31(10):1160-1167.