曾慶堂，起衛(wèi)羅，馬志剛，余慶坤，王恩康，楊銘昌

(1.云南省地震局騰沖地震臺(tái)，云南 騰沖 679100；2.云南地震局元謀地震臺(tái)，云南 元謀 651399；3.云南省地震局應(yīng)急保障中心，云南 昆明 650224）

MATLAB消除騰沖臺(tái)數(shù)字地震記錄中干擾波的應(yīng)用

曾慶堂1，起衛(wèi)羅2，馬志剛1，余慶坤3，王恩康1，楊銘昌1

(1.云南省地震局騰沖地震臺(tái)，云南 騰沖 679100；2.云南地震局元謀地震臺(tái)，云南 元謀 651399；3.云南省地震局應(yīng)急保障中心，云南 昆明 650224）

數(shù)字化地震記錄資料中包含了大量的干擾信號(hào)，影響資料分析處理，有時(shí)還造成波形不完整或嚴(yán)重變形，很難進(jìn)行深入的分析。利用MATLAB設(shè)計(jì)數(shù)字濾波器對(duì)不同記錄中的干擾波進(jìn)行濾波，達(dá)到消除干擾的目的。實(shí)例表明，此方法具有良好的效果。

數(shù)字地震記錄；MATLAB；消除干擾；應(yīng)用

0 引言

“十五”以來，地震監(jiān)測發(fā)展了數(shù)字化觀測技術(shù)和網(wǎng)絡(luò)通信技術(shù)，極大地豐富了觀測的信息量，同時(shí)人們又面臨如何去處理分析使用這樣龐大的數(shù)字化資料的問題[1、2]。儀器在進(jìn)行地震數(shù)據(jù)采集的時(shí)候，也會(huì)采集到很多干擾信號(hào)，這些干擾嚴(yán)重影響了觀測資料的質(zhì)量，需要加以排除。排除干擾的方法有很多，數(shù)字濾波器就是其中使用最廣泛的一種。MATLAB是一套用于科學(xué)計(jì)算的可視化高性能的語言和軟件環(huán)境，它集數(shù)值分析、矩陣運(yùn)算、信號(hào)處理和圖形顯示于一體，工具箱中包含各種經(jīng)典和現(xiàn)代的數(shù)字信號(hào)處理技術(shù)，能實(shí)現(xiàn)各種數(shù)字濾波器的設(shè)計(jì)[5～7]。文中給出了IIR數(shù)字濾波器的MATLAB設(shè)計(jì)方法，并給出了用這種方法濾除數(shù)字地震記錄中干擾波的應(yīng)用實(shí)例。

1 濾波器設(shè)計(jì)原理

濾波器是一個(gè)選頻裝置，能傳輸有用信號(hào)，而又能抑制無用信號(hào)。有用信號(hào)和無用信號(hào)往往占有不同的頻帶，信號(hào)能通過濾波器的頻帶稱為通帶(Passband)，信號(hào)被抑制的頻帶稱為阻帶（Stopband）。理想濾波器頻率特性可寫為：

數(shù)字濾波器是數(shù)字信號(hào)處理技術(shù)的重要內(nèi)容。和模擬濾波器一樣，數(shù)字濾波器的主要功能是對(duì)數(shù)字信號(hào)進(jìn)行處理，最常見的處理是保留數(shù)字信號(hào)中有用的頻率成分，去除信號(hào)中無用的頻率成分[4、5]。

IIR濾波器的傳遞函數(shù)為：

h(n)為濾波器的脈沖響應(yīng)，從n=0～∞均有值。 M和N為分解的分子和分母多項(xiàng)式的系數(shù)個(gè)數(shù)[3～5]。

數(shù)字濾波器的設(shè)計(jì)目的是使濾波器的頻率特性達(dá)到所給定的性能指標(biāo)。其性能指標(biāo)包括通帶波紋RP(dB)、阻帶衰減 RS(dB)、通帶邊界頻率（Hz）、阻帶邊界頻率等。

IIR數(shù)字濾波器的設(shè)計(jì)方法:

IIR濾波器的設(shè)計(jì)就是根據(jù)濾波器的性能指標(biāo)要求，設(shè)計(jì)濾波器的分子和分母多項(xiàng)式系數(shù)。它和FIR濾波器相比優(yōu)點(diǎn)是在滿足相同性能指標(biāo)要求的前提下，IIR濾波器的階數(shù)要明顯低于FIR濾波器，但I(xiàn)IR濾波器的相位是非線性的。

IIR經(jīng)典設(shè)計(jì)就是將已設(shè)計(jì)好的模擬濾波器按一定變換原理轉(zhuǎn)換為數(shù)字濾波器。該方法先根據(jù)濾波器的技術(shù)指標(biāo)設(shè)計(jì)出相應(yīng)的模擬濾波器，然后再將設(shè)計(jì)好的模擬濾波器變換成數(shù)字濾波器。在MATLAB中，經(jīng)典法設(shè)計(jì)IIR數(shù)字濾波器主要采用以下步驟：

由圖1可見，IIR濾波器設(shè)計(jì)利用了模擬濾波器的設(shè)計(jì)成果。第一步和第二步完成后，一個(gè)達(dá)到期望性能指標(biāo)的模擬濾波器 （低通、高通、帶通、帶阻等）已經(jīng)設(shè)計(jì)出來。第三步模擬離散化的主要任務(wù)就是把模擬濾波器變換成數(shù)字濾波器，即把模擬濾波器的系統(tǒng)函數(shù)H(s)映射為數(shù)字濾波器的系統(tǒng)函數(shù)H(z)。這樣，數(shù)字濾波器的設(shè)計(jì)工作全部完成。

圖1 IIR濾波器的設(shè)計(jì)步驟Fig.1 The design procedures of IIR filter

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

MATLB中提供了大量設(shè)計(jì)IIR數(shù)字濾波器的函數(shù)，下面以設(shè)計(jì)Butterworth濾波器為例，介紹具體的設(shè)計(jì)步驟和方法。

（1）確定濾波器的技術(shù)指標(biāo)，主要有：采樣頻率FS，通帶臨界頻率FCP，阻帶臨界頻率FCS，通帶衰減RP和阻帶衰減RS。

有時(shí)候幾個(gè)孩子手上的鐵環(huán)沒有控制好方向，滾著滾著就互相碰撞起來，甚至也有的孩子連自己都跟著鐵環(huán)跑到了小路邊干涸的溝渠里，滾得滿頭滿臉都是塵土……

（2）用函數(shù) [N，Wn]=buttord(WP，WS，RP，RS，'hour')求出濾波器的最小階數(shù)N和標(biāo)準(zhǔn)化臨界頻率Wn；用函數(shù)[b，a]=lp2lp(Bap，Aap，Wn)求出濾波器傳遞函數(shù)多項(xiàng)式系數(shù)向量a和b。

（3）用函數(shù)freqz(bz，az，Nn，1/T)分析出所設(shè)計(jì)濾波器的幅頻特性和相頻特性，并用函數(shù)yt=filtfilt(bz，az，xt)完成對(duì)原始信號(hào)xt的濾波得到輸出信號(hào)yt。

（4） 用 subplot(2，1，1)，plot(t，x)，title('Input signal')繪出濾波器的性能圖。用subplot(2，1，2)，plot (t，y)，title('OutputSignal')，xlabel('Time/min')繪出信號(hào)濾波后的效果圖。

實(shí)現(xiàn)濾波功能的詳細(xì)MATLAB代碼如下：

以上代碼中，tct.txt為信號(hào)數(shù)據(jù)文件。

3 數(shù)字地震記錄中干擾波的排除應(yīng)用實(shí)例

數(shù)字地震記錄中記錄到的大部分是近震，由于干擾波或儀器本身問題，經(jīng)常記錄到的波形疊加了低頻或高頻成分，波形產(chǎn)生較大變形，很難分析，為把其中的干擾排除，針對(duì)干擾的不同頻率成分情況，設(shè)計(jì)不同帶通的濾波器進(jìn)行排除。

騰沖臺(tái)民振子臺(tái)記錄到的地震經(jīng)常疊加了低頻干擾，波形嚴(yán)重變形，圖2是所設(shè)計(jì)的濾波器（chebyshevI），通帶為1～5 Hz，過渡帶寬為0.9 Hz、2 Hz，通帶波紋小于1dB，阻帶衰減大于30 dB，采樣頻率為50 Hz。

圖2 所設(shè)計(jì)濾波器Fig.2 The developed filter

騰沖沙壩臺(tái)由于擺的問題，記錄到的地震P波和S波震相不明顯。下面設(shè)計(jì)了一個(gè)帶通chebyshevI型數(shù)字濾波器，通帶為2～5 HZ，通帶波紋小于1 DB，阻帶內(nèi)頻率小于1.3 HZ和大于6 HZ的波不能通過，采樣頻率為50 HZ。用此濾波器濾掉地震中背景噪聲，清晰的反映出P波、S波的震相特征。圖4為所設(shè)計(jì)的濾波器，令Fcp= [25]，F(xiàn)cs=[1.36]；圖5為輸入的波形和濾波后的輸出波形。

還可對(duì)其它數(shù)字化前兆儀器記錄數(shù)據(jù)中的干擾進(jìn)行排除，以下為一段騰沖臺(tái)重力分鐘值原始資料曲線及其濾波后的圖象。 原始資料曲線中，由于儀器故障,使得記錄曲線中在固體潮的低頻成分上疊加了一些高頻成分。為了得到固體潮成分的本來面目，設(shè)計(jì)了一低通濾波器，其幅頻特性和相頻特性見圖6，圖7是濾波前和濾波后的波形圖，經(jīng)運(yùn)行該濾波器程序后，取得良好效果。

以上實(shí)例都是運(yùn)行文中所給的MATLAB代碼，只是對(duì)不同記錄中，根據(jù)干擾頻率特性對(duì)通帶臨界頻率Fcp和阻帶臨界頻率Fcs作相應(yīng)改變，即可達(dá)到消除干擾的目的，方法極具靈活性和實(shí)用性。

從圖3、圖5和圖7中可以看出，所設(shè)計(jì)的濾波器均很好地排除了數(shù)據(jù)資料記錄中的干擾波，信號(hào)經(jīng)重構(gòu)后，得到較為理想的曲線，這樣既提高了觀測資料質(zhì)量，更易于進(jìn)一步提高地震波的分析及其前兆信息的分析處理的準(zhǔn)確性，對(duì)臺(tái)站的資料處理很實(shí)用。

4 結(jié)語

數(shù)字地震記錄中的干擾波有很多，要針對(duì)臺(tái)站的實(shí)際情況對(duì)各種干擾波進(jìn)行分析，找出其干擾頻率和特點(diǎn)，做為臺(tái)站的干擾背景資料，為精確地進(jìn)行數(shù)字化分析研究提供依據(jù)。

圖3 原始曲線圖和濾波后的曲線Fig.3 The original earthquake waveforms(a)and the filtered ones(b)

圖4 所設(shè)計(jì)的濾波器Fig.4 The developed filter

數(shù)字濾波器的應(yīng)用十分廣泛，運(yùn)用MATLAB語言能方便快捷地設(shè)計(jì)出滿足要求的濾波器。在排除數(shù)字地震記錄中的干擾波時(shí)，只要通過對(duì)臺(tái)站的正常背景信號(hào)和干擾背景信號(hào)做分析，準(zhǔn)確選取通帶臨界頻率Fcp和阻帶臨界頻率Fcs，就能將原始資料去其糟粕，取其精華，便于進(jìn)行深入和準(zhǔn)確的分析處理。

此方法不但適用于地震波的消噪處理，還適用于所有數(shù)字化資料的消噪處理，且簡單易行，實(shí)例表明，此方法具有良好的效果。

圖5 輸入波形和濾波后的輸出波形Fig.5 The waveforms before filtering(a)and after filtering(b)

圖6 濾波器的幅頻和相頻特性Fig.6 The amplitude-frequency characteristic(a)and the phase-frequency characteristic of the filter(b)

圖7 騰沖臺(tái)重力資料濾波前與濾波后的圖象Fig.7 The gravity records by Tengchong earthquake station before filtering(a)and after filtering(b)

[1]張少泉.地震波分析與應(yīng)用 [M].北京：地震出版社，1998.

[2]張平，山秀明，毛玉平，等.數(shù)字地磁信號(hào)的小波消噪處理[J].地震研究，2004，27(2)：198-202.

[3]徐東艷，孟曉剛編著.MATLAB函數(shù)庫查詢?cè)~典[M].北京：中國鐵道出版社，2006.

[4]董維國著.深入淺出MATLAB混合編程[M].北京：機(jī)械工業(yè)出版社，2006.

[5]萬永革.MATLAB數(shù)字信號(hào)處理實(shí)例教程[M].北京：科學(xué)出出版社，2004.

[6]張志涌，楊祖櫻，徐彥琴，等.精通MATLAB編程6.5版[M].北京：航空航天大學(xué)出版社，2003.

[7]李敬，甘延鋒，黃友明，等.數(shù)字地震記錄中干擾波的排除[J].防災(zāi)技術(shù)高等?？茖W(xué)校學(xué)報(bào)，2004，6（3）：20-25.

Application of MATLAB Software to the Elimination of the Disturbance Signals in the Digital Seismic Records at Tengchong Earthquake Station

ZENG Qingtan1，QI Weiluo2，MA Zhigang1，YU Qingkun3，WANG Enkang1，YANG Mingchang1

（1.Tenghcong Earthquake Station,Earthquake Administration of Yunnan Province,Tengchong 679100，China; 2.Yuanmou Earthquake Station，Earthquake Administration of Yunnan Province，Yuanmou 651399，China; 3.Earthquake Emergency and Logistics Center，Earthquake Administration of Yunnan Province，Kunming 650224，China）

Lots of disturbance signals such as low-frequency components and high-frequency ones contained in the digital seismic records often lead to the distorted or incomplete seismic waveforms,which are difficult for the seismologists to analyze.On the basis of MATLAB software,the authors developed digital filters to eliminate the disturbance signals in the digital seismic data recorded by Tengchong earthquake station.The filters is proved to be practical.

Digital seismic records；MATLAB software；Disturbance elimination；Application

P315.62

A

1001-8662（2014）01-0058-05

10.13512/j.hndz.2014.01.009

2013-03-08

曾慶堂，男，1968年生，工程師，主要從事地震監(jiān)測和火山監(jiān)測工作.

E-mail：zqingt@163.com.