CN102183790A - 基于时空双变网格的弹性波正演模拟技术 - Google Patents

Abstract

本发明提供基于时空双变网格的弹性波正演模拟技术，针对已知的地震地质模型，根据地下介质的速度参数和弹性参数确定目标区域和背景区域的稳定性条件，从而计算出不同尺度的空间网格大小并对速度场进行离散化处理；然后利用泰勒公式和有限差分技术得到不同尺度网格点上变量的空间任意精度差分公式，即空间常网格算子和空间变网格算子的任意偶数阶精度空间差分公式，利用推导出的空间差分公式计算精细剖分区域、粗糙剖分区域和中间过渡带区域的波场值；在时间层的波场外推过程中，对不同网格剖分区域采用不同的时间采样间隔满足局部稳定性条件，分别得到精细时间层内精细网格点处的波场值以及全局时间层内粗糙网格点处的波场值；在每一个全局时间层内，将精细网格剖分区域的波场值传递到粗糙网格剖分区域，得到整个速度场在此时刻的空间波场值；最后利用上述方法在每个时间层进行波场外推，得到单炮的正演模拟记录；满足了浅层的局部稳定性条件，从而避免了频散噪音的产生。

Description

基于时空双变网格的弹性波正演模拟技术

所属技术领域

本发明属于地震勘探基础应用领域，是一种针对地下复杂介质的高效地震波正演模拟技术。

背景技术

随着地震勘探开发的深入，勘探领域逐步过渡到复杂地表、复杂地质结构的双复杂构造环境。此时地震勘探的成效，在很大程度上取决于所研发技术对复杂实际介质模型的适应性。因此在进行基于实际模型的地震波正演模拟时，模型的地表和深层复杂性是不可避免的。传统有限差分地震波模拟方法基于规则的矩形网格，是处理非均匀介质，特别是连续非均匀介质波动问题的有效方法。但是，当模型存在低降速层或小型非均质体时，网格间距必须取得较小，以保证计算精度和稳定性，但是这种规则的网格离散方法极大的增加了计算量，降低了计算效率。

地震波场正演过程中，模型网格划分的优劣与波场正演的效率密切相关。网格划分可分为规则网格剖分和不规则网格剖分，而简单的不规则网格可以通过沿x轴和z轴取不同的网格间距来实现。早期的变网格算法是在规则的矩形网格基础上实现的，在精细网格和粗糙网格的过渡区域往往要使用插值算法来得到每一时刻的波场值。Jastram和Behle(1992)提出了基于二维声波方程中某一深度变网格步长的算法；Jastram和Tessemer(1994)用了一种类似的方法针对二维弹性波模型，提出了在垂向上网格步长逐渐变化的算法。随后Falk(1996)将这一方法用于模拟井间地震中产生的管波。Arben Pitarka(1999)实现了低阶精度的不规则交错网格正演模拟，并首次分析了网格尺度变化带来的数值频散问题。Wang等(2001)年利用不规则网格实现了对粘弹性介质的模拟，在他的算法中需要在波数域对过渡区进行插值计算，并且只能适应2倍和3倍的空间网格变化。

在正演模拟中时间算子做为一个全局变量，局部变化非常困难。Falk等(1998)、Tessmer(2000)发展了二阶运动方程的变时间算法，但是对一阶速度-应力方程中时间采样在半程上交错计算的情况并不适用。Thomas(2000年)利用交错位移-应力方程实现了大地地震声波模拟中的P波传播。尽管Thomas实现了局部稳定性条件，但是在不同区域边界部分如何实现局部时间采样步长变化并没有给出详细的处理说明。Hayashi等(2001)用3倍的精细网格对近地表的低速区域实现了不规则网格正演模拟。

发明内容

本发明的目的在于针对地表具有低降速带、地下具有薄层介质、小型化构造和小尺度非均质体的情况，提出了一种适用于复杂地震地质模型，并且具有较高计算效率以及模拟精度的基于时空双变网格的弹性波正演模拟技术，对复杂速度模型采用了更为合理的多尺度网格离散方案，局部采用精细网格剖分，从而避免了对整个模型采用精细网格剖分而导致的采样过密带来的计算耗时长，占用内存多等问题。同时，基于采样定理给出的局部稳定性条件，对精细网格进行局部精细时间步长延拓，从根本上提高了计算效率。

本发明所采用的技术方案是：基于时空双变网格的弹性波正演模拟技术，地震波在地下的传播过程可以用弹性波波动方程进行描述，弹性波一阶速度——应力方程求解地震波在给定速度场中的传播，其特征在于其步骤是：

(1)首先针对已知的地震地质模型，根据地下介质的速度参数和弹性参数确定目标区域和背景区域的稳定性条件，从而计算出不同尺度的空间网格大小并对速度场进行离散化处理；

(2)然后利用泰勒公式和有限差分技术得到不同尺度网格点上变量的空间任意精度差分公式，即空间常网格算子和空间变网格算子的任意偶数阶精度空间差分公式，利用推导出的空间差分公式计算精细剖分区域、粗糙剖分区域和中间过渡带区域的波场值；

(3)同时，在时间层的波场外推过程中，对不同网格剖分区域采用不同的时间采样间隔满足局部稳定性条件，分别得到精细时间层内精细网格点处的波场值以及全局时间层内粗糙网格点处的波场值；

(4)在每一个全局时间层内，将精细网格剖分区域的波场值传递到粗糙网格剖分区域，得到整个速度场在此时刻的空间波场值；

(5)最后利用上述方法在每个时间层进行波场外推，得到单炮的正演模拟记录。

