做了大概三个月的波场模拟，本来只是本科老师一门课的课设，却激发了我无限的兴趣。通过这次地震波场的模拟，深深的感觉理论与实践结合是多么的重要。我这里对想从事波长数值模拟的人又找不到比较系统入门的方法做个总结，希望能给那些对数值模拟方面感兴趣的人有所帮助。
    众所周知，地震波的传播规律可以用一双曲线方程表示，我们称之为波动方程，我们需要在已知地下介质参数和物理模型的情况下，得到各个质点的振幅信息。通常在地震勘探中有两种形式表示，二阶的位移-应力方程和一阶速度-应力方程；介于一阶速度-应力方程不需要对弹性参数进行空间微分，并且一阶速度-应力方程的耦合性，利用现在已经非常完善的交错网格技术来进行模拟，对模拟的精度和计算速度，交错网格比常规网格有着更大的优势（当然有比这更好的网格离散型式），所以我决定采用交错网格对声学介质中的声波方程进行离散。
     1、第一步，当然是对要进行求解的波动方程进行离散，在有限差分的思想中，对方程中的微分项我们采用差分来近似代替。而交错网格中核心思想是，1）将所求的变量分别储存在不同的网格系统中，2）交错网格将变量的导数是在相应的变量的网格的半程上计算的。本次一阶速度-声压方程有三个变量，速度的水平分量，速度的垂直分量，声压。我将各个变量的储存位置由如下图所示：


     2、利用差分的思想，对各个方程进行离散化，离散后的交错网格任意阶数方程表达式可以参考文献：《波动方程数值模拟的三种方法及对比》。任意阶交错网格的系数可以参考文献：董良国《一阶弹性波方程交错网格高阶差分解法》。那么最后我们得到了所需要的离散方程组。，对于交错网格在文献：《P-S wave  propagation  in  heterogeneous  media: Velocity-stress  finite-difference  method》有详细的介绍。这里我给出自己对交错网格中自己的对空间导数半个网格点的理解，以便对编制程序有所帮助。
     3、利用交错网格对方程进行离散后，我们通常会得到半个网格点的节点编号，由于各个变量储存在不同的网格系统中，所以对节点的编号的协调性就显得尤为重要。这里我采用在空间中对网格进行平移的思想来进行协调，例如，水平速度分量储存在网格的（I,J+1/2)位置，而声压储存在网格节点上（I,J),那么我们对水平速度分量所在的网格系统整个向左平移半个网格点，那么这个点的速度分量就可以和声压所在的网格点重合了，那么就不存在半个网格点，对垂直速度也是一样的理解，从而在程序的节点编号就可以实现了。记住一点：在循环中，每一次循环，求出的变量值是同一时间同一节点而不同位置的变量值。
     4、对于利用有限差分对波动方程的数值模拟，会产生两个比较大的问题：数值频散和人工边界反射。我们需要单独处理。
1）利用差分代替原微分方程，那么必然会有截断误差，这样造成方程的弹性参数系数变大或变小，造成波在传播的过程中的相速度和群速度不一致，造成伪波动，那么在波场中就会造成波场扩散，会出现很多波前面，这种现象被称为数值频散。对于这个方面得详细讨论，可参考文献：董良国<地震波传播数值模拟中的频散问题>。
2）我们进行波场模拟只能在有限的区域中进行模拟，那么必然带来的人工边界的问题。在计算区域之内，各个物理参数和介质参数不为0，而在计算区域外，介质参数都为0，这样在边界上必然会产生反射，而在现实中波的传播是在无限的介质中，没有这样的反射产生，如果不加以处理，会扭曲波在无限介质中传播的真实性。
       关于这两个问题的解决方法有大量的文献可以参考，以及各种方法的优势与弊端都有很好的说明，我看了大量的文献，选取其中两种效果比较好的方法进行简单的说明。
    5、对于频散方面，造成的因素主要有三个：1）地震波的传播方向；2）空间差分的精度；3）一个波长内离散点数目；如对弹性波动方程，（则还有泊松比的影响）
     对于差分精度，我采用高阶差分与交错网格结合，提高空间差分的精度来降低频散；另外利用FCT（通量传输校正），减小采样间隔和空间步长来增加一个步长中网格点的数目。
对于交错网格高阶差分的精度与稳定性分析可以参考文献：董良国《一阶弹性波方程交错网格高阶差分解法稳定性研究》。对于FCT校正的详细原理可以参考文献：《Elimination of numerical dispersion infinite-difference modeling and migration by flux-corrected transport》。
    6、对于边界反射方面，很多人在这方面研究，并且已经产生了很多消除边界反射的方法，总体上可分为局部人工边界方法和全局人工边界。通过前人大量的研究，PML方法被目前认为是最有效的边界反射处理方法。这里我对PML方法做一简单的说明。PML方法在计算边界处引入一个吸收层，在这个层中，将所求的未知量的各个分量分解为两个辅助分量，一个平行于边界面，在方程中只保留平行于交界面方向的偏导数分量；一个垂直于边界面，只保留垂直于交界面方向的偏导数分量。最后对垂直于边界的分量引入阻尼因子，使界面上不产生反射的特性。由于在界面上，阻尼因子为0，与计算区域的方程相同，故不会产生反射，达到了消除边界反射的目的。在加入PML后，可以不必单独对方程在计算区域和PML区域内单独处理，只需要在计算区域内令阻尼因子为0，那么在全区域都可以由一组方程组来进行求解。实行起来非常方便。关于对PML的详细介绍，可以参考文献：《Application of the perfectly matched absorbing layer model to the linear elastodynamic problemin anisotropic heterogeneous media》
    7、通常我们还需要边界条件（不同于人工边界），通常解决数学模型问题，我们需要偏微分方程+定解条件来描述这样一个数学模型，这样构成一个定解问题，这样方程才会得到特解。定解条件分为初始条件+边界条件。本文采用如下定解条件：
初始条件：T<=0时刻介质处于平衡状态，此时介质中的速度和应力处处为0.
边界条件：1）内部弹性界面：这种边界是由于介质中速度和密度等弹性性质发生变化而形成的两层介质的分界面，在差分计算中，通过公式中弹性参数的改变自然就能体现出来，不需要单独处理；
        2）自由边界：交错网格中有限差分在自由表面可以采取0应力的边界条件。
这里对自由边界，给出两点说明：1）：自由表面上方区域的波场值遵从反对称原则；
                              2）：自由表面上方区域的速度值设置为空气层，即V=0。
     8、下面我对整个模拟过程进行图片说明：
首先是频散的表现：未加FCT

  进行FCT后：

对于边界反射：

加入PML后：注释：顶界面满足自由边界条件，所以会有正常的反射。

以上都是基于自己的理解，难免有些理解错误的地方，望读者能够指点出来，我会加以改正，我会继续往下深化！

希望能给想给做这方面初学者起到一定作用，那我就倍感荣幸了！通过这次小研究，对学习研究有所感悟，通常最重要的就是需要阅读大量的文献资料，这样才对这个方面有比较全面的了解，做研究一定不能急于求成，往往是万事开头难！