亞像素相位相關(guān)法在獲取汶川地震近場(chǎng)形變中的應(yīng)用

劉云華 單新建 屈春燕 宋小剛 張桂芳

(中國(guó)地震局地質(zhì)研究所，地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029)

0 引言

地震是一種以力學(xué)運(yùn)動(dòng)為主體的自然現(xiàn)象。因此，測(cè)量并掌握地殼的力學(xué)狀態(tài)與力學(xué)運(yùn)動(dòng)，是地震學(xué)與地震預(yù)報(bào)研究中最為重要的基礎(chǔ)性工作(國(guó)家地震局科技監(jiān)測(cè)司，1995)。獲取大范圍的地表連續(xù)覆蓋的同震位移形變場(chǎng)圖像，并利用這些地表同震形變場(chǎng)來揭示震源破裂面的斷層位移和運(yùn)動(dòng)速率等的幾何學(xué)和運(yùn)動(dòng)學(xué)特征，對(duì)深入了解地震重復(fù)間隔、潛在震源區(qū)等問題，從而認(rèn)清地震應(yīng)力應(yīng)變的積累與釋放過程，具有重要意義。地殼形變監(jiān)測(cè)是InSAR技術(shù)應(yīng)用較為成功的領(lǐng)域之一，自Massonnet等(1993)在《nature》上發(fā)表了利用ERS數(shù)據(jù)獲取的美國(guó)Landers地震形變場(chǎng)以來，各國(guó)學(xué)者開展了大量的相關(guān)研究工作，成為研究地震形變場(chǎng)較為流行的一種方法。然而該方法也有其弱點(diǎn):1)對(duì)相位失相關(guān)很敏感，因此對(duì)于震中失相關(guān)嚴(yán)重地區(qū)無法提供形變細(xì)節(jié);2)為獲得形變場(chǎng)必須進(jìn)行相位解纏，需要進(jìn)行大量計(jì)算;3)只能監(jiān)測(cè)到雷達(dá)視線方向上的一維形變量，不能完全反映監(jiān)測(cè)地區(qū)的形變(Van Puymbroeck et al.，2000)。

為了解決上述問題，一些學(xué)者開始探索使用光學(xué)影像偏移量法來作為InSAR技術(shù)的補(bǔ)充(Van Puymbroeck et al.，2000)。人們發(fā)現(xiàn)地表形變還可以通過計(jì)算震前與震后的兩景衛(wèi)星圖像之間的偏移量來獲得，這個(gè)偏移量反映了兩幅圖像上地表指定點(diǎn)的位置偏差。在該思想指導(dǎo)下，近年來利用亞像素級(jí)相位相關(guān)性獲取地表形變逐漸成為研究的熱點(diǎn)。并在一些震例(Massonnet et al.，2000;Michel et al.，2002;Dominguez，2003;Binet，2005)中獲得了較好的結(jié)果，能夠獲取震中附近的形變量。

汶川地震后，一些學(xué)者利用ALOS數(shù)據(jù)獲取了地表InSAR同震形變場(chǎng)(Linlin et al.，2008;孫建寶等，2008;單新建等，2009;屈春燕等，2010)，為理解此次地震的形變過程和發(fā)生機(jī)理，提供了重要的觀測(cè)資料。但由于沿?cái)鄬訋巫兲荻容^大，在形變中心帶出現(xiàn)了1條無干涉條紋的非相干條帶，因此無法獲得斷層帶近場(chǎng)形變量，使得對(duì)了解汶川地震總體同震位錯(cuò)量有一定難度。地震后，國(guó)內(nèi)很多研究團(tuán)隊(duì)對(duì)汶川地震地表破裂帶進(jìn)行了野外實(shí)地考察和測(cè)量，取得了大量的第一手寶貴資料(李海兵等，2008;任俊杰等，2008)，但是人們?cè)谕鹞灰?、地表破裂帶長(zhǎng)度等方面還存在認(rèn)識(shí)上的差別，主要有3方面的原因:不同的測(cè)量方法、不同的觀測(cè)地點(diǎn)以及很多觀測(cè)點(diǎn)無法到達(dá)的限制、對(duì)同一現(xiàn)象的認(rèn)識(shí)不同(李海兵等，2008)。其中，最主要的是由于汶川地震造成大量山體滑坡和崩塌，使道路不通，很多有地表破裂的地區(qū)無法到達(dá)，或被深埋，或破裂帶穿過陡峭山坡無法測(cè)量。而遙感技術(shù)以其覆蓋范圍廣、不受地面條件限制等獨(dú)特優(yōu)勢(shì)，可以成為野外實(shí)地測(cè)量的有益補(bǔ)充。因此，本文擬利用光學(xué)影像相位相關(guān)估算偏移量法，來獲取發(fā)震斷裂帶附近的形變量。

