張龍，江在森，武艷強(qiáng)，2＊，鄒鎮(zhèn)宇，3，劉曉霞，3，魏文薪

1中國(guó)地震局地震預(yù)測(cè)研究所地震預(yù)測(cè)重點(diǎn)實(shí)驗(yàn)室，北京 100036

2中國(guó)地震局第一監(jiān)測(cè)中心，天津 300180

3中國(guó)地震局地質(zhì)研究所地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029

1 引言

非連續(xù)變形分析（DDA）方法是巖土工程領(lǐng)域模擬節(jié)理巖體運(yùn)動(dòng)與變形的一種數(shù)值解析方法（Shi，1992）.該方法將不連續(xù)面分割的巖體視為塊體集合，以最小勢(shì)能法建立總體平衡方程，利用單純形積分與開閉迭代方法保證結(jié)果的精確性與計(jì)算的高效性，最終求得各塊體的位移與變形（Shi，2001；Wu et al.，2012）.DDA方法因收斂速度高以及模擬結(jié)果精確等優(yōu)勢(shì)，廣泛應(yīng)用于邊坡破壞（Sitar et al.，2005；Wu，2010；Wu and Chen，2011）、洞室開挖（鄔愛清等，2006；Law and Lam，2003）以及爆破（劉剛等，2003）等動(dòng)力學(xué)模擬研究中.此外，DDA方法在模擬區(qū)域地殼變形與應(yīng)變積累領(lǐng)域也得到了初步應(yīng)用（陳兵等，2000；陳祖安等，2008；秦小軍等，2001；王輝等，2003；Liu et al.，1996）.

DDA方法中塊體單元之間的接觸界面是位移間斷面，界面滑動(dòng)與否由Mohr－Coulomb準(zhǔn)則控制，當(dāng)剪應(yīng)力τ滿足公式（1）時(shí)，接觸狀態(tài)由鎖定變?yōu)槟Σ粱瑒?dòng).

式中σn為接觸邊界上的正應(yīng)力，φ和C分別是摩擦角和黏著力（Shi，1992，2001），摩擦系數(shù)μ（＝tanφ）在模擬過程中保持不變.由于巖土工程主要針對(duì)水壩、邊坡等百米尺度、巖性相對(duì)均一的研究目標(biāo)，采用常摩擦系數(shù)進(jìn)行數(shù)值模擬已足夠精確.然而，針對(duì)更大尺度的地學(xué)問題，實(shí)驗(yàn)室測(cè)量與野外調(diào)查結(jié)果均表明，地震成核、震前滑動(dòng)、同震破裂、余震和震后余滑以及震間期斷層的愈合過程中，斷層面的摩擦系數(shù)并非常數(shù)（Marone，1998）.因此，為合理解釋地震周期各階段發(fā)震斷層行為及剪應(yīng)力變化過程，需要在數(shù)值模擬時(shí)采用更符合實(shí)際的摩擦本構(gòu)定律（張龍等，2013）.

基于實(shí)驗(yàn)室摩擦實(shí)驗(yàn)結(jié)果得到的速度－狀態(tài)摩擦本構(gòu)定律能夠模擬斷層面上的成核過程以及發(fā)生在俯沖帶的周期性慢滑移，結(jié)合數(shù)值模擬方法還可直觀給出地震周期的輪廓 （Lapusta and Rice，2003；Liu and Rice，2007）.該定律中，斷層面的摩擦強(qiáng)度由滑動(dòng)速度V和狀態(tài)變量θ共同描述，其基本形式如下（Dieterich，1979，1981；Ruina，1983）：

為了實(shí)現(xiàn)復(fù)雜介質(zhì)的地震孕育、發(fā)生及震后調(diào)整全過程的數(shù)值模擬研究，首先需要解決速度－狀態(tài)摩擦本構(gòu)定律與動(dòng)力學(xué)數(shù)值模擬方法的融合問題.本文選擇三維DDA方法作為研究對(duì)象，首先推導(dǎo)了基于速度－狀態(tài)摩擦本構(gòu)定律的三維DDA模擬實(shí)用公式，隨后設(shè)計(jì)了滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)和速度步進(jìn)實(shí)驗(yàn)兩算例以驗(yàn)證模擬效果，并分析了數(shù)值模擬結(jié)果與實(shí)驗(yàn)室摩擦實(shí)驗(yàn)結(jié)果存在細(xì)微差別可能的原因.

