一种自动衡量地震数据谐波噪声强弱的振幅谱比值方法与流程

本发明属于地震勘探数据处理领域，特别涉及一种谐波噪声强弱的自动判别方法。

背景技术：

可控震源是不同于爆炸震源的另一种形式的激发源，不但可以降低成本，还可以提高效率。目前，滑动扫描是可控震源的一种高效采集方式。滑动扫描是指下一组震源不必等上一组震源震动结束即可开始震动的一种采集方法。滑动扫描采集技术能够得到高密的炮点、高覆盖次数的地震数据，但采用滑动扫描采集的最大问题是谐波畸变。谐波噪声对采集信号造成了污染，降低资料信噪比，严重时会掩盖掉有效信息，影响对采集到的数据的分析使用。

因此，具有良好处理效果的去谐波方法是滑动扫描采集方式继续推向前进的关键。同时，滑动扫描采集方式每天可以得到近万炮的大量数据，这就要求去谐波方法具有较高的处理效率以及需要对不同类型记录数据要有更广泛的适用性，才能在工业上得到推广使用。在高效率压制谐波噪声过程中，谐波噪声强弱是一个十分重要的指标，影响各种处理参数的选择，因此研究一种可以自动判断地震采集信号中谐波噪声强弱的方法十分重要，具有巨大的理论意义与市场价值。

现有技术1：

经验法。该方法主要是从业者根据地震数据处理经验，通过观察地震数据的纹理走向，来判断确定地震数据谐波噪声的强弱。地震数据有效信号表现为明显的线性纹理，尤其是明显的同相轴位置，但是谐波噪声的纹理表现为剧烈的振荡波纹状。

现有技术1的缺点：

1、该方法判断确定的指标，属于一个主观估计值，不够精确，在准确性方面有所欠缺。

2、该方法针对多道数据确定一个综合强弱指标，无法针对每一道数据确定强弱。

现有技术2：

深浅层能量比值法。该方法首先确定初至波的位置，然后向下延伸一段数据作为深浅部分的分割线，计算深层数据的能量与浅层数据的能量比值作为衡量地震数据中谐波噪声强弱的指标。

现有技术2的缺点：

1、该方法通过时域能量比值的方法确定，对谐波噪声的位置要求较高，适用性不强。

2、该方法针对不同采样长度的数据确定深浅层分界时的参数不同，灵活性较差。

技术实现要素：

本发明的目的在于提供一种自动衡量地震数据谐波噪声强弱的振幅谱比值方法，以解决上述技术问题。本发明通过对单道地震数据做傅里叶变换，在频率域选择合适的分界计算高频能量占总频带能量的比值，并将此比值作为衡量谐波噪声强弱的指标。这种非人工的自动判断方法提高了判断衡量单道数据谐波噪声强弱的精确性，是自适应压制谐波噪声的基础。

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

一种自动衡量地震数据谐波噪声强弱的振幅谱比值方法，包括以下步骤：

步骤01：读取原始地震记录信号的单道数据；

步骤02：对步骤01读取的单道数据做快速傅里叶变换，从时域转换到频率域；

步骤03：确定高低频频谱分界，计算高频能量占总频带能量的比值；

重复步骤01-03直到所有道数据处理完成。

进一步的，步骤02中对地震资料单道数据做快速傅里叶变换，具体包括：

假设x(t)表示待分析信号，X(ω)为x(t)的频谱，使用快速傅里叶变换求取频谱：

X(ω)＝FFT(x(t))。

假设待分析信号为离散信号，长度为N，则离散的频谱长度也为N；由于FFT得到的频谱具有对称性，因此只考虑前N/2个离散点频谱；

假设第i点的频率值为ω_i，则有下式：

其中，dt为采样间隔，则为采样频率；采样频率记为f_N，采样频率一半记为f_N/2。

进一步的，步骤03中高低频频谱分界为40Hz。

进一步的，步骤03中计算高频能量占总频带能量的比值，包括：

选择40Hz作为低频与高频分界，接下来计算高频段40Hz-f_N/2的能量，记为f_high：

其中，i₄₀是40Hz频率点的位置。

计算总频段0-f_N/2的能量，记为f_all：

计算高频段能量与总频段能量的比值，记为α：

