張煥鈞，陳祖斌，李 昊，楊興林

(吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院, 吉林 長春 130026)

0 引言

隧道施工時常發(fā)生塌方、突水等地質(zhì)災(zāi)害事故，需要有效的隧道超前探測手段來對未開采區(qū)域進行預(yù)報，以減少生命和財產(chǎn)的損失。TSP(tunnel seismic prediction)法是一種常用的利用地震波進行隧道探測的有效手段[1]，為了給其提供良好的模擬數(shù)據(jù)，并檢驗相應(yīng)TSP法反演的準確性和穩(wěn)定性，需要準確、高精度的地震數(shù)值正演模擬提供數(shù)據(jù)支撐。

目前，地球物理領(lǐng)域常用的正演方法有射線追蹤法[2]、積分方程法、有限元微分方程法[3]和有限差分微分方程法。其中，射線追蹤法較為直觀，但要求速度模型變化要平緩；積分方程法精度高，但數(shù)值求解困難；有限元微分方程法網(wǎng)格劃分靈活，但運行效率低下；而使用有限差分法求解波動方程模擬波場是計算效率和精度較高的正演模擬手段，非常適合實際工程應(yīng)用。20世紀80年代，Dablain[4]先提出了高階有限差分求解波場的方法；Jastram等[5]又提出了變步長有限差分法，提高了運算效率；之后，國內(nèi)學(xué)者[6-7]提出了交錯網(wǎng)格有限差分法，進一步提高了波場模擬的精度。隨著計算機效率的提高，近幾年學(xué)者們開始對復(fù)雜地質(zhì)結(jié)構(gòu)進行模擬，包括孔隙介質(zhì)[8]、黏彈性介質(zhì)[9]、各向異性介質(zhì)[10]等，還有一些學(xué)者專注于使用GPU加速正演算法的運行[11]，但與實際應(yīng)用場景(如隧道探測、礦井探測等)結(jié)合較少。

本文提出的隧道探測正演模擬服務(wù)于GEI-TSP隧道超前預(yù)報系統(tǒng)，該系統(tǒng)采用少炮多道的采集系統(tǒng)，具有成本低、安全性好的優(yōu)點。本文在該系統(tǒng)的基礎(chǔ)上進行數(shù)值模擬計算，并結(jié)合隧道模型給出相應(yīng)的模擬結(jié)果。

1 數(shù)值模擬

1.1 彈性波方程與交錯網(wǎng)格有限差分

地震波數(shù)值模擬可以采用聲波方程，如一階速度-應(yīng)力彈性波方程、二階位移彈性波方程等。由于隧道探測多為多波多分量地震信號的采集，對介質(zhì)的密度和速度較為敏感，且只需要探測掌子面后方平行或傾斜于掌子面的斷層、裂隙，故本文選擇二維一階速度-應(yīng)力彈性波方程進行模擬。假設(shè)介質(zhì)各向同性，密度不均勻，則其表達式為[12]：

(1)

式中：x為平行掌子面方向；z為垂直掌子面方向；ρ為介質(zhì)密度；vx、vz分別為速度在x、z方向上的分量；τxx、τzz分別為x、z方向的正應(yīng)力；τxz為切應(yīng)力；μ、Λ為拉梅系數(shù)，與傳播介質(zhì)的速度、密度相關(guān)，其表達式見式(2)—(3)。

Λ=ρ(vp2-2vs2)；

(2)

μ=ρvs2。

(3)

式(2)—(3)中：vp為介質(zhì)的縱波速度；vs為介質(zhì)的橫波速度。

對于式(1)，本文使用二維交錯網(wǎng)格有限差分法進行求解。有限差分的原理是用網(wǎng)格剖分的辦法將方程離散，用差分項代替微分項后再代入邊界條件進行求解，得到一個基于網(wǎng)格精度的數(shù)值解。而交錯網(wǎng)格的原理就是將速度、應(yīng)力定義在一個交錯了半個空間步長的網(wǎng)格上，這么做可以在保證精度的同時大大提高計算的效率，減少頻散的產(chǎn)生。交錯網(wǎng)格有限差分示意如圖1所示。

