基于常Q模型的分数阶粘弹介质数值模拟方法

杨佩城， 韩立国， 蔡中正

(吉林大学 地球探测科学与技术学院，长春 130026)

0 引言

地震波在实际地下传播过程中，经历了大地滤波作用后，将导致能量的衰减和相位的错动。能量的衰减将严重影响勘探的深度；相位的错动则将影响构造位置的确定。早期的波场模拟和成像并未考虑这些衰减因素，对油气勘探造成了一定程度的误导。随着近年来勘探精度要求的提高，基于粘滞波动方程的正演和偏移算法也已经得到越来越多的重视。Kjartansson[1]提出的常Q模型，认为与地震勘探频率有关的频段内，Q值不随频率变化。近些年常Q模型得到了广泛地研究和应用。Zhang等[2]提出了一种基于Kjartansson常Q模型的分数阶粘弹性波动方程，该方程中相位校正项与衰减补偿项是解耦的；Deng 等[3]研究了基于标准线性体模型的逆时偏移衰减补偿，此方法能够在地震波传播过程中对全路径进行振幅补偿和相位校正，从而克服了反Q滤波方法的局限性；张平等[4]运用谱元法，对地震波在各向异性粘弹性介质中的传播进行了波场模拟以及分析；孟凡顺等[5]推导了二维非均匀黏滞性介质中地震波传播的二阶显示差分格式，并实现了对任意复杂地质体黏滞性波场进行的正演模拟。此外，为了提高计算效率，王德利等[6]使用MPI进行三维标准线性体波动方程有限差分并行数值模拟。

近年来，基于常Q模型，有一种包含非整数阶的波动方程被提出[7]。有研究表明，分数阶波动方程能更好地刻画地震波在粘弹性介质中的衰减和频散过程[8]。Carcione[9-10]采用Grunwald-Letnikov逼近方法近似求解时间分数阶波动方程，实现了粘滞声波和弹性波的波场模拟。然而，计算时间分数阶偏导数需要使用当前时刻以前所有时刻所对应的波场值，这会导致巨大的内存需求。为了解决这一问题，Carcione[11]又将时间偏导数转换为分数阶拉普拉斯算子，该算子可以通过快速傅里叶变换进行计算。这种算法避免了对当前时刻以前所有波场值进行存储这一问题，同时数值计算也大为简化。在此基础上，Zhu等[12]推导了基于解耦拉普拉斯算子的粘滞波动方程。解耦的拉普拉斯算子含有两项，其中一项表征振幅衰减，另外一项表征相位错动。该解耦的拉普拉斯算子应用于逆时偏移(RTM)时，只需将表征振幅衰减的一项取负号，而保持表征相位错动一项符号不变即可同时实现振幅的补偿和相位的校正。最近Zhu等[13- 14]、Wang等[15]、吴玉等[16]发展了基于这种解耦的分数阶方程的逆时偏移方法,并取得了理想的偏移结果。

然而解耦的分数阶方程在数值求解时存在一定的困难。因为该分数阶的阶数是与Q 有关的随空间变化的的函数，所以分数阶的拉普拉斯算子就是关于波数和空间的混合算子，而这个混合域的算子是不能直接进行傅里叶反变换的。Zhu在计算时，简单地采用了将Q 取空间平均值的办法，这种做法虽然解决了上述问题，但无疑会引入计算误差。为了精确求解该分数阶波动方程，Li等[17]使用最小二乘拟合的办法来近似随空间变化的拉普拉斯算子；Chen 等[18]、 Sun等[19]使用Low-rank近似方法将混合算子分离并获得了较高的波场模拟精度；Yao等[20]采用Hermite分布近似函数(HDAF)来计算分数阶拉普拉斯算子。HDAF方法计算时不依赖于FFT，因此相比于其他方法，具有更好的灵活性和并行性。此外，Chen等还提出了一种基于泰勒展开的方法来近似混合域的拉普拉斯算子，相比于其他方法，这种基于泰勒展开的方法计算量更小，同时具有较高的波场模拟精度。

上述方法都基于分数阶粘滞声波方程，弹性波方程相比于声波方程计算更复杂，但其包含的波场信息也更为丰富。因此，发展分数阶解耦的粘滞弹性波动方程是具有实际意义的。笔者从Zhu所提出的分数阶弹性波动方程出发，利用泰勒近似的方法分离混合域的拉普拉斯算子，最终推导出一组阶数固定的分数阶粘滞弹性波动方程，从而有效地避免了变分数阶拉普拉斯算子的混合域求解问题。本文所提出的方程只包含速度-应力两个分量，而Zhu的方程包含应力-应变-速度三个分量。因此我们的方程计算简单，内存占用小更小。另外，笔者采用伪谱法计算分数阶拉普拉斯算子，可以有效地避免空间频散问题。

1 方法原理

一阶速度-应力分数阶粘滞弹性波动方程可以表示为:

(1)

且

(2)

式中：σxx、σxz和σzz表示应力分量;vx和vz表示速度分量;ρ代表密度;cop以及cos为在参考频率为ω0时的P波和 S波的速度;fx和fz为震源的水平和垂直分量。ηp与ηs为方程中主控相位的参数，τp与τs为方程中主控振幅的参数；Mo和μo分别为P波和 S波的模量；γp,s为分数阶的阶数，由其定义可知，γp,s与空间可变的品质因子Qp,s有关。对于非均匀介质，Zhu采用的取平均值确定γp,s的做法显然会引入计算误差。下面将介绍一种更加准确的近似方法。

