張 宇 汪建軍,2,3 劉 洋,2,3 熊 維,4
1 武漢大學測繪學院,武漢市珞喻路129號,430079
2 武漢大學地球空間環境與大地測量教育部重點實驗室,武漢市珞喻路129號,430079
3 自然資源部地球物理大地測量重點實驗室,武漢市珞喻路129號,430079
4 中國地震局地震研究所,武漢市洪山側路40號,430071
研究大地震的震后時變機制,對認識地震破裂與發震環境的關系,探索地震斷層的震后運動行為、摩擦性質及巖石圈流變結構等具有重要意義[1]。地震斷層發生同震錯動后,斷層面的震后余滑會導致地殼產生形變,許多地震的早期震后形變可以用余滑模型來解釋。目前用于解釋震后余滑形變的模型包括運動學余滑模型和基于速率-狀態摩擦準則的應力驅動余滑模型。2021-05-22青?,敹喟l生MW7.4地震,震中為98.34°E、34.59°N,震源深度約17 km。余震精定位結果顯示,地震發生在巴顏喀拉塊體北部的NW-W向昆侖山口-江錯斷裂帶[2-3]。布設在瑪多地區的GPS連續運行參考站記錄了地震的震后高時空分辨率形變信息[4-5]。學者們利用大地測量觀測數據,通過運動學余滑模型反演震后余滑[3,6-7],同時利用應力驅動余滑模型解釋震后形變時空特征[7]。本文采用瑪多地震震后160 d的GPS形變時間序列數據,利用基于速率-狀態摩擦準則的應力驅動余滑模型[8],正演此次瑪多地震的震后余滑分布,以分析該地震的震后余滑特征及斷層摩擦屬性。
瑪多地震后,學者們使用不同的觀測數據給出各自的斷層滑動分布模型[5-6],所得滑移大小和地震矩基本相同。震后早期余滑(約160 d)形變由位于發震斷層約50 km內的16個連續GPS (cGPS)站記錄得到,其中13個站點在主震后10 d內建立,另外3個在2021年之前建立(圖1)。本文使用的數據為Xiong等[5]剔除16個cGPS站點數據的震間形變場和季節性變形后得到的震后早期形變時間序列,形變表現出對稱特征,符合高傾角斷層及以走滑為主的滑動機制。值得注意的是,4454站點僅記錄到微弱的震后位移,表明該站點附近幾乎沒有發生余滑,這一結論同樣體現在InSAR觀測和運動學反演結果中[6-7]。

圖1 GPS參考站分布
本文采用Jin等[9]利用InSAR觀測數據反演的瑪多地震同震滑動模型作為同震輸入模型,通過應力驅動余滑模型擬合Xiong等[5]處理得到的震后早期形變時間序列數據。
使用Relax程序[10]模擬瑪多地震震后余滑引起的地表形變。應力驅動余滑模型描述的是地震斷層在受同震應力加載和斷層面摩擦參數控制下的震后滑動行為。根據速率-狀態摩擦準則,以同震破裂在斷層面上產生的剪應力變化為應力初始條件,當給定斷層的摩擦參數(a-b)σ和參考速率V0時,可計算斷層面的震后余滑和地表形變隨時間的演化。與運動學模型相比,應力驅動余滑模型顧及了同震滑動產生的應力變化和斷層面的摩擦參數對震后余滑行為的控制作用。在應力驅動余滑模型中,同震庫侖應力的增加促進了斷層繼續滑動[11]。斷層滑動速率可以表示為:
(1)
式中,V為斷層的震后余滑滑動速率,Δτ為同震剪應力變化,V0為參考滑動速率,(a-b)σ為0.1~10 MPa范圍內的本構參數,σ為斷層面上的有效正應力。
本文取剪切模量為30 GPa,泊松比為0.25,摩擦系數為0.6。采用Barbot等[12]的方法計算應力變化的同震滑動模型,即假設同震滑動量較大的區域是速度弱化區域,不參與震后蠕滑。為使同震滑動和余滑之間的重疊區域允許發生余滑,去除同震滑動峰值30%以下的滑動量,以集中同震滑動分布;然后將其余的滑動量放大,以保持地震矩不變[12-13]。
實際操作中,由于V0和(a-b)σ之間通常存在很強的相關性,不可能獨立確定。盡管改變(a-b)σ的值會導致搜索結果不同,但并不能通過數據擬合來區分不同的(a-b)σ的值[14],這種相關性正是Barbot等[12]認為不能賦予V0任何特定物理意義的部分原因。因此,本文遵循以往使用Relax軟件的方法,將(a-b)σ的值固定為5 MPa,以降低搜索的復雜程度[11]。
震后形變機制主要包括同震應力驅動的震后余滑、下地殼和上地幔粘彈性松弛和地殼淺部的孔隙彈性回彈。為評估除震后余滑外的震后形變機制是否參與了震后早期形變,參考Huang等[15]研究巴顏喀拉塊體內斷層相互作用時使用的下地殼和上地幔粘滯系數(下地殼1018~1019Pa·s,上地幔1020Pa·s),使用Relax軟件計算瑪多地震震后160 d粘彈性松弛引起的累積地表形變,預測水平位移小于6 mm,遠小于GPS觀測值。采用Jin等[9]的同震滑動模型,利用地震前后泊松比的變化計算地表形變場,并對其進行差分,得到平衡狀態的孔隙彈性響應形變,預測水平位移不足1 mm,表明孔隙回彈引起的形變可忽略不計。
假定震后斷層余滑面與同震斷層滑動面完全重合。由于瑪多地震以走滑分量為主,震后GPS形變觀測數據的EW分量整體遠大于NS分量,故擬合GPS觀測數據時僅使用EW分量計算殘差,而用NS分量作定性比較。通過改變V0的值對地表震后EW向形變時間序列進行擬合,得到最優擬合參數為V0=1.5 m/a,RMS=0.76 cm,從而得到GPS站點震后形變時間序列(圖2和3)和余滑分布(圖4(a))。可以看出,震后余滑主要分布在同震斷層顯著滑動的下傾和兩側區域,最大滑移量為0.67 m。余滑模型1所得的地表形變在距離斷層跡線較近的GPS站點處均大于觀測數據,而較遠站點處則明顯小于觀測值,且大部分站點具有比模型1更明顯的衰減趨勢。根據前人研究結果[9,16],瑪多地震發震斷層的摩擦屬性在深度上可能存在差異,因此進一步構建在深部和淺部具有不同摩擦屬性的應力驅動余滑模型(模型2)。

