地震波正演数值模拟方法是研究地震波在地球介质中传播规律的重要手段之一，发展和研究高精度高效率的正演数值模拟技术成为地球物理学中最关注的研究内容.简单且使用频率高的一类地震波即声波在地球介质中的传播过程可视为波动方程的正演模拟过程.目前，常见的地震波数值模拟方法包括：有限元法(Chen，1984; 杨顶辉，2002)、射线追踪法(Cerveny and Firbas，1984; Chapman and Shearer，1989)、谱元法(Huang，1992; Komatitsch and Vilotte，1998)和有限差分法(Kelly et al.，1976; Virieux J，1984; 1986; Dablain，1986; 董良国等，2000)等.

当前，有限差分法因其设计简单、存储量小和计算效率高等优点而被广泛使用.但是该方法因采用差分算子离散空间高阶偏导数而导致地震波场模拟过程中产生严重的数值频散，以致大大降低了计算精度(Fei and Larner，1995; Zhang et al.，1999; Yang et al.，2003).为避免数值频散的发生以有效提高数值模拟精度，Yang等于2003年引入近似解析离散算子以修正差分算子从而提出了一种所谓的近似解析离散(NAD)方法(Yang et al.，2003; 2006a; 王磊等，2009; 宋国杰等，2010; 马啸等，2010).该方法通过位移和位移梯度的线性组合来实现对空间高阶偏导数的离散化，从而有效地压制数值频散以达到提高计算效率.但是该方法也存在相比有限差分而言需要更多的存储量和计算效率不很高等不足(卢明和杨顶辉，2005).为此，Yang等人(Yang et al.，2004;2006b; 卢明和杨顶辉，2005)接着又提出了一种优化的NAD方法即ONAD方法，该方法具有存储量小和计算效率高等优势.然而，该方法使用的近似解析离散算子只有四阶的，即ONAD方法在空间离散上只是四阶精度的.

为更有效地压制数值频散以最大限度地提高模拟精度，本文在前人基础上发展了一种求解二维声波方程的高精度数值模拟算法——八阶ONAD方法.该方法使用八阶精度的近似解析离散算子(Tong et al.，2013)离散空间高阶偏导数，而利用四阶截断的泰勒式离散时间偏导数.数值误差分析、计算效率分析和复杂介质地震波数值模拟结果等显示，比起四阶LAX-Wendroff Corrected(LWC)方法和八阶LWC方法而言，八阶ONAD方法具有数值频散低、计算效率高和波场模拟能力强等优势，是一种非常有效的地震波数值正演模拟方法. 1 八阶ONAD方法

我们以二维声波方程为例来阐述八阶ONAD方法求解波动方程的基本思路.对于均匀介质情形，二维声波方程为

为了求解二维声波方程(1)式，我们令

利用四阶截断的泰勒展开式可得位移向前和向后的计算表达式分别为

这里，Uⁿ=U(nΔt).

将(4)式和(5)式相加，可得：

这样，求解二维声波方程(1)式就转化为计算(7)式.(7)式含有Uⁿ的高阶空间偏导数，本文采用文献[Tong et al.，2013]推导出的具有八阶精度的近似解析离散公式来计算Uⁿ的高阶空间偏导数，具体的计算公式如附录A所示.

我们把用来求解声波动方程(1)式的逼近表达式(7)式称为八阶ONAD方法. 2 误差分析

2.1 理论误差分析

对于求解二维声波动方程的八阶ONAD方法，由截断的泰勒展开式得到时间离散公式(7)式，故其时间误差为O(Δt⁴)；由泰勒展式离散方程(7)式右边的二阶和三阶空间偏导数，其计算公式的空间误差(Tong et al.，2013)为O(Δx⁸+Δz⁸).所以求解二维声波动方程的八阶ONAD方法的误差精度为O(Δt⁴+Δx⁸+Δz⁸). 2.2 数值误差分析

为考察八阶ONAD方法的模拟精度，我们选取八阶ONAD方法、四阶LWC方法和八阶LWC方法求解二维声波动方程初值问题，并对它们的计算结果进行数值误差比较.设二维声波动方程的初值问题为

这里，f₀是频率，c是声波波速，θ₀是声波在t=0时刻的入射角度.易知，该初值问题的解析解为

