CN114966837A - 基于卷积型w2距离目标函数的全波形反演方法及装置 - Google Patents

Abstract

本说明书提供了基于卷积型W₂距离目标函数的全波形反演方法及装置。该方法包括：获取初始速度模型，根据初始速度模型计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据；获取观测地震数据；根据初始的模拟地震数据、观测地震数据和卷积型W₂距离目标函数，计算伴随源；根据伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；根据伴随波场和逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；根据卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。基于上述方法能够解决现有方法中存在的周波跳跃的技术问题，实现速度模型的精确构建。

Description

基于卷积型W2距离目标函数的全波形反演方法及装置

技术领域

本说明书涉及地震成像技术领域，特别地，涉及基于卷积型W₂距离目标函数的全波形反演方法及装置。

背景技术

地震成像质量不仅受偏移算法影响，还依赖于速度建模精度。通常速度模型越精确，速度估计越准确，得到偏移剖面的空间结构相对越清晰；同时，偏移剖面的振幅和相位信息也越可靠。

现有的建模方法中，例如，全波形反演(Full Waveform Inversion,FWI)能综合利用地震波走时、振幅、相位等全波场信息，获得高分辨率反演结果，具有广泛应用前景。FWI最优化问题具体又可以通过全局和局部两类算法实现。其中，全局算法可以避免陷入局部极值，然而，其计算成本过高，难以实际推广。在局部算法中，FWI大多采用模拟和观测数据残差的L₁、L₂或两种范数的混合作为反演的目标函数，但在低波数不准确的情况下，基于该方法使用梯度算法的FWI很容易陷入局部极小值，产生周波跳跃问题。

因此，目前亟需一种解决周波跳跃的全波形反演方法。

发明内容

本说明书提供了基于卷积型W₂距离目标函数的全波形反演方法及装置，能够解决现有方法中存在的周波跳跃的技术问题，得到精准的速度模型。

本说明书实施例的目的是提供基于卷积型W₂距离目标函数的全波形反演方法，包括：

获取初始速度模型，根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据；

获取观测地震数据；根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源；

根据所述伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；

根据所述伴随波场和所述逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；

根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。

进一步地，所述方法的另一个实施例中，所述根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据，包括：

在计算区域内，根据初始速度模型，构建震源波动方程；

在计算区域外，加载分裂式完美匹配层(Perfectly Matched Layer，PML)，构建分裂式完美匹配层波动方程；

根据所述震源波动方程、所述分裂式完美匹配层波动方程，计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据。

进一步地，所述方法的另一个实施例中，所述根据所述震源波动方程、所述分裂式完美匹配层波动方程，计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据，包括：

按照以下方式计算当前轮的正传震源波场：

获取与分裂式完美匹配层边界相关的上一轮的关联卷积项；

在计算区域外，根据与分裂式完美匹配层边界相关的上一轮的关联卷积项和分裂式完美匹配层波动方程，得到PML区域震源波场；

在计算区域内，根据震源波动方程和PML区域震源波场，得到当前轮的正传震源波场；

更新与分裂式完美匹配层边界相关的关联卷积项，得到当前轮的关联卷积项；

检测正时间序列上的正传震源波场是否计算完毕，如果没有计算完毕，计算下一轮的正传震源波场；如果计算完毕，根据正传震源波场确定出初始的模拟地震数据。

进一步地，所述方法的另一个实施例中，所述根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源，包括：

根据初始的模拟地震数据，计算归一化后的模拟地震数据；

根据归一化后的模拟地震数据，计算反函数；

根据反函数、归一化后的模拟地震数据、观测地震数据、卷积型W₂距离目标函数，计算卷积型W₂距离目标函数值；

根据卷积型W₂距离目标函数值，计算伴随源。

进一步地，所述方法的另一个实施例中，所述根据反函数、归一化后的模拟地震数据、观测地震数据、卷积型W₂距离目标函数，计算卷积型W₂距离目标函数值，包括：

按照以下算式，计算卷积型W₂距离目标函数值：

其中，W₂为卷积型W₂距离目标函数值，f＝{f_i}(i＝1,2,...,N)为初始的模拟地震数据，d＝{d_i}(i＝1,2,...,N)为观测地震数据，N为最大采样时刻，M为模型的空间网格数，Δt为采样间隔，t_i＝(i-1)Δt(i＝1,2,...,N)为第i个采样点，G_i ^-1表示求反函数，G_i ^-1(F_i)表示函数G_i取值为F_i时所对应的时刻，

为归一化后的观测地震数据，

为归一化后的模拟地震数据。

进一步地，所述方法的另一个实施例中，所述基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场，包括：

按照以下方式进行当前轮的震源波场逆时重构：

获取与分裂式完美匹配层边界相关的上一轮的关联卷积项；

根据局部检查点值，从最大时刻逆时计算PML区域震源波场，得到PML区域震源波场逆时间序列上的结果；

根据PML区域震源波场逆时间序列上的结果，计算得到重构震源波场逆时间序列上的结果；

根据局部检查点值，更新与分裂式完美匹配层边界相关的关联卷积项，得到当前轮的关联卷积项；

检测当前时间段内逆时间序列上的重构震源波场是否计算完毕，以确定是否触发进行下一轮的震源波场逆时重构。

进一步地，所述方法的另一个实施例中，所述根据局部检查点值，从最大时刻逆时计算PML区域震源波场，得到PML区域震源波场逆时间序列上的结果，包括：

按照以下方式，计算得到逆时间序列上的结果中的当前逆时间序列PML区域震源波场数据值：

在计算出当前逆时间序列PML区域震源波场数据值后，将当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点进行比较；

在确定当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点相同的情况下，使用局部检查点中该时间点对应的PML区域震源波场数据值替换当前逆时间序列PML区域震源波场数据值；

