CN105676280A - 基于旋转交错网格的双相介质地质数据获取方法和装置 - Google Patents

Abstract

本发明提供了基于旋转交错网格的双相介质地质数据获取方法和装置，涉及地质数据测量领域。本发明提供的方法采用旋转交错网格对地质数据计算式进行处理的方式，与现有技术中使用一般交错网格对地质数据计算进行处理，导致在计算物理量(如速度和应力)的偏导时，需要通过内插法进行插值处理，才能够求得最终的结果，由于需要进行插值，因而导致最终求得的结果会有一定的偏差相比，其通过旋转交错网格对地质数据计算式进行处理，在后续模拟的时候，避免了场量和模型参数的插值问题，因而使得后续计算出的最终结果更为准确提高了计算精度和稳定性。

Description

基于旋转交错网格的双相介质地质数据获取方法和装置

技术领域

本发明涉及地质数据测量领域，具体而言，涉及基于旋转交错网格的双相介质地质数据获取方法和装置。

背景技术

随着工业技术的发展，人们对能源的需求量越来越大，常见的能源有电能、风能、化石能源等。日常生活中，如汽车中所使用的石油，即为化石能源的一种。不同种类能源的来源是不相同的，如电能通常是靠其他形式的能量转化得来，如风力发电、水力发电均是将动能转化为电能。与电能不同的是，化石能源是通过动植物深埋一定时间之后所转化的，并且需要通过勘探和挖掘来获得。

随着化石能源的开发，当前社会越来越多的使用了化石能源。为了供给充分的化石能源，高效开采化石能源的工作越发受到重视。化石能源的获取可以分为两个步骤，第一个步骤是勘探，来判断当地是否存在化石能源，或者说是了解指定地域的化石能源的情况；如果第一个步骤所得出的结果比较理想的话，就会执行第二个步骤，也就是通过开采的方式来得到化石能源。由此可见，如何准确的了解化石能源的情况是首要的。并不只是化石能源，所有深埋在地下的矿产资源在挖掘之前，都要先进行勘探、分析，进而确定该地域的地址情况，继而确定该地域是否值得进行挖掘。

在进行勘探的时候，通常会采用地震波测量技术，如通过在指定地域中埋设激发，进而通过激发形成地震波，再通过地震仪和检波器准确的对该地域的地震波进行采集，最后通过分析手段来分析出该地域的地质情况。

同时，矿产资源所埋藏的地域，通常是由多种不同的介质组成了储层(矿产资源即是埋藏在储层中的)。一般来说，储层介质多是由固体骨架和孔隙中的流体(油、气或水)组成的双相介质，而非单一的介质。由于不同介质的物理属性差异较大，因此，对于这种介质构造较为复杂的双相介质而言，储层内部多表现为各向异性的性质，这增加了研究储层构造和特性的难度。

考虑到双相介质构造的复杂性，为了更有效的对指定地域的地质情况进行确定，需要探究双相各向异性的介质环境中地震波的传播规律，来进一步探明该地域中储层的性质。

发明内容

本发明的目的在于提供基于旋转交错网格的双相介质地质数据获取方法和装置，以获取指定地域的应力和速度的准确性。

第一方面，本发明实施例提供了基于旋转交错网格的双相介质地质数据获取方法，包括：

根据双相介质运动方程，确定目标区域的地质数据计算式；

使用旋转交错网格对所述地质数据计算式进行离散化处理；

将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，以确定所述目标区域的优化地质数据，所述优化地质数据包括优化速度值和优化应力值。

结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，所述使用旋转交错网格对所述地质数据计算式进行离散化处理包括：

按照将速度和密度置于整网格点处，将弹性系数和应力置于半网格点处的方式，计算所述地质数据计算式中，初步地质数据；所述初步地质数据包括初步速度值和初步应力值。

结合第一方面，本发明实施例提供了第一方面的第二种可能的实施方式，其中，所述按照将速度和密度置于整网格点处，将弹性系数和应力置于半网格点处的方式，计算所述地质数据计算式中，初步速度值和初步应力值包括：

将速度和密度赋值于整网格点处，将弹性系数和应力赋值于半网格点处；

使用所述整网格点沿旋转交错网格的对角线的四个应力分量，通过中心差分计算的方式计算速度对时间的偏导数，即初步速度值；

使用所述半网格点沿旋转交错网格的对角线的四个速度分量，通过中心差分计算的方式计算应力对时间的偏导数，即初步应力值。

结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，所述将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟包括：

按照随地震波传播时间变化的顺序，依次循环迭代调整所述离散化处理后的地质数据计算式，每次迭代后，均保存迭代计算出的速度值和应力值，以确定优化速度值和优化应力值。

