基于蝶形破壞理論的地震能量來源

馬 驥，趙志強，師皓宇,3，郭曉菲，喬建永,2，馬念杰

(1.中國礦業大學(北京) 能源與礦業學院，北京 100083; 2.北京郵電大學 理學院，北京 100876; 3.華北科技學院 安全工程學院，北京 101601)

地震是地殼板塊運動時快速釋放能量過程中造成的急劇震動，并會產生地震波的一種自然現象。地震過程的第1階段，即在地球內部能量激發的起因，至今仍是眾說紛紜，并沒有形成統一的認識。自20世紀初，科學家們就地震的能量來源，相繼提出彈性回跳、相變和巖漿沖擊3個比較有影響的假說[1-2]。彈性回跳假說詳細并系統地總結了地震的成因與能量的來源，即地震是地殼巖體突然斷裂錯動引起的，接著巖體又沿著斷裂面整體彈性回跳到初始狀態，地震能量來源于斷層兩側巖體發生不均勻彈性變形所產生和積累的彈性變形能[3];相變說揭示的能量來源于地下物質由于臨界溫度和壓力作用下，其自身密度增大體積突然變小，發生相變，周圍巖體由于應力狀態發生改變，快速擠壓該相變物質從而激發地震波[4];巖漿沖擊說認為地震能量是地下深處高溫高壓的巖漿涌入地殼，使得地殼巖體導熱不均，部分巖體體積膨脹擠壓周圍巖體，導致周圍巖體發生破裂激發的地震波[5]。

受三大經典假說影響，后來學者針對不同類型地震分別提出了一系列假說，試圖揭示地震能量來源。黏滑說基于彈性回跳假說將斷層變形問題轉化為斷裂兩側巖體的摩擦問題，指出地球上的淺層地震(深度<70 km)為新老斷層滑動過程中的黏滑，地震能量來源于相鄰震區提供給震源體的應變能[6-8];既然因斷層運動導致巖石破裂就能引發地震，研究巖石本身破裂機制十分必要。巖石破裂說提出斷層面上存在障礙體與凹凸體破裂模型，地震發生能量來源于積累高能量載體自身突然破裂[9-10];剪切熔融假說認為地殼巖體在剪切應力持續作用下發生蠕變，在蠕變加速過程中巖體變形并集中在一個具有高速蠕變特征的薄層中，該薄層由于溫度升高產生熔融進而導致沿剪切帶滑移，最終由于耗散產生剪切失穩發生地震(深度>300 km)，此過程描述了地震的能量主要來源于高速蠕變的薄層溫度升高使得局部巖體弱化存儲的能量[11];脫水脆裂假說解釋板塊在俯沖過程中將一定數量的水帶入地球深部，隨深度增加溫度和壓力不斷增大，當達到一定的溫壓條件時引起礦物脫水脆化破裂，以彈性波的形式釋放積聚的彈性變形能進而引發地震(70～300 km)[12];反裂隙斷層假說基于相變說的論述，指出地震能量來源于超高壓作用下，深源巖體內部反向裂隙貫通形成斷層時以彈性波形成釋放的彈性應變能[13];孕震斷層多鎖固段脆性破裂理論認為孕震斷層存在一個或多個鎖固段，強度高的鎖固段承受應力集中，是高能量的載體，發震能量主要集中于此;某個鎖固段發生宏觀破裂后，應力向下一個鎖固段轉移，導致下一個鎖固段承受應力集中，以此類推[14];重力塌陷假說解釋地震能量來源于震源巖體直接釋放出的重力位能，能量的大小取決于作用在震源巖體上的重力所產生的構造應力[15];地球自轉速率的加強或減慢觸發板塊構造運動，使得板塊邊緣與構造聚集帶地震活躍，揭示了地震能量與地球自轉有關系[16];共振假說指出地震能量主要來源于震源區巖體自身的自由振動周期與天體、太陽以及地球自轉產生的某些周期相一致或接近整數倍時產生共振做功所積累的動能[17]。這些假說對人們從宏觀上認識地震的能量來源提供了良好借鑒。