在确定当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点不相同的情况下，利用该时刻之前最近的局部检查点处的PML区域震源波场逆时间序列计算该时刻的PML区域震源波场数据值。

进一步地，所述方法的另一个实施例中，所述根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，包括：

根据卷积型W₂距离目标函数的梯度，确定速度模型的更新方向；并计算出相应的更新步长；

根据所述更新方向和更新步长，更新当前速度模型。

另一方面，本申请提供了基于卷积型W₂距离目标函数全波形反演装置，包括：

第一处理模块，用于获取初始速度模型，根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据；

计算模块，用于获取观测地震数据；根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源；

第二处理模块，用于根据所述伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；

求取模块，用于根据所述伴随波场和所述逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；

更新模块，用于根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。

再一方面，本申请还提供了一种计算机可读存储介质，其上存储有计算机指令，所述计算机可读存储介质执行所述指令时实现上述基于卷积型W₂距离目标函数的全波形反演方法。

本说明书提供的基于卷积型W₂距离目标函数的全波形反演方法及装置，通过获取初始速度模型，根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据；获取观测地震数据；根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源；根据所述伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；根据所述伴随波场和所述逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。依据本说明提供的方法，能够消除子波估计不准确对全波形反演的影响，解决了现有方法中存在的周波跳跃、子波误差的技术问题，实现速度模型的精确构建，还可以拓展全波形反演的适用范围。

并且，在具体更新当前速度模型时，通过计算并利用卷积型W₂距离目标函数的梯度，以根据初始的模拟地震数据、观测地震数据的匹配程度来决定下一次模型参数的更新方向；再通过更新模型参数，对当前速度模型进行更新，最终使得初始的模拟地震数据与观测地震数据达到最佳匹配。

进一步，还通过基于双层检查点进行震源波场逆时重构，使得正传震源波场的递推方式与重构震源波场时序一致，便于进行基于正传波场与重构震源波场的零延时互相关的梯度计算，与传统的边界存储方法相比，该方法只需要四次求解波动方程，能够降低数据处理量。

此外，还通过根据相应的运算规则，够将本说明书提供的基于卷积型W₂距离目标函数的全波形反演方法部署于GPU(Graphics Processing Unit)并行计算平台上实现，从而还能够达到保证卷积型W₂距离目标函数和伴随源的计算不会显著增加反演计算量的同时，实现了高效稳健的全波形反演的技术效果。

附图说明

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

图1是本说明书提供的基于卷积型W₂距离目标函数的全波形反演方法一个实施例的流程示意图；

图2a是本说明书一个实施例中的时移子波与准确子波凸性对比示意图；

图2b是本说明书一个实施例中的一维时移数据不同目标函数的凸性对比示意图；

图3a是本说明书一个实施例中的加常数非负化预处理方法对一维信号时移量的凸性分析结果图；

图3b是本说明书一个实施例中的幂指数非负化预处理方法对一维信号时移量的凸性分析结果图；

图3c是本说明书一个实施例中的平方非负化预处理方法对一维信号时移量的凸性分析结果图；

图4是本说明书一个实施例中的双层检查点震源波场逆时重构示意图；

图5a是本说明书一个实施例中的Marmousi速度准确模型示意图；

图5b是本说明书一个实施例中的Marmousi速度初始模型示意图；

图6a是本说明书一个实施例中的相同主频不准确子波类型Marmousi模型卷积型L₂距离目标函数全波形反演第8次迭代梯度示意图；

图6b是本说明书一个实施例中的相同主频不准确子波类型Marmousi模型卷积型W₂距离目标函数全波形反演第8次迭代梯度示意图；

图7a是本说明书一个实施例中的相同主频不准确子波类型Marmousi模型卷积型L₂距离目标函数第一个频带全波形反演结果示意图；

图7b是本说明书一个实施例中的相同主频不准确子波类型Marmousi模型卷积型W₂距离目标函数第一个频带全波形反演结果示意图；

图7c是本说明书一个实施例中的相同主频不准确子波类型Marmousi模型卷积型W₂距离目标函数最终频带全波形反演结果示意图；

图8a是本说明书一个实施例中的Overthrust速度准确模型示意图；

图8b是本说明书一个实施例中的Overthrust速度初始模型示意图；

图9a是本说明书一个实施例中的无噪观测地震数据示意图；

图9b是本说明书一个实施例中的信噪比为1dB的观测地震数据示意图；

图10a是本说明书一个实施例中的无噪观测地震数据迭代40次的反演结果示意图；

图10b是本说明书一个实施例中的信噪比为1dB的观测地震数据迭代40次的反演结果示意图；

图10c是本说明书一个实施例中的信噪比为2dB的观测地震数据迭代40次的反演结果示意图；

图10d是本说明书一个实施例中的信噪比为4dB的观测地震数据迭代40次的反演结果示意图；

图10e是本说明书一个实施例中的信噪比为20dB的观测地震数据迭代40次的反演结果示意图；

图11a是本说明书一个实施例中的不同信噪比观测地震数据在x＝2000m处反演结果垂直速度对比示意图；

图11b是本说明书一个实施例中的不同信噪比观测地震数据在x＝3500m处反演结果垂直速度对比示意图；

图12a是本说明书一个实施例中的无噪观测地震数据反演40次归一化目标函数值下降值示意图；

图12b是本说明书一个实施例中的信噪比为1dB的观测地震数据反演40次归一化目标函数值下降值示意图；

图12c是本说明书一个实施例中的信噪比为2dB的观测地震数据反演40次归一化目标函数值下降值示意图；

图12d是本说明书一个实施例中的信噪比为4dB的观测地震数据反演40次归一化目标函数值下降值示意图；

图12e是本说明书一个实施例中的信噪比为20dB的观测地震数据反演40次归一化目标函数值下降值示意图；