圖1 交錯網(wǎng)格有限差分示意圖

對于時間2階、空間2N階精度的有限差分法，式(1)中的空間微分項應(yīng)該用差分格式(4)來替代，時間微分項用差分格式(5)來替代。

(4)

(5)

式(4)—(5)中：f為需要進行微分的速度和應(yīng)力分量； Δx、Δz分別為x、z方向的空間步長；t為模擬時間； Δt為時間步長；an為差分系數(shù)，由式(6)確定。

(6)

式(4)的截斷誤差

Δf=C2Nf2N+1xΔx2N+1+O(Δx2N+2)。

(7)

1.2 有限差分數(shù)值模擬中的一些問題

1.2.1 穩(wěn)定性

差分方程的穩(wěn)定性是差分方程求解波動方程數(shù)值解的基本問題。交錯網(wǎng)格差分精度應(yīng)該滿足式(8)，才能使得數(shù)值計算過程穩(wěn)定。

(8)

式中γ為穩(wěn)定性閾值，其值見表1。

表1 γ與空間差分階數(shù)的關(guān)系[13]

1.2.2 頻散

頻散是差分方程近似微分方程過程中出現(xiàn)的不同波數(shù)的分量以不同速度傳播的現(xiàn)象，包括空間頻散和時間頻散。由于穩(wěn)定性(式(8))的要求，Δt與子波周期的比值很小，所以頻散主要由空間頻散造成[14]。交錯網(wǎng)格差分方程求解一階彈性波方程的空間頻散曲線公式如下：

(9)

式中p、q應(yīng)滿足條件見式(10)。

(10)

式中：v0為未發(fā)生頻散情況下的縱波速度；υ為泊松比；λ為波長；θ為傳播方向。

υ取0.25、θ取45°時，不同階數(shù)交錯網(wǎng)格差分方程的空間頻散曲線如圖2所示。可以看出，在頻散問題上，交錯網(wǎng)格法相比于常規(guī)網(wǎng)格有大幅度的改進，并且差分精度的提升也能明顯改善有限差分法的頻散問題。本文選擇100 Hz雷克子波作為震源，在最大縱波速度vp=4 000 m/s的情況下，波長λ為40 m，則Δx/λ=0.025，在這個步長波長條件下，頻散是十分微弱的。高階差分和交錯網(wǎng)格對于頻散問題改進的波場快照示意如圖3所示，可以看出提高階數(shù)和使用交錯網(wǎng)格法可以有效地改善頻散問題。

此外，本文還使用通量傳輸校正法[15]對頻散進行進一步的壓制。該方法假設(shè)波形中所有的極值點均是由數(shù)值頻散所引起的，在迭代求解波長過程中通過計算漫射通量對波形進行平滑處理并壓制頻散；之后，再計算反漫射通量補償平滑處理時的振幅損失。針對一階速度-應(yīng)力彈性波方程的通量傳輸矯正方法，只需要在速度項上進行校正，應(yīng)力項通過速度-應(yīng)力的關(guān)系自然會被矯正，這樣就提高了計算效率。

圖2 不同階數(shù)交錯網(wǎng)格差分方程的空間頻散曲線

1.2.3 邊界條件

由于計算機計算能力有限，只會在掌子面前方一定長度的區(qū)域進行數(shù)值模擬，此時人為劃分的邊界會產(chǎn)生本不存在的虛假反射，所以需要一定手段消除這一反射。Collino等[16]將Berenger提出的PML吸收邊界應(yīng)用在地震波領(lǐng)域，其原理是在波動方程中引入一個變化的衰減函數(shù)d(s)，在有限的吸收區(qū)域?qū)⑻摷俚姆瓷渌p到不會影響到正演區(qū)域的數(shù)值大小。該函數(shù)的系數(shù)部分應(yīng)與介質(zhì)速度成正比，與吸收厚度成反比。整個函數(shù)應(yīng)是一個與正演區(qū)域邊界距離s有關(guān)的凹函數(shù)，因為凹函數(shù)前期增長慢，與正演區(qū)域完美匹配，幾乎不產(chǎn)生反射，后期增長快，衰減效率高。此前學(xué)者[17-18]提出的衰減函數(shù)都依賴于經(jīng)驗，且衰減過程不平穩(wěn)，故本文提出一種新的衰減函數(shù)如下：

