二维频率域声波方程正演模拟

2D Acoustic Wave Equation Forward Modeling in the Frequency Domain

摘要: 频率域正演在地震波数值模拟中占有十分重要的地位。相比于时间域正演，频率域正演具有适合多炮并行运算、无时间频散、频带选取灵活、误差小等优点。频率域正演时，不同频率的系数矩阵相对独立，适合并行计算加速处理，大大提高计算效率。这里针对频率域声波方程的最优化9点差分格式，开展了隐式表达方式和稀疏矩阵的求解研究，正演模拟了地震波场。模型计算验证了方法技术的准确性和有效性。

Abstract: Forward modeling in frequency domain plays an important role in the numerical simulation of seismic waves. Compared with time domain forward modeling, frequency domain forward modeling has many advantages, such as suitable multi shot parallel operation, no time dispersion, flexible frequency band selection and small error. The coefficient matrix of different frequencies is relatively independent in the frequency domain forward modeling, which is suitable for the acceleration of parallel computing and greatly improves the computing efficiency. In this paper, for the optimal 9-point difference scheme of frequency domain acoustic equation, the implicit expression and sparse matrix solution are studied, and the seismic wave field is simulated forward. The accuracy and validity of the method are verified by model calculation.

1. 概述

频率域波场正演相对于时间域数值模拟来说，有其自身的优势。首先，在多炮数值模拟情况下，频率域相对于时间域效率更高，每个频率成分的阻抗矩阵只需要计算一次，加入并行计算后可以极大地提高计算效率；其次，就方程本身而言，频率域的衰减项更容易添加；另外，频率域模拟方法是计算区域所有网格点上的全部时间上的频率切片解，误差将分配到所有参与计算的每一个网格点上，因此不存在类似时间域模拟的累计误差，可以考虑长时间的地震波数值模拟。

Lysmer和Drake用有限元方法实现了地震波的数值模拟，初次提出了频率域正演模拟 [1]。Marfurt指出频率域模拟地震波，能够非常有效地模拟出各向异性介质中波的传播及其衰减效应，而且适合大量炮点的同时模拟 [2]。Pratt和Worthington将频率域正演引入到波形反演中，应用常规二阶差分法(5点有限差分格式)对声波方程进行了离散化处理，结果表明当每个波长取至少13个网格点时，相速度误差就可以控制在1%以内，研究过程中详细推导了频率域声波方程和弹性波方程差分格式，构建了阻抗矩阵，指出频率域正演模拟适合多炮计算(每个频率成分共用一个阻抗矩阵)，便于模拟衰减效应，为频率域波场正演奠定了基础 [3]。由于频散严重和巨大的计算量等问题，频率域波场正演模拟没有得到广泛地应用，Jo和Shin等人在研究频率空间域声波波场正演模拟过程中，首次引入旋转坐标算子，将常规网格需要的5个网格点扩充到9个网格点，实现了频率域粘滞声波的正演模拟，频散得到了很好的压制，每个波长只需求4个网格就可以得到1%的精度要求，这使得频率域正演得到了极大地发展 [4]。Shin在Jo的基础上，由9点差分格式扩充到25点，相对提高了精度，这样在正演过程中单个波长内的网格点可以减少到2.5个，得到了理想的结果 [5]。

进入21世纪后，频率域正演方法上有了更大的提高。Min等人在Shin和Sohn研究的基础上做了进一步的研究，首先提出了一种新的25点有限差分加权最优化差分算子，通过Gauss-Newton法求出最优化加权系数，最后将频率空间域弹性波正演结果和方程的解析解进行了对比，从结果看出，利用新的差分算子计算的结果对于复杂介质来说结果比较理想，随后他又优化了25点算子的频散系数，通过最优化原理确定了新的加权系数，频散效应得到了更大限度的压制，并且减少了单位波长内所需的网格点数 [6]。高精度的时间域正演非常依赖于交错网格法，而频率域方法主要是基于混合网格法，Hustedt等人研究对比了频率域地震波模拟中交错网格法和混合网格法，指出对于二维模拟来说，交错网格法在内存需求和计算效率远不及混合网格法 [7]。Xu等从CPU计算时间、内存需求和硬盘读写速度三个方面研究了频率域声波模拟结果，并对比了LU分解和迭代法的优劣 [8]。

