阅读说明：本技术 地震波模拟方法及系统 (Seismic wave simulation method and system ) 是由 刘炯 刘喜武 刘振峰 张金强 钱恪然 刘志远 霍志周 刘宇巍 张远银 于 2018-06-27 设计创作，主要内容包括：公开了一种地震波模拟方法及系统。该方法可以包括：步骤1：建立地震波波动方程,以g表示其中任意一个变量；步骤2：对三维空间进行网格剖分；步骤3：x方向采用傅里叶方程求解下一时刻点的内部网格点的变量g,z方向采用采用切比雪夫方程求解下一时刻点的内部网格点的变量g；步骤4：针对某一时刻点进行边界方程求解,获得下一时刻点的边界上的变量g；步骤5：随时间步进,重复步骤3-4,直至获得每个时刻点的内部网格点和边界上的变量g。本发明通过傅里叶伪谱法计算水平导数,用切比雪夫伪谱法计算垂直方向导数,在解决自由面模拟的同时,保证了模拟的高精度。(A method and system for simulating seismic waves are disclosed. The method can comprise the following steps: step 1: establishing a seismic wave equation, and expressing any variable in g; step 2: mesh generation is carried out on the three-dimensional space; and step 3: solving a variable g of an internal grid point of the next time point by adopting a Fourier equation in the x direction, and solving a variable g of the internal grid point of the next time point by adopting a Chebyshev equation in the z direction; and 4, step 4: solving a boundary equation aiming at a certain time point to obtain a variable g on the boundary of the next time point; and 5: and repeating the steps 3-4 with the time stepping until the variable g on the internal grid points and the boundary of each time point is obtained. According to the method, the horizontal derivative is calculated by a Fourier pseudo-spectrum method, the vertical derivative is calculated by a Chebyshev pseudo-spectrum method, and high simulation precision is ensured while free surface simulation is solved.)

地震波模拟方法及系统

技术领域

本发明涉及油汽地球物理技术领域，更具体地，涉及一种地震波模拟方法及系统。

背景技术

伪谱法与有限差分法、有限元法并列为地震波数值模拟的三种主要方法。它通过数学变换的方法把物理变量对空间的微分转化为转换域中的代数运算，然后再把结果通过逆变换转换到物理空间，从而求得对应量的空间微分。相对于有限差分法和有限元法，伪谱法具有精度高的特点，运算速度比较快的特点。因此，伪谱法在地震勘探中有着比较广泛的应用。

目前地震领域中的伪谱法主要指傅里叶伪谱法。傅里叶伪谱法中的基函数为三角函数，具有周期特性，因此，傅里叶伪谱法求得的空间导数值也是周期的，这种伪谱法一般用于周期边界条件的情况，而不适合有固定条件的情况。地球介质是有边界的，如地表、固定底面，因此边界条件往往在波场模拟中需要考虑。用传统的傅里叶伪谱法不能很好处理这类含边界的问题。因此，有必要开发一种地震波模拟方法及系统。

公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。

发明内容

本发明提出了一种地震波模拟方法及系统，其能够通过傅里叶伪谱法计算水平导数，用切比雪夫伪谱法计算垂直方向导数，在解决自由面模拟的同时，保证了模拟的高精度。

根据本发明的一方面，提出了一种地震波模拟方法。所述方法可以包括：步骤1：建立含自由表面的地层介质中的地震波波动方程，以g表示其中任意一个变量；步骤2：对三维空间进行网格剖分，其中x方向采用等距网格，z方向采用高斯-切比雪夫-洛巴托点获得空间离散点；步骤3：针对某一时刻点的内部网格点，x方向采用傅里叶方程求解下一时刻点的内部网格点的变量g，z方向采用采用切比雪夫方程求解下一时刻点的内部网格点的变量g；步骤4：施加物理边界条件和数值边界条件，针对某一时刻点进行边界方程求解，获得下一时刻点的边界上的变量g；步骤5：随时间步进，重复步骤3-4，直至获得每个时刻点的内部网格点和边界上的变量g。

优选地，x方向采用傅里叶方程求解变量g的一阶微分包括：

其中，k代表波数，Δk为波数间隔，Δk＝2π/(NΔx)，i是虚数单位，G(lΔk)表示函数g(nΔx)的傅里叶变换，g(nΔx)为g在x方向的离散表达式。

