隧道地震預報波場的有限元數值模擬

王朝令，劉爭平，黃云艷，黃顯彬

1.四川農業大學土木工程學院，成都 611830 2.西南交通大學土木工程學院，成都 610031

隧道地震預報波場的有限元數值模擬

王朝令1，2，劉爭平2，黃云艷1，黃顯彬1

1.四川農業大學土木工程學院，成都 611830 2.西南交通大學土木工程學院，成都 610031

在勘探地球物理領域中，能快速、有效地數值模擬2D和3D的地震全波場，包括同時模擬P波、S波、面波多分量波場特征的軟件并不是很多；但在結構力學領域中，已有多種成熟的多波多分量有限元數值模擬大型商用軟件用于求解彈性本構方程，如ANSYS，PLAXIS等。從方程結構來看，彈性本構方程與完全彈性波動方程相比，只多一項對時間的一階微分項——阻尼項。因此，調整合理的介質參數——Rayleigh系數，令阻尼項近似為0，彈性本構方程就蛻化為完全彈性波動方程。隧道地震波場是在巖體中傳播，具有彈性參數較好、無常規地震勘探中覆蓋層等低速帶干擾等的優點，采用ANSYS軟件對隧道地震波場進行數值模擬，研究了頻散、阻尼與吸收邊界等數值模擬的相關技術。針對介質吸收，比較了數值模擬中有無阻尼的時間記錄和波場快照，對隧道地震預報來說，將阻尼系數設定為0是一種合理的假設；比較了實例中不同網格長度與頻散的關系，當網格長度小于波長的1/π倍時，才能消除頻散；通過由波方程的推導，引入了黏彈性邊界條件，通過實例計算證明它可以有效地吸收邊界反射；最后對隧道地震預報進行了實例計算，算例表明采用ANSYS軟件可以有效地模擬隧道復雜地質條件下全波場的激發傳播過程。

數值模擬；ANSYS；有限元；隧道地震預報；地震波場；阻尼；頻散；邊界條件

0 引言

隧道屬于地下建筑工程，其圍巖地質情況復雜多變［1］，常遇到滑坡、崩坍、巖堆、巖洞、高應力高強度地層、松散地層、軟土等不良地質地段，以及膨脹地層、含水未固結圍巖、溶洞、斷層、巖爆、流沙和瓦斯溢出地層等特殊地質地段，它們具有不同的性質、形狀和規模，且成因和變異條件非常復雜，地質條件具有突發性。設計和勘察所提供的地質資料不可能自始至終符合實際情況，其詳細程度也不能完全達到隧道施工的具體要求，如果沒有可靠的超前地質預報，潛在的地質隱患就可能變成災難。因此，隧道施工中安全和工程事故時有發生，輕則延長工期、增加投資，重則毀壞機械設備，甚至造成人員傷亡，而且事故處理難度大。在隧道地震超前預報方面：20世紀90年代初，瑞士Amberg測量技術有限公司基于地震波反射原理開發研制了一套超前預報系統設備——隧道地震預報（TSP）系統，可探測和預報圍巖軟硬變化，斷層、破碎帶、節理裂隙發育情況，含水、溶洞及圍巖類別等［2］，TSP在瑞士、德國、法國、意大利、奧地利、日本等發達國家的隧道施工中，得到了廣泛的應用；Hanson等［3］采用真反射成像法（true reneetion tomo，TRT），該方法在觀測方式上采用的是空間多點接收和激發系統，檢波器和激發的炮點呈空間分布，以充分獲得空間波場信息，提高對前方不良地質體的定位精度。

地震模擬技術可以分為物理模擬和數值模擬2種［4］，數值模擬技術主要包括波動方程法和射線法兩大類。射線法模擬精度較低；波動方程法實質是求解地震波波動方程，它能比較準確地模擬任意復雜介質中的地震波場，為研究地震波傳播機理和復雜底層的解釋提供更多的數學及物理依據。鑒于隧道地震波場的復雜性，開展波場數值模擬技術是非常有意義的。

