地震波誘發(fā)的深部地層動應(yīng)力機(jī)制及斷裂位移研究?以汶川地震對龍門山前陸盆地的作用為例*

張子平 袁 青 熊齊歡

1)中交第二航務(wù)工程局有限公司，武漢 430040

2)長大橋梁建設(shè)施工技術(shù)交通行業(yè)重點(diǎn)實(shí)驗室，武漢 430040

引言

2008年5月12日四川省汶川縣發(fā)生了8.0級大地震。地震造成了巨大的經(jīng)濟(jì)損失和民眾傷亡。伴隨著地震的發(fā)生，周邊地區(qū)地層產(chǎn)生了眾多的形變構(gòu)造，給后續(xù)的工程重建和人民財產(chǎn)安全帶來了許多隱患。

通過對震后遙感記錄和實(shí)地考察所得資料的調(diào)查研究，發(fā)現(xiàn)汶川地震后誘發(fā)了眾多的滑坡和褶皺斷層。曹俊興等[1-2]發(fā)現(xiàn)在距離汶川地震破裂帶30—40 km的龍門山前陸盆地地腹構(gòu)造處曾存在斷裂和褶皺變形。這些研究表明在汶川地震的影響下，近場至中場區(qū)域的地腹巖層出現(xiàn)了斷裂與褶皺變形。斷層作為一種分布極為廣泛的構(gòu)造形態(tài)，在許多構(gòu)造環(huán)境中均有發(fā)育。它不僅是一個與構(gòu)造地質(zhì)學(xué)、地震學(xué)、地球動力學(xué)有關(guān)的重要地質(zhì)現(xiàn)象，也同礦產(chǎn)資源的形成與分布、大型工程的基礎(chǔ)穩(wěn)定性、地震危險區(qū)劃分等有十分密切的聯(lián)系。

近年來，關(guān)于地震與斷層間的應(yīng)力應(yīng)變研究已有發(fā)展，早在20世紀(jì)中后期，地震學(xué)家們就用庫侖破裂應(yīng)力來描述地震引起的地層應(yīng)力應(yīng)變效果[3-6]。1951年，Anderson[7]根據(jù)斷層形成時3個主應(yīng)力的作用方向與相對大小，提出了斷層作用理論。Papadimitriou等[8]利用中國四川西部地震活動區(qū)域內(nèi)地震數(shù)據(jù)研究了1893年以來庫侖破裂函數(shù)的變化以及近110年該地區(qū)應(yīng)力場的演變，計算出地震引發(fā)的彈性半空間模型中的靜態(tài)位錯。戴俊[9]基于爆炸應(yīng)力波隨距離衰減原理研究沖擊波如何轉(zhuǎn)變成地震波。Aiken等[10]引入流體的影響來研究大范圍地震的動態(tài)庫侖應(yīng)力作用。Okada[11]利用半無限空間中任意傾角的傾向滑動剪切斷層錯動位移量計算出產(chǎn)生的應(yīng)力場與位移場。Wang等[12]利用格林函數(shù)來計算地震區(qū)域內(nèi)的水平分層粘彈性介質(zhì)模型。前人對地震作用下巖石構(gòu)造的破壞分析歷經(jīng)多年的研究，常采用的研究方法有Newmark分析法、動力有限元時程分析法[13]、動力有限元強(qiáng)度折減法[14]等，這些方法或理論分析，或構(gòu)建模型，均綜合考慮了地震影響下巖土體的破壞性問題。

這些研究方法雖然成果豐富，但少有對于深部地層破壞的研究，且沒有過多考慮巖層間波的折射衰減效應(yīng)?；诖?，本文通過波動力學(xué)和爆破力學(xué)衍生出的力學(xué)方程和地震波的層間折射效應(yīng)，提出巖層在地震波作用下的受力機(jī)制。依托汶川地震對龍門山前陸盆地深部地層的具體破壞效應(yīng)，根據(jù)實(shí)際地層的物理力學(xué)參數(shù)，對深部地層進(jìn)行應(yīng)力平衡分析。建立地層與斷層面模型，設(shè)置應(yīng)力約束和阻尼邊界，據(jù)Coulomb準(zhǔn)則和摩爾應(yīng)力圓判斷是否發(fā)生剪切破壞與位移情況。根據(jù)最后的位移量結(jié)果與監(jiān)測主應(yīng)力情況，得出在汶川地震波的沖擊作用下，龍門山前陸盆地深部地層的斷層破壞形變情況。同時證實(shí)此地震波應(yīng)力計算機(jī)制的可行性以及大地震對斷層構(gòu)造的誘發(fā)影響。

