利用MATLAB信号处理工具箱SPTool完成地震资料滤波处理

2013-12-27孙炜

防灾减灾学报 2013年3期

孙炜

（德都地震台，黑龙江五大连池 164155）

0 引言

数字滤波器是一种用来过滤时间离散信号的数字系统，它是通过对抽样数据进行数学处理来达到频域滤波的目的。数字滤波器可以用软件或设计专用的数字处理硬件两种方式来实现。在用软件来实现数字滤波器中，可以根据滤波对象的实际情况，设置所需要的滤波器参数，改变滤波器的性能，达到所需要的滤波效果。在台站日常地震资料处理中，由于各种环境因素的影响，可能使得震相图中出现区别于地震信息的其他干扰信号，如车辆、湖面水波动、爆破等等［1］。因此，对地震资料的滤波处理是一个十分重要的环节，有助于提高震相分析的准确性。要实现对地震资料的滤波处理，必须设计相应的滤波器。对于数字滤波器的软件实现，可以借助高级语言（如Fortran、MATLAB语言）来实现，但是这均需要比较好的计算机编程能力，对于台站年龄较大的计算机基础或信号处理基础较薄弱的同志，很难单独实现，笔者直接使用MATLAB中自带的信号处理工具箱SPTool工具较方便的实现了数字滤波器的设计，可以省略了复杂的滤波器程序编写过程［2］，并完成了地震信号的滤波处理。

1 数字滤波器的基本原理

从数学意义来说，数字滤波器是一个离散的时间系统，输入x（n）是一个时间序列，输出y（n）也是一个时间序列，如数字滤波器脉冲响应为h（n），时间域内存在下列关系：y（n）=x（n）⊗y（n），频率域内输入输出关系为 Y（jw）=X（jw）H（jw）。同时滤波器是一个选频装置。理想的滤波器应该能无失真地传输有用信号，而又能完全的抑制掉无用信号。有用信号和无用信号往往占有不同的频带。信号能够通过滤波器的频带成为通带（passband），信号被抑制的频带称为阻带（stopband）。理想滤波器的频率特征可以写为：［3］

理想滤波器是物理上不可实现的系统，实际滤波器的频率特性只能够 “逼近”理想滤波器。按频率域特性来讲，数字滤波器可以分为低通、高通、带通、带阻等。由于实际的滤波器只能 “逼近”理想滤波器，所以，实际中滤波器的幅频响应不能如（1）、（2）式那样理想，通带内部不是完全平直，而是呈波纹变化；在阻带内，其幅频特性也不为零，而是衰减至某个值；在通带和阻带之间还存在一个过渡带，也并不是突然下降。因此设计滤波器的技术指标就包括了通带波纹Rp、阻带衰减Rs、阻带边界频率ws、通带边界频率wp。滤波器通带波纹Rp为相对于频率响应最大点的下降，因此下降越少越好，越少说明通带越平直，滤波效果越好，阻带衰减Rs也是相对频率响应最大点的下降，下降越多说明信号在阻带里受到越大的抑制，滤波效果越好。因此在实际应用中只能根据具体实际需要来设计滤波器。一般来说数字滤波器分为无限冲激和有限冲激（脉冲）响应数字滤波器两类。两种滤波器各有优缺点，在此不再详述，利用MATLAB的SPTool均能根据滤波器的性能要求实现以上两种滤波器设计，文中后面主要以IIR滤波器设计为例。

2 SPTool实现滤波器设计

一般情况下，如果不使用信号处理工具箱SPTool来设计滤波器，严格的数字滤波器设计，以IIR滤波器为例，在MATLAB中一般IIR数字滤波器设计的步骤如下：（1）将给定的数字滤波器的性能指标进行转换，转换后的频率指标作为模拟滤波器原型设计指标；（2）估计模拟滤波器的最小阶次和边界频率；（3）设计模拟低通滤波器原型；（4）由模拟原型经频率变换得到实际的模拟（低通、高通、带通、带阻）滤波器；（5）将模拟滤波器转换为IIR数字滤波器，可利用脉冲响应不变法或双线性变换法，以上步骤对于熟悉信号处理以及具备计算机高级语言编程能力的人来说，能够单独完成全部过程，但是对于不懂计算机语言编程的人来说，以上过程十分的繁琐复杂，倘若利用MATLAB自带的仿真工具SPTool，就会使得以上繁琐的滤波器设计过程变得十分简单。SPTool工具是MATLAB自带的一个交互式信号处理工具，专门用于完成常用的数字信号处理任务。通过这个工具，只要按SPTool界面要求，基层台站工作人员均能设计好实际工作中需要的滤波器，并输入相应的地震数据文件，就可以完成十分复杂的地震数字信号滤波处理任务，而不需要用户对数字滤波器的设计原理非常熟悉。下面实例说明笔者使用MATLAB中的交互式的信号处理工具SPTool滤波器设计。首先，我们设计一个阻带边界频率为1Hz，通带边界频率为2Hz，通带波纹为1dB，阻带衰减为30dB，采用Butterworth滤波器。点击进入SPTool中（Filter Design）功能模块，如图1所示。

图1 滤波器设计界面Fig.1 The printscreen of filter designing

可以指定滤波器的性能参数（如采样频率、通带频率范围等），也可以选择设计函数即设计方法，如Equiripple FIR，LeastSquares FIR，切比歇夫IIR，椭圆IIR等。滤波器设计完之后，另一个重要的事情是进行滤波器性能分析，SPTool提供了进行滤波器的频域和时域分析，（Filter View）是一个基于图形用户界面的分析工具，可以显示特定滤波器6个不同的特征子图，包括幅值响应、相位响应、群延迟、零极点和脉冲响应、阶跃响应。其中相位响应是滤波器频率响应的角度分量，零极点图用于显示滤波器传递函数的零极点与单位圆的位置。下图2是（Filter View）给出的幅频响应与相位响应。点选图中 group delay、phase、zeros and poles、impulse response、 step response等可以得到群延迟、零极点和脉冲响应、阶跃响应等其他图形。

图2 滤波器特征分析界面Fig.2 The printscreen of filter character analyzing

3 地震实例滤波

台站每日观测得到的震相数据中，除了由于地震自身所引发的之外，还包含了地脉动、海浪干扰等低频干扰，给地震波震相标注带来了麻烦。下面的实例，就给出了如何利用上面SPTool设计的高通滤波器进行低频干扰滤除。数据来源为五大连池地震火山监测站2007年11月9日16时57分一个近地方震，利用SPTool工具中Import功能导入地震数据，导入数据必须先将地震数据格式转换为MATLAB中的MAT文件。SPTool中（Signal Browser）会显示所导入的地震波形，如图3所示。然后在SPTool中点选APPLY利用前述滤波器进行滤波，滤波时可以选择零相位差滤波，消除IIR滤波器的时域相位延迟，如图4所示。滤波之后，SPTool中（Signal Browser）显示出滤波结果，如图5所示。通过滤波结果可以看到，地震波中低频干扰得到了很好的消除，突出了明显的近地方震震相，有助于台站地震分析人员明确其震相。