常用的地震波場數值模擬方法主要包括有限差分法（FDM）、有限元法（FEM）等，目前大多采用基于有限差分的彈性或黏彈性雙程波動方程進行計算。Boore［5］將有限差分法用于SH波非均勻介質地震波的CT成像模擬；Carcione等［6］提出了線性黏彈性介質中地震波傳播的模擬方法；Talezcr等［7］進行了線性黏彈性介質中地震波傳播的方法研究；Robertsson等［8］給出了黏彈性波有限差分模擬方法；Graves［9］給出了三維速度－應力方程交錯網格有限差分法彈性波傳播的模擬方法；Carcione等［10］研究了孔隙彈性介質中Biot縱波的數值模擬問題；Ozdenvar等［11］給出了孔隙彈性介質中地震波的交錯網格有限差分模擬方法；Carcione等［12］提出了孔隙黏彈性介質中地震波傳播的交錯網格有限差分模擬方法；何樵登等［13］研究了橫向各向同性介質中地震波數值模擬方法；王延光等［14］給出了一種適于強橫向變速的高階差分正演模擬方法；裴正林等［15］給出了三維任意各向異性介質中彈性波傳播的交錯網格高階有限差分模擬方法；董良國［16］利用交錯網格高階差分方法對起伏地表情況下的彈性波傳播進行了數值模擬；戴志陽等［17］利用褶積微分算子實現了對地震波場的數值模擬，并利用模型驗證了其有效性和正確性；周曉華等［18］利用交錯網格有限差分方法模擬了微動信號，結果表明采用這種方法得到的實測結果更符合實測結果。

這些方法具有較強的地震波傳播模擬能力，也能很好地描述小尺度不均勻各向異性介質（巖層的孔隙、裂隙等）中地震波的傳播特征。

在隧道波場數值模擬方面：Bohlen等［19－21］通過數值模擬在隧道中開展集成地震成像系統（ISIS）進行隧道超前預報，并通過3D有限差分的方法模擬了隧道超前預報，詳細地模擬了隧道邊墻上橫波與面波的轉換特征；劉晶波等［21］采用非線性有限元軟件LS-DYNA模擬了爆破對隧道結構的影響，并與現場爆破進行了對比，證明模擬方法的有效性；魯光銀等［22］采用高階交錯網格差分方法進行了隧道超前預報的全波場數值模擬計算。

目前，在勘探地球物理領域中，能快速、有效地數值模擬2D和3D地震全波場，包括同時模擬P波、S波、面波多波多分量波場特征的軟件，并不是很多。但在結構力學領域中，已有多種成熟的多波多分量有限元數值模擬軟件用于解彈性本構方程，如ANSYS，PLAXIS等。從方程結構來看，彈性本構方程與完全彈性波動方程相比，只多一項對時間的一階微分項——阻尼項。因此，調整合理的介質參數——Rayleigh系數，令阻尼項近似為0，彈性本構方程就蛻化為完全彈性波動方程。且這些軟件中模擬震源可以用任意自定義數據加載，因此可方便地用現場實驗和物理模擬提取的噪聲震源函數加載模擬，顯然比一般地球物理數值模擬中，主要用Ricker子波進行模擬更接近于實際情況。數值模擬中的關鍵技術是介質Rayleigh系數調整和吸收邊界條件的合理加載。

在本文中，筆者首先引入彈性力學的本構、Cauchy、Navier3個基本方程，這3個方程涉及應力、應變和位移3個物理量，是彈性波傳播理論中的基礎方程；然后以這3個方程構建波動方程；最后基于由以上步驟的理論準備，利用ANSYS有限元軟件實現對地震波場的數值模擬。

1 理論基礎

人工地震中，當介質（即物體）受到非零的外力（如人工激發震源）時，該外力要轉化為介質內部的應力，并使彈性介質內部發生應變和位移，形成地震波（即彈性波）場，介質的應力、應變以及能量的改變都是動態的運動過程。

1）本構方程的矩陣表達形式［23］為

式中：σ6×1和ε6×1分別為單角標形式的應力和應變矩陣；C6×6是物性矩陣，也稱為介質的彈性性質矩陣，它描述了介質的彈性和相關物理性質。

2）Cauchy方程的一維矩陣形式可以寫成

式中：u3×1是位移矩陣；L6×3是偏導數算符矩陣。

3）Naiver方程：

其中：ρ是介質密度；算符L2t表示對時間變量求二階偏導 數；ux，uy，uz表 示x，y，z方 向 的 位 移；［LxLyLz］T＝［?/?x?/?y?/?z］T；［FxFyFz］T表示力矢量。式（3）就是振動介質的“平動運動方程”，即Navier方程。

將Cauchy方程（2）代入本構方程式（1），再將得到的方程式代入Navier方程式（3）得

式中：位移矢量u的矩陣表達式為u3×1＝［uxuyuz］T；¨u3×1為節點的加速度矢量。式（4）即在外力作用下用位函數表示的波動方程。在地震波場中的外力就是震源的激發力，求解波動方程的問題就是波的傳播問題。

在有限元求解中，存在下列矩陣微分方程：

式中：m是質量矩陣；R是阻尼矩陣；K是剛度矩陣；˙u為節點速度矢量。式（5）是將單元剛度系數和單元力按照前面給出的方式相加得到的，而單元力包括外部力、初始應力等，此式相比于波動方程式（4）多了一個阻尼項R˙u，新的矩陣R和m也是按照一般方法由單元子矩陣組裝得到的，即式中：Re與me分別是單元阻尼矩陣和單元質量矩陣；Ni為形狀函數矩陣；μ為材料的黏滯矩陣；Ω為積分區域。

