隧道突涌水防突结构微震响应初至拾取方法研究

李利平，陈彦好，靳昊，成帅，孙子正，胡慧江，陈迪杨

（山东大学 齐鲁交通学院，山东 济南 250002）

0 引言

隧道突涌水灾害是严重威胁我国交通基础设施建设安全的地质灾害之一，具有强破坏性、突发性等特点[1]，利用先进技术实现对地质灾害的时空监测预警是当前研究的热门方向。微震监测技术因其远程无损监测的优势得以在隧道工程、油田服务、矿山安全监测等领域广泛应用[2]。在隧道工程领域，我国岩溶发育的隧址地区，特别是高山险峻的西部，突涌水灾害由于高水压、高地应力因素使破坏力激增，灾变机理更加复杂。利用微震监测技术对防突结构的微震活动进行监测和分析，实现防突结构整体稳定性评价，是突涌水灾害监测的有效方法。因此，有必要开展一系列微震监测技术在突涌水灾害监测方面的应用研究，建立防突结构稳定性评价与微震结果之间的联系，为灾害监测预警奠定基础。

在微震监测技术应用过程中，震源定位和震源机制分析的关键是P 波初至时间拾取。突涌水过程存在水力压裂、通道流体撞击等微小振动信号，需要从中提取清晰的表征微震事件并实现精准定位。理论方法在理想信号或高信噪比信号中效果较好，但在低信噪比信号中具有一定难度，误差较大；算法选择、拾取阈值等人为因素也会导致误差加大。

对此，各国学者开展了大量研究。针对瞬态非平台的复杂微震信号，常规识别检测方法有时域法、频域法、时频域法及综合法[3]。基于频域、时频域的算法也被相继提出[4]，如时间-频率域能量比[5]、时频分析[6]等。此外，运用分形理论[7-8]、互相关[9]、神经网络[10]、数字图像[11]等技术亦可进行地震波初至时间拾取。通过对波形数据的统计分析，获取信号的偏斜度和峰度，进而拾取初至时间的PAI-S/K 法，可获得较准确的初至时间。在时域法——数值P 波初至拾取中，长短时均值比拾取技术是应用最广泛的技术之一[12]，该技术利用信号特征函数STA 和LTA 之比识别P 波初至相位，当STA/LTA 超过预定或动态指定的阈值时触发信号识别。随着现代计算机技术的迅猛发展，人工神经网络[13]、小波变换[14-15]、波偏振分析[16]等也逐渐应用到P 波初至时间拾取中，一定程度上给P波初至检测带来新的突破口，但这些方法需要1 个或多个预设标准，如检测间隔、阈值设置，不同的研究人员对此难以形成共识标准，因此分析结果存在一定差异。

在单自由度振动体系结构抗震研究领域中，已有研究以结构的总地震输入能及滞回耗能占总耗能的比例评估结构滞回耗能水平，相关理论研究在该领域已趋于成熟[17]，并且用于抗震设计的弹塑性单自由度体系的总能量谱研究也有较大进展[18]，许多学者采用以“单”代“多”的方法，研究单自由度体系和多自由度体系输入能的相关性，通过单自由度体系等效叠加得到多自由度体系的相关结论[19]。以上研究成果在微震信号分析方面尚无成熟应用，采用单自由度体系这一相对稳定的结构处理微震信号，利用体系能量谱的形式估计振动输入能，尝试实现对微震信号P 波初至特征的捕捉，提供了一种能量分析拾取算法的新角度。

1 单自由度体系振动规律分析

单自由度体系分为无阻尼和有阻尼2 种，前者振动总是以动能和势能交换为特征，并未考虑体系能量的耗散，即结构体系在振动过程中总能量保持不变，因而与能量大小密切相关的振幅始终不变[20]。实际上，单个微震事件信号振幅逐渐衰减，最终消失于背景噪声中，在单自由度体系中表现为质量块m 逐渐静止在静力平衡位置，即质量块m 从微震波动能量输入后进行的是有阻尼振动。

在单自由度体系中，阻尼主要由2 部分因素构成：一是外部介质的摩擦阻尼力；二是结构内部变形的内耗。阻尼力的性质比较复杂，为了进行阻尼结构的振动计算，采用线性粘滞阻尼理论求解。最简单的单自由度阻尼结构模型见图1。

图1 单自由度阻尼结构模型

假定阻尼力与变形速度成正比，但是方向与速度方向相反，公式为：

式中：D(t)为阻尼力；y为质量块的运行速度；c为结构系统的阻尼系数。结构的阻尼力由多种因素引起，机理复杂，多数情况将阻尼系数理解为实际结构中各种因素的综合阻尼系数，不影响最终计算结论。

设ξ为阻尼比，定义为：

式中：ω为结构自振角频率。

将式（3）代入式（2），设方程解为y=ert，则：

式（4）方程特征根为：

式（6）表达了单自由度含阻尼体系质量块在外界振动能量输入条件下，位移y随时间t的变化规律。当阻尼比设置为0.05，自振频率为1 时，含阻尼自由振动位移-时间曲线见图2。

