基于特征點的SAR影像配準方法

劉智勇 來立永 張耿斌 祁宏昌潘屹峰 彭林才 吳希文

1)廣東電網有限責任公司， 廣州供電局， 廣州 510620 2)中科云圖科技有限公司， 鄭州 450003 3)廣東工業大學， 土木與交通工程學院， 廣州 510006 4)廣東工業大學， 大灣區城市環境安全與綠色發展教育部重點實驗室， 廣州 510006

0 引言

合成孔徑雷達干涉測量技術(Interferometric Synthetic Aperture Radar， InSAR)已被廣泛應用于自然災害的監測工作中， 如地震(Wangetal.， 2007， 2009， 2019)、 地表沉降(王華等， 2014a； Ngetal.， 2018)、 滑坡(Tarchietal.， 2003； Sunetal.， 2015)等。通常， 在形成干涉圖之前需對SAR影像進行配準并重新采樣到與參考影像(主影像)統一的網格上。配準作為InSAR數據處理的關鍵步驟， 其精度直接影響InSAR生成干涉圖的質量。為了不影響生成干涉條紋的質量， 對于Stripmap/ScanSAR成像模式的影像配準誤差必須要控制在0.125個像素以下(王超等， 2002)， 對于方位向漸進式掃描的地形觀測(Terrain Observation with Progressive Scans in azimuth， TOPS)成像模式， 其配準精度需要達到千分之一個像素才不至于影響干涉處理(吳文豪等， 2017)。

SAR影像配準主要用于匹配多時相、 多傳感器或從不同視點拍攝的2個或多個影像(Brown， 1992)。常規的配準方法是在主強度圖上按照一定的間距規則選取控制點， 即規則點(Regular Points， RP)。由于規則點的選取不考慮地面的散射特性， 在具有較高相干性的區域效果較好， 但在實際應用中， 規則點經常位于低相干、 失相干區域(例如大面積的水域、 植被茂密的森林或季節變化的農業區域等)， 導致配準存在計算量大、 誤匹配率高以及可靠性低等缺陷。

為了改善常規配準方法所存在的缺陷， 一些研究人員利用永久散射體(Permanent Scatterer， PS)、 孤立點(Isolated Points， IPs)等特殊點目標替代規則點進行配準。雖然PS點具有很好的穩定性， 但PS點的可靠識別需要大量SAR影像， 而研究中通常只有很少的數據可用， 致使此方法受到較大的限制。此外， IPs通常是孤立的比較明亮的點(Serafino， 2006)， IPs能有效避開失相干的區域(例如水體)， 但由于SAR影像的成像特性， 這些亮點可能會聚集在一定的空間中， 因此可能會導致類似斑狀的偏移估計。此外， SAR影像中存在大量的斑點噪聲， 呈現出高亮度值， 同樣會導致孤立點定位錯誤。

配準的主要目的是獲取主、 輔影像之間的坐標轉換參數， 除了PS點和IPs點之外， 可利用具有明顯特征的點進行配準， 因此， 特征提取是SAR影像配準技術中的關鍵， 準確的特征提取為配準的成功進行提供了保障。Bay等(2006)提出的加速魯棒特征算子(Speeded Up Robust Feature， SURF)算法是用于遙感影像特征檢測與描述的流行方法之一， 其可利用SURF算子對SAR影像進行配準(熊智等， 2011； Durgametal.， 2016)。但由于SAR影像在成像時存在大量噪聲點， 僅通過使用SURF算法進行配準的精度低且誤匹配率高， 難以滿足InSAR的精度要求。

針對上述配準問題， 本文提出了一種結合SURF、 土地覆蓋數據和互相關算法等的配準方法。首先， 對主強度圖進行分塊處理， 并將每個子塊數據的灰度值線性拉伸到0～255； 然后， 通過SURF對每個子塊進行特征點探測， 接著利用土地覆蓋數據剔除不當的特征點， 對保留的特征點進行互相關； 最后， 通過多項式擬合獲取主、 輔影像之間的仿射變換參數。本文通過珠三角地區的SAR影像數據的配準實驗評估該方法的性能。

1 Sentinel-1A雷達數據及配準方法

1.1 數據

本文所用的研究數據覆蓋了廣州、 佛山、 東莞等珠三角地區， 采用2015年12月12日以及2015年12月24日的一對C波段Sentinel-1A雷達影像數據， 其成像模式為TOPS模式， 影像大小為20278像素×68514像素， 在距離向和方位向上的像素大小為2.3m×14.1m， 空間分辨率為5m×23m。

1.2 基于特征點的SAR影像配準

