非等边长网格波动方程有限差分模板优化设计方法

摘要

本发明提供一种非等边长网格波动方程有限差分模板优化设计方法，其中，该方法包括以下步骤：根据时间采样间隔和空间采样间隔对实际地质模型的模拟区域进行网格剖分；根据给定的最大允许误差和波数范围，对不同的声波速度求出对应的算子长度；基于最小二乘优化的时空域有限差分法和所述的对应的算子长度，获取网格的优化有限差分系数；将得到的优化有限差分系数带入差分格式的波动方程，进行波动方程正演模拟。本发明方法将只适用于正方形和正方体网格的时空域有限差分方法扩展到了矩形或长方体网格中，满足实际生产中特定的模拟精度要求和节省计算量的需求。

说明书

技术领域

本发明涉及地震波正演数值模拟技术领域，特别涉及一种非等边长网格波动方程 有限差分模板优化设计方法。

背景技术

地震波正演数值模拟技术是在复杂地质模型(包括各向同性介质、各向异性介质、 Biot多相各向异性介质、随机孔洞介质等)已知的情况下，利用数值计算方法使波在 这种介质中传播，经地下地质构造的多次透射、反射、散射，被地表或地下布置的检 波器接收的过程。利用精确的波动方程数值求解来模拟地下复杂地质构造的地震响 应，为研究地震波传播机理、地震资料的特殊处理方法以及复杂地层的解释等许多方 面提供更为科学的数学物理依据。近年来，波动方程数值模拟方法被广泛应用于逆时 偏移和全波形反演中。

波动方程正演有多种方法，较常见的有：有限差分方法、伪谱法、有限元法、边 界元法、谱元法等等。其中有限差分方法因其计算量小、计算效率高、可以适应较复 杂速度模型而被广泛使用。有限差分法根据不同的标准可以分为：显式有限差分和隐 式有限差分；规则网格有限差分，交错网格有限差分和旋转交错网格有限差分。有限 差分法中，差分系数可以通过泰勒级数展开或最优化方法求得，分别对应以泰勒级数 展开为基础的有限差分和以最优化为基础的有限差分。常规有限差分法中，差分系数 是通过极小化空间域的频散关系得到的。近年来，出现了一种时空域有限差分法，该 方法通过极小化时间域和空间域的频散关系来求取差分系数，具有更高的模拟精度和 更好的稳定性。

目前的时空域有限差分方法要求各方向上的空间采样间隔相等，也就是需要把模 型剖分为正方形或正方体网格。而在实际生产中，为了满足特定的精度要求或为了节 省计算量，我们常常需要将模型剖分为矩形或长方体网格(通常是深度方向的网格间 距不同于水平方向)。适用于正方形和正方体网格的时空域有限差分方法，不能满足 特定的精度要求或不能节省计算量。

发明内容

本发明实施例提供了一种非等边长网格波动方程有限差分模板优化设计方法，将 只适用于正方形和正方体网格的时空域有限差分方法扩展到矩形或长方体网格中，满 足实际生产中特定的模拟精度要求和节省计算量的需求，该方法包括：

根据时间采样间隔和空间采样间隔对实际地质模型的模拟区域进行网格剖分；

根据给定的最大允许误差和波数范围，对不同的声波速度求出对应的算子长度；

基于最小二乘优化的时空域有限差分法和所述的对应的算子长度，获取网格的优 化有限差分系数；

将得到的优化有限差分系数带入差分格式的波动方程，进行波动方程正演模拟。

在一个实施例中，所述根据时间采样间隔和空间采样间隔对实际地质模型的模拟 区域进行差分网格剖分，包括：将实际地质模型的模拟区域剖分成长方体网格或矩形 网格；

基于最小二乘优化的时空域有限差分法和所述的对应的算子长度，获取网格的优 化有限差分系数，包括：获取长方体网格的优化有限差分系数或矩形网格的优化有限 差分系数。

在一个实施例中，所述获取长方体网格的优化有限差分系数，按照如下公式计算：

其中，

其中，b为波数，M为算子长度，a_m为优化后的有限差分系数，θ为平面波传播 方向与水平面的夹角，θ∈[0，π]；φ为平面波传播的方位角，φ∈[0，2π]；V为声波速度，τ为时间采样间隔，h为x、y方向采样间隔；c＝Δz/h，c为参数变 量，Δz为z方向采样间隔；β＝kh，k为参数变量，β为波数范围，β∈[0，b]；m 为参数变量，m为整数，m∈[1，M]；n为参数变量，n为整数，n∈[1，M]。

