張子平 袁 青 熊齊歡

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

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

引言

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

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

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

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

1 計(jì)算原理

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

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

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

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

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

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

假設(shè)層面一共有n+1個(gè)層面（圖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對(duì)應(yīng)為相應(yīng)層面的波速，與其對(duì)應(yīng)的入射角有關(guān)，，C即為射線參數(shù)。將式（3），（4）對(duì)時(shí)間t求導(dǎo)后分別代入式（2），在考慮透射系數(shù)的情況下可以得到P波、S波在通過地層分界層面（兩層巖石組成不同、物理性質(zhì)不同）后對(duì)目的地層的應(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波對(duì)應(yīng)的振幅。ω為波形函數(shù)。由于地震波速多為不規(guī)律的，因而采取某一時(shí)刻的數(shù)據(jù)進(jìn)行計(jì)算，t為對(duì)應(yīng)的時(shí)間。

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

2.1 模型構(gòu)建

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

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

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

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

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

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

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

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

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

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

圖5 加速度時(shí)程信息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)力伴隨時(shí)間完成由0瞬間增加至壓應(yīng)力極值再逐漸遞減的轉(zhuǎn)變，也與地震導(dǎo)致巖石破裂的時(shí)間范圍接近，后200 s的數(shù)據(jù)由于應(yīng)力積累等原因呈現(xiàn)周期性反復(fù)。而由于地震波中橫波應(yīng)力效應(yīng)一般大于縱波，也在圖中體現(xiàn)為P波的峰值大于S波。

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

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

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

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

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

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

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

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

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

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

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

在整個(gè)動(dòng)力過程中每個(gè)監(jiān)測(cè)點(diǎn)的Z向位移變化趨勢(shì)基本相同，初始都為0，隨時(shí)間不斷呈反向遞增趨勢(shì)，即整個(gè)模型體在縱向上不斷發(fā)生負(fù)向位移。整體上，上部地層的Z向位移比下部地層要少，因?yàn)榈卣鸬募铀俣仁侵苯邮┘釉谀Ｐ腕w的底面上。而監(jiān)測(cè)點(diǎn)9至監(jiān)測(cè)點(diǎn)10、監(jiān)測(cè)點(diǎn)11至監(jiān)測(cè)點(diǎn)12兩個(gè)層面間的Z向位移量相差無幾，此處的縱向位移不太明顯。

而對(duì)于中間斜面監(jiān)測(cè)點(diǎn)的X向位移來說，點(diǎn)7：初始位移為0，108 000 step（step為內(nèi)置時(shí)間步數(shù)，代表程序運(yùn)行時(shí)間，基本400 000 step對(duì)應(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)測(cè)點(diǎn)的X向位移均有向負(fù)向遞增的趨勢(shì)。對(duì)于上盤3個(gè)點(diǎn)來說，深度越大，變化幅度越小，峰值越小。而對(duì)于下盤的3個(gè)點(diǎn)來說，深度越大，峰值越小。在斜面上的監(jiān)測(cè)點(diǎn)呈現(xiàn)出與其Z向位移相當(dāng)?shù)奶攸c(diǎn)，從整體上來說，深度越大的層面，其發(fā)生的X向位移量越小。且在110000 step這個(gè)應(yīng)力激發(fā)時(shí)刻7、8、9、10，4個(gè)監(jiān)測(cè)點(diǎn)出現(xiàn)了正向位移。之后4個(gè)點(diǎn)又向負(fù)向位移，最后穩(wěn)定為負(fù)向位移。綜合3個(gè)圖，發(fā)現(xiàn)每個(gè)監(jiān)測(cè)點(diǎn)的X向位移均有所不同，即每一層位都存在X向的位移偏差。

3.2 模型體位移變化

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

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

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

4 討論與結(jié)論

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

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

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

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