结合第一方面，本发明实施例提供了第一方面的第四种可能的实施方式，其中，所述地质数据计算式为：

${\mathrm{\ρ}}_{11}\frac{\∂u_i}{\∂t^2}+{\mathrm{\ρ}}_{12}\frac{\∂U_i}{\∂t^2}={\mathrm{\σ}}_{{\mathrm{ij}}^{\′}j}-b_{ij}(\frac{\∂U_j}{\∂t}-\frac{\∂u_j}{\∂t})$

${\mathrm{\ρ}}_{12}\frac{\∂u_i}{\∂t^2}+{\mathrm{\ρ}}_{22}\frac{\∂U_i}{\∂t^2}=s_{i_{\′}}+b_{ij}(\frac{\∂U_j}{\∂t}-\frac{\∂u_j}{\∂t});$

其中，i和j分别表示水平方向x,y或垂直方向z中的一个，u_j和U_j分别是固相和流相的位移在j方向的分量，b_ij为流体相对固体骨架运动时的耗散系数，σ_ij'j为作用在固相应力分量在j方向上的偏导，s为作用在流体单元侧面上的应力，ρ₁₁和ρ₂₂分别表示介质单位体积内固相和流相部分的有效质量，ρ₁₂为视质量，即流相相对固相运动时的质量耦合系数。

结合第一方面，本发明实施例提供了第一方面的第五种可能的实施方式，其中，所述地质数据计算式为：

其中，ν_i、V_i表示速度，σ_ij表示应力，c_ij为弹性系数，D_i表示关于密度的多项式，b_ij表示耗散系数，Q_i为固相和流相之间体积变化的耦合系数，R是描述孔隙流体的弹性参数。

结合第一方面，本发明实施例提供了第一方面的第六种可能的实施方式，其中，所述将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟包括：

采用如下稳定性条件对所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，

其中，Δt为时间步长，V_max是最大相速度，Δh是空间步长，c_k是交错网格空间差分系数，D是空间维数。

结合第一方面，本发明实施例提供了第一方面的第七种可能的实施方式，其中，所述将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟包括：

采用如下吸收边界条件对所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，

其中，s_x是扩展函数，a_x≥0，k_x≥1，a_x、k_x、d_x均为关于吸收边界宽度的修正参数，i为虚数单位，w为频率变量。

结合第一方面，本发明实施例提供了第一方面的第八种可能的实施方式，其中，所述离散化处理后的地质数据计算式为：

；

其中，

$\begin{array}{c}{W}_x^{-}\left(u{\left(x\right)}_{i,j}\right)=\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(x\right)}_{i+\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}-u{\left(x\right)}_{i-\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}}{2\mathrm{\Δ}x}\\ +\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(x\right)}_{i+\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}-u{\left(x\right)}_{i-\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}}{2\mathrm{\Δ}x};\end{array}$

$\begin{array}{c}{W}_z^{-}\left(u{\left(z\right)}_{i,j}\right)=\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(z\right)}_{i+\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}-u{\left(z\right)}_{i-\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}}{2\mathrm{\Δ}z}\\ +\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(z\right)}_{i+\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}-u{\left(z\right)}_{i-\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}}{2\mathrm{\Δ}z};\end{array}$

a为空间差分系数，u为表达式W示范参数，1≤n≤N，N表示空间差分阶数，Δx、Δz为沿水平方向X和垂直方向Z的空间差分步长，i、j分别表示所对应的网格点的位置。

第二方面，本发明实施例还提供了基于旋转交错网格的双相介质地质数据获取装置，包括：

确定模块，用于根据双相介质运动方程，确定目标区域的地质数据计算式；

离散化处理模块，用于使用旋转交错网格对所述地质数据计算式进行离散化处理；

模拟模块，用于将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，以确定所述目标区域的优化地质数据，所述优化地质数据包括优化速度值和优化应力值。

本发明实施例提供的基于旋转交错网格的双相介质地质数据获取方法，采用旋转交错网格对地质数据计算式进行处理的方式，与现有技术中使用一般交错网格对地质数据计算进行处理，导致在计算物理量(如速度和应力)的偏导时，需要通过内插法进行插值处理，才能够求得最终的结果，由于需要进行插值，因而导致最终求得的结果会有一定的偏差相比，其通过旋转交错网格对地质数据计算式进行处理，在后续模拟的时候，避免了场量和模型参数的插值问题，因而使得后续计算出的最终结果更为准确提高了计算精度。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1示出了传统交错网格的定义方式；

图2示出了使用传统的交错网格方法时，在地质异常体B边界出现了振幅异常的现象的示意图；

图3示出了本发明所提供的基于旋转交错网格的双相介质地质数据获取方法的基本流程图；