2 計(jì)算方法

由公式（2）和公式（3）知，塊體接觸面上每一時(shí)間步的摩擦系數(shù)μ由當(dāng)前時(shí)步的速度V和狀態(tài)變量θ共同決定，其他參數(shù)（μ0，V0，Dc，a和b）均為常數(shù).下面詳細(xì)推導(dǎo)了改進(jìn)的三維DDA算法中計(jì)算摩擦系數(shù)的實(shí)用公式.

2.1 速度

應(yīng)用三維DDA方法求解塊體速度時(shí)，選擇動(dòng)力學(xué)計(jì)算模式，即每一時(shí)步結(jié)束時(shí)的速度為下一時(shí)步開始時(shí)的速度.D＝ （uvwrxryrzεxεyεzγxyγyz γzx）T表示塊體變形量，包含塊體位移、旋轉(zhuǎn)和應(yīng)變.令Δt為時(shí)步長(zhǎng)，Di和Di＋1為t和t＋Δt時(shí)刻的塊體變形量，采用向前差分格式（Shi，1992，2001），則有

在小步長(zhǎng)情況下可假定每一時(shí)步內(nèi)的加速度為常數(shù)，即

令當(dāng)前時(shí)步初始變形量為0，即Di＝0，公式（4）略去高階無(wú)窮小并代入公式（5）得

勻加速直線運(yùn)動(dòng)速度計(jì)算表達(dá)式為

部分公民環(huán)保意識(shí)薄弱，碳信息披露可以促進(jìn)企業(yè)可持續(xù)發(fā)展，但企業(yè)對(duì)此認(rèn)識(shí)不夠深刻，這些都使碳披露難以發(fā)揮其對(duì)經(jīng)營(yíng)管理的指導(dǎo)作用。這就需要社會(huì)廣泛宣傳低碳經(jīng)濟(jì)理念，提高公民的環(huán)保意識(shí)，進(jìn)而逐步上升到企業(yè)管理層面，使管理者提高環(huán)保責(zé)任意識(shí)，加大對(duì)碳信息的重視。通過學(xué)習(xí)，掌握必要的碳披露知識(shí)，推動(dòng)經(jīng)濟(jì)社會(huì)可持續(xù)發(fā)展。

將公式（6）代入公式（7），整理得速度迭代表達(dá)式

令Vi和Vi＋1分別表示，則（8）式可表示為

2.2 狀態(tài)變量

根據(jù)公式（2）或（3）的演化方程計(jì)算得到狀態(tài)變量θ.以公式（2）為例，將演化方程用4階Runge－Kutta方法離散，即

式中θi和θi＋1分別為t和t＋Δt時(shí)刻的狀態(tài)變量，Kj（j＝1，2，3，4）為中間變量.初始時(shí)刻的θ值由公式（2）（或（3））的主方程求得.

根據(jù)公式（9）和（10）分別計(jì)算得到下一時(shí)步的速度V和狀態(tài)變量θ，代入公式（2）（或（3））可計(jì)算得到下一時(shí)步的摩擦系數(shù)μ，完成該時(shí)間步計(jì)算.將上述算法嵌入三維DDA程序中，建立了應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法.

3 數(shù)值模擬

為驗(yàn)證結(jié)合速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法在計(jì)算地震周期中斷層行為問題時(shí)的精確性，本文設(shè)計(jì)了滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)和速度步進(jìn)實(shí)驗(yàn)兩個(gè)算例.

3.1 滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)

研究表明，相比于其他形式的本構(gòu)定律，Dieterich law能夠準(zhǔn)確描述應(yīng)力降與保持時(shí)間的關(guān)系，即靜摩擦的時(shí)間依賴性（Dieterich，1978）.最初，實(shí)驗(yàn)室做滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)使用“三明治”直剪實(shí)驗(yàn)機(jī)（圖1a），Scholz（2002）將其簡(jiǎn)化為圖1b形式.

