基于蝶形破壞理論的地震能量來(lái)源

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

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

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

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

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

1 地震能量來(lái)源的理論分析

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

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

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

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

圖2 軟弱異性體周?chē)芰坑?jì)算分析模型Fig.2 Computational,and analytical,model for energies around a soft anisotropic body

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

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

(1)

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

(3)

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

(4)

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

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

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

(5)

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

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

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

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

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

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

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

1.3.1軟弱異性體周?chē)N(yùn)含能量的變化規(guī)律

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

圖4 軟弱異性體蘊(yùn)含震源能、里氏震級(jí)與水平主應(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)變化與對(duì)應(yīng)震源能密度分布特征。圖6中A,B,C分別對(duì)應(yīng)圖4劃分的階段(I，II，III)中任意取得的3個(gè)特征點(diǎn)。分析可得軟弱異性體圍巖形成的蝶形破壞區(qū)周?chē)辛舜罅康膹椥阅?，是地震發(fā)生時(shí)能量的主要來(lái)源。以擠壓應(yīng)力區(qū)為例，原巖應(yīng)力Px=270 MPa，塑性區(qū)的形狀為圓形。地質(zhì)構(gòu)造劇烈變化，水平應(yīng)力開(kāi)始增加，當(dāng)Px=810 MPa時(shí)(Px/Pz=3)，軟弱異性體周?chē)钠茐漠a(chǎn)生質(zhì)的轉(zhuǎn)變，呈蝶形變化特征。也就是說(shuō)，滿(mǎn)足條件270 MPa

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

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

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

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

依據(jù)前文定義:震源能是指某種破壞狀態(tài)下軟弱異性體周?chē)鷰r體內(nèi)的集中能量，是某一時(shí)刻的狀態(tài)量值，與時(shí)間無(wú)關(guān)。而真實(shí)地震發(fā)生時(shí)記錄到的以地震波形成釋放的能量是以時(shí)間為周期加以衡量，例如里氏震級(jí)(ML)、面波震級(jí)(Ms)及體波震級(jí)(Mb，MB分別為用1 s和5 s左右的地震體波振幅來(lái)量度地震的大小)中能量的計(jì)算都是以波的周期長(zhǎng)度作為時(shí)間單位。為了方便起見(jiàn)，本文定義單位時(shí)間釋放的震源能為地震能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ā)生時(shí)每秒釋放的震源能，J/s;t為時(shí)間變量，s;Px,y,z(t)為地震發(fā)生前時(shí)間t的區(qū)域主應(yīng)力，MPa;Px,y,z(t+Δt)為地震發(fā)生時(shí)的區(qū)域主應(yīng)力，MPa;ΔPx,y,z為經(jīng)歷Δt時(shí)刻的區(qū)域主應(yīng)力增量，MPa/s;D為常量；R為全體實(shí)數(shù)域。

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

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

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

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

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

2 震例分析

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

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

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

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

(7)

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

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

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

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

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

2.2 邊界條件與初始條件

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

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

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

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

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

圖11 特征點(diǎn)集中能量分布對(duì)比Fig.11 Comparison of distributions of energies accumulated at feature points

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

Pregnant Period:蝶形塑性破壞還未形成時(shí)期;Growth Period:蝶形塑性破壞形成初期;Upheaval Period:蝶形塑性破壞擴(kuò)展末期。

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

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

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

4 結(jié) 論

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

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

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