CN113504566A - 基于波动方程的地震反演方法、系统、装置及介质 - Google Patents

基于波动方程的地震反演方法、系统、装置及介质 Download PDF

Info

Publication number
CN113504566A
CN113504566A CN202110609064.1A CN202110609064A CN113504566A CN 113504566 A CN113504566 A CN 113504566A CN 202110609064 A CN202110609064 A CN 202110609064A CN 113504566 A CN113504566 A CN 113504566A
Authority
CN
China
Prior art keywords
seismic
frequency
wave
determining
travel time
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.)
Granted
Application number
CN202110609064.1A
Other languages
English (en)
Other versions
CN113504566B (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.)
Shanghai Qingfeng Zhiyuan Geophysical Geological Exploration Technology Co ltd
CNOOC China Ltd Zhanjiang Branch
Southern Marine Science and Engineering Guangdong Laboratory Zhanjiang
Original Assignee
Shanghai Qingfeng Zhiyuan Geophysical Geological Exploration Technology Co ltd
CNOOC China Ltd Zhanjiang Branch
Southern Marine Science and Engineering Guangdong Laboratory Zhanjiang
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 Shanghai Qingfeng Zhiyuan Geophysical Geological Exploration Technology Co ltd, CNOOC China Ltd Zhanjiang Branch, Southern Marine Science and Engineering Guangdong Laboratory Zhanjiang filed Critical Shanghai Qingfeng Zhiyuan Geophysical Geological Exploration Technology Co ltd
Priority to CN202110609064.1A priority Critical patent/CN113504566B/zh
Publication of CN113504566A publication Critical patent/CN113504566A/zh
Application granted granted Critical
Publication of CN113504566B publication Critical patent/CN113504566B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/303Analysis for determining velocity profiles or travel times
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282Application of seismic models, synthetic seismograms
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/622Velocity, density or impedance
    • G01V2210/6222Velocity; travel time
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/67Wave propagation modeling
    • G01V2210/675Wave equation; Green's functions

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

