一种模拟地震波在ti介质中传播规律的方法
【技术领域】
[0001] 本发明属于地震勘探领域，具体涉及一种模拟地震波在TI介质中传播规律的方 法，可模拟地震波在横向各向同性（TI)介质中传播规律，用于波动方程的正演模拟，来模 拟地震波在介质中的传播规律。
【背景技术】
[0002] 研究地震波在复杂介质中的传播规律，是地球物理学中的一项重要内容，是帮助 人们认识地球介质的有效模拟手段。
[0003] 地震波数值模拟是根据已知的地下结构和物性参数来模拟研究地震波在地下各 种介质中的传播规律，在地震勘探和天然地震领域有广泛的应用，而实际地下介质的各向 异性是普遍存在的.描述地下介质波的传播过程可以用波动方程，因此基于波动方程的 数值解法如有限差分法、有限元法、伪谱法、谱元法等在各向异性介质模拟中均有广泛应 用.这些方法都有各自的优缺点：有限差分法理论基础为基于Taylor展开，计算速度快， 效率高，但对于非规则区域问题较难处理；有限元法基于变分理论，模拟精度高，可以处理 复杂形状的地质体，但计算过程复杂，计算量大；伪谱法也是解决偏微分方程的数值解法之 一，在计算过程中引入了快速Fourier变换，精度可以达到谱的精度，但由于Fourier变换 需要全局信息，伪谱法在处理复杂介质中波场传播问题时显得力不从心；将有限元和伪谱 法优点结合起来的谱元法也存在着譬如边界和新方法难以匹配等问题。
【发明内容】
[0004] 本发明的目的在于解决上述现有技术中存在的难题，提供一种可用于TI介质正 演模拟的褶积微分方法，更好地描述波在介质中的传播特征，并降低地震波数值模拟中数 值频散现象。
[0005] 本发明是通过以下技术方案实现的：
[0006] -种模拟地震波在TI介质中传播规律的方法，包括：
[0007] 步骤1 :根据待求取的问题确定六阶矩阵系数：
[0008] 步骤2 :求解一阶速度--应力方程；
[0009] 步骤3 :求取步骤2中的褶积微分算子；
[0010] 步骤4 :求取窗函数系数，对褶积微分算子进行加窗截断；
[0011] 步骤5 :根据地质体，给定模型参数、震源参数、震源频率，对模型进行离散化；
[0012] 步骤6:进行时间迭代，得到整个区域内的速度和应力值，若设定时间到，则模拟 过程结束，输出波场快照和地震记录。
[0013] 所述步骤1是这样实现的：
[0014] 将TI介质各向异性参数表示为六阶矩阵形式：
[0015]
[0016] 其中C12 = Cn_2C66.式中：C为各向异性介质参数.在各向同性情况下，将该六阶 矩阵弱化为用拉梅常数λ和剪切模量μ来描述介质，此时各参数之间满足如下关系式：C n -C33 -入 +2 P·，C13 -入，C55 - P·。
[0017] 所述步骤2是这样实现的：
[0018] 考虑x-z平面的二维问题，速度-应力方程为
[0024] (6)
[0025] 式中，σ χχ、σ χζ、σ zz为应力分量；vx、Vz分别为速度水平分量和垂直分量，P为密 度，D x、Dz分别为空间沿水平和垂直方向的褶积微分算子。
[0026] 所述步骤3是这样实现的：
[0027] 将一个函数u (X)表示为奇异核函数与元函数的褶积形式：
[0028]
[0029] 用Shannon奇异核函数来逼近⑵式中的δ函数：
[0030]
[0031] 其中，Λχ为空间步长.当Λχ - 0时，所得函数即为Delta函数.
[0032] 对于一个函数的离散序列U(Xn)，通过对奇异核函数求导，得到该函数在全空间的 插值及微分形式：
[0033]
[0034] 其中表示最接近x的格点，m表示离散点的个数，q表示空间微分阶数，W为算 子半宽度.根据公式（2)、（3)、（4)，得到基于Shannon奇异核的一阶褶积微分算子为：
[0035]
[0036]
[0037]
[0038] 公式（5)和（6)得到的即为公式（1)中的离散形式的褶积微分算子D3^P Dz，其中 公式（5)适用于普通网格的地震波数值模拟，公式（6)适用于交错网格情况下的地震波数 值模拟。
[0039] 所述步骤6是这样实现的：按步骤2给出的公式计算应力及速度的离散值：其中 公式⑴中（l-a)、（l-b)、（l-c)为由速度分量v x、vz求取应力分量〇xx、〇xz、〇zz的公式， (Ι-d)、（Ι-e)为由应力分量σ χχ、σ χζ、σ zz求取速度分量vx、Vz的公式，这些变量的初始值 均赋为〇,步骤5中的震源参数与时间有关，随着时间迭代，通过这五个公式的迭代，得到整 个区域内的速度和应力值，若设定时间到，则模拟过程结束，输出波场快照和地震记录。
[0040] 与现有技术相比，本发明的有益效果是：
[0041] (1)更有利于描述波在介质中的传播特征；
[0042] (2)能够降低地震波数值模拟中数值频散现象。
【附图说明】
[0043] 图1是TI介质空间交错网格示意图。
[0044] 图2-1是优化褶积微分法水平分量。
[0045] 图2-2是优化褶积微分法垂直分量。
[0046] 图2-3是有限差分法水平分量。
[0047] 图2_4是有限差分法垂直分量。
[0048] 图3是有限差分方法与褶积微分方法在检波器位于（500m，400m)处的地震记录对 比图（X为水平分量，z为垂直分量）。
[0049] 图4是本方法的步骤框图。
【具体实施方式】
[0050] 下面结合附图对本发明作进一步详细描述：
[0051] 波动方程的空间微分可以采用多种形式，除了以上所述的有限差分、有限元、伪谱 法、谱元法等之外，本发明提出一种褶积微分算子对波动方程的空间进行微分。在20世纪 70年代早期，就有学者提出褶积微分算子的概念，但将其应用到地震波场模拟中比较晚。 Holberg运用Fourier变换，对褶积微分算子进行了设计，并在三维弹性波数值模拟中取得 了较好的结果。Zhou等、张中杰等和戴志阳等也把这一方法进行了深入研究和发展。龙桂 华等在前人的工作基础上，根据函数分布理论，引入奇异核褶积微分算子，通过一组经过加 窗处理的微分算子和波场变量进行褶积，来求取波场变量关于空间微分。褶积微分算子和 伪谱法相比，其优势在于褶积微分算子是一种优化短算子，更能突出空间微分本身的局部 属性；而伪谱法在计算微分时须利用全局信息，就复杂非均匀介质中波的传播问题而言，与 物理因果律相悖。褶积微分算子和有限差分算子相比，精度更高，能够弥补有限差分算子描 述复杂介质的不足。因此褶积微分算子应用于波动方程的空间离散或空间微分近似是非常 有效的。
[0052] 如图4所示，本发明是通过以下技术方案来解决的，以模拟TI介质中波的传播规 律为例，具体包含如下步骤：
[0053] 步骤1 :根据待求取的问题确定六阶矩阵系数：
[0054] TI介质各向异性参数可以表示为六阶矩阵形式：
[0055]
[0056] 其中C12 = Cn_2C66.式中：C为各向异性介质参数.在各向同性情况下，该六阶矩 阵可以弱化为用拉梅常数λ (弹性模量）和μ (剪切模量）来描述介质，此时各参数之间 满足如下关系式：C11 = C33 = λ +2 μ , C13 = λ，C55 = μ ·
[0057] 步骤2 :求解一阶速度--应力方程：
[0058] 考虑x-z平面的二维问题，速度-应力方程为
[0065] 式中，σ χχ、σ χζ、σ ζζ为应力分量；vx、Vz分别为速度水平分量和垂直分量，P为密 度，D X、DZ分别为空间沿水平和垂直方向的褶积微分算子.本发明提出的速度-应力方程与 常规速度-应力方程不同之处在于：利用褶积运算代替了常规速度-应力方程中的求导运 算。以公式（Ι-a)为例（以下公式中类似），D xVx表示对Vx在空间上沿水平方向进行褶积 运算，在常规速度-应力方程中，需要对V x在空间上沿水平方向求导，用数学表达式则记为 CV f。本发明提出的方法避免了空问上的求导运算。 VX
[0066] 公式（1)为应力与速度之间的关系，若使用程序来具体模拟计算时，上述公式中 所有变量均取离散格式。左边均为分量对时间求导，离散格式下可用普通时间差分进行计 算。
[0067] 数值模拟可在普通网格或者交错网格上进行。若采用普通网格形式，这些变量均 取网格点上的数值；若采用交错网格形式，这些变量有些需采用网格点上的数值，有些需采 用半网格点上的数值（即取周围网格点的平均值），本例中TI介质的各参数分布情况如图 1所示，其中实线与实线的交点为网格点，实线与虚线的交点为半网格点（取值为该点前 后两个点的平均值），虚线与虚线的交点也为半网格点（取值为该点周围的四个点的平均 值）。
[0068] 步骤3 :求取步骤2中的褶积微分算子。
[0069] 具体计算方法为：
[0070] 将一个函数u(X)根据函数分布理论，可表示为奇异核函数（可取Delta(S)函 数）与元函数（可取u (X)本身）的褶积形式：
[0071]
[0072]