采用各向異性尺度空間的遙感圖像配準*

李 宏，王 鵬，畢 波，唐錦萍

(1.東北石油大學 a.電氣信息工程學院;b.數學與統計學院，黑龍江 大慶 163318；2.海南醫學院 公共衛生學院，海口 571101；3.黑龍江大學 數據科學與技術學院，哈爾濱 150080)

0 引 言

圖像配準是在不同的時間、不同的成像設備或不同的條件(天氣、光亮、相機位置和視角等)下對采集的兩個或多個圖像進行匹配的過程[1]，它是許多遙感任務(例如變化檢測、圖像拼接和環境監控)必不可少的部分。不同于基于頻域的方法，空間域中的基于特征的方法[2]主要依賴于提取顯著特征。顯著特征是圖像域的子集，通常是獨立的點、連續曲線或連續區域[3-4]。

在基于特征的方法中，最經典的算法是SIFT(Scale Invariant Feature Transform)，而基于SIFT的圖像配準方法主要使用高斯濾波構建尺度空間。因為高斯函數是各向同性，因此待配準圖像的邊緣和細節同噪聲一樣變得模糊，關鍵點很難實現亞像素級別的準確定位，所以特征描述子的區別度低。文獻[5]提出了KAZE算法，首次使用各向異性技術生成非線性尺度空間，檢測和描述性能優于傳統的SIFT，但是使用Hessian算子需要計算每層尺度空間的二階梯度，對噪聲特別敏感。文獻[6]提出了加速KAZE(Accelerated-KAZE,AKAZE)算法，采用快速顯式擴散函數和二進制特征描述子，在保留KAZE優良性能的同時加快了KAZE算法的計算時間。文獻[7]提出了基于多特征點的圖像配準算法，聯合各向異性尺度空間和Harris角點響應的策略，提高了魯棒性和亞像素級別的準確定位。文獻[8]提出了一種結合SURF(Speeded-Up Robust Features)和AKAZE的改進算法，使用了Hessian算子和二進制特征描述子進行特征描述。文獻[9]提出了基于改進AKAZE的配準算法，速度比AKAZE提高了50%且特征點也大幅增加，但是更適用于喀斯特地區的遙感圖像。文獻[10]提出了基于特征的圖像配準的算法，利用改進的復合非線性擴散濾波策略和逐漸降低的多尺度金字塔空間分辨率及二進制代碼串匹配提高了準確性和計算效率，但在計算和匹配二進制代碼串花費時間較長。

對于油田地區，有效掌握地面建筑物的實時位置及狀態信息對于油田生產的安全尤為重要。然而，油田遙感圖像中的建筑物的形狀不均勻，所以圖像配準更加困難。

本文提出了一種基于各向異性尺度空間的遙感圖像配準方法，使用了自適應的側窗濾波來構建各向異性尺度空間；接著選取了一種采用了改進的鄰域分塊思想進行特征描述；最后采用了一種增強的匹配算法進行匹配點過濾，提高了計算速度，得到了更多精確的匹配點。

1 基于多特征點的遙感圖像配準算法

基于多特征點的遙感圖像配準算法[7]主要由各向異性尺度空間生成、尺度空間特征點檢測、特征點描述符提取和特征點匹配[11]四個步驟構成，其流程圖如圖1所示。

圖1 基于多特征點的遙感圖像配準算法流程圖

1.1 梯度計算的原理

為了增加描述符對圖像亮度的顯著非線性差異的魯棒性，基于多特征點的遙感圖像配準方法為高斯尺度空間中的每個像素給了梯度定義。首先，原算法通過Sobel濾波器計算高斯尺度空間圖像的梯度幅度為

(1)

(2)

(3)

1.2 各向異性尺度空間的生成

各向異性尺度空間將圖像的亮度擴散到非線性比例空間中，可以由各向異性擴散方程來描述。Weickert等人[12]引進了加性算子分裂(Additive Operator Splitting,AOS)算法進行求解各向異性擴散方程。

(4)

公式(4)就是AOS算法的根本思想，由于圖像是二維的，所以m=2。

各向異性技術的算法偽代碼如下：

Input：原始輸入圖像I，時間度量尺度值tn，n∈[0,1,…,S]，圖像的高度M和寬度N

Output：各向異性尺度空間圖像，其中Ln,n∈[0,1,…,S]

1對圖像I進行標準差為σ0的高斯濾波，并把濾波結果作為尺度空間第0層圖像L0

2按照下面過程迭代生成各向異性尺度空間：

3 forn=0,1,…,S-1

(3)計算擴散時間步長τ，其中τ=tn+1-tn

fori=0,…,M-1

end for

forj=0,1,…,N-1

end for

(6)由下面公式計算各向異性尺度空間第n+1層圖像Ln：

4 end for

