目前最常使用的地震触发检测算法是短时绝对均值与长时绝对均值之比(即STA/LTA)(彼得·鲍曼,2006)。STA测量信号的短时间变化,随时监视着地震; 而LTA测量台站当前地震噪声的振幅水平。

(1)STA/LTA算法原理

首先分别计算移动的STA和LTA时间窗中振幅的绝对均值,再计算这两个均值的比值(STA/LTA)。程序将此比值与用户设定的触发阈值相比较。若大于触发阈值时,就认为检测到地震并触发报警。 当STA/LTA比值低于解除触发阈值时,就认为地震信号结束并解除报警。STA/LTA解除触发阈值要小于STA/LTA触发阈值(图1)。

STA、LTA都作为滑动均值进行计算:

STA_i=STA_i-1+(|x_i|-STA_i-1)/(N_STA),(1)

LTA_i=LTA_i-1+(|x_i|-LTA_i-1)/(N_LTA),(2)

R_i=(STA_i)/(LTA_i).(3)

式中:x_i是地震信号(经滤波或未经滤波); STA是短时均值; LTA是长时均值; N_STA和N_LTA分别是STA窗和LTA窗中的点数。

图1 STA/LTA触发算法原理
(a)一个连续的(经过滤波)的地震信号;(b)在STA、
LTA两个窗口里平均绝对信号随时间的变化;
(c)STA/LTA随时间的变化
Fig.1 The principle of STA/LTA trigger algorithm (a)a continuous(filtered)seismic signal;(b)mean absolute signals in STA and LTA windows change with time;(c)STA/LTA changes with time

触发参数设置是否合适对能否成功检测地震事件很关键。STA/LTA的基本触发参数有:STA时间窗长、LTA时间窗长、STA/LTA触发阈值、STA/LTA解除触发阈值。对于仿真成短周期的信号,STA用来测量地震信号的“瞬时”振幅,STA窗长一般为0.5 s; LTA测量台站当前地震噪声的振幅水平,通常LTA窗长为50～500 s。安静的台站,STA/LTA触发阈值常设置为4; 人为噪声强的台站,触发阈值可取8以上。STA/LTA解除触发阈值决定何时解除报警,能记录弱地震活动的安静台站,STA/LTA解除触发阈值的取值是2～3; 对噪声大些的台站,该阈值应设置高一点。

(3)改进STA/LTA地震检测算法

为提高地震检测灵敏度,并有效减少地震的漏触发、虚触发和重复触发,本研究采用冻结LTA值和设置触发有效的最小持续时间、相临两次有效的地震检测的最小间隔时间来改进传统的STA/LTA地震检测算法。

冻结LTA:一旦检测到事件后就冻结LTA,从而使背景噪声参考水平不受地震事件信号的影响; 当事件检测结束时,再解冻LTA来跟踪背景噪声变化。这样STA/LTA可以突出地震信号的变化。

检测触发有效的最小持续时间:为了有效减少尖锐脉冲或爆破等引发的虚触发,根据此类脉冲持续时间短的特点,采取一旦达到触发阈值,检测触发有效的最小持续时间,比如持续时间小于10 s可认定为干扰,持续时间超过10 s就认定为地震,从而能排除绝大部分因干扰引发的报警。

相临两次有效的地震检测最小间隔时间:为了避免同一次地震因波形起伏大引起的重复触发,可采取设置相临两次有效地震的最小间隔时间,即后一次地震的触发时刻与前一次地震解除触发时刻之间的时差(比如为240 s),间隔小于最小间隔时间的触发就认定为同一次地震,间隔超过最小间隔时间的触发则认定为另一次地震。

台站观测仪器工作环境好、仪器运行正常,才能产出正常、连续、可靠的地震波形数据,这是提高地震检测准确性的前提和基础。如果因仪器工作环境恶劣以及管理不善而导致仪器长期运行不正常、产出断续甚至畸形的波形,那么即便采用好的地震检测算法也无法显著提高地震报警准确率。

(1)台址选择

