利用MATLAB信號(hào)處理工具箱SPTool完成地震資料濾波處理

孫 煒

（德都地震臺(tái)，黑龍江 五大連池 164155）

0 引言

數(shù)字濾波器是一種用來(lái)過(guò)濾時(shí)間離散信號(hào)的數(shù)字系統(tǒng)，它是通過(guò)對(duì)抽樣數(shù)據(jù)進(jìn)行數(shù)學(xué)處理來(lái)達(dá)到頻域?yàn)V波的目的。數(shù)字濾波器可以用軟件或設(shè)計(jì)專用的數(shù)字處理硬件兩種方式來(lái)實(shí)現(xiàn)。在用軟件來(lái)實(shí)現(xiàn)數(shù)字濾波器中，可以根據(jù)濾波對(duì)象的實(shí)際情況，設(shè)置所需要的濾波器參數(shù)，改變?yōu)V波器的性能，達(dá)到所需要的濾波效果。在臺(tái)站日常地震資料處理中，由于各種環(huán)境因素的影響，可能使得震相圖中出現(xiàn)區(qū)別于地震信息的其他干擾信號(hào)，如車輛、湖面水波動(dòng)、爆破等等［1］。因此，對(duì)地震資料的濾波處理是一個(gè)十分重要的環(huán)節(jié)，有助于提高震相分析的準(zhǔn)確性。要實(shí)現(xiàn)對(duì)地震資料的濾波處理，必須設(shè)計(jì)相應(yīng)的濾波器。對(duì)于數(shù)字濾波器的軟件實(shí)現(xiàn)，可以借助高級(jí)語(yǔ)言（如Fortran、MATLAB語(yǔ)言）來(lái)實(shí)現(xiàn)，但是這均需要比較好的計(jì)算機(jī)編程能力，對(duì)于臺(tái)站年齡較大的計(jì)算機(jī)基礎(chǔ)或信號(hào)處理基礎(chǔ)較薄弱的同志，很難單獨(dú)實(shí)現(xiàn)，筆者直接使用MATLAB中自帶的信號(hào)處理工具箱SPTool工具較方便的實(shí)現(xiàn)了數(shù)字濾波器的設(shè)計(jì)，可以省略了復(fù)雜的濾波器程序編寫(xiě)過(guò)程［2］，并完成了地震信號(hào)的濾波處理。

1 數(shù)字濾波器的基本原理

從數(shù)學(xué)意義來(lái)說(shuō)，數(shù)字濾波器是一個(gè)離散的時(shí)間系統(tǒng)，輸入x（n）是一個(gè)時(shí)間序列，輸出y（n）也是一個(gè)時(shí)間序列，如數(shù)字濾波器脈沖響應(yīng)為h（n），時(shí)間域內(nèi)存在下列關(guān)系：y（n）=x（n）?y（n），頻率域內(nèi)輸入輸出關(guān)系為 Y（jw）=X（jw）H（jw）。同時(shí)濾波器是一個(gè)選頻裝置。理想的濾波器應(yīng)該能無(wú)失真地傳輸有用信號(hào)，而又能完全的抑制掉無(wú)用信號(hào)。有用信號(hào)和無(wú)用信號(hào)往往占有不同的頻帶。信號(hào)能夠通過(guò)濾波器的頻帶成為通帶（passband），信號(hào)被抑制的頻帶稱為阻帶（stopband）。理想濾波器的頻率特征可以寫(xiě)為：［3］

