董載天 王軍民

摘 ?????要： 由于隧道超前探測(cè)的環(huán)境與常規(guī)石油地震勘探的差異，不宜直接將過(guò)去的模擬方法拿到隧道中進(jìn)行應(yīng)用。結(jié)合層狀地層表現(xiàn)出橫向各向同性的特性，推導(dǎo)出在隧道空間適用的二維VTI介質(zhì)一階qP波方程，分析了其穩(wěn)定條件。建立斷層和軟弱夾層的幾種數(shù)值模型，利用PML邊界解決模型邊界反射，采用交錯(cuò)網(wǎng)格進(jìn)行有限差分正演模擬，得到其共炮點(diǎn)道集。結(jié)果表明，qP波有限差分法可以快速地模擬不良的地震響應(yīng)，規(guī)避偽S波和數(shù)值不穩(wěn)定問(wèn)題，為隧道超前探測(cè)資料的解釋提供依據(jù)，提高預(yù)測(cè)的準(zhǔn)確性。

關(guān) ?鍵 ?詞：隧道超前探測(cè);VTI介質(zhì);qP波;PML

中圖分類(lèi)號(hào)：TV221.2???????文獻(xiàn)標(biāo)識(shí)碼：?A ??????文章編號(hào)： 1671-0460（2020）03-0683-07

Forward Modeling of Tunnel Advanced

Detection qP?Wave Based on VTI Medium

DONG Zai-tian， WANG Jun-min

（Yangtze University， Hubei Wuhan?430100， China）

Abstract： ?Due to the difference between the environment of tunnel pre-detection and conventional petroleum seismic exploration， it is not appropriate to directly apply the past simulation method to the tunnel. Combined with the transverse isotropy?characteristic?of layered stratum， the first-order qP wave equation of the two-dimensional VTI medium suitable for tunnel space was?derived， and the stability conditions were?analyzed. Several numerical models of fault and weak interlayer were?established. The boundary reflection of the model was?solved by PML boundary， and the finite difference forward modeling was?carried out by staggered grid to obtain the common shot gather. The results show that the qP wave finite difference method can quickly simulate the bad seismic response， avoid the pseudo S wave and numerical instability problem， and provide a basis for the interpretation of tunnel advance detection data?to?improve the accuracy of prediction.

Key words： ?tunnel advance detection; ?TI media;??qP wave; ?PML

在我國(guó)經(jīng)濟(jì)高速發(fā)展的當(dāng)下，大量公路鐵路工程都有隧道施工。在隧道工程施工過(guò)程中，由于掌子面前方地質(zhì)狀況不詳，面對(duì)破碎帶、裂隙和溶洞的支護(hù)不當(dāng)或支護(hù)不及時(shí)，經(jīng)常出現(xiàn)坍塌、突水突泥、冒頂?shù)戎卮鬄?zāi)害，造成嚴(yán)重的人員傷亡和財(cái)產(chǎn)損失。顯然在隧道工程中，使用超前探測(cè)技術(shù)了解掌子面前方地質(zhì)狀況，對(duì)工程施工安全極其重要。超前探測(cè)方法從原理上可以劃分為：地震波法、電磁法（瞬變電磁、電阻率、激發(fā)極化、核磁共振）、巖體溫度法等，其中地震波法在遠(yuǎn)距離勘探中有較好的分辨率，成為目前隧道超前探測(cè)中的常用方法之一。在隧道地震波法探測(cè)中，常用的方法有VSP、HSP^[1]、TSP^[2]、TVSP^[3]、陸地聲納^[4]、HSP^[5]、TRT^[6，7]、TST^[8]、ISIS^[9]等方法。由于隧道施工現(xiàn)場(chǎng)的空間狹小且工程工期壓力大，所以隧道超前探測(cè)技術(shù)必須滿(mǎn)足狹小場(chǎng)地條件下進(jìn)行大距離、高精度、高效率探測(cè)要求^[10]，這就要求在對(duì)傳統(tǒng)方法技術(shù)（測(cè)線(xiàn)與探測(cè)對(duì)象平行）進(jìn)行改進(jìn)（測(cè)線(xiàn)與探測(cè)對(duì)象垂直）的同時(shí)，進(jìn)一步研究地震波在隧道復(fù)雜條件下的傳播機(jī)理特性^[11]。因此，本文在層狀地層條件下，構(gòu)建斷層、軟弱夾層等不良地質(zhì)災(zāi)害體模型，從一階速度?應(yīng)力qP波方程出發(fā)，采用交錯(cuò)網(wǎng)格高階有限差分進(jìn)行隧道復(fù)雜波場(chǎng)的正演模擬計(jì)算，以充分認(rèn)識(shí)不同條件下地震波響應(yīng)特征。