以式(1)中的第三项为例，假设介质均匀，使用广义傅里叶伪谱方法将其变换到波数域，可以得到如下表达式：

(3)

且

(4)

(5)

(6)

类似地，方程(1)中的其他项也可以运用上述的方法进行近似，最终推导出新的解耦的分数阶粘滞弹性波方程：

(7)

以及，

(8)

其中，

(9)

对比方程(1)与方程(7)和方程(8)可以发现，我们推导的新的分数阶方程的阶数是常数，从而避免了阶数随空间变化的原方程的求解误差。

2 实验结果

2.1 层状模型

首先考虑一个简单的双层介质模型，模型的介质参数如下：上层介质纵波速度为3 000 m/s，横波速度为2 200 m/s，纵横波的品质因子都取为20；下层介质纵波速度为4 000 m/s，横波速度为3 200 m/s，纵横波的品质因子都取为100。理论上，横波的品质因子应该小于纵波品质因子，但为了简便，设置为二者相同。模型尺寸为400×400个网格点，空间采样间隔为10 m，时间步长为1 ms。水平界面位于深度为2 200 m处。主频为25 Hz的雷克子波震源位于模型中心。对于该模型，我们采用三种方法进行求解：Zhu的平均求解方法(方程(1))，本文提出的方法(方程(7)和方程(8))，以及Point-wise求解方法。Point-wise方法即为对该层状模型的上、下两部分分别求取分数阶的拉普拉斯算子，该方法是精确的，可作为参考解。但是由于该方法需要对每一块参数均匀的区域进行计算，当模型复杂时，Point-wise求解方法将面临巨大的计算开销。图1展示的波场快照中：图1(a)是采用Zhu的平均方法，图1(b)是采用本文提出的方法，图1(c)是采用Point-wise方法，图1(d)是图1(a)与图1(c)的差异，图1(e)是图1(b)与图1(c)的差异，图1(f)是使用完全弹性波动方程模拟该模型的结果。所有图片都显示在同一色标范围内。由图1可见：①Zhu的平均求解方法与参考解之间存在较大误差，本文提出的方法与参考解之间的误差则很小；②相比于完全弹性方程，粘滞波动方程模拟的振幅存在明显衰减。

为了进一步对比，图2显示了距离震源1 km处的单道地震记录。由图2可见：①相比于完全弹性方程，吸收衰减介质不仅造成了地震波振幅的衰减，同时也造成了相位的畸变；②与参考解相比，Zhu的平均求解方法同时造成振幅和相位都存在误差，而本文提出的方法则很好地匹配了参考解。

图1 不同方法计算得到的波场快照及其差值Fig.1 Snapshots computed by different methods and their residuals(a)平均方法；(b)本文方法；(c)参考解；(d)为(a)与(c)的差值；(e)为(b)与(c)的差值；(f)完全弹性介质的波场快照

图2 双层介质1 km处单道地震记录对比Fig.2 Traces computed by different methods of two-layer model at 1 km

图3 BP模型Fig.3 BP model(a)P波速度；(b)P波品质因子

2.2 复杂模型

图3展示的是经典的BP气云模型，图3(a)为P波速度Vp模型，S波速度由Vs=Vp/1.73 换算得到；图3(b)为P波品质因子Qp模型，S波品质因子由Qs=Qp/1.3 换算得到。采用主频为20 Hz的雷克子波震源，于(x,z)=(2 000 m，10 m)处激发，空间采样间隔为10 m，时间步长为1 ms。由于该模型比较复杂，上述的Point-wise方法则不再使用，这里使用Low-rank分解方法[21]直接对方程(1)进行计算，并取较小的时间步长0.2 ms，将其模拟结果作为参考解。

图4中：图4(a)为Zhu的平均求解方法(方程(1))，图4(b)为本文提出的方法，以及图4(c)Low-rank方法模拟的单炮记录。图5则是对应的(x,z)=(1 500 m，10 m)处的地震记录，由图4、图5可以发现，对于非均匀介质，平均方法求解的地震记录与参考解之间存在较为明显的误差；而本文提出的新的解耦分数阶粘滞波动方程的模拟结果几乎与参考解完全吻合。

图4 不同方法计算得到的单炮记录Fig.4 The common-shot gathers calculated by different methods(a)平均方法；(b)本文方法；(c)Low-rank方法(参考解)

图5 BP模型1 km处单道地震记录对比Fig.5 Traces computed by different methods of BP model at 1 km

3 结论

笔者在基于常Q模型的解耦分数阶粘滞弹性波动方程基础上，推导了更加适合于非均匀介质模拟的固定分数阶的解耦粘滞弹性波方程。该方程中的拉普拉斯算子的阶数是固定的，从而避免了由于平均阶数方法引入的数值误差。通过数值算例已经说明，无论对于简单的层状模型还是复杂的BP模型，本文推导的新方程都具有较好的精度。同时，本文的一阶方程中，只包含速度-应力两个分量，因此占用更小的内存开销。该方程利用交错网格伪谱法求解，可以有效地避免空间数值频散。本文发展的固定分数阶的解耦粘滞弹性波方程，可用于粘弹性介质的衰减补偿逆时偏移。