图13是本说明书一个实施例中的迭代40次后不同信噪比观测地震数据的反演结果与无噪数据反演结果后误差对比示意图；

图14是本说明书提供的基于卷积型W₂距离目标函数的全波形反演装置一个实施例的模块结构示意图。

具体实施方式

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

考虑到现有的FWI(Full Waveform Inversion,全波形反演)方法大多采用局部算法，即利用目标函数对模型参数的梯度等信息来获得模型参数的更新量，不断迭代以获得最终的反演结果。在获取梯度信息前需要进行波形匹配，在波形匹配过程中，如果初始的模拟地震数据与观测地震数据的波形相位差大于半个周期，将导致不同周期波形的错误匹配，并最终导致反演结果陷入局部极小值，这种现象称为周波跳跃。与传统的L₁、L₂目标函数相比，基于最优传输距离目标函数考虑的是地震信号之间的整体差异，而不是简单地点对点比较。在基于最优传输距离目标函数中，地震波的传播时间和相位等低频信息被更多地使用，所以具有更好的凸性。因此，基于最优传输距离目标函数对周波跳跃现象具有突出的抑制能力。

进一步，还考虑到进行全波形反演时通常假设震源子波已知，只对模型参数进行反演，然而在实际情况中震源子波是难以准确估计的。当子波不准确时，震源引起的误差会错误地归咎于模型参数，导致模型收敛到错误的方向，反演结果精度降低。采用卷积算法对子波进行处理，可以消除误差，有效解决震源子波不准确的问题。

针对现有方法存在的上述问题以及产生上述问题的具体原因，本申请考虑可以引入基于最优传输距离目标函数进行全波形反演以提升周波跳跃抑制能力；同时，可以采用卷积算法消除震源子波误差，以提升全波形反演精度。

基于上述思路，本说明书提出基于卷积型W₂距离目标函数的全波形反演方法。首先，获取初始速度模型，根据初始速度模型计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据；获取观测地震数据；根据初始的模拟地震数据、观测地震数据和卷积型W₂距离目标函数，计算伴随源；根据伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；根据伴随波场和逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；最后，根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。参阅图1所示，本说明书实施例提供了基于卷积型W₂距离目标函数的全波形反演方法。具体实施时，该方法可以包括以下内容。

S101：获取初始速度模型，根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据。

在一些实施例中，上述获取初始速度模型，可以通过以下方式：在勘探区域，布设二维或三维测线，使用炸药震源激发地震波，沿着测线等间距布设多个检波器接收地震波信号；检波器接收到地震波信号后，以等时间间隔离散采样初始速度模型，并以数字形式记录于存储介质中。所述初始速度模型具体可以包括激发参数、排列参数、接收参数。

在一些实施例中，上述根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据，具体实施时，可以包括：

S1：在计算区域内，根据初始速度模型，构建震源波动方程；

S2：在计算区域外，加载分裂式完美匹配层，构建分裂式完美匹配层波动方程；

S3：根据所述震源波动方程、所述分裂式完美匹配层波动方程，计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据。

在一些实施例中，上述完美匹配层(Perfectly Matched Layer，PML)是在计算区域边界设置的一种特殊介质层，地震波可以无反射地穿过分界面而进入完美匹配层。

在一些实施例中，上述在计算区域内，根据初始速度模型，构建震源波动方程，包括：根据初始速度模型提取得到采样间隔、震源坐标、震源子波；再根据采样间隔、震源坐标、震源子波，在计算区域内，构建震源波动方程。

在一些实施例中，上述根据采样间隔、震源坐标、震源子波，在计算区域内，构建震源波动方程，包括：

按照以下算式构建震源波动方程：

式中，

为震源波场，v₀为初始模型速度，s(t)为震源子波，δ为Delta函数，(x_s,z_s)为震源坐标，x为空间任意一点横坐标，z为空间任意一点纵坐标,t为时间。

在一些实施例中，上述在计算区域外，加载分裂式完美匹配层，构建分裂式完美匹配层波动方程，具体实施时，可以包括：

按照以下算式构建分裂式完美匹配层波动方程：

式中，

为PML区域震源波场，d_x和d_z分别为沿x方向和z方向的衰减函数，

和

为PML区域的震源波场分量，C_x和C_z为与分裂式完美匹配层区域有关的关联卷积项，可以表示为：

在一些实施例中，上述根据所述震源波动方程、所述分裂式完美匹配层波动方程，计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据，具体实施时，可以包括：

按照以下方式计算当前轮的正传震源波场：

S1：获取与分裂式完美匹配层边界相关的上一轮的关联卷积项；

S2：在计算区域外，根据与分裂式完美匹配层边界相关的上一轮的关联卷积项和分裂式完美匹配层波动方程，得到PML区域震源波场；

S3：在计算区域内，根据震源波动方程和PML区域震源波场，得到当前轮的正传震源波场；

S4：更新与分裂式完美匹配层边界相关的关联卷积项，得到当前轮的关联卷积项；

S5：检测正时间序列上的正传震源波场是否计算完毕，如果没有计算完毕，计算下一轮的正传震源波场；如果计算完毕，根据正传震源波场确定出初始的模拟地震数据。

在一些实施例中，上述在计算区域外，根据与分裂式完美匹配层边界相关的上一轮的关联卷积项和分裂式完美匹配层波动方程，得到PML区域震源波场，具体实施时，可以包括：

按照以下算式得到PML区域震源波场：

式中，

为n+1时刻的PML区域震源波场，即经过离散化处理的PML区域震源波场，Δt表示采样间隔，i、j为该计算点在二维格网中的位置横、纵坐标，D_xx为x方向的二阶导数算子，D_zz为z方向的二阶导数算子。

