1 引言

逆时偏移是一种基于全波方程的地震偏移方法，由于全波方程能够更加精确地描述地震波在地球介质中各个方向的传播，所以相比于Kirchhoff偏移、射线束偏移和单程波偏移等传统偏移技术，逆时偏移的成像精度大为提高. Whitmore(1983)首次提出逆时偏移方法，同年，Baysal等(1983)详细阐述了逆时偏移的概念和方法，利用伪谱法实现了逆时偏移，McMechan(1983)利用有限差分方法实现了逆时偏移.最初上述各方法用于叠后偏移.Chang和McMechan(1986，1987，1990，1994)先后采用有限差分方法实现了声波、弹性波的叠前逆时偏移；Wu等(1996)分析了高阶有限差分叠前逆时偏移的成像精度；Yoon等(2003)采用高阶有限差分方法和伪谱法实现了三维叠前逆时偏移；Xu等(2010)利用频率域外推方法实现了叠前逆时偏移；Duveneck和Bakker(2011)和Zhang等(2011)利用有限差分方法均实现了各向异性声波叠前逆时偏移.在国内，底青云和王妙月(1997)开展了弹性波有限元逆时偏移研究；杜启振和秦童(2009)、刘红伟等(2010)、严红勇和刘洋(2013)先后利用高阶有限差分方法实现了叠前逆时偏移；Liu等(2011)提出了利用多次波进行逆时偏移.另外，Chattopadhyay和McMeehan(2008)对几种成像条件进行了对比分析；Zhang和Sun(2009)提出了利用拉普拉斯算子滤波方法来压制低频噪声.

求解全波方程是逆时偏移的关键步骤之一，运用快速、高精度的数值算法是逆时偏移中一个非常重要的问题. 有限差分方法因其编程易实现、相同网格点数情况下存储量小、计算速度快、并行性高等优点，在地球物理学中得到广泛应用.但传统的差分格式会因为数值波速与真实波速的不同步，在粗网格情况下产生假波，引起严重的数值频散(Moczo et al., 2000). 为了消除数值频散，目前通常采用的方法有三种: 一是加细网格，二是构造高阶数值格式(Dablain，1986)，三是对差分系数进行优化(Liu，2013； Zhang and Yao， 2013). 加细网格会增加计算量和存储量，这种做法不现实. 而高阶数值格式虽然具有高阶精度，但不一定具有低数值频散，且高精度导致长的差分算子，不利于边界处理和并行的实现. 对差分算子进行优化，能在一定程度上有效降低数值频散. 除了上述三种方法之外，Yang 等(2003)提出一种新的方法，将一种近似解析离散化算子(Kondoh et al., 1994)引入到地球物理领域. 该算子基于截断的泰勒展式以构造空间局部插值函数，并通过局部网格点之间的连接关系，用位移及梯度的线性组合来逼近空间高阶偏导数，包含了更多的波场信息，进而能有效压制数值频散. 在使用该离散算子的同时，结合不同的时间离散方法，可以得到具有低数值频散的数值格式，如NADM(Yang et al., 2003)、ONADM(Yang et al., 2006，2007a)、INADM(Yang et al., 2007b)、WNADM(王磊等，2009)等基于泰勒公式展开进行时间离散的数值方法，以及同时使用了三阶Runge-Kutta和算子级数方法的IRK-DSM方法(Yang et al., 2009)，这些方法的共同特点是均具有良好的压制数值频散的能力，计算效率高.另一方面，波动方程可以表述为哈密尔顿系统的形式，且其时间上具有辛结构. 冯康和秦孟兆(2003)认为，对于哈密尔顿系统的离散格式应当采用保辛的数值格式. 保辛格式具有计算快速稳健、能长时间跟踪等优点(李一琼等，2011；Li et al., 2011，2012). 为了保证在求解哈密尔顿系统时离散方法的保辛性，马啸等(2010)基于声波方程的哈密尔顿表述，时间上采用二阶显式的辛可分Runge-Kutta III A-III B格式，空间上采用近似解析离散的方法，得到了一种时间精度为二阶、空间精度为四阶的新型辛格式(NSPRK).

