CN106680874A - 基于波形形态特征稀疏化建模的谐波噪声压制方法 - Google Patents

Abstract

本发明公开一种基于波形形态特征稀疏化建模的谐波噪声压制方法，包括：1)、根据可控震源滑动扫描采集地震记录中谐波噪声波形形态特征构造Chirplet变换，并与连续小波变换构成超完备字典；2)、Chirplet正、反变换的快速实现；3)、基于相关后数据时频分布特征确定Chirplet变换参数；4)、根据参考扫描信号的起始频率确定谐波噪声的滤波截止频率，实现有效信号与谐波噪声保真分离。本发明解决了可控震源滑动扫描采集地震数据中谐波噪声干扰的问题，达到了提高了地震资料信噪比的目的。本发明方法根据谐波噪声的时频分布特征确定选择Chirplet变换，保证了稀疏性；Chirplet变换的快速实现方式保证了变换的效率，修正系数保证了反变换重构的准确性；根据数据驱动自动确定Chirplet变换的参数，具有很强的适应性，且单道计算，便于并行处理。

Description

基于波形形态特征稀疏化建模的谐波噪声压制方法

技术领域

本发明属于地震勘探中数据处理技术领域，特别涉及一种地震采集数据中谐波噪声压制的方法。

背景技术

随着全球经济的快速发展以及环保意识的提高，传统的炸药震源逐渐被可控震源替代。滑动扫描是可控震源施工中常用的采集方式，提高了采集效率。但滑动扫描存在采集地震记录受到较强谐波噪声干扰的问题，有效信号往往掩藏在谐波噪声中，对数据资料的质量造成了严重影响。目前压制谐波噪声有如下方法：一是纯相移滤波法，但该方法处理相关前地震记录，数据量较大、计算效率低；二是分频异常振幅压制法，但该方法对有效信号会造成损伤，保真性不够；三是地面力信号法，但需要采集地面力信号，采集成本较高，且必须从未受影响的最后一炮数据处理，无法进行单炮集处理。

随着信号稀疏性理论的发展，Starck等人提出了形态成分分析的混合信号分解方法。形态成分分析是指，根据复杂信号的组成成分波形形态特征，将两种具有不同原子特征的变换字典构成超完备字典，实现对复杂信号更稀疏的表示方式以及更有效的信息识别能力，实现两种成分的分离。字典通常是根据经验从已知的数学变换中选择或者构造，而选出的字典是否充分满足形态成分分析的假设，是形态成分分析方法能否成功的关键。使用形态成分分析的方法对信号与谐波噪声实现分离来达到压制谐波噪声的目的，需要根据有效信号与谐波噪声的波形形态特征构造合适的稀疏变换。已有文献研究选择使用连续小波变换来稀疏表示有效信号比较合适。但谐波噪声较为复杂，需要构造较为合适匹配的变换来稀疏表示，并选择确定变换的具体参数，目前相关研究较少。

针对上述问题，本发明利用形态成分分析理论，给出了针对谐波噪声的波形形态特征，构造合适的稀疏表示变换即Chirplet变换，与连续小波变换构成超完备字典，并根据数据特征选择确定Chirplet变换参数，以及使用坐标分块松弛算法方法实现信噪分离达到压制谐波噪声的解决方案。

发明内容

本发明的目的在于提供了一种基于波形形态特征稀疏化建模的谐波噪声压制方法。针对谐波噪声的波形形态特征，构造合适的稀疏表示变换即Chirplet变换，与连续小波变换联合构成超完备字典，并根据数据特征选择确定Chirplet变换参数，以及使用坐标分块松弛算法方法实现信噪分离以达到有效压制可控震源滑动扫描记录的地震数据中的谐波噪声的目的。

为了实现上述目的，本发明采用如下技术方案：

基于波形形态特征稀疏化建模的谐波噪声压制方法，包括以下步骤：

1)、根据可控震源滑动扫描采集地震记录中谐波噪声波形形态特征构造Chirplet变换，并与连续小波变换联合构成超完备字典；

2)、Chirplet正、反变换的快速实现；

3)、基于相关后数据时频分布特征确定Chirplet变换参数；

