一种基于弹性波解耦方程的优化拟解析方法及装置

摘要

本发明涉及一种基于弹性波解耦方程的优化拟解析方法及装置，其中，方法包括：获得等效的P波、S波解耦的弹性波方程，根据等效的P波、S波解耦的弹性波方程分别获得对应地P波、S波优化归一化拟拉普拉斯算子；根据P波、S波优化归一化拟拉普拉斯算子获得P波、S波优化拟微分算子的象征；对P波、S波优化拟微分算子的象征分别进行分解；将P波、S波优化拟微分算子的象征的分解结果代入到弹性波优化拟微分算子的表达式中，获得弹性波优化拟微分算子的最终表达式；利用弹性波优化微分算子的最终表达式分别对等效的P波、S波解耦的弹性波方程进行计算，获得P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。

权利要求

1.一种基于弹性波解耦方程的优化拟解析方法，其特征在于，包括：

对二阶位移常密度弹性波解耦方程进行分解，获得等效的P波、S波解耦的弹性波方程，并利用所述等效的P波、S波解耦的弹性波方程分别获得对应地P波、S波优化归一化拟拉普拉斯算子；

对P波、S波优化归一化拟拉普拉斯算子利用弹性波优化拟微分算子的表达式进行表示，获得P波、S波优化拟微分算子的象征；

对P波、S波优化拟微分算子的象征分别进行分解；

将P波、S波优化拟微分算子的象征的分解结果代入到弹性波优化拟微分算子的表达式中，获得弹性波优化拟微分算子的最终表达式；

利用所述弹性波优化微分算子的最终表达式分别对所述等效的P波、S波解耦的弹性波方程进行计算，获得P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。

2.如权利要求1所述的方法，其特征在于，所述二阶位移常密度弹性波解耦方程的表达式为：其中，x为空间矢量，x、y、z表示空间矢量x中的三个垂直方向的坐标；Vp(x)和Vs(x)分别为纵波速度与横波速度，u(x,t)、v(x,t)、w(x,t)分别为质点在混合波场作用下空间矢量x中的三个垂直方向上产生的位移分量。

3.如权利要求2所述的方法，其特征在于，所述解耦的弹性波方程中的P波方程的表达式为：式中，up(x,t)、vp(x,t)、wp(x,t)分别为质点在P波作用下空间矢量x中的三个垂直方向产生的位移分量；

所述解耦的弹性波方程中的S波方程的表达式为：

式中，us(x,t)、vs(x,t)、ws(x,t)分别为质点在S波作用下空间矢量x中的三个垂直方向产生的位移分量。

4.如权利要求3所述的方法，其特征在于，所述P波、S波优化归一化拟拉普拉斯算子的表达式为：式中，kela是空间波数矢量，其中，kp、ks分别为P波、S波空间波数矢量，ela＝p,s表示P波和S波； 为弹性波优化归一化拟拉普拉斯算子，其中， 为P波优化归一化拟拉普拉斯算子， 为S波优化归一化拟拉普拉斯算子；Δt是时间采样间隔。

5.如权利要求4所述的方法，其特征在于，所述弹性波优化拟微分算子的表达式为：其中，α(x,kela)称为弹性波优化拟微分算子A[E(x,t)]的象征，E(kela,t)为位移矢量E(x,t)的空间傅里叶变换。

6.如权利要求5所述的方法，其特征在于，所述P波、S波优化拟微分算子的象征的表达式为：式中，αopt,ela(x,kela)中ela＝p,s分别表示P波优化拟微分算子的象征和S波优化拟微分算子的象征。

7.如权利要求6所述的方法，其特征在于，所述弹性波优化拟微分算子的最终表达式为：式中，N表示空间矢量的最大个数；n表示空间矢量的个数编号；M表示空间波数矢量的最大个数；m表示空间波数矢量的个数编号；W(x,km)是象征αopt,ela(x,kela)的子矩阵,由与km有关的特定的列元素组成；W(xn,kela)是象征αopt,ela(x,kela)的另一个子矩阵，由与xn有关的特定行元素组成；αmn为系数矩阵。

8.如权利要求7所述的方法，其特征在于，对P波、S波优化拟微分算子的象征分别采用低秩分解近似方法进行分解。

