基于Sentinel-1A的2021年雙湖MW5.7地震同震破裂模型

夏 濤 劉小利 余鵬飛 鄧德貝爾 樂子揚(yáng) 高天琪

1 中國(guó)地震局地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室，武漢市洪山側(cè)路40號(hào)，430071 2 防災(zāi)科技學(xué)院，河北省三河市學(xué)院街465號(hào)，065201

2021-03-30西藏雙湖縣黃水湖西岸發(fā)生MW5.7地震，中國(guó)地震臺(tái)網(wǎng)中心(CENC)給出震中位置為34.38°N、87.68°E，震源深度為10 km。雙湖縣位于青藏高原腹地的羌塘塊體中部，自1970年有地震記錄以來(lái)，羌塘塊體中部共發(fā)生過14次5級(jí)以上地震，本次地震是周邊100 km范圍內(nèi)的首個(gè)5級(jí)以上地震(圖1(a))。由于自然環(huán)境惡劣、交通限制等原因，羌塘塊體中部測(cè)震臺(tái)和GNSS連續(xù)站非常稀疏，地質(zhì)學(xué)資料有限，距離震中200 km范圍內(nèi)除震中西側(cè)約1 km處有一條近SN向、傾角不明的無(wú)名正斷層外(本文稱黃水湖斷層)，尚未有其他已知斷層[1-2]。受限于震中區(qū)域地質(zhì)構(gòu)造、余震和同震GPS等資料較少，發(fā)震斷層的確定較難，CENC未給出此次地震的震源機(jī)制解，而GCMT和USGS給出的震源機(jī)制解差異較大，發(fā)震斷層面選擇不統(tǒng)一。

斷裂數(shù)據(jù)來(lái)源于文獻(xiàn)[1-2]；震源機(jī)制解源于哈佛大學(xué)GCMT；歷史地震源于CENC；GPS速度場(chǎng)源于文獻(xiàn)[3]圖1 羌塘塊體中部歷史地震及斷裂分布Fig.1 The distribution of the faults and historical earthquakes of the central Qiangtang block

InSAR可大范圍、高精度地獲取地震同震形變場(chǎng)，是確定中強(qiáng)地震發(fā)震構(gòu)造的重要手段之一。本文利用Sentinel-1A SAR-C地震前后升軌和降軌4景數(shù)據(jù)獲取2021年雙湖地震同震形變場(chǎng)，基于半空間位錯(cuò)模型確定發(fā)震構(gòu)造幾何學(xué)參數(shù)，計(jì)算不同深度上的同震靜態(tài)庫(kù)侖應(yīng)力變化，并結(jié)合羌塘塊體構(gòu)造背景討論本次地震的發(fā)震機(jī)制及其對(duì)區(qū)域未來(lái)地震危險(xiǎn)性的影響，為深入理解羌塘塊體中部構(gòu)造活動(dòng)性和地殼變形機(jī)制提供參考依據(jù)。

1 區(qū)域構(gòu)造背景

GPS速度場(chǎng)(圖1(b))顯示，青藏高原內(nèi)部各次級(jí)塊體呈現(xiàn)不同的地殼形變機(jī)制[4-5]。班公湖-怒江縫合帶以北的羌塘塊體從西向東滑移速率逐漸減小，羌塘塊體中部至東部滑移速率相差可達(dá)10.3 mm/a[6]，表現(xiàn)出明顯的空間差異性；班公湖-怒江縫合帶以西在上地殼物質(zhì)持續(xù)向東擠出的構(gòu)造作用下，沿羌塘內(nèi)部小型正斷層和共軛走滑斷裂發(fā)生彌散型變形特征，構(gòu)造活動(dòng)以EW向延伸擴(kuò)展為主；班公湖-怒江縫合帶以東則直接向東流出，使得塊體東邊界不斷向東擴(kuò)張，發(fā)生剛性塊體變形，呈現(xiàn)明顯的物質(zhì)快速轉(zhuǎn)移運(yùn)動(dòng)特征[4]。