圖1 （a）“三明治”直接剪切實(shí)驗(yàn)機(jī)；（b）Scholz簡(jiǎn)化的直接剪切實(shí)驗(yàn)機(jī)Fig.1 （a）Schematic diagram of“sandwich”direct shear apparatus；（b）Direct shear apparatus modifid by Scholz

利用改進(jìn)的三維DDA方法數(shù)值模擬該實(shí)驗(yàn)，需要設(shè)計(jì)形狀合適的模型并合理選擇材料的物理參數(shù).考慮到摩擦強(qiáng)度變化與接觸面積大小無(wú)關(guān)，為簡(jiǎn)化模擬，對(duì)模型幾何形狀略作修改，如圖2所示.模型的介質(zhì)參數(shù)、數(shù)值控制參數(shù)及摩擦本構(gòu)參數(shù)見表1.

表1 滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)數(shù)值參數(shù)表Table 1 Numerical parameters of slide－h(huán)old－slide tests for 3DDDA forward modeling simulation

為與實(shí)驗(yàn)室摩擦實(shí)驗(yàn)數(shù)據(jù)（Beeler et al.，1994）進(jìn)行對(duì)比，本文分別做了保持時(shí)間為1s、10s、100s、1000s和10000s的滑動(dòng)－保持－滑動(dòng)數(shù)值實(shí)驗(yàn).限于篇幅，僅以圖3為例給出保持時(shí)間為10s的滑動(dòng)－保持－滑動(dòng)數(shù)值實(shí)驗(yàn)?zāi)M結(jié)果.

圖2 滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)的模型幾何Fig.2 Numerical model for slide－h(huán)old－slide tests

圖3 滑動(dòng)－保持－滑動(dòng)數(shù)值實(shí)驗(yàn)?zāi)M結(jié)果（保持時(shí)間為10s）Fig.3 Numerical results of slide－h(huán)old－slide tests（hold period is 10s）

如圖3所示，前5s滑塊在荷載力和摩擦力的共同作用下做勻速直線運(yùn)動(dòng)，即“滑動(dòng)”過程；5～15s為“保持”階段，即令荷載速度為零，滑塊在摩擦力的作用下做減速運(yùn)動(dòng)，此時(shí)摩擦系數(shù)也發(fā)生相應(yīng)變化；15～20s仍為“滑動(dòng)”階段，15s時(shí)重載初始時(shí)刻荷載速度，可明顯觀察到靜摩擦系數(shù)瞬間達(dá)到峰值，隨著時(shí)間的增加不斷減小，最終恢復(fù)到初始水平.不同保持時(shí)長(zhǎng)的滑動(dòng)－保持－滑動(dòng)數(shù)值實(shí)驗(yàn)?zāi)Σ料禂?shù)演化過程與該結(jié)果類似.

圖4給出了靜摩擦變化量Δμs（峰值摩擦系數(shù)與穩(wěn)態(tài)值的差值）與保持時(shí)間的關(guān)系.實(shí)驗(yàn)室摩擦實(shí)驗(yàn)結(jié)果與地震學(xué)對(duì)斷層愈合速率的估計(jì)均表明，準(zhǔn)靜態(tài)接觸時(shí)摩擦愈合量Δμs與保持時(shí)間的對(duì)數(shù)呈線性關(guān)系（Beeler et al.，1994），即所謂靜摩擦?xí)r間依賴性.圖4表明，模擬結(jié)果與實(shí)驗(yàn)室測(cè)量結(jié)果（Beeler et al.，1994）基本吻合.當(dāng)保持時(shí)間為101～104s時(shí)，模擬數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)對(duì)數(shù)擬合后呈線性變化且斜率相同，擬合優(yōu)度均為0.999，將兩部分?jǐn)?shù)據(jù)混合再進(jìn)行對(duì)數(shù)擬合的擬合優(yōu)度依然高達(dá)0.997.由此得出結(jié)論，應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法能夠比較準(zhǔn)確地描述靜摩擦的時(shí)間依賴性.另外，摩擦愈合速率β（＝Δμs每十年的變化量）.在速度接近零時(shí)可由摩擦本構(gòu)參數(shù)b衡量（Beeler et al.，1994）.換言之，用該方法可求出Dieterich law中本構(gòu)參數(shù)b.圖4中保持時(shí)間為101～104s的數(shù)值模擬結(jié)果斜率為0.00856，與摩擦本構(gòu)參數(shù)b（＝0.009）十分接近，相對(duì)誤差約為5%，表明應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法進(jìn)行數(shù)值模擬可以基本滿足精度要求.