基于波动方程的地震反演方法、系统、装置及介质
技术领域
本发明涉及地震速度建模技术领域,尤其是基于波动方程的地震反演方法、系统、装置及介质。
背景技术
在传统的波动方程走时反演(Luo and Schuster,1991;Tromp et al.,2005)中,目标函数的梯度可以采用伴随状态理论推导,通常可以表示为一个入射场与伴随场的零延迟互相关(Lailly,1983;Liu and Tromp,2006,2008;Fichtne et al.,2006)。考虑到入射场是沿时间正向传播的,而伴随波场逆时传播的。为同时获得某一时刻的入射场和伴随场,需要采取某些策略从而兼顾计算和存储:如最优检查点方法(optimal checkpointingmethod,Symes,2007)、随机边界条件(random boundary condition,Clapp,2009),波场重构策略(Dussaud et al.,2008;Feng et al.,2011)等。上述方法中,波场重构策略在计算精度、重计算比和数据I/O量三个方面取得了较好的平衡,因此广泛应用于波形反演及逆时偏移(RTM)中。
然而,对于三维波动方程初至层析,即使采用了波场重构策略计算梯度,其计算量仍然比射线层析提高了几个数量级,因此严重限制了该技术的实用化。
发明内容
有鉴于此,为至少部分解决上述技术问题之一,本发明实施例目的在于提供一种计算量小,更为轻量的基于波动方程的地震反演方法,以及本申请还提供了能够对应实现该方法的系统、装置及计算机可读的存储介质。
第一方面,本申请的技术方案提供了基于波动方程的地震反演方法,其步骤包括:
获取并根据地震观测参数与初始速度模型参数确定单频敏感度核函数的计算参数;
生成随机边界的速度模型;
根据所述计算参数和所述速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场;
根据所述合成地震记录确定走时伴随震源,根据所述走时伴随震源确定走时目标泛函的单频梯度;
从所述单频梯度提取低波数部分,根据所述低波数进行地震速度反演。
在本申请方案的一种可行的实施例中,所述生成随机边界的速度模型这一步骤,其包括:
构建三维速度模型,将所述三维速度模型的第一表面设置为刚性地标,将除第一表面外的其余表面设置为随机速度模型;
所述随机速度模型为随机速度层,所述随机速度层由若干随机速度网格构成。
在本申请方案的一种可行的实施例中,所述根据所述计算参数和所述速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场这一步骤,其包括:
根据所述计算参数中的地震子波、具有随机边界的所述速度模型以及正演模拟参数,进行波动方程正演;
对所述波动方程正演中时间步的波场进行单频离散傅里叶积分得到所述合成地震记录和所述频率域入射波场。
在本申请方案的一种可行的实施例中,所述根据所述合成地震记录确定走时伴随震源,根据所述走时伴随震源确定走时目标泛函的单频梯度这一步骤,其包括:
根据所述走时伴随震源确定边界条件,通过波动方程逆时传播确定伴随波场;
在所述波动方程逆时传播中对所述伴随波场进行单频离散傅里叶积分,得到频域伴随波场;
根据所述频域伴随波场的入射波场和伴随波场,确定所述走时目标泛函的单频梯度。
在本申请方案的一种可行的实施例中,所述根据所述频域伴随波场的入射波场和伴随波场,确定所述走时目标泛函的单频梯度这一步骤,其包括:
获取若干单炮道集的单频梯度,对若干所述单炮道集的单频梯度进行求和,得到
所述走时目标泛函的单频梯度。
在本申请方案的一种可行的实施例中,所述从所述单频梯度提取低波数部分,根据所述低波数进行地震速度反演这一步骤,其包括:
通过高斯光滑滤波器,对所述单频梯度进行光滑处理,提取得到所述低波数部分。
在本申请方案的一种可行的实施例中,所述获取并根据地震观测参数与初始速度模型参数确定单频敏感度核函数的计算参数,其包括以下步骤至少之一:
确定震源即检波器坐标;
确定所述单频敏感度核函数的频率;
确定地震子波。
第二方面,本发明的技术方案还提供基于波动方程的地震反演系统,其包括:
参数确定模块,用于获取并根据地震观测参数与初始速度模型参数确定单频敏感度核函数的计算参数;
模型生成模块,用于生成随机边界的速度模型;
数据处理模块,用于根据所述计算参数和所述速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场,以及用于根据所述合成地震记录确定走时伴随震源,根据所述走时伴随震源确定走时目标泛函的单频梯度;
地震反演模块,用于从所述单频梯度提取低波数部分,根据所述低波数进行地震速度反演。
第三方面,本发明的技术方案还提供一种基于波动方程的地震反演装置,其包括:
至少一个处理器;
至少一个存储器,用于存储至少一个程序;
当至少一个程序被至少一个处理器执行,使得至少一个处理器运行第一方面中的基于波动方程的地震反演方法。
第四方面,本发明的技术方案还提供了一种存储介质,其中存储有处理器可执行的程序,处理器可执行的程序在由处理器执行时用于运行第一方面中的方法。
本发明的优点和有益效果将在下面的描述中部分给出,其他部分可以通过本发明的具体实施方式了解得到:
本申请技术方案针对波动方程层析中计算和存储问题,提出了基于随机边界条件的单频走时敏感度核函数,并用于走时误差泛函梯度的高效计算,并通过提取单频梯度中的低波数部分用于背景速度反演;本方案相比于基于波场重构框架的梯度计算,提出的单频梯度计算方法对内存需求极低,且计算量至少减少1/3;并且,随机边界的引入使得正演算法大幅简化,因而更适用于GPU加速计算,相对于基于波场重构框架的波动方程走时层析技术,本发明在计算效率、内存占用及算法复杂度方面具有显著优势。
附图说明
为了更清楚地说明本申请实施例中的技术方案,下面将对实施例描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本申请的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的基于波动方程的地震反演方法的步骤流程图;
图2为本发明实施例中含有随机边界的速度模型示意图;
图3为本发明实施例中基于随机边界条件计算的单频走时敏感度核函数的示意图;
图4为本发明实施例中基于单频走时敏感度核函数的走时误差泛函的梯度的示意图;
图5为本发明实施例中基于有限频走时敏感度核函数的走时误差泛函的梯度的示意图;
图6为本发明实施例中经过高斯平滑后的单频梯度的示意图;
图7为本发明实施例中经过高斯平滑后的有限频梯度的示意图。
具体实施方式
下面详细描述本发明的实施例,实施例的示例在附图中示出,其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的,仅用于解释本发明,而不能理解为对本发明的限制。对于以下实施例中的步骤编号,其仅为了便于阐述说明而设置,对步骤之间的顺序不做任何限定,实施例中的各步骤的执行顺序均可根据本领域技术人员的理解来进行适应性调整。
在现有技术中,对于三维波动方程初至层析,即使采用了波场重构策略计算梯度,其计算量仍然比射线层析提高了几个数量级,因此严重限制了该技术的实用化。所以,本申请的技术方案针对波动方程层析中计算和存储问题,创造性地提出了基于随机边界条件的单频走时敏感度核函数,并用于走时误差泛函梯度的高效计算。
第一方面,如图1所示,本申请提供可以基于第一方面中的系统实现的基于波动方程的地震反演方法,其步骤包括S100-S500:
S100、获取并根据地震观测参数与初始速度模型参数确定单频敏感度核函数的计算参数;
其中,地震观测参数是通过地震观测系统采集得到的相关系数,地震观测系统可以描述地震勘探中激发点和接收点排列之间的相对空间位置关系的布置方式。针对不同的地震勘探方法,采用不同的观测系统。得到的相关参数包括但不限于最大偏移距、炮间距、偏移孔径等;初始速度模型参数是指初始速度模型的相关参数,其包括但不限于速度模型采样点数与采样间隔等;所确定的单频敏感度核函数的计算参数,其包括但不限于震源及检波器坐标、核函数的频率、地震子波等。具体地,实施例根据地震观测系统及反演目标的深度及尺度,得到相应的单频敏感度核函数的计算参数作为反演参数,可以理解的是,实施例中,计算参数其还可以包括单频梯度对应的频率,高斯平滑参数等。
S200、生成随机边界的速度模型;
具体地,实施例将生成单炮孔径内包含随机边界的速度模型,以利用输入的地震子波、速度模型及正演参数进行波动方程正演。示例性地,实施例中该速度模型可以选用三维速度模型。
S300、根据计算参数和速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场;
具体地,实施例利用步骤S100中确定的单频敏感度核函数的计算参数,例如输入的地震子波、正演模拟参数结合包含随机边界的速度模型进行波动方程正演。在正演过程中,对每个时间步的波场进行单频离散傅里叶积分。正演结束后得到合成地震记录以及频率域入射波场。
S400、根据合成地震记录确定走时伴随震源,根据走时伴随震源确定走时目标泛函的单频梯度;
具体地,实施例根据步骤S300中得到的合成地震记录进一步计算走时伴随震源,并将走时伴随震源作为边界条件,在进行波动方程逆时传播计算伴随波场,并可以在逆时传播过程中对伴随波场进行单频离散傅里叶积分,得到频域伴随波场。实施例再根据频域的入射波场和伴随波场,计算单炮道集的单频梯度,根据若干单炮道集的单频梯度求和得到走时目标泛函的单频梯度。
S500、从单频梯度提取低波数部分,根据低波数进行地震速度反演;
具体地,实施例利用滤波器,对步骤S400中得到的单频梯度进行光滑处理,提取得到单频梯度中的低波数部分,将其作为波动方程初至层析目标函数的梯度用于地震速度反演。
在一些可选的实施例中,生成随机边界的速度模型这一步骤S200,其可以包括步骤S210-S220:
S210、构建三维速度模型,将三维速度模型的第一表面设置为刚性地标,将除第一表面外的其余表面设置为随机速度模型;
S220、随机速度模型为随机速度层,随机速度层由若干随机速度网格构成;
具体地,实施例中可以选用三维速度模型,对于三维速度模型,实施例将其上表面设置为刚性地表,无随机速度;其余5个面为随机速度模型,随机速度层的厚度为用户输入参数,通常取30-50个网格,随机速度层每个网格上的速度为随机生成的,速度取值范围为[Vmin,Vmax],其中Vmin和Vmax分别是输入速度模型的最小速度和最大速度。
在一些可选的实施例中,根据计算参数和速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场这一步骤S300,可以进一步细分得到步骤S310-S320:
S310、根据计算参数中的地震子波、具有随机边界的速度模型以及正演模拟参数,进行波动方程正演;
S320、对波动方程正演中时间步的波场进行单频离散傅里叶积分得到合成地震记录和频率域入射波场;
具体地,实施例利用输入的地震子波、包含随机边界的速度模型及正演模拟参数进行波动方程正演。在本实施例中,具体采用标量声波方程正演模拟,并对每个时间步的波场进行单频离散傅里叶积分;正演结束后得到合成地震记录以及频率域入射波场。离散傅里叶积分公式如下:
Figure BDA0003095262320000061
在公式(1)中,u(x,kΔt;xs)为波动方程正演模拟得到的时间-空间域全波场,xs为震源位置,x为地下空间任一点的坐标,exp(-i2πfkΔt)为离散Fourier积分核,u(x,f;xs)为频率域入射波场,f为频率。合成地震记录定义为时间-空间域全波场在检波器位置xr处的空间采样,其表示为
Figure BDA0003095262320000062
需要补充说明的是,本实施例中,以二阶常密度标量声波方程的有限差分数值解为例,给出波动方程正演模拟及逆时传播对应的数学描述:
(1)波动方程正演模拟,需要求解如下初值问题:
Figure BDA0003095262320000063
在公式(2)中,x=(x,y,z),f(t)为震源子波,xs为震源坐标,合成地震记录定义为地震波场u(x,t)在检波器位置xr处的采样,表示为:
Figure BDA0003095262320000064
(2)波动方程逆时传播,需要求解如下边值问题(伴随震源作为边值条件):
Figure BDA0003095262320000065
在一些可选的实施例中,根据合成地震记录确定走时伴随震源,根据走时伴随震源确定走时目标泛函的单频梯度这一步骤S400,其包括步骤S410-S430:
S410、根据走时伴随震源确定边界条件,通过波动方程逆时传播确定伴随波场;
S420、在波动方程逆时传播中对伴随波场进行单频离散傅里叶积分,得到频域伴随波场;
S430、根据频域伴随波场的入射波场和伴随波场,确定走时目标泛函的单频梯度;
具体地,实施例根据公式(2),通过合成地震记录计算走时伴随震源:
Figure BDA0003095262320000066
然后,实施例将走时伴随震源作为边界条件,进行波动方程逆时传播计算伴随波场。根据公式(5),在逆时传播过程中对伴随波场进行单频离散傅里叶积分,得到频率域伴随波场λ(x,f):
Figure BDA0003095262320000071
在公式(5)中,λ(x,kΔt)代表走时伴随震源逆时传播得到的伴随波场。随后实施例再利用频率域的入射波场和伴随波场,计算单炮道集的单频梯度:
grad(x;xs,f)=Re{(u(x,f;xs))*λ(x,f)} (6)
在公式(6),Re表示取复数的实部,上标*代表复数共轭。根据公式(6)算得到单炮道集的单频梯度确定走时目标泛函的单频梯度。
在本实施例中,步骤S430、根据频域伴随波场的入射波场和伴随波场,确定走时目标泛函的单频梯度;其还可以包括步骤S430a、获取若干单炮道集的单频梯度,对若干单炮道集的单频梯度进行求和,得到走时目标泛函的单频梯度。
具体地,实施例重复步骤S210-S430,计算观测系统中所有炮的单频梯度grad(x;xs i,f),i=1,…,Nshot,其中,Nshot为炮数,并进行求和得到走时目标泛函的单频梯度:
Figure BDA0003095262320000072
在一些其他的可选实施例中,从单频梯度提取低波数部分,根据低波数进行地震速度反演这一步骤S500,其具体为:通过高斯光滑滤波器,对单频梯度进行光滑处理,提取得到低波数部分。
具体地,实施例采用高斯光滑滤波器对步骤S430a中得到的单频梯度进行光滑处理,提取单频梯度中的低波数部分。实施例中所采用的高斯光滑滤波器如下:
Figure BDA0003095262320000073
在公式(8)中,σx,σy,σz分别为高斯光滑滤波器在x,y,z方向上的平滑半径。通过将高斯光滑滤波器与单频梯度卷积:G(x,σxyz)*grad(x;f),实现了对梯度进行光滑化处理。
更为具体地,实施例的完整实施例过程为:
首先,实施例生成一个包含随机边界的速度模型,如图1所示,上表面设置为刚性地表,无随机速度,其余5个面为随机速度模型,随机速度层的厚度为50个网格。
其次,输入震源、检波器坐标及核函数频率,计算单频走时敏感度核函数,如图2所示,在第一菲涅尔体内,单频核函数的形态与传统的有限频核函数相近;而在第一菲涅尔体外,单频核函数呈明显的空间震荡特性,符合理论预测。结合附图2和附图3,展示了单频核函数的形态。
然后,如图4所示,实施例将单频核函数应用于梯度计算,可以得到走时目标泛函的单频梯度,作为对比,如图5所示,展示了走时目标泛函的有限频梯度。二者均在震源附近的存在较强的高波数噪声。
最后,实施例对单频梯度进行光滑处理得到其低波数部分,如图6所示,可以看出,震源附近的高波数噪声已经完全消除,单频梯度表现出较好的背景速度更新特性。作为对比,如图7所示,展示了高斯光滑后的有限频梯度。对比图6和图7可以看出,光滑后的单频梯度与传统有限频梯度的低波数部分形态相似,因此利用单频梯度进行背景速度反演具有可行性。
第二方面,本申请的技术方案还提供基于波动方程的地震反演系统,其包括:
参数确定模块,用于获取并根据地震观测参数与初始速度模型参数确定单频敏感度核函数的计算参数;
模型生成模块,用于生成随机边界的速度模型;
数据处理模块,用于根据计算参数和速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场,以及用于根据合成地震记录确定走时伴随震源,根据走时伴随震源确定走时目标泛函的单频梯度;
地震反演模块,用于从单频梯度提取低波数部分,根据低波数进行地震速度反演。
第三方面,本申请的技术方案还提供基于波动方程的地震反演装置,其包括至少一个处理器;至少一个存储器,用于存储至少一个程序;当至少一个程序被至少一个处理器执行,使得至少一个处理器运行如第一方面中的基于波动方程的地震反演方法。
本发明实施例还提供了一种存储介质内存储有程序,程序被处理器执行,实现如第一方面中的方法。
从上述具体的实施过程,可以总结出,本发明所提供的技术方案相较于现有技术存在以下优点或优势:
本申请技术方案在计算和存储方面,相比于基于波场重构框架(需要三次波传播且需要处理数值吸收边界)的梯度计算策略,本发明提出的单频梯度计算方法对内存需求极低(相对于波动方程正演,仅需要增加两个单频波场以及一个三维成像体),且计算量至少减少1/3(仅需要两次波传播且无需处理数值吸收边界)。在算法复杂度方面,随机边界的引入使得正演算法大幅简化(无需单独处理吸收边界加角点反射),因而更适用于GPU加速计算。因此,相对于基于波场重构框架的波动方程走时层析技术(即传统方法),本申请的技术方案在计算效率、内存占用及算法复杂度方面具有显著优势。
在一些可选择的实施例中,在方框图中提到的功能/操作可以不按照操作示图提到的顺序发生。例如,取决于所涉及的功能/操作,连续示出的两个方框实际上可以被大体上同时地执行或所述方框有时能以相反顺序被执行。此外,在本发明的流程图中所呈现和描述的实施例以示例的方式被提供,目的在于提供对技术更全面的理解。所公开的方法不限于本文所呈现的操作和逻辑流程。可选择的实施例是可预期的,其中各种操作的顺序被改变以及其中被描述为较大操作的一部分的子操作被独立地执行。
此外,虽然在功能性模块的背景下描述了本发明,但应当理解的是,除非另有相反说明,功能和/或特征中的一个或多个可以被集成在单个物理装置和/或软件模块中,或者一个或多个功能和/或特征可以在单独的物理装置或软件模块中被实现。还可以理解的是,有关每个模块的实际实现的详细讨论对于理解本发明是不必要的。更确切地说,考虑到在本文中公开的装置中各种功能模块的属性、功能和内部关系的情况下,在工程师的常规技术内将会了解该模块的实际实现。因此,本领域技术人员运用普通技术就能够在无需过度试验的情况下实现在权利要求书中所阐明的本发明。还可以理解的是,所公开的特定概念仅仅是说明性的,并不意在限制本发明的范围,本发明的范围由所附权利要求书及其等同方案的全部范围来决定。
在流程图中表示或在此以其他方式描述的逻辑和/或步骤,例如,可以被认为是用于实现逻辑功能的可执行指令的定序列表,可以具体实现在任何计算机可读介质中,以供指令执行系统、装置或设备(如基于计算机的系统、包括处理器的系统或其他可以从指令执行系统、装置或设备取指令并执行指令的系统)使用,或结合这些指令执行系统、装置或设备而使用。
在本说明书的描述中,参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中,对上述术语的示意性表述不一定指的是相同的实施例或示例。而且,描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
尽管已经示出和描述了本发明的实施例,本领域的普通技术人员可以理解:在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型,本发明的范围由权利要求及其等同物限定。
以上是对本发明的较佳实施进行了具体说明,但本发明并不限于上述实施例,熟悉本领域的技术人员在不违背本发明精神的前提下还可做作出种种的等同变形或替换,这些等同的变形或替换均包含在本申请权利要求所限定的范围内。

