許 可，劉瑞瑞，孔繁旭，朱 宏

（天津市地震局，天津 300201）

0 引言

在數(shù)字地震觀測(cè)系統(tǒng)記錄中，除了記錄到真實(shí)的地震信號(hào)外，還記錄到很多干擾信號(hào)。這些干擾信號(hào)嚴(yán)重影響了觀測(cè)資料的質(zhì)量需要加以排除[1]。隨著信息化社會(huì)的到來(lái)，計(jì)算機(jī)和微電子技術(shù)的迅速發(fā)展，數(shù)字信號(hào)處理的理論、算法以及實(shí)現(xiàn)手段都得到了全面的進(jìn)步，數(shù)字濾波器作為數(shù)字信號(hào)處理技術(shù)的一個(gè)重要工具，可用來(lái)過(guò)濾數(shù)字信號(hào)。因此，在處理數(shù)字地震觀測(cè)記錄中，可以利用濾波器來(lái)提取有用信號(hào)并抑制干擾信號(hào)。這是數(shù)字地震資料的分析和處理的重要內(nèi)容。Matlab是一種面向科學(xué)與工程數(shù)值計(jì)算的計(jì)算機(jī)軟件，工具箱中包含了各種經(jīng)典和現(xiàn)代數(shù)字信號(hào)處理技術(shù)，能實(shí)現(xiàn)各種數(shù)字濾波器的設(shè)計(jì)[2]。本文給出了FIR（有限沖激響應(yīng)）數(shù)字濾波器的Matlab設(shè)計(jì)方法，并給出了應(yīng)用這種方法消除天津數(shù)字地震記錄中干擾的應(yīng)用實(shí)例。

1 FIR數(shù)字濾波器的原理與設(shè)計(jì)

1.1 FIR濾波器的原理[3]

有限沖激響應(yīng)FIR數(shù)字濾波器的單位沖激響應(yīng)h（n）， 是有限長(zhǎng)的 0≤n≤N-1 ， 長(zhǎng)度為 N，（階數(shù)為N-1）的FIR系統(tǒng)函數(shù)為：

式（1）中 X （z）、 Y（z）分別為輸入 x（n）和輸出 y（n）的z變換，h（n）為濾波器的脈沖響應(yīng)，該式表示：濾波器的脈沖響應(yīng)h（n）在n=0，1，…，N-1的有限個(gè)點(diǎn)（N個(gè)點(diǎn)）上有值。

可得FIR濾波器的系統(tǒng)差分方程為:

因此，F(xiàn)IR濾波器又稱(chēng)為卷積濾波器，頻率響應(yīng)表達(dá)式為：

濾波器在通帶內(nèi)具有恒定的幅頻特性和線(xiàn)性相位特性。當(dāng)FIR濾波器的系數(shù)滿(mǎn)足下列中心對(duì)稱(chēng)條件： b（n）=b（N-1-n）或 b（n）=-b（N-1-n）時(shí)，濾波器設(shè)計(jì)在逼近平直幅頻特性的同時(shí)，還能獲得嚴(yán)格的線(xiàn)性相位特性。線(xiàn)性相位FIR濾波器的相位滯后和群延遲在整個(gè)品帶上是相等且不變的。對(duì)于一個(gè)N階的線(xiàn)性相位FIR濾波器，群延遲為常數(shù)，即濾波后的信號(hào)簡(jiǎn)單地延遲常數(shù)個(gè)時(shí)間步長(zhǎng)。這以特性使通帶頻率內(nèi)信號(hào)通過(guò)濾波器后仍保持原有波形形狀而無(wú)相位失真。

1.2 FIR濾波器窗函數(shù)方法

FIR濾波器主要的設(shè)計(jì)方法有窗函數(shù)法、最優(yōu)化設(shè)計(jì)法、約束最小二乘逼近法、升余弦函數(shù)法。本文采用窗函數(shù)法設(shè)計(jì)FIR濾波器。如果希望的理想頻率響應(yīng)函數(shù)為Hd（ejω），則其對(duì)應(yīng)的單位脈沖響應(yīng)為：

