CN113050160A - 一种数据阻尼全波形反演方法、装置和计算机设备 - Google Patents

Abstract

本文提供了一种数据阻尼全波形反演方法、装置和计算机设备，其中方法包括：基于观测数据、模拟数据和衰减函数，构建表征观测数据与模拟数据之间数据残差的目标函数；获取震源的正传波场；根据衰减函数获得虚拟震源的反传波场；根据正传波场、反传波场和目标函数关于速度场模型参数的梯度计算表达式，计算速度场模型参数的梯度；根据梯度对速度场模型参数进行迭代更新，直至数据残差达到设定阈值或迭代更新次数达到预设值。本文提供的一种数据阻尼全波形反演方法，使用分段衰减函数作用于反传源以模拟得到反传波场，再分别对正传和反传波场实施波场解耦计算，能够有效地降低反演中高波数噪声的影响，从而得到高精度的反演结果。

Description

一种数据阻尼全波形反演方法、装置和计算机设备

技术领域

本发明涉及石油地地震勘探速度建模技术领域，尤其是一种数据阻尼全波形反演方法、装置和计算机设备。

背景技术

石油和天然气等能源产业支撑着我国经济的快速稳定发展。随着石油、天然气勘探开发程度的加深，一些埋藏浅、构造简单的油气藏已经逐渐开发殆尽。为了实现对埋藏深、构造复杂及岩性油气藏的开发利用，亟需先进的勘探手段和高精度的成像技术。近些年来，主要发展宽频带、宽方位、高密度的地震数据采集技术，通过地震成像技术来实现地下参数的高信噪比、高分辨率以及高保真度的建模。最终依据成像结果，进行地质解释、油藏定位以及井位确定等工作。

全波形反演(Full Waveform Inversion，FWI)以层析反演或者速度分析结果为初始模型，可实现高分辨率的速度建模，结合保幅的偏移成像方法，获得高质量的成像剖面，为地质解释提供支撑和依据，因此全波形反演成为近些年来的研究热点之一。然而，计算效率低下是制约全波形反演进行大规模计算以及实际资料处理的主要因素，现有的一种解决方法是使用震源编码策略，其思想是通过某种震源组合方式将很多单炮形成一个超级震源，从而减少波场模拟的次数。其中被广泛应用的一种为基于频率选择的同时源全波形反演方法，它可以获得无串扰的反演结果，但是在一定程度上存在着反演质量差的现象：当初始模型和真实模型相差较大时，低频模型不能被充分的更新，会导致模型里的高波数信息不能正确地更新，从而产生比较严重的高波数噪声。

因此，如何提供一种新的同时源反演方法，能够克服现有技术存在的反演精度差的问题是本领域技术人员亟需解决的技术难题。

发明内容

针对现有技术的上述问题，本文的目的在于，提供一种数据阻尼全波形反演方法、装置和计算机设备，以解决现有技术中同时源反演产生的高波数噪声问题。

为了解决上述技术问题，本文的具体技术方案如下：

第一方面，本文提供一种数据阻尼全波形反演方法，包括：

基于观测数据、模拟数据和衰减函数，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；

获取震源的正传波场；

根据所述衰减函数获得虚拟震源的反传波场，所述虚拟震源与所述震源相对应；

根据所述正传波场、所述反传波场和所述目标函数关于速度场模型参数的梯度计算表达式，计算所述速度场模型参数的梯度；

根据所述梯度对所述速度场模型参数进行迭代更新，直至所述数据残差达到设定阈值。

具体地，所述衰减函数为：

其中，α_γ为常数，m∈Z，T_p为最小积分区间。

具体地，所述目标函数为：

其中，h(α_γ,t)＝h^′2(α_γ,t)；d^obs(X_r,Ω,t)为在检波点X_r处获得的所述观测数据，u(X_r,Ω,t；m)为当前速度场模型参数m下获得的所述模拟数据；Ω表示角频率，x表示计算区域，并且

m(x)为计算区域xx处的速度场模型参数，X_r为检波点的位置信息。

进一步地，所述震源通过多个位于不同位置的单频子波同时激发获得，所述震源为：

其中，

为位于X_s的单频随机信号；

是表示编码信号中第i个信号的位置，

ω_i表示编码信号中第i个信号的角频率，Ω＝(ω₁,ω₂,…,ω_n)，n_s表示一个震源内编码信号的个数。

具体地，所述获取震源的正传波场包括：

求解正传波场的模拟方程，得到所述正传波场，所述模拟方程为：

其中，

x表示计算区域，X_s震源位置信息，L[·]是声波正演模拟算子。

具体地，所述根据所述衰减函数获得虚拟震源的反传波场，包括：

根据所述衰减函数，计算虚拟震源，所述虚拟震源为：

其中，反传源为：

根据所述虚拟震源，获取伴随方程，所述伴随方程为：

其中，

表示伴随算子，u′(x,Ω,t；m)表示反传波场；

求解所述伴随方程，获得所述虚拟震源的反传波场。

进一步地，所述根据所述正传波场、所述反传波场和所述目标函数关于速度场模型参数的梯度计算表达式，计算速度场模型参数的梯度，包括：

