三维Laplace域声波方程数值模拟方法及装置

摘要

本发明提供了一种三维Laplace域声波方程数值模拟方法及装置，该方法包括：基于三维Laplace域声波方程，利用平均导数方法建立包含多个自由参数的27点差分方程；利用所述包含多个自由参数的27点差分方程，通过Laplace域频散分析确定所述多个自由参数的值；利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图。本发明利用包含多个自由参数的27点差分方程进行数值模拟，可以显著提高模拟效率。

说明书

技术领域

本发明涉及声波数值模拟技术领域，尤其涉及一种三维Laplace域声波方程数值模拟方法及装置。

背景技术

三维频率域声波方程多尺度全波形反演，可以基于简单的初始模型，直接利用所接收的地震数据得到关于地下介质参数的定量信息，是一种有效的地震成像方法。但当地震数据缺乏低频信息时，对于频率域多尺度全波形反演来讲，简单的初始模型则不再满足要求，此时需要其它反演方法提供包含长波长速度分量的初始模型。而且，地震数据通常会缺乏低频信息或很难得到可靠的低频信息，所以频率域多尺度全波形反演方法的实际应用通常存在局限性。

为了解决三维频率域声波方程多尺度全波形反演对地震数据低频信息的依赖问题，人们发展了三维Laplace域声波方程全波形反演方法。这种方法对地震数据进行人为衰减，然后从衰减的地震数据中提取低频信息，再利用波形反演方法来构造包含长波长速度分量的初始速度模型。最后，以该初始速度模型作为输入，利用频率域多尺度全波形反演得到关于地下介质参数的定量信息。

三维Laplace域声波方程数值模拟是三维Laplace域声波方程全波形反演方法的基础，由于要对多个Laplace参数进行反演，三维Laplace域声波方程全波形反演方法计算量大，效率低。而且，三维Laplace域声波方程数值模拟使用的常规27点差分方法，在百分之一的频散误差范围内，每一个最小拟波长需要16个采样点，计算效率比较低。

发明内容

本发明提供一种三维Laplace域声波方程数值模拟方法及装置，以提高三维Laplace域声波方程数值模拟的效率。

本发明提供一种三维Laplace域声波方程数值模拟方法，包括：基于三维Laplace域声波方程，利用平均导数方法建立包含多个自由参数的27点差分方程；利用所述包含多个自由参数的27点差分方程，通过Laplace域频散分析确定所述多个自由参数的值；利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图。

一个实施例中，所述包含多个自由参数的27点差分方程为：

$\left(\begin{array}{c}\frac{{P^{\left(1\right)}}_{m+1,l,n}-2{P^{\left(1\right)}}_{m,l,n}+{P^{\left(1\right)}}_{m-1,l,n}}{{\mathrm{\Delta x}}^2}+\frac{{P^{\left(2\right)}}_{m,l+1,n}-2{P^{\left(2\right)}}_{m,l,n}+{P^{\left(2\right)}}_{m,l-1,n}}{{\mathrm{\Delta y}}^2}\\ +\frac{{P^{\left(3\right)}}_{m,l,n+1}-2{P^{\left(3\right)}}_{m,l,n}+{P^{\left(3\right)}}_{m,l,n-1}}{{\mathrm{\Delta z}}^2}-\frac{s^2}{v_{m,l,n}^2}\Sigma =0,\end{array}\right)$

其中，

P⁽¹⁾_m+j,l,n＝α₁(P_m+j,l+1,n+P_m+j,l,n+1+P_m+j,l-1,n+P_m+j,l,n-1)

+α₂(P_m+j,l+1,n+1+P_m+j,l-1,n+1+P_m+j,l+1,n-1+P_m+j,l-1,n-1)

+(1-4α₁-4α₂)P_m+j,l,n,j＝1,0,-1,

P⁽²⁾_m,l+j,n＝α₃(P_m+1,l+j,n+P_m,l+j,n+1+P_m-1,l+j,n+P_m,l+j,n-1)

+α₄(P_m+1,l+j,n+1+P_m+1,l+j,n-1+P_m-1,l+j,n+1+P_m-1,l+j,n-1)

+(1-4α₃-4α₄)P_m,l+j,n,j＝1,0,-1,

P⁽³⁾_m,l,n+j＝α₅(P_m+1,l,n+j+P_m,l+1,n+j+P_m-1,l,n+j+P_m,l-1,n+j)

+α₆(P_m+1,l+1,n+j+P_m+1,l-1,n+j+P_m-1,l+1,n+j+P_m-1,l-1,n+j)

+(1-4α₅-4α₆)P_m,l,n+j,j＝1,0,-1,

$\left(\begin{array}{c}\Sigma ={\alpha }_7{P}_{m,l,n}+{\alpha }_8(P_{m,l+1,n}+P_{m,l,n+1}+P_{m,l-1,n}+P_{m,l,n-1}+P_{m+1,l,n}+P_{m-1,l,n})\\ +\frac{1}{12}(1-{\alpha }_7-6{\alpha }_8)(P_{m+1,l+1,n}+P_{m+1,l,n+1}+P_{m+1,l-1,n}+P_{m+1,l,n-1}+P_{m-1,l+1,n}+P_{m-1,l,n+1}\\ +P_{m-1,l-1,n}+P_{m-1,l,n-1}+P_{m,l+1,n+1}+P_{m,l-1,n+1}+P_{m,l+1,n-1}+P_{m,l-1,n-1}),\end{array}\right)$

P_l,m,n≈P((m-1)Δx,(l-1)Δy,(n-1)Δz),m＝1,...,M；l＝1,...,L；n＝1,...,N,

其中，P表示压力波场；压力波场P的下标m+j,l+j,n+j表示离散点的三个不同指标；压力波场P的上标(1),(2),(3)分别表示对第一下标m+j、第二下标l+j、第三个下标n+j的离散点处的压力波场进行平均；s为Laplace阻尼常数；v_m,l,n表示三个指标分别为m,l,n的离散点处的波速；Δx,Δy,Δz表示离散点在网格方向x,y,z上的间隔；M,L,N分别表示离散点的三个指标的最大值，M,L,N为正整数；α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈为自由参数。

一个实施例中，利用所述包含多个自由参数的27点差分方程，通过Laplace域频散分析确定所述多个自由参数的值，包括：将设定波场代入所述包含多个自由参数的27点差分方程，计算得到包含所述多个自由参数的波场衰减传播数值速度，所述设定波场的拟波长为波振幅衰减到初始波振幅的1/e^2π时波所经过的空间距离；利用所述波场衰减传播数值速度，建立包含所述多个自由参数的频散误差最小化目标函数；利用所述频散误差最小化目标函数确定使频散误差最小的所述自由参数的值。

一个实施例中，归一化的所述波场衰减传播数值速度为：

$\frac{V_{num}(\theta ,\phi ,\tilde{k};{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6,{\alpha }_7,{\alpha }_8)}{v}=\frac{G}{2\pi }{\left(\frac{ϵ}{F}\right)}^{\frac{1}{2}},$

其中，

$\left(\begin{array}{c}ϵ=2[\cosh \left(\frac{2\pi sin\theta cos\phi }{G}\right)-1]E_x+2g_1^2[\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)-1]E_y\\ +2g_2^2[\cosh \left(\frac{2\pi \cos \theta }{g_{2}G}\right)-1]E_z,\end{array}\right)$

