瀘定M6.8地震同震重力與形變效應模擬研究
——閉鎖剪力孕震模式的新證據

金雯 談洪波 申重陽 申宇彤

1)中國地震局地震研究所,武漢 430071

2)湖北省地震局,武漢 430071

0 引言

據中國地震臺網測定,北京時間2022年9月5日12時52分四川甘孜藏族自治州瀘定縣(29.59°N,102.08°E)發生M6.8地震,此次地震是繼漾濞M6.4、瑪多M7.4、蘆山M6.1、馬爾康M6.0地震之后,近2年來在青藏高原中東部地區發生的第5次M≥6.0地震(李傳友等,2022),也是近40年來四川地區發生的M6.0～6.9地震中人員震亡最多的一次(許娟等,2023)。瀘定M6.8地震位于青藏高原東南緣鮮水河斷裂帶南東段的磨西斷裂帶附近,系鮮水河、安寧河和龍門山斷裂的交匯處(圖1)(單新建等,2023)。鮮水河斷裂帶是巴顏喀拉塊體與川滇菱形塊體的邊界,具有左旋走滑活動特征(郭長寶等,2015; 孫東等,2023)。研究表明,本次地震的發震斷層磨西段是全新世強烈活動的斷裂,其位于貢嘎山東側、大渡河西側,斷裂沿線地質條件和地貌地形較為復雜,歷史地震頻發(熊探宇等,2010; 李彩虹等,2022)。

注: 時間為1900—2022年,紅色五星表示瀘定M6.8地震,下同。

關于瀘定M6.8地震的破裂模型,已有一些地震學研究結果,如中國地震局地震預測研究所給出的結果[注]https://www.ief.ac.cn/kydts/info/2022/69528.html.顯示,瀘定地震為左旋走滑剪切型破裂,破裂面呈NW-SE走向,同時向南北兩側擴展,破裂方向與鮮水河斷裂帶南段磨西斷裂延伸方向基本一致,最大滑動量達1.84m; 中國地震局地球物理研究所[注]https://www.cea-igp.ac.cn/kydt/279423.html.利用全波形數據的迭代反褶積與疊加方法來反演地震破裂過程,得到斷層破裂模型,地震矩震級為6.7級; 王衛民等[注]http://www.itpcas.cas.cn/new＿kycg/new＿kyjz/202209/t20220906＿6509485.html.通過遠場體波波形數據獲得了該地震的震源破裂過程模型,得出本次地震為高傾角走滑型事件,震源深度約為13km,最大滑動量為1.92m。這些地震學結果為更深入地理解瀘定地震機制和理論模擬研究提供了參考。

根據震前重力場變化圖像,中國地震局重力學科組對此次地震進行了較好的預測,瀘定地震正好位于重力學科2021年底劃定的年度重點危險區內[注]湖北省地震局. 2021. 2022年度重力學科地震趨勢會商報告.,其主要判據之一為2019年9月—2020年9月重力測量結果(圖2):顯著正變化呈現NE向條帶狀,正負最大差異達 90×10-8m/s2; 震前重力變化大致以石棉—康定為中心呈現四象限分布特征。這種震前特征是否與同震破裂機制存在聯系,是否可用申重陽等(2011)提出的閉鎖剪力模式來解釋,是值得研究的問題。為此,本文基于彈性半無限空間矩形位錯理論,采用張勇②和王衛民③發布的由遠場地震波反演的瀘定M6.8地震發震斷層模型,從理論上模擬了此次地震引起的地表同震重力和形變效應,并與實際觀測資料進行對比分析,以期理解瀘定地震的發震和孕育機制,為地震預測和地震危險性評估提供實例參考。

圖2 瀘定地震前實測重力場變化

1 同震位錯引起的地表重力和形變變化的理論公式