4)、根据参考扫描信号的起始频率确定谐波噪声的滤波截止频率，实现有效信号与谐波噪声保真分离。

进一步的，根据可控震源滑动扫描采集地震记录中谐波噪声波形形态特征构造Chirplet变换，并与连续小波变换构成超完备字典，包括：

形态成分分析的对象是含有两种具有不同形态特征的成分：

式中：表示待分析信号；表示信号中的两种成分，具有不同的形态特征；形态成分分析的目标是分别提取出两种成分；假设和可以分别由字典A₁和A₂有效的稀疏表示，但是用A₂稀疏表示和用A₁稀疏表示时稀疏性差；

在地震数据处理中，一般选择连续小波变换作为稀疏表示有效信号成分的变换字典，其中连续小波变换为：

式中WT_x(A,τ)为变换系数，A表示尺度因子，x(t)表示待分析信号，ψ(t)表示Morlet母小波；其中t为时间，τ为平移量，*表示共轭；

连续小波变换的反变换为：

式中常数C_Ψ＜∞为其容许条件；

本发明中，构成合适的Chirplet变换作为稀疏表示谐波噪声的字典，其中Chirplet正变换的定义为：

CT_x(t,f,a,c,d)＝∫x(τ)h^*(τ-t,f,a,c,d)dτ，

式中CT_x(t,f,a,c,d)表示待分析信号的Chirplet变换值，x(t)表示待分析信号，a为伸缩参数，c为时间上的线性调频率，f为频率，d为频率上的线性调频率；τ为时间平移量，*表示共轭；h(t,f,a,c,d)为Chirplet变换的核函数，其表达式如下：

式中g(t)为高斯函数；v为平移量；

忽略频率上的线性调频操作，核函数退化为：

根据形态成分分析理论，用上面选定的字典A₁即连续小波变换和A₂即Chirplet变换，构成超完备字典，稀疏表示信号计算稀疏表示系数：

式中：x₁为重构系数中与A₁对应的部分；x₂为重构系数中与A₂对应的部分。为拉格朗日乘子；该优化问题可以通过坐标分块松弛算法求解。坐标分块松弛算法的基本思想是交替迭代的计算x₁和x₂。其主要内容步骤为：

初始化：初始迭代步数k＝0，初始解

迭代：每步迭代k增加1，并计算：

式中，T_λ为硬阈值函数；与A₁构成一对正反变换，与A₂构成一对正反变换；

终止条件：当小于预设的值时，即继续迭代对结果的影响足够小时，迭代终止。

输出：

进一步的，步骤2)即Chirplet正、反变换的快速实现中，具体包括：

步骤101：从原始数据中读出单道时域信号，记为x；

步骤102：计算得到信号x的Chirplet正变换结果，记为y，其实现方法类似于短时傅里叶变换；

步骤103：计算y的Chirplet反变换结果，记为

由于步骤103中Chirplet反变换得到的重构结果与原始信号相比会在幅度上存在一定的差别，因此接下来使用一个模型信号计算修正系数来进行幅度修正；

步骤201：构造一个线性升频信号模型，记为s，其长度与信号x相同；

步骤202：计算得到信号s的Chirplet正变换结果，记为s₀，其实现方法类似于短时傅里叶变换；

步骤203：计算得到s₀的Chirplet反变换结果，记为

步骤204：根据模型反变换结果计算修正系数系数，即其中||.||₂表示二范数；

步骤104：将修正系数r作用到单道时域信号x的Chirplet反变换得到修正的Chirplet反变换结果，记为则

进一步的，步骤3)即相关后数据的时频分布特征计算Chirplet变换的参数包括：

地震资料数据是阵列信号，利用形态成分分析实现有效信号与谐波噪声分离的方法基于单道数据处理，而由于各道数据的谐波噪声强弱不同，因此需要根据数据计算确定Chirplet变换的参数才能保证对谐波噪声稀疏，更好的实现信噪分离。其中具体实现步骤如下：

步骤301：将单道时域信号x转换到频域，记为F；

步骤302：确定频率域的高、低频分界，记为κ；

步骤303：计算高频(频率分界以上)能量占总频带能量的比值，即振幅谱比值，记为b；

