基于非線性時(shí)頻分析的地震和爆破識(shí)別

張 帆， 朱新運(yùn)， 熊 丹， 鐘羽云， 陳曉燕

(1. 浙江省地震局， 杭州 310013； 2. 河北省地震局， 石家莊 050021）

0 引言

自開(kāi)展地震觀測(cè)以來(lái)， 如何區(qū)分地震及爆破就成為地震觀測(cè)研究的基礎(chǔ)課題， 并已經(jīng)產(chǎn)生了大量的科技成果， 隨著數(shù)字地震觀測(cè)、 計(jì)算機(jī)技術(shù)、 數(shù)學(xué)算法等的發(fā)展， 在最新的觀測(cè)及算法條件下， 通過(guò)相關(guān)指標(biāo)定性、 定量的識(shí)別地震及爆破有廣闊的研究空間。 在我國(guó)20 世紀(jì)70年代，可以從波形的圖像對(duì)比法、 振幅與周期關(guān)系法[1]等方法來(lái)區(qū)分地震和爆破。 8O年代以來(lái)， 傅淑芳、張誠(chéng)、 趙榮國(guó)等從波形、 波譜、 P 波初動(dòng)、 振幅比及尾波衰減等方面對(duì)爆破識(shí)別進(jìn)行了研究[2-4]， 得出了大致相同的結(jié)果。 很多工作人員把這些方法應(yīng)用到了工作中， 并利用波譜求出震源參數(shù)來(lái)區(qū)分地震和爆破[5-6]。 遺傳BP 網(wǎng)絡(luò)[7-8]、 分形分析[9]也被應(yīng)用于地震和爆破的識(shí)別。 以前廣泛應(yīng)用的波譜分析方法， 把地震波作為是平穩(wěn)信號(hào)， 利用傅立葉分析處理[10-11]。 隨著數(shù)字信號(hào)處理的發(fā)展， 研究者發(fā)現(xiàn)地震波是一種極為復(fù)雜的隨機(jī)波， 所產(chǎn)生的地震信號(hào)具有短時(shí)、 突變等特點(diǎn)， 是一種典型的非平穩(wěn)隨機(jī)信號(hào)[12]， 因而反映地震波非平穩(wěn)特征的線性時(shí)頻分析方法被應(yīng)用到地震分析及識(shí)別爆破中， 劉希強(qiáng)等提出小波變換能量線性度方法識(shí)別天然地震與爆破或塌方[13]， 還利用高斯線調(diào)頻連續(xù)小波變換的時(shí)頻能量衰減因子方法來(lái)識(shí)別[14]。孫甲寧， 夏愛(ài)國(guó)， 蘇乃秦用小波分析方法得到了烏魯木齊臺(tái)記錄的天然地震和人工爆破在時(shí)頻域的能量分布特征[15]。 筆者對(duì)比了波恩一約旦（Born-Jordan)、 喬伊一威廉斯（Choi-Williams）、 趙-阿特拉斯-馬克斯（Zhao-Atlas-Marks）、 維格納（WV）四種非線性時(shí)頻分析方法[16]， 發(fā)現(xiàn)趙-阿特拉斯-馬克斯（Zhao-Atlas-Marks）時(shí)頻分析方法， 抑制交叉干擾項(xiàng)的效果和時(shí)頻聚集性更好， 更適合于分析地震波。 本文利用趙-阿特拉斯-馬克斯（Zhao-Atlas-Marks）時(shí)頻分析方法分析浙江省的地震和爆破， 得到他們的時(shí)頻譜特征。

1 理論和方法

時(shí)頻分析方法， 主要是研究信號(hào)在時(shí)頻域上的能量密度是如何變化的， 它不僅可以同時(shí)用時(shí)間和頻率描述信號(hào)的能量密度， 而且能計(jì)算出任何密度。

所有的時(shí)頻分布都可以由下面的方程得到：