对所述正传波场和所述反传波场分别进行无串扰解耦，得到波场变量；

根据所述波场变量和所述梯度计算表达式，计算所述速度场模型参数的梯度。

进一步地，所述对所述正传波场和所述反传波场分别进行无串扰解耦，包括：

获取第一参考信号和第二参考信号，所述第一参考信号的振幅和频率均与所述第二参考信号的振幅和频率相同，所述第一参考信号的相位与所述第二参考信号的相位相差90°；

根据所述第一参考信号和所述第二参考信号，对所述正传波场和所述反传波场分别进行无串扰解耦。

第二方面，本文还提供一种数据阻尼全波形反演装置，包括：

目标函数构建模块，用于基于观测数据、模拟数据和衰减函数，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；

正传波场获取模块，用于获取震源的正传波场；

反传波场获取模块，用于根据所述衰减函数获得虚拟震源的反传波场，所述虚拟震源与所述震源相对应；

梯度计算模块，用于根据所述正传波场、所述反传波场和所述目标函数关于速度场模型参数的梯度计算表达式，计算速度场模型参数的梯度；

迭代更新模块，用于根据所述梯度对所述速度场模型参数进行迭代更新，直至所述数据残差达到设定阈值或者迭代更新次数达到预设值。

第三方面，本文还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上述技术方案提供的数据阻尼全波形反演方法。

采用上述技术方案，本文提供的一种数据阻尼全波形反演方法，使用分段衰减函数作用于反传源以模拟得到反传波场，再分别正传和反传波场实施波场解耦计算，能够有效地降低反演中高波数噪声的影响，从而得到高精度的反演结果。

为让本文的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附图式，作详细说明如下。

附图说明

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

图1示出了本文实施例提供的一种数据阻尼全波形反演方法的步骤示意图；

图2示出了获得反传波场的步骤示意图；

图3a至图3d示出了采用相敏检测法对衰减信号进行解耦恢复的示意图；

图4a至图4b示出了Overthrust模型下真实速度场模型参数和初始速度场模型参数的示意图；

图5a至图5c示出了使用分段式衰减函数作用于反传源的前后对比图；

图6a至图6d示出了采用相敏检测方法对未衰减的反传源波场进行解耦得到的其对应的解耦波场和对经分段衰减函数衰减处理的反传源波场进行解耦得到的其对应的解耦波场的对比示意图；

图7a至图7b示出了分别使用未经衰减的反传源信号和经衰减的反传源信号计算得到的梯度的对比示意图；

图8a至图8d示出了本发明实施例提供的数据阻尼全波形反演方法运用到Overthrust模型得到的反演结果与常规方法得到的反演结果的对比示意图；

图9a至图9e示出了本发明实施例提供的数据阻尼全波形反演方法运用到Marmousi模型中得到的反演结果与常规方法得到的反演结果的对比示意图；

图10a至图10g示出了本发明实施例提供的数据阻尼全波形反演方法运用到某工区实际场景下得到的反演结果与常规方法运用于此场景下得到的反演结果的对比示意图；

图11示出了本发明实施例提供的数据阻尼全波形反演装置；

图12示出了一种计算机设备的结构示意图。

附图符号说明：

10、目标函数构建模块；

20、正传波场获取模块；

30、反传波场获取模块；

40、梯度计算模块；

50、迭代更新模块；

1202、计算机设备；

1204、处理器；

1206、存储器；

1208、驱动机构；

1210、输入/输出模块；

1212、输入设备；

1214、输出设备；

1216、呈现设备；

1218、图形用户接口；

1220、网络接口；

1222、通信链路；

1224、通信总线。

具体实施方式

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

需要说明的是，本文的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本文的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、装置、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

石油和天然气等能源产业支撑着我国经济的快速稳定发展。随着石油、天然气勘探开发程度的加深，一些埋藏浅、构造简单的油气藏已经逐渐开发殆尽。为了实现对埋藏深、复杂构造及岩性油气藏的开发利用，亟需先进的勘探手段和高精度的成像技术。近些年来，主要发展宽频带、宽方位、高密度的地震数据采集技术，通过地震成像技术来实现地下参数的高信噪比、高分辨率以及高保真度的建模。最终依据成像结果，进行地质解释、油藏定位以及井位确定等工作。

就成像方法而言，大致可分为：叠加成像技术、偏移成像技术和波形反演技术。其中，波形反演成像技术具有较高的成像精度，主要以最小二乘成像和波形反演为代表，能够实现构造信息的保幅成像和岩性油气藏的成像。对偏移成像来说，建立初始速度场是非常关键的技术环节。只有建立准确的背景速度场才能保证正确的成像结果。而全波形反演(Full Waveform Inversion，FWI)是以层析反演或者速度分析结果为初始模型，利用全波形反演技术实现高分辨率的速度建模，结合保幅的偏移成像方法，能够获得高质量的成像剖面，为地质解释提供支撑和依据，因此成为近些年来的研究热点之一。但计算效率低下是制约全波形反演进行大规模计算以及实际资料处理的主要因素，因此提升全波形反演的效率就显得十分重要。除了使用高性能计算平台加速之外，通常使用算法类的效率优化方案来加速。例如，使用频率域的全波形反演，通过有限的频率反演来节省计算耗时；使用一些收敛性更好的迭代算法，如L-BFGS算法，也可以减少迭代次数，提升计算效率；再如在三维大规模计算中，由于频繁地从硬盘里读取大量的波场数据的操作非常耗时，所以引进有效边界存储技术，可以用额外的正演模拟计算来换取数据存储，从而起到均衡内存和计算效率的作用。另外一类方法是使用震源编码策略，震源编码是指通过某种震源组合方式将很多单个震源形成一个超级震源，从而减少波场模拟的次数,可以同时实现几十个甚至上百个震源的数值模拟，从而极大地减少全波形反演的计算量。