目前，對(duì)羌塘塊體的斷裂活動(dòng)性仍缺乏定量研究，特別是地處青藏高原腹地的塊體中部，僅Wang等[7]利用時(shí)序InSAR技術(shù)發(fā)現(xiàn)羌塘塊體中部存在持續(xù)向東的滑移運(yùn)動(dòng)，速率約為4 mm/a。2021年雙湖發(fā)震發(fā)生于羌塘塊體中部的黃水湖西岸，其發(fā)震構(gòu)造和變形機(jī)制不明，變形特征是呈彌散型還是剛性整體向東擠出均有待進(jìn)一步研究。

2 InSAR數(shù)據(jù)處理

本文利用Sentinel-1A影像(表1)，基于兩軌差分干涉獲取本次雙湖地震的干涉影像，提取相應(yīng)的同震形變場(chǎng)?；陂_源雷達(dá)數(shù)據(jù)處理平臺(tái)GMTSAR[8]進(jìn)行以下數(shù)據(jù)處理：1)采用20∶4的多視比例；2)通過GACOS移除影像中的大氣誤差[9]；3)根據(jù)精密軌道和SRTM提供的30 m DEM去除軌道誤差和平地相位；4)采用Goldstein濾波方法和最小費(fèi)用流解纏方法分別對(duì)干涉圖進(jìn)行濾波和相位解纏處理；5)地理編碼。

表1 干涉對(duì)信息Tab.1 Interferometry pairs information

主震升降軌干涉圖和同震形變場(chǎng)如圖2所示?？梢钥闯?，形變場(chǎng)以近SN向分界呈三瓣式(L1、L2、L3)分布。L1部分在升降軌形變場(chǎng)中最大視線向(line of sight，LOS)形變量分別為0.53 cm、-0.46 cm，二者符號(hào)相反、量值接近；升軌時(shí)遠(yuǎn)離、降軌時(shí)靠近L2部分，指示L1存在NWW向水平運(yùn)動(dòng)。L2部分在升降軌形變場(chǎng)中皆呈現(xiàn)負(fù)值，LOS向形變量分別為-6.78 cm、-4.73 cm，表明該區(qū)域整體為垂直沉降運(yùn)動(dòng)。L3部分在升降軌形變場(chǎng)中LOS形變量分別為-1.65 cm、1.45 cm，符號(hào)相反、量值接近；升軌時(shí)接近、降軌時(shí)遠(yuǎn)離L2部分，指示L3存在近東向水平運(yùn)動(dòng)。圖2(g)顯示，沿AB測(cè)線得到的剖面T114A在從L1向L2過渡時(shí)(L1與L2之間的虛線)，其梯度明顯大于剖面T121D，可能是由于T121D的成像視角與干涉圖長(zhǎng)軸方向近乎平行，導(dǎo)致真實(shí)地表形變?cè)谛巫儓?chǎng)中的投影分量較少，限制了LOS形變的觀測(cè)精度。表2中最佳均勻滑動(dòng)模型的傾角為本次雙湖事件形變場(chǎng)長(zhǎng)軸為近SN走向，主要以垂直向形變?yōu)橹鳌Ｒ騃nSAR對(duì)SN方向形變分量不敏感，因此在忽略SN向形變分量的情況下，利用不同視角的升降軌同震形變場(chǎng)提取震區(qū)EW水平位移場(chǎng)和垂直向位移場(chǎng)[10](圖2(c)、2(f))更適用于本次地震。計(jì)算得出，EW向和垂直向形變場(chǎng)的殘差(RMS)分別為0.4 cm和1.3 cm，在InSAR精度范圍內(nèi)，說(shuō)明2.5D形變場(chǎng)具有一定的可靠性。從圖2(c)水平向形變場(chǎng)來(lái)看，L1表現(xiàn)為西向水平位移，最大位移值為2.5 cm；L2、L3則為東向水平位移，最大位移值為5.8 cm，且出現(xiàn)在L2區(qū)間；從圖2(f)垂直向形變場(chǎng)來(lái)看，最大沉降量達(dá)11.8 cm，出現(xiàn)在L2區(qū)間；L1、L3部分呈微弱抬升，最大抬升量為2 cm。此外，對(duì)比3組形變場(chǎng)來(lái)看，在L1與L2的過渡區(qū)域，與T121D形變場(chǎng)相比，T114A表現(xiàn)出更為明顯的近SN向長(zhǎng)軸走向。由此推斷，造成T114A和T121D沿測(cè)線形變量變化梯度差距較大(圖2(g))的可能原因是T121D成像角度與形變場(chǎng)近SN向長(zhǎng)軸幾乎垂直，對(duì)形變場(chǎng)水平向的位移并不敏感，且L1部分本身水平位移量較小。

