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

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

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

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

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

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

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

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

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

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

(1)

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

(2)

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

Xf=fft(xt,N)

(3)

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

圖1 地震記錄的原始波形

圖2 地震記錄的頻率譜

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

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

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

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

(4)

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

(5)

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

(6)

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

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

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è)計(jì)

濾波器的設(shè)計(jì)本質(zhì)上是尋找一個(gè)既能物理實(shí)現(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ù)個(gè)數(shù)[4]。IIR濾波器在同樣的性能指標(biāo)的前提下，需要的階數(shù)明顯要低于FIR濾波器。所用存儲(chǔ)單元少，計(jì)算量小，效率更高，但其缺點(diǎn)是相位的非線性。根據(jù)頻譜分析的結(jié)果以及干擾信號(hào)的頻率區(qū)間，需要設(shè)計(jì)IIR帶阻濾波器。設(shè)計(jì)思路是,利用已有的模擬濾波器設(shè)計(jì)理論，采用巴特沃斯butterworth模擬濾波器，然后再通過雙線性變換法，完成從模擬到數(shù)字的變換。確定濾波器的技術(shù)指標(biāo)，包括采樣頻率Fs，通帶和阻帶的上、下限歸一化角頻率(wlp、wls、whp、whs)，通帶衰減Rp,阻帶衰減Rs。通過計(jì)算求取頻率預(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ù)多項(xiàng)式系數(shù)向量b,a。相關(guān)程序：

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

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

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

繪出所設(shè)計(jì)濾波器幅幅頻響應(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ù)，對(duì)地震數(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 對(duì)地震數(shù)據(jù)濾波處理實(shí)例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[4] 萬永革. MatLab數(shù)字信號(hào)處理實(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