融合DInSAR和MAI技術(shù)反演九寨溝地震三維形變場(chǎng)

宮熙雯,張正加,陳啟浩,劉修國(guó)

(中國(guó)地質(zhì)大學(xué)(武漢)，湖北 武漢 430070)

2017年8月8日21:19:46，四川省北部阿壩州九寨溝縣發(fā)生7.0級(jí)地震(33.2°N，103.82°E)，震源深度20 km，截至2017年8月18日09:00共記錄到地震總數(shù)為5 020個(gè)(余震總數(shù)5 019個(gè))，其中最大余震為：2017年8月9日10:17四川阿壩州九寨溝縣4.8級(jí)。此次地震，中國(guó)地震局徐錫偉等專家判斷這次地震與巴顏喀拉塊體南東向推擠有關(guān)系，發(fā)震構(gòu)造推測(cè)為塔藏?cái)嗔涯蟼?cè)分支和虎牙斷裂北段[1]。

一些學(xué)者對(duì)該地震進(jìn)行了相關(guān)研究，鄭緒君等采用迭代反褶積和疊加方法反演得到破裂模型，用于確定地震的破裂過(guò)程[2]，汪建軍等采用庫(kù)侖應(yīng)力模型研究地震激發(fā)的庫(kù)侖應(yīng)力變化[3]，文獻(xiàn)[4]采用多步定位法對(duì)九寨溝地震Ms7.0的主震及其余震進(jìn)行重新定位，提出震源機(jī)制解決方案，并通過(guò)CAP法(Cut and Paste Method)反演了九寨溝地震的質(zhì)心深度。除了對(duì)以上地震參數(shù)的反演，獲取地震同震形變場(chǎng)對(duì)揭示地表形變、理解和分析發(fā)震構(gòu)造具有重要意義。由于差分雷達(dá)干涉測(cè)量(DInSAR)可以獲取全天時(shí)全天候、大范圍、連續(xù)空間覆蓋、高精度的地震區(qū)域同震形變圖[5]，在地震監(jiān)測(cè)方面應(yīng)用廣泛，對(duì)于本次九寨溝Ms7.0地震，趙強(qiáng)等利用Sentinel-1A和RadarSat-2升軌數(shù)據(jù)結(jié)合兩軌法DInSAR反演了地震的同震形變場(chǎng)[6]；熊威等使用Sentinel-1A數(shù)據(jù)進(jìn)行同震形變場(chǎng)反演，并確定形變區(qū)的相干系數(shù)[7]；申文豪等采用DInSAR方法利用升、降軌Sentinel-1數(shù)據(jù)反演了同震形變場(chǎng)[8]。Zhao基于InSAR和GPS測(cè)量，反演地震同震滑動(dòng)分布和傾角，并獲得地震視線向的同震形變場(chǎng)[9]。但由于DInSAR方法的設(shè)計(jì)本質(zhì)是一種側(cè)視雷達(dá)測(cè)量技術(shù)，得到的并不是地表真實(shí)形變，而是地表真實(shí)形變矢量通過(guò)一系列坐標(biāo)系轉(zhuǎn)換和幾何轉(zhuǎn)換后在視線向的投影[10]。因此，同震形變的反演僅僅依靠DInSAR技術(shù)獲取視線向形變是不夠的，需要增加觀測(cè)角度以建立地震三維同震形變場(chǎng)。

1999年，Michael 等人提出了一種基于方位向像素偏移(Azimuth Pixel Offset，AZO)的技術(shù)[11]，2002年，F(xiàn)ialko等結(jié)合DInSAR和像元偏移量技術(shù)(Offset Tracking)成功獲取Hector Mine地震的三維位移場(chǎng)[12]。2004年，Wright 等提出通過(guò)不同入射角影像進(jìn)行干涉的多視角InSAR技術(shù)并反演Nenana Mountain地震的三維形變場(chǎng)[13]，王家慶等在此基礎(chǔ)上結(jié)合先驗(yàn)條件的最小二乘迭代逼近法解算改后地震的三維同震形變場(chǎng)[14]，但由于多視角InSAR對(duì)數(shù)據(jù)入射角差異要求較高，獲取滿足條件的數(shù)據(jù)源對(duì)絕大多數(shù)地區(qū)較難實(shí)現(xiàn)。2006年，Bechor等提出了多孔徑干涉技術(shù)(Multiple Aperture InSAR，MAI)以獲取方位向形變，較AZO技術(shù)獲取方位向位移矢量精度更高[15]。2014年，Wang等使用該方法并結(jié)合DInSAR方法，使用L波段PALSAR數(shù)據(jù)和C波段ASAR數(shù)據(jù)反演了玉樹(shù)地震三維形變場(chǎng)，實(shí)驗(yàn)結(jié)果表明，使用DInSAR和MAI方法可有效提高震源參數(shù)的反演精度[16]。

