運用InSAR技術(shù)解算四川長寧MS6.0地震三維形變場及時間序列形變分析

于書媛, 方良好, 宴金旭, 倪紅玉, 陳 靚, 丁 娟

(1. 安徽省地震局, 安徽 合肥 230031; 2. 安徽蒙城地球物理國家野外科學(xué)觀測研究站, 安徽 蒙城 233527;3. 四川省地震局, 四川 成都 610041)

0 引言

北京時間2019年6月17日22時55分43秒,四川長寧發(fā)生6.0級地震。參考震中位于104.90°E,28.34°N,震源深度16 km。重慶、四川、貴州、云南等多省對此次地震有感。截至2019年6月26日8時,共記錄到大于2.0級余震182次,其中大于MS4.0余震7次,最大余震震級5.4級。此次地震最大烈度達(dá)到Ⅷ度,截至2019年6月21日,地震共造成13人死亡,226人受傷。自長寧MS6.0地震發(fā)生后,相關(guān)地震行業(yè)專家從波形資料、地殼形變、重力場、地電場、地磁場、地質(zhì)構(gòu)造等方面對地震構(gòu)造應(yīng)力場特征、孕震機理、構(gòu)造特征進(jìn)行深入探索[1-3]。我國地震多發(fā)生在高原山區(qū),震中附近地質(zhì)構(gòu)造復(fù)雜,植被覆蓋,人不能至,傳統(tǒng)手段無法獲取連續(xù)、宏觀的形變場。合成孔徑雷達(dá)干涉測量技術(shù)(Interferometric Synthetic Aperture Radar,InSAR)自1993年首次應(yīng)用于地震以來,已在全球地震研究中得到快速發(fā)展和成熟應(yīng)用。本文利用覆蓋長寧縣地區(qū)的哨兵1A升降軌SAR數(shù)據(jù)(http://www.asf.allaska.edu/sar-data/palsar/),采用InSAR技術(shù)獲取地震的三維同震形變場,并進(jìn)一步研究震源位置,反演斷層幾何參數(shù)和近6個月的形變累積場,討論其發(fā)震原因與雙河背斜構(gòu)造關(guān)系,彌補該地區(qū)活動構(gòu)造研究的不足,為判定區(qū)域地震形勢提供參考。

1 研究區(qū)及InSAR數(shù)據(jù)

1.1 研究區(qū)概況

本文研究區(qū)以長寧縣地震震中為中心,覆蓋面積約為30 km×30 km,長寧地震震中位于青藏高原東緣、云貴高原北緣地殼運動較為活動的區(qū)域,地貌位置表明西南側(cè)是抬升的高原,北側(cè)是相對沉陷的四川盆地。此次地震發(fā)生于川西南地區(qū),其斷裂活動程度及地震強度雖不能和橫斷山東部及其川滇地震帶相比,但仍然是中強震活動較為頻繁的區(qū)域,研究區(qū)構(gòu)造情況如圖1所示。

1.2 研究數(shù)據(jù)

歐洲航天局(ESA)2014年4月3日發(fā)射的Sentinel-1A衛(wèi)星采用12天重訪周期進(jìn)行全球覆蓋,其獲取了長寧地震震前震后多景SAR影像。具體獲取參數(shù)見表1。采用的外部參考DEM數(shù)據(jù)為NASA SRTM DEM,數(shù)據(jù)精度為30 m。

表1 升降軌差分干涉影像參數(shù)

時間序列研究選取降軌2019-07-26作為超級主影像,其他影像作為副影像,具體獲取參數(shù)見表2。分析12景SLC影像的時間基線距和垂直基線距,在已有SAR影像數(shù)據(jù)集中形成基于不同主影像的時間序列干涉圖子集,利用短基線合成孔徑雷達(dá)干涉測量(Small Baseline Subset-InSAR,SBAS-InSAR)求解近6個月的震間形變累積、削弱信息。

表2 時間序列Sentine-1A影像參數(shù)

2 InSAR數(shù)據(jù)處理及三維同震形變場

常規(guī)單一的一對數(shù)據(jù)進(jìn)行差分干涉測量(Differential InSAR,D-InSAR)處理,獲得的是雷達(dá)視線向的形變,可能存在獲取的差分干涉結(jié)果形變量差異較大的情況,因此單一方向的形變結(jié)果不能準(zhǔn)確反映研究區(qū)的真實三維形變情況。本文結(jié)合升、降軌影像數(shù)據(jù),采用D-InSAR方法獲取沿LOS向形變信息并進(jìn)行三維同震形變場聯(lián)合解算,獲取同震東西、南北、垂直三個方向上的形變量。

2.1 D-InSAR獲取LOS向同震形變

本文運用震前震后各2景升軌、降軌Sentinel-1A幅寬約為250 km數(shù)據(jù)通過截取研究區(qū)經(jīng)差分干涉處理獲取震區(qū)三維同震形變場,其中升軌組成的干涉對為2019-06-09和2019-06-21、2019-06-14和2019-06-26,降軌組成的干涉對為2019-06-16和2019-06-28,時間基線均為12 d。利用商業(yè)軟件ENVI的SARscape模塊對Sentinel-1A數(shù)據(jù)進(jìn)行SAR差分干涉技術(shù)(D-InSAR)處理[4],由于選用的三對干涉對具有較高的相干性,獲取了長寧地震的形變場。從圖2中可看出,升軌的雷達(dá)視線向形變場整體以沉降為主,降軌的視線向形變場整體以抬升為主,且升降軌衛(wèi)星數(shù)據(jù)同區(qū)域形變場觀測值符號相反,其形變場范圍大致為18 km×15 km,形狀似橢圓形(長軸大體呈北西向),區(qū)域地質(zhì)構(gòu)造與長寧背斜構(gòu)造附近規(guī)模不大的次級斷層對應(yīng)。衛(wèi)星LOS向最大沉降量(遠(yuǎn)離衛(wèi)星傳感器)分別是7.3 cm(升軌)和6.26 cm(降軌),最大抬升量(靠近衛(wèi)星傳感器)分別是5.5 cm(升軌)和8.12 cm(降軌)。綜合本次地震的發(fā)震構(gòu)造和升降軌LOS向沉降和抬升之間的關(guān)系,判定本次形變場是由斷層的走滑和逆沖綜合運動形成。根據(jù)形變場抬升量、沉降量分布及區(qū)域地質(zhì)構(gòu)造背景,確定一條北西—南東走向的斷層跡線(圖2)。結(jié)合發(fā)震地區(qū)1∶5萬地質(zhì)圖和地質(zhì)構(gòu)造背景分析,長寧地震發(fā)生在四川盆地南緣盆山轉(zhuǎn)換帶,發(fā)震區(qū)內(nèi)主要地質(zhì)構(gòu)造為褶皺及伴生斷層,出露較少,走向和規(guī)模多變。本文形變場的整體呈NW-SE走向,其分布特征與長寧—雙河復(fù)式大背斜主體構(gòu)造一致,認(rèn)為本次形變場推斷斷層水平位置與雙河鄉(xiāng)處存在的NW向伴生斷層存在重合。圖3為跨斷層沿AA′剖面的升軌和降軌位移剖面圖,其中升軌LOS向形變量約為-4～1 cm,降軌LOS向形變量約為-1～4 cm。

圖2 D-InSAR雷達(dá)視線向形變Fig.2 D-InSAR LOS deformation

圖3 沿圖2剖面AA′的升降軌形變Fig.3 Ascending and descending deformation along the profile AA′ in figure 2

2.2 三維同震形變場解算

采用常規(guī)D-InSAR方法對升、降軌SAR數(shù)據(jù)對提取的雷達(dá)視線向形變量為一維方向的形變量,由于一維的形變量是由東西、南北、垂直三個維度方面投影得到的,因此通過D-InSAR技術(shù)處理得到的形變信息存在視線向模糊問題[5-8]。本文為解算上述3個方向的形變分量,設(shè)定dE、dN、dU為同震形變在東、北及垂直方向上的分量,排除衛(wèi)星軌道誤差、測量誤差等因素,通過地表的空間幾何與衛(wèi)星視線和軌道的關(guān)系得到函數(shù)式,如式(1)。

dU·cosθ-dN·cosα·sinθ-dE·sinαsinθ=dLOS

(1)

式中:θ為雷達(dá)視線方向的入射角;α為方位角;dE、dN、dU分別為分解后同震形變在東南、西北、垂直方向上的分量;dLOS差分干涉測量的同震形變量。同時,為了解算同震形變?nèi)S方向的形變量,原則上需要3個視線方向的形變結(jié)果。因此,本文為解算更可靠的3D分解同震形變結(jié)果,結(jié)合2個升軌LOS向同震形變場和1個降軌LOS向同震形變場,采用加權(quán)最小二乘法反演(權(quán)重來自于LOS位移精度)。其中,LOS位移精度由相干性和波長導(dǎo)出,提供測量精度的估計值(即標(biāo)準(zhǔn)偏差值)。這個估計的標(biāo)準(zhǔn)偏差值越高,測量值精確度越低,精度公式表達(dá)為:

(2)

式中:γ表示干涉相干性;λ表示波長。

本文計算出的三維形變場EW向形變量(向E為正)、SN向形變量(向N為正)、垂直向形變量(向上為正)。長寧地震三維同震形變場解算結(jié)果顯示,如圖4(c)所示,地震形變以向北運動為主,向北最大偏移量達(dá)到17.8 cm。結(jié)合圖4(a)、(b)所示,橢圓抬升區(qū)域及其西南側(cè)大部分沉降區(qū)域向E移動,最大偏移量為8 cm;北東側(cè)沉降區(qū)則向W運動,最大偏移量為7 cm,該同震形變場符合逆斷層特征。由于地震烈度圖可反映地震區(qū)域內(nèi)建筑物的破壞范圍、破壞程度以及地表變化情況[9]。本文將四川地震局發(fā)布的長寧6.0級地震烈度圖與InSAR監(jiān)測結(jié)果進(jìn)行對比可以看出,烈度范圍與InSAR技術(shù)獲取的地表形變區(qū)域空間分布較為一致,發(fā)震斷層總體沿NW向延伸。三維同震形變中UD向分量結(jié)果表明,該區(qū)域存在NW向共軛的2條長度約15 km的斷裂,呈現(xiàn)中心抬升、兩側(cè)下沉的特點,形變場特征走向也為NW向。

