VTI介质声波方程非分裂式PML吸收边界条件研究

张 衡,刘 洪,李 博,丁仁伟,李丽青,李福元

(1.国土资源部海底矿产资源重点实验室,中国地质调查局广州海洋地质调查局,广东广州510075;2.中国科学院地质与地球物理研究所,中国科学院油气资源研究重点实验室,北京100029;3.中国石油化工股份有限公司石油物探技术研究院,江苏南京211103;4.山东科技大学地球科学与工程学院,山东青岛266590)



VTI介质声波方程非分裂式PML吸收边界条件研究

张 衡1,刘 洪2,李 博3,丁仁伟4,李丽青1,李福元1

(1.国土资源部海底矿产资源重点实验室,中国地质调查局广州海洋地质调查局,广东广州510075;2.中国科学院地质与地球物理研究所,中国科学院油气资源研究重点实验室,北京100029;3.中国石油化工股份有限公司石油物探技术研究院,江苏南京211103;4.山东科技大学地球科学与工程学院,山东青岛266590)

针对垂直横向各向同性(VTI)介质声波波动方程的特点,研究了VTI介质非分裂式完全匹配层(Unsplit perfectly matched layer,UPML)吸收边界条件。首先从VTI介质波场传播稳定性和减少伪横波噪声的角度对常见的几种VTI介质声波波动方程进行了归纳和对比,发现VTI介质声波近似方程更适用于VTI介质情形。然后以VTI介质声波近似二阶波动方程为基础,推导出UPML波动方程,并给出求解该波动方程的具体数值实现方法。数值模拟结果表明,UPML吸收边界条件能达到很好的人工边界反射吸收效果,相比于优化海绵吸收边界条件,其人工边界反射吸收效果明显提高。

VTI介质;非分裂式完全匹配层;伪横波噪声;声波近似方程;二阶波动方程

VTI介质正演数值模拟是逆时偏移和全波形反演的基础[1]。VTI介质参数又称Thomsen各向异性参数(包括纵波各向异性参数ε和变异系数δ),由THOMSEN[2]于1986年引入地震偏移领域,因其能够描述地球介质中大多数各向异性情形而得到广泛应用。THOMSEN[2]在对全球范围内已知的地球介质各向异性数据统计的基础上,归纳出弱各向异性、强各向异性和极端各向异性3种情形,并根据统计结果得出大多数各向异性介质为弱各向异性介质即横向各向同性介质(TI介质)的结论。垂直横向各向同性介质(VTI介质)作为一种最基本的弱各向异性介质,其研究具有重要的理论和实际意义。VTI介质可以看作是具有倾斜对称轴的横向各向同性(TTI)介质在对称轴角度为零时的一种特殊情形。

1996年,TSVANKIN[3]基于Thomsen弱各向异性理论推导得到P-SV波VTI介质精确相速度频散方程。2006年,ZHOU等[4]从TSVANKIN相速度频散方程出发,采用横波速度为零的VTI介质声波近似,推导得到了P-SV波耦合的四阶VTI介质qP波波动方程。因为四阶各向异性波动方程难以求解,ZHOU等[4]引入辅助波场项将四阶各向异性波动方程转换为两个二阶的各向异性波动方程,但是其仍然是一个P-SV波耦合的波动方程,因为各向异性波场传播过程中伪横波(SV波)仍然存在[5],在角度剧变情况下波场传播存在严重的不稳定性问题[6]。为了解决TTI介质声波近似各向异性波场传播存在的不稳定性问题,2009年,FLETCHER等[7]引入适当的横波分量来减小角度剧变时产生的不稳定性,但是同时也带来了较强的伪SV波假象。DUVENECK等[8]从Hooke定律和运动学方程出发推导得到新的VTI介质方程,ZHANG等[9-10]以VTI介质弹性波方程为基础,利用自共轭旋转方式推导了自共轭3D VTI波动方程,但这些方程仍然存在伪SV波假象且在角度剧变时不稳定。HESTHOLM[11]将二阶VTI方程变换成6个一阶耦合偏微分方程的形式,研究了VTI介质的高阶有限差分模拟。LIU等[12]研究了一种新的基于时空域频散关系优化的VTI介质数值模拟算法,相比常规算法,该算法明显提高了计算精度。FOWLER等[13]从特征值矩阵分析的角度对各种VTI介质耦合方程的精度和稳定性进行了综合分析。