按照以下算式得到当前轮的与分裂式完美匹配层边界相关的关联卷积项：

表示n时刻与分裂式完美匹配层边界相关的关联卷积项。

在一些实施例中，上述在计算区域内，根据震源波动方程和PML区域震源波场，得到当前轮的正传震源波场，具体实施时，包括：

按照以下算式采用傅里叶伪谱法计算空间偏导数：

式中，

为空间波数，k_x、k_z分别为沿x和z方向的空间波数，FFT与IFFT分别为快速傅里叶变换和其反变换，FTT与IFFT同时也是关于PML区域震源波场的函数。

按照以下算式得到正传震源波场：

式中，

分别表示n-1时刻、n时刻、n+1时刻的正传震源波场，即表示经过离散化处理的震源波场。

在一些实施例中，上述根据正传震源波场确定出初始的模拟地震数据，包括：

按照以下算式得到初始的模拟地震数据：

f为初始的模拟地震数据，N为最大采样时刻。

在计算正传震源波场时，可以根据固定的时间步数存储波场值，存储的波场值用于后续震源波场逆时计算。具体地，在正传震源波场沿时间从最小时刻到最大采样时刻N进行传播的过程中，根据固定时间步数ΔT₁按以下算式存储波场值W_g：

式中，W_g为全局检查点(checkpoint)。

在正传震源波场由N-ΔT₁时刻延拓至N时刻的过程中，即在时间段[N-ΔT₁,N]内，根据固定时间步数ΔT₂按以下算式存储波场值W_l(ΔT₂＜ΔT₁)：

式中，W_l为局部检查点(checkpoint)。

S102：获取观测地震数据；根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源。

在一些实施例中，W₂距离用来表示两个分布的相似程度，经典的基于最小平方度量的全波形反演存在着较为严重的周波跳跃问题，而采用W₂距离计算全波形反演的距离目标函数具有更加良好的凸性，能够有效抑制周波跳跃现象。本说明书中，可以采用卷积、非负化、归一化操作对初始的模拟地震数据与观测地震数据进行预处理，通过构造卷积型W₂距离目标函数，以对初始的模拟地震数据和观测地震数据的差异进行对比，得到较为可靠的卷积型W₂距离目标函数值。

在一些实施例中，伴随波场的震源称作伴随源。

在一些实施例中，上述根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源，具体实施时，可以包括：

S1：根据初始的模拟地震数据，计算归一化后的模拟地震数据；

S2：根据归一化后模拟地震数据，计算反函数；

S3：根据反函数、归一化后的模拟地震数据、观测地震数据、卷积型W₂距离目标函数，计算卷积型W₂距离目标函数值；

S4：根据卷积型W₂距离目标函数值，计算伴随源。

在一些实施例中，上述根据初始的模拟地震数据，计算归一化后的模拟地震数据，具体实施时，可以包括：

S1：根据初始的模拟地震数据，计算卷积后的模拟地震数据；

S2：根据卷积后的模拟地震数据，计算非负化后的模拟地震数据；

S3：根据非负化后的模拟地震数据，计算归一化后的模拟地震数据；

在一些实施例中，上述根据初始的模拟地震数据，计算卷积后的模拟地震数据，包括：

按照以下算式计算卷积后的模拟地震数据：

式中，f_k(k＝1,2,...,i)为初始的模拟地震数据，

为从多道观测地震数据中选取的参考道，

为卷积后的模拟地震数据。

在一些实施例中，上述根据卷积后的模拟地震数据，计算非负化后的模拟地震数据，包括：

按照以下算式计算非负化后的模拟地震数据：

式中，c为非负化系数，d_min为观测地震数据的最小值，

为非负化后的模拟地震数据。

在一些实施例中，上述根据非负化后的模拟地震数据，计算归一化后的模拟地震数据，包括：

按照以下算式计算归一化后的模拟地震数据：

式中，

表示归一化后的模拟地震数据。

在一些实施例中，上述根据归一化后的模拟地震数据，计算反函数，具体操作时，可以包括：

S1：根据观测地震数据，计算卷积后的观测地震数据；

S2：根据卷积后的观测地震数据，计算非负化后的观测地震数据；

S3：根据非负化后的观测地震数据，计算归一化后的观测地震数据；

S4：根据归一化后的观测地震数据和采样间隔，构造第一中间函数；

S5：根据归一化后的模拟地震数据和采样间隔，构造第二中间函数；

S6：根据第二中间函数，利用二分法确定第二中间函数所在的值域区间；并根据第二中间函数所在的值域区间确定第一中间函数的反函数所在的时间区间；

S7：根据第二中间函数所在的值域区间、第一中间函数的反函数所在的时间区间、第二中间函数，计算第一中间函数的反函数；

S8：根据第一中间函数的反函数，计算插值后的观测地震数据。

其中，上述第一中间函数可以记为G_i；上述第二中间函数可以记为F_i；上述第二中间函数所在的值域区间可以记为[G_k,G_k+1]；上述第一中间函数的反函数可以记为G_i ^-1(F_i)；上述第一中间函数的反函数所在的时间区间可以记为[t_k,t_k+1]。

在一些实施例中，上述根据观测地震数据，计算卷积后的观测地震数据，包括：

根据以下算式计算卷积后的观测地震数据：

式中，d_k(k＝1,2,...,i)为观测地震数据，

为从初始的模拟地震数据中选取的参考道，

为卷积后的观测地震数据。

在一些实施例中，上述根据卷积后的观测地震数据，计算非负化后的观测地震数据，包括：

按照以下算式计算非负化后的观测地震数据：

式中，

为非负化后的观测地震数据。

在一些实施例中，上述非负化后的观测地震数据，计算归一化后的观测地震数据，包括：

按照以下算式计算归一化后的观测地震数据：

式中，