窗函數(shù)設(shè)計(jì)法的基本原理是用有限長(zhǎng)單位脈沖響應(yīng)序列逼近 hd（n）。 由于 hd（n）是無(wú)限長(zhǎng)序列，而且是非因果的， 所以用窗函數(shù)w（n）將 hd（n）截?cái)?，得到?/p>

h（n）就作為實(shí)際設(shè)計(jì)的FIR數(shù)字濾波器的單位脈沖響應(yīng)序列，其頻率響應(yīng)函數(shù)H（ejω）為:

式（6）中N為所選窗函數(shù)w（n）的長(zhǎng)度。

1.3 基于窗函數(shù)的FIR濾波器的設(shè)計(jì)步驟

（1）對(duì)濾波器的理想頻域幅值響應(yīng)進(jìn)行Fourier逆變換獲得理想濾波器的單位脈沖響應(yīng)hd（n）。

（2）由濾波器的性能指標(biāo)根據(jù)窗函數(shù)特點(diǎn)，確定滿(mǎn)足阻帶衰減的窗函數(shù)類(lèi)型w（n）。

（3）求實(shí)際濾波器的單位響應(yīng)h（n）。

（4）檢驗(yàn)濾波器的性能。

本文所用的窗函數(shù)是漢寧窗（hanning window），又稱(chēng)余弦窗，主瓣寬是8 π/N，第一旁瓣相對(duì)主瓣衰減-31 dB，它可以看做是3個(gè)矩形時(shí)間窗的頻譜之和，或者說(shuō)是3個(gè)sin（t）型函數(shù)之和，相當(dāng)于一個(gè)譜窗，向左、右各移動(dòng)了π/T，從而使旁瓣互相抵消，消除高頻干擾和漏能?？梢钥闯?，漢寧窗的主瓣加寬并降低，旁瓣則顯著減小（如圖1），從減小泄露的觀點(diǎn)出發(fā)，漢寧窗優(yōu)于矩形窗，但漢寧窗的主瓣加寬，相當(dāng)于分析帶寬加寬，頻率分辨率會(huì)相應(yīng)下降。通常來(lái)講漢寧窗的主瓣有較小的旁瓣和較大的衰減速度，是較為常用的窗函數(shù)。

2 FIR數(shù)字濾波器的Matlab實(shí)現(xiàn)

圖1 漢寧窗函數(shù)的幅頻形狀Fig.1 The amplitude-frequency shape of the hanning window function

Matlab已經(jīng)成為數(shù)字信號(hào)處理應(yīng)用中分析和仿真的主要工具，它的工具箱中提供了大量設(shè)計(jì)FIR數(shù)字濾波器的函數(shù)，這給濾波器的設(shè)計(jì)帶來(lái)了極大的方便。下面以基于漢寧窗函數(shù)的FIR數(shù)字濾波器為例，介紹具體的Matlab實(shí)現(xiàn)過(guò)程。

（1）確定濾波器的技術(shù)指標(biāo)，主要有采樣率Fs，通帶邊界頻率Wp，阻帶邊界頻率Ws，通帶波紋Rp，阻帶衰減Rs，過(guò)渡帶寬Wdelta。

（2）根據(jù)阻帶衰減確定窗函數(shù)類(lèi)型，使過(guò)渡帶寬近似于窗函數(shù)主瓣寬度，則可求得滿(mǎn)足性能指標(biāo)的窗函數(shù)的最小長(zhǎng)度。如果主瓣寬度為a，則可應(yīng)用函數(shù)N=ceil(a/Wdelta)求得窗函數(shù)的最小長(zhǎng)度。