1 ?TI介質(zhì)的波動(dòng)方程

在我國(guó)，層狀地層分布廣泛，約占陸地面積的77.3%^[12]，在隧道施工過(guò)程中大多要面對(duì)這樣的地質(zhì)環(huán)境，而地下層狀在橫向各向同性地層往往表現(xiàn)出橫向各向同性的物理響應(yīng)特征。介于筆者使用波動(dòng)方程法進(jìn)行地震數(shù)值模擬，所以推導(dǎo)出橫向各向同性（TransverseIsotropy）介質(zhì)的波動(dòng)方程尤為重要。

1.1 ?波動(dòng)方程導(dǎo)出

波動(dòng)方程的建立涉及應(yīng)力、應(yīng)變和位移這三個(gè)物理量，需要使用彈性體的本構(gòu)方程、運(yùn)動(dòng)微分方程和柯西方程。由本構(gòu)方程，可知彈性體應(yīng)力與應(yīng)變關(guān)系表達(dá)式：

????????（1）

式中：?— 應(yīng)力張量;

e?— 應(yīng)變張量;

C?— 彈性常數(shù)。

彈性常數(shù)共有81個(gè)，由柯西方程可知，應(yīng)變張量有對(duì)稱(chēng)性：

???????????????（2）

其中：u、v —位移分量;

x，?y?—笛卡爾坐標(biāo)。

而應(yīng)力張量也有此種對(duì)稱(chēng)性，即。因此彈性常數(shù)減少為36個(gè)，同時(shí)彈性常數(shù)張量相對(duì)于主對(duì)角線(xiàn)對(duì)稱(chēng)，有。使用單下標(biāo)替代雙下標(biāo)11→1，22→2，33→3，23（32）→4，13（31）→5，12（21）→6，則可的表達(dá)式：

??????（3）

上式可簡(jiǎn)記為：

?????????????（4）

從彈性體的角度，描述完應(yīng)變與應(yīng)力關(guān)系、應(yīng)變與位移關(guān)系后，通過(guò)牛頓第二定律可建立運(yùn)動(dòng)平衡微分方程：

??????????????（5）

式中：—?彈性體密度，kg/m;

U?— 位移張量;

— 應(yīng)力張量;

t — 時(shí)間，s;

F —體力張量;

L —偏導(dǎo)算子。

偏導(dǎo)算子L為：

???????（6）

利用Cauchy方程可將（5）式表達(dá)為：

?????????（7）

1.2 ?TI介質(zhì)中qP波波動(dòng)方程

地表巖土物理性質(zhì)極其復(fù)雜不便于數(shù)值模擬研究，地球物理學(xué)家為了能借助波動(dòng)物理研究其波傳播特性，根據(jù)各項(xiàng)異性彈性體介質(zhì)的基本對(duì)稱(chēng)性進(jìn)行了分類(lèi)。其中三方、六方各項(xiàng)異性介質(zhì)常叫做橫向各向同性（Transverse Isotropy，簡(jiǎn)稱(chēng)TI）介質(zhì)，是目前實(shí)際應(yīng)用研究中使用最廣泛的模型。TI介質(zhì)又根據(jù)對(duì)稱(chēng)軸的水平和垂直細(xì)分為，VTI介質(zhì)和HTI介質(zhì)（圖1）。

兩介質(zhì)都需要5個(gè)獨(dú)立的彈性參數(shù)來(lái)描述，其彈性矩陣分別為：

分別將兩式帶入公式（7），可得VTI介質(zhì)三維波動(dòng)方程：

由于VTI介質(zhì)與地下實(shí)際介質(zhì)的響應(yīng)特質(zhì)較為一致且研究最為廣泛，所以筆者主要使用VTI介質(zhì)的彈性波波動(dòng)方程，因此此處不再給出HTI介質(zhì)的波動(dòng)方程。

