一种基于局部存储策略的声波全波形反演方法及装置

1.本说明书属于地震探勘技术领域，尤其涉及一种基于局部存储策略的声波全波形反演方法及装置。


背景技术：

2.在现有技术中，通常采用全波形反演法确定地震波的纵波波速。全波形反演是一种获取高分辨率地下速度模型的方法，它使用有限差分技术进行波场模拟，利用速度模型得到模拟波场，并基于观测记录(观测记录来源于观测波场)与模拟记录(模拟记录来源于模拟波场)的目标函数，确定观测波场与模拟波场的差异，在差异较大的情况下利用优化算法更新速度模型，利用更新后的速度模型调整模拟波场，直至模拟记录与观测记录接近，也就是使得模拟波场接近观测波场，从而获得较高精度的速度模型。现有的全波形反演方法需要计算并存储全部时刻的模拟波场，计算机的数据处理能力与数据传输能力难以满足全波形反演的需求，因此现有的全波形反演计算效率低、内存占用大。
3.针对上述技术问题，目前尚未提出有效的解决方案。


技术实现要素：

4.本说明书提供了一种基于局部存储策略的声波全波形反演方法及装置，能够高效进行全波形反演，减少数据处理量，并最终确定高精度的速度模型。
5.本说明书实施例的目的是提供一种基于局部存储策略的声波全波形反演方法，包括：
6.获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；
7.根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；
8.根据所述观测记录和所述模拟记录，计算初始l2目标函数值；
9.检测所述初始l2目标函数值是否大于预设值；
10.在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
11.进一步地，所述方法的另一个实施例中，还包括：在确定所述初始l2目标函数值小于或等于预设值的情况下，利用所述初始速度模型进行地震偏移成像。
12.进一步地，所述方法的另一个实施例中，所述根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出目标时刻的模拟波场，包括：
13.根据所述地震波参数和所述初始速度模型，对所述目标区域进行模拟波场计算，确定全局模拟波场以及地震波到达时刻；
14.根据所述地震波参数，确定采样时间步长、采样点数量；
15.根据所述采样时间步长、所述采样点数量、所述地震波到达时刻，确定目标时刻；
16.根据所述目标时刻，从所述全局模拟波场中，提取出所述目标时刻的模拟波场。
17.进一步地，所述方法的另一个实施例中，所述根据所述采样时间步长、所述采样点数量、所述地震波到达时刻，确定目标时刻，包括：
18.按照如下算式确定目标时刻：
19.t0(x,y,z)+t
ny
(x,y,z)＝t0+a(x,y,z)
·
ts(a(x,y,z)＝1,2,,an(x,y,z))
20.其中，t0(x,y,z)表示网格点x(x,y,z)的地震波到达时刻，t
ny
(x,y,z)表示网格点x(x,y,z)的采样时刻，t0(x,y,z)+t
ny
(x,y,z)表示网格点x(x,y,z)的目标时刻，a(x,y,z)表示网格点x(x,y,z)的采样次数，ts表示采样时间步长，an(x,y,z)表示网格点x(x,y,z)的采样点数量，x表示网格点的x轴坐标值，y表示网格点的y轴坐标值，z表示网格点的z轴坐标值。
21.进一步地，所述方法的另一个实施例中，所述根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型，包括：
22.根据所述模拟记录和所述观测记录的差异，计算当前轮次的伴随波场；
23.根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度；
24.根据所述当前轮次的l2目标函数梯度，利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
25.进一步地，所述方法的另一个实施例中，在得到当前轮次的目标速度模型之后，所述方法还包括：
26.根据所述地震波参数和所述当前轮次的目标速度模型，确定下一轮次的全局模拟波场；并从所述下一轮次的全局模拟波场中，确定出下一轮次的目标时刻的模拟波场、下一轮次的模拟记录；
27.根据所述观测记录和所述下一轮次的模拟记录，计算下一轮次的l2目标函数值；
28.检测所述下一轮次的l2目标函数值是否大于预设值；
29.在确定所述下一轮次的l2目标函数值大于预设值的情况下，根据所述下一轮次的目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述当前轮次的目标速度模型进行更新，得到下一轮次的目标速度模型。
30.进一步地，所述方法的另一个实施例中，所述根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度，包括：
31.按照如下算式计算当前轮次的l2目标函数梯度：
[0032][0033]
其中，t0表示地震波到达时刻，t
ny
(x,y,z)表示采样时刻，表示s2函数对体积模量的梯度，l
p
表示伴随波场，v
x
表示x轴方向的网格点振动速度分量，vy表示表示y轴方
向的网格点振动速度分量，vz表示z轴方向的网格点振动速度分量，t表示时间，表示s2函数对纵波波速的梯度，表示体积模量对纵波波速的梯度，ρ表示密度，v表示纵波波速，k表示体积模量，x、y、z分别为网格点的x轴坐标值、y轴坐标值、z轴坐标值，表示网格点振动速度分量的散度，gn表示l2目标函数对纵波波速的梯度，简称l2目标函数梯度，n表示当前轮次。
[0034]
进一步地，所述方法的另一个实施例中，所述根据当前轮次的l2目标函数梯度，利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型，包括：
[0035]
根据当前轮次的l2目标函数梯度，确定当前轮次的更新方向；
[0036]
根据当前轮次的更新方向，对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0037]
另一方面，本说明书实施例还提供了一种基于局部存储策略的声波全波形反演装置，包括：
[0038]
获取模块，用于获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；
[0039]
第一计算模块，用于根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；
[0040]
第二计算模块，用于根据所述观测记录和所述模拟记录，计算初始l2目标函数值；
[0041]
检测模块，用于检测所述初始l2目标函数值是否大于预设值；
[0042]
更新模块，用于在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0043]
再一方面，本说明书实施例还提供了一种计算机可读存储介质，其上存储有计算机指令，所述计算机可读存储介质执行所述指令时实现上述基于局部存储策略的声波全波形反演方法。
[0044]
本说明书实施例提供的一种基于局部存储策略的声波全波形反演方法，通过获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；根据所述观测记录和所述模拟记录，计算初始l2目标函数值；检测所述初始l2目标函数值是否大于预设值；在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0045]
进一步，可以通过如下方式，得到当前轮次的目标速度模型：根据所述模拟记录和所述观测记录的差异，计算当前轮次的伴随波场；根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度；根据所述当前轮次的l2目标函数梯度，利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
附图说明
[0046]
为了更清楚地说明本说明书实施例，下面将对实施例中所需要使用的附图作简单地介绍，下面描述中的附图仅仅是本说明书中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0047]
图1是本说明书提供的一种基于局部存储策略的声波全波形反演方法一个实施例的流程示意图；
[0048]
图2是本说明书提供的一个具体场景示例中输出目标速度模型的流程示意图；
[0049]
图3a是现有技术中的全时刻互相关方法示意图；
[0050]
图3b是现有技术中的激发振幅时刻互相关方法示意图；
[0051]
图3c是本说明书提供的进行目标时刻互相关计算方法示意图；
[0052]
图4a是基于二维overthrust模型的真实速度模型结果图；
[0053]
图4b是基于二维overthrust模型的初始速度模型结果图；
[0054]
图5a是现有技术中的全时刻互相关方法第一次计算l2目标函数梯度的结果图；
[0055]
图5b是基于本技术提供的目标时刻互相关计算方法第一次计算l2目标函数梯度的结果图；
[0056]
图5c是使用激发振幅时刻互相关方法第一次计算l2目标函数梯度的结果图；
[0057]
图6a是ts＝50ms所得到的l2目标函数梯度的结果图；
[0058]
图6b是ts＝100ms所得到的l2目标函数梯度的结果图；
[0059]
图6c是ts＝150ms所得到的l2目标函数梯度的结果图；
[0060]
图6d是ts＝200ms所得到的l2目标函数梯度的结果图；
[0061]
图7a是基于本技术提供的起跳阈值方法得到的地震波到达时间场结果图；
[0062]
图7b是基于现有技术的激发振幅方法得到的地震波到达时间场结果图；
[0063]
图8a是使用现有技术中的全时刻互相关方法以及无噪声数据进行全波形反演的结果图；
[0064]
图8b是基于本技术提供的目标时刻互相关计算方法以及无噪声数据进行全波形反演的结果图；
[0065]
图8c是使用现有技术中的激发振幅时刻互相关方法以及无噪声数据进行全波形反演的结果图；
[0066]
图9a是使用现有技术中的全时刻互相关方法以及含噪声数据进行全波形反演的结果图；
[0067]
图9b是基于本技术提供的目标时刻互相关计算方法以及含噪声数据进行全波形反演的结果图；
[0068]
图9c是使用现有技术中的激发振幅时刻互相关方法以及含噪声数据进行全波形反演的结果图；
[0069]
图10a是基于三维overthrust模型得到的真实速度模型结果图；
[0070]
图10b是基于三维overthrust模型得到的初始速度模型结果图；
[0071]
图11a是基于三维overthrust模型得到的初始速度模型采用现有技术中的全时刻互相关方法进行全波形反演的结果图；
[0072]
图11b基于三维overthrust模型得到的初始速度模型采用本技术提供的目标时刻
互相关方法进行全波形反演的结果图；
[0073]
图12是本说明书提供的一种基于局部存储策略的声波全波形反演装置一个实施例的结构示意图；
[0074]
图13是本说明书提供的一种服务器的结构组成示意图。
具体实施方式
[0075]
为了使本技术领域的人员更好地理解本说明书中的技术方案，下面将结合本说明书实施例中的附图，对本说明书实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本说明书一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本说明书保护的范围。
[0076]
在现有技术中，通常采用全波形反演法确定地震波的纵波波速。全波形反演是一种获取高分辨率地下速度模型的方法，它使用有限差分技术进行波场模拟，利用速度模型得到模拟波场，并基于观测记录(观测记录来源于观测波场)与模拟记录(模拟记录来源于模拟波场)的目标函数，确定观测波场与模拟波场的差异，在差异较大的情况下利用优化算法更新速度模型，利用更新后的速度模型调整模拟波场，直至模拟记录与观测记录接近，也就是使得模拟波场接近观测波场，从而获得较高精度的速度模型。
[0077]
考虑到现有的全波形反演方法，需要计算并存储全部时刻的模拟波场，然后利用全部时刻的模拟波场进行互相关计算求取梯度，根据梯度对速度模型进行更新；并且，大量的模拟波场数据需要存储在硬盘中。在具体计算时，采用gpu(图形处理器)和mpi(信息传递接口)并行计算方法，利用gpu的分块计算方式提高计算速率，并通过mpi多gpu并行同时计算多炮的梯度，以达到节省反演时间的目的。然而，gpu的存储容量小，无法用于存储数据，并且也无法直接与硬盘进行数据交互。在每一次的互相关计算中，大量的模拟波场数据必须先传输至cpu，再从cpu传输至gpu，进行全波形反演的计算，得到调整后的速度模型与模拟波场数据；调整后的速度模型与模拟波场数据再经过cpu传输至硬盘中进行存储。
[0078]
在现有技术中，针对三维的目标区域，所需要的模拟波场数据的个数为n
x
·
ny·
nz·nt
，n
x
表示目标区域x方向上的网格点数量，ny表示目标区域y方向上的网格点数量，nz表示目标区域z方向上的网格点数量，n
t
表示时间采样点的数量，因此现有方法中所需要的模拟波场数据的数据量是巨大的。因此，计算机的数据处理能力与数据传输能力难以满足全波形反演高速运算的需求，现有的全波形反演存在计算效率低、内存占用大的问题。
[0079]
针对现有方法存在的上述问题以及产生上述问题的具体原因，本技术引入基于局部存储策略的声波全波形反演方法，以减少全波形反演中的数据处理量，提升全波形反演计算效率。
[0080]
基于上述思路，本说明书提出一种基于局部存储策略的声波全波形反演方法。首先，获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出目标时刻的模拟波场、模拟记录；然后，根据所述观测记录和所述模拟记录，计算初始l2目标函数值；检测所述初始l2目标函数值是否大于预设值；最后，在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场利用梯度类优化算法对所述初始速度
模型进行更新，得到当前轮次的目标速度模型。
[0081]
参阅图1所示，本说明书实施例提供了一种基于局部存储策略的声波全波形反演方法。具体实施时，该方法可以包括以下内容。
[0082]
s101：获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数。
[0083]
在一些实施例中，目标区域指的是待全波形反演的地下区域。
[0084]
在一些实施例中，在获取目标区域之后，确定目标区域的三维坐标，并对目标区域进行三维的网格划分，将连续的三维目标区域划分为多个离散的网格点，网格点记为x(x,y,z)，x、y、z分别表示网格点x的x轴坐标值、y轴坐标值、z轴坐标值；并且，针对每个网格点均设定一个振幅阈值p0(x；t)，t表示在模拟波场计算时地震波到达该网格点的时间。
[0085]
在一些实施例中，目标区域的初始速度模型，可以记为m0(v,ρ)，m0是由地震波在目标区域传播的纵波波速v、目标区域的密度ρ构成的函数，为了简化计算，本技术中假定密度是已知的定值，后续也不会对密度进行更新，因此，m0表示地震波在目标区域传播的纵波波速v。
[0086]
在一些实施例中，通过对目标区域进行野外地震实验，在炮点发射出地震波，然后在检波点安置采集仪器，采集到检波点的观测波场(观测波场也就是真实波场)，作为观测记录；其中，检波点的观测波场具体是指：检波点的观测波场的散度。
[0087]
在一些实施例中，可以获取目标区域的地质资料，通过地质资料结合经验公式确定目标区域的初始速度模型与地震波参数；地震波参数包括：子波频率、子波周期；子波指的是能够稳定传播的地震波。例如，可以基于地质资料通过偏移速度分析、层析成像、人工智能模型等方法得到初始速度模型。
[0088]
s102：根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场。
[0089]
在一些实施例中，根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出目标时刻的模拟波场，具体包括：
[0090]
s1：根据所述地震波参数和所述初始速度模型，对所述目标区域进行模拟波场计算，确定全局模拟波场以及地震波到达时刻；
[0091]
s2：根据所述地震波参数，确定采样时间步长、采样点数量；
[0092]
s3：根据所述采样时间步长、所述采样点数量、所述地震波到达时刻，确定目标时刻；
[0093]
s4：根据所述目标时刻，从所述全局模拟波场中，提取出所述目标时刻的模拟波场。
[0094]
在一些实施例中，全局模拟波场指的是：所有网格点在所有计算时刻构成的模拟波场。
[0095]
在一些实施例中，根据所述地震波参数和所述初始速度模型，对所述目标区域进行模拟波场计算，确定全局模拟波场以及地震波到达时刻，具体包括：根据初始速度模型和地震波参数，利用gpu模拟炮点放炮过程进而通过模拟计算得到每个网格点的模拟波场，在模拟波场计算中，按照预设的时间间隔检测网格点的振幅，若满足条件p(x；t)≥p0(x；t)时，则认为地震波到达该网格点，将第一次满足p(x；t)≥p0(x；t)的时间t作为地震波到达
时刻：
[0096]
t0(x)＝t(1)
[0097]
其中，p(x；t)表示网格点x的振幅，t0(x)表示网格点x的地震波到达时刻，t0(x)可以简记为t0。
[0098]
在一些实施例中，通过公式1确定地震波到达时刻的方法也被称为起跳阈值法。
[0099]
通过上述实施例，采用p(x；t)≥p0(x；t)作为地震波到达时刻的判断条件，可以避免现有技术中的激发振幅时刻方法中因为激发振幅时间场计算不准而导致的梯度计算误差。
[0100]
在一些实施例中，可以按照如下算式，计算全局模拟波场：
[0101][0102]
其中，v
x
、vy、vz分别表示沿x轴方向、y轴方向、z轴方向的网格点(网格点视为质点)振动速度分量；p表示压力场数据；f
p
表示应力震源项；k为体积模量，且k满足：k＝ρv2，v表示地震波在目标区域传播的纵波波速，ρ表示密度；t表示时间；x、y、z分别为网格点的x轴坐标值、y轴坐标值、z轴坐标值；表示网格点振动速度分量的散度。
[0103]
通过公式2，可以求出每个网格点振动速度分量的散度，构成了全局模拟波场。
[0104]
在一些实施例中，根据所述地震波参数，确定采样时间步长、采样点数量，具体包括：
[0105]
按照如下公式确定采样时间步长：
[0106][0107]
其中，f
max
表示子波频率的最大值，ts表示采样时间步长；
[0108]
按照如下公式确定采样点数量：
[0109][0110]
其中，an(x,y,z)表示网格点x(x,y,z)对应的采样点数量，an(x,y,z)为正整数，若an(x,y,z)不为正整数则向下取整；tw(x,y,z)表示网格点x(x,y,z)对应的窗口长度，tw(x,y,z)可以取值为子波周期。
[0111]
在一些实施例中，表示的含义是奈奎斯特频率的倒数，公式3的含义是：采样时间步长ts是满足奈奎斯特采样定理的采样间隔。
[0112]
在一些实施例中，当观测波场的信噪比较低，噪声强度较大时，可以适当缩小ts，同时增加采样点数量来增加抗噪性，以取得更好的速度模型计算效果。
[0113]
在一些实施例中，在确定采样时间步长之后，还可以通过如下公式确定采样时刻：
[0114]
t
ny
(x,y,z)＝a(x,y,z)
·
ts(a(x,y,z)＝1,2，，an(x,y,z))(5)
[0115]
其中，t
ny
(x,y,z)表示网格点x对应的采样时刻，可以简记为t
ny
；a(x,y,x)表示采样次数。
[0116]
在一些实施例中，根据所述采样时间步长、所述采样点数量、所述地震波到达时刻，确定目标时刻，具体包括：
[0117]
按照如下算式，确定目标时刻：
[0118]
t0(x,y,z)+t
ny
(x,y,z)＝t0+a(x,y,z)
·
ts(a(x,y,z)＝1,2,,an(x,y,z))(6)
[0119]
其中，t0(x,y,z)表示网格点x(x,y,z)的地震波到达时刻，t
ny
(x,y,z)表示网格点x(x,y,z)的采样时刻，t0(x,y,z)+t
ny
(x,y,z)表示网格点x(x,y,z)的目标时刻，a(x,y,z)表示网格点x(x,y,z)的采样次数，ts表示采样时间步长，an(x,y,z)表示网格点x(x,y,z)的采样点数量，x表示网格点的x轴坐标值，y表示网格点的y轴坐标值，z表示网格点的z轴坐标值。
[0120]
在一些实施例中，a(x,y,z)依次取值为1,2,
…
，an(x,y,z)，相应的，目标时刻t0+t
ny
(x,y,z)分别等于：t0+ts、t0+2ts、
……
、t0+an(x,y,z)
·
ts。
[0121]
在一些实施例中，在确定出确定目标时刻之后，从所有网格点在所有计算时刻构成的全局模拟波场中，提取出所有网格点的目标时刻对应的模拟波场。
[0122]
在一些实施例中，从所有网格点在所有计算时刻构成的全局模拟波场中，提取出检波点对应的全局模拟波场，作为模拟记录。
[0123]
s103：根据所述观测记录和所述模拟记录，计算初始l2目标函数值。
[0124]
在一些实施例中，l2目标函数用于表征模拟波场和观测波场之间的差异；在差异较大的情况下，需要利用优化算法对速度模型进行更新，利用更新后的速度模型得到更加准确、更加接近观测波场的模拟波场。
[0125]
在一些实施例中，通过如下算式，计算初始l2目标函数值：
[0126][0127]
其中，l2(m)表示l2目标函数值，l2目标函数是与m相关的函数，m表示速度模型；t
max
表示模拟波场对应的最大时刻；ω表示目标区域；u
syn
(xr,t；xs)、u
obs
(xr,t；xs)分别表示同一炮、同一位置对应的模拟记录、观测记录，xs表示炮点位置，xr表示检波点位置；δ(x-xr)表示单位脉冲函数。
[0128]
s104：检测所述初始l2目标函数值是否大于预设值。
[0129]
在一些实施例中，预设值可以设定为10-3
。在初始l2目标函数值大于预设值的情况下，说明观测记录和模拟记录之间的差异较大，初始速度模型精度不符合要求，需要对初始速度模型进行更新；在初始l2目标函数值小于或等于预设值的情况下，说明观测记录和模拟记录之间的差异较小，初始速度模型精度符合要求，可以使用初始速度模型进行地震偏移成像。
[0130]
s105：在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的
模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0131]
在一些实施例中，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型，具体包括：
[0132]
s1：根据所述模拟记录和所述观测记录的差异，计算当前轮次的伴随波场；
[0133]
s2：根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度；
[0134]
s3：根据所述当前轮次的l2目标函数梯度，利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0135]
在一些实施例中，伴随波场和模拟波场相对应，模拟波场和伴随波场进行互相关计算，可以得到梯度。
[0136]
在一些实施例中，对伴随波场的计算过程以及梯度的计算过程进行说明：
[0137]
公式2称为各向同性介质的时间域三维速度-应力一阶声波方程或者波动方程，公式2还可以写为如下形式：
[0138]
e(x,t,m,f)u＝0(8)
[0139]
其中，e(x,t,m,f)称为波动方程算子，u＝(p,v
x
,vy,vz)
t
称为模拟波场向量；m表示速度模型，具体表示地震波在目标区域传播的纵波波速；f表示震源向量，f＝(0,0,0,f
p
)。
[0140]
基于拉格朗日伴随方法，将公式8作为约束项加入到l2目标函数(即公式7)的约束中，将添加了约束项的l2目标函数记为s2函数，s2函数的形式如公式9所示：
[0141][0142]
l(x,t)＝(l
x
,ly,lz,l
p
)
t
(10)
[0143]
其中，l(x,t)＝(l
x
,ly,lz,l
p
)
t
表示伴随波场向量，l
x
表示x方向的拉格朗日乘子，ly表示y方向的拉格朗日乘子，lz表示z方向的拉格朗日乘子，l
p
表示伴随波场。
[0144]
将公式9展开可以得到如下形式：
[0145][0146]
根据最优化理论，求取公式9的极小值需要满足以下三个最优化条件：
[0147][0148][0149][0150]
令即为e(x,t,m,f)u＝0，上式是波动方程，因此是成立的。
[0151]
令并结合边界条件和终止条件，可以推导出伴随波场的计算公式：
[0152][0153]
其中，f
lp
表示伴随源，伴随源的具体含义是：u
obs
(xr,t；xs)和u
syn
(xr,t；xs)的残差；伴随源可以表示“模拟记录和观测记录的差异”。
[0154]
利用公式15计算伴随波场需要逆时加载模拟波场然后反向递推。通过公式15，可以计算出伴随波场l
p
。
[0155]
最后，令并且假定密度为已知的定值，可以得到s2函数对纵波波速的梯度的计算公式：
[0156][0157]
其中，表示s2函数对纵波波速的梯度，表示s2函数对体积模量的梯度，表示体积模量对纵波波速的梯度。
[0158]
根据公式16，可以得到：s2函数对纵波波速的梯度是伴随波场与模拟波场质点振动速度的散度零延迟乘积沿时间的积分。由于利用公式15计算伴随波场l
p
，需要逆时加载模拟波场然后反向递推，并且计算伴随源需要预先获得所有时刻的模拟记录和观测记录，因此想要通过公式15的计算结果进一步利用公式16计算梯度，就需要在模拟波场计算的过
程中存储所有时刻的模拟波场。
[0159]
结合公式3、公式4、公式5、公式6、公式16，可以得到基于目标时刻，计算s2函数梯度的公式：
[0160][0161]
其中，t0(x,y,z)表示地震波到达时刻，t
ny
(x,y,z)表示采样时刻，表示s2函数对体积模量的梯度，l
p
表示伴随波场，v
x
表示x轴方向的网格点振动速度分量，vy表示表示y轴方向的网格点振动速度分量，vz表示z轴方向的网格点振动速度分量，t表示时间，表示s2函数对纵波波速的梯度，表示体积模量对纵波波速的梯度，ρ表示密度，v表示纵波波速，k表示体积模量，x、y、z分别为网格点的x轴坐标值、y轴坐标值、z轴坐标值，表示网格点振动速度分量的散度。
[0162]
在一些实施例中，根据拉格朗日乘子方法，s2函数的对纵波波速的梯度就是l2目标函数在波动方程约束下对纵波波速的梯度，将l2目标函数在波动方程约束下对纵波波速的梯度记为gn，则可以得到因此，根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度，包括：
[0163]
按照如下算式计算当前轮次的l2目标函数梯度：
[0164][0165]
其中，gn表示l2目标函数对纵波波速的梯度，简称l2目标函数梯度，n表示当前轮次。
[0166]
在一些实施例中，使用公式18计算梯度也被称为基于目标时刻的互相关计算。使用公式16计算梯度就是现有技术的全时刻互相关计算。通过公式16可知，现有技术的全时刻互相关计算需要用到从0至t
max
时间段内的所有时刻的波场，因此计算量大。通过公式18可知，基于目标时刻的互相关计算，只需要用到目标时刻的波场以及t0时刻的波场，因此可以有效减少计算量，提升数据处理效率与速度。
[0167]
在一些实施例中，根据所述模拟记录和所述观测记录的差异，计算当前轮次的伴随波场，具体包括：根据公式15计算当前轮次的伴随波场；
[0168]
在一些实施例中，根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度，具体包括：根据公式18计算当前轮次的l2目标函数梯度。
[0169]
在一些实施例中，根据当前轮次的l2目标函数梯度，利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型，具体包括：
[0170]
s1：根据当前轮次的l2目标函数梯度，确定当前轮次的更新方向；
[0171]
s2：根据当前轮次的更新方向，对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0172]
在一些实施例中，根据当前轮次的l2目标函数梯度，确定当前轮次的更新方向，具体包括：
[0173]
按照如下算式，计算第一轮次的更新方向：
[0174]
d1＝-g1(19)
[0175]
其中，d1表示第一轮次的更新方向。
[0176]
在一些实施例中，根据当前轮次的更新方向，对所述初始速度模型进行更新，得到当前轮次的目标速度模型，具体包括：
[0177]
按照如下算式，对初始速度模型进行更新，得到第一轮次的目标速度模型：
[0178]
m1＝m0+αd1(20)
[0179]
其中，α表示更新步长，m1表示第一轮次的目标速度模型。
[0180]
在一些实施例中，在得到当前轮次的目标速度模型之后，所述方法还包括：
[0181]
s1：根据所述地震波参数和所述当前轮次的目标速度模型，确定下一轮次的全局模拟波场；并从所述下一轮次的全局模拟波场中，确定出下一轮次的目标时刻的模拟波场、下一轮次的模拟记录；
[0182]
s2：根据所述观测记录和所述下一轮次的模拟记录，计算下一轮次的l2目标函数值；
[0183]
s3：检测所述下一轮次的l2目标函数值是否大于预设值；
[0184]
s4：在确定所述下一轮次的l2目标函数值大于预设值的情况下，根据所述下一轮次的目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述当前轮次的目标速度模型进行更新，得到下一轮次的目标速度模型。
[0185]
在一些实施例中，根据所述地震波参数和所述当前轮次的目标速度模型，确定下一轮次的全局模拟波场，具体包括：根据地震波参数和当前轮次的目标速度模型，通过公式2对目标区域进行下一轮次的模拟波场计算，确定下一轮次的全局模拟波场。
[0186]
在一些实施例中，从所述下一轮次的全局模拟波场中，通过公式1、公式3、公式4、公式5、公式6确定出下一轮次的目标时刻的模拟波场，其具体过程可参阅上述实施例，在此不做赘述。
[0187]
在一些实施例中，根据所述观测记录和所述下一轮次的模拟记录，通过公式7计算下一轮次的l2目标函数值。
[0188]
在一些实施例中，在确定所述下一轮次的l2目标函数值大于预设值的情况下，根据所述下一轮次的目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述当前轮次的目标速度模型进行更新，得到下一轮次的目标速度模型，具体包括：
[0189]
s1：根据下一轮次的模拟记录，通过公式15计算下一轮次的伴随波场；
[0190]
s2：根据下一轮次的伴随波场以及下一轮次的目标时刻的模拟波场，通过公式18计算下一轮次的l2目标函数梯度；
[0191]
s3：根据所述下一轮次的l2目标函数梯度，利用梯度类优化算法对当前轮次的目标速度模型进行更新，得到下一轮次的目标速度模型。
[0192]
在一些实施例中，在确定下一轮次的l2目标函数值小于或等于预设值的情况下，说明当前轮次的目标速度模型的精度符合要求，停止对速度模型的更新，利用当前轮次的目标速度模型进行地震偏移成像。
[0193]
在一些实施例中，根据所述下一轮次的l2目标函数梯度，利用梯度类优化算法对当前轮次的目标速度模型进行更新，得到下一轮次的目标速度模型，具体包括：
[0194]
s1：根据下一轮次的l2目标函数梯度，确定下一轮次的更新方向；
[0195]
s2：根据下一轮次的更新方向，对当前轮次的目标速度模型进行更新，得到下一轮次的目标速度模型。
[0196]
在一些实施例中，可以通过如下算式确定下一轮次的更新方向：
[0197][0198]
其中，n表示当前轮次，n+1表示下一轮次，d
n+1
表示第n+1轮次的更新方向，gn表示第n轮次的l2目标函数梯度，g
n+1
表示第n+1轮次的l2目标函数梯度。
[0199]
在一些实施例中，可以通过如下算式，对当前轮次的目标速度模型进行更新：
[0200]mn+1
＝mn+αd
n+1
(22)
[0201]
其中，mn表示第n轮次的目标速度模型，m
n+1
表示第n+1轮次的目标速度模型。
[0202]
在一些实施例中，参阅图2所示，可以按照本技术记载的方法，对目标速度模型进行反复更新，并利用l2目标函数值来检测目标速度模型的精度，l2目标函数值越小说明目标速度模型的精度越高，当l2目标函数值大于或等于预设值时重复迭代更新，当l2目标函数值小于预设值时停止更新，最终输出符合要求的目标速度模型。
[0203]
在一些实施例中，图3a表示现有技术中的全时刻互相关计算法，其中震源波场就是模拟波场，t
max
表示模拟波场传播的最大时间，t表示模拟波场正向传播的时间，-t表示伴随波场反向传播的时间。图3a中的全时刻互相关计算法需要已知所有时刻的模拟波场并进行互相关计算，图3a中每个圆点表示一个时刻，乘号表示进行互相关计算。和图3a相比，本技术记载的方法避免了现有的全时刻互相关计算中大量无效时刻点的计算，使得波场重构的内存大小与差分精度和地震记录时间长度无关，大幅降低了对计算机的内存需求；同时减少了一次波场重构，即避免了一次额外的求解波动方程，提高了互相关计算梯度的效率；并使得全波形反演摆脱了现有的逆时波场重构限制，将互相关的过程与模拟波场计算过程独立，让其在计算梯度的互相关时无需再考虑不同模拟方法、介质等额外因素带来的震源波场重构问题，极大地简化了梯度的计算过程。
[0204]
在一些实施例中，图3b表示现有技术中的激发振幅方法，t
examp
表示激发振幅时刻点，t
max
表示模拟波场传播的最大时间，t表示模拟波场传播的时间，t对应的波场为模拟波场，-t对应的波场为伴随波场。本技术记载的方法相较于现有技术中的激发振幅方法，并不需要已知精确的激发振幅时刻点，而只需要已知目标时刻。这就降低了因为伴随波场和模拟波场因为激发振幅时刻不同而导致的误差；同时，因为使用了奈奎斯特采样定理进行进一步稀疏采样，能够再次降低存储量并保留信号的原始波形，从而保证互相关计算梯度结果的准确性，并且相较于激发振幅时刻有更好的抗噪性。
[0205]
在一些实施例中，图3c表示本技术中基于目标时刻的互相关计算方法，其中，t表示模拟波场传播的时间，t对应的波场为模拟波场，-t对应的波场为伴随波场，t
max
表示模拟波场传播的最大时间，存储间隔就是ts，起跳时刻表示表示地震波到达时刻t0，表示t0+ts，表示t0+(n-1)ts，表示t0+n
·
ts，此处的n表示目标时刻的个数，有效时刻窗包括所有的目标时刻以及t0。
[0206]
在一个具体的场景示例中，分别使用现有技术中的互相关计算法、本技术记载的方法、现有技术的激发振幅方法，进行梯度计算，结果如图5a至图8c所示。
[0207]
在一些实施例中，图4a表示基于二维overthrust模型的真实速度模型，真实速度模型用于得到观测记录。图4b表示基于二维overthrust模型的初始速度模型。
[0208]
在一些实施例中，图5a表示使用现有技术中的全时刻互相关方法第一次计算l2目标函数梯度的结果图，图5b表示使用本技术提供的目标时刻互相关计算方法第一次计算l2目标函数梯度的结果图，图5c表示使用激发振幅时刻互相关方法第一次计算l2目标函数梯度的结果图，可以看出，在细节上图5b和图5a基本一致，而图5c在细节上与图5a、图5b有一定差距，如箭头所指之处的计算结果存在差异。这主要是由于激发振幅时刻互相关方法仅用了一个时刻的波场计算梯度，而激发振幅时间的计算并不能做到完全精确，因此当激发振幅时间存在误差时，就会直接影响梯度计算的精度。相比之下，本技术记载的方法不依赖激发振幅时刻的计算精度，只要已知目标时刻就可以保证有较好的梯度计算精度。
[0209]
在一些实施例中，图6a表示ts＝50ms所得到的l2目标函数梯度的结果图，图6b表示ts＝100ms所得到的l2目标函数梯度的结果图，图6c表示ts＝150ms所得到的l2目标函数梯度的结果图，图6d表示ts＝200ms所得到的l2目标函数梯度的结果图，从图6a至图6d可以看出，当ts为50ms和100ms时，l2目标函数梯度计算结果基本一致。这说明当地震资料信噪比较高时，缩小ts并不会提高l2目标函数梯度计算精度，反而会增加目标时刻的个数导致内存占用量增大。此外，从图6a至图6d可以看出，当ts过大以至于不符合奈奎斯特采样定理时，就会影响l2目标函数梯度计算结果。因此，只有根据地震资料的信噪比和反演频带合理地设置ts大小，才能较好地平衡梯度计算结果的精度和互相关计算的内存占用量。
[0210]
在一些实施例中，图7a是利用本技术提供的起跳阈值方法得到的地震波到达时间场；图7b是基于激发振幅方法得到的地震波到达时间场结果图。从图7a看出，本技术所使用的地震波到达时间的判断方式避免了传统激发振幅时间场(图7b)中箭头指示处的突变，这种突变会造成梯度计算结果错误，降低反演精度，而本技术记载的方法则很好地避免了这一突变现象，有效改善了目标时刻的计算精度。
[0211]
在一些实施例中，图8a表示使用现有技术中的全时刻互相关方法以及无噪声数据进行全波形反演的结果图，图8b表示使用本技术提供的目标时刻互相关计算方法以及无噪声数据进行全波形反演的结果图，图8c表示使用激发振幅时刻互相关方法以及无噪声数据进行全波形反演的结果图，图8c箭头指示处的结果和图8a存在差异。
[0212]
在一个具体的场景示例中，使用含噪声数据验证本技术方法的抗噪性，并且使用的初始速度模型与真实模型和图4a、图4b一致，结果如图9a至图9c所示。图9a表示使用现有技术中的全时刻互相关方法以及含噪声数据进行全波形反演的结果图，图9b表示使用本技术提供的目标时刻互相关计算方法以及含噪声数据进行全波形反演的结果图，图9c表示使
用激发振幅时刻互相关方法以及含噪声数据进行全波形反演的结果图。从图9c可以看出，由于受到随机噪声的干扰，激发振幅时刻方法的反演结果信噪比较低，推覆体模型浅部的异常体和深部的逆冲构造细节均未反演出来。而图9b受到噪声的干扰明显较小，与图9a基本一致，无论是深部的构造细节还是浅部的异常体均能得到与图9a一致的结果。增加目标时刻的个数可以进一步提高本技术方法的抗噪性，但需要在抗噪性和内存占用量之间做出平衡。因此，相较于激发振幅时刻方法，本技术方法在抗噪性方面表现更好，即使使用信噪比较差的地震数据，也能获得良好的反演结果。
[0213]
在一个具体的场景示例中，图10a表示基于三维overthrust模型得到的真实速度模型，真实速度模型用于得到观测记录。图10b表示基于三维overthrust模型得到的初始速度模型。
[0214]
在一个具体的场景示例中，图11a表示使用图10b初始速度模型采用现有技术中的全时刻互相关方法进行全波形反演的结果图，其中网格点的数量为100
×
115
×
115，且需要存储三个质点振动速度分量，其内存占用量约为3.23gb。
[0215]
在一个具体的场景示例中，图11b表示使用图10b初始速度模型采用本技术提供的目标时刻互相关方法进行全波形反演的结果图，其中网格点的数量为100
×
115
×
115，需要存储100个目标时刻，内存占用量仅为0.50gb左右，约为图11a的15％。即使图11a牺牲一定的重构精度只存1个波场分量p(p表示压力场数据)，内存占用量仍有1.07gb，本技术方法仍然能够多节省一半的内存。并且，本技术方法无需重构波场，减少了1次三维波动方程的正演计算。图11a和图11b对目标区域的反演结果十分接近，两种方法对三维推覆体的反演结果基本一致。
[0216]
基于上述实施例，本技术提供的一种基于局部存储策略的声波全波形反演方法，采用一阶速度-应力声波动方程，在模拟波场计算中根据设定的振幅阈值判断地震波到达时间，并根据奈奎斯特定理给出采样时间步长的取值范围，采样时间步长再根据地震资料信噪比进行调整。随后，在逆时计算伴随波场过程中提取目标时刻模拟波场进行目标时刻互相关计算梯度。在梯度计算中省略掉一次波动方程的求解，也不需要全部时刻的模拟波场，因此减少了计算量，并大幅降低了这个过程的内存占用量，此外还能够获得与常规全局全波形反演梯度几乎一致的结果，保证了计算精度。
[0217]
基于上述地一种基于局部存储策略的声波全波形反演方法，本说明书还提出一种基于局部存储策略的声波全波形反演装置的实施例，参阅图12所示，所述基于局部存储策略的声波全波形反演装置具体包括以下模块：
[0218]
获取模块1201，用于获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；
[0219]
第一计算模块1202，用于根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；
[0220]
第二计算模块1203，用于根据所述观测记录和所述模拟记录，计算初始l2目标函数值；
[0221]
检测模块1204，用于检测所述初始l2目标函数值是否大于预设值；
[0222]
更新模块1205，用于在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行
更新，得到当前轮次的目标速度模型。
[0223]
在一些实施例中，上述第一计算模块1202，具体用于根据所述地震波参数和所述初始速度模型，对所述目标区域进行模拟波场计算，确定全局模拟波场以及地震波到达时刻；根据所述地震波参数，确定采样时间步长、采样点数量；根据所述采样时间步长、所述采样点数量、所述地震波到达时刻，确定目标时刻；根据所述目标时刻，从所述全局模拟波场中，提取出所述目标时刻的模拟波场。
[0224]
在一些实施例中，上述更新模块1205，具体用于根据所述模拟记录所述观测记录的差异，计算当前轮次的伴随波场；根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度；根据所述当前轮次的l2目标函数梯度，利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0225]
需要说明的是，上述实施例阐明的单元、装置或模块等，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本说明书时可以把各模块的功能在同一个或多个软件和/或硬件中实现，也可以将实现同一功能的模块由多个子模块或子单元的组合实现等。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。
[0226]
本说明书实施例还提供一种基于局部存储策略的声波全波形反演方法的计算机存储介质，所述计算机存储介质存储有计算机程序指令，在所述计算机程序指令被执行时实现：获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；根据所述观测记录和所述模拟记录，计算初始l2目标函数值；检测所述初始l2目标函数值是否大于预设值；在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0227]
在本实施例中，上述存储介质包括但不限于随机存取存储器(random access memory,ram)、只读存储器(read-only memory,rom)、缓存(cache)、硬盘(hard disk drive,hdd)或者存储卡(memory card)。所述存储器可以用于存储计算机程序指令。网络通信单元可以是依照通信协议规定的标准设置的，用于进行网络连接通信的接口。
[0228]
在本实施例中，该计算机存储介质存储的程序指令具体实现的功能和效果，可以与其它实施方式对照解释，在此不再赘述。
[0229]
本说明书还提供一种服务器，包括处理器以及用于存储处理器可执行指令的存储器，所述处理器具体实施时可以根据指令执行以下步骤：获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；根据所述观测记录和所述模拟记录，计算初始l2目标函数值；检测所述初始l2目标函数值是否大于预设值；在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟
波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0230]
为了能够更加准确地完成上述指令，参阅图13所示，本说明书实施例还提供了另一种具体的服务器，其中，所述服务器包括网络通信端口1301、处理器1302以及存储器1303，上述结构通过内部线缆相连，以便各个结构可以进行具体的数据交互。
[0231]
其中，所述网络通信端口1301，具体可以用于获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数。
[0232]
所述处理器1302，具体可以用于根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；根据所述观测记录和所述模拟记录，计算初始l2目标函数值；检测所述初始l2目标函数值是否大于预设值；在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。
[0233]
所述存储器1303，具体可以用于存储相应的指令程序。
[0234]
在本实施例中，所述网络通信端口1301可以是与不同的通信协议进行绑定，从而可以发送或接收不同数据的虚拟端口。例如，所述网络通信端口可以是负责进行web数据通信的端口，也可以是负责进行ftp数据通信的端口，还可以是负责进行邮件数据通信的端口。此外，所述网络通信端口还可以是实体的通信接口或者通信芯片。例如，其可以为无线移动网络通信芯片，如gsm、cdma等；其还可以为wifi芯片；其还可以为蓝牙芯片。
[0235]
在本实施例中，所述处理器1302可以按任何适当的方式实现。例如，处理器可以采取例如微处理器或处理器以及存储可由该(微)处理器执行的计算机可读程序代码(例如软件或固件)的计算机可读介质、逻辑门、开关、专用集成电路(application specific integrated circuit，asic)、可编程逻辑控制器和嵌入微控制器的形式等等。本说明书并不作限定。
[0236]
在本实施例中，所述存储器1303可以包括多个层次，在数字系统中，只要能保存二进制数据的都可以是存储器；在集成电路中，一个没有实物形式的具有存储功能的电路也叫存储器，如ram、fifo等；在系统中，具有实物形式的存储设备也叫存储器，如内存条、tf卡等。
[0237]
虽然本说明书提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的手段可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的装置或客户端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境，甚至为分布式数据处理环境)。术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、产品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、产品或者设备所固有的要素。在没有更多限制的情况下，并不排除在包括所述要素的过程、方法、产品或者设备中还存在另外的相同或等同要素。第一，第二等词语用来表示名称，而并不表示任何特定的顺序。
[0238]
本领域技术人员也知道，除了以纯计算机可读程序代码方式实现控制器以外，完
全可以通过将方法步骤进行逻辑编程来使得控制器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种控制器可以被认为是一种硬件部件，而对其内部包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。
[0239]
本说明书可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构、类等等。也可以在分布式计算环境中实践本说明书，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
[0240]
通过以上的实施例的描述可知，本领域的技术人员可以清楚地了解到本说明书可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本说明书的技术方案本质上可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如rom/ram、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，移动终端，服务器，或者网络设备等)执行本说明书各个实施例或者实施例的某些部分所述的方法。
[0241]
本说明书中的各个实施例采用递进的方式描述，各个实施例之间相同或相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。本说明书可用于众多通用或专用的计算机系统环境或配置中。例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备、多处理器系统、基于微处理器的系统、置顶盒、可编程的电子设备、网络pc、小型计算机、大型计算机、包括以上任何系统或设备的分布式计算环境等等。
[0242]
虽然通过实施例描绘了本说明书，本领域普通技术人员知道，本说明书有许多变形和变化而不脱离本说明书的精神，希望所附的权利要求包括这些变形和变化而不脱离本说明书的精神。