自Steketee(1958)將位錯理論引入地震學以來,很多學者研究了半無限空間均勻介質地球模型的同震變形問題。Okada(1985)總結并整理了前人的工作,給出一套完整簡潔實用的同震變形計算公式,適用于計算任何剪切與張性斷層引起的位移、應變和傾斜變形; Okubo(1991、1992)研究了半無限空間介質內剪切和張裂斷層的重力變化問題,利用類似于Okada(1985)的方法導出點源和有限斷層的同震重力位和重力變化的解析表達式; Sun(1992)和Sun等(1993)基于地球曲率和層狀結構,提出球形地球模型位錯理論; Fu(2007)研究了地球內部橫向不均勻結構的影響,對球形位錯理論進行了補充?？紤]到地震破裂的局域性和地震孕育過程的能量積累主要以彈性應力應變形式進行,地震位錯引起的重力及形變變化往往具有以發震斷層為中心、隨距離增大而較快衰減的特性,因此本文選擇彈性平面位錯理論進行研究。

1.1 單條斷層錯動

Steketee(1958)指出,在各向同性的彈性半無限空間介質中,一個沿斷層面的位錯Σ在空間中任意一點(x1,x2,x3)產生的位移場為

(1)

在O-x1x2x3(左旋)直角坐標系下,x1x2為水平地面,x3垂直向下; 設矩形斷層(斷層下盤)的長和寬分別為L和W,δ為斷層面傾角,U1、U2、U3分別對應于任意位錯的走滑、傾滑和引張位錯分量,d為斷層底界深度(圖3),在自由地表某點(x1,x2,0)引起的重力變化Δg和地表變形Δu可以分別表示為(Okubo,1992)

Δg(x1,x2)={ρG[U1Sg(ξ,η)+U2Dg(ξ,η)+U3Tg(ξ,η)]+

ΔρGU3Cg(ξ,η)}‖-βΔh(x1,x2)

(2)

(3)

其中,G為萬有引力常數;ρ為介質密度; Δρ為張裂紋內密度與介質密度之差;Sg(ξ,η)、Dg(ξ,η)、Tg(ξ,η)、Cg(ξ,η)和S(ξ,η)、D(ξ,η)、T(ξ,η)為系數(Okubo,1992;Okada,1985); 自由空氣重力梯度β=0.3086×10-5/s2; Δh為地表高程變化(Okada,1985); 雙豎線符號“‖”可表示為(Okubo,1992)

f(ξ,η)‖=f(x1,p)-f(x1,p-W)-f(x1-L,p)+f(x1-L,p-W)

(4)

其中,p與斷層傾角δ關系式為

p=x2cosδ+dsinδ

(5)

1.2 多斷層錯動

設一地表坐標系O-XYZ,X軸指向正東,Y軸指向正北,Z軸垂直向下,S(sx,sy,0)為斷層頂部地表跡線(或投影線)中點(圖4)。設斷層方位角為α,則由N條斷層位錯引起的地表任意一點(x,y,0)的地表形變u和重力變化ΔG可表示為多條單斷層模型的疊加,即

(6)

圖4 斷層運動地表坐標系示意圖

(7)

由公式(6)、(7)可以得到,斷層尺寸(長、寬)、斷層產狀(傾向、傾角、走向)、斷層位置(埋深、方位、中心投影坐標)以及錯動性質(錯動量的大小、錯動方式)是影響地表形變效應和重力變化的關鍵因素。

2 瀘定M6.8地震同震破裂模型

選取張勇②給出的同震破裂模型(簡稱張勇模型,下同)和王衛民③給出的同震破裂模型(簡稱王衛民模型,下同)分別進行計算分析(圖5)。2個破裂模型包含多個子斷層,每個子斷層的參數包括中心位置(經度、緯度)、中心深度、滑動量、滑動角、方位角和傾角。

圖5 瀘定M6.8地震滑動分布模型

