光學影像匹配技術在地表形變監測中的應用研究

張麗

(上海市測繪院，上海 200129)

1 引 言

地表形變是在地殼表層發生的一種空間位移變化，由于地殼運動、火山運動、地震活動等自然因素，或者是地下水開采、礦山開采、城市建設等人為因素引發。近幾年地震頻發，有專家預測全球開始進入板塊運動活躍期，因此，對地表進行變形監測顯得尤為重要。

長期以來，研究者主要依靠傳統的地面水準測量和GPS監測網來進行監測，其監測精度雖可高達亞毫米級，但不適宜大區域的長期重復監測，且要消耗大量人力物力。近年來，隨著SAR技術的發展，其全天時、全天候工作，且可以大范圍、高精度地從空間獲取地形信息和地表微小的形變信息的優勢彌補了這一缺陷；然而該技術容易受相位去相關影響，在形變過大區域常常無法監測，且只能獲取沿雷達視線向即距離向的地表形變。

而光學影像匹配技術是借助遙感影像的幅度信息，采用互相關技術尋找兩幅圖像的同名點，從而獲得同名點處對應像素在方位向和距離向的偏移量。由于避免了相位運算，該方法能很好地彌補了受相位去相關影響這一不足。因此，本文將該方法應用到具體地表形變監測案例中，并研究其可行性。

2 光學影像匹配

2.1 光學影像匹配技術原理

影像匹配技術的根本就是圖像配準，在光學影像里，圖像配準也稱為幾何精校正，它規避了成像的空間幾何過程，直接利用地面控制點(GCP，Ground Control Points)數據，對衛星影像的幾何畸變過程進行數學模擬，利用已知條件確定模型中的位置參數，建立原始的畸變圖像空間與校正空間之間的相互對應關系，然后利用這種對應關系把畸變圖像空間中的全部元素變換到校正圖像空間中去，從而實現幾何精校正。精校正后，光學遙感影像間的幅度信息仍存在差異，通過傅立葉變換，在頻率域內采用高精度的相位相關技術即可提取偏移量。

2.2 光學影像匹配技術流程

透過原理可以看出，光學影像匹配技術主要分為正射校正和相位相關計算兩步，本文所采用的數據處理流程大致包括：

①利用外部DEM數據及影像的太陽高度角，方位角模擬生成DEM陰影圖；

②從震前影像與DEM陰影圖上選取同名點，對震前影像進行正射糾正；

③從震后影像與已校正的震前影像上選取同名點，對震后影像進行正射糾正；

④對已校正后的震前震后影像做相位相關計算，獲取有效形變值。

(1)選取同名點

在上述兩次配準過程中均涉及同名點的選取，通常選取同名點分為人工選點和自動選點兩種方法。本文采取人工選點法，即通過目視在兩幅待配準影像上，手動選擇若干同名點。同名點選取時需要遵循以下原則：

①至少選3個點，若待糾正影像和一幅色彩較暗的DEM配準時需要選15～30個點；

②同名點要連續選擇，但需遠離形變區域，若整幅圖都有形變就在形變最小處選點；

③點選取時不能太靠近圖像邊界；

④同名點周圍最好有容易識別的地物，最好選擇無須校正的人工物。

(2)影像正射校正與重采樣

正射校正是對圖像空間和幾何畸變進行校正生成多中心投影平面正射圖像的處理過程。它除了能糾正一般系統因素產生的幾何畸變外，還可以消除地形引起的幾何畸變。采用少量的地面控制點與相機或衛星成像模型相結合，確立相機或傳感器平臺、圖像和地面三者之間的空間位置關系，建立正確的校正公式，產生精確的正射影像。

影像重采樣是影像數據重新組織過程中的灰度處理方法，影像采樣是按一定的間隔采集影像灰度數值，當閾值不位于采樣點上的原始函數的數值時，就需要利用已采樣點進行內插，稱為重采樣。本文所采用的重采樣法為三次卷積內插法(也叫辛克插值法)，使用內插點周圍的16個像元值進行距離加權計算柵格值，先在Y方向內插四次(或X方向)，再在X方向(或Y方向)內插四次，最終得到該像元的柵格值。

(3)相位相關法估算偏移量

在經過正射校正和重采樣后，兩幅影像應該具有相同的分辨率，統一的坐標系，但因為存在地表形變，兩幅影像并不能完全重合。因此，還需要利用COSI-CORR軟件模塊提供的高精度相位相關算法進行偏移量估計。該模塊包含了兩種算法，一種是基于頻率域(Frequential)變換的方法，另一種是基于統計學(Statistical)理論的方法。相位相關法是用于配準圖像進行平移變換的典型方法，它主要是依賴于傅立葉偏移理論。

假設有兩幅影像i1和i2，用(△x，△y)來代表兩幅影像上對應點的形變差異，則有公式：

i2(x，y)=i1(x-△x，y-△y)

(1)

用I1和I2分別代表兩幅影像的傅立葉變換，則有以下關系式：

I2(wx，wy)=I1(wx，wy)e-j(wx△x+wy△y)

(2)

其中wx和wy分別代表行和列的頻率變化。定義兩幅影像的互能量譜為如下：

(3)

上式中，C(wx，wy)是i1(x，y)和i2(x，y)之間的互能量譜，*代表復數共軛相乘。

影像間的相對形變量可以通過互能量譜的二維斜率來計算，因此對式(3)進行傅立葉逆變換：

F-1{ej(wx△x+wy△y)}=δ(wx△x+wy△y)

(4)

而在應用中，C(wx，wy)僅為兩幅影像間理論上的互能量譜，由實際計算得到的互能量譜為Q(wx，wy)。定義一個目標函數為φ(△x，△y)：

(5)

