CN103630933B - 基于非线性优化的时空域交错网格有限差分方法和装置 - Google Patents

Abstract

本发明提供了一种基于非线性优化的时空域交错网格有限差分方法和装置，其中，该方法包括：确定有限差分系数；基于时空域频散关系和非线性反演算法对有限差分系数进行优化；利用优化后的有限差分系数进行弹性波正演模拟。本发明解决了现有技术中采用泰勒级数展开和空间域频散关系的有限差分法获得有限差分系数进行弹性波正演模拟而导致的中高频段频散较大，模拟精度较低的技术问题，达到了减小中高频段的频散，提高模拟精度的技术效果。

Description

基于非线性优化的时空域交错网格有限差分方法和装置

技术领域

本发明涉及正演模拟技术领域，特别涉及一种基于非线性优化的时空域交错网格有限差分方法和装置。

背景技术

地震数值模拟技术就是对特定的地质、地球物理问题作适当的简化，以形成简化的数学模型，然后采用数值计算方法获取地震响应的过程。地震数值模拟技术是理解地震波在地下传播特点，帮助解释观测数据的有效手段。地震数值模拟还可以为新技术的提出、可行性分析和应用试验提供高质量的模拟数据；帮助地球物理工作者测试新的算法和处理技术，为地震反演问题提供思路和有效的验证数据。近年来，波动方程数值模拟方法被广泛应用于逆时偏移和全波形反演中。

常用的数值模拟方法主要包括：有限元素法、有限差分法和伪谱法等。其中，有限差分法是偏微分方程的主要数值解法之一，也是最早出现的数值模拟方法，其主要优点是物理意义直观，易于实现，能够较精确地模拟任意非均匀介质中的地震波场。有限差分法根据不同的标准可以分为：显式有限差分、隐式有限差分、规则网格有限差分、交错网格有限差分和旋转交错网格有限差分。在有限差分法中，差分系数可以通过泰勒级数展开或最优化方法求得，分别对应以泰勒级数展开为基础的有限差分和以最优化为基础的有限差分。

然而，使用基于泰勒级数展开和空间域频散关系的有限差分法存在如下问题：在低频段，频散接近于零，然而在中高频段，频散较大，从而导致模拟精度较低。

发明内容

本发明实施例提供了一种基于非线性优化的时空域交错网格有限差分方法，以达到减小中高频段的频散，提高模拟精度的目的，该方法包括：

确定有限差分系数；

基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

利用优化后的有限差分系数进行弹性波正演模拟。

在一个实施例中，所述确定有限差分系数，包括：

按照以下公式确定有限差分系数：

$a_m=\frac{{(-1)}^{m+1}}{2m-1}\underset{1\≤n\≤M,n\≠m}{\mathrm{\Π}}\left|\frac{{(2n-1)}^2}{{(2n-1)}^2-{(2m-1)}^2}\right|$

其中，a_m为有限差分系数，M为算子长度，m为有限差分系数的序号，1≤m≤M，n为连乘的变量。

在一个实施例中，基于时空域频散关系和非线性反演算法对有限差分系数进行优化，包括：

将有限差分系数作为初值确定P波和S波的频散大小；

根据确定的P波和S波的频散计算共轭梯度矢量；

根据所述共轭梯度矢量迭代对所述有限差分系数进行优化。

在一个实施例中，将有限差分系数作为初值确定的P波和S波的频散大小为：

其中， $A={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\cos \mathrm{\θ}\right)\right)}^2,B={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\sin \mathrm{\θ}\right)\right)}^2,$ 为P波的频散大小，为S波的频散大小，M为算子长度，κ=kh，k为波数，h为空间采样间隔，r_P=v_PΔt/h，r_S=v_SΔt/h，v_P为P波的传播速度，v_S为S波的传播速度，Δt为时间采样间隔，θ为波传播方向角度，a_m为有限差分系数；

根据P波和S波的频散大小计算共轭梯度矢量包括：

确定目标函数：

其中， $b1=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_P}h,\mathrm{\π}),b2=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_S}h,\mathrm{\π}),$ f_max为最大频率；

根据所述目标函数计算共轭梯度矢量：

$p_{k+1}=-\▿E_{k+1}+\frac{\▿E_{k+1}^T\▿E_{k+1}}{\▿E_k^T\▿E_k}{p}_k,p_0=-\▿E_0$

其中，p_k表示当前时刻的共轭梯度矢量， $\▿E={\left[\begin{array}{cccc}\frac{\∂E}{\∂a_1}& \frac{\∂E}{\∂a_2}& \·\·\·& \frac{\∂E}{\∂a_M}\end{array}\right]}^T,$ k表示当前时刻，k+1表示当前时刻的下一时刻；

按照以下公式，根据所述共轭梯度矢量迭代对有限差分系数进行优化：

a_k+1=a_k+α_kp_k+1,a=[_a1a₂…a_M]^T，其中，a为不同时刻优化后的有限差分系数向量，a₁,a₂…a_M为优化后的有限差分系数，α_k为迭代步长。

在一个实施例中，基于时空域频散关系和非线性反演算法对有限差分系数进行优化，包括：

基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

对优化后的有限差分系数进行校验；

如果校验结果不满足约束条件，则改变算子长度，再根据改变后的算子长度对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件；

或者，如果校验结果不满足约束条件，则改变最大频率值，再根据改变后的最大频率值对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件。

在一个实施例中，利用优化后的有限差分系数进行弹性波正演模拟，包括：

将优化后的有限差分系数代入二维弹性介质速度应力方程以实现弹性波正演模拟。

在一个实施例中，所述二维弹性介质速度应力方程为：

$\mathrm{\ρ}\frac{v_{\mathrm{xi},j}^{n+1}-v_{\mathrm{xi},j}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xxi}+m-1/2,j}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xxi}-m+1/2,j}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{xzi},j+m-1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi},j-m+1/2}^{n+1/2})\end{array}\right]$

