孟慶筱 黨學會

摘要：以中國地殼運動觀測網(wǎng)絡(luò)2009—2013年GPS觀測數(shù)據(jù)為邊界條件，使用非連續(xù)接觸有限元技術(shù)構(gòu)建九寨溝地區(qū)二維有限元模型，在不確定性分析的基礎(chǔ)之上，計算區(qū)內(nèi)主要斷裂帶現(xiàn)今運動速率。研究結(jié)果表明：在巴顏喀拉塊體整體近于NE向的推擠過程中，九寨溝地區(qū)的塔藏斷裂、虎牙斷裂、樹正斷裂均呈現(xiàn)為較高的左旋走滑兼具擠壓的現(xiàn)今運動特征；岷江斷裂、龍日壩斷裂和龍門山斷裂則呈現(xiàn)為右旋走滑兼有擠壓的運動特性。結(jié)合區(qū)域主應變率計算結(jié)果，發(fā)現(xiàn)九寨溝地區(qū)仍然具有較高的應變積累背景。樹正斷裂作為2017年8月8日九寨溝M7.0地震的發(fā)震斷層，其現(xiàn)今左旋滑動速率為3.0 mm/a，與東昆侖斷裂帶瑪沁—瑪曲段附近的左旋走滑速率4.1mm/a基本匹配，說明該斷裂可能是東昆侖斷裂帶東端分支斷裂之一，而東昆侖斷裂與虎牙斷裂之間的歷史地震空區(qū)可能已被九寨溝地震事件貫通。

關(guān)鍵詞：九寨溝地震；GPS；有限元；接觸單元；樹正斷裂

中圖分類號：P315.725 文獻標識碼：A 文章編號：1000-0666（2018）03-0390-08

0 引言

2017年8月8日四川省阿壩州九寨溝縣發(fā)生M7.0地震，震中位于（33.22°N、103.83°E），震源深度約20 km。該地震是繼2008年汶川M8.0地震和2013年蘆山M7.0地震之后發(fā)生在巴彥喀拉塊體東部的又一強震事件。巴彥喀拉塊體在印度—歐亞板塊碰撞及高原側(cè)向擴展作用下呈現(xiàn)為現(xiàn)今整體NE向的運動特征，其東部受到華南地塊的阻擋，地殼呈近EW向擠壓縮短變形，造成塊體東緣密集分布的活動構(gòu)造和頻繁的強震活動（徐錫偉等，2017；易桂喜等，2017；屈勇，朱航，2017）。隨著區(qū)域構(gòu)造應力的長期積累、集中和加強，應變能最終可能在特定構(gòu)造部位處以中強地震事件的形式突然釋放出來。

中國地殼運動觀測網(wǎng)絡(luò)GPS觀測結(jié)果可以提供高精度、大范圍和準實時的地殼運動定量數(shù)據(jù)，利用復測獲得的地殼運動速度場，通過研究斷裂帶運動學特征，可以為認識塊體內(nèi)部動力學驅(qū)動機制提供重要依據(jù)，對于未來區(qū)域地震危險性判斷具有重要意義（陳長云等，2012）。因此，本文嘗試以中國地殼運動觀測網(wǎng)絡(luò)2009—2013年GPS觀測數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)，在忽略深度對斷層活動差異性影響的基本假設(shè)下，使用非連續(xù)接觸技術(shù)構(gòu)建九寨溝地區(qū)二維有限元模型。在對斷裂帶摩擦系數(shù)等重要參數(shù)的不確定性分析基礎(chǔ)之上，將有限元模擬結(jié)果作為GPS速度場強約束下復雜斷裂帶滑動速率的有效近似，分析九寨溝地區(qū)的強震危險性。

1 模型與方法

1.1 九寨溝地區(qū)二維非連續(xù)接觸有限元模型

考慮到研究區(qū)內(nèi)斷裂構(gòu)造十分發(fā)育，為了提高模型的可分析性，在保證模型可以客觀反映基本問題的前提下對次級斷裂帶進行簡化，納入模型的斷層共11條，如圖1所示。為了清晰地展示九寨溝地區(qū)有限元模型，圖2給出了格網(wǎng)空間步長為5 km的示意圖，而實際后續(xù)研究中筆者使用更加精細的有限元模型，空間步長取0.5 km。在ANSYS有限元環(huán)境下，本文模型中構(gòu)造塊體和斷裂過程帶采用彈性連續(xù)本構(gòu)關(guān)系，使用平面應變單元對研究區(qū)進行單元剖分。