为归一化后的观测地震数据。

在一些实施例中，上述根据归一化后的观测地震数据和采样间隔，构造第一中间函数，包括：

按照以下算式构造第一中间函数：

在一些实施例中，上述根据归一化后的模拟地震数据和采样间隔，构造第二中间函数，包括：

按照以下算式构造第二中间函数：

在一些实施例中，上述根据第二中间函数所在的值域区间、第一中间函数的反函数所在的时间区间、第二中间函数，计算第一中间函数的反函数，包括：

按照以下算式计算第一中间函数的反函数：

式中，G_i ^-1(F_i)为G_i的反函数，它表示函数G_i取值为F_i时所对应的时刻。

在一些实施例中，上述根据第一中间函数的反函数，计算插值后的观测地震数据，包括：

按照如下公式计算插值后的观测地震数据：

式中，

表示向下取整后的第n_i个采样点，

表示插值后的观测数据，插值后的观测地震数据后续用于计算伴随源。

在一些实施例中，上述根据反函数、归一化后的模拟地震数据、观测地震数据、卷积型W₂距离目标函数，计算卷积型W₂距离目标函数值，包括：

按照以算式计算W₂距离目标函数值：

式中，f＝{f_i}(i＝1,2,...,N)为初始的模拟地震数据，d＝{d_i}(i＝1,2,...,N)为观测地震数据，为方便说明，在此处仅考虑单道数据。

在一些实施例中，上述根据卷积型W₂距离目标函数值，计算伴随源，具体实施时，可以包括：

不同的目标函数会产生不同的伴随源。本发明针对卷积型W₂距离目标函数推导新的伴随源表达式。为此，推导卷积型W₂距离目标函数对于当前速度模型v_j的梯度：

式中，

为Fréchet(泛函数)导数，g_j表示空间位置为j处的梯度，

为伴随源，即伴随源等于卷积型W₂距离目标函数对初始的模拟地震数据的导数。具体地，可以将g_j写为如下形式：

分别对以上两项进行整理，对于第一项有：

将式(24)写为矩阵-向量形式：

类似地，对于第二项有：

将式(26)写为矩阵-向量形式：

根据式(22)至式(27)可知，梯度g_j可以表示为：

式中，r＝{r_i}(i＝1,2,...,N)为使用卷积型W₂距离目标函数所对应的伴随源，其表达式为：

式(29)中各项为：

K＝diag(t_i-G^-1(F_i)),(i＝1,2,...,N) (33)

S103：根据所述伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场。

在一些实施例中，上述根据所述伴随源计算伴随波场，包括：

根据以下算式计算伴随波场：

式中，

为伴随波场，x_r表示检波点坐标。

参见图4，在一些实施例中，上述基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场，具体实施时，可以包括：

按照以下方式进行当前轮的震源波场逆时重构：

S1：获取与分裂式完美匹配层边界相关的上一轮的关联卷积项；

S2：根据局部检查点值，从最大时刻逆时计算PML区域震源波场，得到PML区域震源波场逆时间序列上的结果；

S3：根据PML区域震源波场逆时间序列上的结果，计算得到重构震源波场逆时间序列上的结果；

S4：根据局部检查点值，更新与分裂式完美匹配层边界相关的关联卷积项，得到当前轮的关联卷积项；

S5：检测当前时间段内逆时间序列上的重构震源波场是否计算完毕，以确定是否触发进行下一轮的震源波场逆时重构。

在一些实施例中，根据局部检查点值，从最大时刻逆时计算PML区域震源波场，得到PML区域震源波场逆时间序列上的结果，包括：

按照以下方式，计算得到逆时间序列上的结果中的当前逆时间序列PML区域震源波场数据值：

S1：在计算出当前逆时间序列PML区域震源波场数据值后，将当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点进行比较；

S2：在确定当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点相同的情况下，使用局部检查点中该时间点对应的PML区域震源波场数据值替换当前逆时间序列PML区域震源波场数据值；

S3：在确定当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点不相同的情况下，利用该时刻之前最近的局部检查点处的PML区域震源波场逆时间序列计算该时刻的PML区域震源波场数据值。

在一些实施例中，上述根据局部检查点值，从最大时刻逆时计算PML区域震源波场，得到PML区域震源波场逆时间序列上的结果，包括：

按照以下算式计算PML区域震源波场逆时间序列上的结果：

在一些实施例中，上述根据PML区域震源波场逆时间序列上的结果，计算得到重构震源波场逆时间序列上的结果，包括：

按照以下算式得到重构震源波场逆时间序列上的结果：

分别表示n+1时刻、n时刻、n-1时刻的重构震源波场。

在一些实施例中，上述根据局部检查点值，更新与分裂式完美匹配层边界相关的关联卷积项，得到当前轮的关联卷积项，包括：

按照以下算式计算当前轮的关联卷积项：

因为逆时计算的PML区域震源波场是不稳定的，所以要利用局部检查点存储的PML区域震源波场替换逆时重构的PML区域震源波场，即在时间段[N-1-ΔT₁,N-1]内，如果当前时刻为局部检查点所在时刻，将局部检查点存储的PML区域震源波场直接赋给逆时计算的PML区域震源波场。进行替换的意义在于，保证在逆时重构过程中，PML区域震源波场出现不稳定前用准确的PML区域震源波场替换掉不稳定的PML区域震源波场，避免不稳定问题。如前所述，本说明书仅在最后的时间段[N-1-ΔT₁,N-1]内存储了n_l个局部检查点，当逆时传播越过该时间段，局部检查点将处于空闲状态，如果不对它们进行更新，后续逆时重构仍然会出现不稳定。为此，当逆时重构到达N-1-ΔT₁时刻，利用前一个全局检查点存储的数据值作为初始值，重新启动正传震源波场的计算，即从N-1-2ΔT₁计算到N-1-ΔT₁时刻，此时在时间段[N-1-2ΔT₁,N-1-ΔT]内以ΔT₂为间隔重新存储局部检查点的数据，这样就实现了n_l个局部检查点的循环利用。当重启的正传震源波场计算到达N-1-ΔT₁时刻，则继续伴随波场的逆时计算和震源波场的逆时重构，直至最小时刻。

