非均匀黏弹性介质中声波模拟的非规则网格方法

豆辉，张剑锋

1 中国科学院地质与地球物理研究所，中国科学院油气资源研究重点实验室，北京 1000292 中国科学院大学，北京 100049



非均匀黏弹性介质中声波模拟的非规则网格方法

豆辉1,2，张剑锋1

1 中国科学院地质与地球物理研究所，中国科学院油气资源研究重点实验室，北京 1000292 中国科学院大学，北京 100049

本文研究了满足衡Q模型的黏弹性介质中声波模拟的时间域数值解法.通过采用有理式基函数描述频率相关的体积模量，使模型满足地震波黏性吸收的衡Q要求.结合弹性波模拟的非规则、非结构网格方法——格子法，本文提出了一种非均匀黏弹性介质中声波模拟的非规则网格方法.非规则、非结构网格的使用，可以精细地刻画地下介质中复杂速度、Q值的变化，及界面的形状.通过和解析解的比较及复杂模型算例分析，验证了该方法的精度及有效性.

时间域；地震波模拟；非规则网格；黏弹性；衡Q模型；有理函数

1 引言

地震波数值模拟对于定量地分析研究复杂地下介质中的地震波成像非常重要.传统的地震波模拟基于弹性模型假设，可以满足一般的地震勘探精度要求，得到很好的数值模拟结果(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)

格子法的理论核心是在非规则网格下，给出式(1)的基于积分平衡的微分方程弱形式.如图1所示，在节点k的邻域，定义虚线轮廓线1-2-3-…-12-1所围区域为Ω，简记该轮廓线为∂Ω，对式(1)两边作面积分，并对等式右边应用Green定理得到

(2)

式中m为围绕节点k的三角形单元数目，Skl对应节点k周围第l个三角形单元中的虚线段，α与β分别是边界Ω外法线沿x与z轴的方向余弦.利用动力学计算的集中质量原理及三角形线性插值方法可得到式(2)的离散形式：

(3)

式中Mk是节点k各邻域三角形单元中的∬Ω1/ρv2dΩ值总和的三分之一，(ρ)l表示第l个单元的密度，Dxp,Dzp是三角形差分算子，可表述为

(4)

图1 非规则三角网格局部示意图Fig.1 Local irregular triangle grid

式中，Δ为三角形ijk的面积，下标r代表三角形中的各顶点，(bk)l，(ck)l分别表示围绕节点k的第l个单元对应的几何参数，以图1中的三角形单元ijk为例，可简单地表述为bk=(zj-zi)/2,ck=(xj-xi)/2.

最后，利用中心差分格式离散方程式(3)左端中对时间的二阶偏导数，我们就可以实现声波波场的迭代更新.

3 体积模量的有理函数展开和衡Q模型

3.1 体积模量的有理函数展开

黏弹性介质在频率域中有如下应力应变关系：

(5)

其中，M(ω)是与频率有关的复体积模量，σ(ω)是应力，ε(ω)为应变，经由Fourier变换，得到如下褶积关系：

(6)

这意味着在时间域计算时，需要存储历史时刻的应变值，会导致极大的计算负担，是数值计算很不乐见的.所以，我们需要考虑从其他方面对其做简化近似处理.

已知实验测得实际介质的应力松弛响应关系为：在常应变作用下，应力呈指数衰减变化，即实际黏性介质的松弛函数可用一个负指数衰减的函数表示.不同黏性介质的松弛响应对应不同的幅值及衰减系数表示.我们总可以找到这样的一组基函数——由一系列幅值和衰减系数组合表示任意的衰减函数，并用其表示介质的松弛响应.我们定义未松弛模量MU及松弛模量MR，分别表示介质衰减前、后的体积模量，并基于此构建如下负指数衰减的松弛函数：

(7)

至此，可将体积模量表示为

(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)

图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

一般低阶(n≤8)时有2n+1个模型参数，还比较少，我们可以采用蒙特卡洛、模拟退火、遗传算法等全局搜索方法.本文采用的是优化的模拟退火方法(Goffeetal.，1994)，可实现全局的快速搜索，求得最优解.然后将最优解对应的参数结果列表记录下来，便于以后模拟计算时查表使用.最终得到的结果和目标值的拟合效果，可以参见图2所示，图中实线表示目标值Q0=20.0，虚线为Q-1(ω)拟合结果.图2中，分别对应阶数n=2，n=3时Q(ω)拟合得到的曲线结果和目标值的比较.可以看出，在频带3～100 Hz内，n=2时即可得到不错的拟合效果；n=3时，Q(ω)曲线可以更加贴合目标值，拟合效果更优.本文数值模拟的简单算例，采用的就是阶数n=3时的参数结果.

假设δM/MU≪1，则式(10)可近似为

(12)

4 非均匀黏弹性介质中声波模拟的非规则网格方法