9.如权利要求8所述的方法，其特征在于，所述P波、S波优化拟微分算子的象征近似分解之后的表达式为：其中，W(x,km)是象征αopt,ela(x,kela)的子矩阵,其由与km有关的特定的列元素组成；W(xn,kela)是象征αopt,ela(x,kela)的另一个子矩阵，其由与xn有关的特定行元素组成；αmn代表中间的系数矩阵。

10.一种基于弹性波解耦方程的优化拟解析装置，其特征在于，包括：弹性波优化归一化拟拉普拉斯算子计算单元，用于对二阶位移常密度弹性波解耦方程进行分解，获得等效的P波、S波解耦的弹性波方程，并利用所述等效的P波、S波解耦的弹性波方程分别获得对应地P波、S波优化归一化拟拉普拉斯算子；

弹性波优化拟微分算子的象征确定单元，用于对P波、S波优化归一化拟拉普拉斯算子利用弹性波优化拟微分算子表示，获得P波、S波优化拟微分算子的象征；

分解单元，用于对P波、S波优化拟微分算子的象征分别进行分解；

弹性波优化拟微分算子的最终表达式获取单元，用于将P波、S波优化拟微分算子象征的分解结果代入到弹性波优化拟微分算子表达式中，获得弹性波优化拟微分算子的最终表达式；

波场分量计算单元，用于利用所述弹性波优化微分算子的最终表达式分别对所述等效的P波、S波解耦的弹性波方程进行计算，获得P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。

说明书

一种基于弹性波解耦方程的优化拟解析方法及装置

技术领域

[0001] 本发明涉及地震波波场模拟技术领域，特别涉及一种基于弹性波解耦方程的优化拟解析。

背景技术

[0002] 地震波波场模拟(地震正演)和成像(逆时偏移)是勘探地球物理领域两个重要研究方向。高效高精度地震正演算子是地震成像和反演的基础和关键。常用的地震波波场模拟方法可以分为：有限差分法、伪谱法等。其中，常规有限差分正演算子具有精度低且不稳定等缺点，为了抑制数值频散，空间网格的大小通常比地震波长小8倍～10倍。而为了满足稳定性条件，时间采样间隔需取得较小。因此，为了达到地震波波场模拟精度，基于泰勒级数展开的有限差分方法，由于空间和时间采样不能太大，导致其计算量很大。对于常规伪谱法来说，在时间方向上二阶差分格式精度较低，且对于大步长时间采样间隔，伪谱法不稳定；由于伪谱法多次利用快速傅立叶正变换和反变换，其计算量相对于有限差分有所增加。

[0003] 拟解析法首先由Etgen和Brandsberg-Dahl在2009年提出，类似于Lax-Wendroff方法的补偿思想：利用空间导数项，近似补偿时间二阶导数离散产生的误差。Lax-Wendroff方法利用空间的高阶导数项来近似代替高阶时间差分，补偿时间二阶差分近似所造成的误差，从而提高了时间精度。拟解析法可以在波数域修正拉普拉斯算子。在波数-空间域，利用拟拉普拉斯算子，可以补偿因时间二阶差分而造成的解析误差，从而压制时间方向二阶差分格式造成的误差。在速度变化缓慢的情况下，拟解析法可以有效地消除时间方向的数值频散。对于速度变化不是太剧烈的模型，拟解析法可以通过级联几个恒定补偿速度的拟拉普拉斯算子，补偿由于二阶时间差分引起的误差。但当速度变化较为剧烈时，Etgen和Brandsberg-Dahl指出可采用多个常速下的拟拉普拉斯算子进行速度插值，近似变速的拟拉普拉斯算子。虽然采用速度插值策略(3个～5个参考补偿速度)对于各向同性介质可取得较好效果，但对于各向异性介质(VTI或TTI)效果较差，因拟拉普拉斯算子不仅包括速度，而且还包含各向异性等参数，其需要组合速度和各向异性参数对拟拉普拉斯算子进行插值，不仅计算量增大，且拟解析法在时间和空间上均会产生较大误差。

发明内容

[0004] 本发明实施例的主要目的在于提出一种基于弹性波解耦方程的优化拟解析方法及装置，将声波归一化拟拉普拉斯算子改进为弹性波优化归一化拟拉普拉斯算子，利用弹性波优化拟微分算子表示弹性波优化归一化拟拉普拉斯算子，该算子不仅包括原始拟微分算子的谱估计，而且还包含一个时间补偿项，将在时间频散和空间频散方面，优于传统的伪谱法和拟解析法。

