利用HHT方法浅析台站地震前地脉动信息变化特征

2018-04-27任建辉李继业武晓军孙强秦丽岩孙鹏宇张彦吉吴海波杜天娇

防灾科技学院学报 2018年1期

任建辉李继业武晓军孙强秦丽岩孙鹏宇张彦吉吴海波杜天娇

（1.哈尔滨市防震减灾技术中心，黑龙江哈尔滨 150021；2.黑龙江省地震局，黑龙江哈尔滨 150090；3.哈尔滨市通河地震台，黑龙江哈尔滨 150900；4.哈尔滨市阿城地震台，黑龙江哈尔滨 150300；5.哈尔滨市城乡建设委员会，黑龙江哈尔滨 150000；6.哈尔滨市木兰地震台，黑龙江哈尔滨 151900）

0 引言

从地震观测的角度可以把地脉动分为两类，一类是短周期地脉动，一类是长周期地脉动。地脉动也就是地震噪声，因为它干扰了小地震或远距离地震的观测，许多人对地脉动进行了研究。19世纪以来，微弱的震动，与地震、爆炸或当地的噪声源无关的地震波的持续振荡，已经在地震波记录中观测到了（bertelli，1872）。自 Omori（1908）用高倍测振仪器观测到地脉动，并发表“论地脉动”至今已有一百多年，国内外学者对地脉动的研究和应用已做了许多的工作［1］。早期的工作主要是对地脉动的成因、传播、定位等的研究，这些研究中主要是以长周期地脉动为对象，长周期的脉动源是如雨、海浪、地质内力作用等的自然因素（Douze， 1967； Capon， 1973）［2］； Klaus Strobach（1965）分析指出地脉动不仅包含了基本的瑞利面波，还有勒夫面波，两者在不同的台站所占有的比例是不同的；Astiz（1997）画出了全球台网的100个不同台站的典型噪声水平，发现脉动峰值在5～8s的周期出现，脉动峰值的源涉及在海洋产生的地震波［3］；Longuet-Higgins认为主脉动的峰值是由在开阔的海洋里沿不同方向传播的海洋波之间干涉产生的驻波所造成的［4］。

在临震预报方面，通过孕震介质转播的地脉动的频谱特征量的明显变化，为地震预报提供了有用的信息。牛志仁等分析海城、唐山地震前后地脉动的变化，发现震前同样的脉动幅值减小现象［5］；陈化然等利用功率谱和Sompi谱分析方法，对大震前的地脉动及特殊波形进行研究表明地脉动特征量可为地震预报提供有用的信息［6］；范叶萍等介绍了辽宁数字台网运行的脉动异常跟踪系统在2004年印尼苏门答腊西北近海8.7级地震前，营口、锦州台地脉动波形的幅度和周期明显异常情况，并展现了该系统的研究及应用前景［7］；之后，许多学者在长周期事件中发现了慢地震事件，并在一些大震前几天，发现了长周期形变波［8～9］；在此基础上，尹亮进行了进一步的震例研究，并从重力等角度进行了不同地球物理验证研究［10］；姚家俊则从缩短计算地脉动序列时间长度、频谱加窗处理等角度进行相关方法和实用技术的探索研究［11］；冯少孔等利用高密度面波勘探和小型台阵地脉动观测在西部地区地震防灾应用中发挥积极作用［12］。此外，还有许多科技人员利用地脉动做震源参数计算、地下速度结构反演等大量的研究工作［13］。

随着我国测震台网的建设，已经积累了巨量的大动态、高分辨率、高精度的波形资源，这些波形既包括地震事件波形，而更多的是无震的地脉动数据。黑龙江省数字地震监测台网数据从2002年产出，共有35个专业台站。2008年完成十五项目改造后，全部实现了数字化，目前测震台网包含省内及邻省共计50余个台站。本文旨在对数字地震观测台站积累的地脉动数据进行一次观测应用和预测应用的尝试研究，使数据资源发挥更大的作用。

