CN114578452B - 一种定量计算地下反射系数振幅占比的方法 - Google Patents

一种定量计算地下反射系数振幅占比的方法 Download PDF

Info

Publication number
CN114578452B
CN114578452B CN202210489190.2A CN202210489190A CN114578452B CN 114578452 B CN114578452 B CN 114578452B CN 202210489190 A CN202210489190 A CN 202210489190A CN 114578452 B CN114578452 B CN 114578452B
Authority
CN
China
Prior art keywords
curve
reflection coefficient
syn
peak
trough
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.)
Active
Application number
CN202210489190.2A
Other languages
English (en)
Other versions
CN114578452A (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.)
Dida Huineng Beijing Technology Co ltd
Original Assignee
Dida Huineng Beijing Technology Co ltd
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 Dida Huineng Beijing Technology Co ltd filed Critical Dida Huineng Beijing Technology Co ltd
Priority to CN202210489190.2A priority Critical patent/CN114578452B/zh
Publication of CN114578452A publication Critical patent/CN114578452A/zh
Application granted granted Critical
Publication of CN114578452B publication Critical patent/CN114578452B/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
    • G01V11/00Prospecting or detecting by methods combining techniques covered by two or more of main groups G01V1/00 - G01V9/00

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及石油勘探开采技术领域,具体涉及一种定量计算地下反射系数振幅占比的方法,通过在常规合成记录标定原始地震与合成记录相关性的基础上,还同时标定出井的地层反射系数序列,输出各反射系数对指定的波峰或波谷的振幅占比,实现了地震波形与反射系数组成和振幅占比的定量关系,标定出引起波峰或波谷反射的主要反射系数位置以及主要的干扰层位置,实现复杂储层地震响应特征的定量分析。本方法特别适用于各种岩性组成的薄互层、存在强背景干扰层、弱反射层等地区的地震勘探。

Description

