基于InSAR約束的2022年瀘定MW6.6地震同震滑動分布及庫侖應(yīng)力變化

王 欣 李水平,2,3 宋舜躍

1 合肥工業(yè)大學(xué)土木與水利工程學(xué)院，合肥市屯溪路193號，230009 2 武漢引力與固體潮國家野外科學(xué)觀測研究站，武漢市洪山側(cè)路40號，430071 3 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，武漢市魯磨路388號，430074

據(jù)中國地震臺網(wǎng)測定，2022-09-05四川瀘定縣發(fā)生MS6.8地震，震中位于29.59°N、102.08°E(圖1)，震源深度16 km。瀘定地震是該地區(qū)自1955年康定7.5級地震以來的又一次強震，也是近年來造成人員傷亡最嚴重的一次地震，受到社會的廣泛關(guān)注。GCMT、USGS、IPGP、GFZ、CENC等機構(gòu)發(fā)布了此次地震的震源機制解(表1)，張龍等[1]、韓炳權(quán)等[2]也開展了相關(guān)研究。雖然所用數(shù)據(jù)和方法不同，但研究結(jié)果總體認為瀘定地震斷層的走向為NNW-SSE。

表1 瀘定地震破裂參數(shù)Tab.1 Rupture parameters of the Luding earthquake

此次瀘定地震發(fā)生在鮮水河斷裂磨西段附近，距龍門山斷裂帶、大涼山斷裂帶、安寧河斷裂帶和玉龍希斷裂帶距離較近，其中鮮水河斷裂帶、龍門山斷裂帶和安寧河斷裂帶交會構(gòu)成了“Y”字形斷裂帶。鮮水河斷裂帶地處川西高山區(qū)，是巴顏喀拉塊體的西南邊界和川滇塊體的北邊界，地質(zhì)構(gòu)造復(fù)雜，始于甘孜東谷一帶，向南經(jīng)過爐霍、道孚、康定、磨西一線，至石棉安順場一帶逐漸消失，全長約350 km，總體走向為N320°～330°W，呈略向NE凸出的弧形，是青藏高原東南緣的一條全新世活動強烈的左旋走滑型斷裂帶[3-4]。根據(jù)構(gòu)造地貌學(xué)，晚第四紀以來，鮮水河斷裂磨西段的滑動速率為9.6～13.4 mm/a[5]；根據(jù)時序InSAR方法，鮮水河斷裂帶整體滑動速率為8.8～17.9 mm/a，閉鎖深度為7.6～18.5 km[6]；根據(jù)GPS或水準測量方法估計出的鮮水河斷裂北西段滑動速率為8～15 mm/a，南段滑動速率為5～14 mm/a[7]。18世紀以來，整條鮮水河斷裂帶幾乎都發(fā)生過地表破裂，且復(fù)發(fā)頻繁，有歷史記錄的7級以上地震共8次[4](圖1)。已有研究表明，2008年汶川地震和2013年蘆山地震后，鮮水河斷裂磨西段的庫侖應(yīng)力有一定程度的加載[8]，說明該地區(qū)未來地震危險性較高。

紅色五角星為CENC震中位置，黑色震源球為GCMT、USGS、IPGP和GFZ等機構(gòu)給出的震源機制解，黃色震源球為歷史上鮮水河斷裂帶7級以上地震，深藍色震源球為龍門山斷裂帶上的蘆山和汶川地震，圓圈為余震，藍色箭頭為GPS速度場圖1 瀘定地震區(qū)域構(gòu)造背景Fig.1 Tectonic background of Luding seismic area

瀘定地區(qū)地形復(fù)雜，海拔高差大，容易發(fā)生滑坡等地質(zhì)災(zāi)害，很難進行常規(guī)大地測量以及野外現(xiàn)場考察工作。本文采用D-InSAR獲取2022-09-05瀘定地震的同震形變場，首先利用貝葉斯方法[9]中的均勻滑動模型搜索斷層的先驗幾何參數(shù)，然后采用非負最小二乘方法[10]獲得斷層面的同震滑動分布，最后通過計算同震庫侖應(yīng)力變化對區(qū)域地震危險性進行評估，計算震間應(yīng)變場對發(fā)震斷層構(gòu)造進行討論。

1 瀘定地震InSAR同震形變場

本文使用瀘定地震前后Sentinel-1A升降軌影像數(shù)據(jù)，參數(shù)見表2。使用InSAR處理軟件GMTSAR，采用二軌法進行SAR數(shù)據(jù)差分干涉處理，具體流程見圖2。

