CN107015276A - 一种基于改进霍夫变换的自适应时频峰值滤波消噪方法 - Google Patents

Abstract

本发明涉及一种基于改进霍夫变换的自适应时频峰值滤波消噪方法，属于地球物理技术领域。针对直线型同相轴地震记录提出一种多层分数阶傅里叶霍夫，对地震数据进行同相轴包络识别，建立与同相轴走势相近的滤波轨线方程，实现时频峰值滤波最优轨线的自适应选取，从而充分降低滤波偏差。本发明能够在强噪声条件下快速得出检测结果，具有较高的准确性及较强的抗干扰能力。

Description

一种基于改进霍夫变换的自适应时频峰值滤波消噪方法

技术领域

本发明属于地球物理技术领域，尤其是指针对包含直线型同相轴的地震记录提出了一种基于多层分数阶傅里叶霍夫变换的自适应时频峰值滤波消噪方法，该方法通过一种新型的霍夫变换获取地震数据中的同相轴包络斜率，进而建立具有该斜率的滤波轨线方程并沿着这些轨线对含噪地震数据进行样本重采样，其目的是使得采样后数据中有效信号线性度充分降低，最后对重采样后数据进行时频峰值滤波有效压制噪声，为地震勘探信号消噪处理及地震数据中有效信号识别提供了新方法。

背景技术

地震勘探作为一种重要的物理勘探手段，在油田的开发中起着重要的作用。野外地震数据采集和地震资料的数字处理技术是地震勘探的两个重要环节。只有得到高质量的地震剖面图，才能为石油勘探提供可靠的依据。在利用检波器进行实际的野外数据采集过程中，有效地震信号经常会受到各种噪声的干扰(规则噪声和随机噪声)，这些都直接影响地震解释工作的准确性。相比之下，随机噪声是不规则的、没有规律的，且相邻道间彼此是互不相关的，它们没有固定的频率，几乎分布于整个频带，严重影响着地震记录的信噪比。因此随机噪声压制是地震数据处理中的重点和难点。

近几年来，对于随机噪声压制，除了目前发展较为成熟的几种消噪算法，如自适应滤波、经验模态分解、多项式拟合、f-x预测滤波、小波变换、维纳滤波等，时频峰值滤波算法被成功地应用于陆地地震勘探随机噪声压制，并先后发展出多种改进算法，如局部径向轨线时频峰值滤波、二次轨线时频峰值滤波及分段轨线时频峰值滤波等。时频峰值滤波方法较上述成熟算法在强随机噪声、非平稳信号及先验信息未知等条件下存在明显的优势，但现有时频峰值滤波方法中轨线参数多针对某一地震记录人为设置，灵活性较差。

发明内容

本发明提供一种基于改进霍夫变换的自适应时频峰值滤波消噪方法，以解决现有轨线型时频峰值滤波方法中存在的轨线参数及形状固定，灵活性和适应性较差问题。目的在于针对包含直线型同相轴的地震记录，提出一种改进的霍夫变换，并将其应用与地震记录中的同相轴形态识别，从而实现时频峰值滤波方法中最优轨线的自适应选取。

本发明采取的技术方案是，包括下列步骤：

步骤一：一张含噪的、包含直线型同相轴的地震记录，记录中的反射同相轴就相当于图像中的边缘信息，为了便于对地震数据进行形态识别，首先对地震数据进行灰度处理，将数据值转换到0-255之间，处理后的地震数据为一张图像数据，接着，利用改进的霍夫变换实现地震记录中同相轴形态识别，具体过程为：先利用公式对灰度处理后的图像数据进行多层分数阶傅里叶变换，将地震数据变换到时域和频域之间的一个中间域，将其旋转到从时间轴指向频率轴的一个角度，得到多组不同阶数的傅里叶变换结果；再将变换结果转换成由模值ρ和相角θ构成的极坐标形式，在极坐标内，对数据进行共轭对称，并沿模值方向进行一维傅里叶反变换便可得到一张由多条正弦线构成的正弦谱图；实际上，从多层分数阶傅里叶变换到转换成极坐标，再到进行一维傅里叶反变换，整个过程相当于对地震数据进行了一次霍夫变换，得到的正弦谱图中多条正弦线的交点，即正弦谱矩阵中的数据最大值就对应了地震数据中的同相轴；对于直线型同相轴的地震记录，通过搜索正弦谱图中的峰值便可得到与其对应的地震数据中的某条同相轴在极坐标内的模值ρ_peak和相角θ_peak，再根据极坐标与原图像坐标之间的关系式ρ＝xcosθ+ysinθ便可以得到地震数据内同相轴包络线的斜率k_l＝ctanθ_peak；于是沿着多条斜率为k_l的直线型滤波轨线对含噪地震记录进行样本重采样，有效提高采样后信号的线性度，降低滤波偏差；

