石富強 朱琳 王瑩 丁曉光 邱玉榮 石軍

陜西省地震局，西安市碑林區(qū)水文巷4號 710068

0 引言

巴顏喀拉塊體東北緣位于南北地震帶中段，區(qū)域內分布有龍門山斷裂、東昆侖斷裂、西秦嶺北緣斷裂等多條大型活動斷裂帶，為現(xiàn)今地殼變形顯著、地震活躍的區(qū)域。僅20世紀以來，該區(qū)域就發(fā)生了1933年疊溪7.5級、1976年松潘平武7.2級、2008年汶川8.0級、2017年九寨溝7.0級等多次7級以上強震。對該區(qū)域的地震危險性跟蹤一直是地震科技工作者關注的焦點，M7專項工作組根據(jù)歷史地震活動、大地測量及地球物理場數(shù)據(jù)對該區(qū)域給出了甘青川交界危險區(qū)、西秦嶺中西段危險區(qū)以及六盤山南-西秦嶺東危險區(qū)等強震高風險區(qū)域。汶川 8.0級地震后，一些學者（Xiong et al，2010；Toda et al，2008；Wan et al，2010；單斌等，2009）基于歷史強震演化的斷層庫侖應力變化的研究結果表明，東昆侖斷裂東段以及西秦嶺北緣斷裂西段的地震危險性值得關注。此外，該區(qū)域還存在歷史地震破裂空段、低b值、小震活動圖像、地震矩加速釋放等異常（陳為濤等，2013；M7專項工作組，2012）以及跨斷層水準等地球物理場觀測數(shù)據(jù)異常①中國地震局第二監(jiān)測中心，2016，2017年度地震趨勢會商報告。。

庫侖應力研究作為研究地震輪回中構造應力場時空變化的重要組成部分，已成為研究斷層應力狀態(tài)及區(qū)域地震危險性的重要手段之一。King等（1994）研究指出，1992年美國Landers 7.3級地震在Big Bear產(chǎn)生的0.2MPa左右的應力加載量加速了Big Bear 6.5級地震的發(fā)生。1999年Izm it 7.4級地震發(fā)生前，Stein等（1997）等基于庫侖應力計算指出該區(qū)域地震危險性較高。而2004年Parkfield 6.0級地震的發(fā)生則受1983年Coalinga 6.5級地震及Nu?ez多個中等地震的應力卸載的影響，出現(xiàn)了約11年的時間延遲（Bakun et al，1985）。此次九寨溝地震就發(fā)生在M7專項工作組（2012）給出的甘青川交界危險區(qū)內，基于庫侖應力的強震作用研究表明，區(qū)域內1976年松潘平武7.2級地震、2008年汶川8.0級地震等都對此次九寨溝地震發(fā)震斷層有較強的庫侖應力加載作用（徐晶等，2017；單斌等，2017）。

有限元數(shù)值模擬技術作為探索構造變形等的虛擬實驗手段，在地震危險性分析中得到了廣泛應用（李玉江等，2013；Moreno et al，2010；Hergert et al，2010；陳連旺等，2008；祝愛玉等，2015）。為進一步分析九寨溝地震對巴顏喀拉塊體東北緣主要活動斷裂的加卸載效應以及該區(qū)域未來可能的強震風險，本文構建了巴顏喀拉塊體東北緣三維粘彈性有限元模型，模擬計算了九寨溝地震對周邊活動斷裂的加卸載效應，分析了震后地震活動及其與九寨溝地震間的可能關系以及未來可能的強震風險。

1 模型構建

1.1 幾何模型及物性參數(shù)

本文在詳細調研地震地質資料的基礎上，通過對所調研的地質資料進行假設構建了巴顏喀拉塊體東北緣有限元模型：模型縱向深100km，地殼分層、速度結構以及密度結構根據(jù)Crust1.0速度結構給出（Laske et al，2013）；模型中包含禮縣-羅家堡斷裂、虎牙斷裂、雪山斷裂、東昆侖斷裂、龍門山斷裂、青川斷裂、青川-平武斷裂、臨潭-宕昌斷裂、白龍江斷裂、塔藏斷裂、龍日壩斷裂、岷江斷裂、西秦嶺北緣斷裂以及基于有限斷層反演給出的此次九寨溝地震發(fā)震斷裂等等；區(qū)域孕震層深度（段星北，1997；張國民等，2002），除九寨溝地震發(fā)震斷層外，假定其他斷層面下邊界深度為20km；除九寨溝地震發(fā)震斷層施加位錯荷載外，其余斷層均為摩擦接觸斷層，斷層空間幾何結構簡化為由斷層地表出露沿著傾角向下延伸的空間斜面；下地殼及地幔粘滯性結構根據(jù)邵志剛等（2008）、徐晶等（2017）的研究給出。地殼單元尺寸控制在10km內，地幔單元尺寸控制在20km內，斷層面上單元尺寸控制在3km內，模型共包含359676個單元、67796個節(jié)點，構建的有限元模型及對應的物性參數(shù)分別見圖1和表1。