1 計算原理

1.1 一維情況下地震波應(yīng)力效應(yīng)

天然地震傳播的地震波主要包含縱波和橫波，而由此衍生的應(yīng)力也是由橫波和縱波產(chǎn)生的水平應(yīng)力和垂直應(yīng)力。其應(yīng)力效應(yīng)在一維情況下可由波動力學(xué)公式求解[15]：

式中，σ為地震波產(chǎn)生的動應(yīng)力，縱波影響垂直應(yīng)力，橫波影響水平應(yīng)力；ρ為巖體密度；c為巖體中的波速；μ為巖體中質(zhì)點(diǎn)的振動速度。爆破地震學(xué)[16]中由動量守恒、質(zhì)量守恒和能量守恒三大定律導(dǎo)出的巖石中沖擊波參數(shù)計算公式：

式中，g為重力加速度；U為波陣面上質(zhì)點(diǎn)的運(yùn)動速度；Ub為沖擊波傳播速度（地震波波速）；γ0為比重（震前巖石密度）?？梢钥闯?，兩者的計算公式具有高度的一致性，因為天然地震波是彈性波[17]，而爆炸發(fā)出的沖擊波和應(yīng)力波向外傳播后也會衰減為彈性波。

1.2 地震波在分界面的折射與應(yīng)力衰減

巖層密度和巖體中的波速均可由實(shí)驗或當(dāng)?shù)氐刭|(zhì)情況得到，故巖體中質(zhì)點(diǎn)的振動速度就是所需要測算的關(guān)鍵因子。對均質(zhì)地層分層后進(jìn)行橫波縱波透射后的位移量計算：

假設(shè)層面一共有n+1個層面（圖1），下部介質(zhì)為第1層，上部介質(zhì)為第n+1層，中間n?1層。P波和S波的入射波為Pr、Sr；反射波為Pf、Sf；透射波為Pt、St。其中的入射角為。其中上部介質(zhì)，第n+1層的介質(zhì)的總位移量（包括入射波與反射波）為[18]：

圖1 地震波在厚地層中的反射與透射Fig. 1 Reflection and transmission of seismic waves in thick formations

P波：

S波：

式中，A，B是P波和S波的振幅；j為常數(shù)，ω是波形數(shù)據(jù)。且，這里的αn+1和βn+1對應(yīng)為相應(yīng)層面的波速，與其對應(yīng)的入射角有關(guān)，，C即為射線參數(shù)。將式（3），（4）對時間t求導(dǎo)后分別代入式（2），在考慮透射系數(shù)的情況下可以得到P波、S波在通過地層分界層面（兩層巖石組成不同、物理性質(zhì)不同）后對目的地層的應(yīng)力影響：

式中，j為常數(shù)。ρ1為所要求解層面的密度，cP1為其P波波速，cS1為其S波波速；ρ3則為入射進(jìn)此層面的初始層面的密度，cP3為其P波波速，cS3為其S波波速。iP，iS為P波和S波的入射角。x是目的層位距離震源的水平距離，z是目的層位距離震源的豎直距離。A、B分別為P波、S波對應(yīng)的振幅。ω為波形函數(shù)。由于地震波速多為不規(guī)律的，因而采取某一時刻的數(shù)據(jù)進(jìn)行計算，t為對應(yīng)的時間。

2 模型構(gòu)建與應(yīng)力計算監(jiān)測

2.1 模型構(gòu)建

龍門山前陸盆地位于龍門山?jīng)_斷帶與龍泉山前陸隆起之間，北起安縣秀水，南抵名山、彭山一線，西部已卷入龍門山?jīng)_斷帶[19]。龍門山前陸盆地自古以來就充填了大量的海陸相沉積物[20]，以時間可分為三疊系、侏羅系、白堊系、古近系、新近系和第四系6個層位。根據(jù)本文的研究內(nèi)容，基于前人勘探資料，對巖層的物理力學(xué)參數(shù)進(jìn)行平均化計算。共計建立8個層位，表1為其物理力學(xué)參數(shù)[21-23]。

表1 各地層的物理力學(xué)參數(shù)與圍巖應(yīng)力Table 1 Physical mechanical parameters and surrounding rock stress of each layer

