和 平，李志雄，陸遠忠，邵志剛

（1.中國地震局地震預測研究所，北京 100036；2.中國地震局地殼應力研究所，北京 100085）

0 引言

地震的孕育、發(fā)生與斷層有著密切的關系［1－5］。為了更好地認識地震發(fā)生過程，前人對斷層模型和機制曾進行過大量的研究。Sibson［6］在1977年提出了將斷層劃分為摩擦和準塑性兩個域進行處理，其中摩擦域溫度低，碎裂作用占主導；準塑性域溫度高，主要由造巖礦物發(fā)生塑性變形所致。Strehlau［7］認為由三種構造流變層組成斷層模型，即斷層泥組成的上部摩擦層、復雜半脆性巖石組成的中部轉換層、糜棱質剪切帶組成的下部延性層。Dieterich［8］和 Gomberg等［9］提出斷層滑動具有摩擦狀態(tài)－速率依從機制，即斷層的摩擦系數(shù)與滑動狀態(tài)和歷史有關。劉桂萍等［10］采用Gomberg建立的具有摩擦狀態(tài)—速率依從機制、地震成核斷層函數(shù)組成的數(shù)學物理模型，模擬大地震后由區(qū)域靜應力場變化引起的區(qū)域觸發(fā)地震現(xiàn)象，并探討相關力學機制。

然而，在地球動力學、地震學、地殼形變學科的研究中，由于實際地質構造形態(tài)復雜性、地球介質分布非均勻性、以及介質物理性質多樣性，基于上述模型求解偏微分方程解析解時遇到了很大的困難。20世紀60年代后期，有限單元數(shù)值模擬方法以其獨特的優(yōu)勢——復雜幾何構型適應性和各種物理問題適用性，被引入到地學研究中，并因它的優(yōu)越性而得到迅速推廣［11－15］。然而，在利用有限單元法在對地學問題進行研究時所面臨到的第一個重要難題就是利用基于連續(xù)介質理論的有限單元法如何引入對具有非連續(xù)性間斷面的處理，即對斷層進行數(shù)值模擬。這一問題引起了很多學者關注，眾多學者也曾對此進行過許多有益嘗試，并提出過多種斷層模擬方法和模型。

本文總結歸納了近年來國內外學者用有限單元法進行地學，特別是地震學、地殼形變學科研究時，處理斷層所采用的主流方法?？傮w而言，斷層處理方式可歸納為兩類，一類是將斷層處理為特殊連續(xù)介質體，即弱化帶模型；另一類則將斷層作為非連續(xù)特性的位錯面來模擬，主要有劈節(jié)點模型、接觸模型以及塊體模型。本文就上述四大斷層模擬方法及其基礎上展開的地震學、形變學研究展開比較詳盡的概述，以期為進行地學數(shù)值模擬時處理斷層提供一些借鑒和參考。

1 介質弱化帶模型

弱化帶模型假定斷層具有一定的厚度并被某種物質充填，斷層和周圍地質體的差別僅在于內在物理屬性和介質參數(shù)的不同，介質彼此仍是連續(xù)的，就像沒有“破”一樣。弱化帶模型從某種程度上說是合理的，因為根據(jù)野外地質調查，斷層帶內的介質常常發(fā)育成糜棱巖，斷層泥等系列，和周圍地質體差異較大，實驗室研究也證明斷層帶填充物對斷層的行為有重要影響［16］。

在弱化帶模型中，對斷層的模擬主要通過改變斷層附近介質物性特征或者改變其介質參數(shù)來模擬近似斷層帶附近易變形的特征。部分學者對斷層填充介質層用不同的本構關系去模擬，例如彈塑性模型［17－18］、應變軟化模型［19－20］、粘彈性模型［21］。

王仁等［17－18］用二維有限單元法研究華北地區(qū)地震遷移規(guī)律時，對斷裂帶的處理是將相鄰的同方向的斷層合并在一起作為介質弱化帶，寬度10km左右?？紤]到斷裂帶內巖石體是經過錯動的，因此將彈性模量E和摩擦系數(shù)μ取低；又由于斷裂帶取得較寬，E，μ的降低量不會很大?？紤]的破裂有兩種，即沿斷層方向的剪破裂和拉破裂。其中采用的剪破裂準則為庫侖破裂準則，當斷層發(fā)生錯動時，斷層內應力這時與應變之間服從理想塑性體法則。模擬結果基本重現(xiàn)出了近十二年來的地震遷移規(guī)律，得到了華北地區(qū)應力場，并為未來地震危險區(qū)提出了參考。

殷有泉等［19］從巖石彈塑性一般本構方程出發(fā)，考慮斷層介質摩擦系數(shù)和內聚力因變形而降低的特點，建立了斷層介質的非穩(wěn)定本構關系，并用一個能量表述的非穩(wěn)定準則來判斷包含斷層和圍巖在內的整個系統(tǒng)的非穩(wěn)定性的到來。在此基礎上模擬出了1976年唐山7.8級地震彈性回跳和1971年圣費爾南多6.6級地震，得到的模擬結果與當時觀測的地震位移數(shù)據(jù)在規(guī)律和數(shù)量級上大體一致。