理想濾波器是物理上不可實(shí)現(xiàn)的系統(tǒng)，實(shí)際濾波器的頻率特性只能夠 “逼近”理想濾波器。按頻率域特性來(lái)講，數(shù)字濾波器可以分為低通、高通、帶通、帶阻等。由于實(shí)際的濾波器只能 “逼近”理想濾波器，所以，實(shí)際中濾波器的幅頻響應(yīng)不能如（1）、（2）式那樣理想，通帶內(nèi)部不是完全平直，而是呈波紋變化；在阻帶內(nèi)，其幅頻特性也不為零，而是衰減至某個(gè)值；在通帶和阻帶之間還存在一個(gè)過(guò)渡帶，也并不是突然下降。因此設(shè)計(jì)濾波器的技術(shù)指標(biāo)就包括了通帶波紋Rp、阻帶衰減Rs、阻帶邊界頻率ws、通帶邊界頻率wp。濾波器通帶波紋Rp為相對(duì)于頻率響應(yīng)最大點(diǎn)的下降，因此下降越少越好，越少說(shuō)明通帶越平直，濾波效果越好，阻帶衰減Rs也是相對(duì)頻率響應(yīng)最大點(diǎn)的下降，下降越多說(shuō)明信號(hào)在阻帶里受到越大的抑制，濾波效果越好。因此在實(shí)際應(yīng)用中只能根據(jù)具體實(shí)際需要來(lái)設(shè)計(jì)濾波器。一般來(lái)說(shuō)數(shù)字濾波器分為無(wú)限沖激和有限沖激（脈沖）響應(yīng)數(shù)字濾波器兩類。兩種濾波器各有優(yōu)缺點(diǎn)，在此不再詳述，利用MATLAB的SPTool均能根據(jù)濾波器的性能要求實(shí)現(xiàn)以上兩種濾波器設(shè)計(jì)，文中后面主要以IIR濾波器設(shè)計(jì)為例。

2 SPTool實(shí)現(xiàn)濾波器設(shè)計(jì)

一般情況下，如果不使用信號(hào)處理工具箱SPTool來(lái)設(shè)計(jì)濾波器，嚴(yán)格的數(shù)字濾波器設(shè)計(jì)，以IIR濾波器為例，在MATLAB中一般IIR數(shù)字濾波器設(shè)計(jì)的步驟如下：（1）將給定的數(shù)字濾波器的性能指標(biāo)進(jìn)行轉(zhuǎn)換，轉(zhuǎn)換后的頻率指標(biāo)作為模擬濾波器原型設(shè)計(jì)指標(biāo)；（2）估計(jì)模擬濾波器的最小階次和邊界頻率；（3）設(shè)計(jì)模擬低通濾波器原型；（4）由模擬原型經(jīng)頻率變換得到實(shí)際的模擬（低通、高通、帶通、帶阻）濾波器；（5）將模擬濾波器轉(zhuǎn)換為IIR數(shù)字濾波器，可利用脈沖響應(yīng)不變法或雙線性變換法，以上步驟對(duì)于熟悉信號(hào)處理以及具備計(jì)算機(jī)高級(jí)語(yǔ)言編程能力的人來(lái)說(shuō)，能夠單獨(dú)完成全部過(guò)程，但是對(duì)于不懂計(jì)算機(jī)語(yǔ)言編程的人來(lái)說(shuō)，以上過(guò)程十分的繁瑣復(fù)雜，倘若利用MATLAB自帶的仿真工具SPTool，就會(huì)使得以上繁瑣的濾波器設(shè)計(jì)過(guò)程變得十分簡(jiǎn)單。SPTool工具是MATLAB自帶的一個(gè)交互式信號(hào)處理工具，專門(mén)用于完成常用的數(shù)字信號(hào)處理任務(wù)。通過(guò)這個(gè)工具，只要按SPTool界面要求，基層臺(tái)站工作人員均能設(shè)計(jì)好實(shí)際工作中需要的濾波器，并輸入相應(yīng)的地震數(shù)據(jù)文件，就可以完成十分復(fù)雜的地震數(shù)字信號(hào)濾波處理任務(wù)，而不需要用戶對(duì)數(shù)字濾波器的設(shè)計(jì)原理非常熟悉。下面實(shí)例說(shuō)明筆者使用MATLAB中的交互式的信號(hào)處理工具SPTool濾波器設(shè)計(jì)。首先，我們?cè)O(shè)計(jì)一個(gè)阻帶邊界頻率為1Hz，通帶邊界頻率為2Hz，通帶波紋為1dB，阻帶衰減為30dB，采用Butterworth濾波器。點(diǎn)擊進(jìn)入SPTool中（Filter Design）功能模塊，如圖1所示。

圖1 濾波器設(shè)計(jì)界面Fig.1 The printscreen of filter designing