频域地震波数值模拟需要求解一个大型稀疏矩阵方程，该矩阵的带宽取决于正演算法(有限差分方法)，由于该矩阵具有高稀疏化、非对称性以及非正定性的特点，这些增加了求解的难度 [9] [10]。大型稀疏矩阵的求解方法分为直接法和迭代法。直接法中应用较为普遍的是LU分解方法，对于二维模型来说，LU分解能够得到数值解，多炮情况下计算效率高，但是直接法对内存的需求较大。迭代法对内存的需求小，但是并行计算效率不如直接法。这里我们首先简要叙述了地震波频率域模拟方法的发展及研究现状，继而论述了有限差分数值模拟方法的基本原理，研究了频率域声波方程的最优化9点差分格式的隐式表达方式和稀疏矩阵的求解，最后通过正演算例验证了方法技术的准确性和有效性。

2. 频率域声波方程

时间域声波方程为：

$\frac{{\partial }^{2}u\left(x,z,t\right)}{\partial x^2}+\frac{{\partial }^{2}u\left(x,z,t\right)}{\partial z^2}-\frac{1}{v{\left(x,z\right)}^2}\frac{{\partial }^{2}u\left(x,z,t\right)}{\partial t^2}=-f\left(x,z,t\right)$

式中的v是介质的传播相速度，非均匀介质中为空间位置的函数，f为震源函数，正演的时候通常选用雷克子波。将上式变换到频率域，即对t做傅里叶变换， $u\left(x,y,t\right)$ 变为 $u\left(x,y,\omega \right)$ ， $f\left(x,y,t\right)$ 变为 $f\left(x,y,\omega \right)$ 。便有了下面的声波方程：

$\frac{{\partial }^{2}u\left(x,z,\omega \right)}{\partial x^2}+\frac{{\partial }^{2}u\left(x,y,\omega \right)}{\partial z^2}-\frac{{\omega }^2}{v^2}u\left(x,z,\omega \right)=-f\left(x,z,\omega \right)$

3. 波动方程正演模拟

有限差分数值模拟方法的原理是以离散代替连续，以差分代替微分，从而将波动方程差分化，这样就将一个微分问题转化为离散的数学问题 [11] [12] [13]。首先，将研究区域进行网格剖分；其次，采用有限差分算子将连续性方程在每个网格点处离散化；最后，求解差分离散后得到的差分方程组，得到每个网格点处解的近似值。这里我们采用Jo和Shin等人提出的最优化9点差分格式进行计算。根据最优化9点差分格式，可以构造出一个庞大的矩阵方程，该方程用矩阵运算简单表示为： $A\left(v,\omega \right)P\left(\omega \right)=-B\left(\omega \right)$ 。

式中， $B\left(\omega \right)$ 为模拟角频率为 $\omega$ 的雷克子波的傅里叶变换结果，矩阵 $B\left(\omega \right)$ 只在震源处有值，其他处都为零。其中 $A\left(v,\omega \right)$ 是一个与模型参量、有限差分离散格式、频率和边界条件有关的矩阵，称为阻抗矩阵(Impedance Matrix)。阻抗矩阵的大小取决于我们参与计算的网格大小，实际模拟过程中，阻抗矩阵往往是几万阶的矩阵，零元素占据大部分空间，所以我们称 $A\left(v,\omega \right)$ 为大型稀疏矩阵。

求解此方程通常采用以下两种方法。一种是直接法，此方法通常是对系数矩阵做LU分解，进而用分解结果回代求解方程。同频率的多炮正演共用一个分解结果，在二维情况下，计算效率较高。系数矩阵虽然是一个稀疏矩阵，但是分解后的结果会有大量非零注入元，存储分解结果需要很大的存储空间，以致三维的存储变得不现实。另一种方法为间接法，此方法通过迭代来求解方程的解，常用的方法有Jacobi迭代、Gauss-Siedel迭代、松弛迭代等；预条件对加速收敛的意义重大。此方法避免了直接法的大矩阵存储问题，节约内存，但是计算效率偏低，失去了多炮共享分解结果的优势。这里我们基于直接法进行求解。

