基于STA/LTA算法的地脉动有效信号自动识别研究

朱 胜，李 平，李玉影，刘应慈，周 楷

（1.防灾科技学院地质工程学院，河北三河 065201；2.河北省地震灾害防御与风险评价重点实验室，河北三河 065201）

引言

地脉动是地面上一种平稳的非重复性随机波动，主要是由风、海浪活动等自然振源或场地周围的机械振动、交通工具等人工振源激发，经过在不同性质的岩土层界多次反射和折射，叠加后传播到地面的振动。地脉动具有频率低、振幅小等特点，其频率范围一般为0.1～10 Hz，而振幅一般在0.1～1μm。地脉动的不同振幅和频率的变化，既能反映测试场地的土层特征，又能反映工程场地的动力特性。早在1907年，学者Omori［1］观测到了一种振幅很小的振动，他将这种微小的振动称之为地动脉（Microtremor），并发表了《论地脉动》一文，在文中初步描述了地脉动的振动特征，此后日本和美国的学者对地脉动这一热点课题进行了大量的研究。但是很多学者对地脉动产生的机理不明确，所以地脉动的研究一度陷入沉寂，直到19世纪60年代日本学者Kiyoshi等［2］提出将地脉动测试应用于场地评价中去，日本学者Nakamura提出用单点谱比法处理地脉动数据，地脉动研究又重新引起国内外学者的兴趣，随后美国、欧盟、中国等各国学者也相继开展大量的地脉动测试研究，促进了地脉动的在工程中的应用和发展；但是Lermo等［3］一直认为地脉动产生的机理不明确，对地脉动用于场地评价仍一直持有否定的态度；地脉动的研究至今已有百年的历史，我国学者对地脉动的研究起步于20世纪60年代，郭明珠等［4］，师黎静等［5］，刘红彪等［6］，分别地脉动测试应用有在震害分析、建筑结构设计、场地类别划分和地基场地处理效果评价等工作中取得了一些的工程应用成果；师黎静等［7］通过选取46处场地的地脉动数据，采用单点谱比法计算出场地卓越频率，探究地脉动测试法应用于场地特征参数测试的适用性，并且给出了覆盖层厚度、不同深度等效剪切波速以及特征周期与卓越频率之间的拟合曲线，验证其相关性；王伟君等［8］利用在北京五棵松地区进行地脉动测试，采用单点谱比法分析场地卓越频率，台阵地脉动测试F-K法，反演土层的速度结构。研究成果表明：地脉动能够反映出场地土层的构造信息，相比于目前的场地动力特性测试方法，具有经济、便捷和高效的特点，有广阔的工程应用前景。就目前而言，学者们已经对地脉动测试的现场采集方法达成了共识，但是在后期数据的处理仍存在争议。鉴于此，文中以松原5.7级地震砂土液化场地的地脉动测试数据为研究对象，施加反STA/LTA算法，探究数据处理中有效信号时间窗的选取问题，为以后规范化处理数据提供依据。

1 地脉动数据来源

2018年5月28日1点50分，在吉林松原市宁江区发生了5.7级地震，震源深度为13 km。地震发生后，防灾科技学院迅速组织人员，对震中受灾区域开展震害调查工作。选取了26处典型的场地进行地脉动测试数据采集，其中液化场地为19处，未发生液化场地7处，先后2年时间一共采集108条地脉动数据记录。

本次测试设备采用美国Kinemet-rics公司生产的Basalt型一体式强震仪，其内部设置三通道加速度计，能够同时记录2个水平方向和垂直方向的加速度，仪器的动态范围为155 dB。仪器布置时利用地质罗盘定位，按南北、东西和垂直3个正交方向放置。采样率设为200 sps，电压灵敏度为2.5 V/g，每个测点每次测试连续采样300 s，每个测点数据采样3次［19］，地脉动的采集结果如图1所示。

图1 测点M1地脉动三分量记录结果Fig.1 Microtremors 3-components time-histories at site M1

地脉动测试是一种在无强震源的情况下，仅仅是依靠自然震源或人为震源产生的地表微弱的振动，需要使用灵敏度较高的加速度传感器进行监测得到。因为振动幅值比较小，所以在数据采集的过程中，需在观测环境中在一定范围内无特定的振动源，如汽车、机械运转的影响等。地脉动的测点一般选择自由场地，并且远离道路。选择的数据采集的时间一般选择晚上，但是考虑到实际的测试情况进而可以考虑选择白天进行地脉动测试，从早上8点到下午6点，进行地脉动数据观测。数据采集时一般在风和日丽的天气进行。刮风下雨天气对测试结果也有一定的影响，下雨之后，雨水渗透作用导致覆盖土层抗剪强度下降，对高频波有抑制作用，对低频波有放大效应。

