石富強 朱琳 李玉江 丁曉光

摘要：總結了“介質弱化帶”和“摩擦接觸”這2類常見斷層模型在地殼動力學數(shù)值模擬研究中的應用研究進展，再利用2者各自的特點，引入復合材料力學中各向異性理論，結合斷層產(chǎn)狀給出了針對不同走向、傾角的斷層各向異性本構方程。數(shù)值模擬試驗對比研究顯示：基于正交各向異性本構關系的斷層模型不但能夠避免斷層垂直方向上的形變場畸變現(xiàn)象，而且能夠有效地刻畫斷層兩盤相對滑動，模擬結果與基于接觸摩擦斷層模型的模擬結果吻合良好。該模型可以在研究區(qū)域構造復雜，接觸計算收斂困難的情況下替代摩擦接觸模型。

關鍵詞：斷層變形；數(shù)值模擬；介質弱化模型；接觸摩擦模型；正交各向異性理論

中圖分類號：P313.5 文獻標識碼：A 文章編號：1000-0666（2018）01-0064-09

0 引言

地殼形變是地殼在力學、熱學等復雜環(huán)境作用下表現(xiàn)出來的綜合響應，研究地殼形變特征及其演化規(guī)律是認識地殼運動及其伴生地質災害的一個重要途徑。隨著計算機科學的日益發(fā)展，數(shù)值模擬技術已經(jīng)成為地殼動力學機理研究的重要手段（Moreno et al，2010；Hergert，Heidbach，2010；Li et al，2014；Zhu，Zhang，2010；He et al，2013；Lu et al，2011；陳連旺等，1999，2008；邵志剛等，2008；李玉江等，2014）。在復雜地質構造作用下，斷層作為巖石圈巖石破裂形成的特殊產(chǎn)物，在現(xiàn)今地殼形變特征以及地震的孕育、發(fā)生中有重要的作用。在探索地震發(fā)生孕育的過程中，國內外科學家從不同角度闡述了斷層物理機制（Sib-son，1977；馬勝利，本利彥，1995）。但由于地質構造運動的復雜性，現(xiàn)在還沒有一種行之有效的辦法可以模擬真實地殼環(huán)境中斷層的形成、發(fā)展變化。因此通過數(shù)值模擬解釋現(xiàn)今地殼形變特征時，往往需要結合野外考察、深地震反射、地震精定位、震源機制解等多學科方法給出的斷層產(chǎn)狀及斷裂帶斷層泥等做出一定程度的簡化，進而構建斷層物理模型。

在長期的實踐應用過程中，根據(jù)物理模型及具體問題的研究，這些簡化的斷層模型可以歸納為：劈節(jié)點模型、塊體模型、介質弱化模型以及摩擦接觸模型（和平等，2011）。劈節(jié)點模型本質上是在預設斷層位錯的情況下模擬斷層，因而在同震或震后效應研究中具有極大的優(yōu)勢，得到了廣泛應用；塊體模型本質基礎是Shi（2007）提出的DDA（Discontinuous Deformation Analysis）算法，它將地球動力過程看作為一系列塊體相互作用的過程，斷層其實就是塊體與塊體之間的邊界。介質弱化模型和摩擦接觸模型作為2種原理簡單且便于在一般軟件中設置的模型，在區(qū)域地殼動力學形變特征模擬及斷層活動狀態(tài)模擬中得到大量應用。由于斷層的介質力學參數(shù)對區(qū)域應力場以及形變場均有一定程度的影響（李玉江等，2010），因此在模擬過程中，介質弱化模型是將斷層視為有一定厚度的弱介質層，通常通過降低該介質層彈性模量來模擬斷層對區(qū)域地殼動力過程的影響（Lu et al， 2011，2012；祝愛玉等，2015a；胡勐乾等，2014；楊興悅等，2013）；摩擦接觸模型則是直接在創(chuàng)建好的模型基礎上給斷層面賦摩擦系數(shù)，進而實現(xiàn)斷層兩盤的相對運動（Hergert，Heidbach，2010； Zhu，Zhang，2010； He et al，2013；李玉江等，2014；祝愛玉等，2015b）。但也存在一些問題：①對于介質弱化模型來說，其本質上是通過降低斷層介質彈性模量來近似模擬斷層兩盤的非連續(xù)變形。盡管野外觀測到斷層帶內有糜棱巖或斷層泥等弱力學屬性材料填充，但這種弱力學屬性的填充物已經(jīng)在長期的構造作用下與周邊介質處于力平衡狀態(tài)。因此，在構建有限元模型時直接賦予該填充物一個較低的彈性模量，以期獲得斷層兩盤的運動差異和斷層附近應力狀態(tài)的同時，往往忽略了降低彈性模量，也使得該填充層內部、垂直于斷層走向的力平衡破壞，在斷層內部垂直于斷層方向產(chǎn)生新的附加力。這可能影響模擬結果的可靠性。此外，實驗室關于斷層泥的弱力學屬性對斷層行為的影響也表現(xiàn)在剪切摩擦行為上（馬勝利，1986）。在各向同性彈性理論前提下還可能會出現(xiàn)不合理的負泊松比的情況（董培育，石耀霖，2013）。②摩擦接觸模型可以很好地刻畫斷層兩盤相對運動的非連續(xù)變形，相比介質弱化模型更接近于真實斷層狀態(tài)。但摩擦接觸計算是個強非線性問題，涉及到收斂與否以及計算量等問題（李玉江等，2009），特別是在研究區(qū)域內斷層較多的情況下，如果將每條斷層都構建為摩擦基礎模型，收斂問題將變得更加緊迫，為此研究中還可以將重點斷層考慮為摩擦接觸，而將次要斷層處理為弱化帶模型（Li et al，2015），這一簡化處理提高了計算效率，但又增加了弱化介質模型自身的間題。