暖色代表接近衛(wèi)星方向或抬升運(yùn)動(dòng)，冷色代表遠(yuǎn)離衛(wèi)星方向或沉降運(yùn)動(dòng)圖2 雙湖地震同震形變場(chǎng)Fig.2 Coseismic deformation field of the Shuanghu earthquake

3 斷層滑動(dòng)分布反演

為精細(xì)地確定發(fā)震斷層位錯(cuò)滑移特征，采用以下步驟進(jìn)行計(jì)算：1)以同震形變場(chǎng)為輸入數(shù)據(jù)，利用蒙特卡洛算法搜索最佳斷層參數(shù)，獲得發(fā)震斷層的均勻滑動(dòng)模型；2)以最佳均勻滑動(dòng)模型參數(shù)為先驗(yàn)條件約束，利用Okada位錯(cuò)模型反演同震滑動(dòng)分布。

3.1 均勻滑動(dòng)反演

為提高最佳斷層參數(shù)的搜索效率，確保有效收斂，使用四叉樹采樣算法[11]對(duì)升降軌同震形變場(chǎng)進(jìn)行降采樣處理。處理后得到494個(gè)升軌、422個(gè)降軌形變點(diǎn)位，以此為蒙特卡洛搜索約束[12]，獲得收斂性較好的馬爾科夫鏈(圖3)，并統(tǒng)計(jì)特征最佳的均勻滑動(dòng)模型斷層參數(shù)(表2)。

紅線表示最佳斷層模型參數(shù)圖3 蒙特卡洛-馬爾科夫鏈聯(lián)合概率密度分布Fig.3 Distribution of joint probability density of Monte Carlo-Markov chain

表2 均勻滑動(dòng)模型Tab.2 Uniform slip model

表2中最佳均勻滑動(dòng)模型的傾角為51°，走向?yàn)?6.6°，走向與主震干涉圖NNE向長(zhǎng)軸走向一致；斷層走滑量為0.39 m，傾滑量為-0.61 m，指示雙湖地震可能是一次正斷為主、兼右旋走滑的地震，但已有研究并未表明該區(qū)域存在右旋走滑斷層[13]。因此，有必要構(gòu)建滑動(dòng)分布模型，進(jìn)一步正演精細(xì)的三維形變特征以確定發(fā)震斷層幾何學(xué)參數(shù)。

3.2 斷層滑動(dòng)分布反演

以均勻滑動(dòng)模型得到的最佳斷層模型參數(shù)為先驗(yàn)條件，基于Okada位錯(cuò)模型[14]進(jìn)行同震斷層滑動(dòng)分布反演[15]，使用Crust1.0模型，反演約束方程為：

式中，G為聯(lián)系觀測(cè)值與彈性位錯(cuò)模型的格林函數(shù)矩陣，d為觀測(cè)值，?2為拉普拉斯算子，γ2為平滑因子。

根據(jù)均勻滑動(dòng)模型，設(shè)定以下初始參數(shù)：斷層長(zhǎng)18 km、寬20 km、傾角50°、走向30°、滑動(dòng)角范圍-130°～-50°。將斷層按1 km×1 km劃分為360個(gè)子斷層，并且設(shè)定平滑因子γ2為0.07。最終得到的斷層滑動(dòng)模型為：震中34.37°N、87.71°E，震源深度6.51 km，走向33.0°、傾角50°、平均滑動(dòng)角-74°、最大滑動(dòng)量0.26 m(圖4)。從圖4(a)～4(c)可見，破裂集中在沿走向8～16 km(滑動(dòng)量大于10 cm)、地下深度5～15 km處，斷層滑動(dòng)以傾滑為主，兼有微弱左旋分量。假定剪切模量為30 GPa，計(jì)算得到地震矩張量為4.5×1017Nm，相當(dāng)于震級(jí)MW5.74。