相比于传统的多震源模拟(震源编码)方法，比如随机极性编码、随机时移、平面波编码以及随机优化算法，基于频率选择的多震源编码方法由于其不仅能够保证计算效率，而且可以获得不受编码串扰噪声影响的反演结果而得到广泛应用。在此基础上，又发展得到了时间–频率域(混合域)的多震源反演方法。这种方法是指，在混合域多震源模拟方法中，将一系列不同频率的单频源组合起来形成了一个编码源，并在时间域完成波场的模拟。超级震源的梯度计算则是在频率域完成的，通过计算每一个单频子波对应的梯度然后叠加形成编码信号对应的梯度。这种方法的关键点是如何准确地解耦混叠波场。尽管基于频率选择的同时源全波形反演方法可以获得无串扰的反演结果，但是基于这种方法得到的反演结果，在一定程度上存在着反演质量差的现象：当初始模型和真实模型相差比较大时，低频模型不能被充分地更新，会导致模型里的高波数信息不能正确地更新，会产生比较严重的高波数噪声。为了解决同时源全波形反演方法存在的质量差的问题，本文实施例提供了一种数据阻尼全波形反演方法。图1是本文实施例提供的一种数据阻尼全波形反演方法的步骤示意图，本说明书提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的劳动可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的系统或装置产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行。

为帮助本领取技术人员更好地理解本文的技术方案，本说明书对一种无串扰的同时源全波形反演方法进行如下解释说明，该方法包括：

S1：建立表征观测数据与模拟数据之间数据残差的目标函数；

所述目标函数采用最小化位于检波点X_r的观测数据d^obs(X_r,Ω,t)和模拟数据u(X_r,Ω,t；m)的残差，具体为：

其中，x表示计算区域，m(x)表示计算区域x处的速度场模型参数；

X_r为检波点的位置信息；Ω表示角频率，t表示时间；m表示速度场模型参数，u(X_r,Ω,t；m)表示在当前速度场模型参数m下获得的模拟数据。

S2：通过模拟，获取震源的正传波场和反传源的反传波场；

所述震源为通过多个位于不同位置的单频子波同时激发获得，即所述震源为同时源模拟波场。

给定一系列位于X_s的单频随机信号s(X_s,Ω)，则所述震源可表达为：

其中，X_s为震源位置信息，

n_s表示一个震源内编码信号的个数；

是表示震源编码信号中第i个信号的位置，

ω_i表示编码信号中第i个信号的角频率，Ω＝(ω₁,ω₂,…,ω_n)。

则震源

对应的正传波场模拟可以表示为：

其中，

u(x,Ω,t；m)表示正传波场，L[·]为声波正演模拟算子。

相应地，伴随方程为：

式(4)中，

表示伴随算子，u′(x,Ω,t；m)表示反传波场；

为所述反传源，则所述反传源可表达为：

式(5)中，

表示编码信号中第i个信号对应检波点的位置。求解式(4)和式(5)，即可模拟得到震源的正传波场和虚拟震源的反传波场。

通常情况下，由于观测数据和模拟数据的观测系统不一致，所以需要分别从观测和模拟数据对应的单个信号数据中提取相同频率的信号，在观测数据的观测系统约束下求取数据的残差，以此解决观测系统不一致的问题。从而可以适用于海上拖缆或者陆上滚动采集观测系统，消除了观测系统不一致对反演的影响。

以及，S3使用正传波场和反传波场进行梯度计算，迭代计算目标函数数据残差以更新速度场模型参数。

具体地，包括如下步骤：

S31：获得目标函数对速度场模型参数的梯度计算表达式；

S32：根据同时源模拟进行正演和反传计算，获得解耦的波场变量；

本说明书中，梯度计算表达式为正传波场与反传波场的特定形式的互相关，可简写为：

梯度简写形式G(x,m)表达式中，ψ_s和ψ_r分别表示正传波场与反传波场。从而，由(6)式的第二个等式可以看到，等式右端的第一项表示每一个编码信号对应波场数据的自相关运算；第二项表示当前编码信号与其它编码信号的波场互相关运算，即串扰。由(6)式可知，只有从混叠的波场中解耦出所有编码信号的波场，然后只计算当前编码信号对应的正反传波场的自相关，才能从根本上消除串扰噪声。