基于研究地震波引發(fā)動力對深部地層造成的破壞效應(yīng)，將地層看作一個整體構(gòu)造，建立15 km×5 km×10 km的立方體模型（圖2），內(nèi)部設(shè)置斜面模型觀測地震波動力效果下斷層面上的應(yīng)力應(yīng)變情況。由于采集數(shù)據(jù)的地震臺站所處區(qū)域的斷層帶多為30°傾角的逆斷層，故設(shè)計此結(jié)構(gòu)對其進(jìn)行分析，通過多次研究發(fā)現(xiàn)，20°—45°傾角范圍內(nèi)整個地層結(jié)構(gòu)體的受力和位移情況相差不大，其他傾角情況下，傾角越小，越容易產(chǎn)生錯斷位移情況。

圖2 主體地層模型Fig. 2 The model of main stratigraphic

基于莫爾—庫侖失效準(zhǔn)則對于斷層失穩(wěn)的判斷分析，對斷層前后層面周邊的地層情況進(jìn)行數(shù)值模擬，通過監(jiān)測巖層的應(yīng)力與位移形變信息來具體分析。設(shè)置在斜面上和周圍區(qū)域共13個監(jiān)測點(diǎn)（圖3）。

圖3 監(jiān)測點(diǎn)在XOZ平面上的投影Fig. 3 Projection of the monitoring point on the XOZ plane

2.2 應(yīng)力計算與監(jiān)測

震源距離所取數(shù)據(jù)的地震臺站水平距離36 km，如圖4所示分為9等份，每4 km取一個質(zhì)點(diǎn)。z對應(yīng)每層的厚度（圖4）。采集地震臺站（江油市含增鎮(zhèn)編號51JYH）記錄的加速度時程信息作為初始條件，模擬地震波影響下的動力阻尼。

圖4 地下分層情況與質(zhì)點(diǎn)分布Fig. 4 Underground stratification and particle distribution

對巖石地基，在底部施加加速度轉(zhuǎn)化的應(yīng)力荷載。采取以0.02 s為間隔的12500個樣本時間數(shù)據(jù)（圖5）。根據(jù)式（5）、（6）以及表1中的巖石物理參數(shù)計算對應(yīng)的各地層質(zhì)點(diǎn)P波、S波在通過地層分界層面（兩層巖石組成不同、物理性質(zhì)不同）后對目的地層的應(yīng)力影響（以第7層質(zhì)點(diǎn)2為例）（圖6）。

圖5 加速度時程信息Fig. 5 Acceleration time history information

圖6 第7層質(zhì)點(diǎn)2處的P波和S波應(yīng)力影響Fig. 6 The influence of P-wave and S-wave stress at the particle #2 of the 7th layer

可以看出，P波和S波的應(yīng)力在（?100，100） MPa區(qū)間內(nèi)。S波在初始的10 s內(nèi)峰值達(dá)到?36 MPa左右，而P波在初始的10 s內(nèi)峰值達(dá)到?50 MPa。其中正數(shù)值代表拉應(yīng)力，負(fù)數(shù)值代表壓應(yīng)力。前50 s內(nèi)的應(yīng)力數(shù)據(jù)具有代表性，應(yīng)力伴隨時間完成由0瞬間增加至壓應(yīng)力極值再逐漸遞減的轉(zhuǎn)變，也與地震導(dǎo)致巖石破裂的時間范圍接近，后200 s的數(shù)據(jù)由于應(yīng)力積累等原因呈現(xiàn)周期性反復(fù)。而由于地震波中橫波應(yīng)力效應(yīng)一般大于縱波，也在圖中體現(xiàn)為P波的峰值大于S波。

通過將式（5）、（6）計算好的各地層地震波P波和S波動應(yīng)力影響施加于建立好的模型體上，發(fā)現(xiàn)隨著地震波的震蕩變化，地震應(yīng)力和不平衡力變化（圖7）也呈現(xiàn)上下波動的趨勢。

圖7 動力過程中的最大不平衡力Fig. 7 Maximum unbalanced force during dynamic process