標準的配準方法首先在主強度圖上按照等間距選取控制點， 然后采用相關系數法進行像素級配準， 接著進行過采樣處理， 并再進行一次互相關以實現亞像素級配準， 最終進行多項式擬合計算主、 輔影像之間的轉換參數(靳國旺， 2007)。本文提出的SAR影像配準方法首先利用SURF算子探測主影像的特征點(Feature Points， FP)， 然后進行水體掩膜， 剔除無效特征點， 并對保留的特征點進行互相關， 最后進行多項式擬合獲取主、 輔影像之間的轉換參數。

1.2.1 利用SURF提取特征點

SURF算法是由Bay等(2006)提出的尺度與旋轉不變的興趣點檢測和描述符， 是改進尺度不變特征變換算子(Scaleinvariant feature transform， SIFT)(Low， 2004)的一種特征檢測算法， 在獨特性和魯棒性方面均近似甚至優于SIFT， 且計算速度快于SIFT。其特征點提取主要流程是：

(1)構造Hessian矩陣， 并通過Hessian矩陣的行列式初步判定特征點位置

當行列式值為0時， 可判定該點比鄰域內的其他點更加明亮或者暗淡。為提高計算性能， SURF使用盒式濾波器模板代替高斯二階微分模板。圖 1 為9×9的盒式濾波器， 左側2張圖分別是y方向和xy的二階高斯微分模板， 右側2張圖是采用盒子濾波器對其近似的結果， 模板中的數字(-1、 1、 -2)表示對應顏色區域的權值， 灰色區域權值為0。因為盒式濾波器在積分影像上的計算速度非常快， 故SURF與SIFT相比在性能上有了很大的提升。

圖 1 盒式濾波器模板示意圖Fig. 1 Schematic diagram of box filter template.

(2)建立影像金字塔

SIFT通過保持濾波器模板不變， 對原始影像進行降采樣獲得影像金字塔(圖 2a)。SURF通過改變濾波器的尺寸而影像不變來保持尺度不變性(圖 2b)。與SIFT相比， SURF具有更快的速度和更高的精度。

圖 2 影像的金字塔Fig. 2 The image pyramid.

(3)定位特征點

經Hessian矩陣處理后的像素點須與相同或相鄰尺度空間中的26個像素進行比較(圖 3)， 濾除響應值較弱及錯誤定位的像素點， 并選定穩定像素點作為特征點。

圖 3 特征點定位Fig. 3 Location of feature points.

SAR影像的成像過程是由不同散射體的回波路徑相干疊加而成， 導致SAR影像中存在大量或明或暗的斑點， 即相干斑(Eom， 2011)， 這降低了SAR影像的目標識別能力。為了抑制SAR影像的相干斑， 可進行多視處理， 但該處理同時會降低SAR影像的分辨率。本文中， 為在確保不犧牲SAR影像分辨率的前提下提高特征點識別的速度并抑制噪聲， 首先對每個子塊的灰度值進行初步處理， 以保證整體影像不會太明亮或太過暗淡。在0.5%≤灰度值 ≤99.5% 的范圍內將其進行線性拉伸到0～255， 對于灰度值 >99.5% 或 <0.5% 的像素則分別設置其值為255和0。對主強度圖進行初步處理后， 使用SURF算子探測特征點。SURF可處理大塊數據， 具體窗口大小可根據計算機的物理內存進行配置。本文為加快計算速度并防止整個檢測窗口位于水體中， 統一將特征點的探測窗口大小設置為5000像素×5000像素。

從圖 4 中可以比較明顯地看出水體與陸地之間的灰度值差異， 這對陸地特征點的準確定位有所幫助。設置Hessian的閾值為25000， Octaves為4， Octaves layers為3， 一共探測到了8693個特征點， 大部分都集中于陸地(圖4a中的紅色圓圈)， 但有一部分特征點仍落入水體中(圖4a中的藍色圓圈)， 其主要原因是在珠江流域及入海口停泊了大量船只。圖4a黃色方框區域的放大圖如圖4b、 c所示， 從圖4b中可以明顯看出這些船只在雷達強度圖中相對于周圍水體具有非常明顯的差異， 導致SURF將其識別為特征點(圖4b右側)。但經過一段時間后在相同區域再次進行檢測， 探測到的特征點位置隨船只的移動具有很大的差異(圖4c右側)， 而陸地上探測到的特征點仍然具有非常穩定的特性， 其相應的位置基本保持不變。因此， 在進行配準之前， 需利用外部土地覆蓋數據剔除落入水中的特征點， 以剔除不穩定的特征點。

圖 4 特征點分布圖Fig. 4 Distribution map of feature points.a 主強度圖的特征點； b 主影像局部特征點； c輔影像局部特征點

1.2.2 水體掩膜