图4示出了本发明所提供的基于旋转交错网格的双相介质地质数据获取方法中，旋转交错网格的示意图；

图5示出了本发明所提供的基于旋转交错网格的双相介质地质数据获取方法中，均匀双相VTI介质旋转交错网格数值模拟波场快照；

图6示出了本发明所提供的基于旋转交错网格的双相介质地质数据获取方法中，含有CPML吸收边界前后波场快照对比示意图；

图7示出了本发明所提供的基于旋转交错网格的双相介质地质数据获取方法中，非均匀双相各向同性介质模型的示意图；

图8示出了本发明所提供的基于旋转交错网格的双相介质地质数据获取方法中，非均匀双相各向同性介质传统交错网格数值模拟波场快照；

图9示出了本发明所提供的基于旋转交错网格的双相介质地质数据获取方法中，非均匀双相各向同性介质旋转交错网格波场快照。

具体实施方式

下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

为了准确的确定双相介质构造的储层的性质，通常会采用地震波正演模拟的方式来确定目标参数。下面以常用的传统交错网格方法的解算过程来说明相关技术中如何确定目标参数。

传统的交错网格方法是将网格点分为整网格点和半网格点，在相邻网格的中间位置(即半网格点处)，计算导数值，将速度和应力分量分别定义在两个不同的时间层上，并且相邻两个时间层上的物理量在空间分布上恰好交错半个网格，从而达到了时间和空间上的交错的目的。

如图1所示，示出了传统交错网格的定义方式。其中，t表示时间，Δt为时间间隔，ν_x、ν_z分别表示速度沿X和Z方向的速度分量，σ_ij表示沿i方向的应力在j方向的分量，i、j分别表示水平方向X，Y和垂直方向Z不同分量中的一种。按照地震波正常传播规律，地震波在均匀情况下会以同样的速度向四周扩散，呈圆形对外传播；遇到岩性分界面时，会出现反射、透射的现象，进而产生反射波，透射波；遇到地层剧烈变化的地方，如断层断点、断棱、地层尖灭点等会产生绕射波。在使用传统交错网格方法进行非均匀双相各向同性介质波场模拟时，会出现振幅异常现象，不符合地震波动的正常传播规律。如图2所示，使用传统的交错网格方法时，在地质异常体B边界出现了振幅异常的现象(方框圈出的位置)，且没有绕射波产生，很明显违背了地震波的传播规律。如此，在特殊情况下不能适用，稳定性不足。下面以差分格式进行说明，参见公式1和公式2。

$\frac{\∂{\mathrm{\υ}}_y}{\∂t}=(D_2+D_3)b_{22}({\mathrm{\υ}}_y-V_y)-D_3(\frac{\∂{\mathrm{\σ}}_{xy}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{yz}}{\∂z})$ 公式1；

$\frac{\∂{\mathrm{\σ}}_{xy}}{\∂t}=c_{16}\frac{\∂{\mathrm{\υ}}_x}{\∂x}+c_{36}\frac{\∂{\mathrm{\υ}}_z}{\∂z}+c_{46}\frac{\∂{\mathrm{\υ}}_y}{\∂z}+c_{56}(\frac{\∂{\mathrm{\υ}}_x}{\∂z}+\frac{\∂{\mathrm{\υ}}_z}{\∂x})+c_{66}\frac{\∂{\mathrm{\υ}}_y}{\∂x}$ 公式2；

其中，ν_i、V_i表示速度，σ_ij表示应力，c_ij为弹性系数，D_i表示关于密度的多项式，b_ij表示耗散系数，Q_i为固相和流相之间体积变化的耦合系数。

在二维三分量双相介质波场模拟时，由速度—应力方程中速度Y分量方程(1)、(2)可知，按照传统交错网格定义方式，计算速度分量导数时，由应力速度相对空间位置关系可知，需要计算应力σ_xy对水平方向导数，故需对σ_xy进行插值；计算应力分量导数时，由应力速度相对空间位置关系可知，需要计算速度v_x,v_y,v_z对水平和垂直两个方向的导数，故需对v_x,v_y,v_z进行插值(其它分量可按照该方式，类似得出)。

由于上述使用传统交错网格进行求取时，会受到插值因素的影响，导致传统双相介质波场模拟的精度大大降低，波场数据保真度受到影响。有鉴于此，本申请提供了基于旋转交错网格的双相介质地质数据获取方法，本申请的方法核心是利用旋转交错网格替代了传统的交错网格，再进行二维三分量双相介质波场模拟，以获得高保真的波场数据(速度和应力)。