1 研究区地质构造

依据地形地貌来分析，东北地块长期以来受到西北方向力和北东方向力的作用。西北方向的作用力主要来源于太平洋板块挤压扩张和抬升作用，形成了大兴安岭隆起和盆地，东北方向的力来源于北美板块，形成小兴安岭，将盆地分隔成结雅盆地和松嫩盆地。东北地区的地震可分为日本海―长白山深震带、郯庐断裂带北段地震区、松辽盆地内部地震带和松辽盆地西缘地震区［14］。松辽盆地是东北地区最大的中新生代陆相沉积盆地，盆地东、西两侧的依舒断裂带和嫩江断裂带的活动较弱，但盆地中的断裂仍有一定程度的活动。松辽盆地也是东北地区发生浅源中强地震的场所，地震主要分布在盆地的中部。考虑到数据获取和掌握情况，以及地震活跃程度和地域研究的必要性，研究区域主要集中在对黑龙江省防震减灾工作影响更具有意义的松辽盆地及其周边邻近地区。

2 资料选取

考虑到研究的统一性和未来的可持续应用性，选取“十五”项目改造之后的数据和震例作为研究对象。本次研究选用的数据是垂直分量UD向分钟值（在一般情况下，对于三个分向的数据选取没有特殊要求），地震选取分别为松辽盆地内部2009年5月10日安达地震、2013年10月31日前郭地震和松辽盆地附近地区2013年11月20日桦南地震。研究区内显著断裂展布及地震的空间分布见图1，地震具体参数信息见表1。

图1 研究区内主要断裂展布及3次M≥4.5地震空间分布图（F1.敦密断裂；F2.依舒断裂；F3.扶余-肇东断裂；F4.滨州断裂；F5.勃利-北安断裂；F6.嫩江断裂）Fig.1 Spatial distribution of fault zones and M≥4.5 earthquakes in study area

根据以往学者的研究成果，地脉动在短临预报方面有较好的应用，有成功的预测经历，具有数小时至7天左右的地震短临预报指标意义［15-16］。选取震前20天，震后10天的地脉动数据总结震例的短临预测指标，研究使用的数据均来自黑龙江省地震局监测中心，处理数据累计9264小时（约合555840分钟），涉及台站12个，其中黑龙江省台网10个，包括庆安台、北林台、碾子山台、五常台、肇源台、牡丹江台、双鸭山台、依兰台、密山台和佳木斯台；其他台网2个，为乾安台和榆树台。台站与震中分布均属同一构造块体、震中距250km范围内。

表1 研究地震信息表Tab.1 The information of earthquakes in study area

3 HHT方法原理

对于非平稳信号，如地震数据，其频谱是时间的函数，所以要求分析方法能够准确反映出信号的局部时变频谱特性，希尔伯特-黄变换（HHT）能够进行非线性、非平稳信号的线性化和平稳化处理，被认为是近年来对以傅立叶变换为基础的线性和稳态谱分析的一个重大突破［17］。它由经验模态分解（Empirical Mode Decomposition，简称EMD）方法和 Hilbert变换（Hilbert Transform）两部分组成，其核心是EMD分解。分解的过程中保留了数据本身的特性，再对各 IMF分量进行Hilbert变换，得到 Hilbert谱，算法流程如图2。Hilbert谱表示的是信号幅值在整个频率段上随时间和频率的变化规律，Hilbert边际谱表示信号幅值在整个频率段上随频率的变化情况，是通过对Hilbert谱积分得到的，Hilbert边际谱是具有一定的概率意义，可以看作是一种加权的联合幅值-频率-时间分布，而赋予每个时间-频率单元的权重即为局部幅值，而不是其信号上的幅值，只是通过其判断哪个具体时刻出现振动的可能性。

图2 HHT算法流程图（据徐克全，2011）Fig.2 HHT algorithm flowchart（based on Xu，2011）

