CN113740901B - 基于复杂起伏地表的陆上地震数据全波形反演方法及装置 - Google Patents

基于复杂起伏地表的陆上地震数据全波形反演方法及装置 Download PDF

Info

Publication number
CN113740901B
CN113740901B CN202010478692.6A CN202010478692A CN113740901B CN 113740901 B CN113740901 B CN 113740901B CN 202010478692 A CN202010478692 A CN 202010478692A CN 113740901 B CN113740901 B CN 113740901B
Authority
CN
China
Prior art keywords
seismic
data
work area
model
target work
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
CN202010478692.6A
Other languages
English (en)
Other versions
CN113740901A (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.)
Petrochina Co Ltd
Original Assignee
Petrochina Co Ltd
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 Petrochina Co Ltd filed Critical Petrochina Co Ltd
Priority to CN202010478692.6A priority Critical patent/CN113740901B/zh
Publication of CN113740901A publication Critical patent/CN113740901A/zh
Application granted granted Critical
Publication of CN113740901B publication Critical patent/CN113740901B/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
    • 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/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616—Data from specific type of measurement
    • G01V2210/6161—Seismic or acoustic, e.g. land or sea measurements
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/66—Subsurface modeling
    • G01V2210/665—Subsurface modeling using geostatistical modeling
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/67—Wave propagation modeling
    • G01V2210/673—Finite-element; Finite-difference

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)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种基于复杂起伏地表的陆上地震数据全波形反演方法及装置，该方法包括：对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；根据所述多尺度反向传播波场更新所述速度模型。本发明实施例可以提供一种适用于复杂起伏地表的陆上地震资料全波形反演方法，可提高陆上双复杂探区深度域速度建模精度，从而改善成像质量。

Description