S104：根据所述伴随波场和所述逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度。

全波形反演更新速度模型需要计算梯度，可以采用伴随状态法计算重构震源波场和伴随波场的零延迟互相关求取梯度。

具体的，可以按照以下算式采用伴随状态法计算重构震源波场和伴随波场的零延迟互相关求取梯度：

S105：根据卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。

在一些实施例中，上述根据卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型，具体实施时，包括：

S1：根据卷积型W₂距离目标函数的梯度，确定速度模型的更新方向；并计算出相应的更新步长；

S2：根据所述更新方向和更新步长，更新当前速度模型。

在一些实施例中，上述根据卷积型W₂距离目标函数的梯度，确定速度模型的更新方向；并计算出相应的更新步长，包括：

S1：利用拟牛顿算法计算目标速度模型的更新方向；

S2：利用wolfe准则确定更新步长；

在一些实施例中，wolfe准则是一种非精确线搜索方法，它在每一次迭代中不要求目标函数达到精确最小，只需要其有一定的下降量，虽然迭代次数有所增加，但是收敛速度较快。

在一些实施例中，上述根据所述更新方向和更新步长，更新当前速度模型，包括：

按照以下算式计算更新后的目标速度模型：

v_l+1＝v_l+α_lΔv_l (40)

式中，l为迭代次数，Δv_l为目标速度模型更新方向，α_l为满足wolfe条件的更新步长，v_l为更新前的目标速度模型，v_l+1为更新后的目标速度模型。

在一些实施例中，依据式(40)在每一次循环计算过程中均对目标速度模型进行更新，直到满足迭代终止条件。上述迭代终止条件，可以为超过设定的最大迭代次数和卷积型W₂距离目标函数下降率小于设定比例。

在一些实施例中，在得到目标速度模型之后，可以利用该目标速度模型对目标区域进行地震成像。

在一些实施例中，上述计算过程均在GPU并行计算平台上实现，它可以针对大量数据进行复杂的密集同类计算。

在一个具体的场景示例中，针对一维时移数据分别采用传统型L₂距离目标函数、卷积型L₂距离目标函数、卷积型W₂距离目标函数进行处理，得到实验结果。其中，图2a表示时移子波和准确子波的波形，虚线表示时移子波，实线表示准确子波。图2b表示分别采用传统型L₂距离目标函数、卷积型L₂距离目标函数、卷积型W₂距离目标函数对时移子波进行计算后得到的凸性对比结果，由图2b可知，卷积型W₂距离目标函数对时移子波进行计算后能得到较好的结果。

在一个具体的场景示例中，针对一维信号时移量分别采用传统型L₂距离目标函数、卷积型L₂距离目标函数、卷积型W₂距离目标函数进行处理，得到实验结果。其中，图3a表示采用加常数非负化预处理方法对一维信号时移量的凸性分析结果，图3b表示采用幂指数非负化预处理方法对一维信号时移量的凸性分析结果，图3c表示采用平方非负化预处理方法对一维信号时移量的凸性分析结果。由图3a、图3b、图3c可知，采用加常数非负化预处理方法能产生较好的结果。

在一个具体的场景示例中，使用Marmousi模型比较传统卷积型L₂距离目标函数与卷积型W₂距离目标函数的反演效果。设定准确模型网格大小为234m×663m，网格间距为10m。设置共60炮的双边变排列观测系统，炮间距100m，最大炮检距4000m。准确的震源子波为雷克子波的一阶导数，记录时长为3s，采样间隔为0.8ms。用于反演的初始速度模型深度小于60m的部分为真实模型，随后速度随深度线性增加。图5a、图5b分别为准确的Marmousi速度模型和用于反演的初始速度模型，该初始模型只保留了准确模型纵向速度变化趋势。在反演过程中计算预测数据的子波为雷克子波，两种震源子波的主频相同。首先从主频为3Hz子波进行反演，图6a、图6b分别为基于卷积型L₂距离目标函数和卷积型W₂距离目标函数全波形反演的第一次迭代梯度。在初始模型较差的情况下，图6a、图6b表明基于卷积型W₂距离目标函数计算的梯度包含更丰富的低波数更新，这说明由卷积型W₂距离目标函数计算的伴随源的低频成分较传统卷积型L₂距离目标函数更丰富。随着迭代次数的增加，基于卷积型L₂距离目标函数计算的梯度呈现越来越多的高波数成分，在第8次迭代时陷入局部极小，如图7a所示。相比之下，基于卷积型W₂距离目标函数反演得到的速度模型没有出现明显的局部极值问题，如图7b所示。图7c为反演最后频带数据卷积型W₂距离目标函数的结果，该结果与真实模型十分接近，证实了卷积型W₂距离目标函数的稳健性。在本场景示例中对多频带数据依次反演实现多尺度反演，第一个频带数据主频为3Hz，迭代30次，第二个频带主频5Hz，迭代30次，第三个频带数据的主频为10Hz，迭代15次。