1.2 邊界約束及模擬方案

地震是構造應力作用下斷層應力狀態(tài)達到一定程度而斷層破壞形成的自然災害。在數(shù)值模擬中，由于我們并不知道斷層當前絕對應力狀態(tài)以及斷層自身的承載極限，因此，很難做到對地震自發(fā)的模擬。為此，本文首先對模型施加以GPS約束的邊界荷載來模擬初始構造應力場；然后，再施加地震同震位錯，模擬計算給出九寨溝地震后區(qū)域的模擬構造應力場；最后，從震后模擬應力場中扣除模擬初始應力場，給出九寨溝地震的發(fā)生對區(qū)域主要活動斷裂的加卸載效應。其中，九寨溝地震位錯模型參考波形反演結果（http：//www.cea-igp.ac.cn/tpxw275883.html）給出。所有模擬計算過程均在有限元開源代碼 Pylith（https：//geodynamics.org/cig/software/pylith/）平臺下進行。

圖1 本文構建的巴顏喀拉塊體東北緣有限元模型

表1 本文有限元模型力學參數(shù)

1.3 庫侖應力計算

三維有限元模擬軟件并不能直接給出斷層面上的庫侖應力，因此，需要通過計算給出的應力張量來計算斷層面上的正應力及剪應力，進而計算斷層面庫侖應力變化。圖2為斷層面主方向（）與模型坐標系（）相對關系，其中，分別為斷層走向（面向走向，右手在上盤）、傾向（逆沖為正）以及張拉方向（張為正）；分別為地理E、N向以及垂直地表向上的方向。α為斷層走向；β為斷層傾角。將斷層面應力分量在模型坐標系（、）下表示為

其中，T為矩陣轉置符號。

圖2 斷層變形模式示意圖（a）以及斷層面主方向與模型坐標間的關系（b）

由圖2（b）可知，將模型坐標系下的應力張量繞z軸逆時針旋轉90°－α，然后，再繞主軸1逆時針旋轉β，便可得到主方向下的應力張量。則

由式（1）、（2）給出的新的應力分量σ33即為斷層面上的正應力變化 Δσn，σ23、σ13分別為斷層面上沿傾向和走向的剪應力分量。令斷層面優(yōu)勢破裂方向的滑動角為θ，則斷層面優(yōu)勢破裂滑動方向的剪應力變化為

根據(jù)斷層面庫侖應力變化的定義（King et al，1994），給出斷層面優(yōu)勢破裂方向的庫侖應力變化為

其中，μ為斷層面有效摩擦系數(shù)。

2 模擬結果分析

2.1 模擬初始應力場

以1999～2007年GPS觀測速度為約束，通過本文構建的有限元動力學模型計算了巴顏喀拉塊體東北緣的主應力場（圖3）。模擬結果與基于小震震源機制反演的區(qū)域構造應力場（王曉山等，2015）以及基于GPS推算的區(qū)域應變率場（陳長云等，2013）基本一致。區(qū)域主要處于NW-SE拉張、NE-SW西向擠壓的動力環(huán)境。九寨溝地震震中南側至龍門山斷裂，主張應力逐漸減小，主壓方向逐漸偏轉并呈現(xiàn)出近EW向擠壓態(tài)勢。九寨溝地震震中附近主應力方向與 USGS給出的震源機制解（https：//earthquake.usgs.gov/earthquakes/eventpage/us2000a5x1＃moment-tensor）的 P軸（方位角 104°）和 T軸（方位角 204°）基本一致。

2.2 九寨溝地震及其對周邊斷裂的加卸載效應

以同震位錯為荷載模擬給出的九寨溝地震震中附近的同震位移場如圖4所示。由圖4可見，震中最大同震位移約為9cm，九寨溝地震的發(fā)生對東昆侖斷裂、塔藏斷裂、西秦嶺北緣斷裂東段以及龍門山斷裂均為左旋作用，同震位移場的左旋作用有利于左旋活動為主的東昆侖斷裂、塔藏斷裂以及西秦嶺北緣斷裂東段斷層面的應力積累；而對于具有右旋逆沖活動的龍門山斷裂則為卸載作用。同理，對西秦嶺北緣斷裂西段應為卸載影響。為進一步定量分析九寨溝地震對各斷裂的加卸載效應，我們基于本文模擬結果，根據(jù)式（2）、（4）計算了九寨溝地震在這些斷裂斷層面上引起的庫侖應力變化（圖5）。

圖3 本文模擬給出的巴顏喀拉塊體東北緣初始主應力場