圖4 滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)的數(shù)值結(jié)果與實(shí)驗(yàn)室結(jié)果圖中三角形表示三維DDA的數(shù)值實(shí)驗(yàn)結(jié)果，圓點(diǎn)和菱形表示Beeler等（1994）實(shí)驗(yàn)室摩擦實(shí)驗(yàn)結(jié)果.實(shí)線代表用Dieterich law對(duì)滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)數(shù)值模擬的結(jié)果.Fig.4 Comparison between numerical results and laboratory measurements for slide－h(huán)old－slide tests In this figure，numerical results performed by 3DDDA are showed as triangles.Circles and diamonds represent the laboratory data measured by Beeler et al.（1994）.The curve represents numerical results performed by traditional integral method with Dieterich law.

3.2 速度步進(jìn)實(shí)驗(yàn)

速度步進(jìn)實(shí)驗(yàn)可用來(lái)研究擴(kuò)容作用及斷層泥對(duì)動(dòng)摩擦速度依賴性的影響 （Dieterich and Kilgore，1994）.研究表明，Ruina law能夠準(zhǔn)確描述摩擦系數(shù)對(duì)速度變化的響應(yīng)，即動(dòng)摩擦的速度依賴性（何昌榮，1999；Dieterich and Kilgore，1994），因而在該算例中選擇公式 （3）計(jì)算摩擦系數(shù).本實(shí)驗(yàn)選擇與滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)相同的幾何模型（圖2），介質(zhì)參數(shù)、數(shù)值控制參數(shù)及摩擦本構(gòu)參數(shù)見表2.

模擬結(jié)果如圖5b所示，荷載點(diǎn)速度為10μm／s，滑塊穩(wěn)定滑動(dòng)；運(yùn)動(dòng)25μm后將荷載點(diǎn)速度突然降至1μm／s，運(yùn)動(dòng)一段距離后滑塊恢復(fù)穩(wěn)定滑動(dòng)狀態(tài)，重復(fù)上述過程一次.

圖5表明，應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法能夠有效模擬動(dòng)摩擦的速度依賴性.選擇與摩擦實(shí)驗(yàn)相同的加載方式和加載條件，三維DDA方法能夠比較準(zhǔn)確地描述摩擦系數(shù)對(duì)速度由高值到低值再到高值的響應(yīng)，曲線形態(tài)與實(shí)驗(yàn)結(jié)果基本吻合.值得注意的是，數(shù)值模擬曲線中摩擦系數(shù)在速度突變時(shí)的響應(yīng)較實(shí)驗(yàn)曲線銳利，通過設(shè)置不同的數(shù)值彈簧剛度及介質(zhì)材料參數(shù)，并未使得模擬曲線峰尖變緩.分析表明，物理實(shí)驗(yàn)的加載（速度變化）做不到數(shù)值模擬那樣的瞬時(shí)響應(yīng)，因而存在這種差別是合理的.

表2 速度步進(jìn)實(shí)驗(yàn)數(shù)值參數(shù)表Table 2 Numerical parameters of velocity stepping tests for 3DDDA forward modeling simulation

圖5 速度步進(jìn)摩擦實(shí)驗(yàn)結(jié)果與數(shù)值模擬結(jié)果Fig.5 Comparison between numerical results and laboratory measurements for velocity stepping tests

4 結(jié)論

本文通過公式推導(dǎo)、模型驗(yàn)證等過程，拓展了基于速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法，取得如下認(rèn)識(shí)和結(jié)論：

（1）從速度－狀態(tài)摩擦本構(gòu)定律的基本形式出發(fā)，通過推導(dǎo)給出了融合該定律的三維DDA方法計(jì)算摩擦系數(shù)的實(shí)用公式，拓展了三維DDA方法.

（2）滑動(dòng)－保持－滑動(dòng)實(shí)驗(yàn)和速度步進(jìn)實(shí)驗(yàn)?zāi)M結(jié)果表明，應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法可以比較準(zhǔn)確地模擬實(shí)驗(yàn)中觀察到的靜摩擦?xí)r間依賴性和動(dòng)摩擦速度依賴性，進(jìn)而表明該方法模擬接觸面剪應(yīng)力動(dòng)態(tài)變化是有效的.