一般將式（5）定義的質量矩陣叫作一致性質量矩陣，是由Zienkiewicz和Cheung提出的［24］，其中的矩陣Re和R也稱為一致性阻尼矩陣。在計算過程中，集中質量矩陣更方便、更高效。實際上，阻尼矩陣R很難確定，通常將阻尼矩陣定義為剛度矩陣和質量矩陣的線性組合，即

式中，質量阻尼系數α和阻尼剛度系數β可通過實驗確定。這樣的阻尼稱為Rayleigh阻尼。

與波動方程式（4）形式類比，有限元問題方程（5）除去一阻尼項外，二者具有完全相同的形式，表明二者求解的彈性波場問題是等價的。

2 數值模擬實現策略

基于上述分析可知，瞬態動力學分析的基本方程式如式（5）所示，加入時間變量后變換為如下形式：

相比于式（5），將其中的力矩陣F變為荷載矢量F（t）。

在ANSYS軟件中，采用Newmark時間積分方法（包括改進的算法HHT）和中心差分時間積分法求解瞬態動力學基本方程。

在采用ANSYS軟件實現彈性波場模擬中，主要牽涉以下幾個方面的問題：1）單元選擇；2）介質吸收；3）頻散效應；4）邊界條件。下面分別進行討論。

2.1 單元選擇

由于計算機能力的限制，筆者首先研究2D的數值模擬，ANSYS中針對結構的二維單元包括Plane25、Plane42、Plane83、Plane182、Plane183。其中：Plane183是高階二維八節點單元，具有二次位移項，適于生成不規則網格模型，因此不選用此單元；Plane25、Plane83均是諧響應對稱單元，在本文研究的模型中，由于震源加載的限制，不能滿足對此類單元加載載荷的要求，因此不選擇此2種單元。比較之下，可供選擇的平面單元只有2個：Plane42和Plane182。

Plane42和Plane182單元非常相近。Plane182單元具有選擇簡化積分這一選項，這對于不可壓縮情況有助于防止網格被鎖死；此外Plane182單元具有4個用戶不可調節的自由度，相比于Plane42單元，在計算中具有較高的精確度。在本文模擬中，選擇Plane182單元進行彈性波場的數值模擬。圖1是在分別基于Plane42和Plane182單元采用相同震源激發等條件時接收到的地震記錄，Plane182單元計算結果更好，但這種差異是比較微小的。

圖1 不同單元數值模擬Fig.1 Modeling in different elements

2.2 介質吸收

ANSYS可以進行加載阻尼的運算。雖然筆者在本文中只研究彈性介質，不需要考慮阻尼因素，但當進行地震波場模擬時，需要考慮Rayleigh阻尼所造成的介質吸收問題，否則會出現意想不到的結果。據式（8）知，要計算Rayleigh阻尼，必須定義α和β；但一般情況下，α和β不能直接獲得，通常由模態阻尼比ξ計算得出。模態阻尼比是特定振動模式下實際阻尼與臨界阻尼的比值。若ωi是第i階的自然圓頻率，那么α和β滿足如下關系式：

關于模態阻尼比ξ，不少文獻中的研究表明取0.02～0.05是比較合適的；在一次模擬中，往往取一個常數；對于巖石等，一般取0.05。在許多實際問題中，往往忽略質量阻尼系數α。因此，將ωi＝2πfi代入式（11）可以得到

式中：f為頻率。

若加載主頻為100Hz，則由式（13）可以得到剛度阻尼系數β為1.59×10－4。大部分體系的剛度阻尼系數都是落在0～0.05這一范圍的，在此范圍內進行動力學計算時，位移的變化不超過5%；β為零時，系統的位移最大。因此，出于保守的目的，當系統承受沖擊荷載時，取β為0是很合理的逼近。

下面通過數值計算來研究在加載阻尼系數β為1.59×10－4和0（即不加載阻尼系數）時的模擬效果。模型如圖2所示，網格長度為0.5m，震源主頻為300Hz，提取偏移距為10m、道間距1m的時間記錄。

圖2 阻尼加載模型圖Fig.2 Model figure of damping loading

圖3是加載阻尼和無阻尼情況下的波場快照。可以看出：2種情況下傳播時的差異比較小，且相同時刻波場也非常接近；圖3a中由于施加了阻尼，在波場擴散中能量不斷減小，對比50ms時刻的波場快照尤其明顯。

圖4為提取的時間記錄，用來比較加載前后的影響。由圖4可以看出，加載阻尼和無阻尼的時間記錄非常接近，無阻尼時得到的剖面信噪比更高，反射波的同相軸更清晰。