[0005] 为实现上述目的，本发明实施例提供了一种基于弹性波解耦方程的优化拟解析方法，包括：

[0006] 对二阶位移常密度弹性波解耦方程进行分解，获得等效的P波、S波解耦的弹性波方程，并利用所述等效的P波、S波解耦的弹性波方程分别获得对应地P波、S波优化归一化拟拉普拉斯算子；

[0007] 对P波、S波优化归一化拟拉普拉斯算子利用弹性波优化拟微分算子的表达式进行表示，获得P波、S波优化拟微分算子的象征；

[0008] 对P波、S波优化拟微分算子的象征分别进行分解；

[0009] 将P波、S波优化拟微分算子象征的分解结果代入到弹性波优化拟微分算子的表达式中，获得弹性波优化拟微分算子的最终表达式；

[0010] 利用所述弹性波优化微分算子的最终表达式分别对所述等效的P波、S波解耦的弹性波方程进行计算，获得P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。

[0011] 可选的，在本发明一实施例中，所述二阶位移常密度弹性波解耦方程的表达式为：

[0012]

[0013] 其中，x为空间矢量，x、y、z表示空间矢量x中的三个垂直方向的坐标；Vp(x)和Vs(x)分别为纵波速度与横波速度，u(x,t)、v(x,t)、w(x,t)分别为质点在混合波场作用下在空间矢量x中的三个垂直方向上产生的位移分量。

[0014] 可选的，在本发明一实施例中，所述解耦的弹性波方程中的P波方程的表达式为：

[0015]

[0016] 式中，up(x,t)、vp(x,t)、wp(x,t)分别为质点在P波作用下空间矢量x中的三个垂直方向产生的位移分量；

[0017] 所述解耦的弹性波方程中的S波方程的表达式为：

[0018]

[0019] 式中，us(x,t)、vs(x,t)、ws(x,t)分别为质点在S波作用下空间矢量x中的三个垂直方向产生的位移分量。

[0020] 可选的，在本发明一实施例中，所述P波、S波优化归一化拟拉普拉斯算子的表达式为：

[0021]

[0022] 式中，kela是空间波数矢量，其中，kp、ks分别为P波、S波空间波数矢量，ela＝p,s表示P波和S波； 为弹性波优化归一化拟拉普拉斯算子，其中， 为P波优化归一化拟拉普拉斯算子， 为S波优化归一化拟拉普拉斯算子；Δt是时间采样间隔。

[0023] 可选的，在本发明一实施例中，所述弹性波优化拟微分算子的表达式为：

[0024]

[0025] 其中，α(x,kela)称为弹性波优化拟微分算子A[E(x,t)]的象征，E(kela,t)为位移矢量E(x,t)的空间傅里叶变换。

[0026] 可选的，在本发明一实施例中，所述P波、S波优化拟微分算子的象征的表达式为：

[0027]

[0028] 式中，αopt,ela(x,kela)中ela＝p,s分别表示P波优化拟微分算子的象征和S波优化拟微分算子的象征。

[0029] 可选的，在本发明一实施例中，所述弹性波优化拟微分算子的最终表达式为：

[0030]

[0031] 式中，N表示空间矢量的最大个数；n表示空间矢量的个数编号；M表示空间波数矢量的最大个数；m表示空间波数矢量的个数编号。

[0032] 可选的，在本发明一实施例中，对P波、S波优化拟微分算子的象征分别采用低秩分解近似方法进行分解。

[0033] 可选的，在本发明一实施例中，所述P波、S波优化拟微分算子的象征近似分解之后的表达式为：

[0034]

[0035] 其中，W(x,km)是象征αopt,ela(x,kela)的子矩阵,其由与km有关的特定的列元素组成；W(xn,kela)是象征αopt,ela(x,kela)的另一个子矩阵，其由与xn有关的特定行元素组成；αmn代表中间的系数矩阵。

[0036] 对应地，为实现上述目的，本发明实施例提供了一种基于弹性波解耦方程的优化拟解析装置，包括：

[0037] 弹性波优化归一化拟拉普拉斯算子计算单元，用于对二阶位移常密度弹性波解耦方程进行分解，获得等效的P波、S波解耦的弹性波方程，并利用所述等效的P波、S波解耦的弹性波方程分别获得对应地P波、S波优化归一化拟拉普拉斯算子；