在利用此方程進(jìn)行TI介質(zhì)的彈性波正演模擬中，波動(dòng)方程的包含三個(gè)分量，占用計(jì)算機(jī)大量?jī)?nèi)存，同時(shí)由于縱波速度快于橫波，而網(wǎng)格的空間步長(zhǎng)和采樣時(shí)間步長(zhǎng)主要受最波的小速度影響，為穩(wěn)定計(jì)算，只能減小空間和采樣的間隔，使得計(jì)算的時(shí)間點(diǎn)數(shù)與空間點(diǎn)數(shù)大幅提高，進(jìn)一步減慢了計(jì)算速度，不利于計(jì)算機(jī)模擬。且目前實(shí)際勘探仍以縱波為主，多分量數(shù)據(jù)中進(jìn)行橫波分離仍然是個(gè)難點(diǎn)問(wèn)題。

因此筆者借用Alkhalifah提出的聲學(xué)假設(shè)條件，直接將橫波速度設(shè)為零，從而在復(fù)雜的 VTI 介質(zhì)彈性波方程中分裂出qP波方程^[14]。該方程不產(chǎn)生橫波，雖然現(xiàn)實(shí)中上不可實(shí)現(xiàn)，但在運(yùn)動(dòng)學(xué)上能對(duì)彈性波方程進(jìn)行有效地近似，便于研究TI介質(zhì)中qP波的傳播規(guī)律。

根據(jù)Alkhalifah的假設(shè)導(dǎo)出的三維 VTI 介質(zhì)中的qP波方程為：

???（11）

式中：F?— 為波場(chǎng)函數(shù);

t?— 為時(shí)間，s;

x、y、z?— 為三維空間坐標(biāo);

v 、v_v?— qP波的動(dòng)校正速度和垂向速度，m/s;

η?— 各向異性參數(shù)。

已知：e、d 為 Thomsen 參數(shù)^[15]。

由于隧道超前探測(cè)觀測(cè)方向?yàn)檎谱用嬲胺?，而常?guī)地震勘探觀測(cè)方向?yàn)榈孛娲瓜颍ㄈ鐖D2），所以筆者后面推導(dǎo)VTI介質(zhì)中的波動(dòng)方程時(shí)只保留xy分量。

因此，由式（11）可得二維VTI介質(zhì)中qP波方程：

直接求解上式需要消耗大量?jī)?nèi)存和時(shí)間，所以利用，可將上式轉(zhuǎn)化為：

其中為應(yīng)力。借鑒Hestholm推導(dǎo)二維VTI介質(zhì)中qP波方程的一階應(yīng)力-速度方程形式的方法，可得：

式中：— 彈性體密度，kg/m;

v_x、v?_y— x、y方向速度，m/s;

?— 計(jì)算過(guò)程中的中間變量。

2 ?交錯(cuò)網(wǎng)格

在實(shí)際計(jì)算過(guò)程中必須對(duì)介質(zhì)模型進(jìn)行離散，將介質(zhì)模型網(wǎng)格化。實(shí)際應(yīng)用過(guò)程中，交錯(cuò)網(wǎng)格相比規(guī)則網(wǎng)格有更好的差分算子局部特性，更高的精度和更好的收斂性。所以筆者采用交錯(cuò)網(wǎng)格（如圖3）將放置在網(wǎng)格節(jié)點(diǎn)。

利用Taylor級(jí)數(shù)展開(kāi)法，用網(wǎng)格節(jié)點(diǎn)上數(shù)值的差商近似代替前式中的導(dǎo)數(shù)。推導(dǎo)出了時(shí)間域 2 階、空間域2N階精度的交錯(cuò)網(wǎng)格離散差分格式。

式中：?— 空間間隔，m;

— 時(shí)間間隔，s;

i、j、k— 空間和時(shí)間節(jié)點(diǎn);

U、V— 速度分量v_x和v_y的離散值，m/s;

P、Q、R— 應(yīng)力和的離散值。

3 ?有限差分模擬中的重要問(wèn)題

3.1 ?PML邊界條件

針對(duì)隧道超前探測(cè)的彈性波數(shù)值模擬而設(shè)計(jì)的數(shù)據(jù)模型，僅考慮深埋隧道時(shí)，隧道空間為自由邊界，計(jì)算區(qū)域外部邊界則都是人工截?cái)噙吔?。人工邊界?huì)使向外傳播的彈性波完全反射回去，從而使正演模擬計(jì)算結(jié)果與實(shí)際情況大不符。為了盡量降低人工截?cái)噙吔缢a(chǎn)生的影響，引入Berenger針對(duì)電磁波傳播提出了完美匹配層（Perfect Matched Layer，簡(jiǎn)稱(chēng) PML）吸收邊界條件，削弱人工截?cái)噙吔绲姆瓷洌纳普菽M效果。

