CN111208568A - 一种时间域多尺度全波形反演方法及*** - Google Patents
一种时间域多尺度全波形反演方法及*** Download PDFInfo
- Publication number
- CN111208568A CN111208568A CN202010045893.7A CN202010045893A CN111208568A CN 111208568 A CN111208568 A CN 111208568A CN 202010045893 A CN202010045893 A CN 202010045893A CN 111208568 A CN111208568 A CN 111208568A
- Authority
- CN
- China
- Prior art keywords
- velocity model
- observation data
- model
- seismic observation
- iteration
- 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 80
- 238000003384 imaging method Methods 0.000 claims abstract description 19
- 238000001914 filtration Methods 0.000 claims abstract description 12
- 238000003325 tomography Methods 0.000 claims abstract description 7
- 230000006870 function Effects 0.000 claims description 93
- 238000004088 simulation Methods 0.000 claims description 20
- 238000010276 construction Methods 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 abstract description 13
- 230000036039 immunity Effects 0.000 abstract description 6
- 230000002708 enhancing effect Effects 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 17
- 230000008569 process Effects 0.000 description 10
- 239000011159 matrix material Substances 0.000 description 9
- 108010001267 Protein Subunits Proteins 0.000 description 6
- 230000005540 biological transmission Effects 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 208000010392 Bone Fractures Diseases 0.000 description 2
- 206010017076 Fracture Diseases 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 241001482237 Pica Species 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012795 verification Methods 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- 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/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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
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
本发明公开一种时间域多尺度全波形反演方法及***，所述方法包括：步骤S1：利用检波器获取地震观测数据；步骤S2：利用层析成像方法确定初始速度模型；步骤S3：利用维纳低通滤波器将所述地震观测数据滤波到不同频段；步骤S4：根据所述初始速度模型和不同频段的地震观测数据确定最优速度模型；步骤S5：基于所述最优速度模型进行地震成像，提高全波形反演的收敛速度和计算效率，增强反演的抗噪性。
Description
技术领域
本发明涉及地震勘探技术领域，特别是涉及一种时间域多尺度全波形反演方法及***。
背景技术
随着社会工业的快速发展，对资源需求不断增加，所以对厚覆盖层区、剧烈起伏地形(陡峭山地)、复杂地下地质结构(地下断裂发育、结构破碎、地层变形严重、岩性复杂)区进行深部资源探测需求越来越强烈。地震勘探是一个***工程，地震方法能否成功用于资源勘探，很大程度上取决于地震资料处理技术。地震成像是地震资料处理中的核心环节，是最终成果和效益的集中体现。随着计算机硬件的发展，以逆时偏移成像技术(RTM)为代表的地震成像技术已经得到了广泛的发展和应用。而以全波场信息拟合介质模型为手段的全波形反演(Full waveform inversion)技术可以为构造最为复杂地区成像提供高分辨率、高精度速度模型。
全波形反演从波动方程出发，综合利用叠前地震波场中的振幅、相位等信息，通过拟合实际波场和预测波场来修正初始假定的速度模型，以获取高精度的地下速度信息。全波形反演为深部大尺度构造演化分析、精确的地震成像及速度建模等方面提供可靠依据，因此，其作为一种高精度速度模型构建方法受到重视，是目前国际地球物理学界的研究热点。全波形反演是一个高度非线性的问题，对实际观测数据的质量(偏移距大小、有效低频以及噪音情况)、反演过程中采用的技术手段(反演算法、步长搜索方法、目标函数、反演策略)和初始输入模型的准确度都比较敏感。同时，这些因素都不同程度地影响着全波形反演在实际地震资料中的应用。
全波形反演的主要目标是最小化观测数据与合成地震数据构成的目标函数，通过伴随算子计算梯度，利用一定的搜索方法计算步长，结合合适的反演算法构造搜索方向，不断地对模型参数进行更新，最终获得高分辨率参数模型。目标函数和步长搜索方法是全波形反演中至关重要的两个因素。
目标函数决定着全波形反演的非线性程度和抗噪性。传统的全波形反演的目标函数是基于计算数据和观测数据残差的L2范数，但是L2范数存在对非高斯噪音敏感的问题，在面对实际数据中存在的各种强噪声，全波形反演容易陷入局部极小值。
优化步长搜索方法可以加快全波形反演的收敛到极小值的速度，并且涉及很少的额外正演，提高全波形反演的计算效率，但现有Pica等(1990)提出一种从目标函数出发，直接求解最优步长的方法，但这种方法仅适用于抗噪性较差的L2范数。Vigh等(2009)提出使用抛物线方法，利用三个满足关系的步长值和其对应的目标函数值求解抛物线系数，获得最优步长值，但此方法至少需要额外两次正演，降低了计算效率。非精确线性搜索方法保证函数估计次数较小，但是搜索的过程通常需要多次额外的正演计算，计算成本较高。
发明内容
基于此，本发明的目的是提供一种时间域多尺度全波形反演方法及***，以提高全波形反演的抗噪性、收敛速度和计算效率。
为实现上述目的，本发明提供一种时间域多尺度全波形反演方法，所述方法包括：
步骤S1：利用检波器获取地震观测数据；
步骤S2：利用层析成像方法确定初始速度模型；
步骤S3：利用维纳低通滤波器将所述地震观测数据滤波到不同频段，分别为低频带地震观测数据、中频带地震观测数据和高频带地震观测数据；
步骤S4：根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型；
步骤S5：基于所述最优速度模型进行地震成像。
可选的，所述根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型，包括：
步骤S41：根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型；
步骤S42：根据所述中频带地震观测数据和所述第一速度模型进行反演，获得第二速度模型；
步骤S43：根据所述高频带地震观测数据和所述第二速度模型进行反演，获得最优速度模型。
可选的，所述根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型，包括：
步骤S411：计算所述初始速度模型对应的模拟波场数据；
步骤S412：基于所述模拟波场数据和所述地震观测数据构建目标函数；
步骤S413：采用伴随状态法计算所述目标函数的梯度；
步骤S414：根据所述梯度确定所述目标函数的搜索方向；
步骤S415：确定所述目标函数的最优步长；
步骤S416：根据所述搜索方向和所述最优步长确定下次迭代的速度模型；
步骤S417：判断迭代次数是否大于或等于第一迭代次数阈值；如果迭代次数小于第一迭代次数阈值，则将下次迭代的速度模型作为所述初始速度模型，并返回“步骤S411”；如果迭代次数大于或等于第一迭代次数阈值，则将下次迭代的速度模型作为第一速度模型，并执行“步骤S42”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于第一设定系数，则将下次迭代的速度模型作为所述初始速度模型，并返回“步骤S411”；如果所述连续目标函数相对差值小于第一设定系数，则将下次迭代的速度模型作为第一速度模型，并执行“步骤S42”。
可选的，所述计算所述初始速度模型对应的模拟波场数据，包括：
采用有限差分法，利用声波方程计算所述初始速度模型对应的模拟波场数据。
本发明还提供一种时间域多尺度全波形反演***，所述***包括：
地震观测数据获取模块，用于利用检波器获取地震观测数据；
初始速度模型确定模块，用于利用层析成像方法确定初始速度模型；
滤波模块，用于利用维纳低通滤波器将所述地震观测数据滤波到不同频段，分别为低频带地震观测数据、中频带地震观测数据和高频带地震观测数据；
最优速度模型确定模块，用于根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型；
地震成像模块，用于基于所述最优速度模型进行地震成像。
可选的，所述最优速度模型确定模块，包括：
第一速度模型确定单元，用于根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型；
第二速度模型确定单元，用于根据所述中频带地震观测数据和所述第一速度模型进行反演，获得第二速度模型；
最优速度模型确定单元，用于根据所述高频带地震观测数据和所述第二速度模型进行反演，获得最优速度模型。
可选的，所述第一速度模型确定单元，包括：
第一模拟波场数据确定子单元，用于计算所述初始速度模型对应的模拟波场数据；
目标函数构建子单元，用于基于所述模拟波场数据和所述地震观测数据构建目标函数；
梯度确定子单元，用于采用伴随状态法计算所述目标函数的梯度；
搜索方向确定子单元，用于根据所述梯度确定所述目标函数的搜索方向；
最优步长确定子单元，用于确定所述目标函数的最优步长；
下次迭代的速度模型确定子单元，用于根据所述搜索方向和所述最优步长确定下次迭代的速度模型；
第一判断子单元，用于判断迭代次数是否大于或等于第一迭代次数阈值；如果迭代次数小于第一迭代次数阈值，则将下次迭代的速度模型作为所述初始速度模型，并返回“第一模拟波场数据确定子单元”；如果迭代次数大于或等于第一迭代次数阈值，则将下次迭代的速度模型作为第一速度模型，并执行“第二速度模型确定单元”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于第一设定系数，则将下次迭代的速度模型作为所述初始速度模型，并返回“第一模拟波场数据确定子单元”；如果所述连续目标函数相对差值小于第一设定系数，则将下次迭代的速度模型作为第一速度模型，并执行“第二速度模型确定单元”。
可选的，所述第一模拟波场数据确定子单元，具体包括：
用于采用有限差分法，利用声波方程计算所述初始速度模型对应的模拟波场数据。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明公开一种时间域多尺度全波形反演方法及***，所述方法包括：步骤S1：利用检波器获取地震观测数据；步骤S2：利用层析成像方法确定初始速度模型；步骤S3：利用维纳低通滤波器将所述地震观测数据滤波到不同频段；步骤S4：根据所述初始速度模型和不同频段的地震观测数据确定最优速度模型；步骤S5：基于所述最优速度模型进行地震成像，提高全波形反演的收敛速度和计算效率，增强反演的抗噪性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例时间域多尺度全波形反演方法流程图；
图2为本发明实施例初始速度模型示意图；
图3为本发明实施例反演获得的最优速度模型；
图4为本发明实施例反演获得的最优速度模型深度剖面分别在水平位置示意图；
图5为本发明实施例反演过程的目标函数收敛曲线图；
图6为本发明实施例地震观测数据分布示意图；
图7为本发明实施例加入噪音后反演获得的最优速度模型；
图8为本发明实施例加入噪音后反演获得的最优速度模型深度剖面分别在水平位置示意图；
图9为本发明实施例加入噪音后反演过程的目标函数收敛曲线图；
图10为本发明实施例时间域多尺度全波形反演***结构图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种时间域多尺度全波形反演方法及***，以提高全波形反演的抗噪性、收敛速度和计算效率。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明实施例时间域多尺度全波形反演方法流程图，如图1所示，本发明提供一种时间域多尺度全波形反演方法，所述方法包括：
步骤S1：利用检波器获取地震观测数据。
步骤S2：利用层析成像方法确定初始速度模型。
步骤S3：利用维纳低通滤波器将所述地震观测数据滤波到不同频段，分别为低频带地震观测数据、中频带地震观测数据和高频带地震观测数据。
步骤S4：根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型。
步骤S5：基于所述最优速度模型进行地震成像。
下面对各个步骤进行详细论述：
步骤S1：利用检波器获取地震观测数据。
在野外激发震源，利用检波器获取地震观测数据，所述地震观测数据既可以是加入噪声的地震观测数据，也可以不加入噪声的地震观测数据。
步骤S2：利用层析成像方法确定初始速度模型。
全波形反演中数据残差和模型参数的强非线性关系，导致反演过程容易陷入局部极小值，可见初始速度模型直接决定反演成功与否，因此本申请利用层析成像方法确定初始速度模型，防止在反演过程中陷入局部极小值的现象。
步骤S3：利用维纳低通滤波器将所述地震观测数据滤波到不同频段，分别为低频带地震观测数据、中频带地震观测数据和高频带地震观测数据。
本发明中各个频带的范围是根据实际数据的频谱分布来确定，不同地区或者不同数据的低频带和高频带会差别，没有统一标准。
本发明低频带地震观测数据可以构造出宏观构造的速度模型，然后将低频带地震观测数据输入到中频带反演流程后再输入到高频带反演流程，随着反演的频带不断提高，可以对模型进行精细刻画，进而获得精确度较高速度模型。
具体步骤如下：
步骤S4：所述根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型，包括：
步骤S41：根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型。
步骤S42：根据所述中频带地震观测数据和所述第一速度模型进行反演，获得第二速度模型。
步骤S43：根据所述高频带地震观测数据和所述第二速度模型进行反演，获得最优速度模型。
其中，步骤S41：根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型，包括：
步骤S411：计算所述初始速度模型对应的模拟波场数据。
本发明使用有限差分方法求解方程(1)，其是一种广泛使用的数值模拟方法，利用局部算子，求解简单、容易编程实现且所需内存小，计算效率高。本发明方法所涉及的震源正传、残差反传均采用时间2阶、空间10阶的有限差分格式进行，数值模拟中采用完全匹配层吸收边界条件(PML)对人为边界产生的强烈反射进行吸收，完成震源正传，具体步骤如下：
采用有限差分法，利用声波方程计算所述初始速度模型对应的模拟波场数据。
所述声波方程为：
其中，t为波场传播时间，x、z分别为波场的水平坐标和垂直坐标，
为第k次迭代的模拟波场数据，vk为第k次迭代的速度模型，s(x,z,t)为震源项。
步骤S412：基于所述模拟波场数据和所述地震观测数据构建目标函数；所述目标函数为L1范数型目标函数；所述目标函数为：
为方便计算，在公式(2)中令：
其中，s为炮数，l为检波器数，
为第k次迭代的模拟波场数据，pk obs(x,z,t)为第k次迭代的地震观测数据。
步骤S413：采用伴随状态法计算所述目标函数的梯度；所述梯度公式为：
其中，gk为第k次迭代的梯度，T为最大计算时间，λ为检波器处残差反传波场，δpk为第k次迭代的反传残差震源。
拟牛顿算法的特点是不需要直接计算海森矩阵的逆(计算成本较高)，而是采用的一定的手段来近似构造海森矩阵的逆矩阵。有限内存牛顿法(L-BFGS)避免了存储近似海森矩阵逆矩阵，而是通过存储有限个梯度变化量和速度模型更新量(3～20)构造搜索方向，因此公开了以下步骤：
步骤S414：根据所述梯度确定所述目标函数的搜索方向；具体的，采用双循环迭代法，根据所述梯度确定所述目标函数的搜索方向；所述搜索方向公式为：
dk＝-Hkgk (5)
其中，gk为第k次迭代的梯度，dk为第k次迭代的搜索方向，Hk为第k次迭代中海森矩阵的逆矩阵，Hk的表达式如下：
其中，
I为单位矩阵，
yk为第k次迭代的梯度变化量，yk＝gk+1-gk，sk为第k次迭代的速度更新量，
为近似海森逆矩阵的初始矩阵，m为存储的梯度变化量的数量，通常为3～20。
步骤S415：确定所述目标函数的最优步长。
对模拟波场数据进行泰勒展开，可以得到：
其中，pcal(vk)为模拟观测波场，αk为第k次迭代的最优步长，vk代表第k次迭代的速度模型，即最优速度模型，
为拉普拉斯算子。
因此目标函数(2)可以写为：
其中，pobs(vk)为实际观测波场；
为了求得最优目标函数，使目标函数对变量步长进行求导，并使导数等于0，即：
通过上面的一系列公式计算和推导，可以得到针对L1范数型目标函数的最优步长计算公式为：
其中，αt为测试步长，c1为一个权重数值，
步骤S416：根据所述搜索方向和所述最优步长确定下次迭代的速度模型，具体公式为：
vk+1＝vk+αkdk (10)
其中，vk+1为第k+1次迭代的速度模型，vk为第k次迭代的速度模型，αk为第k次迭代的最优步长，dk为第k次迭代的搜索方向。
步骤S417：判断迭代次数是否大于或等于第一迭代次数阈值；如果迭代次数小于第一迭代次数阈值，则将下次迭代的速度模型作为所述初始速度模型，并返回“步骤S411”；如果迭代次数大于或等于第一迭代次数阈值，则将下次迭代的速度模型作为第一速度模型，并执行“步骤S42”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于第一设定系数，则将下次迭代的速度模型作为所述初始速度模型，并返回“步骤S411”；如果所述连续目标函数相对差值小于第一设定系数，则将下次迭代的速度模型作为第一速度模型，并执行“步骤S42”。
步骤S42：根据所述中频带地震观测数据和所述第一速度模型进行反演，获得第二速度模型。
步骤S421：计算所述第一速度模型对应的模拟波场数据，具体包括：
采用有限差分法，利用声波方程计算所述第一速度模型对应的模拟波场数据。
步骤S422～步骤S426分别与步骤S412～步骤S416的步骤完全相同，在此不再一一赘述。
步骤S427：判断迭代次数是否大于或等于第二迭代次数阈值；如果迭代次数小于第二迭代次数阈值，则将下次迭代的速度模型作为所述第一速度模型，并返回“步骤S421”；如果迭代次数大于或等于第二迭代次数阈值，则将下次迭代的速度模型作为第二速度模型，并执行“步骤S43”。
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第二设定系数；如果所述连续目标函数相对差值大于或等于第二设定系数，则将下次迭代的速度模型作为所述第一速度模型，并返回“步骤S421”；如果所述连续目标函数相对差值小于第二设定系数，则将下次迭代的速度模型作为第二速度模型，并执行“步骤S43”。
步骤S43：根据所述高频带地震观测数据和所述第二速度模型进行反演，获得最优速度模型，包括：
步骤S431：计算所述第二速度模型对应的模拟波场数据，具体包括：
采用有限差分法，利用声波方程计算所述第二速度模型对应的模拟波场数据。
步骤S432～步骤S436分别与步骤S412～步骤S416的步骤完全相同，在此不再一一赘述。
步骤S437：判断迭代次数是否大于或等于第三迭代次数阈值；如果迭代次数小于第三迭代次数阈值，则将下次迭代的速度模型作为所述第二速度模型，并返回“步骤S431”；如果迭代次数大于或等于第三迭代次数阈值，则将下次迭代的速度模型作为最优速度模型，并执行“步骤S5”。
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第三设定系数；如果所述连续目标函数相对差值大于或等于第三设定系数，则将下次迭代的速度模型作为所述第二速度模型，并返回“步骤S431”；如果所述连续目标函数相对差值小于第三设定系数，则将下次迭代的速度模型作为最优速度模型，并执行“步骤S5”。
图10为本发明实施例时间域多尺度全波形反演***结构图，如图10所示，本发明还提供一种时间域多尺度全波形反演***，所述***包括：
地震观测数据获取模块1，用于利用检波器获取地震观测数据；
初始速度模型确定模块2，用于利用层析成像方法确定初始速度模型；
滤波模块3，用于利用维纳低通滤波器将所述地震观测数据滤波到不同频段，分别为低频带地震观测数据、中频带地震观测数据和高频带地震观测数据；
最优速度模型确定模块4，用于根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型；
地震成像模块5，用于基于所述最优速度模型进行地震成像。
作为一种实施方式，本发明所述最优速度模型确定模块4，包括：
第一速度模型确定单元，用于根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型。
第二速度模型确定单元，用于根据所述中频带地震观测数据和所述第一速度模型进行反演，获得第二速度模型。
最优速度模型确定单元，用于根据所述高频带地震观测数据和所述第二速度模型进行反演，获得最优速度模型。
作为一种实施方式，本发明所述第一速度模型确定单元，包括：
第一模拟波场数据确定子单元，用于计算所述初始速度模型对应的模拟波场数据；
目标函数构建子单元，用于基于所述模拟波场数据和所述地震观测数据构建目标函数；
梯度确定子单元，用于采用伴随状态法计算所述目标函数的梯度；
搜索方向确定子单元，用于根据所述梯度确定所述目标函数的搜索方向；
最优步长确定子单元，用于确定所述目标函数的最优步长；
下次迭代的速度模型确定子单元，用于根据所述搜索方向和所述最优步长确定下次迭代的速度模型；
第一判断子单元，用于判断迭代次数是否大于或等于第一迭代次数阈值；如果迭代次数小于第一迭代次数阈值，则将下次迭代的速度模型作为所述初始速度模型，并返回“第一模拟波场数据确定子单元”；如果迭代次数大于或等于第一迭代次数阈值，则将下次迭代的速度模型作为第一速度模型，并执行“第二速度模型确定单元”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于第一设定系数，则将下次迭代的速度模型作为所述初始速度模型，并返回“第一模拟波场数据确定子单元”；如果所述连续目标函数相对差值小于第一设定系数，则将下次迭代的速度模型作为第一速度模型，并执行“第二速度模型确定单元”。
作为一种实施方式，本发明所述第一模拟波场数据确定子单元，具体包括：
采用有限差分法，利用声波方程计算所述初始速度模型对应的模拟波场数据。
作为一种实施方式，本发明所述第二速度模型确定单元，包括：
第二模拟波场数据确定子单元，用于计算所述第一速度模型对应的模拟波场数据；
第二判断子单元，用于判断迭代次数是否大于或等于第二迭代次数阈值；如果迭代次数小于第二迭代次数阈值，则将下次迭代的速度模型作为所述第一速度模型，并返回“第二模拟波场数据确定子单元”；如果迭代次数大于或等于第二迭代次数阈值，则将下次迭代的速度模型作为第二速度模型，并执行“第三速度模型确定单元”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第二设定系数；如果所述连续目标函数相对差值大于或等于第二设定系数，则将下次迭代的速度模型作为所述第一速度模型，并返回“第二模拟波场数据确定子单元”；如果所述连续目标函数相对差值小于第二设定系数，则将下次迭代的速度模型作为第二速度模型，并执行“第三速度模型确定单元”。
所述第二速度模型确定单元中的剩下的子单元与所述第一速度模型确定单元中目标函数构建子单元、梯度确定子单元、搜索方向确定子单元、最优步长确定子单元和下次迭代的速度模型确定子单元完全相同，在此不再一一赘述。
作为一种实施方式，本发明所述第二模拟波场数据确定子单元，具体包括：
采用有限差分法，利用声波方程计算所述第一速度模型对应的模拟波场数据。
作为一种实施方式，本发明所述最优速度模型确定单元，包括：
第三模拟波场数据确定子单元，用于计算所述第二速度模型对应的模拟波场数据；
第三判断子单元，用于判断迭代次数是否大于或等于第三迭代次数阈值；如果迭代次数小于第三迭代次数阈值，则将下次迭代的速度模型作为所述第二速度模型，并返回“第三模拟波场数据确定子单元”；如果迭代次数大于或等于第三迭代次数阈值，则将下次迭代的速度模型作为最优速度模型，并执行“地震成像模块”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第三设定系数；如果所述连续目标函数相对差值大于或等于第三设定系数，则将下次迭代的速度模型作为所述第二速度模型，并返回“第三模拟波场数据确定子单元”；如果所述连续目标函数相对差值小于第三设定系数，则将下次迭代的速度模型作为最优速度模型，并执行“地震成像模块”。
所述最优速度模型确定单元中的剩下的子单元与所述第一速度模型确定单元中目标函数构建子单元、梯度确定子单元、搜索方向确定子单元、最优步长确定子单元和下次迭代的速度模型确定子单元完全，在此不再一一赘述。
作为一种实施方式，本发明所述第三模拟波场数据确定子单元，具体包括：
采用有限差分法，利用声波方程计算所述第二速度模型对应的模拟波场数据。
采用本发明公开的时间域多尺度全波形反演方法及***进行试验验证，图2为本发明实施例初始速度模型示意图；图2(a)和图2(b)分别是Overthrust真实模型和光滑初始模型；图3为本发明实施例反演获得的最优速度模型，可以发现本发明反演获得的速度模型结构清晰，推覆体和和层状模型都很清楚；图4为本发明实施例反演获得的最优速度模型任意抽取的三条深度剖面分别在水平位置示意图；(a)为反演获得的最优速度模型深度剖面在水平位置3.0km处的示意图，(b)为反演获得的最优速度模型深度剖面在水平位置4.8km处的示意图，(c)为反演获得的最优速度模型深度剖面在水平位置6.0km处的示意图，可以发现本发明方法获得的深度剖面与真实的深度剖面匹配度很好；图5为本发明实施例反演过程的目标函数收敛曲线图；图6为本发明实施例地震观测数据分布示意图，可以发现本发明对应的目标函数收敛的很快，在较少的迭代次数情况下，反演误差小，精度高；图6中(a)和(b)分别为Overthust模型的第30炮的地震观测数据和加入强噪音的地震观测数据分布示意图；图7为本发明实施例加入噪音后反演获得的最优速度模型，可以发现在数据中掺有强噪音的情况下，虽然与无噪音情况相比，反演精度下降，但是本发明依旧可以获得高精度的速度模型；图8为本发明实施例加入噪音后反演获得的最优速度模型任意抽取的三条深度剖面分别在水平位置示意图；(a)为加入噪音后反演获得的最优速度模型深度剖面在水平位置3.0km处的示意图，(b)为加入噪音后反演获得的最优速度模型深度剖面在水平位置4.8km处的示意图，(c)为加入噪音后反演获得的最优速度模型深度剖面在水平位置6.0km处的示意图，可以发现，和真实速度模型整体匹配度很好，深层和浅层的速度都得到了很好的恢复重建；图9为本发明实施例加入噪音后反演过程的目标函数收敛曲线图，可以发现，在加入强噪音的情况下，目标函数曲线在保证精度的前提下，收敛速度较快；图2-图9证明了本发明方法的稳健性。
本发明提出了一种时间域多尺度全波形反演方法及***，在全波形反演中使用L1范数型目标函数，提高反演的抗噪性，在地震数据含有强噪音的情况下，依然可以获得满意的反演结果。针对L1范数型目标函数，提出一种适用的最优步长搜索方法，从而提高全波形反演的收敛速度和计算效率，增强反演的抗噪性。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
Claims (8)
1.一种时间域多尺度全波形反演方法，其特征在于，所述方法包括：
步骤S1：利用检波器获取地震观测数据；
步骤S2：利用层析成像方法确定初始速度模型；
步骤S3：利用维纳低通滤波器将所述地震观测数据滤波到不同频段，分别为低频带地震观测数据、中频带地震观测数据和高频带地震观测数据；
步骤S4：根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型；
步骤S5：基于所述最优速度模型进行地震成像。
2.根据权利要求1所述的时间域多尺度全波形反演方法，其特征在于，所述根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型，包括：
步骤S41：根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型；
步骤S42：根据所述中频带地震观测数据和所述第一速度模型进行反演，获得第二速度模型；
步骤S43：根据所述高频带地震观测数据和所述第二速度模型进行反演，获得最优速度模型。
3.根据权利要求2所述的时间域多尺度全波形反演方法，其特征在于，所述根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型，包括：
步骤S411：计算所述初始速度模型对应的模拟波场数据；
步骤S412：基于所述模拟波场数据和所述地震观测数据构建目标函数；
步骤S413：采用伴随状态法计算所述目标函数的梯度；
步骤S414：根据所述梯度确定所述目标函数的搜索方向；
步骤S415：确定所述目标函数的最优步长；
步骤S416：根据所述搜索方向和所述最优步长确定下次迭代的速度模型；
步骤S417：判断迭代次数是否大于或等于第一迭代次数阈值；如果迭代次数小于第一迭代次数阈值，则将下次迭代的速度模型作为所述初始速度模型，并返回“步骤S411”；如果迭代次数大于或等于第一迭代次数阈值，则将下次迭代的速度模型作为第一速度模型，并执行“步骤S42”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于第一设定系数，则将下次迭代的速度模型作为所述初始速度模型，并返回“步骤S411”；如果所述连续目标函数相对差值小于第一设定系数，则将下次迭代的速度模型作为第一速度模型，并执行“步骤S42”。
4.根据权利要求3所述的时间域多尺度全波形反演方法，其特征在于，所述计算所述初始速度模型对应的模拟波场数据，包括：
采用有限差分法，利用声波方程计算所述初始速度模型对应的模拟波场数据。
5.一种时间域多尺度全波形反演***，其特征在于，所述***包括：
地震观测数据获取模块，用于利用检波器获取地震观测数据；
初始速度模型确定模块，用于利用层析成像方法确定初始速度模型；
滤波模块，用于利用维纳低通滤波器将所述地震观测数据滤波到不同频段，分别为低频带地震观测数据、中频带地震观测数据和高频带地震观测数据；
最优速度模型确定模块，用于根据所述初始速度模型、低频带地震观测数据、中频带地震观测数据和高频带地震观测数据确定最优速度模型；
地震成像模块，用于基于所述最优速度模型进行地震成像。
6.根据权利要求5所述的时间域多尺度全波形反演***，其特征在于，所述最优速度模型确定模块，包括：
第一速度模型确定单元，用于根据所述低频带地震观测数据和所述初始速度模型进行反演，获得第一速度模型；
第二速度模型确定单元，用于根据所述中频带地震观测数据和所述第一速度模型进行反演，获得第二速度模型；
最优速度模型确定单元，用于根据所述高频带地震观测数据和所述第二速度模型进行反演，获得最优速度模型。
7.根据权利要求5所述的时间域多尺度全波形反演***，其特征在于，所述第一速度模型确定单元，包括：
第一模拟波场数据确定子单元，用于计算所述初始速度模型对应的模拟波场数据；
目标函数构建子单元，用于基于所述模拟波场数据和所述地震观测数据构建目标函数；
梯度确定子单元，用于采用伴随状态法计算所述目标函数的梯度；
搜索方向确定子单元，用于根据所述梯度确定所述目标函数的搜索方向；
最优步长确定子单元，用于确定所述目标函数的最优步长；
下次迭代的速度模型确定子单元，用于根据所述搜索方向和所述最优步长确定下次迭代的速度模型；
第一判断子单元，用于判断迭代次数是否大于或等于第一迭代次数阈值；如果迭代次数小于第一迭代次数阈值，则将下次迭代的速度模型作为所述初始速度模型，并返回“第一模拟波场数据确定子单元”；如果迭代次数大于或等于第一迭代次数阈值，则将下次迭代的速度模型作为第一速度模型，并执行“第二速度模型确定单元”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于第一设定系数，则将下次迭代的速度模型作为所述初始速度模型，并返回“第一模拟波场数据确定子单元”；如果所述连续目标函数相对差值小于第一设定系数，则将下次迭代的速度模型作为第一速度模型，并执行“第二速度模型确定单元”。
8.根据权利要求7所述的时间域多尺度全波形反演***，其特征在于，所述第一模拟波场数据确定子单元，包括：
用于采用有限差分法，利用声波方程计算所述初始速度模型对应的模拟波场数据。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010045893.7A CN111208568B (zh) | 2020-01-16 | 2020-01-16 | 一种时间域多尺度全波形反演方法及*** |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010045893.7A CN111208568B (zh) | 2020-01-16 | 2020-01-16 | 一种时间域多尺度全波形反演方法及*** |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111208568A true CN111208568A (zh) | 2020-05-29 |
| CN111208568B CN111208568B (zh) | 2021-04-09 |
Family
ID=70787266
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010045893.7A Expired - Fee Related CN111208568B (zh) | 2020-01-16 | 2020-01-16 | 一种时间域多尺度全波形反演方法及*** |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111208568B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112379415A (zh) * | 2020-11-06 | 2021-02-19 | 中国科学技术大学 | 基于减采样的重构低频数据多尺度全波形反演方法及装置 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105005076A (zh) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | 基于最小二乘梯度更新速度模型的地震波全波形反演方法 |
| CN105319581A (zh) * | 2014-07-31 | 2016-02-10 | 中国石油化工股份有限公司 | 一种高效的时间域全波形反演方法 |
| CN105549079A (zh) * | 2016-01-12 | 2016-05-04 | 中国矿业大学(北京) | 一种地球物理参数的全波形反演模型的建立方法和装置 |
| WO2016165064A1 (zh) * | 2015-04-14 | 2016-10-20 | 中国科学院自动化研究所 | 基于多视角学习的鲁棒性前景检测方法 |
| CN106501852A (zh) * | 2016-10-21 | 2017-03-15 | 中国科学院地质与地球物理研究所 | 一种三维声波方程任意域多尺度全波形反演方法及装置 |
| US20180052245A1 (en) * | 2016-08-17 | 2018-02-22 | Pgs Geophysical As | Marine Vibrator Source Acceleration and Pressure |
| CN107765302A (zh) * | 2017-10-20 | 2018-03-06 | 吉林大学 | 不依赖震源子波的时间域单频波形走时反演方法 |
| CN107894618A (zh) * | 2017-11-10 | 2018-04-10 | 中国海洋大学 | 一种基于模型平滑算法的全波形反演梯度预处理方法 |
-
2020
- 2020-01-16 CN CN202010045893.7A patent/CN111208568B/zh not_active Expired - Fee Related
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105319581A (zh) * | 2014-07-31 | 2016-02-10 | 中国石油化工股份有限公司 | 一种高效的时间域全波形反演方法 |
| WO2016165064A1 (zh) * | 2015-04-14 | 2016-10-20 | 中国科学院自动化研究所 | 基于多视角学习的鲁棒性前景检测方法 |
| CN105005076A (zh) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | 基于最小二乘梯度更新速度模型的地震波全波形反演方法 |
| CN105549079A (zh) * | 2016-01-12 | 2016-05-04 | 中国矿业大学(北京) | 一种地球物理参数的全波形反演模型的建立方法和装置 |
| US20180052245A1 (en) * | 2016-08-17 | 2018-02-22 | Pgs Geophysical As | Marine Vibrator Source Acceleration and Pressure |
| CN106501852A (zh) * | 2016-10-21 | 2017-03-15 | 中国科学院地质与地球物理研究所 | 一种三维声波方程任意域多尺度全波形反演方法及装置 |
| CN107765302A (zh) * | 2017-10-20 | 2018-03-06 | 吉林大学 | 不依赖震源子波的时间域单频波形走时反演方法 |
| CN107894618A (zh) * | 2017-11-10 | 2018-04-10 | 中国海洋大学 | 一种基于模型平滑算法的全波形反演梯度预处理方法 |
Non-Patent Citations (3)
| Title |
|---|
| 崔永福 等: "全波形反演在缝洞型储层速度建模中的应用", 《地球物理学报》 * |
| 成景旺 等: "基于三次卷积插值的时间域多尺度全波形反演", 《煤炭学报》 * |
| 苗永康: "基于L-BFGS算法的时间域全波形反演", 《石油地球物理勘探》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112379415A (zh) * | 2020-11-06 | 2021-02-19 | 中国科学技术大学 | 基于减采样的重构低频数据多尺度全波形反演方法及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111208568B (zh) | 2021-04-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106526674B (zh) | 一种三维全波形反演能量加权梯度预处理方法 | |
| CA2726462C (en) | Estimation of soil properties using waveforms of seismic surface waves | |
| KR101549388B1 (ko) | 탄성파 다성분 자료에 대한 중합전 ｅｇｓ 구조보정 방법 | |
| Vigh et al. | Developing earth models with full waveform inversion | |
| CN107843925B (zh) | 一种基于修正相位的反射波波形反演方法 | |
| US8990053B2 (en) | Method of wavelet estimation and multiple prediction in full wavefield inversion | |
| Raknes et al. | Time-lapse full-waveform inversion of limited-offset seismic data using a local migration regularization | |
| Li et al. | Elastic reflection waveform inversion with variable density | |
| CN113740901B (zh) | 基于复杂起伏地表的陆上地震数据全波形反演方法及装置 | |
| EA032186B1 (ru) | Сейсмическая адаптивная фокусировка | |
| KR20170118185A (ko) | 다중반사파 없는 데이터 세트를 생성하는 다단식 전 파동장 역산 프로세스 | |
| CN111580163B (zh) | 一种基于非单调搜索技术的全波形反演方法及*** | |
| CN109655890B (zh) | 一种深度域浅中深层联合层析反演速度建模方法及*** | |
| CN109946742A (zh) | 一种TTI介质中纯qP波地震数据模拟方法 | |
| CN104749625B (zh) | 一种基于正则化技术的地震数据倾角估计方法及装置 | |
| Hu et al. | Ray-illumination compensation for adjoint-state first-arrival traveltime tomography | |
| Gao et al. | An efficient vector elastic reverse time migration method in the hybrid time and frequency domain for anisotropic media | |
| CN111208568B (zh) | 一种时间域多尺度全波形反演方法及*** | |
| Li et al. | Waveform inversion of seismic first arrivals acquired on irregular surface | |
| Gong et al. | Combined migration velocity model-building and its application in tunnel seismic prediction | |
| CN111435174B (zh) | 强反射地区地震资料振幅补偿方法及装置 | |
| CN108680957B (zh) | 基于加权的局部互相关时频域相位反演方法 | |
| CN114814944B (zh) | 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 | |
| Huang et al. | Data-assimilated time-lapse visco-acoustic full-waveform inversion: Theory and application for injected CO2 plume monitoring | |
| Wang et al. | Rock Fracture Monitoring Based on High‐Precision Microseismic Event Location Using 3D Multiscale Waveform Inversion |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210409 Termination date: 20220116 |