地震波正演模拟中对速度场进行网格剖分，基于可变网格和不规则网格的地震波数值模拟方法对地震地质模型进行离散化，根据地震地质模型的背景速度参数和目标区域的弹性参数，将整体速度场剖分为小尺度网格区域、大尺度网格区域和过渡区域三个主要部分。

对各向同性非均匀介质的弹性波一阶速度-应力方程(方程1)进行离散求解，采用交错网格有限差分法，时间域二阶精度的差分格式对各向同性非均匀介质的弹性波一阶速度-应力方程进行离散求解：

$\mathrm{\ρ}\frac{\∂{\mathrm{\υ}}_x}{\∂t}=\frac{\∂{\mathrm{\τ}}_{\mathrm{xx}}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂z}$

$\mathrm{\ρ}\frac{\∂{\mathrm{\υ}}_z}{\∂t}=\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{\mathrm{zz}}}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{xx}}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂{\mathrm{\υ}}_x}{\∂x}+\mathrm{\λ}\frac{\∂{\mathrm{\υ}}_z}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{zz}}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂{\mathrm{\υ}}_z}{\∂z}+\mathrm{\λ}\frac{\∂{\mathrm{\υ}}_x}{\∂x}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂t}=\mathrm{\μ}(\frac{\∂{\mathrm{\υ}}_x}{\∂z}+\frac{\∂{\mathrm{\υ}}_z}{\∂x})---\left(1\right)$

将运动方程(1)改写为离散形式：

${{\mathrm{\υ}}_x^{n+\frac{1}{2}}}_{i,j}={{\mathrm{\υ}}_x^{n-\frac{1}{2}}}_{i,j}+\mathrm{\Δt\ρ}(D_x{\mathrm{\τ}}_{\mathrm{xx}}+D_z{\mathrm{\τ}}_{\mathrm{xz}}){|}_{i,j}^n$

${{\mathrm{\υ}}_z^{n+\frac{1}{2}}}_{i+\frac{1}{2},j+\frac{1}{2}}={{\mathrm{\υ}}_z^{n-\frac{1}{2}}}_{i+\frac{1}{2},j+\frac{1}{2}}+\mathrm{\Δt\ρ}(D_x{\mathrm{\τ}}_{\mathrm{xz}}+D_z{\mathrm{\τ}}_{\mathrm{zz}}){|}_{i+\frac{1}{2},j+\frac{1}{2}}^n$

${{\mathrm{\τ}}_{\mathrm{xx}}^{n+1}}_{i+\frac{1}{2},j}={\mathrm{\τ}}_{\mathrm{xxi}+\frac{1}{2},j}^n+\mathrm{\Δt}{[(\mathrm{\λ}+2\mathrm{\μ})D_x{\mathrm{\υ}}_x+\mathrm{\λ}{D}_z{\mathrm{\υ}}_z]|}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}$

${{\mathrm{\τ}}_{\mathrm{zz}}^{n+1}}_{i+\frac{1}{2},j}={\mathrm{\τ}}_{\mathrm{zzi}+\frac{1}{2},j}^n+\mathrm{\Δt}{[(\mathrm{\λ}+2\mathrm{\μ})D_z{\mathrm{\υ}}_z+\mathrm{\λ}{D}_x{\mathrm{\υ}}_x]|}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}$

${{\mathrm{\τ}}_{\mathrm{xz}}^{n+1}}_{i,j+\frac{1}{2}}={\mathrm{\τ}}_{\mathrm{xzi},j+\frac{1}{2}}^n+\mathrm{\Δt}{\left[\mathrm{\μ}(D_z{\mathrm{\υ}}_x+D_x{\mathrm{\υ}}_z)\right]|}_{i,j+\frac{1}{2}}^{n+\frac{1}{n}}---\left(2\right)$

式中，D_x，D_z分别表示对x，z的一阶微分算子。用变量g来表示速度υ_x、υ_z或应力张量τ_xx、τ_zz、τ_xz，其中，传统的常步长任意偶数阶精度空间差分格式为(Igel，1992)：