圖4 三維同震形變Fig.4 3D coseismic deformation

綜合分析升降軌LOS向形變和三維同震形變特征可知,震中參考位置位于橢圓形變場的西南區(qū)域,結(jié)合研究區(qū)行政區(qū)劃圖和地震烈度圖可知,橢圓形形變場區(qū)域為受到長寧地震影響最大的區(qū)域,該區(qū)域分布在長寧縣雙河鎮(zhèn)、富興鄉(xiāng),其UD向形變量最大約為8 mm。橢圓形抬升區(qū)域兩側(cè)為地震造成的下沉區(qū)域,主要在長寧雙河鎮(zhèn)附近,結(jié)合相關(guān)應(yīng)急救援資料可知,長寧地震造成的雙河鎮(zhèn)傷亡人員數(shù)量、人工建(構(gòu))筑物的損壞,這與InSAR監(jiān)測結(jié)果吻合。

3 同震滑動分布反演

本文運用兩軌升降軌同震形變場等權(quán)為反演約束條件,基于Okada彈性半空間位錯模型[10],反演地震斷層滑動分布。同震滑動震源參數(shù)反演步驟包括:

(1)對升軌2019-06-09和2019-06-21、降軌2019-06-16和2019-06-28得到的D-InSAR形變場通過繪制矢量多邊形圖層,劃定形變場和反演模型范圍,形變場采樣間隔是1 km,反演模型范圍采樣間隔為2 km,升降軌形變場分別采樣得到1 320、976個數(shù)據(jù)點,兩軌數(shù)據(jù)等權(quán)重值作為反演斷層滑動的約束條件,震中形變場的降采樣結(jié)果、模擬形變場和殘差分布如圖5所示。圖5(a)、圖5(b)、圖5(c)和圖5(d)、圖5(e)、圖5(f)分別為升軌和降軌的觀測值、模擬值和殘差,升降軌誤差集中在-1～1 cm區(qū)間。