图2 典型含阻尼自由振动位移-时间曲线

为了理解微震信号能量的输入过程，并理解体系内能量转化，对以下关键参量进行说明：

（1）振幅衰减。

由式（6）可知，含阻尼自由振动的振幅为Ae-ξωt，且振幅值按照指数e-ξωt的规律迅速衰减。如图2 所示，对于每1 个周期振幅的衰减情况，采用衰减率η表示振幅的衰减：

由此可见，阻尼造成振幅不断衰减，从初始时刻开始外界振动输入能逐渐耗散。

（2）频率减小。

为进一步表明结构体固有阻尼参数和自振周期对输入能衰减过程不同的反应，基于Matlab 完成相同初始位移、初始速度及时间步长条件下不同阻尼比的单自由度体系位移-时间曲线计算（见图3）。

图3 不同阻尼比的单自由度体系位移-时间曲线

如图3 所示，在输入体系能量一致的情况下，含有不同阻尼比的体系对输入能的振动反应有所不同。阻尼比越大，原始振动振幅和频率衰减越明显，最终趋于静止，阻尼比对输入能的耗散作用明显，同时也反映了单自由度模型对复杂振动过程的表征能力。本能量分析拾取法即利用单自由度体系的固有阻尼对微震信号输入能的耗散作用来表征有效微震信号的能量变化过程，从而拾取P 波初至时间。

2 单自由度体系振动能量反应

研究单自由度体系的微震能量输入，通过结构弹塑性地震反应的时程分析法，计算结构在微震作用下的振动输入能、阻尼耗能和滞回耗能时程曲线[21]，为微震事件信号中P 波初至时间拾取提供依据。单自由度体系在微震作用下的动力方程为：

将式（7）进行从y(0)到y(t0)的积分：

上式简化为：

式中：EK为单自由度体系的相对动能；ED为单自由度体系的阻尼耗能；EA为单自由度体系的变形能，由2 部分组成，一是可恢复的弹性应变能，另一部分是不可恢复的塑性累计滞回耗能；EI为微震事件输入能。事实上，如将包含微震事件信号的加速度振动谱作为输入能，其因阻尼逐渐发生能量耗散，直至结构停止振动，即体系EK=0。

3 微震信号能量辨识与初至拾取

计算截取2 种典型微震监测的现场数据作为输入能，基于Matlab 计算结构体系相对于地面的振动输入能、动能、阻尼耗能和变形能等各种能量时程曲线，解析微震信号时间序列上的能量变化过程。

将微震弹性波能量作为单自由度体系输入的总能量，受体系内阻尼力和变形作用，总能量以阻尼能和变形能的形式耗散，转化过程可由Matlab 计算实现。单自由度体系能量响应计算流程见图4。

图4 单自由度体系能量响应计算流程

采集某隧道微震监测现场数据，微震信号加速度时间序列见图5，输入能时程曲线见图6。

图5 微震信号加速度时间序列

图6 输入能时程曲线

由图5 可见，高信噪比微震信号P 波振幅突跳明显，波形轮廓清晰；低信噪比微震信号P 波突跳振幅不明显，P 波波至前存在噪声干扰。针对图5 的2 种现场典型微震信号进行分析，可得微震信号计算处理结果，输入能在体系能转换后的动能响应时程曲线见图7，阻尼耗能时程曲线见图8。

图7 动能响应时程曲线

图8 阻尼耗能时程曲线

在图6 所示各项能量中，输入能和阻尼耗能的表征曲线整体呈递增状。如图5（a）所示，动能曲线在200 ms 前递增明显，但维持时间较短，到达顶峰值后迅速衰减，在400 ms 后微震信号持续时间内在零线之上微小波动。从微震信号原始加速度波形分析，体系质量块受微震输入能影响，质点从静止开始振动，阻尼能量耗散作用开始，相应动能迅速增加，与微震信号转化为单自由度能量响应域的理论预期相符。如图6（a）所示，较高信噪比微震信号可清晰识别P 波初至振幅突跳位置，在输入能时程曲线相应时刻出现陡增。阻尼耗能随体系输入能一开始即发生能量耗散作用，所以阻尼耗能与输入能时程曲线变化趋势相对一致。动能是体系质量块运动状态的直接反映，体系动能曲线相对阻尼耗能变化复杂。在变化趋势上，动能时程曲线与微震信号加速度谱几乎同步达到峰值，然后迅速衰减。

低信噪比的微震信号处理是研究难点，微震事件的初至受噪声影响容易被掩盖。如图7 所示，微震信号来自同一次微震事件的加速度波形，但图7（b）相对于图7（a），有效信号在抵达之前存在一个噪声段，幅值相当于微震事件信号段峰值的1/2，容易造成算法识别错误。

