高二濤，蘭艷萍，黃小梅，李 豪，雍 琦

(1.桂林理工大學(xué)測(cè)繪地理信息學(xué)院，桂林 541006；2.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都 611756；3.山西省工業(yè)設(shè)備安裝集團(tuán)有限公司，太原 030000)

2017年8月8日21點(diǎn)19分46秒，四川省阿壩州九寨溝縣發(fā)生7.0級(jí)地震，震中位于距離九寨溝核心景區(qū)西側(cè)約5 km處的比芒村(33°12′N，103°49′E)，震源深度20 km。距離震中200 km的范圍以內(nèi)，近5年期間發(fā)生過(guò)142次3級(jí)以上的地震，此次地震是震級(jí)最大、破壞最嚴(yán)重的一次。

隨著InSAR(interferometric synthetic aperture radar)技術(shù)的深入研究，其在地震形變場(chǎng)提取、火山運(yùn)動(dòng)監(jiān)測(cè)、城市地表沉降和滑坡監(jiān)測(cè)等方面表現(xiàn)出了良好的應(yīng)用前景[1]。洪順英等人利用3種不同視線(line of sight,LOS)向的ENVISAT ASAR數(shù)據(jù)進(jìn)行干涉處理，獲取多視線向(multi-LOS)的同震形變場(chǎng)[2]。2018年，曹海坤等[3]利用InSAR觀測(cè)技術(shù)融合GPS(global positioning system)數(shù)據(jù)，通過(guò)添加一個(gè)跟點(diǎn)位有關(guān)的系統(tǒng)誤差函數(shù)，對(duì)系統(tǒng)誤差進(jìn)行修正，得到更為準(zhǔn)確的三維形變場(chǎng)。在地震形變場(chǎng)觀測(cè)中，D-InSAR(differential synthetic aperture radar interferometric)技術(shù)具有不受時(shí)相和氣象條件限制、可以全天時(shí)全天候的探測(cè)地球表面細(xì)微的變化[4]、覆蓋范圍廣、觀測(cè)精度高等獨(dú)特優(yōu)勢(shì)。

由于地震發(fā)生在九寨溝景區(qū)內(nèi)，植被、森林密布給GPS、水準(zhǔn)測(cè)量等常規(guī)方法帶來(lái)了較大的困難，InSAR技術(shù)無(wú)需到實(shí)地測(cè)量，并且可以獲取大范圍監(jiān)測(cè)信息，具有較大的優(yōu)勢(shì)。基于歐空局Sentienl-1A衛(wèi)星獲取的九寨溝震前和震后的升降軌雷達(dá)影像，結(jié)合精度較高的AW3D30(ALOS World 3D-30 m)數(shù)字高程模型，利用D-InSAR技術(shù)分別對(duì)升降影像進(jìn)行數(shù)據(jù)處理，為了獲取九寨溝地震的精確地表同震形變場(chǎng)，進(jìn)行了升降軌結(jié)果的交叉驗(yàn)證分析。結(jié)合MATLAB和Surfer繪圖軟件實(shí)現(xiàn)對(duì)形變場(chǎng)的精度分析，最后導(dǎo)入Google Earth驗(yàn)證形變場(chǎng)位置的準(zhǔn)確性。

1 差分干涉測(cè)量基本原理

差分干涉測(cè)量基本原理是基于地震前后SAR影像獲得地表的形變信息。即利用兩幅跨越形變周期的SAR圖像對(duì)，共軛相乘獲取干涉相位信息[5]。通過(guò)衛(wèi)星觀測(cè)幾何參數(shù)和軌道信息去除參考橢球相位，使用外部DEM(digital elevation model)去除地形相位，扣除大氣和噪聲相位，最后得到地表形變相位信息。

合成孔徑雷達(dá)差分干涉測(cè)量工作模式通常有兩種，即通過(guò)兩個(gè)傳感器同時(shí)對(duì)同一區(qū)域進(jìn)行觀測(cè)的單軌模式和單一傳感器重復(fù)對(duì)同一區(qū)域進(jìn)行觀測(cè)的重復(fù)軌道模式[6](圖1)。

圖1 差分測(cè)量原理

S1為衛(wèi)星第一次成像的位置，S2為衛(wèi)星第二次成像的位置，R1和R2分別對(duì)應(yīng)兩次成像到地面目標(biāo)點(diǎn)的視線距離，令R2=R1+Δr。B是衛(wèi)星兩次成像的空間基線，α是基線與水平面的夾角，S1位置時(shí)衛(wèi)星的側(cè)視角為θ，S1與參考平面的高度為H，h為目標(biāo)點(diǎn)P點(diǎn)與參考平面之間的高度。干涉相位可以表示為