步骤304：根据振幅谱比值计算Chirplet变换的伸缩参数，记为a，线性调频参数，记为c以及阈值权系数参数，记为λ：

a＝k₁×b，λ＝k₂-b，c＝-k₃×b×a²，

其中，b为振幅谱比值；k₁,k₂,k₃分别为参数a、c与λ的修正系数；高低频分界以及修正系数需在分析地震资料数据的特点之后合理给出。

进一步的，步骤4)即根据参考扫描信号的起始频率确谐波噪声的滤波截止频率，实现有效信号与谐波噪声保真分离具体包括：

步骤401：根据参考扫描信号的起始频率确定谐波噪声的起始频率f_b；

步骤402：谐波噪声做截止频率为f_b的高通滤波；

步骤403：原始数据减去滤波后的谐波噪声得到有效信号，实现有效信号与谐波噪声的保真分离。

相对于现有技术，本发明具有以下有益效果：

1)本发明方法根据谐波噪声的时频分布特征构造合适的变换，保证了变换对谐波噪声的稀疏性；

2)Chirplet正、反变换的快速实现方式保证了变换的计算效率，使其可用于海量实际资料处理；

3)Chirplet反变换中使用幅度修正系数使得反变换重构的准确性得到保证；

4)根据数据驱动自动确定Chirplet变换的参数，具有很强的适应性，且单道计算，易于并行实现；

5)本发明对有效信号具有高保真性，能够有效的保护有效波的高、低频成分。

基于本发明的时频域稀疏优化去噪方法解决了可控震源滑动扫描采集地震数据中谐波噪声干扰的问题，达到了提高了地震资料信噪比的目的。

附图说明

图1为谐波噪声的时域波形图；

图2为Chirplet变换原子示意图；

图3为Chirplet正、反变换的快速实现流程图；

图4为根据相关后数据时频分布特征确定Chirplet变换参数流程图；

图5为基于参考扫描信号起始频率的有效信号与谐波噪声保真分离流程图；

图6A为未受谐波干扰模拟数据；图6B为受谐波干扰模拟数据；

图7A本发明方法从图6B中分离的有效信号；图7B本发明方法从图6B中分离的谐波噪声；

图8A为实际滑动扫描地震记录图；

图8B为本发明方法压制图8A谐波噪声后剖面；图8C为本发明方法获得的图8A谐波噪声剖面；

图9A为图8A局部放大显示结果；图9B为图8B局部放大显示结果；图9C为图8C局部放大显示结果。

具体实施方式

为了使本发明的目的、技术方案更加清楚明白，下面结合附图和具体实施方式对本发明做进一步详细的说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为本发明的限定。

在本发明实施例中，提出了一种基于波形形态特征稀疏化建模的谐波噪声压制方法，包括以下内容：

1)、根据可控震源滑动扫描采集地震记录中谐波噪声波形形态特征构造Chirplet变换，并与连续小波变换构成超完备字典；

2)、Chirplet正、反变换的快速实现方法；

3)、基于相关后谐波噪声时频分布特征确定Chirplet变换参数；

4)、根据参考扫描信号的起始频率确定谐波噪声的滤波截止频率，实现有效信号与谐波噪声保真分离。

步骤1)中根据可控震源滑动扫描采集地震记录中谐波噪声波形形态特征构造Chirplet变换，并与连续小波变换构成超完备字典包括：

选择连续小波变换作为稀疏表示有效信号成分的变换字典：

式中WT_x(A,τ)为变换系数，A表示尺度因子，x(t)表示待分析信号，ψ(t)表示Morlet母小波；其中t为时间，τ为平移量，*表示共轭。

连续小波变换的反变换为：

式中常数C_Ψ＜∞为其容许条件。

Chirplet正变换的定义为：

CT_x(t,f,a,c,d)＝∫x(τ)h^*(τ-t,f,a,c,d)dτ，

式中CT_x(t,f,a,c,d)表示待分析信号的Chirplet变换值，x(t)表示待分析信号，a为伸缩参数，c为时间上的线性调频率，f为频率，d为频率上的线性调频率；τ为时间平移量，*表示共轭；h(t,f,a,c,d)为Chirplet变换的核函数，其表达式如下：

式中g(t)为高斯函数；v为平移量。