近期出現的蝶形破壞理論[18-20]可以從微觀角度認識地震的發生機理和能量來源，它較好的解釋了煤礦巷道沖擊地壓發生機理[21-22]與掘進巷道煤與瓦斯突出的機理[23]，為地震能量來源的研究提供了一種新方法。筆者首先從軟弱異性體周圍巖體發生塑性破壞引發能量改變角度出發，構建地震發生時以軟弱異性體為中心的能量分析模型，提出地震震源能與地震能的計算方法;然后，采用FLAC3D數值模擬手段理論分析地震發生過程中能量積聚與釋放的變化特征;最后，對比分析真實地殼應力環境下有無軟弱異性體存在，地震震源能分布特征與觸發地震能的變化特征，闡明不同類型(形狀、尺寸與性質)軟弱異性體存在所揭示規律的一致性，及能量分析模型與計算方法的合理性與適用性。

1 地震能量來源的理論分析

1.1 均勻彈塑性介質中軟弱異性體圍巖的能量分析模型

地震時會釋放大量的能量，研究地震發生機理，首先需要弄清地震的能量來源。受蝶形破壞理論的啟發[18-20]，考慮從軟弱異性體周圍能量集中現象進行分析，如圖1所示，這里的軟弱異性體是指真實地殼巖體中存在著的強度較低、具有一定尺寸，任意形態特征，可能充斥有氣態、液態物質的軟弱巖體、破碎固體等，是一個廣義的概念。

圖1 自然界中真實存在的軟弱異性體Fig.1 Soft anisotropic bodies in nature

由于軟弱異性體的存在，周圍巖體應力將重新分布，形成圍繞軟弱異性體周邊的應力集中。應力的變化也帶來了軟弱異性體周圍巖體中能量分布的變化，即能量集中。由彈性力學理論可知[24-25]，集中在軟弱異性體周圍巖體內的彈性能可以通過原巖應力和周圍巖體的物理力學參數計算出來。為此，構建了軟弱異性體存在發生塑性破壞集中與釋放能量的計算分析模型如圖2所示。

圖2 軟弱異性體周圍能量計算分析模型Fig.2 Computational,and analytical,model for energies around a soft anisotropic body

1.2 含軟弱異性體圍巖集中彈性能的計算

板塊構造學說認為:劃分板塊的巖石圈具有較高的剛性和彈性(均厚70～80 km)，漂浮在低密度、塑性軟流圈之上作大規模運動[26]。因此，對地震過程有意義的彈性應變能可等效為引起巖石圈體積發生變化的體變能。當應力和應變滿足線性關系時，巖石圈任意微單元體均滿足虛功原理[27]，則微單元單位體積應變能，即應變能密度可表示為

(1)

其中，σ1,σ2,σ3分別為區域應力場主導下的最大，中間與最小主應力；E為均質巖體的彈性模量；μ為巖體的泊松比。如圖2(a)所示，假設大尺度的地質體為理想彈性體，則其空間內任意一點的應變能密度可由式(1)計算(軟弱異性體及其周圍巖體構成空間閉區域任意單元的應變能密度)。當含軟弱異性體圍巖出現破壞后，其塑性破壞區的微小單元體仍看作處于彈性狀態下，則該區域彈性應變能密度可等效用式(2)計算(此時的彈性模量與泊松比為塑性狀態下的量值)。

σp2σp3+σp1σp3)]

(2)

式中，σp1,σp2,σp3分別為塑性破壞區域微小單元受到的最大，中間與最小主應力；Ep為該區域巖體的彈性模量；μp為塑性破壞區巖體的泊松比。軟弱異形體與其圍巖共同構成空間閉區域Ω(規則圓形軟弱異形體半徑a的5倍以上范圍)，當軟弱異性體圍巖沒有出現任何破壞時其彈性應變能最大，表示為

(3)

塑性破壞區出現后，必定消耗能量，這就意味著圖2的能量計算模型存在著彈性Ωe和塑性Ωp(應力極限平衡)兩種不同的區域，其中Ωe∪Ωp=Ω。則在軟弱異性體周圍部分巖體出現破壞后，含軟弱異性巖體的全部彈性應變能用Uep表示，即

(4)

