[发明专利]基于时空双变网格的弹性波正演模拟技术

权利要求书

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

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

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

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

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

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

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

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

ρ∂υx∂t=∂τxx∂x+∂τxz∂z]]>

ρ∂υz∂t=∂τxz∂x+∂τzz∂z]]>

∂τxx∂t=(λ+2μ)∂υx∂x+λ∂υz∂z]]>

∂τzz∂t=(λ+2μ)∂υz∂z+λ∂υx∂x]]>

∂τxz∂t=μ(∂υx∂z+∂υz∂x)---(1)]]>

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

υxn+12i,j=υxn-12i,j+Δtρ(Dxτxx+Dzτxz)|i,jn]]>

υzn+12i+12,j+12=υzn-12i+12,j+12+Δtρ(Dxτxz+Dzτzz)|i+12,j+12n]]>

τxxn+1i+12,j=τxxi+12,jn+Δt[(λ+2μ)Dxυx+λDzυz]|i+12,jn+12]]>

τzzn+1i+12,j=τzzi+12,jn+Δt[(λ+2μ)Dzυz+λDxυx]|i+12,jn+12]]>

τxzn+1i,j+12=τxzi,j+12n+Δt[μ(Dzυx+Dxυz)]|i,j+12n+1n---(2)]]>

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

Dxg(x,z)=1ΔxΣi=1nCin{g[x+Δx2i-1,z]-g[x-Δx2i-1,z]}---(3)]]>

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

13···2n-1133···(2n-1)3············132n-1···(2n-1)2n-1C1C2···Cn=10···0---(4)]]>

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

Dxg(x,z)=Σi-1n[c2i-1g(x+Δ2i-1,z)+c2ig(x-Δ2i,z)]---(5)]]>

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

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

Δn-1=Σk=1n-1dxi+k+dxi/2]]>Δn=Σk=1n-1dxi-k+dxi/2]]>

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

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

ik=Σi-1n[c2i-1exp(+ikΔ2i-1)+c2iexp(-ikΔ2i)]---(6)]]>

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

exp(+ikΔi)≈1+ikΔi+12(ik)2Δi2+......+(i)2n-1(2n-1)!(kΔi)2n-1+o(Δi2n)---(7)]]>

把式(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(i)=(nk-i)F0+iF1nk]]>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。