為了獲得九寨溝地區(qū)斷裂帶的現(xiàn)今形變特征，使用接觸單元構(gòu)建斷裂帶。在對非連續(xù)接觸的摩擦方面進行控制時，為了處理斷裂帶兩盤的接觸和相對滑動問題，采用基于接觸力的對稱罰函數(shù)法，并將接觸摩擦力考慮為庫侖摩擦，使用ANSYS中Conta175和Targe169接觸單元進行建模（安關(guān)峰等，2001）。

由于斷裂帶內(nèi)介質(zhì)常常發(fā)育有糜棱巖、斷層泥等軟弱介質(zhì)，其物理力學性質(zhì)與周圍地質(zhì)體具有較大差異，為了更加真實地反映斷裂帶附近破碎巖體較圍巖強度更低的實際情況，本文通過介質(zhì)物性參數(shù)弱化的方法在模型中構(gòu)建了斷裂過程帶，如圖3所示（Vermilye，Scholz，1998）。

1.2 介質(zhì)物理力學參數(shù)

為了盡可能降低有限元模型的參數(shù)不確定性，筆者使用中國地震局工程力學研究所國家強震動臺網(wǎng)中心提供的2017年8月8日九寨溝M7.0地震觀測波形計算了研究區(qū)域內(nèi)的P波與S波波速（Kenichi et al，2010），用于有限元彈性介質(zhì)的物理力學參數(shù)設(shè)置。圖4給出了P波與S波波速計算結(jié)果及其所用強震觀測時程。

根據(jù)圖4，將研究區(qū)域基巖塊體的P波和S波波速取為6.1 km/s和3.5 km/s。參考徐晶等（2017）給出的巴彥喀拉塊體東部分層介質(zhì)模型參數(shù)，將介質(zhì)密度取作2.8×103 kg/m3。根據(jù)Kenichi等（2010）對2000年日本鳥取MW6.6地震事件的研究結(jié)果，破裂過程帶中存在的大量裂隙將會導致S波波速顯著降低，可以將破裂過程帶的S波波速設(shè)置為基巖塊體區(qū)域的構(gòu)建破裂過程帶。因此，本文將九寨溝地區(qū)斷裂帶的破裂過程區(qū)的S波波速設(shè)置為基巖塊體的0.5倍，即1.75 km/s。

斷裂過程帶寬度設(shè)置方面，Vermilye和Scholz（1998）將斷裂帶最大裂隙密度作為斷裂過程帶寬度的閾值標記，給出了破裂過程帶寬度與斷層長度之間的統(tǒng)計分布，如圖5所示。據(jù)此，本文將九寨溝地區(qū)斷層破裂過程帶寬度設(shè)置為斷裂帶長度的1%。進一步參考雷顯權(quán)等（2011）的研究，將斷裂帶接觸單元法向接觸剛度取為75 GPa/m，切向接觸剛度取為3.5 GPa/m。根據(jù)劉曉紅等（1986）使用雙剪法摩擦實驗得到的斷層泥摩擦系數(shù)平均值，將區(qū)內(nèi)各斷裂的摩擦系數(shù)統(tǒng)一為0.7。

1.3 邊界條件

本文采用中國地殼運動觀測網(wǎng)絡(luò)2009—2013年GPS觀測數(shù)據(jù)，利用GAMIT/GLOBK和QOCA軟件進行處理，得到了相應的速度場結(jié)果，其精度約為1 mm/a。為了盡可能剔除2008年汶川M8.0地震同震形變的影響，解算過程中通過插值獲得各個GPS流動站處汶川M8.0地震的同震位移場，進而對同震位移大于3 mm/a的測站加以改正（國家重大科學工程“中國地殼運動觀測網(wǎng)絡(luò)”項目組，2008；江在森等，2009）。通過坐標旋轉(zhuǎn)，得到區(qū)域無旋轉(zhuǎn)坐標框架下GPS速度場（陳長云等，2012），如圖6所示。將靠近模型矩形區(qū)域邊界處的GPS測點的實際觀測速度作為距離最近處節(jié)點的位移約束，通過雙線性插值方法得到關(guān)鍵控制點之間其他節(jié)點的邊界條件。

本文研究區(qū)域距離汶川地震震中較近，雖然在解算過程中對汶川地震同震形變進行了剔除，但是考慮到震后影響的復雜性，未能找到有效模型對汶川地震震后效應進行改正。因此本文所使用的GPS觀測數(shù)據(jù)既包含了汶川地震后的變形，也包含了一部分汶川地震事件中尚未扣除干凈的同震位移，這可能會導致最終斷裂帶滑動速率計算結(jié)果偏大（武艷強等，2012）。