因此，如何做到既保證計算效率、避免介質弱化模型可能帶來的不真實信息，又能快速合理模擬斷層兩盤相對運動是地殼動力學數(shù)值模擬的一個基礎性工作。本文在介質弱化帶模型基礎上引入正交各向異性理論，用正交各向異性本構方程替代傳統(tǒng)弱化帶中的線彈性本構方程，實現(xiàn)了彈性模量與剪切模量相互獨立；同時考慮到實際建模過程中可能遇到不同斷層走向、傾角各不相同的情況，通過坐標轉換給出了針對任意走向、傾角的斷層本構關系。最后通過數(shù)值試驗對比分析了介質弱化模型、摩擦接觸模型以及正交各向異性斷層模型。

1 正交各向異性本構關系

假定斷層帶為夾在塊體之間的類似于層合板的一類特殊介質：其在垂直于斷層方向力學性質與周邊塊體相近，不易產(chǎn)生大擠壓或拉張變形；而在平行于斷層運動方向上的剪切模量較小，能夠產(chǎn)生較大的剪切變形。則根據(jù)復合材料力學（沈觀林，胡更開，2006）基礎知識，斷層內應力應變關系在其3個主方向上滿足如下關系：式中：ε=[e11 e22 e33 e12 e23 e13]T，σ=[s11；s22 s33s12 S23 S13]T分別為3個主方向的應變和應力分量；x為柔度矩陣；E，，E2，E3，v12，V23，v13，G12，G23，G13分別為其3個主方向的工程彈性常數(shù)。K可表示為：

由于本文重點考慮斷層介質在走向和傾向的剪切行為，其他方向上的介質力學性質與周邊介質一致。為此將式（2）進一步簡化為：式中：E，v分別為周邊介質彈性模量和泊松比；Gd；P為斷層在剪切方向上的剪切模量，獨立于周邊介質參數(shù)E，v。

在實際操作過程中，由于區(qū)域斷層較多且其走向和傾角各不相同，因此需要將其主方向下的應力應變矩陣旋轉變換到地理坐標（x，y，z）下。圖1為斷層面主方向（1，2，3）以及地理坐標系（x，y，z）相對關系，其中：1，2，3分別為斷層走向（面向走向，右手在上盤）、傾向（逆沖為正）以及張拉方向（張為正）；x，y，z分別為地理東、北向以及垂直地表向上的方向。α為斷層走向，β為斷層傾角。令斷層面應變和應力分量在地理坐標系（x，y，z）下分別為：

e=[exx eyy ezz exy eyz exz]T

s=[sxx syy szz sxy syz sxz]（4）

由圖1b可知，將地理坐標系下的應力張量繞z軸逆時針旋轉，然后再繞主軸1逆時針旋轉β便可得到主方向下的應力張量。則：式中：T=A*B*結合式（1）和（4），式（5）可以轉化為：式中：F為坐標轉換矩陣，可表示為：同理可得：將式（6）和式（7）帶入式（1），可得：

