应用配景

有限差分数值仿真技术是声波或弹性波场数值模拟中最为流行的方法之一，然而传统的有 限差分方法在求解颠簸方程时，会产生不期望的数值频散或称网格频散，导致了数值模拟结果 分辨率的降低。究其根源，是因为基于颠簸方程的有限差分求解过程，通常是利用一离散化的 有限差分方程去迫近颠簸方程，从而使得相速度酿成了离散空间间隔的函数。因此，当每一波长内空间采样太少（即空间网格太粗）时，就会产生数值频散。

所谓数值频散实质上是一种因离散化求解颠簸方程而产生的伪颠簸。这种频散既差别于颠簸方程自己引起的频散，也差别于因波流传的速度、频率和角度而引起的频散，它是有限差分 方法求解颠簸方程时所固有的本质特征，无法制止。

为了消除这种数值频散，人们进行了大量的研究[1][2][3]，他们的结论是基本一致的，即 为了消除数值频散，在使用二阶有限差分方法时，每个二分之一功率对应的波长至少必须使用 nn个网格点，而四阶有限差分则可用二阶差分网格点数的一半。王才经 (1990) 也得到了基本 一致的结论[4]。从他们的研究结论中，我们似乎可以得出，提高有限差分格式的精度，可充实 减少每一个波长内的采样点。但事实上并非如此，因为再提高差分格式的精度，将会大大增加 计算量。因此，通过提高有限差分格式的精度以便在粗网格计算条件下消除网格频散的方法并 不可取。

在流体动力学连续方程的求解中，Boris和Book (1973)[5]、Book等 (1975)[6]发展了一种通量校正传输方法 (FCT)，并将他们的FCT方法应用于声波方程的求解中，有效地压制了在粗网格情况下差分计算产生的数值频散。杨顶辉和滕吉文（1997）将这种FCT技术与求解各向异性颠簸方程组的有限差分方法结合，获得了一种适用于求解各向异性介质中二阶声波和弹性波方程的FCT有限差分算法。其过程主要包罗三步：有限差分计算、平滑方程的解和反漫射处置惩罚。

正如上面所述，在使用有限差分求解颠簸方程时，如果每一波长的采样点太少，会产生 数值频散。因此，为了消除差分计算引起的频散，必须作漫射通量的计算，通过漫射通量来 修改差分方程的解 ( 即平滑方程的解 )。由于这种平滑处置惩罚是对计算区域内的每一个网格点实施的，从而一定水平地压制了真实的颠簸，造成振幅精度的部门损失。固然这种精度的损失并 不影响FCT技术的推广，因为损失是不明显的。漫射是以守恒方式实现的，即任何时候在任一 点去掉一部门波场时，就必须有相同数量加回到某一地方。因此，一旦解发生漫射，就必须在 任何不需要漫射的地方引入反漫射以抵消漫射，以使真实波场得以恢复。因此对解作反漫射修正，比漫射计算更复杂，通常是通过一非线性局部搜索来实现。

研究进展

基本理论

通量传输校正技术虽然能有效的压制数值频散，但存在较难选取合适参数的问题，从而影响校正效果。王婷等[10]将接纳改进后的FCT与旋转交错网格相结合，借助旋转交错网格来实现对弹性波方程的有限差分数值模拟，在增加少量计算机内存的情况下有效压制粗网格造成的数值频散。

FCT校正[9]是在假设所有的极值点都是由数值频散引起的前提下，对到场计算的所有网格点进行扩散通量校正处置惩罚，再对非局部极值点进行逆扩散通量校正来赔偿损耗。FCT旋转交错网格有限差分方法的实现分两个阶段：旋转交错网格有限差分阶段和通量校正阶段。通量校正阶段又分为两步：平滑方程的解和反漫射矫正。

其具体步调如下：

（1）用旋转交错网格有限差分格式表现弹性波颠簸方程。求解 n+1n+1 时刻和 nn时刻的速度分量。

首先建立应力-速度方程的旋转交错网格有限差分格式，进行计算得到 n+1n+1 时刻的速度分量 vxn+1(i,j),vzn+1(i,j)v_{x}^{n+1}(i, j),v_{z}^{n+1}(i, j) 和各应力分量 τxxn+1(i,j),τzzn+1(i,j),τxzn+1(i,j)\tau_{x x}^{n+1}(i, j),\tau_{zz}^{n+1}(i, j),\tau_{xz}^{n+1}(i, j) 的值。为了节省时间，提高计算效率，通量校正运算都是只针对速度分量 vx,vzv_{x},v_{z}进行计算，而无须校正应力分量，因为一旦速度分量被校正完毕，对应力分量的校正可通过弹性波方程来完成。

（2）以速度vxv_{x}为例，计算 nn 时刻的弥散通量 PP 和 QQ 。

以速度分量 vxn+1(i,j)v_{x}^{n+1}(i, j) 为例，计算 n−1n-1 时刻的漫射通量

，Pi+1/2,jn−1=η1((vx)i+1,jn−1−(vx)i,jn−1)Qi,j+1/2n−1=η1((vx)i,j+1n−1−(vx)i,jn−1)，(1)\begin{aligned} P_{i+1 / 2, j}^{n-1} &=\eta_{1}\left((v_{x})_{i+1, j}^{n-1}-(v_{x})_{i, j}^{n-1}\right) \\ Q_{i, j+1 / 2}^{n-1} &=\eta_{1}\left((v_{x})_{i, j+1}^{n-1}-(v_{x})_{i, j}^{n-1}\right) \end{aligned}，(1)

