楊士超，查顯杰，郭曉丹，張志宏，張 琪

（1.遼寧省地震局，遼寧 沈陽 110034；2.中國科學(xué)技術(shù)大學(xué)，安徽 合肥 230026）

0 引言

大地震地表破裂的觀測資料對于研究地震力學(xué)和斷層區(qū)域結(jié)構(gòu)性質(zhì)都具有重要意義。斷層幾何結(jié)構(gòu)和相關(guān)滑移分布也是決定地震動力學(xué)破裂過程的關(guān)鍵所在[1]。利用地表形變位移場可以反演斷層的滑移分布，研究斷層結(jié)構(gòu)和斷層破裂狀態(tài)；同時也可以將地表位移場作為邊界條件，對地震斷裂運動進(jìn)行約束。因此對于研究地震動力學(xué)和運動學(xué)問題來說，獲取大地震同震位移場都是至關(guān)重要的科學(xué)問題。

目前，地表形變監(jiān)測技術(shù)主要包括野外地質(zhì)調(diào)查[2]、GPS 臺站測量[3]、合成孔徑雷達(dá)干涉（InSAR）測量[4]、子像素偏移追蹤測量[5-6]。對于廣泛區(qū)域的偏移場獲取，野外地質(zhì)調(diào)查需要大量的人力和物力，而偏遠(yuǎn)地區(qū)GPS 臺站覆蓋有限，GPS 技術(shù)很難有效發(fā)揮作用。InSAR 是一種理想的大地震近場形變監(jiān)測技術(shù)。然而，大地震的近場形變變化梯往往較大，斷層附近的InSAR 形變干涉圖條紋常常非常密集（即出現(xiàn)信號飽和現(xiàn)象），因此InSAR 方法獲取的斷層附近區(qū)域的位移場精度會顯著降低。子像素偏移追蹤算法主要是通過處理地震前和地震后獲取的兩景光學(xué)遙感影像，采用互相關(guān)技術(shù)得到大地震的水平位移場[7-9]。該方法不存在形變信號飽和的問題，因此它非常適合監(jiān)測大地震的地表形變。本研究小組在強(qiáng)度互相關(guān)、雙線性插值和二維高斯回歸三種技術(shù)基礎(chǔ)上新發(fā)展了一套子像素偏移追蹤算法[10]。與現(xiàn)有算法相比，該算法能以更高的精度從遙感影像對中獲取位移場。為了驗證算法的精度和可靠性，本文利用新發(fā)展的子像素偏移追蹤算法對2014 年2 月12 日新疆于田地震和2013 年9 月24 日巴基斯坦地震觀測資料進(jìn)行處理，獲取這兩次地震的同震形變位移場。

1 算法原理與實現(xiàn)

1.1 算法原理

利用光學(xué)遙感影像求解形變位移場的精度取決于影像的分辨率和子像素偏移追蹤算法的精度。由于遙感成像技術(shù)的限制，所獲取的光學(xué)影像分辨率都是一定的（通常商用衛(wèi)星數(shù)據(jù)空間分辨率為米級或亞米級），因此需要發(fā)展一種能獲取高精度位移場的子像素偏移追蹤算法。目前，子像素偏移追蹤算法主要有兩種：一種是頻域互相關(guān)算法，例如COSI-Corr 程序包[11]；另一種是空間域互相關(guān)算法，例如ImCorr 程序包[12]。

1.1.1 頻域互相關(guān)算法

頻域互相關(guān)又稱作相位互相關(guān)，這種算法主要基于傅里葉偏移原理[13]?？臻g域內(nèi)兩個函數(shù)存在的偏移值可以轉(zhuǎn)變?yōu)轭l域內(nèi)線性的相位差?？梢杂萌缦鹿奖磉_(dá)：

假設(shè)f1（x，y）和f2（x，y）分別表示兩個圖像的函數(shù)，（x0，y0）表示f1和f2之間的相對偏移值，那么

根據(jù)傅里葉偏移原理可得：

其中，F(xiàn)1（u，v）和F2（u，v）分別表示f1（x，y）和f2（x，y）的傅里葉變換。因此，標(biāo)準(zhǔn)互功率譜可以表示為：

其中*表示復(fù)共軛。