$D_{x}g(x,z)=\frac{1}{\mathrm{\Δx}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_i^n\left\{g\right[x+\mathrm{\Δ}{x}_{2i-1},z]-g[x-\mathrm{\Δ}{x}_{2i-1},z\left]\right\}---\left(3\right)$

任意2L阶精度中心有限差分系数由在x处将g[x+Δx_2i-1，z]和g[x-Δx_2i-1，z]作Taylor展开而得到，求解差分系数的矩阵如下：

$\left(\begin{array}{cccc}1& 3& \·\·\·& 2n-1\\ 1& 3^3& \·\·\·& {(2n-1)}^3\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ 1& 3^{2n-1}& \·\·\·& {(2n-1)}^{2n-1}\end{array}\right)\left(\begin{array}{c}{C}_1\\ C_2\\ \·\\ \·\\ \·\\ C_n\end{array}\right)=\left(\begin{array}{c}1\\ 0\\ \·\\ \·\\ \·\\ 0\end{array}\right)---\left(4\right)$

变步长交错网格算法的任意偶数阶精度差分近似式为：

$D_{x}g(x,z)=\underset{i-1}{\overset{n}{\mathrm{\Σ}}}[c_{2i-1}g(x+{\mathrm{\Δ}}_{2i-1},z)+c_{2i}g(x-{\mathrm{\Δ}}_{2i},z)]---\left(5\right)$

其中，c_i为待求的差分系数，Δ_i是空间差分算子，它是网格步长dx的函数，在交错网格技术中差分算子有两个对称点：x_i或x_i+1/2，以对称点x_i+1/2为例求取Δ_i，

Δ₁＝dx_i/2         Δ₂＝dx_i/2

${\mathrm{\Δ}}_{n-1}=\underset{k=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{dx}}_{i+k}+{\mathrm{dx}}_i/2$ ${\mathrm{\Δ}}_n=\underset{k=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{dx}}_{i-k}+{\mathrm{dx}}_i/2$

以x_i为对称点的差分算子为

令g(x，z)＝gz exp(ikx)，(5)式可改写为：

$\mathrm{ik}=\underset{i-1}{\overset{n}{\mathrm{\Σ}}}[c_{2i-1}\exp (+\mathrm{ik}{\mathrm{\Δ}}_{2i-1})+c_{2i}\exp (-\mathrm{ik}{\mathrm{\Δ}}_{2i})]---\left(6\right)$

对(6)式中的指数项进行泰勒展开，它的2n阶泰勒展开式为：

$\exp (+\mathrm{ik}{\mathrm{\Δ}}_i)\≈1+\mathrm{ik}{\mathrm{\Δ}}_i+\frac{1}{2}{\left(\mathrm{ik}\right)}^2{\mathrm{\Δ}}_i^2+......+\frac{{\left(i\right)}^{2n-1}}{(2n-1)!}{\left(k{\mathrm{\Δ}}_i\right)}^{2n-1}+o\left({\mathrm{\Δ}}_i^{2n}\right)---\left(7\right)$

把式(7)代入(6)式并写成矩阵的形式，差分系数由以下方程确定：

通过解方程组(8)，得到有限差分算子D_x的系数c_i。

空间网格实现步长变化的同时，时间采样满足精细空间网格对应的稳定性条件，根据速度场为二维，给出假设：(1)区域B的空间网格步长ΔH为区域A的空间网格步长Δh的2n+1(n＝0，1，2，……)倍，取值为5倍，则对应的区域B的时间采样步长ΔT为区域A的时间采样步长Δt的5倍；(2)下标k代表时间采样间隔为5Δt，区域B的速度分量初始时刻为t_k-1/2，应力分量初始时刻为t_k，区域A的速度分量初始时刻为t_k+3/10，应力分量初始时刻为t_k+4/10；(3)空间差分为四阶精度，且已知t_i-1/2时刻V_x，V_z的波场值和t_i时刻τ_xx，τ_zz，τ_xz的波场值，已知t＝t_k+3/10时刻v_x，v_z的波场值和t＝t_k+4/10时刻t_xx，t_zz，t_xz的波场值；其具体实施步骤和技术路线如下：

(1)首先计算t_i+1/2时刻区域

的V_x，V_z值；计算t＝t_k+5/10时刻区域

时v_x，v_z波场值；然后利用交错网格变步长的空间差分公式(5)计算t＝t_k+5/10时刻过渡区域z_j-1≤z≤z_j+1内的v_x，v_z；

(2)将过渡区域z_j-1≤z≤z_j+1的v_x，v_z赋值给区域B对应点的V_x，V_z。例如：

(3)用和第一步相同的方法求解t_i+1时刻z≥z_j的τ_xx，τ_zz，τ_xz值和t_i+3/2时刻z≥z_j+2时V_x，V_z值，然后根据已知的V_x，V_z值求解v_x，v_z值；

(4)利用双线性插值公式(9)求解精细时间采样平面(t＝t_k+6/10，t_k+7/10，t_k+8/10，t_k+9/10的边界值

$f\left(i\right)=\frac{(\mathrm{nk}-i)F_0+i{F}_1}{\mathrm{nk}}$ i＝0，1，2，.....nk     (9)

其中F₀，F₁表示待求点f(i)的前一全局时刻和下一全局时刻的波场值，相同边界处的应力值t_xx，t_zz，t_xz用同样的方法利用相邻两个全局时刻的τ_xx，τ_zz，τ_xz值计算得到；

(5)利用公式(3)计算精细采样平面内部z ≤z_j-2的的v_x，v_z值和t_xx，t_zz，t_xz值和t＝t_k+10/10时刻z≤z_j-1的t_xx，t_zz，t_xz值，并将过渡区域z_j-1≤z≤z_j+1的t_xx，t_zz，t_xz值赋予区域B对应点的τ_xx，τ_zz，τ_xz。

与现有技术相比，本发明的有益效果是：满足了浅层的局部稳定性条件，从而避免了频散噪音的产生，比传统算法节省了计算时间，提高了计算效率，避免了空间网格变化可能会导致的计算不稳定和人工数值反射现象，实现了精细插值算法处理，提高了算法的精度和可靠性，且模拟过程较为稳定。

附图说明

图1为速度场离散方式(1)。

图2为速度场离散方式(2)。

图3为空间网格变化图。

图4为局部时间变化图。

图5为中浅层速度模型。

图6为基于传统算法的x分量模拟记录。

图7为基于时空变网格算法的x分量模拟记录。

图8为陡坡带砂砾岩体速度模型。

图9为模拟耗时对比图。

图10为时空双变网格算法流程图。

图11为裂缝发育带位于灰岩层。

图12为含裂缝模型叠加剖面。

图13为无裂缝模型叠加剖面。

图14为炮点位于轴部区域时的叠前单炮记录局部放大图，其中(a)含裂缝模型单炮记录(b)无裂缝模型单炮记录。

具体实施方式

如图10所示：基于时空双变网格的弹性波正演模拟技术，地震波在地下的传播过程可以用弹性波波动方程进行描述，弹性波一阶速度——应力方程求解地震波在给定速度场中的传播，其步骤是：

(1)首先针对已知的地震地质模型，根据地下介质的速度参数和弹性参数确定目标区域和背景区域的稳定性条件，从而计算出不同尺度的空间网格大小并对速度场进行离散化处理；

(2)然后利用泰勒公式和有限差分技术得到不同尺度网格点上变量的空间任意精度差分公式，即空间常网格算子和空间变网格算子的任意偶数阶精度空间差分公式，利用推导出的空间差分公式计算精细剖分区域、粗糙剖分区域和中间过渡带区域的波场值；

(3)同时，在时间层的波场外推过程中，对不同网格剖分区域采用不同的时间采样间隔满足局部稳定性条件，分别得到精细时间层内精细网格点处的波场值以及全局时间层内粗糙网格点处的波场值；

(4)在每一个全局时间层内，将精细网格剖分区域的波场值传递到粗糙网格剖分区域，得到整个速度场在此时刻的空间波场值；

(5)最后利用上述方法在每个时间层进行波场外推，得到单炮的正演模拟记录。

地震波正演模拟中对速度场进行网格剖分，基于可变网格和不规则网格的地震波数值模拟方法对地震地质模型进行离散化，根据地震地质模型的背景速度参数和目标区域的弹性参数，将整体速度场剖分为小尺度网格区域、大尺度网格区域和过渡区域三个主要部分。

若对速度场离散的不同尺度网格均接近于地震勘探的尺度范围(米级)，则不需要过渡区域，两种离散方式的模型剖面如图(1，2)所示。

图1为速度模型不规则剖分时，局部储层尺度较大时，不需要过渡区域的情况；而当我们需要精细描述储层特征时，比如储层描述尺度为毫米级时，在剖分速度模型时，需要加入速度过渡区域，从而实现速度模型的自然过渡，波场正演过程也较为稳定。

采用交错网格有限差分法，时间域二阶精度的差分格式对各向同性非均匀介质的弹性波一阶速度-应力方程进行离散求解：

$\mathrm{\ρ}\frac{\∂{\mathrm{\υ}}_x}{\∂t}=\frac{\∂{\mathrm{\τ}}_{\mathrm{xx}}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂z}$

$\mathrm{\ρ}\frac{\∂{\mathrm{\υ}}_z}{\∂t}=\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{\mathrm{zz}}}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{xx}}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂{\mathrm{\υ}}_x}{\∂x}+\mathrm{\λ}\frac{\∂{\mathrm{\υ}}_z}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{zz}}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂{\mathrm{\υ}}_z}{\∂z}+\mathrm{\λ}\frac{\∂{\mathrm{\υ}}_x}{\∂x}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂t}=\mathrm{\μ}(\frac{\∂{\mathrm{\υ}}_x}{\∂z}+\frac{\∂{\mathrm{\υ}}_z}{\∂x})---\left(1\right)$