Φ=Φflat+Φtopo+Φdef+Φatm+Φnoise

(1)

式(1)中：Φflat為平地相位；Φtopo為地形相位；Φdef為要獲取的地表形變相位；Φatm、Φnoise分別為大氣和噪聲相位。在去除大氣和噪聲相位后，式(1)可表達(dá)為

(2)

(3)

式中：-Δr表示衛(wèi)星兩次對(duì)同一目標(biāo)成像的路程差；λ表示雷達(dá)波長(zhǎng)。其中，Φflat可以用傳感器本身的頭文件等相關(guān)參數(shù)計(jì)算后消除，Φtopo可以用外部DEM模擬出地形相位，并從中扣除。這樣，地表形變相位可以通過(guò)式(3)獲得。

2 九寨溝地質(zhì)背景構(gòu)造

根據(jù)九寨溝地震震中位置(圖2)，結(jié)合已有的地質(zhì)資料可知，此次地震位于巴顏喀拉塊體東緣岷江斷裂、塔藏?cái)嗔押突⒀罃嗔阎車前皖伩瓑K體發(fā)生激烈碰撞的結(jié)果。東昆侖斷裂向東部延伸形成塔藏?cái)嗔眩诎皖伩w塊東北方向上，大體上呈現(xiàn)出西北走向的趨勢(shì)，它涵蓋了東北村、塔藏、九寨溝、沙尕里等廣泛地區(qū)[7]。對(duì)塔藏?cái)嗔盐鞫蔚难芯?，主要以左旋剪切走滑為主兼擠壓活動(dòng)為主，滑動(dòng)速率為2.43～2.89 mm/a[8-9]。

圖2 地震地質(zhì)背景

塔藏塊體東北部在青藏高原的東北方向上，是柴達(dá)木地塊、揚(yáng)子地塊、巴顏喀拉地塊的交匯地區(qū)，是我國(guó)南北地震帶地震中段。此處地形極其復(fù)雜，板塊之間的運(yùn)動(dòng)劇烈，發(fā)生過(guò)許多強(qiáng)震。岷江斷裂北起九寨溝縣，南至松潘縣，斷層走向主要以左旋為主，傾向?yàn)槲鞅狈较?，傾角在60°～70°之間，顯示出該斷裂強(qiáng)烈活動(dòng)的現(xiàn)象。

虎牙斷裂是這次地震的發(fā)震斷裂，呈現(xiàn)出由西北向南北方向的延伸，為上沖兼左旋走滑方式[10]。斷裂的走向方向表現(xiàn)出上部和下部不同的現(xiàn)象，在上部，發(fā)生從北西方向向南北方向過(guò)渡，向東傾斜80°左右；在下部由南北方向向南東方向過(guò)渡，傾向向西，傾角在30°～70°范圍變化。在這樣的情況和地勢(shì)下，使得巴顏喀拉塊體向東部運(yùn)動(dòng)激烈，導(dǎo)致九寨溝地震的發(fā)生。

3 九寨溝地震形變場(chǎng)提取

3.1 實(shí)驗(yàn)數(shù)據(jù)

利用Sentinel-1A衛(wèi)星獲取LOS方向形變場(chǎng)。其中使用的4幅影像參數(shù)如表1所示。

表1 Sentinel-1A升降軌影像參數(shù)

DEM數(shù)據(jù)使用的是AW3D30。它是日本宇航局對(duì)地觀測(cè)中心于2016年5月發(fā)布的全球數(shù)字高程模型，它是由ALOS(advanced land observation satellite)衛(wèi)星上搭載的全色立體遙感測(cè)繪儀PRISM(panchromatic remote-sensing instrument for stereo mapping)生成，覆蓋了北緯82°到南緯82°之間的全球陸地面積，空間分辨率為30 m，高程精度為5 m。AW3D30是當(dāng)今用戶獲取的覆蓋范圍最廣，精度最高的開(kāi)源數(shù)字高程模型之一[11]。通過(guò)MATLAB軟件對(duì)DEM數(shù)據(jù)進(jìn)行拼接、裁剪得到研究區(qū)DEM，運(yùn)用其去地形相位處理。

3.2 數(shù)據(jù)處理過(guò)程

