一种基于弥散黏滞性波动方程的吸收边界反射方法
【技术领域】
[0001] 本发明属于地球物理勘探领域，涉及一种地震波数值模拟中的吸收边界处理方 法，特别涉及一种基于弥散黏滞性波动方程的吸收边界反射方法。
【背景技术】
[0002] 地震波数值模拟技术在油气勘探中有着重要的作用，它对于认识地下介质中波的 传播规律、验证所研宄的各种方法技术起关键性的作用。由于实际的地下介质是无限大 的，而数值模拟的计算区域是有限的，这样会造成由于计算区域被截断而产生人工边界，当 地震波传到该边界时会产生人工边界反射波。因此，需要在计算区域的边界处使用合理的 吸收边界条件进行处理，才能模拟实际地下介质中波的传播。关于人工边界反射的问题， 学者们已经发展了很多种方法，这些方法总体分为三类：第一类是预测法一利用波动方程 的不同形式如单程波方程、波场外推等近似计算边界处的波场值；第二类是阻尼衰减法一 利用合适的衰减函数将临近边界处的波场值进行衰减，在非均匀介质中由于波以不同的传 播速度沿着各个方向传播，因此，这类方法中衰减函数比较难确定；第三类是完全匹配层法 (PML一PerfectMatchedLayer)一计算区域的内部利用波动方程计算，而在计算区域的外 部需要加入适当的吸收层，并利用相应的带衰减函数的吸收边界控制方程进行计算，只要 吸收层内衰减函数选择的合理，这种方法就能很好地吸收各种频率来自各个方向的波。每 一种吸收边界反射的方法都有各自的优点和缺点，在第一类方法中，如Clayton-Engquist 方法和Higdon方法的优点是需要相对较少的存储量，但是这些方法只在有限的入射角范 围能够很好的吸收边界反射波；Lindman方法在较大的角度范围能达到很好的吸收效果， 但它只在边界是均匀的介质中有效；Liao方法能够在较宽的角度范围和非均匀介质边界 处有很好的效果，但它需要双精度实现来保证其稳定性，也需要更多的计算存储。第二类方 法中，如Kosloff方法在法向入射时比较有效，但对于非法向入射会产生明显的边界反射， 同时它也需要更大的存储空间。第三类PML方法是由Berenger在电磁波模拟中提出，它是 目前被认为效果最好的吸收边界条件，它的实现比较简单、稳定，能够适用于具有大范围入 射角和宽频带的入射波，广泛应用于很多领域，如线性Euler方程、声学介质波场模拟、弹 性波的传播、黏弹性介质中波的传播和孔隙弹性介质中波的传播等。
[0003] 弥散黏滞性波动方程是为了解释实际和实验中观察到的依赖于频率变化的反射 现象（如：低频伴影）而提出的，这种波动方程还可用于研宄在含有流体介质中地震波的衰 减和频散问题。因而，弥散黏滞性波动方程的求解在地震勘探领域有着重要的意义。利用弥 散黏滞性波动方程模拟地震波在含有流体介质中的传播过程也会遇到人工边界处理问题。
【发明内容】
[0004] 本发明的目的在于针对弥散黏滞性波动方程数值模拟中的人工边界反射问题，提 出一种基于弥散黏滞性波动方程的吸收边界反射方法，即在计算区域的外部加入适当的吸 收层，并利用相应的带衰减函数的吸收边界控制方程进行计算，计算区域的内部利用波动 方程计算。由于声波方程和Stokes方程是弥散黏滞性波动方程的特例，所以该方法也适用 于声波方程和Stokes方程。
[0005] 为了实现上述目的，本发明所采用的技术方案包括以下步骤：
[0006] 1)根据求解一维弥散性黏滞性波动方程的吸收边界控制方程，推广确定二维弥散 性黏滞性波动方程的吸收边界反射控制方程；其中，一维弥散性黏滞性波动方程为：
【主权项】
1. 一种基于弥散黏滞性波动方程的吸收边界反射方法，其特征在于，包括以下步骤： 1) 根据求解一维弥散性黏滞性波动方程的吸收边界控制方程，推广确定二维弥散性黏 滞性波动方程的吸收边界反射控制方程；其中，一维弥散性黏滞性波动方程为：
式中为波场函数；γ，Tl分别为弥散和黏滞性衰减系数，它们是岩石的孔隙度、渗透 率以及流体的密度、粘度等参数的函数；u为非频散介质中波的传播速度；x，t分别为空间 和时间变量；式（1)和式（2)中左端第一项表示惯性项，第二项为弥散耗损力即扩散项，第 三项表示黏滞性阻尼，第四项为波动方程的弹性部分； 二维弥散性黏滞性波动方程的吸收边界控制方程为：
式中M1(X), d2 (Z)分别为沿X，z方向的衰减函数； 2) 建立介质模型，在其周围加入吸收边界层，吸收边界层的厚度取值至少为一个波 长； 3) 选择吸收边界中的吸收衰减函数，其中吸收衰减函数包括指数函数、对数函数、正/ 余弦函数或其中的一种或几种的复合函数； 4) 对所计算的介质模型进行网格离散； 5) 利用有限差分法对二维弥散黏滞性波动方程进行数值模拟，在计算区域内部利用式 (4)弥散黏滞性波动方程的有限差分格式计算，在吸收边界层内利用式（3)计算，以消除人 工边界反射波；其中，式（4)如下：
C表示第η时间步在网格点（Xj, zm)处的波场值；η为时 间采样点，j为空间X方向的采样点，m为空间ζ方向的采样点； 震源函数采用爆炸源，在空间上采用高斯函数，时间上采用Ricker子波，震源函数的 形式为：
式中：&表示Ricker子波的中心频率，模型计算中β为常数；(X(l，Ztl)表示震源的空间 位置。
2.根据权利要求1所述的基于弥散黏滞性波动方程的吸收边界反射方法，其特征在 于：所述步骤1)中，根据求解一维弥散性黏滞性波动方程的吸收边界算法，推广确定二维 弥散性黏滞性波动方程的吸收边界反射控制方程的具体方法如下： 1-1) 一维弥散黏滞性波动方程的求解： 对式（1) 一维弥散黏滞性波动方程关于时间t做Fourier变换，整理后得到：
式中：是u(x, t)关于时间t的Fourier变换；ω为角频率； 式（5)办'/(Aj)关于X的常微分方程，其特征值λ为：
式中：
^为Xk的模；Θ为相位角，由arctan Θ = I/R确定；λ k的实部R和 虚部I分别为：
所以式⑴解的形式为：
式中C2为常数，由初始条件确定； 当衰减参数γ为零时，式（1)是即为Stokes方程，而这种情况下式（10)是即为Stokes 方程的解；当衰减参数γ，η均为零时，式（1〇)式即为一维声波方程的解； 将式（1)写成等价的偏微分方程组，然后构造与其相应的辅助方程组并引入坐标变 换，得到在新坐标系下的方程组，进而得到式（1)的解是指数衰减的，从而利用构建带衰减 函数的方程组用于吸收人工边界反射波； 与式（1)等价的一阶偏微分方程组为：
构建与式（1)和式（11)相对应的含有衰减函数的吸收边界控制方程组为：
氺 式中 ^为中间变量；sgn(t)为符号函数；表示时间域卷 积； 关于空间变量X引入坐标变换
将式（13)代入式（12)，再关于t做Fourier变换，得到：
式（15)和式（11)的形式完全相同，不同在于空间坐标JE与X，前者为复数，后者为实 数；由式⑴的解可知，式（15)的解为：
式中：c3, C4为常数，它们由初始条件确定； 从式（16)可以明显看出，方程组（15)的解沿X方向指数衰减，衰减形态和衰减大小由 d(x)决定；因此，只要恰当的选择该函数使得式（13)在计算区域内部为式（1)式的解，而 在边界处波函数呈现快速衰减，从而达到吸收边界反射波的目的； 1-2)二维弥散黏滞性波动方程吸收边界算法原理 与一维情况相同，引入中间变量Sx，Sz，Tx，Tz，R x，Rz，ux，Uz，并令
，则与式⑵二维弥散黏 滞性波动方程等价的偏微分方程组的形式为：
3.根据权利要求1所述的基于弥散黏滞性波动方程的吸收边界反射方法，其特征在 于：所述步骤3)中，吸收衰减函数采用正余弦函数，其形式为： Ci1(X) =M[1_cos(2jix/1)] d2(z) = M[l-cos (2 π z/1)] 式中：M为幅值常量，I为吸收边界层的厚度。
【专利摘要】本发明公开了一种基于弥散黏滞性波动方程的吸收边界反射的方法，为利用该方程进行复杂介质中地震波数值模拟提供了一种处理人工边界反射波的有效途径。首先从理论上阐述了一维和二维情况下该方法用于处理人工边界反射的原理；其次，在均匀介质和复杂储层介质模型中利用有限差分法对弥散黏滞性波动方程进行数值模拟，并采用所提的吸收边界反射方法进行处理，数值模拟结果验证了该方法的有效性。该方法也适用于声波方程和Stokes方程。此外，该技术方案易于实现，可操作性强。
【IPC分类】G01V1-28
【公开号】CN104749628
【申请号】CN201510145189
【发明人】高静怀, 赵海霞
【申请人】西安交通大学
【公开日】2015年7月1日
【申请日】2015年3月30日