将运动方程(1)改写为离散形式：

${{\mathrm{\υ}}_x^{n+\frac{1}{2}}}_{i,j}={{\mathrm{\υ}}_x^{n-\frac{1}{2}}}_{i,j}+\mathrm{\Δt\ρ}(D_x{\mathrm{\τ}}_{\mathrm{xx}}+D_z{\mathrm{\τ}}_{\mathrm{xz}}){|}_{i,j}^n$

${{\mathrm{\υ}}_z^{n+\frac{1}{2}}}_{i+\frac{1}{2},j+\frac{1}{2}}={{\mathrm{\υ}}_z^{n-\frac{1}{2}}}_{i+\frac{1}{2},j+\frac{1}{2}}+\mathrm{\Δt\ρ}(D_x{\mathrm{\τ}}_{\mathrm{xz}}+D_z{\mathrm{\τ}}_{\mathrm{zz}}){|}_{i+\frac{1}{2},j+\frac{1}{2}}^n$

${{\mathrm{\τ}}_{\mathrm{xx}}^{n+1}}_{i+\frac{1}{2},j}={\mathrm{\τ}}_{\mathrm{xxi}+\frac{1}{2},j}^n+\mathrm{\Δt}{[(\mathrm{\λ}+2\mathrm{\μ})D_x{\mathrm{\υ}}_x+\mathrm{\λ}{D}_z{\mathrm{\υ}}_z]|}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}$

${{\mathrm{\τ}}_{\mathrm{zz}}^{n+1}}_{i+\frac{1}{2},j}={\mathrm{\τ}}_{\mathrm{zzi}+\frac{1}{2},j}^n+\mathrm{\Δt}{[(\mathrm{\λ}+2\mathrm{\μ})D_z{\mathrm{\υ}}_z+\mathrm{\λ}{D}_x{\mathrm{\υ}}_x]|}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}$

${{\mathrm{\τ}}_{\mathrm{xz}}^{n+1}}_{i,j+\frac{1}{2}}={\mathrm{\τ}}_{\mathrm{xzi},j+\frac{1}{2}}^n+\mathrm{\Δt}{\left[\mathrm{\μ}(D_z{\mathrm{\υ}}_x+D_x{\mathrm{\υ}}_z)\right]|}_{i,j+\frac{1}{2}}^{n+\frac{1}{n}}---\left(2\right)$

式中，D_x，D_z分别表示对x，z的一阶微分算子。用变量g来表示速度υ_x、υ_z或应力张量τ_xx、τ_zz、τ_xz，其中，传统的常步长任意偶数阶精度空间差分格式为(Igel，1992)：