步骤二：根据检测到的同相轴包络斜率k_l建立与同相轴走势接近的滤波轨线方程y_l＝k_lx_l+b_l，其中x_l为地震数据点的横坐标，代表炮检距方向，y_l为地震数据点的纵坐标，代表时间方向，k_l为轨线斜率，b_l为截距，决定轨线与地震数据纵坐标的交点位置，改变截距b_l的数值便会生成多条具有相同斜率的直线轨线，沿着这些轨线对地震数据进行样本重采样，便可实现对地震数据沿同相轴形态方向的数据提取；

步骤三：对采样后的数据样本，利用时频峰值滤波方法进行噪声压制，由于对地震数据进行样本重采样时所采用的轨线斜率是通过检测同相轴形态得出的，因此轨线与同相轴趋势接近，沿着这些轨线重采样后的地震数据有效信号线性度充分提高，主频充分降低，因此对其进行时频峰值滤波时所采用的滤波窗长比传统时频峰值滤波方法中长很多，较长的窗长在完整保持有效成分的同时可以有效地压制随机噪声。

本发明所述步骤一中对于直线型同相轴，地震数据中的某一点(x,y)与极坐标内的参数对(ρ，θ)满足如下关系：ρ＝xcosθ+ysinθ,θ∈[0,π]，其中θ为直线法线与x轴的夹角，ρ为图像空间中原点到直线的距离。

本发明以提高时频峰值滤波方法的适用性和灵活性，有效压制地震记录中随机噪声为目的，针对包含直线型同相轴地震记录，基于多层分数阶傅里叶提出了一种改进的霍夫变换，并将其应用于同相轴形态检测，实现了时频峰值滤波方法中最优滤波轨线的自适应选取。改进的霍夫变换与现有霍夫变换相比具有较高的分辨率，没有梯度信息取向的损失，忽略了从笛卡尔到极坐标映射的零填充，在同相轴形态检测方面计算速度快，准确性高。

多层分数阶傅里叶变换在频域内具有较高的分辨率，没有梯度信息取向的损失，忽略了从笛卡尔到极坐标映射的零填充，计算速度明显提高。而现有的随机霍夫变换有时不能检测到短线段；线段检测器对噪声比较敏感，在嘈杂的图像中检测性能明显下降。

在重采样过程中由于位移-时间域内的样本格点与新作用域内的样本格点不一致，部分样本点不能落到具体的地震道上。如果丢弃这些“脱道”的点，采样得到的序列中就会出现一些断点和尖峰脉冲。本发明通过在采样前在相邻地震道间插入一些虚拟地震道的方式对地震数据进行扩充，使得原本“脱道”的点能够落在具体的地震道上，从而有效降低采样误差。

本发明的有益效果是提出了一种基于多层分数阶傅里叶变换实现的改进霍夫变换方法，并将其应用于直线型地震记录中同相轴形态检测，获取同相轴斜率。改进的霍夫变换与现有霍夫变换相比具有较高的分辨率，没有梯度信息取向的损失，忽略了从笛卡尔到极坐标映射的零填充，在同相轴形态检测方面计算速度快，准确性高。同相轴形态的准确识别轨线型时频峰值滤波方法中最优滤波轨线的自适应选取为提供了先决条件。此外，沿滤波轨线进行的插值重采样方案通过在相邻地震道间插入虚拟地震道的方式，有效避免了采样点的脱道现象，降低了重采样过程造成的误差。以上方面均对充分降低滤波偏差，完整保持有效信号具有重要意义。

附图说明

图1(a)图像空间内直线各参数情况，其中点P₁、P₂和P₃为直线上的三个点，直线与坐标原点距离为ρ，直线与y轴夹角为θ；