王繼存等［20］在應變軟化模型基礎上研究了唐山地震斷層破壞及其構造應力場，探討了唐山地震從孕育到發(fā)震的力學過程，分析了唐山發(fā)震斷層破壞的力學機制，并利用增量分析的方法，按時間序列近似地模擬出了發(fā)震斷層和周圍巖石從孕震到發(fā)震的過程。模擬中，唐山地區(qū)斷層帶首先進入塑性階段，并在閉鎖段出現(xiàn)應力集中，相當于唐山斷層處于臨震狀態(tài)；然后在閉鎖段兩側施以方向相反的指定位移，模擬發(fā)震時斷層的突發(fā)走滑，相當于斷層閉鎖段被剪破裂連通而發(fā)震。通過調整介質內摩擦系數(shù)和內聚力的方法來模擬震源區(qū)的應力釋放，模擬得到的塑性區(qū)分布與各時段震余震分布有著較好的吻合。

Mian Liu等［21］用三維粘彈性模型研究了區(qū)域壓應力場背景下青藏高原擴張運動的機理，并通過調整歐亞板塊與印度板塊間的斷裂帶介質粘彈性來近似模擬出兩板塊間的板塊擠壓作用、重力側向力作用機制以及地面速度變化。研究結果表明，若青藏高原高度只為本身高程百分之五十時，則該區(qū)域走滑斷層與逆斷層為主要斷層，且東西向無擴張運動趨勢。當今高原地殼擴張狀態(tài)只有當青藏高原達當今75%高程時才會出現(xiàn)。

陸遠忠等［22］建立了含有母體巖石、硬包體和隨機分布的小裂紋的三維有限元模型，計算了包體和各層實體中的應力分布，模擬強震前巖石的破裂和小震的空間分布特征。結果表明，強震前在孕震體即包體附近均出現(xiàn)了高應力集中單元，它們是形成小震空區(qū)、條帶和地震空間叢集圖像的基礎。隨機裂紋的存在有利于應力在孕震體（包體）外的裂紋端部集中，發(fā)生小震，形成包圍孕震包體的前兆地震活動圖像。而包體中的應力逐漸增加，為發(fā)生強震提供了條件。

陳連旺等［23－24］在研究華北地區(qū)構造應力場及其演化時，根據(jù)華北地區(qū)地殼波速三維結構構建了縱向分層、橫向分區(qū)的三維有限元模型，不同層、分區(qū)內的物性參數(shù)不同，對斷層帶的處理主要是增加其介質泊松比和降低楊氏模量來加以實現(xiàn)。計算結果表明，華北地區(qū)大多數(shù)區(qū)域，最大主壓應力與最小主壓應力（最大主張應力）均位于接近水平的平面內；中等主應力位于接近垂直的平面內即處于走滑應力狀態(tài)。之后，陳連旺等［25］又在上述弱化帶處理的基礎上建立了三維粘彈性有限元模型，同時引入了區(qū)域背景應力場和庫侖破裂準則，模擬了1966年邢臺地震引起的華北地區(qū)應力場動態(tài)演化過程，得到了強地震對華北地區(qū)應力場的震時擾動作用及其引起的一年時間尺度的庫侖破裂應力的動態(tài)變化速率和對下一次強震觸發(fā)作用。

陳化然等［26］建立了包括上地殼、下地殼、上地幔等4層結構的三維有限元模型，計算了川滇地區(qū)背景應力場、斷層蠕動和強震觸發(fā)的應力場及它們的動態(tài)變化。結果表明：川滇地區(qū)后續(xù)地震大多發(fā)生在前面地震引發(fā)的庫侖破裂應力正值區(qū)，前一地震對后續(xù)地震有一定的觸發(fā)作用，強震是在較高的應力背景下成組發(fā)生的。

李平等［27］用三維有限元數(shù)值模擬研究了遼寧及鄰近地殼構造應力場及其與地震活動的關系，模型考慮了地殼上地幔三維地震波速結構、地質構造及地殼深部構造環(huán)境等，其中斷層以寬約1km的條帶模擬，并結合震源機制解與實際地應力測量資料確定了隨深度變化的水平邊界力。研究表明，地質構造、斷層分布與介質力學參數(shù)直接影響到應力場分布，因此對這些參數(shù)研究的精細程度將直接影響到應力場分布的可靠程度。

總體而言，弱化帶模型能很好地描述出了斷層帶內介質的屬性特征以及填充物對斷層行為的影響，且模擬出斷層帶更容易發(fā)生變形的特點。然而其不足之處在于模型計算過程中斷層上單元的共同節(jié)點有且僅有一個方向的位移，因此該模型無法模擬出斷層和周圍介質的交界發(fā)生的位錯，更難模擬出斷層上存在的摩擦機制。

2 劈節(jié)點模型

劈節(jié)點（split node technique）模型由 Melosh等［28－29］在1981年正式提出，是對Jungels等［30－31］擬斷層滑動方法的一種修正和改進。其核心思想是通過引入一個新載荷向量，使得模型中斷層處所對應的節(jié)點在毗鄰的單元上產生大小相同而方向相反的位移即位錯，就像斷層面上的節(jié)點被“劈”開了一樣，斷層則是由這些“劈”節(jié)點所組成。

劈節(jié)點模型主要是通過初始位移的形式將斷層位錯引入有限單元數(shù)值模擬中。以一維的兩個單元為例，以下是該模型的處理思路。單元1，2分別位于斷層的左右兩側，U為位移，上標表示單元號，下標表示單元內的節(jié)點號；ΔU表示位錯量；則有

其中分裂節(jié)點上位移平均值為ˉU12＝ˉU21，ΔU12＝－ΔU21。其示意如圖1。

圖1 一維劈節(jié)點模型示意圖Fig.1 Sketch of One Dimensional Split Node Model.

利用以上關系，可以將單元局部剛度矩陣合成以下全局剛度矩陣

其中Ui為位移；Fi為位移；i為全局節(jié)點號。