在一个具体的场景示例中，使用Overthrust模型证实本说明书提出的卷积型W₂距离目标函数的稳健性，其中图8a和图8b分别为准确的Overthrust速度模型和用于反演的初始速度模型。图9a、图9b为基于准确速度模型计算得到的观测地震数据，图9a是无噪观测地震数据，图9b是含随机噪声数据(信噪比1dB)的观测地震数据。由于信噪比低，有效的反射波形几乎被噪声全部淹没。尽管如此，基于卷积型W₂距离目标函数的全波形反演方法在迭代40次后，不含噪和含噪数据(信噪比分别为1dB、2dB、4dB和20dB)的反演结果几乎完全一致，如图10a、图10b、图10c、图10d、图10e所示，这说明卷积型W₂距离目标函数具有很好的抗随机噪声能力。提取不同信噪比地震数据反演结果的垂直速度曲线并在图11a、图11b中对比，通过曲线的高度重合，进一步证实了随机噪声对反演结果的影响微弱。图12a、图12b、图12c、图12d、图12e对比了反演不同信噪比数据时卷积型W₂距离目标函数的下降趋势，由图12a、图12b、图12c、图12d、图12e可知在不含噪和含噪情况下卷积型W₂距离目标函数的下降速度均较快，这也说明卷积型W₂距离目标函数具有良好的抗噪性。图13进一步对比了含噪情况下反演的速度模型与准确速度模型之间的均方根误差，对比表明，即使观测地震数据受随机噪声严重污染(信噪比1dB)，其对应反演结果的均方根误差也小于0.8％。随着信噪比提高，反演结果的准确度随之提高。

虽然本说明书提供了如下述实施例或附图14所示的方法操作步骤或装置结构，但基于常规或者无需创造性的劳动在所述方法或装置中可以包括更多或者部分合并后更少的操作步骤或模块单元。在逻辑性上不存在必要因果关系的步骤或结构中，这些步骤的执行顺序或装置的模块结构不限于本说明书实施例或附图所示的执行顺序或模块结构。所述的方法或模块结构的在实际中的装置、服务器或终端产品应用时，可以按照实施例或者附图所示的方法或模块结构进行顺序执行或者并行执行(例如并行处理器或者多线程处理的环境、甚至包括分布式处理、服务器集群的实施环境)。

基于上述基于卷积型W₂距离目标函数的全波形反演方法，本说明书还提出基于卷积型W₂距离目标函数的全波形反演装置的实施例。如图14所示，所述基于卷积型W₂距离目标函数的全波形反演装置具体包括以下模块：

第一处理模块1401，用于获取初始速度模型，根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据；

计算模块1402，用于获取观测地震数据；根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源；

第二处理模块1403，用于根据所述伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；

求取模块1404，用于根据所述伴随波场和所述逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；

更新模块1405，用于根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。

在一些实施例中，上述第一处理模块1401具体可以用于在计算区域内，根据初始速度模型，构建震源波动方程；在计算区域外，加载分裂式完美匹配层，构建分裂式完美匹配层波动方程；根据所述震源波动方程、所述分裂式完美匹配层波动方程，计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据。

在一些实施例中，上述计算模块1402具体可以用于获取观测地震数据；根据初始的模拟地震数据，计算归一化后的模拟地震数据；根据归一化后的模拟地震数据，计算反函数；根据反函数、归一化后的模拟地震数据、观测地震数据、卷积型W₂距离目标函数，计算卷积型W₂距离目标函数值；根据卷积型W₂距离目标函数值，计算伴随源。

在一些实施例中，上述第二处理模块1403具体可以用于根据伴随源计算伴随波场；并按照以下方式计算当前轮的正传震源波场：获取与分裂式完美匹配层边界相关的上一轮的关联卷积项；在计算区域外，根据与分裂式完美匹配层边界相关的上一轮的关联卷积项和分裂式完美匹配层波动方程，得到PML区域震源波场；在计算区域内，根据震源波动方程和PML区域震源波场，得到当前轮的正传震源波场；更新与分裂式完美匹配层边界相关的关联卷积项，得到当前轮的关联卷积项；检测正时间序列上的正传震源波场是否计算完毕，如果没有计算完毕，计算下一轮的正传震源波场；如果计算完毕，根据正传震源波场确定出初始的模拟地震数据。

在一些实施例中，上述求取模块1404具体可以用于依据伴随状态法计算逆时间序列上的重构震源波场和伴随波场的零延迟互相关求取梯度。

在一些实施例中，上述更新模块1405具体可以用于根据卷积型W₂距离目标函数的梯度，确定速度模型的更新方向；并计算出相应的更新步长；根据所述更新方向和更新步长，更新当前速度模型；每一次循环计算过程中均对目标速度模型进行更新，直到满足迭代终止条件。

本说明书实施例还提供了一种基于上述基于卷积型W₂距离目标函数的全波形反演方法的计算机存储介质，所述计算机存储介质存储有计算机程序指令，在所述计算机程序指令被执行时实现：获取初始速度模型，根据初始速度模型计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据；获取观测地震数据；根据初始的模拟地震数据、观测地震数据和卷积型W₂距离目标函数，计算伴随源；根据伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；根据伴随波场和逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；根据卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。

在本实施例中，上述存储介质包括但不限于随机存取存储器(RandomAccessMemory,RAM)、只读存储器(Read-Only Memory,ROM)、缓存(Cache)、硬盘(Hard DiskDrive,HDD)或者存储卡(Memory Card)。所述存储器可以用于存储计算机程序指令。网络通信单元可以是依照通信协议规定的标准设置的，用于进行网络连接通信的接口。

虽然本说明书提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的手段可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的装置或客户端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境，甚至为分布式数据处理环境)。术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、产品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、产品或者设备所固有的要素。在没有更多限制的情况下，并不排除在包括所述要素的过程、方法、产品或者设备中还存在另外的相同或等同要素。第一，第二等词语用来表示名称，而并不表示任何特定的顺序。

本说明书可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构、类等等。也可以在分布式计算环境中实践本说明书，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。

通过以上的实施例的描述可知，本领域的技术人员可以清楚地了解到本说明书可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本说明书的技术方案本质上可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，移动终端，服务器，或者网络设备等)执行本说明书各个实施例或者实施例的某些部分所述的方法。