主要使用ENVI軟件中的SARscape模塊，結(jié)合AW3D30 DEM數(shù)據(jù)分別對(duì)Sentinel-1A升降軌影像對(duì)進(jìn)行差分干涉處理。技術(shù)路線如圖3所示。

圖3 D-InSAR技術(shù)路線

3.2.1 影像配準(zhǔn)

首先找出兩景影像的同名點(diǎn)，計(jì)算它們的坐標(biāo)函數(shù)關(guān)系，按照函數(shù)關(guān)系把從影像經(jīng)過(guò)簡(jiǎn)單的平移，使兩幅影像有統(tǒng)一的分辨單元。配準(zhǔn)的步驟分為粗配準(zhǔn)和精配準(zhǔn)[6]，得到配準(zhǔn)信息如表2所示。

表2 主從影像配準(zhǔn)信息

從表2可以看出，圖像配準(zhǔn)精度較好，距離向、方位向、信噪比和相關(guān)系數(shù)都比較穩(wěn)定。

3.2.2 生成干涉圖

對(duì)配準(zhǔn)后的主影像和從影像進(jìn)行共軛乘法運(yùn)算來(lái)生成干涉圖，設(shè)置距離向和方位向的多視比例為19∶5。評(píng)估所生成的干涉圖的質(zhì)量，通常用的方法是使用相干系數(shù)對(duì)干涉圖進(jìn)行評(píng)價(jià)[12]。相關(guān)系數(shù)的范圍為[0,1]，一般相干系數(shù)大于0.3，說(shuō)明干涉圖質(zhì)量較好。

對(duì)于平地相位使用傳感器本身的頭文件攜帶的軌道參數(shù)計(jì)算并減去；利用AW3D30 DEM模擬地形相位信息，并從干涉結(jié)果中去除。經(jīng)過(guò)以上兩步處理后，利用MATLAB對(duì)干涉圖進(jìn)行裁剪，去平、去地形相位后的干涉圖如圖4所示。

圖4 扣除地形相位后的干涉圖

從圖4可以看出地震的干涉圖條紋不是很清晰，因?yàn)榇嬖谠肼?、近?chǎng)失相關(guān)等的影響，需對(duì)其進(jìn)行濾波處理，如圖5所示為自適應(yīng)濾波并編碼的結(jié)果，通過(guò)濾波處理之后，干涉條紋信息明顯清晰了很多。

圖5 D-InSAR差分干涉圖濾波

3.2.3 相位解纏

相位解纏的本質(zhì)是恢復(fù)整周相位周期的過(guò)程。九寨溝景區(qū)，植被密度較高，使用最小費(fèi)用流(minimum cost flow, MCF)的方法進(jìn)行解纏可以提高解纏效率。這種方法對(duì)像元進(jìn)行了充分的考慮，對(duì)于一些相干性系數(shù)小的像元，進(jìn)行了掩膜處理，提高了解纏的精度。如圖6所示，利用MATLAB裁剪出的相位解纏結(jié)果。

圖6 DInSAR相位解纏

3.2.4 軌道精煉和重去平

軌道精煉是通過(guò)攝取控制點(diǎn)來(lái)校正細(xì)小的軌道誤差，從而提高精度。在解纏和干涉圖平坦的地方選點(diǎn)，震中區(qū)域不可以選點(diǎn)，由于震中失相關(guān)，會(huì)給重去平造成很大的誤差。

3.2.5 相位轉(zhuǎn)形變和地理編碼

把軌道精煉后的結(jié)果轉(zhuǎn)換成形變量，它的本質(zhì)是把雷達(dá)的形變相位轉(zhuǎn)換成一個(gè)形變量(以米為單位)，即乘以λ/(4π)以獲得LOS方向形變場(chǎng)，并與外部DEM結(jié)合通過(guò)地理編碼投影到WGS84(World Geodetic System 1984)坐標(biāo)系。

利用MATLAB軟件裁剪形變場(chǎng)，使用freadbk命令讀取結(jié)果并裁剪出所需要的形變場(chǎng)，并且通過(guò)格式轉(zhuǎn)換，導(dǎo)入Surfer軟件，如圖7所示，實(shí)驗(yàn)獲取的升降軌形變場(chǎng)。

圖7 D-InSAR獲取的地表形變場(chǎng)

為了驗(yàn)證形變場(chǎng)位置的準(zhǔn)確性，在Surfer軟件中，把形變場(chǎng)轉(zhuǎn)化為Google Earth能識(shí)別的KML(keyhole markup language)文件，將實(shí)驗(yàn)得到的升降軌形變場(chǎng)疊加顯示在Google Earth地圖上。從圖8中可以看出形變出現(xiàn)的區(qū)域與實(shí)際地震發(fā)生的區(qū)域吻合度較高。