[0038] 弹性波优化拟微分算子的象征确定单元，用于对P波、S波优化归一化拟拉普拉斯算子利用弹性波优化拟微分算子表示，获得P波、S波优化拟微分算子的象征；

[0039] 分解单元，用于对P波、S波优化拟微分算子的象征分别进行分解；

[0040] 弹性波优化拟微分算子的最终表达式获取单元，用于将P波、S波优化拟微分算子象征的分解结果代入到弹性波优化拟微分算子表达式中，获得弹性波优化拟微分算子的最终表达式；

[0041] 波场分量计算单元，用于利用所述弹性波优化微分算子的最终表达式分别对所述等效的P波、S波解耦的弹性波方程进行计算，获得P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。

[0042] 上述技术方案具有如下有益效果：

[0043] 在声波拟解析法的基础上，本技术方案将声波拟拉普拉斯算子改进为弹性波优化归一化拟拉普拉斯算子，该弹性波优化归一化拟拉普拉斯算子，将介质恒定补偿速度，替换为随空间变换的真实介质的纵横波速度。利用弹性波优化拟微分算子表示弹性波优化归一化拟拉普拉斯算子，相对于传统的伪谱法和拟解析法，该弹性波优化拟微分算子不仅包括原始拟微分算子的谱估计，而且还包含一个时间补偿项。因此，弹性波优化拟微分算子可以在空间每个网格点上，精确地补偿二阶时间差分引起的误差。在时间采样间隔较小的情况下，尽管传统的伪谱法仍然具有很高的空间精度，但有明显的时间频散。当采样时间间隔超过了CFL(收敛条件)数时，伪谱法不稳定。因此，在时间采样间隔变大时，虽然伪谱法的空间精度很高，但是容易出现时间频散和不稳定现象。虽然，拟解析法可以有效地补偿时间误差，减少时间频散，但由于补偿速度大于实际介质速度，空间误差会明显增加。在速度变化剧烈的情况下，即使级联多个恒速补偿拟拉普拉斯算子，拟解析法也会产生空间误差。而利用本技术方案提出的弹性波优化拟微分算子，将在时间频散和空间频散方面，优于传统的伪谱法和拟解析法。

附图说明

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

[0045] 图1为本实施例提出一种基于弹性波解耦方程的优化拟解析方法流程图；

[0046] 图2为本实施例提出一种基于弹性波解耦方程的优化拟解析装置框图；

[0047] 图3a为实施例一中伪谱法得到的X分量示意图；

[0048] 图3b为实施例一中伪谱法得到的Z分量示意图；

[0049] 图3c为实施例一中拟解析法得到的X分量示意图；

[0050] 图3d为实施例一中拟解析法得到的Z分量示意图；

[0051] 图3e为实施例一中优化拟解析法得到的X分量示意图；

[0052] 图3f为实施例一中优化拟解析法得到的Z分量示意图；

[0053] 图4为实施例二中倾斜层状模型示意图；

[0054] 图5a为实施例二中优化拟解析法得到的X分量示意图；

[0055] 图5b为实施例二中优化拟解析法得到的Z分量示意图；

[0056] 图5c为实施例二中拟解析法得到的X分量示意图；

[0057] 图5d为实施例二中拟解析法得到的Z分量示意图；

[0058] 图5e为实施例二中伪谱法得到的X分量示意图；

[0059] 图5f为实施例二中伪谱法得到的Z分量示意图；

[0060] 图6a为实施例二中伪谱法得到的水平分量单炮记录图中黑色方框①②部分需局部放大示意图；

[0061] 图6b为实施例二中拟解析法得到的水平分量单炮记录图中黑色方框①②部分需局部放大示意图；

[0062] 图6c为实施例二中优化拟解析法得到的水平分量单炮记录图中黑色方框①②部分需局部放大示意图；

[0063] 图7a为实施例二中伪谱法得到的水平分量单炮记录图中黑色方框①②部分局部放大示意图；

[0064] 图7b为实施例二中拟解析法得到的水平分量单炮记录图中黑色方框①②部分局部放大示意图；

[0065] 图7c为实施例二中优化拟解析法得到的水平分量单炮记录图中黑色方框①②部分局部放大示意图。