图1(b)霍夫变换参数空间内正弦曲线谱图，其中三条正弦曲线与点P₁、P₂和P₃三点对应，三条线相交于一点，该点的模值和相角即对应了图1a中直线与坐标原点的距离ρ及其与y轴夹角θ；

图2(a)具有不同噪声强度的含噪地震记录，该记录中包含三条直线型同相轴，其中区域R₁信噪比较低，有效信号被淹没其中难以辨识；区域R₂包含较强的脉冲噪声；区域E₁存在交叉轴；

图2(b)变换后得到的参数空间的正弦谱图，其中存在三个明显的交点(黑色方框标出)，分别对应原地震数据中的三条直线型同相轴；

图2(c)变换后得到的参数空间的正弦谱图的三维图，其中的三个尖峰就是图2b正弦谱图中的三个交点，通过搜索峰值便可找到正弦谱图中的交点位置；

图2(d)含噪记录及同相轴检测结果，图中黑色线给出了三条同相轴的检测结果，可以看到检测结果与同相轴几乎完全吻合；

图3重采样过程中虚拟地震道插值示意图，其中x方向为地震道方向，x_i-1，x_i，x_i+1，x_i+2为实际的地震道，其间虚线为插入的虚拟地震道；t方向为时间方向，点P₁、P₂和P₃为与时刻t_i-1，t_i，t_i+1对应的某一滤波轨线上的三个点，可见点P₂没有落到具体的地震道上，而点P₄为对点P₂横坐标四舍五入对应的点，插入虚拟地震道后，点P₂可以落到某一具体的虚拟地震道上，有效提高了数据采样的准确性；

图4(a)实际的含噪地震记录，中包含的同相轴呈现直线形态；

图4(b)传统时频峰值滤波消噪结果，其滤波方向沿着时间方向，选取的窗长为9，滤波后同相轴可见；

图4(c)基于形态检测的自适应时频峰值滤波消噪结果，其滤波方向沿着检测到的同相轴方向，滤波窗长为32，滤波后随机噪声得到了有效压制，同相轴清晰连贯。

具体实施方式

包括下列步骤：

步骤一：一张含噪的地震记录可以看成是一张图像数据，记录中的反射同相轴就相当于图像中的边缘信息。为了便于对地震数据进行形态识别，首先对地震数据进行灰度处理，将数据值转换到0-255之间，处理后的地震数据可看成是一张图像数据；接着利用改进霍夫变换对灰度化的图像数据中的同相轴形态进行识别，首先，利用式(1)中的分数阶傅里叶变换公式将地震数据变换到时域和频域之间的一个中间域，将其旋转到从时间轴指向频率轴的一个角度；

其中F^α(k₁,k₂)的N×N个频率分散在[-απ,απ]×[-απ,απ]区间内，给定一个分数系数α，就会产生一系列的分数阶傅里叶频率，于是有

P_i{(α_ik₁,α_ik₂)|-N/2≤k₁,k₂≤N/2}， (2)

其中i＝1,2,…,L，L是层的个数，整个频域内的数据点为L层数据点的合并，

L越大，得到的频率分辨率越高，对含噪地震数据图像进行多层分数阶傅里叶变换后，得到多个不同阶数的傅里叶变换频域散点图；

接着，将变换得到的散点图转换成由模值ρ和相角θ构成的极坐标形式，在极坐标内，对数据进行共轭对称，并沿模值方向进行一维傅里叶反变换便可得到一张由多条正弦线构成的正弦谱图；实际上，从多层分数阶傅里叶变换到转换成极坐标，再到进行一维傅里叶反变换，整个过程相当于对地震数据进行了一次霍夫变换过程；

得到的正弦谱图中多条正弦线的交点(即正弦谱矩阵中的数据最大值)就对应了地震数据中的同相轴；于是，对于直线型同相轴，通过搜索正弦谱图中峰值便可得到与其对应的地震数据中的某一线条在极坐标内的模值ρ_peak和相角θ_peak，再根据极坐标与原图像坐标之间的关系式ρ＝xcosθ+ysinθ便可以得到地震数据内同相轴包络线的斜率k_l＝ctgθ_peak，于是就可以沿着多条斜率为k_l的直线对含噪地震记录进行样本重采样，使得采样轨迹充分逼近同相轴形态，从而有效降低采样后得到信号的线性度，降低滤波偏差；