$G=\frac{2\pi }{{\mathrm{k\Delta x}}^{\prime }},{\mathrm{\Delta x}}^{\prime }=max\{\Delta x,\Delta y,\Delta z\},\tilde{k}=\frac{1}{G},g_1=\frac{\Delta x}{\Delta y},g_2=\frac{\Delta x}{\Delta z},$

$\left(\begin{array}{c}{E}_x=2{\alpha }_1\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{g_{1}G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)+4{\alpha }_2\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +(1-4{\alpha }_1-4{\alpha }_2),\end{array}\right)$

$\left(\begin{array}{c}{E}_y=2{\alpha }_3\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)+4{\alpha }_4\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +(1-4{\alpha }_3-4{\alpha }_4),\end{array}\right)$

$\left(\begin{array}{c}{E}_z=2{\alpha }_5\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{G})+\cosh (\frac{2\pi sin\theta sin\phi }{g_{1}G}\left)\right)+4{\alpha }_6\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\\ +(1-4{\alpha }_5-4{\alpha }_6),\end{array}\right)$

$\left(\begin{array}{c}F={\alpha }_7+2{\alpha }_8\left(\cosh \right(\frac{2\pi sin\theta cos\phi }{G})+\cosh (\frac{2\pi sin\theta cos\phi }{g_{1}G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)\\ +\frac{1-{\alpha }_7-6{\alpha }_8}{12}\left(\begin{array}{c}\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)+\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\end{array}\right),\end{array}\right)$

其中，θ为传播角；φ为方位角；v为离散点上的波速；α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈为自由参数；为归一化的波场衰减传播数值速度；k为拟波长对应的拟波数；Δx,Δy,Δz表示离散点在网格方向x,y,z上的间隔；

所述频散误差最小化目标函数为：

$E=\int \int \int {[1-\frac{V_{num}(\theta ,\phi ,\tilde{k};{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6,{\alpha }_7,{\alpha }_8)}{v}]}^{2}d\tilde{k}d\theta d\phi ,$

其中，E为衰减传播数值速度频散误差。

一个实施例中，利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图，包括：根据离散速度模型的大小，利用确定所述多个自由参数的值后的27点差分方程建立离散波场线性方程组；将阻尼常数、源子波、所述离散速度模型的网格间隔及所述离散速度模型上各离散点上的波速输入所述离散波场线性方程组，并利用所述离散速度模型的吸收边界条件，依据所述离散波场线性方程组中阻抗矩阵的大小选择求解方法求解所述离散波场线性方程组，得到波场在所述离散速度模型的各离散点上的值。

一个实施例中，利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图，还包括：根据所述离散速度模型的各离散点上的值和检波器在所述离散速度模型中的位置输出所述检波器检测的震波图。

本发明还提供一种三维Laplace域声波方程数值模拟装置，包括：差分方程建立单元，用于执行：基于三维Laplace域声波方程，利用平均导数方法建立包含多个自由参数的27点差分方程；自由参数确定单元，用于执行：利用所述包含多个自由参数的27点差分方程，通过Laplace域频散分析确定所述多个自由参数的值；数值模拟单元，用于执行：利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图。

一个实施例中，所述差分方程建立单元，还用于执行：所述包含多个自由参数的27点差分方程为：

$\left(\begin{array}{c}\frac{{P^{\left(1\right)}}_{m+1,l,n}-2{P^{\left(1\right)}}_{m,l,n}+{P^{\left(1\right)}}_{m-1,l,n}}{{\mathrm{\Delta x}}^2}+\frac{{P^{\left(2\right)}}_{m,l+1,n}-2{P^{\left(2\right)}}_{m,l,n}+{P^{\left(2\right)}}_{m,l-1,n}}{{\mathrm{\Delta y}}^2}\\ +\frac{{P^{\left(3\right)}}_{m,l,n+1}-2{P^{\left(3\right)}}_{m,l,n}+{P^{\left(3\right)}}_{m,l,n-1}}{{\mathrm{\Delta z}}^2}-\frac{s^2}{v_{m,l,n}^2}\Sigma =0,\end{array}\right)$

其中，

P⁽¹⁾_m+j,l,n＝α₁(P_m+j,l+1,n+P_m+j,l,n+1+P_m+j,l-1,n+P_m+j,l,n-1)

+α₂(P_m+j,l+1,n+1+P_m+j,l-1,n+1+P_m+j,l+1,n-1+P_m+j,l-1,n-1)

+(1-4α₁-4α₂)P_m+j,l,n,j＝1,0,-1,

P⁽²⁾_m,l+j,n＝α₃(P_m+1,l+j,n+P_m,l+j,n+1+P_m-1,l+j,n+P_m,l+j,n-1)

+α₄(P_m+1,l+j,n+1+P_m+1,l+j,n-1+P_m-1,l+j,n+1+P_m-1,l+j,n-1)

+(1-4α₃-4α₄)P_m,l+j,n,j＝1,0,-1,

P⁽³⁾_m,l,n+j＝α₅(P_m+1,l,n+j+P_m,l+1,n+j+P_m-1,l,n+j+P_m,l-1,n+j)

+α₆(P_m+1,l+1,n+j+P_m+1,l-1,n+j+P_m-1,l+1,n+j+P_m-1,l-1,n+j)

+(1-4α₅-4α₆)P_m,l,n+j,j＝1,0,-1,

$\left(\begin{array}{c}\Sigma ={\alpha }_7{P}_{m,l,n}+{\alpha }_8(P_{m,l+1,n}+P_{m,l,n+1}+P_{m,l-1,n}+P_{m,l,n-1}+P_{m+1,l,n}+P_{m-1,l,n})\\ +\frac{1}{12}(1-{\alpha }_7-6{\alpha }_8)(P_{m+1,l+1,n}+P_{m+1,l,n+1}+P_{m+1,l-1,n}+P_{m+1,l,n-1}+P_{m-1,l+1,n}+P_{m-1,l,n+1}\\ +P_{m-1,l-1,n}+P_{m-1,l,n-1}+P_{m,l+1,n+1}+P_{m,l-1,n+1}+P_{m,l+1,n-1}+P_{m,l-1,n-1}),\end{array}\right)$

P_l,m,n≈P((m-1)Δx,(l-1)Δy,(n-1)Δz),m＝1,...,M；l＝1,...,L；n＝1,...,N,

其中，P表示压力波场；压力波场P的下标m+j,l+j,n+j表示离散点的三个不同指标；压力波场P的上标(1),(2),(3)分别表示对第一下标m+j、第二下标l+j、第三个下标n+j的离散点处的压力波场进行平均；s为Laplace阻尼常数；v_m,l,n表示三个指标分别为m,l,n的离散点处的波速；Δx,Δy,Δz表示离散点在网格方向x,y,z上的间隔；M,L,N分别表示离散点的三个指标的最大值，M,L,N为正整数；α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈为自由参数。

一个实施例中，所述自由参数确定单元，包括：数值速度获取模块，用于执行：将设定波场代入所述包含多个自由参数的27点差分方程，计算得到包含所述多个自由参数的波场衰减传播数值速度，所述设定波场的拟波长为波振幅衰减到初始波振幅的1/e^2π时波所经过的空间距离；目标函数建立模块，用于执行：利用所述波场衰减传播数值速度，建立包含所述多个自由参数的频散误差最小化目标函数；自由参数确定模块，用于执行：利用所述频散误差最小化目标函数确定使频散误差最小的所述自由参数的值。

一个实施例中，所述数值速度获取模块，还用于执行：归一化的所述波场衰减传播数值速度为：

$\frac{V_{num}(\theta ,\phi ,\tilde{k};{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6,{\alpha }_7,{\alpha }_8)}{v}=\frac{G}{2\pi }{\left(\frac{ϵ}{F}\right)}^{\frac{1}{2}},$

其中，

$\left(\begin{array}{c}ϵ=2[\cosh \left(\frac{2\pi sin\theta cos\phi }{G}\right)-1]E_x+2g_1^2[\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)-1]E_y\\ +2g_2^2[\cosh \left(\frac{2\pi \cos \theta }{g_{2}G}\right)-1]E_z,\end{array}\right)$

$G=\frac{2\pi }{{\mathrm{k\Delta x}}^{\prime }},{\mathrm{\Delta x}}^{\prime }=max\{\Delta x,\Delta y,\Delta z\},\tilde{k}=\frac{1}{G},g_1=\frac{\Delta x}{\Delta y},g_2=\frac{\Delta x}{\Delta z},$

$\left(\begin{array}{c}{E}_x=2{\alpha }_1\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{g_{1}G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)+4{\alpha }_2\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +(1-4{\alpha }_1-4{\alpha }_2),\end{array}\right)$

$\left(\begin{array}{c}{E}_y=2{\alpha }_3\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)+4{\alpha }_4\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +(1-4{\alpha }_3-4{\alpha }_4),\end{array}\right)$

$\left(\begin{array}{c}{E}_z=2{\alpha }_5\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{G})+\cosh (\frac{2\pi sin\theta sin\phi }{g_{1}G}\left)\right)+4{\alpha }_6\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\\ +(1-4{\alpha }_5-4{\alpha }_6),\end{array}\right)$