5 輸出各向異性尺度空間圖像Ln,n∈[0,1,…,S]

1.3 局限性

由于大慶油田分布錯綜復雜，而油田遙感圖像中的建筑物形狀不均勻，因此，從任一種遙感影像數據中很難將復雜的油田圖像準確、完整地提取出來。另外，遙感圖像中地物復雜，干擾因素較多，而且基于高分辨率遙感影像進行圖像配準技術環節，每道工序對檢測的最后結果都有著重要的影響。基于多特征點的遙感圖像配準方法應用在無人機拍攝的油田遙感圖像上可以得到很好的精度，但是得到的特征點匹配對少，而且花費時間也較長。

2 基于各向異性尺度空間的遙感圖像配準

本文對原算法的改進主要分以下幾個部分：

(1)基于各向異性技術，采用自適應的側窗濾波思想，更好地保留了圖像邊緣，得到了更加穩定的特征點；

(2)基于極坐標的平面系統，采用環形和扇形結合的分塊思想，有效地降低了原算法中特征描述子的維度，縮短了在特征描述和特征匹配階段所花費的時間；

(3)基于RANSAC算法，提出了增強的匹配算法，增加了特征點匹配的數量，減少了運行時間。

2.1 自適應的側窗濾波

基于窗口的操作是圖像處理中最常用的操作之一，傳統做法幾乎總是將窗口的中心與要處理的像素對齊，這種傳統觀點并非普遍適用。Yin等人[13]提出了新的側窗濾波方法(Side Window Filter,SWF)，可以將窗口的側面或角落與要處理的像素對齊。在圖像平滑、去噪、增強和保留結構的紋理去除等應用中，實施SWF方法可以顯著提高其邊緣保留能力并達到最好的性能。

SWF方法首先在自變量θ和r連續的情況下定義了側窗，如圖2(a)所示，θ為旋轉的窗口與水平線之間的角度，r是窗口的寬度，p∈{0,r},(x,y)是目標像素i的位置。然后，通過改變θ并固定(x,y)，可以改變窗口的方向，同時將其側面與像素i對齊。

圖2 側窗的定義

(5)

SWF方法保留邊緣通過最小化目標函數來實現，其目標函數為

(6)

基于SWF方法實現了許多流行的線性和非線性濾波算法，所以可以顯著提高其邊緣保留能力。考慮到SWF方法應用在非線性濾波算法的顯著性和靈活性，改進的方法是，側窗的半徑采用每層的尺度計算的擴散時間步長，而非固定值。

SIFT方法中尺度空間有O組，每組S層，尺度空間以金字塔的形式存在。為了簡單起見，改進算法未對圖像執行下采樣操作，即各向異性尺度空間中圖像的分辨率是一樣的。設尺度空間共S+1層，則各層圖像的尺度為

σn=σ02n/3,n∈[0,1,…,S]。

(7)

式中：σ0為尺度空間中第一層圖像的尺度。從公式(7)可以看出，和原SIFT算法相同，相鄰層圖像的尺度比是21/3。因為非線性擴散方程以時間進行的定義，因此把尺度值轉換為時間值，轉換公式為

(8)

各向異性尺度空間的構建中采用了自適應的窗口半徑，可以顯著提高其邊緣保留能力并達到最好的性能。側窗的半徑即時間步長，轉換公式為

τ=tn+1-tn。

(9)

本文采用自適應的SWF方法構建了各向異性擴散的尺度空間，相對于原算法對尺度空間中第n層圖像Ln做標準差σ恒為1的高斯濾波處理，更好地保留了圖像邊緣。針對特征檢測中得到的特征點多且重復度高的問題，改進的方法是在各向異性尺度空間提取到特征點后使用非極大值抑制。

2.2 改進的描述符

馬曉娟[14]總結了對三維模型采用網格劃分的特征提取算法，可以實現對三維模型進行不同形式的合理劃分。本文將三維模型網格劃分區域的思想應用到二維鄰域劃分上，如圖3所示。

圖3 二維鄰域劃分圖

改進的方法是獲取關鍵點為圓心、半徑為R1=12σ的圓形鄰域，然后把取得的鄰域劃分為3個區域和9個子區域的對數極坐標網格生成特征描述子。特征點的圓形鄰域的區域劃分如圖4所示，其中R2/R1=0.73，R3/R1=0.25,此時圓形鄰域的9個子區域有近似一致的面積。每個小型區域形成一個8維的梯度幅度和方向直方圖，因此最后改進的描述子是72維。

圖4 新的描述符的鄰域圖

改進算法的特征描述子的維度由原來的136維下降到了現在的72維，明顯降低了描述子的維度。由于遙感圖像一般有很多局部一致性，所以相較整幅圖像，區域劃分的思想對視角變換有較好的優勢：增加的鄰域信息使得特征描述符更穩定，為后期圖像匹配的準確性提供了保證。