优选地，g(nΔx)的表达式为：

式中，n,l为节点索引，N为傅里叶多项式的最高阶数。

优选地，z方向的高斯-切比雪夫-洛巴托离散点的表达式为：

其中，计算区间z被限制在[-1,1]之间，端点z₀＝1，z_M＝-1，M是运算中所取切比雪夫多项式的阶数。

优选地，z方向采用切比雪夫方程求解变量g的一阶微分包括：

其中，

的离散表达式为

其中，D_c为(M+1)×(M+1)的矩阵，

根据本发明的另一方面，提出了一种地震波模拟系统，其上存储有计算机程序，其特征在于，所述程序被处理器执行时实现以下步骤：步骤1：建立含自由表面的地层介质中的地震波波动方程，以g表示其中任意一个变量；步骤2：对三维空间进行网格剖分，其中x方向采用等距网格，z方向采用高斯-切比雪夫-洛巴托点获得空间离散点；步骤3：针对某一时刻点的内部网格点，x方向采用傅里叶方程求解下一时刻点的内部网格点的变量g，z方向采用采用切比雪夫方程求解下一时刻点的内部网格点的变量g；步骤4：施加物理边界条件和数值边界条件，针对某一时刻点进行边界方程求解，获得下一时刻点的边界上的变量g；步骤5：随时间步进，重复步骤3-4，直至获得每个时刻点的内部网格点和边界上的变量g。

优选地，x方向采用傅里叶方程求解变量g的一阶微分包括：

其中，k代表波数，Δk为波数间隔，Δk＝2π/(NΔx)，i是虚数单位，G(lΔk)表示函数g(nΔx)的傅里叶变换，g(nΔx)为g在x方向的离散表达式。

优选地，g(nΔx)的表达式为：

式中，n,l为节点索引，N为傅里叶多项式的最高阶数。

优选地，z方向的高斯-切比雪夫-洛巴托离散点的表达式为：

其中，计算区间z被限制在[-1,1]之间，端点z₀＝1，z_M＝-1，M是运算中所取切比雪夫多项式的阶数。

优选地，z方向采用切比雪夫方程求解变量g的一阶微分包括：

其中，

的离散表达式为

其中，D_c为(M+1)×(M+1)的矩阵，

本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的

具体实施方式

中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述，这些附图和具体实施方式共同用于解释本发明的特定原理。

附图说明

通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施例中，相同的参考标号通常代表相同部件。

图1示出了根据本发明的地震波模拟方法的步骤的流程图。

图2示出了根据本发明的一个实施例的网格剖分的示意图。

图3示出了根据本发明的一个实施例的二维x-z空间模型的示意图。

图4示出了根据本发明的一个实施例的震源子波在0.35秒时刻的波场快照的示意图。

图5a、5b、5c分别示出了根据本发明的一个实施例的模拟得到的检波器R1、R2、R3的位移记录和理论结果的对比图。

具体实施方式

下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。

图1示出了根据本发明的地震波模拟方法的步骤的流程图。

在该实施例中，根据本发明的地震波模拟方法可以包括：步骤1：建立含自由表面的地层介质中的地震波波动方程，以g表示其中任意一个变量；步骤2：对三维空间进行网格剖分，其中x方向采用等距网格，z方向采用高斯-切比雪夫-洛巴托点获得空间离散点；步骤3：针对某一时刻点的内部网格点，x方向采用傅里叶方程求解下一时刻点的内部网格点的变量g，z方向采用采用切比雪夫方程求解下一时刻点的内部网格点的变量g；步骤4：施加物理边界条件和数值边界条件，针对某一时刻点进行边界方程求解，获得下一时刻点的边界上的变量g；步骤5：随时间步进，重复步骤3-4，直至获得每个时刻点的内部网格点和边界上的变量g。

在一个示例中，x方向采用傅里叶方程求解变量g的一阶微分包括：

其中，k代表波数，Δk为波数间隔，Δk＝2π/(NΔx)，i是虚数单位，G(lΔk)表示函数g(nΔx)的傅里叶变换，g(nΔx)为g在x方向的离散表达式。

在一个示例中，g(nΔx)的表达式为：

式中，n,l为节点索引，N为傅里叶多项式的最高阶数。