即，有效地解耦混叠波场是解决串扰噪声的关键，本说明书实施例引入相敏检测(PSD)法来解决这个问题。PSD的主要思想是：利用两个与未知信号具有相同频率但彼此相位相差90°的参考信号进行时间积分来提取振幅为E相位为θ的未知信号(即目标信号)。

本说明书实施例中，取E_ref sin(ω_it+θ_ref)为第一参考信号为，取E_ref cos(ω_it+θ_ref)为第二参考信号；其中，E_ref为第一参考信号和第二参考信号的振幅；θ_ref为所述第一参考信号和第二参考信号的相位。由此可以看出，所述第一参考信号的振幅和频率均与所述第二参考信号的振幅和频率相同，所述第一参考信号的相位与所述第二参考信号的相位相差90°

由于空间某一点的混叠波场，均可以表示为：

其中，

分别为每个单频编码信号的振幅，

分别为每个单频编码信号的频率；以及

分别为每个单频编码信号的相位。

使|ω_i-ω_j|≥pΔω(i,j∈n_s,p∈Z)，且Δω＝2π/T。这里的T是地震记录的长度，p是最小的整数，pΔω表示编码信号之间的最小频率间隔。当定义最小的积分区间为T_p＝T/p，对于任何q(p∈[1,p])倍的T_p区间qT_p都可以作为合理的积分区间。通常，需要一定的时间等待波场达到稳态状况，所以q比p要小。

在积分区间qT_p上，将所述第一参考信号与混叠波场的信号相乘并积分，得到：

T_p为最小积分区间。

以及，将所述第二参考信号与所述混叠波场的信号相乘并积分，得到：

从而，与所述第一参考信号和第二参考信号频率相同的编码信号的对应波场便可以解耦出来，其计算公式为：

上述积分过程主要利用了三角函数正交性，只有与第一参考信号和第二参考信号频率相同的信号其得到积分结果不为零，其他与第一参考信号和第二参考信号频率不同的信号其得到的积分结果为零，从而实现对空间中任意一混叠波场中各频率编码信号其对应波场的解耦。

S33：利用目标函数对速度场模型参数的梯度计算表达式和解耦的波场变量，获取模型参数梯度；

其中速度模型的梯度计算式为：

其中，Re表示取向量的实数部分；

为正传波场中，频率为ω_i的编码信号解耦得到的其对应的波场变量，

为反传波场中，频率为ω_i的编码信号经解耦得到其对应的波场变量。式(11)即是对频率均为ω_i的两个波场变量进行正相关运算，从而实现无串扰的梯度求解。

S34：采用共轭梯度方法求取更新方向，对模型参数梯度进行迭代更新，直至残差达到设定值或者迭代更新次数达到预设值，输出速度模型反演参数，确定速度模型。

采用共轭梯度方法求取更新方向Δm_k+1，

式(12)中的下标k表示迭代次数，s_k是用于计算共轭方向的一个中间变量，上标T表示矩阵转置，更新后的速度场模型为：

m_k+1＝m_k+t_kΔm_k+1 (13)

即将式(12)中计算得到的Δm_k+1带入到式(13)中，式(13)中的t_k表示选取的更新步长，t_kΔm_k+1即为此次迭代过程中速度场模型参数的更新量。

在上述无串扰的同时源反演方法中，使用一系列单频子波同时激发形成混叠波场；并采用了相敏检测法从混叠波场中解耦出每一个编码信号对应的波场变量，实现了无串扰的梯度计算，从而有利于获得准确的速度场模型参数的更新方向。

本说明书实施例基于无串扰的同时源模拟及反演的基础上，提供的一种数据阻尼全波形反演方法，具体的如图1所示，所述方法包括如下步骤：

S100：基于观测数据、模拟数据和衰减函数，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；

所述目标函数为：

其中，h′(α_γ,t)是衰减函数，α_γ为常数，定义α_γ ^-1为衰减因子，本说明书实施例中优选地，衰减因子的一般取值范围为[1,10]，其他变量的含义与式(1)中一致。

S200：获取震源的正传波场；

本说明书实施例中，所述震源通过多个位于不同位置的单频子波同时激发获得，即所述震源为同时源模拟波场。

所述震源为：

其中，

为位于X_s的单频随机信号；X_s为震源位置信息，下标s表示震源source，

n_s表示一个震源内编码信号的个数；

是表示编码信号中第i个信号的位置，

ω_i表示编码信号中第i个信号的角频率，Ω＝(ω₁,ω₂,…,ω_n)。

则震源

的对应的正传波场模拟可以表示为：

其中，

x表示计算区域，L[·]为声波正演模拟算子。

较为优选地，本说明书实施例中，采用一阶速度—应力波动方程求解声波正演模拟算子，从而求解上式可模拟得到所述震源的正传波场。即在本说明书实施例提供的一种数据阻尼全波形反演方法中，震源及所述震源的正传波场分别与无串扰的同时源反演方法中的震源和震源的正传波场相同。

S300：根据所述衰减函数获得虚拟震源的反传波场，所述虚拟震源与所述震源相对应；

具体地，如图2所示，步骤S300可以包括以下步骤：

S310：根据所述衰减函数，计算虚拟震源；

所述虚拟震源为：

其中，反传源为：