一种定量计算地下反射系数振幅占比的方法
技术领域
本发明涉及石油勘探开采技术领域,尤其涉及一种定量计算地下反射系数振幅占比的方法。
背景技术
合成地震记录是地震模型技术中应用非常广泛的一种,其精度直接影响到地质层位的准确标定和地质体精细解释。
油气勘探工作越来越向隐蔽性油气藏发展,目标尺度越来越小,对合成地震记录提出了更高的要求。随着油气勘探目标逐步转向低幅度构造和薄层,地震解释人员通过地震剖面中的同相轴认识地下构造层的难度也随之增加。对于精细保幅地震处理、构造解释和精细岩性反演而言,同相轴与构造层的错误匹配将导致不正确甚至错误的地质认识和结论,因此,地震正演模型的研究工作就显得尤为基础和重要。
传统的合成地震记录的制作是一个简化的一维正演的过程,由声波和密度测井曲线计算得到反射系数,将反射系数与提取的地震子波进行褶积得到初始合成地震记录。根据较精确的速度场对初始合成地震记录进行校正,再与井旁地震道匹配调整,得到最终合成地震记录。
上述合成记录制作方法存在的主要问题是不能确定反射系数与地震反射的定量关系,一般只用于时深转换。老油区随着勘探的深入,勘探目标的会尺度越来越小,最容易出现的问题是错误的将一些薄层砂体只通过时深对应就认为某某同相轴是该砂体的反射而进行追踪解释,没有考虑到地震波形通常是地下多层反射的共同响应造成的,对合成地震记录提出了更高的要求,如何定量分析反射界面组成与地震波形的对应关系,需要新的技术方法。
发明内容
有鉴于此,本发明的目的在于提出一种定量计算地下反射系数振幅占比的方法,以解决需要定量分析反射界面组成与地震波形的对应关系的问题。
基于上述目的,本发明提供了一种定量计算地下反射系数振幅占比的方法,包括以下步骤:
获得研究区过探井位置的地震数据,以及该探井的声波时差测井曲线和密度测井曲线;
根据声波时差测井曲线和密度测井曲线计算波阻抗深度曲线;
对声波时差测井曲线、密度测井曲线和波阻抗时间曲线进行深度域到时间域的转换,得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线;
提取波阻抗的时间域曲线的反射系数序列;
基于过探井位置的地震数据得到雷克子波曲线,将雷克子波曲线与反射系数序列褶积得到合成地震记录;
计算合成地震记录的波峰波谷位置,并进行排序;
选择指定的波峰或波谷,根据反射系数序列获取其反射系数构成,计算该波峰或波谷前后设定范围内的各反射系数在该波峰或波谷处的振幅响应值;
计算各反射系数在指定的波峰或波谷时刻的振幅响应值与该波峰或波谷的总振幅的比值,输出为各反射系数对指定的波峰或波谷的振幅占比。
优选地,根据声波时差测井曲线和密度测井曲线计算波阻抗深度曲线包括:
将声波时差测井曲线和密度测井曲线的每一个深度采样点换算得到声波速度曲线vi和密度曲线ρi,利用公式Zsi=ρi×vi,(i=1,2……,nn)计算得到波阻抗深度曲线,式中Zsi表示波阻抗深度曲线,i表示采样点,nn表示测井曲线有nn个采样点。
优选地,对声波时差测井曲线、密度测井曲线和波阻抗时间曲线进行深度域到时间域的转换,得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线包括:
利用声波速度曲线并进行时间重采样,得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线Zi(i=1,2……n)。
优选地,提取波阻抗的时间域曲线的反射系数序列包括:
利用反射系数计算公式Ri=(Zi+1-Zi)/(Zi+1+Zi),(i=1,2……,n)将波阻抗的时间域曲线Zi换算得到时间域反射系数曲线Ri,从而获得反射系数序列,式中Ri表示反射系数曲线,Zi表示波阻抗的时间域曲线,i表示采样点。
优选地,本方法还包括:
设置阈值ε1,拾取大于ε1的Ri得到主要反射系数Rk(k=1,2……,m),并按照局部极大和极小处理进行方波化。
优选地,过探井位置的地震数据包括井旁地震道,基于过探井位置的地震数据得到雷克子波曲线,将雷克子波曲线与反射系数序列褶积得到合成地震记录包括:
输入与井旁地震道seisdata一致的子波频率、长度、极性和地震采样率,采用28HZ雷克子波,得到雷克子波曲线wi(i=1,2……n1),并进行显示,式中i为采样点;
采用褶积公式得到合成地震记录syn=wi*Ri,并与井旁地震道seisdata进行相关性比对,修改子波参数,直到相关性匹配结果达到设定值,其中Ri为时间域反射系数曲线。
优选地,计算合成地震记录的波峰波谷位置,并进行排序包括:
依次拾取syn(i-1)、syn(i)和syn(i+1),如果syn(i)>syn(i-1),且syn(i)>syn(i+1),且|syn(i)-syn(i-1)|>ε2,则标定为peak(1),依次拾取到合成记录最有一个采样点,即得到离散的自上而下波峰振幅值peak(kk)(kk=1,2……M);
依次拾取syn(i-1)、syn(i)和syn(i+1),如果syn(i)<syn(i-1),且syn(i)<syn(i+1),且|syn(i)-syn(i-1)|>ε2,得到离散的波谷振幅值lough(kk);
式中,kk表示波峰数或波谷数,ε2为阈值。
优选地,选择指定的波峰或波谷,根据反射系数序列获取其反射系数构成,计算该波峰或波谷前后设定范围内的各反射系数在该波峰或波谷处的振幅响应值包括:
输入指定的波峰或波谷的序列号kk,拾取该波峰或波谷前后半个子波长度范围内的反射系数Rk'(k=1,2……N1),N1为该范围内的反射系数个数,分别对Rk'与输入的地震子波进行褶积,并计算在指定的波峰或波谷对应的时间处其振幅响应值,得到ZZ(i=1,2,……N1),ZZ为各反射系数在该处的振幅响应值。
优选地,振幅占比的计算公式为:
fsxszb(i)=ZZi/∑ZZi
其中,fsxszb表示振幅占比,ZZi为各反射系数在指定波峰或波谷处的振幅响应值,i=1,2……N1。
优选地,本方法还包括:
输出各反射系数对指定的波峰或波谷的振幅占比图。
本发明的有益效果:通过在常规合成记录标定原始地震与合成记录相关性的基础上,还同时标定出井的地层反射系数序列,输出各反射系数对指定的波峰或波谷的振幅占比,实现了地震波形与反射系数组成和振幅占比的定量关系,标定出引起波峰或波谷反射的主要反射系数位置以及主要的干扰层位置,实现复杂储层地震响应特征的定量分析。本方法特别适用于各种岩性组成的薄互层、存在强背景干扰层、弱反射层等地区的地震勘探。
附图说明
为了更清楚地说明本发明或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例的定量计算地下反射系数振幅占比的方法流程示意图;
图2为本发明实施例的过草322井的南北向地震剖面图;
图3为本发明实施例的原始声波时差曲线、原始密度测井曲线、声波时差的时间域曲线和密度的时间域示意图;
图4为本发明实施例的全井段反射系数序列及合成地震记录图;
图5为本发明实施例的全井段合成地震记录及对波峰响应占比最大的反射系数构成图;
图6为本发明实施例的第1到5号波峰的反射系数振幅占比图;
图7为本发明实施例的1-5个地震同相轴的反射系数振幅占比图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白,以下结合具体实施例,对本发明进一步详细说明。
需要说明的是,除非另外定义,本发明使用的技术术语或者科学术语应当为本发明所属领域内具有一般技能的人士所理解的通常意义。本发明中使用的“第一”、“第二”以及类似的词语并不表示任何顺序、数量或者重要性,而只是用来区分不同的组成部分。“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同,而不排除其他元件或者物件。“连接”或者“相连”等类似的词语并非限定于物理的或者机械的连接,而是可以包括电性的连接,不管是直接的还是间接的。“上”、“下”、“左”、“右”等仅用于表示相对位置关系,当被描述对象的绝对位置改变后,则该相对位置关系也可能相应地改变。
如图1所示,本说明书实施例提供一种定量计算地下反射系数振幅占比的方法,包括以下步骤:
S101获得研究区过探井位置的地震数据,以及该探井的声波时差测井曲线和密度测井曲线。
S102根据声波时差测井曲线和密度测井曲线计算波阻抗深度曲线;
举例来说,测井曲线含有nn个采样点,深度采样波长为dz,该步骤具体包括将声波时差测井曲线和密度测井曲线的每一个深度采样点换算得到声波速度曲线vi和密度曲线ρi,其中声波速度vi即为声波时差的倒数vi=1/dti,利用公式Zsi=ρi×vi,(i=1,2……,nn)计算得到波阻抗深度曲线Zsi(i=1,2……,nn),式中Zsi表示波阻抗深度曲线,i表示采样点,此处需要给定测井曲线的记录的格式,选择米制还是英制。
S103对声波时差测井曲线、密度测井曲线和波阻抗时间曲线进行深度域到时间域的转换,得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线。
具体来说是利用声波速度曲线并进行时间重采样,得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线Zi(i=1,2……n)。
S104提取波阻抗的时间域曲线的反射系数序列。
具体来说是利用反射系数计算公式Ri=(Zi+1-Zi)/(Zi+1+Zi),(i=1,2……,n)将波阻抗的时间域曲线Zi换算得到时间域反射系数曲线Ri,从而获得反射系数序列,式中Ri表示反射系数曲线,Zi表示波阻抗的时间域曲线。
举例来说,本方法在该步骤后还包括:
设置阈值ε1,拾取大于ε1的Ri得到主要反射系数Rk(k=1,2……,m),并按照局部极大和极小处理进行方波化。
S105基于过探井位置的地震数据得到雷克子波曲线,将雷克子波曲线与反射系数序列褶积得到合成地震记录。
举例来说,本步骤具体包括:
输入与井旁地震道seisdata一致的子波频率、长度、极性和地震采样率,采用28HZ雷克子波,得到雷克子波曲线wi(i=1,2……n1),并进行显示,式中i为采样点;
采用褶积公式得到合成地震记录syn=wi*Ri,并与井旁地震道seisdata进行相关性比对,修改子波参数,直到相关性匹配结果达到设定值,即匹配到满意为止。
S106计算合成地震记录的波峰波谷位置,并进行排序。
举例来说,本步骤进一步包括:
依次拾取syn(i-1)、syn(i)和syn(i+1),如果syn(i)>syn(i-1),且syn(i)>syn(i+1),且|syn(i)-syn(i-1)|>ε2,则标定为peak(1),依次拾取到合成记录最后一个采样点,即得到离散的自上而下波峰振幅值peak(kk)(kk=1,2……M);
依次拾取syn(i-1)、syn(i)和syn(i+1),如果syn(i)<syn(i-1),且syn(i)<syn(i+1),且|syn(i)-syn(i-1)|>ε2,得到离散的波谷振幅值lough(kk);
式中,kk表示波峰数或波谷数,ε2为阈值,用于剔除微小的扰动。
S107选择指定的波峰或波谷,即需要计算的波峰或波谷,根据反射系数序列获取其反射系数构成,计算该波峰或波谷前后设定范围内的各反射系数在该波峰或波谷处的振幅响应值。
举例来说,本步骤具体包括:
输入指定的波峰或波谷的序列号kk,拾取该波峰或波谷前后半个子波长度范围内的反射系数Rk'(k=1,2……N1),N1为该范围内的反射系数个数,分别对Rk'与输入的地震子波进行褶积,并计算在指定的波峰或波谷对应的时间处其振幅响应值,得到ZZ(i=1,2,……N1),ZZ为各反射系数在该处的振幅响应值。
如假设kk=2.则计算在kk=2时波峰对应的时间处其振幅响应值,得到ZZ(i=1,2,……N1),ZZ为各反射系数在kk=2处的振幅响应值。
S108计算各反射系数在指定的波峰或波谷时刻的振幅响应值与该波峰或波谷的总振幅的比值,输出为各反射系数对指定的波峰或波谷的振幅占比。
举例来说,是通过计算fsxszb(i)=ZZi/∑ZZi,fsxszb即各反射系数对波峰的振幅占比,用于指示各反射系数对波峰形成的贡献,式中i=1,2……N1。
通过反射系数振幅占比计算。ZZ有正有负,对于波峰,取Zmax_location=t(max(ZZ)),即最大振幅占比对应的反射系数位置,用来指示对该地震波形贡献最大的地层;Zmin_location=t(min(ZZ)),即为最大负振幅占比对应的反射系数位置,通常用来指示最大的干扰地层。也可以指定任意波峰或波谷,输出各反射系数对波峰或波谷的振幅占比图。
下面以以济阳坳陷广饶地区为例,介绍本方法的实现过程和取得的地质认识。
济阳坳陷广饶地区在新近系馆陶组底部广泛分布有厚度在10-20m的底砾岩,底砾岩之上上伏有厚度不等的泥岩和多层系的火山玄武岩。地震反射剖面上底砾岩呈连续强反射,底砾岩之上的火山岩也呈强反射,在两个强反射之间,发育有短轴透镜状中-弱反射,如图2所示。以往认为这些反射可能是河道砂岩的反射,为此专门部署了探井草322井失利,未见河道砂。为开展失利井分析,采用提出的基于计算地下反射界面地震响应占比的合成地震记录制作方法,对该地区各种反射的成因进行了分析,解决了储层描述不合理的问题。
具体实现过程以草322井为例可以总结为以下几个步骤:
1.首先实现传统合成记录制作过程:
(1)数据加载:读取研究区过草322井地震成像数据道,以及草322井的声波时差测井曲线和密度测井曲线。
(2)计算波阻抗深度曲线:这里测井曲线含有nn个测量点,深度采样步长为dz=0.125;并将声波时差测井曲线和密度测井曲线每一个深度样点换算得到声波速度曲线vi和密度曲线ρi,利用公式Zsi=ρi×vi。得到波阻抗深度曲线,Zsi(i=1,2……,nn)。
(3)进行深度域到时间域的转换,方法是利用速度曲线并进行时间重采样,本次采用1ms的采样间隔,可以得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线Zi(i=1,2……n)。如图3所示分别为原始声波时差测井曲线和密度测井曲线和采用上述方法得到的对应的时间域曲线。
(4)提取反射系数序列:利用反射系数计算公式Ri=(Zi+1-Zi)/(Zi+1+Zi),(i=1,2……,n)计算反射系数,为避免一些极弱的可以忽略不计的反射系数影响,设置阈值ε1=0.001,拾取大于ε1的Ri得到主要反射系数Rk(k=1,2……,m)。同时按照局部极大和极小判识法进行方波化。
(5)合成记录制作:输入与井旁地震道一致的子波,设置频率、长度、极性和地震采样率等参数,本次根据广饶地区实际,采用28Hz,长度60ms、正极性、0.001ms采样率的雷克子波,得到雷克子波曲线wi(i=1,2……n1);与反射系数进行褶积得到合成地震记录syn。并与井旁地震道进行相关性比对,如相关差,可以修改子波参数,直到匹配满意为止。如图4为得到的该井反射系数曲线与合成记录,与原始地震相关性达到85%。
2.实现地层反射系数、最大振幅占比的反射系数、各反射系数振幅占比的标定:
(1)自动计算合成记录波峰或波谷的位置,并排序。
实际工作中,需要对所有感兴趣的地震同相轴的进行分析,为此需要对地震的波峰或波谷进行标识,以实现机器的自动识别。
方法是采用极大极小判识法,依次拾取syn(i-1)、syn(i)和syn(i+1),如果syn(i)>syn(i-1),且syn(i)>syn(i+1),且|syn(i)-syn(i-1)|>ε22为阈值,用于剔除微小的扰动,这里取值0.00001),则可以标定为peak(1),依次拾取到合成记录最后一个样点,即可以得到离散的自上而下波峰振幅值peak(kk)(kk=1,2……M),kk为总的波峰数。同理依次拾取syn(i-1)、syn(i)和syn(i+1),如果syn(i)<syn(i-1),且syn(i)<syn(i+1),且|syn(i)-syn(i-1)|>ε2,可以得到离散的波谷振幅值lough(kk)。本次工作只对波峰进行了标识。
(2)离散合成记录的制作与振幅占比的提取。
计算每一个反射系数在对应波峰处的振幅响应值并计算振幅占比,是本方法创新的核心和基础。
选择要计算的波峰或波谷,观察其反射系数构成。方法是输入kk,以第9个波峰为例,输入kk=9,即获得第9个波峰处地震响应的构成。拾取在第9个波峰前后半个子波长度范围内的反射系数,因为只有在此范围内的反射系数才能对波峰的形成有影响,设定为Rk'(k=1,2……N1),N1为该范围内的反射系数个数,分别对Ri采用褶积公式并计算在kk=9时波峰对应的时间处其振幅响应值,得到ZZ(i=1,2,……N1),ZZ为各反射系数在kk=9处的振幅响应值。如图5展示了反射系数的振幅取值方法,自上而下通过对每一个反射系数Rk(k=1,2……N1)与子波的褶积,得到地震波形,并计算该第二个波峰对应时间的振幅取值,得到ZZk(k=1,2,……N1)。
进一步,计算每一个反射系数在指定波峰时刻的振幅与波峰总振幅的比值,即fsxszb(i)=ZZi/∑ZZi,这里定义为振幅占比,它反映了反射系数对地震波形振幅贡献的大小,式中i=1,2……N1。
以第9个波峰为例,如表1是振幅占比的程序计算结果,第一列为反射系数的位置,第二列为第9个波峰,时间位置为235ms出的振幅响应,第三列为总的波峰振幅,为恒定值,相除,得到振幅占比。
表1振幅占比的程序计算结果
Figure GDA0003771531270000111
Figure GDA0003771531270000121
(3)最大反射系数振幅占比计算。
通过拾取振幅占比fsxszb的最大值即可获得对波峰影响最大的反射系数和位置。实际拾取过程中,fsxszb有正有负,对于波峰,取Zmax_location=t(max(fsxszb)),即最大振幅占比对应的反射系数位置,用来指示对地震波形贡献最大的地层反射系数;取Zmin_location=t(min(fsxszb)),即为最大负振幅占比对应的反射系数位置,通常用来指示最大的干扰地层。如上表中,对波峰9起最大作用的233ms处的反射系数,振幅占比值为1.5826。
如图6展示了该井各波峰起决定作用的反射系数及其位置,自左向右,分别为综合录井、深度域声波曲线、时间域声波曲线、反射系数序列和最大振幅占比反射系数(带有数值的反射系数)、合成地震记录和原始地震。图中,合成记录上的数字序号为计算得到的波峰序列,带有数值的反射系数为对应这些波峰形成的振幅贡献最大的反射系数,其旁的数值为占整个波峰的振幅占比,最大振幅占比的反射系数位置同时也以数字序号的形式标定在时间域和深度域曲线上。从图中可以看出原本认为是河道砂体的(4)号波峰,真实的情况是由第四套火成岩底界面而形成,非河道砂体的反射,从而避免了下一步的认识错误。
第四例带有数值的反射系数为振幅占比最大的。1到5号波峰,起最大作用的反射系数分别在42、65、90、113、139ms处,占比值分别为0.8、3.1、1.1、38.9和1。
从图上可以看出,对于波峰振幅的大小,并不是反射系数大的起的作用就大,还与地质体的厚度有关。波峰3是多层高速火成岩的综合响应;呈透镜状反射的波峰4,其信号产生的根源是火成岩第四套火成岩的反射系数形成,并不是河道砂体的反射,从而避免了进一步勘探的失误;而波峰5更多的是底砾岩底界面负反射形成的旁瓣而形成。
(4)各反射系数对波峰影响的振幅占比的计算
如果对指定的地震同相轴进行分析,可以输入kk,如kk=2,表示自上而下的第2个波峰处,输出fsxszb即各反射系数对波峰的振幅占比,可以用来指示对地震波形影响的相对大小。
借助fsxszb,可以对形成地震波峰的原因和信号的类型进行定量的分析。
如图7分别为1-5个地震同相轴的反射系数振幅占比。自上而下分别为合成地震记录(其上的数字为波峰自上而下的序列号)、第1至第5波峰反射系数的振幅占比图。反射系数定量的振幅占比图的横坐标为时间,纵坐标为各反射系数的振幅占比,大于0为贡献率为正,小于0贡献率为负,高低指示了占比值的大小,从反射系数占比图中可以看出引起地震响应的信号复杂程度,1号波峰为相对简单的正反射系数单信号源组成,2和3号波峰表现为多个正反射系数组成的复杂信号源,而4-5号波峰,负反射系数占比很大,代表了存在强负反射干扰的信号源。
在占比图上纵坐标为占比值,向上(值为正)表示对波峰起正作用,向下(值为负)表示对波峰起负作用。从其反射系数对波峰振幅占比可以看出,1号波峰存在一个占绝对主力的反射层,可以认为是相对简单的地震信号源;3号波峰,表现为一个较强和几个相对较弱的组合,可以认为是多层火成岩组成的调谐信号源;4号波峰的反射系数振幅占比图最为复杂,没有起决定作用的反射系数,起负作用的反射系数作用最大,说明调谐作用最为复杂,而本身波峰振幅也弱,地震解释起来难度最大,需要仔细开展信号源的分析。通过上述分析,反射系数振幅占比图可以为岩性组合及其反射特征进行定量的分析,并在随后的处理中开展针对性去干扰或弱信号增强处理。
3.成果输出
在传统的合成记录,即合成记录+地震剖面基础上,同时输出最大振幅占比反射系数序列(如图6),指定波峰或波谷时的反射系数构成和振幅占比图(如图7)。
综上,本方法在传统标定的基础上,同时标定出地层反射系数序列、通过振幅占比标定对地震波形影响最大的反射系数、以及各反射系数对地震波形影响的振幅占比图。利用占比最大的反射系数,便于我们准确的追踪解释储层,利用反射系数振幅占比图,可以有效分析信号源的复杂程度。
所属领域的普通技术人员应当理解:以上任何实施例的讨论仅为示例性的,并非旨在暗示本发明的范围(包括权利要求)被限于这些例子;在本发明的思路下,以上实施例或者不同实施例中的技术特征之间也可以进行组合,步骤可以以任意顺序实现,并存在如上所述的本发明的不同方面的许多其它变化,为了简明它们没有在细节中提供。
本发明旨在涵盖落入所附权利要求的宽泛范围之内的所有这样的替换、修改和变型。因此,凡在本发明的精神和原则之内,所做的任何省略、修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (8)