由此可見，ANSYS可以很好地模擬阻尼介質的波場。在地震超前預報的數值模擬中，將質量阻尼系數α和阻尼剛度系數β都設定為0。

2.3 網格剖分的頻散效應

對模型進行網格劃分時，由于要進行瞬態分析來模擬波場傳播，必須考慮由于波長與網格長度所產生的頻散效應，這主要取決于網格長度和時間步長的大小。從波動力學帶寬定理［25］得到

圖3 阻尼波場快照Fig.3 Snapshot of damping

式中：k為波數；Δt為時間步長；Δx為網格長度。式（15）表明，能量所在空間或時間越集中，其頻譜就越寬。

國內外有不少學者從不同角度提出離散化準則，如 Kausel［26］，Costantino等［27］和 Miller等［28］分別給出：

式中：λ為波長；λmin為最小波長。

圖4 阻尼加載時間記錄Fig.4 Time recording of damping loading

據宗福開［25］的研究結果，由λmin＝v/fmax（v為速度，fmax為最高頻率），當波在傳播方向劃分的網格數n大于10時，則得出有限元離散化準則為

式（18）具有明確的物理意義，即要求網格長度小于波長的1/π倍。應當指出，根據波傳播邊界條件的性質，通常Δx≤（1/π～1/8）λmin。

據以上結論，筆者進行了網格長度不同的有限元數值模擬，模型示意圖如圖5所示。為減小邊界對模擬結果的影響，模型4個邊界均加載黏彈性吸收邊界。模型背景設定橫波波速為2 000m/s，縱波波速為4 000m/s，密度為1 800kg/m3。由實際資料分析，震源主頻取300Hz。由于橫波的波長較縱波的小，故模型的λmin取為橫波的主波長，約為20 m。分別設定網格長度為10.0、5.0、2.5、1.0和0.5 m，分析波場數值模擬的頻散效應。

圖6是不同網格長度、同一節點的位移數值。根據公式（15），當網格長度為10.0m＞20/π＝6.366m時，波場的位移數值解震蕩很嚴重。隨著網格長度越來越短，波的震蕩現象逐漸減弱。當網格長度為5.0m時，相比于10.0m的計算值，波的震蕩變小了；當網格長度為2.5m時，仍然存在輕微的震蕩；當網格長度由1.0m進一步變小到0.5 m時，模擬計算結果并沒有太大的改進。考慮到計算量的問題，對于實際模型模擬時，采用最大網格長度為1.0m的剖分網格就可以最大程度地減小頻散效應，很好地模擬波場的傳播了。因此，在進行模擬時，按照式（18）所要求的網格長度進行劃分，即可消除頻散的影響。

圖5 網格剖分模型圖Fig.5 Model figure of grid meshing

圖6 采樣點位移隨網格長度變化圖Fig.6 Displacement of sampling point varied with meshing

2.4 邊界條件

采用有限元模擬時，通常采用時間域的吸收邊界，這種邊界由一維或二維條件導出，并使得沒有能量從邊界反射。其中：應用最廣泛的是Lysmer等［29］引入的黏彈性邊界，它采用黏彈性阻尼來代替遠場邊界；Novak等［30］提出了平面應變邊界，這種邊界主要應用到嵌入基礎和樁基中。但是，平面應變邊界包含依賴于頻率的項，這在加載瞬態載荷的瞬態響應分析中會變得比較復雜；相比之下黏彈性邊界應用更廣泛，加載簡單，且對整個計算所增加的計算量可以忽略不計，在采用黏彈性邊界時，極少會產生這種反射。

當波穿過介質時，波前的應力產生在法向和切向2個分量上。假定在平面應力條件下，徑向位移的波動方程為［31］

其中：ur表示徑向位移；G為壓縮模量；r為半徑；Λ為Lame常數。

式中：E為Young模量；v是Possion比。

將徑向位移ur用位移勢函數φ表示：

則式（19）可表示為

兩邊都對徑向r進行積分，則波方程用φ表示為

可以得到估計表達式：

用f′表示頻率f的導數，則徑向位移ur、徑向應力σr、徑向應變εr、橫向應變εθ可以用頻率f和半徑r表示：

對于任意的半徑rb，對頻率f的任意時刻導數都是其對幅角導數的相反值，則半徑rb處的位移值、速度值和加速度值可以表示為

在半徑rb處二階時間導數為

將式（33）代入式（29）可以得到f（rb，t）在邊界的應力：

對式（34）做時間微分：

聯立式（34）和式（35），可以消去未知函數f：

此種在有限元中集合時間的邊界，在實現中，若采用彈簧形式將會非常方便。Wolf［32］提供了如圖7所示的彈簧。

圖7 彈簧模型Fig.7 Spring model

圖7中彈簧模型的波動方程如下：

由式（37）可以得到

式（39）對時間求導：

