引言

因介质的内摩擦力、含流体孔隙及非均匀性引起散射衰减等原因，地下介质通常表现为非完全弹性，因此地震波在地下介质中传播时发生衰减并对其振幅、相位等产生明显影响，从而降低了地震资料的分辨率.随着对地震勘探精度的要求越来越高，使用黏弹性模型可以更接近实际地描述地下介质的黏滞性，通过对黏弹性介质的地震波场进行数值模拟来研究和分析地震波在传播过程中的衰减特征，对指导地震资料处理及提高地震资料的分辨率具有重要意义.

在地震勘探中，通常使用品质因子Q来度量介质对地震波能量的吸收衰减.实验观测数据表明，在地震勘探频带范围内，地层的品质因子基本不随频率变化，此时衰减系数与频率成正比(Aki，Richards，2002).近三十年来，许多学者在黏弹性理论模型的建立及地震波衰减模拟等方面开展了大量研究工作.目前通常使用机械元件的组合来构建黏弹性模型，例如开尔文-佛格特(Kelvin-Voigt)模型、麦克斯韦(Maxwell)模型，但上述两种基本模型的品质因子与频率分别成反比和正比关系，其对高频成分的衰减分别存在衰减过量和衰减不足的问题，均不能真实地表达黏弹性介质的物理和力学特征(孙成禹，2007).将多个标准线性体并联得到的广义标准线性固体模型(Liu et al，1976；Day，Minster，1984；Carcione et al，1988a, b；Robertsson et al，1994；Zhu et al，2013)或者由多个麦克斯韦模型组合得到的广义麦克斯韦模型(Emmerich，Korn，1987)，可以在给定频带较好地拟合品质因子Q不随频率变化的特征，因此被广泛地用于构建黏弹性模型，被称为近似常Q模型(或称广义流变模型). Cao和Yin(2014)分析证明了不同流变模型之间的等价性. 孙成禹和印兴耀(2007)及孙成禹等(2010)提出了一种使用三参数来构造常Q黏弹性模型的方法，并对开尔文-佛格特模型、麦克斯韦模型的黏弹性介质有限差分正演进行了稳定性分析.由于这些流变模型引入了与应力、应变松弛时间有关的记忆变量，形式复杂，考虑到计算效率，在实际正演模拟和偏移成像等方面较少使用. Kjartansson(1979)建立了一种基于数学模式的完全常Q模型，可以定量地描述非弹性介质的地震波能量吸收和速度频散.依据该理论，Carcione等(2002)和Carcione(2009)建立了常Q黏声和黏弹性波动方程，并使用伪谱法进行了地震波场模拟.需要注意的是，该模型的松弛函数含有与时间相关的幂律形式，在将本构方程中的时间域褶积转化为乘积形式的过程中会引入时间分数阶导数.对于时间分数阶导数的差分近似，需要利用过去所有时刻的波场，尽管数值模拟中通常使用短时记忆方法截断差分算子(Carcione et al，2002)，但依然要耗费很大的计算量和内存资源.

基于波动方程理论的正演模拟方法主要有有限差分法、有限元法、伪谱法、边界元法等(孙成禹，2007).相较于其它方法，有限差分方法实现起来简单快速，占用内存小，可以较好地适应剧烈变化的复杂地下介质.然而由于有限差分法在对时间和空间导数离散近似的过程中存在误差，使得其在地震波场模拟中容易产生数值频散.交错网格有限差分方法可以有效地提高数值模拟的精度和稳定性(Virieux，1986)，但仍然会不可避免地产生数值频散.另外，黏弹性介质本身亦固有频散特征，因此在数值模拟过程中控制数值频散变得尤为重要.为压制数值频散，前人提出了很多优化交错网格有限差分系数的方法，主要分为基于最优化方法和基于窗函数截断逼近方法两类.例如，Liu(2014)提出了基于最小二乘法的交错网格差分系数优化方法，Wang等(2014)基于声波方程时空域频散关系，采用正则化方法优化了交错网格差分算子.相较于传统的通过泰勒级数展开得到差分系数的方法，这些方法均在一定程度上提高了差分精度，降低了数值频散效应，但由于涉及到控制参数，在实际应用中难以确定.与最优化方法相比，Chu和Stoffa (2012)提出了基于改进二项式窗函数截断的优化有限差分算子，该方法计算简单，并证明了有限差分法与伪谱法的统一性.