在一个实施例中，所述长方体网格的优化有限差分系数对应的最大误差满足如下 约束条件：

ξ_1max<η；

其中，ξ_1max为长方体网格的优化有限差分系数对应的最大误差；η是最大允许误 5差。

在一个实施例中，所述长方体网格的优化有限差分系数对应的最大误差ξ_1max按如 下公式计算：

${\xi }_{1\max }=\underset{\beta \in [0,b],\theta \in \left[\mathrm{0,2}\pi \right]}{\max }\left|{\xi }_1(\beta ,\theta )\right|;$

其中，

${\xi }_1(\beta ,\theta )=\frac{2}{\mathrm{r\beta }}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\Sigma }}{a}_m\left(\begin{array}{c}\left(\begin{array}{c}\mathrm{si}{n}^2\left(\frac{\mathrm{m\beta }}{2}\cos \theta \cos \phi \right)-\mathrm{si}{n}^2\left(\frac{\mathrm{m\beta }}{2}\cos \theta \sin \phi \right)\\ -\frac{1}{c^2}\mathrm{si}{n}^2\left(\frac{\mathrm{mc\beta }}{2}\sin \theta \right)\end{array}\right)\end{array}\right)}-1.$

在一个实施例中，所述将得到的优化有限差分系数带入差分格式的波动方程，进 行波动方程正演模拟，包括：将得到的长方体网格的优化有限差分系数代入差分格式 的三维声波波动方程中，进行三维声波波动方程正演模拟；所述差分格式的三维声波 波动方程为：

$\left(\begin{array}{c}\frac{1}{V^2{\tau }^2}(P_{x,y,z}^{t-1}-2P_{x,y,z}^t+P_{x,y,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\Sigma }}{a}_m(P_{x,y,z-m}^t+P_{x,y,z+m}^t)]\\ +\frac{1}{h^2}[2a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\Sigma }}{a}_m(P_{x-m,y,z}^t+P_{x+m,y,z}^t+P_{x,y-m,z}^t+P_{x,y+m,z}^t)]\end{array}\right);$

其中，P为声压。所述获取矩形网格的优化有限差分系数，按照如下公式计算：

其中，

其中，b为波数，M为算子长度，a_m为优化后的有限差分系数，θ为平面波传播 方向与水平面的夹角，θ∈[0，π]；V为声波速度，τ为时间采样间隔，h为 x方向采样间隔；c＝Δz/h，c为参数变量，Δz为z方向采样间隔；β＝kh，k为参 数变量，β为波数范围，β∈[0，b]；m为参数变量，m为整数，m∈[1，M]；n为 参数变量，n为整数，n∈[1，M]。

在一个实施例中，所述矩形网格的优化有限差分系数对应的最大误差满足如下约 束条件：

ξ_2max<η；

其中，ξ_2max为矩形网格的优化有限差分系数对应的最大误差；η是最大允许误差。

在一个实施例中，所述矩形网格的优化有限差分系数对应的最大误差ξ_2max按如 下公式计算：

${\xi }_{2\max }=\underset{\beta \in [0,b],\theta \in \left[\mathrm{0,2}\pi \right]}{\max }\left|{\xi }_2(\beta ,\theta )\right|;$

其中，${\xi }_2(\beta ,\theta )=\frac{2}{\mathrm{r\beta }}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\Sigma }}{a}_m[\mathrm{si}{n}^2\left(\frac{\mathrm{m\beta }}{2}\cos \theta \right)+\frac{1}{c^2}\mathrm{si}{n}^2\left(\frac{\mathrm{mc\beta }}{2}\sin \theta \right)]}-1.$

在一个实施例中，所述将得到的优化有限差分系数带入差分格式的波动方程中， 进行波动方程正演模拟，还包括：将得到的矩形网格的优化有限差分系数代入差分格 式的二维声波波动方程中，进行二维声波波动方程正演模拟；所述差分格式的二维声 波波动方程为：