根據(jù)PML的推導(dǎo)思路，可給出式（14）的吸收邊界方程：

式中：d（x）、d（y） — 阻尼因子;

R?— 反射系數(shù);

?— 完美匹配層吸收邊界的厚度，m;

h?— 完美匹配層吸收邊界的深度，m;

v?— 匹配層中的彈性波波速，m/s。

3.2??穩(wěn)定性

波動(dòng)方程的差分格式按照時(shí)間進(jìn)行遞推，前一時(shí)刻的波函數(shù)數(shù)值解的舍入誤差必然會(huì)對(duì)下一時(shí)刻的波場(chǎng)產(chǎn)生影響。這就需要對(duì)誤差傳播情況進(jìn)行估計(jì)，使其不會(huì)隨著時(shí)間推移而快速增長(zhǎng)，以致影響數(shù)值模擬的效果進(jìn)而破壞計(jì)算，使模擬無(wú)法繼續(xù)進(jìn)行。筆者此處借鑒Virieux在1986年提出的推導(dǎo)方法，得到二維VTI介質(zhì)中二階時(shí)間差分精度、2N階空間差分精度的穩(wěn)定條件：

??????（17）

當(dāng)中表示模型中最大速度的平方。當(dāng)模型橫縱空間步長(zhǎng)一致時(shí)，亦可將方程改寫(xiě)為：

?????????????（18）

3.3??震源加載

此部分還有一個(gè)重要問(wèn)題就是加載震源。本文為保證成像的分辨率，選擇使用零相位雷克子波作為震源。其時(shí)間域表達(dá)形式為：

????????（19）

式中：t₀?—初始時(shí)刻，ms;

f?—主頻， Hz。

4 ?數(shù)值模擬實(shí)例

建立隧道空間的模型（如圖4），正演模型的尺寸為300 m×100 m，模型四周使用PML吸收邊界。

隧道尺寸為60 m×8 m，空腔中充填空氣，波速340?m/s;?圍巖波速3 600 m/s。采樣間隔為0.2 ms，采樣時(shí)長(zhǎng)為0.1 s。

觀測(cè)方式如所示，掌子面中心單震源，左右兩側(cè)墻壁各布置4個(gè)檢波器，間隔2 m，水平最小炮間距10 m。震源與檢波器皆置于墻內(nèi)，埋深2 m。

通過(guò)正演模擬可得到地震響應(yīng)記錄（圖6），為共炮點(diǎn)道集（CSP）?？梢钥闯鲋边_(dá)波能量很強(qiáng)，無(wú)邊界反射和偽S波干擾。根據(jù)費(fèi)馬原理計(jì)算的初至?xí)r間在4 ms左右，基本與記錄吻合。

4.1 ?不良地質(zhì)體模型及其地震響應(yīng)

隧道施工過(guò)程中會(huì)遇到很多危害施工安全的不良地質(zhì)體，常見(jiàn)的有：斷層破碎帶、軟弱夾層等。由于二維平面不存在傾角，所以后續(xù)模型中需要關(guān)注不良地質(zhì)體的走向。

4.1.1 ?斷層模型及響應(yīng)

斷層破碎帶是由于斷層上下盤(pán)搓動(dòng)引起斷面附近的巖體發(fā)生破碎，形成破碎帶。其穩(wěn)定性較差，容易引發(fā)隧道塌方;在地下水富集區(qū)域，破碎帶變成了一條良好水通道，極易引發(fā)突水突泥災(zāi)害。由于破碎帶與斷層的伴生特性且破碎帶本身理化性質(zhì)過(guò)于復(fù)雜，筆者此處將斷層破碎帶的模型設(shè)置為斷面兩側(cè)阻抗不同，破碎帶并不在模型當(dāng)中直接體現(xiàn)。數(shù)值模型構(gòu)建時(shí)，首先選擇走向與隧道方向垂直，距掌子面80和120 m的斷層，斷層右盤(pán)波速為4 000 m/s（圖7）。