α＝f_high/f_all。

进一步的，步骤03中在计算高频段能量时将100Hz以上的频率点舍弃掉；此时，f_high、f_all的改进形式为：

其中，

本发明称根据振幅谱比值来衡量可控震源滑动扫描记录地震资料数据中谐波噪声强弱的方法为“振幅谱比值法”。

相对于现有技术，本发明具有以下有益效果：利用本发明的自动衡量谐波噪声强弱的方法分别对强谐波干扰与弱谐波干扰数据进行处理，得到的振幅谱比值指标能够很好地衡量谐波噪声的强弱。同时在振幅谱比值的基础上，根据此振幅谱比值可以确定稀疏优化方法的字典参数，取得了较好的结果。

附图说明

图1A为含有谐波噪声的实际地震资料数据；图1B为抽取图1A第300道数据的时频图；

图2A为单道待分析信号的振幅谱；图2B为单道待分析信号的振幅谱左半部分；

图3A为实际地震资料数据；图3B为图3A所示数据的各道深浅层能量比值；图3C为图3A所示数据的各道振幅谱比值；

图4为实际地震数据；

图5A为图4所示数据的局部区域放大图；图5B为基于深浅层能量比值法压制图5A所示数据中谐波噪声的结果；图5C为基于振幅谱比值法压制图5A所示数据中谐波噪声的结果；

图6为图4所示数据抽取第240道信号基于深浅层能量比值法去噪后与基于振幅谱比值法去噪后的时域波形图；

图7A为图4所示数据抽取第240道信号基于深浅层能量比值法得到的有效信号时频图；图7B为图4所示数据抽取第240道信号基于振幅谱比值法得到的有效信号时频图；

图8为本发明流程图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步详细的说明。

本发明为一种自动衡量地震数据谐波噪声强弱的振幅谱比值方法，通过将单道地震数据从时间域转换到频率域，在频率域选择合适的频率分界计算高频能量占总频带能量的比值，并将此比值作为衡量谐波噪声强弱的指标。

请参阅图8所示，本发明为一种自动衡量地震数据谐波噪声强弱的振幅谱比值方法，对地震记录信号实施步骤01-步骤03，包括：

步骤01：读取原始地震记录信号的单道数据；

步骤02：对步骤01读取的单道数据做快速傅里叶变换，从时域转换到频率域；

步骤03：确定高低频频谱分界，计算高频能量占总频带能量的比值；

重复步骤01-03直到所有道数据处理完成。

步骤02中对地震资料单道数据做快速傅里叶变换，具体包括：

假设x(t)表示待分析信号，X(ω)为x(t)的频谱，使用快速傅里叶变换求取频谱：

X(ω)＝FFT(x(t))；

假设待分析信号为离散信号，长度为N，则离散的频谱长度也为N；由于FFT得到的频谱具有对称性，因此只考虑前N/2个离散点频谱；

假设第i点的频率值为ω_i，则有下式：

其中，dt为采样间隔，则为采样频率；采样频率记为f_N，采样频率一半记为f_N/2。

步骤03中在频率域选择高低频频谱分界：

谐波噪声主要相对表现为高频能量，有效信号主要相对表现为低频能量；比如，二次谐波的时频曲线是有效信号(基波)的二倍，三次谐波的时频曲线是有效信号(基波)的三倍，更高次的谐波的频率会表现更高；

一般的可控震源滑动扫描带来的谐波噪声的频率都比有效信号频率高，虽然二者的频带存在交叉，这也是无法直接使用滤波器压制谐波噪声的原因，但是，通过在频率域选择一个分界点，可以近似认为谐波噪声能量主要集中在分界点以上，有效信号能量主要集中在分界点以下。

步骤03中计算高频能量占总频带能量的比值，包括：

本发明中选择40Hz作为低频与高频分界，接下来计算高频段(40Hz-f_N/2)的能量，记为f_high：

其中，i₄₀是40Hz频率点的位置。

计算总频段(0-f_N/2)的能量，记为f_all：

计算高频段能量与总频段能量的比值，记为α：

α＝f_high/f_all。