通過公式（3）求解偏移值（x0，y0）有兩種方式：一種方式是直接在頻率域內(nèi)求解，擬合相位平面，相位平面的周期性使得這種方法很難實現(xiàn)；另一種方式是將標(biāo)準(zhǔn)互功率譜進(jìn)行傅里葉反變換，便可以得到δ（x-x0，y-y0），通過找到狄拉克δ 函數(shù)的峰值位置來確定偏移值（x0，y0）[14-15]。后一種方法只能搜索到整數(shù)偏移值，要求得更高精度的偏移值還需要借助曲面擬合算法來求解。

1.1.2 空間域互相關(guān)算法

空間域互相關(guān)算法又稱作強(qiáng)度互相關(guān)算法，通過查找主像與從像的相關(guān)系數(shù)矩陣的峰值位置來確定兩幅影像的相對偏移值。相關(guān)系數(shù)的表達(dá)方式有很多種，因為零均值標(biāo)準(zhǔn)化相關(guān)系數(shù)（ZNCC）具有良好的抗噪性，對于偏移和亮度的線性變化不敏感[16]，本文選擇ZNCC 算法計算相關(guān)系數(shù)。ZNCC 相關(guān)系數(shù)表示為：

其中|*|表示絕對值，m（i，j）和s（i，j）分別表示主像與從像第i 行j 列的像素灰度值，M和N 表示所選像對的大小，分別表示m（i，j）和s（i，j）的均值。

利用以上算法可以求得一個相關(guān)系數(shù)矩陣，通過搜索峰值位置可以找到整數(shù)倍偏移值。接下來需要通過插值和曲面擬合來計算得到小數(shù)倍偏移值。一般插值算法可以提高影像的分辨率，但是高倍插值計算會耗費大量內(nèi)存空間和運算時間。本文采用了一種基于雙線性插值的相關(guān)系數(shù)插值方法[17]。這種算法可以將估算子像素偏移問題轉(zhuǎn)化成一個關(guān)于偏移值的非線性反演問題，通過尋找最優(yōu)解來確定偏移值。為了獲取更高精度的偏移，本文算法進(jìn)一步選取相關(guān)系數(shù)矩陣峰值周圍的9 個點，對這9 個點進(jìn)行9 點高斯回歸求取微小偏移值[18]：

其中，C00，C01，C10，C11，C02和C20為通過回歸分析推導(dǎo)出的6 個回歸系數(shù)。

1.2 本文算法的實現(xiàn)

本文所驗證的子像素偏移追蹤算法主要基于上面提到的強(qiáng)度互相關(guān)算法、雙線性插值算法和二維高斯回歸算法，把這3 種算法巧妙的結(jié)合到一起來進(jìn)行偏移計算。該算法的實現(xiàn)由以下5 步組成：

（1）利用兩個大小一樣的滑動窗在地震前和地震后獲取的兩景影像相應(yīng)位置選取一組像片對，再用一個較小的查找窗分別放在主像片中心，用另一個一樣大小的查找窗在從像片中滑動，這樣可以得到像對矩陣。采用強(qiáng)度互相關(guān)算法計算出一個相關(guān)系數(shù)矩陣；

（2）在相關(guān)系數(shù)矩陣中查找峰值位置，即整數(shù)像素的偏移；

（3）再將從像片中查找窗移至峰值點位置，并利用相關(guān)系數(shù)插值算法在峰值周圍對相關(guān)系數(shù)平面進(jìn)行插值；

（4）在插值的相關(guān)系數(shù)矩陣中查找峰值位置，利用高斯回歸方法獲取子像素偏移量，即子像素偏移；

（5）將整數(shù)像素偏移和子像素偏移加在一起，構(gòu)成位移場。

由于兩景影像的成像時間存在差異，大氣和光照的影響會使兩景影像視覺效果上存在一定差異，降低影像對的相關(guān)性。并且在沒有高精度的DEM 數(shù)據(jù)的情況下，地理坐標(biāo)系的定位也可能存在差異，所以有必要在進(jìn)行偏移場計算之前，對影像對進(jìn)行預(yù)處理來消除這些影響。由于云、冰蓋等因素的影響，部分區(qū)域的位移值可信度較低，需要利用峰值相關(guān)系數(shù)和相關(guān)系數(shù)比值來除去這些可靠性低的點。本文利用峰值相關(guān)系數(shù)和相關(guān)系數(shù)比值設(shè)置閾值，對云或冰蓋覆蓋區(qū)域進(jìn)行掩膜處理。經(jīng)過掩膜處理后的位移場仍然存在著由于成像過程和預(yù)處理產(chǎn)生的長波長畸變，長波長畸變需要采用適當(dāng)?shù)暮筇幚硭惴▉硐?。?jīng)過以上處理便可得到平滑、正確的位移場。