本文基于Kjartansson常Q模型，首先推导了含有分数阶导数的常Q黏声及黏弹性波动方程的速度-应力形式，在时间域采用基于改进二项式窗函数的优化交错网格有限差分法求解该方程，并结合卷积完全匹配层(convolutional perfectly matched layer，简写为CPML)边界条件(Komatitsch，Martin，2007)实现了对衰减介质的地震波场正演模拟.针对使用短时记忆方法计算时间分数阶导数效率低的问题，本文引入了自适应时间记忆步长方法，提高了波动方程的离散化精度和计算效率.通过对均匀模型的数值解与解析解的比较证明了该算法的精确性，并详细分析了子波在不同品质因子下的频散与衰减特征.最后对BP盐丘模型进行数值模拟并分析其衰减特征，验证了本文算法对复杂介质的适用性.

1.   时间分数阶波动方程

1.1   常Q模型

Kjartansson(1979)提出的常Q模型由两个参数所限定，即参考频率处的相速度和品质因子Q，且Q完全与频率无关.其松弛函数准确地描述了地下介质的线性衰减特征.松弛函数的形式为

式中：M₀为参考体积模量，M₀=ρc₀²cos²(πγ/2)，其中ρ为介质密度，c₀为参考速度；γ为与品质因子有关的无量纲量，γ=arctan(1/Q)/π (0 < γ < 0.5)；Γ为伽马函数；t为时间；t₀为参考时间.

根据Caputo分数阶导数的定义(Caputo，1967)，线性黏弹性介质中本构关系的褶积形式可转化为乘积形式，即

式中，σ(t)为应力，*为褶积运算，ə^2γε(t)/ət^2γ为应变ε(t)对时间的2γ阶分数阶导数.对式(2) 进行傅里叶变换，得到线性黏弹性介质在频率域的本构关系为

式中，M(ω)为复模量，M(ω)=M₀(iω/ω₀)^2γ，i为虚数单位，ω为角频率，ω₀为参考频率.

品质因子定义为

式中，Re和Im分别为复数的实部和虚部，由式(4) 可知，品质因子与频率无关.

1.2   分数阶黏声波动方程

声波在二维均匀介质中传播时的质点运动平衡方程为

式中，p为声压，u_x和u_z分别为质点在x，z轴方向的位移分量.

结合式(2)，(5) 和应变-位移关系(几何方程)可得

式中，c²=c₀²cos²(πγ/2)，▽²为拉普拉斯算子.显然，当γ→0(Q→∞)时，式(6) 趋于完全弹性介质下的二阶声波波动方程.

由式(5)，(6) 可推导出常Q衰减介质的声波速度-应力方程为

式中，v_x，v_z表示质点在x，z轴方向的振动速度，D^1－2γ为时间的1－2γ阶导数.

声波相速度频散关系和衰减系数分别为

式中，c_P为相速度，α为衰减系数.可以看出，两者均与频率有关：前者表明常Q衰减介质中波的传播速度是频率的函数，即存在频散现象；后者说明介质对地震波能量的吸收具有频率选择性.

1.3   分数阶黏弹性波动方程

常Q黏弹性各向同性介质的本构方程可表示为(Carcione，2009)

式中，σ_ij, ε_ij分别为应力、应变张量的分量， δ_ij为克罗内克符号，D^2γ_P和D^2γ_S分别表示P波和S波的时间分数阶导数，C_λ=M₀ω₀^－2γ_P，C_μ=μ₀ω₀^－2γ_S，M₀和μ₀分别为P波和S波的参考模量，表示为

式中，c_P0和c_S0分别为P波、S波在参考频率ω₀处的相速度. γ_P和γ_S与P波、S波的品质因子有关，表示为

弹性体运动平衡方程为

式中，f_i为体力分量，u_i分别为x，y，z方向的位移分量.

考虑二维情况，结合式(10)，(13) 及几何方程，得到常Q衰减介质的弹性波速度-应力方程为

明显地，当Q_{P, S}→∞(γ_{P, S}→0) 时，式(14) 趋于完全弹性一阶速度-应力波动方程.