首先对本申请所提供方法的特点进行说明，在交错网格波场模拟中，对交错网格进行了重新定义，即旋转交错网格方法。如图4所示，图4中i，j分别表示x，y和z不同分量中的一个，ν表示速度，ρ表示密度，C表示弹性系数。旋转交错网格是在交错网格的基础上对网格进行了旋转，将同一物理量定义在同一个网格点上。其主要思想是：首先计算沿网格对角线物理量的差分，然后将计算结果通过坐标变换得到沿坐标轴的差分。旋转交错网格有效地避免了场量及模型参数的插值问题，定义方式更加灵活、有效。

下面简要说明本申请的三个基本步骤，如图3所示，包括如下三个步骤：

S101，根据双相介质运动方程，确定目标区域的地质数据计算式；

S102，使用旋转交错网格对地质数据计算式进行离散化处理；

S103，将离散化处理后的地质数据计算式在目标区域所对应的模型下进行数值模拟，以确定目标区域的优化地质数据，优化地质数据包括优化速度值和优化应力值。

其中，步骤S101，主要是基于Biot双相介质理论，得到双相各向异性介质的运动方程，根据双相介质的运动方程，应用几何方程可进行化简，得到双相各向异性介质速度-应力方程表达式(二维三分量差分方程，即，地质数据计算式)。如公式3所示。

${\mathrm{\ρ}}_{11}\frac{\∂u_i}{\∂t^2}+{\mathrm{\ρ}}_{12}\frac{\∂U_i}{\∂t^2}={\mathrm{\σ}}_{{\mathrm{ij}}^{\′}j}-b_{ij}(\frac{\∂U_j}{\∂t}-\frac{\∂u_j}{\∂t})$

${\mathrm{\ρ}}_{12}\frac{\∂u_i}{\∂t^2}+{\mathrm{\ρ}}_{22}\frac{\∂U_i}{\∂t^2}=s_{i_{\′}}+b_{ij}(\frac{\∂U_j}{\∂t}-\frac{\∂u_j}{\∂t})$ 公式3；

其中，i和j分别表示水平方向x,y或垂直方向z中的一个，u_j和U_j分别是固相和流相的位移在j方向的分量，b_ij为流体相对固体骨架运动时的耗散系数，σ_ij'j为作用在固相应力分量在j方向上的偏导，s为作用在流体单元侧面上的应力，ρ₁₁和ρ₂₂分别表示介质单位体积内固相和流相部分的有效质量，ρ₁₂为视质量，即流相相对固相运动时的质量耦合系数。

步骤S102，采用了新的网格定义方式替代了传统的交错网格，新的网格定义方式即是旋转交错网格。网格定义参见图4所示。具体过程是：在双相介质速度-应力方程式的基础上，使用旋转交错网格对速度-应力方程进行离散化处理；将方程式和网格定义相互对应，将速度—应力方程中的速度和密度置于整网格点处，弹性系数和应力置于半网格点处；按照旋转交错网格的定义方式，计算整网格点处速度对时间的偏导数时，使用该点沿对角线的四个应力分量，通过中心差分的方式计算空间四阶导数，即速度对时间的偏导数；计算半网格点处应力对时间的偏导数时，使用该点沿对角线的四个速度分量，通过中心差分的方式计算空间四阶导数，即应力对时间的偏导数。由以上过程确定整个速度—应力方程(地质数据计算式)的速度值(即，初步速度值)和应力值(初步应力值)。

步骤S103，利用离散化后得到的双相介质速度-应力方程在特定模型下进行数值模拟，进而获得最终的优化地质数据(包括优化速度值和优化应力值)。具体过程是：通过速度和应力的空间四阶导数，转化成速度分量和应力分量的求取；随着地震波传播时间的延续，依次循环迭代速度和应力的计算公式，迭代的过程中，先由应力值求取速度，再由速度值求取应力，达到时间和空间相互交错的目的；以速度场值的分布表示地震波场分布，将每个时刻迭代计算的速度值和应力值保存下来，获得高保真地地震波场数据(即优化地质数据)。

下面以就本申请所提供的方法进行整体性说明。

首先说明基本的地质数据计算式的获取过程：

由Biot给出的饱和流体孔隙介质地震波传播理论，可以得到双相TTI介质的运动方程：

${\mathrm{\ρ}}_{11}\frac{\∂u_i}{\∂t^2}+{\mathrm{\ρ}}_{12}\frac{\∂U_i}{\∂t^2}={\mathrm{\σ}}_{{\mathrm{ij}}^{\′}j}-b_{ij}(\frac{\∂U_j}{\∂t}-\frac{\∂u_j}{\∂t})$

${\mathrm{\ρ}}_{12}\frac{\∂u_i}{\∂t^2}+{\mathrm{\ρ}}_{22}\frac{\∂U_i}{\∂t^2}=s_{i_{\′}}+b_{ij}(\frac{\∂U_j}{\∂t}-\frac{\∂u_j}{\∂t})$ 公式4；