隨著地震波加速度數(shù)據(jù)的正負(fù)轉(zhuǎn)化，與其成正相關(guān)關(guān)系的地震應(yīng)力也隨之發(fā)生變化，導(dǎo)致最大不平衡力一直在0與最大值之間波動，且由于大尺度模型的自由邊界問題，未出現(xiàn)平衡收斂。從這也可以推測出，地層結(jié)構(gòu)已出現(xiàn)了破壞。對靜力平衡下及動力施加過程中上盤3個點(diǎn)的主應(yīng)力情況（圖8）進(jìn)行統(tǒng)計。

圖8 動力過程中的主應(yīng)力變化Fig. 8 The change of principal stress during dynamic process

根據(jù)表2的數(shù)據(jù)繪制摩爾應(yīng)力圓和大小主應(yīng)力的變化曲線，藍(lán)色摩爾圓由靜力下的最大最小主應(yīng)力峰值所繪制，選取此情況為極限平衡條件的摩爾圓。從動力施加下的主應(yīng)力情況來看，110000動態(tài)時間步后主應(yīng)力趨于穩(wěn)定（表2）。各點(diǎn)上最大主應(yīng)力的變化情況均大于最小主應(yīng)力的變化，判斷紅色主應(yīng)力圓會相交于同監(jiān)測點(diǎn)藍(lán)色穩(wěn)定情況的摩爾圓（圖9），即此位置的層面已發(fā)現(xiàn)了斷裂破壞。初步判斷在汶川地震影響下龍門山前陸盆地深部地層出現(xiàn)了斷裂構(gòu)造。

圖9 監(jiān)測點(diǎn)的主應(yīng)力變化Fig. 9 The change of principal stress at monitoring points

表2 監(jiān)測點(diǎn)的主應(yīng)力變化Table 2 The change of principal stress at monitoring points

3 監(jiān)測點(diǎn)與模型體的位移分析

3.1 監(jiān)測點(diǎn)位移變化

主應(yīng)力可初步判斷是否發(fā)生斷層錯位，監(jiān)測各點(diǎn)各層位的位移變化可以驗證每層的錯斷與否并統(tǒng)計其位移變化情況。圖10為13個監(jiān)測點(diǎn)的Z向位移和中間斜面監(jiān)測點(diǎn)的X向位移記錄。

圖10 動力過程中的監(jiān)測點(diǎn)位移Fig. 10 Displacement of monitoring point during dynamic process

在整個動力過程中每個監(jiān)測點(diǎn)的Z向位移變化趨勢基本相同，初始都為0，隨時間不斷呈反向遞增趨勢，即整個模型體在縱向上不斷發(fā)生負(fù)向位移。整體上，上部地層的Z向位移比下部地層要少，因為地震的加速度是直接施加在模型體的底面上。而監(jiān)測點(diǎn)9至監(jiān)測點(diǎn)10、監(jiān)測點(diǎn)11至監(jiān)測點(diǎn)12兩個層面間的Z向位移量相差無幾，此處的縱向位移不太明顯。

而對于中間斜面監(jiān)測點(diǎn)的X向位移來說，點(diǎn)7：初始位移為0，108 000 step（step為內(nèi)置時間步數(shù)，代表程序運(yùn)行時間，基本400 000 step對應(yīng)250 s）之后向X正向遞增到7 m，再負(fù)向位移，158 000 step回到初始位置，繼續(xù)負(fù)向位移至?28 m左右。點(diǎn)8：初始位移為0，104 000 step之后向X正向遞增到7 m，再負(fù)向位移，140 000 step回到初始位置，繼續(xù)負(fù)向位移至?47 m左右。點(diǎn)9：初始位移為0，100 000 step之后向X正向遞增到6 m，再負(fù)向位移，124 000 step回到初始位置，繼續(xù)負(fù)向位移至?70 m左右。點(diǎn)10：初始位移為0，100 000 step之后向X正向遞增到4 m，再負(fù)向位移，117 000 step回到初始位置，繼續(xù)負(fù)向位移至?84 m左右。點(diǎn)11：初始位移為0，100 000 step之后向X正向遞增到3 m，再負(fù)向位移，105 000 step回到初始位置，繼續(xù)負(fù)向位移至?118 m左右。點(diǎn)12：初始位移為0，基本一直負(fù)向遞增位移，最后達(dá)到?136 m左右。點(diǎn)13：初始位移為0，基本一直負(fù)向遞增位移，最后達(dá)到?168 m左右。