$\left(\begin{array}{c}\frac{1}{V^2{\tau }^2}(P_{x,z}^{t-1}-2P_{x,z}^t+P_{x,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\Sigma }}{a}_m(P_{x,z-m}^t+P_{x,z+m}^t)]\\ +\frac{1}{h^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\Sigma }}{a}_m(P_{x-m,z}^t+P_{x+m,z}^t)]\end{array}\right);$

其中，P为声压。

在本发明实施例中，提出了一种非等边长网格波动方程有限差分模板优化设计方 法，该方法将只适用于正方形和正方体网格的时空域有限差分方法扩展到了矩形或长 方体网格中，满足实际生产中特定的模拟精度要求和节省计算量的需求。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不 构成对本发明的限定。在附图中：

图1是本发明实施例提供的一种非等边长网格波动方程有限差分模板优化设计 方法流程图；

图2是本发明实施例提供的由传统的有限差分方法获得的误差曲线随传播方向 的变化规律图；

图3是本发明实施例提供的由非等边长网格波动方程有限差分模板优化设计方 法获得的误差曲线随传播方向的变化规律图；

图4是本发明实施例提供的由传统的有限差分方法获得的误差曲线随算子长度 的变化规律图；

图5是本发明实施例提供的由非等边长网格波动方程有限差分模板优化设计方 法获得的误差曲线随算子长度的变化规律图；

图6是本发明实施例提供的通过传统的有限差分方法和非等边长网格波动方程 有限差分模板优化设计方法分别获得的均匀声波介质在1.0s时刻的波场快照对比图；

图7是图6中虚线所在位置的波形对比图；

图8是本发明实施例提供的通过传统的有限差分方法和非等边长网格波动方程 有限差分模板优化设计方法分别获得的Marmousi模型在1.0s时刻的波场快照对比 图；

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图， 对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发 明，但并不作为对本发明的限定。

发明人发现，目前的时空域有限差分方法要求各方向上的空间采样间隔相等，也 就是需要把模型剖分为正方形或正方体网格，但是该种方法不能够满足特定的精度要 求或不能够节省计算量。如果将模型剖分为矩形或长方体网格，则可以满足特定的精 度要求或节省计算量。基于此，本发明提出一种非等边长网格波动方程有限差分模板 优化设计方法。

图1是本发明实施例提供的一种非等边长网格波动方程有限差分模板优化设计 方法流程图，如图1所示，该方法包括：

步骤101：根据时间采样间隔和空间采样间隔对实际地质模型的模拟区域进行网 格剖分；

步骤102：根据给定的最大允许误差和波数范围，对不同的声波速度求出对应的 算子长度；

步骤103：基于最小二乘优化的时空域有限差分法和所述的对应的算子长度，获 取网格的优化有限差分系数；

步骤104：将得到的优化有限差分系数带入差分格式的波动方程，进行波动方程 正演模拟。

具体实施时，需要将实际地质模型的模拟区域(或计算区域)剖分为矩形网格或 长方体网格(通常是深度方向的网格间距不同于水平方向)。

三维声波波动方程如下：

$\frac{{\partial }^{2}P}{\partial x^2}+\frac{{\partial }^{2}P}{\partial y^2}+\frac{{\partial }^{2}P}{\partial z^2}=\frac{1}{V^2}\frac{{\partial }^{2}P}{\partial t^2}---\left(1\right)$

式中，P代表声压，V代表声波速度。

以z方向网格间距不同于x、y方向的情况为例，具体的长方体网格的优化有限 差分格式如下：

$\left(\begin{array}{c}\frac{1}{V^2{\tau }^2}(P_{x,y,z}^{t-1}-2P_{x,y,z}^t+P_{x,y,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\Sigma }}{a}_m(P_{x,y,z-m}^t+P_{x,y,z+m}^t)]\\ +\frac{1}{h^2}[2a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\Sigma }}{a}_m(P_{x-m,y,z}^t+P_{x+m,y,z}^t+P_{x,y-m,z}^t+P_{x,y+m,z}^t)]\end{array}\right)---\left(2\right)$

其中，τ为时间采样间隔；h为x、y方向采样间隔；c＝Δz/h，c为参数变量， Δz为z方向采样间隔；a_m为优化后的有限差分系数，M表示算子长度；m为参数变 量，m为整数，m∈[1，M]。

为了简化符号，我们将简写为根据平面波理论，我们令：