在一个示例中，z方向的高斯-切比雪夫-洛巴托离散点的表达式为：

其中，计算区间z被限制在[-1,1]之间，端点z₀＝1，z_M＝-1，M是运算中所取切比雪夫多项式的阶数。

在一个示例中，z方向采用切比雪夫方程求解变量g的一阶微分包括：

其中，

的离散表达式为

其中，D_c为(M+1)×(M+1)的矩阵，

具体地，以一阶SH波方程为例，根据本发明的地震波模拟方法可以包括：

对三维空间进行网格剖分，其中x方向采用等距网格，z方向采用高斯-切比雪夫-洛巴托点离散。根据弹性力学原理，针对二维x-z空间建立一阶SH波方程为公式(6)：

其中，V表示介质沿垂直x-z平面的y方向的速度，τ_xy、τ_yz是应力分量，f_y是单位质量上沿y方向的外力，ρ代表介质密度，μ是介质的剪切模量。

针对某一时刻点，对二维x-z空间的内部网格点进行一阶SH波方程求解，获得下一时刻点的内部网格点的速度、应力，其中，x方向采用傅里叶方程求解，z方向采用切比雪夫方程求解。

对于x方向，采用等距网格，得到对应的函数离散值τ_xy、V，然后将离散后的函数变换至波数域分别为为公式(7)、(8)：

其中，

表示函数τ_xy、V的傅里叶变换。根据函数微分的傅里叶变换和函数的傅里叶变换的关系，将和ilΔk相乘，并将结果逆傅里叶变换回空间域，得到函数的1阶空间微分结果，即通过公式(9)求解x方向的应力分量的一阶微分：

其中，τ_xy是应力分量，k代表波数，波数间隔Δk＝2π/(NΔx)，i是虚数单位，通过公式(10)求解速度的一阶微分：

其中，V表示介质沿垂直x-z平面的y方向的速度，k代表波数，波数间隔Δk＝2π/(NΔx)。

对于z方向，z方向采用高斯-切比雪夫-洛巴托点离散，其表达形式为公式(3)，其中计算区间z被限制在[-1,1]之间，端点z₀＝1，z_M＝-1，[a,b]为z方向的物理区间(a、b可以为有界的任意实数)，通过坐标转换映射至计算坐标z’，其中z’属于[-1,1]，然后通过公式(11)求解z方向的应力分量的一阶微分：

其中

的离散表达式为

其中，

公式(11)中

[a,b]为z方向的物理区间。

通过公式(13)求解速度在z方向的一阶微分：

其中

的离散表达式为

其中，

公式(14)中[a,b]为z方向的物理区间。

求得和

后，公式(6)右边项的值就可以得到了(f_y，ρ和μ是已知的)，于是公式(6)由偏微分方程组变为了常微分方程组，再利用数值积分技术(如欧拉法，4阶龙格库塔法等)求得下一时刻的V、τ_xy和τ_yz。

施加物理边界条件和数值边界条件，针对某一时刻点，对二维x-z空间的边界进行边界方程(15)求解：

τ_yz＝0,当z＝0 (15)

获得下一时刻点的边界上的位移和应力。

随时间步进，获得每个时刻点的内部网格点的速度、应力和边界上的位移、应力。根据每个时刻点的内部网格点的速度、应力和边界上的位移、应力，获得SH波的波场记录，得到空间各场量移随时间的变化关系。

本方法通过傅里叶伪谱法计算水平导数，用切比雪夫伪谱法计算垂直方向导数，在解决自由面模拟的同时，保证了模拟的高精度。

应用示例

为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。

图2示出了根据本发明的一个实施例的网格剖分的示意图。

图3示出了根据本发明的一个实施例的二维x-z空间模型的示意图。

图4示出了根据本发明的一个实施例的震源子波在0.35秒时刻的波场快照的示意图。

根据本发明的地震波模拟方法包括：

对三维空间进行网格剖分，其中x方向采用等距网格，z方向采用高斯-切比雪夫-洛巴托点离散，如图2所示，对如图2所示的二维x-z空间进行SH波数值模拟，安置了3个虚拟检波器R1、R2、R3，来记录空间点的位移，如图3所示，模拟中震源采用50Hz的Ricker子波，如图4所示。