$\left(\begin{array}{c}F={\alpha }_7+2{\alpha }_8\left(\cosh \right(\frac{2\pi sin\theta cos\phi }{G})+\cosh (\frac{2\pi sin\theta cos\phi }{g_{1}G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)\\ +\frac{1-{\alpha }_7-6{\alpha }_8}{12}\left(\begin{array}{c}\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)+\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\end{array}\right),\end{array}\right)$

其中，θ为传播角；φ为方位角；v为离散点上的波速；α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈为自由参数；为归一化的波场衰减传播数值速度；k为拟波长对应的拟波数；Δx,Δy,Δz表示离散点在网格方向x,y,z上的间隔；

所述频散误差最小化目标函数为：

$E=\int \int \int {[1-\frac{V_{num}(\theta ,\phi ,\tilde{k};{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6,{\alpha }_7,{\alpha }_8)}{v}]}^{2}d\tilde{k}d\theta d\phi ,$

其中，E为衰减传播数值速度频散误差。

一个实施例中，所述数值模拟单元，包括：线性方程组建立模块，用于执行：根据离散速度模型的大小，利用确定所述多个自由参数的值后的27点差分方程建立离散波场线性方程组；离散点波场获取模块，用于执行：将阻尼常数、源子波、所述离散速度模型的网格间隔及所述离散速度模型上各离散点上的波速输入所述离散波场线性方程组，并利用所述离散速度模型的吸收边界条件，依据所述离散波场线性方程组中阻抗矩阵的大小选择求解方法求解所述离散波场线性方程组，得到波场在所述离散速度模型的各离散点上的值。

一个实施例中，所述数值模拟单元，还包括：震波图输出模块，用于执行：根据所述离散速度模型的各离散点上的值和检波器在所述离散速度模型中的位置输出所述检波器检测的震波图。

本发明的三维Laplace域声波方程数值模拟方法及装置，利用平均导数技巧成功实现构造带自由参数的27点差分方法，适应于任意的网格方向间隔比，通过使用自由参数，可以得到优化方法。进一步，本发明通过定义拟波长，波场衰减传播数值速度除以波场衰减传播的实际速度就能够得到归一化的波场衰减传播速度。归一化的波场衰减传播速度与自由参数相关，与Laplace参数、具体的速度、具体的采样间隔等具体量无关，所以是三维Laplace域频散分析的一般方法，克服了现有技术中频散关系依赖太多具体量而导致很难得到一般性的结论的问题。本发明优化的27点差分方法可以包含8个自由参数，再利用频散误差最小化目标函数得到自由参数的最优值，以此进行数值模拟，可使差分方法在百分之一的频散误差范围内，每一个最小拟波长只需要4个采样点，与现有方法相比，每个拟波长内的采样点从16个减少到4个，在同样模拟精度下，所需要的总的采样点仅为现有方法的1/64，计算效率至少提高64倍，数值计算效率极大提高，进而可提高三维Laplace域声波方程全波形反演的效率，大大减少了机时和能耗。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：

图1是本发明一实施例的三维Laplace域声波方程数值模拟方法的流程示意图；

图2是本发明一实施例中确定自由参数值的方法流程示意图；

图3至图5分别是利用现有27点差分方法在1/G的值分别为0.1、0.2及0.251时的频散相对误差结果；

图6至图8分别是利用本发明实施例的27点差分方法在1/G的值分别为0.1、0.2及0.251时的频散相对误差结果；

图9是本发明一实施例中利用27点差分方程进行数值模拟的方法流程示意图；

图10是本发明另一实施例中利用27点差分方程进行数值模拟的方法流程示意图；

图11是本发明另一实施例的三维Laplace域声波方程数值模拟方法的流程示意图；

图12是利用本发明实施例的27点差分方法、解析方法及现有27点差分方法进行数值模拟得到的震波图；

图13是本发明一实施例中三维Laplace域声波方程数值模拟装置的结构示意图；

图14是本发明一实施例中自由参数确定单元的结构示意图；

图15是本发明一实施例中数值模拟单元的结构示意图；

图16是本发明另一实施例中数值模拟单元的结构示意图。

具体实施方式

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

针对现有三维Laplace域声波方程全波形反演方法计算量大、效率低的问题，发明人考虑到要提高反演的效率，关键是提高三维Laplace域声波方程数值模拟的效率。发明人进一步分析现有问题的原因，发现现有三维Laplace域声波方程数值模方法是使用常规的27点差分方法，需要比较细的网格剖分来保持计算的精度，所以导致数目巨大的网格点数，计算量非常大。当应用于三维Laplace域声波方程全波形反演时，很容易导致计算效率低下，阻碍反演方法的应用。发展高效率的三维Laplace-Fourier域声波正演算法存在困难，主要是因为缺乏Laplace域正演算法的数值分析理论。首先遇到的困难是如何发展Laplace域数值分析方法以及Laplace参数所对应的拟波长如何定义的问题。现有的Laplace域数值分析的方法，将Laplace域频散关系表示为数值特征值和解析特征值之比。但这种频散关系依赖于Laplace参数、具体的速度、具体的采样间隔以及传播角和方位角，由于依赖于太多的具体的量，利用这种频散关系，很难得到一般性的结论，也无法对算法进行优化。基于上述一系列创造性思维过程，本发明提出一种优化的三维Laplace域声波方程数值模拟方法，能够解决现有技术遇到的问题。

图1是本发明一实施例的三维Laplace域声波方程数值模拟方法的流程示意图。如图1所示，本发明实施例的三维Laplace域声波方程数值模拟方法，可包括步骤：

S110：基于三维Laplace域声波方程，利用平均导数方法建立包含多个自由参数的27点差分方程；

S120：利用所述包含多个自由参数的27点差分方程，通过Laplace域频散分析确定所述多个自由参数的值；

S130：利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图。

在上述步骤S110中，该三维Laplace域声波方程可为：

$\frac{{\partial }^{2}P}{\partial x^2}+\frac{{\partial }^{2}P}{\partial y^2}+\frac{{\partial }^{2}P}{\partial z^2}-\frac{s^2}{v^2}P=f,---\left(1\right)$

其中，P可表示为P(x,y,z；s)，P表示压力波场；s表示Laplace阻尼常数，s为实数；v可表示为v(x,y,z)，v表示压力波场P的传播速度；f为源项，可表示源子波。

利用数值方法将微分方程(1)离散，使之变为差分方程。这一步是决定三维Laplace域声波方程数值模拟效率的关键。在上述步骤S110和步骤S120中，给出了一种优化的27点差分方法。具体地，在上述步骤S110中，基于三维Laplace域声波方程，利用平均导数方法建立包含多个自由参数的27点差分方程，得到了带有自由参数的差分方法，通过使用自由参数，可以得到优化方法。在上述步骤S120中，通过Laplace域频散分析确定上述多个自由参数的值，可以得到数值模拟所需的高效率的数值格式。

其中，三维Laplace域声波方程包含三个方向，利用平均导数方法建立包含多个自由参数的27点差分方程具体实施方式可为：在对一个方向上的导数进行离散时，同时利用其他两个方向上的波场的平均值，建立包含多个自由参数的27点差分方程。

在上述步骤S130中，利用优化后的27点差分方程进行数值模拟，可以降低计算量，显著提高数值模拟的效率。该种数值模拟方法用于全波形反演可提高大大提高反演效率。

本发明实施例，由于构建了带有自由参数的27点差分方程(格式)，可以对方法进行优化，适用于任意的方向网格间隔，应用范围广，可以提高数值模拟效率高。

一个实施例中，上述包含多个自由参数的27点差分方程可为：

$\begin{array}{c}\frac{{P^{\left(1\right)}}_{m+1,l,n}-2{P^{\left(1\right)}}_{m,l,n}+{P^{\left(1\right)}}_{m-1,l,n}}{{\mathrm{\Delta x}}^2}+\frac{{P^{\left(2\right)}}_{m,l+1,n}-2{P^{\left(2\right)}}_{m,l,n}+{P^{\left(2\right)}}_{m,l-1,n}}{{\mathrm{\Delta y}}^2}\\ +\frac{{P^{\left(3\right)}}_{m,l,n+1}-2{P^{\left(3\right)}}_{m,l,n}+{P^{\left(3\right)}}_{m,l,n-1}}{{\mathrm{\Delta z}}^2}-\frac{s^2}{v_{m,l,n}^2}\Sigma =0,\end{array}---\left(2\right)$

其中，

P⁽¹⁾_m+j,l,n＝α₁(P_m+j,l+1,n+P_m+j,l,n+1+P_m+j,l-1,n+P_m+j,l,n-1)

+α₂(P_m+j,l+1,n+1+P_m+j,l-1,n+1+P_m+j,l+1,n-1+P_m+j,l-1,n-1)

+(1-4α₁-4α₂)P_m+j,l,n,j＝1,0,-1,

P⁽²⁾_m,l+j,n＝α₃(P_m+1,l+j,n+P_m,l+j,n+1+P_m-1,l+j,n+P_m,l+j,n-1)

+α₄(P_m+1,l+j,n+1+P_m+1,l+j,n-1+P_m-1,l+j,n+1+P_m-1,l+j,n-1)

+(1-4α₃-4α₄)P_m,l+j,n,j＝1,0,-1,

P⁽³⁾_m,l,n+j＝α₅(P_m+1,l,n+j+P_m,l+1,n+j+P_m-1,l,n+j+P_m,l-1,n+j)

+α₆(P_m+1,l+1,n+j+P_m+1,l-1,n+j+P_m-1,l+1,n+j+P_m-1,l-1,n+j)

+(1-4α₅-4α₆)P_m,l,n+j,j＝1,0,-1,

$\left(\begin{array}{c}\Sigma ={\alpha }_7{P}_{m,l,n}+{\alpha }_8(P_{m,l+1,n}+P_{m,l,n+1}+P_{m,l-1,n}+P_{m,l,n-1}+P_{m+1,l,n}+P_{m-1,l,n})\\ +\frac{1}{12}(1-{\alpha }_7-6{\alpha }_8)(P_{m+1,l+1,n}+P_{m+1,l,n+1}+P_{m+1,l-1,n}+P_{m+1,l,n-1}+P_{m-1,l+1,n}+P_{m-1,l,n+1}\\ +P_{m-1,l-1,n}+P_{m-1,l,n-1}+P_{m,l+1,n+1}+P_{m,l-1,n+1}+P_{m,l+1,n-1}+P_{m,l-1,n-1}),\end{array}\right)$

P_l,m,n≈P((m-1)Δx,(l-1)Δy,(n-1)Δz),m＝1,...,M；l＝1,...,L；n＝1,...,N,

其中，P表示压力波场；压力波场P的下标m+j,l+j,n+j表示离散点的三个不同指标；压力波场P的上标(1),(2),(3)分别表示对第一下标m+j、第二下标l+j、第三个下标n+j的离散点处的压力波场进行平均；s为Laplace阻尼常数；v_m,l,n表示三个指标分别为m,l,n的离散点处的波速；Δx,Δy,Δz表示离散点在网格方向x,y,z上的间隔；M,L,N分别表示离散点的三个指标的最大值，M,L,N为正整数；α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈为自由参数。具体符号的意义，本领域技术人员能够结合上面标号的说明得以明确。

图2是本发明一实施例中确定自由参数值的方法流程示意图。如图2所示，在上述步骤S120中，利用所述包含多个自由参数的27点差分方程，通过Laplace域频散分析确定所述多个自由参数的值的方法，可包括步骤：

S121：将设定波场代入所述包含多个自由参数的27点差分方程，计算得到包含所述多个自由参数的波场衰减传播数值速度，所述设定波场的拟波长为波振幅衰减到初始波振幅的1/e^2π时波所经过的空间距离；

S122：利用所述波场衰减传播数值速度，建立包含所述多个自由参数的频散误差最小化目标函数；

S123：利用所述频散误差最小化目标函数确定使频散误差最小的所述自由参数的值。

为了确定这自由参数的值，以得到满足要求的优化算法，在上述步骤S121～S123中，建立了三维Laplace域声波方程数值分析的一般方法。具体地，在上述步骤S121中，定义了拟波长λ为波的振幅衰减到原来的1/e^2π时所经过的空间距离，而拟波数为k＝2π/λ，根据拟波长λ可以得到拟波数，进而根据拟波数k可以得到波场衰减传播数值速度。如此一来，定义了拟波长λ为波的振幅衰减到原来的1/e^2π时所经过的空间距离，可以解决现有技中不知如何对Laplace参数所对应的拟波长进行定义的问题。波场衰减传播数值速度除以波场衰减传播的实际速度就能够得到归一化的波场衰减传播速度。归一化的波场衰减传播速度与自由参数相关，与Laplace参数、具体的速度、具体的采样间隔等具体量无关，所以是三维Laplace域频散分析的一般方法，克服了现有技术中频散关系依赖太多具体量而导致很难得到一般性的结论的问题。在上述步骤S122和步骤S123中，利用所述波场衰减传播数值速度，建立包含所述多个自由参数的频散误差最小化目标函数，利用最优化理论来确定自由参数的最优值，可以针对不同的网格方向间隔比得到自由参数的最优值，从而可使27点差分方法的频散误差保持较小情况小，减少一个拟波长内所需要的采样点，从而可以在保持同等数值计算精确的情况下，无需要求太细的网格，所以能够减少数值模拟计算量，提高模拟效率。

一个实施例中，波场P(k,x,y,z)可为：

P(k,x,y,z)＝P₀e^{-k(sinθcosφx+sinθsinφy+cosθz)},(3)

其中，P₀为波场振幅，θ为传播角，而φ为方位角，k为拟波数，x,y,z为网格的三个方向。将定义的拟波长得到的拟波数代入公式(3)可以得到具体波场P(k,x,y,z)表达式。

一个实施例中，归一化的所述波场衰减传播数值速度可为：

$\frac{V_{num}(\theta ,\phi ,\tilde{k};{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6,{\alpha }_7,{\alpha }_8)}{v}=\frac{G}{2\pi }{\left(\frac{ϵ}{F}\right)}^{\frac{1}{2}},---\left(4\right)$

其中，

$\left(\begin{array}{c}ϵ=2[\cosh \left(\frac{2\pi sin\theta cos\phi }{G}\right)-1]E_x+2g_1^2[\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)-1]E_y\\ +2g_2^2[\cosh \left(\frac{2\pi \cos \theta }{g_{2}G}\right)-1]E_z,\end{array}\right)$

$G=\frac{2\pi }{{\mathrm{k\Delta x}}^{\prime }},{\mathrm{\Delta x}}^{\prime }=max\{\Delta x,\Delta y,\Delta z\},\tilde{k}=\frac{1}{G},g_1=\frac{\Delta x}{\Delta y},g_2=\frac{\Delta x}{\Delta z},$

$\left(\begin{array}{c}{E}_x=2{\alpha }_1\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{g_{1}G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)+4{\alpha }_2\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +(1-4{\alpha }_1-4{\alpha }_2),\end{array}\right)$

$\left(\begin{array}{c}{E}_y=2{\alpha }_3\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)+4{\alpha }_4\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +(1-4{\alpha }_3-4{\alpha }_4),\end{array}\right)$

$\left(\begin{array}{c}{E}_z=2{\alpha }_5\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{G})+\cosh (\frac{2\pi sin\theta sin\phi }{g_{1}G}\left)\right)+4{\alpha }_6\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\\ +(1-4{\alpha }_5-4{\alpha }_6),\end{array}\right)$

