MatLab在成都地震臺數(shù)字地震記錄中的應(yīng)用

何思源，劉華嬌，徐文海，鐘李彬，董立杰

(成都地震基準(zhǔn)臺 四川省地震局，四川 成都 611730)

MatLab在成都地震臺數(shù)字地震記錄中的應(yīng)用

何思源，劉華嬌，徐文海，鐘李彬，董立杰

(成都地震基準(zhǔn)臺 四川省地震局，四川 成都 611730)

成都地震基準(zhǔn)臺數(shù)字臺網(wǎng)中記錄到的小震級地震數(shù)據(jù)中有的包含較為嚴(yán)重的干擾，本文利用MatLab計算軟件對數(shù)字地震記錄進(jìn)行頻譜分析，設(shè)計FIR、IIR兩種類型濾波器對這些干擾數(shù)據(jù)進(jìn)行濾波，對兩種不同類型的濾波器的應(yīng)用結(jié)果進(jìn)行對比分析。通過濾波器處理后干擾得到有效消除，數(shù)據(jù)質(zhì)量得到顯著提高，分析認(rèn)為2～4 Hz的干擾信息是成都地震臺地震記錄質(zhì)量不佳的重要原因，這對于今后的地震數(shù)據(jù)分析工作具有一定的參考價值。

MatLab；數(shù)字地震記錄；干擾信息；濾波器

成都地震基準(zhǔn)臺目前使用了JCZ-1型和JCZ-1T型兩套甚寬頻帶地震計[1]。其中JCZ-1是分體裝置，由一個垂向和兩個水平向地震計構(gòu)成，頻帶寬為20 Hz-DC，JCZ-1T型是JCZ-1型的改進(jìn)型，是三分量一體機，地震計頻帶寬為50 Hz-DC。地震計數(shù)字化極大豐富了觀測的信息量，同時人們又面臨如何去處理和使用這樣規(guī)模龐大的數(shù)字化觀測資料的問題。地震計儀器在數(shù)據(jù)采集的時候，也會采集到很多干擾信號，這些干擾嚴(yán)重影響了觀測地震波數(shù)據(jù)的質(zhì)量，因此需要對這些干擾加以排除[2]。MatLab軟件是一套進(jìn)行科學(xué)計算的高性能軟件，興起于20世紀(jì)70年代，在傳統(tǒng)的C語言和FORTRAN計算語言的基礎(chǔ)之上，MatLab語言以更接近數(shù)學(xué)公式的表達(dá)方式，給用戶提供了最簡潔最直觀的數(shù)字計算開發(fā)環(huán)境[3]。本文介紹在MatLab軟件中利用快速傅里葉變換(FFT)方法實現(xiàn)有限序列長度的地震波形記錄頻譜分析，以及FIR、IIR兩種濾波器的設(shè)計原理和應(yīng)用于成都臺測震數(shù)據(jù)后的對比分析。

1 快速傅里葉變換(FFT)原理及頻譜分析MatLab實現(xiàn)

傅里葉變換在信號處理中具有十分重要的作用，但是離散時間的傅里葉變換具有很大的時間復(fù)雜度。非周期性連續(xù)時間信號x(t)的傅里葉變換可以表示為：

(1)

公式中計算出來的是信號x(t)的連續(xù)頻譜。但是，在實際的控制系統(tǒng)中能夠得到的是連續(xù)信號x(t)的離散采樣值x(nT)。因此需要利用離散信號x(nT)來計算信號x(t)的頻譜。有限長離散信號x(n)，n=0，1，…，N-1的DFT定義為：

(2)

這里，k=0,1,…,N-1?？梢钥闯?，DFT需要計算大約N2次乘法和N2次加法。當(dāng)N較大時，這個計算量是很大的。而快速傅里葉變換(FFT)減少了DFT的運算次數(shù),通過在時域?qū)⑿蛄兄鸫畏纸鉃橐唤M子序列，然后利用子序列的DFT來實現(xiàn)整個序列的DFT，從而減少了離散傅里葉變換的運算量，提高了計算效率。在MatLab中對輸入信號實現(xiàn)快速傅里葉變換的命令如下[4]：

Xf=fft(xt,N)

(3)