本文的目的主要是基于四阶NSPRK数值格式的思想，发展新的具有六阶空间精度的格式，以提高现在NSPRK方法的精度，在此基础上进一步研究新格式稳定性条件和数值频散关系，并进行计算效率比较和波场模拟，然后和声波叠前逆时偏移理论相结合，给出基于六阶NSPRK方法的声波叠前逆时偏移方法，最后利用该方法对Sigsbee2B模型进行了偏移成像. 基本步骤是，在正演过程中采用六阶NSPRK方法，而边界处采用PML吸收边界条件(Komatitsch and Tromp， 2003)，并利用归一化互相关成像条件进行成像，最后利用拉普拉斯算子滤波(Zhang and Sun， 2009)，进而得到成像结果，同时和四阶NSPRK方法以及传统的有限差分方法如六阶差分方法、四阶Lax-Wendroff correction方法(简称LWC)(Dablain，1986)的偏移结果进行比较.数值结果表明，基于六阶NSPRK方法的声波叠前逆时偏移能得到更好的偏移结果，是一种优于四阶NSPRK方法、六阶传统差分方法、四阶LWC方法的叠前逆时偏移方法的新算法，尤其是在粗网格情况下具有更加明显的优越性. 2 正演算法 2.1 保辛格式

二维声波方程可写为：

其中，偏微分算子若令v=，则可将方程(1)写成如下线性可分的哈密尔顿系统:

针对可分的哈密尔顿系统(2)，时间上采用孙耿(1997)提出的二阶对称辛可分Runge-Kutta III A-III B格式，得到:

其中uⁿ、vⁿ分别表示第n时间层的位移和速度，该格式具有保辛的性质.

令，则方程(2)可变成:

其中.针对哈密尔顿系统(4)，采用近似解析离散化方法(Yang et al., 2003)离散方程(4)中的高阶空间偏导数，于是得到时间精度二阶、空间精度四阶的NSPRK(马啸等，2010)正演数值算法. 为了提高保辛算法的空间精度，我们利用四阶近似解析离散化算子的基本思想(Yang et al., 2003，2012)，发展了具有六阶空间精度的计算公式，然后参照马啸等(2010)的做法，将时间保辛格式与六阶近似解析离散化算子相结合，于是得到时间精度二阶、空间精度六阶的NSPRK算法. 由于这种算法在时间上采用了辛格式，因此它在时间上具有保辛性，另一方面，由于在空间上采用六阶近似解析离散化方法，因此它具有低数值频散的特点. 在后面，我们将通过数值算例来说明这一点.

参照原四阶近似解析离散化算子(Yang et al., 2003，2012)的基本思想，本文详细推导了具有六阶精度的高阶空间偏导数计算公式，为了方便起见，以为例，新的六阶精度计算公式如下：

其中Δx，Δz分别是空间两个方向的网格步长，u为位移. 用v代替u可得关于v的六阶精度高阶空间 偏导数的近似计算公式，u和v其余的六阶精度高 阶空间偏导数的近似计算公式及推导过程详见附录A.

的四阶精度近似解析离散化方法(Yang et al., 2003)计算公式如下:

u和v其余的高阶空间偏导数的近似计算公式见参考文献Yang 等(2003).

2.2 传统差分格式

对于传统的有限差分方法，在离散高阶空间偏导数时，只利用函数值.

例如四阶LWC方法(Dablain，1986)，首先通过泰勒展开和时空转换得到波动方程的差分格式如下：

然后对格式(7)中的空间偏导数离散，其中二阶空间偏导数以为例，其四阶精度的计算公式如下：

四阶空间偏导数通过对二阶空间偏导数进行中心差分计算(Dablain，1986)，以为例，计算公式如下：

对于时间精度二阶、空间精度六阶的传统差分方法(以下简称六阶CFD)，其格式如下：

六阶CFD在离散高阶空间偏导数时，只利用函数值，以为例，其计算公式如下：

2.3 稳定性条件

对于显式的差分方法，需要满足一定的稳定性条件. 仿照Yang等(2006)的稳定性分析方法，可得到本文六阶NSPRK方法在二维情况下满足的稳定性条件如下：

其中h=Δx=Δz为网格步长，C为库郎数. 稳定性条件(12)的推导过程详见附录B.

2.4 数值频散分析

数值算法会带来数值波速和真实波速之间的差距，称之为数值频散. 数值频散的大小是衡量数值算法好坏的重要指标. 本节对我们发展的六阶NSPRK方法进行数值频散分析，并和原来的四阶NSPRK方法、四阶LWC方法、六阶CFD方法进行对比.

首先定义数值频散比率如下(Moczo et al., 2000)：

