基于MATLAB的地震数据的分析

基于MATLAB的地震数据的分析孙玉柱冯光房桂梅摘要：地震波原始数据中存在的干扰信号，会影响震相分析的准确性。

为了滤除干扰信号，对地震波原始信号进行了频谱分析，给出了一种基于MATLAB的FIR数字滤波器的优化设计方案，将其用于地震波数据的分析中，并进行了仿真分析。

仿真结果表明，FIR数字滤波器对地震波原始信号进行滤波处理后，提高了震相分析的准确性，得到了理想的效果，达到了预期的目的。

关键词：MATLAB；FIR数字滤波器；优化；滤波the Analysis of Earthquake Data Based on MATLAB SUN Yuzhu,FENG Guang,FANG Guimei Abstract: The interference that existed in the earthquake data will affect the accuracy of the seismic phase analysis. In order to filter the disturbance signal, this paper carries out spectrum analysis of the earthquake data, proposes an optimum design method for FIR digital filter based on MATLAB and applies it to the analysis of earthquake data. After the filter of the noise jamming, the true information of the earthquake wave is clearly reflected. The simulation results manifest that it canimprove the accuracy of seismic phase analysis and arrive at the purpose desired.Key words: MATLAB；FIR digital filter；optimization；filter1 引言地震带给人类的损失是巨大的，汶川大地震依旧在我们的记忆深处清晰存在。

大地震的每次不约而至，都对国家和人民造成了巨大的损失。

地震预测是世界性难题，全世界的地震科学家不断在探索，尽最大努力减少其破坏性。

地震台站提供的地震观测资料的可靠性和准确性，是地震学家进行地震预测的基础[1]。

但是地震波信号变化的不平稳性和复杂性给地震的分析和预测带来了很大困难，并且在地震波的原始记录中往往还掺杂着来自外界的各种干扰，如仪器、环境噪声、爆破、采矿、火车的震动等，这些都给地震波的分析带来了严重影响，甚至导致分析结果的错误。

为保证地震分析的准确性，可先对原始记录进行频谱分析，选择性能优良的滤波器对其进行优化处理，把干扰信号尽量滤除，然后再对处理后的地震波数据进行分析处理，则会得到良好的效果。

地震记录的数字化使得利用计算机对地震信号进行分析处理得以实现。

在数字信号的分析处理中，Fourier变换和数字滤波器的应用极为普遍，语音、雷达、地震、图像、机械振动、地质勘探等众多领域都广泛采用数字滤波器。

MATLAB是一种集数值分析、矩阵运算、信号处理和图像显示于一体、功能极其强大的高性能软件，其工具箱中包含了各种经典和现代数字信号处理技术，很容易实现Fourier 变换和各种数字滤波器的设计，在地震数据的分析处理中起着重要作用。

本文首先介绍了数字滤波器，给出了基于MATLAB 的FIR 数字滤波器的一种优化设计方案，最后将其用于地震波数据的处理中，并进行了仿真分析。

2快速Fourier 变换（FFT ）及频谱分析在地震波的原始记录数据中往往夹杂不同频率范围的噪声干扰信号，为了显示出地震波数据中的优势频率和干扰频率，保证地震分析的准确性，应首先采用频谱分析，再针对干扰波的频率范围，设计合适的滤波器参数。

在对有限长信号序列进行频谱分析时，离散Fourier 变换（DFT ）应用非常广泛，它可以很好地反映序列的频谱特性[2]。

设()x n 是一个长度为N 的有限长序列，则()x n 的N 点离散傅里叶变换为[3]()120()(), 0,1,...,1N j N kn n X k x n e k N π--===-∑ （1）DFT 是信号分析与处理中的一种重要变换，但是当N 较大时，其计算量太大。

快速Fourier 变换（FFT ）是减少 DFT 运算次数的一种快速算法，通过在时域将序列逐次分解为一组子序列，然后利用子序列的DFT 来实现整个序列的DFT ，从而减少离散Fourier 变换的运算量，提高了计算效率。

在MATLAB中可调用函数()X=fft x,N来进行快速傅里叶变换[4]。

f t其中，x为时域内的输入信号序列，N为序列长度，f X为频率域的t输出信号，即x的频谱特征。

t3 FIR数字滤波器的优化设计3.1 数字滤波器的选择在地震分析中必须要先对原始信号做滤波处理，滤波的目的是为了去除噪声，使原始信号通过滤波器后能够清晰地显示出优势频率，为更好的分析地震信号（比如震相等）做准备。

滤波器包括模拟滤波器和数字滤波器，模拟滤波器又可分为无源滤波器（主要由R、C和L构成）和有源滤波器（主要由集成运放和R、C元件构成），数字滤波器可用计算机软件或大规模集成数字硬件实现。

模拟滤波器存在电压漂移、温度漂移和噪声等问题，而数字滤波器不存在这些问题，可以达到很高的稳定度和精度。

根据实现的网络结构不同，数字滤波器可分为无限脉冲响应（IIR）滤波器和有限脉冲响应（FIR）滤波器。