忽略频率上的线性调频操作，核函数退化为：

选择构造的Chirplet变换作为稀疏表示谐波噪声的字典。

如图1与图2所示，退化后的Chirplet变换原子与谐波噪声波形形态特征比较匹配，能更有效表示谐波噪声。

根据形态成分分析理论，用上面选定的字典A₁即连续小波变换和A₂即Chirplet变换构成超完备字典稀疏表示信号计算稀疏表示系数：

式中：x₁为重构系数中与A₁对应的部分；x₂为重构系数中与A₂对应的部分。λ为拉格朗日乘子。该优化问题可以通过坐标分块松弛算法求解。

如图3所示，步骤2)中Chirplet正、反变换的快速实现流程包括如下步骤：

步骤101：从原始数据中读出单道时域信号，记为x；

步骤102：计算得到信号x的Chirplet正变换结果，记为y，其实现方法类似于短时傅里叶变换；

步骤103：计算y的Chirplet反变换结果，记为

由于步骤103中Chirplet反变换得到的重构结果与原始信号相比会在幅度上存在一定的差别，因此接下来使用一个模型信号计算修正系数来进行幅度修正。

步骤201：构造一个线性升频信号模型，记为s，其长度与信号x相同；

步骤202：计算得到信号s的Chirplet正变换结果，记为s₀，其实现方法类似于短时傅里叶变换；

步骤203：计算得到s₀的Chirplet反变换结果，记为

步骤204：根据模型反变换结果计算修正系数系数，即其中||.||₂表示二范数；

步骤104：将修正系数r作用到单道时域信号x的Chirplet反变换x，得到修正的Chirplet反变换结果，记为则

通过本发明的Chirplet快速正反变换，计算效率较高，且修正系数使得重构精确度更高。

如图4所示，步骤3)中基于相关后数据时频分布特征确定Chirplet变换参数的方法包括如下步骤：

步骤301：将单道时域信号x转换到频域，记为F；

步骤302：确定频域的高低频分界，记为κ；

步骤303：计算高频(频率分界以上)能量占总频带能量的比值，即振幅谱比值，记为b；

步骤304：根据振幅谱比值计算Chirplet变换的伸缩参数，记为a，线性调频参数，记为c以及阈值权系数参数，记为λ：

a＝k₁×b，λ＝k₂-b，c＝-k₃×b×a²，

其中，b为振幅谱比值；k₁,k₂,k₃分别为参数a、c与λ的修正系数。

本发明实施例中，高、频分界选择为40Hz，k₁,k₂,k₃分别取：k₁＝2.0,k₂＝1.5,k₃＝5。

通过本发明的确定Chirplet变换参数的方法，能够各道数据得到不同的参数，从而具有较强的适应性，可以处理含不同强弱程度谐波噪声的地震记录数据。

如图5所示，步骤4)中基于参考扫描信号起始频率的谐波噪声保真分离方法包括如下步骤：

步骤401：根据参考扫描信号的起始频率确定谐波噪声的起始频率，记为f_b；

步骤402：谐波噪声做截止频率为f_b的高通滤波；

步骤403：原始数据减去滤波后的谐波噪声得到有效信号，实现有效信号与谐波噪声的保真分离。

使用基于本发明的保真分离方法，能够更加有效地保护有效信号，降低对有效信号的损伤。

在本例中，提出了一种基于波形形态特征稀疏化建模的谐波噪声压制方法，从而达到了有效压制可控震源滑动扫描记录地震数据中谐波噪声的目的。本发明具有如下有益效果：

1)本发明方法根据谐波噪声的时频分布特征构造合适的变换，保证了稀疏性；

2)Chirplet正、反变换的快速实现方式保证了变换的效率，使其可用于海量实际资料处理；

3)Chirplet反变换使用修正系数保证了结果的准确性；

4)根据数据驱动自动确定Chirplet变换的参数，具有很强的适应性，且单道计算，便于并行处理；

5)本发明对有效信号具有高保真性，能够有效的保护有效波的高、低频成分。

下面结合一个具体实施例对上述变换构造及参数确定方法进行具体说明，且该实施例仅是为了更好说明本发明，并不构成对本发明的不当限定。

