地震聲波數值模擬中人工邊界條件的差別與組合

韓復興，王若雯，孫章慶，高正輝

吉林大學地球探測科學與技術學院，長春 130026

0 引言

在地震聲波數值模擬計算過程中,選取的有限計算區域會產生邊界反射，干擾正常波場計算，因此需要引入人工邊界條件降低邊界反射的影響。人工邊界條件按方法原理可以分為吸收類邊界條件和衰減類邊界條件兩大類。吸收類邊界條件通過在邊界處使用單程波動方程模擬地震波的傳播，使入射波只向計算區域外傳播而不會產生邊界反射；基于傍軸近似原理的Clayton-Engquist(CE)邊界條件是吸收類邊界條件中最經典的一種邊界條件。衰減類邊界條件的方法原理是在計算區域外引入一層衰減區域，入射波在衰減層內傳播時，入射波能量逐漸衰減從而不會產生邊界反射；其中完全匹配層(PML)邊界條件應用效果最優，近年來得到廣泛應用[1]。

近年來許多相關專家學者都對CE邊界條件和PML邊界條件的應用效果進行了比較，包括定性比較和定量比較。定性比較主要通過觀察波場快照或合成地震記錄的邊界反射細節來比較邊界條件的吸收效果。如：王守東[2]導出了完全匹配層法的聲波方程表達形式，并通過增強地震記錄振幅的顯示方式得出PML邊界條件吸收效果優于CE邊界條件的結論；邢麗[3]通過數值算例分析認為，CE邊界條件在程序方面易于實現，但1階CE邊界條件精度低，2階CE邊界條件有時又會出現不穩定的現象，PML邊界條件雖然程序實現較為繁瑣，但計算過程穩定，選取適當的參數時，幾乎可以吸收全部的邊界反射；王永剛等[4]分別采用井間和地面地質模型進行數值模擬，結果表明PML邊界條件相對于旁軸近似法具有優越性；裴俊勇等[5]利用CE邊界條件和PML邊界條件進行有限差分正演模擬，模擬圖像觀察結果說明PML邊界條件吸收邊界效果較好。

定量比較通過量化邊界反射數值進行，主要方法為根據反射波峰值高低比較邊界條件吸收效果強弱。如，付小波等[6]對CE邊界條件和PML邊界條件進行了精細的定量比較性研究，發現CE邊界條件相對PML邊界條件受入射波角度限制更大，當入射角增大時，CE邊界條件吸收效果大大降低，PML邊界條件則不受入射波角度的限制。

本文通過對PML邊界條件和CE邊界條件進行詳細的分析研究，考慮到其各自的優缺點，將CE邊界條件和PML邊界條件組成新的組合邊界，并以數值模擬結果驗證算法的有效性，以期達到在減少衰減帶厚度的同時提高計算效率的目的。

1 PML邊界條件兩種差分形式對比

PML衰減介質中的波動方程可以看作是常規波動方程的推廣，入射波傳播到介質中振幅呈指數衰減。當衰減因子和衰減帶厚度選擇適當時，理論上入射波不會再反射回計算區域[7]。

PML的實現形式有分裂式(SPML)和非分裂式(NPML)兩種[8-9]。SPML是原始的PML形式，其是在分裂原始波場分量的過程中加入PML衰減因子得出基于PML邊界的波動方程。在計算區域和衰減區域都使用基于PML邊界的波動方程，將計算區域的PML衰減因子設置為0，在衰減區域設置衰減因子，這種方式稱為全局SPML。全局SPML雖然數據存儲量增大，但相對而言編程實現簡單[10]，也是本文考慮和使用的PML形式。在計算區域使用加入PML衰減因子前的普通波動方程、在衰減區域使用基于PML邊界波動方程的SPML形式則為局部SPML。這種形式節省了大量的數據存儲量，但要考慮不同區域的PML波動方程，編程難度也有所增加。NPML不需要分裂波場，有一定的優勢，但在計算過程中要進行復雜的卷積計算，計算效率低[11-12]。

本文將討論在地震聲波有限差分數值模擬過程中實現全局SPML的非解耦與解耦差分形式的基本原理和實現效果。