Yoshioka等［32］通過橫向異性的三維有限元模型模擬了1993年Kushiro－oki地震的同震和震后位移以及應力場的變化，模擬中斷層位錯量大小通過近場地震波反演來給出，研究認為同震形變在斷層的上下面較大而震后形變主要位于俯沖面與周圍軟流圈附近，且同震狀態(tài)下應力將在未來20年內保持不變。

Hearn等［33］研究了震后余滑作用和粘彈性松弛所導致的震后形變的區(qū)別，以及如何利用這些區(qū)別來辨別出地殼內部及地幔的分布和流變特性。研究認為一個設計良好的GPS流動測點觀測臺網(wǎng)可足夠將震后余滑同線性粘彈性松弛區(qū)分開來，同時還可得到粘彈性層的厚度和粘性參數(shù)。研究指出為最大限度利用GPS觀測資料，部分GPS站點應沿地表破裂帶布設。

Sheu等［34］建立了三維有限元模型對1999年臺灣集集7.6級地震震后形變的機理進行了研究，通過對比分析集集地震后97天內的GPS資料，認為集集地震震后形變是震后余滑和應力松弛共同引起的，其中余滑起主要作用。Buiter等［35］等用二維數(shù)值模型研究了在百萬年時間尺度下，不同俯沖帶板系統(tǒng)（重力拖拽力、板塊運動速度、俯沖帶交界面摩擦系數(shù)等）對板塊接觸邊緣帶垂直位移的影響。結果表明板塊接觸邊緣帶垂直位移主要受俯沖斷層摩擦系數(shù)和傾角參數(shù)等影響，對板塊邊緣垂直影響幅值介于2km到4km之間。

Hu等［36］用Maxwell粘彈性體模型研究了1960年智利地震的震后形變，探討了介質的流變性對地殼形變的影響，其中俯沖板塊，海洋及大陸板塊為彈性介質，大陸及海洋地幔為粘彈性介質。模型中以拋物線函數(shù)擬合了俯沖板塊剖面，由二維延展至三維，以賦予斷層面一個位錯量來模擬地震的發(fā)生。模擬結果認為1960年智利地震35年后GPS觀測得到地面水平位移特征可能是由地幔應力松弛導致的。

Hyodo等［37］進行了日本中部 Niigata－Kobe構造帶上震間應變積累研究。研究表明當板間地震周期比較長時，Niigata－Kobe構造地殼相對較薄的部位存在著較高的應變積累。Zhao等［38］用三維模型探討了理論上三類斷層中地殼分層結構對地表形變和應力的影響，結果表明地殼分層與各項同性地殼結構可導致至少有20%的地表形變差異和15%的應力差異，因此當在進行數(shù)值模擬時必須考慮地殼縱向分層、橫向分塊的結構對模擬結果的影響。

Dixon等［39］研究了加州Baja半島的斷層滑動率受地震輪回和巖石圈流變性的影響，并對當?shù)匦聵嬙旎顒雍偷卣鹞ｋU性做出了評估。模擬結果表明Agua Blanca和San Miguel斷層均在一定程度上受上述兩個因素的影響均比較明顯，并導致約2 mm／yr斷層活動率的變化。此外，Dixon通過計算認為Agua Blanca斷層目前可能處于地震輪回的晚期階段，同時具備著每隔百年孕育6.1至7.0級地震的能力。

Aagaard等［40］利用 Andrews等［41］提出的拉拽劈節(jié)點（TSN）斷層處理方法研究了在同質和分層半空間的情況下走滑和逆沖斷層的動態(tài)斷裂過程和由此引起的地震動，研究表明介質摩擦系數(shù)，剪切模量以及斷層的深度對破裂過程有著一定的影響，且與斷裂的速度和傳播方向有著直接的聯(lián)系。

劈節(jié)點模型簡單有效地以初始位移的形式將斷層近似引入有限單元數(shù)值模擬中，同時通過引入一個載荷向量模擬地震引起的位錯。然而，由于這種模型只是通過模擬斷層的位錯行為來描述斷層的存在，從本質上來說并沒有對客觀原本存在的斷層進行模擬。因此劈節(jié)點模型在地學中的應用主要局限于地震同震和震后效應等相關研究。

3 接觸單元模型

接觸模型主要是通過引入一對接觸面，并在這對接觸面的表面“焊接”上特殊面單元（零厚度的接觸單元）來模擬地質體中斷層的存在和其一定摩擦機制下的滑動特征的。實質是將斷層視為不同地質體之間所接觸的部分，斷層的活動則通過不同地質體間的接觸滑動來模擬。

接觸模型的分析過程比較復雜，有限元方程組無法通過一次求解而實現(xiàn)，而須采取增量求解法對于每個增量必須使用“試探—校核”的迭代方法進行求解。首先根據(jù)前一步結果和本步給定的載荷條件，并通過接觸條件的檢查和搜尋，設定此步第一次迭代求解時的接觸面的區(qū)域和狀態(tài)，然后再根據(jù)上述接觸面區(qū)域及狀態(tài)的假設對于接觸面上的每一點將接觸界面約束條件改為等式約束作為定解條件引入方程組并進行方程的求解，最后利用接觸面上于等式約束對應的不等式約束條件作為校核條件對解的結果進行檢查。如果兩者相符則轉入下一增量步計算，如兩者不相符合則重新進行搜尋及迭代求解，直至滿足校核條件而轉入下一增量步計算。