其中c_num和c₀分别是数值波速和真实波速，则R和1越接近，表明数值方法的数值频散越小，反之，数值频散越严重. 定义采样率，其中h和λ^grid分别为网格步长和数值波长. 根据Nyquist-Shannon采样定理，Sp≤0.5才有意义(Dablain，1986). 经过简单推导，我们有，其中w_num为数值角频率，Δt为时间步长，C为库郎数，因此，只要推出w_numΔt的表达式，就可以得到数值频散关系. 对于本文发展的六阶NSPRK方法，其数值频散关系及推导过程详见附录B.
图 1给出了库郎数C=0.3时，在不同方向的六阶NSPRK方法、四阶NSPRK方法、六阶CFD方法、四阶LWC方法所对应的数值频散关系. 由图 1可以看出，六阶NSPRK方法、四阶NSPRK方法、六阶CFD方法、四阶LWC方法的最大数值 频散误差分别为：1.61%、7.02%、18.19%、25.82%. 这表明，相比于传统的差分方法，六阶NSPRK方法和四阶NSPRK方法由于采用了近似解析离散化算子离散空间高阶偏导数，均能有效压制数值频散，特别地，在这四种方法中，六阶NSPRK方法的数值频散最小. 另一方面，由图 1可以看出，方向不同，四种方法的数值频散误差也不同，这种现象称为数值频散各向异性，其中六阶NSPRK方法的数值频散各向异性最小. 2.5 计算效率比较

下面研究我们发展的六阶NSPRK方法的计算效率. 为此，选取如下的波动方程：

其中(x_s，z_s)表示震源位置. 选取计算区域0≤x≤18 km，0≤z≤18 km，声波速度c₀=4.0 km·s^-1，计算时库郎数C=c₀Δt/Δx=0.32，震源位于计算区域的中心，震源是主频为f₀=20 Hz的Ricker子波，其随时间变化的函数为：

图 2给出了T=1.8 s时刻的波场快照，其中图 2a、图 2b、图 2c、图 2d分别由六阶NSPRK方法、四阶NSPRK方法、六阶CFD方法、四阶LWC方法这四种方法在粗网格(Δx=Δz=75 m)情况下计算所得. 对比这四个波场快照图，可以看出，在相同条件下，六阶CFD方法和四阶LWC方法有严重的数值频散，四阶NSPRK方法数值频散较弱，而六阶NSPRK方法没有可见的数值频散. 由此可见，相比于六阶CFD方法和四阶LWC方法这两种传统的差分方法，六阶NSPRK方法和四阶NSPRK方法能有效地压制数值频散，同时，六阶NSPRK方法比四阶NSPRK方法压制数值频散的能力更强.

对于四阶NSPRK方法、六阶CFD方法和四阶LWC方法，为了消除数值频散，需要加细网格. 当网格步长分别加细到Δx=Δz=50 m、Δx=Δz=25 m、Δx=Δz=18 m时，四阶NSPRK方法、六阶CFD方法和四阶LWC方法即可消除数值频散，得到清晰的波场快照(图 3a、图 3b和图 3c).比较图 2a和图 3a、图 3b、图 3c，它们都没有可见的数值频散，但是计算得到这些波场快照所需要的CPU时 间以及存储量是不一样的. 表 1给出了六阶NSPRK方法、 四阶NSPRK方法、六阶CFD方法和四阶LWC方法分别计算得到图 2a、图 3a、图 3b、图 3c所需要的计算时间和存储量. 由表 1知，在消除数值频散方面达到相同效果的情况下，六阶NSPRK方法需要的计算时间是7 s，而四阶NSPRK方法、六阶CFD方法、四阶LWC方法所需的计算时间分别为16 s、28 s、173 s，这表明，六阶NSPRK方法的计算速度大约分别是四阶NSPRK方法的2.3倍、六阶CFD方法的4.0倍、四阶LWC方法的24.7倍. 从存储量上看，六阶NSPRK方法需要241×241的12个存储数组，而四阶NSPRK方法、六阶CFD方法、四阶LWC方法分别需要361×361的12个、721×721的3个、1001×1001的3个存储数组. 这表明，六阶NSPRK方法所需要的存储量大约分别是四阶 NSPRK方法、六阶CFD方法、四阶LWC方法的 44.6%、44.7%、23.2%.因此，相比于四阶NSPRK方法、 六阶CFD方法和四阶LWC方法，六阶NSPRK方法是一种计算速度快、低存储、能有效压制数值频 散的保辛差分格式. 同时，六阶NSPRK方法和四阶NSPRK方法一样，与传统方法相比，可以提供除了位移之外更多的信息，包括位移的梯度场、粒子速度场、粒子速度的梯度场等信息.

2.6 波场模拟

为了进一步说明六阶NSPRK方法的计算效 果，本节选取了两个不同的模型，分别进行波场模拟. 2.6.1 双层介质模型