黏弹性介质中，P波和S波的相速度分别为

P波、S波的衰减系数分别为

2.   时间分数阶导数自适应时间步长记忆离散方法

对于任意函数f(t)，其时间分数阶导数的中心差分可近似为(Carcione et al，2002)

式中，D^β为时间分数阶导数，β为阶数，t为计算时间，τ为时间步长，K=t/τ－1为时间点数，φ(β, k)为二项式系数，表示为

图 1显示了当β=0.98时，二项式系数的十进对数随k的变化趋势.可以看出，相对于当前时刻，过去不同时间点对式(17) 求和产生的影响不同，接近当前时刻的时间点影响较大，而随着时间间隔的增大，影响越来越小.利用这一性质，本文使用自适应时间步长记忆方法对整个时间历程通过加权采样进行离散，确保了离散精度，且该方法由于使用更少的时间节点参与计算，提高了计算效率.

图  1  二项式系数的十进对数随k的变化关系

Figure  1.  Decimal logarithm of binomial coefficients versus summation index k

下载: 全尺寸图片 幻灯片

图 2为时间分数阶导数两种离散方法的原理示意图.对于任意数值a，设初始区间为[0, a]，该区间内的所有时间节点均参与计算，且将过去时间历程分割为不同时间长度的区间，区间集合定义为

图  2  短时记忆方法(a)和自适应时间步长记忆方法(b)离散示意图(引自MacDonald等，2015)

k=0表示当前时刻；k为过去的时间节点数，即迭代次数；L为短时记忆长度；黑线为过去所有的时间节点；红线为参与计算的时间节点

Figure  2.  Short memory method (a) and adaptive time step memory method (b) for the discretization (after MacDonald et al, 2015)

k=0 represents the current moment; k represents the number of iterations; L represents short-term memory length; the black lines represent all the time nodes in the past; the red lines represent the time points included in the computation

下载: 全尺寸图片 幻灯片

根据加权采样的思想，每个区间内采样点的权重为2i-1，相应的时间点采样间隔d为

根据上述定义，使用自适应时间步长记忆方法将分数阶导数离散为

对于式(7) 中的时间分数阶导数D^1－2γp，在t=(n+1/2)Δt时刻的短时记忆中心差分离散形式为

同理，使用自适应时间步长记忆方法，在t=(n+1/2)Δt时刻将其离散化为

式中p_{i, j}ⁿ表示在空间节点(i, j)处时间节点n的声压.

3.   优化交错网格有限差分方法

使用交错网格差分方法求解常Q声波和弹性波速度-应力方程，是在不同网格节点上计算应力和速度分量，其在二维空间中的分布如图 3所示.以声波方程式(7) 中x方向的声压为例，空间域2N阶精度交错网格差分格式为

图  3  常Q声波和弹性波方程变量空间域交错网格分布

Figure  3.  Variables of constant-Q acoustic and elastic wave equations located in staggered grids in 2D space domain

下载: 全尺寸图片 幻灯片

式中，a_n为交错网格差分系数，p_{i, j}=p(x+ih, z+jh).将平面波p=p₀e^{i(ωt－k_xx)}带入上式有

可以看到，在波数域，有限差分方法的三角多项式逼近等效于伪谱法. Chu和Stoffa(2012)提出，交错网格有限差分系数a_n可以由窗函数截断逼近微分算子得到，即

式中w_n^N+M为改进的二项式窗函数，表示为

式中M为任意非负偶数.当M=0时，上式等同于常规交错网格有限差分算子，即传统的泰勒级数展开方法(Taylor series expansion method, 简写为TEM).

由式(25) 得到频散误差函数E为

图 4a，b比较了不同二项式窗函数(M=0, 4, 8, 16, 32) 的幅值响应及对应窗函数截断逼近得到的交错网格算子精度误差曲线，可以看出, 随着M的增大，频散误差有较大的波数覆盖范围，说明采用优化差分算子可达到更高的精度，但是，当M过大(M>16) 时，精度误差波动较大，稳定性差，对频散效果影响很大.

图  4  改进二项式窗函数(M=0, 4, 8, 16, 32) 的幅值响应(a)及截断逼近的交错网格有限差分算子精度误差曲线(b)