將式（39）和式（40）代入到式（38）可以得到第一個自由度的波動方程：

聯立式（36）—式（41）可以得到等價的剛度、阻尼和質量：

可以將式（42）（43）（44）所得的值賦予彈簧，來實現黏彈性邊界。

歸納二維人工邊界等效物理系統的彈簧系數K和阻尼系數C在針對不同邊界時分別如下。

切向邊界：

法向邊界：

式中：KBN、KBT分別為彈簧的法向與切向剛度；CBT、CBN分別為彈簧的法向與切向阻尼；rb為波源至人工邊界點的距離；γT和γN分別為切向與法向黏彈性人工邊界參數。根據劉晶波等［33］的研究結果，人工邊界比較合適的參數γT取值范圍為［0.35，0.65］，γN的取值范圍為［0.8，1.2］。

在ANSYS軟件中，適合用來做黏彈性邊界的單元選用Combin14，它在1維、2維或3維應用中具有軸向或扭轉的性能。縱向阻尼彈簧是單軸壓縮張力的單元，在每個節點有3個自由度：x、y和z軸。不考慮扭轉和彎曲，另外彈簧阻尼單元沒有質量矩陣，圖8為Combin14單元的示意圖和此單元在ANSYS數值模擬計算中的加載方式。

圖8 Combin14單元示意圖（a）和邊界加載方式（b）Fig.8 Combin14unit schematic（a）and boundary loading mode（b）

圖9 邊界比較模型圖Fig.9 Model figure of boundary comparing

按照式（19）和式（21）將Combin14在模型周邊布設，利用如圖9所示的模型計算四周是Dirichlet邊界和黏彈性邊界的數值模擬。震源為100Hz的Ricker子波，加載方向豎直向下，網格長度0.5m，觀測系統以震源為中心，兩邊等間距，道間距為2 m，道數為48道，分別取得震源兩邊的時間記錄，如圖10所示。

在模型四周施加了黏彈性邊界時，直達P波和直達S波都受到了抑制，且相較之下，S波受影響更大。這是由于黏彈性邊界彈簧的存在使除直達P、S波之外的邊界反射被吸收的緣故，證明了黏彈性邊界可以有效地吸收反射。在默認位移為0的Dirichlet邊界所得到的記錄中（圖11），各邊界的反射波返回到接收排列，使得對波的辨識和同相軸的判別都變得比較困難。由圖10時間記錄和圖11波場快照的相同時刻的對比可以看出，Dirichlet邊界的波場快照中各種反射波夾雜在一起，使得波場變得非常復雜，不易辨識，增加了波場分析的難度。

圖10 2種邊界的時間記錄Fig.10 Time record of two boundaries

圖11 2種邊界的波場快照Fig.11 Wavefield snapshot of two boundaries

3 計算實例

圖12為邊墻激發二維隧道數值模型。模型隧道開挖工作面前方具有一模擬斷層垂直傾角的低波阻抗反射帶。隧道模型區域I中的縱、橫波速度分別設為4 000m/s和2 000m/s，低波阻抗區域II的縱、橫波速度為2 000m/s和1 000m/s，另一側區域III的縱、橫波速度設為3 000m/s和1 500m/s；其觀測系統中激發震源位于邊墻上，接收排列的48道檢波器也設置在邊墻上，并與震源在同一直線上，最小偏移距10m，道間距1m。由互換原理可知，若將接收與激發換位，即為TSP觀測系統。模擬震源為主頻300Hz的Ricker。采樣間隔為0.1ms，采樣視窗長0.1s，網格長度0.5m。模型邊界均采用黏彈性邊界條件。

圖12 邊墻激發數值模型Fig.12 Numerical model of sidewall excited

圖13 實例波場快照Fig.13 Snapshot of the example