圖2 D-InSAR處理流程Fig.2 D-InSAR processing flow chart

表2 Sentinel-1A SAR影像干涉對參數(shù)Tab.2 Parameters of Sentinel-1A SAR image interference

圖3為升降軌干涉圖和同震形變場。由圖3(d)可見，2022-08-26～09-19升軌T26同震形變場最大視線向沉降量為15 cm，最大視線向抬升量為13.8 cm；由圖3(e)可見，降軌T135同震形變場最大視線向沉降量為7.4 cm，最大視線向抬升量為12.7 cm；由圖3(f)可見，2022-08-26～09-07升軌T26同震形變最大視線向沉降量為7.7 cm，最大視線向抬升量為14.9 cm。升降軌沉降量和抬升量存在差異，主要是升降軌觀測角度不同所致。本文結(jié)果與張龍等[1]和韓炳權(quán)等[2]的結(jié)果大致相同。需要說明的是，相比于2022-08-26～09-19升軌T26圖像，2022-08-26～09-07升軌T26圖像存在更大的觀測噪聲，且存在將一些未發(fā)生形變區(qū)域錯誤解纏成有相位區(qū)域的情況，因此其可靠性相對較差。為提高后續(xù)斷層滑動分布反演效率，本文對升降軌同震形變場進行四叉樹降采樣處理[11]，最終得到1 355 個2022-08-26～09-19升軌T26視線向形變觀測數(shù)據(jù)、1 567個2022-08-26～09-07升軌T26視線向形變觀測數(shù)據(jù)和1 382個降軌T135視線向數(shù)據(jù)。

圖3 InSAR升降軌干涉圖和同震形變場Fig.3 InSAR interferogram of ascending and descending and coseismic deformation field

2 斷層滑動分布反演

2.1 均勻滑動分布

由于瀘定地震缺乏發(fā)震斷層的先驗幾何參數(shù)，因此在進行分布式滑動分布反演前，需要使用均勻滑動模型獲得斷層傾角、走向、震源深度等參數(shù)。本文采用Bagnardi等[9]開發(fā)的GBIS MATLAB軟件包，基于升降軌InSAR數(shù)據(jù)，使用貝葉斯算法進行斷層參數(shù)估計。在參數(shù)估計過程中，對初始模型參數(shù)空間不加以嚴格約束，如將斷層走向設(shè)為0°～360°、傾角設(shè)為-90°～90°，可以保證搜索結(jié)果的可靠性。迭代106次后，最終得到斷層模型參數(shù)的后驗概率分布見圖4。由圖可見，斷層最優(yōu)長度為19.5±0.06 km、寬度為15.29±0.15 km、深度為15.61±0.13 km、走向為166.86°±0.09°、傾角為74.23°±0.26°、走滑分量為0.905±0.005 m、傾滑分量為-0.095±0.001 m，表明本次地震以左旋走滑為主。這與韓炳權(quán)等[2]GCMT和IPGP等給出的結(jié)果基本一致(表1)。

圖4 斷層模型參數(shù)的后驗概率分布Fig.4 The posterior probability distribution of fault model parameters

2.2 分布式滑動分布

由于均勻滑動分布無法得到精細的分布式滑動斷層模型，因此采用非負最小二乘法求解斷層滑動與地表形變的關(guān)系[12]，反演的目標函數(shù)為：

‖W(Gs-d)‖2+k2‖Ls‖2=min

(1)

式中，W為觀測值權(quán)陣，G為格林函數(shù)，s為每個斷層單元上的滑移量，d為地表觀測值，k2為平滑因子，L為拉普拉斯二階有限差分算子。

由于平滑因子對分布式滑動分布反演結(jié)果影響較大，因此固定其他參數(shù)以測試平滑因子的大小，最終確定平滑因子為35。為構(gòu)建斷層面的精細同震破裂模型，對斷層長度和深度分別進行適度擴展，以消除反演結(jié)果的邊緣效應(yīng)，并將斷層面劃分為1 km×1 km的矩形單元。

建立分布式斷層滑動模型是同震反演的關(guān)鍵步驟。由于2022-09-07的升軌數(shù)據(jù)質(zhì)量存在問題，為測試其對反演結(jié)果的影響，對本次地震升降軌數(shù)據(jù)進行以下2種實驗：1)使用2022-08-26～09-07升軌T26和降軌視線向數(shù)據(jù)進行斷層滑動分布反演；2)使用2022-08-26～09-19升軌T26和降軌視線向數(shù)據(jù)進行斷層滑動分布反演?；?種實驗方式分別構(gòu)建分布式滑動分布模型，然后根據(jù)數(shù)據(jù)模擬效果評估反演結(jié)果。

