0 引言

地震波波场正演模拟始终是地球物理领域研究的核心内容，寻求高精确的波场正演算子是地震波成像和反演的基础.

有限差分法因其高效性是求解波动方程的常用算法，但是它本质上还是在零频处展开得到的差分系数(Dablain, 1986; Kindelan et al., 1990)，对于高频波的数值模拟，其存在的时间和空间频散限制了这种算法的精度，尤其是在大网格计算的情况下.为了解决有限差分的数值频散问题，许多科学家采用不同的优化算法来提高差分算子的精度(Zhou et al., 2008; Liu and Sen, 2009; Liu et al., 2009; Kosloff et al., 2010; Zhou and Zhang, 2011), 这些方法大都是通过拟合精确的频散关系得到优化的差分系数来提高计算精度和稳定性，但是没有从根本上解决这一问题.

实现时间域波场延拓的另一种方法是谱方法(Fowler et al., 2010; Du et al., 2014).最常用的伪谱法(Kosloff and Baysal, 1982)在波数域求解空间导数，可以避免空间频散，但是由于对于时间微分项是用简单的差分近似，因此无法解决时间频散，且是条件稳定的.理论上，时间域波场延拓可以归结为对一个波数-空间混合域象征进行处理(Wards et al., 2008).关于如何近似该混合域象征，涌现了一批算法：Zhang等(2007)、Zhang Y和Zhang G Q (2009)的单步延拓法利用可分近似的方法，将该混合域象征进行展开，展开过程涉及到特征向量的计算，这种算法与单程波方法类似(Chen and Liu, 2004)，计算量非常大，但是不存在稳定性问题以及数值频散；Soubaras和Zhang(2008)的双步法是利用多项式展开来逼近该混合域象征，然后利用Remez最优化方法来求取逼近系数，这种方法能够使用大时间步长进行波场延拓.快速展开法(Tal-Ezer et al., 1987)是求解波动方程的解析算法，该算法利用Chebyshev多项式来逼近混合域象征，而且没有数值频散；目前，已有科学家(Pestana and Stoffa, 2009, 2010)将这一算法运用于逆时偏移中，取得了较好的效果，但是该方法的缺点在于计算效率太低，原因在于需要进行大量的傅里叶变换.Etgen和Brandsberg-Dahl(2009)提出用拟拉普拉斯算子代替拉普拉斯算子对声波方程进行模拟，补偿因时间差分造成的误差，但是对于速度变化剧烈的速度模型，这种拟解析法在时间和空间上都具有较大的误差.Mast等(2001)提出了波数-空间混合域算法处理超声波在时间上的延拓，Tabei等(2002)将该方法扩展到耦合的一阶速度-应力方程，有效地处理了弱非均匀介质中的波场延拓问题.Liu(1995)最先将波数-空间域方法应用于弹性波，并且基于并矢格林函数解析解，推导了二阶弹性波方程的波数-空间域象征.

在求解波数-空间混合域的谱方法中，二阶时间差分算子的频散误差可以通过波数域中的谱算子进行补偿，这种补偿对于常速介质是精确的，对于一般的变速介质，该混合域象征的计算量非常大，Fomel等(2010)提出了低秩分解来近似该混合域象征算子，这种方法在每一时间步上仅需少量的几次快速傅里叶变换，这是一种计算量与网格点数呈线性关系的波场延拓方法.Song和Fomel(2011)提出低秩有限差分法解决声波方程.Fang等(2014)将低秩有限差分方法应用于交错网格上，提出交错网格低秩有限差分法，并且用来解决一阶声波方程.Han等(2016)将交错网格低秩有限差分法用来解决一阶解耦的速度应力方程.

本文致力于解决解耦的二阶弹性波位移方程，将交错网格低秩方法结合有限差分法，提出交错网格低秩有限差分法，试图改善常规交错网格有限差分法的精度.最后，数值模拟结果证明，我们提出的交错网格低秩方法和交错网格低秩有限差分法具有较高的数值模拟精度.

1 理论 1.1 解耦的弹性波方程

在均匀各向同性完全弹性介质中，弹性波的传播由下述拉梅方程描述(Aki and Richards, 2009)：

(1)

式中，u =(u, v, w)为位移场，λ，μ是拉梅系数，ρ是介质的密度.▽是梯度算子.

对于拉普拉斯算符▽²u可有：

(2)

于是，(1)式还可以写成：

(3)

其中：，分别代表纵、横波速度.

用位移向量的三个分量u，v，w描述的拉梅方程为：

(4)

式中，x =(x, y, z)代表三维空间位矢.

在三维介质中，等价的解耦弹性波方程可以写为(马德堂和朱光明，2003)：

(5)

在二维介质中，可以写为：

(6)

(6) 式中，空间二阶微分算子可以分解成一对一阶微分算子.为了提高模拟精度，我们可以采用交错网格有限差分算子来近似一阶微分算子(陈可洋等，2009)：