式中，u′e，up分別為彈性區Ωe與塑性區Ωp的微單元單位體積應變能密度；dVe，dVp為對應微單元的單位體積。

含軟弱異性體的最大彈性應變能Uemax大于軟弱異性及其周圍巖體破壞后的彈性應變能Uep，這一差值伴隨著軟弱異性體圍巖的破壞而消失。其中部分用于巖石內結晶晶格錯位，部分產生熱量，還有一部分會引起巖體的震動并以地震波的形式傳播出去。定義這種伴隨軟弱異性體圍巖破壞而引發巖體震動，并以地震波的形式傳播出去的能量為震源能We，其數學表達式為

We[Px,y,z(t)]=β(Uemax-Uep)

(5)

其中，We為震源能，106J；β為震動能因子，0<β<1;Px,y,z(t)為以軟弱異性體為中心形成的整個空間閉區域巖體受到3個方向的主應力;We為一個以Px,y,z(t)為自變量的復合函數。

1.3 基于能量分析模型的數值分析

板塊構造學說的提出，合理解釋了地殼運動相對劇烈的板塊邊界附近是地震活動的頻發區域[28]:由于板塊之間的相互運動，必然造成板塊巖體應力的變化，通過“世界應力圖”發現全球大部分地區的最大水平主應力方向與板塊絕對運動跡線保持較好的一致性，反映出構造應力與板塊運動的關系密切[29];被人們稱之為“世界屋脊”的青藏高原是印度板塊與歐亞板塊間陸陸碰撞的結果，如今印度次大陸仍向北運動，使得青藏高原受南北向巨大擠壓應力場作用，表現為強烈擠壓沖斷,大規模的走滑與剪切、正斷與拉伸等構造特征[30];有“地球上的傷疤”之稱的東非大裂谷是由于非洲板塊SW向運動和印度洋板塊NE向運動拉伸張裂形成的，地殼下面的地幔物質上升分流，產生巨大的張力，正是在這種張力的作用之下，地殼發生大斷裂，從而形成裂谷[31];總之，板塊運動會產生巨大的壓力與張力(構造應力場變化)，其應力方向與最大水平主應力方向趨于一致，由于地殼局部區域較強的剛性，致使彈性應變能積累，當這個應變能積累到超過局部巖體強度極限時勢必產生塑性破壞，并伴隨能量釋放。

圖3 含軟弱異性體圍巖體模型Fig.3 Model of the rocks surrounding a soft anisotropic body

數值模擬方法是近似直觀反演地震發生物理過程的有效手段之一，本文應用FLAC3D數值模擬構造應力場變化產生局部塑性破壞，伴隨能量釋放的變化過程。構建以軟弱異性體為中心的圍巖體數值模型，如圖3所示?；谑ゾS南原理[32]，在小邊界上進行面力的靜力等效變換后，只影響附近局部區域的應力，對絕大部分彈性體區域的應力沒有明顯影響。因此假設軟弱異性體(充斥著各種氣態物質的圓形孔洞構造)為規則的圓形，半徑為10 m，長度為10 km，位于地下深處10 km，取軟弱異性體周圍500 m半徑區域范圍內的巖體為均質變質巖或火成巖(參考了花崗巖的力學參數取值)，具體巖石力學參數見表1。模型x軸、y軸邊界水平位移與z軸上下邊界垂直位移均固定。

表1 均質巖體物理力學參數Table 1 Physical and mechanical parameters of homogeneous rock masses

設軟弱異性體及其周圍巖體處于相對穩定的等壓應力狀態下，即Px(t)=Py(t)=Pz(t)=270 MPa，區域圍壓均為時間t的函數;結合彈塑性力學平面應變問題解法，假設垂直向主應力為軟弱異性體上覆巖層自重，即Pz(t)=270 MPa，則水平向區域主應力Px(t)為惟一自變量。

1.3.1軟弱異性體周圍蘊含能量的變化規律

