一種適用于任意高階間斷有限元的高精度非分裂完全匹配層吸收邊界方法*

何洋洋 翁 斌 張金淼

(1.中海油研究總院 北京 100028； 2.中國(guó)石油大學(xué)(北京) 北京 102249)

一種適用于任意高階間斷有限元的高精度非分裂完全匹配層吸收邊界方法*

何洋洋1,2翁 斌1張金淼1

(1.中海油研究總院 北京 100028； 2.中國(guó)石油大學(xué)(北京) 北京 102249)

任意高階間斷有限元法常用的吸收邊界條件存在對(duì)切向入射波吸收效果不佳等問(wèn)題。本文提出了一種適用于任意高階間斷有限元的高精度非分裂完全匹配層吸收邊界方法，將非分裂完全匹配層吸收邊界應(yīng)用于任意高階間斷有限元地震波數(shù)值模擬方法中，建立了完全匹配層內(nèi)新的波動(dòng)方程，推導(dǎo)了其求解過(guò)程及在三角形單元內(nèi)的表達(dá)形式，最后給出了離散化格式。該方法在完全匹配層也可以求得任意高階時(shí)間精度和空間精度的數(shù)值解，與計(jì)算區(qū)域的精度一致，減少了邊界反射波的能量。數(shù)值算例模擬結(jié)果表明，本文方法對(duì)不同角度的入射波具有較好的吸收效果，適用于高泊松比介質(zhì)。

非分裂完全匹配層；任意高階間斷有限元；高精度；不同角度入射波；高泊松比

地震波場(chǎng)數(shù)值模擬是勘探地球物理的重要研究手段，實(shí)際應(yīng)用中，可以通過(guò)求解彈性波動(dòng)方程來(lái)實(shí)現(xiàn)[1-2]。針對(duì)具有復(fù)雜幾何結(jié)構(gòu)介質(zhì)中的地震波傳播進(jìn)行高精度模擬是一項(xiàng)具有挑戰(zhàn)性的工作，近年來(lái)已提出許多新的方法，比如有限體積法[3]、混合有限元法[4]、邊界元法[5]等。K?ser和Dumbser[6-7]提出了彈性介質(zhì)中的任意高階間斷有限元地震波數(shù)值模擬方法，該方法將間斷有限元法[8]和任意高階時(shí)間積分法[9]相結(jié)合，從而得到空間和時(shí)間上任意高階精度的數(shù)值解，能夠很好地處理具有局部小尺度結(jié)構(gòu)的復(fù)雜介質(zhì)模型，應(yīng)用范圍廣泛。實(shí)際應(yīng)用中，由于地震波數(shù)值模擬每次只能計(jì)算有限區(qū)域中的波傳播，所以必須在計(jì)算區(qū)域的人工邊界處使用吸收邊界才能模擬無(wú)限大的地下介質(zhì)。K?ser和Dumbser[6]給出了一種簡(jiǎn)單的吸收邊界，可以很好地吸收法向方向入射波，是任意高階間斷有限元法常用的吸收邊界處理方法，但在某些情況下，比如切向入射波，或者模型的角點(diǎn)處，這種吸收邊界不能有效地對(duì)入射波進(jìn)行吸收。

Bérenger[10]在電磁波場(chǎng)模擬中提出了完全匹配層吸收邊界條件(PML)，這是目前已知的效果最好的吸收邊界條件。Chew和Liu[11]以及Hastings等[12]將完全匹配層用于彈性波模擬，取得了好的吸收效果。目前完全匹配層普遍使用波場(chǎng)分裂的方法，以避免計(jì)算中出現(xiàn)時(shí)間域卷積，但是當(dāng)?shù)卣鸩ㄇ邢蛉肷鋾r(shí)該方法的邊界反射系數(shù)會(huì)變大，產(chǎn)生強(qiáng)能量的反射波，對(duì)波場(chǎng)模擬結(jié)果造成干擾。

Komatitsch和Martin[13]，Drossaert和Giannopoulos[14]等分別提出了一種非分裂形式的完全匹配層算法(CPML)，采用遞歸卷積技術(shù)[15]，在提高吸收效果的同時(shí)并沒(méi)有增加運(yùn)算量。該方法在完全匹配層內(nèi)添加了一個(gè)濾波器，地震波切向入射時(shí)會(huì)在匹配層內(nèi)傳播很長(zhǎng)的距離，濾波器會(huì)對(duì)這種波產(chǎn)生強(qiáng)的吸收衰減，吸收效果優(yōu)于傳統(tǒng)的完全匹配層算法。