具体实施方式

[0066] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0067] 本技术方案的工作原理为：从现有技术我们了解到，利用恒定补偿速度求取拟微分算子的拟解析法方案，最关键的技术是恒定速度的选取以及拟微分算子的计算。因为速度的选取和拟拉普拉斯算子的计算结果，间接和直接决定了正演的时间和空间精度。对于速度变化不同的速度模型，拟解析法可采用不同的方案来求解拟拉普拉斯算子，均可得到较好模拟效果。但是，最后的模拟效果仍然有较大提升空间。现有的拟解析法是基于常密度声波方程，并且针对各向同性介质，并不符合实际地质情况。因此，对于各向异性介质，拟解析法模拟效果较差，此时，我们采用弹性波方程来进行波场模拟。因此，本技术方案将声波拟拉普拉斯算子改进为弹性波优化归一化拟拉普拉斯算子，利用弹性波优化拟微分算子表示弹性波优化归一化拟拉普拉斯算子，该拟微分算子不仅包括原始拟微分算子的谱估计，而且还包含一个时间补偿项，补偿二阶时间差分引起的误差。弹性波优化拟微分算子不能直接应用空间快速傅里叶变换进行计算，若直接计算，其计算量较大。为了减少计算量，本技术方案对弹性波优化拟微分算子中的象征采用低秩分解近似方法进行分解。最终达到以下两方面的效果：当时间采样间隔较小时，本技术方案正演的时间和空间模拟精度均优于传统的伪谱法和拟解析法，但是计算量又没有提高；当时间采样间隔较大时，本技术方案在稳定性上优于伪谱法，并且时间和空间精度都要优于伪谱法和拟解析法。

[0068] 二阶位移常密度弹性波解耦方程的表达式为：

[0069]

[0070] 其中，x为空间矢量，x、y、z表示空间矢量x中的三个垂直方向的坐标；Vp(x)和Vs(x)分别为纵波速度与横波速度，u(x,t)、v(x,t)、w(x,t)分别为质点在混合波场作用下空间矢量x中的三个垂直方向上产生的位移分量。

[0071] 在式(1)中，位移分量时间二阶差分的计算，不仅与纵波速度Vp(x)有关，还与横波速度Vs(x)有关，使得我们不能直接应用下文所提出的弹性波优化拟微分算子，所以将式(1)分解为与其等效的P波、S波解耦的弹性波方程：

[0072]

[0073]

[0074] 其中，式(2)为解耦的弹性波方程中的P波方程，式(3)为解耦的弹性波方程中的S波方程。up(x,t)、vp(x,t)、wp(x,t)分别为质点在P波作用下空间矢量x中的三个垂直方向产生的位移分量；us(x,t)、vs(x,t)、ws(x,t)分别为质点在S波作用下空间矢量x中的三个垂直方向产生的位移分量。那么，质点在混合波场作用下空间矢量x中的三个垂直方向的位移以及P波、S波波场可以表示为：

[0075]

[0076] 在式4中，u(x,t)、v(x,t)、w(x,t)质点在混合波场作用下空间矢量x中的三个垂直方向的位移；p(x,t)、s(x,t)分别表示P波位移分量和S波位移分量。其中，P波位移分量即P波波场，S波位移分量即S波波场。

[0077] 空间位移矢量为：

[0078] E(x,t)＝p(x,t)+s(x,t)  (5)

[0079] 上述解耦的弹性波方程可以在延拓矢量波场的同时，分解并延拓纯纵波和纯横波波场。

[0080] 对于本技术方案来说，P波波场水平分量为：up(x,t)和vp(x,t)，P波波场垂直分量为：wp(x,t)，S波波场水平分量为：us(x,t)和vs(x,t)，S波波场垂直分量为ws(x,t)；在混合波场作用下，u(x,t)、v(x,t)为水平分量，w(x,t)为垂直分量。对公式(2)、公式(3)进行求解，便可求出P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。然后根据公式(4)，便可求出质点在混合波场作用下空间矢量x中的三个垂直方向的位移u(x,t)、v(x,t)、w(x,t)以及P波位移分量p(x,t)、S波位移分量s(x,t)，根据公式(5)便可求出混合波场的的位移分量E(x,t)。

[0081] 各向同性常密度声波方程可写为：

[0082]

