2017, Vol. 38
表1(Table 1)
非平稳信号的时频分析就是在时间和频率的二维域内展示信号频率随时间的变化情况。时频分析方法一般有线性和非线性两类，其中典型线性分析方法有短时傅里叶变换、小波变换和Gabor变换（崔鑫等，2016）。为了更好地理解语音信号，Potter等在1947年首次提出一种实用的时频分析方法——短时傅里叶变换（王新刚等，2013）。短时傅里叶变换计算速度快，曹妃甸测震台网采用该时频分析方法分析地震波和爆破，具有较好的分辨率。
唐山迁安地区矿产丰富，爆破频繁，曹妃甸测震台网每天均能记录到矿区爆破波形，还监测到2015年8月12日天津滨海化学爆破及2016年1月6日朝鲜核试验爆破。分辨地震和爆破的常用手段有很大主观性和经验性，判别依据有:①初动方向:爆破初动方向一般向上，地震则不确定；②振幅:爆破的P波振幅大于S波，地震与其相反；③面波发育:地震面波一般不发育，震中距较大的爆破会发育明显面波。本文运用时频分析方法，研究曹妃甸测震台网记录的核爆、化爆、矿爆和天然地震，综合分析其时频分布特点，通过对比分析，寻找规律，以提高该台网识别爆破与地震的能力。
1 基本理论短时傅里叶变换又称为窗式傅里叶变换，是常用时频分析方法，基本思想是将非平稳信号分解成一系列短时平稳信号。具体实现方法为:用一个窗函数乘时间信号，对窗函数内的信号进行傅里叶变换，窗函数在时间轴上滑动，即可计算整个非平稳信号频率随时间的变化。即
| $ {\rm{STF}}{{\rm{T}}_x}\left({n, \omega } \right) = \int\limits_{ - \infty }^\infty {x\left(m \right)} w\left({n - m} \right){{\rm{e}}^{ - j\omega m}}{\rm{d}}m $ |
或
| $ {\rm{STF}}{{\rm{T}}_x}\left({n, \omega } \right) = \sum\limits_{m = - \infty }^\infty {x\left(m \right)w\left({n - m} \right){{\rm{e}}^{ - j\omega m}}} $ |
式中:w (n-m)为采用的时间窗，n为时间点，ω为频率点，m为变换参数。
2 资料选取河北省地震局曹妃甸地震台网由34个网络台和3个子台（曹妃甸台、南堡台、月岛台）组成，始建于2010年，2013年建成投入运行，主要观测环渤海地区地震活动。曹妃甸台附近构造主要有西南庄断裂、高柳断裂、柏各庄断裂、渤中断裂、唐山断裂等（赵建明等，2015）。曹妃甸测震台网台站分布见图 1。
|
图 1 所选地震及爆破事件与台站分布 Fig.1 The distribution of the seismic and blasting events and the stations |
本文选取曹妃甸测震台网在2014—2015年监测的10次M 0.8—4.2天然地震及10次典型爆破（朝鲜丰溪里核爆炸发生在2016年1月6日）进行时频分析。所选地震及爆破分布见图 1（2016年1月6日朝鲜丰溪里核爆炸因距离远未标注），统计数据见表 1。研究包括天津滨海化学爆破和朝鲜核试验在内的各类爆破时频图，可丰富研究的爆破类型，利于获取科学全面的研究结果。
|
表 1 天然地震与爆破事件统计 Tab.1 Statistics of earthquake and blasting events |
为保证研究的科学性和严谨性，对比分析的地震和爆破事件应发生在同一地点且震级相差不超过M 0.5，但实际几乎不存在此类事件，因此在震级相近的基础上，如果地震和爆破的震中到地震台站的距离小于二者到地震台站平均距离的30%，就认定二者发生在同一地点。此次选取的10次地震和10次爆破波形均在垂直向选取，截取的原始波形从P波开始至整个波形结束，即只截取全波。为排除噪声干扰，对地震和爆破事件的原始数据进行预处理，即计算而二者的幅频图，确定优势频率，通过带通滤波器，滤除干扰频率，保留有用信号。
3.2 时频分析时频分析结果好坏取决于窗函数类型。为获得清晰的时频分析图像，笔者试验了矩形窗、三角窗、海明窗和高斯窗，结果发现，用4种窗函数获得的时频分析结果本质一致，细节处稍有不同，因海明窗计算迅速，优点突出，本研究选用海明窗。原始数据采样率为50 Hz，设滑动步长50，窗口长度421（万永革，2007），利用Matlab软件编程，实现短时傅里叶变换，绘制地震波形的时幅图和时频图。统计10组地震与爆破的时频分析结果，并以2015年6月11日河北迁安ML 2.1矿爆炸、2015年8月12日天津滨海地区ML 2.9化学爆破、2015年9月14日河北昌黎M 4.1天然地震及2016年1月6日朝鲜丰溪里ML 4.9核爆炸为例，绘制地震波形的时幅图和时频图，见图 2。
|
图 2 典型地震及爆破时频分析 (a) 2015年6月11日河北迁安ML 2.1矿爆炸；(b)2015年8月12日天津滨海地区ML 2.9化学爆破；2015年9月14日河北昌黎M 4.1天然地震；(d)2016年1月6日朝鲜丰溪里ML 4.9核爆炸 Fig.2 Typical seismic and blasting time -frequency analysis |
由图 2可知:①地震事件频率分布在4—15 Hz范围内，其中6次地震频率范围在8—15 Hz，天然地震的瞬态谱散布在低频到高频的较大范围（韩绍卿等，2010）；爆破事件的频率范围在0—8 Hz，其中7次爆破的频率范围是2—4 Hz，且主频率集中在4 Hz；②地震事件有多个频率聚集点，且频率成分出现在S波段，而9次爆破的频率成分均出现在P波段，且只有一个频率聚集点。
图 2(d)是曹妃甸测震台网接收到的朝鲜核爆炸的时频分析图，从图中可见此次核爆炸的频率范围是0—0.5 Hz，频率成分只出现在S波段及后续波形中。正常情况下，地下核爆炸的力学机制是简单的球对称压缩，其震源时间函数表现为单脉冲形式且能量为高速释放的物理过程，源激发出各种频率的波，其中有的具有极高频率和极短波长（安镇文等，2011）。但曹妃甸测震台网接收的核爆炸频率较小，且P波段几乎无频率成分。出现该反常现象的原因在于，核爆炸的震中位置距曹妃甸测震台网较远，爆破源位于强度较低的介质中，爆破波比地震波穿透的地层浅，因而高频成分容易被吸收（刘芳彤等，2015）。由于传输介质密度低且传输距离大，核爆炸波在地表岩层传输过程中高频成分损耗，导致曹妃甸测震台网只接收到核爆炸波的低频成分。
4 结论通过对选取的10次地震和10次爆破波形进行时频分析，可以看出二者在时频域分布上存在明显差异，不同种类的爆破时频分布也不同。
（1）天然地震频率从低频到高频均有分布，选取的10次地震频率集中在几到十几Hz，无固定主频；爆破频率只有几Hz，主频一般出现在2—4 Hz。
（2）地震时频谱较爆破“分散”，而爆破的时频谱则相对“集中”（张帆等，2014）。地震时频图复杂且有多个频率聚集点，即出现“多峰”特征，而爆破时频图简单且只存在一个频率聚集点。由此可知，爆破时频图比天然地震更聚集。
（3）地震的频率成分小部分集中在P波，剩余大部分集中在S波，主要能量在S波阶段释放，与地震发生的实际情况相符。爆破的频率成分几乎全部集中在P波上，表明爆破能量集中在P波段释放。
震源性质决定了时频分布的特点，天然地震与爆破的震源不同，因此二者的时频分布特征存在差异。地震震源比爆破的震源深，地震波需要通过多层介质的折射和反射才能被接收，爆破的震源相对较浅，造成地震波和爆破波的传播路径不同。另一方面，地震的震源机制相对复杂，若干微小破裂组成地震震源，每个微小破裂均会反映在时频图上，爆破则是大膨胀后产生振动，爆破的震源过程近似于单脉冲函数，因此天然地震的频率散布在高频和低频中，且其时频图“多峰”复杂，而爆破的频率集中在低频段，其时频图“少峰”简单。
本文仅就爆破与地震的时频特征进行初步的定量研究，尚需大量观测资料进行深入探索。
| 安镇文, 魏富胜. 震源性质的时频分析与事件识别[J]. 地震地磁观测与研究, 2011, 32(1): 1-9. | |
| 崔鑫, 许力生, 许忠淮, 等. 小地震与人工爆破记录的时频分析[J]. 地震工程学报, 2016, 38(1): 71-78. | |
| 韩绍卿, 李夕海, 安跃文, 等. 核爆、化爆、地震识别研究综述[J]. 地球物理学进展, 2010, 25(4): 1206-1218. | |
| 刘芳彤, 杜春鸿. 朝鲜地区地震与爆破波形识别研究[J]. 防灾减灾学报, 2015, 31(2): 58-66. | |
| 王新刚, 郭琦南. 时频分析方法在探测信号处理中的应用[J]. 电脑知识与技术, 2013, 9(26): 5973-5975. | |
| 万永革. 数字信号处理的Matlab实现[M]. 北京: 科学出版社, 2007, 406-411. | |
| 张帆, 朱新运, 熊丹, 等. 基于非线性时频分析的地震和爆破识别[J]. 华南地震, 2014, 34(2): 56-63. | |
| 赵建明, 任佳, 符泽宇, 等. 曹妃甸地震台电磁扰动异常分析[J]. 地震地磁观测与研究, 2015, 36(4): 26-29. |