通過零厚度的單元來模擬非連續(xù)面的行為的思想最早由Goodman等1968年［42］提出，起初為節(jié)理單元法，后經過 Hallquist、Cescotto、Peric等多人［43－45］陸續(xù)調整修改逐漸形成了一套較完善的接觸—碰撞理論，并被廣泛應用于巖土工程、建筑等各個領域，后在地學中也得到了越來越多的應用。

Tom Parsons［46］最早利用含有接觸面單元的三維粘彈性模型模擬了1906年舊金山7.8級地震后圣安地列斯斷裂帶的應力恢復過程，模型中的間斷面代表舊金山斷裂帶上三條的走滑斷裂。斷層面通過零厚度接觸面單元來模擬，接觸單元焊接在周圍粘彈性單元的表面上，斷層可發(fā)生形變且遵循庫倫摩擦準則。此外Parsons等［47］還利用含地表高程，不同厚度地殼層和蠕滑斷層的復雜三維模型并結合GPS觀測資料研究了加利福尼亞地區(qū)的構造應力情況，研究發(fā)現(xiàn)部分地殼運動并非板塊構造所驅動，而是由于地形下降所導致，Garlock地震并非由于高構造應力導致，而是與圣安地列斯斷層的滑動引起的。

Masterlark［48］以1985年墨西哥Jalisco－Colima的Mw8.0俯沖帶地震為例，結合GPS觀測資料，通過對比位錯模型解析解與三維有限元數(shù)值計算結果，研究了位錯模型中均勻各項同性半空間泊松體（HIPSHS）中的每一種假設對地形變的影響。研究表明，各假設對地形變模擬影響敏感度由低到高分別為半空間、泊松體、彈性各向同性和同質體，且除半空間外，各假設所導致的敏感度大小遠大于GPS數(shù)據(jù)的不確定度。

Ellis等［49］用含有摩擦機制的斷層面模擬了新西蘭阿爾卑斯斷裂帶的彈性和非彈性應變積累，地殼體中斷層滑動機制和介質行為由巖石實驗和地球物理測量所確定。研究認為發(fā)生地震的粘滑斷層之下存在著一個蠕動延展帶，且控制著地面水平位移，所模擬得到的地面水平速度場與GPS觀測結果取得了很好的一致性。

Ganas等［50］模擬了萬年尺度下希臘 Hellenic Arc地區(qū)地形變、強地震活動性以及Crete地區(qū)隆起的機理，模型表面由GPS地面水平位移約束。研究表明Aegean地形變、Crete隆起主要由一個33 mm／年快速移動的上沖板塊與一個相對靜止（大概5mm／年）下俯沖板塊導致的，因此Hellenic Arc的構造更接近于是陸地俯沖板塊。

劉峽等［51］在研究華北地區(qū)現(xiàn)今地殼運動及形變動力學數(shù)值模擬時分別采用了三維連續(xù)弱化帶模型和二維滑動摩擦模型。在二維模型中采用了接觸單元法來實現(xiàn)對非連續(xù)變形斷層的模擬，其相對運動由接觸面的正壓力、摩擦系數(shù)和抗摩擦強度來控制，一旦接觸面的剪切力超過摩擦強度或靜摩擦力，則相對運動發(fā)生。研究表明，接觸模型模擬結果與觀測到的GPS數(shù)據(jù)平均離散度高于三維介質弱化模型，這可能是由于接觸模型忽略了華北地區(qū)地殼縱向不均勻性導致的。

李紅等［52］根據(jù)廣州地區(qū)地震地質、地球物理、震源機制解以及GPS觀測等相關資料，用三維有限元模型模擬計算了廣州地區(qū)四條斷裂的現(xiàn)今活動方式以及庫侖破裂應力年累積速率。其中，對廣州地區(qū)的主要活動斷裂主要采用接觸單元表示的非連續(xù)變形的斷層間斷面來模擬，間斷面之間遵從庫侖滑動摩擦準則，模型單元的楊氏模量和泊松比依據(jù)廣州地區(qū)波速結構確定。模擬結果表明，廣州—從化斷裂和瘦狗嶺斷裂東段斷裂面上的庫侖破裂應力年累積速率明顯高于瘦狗嶺斷裂西段和珠江口斷裂西支、東支，數(shù)值模擬計算結果與該地區(qū)地震活動特征基本一致。

王凱英［53］對川滇地區(qū)應力場形成機制及斷層相互作用進行數(shù)值模擬，模型中考慮了主要活動斷裂的存在，并將斷層處理為接觸摩擦面，且斷層帶由具有一定寬度的相對軟弱介質組成。分析表明川滇地區(qū)斷層相互作用現(xiàn)象是塊體非均勻運動過程中應力場調整的反映，是塊體運動的結果。此外王凱英［54］以青藏高原北部地區(qū)幾次強震為例研究了介質非連續(xù)的斷—塊構造模型中地震永久位移造成的區(qū)域斷層靜庫侖應力變化。研究表明在斷塊構造模型中走滑型為主的地震位錯引起的區(qū)域斷層庫侖應力變化在地震破裂面以外隨距離衰減很快，影響范圍僅為幾十公里，和余震分布尺度相當；對青藏高原北部的多次地震的模擬結果顯示，先前地震的斷層位移在后續(xù)地震斷層面所造成的庫侖應力甚微，預示了遠距離地震之間的靜應力觸發(fā)效應不明顯。

接觸模型通過引入一種特殊的零厚度單元來進行模擬斷層存在，斷層被近似為由接觸單元對所組成的接觸面。這種模型最大的優(yōu)點在于可以充分考慮斷層面的摩擦機制對模擬結果的影響。不足之處在于接觸模型的計算分析過程十分復雜，而且接觸剛度的選取困難。為避免接觸面互相穿透而影響計算精度，應選取較大的接觸剛度，然而太大的接觸剛度會產生計算收斂困難，接觸表面互相跳開等。此外接觸模型也沒有考慮斷層帶介質參數(shù)及斷層內填充物（斷層泥）對模擬結果的影響。