[0083] 其中，U(x,t)为地震声压波场，v(x)为地震波波速， 为拉普拉斯算子。

[0084] 在v(x)＝v(常数)的假设条件下，拟拉普拉斯算子(二阶波数-空间算子)F(k)为：

[0085]

[0086] 其中，k为空间波数矢量，Δt是时间采样间隔。

[0087] 由式(6)和式(7)可得波数-空间域拟解析法表达式为：

[0088]

[0089] 在时间-空间域拟解析法表达式为：

[0090]

[0091] 其中，FFT表示空间快速傅里叶正变换，FFT-1表示空间快速傅里叶反变换。

[0092] 对于均匀介质，式(9)中的空间快速傅里叶反变换可直接应用，并且式(9)可提供解析解或近似解析解；对于非均匀介质，速度随空间变化时，不能直接应用式(9)中空间快速傅里叶反变换。

[0093] 为了更为精确地补偿二阶有限差分引起的时间误差，将式(7)中拟拉普拉斯算子的恒定补偿速度v替换为随空间变换的真实速度v(x)。那么，拟拉普拉斯算子可改写为：

[0094]

[0095] 其中，Fopt(k)是空间波数矢量k的函数，且与空间网格坐标有关，即与速度有关。本技术方案将Fopt(k)称为优化拟拉普拉斯算子。

[0096] 将式(10)中的余弦函数使用泰勒展开，那么优化拟拉普拉斯算子可改写为：

[0097]

[0098] 将式(11)替换式(8)中F(k)，相应地，频率-空间域优化拟解析法可改写为：

[0099]

[0100] 在时间-空间域，式(12)可写为任意阶Lax-Wendroff格式：

[0101]

[0102] 其中，

[0103]

[0104] 很明显，Lax-Wendroff格式利用空间导数项来近似补偿时间二阶导数的离散误差。

[0105] 根据Chu和Stoffa提出的归一化拟拉普拉斯算子：

[0106]

[0107] 将式(10)优化拟拉普拉斯算子改进为优化归一化拟拉普拉斯算子：

[0108]

[0109] 由于式(15)只能应用到标量声波方程，而本技术方案应用到地震波方程为弹性波解耦方程，那么弹性波优化归一化拟拉普拉斯算子可表示为：

[0110]

[0111] 其中，kp、ks分别为P波、S波空间波数矢量，ela＝p,s表示P波和S波。不仅是二阶微分算子，而且还能在波数-空间域补偿二阶时间有限差分离散引起的误差。此时，补偿速度分别取纵波速度Vp(x)和横波速度Vs(x)。

[0112] 拟微分算子可以看作是二阶微分算子的推广和一般形式，可以用傅里叶积分形式表示为：

[0113]

[0114] 其中，α(x,kela)称为弹性波优化拟微分算子A[E(x,t)]的象征，E(kela,t)为位移矢量E(x,t)的空间傅里叶变换。

[0115] 对于二阶微分算子和弹性波优化拟微分算子，本技术方案将象征α(x,kela)改写为：

[0116]

[0117] 其中， 在式(16)中已定义。相应地，将公式(18)中的αopt,ela(x,kela)替换公式(17)中的α(x,kela)，弹性波优化拟微分算子可以写为:

[0118]

[0119] 该弹性波优化拟微分算子不仅包括原始拟微分算子的谱估计，而且还包含一个时间补偿项。该时间补偿项可在波数-空间域精确地补偿二阶时间差分引起的误差。补偿速度在式(16)中为介质实际速度，因此，弹性波优化拟微分算子可以在空间每个网格点上，精确地补偿二阶时间差分引起的误差。

[0120] 我们将Aoptela[E(x,t)]简记为：

[0121] Aopt,ela[E(x,t)]＝FT-1{αopt,ela(x,kela)FFT[E(x,t)]}  (20)

[0122] 其中，FT-1表示空间傅里叶反变换。

[0123] 那么，在式(2)(ela＝p时)与式(3)(ela＝s时)中，二阶微分算子与象征αopt,ela(x,kela)之间的关系可表示为：

[0124]

[0125] 由于弹性波优化拟微分算子中包含随空间坐标变化的P波、S波补偿速度vp(x)、vs(x)，所以方程(20)不能直接应用空间快速傅里叶变换进行计算。若直接计算方程(20)，其计算量为 其中Nx为3维网格总数。为了减少计算量，本技术方案对弹性波优化拟微分算子中的象征，采用低秩分解近似方法。

