基于Butterworth 滤波器的线性噪声压制方法研究及实现

李子伟，邹明亮，曹成寅，吴曲波

（1. 核工业北京地质研究院, 北京100029； 2. 核工业二三〇研究所, 长沙 410007）

作为发展核电和巩固国防的战略资源,铀矿有着不同于其他矿产的特殊地位。 21 世纪初我国铀矿战略转移至北方, 此后的十余年落实和扩大了一批砂岩型铀矿床, 资源增长量很快, 在新增探明的资源量中约占4/5,也使得砂岩型铀矿在我国铀资源总量中所占的比例超过了40%, 因此我国砂岩型铀矿存在巨大的发展潜力[1-2]。

砂岩型铀矿勘查实践表明利用地球物理方法查明地层结构、 构造和砂体分布规律对于砂岩型铀矿的研究具有重要意义[2-6]。 在众多地球物理方法中, 地震勘探方法具有较大的探测深度和较高的探测精度[5-6], 是未来砂岩型铀矿地球物理勘探中不可或缺的主要方法之一。 随着勘探深度的增加, 砂岩型铀矿的地质条件越来越复杂, 这对地震资料的处理也提出了更高的要求, 其中地震资料的信噪比是保证地震资料品质的基础, 因此提高信噪比是地震数据处理的首要任务。

渥·伊尔马兹等人将地震资料中的噪声分为两类: 随机噪声和相干噪声, 其中相干噪声又包括线性噪声、 交混回响和多次波[7-9]。相干线性噪声在陆上资料中以频散瑞雷波的形式存在, 通常把它称为面波, 这些噪声的存在直接影响地震资料叠加剖面的成像效果,从而增加了地震资料解释的难度。 为满足当前砂岩型铀矿勘查的要求, 研究与开发压制线性噪声的技术方法显得十分重要。 笔者基于Butterworth 滤波器, 开展了频率-波数域（f-k 域） 的线性噪声压制方法研究, 并利用Fortran 语言编制了计算机程序, 实现了地震数据线性噪声的压制, 新疆准噶尔盆地的实际地震数据处理结果表明, 该方法能够较好地压制线性噪声, 有效提高地震资料的信噪比。

1 频率-波数域的数字滤波

在地震勘探中, 用数字仪记录地震波时,为了保持更多的地震波的特征, 通常利用宽频带进行数据记录, 因此, 在宽频带范围记录了各种有效反射波的同时, 也记录了各种干扰波。 有效波和干扰波的差异表现在多个方面, 例如频谱、 传播方向、 能量等, 在地震数字处理中利用频谱特征的不同来压制干扰波, 以突出有效波的方法就是数字滤波,是地震数据处理方法的基础[9]。

图1 频率-波数域（f-k 域）示意图Fig. 1 Schematic map of f-k domain

对于线性干扰而言, 它与有效信号的差异主要体现在频率（f）、 视速度（v）方面, 该差异在频率-波数域（f-k 域）易于区分, 如图1所示, 通过原点的两条直线的斜率就是视速度, I 区为高速干扰区, II 区是有效信号区,III 区是低速干扰区, 常见的面波干扰即位于该区域。 可以看出, 视速度存在差异的有效信号和干扰信号在f-k 域能有效区分, 因此利用f-k 域滤波可以压制各种线性干扰。

用X（t, x）、 Y（t, x）表示滤波前后的地震记录, 用H（t, x）表示滤波因子, 对地震记录X（t, x）施加滤波因子H（t, x）, 经处理后即可得到滤波后的地震记录Y（t, x）, 该时空域的滤波过程可以用图2 简要表示。

图2 时间-空间域滤波过程示意图Fig. 2 Schematic map of the time-space domain filtering process

以上时间-空间域滤波过程可以用二维褶积公式（1）表示:

根据褶积原理, 两个时间函数在时间域褶积的频谱等于两个函数的频谱在频率域的乘积, 因此式（1）的时间-空间域的二维褶积运算可以在频率-波数域完成, 即输出信号的频-波谱是输入信号的频-波谱与滤波器的频-波谱的乘积。

在进行滤波处理时, 通过二维傅里叶变换将时空域的地震数据变换到f-k 域,得到地震数据在f-k 域的频-波谱, 与构建的滤波器的频-波谱进行乘积运算, 即可得到输出信号的频-波谱, 再利用二维反傅里叶变换, 获得滤波后的时空域地震数据。

2 Butterworth 滤波器

在f-k 域数字滤波过程中, 滤波器选取的好坏, 直接决定了滤波的效果。 常规的扇形滤波器, 由于通带和阻带间尖锐边界的存在,会产生严重的吉布斯效应[8-9], 从而影响滤波后的地震信号的保真度。 为最大程度地消除吉布斯效应, Hale & Claerbout 提出了一种在时空域实现的倾角滤波方法---Butterworth 类倾角滤波方法[10], 该方法采用的Butterworth类倾角滤波器在通带与阻带间变化缓慢, 其算子带内、 带外幅频特性最佳平稳, 可以避免滤波后因吉布斯效应产生的断断续续的假同相轴现象。

考虑到砂岩型铀矿勘查中的地震资料的干扰主要以面波和声波为主, 在地震剖面上表现为较低视速度、 较高倾角的同相轴, 因此本文在Hale & Claerbout 提出的Butterworth滤波器基础上进行了简化, 仅采用高通倾角滤波器滤压制低视速度的噪声[10-11]。

Hale & Claerbout 提出的Butterworth 高通倾角滤波器在f-k 域的频谱响应可以用公式（3）表示:

式中: ω-弧度, ω＝2πf。

笔者对公式（3）进行简化后可用公式（4）表示:

根据LI Jianchao 等提出的滤波参数选择结果[12], 令α＝fNνN, 代入公式（4）, 可得到简化后的Butterworth 高通倾角滤波器在f-k域的频谱响应, 用公式（5）表示:

基于公式（5）表示的Butterworth 滤波器,利用Fortran 语言进行计算机实现, 完成f-k域的线性噪音压制。

3 线性噪声压制方法的实现

在前述数字滤波和Butterworth 滤波器的数学模型基础上, 构建了f-k 域的Butterworth滤波器, 在f-k 域通过施加该滤波器进行线性噪声的滤波处理, 并采用Fortran 语言编写了以上滤波程序的核心算法, 实现了线性噪声的压制, 达到了提高地震数据信噪比的目的,具体的计算机实现流程见图3。

图3 线性噪声压制方法的计算机实现流程及函数调用Fig. 3 Computer implementation process and function call of linear noise suppression method

在时空域和f-k 域之间的数据转换依赖常规的傅里叶变换算法, 此处不再赘述。 本算法的核心为Butterworth 滤波器在频率-波数域的应用, 此处简要介绍调用函数的基本情况。

函 数 定 义 -SUBROUTINE BUTDIP（IHEAD, GATHER, IPAM, TABLE, ISTAT）；函数功能-核心算法主子程序, 执行在频率-波数域的倾角滤波； 返回值-IHEAD-地震数据道集道头, GATHER-地震数据道集数据体,IPAM-滤波参数信息, TABLE-FFT 用余铉函数表, ISTAT-程序运行结果指示信息； IPAM滤波参数信息: fN-压制截频, vN-压制视速度, dx-道间距。

利用以上核心算法程序, 对读取的地震数据开展傅里叶变换, 变换至f-k 域后施加滤波器, 然后再利用二维反傅里叶变换, 获取滤波后的时空域地震数据。 根据线性噪声与有效信号之间的视速度差异, 确定输入压制截频参数（fN）和压制视速度（vN）, 调整Butterworth 滤波器参数, 压制指定速度的线性噪声, 进而实现不同类型线性噪声的压制,达到提高地震数据信噪比的目的。

4 线性噪声压制方法的应用效果对比

采用本文开发的地震数据处理程序, 对新疆准噶尔盆地的实际地震数据进行了处理测试, 经处理后的单炮记录和时间剖面的线性干扰得到了较好的压制, 突出了有效信号,提高了地震资料的信噪比, 这对后续的地震解疑具有重要的意义。

图4 展示了开发的处理程序在准噶尔盆地的单炮地震记录上的应用效果。 图4a 为原始的单炮地震记录, 可以看到, 在原始单炮记录上线性干扰严重, 双曲线特征的有效信号被淹没在线性干扰中。 采用同一组滤波参数（带通滤波15～150 Hz, 压制视速度2 000 m/s）,利用某商业软件的f-k 滤波模块和研发的滤波程序分别进行滤波处理, 结果分别如图4b 和图4c 所示。 图4b 为某商业软件的f-k 滤波模块的处理效果, 可以看到线性干扰受到一定程度的压制, 有效信号凸显, 但线性干扰仍然有所残留； 图4c 为开发的滤波程序的处理效果, 可以看到, 相对某商业软件的f-k 滤波效果, 开发的滤波程序对线性噪声的压制效果更佳, 经处理后的单炮记录有效信号更加明显。

图5 展示开发的处理程序在新疆准噶尔盆地地震叠加剖面上的应用效果。 图5a 为原始的地震叠加剖面, 可以看到, 在原始叠加剖面上线性干扰严重, 影响了有效信号的成像。 利用某商业软件的f-k 滤波模块和研发的滤波程序分别进行滤波处理, 结果如图5b 和图5c 所示。 图5b 为某商业软件的f-k 滤波模块的处理效果, 可以看到线性干扰受到一定程度的压制, 但仍然存在部分残留, 尤其是在剖面的右侧, 线性干扰仍然比较严重, 有效信号的成像受到了较大的影响； 图5c 为开发的滤波程序的处理效果, 可以看到, 相对某商业软件f-k 滤波效果, 开发的滤波程序对线性噪声的压制效果更佳, 经处理后的叠加剖面成像效果更好, 尤其是在剖面的右侧,有效信号的成像效果明显增强。

通过对比分析开发的滤波方法和某商业软件的f-k 滤波模块的去噪效果, 结果表明开发的基于Butterworth 滤波器的去噪方法能够有效地压制线性噪声, 效果优于常规商业软件的f-k 滤波模块。

图4 准噶尔盆地去噪前后的地震单炮记录对比Fig. 4 Comparison of seismic single-shot record before and after filtering in Junggar Basin

5 结论

开展了Butterworth 滤波器的研究, 并基于该滤波器开展了f-k 域的数字滤波研究, 利用Fortran 语言编制了计算机程序, 实现了线性噪声的压制。 利用新疆准噶尔盆地的地震数据进行了实际测试, 单炮记录和叠加剖面的测试结果表明, 地震数据的线性噪声得到了较好的压制, 突出了有效信号, 成像效果得到了较大改善, 有效地提高了地震资料的信噪比。

研究成果可推广应用到砂岩型铀矿地震资料的去噪处理中, 可以有效提高地震剖面的信噪比和成像效果, 有利于降低地震解释的难度, 保障地质解释的可靠性。