本次地脉动的测试采集应为是地震应急科考期间时间比较急迫，所以选择白天开展地脉动测试。在地脉动数据处理的过程中，发现了地脉动数据有一定的“毛刺”和脉冲信号干扰，采用带通滤波仍然无法有效的去除，如图2所示。因此，文中提出对地脉动信号施加反STA/LTA算法，提取地脉动的有效信号段，剔除“毛刺”和脉冲信号段以达到滤波的效果。

图2 测点M 13地脉动三分量记录结果Fig.2 Microtremors 3-components time-histories at site M13

2 STA/LTA法简介

震相识别是地震学研究中一个重要课题分支，它为地球内部构造的研究、振源的定位、震源机制的推演、地震速报和预警等一系列研究工作提供了基础性的研究资料。根据识别的方法，可将震相识别划分为人工识别和自动识别。传统的人工识别方法主要有：直观检验法、走时表检验法、综合检验法和合成地震图检验法等等，这些方法都是基于识别者的工作经验来加以判断、辨识，主要的缺点在于判别效率低，并且由于人为的误差、读数错误，可能导致不同的人识别的精度相差较大［9-10］。随着2008年以后我国数字化地震台网的大规模建立，地震观测数据量巨大，单单依靠人工识别是不切实际的，因此应运而生了稳定、快速和高效的自动识别方法。目前，自动识别方法主要分为：时域分析法、频域分析法、时频分析法、模式识别法和综合分析法［9-10］。文中着重介绍时域分析法中的STA/LTA法。

STA/LTA法是由Stevenson［11］最早提出，后经过Allen等［12］的不断改进，用来检测地震事件和拾取震相。其原理是在原信号上施加两个时间窗，分别为短时窗和长时窗，用短时窗的平均值（short-term average，简称STA）和长时窗的时长（long-term average，简称LTA）的比值（R）来反映信号能量的变化程度，当地震信号到达时，STA值比LTA值变化的快，STA值快速变大，R值会有一个显著的增大趋势，当R值超过设定的阈值（THR）时，则此点被判断为初动，地震记录系统开始记录有效信号。

式中：X(i),(i=1,2,3···,N)，表示短时间窗内记录数据；Y(j),(j=1,2,3···,M)，表示长时窗内记录数据，M和N分别表示长短窗内的样本个数。

STA/LTA方法计算量小、运算程序简单，但是在实际应用过程中，其计算结果常常出现计算偏差，导致误判的现象。为了提高信号识别的准确率，增加信号到达时的灵敏度，所以学者Allen［10］提出基于原始的地震记录，重新构造能够灵敏反应信号时间序列CF（t），并将新的时间序列作为计算的输入参数。新的时间序列C F（t）被称为特征函数。常用的特征函数有3种：

随着研究的不断深入，很多学者提出改进的特征函数形式，国内学者武东坡和高淑芳［9，17］利用三角函数推导，提出了一种改进的STA/LTA震相自动识别的特征函数：

刘晗通过大量微震监测信号分析发现，STA/LTA算法特征函数取式（3）和式（5），对微震记录信号识别响应比较灵敏，并且式（5）对频率变化也比较敏感［13］。根据短时窗与长时窗的关系，可以将STA/LTA计算方法分为标准STA/LTA法和延迟STA/LTA法。标准STA/LTA法是短时窗与长时窗之间没有重叠或延迟，进行截取计算；而延迟STA/LTA法，顾名思义即在短时窗和长时窗之间有一段延迟。实际应用中，为了提高计算效率，在标准法的基础上，逐渐发展了递归STA/LTA法。递归法相比于标准法计算更加快速，拾取信号的灵敏度更高，对脉冲信号有一定的抑制作用，减小误触的风险。

递归法的计算公式如下：

式中：STAi和LTAi分别为信号在i时刻短时窗与长时窗的平均值；C F（i）为信号在i时刻特征函数值；Nsta和Nlta分别为信号短时窗和长时窗的信号记录点数。

3 基于地脉动数据有效信号自动识别

目前随着各个学科之间不断交融与发展，STA/LTA算法的应用不仅仅局限于震相识别领域的，在其他信号“初动”识别领域也被广泛应用，如煤矿冲击地压的预测、煤与瓦斯突出的预防、深部矿山评估岩体稳定性、预测岩爆等动力型破坏等方面［13-14］。

