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

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

0 引言

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

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

1.1 FIR濾波器的原理[3]

有限沖激響應(yīng)FIR數(shù)字濾波器的單位沖激響應(yīng)h（n）， 是有限長的 0≤n≤N-1 ， 長度為 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的有限個點（N個點）上有值。

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

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

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

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

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

窗函數(shù)設(shè)計法的基本原理是用有限長單位脈沖響應(yīng)序列逼近 hd（n）。 由于 hd（n）是無限長序列，而且是非因果的， 所以用窗函數(shù)w（n）將 hd（n）截斷，得到：

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

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

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

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

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

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

（4）檢驗濾波器的性能。

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

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

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

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

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

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

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

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

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

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

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

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

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

根據(jù)（5）使信號通過濾波器進行濾波。結(jié)果顯示信號通過濾波器后濾除了高頻振動干擾，濾波之前信號被高頻干擾所掩蓋，濾波之后可以很容易的分辨出地震波的震相。并且濾波后的信號相位無失真，濾波器達到了預(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

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

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

圖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é)語

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

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

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

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

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

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

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

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