張曉波， 宋　鵬,2,3??， 李金山,2,3， 譚　軍,2,3， 劉照侖， 夏冬明,2,3， 姜秀萍,2,， 趙　波,2

(1.中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 山東 青島 266100；2.中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室， 山東 青島 266100；3.青島海洋科學(xué)與技術(shù)國家實(shí)驗(yàn)室， 山東 青島 266100)

?

理論反射系數(shù)對于完全匹配層邊界的吸收效果影響分析?

張曉波1， 宋鵬1,2,3??， 李金山1,2,3， 譚軍1,2,3， 劉照侖1， 夏冬明1,2,3， 姜秀萍1,2,， 趙波1,2

(1.中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 山東 青島 266100；2.中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室， 山東 青島 266100；3.青島海洋科學(xué)與技術(shù)國家實(shí)驗(yàn)室， 山東 青島 266100)

摘要：本文基于平面波理論和空間二階精度差分格式推導(dǎo)出聲波方程數(shù)值模擬中的完全匹配層邊界條件反射系數(shù)公式，并基于該公式系統(tǒng)分析了理論反射系數(shù)對于完全匹配層吸收效果的影響。實(shí)驗(yàn)結(jié)果表明，對于小入射角度(小于30°)的入射波，當(dāng)理論反射系數(shù)R取10-4時(shí)，完全匹配層有最優(yōu)的吸收效果；對于入射角度為40°、50°和60°的入射波，當(dāng)理論反射系數(shù)R分別取10-5、10-6和10-8時(shí)完全匹配層邊界反射系數(shù)最小；而對更大角度的入射波，需要更小的理論反射系數(shù)R才能使得完全匹配層達(dá)到最優(yōu)的吸收效果。

關(guān)鍵詞：聲波方程；有限差分?jǐn)?shù)值模擬；完全匹配層；理論反射系數(shù)

引用格式：張曉波，宋鵬，李金山，等.理論反射系數(shù)對于完全匹配層邊界的吸收效果影響分析[J]. 中國海洋大學(xué)學(xué)報(bào)(自然科學(xué)版)， 2016,46(6)： 84-89.

ZHANG Xiao-Bo,SONG Peng,LI Jin-Shan,et al. Analysis of the impactof the theoretical reflection coefficient on the perfectly matched layer absorbing boundary[J]. Periodical of Ocean University of China, 2016，46(6)： 84-89.

數(shù)值模擬是研究地震波傳播規(guī)律的有效手段，在地震勘探的采集、處理以及反演等環(huán)節(jié)中有著廣泛應(yīng)用。地震波數(shù)值模擬中，必須用有限的計(jì)算空間來模擬實(shí)際的無限空間，這就引入了人工邊界問題。如果人工邊界處理不當(dāng)，不僅會產(chǎn)生污染中心波場的虛假反射，而且會使地震波的運(yùn)動學(xué)和動力學(xué)特征產(chǎn)生畸變，影響模擬精度，因此，人工邊界的正確處理是地震波高精度數(shù)值模擬的前提。

目前常用的人工邊界處理方法主要有吸收邊界條件(ABC)方法和完全匹配層(PML)邊界條件兩種。ABC方法[1-5]主要通過在邊界處應(yīng)用單程外行波方程來消除邊界反射，由于單程波方程往往具有小角度入射的假設(shè)，這種方法很難對大角度入射波取得滿意的吸收效果。PML方法主要通過在中心波場計(jì)算區(qū)域外加入吸收層，并在吸收層內(nèi)引入衰減因子來達(dá)到消除邊界反射的目的，理論上PML方法不存在對地震波入射角和頻率的假設(shè)，其具有更好的吸收效果，因此PML方法已成為目前應(yīng)用最為廣泛的邊界處理手段[6-11]。

對于PML邊界條件來說，當(dāng)吸收層數(shù)相同時(shí)，其吸收效果取決于PML中的吸收衰減因子和理論反射系數(shù)，衰減因子的設(shè)定及其特性已有大量文獻(xiàn)發(fā)表[12-16]，本文不再贅述。在PML理論反射系數(shù)的研究方面，相關(guān)文獻(xiàn)或討論很少，目前較為常用的方法或思路是：在邊界處理過程中直接給定一個(gè)相對較小的常數(shù)作為理論反射系數(shù)，然后依據(jù)波場快照和合成記錄判斷PML吸收效果。事實(shí)上理論反射系數(shù)的設(shè)置對于PML的吸收效果有重要影響，該系數(shù)的設(shè)定需要理論依據(jù)作為支撐，而這恰好是當(dāng)前國內(nèi)外在PML研究領(lǐng)域的盲區(qū)。

本文基于平面波理論并采用二階差分精度格式，推導(dǎo)出了離散情況下頻率域PML的邊界反射系數(shù)的計(jì)算公式，并以聲波方程數(shù)值模擬為例，分析了理論反射系數(shù)對于PML吸收效果的影響規(guī)律。