對(duì)比兩模型的CSP道集（如圖8）可以看出，當(dāng)斷面與掌子面距離減小時(shí)，阻抗界面的反射波逐漸靠近強(qiáng)振幅的直達(dá)波，初至有可能被遮蔽，此時(shí)僅靠時(shí)域記錄是難以判斷前方斷面距離的;當(dāng)兩者間距拉大，反射波整體向后時(shí)移，初至明顯但振幅減弱，現(xiàn)場(chǎng)施工噪音較強(qiáng)時(shí)，會(huì)難以見(jiàn)到完整波形。

4.1.2??軟弱夾層及響應(yīng)

軟弱夾層成因復(fù)雜且尚無(wú)統(tǒng)一的分類(lèi)原則，此處只針對(duì)沉積型軟弱夾層進(jìn)行分析。此種軟弱夾層往往形成與河流沉積環(huán)境，層中粉砂巖、泥巖含量偏高，強(qiáng)度低質(zhì)地較軟弱^[17]。當(dāng)構(gòu)造作用使夾層被剪切、破碎時(shí)，無(wú)論區(qū)域是否富水，都會(huì)產(chǎn)生低速低密度的軟弱夾層。軟弱夾層的數(shù)值模型參考斷層模型，將夾層設(shè)置在掌子面前方80 m處且走向與隧道方向垂直，此時(shí)一次反射波不會(huì)被直達(dá)波遮蔽。分別構(gòu)建10、20 m厚的夾層，波速為800 m/s（如圖9）。

此時(shí)獲得CSP道集（如圖10）。由于此時(shí)夾層與圍巖阻抗差異過(guò)大，且檢波器內(nèi)置在墻壁中，導(dǎo)致阻抗界面反射波振幅較強(qiáng)。對(duì)于薄夾層來(lái)說(shuō)，子波在夾層中旅行時(shí)較短，夾層右側(cè)依舊為高阻抗差界面，所以10 m薄夾層的地震反射記錄上出現(xiàn)多次波。而層厚為20 m的夾層，子波旅行時(shí)被拉長(zhǎng)此時(shí)反射波沒(méi)有在記錄時(shí)間內(nèi)返回到檢波器，故道集上無(wú)明顯反射波。

將軟弱夾層與掌子面之間距離拉大至120 m（模型如圖11）。

從兩模型CSP道集（如圖12）可以看出，此時(shí)兩種模型的地震響應(yīng)頗為一致，顯然夾層右界面的反射波由于旅行時(shí)較長(zhǎng)沒(méi)有在記錄時(shí)長(zhǎng)內(nèi)返回到檢波器位置。

構(gòu)建一個(gè)右盤(pán)為低速巖體的斷層（如圖13），此時(shí)獲得其CSP道集（如圖14），其與掌子面距離為120 m的軟弱夾層的地震響應(yīng)基本一致。可見(jiàn)固定采樣時(shí)長(zhǎng)有時(shí)無(wú)法分辨的不同的不良地質(zhì)體，此時(shí)需要增加采樣時(shí)長(zhǎng)或者利用其他地質(zhì)勘探資料來(lái)確定不良地質(zhì)體形態(tài)。

5??結(jié) 論

隧道超前探測(cè)尚處于發(fā)展階段，由于其探測(cè)空間的復(fù)雜以及其與傳統(tǒng)地震勘探的差異，使得隧道超前探測(cè)在實(shí)際應(yīng)用過(guò)程中存在一些問(wèn)題。筆者從彈性體基本理論出發(fā)，逐步推導(dǎo)出適用于隧道超前探測(cè)二維正演模擬公式，并模擬了斷層、軟弱夾層的地震響應(yīng)。從中得到一些結(jié)論：

（1）利用三維VTI介質(zhì)波動(dòng)方程，根據(jù)隧道超前探測(cè)的觀測(cè)方式，推出了適用于隧道勘探正演模擬的二維VTI介質(zhì)qP方程，后導(dǎo)出其1階應(yīng)力-速度方程，進(jìn)行時(shí)間2階空間2N階的有限差分正演模擬。其計(jì)算迅速，規(guī)避了偽S波干擾和計(jì)算不穩(wěn)定性，便于快速得到不同地質(zhì)模型的地震響應(yīng)特征。