4 塊體模型

塊體模型是指將地學問題看成一個多地塊相互作用的動力學過程，各塊體本身同時具有滑移和變形的性質，模型中的斷層實質就是塊體系統(tǒng)中的地塊邊界。通過多塊體作用體系來近似模擬地學現(xiàn)象可能最早起源于活動亞板塊和構造塊體相對運動假說，而這一思想能被成功應用于有限元數(shù)值模擬更大程度上歸功于美籍華人石根華提出的非連續(xù)變形分 析 方 法 （Discontinuous Deformation Analysis）［55－56］。

石根華提出的DDA方法假定每個塊體處處具有常應力，常應變，塊體的運動和變形由12個獨立的變形參數(shù)確定，如下：

［D］＝ ［u0，v0，w0，α0，β0，γ0，εx，εy，εz，γxy，γyz，γzx］T（3）其中 （u0，v0，w0）指塊體質心（x0，y0，z0）的剛體平移；（α0，β0，γ0）指塊體繞x，y，z軸的轉角；（εx，εy，εz，γxy，γyz，γzx）指塊體的3個正應變和3個剪應變。

由于式（3）的常應變假定限制了塊體復雜變形的能力，為更好地求解塊體內部的位移場和應力分布，DDA的塊體單元被離散成mb個四面體單元，內嵌有有限元網(wǎng)格的塊體單元勢能為

式中，［ke］為單元的剛度矩陣；［pe］為單元等效節(jié)點載荷矩陣；［ae］為節(jié)點位移。

對于離散模型，塊體單元間的接觸界面是位移間斷面，界面的滑動和分離是由Mohr－Coulomb準則和無拉應力準則控制的。DDA采用罰函數(shù)求解接觸問題，等價于在接觸點施加剛度很大的法向和切向彈簧，彈簧的加設和移開稱之為開—閉。迭代假設所定義的塊體系統(tǒng)有n個塊體和k個接觸，則塊體系統(tǒng)的總勢能包括各塊體單元勢能及接觸彈簧的應變能：

式中，［kcij］為接觸彈簧矩陣，根據(jù)變分原理δ∏p＝0，經推導可得塊體系統(tǒng)的總體平衡方程：［K］［a］＝［P］，其中［K］為塊體系統(tǒng)整體剛度矩陣，由有限單元剛度矩陣和接觸矩陣集合而成；［P］為塊體系統(tǒng)整體荷載列陣，由有限元等效節(jié)點荷載列陣集合而成。鑒于塊體模型以塊體為研究重要組成部分，因此塊體模型主要研究大尺度的的地學問題。

Liu［57］基于華北地區(qū)的地質構造的斷塊特征，用非連續(xù)變形方法模擬了該地區(qū)的長期形變，包括無震斷層滑動和斷層形變。研究發(fā)現(xiàn)對一條斷層，由相鄰斷塊內塊體旋轉引起的斷層滑動方向可能與區(qū)域構造擠壓引起的滑動方向相反，但由于由構造擠壓引起的滑動量級遠大于由斷塊旋轉引起的斷層滑動，因此，斷層滑動方向基本與一個地區(qū)構造擠壓方向一致。

蔡永恩等［58］用新的 LDDA（Lagrangian Discontinuous Deformation Analysis）方法模擬了唐山斷層破裂、錯動和應力釋放的整個過程。LDDA采用了DDA方法中的接觸判斷準則，采用區(qū)域分解方法求解出接觸面上滿足的約束條件邊界力，再以接觸力為邊界條件求出每個塊體內部的位移和應力。研究發(fā)現(xiàn)唐山地震斷層的破裂速度和應力將與斷層上的初始剪應力有關，且唐山發(fā)震斷層的最大動態(tài)、準靜態(tài)位錯和剪應力降均發(fā)生在中間部位，斷層破裂的傳播速度從震中向東南和西北方向分別為3.08km／s和1.18km／s。

白武明等［59］用DDA＋FEM模型研究了華北地區(qū)各地塊相互制約的構造環(huán)境中發(fā)生1976年唐山地震的動力學過程，并探討了大震引起的華北地區(qū)各地塊特別是鄂爾多斯地塊運動變形及各地塊邊界斷層上應力狀態(tài)變化的特征。研究認為唐山地震的發(fā)生引起了華北地區(qū)各地塊邊界斷層上的應力狀態(tài)發(fā)生一定程度的變化，其中通過唐山震區(qū)的張家口—蓬萊北西走向的斷裂帶及鄂爾多斯地塊北部邊界斷層的主要地帶上剪應力增加而法向應力減小，使得這些斷裂帶的地震危險性增加，這與唐山地震后鄂爾多斯地塊北部邊界斷層發(fā)生一系列6級強震事實基本相符。

張瑞青等［60］通過對1999年岫巖5.4級地震前震及余震分布圖像研究和前人對海城地震的研究，提出了海城、岫巖地震發(fā)震構造模型，并利用DDA＋FEM方法模擬了2次地震釋放的主應力變化、最大剪應力變化等值線圖，地震前后位移變化矢量圖，及發(fā)震斷層滑移隨時間的變化。模擬結果分別與相應地震的震源機制、宏觀等震線、發(fā)震斷層的走向性質等基本一致。