其中Xt是輸入信號序列，N為序列的長度，Xf是輸出的信號序列，也就是Xt的頻譜特征。圖1是成都臺JCZ-1T于2016年6月18日09時23分記錄到的Unimak Island Region, Alaska地區(qū)的地震，震級M4.7，震中距67°，震源深度10 km。通過對該地震事件進(jìn)行快速傅里葉(FFT)變換，得到了地震的優(yōu)勢頻率和干擾信號的優(yōu)勢頻率。由圖1中可以看出，經(jīng)過快速傅里葉(FFT)變換后得到的地震記錄的頻譜中，地震信號的優(yōu)勢頻率是在0-1 Hz的區(qū)間，干擾信號的優(yōu)勢頻率是在2～5 Hz的區(qū)間。在設(shè)計濾波器對干擾部分進(jìn)行“過濾”的同時，應(yīng)當(dāng)對于包含地震信息的頻率部分盡可能保護(hù)，這一點在設(shè)計濾波器時需要格外注意。

圖1 地震記錄的原始波形

圖2 地震記錄的頻率譜

2 FIR與IIR數(shù)字濾波器的設(shè)計

在地震分析中必須要先對原始信號進(jìn)行濾波處理，濾波的目的是為了去除噪聲，使原始信號通過濾波器后能夠清晰地顯示出優(yōu)勢頻率。數(shù)字濾波器的主要功能是對數(shù)字信號進(jìn)行處理時最常見的方法是保留數(shù)字信號中有用的頻率成分，去除信號中無用的頻率成分[4]。

2.1 FIR數(shù)字濾波器原理及設(shè)計

FIR濾波器的單位沖擊響應(yīng)h(n)是有限長(0≤n≤N-1)，其Z變換為：

(4)

在有限Z平面有(N-1)個零點，而它的(N-1)個極點均位于原點z=0處。FIR濾波器的系統(tǒng)差分方程為：

(5)

h(n)的頻率響應(yīng)H(ejw)可表示為：

(6)

信號通過FIR不失真條件是濾波器在通帶內(nèi)具有恒定的幅頻特性和線性相位特性。線性相位FIR濾波器的相位滯后和群延遲在整個頻帶上是相等且不變的，所以，F(xiàn)IR濾波器雖然有相位延遲，但是通過濾波器的波形仍然保持原有形狀且不會相位失真。FIR帶阻濾波器的設(shè)計思路大致如下：確定帶阻濾波器性能參數(shù)，如采樣頻率Fs，通帶上、下截止頻率flp,fhp,阻帶下、上截止頻率fls,fhs。通過“wlp=2*flp*pi/Fs”公式依次計算阻帶、通帶上下限歸一化角頻率wls,whs,whp,程序：wlp=2*flp*pi/Fs;whp=2*fhp*pi/Fs;wls=2*fls*pi/Fs;whs=2*fhs*pi/Fs(歸一化角頻率)。Fir1是MatLab信號處理工具箱中采用經(jīng)典窗函數(shù)法設(shè)計線性FIR數(shù)字濾波器的函數(shù)，調(diào)用格式：b=fir1(N,wn,'type',window)，其中b為濾波器系數(shù)；N為濾波器階數(shù)；wn為濾波器截止頻率；‘type’為設(shè)計濾波器的類型，如帶阻，type=stop，高通則type=high；window為窗函數(shù)類型。計算濾波器截止頻率wn，wn=[(wlp+wls)/(2*pi),(whp+whs)/(2*pi)]。計算濾波器階數(shù)N，程序如下:k1=wls-wlp;k2=whp-whs;kw=min(k1,k2);N=ceil(1.8*pi/kw);N=N+rem(N,2);使用窗函數(shù);window=Boxcar(N+1);設(shè)置濾波器類型：type=stop;得到濾波器系數(shù)：b=fir1(N,wn,'stop',window);

調(diào)入地震原始數(shù)據(jù)，對地震數(shù)據(jù)進(jìn)行濾波并輸出濾波前后波形圖。程序：

subplot(2,1,1),plot(t,a),xlabel(′t/s′),ylabel(′magnitude/count′),title(‘濾波前地震波形’),

subplot(2,1,2),plot(t,sf),xlabel(′t/s′),ylabel(′magnitude/count′),title(‘濾波后地震波形’);

512個點繪制濾波器頻譜，程序：

[h1,w]=freqz(b,1,512); %使用512個點繪制濾波器幅頻特性