$\mathrm{\ρ}\frac{v_{\mathrm{zi}+1/2,j+1/2}^{n+1}-v_{\mathrm{zi}+1/2,j+1/2}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xzi}+m,j+1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi}-m+1,j+1/2}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{zzi}+1/2,j+m}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j-m+1}^{n+1/2})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +\mathrm{\λ}(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\λ}(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\μ}(v_{\mathrm{zi}+m-1/2,j}^{n+1}-v_{\mathrm{zi}-m+1/2,j}^{n+1})\\ +\mathrm{\μ}(v_{\mathrm{xi},j+m}^{n+1}-v_{\mathrm{xi},j-m+1}^{n+1})\end{array}\right]$

其中，(v_x,v_z)为速度矢量，(τ_xx,τ_zz,τ_xz)为应力的三个分量，λ(x,z)和μ(x,z)为拉梅系数，ρ(x,z)为密度，Δt为时间采样间隔，h为空间采样间隔，a_m为优化后的有限差分系数，M为算子长度，x为沿水平方向的坐标，z为沿铅垂方向的坐标，i为沿水平方向的采样序号，j为沿铅垂方向的采样序号，n为沿时间方向的采样序号，m为有限差分系数的序号，1≤m≤M。

本发明实施例还提供了一种基于非线性优化的时空域交错网格有限差分装置，以达到减小中高频段的频散，提高模拟精度的目的，该装置包括：

确定模块，用于确定有限差分系数；

优化模块，用于基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

模拟模块，用于利用优化后的有限差分系数进行弹性波正演模拟。

在一个实施例中，所述确定模块具体用于按照以下公式确定有限差分系数：

$a_m=\frac{{(-1)}^{m+1}}{2m-1}\underset{1\≤n\≤M,n\≠m}{\mathrm{\Π}}\left|\frac{{(2n-1)}^2}{{(2n-1)}^2-{(2m-1)}^2}\right|$

其中，a_m为有限差分系数，M为算子长度，m为有限差分系数的序号，1≤m≤M，n为连乘的变量。

在一个实施例中，所述优化模块包括：

确定单元，用于将有限差分系数作为初值确定P波和S波的频散大小；

计算单元，用于根据确定的P波和S波的频散计算共轭梯度矢量；

优化单元，用于根据所述共轭梯度矢量迭代对所述有限差分系数进行优化。

在一个实施例中，所述确定单元具体用于按照以下公式确定P波和S波的频散大小：

其中， $A={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\cos \mathrm{\θ}\right)\right)}^2,B={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\sin \mathrm{\θ}\right)\right)}^2,$ 为P波的频散大小，为S波的频散大小，M为算子长度，κ=kh，k为波数，h为空间采样间隔，r_P=v_PΔt/h，r_S=v_SΔt/h，v_P为P波的传播速度，v_S为S波的传播速度，Δt为时间采样间隔，θ为波传播方向角度，a_m为有限差分系数；

所述计算单元具体按照以下方式计算共轭梯度矢量：

确定目标函数：

其中， $b1=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_P}h,\mathrm{\π}),b2=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_S}h,\mathrm{\π}),$ f_max为最大频率；

根据所述目标函数计算共轭梯度矢量：

$p_{k+1}=-\▿E_{k+1}+\frac{\▿E_{k+1}^T\▿E_{k+1}}{\▿E_k^T\▿E_k}{p}_k,p_0=-\▿E_0$

其中，p_k表示当前时刻的共轭梯度矢量， $\▿E={\left[\begin{array}{cccc}\frac{\∂E}{\∂a_1}& \frac{\∂E}{\∂a_2}& \·\·\·& \frac{\∂E}{\∂a_M}\end{array}\right]}^T,$ k表示当前时刻，k+1表示当前时刻的下一时刻；

所述优化单元具体用于按照以下公式对有限差分系数进行优化：

a_k+1=a_k+α_kp_k+1,a=[a₁a₂…a_M]^T，其中，a为不同时刻优化后的有限差分系数向量，a₁,a₂…a_M为优化后的有限差分系数，α_k为迭代步长。

在一个实施例中，优化模块包括：

系数优化单元，用于基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

校验单元，用于对优化后的有限差分系数进行校验，如果校验结果不满足约束条件，则改变算子长度，再根据改变后的算子长度对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件；或者，如果校验结果不满足约束条件，则改变最大频率值，再根据改变后的最大频率值对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件。

在一个实施例中，所述模拟模块具体用于将优化后的有限差分系数代入二维弹性介质速度应力方程以实现弹性波正演模拟。

在一个实施例中，所述二维弹性介质速度应力方程为：

$\mathrm{\ρ}\frac{v_{\mathrm{xi},j}^{n+1}-v_{\mathrm{xi},j}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xxi}+m-1/2,j}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xxi}-m+1/2,j}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{xzi},j+m-1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi},j-m+1/2}^{n+1/2})\end{array}\right]$

$\mathrm{\ρ}\frac{v_{\mathrm{zi}+1/2,j+1/2}^{n+1}-v_{\mathrm{zi}+1/2,j+1/2}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xzi}+m,j+1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi}-m+1,j+1/2}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{zzi}+1/2,j+m}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j-m+1}^{n+1/2})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +\mathrm{\λ}(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\λ}(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\μ}(v_{\mathrm{zi}+m-1/2,j}^{n+1}-v_{\mathrm{zi}-m+1/2,j}^{n+1})\\ +\mathrm{\μ}(v_{\mathrm{xi},j+m}^{n+1}-v_{\mathrm{xi},j-m+1}^{n+1})\end{array}\right]$

其中，(v_x,v_z)为速度矢量，(τ_xx,τ_zz,τ_xz)为应力的三个分量，λ(x,z)和μ(x,z)为拉梅系数，ρ(x,z)为密度，Δt为时间采样间隔，h为空间采样间隔，a_m为优化后的有限差分系数，M为算子长度，x为沿水平方向的坐标，z为沿铅垂方向的坐标，i为沿水平方向的采样序号，j为沿铅垂方向的采样序号，n为沿时间方向的采样序号，m为有限差分系数的序号，1≤m≤M。