在本数值试验中，我们选取双层介质模型，其 上、下层介质的波速分别为4.0 km·s^-1、6.0 km·s^-1，选取计算区域0≤x≤12 km，0≤z≤12 km，速度的间断面位于竖直方向的中心，空间步长为Δx=Δz= 60 m，时间步长为Δt=3 ms，震源位于(6 km，5.1 km)处，震源是主频为f₀=20 Hz的Ricker子波，震源随时间变化的函数如(15)式.
图 4给出了由四种方法计算得到的在T=1.0 s 时刻的波场快照图. 由图 4可以看出，六阶CFD方法和四阶LWC方法具有严重的数值频散，四阶NSPRK方法具有较弱的数值频散，而六阶NSPRK方法没有可见的数值频散. 这表明，六阶NSPRK方法在较粗的网格情况下，对于强间断层状介质中的波场模拟也十分有效.

2.6.2 Sigsbee2B模型

在接下来的试验中，我们选取Sigsbee2B模型，来研究六阶NSPRK方法在非均匀复杂介质中的计 算效果.Sigsbee2B模型的速度结构如图 5所示，它具有很强的非均匀性.

选取空间步长为Δx=Δz=75 ft，时间步长为Δt=2 ms，震源位于计算区域的(37500 ft，25 ft)处，震源是主频为f₀=15 Hz的Ricker子波，震源随时间变化的函数如(15)式. 图 6给出由六阶NSPRK方法计算所得各个时刻的波场快照图. 由这些波场快照图可以看出，六阶NSPRK方法可以计算得到清晰的波场快照图，能有效模拟声波在复杂介质中传播. 图 7为由六阶NSPRK方法计算得到的0 s至9.6 s时刻的地表地震记录，其中图 7a 使用了PML(Perfectly Matched Layer，完全匹配层)吸收边界条件(Komatitsch and Tromp， 2003)，图 7b使用了刚性边界条件. 对比图 7a和图 7b可以看出，六阶NSPRK方法能与PML吸收边界条件有效结合，吸收来自人为边界的反射波.

3 声波叠前逆时偏移理论

叠前逆时偏移实现的基本步骤为：首先，正演计算震源波场，计算时间从T₀时刻开始到T时刻结束，各个时刻的震源波场记为s(t，x，z)；然后反向逆时外推记录波场，计算时间从T时刻开始反推到T₀时刻结束，各个时刻的记录波场记为r(t，x，z)；最后，利用合适的成像条件以及震源波场 s(t，x，z)和记录波场r(t，x，z)，计算得到偏移成像结果，偏移结果记为I(x，z).其中计算震源波场s(t，x，z)和记录波场r(t，x，z)的过程就是数值求解全波方程的过程，该过程利用本文发展的六阶NSPRK方法来计算.

叠前逆时偏移成像条件的选取和计算也是叠前逆时偏移中的核心问题之一，成像条件的好坏会直接影响最终成像结果的质量. 目前叠前逆时偏移中常用的成像条件主要有三种：激发时刻成像条件，互相关成像条件，振幅比成像条件. Chattopadhyay和McMeehan(2008)对几种成像条件进行了对比分 析，认为震源归一化互相关成像条件能够反映地下 界面的反射系数，提供较为准确的成像结果，因此，本文利用震源归一化互相关成像条件进行叠前逆时偏移成像，其表达式如下：

其中的ε表示很小的正数，它的选取是为了使得成像条件表达式的分母不为零，便于计算.

在叠前逆时偏移的过程中往往会产生低频噪声，这些噪声的产生是由潜水波、首波和绕射波等引起的，噪声的存在会影响到成像的结果，所以要得到最终清晰的成像结果，我们需要抑制噪声，对偏移产生的最初成像结果进行去噪处理. Zhang和Sun(2009)提出了利用拉普拉斯算子滤波方法来压制低频噪声，拉普拉斯算子滤波相当于对波场进行角度域衰减，该方法简单易操作被广泛使用. 本文利用拉普拉斯算子滤波方法对成像结果进行去噪处理. 4 声波叠前逆时偏移算例

我们选取被国际广为认可的Sigsbee2B模型进行偏移成像并进行对比研究. Sigsbee2B速度模型如图 5所示. 该模型的合成数据共496炮，其中每一炮最大的道集数为348，最小的道集数为65，而对Sigsbee2B模型进行偏移成像时，选取其中的248炮数据，同时选取两套网格，粗网格 情形下网格步长为Δx=Δz=75 ft，时间步长为Δt=2 ms，细网格情形下网格步长为Δx=Δz=37.5 ft，时间步长为Δt=1 ms. 震源选择主频为f₀=15 Hz的Ricker子波，其随时间变化的函数如(15)式.