VTI介质正演时,吸收边界条件至关重要,因为吸收边界条件的效果会极大影响正演模拟效果。各向异性地震波场正演数值模拟时计算区域有限,必须研究有效的吸收边界条件对强人工边界反射进行有效吸收。目前常用的吸收边界条件包括单程波旁轴近似吸收边界条件[14]和衰减吸收边界条件[15-16]两大类。CLAYTON等[14]提出的单程波旁轴近似吸收边界条件应用广泛,但是其对大角度入射的波吸收效果较差。CERJAN等[15]提出了添加阻尼层的海绵吸收边界条件,通过添加阻尼层达到衰减效果,但是边界吸收效果仍较差。BORDING[17]对CERJAN方法进行改进,提出了优化的海绵吸收边界条件,但是计算效率相对CERJAN方法有所降低。完全匹配层(PML)吸收边界条件最早于1994年由BERENGER[16]在电磁学领域提出,并于1996年由CHEW等[18]将其应用于地震波数值模拟,取得了很好的应用效果。BÉCACHE等[19]指出了各向异性介质应用PML吸收边界条件时的不稳定性。PML吸收边界条件是目前应用效果最好的吸收边界条件,但是其公式推导和编程实现复杂,需要对不同的波动方程推导不同的PML方程形式。PML吸收边界条件技术的研究代表吸收边界技术研究的前沿发展方向。

作为PML吸收边界条件技术研究的一个重要方面,UPML吸收边界条件已经得到了广泛应用。WANG等[20]在弹性波介质模拟中提出了一种新的UPML算法。DROSSAERT等[21]提出了一种基于递归积分复频移UPML算法。QIN等[22]推导得到一种简洁有效的非褶积实现的UPML算法,该算法既不需要对波场分量进行分裂,也不需要做复杂的褶积运算,且具有良好的吸收衰减性能。

本文分析了常见的几种VTI介质波动方程,并从波场传播稳定性和减少伪横波噪声的角度对VTI介质波动方程进行了对比,得出了VTI介质声波近似方程更适用于VTI介质情形的结论。然后从VTI介质声波近似二阶波动方程出发推导得到2D VTI介质时间域UPML波动方程,并给出了2D VTI介质时间域UPML波动方程高阶规则网格有限差分算法求解的具体实现方法。最后利用简单均匀2D VTI介质模型和复杂BP2007 2D VTI介质海洋标准模型对UPML波动方程进行了测试,结果表明本文研究的VTI介质时间域UPML吸收边界条件能很好地吸收边界反射,能够实现VTI介质高精度数值模拟。

1 方法原理

1.1 VTI介质耦合qP波波动方程

P-SV波VTI介质相速度频散方程[3]为:

(1)

式中:vP(θ)表示精确相速度;θ表示相速度传播角度;vSZ表示横波速度;vPZ表示纵波速度;ε和δ表示Thomsen各向异性参数[1]。

直接引入VTI介质声波近似(VTI介质声波近似情况下直接令横波速度vSZ=0)[23],推导得到的2D VTI介质波动方程为[4]:

(2)

本文将采用直接令横波为零的VTI介质声波近似思想推导得到的VTI介质方程称为基于VTI介质声波近似的VTI介质声波近似方程。但是这种方程在实际计算中往往会产生数值不稳定问题,尤其是在角度剧变的情况下[6]。FLETCHER等[7]完全放弃了采用VTI介质声波近似的思路,提出在方程推导中保留横波分量,通过引入横波分量的方式能有效消除VTI介质声波近似产生的不稳定问题,FLETCHER等推导的2D VTI介质纵横波耦合qP波方程为:

(3)

一般情况下α=1,本文σ取经验常数0.75[7]。

本文将FLETCHER等[7]和李博等[24]通过引入横波分量的方式推导得到的VTI介质方程称为VTI介质纵横波耦合方程。

VTI介质波场模拟首先必须满足物理稳定性条件,即P-SV波VTI介质精确相速度频散方程((1)式)中相速度的平方不能小于0[25]:

(4)

由(4)式得到TI介质波场模拟中的物理稳定性条件:

(5)

特别地,在VTI介质声波近似情况下(vSZ=0),即f=1,由(5)式可推导出TI介质声波近似情况下的物理稳定性条件为:

(6)

本文设计了一个楔形2D VTI介质模型进行算法测试,其主要特点是模型左部有一个楔形构造,上层为各向同性介质(图1)。模型参数中横向网格数和纵向网格数都为300,模型横向网格间距和纵向网格间距都为10m。因为本例震源置于各向同性的海水层中,VTI数值模拟时不会出现影响qP波场传播的伪SV波噪声,此种情况下不需要采取在震源处添加各向同性薄层消除伪SV波噪声的策略[5,26]。

数值模拟时,雷克震源子波主频取25Hz,传播时间为3s,时间采样间隔为0.5ms,吸收边界条件采用优化海绵吸收型[17],分别采用VTI介质声波近似方程(ZHOU等的方程)和引入横波分量的VTI介质纵横波耦合方程(FLETCHER方程)。

对波场进行分析,从t=1.0s和t=1.5s时刻的波场快照(图2,图3)看出,ZHOU等的方程都能保持稳定传播而且没有出现影响qP波场传播的伪SV噪声,而FLETCHER等的方程由于引入了横波分量解决了地震波在各向异性介质中传播的不稳定性问题,但是也带来了较多的伪SV噪声(图2b和图3b中红色箭头处)。因此VTI介质情况下更适合采用基于VTI介质声波近似方程,因此在二维和三维VTI介质正演中,我们将采用稳定的ZHOU等的方程。原因是VTI介质的对称轴角度为零,在VTI介质声波近似假设下基本不存在稳定性问题,另外VTI介质声波近似方程相比VTI介质纵横波耦合方程而言伪SV噪声也更少。

综上所述,基于VTI介质精确相速度频散方程推导出的VTI介质波动方程可以分为两类,一类为基于VTI介质声波近似的VTI介质声波近似方程,另一类为考虑横波分量的VTI介质纵横波耦合方程,其中VTI介质声波近似方程可以视为VTI介质纵横波耦合方程在横波速度为零时的特例。通过研究发现,VTI介质声波近似方程更适用于VTI介质。

图1 楔形2D VTI介质模型各向异性参数a 纵波速度vP(km/s); b 纵波各向异性参数ε; c 变异系数δ

图2 楔形2D VTI介质模型t=1.0s时的波场快照a ZHOU等的方程; b FLETCHER方程

图3 楔形2D VTI介质模型t=1.5s时的波场快照a ZHOU等的方程; b FLETCHER方程

1.2 2D VTI介质时间域UPML波动方程推导

VTI介质波动方程的特点是不含比较复杂的交叉导数项,因此可以采用计算效率更高的UPML吸收边界条件来实现。针对VTI介质波动方程上述特点,我们以ZHOU等的二阶2D VTI介质声波近似时间域波动方程(公式(2))为基础,推导2D VTI介质时间域UPML波动方程。

首先对二阶2D VTI介质时间域波动方程(公式(2))做傅里叶变换,得到其频率域波动方程的形式:

(7)

再引入X和Z方向的频率域PML拉伸函数ξx,ξz:

(8)

其中,dx,dz为PML衰减因子,其表达式为[27]:

(9)

式中:LPML是PML边界层的厚度;θ是计算点距PML边界的距离;Rcoeff是理论反射系数值,本文Rcoeff取0.0001。

(10)