在进行数值实验时，选取空间网格点数为N=M=201，频率f₀=25 Hz，声波波速c=4 km/s，入射角度θ₀=π/4，设h=Δx=Δz为空间步长.数值解的相对误差定义为

其中，u(t_n，x_i，z_j)为精确解，uⁿ_i，j为数值解.

表 1是八阶ONAD方法、四阶LWC方法和八阶LWC方法在三种不同情形下的最大相对误差，图 1是该三种方法在三种不同情形下的相对误差曲线.从表 1和图 1可知，八阶ONAD方法、四阶LWC方法和八阶LWC方法在三种情形下的相对误差均随积分时间的增加成递增趋势，但相同条件下八阶ONAD方法的相对误差大大地小于四阶LWC方法，也小于八阶LWC方法.

表 1(Table 1)

表 1 不同方法在不同参数情形下位移的最大相对误差比较 Table 1 The comparision of the maximum relative errors on the displacement，respectively computed by different methods for three cases 计算方法 八阶ONAD方法 四阶LWC方法 八阶LWC方法 情形a: h=30 m，Δt=0.5 ms 0.143529 42.4741 0.591174 情形b: h=40 m，Δt=0.8 ms 1.9502 129.7265 5.6062 情形c: h=50 m，Δt=1 ms 13.4508 205.8764 29.1554

表 1 不同方法在不同参数情形下位移的最大相对误差比较 Table 1 The comparision of the maximum relative errors on the displacement，respectively computed by different methods for three cases

图 1

Fig. 1

图 1 八阶ONAD方法与四阶LWC方法、八阶LWC方法在求解初值问题(8)时，在不同网格情况下位移在不同时刻的相对误差 (a)h=30 m，Δt=0.5 ms；(b)h=40 m，Δt=0.8 ms；(c)h=50 m，Δt=1 ms. Fig. 1 The relative errors of the displacement，respectively computed by the eighth-order ONAD method，the fourth-order LWC method，and the eighth-order LWC method，are shown in a semi-log scale for the 2D initial-value problem(8)，corresponding to the three different cases: (a)h=30 m，Δt=0.5 ms；(b)h=40 m，Δt=0.8ms；(c)h=50 m，Δt=1 ms.

本节，我们来研究八阶ONAD方法在均匀各向同性介质中二维声波方程的波场模拟，分别使用八阶ONAD格式、四阶LWC格式和八阶LWC格式进行计算，通过比较验证八阶ONAD方法的计算效率.

在数值实验中，我们取定Courant数α=0.16，声波速度c=4.0 km/s，计算区域为0≤x，z≤10 km，震源位于区域中心，频率f₀=25 Hz，空间步长为Δx=Δz=h，两个接收器分别为R₁(5 km，3 km)和R₂(6 km，3.27 km).接收器R₁在震源上方2 km处，R₂距震源2 km处，且与震源的连线与水平方向成60°的夹角.

图 2a~c分别是八阶ONAD方法、四阶LWC方法和八阶LWC方法模拟得到的粗网格(h=50 m)情况下T=1.0 s时刻的波场快照，图 3a~c分别由这三种方法计算得到的在接收器R₁处的波形记录，图 4a~c分别由这三种方法模拟得到的在接收器R₂处的波形图.由快照和波形图可看到，三种方法模拟得到的数值波场几乎一致，但是四阶LWC方法和八阶LWC方法有明显的数值频散，而八阶ONAD方法模拟得到的波形光滑清晰，没有可见数值频散.显然，在粗网格下八阶ONAD方法在压制数值频散上明显优越于四阶LWC方法，也好于八阶LWC方法.另外，从图形上可以看到，四阶LWC方法和八阶LWC方法在R₁处的波形数值频散很严重，而在R₂处的波形图数值频散较小，这验证了该两种方法的数值频散各向异性比较明显，而八阶ONAD方法几乎没有数值频散各向异性.

图 2

Fig. 2

图 2 粗网格(h=50 m)条件下T=1.0 s时刻的声波场瞬时快照 (a)八阶ONAD方法；(b)四阶LWC方法；(c)八阶LWC方法.Fig. 2 The snapshots of wave fields at time T=1.0 s on the coarse grid(h=50 m)，respectively generated by (a)The eighth-order ONAD;(b)The fourth-order LWC;(c)The eighth-order LWC.

图 3

Fig. 3

