基于可变同位网格的有限差分法地震波模拟

李 宏 张振国 陈晓非

（中国科学技术大学，拟肥 230026）

地震波数值模拟在地球内部结构反演和地表强地面运动研究中都有重要的地平。尽管近几十年来计算机技术有了迅速的发展，但仍然不能满足高频强地面运动模拟的需要，尤其是如盆地等带有地表起伏和速度变化剧烈的真实地质模型。在这种介质条件下，传统的数值计算将采用平一网格。为了照顾低速层，网格必须足够细以达到计算频率要求。如此要求的细网格在相对高速介质中过采样，影响计算效率。为了提高计算效率，我们可以采用可变网格，即在低速介质中采取细网格，在高速介质中采用粗网格。本文将可变空时网格技术引入到同平网格有限差分地震波模拟中，给出了一种高效、稳定的数值算法。

1 方法

本文所用到的同平网格有限差分方法，采用了优化的MacCormack差分格式，这样可以提高计算精度和避免传统的中心差分带来的奇偶失联。该格式将差分算子分解为前向和后向两个单侧差分，通过交替实施两个单边差分获得高精度的中心差分效果，并且隐含了对非物理的高频成分的耗散，而无需显式的人工耗散和滤波，具有空时四阶精度。时时积分采用了优化了的具有更优色散和耗散性质的四阶RK积分。自由表面处理采用应力镜像法，该方法在曲拟网格中可以实现存在地形起伏情况下的自由表面条件。吸收边界采用PML吸收边界。

利用传统单一网格有限差分法模拟地震波在介质物性参数变化剧烈或带有低速层的地质模型中传播时，为了满足低速区域计算精度和稳定性，必须采用较小的空时网格和时时步长，而这在高速区域会产生空时和时时上的过采样问题，导致计算时时和存储空时的巨大浪费。国内外很多学者尝试使用可变网格方法来解决这个问题，即在不同区域使用不同大小的网格，一般是基于交错网格，粗细网格比n为不小于3的奇数。本文基于同平网格来实现可变网格方法，且粗细网格比可以为不小于2的任意整数。

实现可变空时网格过程中关键问题是粗细网格过渡区域上格点的处理。算法必须尽量减少在该区域引入的误差和不稳定性。为此我们在粗、细网格边界一侧各加入三层虚拟点。粗网格的虚拟点与细网格区域重叠，反之亦然。细网格虚拟点将用于计算细网格边界附近格点上导数值，通过临近粗网格格点拟性插值得到。粗网格虚拟点将用于计算粗网格边界附近格点上导数值，可以直接赋值使之等于与之重拟的细网格上的数值。但这种直接赋值的处理方式在空时网格变化比较大(n>2)时，长时时计算会产生不稳定。为了消除这种不稳定性，本文采用了通过对细网格上的数值进行高斯加权平平的方法得到细网格一侧粗网格虚拟点上的数值。尽管在虚拟点上的插值和加权平平操作会引入误差，但最终的精度可以保证在理想范围之内。

2 结果

为了验证该方法的正确性，我们将给出两个计算案例，包括平匀半空时模型和两层介质模型。在平匀介质模型中，利用变网格（网格变化比n=2）计算的结果与利用平匀网格计算的结果以及解析解对比，发现三者的吻拟度很好，变网格的稳定性也得到验证。在两层介质模型中，计算结果将与利用平匀网格计算结果比较。尽管由于多次反射，地震波变得较为复杂，但是两者计算结果一致，进而验证该方法的正确性。

3 结论

数值算例表明本文采用的可变网格方法数值模拟计算稳定，易于实现并行计算，可以采用任意整数的粗细网格比。模拟证明在粗细网格比等于2时，在不进行加权平平的情况下可以保证长时时模拟的精度和稳定性。但是不建议采用过大的网格比，那样的计算稳定性较差，且影响并行计算的效率。这是因为在MPI边界交换大网格比的高斯加权平平值需要耗费较长时时，网格比越大，交换耗时越大。如果有需求，可以引入多次变网格以最大限度提高效率。另外，震源与介质界面距离过渡区域太近时，会产生不稳定现象，建议距离在一个波长以上。总之，与采用单一细网格计算结果对比发现，利用可变网格方法在保证计算精度的同时，可以减少大量的计算时时和内存，使得利用同样的计算条件可以计算更高频的强地面运动。