2 應(yīng)用實例

2.1 2014 年新疆于田地震

2014 年2 月12 日17 時19 分，我國新疆維吾爾自治區(qū)和田地區(qū)于田縣發(fā)生強(qiáng)烈地震，據(jù)中國地震臺網(wǎng)測定，震中位于36.1°N、82.5°E，震源深度12 km，震級Ms7.3（CENC）。此次于田地震震中位于阿爾金斷裂西南段尾端的分支斷裂的交匯處，為拉張與擠壓變形的過渡帶（圖1），地質(zhì)構(gòu)造復(fù)雜，具有發(fā)生強(qiáng)震的背景[19]。經(jīng)野外調(diào)查顯示，這次地震在高海拔地區(qū)形成了由一系列張裂隙、張剪裂隙、剪切裂隙以及擠壓鼓包和裂陷等斜列狀組合而成的地表破裂帶,整體呈NEE—SWW 走向，全長約25 km，分為不連續(xù)的兩段，顯示出左旋走滑伴隨有正滑分量的特征，最大左旋位移約1 m[20]。

圖1 2014 年2 月12 日于田地震構(gòu)造背景及震中位置Fig.1 Construction background and epicenter location of Yutian earthquake on Feb.12，2014

圖2 2014 年2 月12 日于田地震預(yù)處理后的光學(xué)遙感影像對Fig.2 Preprocessed optical remote sensing image pair of Yutian earthquake on Feb.12，2014

圖3 影像對的峰值相關(guān)系數(shù)圖Fig.3 Peak correlation coefficient map of image pairs

本文從美國地質(zhì)調(diào)查局USGS 上選取了2013 年5 月31 日和2014 年5 月18 日的兩景Landsat TM8 遙感數(shù)據(jù)分別作為監(jiān)測2014 年2月12 日于田地震同震位移場的震前圖像和震后圖像，并選用圖像分辨率為30 m 的第七波段。利用ENVI 軟件，將兩景圖像進(jìn)行影像增強(qiáng)、初配準(zhǔn)、紋理增強(qiáng)、影像剪裁和灰度值匹配的預(yù)處理，獲得一對7000×7000 像素的影像。圖2a 和2b 分別為預(yù)處理之后的震前圖像和震后圖像。之后我們選用256×256、300×300、360×360、512×512、600×600 等像素大小的查找窗口，在行和列方向的采樣間隔均為100 像素，利用上述子像素偏移追蹤算法對預(yù)處理后的影像對進(jìn)行計算得到位移場。圖3 為峰值相關(guān)系數(shù)圖。本文設(shè)置相關(guān)系數(shù)閾值和相關(guān)系數(shù)比值對偏移場掩膜，獲取了2014 年于田地震同震偏移場（圖4）。本文利用遠(yuǎn)場選取了部分位移場數(shù)據(jù)分析了算法的統(tǒng)計誤差。經(jīng)計算，東西向位移場誤差為0.4170 m，南北向誤差為0.4073 m，精度為百分之幾個像素（圖5）。

圖4 2014 年2 月12 日于田地震同震位移場Fig.4 Displacement fields of 2014 Yutian earthquake

從我們得到的同震偏移場（圖4）可以看出，兩條黑色線段上部呈北西向運動，黑色線段下部呈南東向運動，兩個條帶存在左旋走滑趨勢的斷層，分別為北硝爾庫勒斷裂和阿什庫勒—硝爾庫勒斷裂，斷裂以左旋走滑為主，帶有一定量的正滑分量。北硝爾庫勒斷裂和阿什庫勒—硝爾庫勒斷裂全長約20 km，整體走向約N60°E，最大偏移量可達(dá)到約1.5 m，與李海兵等（2014）野外地質(zhì)調(diào)查結(jié)果一致。震源南方地貌復(fù)雜，不易調(diào)查，本文得到的位移場表明該區(qū)域存在明顯的地震地表破裂。

