2019年3月—8月,河南平顶山市平煤矿区附近小震持续活动,河南测震台网共记录到21次震动事件,最大为2019年3月15日河南宝丰ML 2.9地震。由于矿区复杂的构造环境、井下爆破施工和采空区塌陷等因素影响,事件类型成为地震危险性研判的关键依据,判定方法主要有:时间域初动方向、振幅和周期测量(赵永等,1995),频率域频谱特征(沈萍等,1999;黄汉明等,2010),倒谱域判据(魏富胜等,2003;郭祥云等,2010)以及时域、频域信息与人工神经网络结合(Reynen et al,2017;隗永刚等,2019)等。
震源机制不同,信号频段组成有所差别(张萍等,2005),而时频分析方法能够描述信号频率随时间的变化特征,是刻画地震记录等非平稳信号在时间和频率轴上能量强度分布的有效手段(郑建常等,2014)。近年来,广泛应用于地震波初至检测(詹毅等,2004;邹文等,2006;刘希康等,2016)、地震资料分辨率提高(高静怀等,1996;高静怀,1997;柳建新等,2006)、地震波能量补偿(焦叙明,2007)、窗函数选择(马长征等,1999;章兰英等,2005;肖瑛等,2010)、事件类型识别(张帆等,2014;崔鑫等,2016;荣伟健等,2017;李丽等,2018)等。
本文在传统时域初动方向、振幅和周期等波形特征判定事件类型的基础上,采用短时傅里叶变换(STFT)时频分析方法,对2019年3月—8月平煤矿区发生的ML 2.0—2.9天然地震、爆破、塌陷等9次事件,进行时频波谱对比,梳理不同类型震动事件识别的有效定量指标,为研究区事件类型判断提供参考依据。
1 构造背景和地震活动平顶山煤田是华北晚古生代聚煤盆地的一部分,处于豫西断隆、华北断拗和北秦岭褶皱带的衔接部位。喜山运动使前期形成的凹陷进一步发展,接受厚度达千米以上的新生代沉积,形成目前煤田的地质景观。区域构造形态为一系列NW向褶皱,以李口向斜和灵武山向斜为主,为东南端收敛西北端呈扇型展开的宽缓褶曲,向NW倾伏(中国科学院地质研究所,1980)。区内断层主要是1组平行褶曲轴向的正断层及逆断层,其中凤凰岭逆断层位于整个矿区西部,走向近NS,倾向西,破坏了煤层向西的延伸;锅底山—小山正断层为研究区中部较大的断裂构造,呈NW—SE向展布,与李口向斜大致平行,倾向SW,为高角度阻水正断层(王志铄,2017)。矿区10余座矿井主要分布在平顶山市区以北与宝丰县交界地区(图 1)。
平顶山地区深部结构复杂,湖北随州至安阳地学断面揭示该地区存在深部地球物理变异带,深地震测深结果显示,断裂主要在上地壳展布,20 km以下的中下地壳存在高速异常体(胡鸿翔等,1986)。该地区存在发生构造地震的深部背景和孕震环境,但断裂构造以早—中更新世为主,历史记录中未发生过5级以上地震;2008年地震监测台网数字化改造以来仅记录到ML 3.0以下地震事件45次,部分地震沿区域内主要断裂分布(图 1)。
2 方法和数据采用短时傅里叶变换(short-time Fourier transform,STFT)(Gabor,1946)方法,研究平煤矿区不同类型事件的时频特征。STFT是一种线性时频分析方法,基本思想是,用窗函数来截取信号,并假定信号在窗内是平稳的,采用傅里叶变换分析移动窗口内信号,得到信号频率分布随时间的变化关系。
时间序列x(n)的短时傅里叶变换表示为
$ {\rm{STFT}}x(n, \omega) = \int\limits_{ - \infty }^\infty {x(m)\omega (n - m){e^{ - j\omega m}}{\rm{d}}} m $ |
或
$ {\rm{STFT}}x(n, \omega) = \sum\limits_{ - \infty }^\infty {x(m)\omega (n - m){e^{ - j\omega m}}} $ |
式中:ω(n-m)表示时间窗函数,n表示时间点,ω表示频率点,m为变换参数。
STFT对计算机性能要求不高,具有计算速度快、分辨率较好等优点(崔鑫等,2016;荣伟健等,2017),且对非平稳信号无交叉干扰项(刘希康等,2016)。因此选择STFT方法开展时频计算。
选取2019年3月—8月平煤矿区发生的ML 2.0—2.9天然地震、爆破、塌陷等9次震动事件作为研究资料,其中天然地震6次,爆破2次,塌陷1次,地震基本参数见表 1。选取其中震级相当的3次不同类型的震动事件,利用距平煤矿区最近的4个地震台站(台站代码PDS、LUS、JS、ZMD)记录,分析天然地震、爆破、塌陷的初动方向、周期、振幅比等特征参数,结果见表 2。
统计结果显示,天然地震垂直向初动一般按四象限分布,P波周期略小于S波,振幅比APg/ASg<0.5,近台面波明显但远台无面波;与天然地震相比,爆破垂直向初动均向上而塌陷均向下,非天然地震振幅比APg/ASg>1.0,面波发育。这些波形特征是事件类型判断的主要依据。
截取距平煤矿区最近的平顶山(PDS)台(距矿区中心约12 km)垂直分量记录初至波到时前5 s、长度30 s的完整波形作为样本,采用最大值归一化方法预处理,经傅里叶变换计算幅频图,确定以优势频率为0.5—20 Hz的带通滤波,过滤掉地震信号中的环境噪声和低频干扰。
3 窗函数及窗长选取万永革(2012)研究发现,短时傅里叶变换引起的频谱泄露与窗函数频谱的两侧旁瓣有关,若两侧旁瓣高度趋于零,而能量相对集中在主瓣,则认为接近真实频谱。矩形窗、汉宁窗、海明窗3种窗函数的时域和幅频形状显示(图 2,图中n表示离散化的时间序列;w(n)为返回的窗函数序列),各窗函数均有明显的主瓣和旁瓣,其中矩形窗的主瓣最窄,旁瓣峰值最大;海明窗和汉宁窗主瓣具有较小的旁瓣和较大的衰减速度,而海明窗第一旁瓣相对主瓣衰减更快。因此,选择海明窗作为窗函数进行时频计算。
短时傅里叶变换通过滑动窗来计算频谱,短窗口时间分辨率较高,但频率分辨能力差;长窗口频率分辨率较高,但时间分辨能力弱(陈雨红等,2006;万永革, 2007, 2012;董建华等,2007)。原始数据采样率为100 Hz,即每秒100个采样点,分别以51、171、391个数据为海明窗窗口长度进行分辨率测试,结果见图 3,可知:当窗长(采样率为100 Hz时离散数据点个数)为51时,时间分辨率较高,频率分辨率较差;当窗长为391时,频率分辨率较好,时间分辨率较差;当窗长为171时,频率和时间分辨率均较为适中。因此,时频计算时采用海明窗窗口长度为171。
采用STFT时频分析方法,计算表 1所列9次震动事件时频分布。为了便于提取量化识别指标,将滤波后的垂直向记录振幅与时频图能量进行归一化处理。
前人对不同类型事件的时频特征研究发现,天然地震由于震源深、衰减小,频率成分较为丰富(沈萍等,1999),一般均匀分布在几到十几Hz之间(张帆等,2014;崔鑫等,2016),S波能量较大且存在“多峰”特点(荣伟建等,2017;李丽等,2018);爆破和塌陷相对天然地震发震机制比较简单,由于震源浅,高频成分多被松散沉积层吸收,短周期面波发育(张萍等,2005),主频较低、频率成分较单一,相对“少峰”(崔鑫等,2016),时频分布集中性明显高于天然地震(张帆等,2014);塌陷事件的震源一般在地表以下几百米范围内,相对于爆破,塌陷事件初至波时频谱无明显峰值,面波能量更为突出(李丽等,2018)。
4.1 天然地震波谱特征图 4给出6次天然地震垂直向波形和时频分布。在波形记录上,P波、S波、面波清晰并依次出现,P波振幅小于S波,面波周期最大。时频计算结果显示,天然地震频率成分丰富,在10 Hz以内高、低频分布均匀,无明显主频;P波部分在约3 Hz和8 Hz处存在2个峰值;S波部分能量最大,存在多个峰值且随时间衰减较快;面波频率最低,多分布在3 Hz以内,随时间衰减较慢。
图 5给出2次爆破事件垂直向波形和时频分布,可见垂直向波形初动均向上,P波振幅最大、周期最小,面波发育且周期最大,S波淹没在面波中不易识别,且较天然地震衰减快。爆破事件时频谱相对集中,以低频为主,高于5 Hz的频率成分能量较小;初至波频率最大,且峰值约5 Hz,信号主频随时间大致呈线性变化,逐渐降低至1—2 Hz;面波频率最低;能量主要集中在P波和面波上。
塌陷事件记录波形垂直向初动向下,频率成分较为简单,面波发育且周期、振幅最大。图 6给出塌陷事件垂直向波形和时频图,可见:频率成分以4 Hz以下的低频为主;P波部分无峰值,频率成分单一;面波部分频率分布较均匀,主频出现在约2 Hz处,低频面波发育且频率随时间平缓衰减。
平煤矿区天然地震、爆破、塌陷事件波形特征和时频特征同震源性质、地质构造、深部背景相关。区域断裂主要在上地壳展布,天然地震可用双力偶模型解释初动四象限分布等波形特征,良好的传播路径使震相丰富、频率成分分布均匀;平煤矿区人工爆破主要集中在约600 m深的第三开采层位,膨胀源性质决定垂直向记录初动向上、初至波能量大,高频成分多被浅表岩层吸收,故波形记录以低频为主;平煤矿区采空区处理以自然垮塌为主,塌陷岩体向底部冲击,故近台地震波初动向下,频率成分简单且衰减较快,大周期面波发育。
时频计算结果显示,不同类型震动事件的时频特征具有明显差异,事件类型时频判断与经验性判断结果一致:天然地震频率成分均匀分布在10 Hz以内,P波波形在约3 Hz和8 Hz处存在2个峰值,S波波形的频率成分存在多个峰值,且主频随时间推移逐渐降低;爆破事件频率成分较为集中,一般低于5 Hz,初至波能量最大,频率成分峰值约5 Hz,信号主频随时间大致呈线性变化,快速衰减至1—2 Hz;塌陷事件频率成分以4 Hz以下的低频为主,初至波无明显峰值,主频出现约2 Hz处,低频面波发育,且频率随时间平缓线性衰减。时频分析方法判断事件类型,可为区域震情趋势研判、矿区非天然地震监测等提供参考依据。
陈雨红, 杨长春, 曹齐放, 等. 2006. 几种时频分析方法比较[J]. 地球物理学进展, 21(4): 1180-1185. DOI:10.3969/j.issn.1004-2903.2006.04.019 |
崔鑫, 许力生, 许忠淮, 等. 2016. 小地震与人工爆破记录的时频分析[J]. 地震工程学报, 38(1): 71-78. DOI:10.3969/j.issn.1000-0844.2016.01.0071 |
董建华, 顾汉明, 张星. 2007. 几种时频分析方法的比较及应用[J]. 工程地球物理学报, 4(4): 312-316. DOI:10.3969/j.issn.1672-7940.2007.04.009 |
高静怀, 汪文秉, 朱光明, 等. 1996. 地震资料处理中小波函数的选取研究[J]. 地球物理学报, 39(3): 392-400. DOI:10.3321/j.issn:0001-5733.1996.03.013 |
高静怀. 小波变换用于高分辨率地震资料处理及瞬时特征分析的理论与方法研究[D]. 西安:西安交通大学,1997.
|
郭祥云, 王淑辉, 魏富胜. 2010. 小震级地震事件的倒谱差异[J]. 地震地磁观测与研究, 31(1): 12-16. DOI:10.3969/j.issn.1003-3246.2010.01.003 |
胡鸿翔, 陈学波, 张碧秀, 等. 1986. 中国中原地区随县-安阳剖面深地震测深资料的解释[J]. 地震学报, 8(1): 37-49. |
黄汉明, 边银菊, 卢世军, 等. 2010. 天然地震与人工爆破的波形小波特征研究[J]. 地震学报, 32(3): 270-276. DOI:10.3969/j.issn.0253-3782.2010.03.002 |
焦叙明. 时频分析及其在地震资料处理分析中的应用[D]. 青岛:中国海洋大学,2007.
|
李丽, 刘剑, 李自红, 等. 2018. 山西运城振动事件S变换时频分析[J]. 地震地磁观测与研究, 39(2): 77-83. DOI:10.3969/j.issn.1003-3246.2018.02.011 |
柳建新, 韩世礼, 马捷. 2006. 小波分析在地震资料去噪中的应用[J]. 地球物理学进展, 21(2): 541-545. DOI:10.3969/j.issn.1004-2903.2006.02.032 |
刘希康, 丁志峰, 李媛, 等. 2016. STFT与WVD在地震波时频分析中的应用[J]. 地震地磁观测与研究, 37(1): 15-21. |
马长征, 张守宏, 焦李成. 1999. 短时傅里叶变换和拟Wigner分布最佳窗函数[J]. 电子科学学刊, 21(4): 467-472. |
荣伟健, 赵建明, 董祎玮, 等. 2017. 曹妃甸地震台网天然地震与人工爆破时频分析[J]. 地震地磁观测与研究, 38(5): 39-43. DOI:10.3969/j.issn.1003-3246.2017.05.007 |
沈萍, 郑治真. 1999. 瞬态谱在地震与核爆识别中的应用[J]. 地球物理学报, 42(2): 233-240. DOI:10.3321/j.issn:0001-5733.1999.02.011 |
万永革. 2007. 数字信号处理的Matlab实现[M]. 北京: 科学出版社, 188-201.
|
万永革. 2012. 数字信号处理的Matlab实现[M]. 2版. 北京: 科学出版社, 384-395.
|
王志铄. 2017. 河南省地震构造特征[M]. 北京: 地震出版社, 1-23.
|
魏富胜, 黎明. 2003. 震源性质的倒谱分析[J]. 地震学报, 25(1): 47-54. DOI:10.3321/j.issn:0253-3782.2003.01.006 |
隗永刚, 杨千里, 王婷婷, 等. 2019. 基于深度学习残差网络模型的地震和爆破识别[J]. 地震学报, 41(5): 646-657. |
肖瑛, 冯长建. 2010. 组合窗函数的短时傅里叶变换时频表示方法[J]. 探测与控制学报, 32(3): 43-47. DOI:10.3969/j.issn.1008-1194.2010.03.011 |
詹毅, 钟本善. 2004. 利用小波变换提高地震波初至拾取的精确度[J]. 成都理工大学学报(自然科学版), 31(6): 703-707. DOI:10.3969/j.issn.1671-9727.2004.06.027 |
张帆, 朱新运, 熊丹, 等. 2014. 基于非线性时频分析的地震和爆破识别[J]. 华南地震, 34(2): 56-63. |
章兰英, 侯孝民, 郑海昕. 2005. 短时傅里叶变换软件解调中窗函数影响分析[J]. 装备指挥技术学院学报, 16(6): 98-101. DOI:10.3783/j.issn.1673-0127.2005.06.023 |
张萍, 蒋秀琴, 苗春兰, 等. 2005. 爆破、矿震与地震的波谱差异[J]. 地震地磁观测与研究, 26(3): 24-34. DOI:10.3969/j.issn.1003-3246.2005.03.004 |
赵永, 刘卫红, 高艳玲. 1995. 北京地区地震、爆破和矿震的记录图识别[J]. 地震地磁观测与研究, 16(4): 48-54. |
郑建常, 徐长朋, 赵金花, 等. 2014. 基于S变换的矿山地区地震波形时频分析[J]. 地震地磁观测与研究, 35(5/6): 8-14. |
中国科学院地质研究所. 1980. 华北断块区的形成与发展[M]. 北京: 科学出版社.
|
邹文, 陈爱萍, 顾汉明, 等. 2006. S-变换初至拾取技术及其在复杂山区中的应用[J]. 石油天然气学报(江汉石油学院学报), 28(2): 63-65. |
Gabor D. 1946. Theory of communication[J]. J Inst Electr Eng, 93: 429-457. |
Reynen A, Audet P. 2017. Supervised machine learning on a network scale:Application to seismic event classification and detection[J]. Geophys J Int, 210(3): 1394-1409. |