F·e=K·F·s（8）即：

e=D·s=F-1·K·F·s（9）式中：D為正交各向異性材料在地理坐標系（x、y、z）下的柔度矩陣，且：式中：柔度矩陣D為21個獨立參數(shù)的對稱矩陣，因此在實際模擬計算過程中，只需修改正交各向異性本夠關系柔度矩陣K在其相應主方向下對應滑動方向的剪切模量（式（3）中的Gslip），結合坐標轉換矩陣F，代入式（10）便可生成斷層有限元模型的柔度矩陣，進而可以構建有限元計算的整體剛度矩陣，實現(xiàn)垂直斷層方向無大變形，平行斷層滑動方向有大變形滑動的斷層運動狀態(tài)。為方便控制，本文令斷層走向和傾向剪切模量與周邊彈性介質剪切模量之比為λ∈[0，1]，即Gslip=λ

2 數(shù)值試驗對比研究

2.1 模型構建

為檢驗正交各向異性本構關系模擬斷層的可行性，本文設計1個虛擬理想彈性地殼模型，模型為邊長L=6000m的正方形，縱向（z）深2000m且上下地殼各1000m厚，模型中部沿x軸走向有一傾角90°的右旋走滑斷裂帶，長2000m，寬50m。模型施加狄里克萊邊界條件：-y邊界約束y方向位移為0，+y邊界只有y方向位移，Uy=-0.1m；-x邊界上施加一反正切函數(shù)的x方向剪切位移，如圖2a所示；模型-z底面約束z方向位移為0。假定模型彈性模量E=84GPa，泊松比v=0.25。然后針對常用的介質弱化模型，摩擦接觸模型以及本文正交各向異性模型分別設計2組工況：第一組工況（Case 1）假定上地殼斷層閉鎖，下地殼可以剪切滑動；第二組工況（Case 2）認為整個斷層處于非閉鎖狀態(tài)，斷層上下均可以相對自由滑動，如圖2b所示。

2.2 結果對比

在一個地震輪回過程中，斷層在地震發(fā)生時由于快速滑動而破裂，使得其力學屬性快速降低；在震后經(jīng)歷長時間的愈合，在下一次地震來臨前，斷層力學屬性又增加到相對較強的狀態(tài)。根據(jù)地震輪回的理論模式（Savage，Burford，1973；Scholz，1998），借鑒Hergert和Heidbach（2010）在馬爾馬拉海的工作思路，本文首先用介質弱化模型、摩擦接觸模型以及本文給出的正交各向異性本構模型分別模擬一組震間地殼運動位移場圖像。假定上地殼己經(jīng)經(jīng)歷了長久的愈合，斷層力學屬性與周邊地殼介質基本一致，而深部斷層強度低，依然可以自由剪切，斷層模型參數(shù)設置如圖2b（Case 1），模擬結果如圖3所示。

基于3種斷層模型地殼運動位移場圖像的空間分布形態(tài)基本一致，但一般文獻多用的介質弱化模型在斷層附近局部，x方向的位移場存在一定程度的畸變，但并未影響Ux沿y的反正切漸變分布特征。另外在空間分布圖像上，基于介質弱化斷層模型的Ux和基于摩擦接觸模型及正交各向異性模型的Ux結果在定量上有所差異。而基于摩擦接觸斷層模型和各向異性斷層模型的對比結果顯示，2者在位移場圖像上高度一致，這與本文構建模型的初衷一致，因此在考慮斷層較多，摩擦接觸建模復雜，計算收斂困難的時候可以考慮用各向異性本構關系來模擬斷層，其模型構建與計算量與常用的介質弱化模型完全一致，計算中也不會涉及到非線性收斂困難等情況，不同的是有限元方程中的柔度矩陣（或剛度矩陣）內的具體數(shù)值有所不同。