其中，W(wx，wy)表示加權矩陣。目標函數φ(△x，△y)在其他位置都為零，只有在平移位置處不為零，因此這個位置就是影像間的相對偏移，當目標函數φ(△x，△y)達到最小值時的(△x，△y)就是所求的偏移量。

3 實驗數據處理

北京時間2010年4月14日，在我國青海省玉樹縣發生了震級為里氏7.1級的大地震。此次地震的發震構造是青藏高原東部規模巨大的甘孜-玉樹斷裂帶，該斷裂帶長度約 500 km，總體走向呈285°～315°弧形展布，是一條大型左旋走滑斷裂帶。

3.1 實驗數據參數

本實驗選用了兩幅SPOT5衛星影像。影像數據采集的時間分別為2009年11月5日和2010年4月15日，影像數據具體參數如表1所示。選用時間相距較近的衛星影像可以最大限度排除人工建筑等因素的影響，協調時間、太陽方位角、入射角的一致，確保了不會由于角度問題導致正射校正的偏差。

實驗選用的SPOT5衛星影像參數 表1

續表1

選用的SRTM數據絕對高程精度為 16 m，相對高程精度為 10 m左右，坐標系為WGS84地理坐標系，投影方式為UTM投影。

3.2 實驗結果分析

本文利用基于ENVI平臺開發的COSI-CORR軟件模塊，按照上述技術路線進行數據處理。

(1)生成DEM陰影圖

由于已有SRTM數據與SPOT5影像分辨率存在較大差異，數據處理的第一步是先將DEM插值到空間分辨率 10 m，然后根據太陽高度角和太陽方位角生成一幅DEM陰影圖，將其作為參考影像，用于震前影像的正射校正中，如圖1所示。

(2)SPOT5影像匹配

兩次配準時，在SPOT影像和DEM陰影圖上，通過人工選點法分別選取了18個和19個均勻分布的同名點。在選取過4個及以上同名點后，軟件模塊會自動預測后續同名點位置，為提高配準精度，在經過重新選取和刪除誤差較大的點后，保證兩次的均方根誤為分別為0.48和0.43，均小于1個像元。

圖1 生成的DEM陰影圖

選取同名點后，需將其轉換為GCP點并進行優化，才能用于后續校正與重采樣。而COSI-CORR軟件模塊提供兩種優化方法(相位相關法和統計相關法)，可以使同名點間配準誤差值達到分米級甚至是厘米級。第一次正射校正時，由于參考影像為DEM陰影圖，其影像特點與SPOT影像存在很大差異，因而選用統計相關法。而第二次正射校正時，由于參考影像為校正過后的震前影像，兩幅影像特點相一致，因此選用相位相關法來提高計算精度。兩次優化的結果表明，選取5次迭代計算能夠使誤配準值即殘差偏移量達到收斂。點優化過程的迭代結果如圖2所示。

圖2 5次迭代后GCP點處的殘差偏移量均值和標準差

(3)影像校正與重采樣

當GCP點優化結果達到收斂且滿足精度要求后，影像進行正射校正與重采樣。此模塊采用的正射校正模型為嚴格軌道模型(Pushroom Sensor)，校正過程利用頭文件中衛星狀態參數提高幾何校正的精度。在選用精確度最高的辛克插值法進行重采樣后，得到校正后的震前、震后影像如圖3和圖4所示。

圖3 校正后的震前影像

圖4 校正后的震后影像

(4)相位相關計算偏移量

本文經過反復對比，進行相位相關法計算時，選用滑動窗口為64×64的單一型窗口，即X和Y方向的Initial和final都選擇64；滑動窗口在影像X和Y方向移動的像素數都采用默認值16；迭代次數選擇2，掩模閾值選擇0.9。獲得了此次地震東西向和南北向的同震位移量，結果如圖5所示。此一步可輸出偏移量文件，共包含三個波段數據，分別為地表東西向形變值，南北向形變值及信噪比SNR值。

圖5玉樹地震同震位移圖(左為東西向形變，右為南北向形變)

圖6 去噪后玉樹地震局部同震位移分布圖

(5)形變結果優化

通過對結果的SNR值(0效果差，1效果優)分析，得到了同震位移場的質量評價。圖6中，上圖為震前SPOT5衛星影像，紅線為同震地表破裂帶，中圖為東西向的同震位移分布圖，下圖為南北向的同震位移分布圖。從圖中可以看到，雖有一些效果比較差的地方，但主要分布在遠離地表破裂帶處，在破裂帶展布區域分布很少。這是由于玉樹地處高原地區，地形起伏比較大，在一些地區得到的數據結果不是很理想，而且由于震后SPOT5衛星影像為第二天應急拍攝，沒有考慮云量，導致震后影像云量比較多，但幸運的是云量基本沒有覆蓋地表破裂帶的范圍，而且地形起伏導致的不理想數據結果也很少在地表破裂帶上，這保證了得到的同震位移場至少在近地表破裂帶處是比較可信的。

從圖中可以明顯地看出地表破裂帶的展布以及同震位移量的分布。此次玉樹地震是左旋走滑型地震，產生的同震位移主要為水平位移，垂直位移分量比較小。其最大位移不超過 2 m，地表裂長度為 30 km左右。

4 結 語

本文對光學影像匹配技術進行了初步研究，利用光學影像空間分辨率高的特點，獲取了玉樹地震水平方向的形變場。根據野外地表破裂帶調查，此次地震較清晰的地表破裂帶由3條主破裂左階組成，總長約 31 km，左旋走滑性質，最大左旋位移量約為 1.8 m。這與本文得到的破裂帶長度和最大同震位移量基本一致，這也從另一個角度說明了實驗得到的同震位移場的合理性和正確性，顯示了光學影像匹配技術在監測地表形變方面的巨大應用價值。