（3）用函數(shù) b=fir1 （N， Wn， ‘ftype’， window）求出濾波器的傳遞函數(shù)多項(xiàng)式系數(shù)向量b。N為濾波器的階數(shù)，window為窗函數(shù)的列向量，其長(zhǎng)度為N+1。Wn為濾波器的截至頻率。對(duì)于帶通、帶阻濾波器Wn=[w1，w2]。w1和w2分別為帶通、帶阻濾波器的邊界頻率。對(duì)于低通、高通濾波器Wn=(Wp+Ws)/2。 ‘ftype’為濾波器的類(lèi)型。

（4）用函數(shù) [H， f]=freqz （b， 1， 512， Fs）分析出所設(shè)計(jì)的濾波器的幅頻特性和相頻特性。

（5）用函數(shù) y=filtfilt（b，1）完成對(duì)原始信號(hào) x 的濾波得到輸出信號(hào)y，并進(jìn)行Fourier變換與原信號(hào)進(jìn)行比較。

3 消除天津數(shù)字地震記錄中的干擾

天津位于華北平原東北部，與河北省和北京市為鄰，東鄰渤海，全境絕大部分屬華北平原，數(shù)字地震記錄到的大部分是近震。天津數(shù)字地震臺(tái)網(wǎng)各子臺(tái)分布在天津18個(gè)區(qū)縣的學(xué)校、醫(yī)院、工廠、村莊等。由于天津東部臨近渤海海域，南部是經(jīng)濟(jì)和工業(yè)發(fā)達(dá)的地區(qū)，人口稠密。因此，經(jīng)常記錄到的波形疊加了高頻或低頻的成分。針對(duì)干擾的不同頻率成分，可設(shè)計(jì)基于窗函數(shù)的FIR數(shù)字濾波器進(jìn)行濾波，消除干擾。

2014年3月30日23時(shí)11分，天津漢沽發(fā)生ML2.3級(jí)地震。天津沙井子臺(tái)北南向，記錄到的原始信號(hào)中夾雜著多種干擾，掩蓋了地震信息。通過(guò)快速Fourier變換（fft）分析可知，波形11～45 Hz頻段內(nèi)存在許多高頻振動(dòng)干擾，主要是受人為噪聲干擾影響。沙井子地震臺(tái)位于天津最南部的大港工業(yè)區(qū)，距離沙井子臺(tái)1 km處有一個(gè)風(fēng)力發(fā)電廠，每天有許多風(fēng)力發(fā)電機(jī)運(yùn)行，地震波形的Fourier變換上有許多比較有規(guī)律的高頻干擾，如圖2。

根據(jù)波形干擾的頻譜范圍特點(diǎn)，可設(shè)計(jì)帶阻FIR濾波器，通帶的起始頻率w1=11 Hz，截至頻率 w2=45 Hz，阻帶衰減 Rs=30 dB，過(guò)渡帶寬Wdelta=0.5 Hz，采樣率 Fs=100 Hz。根據(jù)上述Matlab實(shí)現(xiàn)過(guò)程（2），可確定選取的窗函數(shù)為漢寧窗，F(xiàn)IR濾波器階數(shù)800。由（4）可得到帶阻濾波器的幅頻特性和相頻特性。如圖4。

根據(jù)（5）使信號(hào)通過(guò)濾波器進(jìn)行濾波。結(jié)果顯示信號(hào)通過(guò)濾波器后濾除了高頻振動(dòng)干擾，濾波之前信號(hào)被高頻干擾所掩蓋，濾波之后可以很容易的分辨出地震波的震相。并且濾波后的信號(hào)相位無(wú)失真，濾波器達(dá)到了預(yù)期的要求。采用Fourier變換分析濾波后的振幅譜，其中阻帶內(nèi)的頻率成分被濾除掉。如圖5。

圖2 原始波形及Fourier變換Fig.2 The original waveforms and the Fourier transform

圖3 帶阻濾波器的幅頻和相頻特性Fig.3 The amplitude frequency and phase frequency characteristics of Band-stop filter

圖4 濾波后的波形及Fourier變換Fig.4 The filtered waveforms and the Fourier transform