圖13是邊墻激發所產生的波場快照。由圖13可見，震源激發的橫波能量遠大于縱波能量，所以波場快照中縱波傳播和反射特征表現不如橫波那么明顯。由圖13分析，P波和S波的隧道波場特征有很大的差異，因此對二者的特征分述如下。首先，當t＝5ms時，傳播中的P、S波還沒有完全分離；t＝10ms、t＝15ms時，由于P波比S波速度快，二者已經完全分離。對于P波：當t＝23ms時，P波傳播到掌子面；t＝30ms時，P波經過掌子面，并繼續向前方介質傳播；t＝37ms時，P波到達第一波阻抗變化界面，部分能量反射回來，另一部分則透射過界面繼續向前傳播，反射回來的P波將會被沿邊墻布置的檢波器觀測排列接收。對于S波，激發的S波一部分能量以體波的形式向邊墻外的介質體內傳播，另一部分能量則轉換為瑞利（Rayleigh）面波，沿邊墻向前傳播。t＝40ms時，S波到達掌子面。t＝46ms、t＝50ms時，瑞利面波的部分能量在掌子面轉換為S波，新產生的S波以體波的波陣面繼續向前傳播［20］，部分能量則繼續以面波的波陣面沿邊墻向后傳播，可被沿邊墻布置的檢波器觀測排列接收。這類來自于掌子面的反射面波被稱為SRR波，顯然對于提取掌子面前方介質的反射S波信息來看，是干擾波。t＝55ms、t＝58ms時，轉換S波繼續向前傳播，直到t＝75ms時刻，轉換S波到達第一反射界面，部分能量發生反射回傳，部分發生透射繼續向前傳播。反射回傳的S波到達隧道邊墻時，部分能量將再次轉換為面波，并會被沿邊墻布置的檢波器觀測排列接收。分析表明，隧道邊墻激發條件下，觀測排列記錄的反射S波實際上是由S波轉換為面波，面波再轉換為S波，最后再轉換為面波的復雜轉換波，被稱為SRSR波。波場快照還顯示，SRSR波的能量遠大于反射P波。

圖14是邊墻激發數值模型位移x分量、y分量的時間記錄。圖14中的同相軸主要包括反射縱波、反射橫波和面波轉換橫波，與圖13的波場快照相對應，圖14中的負視速度同相軸就是面波在掌子面轉換為橫波傳播到接收排列所形成的同相軸SRSR波。

圖14 實例時間記錄Fig.14 Profiles of the example

以上分析表明，有限元數值方法可以有效地模擬隧道復雜地質條件下的全波場的激發傳播過程。

4 結論

1）由于隧道地震預報的波場傳播主要是在巖體中進行，相比于常規地面地震勘探，具有更好的彈性參數。對隧道地震波場的模擬而言，將質量阻尼系數α和阻尼剛度系數β都設置為0是比較合理的近似，能滿足數值模擬的需要；若模擬地表地震勘探，則需要考慮介質阻尼。

2）針對網格長度，為滿足計算要求，根據網格長度不大于（1/π～1/8）λmin的原則，針對主頻為300 Hz的震源而言，取網格長度為1m就可以消除頻散效應的影響。

3）加載黏彈性邊界條件可以很好地吸收邊界反射，并且其有良好的魯棒性。

4）算例表明，采用大型有限元軟件可以有效地模擬隧道復雜地質條件下的全波場的激發傳播過程。

5）本文只是對二維隧道地震波場進行了數值模擬研究，相比于實際需求，開展三維隧道地震波場的數值模擬研究更具實際意義。

本文采用數值模擬方法雖然只研究了隧道地震反射波場的傳播規律，根據不同模型，它也可以模擬地震勘探領域中其他地震波場的傳播規律，這為地震波場的數值模擬提供了一種解決手段。

（References）：

［1］魯光銀.隧道地質災害反射波法探測數值模擬及圍巖f－Ahp分級研究［D］.長沙：中南大學，2009.

Lu Guangyin.Numerical Simulation for Tunnel Geological Hazard Detection by Reflected Wave Method and Rock Dynmanic Classification by f-Ahp Method Research［D］.Changsha：South China University，2009.

［2］劉志剛，劉秀峰.TSP隧道地震勘探在隧道隧洞超前預報中的應用與發展［J］.巖石力學與工程學報，2003，22（8）：1399-1402.

Liu Zhigang，Liu Xiufeng.TSP Application and Development in Tunnel Lead Forecast［J］.Chinese Journal of Rock Mechanics and Engineering，2003，22（8）：1399-1402.

［3］Hanson D R，Haramy K Y，Neil D M.Seismic Tomography Applied to Site Characterization［M］.Denver：American Society of Civil Engineers，2000.

［4］佘德平.波場數值模擬技術［J］.勘探地球物理進展，2004，27（1）：16-21.She Deping.On Wave Field Numerical Modeling Techniques［J］.Progress in Exploration Geophysics，2004，27（1）：16-21.

［5］Boore D M A.Note on the Effect of Simple Topography on Seismic SH Waves［J］.Bulletin of the Seismological Society of America，1972，62（1）：275-284.

［6］Carcione M，Kosloff，Kosloff R.Wave Propagation Simulation in a Linear Viscoelastic Medium［J］.Geophysical Journal，1988，95（3）：597-611.

［7］Talezer H，Carcione J M，Kosloff D.An Accurate and Efficient Scheme for Wave-Propagation in Linear Viscoelastic Media［J］.Geophysics，1990，55（10）：1366-1379.

［8］Robertsson.Numerical Free-Surface Condition for Elastic Viscoelastic Finite-Difference Modeling in the Presence of Topography［J］.Geophysics，1996，61（6）：1921-1934.