基于复杂起伏地表的陆上地震数据全波形反演方法及装置
技术领域
本发明涉及地质勘探与油气勘探技术领域，特别是关于石油地球物理勘探地震数据处理技术领域，具体涉及一种基于复杂起伏地表的陆上地震数据全波形反演方法及装置。
背景技术
我国陆上中西部地区具有广阔的山前推覆构造带，成藏条件优越，赋存大量油气资源，整体探明程度低，是油气勘探突破的重点领域。山前带勘探区通常位于盆山过渡带，具有复杂地表、复杂地下构造的双复杂地质结构特征，地震采集难度大，各类干扰波广泛发育，地震数据品质低，地震成像精度低，提高双复杂探区地震成像精度是当前地震成像处理技术攻关与创新研究的主要方向。实现复杂构造成像的核心方法是叠前深度偏移技术。叠前深度偏移成像的准确性极大的依赖于深度域速度建模的可靠性。但是，双复杂探区存在地形起伏剧烈，地表出露岩性变化快，近地表厚度、速度横向变化大，地下构造破碎，断层发育，地层倾角大，产状变化多，速度横向变化剧烈，局部垂向速度反转等复杂地震地质特征，现有速度建模技术对复杂地质条件适应性差，难以为双复杂探区构建可靠的深度域速度模型，亟待发展适应陆上双复杂构造特征的速度反演建模新技术以提高深度域速度建模的可靠性和成像质量。
发明内容
针对现有技术中的问题，本发明提供的基于复杂起伏地表的陆上地震数据全波形反演方法及装置，能够有效适应复杂起伏的地表条件，有效反演复杂变化的速度结构特征。
为解决上述技术问题，本发明提供以下技术方案：
第一方面，本发明提供一种基于复杂起伏地表的陆上地震数据全波形反演方法，包括：
对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；
根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；
根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；
根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；
根据所述多尺度反向传播波场更新所述速度模型。
一实施例中，所述对目标工区进行非规则网格离散，以生成所述目标工区的几何模型包括：
解析所述目的工区的地震数据，以获取所述地震数据的局部坐标采集观测***；所述局部坐标采集观测***包括：震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程；
根据所述震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程，计算所述目标工区的地表坐标；
利用非规则网格离散方法，根据所述地表坐标以及工区范围生成所述目标工区的几何模型。
一实施例中，所述根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据，包括：
根据所述地震数据求取地震子波；
利用阻尼间断伽辽金数值模拟方法，根据所述地震子波、所述速度模型、所述目标工区非规则网格下的密度模型、非规则网格节点、有限元数据以及所述地震子波对所述目标工区进行正演模拟，以生成多尺度正演模拟数据。
一实施例中，所述根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场，包括：
利用所述阻尼间断伽辽金数值模拟方法以及频率域地震波数学物理方程，根据速度模型以及伴随震源场求解所述多尺度反向传播波场。
一实施例中，所述根据所述多尺度反向传播波场更新所述速度模型包括：
利用伴随状态法，根据所述多尺度反向传播波场以及多尺度正演模拟数据，生成所述速度模型的梯度以及对角线海森矩阵；
利用共轭梯度法优化方法以及递进尺度方式，根据所述梯度以及对角线海森矩阵更新所述速度模型。
第二方面，本发明提供一种基于复杂起伏地表的陆上地震数据全波形反演装置，该装置包括：
几何模型生成单元，用于对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；
模拟数据生成单元，用于根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；
震源场生成单元，用于根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；
传播波场生成单元，用于根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；
模型更新单元，用于根据所述多尺度反向传播波场更新所述速度模型。
一实施例中，所述几何模型生成单元包括：
地震数据解析模块，用于解析所述目的工区的地震数据，以获取所述地震数据的局部坐标采集观测***；所述局部坐标采集观测***包括：震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程；
地表坐标计算模块，用于根据所述震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程，计算所述目标工区的地表坐标；
几何模型生成模块，用于利用非规则网格离散方法，根据所述地表坐标以及工区范围生成所述目标工区的几何模型。
一实施例中，所述模拟数据生成单元包括：
子波求取模块，用于根据所述地震数据求取地震子波；
模拟数据生成模块，用于利用阻尼间断伽辽金数值模拟方法，根据所述地震子波、所述速度模型、所述目标工区非规则网格下的密度模型、非规则网格节点、有限元数据以及所述地震子波对所述目标工区进行正演模拟，以生成多尺度正演模拟数据。
一实施例中，所述传播波场生成单元具体用于利用所述阻尼间断伽辽金数值模拟方法以及频率域地震波数学物理方程，根据速度模型以及伴随震源场求解所述多尺度反向传播波场。
一实施例中，所述模型更新单元包括：
梯度生成模块，用于利用伴随状态法，根据所述多尺度反向传播波场以及多尺度正演模拟数据，生成所述速度模型的梯度以及对角线海森矩阵；
模型更新模块，用于利用共轭梯度法优化方法以及递进尺度方式，根据所述梯度以及对角线海森矩阵更新所述速度模型。
第三方面，本发明提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行程序时实现基于复杂起伏地表的陆上地震数据全波形反演方法的步骤。
第四方面，本发明提供一种计算机可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现基于复杂起伏地表的陆上地震数据全波形反演方法的步骤。
从上述描述可知，本发明实施例提供基于复杂起伏地表的陆上地震数据全波形反演方法及装置，首先对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；本发明实施例可以提供一种适用于复杂起伏地表的陆上地震资料全波形反演方法，可提高陆上双复杂探区深度域速度建模精度，从而改善成像质量。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明的实施例中的基于复杂起伏地表的陆上地震数据全波形反演方法的流程示意图；
图2为本发明的实施例中步骤100的流程示意图；
图3为本发明的实施例中步骤200的流程示意图；
图4为本发明的实施例中步骤400的流程示意图；
图5为本发明的实施例中步骤500的流程示意图；
图6为本发明的具体应用实例中基于复杂起伏地表的陆上地震数据全波形反演方法的流程示意图；
图7为本发明的具体应用实例中基于复杂起伏地表的陆上地震数据全波形反演方法的思维导图；
图8为本发明的具体应用实例中包含起伏地表的几何模型示意图
图9为本发明的具体应用实例中几何模型对应的速度场；
图10为本发明的具体应用实例中几何模型对应的密度场；
图11为本发明的具体应用实例中几何模型的剖分示意图(不依赖背景场)；
图12为本发明的具体应用实例中几何模型的剖分示意图(依赖背景场)；
图13为本发明的具体应用实例中Foothill速度模型示意图；
图14为本发明的具体应用实例中用于起伏地表波形层析的初始速度模型示意图；
图15为本发明的具体应用实例中Foothill几何模型示意图；
图16为本发明的具体应用实例中Foothill网格数据示意图；
图17为本发明的具体应用实例中地震记录数据示意图(炮号：98-103)；
图18为本发明的具体应用实例中用于起伏地表波形层析反演所得的速度模型示意图(第一个尺度(1Hz～6Hz)第10次迭代结果)；
图19为本发明的具体应用实例中用于起伏地表波形层析反演所得的速度模型示意图(第二个尺度(6Hz～12Hz)第4次迭代结果)；
图20为本发明的实施例中基于复杂起伏地表的陆上地震数据全波形反演装置的结构示意图；
图21为本发明的实施例中几何模型生成单元的结构框图；
图22为本发明的实施例中模拟数据生成单元的结构框图；
图23为本发明的实施例中模型更新单元的结构框图；
图24为本发明的实施例中的电子设备的结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
地球物理方法是资源勘探的重要工具，包括重力、磁场、电场、地热及地震等方法，其中，地震勘探方法是石油勘探应用最广泛的方法，该方法包括地震资料采集、处理和解释三大步骤。根据勘探环境的不同可以分为陆上地震勘探和海上地震勘探，它们涉及的主要问题差别明显，尤其是地震资料处理领域，陆上资料受复杂地表的影响，资料信噪比差，难以得到复杂的近地表速度模型。
全波形反演方法是目前地震资料处理解释领域最先进的基础理论方法，传统的地震资料处理流程以都可归于全波形反演理论框架内。全波形反演有两部分的主要内容：地震波正演数值模拟方法和反演优化算法。复杂地表的问题既涉及地表起伏又涉及到复杂变化的近地表速度，因此传统有限差分正演模拟数值方法不能很好地解决正演问题，因而反演问题也无法有效地解决。本发明所提适用于复杂起伏地表的波形匹配全波形反演方法,使用非规则网格表示复杂起伏地表模型和改进的有限元数值模拟方法高效地实现地震波数学物理方程的数值模拟，以及使用多尺度波形匹配的反演方法实现速度模型更新。
陆上地震勘探的难点是众所周知的“双复杂”问题，而复杂地表问题更为突出，不解决复杂地表问题去考虑解决复杂地下问题也就无从谈起。复杂地表问题的表现有两方面：一是“起伏”问题，即地表高程变化大，西部山区变化还特别剧烈；二是“变速”问题，表层岩石塑性严重、高速岩体出露、横向速度变化大以及纵向速度差异明显。这些近地表问题直接导致地震波传播复杂，地震资料信噪比差，给后续的处理流程包括地震速度分析和成像带了巨大的困难。解决近地表问题，常规的处理流程中主要考虑地表起伏对旅行时差的影响而进行静态时移校正，最多使用的是基于地表一致性假设条件的各类静校正方法；进一步的处理方法是波动方程或者克西霍夫积分法基准面校正(Berryhill,1984和Wiggins,1984)，这类方法需要人为确定填充速度，有时不合理的填充速度反而给后续的成像处理引入新的旅行时误差；直接在起伏地表进行成像也是解决近地表问题的途径之一，但前提需要准确的近地表模型。实际上，准确的近地表速度模型，不但是起伏地表成像方法的依赖条件，也能够提高常规处理流程步骤的精度。因此，如何建立准确的近地表模型是中国陆上勘探迫切需要解决的难题之一。近地表建模的基础是表层速度分析，这些都有一些经典方案，包括小折射、微测井、面波反演、初至波反演等方法；近地表建模除了表层的速度还包括浅层速度模型，广泛应用的方法仍然是基于初至波旅行时或者基于波形信息的层析技术或者波形反演方法更新速度场。
建模问题可以归类为“反问题”。在石油勘探领域，地震勘探方法使用人工震源产生的地震波沿地球介质传播，经过波动传播理论的透射、反射、折射等传播机制，把地球分层、岩性等信息携带到地面并通过地面传感器转化成数字信号，最后信号处理工程师、地球物理学家和地质学家对这些信号进行一系列处理和解释得到地球地下介质可视化图像和可识别的岩性资料。这种识别地下构造和岩性的地震勘探过程在数学上可以抽象归类为“反问题”，地震勘探最终的目标是建立地下介质岩性和物性的“可视化”地震像，地震建模可以看成是整个勘探反问题的一个中间环节，但其本身就是典型的“反问题”。求解该反问题的过程首先是通过波动力学的数学物理方程建立波场和速度模型正确的映射关系模型，然后选择合适的线性或者反线性，全局或者局部优化方案求解反演问题。在20世纪80年代首次提出的全波形反演方法(Lailly，1983；Tarantola，1984)与传统地震资料处理建立反问题的流程相比，就是对地震勘探反问题最直接的和最全面的定义。然而该反演问题的非线性、初始模型偏差过大、采集数据缺乏大偏移距和低频信息、实际地震子波的复杂性、噪音干扰以及数据物理方程描述不准确等众多因素的影响，全波形反演方法在实际应用中很难得到期望的高分辨率的岩石物理参数模型。面对低波数信息的缺失，也缺少精确表达各种物理岩性参数的数学物理方程，实际应用中全波形反演仍主要作为一个速度建模方法被广泛研究，即使如此，非线性和缺少低频大偏移数据等问题仍然存在。
由于全波形反演方法在理论上仍然肩负着地球物理学家对提高实际地震介质高精度地震像的期望，上述诸多问题近几年通过地球物理学家提出的不同方法和策略逐步解决或者减缓(Luo&Schuster,1991；Pratt,1999；Ma&Hale,2013；Engquist&Froese,2014；Metivier,et al.2016；Yang et al.,2018)。总的来说，上述问题在反演求解中多数表现为优化方法适用性问题，如非线性特性和周期跳跃造成因全波形常用的最小二乘算法常常收敛缓慢或者收敛于局部极小值。解决上述问题的思路，会涉及到精确的数学物理模型(方程)建立、数学物理方程的数值模拟、反演优化方法选择和实施，一般会分为正演、反演两个问题来研究。在起伏地表情况下，速度模型复杂构造的精确表达也是需要关注的重点问题。基于上述原因，本发明的实施例提供一种基于复杂起伏地表的陆上地震数据全波形反演方法的具体实施方式，参见图1，该方法具体包括如下内容：
步骤100：对目标工区进行非规则网格离散，以生成所述目标工区的几何模型。
步骤100在实施时，具体为：基于复杂起伏地表初始模型对模型空间进行非规则网格离散，其有益效果是解决了复杂速度模型及复杂几何形态的高精度离散表达问题。
步骤200：根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据。
具体地，基于从地震记录数据求取的地震子波、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法实现地震波的多炮正演模拟，可以理解的是，步骤200实现复杂速度模型的多尺度波场高效正演模拟以及解决频率域震源周期效应。
步骤300：根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场。
具体地，基于多尺度正演模拟数据和观测地震记录数据的多尺度表示，使用匹配波场的目标函数求解函数值以及伴随震源，可以实现观测地震记录数据和正演模拟数据的有效匹配，并减少陆上勘探地震资料波形层析反演的非线性。
步骤400：根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场。
步骤400在实施时，具体为：基于伴随波场、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法求解伴随震源场的多尺度反向传播波场，可以高效地实现适应起伏地表的反向地震波场传播。
步骤500：根据所述多尺度反向传播波场更新所述速度模型。
具体地，根据多尺度正演模拟数据和多尺度反向传播波场，使用伴随状态法得到梯度和对角线海森矩阵，并进一步使用共轭梯度法优化方法和递进尺度方式更新速度场，以实现最终复杂起伏地表速度模型的有效更新。
从上述描述可知，本发明实施例提供基于复杂起伏地表的陆上地震数据全波形反演方法，首先对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；本发明实施例可以提供一种适用于复杂起伏地表的陆上地震资料全波形反演方法，可提高陆上双复杂探区深度域速度建模精度，从而改善成像质量。
本发明实施例所提供的方法是为了解决陆上地震勘探面临的复杂起伏地表速度建模问题。首先，基于复杂起伏地表初始模型对模型空间进行非规则网格离散，目的是解决复杂速度模型及复杂几何形态的高精度离散表达问题；然后，基于从地震记录数据求取的地震子波、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法实现地震波的多炮正演模拟，目的是实现复杂速度模型的多尺度波场高效正演模拟以及解决频率域震源周期效应；其次，基于多尺度正演模拟数据和观测地震记录数据的多尺度表示，使用匹配波场的目标函数求解函数值以及伴随震源，目的是实现观测地震记录数据和正演模拟数据的有效匹配，减少陆上勘探地震资料波形层析反演的非线性；再次，基于伴随波场、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法求解伴随震源场的多尺度反向传播波场，目的是高效地实现适应起伏地表的反向地震波场传播；最后，基于多尺度正演模拟数据和多尺度反向传播波场，使用伴随状态法得到梯度和对角线海森矩阵，进一步使用共轭梯度法优化方法和递进尺度方式更新速度场，目的是实现最终复杂起伏地表速度模型的有效更新。
一实施例中，参见图2，步骤100包括：
步骤101：解析所述目的工区的地震数据，以获取所述地震数据的局部坐标采集观测***。
在步骤101中，局部坐标采集观测***包括：震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程。
步骤102：根据所述震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程，计算所述目标工区的地表坐标。
步骤103：利用非规则网格离散方法，根据所述地表坐标以及工区范围生成所述目标工区的几何模型。
在步骤101至步骤103中，从地震记录数据道头信息中得到局部坐标采集观测***，其中包括地震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标、检波点高程等信息。根据采集观测***中的震源点坐标、震源点高程、震源点深度、检波点坐标、检波点高程进行二次样条插值以及外推得到数据所涉及工区的地表坐标(起伏地表坐标)，通过起伏地表坐标和工区范围组建所涉及工区的几何模型。进一步地，根据局部坐标采集观测***裁剪初始速度模型、密度模型，得到数据所涉及工区的局部初始速度模型、密度模型以及对应的坐标范围，进一步构造出几何模型对应的背景速度场和密度场，接着，依据背景速度场、密度场描述的地下复杂构造和几何模型包含的起伏地表复杂构造，对所涉及工区的几何模型使用德劳内方法进行三角网格剖分，获得包括元素、节点、空间位置等信息在内的网格数据，进一步抽取本发明需要的元素和节点的数据。节点数据主要包括节点的编号、坐标，元素数据主要包括元素的编号、构成该元素的节点数目和编号，元素和节点信息反映出剖分网格的布局和大小。
一实施例中，参见图3，步骤200包括：
步骤201：根据所述地震数据求取地震子波。
可以理解的是，在步骤201之前还需要将地震记录数据炮点坐标和检波点坐标与局部模型坐标进行匹配，筛选工区范围内的地震数据，对所选地震数据进行预处理操作主要操作包括：去除随机噪声，衰减面波等强能量的干扰波，剔除坏道，切除早至波之前的异常干扰，并进行子波整形预处理。
步骤202：利用阻尼间断伽辽金数值模拟方法，根据所述地震子波、所述速度模型、所述目标工区非规则网格下的密度模型、非规则网格节点、有限元数据以及所述地震子波对所述目标工区进行正演模拟，以生成多尺度正演模拟数据。
具体地，根据元素和节点数据以及初始速度模型和密度模型，使用线性插值可以得到非规则网格点对应的速度值和密度值，构造出非规则网格下速度模型和密度模型。接着，对预处理后的地震记录数据进行傅里叶变换，将其由时间域变换到频率域中，并根据尺度个数，尺度范围，得到不同尺度下的地震数据分量。对不同尺度的地震数据进行线性反演得到频率域多尺度子波。
从非规则网格离散的速度模型、密度模型、非规则网格节点和有限元数据以及提取所得的多尺度地震子波，使用阻尼间断伽辽金多炮数值模拟方法进行数值模拟。阻尼间断伽辽金多炮数值模拟方法与传统间断伽辽金有限元的区别在于为了解决频率域模拟的周期效应，在传统间断伽辽金有限元的基础上引入阻尼项。
一实施例中，参见图4，步骤400包括：
步骤401：利用所述阻尼间断伽辽金数值模拟方法以及频率域地震波数学物理方程，根据速度模型以及伴随震源场求解所述多尺度反向传播波场。
具体地，基于伴随波场、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法求解伴随震源场的多尺度反向传播波场。
一实施例中，参见图5，步骤500包括：
步骤501：利用伴随状态法，根据所述多尺度反向传播波场以及多尺度正演模拟数据，生成所述速度模型的梯度以及对角线海森矩阵。
步骤502：利用共轭梯度法优化方法以及递进尺度方式，根据所述梯度以及对角线海森矩阵更新所述速度模型。
在步骤501以及步骤502中，基于正向传播的波场数据和伴随波场数据，使用伴随状态方法求取目标函数相对于速度的梯度，进一步使用共轭梯度优化方法及输入参数初始步长，求取速度更新量，并更新速度值。
递进尺度方式速度场更新：根据尺度序号、对应的尺度范围和最大迭代次数，按照步骤200至步骤400更新速度值，然后进入到下一个尺度，依次递进尺度方式更新速度场直到完成所有尺度计算。
从上述描述可知，本发明实施例提供基于复杂起伏地表的陆上地震数据全波形反演方法，基于复杂起伏地表初始模型对模型空间进行非规则网格离散；接着，基于从地震记录数据求取的地震子波、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法实现地震波的多炮正演模拟；根据多尺度正演模拟数据和观测地震记录数据的多尺度表示，使用匹配波场的目标函数求解函数值以及伴随震源；并基于伴随波场、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法求解伴随震源场的多尺度反向传播波场；最后基于多尺度正演模拟数据和多尺度反向传播波场，使用伴随状态法得到梯度和对角线海森矩阵，进一步使用共轭梯度法优化方法和递进尺度方式更新速度场。
本发明的有益效果体现在：(1)在复杂起伏地表处没有网格离散阶梯效应，模拟结果没有“离散阶梯”数值噪声；(2)波形匹配目标函数求残差具有强鲁棒性特征，没有“周波跳跃”现象，能正确求取模拟结果与观测数据的残差，基于波形匹配残差的反演更新能获得正确的速度模型；(3)多尺度反演能够适应复杂起伏地表变化，能够正确更新速度细节特征；(4)基于本发明方法更新的速度模型，具有正确的速度结构和数值，能够显著提升成像质量。
为进一步地说明本方案，本发明采用工业和学术界标准的复杂起伏地表和复杂地下构造模型Foothill模型，其用来测试基于复杂起伏地表的陆上地震数据全波形反演方法，以此提供基于复杂起伏地表的陆上地震数据全波形反演方法的具体应用实例，该具体应用实例具体包括如下内容，参见图6以及图7。
复杂起伏地表是导致陆上双复杂探区深度域速度建模精度低的主要原因。本具体应用实例的目的是提供一种适用于复杂起伏地表的陆上地震资料全波形反演方法，用以提高陆上双复杂探区深度域速度建模精度，改善成像质量。
可以理解的是，全波形反演是当前精度最高且最先进的速度建模技术。全波形反演技术实施的关键环节有两个，地震波正演数值模拟求残差和反演优化速度更新。复杂地表问题涉及地表起伏和近地表速度复杂变化两方面，传统有限差分正演模拟数值方法难以刻画复杂变化的自由边界条件，不能反映复杂变化的速度特征，难以提供逼真的正演模拟结果，从而不能获取正确的数值残差，导致反演失败。传统有限差分不适用于复杂地区和复杂地质体模拟和反演。
在本具体应用实例中，复杂地下构造模型Foothill模型密度为常量，速度模型的大小：横向25005米，纵向9990米；原始规则网格速度模型横向采样间隔为15米，横向网格点数为1668个，纵向采样间隔为10米，纵向网格点数为1000个；模型的地表比较复杂，起伏陡峭，高差有1200多米，地表岩性变化较大，横向速度差异变化明显，速度从3600米/秒～4000米/秒，近地表纵向速度分层也很明显，因此近地表建模困难；地下构造复杂，断层发育，褶皱严重，其地震波传播不但受地表复杂性的影响，本身成像也具有一定难度。
S1：数据输入和参数准备。
具体地，获取地震记录数据，在本具体应用实例中，数据使用标准的SEGY格式或者包含炮号、震源点坐标、震源点高程、震源点深度、检波点坐标、检波点高程、初至旅行时时间、地震道道数、时间采样间隔和采样点数信息的数据格式，本发明中以下叙述以SEGY格式为例；地震记录数据SEGY数据格式包括三部分，分别为EBCDIC文件头、二进制文件头以及地震道，每条地震道均包含道头信息和地震道数据；所述道头信息内保存所述地震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标、检波点高程、初至旅行时时间、地震道道数、时间采样间隔dt和采样样点数nt。
获取初始速度模型，数据使用标准的SEGY格式或者包含初始速度模型中心点位置坐标、地震道道数、纵向采样间隔和采样样点数信息的数据格式。初始速度模型数据SEGY格式包括三部分，分别为EBCDIC文件头、二进制文件头以及地震道，每条地震道均包含道头信息和地震道数据，本发明中以下叙述以SEGY格式为例；所述道头信息内保存所述地震道中心点位置坐标、地震道道数、纵向采样间隔和采样样点数信息的数据格式.
获取偏移距加权系数数据，数据格式为二进制格式，可表示为向量wh。获取主要参数：阻尼系数β；尺度个数nf，尺度范围[f1,f2,...,fn+1]；共轭梯度优化方法初始步长α；偏移距个数noff，偏移距权值系数文件偏移距个数doff；最大迭代次数m；依赖背景场(速度场或者密度场)的不规则网格细化比例系数γ。
S2：基于复杂起伏地表初始速度模型对模型空间进行非规则网格剖分。
从地震记录数据SEGY道头信息中得到局部坐标采集观测***，其中包括地震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标、检波点高程等信息。根据采集观测***中的震源点坐标、震源点高程、震源点深度、检波点坐标、检波点高程进行二次样条插值以及外推得到数据所涉及工区的起伏地表坐标，通过起伏地表坐标和工区范围组建所涉及工区的几何模型，参见图8。
根据局部坐标采集观测***裁剪初始速度模型、密度模型得到数据所涉及工区的局部初始速度模型、密度模型以及对应的坐标范围，进一步构造出几何模型对应的背景速度场和密度场，参见图9以及图10。
依据背景速度场、密度场描述的地下复杂构造和几何模型包含的起伏地表复杂构造，对所涉及工区的几何模型使用德劳内方法进行三角网格剖分，获得包括元素、节点、空间位置等信息在内的网格数据，进一步抽取本发明需要的元素和节点的数据。节点数据主要包括节点的编号、坐标，元素数据主要包括元素的编号、构成该元素的节点数目和编号，元素和节点信息反映出剖分网格的布局和大小。上述三角剖分方法，本发明中称之为尺度依赖背景场(包含复杂构造信息)的德劳内三角剖分方法(如图11以及图12所示)，该方法的实现算法，是基于传统的德劳内三角剖分算法，在网格细化时网格的尺度l定义如下：
其中，γ为不规则网格细化比例系数,v(x)为空间位置x处的速度值，vmax和vmin分别时速度最大值和最小值。
S3：利用阻尼间断伽辽金多炮数值模拟方法对目标工区进行正演模拟。
步骤S3在实施时，具体地，根据元素和节点数据以及初始速度模型和密度模型使用线性插值可以得到非规则网格点对应的速度值和密度值，构造出非规则网格下速度模型和密度模型。
地震记录数据炮点坐标和检波点坐标与局部模型坐标进行匹配，筛选工区范围内的地震数据，对所选地震数据进行预处理操作主要操作包括：去除随机噪声，衰减面波等强能量的干扰波，剔除坏道，切除早至波之前的异常干扰，并进行子波整形预处理。
对预处理后的地震记录数据进行傅里叶变换，将其由时间域变换到频率域中，并根据尺度个数nf，尺度范围[f1,f2,...,fn+1]，得到不同尺度下的地震数据分量，表示为uobs。对不同尺度的地震数据进行线性反演得到频率域多尺度子波。
从非规则网格离散的速度模型、密度模型、非规则网格节点和有限元数据以及提取所得的多尺度地震子波，使用阻尼间断伽辽金多炮数值模拟方法进行数值模拟。阻尼间断伽辽金多炮数值模拟方法与传统间断伽辽金有限元的区别在于为了解决频率域模拟的周期效应，在传统间断伽辽金有限元的基础上引入阻尼项。频率域的变密度声波方程为
其中，ρ(x)为密度，κ(x)＝ρ(x)v2(x)为体变模量，v(x)为速度，t为旅行时间，x为空间位置点，xs为震源点空间位置点，▽是梯度符号，ω为角频率，u(x,ω)为单频波场，s(xs,ω)为频率域的震源函数。
本具体应用实例使用间断伽辽金有限元方法求解方程(2)，首先，用近似的方法来计算精确的解其中φk是线性独立的基函数(有限元中称为形函数或插值函数)，uk取节点处的场值。取与基函数相同的权函数φj，经推导可以得到伽辽金有限元求解的弱形式方程，如下
其中,表示计算内部区域边界的法向。方程(3)是声波方程的弱形式，通常可以表示成如下简单的线性方程形式：
Lu＝s (4)
其中，L为阻抗矩阵或者微分算子矩阵，依赖于密度、体积模量以及频率；u和s分别是波场和震源的矩阵表示。与方程(4)类比，改进的阻尼方法方程表示如下，
Lu′＝s′ (5)
阻尼方法方程(5)中的u′和s′分别为，
和
其中，β阻尼系数，t地震记录数据的纵向坐标旅行时，dt和nt分别时地震记录数据的采样间隔和样点长度。
方程(4)和(5)右侧都可以引入多个炮点记录，如此公式，根据非规则速度模型和密度模型、多尺度地震子波可以求解多炮多尺度地震合成记录并记录正向传播的波场数据，下面表示为usyn。
S4：使用匹配波场的目标函数求解函数值以及伴随震源。
根据输入的偏移距权值系数wh，多尺度地震记录数据uobs和多尺度地震合成记录usyn定义如下目标函数：
其中，wh偏移距权值系数，因此残差波场可以表示为
δu＝ln(wh·usyn)-ln(wh·uobs) (9)
目标函数的梯度：
方程(10)中，为伴随震源。
S5：求取多尺度反向传播波场。
根据多尺度地震数据记录和多尺度合成数据记录求取对数波场残差(方程9)，然后根据对数波场参数和多尺度合成数据记录构建伴随震源，使用进行伴随震源的反向传播，即使用方程(5)进行波场求解。
S6：更新速度场。
基于正向传播的波场数据和伴随波场数据，使用伴随状态方法和方程(10)求取目标函数相对于速度的梯度，进一步使用共轭梯度优化方法及输入参数初始步长α，求取速度更新量，并更新速度值。
递进尺度方式速度场更新：根据尺度序号i、对应的尺度范围[fi,fi+1]和最大迭代次数m，利用步骤S3至S5更新速度值，然后进入到下一个尺度需要将i加1(i+1)，依次递进尺度方式更新速度场直到完成所有尺度计算。
具体实施方式以及发明的效果：
地震记录数据是通过实际速度模型(图13)正演模拟获取，一共278炮，炮间距为90米，每炮240个检波点，检波距为15米，最大偏移距为3600米，时间采样间隔为4毫秒，道长5秒。图14用于起伏地表波形层析的初始速度模型是流程中需要输入的初始模型。通过地震记录数据的道头信息提取起伏地表坐标；图14中的速度模型可以通过边缘检测手段提取出更多的构造信息帮助构建更完整的Foothill几何模型(如图15)。由于密度场为常量，这里基于速度场(图14)和几何模型(图15)生成依赖于背景场的网格数据(如图16)，网格的大小与速度场分布一致。
基于实际速度模型(图13)获取的地震数据(图17，第98-103炮数据)，以及网格数据抽取的节点和有限元数据(图15)得到的非规则网格初始速度模型(图14)，并设置本发明适应起伏地表的波形层析方法相关参数：阻尼系数β＝50；尺度个数nf＝2，尺度范围[1.0Hz，6.0Hz，12.0Hz]；共轭梯度优化方法初始步长α＝0.01，偏移距权值系数为常系数1，偏移距个数noff＝23、偏移距权值系数文件偏移距个数doff＝200，最大迭代次数设置为m＝10。
经过本发明适用于复杂起伏地表的波形匹配全波形反演方法反演更新初始速度，得到多尺度反演结果图18和图19。图18是第一个尺度(1-6Hz)第10次迭代地震记录数据反演更新速度结果，图中可以观察到在初始速度(图14)上更高波数的速度模型成分被恢复；图19是第二个尺度(6-12Hz)第四次迭代地震记录数据反演更新速度结果，同样可以观察到明显高波数速度成分被恢复，与真实速度模型(图13)比较起伏地表近地表的构造已经相当清晰，得到理想的结果。
从上述描述可知，本发明实施例提供基于复杂起伏地表的陆上地震数据全波形反演方法，首先使用非规则网格表示复杂起伏地表模型和改进的有限元数值模拟方法高效地实现地震波数学物理方程的数值模拟，以及使用多尺度波形匹配的反演方法实现速度模型更新。关键步骤包括：一、基于复杂起伏地表初始模型对模型空间进行非规则网格离散；二、基于从地震记录数据求取的地震子波、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法实现地震波的多炮正演模拟；三、基于多尺度正演模拟数据和观测地震记录数据的多尺度表示，使用匹配波场的目标函数求解函数值以及伴随震源；四、基于伴随波场、频率域地震波数学物理方程及非规则网格离散的速度模型，使用阻尼间断伽辽金数值模拟方法求解伴随震源场的多尺度反向传播波场；五、基于多尺度正演模拟数据和多尺度反向传播波场，使用伴随状态法得到梯度和对角线海森矩阵，进一步使用共轭梯度法优化方法和递进尺度方式更新速度场。
相对地，现有技术特征包括：基于矩形网格或非规则网格离散速度模型，使用有限差分方法实现地震波的多炮多尺度正演模拟，使用减法目标函数求取模拟数据与观测数据的残差，使用共轭梯度法优化方法更新速度场。
现有技术效果包括：(1)在复杂起伏地表处存在网格离散阶梯效应，模拟结果存在“离散阶梯”数值噪声；(2)减法目标函数求残差导致严重的“周波跳跃”现象，不能正确求取模拟结果与观测数据的残差，基于减法残差的反演更新不能获得正确的速度模型；(3)多尺度反演无法适应复杂起伏地表变化，不能正确更新速度细节特征；(4)基于现有技术更新的速度模型，因更新数值不正确，故无法显著提高成像质量。
模型试验和实际数据测试分析表明，本方法适用于具有复杂起伏地表、复杂地下速度结构的陆上双复杂探区叠前深度域速度反演，能够有效适应复杂起伏的地表条件，能够有效反演复杂变化的速度结构特征与细节特征，有效提升陆上双复杂探区叠前深度域速度建模精度，改善成像质量。
基于同一发明构思，本申请实施例还提供了基于复杂起伏地表的陆上地震数据全波形反演装置，可以用于实现上述实施例所描述的方法，如下面的实施例。由于基于复杂起伏地表的陆上地震数据全波形反演装置解决问题的原理与基于复杂起伏地表的陆上地震数据全波形反演方法相似，因此基于复杂起伏地表的陆上地震数据全波形反演装置的实施可以参见基于复杂起伏地表的陆上地震数据全波形反演方法实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的***较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。
本发明的实施例提供一种能够实现基于复杂起伏地表的陆上地震数据全波形反演方法的基于复杂起伏地表的陆上地震数据全波形反演装置的具体实施方式，参见图20，基于复杂起伏地表的陆上地震数据全波形反演装置具体包括如下内容：
几何模型生成单元10，用于对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；
模拟数据生成单元20，用于根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；
震源场生成单元30，用于根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；
传播波场生成单元40，用于根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；
模型更新单元50，用于根据所述多尺度反向传播波场更新所述速度模型。
一实施例中，参见图21，所述几何模型生成单元10包括：
地震数据解析模块101，用于解析所述目的工区的地震数据，以获取所述地震数据的局部坐标采集观测***；所述局部坐标采集观测***包括：震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程；
地表坐标计算模块102，用于根据所述震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程，计算所述目标工区的地表坐标；
几何模型生成模块103，用于利用非规则网格离散方法，根据所述地表坐标以及工区范围生成所述目标工区的几何模型。
一实施例中，参见图22，所述模拟数据生成单元20包括：
子波求取模块201，用于根据所述地震数据求取地震子波；
模拟数据生成模块202，用于利用阻尼间断伽辽金数值模拟方法，根据所述地震子波、所述速度模型、所述目标工区非规则网格下的密度模型、非规则网格节点、有限元数据以及所述地震子波对所述目标工区进行正演模拟，以生成多尺度正演模拟数据。
一实施例中，所述传播波场生成单元具体用于利用所述阻尼间断伽辽金数值模拟方法以及频率域地震波数学物理方程，根据速度模型以及伴随震源场求解所述多尺度反向传播波场。
一实施例中，参见图23，所述模型更新单元50包括：
梯度生成模块501，用于利用伴随状态法，根据所述多尺度反向传播波场以及多尺度正演模拟数据，生成所述速度模型的梯度以及对角线海森矩阵；
模型更新模块502，用于利用共轭梯度法优化方法以及递进尺度方式，根据所述梯度以及对角线海森矩阵更新所述速度模型。
从上述描述可知，本发明实施例提供基于复杂起伏地表的陆上地震数据全波形反演装置，首先对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；本发明实施例可以提供一种适用于复杂起伏地表的陆上地震资料全波形反演方法，可提高陆上双复杂探区深度域速度建模精度，从而改善成像质量。
本申请的实施例还提供能够实现上述实施例中的基于复杂起伏地表的陆上地震数据全波形反演方法中全部步骤的一种电子设备的具体实施方式，参见图24，电子设备具体包括如下内容：
处理器(processor)1201、存储器(memory)1202、通信接口(CommunicationsInterface)1203和总线1204；
其中，处理器1201、存储器1202、通信接口1203通过总线1204完成相互间的通信；通信接口1203用于实现服务器端设备、采集设备以及用户端设备等相关设备之间的信息传输。
处理器1201用于调用存储器1202中的计算机程序，处理器执行计算机程序时实现上述实施例中的基于复杂起伏地表的陆上地震数据全波形反演方法中的全部步骤，例如，处理器执行计算机程序时实现下述步骤：
步骤100：对目标工区进行非规则网格离散，以生成所述目标工区的几何模型。
步骤200：根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据。
步骤300：根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场。
步骤400：根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场。
步骤500：根据所述多尺度反向传播波场更新所述速度模型。
本申请的实施例还提供能够实现上述实施例中的基于复杂起伏地表的陆上地震数据全波形反演方法中全部步骤的一种计算机可读存储介质，计算机可读存储介质上存储有计算机程序，该计算机程序被处理器执行时实现上述实施例中的基于复杂起伏地表的陆上地震数据全波形反演方法的全部步骤，例如，处理器执行计算机程序时实现下述步骤：
步骤100：对目标工区进行非规则网格离散，以生成所述目标工区的几何模型。
步骤200：根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据。
步骤300：根据所述正演模拟数据求解目的工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场。
步骤400：根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场。
步骤500：根据所述多尺度反向传播波场更新所述速度模型。
本领域内的技术人员应明白，本发明的实施例可提供为方法、***、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于复杂起伏地表的陆上地震数据全波形反演方法，其特征在于，包括：
对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；
根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；
根据所述正演模拟数据求解目标工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；
根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；
根据所述多尺度反向传播波场更新所述速度模型；
所述根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据，包括：
根据所述地震数据求取地震子波；
利用阻尼间断伽辽金数值模拟方法，根据所述地震子波、所述速度模型、所述目标工区非规则网格下的密度模型、非规则网格节点、有限元数据以及所述地震子波对所述目标工区进行正演模拟，以生成多尺度正演模拟数据；
所述利用阻尼间断伽辽金数值模拟方法，根据所述地震子波、所述速度模型、所述目标工区非规则网格下的密度模型、非规则网格节点、有限元数据以及所述地震子波对所述目标工区进行正演模拟，以生成多尺度正演模拟数据，包括：根据元素、节点数据、初始速度模型和密度模型，使用线性插值得到非规则网格点对应的速度值和密度值，构造出非规则网格下速度模型和密度模型；
对预处理后的地震记录数据进行傅里叶变换，将其由时间域变换到频率域中，并根据尺度个数，尺度范围，得到不同尺度下的地震数据分量；
对不同尺度的地震数据分量进行线性反演得到频率域多尺度子波。
2.根据权利要求1所述的陆上地震数据全波形反演方法，其特征在于，所述对目标工区进行非规则网格离散，以生成所述目标工区的几何模型包括：
解析所述目标工区的地震数据，以获取所述地震数据的局部坐标采集观测***；所述局部坐标采集观测***包括：震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程；
根据所述震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程，计算所述目标工区的地表坐标；
利用非规则网格离散方法，根据所述地表坐标以及工区范围生成所述目标工区的几何模型。
3.根据权利要求2所述的陆上地震数据全波形反演方法，其特征在于，所述根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场，包括：
利用所述阻尼间断伽辽金数值模拟方法以及频率域地震波数学物理方程，根据速度模型以及伴随震源场求解所述多尺度反向传播波场。
4.根据权利要求3所述的陆上地震数据全波形反演方法，其特征在于，所述根据所述多尺度反向传播波场更新所述速度模型包括：
利用伴随状态法，根据所述多尺度反向传播波场以及多尺度正演模拟数据，生成所述速度模型的梯度以及对角线海森矩阵；
利用共轭梯度法优化方法以及递进尺度方式，根据所述梯度以及对角线海森矩阵更新所述速度模型。
5.一种基于复杂起伏地表的陆上地震数据全波形反演装置，其特征在于，包括：
几何模型生成单元，用于对目标工区进行非规则网格离散，以生成所述目标工区的几何模型；
模拟数据生成单元，用于根据所述几何模型对所述目标工区进行正演模拟，以生成正演模拟数据；
震源场生成单元，用于根据所述正演模拟数据求解目标工区的匹配波场的目标函数以及伴随震源，以生成伴随震源场；
传播波场生成单元，用于根据所述目标工区非规则网格下的速度模型以及所述伴随震源场生成所述伴随震源场的多尺度反向传播波场；
模型更新单元，用于根据所述多尺度反向传播波场更新所述速度模型；
所述模拟数据生成单元包括：
子波求取模块，用于根据所述地震数据求取地震子波；
模拟数据生成模块，用于利用阻尼间断伽辽金数值模拟方法，根据所述地震子波、所述速度模型、所述目标工区非规则网格下的密度模型、非规则网格节点、有限元数据以及所述地震子波对所述目标工区进行正演模拟，以生成多尺度正演模拟数据；
所述模拟数据生成单元具体用于：
根据元素、节点数据、初始速度模型和密度模型，使用线性插值得到非规则网格点对应的速度值和密度值，构造出非规则网格下速度模型和密度模型；
对预处理后的地震记录数据进行傅里叶变换，将其由时间域变换到频率域中，并根据尺度个数，尺度范围，得到不同尺度下的地震数据分量；
对不同尺度的地震数据分量进行线性反演得到频率域多尺度子波。
6.根据权利要求5所述的陆上地震数据全波形反演装置，其特征在于，所述几何模型生成单元包括：
地震数据解析模块，用于解析所述目标工区的地震数据，以获取所述地震数据的局部坐标采集观测***；所述局部坐标采集观测***包括：震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程；
地表坐标计算模块，用于根据所述震道对应的炮号、震源点坐标、震源点高程、震源点深度、检波点坐标以及检波点高程，计算所述目标工区的地表坐标；
几何模型生成模块，用于利用非规则网格离散方法，根据所述地表坐标以及工区范围生成所述目标工区的几何模型。
7.根据权利要求6所述的陆上地震数据全波形反演装置，其特征在于，所述传播波场生成单元具体用于利用所述阻尼间断伽辽金数值模拟方法以及频率域地震波数学物理方程，根据速度模型以及伴随震源场求解所述多尺度反向传播波场。
8.根据权利要求7所述的陆上地震数据全波形反演装置，其特征在于，所述模型更新单元包括：
梯度生成模块，用于利用伴随状态法，根据所述多尺度反向传播波场以及多尺度正演模拟数据，生成所述速度模型的梯度以及对角线海森矩阵；
模型更新模块，用于利用共轭梯度法优化方法以及递进尺度方式，根据所述梯度以及对角线海森矩阵更新所述速度模型。
9.一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现权利要求1至4任一项所述基于复杂起伏地表的陆上地震数据全波形反演方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该计算机程序被处理器执行时实现权利要求1至4任一项所述基于复杂起伏地表的陆上地震数据全波形反演方法的步骤。
CN202010478692.6A 2020-05-29 2020-05-29 基于复杂起伏地表的陆上地震数据全波形反演方法及装置 Active CN113740901B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010478692.6A CN113740901B (zh) 2020-05-29 2020-05-29 基于复杂起伏地表的陆上地震数据全波形反演方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010478692.6A CN113740901B (zh) 2020-05-29 2020-05-29 基于复杂起伏地表的陆上地震数据全波形反演方法及装置

