中强震孕育过程中通常伴随着比较明显的应力变化,监测地震前后应力变化过程和规律是深入了解和研究地震的一种重要手段(陈颙等,2009),而直接测量深部介质的应力难以实现,往往通过利用重复震源测量介质波速变化的方式代替(杨润海等,2011)。气枪震源是一种目前应用较广的重复震源,与传统人工震源相比,它既不会污染和破坏生态环境,又能传播较远的距离; 与背景噪声等被动源相比,有明显的时间分辨率优势; 与天然地震相比,有人工可控、震源已知等优点,因此被认为是目前探索地球介质结构和介质变化的最佳人工震源(陈颙等,2007; Chen et al,2008,2017; Wang et al,2016)。气枪子波分为压力脉冲和气泡脉冲,压力脉冲频率高,传播距离短,气枪信号的远距离传播得益于其低频的气泡脉冲,但其能量相对于压力脉冲较小,因此较远台站接收信号信噪比较低,甚至难以辨认。为了提高信噪比,通常会对信号进行滤波和叠加等处理,但部分台站由于随机噪声干扰太大,处理后的信号信噪比仍然较低。而根据克拉默—拉奥下界(Cramer-Rao Lower Bound)(Carter,1987; 王伟涛,2009; 杨微,2013),在波速变化的计算过程中,信噪比是最重要的影响因素之一,因此在计算气枪信号波速变化(或走时差变化)时,较低的信噪比会使计算结果出现较大的误差,甚至不能反映介质的真实变化情况。

二次相关时延检测方法(唐娟,行鸿彦,2007)是一种简单易实现的运算方法,利用有效信号与噪声之间、噪声与噪声之间的弱相关特征,将自相关和互相关相结合,可在较低的信噪比信号中获得较高的时延估计精度。目前,该方法在声信号时延估计、闪电辐射源定位及卫星干扰源定位等领域有较为广泛的应用(杜鹃,程擂,2010; 周康辉等,2013; 窦慧晶等,2016)。本文将该方法应用于地震信号的时延检测,首先利用模拟信号检测二次相关时延检测方法应用于地震信号的可行性,随后利用该方法计算宾川较远台站气枪信号的走时变化。

相关函数是信号间相似度的一种表现形式。自相关函数表示相同信号间的相关程度,互相关函数表示不同信号间的相关程度。

互相关函数在背景噪声中已经有相当广泛的应用(Aki,1957; Shapiro,Campillo,2004),背景噪声中隐藏着微弱的周期性地脉动信号,通过互相关创建经验格林函数可以将背景噪声与地下结构联系起来(Friedrich et al,1998; Hasselmann,1963; Longuet-Higgins,1950)。地震背景噪声周期较长、振幅较小,单次互相关所得信号信噪比低(王伟涛等,2011)。本文处理信号的时间尺度较小,噪声可认为是随机噪声,并假设有效信号与噪声之间、噪声与噪声之间不相关。因此经互相关处理后,有效信号与噪声之间、噪声与噪声之间互相关函数值被压制,有效信号之间函数值被放大。

将两道信号f₁(t)和f₂(t)分别定义为:

f₁(t)=x(t)+n₁(t)(1)

f₂(t)=x(t+D)+n₂(t)(2)

式中:x(t)为信号源辐射的有效信号; D为两次信号的相对延迟时间; n₁(t)和n₂(t)分别为两道信号中的加性高斯白噪声。

单道信号f₁(t)的自相关函数表示为:

R₁₁(τ)=E[f₁(t)f₁(t+τ)]=R_xx(τ)+R_xn1(τ)+

R_n1x(τ)+R_n1n1(τ)(3)

式中:R_xx为有效信号之间的互相关函数; R_xn1和R_n1x为有效信号与噪声之间的互相关函数; R_n1n1为噪声之间的互相关函数。

假设噪声与信号之间不相关,式(3)可以简化为:

R₁₁(τ)=R_xx(τ)+R_n1n1(τ)(4)

两道信号f₁(t)和f₂(t)的互相关函数为:

R₁₂(τ)=E[f₁(t)f₂(t+τ)]=R_xx(τ+D)+

R_xn2(τ)+R_n1x(τ+D)+R_n1n2(τ)(5)

同样忽略有效信号与噪声之间的互相关函数,可以将式(5)简化为:

