基于MATLAB的地震数据的分析.docx
《基于MATLAB的地震数据的分析.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的地震数据的分析.docx（17页珍藏版）》请在冰点文库上搜索。
基于MATLAB的地震数据的分析
基于MATLA啲地震数据的分析
孙玉柱冯光房桂梅摘要：
地震波原始数据中存在的干扰信号，会影响震相分析的准确性。
为了滤除干扰信号，对地震波原始信号进行了频谱分析，给出了一种基于MATLAB勺FIR数字滤波器的优化设计方案，将其用于地震波数据的分析中，并进行了仿真分析。
仿真结果表明，FIR数字滤波器对地震波原始信号进行滤波处理后，提高了震相分析的准确性，得到了理想的效果，达到了预期的目的。
关键词：
MATLABFIR数字滤波器；优化；滤波
theAnalysisofEarthquakeDataBasedonMATLAB
SUNYuzhu,FENGGuang,FANGGuimei
Abstract:
Theinterferencethatexistedintheearthquakedatawillaffecttheaccuracyoftheseismicphaseanalysis.Inordertofilterthedisturbancesignal,thispapercarriesoutspectrumanalysisoftheearthquakedata,proposesanoptimumdesignmethodforFIRdigitalfilterbasedonMATLABandappliesittotheanalysisofearthquakedata.Afterthefilterofthenoisejamming,thetrueinformationoftheearthquakewaveisclearlyreflected.Thesimulationresultsmanifestthatitcan
improvetheaccuracyofseismicphaseanalysisandarriveatthepurposedesired.
Keywords:
MATLA；BFIRdigitalfilter；optimization；filter
1引言
地震带给人类的损失是巨大的，汶川大地震依旧在我们的记忆深处清晰存在。
大地震的每次不约而至，都对国家和人民造成了巨大的损失。
地震预测是世界性难题，全世界的地震科学家不断在探索，尽最大努力减少其破坏性。
地震台站提供的地震观测资料的可靠性和准确性，是地震学家进行地震预测的基础[1]。
但是地震波信号变化的不平稳性和复杂性给地震的分析和预测带来了很大困难，并且在地震波的原始记录中往往还掺杂着来自外界的各种干扰，如仪器、环境噪声、爆破、采矿、火车的震动等，这些都给地震波的分析带来了严重影响，甚至导致分析结果的错误。
为保证地震分析的准确性，可先对原始记录进行频谱分析，选择性能优良的滤波器对其进行优化处理，把干扰信号尽量滤除，然后再对处理后的地震波数据进行分析处理，则会得到良好的效果。
地震记录的数字化使得利用计算机对地震信号进行分析处理得以实现。
在数字信号的分析处理中，Fourier变换和数字滤波器的应用极为普遍，语音、雷达、地震、图像、机械振动、地质勘探等众多领域都广泛采用数字滤波器。
MATLA是一种集数值分析、矩阵运算、信号处理和图像显示于一体、功能极其强大的高性能软件，其工具箱中包含了各种经典和现代数字信号处理技术，很容易实现Fourier变换和各种数字滤波器的设计，在地震数据的分析处理中起着重要作用。
本文首先介绍了数字滤波器，给出了基于MATLAB勺FIR数字滤波器的一种优化设计方案，最后将其用于地震波数据的处理中，并进行了仿真分析。
2快速Fourier变换（FFT）及频谱分析
在地震波的原始记录数据中往往夹杂不同频率范围的噪声干扰信号，为了显示出地震波数据中的优势频率和干扰频率，保证地震分
析的准确性，应首先采用频谱分析，再针对干扰波的频率范围，设计合适的滤波器参数。
在对有限长信号序列进行频谱分析时，离散Fourier变换（DFT应用非常广泛，它可以很好地反映序列的频谱特性[2]。
设x（n）是一个长度为N的有限长序列，则x（n）的N点离散傅里叶变换为⑶
N-1
X（k）八x（n）/2二Nkn,k=0,1，…，N-1
（1）
n=0
DFT是信号分析与处理中的一种重要变换，但是当N较大时，其计算量太大。
快速Fourier变换（FFT）是减少DFT运算次数的一种快速算法，通过在时域将序列逐次分解为一组子序列，然后利用子序列的DFT来实现整个序列的DFT从而减少离散Fourier变换的运算量，提高了计算效率。
在MATLAB^可调用函数Xf=fftxt,N来进行快速傅里叶变换⑷。
其中，Xt为时域内的输入信号序列，N为序列长度，Xf为频率域的输出信号，即Xt的频谱特征。
3FIR数字滤波器的优化设计
3.1数字滤波器的选择
在地震分析中必须要先对原始信号做滤波处理，滤波的目的是为了去除噪声，使原始信号通过滤波器后能够清晰地显示出优势频率，为更好的分析地震信号（比如震相等）做准备。
滤波器包括模拟滤波器和数字滤波器，模拟滤波器又可分为无源滤波器（主要由R、C和
L构成）和有源滤波器（主要由集成运放和R、C元件构成），数字滤波器可用计算机软件或大规模集成数字硬件实现。
模拟滤波器存在电压漂移、温度漂移和噪声等问题，而数字滤波器不存在这些问题，可以达到很高的稳定度和精度。
根据实现的网络结构不同，数字滤波器可分为无限脉冲响应（IIR）滤波器和有限脉冲响应（FIR）滤波器。
1.IIR数字滤波器
IIR数字滤波器最大的优点是可取得非常好的通带和阻带衰减，还可得到准确的通带与阻带的边缘频率，而且滤波时需要的计算量小；缺点是滤波器的响应灵活性较差，不具有线性相位（若得到线性相位则需要用全通系统进行相位补偿）且存在稳定性问题，而这些在实际应用（如图像信号处理、地震信号处理、数据通信等领域）中又是非常重要的。
2.FIR数字滤波器
FIR数字滤波器又称为卷积滤波器，通常采用迭代算法来满足设
定的技术指标，由于FIR系统是全零点系统，其单位抽样响应h（n）为有限长，因此容易实现某种对称性，从而获得在设计任意幅频特性的同时保证严格的线性相位特性；由于采用非递归结构，则不存在稳定性问题，运算误差也较少，并且能容易实现多通带、多阻带的滤波器的设计；又因为FIR数字滤波器在通带内具有恒定的幅频特性和线性相位特性，从而使信号通过时不失真，有时还可实现零相位滤波。
考虑到地震波数据的特点，本文选用FIR数字滤波器。
3.2FIR数字滤波器的MATLA优化设计
其中，h（n）是FIR滤波器的单位脉冲响应，其长度为N，非零区间为
[0,N-1]。
（1）式表明H（z）是z」的（N-1）次多项式，它在z平面上有（N-1）个零点。
原点z=0是（N-1）阶重极点。
因此，H（z）肯定是稳定的。
由FIR滤波器的传递函数，可得到其频率响应表达式
N-1
H©）八h（n）e_n（3）
n=0
如果h（n）是实序列，并且满足下列中心对称条件⑹
h（n）=h（N-n-1）或h（n）=-h（N-nT）（4）
则FIR滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。
FIR滤波器虽然具有相位滞后的缺点，但是其相位滞后和群延迟在整个频带上是相等且不变的。
一个N阶的线性相位FIR滤波器群延迟为N2，即滤波后的信号简单地延长常数个时间步长，这一特性使通带频率内信号通过滤波器后仍然保持原有波形形状而无相位失真。
FIR数字滤波器的设计方法很多，本文采用目前最优秀的设计方法一一切比雪夫逼近法。
切比雪夫逼近法是一种等波纹逼近法，它使
误差在整个频带均匀分布，对同样的技术指标，这种逼近法需要的滤波器阶数低，而对同样的滤波器阶数，这种逼近法的最大误差最小。
设希望设计的滤波器幅度特性为Hd（），实际设计的滤波器幅度特性为Hg（J，误差加权系数为WC），则切比雪夫最佳一致逼近准则是使加权函数⑺
E（）=W（）||HdCJ-Hg（）（5）
的最大值达到最小，即满足
minmaxE（时）I（6）
切比雪夫逼近理论可以解决EC）的存在性、唯一性及如何构造等一系列问题，克服了通带和阻带的边缘不易精确确定和Gibbs现象等多种缺点，使设计出的滤波器具有明显优点。
因为在一定意义上对HdC）做最佳逼近，可获取较好的通带和阻带性能，并能准确地指定通带和阻带的边缘频率，此时又具有了IIR滤波器的优点，所以这是一种很好的设计方法。
MATLAB提供了大量设计FIR数字滤波器的函数，下面针对切比雪夫逼近法，介绍其设计步骤和方法。
（1）根据得出的地震波原始记录的频谱分析波形，确定FIR数字滤波器的技术指标;
（2）用函数M,F0,A0,Wl=remezordf,a,dev,Fs问估算设计的等波纹逼近法的参数：
最低滤波器阶数M、频率向量Fo、幅度向量Ao和加权向量W。
其中，Fs为采样频率；f可以是模拟频率或归一化频率，但必须以0开始，以Fs2（用归一化频率为1）结束，而且省略了0和Fs2两个频点；dev为各逼近段允许的幅度响应偏差（波纹振幅）；a为滤波器在各个频段上的幅度值，一般对通带取值为1,对阻带取值为0。
（3）用函数h=remez（M,F0,A0,W）问完成等纹波FIR滤波器的设计，并用函数y=filterh,1,x或y=filtfilth,1,x完成对输入信号x的滤波，得到输出信号y。
4基于MATLA的地震数据的分析
下面以汕头台受到距台站300m处的汽车干扰的波形记录图的数字地震记录资料为例，进行基于MATLAB勺地震数据的分析。
该地震台采用的是频带范围为0.05~20Hz的FBS-3A型宽频带数字地震记录仪和EDAS-C2鰹数据采集器，系统的采样频率为50Hz
将该地震波原始数据保存在MATLAB勺work文件夹中，并对其进行地频谱分析，如图1所示。
从图1中可以很清楚的看出地震波原始记录中地震信号的优势频率为0.25Hz，频段范围在0~1Hz之间；干扰的优势频率为12.5Hz,频段范围在10~15Hz之间。
时域波形图
图1地震波原始数据的时域波形图和频谱分析图
Figure1theoriginaltime-domainandspectrumanalysisofthe
earthquakedata
为了滤除干扰信号，最大限度的保留其中的有用信号，选取的
FIR滤波器为带阻滤波器，其参数选定为：
通带上截止频率Fp1=7,
阻带下截止频率Fs1=7.1，阻带上截止频率Fs1=18.9，通带下截止频率Fp2=19,通带波纹峰值dp=0.01，阻带波纹峰值ds=0.01。
则设计的带阻FIR滤波器的频率特性曲线如图2所示。
图2带阻FIR滤波器的频率特性曲线
Figure2thefrequencycharacteristicofband-stopFIRfilter
由图2可以看出设计的FIR滤波器的特性满足设计要求。
利用该滤波器对地震波原始记录信号进行滤波处理，滤波前后的时域波形图
如图3所示，图4为滤波前后频域图的对比。
图3滤波前后地震波时域波形图的比较
Figure3thecomparisonoftime-domainbetweenbeforeandafterfilter
processing
图4滤波前后频域波形图的比较
Figure4thecomparisonoffrequency-domainbetweenbeforeandafter
filterprocessing
由图3和图4可以看出，地震波原始记录信号经FIR滤波器处理后，干扰波信号被滤除了，而地震波信号很好的显示了出来。
由于FIR滤波器具有相位滞后的特点，所以滤波后的信号和原
信号相比有了一定程度的相位延迟。
为了不产生相位延迟，使滤波后
行修改:
第38行后加m=（M-1）/2/Fs;%计算相位延迟
第50行后加t=t-m;%计算相位延迟
第53行后加xlim（[0,max（t）]）;
然后再对地震波原始记录信号进行处理，则滤波前后地震波时域
波形图如图5所示，该图形非常明显的反映出了滤波前后地震波的特
征，得到了理想的效果，达到了预期设计的目的
图5滤波前后时域波形图的比较
Figure5thecomparisonoftime-domainbetweenbeforeandafterfilter
processing
5结束语
地震数据中包含了很多干扰信号，直接影响分析的准确性。
干扰信号的处理虽然已经有很多方法，但本文采用目前非常流行的MATLAB软件，使用数字信号处理中的快速Fourier变换和最优滤波器的设计方法，对采集到的地震数据进行频谱分析和滤波处理，去除干扰并最
大限度地保留有用的频率成份，并使信号无失真，提高了震相分析的准确度。
该方法可用于结构地震动力分析、地震台等领域，对地震的观测、分析、预报和研究有着极其重要的作用。
参考文献
[1]宋建锁.滤波在地震分析中的应用[J].防灾技术高等专科学校学报,2006,8
（1）:
75-79
[2]李敬,甘延锋,黄友明.数字地震记录中干扰波的排除[J].防灾技术高等专科学校学报,2004,6（3）:
20-25
[3]丁玉美,高西金.数字信号处理[M].西安:
西安电子科技大学出版社,2004
[4]万永革.数字信号处理的MATLAB实现[M].北京：
科学出版社,2007
[5]胡广书.数字信号处理[M].北京：
清华大学出版社,2003
[6]陈鲁刚,王锐锋.用MATLA实现滤波器对数字地震资料处理初探
[J].防灾技术高等专科学校学报,2006,8
（1）：
62-65
[7]陈后金.数字信号处理[M].北京：
高等教育出版社,2004
[8]唐向红,岳恒立,郑雪峰.MATLAB及在电子信息类课程中的应用[M].北京：
电子工业出版社,2006
附件
程序1对地震记录数据进行频谱分析的实现程序
loadi.txt%
调入原始数据文件（加载地震波数据）
Xt=i;
%
得到原始信号序列
Fs=50;
%
采样频率为50Hz
dt=1/Fs;
%
采样间隔（单位为s）
N=length（Xt）;
%
原始信号序列长度
Xf=fft（Xt）;%
对原始波形数据进行快速Fourier变换
n=0:
N-1;
t=n*dt;
%
得到时间序列
f=n/（N*dt）;
%
得到频率序列
subplot（3,1,1）;
%
时域坐标方框图
plot（t,Xt）;
%
画出时域中的原始波形图
xlabel（'时间/s'）;
%X
轴标示
ylabel（'振幅/counts'）;%Y
轴标示
title（'时域波形图'）;
%
加注标题
gridon
subplot（3,1,2）;
%
频域坐标方框图
plot（f,abs（Xf））;
%画出频域中的FFT波形图
xlabel（'频率/Hz'）;
%X
轴标示
ylabel（'频谱密度'）;
%Y
轴标示
加注标题
频域之画出采样频率的一半
%
频域坐标方框图
画出频域中的
FFT波形图，丫轴为对数
%X
轴标示
%Y
轴标示
%
加注标题
频域之画出采样频率的一半
title（'频域波形图'）;xlim（[0Fs/2]）;%gridonsubplot（3,1,3）;semilogy（f,abs（Xf））;%xlabel（'频率/Hz'）;ylabel（'频谱密度取对数'）;title（'频域波形图'）;xlim（[0Fs/2]）;%gridon
loadi.txt
%
Xt=i;
%
Fs=50;
%
dt=1/Fs;
%
N=length（Xt）;%Xf=fft（Xt）;%n=0:
N-1;
t=n*dt;%f=n/（N*dt）;%figure
（1）subplot（2,1,1）;%plot（t,Xt）;%
程序2用FIR数字滤波器实现的滤波程序
调入原始数据文件（加载地震波数据）
得到原始信号序列
采样频率为50Hz
采样间隔（单位为s）
原始信号序列长度
对原始波形数据进行快速Fourier变换
得到时间序列
得到频率序列
时域坐标方框图
xlabel（'
时间/s'）;
%X
轴标示
ylabel（'
振幅/counts'）;
%Y
轴标示
title（'
时域波形图'）;
%
加注标题
gridon
subplot（2,1,2）;
%
频域坐标方框图
plot（f,abs（Xf））;
%
画出频域中的FFT波形图
xlabel（'
频率/Hz'）;
%X
轴标示
ylabel（'
频谱密度'）;
%Y
轴标示
title（'
频域波形图'）;
%
加注标题
画出时域中的原始波形图
figure（3）;
subplot（2,1,1）;
plot（t,Xt）;%
画出时域中的原始波形图
xlabel（'时间/s'）;
%X
轴标示
ylabel（'振幅/counts'）;
%Y
轴标示
title（'滤波前时域波形图'）;
%
加注标题
gridon
z=filter（h,1,Xt）;%
对原始信号进行滤波
subplot（2,1,2）;
plot（t,z）;%
画出滤波后的时域波形图
xlabel（'时间/s'）;
%X
轴标示
ylabel（'振幅/counts'）;
%Y
轴标示
title（'滤波后时域波形图'）;
%
加注标题
gridon
Q=length（z）;%
求滤波后信号序列长度
figure（4）;
subplot（2,1,1）;
plot（f,abs（Xf））;%
画出频域中的原始FFT波形图
xlabel（'频率/Hz'）;
%X
轴标示
ylabel（'频谱密度'）;
%Y
轴标示
title（'滤波前频域波形图'）;
%
加注标题
xlim（[0Fs/2]）;%
频域直画出采样频率的一半
gridon
subplot（2,1,2）;
plot（[0:
Q-1]/（Q*dt）,abs（p））;%
画出频域中的滤波后信号的FFT
波形图
xlabel（'
频率/Hz'）;
%X
轴标示
ylabel（'
频谱密度'）;
%Y
轴标示
title（'
滤波后频域波形图'）;
%
加注标题
xlim（[0Fs/2]）;%
频域直画出采样频率的一半
p=fft（z）;
对滤波后的波形数据进行快速Fourier变换
gridon
升级会员 