2 計算結(jié)果與討論

2.1 不確定性分析

為了分析參數(shù)不確定性對有限元模擬結(jié)果的影響，圖7給出了研究區(qū)內(nèi)21個GPS測站有限元模擬結(jié)果和GPS觀測結(jié)果在EW分量和SN分量上的誤差對比情況。由圖可見，邊界條件附近兩者一致性較好，而區(qū)域中心位置附近兩者的差異性則相對有所增大，但是模擬值總體上處于GPS觀測結(jié)果解算誤差范圍之內(nèi)，說明當前參數(shù)條件下的模擬結(jié)果具有較好的可靠度，可以用于進一步的分析。

本文嘗試計算了斷裂帶接觸單元摩擦系數(shù)取0.5，0.6，0.7和0.8這4種情況下的地殼運動速度場的模擬結(jié)果，并檢測與實際觀測結(jié)果的一致性。結(jié)果表明，摩擦系數(shù)取0.7時，模擬結(jié)果與觀測結(jié)果的一致性最好，但是摩擦系數(shù)取值對計算結(jié)果的影響作用較小，其原因可能在于本文模擬過程未能進行斷層閉鎖到滑動的動態(tài)模擬。

2.2 區(qū)內(nèi)斷裂運動速度分析

本文利用九寨溝地區(qū)二維有限元模型中斷裂帶兩盤接觸單元的計算得到斷裂現(xiàn)今運動特征與運動速率，主要包括平行于斷層走向的走滑分量（右旋/左旋）Δvτ以及垂直于斷層走向的壓張分量（擠壓/拉張）Δvn。

圖8給出了九寨溝地區(qū)斷裂帶現(xiàn)今走滑速率及壓張速率，其中線寬代表速率大小，斷裂線條越寬則表示斷裂此處運動速率越大。按照圖8的現(xiàn)今運動性質(zhì)計算結(jié)果，研究區(qū)內(nèi)斷裂帶可以分為左旋走滑兼有擠壓、左旋走滑兼有輕微拉張、右旋走滑兼有擠壓3組。其中，東昆侖斷裂、白龍江斷裂、光蓋山—迭山斷裂、哈南—青山灣—稻畦子斷裂、塔藏斷裂、虎牙斷裂、樹正斷裂為左旋走滑運動兼有擠壓的特征，左旋走滑速率分別為4.1，0.8，1.5，2.5，4.0，2.5，3.0 mm/a，垂直于斷裂走向的擠壓速率分別為0.5，2.0，0.5，2.4，2.3，2.8，0.2 mm/a；雪山梁子斷裂為左旋走滑兼有輕微拉張的現(xiàn)今運動特性，其左旋走滑速率為0.8 mm/a，拉張運動速率為0.2 mm/a；岷江斷裂、龍日壩斷裂和龍門山斷裂則呈現(xiàn)為右旋走滑兼有擠壓的運動特性，右旋走滑速率分別為2.1，4.8，1.5 mm/a，垂直于斷裂走向的擠壓速率分別為0.5，1.5，3.1 mm/a（表2）。

本文關(guān)于東昆侖斷裂、白龍江斷裂、光蓋山—迭山斷裂、哈南—青山灣—稻畦子斷裂、塔藏斷裂、虎牙斷裂、龍日壩斷裂、龍門山斷裂無論是現(xiàn)今運動性質(zhì)還是運動速率的分析，均與前人使用最鄰近斷層兩側(cè)的GPS測站位移速率計算方法以及斷裂位錯模型得到研究結(jié)果具有較好的一致性（程佳等，2009；陳長云等，2012；易桂喜等，2017），符合前人研究中巴彥喀拉塊體向東方向運動過程中受到華南塊體阻擋而強烈變形的相關(guān)認識。岷江斷裂的運動特性存在爭議，本文結(jié)果與杜方等（2009）和陳長云等（2012）根據(jù)GPS觀測數(shù)據(jù)得到的岷江斷裂右旋走滑速率分別為2.0，2.5 mm/a一致，而與趙小麟等（1994）利用地質(zhì)、地貌資料得到的左旋走滑運動特征相左?？紤]到地震地質(zhì)資料反映的是斷裂活動的長期運動特征，而GPS觀測數(shù)據(jù)獲得的是斷裂帶現(xiàn)今運動速率，該差異性說明岷江斷裂除具有擠壓特性以外，自形成以來可能經(jīng)歷了左旋與右旋走滑兼有的多期運動（陳長云等，2012；徐錫偉等，2017）。

2.3 樹正斷裂與九寨溝地震