对于直线型同相轴，地震数据中的某一点(x,y)与极坐标内的参数对(ρ，θ)满足如下关系：

ρ＝xcosθ+ysinθ,θ∈[0,π] (4)

其中θ为直线法线与x轴的夹角，ρ为图像空间中原点到直线的距离，图1给出了图像空间与参数空间的对应关系示意图。图1(a)为图像数据空间中某一直线上的点P₁、P₂和P₃，图1(b)为这三个点对应的极坐标空间的三条正弦曲线，这三条正弦线交于坐标为(ρ₁,θ₁)的点P，通过点P的模值和相角信息便可得到原图像数据中的直线(即同相轴)的斜率；

为了验证改进的多层分数阶傅里叶霍夫变换方法在地震记录同相轴形态检测中的有效性，图2(a)给出了一张40道，每道1024个采样点的人工合成地震记录，各道加入不同强度的随机噪声得到信噪比不同的含噪记录，信噪比范围由-8dB到0dB。该记录中包含三条直线型同相轴，其中两条相交于E1区域，在R2区域存在一段突变噪声，在R1区域有效信号淹没于强随机噪声当中，信噪比较低。图2(b)中的蝶型图中存在三个明显的蝶型中心点，分别对应地震记录中的三条同相轴。图2(c)是蝶型图的一个三维图，图中的三个明显尖峰就对应蝶型图中的三个蝶型中心点。从参数空间的三维立体图可以看出，累加器矩阵存在三个明显峰值，选取其中峰值最大的前三个，即可恢复出图像域内相应的三条同相轴包络线，如图2(d)中黑色线所示。

步骤二：根据检测到的同相轴包络斜率k_l建立与同相轴形态接近的滤波轨线方程

y_l＝k_lx_l+b_l, (5)

其中x_l为地震数据点的横坐标，代表炮检距方向，y_l为地震数据点的纵坐标，代表时间方向，k_l为轨线斜率，b_l为截距，决定轨线与地震数据纵坐标的交点位置；不断改变截距b_l的数值便会生成多条具有相同斜率的直线轨线，沿着这些轨线对地震数据进行样本重采样，便可实现对地震数据沿同相轴形态方向数据提取；

在重采样过程中由于位移-时间域内的样本格点与新作用域内的样本格点不一致，部分样本点不能落到具体的地震道上，如果丢弃这些“脱道”的点，采样得到的序列中就会出现一些断点和尖峰脉冲，本发明通过在采样前在相邻地震道间插入一些虚拟地震道的方式对地震数据进行扩充，使得原本“脱道”的点能够落在具体的地震道上，从而有效降低采样误差，如图3所示，图中虚线是滤波轨线，横坐标为地震道方向，纵坐标为时间方向，点P₁、P₃为时刻点t_i-1和t_i+1对应的轨线上的样本点，它们分别落在地震道x_i-1和x_i+2上，从图中可以看出，菱形点P₂在插值后可以落到了某一虚拟道上，在采样过程中它就不会被丢失，也不再利用其它点进行代替，使得采样得到的序列更准确；

采样过程首先在地震数据中选取一个样本点，设其纵坐标为将其代入式(5)中的轨线方程计算出相应的横坐标即可提取一个坐标为的数据点；再选取纵坐标代入式(5)中的轨线方程计算出相应的横坐标x_li+1，提取第二个样本点依次类推。直到计算的横坐标达到最大的地震道数为止，这样就可提取得到一个样本序列。接着依次提取第二个、第三个…样本序列。重采样后得到的是一个包含多条序列的数据矩阵；

步骤三：选取合适窗长，利用时频峰值滤波方法对沿多条轨线采样后的重采样数据进行噪声压制，由于对地震数据进行样本重采样时所采用的轨线斜率是通过检测同相轴形态得出的，因此轨线与同相轴趋势接近，沿着这些轨线重采样的地震数据中有效信号线性度充分提高，主频充分降低，因此对其进行时频峰值滤波时所采用的滤波窗长比传统时频峰值滤波方法中要长很多，较长的窗长在完整保持有效成分的同时可以有效地压制随机噪声。选取较大窗长对重采样后数据序列进行时频峰值滤波，实现对含噪地震记录的消噪处理。

