1 引言

实际的地球介质不是完全弹性的，是黏声介质.传播过程中地震波发生频散，导致能量逐渐衰减，相位畸变.基于声波方程的常规地震模拟与成像方法忽略了黏声介质的这种吸收衰减作用.为了更准确描述波在实际介质中的传播规律，需要在波场模拟中考虑介质的吸收衰减作用 (傅旦丹等，1993；张智等，2005；Dvorkin and Mavko, 2006；单启铜和乐友喜，2007).长期以来，对黏声介质的吸收衰减机制进行了广泛的研究，相关衰减模型可以大致分为两类：(1) 在频域引入复速度 (Aki and Richards, 2002; Liao and McMechan, 1996; Štekl and Pratt, 1998); (2) 在时域波动方程中引入品质因子Q.其中，时域的方法通常用一系列黏性参数 (即标准线性体) 的组合来描述介质的黏性，其相关模型被称为近似的常Q模型 (Liu et al., 1976; Emmerich and Korn, 1987; Carcione, 2007; Zhu et al., 2013).Kjartansson (1979)提出的常Q模型，衰减因子是频率的线性函数，其Q在常规地震勘探的频带内与频率基本无关.基于常Q模型推导的波动方程含有分数阶时间导数 (Caputo and Mainardi, 1971; Mainardi and Tomirotti, 2010)，Carcione采用Grünwald-Letnikov (GL) 逼近来近似求解分数阶时间导数，实现了黏声波方程和黏弹性波方程的波场模拟 (Carcione et al., 2002; Carcione, 2008).因为一个变量的分数阶时间导数依赖于该变量此前时刻的所有值，因此，直接数值求解黏声波方程需要保存每一时刻的应力-应变值，通常需要很大的内存，计算效率较低.可以对算子进行截断 (Podlubny, 1999; Carcione et al., 2002)，适当降低精度来减少内存需求.即使如此，内存需求问题依然不能忽略，特别是黏声介质三维地震模拟.

为了解决内存需求问题，Lu和Hanyga (2004)采用一系列中间变量来求解分数阶导数，但是这种方法以增加计算量为代价.Chen和Holm (2004)采用分数阶拉普拉斯算子计算黏声介质中的波场，这种分数阶拉普拉斯算子可以通过一种改进的伪谱法在波数域中方便求取，该算子把分数阶时间导数转化为分数阶空间导数，提高了计算效率.Carcione (2010)直接采用改进的伪谱法进行黏声介质中的地震波场模拟，由于衰减项和频散项是耦合在一起的，数值求解需要大量内存.Treeby和Cox (2010)进一步将频散项和衰减项关联方程进行显式分解，极大地提高了数值求解的计算效率，有利于反问题的求解.特别是在地球物理领域，黏声介质地震偏移，实现了在波场延拓过程中同时进行频散校正和衰减补偿，提高成像质量 (Zhu, 2014; Zhu et al., 2014).

本文从Kjartansson常Q应力-应变关系出发，结合声波介质的运动方程和动力学方程，推导了含分数阶时间导数的时间域常Q黏声波方程.该方程的求解由于内存需求大，不利于逆时偏移的实施.本文基于Treeby-Cox分解方法 (Treeby and Cox, 2010)，通过一系列近似，推导了频散项和衰减项解耦的分数阶拉普拉斯算子黏声波方程.通过交错网格有限差分逼近时间导数，改进的伪谱法求解分数阶拉普拉斯算子，本文采用PML边界条件，利用含分数阶拉普拉斯算子黏声波方程进行黏声介质一维和二维地震波场正演模拟.由于具有较高的计算效率，该方程主要用于黏声介质的地震数据逆时偏移，改善深层构造的成像精度.黏声介质逆时偏移需要解决的问题是波场衰减补偿后延拓过程的稳定性.在延拓过程中，波场衰减补偿的同时，高频噪声能量得到过度放大，导致数值延拓发散.因此，需要在每步波场延拓后进行带通滤波，消除高频噪声.模型试算的结果表明本文方法可以进行较高精度的波场正演模拟，逆时偏移可以有效地进行频散校正和衰减补偿，显著提高成像质量.

2 常Q模型和分数阶拉普拉斯算子 2.1 常Q模型 2.1.1 应力-应变关系

尽管有实验数据表明在超声频段Q与频率密切相关 (Kan et al., 1983; Sams et al., 1997; Batzle et al., 2006)，但对于低频段的地震勘探而言，常Q模型是对实际衰减介质的一个很好近似 (McDonal et al., 1958; Aki and Richards, 2002).Kjartansson (1979)给出了如下的松弛函数来模拟常Q介质的线性衰减特性：