在本发明实施例中，通过时空域频散关系和非线性反演算法对有限差分系数进行优化，并利用优化后的有限差分系数进行弹性波正演模拟，从而解决了现有技术中采用泰勒级数展开和空间域频散关系的有限差分法获得有限差分系数进行弹性波正演模拟而导致的中高频段频散较大，模拟精度较低的技术问题，达到了减小中高频段的频散，提高模拟精度的技术效果。

附图说明

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

图1是本发明实施例的基于非线性优化的时空域交错网格有限差分方法的流程图；

图2是本发明实施例的进行有限差分系数优化的方法流程图；

图3是本发明实施例的基于非线性优化的时空域交错网格有限差分装置的结构框图；

图4是本发明实施例的基于非线性优化的时空域交错网格有限差分方法的具体流程图；

图5是常规方法确定的不同传播角度时δ随频率的变化曲线示意图；

图6是本发明实施例确定的不同传播角度时δ随频率的变化曲线示意图；

图7是常规方法确定的不同算子长度时δ随频率的变化曲线；

图8本发明实施例确定的不同算子长度时δ随频率的变化曲线。

具体实施方式

发明人发现，在常规的有限差分方法中，差分系数一般是通过极小化空间域的频散关系得到的，这种方式在低频段，频散误差接近于零，然而在中高频段，频散较大，不能很好地描述地震波在时空域传播的规律。为此，在本实施例中提出了一种基于非线性优化的时空域交错网格有限差分方法来进行弹性波正演模拟，该方法通过极小化时间域和空间域的频散关系，以及非线性反演算法来求取最优的差分系数，以达到减小在中高频段的频散，提高模拟精度的目的。

在本发明实施例中，提出了一种基于非线性优化的时空域交错网格有限差分方法，如图1所示，包括以下步骤：

步骤101：确定有限差分系数；

步骤102：基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

步骤103：利用优化后的有限差分系数进行弹性波正演模拟。

在上述实施例中，通过时空域频散关系和非线性反演算法对有限差分系数进行优化，并利用优化后的有限差分系数进行弹性波正演模拟，从而解决了现有技术中采用泰勒级数展开和空间域频散关系的有限差分法获得有限差分系数进行弹性波正演模拟而导致的中高频段频散较大，模拟精度较低的技术问题，达到了减小中高频段的频散，提高模拟精度的技术效果。

在上述步骤101中，可以采用常规的有限差分系数求取方法获得有限差分系数，例如可以按照以下公式计算有限差分系数：

$a_m=\frac{{(-1)}^{m+1}}{2m-1}\underset{1\≤n\≤M,n\≠m}{\mathrm{\Π}}\left|\frac{{(2n-1)}^2}{{(2n-1)}^2-{(2m-1)}^2}\right|$

其中，a_m为有限差分系数，M为算子长度，m为有限差分系数的序号，1≤m≤M，n为连乘的变量。

具体的，上述步骤102基于时空域频散关系和非线性反演算法对有限差分系数进行优化可以如图2所示，包括以下步骤：

步骤201：将有限差分系数作为初值确定P波和S波的频散大小，其中，P波是指纵波（Primary/Compressional wave），S波是指横波（Second/Shear wave）：

其中， $A={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\cos \mathrm{\θ}\right)\right)}^2,B={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\sin \mathrm{\θ}\right)\right)}^2,$ 为P波的频散大小，为S波的频散大小，M为算子长度，κ=kh，k为波数，h为空间采样间隔，r_P=v_PΔt/h，r_S=v_SΔt/h，v_P为P波的传播速度，v_S为S波的传播速度，Δt为时间采样间隔，θ为波传播方向角度，a_m表示初始有限差分系数；

步骤202：根据确定的P波和S波的频散计算共轭梯度矢量；

先确定目标函数：

其中， $b1=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_P}h,\mathrm{\π}),b2=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_S}h,\mathrm{\π}),$ f_max为最大频率；

然后根据所述目标函数计算共轭梯度矢量：

$p_{k+1}=-\▿E_{k+1}+\frac{\▿E_{k+1}^T\▿E_{k+1}}{\▿E_k^T\▿E_k}{p}_k,p_0=-\▿E_0$

其中，p_k表示当前时刻的共轭梯度矢量， $\▿E={\left[\begin{array}{cccc}\frac{\∂E}{\∂a_1}& \frac{\∂E}{\∂a_2}& \·\·\·& \frac{\∂E}{\∂a_M}\end{array}\right]}^T,$ k表示当前时刻，k+1表示当前时刻的下一时刻；

步骤203：按照以下公式，根据所述共轭梯度矢量迭代对有限差分系数进行优化：

a_k+1=a_k+α_kp_k+1,a=[a₁a₂…a_M]^T，其中，a表示不同时刻优化后的有限差分系数向量，a₁,a₂…a_M为优化后的有限差分系数，α_k表示迭代步长。

考虑到算子长度M和最大频率在对有限差分系数的优化过程中起着重要的作用，可以通过对优化后的有限差分系数进行校验，判断其是否满足约束条件，如果不满足，可以更改算子长度M的值或者改变最大频率f_max的值，然后重新确定最佳有限差分系数。

具体的可以包括：

步骤1：基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

步骤2：对优化后的有限差分系数进行校验；

步骤3：如果校验结果不满足约束条件，则改变算子长度，再根据改变后的算子长度对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件；或者，如果校验结果不满足约束条件，则改变最大频率值，再根据改变后的最大频率值对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件。

在上述各个实施例中，上述步骤103利用优化后的有限差分系数进行弹性波正演模拟，可以包括：将优化后的有限差分系数代入二维弹性介质速度应力方程以实现弹性波正演模拟。

其中，该二维弹性介质速度应力方程为：

$\mathrm{\ρ}\frac{v_{\mathrm{xi},j}^{n+1}-v_{\mathrm{xi},j}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xxi}+m-1/2,j}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xxi}-m+1/2,j}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{xzi},j+m-1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi},j-m+1/2}^{n+1/2})\end{array}\right]$