(2)將得到的升、降軌InSAR視線向形變場結(jié)果,采用非線性方法反演均勻滑動斷層幾何參數(shù)。為了準(zhǔn)確進(jìn)行反演,加入6個軌道參數(shù)進(jìn)行線性估計軌道誤差,Okada模型參數(shù)設(shè)置采用GCMT的震源機制解進(jìn)行約束。其中,采用Levemberg-Marquardt最小二乘優(yōu)化算法迭代,進(jìn)行9個幾何參數(shù)和6個軌道參數(shù)求解。

(3)最后基于上一步反演的斷層參數(shù)采用線性反演方法并結(jié)合非負(fù)最小二乘算法求解活動斷層精細(xì)滑動分布。斷層參數(shù)反演結(jié)果如表3所列,發(fā)震斷層二三維模型如圖6所示。圖6(a)、圖6(b)分別為坐標(biāo)系下二維同震滑動分布、二維斷層滑動分布。通過斷層幾何參數(shù)反演結(jié)果表明本次地震的宏觀震中為104.87°E,28.39°N,矩震級為MW5.8,破裂尺度模型斷層長度20 km,寬度5 km,破裂深度主要集中在2～13 km范圍內(nèi),震源深度10.2 km;斷層傾角為56.42°,平均滑動角為44.37°,斷層走向319.3°。此次地震時兼具逆沖和走滑的破裂模式,對比長寧地震烈度圖可看出反演出斷層的走向與烈度圈的方向一致。此外,如圖6所示,基于DInSAR形變結(jié)果反演的地震形變宏觀震中與多家研究結(jié)構(gòu)基于地震波確定的震中位置具有比較大的差異,分析原因可能與臺站的分布情況、觀測手段及精度對研究機構(gòu)的速報結(jié)果精度均有影響。同震形變場得到的宏觀震中更能真實地反映中小震源發(fā)震位置,為地震監(jiān)測預(yù)報中心提供更加準(zhǔn)確的斷層位置和幾何信息。

圖5 形變場采樣及反演殘差Fig.5 Sampling results and inversion residuals of deformation field

圖6 長寧地震同震破裂模型的滑動分布Fig.6 Co-seismic slip distribution of rupture model in Changning earthquake

圖7 震中位置Fig.7 Epicenter location

4 SBAS-InSAR震間形變分析

4.1 時間序列形變

時間序列InSAR觀測結(jié)果可從形成場時空分布特征和斷層活動狀態(tài)進(jìn)行活動斷裂帶地殼運動的精細(xì)結(jié)構(gòu)研究,并可從斷裂帶及近場區(qū)域的宏觀面狀動態(tài)形變圖像研究斷層震前-震間-震后累積形變量、活動速率等[10-11]。本文結(jié)合研究區(qū)域的相干性和使用的數(shù)據(jù)量,選擇SBAS-InSAR方法[12],計算平均形變結(jié)果。

4.2 連續(xù)空間形變場模擬