利用本发明提供的分析处理方法应用到合成模拟与实际滑动扫描地震记录数据中验证本发明构造Chirplet变换及参数的方法压制谐波噪声的有效性。本发明的方法不仅能够有效的压制谐波噪声，而且有效信号具有较高保真性。

图6A-图6B为针对滑动扫描地震记录特性，给出的未受谐波干扰的合成记录以及受谐波干扰的合成记录。该模型包含三个反射层，所用扫描信号为线性升频扫描信号，其最低频率为3Hz，最高频率为90Hz，扫描信号时长16s，有效信号包含直达波以及反射波，谐波噪声包含二次谐波和三次谐波。有效信号受到很强的谐波噪声干扰，尤其是深层反射波受到谐波噪声干扰更为突出。

如图7A-图7B所示，本发明方法对该含谐波噪声的炮记录(图6B)进行处理，得到对应的有效信号剖面和对应的谐波噪声剖面。本发明方法可以有效的对谐波噪声进行压制，谐波噪声中几乎没有残留有效信号能量，表明本发明的方法能够有效压制谐波噪声，而且对有效信号具有高保真性。

图8A为我国西部某油田基于滑动扫描观测到的相关后炮集数据，接下来验证本发明方法处理实际地震资料的有效性。该炮集数据总共401道，采样时间间隔为2ms，记录长度为6s。从剖面上可以看出，该记录中的信号受到很强的谐波噪声干扰，导致有效信号被谐波噪声覆盖，对地震资料的分析以及解释造成严重影响。

如图8B-图8C所示，应用本发明的方法对该炮集数据进行处理获得其有效信号和谐波噪声，本发明方法有效地压制了地震记录中的谐波噪声。并且如图8B中方框区域指示处被谐波噪声覆盖的有效信号在有效信号剖面中清晰地显示。

对矩形框区域进行放大显示，如图9A-9C所示，在谐波噪声剖面中几乎不含有效信号能量，表明本发明的方法对有效信号具有较高保真性。

以上的模型和实际资料算例中，利用本发明的方法构造谐波噪声的稀疏表示变换及确定变换的参数，对地震资料数据进行谐波噪声压制，不仅能够有效的压制谐波噪声，而且有效信号具有较高保真性，为后续资料的分析奠定基础。

最后需要说明的是，以上模型和实际资料算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围。凡在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。

Claims (5)

1.基于波形形态特征稀疏化建模的谐波噪声压制方法，其特征在于，包括以下步骤：

1)、根据可控震源滑动扫描采集地震记录中谐波噪声的波形形态特征构造Chirplet变换，并与连续小波变换构成超完备字典；

2)、Chirplet正、反变换的快速实现；

3)、基于相关后数据时频分布特征确定Chirplet变换参数；

4)、根据参考扫描信号的起始频率确定谐波噪声滤波截止频率，实现有效信号与谐波噪声保真分离。

2.如权利要求1所述的基于波形形态特征稀疏化建模的谐波噪声压制方法，其特征在于，步骤1)，包括：

形态成分分析的对象是含有两种具有不同形态特征的成分：

$\stackrel{\·}{s}={\stackrel{\·}{s}}_1+{\stackrel{\·}{s}}_2,$

式中：表示待分析信号；表示信号中的两种成分，具有不同的形态特征；形态成分分析的目标是分别提取出两种成分；假设和分别由字典A₁和A₂有效的稀疏表示，但是用A₂稀疏表示和用A₁稀疏表示时稀疏性差；

在地震数据处理中，一般选择连续小波变换作为稀疏表示有效信号成分的变换字典，其中连续小波变换为：

${\mathrm{WT}}_x(A,\mathrm{\τ})=\frac{1}{\sqrt{A}}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}x\left(t\right){\mathrm{\ψ}}^{*}\left(\frac{t-\mathrm{\τ}}{A}\right)dt,$

式中WT_x(A,τ)为变换系数，A表示尺度因子，x(t)表示待分析信号，ψ(t)表示Morlet母小波；其中t为时间，τ为平移量，*表示共轭；

连续小波变换的反变换为：

