張 帆， 朱新運， 熊 丹， 鐘羽云， 陳曉燕

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

0 引言

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

1 理論和方法

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

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

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

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

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

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

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

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

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

在頻率域中，

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

2 資料選取

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

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

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

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

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

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

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

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

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

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

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

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

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

爆破和地震的平均時間tm和平均頻率fm有明顯的不同。 對于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 波的平均頻率比地震的平均頻率低很多， 甚至不是一個數(shù)量級的， 爆破P 波平均頻率fm為1.7～4 Hz， 地震S 波平均頻率為6～12 Hz。

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

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

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

4 結(jié)論與討論

圖7 爆破和地震位移記錄S 波的平均頻率和平均時間（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 和平均時間tm 對比表Tabel 3 The comparison table of fm and tm of S wave between explosions and earthquakes

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

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

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

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

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

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

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

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

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

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

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

[4] 時振粱， 強少泉， 趙榮國， 等． 地震工作手冊[M]. 北京： 地震出版杜， 1990.

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

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

[7] 邊銀菊. 遺傳BP 網(wǎng)絡(luò)在地震和爆破識別中的應(yīng)用[J]. 地震學(xué)報， 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] 劉代志， 鄒紅星， 韋蔭康. 分形分析與核爆地震模式識別[J]. 模式識別與人工智能， 1997， 10 (2)： 153-158.

[10] 林秀英， 張志呈. 爆破地震波的頻譜分析[J]. 中國礦業(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] 劉希強， 沈萍， 張玲， 等. 用小波變換能量線性度方法識別天然地震與爆破或塌方[J]. 西北地震學(xué)報，2003， 25（3）： 204-209.

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

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

[16] 張 帆， 鐘羽云. 時頻分析方法及在地震波譜研究中的應(yīng)用[J]. 地震地磁觀測與研究， 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] 陳運泰， 吳忠良， 王培德， 等. 數(shù)字地震學(xué)[M]. 北京：地震出版社， 2000.

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

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

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

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