圖2 應力驅動余滑模型模擬的震后地表EW向形變與GPS觀測數據對比

圖3 應力驅動余滑模型模擬的震后地表NS向形變與GPS觀測數據對比

圖4 最佳擬合余滑分布
基于余滑模型1的結果,通過大量正演實驗測試余滑界面位置和摩擦屬性分層深度的影響,最終將震后余滑對應的不同摩擦屬性深度分界設定為13 km(圖5)。余滑模型2以該深度為界,對不同深度的余滑界面設置不同的參考滑動速率V0[17],并參考運動學余滑結果[5-7]約束斷層兩端7 km深度以上的淺部余滑。

圖5 不同摩擦屬性分界深度對應的最小殘差


圖6 余滑模型2對應的擬合殘差
為測試不同同震模型對計算結果的影響,在其余參數不變的情況下,分別使用Wang等[18]和Zhao等[7]的同震滑動模型計算余滑模型2的震后地表LOS向形變,并與本文結果進行對比。由圖7可知,不同模型結果的主要特征比較接近,均沿走向分布在斷層跡線兩側,但地表形變分布的集中程度不同,差異主要在斷層附近。本文得到的模擬結果在形變分布和量級上介于Wang等[18]和Zhao等[7]的結果之間。

圖7 不同同震模型的LOS向形變模擬結果

Xiong等[5]和He等[6]利用運動學余滑模型反演得到的滑動模型均能較好地解釋瑪多地震的震后短期大地測量觀測數據。Xiong等[5]認為,余滑主要發生在同震破裂下側,同震破裂區域基本沒有發生滑移,最大滑動量約為0.3 m,約位于14 km深度處。He等[6]認為,余滑主要發生在同震破裂下傾區域,余滑峰值為0.3 m,約位于20 km深度處。運動學震后余滑特征的差異可能與采用的數據不同有關。Zhao等[7]采用InSAR形變約束的應力驅動余滑結果表明,余滑主要分布在發震斷層的深部,淺部幾乎沒有發生余滑,滑動峰值為0.8 m,位于約10 km深度處,與本文余滑集中區域一致,量級接近余滑模型1的結果。

本文采用應力驅動余滑模型對2021年青?,敹郙W7.4地震的震后余滑進行模擬,較好地解釋了瑪多地震震后160 d的GPS形變觀測數據,獲得了地震的震后余滑特征?,敹嗟卣鸬恼鸷笥嗷饕l生在同震斷層顯著滑動的下傾和兩側區域,集中分布在5~15 km深度(占總量80%),最大滑移量達1.24 m,累積釋放地震矩1.08×1019Nm,相當于MW6.62地震。以13 km為界,淺部和深部在最優擬合下對應的V0分別為6.50 m/a和0.15 m/a,表明斷層摩擦屬性存在深度差異。