图 3 粗网格(h=50 m)条件下在接收点R1(5 km，3 km)处的波形记录 (a)八阶ONAD方法；(b)四阶LWC方法；(c)八阶LWC方法.Fig. 3 The waveforms at the receiver R1(5 km，3 km)on the coarse grid(h=50 m)，respectively generated by (a)The eighth-order ONAD;(b)The fourth-order LWC;(c)The eighth-order LWC.

图 4

Fig. 4

图 4 粗网格(h=50 m)条件下在接收点R2(6 km，3.27 km)处的波形记录 (a)八阶ONAD方法；(b)四阶LWC方法；(c)八阶LWC方法.Fig. 4 The waveforms at the receiver R2(6 km，3.27 km)on the coarse grid(h=50 m)，respectively generated by (a)The eighth-order ONAD;(b)The fourth-order LWC;(c)The eighth-order LWC.

为了分析该三种方法在消除数值频散条件下的计算效率和计算存储量，在Courant数α=0.16下，逐渐减小四阶LWC方法和八阶LWC方法的空间步长.当分别取h=22 m和h=30 m时，四阶LWC方法和八阶LWC方法可以有效消除地震波场模拟中的数值频散，此时波场模拟均不见数值频散.但在Intel(R)Core(TM)2Duo CPU和内存为2 G的计算机环境下，三种方法为消除数值频散所付出的存储量和CPU时间却明显不同.四阶LWC方法和八阶LWC方法为消除数值频散所消耗的CPU时间分别为94 s和67 s，而用八阶ONAD方法模拟得到图 2a所花的时间仅为19 s.这意味着八阶ONAD方法的计算速度约为四阶LWC方法的4.95倍、为八阶LWC方法的3.53倍.在存储量上，四阶LWC方法需要的存储规模为455×455；八阶LWC方法需要的存储规模为334×334；但是八阶ONAD方法需要的存储规模仅为200×200.这意味着八阶ONAD方法需要的存储量是四阶LWC方法的19.32%、是八阶LWC方法的35.86%. 4 数值模拟

本节将通过具体的三个实验来研究八阶ONAD方法在二维声波复杂介质中的数值模拟效果和模拟能力，并同四阶LWC方法和八阶LWC方法进行比较，以体现八阶ONAD方法的优势. 4.1 双层介质模型

首先，考察八阶ONAD方法在双层介质模型中模拟声波传播的效果.模型计算区域为0≤x，z≤20 km，上层速度为c₁=2.4 km/s，下层速度为c₂=4.8 km/s，层间分界面为z=8 km.震源位于(10 km，7 km)，频率为f₀=15 Hz.取空间步长为h=Δx=Δz=50 m，时间步长为Δt=1×10^－3 s.

图 5a~c分别是八阶ONAD方法、四阶LWC方法和八阶LWC方法计算得到在T=2.0 s时刻的波场快照.从图 5a可看出，波在水平面断界面上反射波和透射波的快照图都非常清晰，在间断附近没有可见数值频散，这说明八阶ONAD方法可以很好地模拟层状模型中的波传播规律.由四阶LWC方法和八阶LWC方法分别计算得到的图 5b和图 5c有着明显的数值频散.这进一步揭示了对于层间速度反差较大模型的波场模拟八阶ONAD方法比起其他两种有限差分方法更有优势，有更好的适用能力.

图 5

Fig. 5

图 5 双层介质模型中粗网格(h=50 m)下T=2.0 s时刻的波场瞬时快照(a)八阶ONAD方法；(b)四阶LWC方法；(c)八阶LWC方法.Fig. 5 Snapshot of wave-field at time T=2.0 s on the coarse grid(h=50 m)，respectively generated by (a)The eighth-order ONAD;(b)The fourth-order LWC;(c)The eighth-order LWC.

接下来，我们研究八阶ONAD方法模拟声波在更复杂的三层非均匀介质模型(见图 6)中的传播能力.模型计算区域大小、层间界面和各层速度如图 6所示.取空间步长和时间步长分别为h=Δx=Δz=30 m、Δt=1×10^－3 s.震源位于模型区域中心，频散为f₀=20 Hz.

图 6

Fig. 6

图 6 三层介质模型Fig. 6 Three-layer model