台址的选择,对地震波形记录质量尤为关键。必须按《地震台站观测规范》(国家地震局,1990)、《数字地震观测技术》(中国地震局监测预报司,2003)和《地震学与地震观测》(中国地震局监测预报司,2007)来选择。选址的注意事项主要有以下几条:测震台的选址,必须远离各种振动干扰源; 台基应选择在岩性坚硬致密、完整的、并大面积出露的基岩上; 台址的地势起伏要小,应避开风口; 台站观测环境地噪声水平要符合要求; 摆房要防潮保温。

(2)改善仪器运行环境

大多数测震仪器要在干燥的环境中才能正常工作。对潮湿的观测山洞,可以采取综合措施来降低潮湿程度,如采取工程措施对洞室内壁作防水防潮处理; 用密闭的仪器罩罩住地震计,罩内放置能吸湿的硅胶干燥剂; 不定期用除湿机对山洞除湿。对于室外摆线的铺设或架设,要采取加装外套或穿管埋地铺设等防护措施,防止鼠类动物啃咬线路。对于雨季雷患严重的地区,要制作良好的接地、电源防雷和GPS天线防雷,防止雷电串入干扰甚至烧毁设备。稳定的供电也是保证仪器正常运行的一个关键因素。

实践表明,如果将数字测震仪波形仿真成模拟短周期测震仪波形,就能滤除干扰的主要成分,突出每个地震的纵波,从而把各个地震区分出来。本研究采用无延时的实时仿真算法来仿真原始记录波形,原理如下:

结合宽频带数字测震仪的采样频率和模拟测震仪的传递函数,使用双线性变换法(王洪体等,2006; 李万金,邓存华,2011)将模拟测震仪器的传递函数滤波器变换为不同采样率下的数字滤波器传递函数,从而构造出递归仿真滤波器; 再用滤波的方法便可将数字地震波形仿真成模拟仪器记录波形,从而实现实时波形仿真。

由于模拟地震仪记录的是地动位移信号,而数字地震仪记录的是地动速度或地动加速度信号,因此需要将数字波形信号先转换成位移量。要将速度型数字地震仪波形转换成位移量,需要积分一次。

一般的做法是:先对原始数字地震记录数据进行积分,再用模拟仪器的传递函数来仿真。本研究的做法是:将速度型数字地震仪波形转换成位移量,需要积分一次,只需要将模拟地震仪的传递函数乘1/S后得到的传递函数来仿真; 将加速度型数字地震仪波形转换成位移量,要积分两次,只需要将模拟地震仪的传递函数乘1/S²后得到的传递函数来仿真 。这样可以大幅提高波形仿真的速度。

将减少了零点个数的模拟地震仪器传递函数变换为不同采样率下的数字滤波器传递函数,从而构造出递归仿真滤波器,可实现实时波形仿真。

使用MATLAB软件构造递归仿真滤波器(万永革,2007)的步骤:

假设模拟仪器的传递函数形式为 H(s)=f(s)/g(s),则

第一步,使用MATLAB函数f=poly(z)*k和g=poly(p)便可将仪器剩下的零极点转换成传递函数形式。

第二步,使用MATLAB双线性变换函数[a,b]=bilinear(f,g,F_s)将模拟仪器传递函数离散成数字滤波器。式中:f、g为模拟滤波器传递函数分子和分母多项式向量系数; F_s为采样频率,单位为Hz; a、b为数字滤波器传递函数分子和分母多项式向量系数。

如:将DD-1型模拟地震仪传递函数转换成采样率为100 Hz的速度信号的数字仿真滤波器:

a=0.2440 -0.8741 0.8128 0.5086 -1.2196 0.4672 0.1628 -0.1,

b=1.0000 -5.3974 12.2810 -15.3370 11.4941 -5.3084 1.4621 -0.1944.

第三步,将构造出的递归仿真滤波器系数代入下式的滤波算法,即可对波形进行实时仿真:

y(k)=1/(α₀)*(∑ⁿ_i=0b_ix(k-i)-∑ⁿ_i=1a_iy(k-i)).(4)

