傅里叶变换是一种常用的信号频率的分析方

法,但傅里叶变换只能反映给定信号的整体频率特征。S变换是由Stockwell等(1996)提出来的一种无损可逆的加时窗傅里叶变换,变换过程中时窗宽度与频率呈反向变化,改善了傅里叶变换不能反映频率随时间变化的不足,能够很好的进行信号时频局部分析(李玲莉等,2012; 郑成龙,王宝善,2015)。信号h(t)的S变换表示为:

S(τ,f)=∫^+∞∫_-∞h(t)(〖JB<1|〗f〖JB>1|〗)/((2π)^1/2)exp[(-(τ-t)² f ²)/2] exp(-i2πft)dt(1)

式中:h(t)是台站记录到的时域信号,S(τ,f)是经S变换后的结果,f为瞬时频率,τ为时间窗在信号上的中点,i为虚数单位。高斯窗函数定义为:

g_f(t)=(〖JB<1|〗f〖JB>1|〗)/((2π)^1/2)exp(-(t² f ²)/2)(2)

基本小波定义为:

w_f(t)=(〖JB<1|〗f〖JB>1|〗)/((2π)^1/2)exp(-(t² f ²)/2-i2πft))

=g_f(t)exp(-i2πft)(3)

为对比研究气枪主动源地震波信号与同当量的天然地震波信号,我们选取北京时间2013年

图1 53274台记录的ML0.7天然地震信号(a)和单次气枪激发信号(b)
Fig.1 Waveform of the ML0.7natural earthquake(a)and single air-gun excitation signals(b)recorded by the Station 53274

11月22日10时37分13.6秒发生在云南省宾川县(100.42°E,25.73°N)的一次M_L0.7地震作为对比研究对象。此次地震震源深度5 km,距气枪主动源发射台11 km。53274台距气枪源8.8 km、距天然地震的震中距7.7 km,我们用53274台记录的地震信号和气枪信号进行对比研究。在53274台的记录数据中,分别截取此次地震的记录和一次气枪激发(激发时间:2016年2月22日7时14分59秒; 沉放深度:10 m; 工作压力:15 MPa)信号(图1),应用S变换分析2种地震信号的时频特征。

将53274台接收到的气枪激发信号和天然地震信号的东西(SHE)、南北(SHN)和垂直(SHZ)三分量分别进行S变换,得到2种地震信号三分量的时频特征。如图2所示,通过S变换,可以很清楚地分析各地震信号频率随时间的变化情况。频率特征上,天然地震信号在传播7.7 km时三分量的频率范围为1～40 Hz,主要频率为10～20 Hz; 气枪信号在传播8.8 km时三分量的频率范围为1～10 Hz,主要频率为3～6 Hz。在信号震动时间上,天然地震持续了约6 s,而气枪震源持续了约9 s。在信号振幅上,天然地震信号SHE、SHN、SHZ三分量最大振幅分别为2.098 5×10⁴ count、3.328 0×10⁴ count、3.318 6×10⁴ count,气枪信号三分量最大振幅分别为2.061 0×10³ count、1.632 3×10³ count、2.216 2×10⁴ count,天然地震波的振幅比气枪震源地震波的振幅大了一个数量级。综合对比分析可知,天然地震信号相对气枪信号振幅大、频带宽,震动时间短,能量主要集中在高频部分。

图2 53274台天然地震信号时频图(a)和气枪信号时频图(b)
Fig.2 Time-frequency spectrum of natural earthquakes(a)and single air-gun excitation signals(b)of the Station 53274

为了得到天然地震信号和气枪激发信号的衰减特性,我们对比了2种地震信号不同震中距Z分量的最大振幅值。图3展示了2种地震信号的最大振幅衰减特征,其中实心圆和正方形分别表示天然地震信号和气枪激发信号幅值,实线和虚线分别是根据各振幅散点拟合的曲线。可见,天然地震信号的振幅衰减速度快于气枪激发信号。

分别对不同台站记录到的气枪信号进行2～7 Hz滤波,天然地震信号进行1～40 Hz滤波,得到2种震源激发的地震波在不同距离上的记录(图4)。图4a为气枪单次激发归一化信号,图4b为M_L0.7天然地震信号归一化信号。在已有的流动观测台中,距离气枪源151 km的台站可以接收到气枪单次激发信号,而同当量的天然地震信号仅能传播约60 km。由此可见低频的气枪信号可以传播得更远,同当量天然地震信号则在传播过程中衰减较快,传播距离要短得多,且由于气枪信号传播距离较远,震相发育也更丰富。

图3 天然地震信号和气枪信号的最大振幅衰减特性
Fig.3 Characteristics of the maximum amplitude attenuation between natural earthquakes and single air-gun excitation signals

图4 单次气枪激发信号(a)和天然地震信号(b)传播对比
Fig.4 Comparison of the signals of air-gun excitation(a)and natural earthquakes(b)

