賀同江,劉紅艷,李小凡

(1.天津市地震局,天津 300201;2.中國科學(xué)院 地質(zhì)與地球物理研究所,北京 100029)

褶積微分算子法數(shù)值模擬的邊界問題研究

賀同江1,劉紅艷1,李小凡2

(1.天津市地震局,天津 300201;2.中國科學(xué)院 地質(zhì)與地球物理研究所,北京 100029)

褶積微分算子法是一種全新的數(shù)值模擬方法,已被廣泛應(yīng)用于復(fù)雜介質(zhì)的地震波場數(shù)值模擬,但是其邊界反射問題一直沒有解決。這里將最佳匹配層(PML)吸收邊界條件引入到褶積微分算子法中,此方法是在研究區(qū)域的邊界上加入吸收層,使邊界上傳入吸收層的波,隨傳播距離按指數(shù)規(guī)律衰減,不產(chǎn)生任何反射,以達(dá)到消除邊界反射的目的。構(gòu)造不同的模型,通過對比分析證明,PML吸收邊界條件能比較好地解決褶積微分算子法的邊界問題,從而驗(yàn)證了完全匹配層吸收效果的優(yōu)越性。

褶積微分算子;數(shù)值模擬;PML吸收邊界

邊界吸收的方法很多,可謂舉不勝舉,但是真正能有效的用于作者在本文所構(gòu)造的褶積微分算子法數(shù)值模擬的卻少之又少。通過模型試算分析可知,C layton吸收邊界條件適用于有限差分?jǐn)?shù)值模擬,但是對作者所構(gòu)造的廣義正交多項(xiàng)式褶積微分算子法的數(shù)值模擬的邊界反射無效。Cerjan吸收邊界條件適用于偽譜法數(shù)值模擬,理論上也能適用于作者在本文所構(gòu)造的基于Fo rsyte廣義正交多項(xiàng)式褶積微分算子法的數(shù)值模擬,但是采用Cerjan吸收邊界條件并不能完全適用本文方法。因此作者在本文中采用PML吸收邊界條件,通過模擬實(shí)例,證明了PML吸收邊界能比較好地適用于褶積微分算子法。

1 Fo rsyte廣義正交多項(xiàng)式褶積微分算子求導(dǎo)原理

作者在本文中利用錯格有限差分褶積微分算子法來計算彈性介質(zhì)中的地震波場,其主導(dǎo)思想是:利用基于Forsyte廣義正交多項(xiàng)式褶積微分算子,有效表示計算波場對空間的偏導(dǎo)數(shù),采用錯格有限差分法計算對時間的偏導(dǎo)數(shù),其計算思路類似于偽譜法[1]。

首先給出Fo rsyte多項(xiàng)式微分算子。Forsyte多項(xiàng)式是一個廣義正交多項(xiàng)式,其插值函數(shù)為:

其中 P0=1

f(xi)為被插值函數(shù)f(x)在點(diǎn)xi處的值。在式(1)中的P0(x)、…、Pj+1(x)定義為Forsyte多項(xiàng)式系數(shù)。

對式(1)中的x求導(dǎo),可得:

其中

Forsyte多項(xiàng)式微分算子可寫為式(3)。

將上述微分算子式(3)離散化,可得式(4)。

其中 i為采樣指標(biāo);Δx為沿著x軸的采樣間隔。

就實(shí)際應(yīng)用而言,須將微分算子截成短算子,這勢必引起Gibbs現(xiàn)象。另外,多項(xiàng)式的引入還將引起Runge現(xiàn)象,為了消除這些現(xiàn)象,必須采用窗函數(shù),以截斷長微分算子。作者在本文采用的是下列Gaussian窗函數(shù):

其中 m x為單邊截斷長度的采樣數(shù);c為常數(shù);a(0.1≤a≤0.75)為衰減因子。

將微分算子式(4)用式(5)截斷并鋸齒化后,可得如下實(shí)用的一階褶積微分算子:

通過對算子長度的調(diào)節(jié)及算子系數(shù)的優(yōu)化,可同時兼顧波場解的全局信息與局部信息。作者運(yùn)用算子長度為9點(diǎn)的一階褶積微分算子,求解波動方程。而9點(diǎn)褶積算子的最優(yōu)權(quán)系數(shù)為:

可以看出,算子的系數(shù)是反對稱的。

2 二維各向同性介質(zhì)一階速度(應(yīng)力彈性波動方程及PML吸收邊界)

利用完全匹配層作為吸收邊界的基本做法,是在所研究區(qū)域的四周引入完全匹配層。如圖1所示,區(qū)域ABCD為要研究的區(qū)域,即我們要在此區(qū)域中研究波的傳播問題。在區(qū)域的周圍加上完全匹配層,在區(qū)域“1”中,令d(x)≠0;d(z)≠0,速度V都等于角點(diǎn)的速度。在區(qū)域“2”中,令d(x)=0;d(z)≠0,速度V在z方向?yàn)槌?shù),在x方向和邊界的速度相等。在區(qū)域“3”中,令d(x)≠0;d(z)=0,速度V在x方向?yàn)槌?shù),在z方向和邊界的速度相等。這樣在計算邊界的周圍都有完全匹配層吸收介質(zhì),波由區(qū)域內(nèi)通過邊界傳播到完全匹配層時,就不會產(chǎn)生任何反射。波在完全匹配層中傳播時,也不會產(chǎn)生反射,并且按傳播距離的指數(shù)規(guī)律衰減。當(dāng)波傳播到完全匹配層的邊界時,波場近似為零,也不會產(chǎn)生反射。

圖1 PML吸收邊界示意圖Fig.1 Schem atic diagram of the PML abso rbing boundary

對于二維非均勻各向同性介質(zhì),一階速度,即應(yīng)力彈性波動方程(假定體力為零)為式(7):

式中 vx=?ux/?t,vz=?uz/?t為速度的x、z分量;σxx、σzz為正應(yīng)力;σxz為剪應(yīng)力;λ、μ為拉梅常數(shù);λ=ρ(V2p-2V2S);μ=ρVS2;VP、VS分別為介質(zhì)縱波與橫波速度;ρ為介質(zhì)密度。

將一階應(yīng)力,即速度波動方程的波場分為與x方向和與z方向有關(guān)的二個子分量:

在式(8)中,上標(biāo)x和z代表該項(xiàng)只與相應(yīng)的空間導(dǎo)數(shù)有關(guān)。

根據(jù)二維非均勻各向同性介質(zhì)波動方程組,可以分裂得到二維非均勻各向同性介質(zhì)波動方程帶有衰減因子的PML吸收邊界系統(tǒng)方程組:

其中

式中 dPML是PML吸收層的厚度,作者取二十個網(wǎng)格點(diǎn)數(shù);x、z是PML內(nèi)計算點(diǎn)到內(nèi)部區(qū)域和PML層交界面的距離;R為理論反射系數(shù)(作者在數(shù)值模擬中取為10-6),式(10)的解是衰減的。

d(x)和d(z)分別為x方向和z方向的衰減系數(shù),也就是說d(x)起到衰減x方向傳播的波,d(z)起到衰減z方向傳播的波的作用。對于任意方向傳播的波,可以通過矢量分解,分解成x方向和z方向傳播的波,分別進(jìn)行衰減。波場是按傳播距離的指數(shù)規(guī)律衰減,衰減速度很快。當(dāng)衰減系數(shù)d(x)、d(z)隨空間位置變化時,不會在介質(zhì)中產(chǎn)生任何反射。

3 二維一階速度(應(yīng)力彈性波動方程PML吸收邊界條件下的數(shù)值解格式)