（1）式中φ (θ，τ）是二維函數(shù)叫做核[17]， 式中*表示解析信號(hào)復(fù)共軛。 取不同的核函數(shù)可得到不同的時(shí)頻分布。

由于時(shí)頻分布是信號(hào)的雙線性函數(shù)，兩個(gè)信號(hào)和的時(shí)頻分布并不是每個(gè)信號(hào)的時(shí)頻分布之和，而多了兩個(gè)附加項(xiàng)叫做交叉干擾項(xiàng)，（2）式中C11（t，ω）， C22（t， ω）， C12（t， ω）， C21（t， ω）分別是s1（t）， s2（t）的自項(xiàng)和他們之間的交叉干擾項(xiàng)， 且C12（t， ω）=C21（t， ω）。

（3）式中ω12=（ω1+ω2）/2， 通過(guò)交換角標(biāo)1 和2， 就可以得到C22和C21。 那么K（ω）的性質(zhì)決定了交叉干擾項(xiàng)的大小和分布[18]， 而核φ（0，τ）又唯一確定了K（ω）。 對(duì)于趙-阿特拉斯-馬克斯（Zhao-Atlas-Marks）時(shí) 頻 分 布， 取 核在趙一阿特拉斯一馬克斯的原著中， g（τ）取1， α 取1/2， 則

（4）式中可以看出K（ω）是ω－ω1和ω－ω2的函數(shù)，這些地方正好是自項(xiàng)的地方， 即交叉干擾項(xiàng)完全位于自項(xiàng)中間， 而其它時(shí)頻分布的交叉干擾項(xiàng)均未完全位于自項(xiàng)中間。 將（4）式代入（3）式， 求得K（ω）后， 代入（2）式， 即可得到趙-阿特拉斯-馬克斯（Zhao-Atlas-Marks）時(shí)頻分布。

此外， 對(duì)于任何的統(tǒng)計(jì)量x， 如果知道了它的分布P (x) 后， 可以求出該統(tǒng)計(jì)量的n 階矩， 即xn的平均值， 也就是

矩具有很多方面的重要意義： 首先， 一階矩給出了分布的基本性質(zhì)的表示； 其次， 知道的矩階數(shù)越高， 關(guān)于分布的了解就越多； 第三， 對(duì)于性能良好的分布， 矩唯一地確定了分布。 可以根據(jù)矩的理論求出信號(hào)和信號(hào)傅立葉變換的矩， 得到信號(hào)的時(shí)間局域特性和頻率局域特性。 在時(shí)間域中，

（6）式中tm為平均時(shí)間， T 為時(shí)間散布。

在頻率域中，

（7）式中X（ν）為x（t）的傅立葉變換， 其中fm為平均頻率， F 為頻率散布。

2 資料選取

選取浙江省數(shù)字化地震臺(tái)網(wǎng)的觀測(cè)資料來(lái)分析。 為了減小地震波在傳播過(guò)程中介質(zhì)的影響，盡量選取震中距最小的臺(tái)站記錄的爆破和地震進(jìn)行分析。 地震儀的頻帶不同所記錄到的地震波的波形也有所不同， 選取的地震波都是短周期地震儀記錄到的波形。 此外， 由于爆破的震級(jí)大部分都在ML3 級(jí)以下， 在挑選地震的時(shí)候， 也應(yīng)該挑選震級(jí)在ML3 級(jí)以下的地震來(lái)進(jìn)行對(duì)比分析。 根據(jù)以上的原則， 挑選了2000年3月至2006年5月間浙江省臺(tái)網(wǎng)記錄到的安吉爆破、 鄞縣地震、文成地震、 嵊州地震、 桐廬地震共51個(gè)事件進(jìn)行了分析和對(duì)比， 其中爆破11個(gè)， 地震40個(gè)， 震級(jí)范圍在ML1.5 和ML3.1 之間。 所選事件的震中分布和浙江省的臺(tái)站分布如圖1 所示。

選擇離震中最近的臺(tái)站的地震記錄來(lái)分析，臺(tái)站的震中距大部分都在20 km 以內(nèi)。 事件及所選臺(tái)站儀器參數(shù)見(jiàn)表1。

圖1 震中分布和浙江省臺(tái)站分布圖Fig.1 Distribution map of epicenters and stations

表1 所選臺(tái)站儀器參數(shù)表Table 1 The instrument parameter of selected stations

3 計(jì)算方法和結(jié)果

對(duì)每一個(gè)事件波形的三個(gè)分量分別經(jīng)過(guò)如下處理得到時(shí)頻分布： 第一步， 利用Matlab 中的detrend 函數(shù)去傾， 再讀入儀器零極點(diǎn)文件得到傳遞函數(shù)， 用簡(jiǎn)單易行的傅里葉分析法即用地震記錄的頻譜除以傳遞函數(shù)扣除儀器響應(yīng)后經(jīng)過(guò)反傅立葉變換[19]恢復(fù)地動(dòng)速度； 第二步， 根據(jù)儀器的頻帶范圍有選擇地進(jìn)行帶通濾波（本文中所選擇的濾波范圍為0.5～20 Hz）后， 按照累計(jì)梯形積分法，利用Matlab 中提供的cumtrapz 函數(shù)[20]對(duì)記錄進(jìn)行積分得到了位移記錄； 第三步， 從位移記錄中分別截取P 波和S 波進(jìn)行Hilbert 變換得到它們的解析信號(hào)后， 求出趙-阿特拉斯-馬克斯（Zhao-Atlas-Marks）時(shí)頻分布， 以及P 波和S 波的平均頻率fm、頻率散布F， 平均時(shí)間tm、 時(shí)間散布T。 最終， 時(shí)頻分布以等值線的形式給出。 典型的爆破和地震P波的時(shí)頻譜如圖2 和圖3 所示。

圖2 爆破位移記錄P 波時(shí)頻譜圖Fig.2 The P wave time-frequency spectrum of explosion displacement record

圖3 地震位移記錄P 波時(shí)頻譜圖Fig.3 The P wave time-frequency spectrum of earthquake displacement record

從圖2 和圖3 中可以看出， 相對(duì)于爆破P 波的時(shí)頻譜圖， 地震的時(shí)頻譜要分散一些。 爆破的P波的時(shí)頻譜峰值主要集中在1～8 Hz， 地震的P 波的時(shí)頻譜峰值主要集中在6～12 Hz。 爆破和地震S波的時(shí)頻譜圖分別如圖4 和5 所示。

圖4 爆破位移記錄S 波時(shí)頻譜圖Fig. 4 The S wave time-frequency spectrum of explosion displacement record

圖5 地震位移記錄S 波時(shí)頻譜圖Fig.5 The S wave time-frequency spectrum of earthquake displacement record

從圖4 和圖5 中同樣可以看出， 相對(duì)于爆破S波的時(shí)頻譜圖， 地震的時(shí)頻譜要分散。 爆破的S波的時(shí)頻譜峰值主要集中在2 Hz 左右。 地震的S波時(shí)頻譜峰值主要集中在2～12 Hz。 爆破的P 波時(shí)頻譜峰值主要分布在2 Hz 左右的區(qū)域里， 而地震的時(shí)頻譜峰值， 在頻率高和低的區(qū)域都有分布。 爆破的S 波時(shí)頻譜的最大峰值主要分布在較低的區(qū)域， 大約在1～2 Hz， 同樣地震的S 波多數(shù)時(shí)頻譜相對(duì)來(lái)說(shuō)有很多個(gè)峰值， 并且峰值出現(xiàn)的地方不如爆破有規(guī)律性， 有的峰值出現(xiàn)在頻率較高的區(qū)域， 有的峰值出現(xiàn)在頻率較低的區(qū)域。

地震的時(shí)頻譜相對(duì)與爆破來(lái)說(shuō)比較 “分散”，而爆破的時(shí)頻譜相對(duì)于地震來(lái)說(shuō)比較 “集中”。 S波最大振幅為S 波到達(dá)1～2 s 后， 主要頻率成分為2 Hz 左右， 這是由于S 波到達(dá)后， 短周期瑞利面波在此時(shí)出現(xiàn)。

爆破和地震的平均時(shí)間tm和平均頻率fm有明顯的不同。 對(duì)于P 波我們僅討論平均頻率。 爆破和地震的P 波平均頻率如圖6 所示。

圖6 爆破和地震位移記錄P 波的平均頻率（a：UD；b：EW；c：NS）Fig.6 The average frequency of the P wave of explosions and earthquakes displacement record (a：UD，b： EW，c：NS)

從表2 中看出， 爆破P 波的平均頻率比地震的平均頻率低很多， 甚至不是一個(gè)數(shù)量級(jí)的， 爆破P 波平均頻率fm為1.7～4 Hz， 地震S 波平均頻率為6～12 Hz。

分析S 波時(shí)， 同時(shí)考慮S 波的平均時(shí)間和平均頻率， 如圖7 所示。

表2 爆破和地震P 波平均頻率fm 對(duì)比表Table 2 The fm comparison table of P wave between explosions and earthquakes

從表3 看出， 爆破S 波三個(gè)分量的平均頻率fm為1～2 Hz， 不同分量平均時(shí)間tm有所不同， 約在1～11 s 之間， 地震S 波三個(gè)分量的平均頻率fm有所不同， 頻率范圍比較廣， 約在3～11 Hz 之間，平均時(shí)間tm在0.4～2.5 s 之間， 甚至有小于1 s 的事件， 而爆破都是大于1 s 的。 無(wú)論哪個(gè)分量， 爆破的時(shí)頻中心相對(duì)于地震分布在時(shí)間—頻率面的右下角， 而地震的時(shí)頻中心相對(duì)于爆破分布在時(shí)間—頻率面的左上角， 這種說(shuō)明了地震的頻率比爆破的頻率高， 地震的振幅比爆破的振幅先達(dá)到最大值。

4 結(jié)論與討論

圖7 爆破和地震位移記錄S 波的平均頻率和平均時(shí)間（a：UD；b：EW；c：NS）Fig.7 The average frequency and time of S wave of explosions and earthquakes displacement record (a：UD, b: EW, c: NS)

表3 爆破和地震S 波平均頻率fm 和平均時(shí)間tm 對(duì)比表Tabel 3 The comparison table of fm and tm of S wave between explosions and earthquakes

選取了浙江省數(shù)字化地震臺(tái)網(wǎng)提供的一些爆破和地震記錄共51個(gè)事件（其中爆破11個(gè)， 地震40個(gè)）， 經(jīng)過(guò)預(yù)處理、 扣除儀器響應(yīng)后， 分別求出位移記錄的ZAM 時(shí)頻分布后， 得到它們的時(shí)頻譜的差異如下：

（1）地震P 波的平均頻率大于爆破P 波的平均頻率。

（2）地震S 波的平均頻率大于爆破S 波的平均頻率， 地震S 波的平均時(shí)間小于爆破S 波的平均時(shí)間， 表現(xiàn)為地震的平均時(shí)頻中心相對(duì)于爆破分布在時(shí)間—頻率面上左上角。

（3）無(wú)論是從P 波還是從S 波的時(shí)頻譜看，地震的時(shí)頻譜有多個(gè)峰值， 而爆破的時(shí)頻譜僅有一個(gè)或者兩個(gè)峰值。

（4）爆破的時(shí)頻譜的時(shí)頻聚集性比地震的時(shí)頻聚集性好。

爆破和地震時(shí)頻譜這些差異和它們倆之間的震源、 介質(zhì)、 傳播過(guò)程等的不同是密切相關(guān)的。理想情況下， 爆炸源只產(chǎn)生P 波和次生瑞利面波（國(guó)外稱為Rg 波）， 爆炸點(diǎn)附近巖石的不均勻性導(dǎo)致比淺源地震弱的S 波[21]。 這種Rg 波在震中距7～90 km 的范圍內(nèi)其最大振幅明顯大于S 波最大振幅， 在爆破記錄中以優(yōu)勢(shì)波出現(xiàn)[22]， 而地震波形的各道記錄中卻沒(méi)有發(fā)現(xiàn)這種波[23]， Rg 波周期比S波周期大。 而在分析的地震和爆破中， 該震相在S波到達(dá)后的1～2 s 內(nèi)到達(dá)， 導(dǎo)致爆破S 波時(shí)頻譜的峰值出現(xiàn)在S 波到達(dá)后的1～2 s 內(nèi)， 頻率為1～2 Hz。 在時(shí)頻面上表現(xiàn)為爆破的平均時(shí)頻中心相對(duì)于地震來(lái)說(shuō)處位于時(shí)間—頻率平面右下角的區(qū)域。

地震震源較深， 頻率成分復(fù)雜， 而爆破的震源淺， 經(jīng)過(guò)松散地層， 高頻成分被吸收的多， 因此， 地震的P 波的平均頻率和S 波的平均頻率都比爆破的要高。 爆破的震源是膨脹源， 相對(duì)于地震來(lái)說(shuō)破裂和震源機(jī)制比較簡(jiǎn)單， 爆破所產(chǎn)生的波也相對(duì)于地震來(lái)說(shuō)比較簡(jiǎn)單， 頻率成分也比較單一， 時(shí)頻譜相對(duì)地震來(lái)說(shuō)聚集性較好。

時(shí)頻分析方法將地震和爆破看作一種非平穩(wěn)信號(hào)來(lái)處理， 可以很好地描述地震波的頻率隨時(shí)間的變化特征， 同時(shí)在時(shí)間尺度和頻率尺度描述地震波的能量分布， 而傅立葉分析僅僅只能在頻率域研究不同頻率的能量大小。 因此該方法除了可以識(shí)別爆破和地震以外， 還可以對(duì)地震和爆破的復(fù)雜性經(jīng)行研究。 但該方法的計(jì)算過(guò)程復(fù)雜，運(yùn)算速度有待提高。

[1] 張少泉． 近震分析[M]． 北京： 地震出版社， 1977．

[2] 傅淑芳， 劉寶誠(chéng)， 李文藝. 地震學(xué)教程(下冊(cè))[M]． 北京：地震出版社， 1980．

[3] 張誠(chéng)． 地震分析基礎(chǔ)[M]． 北京： 地震出版社， 1988.

[4] 時(shí)振粱， 強(qiáng)少泉， 趙榮國(guó)， 等． 地震工作手冊(cè)[M]. 北京： 地震出版杜， 1990.

[5] 張萍， 蔣秀琴， 苗春蘭， 等. 爆破、 礦震與地震的波譜差異[J]. 地震地磁觀測(cè)與研究， 2005， 26（3）： 24-34.

[6] 許健生， 尹志文. 地震與爆破的波譜差異[J]. 地震地磁觀測(cè)與研究， 1999， 20（3）： 18-23.

[7] 邊銀菊. 遺傳BP 網(wǎng)絡(luò)在地震和爆破識(shí)別中的應(yīng)用[J]. 地震學(xué)報(bào)， 2002， 24（5）： 516-524.

[8] Dowla F U ，Taylor S R ，Anderson R W. Seismic discrimination with artificial neural networks：Preliminary results with regional spectral data[J]． Bull Seism Soc A m ，1990， 80： 1346-1373.

[9] 劉代志， 鄒紅星， 韋蔭康. 分形分析與核爆地震模式識(shí)別[J]. 模式識(shí)別與人工智能， 1997， 10 (2)： 153-158.

[10] 林秀英， 張志呈. 爆破地震波的頻譜分析[J]. 中國(guó)礦業(yè)， 2000， 9（6）： 77-80.

[11] 郭學(xué)彬， 肖正學(xué)， 張志呈， 等. 爆破地震波的頻譜特性研究[J]. 化工礦物與加工， 1999， 7： 18-20.

[12] Dowding C H．Monitoring and control of Explosioning effects [M]. Mining Engineering Handbook. Prentice Hal1．1985， 746-760.

[13] 劉希強(qiáng)， 沈萍， 張玲， 等. 用小波變換能量線性度方法識(shí)別天然地震與爆破或塌方[J]. 西北地震學(xué)報(bào)，2003， 25（3）： 204-209.

[14] 劉希強(qiáng)， 沈 萍， 李 紅， 等. 高斯線調(diào)頻連續(xù)小波變換的時(shí)頻能量衰減因子方法及其應(yīng)用[J]. 中國(guó)地震，2003， 19（3）： 225-235.

[15] 孫甲寧， 夏愛(ài)國(guó)， 蘇乃秦. 地震和爆破時(shí)頻域能量分布特征的對(duì)比分析[J]. 華南地震， 2005， 25（2）： 68-74.

[16] 張 帆， 鐘羽云. 時(shí)頻分析方法及在地震波譜研究中的應(yīng)用[J]. 地震地磁觀測(cè)與研究， 2006， 27（4）： 17-22.

[17] T.A.C.M.Claasen and W.F.G. mecklenbr?uker. The Wigner distribution－a tool for time-frequency signal analysis－PartⅢ ：Relation with other time -frequency signal transformations [J]. Phlips Jour. Research，1998，35：372-389.

[18] L. Cohen. Time -Frequency Analysis: Theory and Application[M]. Englewood Cliffs：Prentice Hall， 1995.

[19] 陳運(yùn)泰， 吳忠良， 王培德， 等. 數(shù)字地震學(xué)[M]. 北京：地震出版社， 2000.

[20] 李海濤， 鄧櫻． Matlab 6.1 基礎(chǔ)及應(yīng)用技巧[M]． 北京：國(guó)防工業(yè)出版社， 2002．

[21] 趙永， 劉衛(wèi)虹， 高艷玲. 北京地區(qū)地震、 爆波和礦震的記錄圖識(shí)別[J]. 地震地觀測(cè)與研究， 1995， 16（4）：48-54.

[22] 山長(zhǎng)侖， 張玲， 李永紅， 等. 爆破短周期面波的分析[J]. 地震地磁觀測(cè)與研究， 2001， 22（6）： 43-48.

[23] 中國(guó)科學(xué)院地球物理研究所. 近震分析[M]. 北京： 地震出版杜， 1978.