（3）同原有方法相比，改進(jìn)后的三維DDA方法解決了應(yīng)用于大尺度構(gòu)造應(yīng)力演化和區(qū)域應(yīng)變積累的基本問題，即Coulomb定律中的常摩擦系數(shù)不能反映地震周期中斷層面剪應(yīng)力變化.然而，改進(jìn)的三維DDA方法在模擬中依然存在由于不連續(xù)面剛度的不確定性，發(fā)生嵌入深淺不一的問題，往往會(huì)影響求解接觸力或應(yīng)力的精度.相信隨著三維DDA凹?jí)K體接觸算法不斷完善，本文這種改進(jìn)將有助于今后與觀測(cè)資料（GPS等）相結(jié)合，定量描述區(qū)域地殼變形特征.

致謝感謝石根華先生、九州大學(xué)陳光齊教授提供最新版本的三維DDA代碼，感謝中國(guó)地震局地質(zhì)研究所何昌榮研究員關(guān)于數(shù)值模擬技術(shù)的有益指導(dǎo)，感謝兩位匿名審稿專家為完善本文提出的寶貴意見.

Beeler N M，Tullis T E，Weeks J D.1994.The roles of time and displacement in the evolution effect in rock friction.Geophys.Res.Lett.，21（18）：1987－1990.

Chen B，Jiang Z S，Wang S X，et al.2000.Preliminary study on block movement and its stress field by discontinuous deformation analysis.CrustalDeformationandEarthquake（in Chinese），20（1）：38－42.

Chen Z A，Lin B H，Bai W M.2008.3－D numerical simulation on influence of 1997Mani earthquake occurrence to stability of tectonic blocks system in Qingzang and Chuandian zone using DDA＋FEM method.ChineseJ.Geophys.（in Chinese），51（5）：1422－1430，doi：10.3321／j.issn：0001－5733.2008.05.015.Dieterich J H.1978.Time－dependence friction and the mechanics of stick－slip.PureAppl.Geophys.，116（4－5）：790－806.

Dieterich J H.1979.Modeling of rock friction：1.Experimental results and constitutive equations.J.Geophys.Res.，84（B5）：2161－2168.

Dieterich J H.1981.Constitutive properties of faults with simulated gouge.／／Carter N L，F(xiàn)riedman M，Logan J M eds.Mechanical Behavior of Crustal Rocks：The Handin Volume.Washington，DC：Am.Geophys.Union：24：103－120.

Dieterich J H，Kilgore B D.1994.Direct observation of frictional contacts：new insights for state－dependent properites.Pure Appl.Geophys.，143（1－3）：283－302.

He C R.1999.Comparing two types of rate and state dependent friction laws.SeismologyandGeology（in Chinese），21（2）：137－146.

Lapusta N，Rice J R.2003.Nucleation and early seismic propagation of small and large events in a crustal earthquake model.J.Geophys.Res.，108（B4），doi：10.1029／2001JB00793.Law H K，Lam I P.2003.Evaluation of seismic performance for tunnel retrofit project.J.Geotech.Geoenviron.Eng.，129（7）：575－589.

Liu G，Shu D Q，Jiang Q H.2003.Application of discontinuous deformation analysis method in analyzing shaft well stabilization.Blasting（in Chinese），20（4）：38－44.

Liu L B，Linde A T，Sacks I S，et al.1996.A seismic fault slip and block deformation in north China.PureAppl.Geophys.，146（3－4）：717－740.

Liu Y J，Rice J R.2007.Slow slip predictions based on granite and gabbro friction data compared to GPS measurements in northern Cascadia.J.Geophys.Res.，114（B9）：B09407，doi：10.1029／2008JB006142.

Marone C.1998. Laboratory－derived friction laws and their application to seismic faulting.Annu.Rev.EarthPlanet.Sci.，26：643－696.

Qin X J，Zhou S Y，Zhao L Q，et al.2001.Research on horizontal crustal movement from 1994to 1998in Jiashi and the northeast to Pamir ＆its relation to strong earthquake swarm on the basis of GPS data.CrustalDeformationandEarthquake（in Chinese），21（4）：1－7.