地脉动是一种平稳、随机波动，其最大的特点就是振幅小、频率低。地脉动测试受环境影响较大，因此测试最理想的时间段应选择在夜深人静的深夜，但是在实际测量中却做不到这一点。而白天地脉动测试时，拾震器往往接收到一些振幅较大的干扰信号，干扰信号直接影响数据处理的结果，对结果造成不利的影响，因此需要对这些干扰信号进行筛选剔除。文中基于震相自动识别思想，对地脉动数据施加反STA/LTA算法，快速有效地对地脉动数据进行筛选，自动识别出异常信号，获得有价值的信号，为后期数据处理奠定了基础。

3.1 反STA/LTA算法识别流程

地脉动数据施加反STA/LTA算法的识别过程为：信号输入→计算→判断→拾取。输入参数有地脉动测试信号Y（t）、数据采样频率fs、短时窗的时长STA、长时窗的时长LTA、对应的判断阈值Rd和切片时间窗时长Tw。识别算法具体流程如下：

（1）计算地脉动数据t时刻的C F（t）值；

（2）计算STA/LTA法的t时刻的R（t）值；

（3）判别t时刻的R（t）值是否小于Rd值，若判断结果为是，则进行下一步判断；若判断结果为否，则进入下一时刻并返回（2）；

（4）将地脉动时程进行切片为t时刻到t+Tw时刻，分别计算出切片段的R（t）值，判别切片段的最大值是否小于Rd值，若判断结果为是，则进行下一步；若判断结果为否，则进入下一时刻并返回（2）；

（5）记录下T0=t时刻和T1=t+Tw时刻值，生成时间窗矩阵，并进行下一步判别；

（6）判别t时刻是否到达结束时刻t0，若判断结果为是，则绘制时间窗图形；若判断结果为否，则将进入t+Tw时刻，并返回（2）。

（7）如图3所示，为算法程序的流程图。

图3 STA/LTA算法流程图Fig.3 Flow chart of STA/LTA algorithm

3.2 反STA/LTA算法的影响因素探究

由于STA/LTA算法常常被用于地震震相的识别，震相识别与地脉动干扰信号识别的相似之处在于：两者都是在平稳信号中，识别突变的异常信号；而地震信号与地脉动信号的最大的区别在于：当地震信号来临时，振幅迅速增大，短时窗平均值会迅速变化，STA/LTA值短时间内增大，超过阈值，地震记录系统开始工作，记录地震信号；而地脉动最大的特点是振幅小、频率低，并且干扰信号的幅值，相较于地脉动信号的幅值相差不大。文中将从STA/LTA计算方法、特征函数和长短时窗的时长的选取等方面，讨论STA/LTA法对于地脉动测试数据处理的适用性。文中选取了一条在2018年5月27日松原发生5.7级地震后采集的典型地脉动测试记录作为研究对象。此地脉动记录点为姜家围子村委会前测试点，截取的时程为5 min，如图4所示。

图4 测点M1地脉动三分量记录结果Fig.4 Microtremors 3-components time-histories at site M1

3.2.1 不同的计算方法对算法的影响

如图5（a）和图5（b）所示，STA/LTA的算法分别选择的是递归法和标准法，递归法选择的是式（6）、式（7），特征函数选择的是式（5），短时窗的时长为4 s，长时窗的时长为20.48 s，阈值设定为3，时间窗长为20.48 s。在数据处理的过程中发现，地脉动的三分量中的垂直分量受脉冲影响较大，在施加STA/LTA算法的计算过程中，以垂直分量作为计算的首选计算样本，截取出数据窗之后再检验水平分量的R值是否超过阈值，若超过阈值则将该数据窗舍去，反之没有超过则进行下一步的计算。为了减小之后对地脉动数据分析处理的误差，相邻的时间窗设定重合50%，对比2种方法，结果表明：采用递归法，对干扰信号的识别的灵敏度较高，而传统的标准法，对干扰信号计算R值有一定的平滑作用，灵敏度较低。

如图5（a）、（c）和（d）所示，探究特征函数的选取对STA/LTA算法的影响。特征函数分别选择的是式（5）、式（2）和式（3），STA/LTA计算方法选择递归法，其他参数都设置一样的，短时窗的时长为4 s，长时窗的时长为20.48 s，阈值设定为3，时间窗长为20.48 s，相邻的时间窗重合50%；对比3种不同特征函数，如图5中的相同段e与f段，结果表明：特征函数选择平方法和高淑芳提出的改进法，对于干扰信号识别的灵敏度较高，而采用绝对值法的对干扰信号的识别较差。其中原因可能是地脉动的幅值较小，利用绝对值法计算出的特征函数，对R值的计算结果的放大作用不明显。

图5 不同特征函数选取对时间窗的影响Fig.5 The influence of different characteristic function selection on time window

3.2.2 不同参数选取对算法的影响