$x\left(t\right)=\frac{1}{C_{\mathrm{\ψ}}}\underset{0}{\overset{+\mathrm{\∞}}{\∫}}\frac{dA}{A}\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{\mathrm{WT}}_x(A,\mathrm{\τ})\frac{1}{A^2}\mathrm{\ψ}\left(\frac{t-\mathrm{\τ}}{A}\right)d\mathrm{\τ},$

式中常数C_Ψ＜∞为其容许条件；

构成合适的Chirplet变换作为稀疏表示谐波噪声的字典，其中Chirplet正变换的定义为：

CT_x(t,f,a,c,d)＝∫x(τ)h^*(τ-t,f,a,c,d)dτ，

式中CT_x(t,f,a,c,d)表示待分析信号的Chirplet变换系数，x(t)表示待分析信号，a为伸缩参数，c为时间上的线性调频率，f为频率，d为频率上的线性调频率；τ为时间平移量，*表示共轭；h(t,f,a,c,d)为Chirplet变换的核函数，其表达式如下：

$h(\mathrm{\τ}-t,f,a,c,d)={(-jda)}^{\frac{1}{2}}\∫g(\frac{\mathrm{\τ}-t}{a}-v)e^{j\mathrm{\π}(2f\mathrm{\τ}+c{\left(\frac{\mathrm{\τ}-t}{a}-v\right)}^2+\frac{v^2}{d})}dv,$

式中g(t)为高斯函数；v为平移量；忽略频率上的线性调频操作，核函数退化为：

$h(t,f,a,c)=g\left(\frac{t}{a}\right)e^{j\mathrm{\π}2ft+j\mathrm{\π}c{\left(\frac{t}{a}\right)}^2},$

根据形态成分分析理论，用上面选定的字典A₁和A₂，构成超完备字典，稀疏表示信号计算稀疏表示系数：

$\{x_1^{*},x_2^{*}\}=\arg \underset{x_1,x_2}{min}\stackrel{\·}{\mathrm{\λ}}\left(\right|\left|x_1\right|{|}_1+\left|\right|x_2\left|{|}_1\right)+\frac{1}{2}\left|\right|\stackrel{\·}{s}-A_1{x}_1-A_2{x}_2|{|}_2^2,$

式中：x₁为重构系数中与A₁对应的部分；x₂为重构系数中与A₂对应的部分；为拉格朗日乘子；该优化问题通过坐标分块松弛算法求解。

3.如权利要求1所述的基于波形形态特征稀疏化建模的谐波噪声压制方法，其特征在于，步骤2)具体包括：

步骤101：从原始数据中读出单道时域信号，记为x；

步骤102：计算得到信号x的Chirplet正变换的结果，记为y；

步骤103：计算得到y的Chirplet反变换结果,记为

步骤201：构造一个线性升频信号模型，记为s，其长度与信号x相同；

步骤202：计算得到信号s的Chirplet正变换结果，记为s₀；

步骤203：计算得到s₀的Chirplet反变换结果，记为

步骤204：根据模型反变换结果计算修正系数，即其中||.||₂表示二范数；

步骤104：将修正系数r作用到单道时域信号x的Chirplet反变换得到修正的Chirplet反变换结果，记为则

4.如权利要求1所述的基于波形形态特征稀疏化建模的谐波噪声压制方法，其特征在于，步骤3)包括：

步骤301：将单道时域信号x转换到频域，记为F；

步骤302：确定频域的高、低频分界，记为κ；

步骤303：计算频率分界以上的高频能量占总频带能量的比值，即振幅谱比值，记为b；

步骤304：根据振幅谱比值计算Chirplet变换的伸缩参数记为a，线性调频参数记为c以及阈值权系数参数记为λ：

a＝k₁×b,λ＝k₂-b,c＝-k₃×b×a²，

其中，b为振幅谱比值；k₁,k₂,k₃分别为参数a、c与λ的修正系数。

5.如权利要求1所述的基于波形形态特征稀疏化建模的谐波噪声压制方法，其特征在于，步骤4)具体包括：

步骤401：根据参考扫描信号的起始频率确定谐波噪声的起始频率f_b；

步骤402：谐波噪声做截止频率为f_b的高通滤波；

步骤403：原始数据减去滤波后的谐波噪声得到有效信号，实现有效信号与谐波噪声的保真分离。