$\mathrm{\ρ}\frac{v_{\mathrm{zi}+1/2,j+1/2}^{n+1}-v_{\mathrm{zi}+1/2,j+1/2}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xzi}+m,j+1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi}-m+1,j+1/2}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{zzi}+1/2,j+m}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j-m+1}^{n+1/2})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +\mathrm{\λ}(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\λ}(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\μ}(v_{\mathrm{zi}+m-1/2,j}^{n+1}-v_{\mathrm{zi}-m+1/2,j}^{n+1})\\ +\mathrm{\μ}(v_{\mathrm{xi},j+m}^{n+1}-v_{\mathrm{xi},j-m+1}^{n+1})\end{array}\right]$

其中，(v_x,v_z)为速度矢量，(τ_xx,τ_zz,τ_xz)为应力的三个分量，λ(x,z)和μ(x,z)为拉梅系数，ρ(x,z)为密度，Δt为时间采样间隔，h为空间采样间隔，a_m为优化的有限差分系数,M为算子长度，x为沿水平方向的坐标，z为沿铅垂方向的坐标，i为沿水平方向的采样序号，j为沿铅垂方向的采样序号，n为沿时间方向的采样序号，m为有限差分系数的序号，1≤m≤M。

基于同一发明构思，本发明实施例中还提供了一种基于非线性优化的时空域交错网格有限差分装置，如下面的实施例所述。由于基于非线性优化的时空域交错网格有限差分装置解决问题的原理与基于非线性优化的时空域交错网格有限差分方法相似，因此基于非线性优化的时空域交错网格有限差分装置的实施可以参见基于非线性优化的时空域交错网格有限差分方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。图3是本发明实施例的基于非线性优化的时空域交错网格有限差分装置的一种结构框图，如图3所示，包括：确定模块301、优化模块302和模拟模块303，下面对该结构进行具体说明。

确定模块301，用于确定有限差分系数；

优化模块302，用于基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

模拟模块303，用于利用优化后的有限差分系数进行弹性波正演模拟。

在一个实施例中，确定模块301具体用于按照以下公式确定有限差分系数：

$a_m=\frac{{(-1)}^{m+1}}{2m-1}\underset{1\≤n\≤M,n\≠m}{\mathrm{\Π}}\left|\frac{{(2n-1)}^2}{{(2n-1)}^2-{(2m-1)}^2}\right|$

其中，a_m为有限差分系数，M为算子长度，m为有限差分系数的序号，1≤m≤M，n为连乘的变量。

在一个实施例中，优化模块302包括：确定单元，用于将有限差分系数作为初值确定P波和S波的频散大小；计算单元，用于根据确定的P波和S波的频散计算共轭梯度矢量；优化单元，用于根据所述共轭梯度矢量迭代对所述有限差分系数进行优化。

在一个实施例中，确定单元具体用于按照以下公式确定P波和S波的频散大小：

其中， $A={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\cos \mathrm{\θ}\right)\right)}^2,B={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\sin \mathrm{\θ}\right)\right)}^2,$ 为P波的频散大小，为S波的频散大小，M为算子长度，κ=kh，k为波数，h为空间采样间隔，r_P=v_PΔt/h，r_S=v_SΔt/h，v_P为P波的传播速度，v_S为S波的传播速度，Δt为时间采样间隔，θ为波传播方向角度，a_m为有限差分系数；

计算单元具体按照以下方式计算共轭梯度矢量：确定目标函数：

其中， $b1=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_P}h,\mathrm{\π}),b2=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_S}h,\mathrm{\π}),$ f_max为最大频率；

根据所述目标函数计算共轭梯度矢量：

$p_{k+1}=-\▿E_{k+1}+\frac{\▿E_{k+1}^T\▿E_{k+1}}{\▿E_k^T\▿E_k}{p}_k,p_0=-\▿E_0$

其中，p_k表示当前时刻的共轭梯度矢量， $\▿E={\left[\begin{array}{cccc}\frac{\∂E}{\∂a_1}& \frac{\∂E}{\∂a_2}& \·\·\·& \frac{\∂E}{\∂a_M}\end{array}\right]}^T,$ k表示当前时刻，k+1表示当前时刻的下一时刻；

优化单元具体用于按照以下公式对有限差分系数进行优化：

a_k+1=a_k+α_kp_k+1,a=[a₁a₂…a_M]^T，其中，a为不同时刻优化后的有限差分系数向量，a₁,a₂…a_M为优化后的有限差分系数，α_k为迭代步长。

在一个实施例中，优化模块302包括：系数优化单元，用于基于时空域频散关系和非线性反演算法对有限差分系数进行优化；校验单元，用于对优化后的有限差分系数进行校验；如果校验结果不满足约束条件，则改变算子长度，再根据改变后的算子长度对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件；或者，如果校验结果不满足约束条件，则改变最大频率值，再根据改变后的最大频率值对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件。

在一个实施例中，模拟模块303具体用于将优化后的有限差分系数代入二维弹性介质速度应力方程以实现弹性波正演模拟。

在一个实施例中，二维弹性介质速度应力方程为：

$\mathrm{\ρ}\frac{v_{\mathrm{xi},j}^{n+1}-v_{\mathrm{xi},j}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xxi}+m-1/2,j}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xxi}-m+1/2,j}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{xzi},j+m-1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi},j-m+1/2}^{n+1/2})\end{array}\right]$

$\mathrm{\ρ}\frac{v_{\mathrm{zi}+1/2,j+1/2}^{n+1}-v_{\mathrm{zi}+1/2,j+1/2}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xzi}+m,j+1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi}-m+1,j+1/2}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{zzi}+1/2,j+m}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j-m+1}^{n+1/2})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +\mathrm{\λ}(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\λ}(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\μ}(v_{\mathrm{zi}+m-1/2,j}^{n+1}-v_{\mathrm{zi}-m+1/2,j}^{n+1})\\ +\mathrm{\μ}(v_{\mathrm{xi},j+m}^{n+1}-v_{\mathrm{xi},j-m+1}^{n+1})\end{array}\right]$