以2022-08-26～09-07升軌T26和降軌T135同震形變場為約束進行滑動分布反演，將斷層走向設(shè)置為NNW-SSE，對應(yīng)的數(shù)據(jù)模擬結(jié)果見圖5。由圖可見，升軌的觀測值和模擬值顯著不同，存在明顯的殘差，殘差主要集中在震中東側(cè)；降軌的觀測值和模擬值吻合較好。升軌數(shù)據(jù)在震中東側(cè)產(chǎn)生明顯殘差的原因可能是該地區(qū)植被茂盛，且地形梯度較大，使得升軌圖像中存在較為明顯的大氣誤差。

圖5 升降軌同震形變場觀測值、模擬值和殘差分布Fig.5 Observed and simulated values and residual distribution of the coseismic deformation of ascending and descending

利用2022-08-26～09-19升軌T26和降軌T135同震形變場作為約束，反演同震滑動分布，數(shù)據(jù)模擬結(jié)果見圖6。由圖可見，升軌平均殘差為1.2 cm，降軌平均殘差為1.8 cm，殘差圖中不存在明顯的系統(tǒng)誤差，表明同震滑動模型能較好地模擬地表InSAR觀測數(shù)據(jù)。

圖6 升降軌同震形變場觀測值、模擬值和殘差分布Fig.6 Observed and simulated values and residual distribution of the coseismic deformation of ascending and descending

圖6中對應(yīng)的二維和三維同震滑動分布模型見圖7。反演結(jié)果顯示，斷層滑動主要以走滑為主，兼少量逆沖分量，此結(jié)果與GCMT、USGS等給出的結(jié)果大致相同。地震破裂深度主要集中在0～17 km，破裂沿走向延伸約55 km，大于0.8 m的較大量級滑動主要發(fā)生在0～7 km深度處，斷層最大滑動量為1.12 m，對應(yīng)深度為1 km，釋放的總地震矩為1.02×1019Nm，相當于MW6.64，此結(jié)果與GCMT計算的結(jié)果大致相同。由表1可見，不同類型觀測資料約束的地震震級和發(fā)震斷層參數(shù)基本一致。余震活動主要發(fā)生在同震破裂區(qū)的中部、深部和邊緣擴展區(qū)，總體上看，余震分布與同震破裂具有較好的空間互補關(guān)系，說明瀘定地區(qū)余震活動是由地震同震破裂觸發(fā)的應(yīng)力加載所致。

圖7 瀘定地震滑動分布模型Fig.7 Slip distribution model of Luding earthquake

3 同震庫侖應(yīng)力變化

地震的破裂滑動會使區(qū)域應(yīng)力場發(fā)生改變，進而影響余震活動，因此余震活動分布可以驗證破裂滑動分布結(jié)果的正確性。分析震中CFS的變化是預(yù)測震中附近發(fā)生地震可能性的重要方法之一。利用反演的斷層滑動分布參數(shù)，運用Coulomb3.3軟件包計算瀘定地震附近地區(qū)的庫侖應(yīng)力變化。

CFS=Δτ+μ′Δσn

(2)

式中，Δτ為剪切應(yīng)力變化，μ′為有效摩擦系數(shù)，Δσn為法向應(yīng)力變化。

將泊松比設(shè)為0.25，有效摩擦系數(shù)設(shè)為0.4，分別計算5 km、10 km、15 km和20 km深度處的同震庫侖應(yīng)力變化，結(jié)果見圖8。由圖可見，隨著深度增加，庫侖應(yīng)力在震中附近及靠南地區(qū)發(fā)生顯著變化，5～15 km深度范圍內(nèi)大部分余震分布在庫侖應(yīng)力加載區(qū)域，符合地震應(yīng)力傳輸和同震庫侖應(yīng)力觸發(fā)機制；小部分余震活動分布在應(yīng)力卸載區(qū)域，表明此次地震未能釋放斷層面積累的全部應(yīng)力。余震分布趨勢為NNW-SSE向，與斷層滑動分布走向大致相同。瀘定地震發(fā)生后，庫侖應(yīng)力的增加主要集中在破裂斷層的西北段、東南段以及西南段，位于鮮水河斷裂帶、安寧河斷裂帶北段、玉龍希斷裂中北段。上述地區(qū)未來發(fā)生地震的危險性很大，尤其需要注意安寧河斷裂帶北段，該斷裂帶在1480年和1536年發(fā)生7.5級地震、在1952年發(fā)生6.75級地震，是需要重點關(guān)注的地震空區(qū)[13]。