(11)

式中：L為吸收區(qū)域總層數(shù)；s為距正演區(qū)域的層數(shù)；R=0.000 1。

為了體現(xiàn)上述衰減函數(shù)的吸收效果，本文將其與目前廣泛使用的余弦衰減函數(shù)、指數(shù)衰減函數(shù)和二次冪衰減函數(shù)做對比測試。測試的區(qū)域是邊長為150 m的正方形，震源位置在區(qū)域的中心。4種衰減函數(shù)的波場快照如圖4所示，其中對振幅求取了絕對值，便于比較吸收效果。t=80 ms時4種衰減函數(shù)的波場快照中經(jīng)過吸收后剩余峰值的最大值和計算后的反射率見表2。

(a) 4階普通網(wǎng)格 (b) 10階普通網(wǎng)格

表2 4種衰減函數(shù)吸收效果對比

圖4所示的是中心震源激發(fā)后經(jīng)過4個邊和4個角吸收后反射回來的波。結(jié)合表2數(shù)據(jù)可以看出，幾種衰減函數(shù)都有效壓制了人為邊界造成的反射，但余弦函數(shù)和冪函數(shù)在衰減前中部反射過大，指數(shù)函數(shù)在衰減末端反射過大。而本文提出的衰減函數(shù)，衰減更為平緩，將進入吸收邊界的波的能量均勻吸收和反射，衰減后峰值最小，衰減效果最好。

1.2.4 計算效率

計算效率是在實際數(shù)值計算中不得不考慮的問題，尤其是在需要實時數(shù)據(jù)處理或者程序運行環(huán)境不佳的情況下。在正演模擬中，模擬區(qū)域的范圍(包括吸收邊界)與模擬的時長往往是固定的，此時數(shù)值計算效率就只與網(wǎng)格剖分精度Δx、Δz、Δt和差分階數(shù)2N(空間)、2M(時間)有關(guān)。

假設(shè)模擬區(qū)域?qū)挒閄、長為Z，模擬總時長為T，則數(shù)值計算的空間復(fù)雜度為O(n)，時間復(fù)雜度為O(m)，其中n滿足式(12)，m滿足式(13)，式中常數(shù)10代表了速度和應(yīng)力5個變量，而每個變量都在x和z方向上有分量。

(12)

(13)

雖然式(13)中數(shù)值計算的時間復(fù)雜度與差分階數(shù)和網(wǎng)格精度呈比例關(guān)系，但實際上由于穩(wěn)定性的要求，提高空間差分階數(shù)N或者縮小網(wǎng)格Δx、Δz都會進一步降低Δt，因此時間復(fù)雜度會進一步提高。這意味著過于追求計算精度可能會損失更多的計算效率。

根據(jù)以上4點問題的論述，本文算例采用時間2階空間10階交錯網(wǎng)格差分法。根據(jù)式(6)得到的空間10階精度下差分系數(shù)an見表3；根據(jù)穩(wěn)定性條件(式(8))，本文算例選擇空間網(wǎng)格剖分精度Δx=Δz=1 m，時間網(wǎng)格精度Δt=0.1 ms，模型最大縱波速度vp=4 000 m/s。

表3 10階精度差分系數(shù)

2 隧道超前探測正演算例

本文按照GEI-TSP隧道超前預(yù)報系統(tǒng)的要求，在掌子面附近安放炸藥震源，在隧道震源一側(cè)安放16道檢波器，最小炮檢距為15 m，道間距為2 m。單炮多道采集示意如圖5所示。本文以此為基準，分別建立層狀模型和空洞模型來驗證正演數(shù)值模擬的有效性。