圖(c)中，視線方位角為150°，高度角為20°，加粗實(shí)線為斷層上盤；圖(d)中，箭頭表示水平運(yùn)動(dòng)方向圖4 斷層滑動(dòng)分布模型及模擬三維形變場(chǎng)Fig.4 Fault slip distribution model and the simulated 3D deformation field

圖4(d)是基于滑動(dòng)分布模型重建的矢量形變場(chǎng)。水平向位移顯示，斷層西側(cè)向西移動(dòng)，東側(cè)則微弱東移，呈左旋走滑特征；垂直向形變場(chǎng)存在明顯的NNE向分界線，以沉降運(yùn)動(dòng)為主。顯然，模擬的矢量形變場(chǎng)整體形變特征與2.5D形變場(chǎng)基本一致，且更加精細(xì)地展示了斷層在不同維度上的滑動(dòng)量值。

圖5是基于均勻滑動(dòng)模型、滑動(dòng)分布模型得到的正演模型及其擬合殘差。可以看到，T114A的2個(gè)模型與觀測(cè)值擬合度更好，RMS均為0.3 cm，殘差主要出現(xiàn)在形變場(chǎng)沉降區(qū)，主要是因黃水湖干擾而缺乏部分近場(chǎng)觀測(cè)信號(hào)導(dǎo)致的，滑動(dòng)分布模型以清晰的NNE向分界線指示了發(fā)震構(gòu)造行跡。T121D的均勻滑動(dòng)模型和滑動(dòng)分布模型擬合的RMS分別為0.3 cm和0.6 cm，但均勻滑動(dòng)模型的指示意義相對(duì)不佳。聯(lián)合T114A得到的斷層走向參數(shù)約束后，可以改善T121D滑動(dòng)分布模型在斷層走向的擬合效果(圖5(g)、5(i))，進(jìn)一步表明T121D近乎平行于斷層走向的飛行方向?qū)е略摲较蛏系男巫儷@取受限，進(jìn)而使得在該方向上的形變量低于T114A。

黑色矩形框表示滑動(dòng)分布模型斷層地表投影，加粗實(shí)線表示斷層上盤，(a)～(e)和(f)～(j)分別對(duì)應(yīng)T114A和T121D圖5 同震形變場(chǎng)觀測(cè)值、模型值與殘差Fig.5 The observation values, model values and residual of coseismic deformation field

4 討 論

均勻滑動(dòng)模型和滑動(dòng)分布模型獲得的震源機(jī)制參數(shù)見表3?？梢钥闯?，InSAR反演的模型反映的是地震破裂面的總體特征，與USGS、GCMT給出的遠(yuǎn)震體波模型結(jié)果基本一致，僅震源深度相差較大，這是由于遠(yuǎn)震體波模型地震矩張量反演中波形擬合誤差隨深度變化，深度結(jié)果誤差相對(duì)較大[16-17]。

表3 2021年雙湖MW5.7地震震源參數(shù)Tab.3 Focal parameters of the 2021 Shuanghu MW5.7 earthquake

測(cè)震學(xué)和大地測(cè)量學(xué)結(jié)果都顯示，2021年雙湖MW5.7地震為一次正斷為主、兼具走滑分量的事件。但受特定方位的遠(yuǎn)震地震資料信噪比低等影響，USGS和GCMT給出的節(jié)面I的滑動(dòng)角以及InSAR均勻滑動(dòng)模型和滑動(dòng)分布模型的滑動(dòng)角分別表現(xiàn)出右旋或左旋走滑，干擾了走滑方向的確定。已有地質(zhì)學(xué)資料顯示，羌塘塊體中部已知斷層未見右旋走滑性質(zhì)[13]。由于震中區(qū)域缺少GPS連續(xù)觀測(cè)站，測(cè)震臺(tái)也非常稀疏(圖1)，主震后2個(gè)月內(nèi)僅發(fā)生8次余震(表4，數(shù)據(jù)來(lái)源于CENC地震目錄)，無(wú)法提供更多資料約束發(fā)震節(jié)面。

