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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

σp2σp3+σp1σp3)]

(2)

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

(3)

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

(4)

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

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

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

(5)

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

1.3 基于能量分析模型的數(shù)值分析

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

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

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

表1 均質(zhì)巖體物理力學(xué)參數(shù)Table 1 Physical and mechanical parameters of homogeneous rock masses

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

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

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

圖4 軟弱異性體蘊含震源能、里氏震級與水平主應(yīng)力關(guān)系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提取出不同區(qū)域主應(yīng)力Px作用下，塑性破壞區(qū)形態(tài)變化與對應(yīng)震源能密度分布特征。圖6中A,B,C分別對應(yīng)圖4劃分的階段(I，II，III)中任意取得的3個特征點。分析可得軟弱異性體圍巖形成的蝶形破壞區(qū)周圍集中了大量的彈性能，是地震發(fā)生時能量的主要來源。以擠壓應(yīng)力區(qū)為例，原巖應(yīng)力Px=270 MPa，塑性區(qū)的形狀為圓形。地質(zhì)構(gòu)造劇烈變化，水平應(yīng)力開始增加，當(dāng)Px=810 MPa時(Px/Pz=3)，軟弱異性體周圍的破壞產(chǎn)生質(zhì)的轉(zhuǎn)變，呈蝶形變化特征。也就是說，滿足條件270 MPa

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

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

1.3.3 自然地震的觸發(fā)

(1)地震能與震源能的關(guān)系

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

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為地震能，即地震發(fā)生時每秒釋放的震源能，J/s;t為時間變量，s;Px,y,z(t)為地震發(fā)生前時間t的區(qū)域主應(yīng)力，MPa;Px,y,z(t+Δt)為地震發(fā)生時的區(qū)域主應(yīng)力，MPa;ΔPx,y,z為經(jīng)歷Δt時刻的區(qū)域主應(yīng)力增量，MPa/s;D為常量；R為全體實數(shù)域。

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

(2)區(qū)域主應(yīng)力微小變化觸發(fā)地震的可能性分析

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

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

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

2 震例分析

2.1 區(qū)域構(gòu)造背景及有限差分模型

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

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

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

(7)

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

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

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

圖8 西南地區(qū)構(gòu)造簡圖[42]Fig.8 Geological structures of southwest China[42]

表2 模型物理力學(xué)參數(shù)Table 2 Physical and mechanical parameters of the model

2.2 邊界條件與初始條件

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

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

2.3 數(shù)值模擬結(jié)果分析

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

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

2.3.2不同應(yīng)力狀態(tài)下震源能的分布特征分析

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

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

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

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

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

2.3.3不同圍巖塑性區(qū)狀態(tài)下地震的觸發(fā)

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

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

4 結(jié) 論

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

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

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