(7)

式中，L_x⁺和L_z⁺分别表示沿着x轴和z轴的前向差分算子，L_x^-和L_z^-分别表示沿着x轴和z轴的后向差分算子.本文中各位移分量在交错网格节点配置的方法如图 1所示：u, u_p, u_s定义在主网格节点，w, w_p, w_s定义在半网格节点.

时间微分算子通常由二阶中心差分算子L_tt近似：我们可以给出方程的离散形式：

(8)

利用傅里叶平移变换的性质，可以得到一阶微分算子的一种等价表示：

(9)

其中，F和F ^－1分别代表傅里叶正、反变换，k =(k_x, k_z)表示波数向量.对于地震波场这样的有限带宽信号，对波场应用傅里叶变换可以准确求得波场在空间上的偏导数，可以看到，谱方法计算空间微分，由于没有任何近似，具有很高的精度.然而，我们仍然需要采用中心二阶差分计算时间微分项，受到有限差分固有的稳定性条件限制，对于大时间步长采样，会遇到不稳定问题；并且时间频散很严重(Firouzi et al., 2012).

1.2 拟微分算子与混合域象征

在提出本文的方法之前，先介绍拟微分算子以及它的象征.

拟微分算子可以视为微分算子的推广，我们知道，微分与乘法对偶(齐民友等，2005)，定义拟微分算子如下：

(10)

其中，A称为拟微分算子，而a(x, ξ)称为拟微分算子的象征. 是u(x)的傅里叶变换，n代表空间的维数，本文中，n=2.

1.2.1 低秩分解近似一阶混合域象征

针对大步长采样不稳定问题，我们用波数-空间混合域象征和 (Firouzi et al., 2012)修正谱算子ik_j, j=x, z, 得到了一阶混合域象征，这种谱算子的修正能补偿由于时间项差分离散带来的误差，并且对于均匀介质，这种补偿是精确的.根据这些混合域象征，引入对应的拟微分算子：

(11)

其中，j=p, s.对于不均匀介质，速度随着空间位置而变化，直接采用(10)式计算，计算量是O(N_xlogN_x)，其中N_x是总的采样点数.为了提高计算效率，Fomel等(2013), Engquist和Ying(2009)研究给出了混合域象征的低秩分解，以W_{p, x}=k_xsinc(v_p(x)Δtk/2)为例，它可以有如下的低秩分解：

(12)

其中，w_px(x, k_m)和w_px(x_n, k)分别是从W_{p, x}选取M列和N行组成的子矩阵，M和N是依赖于矩阵的整数，一般是2或者3.利用(11)式，我们可以将拟微分算子A_{x, p}⁺写成：

(13)

同样的分解也可以用于其余三个拟微分算子.可以看到式(12)仅需N次反傅里叶变换，而直接计算A_{x, p}⁺需要N_x次反傅里叶变换，显著地减少了计算量，本文称此方法为交错低秩方法.

1.2.2 低秩有限差分近似拟微分算子

注意到w_px(x_n, k)是一个只与波数有关的矩阵，它可以进一步被分解为(Song et al., 2011)：

(14)

其中，B是L×N_k的矩阵，N_k是总的波数采样点数，L是有限差分算子的长度.它的元素形如sin((2l－1)k_xΔx/2)，C_px(x_n, l)是w_px(x_n, k)和B的广义逆的乘积.

定义新的N_x×L矩阵：

(15)

因此，可以将A_{x, p}⁺u(x, t)写成：

(16)

利用傅里叶平移变换的移位性质，(16)式可以重新写成

(17)

其中，x_R=(x+lΔx, y)，x_L=(x－(l－1)Δx, y)，l=1, 2, …, L.类似地可以推导其余三个拟微分算子的交错低秩有限差分格式.与常规交错有限差分格式相比，差分系数G_px是依赖于空间位置的，它随着模型参数的改变而自适应地变化，因而具有较高的模拟精度.本文称此方法为交错网格低秩有限差分法.接下来用平面波理论分析交错网格低秩有限差分法的数值频散特性，将p=e^i(kx－wt)代入公式(17)，并且将ξ定义为相速度的相对误差，根据频散理论，w= |k| v，以一维为例，可以得到：

(18)

利用(18)式，可以计算交错网格低秩有限差分法的频散特性.ξ如果等于0，表明不存在数值频散，ξ偏离0越大，表明数值频散越严重.

图 2对比了不同差分阶数(4，6，8，16阶)的交错网格低秩有限差分方法和常规交错网格有限差分方法的频散曲线，可以看出两种方法都随着差分阶数的提高频散减弱，注意到常规交错网格有限差分方法随着差分阶数的提高，能够精确拟合频散关系的频率范围并未增加，而对于交错网格低秩有限差分，随着有限差分阶数的增加，其精确拟合频散关系的频率范围变的更宽.