1 原理介紹

利用光學(xué)影像偏移量來獲取斷裂帶附近的形變量需要圖像配準(zhǔn)達(dá)到亞像素級(jí)，這是由于地震形變的量級(jí)一般小于遙感影像的分辨率(少數(shù)大地震可以達(dá)到m級(jí))。并且由于衛(wèi)星的軌道運(yùn)動(dòng)、相機(jī)的掃描運(yùn)動(dòng)及數(shù)字高程模型的誤差等都給精確配準(zhǔn)帶來了困難，需要用嚴(yán)格成像模型來修正原始遙感影像所存在的幾何變形并投影到合適的參考系統(tǒng)中。成像幾何模型反映了地面點(diǎn)3維空間坐標(biāo)與相應(yīng)像點(diǎn)在像平面坐標(biāo)系的2維坐標(biāo)之間的數(shù)學(xué)關(guān)系。另外，地面形變測(cè)量的精度，還依賴于圖像重采樣及相關(guān)性估計(jì)技術(shù)。在圖像配準(zhǔn)前，遙感影像需要投影并重采樣到同一參考系下，一般有2種做法:方法1為選一主圖像為參考基準(zhǔn)，將從影像配準(zhǔn)并投影到主影像的參考系中，雷達(dá)干涉測(cè)量中采用這種方法。方法2為將主、從影像均投影并重采樣到獨(dú)立于衛(wèi)星成像幾何的地面參考系中，之后對(duì)影像的變化檢測(cè)均是基于地面參考系的。本文將采用方法2對(duì)SPOT數(shù)據(jù)進(jìn)行處理來獲得地表形變，整個(gè)過程涉及到SPOT正射糾正、影像重采樣、相關(guān)性估計(jì)等步驟，下面從這3個(gè)方面來進(jìn)行介紹。

1.1 SPOT影像正射糾正

SPOT衛(wèi)星為推掃式成像系統(tǒng)，即光學(xué)成像系統(tǒng)固定，隨著衛(wèi)星平臺(tái)的向前移動(dòng)，逐行以時(shí)序方式獲取沿軌道的連續(xù)影像條帶。SPOT 1A級(jí)產(chǎn)品是提供給用戶的經(jīng)過最低預(yù)處理級(jí)的影像，用戶可在此基礎(chǔ)上根據(jù)需要進(jìn)行嚴(yán)格的集合校正，1A及數(shù)據(jù)各掃描行的像素與相應(yīng)地面行的地面點(diǎn)之間保持嚴(yán)格的中心投影關(guān)系。依據(jù)SPOT提供的高精度衛(wèi)星星歷數(shù)據(jù)、姿態(tài)數(shù)據(jù)及時(shí)間數(shù)據(jù)，建立地面點(diǎn)坐標(biāo)和像點(diǎn)坐標(biāo)之間的幾何關(guān)系進(jìn)行正射糾正，這種幾何關(guān)系可以通過正解法與反解法2種方式建立。

正解法最終的目的為求解像點(diǎn)P(c，r)在地心坐標(biāo)系的視向矢量→u3與高程為h的地球橢球面的交點(diǎn)，期間要經(jīng)過本體坐標(biāo)系到軌道坐標(biāo)系、再由軌道坐標(biāo)系到地心坐標(biāo)系的轉(zhuǎn)換，用到的坐標(biāo)系統(tǒng)如圖1所示(Riazanoff，2002)。

正解法將像平面中規(guī)則的像素網(wǎng)格投影到不規(guī)則的地球表面上，在大范圍尺度上，會(huì)由于地球自轉(zhuǎn)導(dǎo)致畸變，小尺度范圍內(nèi)也會(huì)因?yàn)槌上裣到y(tǒng)姿態(tài)變換及地形起伏而使圖像變形。造成糾正影像上所得之像點(diǎn)非規(guī)則排列，有的像元素內(nèi)可能“空白”(無像點(diǎn))，有的可能重復(fù)(多個(gè)像點(diǎn))，就難以實(shí)現(xiàn)灰度內(nèi)插(張永生，2000)。由于受到圖像分辨率的限制，要想測(cè)量地震形變，圖像的重采樣必須要保持來自原始數(shù)據(jù)中的亞像素級(jí)信息，因此以反解法為宜。

