2016, Vol. 59
2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
地震波数值模拟对于定量地分析研究复杂地下介质中的地震波成像非常重要.传统的地震波模拟基于弹性模型假设，可以满足一般的地震勘探精度要求，得到很好的数值模拟结果(Gazdag，1981).但近年来勘探技术的发展对数值模拟的精度提出了更高的要求，需要考虑地下介质黏性对数值模拟结果的影响.同时，计算机技术的飞速发展为高精度黏弹性数值计算提供了可能.前人关于黏弹性数值模拟的研究分为频率域和时间域两类.频率域数值模拟可以适用任意的衰减准则描述介质的衰减特性，得到广泛的应用(Brossier et al.，2010；廖建平等，2011).但该方法在处理大型数据时可能面临着内存需求过大的问题.而时间域数值模拟时，面临着黏性应力-应变之间的褶积运算会增加数值计算难度的问题.为改善时间域数值计算的困难，前人们提出了很多描述黏性衰减的策略，如引入流变模型或改变波动方程表示等.
目前描述衰减问题常用的流变模型，是通过不同形式的串并联线性弹簧体和阻尼器得到的一系列简单的线性黏弹性固体模型，如Kelvin-Voigt固体模型，Maxwell固体模型等(Blake et al.，1982).这些流变模型基于相应的物理模型，描述介质的蠕变和松弛响应，及介质的应力应变关系.但也存在一定的局限性.首先，和实验所测材料的应力应变关系相比，Maxwell模型和Kelvin-Voigt模型都只能部分(蠕变或松弛过程)和实验结果相吻合.其次，该类方法没有考虑模型的Q(ω)函数是否符合衡Q要求，不能得到地震波场的真实振幅，与实际地球介质属性不相符.
另一种描述衰减的策略是对弹性波方程做变动(吴玉等，2014)，以克服传统计算中的内存需求过大的问题.Zhu和Harris(2013，2014)、Zhu(2014)基于Kjartansson(1979)的衡Q模型和Chen和Holm(2004)的分数式Laplace算子方法，通过在弹性波波动方程上加上黏性近似项，实现黏弹性介质的时间域数值模拟.该方法可以描述近乎衡Q的衰减频散响应，但所加的黏性项是近似项，容易造成算法的不稳定，导致计算精度受限.
还有一种方法是Liu等(1976)提出的用松弛机制谱的概念定义本构关系.Liu认为黏性Q在地震波频带范围内近似为常数，可由多个松弛机制的组合近似表示.Day和Minster(1984)更进一步提出，借助Padé近似将时间域的褶积算子化为一系列微分算子，进而由一个n阶有理函数近似表示黏弹性模量，并解析地求出其中的参数.Emmerich和Korn(1987)在此基础上，赋予了该黏弹性模量相应的物理意义--可由广义的Maxwell流变模型表示.该广义模型由一个弹簧体和n个经典的Maxwell体并联组成，n个经典的Maxwell体分别对应有理式中的n个叠加项.该类方法可以近似地满足衡Q模型要求，但需要寻得尽可能低阶的n，以减少计算时对内部变量的存储需求，保证计算效率.
本文描述衰减问题的方法就是基于上述策略展开的.为了快速地搜索到低阶的又能满足衡Q要求的模型参数，本文以无穷范数构建目标函数，结合全局搜索的模拟退火方法寻求最优解.最终确定，当阶数n=2或3时，我们就能找到满足要求的衡Q模型，从而极大地减少了数值计算量和所需的存储空间.另一方面，为进一步提高计算效率和精度，本文结合了非规则非结构网格的格子法(Zhang and Liu，1999)作为黏弹性介质的数值模拟算法.相比于常见的有限差分法(Robertsson et al.，1994；Saenger and Bohlen，2004；Operto et al.，2007；王忠仁等，2014)、有限元法(Marfurt，1984)和伪谱法(Carcione，2010)等，格子法可以自然地满足自由边界条件，更精细地刻画复杂介质的速度变化及Q值变化.为验证本文算法的精度和有效性，文中给出了简单模型算例中数值解和解析解的比较，测试了算法在气云构造和Marmousi模型中的模拟效果，都取得了很好的结果.
2 非规则、非结构网格声波模拟的格子法本文数值解法是基于声波模拟的格子法(Zhang and Liu，1999，2002)展开工作.格子法采用非规则、非结构化网格离散速度模型，既可以精细地刻画复杂边界条件、速度及Q值的变化，也具有与有限差分法相当的计算效率.
这里简单地介绍一下声波模拟下格子法的基本原理.考虑非均匀介质中的无源声波方程(徐义和张剑锋，2008)：
|
(1) |
式中 p(x，z，t) 表示波场值， ρ(x，z)，v(x，z) 分别表示介质的密度和速度， Δ 表示梯度.
格子法的理论核心是在非规则网格下，给出式(1)的基于积分平衡的微分方程弱形式.如图 1所示，在节点 k 的邻域，定义虚线轮廓线1-2-3-…-12-1所围区域为Ω，简记该轮廓线为Ω，对式(1)两边作面积分，并对等式右边应用Green定理得到
|
(2) |
式中m为围绕节点k的三角形单元数目，Skl对应节点k周围第l个三角形单元中的虚线段，α与β分别是边界Ω外法线沿x与z轴的方向余弦.利用动力学计算的集中质量原理及三角形线性插值方法可得到式(2)的离散形式：
|
(3) |
式中 Mk 是节点 k 各邻域三角形单元中的 
|
(4) |
式中，Δ为三角形ijk的面积，下标r代表三角形中的各顶点， (bk)l，(ck)l 分别表示围绕节点k的第l个单元对应的几何参数，以图 1中的三角形单元ijk为例，可简单地表述为bk=(zj-zi)/2，ck=(xj-xi)/2.
|
图 1 非规则三角网格局部示意图 Fig. 1 Local irregular triangle grid |
最后，利用中心差分格式离散方程式(3)左端中对时间的二阶偏导数，我们就可以实现声波波场的迭代更新.
3 体积模量的有理函数展开和衡Q模型 3.1 体积模量的有理函数展开黏弹性介质在频率域中有如下应力应变关系：
|
(5) |
其中， M(ω) 是与频率有关的复体积模量， σ(ω) 是应力， ε(ω) 为应变，经由Fourier变换，得到如下褶积关系：
|
(6) |
这意味着在时间域计算时，需要存储历史时刻的应变值，会导致极大的计算负担，是数值计算很不乐见的.所以，我们需要考虑从其他方面对其做简化近似处理.
已知实验测得实际介质的应力松弛响应关系为：在常应变作用下，应力呈指数衰减变化，即实际黏性介质的松弛函数可用一个负指数衰减的函数表示.不同黏性介质的松弛响应对应不同的幅值及衰减系数表示.我们总可以找到这样的一组基函数--由一系列幅值和衰减系数组合表示任意的衰减函数，并用其表示介质的松弛响应.我们定义未松弛模量MU及松弛模量MR，分别表示介质衰减前、后的体积模量，并基于此构建如下负指数衰减的松弛函数：
|
(7) |
式中， H(t) 为Heaviside阶梯函数，权重 aj满足
至此，可将体积模量表示为
|
(8) |
我们通过Fourier变换，将其变换到频率域 (e-ωjt→FFT→1/(iω+ωj))， 有
|
(9) |
黏弹性介质中，频率描述的体积模量可以由n个松弛频率的有理式函数展开表示.
3.2 衡Q模型实际地球介质中，黏性 Q(ω) 是随频率变化的函数.但是在地震频带范围内Q可近似为常数，或随频率微弱变化.地球物理勘探中的信号频带一般在3~100 Hz内，我们可以采用在该频带范围内满足衡Q要求的模型描述黏弹性介质的高精度勘探问题.在介质黏性吸收满足衡Q要求时，我们就可以确定体积模量Mn(ω)中有理式函数的系数.
针对式(9)中有理函数表示的体积模量Mn(ω)，品质因子Q可表示为
|
(10) |
为寻求满足衡Q模型要求的最优解，我们用Q0表示目标Q值，以无穷范数构建目标函数，得到下列表示：
|
(11) |
一般低阶(n≤8)时有2n+1个模型参数，还比较少，我们可以采用蒙特卡洛、模拟退火、遗传算法等全局搜索方法.本文采用的是优化的模拟退火方法(Goffe et al.，1994)，可实现全局的快速搜索，求得最优解.然后将最优解对应的参数结果列表记录下来，便于以后模拟计算时查表使用.最终得到的结果和目标值的拟合效果，可以参见图 2所示，图中实线表示目标值Q0=20.0，虚线为Q-1(ω)拟合结果.图 2中，分别对应阶数n=2，n=3时Q(ω)拟合得到的曲线结果和目标值的比较.可以看出，在频带3~100 Hz内，n=2时即可得到不错的拟合效果；n=3时，Q(ω)曲线可以更加贴合目标值，拟合效果更优.本文数值模拟的简单算例，采用的就是阶数n=3时的参数结果.
|
图 2 阶数(a)n=2，(b)n=3时，Q-1=0.05的目标值(实线)和Q-1(ω)曲线拟合结果(虚线) Fig. 2 Target values(solid lines)and curve fitting results(dash lines)of the Q-1(ω)， with different order:(a)n=2，(b)n=3 |
假设 δM/MU<<1， 则式(10)可近似为
|
(12) |
观察式(11)， 
将式(9)表示的体积模量 Mn(ω) 代入式(5)中，并引入中间变量 φj(ω)， 得到
|
(13) |
式中，中间变量 φj(ω) 满足关系式
|
(14) |
对式(13)和式(14)做Fourier反变换，得到时间域黏弹性介质的波动方程：
|
(15) |
其中，
由此，我们得到如下非均匀黏弹性介质中的声波方程组：
|
(16) |
|
(17) |
|
(18) |
式中的系数 aj，ωj，及δM 可由本文3.2节中的计算列表查询得到.当 δM=0 时，介质不含黏性衰减，则式(16)-(18)可合并成式(1)所示的声波方程.
为便于模拟复杂介质下波传播情况，本文采用格子法中的三角网格单元，即非规则、非结构网格离散速度模型，并延用格子法原理到黏声波方程组(式(16)-(18))，最终得到非规则网格的黏声波数值模拟解法.
仍以图 1为例，延用声波格子法的积分近似推导，在虚线包围的k节点邻域对式(17)做面积分，并利用Green定理，得到
|
(19) |
应用前面提到的集中质量原理和三角差分算子表示，得到
|
(20) |
式中， ck，bk 和式(4)中三角单元的几何参数有相同的意义，Δ为三角单元的面积.最终，对式(16)、(20)和(18)使用时间域的中心差分近似，可得
|
(21) |
|
(22) |
|
(23) |
式中上标 q 代表时间，下标 l 表示节点 k 周围第 l 个三角形单元.式(21)-(23)构成了非规则网格黏声波方法的应变和压强更新基础.由 q 时刻的压强 p 代入式(22)可求得 q+0.5 时刻各节点处的形变 ε ，再结合式(23)求得 q+0.5 时刻各节点处的中间变量 φj， 最后由式(21)求得 q+1 时刻的压强波场 p， 完成对压强 p 在时间上的更新.
相比于声波方程，该黏声波方程需要额外计算n个一阶时间导方程，储存n+1个中间变量 ε(x，z，t)，φj(x，z，t)(j=1，…，n). 为保证计算效率，我们需要使用尽可能低阶的阶数n，同时又要满足地震波黏性吸收的衡Q要求.本文3.2节已经确定在n=2或3时找到很好的满足衡Q要求的参数.另一方面，格子法中的非规则网格的使用，允许在数值计算时，黏弹性介质区域参与黏声波方程的计算，弹性介质仍按声波方程模拟计算，也可大大缩减计算存储需求，提高计算效率.
数值模拟时四周的边界条件，延用徐义和张剑锋(2008)的非规则PML吸收方法.在吸收带内，本文采用随边界位置变化的局部坐标系，见图 1中的 (w，v) 坐标系，结合引入的中间变量φj，构建基于局部坐标系的黏声波PML分裂方程，这里不再详述.该方法可以很好地吸收或透射不同方向的入射波，避免边界反射问题.
根据Zhang和Liu(1999)对弹性波格子法的数值分析，本文讨论以等边三角形单元为剖分网格时黏声波数值计算的稳定性条件及频散关系问题.在空间采样方面，根据不同介质中的最小速度cmin，结合频散关系公式，确定网格尺寸 h≤cmin/nf，f 为震源主频， n 与具体的算法有关，格子法中取20以上即可.时间的采样间隔由 Δt≤h/
为验证本文提出的非规则网格黏声波数值模拟解法的正确性和精度，首先将其应用于均匀速度模型(对此模型可容易得到黏声波传播的解析解)，并将其计算结果与解析解对比分析.然后，给出一个强Q衰减扰动情况下的数值模拟解算例，以验证该计算方法的稳定性.最后应用到Marmousi模型考察该方法对复杂横向速度、Q值变化情况的适应能力.
5.1 均匀介质模型为了评估本文解法的精度，我们分别从该算法对炮检距和Q值的响应这两个方面讨论、比较数值模拟得到的数值解与解析解.用于对比的解析解，是根据对应原理(Carcione et al.，1988)，在频率域中用其中的黏弹性参数代替弹性解析解中的弹性参数，再反变换到时间域得到.数值解方面，我们设计了一个横向4 km，纵向2 km的各向同性均匀黏弹性介质模型，密度、速度分别为ρ=2.0 g·cm-3，v=2000 m·s-1，炮点位于模型的中心位置.文中我们采用主频为30Hz的Ricker子波作为震源，时间采样步长为0.33 ms，利用非规则网格离散该速度模型共得到4527531个网格单元参与数值计算.
|
图 3 不同炮检距时，数值解(虚线)和解析解(实线)的比较.检波器距离震源距离依次为(a)600，(b)1600，(c)3200 m.二者结果吻合较好 Fig. 3 Comparison of numerical solutions(dash line)with the analytical solutions(solid line). The receivers are at(a)600，(b)1600，(c)3200 m from the source. The two results are matched with each other well |
本文首先研究Q=20时，数值解对炮检距的敏感度--能否在远炮检距处保证相应的计算精度.我们分别在距震源600 m，1600 m，3200 m处，放置一个检波器记录波场值，并与相同位置的解析解对比，最终结果如图 3所示(图中实线表示解析解，虚线表示数值解).从图 3中可以看出，该算法的数值解与解析解都能较好的吻合.为进一步定量地评估本文算法的精度，本文利用数值解与解析解的均方根误差来衡量，并定义均方根误差E：
|
(24) |
式中N表示地震道采样的时间数， djo，dja 分别为j采样点时得到的数值解和解析解.经过式(24)的计算，我们在炮检距600 m，1600 m，3200 m处算得的均方根误差依此为6.9×10-3，1.3×10-2，8.1×10-3.进一步说明，炮检距越近，数值解与解析解吻合的越好；远炮检距处的数值解与解析解仍可以很好地吻合.本文算法对远距离波场传播仍能满足较高的精度要求.
5.1.2 不同Q值时在这部分，我们讨论数值模拟结果对Q值的敏感度，尤其在低Q时，模拟结果的精度问题.仍采用上面的均匀速度模型和震源，在点(2.6 km，1.0 km)处对应的数值解(虚线)和解析解(实线)的比较.它们的结果都能很好的吻合
放置一检波器，记录四个不同Q值时波传播情况.图 4为不同Q值时，非规则网格数值模拟结果在T=200 ms时刻的波场快照图：从左到右，从上到下，Q值分别为1000，100，30，10.从图 4中可以看到，Q值的减小伴随着波场幅值更强的衰减及相位延迟.为了更详细地了解模拟结果对Q的敏感度，我们取不同Q值在同一检波器处的波场记录，并与相应的解析解对比分析，如图 5所示，图中实线表示解析解，虚线为数值解.我们仍利用式(24)得到相应的均方根误差分别为3.8×10-3，3.9×10-3，5.4×10-3，1.3×10-2，说明Q值很小时，本文解法仍可得到稳定的高精度的计算结果.
以上两个数值算例结果都验证了本文数值解法在远炮检距、低Q处均具有稳定的高精度解.
|
图 4 均匀衰减介质中，T=200 ms时刻的波场快照图，图中的四个部分分别对应不同的Q值：(a)1000，(b)100，(c)30，(d)10.震源位于模型的中间位置.随着Q变小，地震波波场有更强的幅值衰减和相位延迟 Fig. 4 Four snapshot parts corresponding to four Q values:(a)1000，(b)100，(c)30，(d)10. A point source located at the center of the model. Snapshots are recorded at 200 ms in homogeneous attenuating media.The amplitude of wavefields decreases rapidly with the Q values decline |
|
图 5 震源位于点(2.0，1.0)km，检波器位于(2.6，1.0)km处，不同的Q值：(a)1000，(b)100，(c)30，(d)10时， Fig. 5 Comparison between simulated solutions(dash lines)and analytical solutions(solid lines)corresponding to four Q values: (a)1000，(b)100，(c)30，and(d)10 at point(2.6，1.0)km from a source(2.0，1.0)km . The results show an excellent agreement |
具有强衰减(即低Q)特征的气云构造常给地震成像带来很大的困扰，应用本文数值算法到气云模型，可以验证算法的稳定性.模型设计如下：模型大小为2000 m×1000 m，地下介质为速度v=2000 m·s-1的黏弹性介质.在点(600 m，600 m)处有一半径为50 m的圆形气云(图 6中的白色圆圈)，气云内Q=10，气云外为Q=∞(即弹性介质).数值模拟时，取主频为30 Hz的Ricker子波作为爆炸源置于点(1000 m，10 m)处(图 6中的星号)，采样时间步长为0.333 ms.图 6a和6b分别对应传播时间为370 ms和450 ms时的地震波场快照.从图中可以看出，在气云界面有弱的反射波产生，透过气云的波场存在幅值的衰减和相位的略微滞后现象.说明本文方法在强Q变化模型中可得到稳定的计算结果.另一方面，鉴于非规则、非结构网格的优势，计算时气云内使用黏声波方程，气云外的区域按弹性波方程计算，可以有效节省黏声波计算时的存储空间，提高计算效率，是其他规则网格方法所不具备的.
|
图 6 (a) 370 ms，(b)450 ms时刻的波场快照. 在图中白色圆圈代表的气云周围有反射产生 Fig. 6 Snapshots of pressures at(a)370 ms and(b)450 ms propagation time. There are reflections at interface of the gas-pocket which is represented by the white circle in the figures |
为进一步验证本文方法对复杂横向速度及Q值变化的适应能力，将其应用于SEG 的Marmousi模型，该模型存在断层、背斜、尖灭等复杂构造，同时具有较强的横向介质变化.图 7显示了其速度模型和根据速度计算得到的Q值模型.其中，Q值的计算参照经验公式(Li，1993)： Q=3.516×v2.2×10-6， 得到的Q值范围为34.2~595.5.震源置于点(4.8 km，0.009 km)处，为主频30 Hz的Ricker子波，采样时间间隔为0.267 ms，剖分模型得到11931119个三角形网格单元.为简便计算，本文取Q=50时对应的三阶Q值模型参数为参考值，利用式(11)求取其余Q值时的参数，即对应的黏声波方程系数.应用本文算法模拟Marmousi模型的黏声波传播，得到传播时间为1.2 s时的波场快照，如图 8a所示.为对比验证黏性的引入对波场的影响，本文也模拟了同等情况下，声波在Marmousi模型中的波场传播.取同在水平剖面0，4.6 km内，传播时刻均为1.2 s时的声波模拟结果和黏声波模拟结果的波场快照结果对称地显示在同一张图里，比较分析其传播情况(图 8b).在图 8b所示的图中，左半部分为声波模拟结果，右半部分为其同一位置同一时刻的黏声波模拟结果.图 9为二者的频谱图比较，图中实线为声波频谱结果，虚线为黏声波频谱.对比发现：首先，黏声波模拟和声波模拟的波场形态基本是一致的，证明了本文黏声波算法的正确性.其次，观察二者的走时变化，可以看到黏性Q的引入，使得波传播速度变慢，波场传播存在相位延迟现象.然后，由二者幅值的比较可以看出，黏性Q会导致地震波幅值的衰减.最后，频谱分析说明黏性介质使得地震波能量的峰值由主频30 Hz向低频方向移动，对地震子波中的高频成分吸收强烈，进一步验证了本文黏声波算法的正确性与合理性.
|
图 7 (a) Marmousi速度模型;(b)根据速度计算得到的Q值模型 Fig. 7 (a) Marmousi velocity model；(b)Q model derived from the velocity model |
|
图 8 传播时刻T=1.2 s时的波场快照 (a)黏声波结果;(b)声波结果(左半部分)与黏声波结果(右半部分)在同一时刻的对比.黏声波波场的幅值有明显的衰减. Fig. 8 The snapshots of pressure at T=1.2 s (a)The viscoacoustic wave;(b)The comparisons of the acoustic wave(the left part)and the viscoacoustic wave(the right part) at the same time. The amplitude of viscoacoustic wave field is significant reduced by the attenuation. |
|
图 9 Marmousi模型模拟结果的声波频谱与黏声波频谱对比图.图中实线代表声波频谱，虚线代表黏声波频谱.黏性介质下地震波能量向低频方向移动，高频成分被强烈吸收 Fig. 9 The comparison of acoustic spectrum and viscoacoustic spectrum derived from the Marmousi Model. The solid line represents the acoustic spectrum，and the dash line represents the viscoacoustic spectrum. In the attenuation medium，the seismic energy moves toward to the low frequencies，and strongly absorbs the high frequency components |
本文发展了基于非规则化网格的时间域黏声波数值模拟算法.文中采用有理式基函数描述频率相关的体积模量，通过无穷范数构建目标函数，并结合模拟退火方法，实现了满足衡Q要求的低阶基函数参数的快速搜索，降低了计算中的内存需求.同时，由引入n个关于中间变量的一阶偏微分方程，替代常规时间域黏声波数值模拟中的褶积运算，极大地缓和了计算存储问题，提高了计算效率.低阶有理式表示的体积模量和格子法的结合，使得本文的数值模拟算法，既可精细地刻画复杂边界条件，也能兼顾计算精度和效率.简单模型中和理论解的分析比较证明了方法的精度及稳定性.二维Marmousi模型算例表明，本文算法对复杂介质模型具有很好适用性.本文方法也可延用到黏弹性波模拟及成像，进一步提高地震勘探精度，有很好的应用前景.
| Blake R J, Bond L J, Downie A L. 1982. Advances in numerical studies of elastic wave propagation and scattering.//Reviewof Progress in Quantitative Nondestructive Evaluation.New York:Plenum Press, 157-166. | |
| Brossier R, Operto S, Virieux J, et al. 2010. Frequency-domain numerical modelling of visco-acousticwaveswith finite-difference and finite-element discontinuous Galerkinmethods.//Dissanayake DW. Acoustic Waves.SCIYO. | |
| Carcione J M, Kosloff D, Kosloff R. 1988. Wave propagation simulation in a linear viscoelastic medium. Geophysical Journal International , 95 (3) : 597-611. DOI:10.1111/j.1365-246X.1988.tb06706.x | |
| Carcione J M. 2010. A generalization of the Fourier pseudo spectral method. Geophysics , 75 (6) : A53-A56. DOI:10.1190/1.3509472 | |
| Chen W, Holm S. 2004. Fractional Laplacian time-space models for linear and nonlinear lossymedia exhibiting arbitrary frequency power-law dependency. The Journal of the Acoustical Society of America , 115 (4) : 1424-1430. DOI:10.1121/1.1646399 | |
| Day S M, Minster J B. 1984. Numerical simulation of attenuated wavefields using a Padéapproximant method. Geophysical Journal International , 78 (1) : 105-118. DOI:10.1111/j.1365-246X.1984.tb06474.x | |
| Emmerich H, Korn M. 1987. Incorporation of attenuation into time-domain computations of seismic wave fields. Geophysics , 52 (9) : 1252-1264. DOI:10.1190/1.1442386 | |
| Gazdag J. 1981. Modeling of the acoustic wave equation with transform methods. Geophysics , 46 (6) : 854-859. DOI:10.1190/1.1441223 | |
| Goffe W L, FerrierG D, Rogers J. 1994. Global optimization of statistical functions with simulated annealing. Journal of Econometrics , 60 (1-2) : 65-99. DOI:10.1016/0304-4076(94)90038-8 | |
| Kjartansson E. 1979. Constant Q-wave propagation and attenuation. Journal of Geophysical Research:Solid Earth , 84 (B9) : 4737-4748. DOI:10.1029/JB084iB09p04737 | |
| Li Q Z. High-resolution Seismic Data Processing (in Chinese). (in Chinese) Beijing: Petroleum Industry Press, 1993 : 38 . | |
| Liao J P, Liu H X, Wang H Z, et al. 2011. Two dimensional visco-acoustic wave modeling research in frequency-space domain. Progress in Geophysics (in Chinese) , 26 (6) : 2075-2081. DOI:10.3969/j.issn.1004-2903.2011.06.023 | |
| Liu H P, Anderson D L, Kanamori H. 1976. Velocity dispersion due to anelasticity:Implications for seismology and mantle composition. Geophysical Journal International , 47 (1) : 41-58. DOI:10.1111/j.1365-246X.1976.tb01261.x | |
| Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics , 49 (5) : 533-549. DOI:10.1190/1.1441689 | |
| Operto S, Virieux J, Amestoy P, et al. 2007. 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver:A feasibility study. Geophysics , 72 (5) : SM195-SM211. DOI:10.1190/1.2759835 | |
| Robertsson J O, Blanch J O, Symes W W. 1994. Viscoelastic finite-difference modeling. Geophysics , 59 (9) : 1444-1456. DOI:10.1190/1.1443701 | |
| Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid. Geophysics , 69 (2) : 583-591. DOI:10.1190/1.1707078 | |
| Wang R Z, Liu R, Zhao B X. 2014. Numerical simulation of wave equation with segmented smooth curve boundaries. Chinese Journal of Geophysics (in Chinese) , 57 (4) : 1199-1208. DOI:10.6038/cjg20140417 | |
| Xu Y, Zhang J F. 2008. An irregular-grid perfectly matched layer absorbing boundary for seismic wave modeling. Chinese Journal of Geophysics (in Chinese) , 51 (5) : 1520-1526. DOI:10.3321/j.issn:0001-5733.2008.05.026 | |
| Zhang J F, Liu T L. 1999. P-SV-wave propagation in heterogeneous media:Grid method. Geophysical Journal International , 136 (2) : 431-438. DOI:10.1111/j.1365-246X.1999.tb07129.x | |
| Zhang J F, Liu T L. 2002. Elastic wave modelling in 3-D heterogeneous media:3D grid method. Geophysical Journal International , 150 (3) : 780-799. DOI:10.1046/j.1365-246X.2002.01743.x | |
| Zhu T Y, Harris J M. 2013. A constant-Q time-domain wave equation using the fractional Laplacian.//SEG Technical Program Expanded Abstracts 2013.Society of Exploration Geophysicists, 3417-3422. | |
| Zhu T Y. 2014. Time-reverse modelling of acoustic wave propagation in attenuating media. Geophysical Journal International , 197 (1) : 483-494. DOI:10.1093/gji/ggt519 | |
| Zhu T Y, Harris J M. 2014. Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians. Geophysics , 79 (3) : T105-T116. DOI:10.1190/geo2013-0245.1 | |
| 李庆忠. 走向精确勘探的道路. 北京: 石油工业出版社, 1993 : 38 . | |
| 廖建平, 刘和秀, 王华忠, 等. 2011. 二维频率空间域黏声波正演模拟研究. 地球物理学进展 , 26 (6) : 2075–2081. DOI:10.3969/j.issn.1004-2903.2011.06.023 | |
| 王忠仁, 刘瑞, 赵博雄. 2014. 分段光滑曲线边界波动方程数值模拟研究. 地球物理学报 , 57 (4) : 1199–1208. DOI:10.6038/cjg20140417 | |
| 吴玉, 符力耘, 陈高祥. 2014. 基于分数阶拉普拉斯算子解耦的粘声介质地震正演模拟与逆时偏移.//2014年中国地球科学联合学术年会-专题19:地震波传播与成像论文集. 北京:中国地球物理学会. http://cpfd.cnki.com.cn/article/cpfdtotal-zgdw201410024117.htm | |
| 徐义, 张剑锋. 2008. 地震波数值模拟的非规则网格PML吸收边界. 地球物理学报 , 51 (5) : 1520–1526. DOI:10.3321/j.issn:0001-5733.2008.05.026 | |