地震matlab频域分析,基于matlab的地震数据的分析.doc

基于matlab的地震数据的分析.doc

基于MATLAB的地震数据的分析孙玉柱冯光房桂梅摘要地震波原始数据中存在的干扰信号，会影响震相分析的准确性。为了滤除干扰信号，对地震波原始信号进行了频谱分析，给出了一种基于MATLAB的FIR数字滤波器的优化设计方案，将其用于地震波数据的分析中，并进行了仿真分析。仿真结果表明，FIR数字滤波器对地震波原始信号进行滤波处理后，提高了震相分析的准确性，得到了理想的效果，达到了预期的目的。4基于MATLAB的地震数据的分析下面以汕头台受到距台站300M处的汽车干扰的波形记录图的数字地震记录资料为例，进行基于MATLAB的地震数据的分析。该地震台采用的是频带范围为00520HZ的FBS3A型宽频带数字地震记录仪和EDASC24型数据采集器，系统的采样频率为50HZ。将该地震波原始数据保存在MATLAB的WORK文件夹中，并对其进行地频谱分析，如图1所示。从图1中可以很清楚的看出地震波原始记录中地震信号的优势频率为025HZ，频段范围在01HZ之间；干扰的优势频率为125HZ，频段范围在1015HZ之间。050100150100001000前前/S前前/COUNTS前前前前前0510152025012X106前前/HZ前前前前前前前前前05101520251001051010前前/HZ前前前前前前前前前前前前图1地震波原始数据的时域波形图和频谱分析图FIGURE1THEORIGINALTIMEDOMAINANDSPECTRUMANALYSISOFTHEEARTHQUAKEDATA为了滤除干扰信号，最大限度的保留其中的有用信号，选取的FIR滤波器为带阻滤波器，其参数选定为通带上截止频率FP17，阻带下截止频率FS171，阻带上截止频率FS1189，通带下截止频率FP219，通带波纹峰值DP001，阻带波纹峰值DS001。则设计的带阻FIR滤波器的频率特性曲线如图2所示。05101520256420X104FREQUENCYHZPHASEDEGREES051015202515010050050FREQUENCYHZMAGNITUDEDB图2带阻FIR滤波器的频率特性曲线FIGURE2THEFREQUENCYCHARACTERISTICOFBANDSTOPFIRFILTER由图2可以看出设计的FIR滤波器的特性满足设计要求。利用该滤波器对地震波原始记录信号进行滤波处理，滤波前后的时域波形图如图3所示，图4为滤波前后频域图的对比。05010015010005000500前前/S前前/COUNTS前前前前前前前前0501001506004002000200前前/S前前/COUNTS前前前前前前前前图3滤波前后地震波时域波形图的比较FIGURE3THECOMPARISONOFTIMEDOMAINBETWEENBEFOREANDAFTERFILTERPROCESSING0510152025051015X105前前/HZ前前前前前前前前前前前前0510152025051015X105前前/HZ前前前前前前前前前前前前图4滤波前后频域波形图的比较FIGURE4THECOMPARISONOFFREQUENCYDOMAINBETWEENBEFOREANDAFTERFILTERPROCESSING由图3和图4可以看出，地震波原始记录信号经FIR滤波器处理后，干扰波信号被滤除了，而地震波信号很好的显示了出来。由于FIR滤波器具有相位滞后的特点，所以滤波后的信号和原信号相比有了一定程度的相位延迟。为了不产生相位延迟，使滤波后的波形图能更好的与原始信号进行比较，可以对程序2(见附件)进行修改第38行后加MM1/2/FS计算相位延迟第50行后加TTM计算相位延迟第53行后加XLIM0,MAXT然后再对地震波原始记录信号进行处理，则滤波前后地震波时域波形图如图5所示，该图形非常明显的反映出了滤波前后地震波的特征，得到了理想的效果，达到了预期设计的目的。05010015010005000500前前/S前前/COUNTS前前前前前前前前0204060801001201406004002000200前前/S前前/COUNTS前前前前前前前前图5滤波前后时域波形图的比较FIGURE5THECOMPARISONOFTIMEDOMAINBETWEENBEFOREANDAFTERFILTERPROCESSING5结束语地震数据中包含了很多干扰信号，直接影响分析的准确性。干扰信号的处理虽然已经有很多方法，但本文采用目前非常流行的MATLAB软件，使用数字信号处理中的快速FOURIER变换和最优滤波器的设计方法，对采集到的地震数据进行频谱分析和滤波处理，去除干扰并最大限度地保留有用的频率成份，并使信号无失真，提高了震相分析的准确度。该方法可用于结构地震动力分析、地震台等领域，对地震的观测、分析、预报和研究有着极其重要的作用。附件程序1对地震记录数据进行频谱分析的实现程序LOADITXT调入原始数据文件(加载地震波数据)XTI得到原始信号序列FS50采样频率为50HZDT1/FS采样间隔(单位为S)NLENGTHXT原始信号序列长度XFFFTXT对原始波形数据进行快速FOURIER变换N0N1TNDT得到时间序列FN/NDT得到频率序列SUBPLOT3,1,1时域坐标方框图PLOTT,XT画出时域中的原始波形图XLABEL 时间/S X轴标示YLABEL 振幅/COUNTS Y轴标示TITLE 时域波形图 加注标题GRIDONSUBPLOT3,1,2频域坐标方框图PLOTF,ABSXF画出频域中的FFT波形图XLABEL 频率/HZ X轴标示YLABEL 频谱密度 Y轴标示TITLE 频域波形图 加注标题XLIM0FS/2频域之画出采样频率的一半GRIDONSUBPLOT3,1,3频域坐标方框图SEMILOGYF,ABSXF画出频域中的FFT波形图，Y轴为对数XLABEL 频率/HZ X轴标示YLABEL 频谱密度取对数 Y轴标示TITLE 频域波形图 加注标题XLIM0FS/2频域之画出采样频率的一半GRIDON程序2用FIR数字滤波器实现的滤波程序LOADITXT调入原始数据文件(加载地震波数据)XTI得到原始信号序列FS50采样频率为50HZDT1/FS采样间隔(单位为S)NLENGTHXT原始信号序列长度XFFFTXT对原始波形数据进行快速FOURIER变换N0N1TNDT得到时间序列FN/NDT得到频率序列FIGURE1SUBPLOT2,1,1时域坐标方框图PLOTT,XT画出时域中的原始波形图XLABEL 时间/S X轴标示YLABEL 振幅/COUNTS Y轴标示TITLE 时域波形图 加注标题GRIDONSUBPLOT2,1,2频域坐标方框图PLOTF,ABSXF画出频域中的FFT波形图XLABEL 频率/HZ X轴标示YLABEL 频谱密度 Y轴标示TITLE 频域波形图 加注标题XLIM0FS/2频域之画出采样频率的一半GRIDONFP17带阻滤波器通带上截止频率FS171带阻滤波器阻带下截止频率FS2189带阻滤波器阻带上截止频率FP219带阻滤波器通带下截止频率DP001带阻滤波器通带波纹峰值DS001带阻滤波器阻带波纹峰值W