CN116626760B - 基于自适应高阶最大熵wvd的地层不连续性检测方法及装置 - Google Patents

基于自适应高阶最大熵wvd的地层不连续性检测方法及装置 Download PDF

Info

Publication number
CN116626760B
CN116626760B CN202310649625.XA CN202310649625A CN116626760B CN 116626760 B CN116626760 B CN 116626760B CN 202310649625 A CN202310649625 A CN 202310649625A CN 116626760 B CN116626760 B CN 116626760B
Authority
CN
China
Prior art keywords
data
order
seismic
maximum entropy
wvd
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
CN202310649625.XA
Other languages
English (en)
Other versions
CN116626760A (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.)
Sichuan Natural Resources Investment Group Geophysical Exploration Institute Co ltd
Original Assignee
Sichuan Natural Resources Investment Group Geophysical Exploration Institute 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 Sichuan Natural Resources Investment Group Geophysical Exploration Institute Co ltd filed Critical Sichuan Natural Resources Investment Group Geophysical Exploration Institute Co ltd
Priority to CN202310649625.XA priority Critical patent/CN116626760B/zh
Publication of CN116626760A publication Critical patent/CN116626760A/zh
Application granted granted Critical
Publication of CN116626760B publication Critical patent/CN116626760B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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/301—Analysis for determining seismic cross-sections or geostructures
    • 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
    • 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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/10—Pre-processing; Data cleansing
    • 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/63—Seismic attributes, e.g. amplitude, polarity, instant phase
    • 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/64—Geostructures, e.g. in 3D data cubes
    • G01V2210/641—Continuity of geobodies
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Theoretical Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及地球物理勘探地震资料解译领域，具体为一种基于自适应高阶最大熵WVD的地层不连续性检测方法及装置。检测方法包括获取叠后地震数据；所述叠后地震数据进行滤波处理，得到滤波后地震数据；对所述滤波后地震数据做自适应高阶最大熵WVD变换，得到地震时频分布数据；由所述地震时频分布数据提取地震瞬时频率数据；对所述地震瞬时频率数据计算相干属性，得到相干剖面；识别所述相干剖面的异常区，输出地层不连续区。本方法利用从频率域计算的相干属性检测地层不连续性，具有更优的计算效率以及更可靠的检测结果。

Description