圖1 參考坐標(biāo)系Fig.1 Reference coordinate system.

圖2 反解法示意圖Fig.2 Schematic diagram of indirect image formation.

反解法是從空白的輸出圖像陣列(標(biāo)準(zhǔn)空間)出發(fā)，亦按行列的順序依次對(duì)每個(gè)輸出像元點(diǎn)位反求其原始圖像坐標(biāo)(畸變空間)中的位置。如圖2，首先在地面建立規(guī)則的格網(wǎng)，因?yàn)橐雽?duì)配準(zhǔn)影像進(jìn)行比對(duì)，則需要將它們都校正到相同格網(wǎng)上，這樣2幅影像就具有相同的分辨率，不會(huì)在相關(guān)計(jì)算中出現(xiàn)影像像元錯(cuò)位。

為避免投影后造成像點(diǎn)不規(guī)則問題，Leprince等(2007)提出了一種反變換模型，以滿足嚴(yán)格的重采樣需求。下面對(duì)這一模型進(jìn)行介紹。

地面點(diǎn)的位置僅與成像方向→u3和高程h有關(guān)，若把h確定下來，則成像方向與橢球的交點(diǎn)必定位于高程為h的橢球上，如果能通過像點(diǎn)坐標(biāo)的初始值(c，r)初步確定→u3，然后通過迭代過程逐步修正(c，r)，即可獲得實(shí)際的成像方向。

模型假設(shè)地面任一給定點(diǎn)都能在像平面上找到一個(gè)并且是惟一一個(gè)對(duì)應(yīng)像點(diǎn)。如圖3所示，M為地面一給定點(diǎn)，該點(diǎn)的高程由DEM數(shù)據(jù)經(jīng)3次線性內(nèi)插得到，并且為WGS84坐標(biāo)系統(tǒng)。

考慮視線矢量U3(c，r)對(duì)所有的c，r=1，…，N，拓展到連續(xù)的二維實(shí)數(shù)空間，即u3(x，y)和(x，y)∈R2。那么在像空間中尋找像點(diǎn)坐標(biāo)(x，y)即相當(dāng)于求點(diǎn)(x，y)∈R2使得函數(shù)Ф(x，y)取得最小值:

式(1)中M'(x，y)為從視向矢量→u3所看到的地面點(diǎn)，另為視向矢量的衛(wèi)星位矢，假設(shè)光線在大氣層中為直線傳播，那么的光線延長(zhǎng)矢量，如果 M'位于 s→和地形面的交點(diǎn)，那么由于地形面的非線性特性而使得計(jì)算M'坐標(biāo)變得相當(dāng)繁瑣，為了簡(jiǎn)化問題，為每一個(gè)點(diǎn)M構(gòu)筑一投影平面P(M)，使得M'點(diǎn)實(shí)際位于該平面上，該投影平面P(M)經(jīng)過點(diǎn)M并且垂直于(圖3)。而M∈P(M)，使得方程(1)的解并沒有變。所有位于投影平面P(M)的點(diǎn)M'(α，β，γ)必須滿足，因此該投影平面的方程可以由下式?jīng)Q定:

圖3 反變換模型(據(jù)Leprince et al.，2007)Fig.3 Inverse orthorectification model(after Leprince et al.，2007).

s和投影平面P(M)相交于

其中

則反變換的解可以通過使下式取得極小值得到

其中

這個(gè)過程是一個(gè)反復(fù)迭代的過程，經(jīng)測(cè)試兩點(diǎn)梯度法(TPSS)在穩(wěn)定性及效率方面均比準(zhǔn)牛頓法、最速下降法要好(Barzilai et al.，1988)。

計(jì)算得到的最小值存儲(chǔ)在2個(gè)與地面格網(wǎng)相同維數(shù)大小的矩陣中，如果格網(wǎng)左上角頂點(diǎn)坐標(biāo)為(E0，N0)，格網(wǎng)的分辨率為r，在像平面中的像點(diǎn)［X(i，j)，Y(i，j)］即被投影到地面對(duì)應(yīng)的(E0+i·r，N0－ j·r)點(diǎn)上。