張勇模型斷層面上的靜態滑動分布如圖5(a)所示,結果顯示震中29.59°N、102.08°E,震源深度16km,越靠近最佳矩心震源,滑動量越大,最大靜態位錯約為1.6m。其發震斷層的參數為:走向163°、傾角80°、滑動角8°,該斷層分別沿走向和傾向方向均勻地分成13×7塊子斷層,每個子斷層的尺度為3km×3km,并給出了每個子斷層的滑動角和滑動量。

王衛民模型斷層面上的靜態滑動分布如圖5(b)所示,結果顯示震中29.59°N、102.08°E,震源深度13km,最大靜態位錯約2m。發震斷層參數為:走向166°、傾角78°、滑動角約0°,該斷層分別沿走向和傾向方向均勻地分成17×9塊子斷層,每個子斷層的尺度為3km×2km。

這2種模型均是基于瀘定地震的遠場地震波形數據反演得到的,盡管其發震斷層均表現為左旋走滑特征,但在斷層破裂長度、寬度、震源深度、滑移量等方面具有明顯差異。

3 地表同震重力變化模擬結果

基于上述理論公式和2種斷層破裂模型,模擬計算瀘定M6.8地震引起的地表同震重力變化。根據同震效應影響的范圍,選取計算區域為28°N～31°N、100.5°E～103.5°E,單元網格尺寸為0.025°×0.025°。

根據張勇模型的模擬結果(圖6(a)),地表重力變化圖像呈現明顯的四象限分布特征,以發震斷層為界,近震區(距發震斷層小于50km)以東區域呈現正的重力變化,最大正變化量約7×10-8m/s2,震中以西區域呈現負的重力變化,最大負變化量約-14 ×10-8m/s2; 遠震區(距發震斷層大于50km)的北東區域和西南區域呈現正重力變化,而西北區域和東南區域呈負重力變化。重力變化峰值出現在鮮水河斷裂南東端和安寧河斷裂北端交匯處。

圖6 瀘定M6.8地震斷層位錯引起的同震重力變化

根據王為民模型的模擬結果(圖6(b)),其遠震區的重力變化與張勇模型模擬結果基本特征一致; 近震區存在差異,主要體現在以發震斷層為界,震中EW向呈現正向變化,最大正變化量約6×10-8m/s2,SN向則呈現負向變化,最大負變化量約-8 ×10-8m/s2。

從圖6可以看出,2種模型計算的同震重力變化具有共性和差異性。共性主要表現在遠場,2種模型的模擬結果均具有相似的四象限特征,顯示出斷層源以左旋為主的走滑特征。差異性主要表現在震中近場細節,張勇模型結果顯示震中以東為正變化,以西方向表現出負向變化; 王衛民模型結果展示了震中SN向呈現負向變化,EW向則呈現正向變化。這種近場差異可以歸因于2種斷層模型的細節差異,首先,可能是因為張勇模型聚焦在單一破裂集中區,而王衛民模型涵蓋了2個破裂集中區; 其次,張勇模型主要呈現走滑運動,并帶有逆沖分量(滑動角為8°),而王衛民模型則主要表現為純走滑特征(滑動角為0°)。圖6 震中近場主要呈現局部“高頻”變化,目前還難以用重力觀測等手段予以驗證。同震重力變化峰值位置均出現在鮮水河斷裂南東端和安寧河斷裂交匯處; 同震造成震中東北和西南區域重力值繼續增大,表明這些區域受到的擠壓效應在地震過程中不斷增強,處于物質累積狀態。此外,發現這種擠壓效應的增強與震前重力正變化區域的位置高度重合,這種疊加效應可能會使相關區域地震危險性增強。

4 地表同震位移模擬結果

模擬計算瀘定M6.8地震引起的區域同震位移,結果如圖7所示。

圖7 瀘定M6.8地震引起的同震地表變形