在某些情况下，频域分析比时域分析更有优势。地震记录中蕴含着多种频率成分，就地震波本身来说，P波、S波、面波等具有不同的优势频率分布，而诸如环境变化、气象变化、爆破等干扰因素，也常常具有各自不同的频率特征。这些不同来源、不同特性的波动叠加在地震记录中，仅仅从地震波形特征的角度常常难以分辨；但是根据不同来源波动的频率差异，在频率域内进行分离，则可以较好地分辨出不同的影响因素及其变化。本文为了避免单独事件对数据的干扰，应用HHT方法边际谱对地脉动的每日分钟值数据进行分析，找出每日频率叠加幅值最大时段，这样可以有效地分离了单个事件（前震、爆破、余震、其它干扰等）频率对整体数据分析的影响。计算中，采用Matlab提供的HHT算法对分钟地震记录进行边际谱计算。

4 数据分析

利用HHT方法对地脉动UD分量每日分钟值进行频谱信号计算与提取（图3）。

通过对3次地震前，多台站、多时段的边际谱分析发现，边际谱在20～49μHz脉动的幅值出现快速、持续地增大，而此前的相应结果比较稳定，仅仅在一定的范围内波动。下面进行具体分析。

安达地震研究选取北林、碾子山和庆安3个台站（图4a），数据截取2009年4月20日至5月20日；前郭地震研究选取乾安、五常、榆树和肇源4个台站（图4b），数据截取2013年10月11日至11月10日；桦南地震研究选取牡丹江、双鸭山、依兰和密山4个台站（图4c），数据截取2013年10月31日至11月30日。分别利用HHT方法分别对每台、每日UD分量进行边际谱分析，具体结果见表2。

续表2

图3 HHT方法计算地脉动UD分量24小时分钟值Fig.3 24 hours minute value of microtremor UD component by HHT method

图4 选取台站与地震分布图Fig.4 Distribution map of selected seismic stations and earthquakes

图5为安达地震地脉动边际普幅值观测序列。其中，北林台、庆安台在4月21日出现幅度值随时间单调增大的现象，碾子山台则在4月22日出现类似情况，北林台在5月2日后恢复到平静值，庆安台在5月3日恢复，碾子山台则在6日恢复。距离震中最近的北林台持续时间为12天，庆安台持续13天，最远的碾子山台持续了15天，距离越远持续的时间相对越长；从出现幅值增大的时间来看，距离震中位置相对较近的北林台（140km）和庆安台（180km）最先出现，且最大幅度值也较大（分别为15.2×106A、29.6×106A，A表示幅值，单位counts，下同），距离震中200km的碾子山台则出现的较晚，幅度值相对也较小。优势频率集中在20～38μHz，其中优势幅值频率为26μHz。本文又分别对前郭地震和桦南地震的地脉动边际普幅值序列进行了时间扫描，由于篇幅的关系，仅把扫描结果进行说明，具体情况见表3。

表3显示，前郭地震前，乾安台、五常台和榆树台在10月12日出现幅度值随时间单调增大的现象，肇源台则在10月14日出现类似情况，榆树台在10月26日后恢复到平静值，乾安台、五常台在5月28日恢复，肇源台则在31日地震当日恢复。距离震中最近的乾安台持续时间为17天，肇源台持续18天，相对较远的榆树台和五常台则分别持续了15天、17天，距离越远持续的时间相对越短；从幅值增大的时间来看，距离震中位置相对较近的乾安台（50km）最先出现，且最大幅度值也较大（19.0×106A），肇源台（110km）虽然出现幅值增大的时间相对较晚，但最大幅值大于距离震中最远的五常台。该震例中存在一个特殊情况就是，最大幅值出现在榆树台（230km），达到了71.1×106A，分析认为，这可能与前郭地震类型有关。众所周知，前郭地震是震群型，最大震级达到5.8级，且时间跨度近一个月。此外，根据作者以往的研究，一个中强地震可能与一个或者几个典型场有关，一个典型场也可能与一个或几个中强地震对应［18］。10月31日地震后，仍有部分台站出现幅值增大的现象，可能与震兆场对应多个地震现象有关。前郭地震的优势频率集中在20～49μHz，其中优势幅值频率为32μHz。