对公式(10)求导并结合公式(8)推导得到:

(11)

在正常计算区域内,加PML衰减项的声波波动方程与不加PML衰减项的声波波动方程是一致的,因为在正常计算区域内di=0。如果PML计算区域内di>0,则在PML区域内声波波动方程的平面波解为指数衰减的形式:

(12)

当波从内部区域过渡到PML边界层后波场逐渐衰减。

将公式(11)代入2D VTI介质频率域波动方程(公式(7)),并引入6个辅助波场项Epx,Epz,Eqx,Hpx,Hpz,Hqx进行推导,最后得到2D VTI介质频率域UPML波动方程,将此方程反傅里叶变换到时间域即可得到2D VTI介质时间域UPML波动方程:

(13)

其中6个辅助波场项Epx,Epz,Eqx,Hpx,Hpz,Hqx对应的控制方程为:

公式(13)的右端项由两部分组成,一部分是正常项(不含辅助波场项),另一部分是衰减项(含辅助波场项)。如果衰减项为零,则UPML波动方程就退化为常规的2D VTI介质波动方程(公式(2))。

通过一系列公式推导可以发现不同2D VTI介质方程的UPML方程的形式不一样,差分离散形式和辅助波场项的表达形式也不一致,其具体形式根据不同的VTI介质方程形式而定。本文只讨论二维VTI介质情形,但是这种方法也可以扩展到三维VTI介质情形(详细推导见附录A),在此不再展开讨论。

1.3 2D VTI介质UPML波动方程数值实现

从二维UPML波动方程数值实现示意图(图4)可以看出,UPML具体实现的时候不再有如分裂式完全匹配层(SPML)波动方程一样的边和角的概念,只有方向的概念,我们只需要考虑X和Z2个方向即上下左右4个PML区域波场的衰减,而且UPML不需要对波场项进行方向分裂,相比二维SPML需要考虑8种情形并对波场项进行方向分裂而言,UPML简化了计算量,也减少了存储空间,编程工作量大大减少。

图4 二维UPML数值实现示意

2D VTI介质UPML波动方程具体数值实现过程如下:

1) 在所有区域进行波场计算。采用12阶高阶规则网格有限差分算法对公式(13)的正常项进行波场计算。

2) 在PML区域内进行波场更新。在PML计算区域内对公式(13)的衰减项进行有限差分计算并更新PML计算区域内的波场,共考虑4个PML区域情形,其中4个角域需要计算2次。

3) 在PML区域内更新辅助波场项Epx,Epz,Eqx,Hpx,Hpz,Hqx。在PML计算区域内对公式(13)的6个辅助波场项的控制方程进行有限差分计算并更新PML计算区域内的辅助波场项,共考虑4个PML区域情形,其中4个角域需要计算2次。

2 数值计算实例

2.1 简单均匀2D VTI介质模型

为了验证本文提出的VTI介质UPML吸收边界条件的有效性,首先采用简单均匀2D VTI介质模型进行UPML点源响应测试。模型参数中横向网格数和纵向网格数都为256,模型横向网格间距和纵向网格间距都为10m,速度为常速3000m/s。VTI介质的各向异性参数ε=0.3,δ=0.1。数值模拟时,雷克震源子波主频取25Hz,震源位置设置在模型的中心,传播时间为1.5s,采样间隔为1.5ms,采用12阶高阶规则网格有限差分格式进行模拟。采用UPML吸收边界条件,PML边界吸收层数分别设为10层、20层和30层。图5为3种不同PML吸收层数情形下,t=0.54s时刻的正演波场快照,可以看出,PML吸收层数为10层时边界反射吸收效果不理想(图5a);PML吸收层数为20层时边界反射吸收效果明显改善,但是仍然存在着微弱的边界反射(图5b);当PML吸收层数为30层时边界反射已经基本被消除(图5c)。因此本例选定PML层数为30层。图6a到图6f分别是UPML波动方程在各个时刻的波场快照,由图可见UPML对边界反射吸收效果明显。