Claims (10)

1.基于波动方程的地震反演方法,其特征在于,
获取并根据地震观测参数与初始速度模型参数确定单频敏感度核函数的计算参数;
生成随机边界的速度模型;
根据所述计算参数和所述速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场;
根据所述合成地震记录确定走时伴随震源,根据所述走时伴随震源确定走时目标泛函的单频梯度;
从所述单频梯度提取低波数部分,根据所述低波数进行地震速度反演。
2.根据权利要求1所述的基于波动方程的地震反演方法,其特征在于,所述生成随机边界的速度模型这一步骤,其包括:
构建三维速度模型,将所述三维速度模型的第一表面设置为刚性地标,将除第一表面外的其余表面设置为随机速度模型;
所述随机速度模型为随机速度层,所述随机速度层由若干随机速度网格构成。
3.根据权利要求2所述的基于波动方程的地震反演方法,其特征在于,所述根据所述计算参数和所述速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场这一步骤,其包括:
根据所述计算参数中的地震子波、具有随机边界的所述速度模型以及正演模拟参数,进行波动方程正演;
对所述波动方程正演中时间步的波场进行单频离散傅里叶积分得到所述合成地震记录和所述频率域入射波场。
4.根据权利要求1所述的基于波动方程的地震反演方法,其特征在于,所述根据所述合成地震记录确定走时伴随震源,根据所述走时伴随震源确定走时目标泛函的单频梯度这一步骤,其包括:
根据所述走时伴随震源确定边界条件,通过波动方程逆时传播确定伴随波场;
在所述波动方程逆时传播中对所述伴随波场进行单频离散傅里叶积分,得到频域伴随波场;
根据所述频域伴随波场的入射波场和伴随波场,确定所述走时目标泛函的单频梯度。
5.根据权利要求4所述的基于波动方程的地震反演方法,其特征在于,所述根据所述频域伴随波场的入射波场和伴随波场,确定所述走时目标泛函的单频梯度这一步骤,其包括:获取若干单炮道集的单频梯度,对若干所述单炮道集的单频梯度进行求和,得到所述走时目标泛函的单频梯度。
6.根据权利要求1所述的基于波动方程的地震反演方法,其特征在于,所述从所述单频梯度提取低波数部分,根据所述低波数进行地震速度反演这一步骤,其包括:
通过高斯光滑滤波器,对所述单频梯度进行光滑处理,提取得到所述低波数部分。
7.根据权利要求1-6任一项所述的基于波动方程的地震反演方法,其特征在于,所述获取并根据地震观测参数与初始速度模型参数确定单频敏感度核函数的计算参数,其包括以下步骤至少之一:
确定震源即检波器坐标;
确定所述单频敏感度核函数的频率;
确定地震子波。
8.基于波动方程的地震反演系统,其特征在于,包括:
参数确定模块,用于获取并根据地震观测参数与初始速度模型参数确定单频敏感度核函数的计算参数;
模型生成模块,用于生成随机边界的速度模型;
数据处理模块,用于根据所述计算参数和所述速度模型进行波动方程正演,得到合成地震记录以及频率域入射波场,以及用于根据所述合成地震记录确定走时伴随震源,根据所述走时伴随震源确定走时目标泛函的单频梯度;
地震反演模块,用于从所述单频梯度提取低波数部分,根据所述低波数进行地震速度反演。
9.基于波动方程的地震反演装置,其特征在于,包括:
至少一个处理器;
至少一个存储器,用于存储至少一个程序;
当所述至少一个程序被所述至少一个处理器执行,使得所述至少一个处理器运行如权利要求1-7任一项所述的基于波动方程的地震反演方法。
10.一种存储介质,其中存储有处理器可执行的程序,其特征在于,所述处理器可执行的程序在由处理器执行时用于运行如权利要求1-7中任一项所述的基于波动方程的地震反演方法。
CN202110609064.1A 2021-06-01 2021-06-01 基于波动方程的地震反演方法、系统、装置及介质 Active CN113504566B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110609064.1A CN113504566B (zh) 2021-06-01 2021-06-01 基于波动方程的地震反演方法、系统、装置及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110609064.1A CN113504566B (zh) 2021-06-01 2021-06-01 基于波动方程的地震反演方法、系统、装置及介质

