基于交错网格有限差分算法的波数域保幅波场分离方法

摘要

本发明涉及一种基于交错网格有限差分算法的波数域保幅波场分离方法，包括以下步骤：利用交错网格有限差分算法进行弹性波场数值模拟；利用傅里叶变换将波场快照从空间域变换到波数域；利用波数域插值算子将速度分量插值到相同的网格节点；利用归一化波数进行保幅波场分离；利用傅里叶反变换将分离结果变换到空间域获得保幅的P波和S波。本发明利用交错网格有限差分算法可实现较高精度的地震波场数值模拟；波数域的插值算子具有很高的插值精度，因此可以很好的估算相同网格节点的参数值；利用归一化波数可实现保幅波场分离，获得保幅的P波和S波波场；由于波场分离是在波数域进行的，相对于正反傅里叶变换，波数域插值并不会显著增加计算量。

主权项

一种基于交错网格有限差分算法的波数域保幅波场分离方法，其特征在于，包括以下步骤：步骤1：利用交错网格有限差分算法进行弹性波场数值模拟；其中，所述步骤1为：利用交错网格有限差分格式对弹性波一阶速度‑应力方程进行离散，获得弹性波传播算子；加载给定震源子波，利用弹性波传播算子进行弹性波场数值模拟；具体方法如下：利用弹性动力学的三个基本方程：描述应力应变关系的本构方程、描述应力位移关系的运动平衡微分方程以及描述位移应变关系的几何方程得到弹性波一阶速度‑应力方程：$\frac{\partial V(x,z,t)}{\partial t}=(A^x\frac{\partial }{\partial x}+A^z\frac{\partial }{\partial z})V(x,z,t),$其中，$V=\left[\begin{array}{c}{v}_x\\ v_z\\ {\tau }_{xx}\\ {\tau }_{zz}\\ {\tau }_{xz}\end{array}\right],A^x=\left[\begin{array}{ccccc}0& 0& \frac{1}{\rho }& 0& 0\\ 0& 0& 0& 0& \frac{1}{\rho }\\ C_{11}& 0& 0& 0& 0\\ C_{13}& 0& 0& 0& 0\\ 0& C_{44}& 0& 0& 0\end{array}\right],A^z=\left[\begin{array}{ccccc}0& 0& 0& 0& \frac{1}{\rho }\\ 0& 0& 0& \frac{1}{\rho }& 0\\ 0& C_{13}& 0& 0& 0\\ 0& C_{33}& 0& 0& 0\\ C_{44}& 0& 0& 0& 0\end{array}\right]$v_x、v_z分别为质点振动速度的水平分量和垂直分量，τ_ij为应力分量，C_ij为弹性常数，ρ为介质密度，t表示时间；利用交错网格有限差分格式对弹性波一阶速度‑应力方程进行离散，得到弹性波传播算子：$V_{i,j}^{n+1}=V_{i,j}^n+{\mathrm{\Delta tA}}^x{L}_x\left(V_{ai,j}^{n+\frac{1}{2}}\right)+{\mathrm{\Delta tA}}^z{L}_z\left(V_{bi,j}^{n+\frac{1}{2}}\right),$其中，$\left\{\begin{array}{c}{V}_{ai,j}^n={\left[v_{xi+\frac{1}{2},j}^n,v_{zi,j+\frac{1}{2}}^n,{\tau }_{xxi,j}^n,{\tau }_{zzi,j}^n,{\tau }_{xzi+\frac{1}{2},j+\frac{1}{2}}^n\right]}^T\\ V_{bi,j}^n={\left[v_{xi,j+\frac{1}{2}}^n,v_{zi+\frac{1}{2},j}^n,{\tau }_{xxi,j}^n,{\tau }_{zzi+\frac{1}{2},j+\frac{1}{2}}^n,{\tau }_{xzi,j}^n\right]}^T\end{array}\right.$L_x、L_z为空间2L阶有限差分算子；由于交错网格一阶导数的2L阶精度差分可近似表示为：${\mathrm{\Delta xf}}^{\left(1\right)}\left(x_0\right)=\underset{m=1}{\overset{L}{\Sigma }}{a}_m[f(x_0+\frac{2m-1}{2}\Delta x)-f(x_0-\frac{2m-1}{2}\Delta x)],$其中，Δx为网格间距，a_m为高阶交错网格差分系数；因此，L_x、L_z表示如下：$\left\{\begin{array}{c}{L}_x\left(V_{i,j}^n\right)=\frac{1}{\Delta x}\underset{m=1}{\overset{L}{\Sigma }}{a}_m\left(V_{i+\frac{2m-1}{2},j}^n-V_{i-\frac{2m-1}{2},j}^n\right)\\ L_z\left(V_{i,j}^n\right)=\frac{1}{\Delta z}\underset{m=1}{\overset{L}{\Sigma }}{a}_m\left(V_{i,j+\frac{2m-1}{2}}^n-V_{i,j-\frac{2m-1}{2}}^n\right)\end{array}\right..$其中，a_m通过以下方程确定：加载给定的震源子波，利用弹性波传播算子，即可实现弹性波场数值模拟，构建地下弹性矢量波场；步骤2：利用傅里叶变换将波场快照从空间域变换到波数域；所述步骤2为：利用傅里叶变换将步骤1所得的每一时间步长对应的波场快照从空间域变换到波数域；具体方法如下：2D傅里叶正变换为$\tilde{V}(k_x,k_z)=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}v(x,z)e^{-i2\pi (k_{x}x+k_{z}z)}dxdz,$2D傅里叶反变换为$v(x,z)=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\tilde{V}(k_x,k_z)e^{i2\pi (k_{x}x+k_{z}z)}{\mathrm{dk}}_x{\mathrm{dk}}_z,$其中，v为空间域质点振动速度参数，为其对应的波数域结果，x、z分别为空间域水平方向与垂直方向的坐标，k_x、k_z分别为波数域水平方向与垂直方向的波数；利用2D傅里叶正变换将每一时间步长对应的波场快照从空间域变换到波数域；步骤3：利用波数域插值算子将速度分量插值到相同的网格节点；所述步骤3为：将步骤2所得的每一时间步长对应的波数域结果乘以对应的插值算子，将波数域结果插值到相同的网格节点上；具体方法如下：由于交错网格中v_x和v_z分量定义在不同的网格节点上，网格节点O处的v_x分量用带水平方向插值算子的2D傅里叶反变换进行估计：$v_x(x+\Delta x/2,z)=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}{\tilde{V}}_x(k_x,k_z)e^{i2\pi [k_x(x+\Delta x/2)+k_{z}z]}{\mathrm{dk}}_x{\mathrm{dk}}_z$网格节点O处的v_z分量用带垂直方向插值算子的2D傅里叶反变换进行估计：$v_z(x+\Delta x/2,z)=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}{\tilde{V}}_z(k_x,k_z)e^{i2\pi [k_{x}x+k_z(z-\Delta z/2)]}{\mathrm{dk}}_x{\mathrm{dk}}_z,$其中，和分别与v_x(x,z)和v_z(x+Δx/2,z+Δz/2)对应；在波数域，将乘以插值算子和将乘以插值算子可将它们插值到相同位置；步骤4：利用归一化波数进行保幅波场分离；所述步骤4为：将步骤3所得插值后的波数域波场分量利用归一化波数在波数域实现保幅波场分离；具体方法如下：对空间域波场快照V求散度得P波$P=\frac{\partial V_x}{\partial x}+\frac{\partial V_z}{\partial z},s$对空间域波场快照V求旋度得S波$S=(\frac{\partial V_x}{\partial z}-\frac{\partial V_z}{\partial x})a_y,$其中，V_x、V_z分别为波场V的水平分量和垂直分量，a_y为y方向的单位矢量，此时S波可视为标量波；将对波场快照求散度所得P波波场进行2D傅里叶正变换得对应波数域结果$\tilde{P}={\mathrm{ik}}_x{\tilde{V}}_x+{\mathrm{ik}}_z{\tilde{V}}_z,$将对波场快照求旋度所得S波波场进行2D傅里叶正变换得对应波数域结果$\tilde{S}=({\mathrm{ik}}_z{\tilde{V}}_x-{\mathrm{ik}}_x{\tilde{V}}_z)a_y,$其中，和分别为空间域波场的x与z分量对应的波数域结果，k_x与k_z分别为x方向与z方向的波数；在空间域求偏导数等价于在波数域乘以对应的波数，这会使得分离后的波场相对于原始波场在振幅与单位等方面均发生改变；为了校正求偏导数对波场强度与单位的影响，将和带入P波保幅波场分离公式$\tilde{P}=i\frac{k_x}{k}{\tilde{V}}_x+i\frac{k_z}{k}{\tilde{V}}_z$得保幅的P波结果；将和带入S波保幅波场分离公式$\tilde{S}=i(\frac{k_z}{k}{\tilde{V}}_x-\frac{k_x}{k}{\tilde{V}}_z)a_y$得保幅的S波结果；其中，为波数的绝对值；步骤5：利用傅里叶反变换将分离结果变换到空间域获得保幅的P波和S波；所述步骤5为：将步骤4所得波数域分离的结果利用2D傅里叶反变换变到空间域，从而获得保幅的P波和S波波场。