(1)

其中M₀=ρ₀c₀²cos (πγ/2) 是体积模量, ρ₀是介质密度, c₀是参考频率ω₀处的相速度, t₀是参考时间, H(t) 是赫维赛德阶跃函数, Γ是伽马函数.参数γ=arctan (1/Q)/π是一个无量纲量，因为Q>0, 所以0 < γ < 0.5

黏弹性介质的应力-应变关系可以表示为松弛函数的卷积：

(2)

其中*表示卷积，公式 (2) 包含关于时间的分数阶导数.当Q→∞, γ→0, 公式 (2) 退化为完全弹性介质的本构方程σ=M₀ε.当0 < Q < ∞, 分数阶时间导数项 (∂^2γε/∂t^2γ) 的求解需要知道以前时刻的应变，即当下时刻的应力与之前的应变历史有关.因此，数值求解公式 (2)，内存需求量大，计算效率低.

2.1.2 相速度与衰减因子

对公式 (2) 进行傅里叶变换，可以得到：

(3)

其中复模量.对于声波介质，相速度c_p是角频率ω和实波数k的函数.对于黏声介质，假定k_R是复波数k的实部，则有：

(4)

其中.根据复模量M与实模量M₀的关系，可得：

(5)

其中.假定k_I是复波数k的虚部，则衰减因子α可表示为：

(6)

可见，衰减因子α与角频率满足线性关系.据此，定义黏声介质的品质因子为：

(7)

2.2 含分数阶拉普拉斯算子的常Q黏声波方程

根据声波介质的动力学方程：

(8)

和运动方程：

(9)

其中σ是压力场，ε是应变场，v=(v_x, v_y, v_z) 是质点的振动速度，代入公式 (2) 可得黏声波方程：

(10)

其中∇ ²是拉普拉斯算子，当Q→∞时，γ→0，方程 (10) 退化为经典的声波方程.数值求解方程 (10) 需要保存各个时刻的波场σ值，极大降低计算效率.为了解决这个问题，下面推导方程 (10) 的另一种形式.

把平面波exp (iωt-ik · r) 代入 (10) 式，得到相应的频散方程：

(11)

其中复波数k=k_R+ik_I.把i^2γ=cos (πγ)+isin (πγ) 代入 (11) 式，

(12)

在 (12) 式右边乘以c₀^2γ/c₀^2γ，当γ很小，可采用近似，整理得到：

(13)

，所以，代入 (13) 式得到：

(14)

其中的，(14) 式可近似为：

(15)

(15) 式为新得到的频散关系.

在均匀介质中，c₀和Q都为常数，对式 (15) 进行空间二维傅里叶逆变换和时间一维傅里叶逆变换，并考虑到如下的对应关系：k^2γ+2→(-Δ ²)^γ+1, k^2γ+1→(-∇²)^γ+1/2, (iω)²→ ∂² ∂t², 可得

(16)

其中η=c₀^2γω₀^-2γcos (πγ)，τ=c₀^2γ-1ω₀^-2γ(iω) sin (πγ).

方程 (16) 是新得到的黏声波方程，含有分数阶拉普拉斯算子，和方程 (10) 相比，把时间分数阶导数转移到对空间求分数阶导数.空间分数阶导数可以在波数域方便求取.数值求解方程 (16) 只需保存前一时刻的波场 (一阶时间导数项)，不需要把每一步的波场都保存下来，节省了内存，也提高了计算效率.方程 (16) 还有一个优越的特性是把介质的频散效应和衰减效应显式的区分开来.方程右边第一项主要和频散有关，而第二项主要与衰减有关，这两项的不同作用将在后续的数值试验中得到解释.

2.3 分数阶拉普拉斯算子计算

直接采用差分计算分数阶拉普拉斯算子比较复杂，如果把波场变换到波数域，在波数域求解分数阶拉普拉斯算子将变得非常方便.根据空间微分算子的波数域表示∂_x→ik_x，Carcione (2010)将空间微分算子从整数阶扩展到分数阶，通过如下的傅里叶变换过程：u(x_j)→FFT→ u_r→(ik)^βu_r→FFT^-1→∂_x^βu(x_j)，波数域的分数阶拉普拉斯算子可表示为：

(17)

其中F表示傅里叶变换.

3 正演模拟和逆时偏移 3.1 正演模拟

