0 引言

波形反演作为重建地下速度结构的重要方法受到学界广泛关注.地震波数值模拟是认识地震波在地下复杂介质中传播特征的重要手段，也是波形反演的基础.当前，前人已经提出了多种地震波正演方法，包括有限差分法、有限元法、伪谱法、体积元法等方法.其中有限差分方法因其简单、灵活以及通用性强，容易在计算机上实现等特点，得到了深入研究和广泛应(Alterman and Karal, 1968; Virieux, 1986; Levander, 1988; Zhang, 1997, 1999；Yang et al., 2002, 2003; Operto et al., 2007, 2009；吴国忱等，2007；李桂花等，2011；Liu and Sen, 2011，2013).在地震波场数值模拟研究中，按计算域的不同可以分为时间域正演和频率域正演.时间域正演求解较为直接，一般需要迭代求解，计算量较大，也容易造成误差的累积.相比之下，频率域正演由于其各频率计算是相互独立的，易于并行计算，特别是对于多震源激发的情况可以大大缩短正演计算时间.另外，频率域正演在单频计算过程中数值误差没有累积效应.

随着全波形反演研究的兴起，频率域正演得到了重要应用.全波形反演的一个重要难点在于其巨大的计算量，而频率域正演易于并行计算的特点可以显著缩短计算时间.频率域正演的另一个重要性在于，可以服务于灵活的分频多尺度反演，从而提高反演的稳定性和可靠性.频率域波形正演作为频率域波形反演的重要组成部分，前人对其做了许多研究(Jo et al., 1996；Shin and Shon, 1998；Štekl and Pratt, 1998；Hustedt et al., 2004；Chen, 2012, 2013; Min et al., 2012；Chen and Cao, 2016).当然其也有一定的不足之处，例如，需要巨大的存储量.但总的来说，频率域正演作为一种有效的波场模拟方法具有十分重要的研究价值.

有限差分方法作为频率域波场模拟中的重要方法，有两类方法得到了广泛应用，一类是旋转法(Jo et al., 1996；Shin and Shon, 1998；Štekl and Pratt, 1998; 曹书红和陈景波, 2012)，其基本原理是将水平方向上的差分算子旋转到对角线上，其主要思想是利用对角线上点来构建差分算子，然后对这两种算子进行加权平均获得更高的数值精度.另一类是平均导数算法(Chen, 2012, 2013; 唐祥德等, 2015；张衡等, 2015；Chen and Cao, 2016)，其基本原理是利用单一水平方向和垂直方向上不同层的有限差分算子进行加权平均然后获得更高的精度.传统的旋转法仅适用于正方形网格(Chen, 2012, 2013), 因为其将Laplace算子看做一个整体进行差分，然后再利用坐标旋转的原理构建对角线上的算子，当网格间距不等时，其旋转方式不成立.所以这种方法虽然提高了有限差分格式的数值精度，但限制了该方法的使用范.因此Chen(2012)提出了平均导数算法，但此种方法也有其不足之处，例如其失去了旋转法所具有的良好的几何旋转性质(Chen，2013)，此外纯粹的平均导数算法并不能离散具有混合偏导数的波动方程，因为离散混合偏导数需要用到对角线上的点，而平均导数算法无法对对角线上的点进行平均.例如，弹性波方程.若要求解此类方程，其需要引入旋转法的思想(Chen and Cao, 2016).而纯粹的旋转法却可以求解此类方程(Štelk and Pratt, 1998).因此，深入研究旋转法对于求解各类复杂的波动方程具有十分重要的意义.

本文在前人旋转法的基础上，将Laplace算子分开成x方向二阶偏导数和z方向二阶偏导数来看，先用坐标轴方向上的点构建二阶差分算子，然后利用旋转的思想分别将其旋转到对角线上，构建出新的差分算子.为了获得更高的精度，本文方法采用17点格式，其可以看成两个九点格式的组合，对于二阶偏导数单个九点格式可以构造出三个差分算子，内外层两个九点格式共用中间格点，总共可以构造出6个差分算子，然后对其加权平均，其具体形式后文将会给出.本文将17点差分格式拆分成两个九点格式分开来看的思想具有十分重要的意义，内外两个九点，相互独立但又统一，基于此，无论是传统旋转法九点格式，平均导数方法九点格式，均可以用本文的方法和思想推广到双九点格式，此结论将在附录中借助平均导数算法(Chen, 2012)予以证明.最为重要的是其可以简单的借助传统九点法构造的思路，推广到17点格式用于高精度求解不同的波动方程.由于传统九点法数值精度一般为最小波长网格点数即为4个左右(Jo et al., 1996; Štekl and Pratt, 1998; Chen, 2012, 2013)，对于近地表低速地区，其差分网格间距需要取较小值才能满足相应的精度要求，而17点格式大大提高了相应的数值精度，虽然对于相同网格大小的模型其计算量和存储量有所增加，但其使得差分网格间距的选择范围大大提高.本文双九点格式的提出将为后人借助传统九点法求解各种复杂波动方程提供简单易懂的思路.