樹正斷裂作為2017年8月8日九寨溝M7.0地震的發(fā)震斷層，是位于塔藏斷裂與岷江斷裂之間的一條NW—SE向無名斷裂，九寨溝地震之前該斷裂無相關(guān)研究，震后被正式命名為樹正斷裂（易桂喜等，2017）。本文計算得到的該斷裂左旋走滑速率為3.0 mm/a，同時伴有擠壓作用。圖9給出了九寨溝地區(qū)區(qū)域水平主應變場的有限元模擬結(jié)果。

由圖9可知，2009—2013年九寨溝地區(qū)的水平主壓應變率約為（2.0～9.0）×10-8/a，水平主張應變率約為（2.5～6.5）×10-8/a，該結(jié)果與陳長云等（2012）利用GPS觀測獲得的結(jié)果基本一致，說明巴彥喀拉塊體整體處于近于NE向的推擠過程中，九寨溝地區(qū)仍然具有較高的應變積累背景。該區(qū)域構(gòu)造背景下，2017年九寨溝M7.0地震的震源機制解為純走滑型，且與余震帶長軸方向一致的節(jié)面走向均值約為152.3°（易桂喜等，2017），與樹正斷裂較高的左旋走滑速率以及微弱的擠壓速率特征一致。

徐錫偉等（2017）認為其成因機制主要在于樹正斷裂北段走向偏西而南段偏北，在巴彥喀拉塊體整體近于向東推擠過程中，NNW向的北段左旋運動轉(zhuǎn)化為近EW向的地殼縮短和岷山垂直隆升作用。同時考慮到塔藏斷裂、虎牙斷裂以及本文得到的樹正斷裂左旋滑動速率3.0 mm/a與東昆侖斷裂帶瑪沁—瑪曲段附近左旋走滑速率4.1 mm/a基本匹配，說明塔藏斷裂、虎牙斷裂以及樹正斷裂可能是東昆侖斷裂帶東端分支斷層之一。與此同時，2017年九寨溝M7.0地震及其余震密集帶恰好填補了虎牙斷裂帶西北段歷史地震空區(qū)，故樹正斷裂可能是虎牙斷裂的NW延伸的隱伏部分，而九寨溝地震的發(fā)生則可能貫通了虎牙斷裂與東昆侖斷裂（季靈運等，2017；徐錫偉等，2017），其區(qū)域地震危險性值得進一步關(guān)注。

3 結(jié)論

本文利用中國地殼運動觀測網(wǎng)絡(luò)2009—2013年GPS觀測數(shù)據(jù)，使用二維非連續(xù)接觸數(shù)值模擬方法計算了九寨溝地區(qū)斷裂構(gòu)造現(xiàn)今運動形變特征，得到以下主要結(jié)論：

（1）東昆侖斷裂、白龍江斷裂、光蓋山—迭山斷裂、哈南—青山灣—稻畦子斷裂、塔藏斷裂、虎牙斷裂、樹正斷裂為左旋走滑兼有擠壓的運動特征，左旋走滑速率分別為4.1，0.8，1.5，2.5，4.0，2.5，3.0 mm/a，垂直于斷裂走向的擠壓速率分別為0.5，2.0，0.5，2.4，2.3，2.8，0.2 mm/a；岷江斷裂、龍日壩斷裂和龍門山斷裂則呈現(xiàn)為右旋走滑兼有擠壓的運動特性，右旋走滑速率分別為2.1，4.8，1.5 mm/a，垂直于斷裂走向的擠壓速率分別為0.5，1.5，3.1 mm/a。結(jié)合區(qū)域主應變率計算結(jié)果，說明九寨溝地區(qū)仍然具有較高的應變積累背景。

（2）樹正斷裂作為2017年8月8日九寨溝M7.0地震的發(fā)震斷層，其現(xiàn)今左旋滑動速率約為3.0 mm/a，與東昆侖斷裂帶瑪沁—瑪曲段附近的左旋走滑速率4.1 mm/a基本匹配，說明該斷裂可能是東昆侖斷裂帶東端分支斷裂之一，可能是虎牙斷裂北西段的隱伏部分，九寨溝地震事件可能已將東昆侖斷裂與虎牙斷裂貫通，其區(qū)域地震危險性值得進一步關(guān)注。

本文在撰寫過程中得到中國地震局第一監(jiān)測中心陳長云博士和占偉博士的幫助，在此向他們表示衷心感謝。

參考文獻：

安關(guān)峰，高大釗，趙炯洋.2001.接觸面彈粘塑性本構(gòu)關(guān)系研究[J].土木工程學報，34（1）：88-91.

