利用升降軌Sentinel-1A影像提取2017年九寨溝Ms 7.0地震同震形變場

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

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

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

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

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

1 差分干涉測量基本原理

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

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

圖1 差分測量原理

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

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

(1)

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

(2)

(3)

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

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

根據(jù)九寨溝地震震中位置(圖2)，結(jié)合已有的地質(zhì)資料可知，此次地震位于巴顏喀拉塊體東緣岷江斷裂、塔藏斷裂和虎牙斷裂周圍，是巴顏喀拉塊體發(fā)生激烈碰撞的結(jié)果。東昆侖斷裂向東部延伸形成塔藏斷裂，在巴顏喀拉體塊東北方向上，大體上呈現(xiàn)出西北走向的趨勢，它涵蓋了東北村、塔藏、九寨溝、沙尕里等廣泛地區(qū)[7]。對塔藏斷裂西段的研究，主要以左旋剪切走滑為主兼擠壓活動為主，滑動速率為2.43～2.89 mm/a[8-9]。

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

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

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

3 九寨溝地震形變場提取

3.1 實驗數(shù)據(jù)

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

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

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

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

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

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

3.2.1 影像配準(zhǔn)

首先找出兩景影像的同名點，計算它們的坐標(biāo)函數(shù)關(guān)系，按照函數(shù)關(guān)系把從影像經(jīng)過簡單的平移，使兩幅影像有統(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 生成干涉圖

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

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

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

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

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

3.2.3 相位解纏

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

圖6 DInSAR相位解纏

3.2.4 軌道精煉和重去平

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

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

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

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

圖7 D-InSAR獲取的地表形變場

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

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

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

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

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

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

5 結(jié)論

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

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