克里金插值是基于包含自相關(guān)(測量點之間的統(tǒng)計關(guān)系)的統(tǒng)計模型[13]。本文采用克里金插值方法對離散的形變場矢量點進(jìn)行插值,圖8為由SBAS時間序列方法解算和克里金插值方法得到的形變場結(jié)果。將2019-03-12作為未發(fā)生形變的時間基準(zhǔn),以11個時間節(jié)點的累積形變量來展示此次長寧地震震前、震后的地表LOS向形變動態(tài)演化過程。分析結(jié)果表明,2019-03-24—2019-06-16時間段內(nèi),該區(qū)域范圍基本處于較為微小的形變量,形變場范圍累積形變量較小,震中近場附近最大累積形變量約為±8 mm。2019-06-16—2019-06-28時間為發(fā)震前后時刻(2019-06-17)。由圖8(h)可知,該時間震中附近累積形變量突增,最大值超過±20 mm。2019-06-28—2019-07-22期間,形變場累積形變量略有衰減。2019-08-03—2019-08-15期間,震中附近累積形變場又出現(xiàn)增長,最大抬升值超過30 mm,最小沉降值超過20 mm,分析可能由于本次地震的余震較為豐富是的震后地表形變處于一定時間的不穩(wěn)定狀態(tài)。

圖8 研究區(qū)形變累積時間序列Fig.8 Accumulative time series of deformation in the study area

為較為直觀地分析震間形變過程,選取圖8(i)中具有典型形變特征的A、B兩點分析形變時間序列,A點位于抬升區(qū)域,B點位于沉降區(qū)域。圖8中紅線為長寧地震的發(fā)震時刻(2019-06-17),從圖9可看出,紅線所處的時間段為形變速率最大的一段。地震發(fā)生前A點的形變波動表現(xiàn)為緩慢抬升,沉降量略大于抬升量,而B點沉降量始終大于抬升量,總體呈現(xiàn)沉降趨勢,震中近場其他點呈現(xiàn)類似的形變趨勢。

圖9 圖8中A點和B點形變時間序列Fig.9 Deformation time series at points A and B in figure 8

5 討論

本次地震為典型的構(gòu)造地震活動,區(qū)內(nèi)僅有少量斷層出露地表,地震發(fā)生于四川盆地邊緣的長寧背斜構(gòu)造規(guī)模不大的次級斷層上,在已有的活動斷層數(shù)據(jù)庫中并無標(biāo)識,給發(fā)震構(gòu)造判識帶來很大困難[14],是一次典型的盲斷層破裂事件。發(fā)震區(qū)域位于斷裂活動性較弱、分布式構(gòu)造變形的川西南區(qū)域,西側(cè)為長800 km的華鎣山褶皺帶,形成北西—南東向的現(xiàn)今應(yīng)力場[15],光學(xué)影像上沒有顯著的斷層構(gòu)造特征。DInSAR同震形變場表明,本次地震產(chǎn)生了18 m×15 km的形變場范圍,衛(wèi)星視線向最大沉降量是7.3 cm(升軌)、最大抬升量是8.12 cm(降軌),斷層形變場走向為NW向,最大偏移量達(dá)到17.8 cm(向北)。反演結(jié)果表明,滑動量在斷層深度2～6 km處達(dá)到最大,以逆斷走滑為主,最大滑動量約0.06 m,位于沿斷層面傾斜陷落深度3 km處,震級為MW5.8。

本研究得到的發(fā)震斷層為逆斷層,根據(jù)形變上升的主動盤位于西側(cè)山脊區(qū)域,而被動盤位于東側(cè)海拔較低的低洼區(qū)域,同震形變與地貌特征耦合。根據(jù)郭志[16]的研究顯示長寧地震主震級幾次震級較大余震的震源機制結(jié)果表明,發(fā)震斷層以逆沖為主,結(jié)果與長寧—雙河背斜的地質(zhì)構(gòu)造及應(yīng)力背景較為一致,斷層傾角為59°,略大于本研究中得到的發(fā)震斷層的傾角(45°)。