图 3对比了不同介质速度下，常规16阶交错网格有限差分和交错网格低秩有限差分的频散曲线，其中时间采样间隔是Δt=1 ms，空间采样间隔是Δx=10 m.对于一维的16阶交错网格低秩有限差分方法，L=8.图中可以看出，常规的交错网格有限差分只在零频率(波数)附近是精确的，随着频率增大，常规有限差分的频散越来越严重，而交错网格低秩有限差分能够在较宽的波数(大于70%的Nyquist波数)范围内精确地刻画频散关系，使得该方法能够较好地压制数值频散，具有较高的精度.

2 数值实验

本文中数值实验均采用海绵吸收边界条件，震源皆采用主频为20 Hz(最高频率60 Hz)的Ricker子波.

2.1 均匀模型点源测试

首先，采样2D均匀模型，测试本文提出的交错网格低秩方法与交错网格低秩有限差分法，并与常规的交错网格有限差分方法所得结果进行对比.2D模型大小为256×256，X方向和Z方向采样间隔均为10 m.横波和纵波速度分别为1500 m·s^-1、3000 m·s^-1. Ricker子波的主频为20 Hz，时间采样步长为1 ms.震源位于模型中心，并同时产生横波和纵波.根据Song等(2013), 我们仍然采用CFL(Courant-Friedrichs-Lewy condition)数作为稳定性衡量因子，采用作为频散衡量因子，本例中α=0.33，β=2.5.图 4a—4c分别为4阶交错网格有限差分、4阶交错网格低秩有限差分、交错网格低秩方法得到的X分量波场快照，从图中蓝色箭头所指处可以看出：本文提出的方法数值频散更弱，交错网格低秩方法几乎无频散，因而精度更高.图 4d—4f为图 4a—4c在X=1280 m处选取的相应的垂直道.

2.2 层状模型测试

我们用层状弹性介质对比常规交错网格有限差分法和本文提出的方法.采用的模型在X方向和Z方向的空间间隔均为15 m.上层介质V_p=2500 m·s^-1，V_s=1800 m·s^-1.下层介质V_p=3000 m·s^-1，V_s=2000 m·s^-1.分界面位于Z=2100 m处.震源为主频为20 Hz的Ricker子波，位于模型(X, Z)=(2250 m, 1500 m)处.时间采样步长为2 ms.稳定性因子α=0.4，频散因子β=2.图 5a—5c显示了常规6阶交错网格有限差分法、6阶交错网格低秩有限差分法、交错网格低秩方法，X分量波场快照对比.图 5d—5f显示了常规6阶交错网格有限差分法、6阶低秩有限差分法、交错网格低秩方法，Z分量波场快照对比.从图中可以看出6阶交错网格低秩有限差分法和交错网格低秩方法得到的波场快照比常规6阶交错网格有限差分方法得到的波场快照精度更高(数值频散更弱)，如蓝色箭头所示.

2.3 Marmousi模型测试

为了进一步验证本文方法的有效性，我们选择了部分Marmousi模型进行实验.图 6显示了Marmousi的纵波速度模型.横波速度与纵波速度的关系为：模型的大小是376×376，X方向和Z方向网格间距均为8 m.Ricker子波的主频为20 Hz.震源位于模型(X, Z)=(1504 m, 0 m)处，时间采样步长为0.5 ms.检波器设置在模型表面，总的时间记录长度为1 s.

图 7a、7c、7e分别为使用常规交错网格有限差分法、交错网格低秩有限差分法、交错网格低秩方法正演得到的横波X分量地震记录.图 7b、7d、7f分别为使用常规交错网格有限差分法、交错网格低秩有限差分法、交错网格低秩方法正演得到的横波Z分量地震记录.

图 8为图 7中地震记录黑色虚线部分的放大图，可以看出，本文提出的交错网格低秩有限差分法和交错网格低秩方法频散更弱，精度更高.

3 结论与建议

在本文中，我们提出了交错网格低秩有限差分方法以及交错网格低秩方法求解解耦的二阶弹性波位移方程.交错网格低秩有限差分法的差分系数，是通过对波数-空间混合域的算子作低秩近似而得到的，它随着空间位置的变化而变化，相比于泰勒级数在零波数处的展开而得到传统的有限差分系数，低秩有限差分法能在高波数处拟合频散关系曲线；而交错网格低秩方法，对于均匀介质，它的解实质上就是解析解(如果不考虑矩阵低秩近似带来的微弱误差).数值结果表明：本文提出的两种方法相对于常规的交错网格有限差分法具有更高的精度，并且可以适用于复杂模型.本文的方法可以扩展到三维弹性波模型正演模拟，并且进一步应用于弹性波逆时偏移成像.

致谢

感谢审稿专家和编辑提出的宝贵意见.