图5 安达地震地脉动UD向日值边际谱幅值观测序列Fig.5 Observation sequence of marginal spectrum amplitude of microtremor UD component of Anda earthquake

表3 边际谱幅值信息与地震的关系Tab.3 The relationship between the marginal spectrum amplitude information and earthquakes

桦南地震前，依兰台、密山台在11月6日出现幅度值随时间单调增大的现象，双鸭山台在11月7日随即出现增大现象，牡丹江台则在11月9日出现类似情况，双鸭山台在11月19日后恢复到平静值，牡丹江台、依兰台、密山台在11月20日同时恢复。距离震中最近的双鸭山台持续时间为13天，依兰台、密山台各持续15天，相对较远的牡丹江台则持续了12天，距离越远持续的时间相对越短；从出现幅值增大的时间来看，距离震中位置相对较近的依兰台（103km）和密山台（113km）最先出现，且最大幅度值也相对较大（36.8×106A和33.0×106A），然后是距离震中最近的双鸭山台（76km），幅度值达到了41.0×106A，距离最远的牡丹江台（180km）出现幅值增大的时间最晚，最大幅值也最小。桦南地震的优势频率集中在20～49μHz，其中优势幅值频率为38μHz。

综上所述，地震发生前，一般表现在15天左右，连续出现5次或以上频率幅度高值大于5.00×106A的现象。为了考证这种现象是否与地震前应力场变化引起的脉动异常有一定关联性，本文还对佳木斯、牡丹江和依兰台在无震平静时段的地脉动UD分量边际普幅值观测序列进行时间扫描，15天内未发现5次以上高于5.00×106A的频率幅度值，且频率幅度都相对偏小（图6）。充的作用。由于篇幅的关系，详细结果见表4，图7为前郭地震台站间地脉动互相关时间扫描情况。

图6 平静时段地脉动UD向日值观测序列边际谱（2016.4.1～4.15）Fig.6 Marginal spectrum of daily value observation sequence of microtremor UD component in quiet period

表4 地震与异常时间的关系表Tab.4 The relationship between the abnormal time and earthquakes

图7 前郭地震地脉动UD向整点值互相关曲线Fig.7 Cross-correlation curve of whole value of microtremor UD component of Qianguo earthquake

表4显示，安达地震前，互相关曲线在4月21日出现明显离散现象，在26日达到最大，5月2日逐渐恢复，恢复一周后发生地震。前郭地震前，肇源／乾安互相关曲线在10月12日出现明显离散现象，在16日达到最大，10月27日逐渐恢复，恢复一周内发生地震；肇源／榆树互相关曲线在10月16日出现离散现象，在24日达到最大，10月25日逐渐恢复，恢复一周内发生地震；五常／乾安互相关曲线在10月12日出现离散现象，在16日达到最大，10月27日逐渐恢复，恢复一周内发生地震；榆树／乾安互相关曲线在10月12日出现明显离散现象，在13日达到最大，10月30日恢复后一天发生地震。桦南地震前，密山／双鸭山互相关曲线在11月7日出现离散现象，在18日达到最大，10月20日地震发生后逐渐恢复；密山／依兰互相关曲线在11月6日出现明显离散现象，在8日达到最大，11月19日恢复后一天发生地震；双鸭山／依兰互相关曲线在11月6日出现明显离散现象，在7日达到最大，11月18日逐渐恢复，恢复一周内发生地震。

5 结论与讨论

本文通过HHT方法计算中强地震前地脉动数据的边际谱，总结出地震前边际谱的优势频率范围，发现地脉动幅值在地震前随时间快速、持续地增大，并在震前恢复平静的完整过程。研究的初步结论如下：