将式(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定理，得到

(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分裂方程，这里不再详述.该方法可以很好地吸收或透射不同方向的入射波，避免边界反射问题.

5 数值算例

为验证本文提出的非规则网格黏声波数值模拟解法的正确性和精度，首先将其应用于均匀速度模型(对此模型可容易得到黏声波传播的解析解)，并将其计算结果与解析解对比分析.然后，给出一个强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

5.1.1 不同炮检距时

本文首先研究Q=20时，数值解对炮检距的敏感度——能否在远炮检距处保证相应的计算精度.我们分别在距震源600 m，1600 m，3200 m处，放置一个检波器记录波场值，并与相同位置的解析解对比，最终结果如图3所示(图中实线表示解析解，虚线表示数值解).从图3中可以看出，该算法的数值解与解析解都能较好的吻合.为进一步定量地评估本文算法的精度，本文利用数值解与解析解的均方根误差来衡量，并定义均方根误差E：

(24)

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值很小时，本文解法仍可得到稳定的高精度的计算结果.

图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处均具有稳定的高精度解.

图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

5.2 气云模型

具有强衰减(即低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变化模型中可得到稳定的计算结果.另一方面，鉴于非规则、非结构网格的优势，计算时气云内使用黏声波方程，气云外的区域按弹性波方程计算，可以有效节省黏声波计算时的存储空间，提高计算效率，是其他规则网格方法所不具备的.

5.3 Marmousi模型

图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

6 结论

本文发展了基于非规则化网格的时间域黏声波数值模拟算法.文中采用有理式基函数描述频率相关的体积模量，通过无穷范数构建目标函数，并结合模拟退火方法，实现了满足衡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.

Carcione J M.2010.A generalization of the Fourier pseudo spectral method.Geophysics,75(6):A53-A56.

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.

Day S M,Minster J B.1984.Numerical simulation of attenuated wavefields using a Padéapproximant method.Geophysical Journal International,78(1):105-118.

Emmerich H,Korn M.1987.Incorporation of attenuation into time-domain computations of seismic wave fields.Geophysics,52(9):1252-1264.

Gazdag J.1981.Modeling of the acoustic wave equation with transform methods.Geophysics,46(6):854-859.

Goffe W L,FerrierG D,Rogers J.1994.Global optimization of statistical functions with simulated annealing.Journal of Econometrics,60(1-2):65-99.

Kjartansson E.1979.Constant Q-wave propagation and attenuation.Journal of Geophysical Research:Solid Earth,84(B9):4737-4748.

Li Q Z.1993.High-resolution Seismic Data Processing (in Chinese).Beijing:Petroleum Industry Press,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.

Marfurt K J.1984.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations.Geophysics,49(5):533-549.

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.

Robertsson J O,Blanch J O,Symes W W.1994.Viscoelastic finite-difference modeling.Geophysics,59(9):1444-1456.

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.

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.

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.

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.

Zhu T Y,Harris J M.2014.Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians.Geophysics,79(3):T105-T116.

附中文参考文献

李庆忠.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:地震波传播与成像论文集.北京:中国地球物理学会.

徐义,张剑锋.2008.地震波数值模拟的非规则网格PML吸收边界.地球物理学报,51(5):1520-1526,doi:10.3321/j.issn:0001-5733.2008.05.026.

(本文编辑 胡素芳)

An irregular grid method for acoustic modeling in inhomogeneous viscoelastic medium

DOU Hui1,2，ZHANG Jian-Feng1

1 Key Laboratory of Petroleum Resources Research，Institute of Geology and Geophysics， Chinese Academy of Sciences，Beijing 100029，China2 University of Chinese Academy of Sciences，Beijing 100049，China

Absorption and dispersion in earth media is an issue we should consider for wavefield simulation.In this paper,we present a time-domain numerical modeling algorithm for acoustic wave in inhomogeneous viscoelastic medium.On one hand,to describe the attenuation effect we use the constant-Q model,which can be frequency-independent.On the other hand for the wavefield modeling engine,we utilize the grid method,which uses irregular,unstructured grid to simulate acoustic wave propagation in non-uniform viscoelatsic media.

Time-domain；Seismic modeling；Irregular grid；Viscoelasticity；Constant-Q model；Rational function

豆辉，张剑锋.2016.非均匀黏弹性介质中声波模拟的非规则网格方法.地球物理学报，59(11):4212-4222,

10.6038/cjg20161123.

Dou H，Zhang J F.2016.An irregular grid method for acoustic modeling in inhomogeneous viscoelastic medium.Chinese J.Geophys.(in Chinese),59(11):4212-4222,doi:10.6038/cjg20161123.

国家自然科学基金面上项目(41574135)和国家自然科学基金重点项目(41330316)联合资助.

豆辉，女，1989年生，中国科学院地质与地球物理研究所博士研究生，主要从事时间域数值模拟研究.E-mail:douhui@mail.iggcas.ac.cn

10.6038/cjg20161123

P631

2016-03-17，2016-09-19收修定稿

Specifically to resolve the convolution problem between strain and stress in the time-domain,we approximate the bulk modulus in the frequency domain with a set of rational functions,which meets the constant-Q requirement.To confirm the model parameters for constant-Q,we apply the L-∞ Norm to reconstruct the objective function,and choose the simulated annealing (SA) method to finish the searching.As a result,we find the optimum solution with low order,like n=2 or 3,which simultaneously meets the constant-Q requirement very well,and reduces the memory demand obviously.In addition,the irregular,unstructured grid,which is used to discretize the velocity model and corresponding Q model,can give detailed descriptions and differentiations when complex interfaces are present,without sacrificing the accuracy and efficiency.

As for the numerical examples,we first validate the scheme proposed using simple model and compare the results with analytic solutions.Result with good consistency is obtained regarding the attenuation effect.Moreover the Marmousi model validation proves that the proposed irregular,unstructured gird based simulation method together with the constant-Q model is effective and efficient to describe the acoustic wave propagation in inhomogeneous viscoelastic medium.