如计算结果所见，体系输入能和阻尼耗能随时间累积，整体保持递增趋势。阻尼耗能激增点位置见图9。动能时程曲线随微震信号时间序列出现噪声信号段平稳、微震信号段内激增后迅速衰减的规律。如图7（b）所示，广泛应用的STA/LTA 算法拾取结果受噪声干扰，初至识别位置有所提前。这是因为算法触发机制设定在长短时均值比值超过预设阈值时发生拾取，而一般噪声在加速度时间序列内有着和微震事件相似的地震属性，存在误判的可能。当微震加速度谱转入单自由度体系作为能量变化考虑时，由于噪声信号具有随机特性，频散严重，计算的能量累积处于稳定水平，如图7（b）噪声段的动能时程曲线部分，但微震有效信号具有特定主频段和震源极性，所以噪声和微震事件信号在体系各部分能量时程曲线（见图8（b）、图9）变化上有清晰的辨识特征。

图9 阻尼耗能激增点位置

目前，常用的STA/LTA、AIC 等算法对微震初至的拾取，要求设定检测间隔或触发门槛，一定程度上增加了研究人员的主观因素[22-23]。将微震引起的质点运动加速度记录转化为体系能量域，使得在时间序列中的波至特征在体系输入能、阻尼耗能及弹性动能序列上得到凸显，以能量出现特征反应的时刻作为震波初至时间，进行定位计算。

通过上述计算结果分析可知，阻尼耗能在P 波相位到达之前为0 或接近0，并在P 波相位之后迅速积累。在阻尼作用下，随时间推移产生光滑的能量耗散累积曲线。因为阻尼能函数在低信噪比信号中反应仍然敏感，所以阻尼能变化过程可作为追踪和检测P 波初至时间的指标。根据阻尼耗能变化曲线，取耗能激增点位置0.16 s 为微震初至时刻，比采用长短时均值比检测延后了0.03 s（见图10（a））。针对低信噪比的微震信号，微震信号波至在阻尼耗能曲线上会出现第2 次增幅更大的跳跃，由此将微震初至定为1.85 s，与长短时均值比检测结果相差较大（见图10（b））。这是因为长短时均值比的检测方式一旦检测阈值确定，对于低信噪比的微震信号来说，其容易受噪声干扰，这就导致误将噪声识别为有效波至。

4 拾取误差分析

为了验证微震初至拾取方法的精度，通过建立震源定位双曲线控制的定位误差计算模型[24-25]，基于核密度估计函数得到能量分析拾取法与常用STA/LTA 在同一微震定位模型下的定位误差分布密度函数，从而判断2 种方法的拾取精度。以震波从震源到2 个传感器i和j为例，波至到时差曲线模型见图11。

图10 拾取结果对比

图11 波至到时差曲线模型

微震弹性波从震源到达传感器i和j之间的到时差（二维平面）表示为：

式中：vP表示弹性波P 波波速。式（10）是震源定位双曲线（平面范围）控制方程，形成的曲线为ti-tj到时差轨迹线，线上每点到2 个传感器的距离为vP（ti-tj）。在概率论方法中，概率密度函数是一个描述此变量的输出值，在某个确定取值点附近的可能性的函数。对模型台网内所有台站添加了满足相同的正态分布ξ～N（0，σ2I）的旅行时随机误差，I为单位矩阵，σ为随机误差的方差。利用能量分析拾取法与STA/LTA 计算同一观测系统下的定位双曲线，检验2 种算法在含随机误差情况下的拾取结果，计算定位误差概率分布。通过获得的误差概率密度分布函数判断不同拾取结果下的定位精度，从而实现对2 种方法P 波初至拾取精度的比较。经Matlab 计算，进行2 种方法的反演定位结果与精确值对比，可得P 波初至拾取结果误差分析（见图12）及各自定位误差概率密度分布（见图13）。

图12 拾取结果误差分析

图13 水平向定位误差概率密度分布曲线

如图13 所示，能量分析拾取法的结果在水平向上误差控制更好，误差分布带宽更窄，定位误差在7.8 m范围内波动，STA/LTA 的定位误差则在13.3 m 范围内波动。较小的误差分布带宽，表示在零误差范围内具有较高集中度，说明基于能量分析的拾取方法在反演计算微震事件位置时，对于定位误差具有更好的控制力。

5 结束语

基于单自由度振动体系的能量分析过程实现有效波动初至的拾取，并采用双曲线定位模型进行验证，研究结果表明：

（1）不同信噪比的微震信号在单自由度振动体系能量分析中可以捕获明显特征点。高信噪比信号段在阻尼耗能、动能曲线上激增点凸显，输入能曲线激增幅度大。低信噪比信号在阻尼耗能曲线上出现阶段性递增，为微震事件波至拾取提供参考依据。

（2）由于能量特征点的检测不需要预设阈值，避免了其他预设阈值方法中的人为误差，具有更好的误差控制，通过数值计算验证，能量分析拾取法可获得微震反演更窄的定位误差分布。

隧道工程监测领域所研究的微震监测对事件定位精度比天然地震监测的精度要求更高，众多研究旨在解决这项技术在尺度变化下的适用性问题。本研究的微震事件初至拾取作为复杂系统的一环，力求提供一种新方法和思路。要做好定位精度的控制仍需考虑前后环节的衔接处理和相互影响，如数字滤波技术和定位计算，未来很多相关研究工作必定继续开展。