其中，(v_x,v_z)为速度矢量，(τ_xx,τ_zz,τ_xz)为应力的三个分量，λ(x,z)和μ(x,z)为拉梅系数，ρ(x,z)为密度，Δt为时间采样间隔，h为空间采样间隔，a_m为优化后的有限差分系数，M为算子长度，x为沿水平方向的坐标，z为沿铅垂方向的坐标，i为沿水平方向的采样序号，j为沿铅垂方向的采样序号，n为沿时间方向的采样序号，m为有限差分系数的序号，1≤m≤M。

为了对上述弹性波正演模拟方法进行更为清楚的解释，下面结合一个具体的实施例来进行说明，然而值得注意的是该实施例仅是为了更好地说明本发明，并不构成对本发明不当的限定。

本发明在分析了常规有限差分法的不足的情况下、以时空域频散关系和非线性反演算法作为基础上提出了一种基于非线性优化的时空域交错网格有限差分法，以减小中高频段数值频散的大小、提高弹性波方程正演模拟的精度，该方法适用于弹性波方程数值模拟。

如图4所示，是本发明实施例的基于非线性优化的时空域交错网格有限差分方法流程图，包括以下步骤：

步骤401：根据有关参数对计算区域进行网格剖分；

步骤402：求取常规有限差分法的差分系数：

$a_m=\frac{{(-1)}^{m+1}}{2m-1}\underset{1\≤n\≤M,n\≠m}{\mathrm{\Π}}\left|\frac{{(2n-1)}^2}{{(2n-1)}^2-{(2m-1)}^2}\right|$

步骤403：以常规有限差分系数为初值，基于时空域频散关系和非线性反演算法求取最佳有限差分系数（即优化后的有限差分系数）；

步骤404：将得到的差分系数代入二维弹性介质速度-应力方程，以实现弹性波正演模拟。

其中，二维弹性介质速度-应力方程表示如下：

$\mathrm{\ρ}\frac{\∂v_x}{\∂t}=\frac{\∂{\mathrm{\τ}}_{\mathrm{xx}}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂z}$

$\mathrm{\ρ}\frac{\∂v_z}{\∂t}=\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{\mathrm{zz}}}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{xx}}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_x}{\∂x}+\mathrm{\λ}\frac{\∂v_z}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{zz}}}{\∂t}=\mathrm{\λ}\frac{\∂v_x}{\∂x}+(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_z}{\∂z}$

$\frac{\∂{\mathrm{\τ}}_{\mathrm{xz}}}{\∂t}=\mathrm{\μ}(\frac{\∂v_z}{\∂x}+\frac{\∂v_x}{\∂z})$

其中，(v_x,v_z)为速度矢量，(τ_xx,τ_zz,τ_xz)为应力的三个分量，λ(x,z)和μ(x,z)为拉梅系数，ρ(x,z)为密度。

$\mathrm{\ρ}\frac{v_{\mathrm{xi},j}^{n+1}-v_{\mathrm{xi},j}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xxi}+m-1/2,j}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xxi}-m+1/2,j}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{xzi},j+m-1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi},j-m+1/2}^{n+1/2})\end{array}\right]$

$\mathrm{\ρ}\frac{v_{\mathrm{zi}+1/2,j+1/2}^{n+1}-v_{\mathrm{zi}+1/2,j+1/2}^n}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}({\mathrm{\τ}}_{\mathrm{xzi}+m,j+1/2}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{xzi}-m+1,j+1/2}^{n+1/2})\\ +({\mathrm{\τ}}_{\mathrm{zzi}+1/2,j+m}^{n+1/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j-m+1}^{n+1/2})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xxi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +\mathrm{\λ}(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{zzi}+1/2,j}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\λ}(v_{\mathrm{xi}+m,j}^{n+1}-v_{\mathrm{xi}-m+1,j}^{n+1})\\ +(\mathrm{\λ}+2\mathrm{\μ})(v_{\mathrm{zi}+1/2,j+m-1/2}^{n+1}-v_{\mathrm{zi}+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+3/2}-{\mathrm{\τ}}_{\mathrm{xzi},j+1/2}^{n+1/2}}{\mathrm{\Δt}}=\frac{1}{h}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\left[\begin{array}{c}\mathrm{\μ}(v_{\mathrm{zi}+m-1/2,j}^{n+1}-v_{\mathrm{zi}-m+1/2,j}^{n+1})\\ +\mathrm{\μ}(v_{\mathrm{xi},j+m}^{n+1}-v_{\mathrm{xi},j-m+1}^{n+1})\end{array}\right]$

其中，Δt是时间采样间隔，h为空间采样间隔，a_m为优化后的最佳有限差分系数,M为算子长度。

下面对图4所示的弹性波方程数值模拟流程进行具体描述。

弹性波方程交错网格模拟，时空域频散关系如下：

$A+B\≈r_P^{-2}{\sin }^2\left(0.5{\mathrm{\κr}}_P\right)$

$A+B\≈r_S^{-2}{\sin }^2\left(0.5{\mathrm{\κr}}_S\right)$

其中：

$A={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\cos \mathrm{\θ}\right)\right)}^2$

$B={\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m\sin \left((m-0.5)\mathrm{\κ}\sin \mathrm{\θ}\right)\right)}^2$

其中，κ=kh，k为波数，r_P=v_PΔt/h，r_S=v_SΔt/h，v_P和v_S分别为P波和S波的传播速度。

定义两个新函数：

其中，和为上述公式两边的相对误差，用于分别表示P波和S波的频散大小，如果和都为0，则表示没有频散，通过极小化和可以求得最佳的交错网格的有限差分系数。

定义目标函数为：

对于给定的频率范围（0≤f≤f_max)，P波和S波的κ是不同的，计算公式如下：

$b1=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_P}h,\mathrm{\π}),b2=\max (\frac{{2\mathrm{\πf}}_{\max }}{v_S}h,\mathrm{\π}),$

可以采用非线性共轭梯度法来求解该优化问题，共轭梯度矢量为：