1.2 影像重采樣

由于計(jì)算后的點(diǎn)位置多半不在原圖的像元中心處，即在原始影像中的相應(yīng)位置的坐標(biāo)值不為整數(shù)，因此必須利用在原始影像中該點(diǎn)附近的若干像元的灰度，利用重采樣進(jìn)行賦值。根據(jù)奈奎斯特采樣定律，采樣頻率應(yīng)＞2倍信號(hào)最高頻率，理想的重采樣函數(shù)為辛克(SINC)函數(shù)，在影像處理過程中經(jīng)常用到的最近鄰法、雙線性內(nèi)插及3次卷積等方法，在重采樣過程中會(huì)造成一定程度的失真，會(huì)影響后續(xù)的相關(guān)性處理精度。為此，Leprince等(2007)構(gòu)造了一個(gè)二維重采樣核，具有理想的低通濾波形式:

d x和d y為重采樣距離，代表分別在x和y方向上的鄰采樣的最大間隔。

1.3 相位相關(guān)法估算位錯(cuò)量

相位相關(guān)法的理論依據(jù)是傅立葉變換的相移定理，相移定理是指空間域內(nèi)函數(shù)的位移將會(huì)引起頻域內(nèi)變換函數(shù)的相移，因此一對(duì)配準(zhǔn)圖像之間的相對(duì)位移可以從它們傅里葉變換后的相位差得到。

假設(shè)兩幅影像i1和i2之間只存在位移關(guān)系，平移量為(x0，y0)，即則i1，i2對(duì)應(yīng)的傅里葉變換I1，I2之間的關(guān)系為對(duì)應(yīng)的頻域中2個(gè)圖像的互功率譜為

此函數(shù)在偏移位置處有明顯的尖銳峰值，其他位置的值接近于零，估算峰值處的坐標(biāo)就可找到兩幅圖像間的偏移量。

2 數(shù)據(jù)處理流程

2008年5月12日發(fā)生于龍門山斷裂帶的汶川地震舉世震驚，造成了巨大損失，地表破裂帶長(zhǎng)約240km，汶川地震復(fù)雜的滑動(dòng)分布可能表明幾條活動(dòng)斷層同時(shí)參與了地震破裂過程(Xu et al.，2009)。震后不久葛林林教授便發(fā)表了用InSAR技術(shù)測(cè)量的地表形變初步結(jié)果(Linlin et al.，2008)，體現(xiàn)了遙感技術(shù)在時(shí)效性及覆蓋范圍方面的優(yōu)勢(shì)。然而由于四川地區(qū)多云天氣較多，獲得無云的可見光數(shù)據(jù)比較困難，因此一直未見有利用光學(xué)影像估算地表形變的文章發(fā)表。

本文選取了震前與震后各一景SPOT數(shù)據(jù)作為研究對(duì)象，級(jí)別為L(zhǎng)1A級(jí)。基本參數(shù)見表1。

兩幅影像時(shí)間上相距22個(gè)月，但均在冬季成像，這保證了地物差別不會(huì)很大。選取的1A級(jí)產(chǎn)品是SPOT數(shù)據(jù)經(jīng)輻射校正處理后的產(chǎn)品，包含了用以進(jìn)行后續(xù)的幾何校正處理的輔助數(shù)據(jù)。正射糾正中用到的DEM數(shù)據(jù)為最新發(fā)布的ASTER GDEM version1，地面分辨率為1″(約30m)，在全球范圍內(nèi)滿足垂直精度為20m的置信度為95%(NASA，2009)。下面介紹數(shù)據(jù)處理流程。

(1)選取1/10000的已經(jīng)進(jìn)行正射糾正的SPOT數(shù)據(jù)作為參考影像，對(duì)震前數(shù)據(jù)image1進(jìn)行正射糾正，從影像上選取了12個(gè)GCP點(diǎn)，注意GCP點(diǎn)的選取盡量均勻分布在圖像中間并避開斷裂帶，以避免后續(xù)的衛(wèi)星視向矢量?jī)?yōu)化時(shí)將地震形變量抵消。

表1 SPOT數(shù)據(jù)相關(guān)參數(shù)Table 1 Relevant parameters for the SPOT scenes

(2)從1A級(jí)數(shù)據(jù)的頭文件中提取輔助數(shù)據(jù)，并按照前面介紹的反變換模型對(duì)所選取的GCP點(diǎn)進(jìn)行視向矢量?jī)?yōu)化。