對模型區(qū)間離散后,設(shè)n、m、k分別是沿著空間x軸、z軸、時間t軸的采樣點(diǎn)數(shù),Δx、Δz、Δt分別是沿x軸、z軸、t軸的采樣間隔,m x、m z是沿x軸、z軸采樣數(shù)的半算子長度,則式(10)二維非均勻各向同性介質(zhì)完全匹配層彈性波一階速度——應(yīng)力方程,其離散化的時間錯格差分褶積微分算子法格式為(體力為零)式(11)。

4 數(shù)值算例

模型網(wǎng)格點(diǎn)數(shù)為256×256,震源位于模型正中,震源為45°傾斜集中力源,震源R icker子波主頻為20 Hz。采用不同的吸收邊界條件進(jìn)行對比,模型參數(shù)如表1所示。

(1)圖2(見下頁)不加邊界條件時,邊界反射很強(qiáng),并且在邊界條件處,由于邊界反射引起很嚴(yán)重的數(shù)值頻散現(xiàn)象,使模擬精度降低。

(2)圖3(見下頁)是加入Cerjon吸收邊界條件的模擬結(jié)果,相對圖2邊界反射引起的數(shù)值頻散已幾乎消失,但是邊界反射還很嚴(yán)重。從圖3可以看出,波遇到邊界產(chǎn)生的反射波。

(3)圖4(見下頁)是加入PML吸收邊界條件的模擬結(jié)果,看不到任何邊界反射。從模擬結(jié)果看,邊界條件是數(shù)值模擬中需解決的關(guān)鍵問題,作者在本文中所構(gòu)造的方法,Cerjon吸收邊界條件不能完全適用,但是PML吸收邊界條件卻能取得較好的結(jié)果,適于本文中所構(gòu)造的方法。

我們利用一個比較復(fù)雜的非均勻介質(zhì)模型,對本文方法加PML邊界和Cerjon吸收邊界作比較。復(fù)雜介質(zhì)模型結(jié)構(gòu)如后面圖5所示。模型參數(shù)如后面表2所示。模型的網(wǎng)格大小為256×256,網(wǎng)格間距為d x=d z=10m,時間步長為1m s,震源放置于介質(zhì)①中,坐標(biāo)為(128,50),震源為45°傾斜集中力源,震源R icker子波主頻為20 Hz。

表1 均勻各向同性介質(zhì)模型參數(shù)Tab.1 M odelparam etersof homogeneous and isotrop icm edia

圖2 均勻各向同性介質(zhì)模型,未加入吸收邊界條件400m s時刻波場快照Fig.2 W ave-field snapshots in homogeneous and isotrop icmodelw ith non-absorbing boundary(t=400m s,left:x component,right:z component)

圖3 均勻各向同性介質(zhì)模型,加入Cerjon吸收邊界條件400m s時刻波場快照Fig.3 W ave-field snapsho ts in homogeneous and iso trop icmodelw ith Cerjon abso rbing boundary(t=400m s,left:x component,right:z component)

圖4 均勻各向同性介質(zhì)模型,加入PML吸收邊界條件400m s時刻波場快照Fig.4 W ave-field snap shots in homogeneous and isotrop icmodelw ith PML absorbing boundary(t=400m s,left:x component,right:z component)

表2 復(fù)雜非均勻介質(zhì)模型參數(shù)Tab.2 M odelparam etersof comp lex and inhomogeneousm edia

圖5 復(fù)雜非均勻介質(zhì)模型Fig.5 Com p lex and inhomogeneousmodel

從圖6、圖7(見下頁)可以看出,作者在本文所構(gòu)造的基于Forsyte廣義正交多項(xiàng)式褶積微分算子法,能較好地模擬特別復(fù)雜的非均勻介質(zhì)中的地震波場。采用Cerjon吸收邊界條件,不能有效地消除邊界反射,邊界反射比較明顯,影響了模擬的精度,因而不能完全適用于作者在本文所構(gòu)造的,基于Forsyte廣義正交多項(xiàng)式褶積微分算子法的波場數(shù)值模擬中。而采用PML吸收邊界,不會產(chǎn)生任何反射,可達(dá)到消除邊界反射的目的。因此采用PML吸收邊界條件,能夠完全適應(yīng)本文方法。