$D_{x}g(x,z)=\frac{1}{\mathrm{\Δx}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_i^n\left\{g\right[x+\mathrm{\Δ}{x}_{2i-1},z]-g[x-\mathrm{\Δ}{x}_{2i-1},z\left]\right\}---\left(3\right)$

任意2L阶精度中心有限差分系数

由在x处将g[x+Δx_2i-1，z]和g[x-Δx_2i-1，z]作Taylor展开而得到，求解差分系数的矩阵如下：

$\left(\begin{array}{cccc}1& 3& \·\·\·& 2n-1\\ 1& 3^3& \·\·\·& {(2n-1)}^3\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ 1& 3^{2n-1}& \·\·\·& {(2n-1)}^{2n-1}\end{array}\right)\left(\begin{array}{c}{C}_1\\ C_2\\ \·\\ \·\\ \·\\ C_n\end{array}\right)=\left(\begin{array}{c}1\\ 0\\ \·\\ \·\\ \·\\ 0\end{array}\right)---\left(4\right)$

变步长交错网格算法的任意偶数阶精度差分近似式为：

$D_{x}g(x,z)=\underset{i-1}{\overset{n}{\mathrm{\Σ}}}[c_{2i-1}g(x+{\mathrm{\Δ}}_{2i-1},z)+c_{2i}g(x-{\mathrm{\Δ}}_{2i},z)]---\left(5\right)$

其中，c_i为待求的差分系数，Δ_i是空间差分算子，它是网格步长dx的函数，在交错网格技术中差分算子有两个对称点：x_i或x_i+1/2，以对称点x_i+1/2为例求取Δ_i，

Δ₁＝dx_i/2       Δ₂＝dx_i/2

${\mathrm{\Δ}}_{n-1}=\underset{k=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{dx}}_{i+k}+{\mathrm{dx}}_i/2$ ${\mathrm{\Δ}}_n=\underset{k=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{dx}}_{i-k}+{\mathrm{dx}}_i/2$

以xi为对称点的差分算子为

令g(x，z)＝g_z exp(ikx)，(5)式可改写为：

$\mathrm{ik}=\underset{i-1}{\overset{n}{\mathrm{\Σ}}}[c_{2i-1}\exp (+\mathrm{ik}{\mathrm{\Δ}}_{2i-1})+c_{2i}\exp (-\mathrm{ik}{\mathrm{\Δ}}_{2i})]---\left(6\right)$

对(6)式中的指数项进行泰勒展开，它的2n阶泰勒展开式为：

$\exp (+\mathrm{ik}{\mathrm{\Δ}}_i)\≈1+\mathrm{ik}{\mathrm{\Δ}}_i+\frac{1}{2}{\left(\mathrm{ik}\right)}^2{\mathrm{\Δ}}_i^2+......+\frac{{\left(i\right)}^{2n-1}}{(2n-1)!}{\left(k{\mathrm{\Δ}}_i\right)}^{2n-1}+o\left({\mathrm{\Δ}}_i^{2n}\right)---\left(7\right)$

把式(7)代入(6)式并写成矩阵的形式，差分系数由以下方程确定：

通过解方程组(8)，得到有限差分算子D_x的系数c_i。

在空间网格实现步长变化的同时，由于稳定性条件的限制，时间采样需要满足精细空间网格对应的稳定性条件，因此局部稳定性条件的实施使得时间采样变得宽松，才使得计算效率提高从根本上得以解决。图3和图4分别为不同网格步长和不同时间步长的表征方式。

空间网格实现步长变化的同时，时间采样满足精细空间网格对应的稳定性条件，因为速度场为二维，给出假设：(1)区域B的空间网格步长ΔH为区域A的空间网格步长Δh的2n+1(n＝0，1，2，……)倍，取值为5倍，则对应的区域B的时间采样步长ΔT为区域A的时间采样步长Δt的5倍；(2)下标k代表时间采样间隔为5Δt，区域B的速度分量初始时刻为t_k-1/2，应力分量初始时刻为t_k，区域A的速度分量初始时刻为t_k+3/10，应力分量初始时刻为t_k+4/10；(3)空间差分为四阶精度，且已知t_i-1/2时刻V_x，V_z的波场值和t_i时刻τ_xx，τ_zz，τ_xz的波场值，已知t＝t_k+3/10时刻v_x，v_z的波场值和t＝t_k+4/10时刻t_xx，t_zz，t_xz的波场值；其具体实施步骤和技术路线如下：

(1)首先计算t_i+1/2时刻区域

的V_x，V_z值；计算t＝t_k+5/10时刻区域

时v_x，v_z波场值；然后利用交错网格变步长的空间差分公式(5)计算t＝t_k+5/10时刻过渡区域z_j-1≤z≤z_j+1内的v_x，v_z；

(2)将过渡区域z_j-1≤z≤z_j+1的v_x，v_z赋值给区域B对应点的V_x，V_z。例如：

(3)用和第一步相同的方法求解t_i+1时刻z≥z_j的τ_xx，τ_zz，τ_xz值和t_i+3/2时刻z≥z_j+2时V_x，V_z值，然后根据已知的V_x，V_z值求解v_x，v_z值；

(4)利用双线性插值公式(9)求解精细时间采样平面(t＝t_k+6/10，t_k+7/10，t_k+8/10，t_k+9/10的边界值

$f\left(i\right)=\frac{(\mathrm{nk}-i)F_0+i{F}_1}{\mathrm{nk}}$ i＝0，1，2，.....nk    (9)

其中F₀，F₁表示待求点f(i)的前一全局时刻和下一全局时刻的波场值，相同边界处的应力值t_xx，t_zz，t_xz用同样的方法利用相邻两个全局时刻的τ_xx，τ_zz，τ_xz值计算得到；

(5)利用公式(3)计算精细采样平面内部z≤z_j-2的的v_x，v_z值和t_xx，t_zz，t_xz值和t＝t_k+10/10时刻z ≤z_j-1的t_xx，t_zz，t_xz值，并将过渡区域z_j-1≤z≤z_j+1的t_xx，t_zz，t_xz值赋予区域B对应点的τ_xx，τ_zz，τ_xz。

为了展示本专利技术的发明效果，我们给出了两个模型进行说明：(1)表层低降速带模型算例，(2)陡坡带砂砾岩体速度模型算例；

算例(1)：利用时空双变网格正演模拟技术对浅层含有低速带的模型进行了试算，并与传统交错网格算法进行了对比，验证其计算精度和计算效率。图5为浅层速度模型，图6为基于传统算法的x分量模拟记录，图7为基于时空双变网格算法的x分量模拟记录。基于传统算法的模拟记录中出现了大量的频散噪音(如图7中箭头所示区域)，湮灭了有效反射波能量。而采用时空双变网格正演算法利用精细网格满足了浅层的局部稳定性条件，从而避免了频散噪音的产生(如图8中箭头所示区域)。

(2)陡坡带砂砾岩体速度模型算例

图8是某地区的陡坡带砂砾岩体速度模型。速度模型中部层位较多，层厚度较薄。因此对模型中部进行精细网格剖分，并采用局部变时间算法进行正演模拟(矩形方框区域)。模型参数和模拟参数如表1。在大型集群Cluster下对两种算法的耗时进行了对比，硬件环境为CPU Speed 2500hz，Memory Total 8166004KB，Swap Space Total 16579072KB。图9为模拟耗时对比图，其中采用时空精细网格剖分的传统算法、采用时间精细剖分的空间变网格算法和采用双变算法得到的单炮模拟记录精度一致，但是采用传统算法(全局采用精细剖分)单炮耗时为16613s，采用空间变网格剖分算法单炮耗时为8326s，采用双变算法单炮耗时为3150s，相比传统算法节省了81.3％的计算时间，相比仅采用空间变网格算法节省了62.6％的计算时间。

表1模型参数和正演模拟参数

本发明在在不影响计算精度的基础上，相比传统算法可观的提高了计算效率，可以极大节约正演计算时间。对大规模实际速度模型的处理和海量速度模型的正演模拟研究中，更能够体现本专利技术的优势。

基于上述方法原理推导及技术路线，根据某地区数据记录及裂缝发育模式，对裂缝发育带位于灰岩层的实际含裂缝速度模型进行了正演模拟。模型参数如下：横测线长度为4500m，深度为5500m。裂缝主要发育区域的岩性为灰岩和白云岩。图10为灰岩中含裂缝模型，裂缝紧附于灰岩顶界面，宽度为2mm、长度为20m、裂缝密度为0.1、油充填。图13为含裂缝模型的叠加剖面，图14为无裂缝模型的叠加剖面。红色方框标出部分为2.5s-3.0s的模拟记录部分，裂缝的反射信息多集中在此区域。

对比计算结果：采用规则交错网格技术对裂缝介质描述需要2.475*10⁸个空间网格点，同时根据稳定性条件，时间层采样为0.1μs，也就是说即使介质在空间可以描述，计算机也无法承受如此巨大的时间计算量和低下的计算效率；而利用空间变网格技术则只需要9.394*10⁵个空间网格点，网格数量相比减少了62.04％，从而极大的减小了模型的计算量。并且在时间层延拓中全局仍采用粗糙网格控制的时间步长(0.7ms)，只在局部进行时间层的细分，使一个全局时间层内只有较少的空间网格需要几千次的精细迭代，从而使得对碳酸盐岩裂缝性储层速度模型的正演模拟成为可能。

在叠加剖面中可以观测到微弱的裂缝反射信息，沿地层走向(倾末端)表现为不连续的块状反射轴，构造起伏部位(左翼、轴部)表现为纵向上相邻的多条弱绕射波，受裂缝群影响，裂缝下伏层反射波横向上连续性变差。

在单炮记录中，(1)非连续“段状”反射多表现为平缓构造影响下的裂缝反射信息(如圆圈1所示，为倾末端区域的裂缝反射信息)；(2)起伏构造和埋藏较深的白云岩中，多表现为杂乱反射，或横向极不连续的小块状反射(如中间圆圈所示，为轴部和翼部构造区域的裂缝反射信息)。

我国西部探区碳酸盐岩裂缝型储层是时空双变网格技术的应用领域之一，本技术在此应用中主要体现了两个创新点：

(1)裂缝型储层的正演模拟技术一直是个难点，这是因为地下实际介质中的裂缝开度一般较小(多为毫米量级甚至微米量级，无法用正常尺度的空间网格进行剖分，而精细化网格带来的计算量和计算时间是目前的计算机无法承受的)，埋藏较深(只能在大型速度场中进行模拟)，具有较强的非均质性(传统的正演算法存在技术上的局限性)。本技术利用空间平滑渐进式的剖分网格对速度模型中微裂缝进行了表征，解决了传统正演算法无法表征裂缝的众多弹性参数的缺点，同时采用了局部稳定性条件，极大的提高了计算效率。本专利技术使得直接求解方程得到裂缝型储层的地震波传播特征的方法，从不可能变为可能。

(2)同时，与裂缝储层各向异性介质的等效模拟相比，本专利技术通过直接求解弹性波波动方程，能更直观的体现裂缝的波场特征(如导波、尾波和首波等)，不再是一个平均和模糊处理后的裂缝特征表现。

除了碳酸盐岩裂缝型储层之外，本专利技术还可以广泛应用于孔洞介质、油藏薄层模型、低降速带速度模型和强纵横向速度变化的速度模型。由于采用了变时间变网格的计算策略，本专利技术相对于传统的固定网格正演模拟方法，具有更强的速度模型适应性，同时具有更高的计算效率，使得毫米级尺度的裂缝储层的弹性波正演模拟成为可能，极大提高了弹性波波动方程正演模拟的模型实用范围，使得模拟技术更趋于精细化、实用化。

Claims (4)

1.基于时空双变网格的弹性波正演模拟技术，地震波在地下的传播过程可以用弹性波波动方程进行描述，弹性波一阶速度——应力方程求解地震波在给定速度场中的传播，其特征在于其步骤是：

(1)首先针对已知的地震地质模型，根据地下介质的速度参数和弹性参数确定目标区域和背景区域的稳定性条件，从而计算出不同尺度的空间网格大小并对速度场进行离散化处理；

(2)然后利用泰勒公式和有限差分技术得到不同尺度网格点上变量的空间任意精度差分公式，即空间常网格算子和空间变网格算子的任意偶数阶精度空间差分公式，利用推导出的空间差分公式计算精细剖分区域、粗糙剖分区域和中间过渡带区域的波场值；

(3)同时，在时间层的波场外推过程中，对不同网格剖分区域采用不同的时间采样间隔满足局部稳定性条件，分别得到精细时间层内精细网格点处的波场值以及全局时间层内粗糙网格点处的波场值；

(4)在每一个全局时间层内，将精细网格剖分区域的波场值传递到粗糙网格剖分区域，得到整个速度场在此时刻的空间波场值；

(5)最后利用上述方法在每个时间层进行波场外推，得到单炮的正演模拟记录。

2.根据权利要求1所述的基于时空双变网格的弹性波正演模拟技术，其特征在于地震波正演模拟中对速度场进行网格剖分，基于可变网格和不规则网格的地震波数值模拟方法对地震地质模型进行离散化，根据地震地质模型的背景速度参数和目标区域的弹性参数，将整体速度场剖分为小尺度网格区域、大尺度网格区域和过渡区域三个主要部分。

3.根据权利要求1所述的基于时空双变网格的弹性波正演模拟技术，其特征在于对各向同性非均匀介质的弹性波一阶速度-应力方程(方程1)进行离散求解，采用交错网格有限差分法，时间域二阶精度的差分格式对各向同性非均匀介质的弹性波一阶速度-应力方程进行离散求解：

$\mathrm{\ρ}\frac{\∂{\mathrm{\υ}}_x}{\∂t}=\frac{\∂{\mathrm{\τ}}_{\mathrm{xx}}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂z}$

$\mathrm{\ρ}\frac{\∂{\mathrm{\υ}}_z}{\∂t}=\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{\mathrm{zz}}}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{xx}}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂{\mathrm{\υ}}_x}{\∂x}+\mathrm{\λ}\frac{\∂{\mathrm{\υ}}_z}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{zz}}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂{\mathrm{\υ}}_z}{\∂z}+\mathrm{\λ}\frac{\∂{\mathrm{\υ}}_x}{\∂x}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂t}=\mathrm{\μ}(\frac{\∂{\mathrm{\υ}}_x}{\∂z}+\frac{\∂{\mathrm{\υ}}_z}{\∂x})---\left(1\right)$