根据弹性力学原理，针对二维x-z空间建立一阶SH波方程为公式(6)。

针对某一时刻点，对二维x-z空间的内部网格点进行一阶SH波方程求解，获得下一时刻点的内部网格点的速度、应力，其中，x方向采用傅里叶方程求解，z方向采用切比雪夫方程求解。

对于x方向，采用等距网格，得到对应的函数离散值τ_xy、V，然后将离散后的函数变换至波数域为公式(7)、(8)。根据函数微分的傅里叶变换和函数的傅里叶变换的关系，将

和ilΔk相乘，并将结果逆傅里叶变换回空间域，得到函数的1阶空间微分结果，即通过公式(9)求解x方向的应力分量的一阶微分，通过公式(10)求解速度的一阶微分。

对于z方向，z方向采用高斯-切比雪夫-洛巴托点离散，其中计算区间z被限制在[-1,1]之间，端点z₀＝1，z_M＝-1，[a,b]为z方向的物理区间(a、b可以为有界的任意实数)，通过坐标转换

映射至计算坐标z’，其中z’属于[-1,1]，然后通过公式(11)求解z方向的应力分量的一阶微分，通过公式(13)求解速度在z方向的一阶微分，其中，公式(12)、(14)表示

的离散表达式。

求得

和

后，公式(6)右边项的值就可以得到了(f_y，ρ和μ是已知的)，于是公式(6)由偏微分方程组变为了常微分方程组，再利用数值积分技术(如欧拉法，4阶龙格库塔法等)求得下一时刻的V、τ_xy和τ_yz。

施加物理边界条件和数值边界条件，针对某一时刻点，对二维x-z空间的边界进行边界方程(15)求解，获得下一时刻点的边界上的位移和应力。

随时间步进，获得每个时刻点的内部网格点的速度、应力和边界上的位移、应力。根据每个时刻点的内部网格点的速度、应力和边界上的位移、应力，获得SH波的波场记录，得到空间各场量移随时间的变化关系。

图5a、5b、5c分别示出了根据本发明的一个实施例的模拟得到的检波器R1、R2、R3的位移记录和理论结果的对比图，由图可知，数值解和理论解符合得很好，证明本方法具有很高的模拟精度。

综上所述，本发明通过傅里叶伪谱法计算水平导数，用切比雪夫伪谱法计算垂直方向导数，在解决自由面模拟的同时，保证了模拟的高精度。

本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。

根据本发明的一种地震波模拟系统，其上存储有计算机程序，所述程序被处理器执行时实现以下步骤：步骤1：建立含自由表面的地层介质中的地震波波动方程，以g表示其中任意一个变量；步骤2：对三维空间进行网格剖分，其中x方向采用等距网格，z方向采用高斯-切比雪夫-洛巴托点获得空间离散点；步骤3：针对某一时刻点的内部网格点，x方向采用傅里叶方程求解下一时刻点的内部网格点的变量g，z方向采用采用切比雪夫方程求解下一时刻点的内部网格点的变量g；步骤4：施加物理边界条件和数值边界条件，针对某一时刻点进行边界方程求解，获得下一时刻点的边界上的变量g；步骤5：随时间步进，重复步骤3-4，直至获得每个时刻点的内部网格点和边界上的变量g。

在一个示例中，x方向采用傅里叶方程求解变量g的一阶微分包括：

其中，k代表波数，Δk为波数间隔，Δk＝2π/(NΔx)，i是虚数单位，G(lΔk)表示函数g(nΔx)的傅里叶变换，g(nΔx)为g在x方向的离散表达式。

在一个示例中，g(nΔx)的表达式为：

式中，n,l为节点索引，N为傅里叶多项式的最高阶数。

在一个示例中，z方向的高斯-切比雪夫-洛巴托离散点的表达式为：

其中，计算区间z被限制在[-1,1]之间，端点z₀＝1，z_M＝-1，M是运算中所取切比雪夫多项式的阶数。

在一个示例中，z方向采用切比雪夫方程求解变量g的一阶微分包括：

其中，

的离散表达式为

其中，D_c为(M+1)×(M+1)的矩阵，

本系统通过傅里叶伪谱法计算水平导数，用切比雪夫伪谱法计算垂直方向导数，在解决自由面模拟的同时，保证了模拟的高精度。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。