基于特征點的SAR影像偏移跟蹤方法研究

彭林才，王華，吳希文

(廣東工業(yè)大學測繪工程系，廣東 廣州 510006)

1 引 言

合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar，InSAR)由于具有覆蓋范圍大、全天候、精度高且不受云層影響的特點，廣泛應用于自然災害的監(jiān)測[1，2]。InSAR的精度可以達到毫米級，但是其最大局限于SAR影像之間的失相干，特別是在大梯度的地表形變區(qū)域中易造成失相干而無法獲取近場形變，如地震[3～5]。偏移跟蹤作為InSAR的重要補充，它對相干性相對不敏感，并且能夠處理更大的可測量梯度。

標準偏移跟蹤，即規(guī)則點偏移跟蹤(Regular Points Offset Tracking，RPOT)是在整個場景中以一定的間隔規(guī)則地選擇同名點，不考慮地面的散射特性，然后進行互相關計算主輔影像的偏移量，最后進行多項式擬合來估計地表位移。該方法在具有較高相干性的區(qū)域效果較好，但在實際應用中，規(guī)則點(Regular Points，RP)的互相關遠非最佳。RPOT具有兩個明顯的缺點：①計算量大；②同名點可能位于失相關區(qū)域。

為了改善這種情況，可以使用點目標進行偏移跟蹤[6，7]，例如永久散射體、孤立點等。盡管永久散射體相對穩(wěn)定，但需要一定數量的SAR影像才能可靠識別，而許多研究通常只有在事件發(fā)生前后的數據可用。另一方面，孤立點是以是孤立的亮點作為同名點[7]，雖然可以有效避免同名點位于失相關的區(qū)域(例如水體)，但由于亮點通常會聚集在一定的空間中，因此可能會導致類似斑狀的偏移估計。

由于偏移跟蹤的主要目的是在主輔影像之間找到相同的點，因此利用具有明顯特征的點將對偏移跟蹤很有幫助。Speeded Up Robust Feature(SURF)是用于遙感影像特征檢測和描述的最流行方法之一[8]。已有研究人員利用SURF進行SAR影像配準[9，10]，但是，單一使用SURF算法進行配準的精度低且誤匹配率高，難以滿足InSAR配準和偏移跟蹤的要求。

本文提出了一種結合SURF、外部地理數據和偏移跟蹤的基于特征點的偏移跟蹤方法(Feature Points Offset Tracking，FPOT)。首先利用SURF算子對主影像的振幅圖進行特征點探測，然后利用外部地理信息(如水體數據)剔除不當的特征點，并利用剩下的特征點進行互相關計算，然后進行多項式擬合以及系統(tǒng)誤差的剔除，最終獲取地表形變。下面通過模擬數據分析以及尼泊爾地震來評估該方法的性能。

2 模擬實驗

2.1 數據

本文選擇2016年11月3日新西蘭凱庫拉地區(qū)獲取的一幅C波段Sentinel-1A的單視復數(Single Look Complex，SLC)影像，該影像距離向和方位向的像素大小分別為 2.33 m×14.13 m，對應于地面距離和方位向的空間分辨率分別為 5 m×23 m[11，12]。本文模擬實驗選擇局部的 5 000×5 000像素作為主影像，并將主影像添加40%的高斯噪聲作為輔影像，進行模擬數據的偏移跟蹤實驗。

2.2 基于特征點的SAR影像偏移跟蹤

本文提出的FPOT方法，首先使用SURF提取主影像的特征點(Feature Points，FP)，然后使用外部土地覆蓋數據剔除無效FP，并對保留的FP進行互相關以獲得偏移量，最后進行多項式擬合以及系統(tǒng)誤差的剔除，獲取地表形變。

(1)SURF提取FP

SURF算法[8]是改進ScaleInvariant Feature Transform(SIFT)[13]的一種特征檢測算法。相比SIFT，SURF具有更快的速度以及更好的魯棒性，其特征點提取主要流程是：

①構造Hessian矩陣

在影像I中的任意一點P=(x，y)，空間尺度為σ，其Hessian定義為：

(1)

(2)

通過Hessian矩陣的行列式可初步判定特征點位置。當行列式值為0時，可判定該點比鄰域內的其他點更加明亮或者暗淡。

②建立影像金字塔

SIFT是保持濾波器模板不變，對原始影像進行降采樣獲得影像金字塔，而SURF保持影像不變，改變盒式濾波器大小以保持尺度不變性，相比較SIFT，具有更快的速度和更高的精度。

③定位特征點

將Hessian矩陣處理的像素點，需與同尺度空間的8個領域像素和相鄰尺度空間的18個像素進行比較，然后去除能量較弱及錯誤定位的特征點，最終選定穩(wěn)定特征點。

為了提高FP識別的速度并抑制噪聲，我們對強度圖的灰度值進行排序，對灰度值大于2.5%且小于97.5%的值線性拉伸到0～255，超出上限(97.5%)和下限(2.5%)的像素分別設置為255和0。主影像經過初步處理后，使用OpenCV的SURF算子來探測FP。為了加快計算速度并防止整個檢測窗口位于水體中，我們將窗口大小設置為 5 000×5 000。