将运动方程(1)改写为离散形式：

${{\mathrm{\υ}}_x^{n+\frac{1}{2}}}_{i,j}={{\mathrm{\υ}}_x^{n-\frac{1}{2}}}_{i,j}+\mathrm{\Δt\ρ}(D_x{\mathrm{\τ}}_{\mathrm{xx}}+D_z{\mathrm{\τ}}_{\mathrm{xz}}){|}_{i,j}^n$

${{\mathrm{\υ}}_z^{n+\frac{1}{2}}}_{i+\frac{1}{2},j+\frac{1}{2}}={{\mathrm{\υ}}_z^{n-\frac{1}{2}}}_{i+\frac{1}{2},j+\frac{1}{2}}+\mathrm{\Δt\ρ}(D_x{\mathrm{\τ}}_{\mathrm{xz}}+D_z{\mathrm{\τ}}_{\mathrm{zz}}){|}_{i+\frac{1}{2},j+\frac{1}{2}}^n$

${{\mathrm{\τ}}_{\mathrm{xx}}^{n+1}}_{i+\frac{1}{2},j}={\mathrm{\τ}}_{\mathrm{xxi}+\frac{1}{2},j}^n+\mathrm{\Δt}{[(\mathrm{\λ}+2\mathrm{\μ})D_x{\mathrm{\υ}}_x+\mathrm{\λ}{D}_z{\mathrm{\υ}}_z]|}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}$