2.3 增強的特征匹配

在特征匹配過程中，諸如規則化最小二乘之類的核方法基函數的數量等于訓練樣本的數量N。因此，匹配過程在時間和空間上復雜度分別為O(N3)和O(N2)，當要匹配的特征點較多時，增加了計算量。周榮榮等人[15]提出了SAR配準算法，把向量場共識(Vector Field Consensus,VFC)算法[16]在模擬和真實的山地SAR圖像的關鍵點配準，獲得了很高的匹配正確率。

稀疏向量場共識(Sparse VFC)算法[17]由期望最大化算法(Expectation Maximization,EM)可以得到完全對數后驗Q(θ,θold)，即

(10)

Sparse VFC算法偽代碼如下：

Input：訓練集S={(xn,yn),x∈RP,y∈RD,n∈NN}，N為訓練集的數量，核矩陣Γ，正則化常數λ，基函數數量M，潛在變量zn

Output：向量場f，內在點對集合T

1初始化a,β,γ,V=ODN×1,P=IM×N,σ2,ε

Γ(xi,xj)=e-β‖xi,xj‖2(ωID×D+(1-ωD)ID×D),

3 Repeat

(1)EM中E-step

通過貝葉斯定理計算后驗概率分布pn=P(zn=1|xn,yn,θold)，計算對角矩陣P=diag(p1,p2,…,pN)，pn表示樣本n與當前向量場f相一致的程度:

(2)M-step

①通過下面的公式求解線性系統更新C，計算向量場f:

②使用vn=(f(x1)T,…,f(xN)T)T更新V:

③使用稀疏近似最大化對數似然函數估計θnew=argmaxθQ(θ,θold)更新σ2和γ。對Q(θ)分別求解σ2和γσ2的偏導，并令其結果為0，可求得σ2和γ，tr(·)表示矩陣的跡：

Q(θ)=log(p(Y|X,θ)p(f))

4 Until目標函數Q(θ)收斂，收斂準則為Q(θ)<ε

5 向量場f由②中的公式決定，內在點對集合是T={n:pn>τ,n∈RN}，其中τ是預定義的閾值

提出的增強的特征點匹配步驟如下：

Step1 首次匹配：特征點通過對應描述符的最小和次小的反余弦之比進行初始匹配，首先獲得初始匹配點對集P0。

Step2 再次匹配：Sparse VFC算法用于從點對集P0中刪除錯誤匹配點，然后得到特征點對集P1。

Step3 三次匹配：RANSAC算法用于從點對集P1中刪除錯誤匹配點并計算初始變換參數μ，最后獲得特征點對集合P2。

3 實驗與分析

3.1 實驗數據及預處理

本文進行兩組實驗。

遙感圖像對A為多光譜遙感圖像對P-A[18-19]；參考圖像是從Landsat-7 ETM+衛星于2000年7月24日拍攝的多光譜遙感圖像中提取的第5頻段圖像，分辨率是614 pixel×611 pixel；待配準圖像是從Landsat 4-5TM衛星在1996年6月28日拍攝的多光譜圖像中提取的第三頻段圖像，分辨率為614 pixel×611 pixel。為了增加圖像的復雜度，把待配準圖像順時針旋轉90°。由于兩幅圖像來自不同的傳感器和不同的頻段，因此圖像中相同區域的像素灰度存在顯著的非線性差異。

遙感圖像對B使用無人機拍攝獲取，無人機采用索尼DSC-RX1RM2相機(焦距為35 mm)，拍攝的單張圖像為7 952 pixel×5 304 pixel。遙感區域為大慶油田喇嘛甸地區，并對采集得到的照片進行處理，最后得到了管線正射原始數據的圖像。由于無人機拍攝的油田遙感圖像分辨率太大，使得程序無法直接處理，所以，我們先進行雙三次插值把分辨率縮小為900 pixel×600 pixel。B組中參考圖像和待配準圖像的建筑物形狀不均勻且灰度存在顯著差異。

3.2 評估標準

主觀評價通過視覺對匹配參數精度進行評估，容易受到人為因素的影響；而客觀評價標準通過數字對匹配參數進行評估，因此會更具有說服力，通常選用正確匹配的點對個數、正確率和均方根誤差(Root Mean Square Error,RMSE)對參數精度進行衡量。

3.2.1 正確匹配的點對個數和正確率

正確匹配的點對個數越多，匹配算法計算得到的相似變換的參數就越準確，因此可以把正確匹配點對個數作為算法性能評價標準[20]。正確匹配個數和算法初始檢測到的特征點個數有關，而一般情況下增加初始檢測到的點個數會增加最后正確匹配點對個數，但是正確率會相應下降，因此為了評價標準的有效性，需要對比算法和本文提出的算法正確匹配點對個數和正確率一起比較。