S320：根据所述虚拟震源，获得伴随方程；所述伴随方程为：

其中，

表示伴随算子。

S330：求解所述伴随方程，获得所述虚拟震源的反传波场。同样地，本说明书实施例中，优选地，采用一阶速度—应力波动方程求解上式以模拟得到反传波场。

S400：根据所述正传波场、所述反传波场和所述目标函数关于速度场模型参数的梯度计算表达式，计算所述速度场模型参数的梯度；

S410：对所述正传波场和所述反传波场分别进行无串扰解耦，得到波场变量；

当衰减函数是指数衰减函数并将其直接运用到多震源数据衰减中时，由于PSD积分的缘故，将导致无法完全分离混叠波场，不能实现无串扰的同时源反演方法。以式(9)为例，首先将h(α_γ,t)代入方程，得到：

式(9)中，要实现波场的解耦，必须保证qT_p周期上与目标信号无关的积分项的积分值为零。如果直接将指数衰减函数代入，则无法保证在最大积分区间qT_p上不同频率的三角函数的积分值为零。由于在最小积分区间T_p上保证除了与频率相关的积分项积分都为零，而在q个不同的T_p之间进行数据衰减，从而有效地衰减波场振幅信息，同时又保证波场相位信息不变。为了实现这一目的，本说明书实施例中，提供了一种新的衰减函数：

其中，α_γ ^-1为衰减因子，其取值范围一般为[1,10]；m∈Z,m∈[1,q]。即本说明书实施例中，所述衰减函数为分段函数，且需要说明的是，分段函数在每个T_p周期上为一个常数，根据这个条件简化公式(17)，可以得到：

其中，

同理可得：

将式(18)和式(19)带入到式(10)中，可以得到衰减后，各频率编码信号对介质响应，即得到每个频率编码信号对应的波场，其振幅和相位与未衰减项之间的关系为E′_i＝E_iA_α,θ′_i＝θ_i。由此实现了无串扰解耦。

如图3b所示，为衰减函数的示意图，其中，光顺曲线为衰减因子α_γ ^-1为1、2和5的指数函数，阶梯状的折线为本说明书实施例提供的衰减因子α_γ ^-1为1、2和5的分段式衰减函数。

S420：根据所述波场变量和所述梯度计算表达式，计算所述速度场模型参数的梯度。

所述目标函数关于速度场模型参数的梯度计算表达式为：

其中，Re表示取向量的实数部分，

为正传波场中，频率为ω_i的编码信号经解耦得到的其对应的波场变量；

为反传波场中，频率为ω_i的编码信号经解耦得到其对应的波场变量。

通过正演和反传计算出相应的波场变量，利用梯度表达式求取梯度，采用共轭梯度方法求取更新方向：

其中，下标k表示迭代次数，上标T表示矩阵转置。

S500：所述梯度对所述速度场模型参数进行迭代更新，直至所述数据残差达到设定阈值或者迭代更新次数达到预设值。

m_k+1＝m_k+t_kΔm_k+1 (13)

t_k表示选取计算更新。

获得更新后的速度场模型参数后，将其带入式(14)中重新计算残差判断是否满足预设的优化条件或者判断当前迭代更新次数是否达到预设值。若不满足，则继续进行波场正传、反传，以及无串扰解耦、梯度计算全过程，循环执行直至速度场模型参数足够优化或者直至满足预设的迭代的次数。

需要说明的是，在第k+1次迭代过程中，对反传源进行衰减的衰减函数，其衰减因子可以与第k次迭代过程中所采用的衰减函数的衰减因子相同，也可以不同。

本说明书实施例在无串扰的同时源模拟及反演方法的基础上，提供了一种数据阻尼全波形反演方法，采用了分段式的衰减函数来对反传源进行衰减，进而根据一阶速度—应力波动方程模拟得到衰减后的反传波场。在使用PSD法对这种波场进行解耦时，不会对目标波场的相位信息破坏，同时对该波场的振幅信息只是乘以一个系数项，即对反演中所使用的数据的振幅信息进行调节。通过本说明书实施例中对反传源的衰减方法，不仅仍可以用PSD方法解耦混叠波场以实现无串扰的同时源反演；并且由于只作用于反传源，因此不会显著增加计算量。

除此之外，本说明书实施例提供的一种数据阻尼全波形反演方法，还充分利用了波形反演的优势，并综合使用了振幅和相位信息。在开始反演时，所用衰减函数的衰减因子使数据的振幅信息衰减程度最大，保证反演尽可能匹配相位信息；随着迭代次数的增加，可以逐渐加入更多的振幅信息，从而使反演算法的鲁棒性增加。同时在反演过程中降低振幅的影响，可以有效地降低高频噪声，提高反演精度。

综上，本发明实施例提供的一种数据阻尼全波形反演方法，通过多个单频编码信号同时激发实现同时源模拟；并使用衰减函数对反传源进行衰减处理，采用一阶速度—应力声波动方程获得震源的正传波场和虚拟震源的反传波场；采用相敏检测方法解耦波场以实现无串扰噪声的梯度计算，不仅能够保证相位信息，还可以衰减部分振幅信息起到降低数据量的作用。