為了揭示軟弱異性體周圍蘊含能量的變化規律，建立起區域水平主應力與震源能、里氏震級間的對應關系，如圖4所示，軟弱異性體周圍蘊含能量可以通過式(5)(這里β取1)定量計算得出。以區域水平主應力為自變量，震源能與里氏震級作為因變量表現出幾乎一致的特征:隨著水平主應力Px的改變(以Px=Pz=270 MPa為中心，曲線左側水平主應力減小，右側主應力增大)，震源能與地震震級在曲線的兩側出現了近似指數型增長的變化規律。軟弱異性體周圍積聚的能量最大可達1.80×1016J，如果這部分能量一下釋放，則相當于里氏7.6級地震。以軟弱異性體及其周圍巖體的等壓狀態為中心，隨著水平主應力Px逐漸增加，使得整個巖體系統的總能量增加，所以軟弱異性體周圍巖體的能量也增大，可能引發地震的級別也越大;然而圖4左側函數變化關系顯示:隨著水平主應力Px的減小，系統能量總體減小的情況下軟弱異性體周圍巖體的能量反而增大，對應地震的級別也是越來越大。正是這種隨著系統總能量的減小地震級別越來越高的反常現象，很好的解釋了受構造應力作用強烈的區域附近往往是地震多發帶的自然規律，即構造地應力場的劇烈變化，導致地球內部巖體突然的破裂，從而引發不同級別的地震;于是，以等壓力狀態為中心，將水平主應力減小軟弱異性體周圍巖體彈性能增長的區域劃分為張拉破壞區TFZ，水平主應力增加能量增加的區域劃分為擠壓區CFZ(圖4)。依據軟弱異性體周圍巖體塑性破壞形態的變化特征(圖5,6)，將擠壓應力區地震活動分為3個階段(圖4):I為軟弱異性體周圍巖體塑性區為圓形;II為蝶形塑性區形成初期;III為蝶形塑性區發展末期(張拉應力區具有同樣特征，本文僅以擠壓應力區能量變化特征為例進行分析)。

圖4 軟弱異性體蘊含震源能、里氏震級與水平主應力關系Fig.4 Relationship between SEP-SS in rock masses surrounding the anisotropic body,Richter magnitude,and the horizontal principal stress

1.3.2軟弱異性體圍巖能量的分布特征

圖5,6提取出不同區域主應力Px作用下，塑性破壞區形態變化與對應震源能密度分布特征。圖6中A,B,C分別對應圖4劃分的階段(I，II，III)中任意取得的3個特征點。分析可得軟弱異性體圍巖形成的蝶形破壞區周圍集中了大量的彈性能，是地震發生時能量的主要來源。以擠壓應力區為例，原巖應力Px=270 MPa，塑性區的形狀為圓形。地質構造劇烈變化，水平應力開始增加，當Px=810 MPa時(Px/Pz=3)，軟弱異性體周圍的破壞產生質的轉變，呈蝶形變化特征。也就是說，滿足條件270 MPa

圖5 張拉破壞區的蝶形的演化過程與能量分布Fig.5 Evolutionary process and energy distribution of butterfly failure in TFZ

圖6 擠壓破壞區的蝶形的演化過程與能量分布Fig.6 Evolutionary process and energy distribution of butterfly failure in CFZ

1.3.3 自然地震的觸發

(1)地震能與震源能的關系

依據前文定義:震源能是指某種破壞狀態下軟弱異性體周圍巖體內的集中能量，是某一時刻的狀態量值，與時間無關。而真實地震發生時記錄到的以地震波形成釋放的能量是以時間為周期加以衡量，例如里氏震級(ML)、面波震級(Ms)及體波震級(Mb，MB分別為用1 s和5 s左右的地震體波振幅來量度地震的大小)中能量的計算都是以波的周期長度作為時間單位。為了方便起見，本文定義單位時間釋放的震源能為地震能Wm，如圖2(b)所示。結合式(5)構建地震能Wm與震源能We的關系:

Wm=We[Px,y,z(t+Δt)]-We[Px,y,z(t)]=

We(Px,y,z+ΔPx,y,z)-We(Px,y,z)

(6)

式中，Wm為地震能，即地震發生時每秒釋放的震源能，J/s;t為時間變量，s;Px,y,z(t)為地震發生前時間t的區域主應力，MPa;Px,y,z(t+Δt)為地震發生時的區域主應力，MPa;ΔPx,y,z為經歷Δt時刻的區域主應力增量，MPa/s;D為常量；R為全體實數域。

