波动方程有限差分法频散衰减方法研究

2012-08-20成都理工大学地球物院四川成都610059

石油天然气学报 2012年12期

刘洋，唐湘蓉（成都理工大学油气藏地质及开发工程国家重点实验室成都理工大学地球物院，四川成都610059）

张玺华（成都理工大学沉积地质研究院，四川成都610059）

近年来，相对于其他方法，有限差分法显示了其数值方法的简洁性和生命力。这里通过修改波动方程式来压制有限差分中的频散，新的解法是利用原始波动方程式与一个衰减因子相乘。因此，数值频散在大波数情况下成指数方式衰减，使得有限差分法的效果像高阶一样，但是计算成本要小得多。对模型和真实数据的试验表明了该方法的有效性。由于有限差分法的简单性和可扩展性，有限差分法已经广泛运用于正演和偏移中波场的数值外推，但是必须特别注意的是一种有限差分法固有的误差，即数值频散。为此，笔者提出了一种简单和低计算成本的改进的波动方程式用来压制数值频散。

1 频散的理论分析

在声学中，三维波动方程式如下：

式中，v＝v（x，y，z）为空间速度，m/s；x、y、z为空间坐标；u为位移，m；t为时间，s。

在有限差分法中，数值频散是一个不希望出现的现象。先讨论一维上行波动方程式：

对空间Z方向进行离散，并且用1阶中心差分代替微分：

式中，Δz为Z方向的空间采样间隔，m。

以u＝exp [i （kz＋ωt） ]（其中，kz为Z方向的波数；ω为角频率，rad/s）代入式（3）得：

空间离散引起的相速度与群速度之比为：

式中，O（）为高阶无穷小。

由式（5）可见，空间离散使得相速度低于群速度，并且波数越大的分量滞后得越多。而且当时间和空间采样间隔相对于所使用的有限差分方法过于粗糙，数值频散就会发生。从物理意义上看，这种现象的产生是由于不同的频率具有不同的相速度产生的。而体现在以式（1）为基础的有限差分法正演模拟上，数值频散主要是由下面的原因造成的：①网格太粗糙；②差分算法不精确。

图1（a）显示了一个由式（1）推导出的有限差分方法正演结果，图1（a）使用的是时间2阶精度差分和空间4阶精度差分。虽然在图1（a）中可以看到比较清晰的圆形波场，但是频散几乎覆盖了环形波场内部的全部区域，其主要原因是地震波的算法采用的是低阶（2阶）精度差分方法和太粗糙的差分计算网格。

图1 在各向同性介质中250ms处的波场快照

克服上述问题的一个通常的做法是使用较长的差分算子（高阶有限差分法），另一个常用的方法是减小计算网格的尺寸。由于要达到一定的精确度，每波长上网格点的最小数须足以表示一个波形，算子越长，需要的网格点数越少，反之亦然。但必须注意，每波长的网格点数不能少于2个。

图1显示的是各种不同空间精度有限差分法合成的波场快照，模型是各向同性介质，速度是3500m/s，时间步长是0.001s，X方向采样间隔和Z方向采样间隔均是10m，有限差分计算网格是201×201，震源函数使用的是雷克子波。在图1（b）与图1（a）相比，图1（b）较图1（a）的频散现象已经明显改善。当更进一步，使用8阶精度有限差分法计算同样的模型，其他参数保持不变，结果如图1（c），频散现象较图1（a）和图1（b）进一步改善。图1（d）使用了2阶精度有限差分算法，但是计算所使用的网格大小是图1（a）～（c）所使用的1/4（即网格点数图1（d）是图1（a）～（c）的4倍），其他参数保持不变。综合比较以上4种效果图，图1（b）～（d）展示了比图1（a）好得多的效果。

在三维波动方程式正演中，每一时间计算间隔，减小网格会使计算量呈3次方增长；另外，增加有限差分算子的长度，即用更高阶精度的有限差分算法，将使计算量线性增长。对比这两个方法，后者具有明显的优势。但是，有限差分算子越长有限差分式的稳定性越差。再讨论有限差分法的稳定条件，根据柯朗－弗里德里希斯－列维条件，在三维正演模拟中，每个时间步长必须满足下式：