公式4中，i、j分别表示x,y,z三个不同方向分量中的一个，u_j和U_j分别是固相和流相的位移在j方向的分量，b_ij为流体相对固体骨架运动时的耗散系数，σ_ij'j为作用在固相应力分量在j方向上的偏导，s为作用在流体单元侧面上的应力。假设φ为孔隙度，ρ₁₁和ρ₂₂分别表示介质单位体积内固相和流相部分的有效质量，ρ₁₂为视质量，即流相相对固相运动时的质量耦合系数。三者与固体密度ρ_s和流体密度ρ_f之间满足如下公式5，

ρ₁₁+ρ₁₂＝(1-φ)ρ_s

ρ₁₂+ρ₂₂＝φρ_f公式5；

其中，φ为孔隙度，ρ₁₁和ρ₂₂分别表示介质单位体积内固相和流相部分的有效质量，ρ₁₂为视质量，即流相相对固相运动时的质量耦合系数。

结合双相介质的柯西方程、本构方程和运动微分方程，可以得到双相TTI介质二维三分量一阶速度—应力方程。设固相介质速度向量为ν＝(υ_x,υ_y,υ_z)^Τ,流相介质速度向量为V＝(V_x,V_y,V_z)^Τ,令,则速度和应力分量分别表示为：

公式6

公式6中， $D_1={\mathrm{\ρ}}_{11}/({\mathrm{\ρ}}_{12}^2-{\mathrm{\ρ}}_{11}\·{\mathrm{\ρ}}_{22}),D_2={\mathrm{\ρ}}_{12}/({\mathrm{\ρ}}_{12}^2-{\mathrm{\ρ}}_{11}\·{\mathrm{\ρ}}_{22}),D_3={\mathrm{\ρ}}_{22}/({\mathrm{\ρ}}_{12}^2-{\mathrm{\ρ}}_{11}\·{\mathrm{\ρ}}_{22}),$ ν_i、V_i表示速度，Q_i表示固相和流相之间体积变化的耦合系数，R是描述孔隙流体的弹性参数，表示强迫一定体积流体流入双相介质体积元，而又保持总体积恒定不变时需要作用在流体上的一种力的度量，c_ij表示虎克定律中弹性系数。需要说明的是，公式3是公式6的简化形式。

其次，需要说明旋转交错网格有限差分方法及差分格式。

旋转交错网格有限差分技术是在传统的交错网格基础上发展来的。传统交错网格有限差分技术是将网格点分为整网格点和半网格点，在相邻网格的中间位置即半网格点处，计算导数值，将速度和应力分量分别定义在两个不同的时间层上，并且相邻两个时间层上的物理量在空间分布上恰好交错半个网格，从而达到了时间和空间上的交错的目的。如图1所示，实线表示空间分布上的整网格，虚线表示半网格。

交错网格有限差分技术已经广泛应用于地震波正演模拟中，但当地质体中存在密度、速度极小的异常体时，其稳定性将不再满足。而旋转交错网格方法在各向异性介质正演模拟时，无需对弹性模量进行插值，避免了由于插值带来的误差，一定程度上降低了不稳定性。

旋转交错网格是在交错网格的基础上对网格进行了旋转，将同一物理量定义在同一个网格点上。其主要思想是：首先计算沿网格对角线物理量的差分，然后将计算结果通过坐标变换得到沿坐标轴的差分。旋转交错网格有效地避免了场量及模型参数的插值问题，定义方式更加灵活、有效。如图4所示，i,j分别表示x,y,z不同分量中的一个。

由双相TTI介质的速度—应力方程(公式6)，并结合旋转交错网格有限差分方法可得双相TTI介质弹性波二维三分量旋转交错网格有限差分格式如下：

公式7；

公式7中，

$\begin{array}{c}{W}_x^{-}\left(u{\left(x\right)}_{i,j}\right)=\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(x\right)}_{i+\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}-u{\left(x\right)}_{i-\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}}{2\mathrm{\Δ}x}\\ +\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(x\right)}_{i+\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}-u{\left(x\right)}_{i-\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}}{2\mathrm{\Δ}x};\end{array}$

$\begin{array}{c}{W}_z^{-}\left(u{\left(z\right)}_{i,j}\right)=\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(z\right)}_{i+\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}-u{\left(z\right)}_{i-\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}}{2\mathrm{\Δ}z}\\ +\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(z\right)}_{i+\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}-u{\left(z\right)}_{i-\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}}{2\mathrm{\Δ}z};\end{array}$

a为空间差分系数，u为表达式W示范参数，1≤n≤N，N表示空间差分阶数，Δx、Δz为沿水平方向X和垂直方向Z的空间差分步长，i、j分别表示所对应的网格点的位置。