（1）地震发生前，频率幅值存在随时间上升—下降—恢复的完整动力过程变化，频率集中在 20 ～ 49μHz，其中优势幅值频率为 26μHz、32μHz和 38μHz。

（2）地震前，幅值增高变化在时域上一般表现为15天左右，且连续出现5次或以上频率幅度高值大于5.00×106A的现象。

（3）频率幅值与震级存在一定的正比例关系，震级越大，幅值相对越大。震中距与幅值存在一定的反比例关系，台站距离震中越远，出现幅值相对越小。异常出现时间与震中距存在一定的正比例关系，距离震中越近的台站，可能越早出现异常。

（4）对出现频率幅度高值的台站进行互相关分析，发现相关性曲线存在明显的离散现象，且这种现象在地震发生前逐渐减弱或消失。

（5）地震一般在异常消失后的一周内发生，具有短临预测意义。

通过研究也发现该项工作也存在着一些缺点和需要继续完善的方面：

（1）算法复杂且无法快速的进行计算。数据处理量较大，研究时段长度受制约，适合短临预测以及追溯性检验、总结工作。如将该方法应用到日常工作当中，每日处理数据，可适用于长期震情跟踪工作。

（2）在处理干扰过程中由于仪器不同、工作周期不同、台站受到的干扰因素不同，一些干扰只能有针对的定向去分析，对于干扰严重的数据，可能无法完全消除干扰信息，只是相对提高数据质量。

（3）受到东北地区地震震例少，震级小的限制，在分析提取地震前兆异常中，仅对3个震例进行了分析研究，且台站选取距离震中相对较近。未来需要利用该方法继续开展同构造块体内远场信号的提取与识别工作，尝试降低震级下限探索预测指标的适用性，并逐渐完善平静期背景指标，丰富震例总结数量，拓宽检验时间长度等工作。考虑到地震活动受到许多因素的影响，如太阳黑子活动、地球自转速率变化以及活动断裂运动等等，本次研究并未予以考虑，因而研究结果具有局限性［19］。此外，本文仅对3次震例进行了现象描述和较浅的特征分析，预测指标也需要日后更多的震例研究进行检验与完善，并结合其它方法进行辅助预测。

［1］董蕾.大震前地脉动位移及其频谱特征研究［D］.兰州：中国地震局兰州地震研究所，2010.

［2］ Robert K.Cessaro.Sources of Primary and Secondary Micmseisms［J］.Bulletin of the Seismological Society of America，1994，84（1）：142-148.

［3］ Peter M.Shearer.地震学引论［M］陈章立译.北京：地震出版社，2008.

［4］ Longuet-Higgins M S.A theory of the origin of microseisms［J］.Phil.Trans.R.Soc，1950，London.A 243：1-35.

［5］牛志仁，朱传镇.脉动与地震关系的初步讨论（二）［J］.地球物理学报，1978，21（4）：325-331.

［6］陈化然，郭瑞芝，冯德益.地脉动及特殊波形的频谱特性在地震预报中的应用［J］.地震，1998，18（1）：75-82.

［7］范叶平.数字地震台网大震前地脉动信息变化特征研究［J］.国际地震动态，2004，305 （5）：13-18.

［8］王迪晋，李正嫒，吕品姬.慢地震研究综述［J］.大地测量与地球动力学，2007（27）：21-25.

［9］杨又陵，赵根模，高国英，等.2001年11月14日昆仑山口西M8.1地震前的缓慢地震事件［J］.国际地震动态，2003，297（9）：1-4.

［10］尹亮.宽频带地脉动频谱异常震例分析及多角度地球物理观测的验证研究［D］.兰州：中国地震局兰州地震研究所，2010.

［11］姚家俊.数字化宽频带地脉动资料处理技术和频谱特征研究［D］.中国地震局兰州地震研究所，2010.

［12］冯少孔，车爱兰，吴志坚，等.利用高密度面波勘探和小型台阵地脉动观测在西部地区地震防灾应用［J］. 地震工程学报，2013，35 （4）：780-788.