声波有限差分数值模拟课件

声波有限差分数值模拟

——基于完全匹配层吸收边界条件以下讨论声波模拟的有限差分方法地震资料的数值模拟作用：

1）模拟地震记录，检验处理结果的好坏、处理方法的有效性；

2）正演模拟可以作为反演研究的基础。固体弹性介质简化为声学介质：研究地震波传播问题及地震成像方法时，为了方便求解，只研究纵波的波场特征及成像方法。这种做法是对实际问题的良好近似。因为地表附近存在低速带，地震反射记录中横波信息的能量非常微弱。内容提要一、声波方程二、声波方程的差分近似三、吸收边界条件

一、声波方程声波方程：通常，介质密度相对于其速度变化很小，可以近似地将密度看作常数：即其中：为声压为密度为速度为震源为梯度算子二、声波方程的差分近似三维常密度声波方程：其中：为频带有限的震源子波，震源位于处。对微分方程进行差分近似时，主要利用了Taylor级数展开理论：（2-1）二、声波方程的差分近似设是一个多元函数，是它的一个自变量的差分近似时，先将记为为了表述方便，当研究（2-2）（2-3）由（2-2）和（2-3）得（2-4）二、声波方程的差分近似（2-4）中分别取并做误差截断得：令（2-5）二、声波方程的差分近似（2-5）可以写成（2-6）求解方程组（2-6）得（2-7）所以的六阶误差截断表达式：用类似的方法可以求得的四阶误差截断表达式：（2-8）的二阶误差截断表达式：（2-9）（2-10）二、声波方程的差分近似二、声波方程的差分近似取x、y、z、t方向的差分网格间距分别为、、、令，在点考虑差分近似以时间二阶空间四阶误差截断为例，有：将上面式子代入声波方程（2-1）就可以得到对应的差分方程：，震源项取为并将速度取为二、声波方程的差分近似（2-11）二维情况下，上式变为：（2-12）其中，整理（2-12）式，得在边界处需要降阶二、声波方程的差分近似三、吸收边界条件（图3-1）（图3-2）在用有限差分法进行声波数值模拟时，若采用自然边界条件，即令各个时间层边界处的声压值为零，在波场没有到达人为边界时计算结果是正确的（图3-1）；而当波场超出人为边界时，计算结果中混入了人为边界的反射（图3-2）。（图3-2）（图3-3）工程中的许多波动问题，其计算求解域往往很大，有时甚至是无界的。但实际情况下，我们总是在一个有限的区域内进行求解，因此，需要在所限定的区域的边界上引入吸收边界条件，从而最大限度地降低由于人为划定的边界而造成边界反射。吸收边界的建立对于波动方程的数值模拟起到致关重要的作用，它将直接影响的计算结果的稳定性和准确性。图3-2没有加吸收边界；图3-3加入了吸收边界。三、吸收边界条件

一些用于消除人工边界反射的方法包括：

Smith于1974年提出多次反射法；

Clayton及Engquist于1977年提出的旁轴近似法；

近年来，Berenger针对电磁波传播情况，给出了一种高效的完全匹配层（PML）吸收边界条件，并在理论上证明该方法可以完全吸收来自各个方向、各种频率的电磁波，而不产生任何反射。目前已将完全匹配层的思想引入到了弹性介质中，并建立了弹性介质中的完全匹配层吸收边界条件。以下利用声波方程，导出声波的完全匹配层的控制方程，并作为吸收边界条件应用于数值模拟。实际结果表明这是一种非常有效的方法。声波方程：其中（3-1）等价运动方程连续性方程（3-2）（3-3）质点速度等价（3-6）（3-5）（3-4）三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（3-9）（3-8）（3-7）完全匹配层边界条件：当时，上面方程就是声波方程；可以证明，当时，在方向上起到的衰减作用，以下给出证明。三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

将（3-7）、（3-8）、（3-9）变换到频域（3-12）（3-11）（3-10）对x、y、z做下面坐标拉伸变换：其中同理（3-13）（3-14）三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（3-13）、（3-14）代入（3-10）、（3-11）得（3-12）（3-11）（3-10）（3-10）、（3-11）、（3-12）变回到时域：（3-15）（3-14）（3-13）（3-9）（3-8）（3-7）（3-15）（3-14）（3-13）（3-6）（3-5）（3-4）两组式子形式一样坐标变换因此，若声波方程的解为PML方程的解应该为三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

有了上面的推导作为基础，下面以二维情况为例，说明PML的衰减作用。设二维声波方程的一个平面简谐波解为：若此简谐波在PML介质中传播，则对应的声压解应该为：坐标变换两个解的相位完全相同，且振幅比为：（3-16）三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（3-16）（3-16）可以解释PMLABC效用和其参数的效用：当时，，振幅不发生衰减。当时，可以使，振幅发生衰减，且和分别控制着各自方向上的衰减。时，振幅无衰减，随着它们的增大，振幅呈指数衰减。以下给出（3-7）（3-8）（3-9）交错网格离散表达式：（3-7）以为例，在进行离散（3-8）(3-17)(3-18)对（3-17）和（3-18）做以下变量替换，替换时需要注意在（3-17）和（3-18）中，为了得到从到的声压波场的递推公式，(3-17)(3-19)(3-18)(3-20)(3-9)(3-19)(3-20)交错网格时间层递推以下以二维情况为例说明PML吸收边界的实现方法：三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（图3-4）如图3-4，中间的部分为进行有限差分模拟的区域，为了加入PML吸收边界条件，需要在它的周围加入PML介质，图中，L代表左边界，R代表右边界，U代表上边界，D代表下边界，RD代表右下角边界，LD代表左下角边界，RU代表右上角边界，LU代表左上角边界。垂直方向设置为有衰减，这时可以让从边界处为零值开始向外逐渐增大到某个值。水平方向设置为有衰减，这时可以让从边界处为零值开始向外逐渐增大到某个值。在D和U处，水平方向设置为无衰减，即在L和R处，垂直方向设置为无衰减，即在RD、LD、RU、LU处，水平和垂直方向均设置为有衰减。

图3-4中的细实线，代表每个数组的次边界，粗实线代表每个数组的边界，数值模拟时，每完成一次时间层递推，一个数组里的边界值要和相临数组的次边界值发生交换。三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

如图3-5所示，下面的两个“回”字矩形代表 时刻时间层上的数据，一个数组里的边界值要和相临数组的次边界值发生交换，上面两个矩形代表由 时刻的时间层数据递推得到的 时刻的时间层数据，完成计算后，再进行边界和次边界的数据交换，生成新的数组边界数据，依此类推完成整个时间层的计算。（图3-5）三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

如图3-6，为了程序实现时便于操作数组，将D和U合并在一起存放的一个数组里，将L和R合并在一起存放的一个数组里，将RD、LD、RU、LU合并在一起存放的一个数组里。这种做法的另一好处就是，它将图1-3-4中周围的PML介质四周也连接起来，这样不需要在每次计算完一个时间层时对图1-2-4的最外边界赋零值，而是让波场继续传播，多通过一倍的PML介质，衰减更充分，却没有增加任何计算量。（图3-6）（图3-4）三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

下面两个图说明了PML介质中速度的一种选取方法。以下给出一些使用了PMLABC的二维声波数值模拟的算例时间切片，震源坐标为震源处使用的是一个峰值频率为，峰值振幅为1，的Ricker子波。

介质速度为延时为图3-7、图3-8、图3-9、图3-10，给出了四张时间分别为三、吸收边界条件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（图3-7）t=0.1s（图3-8）t=0.25s（图3-9）t=0.04s（图3-10）t=0.55s三、吸收边界条件