图3a至图3d为本发明实施例提供的采用相敏检测方法对衰减信号进行解耦恢复的示意图。具体地，图3a从上到下分别是频率为6.25Hz的单频信号A、频率为7.5Hz的单频信号B、频率为8.75Hz的单频信号C以及根据这三个单频信号得到的叠加信号D，信号A、信号B和信号C的信号长度均为4s，时间采样步长均为1ms。图3b为衰减函数的示意图，其中，光顺曲线分别为衰减因子为1、2和5的指数函数，阶梯状的折线分别为衰减因子为1、2和5的分段式衰减函数。图3c为通过相敏检测方法解耦后的叠加得到的信号,从上到下分别为：叠加信号D经衰减因子为1的指数型衰减函数阻尼衰减后得到的信号、叠加信号D经衰减因子为1的分段式衰减函数阻尼衰减后得到的信号、对经衰减因子为1的指数型衰减函数衰减后的信号进行解耦并叠加得到的信号D’以及对经衰减因子为1的分段式衰减函数衰减后的信号进行解耦并叠加得到的信号D”。图3d为相敏检测方法解耦后叠加得到的信号与原叠加信号的误差，从上到下分别为：信号D’与原叠加信号D之差、信号D”与原叠加信号D之差、振幅归一化后的信号D’与叠加信号D之差的误差以及振幅归一化后的信号D”与叠加信号D之间的误差。

由图3a至图3d可以看出，指数型衰减函数和分段式衰减函数均可以衰减信号的振幅，但是指数型衰减函数数据阻尼的方式破坏了信号的相位信息，因此无法完全解耦信号。相反地，分段式衰减函数数据阻尼的方式可以保证信号相位信息不被改变，即图3d中最下方的曲线，两个信号之差几乎为0。

图4a至图4b为本发明实施例所使用的Overthrust模型真实速度场模型参数和初始速度场模型参数的示意图。其中图4a为真实的速度模型，图4b为反演所用的初始速度模型。

图5a至图5c为本发明实施例一种使用分段式衰减函数作用于反传源的前后对比图。具体地，图5a表示未衰减的反传源；图5b表示所使用的分段式衰减函数；图5c表示经衰减处理后的反传源，即将图5b中的分段式衰减函数作用于图5a中的反传源所得到的数据。从图5c可以看出，经分段式衰减函数数据阻尼处理后，反传源振幅信息沿着时间方向发生了衰减(图5a中有4个完整的周期信号，而经衰减后图5c中，仅有一个振幅完整的周期性)，而相位信息不变(衰减前后，同相轴沿时间方向的周期不变)。因此，对比图5a和图5c，可知使用本说明书实施例提供的分段式衰减函数数据阻尼的策略，可以有效地减少反传数据量。

图6a至图6d为采用相敏检测方法对未衰减的反传源波场进行解耦得到的其对应的解耦波场和对经分段衰减函数衰减的反传源波场进行解耦得到的其对应的解耦波场的对比示意图。图6a和图6b分别表示未经衰减函数衰减处理的伴随源对应波场的振幅和相位信息；图6c和图6d分别表示经衰减函数衰减处理后的伴随源对应波场的振幅和相位信息。对比图6b与图6d可以看出对反传源使用分段衰减函数阻尼策略，其相位信息不变；而对比图6c与图6a可以看出，分段衰减函数阻尼策略使得振幅信息变得光滑并且能量更加均匀。

图7a至图7b为本发明实施例中分别使用未经衰减的反传源信号和经衰减的反传源信号计算得到的梯度的对比示意图。图7a是使用未经衰减处理的反传源信号进行同时源全波形反演求取的梯度；图7b是使用经分段衰减函数衰减处理的反传源信号进行数据阻尼同时源全波形反演求取的梯度。对比图7a和图7b可以看出，对反传源使用阻尼策略对其进行衰减后，其目标函数对应模型的梯度结果(黑线线框中)连续性更好，深层信息更加突出，从而对速度场模型参数的更新方向更加准确，也就是说本说明书实施例提供的采用分段式衰减函数处理反传源的策略，更有利于得到准确的速度场模型参数。

图8a至图8d为本发明实施例提供的一种数据阻尼全波形反演方法运用到Overthrust模型得到的反演结果与常规方法得到的反演结果的对比示意图。图8a是未经衰减处理得到的同时源反演结果，即同时源全波形反演方法得到的结果；图8b是本说明书实施例提供的数据阻尼(同时源)全波形反演方法得到的结果；图8c是采用常规的序列源反演方法得到的结果，序列源反演方法中震源为多个雷克子波逐个激发形式，该方法能够得到较好的反演结果但其效率低下；图8d由左至右、由上至下分别是准确的速度场模型参数(即图4a)和图8a、图8b以及图8c反演结果的局部放大示意图。对比图8b与图8a，可以看出本发明使用的数据阻尼反演方法相比于传统的同时源反演结果，在连续性、层内光滑性和均匀性方面有明显的改善；图8b与图8c相比，可以看出，本发明实施例提供的数据阻尼全波形反演方法其结果与标准的序列源反演结果精度接近，表明数据阻尼全波形反演方法具体非常好有效性和可靠性。