式(6)表明，對于軟弱異性體而言，發生地震時的震級和釋放的能量僅取決于區域主應力Px,y,z(t)及其增量ΔPx,y,z。

(2)區域主應力微小變化觸發地震的可能性分析

與地球提供給板塊運動的能量相比，地震所釋放的能量僅僅占很小的部分[33];地球板塊運動仿佛一直處在不穩定的邊緣，而地震似乎可以看成是圍繞這一“臨界狀態”的“漲落”[34]。固體潮汐應力、地球自由震蕩、已有斷層的突然破裂、采礦活動等微小應力亦可“觸發”大級別地震的自然現象不斷被發現，甚至極其微小的遠處大地震的面波通過(引起的水平張應力僅為0.025 MPa)也能引發地震[35-39]。由此看來，“觸發”地震似乎無需非常大的應力變化，這一說法令人難以置信。

這種微小應力變化“觸發”地震的自然現象可通過式(6)的計算結果很好的定量解釋。假設軟弱異性體位于擠壓構造應力場中，則在受到1，0.1與0.01 MPa/s的觸發應力作用時，不同震前狀態(圖4，A點540 MPa;B點1 180 MPa;C點1 572 MPa)下的地震能和震級如圖7所示，受到相同觸發應力增量作用時，圍巖塑性破壞形態呈圓形與蝶形塑性破壞形成初期(特征點A+B)的地震能和地震級較小，蝶形塑性破壞發展末期(特征點C)最大，震級可達6.6級。此時，水平主應力僅增加0.01 MPa/s，即可觸發里氏4.3級左右地震，說明地震誘因存在著“壓倒駱駝的最后一根稻草”現象;處于蝶形塑性破壞急劇擴展階段，地震對地應力分布具有敏感依賴性，區域主應力的任何微小增量ΔPx,y,z都有可能引發能量系統的災變，可以推測微小、甚至極微小的應力變化都可“觸發”較大級別的地震。這里需要指出:① 由于應力是矢量，以上表述的觸發應力方向是與區域水平主應力的方向保持一致的。如果二者的作用方向不一致，其結果也會完全不一樣，尤其當觸發應力和區域水平主應力的作用方向完全相反時，即使再大的應力變化也不會“觸發”地震;② 區域主應力(x,y,z方向上)盡管有時存在相互作用關系，但它們之間也可看成獨立的自變量，即以上所有過程基于固定區域垂直主應力，僅分析水平主應力變化引發地震的自然現象，也可用于解釋不同方向上區域主應力變化引起的地震能量改變。

圖7 應力增量對應能量與震級變化Fig.7 Changes in energies and magnitudes corresponding to different stress increments

2 震例分析

2.1 區域構造背景及有限差分模型

自然界中的軟弱異性體形形色色，有各種類型。軟弱異性體極少是圓的，存在各種各樣的形狀;軟弱異性體很多情況下會是被堅硬巖體包圍著的局部相對軟弱的地質軟弱構造。

中國大陸西南地區是印度板塊與歐亞板塊碰撞形成的“擠出構造”，地處南北地震帶中南部的川滇菱形構造塊體及其邊緣區域，構造應力環境極其復雜。研究區域具體的空間位置為:26°～28°N，103°～104°E，其中包含了魯甸6.5級(2014年中國云南省)地震震中所在位置。考慮到魯甸震源區地殼厚度變化較大，從44.5 km增加到59.0 km的特點[40]，采用基于拉格朗日連續介質法的FLAC3D5.0數值模擬軟件[41]，構建了真實地殼的三維有限差分模型，如圖8所示，東西向為模型的走向，其長度為60 km，高度取為30 km，厚度為100 m(南北向表示厚度方向)。為了方便模型的構建和突出所研究的中心區域，建模中進行了必要的簡化，由于本文主要研究地震震源能量來源，對于震中地表及走向范圍內的地勢差異影響未作考慮;由于模型所建高度未達到巖石圈底部Moho面，故真實巖石圈底部的起伏變化也加以忽略。在距模型地表12 km深度處設置一軟弱異性體(長方體單元)(圖8)，其尺寸為500 m×500 m×100 m(x×z×y)。