本文將非分裂完全匹配層吸收邊界應(yīng)用于任意高階間斷有限元地震波數(shù)值模擬方法中，建立了完全匹配層內(nèi)新的波動(dòng)方程，推導(dǎo)了其求解過(guò)程及在三角形單元內(nèi)的表達(dá)形式，最后給出了離散化格式。通過(guò)數(shù)值算例分析，研究了該算法對(duì)不同角度入射波的吸收效果以及該算法在高泊松比介質(zhì)中波場(chǎng)模擬的適用性。

1 高精度非分裂完全匹配層算法

1.1 完全匹配層內(nèi)波動(dòng)方程的建立

二維各向同性彈性介質(zhì)中的波動(dòng)方程為

(1)

(2)

(3)

式(1)～(3)使用張量表示法。其中，u=(σxx,σzz,σxz,v,w)T為5個(gè)未知量組成的向量，up、uq分別為向量中的第p、q個(gè)分量，p=1，…，5，q=1，…，5為自由指標(biāo)；時(shí)間t∈R，x=(x,z)T∈R2；σxx=σxx(x,t)、σzz=σzz(x,t)是法向應(yīng)力分量，σxz=σxz(x,t)是切向應(yīng)力分量；v=v(x,t)和w=w(x,t)是質(zhì)點(diǎn)速度x方向和z方向分量。Apq=Apq(x)和Bpq=Bpq(x)是大小為5×5的Jacobian矩陣，λ=λ(x)和μ=μ(x)是拉梅常數(shù)，ρ=ρ(x)為密度。Sp(x,t)=(S1(x,t),S2(x,t),S3(x,t),S4(x,t),S5(x,t))T是震源向量。

在時(shí)間域內(nèi)，非分裂完全匹配層內(nèi)的波動(dòng)方程可以表示為[13-14]

(4)

(5)

(6)

(7)

式(4)～(7)中：dj(x)為衰減系數(shù)；aj(x)和κj(x)用來(lái)控制能量衰減的速度，地震波模擬中可取κx(x)=κz(x)=1[14]。H(t)為階躍函數(shù)。L為完全匹配層的厚度，可選為1/3或1/2波長(zhǎng)，當(dāng)吸收效果不能滿足要求時(shí)，則適當(dāng)增加完全匹配層的厚度。衰減系數(shù)在完全匹配層的每個(gè)三角形單元內(nèi)是隨空間位置而變化的。對(duì)于一般情況下的波場(chǎng)正演，可令dmax=-3vPlgRc/(2L),αmax=πf0，其中vP是介質(zhì)的縱波速度，f0是震源子波的主頻，Rc為數(shù)值離散后的反射系數(shù)，通常取0.1%[13-14]。

本文對(duì)式(4)進(jìn)行改寫，得到

(8)

(9)

式(9)中：S1、S2、S3、S4、S5是式(1)中的震源項(xiàng)，因?yàn)橥耆ヅ鋵觾?nèi)沒(méi)有震源，所以此處S1=S2=S3=S4=S5=0。

1.2 完全匹配層內(nèi)波動(dòng)方程的求解

(10)

(11)

(12)

(13)

對(duì)于計(jì)算區(qū)域，K?ser和Dumbser[6]在使用任意高階間斷有限元法求解式(1)而得到up的k階時(shí)間導(dǎo)數(shù)時(shí)使用了式(14)，即

-1)k(Ap q+Bp q)kuq(x,t)+

(14)

式(14)中，Sq表示震源向量中的第q個(gè)分量，其k階時(shí)間導(dǎo)數(shù)可以由震源的時(shí)間域表達(dá)式直接計(jì)算得到。在求解t時(shí)刻方程時(shí)，可先求得uq的數(shù)值解，然后通過(guò)式(14)求得up的第k階時(shí)間導(dǎo)數(shù)，從而可以得到k階精度的數(shù)值解。

圖1 完全匹配層內(nèi)up第k階時(shí)間導(dǎo)數(shù)的計(jì)算步驟

1.3 三角形單元內(nèi)的數(shù)值解

(15)

(16)