表4 2021年雙湖MW5.7主震及余震信息Tab.4 2021 Shuanghu MW5.7 main shock and aftershocks information

為了解決上述問題，采用靜態(tài)庫(kù)侖應(yīng)力轉(zhuǎn)移與同震形變場(chǎng)、余震的對(duì)應(yīng)關(guān)系來(lái)討論不同性質(zhì)發(fā)震節(jié)面的差異，進(jìn)而判定發(fā)震構(gòu)造性質(zhì)，評(píng)估本次地震對(duì)區(qū)域未來(lái)地震危險(xiǎn)性的影響。考慮到InSAR滑動(dòng)分布模型給定的左旋節(jié)面解與GCMT節(jié)面I參數(shù)更接近，為對(duì)比分析，采用Crust1.0模型、摩擦系數(shù)0.4、楊氏模量30 GPa、泊松比0.25，分別以InSAR左旋節(jié)面解(走向33°、傾角50°，滑動(dòng)角-74°)和GCMT右旋節(jié)面解(走向179°、傾角54°，滑動(dòng)角-112°)為接收面計(jì)算不同深度上的靜態(tài)庫(kù)侖應(yīng)力，結(jié)果分別見圖6(a)～6(g)、6(h)～6(n)。

從圖6(a)～6(g)看出，5 km深度處，庫(kù)侖應(yīng)力呈現(xiàn)明顯的對(duì)稱特征，應(yīng)力卸載區(qū)呈EW向，應(yīng)力加載區(qū)呈SN向，東側(cè)卸載區(qū)與南側(cè)加載區(qū)分界明顯，呈NNE向，與左旋節(jié)面走向一致，左旋節(jié)面給定的主震即落于該分界處；最大應(yīng)力增加為0.6 bar，超過應(yīng)力觸發(fā)閾值0.1 bar，存在較大地震風(fēng)險(xiǎn)。7 km深度處，應(yīng)力卸載區(qū)主要沿左旋節(jié)面走向進(jìn)一步擴(kuò)張，與左旋節(jié)面給定的主震位置和SAR觀測(cè)到的主要形變區(qū)范圍(圖5)基本一致，同時(shí)南側(cè)應(yīng)力加載明顯?？紤]到InSAR左旋節(jié)面給定的主震深度為6.51 km，且在7 km處釋放了28.9%的力矩，推斷該深度為主震應(yīng)力釋放的主要深度。10 km深度處，SN向應(yīng)力加載區(qū)范圍增大，EW向應(yīng)力卸載區(qū)范圍減小，且呈東西兩瓣沿左旋節(jié)面方向反向?qū)ΨQ；最大應(yīng)力增加達(dá)0.9 bar，遠(yuǎn)大于應(yīng)力觸發(fā)閾值0.1 bar，發(fā)生于2021-04-02、震源深度為9 km的4級(jí)余震即落于東側(cè)卸載區(qū)與北側(cè)加載區(qū)的分界處；南側(cè)應(yīng)力加載也較明顯，仍存在較大地震風(fēng)險(xiǎn)。15 km深度處，SN向應(yīng)力加載區(qū)進(jìn)一步擴(kuò)張并連通，主震所在區(qū)域應(yīng)力持續(xù)卸載。在20～30 km深度處仍以應(yīng)力加載為主，最大應(yīng)力在20 km深度時(shí)增加到0.4 bar，而其他7次余震都集中發(fā)生于該深度區(qū)間，且均落于庫(kù)侖應(yīng)力大于0.1 bar的應(yīng)力加載區(qū)或邊界。

(a)～(g)為InSAR結(jié)果；(h)～(n)為GCMT結(jié)果；綠色、黃色小球分別表示不同震級(jí)的余震分布,序號(hào)與表4對(duì)應(yīng)圖6 不同接收面對(duì)應(yīng)的不同深度上同震靜態(tài)庫(kù)侖應(yīng)力變化與主余震分布Fig.6 Coseismic static Coulomb stress changes at different depths corresponding to different receivers and distribution of main shock and aftershocks