图 7a~c分别由八阶ONAD方法计算得到在T=0.5 s、 T=1.0 s和T=1.5 s时刻的波场快照.图 7a中波的传播还未穿过区域Ⅲ和区域Ⅱ的水平和垂直间断界面，即波在区域Ⅲ的均匀各向同性介质中传播，因此，波场快照显示的是以震源为中心的同心圆，波形非常清晰，没有可见数值频散.图 7b中波的传播已经穿过区域Ⅲ和区域Ⅱ的界面但还没有穿过区域Ⅰ和区域Ⅱ的交界面，图 7c中波的传播已经穿过穿过区域Ⅰ和区域Ⅱ的交界面但还没有出模型区域.从这两个图中可以清晰地看到波穿过水平和垂直间断界面后的反射波和透射波，并且在间断面附近没有可见的数值模拟，这说明八阶ONAD方法对于复杂介质中的地震波模拟也是非常有效的.

图 7

Fig. 7

图 7 八阶ONAD方法模拟得到的声波在复杂介质中不同时刻的波场快照 (a)T=0.5 s；(b)T=1.0 s；(c)T=1.5 s.Fig. 7 Snapshot of wave-field generated by the eighth-order ONAD method in complex media at time (a)T=0.5 s，(b)T=1.0 s，(c)T=1.5 s，respectively.

最后一个实验，我们利用八阶ONAD方法来对Marmousi模型(见图 8)进行波场数值模拟.设模型计算区域为0≤x≤9.192 km和0≤z≤2.904 km，取空间步长和时间步长分别为Δx=Δz=24 m、Δt=0.8×10^－3 s，则Marmousi 模型的网格数为384×122.震源位于(4.596 km，0.024 km)，取频率f₀=15 Hz.边界条件处理采用Yang等(2002)提出的二次吸收边界条件.

图 8

Fig. 8

图 8 armousi模型(声波速度变化范围从1.5 km/s到5.5 km/s)Fig. 8 The Marmousi model(The velocity varies from 1.5 km/s to 5.5 km/s)

图 9给出了八阶ONAD方法计算得到从时间T=0 s至1.5 s的地表地震记录图，可见地震记录图非常清晰，无可见数值频散.这进一步验证了八阶ONAD方法在复杂非均匀介质波场模拟中能有效地压制数值频散.同时，这个实验表明八阶ONAD方法可以有效地结合Yang等(2002)提出的二次吸收边界条件开展复杂地震波场模拟.

图 9

Fig. 9

图 9 由八阶ONAD方法计算得到的Marmousi模型地表合成记录Fig. 9 Synthetic seismogram on the surface，generated by the eighth-order ONAD method

有限差分方法是目前使用最广泛的数值模拟方法之一，但是该方法在利用差分替代微分来实现对空间高阶偏导数的离散，从而导致地震波场模拟过程中出现伪波动即所谓的数值频散，进而严重影响了地震波数值模拟精度.为有效地压制地震波数值模拟中的数值频散，基于二维声波方程，在Yang等人的基础上发展了一种高精度的数值模拟格式——八阶ONAD方法.该方法使用四阶截断的泰勒展开式对时间偏导数进行离散化，利用八阶精度的近似解析离散算子对空间高阶偏导数进行离散化.数值误差分析表明，八阶ONAD方法的相对误差明显小于四阶LWC方法和八阶LWC方法.计算效率调查揭示，八阶ONAD方法几乎没有数值频散各向异性，而四阶LWC方法和八阶LWC方法的数值频散各向异性比较明显.另外，计算效率比较显示，八阶ONAD方法的计算速度为四阶LWC方法的4.95倍、为八阶LWC方法的3.53倍；其存储量是四阶LWC方法的19.32%、是八阶LWC方法的35.86%.最后，我们将八阶ONAD方法应用到层间速度反差较大的双层介质模型、三层非均匀复杂介质模型和经典的Marmousi模型中进行声波场数值模拟，结果表明，该方法能有效地压制数值频散，且对于复杂介质中的地震波模拟也是非常有效的.总之，八阶ONAD方法是一种有望在地震波模拟得到应用的十分有效的正演数值模拟技术. 附录A：空间高阶偏导数的逼近公式

利用局部插值近似方法，可以得到位移和位移梯度关于空间的二阶和三阶偏导数的逼近公式，具体推导过程见参考文献[Tong et al.，2013].这里，只给出了关于位移的高阶偏导数的计算公式：

把(A1)~(A6)式中的u换成v就得到关于位移梯度的二阶和三阶偏导数的逼近公式.