R₁₂(τ)=R_xx(τ+D)+R_n1n2(τ)(6)

此时,我们可以将自相关函数R₁₁和互相关函数R₁₂看成两道新的关于时间t的信号,对两道新的信号进行互相关运算,得到二次相关函数:

R_RR(τ)=E[R₁₁(t)R₁₂(t+τ)]=R_RX(τ+D)+E[R_xx(τ)R_n1n2(t+τ)]+E[R_n1n2(t)R_xx(t+D+τ)]+R_RN(τ)(7)

忽略信号与噪声之间的二次相关函数,式(7)可简化为:

R_RR(τ)=R_RX(τ+D)+R_RN(τ)(8)

式中:R_RX为有效信号的二次相关函数; R_RN为噪声的二次相关函数。

如果加性噪声之间是不相关的,则忽略噪声的二次相关函数,此时二次相关函数变为:

R_RR(τ)=R_RX(τ+D)(9)

由式(9)可知,当时间τ=-D时,二次相关函数取得最大值R_RR(-D)。通过对不同信号二次相关函数最大值对应时间的拾取,可以获取不同信号间的时间延迟变化。前述假设噪声与信号之间不相关,所以噪声与信号之间的互相关函数可以忽略。实际情况下,噪声并不一定是高斯白噪声,噪声与信号之间的互相关函数并不完全等于0(王彬,2009),但与二次相关运算之前的信号相比,噪声相对于有效信号的幅度要小很多(唐娟,行鸿彦,2007)。

二次相关时延检测方法在低信噪比信号的基础上构建了信噪比更高的信号,为了验证该方法应用于地震信号波速变化检测的可行性,本文首先对理论地震数据进行二次相关检测实验。波速变化与走时变化存在如下关系:

δv/v=-δt/t(10)

式中:δv,δt为相对波速、走时变化; v,t为原始路径的地震波速度和走时(Schaff,Beroza,2004)。

本文通过测量走时变化来代替波速变化。已知的理论信号长度为1 400个采样点,信号带宽为2～8 Hz,采样率100 Hz,利用移动窗口压缩—拉伸法(Meier et al,2010; 刘志坤,2010)生成的波速变化率不超过±5‰,且具有2个正弦周期的1 000道信号,然后在每道信号上随机加入高斯白噪声,最后将生成的信噪比为3 dB的1 000道信号进行线性叠加作为参考模板信号c_f,如图1a所示。图中黑色波形为叠加的参考信号,其他5个波形为1 000道信号中随机选取的加入噪声和波速变化信号。从图1a中可以看出,加入噪声的信号相对叠加信号已经基本看不出信号的有效相位。

对于加入噪声的1 000道信号,我们从两个方面进行比较:直接计算走时变化和经二次相关后计算走时变化。直接计算走时变化时,利用基于移动窗口互相关的时延检测方法(罗桂纯等,2008; Wang et al,2008; 刘志坤,2010; 刘自凤等,2015),将这1 000道信号与叠加信号c_f进行互相关时延检测,然后选取大多数信号走时变化比较平稳的相位提取走时变化。二次相关及计算走时变化的处理步骤如下:①对叠加的参考信号c_f进行自相关运算,得到自相关函数c₀; ②将加入噪声的每一道信号分别与参考信号c_f进行互相关运算,得到1 000条互相关函数C₁; ③将C₁里的每一条互相关函数与c₀进行互相关运算,得到1 000条二次相关函数C₂; ④将C₂里的所有二次相关函数进行线性叠加,得到二次相关的参考函数c_s; ⑤将C₂里所有二次相关函数分别与c_s进行互相关时延检测,选取合适的相位提取走时变化,并对检测结果进行插值以提高检测精度。二次相关方法本身具有时延检测的功能,但考虑到精度方面的要求,本文利用互相关时延检测方法对二次相关后的函数进行走时变化的提取。图1b为经二次相关后的信号波形,与图1a中原始信号波形相比,有效信号部分被凸显。

图1 模拟地震信号(a)及二次相关后的信号(b)
Fig.1 The simulated earthquake signals(a)and their second correlation signals(b)