文中在九寨溝Ms7.0地震同震形變場(chǎng)反演過(guò)程中，融合升降軌Sentinel-1A數(shù)據(jù)，使用DInSAR和MAI組合方法獲取雷達(dá)視線向和方位向的形變，構(gòu)建地震三維同震形變場(chǎng)，獲取本次地震在豎直方向、南北方向、東西方向的地表位移三分量，較單一雷達(dá)視線向更充分地解釋了地表真實(shí)形變情況。

1 實(shí)驗(yàn)區(qū)與數(shù)據(jù)源

1.1 實(shí)驗(yàn)區(qū)概況

在四川省的北部，阿壩州九寨溝縣于2017年8月8日發(fā)生了一次Mw7.0級(jí)地震，震源中心為(33.20°N，103.82°E)。地震發(fā)生在青藏高原與四川盆地兩大地貌單元過(guò)渡的深切割高山峽谷地帶，地處岷山山脈北段的復(fù)背斜上，地勢(shì)西北高，東南低，平均海拔在4 000 m以上，受災(zāi)面積約為840 km2(35 km×24 km)。研究區(qū)內(nèi)地形復(fù)雜，有多條斷裂帶縱橫交錯(cuò)，震源中心位于塔藏?cái)嗔涯蟼?cè)，虎牙斷裂延長(zhǎng)線北部，與岷江斷裂相交[1,17]，如圖1所示。

圖1 研究區(qū)概況

1.2 數(shù)據(jù)選取

哨兵1號(hào)衛(wèi)星是歐洲航天局哥白尼計(jì)劃(GMES)中的地球觀測(cè)衛(wèi)星，為雙星系統(tǒng)，搭載C波段合成孔徑雷達(dá)，可全天時(shí)、全天候不受天氣條件影響地提供連續(xù)影像。本研究使用數(shù)據(jù)為Sentinel-1A衛(wèi)星干涉寬幅模式的斜距單視復(fù)數(shù)產(chǎn)品(IW SLC)，空間分辨率為5 m×20 m，極化方式為VV。選取了震前震后的升降軌SAR數(shù)據(jù)共4景用于本次地震同震形變反演，影像基本信息如表1所示。本研究所采用的參考DEM為覆蓋震區(qū)的SRTM-3 version4 DEM數(shù)據(jù)。本次實(shí)驗(yàn)選取的升降軌震前震后兩個(gè)影像對(duì)的時(shí)間基線相近均為12 d，空間基線也相對(duì)較小，分別為36.7 m和65.4 m，盡可能減少了時(shí)空失相關(guān)對(duì)最終結(jié)果造成的誤差。

表1 Sentinel-1A影像參數(shù)信息

2 研究方法

文中使用DInSAR和MAI結(jié)合的方法獲取地表位移分量。由于MAI所獲取方位向位移的精度在方位向形變獲取中精度較高，且和DInSAR方法在一個(gè)數(shù)量級(jí)[18]，因此可結(jié)合上述數(shù)據(jù)處理得到的形變圖數(shù)據(jù)，對(duì)九寨溝地區(qū)的地表形變?nèi)至窟M(jìn)行解算，方法流程圖如圖2所示。

圖2 方法流程圖

2.1 DInSAR