基于自适应高阶最大熵WVD的地层不连续性检测方法及装置
技术领域
本发明涉及地球物理勘探地震资料解译领域，具体为一种基于自适应高阶最大熵WVD的地层不连续性检测方法及装置。
背景技术
不连续性检测是地震勘探资料解译的重要环节。目前针对河道砂体、构造识别等不连续性检测技术主要包括：蚂蚁追踪、时频分析、相干属性及曲率属性。其中，相干技术和时频分析技术是不连续性检测的两种常见且有效的手段。
基于相干技术地层不连续性检测通常利用时间域地震信号从振幅出发计算相干属性。基于互相关原理的第一代相干算法（C1）计算量小，但必须要求地震资料具有极高的信噪比。基于多道相似技术的第二代相干算法（C2）在一定程度上摆脱了地震资料高信噪比的限制，但横向分辨率欠佳。基于本征值结构的第三代相干体计算（C3）提高了横向分辨率，但计算量较大，且对倾角变化不敏感。3种相干算法严格上并没有绝对的优劣，只是适用的条件不同，实际使用时通常根据使用条件和使用目的选取最为合适的相干算法。
时频分析方法作为不连续性检测的另一重要手段，主要包含线性、双线性（二次型）和非线性3种。线性方法包括短时傅里叶变换（STFT）、连续小波变换（CWT）、S变换（ST）和广义S变换（GST）等，优点在于计算效率高，缺陷在于时间分辨率和频率分辨率无法同时提高。双线性方法的代表为Wigner-Ville分布（WVD），其优点在于时频比较聚焦，缺点在于由于二次型分布的特性，不可避免的会产生虚假项，使得结果难以应用。非线性方法包括：经验模态分解（EMD）类和同步挤压变换（SST）类，优势在于时频分辨率高，不足之处是受信噪比影响较大。
综上，相干和时频分析作为不连续检测技术的两种重要手段，在各自的领域均可以实现对地层构造等不连续性的检测，但又都存在各自的局限。本申请发明人发现：对于地震信号，其波形的变化势必带来频率的变化，如何将相干技术与时频分析有机结合，充分发挥时频分析和相干两种技术的优势，是获得更加可靠的地层不连续性检测结果的关键。
发明内容
本发明的目的在于以高精度时频分析技术为基础，结合相干算法应用，从频率域进行地层不连续性检测，提高地层不连续性检测结果可靠性。
为了实现上述发明目的，本发明提供了以下技术方案：
基于自适应高阶最大熵WVD的地层不连续性检测方法，包括：
S1,获取叠后地震数据；
S2,对所述叠后地震数据进行滤波处理，得到滤波后地震数据；
S3,对所述滤波后地震数据做自适应高阶最大熵WVD变换，得到地震时频分布数据；
S4,由所述地震时频分布数据提取地震瞬时频率数据；
S5,对所述地震瞬时频率数据计算相干属性，得到相干剖面；
S6,识别所述相干剖面的异常区，输出地层不连续区。
地震剖面由多道地震数据组成，每道数据通常由道头文件和实际数据文件两部分组成。本方法既可以采用单核处理器逐道计算以实现，亦可以采用多核处理器并行计算或者多计算设备分布计算方式实现。
其中，S2对所述叠后地震数据进行滤波处理包括：
S21,采用短时傅里叶变换将所述叠后地震数据变换至时频域；
S22,在时频域对最小截止频率和最大截止频率之外的频点赋0值；
S23,对赋值后的时频域数据进行短时傅里叶反变换，得到滤波后单道地震数据。
S3对所述滤波后地震数据做自适应高阶最大熵WVD变换包括：
S31,以奇数长度窗口L对包含N个数据点的当前道地震数据进行端点延拓；
S32,对延拓后的数据进行希尔伯特变换，得到滤波后数据的复数道；
S33,以所述奇数长度窗口L对所述复数道信号进行逐点截取，得到N×L的矩阵数据；
S34,以N×L的矩阵数据的每一行与其共轭转置相乘，得到N×1的矩阵，作为求解最大熵自回归模型功率谱的预测误差初值E 0 ;
S35,对于N×L的矩阵数据计算初始阶数m=1的最大熵自回归模型功率谱的主自回归系数c m 和预测误差能量E m ，再令m=m+1，直至求得最高阶数mp下的每一阶的自回归系数c m 和预测误差能量E m ；
可知，E 0 是根据原始数据给定的，可等价于0阶，阶数m从1开始计算。
S36,利用莱文森递归方程计算剩余的自回归系数；
S37,根据自回归系数c与核函数的关系计算自适应高阶最大熵WVD的核函数，对所述核函数做傅里叶变换。
进一步地，S33中以所述奇数长度窗口L对所述复数道信号进行逐点截取具体为以N个数据的第一个数据点开始，以当前数据点为所述窗口L的中心截取其前(L-1)/2至后(L-1)/2的数据点作为当N×L矩阵的当前行，截取直至第N个数据点。
进一步地，采用伯格（Burg）前后向误差递推算法对于N×L的矩阵数据计算初始阶数m=1的最大熵自回归模型功率谱的主自回归系数c m 和预测误差能量E m 包括：
按行处理N×L的矩阵数据，设当前行的数据为z(n)，n∈[1,L]，设定初始前向预测误差efp和后向预测误差ebp:
计算主自回归系数c m ：
根据c m ，更新efp和ebp:
计算预测误差能量E m ：
式中，*表示共轭转置，下同。采用伯格算法计算的主自回归系数c m （等价于c i,j ，i= j=m）与后续莱文森递归方程计算剩余的自回归系数c m,k （等价于c i,j ，i≠j）共同构成的自回归系数c。
其中，所述最高阶数mp通过以下方式求得：
根据滤波后、逐点截取前的N×1数据计算初始阶数m=1的最大熵自回归模型功率谱的主自回归系数c m 和预测误差能量E m ；
令m=m+1,直至求得满足预设条件的最大阶数mp=m；
所述预设条件为预测误差能量E m 阈值T E 和阶数m阈值T m ，当Em<T E 且m≤T m 时，mp=m；当m≤T m 时不存在Em<T E ，则mp=T m 或经验值。经验值的选择综合考虑基本的数据处理要求（可略低于所述预设预测误差能量E m 阈值T E ）以及计算效率，例如若直至m=T m =10仍然不满足Em< T E ，则可指定mp=7。
本领域技术人员应当知晓，对N×1数据计算前向误差和后向误差与对N×L的矩阵数据的过程是一致的。区别仅在于N×L的矩阵数据初始长度是窗口截取后的L；而N×1数据的初始长度是N（单道地震数据长度）其E 0是该列数据的转置乘以该列数据，再除以N，得到所述E 0的值。
进一步地，S36中利用莱文森递归方程计算剩余的自回归系数包括：
。
本领域技术人员应当知晓：时间域信号x(n)做希尔伯特变换可得解析信号：
z(n)= x(n)+jH[x(n)]
式中，n∈[0, N-1]为采样点序号，N为最大采样点数。利用解析信号z(n)可生成N×N阶的自相关系数协方差矩阵C=zz T,上标T表示复共轭转置。矩阵C沿每条交叉对角线的序列构成了离散WVD的核函数K(n)，离散WVD核函数K(n)的累积是输入信号自相关函数的偶数项，表示为：
其中，，/>，k n (0)表示解析信号z(n)的自相关系数协方差矩阵C的主对角线元素；/>；定义k n (l)为：
且仅考虑l取大于等于0的值，离散WVD核函数K(n)的累积是输入信号自相关函数的偶数项，可改写为如下形式：
上式表明，对于单独的每个n进行取值，总有范围[0:2：N-1]的l与其对应，即：
若考虑l取负值，则R n (l)可以进一步改写为：
同时，结合伯格算法在莱文森递归方程中推导出的自回归系数c与核函数的关系：
可以推导出S37中自适应高阶最大熵WVD的核函数为：
式中，kn即为利用最大熵（AR）模型求取的WVD的核函数的右半部分。根据协方差矩阵的共轭性质，kn与WVD核函数的求取方法一致，只需求出右半部分，对其求共轭可得最大熵WVD核函数的左半部分；再对完整的最大熵WVD核函数做快速傅里叶变换（FFT），可求出WVD的最大熵瞬时功率谱，即WVD-MEM时频分布。
S4中利用所述自适应高阶最大熵WVD核函数求取瞬时频率f：
式中，f(n)表示瞬时频率，n为时间采样点，范围从1到N；q表示频率设置，范围从1到N；k n1表示最大熵核函数k n 的虚部；T表示采样周期；k n (0)等价于误差初值E 0，r表示校正系数。
S5采用基于多道相似原理的C2相干算法对所述地震瞬时频率数据计算相干属性。
S6识别所述相干剖面的异常区具体包括：遍历相干剖面筛选相干值小于预设相干值阈值的数据点作为地层不连续区的指示输出。
本申请还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述的基于自适应高阶最大熵WVD的地层不连续性检测方法的步骤。
本申请还提供一种计算机设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述的基于自适应高阶最大熵WVD的地层不连续性检测方法的步骤。
与现有技术相比，本发明的有益效果：
本发明从频率域计算地震数据相干属性，不同于利用时间域从振幅上计算相干属性，其计算的相干属性能够更好地反映地层的不连续性。地震数据的幅值变化较大，一般为从几百至几万的变化范围，难免会存在高频或低频的干扰，而地震数据的频率变化范围一般集中在5-80Hz之间，且频率的变化是与波形的变化紧密联系的，从频率域计算的相干属性具有更可靠的地层不连续性检测结果。
本发明采用自适应高阶最大熵WVD的时频分析方法，利用最大熵模型的特性消除了WVD存在的干扰项，具有高度的时频聚焦性和分辨率。
进一步地，本发明并未使用常规的一阶中心矩来提取瞬时频率，而是利用了扩展核函数来提取瞬时频率，跳过了傅里叶变换这一环节，提取的瞬时频率相比一阶中心矩，噪声更少，波形更光滑，计算效率更高。
此外，本发明针对单道数据利用预测误差自动定阶，使得每一道数据均可按照数据本身特性自动确定最合理阶数，避免了对整个数据体赋值相同的阶数，因而时频分析结果也更加合理可靠，也为提取更加合理可靠的瞬时频率奠定了基础。
附图说明
图1为实施例1地层不连续性检测方法流程图；
图2为实施例1滤波前后地震信号对比图（单道）；
图3为实施例1自适应高阶最大熵WVD与其他时频分析方法结果对比图；
图4为实施例1基于自适应高阶最大熵WVD提取的瞬时频率与其他瞬时频率提取方法结果对比图；
图5为实施例1基于自适应高阶最大熵WVD提取的瞬时频率相干剖面与基于地震振幅的相干剖面结果对比图一；
图6为实施例1基于自适应高阶最大熵WVD提取的瞬时频率相干剖面与基于地震振幅的相干剖面结果对比图二。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图，对本发明实施例中的技术方案进行清楚、完整的描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。
因此，以下对本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的部分实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征和技术方案可以相互组合。
实施例1提供一种基于自适应高阶最大熵WVD的地层不连续性检测方法，如图1所示包括：
S1,获取叠后地震数据；
S2,对所述叠后地震数据进行滤波处理，得到滤波后地震数据；
S3,对所述滤波后地震数据做自适应高阶最大熵WVD变换，得到地震时频分布数据；
S4,由所述地震时频分布数据提取地震瞬时频率数据；
S5,对所述地震瞬时频率数据计算相干属性，得到相干剖面；
S6,识别所述相干剖面的异常区，输出地层不连续区。
以下详述具体各个步骤。
获取叠后地震数据
叠后地震数据是石油工业中最常用的勘探数据，通常由道头文件和实际数据文件两部分组成。道头文件包含线道号、采样率、采样时间范围等信息；实际数据文件为一个矩阵，横轴代表道号，纵轴代表时间，矩阵的值为时间域的振幅属性。
滤波处理
本实施例采用短时傅里叶变换（STFT）及其反变换技术进行时频域滤波：
1）生成短时傅里叶变换窗口，设置窗口类型为高斯窗，设置窗口长度；
2）设定滤波最小截止频率和最大截止频率；
3）对地震数据单道按设定的窗口做短时傅里叶变换，获得时频域矩阵；
4）对短时傅里叶变换后的矩阵，在时频域对每个时间采样点小于最小截止频率和大于最大截止频率的成分赋0值，消除高频和低频的干扰；
5）对赋值后的矩阵进行短时傅里叶反变换，获得滤波后的单道地震数据；
6）对剩余地震道数据按照上面3）~5）进行逐道处理，获得滤波后的地震数据剖面（矩阵）。
短时傅里叶变换（STFT）及其反变换具有计算效率高的特点，可以较好地消除背景干扰噪声。图2展示了原始单道地震信号和利用短时傅里叶变换及其反变换滤波后的地震信号。可见，原始信号噪声干扰非常严重，波形较毛糙、不清晰；经过滤波后的数据，波形光滑且清晰，证明已消除了低频和高频噪声的干扰，为后续可靠的时频分析提供了基础。
自适应高阶最大熵WVD变换
1）设定奇数窗口长度L，选择滤波后数据的第1道，利用自回归模型（Autoregressive model，简称AR模型）对该道数据端点延拓，以消除端点效应；
2）对延拓后的数据进行希尔伯特变换，获得滤波后数据的复数道（复数信号，也叫做解析信号），复数道包含实部和虚部两部分，实部信号与原始信号一致；
3）以奇数窗口L对复数道信号进行逐点截取，获得N×L规模的矩阵数据；
对逐点截取的具体过程进行说明：假设原始数据为1/2/3/4/5/6/7/8/9共9个数字，那么N为9。假设L为3，且进行了延拓，假设延拓后的数据为0/0/0/1/2/3/4/5/6/7/8/9/0/0/0，相比原始数据，两端各多了3个数据。逐点截取后，最终获得N×L（9×3）矩阵中，第一行则为0/1/2，第2行则为1/2/3，依次类推，共9行，即从N个数据的第一个数据点开始，以当前数据点为所述窗口L的中心截取其前(L-1)/2至后(L-1)/2的数据点作为当N×L矩阵的当前行，截取直至第N个数据点；
4）计算预测误差初值E 0：
预测误差初值E 0为一个规模为N×1的矩阵（或者叫一列数），其根据上一步骤3）中的N×L矩阵计算，对于N×L矩阵的每一行，均利用该行数据乘以它的共轭转置，计算出一个数据结果，将每一行的数据结果组合，即得规模为N×1的矩阵（预测误差初值E 0）；
5）对于阶数m，给定初始阶数为1。利用伯格算法和莱文森递归方程中的预测误差能量计算公式，计算自回归系数c m 和预测误差能量E m ；
伯格算法步骤如下：
①对于3）中N×L的矩阵数据，按行处理，首先对于第一行数据，假设其为z(n)，n的取值范围为[1,L]。给定前后向预测误差分别为efp和ebp:
②计算主自回归系数c m ：
③根据c m ，更新efp和ebp:
④计算预测误差能量E m ：
6）根据滤波后、逐点截取前的N×1的一列数据，设定m=1,利用步骤1）~4），重新计算整个序列的预测误差能量E m （m=1），计算完成后，令m = m+1,再次利用步骤4）计算高阶总预测误差能量E m 。设定最大阶数为10，自适应判断合适的阶数mp。当E m <0.01，且m<=10，mp=m；当m>10，令mp=7。
需要说明的是，最大阶数mp的求取相对于N×L矩阵数据的处理的顺序不是固定的，即既可以置于N×L矩阵数据开始求取自回归系数c m 和预测误差能量E m 之后，也可以先行根据N×1滤波后列数据求取最大阶数mp之后再对N×L矩阵数据求取主自回归系数c m 和预测误差能量E m 。
以整道地震数据的预测误差能量为基准自动确定每一道地震数据的合适阶数，不需要再对所有道的地震数据赋予相同的阶数，因而变换结果不仅保留了最大熵WVD的高分辨特性，同时更符合实际信号分解规律。其本质上属于数据驱动的分解方法，从1阶开始，随着阶数上升，其预测误差能量下降较快的，其最终阶数越小，反之其预测误差能量下降较慢的，其最终阶数较大。由信号本身根据自身的预测误差能量来确定阶数，从而最终确定时频分布，其时频分布结果也将更为合理。
7）当利用整个序列的预测误差能量E m 确定好最佳mp后，对于N×L的矩阵数据，令m= m+1,按步骤1）~4）计算每一阶的主自回归系数c m 和预测误差能量E m ，一直到m = mp；
8）当主自回归系数c m 求出后，按照莱文森递归方程，可求出剩余的回归系数；
莱文森递归方程如下：
9）自回归系数c和核函数的关系如下：
因而，可以推导自适应高阶最大熵WVD的核函数为：
上式中，k n 即为最终求出的自适应高阶最大熵WVD的核函数，对其做傅里叶变换（或者快速傅里叶变换），即可得到自适应高阶最大熵WVD时频分布；
10）当单道数据求出后，对滤波后地震数据逐道按步骤1）~9）求取，即可得到所有地震道的时频分布。
图3所示，从左到右依次为：滤波后的单道实测地震信号、广义S变换（GST）、Wigner-Ville分布（WVD）、平滑伪Wigner-Ville分布（SPWVD）、最大熵WVD以及自适应高阶最大熵WVD。如蓝色虚线框内所示，GST由于其为线性时频分析方法，时频较模糊，能量不聚焦，分辨率较低；WVD时频谱比较聚焦，但交叉项干扰严重，并不适用地震数据谱分解；SPWVD相比WVD消除了部分交叉项，消除不彻底，且降低了分辨率；最大熵WVD的时频聚焦性和时频分辨率都很高；自适应高阶最大熵WVD相比最大熵WVD的优势在于，不仅聚焦性和分辨率更高，且不需要设定阶数，其时频结果中每个时间均对应多个频率成分，尤其对复杂的地震信号而言，更符合地震波形规律。
提取地震瞬时频率
无论是峰值检测法还是一阶中心矩法，求瞬时频率的原理均利用了时频分布这个概念，这也带来了时间上的额外消耗。
本申请直接利用自适应高阶最大熵WVD核函数通过积分求取瞬时频率（扩展核函数法），绕过傅里叶变换这一环节，从而能极大地节省时间和提高计算效率。
算法公式如下：
式中，f(n)表示瞬时频率，n为时间采样点，范围从1到N；q表示频率设置，范围从1到N；k n1表示最大熵核函数k n 的虚部；T表示采样周期；k n (0)等价于误差初值E 0，r表示校正系数。
图4从上往下分别为：单道实测地震信号以及峰值检测法、一阶中心矩法、扩展核函数法提取的瞬时频率。以蓝色虚线框内所示，3种方法提取结果基本一致。其中峰值检测法和扩展核函数法计算的瞬时频率更为更滑，一阶中心矩法提取的瞬时频率明显抗噪性稍弱。
表1 3种方法提取单道实测地震信号瞬时频率耗时
表1为利用峰值检测法、一阶中心矩法以及扩展核函数法3种方法提取的单道实测地震信号耗时表，对每种方法分别测试5次，取平均值。可见，峰值检测法平均用时0.0732852s，一阶中心矩法平均用时0.1934980s，扩展核函数法平均用时0.0249848s。扩展核函数法相比峰值检测法计算效率提高了293.32%，相比一阶中心矩法计算效率提高了774.46%。因此，在提取结果基本一致的情况下，计算效率成了关键。扩展核函数法仅需要核函数的参与，跳过了计算傅里叶变换这一环节，计算效率提升极大。
计算相干属性
本实施例采用C2算法计算相干属性。由于计算的是利用时频分布提取的瞬时频率而不是原始的振幅值，对于地震数据，其波形的变化必然导致频率的变化，利用C2相干算法计算频率的相干属性，其检测结果可以直接反映地震波形的变化，且由于避过了时域波形，而利用了频域数据，其计算出的相干属性也更加稳定可靠。
此外，C2算法计算效率高，其输入量为瞬时频率，瞬时频率对应了地震波形，且瞬时频率的相似程度（相干性）反映了地震波形（地震数据）的相似度（相干性），又因为地震波形的连续性与地层的连续性是一致的，因此，瞬时频率的相似程度（相干性）能够反映了地层的连续性，能够很好地用于地层不连续性检测。
检测地层不连续性
利用相干剖面，可以检测地层界面的不连续性。相干技术的基本原理是对地震数据空间连续性的表征，并对异常体突出共性，将原振幅数据中多线、多道、多点的信息浓缩到一个点上来反映该点地质特征，因此，相干技术对不连续性十分敏感。
本实施例中，根据相干剖面识别地层不连续性既可以通过地震资料解译人员读图形式进行。本领域人员可以理解，为进一步提高自动化只能辅助地震资料解译过程，还可以通过算法遍历相干剖面筛选相干值小于预设相干值阈值的数据点作为地层不连续区的指示输出。
图5和图6为基于自适应高阶最大熵WVD提取的瞬时频率和C2相干算法在地层不连续性检测的实际应用效果图。4幅子图从左往右、从上往下依次是：原始数据滤波后的剖面、频率剖面、振幅相干剖面、基于自适应高阶最大熵WVD和C2算法的频率相干剖面。其中，白色表示低相干值，黑色表示高相干值。相干属性对应地震数据的连续性，地震数据的连续性又与地质体的连续性所对应。当地层连续时，地震道之间具有很高的相似性，表现为高相干值；当地层不连续时，地震道之间相似性不高，表现为低相干值。因此，相干体与地层对应，可直接反映地层的不连续性。
可见，利用本发明计算的相干剖面，由于利用了自适应高阶最大熵WVD的自动确定阶数和高分辨的优势，同时在频率域求取相干，结果更稳定，剔除了部分噪声干扰，纵向上对地层不连续性的分辨能力更强（红色矩形框内所示）。
综上所述，本发明采用自适应高阶最大熵WVD的时频分析方法，利用最大熵模型的特性，消除了WVD存在的干扰项，具有高度的时频聚焦性和分辨率，同时，针对各单道利用预测误差来自动确定最合理阶数，时频分析结果也更加合理可靠，使最终提取的瞬时频率也更为合理可靠；利用了扩展核函数来提取瞬时频率，跳过了傅里叶变换这一环节，提取的瞬时频率相比一阶时间矩，噪声更少，波形更光滑，计算效率更高；从频率域计算地震数据相干属性，能够更好地反应地层的不连续性。
以上实施例仅用以说明本发明而并非限制本发明所描述的技术方案，尽管本说明书参照上述的各个实施例对本发明已进行了详细的说明，但本发明不局限于上述具体实施方式，任何对本发明进行修改或等同替换的而不脱离发明的精神和范围的技术方案及其改进，其均涵盖在本发明的权利要求范围当中。