陳祖安等［61］用三維DDA＋FEM 方法以GPS資料做位移速率邊界和震源機制約束，模擬了1997年瑪尼7.9級大震發(fā)生對塊體邊界斷層應力狀態(tài)的影響。此外，陳祖安［62］還研究了2008年汶川8.0級地震孕震機理。研究發(fā)現(xiàn)，龍門山斷裂帶東西兩側地勢、地殼厚度、分層以及物性明顯變化對汶川大地震的孕育和發(fā)生均起了關鍵的作用。其中高原地殼物質向東水平運動，受到龍門山斷裂帶東側介質剛性較大的四川盆地阻擋，使得汶川大地震發(fā)震斷層在大震前已積累了相當水平的應變能，處于不穩(wěn)定狀態(tài)。

塊體模型將剛體位移和塊體變形采用統(tǒng)一的有限元格式求解，所有塊體不僅允許自身有位移和變形，而且允許塊體間有滑動、轉動、張開等運動形式，能夠得到塊體系統(tǒng)的大位移和大變形解。塊體模型結合有限元法和離散元法的優(yōu)點，能夠較好的反映出塊體的運動趨勢，因此在地學中取得了廣泛應用。然而，在塊體模型中對斷層的模擬實質是作為次級塊體的邊界帶，故在二維和三維模型中大多數(shù)斷層往往被簡化為直立斷層，同時也無法考慮塊體內部斷層存在和分布對模擬結果的影響。

5 小結與展望

本文概述了近30年來國內外學者用有限單元法研究地學問題時就具有非連續(xù)性的斷層所采用的幾大主流處理方法，并扼要敘述了在各個斷層模型基礎上展開的地震學和斷層動力學方面的研究。

總體來講，有限單元法關于斷層的處理主要有兩大分流：①20世紀80年代以王仁等人為代表提出的連續(xù)介質弱化帶模型，主要通過改變斷層帶物性特征或其介質參數(shù)來得到斷層帶易變形的特征。這種方法在一定程度上是科學的，因為它考慮到了斷層帶介質與周圍巖體性質的差異性，改變本構方程或介質參數(shù)，使其在一定程度上更容易發(fā)生變形來近似模擬斷層帶附近易變形的特征；②20世紀60年代以來針對處理斷層非連續(xù)特點所提出的非連續(xù)介質間斷面模型。主要代表有1968年Goodman等人提出節(jié)理單元、1981年Melosh提出的劈節(jié)點處理方法、1988年石根華提出的非連續(xù)變形分析方法、接觸單元模型。這些模型的共同點均是將斷層處理為非連續(xù)間斷面，其與弱化帶法相比最直接的優(yōu)點在于能夠更好的模擬出斷層兩側區(qū)域存在的速度和位移場的躍變性以及斷層受到摩擦滑動機制，可相對更真實的表現(xiàn)客觀斷層活動特點，然而，卻很難體現(xiàn)出斷層帶內介質斷層泥等對斷層滑動的影響特征等。

以上斷層處理方法在地學問題模擬分析中得到了廣泛的應用，在對斷層進行處理時，每一類模擬方法也均存在著它自身的優(yōu)勢和局限性。為更準確地認識斷層與地殼形變及地表應力分布關系、深化斷層活動與地震孕育的關系研究，并為未來地震預測提供更加真實可靠的依據(jù)，本文提出今后對斷層的模擬方法可能將會有以下幾個發(fā)展趨勢：

（1）在同一研究區(qū)域中，斷層處理手段的多元化。即根據(jù)不同地質構造選取不同的處理模型，如新發(fā)生的斷層面易用接觸單元來模擬，有斷層泥發(fā)育的斷裂可采用介質弱化模型；根據(jù)不同的研究區(qū)域尺度選取不同的模型，如研究區(qū)域為中國大陸時宜先采用塊體模型進行處理，再對塊體內部的斷層采用接觸或弱化帶模型。

（2）同一斷層面上分層化處理。如將斷層面由淺至深劃分為摩擦和準塑性兩個不同的域，其中摩擦域用接觸模型來模擬，準塑性域則采用應變軟化法來使得巖礦物能夠塑性變形；或可將同一斷層處理為三個構造和流變層組成的斷層模型，上部摩擦層、中部轉換層和下部延性層，以探討各斷層模型對構造應力場和形變場的影響。

（3）過去研究大尺度區(qū)域的地學問題時，相對小的斷裂帶往往都被視為單一平面。然而野外實際調查表明即使是規(guī)模相對較小的斷裂帶其產狀分布也往往不是一簡單平面，斷層的不同部分傾角、形態(tài)、深度均可能會存在較大差別。已有研究表明，斷層產狀和分布的選取在地學中起著重要的角色。因此隨著地球物理勘探技術、地質力學、地球化學和觀測技術的發(fā)展，在對地學問題進行數(shù)值模擬時，斷層產狀、分布及其介質參數(shù)的選取也必將會趨近真實化。

（4）斷層帶物性實測和模擬方法研究。根據(jù)不同條件下的大量實驗結果［16］，斷層泥可大概分為三類：碎屑類、碳酸鹽類、層狀硅酸鹽類。斷層泥的礦物成分又很大程度地影響著斷層的力學行為。因此在對斷層帶進行數(shù)值模擬時，應結合考慮實際的場地構造環(huán)境因素如溫度、圍壓、孔隙壓以及斷層泥礦物成分類型來對斷層的摩擦滑動特征、滑動方式和斷層泥變形機制等加以模擬。

［1］鄧起東.活動斷裂研究的進展與方向［M］.北京：地震出版社，1991.