Publications (2)

Publication Number Publication Date
CN113740901A CN113740901A (zh) 2021-12-03
CN113740901B true CN113740901B (zh) 2024-01-30

Family

ID=78724973

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010478692.6A Active CN113740901B (zh) 2020-05-29 2020-05-29 基于复杂起伏地表的陆上地震数据全波形反演方法及装置

Country Status (1)

Country Link
CN (1) CN113740901B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115220091B (zh) * 2022-02-22 2023-04-28 中国科学院地质与地球物理研究所 一种地质导向的非规则观测***确定方法及***
CN114935319B (zh) * 2022-04-08 2023-02-21 北京大学 多偏移距震电频谱比值获取方法及用于监测潜水面的方法
CN114994749A (zh) * 2022-05-30 2022-09-02 佟小龙 陆地地球物理勘探方法、电子设备及可读存储介质
CN117420596A (zh) * 2022-07-19 2024-01-19 中国石油天然气集团有限公司 地下介质速度的确定方法、装置、设备
CN116819602B (zh) * 2023-07-12 2024-02-09 中国矿业大学 一种深度学习优化的变密度声波方程全波形反演方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107765302A (zh) * 2017-10-20 2018-03-06 吉林大学 不依赖震源子波的时间域单频波形走时反演方法
CN109188512A (zh) * 2018-09-17 2019-01-11 中国石油大学(华东) 基于非规则扇形网格剖分的起伏隧道空间正演模拟***及方法
DE102018108544A1 (de) * 2018-01-31 2019-08-01 Bauhaus-Universität Weimar Detektion und Lokalisierung von Bauwerksschäden mit vollständiger Wellenformumkehrung

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107765302A (zh) * 2017-10-20 2018-03-06 吉林大学 不依赖震源子波的时间域单频波形走时反演方法
DE102018108544A1 (de) * 2018-01-31 2019-08-01 Bauhaus-Universität Weimar Detektion und Lokalisierung von Bauwerksschäden mit vollständiger Wellenformumkehrung
CN109188512A (zh) * 2018-09-17 2019-01-11 中国石油大学(华东) 基于非规则扇形网格剖分的起伏隧道空间正演模拟***及方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
全波形反演中三种优化方法对比研究;刘隼;《中国优秀硕士学位论文全文数据库 基础科学辑》(第第1期期);第A011-712页 *
基于非规则网格声波正演的时间域全波形反演;魏哲枫 等;《地球物理学报》;第第57卷卷(第第2期期);第586-594页 *