Claims (9)

1.基于自适应高阶最大熵WVD的地层不连续性检测方法，其特征在于，包括：
S1,获取叠后地震数据；
S2,对所述叠后地震数据进行滤波处理，得到滤波后地震数据；
S3,对所述滤波后地震数据做自适应高阶最大熵WVD变换，得到地震时频分布数据；
S4,由所述地震时频分布数据提取地震瞬时频率数据；
S5,对所述地震瞬时频率数据计算相干属性，得到相干剖面；
S6,识别所述相干剖面的异常区，输出地层不连续区；其中，
S3对所述滤波后地震数据做自适应高阶最大熵WVD变换包括：
S31,以奇数长度窗口L对包含N个数据点的当前道地震数据进行端点延拓；
S32,对延拓后的数据进行希尔伯特变换，得到滤波后数据的复数道；
S33,以所述奇数长度窗口L对所述复数道信号进行逐点截取，得到N×L的矩阵数据；
S34,以N×L的矩阵数据的每一行与其共轭转置相乘，得到N×1的矩阵，作为求解最大熵自回归模型功率谱的预测误差初值E0；
S35,对于N×L的矩阵数据计算初始阶数m＝1的最大熵自回归模型功率谱的主自回归系数cm和预测误差能量Em，再令m＝m+1，直至求得最高阶数mp下的每一阶的自回归系数cm和预测误差能量Em；
S36,利用莱文森递归方程计算剩余的自回归系数；
S37,根据自回归系数c与核函数的关系计算自适应高阶最大熵WVD的核函数，对所述核函数做傅里叶变换。
2.根据权利要求1所述的方法，其特征在于，所述最高阶数mp通过以下方式求得：
根据滤波后、逐点截取前的N×1数据计算初始阶数m＝1的最大熵自回归模型功率谱的主自回归系数cm和预测误差能量Em；
令m＝m+1,直至求得满足预设条件的最大阶数mp＝m；
所述预设条件为预测误差能量Em阈值TE和阶数m阈值Tm，当Em<TE且m≤Tm时，mp＝m；当m≤Tm时不存在Em<TE，则mp＝Tm或经验值。
3.根据权利要求1所述的方法，其特征在于，对于N×L的矩阵数据计算初始阶数m＝1的最大熵自回归模型功率谱的主自回归系数cm和预测误差能量Em包括：
按行处理N×L的矩阵数据，设当前行的数据为z(n')，n'∈[1,L]，设定初始前向预测误差efp和后向预测误差ebp:
计算主自回归系数cm：
根据cm，更新efp和ebp:
计算预测误差能量Em：
4.根据权利要求1所述的方法，其特征在于，S36中利用莱文森递归方程计算剩余的自回归系数包括：
5.根据权利要求1所述的方法，其特征在于，S37中自适应高阶最大熵WVD的核函数为：
6.根据权利要求5所述的方法，其特征在于，S4中利用所述自适应高阶最大熵WVD核函数求取瞬时频率f：
式中，f(n)表示瞬时频率，n为时间采样点，范围从1到N；q表示频率设置，范围从1到N；kn1表示最大熵核函数kn的虚部；T表示采样周期；kn(0)等价于误差初值E0，r表示校正系数。
7.根据权利要求1所述的方法，其特征在于，S5采用基于多道相似原理的C2相干算法对所述地震瞬时频率数据计算相干属性；
S6识别所述相干剖面的异常区具体包括：遍历相干剖面筛选相干值小于预设相干值阈值的数据点作为地层不连续区的指示输出。
8.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1-7任一项所述方法的步骤。
9.一种计算机设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1-7任一项所述方法的步骤。
CN202310649625.XA 2023-06-02 2023-06-02 基于自适应高阶最大熵wvd的地层不连续性检测方法及装置 Active CN116626760B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310649625.XA CN116626760B (zh) 2023-06-02 2023-06-02 基于自适应高阶最大熵wvd的地层不连续性检测方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310649625.XA CN116626760B (zh) 2023-06-02 2023-06-02 基于自适应高阶最大熵wvd的地层不连续性检测方法及装置

Publications (2)

Publication Number Publication Date
CN116626760A CN116626760A (zh) 2023-08-22
CN116626760B true CN116626760B (zh) 2024-05-10

Family

ID=87638017

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310649625.XA Active CN116626760B (zh) 2023-06-02 2023-06-02 基于自适应高阶最大熵wvd的地层不连续性检测方法及装置

Country Status (1)

Country Link
CN (1) CN116626760B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009011735A1 (en) * 2007-07-16 2009-01-22 Exxonmobil Upstream Research Company Geologic features from curvelet based seismic attributes
CN103235339A (zh) * 2013-04-09 2013-08-07 中国石油大学(北京) 一种时频分解地震流体识别方法
CN104316958A (zh) * 2014-10-20 2015-01-28 中国石油天然气集团公司 一种识别不同尺度地层断裂的相干处理方法
CN109375265A (zh) * 2018-08-22 2019-02-22 中国地质大学（武汉） 一种基于变相位雷克子波匹配追踪的理想地震谱分解方法
CN111708084A (zh) * 2020-04-13 2020-09-25 中铁二院工程集团有限责任公司 一种基于地震映像数据的高精度岩溶识别方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7649805B2 (en) * 2007-09-12 2010-01-19 Schlumberger Technology Corporation Dispersion extraction for acoustic data using time frequency analysis

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009011735A1 (en) * 2007-07-16 2009-01-22 Exxonmobil Upstream Research Company Geologic features from curvelet based seismic attributes
CN103235339A (zh) * 2013-04-09 2013-08-07 中国石油大学(北京) 一种时频分解地震流体识别方法
CN104316958A (zh) * 2014-10-20 2015-01-28 中国石油天然气集团公司 一种识别不同尺度地层断裂的相干处理方法
CN109375265A (zh) * 2018-08-22 2019-02-22 中国地质大学（武汉） 一种基于变相位雷克子波匹配追踪的理想地震谱分解方法
CN111708084A (zh) * 2020-04-13 2020-09-25 中铁二院工程集团有限责任公司 一种基于地震映像数据的高精度岩溶识别方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于最大熵准则的Wigne...与微型古河道刻画方法及应用;徐天吉 等;石油勘探与开发;20211231;第48卷(第6期);1175-1186 *

Also Published As

Publication number Publication date
CN116626760A (zh) 2023-08-22

Similar Documents

Publication Publication Date Title
CN110361778B (zh) 一种基于生成对抗网络的地震数据重建方法
CN111723329B (zh) 一种基于全卷积神经网络的震相特征识别波形反演方法
CN111596366B (zh) 一种基于地震信号优化处理的波阻抗反演方法
CN113093272A (zh) 基于卷积编码的时间域全波形反演方法
CN111766628A (zh) 一种预条件的时间域弹性介质多参数全波形反演方法
CN113687433A (zh) 一种基于Bi-LSTM的大地电磁信号去噪方法及系统
CN113887398A (zh) 一种基于变分模态分解和奇异谱分析的gpr信号去噪方法
CN113777650A (zh) 一种基于混合范数和小波变换的稀疏时频谱分解方法、装置、设备及存储介质
CN113608259A (zh) 一种基于iceemdan约束广义s变换的地震薄层检测方法
Kimiaefar et al. Random noise attenuation by Wiener-ANFIS filtering
CN112183407B (zh) 一种基于时频域谱减法的隧道地震波数据去噪方法及系统
CN102289715A (zh) 基于前向线性预测的自适应小波神经网络去噪建模方法
CN112526599A (zh) 基于加权l1范数稀疏准则的子波相位估计方法及系统
CN116626760B (zh) 基于自适应高阶最大熵wvd的地层不连续性检测方法及装置
CN113655522A (zh) 频率域地震弱信号的增强方法
Jiang et al. Seismic wavefield information extraction method based on adaptive local singular value decomposition
CN109782346B (zh) 一种基于形态成分分析的采集脚印压制方法
CN112817056B (zh) 一种大地电磁信号去噪方法及系统
CN111856559B (zh) 基于稀疏贝叶斯学习理论的多道地震谱反演方法及系统
CN111239805B (zh) 基于反射率法的块约束时移地震差异反演方法及系统
CN112363217A (zh) 一种地震数据随机噪声压制方法及系统
Nose-Filho et al. Algorithms for sparse multichannel blind deconvolution
CN107315713B (zh) 一种基于非局部相似性的一维信号去噪增强方法
Ursin Seismic signal detection and parameter estimation
CN110687605A (zh) 基于改进的k-svd算法在地震信号处理的算法分析应用

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
CB02 Change of applicant information
CB02 Change of applicant information

Address after: No.39 Qinghua Road, Qingyang District, Chengdu, Sichuan 610072

Applicant after: Sichuan Natural Resources Investment Group Geophysical Exploration Institute Co.,Ltd.

Address before: No.39 Qinghua Road, Qingyang District, Chengdu, Sichuan 610072

Applicant before: SICHUAN ZHONGCHENG COALFIELD GEOPHYSICAL ENGINEERING INSTITUTE CO.,LTD.

GR01 Patent grant
GR01 Patent grant