式(15)～(16)使用張量表示法。其中，l=1，…，N，為自由指標(biāo)；N為多項(xiàng)式基函數(shù)的個(gè)數(shù)，取決于多項(xiàng)式基函數(shù)的階數(shù)；Φl(x)的階數(shù)決定了任意高階間斷有限元法的計(jì)算精度，階數(shù)越高，精度越高。

1.4 離散化

對(duì)于式(13)和式(14)，兩邊同時(shí)乘以試驗(yàn)函數(shù)Φn(x)，并在x-z坐標(biāo)系下的三角形單元T(i)內(nèi)進(jìn)行積分，然后將積分式映射到ξ-η坐標(biāo)系下的標(biāo)準(zhǔn)三角形單元TE中?？紤]到多項(xiàng)式基函數(shù)為正交多項(xiàng)式，式(13)和式(14)在第n個(gè)時(shí)間步的離散計(jì)算格式分別為

]n=

(17)

(18)

式(17)～(18)中：l=1,…，N，m=1，…，N，為不同的自由指標(biāo)；〈·〉表示在標(biāo)準(zhǔn)三角形單元TE中的積分，積分項(xiàng)O、P都可以在TE中提前計(jì)算，這樣可以減少正演模擬的計(jì)算量。

2 數(shù)值算例分析

2.1 均勻彈性介質(zhì)

均勻彈性介質(zhì)模型(圖2)大小為600 m×600 m的正方形，介質(zhì)中P波速度為2 500 m/s，S波速度為1 225 m/s，泊松比為0.34，密度為2 000 kg/m3。本模型所用震源有2種:第1種是脹縮源，只產(chǎn)生P波;第2種是剪切源，只產(chǎn)生S波。震源位于(150 m,150 m)處，震源子波為主頻f0=20 Hz的Ricker子波，子波時(shí)延為t0=0.05 s。波場(chǎng)模擬總時(shí)長(zhǎng)為T=0.75 s，時(shí)間步長(zhǎng)為Δt=0.000 3 s。

圖2 均勻彈性介質(zhì)模型

分別使用高精度非分裂完全匹配層算法及常用的吸收邊界處理方法[6]對(duì)均勻彈性介質(zhì)模型進(jìn)行波場(chǎng)模擬，為了衡量吸收邊界對(duì)入射波的吸收效果，再按照Drossaert和Giannopoulos的方法[14]計(jì)算波場(chǎng)的相對(duì)誤差(誤差越小，說(shuō)明吸收邊界的效果越好)，其結(jié)果見圖3?？梢钥闯?，無(wú)論是P波還是S波，本文方法對(duì)不同角度入射波的吸收效果都優(yōu)于常用方法，地震波在x=150 m處為法向入射，在x=600 m處為切向入射，此間吸收邊界與震源的距離不斷增大，而且隨著距離的增加，地震波的入射角變大，常用方法的吸收邊界產(chǎn)生了較強(qiáng)能量的反射，而本文方法仍然可以有效吸收P波和S波。對(duì)于本文方法，可以得出如下認(rèn)識(shí)：①對(duì)P波的吸收效果優(yōu)于對(duì)S波的吸收效果，因?yàn)镻波單位波長(zhǎng)內(nèi)所含的三角形單元個(gè)數(shù)多于S波單位波長(zhǎng)內(nèi)所含的三角形單元個(gè)數(shù)，精度相對(duì)較高；②使用4階精度的吸收效果優(yōu)于使用3階精度的吸收效果；③完全匹配層越厚，邊界吸收效果越好，但是相應(yīng)的計(jì)算量也增加了。

圖3 均勻彈性介質(zhì)中波場(chǎng)模擬的相對(duì)誤差

2.2 兩層彈性介質(zhì)

在高泊松比情況下，一些數(shù)值模擬方法有可能出現(xiàn)不穩(wěn)定的現(xiàn)象，因此通過(guò)模擬高泊松比兩層彈性介質(zhì)中波的傳播檢驗(yàn)本文方法的有效性。兩層彈性介質(zhì)模型(圖4)大小為1 000 m×700 m。第一層介質(zhì)中，P波速度為1 500 m/s，S波速度為0.1 m/s，密度為1 000 kg/m3，泊松比接近0.5，只能傳播P波。第二層介質(zhì)中，P波速度為2 400 m/s，S波速度為1 100 m/s，密度為1 800 kg/m3。本模型所用震源為脹縮源，只產(chǎn)生P波，震源位于(200 m,200 m)處。波場(chǎng)模擬總時(shí)長(zhǎng)為1 s，時(shí)間步長(zhǎng)為Δt=0.000 1 s。