圖4 數(shù)值模擬給出的九寨溝地震震中附近同震位移場

由圖5可見，九寨溝地震對塔藏斷裂西段、虎牙斷裂加載效應顯著，達7056Pa；東昆侖斷裂東段、龍日壩斷裂、迭部-白龍江斷裂西端和東端、臨潭-宕昌斷裂東段以及西秦嶺北緣斷裂東段也均受到不同程度的加載作用；塔藏斷裂東段以及岷江斷裂受九寨溝地震的影響，斷層面庫侖應力變化表現(xiàn)為卸載，最大卸載量達73984Pa。位于此次九寨溝地震破裂面正前方的雪山斷裂斷層面庫侖應力變化較為復雜：由圖4模擬同震位移場可知，九寨溝地震對雪山斷裂西段的影響表現(xiàn)為擠壓，對東段的影響表現(xiàn)為拉張，根據(jù)庫侖應力計算公式可知，當斷層面應力變化主要表現(xiàn)為垂直于斷層面的擠壓/拉張時，有效摩擦系數(shù)的取值對庫侖應力變化影響顯著，甚至可以影響其極性變化。Wan等（2010）在分析汶川地震對周邊斷裂影響時也給出同樣的結論。從本文模擬結果看，當有效摩擦系數(shù)較?。é?0.1）時，九寨溝地震對雪山斷裂西段表現(xiàn)為庫侖應力加載，對東段表現(xiàn)為卸載；當有效摩擦系數(shù)較大（μ=0.7）時，對西段表現(xiàn)為卸載作用，而對東段表現(xiàn)為加載（圖5）。

此外，此次九寨溝地震還對汶川地震破裂區(qū)東北端的青川斷裂以及青川-平武斷裂的西南段有一定的庫侖應力加載，最大加載量可達1139Pa。邵志剛等（2010）研究汶川地震的震后影響發(fā)現(xiàn)，青川斷裂庫侖應力加載非常顯著，實際上青川地區(qū)也是汶川地震余震活動非?；钴S的區(qū)域（石軍，2008），但從2013年開始，該區(qū)域地震活動水平明顯減弱。此次九寨溝地震后不到2個月，9月30日青川（32.27°N，105°E）再次發(fā)生5.4級地震。分析認為，此次青川地震可能是由九寨溝地震的觸發(fā)作用所致。

圖5 有限元模擬給出的不同有效摩擦系數(shù)下斷層面庫侖應力變化

3 結論與討論

基于構建的三維有限元模型，以同震位錯為荷載模擬了九寨溝地震的發(fā)生及其對巴顏喀拉塊體東北緣主要活動斷裂的加卸載效應。數(shù)值模擬結果顯示，九寨溝地震的發(fā)生對處于甘青川交界危險區(qū)內的東昆侖斷裂東段至塔藏斷裂西段庫侖應力加載效應顯著，最大可達7056Pa；對六盤山南-西秦嶺東危險區(qū)的西秦嶺北緣斷裂東段也有一定程度的加載效應，約為346Pa。此外，計算結果還顯示，九寨溝地震對龍日壩斷裂、虎牙斷裂、青川-平武斷裂西段、迭部-白龍江斷裂西段和東段、臨潭-宕昌斷裂東段等也均有不同程度的庫侖應力加載效應；對岷江斷裂、塔藏斷裂東段庫侖應力卸載效應顯著，最大卸載量可達73984Pa。

地震是區(qū)域構造應力積累到一定程度且超越斷層面承載極限而引起斷層剪切破壞的動力過程。但實際上我們很難獲取斷層面的絕對應力狀態(tài)和斷層面的承載極限。庫侖應力作為由地震發(fā)生而引起的構造應力場變化在斷層面優(yōu)勢破裂方向的體現(xiàn)，在地震危險性定性分析判斷中有著廣泛的應用。汶川地震后，許多學者通過庫侖應力計算分析了汶川地震對區(qū)域活動斷層的加卸載效應及地震危險性（Toda et al，2008；單斌等，2009；邵志剛等，2010；李玉江等，2013）。隨后的2013年蘆山7.0級地震、岷縣漳縣6.6級地震以及此次九寨溝7.0級地震均發(fā)生于文獻報道的庫侖應力加載區(qū)域。這些強地震的時空叢集活動也表明，巴顏喀拉地塊東段在未來一段時間內可能繼續(xù)保持這種活躍態(tài)勢。從庫侖應力變化的角度來看，此次九寨溝地震對甘青川交界危險區(qū)內的主要活動斷裂（東昆侖斷裂東段和塔藏斷裂）以及位于六盤山南-西秦嶺東危險區(qū)的主要活動斷裂（西秦嶺北緣斷裂東段），均具有不同程度的庫侖應力加載效應。