为了应用PML吸收边界，需要将方程 (16) 转化为如下的一阶方程组：

(18)

式中的一阶空间导数采用伪谱法在波数域计算，一阶时间导数采用有限差分计算，分数阶拉普拉斯算子按2.3节的方法进行计算.数值离散后的方程为：

(19)

3.2 逆时偏移

逆时偏移和正演模拟依据的方程是相同的，不同的是正演模拟是已知时刻n的波场σⁿ，求时刻n+1的波场σⁿ⁺¹，逆时偏移刚好反过来，所以离散方程为：

(20)

逆时偏移的主要问题是保持波场延拓的稳定性.完全弹性介质中，正演模拟和逆时偏移在参数选取恰当的情况下都是稳定的.在黏性介质中，波场正向传播能量是逐渐衰减的，所以是稳定的.但在逆时延拓的过程中，延拓波场逐层得到补偿，能量逐渐增大，需要数值稳定性特殊处理.波场延拓过程中，能量补偿项是，补偿能量与频率成正比，导致高频噪声能量得到极大增强，破坏波场延拓的稳定性，所以需要进行滤波处理，把高频噪声去除掉.

4 地震正演模拟数值算例 4.1 1D正演数值模拟

为了说明本文正演方法的可行性和准确性，首先进行一维情况下的数值解和解析解的对比.一维声波方程初值问题的解析解为, 其中g(x) 为初值条件.假设g(x)=(1-2π²k₀²x²) e^-π2k₀2x2, 取k₀=0.0462 m^-1.数值模拟参数为：c₀=2164 m·s^-1，ρ₀=2200 kg·m^-3，Q=32.5, nx=1204, dx=1 m, dt=1.15×10^-4s, t_max=0.4 s，接收距离为400 m.图 1是1D情况下数值解和解析解的对比，二者吻合得很好，说明本文方法的可靠性.为进一步验证数值模拟的精度，我们分别计算了上述模拟波场的衰减因子和相速度.衰减因子和相速度的计算公式为：

(21)

(22)

其中A₁，A₂是两个不同位置的接收器接收波场的振幅谱，ϕ₁，ϕ₂是对应的相位谱，d是两个接收器之间的距离.波场模拟参数同上例，d=50 m.计算的衰减因子和相速度见图 2, 图中实线是相应的理论值，根据公式 (5) 和 (6) 计算得到.从图 2可以看出，理论值和数值模拟的结果吻合很好，只是相速度在低频和高频端出现一些偏差.高频端偏差是因为高频衰减比较严重，被数值噪声污染，影响相位谱计算的精度，后续的例子分析进一步印证此结论.低频端偏差是因为本文算法对低频的频散衰减存在低估.

本文黏声介质正演模拟算法的一个显著特点是将频散效应和衰减效应显式分开来处理.为了研究两种效应的不同作用，我们对频散效应和衰减效应分别进行了数值模拟.图 3是仅考虑频散项得到的衰减因子和相速度, 因为没有考虑衰减项，所以衰减因子为0.值得注意的是，由于没有考虑衰减，相速度的数值解与理论值吻合较好 (特别是高频端).图 4是仅考虑衰减项的情形，从相速度上可以看出几乎没有频散.

4.2 2D正演数值模拟

为了验证本文正演方法在二维情况下的可行性和精度，图 5是四种介质二维数值模拟110 ms的波场快照.数值模拟的参数为：黏声介质Q=32.5，网格大小为512×512，网格尺寸为1 m×1 m, 震源位于模型中心，主频为100 Hz, 其他参数同上例.从图 5可见，四种介质的波场特征与理论解释基本吻合.图 5b为黏声介质仅考虑衰减，其波场能量比图 5a的声波介质明显衰减，但二者波前面走时一致，表明波速没有发生频散.图 5c为黏声介质仅考虑频散，其波场能量与图 5a的声波介质基本一致，没有衰减，但其波前面明显延后，表明存在波速频散.图 5d为黏声介质衰减与频散同时考虑，其波形既有能量衰减又有波速频散.

本文算法和直接求解分数阶黏声波方程 (10) 相比，把时间分数阶导数转移到对空间求分数阶导数.空间分数阶导数可以在波数域方便求取.数值求解方程 (16) 只需保存前一时刻的波场 (一阶时间导数项)，不需要把每一步的波场都保存下来，节省了内存，也提高了计算效率.为了定量说明本文算法计算效率的提升，沿用图 5的模型，分别采用本文方法和直接求解方程 (10) 进行波场模拟.从图 6可以看出，两种方法得到的结果相似，表明两种方法都可以准确用于黏声介质的地震波模拟.直接数值求解方程 (10) 用了2334 s, 本文方法只用了276 s, 计算效率差不多提高了8.5倍.图 7比较了黏声介质Q=100, 30, 10, 4时二维数值模拟110 ms的波场快照.其他数值模拟参数与图 5相同.从图中可以看出，Q值越小，波前面越延后，波场能量越弱，表明衰减和频散越严重.