在水域中， 干涉圖的相干性均非常低(接近0)(王華等， 2014b)。同樣， 在大面積水域中進行互相關， 相關度極低且偏移量極不穩定， 這不僅大大降低了配準的計算效率， 同時還導致錯誤的偏移量估計， 使得最終計算的主、 輔影像轉換模型不可靠， 故必須對水體進行掩膜處理。盡管SURF可有效避開水體， 但當Hessian閾值太低、 探測窗口太小或整個探測窗口都在水中時， SURF探測到的特征點可能位于水中。因此， 針對在主強度圖中探測到的特征點需要進行水體掩膜處理， 以剔除位于水體中的特征點。

本文利用SWBD(Shuttle Radar Topography Mission Water Body Dataset)數據濾除SAR強度圖中不適當的特征點， SWBD數據可從相關網站(1)http: ∥geodata.policy support.org/swbd。下載。SWBD數據是僅有-1和0的二值二進制文件， -1為水體， 0為陸地。由于SWBD是地理坐標， 而強度圖是雷達坐標， 故需要將SWBD數據轉換到雷達坐標系下， 判斷特征點所在位置的SWBD值是否為0， 否則舍棄。

1.2.3 互相關計算

經水體掩膜后， 進行互相關計算。互相關可衡量2幅影像的相關程度， 計算公式為

(1)

式中，M、S表示主、 輔影像的強度，E為求平均。當|ρ|=0時， 表示主、 輔影像完全不相關； 當|ρ|=1時， 表示主、 輔影像完全相關。

圖 5 左側是主強度圖的參考窗口m，m的窗口大小一般取2的倍數， 如32×32、 64×64、 128×128等， 中間是輔強度圖的搜索窗口s，s的窗口尺寸必須>m， 紅框是與m窗口大小相等的輔匹配窗口s0， 右側為相關系數ρ， 計算的過程為： 保持m不變，s0逐行逐列地在s內移動， 每移動1次， 就利用式(1)計算1次m與s0之間的相關系數ρ， 并記錄此時計算得到的相關系數、 信噪比以及方位向和距離向的偏移量， 依次遍歷整個搜索窗口。最終， 相關系數極值所在的位置對應的偏移量， 即為2幅影像的最佳偏移。

圖 5 互相關計算示意圖Fig. 5 Schematic diagram of cross-correlation calculation.

1.2.4 多項式擬合

通過互相關計算得到所有控制點的偏移量后， 再進行多項式擬合， 計算模型如式(2)所示， 可初步通過設置閾值剔除粗差， 如設置SNR>5或相關系數>0.45， 并對剩下的數據通過最小二乘反演主、 輔影像之間的幾何轉換參數， 刪除超限值， 然后重新計算轉換參數， 不斷迭代計算， 直到剩下的所有數據滿足條件為止。對于主、 輔影像而言， 只有1組轉換參數， 即旋轉參數(m1，m2，m3，m4)與平移參數(m5，m6)， 而(x1，y1)、 (x2，y2)分別為主、 輔影像的像素坐標。

(2)

2 標準配準方法與基于特征點的配準方法的比較

我們首先對特征點和規則點分別在距離向和方位向的偏移量進行了比較。圖6a與6b是經過水體掩膜后整個SAR影像偏移量的統計直方圖。規則點的取樣間隔為256像素×512像素， 為了更加客觀地比較， 設置Hessian閾值為25000， 獲取的特征點數量與規則點的數量接近。由于本文所選的雷達數據不存在變形或偏移量小于配準所能識別精度的情況， 偏移量一般呈現為常數或按照一定規律變化， 本文通過偏移量與0之間的偏差評估算法的準確性。我們以十分之一個像素作為評估指標， 經過水體掩膜后， 在距離向上基于特征點的配準方法獲取的偏移量約有80%介于-0.1～0.1像素之間， 而標準的配準方法僅有約54%； 在方位向上， 分別有82%、 62%的偏移量介于-0.1～0.1像素之間， 可見本文提出的方法所獲得的偏移量更加穩定可靠。

圖 6 特征點和規則點的偏移量及相關系數比較Fig. 6 Comparison of offset and cross correlation between feature points and regular points.a 距離向偏移量統計； b 方位向偏移量統計； c 相關系數統計

互相關系數是互相關質量的重要指標， 因此， 接下來將對規則點以及不同Hessian閾值的特征點所獲取的相關系數進行對比。圖6c為互相關系數的概率分布， 互相關系數越高表示互相關越穩健， 偏移量將更精確。標準配準方法的峰值約為0.4， 本文提出的配準方法的峰值約為0.9， 遠遠高于前者。通常認為相關系數>0.45的配準是相對可靠的， 我們統計了這2種方法相關系數>0.45的同名點數量， 其中標準方法僅有約43%， 而本文提出的方法 >94% ， 具有非常明顯的優勢。此外， 隨著Hessian的閾值從25000增加到34000， 特征點的數量減少， 通常獲取的相關系數整體將提高， 但Hessian增加到達一定的程度后， 相關系數的峰值就不再增大， 所求得的相關系數已無法進一步提高， 而是更多集中于高相關性區間， 這與獲取的SAR影像數據的質量相關。