由水平經向位移等值線(圖7(a))可以看出,水平經向位移場具有北正南負的四象限特征。以發震斷層為界,西北和東北區域為正變化,西南和東南區域為負變化。在近震區,該位移場的幅值衰減速度較快(張勇模型結果:0.3～1.3cm,-0.1～-1.1cm; 王衛民模型結果:0.4～1.6cm,-0.2～-1cm),而在遠震區位移場的幅值衰減速度則相對緩慢。震中WN向出現正峰值,張勇模型模擬結果約為1.3cm,王衛民模型模擬結果約為1.6cm; 震中ES向則出現負峰值,張勇模型模擬結果約為-1.1cm,王衛民模型模擬結果約為-1cm。在震中附近的地區,漢源縣表現出西向位移(張勇模型結果:1cm; 王衛民模型結果:0.9cm),石棉縣也表現出西向位移(張勇模型結果:0.5cm; 王衛民模型結果:0.5cm),而瀘定縣則出現東向位移(張勇模型結果:0.4cm; 王衛民模型結果:0.5cm)。

由水平緯向位移等值線(圖7(b))可以看出,水平緯向位移場具有東正西負的四象限特征。以發震破裂斷層為界,北東、南東區域為正變化,而西南、西北區域則為負變化。近震區的幅值也顯示出快速衰減的特點(張勇模型結果:0.1～0.9cm,-0.3～-1.7cm; 王衛民模型結果:0.2～1cm,-0.4～-1.6cm)。震中EN向出現正峰值,張勇模型結果約為0.9cm,王衛民模型結果約為1cm; 而在WS向則出現負峰值,張勇模型結果約為-1.7cm,王衛民模型結果約為-1.6cm。同時,漢源縣出現北向位移(張勇模型結果:0.4cm; 王衛民模型結果:0.5cm),石棉縣也出現北向位移(張勇模型結果:0.3cm; 王衛民模型結果:0.3cm),而瀘定縣同樣出現北向位移(張勇模型結果:0.8cm; 王衛民模型結果:0.9cm)。

由垂直位移等值線(圖7(c))可以看出,垂直位移場具有以震中為中心的正負四象限對稱分布特征,垂直位移主要分布在發震斷層的近震區,東北、西南區域為正變化,西北、東南區域為負變化,近場幅值變化較大(張勇模型結果:0.2～1.8cm,-0.2～-1cm; 王衛民模型結果:0.3～1.5cm,-0.1～-1.5cm)。其中,漢源縣出現沉降位移(張勇模型結果:0.3cm; 王衛民模型結果:0.3cm),石棉縣也出現沉降位移(張勇模型結果:0.2cm; 王衛民模型結果:0.2cm),而瀘定縣出現隆起位移(張勇模型結果:0.3cm; 王衛民模型結果:0.5cm)。

水平位移以及垂直位移圖像的對稱性主要體現了同震位錯的走滑型特征,不同震源模型造成了近震區和遠震區移幅值的差異。地表同震重力變化與垂直位移圖像在遠場區域呈現正相關,擠壓隆起區重力值增大,拉張凹陷區重力值減少,說明該區域同震重力變化主要受斷層位錯引起的地殼密度變化的影響; 地表同震重力變化與垂直位移圖像在近場區域呈現負相關,表明地表垂直位移是該區重力變化的主要因素。

5 討論與結論

本文基于Okubo平面矩形彈性位錯理論,采用已有的由地震波反演獲得的同震破裂模型,模擬研究了瀘定M6.8地震產生的地表同震重力變化、垂直位移和水平位移。結果表明:

(1)瀘定M6.8地震同震重力變化和形變分布形態與走滑型斷層位錯引起的地表重力與位移變化的理論結果(談洪波等,2009)具有較好的一致性,表明瀘定地震發震斷層主要呈現左旋走滑特征。