4 形變結(jié)果及分析

基于Sentinel-1A影像，利用InSAR技術(shù)，提取九寨溝地震升降軌形變場(chǎng)，對(duì)形變場(chǎng)的范圍、時(shí)空分布以及形變大小進(jìn)行分析，為地震斷層反演打下良好基礎(chǔ)。

首先，從圖7可以看出，形變場(chǎng)的分布不對(duì)稱，左邊(西邊)的形變程度較嚴(yán)重，右邊(東邊)形變相對(duì)較弱。其主要原因是該地發(fā)震斷層處于塔藏?cái)鄬印⒒⒀罃鄬雍歪航瓟鄬舆@三個(gè)斷層交匯的地震帶上，斷層結(jié)構(gòu)相對(duì)復(fù)雜所致。斷層活動(dòng)表現(xiàn)出不確定特性，結(jié)合歷史資料，該地震帶上發(fā)生過(guò)很多次大地震，顯示出地震群體的特征。斷層活動(dòng)的變化發(fā)生的方式不同，變形場(chǎng)的分布特征也遵循不同的模式。2003年巴姆地震是一次右旋走滑形式，它在LOS形變上呈現(xiàn)出正四面體分布[13]；2008年汶川地震也是一次右旋走滑形式，它在LOS形變分布不明顯[6]。九寨溝地形復(fù)雜，在提取地震形變場(chǎng)時(shí)，需要考慮斷裂山脈的滑坡影響[14]。九寨溝地震形變場(chǎng)與巴姆地震形變場(chǎng)有類似之處，都處于四象限分布，不同之處在于九寨溝地震形變場(chǎng)是不對(duì)稱分布，而巴姆地震形變場(chǎng)是對(duì)稱分布。

其次，從圖5可以看出，震中區(qū)域的相干性較低。其主要原因有兩個(gè)，一是因?yàn)榫耪瘻暇皡^(qū)植被森林密集和地形復(fù)雜所致。植被森林密集在雷達(dá)反射過(guò)程中，由于綠色植被中水分含量較多，吸收率變大，因此反射率就大大下降了。雷達(dá)波接收不到反射光信號(hào)，在干涉處理過(guò)程中就出現(xiàn)了干涉失相關(guān)的現(xiàn)象。此外，震中造成的破壞力度強(qiáng)，導(dǎo)致干涉失相關(guān)。西北和東南方向的干涉條紋顯示出不同形狀。在西北方向干涉條紋表現(xiàn)出彩虹漸變色變化趨勢(shì)；而在東南方向上變化較弱。從另一個(gè)層面上說(shuō)明地形、地物復(fù)雜，變形對(duì)該地的破壞程度不同。

最后，從整體上看，它的總體形狀是一個(gè)“核仁”狀，長(zhǎng)度大約有24 km，寬度大約有29 km，并且呈現(xiàn)出西南向東北走向。從MATLAB中精確的讀取到九寨溝地震LOS向形變的相關(guān)參數(shù)，形變場(chǎng)的面積大約為696 km2，上升的最大形變值達(dá)到22.4 cm，約有8個(gè)干涉條紋；下降極值約為11.2 cm，約有4個(gè)干涉條紋。說(shuō)明該地區(qū)變形較大，地震破壞程度較大。

5 結(jié)論

利用升降軌Sentinel-1A數(shù)據(jù)，基于D-InSAR技術(shù)進(jìn)行SAR數(shù)據(jù)處理，通過(guò)圖像配準(zhǔn)等一系列數(shù)據(jù)處理，結(jié)合已有的資料和谷歌地球等數(shù)據(jù)，利用MATLAB、Surfer軟件實(shí)現(xiàn)對(duì)形變場(chǎng)的裁剪和顯示，得到了2017年九寨溝7.0級(jí)地震精確同震形變場(chǎng)。最后本文對(duì)形變場(chǎng)的位置精度、形變量級(jí)大小等進(jìn)行了系統(tǒng)的分析，為下一步地震斷層位置反演打好了良好的前期基礎(chǔ)。

致謝:感謝歐空太空局在科學(xué)數(shù)據(jù)中心網(wǎng)站公布的Sentinel-1A衛(wèi)星SAR影像；感謝日本宇宙航空航天局地球觀測(cè)研究中心提供的AW3D30 DEM數(shù)據(jù)。