CN109407152A - 基于零均值归一化互相关目标函数的时间域全波形反演方法 - Google Patents
基于零均值归一化互相关目标函数的时间域全波形反演方法 Download PDFInfo
- Publication number
- CN109407152A CN109407152A CN201811550994.9A CN201811550994A CN109407152A CN 109407152 A CN109407152 A CN 109407152A CN 201811550994 A CN201811550994 A CN 201811550994A CN 109407152 A CN109407152 A CN 109407152A
- Authority
- CN
- China
- Prior art keywords
- model
- low
- objective function
- record
- full waveform
- 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 47
- 238000001914 filtration Methods 0.000 claims description 14
- 230000008569 process Effects 0.000 claims description 12
- 238000005457 optimization Methods 0.000 claims description 10
- 238000005070 sampling Methods 0.000 claims description 10
- 238000004422 calculation algorithm Methods 0.000 claims description 9
- 238000009795 derivation Methods 0.000 claims description 4
- 238000007781 pre-processing Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 2
- 238000001514 detection method Methods 0.000 claims 1
- 238000004088 simulation Methods 0.000 abstract description 7
- 238000010606 normalization Methods 0.000 abstract description 5
- 230000006870 function Effects 0.000 description 46
- 239000011159 matrix material Substances 0.000 description 12
- 239000013598 vector Substances 0.000 description 8
- 230000002159 abnormal effect Effects 0.000 description 5
- 238000011161 development Methods 0.000 description 5
- 235000013399 edible fruits Nutrition 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 230000005284 excitation Effects 0.000 description 3
- 238000004513 sizing Methods 0.000 description 3
- 238000009825 accumulation Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000011109 contamination Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 230000002969 morbid Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/364—Seismic filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
- G01V2210/324—Filtering
-
- 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
本发明涉及一种基于零均值归一化互相关目标函数的时间域全波形反演方法,通过计算模拟和观测记录之间的零均值归一化互相关代替传统全波形反演中最小二乘目标函数的求取。通过归一化互相关使模拟数据与观测数据振幅的绝对值在‑1到1之间,有效的减小了因振幅误差导致的反演结果错误。且由于振幅能量的归一化,提高了地震记录中远偏移距的能量权重,因远偏移距数据携带了模型中大尺度构造的信息,该方法能减小跳周的发生。因实际地震记录易受低频噪声的影响,且通常为非均值噪声,引入了去均值项,通过对模拟和观测记录的去均值计算减小了低频噪声的干扰。该方法能在缺失低频信息或者低频信息被噪声污染的情况下准确的反演出地下介质的速度参数。
Description
技术领域
本发明属于地震勘探技术领域,涉及一种地震勘探领域的地下介质成像和介质速度计算方法,具体涉及一种时间域全波形反演方法,特别涉及一种基于零均值归一化互相关目标函数的时间域全波形反演方法。
背景技术
石油矿产等不可再生资源在人类生活的方方面面起着至关重要的作用。因此,相应的资源勘探技术也在不断革新和发展。近年来,随着陆地和海洋油气开发的规模越来越大,勘探目标由浅层简单构造的地质区域向深部复杂构造区域以及隐藏油气藏变化,勘探难度日益增大。为了进一步解决深部能源的勘探与开采问题,地震勘探领域的全波形反演方法快速发展起来。全波形反演方法利用了叠前地震资料的全部信息包括走时、振幅、相位等,使模拟的地震记录逐渐拟合实际勘探采集到的观测记录,最终将地下介质的速度、密度、各向异性参数等特征信息计算出来。全波形反演以其对地下介质参数计算的准确性之高,吸引了全球范围内众多地球物理学家的广泛研究。
全波形反演方法最初在20世纪80年代由Tarantola(1984)在时间域实现。其将地震勘探领域的地下介质参数反演问题,看成模拟地震记录与实际采集地震记录残差以最小二乘目标函数为目标函数的最优化问题。由于该问题是高度非线性以及病态的。因此,在Born近似的理论框架下,该优化问题可以采用局部优化的方法近似求解。反演的梯度用伴随状态法求取,即将目标函数对模型参数求导通过以记录残差为震源的反传波场与正传波场的零延迟互相关得到。伴随状态法的提出避免了直接计算Jacobi矩阵,即一次梯度的求取只需要计算两次波场模拟,从而大大减少计算量,奠定了全波形反演的基础。但是受计算机性能的限制以及反演过程需要很多次迭代,计算效率的问题很长时间内制约了该方法的发展。到了20世纪90年代,Pratt等将全波形反演方法在频率域实现,即采用从低到高的几个离散的频率就可以实现全波形反演,同时从低频到高频的多尺度策略避免了反演陷入局部极小值。这大大推动了之后学者对全波形反演的研究,使该方法发展的越来越快。
全波形反演方法面临的主要问题就是跳周问题,即模拟记录和观测记录波形相位相差半个周期以上,导致波形匹配错误,从而导致反演结果出错。产生跳周的原因主要有两方面,一是反演的初始模型与真实地下模型相差太远,导致模拟记录波形与观测记录波形相差太大,从而导致跳周;二是缺少低频信息,地震记录中的低频信号携带了大尺度构造的信息且不易跳周。然而,实际观测记录中往往缺失低频信息,造成缺失低频信息的原因有很多种,如震源激发的地震子波缺失低频信息、检波器接收不到低频信息或者低频信息受噪声污染被人为切除等。对于传统的基于最小二乘目标函数的全波形反演对振幅信息的要求比较高,由于震源激发或者地震记录采集过程中造成的波形振幅错误对反演过程的影响极大。
发明内容
本发明的目的就在于针对上述现有技术的不足,提供一种基于零均值归一化互相关目标函数的时间域全波形反演方法。
本发明的目是通过以下技术方案实现的:
本发明基于零均值归一化互相关目标函数的时间域全波形反演方法,核心是提出一种零均值归一化互相关目标函数,该方法通过计算模拟记录和观测记录之间的零均值归一化互相关代替传统全波形反演最小二乘目标函数中波场残差的求取。通过归一化互相关使模拟数据与观测数据振幅的绝对值在-1到1之间,有效的减小了因振幅误差导致的反演结果错误。同时,由于振幅能量的归一化,提高了地震记录中远偏移距能量的权重,因为远偏移距数据携带了模型中大尺度构造的信息,因此该方法能在一定程度上减小跳周的发生。由于实际地震记录易受低频噪声的影响,且通常为非均值噪声,故引入了去均值项,通过对模拟记录和观测记录的去均值计算减小了非均值低频噪声的干扰。因此,该方法能在缺失低频信息或者低频信息被噪声污染的情况下准确的反演出地下介质的速度参数。
一种基于零均值归一化互相关目标函数的时间域全波形反演方法,包括以下步骤:
a、对实际地震观测记录进行预处理,对震源子波进行估计;
b、首先在预估速度范围建立线性递增初始模型,根据要求设定时间域全波形反演相关参数,包括地震子波主频f,低通滤波截断频率fluc,模型大小nz×nx,网格距dx,dz,采样总时间T,时间采样间隔dt,每个频段最大迭代次数itermax,最优化算法的迭代步长s,目标函数要求精度tol,模型速度估计的最大值vmax与最小值vmin;
c、用子波在初始模型上进行正演,得到模拟记录,对模拟记录和观测记录做低通滤波处理,得到低频段信号;
d、根据零均值归一化互相关原理构造目标函数:
其中J为目标函数,v为模型速度,d为实际采集的观测记录,u为对速度模型正演得到的模拟记录,上标“-”表示取平均值,为观测记录的平均值,为模拟记录的平均值,x为检波器位置,t为采样时间,令:
目标函数两端对模型速度v求导,可得得到梯度表达式:
Pf为时间域正传波场,L-1λ表示反传波场,其中L-1为反传算子,λ为伴随源,其表达式为:
e、利用L-BFGS优化算法对速度模型进行迭代更新,先反演出模型的大尺度构造,再以低频段反演结果做为初始模型进行全频段反演,将模型的细节构造反演出来,最终得高精度的地下速度模型。
步骤a,所述预处理过程包括:面波切除、多次波衰减、缺失地震道修复以及振幅补偿。
与现有技术相比,本发明的有益效果在于:本发明成功地将零均值归一化目标函数用到了全波形反演中,有效的解决了观测记录中低频信息缺失造成的跳周问题以及低频噪声污染导致的反演结果不准确等问题。
通过预先估计地下速度参数的取值范围建立线性递增的初始模型,采用混合震源进行正演模拟并采用多尺度策略选择适当的低通滤波截断频率,进行基于零均值归一化互相关目标函数的全波形反演,之后以低频带的反演结果做为初始模型进行全频带基于零均值归一化互相关目标函数的全波形反演。与常规基于最小二乘目标函数的时间域全波形反演相比,本发明解决了以下问题:
1、利用了归一化互相关的思想,提高了观测记录和模拟记录中远偏移距信号的权重,而远偏移距信号往往携带了地下介质大尺度构造的信息,相比于传统基于最小二乘目标函数的全波形反演,一定程度上缓解了跳周现象的发生,提高了反演的准确性。
2、归一化互相关的思想使观测记录和模拟记录的振幅绝对值在-1到1之间,从而减小了由于震源激发或者检波器接受故障导致的振幅能量错误的影响,可以在观测记录出现振幅错误的情况下准确的将地下介质参数反演出来,提高了反演对观测记录的容错率。
3、采用了零均值思想去除观测记录中非均值低频噪声的影响。由于实际地震勘探采集到的观测记录易受低频噪声的影响,这些噪声往往是非均值的有色噪声,使得低频段的全波形反演精度受到很大影响。零均值方法的实现对消除低频非均值噪声有很好的效果,大大提高了对受低频噪声污染数据的反演精度,提高了分辨率。
基于零均值归一化互相关目标函数的时间域全波形反演方法实现简单,计算过程与计算效率和传统基于最小二乘目标函数的全波形反演类似,在提高反演效果的同时不会增加额外的计算量,非常适合应用于实际地震勘探领域。
附图说明
图1基于零均值归一化互相关目标函数的时间域全波形反演方法流程图;
图2a 20Hz主频的雷克子波;
图2b雷克子波频谱;
图3a去除低频成分9Hz以下信息的雷克子波图;
图3b去除低频成分9Hz以下信息的雷克子波频谱;
图4真实模型图;
图5线性递增初始模型图;
图6a基于最小二乘目标函数的20Hz低通滤波全波形反演结果图;
图6b基于零均值归一化互相关目标函数的20Hz低通滤波全波形反演结果图;
图7a以图6a为初始模型进行基于最小二乘目标函数的全频带全波形反演结果图;
图7b以图6b为初始模型进行基于零均值归一化互相关目标函数的全频带全波形反演结果图;
图8a观测记录振幅乘0.03后基于最小二乘目标函数的20Hz低通滤波全波形反演结果图;
图8b观测记录振幅乘0.03后基于零均值归一化互相关目标函数的20Hz低通滤波全波形反演结果图;
图9含强低频噪声的观测记录;
图10a基于最小二乘目标函数的20Hz低通滤波全波形反演结果图;
图10b基于互相关目标函数的20Hz低通滤波全波形反演结果图;
图10c基于零均值归一化互相关目标函数的20Hz低通滤波全波形反演结果图;
图11a以图10b为初始模型基于互相关目标函数的全频带全波形反演结果图;
图11b以图10c为初始模型基于零均值归一化互相关目标函数的全频带全波形反演结果图。
具体实施方式
下面结合附图和实例对本发明进一步的详细说明:
如图1所示,本发明基于零均值归一化互相关目标函数的时间域全波形反演方法,包括以下步骤:
a、实际地震观测记录首需进行预处理,对震源子波进行估计,其中预处理过程包括:面波切除、多次波衰减、缺失地震道修复、振幅补偿等。
b、根据预估计的地下介质速度范围建立线性递增的初始模型。根据要求设定时间域全波形反演相关参数,包括地震子波主频f,低通滤波截断频率fluc,模型大小nz×nx,网格距dx,dz,采样总时间T,时间采样间隔dt,每个频段最大迭代次数itermax,最优化算法的迭代步长s,目标函数要求精度tol,模型速度估计的最大值vmax与最小值vmin。
c、将图2a所示子波作为震源混合编码后在初始模型上进行正演计算得到模拟记录,对模拟记录和观测记录做低通滤波处理,得到低频段信号;
d、根据零均值归一化互相关原理构造目标函数:
其中J为目标函数,v为模型速度,d为实际采集的观测记录,u为对速度模型正演得到的模拟记录,上标“-”表示取平均值,为观测记录的平均值,为模拟记录的平均值,x为检波器位置,t为采样时间,令:
目标函数两端对模型速度v求导,可得梯度表达式:
简化后可表示为:
其中λ为伴随源,其表达式为:
通过链式法则可以得到基于伴随状态法的时间域全波形反演梯度公式:
其中Pf为时间域正传波场,L-1λ表示反传波场,其中L-1为反传算子。
e、反演基于L-BFGS优化算法进行迭代更新,并通过Wolfe收敛准则寻找合适的步长;其迭代公式表示为:
mk+1=mk-αkHkgk
其中,Hk为近似Hessian矩阵的逆矩阵,mk为模型更新参数,gk为模型更新梯度,αk为Wolfe线性搜索得到的步长,k表示迭代次数。
L-BFGS优化算法在迭代计算过程中需要保存的矩阵个数很少,对计算机内存要求小,同时,对用于更新Hessian矩阵,其迭代公式如下:
Hk+1=Vk THkVk+ρksksk T
sk=mk+1-mk,yk=gk+1-gk
其中,Vk为计算过程的中间变量,sk为计算过程的中间变量,ρk为中间变量,I表示单位矩阵,Hk+1是根据向量对{sk,yk}和Hk计算得到,只储存n个向量对来隐式表达Hessian矩阵的逆矩阵。Hkgk的乘积可以通过梯度gk与向量对{sk,yk}之间一系列向量的内积与向量的和来获得的。在每一次更新完成后,上一步向量对将被当前新向量对{sk+1,yk+1}取代。因此,向量对集合中包含最近n次迭代的曲率信号。在时实际中,当n≥5时都能获得较满意的结果。L-BFGS优化算法的初始近似Hessian矩阵每一次迭代中都不同。近似Hessian矩阵的逆矩阵Hk需满足以下更新公式:
模型的更新方向,通过以下方法实现:
(1)令则q=q-αiyi,其中,i=k-1,k-2,…k-n;αi为计算过程中的中间变量;
(2)令则r=r-si(αi-β),其中,i=k-n,k-n+1,…,k-1;为初始近似Hessian矩阵;
(3)通过上述步骤得到更新量Hkgk=r。
通过上述方法求得模型的更新量Hkgk,然后再通过Wolfe线性搜索获得合适的步长αk进行更新迭代。L-BFGS优化算法有效的克服了显式求取Hessian矩阵及其逆的困难,其具有超线性收敛速度,计算效率高,占用内存小,精度高等优点,较适合求解大规模非线性优化问题。
最后,判断反演结果是否满足设置的终止条件,即反演结果与真实模型相差很小,若果是则停止迭代,如果不是,则继续计算直到满足终止条件为止,结束计算。
本发明程序是在MATLAB2016b软件框架下编写完成,根据相应的计算要求,搭建MATLAB工作库的安装环境,并安装MATLAB并行计算工具箱(Parallel ComputingToolbox)。
实施例1
根据勘探要求,将Parallel Computing Toolbox和MATLAB DistributedComputing Server(R2016b)在Windows 10专业版系统下进行安装,进行MATLAB并行平台的搭建。利用Marmousi进行缺失低频信息的全波形反演测试,真实模型(图4)和初始模型(图5)。
模型参数如下:
表1缺失低频信息全波形反演测试参数
网格大小 | 网格距 | 测线长度 | 纵向深度 | fluc | 频带宽度 |
69×192 | 12.5m | 2400m | 862.5m | 20Hz | 9~20Hz |
模型网格大小为69×192,网格距12.5,横向距离为2400m,纵向深度为862.5m,模型中地震波速度范围从1500m/s到4000m/s,地震检波器安置在模型表面,每一个网格点都是一个检波器,且检波器之间的间隔为12.5m,一个混合震源中含有12个震源。震源选用缺失9Hz以下信息20Hz主频的雷克子波(图3a),采样间隔为0.001s,实际采样总长度为2s,频率范围从9Hz到20Hz。
图6a和图6b为缺失低频信息的低频带反演结果对比图,传统基于最小二乘目标函数的全波形反演结果(图6a)在模型左上方有不同于真实模型的区域,产生该区域的原因是反演过程中产生了跳周,使该区域的速度更新量错误。而基于零均值归一化互相关目标函数的全波形反演结果(图6b)大致反演出了真实模型的大尺度构造,且没有异常区域。以图6为初始模型,继续进行全波形反演,得到全频带反演结果(图7a和图7b)。基于最小二乘目标函数的全频带反演结果(图7a)中出现了明显的异常区域,这是由于跳周现象得不到解决,导致该区域错误的速度更新量逐渐累加。基于本发明的全波形反演结果(图7b)比较精确的将真实模型反演出来,没有出现速度异常区域,从而使反演结果有很高的可信度。
实施例2
根据勘探要求,将Parallel Computing Toolbox和MATLAB DistributedComputing Server(R2016b)在Windows 10专业版系统下进行安装,进行MATLAB并行平台的搭建。利用Marmousi进行振幅错误时的全波形反演测试,真实模型(图4)和初始模型(图5)。
模型参数如下:
表2观测记录振幅错误低频段全波形反演测试参数
网格大小 | 网格距 | 测线长度 | 纵向深度 | fluc | 频带宽度 |
69×192 | 12.5m | 2400m | 862.5m | 20Hz | 0~20Hz |
模型网格大小为69×192,网格距12.5,横向距离为2400m,纵向深度为862.5m,模型中地震波速度范围从1500m/s到4000m/s,地震检波器安置在模型表面,每一个网格点都是一个检波器,且检波器之间的间隔为12.5m,一个混合震源中含有12个震源。震源选用20Hz主频的雷克子波(图2a),采样间隔为0.001s,实际采样总长度为2s,频率范围从0Hz到20Hz。
反演结果如图8所示,传统基于最小二乘目标函数的全波形反演结果(图8a)完全偏离真实模型。由于将观测记录振幅整体乘以系数0.03,导致观测记录与模拟记录振幅差异异常增大,使模型各个位置的速度更新量都出现了错误值,经过不断迭代计算,错误的更新量不断累积最终导致反演完全失败。而基于零均值归一化互相关目标函数将观测数据与模拟数据的振幅绝对值调整到-1到1之间,使错误振幅产生的影响降到很低,因而其全波形反演结果(图8b)反演出了真实模型的大尺度构造,且没有异常区域。
实施例3
根据勘探要求,将Parallel Computing Toolbox和MATLAB DistributedComputing Server(R2016b)在Windows 10专业版系统下进行安装,进行MATLAB并行平台的搭建。利用Marmousi进行观测记录含强低频噪声时的全波形反演测试,加噪观测记录(图9),真实模型(图4)和初始模型(图5)。。
模型参数如下:
表3含强低频噪声低频段全波形反演测试参数
网格大小 | 网格距 | 测线长度 | 纵向深度 | fluc | 频带宽度 |
69×192 | 12.5m | 2400m | 862.5m | 20Hz | 0~20Hz |
模型网格大小为69×192,网格距12.5,横向距离为2400m,纵向深度为862.5m,模型中地震波速度范围从1500m/s到4000m/s,地震检波器安置在模型表面,每一个网格点都是一个检波器,且检波器之间的间隔为12.5m,一个混合震源中含有12个震源。震源选用20Hz主频的雷克子波(图2a),采样间隔为0.001s,实际采样总长度为2s,频率范围从0Hz到20Hz。
图10a-图10c为含强含强低频噪声观测记录的低频带反演结果对比图。基于最小二乘目标函数的低频段全波形反演结果(图10a)完全偏离真实模型,这是由于低频噪声能量过强湮没了观测记录中的有效信号,导致观测记录与模拟记录振幅、相位等信息相差过大,造成速度更量错误累计,最终偏离真实模型。基于互相关目标函数的低频段全波形反演结果(图10b)相比图10a有了明显的改善,互相关思想的使用使伴随源中有效信号能量权重提高,噪声被一定程度上压制了,大致反演出模型的大尺度构造。但是该反演结果依然受噪声影响比较严重,导致部分区域速度更新量错误,不能够准确反映真实模型速度。基于零均值互相关目标函数的低频段全波形反演结果(图10c)基本准确的将地下速度模型的大尺度构造反演出来,并且没有异常区域。由于伴随源中,观测记录和模拟记录都减去了各自的平均值,使记录中的低频非均值噪声被有效压制,从而将噪声的影响减到最小,得到的反演结果趋于真实模型。
以图10a和图10b为初始模型继续进行全频带全波形反演结果如图11a和图11b所示。由于低频段受噪声干扰较强,基于互相关目标函数的全频带反演结果(图11a)与真实速度模型偏离较大。而基于零均值归一化互相关目标函数的全波形反演结果(图11b)逼近真实模型,没有出现速度异常值,反演精度高。
Claims (2)
1.一种基于零均值归一化互相关目标函数的时间域全波形反演方法,其特征在于,包括以下步骤:
a、对实际地震观测记录进行预处理,对震源子波进行估计;
b、首先在预估速度范围建立线性递增初始模型,根据要求设定时间域全波形反演相关参数,包括地震子波主频f,低通滤波截断频率fluc,模型大小nz×nx,网格距dx,dz,采样总时间T,时间采样间隔dt,每个频段最大迭代次数itermax,最优化算法的迭代步长s,目标函数要求精度tol,模型速度估计的最大值vmax与最小值vmin;
c、用子波在初始模型上进行正演,得到模拟记录,对模拟记录和观测记录做低通滤波处理,得到低频段信号;
d、根据零均值归一化互相关原理构造目标函数:
其中J为目标函数,v为模型速度,d为实际采集的观测记录,u为对速度模型正演得到的模拟记录,上标“-”表示取平均值,为观测记录的平均值,为模拟记录的平均值,x为检波器位置,t为采样时间,令:
目标函数两端对模型速度v求导,可得得到梯度表达式:
Pf为时间域正传波场,L-1λ表示反传波场,其中L-1为反传算子,λ为伴随源,其表达式为:
e、利用L-BFGS优化算法对速度模型进行迭代更新,先反演出模型的大尺度构造,再以低频段反演结果做为初始模型进行全频段反演,将模型的细节构造反演出来,最终得高精度的地下速度模型。
2.根据权利要求1所述的一种基于零均值归一化互相关目标函数的时间域全波形反演方法,其特征在于,步骤a,所述预处理过程包括:面波切除、多次波衰减、缺失地震道修复以及振幅补偿。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811550994.9A CN109407152B (zh) | 2018-12-18 | 2018-12-18 | 基于零均值归一化互相关目标函数的时间域全波形反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811550994.9A CN109407152B (zh) | 2018-12-18 | 2018-12-18 | 基于零均值归一化互相关目标函数的时间域全波形反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109407152A true CN109407152A (zh) | 2019-03-01 |
CN109407152B CN109407152B (zh) | 2019-11-22 |
Family
ID=65459855
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811550994.9A Expired - Fee Related CN109407152B (zh) | 2018-12-18 | 2018-12-18 | 基于零均值归一化互相关目标函数的时间域全波形反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109407152B (zh) |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111239806A (zh) * | 2020-03-09 | 2020-06-05 | 吉林大学 | 基于振幅增量编码的时间域全波形反演方法 |
CN112198550A (zh) * | 2019-07-08 | 2021-01-08 | 中国石油天然气集团有限公司 | 一种时移地震数据可重复性度量方法及装置 |
CN113064203A (zh) * | 2021-03-26 | 2021-07-02 | 中国海洋大学 | 共轭梯度归一化lsrtm方法、系统、存储介质及应用 |
CN113156493A (zh) * | 2021-05-06 | 2021-07-23 | 中国矿业大学 | 一种使用归一化震源的时频域全波形反演方法及装置 |
CN116070706A (zh) * | 2023-01-16 | 2023-05-05 | 华东计算技术研究所(中国电子科技集团公司第三十二研究所) | 用于超导量子芯片标定的自动化处理系统及方法 |
CN117849879A (zh) * | 2023-12-20 | 2024-04-09 | 京全品质能源技术(北京)有限公司 | 波动方程反演方法及装置、电子设备以及存储介质 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140142860A1 (en) * | 2012-11-20 | 2014-05-22 | International Business Machines Corporation | Efficient wavefield compression in seismic imaging |
CN107422379A (zh) * | 2017-07-27 | 2017-12-01 | 中国海洋石油总公司 | 基于局部自适应凸化方法的多尺度地震全波形反演方法 |
CN107765308A (zh) * | 2017-10-12 | 2018-03-06 | 吉林大学 | 基于褶积思想与精确震源的重构低频数据频域全波形反演方法 |
CN108345031A (zh) * | 2018-01-11 | 2018-07-31 | 吉林大学 | 一种弹性介质主动源和被动源混采地震数据全波形反演方法 |
CN108549100A (zh) * | 2018-01-11 | 2018-09-18 | 吉林大学 | 基于非线性高次拓频的时间域多尺度全波形反演方法 |
-
2018
- 2018-12-18 CN CN201811550994.9A patent/CN109407152B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140142860A1 (en) * | 2012-11-20 | 2014-05-22 | International Business Machines Corporation | Efficient wavefield compression in seismic imaging |
CN107422379A (zh) * | 2017-07-27 | 2017-12-01 | 中国海洋石油总公司 | 基于局部自适应凸化方法的多尺度地震全波形反演方法 |
CN107765308A (zh) * | 2017-10-12 | 2018-03-06 | 吉林大学 | 基于褶积思想与精确震源的重构低频数据频域全波形反演方法 |
CN108345031A (zh) * | 2018-01-11 | 2018-07-31 | 吉林大学 | 一种弹性介质主动源和被动源混采地震数据全波形反演方法 |
CN108549100A (zh) * | 2018-01-11 | 2018-09-18 | 吉林大学 | 基于非线性高次拓频的时间域多尺度全波形反演方法 |
Non-Patent Citations (2)
Title |
---|
杨贺龙: "基于波场优化匹配的高精度全波形反演方法研究", 《中国博士学位论文全文数据库 基础科学辑》 * |
胡勇 等: "基于自适应非稳态相位校正的时频域多尺度全波形反演", 《地球物理学报》 * |
Cited By (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112198550A (zh) * | 2019-07-08 | 2021-01-08 | 中国石油天然气集团有限公司 | 一种时移地震数据可重复性度量方法及装置 |
CN112198550B (zh) * | 2019-07-08 | 2023-09-26 | 中国石油天然气集团有限公司 | 一种时移地震数据可重复性度量方法及装置 |
CN111239806A (zh) * | 2020-03-09 | 2020-06-05 | 吉林大学 | 基于振幅增量编码的时间域全波形反演方法 |
CN111239806B (zh) * | 2020-03-09 | 2021-06-22 | 吉林大学 | 基于振幅增量编码的时间域全波形反演方法 |
CN113064203A (zh) * | 2021-03-26 | 2021-07-02 | 中国海洋大学 | 共轭梯度归一化lsrtm方法、系统、存储介质及应用 |
CN113156493A (zh) * | 2021-05-06 | 2021-07-23 | 中国矿业大学 | 一种使用归一化震源的时频域全波形反演方法及装置 |
CN113156493B (zh) * | 2021-05-06 | 2022-02-18 | 中国矿业大学 | 一种使用归一化震源的时频域全波形反演方法及装置 |
CN116070706A (zh) * | 2023-01-16 | 2023-05-05 | 华东计算技术研究所(中国电子科技集团公司第三十二研究所) | 用于超导量子芯片标定的自动化处理系统及方法 |
CN117849879A (zh) * | 2023-12-20 | 2024-04-09 | 京全品质能源技术(北京)有限公司 | 波动方程反演方法及装置、电子设备以及存储介质 |
Also Published As
Publication number | Publication date |
---|---|
CN109407152B (zh) | 2019-11-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109407152B (zh) | 基于零均值归一化互相关目标函数的时间域全波形反演方法 | |
CN109407151B (zh) | 基于波场局部相关时移的时间域全波形反演方法 | |
CN109459789B (zh) | 基于振幅衰减与线性插值的时间域全波形反演方法 | |
Zhang et al. | Parameter prediction of hydraulic fracture for tight reservoir based on micro-seismic and history matching | |
CN104502997B (zh) | 一种利用裂缝密度曲线预测裂缝密度体的方法 | |
CN101630018B (zh) | 一种控制全声波方程反演过程的地震勘探数据处理方法 | |
CN106772577B (zh) | 基于微地震数据和spsa优化算法的震源反演方法 | |
CN102854528B (zh) | 粒子群优化算法叠前非线性反演方法 | |
CN103913774B (zh) | 基于微地震事件的储层地质力学参数反演方法 | |
RU2620785C1 (ru) | Способ определения местоположения очага микросейсмического события | |
MXPA04012989A (es) | Metodo para monitorear eventos sismicos. | |
CN110058302A (zh) | 一种基于预条件共轭梯度加速算法的全波形反演方法 | |
CN105277978A (zh) | 一种确定近地表速度模型的方法及装置 | |
CN109188519B (zh) | 一种极坐标下的弹性波纵横波速度反演系统及方法 | |
CN103513277B (zh) | 一种地震地层裂隙裂缝密度反演方法及系统 | |
CN105388518A (zh) | 一种质心频率与频谱比联合的井中地震品质因子反演方法 | |
CN101021568A (zh) | 基于最大能量旅行时计算的三维积分叠前深度偏移方法 | |
CN104570101A (zh) | 一种基于粒子群算法的avo三参数反演方法 | |
CN113740901B (zh) | 基于复杂起伏地表的陆上地震数据全波形反演方法及装置 | |
CN110579795B (zh) | 基于被动源地震波形及其逆时成像的联合速度反演方法 | |
CN101545986A (zh) | 基于最大能量旅行时计算的三维积分叠前深度偏移方法 | |
CN110007340A (zh) | 基于角度域直接包络反演的盐丘速度密度估计方法 | |
CN109633749A (zh) | 基于散射积分法的非线性菲涅尔体地震走时层析成像方法 | |
WO2008056267A2 (en) | System and method for determing seismic event location | |
CN104199088B (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 | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20191122 Termination date: 20211218 |
|
CF01 | Termination of patent right due to non-payment of annual fee |