Figure  4.  Amplitude response of scaled binomial window functions (M=0, 4, 8, 16, 32) (a) and the accuracy error curves of finite-difference operators in truncated eight-order (2N=8) staggered grids (b)

下载: 全尺寸图片 幻灯片

4.   数值算例与波场分析

4.1   频散与衰减

McDonal等(1958)通过物理实验表明了皮埃尔页岩介质的常Q衰减特征，该模型纵、横波品质因子分别为Q_P≈32和Q_S≈10，纵、横波速度分别为v_P=2164 m/s和v_S=802 m/s(参考频率为100 Hz)，页岩密度约为2.2 g/cm³. 图 5a对使用短时记忆方法(L=30) 和自适应时间步长记忆方法(a=4) 近似得到的频散曲线分别与理论值进行了对比，图 5b将采用两种方法近似得到的衰减曲线分别与理论值进行对比，时间步长为0.2 ms.从图中可以看出，两种方法得到的频散曲线和衰减曲线与理论值均基本一致，且使用自适应时间步长记忆算法得到的结果与理论值更加吻合，尤其在低频处.

图  5  相速度(a)和衰减系数(b)随频率f的变化(a为初始区间长度，下同)

Figure  5.  Phase velocity (a) and attenuation coefficient (b) versus frequency f(a is the initial interval length, the same below)

下载: 全尺寸图片 幻灯片

4.2   均匀介质模型

以皮埃尔页岩为例，考虑二维均匀常Q衰减介质中波的传播情况.设模型网格数为301×301，网格间距为2 m，模型其它参数设置参考4.1节.选用主频为50 Hz的雷克子波，震源位于模型中心位置，时间步长为0.2 ms，记录长度为0.2 s.分别使用基于短时记忆(L=30) 和自适应时间步长记忆(a=4) 的时间分数阶导数离散方法及空间八阶(2N=8) 精度交错网格有限差分方法开展数值模拟. 图 6比较了接收点(231，151) 记录到的波形与解析解，后者由常Q格林函数和震源函数通过卷积运算得到(Carcione, 2007). 表 1比较了使用短时记忆方法和自适应时间步长记忆方法时的程序运行耗时及数值解误差，可以看到，采用自适应时间步长记忆方法在减小误差的同时还可以提高运算效率.数值解的均方根误差定义为

式中，nt表示时间采样点数，d_iⁿ和d_i^a分别为采样点i的计算值和解析值.

图  6  短时记忆方法和自适应时间步长记忆方法的数值解与解析解的比较

Figure  6.  Comparison between analytical and numerical solutions using the short memory method and the adaptive time step memory method

下载: 全尺寸图片 幻灯片

表  1  短时记忆方法和自适应时间步长记忆方法的程序运行耗时及相对误差

Table  1.  Execution time of program and relative error using the short memory method and the adaptive memory method

	计算时间/s	均方根误差
短时记忆方法(L=30)	52.0	0.0130
自适应时间步长记忆 方法(a=4)	26.8	0.0022

下载: 导出CSV 

| 显示表格

保持其它参数不变，研究不同品质因子在均匀介质中的地震波衰减特征. 图 7给出了4种情况(Q=∞，32，10，4) 下的频散和衰减特征对比图.可以明显看出：Q值越小，相速度的频散和衰减越严重，且相速度和衰减系数均与频率成正比；当震源子波的主频取50 Hz时，Q值越小，其相速度越小，且均小于声波速度. 图 8显示了不同Q值下140 ms时刻的波场快照，由于Q取不同值时子波相速度不同，造成波前不一致. 图 9为接收点(231，151) 所记录到的波形和相应的振幅谱.由图 9a可以看出，随着Q值的减小，子波振幅衰减越来越严重，相位越来越往后延迟，波形逐渐变宽；由图 9b可以看出，高频成分比低频成分衰减得快，有效频带变窄. 图 10为弹性波和黏弹性波模拟得到的波场快照，可见黏弹性情况下波场能量较弱. 图 11为检波点(133, 125) 记录的波形，震源选用主频为50 Hz的垂直集中力源.可以看出，横波较纵波衰减严重.

图  7  不同Q值下的频散(a)和衰减(b)曲线

Figure  7.  Dispersion (a) and attenuation (b) curves with different Q values