Publications (2)

Publication Number Publication Date
CN113504566A true CN113504566A (zh) 2021-10-15
CN113504566B CN113504566B (zh) 2024-04-30

Family

ID=78008729

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110609064.1A Active CN113504566B (zh) 2021-06-01 2021-06-01 基于波动方程的地震反演方法、系统、装置及介质

Country Status (1)

Country Link
CN (1) CN113504566B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116626751A (zh) * 2023-05-26 2023-08-22 北京中矿大地地球探测工程技术有限公司 基于多目标函数的黏弹性参数同步反演方法、装置和设备
CN116699695A (zh) * 2023-08-07 2023-09-05 北京中矿大地地球探测工程技术有限公司 一种基于衰减矫正的反演方法、装置和设备

Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5850622A (en) * 1996-11-08 1998-12-15 Amoco Corporation Time-frequency processing and analysis of seismic data using very short-time fourier transforms
US20130077439A1 (en) * 2011-05-27 2013-03-28 Conocophillips Company Reciprocal method two-way wave equation targeted data selection for seismic acquisition of complex geologic structures
CN103135132A (zh) * 2013-01-15 2013-06-05 中国科学院地质与地球物理研究所 Cpu/gpu协同并行计算的混合域全波形反演方法
US20130155814A1 (en) * 2011-12-20 2013-06-20 Conocophillips Company Critical reflection illuminations analysis
CN103675908A (zh) * 2012-09-21 2014-03-26 中国石油化工股份有限公司 一种海量数据图形处理器的波动方程逆时偏移成像方法
CN104297791A (zh) * 2014-09-25 2015-01-21 中国石油天然气股份有限公司 一种基于地震优势频率的反演方法及系统
WO2015106879A1 (en) * 2014-01-14 2015-07-23 Statoil Petroleum As Full wave reverse time migration
CN105005076A (zh) * 2015-06-02 2015-10-28 中国海洋石油总公司 基于最小二乘梯度更新速度模型的地震波全波形反演方法
CN105572734A (zh) * 2014-10-16 2016-05-11 中国石油化工股份有限公司 一种以逆时偏移算法为引擎的波动方程初至走时层析方法
US20160216389A1 (en) * 2015-01-23 2016-07-28 Advanced Geophysical Technology Inc. Beat tone full waveform inversion
CN107765302A (zh) * 2017-10-20 2018-03-06 吉林大学 不依赖震源子波的时间域单频波形走时反演方法
US20190302293A1 (en) * 2018-03-30 2019-10-03 Cgg Services Sas Methods using travel-time full waveform inversion for imaging subsurface formations with salt bodies
CN110579795A (zh) * 2018-06-08 2019-12-17 中国海洋大学 基于被动源地震波形及其逆时成像的联合速度反演方法
CN110858000A (zh) * 2018-08-24 2020-03-03 中国石油天然气股份有限公司 地震数据重构方法、装置、计算机设备及存储介质
CN110888158A (zh) * 2018-09-10 2020-03-17 中国石油化工股份有限公司 一种基于rtm约束的全波形反演方法