圖1(a)可以明顯觀察到水與陸地之間的幅度差異，這有助于特征點進行準確定位在陸地上。設置Hessian閾值為 1 000和 5 000時，FP的數量和位置存在很大的差異，Hessian閾值為 1 000時，仍有一部分落入水體中(圖1藍色星號所示)，而提高Hessian至 5 000時，其數量明顯少于前者的，但大部分集中于城市區(qū)域(圖1中下部)，僅有極少數的FP落入水體中，因此提高Hessian閾值，可以有效避免FP落入水體中。圖1(b)是在圖1(a)右下角黑框表示的區(qū)域放大圖，紅色圓圈標注了FP的位置，這些FP一般位于拐角點、城市建筑等有明顯特征的地方，很少出現在森林、農田、水域等低相關區(qū)域。圖1(c)是與圖1(b)相對應的相關系數圖，底圖是步長為4×4像素的相關系數。顏色越深相關系數越低，顯而易見，水體中的相關系數基本小于0.1，而FP大多數都位于高度相關的區(qū)域中。

圖1 模擬數據的特征點

(2)利用FP進行偏移跟蹤

在大面積水域中進行互相關，相關度極低且偏移量非常不穩(wěn)定，這不僅大大降低了計算效率，而且還導致錯誤的偏移估計。盡管SURF可以很好地避開水體，但是當Hessian閾值設置得太低或探測窗口太小，或者整個探測窗口都在水中時，SURF探測到的FP可能位于水中。因此，在進行偏移跟蹤之前，需利用外部土地覆蓋數據剔除落入水中的FP。

FPOT與RPOT類似，只是選擇同名點的方法不同。在模擬實驗中，對于特征點探測，首先設置Hessian閾值為 3 000，Octaves為4，Octaves layers為3，使用SURF算子獲得了 58 832個特征點。對于規(guī)則點，為了減少計算時間，距離向和方位向設置步長為16×16像素，獲得有效配準點 96 568個，RP的數量幾乎比FP的多一倍，因此其計算量也隨之成倍增加。設置主影像匹配窗口為64×64像素，搜索窗口為84×84像素，過采樣因子為128。不使用任何外部數據進行掩膜，經過互相關計算，RPOT獲得的距離向偏移量如圖2(a)，FPOT獲取的距離向偏移量如圖2(b)。在圖2(a)中，由于SAR影像數據所對應的區(qū)域含大量水體，導致圖上側和右側含有大量噪聲點，因此，RPOT必須經過水體掩膜以剔除粗差。在水體掩膜后，RP剩下 64 913個，濾除了約33%位于水中的RP，FP則保留了 58 299，僅濾除了約0.9%位于水中的FP，在不利用水體數據的情況下也能獲得良好的結果，且噪聲點要遠小于規(guī)則點。因此，使用特征點進行偏移跟蹤，能有效避開水體，提高偏移跟蹤的效率和魯棒性。

圖2 模擬數據距離向偏移量

(3)FPOT和RPOT的比較

互相關系數是互相關質量的重要指標，我們首先對FP和RP的偏移量以及相關系數進行了比較。圖3(a)和圖3(b)是在水體掩膜之后整個SAR影像的偏移量統(tǒng)計直方圖，而偏移量與0之間的偏差可以評估算法的準確性。我們以二十分之一像素作為評估指標。在距離向上，FPOT和RPOT分別有97%、91%的偏移量介于-0.05～0.05像素之間；在方位向上，分別有94%、85%的偏移量介于-0.05～0.05像素之間。FPOT和RPOT在距離向上的標準偏差分別為0.022和0.028像素，在方位向上的標準偏差分別為0.026和0.035像素，FPOT的標準偏差優(yōu)于RPOT。圖3(c)為互相關系數的概率分布，RPOT的峰值約為0.5，FPOT的峰值約為0.6。隨著Hessian閾值的增加，FP的數量減少。當Hessian閾值由 1 000增加到 3 000時，相關系數有了整體的提高，但是再增加Hessian閾值，相關系數的峰值不再改變，因為模擬的數據中添加了大量的高斯噪聲，所求得的相關系數已無法進一步提高，而是更多集中于高相關性區(qū)間。

圖3 FP和RP的偏移量及相關系數比較

此外，本文測試了FP的配準的性能，分別評估了使用FP和RP進行多項式擬合的準確性。我們使用最小二乘估計進行迭代計算。如果殘差大于后驗標準偏差的1.5倍，則刪除異常值，再次使用剩下的點進行最小二乘估計迭代計算，直到收斂。

我們在整個場景中設置步長為128×128像素，獲得了 1 443個RP，并選擇相同數量的響應值最高的FP。在迭代反演中，FP和RP的分別剔除了為26%和55%。其次，我們再次選取響應值最高的前360個FP，經計算得到與 1 443個RP的精度接近，對仿射矩陣的約束同樣好。這說明在配準的過程中RP互相關浪費了大量時間，并且在同等精度下，FP的數量和計算時間僅約為RP的1/4，這可以進一步提高配準效率。

