基于小波分析对形变扰动信号的触发检测方法

刘明辉 李 江 沙成宁 周银兴 陈 阳 崔仁胜 朱小毅 林 湛

1 中国地震局地震预测研究所地震预测重点实验室,北京市复兴路63号,100036 2 青海省地震局,西宁市兴海路1号,810001

震源区应力场在地震前的短临发震阶段处于变动中,越接近发震,发震断层在断裂前扩展具有的脉冲性、间隙性表现越突出,将引起以突跳性、突发性和几起几落为主要特征的前兆现象[1],常被称为前驱波、慢地震、静地震等。这种前兆现象一般在震前7 d至地震发生前出现,周期为几十秒到几十分钟[2-3],可能会被钻孔应变仪、重力仪、倾斜仪、超长周期地震仪等仪器记录下来[4-7],但仪器同时也会记录到地球固体潮和干扰信号。地震前驱波和干扰信号相对于地球固体潮来说周期短、能量弱,叠加在固体潮记录上形成颤抖扰动,不易直观提取和判别,给地震研究带来较大困难。

钻孔应变仪的记录频带范围为DC～10 Hz,能够记录到清晰完整的固体潮[8]。通过计算钻孔应变观测的4个分量之间是否满足自洽方程,可检查观测数据的稳定性和可靠性[9]。小波分析是地壳形变资料分析中一种常用的方法[10],STA/LTA触发检测算法是一种经典的大信号检测方法,广泛应用于日常地震监测、地震烈度速报、地震预警和地震科学研究[11-12]。本文以2008年汶川8.0级地震前玉树地震台钻孔应变观测数据为例,用小波分析从钻孔应变观测数据中提取地震前的扰动信号,用STA/LTA触发检测算法进行扰动信号触发检测。

1 数据处理方法

1.1 自洽检验

钻孔应变仪的4个观测分量(S1、S2、S3和S4)按照顺时针方向相差45°依次分布,如图1所示。其中,S1与S3、S2与S4之间呈90°正交,分别合成分量S13、S24,自洽方程为[9]:

图1 钻孔应变观测分量分布示意图Fig.1 Observation components distribution sketch of strainmeter

S1+S3=S2+S4

(1)

式(1)成立时,4个分量的观测数据可靠,能真实地反映地层应变变化。本文通过计算S13与S24之间的相关系数进行自洽检验。

1.2 相干分析

利用相干分析得到S13和S24信号频率相干程度。2个信号x(t)和y(t)间相干函数为[13]:

(2)

式中,Pxy(ω)为x(t)和y(t)的互功率谱密度,Pxx(ω)、Pyy(ω)分别为x(t)和y(t)的频率域内的自功率谱密度。Cxy(ω)的取值范围为[0,1],若y(t)与x(t)完全不相关,则Cxy(ω)=0;若y(t)为x(t)的线性响应,则Cxy(ω)=1。

1.3 小波分析

小波分析将信号分解到尺度域,通过多分辨率的分解使原始信号中的弱信号成分变得突出,具有优良的时频局部能力。给定一个基本函数[14]:

(3)

式中,a、b为常数,a>0,φa,b(t)由基本函数φ(t)先移位再伸缩后得到,若a、b不断变化,可得到一组φa,b(t)。信号x(t)的小波分析为:

(4)

式中,*代表共轭,a、b和t为连续变量,式(4)的频率域为:

WTx(a,b)=

(5)

式中,X(Ω)为x(t)的傅里叶变换,如果φ(t)是幅频特性比较集中的带通函数,则小波分析就具有表征待分析信号在频域上局部性质的能力。当a较小时,时域观测范围较小,在频域上相当于用高频小波进行细致观察;当a较大时,时域观察范围较大,在频率域上相当于用低频小波进行概貌观察。这为选择小波分析的参数及提取扰动信号的目标层提供了参考。

1.4 时间-频率-能量分析

对小波分析提取出来的扰动信号进行时间-频率-能量分析,得到时间-频率变化关系及能量变化趋势。短时傅里叶变换以固定的滑动窗口对信号进行分析,随着窗函数的滑动,可表征信号的局域频率特性。令g(t)为时间宽度很窄的窗函数,并沿时间轴滑动,信号z(t)的短时傅里叶变换(STFT)定义为[15]:

(6)

由于信号z(t)乘以一个相当窄的窗函数g(u-t)等价于取出信号在分析时间点t附近的一个切片,所以STFT(t,f)可以理解为信号z(t)在分析时间t附近的局部频谱。

1.5 STA/LTA检测算法

STA/LTA检测算法中,STA是短时间窗口内信号幅度的绝对均值,目的是捕获突发大信号;LTA是长时间窗口内信号幅度的绝对均值,目的是得到突发大信号到达之前的平均背景噪声水平,作为突发大信号触发检测的参考基础。STA/LTA算法的步骤为:首先根据检测目标的需要选择适当长度的短时间和长时间窗口,分别计算对应时间窗口内信号的绝对均值STA和LTA;其次计算这2个均值的比值,即STA/LTA;最后将比值与启动触发阈值进行比较。若比值大于启动触发阈值,认为检测到大信号;当比值低于停止触发阈值时,认为大信号结束。

2 数据处理结果

2.1 数据及其自洽检查