传统的观念认为，四阶精度所得的波场模拟效果要好于二阶精度，在时间域中一般来说确实如此，但在频率域中情况则大不相同，因为其不仅是一个提高单个差分算子的问题，更是一个对目标函数进行数值优化的问题.对于单个差分算子的优化问题，前人四阶精度对应的优化系数确实为最优，但由于对方程进行优化时，其多了自由项P的优化，优化问题不一样，导致四阶精度差分算子的优化系数不一定能使目标函数相速度误差最小，其次四阶精度差分算子是由二阶差分算子优化得来的，其仅为二阶差分算子取特殊值时的特殊情况.因此本文在前人研究基础之上基于广义旋转法提出了一种全新的二阶精度17点格式，传统四阶精度17点格式(曹书红和陈景波, 2012)仅为本文格式的特殊形式.本文格式在保留传统旋转法良好对称性的同时成功将有限差分网格推广到矩形区域，和基于平均导数算法的四阶精度17点格式(唐祥德等, 2015)相比，本文方法数值精度更高，以上结论下文将一一证明.

1 差分格式基本原理

已知2D频率域声波方程为

(1)

其中P是波场压力，w是角频率，v(x, z)是地下介质声波速度.本文双九点格式基本思想是将传统17点格式拆分成内外两个九点格式，并分别用对角线上的点构造精度为二阶的有限差分算子，再将其进行加权平均.

对于x方向二阶偏导数有：

(2)

对于z方向二阶偏导数有：

(3)

对于自由项的近似，可根据各点距P_{m, n}点的距离来进行加权近似，公式为

(4)

其中优化系数之间存在关系为

(5)

由式(3)、(4)可以看出，用二阶精度来构造17点格式时可以离散出6个二阶差分项，与前文所述一致.当x与z方向间距相等时，x方向与z方向有限差分系数应该有相等关系.

现在由式(2)、(3)、(4)经过两次简单的赋值来证明前人四阶精度17点格式(曹书红和陈景波, 2012)仅为本文差分格式特殊形式的特殊形式.

对于原式中的差分系数经过简单的赋值，可得到广义四阶精度17点格式，然后当数值差分网格间距相等时，其等价于四阶精度17点格式(曹书红和陈景波, 2012).

当原式中差分系数有如下关系时：

x方向与z方向二阶偏导数变为：

(6)

(7)

式(6)、(7)即为一次近似所得到的广义四阶精度17点格式，由公式的推导可知其仅为本文格式变量取特殊值时的特殊情况.接下来将在等网格间距情形下进行二次近似来证明传统四阶精度17点格式仅为本文格式特殊情况的特殊情况.

当网格间距Δx=Δz=Δ时，由对称性x方向与z方向优化系数有相等关系，即a=b; 则式(6)和(7)化简可得：

(8)

式(8)即为传统四阶精度17点格式(曹书红和陈景波, 2012).

此外本文双九点格式当相应的优化系数a_2i=b_2i=0，与外层九点格式相对应的点系数也为零时，其可以退化为一般的九点格式.而对于差分算子结构已经固定的四阶精度17点格式而言不具备此种性质.

2 优化系数的求取和数值频散分析

在这一部分，对本文提出的差分格式进行优化系数的求取和数值频散分析.

将P(x, z, w)=P₀e^{-i(k_xx+k_zz)}代入(2)、(3)、(4)式整理可得归一化的相速度为

(9)

其中:

其中：，λ为波长，优化系数a_ij, b_ij, c, d₁, d₂, e₁, e₂, f₁可以由目标函数最小化求得，公式为

(10)

其中，θ取值范围为[0，π/2].当r>1时，其中r=Δx/Δz为网格比，1/G_x取值范围为[0, 0.45]，1/G_z=1/(rG_x).当r < 1时, 1/G_z的取值范围为[0, 0.45]，1/G_x=r/G_z.表 1为不同网格比情况下的优化系数.

从表 1可以看出，当r=1时，即x和z方向有限差分网格间距相等，格式无方向差异性，x方向与z方向相应的优化系数具有近似相等的关系，这一方面符合事实规律，另一方面验证了前文所述.为了说明本文优化方法的正确性和高精度，下面比较本文方法与传统未经优化的5点法所得相速度频散曲线.