1一階速度-應(yīng)力聲波方程及其PML邊界處理

二維各向同性介質(zhì)中的聲波方程為[13]：

(1)

數(shù)值求解式(1)時(shí)，常將模擬空間劃分為3個(gè)區(qū)域：中心波場區(qū)域、PML區(qū)域和角點(diǎn)區(qū)域(見圖1)。

PML區(qū)域的控制方程為[13]：

(2)

其中：Px和Pz分別為P在x和z方向上的分量；dx和dz分別為x和z方向上的衰減因子；其常用表達(dá)式為：

其中：m為PML層內(nèi)的計(jì)算點(diǎn)到中心波場邊界的距離；R為理論反射系數(shù)；L是PML層厚度。實(shí)際計(jì)算時(shí)，在中心波場區(qū)域，dx=0，dz=0。在PML區(qū)域1，dx=0，dz=d(m)。在PML區(qū)域2，dx=d(m)，dz=0；在角點(diǎn)區(qū)域dx=d(m)，dz=d(m)。本文主要以式(3)為基礎(chǔ)討論理論反射系數(shù)R對于PML吸收效果的影響。

2離散條件下頻率域PML的邊界反射系數(shù)計(jì)算公式推導(dǎo)

式(1)在均勻介質(zhì)中有如下形式的解[12]：

(4)

其中：c是縱波波速；θ為入射角；u0為波動幅值；t為時(shí)間；ω為頻率；x和y表示空間坐標(biāo)。以x方向的吸收衰減為例，PML層中聲波方程的解[12]為

(5)

(6)

其中：i為虛數(shù)單位；x0為x方向上與PML層鄰接的中心波場邊界的位置。將式(6)代入式(5)中，然后結(jié)合式(4)可得PML層與中心波場的振幅比為

(7)

當(dāng)波在PML層中傳播并反射回來之后，可以很容易地得到經(jīng)PML層吸收衰減之后的反射系數(shù)為

(8)

將吸收衰減因子d(m)(以(3)式為例)代入式(8)得：

(9)

從式(9)中可以看出，在介質(zhì)連續(xù)的情況下，反射系數(shù)僅與入射角度和理論反射系數(shù)的值有關(guān)，當(dāng)垂直入射時(shí)，反射系數(shù)Rm與理論反射系數(shù)R相等，并且對于任意角度的入射波，R的取值越小，PML邊界的反射系數(shù)越小，即PML邊界的吸收效果最好。但在實(shí)際模擬時(shí)，由于有限差分需要對空間進(jìn)行離散化采樣，空間離散化會對反射系數(shù)產(chǎn)生影響，因此當(dāng)理論反射系數(shù)不斷減小時(shí)，實(shí)際反射系數(shù)并不隨之減小[17]。若要深入分析實(shí)際模擬時(shí)理論反射系數(shù)對于PML邊界的吸收效果影響，需推導(dǎo)出離散條件下的PML邊界反射系數(shù)公式。

因?yàn)镻ML吸收邊界條件在各個(gè)方向上的吸收原理具有一致性，因此本文以x方向的吸收衰減為例推導(dǎo)離散條件下的PML邊界反射系數(shù)公式。令式(2)中的衰減因子dx=d(m)，dz=0，則PML控制方程寫為

(10)

(11)

(12)

為了構(gòu)造一種非分裂的PML表達(dá)式，這里引入Bérenger變量[18]

(13)

將其代入式(12)中得

(14)

(15)

令Δx=Δz=h，將上式寫成空間二階精度有限差分的格式

(16)

(17)

為了方便計(jì)算，這里將Pl,j寫作Pl,j=Ple-i(kzjh)，其中kz=ksinθ，這里k表示波數(shù)，θ表示入射角。然后將其代入式(17)中可得

(18)

(19)

(20)

1)當(dāng)l=-1時(shí)，γ-1=γ-2=1，則有

(21)

2)當(dāng)l=0時(shí)，γ-1=1，則有

(22)

3)當(dāng)l=nl-2時(shí)，則有

(23)

(24)

式中

利用上述方程組，代入式(25)所示的反射系數(shù)的表達(dá)式[18]求出反射系數(shù)，從而得到網(wǎng)格離散情況下理論反射系數(shù)與實(shí)際反射系數(shù)的關(guān)系。

(25)

其中：kx=kcosθ，這里k表示波數(shù)；θ表示入射角。

3理論反射系數(shù)對于PML邊界吸收效果影響分析