2008年汶川8.0级地震前,玉树地震台钻孔应变仪记录到的同震效应信号幅度非常大。因此本文以汶川地震前5 d的记录数据为例,用小波分析提取其扰动信号,进行仪器触发检测。利用Tsoft软件对玉树地震台钻孔应变记录到的4个分量进行去零点漂移等预处理[16],波形如图2(a)所示,根据式(1)计算分量S13和S24,波形如图2(b)所示。

图2 玉树地震台钻孔应变仪记录Fig.2 The records of borehole strainmeter at Yushu seismic station

计算可得,S13与S24的相关系数为0.999,满足自洽方程,说明玉树地震台的钻孔应变仪工作正常。由于固体潮的周期长、能量强,而扰动信号的周期短、能量弱,从图2(a)和2(b)中无法直接分辨出明显的扰动信号。

2.2 相干分析

根据式(2)计算合成分量S13和S24的相干系数,结果如图3所示。图3显示,S13和S24中不相干的频率主要分布在120～130 s(2 min左右)和1 000～6 000 s(16～100 min),也是前驱波、慢地震、静地震等的频带范围[17-18]。S13和S24在固体潮频段的信号是完全相干的,说明玉树台钻孔应变仪工作稳定可靠。

2.3 小波分析和时间-频率-能量分析

小波分析将原始信号分成不同的层,使其含有不同频率的信号成分[19]。研究表明,地震前兆异常的判据不仅包括幅度变化,还包括频率变化等[20]。根据式(5)选择的时域观察范围较大,而在频率域上相当于用低频小波作概貌观察。利用db小波分析原始数据,将第1层信号提取出来,并进行时间-频率-能量分析,结果如图4所示。从图4(a)中可以看出,50～60 h之间存在一个扰动信号,其在图4(b)中形成了一个大峰值,在图 4(c)中显示的能量最强,4个分量中尖峰的信号频率各有不同,能量也各有不同,S2分量中的频率多于其他3个分量。从图4(a)中无法直接看到汶川8.0级地震前的扰动信号,但从图4(b)和4(c)中可以看到地震前的扰动信号。图4(c)显示,峰值之后信号能量逐渐增强。

为进一步分析,将汶川地震前24 h的信号作归一化处理后再进行时间-频率-能量分析,结果如图5所示。从图5(a)看出,汶川8.0级地震发生前2 h,4个分量提取的信号中都有明显的扰动成分,S1与S3的扰动信号形态大致相同,S2与S4的扰动信号形态大致相同,4个分量都具有间歇性的特征。图5(b)显示,不同频率的扰动信号的能量有差异。

图5 归一化原始信号、扰动信号及扰动信号时间-频率-能量分布Fig.5 Normalized raw signal, disturbance signal and distribution of time-frequency-energy of disturbance signal

2.4 触发检测

STA/LTA触发检测算法有4个基本触发参数:STA和LTA的时间窗口长度、启动和停止触发阈值。图4(c)和5(b)显示,汶川8.0级地震的扰动信号中能量较高的频率范围主要为2 000～6 000 s。将STA和LTA时间窗口长度分别设置为1 000 s和9 600 s,启动和停止触发阈值分别设置为1.6和0.8,对6个分量(S1、S2、S3、S4、S13和S24)中利用小波分析提取出来的扰动信号进行触发检测。结果表明,6个分量的检测效果都很好。

汶川地震前5 d及12 h的S1分量触发检测效果分别如图6(a)和6(b)所示,图中灰色为提取的扰动信号,红色为触发启动标识,绿色为触发停止标识。图6显示,STA/LTA检测算法检测出了S1分量中利用小波分析提取的扰动信号的开始时间和触发时间。如果用计算机进行实时在线触发检测,一旦检测到扰动信号就可自动启动声、光和图形报警,并立即通过电话、微信等方式通知值班人员进行分析处理[21]。

图6 S1分量触发检测效果Fig.6 Trigger detection results of S1 component

3 结 语

地形变观测过程中往往存在干扰信号和地震前兆信号叠加在固体潮上形成颤抖扰动的情况,给地震研究带来较大的困难。本文首先利用钻孔应变自洽方程检验玉树台钻孔应变仪的工作状态。结果表明,汶川8.0级地震前玉树台钻孔应变仪的观测数据满足自洽方程,仪器工作正常,长周期信号记录稳定可靠。然后利用互相干分析得到玉树台钻孔应变仪S13和S24两个合成分量中信号的相干频率,再通过小波分析分别从S1、S2、S3和S4分量中提取扰动信号,并进行时间-频率-能量分析,得到地震前扰动信号的时间-频率-能量分布。最后利用STA/LTA检测算法对提取出来的扰动信号进行触发检测。结果表明,选择合适的时间窗口长度和触发阈值能够有效检测出扰动信号的开始时间和结束时间。汶川8.0级地震前2 h,玉树台钻孔应变仪4个分量都检测到扰动信号。

合理设置STA/LTA触发检测参数是检测扰动信号的关键。由于触发检测算法本身的局限性,会不可避免地造成误触发和漏触发,需要运行维护人员掌握触发检测算法原理及相应参数的意义,并结合钻孔应变观测台站的实际观测情况不断调整STA/LTA触发检测参数,以达到最优的检测效果。

不同形变台站观测到的扰动信号的分布形态、持续时间、振动频率和能量分布都可能存在差异,需要提取小波分析的不同层来进行触发检测,也可选择多层分别进行触发检测。如果多个台站观测到类似的扰动信号,在排除自然界和人工干扰后,可判断为前驱波。

致谢：感谢童汪练研究员提供帮助!