1.一种定量计算地下反射系数振幅占比的方法,其特征在于,包括以下步骤:
获得研究区过探井位置的地震数据,以及该探井的声波时差测井曲线和密度测井曲线;
根据声波时差测井曲线和密度测井曲线计算波阻抗深度曲线;
对声波时差测井曲线、密度测井曲线和波阻抗时间曲线进行深度域到时间域的转换,得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线;
提取波阻抗时间域曲线的反射系数序列;
基于过探井位置的地震数据得到雷克子波曲线,将雷克子波曲线与反射系数序列褶积得到合成地震记录;
计算合成地震记录的波峰波谷位置,并进行排序;
选择指定的波峰或波谷,根据反射系数序列获取其反射系数构成,计算该波峰或波谷前后设定范围内的各反射系数在该波峰或波谷处的振幅响应值;
计算各反射系数在指定的波峰或波谷时刻的振幅响应值与该波峰或波谷的总振幅的比值,输出为各反射系数对指定的波峰或波谷的振幅占比;
所述过探井位置的地震数据为井旁地震道,所述基于过探井位置的地震数据得到雷克子波曲线,将雷克子波曲线与反射系数序列褶积得到合成地震记录包括:
输入与井旁地震道seisdata一致的子波频率、长度、极性和地震采样率,采用28HZ雷克子波,得到雷克子波曲线wi,并进行显示,式中i为采样点,i=1,2……n1;
采用褶积公式得到合成地震记录syn=wi*Ri,并与井旁地震道seisdata进行相关性比对,修改子波参数,直到相关性匹配结果达到设定值,其中Ri为时间域反射系数曲线;
所述计算合成地震记录的波峰波谷位置,并进行排序包括:
依次拾取syn(i-1)、syn(i)和syn(i+1),如果syn(i)>syn(i-1),且syn(i)>syn(i+1),且|syn(i)-syn(i-1)|>ε2,则标定为peak(1),依次拾取到合成记录最后一个样点,即得到离散的自上而下波峰振幅值peak(kk);
依次拾取syn(i-1)、syn(i)和syn(i+1),如果syn(i)<syn(i-1),且syn(i)<syn(i+1),且|syn(i)-syn(i-1)|>ε2,得到离散的波谷振幅值lough(kk);
式中,kk表示波峰数或波谷数,kk=1,2……M,ε2为阈值。
2.根据权利要求1所述的定量计算地下反射系数振幅占比的方法,其特征在于,所述根据声波时差测井曲线和密度测井曲线计算波阻抗深度曲线包括:
将声波时差测井曲线和密度测井曲线的每一个深度采样点换算得到声波速度曲线vi和密度曲线ρi,利用公式Zsi=ρi×vi,计算得到波阻抗深度曲线,式中Zsi表示波阻抗深度曲线,i表示采样点,i=1,2……,nn,nn表示测井曲线有nn个采样点。
3.根据权利要求2所述的定量计算地下反射系数振幅占比的方法,其特征在于,所述对声波时差测井曲线、密度测井曲线和波阻抗时间曲线进行深度域到时间域的转换,得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线包括:
利用声波速度曲线并进行时间重采样,得到声波时差的时间域曲线、密度的时间域曲线和波阻抗的时间域曲线Zi,其中i=1,2……n。
4.根据权利要求1所述的定量计算地下反射系数振幅占比的方法,其特征在于,所述提取波阻抗的时间域曲线的反射系数序列包括:
利用反射系数计算公式Ri=(Zi+1-Zi)/(Zi+1+Zi)将波阻抗的时间域曲线Zi换算得到时间域反射系数曲线Ri,从而获得反射系数序列, i表示采样点,i=1,2……,n。
5.根据权利要求4所述的定量计算地下反射系数振幅占比的方法,其特征在于,所述方法还包括:
设置阈值ε1,拾取大于ε1的Ri得到主要反射系数Rk,并按照局部极大和极小处理进行方波化,其中k=1,2……,m。
6.根据权利要求1所述的定量计算地下反射系数振幅占比的方法,其特征在于,所述选择指定的波峰或波谷,根据反射系数序列获取其反射系数构成,计算该波峰或波谷前后设定范围内的各反射系数在该波峰或波谷处的振幅响应值包括:
输入指定的波峰或波谷的序列号kk,拾取该波峰或波谷前后半个子波长度范围内的反射系数Rk’,其中k=1,2……N1,N1为该范围内的反射系数个数,分别对Rk’与输入的地震子波进行褶积,并计算在指定的波峰或波谷对应的时间处其振幅响应值,得到ZZi,ZZi为各反射系数在指定波峰或波谷处的振幅响应值;
其中i=1,2,……N1。
7.根据权利要求1所述的定量计算地下反射系数振幅占比的方法,其特征在于,所述振幅占比的计算公式为:
fsxszb(i)=ZZi/∑ZZi其中,fsxszb表示振幅占比,ZZi为各反射系数在指定波峰或波谷处的振幅响应值,i=1,2……N1。
8.根据权利要求1所述的定量计算地下反射系数振幅占比的方法,其特征在于,所述方法还包括:
输出各反射系数对指定的波峰或波谷的振幅占比图。
CN202210489190.2A 2022-05-07 2022-05-07 一种定量计算地下反射系数振幅占比的方法 Active CN114578452B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210489190.2A CN114578452B (zh) 2022-05-07 2022-05-07 一种定量计算地下反射系数振幅占比的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210489190.2A CN114578452B (zh) 2022-05-07 2022-05-07 一种定量计算地下反射系数振幅占比的方法