震后短期內，斷層沒有完全愈合，其介質力學屬性依然較低。在分析這些介質力學屬性本來偏低的斷層時，本文將整個斷層的力學參數(shù)降低，具體如圖2b（Case 2）所示。結果顯示基于3種斷層模型地殼運動位移場圖像空間分布不同程度出現(xiàn)差異，常用的介質弱化模型在水平和垂直位移場上明顯與其他2種模擬不一致，在斷層附近位移場圖像出現(xiàn)巨大的畸變，這與斷層彈性模量的降低引起垂直于斷層方向發(fā)生大幅度的擠壓變形有關。實際上由于長期的構造作用，斷層兩盤的相對運動一般只發(fā)生在其走向和傾向滑動方向，并不會在斷層上存在這樣垂直于斷層的大幅擠壓變形，因此在應用簡單的介質弱化斷層模型模擬地殼形變時，隨著彈性模量的降低，斷層附近的模擬位移場圖像將出現(xiàn)不同程度畸變，進而可能會影響到斷層附近的模擬應力狀態(tài)（董培育，石耀霖，2013）等。對比基于摩擦接觸模型和本文各向異性模型的模擬結果，盡管2者在斷層附近有一定程度的差異，但整體上2者空間分布形態(tài)完全相似，可以清晰地給出斷層兩盤的相對滑動，所不同的是由于摩擦接觸模型給斷層賦予的是摩擦系數(shù)，而本文只是近似給的剪切模量，2者雖然在阻礙變形的機制上相似，但其背后的物理方程不同，因此造成2者在斷層兩盤相對滑動量的差異。后續(xù)的研究需要基于地殼運動位移場圖像討論2者之間的統(tǒng)計關系，這樣會進一步提高本文各向異性斷層模型在模擬中的應用和理解。

為進一步分析3類斷層模型的異同，從圖3和圖4中切出x=0剖面上的Ux和Uy做對比分析，如圖5所示。結果顯示當在Case 1（淺部斷層愈合，深部斷層可以滑動）情況下，總體上這3類斷層模型給出的結果基本一致，但從更細致的角度看，基于摩擦接觸斷層模型和正交各向異性斷層模型給出的結果吻合一致。但介質弱化模型給出的結果與其他2種斷層模型給出的模擬結果略有差異，并且這種差異在空間上呈現(xiàn)近斷層處相對較小，遠離斷層處增大。這是由于淺部斷層與周邊介質屬性一致，保證了表層形變場的連續(xù)性，但基于介質弱化的深部斷層模型破壞了垂直于斷層方向的平衡狀態(tài)，產(chǎn)生附加的形變場畸變信息。進一步降低表層斷層的介質力學屬性（即Case 2），這種畸變信息對表層形變場的影響更加顯著，特別是在垂直于斷層方向上（圖Sd）。但在平行于斷層方向上，3類斷層模型給出的模擬結果均顯示出不同程度的差異，這是因為本文針對3類斷層模型的力學參數(shù)都是隨機賦值，沒有探討3者力學參數(shù)之間的等效關系。但3者均可以表征本文預設的斷層走滑行為。對于垂直于斷層方向，圖Sd顯示基于正交各向異性理論的斷層模型與摩擦接觸斷層模型的模擬結果一致，表明正交各向異性斷層模型相比介質弱化模型，更接近于實際情況。

為進一步分析不同斷層模型對斷層附近模擬應力狀態(tài)的影響，考慮到圖2設計的模型表現(xiàn)為剪切變形行為，以水平剪應力（τyx）的空間分布為討論對象。在圖2的基礎上繼續(xù)分析了不同1況下，不同斷層模型對應的τyx模擬結果，如圖6所示。結果顯示：較介質弱化模型給出的τyx空間分布（圖6a，d），正交各向異性斷層模型的模擬結果（圖6b，e）與摩擦接觸模型的結果（圖6c，f）更為相近。在上地殼斷層閉鎖的情況下（Case 1），水平剪應力τyx云圖在空間上相對y=0軸對稱（圖6b、c），而介質弱化模型給出的結果不同于其他兩個模型的模擬結果，水平剪應力呈現(xiàn)出不對稱形態(tài)（圖6a）；在整個斷層均貫通非閉鎖的情況下（Case 2），盡管正交各向異性斷層模型的模擬結果（圖6e）與摩擦接觸模型（圖6f）的模擬結果有一定的差異，但2者基本保持相似的空間分布特征，而基于介質弱化模型給出的模擬結果（圖6d）則在空間分布上表現(xiàn)出不規(guī)則的分布形態(tài)。這些都是由于在數(shù)值模擬中，介質弱化模型由于其本構關系的原因，不可避免地在垂直于斷層方向上引入了較大的擠壓變形（圖4d），這種不合理的擠壓變形不但影響著位移場的空間分布特征，同時也影響著應力場的空間分布。

3 討論