最后，说明两个约束条件，第一个是稳定性条件。

稳定性是所有数值模拟方法所必须考虑的问题。在旋转交错网格有限差分技术中，Saenger在Neumann稳定性条件下，给出了在空间步长相等时，时间二阶、空间2N阶的旋转交错网格有限差分的稳定性条件：

$\frac{{\mathrm{\ΔtV}}_{max}}{\mathrm{\Δ}h}\≤1/\left(\underset{k=1}{\overset{n}{\Σ}}\right|c_k\left|\right)$ 公式8；

式中，Δt为时间步长，V_max是最大相速度，Δh是空间步长，c_k是交错网格空间差分系数。

对于传统交错网格有限差分技术，在步长相等的情况下，其时间二阶、空间2N阶的传统交错网格有限差分的稳定性条件：

$\frac{{\mathrm{\ΔtV}}_{max}}{\mathrm{\Δ}h}\≤1/\left(\sqrt{D}\underset{k=1}{\overset{n}{\Σ}}\right|c_k\left|\right)$ 公式9；

公式9中，Δt为时间步长，V_max是最大相速度，Δh是空间步长，c_k是交错网格空间差分系数，D是空间维数，通过对比公式8和9可知，旋转交错网格相较于传统交错网格，理论上有更宽松的稳定性条件，这也是旋转交错网格的一个优点。

第二个是吸收边界条件。

为了更准确的模拟半空间介质，吸收边界条件是一个必须要面临的问题。本文将采用CPML进行边界处理。CPML技术是在完全匹配层的基础上发展而来的，通过引入记忆变量和时间上的卷积来达到不分裂的目的，从而避免了常规PML分裂带来的复杂性。

CPML对扩展函数s_x的形式进行了如下改造：

$s_x=k_x+\frac{d_x}{a_x+iw}$ 公式10；

其中，a_x≥0，k_x≥1。当a_x＝0，k_x＝1，s_x即为常规的PML标准形式(扩展函数)，，a_x、k_x、d_x均为关于吸收边界宽度的修正参数，i为虚数单位，w为频率变量。通过对空间导数做时间域上的卷积变换来推导CPML吸收边界公式。

下面对均匀双相VTI介质旋转交错网格数值模拟进行说明：

设定均匀VTI介质模型，模型大小为300m×300m，网格间距dx＝dz＝1m，采样间隔为0.1ms，震源采用固相x方向水平偏振波源激发的Ricker子波，位于中心位置，震源主频为80Hz。物性参数见表1，介质极化角为0°，方位角为0°，取t＝40ms时的波场快照，如图5所示。

表1均匀双相VTI介质物性参数

从图5中可以清晰地观察到，双相VTI介质中存在三种弹性波：快纵波qP¹、慢纵波qP²、快横波qS¹；从图5b、5d中可以明显看到，慢纵波在流相中相对固相能量更强，更容易观察，且相对于流相中其它波场，能量亦更强；在VTI介质数值模拟中，由于固相的各向异性，Y分量没有产生波场；固相为VTI介质时，观察不到横波分裂现象，仅有快横波存在。

为验证不分裂完全匹配层的吸收效果，选取均匀双相HTI介质进行波场模拟。设定均匀HTI介质模型，计算模型大小为260m×260m，其中四周吸收边界宽度均为30m，网格间距dx＝dz＝1m，采样间隔为0.1ms，震源同样采用固相x方向水平偏振波源激发的Ricker子波，位于中心位置，震源主频为80Hz。物性参数同表2，取t＝60ms时的波场快照，吸收边界使用前后波场快照如图6所示。

表2均匀双相TTI介质物性参数

下面对传统的交错网格和旋转交错网格对比说明(主要涉及波场模拟的精度和稳定性)：

一、稳定性方面：

设定非均匀双相各向同性介质模型，模型大小为256m×256m，网格间距dx＝dz＝1m，模型包含A、B两块区域，其中B为地质异常体，是一条密度和速度极小的裂缝，如图7所示。采样间隔为0.1ms，震源采用固相x方向水平偏振波源激发的Ricker子波，位于中心位置，震源主频为120Hz。物性参数见表3，取t＝40ms时的波场快照，如图8、9所示。

表3非均匀双相各向同性介质物性参数

从图8中可以观察到，当使用传统交错网格有限差分技术对非均匀极端介质模拟时，由于地质异常体的存在，在B的边界经过多次插值后，导致地震波场在地质异常体B的两端产生振幅异常现象，不符合地震波传播规律，不能正确反映地震波传播情况，故对于非均匀的极端介质，传统交错网格有限差分技术不能适用。从图9中可以观察到，当使用旋转交错网格有限差分技术对同样的介质模型模拟时，可以清晰的观察到：在地质异常体B的两端产生明显的绕射波，两侧出现了反射波，得到的波场快照符合地震波传播规律。