圖4 兩層彈性介質(zhì)模型

分別使用3階精度非分裂完全匹配層算法及常用的吸收邊界處理方法[6]對(duì)兩層彈性介質(zhì)模型進(jìn)行波場(chǎng)模擬,其結(jié)果見圖5、6。從圖5可以看出：常用方法對(duì)法向入射波具有非常好的吸收效果，垂直入射上邊界的地震波能量完全被吸收，沒(méi)有產(chǎn)生邊界反射波；非垂直入射的地震波能量不能完全被吸收，產(chǎn)生了邊界反射波，而且入射角越大，邊界反射波的能量越強(qiáng)。由傾斜界面反射產(chǎn)生的上行P波在左邊界產(chǎn)生了較強(qiáng)的邊界反射波，由傾斜界面透射產(chǎn)生的下行P波和S波在左邊界和下邊界都產(chǎn)生了較強(qiáng)的邊界反射波。而從圖6可以看出：本文方法對(duì)不同角度入射地震波的吸收效果明顯，不會(huì)產(chǎn)生強(qiáng)能量的反射波。

如果不考慮人工邊界產(chǎn)生的反射波，檢波器應(yīng)該只接收到從兩層介質(zhì)間界面處反射回來(lái)的波，由于第一層介質(zhì)的S波速度近似為0，只有P波能在第一層介質(zhì)中傳播，檢波器只接收到一個(gè)P波波形才是正確的模擬結(jié)果。從圖7可以看出：在常用方法計(jì)算結(jié)果中，3處檢波器在P波前方和后方均接收到了不應(yīng)存在的干擾波波形；而本文方法在高泊松比條件下仍然能夠有效地吸收入射波，3處檢波器均只接收到了一個(gè)P波波形，沒(méi)有受到其他邊界反射波產(chǎn)生的干擾。

圖5 K?ser和Dumbser方法[6]計(jì)算的兩層彈性介質(zhì)波場(chǎng)快照(質(zhì)點(diǎn)速度z方向分量)

圖6 高精度非分裂完全匹配層算法計(jì)算的的兩層彈性介質(zhì)波場(chǎng)快照(質(zhì)點(diǎn)速度z方向分量)

圖7 兩層彈性介質(zhì)中波場(chǎng)模擬檢波器接收到的波形(質(zhì)點(diǎn)速度z方向分量)

3 結(jié)論

本文提出了一種高精度非分裂完全匹配層算法，該算法將非分裂完全匹配層吸收邊界用于任意高階間斷有限元地震波數(shù)值模擬中，使完全匹配層和計(jì)算區(qū)域解的精度一致，能夠達(dá)到任意高階的空間精度和時(shí)間精度。均勻彈性介質(zhì)中波場(chǎng)數(shù)值模擬表明本文算法對(duì)不同角度的入射波具有較好的吸收效果，而兩層彈性介質(zhì)中波場(chǎng)數(shù)值模擬表明本文算法適用于高泊松比介質(zhì)。

[1] 陳可洋.標(biāo)量聲波波動(dòng)方程高階交錯(cuò)網(wǎng)格有限差分法[J].中國(guó)海上油氣,2009,21(4):232-236.

Chen Keyang.High-order staggered-grid finite difference scheme for scalar acoustic wave equation[J].China Offshore Oil and Gas,2009,21(4):232-236.

[2] 范廷恩,余連勇,楊飛龍,等.斜井VSP高斯射線束正演方法[J].中國(guó)海上油氣，2014，26(5):30-35.

Fan Tingen,Yu Lianyong,Yang Feilong,et al.A method of Gaussian beam forward modeling in deviated-well VSP[J].China Offshore Oil and Gas,2014,26(5):30-35.

[3] DORMY E,TARANTOLA A.Numerical simulation of elastic wave propagation using a finite volume method[J].Journal of Geophysical Research:Solid Earth (1978—2012),1995,100(2):2123-2133.

[4] COHEN G,FAUQUEUX S.Mixed finite elements with mass-lumping for the transient wave equation[J].Journal of Computational Acoustics,2000,8(1):171-188.