模型巖石力學參數選取:依據摩爾-庫侖破壞準則，涉及彈塑性介質的5個物理力學參數分別為:彈性模量E、泊松比μ、黏聚力c、內摩擦角φ和抗拉強度σc。由PREM地球模型得出地震P波、S波在地球內部傳播速度變化特征:隨深度增加，地震波波速增加，巖石更加致密，地殼巖體的彈性模量也同等增加。在30 km深度中國大陸西南地區的平均S波波速約為3.74 km/s[43]，平均P波波速為6.45 km/s[44-45]，隨深度變化的密度取值在2.35～2.85 g/cm3內呈梯度增長[46]。于是，E與μ依據式(7)可求出[47]，即

(7)

式中，α與β分別為P波和S波的傳播速度，km/s;ρ為已知巖體的密度，kg/m3。

將計算所得彈性模量按梯度均勻賦值于模型，同時參照花崗巖的巖石力學性質[48]，確定未給出的物理力學參數，具體取值見表2(軟弱異性體參照煤的物理力學參數取值)。

為了方便對比分析，將已構建的含軟弱異性體模型定義為模型1。同時建立無軟弱異性體模型2，具體參數取值參照表2。用六面體單元對2個模型進行網格劃分，對于模型1進行局部網格加密處理，單元總數為246 800個，節點總數為274 349個。

圖8 西南地區構造簡圖[42]Fig.8 Geological structures of southwest China[42]

表2 模型物理力學參數Table 2 Physical and mechanical parameters of the model

2.2 邊界條件與初始條件

板塊運動觀測網絡給出了中國大陸地區地表GPS速度場結果，模型所在區域(圖8)內部魯甸地震臨近地區的GPS監測資料真實反映了包谷垴—小河斷裂帶(1991—2013)東西兩側運動存在明顯差異，即西側運動量值約為10 mm/a左右，東側為6 mm/a[49-50]。本文所建模型的左側邊界相對于右側可以移動，移動速度取2個邊界實際移動速度的差值，并將該相對速度值插值于模型左側邊界面，作為水平向板塊運動位移移進量，且假定從地表到30 km深度保持一致[51-52]，速度取值為4 mm/a(每年4 mm)，在FLAC3D數值模擬軟件中表現為1 a(年)相當于1個step;垂直方向位移保持自由;右側界面水平向固定，等同于板塊運動相對靜止的邊界面，其他方向自由;上表面為自由邊界，即法向應力和剪應力均為0，對于底部邊界，將底面垂直方向位移約束為0，而水平方向自由,具體如圖9所示。整個模型范圍內巖體初始應力環境設定為自重應力場，鉛直應力變化隨深度成正比遞增，由此構建了模型1與模型2的初始應力環境。

圖9 約束條件與初始位移速度加載模型Fig.9 Constraint conditions and loading model for initial displacement speed

2.3 數值模擬結果分析

2.3.1在有無軟弱異性體條件下，地震震源能變化規律的對比分析

軟弱異性體總體積只占整個地殼模型總體積的14‰，為了便于分析軟弱異性體周圍巖體能量的分布規律，選取以上兩數值模型單元體ID=19 902 071的區域水平應力為自變量，10 000 m×5 000 m×100 m(x×z×y)內(圖8)蘊含能量值為因變量繪制曲線圖，具體如圖10所示，能量值由式(5)計算得出。在初始自重應力場影響下，隨著模型邊界位移移進量逐漸增加，該區域水平應力Px(t)逐漸增大，與因變量震源能We呈正指數型增長關系;當區域水平應力達到某一極限值Pxmax=1.414 GPa時，軟弱異性體集中能量具體量值可達2.25×1016J。為了對比分析相同力學條件下無軟弱體圍巖能量分布特征，建立數值分析模型2。顯然，無軟弱異性體模型2集中震源能We與區域水平主應力之間呈線性增長關系;考慮塑性區形態的變化特征與積聚能量的分布特征，將能量變化曲線分為3個階段:蝶形塑性破壞還未形成0～150 000 a(僅軟弱異性體圍巖發生塑性破壞)、蝶形塑性破壞形成初期150 000～320 000 a(蝶形塑性破壞形成，并穩速擴展)與蝶形塑性破壞擴展末期320 000～370 000 a(蝶形塑性破壞急速擴展)。

