2020, Vol. 63
2. 海洋国家实验室海洋矿产资源评价与探测技术功能实验室, 青岛 266071;
3. 中海石油(中国)有限公司深圳分公司, 深圳 518067
2. Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China;
3. Shenzhen Branch of CNOOC(China) Limited, Shenzhen 518067, China
广义的非均匀介质是指介质的弹性参数、密度参数是空间坐标的函数，即介质单元体的尺度大于地震波能够识别的尺度，其性质随空间位置变化并且能够被区分.如果这种介质的弹性参数可随空间坐标连续变化，则称为连续介质.利用弹性力学广义胡克定律方程、柯西方程和牛顿定律的微分方程建立适应弹性参数、密度参数空变的弹性波方程，在非均匀介质中弹性波方程是表达纵、横波耦合的波动方程.由于弹性波场是矢量，因此使用有限差分法模拟波传播时需要计算波场所有分量，在计算效率和内存上比较“昂贵”.此外，时间步长和网格间距选择须参考模型中的最慢速度，即横波速度(Alkhalifah, 2000).为了降低计算成本，本文在充分考虑非均匀介质的纵波速度、密度空间变化率的前提下，尝试用标量波方程近似非均匀介质中纵波传播，突破传统常速度、常密度标量波方程的假设.
地震波在非均匀介质中传播是耦合传播的，根据牛顿第二定律的微分形式，要考虑应力的空间变化率，进而涉及到弹性模量(纵、横波速度参数)、密度参数的空间变化率，不同于均匀介质条件下弹性波传播.在均匀介质中，速度、密度参数是常数(不随空间坐标变化)，不需要考虑其空间变化率，利用涡流定理，定义纵波位移位和横波位移位，根据规范不变性条件得到分离的纵横波位移位方程(孙成禹，2007).Hook(1961)对非均匀弹性波方程进行解耦，通过先寻找能够数学上分离非均质弹性波位移的新的标量位和矢量位，然后进行解耦，且只在特定条件下成立.Richards(1974)针对球对称弹性非均匀介质，给出P-SV-SH波的势函数，并在特定情况下实现纵、横波解耦，其结果与Hook工作类似.对于非均匀弹性波位移的解耦，前人做了大量的工作，但是都未能得出理想的解耦结果，根本原因在于均匀介质假设下，基于涡流定理分离得到的纵、横波位移位在非均匀介质假设下不是无散场和无旋场，即位移位不再代表纵波和横波.
Alkhalifah(2000, 2003)在研究各向异性介质弹性波方程解耦时提出声学近似的思想:在弹性波耦合传播的各向异性介质中，另横波不传播(令横波速度为零)，分离出描述纵波独立传播的高精度波动方程，简化了各向异性介质纵波传播的求解问题.本文引用Alkhalifah声学近似思想对非均匀介质弹性波方程进行解耦，建立非均匀介质纵波速度、密度等空间双变参数的标量波传播方程，为双变参数标量波全波形反演奠定好正演基础.
当前地震波场数值模拟算法数值模拟广泛应用于地球物理和工程应用，有限差分法是最常用的方法之一(Kelly et al., 1976; Robertsson et al., 1994a, b)，其具有计算效率高、内存占用小、方便实现等优点.Boore(1972)提出非均匀介质二阶弹性波有限差分法，Kelly等(1976)对这一方法进行了改进. Madariaga(1976)提出一阶速度-应力弹性波方程交错网格有限差分法，Virieux(1984, 1986)和Levander(1988)将其进行推广并实现四阶差分精度.刘洋等(1998)提出任意偶数阶精度差分格式.董良国等(2000)基于交错网格提出时间四阶、空间高阶有限差分法.殷文等(2006)在频率域实现了有限差分法.
本文从弹性力学中建立弹性波动方程的三个基本方程出发，考虑非均匀介质弹性参数的空变特征，建立非均匀介质的弹性波动方程，利用声学近似思想建立位移表征的纵波波动方程，利用本征值问题求解方法建立标量波频率-波数域传播算子，从而建立描述纵波传播的标量波方程.随后推导含PML边界波动方程差分格式并建立不同模型数值模拟进行数值试算，与均匀假设标量波方程和变密度方程对比证明本方法的准确性和稳定性.
1 非均质双变参数标量纵波方程建立从弹性动力学中的几何方程、本构方程和牛顿运动微分方程三个基本方程(Aki and Richards, 2002)出发建立非均匀介质的弹性波方程.弹性动力学中表示微观弹性应变和位移关系的几何方程又称柯西方程，公式为
|
(1) |
各向同性介质的本构方程，又称广义虎克定律，公式为
|
(2) |
牛顿运动微分方程，又称纳维尔方程，公式为
|
(3) |
公式(1)、(2)和(3)为弹性动力学中建立波动方程所需要的三个基本方程，式中εij为应变张量，σij为应力张量，ui为位移分量，üi为位移分量时间二阶导数，ui, j位移分量ui对空间xj分量一阶导数.结合弹性波速度与弹性模量、密度之间的关系，利用弹性动力学的三个基本方程推导出非均匀介质的无源弹性波方程，公式为
|
(4) |
方程(4)是速度、密度参数及位移分量表征的非均匀介质弹性波动方程，波动方程的系数决定波传播的相速度，本方程考虑了速度、密度参数的空间变化率.接下来利用Alkhalifah声学近似思想解耦非均匀介质弹性波动方程为标量波方程，为此令VS=0，由式(4)推导成为非均匀介质只有纵波位移分量表征的波动方程，公式为
|
(5) |
公式(5)描述的是纵波在非均匀介质中传播的波动方程，这时设想用标量方程解决纵波在非均匀介质中传播问题，相比较均匀介质的波动方程，本方程传播算子
|
(6) |
其中 u =(u1, u2, u3)，矩阵Γ为Kelvin-Christoffel矩阵，表达式为
|
(7) |
其中，系数
|
(8) |
将式(8)算子作用于波场P，进行傅氏反变换，建立了描述非均匀介质纵波传播的双变参数标量波动方程，公式为
|
(9) |
均匀介质标量波方程表达为
|
(10) |
非均匀介质标量波方程(式(9))与均匀介质标量波方程(式(10))相比，增加了速度、密度的空间导数项，这是因为综合考虑了介质参数的空间变化率，基于介质均匀假设条件下，二者能够实现相互转化，二者震源项关系为ρF=f.
2 完全匹配层吸收边界条件在地震数值模拟的过程中，由于研究区域是有限的，不可避免会产生边界反射，因此合理有效的处理边界反射是数值模拟的重心之一.完全匹配层(PML)吸收边界条件(Berenger，1994)已经成为时间域有限差分地震波数值模拟主流的吸收边界条件.它的优势在于可以高效吸收任意入射角度任意频率入射的波并且几乎不产生任何反射(张岩和吴国忱，2013).PML相对于其他吸收边界条件是数值稳定的且所需层数相对较少(Wang and Tang, 2003)，吸收效果也优于其他类的吸收条件(Higdon，1987).目前，PML吸收边界条件已经广泛的应用于声波和弹性波(Liu and Tao, 1997；刘有山等，2012)以及双相介质(Zeng and Richards, 2001)的地震波数值模拟中，并且取得了良好的效果.李振春等(2013)对PML边界非分裂形式进行了研究，进一步减少了计算内存.李青阳等(2018)基于PML吸收边界建立时间四阶精度的弹性波数值模拟方法.为了对文中提出的标量波方程进行更精确的数值模拟，我们通过推导将PML吸收边界条件引入此方程中.
将方程(9)去掉一维得频率域二维双变参数标量方程，公式为
|
(11) |
其中，ω为圆频率.
将空间坐标系变换到复数坐标系，复数坐标系下空间变量表示为(王守东, 2003; 王永刚等, 2007):
|
(12) |
其中，n表示空间坐标x、z；
由式(12)可以得到二阶偏导数新旧坐标系下的对应关系为
|
(13) |
对式(13)在时间上做傅里叶反变换，这样我们就得到了时间域完全匹配层控制方程.但直接做傅里叶反变换时需要做时间域的褶积.因此对波场分裂后进行傅里叶反变换得到时间域的完全匹配层控制方程为
|
(14) |
其中dx(x)和dz(z)分别表示沿x、z方向的衰减因子.
因此方程(14)的截断误差为O(Δx2N, Δz2N, Δt2)的高阶差分方程可表示为
|
(15) |
本节主要从震源加载、数值离散以及稳定性条件对新方程进行分析，本文对方程(9)数值模拟，以及后面涉及到的常密度标量波方程(方程(10))和变密度标量波方程(方程(16))，为统一比较，均选择规则网格离散.在震源加载过程中，方程(9)震源项与其他方程(方程(10)、(16))差一个密度值，加载方式均在一个点上.方程(16)为
|
(16) |
在实际模拟过程中方程(9)和方程(10)标量场P与密度ρ重合，以x方向为例，对于∂2x差分权重落在相邻网格点中间，而对于∂x的差分权重落在网格点本身.方程(16)在空间差分中存在
接下来分析方程(9)的稳定性条件，通过将方程(9)离散后变换到频率域，然后建立差分格式的增长矩阵的特征方程，令方程的模小于等于1，即得到稳定性条件:
|
(17) |
本文建立方程针对的目标为连续介质，等式右端项中速度参数和密度参数的梯度远小于其本身，当介质均匀情况下，方程(17)中第二项考虑空间的Nyquist频率，取m=1，则得到与常规标量波方程一致的稳定性条件:
|
(18) |
事实上方程(17)中系数
为了分析新方程的波场特征，本节选取传统标量波方程和变密度标量波方程进行对比.并通过均匀介质模型、垂向连续非均匀介质模型和弹性间断的层状介质模型三个代表性模型一一展开对比讨论.
4.1 均匀介质模型均匀介质模型网格点数为200×200，空间步长10 m，时间采样间隔0.8 ms，震源采用主频20 Hz的雷克子波，坐标为(100, 100).模型参数包括纵波速度和密度，其中纵波速度为3000 m·s-1，密度为1500 kg·m-3.
前文指出当介质均匀(速度和密度均为常数)时，双变参数标量波方程(Two Variable Parameter Scalar Wave Equation，简写为TVSE)即退化为传统标量波方程(Traditional Scalar Wave Equation，简写为TE)，这一点变密度标量波方程(Density Variable Scalar Wave Equation，简写为DVE)同样适用.
图 1显示了分别采用TE、DVE和TVSE正演模拟并传播0.28s的波场快照，通过对比显然在介质均匀条件下，三个方程模拟下波传播特征基本相同.对波场进行一维采样，分别抽取横向x=1000 m和纵向z=1000 m的两条线记录对比，如图 2所示，黑色实线为TE模拟波场，红色虚线为DVE模拟波场，蓝色圆圈表示TVSE模拟波场，对比结果也证实该结论成立.
|
图 1 均匀介质模型波场快照对比 Fig. 1 Wavefield snapshot comparison of homogeneous media model (a) TE; (b) DVE; (c) TVSE. |
|
图 2 波场一维采样对比(a)波场在x方向的采样(z=1000 m)；(b)波场在z方向的采样(x=1000 m). Fig. 2 One-dimensional sampling of wavefield (a) x direction (z=1000 m); (b) z direction (x=1000 m). |
连续介质设定垂向连续非均匀，横向均匀模型，网格点数为200×200，空间步长10 m，时间采样间隔0.8 ms，震源采用主频20 Hz的雷克子波，坐标为(100, 100).模型参数如图 3所示，纵波速度(图 3a)为2×Depth+2000 m·s-1，密度(图 3b)为4×Depth+1000 kg·m-3.
|
图 3 连续介质模型(a)速度；(b)密度. Fig. 3 Continuous media model (a) Velocity; (b) Density. |
首先固定速度为常数3000 m·s-1，而密度垂向随空间变化，线性增长函数为: 4×Depth+1000 kg·m-3.图 4是分别采用TE、DVE和TVSE正演模拟并传播0.28 s的波场快照，显然TE波场值相比于DVE和TVSE偏大，而且密度的变化没有引起波走时信息的变化.图 5对三个波场进行横向(x=1000 m)和纵向(z=1000 m)一维采样，分别抽取两条线记录对比，图 5a是水平抽样，由于水平均匀性，波场信息基本一致；图 5b是垂直采样，显然在1700 m处，DVE和TE的振幅要高于TVSE，而300 m处则DVE和TE的振幅低于TVSE，反观TE结果在1700 m和300 m处振幅值基本相同.
|
图 4 速度固定，密度变化模型波场快照对比 Fig. 4 Wavefield snapshot comparison for constant velocity and variable density model (a) TE; (b) DVE; (c) TVSE. |
|
图 5 波场一维采样对比(a)波场在x方向的采样(z=1000 m)；(b)波场在z方向的采样(x=1000 m). Fig. 5 One-dimensional sampling of wavefield (a) x direction (z=1000 m); (b) z direction (x=1000 m). |
对于这一现象，可以从垂直入射的Zoeppritz方程分析，首先给出平面波纵波入射(如图 6所示)下的Zoeppritz方程:
|
图 6 平面纵波入射的反射和透射 Fig. 6 Reflection and transmission of plane P-wave incidence |
|
(19) |
公式(19)为利用位移振幅比表示的反射透射系数方程，在纵波垂直入射情况下，入射角α=0，根据snell定律，有
|
(20) |
事实上，连续介质可以近似为参数逐层均匀递增的层状介质，而速度固定情况下，式(20)退化为
|
(21) |
当密度不变情况下，constant=0，Tpp=1；当密度随深度逐渐增大，constant>0，Tpp < 1；当密度随深度逐渐增大，constant < 0，Tpp>1，根据连续介质模型(图 3)计算得到震源在(1000 m，1000 m)时零偏移距两侧的反射透射系数，如图 7所示.
|
图 7 零偏移距(垂直入射)反射透射系数对比(a)上行波透射系数; (b)上行波反射系数; (c)下行波透射系数; (d)下行波反射系数. Fig. 7 Comparison of zero offset (vertical incidence) reflection transmission coefficients (a) Rpp of up-going wave; (b) Tpp of up-going wave; (c) Rpp of down-going wave; (d) Tpp of down-going wave. |
反射透射系数是基于连续介质的层状假设计算的，连续介质并不会产生反射波，这里只关注透射系数对于波传播过程中振幅的影响.根据透射系数可以得到不同深度波的最大振幅.如图 8所示，向上传播过程中振幅是逐渐增大的(Amplitude>1)，向下传播波动振幅逐渐减小(Amplitude < 1).在0.28 s时刻的垂向零偏移距抽样中(图 5)，300 m处TVSE振幅大于TE振幅，而TE振幅大于DVE振幅，三者振幅值均小于1，这是因为波在2D介质传播过程中本身就存在球面扩散.而TE传播介质是均匀的，相当于Tpp=1的特殊层状假设，对比可知，TVSE模拟下波的传播振幅特征与Zoeppritz方程一致，相反，DVE方程结果不正确.在1800 m深度处TVSE振幅小于TE振幅，与垂直入射的Zoeppritz方程结果吻合，而DVE方程同样产生错误的振幅信息.
|
图 8 零偏移距(垂直入射)透射振幅对比(a)上行波振幅随深度变化; (b)下行波振幅随深度变化. Fig. 8 Comparison of zero offset (vertical incidence) transmission coefficient comparison (a) Amplitude of up-going wave varying with depth; (b) Amplitude of down-going wave varying with depth. |
DVE方程本质上是通过体积模量K和密度参数化表征的，习惯上会将体积模量表达为速度和密度的形式，显然这种表达方式在速度固定，密度连续情况下，波的传播特征与Zoeppritz方程结果不相符.
4.2.2 密度固定，速度连续固定密度为常数1500 kg·m-3，纵波速度沿深度递增20×Depth+2000 m·s-1.分别采用TE、DVE和TVSE对该介质进行正演模拟，图 9为传播了0.28 s的波场快照.图 10是对波场进行横向(x=1000 m)和纵向(z=1000 m)采样结果，同样由于水平方向是均匀的，图 10a(水平采样)的振幅信息一致；图 10b(垂直采样)，显然在1700 m处，TVSE的振幅要低于TE和DVE，而300 m处TVSE的振幅高于TE和DVE，这其中TE和DVE的振幅、相位信息完全一致，二者在500 m和1700 m两处振幅极大值几乎相同，大约0.41左右.
|
图 9 密度固定，速度变化介质模型波场快照对比 Fig. 9 Wavefield snapshot comparison of constant density and variable velocity model (a) TE; (b) DVE; (c) TVSE. |
|
图 10 波场一维采样对比(a)波场在x方向的采样(z=1000 m)；(b)波场在z方向的采样(x=1000 m). Fig. 10 One-dimensional sampling of wavefield (a) x direction (z=1000 m); (b) z direction (x=1000 m). |
同样，我们通过垂直入射的Zoeppritz方程来分析透射波的振幅，当密度固定情况下，垂直入射的Zoeppritz方程(公式(20))退化为
|
(22) |
根据速度连续介质模型(图 3)计算得到震源在(1000 m，1000 m)时零偏移距两侧的反射透射系数，如图 11所示，波在向上传播过程中，其透射系数大于1，而向下传播过程中透射系数始终小于1.根据透射系数计算不同深度波的最大振幅(图 12)，其振幅特征与速度固定，密度连续介质类似.对比图 10b，500 m处TVSE振幅高于TE和DVE振幅，而在1700 m处刚好相反，这也与垂直入射的Zoeppritz方程获得的透射波振幅信息相吻合.
|
图 11 零偏移距(垂直入射)反射透射系数对比(a)上行波透射系数; (b)上行波反射系数; (c)下行波透射系数; (d)下行波反射系数. Fig. 11 Comparison of zero offset (vertical incidence) reflection transmission coefficient (a) Rpp of up-going wave; (b) Tpp of up-going wave; (c) Rpp of down-going wave; (d) Tpp of down-going wave. |
|
图 12 零偏移距(垂直入射)透射振幅对比(a)上行波振幅随深度变化; (b)下行波振幅随深度变化. Fig. 12 Comparison of zero offset (vertical incidence) transmission amplitude (a) Amplitude of up-going wave varying with depth; (b) Amplitude of down-going wave varying with depth. |
前面构建速度和密度一个参数固定，另一个参数连续的模型，数值试算后并与平面纵波垂直入射下的Zoeppritz方程对比分析，验证了模拟结果的准确性.但是实际地下介质并不存在这种某一参数恒定不变的情况，因此，这里构建速度、密度同时垂向连续的介质，其参数如图 3所示，纵波速度为2×Depth+2000 m·s-1，密度为4×Depth+1000 kg·m-3.
通过式(20)计算得到垂直入射振幅随深度的变化，如图 15所示，其振幅特征与前面讨论的两种情况均类似，在两个参数共同作用，波在向上传播过程中，振幅在增大，反而亦之.这一点在波场(图 13)的垂向一维采样(图 14)中可以看出，在深度500 m处，TVSE振幅极值大约为0.8，高于TE和DVE振幅，而1700 m处，TVSE振幅小于TE和DVE.在这个模型中，DVE依然存在振幅异常情况，这一点在4.2.1节中已经说明，在此不做分析.
|
图 13 速度、密度连续介质模型波场快照对比 Fig. 13 Wavefield snapshot comparison of variable velocity and density model (a) TE; (b) DVE; (c) TVSE. |
|
图 14 波场一维采样对比(a)波场在x方向的采样(z=1000 m)；(b)波场在z方向的采样(x=1000 m). Fig. 14 One-dimensional sampling of wavefield (a) x direction (z=1000 m); (b) z direction (x=1000 m). |
|
图 15 零偏移距(垂直入射)透射振幅对比(a)上行波振幅随深度变化; (b)下行波振幅随深度变化. Fig. 15 Comparison of Zero offset (vertical incidence) transmission amplitude (a) Amplitude of up-going wave varying with depth; (b) Amplitude of down-going wave varying with depth. |
前面讨论的均匀介质和连续介质均属于特殊的非均匀介质，而介质的非均匀性随处可见.介质中弹性参数可以连续，也可以间断，常见的层状介质在弹性分界面处是不连续的，即不可导.本文的TVSE存在速度、密度的空间导数，因此在数值模拟过程中波遇到弹性间断面，势必存在一定的问题.
建立一个简单双层模型，网格点数为200×200，空间步长5 m，时间采样间隔0.5 ms，震源采用主频20 Hz的雷克子波，坐标(500 m, 0 m).模型分为三组进行测试:
(1) 固定密度1000 kg·m-3，纵波速度上层为2000 m·s-1，下层为3000 m·s-1；
(2) 固定速度2000 m·s-1，密度上层为1000 kg·m-3；下层为2000 kg·m-3；
(3) 纵波速度上层为2000 m·s-1，下层为3000 m·s-1，密度上层为1000 kg·m-3；下层为2000 kg·m-3.
图 16是TE、DVE和TVSE对第一组模型数值模拟至0.4 s时刻的波场快照，TVSE在弹性界面的反射波存在相位翻转现象.图 17是第二组模型模拟结果，因为TE不含密度项，所以不产生反射，而TVSE模拟结果依然存在反射波的相位翻转.图 18是第三组模型模拟结果，TVSE模拟的反射波能量明显高于透射波能量.对比前两组模型，这一现象是由于密度界面引起的，即对于TVSE数值模拟，当介质的速度变化不剧烈的时候(均匀或者连续)，强烈的速度、密度分界面(或者速度、密度扰动)对反射影响更大.
|
图 16 第一组模型数值模拟波场快照 Fig. 16 Wavefield snapshots simulated by first set of models (a) TE; (b) DVE; (c) TVSE. |
|
图 17 第二组模型数值模拟波场快照 Fig. 17 Wavefield snapshot simulated by second set of models (a) TE; (b) DVE; (c) TVSE. |
|
图 18 第三组模型数值模拟波场快照 Fig. 18 Wavefield snapshot simulated by third set of models (a) TE; (b) DVE; (c) TVSE. |
当速度、密度变化不剧烈(或扰动较小)的时候，层状介质可以近似为连续介质(比如4.2节的连续介质)，用TVSE数值模拟是可行的，但是当速度、密度变化剧烈，由于TVSE方程中存在弹性参数的直接求导项，会造成透射能量异常低值，图 19即是固定密度，速度上层速度为2000 m·s-1，下层分别为3000 m·s-1、4000 m·s-1和5000 m·s-1三组模型选取TVSE模拟得到的0.4 s时刻波场图，尽管界面差异增大会造成透射能量减小，但是由于上下界面的速度差增大，透射能量产生明显损失，非常微弱，显然对于存在强反射层的介质，TVSE对于强反射层下的反射波的“反馈”是不理想的，这与方程(9)中对速度和密度进行了空间求导有关，4.2节测试表明TVSE在连续介质中的波形更准确.
|
图 19 不同层状模型TVSE模拟波场快照 Fig. 19 Wavefield snapshots simulated by different layered models using TVSE |
为了验证本文方法对于大型模型数值模拟的稳定性以及PML边界吸收效果，选取Marmoursi-2模型，模型网格数为2401×1201，网格间距5 m，时间采样间隔0.5 ms，震源采用主频为20 Hz的雷克子波，位于坐标(6005 m, 10 m)处，检波器布于水面.具体速度、密度数值见图 20.图 21是通过Seismic Unix包中函数smooth2横纵坐标均平滑20次后得到Marmoursi-2连续介质模型.
|
图 20 Marmousi-2模型(a)速度；(b)密度. Fig. 20 Marmousi-2 model (a) Velocity; (b) Density. |
|
图 21 armousi-2连续介质模型(a)速度；(b)密度. Fig. 21 Marmousi-2 continuous model (a) Velocity; (b) Density. |
TE不考虑密度变化，Marmousi-2模型选用DVE和TVSE进行数值模拟.图 22是分别使用两个方程模拟到2.5 s时刻波场快照，因为模型存在比较强烈的弹性间断面，因此TVSE模拟的波场在深度4000 m至6000 m波的能量更弱，相反DVE对于深层介质的物性差振幅响应更明显，这也印证了4.3节的结论，因此，TVSE不适合介质存在强弹性间断面的模拟.
|
图 22 Marmousi-2模型模拟波场快照 Fig. 22 Wavefield snapshots simulated by Marmousi-2 model (a) DVE; (b) TVSE. |
对连续介质的Marmousi-2模型采用TVSE进行正演模拟，如图 23所示，在介质连续情况下，深层的波场传播稳定，深层透射能量充分，能够准确对地下参数变化产生响应.TVSE模拟过程中才用PML吸收边界的差分方程(式(15))，从图 23可以看出波传播到上边界时有良好的吸收效果.因此可以证明利用此PML条件，在进行连续介质数值模拟的过程中可以取得理想的效果.
|
图 23 Marmousi-2连续介质模型TVSE模拟波场快照 Fig. 23 TVSE wavefield snapshot simulated by Marmousi-2 continuous model |
介质的非均匀是一个相对的概念，本质在于地震波的一个波长能否识别介质的单元尺度，用常用的模型正演模拟过程中，当子波主频较高的时候，地震波的波长小于一个层的厚度，那么非均匀性只存在于弹性间断面，而在地层中的介质是均匀的，因此用常规的均匀假设条件下的标量波方程进行数值模拟是可行的，其非均匀性只存在于界面，不会影响地震波在地层内的传播，界面的非均匀主要影响不明显，这也是为什么可以用生产中(例如，全波形反演、逆时偏移、波动方程走时反演等)可以用均匀假设下的标量波方程进行正演.而当子波主频较低，即一个地震波长大于或者接近地层的厚度(如图 3所示模型)，常规方程不能够真实反演地震波在地下传播的振幅变化，对比图 5、图 10和图 14可知，常密度标量波方程在非均匀介质中传播，其振幅基本相等，而变密度标量波方程能够通过振幅变化指示介质的非均匀性，但是其振幅变化规律与物理意义不符.综上所述，TVSE对于地下介质尺度较小，或者是子波主频较低的情况下正演能够得到更准确反映地下介质的非均匀性的波形特征，这对于反问题中参数模型长波长信息的恢复提供了基础.
5 结论本文借鉴Alkhalifah声学近似思想，充分考虑非均匀介质弹性参数的空变特征，从非均匀介质的弹性波动方程出发，利用本征值问题求解方法建立标量波频率-波数域传播算子，从而建立描述纵波传播的双变参数标量波方程(TVSE)，并推导了PML边界下波动方程差分格式.结合数值示例得出以下几点结论:
(1) 本文通过建立二维双变参数标量波方程来模拟非均匀介质，该方法可以推广到高阶有限差分建模以及三维，采用PML边界处理人工反射得到良好的效果，通过大型Marmousi-2模型证明算法的稳定.
(2) 存在剧烈弹性间断的非均匀模型不宜选取TVSE进行数值模拟，TVSE更适用于连续介质或者近似连续介质(弹性间断相对平缓)，并且相比于TE和DVE，其模拟波场特征(振幅、相位信息)更准确，DVE在传播过程中表现出相反的振幅特征，这与实际并不相符.这一点通过与Zoeppritz方程的讨论中得到证明.
(3) 本文方程可以应用于地球物理反问题，如全波形反演，TVSE能够解决大尺度介质参数的反演，其综合考虑弹性参数的空变特征，有利于更准确解决大尺度背景描述的问题.
致谢 感谢评审专家为本文完善提出的宝贵意见.
Aki K, Richards P G. 2002. Quantitative Seismology. .
|
Alkhalifah T. 2000. An acoustic wave equation for anisotropic media. Geophysics, 65(4): 1239-1250. DOI:10.1190/1.1444815 |
Alkhalifah T. 2003. An acoustic wave equation for orthorhombic anisotropy. Geophysics, 68(4): 1169-1172. DOI:10.1190/1.1598109 |
Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2): 185-200. DOI:10.1006/jcph.1994.1159 |
Boore D M. 1972. Finite difference methods for seismic wave propagation in heterogeneous materials. Methods in Computational Physics:Advances in Research and Applications, 11: 1-37. |
Di Bartolo L, Dors C, Mansur W J. 2012. A new family of finite-difference schemes to solve the heterogeneous acoustic wave equation. Geophysics, 77(5): T187-T199. DOI:10.1190/geo2011-0345.1 |
Dong L G, Ma Z T, Cao J Z. 2000. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation. Chinese Journal of Geophysics (in Chinese), 43(6): 856-864. |
Duan P R, Li Q Y, Zhao Z Q, et al. 2019. High-order finite-difference method based on equivalent staggered grid scheme for scalar wavefield simulation. Progress in Geophysics (in Chinese), 34(3): 1032-1040. DOI:10.6038/pg2019CC0064 |
Higdon R L. 1987. Numerical absorbing boundary conditions for the wave equation. Mathematics of Computation, 49(179): 65-90. DOI:10.1090/S0025-5718-1987-0890254-1 |
Hook J F. 1961. Separation of the vector wave equation of elasticity for certain types of inhomogeneous, Isotropic Media. The Journal of the Acoustical Society of America, 33(3): 302. |
Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms:A finite-difference approach. Geophysics, 41(1): 2-27. DOI:10.1190/1.1440605 |
Levander A R. 1988. Fourth-order finite-difference P-SV seismograms. Geophysics, 53(11): 1425-1436. DOI:10.1190/1.1442422 |
Li Q Y, Wu G C, Liang Z Y. 2018. Time domain high-order pseudo spectral method based on PML boundary for elastic wavefield simulation. Progress in Geophysics (in Chinese), 33(1): 228-235. DOI:10.6038/pg2018BB0014 |
Li Z C, Tian K, Huang J P, et al. 2013. An unsplit implementation of the multiaxial perfectly matched layer. Progress in Geophysics (in Chinese), 28(6): 2984-2992. |
Liu Q H, Tao J. 1997. The perfectly matched layer for acoustic waves in absorptive media. The Journal of the Acoustical Society of America, 102(4): 2072-2082. DOI:10.1121/1.419657 |
Liu Y, Li C C, Mou Y G. 1998. Finite-difference numerical modeling of any even-order accuracy. Oil Geophysical Prospecting, 33(1): 1-10. |
Liu Y S, Liu S L, Zhang M G, et al. 2012. An improved perfectly matched layer absorbing boundary condition for second order elastic wave equation. Progress in Geophysics (in Chinese), 27(5): 283-292. |
Madariaga R. 1976. Dynamics of an expanding circular fault. Bulletin of the Seismological Society of America, 66(3): 639-666. |
Richards P G. 1974. Weakly coupled potentials for high-frequency elastic waves in continuously stratified media. Chinese Medical Sciences Journal, 64(2): 90-90. |
Robertsson J O A, Blanch J O, Symes W W, et al. 1994a. Galerkin-wavelet modeling of wave propagation:Optimal finite-difference stencil design. Mathematical and Computer Modelling, 19(1): 31-38. |
Robertsson J O A, Blanch J O, Symes W W. 1994b. Viscoelastic finite-difference modeling. Geophysics, 59(9): 1444-1456. DOI:10.1190/1.1443701 |
Sun C Y. 2007. Theory and Method of Seismic Waves. .
|
Virieux J. 1984. SH-wave propagation in heterogeneous media; Velocity-stress finite-difference method. Geophysics, 49(11): 1933-1942. DOI:10.1190/1.1441605 |
Virieux J. 1986. P-SV wave propagation in heterogeneous media; velocity-stress finite-difference method. Geophysics, 51(4): 889-901. DOI:10.1190/1.1442147 |
Wang S D. 2003. Absorbing boundary condition for acoustic wave equation by perfectly matched layer. Oil Geophysical Prospecting, 38(1): 31-34. |
Wang T, Tang X M. 2003. Finite-difference modeling of elastic wave propagation:A nonsplitting perfectly matched layer approach. Geophysics, 68(5): 1749-1755. DOI:10.1190/1.1620648 |
Wang Y G, Xin W J, Xie W X., et al. 2007. Study of absorbing boundary condition by perfectly matched layer. Journal of China University of Petroleum (edition of natural science), 31(1): 19-24. |
Yin W, Yin X Y, Wu G C, et al. 2006. The method of finite difference of high precision elastic wave equations in the frequency domain and wave_field simulation. Chinese Journal of Geophysics (in Chinese), 49(2): 561-568. |
Zeng Y Q, Liu Q H. 2001. A staggered-grid finite-difference method with perfectly matched layers for poroelastic wave equations. The Journal of the Acoustical Society of America, 109(6): 2571-2580. DOI:10.1121/1.1369783 |
Zhang Y, Wu G C. 2013. Methods of removing pseudo SV-wave artifacts in TTI media qP-wave reverse-time migration. Chinese Journal of Geophysics (in Chinese), 56(6): 2065-2076. DOI:10.6038/cjg20130627 |
董良国, 马在田, 曹景忠. 2000. 一阶弹性波方程交错网格高阶差分解法稳定性研究. 地球物理学报, 43(6): 856-864. DOI:10.3321/j.issn:0001-5733.2000.06.015 |
段沛然, 李青阳, 赵志强, 等. 2019. 等效交错网格高阶有限差分法标量波波场模拟. 地球物理学进展, 34(3): 1032-1040. DOI:10.6038/pg2019CC0064 |
李青阳, 吴国忱, 梁展源. 2018. 基于PML边界的时间高阶伪谱法弹性波场模拟. 地球物理学进展, 33(1): 228-235. DOI:10.6038/pg2018BB0014 |
李振春, 田坤, 黄建平, 等. 2013. 多轴完全匹配层的非分裂实现. 地球物理学进展, 28(6): 2984-2992. |
刘洋, 李承楚, 牟永光. 1998. 任意偶数阶精度有限差分法数值模拟. 石油地球物理勘探, 33(1): 1-10. DOI:10.3321/j.issn:1000-7210.1998.01.001 |
刘有山, 刘少林, 张美根, 等. 2012. 一种改进的二阶弹性波动方程的最佳匹配层吸收边界条件. 地球物理学进展, 27(5): 283-292. |
孙成禹. 2007. 地震波理论与方法. .
|
王守东. 2003. 声波方程完全匹配层吸收边界. 石油地球物理勘探, 38(1): 31-34. DOI:10.3321/j.issn:1000-7210.2003.01.007 |
王永刚, 邢文军, 谢万学, 等. 2007. 完全匹配层吸收边界条件的研究. 中国石油大学学报(自然科学版), 31(1): 19-24. DOI:10.3321/j.issn:1000-5870.2007.01.004 |
殷文, 印兴耀, 吴国忱, 等. 2006. 高精度频率域弹性波方程有限差分方法及波场模拟. 地球物理学报, 49(2): 561-568. DOI:10.3321/j.issn:0001-5733.2006.02.032 |
张岩, 吴国忱. 2013. TTI介质qP波逆时偏移中伪横波噪声压制方法. 地球物理学报, 56(6): 2065-2076. DOI:10.6038/cjg20130627 |