DInSAR技術(shù)是通過(guò)對(duì)同一區(qū)域不同時(shí)相兩景影像產(chǎn)生的干涉圖差分處理，借助外部DEM消除地形影響而獲得形變信息一種技術(shù)手段[19]。DInSAR獲取視線向形變的主要步驟為：選取主、輔影像生成干涉圖，并借助外部DEM進(jìn)行自動(dòng)配準(zhǔn)并進(jìn)行去平地效應(yīng)；采用Goldstein方法進(jìn)行前置濾波以提高干涉條紋的清晰度；采用Delaunay MCF方法進(jìn)行解纏以消除低相干像元的影響，只對(duì)高相干像元進(jìn)行解纏，最小化相位突變的影響，減少“孤島”情況的存在；對(duì)于山區(qū)可能存在的大氣延遲，主要表現(xiàn)為與地形相關(guān)的相位斜坡，用線性模型進(jìn)行校正[20]，對(duì)于軌道誤差引起的相位斜坡，用雙二次模型進(jìn)行校正[21]，后經(jīng)軌道精煉、重去平及地理編碼等得到視線向形變。

2.2 MAI

MAI技術(shù)是基于分光束InSAR進(jìn)行數(shù)據(jù)處理，生成前向和后向干涉圖，根據(jù)兩幅干涉圖上的相位差獲取軌道方位向的位移分量的一種干涉測(cè)量方法[22]。MAI技術(shù)獲取方位向形變的主要步驟為：選取主、輔影像生成干涉圖，主、輔影像分別通過(guò)帶通濾波，得到各自的前視和后視的單視復(fù)數(shù)影像，主輔影像的前視和后視數(shù)據(jù)對(duì)應(yīng)進(jìn)行干涉處理，分別得到前視干涉圖和后視干涉圖，并將二者進(jìn)行干涉處理得到MAI干涉圖，通過(guò)Goldstein濾波、區(qū)域增長(zhǎng)法解纏，對(duì)干涉圖中存在的軌道誤差和大氣延遲進(jìn)行校正后，經(jīng)軌道精煉、重去平及地理編碼后得到方位向形變。

2.3 三維形變場(chǎng)

文中使用的升軌數(shù)據(jù)的軌道運(yùn)行方向?yàn)闁|南—西北，降軌數(shù)據(jù)的軌道運(yùn)行方向?yàn)闁|北—西南，由于Sentinel-1衛(wèi)星為右視雷達(dá)，升降軌數(shù)據(jù)與真實(shí)地表形變?nèi)S分量的幾何關(guān)系分別如圖3、圖4所示。

圖3 升軌與地表三維向幾何關(guān)系

圖4 降軌與地表三維向幾何關(guān)系

圖中dU，dN，dE分別為豎直方向、南北方向、東西方向上的形變，規(guī)定向上、正北、正東為正；dlos，daz分別為視線向、方位向上的形變，規(guī)定沿衛(wèi)星軌道方位向?yàn)檎?，統(tǒng)一正向可得到地震三維形變的觀測(cè)方程為：

(1)

式中:θ為雷達(dá)波入射角，α為軌道方位角，下標(biāo)A，D分別表示升軌和降軌。觀測(cè)矢量D:

D=(dlos,A,dlos,D,daz,A,daz,D)T.

(2)

得到誤差方程：

V=CX-D.

(3)

其中，向量V為每次觀測(cè)的誤差，即：

V=(Vlos,A,Vlos,D,Vaz,A,Vaz,D)T.

(4)

由加權(quán)最小二乘原理[23]可計(jì)算三維位移場(chǎng)：

X=(ATPA)-1ATPD.

(5)

式中:C為觀測(cè)方程的系數(shù)矩陣，P為觀測(cè)權(quán)陣，由于DInSAR測(cè)量誤差和MAI的方位向形變誤差都在cm級(jí)，近似認(rèn)為4次觀測(cè)的精度是相等的，即權(quán)陣P為單位矩陣。

3 結(jié)果分析

3.1 InSAR處理結(jié)果

融合升降軌SAR數(shù)據(jù)使用DInSAR和MAI技術(shù)反演的地震升降軌視線向、方位向同震形變?nèi)鐖D5所示，同震位移剖面如圖6所示。