（2）利用PML吸收邊界，處理了隧道后方人工邊界引起的反射。采用交錯(cuò)網(wǎng)格和，使方程能較好地模擬斷層和軟弱夾層的地震響應(yīng)，減少頻散。

（3）二維平面模擬中，采樣時(shí)長(zhǎng)與斷層、軟弱夾層的物性參數(shù)組合不當(dāng)，會(huì)導(dǎo)致兩種不同的不良地質(zhì)體產(chǎn)生相同的地震響應(yīng)。在數(shù)值模擬時(shí)，可改變模型參數(shù)，使其顯現(xiàn)不同的地震響應(yīng)。而在現(xiàn)場(chǎng)實(shí)際情況來(lái)說(shuō)，需要現(xiàn)場(chǎng)工程師結(jié)合其他地質(zhì)勘探資料綜合分析。

參考文獻(xiàn)：

[1]柳楊春. HSP地質(zhì)超前預(yù)報(bào)技術(shù)及其應(yīng)用[J]. 西部探礦工程， 1997 （05）： 40-42.

[2]劉志剛， 劉秀峰. TSP（隧道地震勘探）在隧道隧洞超前預(yù)報(bào)中的應(yīng)用與發(fā)展[J]. 巖石力學(xué)與工程學(xué)報(bào)， 2003， 22?（8）：?1399-1402.

[3]曾昭璜. 隧道地震反射法超前預(yù)報(bào)[J]. 地球物理學(xué)報(bào)， 1994 （2）.

[4]鐘世航. 陸地聲納法的原理及其在鐵路地質(zhì)勘測(cè)和隧道施工中的應(yīng)用[J]. 中國(guó)鐵道科學(xué)， 1995 （4）： 48-55.

[5]Inazaki T， Isahai H， Kawamura S， et al. Stepwise application of horizontal seismic profiling for tunnel prediction ahead of the face[J]. The Leading Edge， 1999， 18（12）： 1429?1431.

[6]Otto R， Button E， Bretterebner H， Schwab P.The application of TRT-trae reflection tomography-at the Unterwald Tunnel[J]. Geophysics，?2002， 20 （2）： 51-56.

[7]Yamamoto T，?Shirasagi S，?Yokota Y，?Koizumi Y.?Imaging geological conditions ahead of a tunnel face using three-dimensional seismic reflector tracing system[J]. International Journal of the JCRM， 2011， 6（1）： 23-31.

[8]趙永貴. TST隧道超前預(yù)報(bào)技術(shù)及應(yīng)用[C].中國(guó)地球物理學(xué)會(huì)第22屆年會(huì)論文集，2006.

[9]Borm G， Giese R， Klose C， et al. ISIS integrated seismic imaging system for the geologic prediction ahead in underground construction[C].65th Annual Conference and Exhibition， EAGE， Extended Abstracts. Norway， 2003： 887?899.

[10]宋先海， 顧漢明， 肖柏勛. 我國(guó)隧道地質(zhì)超前預(yù)報(bào)技術(shù)述評(píng)[J]. 地球物理學(xué)進(jìn)展， 2006 （02）： 605-613.

[11]張平松，吳健生.中國(guó)隧道及井巷地震波法超前探測(cè)技術(shù)研究分析[J]. 地球科學(xué)進(jìn)展， 2006 （10）： 1033-1038.

[12]張學(xué)民. 巖石材料各向異性特征及其對(duì)隧道圍巖穩(wěn)定性影響研究[D]. 中南大學(xué)， 2007.

[13]程玖兵， 康瑋， 王騰飛. 各向異性介質(zhì)qP波傳播描述Ⅰ： 偽純模式波動(dòng)方程[J]. 地球物理學(xué)報(bào)， 2013， 56 （10）： 3474-3486.

[14]Alkhalifah， Tariq. An acoustic wave equation for anisotropic media[J]. Geophysics， 2000， 65（4）： 44-48.

[15]Thomsen， Leon. Weak Elastic Anisotropy[J]. Geophysics， 1986， 51：1954-1966.

[16]Clayton R， Engquist B. Absorbing boundary conditions for acoustic and elastic wave equations[J]. Bull. seismol. soc. amer， 1977， 67?（6）：?1529-1540.

[17]胡濤，任光明，聶德鑫， 等. 沉積型軟弱夾層成因及強(qiáng)度特征[J]. 中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào)，2004， 15 （1）： 124-128.