技术领域
[0001] 本
发明属于油气勘探领域，涉及一种基于PML吸收边界的三维声
波数值模拟方法。
背景技术
[0002] 基于PML吸收边界的三维声波模拟中，将计算区域分成两部分：一部分是内部计算区域，另一部分是PML吸收区域，如图2所示，其中边界灰色阴影部分表示PML吸收区域，内部空白为内部计算区域，①表示PML吸收区域的
角点边界，②表示PML吸收区域的棱边界。
[0003] 以N×N×N
三维网格，PML边界厚度为M个网格点的模型为例，其模拟方式有两种：
[0004] 方式一：采用规则网格，即所有物理量及压
力值都定义在网格
节点上，如图3所示。
[0005] 在内部计算区域利用二阶微分方程计算压力P的值，需要分配的内存为N3个浮点类型；在PML吸收区域，计算压力P的值每个网格点还需要引入6个辅助变量，利用PML吸收控制方程模拟
地震波在PML吸收区域的传播，消除人工边界的反射，需要分配的内存为N2×M×6×7，因此方式一的内存需要存储N3+42×M×N2个浮点类型的变量。
[0006] 方式二：采用交错网格，即在内部计算区域和PML吸收区域均利用一组一阶的
应力-速度方程在交错网格上模拟，如图4所示。
[0007] 在图4中，①表示压力P值，②、③、④分别表示质点的速度分量vx,vy,vz。
[0008] 虽然在PML吸收区域不需要引入辅助变量，但是每个整数网格点需要为压力值共分配(N+2×M)3个浮点类型变量，在半数网格点(即整数网格点的中点)上还需要为质点的速度分量分配3×(N+2×M)3个浮点类型变量，因此共需要分配4×(N+2×M)3个浮点类型变量。
[0009] 由此可见，上述方式一存在的技术问题是引入的辅助变量过多，PML边界处理较为复杂；
[0010] 而上述方式二虽然利于PML边界的处理，然而需要同时在整个网格定义压力值和质点速度，内存使用量明显增加。
发明内容
[0011] 针对
现有技术中存在的上述技术问题，本发明提出了一种优化的基于PML吸收边界的三维声波数值模拟方法，在简化PML边界处理的同时，有效减少了内存占用量。
[0012] 为了实现上述目的，本发明采用如下技术方案：
[0013] 基于PML吸收边界的三维声波数值模拟方法，包括如下步骤：
[0014] s1、定义混合网格
[0015] 在内部计算区域应用规则网格，仅在规则网格整数节点上定义波场压力；
[0016] 在PML吸收区域应用交错网格，在交错网格整数节点上定义波场压力，在交错网格半数节点上定义质点的速度分量和介质
密度；
[0017] s2、不同区域波场的计算
[0018] 在内部计算区域采用二阶双曲型偏微分
波动方程更新波场压力值；
[0019] 在PML吸收区域采用一阶应力-速度方程计算质点的速度分量和波场压力值。
[0020] 在步骤s2中，首先结合内部计算区域最外层的波场压力值，计算PML吸收区域的质点的速度分量值，然后计算PML吸收区域的波场压力值。
[0021] 在步骤s2中，在PML吸收区域内最靠近内部计算区域的波场压力值，用于更新下一时刻内部计算区域的波场压力值。
[0022] 本发明具有如下优点：
[0023] 本发明针对现有技术中规则网格PML边界处理需要引入大量辅助变量和交错网格内存需求大的问题，提供了一种基于混合网格的模拟方法，即通过在内部计算区域采用规则网格和二阶双曲型偏微分波动方程求解，在PML吸收区域则采用交错网格和一阶应力-速度方程求解，既减少了内存需求量，扩大了可模拟的模型规模，同时又使得PML边界处理更为简洁。本发明方法尤其适用于油气勘探领域中的
地震波场数值模拟或逆时偏移。
附图说明
[0024] 图1为本发明中基于PML吸收边界的三维声波数值模拟方法的流程示意图；
[0025] 图2为三维介质内部计算区域和PML吸收区域的示意图；
[0026] 图3为三维声波模拟的规则网格示意图；
[0027] 图4为三维声波模拟的交错网格示意图；
[0028] 图5为本发明中三维声波模拟的交错网格示意图，图中显示了XOY平面的变量分布。
具体实施方式
[0029] 下面结合附图以及具体实施方式对本发明作进一步详细说明：
[0030] 基于PML吸收边界的三维声波数值模拟方法，如图1所示，主要包括如下步骤：
[0031] 1、定义混合网格
[0032] 假设速度模型可离散为N×N×N的三维网格，各面的PML吸收层有M层。
[0033] 1.1在内部计算区域应用规则网格，波场压力P仅定义在规则网格整数节点上，需要N3个浮点类型变量；
[0034] 1.2在PML吸收区域应用交错网格，波场压力P定义在交错网格的整数节点上，质点的速度分量vx,vy,vz和介质密度，分别定义在X、Y、Z方向的半数节点上，需要6×M×N2×3个浮点类型数据。
[0035] 因此，混合网格共需要N3+18×M×N2个浮点类型变量。
[0036] 而完全采用规则网格需要的内存量为N3+42×M×N2个浮点数据，完全采用交错网3
格所需的内存量为4×(N+2×M) 个浮点数据。
[0037] 由此可见，本发明方法相对于上述两种传统方式，分别减少约24×M×N2和3×(2×M+1)×N3个浮点类型变量的存储空间，即内存需求量明显降低。
[0038] 以N＝500，M＝50为例，完全采用规则网格、完全采用交错网格及本发明方法的内存使用量分别为2.42G，3.22G，1.30G。通过对比可以发现，
[0039] 本发明方法相较于完全采用规则网格减少了约1.11G内存，相较于完全采用交错网格减少了约1.9G内存，相对于普通计算机的内存容量(2G-8G)，改善效果显著。
[0040] 2、不同区域波场的计算
[0041] 2.1、内部计算区域压力P值的计算
[0042] 利用二阶双曲型偏微分波动方程的有限差分离散格式，计算波场压力P值。
[0043] 2.2、PML吸收区域压力P值的计算
[0044] 在PML吸收区域采用三维声波的一阶应力-速度方程计算质点的速度分量vx,vy,vz和波场压力P值。具体的，
[0045] 首先结合内部计算区域最外层的波场压力P值，计算PML吸收区域的质点的速度分量vx,vy,vz值，然后进行PML吸收区域内波场压力P值的计算。
[0046] 至此，整个交错网格的一个时间步长的波场更新完成。另外，在PML吸收区域内最靠近内部计算区域的波场压力P值，可以用于更新下一时刻内部计算区域的波场压力P值。
[0047] 本发明方法简化了规则网格的PML边界处理，并且相对于交错网格，有效地减少了内存占用量，扩大了可模拟的模型规模。
[0048] 下面给出了本发明中三维声波数值模拟方法更为详细的实施过程：
[0049] 1定义混合网格
[0050] 首先根据给定的速度模型定义三维规则网格，为规则网格结点的波场压力P值分配三维存储空间；
[0051] 然后根据给定的PML吸收层数，在三维模型的周围六个面上为PML吸收区域，在其内定义交错网格；其中，交错网格的整数节点仍然存储波场压力P值，半数节点存储质点的速度分量vx,vy,vz和介质密度，如图5显示了XOY平面的变量分布：
[0052] 其中，粗线框内为规则网格区域，只需计算更新波场压力P值，粗线框外为PML吸收区域采用交错网格；位于网格半数节点的
水平方向的矩形格子表示质点速度的vx分量，竖直方向的矩形格子表示质点速度的vy分量。
[0053] 2计算模拟参数
[0054] 2.1定义子波主频，模拟时长；
[0055] 2.2根据频散要求，定义三个维度网格空间间隔；
[0056] 2.3根据
稳定性条件和空间差分
精度，确定模拟的时间步长，稳定性条件为：
[0057]
[0058] 其中，Δx、Δy、Δz为网格空间间隔，am为差分系数。
[0059] 3开始时间步长的循环
[0060] 3.1内部计算区域波场压力P值的计算
[0061] 此部分计算依据规则网格(如图3)，三维声波方程可表示为：
[0062]
[0063] 利用2N阶精度的差分方程近似微分方程可得到三维声波的离散形式，并
整理得到用于计算内部压力值的更新格式：
[0064] Pn+1(i,j,k)＝2Pn+1(i,j,k)-Pn+1(i,j,k)+Δt2v2{L2x[P(i,j,k)]+L2y[P(i,j,k)]+L2z[P(i,j,k)]}+s(is,js,ks,n+1)
[0065] 其中，
[0066]
[0067]
[0068]
[0069] 根据给定的
震源位置，将对应时刻的子波值加到相应的压力值上，即：
[0070] Pn+1(is,js,ks)＝P(is,js,ks)+s(is,js,ks,n+1)；
[0071] 其中，is,js,ks为震源所处的网格节点索引。
[0072] 3.2PML吸收区域质点的速度分量vx,vy,vz和波场压力P值的计算
[0073] 此部分计算基于交错网格(如图4)，在PML吸收区域，对应的一阶应力-速度方程为：
[0074]
[0075] 引入PML吸收条件后的控制方程可表示为：
[0076]
[0077] 其中,dx,dy,dz分别为X,Y,Z方向的衰减因子。
[0078] 利用2N阶差分离散代替微分，可得到PML吸收区域质点的速度分量vx,vy,vz和波场压力P值的更新格式，即：
[0079]
[0080]
[0081]
[0082] Pn+1/2(i+1/2,j,k)＝Pn-1/2(i+1/2,j,k)-ρv2[Lx(vx)+Ly(vy)+Lz(vz)]，[0083] 其中：
[0084]
[0085] 至此，整个交错网格的一个时间步长的波场更新完成。
[0086] 4输出波场值
[0087] 根据实际需要或者要求，输出某一时刻三维立体波场值或者某一平面切片值。
[0088] 5判断循环条件
[0089] 判断模拟时长是否满足循环要求，进入下一时刻的波场计算，重复步骤3；否则，退出程序，完成整个波场的计算。
[0090] 当然，以上说明仅仅为本发明的较佳
实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本
说明书的教导下，所做出的所有等同替代、明显
变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。