虽然通过实施例描绘了本说明书，本领域普通技术人员知道，本说明书有许多变形和变化而不脱离本说明书的精神，希望所附的权利要求包括这些变形和变化而不脱离本说明书的精神。

Claims (10)

1.一种基于卷积型W₂距离目标函数的全波形反演方法，其特征在于，包括：

获取初始速度模型，根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据；

获取观测地震数据；根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源；

根据所述伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；

根据所述伴随波场和所述逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；

根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。

2.根据权利要求1所述的方法，其特征在于，根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据，包括：

在计算区域内，根据初始速度模型，构建震源波动方程；

在计算区域外，加载分裂式完美匹配层，构建分裂式完美匹配层波动方程；

根据所述震源波动方程、所述分裂式完美匹配层波动方程，计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据。

3.根据权利要求2所述的方法，其特征在于，根据所述震源波动方程、所述分裂式完美匹配层波动方程，计算正传震源波场，根据正传震源波场确定出初始的模拟地震数据，包括：

按照以下方式计算当前轮的正传震源波场：

获取与分裂式完美匹配层边界相关的上一轮的关联卷积项；

在计算区域外，根据与分裂式完美匹配层边界相关的上一轮的关联卷积项和分裂式完美匹配层波动方程，得到PML区域震源波场；

在计算区域内，根据震源波动方程和PML区域震源波场，得到当前轮的正传震源波场；

更新与分裂式完美匹配层边界相关的关联卷积项，得到当前轮的关联卷积项；

检测正时间序列上的正传震源波场是否计算完毕，如果没有计算完毕，计算下一轮的正传震源波场；如果计算完毕，根据正传震源波场确定出初始的模拟地震数据。

4.根据权利要求1所述的方法，其特征在于，根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源，包括：

根据初始的模拟地震数据，计算归一化后的模拟地震数据；

根据归一化后的模拟地震数据，计算反函数；

根据反函数、归一化后的模拟地震数据、观测地震数据、卷积型W₂距离目标函数，计算卷积型W₂距离目标函数值；

根据卷积型W₂距离目标函数值，计算伴随源。

5.根据权利要求4所述的方法，其特征在于，根据反函数、归一化后的模拟地震数据、观测地震数据、卷积型W₂距离目标函数，计算卷积型W₂距离目标函数值，包括：

按照以下算式，计算卷积型W₂距离目标函数值：

其中，W₂为卷积型W₂距离目标函数值，f＝{f_i}(i＝1,2,...,N)为初始的模拟地震数据，d＝{d_i}(i＝1,2,...,N)为观测地震数据，N为最大采样时刻，M为模型的空间网格数，Δt为采样间隔，t_i＝(i-1)Δt(i＝1,2,...,N)为第i个采样点，G_i ^-1表示求反函数，G_i ^-1(F_i)表示函数G_i取值为F_i时所对应的时刻，

为归一化后的观测地震数据，

为归一化后的模拟地震数据。

6.根据权利要求1所述的方法，其特征在于，基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场，包括：

按照以下方式进行当前轮的震源波场逆时重构：

获取与分裂式完美匹配层边界相关的上一轮的关联卷积项；

根据局部检查点值，从最大时刻逆时计算PML区域震源波场，得到PML区域震源波场逆时间序列上的结果；

根据PML区域震源波场逆时间序列上的结果，计算得到重构震源波场逆时间序列上的结果；

根据局部检查点值，更新与分裂式完美匹配层边界相关的关联卷积项，得到当前轮的关联卷积项；

检测当前时间段内逆时间序列上的重构震源波场是否计算完毕，以确定是否触发进行下一轮的震源波场逆时重构。

7.根据权利要求6所述的方法，其特征在于，根据局部检查点值，从最大时刻逆时计算PML区域震源波场，得到PML区域震源波场逆时间序列上的结果，包括：

按照以下方式，计算得到逆时间序列上的结果中的当前逆时间序列PML区域震源波场数据值：

在计算出当前逆时间序列PML区域震源波场数据值后，将当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点进行比较；

在确定当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点相同的情况下，使用局部检查点中该时间点对应的PML区域震源波场数据值替换当前逆时间序列PML区域震源波场数据值；

在确定当前逆时间序列PML区域震源波场数据值所对应的时间点与局部检查点所在时间点不相同的情况下，利用该时刻之前最近的局部检查点处的PML区域震源波场逆时间序列计算该时刻的PML区域震源波场数据值。

8.根据权利要求1所述的方法，其特征在于，根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，包括：

根据卷积型W₂距离目标函数的梯度，确定速度模型的更新方向；并计算出相应的更新步长；

根据所述更新方向和更新步长，更新当前速度模型。

9.一种基于卷积型W₂距离目标函数的全波形反演装置，其特征在于，包括：

第一处理模块，用于获取初始速度模型，根据所述初始速度模型计算正传震源波场，根据所述正传震源波场确定出初始的模拟地震数据；

计算模块，用于获取观测地震数据；根据所述初始的模拟地震数据、所述观测地震数据和卷积型W₂距离目标函数，计算伴随源；

第二处理模块，用于根据所述伴随源计算伴随波场；并基于双层检查点进行震源波场逆时重构，得到逆时间序列上的重构震源波场；

求取模块，用于根据所述伴随波场和所述逆时间序列上的重构震源波场，求取卷积型W₂距离目标函数的梯度；

更新模块，用于根据所述卷积型W₂距离目标函数的梯度，更新当前速度模型，得到目标速度模型。

10.一种计算机可读存储介质，其特征在于，其上存储有计算机指令，所述指令被处理器执行时实现权利要求1至8中任一项所述方法的步骤。

Images