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

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

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

0 引言

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

1 研究區及InSAR數據

1.1 研究區概況

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

1.2 研究數據

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

表1 升降軌差分干涉影像參數

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

表2 時間序列Sentine-1A影像參數

2 InSAR數據處理及三維同震形變場

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

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

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

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

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

2.2 三維同震形變場解算

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

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

(1)

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

(2)

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

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

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

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

3 同震滑動分布反演

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

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

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

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

圖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觀測結果可從形成場時空分布特征和斷層活動狀態進行活動斷裂帶地殼運動的精細結構研究,并可從斷裂帶及近場區域的宏觀面狀動態形變圖像研究斷層震前-震間-震后累積形變量、活動速率等[10-11]。本文結合研究區域的相干性和使用的數據量,選擇SBAS-InSAR方法[12],計算平均形變結果。

4.2 連續空間形變場模擬

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

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

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

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

5 討論

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

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

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

6 結論

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