1. IIR数字滤波器IIR数字滤波器最大的优点是可取得非常好的通带和阻带衰减，还可得到准确的通带与阻带的边缘频率，而且滤波时需要的计算量小；缺点是滤波器的响应灵活性较差，不具有线性相位（若得到线性相位则需要用全通系统进行相位补偿）且存在稳定性问题，而这些在实际应用（如图像信号处理、地震信号处理、数据通信等领域）中又是非常重要的。

2. FIR 数字滤波器FIR 数字滤波器又称为卷积滤波器，通常采用迭代算法来满足设定的技术指标，由于FIR 系统是全零点系统，其单位抽样响应h(n)为有限长，因此容易实现某种对称性，从而获得在设计任意幅频特性的同时保证严格的线性相位特性；由于采用非递归结构，则不存在稳定性问题，运算误差也较少，并且能容易实现多通带、多阻带的滤波器的设计；又因为FIR 数字滤波器在通带内具有恒定的幅频特性和线性相位特性，从而使信号通过时不失真，有时还可实现零相位滤波。

考虑到地震波数据的特点，本文选用FIR 数字滤波器。

3.2 FIR 数字滤波器的MATLAB 优化设计FIR 数字滤波器的传递函数为[5]1()()()()N n n Y z H z h n z X z --===∑ （2） 其中，()h n 是FIR 滤波器的单位脉冲响应，其长度为N ，非零区间为[0,1]N -。

（1）式表明()H z 是1z -的(1)N -次多项式，它在z 平面上有(1)N -个零点。

原点0z =是(1)N -阶重极点。

因此，()H z 肯定是稳定的。

由FIR 滤波器的传递函数，可得到其频率响应表达式10()()N j j n n H e h n e ωω--==∑ （3） 如果()h n 是实序列，并且满足下列中心对称条件[6]()(1)h n h N n =--或()(1)h n h N n =---（4）则FIR 滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。

FIR 滤波器虽然具有相位滞后的缺点，但是其相位滞后和群延迟在整个频带上是相等且不变的。

一个N 阶的线性相位FIR 滤波器群延迟为2N ，即滤波后的信号简单地延长常数个时间步长，这一特性使通带频率内信号通过滤波器后仍然保持原有波形形状而无相位失真。

FIR 数字滤波器的设计方法很多，本文采用目前最优秀的设计方法——切比雪夫逼近法。

切比雪夫逼近法是一种等波纹逼近法，它使误差在整个频带均匀分布，对同样的技术指标，这种逼近法需要的滤波器阶数低，而对同样的滤波器阶数，这种逼近法的最大误差最小。

设希望设计的滤波器幅度特性为()d H ω，实际设计的滤波器幅度特性为()g H ω，误差加权系数为()W ω，则切比雪夫最佳一致逼近准则是使加权函数[7]()()()()d g E W H H ωωωω⎡⎤=-⎣⎦ （5） 的最大值达到最小，即满足min max ()A E ωω∈⎡⎤⎣⎦（6） 切比雪夫逼近理论可以解决()E ω的存在性、唯一性及如何构造等一系列问题，克服了通带和阻带的边缘不易精确确定和Gibbs 现象等多种缺点，使设计出的滤波器具有明显优点。

因为在一定意义上对()d H ω做最佳逼近，可获取较好的通带和阻带性能，并能准确地指定通带和阻带的边缘频率，此时又具有了IIR 滤波器的优点，所以这是一种很好的设计方法。

MATLAB 提供了大量设计FIR 数字滤波器的函数，下面针对切比雪夫逼近法，介绍其设计步骤和方法。

（1）根据得出的地震波原始记录的频谱分析波形，确定FIR数字滤波器的技术指标；（2）用函数[]()M,F,A,W=remezord f,a,dev,Fs[8]估算设计的等波纹逼00近法的参数：最低滤波器阶数M、频率向量F、幅度向量0A和加权向量W。

其中，Fs为采样频率；f可以是模拟频率或归一化频率，但必须以0开始，以Fs2（用归一化频率为1）结束，而且省略了0和Fs2两个频点；dev为各逼近段允许的幅度响应偏差（波纹振幅）；a 为滤波器在各个频段上的幅度值，一般对通带取值为1，对阻带取值为0。

（3）用函数()h=remez M,F,A,W[8]完成等纹波FIR滤波器的设计，00并用函数()y=filtfilt h,1,x完成对输入信号x的滤波，得到y=filter h,1,x或()输出信号y。

4 基于MATLAB的地震数据的分析下面以汕头台受到距台站300m处的汽车干扰的波形记录图的数字地震记录资料为例，进行基于MATLAB的地震数据的分析。

该地震台采用的是频带范围为0.05~20Hz 的FBS-3A型宽频带数字地震记录仪和EDAS-C24型数据采集器，系统的采样频率为50Hz。

将该地震波原始数据保存在MATLAB的work文件夹中，并对其进行地频谱分析，如图1所示。

从图1中可以很清楚的看出地震波原始记录中地震信号的优势频率为0.25Hz，频段范围在0~1Hz之间；干扰的优势频率为12.5Hz，频段范围在10~15Hz之间。