Ruina A L.1983.Slip instability and state variable friction law.J.Geophys.Res.，88（B12）：10359－10370.

Scholz C H.2002.The Mechanics of Earthquakes and Faulting.2nd ed.New York：Cambridge University Press.

Shi G H.1992.Discontinuous deformation analysis：A new numerical model for the statistics and dynamics of deformable block structures.Engng.Comp.，9（2）：157－168.

Shi G H.2001. Three－dimensional discontinuous deformation analysis.／／Elsworth D ed.Proceedings of the 38th U.S.Rock Mechanics Symposium.Taylor ＆Francis，Washington，D.C.，USA，1421－1428.

Sitar N，MacLaughlin M M，Doolin D M.2005.Influence of kinematics on landslide mobility and failure mode.J.Geotech.Geoenviron.Eng.，131（6）：716－728.

Wang H，Zhang G M，Wu Y，et al.2003.The deformation of active crustal－blocks on the Chinese Mainland and its relation with seismic activity.EarthquakeResearchinChina（in Chinese），19（3）：243－254.

Wu A Q，Ding X L，Chen S H，et al.2006.Researches on deformation and failure characteristics of an underground powerhowse with complicated gelolgical conditions by DDA method.ChineseJ.RockMech.Engng.（in Chinese），25（1）：1－8.

Wu J H.2010.Seismic landslide simulations in discontinuous deformation analysis.Comp.Geotech.，37（5）：594－601.

Wu J H，Chen C H.2011.Application of DDA to simulate characteristics of the Tsaoling landslide.Comp.Geotech.，38（5）：741－750.

Wu Y Q，Chen G Q，Jiang Z S，et al.2012.The algorithm of simplex integration in three－dimension and its characteristic analysis.Int.J.Adv.Comp.Tech.，4（10）：246－256.

Zhang L，Jiang Z S，Wu Y Q.2013.Review of rate－and statedependent friction laws and their applications to seismic faulting.ProgressinGeophysics（in Chinese），28（5）：2352－2362，doi：10.6038／pg20130517.

附中文參考文獻(xiàn)

陳兵，江在森，王雙緒等.2000.利用非連續(xù)變形數(shù)值方法研究塊體運(yùn)動(dòng)及其應(yīng)力場(chǎng)初探.地殼形變與地震，20（1）：38－42.

陳祖安，林邦慧，白武明.2008.1997年瑪尼地震對(duì)青藏川滇地區(qū)構(gòu)造塊體系統(tǒng)穩(wěn)定性影響的三維DDA＋FEM方法數(shù)值模擬.地球物理學(xué)報(bào)，51（5）：1422－1430，doi：10.3321／j.issn：0001－5733.2008.05.015.

何昌榮.1999.兩種摩擦本構(gòu)關(guān)系的對(duì)比研究.地震地質(zhì)，21（2）：137－146.

劉剛，舒大強(qiáng)，姜清輝.2003.應(yīng)用DDA方法分析豎井爆振穩(wěn)定問題.爆破，20（4）：38－44.

秦小軍，周碩愚，趙齊樂等.2001.依據(jù)GPS數(shù)據(jù)研究伽師及帕米爾東北側(cè)1994—1998年地殼水平運(yùn)動(dòng)及其與強(qiáng)震群的關(guān)系.地殼形變與地震，21（4）：1－7.

王輝，張國(guó)民，吳云等.2003.中國(guó)大陸活動(dòng)地塊變形與地震活動(dòng)的關(guān)系.中國(guó)地震，19（3）：243－254.

鄔愛清，丁秀麗，陳勝宏等.2006.DDA方法在復(fù)雜地質(zhì)條件下地下廠房圍巖變形與破壞特征分析中的應(yīng)用研究.巖石力學(xué)與工程學(xué)報(bào)，25（1）：1－8.

張龍，江在森，武艷強(qiáng).2013.速度－狀態(tài)摩擦本構(gòu)定律及其在地震斷層中的應(yīng)用研究進(jìn)展.地球物理學(xué)進(jìn)展，28（5）：2352－2362，doi：10.6038／pg20130517.