［2］聞學澤.活動斷裂地震潛勢的定量評價［M］.北京：地震出版社，1996.

［3］薄萬舉.地殼形變與地震預測［M］.北京：地震出版社，2001.

［4］李興才.斷層的震前滑動對唐山地震發(fā)生的影響［J］.地震學報，1992，14（3）：304－308.

［5］周碩愚.斷層形變測量與地震預報［J］.地殼形變與地震，1994，14（4）：90－97.

［6］Sibson R H.Fault rocks and fault mechanisms［J］.J.Geol.，1977，133：191－213.

［7］Strehlau J.A discussion of the depth extent of rupture in large continental earthquakes［J］.Geophys.Monogr.，1986，37：131－145.

［8］Dieterich J H.Earthquake on faults with rate and state－de－pendent friction［J］.Tectonophysics，1992，11：115－134.

［9］Gomberg J.Beeler N M.Earthquake triggering by transient and static deformation［J］.J.Geophys.Res.，1998，103（B10）：24411－24426.

［10］劉桂萍，傅征祥.摩擦狀態(tài)－速率依從的區(qū)域地震觸發(fā)模型研究［J］.地震，2004，24（1）：176－183.

［11］王仁.有限單元等數(shù)值方法在我國地球科學中的應用和發(fā)展［J］.地球物理學報，1994，37（增刊）：128－138.

［12］王妙月，底青云，張美根.地震孕育、發(fā)生、發(fā)展動態(tài)過程的三維有限元數(shù)值模擬［J］.地球物理學報，1999，42（2）：218－227.

［13］王愛國，袁道陽，梁明劍.蘭州盆地最大潛在地震變形數(shù)值模擬［J］.西北地震學報，2008，30（3）：232－238.

［14］章純.中國東部地區(qū)地震活動與構造應力場關系的有限元數(shù)值模擬［J］.西北地震學報，2007，29（3）：230－234.

［15］張彬，楊選輝，陸遠忠.地震動態(tài)應力觸發(fā)研究進展［J］.西北地震學報，2008，30（3）：298－303.

［16］馬勝利.模擬斷層帶摩擦滑動性狀與變形特征［J］.中國地震，1986，2（2）：72－78.

［17］王仁，何國琦，殷有泉.華北地區(qū)地震遷移規(guī)律的數(shù)學模擬［J］.地震學報，1980，2（1）：32－42.

［18］王仁，孫茍英，蔡永恩.華北地區(qū)近700年地震序列的數(shù)學模擬［J］.中國科學（B輯），1982，8：745－753.

［19］殷有泉，張宏.模擬地震的應變軟化的數(shù)學模型地球物理學報［J］.1982，25（5）：414－423.

［20］王繼存，續(xù)春榮.唐山地震斷層破壞及其力學過程的數(shù)值模擬地震地質［J］.1989，11（4）：71－76.

［21］Liu M，Yang Y.Extensional collapse of the Tibetan Plateau：Results of three－dimensional finite element modeling［J］.J.Geophys.Res.，2003，108 （8）：2361.DOI：10.1029／2002JB002248.

［22］陸遠忠，葉金鐸，蔣淳.中國強震前兆地震活動圖像機理的三維數(shù)值模擬研究［J］.地球物理學報，2007，50（2）：499－508.

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

［24］陳連旺，陸遠忠，郭若眉，等.華北地區(qū)斷層運動與三維構造應力場的演化［J］.地震學報，2001，23（4）：249－361.

［25］陳連旺，陸遠忠，劉杰，等.1966年邢臺地震引起的華北地區(qū)應力場動態(tài)演化過程的三維粘彈性模擬［J］.地震學報，2001，23（5）：480－491.

［26］陳化然，陳連旺，馬宏生，等.川滇地區(qū)應力場演化與強震間相互作用的三維有限元模擬［J］.地震學報，2004，26（6）：567－575.

［27］李平，盧良玉，盧造勛，等.遼寧及鄰區(qū)地殼構造應力場及其與地震活動關系的三維有限元數(shù)值模擬研究［J］，地震學報，2001，23（1）：24－35.

［28］Melosh H J，Raefsky A.A simple and efficient method for introducing faults into finite element computations［J］.Bull.Seism.Soc.Am.，1981，71（5）：1391－1400.

［29］Melosh H J，Williams C.A Mechanics of graben formation in crustal rocks：A finite element analysis［J］.J.Geophys.Res.，1989，94（B10）：13961－13973.

［30］Jungels P H.Models of tectonic processes associated with earthquakes［D］.Pasadena：California Institute Technology，California，1973.

［31］Jungels P H，F(xiàn)razier G A.Finite element analysis of the residual displacements for an earthquake rupture：source parameters for the San Fernando earthquake［J］.J.Geophys.Res.1973，78：5062－5083.

［32］Yoshioka S，Tokunaga Y O.Numerical Simulation of Displacement and Stress Fields Associated with the 1993Kushiro－oki，Japan，Earthquake［J］.Pure appl.geophys.，1998，152：443－464.

［33］Hearn E H.What can GPS data tell us about the dynamics of post－seismic deformation［J］.Geophys.J.Int.，2003，155：753－777.

［34］Sheu S Y，Shieh C F.Viscoelastic－afterslip concurrence：a possible mechanism in the early post－seismic deformation of the Mw7.6，1999Chi－Chi（Taiwan）earthquake［J］.Geophys.J.Int.，2004，159（3）：1112－1124.

［35］Buiter J H，Govers Rob，Wortel M J.A modeling study of vertical displacements at convergent plate margins［J］.Geophys.J.Int.2001（147）：415－427.