式中，Δt为时间采样间隔，s；vmax为最大模型速度，m/s；Δx、Δy分别为X、Y 方向的空间采样间隔，m。

由于有式（6）的影响，在二维波动方程式正演下，当网格变小，有限差分总计算量将会成3次方增长；在二维波动方程式正演下，有限差分总计算量将会成4次方增长。图1是使用标准波动方程式在各向同性介质中250ms处的波场快照，震源使用的均是雷克子波，其中图1（a）、（b）、（c）使用的计算网格是201×201。图1（d）使用的是401×401，但是图1（d）在X方向和Z方向的采样间隔是图1（a）～（c）的1/2，网格大小是图1（a）～（c）所使用的1/4。因此，尽管图1（d）比图1（a）的效果要好得多，但是图1（d）需要的计算量大约是图1（a）的8倍。如果能够在2阶精度的有限差分法的基础上压制频散，那么就能够使用比较低的计算量而得到更好效果的结果。

2 声波波动方程式优化

笔者将式（1）的解用下式表示：

式中，r（x，y，z）为随空间变化的非负实数。表面上看，式（9）比式（1）要复杂得多；然而，式（9）等号右边多出的部分只是对式（1）简单的修改。式（9）能够有效地利用低阶（2阶或者4阶）有限差分法。

图2验证了不同空间精度有限差分法计算的单炮正演模型，图2（a）与图2（b）均使用的是空间2阶有限差分法，区别在于图2（a）使用的是式（1），而图2（b）使用的是式（9），图2（b）中的频散比图2（a）中的频散现象明显少得多。图2（c）和图2（d）均使用式（1）计算，但是前者使用的是空间8阶精度有限差分，而后者使用的是空间2阶精度有限差分，且后者使用的网格在X、Z方向的采样间隔均是前者的1/2。图2（b）的效果几乎可以达到图2（c）、图2（d）的效果，在单炮正演记录上基本观测不到频散的影响；而图2（b）的计算量比图2（d）的计算量要小得多。

为了保证式（9）的解有意义，r（x，y，z）必须满足下式：

还有一个简单的处理r（x，y，z）值的方法是把r（x，y，z）设置为一个常数R，R的值与差分方法的阶数和差分法计算所用的网格大小有关，在该次正演模拟中，R值为0.5。

图3给出了修改后的波动方程式空间2阶精度有限差分（图2（b））的正演频谱图（图3（a））和标准波动方程式空间2阶精度有限差分（图2（d））的正演频谱图（图3（b））。比较这2幅图，能够证明修改后的波动方程式能量在频率的分布上几乎没有变化。

式中，X ＝（x，y，z）表示空间位置；p为式（1）的一个解；K ＝（kx，ky，kz）表示空间波数；kx、ky、kz分别为X、Y、Z方向的波数。

上面已经证明了在用有限差分法模拟波场时，大波数的相速度与小波速的相速度不同，而且波数越大的，相比滞后群速度得越多。因此，给出下式：

式中，g（K）为阻尼因子。相比式（7），式（8）多出一个阻尼因子g（K），其是平面波解的阻尼函数，g（K）将在大波数时压制频散噪声。下面给出具有式（1）形式的解：

图2 对由3层不同介质（速度从上到下分别是2700、3100、3500m/s）组成的模型制作的单炮正演记录

3 结语

在波动方程式有限差分解法中，数值频散是固有的一种非物理数值噪声。数值频散可以简单地通过大幅度的计算量来压制。笔者从分析空间差分引起的有限差分法数值频散出发，提出了一个修改的波动方程式，在这个新式中添加了一个阻尼因子，使得新式的解能够在大波数时高速衰减，从而实现低阶精度有限差分法和较粗糙的网格计算得到高阶精度有限差分算法和精细网格才能得到的效果。