$p_{k+1}=-\▿E_{k+1}+\frac{\▿E_{k+1}^T\▿E_{k+1}}{\▿E_k^T\▿E_k}{p}_k,p_0=-\▿E_0$

其中：

$\▿E={\left[\begin{array}{cccc}\frac{\∂E}{\∂a_1}& \frac{\∂E}{\∂a_2}& \·\·\·& \frac{\∂E}{\∂a_M}\end{array}\right]}^T,$ k和k+1分别表示当前时刻和下一时刻。

迭代公式为：

a_k+1=a_k+α_kp_k+1,a=[a₁a₂…a_M]^T，其中,α_k为迭代步长。

考虑到算子长度M在弹性波正演模拟方法中的最佳有限差分系数的确定过程中起着重要的作用，可以通过对求取的最佳有限差分系数进行校验，使其满足约束条件，如果不满足，可以更改算子长度M的值，然后重新确定最佳有限差分系数。具体的可以采用下面的公式来度量数值频散的大小：

${\mathrm{\δ}}_P\left(\mathrm{\θ}\right)=\frac{v_{\mathrm{PFD}}}{v_P}=\frac{2}{r_P\mathrm{\κ}}\arcsin \left(r_P\sqrt{A+B}\right)$

${\mathrm{\δ}}_S\left(\mathrm{\θ}\right)=\frac{v_{\mathrm{SFD}}}{v_S}=\frac{2}{r_P\mathrm{\κ}}\arcsin \left(r_S\sqrt{A+B}\right)$

其中，δ_P(θ)表示P波的相速度误差，δ_S(θ)表示S波的相速度误差，δ_P和δ_S越接近于1表示频散越小。

在本例中，定义以下公式来平衡纵波和横波的频散误差：

δ(f,θ,M)=0.5(δP-1)+0.5(δ_S-1)

该公式可以近似表示纵波和横波相速度的相对误差，δ越接近于0表示频散越小。

最大频散误差δ_max如下：

${\mathrm{\δ}}_{\max }(M,f_{\max })=\underset{f\∈[0,f_{\max }],\mathrm{\θ}\∈\left[\mathrm{0,2}\mathrm{\π}\right]}{\max }\left|\mathrm{\δ}(f,\mathrm{\θ},M)\right|$

为了保证本文方法的有效性，在前面的最优化过程中加入约束条件：

δ_max<η，其中，η为最大允许误差。

其中，δ_max只与f_max和M有关。当M固定时，δ_max决定于f_max，δ_max随着f_max增加而变大。因此，若初始的f_max不满足上述约束条件，则通过逐渐减小f_max直至δ_max<η来获得最佳f_max。另外，当f_max固定时，δ_max只与M有关，δ_max随着M增加而变小，因此如果初始的M不满足上述约束条件，则通过逐渐增大M直到δ_max<η来获得最佳M。

最终将确定的最佳有限差分系数代入到上述二维弹性介质速度-应力方程中就实现了弹性波数值模拟。

以一个均匀弹性介质为例来说明本发明实施例的优点，在本例中，纵波速度为2800m/s，横波速度1500m/s，时间采样间隔为1ms，空间采用间隔为h=20m，算子长度M=14。

图5和图6表示在不同差分算法被使用时数值频散随传播方向的变化规律，图7和图8表示在不同差分算法被使用时数值频散随算子长度的变化规律。由图5至图8可见，与常规交错网格有限差分相比，基于非线性最优化的时空域交错网格有限差分法具有更宽的有效频带及更小的数值频散。通过当算子长度相同时，基于非线性最优化的时空域交错网格有限差分法模拟的精度更高。

本实施例中减小了中高频段的数值频散，以时空域频散关系为基础更符合实际，模拟精度更高，以常规有限差分法的差分系数为初值，结合非线性共轭梯度算法可以在较少的迭代次数下得到满意的结果。该方法可以应用到所有涉及到数值模拟的地球物理研究过程中，例如：逆时偏移、全波形反演，可以大大减少中高频段数值的频散，提高弹性波方程正演模拟的精度。

在另外一个实施例中，还提供了一种软件，该软件用于执行上述实施例及优选实施方式中描述的技术方案。

在另外一个实施例中，还提供了一种存储介质，该存储介质中存储有上述软件，该存储介质包括但不限于：光盘、软盘、硬盘、可擦写存储器等。

从以上的描述中，可以看出，本发明实施例实现了如下技术效果：通过时空域频散关系和非线性反演算法对有限差分系数进行优化，并利用优化后的有限差分系数进行弹性波正演模拟，从而解决了现有技术中采用泰勒级数展开和空间域频散关系的有限差分法获得有限差分系数进行弹性波正演模拟而导致的中高频段频散较大，模拟精度较低的技术问题，达到了减小中高频段的频散，提高模拟精度的技术效果。

显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例不限制于任何特定的硬件和软件结合。

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

Claims (12)

1.一种基于非线性优化的时空域交错网格有限差分方法，其特征在于，包括：

确定有限差分系数；

基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

利用优化后的有限差分系数进行弹性波正演模拟；

其中，基于时空域频散关系和非线性反演算法对有限差分系数进行优化，包括：

将有限差分系数作为初值确定P波和S波的频散大小；

根据确定的P波和S波的频散计算共轭梯度矢量；

根据所述共轭梯度矢量迭代对所述有限差分系数进行优化。

2.如权利要求1所述的方法，其特征在于，所述确定有限差分系数，包括：

按照以下公式确定有限差分系数：

$a_m=\frac{{(-1)}^{m+1}}{2m-1}\underset{1\&le;n\&le;M,n\&NotEqual;m}{\&Pi;}\left|\frac{{(2n-1)}^2}{{(2n-1)}^2-{(2m-1)}^2}\right|$

其中，a_m为有限差分系数，M为算子长度，m为有限差分系数的序号，1≤m≤M，n为连乘的变量。

3.如权利要求1所述的方法，其特征在于：

将有限差分系数作为初值确定的P波和S波的频散大小为：