2.2 复杂BP2007 2D VTI介质海洋标准模型

以复杂BP2007 2D VTI介质海洋标准模型为例,并且与优化海绵吸收边界条件[17]进行对比，进一步说明VTI介质UPML吸收边界条件的应用效果。数值模拟时采用截取的部分BP2007 2D VTI介质海洋标准模型(图7),该部分主要模拟近海岸的地质构造,其特点是左边有一个被VTI各向异性地层包围的高速各向同性盐丘。模型参数包括纵波速度vP,Thomsen纵波各向异性参数ε和Thomsen变异系数δ,模型区域大小为:横向网格数nx=601,纵向网格数nz=361,模型横向网格间距和纵向网格间距都为6.25m,其中海水层为各向同性介质。数值模拟时,雷克震源子波主频取25Hz,采用UPML吸收边界条件,PML边界设为30层,传播时间为2s,时间采样间隔为0.4ms,采用12阶高阶规则网格有限差分格式进行模拟。因为本例震源置于各向同性的海水层中,VTI数值模拟时不会出现影响qP波场传播的伪SV波噪声,此种情况也不需要采取在震源处添加各向同性薄层来消除伪SV波噪声的策略[5,26]。从BP2007 2D VTI介质海洋标准模型的正演单炮记录来看,UPML消除边界反射效果良好(图8a),而优化海绵吸收边界条件仍然有着比较强的边界反射(图8b 红色箭头标示),UPML吸收边界条件明显优于优化海绵吸收边界条件的边界反射吸收效果。

图5 二维UPML波动方程t=0.54s时刻不同PML吸收层数情况下正演波场快照a PML吸收层数为10层; b PML吸收层数为20层; c PML吸收层数为30层

图6 均匀2D VTI介质模型采用二维UPML波动方程不同时刻的正演波场快照a t=0.39s; b t=0.45s; c t=0.51s; d t=0.57s; e t=0.90s; f t=1.35s

图7 截取的部分BP2007 2D VTI介质海洋标准模型a 纵波速度vP; b 纵波各向异性参数ε; c 变异系数δ

图8 BP2007 2D VTI介质海洋标准模型正演单炮记录a UPML吸收边界条件方法; b 优化海绵吸收边界条件方法

3 结束语

本文对常见的几种VTI介质波动方程进行了研究,发现基于VTI介质精确相速度频散方程推导出的VTI介质波动方程可以分为基于VTI介质声波近似的VTI介质声波近似方程和考虑横波分量的VTI介质纵横波耦合方程,就波场传播稳定性和减少伪横波噪声而言,VTI介质声波近似方程更适用于VTI介质情形。从基于VTI介质声波近似的二阶波动方程出发,推导得到了VTI介质UPML波动方程。研究发现不同2D VTI介质UPML方程的形式不一样,差分离散形式和辅助波场项的表达形式也不一致,其具体形式根据不同的VTI方程形式而定。利用简单均匀2D VTI介质模型和复杂BP2007 2D VTI介质海洋标准模型对2D VTI介质UPML方程进行了数值模拟,结果表明,本文研究的UPML吸收边界条件能取得很好的人工边界反射吸收效果,相比优化海绵吸收边界条件,其人工边界反射吸收效果明显提高。在附录中,我们给出了3D VTI介质时间域UPML波动方程的推导结果,从而可以将本文方法进一步推广到三维VTI介质的高精度模拟中。

[1] WARNER M,RATCLIFE A,NANGOO T,et al.Anisotropic 3D full-waveform inversion[J].Geophysics,2013,78(2):R59-R80

[2] THOMSEN L.Weak elastic anisotropy[J].Geophysics,1986,51(10):1954-1966

[3] TSVANKIN I.P-wave signatures and notation for transversely isotropic media:an overview[J].Geophysics,1996,61(2):467-483

[4] ZHOU H B,ZHANG G Q,BLOOR R.An anisotropic acoustic wave equation for VTI media[J].Expanded Abstracts of 68thEAGE Annual Conference,2006:H033