Publications (2)

Publication Number Publication Date
CN114578452A CN114578452A (zh) 2022-06-03
CN114578452B true CN114578452B (zh) 2023-01-06

Family

ID=81769006

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210489190.2A Active CN114578452B (zh) 2022-05-07 2022-05-07 一种定量计算地下反射系数振幅占比的方法

Country Status (1)

Country Link
CN (1) CN114578452B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106526670A (zh) * 2016-09-21 2017-03-22 中石化石油工程技术服务有限公司 一种碎屑岩储层中地震属性砂体空间分布描述及评价的方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7230541B2 (en) * 2003-11-19 2007-06-12 Baker Hughes Incorporated High speed communication for measurement while drilling
CN105005080B (zh) * 2015-06-16 2018-04-13 中国石油化工股份有限公司 一种利用振幅比率属性识别地层圈闭尖灭线的方法
CN105022097B (zh) * 2015-07-09 2018-12-18 广西大学 一种土质边坡滑动面综合预报方法
CN109782337B (zh) * 2018-12-28 2021-05-28 中国石油化工股份有限公司 一种基于地震属性图的有利储层边界自动拾取方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106526670A (zh) * 2016-09-21 2017-03-22 中石化石油工程技术服务有限公司 一种碎屑岩储层中地震属性砂体空间分布描述及评价的方法