其中，为P波的频散大小，为S波的频散大小，M为算子长度，κ＝kh，k为波数，h为空间采样间隔，r_P＝v_PΔt/h，r_S＝v_SΔt/h，v_P为P波的传播速度，v_S为S波的传播速度，Δt为时间采样间隔，θ为波传播方向角度，a_m为有限差分系数，其中，m为有限差分系数的序号，1≤m≤M；

根据P波和S波的频散大小计算共轭梯度矢量包括：

确定目标函数：

其中，f_max为最大频率；

根据所述目标函数计算共轭梯度矢量：

$p_{k+1}=-\&dtri;E_{k+1}+\frac{\&dtri;E_{k+1}^T\&dtri;E_{k+1}}{\&dtri;E_k^T\&dtri;E_k}{p}_k,p_0=-\&dtri;E_0$

其中，p_k表示当前时刻的共轭梯度矢量，k表示当前时刻，k+1表示当前时刻的下一时刻；

按照以下公式，根据所述共轭梯度矢量迭代对有限差分系数进行优化：

a_k+1＝a_k+α_kp_k+1,a＝[a₁ a₂ … a_M]^T，其中，a为不同时刻优化后的有限差分系数向量，a₁,a₂,…,a_M为优化后的有限差分系数，α_k为迭代步长。

4.如权利要求3所述的方法，其特征在于，基于时空域频散关系和非线性反演算法对有限差分系数进行优化，包括：

基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

对优化后的有限差分系数进行校验；

如果校验结果不满足约束条件，则改变算子长度，再根据改变后的算子长度对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件；

或者，如果校验结果不满足约束条件，则改变最大频率值，再根据改变后的最大频率值对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件。

5.如权利要求1至4中任一项所述的方法，利用优化后的有限差分系数进行弹性波正演模拟，包括：

将优化后的有限差分系数代入二维弹性介质速度应力方程以实现弹性波正演模拟。

6.如权利要求5所述的方法，其特征在于，所述二维弹性介质速度应力方程为：

$\mathrm{\&rho;}\frac{{v_x}_{i,j}^{n+1}-{v_x}_{i,j}^n}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}({{\mathrm{\&tau;}}_{xx}}_{i+m-1/2,j}^{n+1/2}-{{\mathrm{\&tau;}}_{xx}}_{i-m+1/2,j}^{n+1/2})\\ +({{\mathrm{\&tau;}}_{xz}}_{j+m-1/2}^{n+1/2}-{{\mathrm{\&tau;}}_{xz}}_{i,j-m+1/2}^{n+1/2})\end{array}\right]$

$\mathrm{\&rho;}\frac{{v_z}_{i+1/2,j+1/2}^{n+1}-{v_z}_{i+1/2,j+1/2}^n}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}({{\mathrm{\&tau;}}_{xz}}_{i+m,j+1/2}^{n+1/2}-{{\mathrm{\&tau;}}_{xz}}_{i-m+1,j+1/2}^{n+1/2})\\ +({{\mathrm{\&tau;}}_{zz}}_{i+1/2,j+m}^{n+1/2}-{{\mathrm{\&tau;}}_{zz}}_{i+1/2,j-m+1}^{n+1/2})\end{array}\right]$

$\frac{{{\mathrm{\&tau;}}_{xx}}_{i+1/2,j}^{n+3/2}-{{\mathrm{\&tau;}}_{xx}}_{i+1/2,j}^{n+1/2}}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}(\mathrm{\&lambda;}+2\mathrm{\&mu;})({v_x}_{i+m,j}^{n+1}-{v_x}_{i-m+1,j}^{n+1})\\ +({v_z}_{i+1/2,j+m-1/2}^{n+1}-{v_z}_{i+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{{\mathrm{\&tau;}}_{zz}}_{i+1/2,j}^{n+3/2}-{{\mathrm{\&tau;}}_{zz}}_{i+1/2,j}^{n+1/2}}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}\mathrm{\&lambda;}({v_x}_{i+m,j}^{n+1}-{v_x}_{i-m+1,j}^{n+1})\\ +(\mathrm{\&lambda;}+2\mathrm{\&mu;})({v_z}_{i+1/2,j+m-1/2}^{n+1}-{v_z}_{i+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{{\mathrm{\&tau;}}_{xz}}_{i,j+1/2}^{n+3/2}-{{\mathrm{\&tau;}}_{xz}}_{i,j+1/2}^{n+1/2}}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}\mathrm{\&mu;}({v_z}_{i+m-1/2,j}^{n+1}-{v_z}_{i-m+1/2,j}^{n+1})\\ +\mathrm{\&mu;}({v_x}_{i,j+m}^{n+1}-{v_x}_{i,j-m+1}^{n+1})\end{array}\right]$

其中，(v_x,v_z)为速度矢量，(τ_xx,τ_zz,τ_xz)为应力的三个分量，λ(x,z)和μ(x,z)为拉梅系数，ρ(x,z)为密度，Δt为时间采样间隔，h为空间采样间隔，a_m为优化后的有限差分系数，M为算子长度，x为沿水平方向的坐标，z为沿铅垂方向的坐标，i为沿水平方向的采样序号，j为沿铅垂方向的采样序号，n为沿时间方向的采样序号，m为有限差分系数的序号，1≤m≤M。

7.一种基于非线性优化的时空域交错网格有限差分装置，其特征在于，包括：

确定模块，用于确定有限差分系数；

优化模块，用于基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

模拟模块，用于利用优化后的有限差分系数进行弹性波正演模拟；

其中，所述优化模块包括：

确定单元，用于将有限差分系数作为初值确定P波和S波的频散大小；

计算单元，用于根据确定的P波和S波的频散计算共轭梯度矢量；

优化单元，用于根据所述共轭梯度矢量迭代对所述有限差分系数进行优化。

8.如权利要求7所述的装置，其特征在于，所述确定模块具体用于按照以下公式确定有限差分系数：

