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

何洋洋 翁 斌 張金淼

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

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

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

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

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

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

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

Bérenger[10]在電磁波場模擬中提出了完全匹配層吸收邊界條件(PML)，這是目前已知的效果最好的吸收邊界條件。Chew和Liu[11]以及Hastings等[12]將完全匹配層用于彈性波模擬，取得了好的吸收效果。目前完全匹配層普遍使用波場分裂的方法，以避免計算中出現時間域卷積，但是當地震波切向入射時該方法的邊界反射系數會變大，產生強能量的反射波，對波場模擬結果造成干擾。

Komatitsch和Martin[13]，Drossaert和Giannopoulos[14]等分別提出了一種非分裂形式的完全匹配層算法(CPML)，采用遞歸卷積技術[15]，在提高吸收效果的同時并沒有增加運算量。該方法在完全匹配層內添加了一個濾波器，地震波切向入射時會在匹配層內傳播很長的距離，濾波器會對這種波產生強的吸收衰減，吸收效果優于傳統的完全匹配層算法。

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

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

1.1 完全匹配層內波動方程的建立

二維各向同性彈性介質中的波動方程為

(1)

(2)

(3)

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

在時間域內，非分裂完全匹配層內的波動方程可以表示為[13-14]

(4)

(5)

(6)

(7)

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

本文對式(4)進行改寫，得到

(8)

(9)

式(9)中：S1、S2、S3、S4、S5是式(1)中的震源項，因為完全匹配層內沒有震源，所以此處S1=S2=S3=S4=S5=0。

1.2 完全匹配層內波動方程的求解

(10)

(11)

(12)

(13)

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

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

(14)

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

圖1 完全匹配層內up第k階時間導數的計算步驟

1.3 三角形單元內的數值解

(15)

(16)

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

1.4 離散化

對于式(13)和式(14)，兩邊同時乘以試驗函數Φn(x)，并在x-z坐標系下的三角形單元T(i)內進行積分，然后將積分式映射到ξ-η坐標系下的標準三角形單元TE中。考慮到多項式基函數為正交多項式，式(13)和式(14)在第n個時間步的離散計算格式分別為

]n=

(17)

(18)

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

2 數值算例分析

2.1 均勻彈性介質

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

圖2 均勻彈性介質模型

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

圖3 均勻彈性介質中波場模擬的相對誤差

2.2 兩層彈性介質

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

圖4 兩層彈性介質模型

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

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

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

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

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

3 結論

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

[1] 陳可洋.標量聲波波動方程高階交錯網格有限差分法[J].中國海上油氣,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].中國海上油氣，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] 李緒宣,于更新,符力耘,等.應用邊界元模擬方法分析復雜海底地震散射特征[J].中國海上油氣,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年畢業于西安交通大學信號與信息處理專業，獲博士學位，主要從事地震波場正演數值模擬研究工作。地址：北京市朝陽區太陽宮南街6號院海油大廈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

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

何洋洋,翁斌,張金淼.一種適用于任意高階間斷有限元的高精度非分裂完全匹配層吸收邊界方法[J].中國海上油氣，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.