$P_{m,l,j}^n=e^{i[k_x(x+\mathrm{mh})+k_y(y+\mathrm{lh})+k_z(z+\mathrm{jch})-\omega (t+\mathrm{n\tau })]}---\left(3\right)$

其中，

$\left(\begin{array}{c}{k}_x\text{=kcos}\left(\theta \right)\cos \left(\phi \right)\\ k_y=k\cos \left(\theta \right)\sin \left(\phi \right)\\ k_z=k\sin \left(\theta \right)\end{array}\right)---\left(4\right)$

其中，i，l，j，n为参数变量，n为整数，n∈[1，M]；k为参数变量；θ∈[0,π] 为平面波传播方向与水平面的夹角，φ∈[0,2π]为平面波传播的方位角。将方程(3)、 (4)代入方程(2)，并适当化简，可得：

$\left(\begin{array}{c}2r^{-2}[\cos \left(\mathrm{\omega \tau }\right)-1]=\frac{1}{c^2}{a}_0+\underset{m=1}{\overset{M}{\Sigma }}\frac{2}{c^2}{a}_m\left[\cos \left(m{k}_{z}h\right)\right]\\ +2a_0+\underset{m=1}{\overset{M}{\Sigma }}2a_m[\cos \left(m{k}_{x}h\right)+\cos \left(m{k}_{y}h\right)]\end{array}\right)---\left(5\right)$

其中，再根据约束条件；

$a_0+2\underset{m=1}{\overset{M}{\Sigma }}{a}_m=0---\left(6\right)$

可得：

$r^{-2}[1-\cos \left(\mathrm{\omega \tau }\right)]=\underset{m=1}{\overset{M}{\Sigma }}{a}_m\left(\begin{array}{c}2-\cos \left(\mathrm{m\beta }\cos \theta \cos \phi \right)\text{-cos}\left(\mathrm{m\beta }\cos \theta \sin \phi \right)\\ +\frac{1}{c^2}[1-\cos \left(\mathrm{mc\beta }\sin \theta \right)]\end{array}\right)---\left(7\right)$

其中，β＝kh，β为波数范围，β∈[0，b]。优化有限差分算法就是要求出一套 差分系数使得方程(7)两端误差最小。根据最小二乘理论，此问题可转化为求解如 下方程组：

其中，

通过求解线性方程组(8)，我们可以得到长方体网格的优化差分系数。再将差分 系数代入到方程(2)中，就可以进行三维声波波动方程的求解。

为了保证本发明方法的有效性，需要在前面的最优化过程中加入下面的约束条 件，即长方体网格的优化有限差分系数对应的最大误差满足如下约束条件：

ξ_1max<η                    (10)

我们通过下式计算优化差分系数对应的最大误差：

${\xi }_{1\max }=\underset{\beta \in [0,b],\theta \in \left[\mathrm{0,2}\pi \right]}{\max }\left|{\xi }_1(\beta ,\theta )\right|---\left(11\right)$

其中，

${\xi }_1(\beta ,\theta )=\frac{2}{\mathrm{r\beta }}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\Sigma }}{a}_m\left(\begin{array}{c}\left(\begin{array}{c}\mathrm{si}{n}^2\left(\frac{\mathrm{m\beta }}{2}\cos \theta \cos \phi \right)-\mathrm{si}{n}^2\left(\frac{\mathrm{m\beta }}{2}\cos \theta \sin \phi \right)\\ -\frac{1}{c^2}\mathrm{si}{n}^2\left(\frac{\mathrm{mc\beta }}{2}\sin \theta \right)\end{array}\right)\end{array}\right)}-1---\left(12\right)$

其中，η是最大允许误差。

ξ_1max只与b和M有关。当M固定时，ξ_1max决定于b，ξ_1max随着b增加而变大。 因此，若初始的b不满足方程(10)，我们应通过逐渐减小b直到ξ_1max<η来获得最 佳b。另外，当b固定，ξ_1max只与M有关，ξ_1max随着M增加而变小。所以，若初始 的M不满足方程(10)，我们应通过逐渐增大M直到ξ_1max<η来获得最佳M。

对于二维矩形网格，以z方向网格间距不同于x方向的情况为例，具体的优化有 限差分格式如下：