5 地震逆时偏移数值算例

时间域常Q黏声波方程 (16)，由于其频散和衰减解耦的特点，以及较高的计算效率，该方程最佳的应用在于黏声介质的地震数据逆时偏移，实现对波场进行频散校正和衰减补偿，有效改善深层构造的成像精度.由于波在黏声介质中传播高频衰减较快，因此，在波场逆时延拓过程中高频能量补偿较大，弱振幅的高频噪声同时也得到极大增强，破坏了波场逆时延拓的稳定性.为了解决这个问题，需要在波场延拓的每一步，进行滤波去除噪声，避免波场中高频噪声不断累积放大.为了了解黏声介质波场逆时延拓过程中频散校正和衰减补偿的直观效果，我们进行如下一维数值试验，首先将初始波场在黏声介质中正向延拓到一定深度，然后再采用不同的校正进行逆时延拓到零时刻来重建初始波场.数值模拟参数为nx=2000, d x=1 m, c₀=2000 m·s^-1, ρ₀=2200 kg·m^-3, Q=10, t_max=0.6 s, dt=2.5×10^-4s.图 8a比较了既无频散校正又无衰减补偿的逆时延拓结果与初始波场，可见重建波场的相位和振幅与初始波场差距较远.图 8b考虑频散校正进行逆时延拓但无衰减补偿，可见重建波场的相位得到校正，但振幅与初始波场差距较远.图 8c同时考虑频散校正和衰减补偿的逆时延拓，重建波场的相位和振幅得到了精确恢复.

图 9是一个凹陷结构的速度和Q模型，网格点数为201×201，采样间隔为dx=dz=5 m.图 10是黏声波方程正演模拟得到的自激自收地震记录，可见深层的波场能量明显衰减，由于频散而频率降低，波形变宽.图 11是分别采用声波方程 (a) 和黏声波方程 (b) 对图 10地震记录进行逆时偏移成像结果.可见，声波方程偏移结果由于未考虑能量补偿和频散校正，深层波场振幅较弱，波形频散，模型界面的子波起跳点成像误差较大.这些问题通过黏声波方程偏移得到了很好的解决.

6 结论

从Kjartansson常Q应力-应变关系出发，结合声波介质的运动方程和动力学方程，本文推导了含分数阶时间导数的时间域常Q黏声波方程.求解该方程需要保存所有时刻的波场值，内存需求大，计算效率低，不利于地震偏移的实施.基于Treeby和Cox的分解方法，通过一系列近似，本文推导了含分数阶拉普拉斯算子的时间域常Q黏声波方程.该方程的优点表现在两个方面：(1) 把分数阶时间导数转换到空间导数上，利用改进的伪谱法，分数阶空间导数可以在波数域方便快速求取；(2) 将频散效应和衰减效应显式解耦，实现分别模拟介质的频散效应和衰减效应.求解该方程仅需保存前一时刻的波场值，内存需求减小，提高了计算效率.

本文将含分数阶拉普拉斯算子的时间域常Q黏声波方程转为了一阶方程组，通过交错网格有限差分逼近时间导数，改进的伪谱法求解分数阶拉普拉斯算子，采用PML边界条件，进行黏声介质一维和二维地震波场正演模拟.模型试算的结果表明本文方法可以进行较高精度的波场正演模拟，模拟波场和理论波场吻合得较好，模拟得到的频散关系、衰减关系和理论值也吻合得较好.本文进而开展黏声介质的地震数据逆时偏移，通过偏移波场的频散校正和衰减补偿，改善深层构造的成像精度.由于波在黏声介质中传播高频衰减较快，因此，在波场逆时延拓过程中高频能量补偿较大，弱振幅的高频噪声同时也得到极大增强，破坏了波场逆时延拓的稳定性.因此，需要在每步波场延拓后进行带通滤波，消除高频噪声.模型试算结果表明：地震数据的黏声介质逆时偏移，可以实现偏移波场的频散校正和衰减补偿，显著提高成像质量，得到的偏移剖面明显优于常规声波偏移剖面.