綜合分析圖5、圖6可知，本次九寨溝Ms7.0地震，形變主要在西北區(qū)域和東南區(qū)域，存在非對(duì)稱性分布特征。升降軌數(shù)據(jù)的形變特征顯示出不一致性：升軌數(shù)據(jù)形變反演結(jié)果顯示，震中西北側(cè)區(qū)域發(fā)生了沿視線向最大值為0.19 m、沿方位向最大值為0.25 m的沉降形變，震中東南側(cè)區(qū)域發(fā)生了沿視線向最大值為0.10 m、沿方位向最大值為0.11 m的隆升形變；而在降軌數(shù)據(jù)形變反演結(jié)果中，震中西北側(cè)區(qū)域發(fā)生了沿視線向最大值為0.14 m、沿方位向最大值為0.14 m的隆升形變，震中東南側(cè)區(qū)域發(fā)生了沿視線向最大值為0.10 m、沿方位向最大值為0.14 m的沉降形變。這種不一致性主要與升降軌不同的成像幾何特征以及地表形變方向有關(guān)，而且升降軌呈現(xiàn)出相反的形變態(tài)勢(shì)，表明九寨溝Ms7.0地震的地表形變以水平形變?yōu)橹?，與地震學(xué)結(jié)論保持一致[24]。

由方位向升降軌同震形變圖可知，九寨溝Ms7.0地震不僅引發(fā)了地表豎直方向上的位移，同時(shí)在地表平面上也存在位移且形變量和視線向的形變量為同一個(gè)數(shù)量級(jí)，在以水平形變?yōu)橹鞯木耪瘻系卣鹬?，?gòu)建三維同震形變場(chǎng)對(duì)于其發(fā)震機(jī)理及構(gòu)造分析有著重要的意義。

圖5 九寨溝Ms7.0地震同震形變場(chǎng)

圖6 九寨溝Ms7.0地震同震形變剖面

3.2 三維形變場(chǎng)

將升軌雷達(dá)波入射角θA、降軌雷達(dá)波入射角θD、升軌軌道方位角αA、降軌軌道方位角αD帶入觀測(cè)方程的系數(shù)矩陣C：

(6)

進(jìn)而求得三維形變場(chǎng)X的系數(shù)矩陣Z：

(7)

從而得到地表的三維同震形變分量，即：

dU= 0.645 8IlosA+0.646 3IlosD+

0.0934IAzA-0.093 5IAzD.

(8)

dN=-0.000 1IlosA+0.000 1IlosD+

0.512 9IAzA-0.512 9IAzD.

(9)

dE=-0.716 9IlosA+ 0.716 5IlosD-

0.258 6IAzA-0.258 4IAzD.

(10)

式中：IlosA，IlosD，IAzA，IAzD分別表示升軌視線向形變結(jié)果、降軌視線向形變結(jié)果、升軌方位向形變結(jié)果和降軌方位向形變結(jié)果。根據(jù)本次地震的影響范圍和升降軌數(shù)據(jù)的重疊情況對(duì)形變結(jié)果進(jìn)行了裁剪，選取103.49°E—104.11°E、32.93°N—33.55°N的震中區(qū)域，并借助Matlab編程計(jì)算求解九寨溝地震三維同震形變場(chǎng)，得到地表位移三分量，地震三維同震形變圖如圖7所示。實(shí)驗(yàn)結(jié)果表明，本次地震的斷裂帶為西北—東南(NW-SE)走向，位置大致與地震實(shí)際破裂帶相符[25]。在豎直方向的同震形變量較小且震源中心附近上下形變貢獻(xiàn)較為平均，存在震源中心東南側(cè)最大向上0.09 m、西北側(cè)最大向下0.11 m的形變；在南北方向存在較為明顯的位移分界，存在震源中心東南側(cè)最大向北0.09 m、西北側(cè)最大向南0.18 m的形變；在東西方向的同震形變量較大，且存在明顯的位移分界，存在震源中心西北側(cè)最大向東0.25 m、東南側(cè)最大向西0.12 m的形變。綜合地表位移三分量可以判斷，地表在震源中心的西北側(cè)表現(xiàn)為沉降運(yùn)動(dòng)，在震源中心的東南側(cè)表現(xiàn)為隆升運(yùn)動(dòng)，但地表運(yùn)動(dòng)縱向起伏較小，地表平面地表位移運(yùn)動(dòng)較為明顯，與李渝生等在震后實(shí)地考察測(cè)量結(jié)果相符[17]。

圖7 九寨溝Ms7.0地震三維同震形變場(chǎng)

