基于有限差分法的三维tti双相介质地震波场数值模拟方法
【技术领域】
[0001] 本发明涉及地震勘探的技术领域，尤其涉及一种用平均入射角道集进行PP波与 PS波联合AVO反演方法。
【背景技术】
[0002] 石油和天然气是经济建设和人类生活不可缺少的重要能源，油气的供应直接影响 着世界各国的经济发展步伐和人们的生活水平，因此，寻找具有工业意义的油气储集层一 直是地球物理学家和地质学家的首要目标。我国从1993年开始成为石油进口国，并且进口 额逐年增大，如何摆脱国民经济对油气进口的依赖，提高国内油气产量，是我国油气工业所 面临的主要任务。
[0003] 油气产量的提高是以地质勘探查明的油气新储量与剩余油储量的增加为前提的， 这就要求地质工作者必须尽快找到具有工业生产价值的新的油气资源或剩余油气资源。石 油地震勘探的主要任务就是解决以上两个问题，为石油工业的持续发展服务，几十年来，我 国的地球物理工作者进行了大量的卓有成效的工作，在许多油气田的勘探与开发中发挥了 重要作用。但以往利用地震勘探资料找油气都属间接找油气范畴，即由物探工作者将野外 地震勘探资料的处理结果提交给地质工作者，再由地质工作者结合地质资料与其它资料确 定勘探区的油气分布与储量。这种方法的优点是多学科综合利用，信息量较大，缺点是钻探 成功率低，勘探周期长、成本高，经济效益差，并且物探资料中许多有效信息未能得到合理 应用，造成资源浪费。利用地震资料直接寻找油气可以大大缩短勘探周期，降低勘探成本， 提高钻探成功率。自从70年代初墨西哥湾利用地震亮点技术直接找气成功后，许多地球 物理学家投身于该领域相关技术研究并在理论上、方法上和实践中取得了一定突破，如AVO 技术、速度反演技术、地震属性技术等。上述技术在某些特定地区取得了好的效果，但在大 规模推广过程中遇到了困难，主要原因是这些技术都以传统的单相介质理论作为出发点， 没有充分考虑油气藏的双相特征，导致了预测或反演结果的较大误差甚至错误结果。采用 基于双相介质理论的油气检测技术可以克服上述缺陷，提高油气检测精度。
[0004] 双相介质理论是上世纪五十年代开始发展起来的一种新的地震波场理论。该理论 假设地下介质是由固相和流相（或气相）组成的。固相是指组成地下岩石的骨架颗粒，流相 是指充填于岩石孔隙中的流体。传统的地球物理方法常将地层介质加以简化，看成是纯固 体（单相介质）介质，这样做在岩石孔隙度很小或孔隙中只含有束缚流体时是成立的，而当 岩石孔隙度较大，而且孔隙中含有连续可动流体时，弹性理论简化就会有较大偏差，甚至得 不到正确结论。实际上，大多充满流体或气体的介质（双相介质或多相介质），由于固体介 质与流体介质的相互作用，使之与单相介质的物理一力学性质有巨大差异。例如，如果充填 流体是理想液体，则在物相分界面上切向应力为零；如果充填的是气体，则骨架与气体分界 面上的切向应力和法向应力均应消失。含油气地层实际上是具有固体状态与流体状态（气 体状态）的双相（三相）介质。实践发现，致密介质的经典模型不能很好的描述波在油层 中的传播过程，需要对它进行完善。双相介质理论充分考虑了介质的结构、流体与气体的特 殊性质、局部特性与整体效应的关系，双相介质模型更接近于实际，更能准确地描述实际地 层结构和地层性质，更能适应越来越复杂的油气储藏勘探的实际需要。因此，研究基于双相 介质理论的油气检测方法已显得非常必要。
[0005] 总之，双相介质的基本理论研究已趋于成熟，而其实际应用还处于探索阶段，目 前，这一理论在工业应用方面的成功实例还不多见，还没有达到成熟阶段，但理论的创新必 将带来技术上的革命，随着研究工作的进一步深入，双相介质中的弹性波理论必将在工业 生产中得到广泛的应用，并产生巨大的经济效益。
【发明内容】
[0006] 本发明的主要目的在于提供一种基于有限差分法的三维TTI双相介质地震波场 数值模拟方法，以实现物理地震波场的实时传播模拟。
[0007] 为解决上述问题，本发明实施例提供一种基于有限差分法的三维TTI双相介质地 震波场数值模拟方法，包括：取得地震波的固体应力张量、流体应力张量、固体应变张量和 流体应变张量；根据应力与应变的对应关系，将所述固体应力张量、流体应力张量、固体应 变张量和流体应变张量转换为所述地震波的本构方程式；根据应力与位移的对应关系，取 得所述地震波的几何方程式；根据所述本构方程式、所述几何方程式、流体相对于固体的运 动及应力与位移的对应关系，取得所述地震波的运动微分方程式；对运动微分方程两边取 散度，取得所述地震波的第一纵波方程式，并令第一纵波方程式中的耗散系数等于零，以取 得第二纵波方程式；对第一纵波方程式，令对y的偏导数等于零，对空间偏导数采用2N阶精 度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，取得第 一差分方程式，其中N为大于1的正整数；对第二纵波方程式，令对 y的偏导数等于零，对空 间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式 进行差分离散，取得第二差分方程式；对所述第一差分方程式与所述第二差分方程式进行 吸收边界条件处理，以取得对应的地震波场数值。
[0008] 根据本发明的技术方案，通过矩形体剖分、离散化，在时间上进行高阶近似，在边 界条件上使用交错网格的吸收边界条件，实现了固体相与流体相耦合作用下的双相介质数 值方程的迭代求解，即实现物理地震波场的实时传播模拟。
【附图说明】
[0009] 此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发 明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0010] 图Ia和图Ib分別是采用以上差分格式计算的均匀双相各向同性介质中固相波场 和流相波场的快照的示意图；
[0011] 图2a和图2b分別是双相各向同性介质中固相波场和流相波场的单炮记录的示意 图；
[0012] 图3为边界入射波的示意图；
[0013] 图4是吸收边界条件中角点处理方法的示意图；
[0014] 图5a和图5b分别为采用以上边界条件计算的均勾双相各向同性介质中固相波场 和流相波场的快照示意图；
[0015] 图6a和图6b分別为没有采用任何边界条件处理计算得到的均匀双相各向同性介 质中的固相波场和流相波场的快照示意图；
[0016] 图7a和图7b分别是利用上述方法计算的双相各向同性介质中X分量波场的固相 波场和流相波场的快照的示意图；
[0017] 图7c和图7d分别是利用上述方法计算的双相各向同性介质中z分量波场的固相 波场和流相波场的快照的示意图；
[0018] 图8是交错网格示意图；
[0019] 图9a和图9b分别是利用上述方法计算的双相各向同性介质中X分量波场的固相 波场和流相波场的快照的示意图；
[0020] 图9c和图9d分别是利用上述方法计算的双相各向同性介质中z分量波场的固相 波场和流相波场的快照的示意图；
[0021] 图IOa和图IOb分别是模型一的固相波场和流相波场的快照的示意图；
[0022] 图I Ia和图I Ib分别是模型二的固相波场和流相波场的快照的示意图；
[0023] 图12是两层介质的模型的示意图；
[0024] 图13a和图13b分别是模型三的固相波场和流相波场的快照的示意图；
[0025] 图14是交错网格的示意图；
[0026] 图15a和图15b分别是固相波场和流相波场的波场快照的示意图；
[0027] 图16a和图16b分别是固相波场和流相波场的波场快照的示意图；
[0028] 图17a和图17b分别是固相波场和流相波场的波场快照的示意图；
[0029] 图18a和图18b分别是固相波场和流相波场的波场快照的示意图；
[0030] 图19a和图19b分别是固相波场和流相波场的波场快照的示意图；
[0031 ] 图20a和图20b分别是固相波场和流相波场的波场快照的示意图；
[0032] 图2Ia和图2Ib分别是固相波场和流相波场的波场快照的示意图；
[0033] 图22和图23分别是X分量和z分量地震记录的不意图；
[0034] 图24和图25分别是采用二阶时间差分精度、十阶空间差分精度的波场快照的示 意图；
[0035] 图26是根据本发明实施例的基于有限差分法的三维TTI双相介质地震波场数值 模拟方法的流程图。
【具体实施方式】
[0036] 为使本发明的目的、技术方案和优点更加清楚，以下结合附图及具体实施例，对本 发明作进一步地详细说明。
[0037] 双相介质理论认为地下介质由固体和流体两部分组成，固体即为岩石的骨架，流 体为充填于骨架孔隙和裂缝中的液体或气体。当孔隙和裂缝中含有两种以上的流体时也称 为多相介质。地震波在固相与液相或固相与气相组成的岩石中传播时与其在单相介质中的 传播规律不同。由于固体介质与流体介质的相互作用，地震波在双相介质中的传播规律变 得更复杂。为方便，在研究双相介质时，一般忽略固体与流体之间的物理、化学作用以及它 们在高温高压环境下的热弹性效应。
[0038] 双相介质中的基本方程包括三组：本构方程、几何方程和运动微分方程。
[0039] I. 1本构方程
[0040] 双相介质中的应力张量分为两部分：固体应力张量和流体应力张量。固体应力张 量指作用在单位容积立方体每一面的固体部分上的应力，如公式（1. 1)所示：
[0042] 其中，σ χχ、σ yy、σ zz为固体相的正应力，τ xy、τ χζ、τ yx、τ yz、τ ζχ、τ zy为固体相 的切应力，且 Txy= τ yx，τχζ= τ ζχ，Tyz= τ zy〇
[0043] 流体应力张量指作用在该单位容积立方体每一面的流体部分上的应力，如公式 (1. 2)所示：
t
[0045] s满足公式（1. 3)，如下所示：
[0046] s = - Φ p (1. 3)
[0047] 其中，Φ为每单位截面中流体面积所占百分数，相当于有效孔隙率，p为流体压 力，负号表示应力s与流体压力p方向相反。由于流体内无切应力，所以应力张量只含有主 对对角分量。
[0048] 双相介质中的应变张量也分为两部分：固体应变张量和流体应变张量。
[0049] 作用于单元体每一固体截面的应变张量，如公式（1. 4)所示：
[0051] 其中，exx、eyy、ezz表不固相正应变，e xy、exz、eyx、eyz、ezx、e zy表不固相切应变，且e xy -6yX? θχζ - θ ζχ? θγζ - θ zyO
[0052] 作用于单元体每一流体截面的应变张量用下列矩阵表示，如公式（I. 5)所示：
[0054] 其中，ε表示流相的体应变。
[0055] 双相线性弹性介质的应力与应变具有线性关系，用广义Hooke定律表示，如公式 (1. 6)所示：
(L6)
[0057] 其中，Cljkl(i、j、k、l = 1，2,3,4)是固体相的弹性参数，R是流体相的弹性参数，Q =(Q1, Q2, Q3, Q4, Q5, 〇6厂是固相和流相耦合关系的弹性参数。
[0058] 当介质为各向同性介质时，广义Hooke定律变为公式（1. 7)，如下所示：
= C2222= C 3333= A+2N。所以，在双相各向同性介质中只有四个独立的弹性参数。此时应力 与应变的关系又可写为公式（1. 8)，如下所示：
[0062] 式中，Θ表示固相体应变，
u为固相 位移向量，ux、uy、Uz分别为固相位移向量u在X、y和z方向的分量；ε表示流相体应变， 11]为流相位移向量，1]：(、1^、1；分别表示流相位移向量1]在^7
' ? 和ζ方向的分量;A和N相当于单相各向同性弹性波理论中的拉梅系数，其中N = μ ;R表 示使一定体积的流体流入某集合体而又使该集合体保持总体积不变所需施加在流体上的 压力的一种量度；Q反映固体与流体体积变化之间的耦合