(3)利用優(yōu)化后的GCP點(diǎn)對(duì)image1進(jìn)行正射糾正處理;得到經(jīng)過精確糾正影像i1。

(4)重復(fù)上述步驟對(duì)震后數(shù)據(jù)image2進(jìn)行處理，不同之處在于此時(shí)選取已處理好的i1為參考影像，得到精確糾正影像i2，且保證了i1和i2之間的配準(zhǔn)。

(5)按照2.3節(jié)中的方法對(duì)i1和i2進(jìn)行相位相關(guān)計(jì)算，選取的滑動(dòng)窗口大小為32×32，步長(zhǎng)為8。

需要補(bǔ)充說明的是，根據(jù)衛(wèi)星自帶的軌道數(shù)據(jù)和一些手動(dòng)選取的初始控制點(diǎn)，得到1個(gè)由在參考影像和待配準(zhǔn)影像上選取的一系列同名點(diǎn)對(duì)構(gòu)成的改正矩陣，運(yùn)用間接正射校正模型和改正矩陣，將原始影像即待配準(zhǔn)影像投影至地面坐標(biāo)系，采用相位相關(guān)技術(shù)對(duì)投影后的影像進(jìn)行匹配運(yùn)算，得到二者在EW向及SN向上的坐標(biāo)偏移量及SNR值，根據(jù)該坐標(biāo)差值在DEM的基礎(chǔ)上就構(gòu)成了一次迭代后的新的地面控制點(diǎn);返回步驟2，再將上步得到的新的控制點(diǎn)作為初始控制點(diǎn)數(shù)據(jù)，生成新的改正矩陣。如此經(jīng)過反復(fù)迭代直到方差達(dá)到穩(wěn)定，即可得到精確的地面控制點(diǎn)，然后便可利用這些控制點(diǎn)進(jìn)行正射校正和配準(zhǔn)。如果軌道數(shù)據(jù)的精度已可滿足初始正射校正的需要，則手動(dòng)選取控制點(diǎn)的步驟可以省略。

3 結(jié)果分析

圖4a為正射糾正后的震后影像image2，圖4c，d分別為EW向及SN向的位錯(cuò)分布圖，而圖4b為相位相關(guān)計(jì)算過程中的SNR圖，它代表了相位相關(guān)計(jì)算的質(zhì)量，范圍從0(去相關(guān))到1(相關(guān)性很好)。從該圖上可以看出去相關(guān)分布區(qū)域，去相關(guān)即失去相關(guān)性，表現(xiàn)為SNR較低或無值(Nan)，當(dāng)相關(guān)性算子無法覆蓋區(qū)域或位錯(cuò)量太大超過設(shè)定值(20m)會(huì)出現(xiàn)無值的情況。導(dǎo)致去相關(guān)的原因可歸納為3種:1)云層造成的去相關(guān)，在震后影像image2(Fig.4a)龍門山前緣上空漂浮著一些云朵，它們所在的位置和SNR圖中去相關(guān)區(qū)域相對(duì)應(yīng);2)積雪覆蓋及大量滑坡坍塌處，在龍門山中央海拔較高，覆蓋有積雪(圖像中央)，而在龍門山鎮(zhèn)附近景區(qū)有大量坍塌(圖像左下角)，導(dǎo)致了這些區(qū)域去相關(guān)面積較大;3)河流，在影像的右下方有2條河流經(jīng)過，造成了2條線狀去相關(guān)帶。

圖4 計(jì)算結(jié)果圖Fig.4 Calculation result map.