通过图8和图9的对比可知：对于存在地质异常体的非均匀各向同性介质，传统交错网格有限差分技术不能适用，其稳定性明显不足，而旋转交错网格有限差分技术则依然适用，稳定性较好。

二、精度方面：

对于二维双相各向异性介质，两种交错网格及对应的弹性波波场分量和弹性参数的空间位置如图1、表4和表5所示。

表4传统交错网格中弹性波场分量和弹性参数位置

表5旋转交错网格中弹性波场分量和弹性参数位置

结合公式11、12，并通过图1，表4及表5可知：在双相各向异性介质的数值模拟中，当使用传统交错网格有限差分技术时，由于速度和应力的相对空间位置关系，计算时需要对应力σ_xy和速度v_x,v_y,v_z进行插值，从而增加了计算误差，降低了计算精度；当使用旋转交错网格有限差分技术时，由于速度分量和应力分量各定义整网格点和半网格点上，计算空间导数时不需要进行任何插值，从而降低了计算误差，提高了计算精度。下面以速度分量和应力分量方程加以说明。

$\frac{\∂{\mathrm{\υ}}_y}{\∂t}=(D_2+D_3)b_{22}({\mathrm{\υ}}_y-V_y)-D_3(\frac{\∂{\mathrm{\σ}}_{xy}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{yz}}{\∂z})$ 公式11；

$\frac{\∂{\mathrm{\σ}}_{xy}}{\∂t}=c_{16}\frac{\∂{\mathrm{\υ}}_x}{\∂x}+c_{36}\frac{\∂{\mathrm{\υ}}_z}{\∂z}+c_{46}\frac{\∂{\mathrm{\υ}}_y}{\∂z}+c_{56}(\frac{\∂{\mathrm{\υ}}_x}{\∂z}+\frac{\∂{\mathrm{\υ}}_z}{\∂x})+c_{66}\frac{\∂{\mathrm{\υ}}_y}{\∂x}$ 公式12；

结合图1，表4及表5可知：在式11，当使用传统交错网格技术计算速度分量导数时，由应力速度相对空间位置关系可知，需要计算应力σ_xy对水平方向导数，故需对σ_xy进行插值；在式12，当使用传统交错网格技术计算应力分量导数时，由应力速度相对空间位置关系可知，需要计算速度v_x,v_y,v_z对水平和垂直两个方向的导数，故需对v_x,v_y,v_z进行插值。其它分量可类似得出。

由以上两个方面的对比可知：传统交错网格有限差分技术在一定条件下会表现出稳定性不足的情况，使用范围有限，由于插值的需要，会增大计算误差，降低计算精度；旋转交错网格有限差分技术表现出更强的适用性，在计算过程中无需插值，提高了计算精度。

本申请实施例还提供了与上述基于旋转交错网格的双相介质地质数据获取方法相对应的基于旋转交错网格的双相介质地质数据获取装置，包括：

确定模块，用于根据双相介质运动方程，确定目标区域的地质数据计算式；

离散化处理模块，用于使用旋转交错网格对所述地质数据计算式进行离散化处理；

模拟模块，用于将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，以确定所述目标区域的优化地质数据，所述优化地质数据包括优化速度值和优化应力值。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统、装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (10)

1.基于旋转交错网格的双相介质地质数据获取方法，其特征在于，包括：

根据双相介质运动方程，确定目标区域的地质数据计算式；

使用旋转交错网格对所述地质数据计算式进行离散化处理；

将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，以确定所述目标区域的优化地质数据，所述优化地质数据包括优化速度值和优化应力值。

2.根据权利要求1所述的基于旋转交错网格的双相介质地质数据获取方法，其特征在于，所述使用旋转交错网格对所述地质数据计算式进行离散化处理包括：

按照将速度和密度置于整网格点处，将弹性系数和应力置于半网格点处的方式，计算所述地质数据计算式中，初步地质数据；所述初步地质数据包括初步速度值和初步应力值。

3.根据权利要求2所述的基于旋转交错网格的双相介质地质数据获取方法，其特征在于，所述按照将速度和密度置于整网格点处，将弹性系数和应力置于半网格点处的方式，计算所述地质数据计算式中，初步速度值和初步应力值包括：

将速度和密度赋值于整网格点处，将弹性系数和应力赋值于半网格点处；

使用所述整网格点沿旋转交错网格的对角线的四个应力分量，通过中心差分计算的方式计算速度对时间的偏导数，即初步速度值；