Patent Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5850622A (en) * 1996-11-08 1998-12-15 Amoco Corporation Time-frequency processing and analysis of seismic data using very short-time fourier transforms
US20130077439A1 (en) * 2011-05-27 2013-03-28 Conocophillips Company Reciprocal method two-way wave equation targeted data selection for seismic acquisition of complex geologic structures
US20130155814A1 (en) * 2011-12-20 2013-06-20 Conocophillips Company Critical reflection illuminations analysis
CN103675908A (zh) * 2012-09-21 2014-03-26 中国石油化工股份有限公司 一种海量数据图形处理器的波动方程逆时偏移成像方法
CN103135132A (zh) * 2013-01-15 2013-06-05 中国科学院地质与地球物理研究所 Cpu/gpu协同并行计算的混合域全波形反演方法
WO2015106879A1 (en) * 2014-01-14 2015-07-23 Statoil Petroleum As Full wave reverse time migration
CN104297791A (zh) * 2014-09-25 2015-01-21 中国石油天然气股份有限公司 一种基于地震优势频率的反演方法及系统
CN105572734A (zh) * 2014-10-16 2016-05-11 中国石油化工股份有限公司 一种以逆时偏移算法为引擎的波动方程初至走时层析方法
US20160216389A1 (en) * 2015-01-23 2016-07-28 Advanced Geophysical Technology Inc. Beat tone full waveform inversion
CN105005076A (zh) * 2015-06-02 2015-10-28 中国海洋石油总公司 基于最小二乘梯度更新速度模型的地震波全波形反演方法
CN107765302A (zh) * 2017-10-20 2018-03-06 吉林大学 不依赖震源子波的时间域单频波形走时反演方法
US20190302293A1 (en) * 2018-03-30 2019-10-03 Cgg Services Sas Methods using travel-time full waveform inversion for imaging subsurface formations with salt bodies
CN110579795A (zh) * 2018-06-08 2019-12-17 中国海洋大学 基于被动源地震波形及其逆时成像的联合速度反演方法
CN110858000A (zh) * 2018-08-24 2020-03-03 中国石油天然气股份有限公司 地震数据重构方法、装置、计算机设备及存储介质
CN110888158A (zh) * 2018-09-10 2020-03-17 中国石油化工股份有限公司 一种基于rtm约束的全波形反演方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
ZHANG, SJ: "" Estimating the acoustic intensity of typhoon Shanshan at very low frequencies"", 《OCEANS - MTS/IEEE KOBE TECHNO-OCEANS CONFERENCE (OTO)》, 1 January 2018 (2018-01-01), pages 1 - 4 *
王大为: ""与速度无关的叠前时间偏移成像方法"", 《工程地球物理学报》, 30 September 2018 (2018-09-30), pages 586 - 591 *
田亚静: ""弹性波菲涅尔带走时层析方法研究"", 《中国优秀硕士学位论文全文数据库》, 15 June 2018 (2018-06-15), pages 1 - 74 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116626751A (zh) * 2023-05-26 2023-08-22 北京中矿大地地球探测工程技术有限公司 基于多目标函数的黏弹性参数同步反演方法、装置和设备
CN116626751B (zh) * 2023-05-26 2024-03-19 北京中矿大地地球探测工程技术有限公司 基于多目标函数的黏弹性参数同步反演方法、装置和设备
CN116699695A (zh) * 2023-08-07 2023-09-05 北京中矿大地地球探测工程技术有限公司 一种基于衰减矫正的反演方法、装置和设备
CN116699695B (zh) * 2023-08-07 2023-11-03 北京中矿大地地球探测工程技术有限公司 一种基于衰减矫正的反演方法、装置和设备

Also Published As

Publication number Publication date
CN113504566B (zh) 2024-04-30

Similar Documents

Publication Publication Date Title
Kim et al. 3-D traveltime computation using second-order ENO scheme
US9158017B2 (en) Seismic imaging apparatus utilizing macro-velocity model and method for the same
RU2577387C2 (ru) Скорость сходимости инверсии полного волнового поля при использовании формирования спектра
RU2587498C2 (ru) Инверсия одновременных источников для данных сейсмоприемной косы с взаимнокорреляционной целевой функцией
US9625593B2 (en) Seismic data processing
CA2972033C (en) Multistage full wavefield inversion process that generates a multiple free data set
CN110058303B (zh) 声波各向异性逆时偏移混合方法
CN113504566A (zh) 基于波动方程的地震反演方法、系统、装置及介质
CN111766628A (zh) 一种预条件的时间域弹性介质多参数全波形反演方法
CN110954945B (zh) 一种基于动态随机震源编码的全波形反演方法
GB2499096A (en) Simultaneous joint estimation of P-P and P-S residual statics
CN111077567B (zh) 一种基于矩阵乘法的双程波叠前深度偏移的方法
US9075160B2 (en) Inversion using a filtering operator
NO20190489A1 (en) Seismic modeling
CN111665556A (zh) 地层声波传播速度模型构建方法
Jia et al. Superwide-angle one-way wave propagator and its application in imaging steep salt flanks
CN111505714A (zh) 基于岩石物理约束的弹性波直接包络反演方法
CN111175822B (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
Ernst et al. Reduction of near-surface scattering effects in seismic data
CN111208568B (zh) 一种时间域多尺度全波形反演方法及系统
Shi et al. Multiscale full-waveform inversion based on shot subsampling
CN111665549A (zh) 地层声波衰减因子反演方法
Robertsson et al. Efficient snapshot-free reverse time migration and computation of multiparameter gradients in full-waveform inversion
Cho Efficient Seismic Depth Imaging and Full-Waveform Inversion via Generalized Multiscale Finite Element
CN117950045A (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