［36］Hu Y，et al.Three－dimensionalviscoelastic finite element model for postseismic deformation of the great 1960Chilean earthquake［J］.J.Geophys.Res.，2004，109.DOI：10.1029／2004JB003163.

［37］Hyodo M K，Hirahara.A viscoelastic model of interseismic strain concentration in Niigata－Kobe Tectonic Zone of central Japan［J］.Earth Planets Space，2003，55：667－675.

［38］Zhao S，et al.3－D finite element modelling of deformation and stress associated with faulting：effect of inhomogeneous crustal structures［J］.Geophys.J.Int.，2004，157：629－644.

［39］Dixon T J，Decaix F，F(xiàn)arina K.Seismic cycle and rheological effects on estimation of present－day slip rates for the Agua Blanca and San Miguel－Vallecitos faults，northern Baja California，Mexico［J］.J.Geophys.Res.，2002，107（B10），2226，DOI：10.1029／2000JB000099.

［40］Aagaard B T，Hall J F，Heaton T H.Characterization of near－source ground motions with earthquake simulations［J］.Earthquake Spectra，2001，17（2）：177－207.

［41］Andrews D J.Test of two methods for faulting in finite－difference calculations［J］.Bull.Seism.Soc.Am.，1999，89：931－937.

［42］Goodman R E，Taylor R L，Brekke T L.A model for the mechanics of jointed rock［J］.Soil Mech.and Found Proc，1968，94：637－658.

［43］Hallquist D，John O.LS－DYNA3DTheoretical Manual［M］.Livermore：Software Technology Corporation，1994.

［44］Cescotto S，Charilier R.Frictional Contact Finite Elements Based on Mixed Variational Principles［J］.International Journal for Numercial Method in Engineering，1992，（36）：1681－1701.

［45］Peric D，Owen D R.Computational Model for 3－D Contact Problems with Friction Based on the Penalty Method［J］.International Journal for Numercial Method in Engineering，1992，35：1289－1309.

［46］Parsons T.Post－1906stress recovery of the San Andreas fault system calculated from three－dimensional finite element analysis［J］.J.Geophys.Res.，2002，107（B8）：2162，DOI：10.1029／2001JB001051.

［47］Parsons T.Tectonic stressing in California modeled from GPS observations［J］.J.Geophys.Res.，2006，111（B03）407，DOI：10.1029／2005JB003946.

［48］Masterlark T.Finite element model predictions of static deformation from dislocation sources in a subduction zone：Sensitivities to homogeneous，isotropic，Poisson－solid，and halfspace assumptions［J］.J.Geophys.Res.，2003，108（B11）：2540，DOI：10.1029／2002JB002296，.

［49］Ellis S，et al.Simplified models of theAlpine Fault seismic cycle：Stress transfer in the mid－crust［J］.Geophy.J.Int.，2006，DOI：10.1111／j.1365－246X.2006.02917.x.

［50］Ganas A，Parsons T.Three－dimensional model of Hellenic Arc deformation and origin of the Cretan uplift［J］.J.Geophys.Res.，2009，114 （B06）：404，DOI：10.1029／2008JB005599.

［51］劉峽，傅容珊，楊國華.用GPS資料研究華北地區(qū)形變場和構造應力場［J］.大地測量與地球動力學，2006，26（3）：33－39.

［52］李紅，陳連旺，李紅江.廣州地區(qū)活動斷裂的數(shù)值模擬［J］.大地測量與地球動力學，2008，28（2）：39－44.

［53］王凱英，馬瑾.川滇地區(qū)斷層相互作用的地震活動證據(jù)及有限元模擬［J］.地震地質，2004，26（2）：259－272.

［54］王凱英.斷塊模型中走滑型地震應力觸發(fā)研究——以青藏高原北部幾次強震為例［J］.地球物理學報，2009，52（7）：1776－1781.

［55］Shi G H.Discontinuous deform ation analysis：a new numerical model for the statics and dynamics of block systems［D］.Berkeley：Department of Cilvil Engineering，University of California，1988.

［56］白武明，林邦慧，陳祖安.1976年唐山大震發(fā)生對華北地區(qū)各地塊運動與變形影響的數(shù)值模擬研究［J］.中國科學（D輯），2003，33（增刊）：99－107.

［57］Liu L.Modeling aseismic fault slips and block deformation in northern china by DDA［A］∥Proc of the First International Forumon Discontinuous Deformation Analysis（DDA）and Simulations of Discontinuous Media［C］.Albuquerque：TSI Press，1996：373－383.

［58］自武明，林邦慧，陳祖安.1976年唐山大震發(fā)生對華北地區(qū)各地塊運動與變形影響的數(shù)值模擬研究［J］.中國科學（D輯），2003，33（增刊）：99－107.

［59］蔡永恩，何濤，王 仁.1976年唐山地震震源動力過程的數(shù)值模擬［J］.地震學報，1999，21（5）：469－477.

［60］張瑞青，魏富勝，喬成斌.用（DDA＋FEM）方法數(shù)值模擬1975年海城、1999年岫巖地震發(fā)生的過程［J］.地震學報，2005，27（3）：163－170.

［61］陳祖安，林邦慧，白武明，等.1997年瑪尼地震對青藏川滇地區(qū)構造塊體系統(tǒng)穩(wěn)定性影響的三維DDA＋FEM方法數(shù)值模擬［J］.地球物理學報，2008，51（5）：1422－1430.

［62］陳祖安，林邦慧，白武明，等.2008年汶川8.0級地震孕震機理研究［J］.地球物理學報，2009，52（2）：408－417.