3 尼泊爾地震

3.1 實驗數據

2015年4月25日下午2點，尼泊爾Gorkha地區(qū)發(fā)生8.1級地震，震中為北緯28.2°，東經84.7°，這次地震致使 8 000多人死亡(主要發(fā)生在尼泊爾首都Kathmandu和鄰近城市地區(qū))以及2萬多人受傷。尼泊爾這個國家的大部分地區(qū)是位于喜馬拉雅山南麓，這些區(qū)域地形起伏大，海拔高，常年積雪覆蓋，而其首都Kathmandu則坐落于一個較小的斷陷盆地之上，而其整體則處于喜馬拉雅斷裂帶的上盤。這種較為特殊的地形以及地質地貌情況會致使雷達信號削弱，容易導致失相干，使得InSAR難以在該區(qū)域進行形變觀測，因此對該地區(qū)的地震事件前后的SAR影像數據進行偏移跟蹤以獲得地震近場形變。在2015年的尼泊爾地震中，使用了C波段的Sentinel_1A數據，其數據參數如表1所示。

尼泊爾地震Sentinel_1A數據參數 表1

3.2 基于特征點的偏移跟蹤

在進行偏移跟蹤前，選取主影像的振幅圖進行特征點探測。為了在地震近場具有更多的可用數據的同時減輕計算負擔，設置地震近場的Hessian閾值為 6 000，在遠場則設置為 10 000，獲取了 265 090個特征點。從圖4的振幅圖所對應的實際地形地貌中，左上角和右下角是密集覆蓋的森林，由于森林中較少明顯的特征，因此該區(qū)域獲取的特征點非常稀疏，而在喜馬拉雅山脈中，具有非常明顯山脊線的特征，因此該區(qū)域獲取的特征點較為密集，而中部是Kathmandu城市區(qū)域以及裸露的高山地區(qū)，特征點非常密集，而震中恰好在此區(qū)域，這對于特征點的探測以及偏移跟蹤的計算都大有裨益。

圖4 尼泊爾特征點分布圖

在進行水體掩膜后進行互相關計算，匹配窗口和搜索窗口分別設置為128×128像素和148×148像素，過采樣因子為128。FPOT獲取的偏移量如圖5(a)、(b)所示。作為對比，設置RP的步長為64×64像素，獲取的有效配準點數為 311 112，RPOT獲得的獲取的偏移量如圖5(c)、(d)所示。從圖5可看出偏移量的噪聲都很大，雖然RPOT的數量占優(yōu)，但其具有更多且更密集的噪聲點，其主要原因是整個場景覆蓋范圍內含有大面積的森林，由于其散射特性，在植被茂密的區(qū)域進行偏移跟蹤，其相關性和信噪比可能會很低，且獲得偏移量含有大量的隨機噪聲，而基于該測量結果進行建模獲取地表形變是不可靠的，因此需要進一步的濾波降噪處理。

圖5 尼泊爾地震偏移量散點圖

針對偏移跟蹤獲取的偏移量含有大量噪聲點的情況，在粗差探測后，屏蔽形變近場進行多項式擬合并獲取偏移量殘差，最后進行濾波處理。為了能更有效比較FPOT和RPOT性能，在偏移跟蹤數據處理中設置相同的參數，表2展示了經粗差探測及濾波處理后保留的點數占初始總點數的百分比，可以看出，在相同的條件下，FPOT保留的點數多于RPOT，這進一步證實了FPOT在性能和質量上要優(yōu)于RPOT。

尼泊爾地震偏移跟蹤數據統(tǒng)計 表2

3.3 FPOT與GPS對比驗證

FPOT獲取的偏移量采樣間隔設置為 550 m×550 m(對應于經緯度0.005°×0.005°)，圖6顯示了與尼泊爾地震有關FPOT衍生的近場的距離向和方位向位移及其與GPS觀測結果對比的散點圖，黑色圓圈表示GPS測站位置，其位移值以相應的顏色進行填充，圖中顯示震中近場的形變非常明顯?？梢悦黠@看出距離向的最大位移 2 m，而方位向的最大位移約 2 m～ 3 m。將GPS三維形變數據轉換到距離向和方位向后，與FPOT的結果對比，RMS大約為 10 cm，在總體上與FPOT的結果一致，且由于分辨率的影響，距離向的RMS優(yōu)于方位向。

圖6 尼泊爾地震形變場

4 結 論

本文通過模擬實驗以及尼泊爾地震測試了FPOT的性能，實驗結果表明，FPOT要優(yōu)于RPOT。使用SURF探測FP能夠有效避免位于水體等失相關區(qū)域，并且在FP上的振幅相關度高，從而獲得更可靠的偏移估計，獲得的結果與GPS數據進行了比較，顯示出良好的一致性。此外，在配準過程中，利用特征點，僅需要規(guī)則點數量的1/4即可達到相同的精度，而在偏移跟蹤中，FPOT保留的點數約為RPOT的2.5倍。因此，基于FP的可以顯著提高偏移跟蹤的效率和可靠性。

致謝：感謝歐洲空間局Sentinel-1A衛(wèi)星數據的免費使用。