1.1 全局SPML非解耦差分形式

均勻介質中二維聲波波動方程為

(1)

式中：u為聲波波場；(x,z)為空間坐標；v為聲波波速；t為時間。孫林潔等[13]推導出基于PML邊界的三維聲波波動方程，根據其推導過程可以得到基于PML邊界的二維聲波波動方程[14]：

(2)

式中：d(x,z)為任意點(x,z)的衰減因子。對式(2)進行時間2階、空間2m階的有限差分。設有限差分空間采樣步長為Δx、Δz，時間采樣步長為Δt，x=iΔx、z=jΔz，t=kΔt，可以得到

(3)

1.2 全局SPML解耦差分形式

(4)

(5)

式中：ux、uz分別為x、z方向的波場值；dx、dz分別為x、z方向的衰減因子。

將式(4)(5)作向前差分，可以得到有限差分格式：

(6)

(7)

(8)

(9)

(10)

1.3 解耦與非解耦PML差分形式數值模擬對比分析

設定Δx=Δz=6 m，Δt=0.8 ms，總采樣時間0.8 s，這里取衰減帶厚度180 m。主頻為30 Hz，雷克子波為激發震源。建立大小為1800 m×1800 m均勻介質模型(圖1a)、層狀介質模型(圖1b)和復雜介質模型(圖1c)進行解耦與非解耦兩種差分形式的數值模擬，震源坐標分別為(900 m, 900 m)、(900 m, 240 m)、(900 m, 240 m)。

圖1 均勻介質(a)、層狀介質(b)、復雜介質(c)模型

均勻介質模型震源位于模型中心，層狀介質和復雜介質模型震源位于模型上層，為更明確地觀察波場的傳播狀態，均勻介質模型波場快照取0.56 s時刻，層狀模型和復雜模型取0.72 s時刻。圖2、圖3、圖4分別為均勻介質模型、層狀介質模型和復雜介質模型使用時間2階、空間2階或4階有限差分格式，非解耦及解耦PML差分形式進行數值模擬得到的波場快照。

由圖2可以看出，當有限差分階數為2階時，非解耦與解耦PML差分形式數值模擬皆出現了少許的頻散現象(圖2a、b)，使用解耦差分形式時的頻散現象更弱；當空間階數升級為4階時，兩種PML差分形式都得到了較好的數值模擬結果(圖2c、d)。

a.非解耦，時間2階、空間2階；b.解耦，時間2階、空間2階；c.非解耦，時間2階、空間4階；d.解耦，時間2階、空間4階。

從圖3、圖4可以看出，當地質模型變得更加復雜時，使用精度低的差分格式出現的頻散會更加嚴重，解耦算法相對非解耦算法而言數值模擬效果要好一些(圖3a、b，圖4a、b)；空間階數提高，解耦與非解耦算法都可以改善頻散現象(圖3c、d，圖4c、d)。

a.非解耦，時間2階、空間2階；b.解耦，時間2階、空間2階；c.非解耦，時間2階、空間4階；d.解耦，時間2階、空間4階。

a.非解耦，時間2階、空間2階；b.解耦，時間2階、空間2階；c.非解耦，時間2階、空間4階；d.解耦，時間2階、空間4階。

在現實生產試驗中，我們通常使用空間8階或更高空間階數的差分格式進行計算，因此選擇計算和編程更加方便的非解耦PML差分形式進行計算更具優勢。

表1為使用3種模型(圖1a、b、c)時，解耦與非解耦PML差分形式不同差分精度數值模擬的計算時間比較。可以看出，相比解耦PML差分形式，非解耦PML差分形式所用計算時間更少。這是因為解耦PML差分形式比非解耦PML差分形式增加了兩項中間量的差分運算(對Ax、Az的差分計算)。同時根據計算原理可以得知，非解耦PML差分形式在計算上較為簡單，在使用高精度差分格式時，計算和編程都比較容易，計算效率高、易于實現；而解耦PML差分形式在計算和編程上都較為繁瑣。

表1 PML差分形式計算時間比較

2 CE邊界條件與PML邊界條件組合