技术特征：
1.一种基于局部存储策略的声波全波形反演方法，其特征在于，包括：获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；根据所述观测记录和所述模拟记录，计算初始l2目标函数值；检测所述初始l2目标函数值是否大于预设值；在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。2.根据权利要求1所述的方法，其特征在于，所述方法还包括：在确定所述初始l2目标函数值小于或等于预设值的情况下，利用所述初始速度模型进行地震偏移成像。3.根据权利要求1所述的方法，其特征在于，根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出目标时刻的模拟波场，包括：根据所述地震波参数和所述初始速度模型，对所述目标区域进行模拟波场计算，确定全局模拟波场以及地震波到达时刻；根据所述地震波参数，确定采样时间步长、采样点数量；根据所述采样时间步长、所述采样点数量、所述地震波到达时刻，确定目标时刻；根据所述目标时刻，从所述全局模拟波场中，提取出所述目标时刻的模拟波场。4.根据权利要求3所述的方法，其特征在于，根据所述采样时间步长、所述采样点数量、所述地震波到达时刻，确定目标时刻，包括：按照如下算式确定目标时刻：其中，t0(x,y,z)表示网格点x(x,y,z)的地震波到达时刻，表示网格点x(x,y,z)的采样时刻，表示网格点x(x,y,z)的目标时刻，a(x,y,z)表示网格点x(x,y,z)的采样次数，t
s
表示采样时间步长，a
n
(x,y,z)表示网格点x(x,y,z)的采样点数量，x表示网格点的x轴坐标值，y表示网格点的y轴坐标值，z表示网格点的z轴坐标值。5.根据权利要求1所述的方法，其特征在于，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型，包括：根据所述模拟记录和所述观测记录的差异，计算当前轮次的伴随波场；根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度；根据所述当前轮次的l2目标函数梯度，利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。6.根据权利要求5所述的方法，其特征在于，在得到当前轮次的目标速度模型之后，所述方法还包括：根据所述地震波参数和所述当前轮次的目标速度模型，确定下一轮次的全局模拟波
场；并从所述下一轮次的全局模拟波场中，确定出下一轮次的目标时刻的模拟波场、下一轮次的模拟记录；根据所述观测记录和所述下一轮次的模拟记录，计算下一轮次的l2目标函数值；检测所述下一轮次的l2目标函数值是否大于预设值；在确定所述下一轮次的l2目标函数值大于预设值的情况下，根据所述下一轮次的目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述当前轮次的目标速度模型进行更新，得到下一轮次的目标速度模型。7.根据权利要求5所述的方法，其特征在于，根据当前轮次的伴随波场以及所述目标时刻的模拟波场，计算当前轮次的l2目标函数梯度，包括：按照如下算式计算当前轮次的l2目标函数梯度：其中，t0(x,y,z)表示地震波到达时刻，t
ny
(x,y,z)表示采样时刻，表示s2函数对体积模量的梯度，l
p
表示伴随波场，v
x
表示x轴方向的网格点振动速度分量，v
y
表示表示y轴方向的网格点振动速度分量，v
z
表示z轴方向的网格点振动速度分量，t表示时间，表示s2函数对纵波波速的梯度，表示体积模量对纵波波速的梯度，ρ表示密度，v表示纵波波速，k表示体积模量，x、y、z分别为网格点的x轴坐标值、y轴坐标值、z轴坐标值，表示网格点振动速度分量的散度，g
n
表示l2目标函数对纵波波速的梯度，简称l2目标函数梯度，n表示当前轮次。8.根据权利要求5所述的方法，其特征在于，根据当前轮次的l2目标函数梯度，利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型，包括：根据当前轮次的l2目标函数梯度，确定当前轮次的更新方向；根据当前轮次的更新方向，对所述初始速度模型进行更新，得到当前轮次的目标速度模型。9.一种基于局部存储策略的声波全波形反演装置，其特征在于，包括：获取模块，用于获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；第一计算模块，用于根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从所述全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；第二计算模块，用于根据所述观测记录和所述模拟记录，计算初始l2目标函数值；检测模块，用于检测所述初始l2目标函数值是否大于预设值；更新模块，用于在确定所述初始l2目标函数值大于预设值的情况下，根据所述目标时刻
的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。10.一种计算机可读存储介质，其特征在于，其上存储有计算机指令，所述指令被处理器执行时实现权利要求1至8中任一项所述方法的步骤。

技术总结
本说明书提出一种基于局部存储策略的声波全波形反演方法及装置。该方法包括：获取目标区域；并确定所述目标区域的初始速度模型、观测记录、地震波参数；根据所述地震波参数和所述初始速度模型，确定全局模拟波场；并从全局模拟波场中，确定出模拟记录、目标时刻的模拟波场；根据观测记录和模拟记录，计算初始L2目标函数值；检测所述初始L2目标函数值是否大于预设值；在确定所述初始L2目标函数值大于预设值的情况下，根据所述目标时刻的模拟波场进行互相关计算并利用梯度类优化算法对所述初始速度模型进行更新，得到当前轮次的目标速度模型。基于上述方法能够减少全波形反演过程中的数据计算量与数据存储量，提升全波形反演的反演效率。反演效率。反演效率。


技术研发人员：杨瀚 陈汉明 王玲谦 周辉
受保护的技术使用者：中国石油大学（北京）
技术研发日：2023.07.07
技术公布日：2023/10/15