$a_m=\frac{{(-1)}^{m+1}}{2m-1}\underset{1\&le;n\&le;M,n\&NotEqual;m}{\&Pi;}\left|\frac{{(2n-1)}^2}{{(2n-1)}^2-{(2m-1)}^2}\right|$

其中，a_m为有限差分系数，M为算子长度，m为有限差分系数的序号，1≤m≤M，n为连乘的变量。

9.如权利要求7所述的装置，其特征在于：

所述确定单元具体用于按照以下公式确定P波和S波的频散大小：

其中，为P波的频散大小，为S波的频散大小，M为算子长度，κ＝kh，k为波数，h为空间采样间隔，r_P＝v_PΔt/h，r_S＝v_SΔt/h，v_P为P波的传播速度，v_S为S波的传播速度，Δt为时间采样间隔，θ为波传播方向角度，a_m为有限差分系数，其中，m为有限差分系数的序号，1≤m≤M；

所述计算单元具体按照以下方式计算共轭梯度矢量：

确定目标函数：

其中，f_max为最大频率；

根据所述目标函数计算共轭梯度矢量：

$p_{k+1}=-\&dtri;E_{k+1}+\frac{\&dtri;E_{k+1}^T\&dtri;E_{k+1}}{\&dtri;E_k^T\&dtri;E_k}{p}_k,p_0=-\&dtri;E_0$

其中，p_k表示当前时刻的共轭梯度矢量，k表示当前时刻，k+1表示当前时刻的下一时刻；

所述优化单元具体用于按照以下公式对有限差分系数进行优化：

a_k+1＝a_k+α_kp_k+1,a＝[a₁ a₂ … a_M]^T，其中，a为不同时刻优化后的有限差分系数向量，a₁,a₂,…,a_M为优化后的有限差分系数，α_k为迭代步长。

10.如权利要求9所述的装置，其特征在于，所述优化模块包括：

系数优化单元，用于基于时空域频散关系和非线性反演算法对有限差分系数进行优化；

校验单元，用于对优化后的有限差分系数进行校验，如果校验结果不满足约束条件，则改变算子长度，再根据改变后的算子长度对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件；或者，如果校验结果不满足约束条件，则改变最大频率值，再根据改变后的最大频率值对有限差分系数进行优化，直至优化后的有限差分系数的校验结果可以满足所述约束条件。

11.如权利要求7至10中任一项所述的装置，所述模拟模块具体用于将优化后的有限差分系数代入二维弹性介质速度应力方程以实现弹性波正演模拟。

12.如权利要求11所述的装置，其特征在于，所述二维弹性介质速度应力方程为：

$\mathrm{\&rho;}\frac{{v_x}_{i,j}^{n+1}-{v_x}_{i,j}^n}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}({{\mathrm{\&tau;}}_{xx}}_{i+m-1/2,j}^{n+1/2}-{{\mathrm{\&tau;}}_{xx}}_{i-m+1/2,j}^{n+1/2})\\ +({{\mathrm{\&tau;}}_{xz}}_{j+m-1/2}^{n+1/2}-{{\mathrm{\&tau;}}_{xz}}_{i,j-m+1/2}^{n+1/2})\end{array}\right]$

$\mathrm{\&rho;}\frac{{v_z}_{i+1/2,j+1/2}^{n+1}-{v_z}_{i+1/2,j+1/2}^n}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}({{\mathrm{\&tau;}}_{xz}}_{i+m,j+1/2}^{n+1/2}-{{\mathrm{\&tau;}}_{xz}}_{i-m+1,j+1/2}^{n+1/2})\\ +({{\mathrm{\&tau;}}_{zz}}_{i+1/2,j+m}^{n+1/2}-{{\mathrm{\&tau;}}_{zz}}_{i+1/2,j-m+1}^{n+1/2})\end{array}\right]$

$\frac{{{\mathrm{\&tau;}}_{xx}}_{i+1/2,j}^{n+3/2}-{{\mathrm{\&tau;}}_{xx}}_{i+1/2,j}^{n+1/2}}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}(\mathrm{\&lambda;}+2\mathrm{\&mu;})({v_x}_{i+m,j}^{n+1}-{v_x}_{i-m+1,j}^{n+1})\\ +({v_z}_{i+1/2,j+m-1/2}^{n+1}-{v_z}_{i+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{{\mathrm{\&tau;}}_{zz}}_{i+1/2,j}^{n+3/2}-{{\mathrm{\&tau;}}_{zz}}_{i+1/2,j}^{n+1/2}}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}\mathrm{\&lambda;}({v_x}_{i+m,j}^{n+1}-{v_x}_{i-m+1,j}^{n+1})\\ +(\mathrm{\&lambda;}+2\mathrm{\&mu;})({v_z}_{i+1/2,j+m-1/2}^{n+1}-{v_z}_{i+1/2,j-m+1/2}^{n+1})\end{array}\right]$

$\frac{{{\mathrm{\&tau;}}_{xz}}_{i,j+1/2}^{n+3/2}-{{\mathrm{\&tau;}}_{xz}}_{i,j+1/2}^{n+1/2}}{\mathrm{\&Delta;}t}=\frac{1}{h}\underset{m=1}{\overset{M}{\&Sigma;}}{a}_m\left[\begin{array}{c}\mathrm{\&mu;}({v_z}_{i+m-1/2,j}^{n+1}-{v_z}_{i-m+1/2,j}^{n+1})\\ +\mathrm{\&mu;}({v_x}_{i,j+m}^{n+1}-{v_x}_{i,j-m+1}^{n+1})\end{array}\right]$

其中，(v_x,v_z)为速度矢量，(τ_xx,τ_zz,τ_xz)为应力的三个分量，λ(x,z)和μ(x,z)为拉梅系数，ρ(x,z)为密度，Δt为时间采样间隔，h为空间采样间隔，a_m为优化后的有限差分系数，M为算子长度，x为沿水平方向的坐标，z为沿铅垂方向的坐标，i为沿水平方向的采样序号，j为沿铅垂方向的采样序号，n为沿时间方向的采样序号，m为有限差分系数的序号，1≤m≤M。