近年來，隨著GPS地殼速度場觀測數(shù)據(jù)的豐富，基于GPS約束的地殼動力學數(shù)值模擬研究日漸深入（Hergert，Heidbach，2010；He et al，2013；Lei et al，2010； Cho，Kuwahara，2013；龍小剛，朱守彪，2015；劉峽等，2010）。一般情況下，研究人員通常都會通過一些途徑選定有限元模型的邊界條件，模擬計算給出模型內部的速度場模擬結果，然后通過對比模擬結果與觀測結果之間的殘差來確定模型是否可靠，參數(shù)選擇是否合理，在大量比較之后給出合理模型的基礎上分析區(qū)域的動力學機制及應力狀態(tài)等，這就要求我們首先要保證所采取的模型是合理的。介質弱化模型作為最為簡單的斷層模型，其建模方便、不涉及非線性問題且不會存在難于收斂的情況，在一些研究中得到大量的應用。因此在應用介質弱化模型的情況下，在對比模擬地殼運動速度場與GPS觀測速度場殘差的同時，如何評價類似于圖4a、d以及圖5d的形變場畸變對模擬結果的影響也是不容忽視的問題。在不可避免使用介質弱化斷層模型的時候，如能進一步修改其有限元剛度矩陣，使其剪切模量與彈性模量相互獨立，對于提高模擬結果的可靠性具有一定的意義。

4 結論

本文在介質弱化帶模型基礎上，引入正交各向異性理論，同時考慮到實際建模過程中可能遇到不同斷層的走向、傾角各不相同的情況，本文通過坐標轉換給出了針對任意走向、傾角的正交各向異性斷層本構關系。最后通過數(shù)值試驗對比分析了介質弱化模型、摩擦接觸模型以及正交各向異性斷層模型在震間地殼速度長圖像模擬中的差異，結果顯示：（1）無論斷層力學性質強弱，傳統(tǒng)的介質弱化模型都會對地表地殼運動位移場產(chǎn)生一定程度的影響，斷層力學性質（對應于圖2不同斷層模型中斷層的力學參數(shù)，如弱化介質模型的彈性模量、各向異性模型的剪切模量以及摩擦接觸模型的摩擦系數(shù)）越弱，影響越顯著，這主要由于是降低斷層彈性模量，使得斷層內垂直于斷層方向的力平衡被打破，進而影響了斷層附近位移場圖像的空間分布，這在模擬中同時也會對斷層附近的應力場和應變場產(chǎn)生一定程度的影響：（2）從位移場圖像以及應力場圖像的空間分布來看，正交各向異性本構關系與摩擦接觸模型給出的模擬結果基本一致，因此在區(qū)域構造復雜的情況下，可以考慮使用各向異性本構關系來表征斷層的變形機制。

非常感謝陳連旺研究員的悉心指導，同時也感謝中國地震局訪問學者計劃的支持和地殼應力研究所提供良好的辦公環(huán)境和軟件支持！部分圖件采用GMT（Wessel et al，2013）繪制，特此表示感謝！

參考文獻：

陳連旺，陸遠忠，張杰，等.1999.華北地區(qū)三維構造應力場[J].地震學報，21（2）：140-149.

陳連旺，張培震，陸遠忠，等.2008.川滇地區(qū)強震序列庫侖破裂應力加卸載效應的數(shù)值模擬[J].地球物理學報，51（5）：1411-1421.

董培育，石耀霖.2013.關于“用單元降剛法探索中國大陸強震遠距離跳遷及主體活動區(qū)域轉移”的討論——橫向各向同性“殺傷單元”才是更好的途徑[J].地球物理學報，56（6）：2133-2139.和平，李志雄，陸遠忠，等.2011.斷層面的有限單元模擬方法綜述[J].西北地震學報，33（2）：186-194.

胡勐乾，鄧志輝，陸遠忠，等.2014.三維數(shù)值模擬在華北地區(qū)現(xiàn)今構造變形分析中的應用[J].地震地質，36（1）：148-165.

李玉江，陳連旺，劉少峰，等.2014.蘆山地震的發(fā)生對周圍斷層影響的數(shù)值模擬[J].地球學報，35（5）：627-634.

李玉江，陳連旺，葉際陽，等.2010.巖石物性變化對區(qū)域應力場的影響[J].地球物理學進展，25（6）：1941-1946：

李玉江，陳連旺，葉際陽.2009.數(shù)值模擬方法在應力場演化及地震科學中的研究進展[J]地球物理學進展，24（2）：418-431.