使用所述半网格点沿旋转交错网格的对角线的四个速度分量，通过中心差分计算的方式计算应力对时间的偏导数，即初步应力值。

4.根据权利要求1所述的基于旋转交错网格的双相介质地质数据获取方法，其特征在于，所述将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟包括：

按照随地震波传播时间变化的顺序，依次循环迭代调整所述离散化处理后的地质数据计算式，每次迭代后，均保存迭代计算出的速度值和应力值，以确定优化速度值和优化应力值。

5.根据权利要求1所述的基于旋转交错网格的双相介质地质数据获取方法，其特征在于，所述地质数据计算式为：

${\mathrm{\ρ}}_{11}\frac{\∂u_i}{\∂t^2}+{\mathrm{\ρ}}_{12}\frac{\∂U_i}{\∂t^2}={\mathrm{\σ}}_{{\mathrm{ij}}^{\′}j}-b_{ij}(\frac{\∂U_j}{\∂t}-\frac{\∂u_j}{\∂t})$

${\mathrm{\ρ}}_{12}\frac{\∂u_i}{\∂t^2}+{\mathrm{\ρ}}_{22}\frac{\∂U_i}{\∂t^2}=s_{{\′}_i}+b_{ij}(\frac{\∂U_j}{\∂t}-\frac{\∂u_j}{\∂t});$

其中，i和j分别表示水平方向x,y或垂直方向z中的一个，u_j和U_j分别是固相和流相的位移在j方向的分量，b_ij为流体相对固体骨架运动时的耗散系数，σ_ij'j为作用在固相应力分量在j方向上的偏导，s为作用在流体单元侧面上的应力，ρ₁₁和ρ₂₂分别表示介质单位体积内固相和流相部分的有效质量，ρ₁₂为视质量，即流相相对固相运动时的质量耦合系数。

6.根据权利要求1所述的基于旋转交错网格的双相介质地质数据获取方法，其特征在于，所述地质数据计算式为：

其中，ν_i、V_i表示速度，σ_ij表示应力，c_ij为弹性系数，D_i表示关于密度的多项式，b_ij表示耗散系数，Q_i为固相和流相之间体积变化的耦合系数，R是描述孔隙流体的弹性参数。

7.根据权利要求4所述的基于旋转交错网格的双相介质地质数据获取方法，其特征在于，所述将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟包括：

采用如下稳定性条件对所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，

其中，Δt为时间步长，V_max是最大相速度，Δh是空间步长，c_k是交错网格空间差分系数，D是空间维数。

8.根据权利要求4所述的基于旋转交错网格的双相介质地质数据获取方法，其特征在于，所述将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟包括：

采用如下吸收边界条件对所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，

其中，s_x是扩展函数，a_x≥0，k_x≥1，a_x、k_x、d_x均为关于吸收边界宽度的修正参数，i为虚数单位，w为频率变量。

9.根据权利要求1所述的基于旋转交错网格的双相介质地质数据获取方法，其特征在于，所述离散化处理后的地质数据计算式为：

；

其中，

$\begin{array}{c}{W}_x^{-}\left(u{\left(x\right)}_{i,j}\right)=\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(x\right)}_{i+\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}-u{\left(x\right)}_{i-\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}}{2\mathrm{\Δ}x}\\ +\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(x\right)}_{i+\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}-u{\left(x\right)}_{i-\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}}{2\mathrm{\Δ}x}\end{array};$

$\begin{array}{c}{W}_z^{-}\left(u{\left(z\right)}_{i,j}\right)=\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(z\right)}_{i+\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}-u{\left(z\right)}_{i-\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}}{2\mathrm{\Δ}z}\\ +\underset{n=1}{\overset{N}{\Σ}}{a}_n^N\frac{u{\left(z\right)}_{i+\frac{(2n-1)}{2},j-\frac{(2n-1)}{2}}-u{\left(z\right)}_{i-\frac{(2n-1)}{2},j+\frac{(2n-1)}{2}}}{2\mathrm{\Δ}z}\end{array};$

a为空间差分系数，u为表达式W示范参数，1≤n≤N，N表示空间差分阶数，Δx、Δz为沿水平方向X和垂直方向Z的空间差分步长，i、j分别表示所对应的网格点的位置。

10.基于旋转交错网格的双相介质地质数据获取装置，其特征在于，包括：

确定模块，用于根据双相介质运动方程，确定目标区域的地质数据计算式；

离散化处理模块，用于使用旋转交错网格对所述地质数据计算式进行离散化处理；

模拟模块，用于将所述离散化处理后的地质数据计算式在所述目标区域所对应的模型下进行数值模拟，以确定所述目标区域的优化地质数据，所述优化地质数据包括优化速度值和优化应力值。