式中:系数a_i(i=0,1,…,n)、 b_i(i=0,1,…,n), n为递归滤波器阶数; x(k)、 y(k)分别为输入的数字测震仪采集数据和输出的DD-1型模拟短周期测震仪仿真波形数据。

运用上述仿真算法将波形仿真成DD-1型模拟短周期测震仪波形,主要有以下作用:

(1)滤除仪器不正常的波形

测震仪器因放大电路板工作受潮、接地不良或雷电感应等原因引起波形大幅振荡或毛刺脉冲,仿真成短周期波形后可以显著的滤除这些干扰,从而减少地震虚报。

(2)消除仪器因零漂产生的偏移

利用零漂周期长的特点在仿真成短周期波形过程中就能有效地消除数据中的仪器零漂,突出地震信号。

(3)滤除长周期大幅度干扰

每年台风季节,宽频带地震仪的波形上常记录到卓越周期为3～8 s的大幅度台风干扰。若不滤除,STA/LTA经常检测不到地震,导致地震漏报。

(4)从强震重叠波形中识别出余震

强震波形会湮没其中的余震波形,仿真成短周期波形后可以显露出这些余震,从而减少地震漏报。

地震仿真实例:以个旧地震台记录的2011年3月11日日本本州东海岸附近海域9.0级大地震(震中距为36.6°)为例,采用MATLAB编程对地震信号进行仿真以及地震检测判别,从仿真后的波形中可以发现主震之后又发生了多次余震(图2)。

图2 2011年3月11日日本9.0级大地震波形
(a)CTS-1仪器记录的原始地震波形;
(b)仿真成DD-1测震仪的波形
Fig.2 Japan MS9.0 earthquake on Mar.11,2011(a)original seismic waveform recorded by CTS-1 seismograph;(b)simulatied to waveforms of DD-1 seismograph

对近年来个旧地震台记录到的一些地震波形和干扰波形采用上述算法进行重新检测,既要有效排除大部分干扰波形,又要能保证一定级别的地震报警。经过比较,个旧地震台的STA/LTA报警检测参数合理设置为:STA窗长为0.5 s,LTA窗长为50 s,触发阈值是4,解除触发阈值是2,检测有效最小持续时间为8 s,相临两次有效地震检测的最小间隔时间为180 s。 以图2b中的垂直向波形为例进行检测演示,如图3所示。

图3 2011年3月11日日本9.0级大地震的检测过程
(a)仿真后的DD-1波形;(b)STA;(c)LTA;
(d)STA/LTA;(e)地震事件触发及解除触发
Fig.3 The detection process of Japan MS9.0 earthquake on Mar.11,2011 (a)simulation waveforms of DD-1 seismograph;(b)STA; (c)LTA;(d)STA/LTA;(e)seismic event triggering

与采用传统算法(图1)对比,从STA/LTA曲线上可看出改进后的STA/LTA算法,当STA/LTA达到触发阈值后,冻结LTA能有效抑制地震信号对LTA的影响,从而提高STA/LTA值; 检测触发有效的最小持续时间参数能够有效地排除短时干扰; 相临两次有效的地震检测最小间隔时间参数能显著抑制因地震波形起伏较大造成对同一次地震的重复检测。

为了充分验证本算法的可行性和准确性,笔者从个旧地震台数字地震仪近年来记录到的地震中选取5个地震事件进行检验(表1),其中涵盖了近震、远震和极远震,波形中或有尖锐的爆破波形干扰、或有大周期台风振幅干扰、或余震完全被主震波形湮没,对台站需要检测的地震和波形所受的干扰具有显著的代表性。从表1的地震检测效果统计结果可看出,在仿真成短周期的地震波形上,使用改进后的STA/LTA触发检测算法来检测地震效果显著,只要地震间的间隔时间不要过短,完全可以实现一定级别地震的准确检测,并有效减少地震的虚报、漏报和重复触发,只是有时触发持续时间较短。

表1 地震检测效果统计
Tab.1 Statistic results of seismic detection effect