2005年，吳國忱等[15]提出了一種將吸收類邊界條件和衰減類邊界條件進行組合的組合邊界思想，即在衰減類邊界條件衰減區域的最外層使用吸收類邊界條件對邊界反射進行再次吸收，數值模擬結果顯示這種組合方式有較好的吸收效果及較少的計算量。2009年，杜啟振等[16]運用這一組合思想將改進的擴邊衰減邊界與特征分析法組合，組合邊界減少了衰減邊界中衰減帶的厚度，數值模擬結果顯示這種組合方式有較好的吸收效果。本文采用這種組合思想，使用更具優勢的非解耦PML差分形式，先使用PML衰減邊界衰減到達邊界的入射波，再在衰減邊界外區域設置占用內存更小、方法更簡單的2階CE吸收邊界吸收未衰減完全的入射波。因為使用2階CE邊界條件進行二次吸收，在PML衰減帶內的入射波可以不衰減為0，所以組合邊界條件使用較小的衰減帶厚度。

分別使用均勻介質模型(圖1a)、層狀介質模型(圖1b)以及復雜介質模型(圖1c)驗證所提出組合邊界條件的有效性。采用時間2階、空間12階的高階有限差分進行聲波波場數值模擬。為說明組合邊界的衰減效果，分別給出衰減帶厚度為90 m時，以及增加衰減帶厚度為180 m后的波場模擬快照。圖5為均勻介質模型組合邊界條件和PML邊界條件在0.56 s時刻的波場模擬快照。圖6、圖7分別為層狀介質模型和復雜介質模型組合邊界條件和PML邊界條件在0.72 s時刻的波場模擬快照。可以看出，在使用較小的衰減帶厚度時，組合邊界條件可以較好地吸收邊界反射(圖5a、圖6a、圖7a)。對比圖5a與圖5b、圖6a與圖6b、圖7a與圖7b可以看出，在采用相同衰減帶厚度時，組合邊界條件可以取得較好的吸收效果；對比圖5a與圖5c、圖6a與圖6c、圖7a與圖7c可以看出，組合邊界條件在采用較少的衰減帶厚度時也能達到很好的吸收效果。

a.組合邊界，衰減帶厚度90 m；b.PML邊界，衰減帶厚度90 m；c.PML邊界，衰減帶厚度180 m。

a.組合邊界，衰減帶厚度90 m；b.PML邊界，衰減帶厚度90 m；c.PML邊界，衰減帶厚度180 m。

a.組合邊界，衰減帶厚度90 m；b.PML邊界，衰減帶厚度90 m；c.PML邊界，衰減帶厚度180 m。

組合邊界條件和PML邊界條件吸收邊界反射所用時間如表2所示。從表2中可以看出：3種模型組合邊界所用時間與同樣衰減帶厚度下PML邊界所用時間相差不大，這表明組合邊界中所使用的2階CE邊界條件所占內存很小；而組合邊界所用時間小于增加衰減帶厚度后PML邊界所用時間，說明組合邊界條件的應用在保證吸收效果的同時可以有效提高計算效率，驗證了這種組合邊界的有效性。

表2 計算效率比較

許多專家學者如裴正林[17-19]、夏凡[20]、王均[21]等均已分別應用CE邊界條件和PML邊界條件在三維地質模型中進行實驗并取得較好的波場模擬效果，本文所提出的組合方法是CE邊界條件和PML邊界條件的有機結合，因此也可推廣應用在三維的波場數值模擬中。

3 結論

1)不同地質模型地震聲波有限差分數值模擬結果表明：差分階數較低時，解耦PML差分形式對頻散的壓制效果優于非解耦PML差分形式；隨著差分階數的提高，頻散得到有效壓制，非解耦差分形式在計算效率和實現方式上都更有優勢。

2)本文給出了一種2階CE邊界條件與PML邊界條件的組合方式，并通過實際算例驗證了組合邊界的有效性。此邊界條件的組合方式在保證吸收效果的同時可以通過減少衰減帶厚度有效提高計算效率，可以作為一種新的人工邊界條件應用在聲波或彈性波波動方程二維或三維數值模擬中。