图9a至图9e为本发明实施例提供的数据阻尼全波形反演方法运用到Marmousi模型中得到的反演结果与常规方法得到的反演结果的对比示意图。图9a为真实的速度模型示意图；图9b为反演所用的初始速度模型示意图；图9c是无数据阻尼的同时源全波形反演方法得到结果；图9d是本说明书实施例提供的数据阻尼全波形反演方法得到的结果；图9e是常规序列源反演结果。对比图9d和图9c可以看出，本发明使用的数据阻尼同时源全波形反演方法相比于传统的同时源反演结果，反演效果改善明显，尤其在深部断层位置、以及地层的构造信息成像方面有着较大的改善；对比图9d和图9e，可以看出，本发明使用的数据阻尼同时源全波形反演方法相比于常规的序列源反演方法，结果精度接近。综合图8a至图8d和图9a至图9e可以看出，本发明使用的数据阻尼同时源全波形反演方法可以适用于多种场景且均具有良好的精度。

本说明书实施例提供的数据阻尼全波形反演方法还能够应用于实际勘探。图10a至图10g为本发明实施例提供的数据阻尼全波形反演方法运用到某工区实际场景下得到的反演结果与常规方法运用于此场景下得到的反演结果的对比示意图。图10a为反演所用的初始速度模型示意图；图10b是无数据阻尼的同时源全波形反演方法得到的结果；图10c是本说明书实施例提供的数据阻尼的同时源全波形反演的结果；图10d是图10b所对应反演结果的逆时偏移成像剖面结果；图10e是图10c所对应反演结果的逆时偏移成像剖面结果；图10f是图10d所对应三个共成像的道集结果；图10g是图10e所对应的三个共成像的道集结果。对比图10b和图10c可以看出，本发明提供的数据阻尼同时源全波形反演方法相比于无串扰的同时源全波形反演方法，得到的反演结果低波数信息更加丰富，地层界面更加清晰；同时，对比图10d和图10e，也验证了使用本说明书实施例提供的数据阻尼全波形方法对速度场模型参数更新的正确性，即图10e中线框中的部分与图10d线框中的部分相比，连续性更好，层更加清晰；对比图10f和图10g，可以看出图10g中能量更加聚焦，即说明图10e的成像优于图10d的成像并且两种方法的共成像点道集拉平程度上也说明了本发明提出的方法可以获得更加可信的速度场模型参数。

如图11所示，本说明书实施例还提供一种数据阻尼全波形反演装置，包括：

目标函数构建模块10，用于基于观测数据、模拟数据和衰减函数，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；

正传波场获取模块20，用于获取震源的正传波场；

反传波场获取模块30，用于根据所述衰减函数获得所述震源的反传波场；

梯度计算模块40，用于根据所述正传波场、所述反传波场和所述目标函数关于速度场模型参数的梯度计算表达式，计算速度场模型参数的梯度；

迭代更新模块50，用于根据所述梯度对所述速度场模型参数进行迭代更新，直至所述数据残差达到设定阈值或者迭代更新次数达到预设值。

综上，本发明实施例提供的一种数据阻尼(同时源)全波形反演方法和装置，采用一阶速度—应力声波动方程，在一次正演模拟中完成多个震源子波的波场模拟，同时使用分段衰减函数作用于多震源反传源，通过相敏检测方法分别对正传和反传波场实施波场解耦计算，求取无串扰的梯度结果，从而得到高精度的反演结果。

如图12所示，为本文实施例提供的一种计算机设备，所述计算机设备1202可以包括一个或多个处理器1204，诸如一个或多个中央处理单元(CPU)，每个处理单元可以实现一个或多个硬件线程。计算机设备1202还可以包括任何存储器1206，其用于存储诸如代码、设置、数据等之类的任何种类的信息。非限制性的，比如，存储器1206可以包括以下任一项或多种组合：任何类型的RAM，任何类型的ROM，闪存设备，硬盘，光盘等。更一般地，任何存储器都可以使用任何技术来存储信息。进一步地，任何存储器可以提供信息的易失性或非易失性保留。进一步地，任何存储器可以表示计算机设备1202的固定或可移除部件。在一种情况下，当处理器1204执行被存储在任何存储器或存储器的组合中的相关联的指令时，计算机设备1202可以执行相关联指令的任一操作。计算机设备1202还包括用于与任何存储器交互的一个或多个驱动机构1208，诸如硬盘驱动机构、光盘驱动机构等。

计算机设备1202还可以包括输入/输出模块1210(I/O)，其用于接收各种输入(经由输入设备1212)和用于提供各种输出(经由输出设备1214))。一个具体输出机构可以包括呈现设备1216和相关联的图形用户接口(GUI)1218。在其他实施例中，还可以不包括输入/输出模块1210(I/O)、输入设备1212以及输出设备1214，仅作为网络中的一台计算机设备。计算机设备1302还可以包括一个或多个网络接口1220，其用于经由一个或多个通信链路1222与其他设备交换数据。一个或多个通信总线1224将上文所描述的部件耦合在一起。