subplot(2,1,1),plot(w*Fs/(2*pi),20*log(abs(h1)/abs(h1(1)))),xlabel('頻率/ HZ′),

ylabel('幅度/dB′),title(‘FIR帶阻濾波器幅頻特性’)。

2.2 IIR數(shù)字濾波器原理及設(shè)計

濾波器的設(shè)計本質(zhì)上是尋找一個既能物理實現(xiàn)，又能滿足給定頻率特性指標(biāo)要求的系統(tǒng)傳輸函數(shù)。IIR濾波器一般采用遞歸型的結(jié)構(gòu)，系統(tǒng)的輸入與輸出服從N階差分方程：

(7)

而IIR的傳輸函數(shù)為：

(8)

h(n)為濾波器的脈沖響應(yīng)，取值范圍在n∈[0,∞]。M和N為分解的分子和分母的系數(shù)個數(shù)[4]。IIR濾波器在同樣的性能指標(biāo)的前提下，需要的階數(shù)明顯要低于FIR濾波器。所用存儲單元少，計算量小，效率更高，但其缺點是相位的非線性。根據(jù)頻譜分析的結(jié)果以及干擾信號的頻率區(qū)間，需要設(shè)計IIR帶阻濾波器。設(shè)計思路是,利用已有的模擬濾波器設(shè)計理論，采用巴特沃斯butterworth模擬濾波器，然后再通過雙線性變換法，完成從模擬到數(shù)字的變換。確定濾波器的技術(shù)指標(biāo)，包括采樣頻率Fs，通帶和阻帶的上、下限歸一化角頻率(wlp、wls、whp、whs)，通帶衰減Rp,阻帶衰減Rs。通過計算求取頻率預(yù)畸變以及歸一化通、阻帶頻率，程序如下：

T=1/Fs; %采樣間隔

wc1=(2/T)*tan(wlp/2);wc2=(2/T)*tan(whp/2);wr1=(2/T)*tan(wls/2);

wr2=(2/T)*tan(whs/2); %頻率預(yù)畸變

w0=sqrt(wc1*wc2); R=wc2-wc1; wp=1; %歸一化通帶截止頻率

ws=wp*(wr1*R)/(w0^2-wr1^2); %歸一化阻帶截止頻率

用函數(shù)[N,wn]=buttord(wp,ws,Rp,Rs,'s')確定最小階數(shù)N和頻率參數(shù)Wn。再求出濾波器傳輸函數(shù)多項式系數(shù)向量b,a。相關(guān)程序：

[Z,P,K]=buttap(N); [MT,NT]=zp2tf(Z,P,K); %將零極點形式轉(zhuǎn)換為傳輸形式

[M,N]=lp2bs(MT,NT,w0,R); %對低通濾波器進(jìn)行頻率變換，轉(zhuǎn)換為帶阻濾波

[b,a]=bilinear(M,N,Fs); %對模擬濾波器進(jìn)行雙性變換。

繪出所設(shè)計濾波器幅幅頻響應(yīng)。相應(yīng)程序：

[H,W]=freqz(b,a); subplot(2,1,1),plot(W*Fs/(2*pi),abs(H)),xlabel(‘頻率’),ylabel(‘幅值’),title(‘IIR帶阻濾波器幅頻響應(yīng)’); %繪出濾波器幅頻響應(yīng)。調(diào)入地震原始數(shù)據(jù)，對地震數(shù)據(jù)進(jìn)行濾波并輸出濾波前后波形圖。程序如下：

subplot(2,1,1),plot(t,xt),xlabel(′t/s′),ylabel(′magnitude/count′),title(‘濾波前地震波形’),

subplot(2,1,2),plot(t,y),xlabel(′t/s′),ylabel(′magnitude/count′),title(‘濾波后地震波形’)。

3 對地震數(shù)據(jù)濾波處理實例

由于成都臺記錄到的中小震級遠(yuǎn)震受到干擾明顯，以圖1的4.7級遠(yuǎn)震為例，通過頻譜分析(參見圖2)，確定地震的優(yōu)勢頻率在0～1.5 Hz，干擾的優(yōu)勢頻率為2～85 Hz，因此設(shè)計了FIR帶阻濾波器與IIR帶阻濾波器分別處理。為了盡可能的保留有用的地震信息頻率，對兩種濾波器做如下參數(shù)設(shè)置：FIR帶阻濾波器：采樣頻率Fs=100 Hz;通帶上截止頻率flp=2 Hz;通帶下截止頻率fhp=5.5 Hz;阻帶下截止頻率fls=2.1 Hz;阻帶上截止頻率fhs=5 Hz;IIR帶阻濾波器：采樣頻率Fs=100 Hz;通帶上截止歸一化角頻率wlp=0.04*pi;通帶下截止歸一化角頻率whp=0.11*pi;阻帶下截止歸一化角頻率wls=0.042*pi;阻帶上截止歸一化角頻率whs=0.10*pi;通帶衰減Rp=1 dB;阻帶衰減Rs=3 dB;經(jīng)過程序運行后，得到濾波器的幅頻特性(參見圖3)，在頻譜上可以看到干擾的頻率被過濾掉(參見圖4)，波形也得到明顯的改善(參見圖5～6)。由圖4可以知，F(xiàn)IR、IIR濾波器應(yīng)用效果明顯，2～5 Hz的干擾被最大程度過濾，保留了地震信息的優(yōu)勢頻率。在圖5和圖6中可以看出，濾波前的原始地震波形記錄“毛刺”較多，波形受干擾嚴(yán)重，P波的初動方向無法清楚辨別。經(jīng)過濾波器濾波后，地震記錄改善十分明顯，之前的“毛刺”大大減少，波形更加清楚，P波的初動方向可以清晰辨認(rèn)。FIR與IIR兩種濾波器都取的了很好的效果，但對比圖5和圖6可以發(fā)現(xiàn)，圖5中地震的初動在時間上相對于原始波形稍有滯后，這是由于FIR濾波器的相位滯后造成的，但波形仍然保持原有形狀且相位無失真。在成都臺地震記錄中，除了中小震級的遠(yuǎn)震，近震和地方震波形都受到很大程度的干擾，甚至混夾干擾信號難以分辨。因此，本文對近震進(jìn)行了處理與分析，以2016年07月27日19時59分記錄于JCZ-1T的四川平武ML2.5級地震為例。首先通過對四川平武ML2.5級地震譜分析(參見圖7)，確定地震信號的優(yōu)勢頻率是在0～0.5 Hz，屬于0～1 Hz的區(qū)間；干擾的優(yōu)勢頻率在2～4 Hz的區(qū)間。因此，需要設(shè)計符合要求的帶阻濾波器去除干擾。改變技術(shù)參數(shù)使得：flp=2;fls=2.1;fhs=4;fhp=5;wlp=0.04*pi;wls=0.042*pi;whs=0.08*pi;whp=0.10*pi;運行程序后得到圖8。

圖3 FIR濾波幅頻特性及IIR濾波器幅頻響應(yīng)

左為FIR濾波結(jié)果；右為IIR濾波結(jié)果圖4 濾波前后的頻率譜

圖5 FIR濾波器濾波后的地震波形

圖6 IIR濾波器濾波后的地震波形

圖7 四川平武ML2.5級地震記錄的頻率譜

(圖左 FIR濾波結(jié)果 圖右IIR濾波結(jié)果)圖8 濾波前后的頻率譜圖

由圖8可以看出2～4 Hz的干擾頻率已經(jīng)被很大程度過濾掉了。圖9、圖10可知，F(xiàn)IR濾波器和IIR濾波器應(yīng)用效果均明顯，地震波形記錄經(jīng)過濾波器處理后，干擾的頻率被過濾，而地震信息的頻率得到了盡可能的保留。經(jīng)過濾波器處理后的地震波形記錄相較于原始地震記錄有了明顯的改善，干擾“毛刺”大大的減少，原始地震波形記錄中被淹沒的近震的波形也更加清晰。FIR濾波器濾波后的波形相位稍有滯后，但波形仍然保持原有形狀且相位無失真。

圖9 FIR濾波器濾波后的地震波形

圖10 IIR濾波器濾波后的地震波形

4 討論

根據(jù)文中列舉的記錄于JCZ-1T的地震實例的頻譜，可以發(fā)現(xiàn)干擾頻率都包含在2～4 Hz的區(qū)間。為了證實干擾是來自外界環(huán)境還是JCZ-1T儀器本身，本文截取了2016年02月19日21點19分記錄于JCZ-1的地震波形數(shù)據(jù)來作為對比研究。通過對數(shù)據(jù)進(jìn)行頻譜分析，發(fā)現(xiàn)仍然存在2-4 Hz的區(qū)間干擾(參見圖11)，并設(shè)計了相應(yīng)的IIR帶阻濾波器進(jìn)行處理，從處理結(jié)果(參見圖12)來看，波形有了明顯改善，干擾明顯減少。為了減小結(jié)果的偶然性，本文另截取了2016年06月18日19點51分記錄于JCZ-1T的地震波形數(shù)據(jù)作為比較。以得出的頻率譜(參見圖13)以及IIR帶阻濾波器處理前后的地震波形記錄(參見圖14)，證實了來自于外界環(huán)境2-4 Hz的干擾信息是造成成都臺地震波形記錄不佳的重要因素。

圖11 記錄于JCZ-1的地震記錄頻率譜

圖12 IIR濾波器濾波前后記錄于JCZ-1的地震記錄

圖13 記錄于JCZ-1T的地震記錄頻率譜

圖14 IIR濾波器濾波前后記錄于JCZ-1T的地震記錄

處理實例都是通過MatLab程序代碼來實現(xiàn)的，針對不同實例的頻譜對濾波器的相關(guān)參數(shù)做了改變，以達(dá)到過濾干擾和提高觀測資料質(zhì)量的目的。通過對成都地震臺測震數(shù)據(jù)的處理結(jié)果來看，收效甚好，確確實實改善了原來的數(shù)據(jù)質(zhì)量，通過研究過程中可總結(jié)出以下認(rèn)識：來自于外界環(huán)境2～4 Hz的干擾信息是造成成都臺地震波形記錄不佳的重要因素，通過設(shè)計帶阻數(shù)字濾波器可以去除干擾。FIR濾波器、IIR濾波器應(yīng)用于成都臺地震記錄，均能達(dá)到理想的濾波效果，濾除了2～4 Hz頻率區(qū)間的干擾，大大減小了地震波形記錄中的干擾信息，改善了數(shù)據(jù)的質(zhì)量，對成都地震臺中小震級地震的觀測記錄有了很大的幫助。FIR濾波器雖然有相位延遲，但是通過濾波器的波形仍然保持原有形狀且不會相位失真，并且相位線性，即不同頻率分量的信號通過FIR濾波器后時間差不變；IIR濾波器在相同性能指標(biāo)的前提下，相較于FIR濾波器所需階數(shù)更少，計算量更小，效率更高，幅頻特性精度更高，但相位非線性。在通過MatLab的濾波器方法對干擾信息進(jìn)行排除時，應(yīng)首先對干擾信息與背景噪聲的分析，明確需要設(shè)計的濾波器類型，根據(jù)實際需要選取通帶或者阻帶的邊界頻率，并且盡可能保留原始觀測信息。

[1] 田文德，葉建慶，胡俊明.成都地震臺JCZ-1與JCZ-1T甚寬頻帶地震儀對比觀測分析[J].地震研究，2013，3(36):372-378.

[2] 曾慶堂，起衛(wèi)羅，馬志剛，等.MatLab消除騰沖臺數(shù)字地震記錄中干擾波的應(yīng)用[J].華南地震，2014，1(34).

[3] 譚雨文，劉國明.MatLab在地震信號處理中的應(yīng)用實例[J].防災(zāi)減災(zāi)學(xué)報，2011，3(27).

[4] 萬永革. MatLab數(shù)字信號處理實例教程[M].北京：科學(xué)出版社，2012.

ApplicationandResearchofMatLabintheDigitalSeismicRecordsatChengduSeismicStandardStation

HE Siyuan, LIU Huajiao, XU Wenhai, ZHONG Libin, DONG Lijie

(Chengdu Seismic Standard Station, Sichuan Earthquake Agency, Sichuan Chengdu 611730, China)

In Chengdu Seismic Standard Station,the seismic data of the middle and small magnitude events contain more serious interference.In this paper,the software MatLab is used to process the digital seismic records in order to reduce the interference.Two software FIR and IIR filters in the MatLab are two types of wave interference filters.The results processed by these two filters are good.It is found that the interference with the frequency 2～4 Hz is an important cause of poor-quality.That is valuable for the work about analysis and processing of seismic data in the future.

MatLab；digital seismic records；interference wave；filter

2016-09-21

何思源(1991-)，男，四川省安岳縣人，助理工程師，主要從事測震分析.

P315.69

B

1001-8115(2017)04-0005-07

10.13716/j.cnki.1001-8115.2017.04.002