$\left(\begin{array}{c}F={\alpha }_7+2{\alpha }_8\left(\cosh \right(\frac{2\pi sin\theta cos\phi }{G})+\cosh (\frac{2\pi sin\theta cos\phi }{g_{1}G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)\\ +\frac{1-{\alpha }_7-6{\alpha }_8}{12}\left(\begin{array}{c}\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)+\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\end{array}\right),\end{array}\right)$

其中，θ为传播角；φ为方位角；v为离散点上的波速；α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈为自由参数；为归一化的波场衰减传播数值速度；k为拟波长对应的拟波数；Δx,Δy,Δz表示离散点在网格方向x,y,z上的间隔；该公式(4)可以通过将公式(3)代入公式(2)得到。

一个实施例中，所述频散误差最小化目标函数可为：

$E=\int \int \int {[1-\frac{V_{num}(\theta ,\phi ,\tilde{k};{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6,{\alpha }_7,{\alpha }_8)}{v}]}^{2}d\tilde{k}d\theta d\phi ,---\left(5\right)$

其中，E为衰减传播数值速度频散误差。

将公式(4)代入公式(5)可以得到具体的频散误差最小化目标函数。除自由参数外，公式(5)还可包括网格方向间隔比。

一个具体实施例中，27点差分方程包含8个自由参数：α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈，该8个自由参数对于不同的网格方向间隔比的最优值如表1所示。图3至图5分别是利用现有27点差分方法在1/G的值分别为0.1、0.2及0.251时的频散相对误差结果。图6至图8分别是利用本发明实施例的27点差分方法在1/G的值分别为0.1、0.2及0.251时的频散相对误差结果。对比图3和图6，图4和图7，图5和图8，即对比现有27点差分方法的频散误差和本发明实施例的27点差分方法与的频散误差对可知，随着拟波数(k＝1/G)的增加，现有27点差分方法的频散相对误差逐渐增加到18％，而本发明实施例的27点差分方法的频散相对误差始终保持在1％以下。由此可知，在百分之一的频散误差范围内，对于现有27点差分方法，每一个最小拟波长需要16个采样点，而对于本发明实施例的27点差分方法，每一个最小拟波长只需要4个采样点。

表1自由参数对于不同的网格方向间隔比的最优值

图9是本发明一实施例中利用27点差分方程进行数值模拟的方法流程示意图。如图9所示，在上述步骤S130中，利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图的方法，可包括步骤：

S131：根据离散速度模型的大小，利用确定所述多个自由参数的值后的27点差分方程建立离散波场线性方程组；

S132：将阻尼常数、源子波、所述离散速度模型的网格间隔及所述离散速度模型上各离散点上的波速输入所述离散波场线性方程组，并利用所述离散速度模型的吸收边界条件，依据所述离散波场线性方程组中阻抗矩阵的大小选择求解方法求解所述离散波场线性方程组，得到波场在所述离散速度模型的各离散点上的值。

上述27点差分方程(例如公式(2))可以表示定义在q＝M×L×N个离散点上的差分方程。在上述步骤S131，根据离散速度模型的大小lx,ly,lz(三个方向的距离)，确定离散点的个数q，将q个离散点处的27点差分方程组合起来形成一个线性方程组。在上述步骤S132中，将未知参数(阻尼常数、源子波、所述离散速度模型的网格间隔及所述离散速度模型上各离散点上的波速)代入离散波场线性方程组，并利用吸收边界条件求解该离散波场线性方程组，可以得到波场P在所述离散速度模型的各离散点上的值。其中，依据离散波场线性方程组中阻抗矩阵的大小可以选择不同的线性方程组求解方法，例如，当阻抗矩阵较小时，可采用直接方法求解，当阻抗矩阵较大时，可采用共轭梯度迭代法求解。

一个实施例中，上述离散波场线性方程组可为：

其中，S为一个q×q矩阵，表示阻抗矩阵，为波场在离散点上的值构成的列向量，F为源向量。

求解线性方程组公式(6)，可得到波场在离散点上的值根据阻抗矩阵S的大小，可采用不同的方法。当阻抗矩阵S较小时，可采用直接方法求解；当阻抗矩阵S为大型矩阵时，可采用共轭梯度法迭代求解。选择合适的求解方法有利于减少计算量提高并计算结果准确度。

图10是本发明另一实施例中利用27点差分方程进行数值模拟的方法流程示意图。如图10所示，图9所示的利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图的方法，还可包括步骤：

S133：根据所述离散速度模型的各离散点上的值和检波器在所述离散速度模型中的位置输出所述检波器检测的震波图。

在上述步骤S133中，波场在离散点上的值构成的列向量可包含波场P在所有离散点上的值，而震波图为波场P在检波器所在离散点上的值。给定检波器所在的离散点的位置，输出震波图。

图11是本发明另一实施例的三维Laplace域声波方程数值模拟方法的流程示意图。本发明实施例的三维Laplace域声波方程数值模拟的优化方法，可对方程(1)进行数值求解，得到波场P在检波器位置上的值，即震波图。如图11所示，本市实施例的方法可包括如下步骤：

第一步，输入离散速度模型参数，包括：离散模型的大小(在三个方向的距离)、网格方向间隔和在网格点上的速度值，并设定Laplace阻尼常数s和源f。

第二步，利用数值方法将微分方程公式(1)离散，使之变为差分方程，一方面，构造带有自由系数的差分方法，另一方面，建立三维Laplace域频散分析的一般方法，利用最优化理论来确定自由参数的值。

第二步是关键，决定了三维Laplace域声波方程数值模拟的效率，本发明实施例的优化的27点差分方法，求解自由参数后可以得到满足要求的数值格式。在该步骤中，在对一个方向上的导数进行离散时，同时利用其他两个方向上的波场的平均值(平均导数技巧)，可以得到方程公式(2)的带有8个自由系数的一个27点差分方法。为了确定这8个自由系数的值，从而得到满足我们要求的优化算法，可建立三维Laplace域声波方程数值分析的一般方法，定义拟波长λ为波的振幅衰减到原来的1/e^2π时所经过的空间距离，而k＝2π/λ为拟波数，代入公式(3)可以得到波场。

第三步，公式(2)表示的差分方程定义在q个离散点上，将这q个差分方程组合起来，并利用吸收边界条件，形成一个线性方程组；其中，q＝M×L×N。

第四步，求解线性方程组公式(6)，得到波场在离散点上的值根据阻抗矩阵S的大小，当S较小时，采用直接方法求解，当S为大型矩阵时，采用共轭梯度法迭代求解。

第五步，波场在离散点上的值包含波场P在所有离散点上的值，利用给定检波器所在的离散点的位置，输出震波图。

为了将本发明实施例的方法和解析解对比，采用均匀速度模型，模型大小可为4km*4km*2km。方向网格间隔可为Δx＝Δy＝100m，Δz＝50m，速度v＝2000m/s。源f可为一个主频为25Hz的Ricker子波，位于模型的中心。Laplace阻尼常数s可取作10π。检波器可位于地表500m处。图12是利用本发明实施例的27点差分方法、解析方法及现有27点差分方法进行数值模拟得到的震波图。本发明实施例的27点差分方法的结果103相比于现有27点差分方法的结果102更加接近解析解结果101。现有27点差分方法如果要达到同样的精度，必须将采样密度增加到原来的64倍，从而使得计算效率大大降低。换句话说，如果要达到同样的计算精度，本发明实施例的27点差分方法可以使用少得多的采样点，从而实现计算效率的大幅提高。对于复杂的模型，由于没有解析解，可以采用高阶的时间域格式和Laplace变化的方法来说明优化的27点差分方法的高效率。

本发明实施例的三维Laplace域声波方程数值模拟方法，在对一个方向上的导数进行离散时，同时利用其他两个方向上的波场的平均值(平均导数技巧)，成功实现构造带自由参数的27点差分方法，适应于任意的网格方向间隔比；进一步，本发明通过定义拟波长，发展了基于波场衰减传播速度的三维Laplace域数值频散分析方法，结合最优化理论确定自由参数的值，提高了数值模拟效率。本发明优化的27点差分方法包含8个自由参数，利用频散误差最小化目标函数得到自由参数的最优值，以此进行数值模拟，可使差分方法在百分之一的频散误差范围内，每一个最小拟波长只需要4个采样点，与现有方法相比，每个拟波长内的采样点从16个减少到4个，在同样模拟精度下，所需要的总的采样点仅为现有方法的1/64，计算效率至少提高64倍，数值计算效率极大提高，进而可提高三维Laplace域声波方程全波形反演的效率，大大减少了机时和能耗。

基于与图1所示的三维Laplace域声波方程数值模拟方法相同的发明构思，本申请实施例还提供了一种三维Laplace域声波方程数值模拟装置，如下面实施例所述。由于该三维Laplace域声波方程数值模拟装置解决问题的原理与三维Laplace域声波方程数值模拟方法相似，因此该三维Laplace域声波方程数值模拟装置的实施可以参见三维Laplace域声波方程数值模拟方法的实施，重复之处不再赘述。

图13是本发明一实施例中三维Laplace域声波方程数值模拟装置的结构示意图。如图13所示，本发明实施例中三维Laplace域声波方程数值模拟装置的结构示意图，可包括：差分方程建立单元210、自由参数确定单元220及数值模拟单元230，上述各单元顺序连接。

差分方程建立单元210用于执行：基于三维Laplace域声波方程，利用平均导数方法建立包含多个自由参数的27点差分方程。

自由参数确定单元220用于执行：利用所述包含多个自由参数的27点差分方程，通过Laplace域频散分析确定所述多个自由参数的值。

数值模拟单元230用于执行：利用确定所述多个自由参数的值后的27点差分方程进行数值模拟，输出震波图。

其中，三维Laplace域声波方程包含三个方向，利用平均导数方法建立包含多个自由参数的27点差分方程具体实施方式可为：在对一个方向上的导数进行离散时，同时利用其他两个方向上的波场的平均值，建立包含多个自由参数的27点差分方程。

在上述差分方程建立单元210中，该三维Laplace域声波方程可为：

$\frac{{\partial }^{2}P}{\partial x^2}+\frac{{\partial }^{2}P}{\partial y^2}+\frac{{\partial }^{2}P}{\partial z^2}-\frac{s^2}{v^2}P=f,---\left(1\right)$

其中，P可表示为P(x,y,z；s)，P表示压力波场；s表示Laplace阻尼常数，s为实数；v可表示为v(x,y,z)，v表示压力波场P的传播速度；f为源项，可表示源子波。

利用数值方法将微分方程(1)离散，使之变为差分方程。这一步是决定三维Laplace域声波方程数值模拟效率的关键。在差分方程建立单元210和自由参数确定单元220中，给出了一种优化的27点差分方法。具体地，在差分方程建立单元210中，基于三维Laplace域声波方程，利用平均导数方法建立包含多个自由参数的27点差分方程，得到了带有自由系数的差分方法，通过使用自由参数，可以得到优化方法。在自由参数确定单元220中，通过Laplace域频散分析确定上述多个自由参数的最优值，可以得到数值模拟所需的高效率数值格式。

在上述数值模拟单元230中，利用优化后的27点差分方程进行数值模拟，可以降低计算量，显著提高数值模拟的效率。该种数值模拟方法用于全波形反演可提高大大提高反演效率。

本发明实施例，由于构建了带有自由参数的27点差分方程(格式)，可以对方法进行优化，适用于任意的方向网格间隔，应用范围广，可以提高数值模拟效率高。

一个实施例中，所述差分方程建立单元210还用于执行：所述包含多个自由参数的27点差分方程为：

$\begin{array}{c}\frac{{P^{\left(1\right)}}_{m+1,l,n}-2{P^{\left(1\right)}}_{m,l,n}+{P^{\left(1\right)}}_{m-1,l,n}}{{\mathrm{\Delta x}}^2}+\frac{{P^{\left(2\right)}}_{m,l+1,n}-2{P^{\left(2\right)}}_{m,l,n}+{P^{\left(2\right)}}_{m,l-1,n}}{{\mathrm{\Delta y}}^2}\\ +\frac{{P^{\left(3\right)}}_{m,l,n+1}-2{P^{\left(3\right)}}_{m,l,n}+{P^{\left(3\right)}}_{m,l,n-1}}{{\mathrm{\Delta z}}^2}-\frac{s^2}{v_{m,l,n}^2}\Sigma =0,\end{array}---\left(2\right)$

其中，

P⁽¹⁾_m+j,l,n＝α₁(P_m+j,l+1,n+P_m+j,l,n+1+P_m+j,l-1,n+P_m+j,l,n-1)

+α₂(P_m+j,l+1,n+1+P_m+j,l-1,n+1+P_m+j,l+1,n-1+P_m+j,l-1,n-1)

+(1-4α₁-4α₂)P_m+j,l,n,j＝1,0,-1,

P⁽²⁾_m,l+j,n＝α₃(P_m+1,l+j,n+P_m,l+j,n+1+P_m-1,l+j,n+P_m,l+j,n-1)

+α₄(P_m+1,l+j,n+1+P_m+1,l+j,n-1+P_m-1,l+j,n+1+P_m-1,l+j,n-1)

+(1-4α₃-4α₄)P_m,l+j,n,j＝1,0,-1,

P⁽³⁾_m,l,n+j＝α₅(P_m+1,l,n+j+P_m,l+1,n+j+P_m-1,l,n+j+P_m,l-1,n+j)

+α₆(P_m+1,l+1,n+j+P_m+1,l-1,n+j+P_m-1,l+1,n+j+P_m-1,l-1,n+j)

+(1-4α₅-4α₆)P_m,l,n+j,j＝1,0,-1,

$\left(\begin{array}{c}\Sigma ={\alpha }_7{P}_{m,l,n}+{\alpha }_8(P_{m,l+1,n}+P_{m,l,n+1}+P_{m,l-1,n}+P_{m,l,n-1}+P_{m+1,l,n}+P_{m-1,l,n})\\ +\frac{1}{12}(1-{\alpha }_7-6{\alpha }_8)(P_{m+1,l+1,n}+P_{m+1,l,n+1}+P_{m+1,l-1,n}+P_{m+1,l,n-1}+P_{m-1,l+1,n}+P_{m-1,l,n+1}\\ +P_{m-1,l-1,n}+P_{m-1,l,n-1}+P_{m,l+1,n+1}+P_{m,l-1,n+1}+P_{m,l+1,n-1}+P_{m,l-1,n-1}),\end{array}\right)$

P_l,m,n≈P((m-1)Δx,(l-1)Δy,(n-1)Δz),m＝1,...,M；l＝1,...,L；n＝1,...,N,

其中，P表示压力波场；压力波场P的下标m+j,l+j,n+j表示离散点的三个不同指标；压力波场P的上标(1),(2),(3)分别表示对第一下标m+j、第二下标l+j、第三个下标n+j的离散点处的压力波场进行平均；s为Laplace阻尼常数；v_m,l,n表示三个指标分别为m,l,n的离散点处的波速；Δx,Δy,Δz表示离散点在网格方向x,y,z上的间隔；M,L,N分别表示离散点的三个指标的最大值，M,L,N为正整数；α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈为自由参数。

图14是本发明一实施例中自由参数确定单元的结构示意图。如图14所示，所述自由参数确定单元220，可包括：数值速度获取模块221、目标函数建立模块222及自由参数确定模块223，上述各模块顺序连接。

数值速度获取模块221用于执行：将设定波场代入所述包含多个自由参数的27点差分方程，计算得到包含所述多个自由参数的波场衰减传播数值速度，所述设定波场的拟波长为波振幅衰减到初始波振幅的1/e^2π时波所经过的空间距离。

目标函数建立模块222用于执行：利用所述波场衰减传播数值速度，建立包含所述多个自由参数的频散误差最小化目标函数。

自由参数确定模块223用于执行：利用所述频散误差最小化目标函数确定使频散误差最小的所述自由参数的值。

为了确定这自由参数的值，以得到满足要求的优化算法，在上述数值速度获取模块221、目标函数建立模块222及自由参数确定模块223中，建立了三维Laplace域声波方程数值分析的一般方法。具体地，在上述数值速度获取模块221中，定义了拟波长λ为波的振幅衰减到原来的1/e^2π时所经过的空间距离，而拟波数为k＝2π/λ，根据拟波长λ可以得到拟波数，进而根据拟波数k可以得到波场衰减传播数值速度。如此一来，定义了拟波长λ为波的振幅衰减到原来的1/e^2π时所经过的空间距离，可以解决现有技中不知如何对Laplace参数所对应的拟波长进行定义的问题。波场衰减传播数值速度波场衰减传播数值速度除以波场衰减传播的实际速度就得到归一化的波场衰减传播速度。归一化的波场衰减传播速度与自由参数相关，与Laplace参数、具体的速度、具体的采样间隔等具体量无关，所以是三维Laplace域频散分析的一般方法，克服了现有技术中频散关系依赖太多具体量而导致很难得到一般性的结论的问题。在上述目标函数建立模块222及自由参数确定模块223中，利用所述波场衰减传播数值速度，建立包含所述多个自由参数的频散误差最小化目标函数，利用最优化理论来确定自由参数的值，可以针对不同的网格方向间隔比得到自由参数的最优值，从而可使27点差分方法的频散误差在保持较小情况下，减少一个拟波长内所可需要的采样点，从而可以在保持同等数值计算精确的情况下，不需要太细的网格，所以能够减少数值模拟计算量，提高模拟效率。

一个实施例中，波场P(k,x,y,z)可为：

P(k,x,y,z)＝P₀e^{-k(sinθcosφx+sinθsinφy+cosθz)},>

其中，P₀为波场振幅，θ为传播角，而φ为方位角，k为拟波数，x,y,z为网格的三个方向。将定义的拟波长得到的拟波数代入公式(3)可以得到具体波场P(k,x,y,z)表达式。

一个实施例中，所述数值速度获取模块221还用于执行：归一化的所述波场衰减传播数值速度可为：

$\frac{V_{num}(\theta ,\phi ,\tilde{k};{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6,{\alpha }_7,{\alpha }_8)}{v}=\frac{G}{2\pi }{\left(\frac{ϵ}{F}\right)}^{\frac{1}{2}},---\left(4\right)$

其中，

$\left(\begin{array}{c}ϵ=2[\cosh \left(\frac{2\pi sin\theta cos\phi }{G}\right)-1]E_x+2g_1^2[\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)-1]E_y\\ +2g_2^2[\cosh \left(\frac{2\pi \cos \theta }{g_{2}G}\right)-1]E_z,\end{array}\right)$

$G=\frac{2\pi }{{\mathrm{k\Delta x}}^{\prime }},{\mathrm{\Delta x}}^{\prime }=max\{\Delta x,\Delta y,\Delta z\},\tilde{k}=\frac{1}{G},g_1=\frac{\Delta x}{\Delta y},g_2=\frac{\Delta x}{\Delta z},$

$\left(\begin{array}{c}{E}_x=2{\alpha }_1\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{g_{1}G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)+4{\alpha }_2\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +(1-4{\alpha }_1-4{\alpha }_2),\end{array}\right)$

$\left(\begin{array}{c}{E}_y=2{\alpha }_3\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)+4{\alpha }_4\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +(1-4{\alpha }_3-4{\alpha }_4),\end{array}\right)$

$\left(\begin{array}{c}{E}_z=2{\alpha }_5\left(\cosh \right(\frac{2\pi sin\theta sin\phi }{G})+\cosh (\frac{2\pi sin\theta sin\phi }{g_{1}G}\left)\right)+4{\alpha }_6\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\\ +(1-4{\alpha }_5-4{\alpha }_6),\end{array}\right)$

$\left(\begin{array}{c}F={\alpha }_7+2{\alpha }_8\left(\cosh \right(\frac{2\pi sin\theta cos\phi }{G})+\cosh (\frac{2\pi sin\theta cos\phi }{g_{1}G})+\cosh (\frac{2\pi cos\theta }{g_{2}G}\left)\right)\\ +\frac{1-{\alpha }_7-6{\alpha }_8}{12}\left(\begin{array}{c}\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)+\cosh \left(\frac{2\pi sin\theta sin\phi }{G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\\ +\cosh \left(\frac{2\pi sin\theta sin\phi }{g_{1}G}\right)\cosh \left(\frac{2\pi cos\theta }{g_{2}G}\right)\end{array}\right),\end{array}\right)$

其中，θ为传播角；φ为方位角；v为离散点上的波速；α₁,α₂,α₃,α₄,α₅,α₆,α₇,α₈为自由参数；为归一化的波场衰减传播数值速度；k为拟波长对应的拟波数；Δx,Δy,Δz表示离散点在网格方向x,y,z上的间隔；该公式(4)可以通过将公式(3)代入公式(2)得到。

一个实施例中，所述频散误差最小化目标函数可为：

$E=\int \int \int {[1-\frac{V_{num}(\theta ,\phi ,\tilde{k};{\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6,{\alpha }_7,{\alpha }_8)}{v}]}^{2}d\tilde{k}d\theta d\phi ,---\left(5\right)$

其中，E为衰减传播数值速度频散误差。将公式(4)代入公式(5)可以得到具体的频散误差最小化目标函数。除自由参数外，公式(5)还包括网格方向间隔比。

图15是本发明一实施例中数值模拟单元的结构示意图。如图15所示，所述数值模拟单元230，可包括：线性方程组建立模块231和离散点波场获取模块232，二者相互连接。

线性方程组建立模块231用于执行：根据离散速度模型的大小，利用确定所述多个自由参数的值后的27点差分方程建立离散波场线性方程组。

离散点波场获取模块232用于执行：将阻尼常数、源子波、所述离散速度模型的网格间隔及所述离散速度模型上各离散点上的波速输入所述离散波场线性方程组，并利用所述离散速度模型的吸收边界条件，依据所述离散波场线性方程组中阻抗矩阵的大小选择求解方法求解所述离散波场线性方程组，得到波场在所述离散速度模型的各离散点上的值。

上述27点差分方程(例如公式(2))可以表示定义在q＝M×L×N个离散点上的差分方程。在上述线性方程组建立模块231，根据离散速度模型的大小lx,ly,lz(三个方向的距离)，确定离散点的个数q，将q个离散点处的27点差分方程组合起来形成一个线性方程组。在上述离散点波场获取模块232中，将未知参数(阻尼常数、源子波、所述离散速度模型的网格间隔及所述离散速度模型上各离散点上的波速)代入离散波场线性方程组，并利用吸收边界条件求解该离散波场线性方程组，可以得到波场P在所述离散速度模型的各离散点上的值。其中，依据离散波场线性方程组中阻抗矩阵的大小可以选择不同的线性方程组求解方法，例如，当阻抗矩阵较小时，可采用直接方法求解，当阻抗矩阵较大时，可采用共轭梯度迭代法求解。

图16是本发明另一实施例中数值模拟单元的结构示意图。如图16所示，图15所示的数值模拟单元230，还可包括：震波图输出模块233，与离散点波场获取模块232连接。

震波图输出模块233用于执行：根据所述离散速度模型的各离散点上的值和检波器在所述离散速度模型中的位置输出所述检波器检测的震波图。

在上述震波图输出模块233中，波场在离散点上的值构成的列向量可包含波场P在所有离散点上的值，而震波图为波场P在检波器所在离散点上的值。给定检波器所在的离散点的位置，输出震波图。

本发明实施例的三维Laplace域声波方程数值模拟装置，在对一个方向上的导数进行离散时，同时利用其他两个方向上的波场的平均值(平均导数技巧)，成功实现构造带自由参数的27点差分方法，适应于任意的网格方向间隔比；进一步，本发明通过定义拟波长，发展了基于波场衰减传播速度的三维Laplace域数值频散分析方法，结合最优化理论确定自由参数的最优值，提高了数值模拟效率。本发明优化的27点差分方法包含8个自由参数，利用频散误差最小化目标函数得到自由参数的最优值，以此进行数值模拟，可使差分方法在百分之一的频散误差范围内，每一个最小拟波长只需要4个采样点，与现有方法相比，每个拟波长内的采样点从16个减少到4个，在同样模拟精度下，所需要的总的采样点仅为现有方法的1/64，计算效率至少提高64倍，数值计算效率极大提高，进而可提高三维Laplace域声波方程全波形反演的效率，大大减少了机时和能耗。

在本说明书的描述中，参考术语“一个实施例”、“一个具体实施例”、“一些实施例”、“例如”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。各实施例中涉及的步骤顺序用于示意性说明本发明的实施，其中的步骤顺序不作限定，可根据需要作适当调整。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。