5 結(jié)論

(1)采用完全匹配層(PML)吸收邊界條件,是消除邊界反射的理想方法之一。在研究區(qū)域的四周加入很薄的匹配層,以很小的計算代價,可換取很好的吸收效果。

(2)在研究區(qū)域的邊界上加入吸收層,使邊界上傳入吸收層的波,隨傳播距離按指數(shù)規(guī)律衰減,不產(chǎn)生任何反射,可達(dá)到消除邊界反射的目的。

圖6 加入Cerjon吸收邊界條件t=500m s時波場快照Fig.6 W ave-field snapshotsw ith Cerjon absorbing boundary of comp lexmodel(t=500m s,left:x component,right:z component)

(3)基于Forsyte廣義正交多項(xiàng)式褶積微分算子法,能較好地模擬特別復(fù)雜的非均勻介質(zhì)中的地震波場,PML吸收邊界解決了本文方法的邊界問題。作者在本文為褶積微分算子法的后續(xù)研究做了鋪墊工作。可以預(yù)期,褶積微分算子法的推出及后續(xù)研究的成功開展,將為高精度地震波模擬,地震波偏移,地震反演,地震波成像,以及地震波在復(fù)雜非均勻介質(zhì)中的傳播等研究,提供更為廣泛的選擇。

[1]程冰潔,李小凡,龍桂華.基于廣義正交多項(xiàng)式褶積微分算子的地震波場數(shù)值模擬方法[J].地球物理學(xué)報,2008,51(2):531.

[2]程冰潔,李小凡.2.5維地震波場褶積微分算子法數(shù)值模擬[J].地球物理學(xué)進(jìn)展,2008,23(4):1099.

[3]劉紅艷,李小凡.非均勻介質(zhì)中地震波傳播的數(shù)值模擬[J].物探化探計算技術(shù),2008,30(3):173.

[4]李信富,李小凡.復(fù)雜介質(zhì)地震波傳播的褶積微分算子數(shù)值模擬[J].地震學(xué)報,2008,30(4):377.

[5]CLAYTON R,ENGQU ISTB.Absorbing boundary conditions for acoustic and elastic wave equations[J].Bu ll.Seis.Soc.Am.,1997,67:1529.

[6]CERJAN C,KOSLOFFD.A non-reflection boundary condition for discrete acoustic and elasticw ave equation[J].Geophysics,1985,50:705.

[7]BERENGER JP.Three-dim ensionalperfectlym atched layer for the absorp tion electrom agnetic waves[J].Comp.Phys,1996,127:363.

[8]CHW EW C,L IU Q H.Perfectly m atched layers for elastodynam ics:A new absorbing boundary condition[J].J.Comp.A coust.,1996,4(4):72.

[9]FRANC ISCOLL INO,CHRYSOULA TSOGKA.App lication of the perfectlym atched absorbing lamodel to the linear elastodynam ic p rob lem in anisotrop ic heterogeneousm edia[J].Geophysics,2001,66(1):294.

[10]FESTA G,N IELSEN S.PML Absorbing Boundaries[J].Bulletin of the Seismo logical Society Am erica,2003,93(2):891.

[11]BECACHE E,FAUQUEUX S,JOLY P.Stability of perfectly m atched layers,group velocities and anisotrop ic waves[J].Journal of Com pu tational Physics,2003,188:399.

[12]ZENG Y Q,HE JQ,L IU Q H.The app lication of the perfectlym atched layer in num ericalmodeling ofwave p ropagation in po roelastic m edia[J].Geophysics,2001,66(4):1258.

[13]SONGRUOLONG,MA JUN,WANG KEX IE.The App lication of the nonsp litting perfectlym atched layer in num ericalmodeling of wave p ropagation in poroelastic m edia[J].App lied Geophysics,2005,2(4):216.