用以上两种方法计算的走时变化和理论波速变化率对比结果如图2所示。从图中可以发现,直接计算的走时变化因为受到低信噪比等因素的影响,与设定的波速变化相比基本没有趋势变化(图2a); 而二次相关后的走时变化趋势同设定的波速变化趋势较为吻合(图2b)。

从理论地震信号的二次相关时延检测的实验可以发现,自相关和互相关运算的结合突出了有效信号而压制了不相干噪声,从而提高了信号的信噪比和后续时延检测的精度。在实际数据处理中,不同台站信号的噪声可能存在差异,但由于噪声是随机的,根据式(1)～(9),相对于有效信号相同的相位和周期,随机噪声二次相关后信噪比明显降低。而不同信号的高斯白噪声同样是随机的,因此不同信号间高斯白噪声的互相关函数值较小。

宾川气枪信号发射台于2011年在云南大理宾川县大银甸水库建成,位于红河断裂和程海—宾

图2 直接计算(a)和二次相关后(b)的走时变化结果与理论波速变化率对比
Fig.2 Comparison between the theoretical velocity variation rate and the travel-time delay variation by directly calculating(a)and the second correlation(b)

川断裂之间,该地区地震活动性明显(杨润海等,2014; 陈佳等,2017),气枪发射台气枪阵列由4支Bolt公司生产的1 500 LL气枪组成(王彬等,2015,2016)。通常情况下,气枪阵列的工作气压为15 MPa,沉放深度为10 m,单次激发能量相当于一次M_L0.7天然地震(杨微,2013)。气枪激发信号由分布在发射台周围的40个流动观测台站及原有的固定接收台站组成的接收台网接收,最远的流动接收台站为距离气枪发射台151 km的53252台,最近的为发射基地室内的CKT0台,距离气枪阵列约40 m,CKT0台被用作参考台,接收信号被视为震源函数。宾川气枪发射台激发信号最远可以传播300多千米,单次激发信号可以被53252台记录到,但震中距大于20 km台站的记录信号信噪比普遍较低(陈蒙,2014; 向涯等,2017)。

二次相关方法应用于模拟地震信号具有较好的效果,本文将其应用于气枪信号走时变化实际检测中,选择距离气枪震源67 km的53036台记录信号的SHZ分量进行走时变化的计算。原始记录信号中往往携带水位变化等震源信息或数据切取过程中的人为误差等干扰信息,因此在进行二次相关之前需要对53036台信号进行反褶积运算(王宝善等,2012; 翟秋实等,2016)。具体处理步骤如下:①53036台信号时频分析,确定滤波范围,并对53036台原始信号滤波以提升信号信噪比,剔除部分信噪比低或者难以辨认的信号; ②53036台数据对CKT0记录信号进行反褶积运算,得到近似格林函数; ③单天所有的格林函数线性叠加以提升信噪比,得到这天的一条格林函数g_d,再次进行滤波,并将所有格林函数进行叠加,得到参考格林函数g₀; ④参考格林函数g₀进行自相关运算得到自相关函数c₀,同时每天的格林函数g_d与参考格林函数g₀进行互相关运算得到每天的互相关函数c_d,再将c_d与自相关函数c₀进行互相关运算,得到每天的二次相关函数c_2d; ⑤利用互相关时延检测技术分别对格林函数以及二次相关函数进行走时变化计算,并对检测出的走时变化结果进行插值以提高检测精度。

不同深度介质的速度变化敏感程度可能存在差异(赵盼盼等,2012),波速变化的大小可能不同,因此为了计算不同相位的走时变化,对于不同时窗内的信号分别计算二次相关函数。图3为53036台的格林函数和不同相位的信号及其二次相关函数,其中图3a为选取的全部格林函数,黑实线为叠加的格林函数,图3b-1为所有格林函数的12.7～14 s信号(相当于P波信号),图3b-2为相对应的二次相关函数,图3c-1为27～27.5 s信号(相当于S波信号),图3c-2为对应的二次相关函数。图3中,由于53036台震中距较大,格林函数存在较多的干扰。选取的两段信号中,12.7～14 s内格林函数振幅较大、信噪比较高,相位之间同步性较好,二次相关函数与原始格林函数相比信号的同步性更好、信噪比更高; 27～27.5 s内信号振幅相对较小、信噪比较低,相位同步性较12.7～14 s内的信号差,