劉峽，馬瑾，傅容珊，等.2010.華北地區(qū)現(xiàn)今地殼運動動力學初步研究[J].地球物理學報，53（6）：1418-1427.

龍小剛，朱守彪.2015.臺灣碰撞帶現(xiàn)今地殼變形場特征及其動力學成因的有限單元法模擬研究[J].地球物理學報，58（7）：2350-2365：

馬勝利，本利彥.1995.蒙脫石的脫水作用對斷層摩擦本構行為的影響[J].地震地質，17（4）：289-304.

馬勝利.1986.模擬斷層帶摩擦滑動性狀與變形特征[J].中國地震，（2）：74-80.

邵志剛，傅容珊，薛霆虓，等，2008.昆侖山MS8.1級地震震后變形場數(shù)值模擬與成因機理探討[J].地球物理學報，51（3）：805-816.

沈觀林，胡更開.2006.復合材料力學[M].北京：清華大學出版社.

楊興悅，陳連旺，楊立明，等.2013.巴顏喀拉塊體強震動力學過程數(shù)值模擬[J].地震學報，35（3）：304-314.

祝愛玉，張東寧，蔣長勝，等.2015a.川滇地區(qū)地殼應變能密度變化率與強震復發(fā)間隔的數(shù)值模擬[J].地震地質，37（3）：906-927.

祝愛玉，張東寧，蔣長勝.2015b.安寧河—則木河—小江斷裂帶應力狀態(tài)分段特征的數(shù)值模擬研究[J].中國科學：地球科學，45（12）：1839-1852.

Cho I，Kuwahara Y.2013.Numerical simulation of crustal deformation u-sing a three-dimensional viscoelastic crustal structure model for theJapanese islands under east-west compression[J].Earth，Planetsand Space，65（9）：1041-1046.

He J，Lu S，Wang W.2013.Three-dimensional mechanical modeling ofthe GPS velocity field around the northeastern Tibetan plateau andsurrounding regions[J].Tectonophysics，584；257-266.

Hergert T，Heidbach 0.2010.Slip-rate variability and distributed de-formation in the Marmara Sea fault system[J].Nature Geoscience，3（2）：132-135.

Lei X，Chen Y，Zhao J，et al.2010.Modelling of current crustal tectonicdeformation in the Chinese Tianshan orogenic belt constrained byGPS observations[J].Journal of Geophysics and Engineering，7（4）：431.

Li S，Moreno M，Rosenau M，et al.2014.Splay fault triggering by greatsubduction earthquakes inferred from finite element models[J].Geo-physical Research Letters，41（2）：385-391.

Li Y，Chen L，Liu S，et al.2015.Coseismic Coulomb stress changescaused by the MW6.9 Yutian earthquake in 2014 and its correlationto the 2008 MW7.2 Yutian earthquake[J].Journal of Asian EarthSciences，105：468-475.

Lu Y，Yang S，Chen L，et al.2011.Mechanism of the spatial distributionand migration of the strong earthquakes in China inferred from nu-merical simulation[J].Journal of Asian Earth Sciences，40（4）：990-1001.

Lu Y，Yang S，Chen L，et al.2012.Migration trend of strong earthquakesin North China from numerical simulations[J].Journal of AsianEarth Sciences，50：116-127.

Moreno M，Rosenau M，0ncken O.2010.2010 Maule earthquake slip cor-relates with pre-seismic locking of Andean subduction zone[J].Nature，467（7312）：198-202.

Savage J C，Burford R O.1973.Geodetic determination of relative platemotion in central California[J].Journal of Geophysical Research，78（5）：832-845.

Scholz C H.1998.Earthquakes and friction laws[J].Nature，391（6662）：37-42.

Shi G.2007.Applications of discontinnons deformation analysis（DDA）to rock engineering[M].Berlin，Heidelberg，Computational Mechan-ics springer；136-147.

Sibson R H.1977.Fault rocks and fault mechanisms[J].Journal of theGeological Society，133（3）：191-213.

Wessel P，Smith W H，Scharroo R，et al.2013.Generic mapping tools：Im-proved version released[J].Eos，Transactions American GeophysicalUnion，94（45）：409-410.

Zhu S，Zhang P.2010.Numeric modeling of the strain accumulation andrelease of the 2008 Wenchuan，Sichuan，China，Earthquake[J].Bul-letin of the Seismological Society of America，100（5B）：2825-2839.