[0126] 可分低秩分解近似主要思想是：选取N个代表性的空间矢量和M个代表性的空间波数矢量，近似原始波数-空间混合域算子。本技术方案将象征近似分解为如下形式：

[0127]

[0128] 其中，W(x,km)是象征αopt,ela(x,kela)的子矩阵,其由与km有关的特定的列元素组成；W(xn,kela)是象征αopt,ela(x,kela)的另一个子矩阵，其由与xn有关的特定行元素组成；αmn代表中间的系数矩阵。

[0129] 通过公式(22)将αopt,ela(x,kela)进行分解，把分解之后的结果代入公式(20)中，那么Aopt,ela[E(x,t)]可改写为：

[0130]

[0131] 式(23)的计算量相当于应用N×2次空间快速傅里叶变换，其计算量为O(NNx log Nx)×2。N一般为一个很小的整数，其大小与式(23)矩阵分解的秩有关。

[0132] 利用公式(23)可以得到等效解耦弹性波方程中的P波方程(对应公式2)和等效解耦弹性波方程中的S波方程(对应公式3)中二阶微分算子的优化拟微分算子，将得到的P波二阶微分算子的优化拟微分算子分别代入公式(2)，将得到的S波二阶微分算子的拟微分算子代入公式(3)，分别对等效解耦弹性波方程中的P波方程和S波方程进行计算，获得P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。

[0133] 本技术方案对弹性波优化拟微分算子中的象征，采用低秩分解近似方法进行分解，大大减少了直接对弹性波优化拟微分算子应用空间快速傅里叶变换的计算量。对于计算效率，在均匀介质情况下，每一个归一化的拟拉普拉斯算子通过低秩分解近似，其秩为2。这意味着相比较于常规的伪谱法，计算归一化的拟拉普拉斯算子需要额外2次快速傅里叶变换。但是，弹性波优化拟解析法能够适用于比较大的时间采样间隔，这将节约不少计算时间。在时间采样间隔较大时，伪谱法不稳定，拟解析法虽然稳定，但是空间频散严重；本技术方案提出的弹性波优化拟解析法无论时间频散还是空间频散，均优于传统的伪谱法和拟解析法；且此情况下，优化拟解析法正演所需要的时间低于伪谱法和拟解析法。

[0134] 基于上述工作原理，本实施例提出一种基于弹性波解耦方程的优化拟解析方法，如图1所示。包括：

[0135] 步骤101)：对二阶位移常密度弹性波解耦方程进行分解，获得等效的P波、S波解耦的弹性波方程，并利用所述等效的P波、S波解耦的弹性波方程分别获得对应地P波、S波优化归一化拟拉普拉斯算子；

[0136] 步骤102)：对P波、S波优化归一化拟拉普拉斯算子利用弹性波优化拟微分算子表示，获得P波、S波优化拟微分算子的象征；

[0137] 步骤103)：对P波、S波优化拟微分算子象征分别进行分解；

[0138] 步骤104)：将对P波、S波优化拟微分算子象征的分解结果代入到弹性波优化拟微分算子的表达式中，获得弹性波优化拟微分算子的最终表达式；

[0139] 步骤105)：利用所述弹性波优化微分算子的最终表达式分别对所述等效的P波、S波解耦的弹性波方程进行计算，获得P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。

[0140] 本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一般计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。

[0141] 如图2所示，为本实施例提出一种基于弹性波解耦方程的优化拟解析装置框图。包括：

[0142] 弹性波优化归一化拟拉普拉斯算子计算单元201，用于对二阶位移常密度弹性波解耦方程进行分解，获得等效的P波、S波解耦的弹性波方程，并利用所述等效的P波、S波解耦的弹性波方程分别获得对应地P波、S波优化归一化拟拉普拉斯算子；

[0143] 弹性波优化拟微分算子的象征确定单元202，用于对P波、S波优化归一化拟拉普拉斯算子利用弹性波优化拟微分算子表示，获得P波、S波优化拟微分算子的象征；

[0144] 分解单元203，用于对P波、S波优化拟微分算子的象征分别进行分解；

[0145] 弹性波优化拟微分算子的最终表达式获取单元204，用于将P波、S波优化拟微分算子象征的分解结果代入到弹性波优化拟微分算子的表达式中，获得弹性波优化拟微分算子的最终表达式；