以均勻介質(zhì)模型(速度為2 500m/s)為例研究理論反射系數(shù)與實(shí)際反射系數(shù)的關(guān)系，模型網(wǎng)格大小為2 m×2 m，計(jì)算時(shí)頻率取為15 Hz，PML層數(shù)由10層并以10層為間隔依次增至50層，并且將理論反射系數(shù)由10-1并以10-1為比例因子依次降至10-15，得到入射角分別為0°～70°(間隔為10°)情況下的理論反射系數(shù)與PML邊界反射系數(shù)的關(guān)系圖線(關(guān)系曲線如圖2所示，其中紅色實(shí)線、黑色實(shí)線、藍(lán)色實(shí)線、紅色虛線和藍(lán)色虛線分別表示PML層數(shù)為10、20、30、40和50層時(shí)的理論反射系數(shù)與PML邊界反射系數(shù)關(guān)系曲線，綠色虛線指示最小實(shí)際反射系數(shù)所對應(yīng)的理論反射系數(shù)值)。

由圖2可以看出，當(dāng)PML層數(shù)為10、20、30、40或50層時(shí)，由于空間的離散化采樣，并不是理論反射系數(shù)R的取值越小，PML邊界吸收效果越好，而是存在一個(gè)最優(yōu)的R值使得PML邊界的反射系數(shù)達(dá)到極小值，并且不同的入射角度其最優(yōu)的R值不同，對于小入射角度(小于30°)的入射波，當(dāng)理論反射系數(shù)R取10-4時(shí)有最優(yōu)的吸收效果，對于入射角度為40°、50°和60°的入射波，當(dāng)理論反射系數(shù)R分別取10-5、10-6和10-8時(shí)其邊界反射系數(shù)最小，而對于更大角度的入射波，需要更小的理論反射系數(shù)才能達(dá)到最優(yōu)的邊界吸收效果。

4結(jié)論

本文在詳細(xì)論述PML基本原理的基礎(chǔ)上，基于平面波理論和空間二階精度差分格式推導(dǎo)出聲波方程數(shù)值模擬中離散條件下的完全匹配層邊界反射系數(shù)公式，并基于該公式系統(tǒng)分析了理論反射系數(shù)對于PML邊界吸收效果的影響。理論分析和實(shí)驗(yàn)結(jié)果得出結(jié)論如下：

(1)基于聲波方程實(shí)際數(shù)值模擬時(shí)，無論P(yáng)ML層數(shù)為10、20、30、40或50層，由于空間的離散化采樣，并不是理論反射系數(shù)R的取值越小，PML邊界吸收效果越好，而是存在一個(gè)最優(yōu)的R值可以使PML邊界的實(shí)際反射系數(shù)值最小。

(2)對不同的入射角度其最優(yōu)的理論反射系數(shù)R值不同，對于小入射角度(小于30°)的入射波，當(dāng)理論反射系數(shù)R取10-4時(shí)有最優(yōu)的吸收效果，對于入射角度為40°、50°和60°的入射波，當(dāng)理論反射系數(shù)R分別取10-5、10-6和10-8時(shí)其邊界反射系數(shù)最小，而對更大角度的入射波，需要更小的理論反射系數(shù)R才能達(dá)到最優(yōu)的邊界吸收效果。

參考文獻(xiàn)：

[1]Engquist B， Majda A． Absorbing boundary conditions for numerical simulation of waves [J]． Proceedings of the National Academy of Sciences， 1977， 74(5)： 1765-1766．

[2]Engquist B， Majda A． Radiation boundary conditions for acoustic and elastic wave calculations [J]． Communications on Pure and Applied Mathematics， 1979， 32(3)： 313-357．

[3]Higdon R L． Absorbing boundary conditions for difference approximations to the multidimensional wave equation [J]． Mathematics of Computation， 1986， 47(176)： 437-459．

[4]HigdonR L. Numerical absorbing boundary conditions for the wave equation [J]. Mathematics of Computation, 1987, 49(179): 65-90.

[5]宋鵬， 王修田． 優(yōu)化系數(shù)的四階吸收邊界條件 [J]． 中國海洋大學(xué)學(xué)報(bào)(自然科學(xué)版)， 2008， 38(2)： 251-258．

Song Peng， Wang Xiutian． The fourth-order absorbing boundary conditions with optimal parameters [J]． Periodical of Ocean University of China， 2008， 38(2)： 251-258．

[6]Bérenger J P． A perfectly matched layer for the absorption of electromagnetic waves [J]． Journal of Computational Physics， 1994， 114(1)： 185-200．

[7]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．

[8]Wang H X． A Study of damping factors in perfectly matched layers for the numerical simulation of seismic waves [J]． Applied Geophysics， 2013， 10(1)： 63-77．

[9]韓令賀， 何兵壽， 張會星． 含直立裂縫粘彈性介質(zhì)地震波場正演模擬 [J]． 山東科技大學(xué)學(xué)報(bào)(自然科學(xué)版)， 2010， 29(1)： 17-25．