下面给出具体的实验例子来说明本发明的实验效果：

图4(a)为某一地震剖面部分记录，其中包含大量随机噪声，有效同相轴淹没其中，该记录中同相轴大都呈现直线型形态分布。分别利用传统时频峰值滤波算法和基于改进霍夫变换的自适应时频峰值滤波算法对该含噪地震记录进行噪声压制处理，结果分别如图4(b)和图4(c)所示。通过观察两种方法的消噪结果可以看出，两种方法去噪后，随机噪声都得到了压制，同相轴变得清晰。但相比之下，本发明中方法消噪后，同相轴变得更加连贯，噪声压制效果更为理想。该实验结果验证了自适应轨线时频峰值滤波算法的有效性。

Claims (2)

1.一种基于改进霍夫变换的自适应时频峰值滤波消噪方法，其特征在于，包括下列步骤：

步骤一：一张含噪的、包含直线型同相轴的地震记录，记录中的反射同相轴相当于图像中的边缘信息，为了便于对地震数据进行形态识别，首先对地震数据进行灰度处理，将数据值转换到0-255之间，处理后的地震数据为一张图像数据，接着，利用改进的霍夫变换实现地震记录中同相轴形态识别，具体过程为：先利用公式(0≤α≤1)对灰度处理后的图像数据进行多层分数阶傅里叶变换，将地震数据变换到时域和频域之间的一个中间域，将其旋转到从时间轴指向频率轴的一个角度，得到多组不同阶数的傅里叶变换结果；再将变换结果转换成由模值ρ和相角θ构成的极坐标形式，在极坐标内，对数据进行共轭对称，并沿模值方向进行一维傅里叶反变换便可得到一张由多条正弦线构成的正弦谱图；实际上，从多层分数阶傅里叶变换到转换成极坐标，再到进行一维傅里叶反变换，整个过程相当于对地震数据进行了一次霍夫变换，得到的正弦谱图中多条正弦线的交点，即正弦谱矩阵中的数据最大值就对应了地震数据中的同相轴；对于直线型同相轴的地震记录，通过搜索正弦谱图中的峰值便可得到与其对应的地震数据中的某条同相轴在极坐标内的模值ρ_peak和相角θ_peak，再根据极坐标与原图像坐标之间的关系式ρ＝xcosθ+ysinθ便可以得到地震数据内同相轴包络线的斜率k_l＝ctanθ_peak；于是沿着多条斜率为k_l的直线型滤波轨线对含噪地震记录进行样本重采样，有效提高采样后信号的线性度，降低滤波偏差；

步骤二：根据检测到的同相轴包络斜率k_l建立与同相轴走势接近的滤波轨线方程y_l＝k_lx_l+b_l，其中x_l为地震数据点的横坐标，代表炮检距方向，y_l为地震数据点的纵坐标，代表时间方向，k_l为轨线斜率，b_l为截距，决定轨线与地震数据纵坐标的交点位置，改变截距b_l的数值便会生成多条具有相同斜率的直线轨线，沿着这些轨线对地震数据进行样本重采样，便可实现对地震数据沿同相轴形态方向的数据提取；

步骤三：对采样后的数据样本，利用时频峰值滤波方法进行噪声压制，由于对地震数据进行样本重采样时所采用的轨线斜率是通过检测同相轴形态得出的，因此轨线与同相轴趋势接近，沿着这些轨线重采样后的地震数据有效信号线性度充分提高，主频充分降低，因此对其进行时频峰值滤波时所采用的滤波窗长比传统时频峰值滤波方法中长很多，较长的窗长在完整保持有效成分的同时可以有效地压制随机噪声。

2.根据权利要求1所述的一种基于改进霍夫变换的自适应时频峰值滤波消噪方法，其特征在于，步骤一中对于直线型同相轴，地震数据中的某一点(x,y)与极坐标内的参数对(ρ，θ)满足如下关系：ρ＝xcosθ+ysinθ,θ∈[0,π]，其中θ为直线法线与x轴的夹角，ρ为图像空间中原点到直线的距离。