CN113552625A - 一种用于常规陆域地震数据的多尺度全波形反演方法 - Google Patents
一种用于常规陆域地震数据的多尺度全波形反演方法 Download PDFInfo
- Publication number
- CN113552625A CN113552625A CN202110683680.1A CN202110683680A CN113552625A CN 113552625 A CN113552625 A CN 113552625A CN 202110683680 A CN202110683680 A CN 202110683680A CN 113552625 A CN113552625 A CN 113552625A
- Authority
- CN
- China
- Prior art keywords
- waveform inversion
- seismic data
- full
- model
- frequency
- 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 54
- 238000003325 tomography Methods 0.000 claims abstract description 40
- 238000011156 evaluation Methods 0.000 claims abstract description 22
- 238000007781 pre-processing Methods 0.000 claims abstract description 13
- 238000004088 simulation Methods 0.000 claims abstract description 8
- 238000009499 grossing Methods 0.000 claims description 13
- 238000012937 correction Methods 0.000 claims description 10
- 238000010521 absorption reaction Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 5
- 230000003068 static effect Effects 0.000 claims description 4
- 230000005284 excitation Effects 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 7
- 238000013508 migration Methods 0.000 abstract description 7
- 230000005012 migration Effects 0.000 abstract description 7
- 238000004458 analytical method Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000013524 data verification Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000007792 addition Methods 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 239000002360 explosive Substances 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000009776 industrial production Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000006467 substitution reaction Methods 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/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/62—Physical property of subsurface
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
本发明提供了一种用于常规陆域地震数据的多尺度全波形反演方法,具体涉及地震勘探技术领域。该方法利用原始地震数据的层析成像结果进行周跳风险评价,建立全波形反演模型,对原始地震数据进行预处理后与正演模拟地震数据振幅匹配,再由时间域变换至频率域作为多尺度全波形反演数据,设置反演参数,基于层剥离的Laplace‑Fourier域多尺度全波形反演方法进行多尺度全波形反演,通过由近及远遍历所有偏移距、由小到大遍历所有Laplace常数以及由低频到高频的遍历所有频率,获得多尺度全波形反演结果。本发明克服了全波形反演方法存在的周跳难题,充分利用了实测数据中的低频信息,有效提高了陆域地震数据的全波形反演效果,为常规陆域的地质解释提供了依据。
Description
技术领域
本发明涉及地震勘探技术领域,具体涉及一种用于常规陆域地震数据的多尺度全波形反演方法。
背景技术
目前,全波形反演已成功应用于海域地震数据中,但是在常规陆域地震数据的应用上还存在诸多难点。与海域地震数据相比,常规陆域地震数据偏移距小且缺少低频信息,由于常规工业检波器的频率为10Hz,而陆域地震数据的频率往往高于5Hz,所以常规陆域地震数据对于全波形反演方法而言频率太高,采用全波形反演方法处理常规陆域地震数据,无法得到较好的效果。
同时,全波形反演方法的应用受限于梯度类反演方法普遍存在的“周跳”现象。“周跳”现象是指反演时,当初始记录与实际记录之间差值小于1/2周期时,反演收敛得到全局最优值,而当初始记录与实际记录之间差值大于1/2周期时,反演收敛得到局部最优值,无法得到正确的结果。因此,亟需提出一种用于常规陆域地震数据的多尺度全波形反演方法,提高全波形反演结果的准确性。
发明内容
针对现有技术中全波形反演方法受限于梯度类反演方法普遍存在的“周跳”难题,本发明提出了一种用于常规陆域地震数据的多尺度全波形反演方法,以传统频率域多尺度为基础,通过增加层剥离多尺度和偏移距多尺度,实现了“由近到远,由浅到深,由低频到高频”的多尺度反演,提高了全波形反演的精度,有利于常规陆域的地质解释。
为了实现上述目的,本发明采用如下技术方案:
一种用于常规陆域地震数据的多尺度全波形反演方法,具体包括如下步骤:
步骤1,输入原始地震数据,获取原始地震数据的初至旅行时,基于二维射线追踪地震层析成像方法,确定层析成像结果,对层析成像结果进行周跳风险评价,确定全波形反演的速度模型,根据全波形反演模型的速度模型,确定全波形反演的密度模型,建立全波形反演模型;
步骤2,对原始地震数据进行预处理,获得校正后的地震数据,利用全波形反演的速度模型进行正演模拟得到正演模拟地震数据,基于均方根振幅最小二乘拟合方法,对正演模拟地震数据与校正后的地震数据进行振幅匹配,并将校正后的地震数据由时间域变换至频率域,作为多尺度全波形反演数据;
步骤3,输入多尺度全波形反演数据和全波形反演模型,设置频率组、拉普拉斯常数组、偏移距组和迭代次数,基于层剥离的Laplace-Fourier域多尺度全波形反演方法,通过设置由近及远遍历偏移距组内偏移距的内部循环、由小到大遍历拉普拉斯常数组内拉普拉斯常数的中间循环以及由低到高的遍历频率组内频率的外部循环,利用全波形反演模型进行全波形反演,获得多尺度全波形反演结果,根据多尺度全波形反演结果进行地质解释。
优选地,所述步骤1中,具体包括以下步骤:
步骤1.1,输入原始地震数据,拾取原始地震数据的初至旅行时获得实际旅行时,结合原始地震数据选取初至数据,设置层析成像参数;
步骤1.2,根据初至数据,基于二维射线追踪地震层析成像方法,得到层析成像结果,确定层析成像速度模型;
步骤1.3,对层析成像速度模型进行正演模拟,计算各炮点激发检波器位置处的初至旅行时,获得各炮点的正演旅行时,结合各炮点的实际旅行时,计算相同炮点和检波器位置处正演旅行时与实际旅行时之间的差值;
步骤1.4,设置周跳风险评价频率,计算评价周期,进行周跳风险评价,判断层析成像速度模型是否会产生周跳现象,若差值小于1/2个评价周期,则层析成像速度模型不会产生周跳现象,执行步骤1.6;若差值不小于1/2个评价周期,则层析成像速度模型会产生周跳现象,执行步骤1.5;
步骤1.5,选取层析成像速度模型中差值不小于1/2个评价周期的区域,增加炮点数据,继续进行正演模拟,返回步骤1.3;
步骤1.6,输出层析成像速度模型作为全波形反演模型中的速度模型,基于Gardner经验关系,确定全波形反演模型中的密度模型,设置品质因子模型,建立全波形反演模型。
优选地,所述步骤1.3中,基于频率域2D粘弹声波方程进行正演模拟,如式(1)所示:
式中,ω为角频率;κ(x,z)为体积模量;p(x,z,ω)为频率域地震波场;ρ(x,z)为密度;s(x,z,ω)为震源项。
优选地,所述步骤1.6中,基于Gardner经验关系,利用全波形反演模型的速度模型确定全波形反演模型中的密度模型,计算公式为:
式中,ρ为密度,单位为g/cm3;Vp为纵波速度,单位为m/s。
优选地,所述步骤2中,具体包括以下步骤:
步骤2.1,原始地震数据预处理,预处理包括静校正、去噪、去面波、振幅补偿、一次性校正和初至切除,得到校正后的地震数据;
步骤2.2,利用全波形反演的速度模型对校正后的地震数据进行正演模拟,得到正演模拟地震数据,基于均方根振幅最小二乘拟合方法,将正演模拟地震数据的振幅与校正后地震数据的振幅逐炮匹配,使得正演模拟地震数据与校正后地震数据的尺度相同;
正演模拟地震数据与校正后地震数据的振幅匹配公式为:
y=Cx (3)
式中,x为正演模拟地震数据的均方根振幅;y为校正后地震数据的均方根振幅;C为缩放系数;
步骤2.3,将校正后的地震数据由时间域变换至频率域,得到多尺度全波形反演数据;
校正后地震数据的变换公式如下:
式中,ω为角频率;t为时间;f(t)为校正后地震数据的时间域信号;F(ω)为校正后地震数据的频率域信号。
优选地,所述步骤3中,具体包括以下步骤:
步骤3.1,输入多尺度全波形反演数据和全波形反演模型,分别设置频率组、拉普拉斯常数组、偏移距组和迭代次数,其中,频率组内按照频率由低到高的顺序设置有a个频率,拉普拉斯常数组内按照数值由小到大的顺序设置有b个拉普拉斯常数,偏移距组内按照数值由小到大的顺序设置有c个偏移距,分别对频率组、拉普拉斯常数组和偏移距组内各参数进行编号;
步骤3.2,分别设置频率、拉普拉斯常数和偏移距的初始值,其中,将频率组内的最低频率设置为频率的初始值,将拉普拉斯常数设置为拉普拉斯常数组内的最小值,将偏移距设置为偏移距组内的最小值;
步骤3.3,设置梯度吸收因子ε和高斯平滑因子ζ;
步骤3.4,根据当前设置的偏移距和频率从多尺度全波形反演数据中提取地震数据,利用全波形反演模型进行全波形反演,得到全波形反演结果,并将本次迭代计算得到的全波形反演结果设置为全波形反演模型;
步骤3.5,判断是否完成设置的迭代次数;
若未完成设置的迭代次数,则返回步骤3.4继续进行全波形反演;若已完成设置的迭代次数,则执行步骤3.6;
步骤3.6,判断是否遍历偏移距组内全部的偏移距;
若未遍历完,则根据偏移距组对当前偏移距Oi进行更新,i=1,…,c-1,i为当前偏移距在偏移距组内的序号,更新后的偏移距为Oi+1,返回步骤3.4继续进行全波形反演;若已遍历完,则对偏移距进行初始化,设置偏移距为O1,重新遍历偏移距组内的偏移距,执行步骤3.7;
步骤3.7,判断是否遍历完拉普拉斯常数组内全部的拉普拉斯常数;
若未遍历完,则根据拉普拉斯常数组对当前拉普拉斯常数Lj进行更新,j=1,…,b-1,j为当前拉普拉斯常数在拉普拉斯常数组内的序号,更新后的拉普拉斯常数为Lj+1,返回步骤3.4继续进行全波形反演;若已遍历完,则对拉普拉斯常数进行初始化,设置拉普拉斯常数为L1,重新遍历拉普拉斯常数组内的拉普拉斯常数,执行步骤3.8;
步骤3.8,判断是否遍历完频率组内全部的频率;
若未遍历完,则根据频率组对当前频率fn进行更新,n=1,…,a-1,n为当前频率在频率组内的序号,更新后的频率为fn+1,返回步骤3.4继续进行全波形反演;若已遍历完,则将当前全波形反演模型作为反演速度模型,执行步骤3.9;
步骤3.9,输出多尺度全波形反演结果,并根据多尺度全波形反演结果进行地质解释。
优选地,所述步骤3.4中,反演过程中全波形反演模型的参数迭代更新公式为:
式中,dobs为实测地震数据,dcal为正演模拟地震数据;C(m)为加权最小二乘目标函数,为转置操作,Wd为数据预处理算子;mk+1为第k+1次迭代得到的属性模型,mk为第k次迭代得到的属性模型,Δmk为第k次模型修正量;为目标函数的最速下降方向,αk为沿最速下降方向的步长,ε为梯度吸收因子,I为单位矩阵,ζ为高斯平滑因子,为高斯-牛顿海森矩阵的近似对角矩阵。
本发明所带来的有益技术效果:
本发明为了克服现全波形反演方法的周跳难题,利用初始模型控制初始数据与实测数据之间的差值,充分利用实测数据中的低频信息,基于层剥离方法在Laplace-Fourier域进行全波形反演,通过依次优化小偏移距浅层地震数据、大偏移距浅层地震数据、小偏移距深层地震数据和大偏移距深层地震数据,采用“由近到远,由浅到深,由低频到高频”的多尺度反演策略,有效提高了速度反演的精度,为常规陆域的地质解释提供了依据。
附图说明
图1为一种用于常规陆域地震数据的多尺度全波形反演方法流程图。
图2为全波形反演模型的构建流程图。
图3为基于层剥离的Laplace-Fourier多尺度全波形反演方法流程图。
图4为反演速度与地温曲线对比图;图(a)为地温曲线,图(b)为反演速度曲线。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明:
实施例1
本发明提出了一种用于常规陆域地震数据的多尺度全波形反演方法,如图1所示,具体包括如下步骤:
步骤1,建立全波形反演模型,如图2所示,具体包括以下步骤:
步骤1.1,输入原始地震数据,拾取原始地震数据的初至旅行时获得实际旅行时,结合原始地震数据选取初至数据,设置层析成像参数。
步骤1.2,根据初至数据,基于二维射线追踪地震层析成像方法,得到层析成像结果,确定层析成像速度模型。
步骤1.3,对层析成像速度模型进行正演模拟,计算各炮点激发检波器位置处的初至旅行时,获得各炮点的正演旅行时,结合各炮点的实际旅行时,计算相同炮点和检波器位置处正演旅行时与实际旅行时之间的差值。
步骤1.4,设置周跳风险评价频率,计算评价周期,进行周跳风险评价,判断层析成像速度模型是否会产生周跳现象,若差值小于1/2个评价周期,则层析成像速度模型不会产生周跳现象,执行步骤1.6;若差值不小于1/2个评价周期,则层析成像速度模型会产生周跳现象,执行步骤1.5。
步骤1.5,选取层析成像速度模型中差值不小于1/2个评价周期的区域,增加炮点数据,继续进行正演模拟,返回步骤1.3。
步骤1.6,输出层析成像速度模型作为全波形反演模型中的速度模型,基于Gardner经验关系,确定全波形反演模型中的密度模型,设置品质因子模型,建立全波形反演模型。
步骤2,通过对原始地震数据进行预处理,获取多尺度全波形反演数据,具体包括以下步骤:
步骤2.1,原始地震数据预处理,预处理包括静校正、去噪、去面波、振幅补偿、一次性校正和初至切除,得到校正后的地震数据。
步骤2.2,利用全波形反演的速度模型对校正后的地震数据进行正演模拟,得到正演模拟地震数据,基于均方根振幅最小二乘拟合方法,将正演模拟地震数据的振幅与校正后地震数据的振幅逐炮匹配,使得正演模拟地震数据与校正后地震数据的尺度相同。
步骤2.3,将校正后的地震数据由时间域变换至频率域,获得多尺度全波形反演数据。
步骤3,多尺度全波形反演,如图3所示,具体包括以下步骤:
步骤3.1,输入多尺度全波形反演数据和全波形反演模型,分别设置频率组f=[f1,…,fa]、拉普拉斯常数组L=[L1,…,Lb]、偏移距组O=[O1,…,Oc]和迭代次数,其中,频率组内的频率按照由低频到高频的顺序排列,拉普拉斯常数组内的拉普拉斯常数按照由小到大的顺序排列,偏移距组内的偏移距按照由小到大的顺序排列,分别对频率组、拉普拉斯常数组和偏移距组内的参数进行编号。
步骤3.2,分别设置反演频率、拉普拉斯常数和偏移距的初始值,其中,将频率组内的最低频率设置为频率的初始值,将拉普拉斯常数设置为拉普拉斯常数组内的最小值,将偏移距设置为偏移距组内的最小值,即设置频率的初始值为f1、拉普拉斯常数的初始值为L1、偏移距的初始值为O1。
步骤3.3,设置梯度吸收因子ε和高斯平滑因子ζ。
步骤3.4,根据当前设置的偏移距和频率从多尺度全波形反演数据中提取地震数据,利用全波形反演模型进行全波形反演,得到全波形反演结果,并将本次迭代计算得到的全波形反演结果设置为全波形反演模型。
步骤3.5,判断是否完成设置的迭代次数。
若未完成设置的迭代次数,则返回步骤3.4继续进行全波形反演;若已完成设置的迭代次数,则执行步骤3.6。
步骤3.6,判断是否遍历偏移距组内全部的偏移距。
若未遍历完,则根据偏移距组对当前偏移距Oi进行更新,i=1,…,c-1,i为当前偏移距在偏移距组内的序号,更新后的偏移距为Oi+1,返回步骤3.4继续进行全波形反演;若已遍历完,则对偏移距进行初始化,设置偏移距为O1,重新遍历偏移距组内的偏移距,执行步骤3.7。
步骤3.7,判断是否遍历完拉普拉斯常数组内全部的拉普拉斯常数。
若未遍历完,则根据拉普拉斯常数组对当前拉普拉斯常数Lj进行更新,j=1,…,b-1,j为当前拉普拉斯常数在拉普拉斯常数组内的序号,更新后的拉普拉斯常数为Lj+1,返回步骤3.4继续进行全波形反演;若已遍历完,则对拉普拉斯常数进行初始化,设置拉普拉斯常数为L1,重新遍历拉普拉斯常数组内的拉普拉斯常数,执行步骤3.8。
步骤3.8,判断是否遍历完频率组内全部的频率。
若未遍历完,则根据频率组对当前频率fn进行更新,n=1,…,a-1,k为当前频率在频率组内的序号,更新后的频率为fn+1,返回步骤3.4继续进行全波形反演;若已遍历完,则将当前全波形反演模型作为反演速度模型,执行步骤3.9。
步骤3.9,输出多尺度全波形反演结果,并根据多尺度全波形反演结果进行地质解释。
实施例2
将本发明提出的一种用于常规陆域地震数据的多尺度全波形反演方法应用于理想模型中进行测试,取得了理想的反演效果。
根据某地区的地震剖面和岩石物性资料,构建二维纵波速度模型,基于Gardner经验关系,将纵波速度转化为密度,构建密度模型,设置品质因子模型为常数,进行多尺度全波形反演,多尺度全波形反演过程中所采用的参数如表1所示。
多尺度全波形反演过程中,根据该地区实际地震数据的观测系统,设置正演地震数据的观测系统最大偏移距为14.4km,正演地震数据观测系统的最大偏移距与实际地震数据观测系统相同。设置偏移距组内偏移距分别为5km和14.4km,设置Laplace常数组内Laplace常数分别为0.5和1.0;梯度吸收因子分为两组,分别为1e-2和1e-6;高斯平滑因子分为两组,水平方向分别为1.2km和0.6km,垂直方向分别为0.6km和0.3km;频率采用渐进式由低频到高频多尺度策略,从2Hz到20Hz分为16个频组。对实际模型进行高斯平滑,作为反演的初始模型。
表1理论模型多尺度全波形反演所采用的参数
由反演结果可得,本发明方法可以很好的反演出模型的构造和速度值,在水平方向速度变化剧烈的地方,反演效果也较理想。通过抽取模型上不同位置处的一维速度进行对比,可以看出,在深度2.5km以上地层中,反演结果能够很好的与实际数据吻合,反演速度相对于初始模型的速度有了很大的提高。
实施例3
将本发明提出的一种用于常规陆域地震数据的多尺度全波形反演方法应用于实际数据处理中,取得了理想的效果。
选取某地区陆域深反射地震数据40km,据测线900米位置处有一口地热井,井深2500米,作为反演结果的验证井。该地震数据激发震源为炸药震源,检波器为常规工业生产中使用的10Hz模拟检波器,道距20m,炮距100m,接收道数1440道,最小偏移距10m,最大偏移距14.4km。
对实际数据进行多尺度全波形反演,具体包括以下步骤:
步骤1,建立初始模型;
输入原始地震数据,获取原始地震数据的初至旅行时,基于二维射线追踪地震层析成像方法,确定层析成像结果,对层析成像结果进行周跳风险评价,确定全波形反演的速度模型,根据全波形反演模型的速度模型,确定全波形反演的密度模型,建立全波形反演模型。
步骤2,地震数据预处理;
首先,对原始地震进行常规地震数据处理,包括:静校正、去噪、球面振幅补偿、一致性校正、陷波、面波去除、初至以上切除,获得校正后的地震数据;再利用全波形反演的速度模型进行正演模拟得到正演模拟地震数据;由于正演模拟地震数据与校正后的地震数据两者的振幅不在同一尺度内,因此,需要对正演模拟地震数据与校正后的地震数据进行振幅匹配,匹配完成后,为了为Laplace-Fourier域多尺度全波形反演做准备,还需要将校正后的地震数据由时间域变换至频率域,获得多尺度全波形反演数据。
步骤3,多尺度全波形反演;
根据偏移距和频率从多尺度全波形反演数据内抽取数据,作为多尺度全波形反演的输入数据,进行多尺度全波形反演,多尺度全波形反演过程中所使用的参数如表2所示。
表2实际数据多尺度全波形反演所采用的参数
通过分别对三组不同高斯平滑因子进行测试获得反演结果,得到采用数值大的平滑参数进行平滑,反演结果更平滑但是局部缺少细节,采用数值小的平滑参数进行平滑,反演结果细节较多但是深层也加入了更多扰动,通过合理设置平滑参数,可以得到高分辨率的反演结果。
利用正演地震记录对比、叠前深度偏移分析和测井数据验证对多尺度全波形反演结果进行评价。
正演地震记录对比:通过正演炮点位置X=18km和20km处的地震记录,得到相比于初始模型地震记录,反演结果的地震记录在中浅层加入了丰富的反射波场信息,初至旅行时和中浅层反射信息与实际数据较一致。
叠前深度偏移分析:分别利用初始速度模型和全波形反演速度模型进行地震炮集的叠前深度偏移,对比后得到,反演速度模型的叠前深度偏移结果在深部成像效果明显好于初始速度模型的成像效果,并且根据偏移道集上可以得到,利用反演模型偏移的道集同相轴较利用初始模型偏移的道集更平,显示了反演的速度更准确,更接近于实际地层速度。
测井数据验证:利用D13地热井的测井数据验证反演模型,通过对比发现反演结果与测井曲线吻合的很好,尤其是在1100-1200m位置处,即测井解释的地热储层位置,同时,该位置对应电阻率曲线上的低阻异常位置。再对比D13井地温曲线与反演的速度曲线,得到在区域基岩界面以上(约900米),地温梯度为3.4℃/100m,符合该区域正常地温梯度,而在地热储层位置处(1100-1200米),地温梯度降到约1.0℃/100m,对应于速度的低阻异常区,符合该区地热储层物性规律。
因此,通过结合正演地震记录对比、叠前深度偏移分析和测井数据验证验证了本发明方法的应用效果,验证了本申请方法充分利用地震数据里的低频信息,有利于常规陆域地震数据的处理,为常规陆域的地质解释提供了依据。
当然,上述说明并非是对本发明的限制,本发明也并不仅限于上述举例,本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换,也应属于本发明的保护范围。
Claims (7)
1.一种用于常规陆域地震数据的多尺度全波形反演方法,其特征在于,具体包括如下步骤:
步骤1,输入原始地震数据,获取原始地震数据的初至旅行时,基于二维射线追踪地震层析成像方法,确定层析成像结果,对层析成像结果进行周跳风险评价,确定全波形反演的速度模型,根据全波形反演模型的速度模型,确定全波形反演的密度模型,建立全波形反演模型;
步骤2,对原始地震数据进行预处理,获得校正后的地震数据,利用全波形反演的速度模型进行正演模拟得到正演模拟地震数据,基于均方根振幅最小二乘拟合方法,对正演模拟地震数据与校正后的地震数据进行振幅匹配,并将校正后的地震数据由时间域变换至频率域,作为多尺度全波形反演数据;
步骤3,输入多尺度全波形反演数据和全波形反演模型,设置频率组、拉普拉斯常数组、偏移距组和迭代次数,基于层剥离的Laplace-Fourier域多尺度全波形反演方法,通过设置由近及远遍历偏移距组内偏移距的内部循环、由小到大遍历拉普拉斯常数组内拉普拉斯常数的中间循环以及由低到高的遍历频率组内频率的外部循环,利用全波形反演模型进行全波形反演,获得多尺度全波形反演结果,根据多尺度全波形反演结果进行地质解释。
2.如权利要求1所述的一种用于常规陆域地震数据的多尺度全波形反演方法,其特征在于,所述步骤1中,具体包括以下步骤:
步骤1.1,输入原始地震数据,拾取原始地震数据的初至旅行时获得实际旅行时,结合原始地震数据选取初至数据,设置层析成像参数;
步骤1.2,根据初至数据,基于二维射线追踪地震层析成像方法,得到层析成像结果,确定层析成像速度模型;
步骤1.3,对层析成像速度模型进行正演模拟,计算各炮点激发检波器位置处的初至旅行时,获得各炮点的正演旅行时,结合各炮点的实际旅行时,计算相同炮点和检波器位置处正演旅行时与实际旅行时之间的差值;
步骤1.4,设置周跳风险评价频率,计算评价周期,进行周跳风险评价,判断层析成像速度模型是否会产生周跳现象,若差值小于1/2个评价周期,则层析成像速度模型不会产生周跳现象,执行步骤1.6;若差值不小于1/2个评价周期,则层析成像速度模型会产生周跳现象,执行步骤1.5;
步骤1.5,选取层析成像速度模型中差值不小于1/2个评价周期的区域,增加炮点数据,继续进行正演模拟,返回步骤1.3;
步骤1.6,输出层析成像速度模型作为全波形反演模型中的速度模型,基于Gardner经验关系,确定全波形反演模型中的密度模型,设置品质因子模型,建立全波形反演模型。
5.如权利要求1所述的一种用于常规陆域地震数据的多尺度全波形反演方法,其特征在于,所述步骤2中,具体包括以下步骤:
步骤2.1,原始地震数据预处理,预处理包括静校正、去噪、去面波、振幅补偿、一次性校正和初至切除,得到校正后的地震数据;
步骤2.2,利用全波形反演的速度模型对校正后的地震数据进行正演模拟,得到正演模拟地震数据,基于均方根振幅最小二乘拟合方法,将正演模拟地震数据的振幅与校正后地震数据的振幅逐炮匹配,使得正演模拟地震数据与校正后地震数据的尺度相同;
正演模拟地震数据与校正后地震数据的振幅匹配公式为:
y=Cx (3)
式中,x为正演模拟地震数据的均方根振幅;y为校正后地震数据的均方根振幅;C为缩放系数;
步骤2.3,将校正后的地震数据由时间域变换至频率域,得到多尺度全波形反演数据;
校正后地震数据的变换公式如下:
式中,ω为角频率;t为时间;f(t)为校正后地震数据的时间域信号;F(ω)为校正后地震数据的频率域信号。
6.如权利要求1所述的一种用于常规陆域地震数据的多尺度全波形反演方法,其特征在于,所述步骤3中,具体包括以下步骤:
步骤3.1,输入多尺度全波形反演数据和全波形反演模型,分别设置频率组、拉普拉斯常数组、偏移距组和迭代次数,其中,频率组内按照频率由低到高的顺序设置有a个频率,拉普拉斯常数组内按照数值由小到大的顺序设置有b个拉普拉斯常数,偏移距组内按照数值由小到大的顺序设置有c个偏移距,分别对频率组、拉普拉斯常数组和偏移距组内各参数进行编号;
步骤3.2,分别设置频率、拉普拉斯常数和偏移距的初始值,其中,将频率组内的最低频率设置为频率的初始值,将拉普拉斯常数设置为拉普拉斯常数组内的最小值,将偏移距设置为偏移距组内的最小值;
步骤3.3,设置梯度吸收因子ε和高斯平滑因子ζ;
步骤3.4,根据当前设置的偏移距和频率从多尺度全波形反演数据中提取地震数据,利用全波形反演模型进行全波形反演,得到全波形反演结果,并将本次迭代计算得到的全波形反演结果设置为全波形反演模型;
步骤3.5,判断是否完成设置的迭代次数;
若未完成设置的迭代次数,则返回步骤3.4继续进行全波形反演;若已完成设置的迭代次数,则执行步骤3.6;
步骤3.6,判断是否遍历偏移距组内全部的偏移距;
若未遍历完,则根据偏移距组对当前偏移距Oi进行更新,i=1,…,c-1,i为当前偏移距在偏移距组内的序号,更新后的偏移距为Oi+1,返回步骤3.4继续进行全波形反演;若已遍历完,则对偏移距进行初始化,设置偏移距为O1,重新遍历偏移距组内的偏移距,执行步骤3.7;
步骤3.7,判断是否遍历完拉普拉斯常数组内全部的拉普拉斯常数;
若未遍历完,则根据拉普拉斯常数组对当前拉普拉斯常数Lj进行更新,j=1,…,b-1,j为当前拉普拉斯常数在拉普拉斯常数组内的序号,更新后的拉普拉斯常数为Lj+1,返回步骤3.4继续进行全波形反演;若已遍历完,则对拉普拉斯常数进行初始化,设置拉普拉斯常数为L1,重新遍历拉普拉斯常数组内的拉普拉斯常数,执行步骤3.8;
步骤3.8,判断是否遍历完频率组内全部的频率;
若未遍历完,则根据频率组对当前频率fn进行更新,n=1,…,a-1,n为当前频率在频率组内的序号,更新后的频率为fn+1,返回步骤3.4继续进行全波形反演;若已遍历完,则将当前全波形反演模型作为反演速度模型,执行步骤3.9;
步骤3.9,输出多尺度全波形反演结果,并根据多尺度全波形反演结果进行地质解释。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110683680.1A CN113552625B (zh) | 2021-06-21 | 2021-06-21 | 一种用于常规陆域地震数据的多尺度全波形反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110683680.1A CN113552625B (zh) | 2021-06-21 | 2021-06-21 | 一种用于常规陆域地震数据的多尺度全波形反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113552625A true CN113552625A (zh) | 2021-10-26 |
CN113552625B CN113552625B (zh) | 2022-05-13 |
Family
ID=78130747
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110683680.1A Active CN113552625B (zh) | 2021-06-21 | 2021-06-21 | 一种用于常规陆域地震数据的多尺度全波形反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113552625B (zh) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116327250A (zh) * | 2023-02-13 | 2023-06-27 | 中国科学院地质与地球物理研究所 | 一种基于全波形反演技术的乳腺超声三维成像方法 |
CN116540297A (zh) * | 2023-05-06 | 2023-08-04 | 中国科学院地质与地球物理研究所 | 一种弹性波地震数据全波形反演方法、系统及设备 |
CN116699695A (zh) * | 2023-08-07 | 2023-09-05 | 北京中矿大地地球探测工程技术有限公司 | 一种基于衰减矫正的反演方法、装置和设备 |
WO2023230946A1 (zh) * | 2022-05-30 | 2023-12-07 | 佟小龙 | 陆地地球物理勘探方法、电子设备及可读存储介质 |
WO2024078134A1 (zh) * | 2022-10-13 | 2024-04-18 | 安徽理工大学 | 基于多参数约束和结构校正的随掘巷道全波形反演方法 |
Citations (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2012170091A1 (en) * | 2011-06-08 | 2012-12-13 | Chevron U.S.A. Inc. | System and method for seismic data inversion |
CN105549080A (zh) * | 2016-01-20 | 2016-05-04 | 中国石油大学(华东) | 一种基于辅助坐标系的起伏地表波形反演方法 |
CN105676277A (zh) * | 2015-12-30 | 2016-06-15 | 中国石油大学(华东) | 一种提高高陡构造速度反演效率的全波形联合反演方法 |
CN105891888A (zh) * | 2016-03-28 | 2016-08-24 | 吉林大学 | 多域分频并行多尺度全波形反演方法 |
CN106054244A (zh) * | 2016-06-16 | 2016-10-26 | 吉林大学 | 截断时窗的低通滤波多尺度全波形反演方法 |
CN107422379A (zh) * | 2017-07-27 | 2017-12-01 | 中国海洋石油总公司 | 基于局部自适应凸化方法的多尺度地震全波形反演方法 |
CN107765302A (zh) * | 2017-10-20 | 2018-03-06 | 吉林大学 | 不依赖震源子波的时间域单频波形走时反演方法 |
CN108254731A (zh) * | 2018-04-25 | 2018-07-06 | 吉林大学 | 探地雷达数据的多尺度阶梯式层剥离全波形反演方法 |
US20180356548A1 (en) * | 2017-06-12 | 2018-12-13 | Institute Of Geology And Geophysics Chinese Academy Of Sciences | Inversion velocity model, method for establishing the same and method for acquiring images of underground structure |
CN110058302A (zh) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | 一种基于预条件共轭梯度加速算法的全波形反演方法 |
WO2019234469A1 (en) * | 2018-06-08 | 2019-12-12 | Total Sa | Method for generating an image of a subsurface of an area of interest from seismic data |
CN111290016A (zh) * | 2020-03-04 | 2020-06-16 | 中国石油大学(华东) | 一种基于地质模型约束的全波形速度建模反演方法 |
-
2021
- 2021-06-21 CN CN202110683680.1A patent/CN113552625B/zh active Active
Patent Citations (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2012170091A1 (en) * | 2011-06-08 | 2012-12-13 | Chevron U.S.A. Inc. | System and method for seismic data inversion |
CN105676277A (zh) * | 2015-12-30 | 2016-06-15 | 中国石油大学(华东) | 一种提高高陡构造速度反演效率的全波形联合反演方法 |
CN105549080A (zh) * | 2016-01-20 | 2016-05-04 | 中国石油大学(华东) | 一种基于辅助坐标系的起伏地表波形反演方法 |
CN105891888A (zh) * | 2016-03-28 | 2016-08-24 | 吉林大学 | 多域分频并行多尺度全波形反演方法 |
CN106054244A (zh) * | 2016-06-16 | 2016-10-26 | 吉林大学 | 截断时窗的低通滤波多尺度全波形反演方法 |
US20180356548A1 (en) * | 2017-06-12 | 2018-12-13 | Institute Of Geology And Geophysics Chinese Academy Of Sciences | Inversion velocity model, method for establishing the same and method for acquiring images of underground structure |
CN107422379A (zh) * | 2017-07-27 | 2017-12-01 | 中国海洋石油总公司 | 基于局部自适应凸化方法的多尺度地震全波形反演方法 |
CN107765302A (zh) * | 2017-10-20 | 2018-03-06 | 吉林大学 | 不依赖震源子波的时间域单频波形走时反演方法 |
CN108254731A (zh) * | 2018-04-25 | 2018-07-06 | 吉林大学 | 探地雷达数据的多尺度阶梯式层剥离全波形反演方法 |
WO2019234469A1 (en) * | 2018-06-08 | 2019-12-12 | Total Sa | Method for generating an image of a subsurface of an area of interest from seismic data |
CN110058302A (zh) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | 一种基于预条件共轭梯度加速算法的全波形反演方法 |
CN111290016A (zh) * | 2020-03-04 | 2020-06-16 | 中国石油大学(华东) | 一种基于地质模型约束的全波形速度建模反演方法 |
Non-Patent Citations (2)
Title |
---|
包乾宗等: "基于地震数据包络的多尺度全波形反演方法", 《石油物探》 * |
胡英等: "Laplace-Fourier域多尺度高效全波形反演方法", 《石油勘探与开发》 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2023230946A1 (zh) * | 2022-05-30 | 2023-12-07 | 佟小龙 | 陆地地球物理勘探方法、电子设备及可读存储介质 |
WO2024078134A1 (zh) * | 2022-10-13 | 2024-04-18 | 安徽理工大学 | 基于多参数约束和结构校正的随掘巷道全波形反演方法 |
CN116327250A (zh) * | 2023-02-13 | 2023-06-27 | 中国科学院地质与地球物理研究所 | 一种基于全波形反演技术的乳腺超声三维成像方法 |
CN116327250B (zh) * | 2023-02-13 | 2023-08-25 | 中国科学院地质与地球物理研究所 | 一种基于全波形反演技术的乳腺超声三维成像方法 |
CN116540297A (zh) * | 2023-05-06 | 2023-08-04 | 中国科学院地质与地球物理研究所 | 一种弹性波地震数据全波形反演方法、系统及设备 |
CN116540297B (zh) * | 2023-05-06 | 2024-03-08 | 中国科学院地质与地球物理研究所 | 一种弹性波地震数据全波形反演方法、系统及设备 |
CN116699695A (zh) * | 2023-08-07 | 2023-09-05 | 北京中矿大地地球探测工程技术有限公司 | 一种基于衰减矫正的反演方法、装置和设备 |
CN116699695B (zh) * | 2023-08-07 | 2023-11-03 | 北京中矿大地地球探测工程技术有限公司 | 一种基于衰减矫正的反演方法、装置和设备 |
Also Published As
Publication number | Publication date |
---|---|
CN113552625B (zh) | 2022-05-13 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113552625B (zh) | 一种用于常规陆域地震数据的多尺度全波形反演方法 | |
CN106405651B (zh) | 一种基于测井匹配的全波形反演初始速度模型构建方法 | |
CN105277978B (zh) | 一种确定近地表速度模型的方法及装置 | |
CN106526674A (zh) | 一种三维全波形反演能量加权梯度预处理方法 | |
CN106932819B (zh) | 基于各向异性马尔科夫随机域的叠前地震参数反演方法 | |
CN109884710B (zh) | 针对激发井深设计的微测井层析成像方法 | |
CN105319589B (zh) | 一种利用局部同相轴斜率的全自动立体层析反演方法 | |
CN108254780A (zh) | 一种微地震定位及各向异性速度结构层析成像方法 | |
CN109917454B (zh) | 基于双基准面的真地表叠前深度偏移成像方法及装置 | |
CN109581501B (zh) | 用于沙漠区深度域速度建模的方法 | |
CN103616723B (zh) | 基于avo特征的crp道集真振幅恢复方法 | |
CN102901985A (zh) | 一种适用于起伏地表的深度域层速度修正方法 | |
CN110531410A (zh) | 一种基于直达波场的最小二乘逆时偏移梯度预条件方法 | |
CN109884700A (zh) | 多信息融合地震速度建模方法 | |
CN105717538B (zh) | 起伏地表地震数据偏移基准面转换方法及装置 | |
CN111399037B (zh) | 高速顶界面提取的方法和装置 | |
CN109490961B (zh) | 起伏地表无射线追踪回折波层析成像方法 | |
CN109725354B (zh) | 各向异性速度建模方法及系统 | |
CN111175822B (zh) | 改进直接包络反演与扰动分解的强散射介质反演方法 | |
CN112305595B (zh) | 基于折射波分析地质体结构的方法及存储介质 | |
CN109901221B (zh) | 一种基于动校正速度参数的地震资料各向异性建模方法 | |
CN112415601A (zh) | 表层品质因子q值的确定方法及装置 | |
CN112649876A (zh) | 建立地震偏移速度模型的方法和装置 | |
CN113267808B (zh) | 振幅补偿方法及装置 | |
Zhou et al. | Anisotropic model building with well control |
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 |