聚类分析在天然地震方面已经有很多的应用(刁桂苓等,1992; 高朝军等,2016),本文选用了层次聚类分析方法进行研究(王伟涛,王宝善,2012)。聚类过程中首先选取合适长度的信号窗口,然后计算N个样本中任意两次信号j和k的相关系数c_jk。定义这两次信号的相似度距离为:

d_jk=1-c_jk(j,k≤N)(4)

式中:d_jk越小,两次信号越相似。所有信号都与其自身以及其他信号进行互相关运算,互相关运算结果构建一个N×N的相似度矩阵M,矩阵中有M_jk=M_kj和M_j=k=1,所以聚类过程中只取上三角或者下三角矩阵即可。聚类分析之前,数据类别为N,要分为多少种类别是未知的,聚类时首先将相似度距离最小的1类信号聚为1类,此时数据类别减少为N-1,如此循环直到所有信号归为一类。两类信号G_p和G_q计算相似度距离时采用最长距离法,即是两类信号中信号的最大距离为两类信号的距离:

D_pq=max(d_pq,p∈G_p,q∈G_q)(5)

之后根据需要,给定阈值α,就可以将相似度距离小于或者等于α的信号聚成一类。给定的阈值不同,聚成的类别就有可能不同。聚类分析得到的结果,提高了叠加信号的相关性,使用叠加方法提高信号的信噪比更为有效。

取参考台记录的宾川主动源发射台2012年9月至2016年4月底激发的共5 311次信号的Z分量进行分析。由前述可知,气枪信号主频约为3～6 Hz,聚类前先将参考台信号进行去均值、去线性趋势、2～7 Hz滤波处理。经过互相关运算得到相似度矩阵M后,给定阈值α= 0.1,得到聚类结果如图5所示,其中纵坐标表示聚成类别数量,横坐标表示信号的时间,信号上方的数字表示该类信号的数量。从图中可以看出,5311次信号被分成18类,信号数量最多的是第4类,共有2 125次。

图5 2012年9月至2016年4月CKT0台信号聚类结果
Fig.5 The results of the hierarchical clustering of the CKT0's signals from September 2012 to April 2016

线性叠加就是将需要叠加的数据先进行算术相加,然后求其平均值,方法如下:

s(t)=1/N∑^N_j=1s_j(t)(6)

式中:s(t)表示线性叠加后的结果,N表示需要叠加的信号数量,s_j(t)表示第j道信号。经n次叠加后,地震信号强度增强n倍,随机噪声信号强度增强n^1/2倍,因此地震信号强度增长速度是随机噪音的n^1/2倍(熊章强等,2010)。

相位加权叠加(Schimmel,Paulssen,1997; 陆一锋等,2011)是一种非线性叠加方法,其中心思想是在线性叠加结果的基础上,信号的每一个采样点乘以一个对应的权重值,从而得到相位加权叠加结果。相位加权叠加的实现,首先用气枪信号与其希尔伯特变换结果构造复数序列以得到瞬时相位,再将线性叠加结果与相位计算的权重值对应相乘,方法如下:

X(t)=x(t)+iH[x(t)]=A(t)exp[iΦ(t)](7)

s(t)=1/N ∑^N_j=1x_j(t)〖JB<2|〗1/N∑^N_j=1exp(iΦ_j(t))〖JB>2|〗^v(8)

其中:(7)式是处理单道信号得到瞬时相位的理论公式,x(t)是台站记录到的单次信号,H[x(t)]是单次信号的希尔伯特变换,A(t)是X(t)的瞬时振幅,i为虚数单位,Φ(t)是瞬时相位;(8)式是处理所有信号得到叠加结果的理论公式,其中s(t)是相位加权叠加的结果,x_j(t)是台站记录到的第j次信号,Φ_j(t)是第j次信号经过(7)式处理过的瞬时相位,v是指数,v=0时是线性叠加,本文中v=2。噪声的瞬时相位随机不定,而气枪信号的瞬时相位基本一致,因此相位加权叠加是一种削弱噪声、增强气枪信号的方法。

选聚类结果中信号最多的第4类数据进行后续的叠加。其他台站选取相同激发时刻的数据,由于并不是所有台站都记录到了全部的激发数据,因此各个台站用于叠加的数据也不尽相同,叠加次数在1 900～2 125之间。在叠加前各台站也要进行去均值、去线性趋势、滤波处理,叠加后再将结果进行滤波、波形灭尖等处理。叠加结果如图6所示,图6a为线性叠加的信号,图6b为相位加权叠加的信号。经过约2 000次叠加后,线性叠加的结果在225 km处信号的震相还比较清晰,更远距离信号的震相已经不大清晰; 相位加权叠加后,在更远的350 km处台站信号的震相清晰可见。