基于高階交錯網格有限差分的隧道超前探測地震波場模擬

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

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

0 引言

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

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

本文提出的隧道探測正演模擬服務于GEI-TSP隧道超前預報系統，該系統采用少炮多道的采集系統，具有成本低、安全性好的優點。本文在該系統的基礎上進行數值模擬計算，并結合隧道模型給出相應的模擬結果。

1 數值模擬

1.1 彈性波方程與交錯網格有限差分

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

(1)

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

Λ=ρ(vp2-2vs2)；

(2)

μ=ρvs2。

(3)

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

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

圖1 交錯網格有限差分示意圖

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

(4)

(5)

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

(6)

式(4)的截斷誤差

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

(7)

1.2 有限差分數值模擬中的一些問題

1.2.1 穩定性

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

(8)

式中γ為穩定性閾值，其值見表1。

表1 γ與空間差分階數的關系[13]

1.2.2 頻散

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

(9)

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

(10)

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

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

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

圖2 不同階數交錯網格差分方程的空間頻散曲線

1.2.3 邊界條件

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

(11)

式中：L為吸收區域總層數；s為距正演區域的層數；R=0.000 1。

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

(a) 4階普通網格 (b) 10階普通網格

表2 4種衰減函數吸收效果對比

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

1.2.4 計算效率

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

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

(12)

(13)

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

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

表3 10階精度差分系數

2 隧道超前探測正演算例

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

圖5 單炮多道采集示意圖

2.1 層狀模型

表4 層狀模型參數

圖6 層狀隧道模型示意圖

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

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

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

2.2 空洞模型

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

圖9 空洞隧道模型示意圖

空洞位于掌子面正前方105 m的位置，形狀是一個直徑為15 m的圓，其內部介質密度和速度均為0。空洞外部為均勻介質，密度ρ=2.8 g/cm3，縱波速度vp=3 700 m/s，其余參數與層狀模型一致。空洞模型波場快照如圖10所示。

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

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

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

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

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

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

(a) z方向分量

(b) x方向分量

3 結論與討論

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

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

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

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