［9］Graves R W.Simulating Seismic Wave Propagation in 3DElastic Media Using Staggered-Grid Finite Differences［J］.Bulletin of the Seismological Society of America，1996，86（4）：1091-1106.

［10］Carcione J，Quiroga-Goode G，Cavallini F.Wavefronts in Dissipative Anisotropic Media：Comparison of the Plane-Wave Theory with Numerical Modeling［J］.Geophysics，1996，61（3）：857-861.

［11］Ozdenvar T，McMechan G.Algorithms for Staggered-Grid Computations for Poroelastic，Elastic，Acoustic，and Scalar Wave Equations［J］.Geophysical Prospecting，1997，45（3）：403-420.

［12］Carcione J，Helle H.Numerical Solution of the Poroviscoelastic Wave Equation on a Staggered Mesh［J］.Journal of Computational Physics，1999，154（2）：520-527.

［13］何樵登，張中杰.3D橫向各向同性介質中的體波［J］.石油地球物理勘探，1990，25（2）：137-147.

He Qiaodeng，Zhang Zhongjie.3DTransversely Isotropic Medium Body Wave［J］.Oil Geophysics Exploration，1990，25（2）：137-147.

［14］王延光，李振春，穆志平，等.一種適于強橫向變速的高階差分正演模擬方法［J］.石油大學學報：自然科學版，2002，26（5）：37-39.

Wang Yanguang，Li Zhenchun，Mu Zhiping，et al.A High-Order Finite-Difference Forward-Modeling Fit for Strong Lateral Velocity Variation［J］.Journal of the University of Petroleum，China，2002，26（5）：37-39.

［15］裴正林.地震波標量方程的小波自適應網格有限差分法數值模擬［J］.石油地球物理勘探，2005，40（3）：267-272.

Pei Zhenglin.Wavelet Scalar Adaptive Mesh Seismic Wave Equation Finite Difference Numerical Simulation［J］.Oil Geophysics Exploration，2005，40（3）：267-272.

［16］董良國.復雜地表條件下地震波傳播數值模擬［J］.勘探地球物理進展，2005，28（3）：187-194.

Dong Liangguo.Complex Surface Conditions Numerical Simulation of Seismic Wave Propagation［J］.Progress in Exploration Geophysics，2005，28（3）：187-194.

［17］戴志陽，孫建國，查顯杰.地震波場模擬中的褶積微分算子法［J］.吉林大學學報：地球科學版，2006，35（4）：520-524.

Dai Zhiyang，Sun Jianguo，Zha Xianjie.Seismic Wave Field Modeling with Convolutional Differentiator Algorithm［J］.Journal of Jilin University：Earth Science Edition，2006，35（4）：520-524.

［18］周曉華，陳祖斌，曾曉獻，等.交錯網格有限差分法模擬微動信號［J］.吉林大學學報：地球科學版，2012，42（3）：852-857.

Zhou Xiaohua，Chen Zubin，Zeng Xiaoxian，et al.Simulation of Microtremor Using Staggered-Grid Finite Difference Method［J］.Journal of Jilin University：Earth Science Edition，2012，42（3）：852-857.

［19］Bohlen T.Parallel 3-D Viscoelastic Finite Difference Seismic Modelling［J］.Computers ＆ Geosciences，2002，28（8）：887-899.

［20］Bohlen T，Lorang U，Rabbel W，et al.Rayleigh-to-Shear Wave Conversion at the Tunnel Face：From 3D-FD Modeling to Ahead-of-Drill Exploration［J］.Geophysics，2007，72（6）：67-79.

［21］Liu Jingbo，Yan Qiushi，Wu Jun.Analysis of Blast Wave Propagation Inside Tunnel［J］.Transactions of Tianjin University，2008，14（5）：358-362.

［22］Lu Guangyin，Han Xuli，Zhu Ziqiang.A Staggered-Grid High-Order Finite Difference Numerical Simulation for Tunnel Seismic Prediction Ahead of Construction［J］.Engineering Computation，2009，23（2）：225-228.

［23］牛濱華.半空間均勻各向同性單相固體彈性介質與地震波傳播［M］.北京：地質出版社，2005：200-205.

Niu Binhua.Single-Phase Half-Space Homogeneous Isotropic Elastic Solid Medium and Seismic Wave Propagation［M］.Beijing：Geological Publishing House，2005：200-205.

［24］Cheung Y K，Jin W G，Zienkiewicz O C.Direct Solution Procedure for Solution of Harmonic Problems Using Complete，Non-Singular，Trefftz Functions［J］.Communications in Applied Numerical Methods，1989，5（3）：159-169.

［25］宗福開.波傳播問題中有限元分析的頻散特性及離散化準則［J］.爆炸與沖擊，1984，4（4）：16-23.