之后， 為了測試本文提出方法的性能， 我們分別評估了利用特征點與規則點進行多項式擬合的精度與效率。利用式(2)對偏移量數據進行最小二乘估計， 設置后驗標準差為0.05， 如果計算的殘差值超過了后驗標準偏差的1.5倍則將其剔除， 接著利用余下的數據進行迭代計算， 直到所有數據的殘差值符合條件為止。

首先對主強度圖設置256像素×512像素的間隔， 獲取了8821個規則點， 作為比較， 設置Hessian為25000， 獲得8693個特征點。在迭代反演中， 規則點和特征點分別保留了為56%和18%。接下來， 進一步提高Hesssian閾值至34000， 獲取了2164個特征點， 經計算其得到與8821個規則點的精度近似甚至更好， 對主、 輔影像的轉換參數模型的約束同樣好。這說明在配準的過程中， 利用規則點進行互相關計算浪費了大量時間， 且在相同的精度條件下， 特征點的數量和計算時間僅約為規則點的1/4， 可極大提高配準效率。

3 智利地震

3.1 數據

2015年9月16日， 智利Illapel的近海發生了8.3級強烈地震， 且引發了海嘯(2)http://earthquake.usgs.gov/。， 其首都圣地亞哥震感強烈。本文使用的數據是2015年8月26日與2015年9月19日018軌道的Sentinel 1A的SAR影像數據， 通過偏移追蹤研究2015年智利MW8.3地震產生的同震形變場， 并對本文提出的方法進行驗證。

3.2 基于FP的偏移追蹤

配準獲得的偏移量可直接應用于地表形變研究， 即進行密集的配準， 稱為偏移追蹤。常規的偏移追蹤(Regular Points Offset Tracking， RPOT)估計的是均勻分布于2個SAR影像上同名點之間的變形， 不考慮地面的散射特性， 從而導致配準或偏移追蹤的效率低下且結果不可靠。下面對常規偏移追蹤和特征點偏移追蹤(FP Offset Tracking， FPOT)進行比較。RPOT設置規則點間距為32像素×32像素， 并進行偏移量計算， 結果如圖7a、 b所示， FPOT設置Hessian閾值為6000以獲取特征點， 然后進行偏移量計算， 結果如圖7c、 d所示。FPOT和RPOT均屏蔽水體， 通過對比可以發現， RPOT的偏移量結果中含有大量噪點， 特別是圖7a、 b的左下角與右上角， 其中左下角是海洋， 右上角是植被茂密的森林， 獲取的偏移量極不可靠。由于在海洋或森林中極少存在特征點， 因此FPOT能很好地規避這些低相干區域， 可見本文提出的方法適用于植被茂密、 水體覆蓋率大等復雜而又易失相干的區域， 極大提高了計算的效率和可靠性。

圖 7 智利地震初始偏移量散點圖Fig. 7 Scatter plot of initial offset of Chile earthquake.a RPOT距離向偏移量； b RPOT方位向偏移量； c FPOT距離向偏移量； d FPOT方位向偏移量

3.3 實驗結果

為了移除噪點對多項式模型的影響， 將相關系數<0.45或者信噪比<5的同名點予以剔除， 并屏蔽形變近場進行多項式擬合， 獲取偏移量殘差， 最后進行濾波處理。為了能更有效地比較FPOT和RPOT的性能， 在進行水體掩膜后， 對偏移追蹤數據設置相同的參數， 以獲取最終的變形場， 如圖 8 所示。最終， FPOT約保留了75%的數據， 而RPOT僅保留了約30%的數據， 可見前者極大提高了偏移追蹤的計算效率。此外， 與GPS觀測數據進行了對比， 在距離向、 方位向與GPS的均方根RMS分別約為15cm、 19cm， 精度約為 1/15 像素。

4 結論

精確配準是InSAR數據處理的關鍵。標準的配準方法存在冗余計算、 相關性低或可靠性差等缺陷。為了提高標準配準方法的可靠性和效率， 我們提出了一種利用特征點進行配準的方法， 該方法結合了特征點探測、 水體掩膜以及互相關等步驟。通常情況下， 如果地表類型差異大、 強度圖區分明顯， 常規配準算法也能獲得較好的配準效果， 此時配準方法對結果影響較小。然而， 在地表植被茂密或者水體較多的情況下， 配準精度較難提高， 本文旨在解決這方面的問題。通過對珠三角地區的SAR影像數據的配準實驗以及智利地震的偏移追蹤的應用表明， 與標準的配準方法相比， 基于特征點的配準方法適用性廣， 并能實現更高的互相關度， 從而獲得更加可靠的偏移估計， 并能以更小的計算量達到與標準方法相同的精度， 進一步提高配準效率。