(2)瀘定M6.8地震孕震過程符合閉鎖剪力孕震模式。地震孕育發生過程會伴隨重力場的變化(申重陽等,2003),震前重力變化可反映孕震現象(申重陽等,2009、2011)。申重陽等(2011、2018)根據2009年姚安M6.0地震前的典型重力變化與震源破裂機制解的比較分析,提出了地震孕育的閉鎖剪力模式:孕震區在外圍動力作用下應先存在一定大小、優勢方向明顯的剪應力以及其作用下形成與之關聯的剪應變,由于該優勢剪切力及其剪應變未超越巖體破裂閾值,使震前孕震區處于“閉鎖”或相對平衡的狀態。通過對大量震例的總結研究(申重陽等,2009、2010、2011、2020; 祝意青等,2020),認為強震前孕震區重力場變化具有典型的梯度帶或四象限分布特征。瀘定M6.8地震前1年,震中區存在明顯的重力場變化正負四象限分布特征:以瀘定為中心,北東區域和西南區域受力擠壓呈正重力變化,西北區域和東南區域為拉張區域呈負重力變化,且震前重力變化和同震重力變化均具有類似的四象限特征,說明震前孕震源與同震破裂具有一定聯系,這種聯系與姚安M6.0地震前重力變化四象限分布特征類似,可能說明引起地震破裂的剪切力在震前就一直存在,即可用閉鎖剪力模式來解釋,這為孕震源的研究提供了新的證據。重力的震前觀測結果與同震模擬結果存在量值差異,可能反映出孕震源與同震破裂源之間的差異。

(3)地表實測GNSS結果(單新建等,2023)與本文理論模擬結果的對比,如表1所示。模擬結果與GNSS實測結果顯示的變形特性一致,位移整體呈現四象限特征分布:震中西南側具有遠離震中的SW向同震形變,震中東南側具有指向震中的NW向同震形變,震中東北側具有遠離震中的NW向的同震形變,震中西北側具有指向震中的SE向的同震形變。整體來說,實測結果與理論模擬結果在遠場區域吻合較好,而在近場區域差異稍大,這可能與本文所采用的發震斷層模型是基于遠場地震波數據反演獲得,其對近場形變控制能力相對較弱有關。雖然本文選取的斷層模型比較簡單,理論模擬結果與實測結果存在一定的偏差,但其基本特征是一致的,本文結果可為瀘定M6.8地震機理研究提供參考。

表1 模擬結果與GNSS實測結果

(4)瀘定M6.8地震對貢嘎山的隆升具有抑制作用。瀘定地震發生于磨西斷裂帶,該斷裂帶位于貢嘎山東側。同時,貢嘎山位于青藏高原東南緣,其主峰距震中約56.7km,貢嘎山是青藏高原隆升作用的重要組成區域之一(劉方斌等,2021)。有研究認為貢嘎山花崗巖體是鮮水河斷裂至安寧河斷裂間擠壓彎曲段吸收、轉換川滇地塊SE向水平運動導致局部快速隆升的產物,其最南端隆升速率超過(3.3±0.8)mm/a(譚錫斌等,2010)。通過模擬貢嘎山主峰(29°35′44″N, 101°52′44″E)的同震水平形變和垂直形變(表2),2個模型模擬結果均顯示貢嘎山有向東的位移量(張勇模型結果:0.73cm; 王衛民模型結果: 0.47cm)和向南的位移量(張勇模型結果:0.41cm; 王衛民模型結果:0.71cm),垂向位移為負值(張勇模型結果: -0.66cm; 王衛民模型結果: -0.52cm),表明貢嘎山出現沉降。對比沉降量與隆升速率,瀘定地震對貢嘎山的隆升產生了2年左右的延遲效應,這表明該地震對貢嘎山隆升產生了一定的抑制作用。因此,此次瀘定地震可能對貢嘎山及其周邊地區的地表變形和地殼運動產生一定影響。

表2 貢嘎山的同震效應模擬結果

致謝：張勇研究員和王衛民研究員提供了發震斷層模型數據,審稿專家提出了寶貴修改意見,在此一并表示感謝。