3.2.2 均方根誤差

(11)

3.3 實驗過程及結果

為了評估算法的有效性，我們進行了A組和B組實驗，并將改進算法與原算法及其他主流算法(如SAR-SIFT[19]、KAZE[5])進行了比較。SAR-SIFT算法是針對SAR圖像配準設計的，而KAZE算法首次使用各向異性技術進行圖像配準，是比較經典的算法。這些算法在Matlab R2018b下使用Intel Core i3-3220 @ 3.30 GHz 雙核處理器、8 GB 1600 MHz物理內存和AMD Radeon HD 6700 Series(1 GB/藍寶石)顯卡來實現。

按照圖1的算法流程圖開始實驗，步驟如下：

Step1 采用SWF方法生成8層各向異性尺度空間，其中擴散時間步長由公式(9)計算得到。

Step2 在Step 1生成的尺度空間中計算Harris響應，若實際計算的角點響應值R大于常量閾值0.8且大于其附近8鄰域內其他點的角點響應，則該像素點為特征點，把該點保存；接著在各向異性尺度空間同一層使用非極大值抑制，非極大值抑制的區域半徑設置為3。Harris算子和角點響應函數分別見公式(12)和(13):

(12)

R(X,σn)=det(u(X,σn))-D·tr(u(X,σn))2。

(13)

Step3 以Step 2中得到的特征點為鄰域中心，采用新的描述符對Step 1中的尺度空間進行極坐標劃分，分別對9個區域統計梯度幅度和方向直方圖，生成新的72維描述符。

Step4 對參考圖像和待配準圖像經過Step 1～3后得到的兩組特征點進行增強的特征匹配。首次匹配中反余弦比的閾值設置為0.9；再次匹配中Sparse VFC算法參數設置為,a=10，γ=0.9，λ=3，τ=0.75，ε=1×10-5,M=16；三次匹配中，RANSAC算法亞像素精度的誤差閾值設置為2。

為加強不同算法的對比，本文調整了對比算法的首次匹配中歐幾里得距離之比的閾值和改進算法的反余弦之比的閾值，保證了不同的對比算法提取的參考圖像/待匹配圖像特征點數量大致相同。兩組實驗的對比算法的特征點匹配結果以連線的形式顯示，A組實驗結果見圖5，B組實驗結果見圖6。

圖5 A組不同算法的特征匹配效果圖

圖6 B組不同算法的特征匹配效果圖

為進一步定量描述算法，本文統計了不同對比算法的特征點數量、匹配對數量、正確匹配對數量、正確匹配率、RMSE以及算法的運行總時間的數據，A組和B組統計結果見表1。為保證算法的效果，所有算法執行10次，以10次時間的平均值作為最終結果。

表1 兩組不同算法實驗結果

3.4 實驗結果分析

觀察圖5可以看出，在A組實驗中，SAR-SIFT算法效果很差，其他算法得到的匹配點大都分布在圖像的凹凸明顯和黑暗明顯的交界處。通過圖6可以看出，在B組配準中，各種算法得到的匹配點大都分布在建筑物角點附近和顏色突變較大的地方；前三種算法得到的匹配對都分布在待配準圖像的右側，而改進后的算法得到的匹配對左右兩側分布均勻，而且可以清楚看到改進后的算法的匹配線也比前三種算法的更加密集，特別是在大廠房旁邊的角點進行了較好的匹配。因此，結果證明了改進算法可以明顯的改善匹配效果。

從表1中的各項統計指標分析可知，在A組配準中，SAR-SIFT算法配準失敗，改進算法和原算法匹配準確度和精度都較高，明顯高于KAZE算法，但是改進算法與原算法相比，不但增加了正確匹配對的數量，而且明顯減少了運行時間;在B組配準中，與其他三種算法相比，改進算法的匹配正確率和配準精度都最高，而且與原算法相比，不但大大增加了正確匹配對的數量，而且提高了正確匹配率，明顯減少了運行時間。

由以上分析可知，A組和B組的配準結果都驗證了改進算法的有效性。

4 結束語

本文針對基于多特征點的遙感圖像配準算法應用在無人機遙感圖像上存在正確匹配點數量相對較少特別是運行時間長的問題，提出了改進算法。該算法首先使用自適應的側窗濾波技術構建了各向異性尺度空間，顯著地保留了圖像的邊緣；然后使用一種新型的描述符，減少了特征描述子的維度；最后使用一種更加魯棒性的匹配點過濾算法，得到了更多精確的匹配點同時也保證了配準精度。實驗結果表明，該算法在正確匹配對的數量、匹配的精度和時間方面顯示出比原算法和其他主流算法更好的性能。

增強的匹配點過濾雖然明顯改善了原算法，但正確率仍有待提高，這也是今后研究的重點。