總體上說，各監(jiān)測點(diǎn)的X向位移均有向負(fù)向遞增的趨勢。對于上盤3個點(diǎn)來說，深度越大，變化幅度越小，峰值越小。而對于下盤的3個點(diǎn)來說，深度越大，峰值越小。在斜面上的監(jiān)測點(diǎn)呈現(xiàn)出與其Z向位移相當(dāng)?shù)奶攸c(diǎn)，從整體上來說，深度越大的層面，其發(fā)生的X向位移量越小。且在110000 step這個應(yīng)力激發(fā)時刻7、8、9、10，4個監(jiān)測點(diǎn)出現(xiàn)了正向位移。之后4個點(diǎn)又向負(fù)向位移，最后穩(wěn)定為負(fù)向位移。綜合3個圖，發(fā)現(xiàn)每個監(jiān)測點(diǎn)的X向位移均有所不同，即每一層位都存在X向的位移偏差。

3.2 模型體位移變化

對每100 step的時間順序依次截取X-Z面的剪切應(yīng)變變化進(jìn)行分析（圖11）。

圖11 X-Z面監(jiān)測的剪切應(yīng)變變化Fig. 11 Shear strain change monitored by X-Z plane

可以看出，因為施加的地震數(shù)據(jù)是置于底面的，整個模型的錯位是從底層開始的，紅色區(qū)域的層位首先發(fā)生斷裂，然后往上，黃色，綠色，藍(lán)色層位隨著時間的推移相繼發(fā)生斷裂，最后形成斷層構(gòu)造。在整個過程中，下部巖層的剪切形變量均大于上部巖層。2 500 step后全部層位已基本完成錯斷，與1、2、3監(jiān)測點(diǎn)的主應(yīng)力信息相符，表明這些位置已發(fā)生斷層偏移現(xiàn)象。地震加速度釋放的350000 step時刻后，全部層位基本發(fā)生錯位。而與剪切應(yīng)變隨時間變化的特征相類似，作為能代表整體形變特征的體積應(yīng)變變化也是從底面開始，隨時間推移，上部巖層慢慢完成層位的錯斷。伴隨著分界斜面，也就是各圖中的斜線，上面的巖層相對斜線以下巖層位移向上偏移。將斜線看作假想斷層面的話，斷層上盤向上位移，下盤向下位移，即發(fā)生了逆斷層。而整個模型體的X向位移本來應(yīng)該是連續(xù)而平緩的，由于上盤巖層向上位移，下盤巖層向下，導(dǎo)致整體X向位移出現(xiàn)不連續(xù)的斷裂。

4 討論與結(jié)論

本文提出了一種考慮地震波分層折射影響的地震應(yīng)力計算機(jī)制，并對龍門山前陸地腹構(gòu)造受汶川地震影響的應(yīng)力與形變進(jìn)行數(shù)值分析，探究大地震對深部巖層斷層構(gòu)造的誘發(fā)機(jī)制。得到以下結(jié)論：

（1）利用分層介質(zhì)傳播衰減規(guī)律計算出的汶川地震波應(yīng)力影響能對應(yīng)各時間點(diǎn)地震波計算出同時期的應(yīng)力影響，并考慮了波在地層分界面之間的折射（地震波能量衰減或增強(qiáng)）效應(yīng)，使結(jié)果具有時間上的連續(xù)性與動態(tài)效應(yīng)，模擬出巖層破壞的變化趨勢，反映地震波對深部地層的動應(yīng)力影響及位移—時間變化。

（2）通過記錄上盤監(jiān)測點(diǎn)的最大、最小主應(yīng)力變化情況，初步判斷龍門山前陸盆地地腹巖石出現(xiàn)了剪切破壞情況，得出汶川地震確實(shí)引發(fā)了深部地層的斷裂構(gòu)造。跟蹤這些點(diǎn)的X與Z向位移，發(fā)現(xiàn)各層位位移量出現(xiàn)差距，也佐證了這一點(diǎn)。在本次的模擬中，各層的Z向位移差別較小，而X向位移偏差較大，且隨著深度的變化具有較為明顯的差異性，深度越大，位移量越小。

（3）通過對模型X—Z面的剪切形變與體積形變進(jìn)行分析，可以很明顯的發(fā)現(xiàn)各層位都有錯斷結(jié)構(gòu)，且根據(jù)體積應(yīng)變的變化趨勢和最后結(jié)果，可以判斷斷層為逆斷層。對整個模型體在地震影響下的X向位移結(jié)果進(jìn)行分析也可佐證此點(diǎn)。