在实际数据处理中，待分析信号的频带一般不会超过100Hz，因此在这种情况下，在计算高频段能量时可以将100Hz以上的频率点舍弃掉，不影响处理结果，并且会提高计算效率；此时，f_high、f_all的改进形式为：

其中，

根据振幅谱比值α的定义，可知此振幅谱比值α越大，说明该道地震数据中谐波噪声越强。

如图1A所示，该片实际叠后地震数据受到较强谐波噪声的干扰，尤其在图中矩形框区域影响最为明显。谐波噪声覆盖了有效信号的同相轴，对下一步分析数据造成了干扰。

如图1B所示，为抽取图1A所示数据抽取第300道的时频图。可以看到在频率上，有效信号与谐波噪声相比，主要表现为低频，谐波噪声表现为高频。

如图2A所示，对时域待分析信号做快速傅里叶变换，得到的频率域最大频率点500Hz即为采样频率，也就是采样间隔0.002秒的倒数。同时，得到的振幅谱关于采样频率的一半即250Hz点对称，这符合奈斯奎特采样定律。图2B是图2A的左半部分，因此只分析图2B中所包含的频段范围即可。如图2B所示，大于100Hz的频段几乎为零，因此在计算振幅谱比值时可以用40Hz-100Hz来代替40Hz-250Hz，以减少运算量。

如图3A所示，是一片实际地震资料，共480道数据，采样点长度为3501，采样间隔为0.002秒。此地震数据在150至300道4秒至6秒区域受到了强谐波噪声干扰。

如图3B-图3C所示，分别是得到的图3A地震资料每道数据的深浅层能量比值法比值与振幅谱比值。可以看出，两者都能基本反映图3A中各道数据的谐波强弱，但图3B深浅层能量比值法比值曲线有快速上升与快速下降，没有图3C振幅谱比值法曲线平滑，同时，图3B深浅层能量比值法根据0.1作为分界线分为两段，每段分别取得一组字典参数，而振幅谱比值法是根据比值来线性得到强弱指标的，因此振幅谱比值法在反映谐波噪声强弱时更加稳定。

本发明具有如下有益效果：

1)本发明方法针对单道数据计算振幅谱比值，作为衡量谐波噪声强弱的指标，而不是一片数据确定一个指标，因此具有很强的自适应性；

2)本发明方法使用快速傅里叶变换，效率高；

3)本发明方法在计算上，各道数据之间没有联系，可以并行处理，便于在工业上大规模应用；

下面将基于本发明的时频域稀疏优化压噪方法应用到实际滑动扫描地震记录中对有效信号和谐波噪声进行分离。应用结果表明，本发明相比较于其他方法，对谐波噪声强弱的判断更加准确，压制谐波噪声效果也更好。

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

如图5A为图4所示的地震资料为便于进一步观察，在200-300道4秒-6秒截取的一个时空窗。

如图5B所示，通过深浅层能量比值法判断地震资料数据(图5A)中谐波噪声的强弱，并根据该指标进行谐波噪声压制，得到对应的有效信号剖面。该方法压制了浅层的谐波噪声，但同时资料数据深层残留了大量谐波噪声。

如图5C所示，通过本发明方法判断地震资料数据(图5A)中谐波噪声强弱，并根据该指标进行谐波噪声压制，得到对应的有效信号剖面。本发明方法对深层与浅层的谐波噪声可以同时进行有效的压制。

可以看到，基于本发明方法的谐波噪声压制去噪效果更加彻底，得到的有效信号数据更加干净。

如图6所示，提取第240道原始数据，基于深浅层能量比值法的噪声压制得到的有效信号残留有明显的谐波噪声，而基于本发明方法的噪声压制将谐波噪声和有效信号进行了分离，有效地压制了谐波噪声。

如图7A-图7B所示，提取的第240道有效信号的时频图，基于深浅层能量比值法上的噪声压制得到的有效信号残留有较为明显的谐波噪声，而基于本发明方法的噪声压制去除谐波噪声效果更加彻底，有效信号没有谐波噪声残留。

以上的实际资料算例中，利用本发明的振幅谱比值方法可以自动衡量可控震源滑动扫描记录地震资料数据中的谐波噪声强弱，该判断指标为后续资料的分析奠定基础，同时本发明方法具有很强的自适应性。

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