Zong Fukai. Frequency-Dispersion Characteristics and Dicretization of the Finite Element Analysis in Wave Propagation Problems［J］.Explosion and Shock Waves，1984，4（4）：16-23.

［26］Kausel E.Local Transmitting Boundaries［J］.Journal of Engineering Mechanics，1988，114（6）：1011-1027.

［27］Costantino C，Miller C，Lufrano LA.Soil Structure Interaction Parameters from Finite Element Analysis［J］.Nuclear Engineering and Design，1976，38（2）：289-302.

［28］Miller T F，Schmidt F W.Use of a Pressure Weighted Interpolation Method for the Solution of the Incompressible Navier-Stokes Equations on a Nonstaggered Grid System［J］.Numerical Heat Transfer：Part A：Applications，1988，14（2）：213-233.

［29］Lysmer J，Kuhlemeyer R L.Finite Dynamic Model for Infinite Media［J］.Journal of Engineering Mechanics，1969，95（4）：859-877.

［30］Novak M，Hindy A.Seismic Analysis of Underground Tubular Structures［C］//Proceedings of the 7th World Conference on Earthquake Engineering.Istanbul：Ankara Natl Common Earthquake Eng，1980：287-294.

［31］Deeks J，Randolph M.Axisymmetrical Time-Domain Transmitting Boundaries［J］.Journal of Engineering Mechanics-Asce，1994，120（1）：25-42.

［32］Wolf J P.A Comparison of Time-Domain Transmitting Boundaries［J］.Earthquake Engineering ＆Structural Dynamics，1986，14（4）：655-673.

［33］劉晶波，谷音，杜義欣.一致粘彈性人工邊界及粘彈性邊界單 元 ［J］.巖 土 工 程 學 報，2006，28（9）：1070-1076.

Liu Jingbo，Gu Yin，Du Yixin.Consistent Viscous-Spring Artificial Boundaries and Viscous-Spring Boundary Elements［J］.Chinese Journal of Geotechnical Engineering，2006，28（9）：1070-1075.

Wavefield Modeling Based on the Finite Element Method for the Tunnel Seismic Prediction

Wang Zhaoling1，2，Liu Zhengping2，Huang Yunyan1，Huang Xianbin1

1.CollegeofCivilEngineering，SichuanAgricultureUniversity，Chengdu611830，China
2.CollegeofCivilEngineering，SouthwestJiaotongUniversity，Chengdu610031，China

In the field of exploration geophysics，the software packages that can fast and efficiently simulate 2Dand 3Dseismic full-wavefield，including simultaneously simulate P-wave，S-wave，and surface wave are not very much.In the structure mechanics，however，many sophisticated commercial software packages based on the finite element method including ANSYS and PLAXIS for multi-wave and multi-component have been employed for the solutions of the elastic constitutive equations.Compared to the complete elastic wave equations，the elastic ones include a damping term，which is a first-order differential term time.Therefore，by adjusting the parameter of the medium，that is，the Rayleigh coefficient to make the damping to zero，the elastic constitutive equations will degenerate into full elastic wave equation.The wavefields propagating in the rocks of a tunnel are free of the disturbance of low velocity resulting from the cover layer in the conventional seismic exploration.We simulated the wavefields of the tunnel seismic prediction using the ANSYS software，and researched the dispersion，damping and absorbing conditions etc.in the numerical modeling.For medium absorption，we compared time records and snapshots of the wavefields in the numerical simulation for the cases with and without the damping.In the tunnel seismic prediction，setting the damping to zero is a reasonable assumption.By comparing the relation of different grid length and dispersion in the numerical example，we find that the dispersion disappears when the mesh size is smaller than the wavelength of 1/π.By introducing a viscoelastic boundary condition deduced from the wave equation，boundary reflections can be effectively absorbed.Finally，an numerical example of the tunnel seismic prediction shows that the usage of the ANSYS software can effectively simulate the propagation of waves in the complicated geological conditions.

numerical modeling；ANSYS；FEM；tunnel seismic prediction；seismic wavefield；damping；dispersion；boundary condition

10.13278/j.cnki.jjuese.201404305

P631.4

A

王朝令，劉爭平，黃云艷，等.隧道地震預報波場的有限元數值模擬.吉林大學學報：地球科學版，2014，44（4）：1369-1381.

10.13278/j.cnki.jjuese.201404305.

Wang Zhaoling，Liu Zhengping，Huang Yunyan，et al.Wavefield Modeling Based on the Finite Element Method for the Tunnel Seismic Prediction.Journal of Jilin University：Earth Science Edition，2014，44（4）：1369-1381.doi：10.13278/j.cnki.jjuese.201404305.

2013-12-12

國家自然科學基金項目（40874051）

王朝令（1980—，男，講師，博士，主要從事地震波場數值模擬、隧道超前預報預報研究，E-mail：wong81010＠gmail.com。