Also Published As

Publication number Publication date
CN113740901A (zh) 2021-12-03

Similar Documents

Publication Publication Date Title
CN113740901B (zh) 基于复杂起伏地表的陆上地震数据全波形反演方法及装置
KR102020759B1 (ko) Q-보상된 전 파동장 반전
EP2491428B1 (en) Full-waveform inversion in the traveltime domain
EP3063562B1 (en) Methods of subsurface exploration, computer program product and computer-readable storage medium
US10317548B2 (en) Reflection seismic data Q tomography
EP2810101B1 (en) Improving efficiency of pixel-based inversion algorithms
NO347906B1 (en) Stratigraphic function
CN104570082B (zh) 一种基于格林函数表征的全波形反演梯度算子的提取方法
Qu et al. Elastic full-waveform inversion for surface topography
US9952341B2 (en) Systems and methods for aligning a monitor seismic survey with a baseline seismic survey
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
CN113945982B (zh) 用于去除低频和低波数噪声以生成增强图像的方法和***
RU2570827C2 (ru) Гибридный способ для полноволновой инверсии с использованием способа одновременных и последовательных источников
Willemsen et al. An efficient coupled acoustic-elastic local solver applied to phase inversion
Ma et al. Topography-dependent eikonal traveltime tomography for upper crustal structure beneath an irregular surface
Li et al. Waveform inversion of seismic first arrivals acquired on irregular surface
Gibson Jr et al. Modeling and velocity analysis with a wavefront-construction algorithm for anisotropic media
CN115373022A (zh) 一种基于振幅相位校正的弹性波场Helmholtz分解方法
Wang et al. Rock Fracture Monitoring Based on High‐Precision Microseismic Event Location Using 3D Multiscale Waveform Inversion
Song et al. Multi-scale seismic full waveform inversion in the frequency-domain with a multi-grid method
CN111175822A (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
Du et al. Full waveform inversion based on well data constraint
Roberts et al. Investigation into the use of 2D elastic waveform inversion from look‐ahead walk‐away VSP surveys
Majdi et al. Application of finite difference eikonal solver for traveltime computation in forward modeling and migration
Krieger 2D Eleastic Full-waveform Inversion of Land Seismic Data with Topographic Variations

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