由对称性，当r＜1时，将r＞1时x与z方向的优化系数互换即可得到相应的优化系数.图 1为传统5点格式与本文二阶精度17点格式频散曲线对比图.从图 1可以看出，为保证数值误差在1%以内，本文格式在有限差分网格间距相等情形下，需要2.4个，在网格间距不相等的情形下数值精度进一步提高，需要2.2个左右即可.从上面可以得出，为满足相同的精度要求，本文所需网格点数较少，本文方法和前人方法相比数值精度更高.从图 1中还可看出本文方法对数值频散抑制性能较好.

3 数值算例

在这一部分，将给出不同的速度模型来验证本文提出算法的正确性.首先在均匀介质模型不同网格比下用本文算法的结果与解析解相比较，然后再在Marmousi模型下进行数值运算，验证单频切片效果.本文所有结果均是采用完全匹配层(PML)吸收边界条件得出来的.

为了充分说明本文方法的正确性和高精度，将借鉴Liu(2014)文中的比较方法，设置如下图 2激发接收系统，炮点位于(21，21)点处，接收点分别位于坐标点(21，1)、(51，1)、(81，1)处，分别表示为点1、点2和点3.地震子波假设为雷克子波，其子波主频分别取10 Hz、15 Hz和20 Hz，分别编号No.0、No.1、No.2, 并在上述三个不同的地方对其进行观测然后比较结果的差异性.当网格间距相等即r=1时, 取dx=dz=18, 当网格间距不等时，选取r=2情形进行验证，此时取dx=18 m，dz=9 m.

对于均匀介质模型解析解可以由式(11)获得(Alford et al., 1974；Chen，2012)，公式为

(11)

其中F^-1和F表示关于时间t的傅里叶和反傅里叶变换. ，这里(x₀, z₀)表示震源位置，(x, z)表示接收点位置.f(t)表示雷克子波，H₀⁽²⁾表示零阶第二类Hanker函数.

从图 4—图 6可以看出，在子波主频较低情况下，本文方法与前人17点方法所得地震记录随偏移距的增大依然均能较好的与解析解重合.但随着偏移距和子波主频的增加，前人17点方法所得地震记录与解析解表现出较大的差异，而本文方法依然能较好的与解析解重合.

接着在相同的观测系统中，对网格间距不等情形将本文方法所得结果与解析解进行对比验证本文方法的正确性.图 7—图 9是网格间距比为2的情形下所得到的本文结果与解析解的对比图，从图中可以看出本文结果与解析解很好的重合，从而说明本文方法的正确性.

最后对一个实际的Marmousi模型进行分析，模型大小为200×250，炮点位于坐标点(20，80)处，差分网格间距为dx=10 m, dz=8 m.此时网格间距比r=1.25，对应的优化系数如表 1所示，分别在频率为15 Hz和30 Hz时验证单频切片数值频散特性和边界吸收效果，从图 11可以看出单频切片没有明显的数值频散效应.

4 结论

本文基于二阶精度广义旋转法提出了一种全新的双九点有限差分格式用于求解频率域波动方程，本文方法不仅拥有十分良好的几何旋转性质而且为保证数值精度在1%以内，单位波长网格点数需要2.2个左右，前人四阶精度17点格式为本文方法等网格间距情形下的特殊情况.本文方法数值精度和差分格式的频散抑制性能均要优于前人平均导数方法所得17点格式.此外，本文差分方法和将17点格式拆分成两个相互独立而又统一九点格式的思想对于在矩形网格中求解各类复杂的波动方程，以及提高数值精度，均具有十分重要的意义.

附录A

将原文中的方程加入PML吸收边界条件后可以得到：

(A1)

其中：

, a为常值，一般取1.79，L表示PML层的厚度，f₀为子波主频.

将式(A1)带入原式结合本文构造的差分格式整理可得：

(A2)

其中：

附录B

双九点格式在平均导数算法中的应用，本文以传统标量波方程求解做为实例来说明.

对于内层九点有：

(B1)

其中：

对于外层九点格式有：

(B2)

其中：

将外层九点格式与内层九点格式组合在一起即为双九点格式在平均导数算法中的应用，其最终表示形式为

(B3)

当a=1, b=0时，即为平均导数算法九点格式(Chen, 2012).本文格式也是基于此原理得出来的.

附录C

对于二维各向同性介质中弹性波方程有：

(C1)

其中：

同理，对u可以进行同样形式的差分，以上均为考虑矩形差分网格给出，同理通过简单的赋零值便可得到相应的九点格式.