图3 53036台格林函数(a),选取的12.7～14 s格林函数(b-1)及其二次相关函数(b-2),27～27.5 s格林函数(c-1)及其二次相关函数(c-2)
Fig.3 The Green functions of station 53036(a),the selected signal of 12.7～14 s(b-1)and its corresponding second correlation functions(b-2),Green functions of the selected signal of 27～27.5 s(c-1)and its corresponding second correlation functions(c-2)

二次相关同样提升了信号的信噪比。

分别计算选取时间段内格林函数和二次相关函数的走时变化,结果如图4所示,图4a为12.7～14 s内信号走时变化,直接计算和二次相关后计算的走时变化结果趋势相近,但二次相关后计算的走时变化误差相对更小,变化更加稳定。图4b为27～27.5 s内信号走时变化,两种方法计算结果显示的总体变化趋势存在一定差异,但在某些时段的变化趋势相同,这可能跟不同时段噪声水平存在差异有关,且二次相关后的计算结果更加稳定。图4c为直接计算的走时变化对比,12.7～14 s内信号的走时变化总体趋势逐渐减小,即波速逐渐增大; 27～27.5 s信号的走时变化总体趋势是逐渐上升,即波速逐渐减小,两段信号的变化趋势存在差异。图4d为经二次相关后计算的走时变化对比,12.7～14 s信号的走时变化趋势为逐渐减小; 27～27.5 s信号的走时变化受信噪比影响波动较大,但总体变化趋势也是逐渐减小。

二次相关时延检测方法通过一次自相关和两次互相关运算构建新的格林函数,压制了理论地震信号的噪声、加强了有效信号,将该方法应用到气枪信号中有较好的效果。与直接计算的走时变化结果相比,理论地震信号二次相关后得到与设定波速变化较吻合的走时变化,说明二次相关时延检测方法具有一定的可靠性。但同时发现,二次相关后的走时变化与设定的波速变化相比,特别是在500道信号之前,有较为明显的相位差,这可能主要与信号较低的信噪比有关。信号的信噪比为3 dB,从图1中可以看出信号的有效相位

图4 53036台格林函数走时变化检测结果
Fig.4 The results of travel-time delay variation of station 53036's Green functions

难以辨别,二次相关可在一定程度上提升信号的信噪比和运算结果的可靠性,但受限于信号的低信噪比,还未达到最理想的效果,因而出现一定的相位差。

将二次相关应用到计算气枪信号走时变化时,气枪信号二次相关后计算的走时变化更加稳定。信噪比较高的12.7～14 s信号,直接计算的和二次相关后计算的走时变化趋势相同; 信噪比较低的27～27.5 s信号,两种方法的计算结果变化趋势不同,但二次相关的计算结果与12.7～14 s信号的结果趋势相同,只是原始信号信噪比较低,计算结果波动较大。在应用互相关时延检测方法计算气枪信号走时变化时,利用Cramer-Rao Lower Bound计算走时变化的误差下限,方法如下:

σ_τ≥σ^m_τ=(3/(2π²f ³₀T(B³+12B))[1/(ρ²)(1+1/(SNR²))²-1])^1/2(11)

式中:σ_τ为实际算得的标准差; σ^m τ为理论标准差; f₀为信号主要频率; T为互相关窗口的长度; B为频宽比; ρ为待检测信号与模板信号间的相关系数; SNR为信号的信噪比。12.7～14 s信号和27～27.5 s信号的二次相关函数的主要频率分别为5 Hz和4 Hz,频宽比分别为2和3,相关系数都约为1,信噪比分别约为300和90,相关窗口的长度都为0.1 s,分别带入式(10)求得两段信号的理论标准差约为6.50×10^-5和2.16×10^-4,计算结果的标准差均值分别为1.58×10^-4和1.30×10^-3,2段信号二次相关的计算结果误差都大于理论误差下限,因此本文气枪信号走时差测量计算结果较为可靠。

二次相关利用噪声与有效信号之间及噪声与噪声之间的弱相关特征,将自相关与互相关运算相结合。本文通过模拟分析和实际资料处理对比,发现二次相关可在滤波、叠加的基础上进一步提升信号的信噪比,经过二次相关处理后计算得到的走时变化也更加稳定可靠。

本文计算走时变化所用互相关时延检测程序由中国科学技术大学王宝善教授提供,在此表示衷心感谢。