Han Linghe， He Bingshou， Zhang Huixing． Simulation of seismic wave field in viscoelastic media with vertical fractures [J]． Journal of Shandong University of Science and Technology， 2010， 29(1)： 17-25．

[10]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： 3061-3069．

[11]何兵壽， 張會星， 范國苗． 雙程聲波方程逆時(shí)深度偏移 [J]． 物探化探計(jì)算技術(shù)， 2010， 32(2)： 113-119， 107．

He Bingshou， Zhang Huixing，F(xiàn)an Guomiao． Reverse-time depth migration of two way acoustic wave equations [J]． Computing Techniques for Geophysical and Geochemical Exploration，2010，32(2)： 113-119，107．

[12]Collino F and Tsogka C． Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media [J]． Geophysics，2001，66(1)： 294-307．

[13]王守東． 聲波方程完全匹配層吸收邊界 [J]． 石油地球物理勘探，2003，38(1)： 31-34．

Wang Shoudong． Absorbing boundary condition for acoustic waveequation by perfectly matched layer [J]． Oil Geophysical Prospecting，2003，38(1)： 31-34．

[14]陳可洋． 完全匹配層吸收邊界條件研究 [J]． 石油物探， 2010， 49(5)： 473-477．

Chen Keyang． Study on perfectly matched layer absorbing boundary condition [J]． Geophysical Prospecting for Petroleum， 2010，49(5)： 473-477．

[15]李飛，程久龍，楊思通． 基于Cerjan衰減函數(shù)的PML吸收邊界條件 [J]． 物探與化探， 2012， 36(6)： 996-1000．

Li Fei，Cheng Jiulong，Yang Sitong． A study of PML absorbing boundary based on Cerjian attenuation function [J]． Geophysical & Geochemical Exploration， 2012， 36(6)： 996-1000．

[16]李寧． 完美匹配層理論及其在地震波模擬中的應(yīng)用 [D]． 北京： 中國地震局工程力學(xué)研究所，2006．

Li Ning． The theory and the application of PerfectlyMatched Layer in the seismic wave simulation [D]． Beijing： Institute of Engineering Mechanics，China Earthquake Administration，2006．

[17]Wu Z，F(xiàn)ang J． Numerical implementation and performance of perfectly matched layer boundary condition for waveguide structures [J]． IEEE Transactions on Microwave Theory and Techniques，1995， 43(12)： 2676-2683．

[18]Collino F，Monk P B． Optimizing the perfectly matched layer [J]． Computer Methods in Applied Mechanics and Engineering，1998，164(1)： 157-171．

責(zé)任編輯徐環(huán)

Analysis of the Impact of the Theoretical Reflection Coefficient on the Perfectly Matched Layer Absorbing Boundary

ZHANG Xiao-Bo1, SONG Peng1,2,3, LI Jin-Shan1,2,3, TAN Jun1,2,3,LIU Zhao-Lun1, XIA Dong-Ming1,2,3, JIANG Xiu-Ping1,2, ZHAO Bo1,2

(1.College of Marine Geo-science, Ocean University of China, Qingdao 266100, China; 2.Key Lab of Submarine Geosciences and Prospecting Techniques Ministry of Education, Qingdao 266100, China; 3.Qingdao National Laboratory for Marine Science and Technology, Qingdao 266100, China)

Abstract:This paper derives the reflection coefficient formula of the perfectly matched layers (denoted as “PML”)for the numerical simulation of the acoustic wave equation, based on the plane wave theory and the second-order accuracy spatial difference scheme. Then, the impact of the theoretical reflection coefficientRon the absorption effect of the PML is analyzed. The experimental results demonstrate that, for the incident wave with a small incident angle (less than 30°), whenRis equal to 10-4, the PML absorbing boundary condition has the optimal absorption effect.And for the incident wave with the incident angles of 40°, 50° and 60°, when theRis equal to 10-5, 10-6and 10-8respectively, the PML absorbing boundary condition has the smallestactual reflection coefficients.While for the incident wave with a large incident angle, in order to achieve the optimal absorption effect, it is necessary to take a smallerR.

Key words:acoustic wave equation; finite difference numerical simulation; perfectly matched layers; theoretical reflection coefficient

基金項(xiàng)目：? 國家自然科學(xué)基金項(xiàng)目(41574105)資助

收稿日期：2015-05-26；修改日期：2015-06-14

作者簡介：張曉波(1990-)，男，博士生，主要從事地震波正演模擬和屬性分析研究。E-mail：xiaobozhang@outlook.com ??通訊作者：E-mail：pengs@ouc.edu.cn

中圖法分類號：P631.4

文獻(xiàn)標(biāo)志碼：A

文章編號：1672-5174(2006)06-084-06

DOI：10.16441/j.cnki.hdxb.20150190

Supported by the National Natural Science Foundation of China (41574105)