下载: 全尺寸图片 幻灯片

图  8  均匀介质取不同Q值时在140 ms时刻的波场快照

(a) Q无穷大；(b) Q=32；(c) Q=10；(d) Q=4

Figure  8.  Snapshots of wave field recorded at 140 ms in homogeneous media for different Q values

(a) Infinite Q; (b) Q=32; (c) Q=10; (d) Q=4

下载: 全尺寸图片 幻灯片

图  9  距震源160 m的接收点(231, 151) 记录到的子波波形(a)和振幅谱(b)

Figure  9.  Waveforms (a) and amplitude spectrum (b) at the receiver (231, 151) with a distance of 160 m from source

下载: 全尺寸图片 幻灯片

图  10  弹性(a)和黏弹性(b)均匀介质在140 ms时刻的垂直分量波场快照

Figure  10.  Snapshots of vertical component recorded at 140 ms in elastic (a) and viscoelastic (b) homogeneous attenuating media

下载: 全尺寸图片 幻灯片

图  11  距震源约60 m的接收点(133, 125) 记录到的子波波形对比

Figure  11.  Comparison of waveforms recorded at the reference point (133, 125) with a distance of about 60 m from the source

下载: 全尺寸图片 幻灯片

4.3   BP盐丘模型

为了验证本文算法对复杂模型的适应性，以BP盐丘模型为例进行正演模拟.模型速度和密度场如图 12所示，品质因子使用李氏经验公式计算得到(Q=14×v^2.2).模型大小约为7.5 km×11.9 km，离散网格点数为1201×1911，网格间距为Δx=Δz=6.25 m，震源采用主频为30 Hz的雷克子波，时间步长为0.2 ms，炮点位于网格点(600，50) 处，整个地表接收，记录时长为5 s，模型四周采用CPML边界，吸收带宽为20个网格点，分别采用常规交错网格和优化交错网格有限差分方法(时间记忆步长a=4、空间八阶精度)进行声波和黏声波数值模拟，地震记录如图 13所示.对比图 13a，b可以看出，使用二项式窗优化交错网格差分算子可以有效地压制数值频散，改善地震记录的质量.由图 13c可以看出，黏声介质中地震波场的能量明显降低，深层反射振幅减弱，分辨率较低，主要表现为记录到的反射波在时间轴上被拉长，同相轴变宽.

图  12  BP盐丘的速度(a)和密度(b)模型

Figure  12.  Velocity (a) and density (b) of BP salt model

下载: 全尺寸图片 幻灯片

图  13  BP盐丘模型的地震记录

(a)常规交错网格差分算子；(b)基于改进二项式窗优化交错网格差分算子；(c)常Q衰减介质下优化交错网格差分算子

Figure  13.  Seismograms of BP salt models

(a) Conventional staggered grid difference operator; (b) Optimal staggered grid difference operator by scaled binomial windows; (c) Optimal staggered grid difference operator in constant-Q attenuation media

下载: 全尺寸图片 幻灯片

5.   讨论与结论

针对目前表征常Q特征的广义流变模型形式复杂、数值模拟效率低的问题，本文推导了基于Kjartansson常Q模型的本构关系及分数阶波动方程，该模型的波动方程形式简单，各参数均有明确的物理意义.使用优化交错网格有限差分方法求解常Q声波和弹性波速度-应力波动方程，数值模拟结果表明，该方法可以准确地模拟常Q衰减介质中地震波的传播过程，模拟精度较高；将CPML吸收边界条件应用于数值模拟中，可以有效地压制人工边界的反射.针对基于短时记忆方法的离散时间分数阶导数计算耗时长的问题，使用自适应时间步长记忆方法可以有效地提高计算效率和离散化精度.另外，使用改进的二项式窗函数得到的优化交错网格差分算子可以很好地压制空间数值频散.常Q模型可以较好地反映介质的黏滞特征和对地震波的吸收效应.通过对Kjartansson常Q模型地震波衰减特征的研究，为与其它衰减模型及实际介质中地震波衰减特征的对比提供了参考，同时也为进一步研究介质品质因子的反演、衰减补偿偏移等提供了有益的借鉴.另外，可以将该模型推广到各向异性介质，发展常Q黏弹性各向异性理论.