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

韓復(fù)興，王若雯，孫章慶，高正輝

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

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

近年來(lái)許多相關(guān)專(zhuān)家學(xué)者都對(duì)CE邊界條件和PML邊界條件的應(yīng)用效果進(jìn)行了比較，包括定性比較和定量比較。定性比較主要通過(guò)觀察波場(chǎng)快照或合成地震記錄的邊界反射細(xì)節(jié)來(lái)比較邊界條件的吸收效果。如：王守東[2]導(dǎo)出了完全匹配層法的聲波方程表達(dá)形式，并通過(guò)增強(qiáng)地震記錄振幅的顯示方式得出PML邊界條件吸收效果優(yōu)于CE邊界條件的結(jié)論；邢麗[3]通過(guò)數(shù)值算例分析認(rèn)為，CE邊界條件在程序方面易于實(shí)現(xiàn)，但1階CE邊界條件精度低，2階CE邊界條件有時(shí)又會(huì)出現(xiàn)不穩(wěn)定的現(xiàn)象，PML邊界條件雖然程序?qū)崿F(xiàn)較為繁瑣，但計(jì)算過(guò)程穩(wěn)定，選取適當(dāng)?shù)膮?shù)時(shí)，幾乎可以吸收全部的邊界反射；王永剛等[4]分別采用井間和地面地質(zhì)模型進(jìn)行數(shù)值模擬，結(jié)果表明PML邊界條件相對(duì)于旁軸近似法具有優(yōu)越性；裴俊勇等[5]利用CE邊界條件和PML邊界條件進(jìn)行有限差分正演模擬，模擬圖像觀察結(jié)果說(shuō)明PML邊界條件吸收邊界效果較好。

定量比較通過(guò)量化邊界反射數(shù)值進(jìn)行，主要方法為根據(jù)反射波峰值高低比較邊界條件吸收效果強(qiáng)弱。如，付小波等[6]對(duì)CE邊界條件和PML邊界條件進(jìn)行了精細(xì)的定量比較性研究，發(fā)現(xiàn)CE邊界條件相對(duì)PML邊界條件受入射波角度限制更大，當(dāng)入射角增大時(shí)，CE邊界條件吸收效果大大降低，PML邊界條件則不受入射波角度的限制。

本文通過(guò)對(duì)PML邊界條件和CE邊界條件進(jìn)行詳細(xì)的分析研究，考慮到其各自的優(yōu)缺點(diǎn)，將CE邊界條件和PML邊界條件組成新的組合邊界，并以數(shù)值模擬結(jié)果驗(yàn)證算法的有效性，以期達(dá)到在減少衰減帶厚度的同時(shí)提高計(jì)算效率的目的。

1 PML邊界條件兩種差分形式對(duì)比

PML衰減介質(zhì)中的波動(dòng)方程可以看作是常規(guī)波動(dòng)方程的推廣，入射波傳播到介質(zhì)中振幅呈指數(shù)衰減。當(dāng)衰減因子和衰減帶厚度選擇適當(dāng)時(shí)，理論上入射波不會(huì)再反射回計(jì)算區(qū)域[7]。

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

本文將討論在地震聲波有限差分?jǐn)?shù)值模擬過(guò)程中實(shí)現(xiàn)全局SPML的非解耦與解耦差分形式的基本原理和實(shí)現(xiàn)效果。

1.1 全局SPML非解耦差分形式

均勻介質(zhì)中二維聲波波動(dòng)方程為

(1)

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

(2)

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

(3)

1.2 全局SPML解耦差分形式

(4)

(5)

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

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

(6)

(7)

(8)

(9)

(10)

1.3 解耦與非解耦PML差分形式數(shù)值模擬對(duì)比分析

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

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

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

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

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

從圖3、圖4可以看出，當(dāng)?shù)刭|(zhì)模型變得更加復(fù)雜時(shí)，使用精度低的差分格式出現(xiàn)的頻散會(huì)更加嚴(yán)重，解耦算法相對(duì)非解耦算法而言數(shù)值模擬效果要好一些(圖3a、b，圖4a、b)；空間階數(shù)提高，解耦與非解耦算法都可以改善頻散現(xiàn)象(圖3c、d，圖4c、d)。

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

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

在現(xiàn)實(shí)生產(chǎn)試驗(yàn)中，我們通常使用空間8階或更高空間階數(shù)的差分格式進(jìn)行計(jì)算，因此選擇計(jì)算和編程更加方便的非解耦PML差分形式進(jìn)行計(jì)算更具優(yōu)勢(shì)。

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

表1 PML差分形式計(jì)算時(shí)間比較

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

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

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

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邊界條件吸收邊界反射所用時(shí)間如表2所示。從表2中可以看出：3種模型組合邊界所用時(shí)間與同樣衰減帶厚度下PML邊界所用時(shí)間相差不大，這表明組合邊界中所使用的2階CE邊界條件所占內(nèi)存很??；而組合邊界所用時(shí)間小于增加衰減帶厚度后PML邊界所用時(shí)間，說(shuō)明組合邊界條件的應(yīng)用在保證吸收效果的同時(shí)可以有效提高計(jì)算效率，驗(yàn)證了這種組合邊界的有效性。

表2 計(jì)算效率比較

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

3 結(jié)論

1)不同地質(zhì)模型地震聲波有限差分?jǐn)?shù)值模擬結(jié)果表明：差分階數(shù)較低時(shí)，解耦PML差分形式對(duì)頻散的壓制效果優(yōu)于非解耦PML差分形式；隨著差分階數(shù)的提高，頻散得到有效壓制，非解耦差分形式在計(jì)算效率和實(shí)現(xiàn)方式上都更有優(yōu)勢(shì)。

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