陳長云，任金衛(wèi)，孟國杰，等.2012.巴顏喀拉塊體北東緣主要斷裂現(xiàn)今活動性分析[J].大地測量與地球動力學，32（3）：27-30.

陳顒，黃庭芳.2009.巖石物理學[M].合肥：中國科學技術(shù)大學出版社.

程佳，甘衛(wèi)軍，王澤河，等.2009.2001年昆侖山口西MS8.1地震前背景形變場的模擬研究[J].地震地質(zhì)，31（1）：97-111.

杜方，聞學澤，張培震，等.2009.2008年汶川8.0級地震前橫跨龍門山斷裂帶的震間形變[J].地球物理學報，52（11）：2729-2738.

國家重大科學工程“中國地殼運動觀測網(wǎng)絡(luò)”項目組.2008.GPS測定的2008年汶川MS8.0級地震的同震位移場[J].中國科學：地球科學，38（10）：1195-1206.

季靈運，劉傳金，徐晶，等.2017.九寨溝MS7.0地震的Insar觀測及發(fā)震構(gòu)造分析[J].地球物理學報，60（10）：4069-4082.

江在森，方穎，武艷強，等.2009.汶川8.0級地震前區(qū)域地殼運動與變形動態(tài)過程[J].地球物理學報，52（2）：505-518.

雷顯權(quán)，陳運平，趙炯洋.2011.天山現(xiàn)今地殼變形的非連續(xù)接觸模型模擬[J].中南大學學報（自然科學版），42（3）：2754-2762.

劉曉紅，方亞如，蔡戴恩，等.1986.海原斷裂帶斷層泥摩擦特性的研究[J].西北地震學報，8（3）：50-54.

屈勇，朱航.2017.巴顏喀拉塊體東—南邊界強震序列庫侖應力觸發(fā)過程[J].地震研究，40（2）：216-225.

武艷強，江在森，楊國華，等.2012.南北地震帶北段近期地殼變形特征研究[J].武漢大學學報（信息科學版），37（9）：1046-1048.

徐晶，邵志剛，劉靜，等.2017.基于庫侖應力變化分析巴顏喀拉地塊東端的強震相互關(guān)系[J].地球物理學報，60（10）：4056-4068.

徐錫偉，陳桂華，王啟欣，等.2017.九寨溝地震發(fā)震斷層屬性及青藏高原東南緣現(xiàn)今應變狀態(tài)討論[J].地球物理學報，60（10）：4018-4026.

易桂喜，龍鋒，梁明劍，等.2017.2017年8月8日九寨溝M7.0地震及余震震源機制解與發(fā)震構(gòu)造分析[J].地球物理學報，60（10）：4083-4097.

趙小麟，鄧起東，陳社發(fā)，等.1994.岷山隆起的構(gòu)造地貌學研究[J].地震地質(zhì)，16（4）：429-439.

Kenichi F，Yuji K，F(xiàn)usanori M.2010.Influence of fault process zone on ground shaking of inland earthquakes：Verification of Mj=7.3 Western Tottori Prefecture and Mj=7.0 West Off Fukuoka Prefecture earthquakes，southwest Japan[J].Engineering Geology，116（8）：157-165.

Vermilye J M，Scholz C H.1998.The Process Zone：a Microstructural View of Fault Growth[J].Journal of Geophysical Research，103（B6）：12223-12237.

Abstract The interpolated GPS observation data during the 2009-2013 period of the China crustal movement observation network was applied to gain the boundary conditions of 2-D element discontinuous model which was used to simulate the slip rates of the fault system in the Jiuzhaigou area and undertake uncertainty analysis over it.The study turned out that the NE direction movement of Bayan Har block had been hindered by the South China block.The eastern Kunlun fault，Tazang fault，Huya fault and Shuzheng fault showed a left-lateral strike slip and compressional characteristic，when the Minjiang fault，Longriba fault and Longmenshan fault showed a right-lateral strike slip and compressional characteristic.Combined with the calculation results of regional principal strain rate，it shows that the Jiuzhaigou area still has a high background of strain accumulation.As the seismogenic fault of the August 8，2017 Jiuzhaigou M7.0 earthquake，the left-lateral slip rates of Shuzheng fault is about 3 mm/year，which basically matched the slip rate of the eastern Kunlun fault.It is suggested that the fault may be one of the branch of the eastern Kunlun fault zone and the gap between the east Kunlun fault and the Huya fault may have been split by the Jiuzhaigou earthquake.

Keywords：Jiuzhaigou earthquake；GPS finite element；contact node；Shuzheng fault