MatLab在成都地震臺數字地震記錄中的應用

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

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

MatLab在成都地震臺數字地震記錄中的應用

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

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

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

MatLab；數字地震記錄；干擾信息；濾波器

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

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

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

(1)

公式中計算出來的是信號x(t)的連續頻譜。但是，在實際的控制系統中能夠得到的是連續信號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次加法。當N較大時，這個計算量是很大的。而快速傅里葉變換(FFT)減少了DFT的運算次數,通過在時域將序列逐次分解為一組子序列，然后利用子序列的DFT來實現整個序列的DFT，從而減少了離散傅里葉變換的運算量，提高了計算效率。在MatLab中對輸入信號實現快速傅里葉變換的命令如下[4]：

Xf=fft(xt,N)

(3)

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

圖1 地震記錄的原始波形

圖2 地震記錄的頻率譜

2 FIR與IIR數字濾波器的設計

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

2.1 FIR數字濾波器原理及設計

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

(4)

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

(5)

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

(6)

信號通過FIR不失真條件是濾波器在通帶內具有恒定的幅頻特性和線性相位特性。線性相位FIR濾波器的相位滯后和群延遲在整個頻帶上是相等且不變的，所以，FIR濾波器雖然有相位延遲，但是通過濾波器的波形仍然保持原有形狀且不會相位失真。FIR帶阻濾波器的設計思路大致如下：確定帶阻濾波器性能參數，如采樣頻率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信號處理工具箱中采用經典窗函數法設計線性FIR數字濾波器的函數，調用格式：b=fir1(N,wn,'type',window)，其中b為濾波器系數；N為濾波器階數；wn為濾波器截止頻率；‘type’為設計濾波器的類型，如帶阻，type=stop，高通則type=high；window為窗函數類型。計算濾波器截止頻率wn，wn=[(wlp+wls)/(2*pi),(whp+whs)/(2*pi)]。計算濾波器階數N，程序如下:k1=wls-wlp;k2=whp-whs;kw=min(k1,k2);N=ceil(1.8*pi/kw);N=N+rem(N,2);使用窗函數;window=Boxcar(N+1);設置濾波器類型：type=stop;得到濾波器系數：b=fir1(N,wn,'stop',window);

調入地震原始數據，對地震數據進行濾波并輸出濾波前后波形圖。程序：

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數字濾波器原理及設計

濾波器的設計本質上是尋找一個既能物理實現，又能滿足給定頻率特性指標要求的系統傳輸函數。IIR濾波器一般采用遞歸型的結構，系統的輸入與輸出服從N階差分方程：

(7)

而IIR的傳輸函數為：

(8)

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

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); %頻率預畸變

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

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

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

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

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

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

繪出所設計濾波器幅幅頻響應。相應程序：

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

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 對地震數據濾波處理實例

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

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

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

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

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

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

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

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

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

4 討論

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

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

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

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

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

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

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

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

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

[4] 萬永革. MatLab數字信號處理實例教程[M].北京：科學出版社，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