截至2019年7月28日,長寧地震主震后發(fā)生余震6 888次,其中包含6次M>4.5的余震,相對于主震呈北西向條帶分布,總長約20 km,正對應(yīng)構(gòu)造應(yīng)力較為集中的背斜軸部區(qū)域,存在眾多背斜伴生斷層[2,17]。在此次地震震中南側(cè)約15 km處曾發(fā)生2018年12月16日四川興文M5.7和2019年1月3日珙縣M5.3地震。易桂喜等[2]通過長寧6.0級地震序列16次MS≥3.6地震序列震源機制解統(tǒng)計認(rèn)為發(fā)震斷層高傾角,斷層走向以NW-SE和NS走向略占優(yōu)勢,發(fā)震構(gòu)造復(fù)雜,構(gòu)造應(yīng)力環(huán)境有利于斷裂發(fā)生逆沖運動。尹欣欣等[17]利用四川長寧6.0級地震的主余震資料,利用CAP震源機制解方法對14次3級以上的地震序列做了震源機制解計算,認(rèn)為應(yīng)力場反演結(jié)果顯示主震區(qū)域網(wǎng)格(28.3°N,104.8°E)斷層為逆沖類型,破裂長度大致為35 km。孫凱等[18]采用SDM反演InSAR數(shù)據(jù)源,迭代計算發(fā)震斷層面最優(yōu)滑動分布,最大滑動量0.28 m,平均滑動角50°,與本文反演結(jié)果基本一致?？傮w來說,認(rèn)為本次地震發(fā)生在上地殼淺部,發(fā)震斷層滑動面上緣接近地表,推測由主震引起的次級斷層活動觸發(fā)短期內(nèi)強余震頻發(fā),認(rèn)為該區(qū)域未來需要進(jìn)一步加強人口密集地區(qū)的地震風(fēng)險防范工作。

6 結(jié)論

本文基于升降軌Sentinel-1A數(shù)據(jù),采用DInSAR技術(shù)獲取2019年6月17日四川長寧MS6.0地震升降軌同震形變場和加權(quán)最小二乘法反演三維同震形變場。分析結(jié)果表明,長寧縣地震形變升軌LOS向形變場以沉降為主(最大沉降量7.3 cm),降軌LOS向形變場以抬升為主(最大抬升量8.12 cm),且升降軌衛(wèi)星數(shù)據(jù)同區(qū)域形變場觀測值符號相反;三維同震形變場以NS向為主,EW向形變量約為-7～8 cm,UD向形變量約為-4～7 cm;形變場空間分布特征與斷裂走向為NW向,且與地震烈度圖的空間分布特征較為一致。利用升降軌形變場采用均勻滑動模型反演結(jié)果表明斷層幾何參數(shù)反演結(jié)果表明本次地震的破裂尺度模型斷層長度20 km,寬度5 km,破裂深度主要集中在2～13 km范圍內(nèi),震源深度10.2 km;斷層傾角為56.42°,平均滑動角為44.37°,斷層走向319.3°,與CMT、USGS、中國地震局地質(zhì)研究所、孫凱等結(jié)果基本一致[18]。SBAS-InSAR震間形變結(jié)果表明,該區(qū)域近6個月內(nèi)區(qū)域累積形變場變化特征主要由長寧地震引起,形變場動態(tài)演化特征明顯,震前時間段該區(qū)域范圍基本處于較為微小的形變量,形變場范圍累積形變量較小,臨震時間段震中附近累積形變量突增,震后短期一段時間形變場累積形變量略有衰減,但最后一段時間震中附近累積形變場又出現(xiàn)增長,最大抬升值超過30 mm,最小沉降值超過20 mm,分析可能由于本次地震的余震較為豐富是的震后地表形變處于一定時間的不穩(wěn)定狀態(tài)。