The Compensation and Correction of Numerical Dispersion of Seismic Wave Propagation
Li Junfeng；Li Guoming；Li Rushan
（College of Geophysics，Chengdu University of Technology，Chengdu 610059，China）
摘要：波动方程有限差分法正演模拟对地震属性、资料地质解释、地震波传播规律、储层评价等研究均具有重要意义。波动理论正演方法是对建立的模型进行网格化的方法，根据波动方程近似，按时间步长迭代计算波场，获得正演波场模拟结果，但有限差分技术求解波动方程时会带来不必要的震动，即数值频散。文章探讨了引起数值频散的原因，并以高阶交错网格差分和通量校正传输（FCT）相结合的方法为基础消除数值频散。
Abstract: The finite difference offorward modeling of wave equationis of great significance for seismic attributes, the data of geological interpretation, seismic wave propagation and reservoir evaluation studies. The establishment of the model is divided to grid, based on wave equation approximation,calculation by time step iterativeof wave field, access to the results of the forward wave field simulation. However, solving the wave equation finite-difference technique will cause unnecessary vibration, which is numerical dispersion. In this paper, I discusses the reasons caused by the numerical dispersion, and based on a combination of methods of order difference equations and flux corrected transport(FCT) to eliminate numerical dispersion.
关键词：波动方程有限差分交错网格数值频散
Key words: finite-difference；staggered grid；wave equation；numerical dispersion
中图分类号：P315.3+1 文献标识码：A文章编号：1006-4311（2011）19-0295-02
0引言
地震数值模拟是在假定地下介质结构模型和相应物参数已知的情况下，模拟研究地震波在地下各种介质中的传播规律，并计算在地面或地下各观测点所应观测到的数值地震记录的一种地震模拟方法。波在传播过程中，波前形状发生变化，并且逐渐散开。这是因为波传播相速度和群速度不一起所引起的，这种现象称为数值频散现象（王才经，1990）[1]；董良国（2000）提出高阶差分方法（包括交错网格）可以显著地降低数值频散[2]；吴国忱等（2005）针对TI介质波场模拟的交错网格有限差分方法，从空间网格离散、时间网格离散和算子近似等三个方面对其产生的数值频散进行了分析[3]。牟永光，裴正林（2003）对三维复杂介质地震数值模拟进行研究。为消除波场模拟中的数值频散问题，诸多学者从不同的角度对有限差分方程的数值频散进行了分析，并给出了相应的解决办法。因此本文通过交错网格技术与高阶差分方程技术有机结合的方法进行分析。
1有限差分模拟的研究进展
有限差分法是在非均匀介质中模拟声波和弹性波传播最主要的数值方法。有限差分方法实质是将波动方程中的偏微分算子用差分算子代替，因此差分格式决定了差分算法的稳定性，网格的划分决定了差分算法的精度。Alterman（1968）首先将有限差分法应用于层状介质弹性波的模拟中。Claerbout（1970）完成了有限差分法在勘探地震学中最著名的应用和改进。Virieux（1984）发展的交错网格方法使用一阶的速度应力方程，其差分精度为O（Δt2+Δx2），在不增加计算工作量和存储空间的前提下，和常规网格相比局部精度提高了4倍。Levander（1988）这种差分网格的精度提高到O （Δt2+Δx4）。Dablain（1986）提出了求解声波方程的高阶差分方法。Crase（1990）将这一方法运用到求解二阶弹性波方程中，提高了计算效率。董良国（2000）采用将速度（应力）对时间的奇数阶高阶导数转化为应力（速度）对空间的导数，将交错网格和高阶差分法有机结合，求解横向同性（Tl）介质一阶速度-应力弹性波方程[2]。
2有限差分模拟数值
有限差分方法（FDM）是计算机数值模拟最早采用的方法。该方法将求解域划分为差分网，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程。该方法是一种直接将微分问题变为代数问题的近似数值解法。
2.1 构造差分的方法构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。不同的微分方程所对应的差分格式是不同的。一般而言，是先做推导，得到差分格式之后再动手写程序。
2.2 交错网格有限差分法 从精度和收敛性上来考察交错网格有限差分方法。对于交错网格差分和规则网格差分，把两者对函u（x）=eiωx一阶导数近似的误差进行比较，当进行二阶近似时，交错网格差分误差和规则网格差分误差的比为0.25，当四阶近似时，误差系数比为0.141，交错网格差分计算精度明显高于规则网格差分。对于复杂构造、复杂地质体、复杂岩性地震波场计算，交错网格高阶有限差分法综合其在占用内存、计算精度、计算效率和并行算法实现等方面的特点，可构成一种综合性能最好的波场计算方法。
3数值频散压制方法研究
提高差分方程的精度包括两个方面：空间差分精度和时间差分精度，单独提高空间差分或时间差分的精度都不会很好的压制网格数值频散。只要空间网格间距和时间网格间距均满足采样定理，不产生时间和空间假频，随着差分方程阶数提高，对物性参数改变就越小，进而数值频散越少。所以，可以通过提高差分方程阶数的方法来提高波场模拟精度，压制数值频散，但随着方程阶数的增加，必然会增大计算量。因此要采用交错网格技术和高阶差分方程有机相结合，在不增加内存量和较少计算量就可提高波场及正演纪录的模拟精度。稳定性问题，是高阶有限差分数值求解波动方程中需要仔细考虑和研究的基本问题。由于差分计算中数值参数选择不合理，很可能造成模拟结果网格频散严重，影响对模型的分析，严重时会引起数据溢出。因此需要进行参数设置对结果的稳定性影响的分析，给出解法的稳定性条件。地震波场模拟中的数值频散是不可避免的，但当我们掌握了它产生的原因和性质后，可以通过从各方面进行压制。压制频散的方法有以下三种：
3.1 差分算子影响波动方程系数的变化，我们可以通过校正差分算子来降低频散。可在频率－波数域进行差分算子补偿，但有一定的难度实现在时间域中进行补偿递推显格式的差分算法。
3.2 提高时间域和空间域差分精度。通过提高差分精度可有效地压制频散，但随着方程阶数的增加，必然增加计算量和计算机内存，为此要采用交错网格技术和高阶差分方程有机结合。我们可以通过选取空间域差分阶数N的值来提高空间差分阶数。通过Taylor展开我们可以得到时间域的2M阶精度差分近似，通过选取M的值来提高时间差分阶数是一个很有效的办法。但应用到波动方程时会遇到对空间求二阶以上的导数运算的问题，推导过程繁琐复杂，容易出错，为了弥补这一缺点，用Runge-Kutta时间离散格式提高时间差分精度。采用Runge-Kutta时间离散格式几乎压制了所有快纵波的频散，效果非常明显，唯一不足的地方是慢纵波的频散还能看到，但考虑到慢纵波本身所具有的频散特性，以及慢纵波的高衰减性，这一现象是可以接受的。
3.3 通量校正传输（FCT）方法是首先应用在流体力学连续守恒式方程的求解中，发展的通量校正传输方法，而后将FCT法应用于求解声波方程，有效压制了在粗网格情况下差分计算产生的数值频散；守恒格式保证差分方程与物理定律相互协调，它可以很容易的处理交错网格的差分频散校正[4]。在交错网格有限差分法中，我们要求在网格边界上有物理量流入和流出，根据能量守恒关系定律，保证每个网格点在同一时刻的通量相等，如果不等就认为是有数值频散存在，因此当两个网格的差分方程相加时，流过共同边界上的通量可相互消掉，所得的差分方程就和能量守恒关系直接应用于组合网格得到的差分方程相同。FCT方法是基于守恒型方程的差分格式，在构造这种差分格式时，要求网格间交界面处守恒物理量的通量互相抵消，即交界面两侧一边流入和另一边流出的通量相等，这样内部单元边界上的截断误差被抵消，避免误差的积累。用FCT方法进行校正，广义地说主要有三步：交错网格差分计算、通量漫射校正和通量抗漫射校正。在波场模拟中，数值频散是不可避免的，即在网格间波场物理量违反了守恒定律，从而引起波场通量漫射，将FCT方法引入到弹性波正演模拟中来，可消除通量漫射以减小数值频散。FCT方法可以用较低阶的差分方程达到较高精度的正演模拟，大大提高了运算效率，但这种方法需要更多的内存量。由于假设所有的极值点都是由数值频散引起的，传统的FCT方法对所有网格点进行通量校正处理,再对局部极值点进行补偿的逆通量校正优化的FCT，通量校正只用在局部极大值上，节省了大约40%的计算量. 同时FCT方法可以放大时间和空间步长，从而抵消计算FCT带来的计算量的增加。通过各种方法的比较我们可以看出，FCT是很好的一种消除数值频散方法。由于漫射因子的取值需要进行大量试验才能得到，这是该方法的主要缺陷。
4结论
交错网格有限差分技术具有效率高、精度大、便于实现的优点，现在被广泛应用于声波方程、黏弹性介质波动方程、双相介质波动方程、各向异性波动方程等介质模型中。空间和时间的网格离散都会带来数值频散，对时间离散格式提高时间差分精度，取得了较好的效果，而通量校正传输法（FCT），对波场通量漫射进行校正，也可有效地压制数值频散。用通量校正（FCT）与交错网格差分法作为正演模拟的理论基础，不仅保证了计算波场的精度，强压制了数值频散，而且大大提高了计算速度。
参考文献：
[1]王才经．波动方程模拟和偏移的频散分析[J]．计算地球物理研究，1990（33）：521-528．
[2]董良国，马在田，曹景忠等．一阶弹性波方程交错网格高阶差分解法. 地球物理学报，2000，43（3）：411-419．
[3]吴国忱，王华忠．波场模拟中的数值频散分析与校正策略．地球物理学进展，2005，20（1）：58-65．
[4]Tong Fei，Ken Larner．Elimination of numerical dispersion in finite-difference modeling and migration by flux-correctedtransport[J]．GEOPHYSICS，1995，60（6）：1830-4842．