SPOT4影像的掃描帶寬度為60km，因此選取的這景影像僅覆蓋茂縣、什邡地區(qū)，看不到汶川地震破裂帶的全貌，但在圖幅覆蓋范圍內(nèi)，地震造成的2條破裂帶清晰可見。圖中紅線為根據(jù)圖像不連續(xù)性畫出的斷層位置。在位錯(cuò)分布圖上可以看出，地表破裂帶主要分布于龍門山中央斷裂帶(北川-映秀斷裂)及前山斷裂帶(彭灌斷裂)上，在后山斷裂帶(茂汶斷裂)未見明顯的地表破裂帶。橫跨斷層畫了6條剖面，AA'、DD'、EE'、FF'橫跨龍門山中央斷裂帶。AA'位于龍門山鎮(zhèn)附近，該處水平縮短量最大，達(dá)到10m左右;FF'位于岳家山附近，水平縮短量最大可達(dá)6m;EE'位于清平鄉(xiāng)附近，水平縮短量為2～4m;DD'位于高川鄉(xiāng)附近，水平縮短量又有增大的趨勢(shì)，最大可達(dá)5m，并且在該處斷裂帶有一轉(zhuǎn)向。從EW向位錯(cuò)及剖面圖上可看出，中央斷裂帶兩端位錯(cuò)量大，中間略小，并且有分段的趨勢(shì)，在高川附近有一橫向轉(zhuǎn)折，發(fā)震斷層西北盤為抬升盤，南東盤斷層附近仍然表現(xiàn)為隆起區(qū)，顯示龍門山中央斷裂帶以逆沖為主、兼具走滑的斷層性質(zhì)，中央斷裂帶位錯(cuò)量一般達(dá)到4～6m，而龍門山鎮(zhèn)附近的剖面AA'的最大位錯(cuò)達(dá)到10m左右。在龍門山鎮(zhèn)附近顯示有大片位錯(cuò)量較大的分布區(qū)，調(diào)查表明該地區(qū)為龍門山景區(qū)的回龍溝及銀廠溝所在，地震造成自然景觀損毀率達(dá)80%～90%。銀廠溝的大龍?zhí)?、小龍?zhí)妒苣嗍骱捅浪襟w阻隔而徹底破壞，銀廠溝蓋坪周圍的居民點(diǎn)已被夷為平地，所有房屋都已垮塌(賈建中等，2008)。

圖5 EW向位錯(cuò)及跨斷層剖面圖Fig.5 East-west displacement and profiles across fault.

BB'和CC'橫跨龍門山前山斷裂帶，BB'位于鎣華鎮(zhèn)附近，水平縮短量為2m左右，CC'位于漢旺附近，水平縮短量為2～4m，可見前山斷裂帶逆沖為主，在山前形成1條從八角鎮(zhèn)至睢水鎮(zhèn)長(zhǎng)約40km的破裂帶(因影像未能覆蓋白鹿鎮(zhèn)及以南地區(qū)，因此未能反映白鹿及小魚洞地區(qū)的破裂情況)，位錯(cuò)量一般為1～2m，漢旺鎮(zhèn)地區(qū)可達(dá)4m左右。綜合以上幾個(gè)剖面，可以看出中央斷裂帶位錯(cuò)量大于前山斷裂帶，且兩端位錯(cuò)量大于中間。

震后國(guó)家地震應(yīng)急指揮部立即組織科技人員進(jìn)行了科考，開展了地震地表破裂考察，通過對(duì)地形地貌、地質(zhì)體和地面人工建(構(gòu))筑物的錯(cuò)動(dòng)變形觀測(cè)，完成了地震地表破裂帶填圖和沿地震斷層同震位移的測(cè)量等科學(xué)考察(徐錫偉等，2008)，取得了大量寶貴的現(xiàn)場(chǎng)資料，本文從中選取一些現(xiàn)場(chǎng)調(diào)查點(diǎn)對(duì)SPOT影像獲取結(jié)果進(jìn)行驗(yàn)證。按照文獻(xiàn)中的描述，在龍門山鎮(zhèn)東林寺(31.285 50°N，103.818 30°E)抬升2.5m，右旋位移2.2m，按本文方法在該處附近得到的位錯(cuò)量約為2.7m。在龍門山鎮(zhèn)附近(31.294 194 4°N，103.848 472 2°E)抬升5m，本文方法在該處附近為6.8m。在綿竹市九龍鎮(zhèn)沙壩村(31.981 1°N，104.118 36°E)NE向陡坎上測(cè)量到的最大垂直位移(3.5±0.2)m，本文方法在該處附近為4.6m左右。在漢旺清水河?xùn)|岸通往清平鄉(xiāng)的公路上(31.461 56°N，104.166 39°E)，可見廢棄公路在破裂帶上地殼縮短量為(1.5±0.2)m，本文方法在該處在1.8m左右。雖然通過SPOT影像得到的位錯(cuò)量較野外測(cè)量值偏大20%左右，但是地表位錯(cuò)量的分布趨勢(shì)達(dá)到了很好的吻合。遙感圖像由于噪聲的干擾，某一點(diǎn)的值和周邊的值具有跳躍性，因此最好選取一定大小窗口的平均值作為該點(diǎn)的值。野外測(cè)量的值是3維的，而遙感影像計(jì)算得到的值只反映水平方向的位錯(cuò)量，因此在比較時(shí)我們更注重于整體形態(tài)的分布。圖6為該文獻(xiàn)中野外調(diào)查點(diǎn)與本文計(jì)算結(jié)果的疊合圖，除岳家山至高川段因當(dāng)時(shí)條件惡劣難以到達(dá)而沒有取得調(diào)查點(diǎn)外，圖中調(diào)查點(diǎn)都分布在根據(jù)遙感影像計(jì)算得到的地表破裂帶上。表2為野外調(diào)查的滑動(dòng)位錯(cuò)量與根據(jù)遙感影像計(jì)算得到的水平位錯(cuò)在對(duì)應(yīng)區(qū)域的對(duì)照表。相比而言，野外調(diào)查是離散的各個(gè)觀測(cè)點(diǎn)，而利用遙感影像獲得的形變場(chǎng)則是連續(xù)的，更能反映破裂帶的整體情況，包括整個(gè)破裂帶分布的確切位置、沿走向的形態(tài)變化等。