[0146] 波场分量计算单元205，用于利用所述弹性波优化微分算子的最终表达式分别对所述等效的P波、S波解耦的弹性波方程进行计算，获得P波波场的水平分量和垂直分量以及S波波场的水平分量和垂直分量。

[0147] 本领域技术人员还可以了解到本发明实施例列出的各种功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。

[0148] 此外，尽管在上文详细描述中提及了装置的若干单元，但是这种划分仅仅并非强制性的。实际上，根据本发明的实施方式，上文描述的两个或更多单元的特征和功能可以在一个单元中具体化。同样，上文描述的一个单元的特征和功能也可以进一步划分为由多个单元来具体化。

[0149] 实施例

[0150] 为了能够更加直观的描述本发明的特点和工作原理，下文将结合实际运用场景来描述。

[0151] 实施例一

[0152] 本实施例采用2D均匀模型，其在X方向和Z方向采样间隔均为10m。横波和纵波速度分别为1500m/s、2500m/s。Ricker子波的主频为20Hz,最大频率为50Hz。震源位于模型中心，并同时产生纵波和横波。每个波场占最小网格点数为：2个，这已达到Nyquist采样定理极限。本次模型试算不采用吸收边界条件。分别使用伪谱法、拟解析法以及优化拟解析法进行正演模拟，时间采样间隔Δt＝2.5ms。图3a和图3b分别为伪谱法得到的X分量和Z分量，图3c和图3d分别为拟解析法得到的X分量和Z分量，图3e和图3f为优化拟解析法得到的X分量和Z分量。分别对比三种测试方法得到的X分量以及Z分量，观测结果表明：当空间采样间隔为10m，时间采样间隔为2.5ms时，伪谱法已经不稳定；虽然拟解析法、优化拟解析法仍然稳定，但优化拟解析法的空间频散(黑色箭头所示)比拟解析法更少。

[0153] 实施例二

[0154] 本实施例采用二维层状模型，应用本技术方案提出的优化拟解析法，来验证其有效性。并与伪谱法以及拟解析法进行对比。二维层状模型包括上部一个水平反射层和下部一个倾斜反射层，倾斜反射层为阶梯界面。从上到下，每层的纵波速度为：2000m/s、2350m/s、2500m/s，横波速度为：1500m/s、1800m/s、2300m/s。模型区域大小为250×250，水平方向和垂直方向网格间距均为10m，如图4所示。正演模拟采用Richer子波，主频为30HZ，可同时产生纵波和横波。波场模拟的时间采样间隔为1.5ms，地震记录总的时间采样点数为1501个。炮点位于模型网格点(125,5)处，检波器均匀分布于模型整个表面。图5a～图5b为优化拟解析法得到的X分量和Z分量，图5c～图5d为拟解析法得到的X分量和Z分量，图5e～图5f为伪谱法得到的X分量和Z分量。图5a～图5f均为t＝1.2s时刻的X分量和Z分量波场快照图。分别对比三种测试方法得到的X分量以及Z分量，观测结果表明：在时间采样间隔1.5ms时，伪谱法具有明显的时间频散(黑色箭头所示)；拟解析法虽然压制了时间频散(黑色箭头所示)，但是由于补偿速度大于实际速度，会产生空间频散(浅色箭头所示)；优化拟解析法，在压制时间和空间频散上，优于伪谱法和拟解析法。

[0155] 将图6a所示的水平分量单炮记录图中黑色方框①②部分需局部放大，对应图7a中的①和②，将图6b所示的水平分量单炮记录图中黑色方框①②部分需局部放大，对应图7b中的①和②，将图6c所示的水平分量单炮记录图中黑色方框①②部分需局部放大，对应图7c中的①和②。可以观察到方块①中：图6a伪谱法，直达波有明显的时间频散(黑色箭头所示)，图6b拟解析法，很好的压制了时间频散，但是有较严重的空间频散(浅色箭头所示)，图
6c优化拟解析法，在压制时间频散的同时，也消除了拟解析法由于补偿速度选取过大产生的空间频散；方块②中同样如此。

[0156] 倾斜层状模型的数值模拟也证明：优化拟解析法在压制时间频散和空间频散上，相较于传统的伪谱法和拟解析法，更有优势。

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