九寨溝地震的震源中心被塔藏?cái)嗔选⑨航瓟嗔押突⒀罃嗔训难娱L(zhǎng)線3條較大的斷裂帶包圍，圖7中黑色虛線位置為本次地震中地表出露的位置[25]，即斷層所在地表部分的破裂區(qū)域，該斷層兩側(cè)在圖中的位移在3個(gè)方向均較為明顯。地震斷裂面西北側(cè)(下降盤(pán))錯(cuò)動(dòng)位移較大，主要向東南偏東方向滑動(dòng)；東南側(cè)(上升盤(pán))較西北側(cè)錯(cuò)動(dòng)位移小且3個(gè)方向位移的貢獻(xiàn)較為平均，大致向正西北方向滑動(dòng)，總體表現(xiàn)為震中東南方向、西北方向向內(nèi)聚斂，震中西南方向、東北方向向外拉張，與中國(guó)地震局地殼應(yīng)力研究所申文豪等的強(qiáng)震記錄解算同震形變結(jié)果[8]及中國(guó)地震局第二監(jiān)測(cè)中心劉芳等得到的九寨溝Ms7.0地震的重力場(chǎng)結(jié)果保持一致[26]。綜合斷裂帶鄰近區(qū)域的同震形變場(chǎng)及斷層位移形式，可以認(rèn)為引發(fā)本次地震的斷層性質(zhì)為左旋走滑斷層，與單新建[25]、李勇[27]、任俊杰[28]等對(duì)本次地震的研究結(jié)果吻合。

3.3 精度分析

由于震中位于山區(qū)地帶，以地面觀測(cè)站數(shù)據(jù)為基礎(chǔ)的常規(guī)觀測(cè)技術(shù)在本次地震中難以實(shí)施，震中附近觀測(cè)站點(diǎn)分布稀疏，文中研究區(qū)內(nèi)只有GNSS的J416、JB33流動(dòng)觀測(cè)站和GPS的51JZZ、51JZB測(cè)站，點(diǎn)位數(shù)目較少[8,29]。很多學(xué)者在九寨溝地震研究過(guò)程中選擇了InSAR方法獲取視線向的形變數(shù)據(jù)[8,24-25,30-31]，其反演結(jié)果如表2所示。由于選用衛(wèi)星、運(yùn)行方向、時(shí)空基線、數(shù)據(jù)處理方式不同等原因，形變結(jié)果存在一定的差異，整理了地震視線向升降軌結(jié)果得到：升軌數(shù)據(jù)得到的隆升形變區(qū)間為0.07～0.11 m、沉降形變區(qū)間為 0.18～0.22 m；降軌數(shù)據(jù)得到的隆升形變區(qū)間為0.10～0.16 m、沉降形變區(qū)間為0.05～0.10 m，本次實(shí)驗(yàn)結(jié)果均在區(qū)間內(nèi)，可以認(rèn)為由此得到的三維形變結(jié)果是可信的。

表2 不同研究者對(duì)九寨溝地震視線向形變的反演結(jié)果

4 結(jié) 論

以四川省九寨溝7.0級(jí)地震為研究對(duì)象，融合升降軌Sentinel-1A數(shù)據(jù)，使用DInSAR和MAI結(jié)合的方法對(duì)該地震的雷達(dá)視線向和軌道運(yùn)行方位向形變進(jìn)行了提取，構(gòu)建了九寨溝Ms7.0地震三維同震形變場(chǎng)，結(jié)果表明：九寨溝7.0級(jí)地震在震源中心的西北側(cè)表現(xiàn)為沉降運(yùn)動(dòng)，在震源中心的東南側(cè)表現(xiàn)為隆升運(yùn)動(dòng)，地表運(yùn)動(dòng)以水平形變?yōu)橹鳎v向起伏較小，其中，東西方向?qū)Ρ敬蔚卣鸬耐鹦巫冐暙I(xiàn)最大，最大形變值為0.25 m。綜合斷裂帶鄰近區(qū)域的同震形變場(chǎng)及斷層位移形式，可判斷引發(fā)本次九寨溝Ms7.0地震的斷層性質(zhì)為左旋走滑斷層。另外，得到的地表位移三分量，更加準(zhǔn)確地反演了地震真實(shí)地表的同震形變情況，也可為地震發(fā)震機(jī)理和地震烈度判斷、震后重建等后續(xù)工作提供三維數(shù)據(jù)支持，具有一定的實(shí)用價(jià)值。文中將觀測(cè)精度近似認(rèn)為相等并選用了單位矩陣進(jìn)行權(quán)重評(píng)定，這種近似存在一定的誤差，三維形變場(chǎng)的提取模型的優(yōu)化有待進(jìn)一步研究。