從圖6(h)～6(n)看出，應(yīng)力變化分布與左旋節(jié)面解結(jié)果存在較大差異，主要表現(xiàn)在4個(gè)方面: 1) 5～10 km深度，應(yīng)力卸載區(qū)呈NW向擴(kuò)張，與右旋節(jié)面走向不符；2) 最大應(yīng)力卸載區(qū)與SAR觀測(cè)到的主要形變區(qū)范圍(圖5)相差較大，且出現(xiàn)在7 km深度處，與表3中GCMT模型給出的主震深度均不一致；3) GCMT給出的主震深度為16.4 km，但在15 km深度處未出現(xiàn)明顯的應(yīng)力卸載區(qū)，相反在主震東側(cè)應(yīng)力加載明顯，而所有余震均未出現(xiàn)在該區(qū)域范圍；4)除2021-04-02的4級(jí)余震震源深度為10 km左右，其余7次余震震源深度都集中在25～32 km，但在25 km、30 km 深度處余震所在位置絕大部分對(duì)應(yīng)應(yīng)力卸載區(qū)。

綜上分析，由左旋節(jié)面作為接收面獲得的靜態(tài)庫(kù)侖應(yīng)力變化與SAR觀測(cè)得到的形變場(chǎng)、主余震分布更加契合，也符合羌塘塊體區(qū)域已有應(yīng)力場(chǎng)及地質(zhì)學(xué)研究結(jié)果[6,13]?？梢哉J(rèn)為，2021年雙湖MW5.7地震的發(fā)震斷層應(yīng)為正斷兼左旋走滑性質(zhì)。從圖5、圖6(a)～6(g)看出，主震和2021-04-02的4級(jí)余震發(fā)生于10 km深度以上的淺部地殼，且分布于黃水湖斷層南段；其他余震均發(fā)生于25 km以下深度，且大都散落于黃水湖斷層兩側(cè)，集中在斷層面地表投影范圍內(nèi)。分析圖6中不同深度上庫(kù)侖應(yīng)力變化及余震分布，認(rèn)為黃水湖斷層南段短時(shí)間內(nèi)地震風(fēng)險(xiǎn)性較小，北段地震風(fēng)險(xiǎn)性需結(jié)合更多資料進(jìn)一步分析。

本次地震發(fā)生于羌塘塊體中部的NNE向地塹區(qū)，InSAR形變場(chǎng)特征反映這是一次典型的地塹擴(kuò)張活動(dòng)。從區(qū)域狀態(tài)來(lái)看[6]，羌塘塊體整體張、壓、剪應(yīng)變都較強(qiáng)烈，為擴(kuò)張型地塊。從西至東，塊體SN向擴(kuò)張率逐漸增加，最大達(dá)24.40×10-9/a；而EW向擴(kuò)張率逐漸減小，低至3.21×10-9/a；塊體中部的主壓、主張、最大剪切和面應(yīng)變率分別為14.29×10-9/a、25.13×10-9/a、35.41×10-9/a和6.84×10-9/a，表明羌塘塊體持續(xù)向東延展。Yu等[18]指出，羌塘塊體東邊界可能已向東延伸至95°E，而不是Taylor等[2]認(rèn)為的93°E。因此，本次地震進(jìn)一步指示，羌塘塊體仍持續(xù)向東擴(kuò)張，因EW向拉伸作用導(dǎo)致近SN向小規(guī)模正斷層呈彌散型變形和地塹發(fā)育。

5 結(jié) 語(yǔ)

1)本次雙湖地震震中位置為34.37°N、87.71°E，震源深度6.51 km，發(fā)震斷層走向33°、傾角50°、傾向東、滑動(dòng)角-74°，以傾滑為主兼少量左旋走滑分量，最大滑動(dòng)量0.26 m；

2)短時(shí)間內(nèi)，黃水湖斷層南段地震風(fēng)險(xiǎn)性較小，北段則需結(jié)合更多資料進(jìn)一步分析；

3)本次地震是在羌塘塊體持續(xù)向東擴(kuò)張的背景下，受EW向拉伸作用使得黃水湖正斷層發(fā)生的一次彌散型變形活動(dòng)，SN向地塹得到進(jìn)一步擴(kuò)張。