表2 野外調(diào)查與遙感圖像計(jì)算對(duì)比表Table 2 Comparison between field survey and computation from remote sensing image

4 結(jié)論

本文以覆蓋茂縣地區(qū)的SPOT影像為數(shù)據(jù)源，采用亞像素相位相關(guān)位錯(cuò)法計(jì)算了影像覆蓋地區(qū)的水平位錯(cuò)量，取得了以下認(rèn)識(shí):

(1)在形變較大的近斷層地區(qū)，In-SAR失相關(guān)嚴(yán)重，無法獲得斷層近場(chǎng)形變信息。本文的研究表明，利用光學(xué)影像相關(guān)法能夠獲得近斷層位錯(cuò)量，可以成為InSAR手段的重要補(bǔ)充。但該方法對(duì)影像配準(zhǔn)要求較高，需要加入衛(wèi)星成像幾何、姿態(tài)能輔助數(shù)據(jù)，利用反解法對(duì)SPOT數(shù)據(jù)進(jìn)行了精確正射校正，對(duì)校正后影像再采用頻域相關(guān)法計(jì)算。

圖6 地表破裂與地面調(diào)查點(diǎn)對(duì)比圖Fig.6 Surface rupture inferred from offset map compared with field survey points.

(2)汶川地震造成龍門山斷裂帶上至少2條斷裂同時(shí)發(fā)生破裂，形成了主要地表破裂帶(龍門山鎮(zhèn)-高川破裂帶)和次級(jí)地表破裂帶(漢旺破裂帶)，沿龍門山鎮(zhèn)-高川破裂帶平均位移量為4～6m，從跨斷層剖面可看出，AA'水平縮短量最大，達(dá)到10m左右，向北至FF'、EE'開始減小，在北端DD'又開始增大，并且在高川附近有一轉(zhuǎn)折;漢旺破裂帶基本為純逆沖性質(zhì)，平均位移量一般為1～2m;與野外實(shí)際調(diào)查數(shù)據(jù)相比，位錯(cuò)量值偏大20%左右，但地表位錯(cuò)量的分布趨勢(shì)與野外調(diào)查結(jié)果相吻合。

(3)該方法雖具有一定優(yōu)勢(shì)，能夠監(jiān)測(cè)水平方向運(yùn)動(dòng)，如對(duì)于走滑斷層水平位錯(cuò)量，以及對(duì)于逆沖斷層的水平縮短量都是可以監(jiān)測(cè)的，但對(duì)于垂直方向位錯(cuò)很不敏感，且數(shù)據(jù)來源受到天氣情況的制約，云量大會(huì)造成圖像去相關(guān)而無法獲得滿意的結(jié)果。

本研究實(shí)現(xiàn)了亞像素級(jí)光學(xué)影像相位相關(guān)算法獲取地表形變，并選取了汶川震區(qū)中的茂縣、什邡為實(shí)驗(yàn)區(qū)，計(jì)算結(jié)果得到一些重要的地質(zhì)信息，對(duì)于理解汶川地震的發(fā)震機(jī)理具有重要意義。研究結(jié)果表明，該方法具有較好的應(yīng)用前景，為獲取斷層附近大形變分布提供了可能。