Also Published As

Publication number Publication date
CN114578452A (zh) 2022-06-03

Similar Documents

Publication Publication Date Title
CN111239802B (zh) 基于地震反射波形和速度谱的深度学习速度建模方法
CN108802812B (zh) 一种井震融合的地层岩性反演方法
CN1170170C (zh) 用于地震分析的频谱解译
CN109188520B (zh) 薄储层厚度预测方法及装置
US7069150B2 (en) Method for optimizing migration fields using time slice analysis
WO2007087435A2 (en) Method for processing acoustic reflections in array data to image near-borehole geological structure
CN109470187B (zh) 基于地震三属性的储层厚度预测方法
CN107065013B (zh) 一种地震尺度下的层速度确定方法及装置
CN104237945A (zh) 一种地震资料自适应高分辨处理方法
CN105093306A (zh) 一种地球物理勘探中储层自动解释与厚度求取方法
CN103901465A (zh) 全息三维地震勘探观测系统的设计方法
CN107831542A (zh) Ddw高精度深度域井震匹配方法
CN104570076A (zh) 一种基于二分法的地震波初至自动拾取方法
CN108226997A (zh) 一种基于叠前地震数据的地震相划分方法
CN114609675A (zh) 基于高频旋回对碳酸盐岩地层沉积微地貌的定量恢复方法
CN114578452B (zh) 一种定量计算地下反射系数振幅占比的方法
CN109283577B (zh) 一种地震层位标定方法
CN110847887A (zh) 一种细粒沉积陆相页岩裂缝识别评价方法
CN113031070B (zh) 一种深度域合成地震记录的制作方法
CN113589365B (zh) 基于时频域信息的储层尖灭线描述方法
CN111624649B (zh) 一种零偏移距vsp建立横向变速层速度模型的方法和装置
CN110568490B (zh) 一种高速层顶薄储层的识别方法
CN110632660B (zh) 基于地震数据体的薄砂体表征方法及装置
US11927711B2 (en) Enhanced-resolution sonic data processing for formation body wave slowness with full offset waveform data
CN110389381B (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
TA01 Transfer of patent application right
TA01 Transfer of patent application right

Effective date of registration: 20221118

Address after: Room 1805, 18th Floor, Building 2, Yard 1, Shangdi 10th Street, Haidian District, Beijing 100000

Applicant after: DIDA HUINENG (BEIJING) TECHNOLOGY CO.,LTD.

Address before: 257000 No. 165, three West Road, Dongying District, Shandong, Dongying

Applicant before: Dongying Jingchuan Petroleum Technology Co.,Ltd.

GR01 Patent grant
GR01 Patent grant