圖5 2014 年2 月12 日于田地震同震位移場誤差圖東西向誤差Fig.5 Error map of 2014 Yutian earthquake displacement fields

隨著印度板塊的向北俯沖，青藏高原沿著郭扎錯斷裂、阿爾金斷裂中段和東段向北東方向推擠。阿爾金斷裂在西南尾端分為兩支：一支為康西瓦斷裂呈西南走向延伸；另一支為郭扎錯斷裂整體呈現(xiàn)北東向走滑形式，并帶有張裂特征[21]。2008 年于田Ms7.3 級地震的破裂帶主要呈NS—NNE 向，位于2014 年于田地震的SWW 方向約120 km（圖1），該斷裂為左旋張裂破裂帶[22]。2014 年于田地震正處于三條斷裂所夾的三角區(qū)域內(nèi)，可能是在這些影響下發(fā)生了2014 年于田地震。

2.2 2013 年巴基斯坦地震

2013 年9 月24 日的巴基斯坦西南地區(qū)發(fā)生Mw7.7 地震，經(jīng)Global CMT 定位地震震源位于26.70°N、65.04°E，震源深度12 km。地震位于歐亞板塊、印度板塊和阿拉伯板塊的交界處附近，該次地震也是板塊間的相互作用引起的。

圖6 2013 年9 月24 日巴基斯坦地震預(yù)處理后的光學(xué)遙感影像對Fig.6 Preprocessed optical remote sensing image pair of Pakistan earthquake on Sep.24，2013

與上面的例子相同，從美國地質(zhì)調(diào)查局USGS 上選取了2013 年9 月10 日和2013 年9月26 日的兩景Landsat TM8 遙感數(shù)據(jù)分別作為監(jiān)測2013 年9 月24 日于田地震同震偏移場的震前圖像和震后圖像（圖6），并選用圖像分辨率為15 m 的第八波段。由于單景影像無法覆蓋斷裂區(qū)域，我們選取了兩個連續(xù)景的光學(xué)遙感數(shù)據(jù)構(gòu)成了兩個影像對。經(jīng)預(yù)處理，再分別求出位移場，最后將兩個位移場繪制到一起構(gòu)成一幅覆蓋2013 年巴基斯坦地震震中區(qū)域的位移圖。本文得到的2013 年9 月巴基斯坦地震同震偏移場如圖7 所示。同樣，我們在遠(yuǎn)場處選取了部分位移數(shù)據(jù)，統(tǒng)計分析位移計算誤差。結(jié)果表明東西向偏移場誤差為0.5318 m，南北向誤差為0.3931 m，精度為百分之三個像素（圖8）。

圖7 2013 年9 月26 日巴基斯坦地震同震形變位移場Fig.7 Displacement fields of Pakistan earthquake on Sep.26，2013

所得同震偏移場給出，該地震斷裂為左旋走滑型的圓弧形斷裂，基本無正斷分量，平均滑移量6 m，最大可達(dá)到10 m，并且斷層西北側(cè)的滑移量明顯大于東南側(cè)，結(jié)果與Barnhart（2014）、Avouac（2014）的研究結(jié)果一致[23-24]。

圖8 求取的2013 年9 月26 日巴基斯坦地震同震位移場的誤差Fig.8 Error map of Pakistan earthquake displacement fields on Sep.26，2013

3 討論

本文利用一種新子像素偏移追蹤算法計算了2013 年巴基斯坦和2014 年新疆于田兩次地震的同震形變場。該子像素偏移追蹤算法將強(qiáng)度互相關(guān)算法、雙線性插值算法和二維高斯回歸算法結(jié)合到一起。應(yīng)用實例表明，該子像素偏移追蹤算法測量精度可以在百分之三個像素水平上，可以很好的運用到求取中強(qiáng)震的同震形變場中。該算法精度略高于目前被廣泛使用的ImCorr 和COSI-Corr 算法。由于本文算法內(nèi)部使用強(qiáng)度互相關(guān)，因此該算法可以使用任意大小的查找窗口進(jìn)行計算。該算法在后處理、遙感影像的條帶狀殘差、兩幅影像合并過程等方面仍需要改進(jìn)。在中強(qiáng)震的偏移場監(jiān)測中，可以與InSAR、GPS 相結(jié)合，進(jìn)行互補(bǔ)，對該算法進(jìn)行優(yōu)化。