张 源 崔庆谷
1 云南省地震局,昆明市北辰大道148号,650224
大震前定点形变数据中周期为d的短持时变化可能预示着大地震的迫近,是短临预测的重要线索和依据[1]。与来源于局部区域的噪声和干扰相比,构造变动引发的暂态异常信号往往同步出现在同一台站的不同记录中,甚至准同步出现在相邻台站的多种记录数据中[2-3],这一特征为异常信号的识别创造了条件。本文在预处理观测数据的基础上,利用线性叠加及零延迟相乘算法对人工合成数据中的暂态信号进行模拟检测,对比2种算法在暂态信号识别中的效果,选出效果较好的算法用于实际观测数据;再利用同测点多道、多测点记录数据零延时相乘的方法有效压制干扰和噪声,使多台、多测项数据中准同步出现的暂态异常信号放大凸显,实现对构造变动信号的初步检测与识别。
噪声的种类复杂多样,压制噪声的方法也因噪声种类的不同而有所不同。对于符合高斯分布的随机白噪声,数据叠加是压制噪声的有效方法。时域线性叠加是提高信噪比常用的技术,将若干条记录时间对齐后直接叠加,再求平均,即为时域线性叠加[4]。其叠加公式为:
(1)
式中,x(t)为叠加输出结果,xj(t)为第j道的时间记录,n为记录数。
在常规线性叠加中,信号是相干的,叠加后信号的振幅和标准差基本不变;而叠加后噪声的振幅和标准差减小,故信噪比增加[5]。零延迟相乘是一种非线性运算,该运算对于数值大于1.0的相干信号有放大作用,且信号的强度越大,放大的倍数也越大;相反,数值小于1.0的相干信号将受到压制。
由于构造变动引发的暂态异常变化信号在同点、甚至相邻测点记录中具有相干性,如果将多道记录相乘,相乘后的信号为多个信号强度的倍乘。当此类暂态信号的幅度大于背景噪声强度时,相乘后得到的参考通道中暂态异常信号的放大倍数将远大于背景噪声和干扰,暂态异常信号在背景噪声中将被凸显出来,该结论可用人工合成的噪声和信号来证实。用计算机生成3组不相干的随机白噪声信号n1(t)、n2(t)、n3(t),其元素的变化幅度在0~0.2之间,利用永胜台垂直摆NS向记录的一段暂态短持时信号经平滑后作为模拟暂态信号Transient(t),将Transient(t)分别与n1(t)、n2(t)、n3(t)叠加后得到d1(t)、d2(t)、d3(t),将d1(t)、d2(t)、d3(t)线性叠加和零延时相乘后的输出数据进行对比,具体见图1。
图1(a)为原始噪声n1(t),图1(b)为n1(t)与模拟暂态信号Transient(t)叠加后的记录数据d1(t),图1(c)为不同数据相乘后的对比结果。可以看出,随着相乘道数的增加,暂态信号的强度呈指数倍增长。由此可见,强度较大的暂态信号相乘后其强度增长较快。图1(d)为2组数据线性叠加及零延迟相乘后的数据对比,可以看出,数据零延迟相乘可以使强度较大的准同步暂态信号凸显,极大地增加了暂态信号的可识别程度。图1(e)为3组数据线性叠加和零延迟相乘的结果,可以看出,零延迟相乘能够更有效地增大准同步暂态信号的可识别程度。
年变及长周期漂移是前兆地球物理观测中最常见、强度最大的干扰信号,为避免该干扰信号对计算结果产生影响,本文采用高通滤波方法去除年变和长趋势漂移。高通滤波器的低频端截止频率设定为30 d,可保证周期大于30 d的低频长周期漂移被去除(定点形变观测数据的采样率为min),滤波器的Z域表达式为:
(2)
式中,分子项系数numZ=[0.999 9 -1.999 8 0.999 9],分母项系数denZ=[1.000 0 -1.999 8 0.999 8],其频率特性及滤波效果见图2。
图2 截止频率30 d的高通滤波器频率特性及滤波效果Fig.2 Performance and effect of high-pass filter with cut-off frequency of 30 days
图2(a)和2(b)展示了用于去除年变和趋势性漂移的高通滤波器幅频特性和相频特性,其低频端截止频率为30 d(3.858 0×10-7Hz)。将2020-02-01~2022-06-10永胜台水管倾斜仪NS向和EW向数据(图2(c)和2(d))输入上述滤波器,滤波后的数据见图2(e)和2(f)。滤波器消除了年变和趋势性变化,保留数星期、数天及更高频段信息,达到预期效果。
对于同一观测点的多套观测仪器,如果观测的物理量相同,构造变动引起的暂态短持时异常信号将会准同步出现在多套仪器记录数据中,且具有相似的形态。此时,利用2种记录的相同分向(同为EW向或NS向)数据零延迟相乘和相加,可以得到信噪比更高的参考通道。该参考通道中不相干噪声和干扰被压制,相干信号被进一步放大[6],参考通道的信噪比得到提高,有利于准同步暂态异常信号的识别。
永胜台洞体应变、钻孔应变2套仪器安装距离约为200 m,记录的物理量都是应变类,其中2020-12-20~2021-06-10的EW向记录数据如图3(a)和3(b)所示,预处理后的数据如图3(c)和3(d)所示,将预处理后的数据时间对齐后相乘,得到参考通道见图3(e)。可以看出,零延迟相乘能有效提高数据的信噪比,使暂态短持时异常信号凸显出来。
图3 用零延迟相乘算法检测2021-05-21漾濞MS6.4地震前永胜定点形变中的准同步暂态信号Fig.3 Detecting quasi-synchronous transient signals in crustal deformation records of Yongsheng station before Yangbi MS6.4 earthquake on May 21st, 2021 by zero-delay multiplication algorithm
地震是应力应变积累的结果,大震孕育过程中的构造变动一般不限于震源区,震中附近甚至距离较远但构造相关的多个台站往往也能够记录到震前的构造变动信号。当多个台站记录的构造变动信号在时间上有重叠时,利用数据零延迟相乘可以检测到多台准同步信号的时间重叠部分。
2002~2022年云南定点形变观测数据中,川滇菱形块体中部丽江-小金河断裂附近的丽江、永胜定点形变台站(2个台站相距约60 km)多次出现持续时间约1周的准同步信号,其后不到40 d,300 km内有MS≥5.0地震发生。图4(a)和4(b)分别为2013-02-01~12-31永胜台和丽江台洞体应变NS向数据,图4(c)和4(d)为滤波后的结果,其中08-12发生左贡芒康6.1、08-31及09-05分别发生香格里拉MS5.0、5.9地震,图4(e)为2组数据滤波后零延迟相乘的结果。图4(e)中可以看到清晰的暂态短持时信号,其信噪比远高于其他通道。利用零延迟相乘算法梳理2002~2022年云南全境的定点形变数据,共检出暂态短持时异常信号11次,这些信号前后1个月或持续期内,云南境内皆有MS≥5.0地震发生(表1)。
表1 11次暂态短持时异常及后续MS≥5.0地震
图4 利用零延迟相乘法检测2013年永胜、丽江定点形变数据中的准同步暂态信号Fig.4 Detecting quasi-synchronous transient signals in crustal deformation records of Yongsheng and Lijiang stations in 2013 by zero-delay multiplication algorithm
表1中11次暂态短持时信号是利用配对台站洞体应变NS向数据零延迟相乘得到参考通道,再从参考通道中识别出来的。零延迟相乘后得到的部分参考通道见图5,其中图5(a)~5(f)与表1部分时间段的数据对应,可以看出,多个MS5.0地震事件发生前,一些台站的数据中出现了明显的准同步暂态短持时信号。
与随机干扰不同,暂态短持时异常信号来源于真实的构造变动,这种信号能够被同台的不同仪器,甚至相邻台站的不同仪器准同步记录下来,利用零延迟相乘算法可以将同步或准同步出现在多组记录中的暂态短持时异常信号检测出来。数值模拟和观测实验都证实,零延迟相乘算法有利于暂态短持时信号的检测与识别。利用该方法梳理云南地区定点形变长期观测数据,得到11组暂态短持时信号。统计结果表明,这11组暂态短持时信号与后续的中强以上地震在时间和空间上存在关联性,具体现象有待深入研究。
另一方面,在云南多个定点形变数据中,相邻台站零延迟相乘检测出准同步暂态信号的台站主要集中在丽江-小金河断裂附近的丽江、永胜和云龙3个台站,这不仅与台站的距离较近有关,同时可能与丽江-小金河断裂的特殊性有一定的关系。丽江-小金河断裂从川滇菱形块体中部穿过,将川滇菱形块体分为北南2个性质截然不同的大地构造单元,是连接川滇活动块体东西边界的纽带[7]。丽江-小金河断裂北连龙门山断裂,南端与维西-乔后断裂交会,既是地形高程过渡带,也是布格重力异常梯度带,断裂两侧的Rayleigh面波相速度存在明显差异,在高原扩展中起到关键作用[8]。丽江-小金河断裂南段基本处于闭锁状态,而北段处于蠕滑状态[9],川滇菱形块体东边界的构造活动可能通过丽江-小金河北段传递到完全闭锁的丽江-小金河断裂南段,这可能是云龙、丽江、永胜台记录数据出现相关性的构造证据。