图 8给出了在粗网格(Δx=Δz=75 ft)情况下由六阶NSPRK方法、四阶NSPRK方法、六阶CFD方法和四阶LWC方法计算所得的逆时偏移结果.对比图 8，可以看到，在粗网格步长情况下，相比于四阶NSPRK方法、六阶CFD方法和四阶LWC方法，由六阶NSPRK方法进行逆时偏移得到的成像结果更清晰一些，尤其是在四个虚线框标示的区域，这四个小区域包含了盐下的小散射点以及盐丘底部的高速区域. 这是由于粗网格情况下，六阶CFD方法和四阶LWC方法计算波动方程时会有严重数值频散，四阶NSPRK方法具有较弱的数值频散，相比之下，六阶NSPRK方法能有效地压制数值频散. 图 9给出了在细网格(Δx=Δz=37.5 ft)情况下由六阶NSPRK方法、四阶NSPRK方法、六阶CFD方法和四阶LWC方法计算所得的逆时偏移结果.对比图 9，可以看到，在细网格步长情况下，由六阶NSPRK方法、四阶NSPRK方法、六阶CFD方法和四阶LWC方法进行逆时偏移得到的成像结果基本没有差别. 事实上，由六阶NSPRK方法在粗网格情况下的逆时偏移结果和另外三种方法在细网 格情况下的逆时偏移结果，也基本没有差别. 这 也表明，六阶NSPRK方法能有效压制数值频散，可以采用粗网格步长进行偏移成像，从而提高计算效率.

5 结论

本文详细推导了具有六阶精度的近似解析离散化方法，给出了高阶空间偏导数的计算公式，将原四阶NSPRK保辛算法的计算精度提高到了六阶，并对新的六阶NSPRK方法进行了稳定性和数值频散分析，以及计算效率比较和波场模拟，进一步将六阶NSPRK方法应用于RTM，得到一种时间上保辛、同时保持低数值频散的且可以应用大步长进行波场延拓的叠前逆时偏移方法. 数值结果表明，这种六阶NSPRK新方法与原四阶NSPRK方法、传统的六阶差分方法以及四阶LWC方法相比，可以有效压制数值频散，且具有计算速度快、低存储量等优 点. 对Sigsbee2B模型的偏移结果显示，六阶NSPRK 方法可以通过使用较大的空间步长和较大的时间步长来提高计算效率和节约存储量，有望在大尺度波场模拟和逆时偏移成像方面具有广泛应用. 附录A 六阶精度近似解析离散化算子推导

参照Yang 等(2003)的基本思想，首先利用截断的泰勒公式，有

然后以为未知数，求解式(A1)、式(A2)、式(A3)组成的线性方程组，即可得到x方向高阶空间偏导数的表达式，z方向高阶空间偏导数的表达式可类似得到. 由泰勒展式可以看出，得到的计算公式具有六阶精度. 为了方便起见，下面给出本文中所使用到的计算式:

接下来利用方向导数推导混合偏导数的计算公式(Yang et al., 2012). 定义如下两个方向向量 l =(Δx，Δz)T，m =(－Δx，Δz)T，以及标量，我们有关于方向偏导数的式子如下：

由式(A8)、(A9)可得到：

将式(A7)代入式(A10)中，同时将式(A10)的方向偏导数可以看成对应方向的偏导数，有类似于式(A6)、(A7)的表达式，整理后可得<

类似地，有

附录B 六阶保辛算法稳定性及数值频散分析

我们采用Von Neumann的分析方法进行稳定性及数值频散分析(Moczo et al., 2000; Yang et al., 2006). 将声波方程的谐波解代入六阶保辛算法的格式中，即式(3)中，其中为虚数单位，则可以得到如下的迭代关系式：

其中迭代矩阵Ι₃ 为3×3的单位矩阵，B =(b_ij)_3×3的元素如下：

其中h=Δx=Δz为空间步长，d₁=exp(I(2πS_pcosθ))，d₂=exp(－I(2πS_pcosθ))，d₃=exp(I(2πS_psinθ))，d₄=exp(－I(2πS_psinθ))，d₅=exp(2I(2πS_pcosθ))，d₆=exp(－2I(2πS_pcosθ))，d₇=exp(2I(2πS_psinθ))，d₈=exp(－2I(2πS_psinθ)).

为了使得该算法稳定，需要满足ρ(A)≤1，从而可以得到库郎数满足的条件为C=c₀Δt/h≤C_max.同时，有

由式(B11)，可得六阶保辛算法的数值频散关系：