賀同江,劉紅艷,李小凡

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

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

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

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

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

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

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

1 Fo rsyte廣義正交多項式褶積微分算子求導原理

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

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

其中 P0=1

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

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

其中

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

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

其中 i為采樣指標;Δx為沿著x軸的采樣間隔。

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

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

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

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

可以看出,算子的系數是反對稱的。

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

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

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

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

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

將一階應力,即速度波動方程的波場分為與x方向和與z方向有關的二個子分量:

在式(8)中,上標x和z代表該項只與相應的空間導數有關。

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

其中

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

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

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

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

4 數值算例

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

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

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

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

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

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

圖2 均勻各向同性介質模型,未加入吸收邊界條件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 均勻各向同性介質模型,加入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 均勻各向同性介質模型,加入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 復雜非均勻介質模型參數Tab.2 M odelparam etersof comp lex and inhomogeneousm edia

圖5 復雜非均勻介質模型Fig.5 Com p lex and inhomogeneousmodel

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

5 結論

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

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

圖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廣義正交多項式褶積微分算子法,能較好地模擬特別復雜的非均勻介質中的地震波場,PML吸收邊界解決了本文方法的邊界問題。作者在本文為褶積微分算子法的后續(xù)研究做了鋪墊工作?？梢灶A期,褶積微分算子法的推出及后續(xù)研究的成功開展,將為高精度地震波模擬,地震波偏移,地震反演,地震波成像,以及地震波在復雜非均勻介質中的傳播等研究,提供更為廣泛的選擇。

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

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

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

[4]李信富,李小凡.復雜介質地震波傳播的褶積微分算子數值模擬[J].地震學報,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]趙海波,王秀明,王東,等.完全匹配層吸收邊界在孔隙介質彈性波模擬中的應用[J].地球物理學報,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 前言

國家自然科學基金重點項目資助(40437018)

2009-07-28

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

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