其中，参数 η1\eta_{1} 为漫射因子， 0⩽η⩽10 \leqslant \eta \leqslant 1 ，其巨细取决于有限差分阶段频散误差的巨细，可以通过简单数值实验来确定，它可以为一个常数，根据实际计算情况而差别。

（3）运用n+1n+1 时刻计算得到的弥散通量对所有网格点(包罗非局部极值点)上的速度 vxv_{x} 进行修正。

这一步是对波形进行平滑处置惩罚，在压制网格频散的同时，也会造成一定的振幅损失。

，(vx)i,jn~+1=(vx)i,jn+1+(Pi+1/2,jn−1−Pi−1/2,jn−1)+(Qi,j+1/2n−1−Qi,j−1/2n−1)，(2)\begin{aligned}(v_{x})_{i, j}^{\tilde{n}+1} &=(v_{x})_{i, j}^{n+1}+\left(P_{i+1 / 2, j}^{n-1}-P_{i-1 / 2, j}^{n-1}\right) \\ &+\left(Q_{i, j+1 / 2}^{n-1}-Q_{i, j-1 / 2}^{n-1}\right) \end{aligned}，(2)

（4）反弥散通量的计算

1）计算 k+1k+1 时刻的弥散通量 PP 和 QQ ；

，Pi+1/2,jn+1=η2((vx)i+1,jn+1−(vx)i,jn+1)Qi,j+1/2n+1=η2((vx)i,j+1n+1−(vx)i,jn+1)，(3)\begin{aligned} P_{i+1 / 2, j}^{n+1} &=\eta_{2}\left((v_{x})_{i+1, j}^{n+1}-(v_{x})_{i, j}^{n+1}\right) \\ Q_{i, j+1 / 2}^{n+1} &=\eta_{2}\left((v_{x})_{i, j+1}^{n+1}-(v_{x})_{i, j}^{n+1}\right) \end{aligned}，(3)

其中，参数 η2\eta_{2} 的选取与 η1\eta_{1} 类似，但其值要略大于 η1\eta_{1} ，这是因为反漫射运算不但要赔偿认为加入的漫反射的损失，还要赔偿传统有限差分运算带来的振幅损失。

2）计算散射后的 vxv_{x} 差别值，得到反漫射通量：

，Xi+1/2,j=(vx)i+1,jn~+1−(vx)i,jn~+1Zi,j+1/2=(vx)i,j+1n~+1−(vx)i,jn~+1，(4)\begin{aligned} X_{i+1 / 2, j} &=(v_{x})_{i+1, j}^{\tilde{n}+1}-(v_{x})_{i, j}^{\tilde{n}+1} \\ Z_{i, j+1 / 2} &=(v_{x})_{i, j+1}^{\tilde{n}+1}-(v_{x})_{i, j}^{\tilde{n}+1} \end{aligned}，(4)

3）利用反漫射通量修正平滑解 (vx)i+1,jn~+1(v_{x})_{i+1, j}^{\tilde{n}+1} 即可得到消除数值频散后的正确波场值。

，(vx)i,jn+1=(vx)i,jm~+1+(Xi+1/2,jc−Xi−1/2,jc)+(Zi,j+1/2c−Zi,j−1/2c)，(5)(v_{x})_{i, j}^{n+1}=(v_{x})_{i, j}^{\widetilde{m}+1}+\left(X_{i+1 / 2, j}^{c}-X_{i-1 / 2, j}^{c}\right)+\left(Z_{i, j+1 / 2}^{c}-Z_{i, j-1 / 2}^{c}\right)，(5)

其中，

Xi+1/2,jc=sx1⋅max{0.0,min[sx1⋅Xi−1/2,j,abs⁡(Pi+1/2,jn+1),sx1⋅Xi+3/2,j}\begin{aligned} X_{i+1/2, j}^{c} =s_{x}^{1} \cdot \max \left\{0.0,\min \left[s_{x}^{1} \cdot X_{i-1 / 2, j},\right.\right. \left.\operatorname{abs}\left(P_{i+1 / 2, j}^{n+1}\right) , s_{x}^{1} \cdot X_{i+3/2, j} \right\} \end{aligned}

，Zi,j+1/2C=sx2⋅max{0.0,min[sz2⋅Zi,j−1/2，abs⁡(Qi+1/2,jn+1,sz2⋅Zi,j+3/2]}\begin{aligned} Z_{i, j+1 / 2}^{C} &=s_{x}^{2} \cdot \max \left\{0.0, \min \left[s_{z}^{2} \cdot Z_{i, j-1 / 2}，\right.\right.\\ &\left.\operatorname{abs}\left(Q_{i+1 / 2, j}^{n+1}, s_{z}^{2} \cdot Z_{i, j+3 / 2}\right]\right\} \end{aligned}

sx1=sign⁡(Pi+1/2,jn+1)sz2=sign⁡(Qi,j+1/2n+1)\begin{aligned} s_{x}^{1} &=\operatorname{sign}\left(P_{i+1 / 2, j}^{n+1}\right) \\ s_{z}^{2} &=\operatorname{sign}\left(Q_{i, j+1 / 2}^{n+1}\right) \end{aligned}

反弥散通量的计算不是对所有的网格点进行处置惩罚，仅针对进行漫射通量校正前未产生频散的网格点（即非局部极值点），通过修正被平滑过的非局部极值点方程的解，恢复不需要进行漫反射位置（未产生数值频散位置）的真实波场。

步伐实现（MATLAB）

步伐实现（C++）

当空间采样变得很粗时，FCT方法无法完全恢复数值频散造成的分辨率损失（Fei and Larner, 1995; Yang et al., 2002）

参考