[14]SHUO MA,PENGCHENG L IU.M odeling of the PerfectlyM atched Layer absorbing boundary and intrinsic attenuation in exp licit Finite-Elem entm ethods[J].Bu lletin of the Seismo logical Society of Am erica,2006,96(5):1779.

[15]趙海波,王秀明,王東,等.完全匹配層吸收邊界在孔隙介質(zhì)彈性波模擬中的應(yīng)用[J].地球物理學(xué)報,2007,50(2):581.

圖7 加入PML吸收邊界t=500m s時波場快照Fig.7 W ave-field snapshotsw ith PML abso rbing boundary of comp lexmodel(t=500m s,left:x component,right:z component)

P 631.4

A

1001—1749(2010)06—0571—07

0 前言

國家自然科學(xué)基金重點(diǎn)項(xiàng)目資助(40437018)

2009-07-28

賀同江(1979-),男,碩士,現(xiàn)在天津市地震局監(jiān)測預(yù)報中心工作。

褶積微分算子法作為一種全新的數(shù)值模擬方法,已被廣泛應(yīng)用于復(fù)雜介質(zhì)的地震波場數(shù)值模擬中,該方法同時具有廣義正交多項(xiàng)式方法的高精度和短算子有限差分算法的高速度。通過對算子長度的調(diào)節(jié)及算子系數(shù)的優(yōu)化,可同時兼顧波場解的全局信息與局部信息。該算法的計算速度快,計算效率高,能夠直觀、高效地反映介質(zhì)中波場的傳播規(guī)律[1～4]。但是在進(jìn)行波動方程數(shù)值模擬時,考慮到計算機(jī)的有限內(nèi)存及有限計算時間,要對考慮問題的無限區(qū)域進(jìn)行截取,使進(jìn)行中的數(shù)值模擬在有限的區(qū)域內(nèi)完成,即可以在有限的區(qū)域中得到較高精度的無限空間解。為此需要引入人工邊界來達(dá)到此目的,但這樣會在人工邊界處產(chǎn)生人為反射,如不消除或者壓制這種虛假反射,就會影響數(shù)值模擬的結(jié)果和精度。為此,邊界條件問題是必須要解決的,也是較難解決的問題。地震波數(shù)值模擬中,在人工邊界上使用吸收邊界條件,應(yīng)盡可能地降低邊界上所產(chǎn)生的反射波強(qiáng)度,使得反射波在邊界上好像被“吸收”了一樣,從而大大減小對計算區(qū)域內(nèi)精度的影響,以達(dá)到邊界吸收的目的。應(yīng)用較普遍的是C layton、Engquist等人[5]提出的吸收邊界條件(即CE邊界),該條件在旁軸近似理論的基礎(chǔ)上導(dǎo)出,在特定的入射角和頻率范圍內(nèi),具有較好的吸收效果。Cerjan etal.[6]引入了有損吸收層來衰減外行波,但這種方法需要足夠厚的吸收層,才可以得到滿意的效果,因此計算量較大。Berenger[7]對海綿吸收邊界條件作了進(jìn)一步的改善,提出一種新的人工邊界條件,即最佳匹配層法(PML吸收邊界),這種方法幾乎達(dá)到零反射。目前,PML技術(shù)已廣泛應(yīng)用于聲波和彈性波的有關(guān)問題。Chew和L iu.[8]首先證明了PML技術(shù)可應(yīng)用于彈性波的數(shù)值模擬中;Co llino和Tsogka[9]把PML用在各向異性非均勻的介質(zhì)情況;Festa和N ielson[10]針對各向同性彈性介質(zhì)PML情況下,離散波動方程穩(wěn)定性進(jìn)行了詳細(xì)討論,并推導(dǎo)了PML理論反射系數(shù)計算公式;Becache et al.[11]在理論上討論了PML在各向異性情況下穩(wěn)定性問題。另外PML也應(yīng)用到其它問題上,如黏彈性介質(zhì)、雙相孔隙彈性介質(zhì)中[12～15]。