可以指定濾波器的性能參數(shù)（如采樣頻率、通帶頻率范圍等），也可以選擇設(shè)計(jì)函數(shù)即設(shè)計(jì)方法，如Equiripple FIR，LeastSquares FIR，切比歇夫IIR，橢圓IIR等。濾波器設(shè)計(jì)完之后，另一個(gè)重要的事情是進(jìn)行濾波器性能分析，SPTool提供了進(jìn)行濾波器的頻域和時(shí)域分析，（Filter View）是一個(gè)基于圖形用戶界面的分析工具，可以顯示特定濾波器6個(gè)不同的特征子圖，包括幅值響應(yīng)、相位響應(yīng)、群延遲、零極點(diǎn)和脈沖響應(yīng)、階躍響應(yīng)。其中相位響應(yīng)是濾波器頻率響應(yīng)的角度分量，零極點(diǎn)圖用于顯示濾波器傳遞函數(shù)的零極點(diǎn)與單位圓的位置。下圖2是（Filter View）給出的幅頻響應(yīng)與相位響應(yīng)。點(diǎn)選圖中 group delay、phase、zeros and poles、impulse response、 step response等可以得到群延遲、零極點(diǎn)和脈沖響應(yīng)、階躍響應(yīng)等其他圖形。

圖2 濾波器特征分析界面Fig.2 The printscreen of filter character analyzing

3 地震實(shí)例濾波

臺(tái)站每日觀測(cè)得到的震相數(shù)據(jù)中，除了由于地震自身所引發(fā)的之外，還包含了地脈動(dòng)、海浪干擾等低頻干擾，給地震波震相標(biāo)注帶來(lái)了麻煩。下面的實(shí)例，就給出了如何利用上面SPTool設(shè)計(jì)的高通濾波器進(jìn)行低頻干擾濾除。數(shù)據(jù)來(lái)源為五大連池地震火山監(jiān)測(cè)站2007年11月9日16時(shí)57分一個(gè)近地方震，利用SPTool工具中Import功能導(dǎo)入地震數(shù)據(jù)，導(dǎo)入數(shù)據(jù)必須先將地震數(shù)據(jù)格式轉(zhuǎn)換為MATLAB中的MAT文件。SPTool中（Signal Browser）會(huì)顯示所導(dǎo)入的地震波形，如圖3所示。然后在SPTool中點(diǎn)選APPLY利用前述濾波器進(jìn)行濾波，濾波時(shí)可以選擇零相位差濾波，消除IIR濾波器的時(shí)域相位延遲，如圖4所示。濾波之后，SPTool中（Signal Browser）顯示出濾波結(jié)果，如圖5所示。通過(guò)濾波結(jié)果可以看到，地震波中低頻干擾得到了很好的消除，突出了明顯的近地方震震相，有助于臺(tái)站地震分析人員明確其震相。

圖3 導(dǎo)入地震波形信號(hào)顯示Fig.3 Importing seismic wave data and displaying

圖4 應(yīng)用濾波器對(duì)導(dǎo)入的地震信號(hào)濾波Fig.4 The filtering of the importing seismic signal by the wave filter

圖5 地震波濾波結(jié)果顯示Fig.5 The result of filtering seismic wave data

4 結(jié)語(yǔ)和討論

通過(guò)利用MATLAB中信號(hào)處理工具SPTool可以十分容易的設(shè)計(jì)出各種類型的數(shù)字濾波器，而且其良好的圖形顯示界面與簡(jiǎn)單的操作，使具備初級(jí)信號(hào)處理知識(shí)的人也能完成從濾波器設(shè)計(jì)到信號(hào)濾波的整個(gè)流程，使臺(tái)站人員也能夠進(jìn)行各種干擾分析、干擾剔除等工作，一定程度上提高了臺(tái)站震相識(shí)別能力。

［1］中國(guó)地震局監(jiān)測(cè)預(yù)報(bào)司.地震學(xué)與地震觀測(cè)［M］.北京：地震出版社，2007.

［2］丁磊，潘貞存，叢偉，等.基于 MATLAB信號(hào)處理工具箱的數(shù)字濾波器設(shè)計(jì)與仿真繼電器［J］.繼電器，2003，31（9）.

［3］萬(wàn)永革.數(shù)字信號(hào)處理的MATLAB實(shí)現(xiàn)［M］.北京：科學(xué)出版社，2007.