[5] 李緒宣,于更新,符力耘,等.應(yīng)用邊界元模擬方法分析復(fù)雜海底地震散射特征[J].中國(guó)海上油氣,2011,23(6):357-361.

Li Xuxuan,Yu Gengxin,Fu Liyun,et al.Analysing seismic scattering characteristics of complex seabed by using the boundary-element simulation method[J].China Offshore Oil and Gas,2011,23(6):357-361.

[8] YE Ruichao,DE HOOP M V.A 3D discontinuous Galerkin method for the propagation and scattering of acousto-elastic waves[C]∥2014 SEG Annual Meeting,2014:3329-3333.

[9] TORO E F,MILLINGTON R C,NEJAD L A M.Towards very high order Godunov schemes[M].Springer US,2001:907-940.

[10] BERENGER J P.A perfectly matched layer for the absorption of electromagnetic waves[J].Journal of computational physics,1994,114(2):185-200.

[11] CHEW W C,LIU Q H.Perfectly matched layers for elastodynamics:a new absorbing boundary condition[J].Journal of Computational Acoustics,1996,4(4):341-359.

[12] HASTINGS F D,SCHNEIDER J B,BROSCHAT S L.Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation[J].The Journal of the Acoustical Society of America,1996,100(5):3061-3069.

[13] KOMATITSCH D,MARTIN R.An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[J].Geophysics,2007,72(5):155 -167.

[14] DROSSAERT F H,GIANNOPOULOS A.Complex frequency shifted convolution PML for FDTD modeling of elastic waves[J].Wave Motion,2007,44(7):593-604.

[15] RODEN J A,GEDNEY S D.Convolutional PML (CPML):an efficient FDTD implementation of the CFS-PML for arbitrary media[J].Microwave and Optical Technology Letters,2000,27(5):334-338.

(編輯：馮 娜)

A new algorithm of high precision unsplit perfectly matched layer absorbing boundary for the arbitrary high-order discontinuous Galerkin finite element

He Yangyang1,2Weng Bin1Zhang Jinmiao1

(1.CNOOCResearchInstitute,Beijing100028,China; 2.ChinaUniversityofPetroleum,Beijing102249,China)

The common absorbing boundary condition used commonly for arbitrary high-order discontinuous Galerkin finite element method could not efficiently absorb wave energy of grazing incidence. A new algorithm of high precision unsplit perfectly matched layer absorbing boundary is proposed for the arbitrary high-order discontinuous Galerkin finite element method. The unsplit perfectly matched layer is applied to seismic wave modeling with arbitrary high-order discontinuous Galerkin method, the new equation of the wave propagation in perfectly matched layer and its solution in triangle element are derived, and the discrete scheme is given. This scheme can obtain arbitrary high-order solutions in both time and space in the unsplit perfect match layer with the same solutions in the computational region, thus reducing the energy of reflectivity from absorbing boundary. Simulation results indicate that the method shows better absorbing effect on incident waves with different angles, and it is suitable for the material with high Poisson’s ratio.

unsplit perfectly matched layer; arbitrary high-order discontinuous Galerkin finite element; high precision; incident waves with different angles; high Poisson’s ratio

何洋洋，男，1984年生，工程師，2013年畢業(yè)于西安交通大學(xué)信號(hào)與信息處理專業(yè)，獲博士學(xué)位，主要從事地震波場(chǎng)正演數(shù)值模擬研究工作。地址：北京市朝陽(yáng)區(qū)太陽(yáng)宮南街6號(hào)院海油大廈B座(郵編：100028)。E-mail：heyy9@cnooc.com.cn。

1673-1506(2016)01-0041-07

10.11935/j.issn.1673-1506.2016.01.006

TE132.1+4；P631

A

2015-03-26 改回日期：2015-05-21

*“十二五”國(guó)家科技重大專項(xiàng)“近海大中型油氣田地震勘探技術(shù)(編號(hào)：2011ZX05023-005)” 部分研究成果。

何洋洋,翁斌,張金淼.一種適用于任意高階間斷有限元的高精度非分裂完全匹配層吸收邊界方法[J].中國(guó)海上油氣，2016,28(1)：41-47.

He Yangyang,Weng Bin,Zhang Jinmiao.A new algorithm of high precision unsplit perfectly matched layer absorbing boundary for the arbitrary high-order discontinuous Galerkin finite element[J].China Offshore Oil and Gas,2016,28(1):41-47.