CN116626751B - 基于多目标函数的黏弹性参数同步反演方法、装置和设备 - Google Patents
基于多目标函数的黏弹性参数同步反演方法、装置和设备 Download PDFInfo
- Publication number
- CN116626751B CN116626751B CN202310612233.6A CN202310612233A CN116626751B CN 116626751 B CN116626751 B CN 116626751B CN 202310612233 A CN202310612233 A CN 202310612233A CN 116626751 B CN116626751 B CN 116626751B
- Authority
- CN
- China
- Prior art keywords
- seismic data
- target
- source
- objective function
- center frequency
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000006870 function Effects 0.000 title claims abstract description 150
- 238000000034 method Methods 0.000 title claims abstract description 76
- 230000001360 synchronised effect Effects 0.000 title claims abstract description 36
- 238000004088 simulation Methods 0.000 claims description 37
- 230000008859 change Effects 0.000 claims description 24
- 238000004364 calculation method Methods 0.000 claims description 20
- 238000012545 processing Methods 0.000 claims description 6
- 238000001228 spectrum Methods 0.000 claims description 4
- 230000008901 benefit Effects 0.000 abstract description 7
- 230000001808 coupling effect Effects 0.000 abstract description 7
- 238000011161 development Methods 0.000 abstract description 7
- 238000010586 diagram Methods 0.000 description 12
- 238000012795 verification Methods 0.000 description 10
- 239000003245 coal Substances 0.000 description 9
- 230000008569 process Effects 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000003384 imaging method Methods 0.000 description 5
- 230000008878 coupling Effects 0.000 description 4
- 238000010168 coupling process Methods 0.000 description 4
- 238000005859 coupling reaction Methods 0.000 description 4
- 238000005457 optimization Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 3
- 238000003860 storage Methods 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000003491 array Methods 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000005755 formation reaction Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000021715 photosynthesis, light harvesting Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/30—Analysis
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Abstract
本发明涉及一种基于多目标函数的黏弹性参数同步反演方法、装置和设备，属于勘探和开发技术领域，该方法、装置和设备在每次反演迭代过程中，分别使用旅行时目标函数和中心频率目标函数计算速度梯度和Q梯度，从而根据计算得到的梯度进行迭代反演，例如，对速度和Q模型进行迭代反演。与基于波形差目标函数的多参数同步反演方法相比，由于构建了中心频率目标函数来计算品质因子Q梯度，充分考虑了品质因子Q对地震旅行时的影响较弱，而地震中心频率主要受Q影响的特性，可以压制参数串扰噪音，减弱速度和Q的耦合效应，降低黏弹性全波形反演的多解性，从而可以获得更加准确和可靠的速度和Q模型。
Description
技术领域
本发明涉及勘探和开发技术领域，具体涉及一种基于多目标函数的黏弹性参数同步反演方法、装置和设备。
背景技术
在勘探和开发技术领域(例如，石油勘探开发、煤矿勘探开发等)，为了在保证安全的前提下实现高效生产，对地下构造的成像精度和分辨率要求越来越高。以煤矿勘探开发为例，通过地球物理探测技术，建立三维地质模型，达到采煤工作面内部的透明化，可实现煤矿的高效安全生产。传统的煤矿巷道地震勘探方法主要利用槽波进行探测。然而，槽波法无法对煤层顶底板进行高精度成像，并在三维空间精细刻画地质异常体(如夹矸、断层和陷落柱等)的展布、位置和形态，难以满足采煤工作面透明化的需求。
相关技术中，通过基于波动方程的全波形地震反演来实现对地下构造的精准成像。其中，基于波动方程的全波形地震反演是一种高精度、高分辨率地球物理成像技术，近年来，在勘探地球物理领域获得广泛关注。该方法基于对波动方程的数值求解，可精确模拟地震波在复杂介质中的传播规律，通过波形拟合和最优化算法对地下介质模型进行迭代反演，可综合利用全波形和全波场信息，在多尺度空间构建模型，具有突破传统地震探测技术局限性的潜力。理论上，基于黏弹性波动方程，全波形反演可以构建高精度的弹性和衰减系数(品质因子Q)模型，对采煤工作面进行高精度成像，为煤矿勘探、开发和安全生产提供技术支撑和保障。
但是，不同物理参数对地震数据的影响相互叠加，衰减造成地震波能量耗散和速度频散，不准确的Q模型会导致速度反演结果的不准确，速度误差也会形成参数串扰噪音，投影到反演的Q模型中，导致不可靠的反演结果，对地下结构的理解和认识产生误导。在黏弹性介质的多参数同步反演中，速度和Q的耦合效应导致非常强的多解性，是黏弹性全波形反演的关键问题之一。因此，受参数耦合影响，黏弹性全波形反演具有非常强的多解性，难以获得准确的弹性参数和Q模型。
发明内容
有鉴于此，本发明的目的在于提供一种基于多目标函数的黏弹性参数同步反演方法、装置和设备，以克服目前受参数耦合影响，黏弹性全波形反演具有非常强的多解性，难以获得准确的弹性参数和Q模型的问题。
为实现以上目的，本发明采用如下技术方案：
一方面，一种基于多目标函数的黏弹性参数同步反演方法，包括：
获取真实地震数据和模拟地震数据，构建地震数据集；
对所述地震数据集中的数据进行去噪处理，得到目标地震数据集；其中，所述目标地震数据集中包括目标真实地震数据和目标模拟地震数据；
基于旅行时目标函数和目标地震数据集计算第一伴随源，根据第一伴随源进行第一伴随波场模拟；基于伴随状态法，通过第一正演模拟波场和第一模拟伴随波场计算速度梯度；以及，基于中心频率目标函数和目标地震数据集计算第二伴随源，根据第二伴随源进行第二伴随波场模拟；基于伴随状态法，通过第二正演模拟波场和第二模拟伴随波场计算品质因子梯度；
其中，所述旅行时目标函数和所述中心频率目标函数均为预先构建；所述第一模拟伴随波场为根据第一伴随源进行第一伴随波场模拟获得，所述第二模拟伴随波场为根据第二伴随源进行第二伴随波场模拟获得；
基于所述速度梯度和所述品质因子梯度，对速度和品质因子模型进行迭代反演。
可选的，所述旅行时目标函数，为：
其中，φCC为旅行时目标函数，△t为旅行时差，xr为接收器。
可选的，所述基于旅行时目标函数和目标地震数据集计算第一伴随源，包括：
根据所述目标真实地震数据和所述目标模拟地震数据的互相关，计算得到所述旅行时差；
根据所述旅行时差，计算得到第一伴随源。
可选的，所述旅行时差的计算公式为：
其中，τ为时间移动量，d(t+τ)为时间移动量后的模拟地震数据，dobs为真实地震数据。
可选的，所述根据所述旅行时差，计算得到第一伴随源，包括：
根据以下公式计算得到第一伴随源：
其中，为第一伴随源。
可选的，所述中心频率目标函数，为：
其中，φCF为中心频率目标函数，ωc和分别为目标模拟地震数据和目标真实地震数据的中心频率,xr为接收器。
可选的，所述基于中心频率目标函数和目标地震数据集计算第二伴随源，包括：
根据目标地震数据集计算目标模拟地震数据的中心频率ωc，以及，目标真实地震数据的中心频率
根据两个中心频率的差值，计算由品质因子扰动产生的中心频率目标函数的变化；以及，计算由品质因子扰动产生的目标模拟地震数据的中心频率的变化；
根据所述目标模拟地震数据的中心频率的变化，以及，中心频率目标函数的变化，计算得到频率域伴随源；
根据所述频率域伴随源和预设伴随源基准，计算得到所述第二伴随源；其中，所述预设伴随源基准为根据全波形反演方法计算得到的计算品质因子梯度时的伴随源。
可选的，所述目标模拟地震数据的中心频率ωc的计算公式为：
其中，A(ω)为振幅谱；ω为角频率；
所述由品质因子扰动产生的中心频率目标函数的变化ΔφCF的计算公式为：
其中，由品质因子扰动产生的目标模拟地震数据的中心频率的变化Δωc的计算公式为：
其中，R[·]代表实部，为目标模拟地震数据的傅里叶变换；
所述频率域伴随源的表示公式为：
所述根据所述频率域伴随源和预设伴随源基准，计算得到所述第二伴随源，包括：将所述预设伴随源基准中的替换为/>得到所述中心频率目标函数计算品质因子梯度的伴随源公式，作为所述第二伴随源。
又一方面，一种基于多目标函数的黏弹性参数同步反演装置，包括：
获取模块，用于获取真实地震数据和模拟地震数据，构建地震数据集；
数据处理模块，用于对所述地震数据集中的数据进行去噪处理，得到目标地震数据集；其中，所述目标地震数据集中包括目标真实地震数据和目标模拟地震数据；
计算反演模块，用于基于旅行时目标函数和目标地震数据集计算第一伴随源，根据第一伴随源进行第一伴随波场模拟；基于伴随状态法，通过第一正演模拟波场和第一模拟伴随波场计算速度梯度；以及，基于中心频率目标函数和目标地震数据集计算第二伴随源，根据第二伴随源进行第二伴随波场模拟；基于伴随状态法，通过第二正演模拟波场和第二模拟伴随波场计算品质因子梯度；
其中，所述旅行时目标函数和所述中心频率目标函数均为预先构建；所述第一模拟伴随波场为根据第一伴随源进行第一伴随波场模拟获得，所述第二模拟伴随波场为根据第二伴随源进行第二伴随波场模拟获得；
基于所述速度梯度和所述品质因子梯度，对速度和品质因子模型进行迭代反演。
又一方面，一种基于多目标函数的黏弹性参数同步反演设备，包括处理器和存储器，所述处理器与存储器相连：
其中，所述处理器，用于调用并执行所述存储器中存储的程序；
所述存储器，用于存储所述程序，所述程序至少用于执行上述任一项所述的基于多目标函数的黏弹性参数同步反演方法。
本发明提供的技术方案，至少具备如下有益效果：
在每次反演迭代过程中，分别使用旅行时目标函数和中心频率目标函数计算速度梯度和Q梯度，从而根据计算得到的黏弹性参数进行迭代反演，例如，对速度和Q模型进行迭代反演。与基于波形差目标函数的多参数同步反演方法相比，由于构建了中心频率目标函数来计算品质因子Q梯度，充分考虑了品质因子Q对地震旅行时的影响较弱，而地震中心频率主要受Q影响的特性，可以压制参数串扰噪音，减弱速度和Q的耦合效应，降低黏弹性全波形反演的多解性，从而可以以供获得更加准确和可靠的速度和Q模型。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明一实施例提供的一种基于波形差目标函数的多参数同步反演方法的流程示意图；
图2为本发明一实施例提供的一种基于多目标函数的黏弹性参数同步反演方法的流程示意图；
图3为本发明又一实施例提供的一种基于多目标函数的黏弹性参数同步反演方法的流程示意图；
图4为本发明一实施例提供的一种基于多目标函数的黏弹性参数同步反演装置的结构示意图；
图5是本发明一实施例提供的一种基于多目标函数的黏弹性参数同步反演设备的结构示意图；
图6为本验证实施例提供的一种起伏地表条件下的真实黏弹性介质模型；
图7为本验证实施例提供的一种波形差目标函数计算各梯度值示意图；
图8为本验证实施例提供的一种分别使用旅行时和中心频率目标函数计算的速度和品质因子梯度示意图；
图9为本验证实施例提供的一种波形差目标函数反演各模型示意图；
图10为本验证实施例提供的一种分别使用旅行时和中心频率目标函数反演的各模型示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将对本发明的技术方案进行详细的描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施方式，都属于本发明所保护的范围。
如背景技术记载，相关技术中，通过基于波动方程的全波形地震反演来实现对地下构造的精准成像。
基于此，本申请对现有技术基于波形差目标函数的多参数同步反演方法进行说明：
全波形反演方法本质上是一个最优化过程，通过将模拟数据与实际观测数据进行波形匹配，获得最优的地下介质模型。在传统全波形反演算法中，首先建立波形差目标函数，测量模拟数据与观测数据之间波形残差。该目标函数φ可以写为：
其中，m代表模型参数，d和dobs分别为模拟数据和观测数据，T为最大记录时间，xr代表接收器位置。然后采用局部优化算法(如梯度法或拟牛顿法)，对模型参数进行迭代更新，最终获得目标函数的最小值和最优介质模型。梯度计算是全波形反演关键步骤之一，为避免雅克比矩阵的直接计算，一般基于伴随状态法，通过正演模拟波场和伴随波场的互相关计算梯度。在黏弹性介质中，基于标准线性固体模型，体积模量κ和剪切模量μ的梯度公式可以推导为：
其中，u为位移场，为位移场散度，/>和/>为伴随应变张量ε*的迹和偏差部分，P为正演应变偏量，/>表示褶积。为计算伴随波场，伴随源为：
基于Kolsk-Futterman模型和链式法则，体积模量和剪切模量品质因子(Qκ和Qμ)的梯度可以表示为：
与计算模量梯度不同，计算品质因子梯度的伴随波场，伴随源可表示为：
其中，ω和ω0为角频率和参考角频率，sgn为符号函数，为公式(3)的傅里叶变换。基于链式法则，纵波速度VP、横波速度VS、纵波品质因子QP和横波品质因子QS的梯度公式可以表示为：
其中，v＝VP/VS，σ＝QP/QS。基于黏弹性波动方程的数值求解，理论上，全波形反演技术可以构建高精度的速度和Q模型。
图1为本发明一实施例提供的一种基于波形差目标函数的多参数同步反演方法的流程示意图，参阅图1，观测地震数据为实际观测数据，预先设定初始模型，通过初始模型进行正演模拟，得到模拟地震数据(即模拟数据)，根据模拟地震数据和实际观测数据计算得到公式(1)，从而根据公式(2)和公式(3)计算伴随源，在计算得到伴随源后，进行伴随波场模拟，从而根据公式(4)至公式(7)计算得到速度和Q梯度。在计算得到速度和Q梯度后，根据现有技术计算搜索方向，从而线性搜索计算步长，对速度和Q模型进行更新，并迭代判断更新后的速度和Q模型是否收敛，在收敛时，停止反演，得到反演结果；在不收敛时，再次通过初始模型进行正演模拟，得到模拟数据，更新模拟数据，重新进行速度和Q模型的更新。
但是，不同物理参数对地震数据的影响相互叠加，衰减造成地震波能量耗散和速度频散，不准确的Q模型会导致速度反演结果的不准确，速度误差也会形成参数串扰噪音，投影到反演的Q模型中，导致不可靠的反演结果，对地下结构的理解和认识产生误导。在黏弹性介质的多参数同步反演中，速度和Q的耦合效应导致非常强的多解性，是黏弹性全波形反演的关键问题之一。
基于此，本发明实施例提供一种基于多目标函数的黏弹性参数同步反演方法、装置和设备。
图2为本发明一实施例提供的一种基于多目标函数的黏弹性参数同步反演方法的流程示意图，参阅图2，本实施例提供的基于多目标函数的黏弹性参数同步反演方法，可以包括以下步骤：
步骤S21、获取真实地震数据和模拟地震数据，构建地震数据集；
步骤S22、对所述地震数据集中的数据进行去噪处理，得到目标地震数据集；其中，所述目标地震数据集中包括目标真实地震数据和目标模拟地震数据；
步骤S23、基于旅行时目标函数和目标地震数据集计算第一伴随源，根据第一伴随源进行第一伴随波场模拟；基于伴随状态法，通过第一正演模拟波场和第一模拟伴随波场计算速度梯度；以及，基于中心频率目标函数和目标地震数据集计算第二伴随源，根据第二伴随源进行第二伴随波场模拟；基于伴随状态法，通过第二正演模拟波场和第二模拟伴随波场计算品质因子梯度；
其中，所述旅行时目标函数和所述中心频率目标函数均为预先构建；所述第一模拟伴随波场为根据第一伴随源进行第一伴随波场模拟获得，所述第二模拟伴随波场为根据第二伴随源进行第二伴随波场模拟获得；
步骤S24、基于所述速度梯度和所述品质因子梯度，对速度和品质因子模型进行迭代反演。
参阅上述对现有技术的记载，可以预设初始模型，通过初始模型进行黏弹性正演模拟，从而得到模拟地震数据。在获取到真实地震数据和模拟地震数据后，构建地震数据集，对数据集进行去噪处理后，得到目标地震数据集。
从而根据预设的旅行时目标函数和中心频率目标函数进行速度梯度和所述品质因子梯度迭代反演。具体的，可以基于所述速度梯度和所述品质因子梯度，以及最优化算法，对速度和品质因子模型进行迭代反演。
可以理解的是，采用本实施例提供的技术方案，在每次反演迭代过程中，分别使用旅行时目标函数和中心频率目标函数计算速度梯度和Q梯度，从而根据计算得到的黏弹性参数进行迭代反演，例如，对速度和Q模型进行迭代反演。与基于波形差目标函数的多参数同步反演方法相比，由于构建了中心频率目标函数来计算品质因子Q梯度，充分考虑了品质因子Q对地震旅行时的影响较弱，而地震中心频率主要受Q影响的特性，可以压制参数串扰噪音，减弱速度和Q的耦合效应，降低黏弹性全波形反演的多解性，从而可以以供获得更加准确和可靠的速度和Q模型。
在一些实施例中，所述旅行时目标函数，为：
其中，φCC为旅行时目标函数，△t为旅行时差，xr为地震接收器。
在一些实施例中，所述基于旅行时目标函数和目标地震数据集计算第一伴随源，包括：
根据所述目标真实地震数据和所述目标模拟地震数据的互相关，计算得到所述旅行时差；
根据所述旅行时差，计算得到第一伴随源。
在一些实施例中，所述旅行时差的计算公式为：
其中，τ为时间移动量，d(t+τ)为时间移动量后的模拟地震数据，dobs为真实地震数据。
在一些实施例中，所述根据所述旅行时差，计算得到第一伴随源，包括：
根据以下公式计算得到第一伴随源：
其中，为第一伴随源。
在一些实施例中，所述中心频率目标函数，为：
其中，φCF为中心频率目标函数，ωc和分别为目标模拟地震数据和目标真实地震数据的中心频率,xr为地震接收器。
在一些实施例中，所述基于中心频率目标函数和目标地震数据集计算第二伴随源，包括：
根据目标地震数据集计算目标模拟地震数据的中心频率ωc，以及，目标真实地震数据的中心频率
根据两个中心频率的差值，计算由品质因子扰动产生的中心频率目标函数的变化；以及，计算由品质因子扰动产生的目标模拟地震数据的中心频率的变化；
根据所述目标模拟地震数据的中心频率的变化，以及，中心频率目标函数的变化，计算得到频率域伴随源；
根据所述频率域伴随源和预设伴随源基准，计算得到所述第二伴随源；其中，所述预设伴随源基准为根据全波形反演方法计算得到的计算品质因子梯度时的伴随源。
在一些实施例中，所述目标模拟地震数据的中心频率ωc的计算公式为：
其中，A(ω)为振幅谱；ω为角频率；
所述由品质因子扰动产生的中心频率目标函数的变化ΔφCF的计算公式为：
其中，由品质因子扰动产生的目标模拟地震数据的中心频率的变化Δωc的计算公式为：
其中，R[·]代表实部，为目标模拟地震数据的傅里叶变换,/>为模拟地震数据的扰动；
所述频率域伴随源的表示公式为：
所述根据所述频率域伴随源和预设伴随源基准，计算得到所述第二伴随源，包括：将所述预设伴随源基准中的替换为/>得到所述中心频率目标函数计算品质因子梯度的伴随源公式，作为所述第二伴随源。
例如，可以将公式(5)中的替换为/>便可以得到中心频率目标函数计算Q梯度的伴随源公式。
在得到两个伴随源后，分别进行正演波场和伴随波场的模拟，将正演波场和伴随波场进行互相关计算，即可得到速度和Q的梯度。
图3为本发明又一实施例提供的一种基于多目标函数的黏弹性参数同步反演方法的流程示意图，参阅图3，在得到观测地震数据(真实地震数据)和模拟地震数据后，进行数据预处理(去噪等)，从而根据去噪后的数据分别与旅行时目标函数、中心频率目标函数相结合，分别计算第一伴随源、第二伴随源；并分别进行第一伴随波场模拟、第二伴随波场模拟，分别进行速度梯度计算核Q梯度计算，将所述速度梯度和所述品质因子梯度作为黏弹性参数，根据所述黏弹性参数进行迭代反演。在得到黏弹性参数后，进行梯度后处理(如，采用高斯平滑等降低不确定性)；根据拟牛顿最优化算法计算搜索方向，计算步长，更新速度和Q模型的更新，并判断该模型是否收敛，在模型收敛时，反演结束；在模型不收敛时，重新进行黏弹性正演模拟，进行迭代。
基于一个总的发明构思，本发明实施例还提供一种基于多目标函数的黏弹性参数同步反演装置，用于实现上述方法实施例。
图4为本发明一实施例提供的一种基于多目标函数的黏弹性参数同步反演装置的结构示意图，参阅图4，本实施例提供的装置，可以包括以下结构：
获取模块41，用于获取真实地震数据和模拟地震数据，构建地震数据集；
数据处理模块42，用于对所述地震数据集中的数据进行去噪处理，得到目标地震数据集；其中，所述目标地震数据集中包括目标真实地震数据和目标模拟地震数据；
计算反演模块43，用于基于旅行时目标函数和目标地震数据集计算第一伴随源，根据第一伴随源进行第一伴随波场模拟；基于伴随状态法，通过第一正演模拟波场和第一模拟伴随波场计算速度梯度；以及，基于中心频率目标函数和目标地震数据集计算第二伴随源，根据第二伴随源进行第二伴随波场模拟；基于伴随状态法，通过第二正演模拟波场和第二模拟伴随波场计算品质因子梯度；
其中，所述旅行时目标函数和所述中心频率目标函数均为预先构建；所述第一模拟伴随波场为根据第一伴随源进行第一伴随波场模拟获得，所述第二模拟伴随波场为根据第二伴随源进行第二伴随波场模拟获得；
将所述速度梯度和所述品质因子梯度作为黏弹性参数，根据所述黏弹性参数进行迭代反演。
关于上述实施例中的装置，其中各个模块执行操作的具体方式已经在有关该方法的实施例中进行了详细描述，此处将不做详细阐述说明。
可以理解的是，采用本实施例提供的技术方案，在每次反演迭代过程中，分别使用旅行时目标函数和中心频率目标函数计算速度梯度和Q梯度，从而根据计算得到的黏弹性参数进行迭代反演，例如，对速度和Q模型进行迭代反演。与基于波形差目标函数的多参数同步反演方法相比，由于构建了中心频率目标函数来计算品质因子Q梯度，充分考虑了品质因子Q对地震旅行时的影响较弱，而地震中心频率主要受Q影响的特性，可以压制参数串扰噪音，减弱速度和Q的耦合效应，降低黏弹性全波形反演的多解性，从而可以以供获得更加准确和可靠的速度和Q模型。
基于一个总的发明构思，本发明实施例还提供一种基于多目标函数的黏弹性参数同步反演设备，用于实现上述方法实施例。
图5是本发明一实施例提供的一种基于多目标函数的黏弹性参数同步反演设备的结构示意图。如图5所示，本实施例的基于多目标函数的黏弹性参数同步反演设备包括处理器51和存储器52，处理器51与存储器52相连。其中，处理器51用于调用并执行所述存储器52中存储的程序；存储器52用于存储所述程序，所述程序至少用于执行以上实施例中的基于多目标函数的黏弹性参数同步反演方法。
本申请实施例提供的基于多目标函数的黏弹性参数同步反演设备的具体实施方案可以参考以上任意实施例的基于多目标函数的黏弹性参数同步反演方法的实施方式，此处不再赘述。
为了对本实施例提供的技术方案的效果进行说明，本实施例提供验证实施例。该验证实施例为数值模拟实验。
采用本申请上述实施例提供的基于多目标函数的黏弹性参数同步反演方法，设计黏弹性介质全波形反演的数值模拟实验，验证本发明方法的有效性和优势。
图6为本验证实施例提供的一种起伏地表条件下的真实黏弹性介质模型；a)为黏弹性介质速度模型，b)为黏弹性介质品质因子模型，黑色曲线代表起伏地表。VP和VS速度结构的背景数值为2400m/s和800m/s。真实VP和VS模型在浅层包含-5％的速度异常区域，如图6a所示。初始VP和VS模型不包含该浅层速度异常区域。真实QP和QS模型的背景数值为150，在浅层包含数值为20的强衰减区域，如图6b所示。初始QP和QS模型不包含该浅层衰减异常区域。
将24个震源和240个接收器均匀布置于起伏地表，使用真实模型进行正演模拟获得观测地震数据，地震子波为主频30Hz的雷克子波。首先使用传统波形差目标函数对速度和品质因子的梯度进行计算，如图7所示(图7为本验证实施例提供的一种波形差目标函数计算各梯度值示意图，a)为波形差目标函数计算VS梯度，b)为波形差目标函数计算QS梯度，c)为波形差目标函数计算VP梯度，d)为波形差目标函数计算QP梯度)。受参数串扰影响，在计算VS和VP梯度中，存在较弱的参数串扰噪音，而在计算QS和QP梯度中，由速度误差导致的噪音非线明显，如图7中箭头所示。图8为本验证实施例提供的一种分别使用旅行时和中心频率目标函数计算的速度和品质因子梯度示意图，a)为旅行时目标函数计算VS梯度，b)为中心频率目标函数计算QS梯度，c)为旅行时目标函数计算VP梯度，d)为中心频率目标函数计算QP梯度。与图7中的梯度相比，在计算的梯度中，速度和衰减异常的位置更加清晰，受参数耦合影响非常弱，如箭头所示。
图9为本验证实施例提供的一种波形差目标函数反演各模型示意图，a)波形差目标函数反演VS模型，b)波形差目标函数反演QS模型，c)波形差目标函数反演VP模型，d)波形差目标函数反演QP模型。参阅图9，使用传统波形差目标函数对速度和品质因子进行同步反演的结果。受衰减误差的影响，反演的速度结构包含非常强的噪音，如箭头所示。这些噪音会对地下结构的认知产生误导。在反演的品质因子模型中，浅层的衰减异常区域则非常模糊，与真实品质因子模型的数值相差较远。图10为本验证实施例提供的一种分别使用旅行时和中心频率目标函数反演的各模型示意图，a)旅行时目标函数反演VS模型，b)中心频率目标函数反演QS模型，c)旅行时目标函数反演VP模型，d)中心频率目标函数反演QP模型。参阅图10，为使用旅行时和中心频率目标函数对速度和品质因子进行同步反演的结果。与图9所示的速度反演结果相比，参数串扰噪音被压制，速度结构更加准确。在反演的品质因子模型中，浅层的衰减异常区域清晰可见，与真实品质因子模型的数值更加接近。
根据上述数值模拟实验结果和分析，可见基于本发明提出的多目标函数黏弹性参数同步反演方法，可以压制参数耦合的影响，降低反演的非线性，更加准确的构建速度和品质因子模型。
因此，在黏弹性全波形反演中，本申请提供的基于多目标函数的速度和品质因子同步反演方法，最大的优势在于，当分别旅行时和中心频率目标函数计算速度和品质因子的梯度，可以压制由参数耦合效应导致的串扰噪音，降低反演的多解性和非线性，获得更加准确和可靠的速度和品质因子模型。
可以理解的是，上述各实施例中相同或相似部分可以相互参考，在一些实施例中未详细说明的内容可以参见其他实施例中相同或相似的内容。
需要说明的是，在本发明的描述中，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性。此外，在本发明的描述中，除非另有说明，“多个”的含义是指至少两个。
流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本发明的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本发明的实施例所属技术领域的技术人员所理解。
应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。
本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。
此外，在本发明各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。
上述提到的存储介质可以是只读存储器，磁盘或光盘等。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
Claims (7)
1.一种基于多目标函数的黏弹性参数同步反演方法，其特征在于，包括：
获取真实地震数据和模拟地震数据，构建地震数据集；
对所述地震数据集中的数据进行去噪处理，得到目标地震数据集；其中，所述目标地震数据集中包括目标真实地震数据和目标模拟地震数据；
基于旅行时目标函数和目标地震数据集计算第一伴随源，根据第一伴随源进行第一伴随波场模拟；基于伴随状态法，通过第一正演模拟波场和第一模拟伴随波场计算速度梯度；以及，基于中心频率目标函数和目标地震数据集计算第二伴随源，根据第二伴随源进行第二伴随波场模拟；基于伴随状态法，通过第二正演模拟波场和第二模拟伴随波场计算品质因子梯度；
其中，所述旅行时目标函数和所述中心频率目标函数均为预先构建；所述第一模拟伴随波场为根据第一伴随源进行第一伴随波场模拟获得，所述第二模拟伴随波场为根据第二伴随源进行第二伴随波场模拟获得；
基于所述速度梯度和所述品质因子梯度，对速度和品质因子模型进行迭代反演；
其中，所述中心频率目标函数，为：
其中，φCF为中心频率目标函数，ωc和分别为目标模拟地震数据和目标真实地震数据的中心频率,xr为地震接收器；
所述基于中心频率目标函数和目标地震数据集计算第二伴随源，包括：
根据目标地震数据集计算目标模拟地震数据的中心频率ωc，以及，目标真实地震数据的中心频率
根据两个中心频率的差值，计算由品质因子扰动产生的中心频率目标函数的变化；以及，计算由品质因子扰动产生的目标模拟地震数据的中心频率的变化；
根据所述目标模拟地震数据的中心频率的变化，以及，中心频率目标函数的变化，计算得到频率域伴随源；
根据所述频率域伴随源和预设伴随源基准，计算得到所述第二伴随源；其中，所述预设伴随源基准为根据全波形反演方法计算得到的计算品质因子梯度时的伴随源；
所述目标模拟地震数据的中心频率ωc的计算公式为：
其中，A(ω)为振幅谱；ω为角频率；
所述由品质因子扰动产生的中心频率目标函数的变化ΔφCF的计算公式为：
其中，由品质因子扰动产生的目标模拟地震数据的中心频率的变化Δωc的计算公式为：
其中，R[·]代表实部，为目标模拟地震数据的傅里叶变换；
所述频率域伴随源的表示公式为：
所述根据所述频率域伴随源和预设伴随源基准，计算得到所述第二伴随源，包括：将所述预设伴随源基准中的替换为/>得到所述中心频率目标函数计算品质因子梯度的伴随源公式，作为所述第二伴随源。
2.根据权利要求1所述的方法，其特征在于，所述旅行时目标函数，为：
其中，φCC为旅行时目标函数，△t为旅行时差，xr为地震接收器。
3.根据权利要求2所述的方法，其特征在于，所述基于旅行时目标函数和目标地震数据集计算第一伴随源，包括：
根据所述目标真实地震数据和所述目标模拟地震数据的互相关，计算得到所述旅行时差；
根据所述旅行时差，计算得到第一伴随源。
4.根据权利要求3所述的方法，其特征在于，所述旅行时差的计算公式为：
其中，τ为时间移动量，d(t+τ)为时间移动量后的模拟地震数据，dobs为真实地震数据。
5.根据权利要求3所述的方法，其特征在于，所述根据所述旅行时差，计算得到第一伴随源，包括：
根据以下公式计算得到第一伴随源：
其中，为第一伴随源。
6.一种基于多目标函数的黏弹性参数同步反演装置，其特征在于，包括：
获取模块，用于获取真实地震数据和模拟地震数据，构建地震数据集；
数据处理模块，用于对所述地震数据集中的数据进行去噪处理，得到目标地震数据集；其中，所述目标地震数据集中包括目标真实地震数据和目标模拟地震数据；
计算反演模块，用于基于旅行时目标函数和目标地震数据集计算第一伴随源，根据第一伴随源进行第一伴随波场模拟；基于伴随状态法，通过第一正演模拟波场和第一模拟伴随波场计算速度梯度；以及，基于中心频率目标函数和目标地震数据集计算第二伴随源，根据第二伴随源进行第二伴随波场模拟；基于伴随状态法，通过第二正演模拟波场和第二模拟伴随波场计算品质因子梯度；
其中，所述旅行时目标函数和所述中心频率目标函数均为预先构建；所述第一模拟伴随波场为根据第一伴随源进行第一伴随波场模拟获得，所述第二模拟伴随波场为根据第二伴随源进行第二伴随波场模拟获得；所述中心频率目标函数，为：
其中，φCF为中心频率目标函数，ωc和分别为目标模拟地震数据和目标真实地震数据的中心频率,xr为地震接收器；
基于所述速度梯度和所述品质因子梯度，对速度和品质因子模型进行迭代反演；
所述基于中心频率目标函数和目标地震数据集计算第二伴随源，包括：
根据目标地震数据集计算目标模拟地震数据的中心频率ωc，以及，目标真实地震数据的中心频率
根据两个中心频率的差值，计算由品质因子扰动产生的中心频率目标函数的变化；以及，计算由品质因子扰动产生的目标模拟地震数据的中心频率的变化；
根据所述目标模拟地震数据的中心频率的变化，以及，中心频率目标函数的变化，计算得到频率域伴随源；
根据所述频率域伴随源和预设伴随源基准，计算得到所述第二伴随源；其中，所述预设伴随源基准为根据全波形反演方法计算得到的计算品质因子梯度时的伴随源；
所述目标模拟地震数据的中心频率ωc的计算公式为：
其中，A(ω)为振幅谱；ω为角频率；
所述由品质因子扰动产生的中心频率目标函数的变化ΔφCF的计算公式为：
其中，由品质因子扰动产生的目标模拟地震数据的中心频率的变化Δωc的计算公式为：
其中，R[·]代表实部，为目标模拟地震数据的傅里叶变换；
所述频率域伴随源的表示公式为：
所述根据所述频率域伴随源和预设伴随源基准，计算得到所述第二伴随源，包括：将所述预设伴随源基准中的替换为/>得到所述中心频率目标函数计算品质因子梯度的伴随源公式，作为所述第二伴随源。
7.一种基于多目标函数的黏弹性参数同步反演设备，其特征在于，包括处理器和存储器，所述处理器与存储器相连：
其中，所述处理器，用于调用并执行所述存储器中存储的程序；
所述存储器，用于存储所述程序，所述程序至少用于执行权利要求1-5任一项所述的基于多目标函数的黏弹性参数同步反演方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310612233.6A CN116626751B (zh) | 2023-05-26 | 2023-05-26 | 基于多目标函数的黏弹性参数同步反演方法、装置和设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310612233.6A CN116626751B (zh) | 2023-05-26 | 2023-05-26 | 基于多目标函数的黏弹性参数同步反演方法、装置和设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116626751A CN116626751A (zh) | 2023-08-22 |
| CN116626751B true CN116626751B (zh) | 2024-03-19 |
Family
ID=87613014
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310612233.6A Active CN116626751B (zh) | 2023-05-26 | 2023-05-26 | 基于多目标函数的黏弹性参数同步反演方法、装置和设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116626751B (zh) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105467444A (zh) * | 2015-12-10 | 2016-04-06 | 中国石油天然气集团公司 | 一种弹性波全波形反演方法及装置 |
| CN108873066A (zh) * | 2018-06-26 | 2018-11-23 | 中国石油大学(华东) | 弹性介质波动方程反射波旅行时反演方法 |
| CN111045077A (zh) * | 2019-12-20 | 2020-04-21 | 核工业北京地质研究院 | 一种陆地地震数据的全波形反演方法 |
| CN113376695A (zh) * | 2021-06-11 | 2021-09-10 | 中国矿业大学 | 一种适用于煤层底板复杂陷落柱的全波形反演方法 |
| CN113504566A (zh) * | 2021-06-01 | 2021-10-15 | 南方海洋科学与工程广东省实验室(湛江) | 基于波动方程的地震反演方法、系统、装置及介质 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11360224B2 (en) * | 2019-05-03 | 2022-06-14 | Exxonmobil Upstream Research Company | Inversion, migration, and imaging related to isotropic wave-mode- independent attenuation |
-
2023
- 2023-05-26 CN CN202310612233.6A patent/CN116626751B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105467444A (zh) * | 2015-12-10 | 2016-04-06 | 中国石油天然气集团公司 | 一种弹性波全波形反演方法及装置 |
| CN108873066A (zh) * | 2018-06-26 | 2018-11-23 | 中国石油大学(华东) | 弹性介质波动方程反射波旅行时反演方法 |
| CN111045077A (zh) * | 2019-12-20 | 2020-04-21 | 核工业北京地质研究院 | 一种陆地地震数据的全波形反演方法 |
| CN113504566A (zh) * | 2021-06-01 | 2021-10-15 | 南方海洋科学与工程广东省实验室(湛江) | 基于波动方程的地震反演方法、系统、装置及介质 |
| CN113376695A (zh) * | 2021-06-11 | 2021-09-10 | 中国矿业大学 | 一种适用于煤层底板复杂陷落柱的全波形反演方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于多目标函数的黏弹性全波形反演;梁旭 等;《煤田地质与勘探》;第51卷(第4期);152-163 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116626751A (zh) | 2023-08-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| RU2693495C1 (ru) | Полная инверсия волнового поля с компенсацией показателя качества | |
| Taillandier et al. | First-arrival traveltime tomography based on the adjoint-state method | |
| US10002211B2 (en) | Artifact reduction in iterative inversion of geophysical data | |
| EP3073296A1 (en) | Full waveform inversion method for seismic data processing using preserved amplitude reverse time migration | |
| CN103293552B (zh) | 一种叠前地震资料的反演方法及系统 | |
| US10520619B2 (en) | FWI model domain angle stacks with amplitude preservation | |
| EP3259620B1 (en) | Multistage full wavefield inversion process that generates a multiple free data set | |
| Cho et al. | Generalized multiscale finite elements for simulation of elastic-wave propagation in fractured media | |
| CN103713315A (zh) | 一种地震各向异性参数全波形反演方法及装置 | |
| Maurer et al. | Receiver-coupling effects in seismic waveform inversions | |
| EA032186B1 (ru) | Сейсмическая адаптивная фокусировка | |
| EP3243089B1 (en) | Method for obtaining estimates of a model parameter so as to characterise the evolution of a subsurface volume over a period of time | |
| Raknes et al. | A numerical study of 3D elastic time-lapse full-waveform inversion using multicomponent seismic data | |
| Lamert et al. | Imaging disturbance zones ahead of a tunnel by elastic full-waveform inversion: Adjoint gradient based inversion vs. parameter space reduction using a level-set method | |
| Lin et al. | Double-difference traveltime tomography with edge-preserving regularization and a priori interfaces | |
| CN116626751B (zh) | 基于多目标函数的黏弹性参数同步反演方法、装置和设备 | |
| Dong et al. | Correlation‐based reflection waveform inversion by one‐way wave equations | |
| Gheymasi et al. | Robust total-variation based geophysical inversion using split Bregman and proximity operators | |
| Jaimes-Osorio et al. | Amplitude variation with offset inversion using acoustic-elastic local solver | |
| Raknes et al. | Challenges and solutions for performing 3D time-domain elastic full-waveform inversion | |
| AU2016220145A1 (en) | Black hole boundary conditions | |
| Riedel et al. | Elastic waveform inversion in the frequency domain for an application in mechanized tunneling | |
| Li et al. | 2D high-resolution crosswell seismic traveltime tomography | |
| Lin et al. | Target-oriented waveform inversion based on Marchenko redatumed data | |
| Liu et al. | An inverted heterogeneous velocity model for microseismic source location in deep buried tunnels |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |