CN113504566A - 基于波动方程的地震反演方法、系统、装置及介质 - Google Patents
基于波动方程的地震反演方法、系统、装置及介质 Download PDFInfo
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 38
- 238000004364 calculation method Methods 0.000 claims abstract description 41
- 230000006870 function Effects 0.000 claims abstract description 25
- 230000035945 sensitivity Effects 0.000 claims abstract description 20
- 238000009499 grossing Methods 0.000 claims description 15
- 230000010354 integration Effects 0.000 claims description 10
- 230000008901 benefit Effects 0.000 abstract description 7
- 238000004422 calculation algorithm Methods 0.000 abstract description 6
- 238000003325 tomography Methods 0.000 abstract description 4
- 238000010586 diagram Methods 0.000 description 7
- 238000004587 chromatography analysis Methods 0.000 description 6
- 230000008569 process Effects 0.000 description 4
- 238000010521 absorption reaction Methods 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 241001344923 Aulorhynchidae Species 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000001902 propagating effect Effects 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
Images
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/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/282—Application of seismic models, synthetic seismograms
-
- 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/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
-
- 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/675—Wave 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、对波动方程正演中时间步的波场进行单频离散傅里叶积分得到合成地震记录和频率域入射波场;
具体地,实施例利用输入的地震子波、包含随机边界的速度模型及正演模拟参数进行波动方程正演。在本实施例中,具体采用标量声波方程正演模拟,并对每个时间步的波场进行单频离散傅里叶积分;正演结束后得到合成地震记录以及频率域入射波场。离散傅里叶积分公式如下:
在公式(1)中,u(x,kΔt;xs)为波动方程正演模拟得到的时间-空间域全波场,xs为震源位置,x为地下空间任一点的坐标,exp(-i2πfkΔt)为离散Fourier积分核,u(x,f;xs)为频率域入射波场,f为频率。合成地震记录定义为时间-空间域全波场在检波器位置xr处的空间采样,其表示为
需要补充说明的是,本实施例中,以二阶常密度标量声波方程的有限差分数值解为例,给出波动方程正演模拟及逆时传播对应的数学描述:
(1)波动方程正演模拟,需要求解如下初值问题:
(2)波动方程逆时传播,需要求解如下边值问题(伴随震源作为边值条件):
在一些可选的实施例中,根据合成地震记录确定走时伴随震源,根据走时伴随震源确定走时目标泛函的单频梯度这一步骤S400,其包括步骤S410-S430:
S410、根据走时伴随震源确定边界条件,通过波动方程逆时传播确定伴随波场;
S420、在波动方程逆时传播中对伴随波场进行单频离散傅里叶积分,得到频域伴随波场;
S430、根据频域伴随波场的入射波场和伴随波场,确定走时目标泛函的单频梯度;
具体地,实施例根据公式(2),通过合成地震记录计算走时伴随震源:
然后,实施例将走时伴随震源作为边界条件,进行波动方程逆时传播计算伴随波场。根据公式(5),在逆时传播过程中对伴随波场进行单频离散傅里叶积分,得到频率域伴随波场λ(x,f):
在公式(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为炮数,并进行求和得到走时目标泛函的单频梯度:
在一些其他的可选实施例中,从单频梯度提取低波数部分,根据低波数进行地震速度反演这一步骤S500,其具体为:通过高斯光滑滤波器,对单频梯度进行光滑处理,提取得到低波数部分。
具体地,实施例采用高斯光滑滤波器对步骤S430a中得到的单频梯度进行光滑处理,提取单频梯度中的低波数部分。实施例中所采用的高斯光滑滤波器如下:
在公式(8)中,σx,σy,σz分别为高斯光滑滤波器在x,y,z方向上的平滑半径。通过将高斯光滑滤波器与单频梯度卷积:G(x,σx,σy,σz)*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中任一项所述的基于波动方程的地震反演方法。
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)
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)
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约束的全波形反演方法 |
-
2021
- 2021-06-01 CN CN202110609064.1A patent/CN113504566B/zh active Active
Patent Citations (15)
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)
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)
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 |