紅色五角星為震中位置，白色實線為斷層位置，白色虛線為斷層俯視圖，黑色實線為斷層，黃色圓圈為不同深度的余震分布圖8 瀘定地震引起的CFS變化Fig.8 Variation of CFS caused by Luding earthquake

4 震源區(qū)孕震特征

瀘定地震發(fā)生在鮮水河斷裂帶東南端，毗鄰鮮水河斷裂帶、龍門山斷裂帶、安寧河斷裂帶的交會區(qū)。分析震源區(qū)震前長期應(yīng)變特征，可為識別強震孕育階段以及預(yù)測強震提供重要參考。本文利用球面多尺度小波變換方法[14]，以GPS速度場為約束，建立瀘定地震震源區(qū)地殼應(yīng)變場模型。GPS速度場數(shù)據(jù)來源于Wang等[15]給出的鮮水河地區(qū)近20 a的GPS觀測結(jié)果。球面多尺度小波分解時，球面網(wǎng)格的最大分辨率(qmax)對結(jié)果的影響較大，根據(jù)數(shù)據(jù)的分布密度，經(jīng)過多次實驗，取qmax=7。此時，最大網(wǎng)格尺度為87.2 km，分辨率低于該閾值的區(qū)域設(shè)定為不確定區(qū)域，在計算應(yīng)變參數(shù)時，忽略不確定區(qū)域的結(jié)果。圖9(a)為鮮水河斷裂及鄰區(qū)主應(yīng)變率和最大剪應(yīng)變率分布?？梢钥闯?，瀘定地震震源區(qū)以NE-SW向的擠壓應(yīng)變?yōu)橹?，主壓?yīng)變率量值約為(20～30)×10-9/a，明顯小于鮮水河斷裂中段和北段，表明鮮水河斷裂擠壓變形從西北段到東南段逐漸減弱，此結(jié)果與藏東高原整體的東向擠出背景基本一致。最大剪切應(yīng)變率基本沿鮮水河-安寧河斷裂帶展布，并且在鮮水河斷裂帶、龍門山斷裂帶和麗江-小金河斷裂帶的交會區(qū)附近存在一個應(yīng)變率高值過渡區(qū)，該過渡區(qū)的位置與瀘定地震的位置基本一致。圖9(b)為震源區(qū)面膨脹率分布。由圖可見，此次瀘定地震震源北側(cè)和東側(cè)以擠壓應(yīng)變?yōu)橹?，西南?cè)以拉張應(yīng)變?yōu)橹?，震中所在的位置正好位于面?yīng)變率從擠壓應(yīng)變到拉張應(yīng)變的轉(zhuǎn)換區(qū)域，該應(yīng)變轉(zhuǎn)換區(qū)的存在可能與多個不同活動塊體在該地區(qū)交會并產(chǎn)生的相互作用有關(guān)。

圖9 瀘定地震震源區(qū)最大剪切應(yīng)變率和面膨脹率分布Fig.9 Distribution of maximum shear strain rate and surface dilatation rate in the focal area of Luding earthquake

5 結(jié) 語

1)瀘定地震升降軌InSAR同震形變場中，升軌最大視線向形變量為-15.0 cm，降軌最大視線向形變量為12.7 cm。

2)基于視線向同震形變場反演的滑動分布模型顯示，發(fā)震斷層以左旋走滑為主，兼少量逆沖，斷層走向沿NNW-SSE，約為167°，沿走向破裂約55 km，傾角約為74°，斷裂主要集中在0～17 km深度處，最大滑動量約為1.12 m，對應(yīng)深度為1 km。瀘定地震釋放的總地震矩為1.02×1019Nm，對應(yīng)矩震級為MW6.64。

3)同震庫侖應(yīng)力結(jié)果顯示，瀘定地震使鮮水河斷裂帶南東段、安寧河斷裂帶北段、玉龍希斷裂中北段處于應(yīng)力加載狀態(tài)，其中安寧河斷裂帶北段庫侖應(yīng)力加載狀態(tài)明顯，未來發(fā)生地震的風(fēng)險較大，需要開展持續(xù)監(jiān)測和危險性評估。

4)震間應(yīng)變場表明，此次瀘定地震震源區(qū)位于拉張應(yīng)變和擠壓應(yīng)變的轉(zhuǎn)換區(qū)域，該應(yīng)變轉(zhuǎn)換區(qū)的存在可能與不同活動塊體在該地區(qū)的交會有關(guān)。