如图5（e）～（h）所示，探究短时窗和长时窗的时长选取值对STA/LTA算法结果的影响。其中计算方法选择递归法，特征函数选择的是式（5），短时窗的时长分别选取了2、4、8 s，长时窗的时长为10.24、20.48、20.72 s，阈值仍然设定为3，时间窗时长为20.48 s，相邻的时间窗重合50%。经过对比，如图6中的g与g′段，结果表明：随着短时窗的时长的增加，STA/LTA算法对干扰识别的灵敏度是降低，而长时窗的时长变化对识别信号的灵敏度影响较小，尤其是持续信号干扰段，算法的识别能力较差，可以通过减小短时窗时长或适当增大长时窗时长的方法提高精度。

图6 不同参数选取对的时间窗选择Fig.6 The influence of different parameter selection on time window selection

文中没有讨论反STA/LTA算法的阈值（THR）设定的问题，本篇文章中提及的所有对比试验阈值都取值3。其主要原因是：在地脉动分析数据时施加反STA/LTA算法的主要目的是对时程曲线进行快速切片提取、过滤干扰信号，以便于下一步的数据分析处理。在实际分析地脉动数据的工作时，由于地脉动是幅值比较小的是平稳信号，而干扰信号多为脉冲信号，易于判断出干扰信号，从而可根据时程曲线图来初步判断出阈值，并调整阈值的取值。

利用本文的结论，对地脉动采集数据施加反STA/LTA算法快速提取出有效信号段，进而利用单点谱比法求取场地的卓越频率。单点谱比法也称为H/V谱比法，是日本学者Nakamura［18］提出的一种场地周期计算的方法，在文献［17］中对单点谱比法的计算也有探讨，文中不再赘述。将提取出来的时间窗分别计算H/V谱，并计算其平均值和标准差，最终求得卓越频率。由于篇幅有限文中只给出了测点M1的结果，如图7所示。由于干扰信号都为脉冲信号，低频成分丰富，相比于施加反STA/LTA算法之后的结果，卓越频率略有减小。

图7 H/V谱比法结果对比Fig.7 Comparison of results of H/V spectra

对松原地区采集的108条地脉动数据记录施加反STA/LTA算法，参数设置为：STA分别为1 s和2 s，LTA为20.48 s，R分别设置为1.5和2.5。如图8所示，第1组（STA=1 s，LTA=20.48 s，R=1.5）的灵敏度较高，提取有效时窗个数平均值在10个；第2组（STA=1 s，LTA=20.48 s，R=2.5），提取的有效时窗个数平均值在15个；第3组（STA=2 s，LTA=20.48 s，R=2.5），提取的有效时窗个数平均值在20个。经过对比分析，结果表明：（1）当短时窗时长较小，施加反STA/LTA算法敏灵度较高，提取时窗个数较少；（2）当提取有效时窗个数小于10个时，计算得到的场地卓越频率有减小的趋势；（3）也存在当提取有效时窗个数达到15个时，计算得到的场地卓越频率呈现稳定的趋势。

图8 H/V谱比法结果对比Fig.8 Comparison of results of H/V spectra

4 结论与讨论

文中以地脉动测试数据为研究对象，施加反STA/LTA算法，提取地脉动信号的平稳段。通过改变计算方法、不同的特征函数和长短时窗的时长的选取，探究其对于地脉动测试数据分析的可行性。对108条地脉动数据记录分析，得到了如下结论：（1）反STA/LTA算法程序对地脉动测试数据的干扰信号，能够成功的识别；（2）使用递归法STA/LTA，特征函数选用平方法和高淑芳提出的改进法灵敏度较高；（3）短时窗的时长对结果影响较大，而长时窗的时长取值对结果影响较小；（4）利用算法提取的时窗数量越多，计算求得的卓越频率越趋于稳定。因此，推荐使用递归STA/LTA法参数设置为：短时窗（STA）的时长设定为1～2 s，长时窗（LTA）的时长设定为20.48 s，阈值（THR）设定为2.5，对地脉动测试数据进行处理。

在地脉动数据的处理时，需要人工截取数据段进行分析，存在偶然性。文中主要的探讨问题是：如何规范化处理地脉动数据。进而文中论证了将STA/LTA算法应用于地脉动数据处理的可行性，还要在今后的实际工程应用中验证其适用性。对于干扰信号的有效识别取决于算法的灵敏度和阈值的设定值。特征函数的选取也间接的影响算法的灵敏度，所以对特征函数的优化选取将是下一步工作的重点；同时在数据处理过程中，当采集的地脉动数据如图9所示时，测试过程中一直有脉冲荷载干扰时，算法无法有效识别出有效信号，只能通过增加短时窗时长或增大阈值。

图9 测点M13计算值Fig.9 Calculated value of M13