技术领域
[0001] 本
发明涉及
地震成像技术领域，特别涉及一种三维声波方程任意域多尺度全波形反演方法及装置。
背景技术
[0002] 三维声波方程多尺度全波形反演是一种有效的地震成像方法，它可以在时间域、
频率域、Laplace-Fourier域进行。此方法从简单的初始模型开始，直接利用所接收的地震数据直接得到关于地下介质参数的定量信息。这种方法的
基础是三维声波方程正演。因为反演中的梯度计算以及步长的选取都通过正演来实现。正演的
精度和效率直接影响着三维声波方程多尺度全波形反演的精度和效率。如何发展保持效率的三维声波方程高精度正演
算法，是三维声波方程多尺度全波形地震成像方法需要解决的重要问题。
[0003] 现有的三维声波方程多尺度全波形反演中所使用的正演方法在时间离散上采用二阶近似,而在空间近似上采用高阶近似。这种正演方法，在时间步长比较小的时候，可以满足精度要求。但当时间步长较大时，就存在较大的误差，这是因为在时间近似上只是采用了二阶近似。因此，为了得到高精度的正演方法，需要发展高阶的时间近似方法，现有的高阶的时间近似方法包括快速展开方法和Lax-Wendroff方法，虽然精度比较高，但效率较低，特别是计算步数比较大的时候，积累误差比较大。
[0004] 现有的三维声波方程多尺度全波形反演中所使用的最优化方法为有限记忆拟
牛顿法。这种方法对于较小规模的速度模型，非常有效。但对于大规模的模型，这种方法的存储量和计算量都非常大，因此反演效率比较低。
发明内容
[0005] 本发明
实施例的主要目的在于提出一种三维声波方程任意域多尺度全波形反演方法及装置，解决三维声波方程多尺度全波形反演中所使用的正演算法精度比较低以及所使用的最优化方法效率比较低的问题。
[0006] 为实现上述目的，本发明提供了一种三维声波方程任意域多尺度全波形反演方法，包括：
[0007] 给定目标区域的实际采集波场，确定离散模型的大小；
[0008] 根据离散模型的大小确定计算精度和效率，利用所述精度和效率确定正演方法，并采用正演方法和不同域之间的波场变换来确定搜索方向；其中，所述正演方法为4阶辛算法；
[0009] 确定全波形反演的尺度；
[0010] 确定目标区域的初始速度模型；
[0011] 对每一个全波形反演的尺度，利用所述搜索方向对所述初始速度模型进行速度
迭代更新，直至满足迭代终止条件；迭代终止之后输出的速度模型为全波形反演结果。
[0012] 可选的，在本发明一实施例中，所述迭代终止条件包括下降率达不到一定要求或迭代次数超过某一数值。
[0013] 可选的，在本发明一实施例中，所述搜索方向确定步骤包括：
[0014] 利用所述正演方法计算目标函数的梯度；
[0015] 对所述目标函数的梯度进行虚源预条件处理，得到搜索方向。
[0016] 可选的，在本发明一实施例中，所述目标函数的梯度表达式为：
[0017]
[0018] 其中，G为目标函数的梯度，Pd和Pu分别为正传波场和反传残差波场，为实际采集波场变换算子，为实际采集波场的相关运算并乘以相应的因子。
[0019] 可选的，在本发明一实施例中，所述搜索方向的表达式为：
[0020] D＝-G./(S+εmax(S)),
[0021] 其中，D为搜索方向，“./”表示向量间的逐分量相除运算，S为由虚源得到的预条件向量，ε为预条件参数，max(S)为预条件向量的最大值。
[0022] 可选的，在本发明一实施例中，所述速度迭代更新的表达式为：
[0023] Vnew＝Vold+αD,
[0024] 其中，Vnew表示更新后的速度，Vold表示更新前的速度，α为步长因子，通过抛物线拟合方法得到。
[0025] 对应地，为实现上述目的，本发明提供了一种三维声波方程任意域多尺度全波形反演装置，包括：
[0026] 离散模型确定单元，用于给定目标区域的实际采集波场，确定离散模型的大小；
[0027] 搜索方向确定单元，用于根据离散模型的大小确定计算精度和效率，利用所述精度和效率确定正演方法，并采用正演方法和不同域之间的波场变换来确定搜索方向；其中，所述正演方法为4阶 辛算法；
[0028] 反演尺度确定单元，用于确定全波形反演的尺度；
[0029] 初始速度模型确定单元，用于确定目标区域的初始速度模型；
[0030] 迭代计算单元，用于对每一个全波形反演的尺度，利用所述搜索方向对所述初始速度模型进行速度迭代更新，直至满足迭代终止条件；迭代终止之后输出的速度模型为全波形反演结果。
[0031] 可选的，在本发明一实施例中，所述搜索方向确定单元包括：
[0032] 目标函数的梯度计算模
块，用于利用所述正演方法计算目标函数的梯度；
[0033] 虚源预条件处理模块，用于对所述目标函数的梯度进行虚源预条件处理，得到搜索方向。
[0034] 可选的，在本发明一实施例中，所述目标函数的梯度计算模块涉及的目标函数的梯度表达式为：
[0035]
[0036] 其中，G为目标函数的梯度，Pd和Pu分别为正传波场和反传残差波场，为实际采集波场变换算子，为实际采集波场的相关运算并乘以相应的因子。
[0037] 可选的，在本发明一实施例中，所述虚源预条件处理模块涉及的计算表达式为：
[0038] D＝-G./(S+εmax(S)),
[0039] 其中，D为搜索方向，“./”表示向量间的逐分量相除运算，S为由虚源得到的预条件向量，ε为预条件参数，max(S)为预条件向量的最大值。
[0040] 上述技术方案具有如下有益效果：
[0041] 本技术方案提供的多尺度全波形反演方法，在正演方面，基于4阶 辛算法的正演方法,在提高算法的精度的同时,保持了算法的效率,并且随着计算步数的增加,优势更加明显；在最
优化算法方面所提出的适应于任意域的虚源预条件梯度方法,不仅效率高,而且应用范围广。
附图说明
[0042] 为了更清楚地说明本发明实施例或
现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0043] 图1为本发明提出的一种三维声波方程任意域多尺度全波形反演方法
流程图；
[0044] 图2为本发明提出的一种三维声波方程任意域多尺度全波形反演装置
框图；
[0045] 图3为三维声波方程任意域多尺度全波形反演装置中搜索方向确定单元的功能框图；
[0046] 图4为本实施例中初始速度模型示意图；
[0047] 图5为本实施例中利用本技术方案反演得到的最终速度模型示意图。
具体实施方式
[0048] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0049] 本技术方案的工作原理为：本
专利为三维声波方程多尺度全波形反演提供了一种基于4阶 辛算法的三维声波方程正演方法。这种方法在空间离散上利用10阶有限差分方法，在时间离散上利用4阶 辛算法。在提高算法的精度的同时，保持了计算的效率。因为这种方法是一种保结构方法，其积累误差比较少，整体精度比较高。
[0050] 进一步地，三维声波方程多尺度全波形反演是通过对由实际采集波场和计算波场构造的目标函数求极值来实现对地下介质参数成像的。对目标函数求极值所使用的最优化方法包括梯度法、预条件梯度法、共轭梯度法、Gauss-Newton法，拟牛顿法和牛顿法。拟牛顿方法类中的有限记忆拟牛顿法比较常用，但其计算量和存储量都比较大，特别是对于大规模实际模型，这种方法计算量和存储量巨大，效率比较低。相比而言，预条件梯度法计算量和存储量都比较小，效率比较高，特别对于基于虚源预条件的梯度法，其效率非常高。对于大规模实际模型的全波形反演来讲，虚源预条件梯度法是一种非常有效的最优化方法。但对于不同的计算域，还没有对虚源预条件梯度法作统一的处理。
[0051] 三维声波方程为：
[0052]
[0053] 这里p(x,y,z,t)为压
力波场，v(x,y,z)为速度，而f(x,y,z,t)为源。
[0054] 设 为方程(1)的数值解，即通过数值计算得到的波场。这里，i＝1,...,I,j＝1,...,J,k＝1,...,K,n＝1,...,N,而Δx,Δy,Δz和Δt分别为空间间隔和时间步长。数值解依赖于离散速度 另一
方面，设 为实际采集的波场。我们构造下面的目标函数：
[0055]
[0056] 这里， 为
检波器所在的指标集， 为波场变换算子：对于时间域， 为恒等变换；对于频率域，为Fourier变换；对于Laplace-Fourier域，为Laplace-Fourier变换。
[0057] 全波形反演就是求解最优的速度V＝[v1,1,1,...,vI,J,K]T，使得目标函数(2)的值最小。如图1所示，基于4阶 辛算法的三维声波方程多尺度全波形反演包括下面几个步骤：
[0058] 步骤101)：给定目标区域的实际采集波场，确定离散模型的大小；
[0059] 步骤102)：根据离散模型的大小确定计算精度和效率，利用所述精度和效率确定正演方法，并采用正演方法和不同域之间的波场变换来确定搜索方向；其中，所述正演方法为4阶 辛算法；
[0060] 因为正演方法是全波形反演的基础，正演方法的精度和效率在很大程度上决定了全波形反演的精度和效率，因此这一步十分关键。本技术方案提出基于4阶 辛算法的三维声波方程的正演方法。其格式为：
[0061]
[0062]
[0063] 其中，
[0064]
[0065]
[0066]
[0067]
[0068] 此格式的精度为ο(Δt4+Δx10+Δy10+Δz10)。由于它是一个辛格式，所以误差积累比较小，整体精度比较高。
[0069] 利用4阶 辛算法计算目标函数的梯度：
[0070]
[0071] 这里G为目标函数的梯度，Pd和Pu分别为正传波场和反传残差波场， 为前面提到的波场变换算子，而 为波场的相关运算并乘以相应的因子。
[0072] 然后对梯度进行虚源预条件处理，得到搜索方向：
[0073] D＝-G./(S+εmax(S)), (6)
[0074] 这里D为搜索方向，“./”表示向量间的逐分量相除运算，S为由虚源得到的预条件向量，ε为预条件参数，而max(S)为预条件向量的最大值。
[0075] 步骤103)：确定全波形反演的尺度；
[0076] 在本实施例中，对于全波形反演的尺度来说，如果是时间域反演，则确定几个频率段；如果是频率域反演，则确定若干个频率；如果是Laplace-Fourier域反演，则确定若干个Laplace-Fourier频率。然后，从大尺度到小尺度逐个进行反演。
[0077] 步骤104)：确定目标区域的初始速度模型；
[0078] 在本实施例中，初始速度模型一般为简单的模型：均匀模型、线性增加模型或平滑模型。
[0079] 步骤105)：对每一个全波形反演的尺度，利用所述搜索方向对所述初始速度模型进行速度迭代更新，直至满足迭代终止条件；迭代终止之后输出的速度模型为全波形反演结果。
[0080] 首先，判断是否满足迭代终止条件。终止条件包括下降率达不到一定要求或迭代次数超过某一数值。如果满足迭代终止条件，则终止迭代，输出更新后的速度模型。如果不满足迭代终止条件，则将更新后的速度模型作为新的初始速度模型继续迭代计算。
[0081] 如图2所示，为本发明提出的一种三维声波方程任意域多尺度全波形反演装置框图。包括：
[0082] 离散模型确定单元201，用于给定目标区域的实际采集波场，确定离散模型的大小；
[0083] 搜索方向确定单元202，用于根据离散模型的大小确定计算精度和效率，利用所述精度和效率确定正演方法，并采用正演方法和不同域之间的波场变换来确定搜索方向；其中，所述正演方法为4阶 辛算法；
[0084] 反演尺度确定单元203，用于确定全波形反演的尺度；
[0085] 初始速度模型确定单元204，用于确定目标区域的初始速度模型；
[0086] 迭代计算单元205，用于对每一个全波形反演的尺度，利用所述搜索方向对所述初始速度模型进行速度迭代更新，直至满足迭代终止条件；迭代终止之后输出的速度模型为全波形反演结果。
[0087] 如图3所示，为三维声波方程任意域多尺度全波形反演装置中搜索方向确定单元的功能框图。所述搜索方向确定单元202包括：
[0088] 目标函数的梯度计算模块2021，用于利用所述正演方法计算目标函数的梯度；所述目标函数的梯度计算模块2021涉及的目标函数的梯度表达式为：
[0089]
[0090] 其中，G为目标函数的梯度，Pd和Pu分别为正传波场和反传残差波场，为实际采集波场变换算子，为实际采集波场的相关运算并乘以相应的因子。
[0091] 虚源预条件处理模块2022，用于对所述目标函数的梯度进行虚源预条件处理，得到搜索方向。所述虚源预条件处理模块2022涉及的计算表达式为：
[0092] D＝-G./(S+εmax(S)),
[0093] 其中，D为搜索方向，“./”表示向量间的逐分量相除运算，S为由虚源得到的预条件向量，ε为预条件参数，max(S)为预条件向量的最大值。
[0094] 本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，可以通过
计算机程序来指令相关的
硬件来完成，所述的程序可存储于一般计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。
[0095] 本领域技术人员还可以了解到本发明实施例列出的各种功能是通过硬件还是
软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。
[0096] 下面以目标区域的实际采集波场的离散模型为三维Marmousi模型作为实施例说明本发明的具体实施方式，离散模型的大小为L×M×N＝166×39×98(包括吸收边界)，网格间距Δx＝Δy＝Δz＝10m。采集波场为合成波场。我们使用时间域多尺度全波形反演方法。分三个尺度：0～10Hz，0～15Hz，0～50Hz。对于初始速度模型，我们使用一个平滑模型，如图4所示。对于每一个尺度，我们迭代5次。图5为反演得到的最终速度模型。和真实模型相比，取得了很好的成像效果。
[0097] 由本实施例可知，本技术方案为三维声波方程多尺度全波形反演建立高效率且保证了高精确度。
[0098] 以上具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