$\left(\begin{array}{c}\frac{1}{V^2{\tau }^2}(P_{x,z}^{t-1}-2P_{x,z}^t+P_{x,z}^{t+1})=\frac{1}{{\left(\mathrm{ch}\right)}^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\Sigma }}{a}_m(P_{x,z-m}^t+P_{x,z+m}^t)]\\ +\frac{1}{h^2}[a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\Sigma }}{a}_m(P_{x-m,z}^t+P_{x+m,z}^t)]\end{array}\right)---\left(13\right)$

式中，P代表声压，V代表声波速度；其中，τ为时间采样间隔；h为x方向采 样间隔；c＝Δz/h，c为参数变量，Δz为z方向采样间隔；a_m为优化后的有限差分 系数，M表示算子长度；m为参数变量，m为整数，m∈[1，M]。

我们可以通过求解下列线性方程来获取矩形网格的优化的差分系数：

其中，

通过求解线性方程组(15)，我们可以得到矩形网格的优化差分系数。再将差分系 数代入到方程(13)中，就可以进行二维声波波动方程的求解。

同样的，为了保证本发明方法的有效性，需要在前面的最优化过程中加入下面的 约束条件，即矩形网格的优化有限差分系数对应的最大误差满足如下约束条件：

ξ_2max<η                  (16)

我们通过下式计算优化差分系数对应的最大误差：

${\xi }_{2\max }=\underset{\beta \in [0,b],\theta \in \left[\mathrm{0,2}\pi \right]}{\max }\left|{\xi }_2(\beta ,\theta )\right|---\left(17\right)$

其中，

${\xi }_2(\beta ,\theta )=\frac{2}{\mathrm{r\beta }}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\Sigma }}{a}_m[\mathrm{si}{n}^2\left(\frac{\mathrm{m\beta }}{2}\cos \theta \right)+\frac{1}{c^2}\mathrm{si}{n}^2\left(\frac{\mathrm{mc\beta }}{2}\sin \theta \right)]}-1---\left(18\right)$

其中，η是最大允许误差。

下面以一个均匀声波介质为例来说明本发明的优势。纵波速度为1500m/s，时间 采样间隔为1ms，x方向采样间隔为10m，z方向采样间隔为5m，算子长度M∈[2，10]。

规定b＝2.74和M＝8，采用本发明方法(即非等边长网格波动方程有限差分模板 优化设计方法)和传统的有限差分方法获得误差曲线随传播方向的变化规律，分别如 图2和图3所示。规定η＝10^-2.5≈0.003，采用本发明方法和传统的有限差分方法获 得误差曲线随算子长度的变化规律，分别如图4和5所示。由图2至图5可知，与基 于泰勒的时空域有限差分方法(传统的有限差分方法)相比，本发明方法具有更宽的 有效频带及更小的数值频散。规定η＝10^-2.5≈0.003，采用传统的有限差分方法和本 发明方法获得均匀声波介质的波场快照对比图，如图6所示。其中，图6中的左边图 为采用传统的有限差分方法获得均匀声波介质的波场快照，右边图为采用本发明方法 获得的均匀声波介质的波场快照。图7为图6中虚线所在位置的波形对比图。由图6 和图7可知，当算子长度相同时，本发明方法模拟精度更高。

为了进一步说明本发明的有效性，我们用Marmousi模型进行正演。时间采样间 隔为1ms，x方向采样间隔为12m，z方向采样间隔为9m。传统的有限差分法的算 子长度M＝9。时空域有限差分采用变算子长度的算法，根据不同的速度选用不同长 度的算子，M∈[2，9]。图8是本发明实施例提供的通过传统的有限差分方法和本发 明方法分别获得的Marmousi模型在1.0s时刻的波场快照对比图，上边图为通过传统 的有限差分方法获得的Marmousi模型在1.0s时刻的波场快照，下边图为通过本发明 方法分别获得的Marmousi模型在1.0s时刻的波场快照；由图8可知，本发明方法适 用于复杂介质的正演模拟，并且其模拟效果明显优于传统的有限差分法。

综上所述，本发明方法具有以下优点：1.减小了中高频段的数值频散。2.以时 空域频散关系为基础更符合实际，模拟精度更高。3.适用于矩形和长方体网格，可 以满足实际生产中特定的精度要求和节省计算量的需求。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技 术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内， 所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。