${{\mathrm{\τ}}_{\mathrm{zz}}^{n+1}}_{i+\frac{1}{2},j}={\mathrm{\τ}}_{\mathrm{zzi}+\frac{1}{2},j}^n+\mathrm{\Δt}{[(\mathrm{\λ}+2\mathrm{\μ})D_z{\mathrm{\υ}}_z+\mathrm{\λ}{D}_x{\mathrm{\υ}}_x]|}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}$

${{\mathrm{\τ}}_{\mathrm{xz}}^{n+1}}_{i,j+\frac{1}{2}}={\mathrm{\τ}}_{\mathrm{xzi},j+\frac{1}{2}}^n+\mathrm{\Δt}{\left[\mathrm{\μ}(D_z{\mathrm{\υ}}_x+D_x{\mathrm{\υ}}_z)\right]|}_{i,j+\frac{1}{2}}^{n+\frac{1}{n}}---\left(2\right)$

式中，D_x，D_z分别表示对x，z的一阶微分算子。用变量g来表示速度υ_x、υ_z或应力张量τ_xx、τ_zz、τ_xz，其中，传统的常步长任意偶数阶精度空间差分格式为(Igel，1992)：

$D_{x}g(x,z)=\frac{1}{\mathrm{\Δx}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_i^n\left\{g\right[x+\mathrm{\Δ}{x}_{2i-1},z]-g[x-\mathrm{\Δ}{x}_{2i-1},z\left]\right\}---\left(3\right)$

任意2L阶精度中心有限差分系数由在x处将g[x+Δx_2i-1，z]和g[x-Δx_2i-1，z]作Taylor展开而得到，求解差分系数的矩阵如下：

$\left(\begin{array}{cccc}1& 3& \·\·\·& 2n-1\\ 1& 3^3& \·\·\·& {(2n-1)}^3\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ 1& 3^{2n-1}& \·\·\·& {(2n-1)}^{2n-1}\end{array}\right)\left(\begin{array}{c}{C}_1\\ C_2\\ \·\\ \·\\ \·\\ C_n\end{array}\right)=\left(\begin{array}{c}1\\ 0\\ \·\\ \·\\ \·\\ 0\end{array}\right)---\left(4\right)$

变步长交错网格算法的任意偶数阶精度差分近似式为：

$D_{x}g(x,z)=\underset{i-1}{\overset{n}{\mathrm{\Σ}}}[c_{2i-1}g(x+{\mathrm{\Δ}}_{2i-1},z)+c_{2i}g(x-{\mathrm{\Δ}}_{2i},z)]---\left(5\right)$

其中，c_i为待求的差分系数，Δ_i是空间差分算子，它是网格步长dx的函数，在交错网格技术中差分算子有两个对称点：x_i或x_i+1/2，以对称点x_i+1/2为例求取Δ_i，

Δ₁＝dx_i/2        Δ₂＝dx_i/2

${\mathrm{\Δ}}_{n-1}=\underset{k=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{dx}}_{i+k}+{\mathrm{dx}}_i/2$ ${\mathrm{\Δ}}_n=\underset{k=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{dx}}_{i-k}+{\mathrm{dx}}_i/2$

以x_i为对称点的差分算子为

令g(x，z)＝g_z exp(ikx)，(5)式可改写为：

$\mathrm{ik}=\underset{i-1}{\overset{n}{\mathrm{\Σ}}}[c_{2i-1}\exp (+\mathrm{ik}{\mathrm{\Δ}}_{2i-1})+c_{2i}\exp (-\mathrm{ik}{\mathrm{\Δ}}_{2i})]---\left(6\right)$

对(6)式中的指数项进行泰勒展开，它的2n阶泰勒展开式为：

$\exp (+\mathrm{ik}{\mathrm{\Δ}}_i)\≈1+\mathrm{ik}{\mathrm{\Δ}}_i+\frac{1}{2}{\left(\mathrm{ik}\right)}^2{\mathrm{\Δ}}_i^2+......+\frac{{\left(i\right)}^{2n-1}}{(2n-1)!}{\left(k{\mathrm{\Δ}}_i\right)}^{2n-1}+o\left({\mathrm{\Δ}}_i^{2n}\right)---\left(7\right)$

把式(7)代入(6)式并写成矩阵的形式，差分系数由以下方程确定：

通过解方程组(8)，得到有限差分算子D_x的系数c_i。

4.根据权利要求1所述的基于时空双变网格的弹性波正演模拟技术，其特征在于空间网格实现步长变化的同时，时间采样满足精细空间网格对应的稳定性条件，根据速度场为二维，给出假设：(1)区域B的空间网格步长ΔH为区域A的空间网格步长Δh的2n+1(n＝0，1，2，……)倍，取值为5倍，则对应的区域B的时间采样步长ΔT为区域A的时间采样步长Δt的5倍；(2)下标k代表时间采样间隔为5Δt，区域B的速度分量初始时刻为t_k-1/2，应力分量初始时刻为t_k，区域A的速度分量初始时刻为t_k+3/10，应力分量初始时刻为t_k+4/10；(3)空间差分为四阶精度，且已知t_i-1/2时刻V_x，V_z的波场值和t_i时刻τ_xx，τ_zz，τ_xz的波场值，已知t＝t_k+3/10时刻v_x，v_z的波场值和t＝t_k+4/10时刻t_xx，t_zz，t_xz的波场值；其具体实施步骤和技术路线如下：

(1)首先计算t_i+1/2时刻区域

的V_x，V_z值；计算t＝t_k+5/10时刻区域

时v_x，v_z波场值；然后利用交错网格变步长的空间差分公式(5)计算t＝t_k+5/10时刻过渡区域z_j-1≤z ≤z_j+1内的v_x，v_z；

(2)将过渡区域z_j-1≤z≤z_j+1的v_x，v_z赋值给区域B对应点的V_x，V_z。例如：

(3)用和第一步相同的方法求解t_i+1时刻z≥z_j的τ_xx，τ_zz，τ_xz值和t_i+3/2时刻z≥z_j+2时V_x，V_z值，然后根据已知的V_x，V_z值求解v_x，v_z值；

(4)利用双线性插值公式(9)求解精细时间采样平面(t＝t_k+6/10，t_k+7/10，t_k+8/10，t_k+9/10的边界值

$f\left(i\right)=\frac{(\mathrm{nk}-i)F_0+i{F}_1}{\mathrm{nk}}$ i＝0，1，2，.....nk    (9)

其中F₀，F₁表示待求点f(i)的前一全局时刻和下一全局时刻的波场值，相同边界处的应力值t_xx，t_zz，t_xz用同样的方法利用相邻两个全局时刻的τ_xx，τ_zz，τ_xz值计算得到；

(5)利用公式(3)计算精细采样平面内部z≤z_j-2的的v_x，v_z值和t_xx，t_zz，t_xz值和t＝t_k+10/10时刻z≤z_j-1的t_xx，t_zz，t_xz值，并将过渡区域z_j-1≤z≤z_j+1的t_xx，t_zz，t_xz值赋予区域B对应点的τ_xx，τ_zz，_xz。

Images