4. 边界条件

对于频率空间域地震波数值模拟来说，不设边界条件很容易产生明显的边界反射，从而对物理波场产生严重干扰，边界条件的构造非常重要，边界反射的压制效果直接影响着正演的精度。因为频率切片是对全网格上的数值计算，如果边界条件构造的不合理，会产生强烈的边界反射。最佳匹配层(PML)吸收边界条件，最初是在电磁学中提出，其主要思想是在研究区域四周加上由幂函数或余弦函数组成的PML吸收层，使边界上传入吸收层的波随传播距离呈指数衰减，吸收效果相对于阻尼法和旁轴近似法要好。这个方法在地震数值模拟中越来越受到欢迎。这里我们应用此方法对边界进行处理。

5. 模型计算

我们用单层介质模型来做试验，网格大小为200 × 200，网格间距是10 m，震源坐标(1000, 1000)，检波器位于地下，距地面50网格点上，如图1所示。时间采样间隔为10 ms，时间记录长度为4 s，震源为雷克子波，主频为30 Hz，纵波速度为3000 m/s。我们取得频率范围为0~100 Hz，频率间隔为0.25 Hz，共400个频率切片。

Figure 1. Velocity model of homogeneous medium

图1. 均匀介质速度模型

频率切片如图2所示，可以看出，随着频率的增大，波长随之变小，表现在频率切片上就是随着频率的增加，同心圆之间的距离越来越小，单位面积上同心圆个数越来越多。同时我们可以发现30 Hz频率实部切片幅值最大，这是由于我们在震源选择时，雷克子波的主频为30 Hz。由于同频率的虚部切片与实部切片有相似的波场特征，这里不再列出。

Figure 2. Real part of frequency slice ((a)-(c) are 10 Hz, 30 Hz and 50 Hz real part slices)

图2. 频率切片的实部((a)~(c)为频率是10、30、50的切片实部)

图3是反傅里叶变换得到的单炮记录。反傅里叶变换频率选择上，我们分别用400、300、200、100个频率成分进行试验。通过4个不同频率成分做傅里叶反变换的结果发现，随着傅里叶反变换所使用频率成分的增加，单炮记录反射轴越来越连续，分辨率越来越高，这和震源子波的主频有关。

Figure 3. Single shot records in time domain ((a)-(d) are time domain single shot record after 400, 300, 200, 100 frequency components inverse transformation)

图3. 时间域单炮记录((a)~(d)分别是400、300、200、100个频率成分反变换的时间域单炮记录)

图4是时间切片图，图4(a)和图4(b)时刻，震源刚刚开始振动，没有传播到边界，波峰能量较强；图4(c)已经传播到PML内能量衰减严重。我们用400个频率切片做反傅里叶变换，因为频率间隔为0.25 Hz，所以单炮记录上采样时间长度为4 s。所用的最大频率切片为100 Hz，所以单炮记录上的采样间隔为0.01 s。反射同相轴是一条曲线，符合地震波传播的实际情况。从首波到达的时间对比上，震源位于模型区域中心，距离计算边界距离有50个网格点，则它们之间的距离为50 × 1 0= 500 m，正演地震波速度为3000 m/s，则地震波到达中心检波器时间为间为500/3000 = 0.167 s，与实际模型相吻合。

Figure 4. Time slice ((a)-(c) are 0.10 s, 0.20 s and 0.30 s time slice)

图4. 时间切片((a)~(c)为0.1秒、0.2秒和0.3秒时的时间切片)

6. 结论

正演是反演的基础，正演的精度和效率决定了反演的结果。从模拟结果我们可以得出三条结论：

1) 频率切片的实部和虚部从低频到高频变化时，同心圆逐渐加密；

2) 在同频率成分下实部和虚部幅值大小不同；

3) 反傅里叶变换时，频率越多时间域精度越高，但是频率成分超过震源子波有值区域后，即使频率成分增加，精度也不会有太大变化。

基金项目

本文为中国地质调查局科研项目(201100307)资助的成果。

NOTES

^*通讯作者。