[5] GRECHKA V,ZHANG L,RECTOR J W.Shear waves in acoustic anisotropic media[J].Geophysics,2004,69(2):576-582

[6] BAKKER P M,DUVENECK E.Stability analysis for acoustic wave propagation in tilted TI media by finite differences[J].Geophysical Journal International,2011,185(2):911-921

[7] FLETCHER R P,DU X,FOWLER P J.Reverse time migration in tilted transversely isotropic (TTI) media[J].Geophysics,2009,74(6):179-187

[8] DUVENECK E,MILCIK P,BAKKER P M.Acoustic VTI wave equations and their application for anisotropic reverse-time migration[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008:2186-2190

[9] ZHANG Y,ZHANG H Z.A stable TTI reverse time migration and its implementation[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:2794-2798

[10] ZHANG Y,ZHANG H Z,ZHANG G Q.A stable TTI reverse time migration and its implementation[J].Geophysics,2011,76(3):WA3-WA11

[11] HESTHOLM S.Acoustic VTI modeling using high-order finite differences[J].Geophysics,2009,74(5):T67-T73

[12] LIU Y,SEN M K.Acoustic VTI modeling with a time-space domain dispersion-relation-based finite-difference scheme[J].Geophysics,2010,75(3):A11-A17

[13] FOWLER P J,DU X,FLETCHER R P.Coupled equations for reverse time migration in transversely isotropic media[J].Geophysics,2010,75(1):S11-S22

[14] CLAYTON R,ENQUIST B.Absorbing boundary conditions for acoustic and elastic wave equations[J].Bulletin of the Seismological Society of America,1977,67(6):1529-1540

[15] CERJAN C,KOSLOFF D,KOSLOFF R,et al.A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J].Geophysics,1985,50(4):705-708

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

[17] BORDING R P.Finite difference modeling-nearly optimal sponge boundary conditions[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004:1921-1924

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

[20] WANG T,TANG X M.Finite-difference modeling of elastic wave propagation:A nonsplitting perfectly matched layer approach[J].Geophysics,2003,68(5):1749-1755

[21] DROSSAERT F H,GIANNOPOULOS A.A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[J].Geophysics,2007,72(2):T9-T17

[22] QIN Z,LU M H,ZHENG X D,et al.The implementation of an improved NPML absorbing boundary condition in elastic wave modeling[J].Applied Geophysics,2009,6(2):113-121

[23] ALKHALIFAH T.Acoustic approximations for processing in transversely isotropic media[J].Geophysics,1998,63(2):623-631

[24] 李博,李敏,刘红伟,等.TTI介质有限差分逆时偏移的稳定性探讨[J].地球物理学报,2012,55(4):1366-1375 LI B,LI M,LIU H W,et al.Stability of reverse time migration in TTI media[J].Chinese Journal of Geophysics,2012,55(4):1366-1375

[25] 李博.地震偏移成像和GPU超算技术[D].北京:中国科学院大学,2011 LI B.Seismic migration and GPU supercomputing technology[D].Beijing:University of Chinese Academy of Sciences,2011

[26] 张衡,刘洪,唐祥德,等.基于平均导数优化方法的VTI介质频率空间域正演[J].地球物理学报,2015,58(9):3306-3316 ZHANG H,LIU H,TANG X D,et al.Forward modeling of VTI media in the frequency-space domain based on an average-derivative optimal method[J].Chinese Journal of Geophysics,2015,58(9):3306-3316

[27] COLLINO F,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

附录A 3D VTI介质时间域UPML波动方程

3D VTI介质时间域UPML波动方程的推导基于ZHOU等的3D VTI介质声波近似波动方程:

(A1)

从公式(A1)出发,对方程做傅里叶变换推导得到其频率域波动方程的形式:

(A2)

引入X,Y和Z方向的频率域PML拉伸函数ξx,ξy,ξz:

(A3)

其中:dx,dy,dz为PML衰减因子,其表达式为:

(A4)

式中:LPML是PML边界层的厚度;θ是计算点距PML边界的距离;Rcoeff是理论反射系数值,取0.0001。

(A5)

对公式(A5)求导并结合公式(A3)推导得到:

(A6)

将公式(A6)代入3D VTI介质频率域波动方程(公式(A2)),并引入10个辅助波场项Epx,Epy,Epz,Eqx,Eqy,Hpx,Hpy,Hpz,Hqx,Hqy,最后得到3D VTI介质频率域UPML波动方程,将此方程反傅里叶变换到时间域即可得到3D VTI介质时间域UPML波动方程:

(A7)

其中10个辅助波场项Epx,Epy,Epz,Eqx,Eqy,Hpx,Hpy,Hpz,Hqx,Hqy对应的控制方程为:

分析公式(A7)的右端项,其由两部分组成,一部分是正常项(不含辅助波场项),另一部分是衰减项(含辅助波场项)。如果衰减项为零,则UPML波动方程就退化为常规的3D VTI介质波动方程(公式(A1))。

通过一系列公式推导研究可以发现不同3D VTI介质方程的UPML方程的形式不一样,差分离散形式和辅助波场项的表达形式也不一致,其具体形式根据不同的VTI介质方程形式而定。

对于三维UPML数值实现而言,只需要考虑X,Y和Z3个方向即上下左右前后共6个PML区域波场的衰减,而且不需要对波场项进行方向分裂,相比三维SPML需要考虑26种情形并对波场项进行方向分裂而言,UPML简化了计算量,也减少了存储空间,编程工作量大大减少。

(编辑：朱文杰)

The research on unsplit PML absorbing boundary conditions of acoustic equation for VTI media

ZHANG Heng1,LIU Hong2,LI Bo3,DING Renwei4,LI Liqing1,LI Fuyuan1

(1.MLRKeyLaboratoryofMarineMineralResources,GuangzhouMarineGeologicalSurvey,ChinaGeologicalSurvey,Guangzhou510075,China;2.KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China;3.SinopecGeophysicalResearchInstitute,Nanjing211103,China;4.CollegeofEarthScienceandEngineering,ShandongUniversityofScienceandTechnology,Qingdao266590,China)

We study unsplit perfectly matched layer (UPML) absorbing boundary conditions for vertical transversely isotropic (VTI) media based on the characteristics of VTI media acoustic wave equation.We firstly summarize several common VTI media acoustic wave equations and compare these wave equations in detail from the perspective of wavefield propagation stability and pseudo SV-wave noise reduction for VTI media.We conclude that the VTI acoustic approximation equation is applicable for VTI media.Then we derive the UPML wave equation from the VTI acoustic approximation second-order wave equation.Afterwards we give the concrete numerical implementation method.The numerical modeling examples show that the UPML boundary condition can absorb the artificial boundary reflection very well.The boundary reflection absorbing effect of the proposed UPML algorithm is much better compared with the optimal sponge boundary condition.

VTI media,unsplit perfectly matched layer (UPML),pseudo SV-wave noise,acoustic approximation equation,second-order wave equation

2016-02-15;改回日期：2016-05-24。

张衡(1987—),男,博士,工程师,主要从事复杂各向异性介质地震波场正演数值模拟和全波形反演研究工作。

国土资源部海底矿产资源重点实验室开放基金(No.KLMMR-2015-A-14)、国家自然科学基金项目(41604110)、国家高技术研究发展计划(863计划)项目(2013AA092501)和国家科技重大专项(2011ZX05003-003)联合资助。

This research is financially supported by the Open Fund of MLR Key Laboratory of Marine Mineral Resources (Grant No.KLMMR-2015-A-14),the National Natural Science Foundation of China (Grant No.41604110),the National High-tech R & D Program of China (863 Program) (Grant No.2013AA092501) and the National Science and Technology Major Project of China (Grant No.2011ZX05003-003).

P631

A

1000-1441(2016)06-0781-12

10.3969/j.issn.1000-1441.2016.06.002