二维声波频率域正演

3
二维频率域正演相比时间域正演在效率上有着巨大优势，在保证模拟精度的前提下，频率域方法可以利用LU分解的中间结果方便快速的计算多炮的单频波场，并且它能够非常有效的模拟各项异性介质中波的传播以及其衰减效应，可以方便地模拟粘弹性介质中不同频率波场能量的衰减。
但是频率域正演也有自身的缺点，最重要的缺点就是其内存需求比时间域大得多，尤其是在3D模型中[]。另一个缺陷就是，在高性能计算快速发展的今天，单炮情形下的基于LU分解的频率域正演无法适应MPI的并行算法以及GPU加速，但如果采用易于并行的迭代法进行求解，则正演效率会大大降低。频率域正演若要大规模应用于3D模型的全波形反演，需要解决其内存需求过大的问题以及难以并行的问题。
图2-2九点差分格式下的阻抗矩阵结构
2.3
在地震波场的正演模拟中，网格节点的数目往往达到数十万，因而阻抗矩阵A将会是一个数十万阶的矩阵，常规的消元法是无法求解如此大规模的方程的。在数值计算领域，求解大型矩阵方程的方法往往有两种，一种是迭代法，给定一个初始解，通过迭代更新使其不断逼近真实解。另一种是直接法，其代表方法是LU分解法，将系数矩阵分解为一个上三角矩阵和一个下三角随着目前国内外勘探开发的不断深入，人们对高精度成像，精确地震反演等技术提出了更高的要求。全波形反演是近些年发展起来的一种反演方法，它能够充分利用地震资料中的时间，振幅，相位等信息，具有高精度多参数建模的能力[]。全波形反演既可以在时间域实现，也可以在频率域实现，频率域反演因为其天生的多尺度特性而受到广泛关注。频率域反演的核心问题是高效的频率域正演模拟方法，它决定了反演的精度和效率，因而近年来频率域正演方法发展迅速。
图3-17.32Hz单频波场
图3-212.21Hz单频波场
图3-315.87Hz单频波场
可以看到，随着频率的增加，波峰和波谷交替越来越密集，波长越来越小。在0-36.63Hz内以0.49Hz的频率间隔进行各频率的波场模拟，对得到的单频波场进行傅立叶反变换进而得到时间域的波场。选取其中第406ms的波前快照如图（3-4）。利用时间域正演得到相同时刻的波前快照如图（3-5），两者的波场基本是一致的，但是也存在着可以接受的微小误差。
就数值模拟本身来说，频率域数值模拟相对于时间域模拟也有着许多优点。在频率域计算波场的衰减效应比时间域更方便。而且对于二维多炮数值模拟，频率域模拟方法比时间域模拟方法更有效。另外，频率域模拟方法是基于单个频率对所有空间网格点进行求解，误差将分配到各个网格点。由于各个频率间是独立计算，所以不存在误差的累计，不同频率波场间可以方便地进行并行计算，也可以考虑进行长时间地震波数值模拟。
2
2.1
时间域常密度二维声波方程可写为：
对上式作傅里叶变换得到：
其中 是 的傅立叶变换， 是震源函数 的傅立叶变换， 为角频率，数值上满足 。
可以看到，式（2-2）其实是一个关于波场值 的方程组，它可以写成如下形式：
A称为阻抗矩阵，它依赖于地下介质的参数以及特定频率。U是一个列向量，它是把特定频率的波场值写成了列向量的形式，容易知道，当有限差分的网格规模为 时，U和S都是有 个元素的列向量，阻抗矩阵A是一个 阶的方形矩阵。这时，我们就把求取波场值的问题转化为了一个方程组的求解问题，不过这个方程组阶数巨大，需要用特殊方法求解。
图2-1九点差分格式示意图
混合网格（i，j）处的波场值U对空间的二阶导数可以写为：
而波场值 则使用其周围八个节点的波场值及其自身的值作加权平均：
其中，最优化系数a=0.5461，b=1-a,c=0.6248,d=0.09381,e=（1-c-4d）/4。
依据上面的思路，可以将频率域声波方程进行离散，进而写成式（2-3）的形式。
2.2
正演时如果不对边界进行处理，则模拟得到的波场包含有大量反射波，无法得到正确的结果。在模拟时间域的声波传播时，如果不添加吸收边界条件，在波传到边界之前不产生边界反射，只有在初至波到达边界才产生反射，然后与有效波加。但频率域声波模拟不同，它的每一个单频切片包含所有时间，所以如果不加吸收边界条件，即使在零时刻也得不到正确的解[]。PML边界是目前使用最为广泛的边界条件，且PML边界本身就是在频率域进行推导和证明的，因而频率域PML边界的引入较时间域更为简单。
参考文献
频率域声波方程带有PML边界的形式可以写为：
一般认为式中的一阶偏导项影响很小，故将其忽略，即B=0，D=0。式中 ，
利用上节的方法对式（2-6）在（i，j）处进行离散可以得到：
其中 。
在所有网格点处进行空间离散，并将其写成式（2-3）所示的方程组形式。式（2-7）的系数构成了阻抗矩阵A。可以知道，矩阵A是一个大型的稀疏矩阵，频率域模拟的关键问题就是此稀疏矩阵方程组的求解。矩阵A的结构如图（2-2）所示。
LU分解法使得频率域波场模拟相对时间域有着巨大优势，在观测系统和频率不变时，LU分解的中间结果是可以重复利用的，当完成了第一炮的求解后，其余炮可以调用已经分解好的中间结果，从而在极短的时间内完成求解。
LU分解一般可以调用第三方的数学函数库完成，这一类的函数库有UMFPACK，SUPERLU等。但是这一类函数包都是在Linux/Unix系统下编写的，Windows系统调用起来较为麻烦，可喜的是UMFPACK包可以在MATLAB中调用，而MATLAB本身就有强大的矩阵运算功能以及完备的函数库，因而本文选择了MATLAB环境搭配UMFPACK进行模拟。
3
3.1
设计一个2000m*2000m均匀介质模型，空间网格步长 。纵波速度为2000m/s。时间步长 。震源为雷克子波，主频为10Hz，道间距10m，震源位于模型中央。PML层厚度50层，理论反射系数 。分别对频率为7.32Hz，12.2Hz，15.87Hz的分量进行模拟，得到它们的实部和虚部如下图：
3.2
图3-6corner-edge模型
Corner-edge模型的网格规模为200x200，网格间距10m。介质速度为1500m/s和2500m/s。用雷克子波激发，震源位于模型顶层的中间。接收点位于地表，道间距10m，共200道。通过正演获取12.2Hz的单频波场如下图：
在0-36.63Hz之间以间隔0.12Hz对各频率进行正演，对得到的各单频波场进行傅立叶反变换，得到820ms时的波前快照如图（3-8）。将接收点接收到的各单频波场做傅立叶反变换可以得到时间域的地震记录，如图（3-9）。从时间域的波前快照和地震记录中可以明显看到来自地下的反射波和尖灭点的绕射波，也证明了本文的频率域波场模拟的结果是正确的。
频率域求解的一个重要问题是如何将式（2-2）写成式（2-3）的形式，这其中涉及到差分阶数和差分方法的问题。目前使用较多的有五点差分法，九点差分法[]和二十五点差分法[]，对于本文的声波方程，九点差分法具有足够的精度，下面介绍九点差分法的实现。
常规的十字型差分法只使用离散点上下及左右方向的点，九点差分法用到了离散点周围所有点，包含左上，左下，右上，右下等方位（图2-1）。它其实是将常规正交网格和45度旋转网格结合的产物。