圖5 單炮多道采集示意圖

2.1 層狀模型

表4 層狀模型參數(shù)

圖6 層狀隧道模型示意圖

圖6中正演區(qū)域長(z方向)300 m、寬(x方向)100 m；黑色部分為隧道，長45 m、寬10 m；紅色點為震源點，震源使用100 Hz雷克子波，震源有2 m的埋深；黃色區(qū)域安放檢波器共16個，最小炮檢距為15 m，檢波器與檢波器之間的距離為2 m。數(shù)值計算時，隧道的邊界設(shè)置為一個剛性反射界面，震源加載在彈性波方程的速度項上。層狀模型波場快照如圖7所示。

由圖7可知：t=30 ms時，地震波沿均勻介質(zhì)球面擴散，由于炸藥有2 m的埋深，會與隧道界面產(chǎn)生1次反射；t=80 ms時，地震波波前已經(jīng)過第1個傾斜反射界面，產(chǎn)生了反射波和透射波；t=120 ms時，可以清楚地看到第1個界面反射波和透射波均有P波和S波；t=170 ms時，地震波經(jīng)過第2個反射界面。

層狀模型共炮點道集如圖8所示，可以明顯地看到能量較強的直達波。將反射波波分放大后，也可以清楚地觀察到2個界面的反射P波和S波以及最初隧道界面反射在反射界面產(chǎn)生的反射波，尤其是在x分量上，反射S波更加明顯。

2.2 空洞模型

建立空洞隧道模型如圖9所示。這種模型可以代表隧道前方有一區(qū)域內(nèi)介質(zhì)的密度和地震波傳播速度極低，如溶洞、水體等。

圖9 空洞隧道模型示意圖

空洞位于掌子面正前方105 m的位置，形狀是一個直徑為15 m的圓，其內(nèi)部介質(zhì)密度和速度均為0。空洞外部為均勻介質(zhì)，密度ρ=2.8 g/cm3，縱波速度vp=3 700 m/s，其余參數(shù)與層狀模型一致?？斩茨Ｐ筒▓隹煺杖鐖D10所示。

(a) z分量(t=65 ms)

(b) z分量(t=120 ms)

(c) x分量(t=65 ms)

(d) x分量(t=120 ms)

由圖10可知，t=65 ms時，地震波沿均勻介質(zhì)球面擴散后，地震波波前與圓形空洞相接觸。在t=120 ms的波場中可以清晰地看到反射的P波和S波以及繞過空洞繼續(xù)向后傳播的彈性波。

空洞模型共炮點道集如圖11所示。圖中可以明顯地看到反射P波與反射S波，證明本文所建立的正演模型適用于介質(zhì)中含有空洞的情形。

(a) z方向分量

(b) x方向分量

3 結(jié)論與討論

波動方程正演是在多種因素制約下進行的波場模擬，其需要在兼顧穩(wěn)定性、頻散和精度的同時，盡可能地保證計算效率。本文提出的高階有限差分隧道超前探測正演模擬得到了符合預(yù)期的結(jié)果，并得出以下結(jié)論：

1)使用高階交錯網(wǎng)格有限差分求解一階速度-應(yīng)力彈性波方程，對隧道模型進行的模擬效果良好，所獲得的模擬數(shù)據(jù)可以為隧道超前探測的反演提供數(shù)據(jù)和模型支撐。

2)PML吸收邊界的衰減函數(shù)可以借鑒已有函數(shù)，但需要結(jié)合實際參數(shù)來修改衰減函數(shù)，這樣才能有好的吸收效果，從而減小吸收邊界寬度，提高計算效率。

使用有限差分求解波動方程的網(wǎng)格精度和差分階數(shù)應(yīng)該在滿足穩(wěn)定性要求的情況下盡可能地壓制頻散，但又不能過于犧牲計算效率。所以，如何在保證數(shù)值計算結(jié)果真實的情況下盡可能降低計算成本，是進一步研究所應(yīng)該探討的問題。