2.3.2不同應力狀態下震源能的分布特征分析

從圖10劃分的3個階段中任意選取特征點a，b，c，d和a′，b′，c′，d′，繪制塑性區形態與集中能量分布對比圖11。由該圖可知，模型1積聚震源能總體分布特征表現為軟弱異性體中心集中能量最多，并由里向外逐級減弱，尤其隨著蝶形塑性區的發展，軟弱異性體圍巖會形成蝶形能量集中區(圖11中32萬a首次生成)，蝶葉周圍巖體也積聚大量彈性能;無軟弱異性體模型2能量呈層狀分布，其積聚能量分布與線性遞增的曲線能量值呈正相關，符合客觀規律。

圖10 塑性破壞區與震源能分布Fig.10 Distribution of plastic damage zones and seismic source’s energy

圖11 特征點集中能量分布對比Fig.11 Comparison of distributions of energies accumulated at feature points

以上對比分析證實:構造應力主導下的真實地殼巖體，由于軟弱異性體的存在造成了其周圍巖體的能量集中，隨著圍巖蝶形破壞的發展能量逐漸呈蝶形分布特征，且隨著蝶葉的擴展其周圍巖體持續破壞并伴隨能量釋放，發生一系列地震。

Pregnant Period:蝶形塑性破壞還未形成時期;Growth Period:蝶形塑性破壞形成初期;Upheaval Period:蝶形塑性破壞擴展末期。

2.3.3不同圍巖塑性區狀態下地震的觸發

選取模型1中不同塑性區狀態下(a,b,c三點)集中震源能為初始值，取微小應力增量ΔP=0.001 MPa，則地震能對應里氏震級變化關系見表3，地震能由式(6)定量計算得到。處于蝶形塑性區形成初期與穩速發展時期，微小應力可能觸發里氏2.3和2.5級地震，急速擴展末期的蝶形塑性區微小應力會引起碟葉發生急劇擴展，此時釋放地震能可達6.02×1010J，相當于里氏4.0級地震，說明地震發生時，真實地殼圍巖系統對于軟弱異性體周圍塑性區狀態具有敏感依賴性，蝶形塑性區一旦出現，系統由穩態向非穩態能量積聚轉變，處于劇變期的圍巖系統本身積聚震源能量值已經很大，任何微小的擾動即可引發能量釋放，蝶形塑性破壞發生災變產生大級別地震。

表3 微小力觸發地震能與震級變化關系Table 3 Relationship of seismic energies triggered by microstress and corresponding magnitudes

4 結 論

(1)得到了一種計算地震震源圍巖能量的方法，發現了地震震源能、相應里氏震級與水平主應力關系曲線呈指數型變化規律，即以等壓狀態為中心，隨著水平主應力減小(形成張拉破壞區)，震源能與對應震級呈負指數型增長;隨著水平主應力增大(形成擠壓破壞區)，震源能與對應震級呈正指數型增長，且蝶形破壞圍巖形成了以軟弱異性體為中心的蝶形能量集中區，蝶葉周圍巖體集中了大量的彈性能。

(2)揭示地震能量來源的力學本質。由于地殼內部軟弱異性體的存在，使得周圍巖體應力重新分布，形成圍繞軟弱異性體的應力集中。應力的改變也帶來了軟弱異性體周圍巖體能量的集中，由于構造應力場變化而產生的巨大壓力與張力，造成軟弱異性體周圍巖體出現蝶形破壞，積聚的彈性能得以釋放，形成以軟弱異性體圍巖蝶形破壞區為震源的地震。也就是說，地震的能量主要來源于軟弱異性體周圍巖體內的能量集中，并以此構建了地震能量的分析模型。

(3)應用數值分析方法初步論證了微小應力擾動可以引發大級別自然地震的觀點。實際上不同的應力變化可引發不同級別的地震，之所以同等應力變化下地震釋放的能量大小存在巨大差別，是由地震發生時塑性破壞范圍與圍巖積聚能量所處的狀態決定的。