1. 引言
随着全球气候变化的影响,极端天气时有发生,滑坡灾害在山区尤其是湿润气候地区发生的频率显著增加,造成的人员伤亡和财产损失不断攀升[1] [2]。在中国,重庆、四川、云南等地区由于其地质构造复杂、降水量大,成为滑坡灾害的高发区域。地质灾害成因复杂,通常受区域地质构造、地貌演化、水文气候条件及人类工程扰动等多重因素的共同影响,历来是地质学及灾害防治领域关注的重点问题[3]。与此同时,人类的工程活动,如矿山深部开采[4]、地下储气库运营、隧道挖掘爆破等,剧烈扰动岩体原始应力场,极易诱发围岩失稳、冲击地压乃至岩爆等灾害,严重威胁施工安全与工程稳定性[5]。微震信号是岩土体在变形、破裂过程中释放的低频地震波,具有传播距离远、衰减小、穿透能力强的优点。与传统的地表位移监测、应力测量等方法相比,微震监测技术通过布设于地下或者结构内部的传感器阵列,时间分辨率更大、空间覆盖范围更广,可以实时捕捉事件演化前兆,从而实现灾害的早期识别预警,其为地质灾害提供了重要的预警和分析手段[6]。
微震监测的效果依赖对海量低信噪比数据的精准处理,处理流程包括信号预处理、事件检测、震相拾取、事件分类与震源定位等多个环节,其中事件识别(长短时窗法)、到时差定位等方法构成了技术核心,在线监测能力则是实现实时预警的保障[7]。该领域目前面临着软件工具困境。一方面,虽然主流商业软件(如RadExPro [8]、Geogiga Seismic Pro [9]等)功能集成化高,界面友好,解决方案成熟,但是其核心算法被封装为黑箱,用户无法了解数据处理过程。这导致研究者难以验证、复现或者改进其处理结果,严重制约了算法的透明化研究与技术的自主创新。另一方面,学术界常见MATLAB或Python脚本[10]的算法代码虽然可以保证算法性能与可修改性,但普遍存在功能单一、缺乏集成、无图形用户界面等问题。研究人员需要投入大量精力进行二次开发和调试,才能将零散的脚本串联成完整的工作流,用户体验较差,极大地阻碍了科研效率的提升与新方法的推广应用。
本文工作提供了交互式微震信号分析环境,集成机器学习分类器用于事件识别与滤波,实现了一套高精度的到时差定位(Time Difference of Arrival Localization, TDOA)流水线,将以上功能无缝集成,既服务于前沿算法的科学研究,也面向工程现场的实时监测应用。
2. 系统与算法
本研究构建了一套集成自动化监测与交互式分析功能的微震分析系统。系统以微震数据池为核心,底层集成多数据源解析、巴斯沃特滤波、连续小波变换(Continuous Wavelet Transform, CWT)时频分析、长短时窗(Short Term Average/Long Term Average, STA/LTA)事件检测、机器学习分类、P波S波分析、到时差定位及落石反演等核心算法模块,实现从数据处理到参数解译的全链条能力。上层通过自动监测界面与手动分析模式双入口设计,分别提供实时处理流水线与包含交互预处理、信号解译和可视化探索的精细分析环境,有效兼顾监测效率与科研深度。如图1所示。
Figure 1. System architecture diagram
图1. 系统架构图
2.1. 微震信号预处理
趋势分量是数据整体变化的趋势,不包含局部波动或者噪声。在微震信号的处理过程中,消除趋势分量有助于分析信号特征。常见的消除趋势分量的方法有滤波法、小波变换法、最小二乘法和经验模态分解法[11]。滤波法设计简单,可以高效地获取预想的频率区间。如图2所示是一段爆破信号采用滤波法消除趋势分量的对比。
软件提供了低通、高通、带通和带阻四种滤波模式。用户可交互式地设置滤波器的关键参数(截止频率、阶数),并实时观察滤波效果。巴特沃斯滤波器以其平滑的频响特性著称,但其过渡带相对平缓。因此,对于要求极为陡峭截止的场景,用户可通过适当提高滤波器阶数来改善阻带衰减,但这会相应增加计算量和相位延迟[12]。
Figure 2. B1 blasting signal eliminates component trends in seismic stations: (a) Amplitude and time curve; (b) Energy and frequency curve
图2. B1爆破信号在地震台站消除分量趋势:(a) 振幅与时间曲线;(b) 能量与频率曲线
在微震信号处理中,时频分析是揭示非平稳信号动态特性的核心手段。鉴于微震信号的频率成分随时间动态变化,传统的傅里叶变换仅能提供全局频谱的局限性,本软件集成了连续小波变换方法以生成高分辨率的时频谱。该功能通过将信号与一系列经过缩放和平移的复高斯母小波(“cgau8”)进行卷积,精确地绘制出信号在时间–频率二维平面上的分布图,从而直观地揭示微震事件发生的精确时刻及其优势频率范围。如图3所示,即为利用本软件对某次微震事件信号进行时频分析所得的可视化结果。
Figure 3. Time and frequency analysis chart
图3. 时频分析图
图中上部为时频分析图,该图采用连续小波变换生成,横轴为时间序列,左侧纵轴为频率,右侧纵轴为振幅。颜色深浅表示信号在特定时刻和频率上的能量强度,暖色调代表高能量区;下部为傅里叶频谱图,基于傅里叶变换生成,横轴为频率,纵轴为傅里叶系数的振幅,反映了信号在整个分析时间段内的全局平均频率成分。上方的时频图精确定位了微震事件的发生时刻并揭示了其频率成分随时间的变化过程,完美展示了信号的非平稳特性。下方的频谱图则从全局角度验证了事件的优势频率,并提供了信号的总体频谱结构。两者结合,为研究者全面解读微震信号的动态特性和频率属性提供了强有力的证据。软件将此两种视角集成于同一界面,极大地便利了对比分析与深入解读。
2.2. 事件检测
为实现微震事件的自动识别与提取,采用了一种基于长短时平均比值的经典触发算法。该算法通过计算信号瞬时能量与背景噪声能量的比值,有效检测信号中的突发性事件,其计算效率高,非常适合对连续波形数据进行批处理与实时监测。本工作对该算法进行了工程实现与优化,并将其封装于多线程任务中,确保了在处理长时序数据时系统界面的响应能力。
对于一个离散的时间序列信号x(n),其STA和LTA的计算首先基于信号的平方值(即能量) E(i),短时窗和长时窗分别通过在设定的时间窗口内对E(n)进行滑动平均计算得到:
(1)
(2)
其中,
和
分别为短时间窗和长时间窗的样本点数。特征函数R(n)即为两者的比值:
(3)
当信号平稳时,R(n)在值1附近波动;当有事件发生时,STA迅速上升而LTA变化缓慢,导致R(n)显著增大。当R(n)超过预设的触发阈值T1时,算法初步判定事件开始。
为提高检测可靠性,本实现引入了一套双阈值判据与状态管理机制。其核心流程如下:
1) 初步触发:当R(n) > T1时,事件长度计数器开始累积。
2) 事件确认:当事件长度达到预设的最小事件持续时间Emin时,算法回溯检查,要求事件起始点附近的R(n)值也必须超过一个较低的确认阈值T2 (T2 < T1)。此步骤有效避免了因个别噪声尖峰造成的误触发。
3) LTA冻结:事件被确认后,LTA值将被冻结在事件开始前的水平。此举旨在防止事件本身的高能量抬升LTA的计算基准,从而保证事件持续期间触发判据的稳定性,避免算法过早终止。
4) 事件终止:当R(n)值回落至T1以下,且持续时长超过最小事件间隔Imin时,判定事件结束。
如图4所示为某事件检测示例结果。所有算法关键参数均可通过软件界面开放用户配置,以适应不同信噪比环境与信号特征的微震数据。检测到的事件及其时间戳将被记录并输出为结构化的CSV文件,同时通过信号槽机制实时通知主界面更新,从而形成一个完整的自动化事件检测、参数化与可视化流程。
2.3. 事件分类
在微震事件的监测与分析中,信号的识别、分类与定位是至关重要的步骤。通过对微震信号的识别与分类,不同类型的事件得以有效区分,从而为后续的定位及滑动面拟合提供清晰的数据。在2021年5月1日至2022年6月30日期间,对所有监测网络内微震信号识别分类后,筛选出至少由4个基站监测
Figure 4. Event detection example
图4. 事件检测示例
到的事件共195个,删去明显不符合实际情况以及噪声的数据,最后得到138个有效的微震事件定位结果。对大量识别出的微震事件进行时频分析,并结合Langet等[13]斜坡微震信号分类中的经验和研究成果,对138个事件进行详细分类,并根据实际观测数据特点对各类信号进一步细化,微震信号被分为以下6类,每一类的定义和特征以信号的频率、持续时间、波形等物理特性为依据,并结合地质背景进行了合理的推测。
高频率事件(High Frequency Events, HF):高频率微震事件的频率范围为30~60 Hz,持续时间较短(通常小于5 s),并表现为突发性波动。低频率事件(Low Frequency Events, LF):与HF事件相比,低频率事件的频率范围为4~20 Hz,持续时间小于5 s。LF事件可能是岩土体的缓慢破裂或持续滑动过程的体现。系列高频事件(High Frequency Series Events, HFS):HFS微震事件包括频率较高的信号序列,并且通常表现为短时间内的频繁发生,持续时间通常大于5 s。系列低频事件(Low Frequency Series Events, LFS):系列低频事件的特点是低频微震信号的多次连续发生。双频率事件(High-Low Frequency Events, HLF):HLF微震事件具有明显的双频率中心,分别位于10~25 Hz和20~40 Hz。这类事件可能对应了深层滑坡岩体内部更为复杂的破裂模式或滑动机制,可能源自岩土体的深层裂缝或破裂面。双频率系列事件(High-Low Frequency Series Events, HLFS):HLFS微震事件是双频率事件的扩展,在双频率中心事件类型的基础上,持续时间较长。
通过对微震信号进行分类,为后续的事件定位与滑动面拟合提供了数据基础。分类考虑了信号的物理特性,结合了不同事件的地质背景和实际情况,对信号类别进行了合理划分。各类别微震信号的波形和时频图如图5所示[14]。
Figure 5. Event waveforms and time-frequency graphs of each category
图5. 各类别事件波形和时频图
2.4. 震源定位
基于“到时差”的网格搜索方法是一种常用的微震震源定位技术。核心是通过计算理论到时与实际观测的到时之间的差异,在三维地形网格中寻找最优震源位置。首先假设空间为均匀的,将监测区域的三维地形图网格化,获取所有网格节点的坐标,搜索最接近震源位置的网格节点
。已知基站坐标为
(
),基站接收微震信号的到达时间
(
)。给定信号传播速度v,则震源到基站的传播时间为:
(4)
若微震事件的发生时间为t0,信号在基站的到达时间为:
(5)
其中,
为地质环境以及信号传播计算公式造成的信号传播时间的误差。为消除震源发生事件
的未知性,计算网络中所有基站到达时间的平均值:
(6)
其中,
为所有基站信号传播时间的平均值。将信号在基站的到达时间和平均值作差值,可以得到:
(7)
通过该式将定位问题转化为仅依赖空间坐标的方程。求解由多个方程联立而成的方程组以实现微震信号的定位。
残差
是实际到达时间与理论到达时间之差,残差平方和
被用作目标函数来衡量震源位置的合理性[7]。
(8)
通过最小化残差平方和,找到一个震源位置,使得理论传播时间与实际传播时间之间差异最小,从而确定整个空间震源位置最优解。
受限于实际地形与台站垂直覆盖的非均匀性,基于时差法的网格搜索在垂直定位精度上往往逊于平面方向,且常呈现向监测区域底部集中的系统性偏差。研究分析认为,这一方面源于地形尺度差异致使误差函数被平面分量主导,算法在更新时倾向于调整水平坐标而“钝化”了垂向优化;另一方面,基站高程的集中分布削弱了对震源深度的敏感性,从而形成垂向模糊带。
针对传统方法单纯追求残差极小而忽视物理可行性(如定位结果偏离实际地形或落入模型边缘极小点)的弊端,本研究引入高程惩罚机制[15]。通过构建三维掩膜并实施形态学腐蚀策略[16],有效剔除低质量边缘信号区,防止震源被错误定位至地表,确保传播路径的物理真实性。此外,设定100米高程下限并在评分函数中增设高程差异惩罚项,既规避了深部误判,又通过聚焦有效计算区域显著提升了反演精度与效率。定义震源搜索空间中任意网格点i的高程为
,台网布置的中心高程面为
。为防止定位结果非物理地偏离监测层位,引入高程偏离惩罚函数
:
(9)
(10)
其中,
表示N个监测台站的平均高程,该平面反映了监测台网的主体控制层位;
为正则化权重系数,用于平衡数据拟合残差与深度先验约束的权重。
在构建设立网格点的概率密度评分时,将上述惩罚项以指数衰减的形式融入评分体系。传统的基于残差的似然函数
被修正为后验概率形式
:
(11)
式中:
为网格点i处的走时残差统计量;
为控制评分分布形态的尺度因子;第二项
即为高程约束项。
由式(11)可知,该算法本质上是将求解目标从单纯寻找“残差极小值”转化为寻找“残差较小且高程合理的概率最大值”。当两个网格点的走时残差
相近时,距离台网中心高程面更近的点(即更小)将获得更高的评分权重。这种机制在不改变水平定位结果的前提下,通过压制垂直方向的离群解,有效收敛了深度方向的置信区间,减小了因基站垂直分布集中而产生的系统性深度偏差。
在微震信号的定位过程中,由于滑坡区域范围较大,为了提高定位精度与计算效率,可采用多尺度网格搜索的方法。其核心是结合不同尺度网格的优点,逐步细化定位区域,优化定位过程。粗网格设计是整个多尺度网格策略中的第一步,为了在保证数据有效的前提下快速缩小震源的搜索范围,对进一步定位提供引导,需要选择一个合理的网格间隔。网格间隔决定了网格的分辨率以及搜索运算量,如果间隔选择过大,网格太粗,可能由于震源恰好在两个网格点之间或分辨率过低,误差计算偏离实际震源位置致使漏掉真实震源点,也存在无法后续需要盲目搜索较大区域;如果间隔过小,则会增大计算时间与计算量,失去多尺度搜索的优势。在完成粗网格定位并初步确定小区域后,对定位区域进行点扩展,其中东西向、南北向扩展形成2级网格,在2级网格内建立更小尺度、更高分辨率的网格,再次进行网格搜索。此步骤将显著缩小搜索范围,并提升定位精度。重复上述过程,可根据需要继续进行多级搜索。每一级都在上一级确定的最优子区域内,使用更精细的网格进行迭代计算,直至网格尺度满足预设的定位精度要求。如图6所示为多尺度网格搜索流程图。
Figure 6. Multi-grid search flow chart
图6. 多网格搜索流程图
3. 实验分析
3.1. 交互式分析功能展示
为实现对自动化处理结果的深度验证与优化,本研究开发了一套综合性的交互式地震波形分析平台。该平台不仅提供了对原始数据和算法中间产出的全面可视化,更赋予研究者从数据预处理到结果生成全流程的人工干预能力。
研究者可手动选择或排除特定台站的通道(如Z、N、E分量)进行处理,这对于处理部分台站数据缺失或质量不佳的情况至关重要。同时,支持对连续波形进行任意时间窗的裁剪,以聚焦于目标事件或剔除已知的强干扰时段,从源头上提升输入数据质量。平台提供带通、高通、低通等可调参数的实时滤波功能。用户可动态调整滤波频带,并实时观察波形变化,从而为不同地质环境或频段特征的事件(如高频的爆破、低频的远震)确定最优预处理方案。结合时频分析(如短时傅里叶变换频谱图)功能,用户可直观判断信号和噪声的频域分布,为滤波器设置提供科学依据。
研究人员能够超越预定义模型,针对特定区域和科学问题创建专属的分类器。该功能模块允许用户手动创建或导入自定义的分类标签文件,从而定义其研究所需的分类体系,极大增强了系统的灵活性和适用性。经过上述交互流程产生的、经过专家知识校验的优质数据集和优化后的模型,可保存为专属模板。
3.2. 自动化流水线
本研究开发了一套基于生产者–消费者模式的在线自动化处理流水线,实现了从数据采集到定位输出的全流程集成化处理。研究者只需预先导入数字高程模型(DEM)并配置监测工作目录,系统即可自动完成连续波形的在线事件检测、实时事件分类与高精度局部定位。
整个系统通过多实例并行处理架构实现计算资源优化,在保证实时处理性能的同时,确保了从事件触发、分类到地形校正定位的全流程自动化运行。该方案不仅为研究者提供了便捷的一键式处理体验,更重要的是通过地形融合算法实现了局部区域定位精度的本质提升,为复杂环境下的微震监测提供了可靠的技术支撑。
4. 工程验证
为验证本软件的实际应用效能,研究选取重庆市旧县坪滑坡的微震监测数据进行处理分析。旧县坪滑坡位于青龙街道复兴社区5组,建民村1、2、3组。平面形态呈近圆弧状,滑体左侧、右侧均以冲沟为界,后缘以基岩陡壁为界,滑坡体前缘由于三峡库区蓄水部分被长江水淹没,滑坡整体边界条件较为清楚。根据前期资料,滑坡体前缘高程约为95 m,后缘高程385 m (其后缘堆积体滑坡最大高程约515 m),滑坡高差290 m。如图7所示为滑坡场地与剖面图。
软件首先对连续波形数据进行自动事件检测,识别出有效的微震信号;随后,基于信号的频谱、波形等特征对检测到的事件进行精确分类;在此基础上,利用台站布局与到时信息,对分类后的事件进行高精度定位。在获得微震事件精确定位与分类结果的基础上,若判定为已知滑坡灾害事件,可进一步开展滑动面几何形态的反演研究。滑动面作为滑坡体与稳定基岩之间的剪切破坏界面,其空间展布特征对揭示滑坡动力学机制及防灾治理具有决定性意义。
针对微震事件在滑坡体内部分布不均的特征,本研究采用主成分分析(PCA)与趋势面–残差融合拟合相结合的方法,对潜在滑动面进行三维几何建模。具体而言,首先通过PCA确定微震事件点云的主优势方位,初步提取滑动面的宏观趋势特征;继而采用趋势面拟合方法构建滑动面的背景几何形态,并通过残差分析捕捉局部结构调整量;最终通过趋势结构与残差场的融合,实现滑动面三维形态的高精度重建。
Figure 7. Jiuxianping landslide site and profile
图7. 旧县坪滑坡场地与剖面图
该方法有效克服了因微震事件空间分布不均导致的形态刻画困难,能够在稀疏观测条件下恢复滑动面的连续几何特征,为滑坡灾害的稳定性评价与防治措施制定提供重要的结构依据。如图8所示为滑动面拟合过程。
最终叠加趋势面与融合残差场得到可以表达整体趋势与局部扰动特征的滑动面模型。如图9所示为滑动面三视图。
图9所展示的三维DEM可视化图像中,叠加绘制了最终构建的滑动面以及三视图。从图中可以清晰看到,拟合得到的滑动面较为准确地穿越了微震事件密集分布的区域,其几何形态与主要震源活动带之间呈现出明显的空间对应关系,进一步验证了提出方法在滑动面建模方面的适用性与可靠性。滑动面与微震事件在垂向分布上具有较强的耦合性,震源密集区域主要分布在滑动面上下附近,说明该面不仅在几何构型上较为合理,也在一定程度上反映了滑体内部破裂带的实际发育特征。尽管震源定位过程中存在一定误差,并且拟合过程中未显式引入岩体力学参数约束,但整体拟合结果依然为后续滑坡结构划分与不稳定边界识别提供了有力支撑与数据基础。
Figure 8. Sliding surface fitting process
图8. 滑动面拟合过程
Figure 9. Three views of the sliding surface
图9. 滑动面三视图
如图10所示为拟合滑动面三维图与俯视图,拟合滑动面高程范围为100~400米之间,定位点高程范围为124.6~233.6米之间,图中定位点在滑动面上呈离散分布,集中于中下部。深部位移监测基站CZK0201处滑坡体高程在170.0米到180.0米范围均有分布,CZK0202处滑坡体高程在88.8米周围。定位点高程集中区有2处,分别为200米与140米,200米处定位点集中靠近监测基站CZK0201且在滑坡体中的监测基站位于上部,因此该基站测得的滑动面高程可作为实际滑动面范围的参考;140米处定位集中点位于监测基站CZK0201与CZK0202之间,所以该基站测得的滑动面高程也可作为实际滑动面范围的参考。通过对比不同方法得到滑动面的高程范围,这不仅表明基于软件定位结果反演滑动面的方法是可行的,也有力地证明了本软件从事件检测、分类到精确定位的整个处理流程具有高度的可靠性与工程实用价值。软件的输出结果能够有效地服务于滑坡内部结构解译与稳定性评价,为地质灾害的精准研判提供了强有力的技术工具。
Figure 10. Fit the sliding surface
图10. 拟合滑动面
5. 结论
本研究成功研发了一套集成化的微震数据分析软件系统,通过对系统架构、算法性能和应用效果的全面评估,主要获得以下结论:
1) 软件系统成功集成了从原始信号预处理、事件检测分类到多尺度定位的完整处理流程,并创新性地融入了崩塌落石物理参数反演功能。通过生产者–消费者模式的在线处理架构与交互式分析模块的协同工作,系统既满足实时监测需求,又支持科研需求的精细分析,真正实现了从数据处理到物理解释的一体化解决方案。
2) 自动化流水线模式提供了标准化的批量处理能力,保证了工程应用的处理效率和一致性;而交互式分析模式则通过人机协同机制,允许专家介入验证和优化自动处理结果。这种设计既保持了算法处理的客观性,又融入了领域专家的主观知识,显著提升了结果的可靠性和实用性,为科研成果向实际应用的转化提供了有效途径。
本系统在复杂地质条件下仍存在局限,主要表现为对速度模型精度和台阵几何结构的依赖性较强,当前模型对非均质地质和各向异性特征的刻画能力有限。未来将致力于集成各向异性速度模型、开发波形反演震源机制解算功能,并通过云平台部署提升系统在复杂环境下的适应性与计算效能,推动微震监测向智能化方向发展。