通信链路1222可以以任何方式实现，例如，通过局域网、广域网(例如，因特网)、点对点连接等、或其任何组合。通信链路1222可以包括由任何协议或协议组合支配的硬连线链路、无线链路、路由器、网关功能、名称服务器等的任何组合。

对应于图1-图2中的方法，本文实施例还提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器运行时执行上述方法的步骤。

本文实施例还提供一种计算机可读指令，其中当处理器执行所述指令时，其中的程序使得处理器执行如图1至图2所示的方法。

应理解，在本文的各种实施例中，上述各过程的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本文实施例的实施过程构成任何限定。还应理解，在本文实施例中，术语“和/或”仅仅是一种描述关联对象的关联关系，表示可以存在三种关系。例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本文的范围。

所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，上述描述的系统、装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。在本文所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另外，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口、装置或单元的间接耦合或通信连接，也可以是电的，机械的或其它的形式连接。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本文实施例方案的目的。

另外，在本文各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以是两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本文的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本文各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

本文中应用了具体实施例对本文的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本文的方法及其核心思想；同时，对于本领域的一般技术人员，依据本文的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本文的限制。

Claims (10)

1.一种数据阻尼全波形反演方法，其特征在于，包括：

基于观测数据、模拟数据和衰减函数，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；

获取震源的正传波场；

根据所述衰减函数获得虚拟震源的反传波场，所述虚拟震源与所述震源相对应；

根据所述正传波场、所述反传波场和所述目标函数关于速度场模型参数的梯度计算表达式，计算所述速度场模型参数的梯度；

根据所述梯度对所述速度场模型参数进行迭代更新，直至所述数据残差达到设定阈值或者迭代更新次数达到预设值。

2.根据权利要求1所述的一种数据阻尼全波形反演方法，其特征在于，所述衰减函数为：

其中，α_γ为常数；m∈Z；T_p为最小积分区间。

3.根据权利要求2所述的一种数据阻尼全波形反演方法，其特征在于，所述目标函数为：

其中，h(α_γ,t)＝h′²(α_γ,t)；d^obs(X_r,Ω,t)为在检波点X_r处获得的所述观测数据；u(X_r,Ω,t；m)为在当前速度场模型参数m下获得的所述模拟数据；Ω表示角频率，t表示时间；x表示计算区域，并且

m为速度场模型参数；m(x)为计算区域x处的速度场模型参数。

4.根据权利要求2所述的一种数据阻尼全波形反演方法，其特征在于，所述震源通过多个位于不同位置的单频子波同时激发获得，所述震源为：

其中，

是表示编码信号中第i个信号的位置，

ω_i表示编码信号中第i个信号的角频率，Ω＝(ω₁,ω₂,…,ω_n)，n_s表示一个震源内编码信号的个数。

5.根据权利要求4所述的一种数据阻尼全波形反演方法，其特征在于，所述获取震源的正传波场，包括：

求解正传波场的模拟方程，得到所述正传波场，所述模拟方程为：

其中，

x表示计算区域，u(x,Ω,t；m)表示正传波场，X_s表示震源位置信息，L[·]是声波正演模拟算子。

6.根据权利要求5所述的一种数据阻尼全波形反演方法，其特征在于，所述根据所述衰减函数获得虚拟震源的反传波场，包括：

根据所述衰减函数，计算虚拟震源，所述虚拟震源为：

其中，反传源为：

根据所述虚拟震源，获得伴随方程，所述伴随方程为：

其中，

表示伴随算子，u′(x,Ω,t；m)表示反传波场；

求解所述伴随方程，获得所述虚拟震源的反传波场。

7.根据权利要求1所述的一种数据阻尼全波形反演方法，其特征在于，所述根据所述正传波场、所述反传波场和所述目标函数关于速度场模型参数的梯度计算表达式，计算速度场模型参数的梯度，包括：

对所述正传波场和所述反传波场分别进行无串扰解耦，得到波场变量；

根据所述波场变量和所述梯度计算表达式，计算所述速度场模型参数的梯度。

8.根据权利要求7所述的一种数据阻尼全波形反演方法，其特征在于，所述对所述正传波场和所述反传波场分别进行无串扰解耦，包括：

获取第一参考信号和第二参考信号，所述第一参考信号的振幅和频率均与所述第二参考信号的振幅和频率相同，所述第一参考信号的相位与所述第二参考信号的相位相差90°；

根据所述第一参考信号和所述第二参考信号，对所述正传波场和所述反传波场分别进行无串扰解耦。

9.一种数据阻尼全波形反演装置，其特征在于，包括：

目标函数构建模块，用于基于观测数据、模拟数据和衰减函数，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；

正传波场获取模块，用于获取震源的正传波场；

反传波场获取模块，用于根据所述衰减函数获得虚拟震源的反传波场，所述虚拟震源与所述震源相对应；

梯度计算模块，用于根据所述正传波场、所述反传波场和所述目标函数关于速度场模型参数的梯度计算表达式，计算速度场模型参数的梯度；

迭代更新模块，用于根据所述梯度对所述速度场模型参数进行迭代更新，直至所述数据残差达到设定阈值或者迭代更新次数达到预设值。

10.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至8任意一项所述数据阻尼全波形反演方法。

Images