[go: up one dir, main page]

CN119224832B - 基于自伴随vti方程的全波形反演方法、装置、介质及设备 - Google Patents

基于自伴随vti方程的全波形反演方法、装置、介质及设备

Info

Publication number
CN119224832B
CN119224832B CN202411218000.9A CN202411218000A CN119224832B CN 119224832 B CN119224832 B CN 119224832B CN 202411218000 A CN202411218000 A CN 202411218000A CN 119224832 B CN119224832 B CN 119224832B
Authority
CN
China
Prior art keywords
vti
accompanying
self
equation
wave equation
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
Application number
CN202411218000.9A
Other languages
English (en)
Other versions
CN119224832A (zh
Inventor
姚刚
吴博
吴迪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum Beijing
Original Assignee
China University of Petroleum Beijing
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Petroleum Beijing filed Critical China University of Petroleum Beijing
Priority to CN202411218000.9A priority Critical patent/CN119224832B/zh
Publication of CN119224832A publication Critical patent/CN119224832A/zh
Application granted granted Critical
Publication of CN119224832B publication Critical patent/CN119224832B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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. for interpretation or for event detection
    • 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. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G01V1/305—Travel times
    • 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. for interpretation or for event detection
    • G01V1/32—Transforming one recording into another or one representation into another
    • 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. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明涉及一种基于自伴随VTI方程的全波形反演方法、装置、介质及设备，方法包括：拆分VTI声波方程和VTI弹性波方程，得到拆分的VTI声波方程和VTI弹性波方程；构建辅助变量p和q并带入拆分的VTI声波方程和VTI弹性波方程中，得到自伴随的VTI声波方程和自伴随的VTI弹性波方程；求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，计算得到正传波场u和预测数据dpre；利用目标函数和预测数据dpre计算残差dres；求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，反传残差dres计算伴随波场u*；基于目标函数对模型参数的梯度公式，利用正传波场u和伴随波场u*计算梯度计算步长α并通过梯度和迭代更新公式更新模型m，判断目标函数是否收敛。

Description

基于自伴随VTI方程的全波形反演方法、装置、介质及设备
技术领域
本发明涉及一种基于自伴随VTI方程的全波形反演方法、装置、介质及设备，属于地质勘探技术领域。
背景技术
垂直对称轴的横向各向同性(VTI)介质的全波形反演可以对VTI各向异性区域实现有效的地下介质建模。全波形反演通常在一个较好的初始模型上通过局部最优化算法进行迭代更新来实现地下介质的高精度反演。局部最优化算法进行迭代更新需要求取目标函数的梯度作为更新方向，而梯度的计算则需要准确的实施正传波场和伴随波场的数值模拟。
目前，利用VTI的弹性波方程可以实现VTI弹性波场的正传波场模拟。为了模拟VTI声波波场，通常通过设置弹性波方程中横波速度沿对称轴方向为0(Duveneck et al.,2008)，来得到VTI声波各向异性波动方程，准确的实现VTI声波波场的正传波场模拟(图1(a))。在实现VTI各向异性介质的全波形反演过程中，通常使用伴随状态法来求取目标函数关于模型参数的梯度。伴随状态法需要推导伴随波动方程和伴随波场。但是，目前普遍使用的VTI弹性波动方程和VTI声波波动方程都不是自对称的。因此，直接通过转置波动方程算子得到的伴随波动方程存在明显问题。首先，VTI弹性介质的弹性波伴随波动方程虽然可以准确的实现伴随波场的数值模拟，但是直接转置波动方程算子的伴随波动方程的自由边界条件难以有效实施。主要原因是转置后的伴随波动方程的伴随变量的物理意义不明确。因此，目前参考文献中，大多学者在波场模拟中只考虑了吸收边界，不实施自由边界(Ren andLiu,2016)。其次，通过直接转置推导的VTI声波伴随状态方程存在两个问题：1、伴随方程不考虑模型参数随空间变化的影响时，伴随波动方程模拟的波场在各向异性参数的突变位置出会现模拟不稳定的问题(如图1(b))；2、伴随方程考虑模型参数随空间变化的影响时，伴随波场会存在深部振幅值较弱的问题(图1(d))。
发明内容
针对上述技术问题，本发明提供一种基于自伴随VTI方程的全波形反演方法、装置、介质及设备，该方法所推导的伴随方程不会存在伴随波场模拟不稳定和深部振幅值较弱的问题；其次，该方法的伴随波场可以准确的实现地表自由边界条件；最后，该方法可以采用同一套代码实现正传波场和伴随波场的数值模拟，从而降低编程的复杂度。
为实现上述目的，本发明采取以下技术方案：
一种基于自伴随VTI方程的全波形反演方法，包括：
对VTI声波方程和VTI弹性波方程中包含有各向异性参数ε的正应力τxx进行拆分，得到拆分的VTI声波方程和VTI弹性波方程；
构建辅助变量p和q并带入拆分的VTI声波方程和VTI弹性波方程中，得到自伴随的VTI声波方程和自伴随的VTI弹性波方程；
求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，计算得到正传波场u和预测数据dpre；
利用目标函数和预测数据dpre计算残差dres；
求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，反传残差dres计算伴随波场u*；
基于目标函数对模型参数的梯度公式，利用正传波场u和伴随波场u*计算梯度
计算步长α并通过梯度和迭代更新公式更新模型m，判断目标函数是否收敛，若目标函数收敛则输出结果，若目标函数不收敛则进行下一次迭代，直至目标函数收敛。
所述的基于自伴随VTI方程的全波形反演方法，优选地，拆分的VTI声波方程具体公式如下：
其中，τxx＝τxx1+τxx2，τxx和τzz分别表示水平和垂直方向的正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，ε和δ分别表示密度，纵波速度和各向异性参数。
所述的基于自伴随VTI方程的全波形反演方法，优选地，拆分的VTI弹性波方程具体公式如下：
其中，τxx＝τxx1+τxx2，τxx和τzz分别表示水平和垂直方向的正应力；τxz表示剪切方向的切应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，vs，ε和δ分别表示密度，纵波速度，横波速度和各向异性参数。
所述的基于自伴随VTI方程的全波形反演方法，优选地，自伴随的VTI声波方程具体公式如下：
其中，p和q表示辅助变量；τxx2表示拆分后水平方向的部分正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，ε和δ分别表示密度，纵波速度和各向异性参数；s表示震源；和分别表示水平方向的空间一阶导数算子，垂直方向的空间一阶导数算子和时间一阶导数算子。
所述的基于自伴随VTI方程的全波形反演方法，优选地，自伴随的VTI弹性波方程具体公式如下：
其中，p和q表示辅助变量；τxx2表示拆分后水平方向的部分正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，vs，ε和δ分别表示密度，纵波速度，横波速度和各向异性参数；s表示震源；和分别表示水平方向的空间一阶导数算子，垂直方向的空间一阶导数算子和时间一阶导数算子。
所述的基于自伴随VTI方程的全波形反演方法，优选地，目标函数对模型参数的梯度公式如下：
其中，“T”表示矩阵的转置；m表示模型参数；u表示正传波场，A表示正演算子；u*表示伴随波场；x表示地下某点的空间坐标。
所述的基于自伴随VTI方程的全波形反演方法，优选地，目标函数的公式如下：
其中，m表示模型参数，dpre表示通过模型m所产生的预测数据，dobs表示野外采集得到的观测数据。
本发明第二方面提供一种基于自伴随VTI方程的全波形反演装置，包括：
第一处理单元，用于对VTI声波方程和VTI弹性波方程中包含有各向异性参数ε的正应力τxx进行拆分，得到拆分的VTI声波方程和VTI弹性波方程；
第二处理单元，用于构建辅助变量p和q并带入拆分的VTI声波方程和VTI弹性波方程中，得到自伴随的VTI声波方程和自伴随的VTI弹性波方程；
第三处理单元，用于求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，计算得到正传波场u和预测数据dpre；
第四处理单元，用于利用目标函数和预测数据dpre计算残差dres；
第五处理单元，用于求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，反传残差dres计算伴随波场u*；
第六处理单元，用于基于目标函数对模型参数的梯度公式，利用正传波场u和伴随波场u*计算梯度
第七处理单元，用于计算步长α并通过梯度和迭代更新公式更新模型m，判断目标函数是否收敛，若目标函数收敛则输出结果，若目标函数不收敛则进行下一次迭代，直至目标函数收敛。
本发明第三方面提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任意一项所述基于自伴随VTI方程的全波形反演方法的步骤。
本发明第四方面提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述任意一项所述基于自伴随VTI方程的全波形反演方法的步骤。
本发明由于采取以上技术方案，其具有以下优点：
1、本发明在基于自伴随VTI方程的全波形反演中，通过先对正应力τxx进行拆分，然后构建辅助变量p和q来实现VTI方程的自伴随，该技术可以降低VTI介质全波形反演的编程复杂度，解决现有技术中伴随波场模拟不稳定和深部振幅值较弱的问题，并且正传波场和伴随波场可以采用同样的自由边界实施方案。
2、本发明通过对VTI声波方程的正应力τxx进行拆分和构建辅助变量p和q实现VTI声波方程的自伴随，基于自伴随VTI声波方程的全波形反演可以有效的避免伴随波场的模拟不稳定和深部振幅值较弱的问题。
3、本发明通过对VTI弹性波方程的正应力τxx进行拆分和构建辅助变量p和q实现VTI弹性波方程的自伴随，基于自伴随VTI弹性波方程的全波形反演可以有效的避免伴随波场无法有效实施自由边界的问题。
附图说明
图1为本发明一实施例提供的Marmousi2模型在1700毫秒所对应的正传和伴随波场；其中，图(a)和(c)为VTI声波方程(公式(2))所计算的正传波场；图(b)为伴随方程不考虑模型参数随空间变化(公式(8))所计算的伴随波场；图(d)为伴随方程考虑模型参数随空间变化(公式(9))所计算的伴随波场；图(e)为本发明所提出的自伴随VTI声波方程(公式(15))所计算的正传波场，图(f)为本发明所提出的自伴随VTI声波方程(公式(15))所计算的伴随波场。震源位置为6.5km，深度为10m，地表为自由边界；
图2为Marmousi2的真实模型和初始模型；其中，图(a)，(c)，(e)，(g)和(i)分别为真实的纵波速度模型，横波速度模型，密度模型和各向异性参数模型，图(b)，(d)，(f)，(h)和(j)分别为对应的初始模型；
图3为Marmousi2模型的VTI声波全波形反演结果，其中图(a)是现有技术所反演的纵波速度；图(b)本发明所反演的纵波速度；
图4为Marmousi2模型的VTI弹性波全波形反演结果，其中，图(a)和(c)分别是现有技术所反演的纵波速度和横波速度；图(b)和(d)分别是本发明所反演的纵波速度和横波速度；
图5为基于自伴随VTI方程的全波形反演方程的流程图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
除非另外定义，本发明使用的技术术语或者科学术语应当为本发明所属领域内具有一般技能的人士所理解的通常意义。本发明中使用的“第一”、“第二”、“第三”、“第四”以及类似的词语并不表示任何顺序、数量或者重要性，而只是用来区分不同的组成部分。“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同，而不排除其他元件或者物件。“连接”或者“相连”等类似的词语并非限定于物理的或者机械的连接，而是可以包括电性的连接，不管是直接的还是间接的。
为了便于描述，可以在文中使用空间相对关系术语来描述如图中示出的一个元件或者特征相对于另一元件或者特征的关系，这些相对关系术语例如为“内部”、“外部”、“内侧”、“外侧”、“下面”、“上面”等。这种空间相对关系术语意于包括除图中描绘的方位之外的在使用或者操作中装置的不同方位。
现有的全波形反演通常采用最小二乘目标函数拟合地表的观测记录和模拟的预测记录来反演地下的介质参数。它的目标函数可以表示为：
其中，m表示模型参数，dpre表示通过模型m所产生的预测数据，dobs表示野外采集得到的观测数据。
对于VTI介质的声波全波形反演，目标函数(1)受约束于VTI声波方程：
其中，τxx和τzz分别表示水平和垂直方向的正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，ε和δ分别表示密度，纵波速度和各向异性参数；s表示震源，并以纯纵波源的形式进行加载。上式(2)给出到了二维VTI声波方程，另外存在三维VTI声波方程。两者形式相似，存在沿y轴的变量τyy和vy。这里就不再赘述其具体形式。
对于VTI介质的弹性波全波形反演，目标函数(1)受约束于VTI弹性波方程：
其中，ρ，vp，vs，ε和δ分别表示密度，纵波速度，横波速度和各向异性参数；f为一个新定义的变量，可以通过vp和vs计算得到；τxx和τzz分别表示水平和垂直方向的正应力；τxz表示剪切方向的切应力；vx和vz分别表示水平和垂直方向的质点振动速度；s表示震源，并以纯纵波源的形式进行加载。同声波方程，上式(3)给出到了二维VTI弹性波动方程，另外存在三维VTI弹性波动方程。两者形式相似，存在与y轴相关的变量τyy、τxy、τyz和vy。这里就不再赘述其具体形式。
公式(2)和(3)都可以用矩阵的形式表示为：
Au(x,t)＝s(xs,t) (4)
其中，u表示正传波场，其包括多个分量，具体定义在公式(2)和(3)中，x表示地下某点的空间坐标；s表示震源，xs表示震源的空间坐标，t表示地震波传播的时间。A表示正演算子，对于公式(2)所示的VTI声波方程，它的正演算子A可以具体表示为：
对于公式(3)所示的VTI弹性波方程，它的正演算子A可以具体表示为：
VTI介质的全波形反演中正传波场可以通过公式(4)进行数值模拟。伴随波场需要所对应的伴随方程进行数值模拟。通过伴随状态法进行推导，公式(4)对应的伴随方程可以用矩阵的形式表示为：
ATu*(x,t)＝dres(xr,t) (7)
其中，u*表示伴随波场，包含有和四个分量，AT表示伴随波场的伴随算子，其是公式(5)和(6)的转置；dres表示数据残差作为伴随源，xr表示检波器的空间坐标。
通过公式(7)可以得到不考虑模型参数随空间变化的情况下，公式(2)的伴随方程可以表示为：
其中，和表示伴随波场的四个分量；δp表示压强的残差数据作为伴随源。同样，考虑模型参数随空间变化的情况下，公式(2)的伴随方程可以表示为：
公式(8)和(9)的主要差别在于伴随方程的求解过程中是否对模型参数求取对应的空间导数。
通过公式(1)、(4)和(7)，目标函数对模型参数的梯度可以表示为：
其中，“T”表示矩阵的转置；m表示模型参数；u表示正传波场，A表示正演算子；u*表示伴随波场；x表示地下某点的空间坐标。
通过公式(10)所示的梯度和局部最优化算法就可实现对初始模型参数的迭代更新，从而实现地下模型参数的反演，即求出地下模型参数，包括纵横波速度和各向异性参数的具体值。具体的迭代更新公式如下所示：
式中，mi表示第i次迭代的模型参数，mi+1表示更新后的梯度。
通过公式(5)和(6)所示的正演算子可以看出，算子A不是自对称的，正演算子A和其伴随算子AT不相等，因此在全波形反演的实施过程中需要不同的代码来分别实现正传波场和伴随波场的数值模拟。同时，基于VTI声波方程(公式2)的全波形反演所推导的伴随方程存在两个问题：1、当A不考虑模型参数随空间变化时，伴随方程在各向异性参数存在突变的位置会模拟不稳定，如图1(b)所示；2、当AT考虑模型参数随空间变化时，伴随波场会出现深部振幅值较弱的问题，如图1(d)。此外，基于VTI弹性波方程(公式3)的全波形反演所推导的伴随方程自由边界难以准确地实施。
针对现有VTI介质的全波形反演的上述缺点，本发明提出了一种基于自伴随VTI方程的全波形反演方法，该方法可以在自由地表的情况下准确的实现VTI介质的全波形反演中正传波场和伴随波场的稳定正确的数值模拟(图1(e)和1(f))。同时，由于正演算子的自伴随，即自对称，因此正演算子和伴随算子一样，编程中的正传波场和伴随波场的数值模拟只需要一套代码进行执行，可以降低代码实现的复杂度。
如图1所示，本发明所提供的基于自伴随VTI方程的全波形反演方法，包括如下具体步骤：
在本发明中，为了实现正演算子A的自伴随，首先对公式(2)和(3)中包含有各向异性参数ε的正应力τxx进行了拆分。公式(2)所示VTI声波方程拆分后可以表示为：
其中，τxx＝τxx1+τxx2，τxx和τzz分别表示水平和垂直方向的正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，ε和δ分别表示密度，纵波速度和各向异性参数。
公式(3)所示的VTI弹性波方程拆分后可以表示为：
其中，τxx＝τxx1+τxx2，τxx和τzz分别表示水平和垂直方向的正应力；τxz表示剪切方向的切应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，vs，ε和δ分别表示密度，纵波速度，横波速度和各向异性参数。
然后通过构建辅助变量p和q来使得正演算子A自伴随：
和
将公式(14)和(15)带入到公式(12)可以得到自伴随的VTI声波方程：
其中，p和q表示辅助变量；τxx2表示拆分后水平方向的部分正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，ε和δ分别表示密度，纵波速度和各向异性参数；s表示震源；和分别表示水平方向的空间一阶导数算子，垂直方向的空间一阶导数算子和时间一阶导数算子。
将公式(14)和(15)带入到公式(13)可以得到自伴随的VTI弹性波方程：
其中，p和q表示辅助变量；τxx2表示拆分后水平方向的部分正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，vs，ε和δ分别表示密度，纵波速度，横波速度和各向异性参数；s表示震源；和分别表示水平方向的空间一阶导数算子，垂直方向的空间一阶导数算子和时间一阶导数算子。
本发明还包括如下步骤：
通过解自伴随VTI声波(公式(16))或弹性波(公式(17))计算正传波场u和预测数据dpre；
通过目标函数(公式(1))和预测数据dpre计算残差dres；
通过解自伴随VTI声波(公式(16))和弹性波(公式(17))反传残差dres计算伴随波场u*
通过公式(10)利用正传波场u和伴随波场u*计算梯度
计算步长α并通过梯度和公式(11)更新模型m，判断目标函数是否收敛，若目标函数收敛则输出结果，若目标函数不收敛则进行下一次迭代，直至目标函数收敛。
公式(16)和(17)所示的自伴随VTI声波和弹性波方程，它们的正演算子A等于其转置AT。因此，基于公式(16)和公式(17)的全波形反演具有相同的正演算子A和伴随算子AT，在全波形反演的实施过程中可以采用同一套代码实现正传波场和伴随波场的数值模拟，并且伴随波场可以和正传波场采用同样的自由边界实施方案。
下面通过具体实例详细说明本发明能有效的实现VTI介质的声波和弹性波全波形反演。
本发明所采用的例子为Marmousi2模型，真实的纵波速度模型如图2(a)所示，真实的横波速度模型如图2(c)所示，真实的密度模型如图2(e)所示，真实的各向异性参数ε模型如图2(g)所示，真实的各向异性参数δ如图2(i)所示。全波形反演所需要的初始速度模型分别如图1(b)，1(d)和1(f)所示。
本发明在VTI介质的声波全波形反演的实施过程中只展示图2(b)所示的初始纵波速度模型进行更新，密度采用Gardner公式进行生成，各向异性参数采用初始模型在反演的过程中不进行更新。最终通过现有技术所得到的纵波速度反演结果如图3(a)所示，本发明的纵波速度反演结果如图3(b)所示。通过对比可以看出，现有技术由于存在伴随波场深部振幅值较弱的问题，最终反演的结果深部更新效果不如本发明的结果，而本发明的反演结果也更接近真实的纵波速度模型。
本发明在VTI介质的弹性波全波形反演的实施过程中对图2(b)和图2(d)所示的初始纵波速度和横波速度模型进行更新，密度采用Gardner公式进行生成，各向异性参数采用初始模型在反演的过程中不进行更新。最终通过现有技术所得到的纵波速度和横波速度反演结果分别如图4(a)和4(c)所示，本发明所得到的纵波速度和横波速度反演结果如图4(b)和4(d)所示。通过图4(a)和4(c)可以对比看出，现有技术的伴随波场不能有效的实施自由边界，最后的反演结果存在大量的噪音。本发明的反演结果图4(b)和4(d)因为其伴随波场可以准确的实施自由边界，因此最终的反演结果也比现有技术更接近真实速度模型。
本发明在基于自伴随VTI方程的全波形反演中，通过先对正应力τxx进行拆分，然后构建辅助变量p和q来实现VTI方程的自伴随，该技术可以降低VTI介质全波形反演的编程复杂度，解决现有技术中伴随波场模拟不稳定和深部振幅值较弱的问题，并且正传波场和伴随波场可以采用同样的自由边界实施方案。
本发明方法具有如下有益效果：首先，该方法所推导的伴随方程不会存在伴随波场模拟不稳定和深部振幅值较弱的问题；其次，该方法的伴随波场可以准确的实现自由边界；最后，该方法可以采用同一套代码实现正传波场和伴随波场的数值模拟，从而降低编程的复杂度。
在本发明的技术方案中，全波形反演的目标函数是基于最小二乘拟合的目标函数。该目标函数也可以是其他已经公开发表的目标函数，比如最优传输目标函数，归一化目标函数等。但这些其它目标函数仍然需要结合本发明中的自伴随方程才能实现稳定、高效的反演。此外，本发明中为了简化推导只写了二维情况下自伴随的VTI声波方程和弹性波方程的推导，针对三维情况也同样可以通过上述流程实现自伴随。
本发明第二方面提供一种基于自伴随VTI方程的全波形反演装置，包括：
第一处理单元，用于对VTI声波方程和VTI弹性波方程中包含有各向异性参数ε的正应力τxx进行拆分，得到拆分的VTI声波方程和VTI弹性波方程；
第二处理单元，用于构建辅助变量p和q并带入拆分的VTI声波方程和VTI弹性波方程中，得到自伴随的VTI声波方程和自伴随的VTI弹性波方程；
第三处理单元，用于求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，计算得到正传波场u和预测数据dpre；
第四处理单元，用于利用目标函数和预测数据dpre计算残差dres；
第五处理单元，用于求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，反传残差dres计算伴随波场u*；
第六处理单元，用于基于目标函数对模型参数的梯度公式，利用正传波场u和伴随波场u*计算梯度
第七处理单元，用于计算步长α并通过梯度和迭代更新公式更新模型m，判断目标函数是否收敛，若目标函数收敛则输出结果，若目标函数不收敛则进行下一次迭代，直至目标函数收敛。
本发明第三方面提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任意一项所述基于自伴随VTI方程的全波形反演方法的步骤。
本发明第四方面提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述任意一项所述基于自伴随VTI方程的全波形反演方法的步骤。
本发明是根据具体实施方式的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解为可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (10)

1.一种基于自伴随VTI方程的全波形反演方法，其特征在于，包括：
对VTI声波方程和VTI弹性波方程中的正应力τxx进行拆分为τxx1和τxx2，得到拆分的VTI声波方程和VTI弹性波方程；
构建辅助变量p和q并代入拆分的VTI声波方程和VTI弹性波方程中，得到自伴随的VTI声波方程和自伴随的VTI弹性波方程；
求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，计算得到正传波场u和预测数据dpre；
利用目标函数和预测数据dpre计算残差dres；
求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，反传残差dres计算伴随波场u*；
基于目标函数对模型参数的梯度公式，利用正传波场u和伴随波场u*计算梯度计算步长α并通过梯度和迭代更新公式更新模型m，判断目标函数是否收敛，若目标函数收敛则输出结果，若目标函数不收敛则进行下一次迭代，直至目标函数收敛。
2.根据权利要求1所述的基于自伴随VTI方程的全波形反演方法，其特征在于，拆分的VTI声波方程具体公式如下：
其中，τxx＝τxx1+τxx2，τxx和τzz分别表示水平和垂直方向的正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，ε和δ分别表示密度，纵波速度和各向异性参数。
3.根据权利要求2所述的基于自伴随VTI方程的全波形反演方法，其特征在于，拆分的VTI弹性波方程具体公式如下：
其中，τxx＝τxx1+τxx2，τxx和τzz分别表示水平和垂直方向的正应力；τxz表示剪切方向的切应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，vs，ε和δ分别表示密度，纵波速度，横波速度和各向异性参数。
4.根据权利要求3所述的基于自伴随VTI方程的全波形反演方法，其特征在于，自伴随的VTI声波方程具体公式如下：
其中，p和q表示辅助变量；τxx2表示拆分后水平方向的正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，ε和δ分别表示密度，纵波速度和各向异性参数；s表示震源；和分别表示水平方向的空间一阶导数算子，垂直方向的空间一阶导数算子和时间一阶导数算子。
5.根据权利要求4所述的基于自伴随VTI方程的全波形反演方法，其特征在于，自伴随的VTI弹性波方程具体公式如下：
其中，p和q表示辅助变量；τxx2表示拆分后水平方向的正应力；vx和vz分别表示水平和垂直方向的质点振动速度；ρ，vp，vs，ε和δ分别表示密度，纵波速度，横波速度和各向异性参数；s表示震源；和分别表示水平方向的空间一阶导数算子，垂直方向的空间一阶导数算子和时间一阶导数算子。
6.根据权利要求5所述的基于自伴随VTI方程的全波形反演方法，其特征在于，目标函数对模型参数的梯度公式如下：
其中，“T”表示矩阵的转置；m表示模型参数；u表示正传波场，A表示正演算子；u*表示伴随波场；x表示地下某点的空间坐标。
7.根据权利要求6所述的基于自伴随VTI方程的全波形反演方法，其特征在于，目标函数的公式如下：
其中，m表示模型参数，dpre表示通过模型m所产生的预测数据，dobs表示野外采集得到的观测数据。
8.一种基于自伴随VTI方程的全波形反演装置，其特征在于，包括：
第一处理单元，用于对VTI声波方程和VTI弹性波方程中的正应力τxx进行拆分为τxx1和τxx2，得到拆分的VTI声波方程和VTI弹性波方程；
第二处理单元，用于构建辅助变量p和q并代入拆分的VTI声波方程和VTI弹性波方程中，得到自伴随的VTI声波方程和自伴随的VTI弹性波方程；
第三处理单元，用于求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，计算得到正传波场u和预测数据dpre；
第四处理单元，用于利用目标函数和预测数据dpre计算残差dres；
第五处理单元，用于求解自伴随的VTI声波方程和自伴随的VTI弹性波方程，反传残差dres计算伴随波场u*；
第六处理单元，用于基于目标函数对模型参数的梯度公式，利用正传波场u和伴随波场u*计算梯度
第七处理单元，用于计算步长α并通过梯度和迭代更新公式更新模型m，判断目标函数是否收敛，若目标函数收敛则输出结果，若目标函数不收敛则进行下一次迭代，直至目标函数收敛。
9.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1-7任意一项所述基于自伴随VTI方程的全波形反演方法的步骤。
10.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1-7任意一项所述基于自伴随VTI方程的全波形反演方法的步骤。
CN202411218000.9A 2024-09-02 2024-09-02 基于自伴随vti方程的全波形反演方法、装置、介质及设备 Active CN119224832B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202411218000.9A CN119224832B (zh) 2024-09-02 2024-09-02 基于自伴随vti方程的全波形反演方法、装置、介质及设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202411218000.9A CN119224832B (zh) 2024-09-02 2024-09-02 基于自伴随vti方程的全波形反演方法、装置、介质及设备

Publications (2)

Publication Number Publication Date
CN119224832A CN119224832A (zh) 2024-12-31
CN119224832B true CN119224832B (zh) 2025-10-03

Family

ID=93940815

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202411218000.9A Active CN119224832B (zh) 2024-09-02 2024-09-02 基于自伴随vti方程的全波形反演方法、装置、介质及设备

Country Status (1)

Country Link
CN (1) CN119224832B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116068621A (zh) * 2021-11-01 2023-05-05 中国石油天然气股份有限公司 一种基于刚度矩阵分解的各向异性介质正演方法及系统
CN116933586A (zh) * 2023-07-13 2023-10-24 中海石油(中国)有限公司 一种TI介质qP波有限差分数值模拟方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015104059A1 (en) * 2014-01-10 2015-07-16 Statoil Petroleum As Determining a component of a wave field
CN118191921A (zh) * 2024-03-20 2024-06-14 中国石油大学(北京) 一种基于声波能量的全波形反演方法和系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116068621A (zh) * 2021-11-01 2023-05-05 中国石油天然气股份有限公司 一种基于刚度矩阵分解的各向异性介质正演方法及系统
CN116933586A (zh) * 2023-07-13 2023-10-24 中海石油(中国)有限公司 一种TI介质qP波有限差分数值模拟方法

Also Published As

Publication number Publication date
CN119224832A (zh) 2024-12-31

Similar Documents

Publication Publication Date Title
CN103091711B (zh) 基于时间域一阶弹性波动方程的全波形反演方法及装置
CN110058307B (zh) 一种基于快速拟牛顿法的全波形反演方法
KR20130060231A (ko) 지구 물리학적 데이터의 반복 반전의 아티팩트 감소
CN112925012A (zh) 地震全波形反演方法及装置
CN113341465A (zh) 方位各向异性介质的地应力预测方法、装置、介质及设备
WO2017162731A1 (en) Method of operating a data-processing system for the simulation of the acoustic wave propagation in the transversely isotropic media comprising an hydrocarbon reservoir
CN115932960B (zh) 全矩张量震源机制反演方法、系统、设备及存储介质
WO2022153984A1 (ja) 学習データ生成方法、モデル生成方法および学習データ生成装置
CN106569262A (zh) 低频地震数据缺失下的背景速度模型重构方法
US11199641B2 (en) Seismic modeling
CN104597489B (zh) 一种震源子波优化设置方法和装置
CN110376642B (zh) 一种基于锥面波的三维地震速度反演方法
CN119224832B (zh) 基于自伴随vti方程的全波形反演方法、装置、介质及设备
Li et al. Time-domain wavefield reconstruction inversion
CN111175822A (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
Fowler et al. Generalized pseudospectral methods for orthorhombic modeling and reverse-time migration
CN117934748B (zh) 基于深度三视图的重力异常智能反演方法
CN118191921A (zh) 一种基于声波能量的全波形反演方法和系统
WO2011041440A2 (en) Correcting an acoustic simulation for elastic effects
CN117933319A (zh) 深部随机混合介质非平稳模型建模方法、系统及电子设备
CN116933586A (zh) 一种TI介质qP波有限差分数值模拟方法
CN116660981A (zh) 基于包络的各向异性参数反演方法、装置及存储介质
CN114089419B (zh) 优化的变网格地震正演模拟方法
CN120762097B (zh) 弹性波全波形反演方法、系统、装置及介质
CN113589362B (zh) 三维陆上耦合波正演模拟方法

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