波场模拟笔记.doc

1、关于介质的分类：22、波场数值模拟的分类23、数值频散24、人工边界条件25、地震波计算理论基础35.1运动平衡微分方程(位移与应力的关系)35.2几何方程(应变与位移的关系)35.3本构方程 广义胡可定律(应力与应变的关系)45.4弹性各向同性介质45.5横向各向同性介质(VTI)65.6声波方程及其有限差分75.6.1均匀各向同性介质二维声波方程75.6.2震源函数85.6.3非均匀介质中声波方程交错网格高阶有限差分数值解86、兰姆问题模拟97、一阶速度应力弹性波动方程108、有限差分的基本原理131、有限差分法简介132、有限差分的差分格式139、TI介质的Thomsen参数表征149.1Thomsen参数149.2弹性介质的Thomsen参数表征1410、交错网格及差分格式1510.1时间上的2M 阶差分近似格式1610.2空间上2N阶差分近似1710.3交错网格任意偶数阶精度有限差分系数计算公式1810.4差分格式2011、PML吸收边界条件2311.1PML吸收衰减函数2511.2Cerjan衰减边界条件2511.3左莹用的PML2612、水平自由表面边界条件2713、起伏地表自由边界条件281、关于介质的分类：TI介质：(Transversely Isotropy)其速度(沿层面)在垂直与介质对称轴的平面内保持不变，在纵向上为非均匀性。水平对称轴的横向各向同性介质：(Transverse Isotropy with a Horizontal axis of symmetry )HTI 介质TI介质的对称轴为水平垂直对称轴的横向各向同性介质：(Transverse Isotropy with a Vertical axis of symmetry ) VTI 介质TI介质的对称轴为垂向OA介质：VTI和HTI介质结合在一起形成正交各向异性(Orthorhombic Anisotropy) 2、波场数值模拟的分类从大的方面讲，地震正演模拟可以分为物理模拟和数值模拟两大类。其中地震数值模拟方法可归纳为三类：波动方程法、积分法和射线追踪法。波动方程法是建立在以弹性或粘弹性理论和牛顿力学为基础的双曲型偏微分方程-波动方程的理论基础上，对波动方程进行数值求解可以得到地震波场记录。在求解波动方程时，首先要将模型离散成有限个网格点，然后给出全波场的值，故又称为网格法和全波方程法。当网格足够密时，模拟可达到足够的精度，并适用于介质参数空间变化复杂的模型。包括有限差分法（FD）、虚谱法（PS）和有限元法（FE）。积分法是建立在波动方程积分解表达式的基础之上的，波场由点震源激发，其理论基础是惠更斯原理。其数学表达形式为波动方程的格林函数域积分方程表达式和边界积分方程表达式。射线追踪法是建立在以射线理论为基础的波动方程高频近似理论基础上的，其数学表达形式为程函方程和传输方程。3、数值频散数值频散：有限差分数值模拟是将微分方程离散为差分方程进行求解。当控件网格间距或时间网格大小不适合就会引起波形畸变，地震波在传播过程中，波前形状就会发生变化，并且逐渐离散。解决办法：增加一个波长内网格点的个数。利用FCT( flux-corrected transport ,通量校正传输 )技术压制数值频散。当波场变化剧烈甚至不连续时常规有限差分算法会造成严重的频散误差，利用FCT技术能够获得满足精度的解。4、人工边界条件在波场数值计算中，除了地表自由边界外，有限计算区域的其它边界都是人工截断边界。实际地球介质是横向无限纵向半无限的区域，而波场数值模拟总是在一个有限的空间区域进行，这样在有限的空间区域就引入了介质边界，从而产生一些边界反射，严重干扰了有效波的信息，引起计算的局部不稳定，直接影响数值模拟结果的精度。吸收边界分类：透射边界条件。衰减边界条件。透射边界条件：在人工边界上给出一种透射波场的近似表达式，从而使向外传播的波场“透过”边界而不被反射回去。衰减边界条件：在边界部位设一衰减(粘滞)带，使向外传播的波及其反射波在此带中同时被衰减，以至于没有明显的反射波回传到计算区域中。完全匹配层(PML)可能够较好地吸收入射到边界上的地震波。5、地震波计算理论基础弹性波是机械振动在弹性介质中传播的过程，而地震波是机械振动在岩石中的传播过程。在一定的条件下，当温度、压力和作用力等大小适当，岩石具有弹性性质，可视为弹性体，地震波可视为岩石中传播的弹性波。5.1运动平衡微分方程(位移与应力的关系)由弹性波理论知，弹性体的运动微分方程可表示为 (5.1)式(5.1)中，为时间分量，分别表示弹性体中质点在三个方向上的位移分量，它们是有关及的函数，为弹性体密度，为正应力，为切应力，为外力分量()。5.2几何方程(应变与位移的关系)在弹性力学中，表示弹性体在外力作用下体内各质点的应变与位移沿坐标轴分量之间关系的几何方程为： (5.2)5.3 本构方程广义胡可定律(应力与应变的关系)广义胡可定律规定，在固体中，任意一点的六个应力中的每一个应力都是六个应变分量的线性函数。广义胡可定律可表示为： (5.3)式中弹性系数 = 1,2,3,4,5,6 共36个，由改点三维弹性介质的性质所决定，表示与弹性体有关的弹性常数。由于弹性能量是应变的单值函数，所以弹性系数必须等于。这样，要完全描述一个复杂的弹性介质所需要的弹性系数总共减少为21个。也就是说在各向异性介质中，最多共有21个独立的弹性系数。5.4弹性各向同性介质在弹性各向异性介质中，3.7式变为： (5.4.1)在弹性各向同性介质中，独立的弹性系数只有两个，所有非零的弹性系数都可用拉梅系数和表示，即，其余弹性系数都为零，相应的(5.4)式变为 (5.4.2)将(5.4.2)式和(5.2)式代入(5.1)式中得到弹性各向同性介质中的位移方程，即 (5.4.3)不考虑分量，相应的二维弹性波方程为 (5.4.4)在拉梅系数和不随空间点的变化时，方程(5.6)可以简化为： (5.4.5)令，为纵波速度，为横波速度，则方程(5.7)可表示为： (5.4.6)5.5 横向各向同性介质(VTI)在各向异性介质中，最少共有21个独立的弹性系数，若对其求解非常困难。油气勘探中所遇到的各向异性介质主要是横向各向同性介质(VTI)。其弹性系数矩阵为： (5.5.1)其中, 在这种情况下，有五个独立的弹性常数，。在这种介质中，介质的弹性特征沿横向是各向同性的，沿纵向是变化的，它对于，面都是对称的，则有，又因为它在方向完全相同，对于VTI介质，有，及，则只需要5个独立的弹性系数。于是(5.5.1)式可变为： (5.5.2)将(5.5.2)和(5.2)式带入(5.1)式就可以得到在VTI介质中以位移为变量的波动方程，即 (5.5.3)不考虑变量，则(5.5.3)式相应的VTI介质中的位移方程为： (5.5.4)当弹性系数不随空间点的变化而变化时，(5.5.4)式可以进一步简化为： (5.5.5) (5.5.6)上述方程组分成两组：方程(5.5.6)是横向各向同性介质中SH波波动方程，方程(5.5.6)则是准P波和SV波波动方程，它们是耦合在一起的，即在横向各向同性介质中，P波传播可以引起SV波，SV波传播也可以引起P波，它们联合成为横向各向同性介质中准P-SV波动方程，如果给定的各向异性弹性参数张量满足弹性各向同性参数张量的条件，则方程(5.5.5)和(5.5.6)即简化为均匀各向同性介质中的P-SV波和SH波波动方程。5.6 声波方程及其有限差分5.6.1 均匀各向同性介质二维声波方程均匀各向同性介质二维声波方程可表示为： (5.6.1)其中，为地表记录的压力波场或位移波场，为声波速度场。为导出方程的离散差分格式，需要将对其进行地震波模拟的模型区域离散化，令，则时间导数采用二阶中心差分，空间导数为十阶差分精度的二维声波方程的高阶差分格式为： (5.6.2)其中，5.6.2 震源函数考虑到零相位子波具有高分辨率及无相位畸变的特征，模拟过程中采用Ricker子波函数作为激发震源。雷克子波的数学表达式为： (5.6.3)当进行正演模拟时，须在网格节点上对应的炮点位置加入震源影响。难么(5.5.3)式变为： (5.6.4)5.6.3 非均匀介质中声波方程交错网格高阶有限差分数值解非均匀各向同性介质中二维声波方程的一阶应力-速度方程形式可表示为： (5.6.5)其中，是质点速度，是法线应力，是密度，是纵波速度。采用交错网格进行剖分，相应的声波波场分量和弹性参数在空间的位置见下图。下面给出2N阶空间差分精度、二阶时间差分精度交错网格高阶有限差分格式，设，分别是应力，速度，的离散值。则方程(5.6.5)的差分格式为： (5.6.6) (5.6.7) (5.6.8)二维各向同性介质声波方程交错网格示意图其中，表示，表示，表示,，表示方向的网格间距；表示时间步长。6、兰姆问题模拟初始模型示于图6.1，兰姆问题分析见图6.2，合成记录如图6.3和图6.4所示。如图6.1所示，一法向沿轴方向均匀分布，于时，突然开始作用，这种情况下空间、时的波场传播规律问题称为兰姆问题。杜世通在其所著的地震波动力学理论与方法中指出，当地下介质为半无限均匀介质时，波场中将出现纵波和横波，此外还有一种特殊的波(Schmidt).Schmidt波只在一定范围内(即，为临界角，分别为横波和纵波速度)存在。图6.2给出了兰姆问题理论波场分析示意图，图中箭头表示波的推进方向，P波、SV波和Schmidt波分别用P、SV、S标于图中。由于兰姆问题有明确的解析解，所以人们常用兰姆问题来检验弹性波动方程数值解法的正确性和精确度。图6.3是兰姆问题的合成记录。比较图6.2和图6.3可以看出：正确地模拟出了纵波和横波；正确地模拟出了非均匀的Schmidt波，而且该波出现的空间区域与理论分析结果一致。图6.1模型示意图图6.2 兰姆问题波场分析示意图图6.3模型波场快照图6.4模型单炮记录7、一阶速度应力弹性波动方程对弹性体传播的运动平衡微分方程(5.1),为了避免对介质弹性常数进行空间微分，则可以对以位移表达的波动方程进行处理，即如果用分别表示弹性体中质点在三个方向上的速度分量，在没有外力作用下或停止外力作用后，就可以得到由速度和应力表示的一阶弹性波方程。 (7.1)将公式(5.2)和带入方程(5.3),并对时间求导可得 (7.2)应用方程(7.1)和(7.2)就可以表示任意各向异性介质中的一阶速度应力方程。将方程组(5.2)和带入(5.5.1),并对时间求导就可以得到如下一阶微分方程组。 (7.3)联合方程组(7.2)和(7.3)，便可以得到横向各向同性介质中用质点速度和应力表示的一阶弹性波方程。二维情况下的一阶弹性波方程为 (7.4)利用(7.4)式就可以对二维横向各向同性介质进行弹性波场数值模拟。同理，可得，弹性各向同性介质的一阶速度应力波动方程。 (7.5)相应的二维各向同性介质一阶速度应力波动方程为 (7.6)8、有限差分的基本原理1、有限差分法简介有限差分法（Finite Differential Method-FDM）是基于差分原理的一种数值计算方法。其基本思想是：将场域离散为许多小网格，应用差分原理，用差商来近似的代替微商将求解连续函数的问题转化为求解网格节点上的差分方程组的问题。2、有限差分的差分格式一个函数在点上的一阶和二阶微商，可以近似的用它临近的两点上的函数值的差分来表示。如对一个单变量函数，为定义在区间的连续变量。以步长将离散化，我们得到一系列节点：，然后求出在这些节点上的近似值。显然步长越小，近似解的精度越高。与节点相邻的节点有，因此在点可以构造如下形式的差值：,节点的一阶向前差分,节点的一阶向后差分,节点的一阶中心差分与相邻两点的泰勒展开式可以写为：公式2.1公式2.2利用2.1式和2.2式，并忽略的平方和更高阶的项得到一阶微分的中心差分表示： 公式2.3利用2.1和2.2式我们还可以得到一阶微分的向前、向后差分表示：公式2.4公式2.5将2.1式和2.2式相加，忽略的立方及更高阶得到二阶微分的中心差分：公式2.6利用2.12.6式，我们就可以构造 出微分方程的差分格式。9、TI介质的Thomsen参数表征9.1Thomsen参数弹性介质模型的性质是由刚度矩阵确定的，刚度矩阵确定了应力和应变之间的关系，但由其确定弹性波动方程系数的物理意义不明确，也很复杂。为了方便理论研究和实际应用，围绕波传播的相速度公式，展开公式的物理意义，Thomsen提出了一套表征TI介质弹性性质的参数，定义如下： (9.1)以上定义TI介质的Thomsen参数为：、和，共五个参数。这里的、分别为qP波和qS波沿TI介质对称轴的相速度；、和是表示TI介质各向异性强度的三个无量纲因子。其中，是度量qP波各向异性强度的参数，越大，介质的纵波各向异性越大，纵波无各向异性；是连接和之间的一种过渡性参数；可以看成是度量qS波各向异性强度或横波分裂强度的参数，越大，介质的横波各向异性越大，时，横波无各向异性。一般情况下，和的单调性是一致ide，即同时增减或为零。9.2 弹性介质的Thomsen参数表征根据TI介质的Thomsen参数表征，使得介质弹性参数的物理意义更加明显。Thomsen参数(9.1)式使用弹性参数表示的；反之，TI介质弹性参数也可用Thomsen参数表征。VTI介质刚度矩阵元素的Thomsen参数表征如下： (9.2)10、交错网格及差分格式在交错网格技术中，不仅要求空间网格的交错，而且在时间上交错，其主要有点是无须对介质的弹性常数进行空间求导，因此它常被用于复杂介质中地震波场数值模拟。再者由于一阶速度应力波动方程中仅含一阶偏导数，容易提高差分近似的阶数，达到提高数值模拟精度的目的，所以常采用交错网格有限差分法求解此方程。无论是各向同性还是在横向各向同性介质中，弹性波场的交错网格有限差分法数值模拟的基本过程是：第一步，将介质中相应的地震波传播方程表达成一阶速度应力方程，此时，方程中仅含有应力、速度变量关于时间或空间变量的一阶偏导数。第二步、将应力、速度变量交错的分布在模型的网格点上。第三步、将方程中的所有一阶偏导数用二阶或高阶差分近似，得到相应的差分方程。第四步、将时间变量离散，随着时间的步进，根据差分方程循环交替地计算各网格点上的应力或速度。下面具体说明交错网格有限差分法地震波场模拟的实现方法。用交错网格有限差分法求解一阶速度应力波动方程3.10式时，应力、速度等分量在交错网格中的位置设置见图 10.1。图 10.1 速度、应力等变量在模型网格点上的交错分布示意图10.1时间上的2M 阶差分近似格式在使用交错网格法求解一阶速度应力方程时，速度和应力两个分量分别是在和时刻进行计算的，为了提高时间的差分精度，得到2M阶精度的时间差分近似格式，可以将和用Taylor公式在处展开为：将两式相减可得公式3.11其中，为时间步长。当M等于1是，公式3.11 就是传统的二阶精度差分近似。直接计算3.11式中的要涉及过多的时间层，需要较大的内存量。为此，利用方程3.10式可以完全准确的将速度对时间的任意奇数阶高阶导数转嫁到应力对空间的导数上去，将应力对时间的任意奇数阶高阶导数转嫁到速度对空间的导数上去。这样，在计算一个时间层上速度（应力）场时，只需要前一时间层的速度（应力）场，以及之间的应力（速度）场，不需要过多的时间层，从而节省内存。如2M = 4时，公式3.11式变为：公式3.12公式中的第三部分是通过应力和速度之间的关系变换，使得应力（速度）对时间的导数转化成速度（应力）对空间的导数。10.2空间上2N阶差分近似在交错网格技术中，变量的导数是在相应的变量网格节点之间的半程上计算的。为此我们用下式计算方程上面方程中的一阶空间导数。上式中待定系数的准确求取是确保一阶空间导数的2N阶差分精度的关键。将和在处用Taylor公式展开后可以发现，通过下列方程组即可以求得待定系数。10.3交错网格任意偶数阶精度有限差分系数计算公式设有7阶导数，则在处的7阶泰勒展开式分别为于是，可见，上式仅为奇数阶导数项，偶数阶导数项被消除了。同理，可得，的泰勒展开式。由于一阶导数 6阶精度中心差分近似式可表示为 其中，为常数，于是，有+=化简可得系数矩阵，即求解系数方程可得，由以上过程可计算任意2L阶精度中心有限差分系数。首先来推导交错网格一阶导数的6阶精度差分系数。交错网格上一阶导数的6阶精度差分系数计算由离散点，确定。要确定网格点上的一阶导数，取，于是有在的一阶导数由离散点，确定。根据泰勒展开式，有此处需要完善求解系数方程组可得，于是如算子的对称点为，则上式可改写为同理，可以推导出交错网格任意2L阶精度有限差分系数计算公式。设有2L+1阶导数，则，在处的2L+1阶泰勒展开式为由于交错网格一阶导数2L阶精度差分近似式可表示为将上式L个方程带入、化简，有其中，差分系数由以下方程确定下面给出了几种不同空间差分精度的差分系数：当2N=4时；=1.125; =-0.04166667;当2N=10时；=1.211243; =-0.08972168; =0.009570313; =-0.00176566; =0.0001186795;10.4差分格式我们采用精度较高的交错差分网格，对一阶弹性波方程4阶时间差分精度高阶近似方程3.12，给出2N阶空间差分精度的差分格式。差分网格如图2.1所示。（图2.1只是显示了应力和速度在空间上的差分，并没有显示速度和应力在时间上的差分，速度和应力的计算是分别在和时刻进行计算的）设,和分别是速度、与应力、的离散值。为了方便起见，假定，则3.12式的精度为的差分格式为：速度水平分量其中，,速度垂直分量其中，应力Pxx其中，应力Qzz其中，应力pxz其中，以上是一阶应力速度弹性波动方程的空间2N阶时间4阶交错网格数值模拟的差分格式。计算步骤：1、 输入主体网格与交错网格的介质本构弹性参数2、 确定激发源的空间位置3、 计算时刻的速度分量；4、 计算时刻网格上应力和；5、 速度、应力分量的时间递推，交错运算。通过以上步骤，就可以求得整个网格节点上不同时刻的波场快照，获得共炮点记录。11、PML吸收边界条件PML的概念有Berenger在1994年首先提出的（A perfectly match layer for the absorption of electromagnetic waves 1994）。它将弹性波分量在吸收边界区域分裂，并能分别对各个分裂的波场分量（P,S,R波等）赋以不同的耗损。这样，就能在时域有限差分法（FDTD）仿真计算区域截断边界处得到一种非物理的特殊吸收介质，该层介质的波阻抗与相邻介质的波阻抗完全匹配，因而入射波将无反射地穿过分界面进入PML层，并且，由于PML层为有耗介质，进入PML层的入射波将迅速衰减。即使PML为有限宽度，它对于入射波仍有很好的吸收效果。文献指出，PML的反射系数是其它吸收边界条件的1/3000，使用完全匹配层可以使FDTD模拟的最大动态范围达到80db。但是PML方法的编程较复杂，需要存储的信息量也比较大。下图是该方法吸收边界处理示意图.图11.1 吸收边界处理(PML)示意现在定义弹性波入射情况下的PML介质层。定义的关键是将弹性波波速和应力分量分裂成5个子分量，记为和，当确定计算区域进行波场模拟时，PML介质中5个分量又分别分为2个，于是共10个波场分量，它们分别为：，、和、。这些波场分量满足以下方程(弹性各向同性介质)： (11.1) (11.2) (11.3) (11.4) (11.5)公式(11.1)(11.5)中，为吸收函数。用交错网格将匹配层控制方程离散化，匹配层内传播的波是逐渐衰减的，只在内部计算区域与匹配层交界的地方由于模型的离散化而产生很小的边界反射，并随着离散网格宽度的减小而趋于理想的边界反射系数，对内部计算区域影响不大，故可采用低阶格式提高计算效率，则匹配层内的差分格式为：11.1PML吸收衰减函数在完全匹配层(PML)的应用中，衰减函数的选择是关键，如果选择不当将会产生严重的边界反射。Collino等提出的一种以内域截断边界与PML边界距离为指数关系的衰减函数： (11.6)式中为PML吸收层内节点到PML吸收层内边界的距离，分别为方向与方向的衰减因子，为最大纵波速度，为匹配层宽度，为理想的边界反射系数(一般取)。当不等于零时表示衰减，当等于零时表示不衰减。11.2Cerjan衰减边界条件Cerjan等基于各向同性的假设提出了衰减边界条件，其基本原理是在计算区域周围引入衰减带，使向外传播的波及其反射波在此衰减带中同时被衰减，一致无法明显反射到计算区域。在每个时间迭代步的计算中，衰减带内的速度和应力的振幅都会通过乘以下面的衰减因子得到衰减。原文衰减因子为： (11.7)衰减因子使波场值在扩充区域内逐渐衰减，其中衰减系数，取值范围为(01)，是衰减带宽度内的网格点数。但是衰减系数和参数要凭经验选取，如果吸收系数选择不合适，则仍然会有很强的边界反射。针对这一缺点，后人提出了以下两种改进：余弦衰减函数：，用以为周期的余弦函数代替原文中的指数衰减函数，可以消除由于递减剧烈产生的反射。抛物线衰减函数:,对抛物线平方项乘以加权系数，增加了选择的自由度，并且可以使波场值在衰减带内平滑递减。为加权因子，理论取值范围一般为(01),为衰减带的总网格点数，为衰减带内待计算的网格点到衰减带内边界的网格点数，简称衰减距离。余弦衰减函数的周期为，衰减带横跨个周期，同常情况下需要满足，结合衰减量及余弦函数的周期性可以容易的确定的范围。在抛物线衰减函数中，为衰减函数的最小值，可以根据衰减量及衰减带厚度来选取加权因子。显然，改进的离散型吸收边界条件是在吸收区域内使波的振幅迅速衰减，从而不产生边界反射，其优点是适应性强、实施简单，适应各种不同的介质，不需要增加附加条件。该方法可以应用到各向异性介质数值模拟中，其扩充吸收带容易加大计算区域的外面，而且反射波衰减得很好。11.3左莹用的PML此方法的基本思想是在研究区域的边界上加入吸收层，使边界上传入吸收层的波随传播距离按指数规律衰减，不产生任何反射，从而达到消除边界反射的目的，下面介绍一下完全匹配曾吸收边界在高阶交错网格中的应用。从地震波的质点速度应力方程可以看出，二维地震波场的质点速度、应力分量都可以分解为两个方向上的次分量：水平方向和垂直方向。按照这个原则我们可以将U分解为：其中，可分别写成式中的和分别为x方向和z方向的衰减系数，我们选择其为其中，为匹配曾区域与内部区域界面的横向距离，为匹配曾区域与内部区域界面的纵向距离，为最大纵波速度，L为匹配层宽度，系数。式子经过变换可以写成：将上式的计算结果代入迭代计算后的U，经过这样的计算边界就可以被吸收掉了。同理，可以得到W、P、Q、S在PML边界处的公式。12、水平自由表面边界条件在水平自由边界上应力分量必须满足这些条件确定了有限差分算法中的质点振动速度的迭代方法。根据上式，在水平自由边界上有，因此在边界上x方向的正应力随时间的变化可表示为利用镜像法可以使应力满足自由的边界条件，得到相对精确而稳定的数值计算方法。本文采用让自由界面通过切应力的节点，而不是和的节点，由于正应力不在自由表面上，因此不需要考虑上式，其在镜像公式：和的作用下自动满足。其次，当边界的形状形成一个顶角时，这种镜像法只需要剪应力为零，因此提高了效率，物理上也是正确的。在水平自由表面上，竖直速度分量W在计算时牵涉到密度，若和分别表示自由界面上和自由界面以下一个节点处的介质密度，这时，则自由界面存在一个反射界面，从而使得球面波在界面反射产生瑞利波。13、起伏地表自由边界条件在上一章中给出了水平自由表面通过切应力的算法，通过分析可以将这种算法推广到不规则起伏地表的情况。为了将这种方法应用于起伏地表的情况，必须仔细的分析水平表面、竖直表面以及他们形成的顶角情况。不规则起伏地表可以分解为许多水平或竖直的表面单元和它们之间的顶角，Robertsson 仔细分析了各种组合情况，把水平地表的镜像法推广到不规则起伏地表的情况。本文结合robertsson的方法以及只允许自由界面通过网格单元切应力点的算法结合，提出了一种新的起伏地表的实现方法。由于根据Robertsson 的算法将位于起伏地表附近的网格点分为七类需要复杂的判断，本文采取新的实现方法，在新方法的实现中，我们并不对这些网格点作分类，而是对起伏地表的拐点进行分类。下图表示不规则起伏表面分解成许多小的表面单元以及其拐点，黑色实现代表起伏地表，可以看出在我们的新方法下自由界面只通过切应力点。这些拐点可以分为5种不同的情况。如下图所示：凹陷地表拐点分类图凸起地表拐点分布图以上两图的正方向分别是水平向右和竖直向下起伏地表自由边界条件在自由表面上，地震波传播满足应力为零的边界条件，在二维起伏地表模型中，我们可以将自由表面分为水平自由表面和竖直自由表面。在水平自由表面上应力分量为零：(3)根据(3)式，：(4)因此，在水平自由边界上，沿方向的正应力随时间的变化可表示为：(5)在竖直自由表面上应力分量为零：(6)根据(6)式，在竖直自由表面上有：(7)因此在竖直自由边界上沿方向的正应力随时间的变化可以表示为：(8)起伏地表处理算法分析常见起伏地表可分为凸起地表和凹陷地表两类，地表拐点及网格分布如图2，将不规则起伏表面分解成许多水平界面和竖直界面以及二者之间的拐点，黑色粗实线代表起伏地表，本文将拐点共分为5种不同的情况且自由表面只通过切应力的采样点。 将位于水平边界上的单元网格节点称为边界点，通过将起伏自由表面水平边界上的某边界点高程和其前后边界点高程的对比可以将拐点分为五类。对于边界的上升或者下降的定义原则是：从正演模型的左侧向右行进，如果界面上升则定义为上升，如果界面下降则定义为下降，如临近的三个界面点的高程相等则定义为水平。假设分别代表边界点处的高程。这五类拐点定义的说明以及镜像处理时的具体规则如下：H：广义水平拐点类包含以下三种情况 U：上升拐点类 D: 下降拐点类 UH: 上升至水平拐点类 HD: 下降至水平拐点类(b