天津糙甸臺(tái)經(jīng)常記錄到的波形疊加了低頻干擾，波形嚴(yán)重變形，很難分析?？稍O(shè)計(jì)高通FIR數(shù)字濾波器，通帶起始頻率為5 Hz，阻帶截止頻率為3 Hz，過(guò)渡帶寬為4 Hz，阻帶衰減大于30 dB，采樣率為100 Hz。濾波器的幅頻與相頻特性，如圖6。

可以看到，該濾波器小于5 Hz的低頻部分不能通過(guò)，大于5 Hz的高頻部分可以通過(guò)。較好的濾除了小于5 Hz的低頻干擾，清晰的反映出P波、S波的震相特征。圖6為原始波形與濾波后的波形對(duì)比。

圖5 高通濾波器的幅頻和相頻特性Fig.5 The amplitude frequency and phase frequency characteristics of High-pass filter

圖6 原始波形及濾波后波形Fig6 The original waveforms and the filtered ones

4 結(jié)語(yǔ)

隨著信息時(shí)代與數(shù)字技術(shù)的發(fā)展，數(shù)字信號(hào)處理己逐漸發(fā)展成為當(dāng)今極其重要的學(xué)科與技術(shù)領(lǐng)域之一。數(shù)字濾波器是數(shù)字信號(hào)處理的重要基礎(chǔ)，在對(duì)信號(hào)的濾波、檢測(cè)及參數(shù)的估計(jì)等信號(hào)應(yīng)用中，數(shù)字濾波器是使用最為廣泛的一種線(xiàn)性系統(tǒng)。數(shù)字濾波器根據(jù)其單位沖擊響應(yīng)函數(shù)的時(shí)域特性可分為兩類(lèi)：無(wú)限沖擊響應(yīng)（IIR）數(shù)字濾波器和有限沖擊響應(yīng)（FIR）數(shù)字濾波器。與IIR數(shù)字濾波器比，F(xiàn)IR數(shù)字濾波器的實(shí)現(xiàn)是非遞歸的，穩(wěn)定性好，精度高；更重要的是FIR數(shù)字濾波器在滿(mǎn)足幅度響應(yīng)要求的同時(shí)，可以獲得嚴(yán)格的線(xiàn)性相位。因此，它在高保真的信號(hào)處理中可以到廣泛應(yīng)用

本文將FIR數(shù)字濾波器應(yīng)用于排除地震記錄中的干擾信號(hào)，運(yùn)用MATLAB語(yǔ)言能容易地設(shè)計(jì)出滿(mǎn)足要求的FIR數(shù)字濾波器，設(shè)計(jì)簡(jiǎn)單方便，大大減少了工作量。在應(yīng)用中只需對(duì)程序中濾波器的起始頻率、截止頻率、采樣率和窗函數(shù)等參數(shù)進(jìn)行修改就可以實(shí)現(xiàn)需要的濾波功能，達(dá)到較好的濾波效果。綜上所述，F(xiàn)IR數(shù)字濾波器可以很好的排除地震記錄的干擾信號(hào)，信號(hào)經(jīng)重構(gòu)后，可以得到較為理想的曲線(xiàn)，這樣既提高了觀測(cè)資料質(zhì)量，又易于進(jìn)一步對(duì)資料進(jìn)行深入準(zhǔn)確的處理[4]，為今后數(shù)字地震資料的分析與應(yīng)用奠定基礎(chǔ)。

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

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

[3]陳杰.MATLAB典[M].北京：電子工業(yè)出版社，2004.

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

[5]魯權(quán)，邢西淳，李西京，等.涇陽(yáng)臺(tái)數(shù)字地磁信號(hào)的干擾分析及去噪處理 [J].華南地震，2013，33（1）：49-54.

[6]肖攀，宋國(guó)華，李露露，等.蚌埠地震臺(tái)測(cè)震干擾分析及處理[J].華南地震，2013，33（1）：86-92.