中国地震  2016, Vol. 32 Issue (2): 356-378
地学长江计划安徽实验中低频可控震源地震波信号提取方法评估
顾庙元1, 姚佳琪1, 张伟1,2, 陶知非3, 张汝杰3     
1. 中国科学技术大学地球与空间科学学院地震学与地球内部物理实验室, 合肥市金寨路96号 230026;
2. 蒙城地球物理国家野外科学观测研究站, 安徽亳州 233527;
3. 中国石油集团东方地球物理公司, 河北涿州 072750
摘要:可控震源的激发过程精确可控、可重复性高。本文评估了重复激发形成的多道记录的信号增强方法,包括叠加方法和有效信号直接分离方法。相比线性叠加,加权叠加充分利用了有效信号的一致性,能进一步压制噪声,提高叠加记录的信噪比。相比时间域相位加权叠加,Semblance加权叠加更加稳定。时频域相位加权叠加能在时频域压制随机分散分布的噪声信号,叠加信号更加清晰。本文中,采用标准SVD分离获得的信号同线性叠加结果类似。本文利用密集分布的临时台阵评价地学长江计划安徽实验中低频可控震源在庐江第2激发点的激发效果。对比不同的处理流程发现,先互相关检测信号再进行时频域相位加权叠加能够取得最佳的处理效果,通过200多次重复激发,可以在50km范围内提取到较清晰的有效信号。
关键词低频可控震源    信号检测    多道记录    信号增强    加权叠加    地震震相    
Evaluating different methods for retrieving seismic wave signals from low frequency vibrator in the Anhui experiment of "Yangtze River Geoscience Project"
Gu Miaoyuan1, Yao Jiaqi1, Zhang Wei1,2, Tao Zhifei3, Zhang Rujie3     
1. Laboratory of Seismology and Physics of Earth's Interior, School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China;
2. Mengcheng National Geophysical Observatory, Hefei 230026, China;
3. BGP Inc., China National Petroleum Corporation, Zhuozhou 072750, China
Abstract: The signal exciting process of vibrator is well controlled and resulting signals are of high repeatability. In recent years, low frequency vibrator, for which the lowest frequency could be 1.5Hz, has been developed and shows its potential for subsurface structure exploring in regional scale. Detecting and enhancing extremely weak signal is the key for processing vibrator seismic signal at long distance. We evaluate different signal enhancing methods for multichannel seismic records excited by repeating low frequency vibrator, including stacking methods and separating methods. Compared with linear stacking, weighted stacking utilizes coherence of signal in multichannel records to suppress noise and enhance the signal-to-noise ratio(SNR)of stacking result. Semblance weighted stacking is more stable and effective than time-domain phase weighted stacking. Compared with time domain weighted stacking, time-frequency domain phase weighted stacking is capable of protecting signal frequency-band and the stacked signal is clearer. Signal value decomposition(SVD)can separate signal from multichannel records. In this experiment, the result of linear standard SVD separation is similar to linear stacking result. During the Anhui experiment of "Yangtze River Geoscience Project", the low frequency vibrator excited repeatedly at two locations. We analyze different processing sequences after preprocess. The best sequence is first to detect signal by correlation for each excitation then enhance signal by time-frequency domain phase weighted stacking. Using repeated signal more than 200 times, clear Pg and Sg phase can be retrieved within 50km.
Key words: Low frequency vibrator     Signal detection     Multichannel records     Signal enhancement     Weighted Stacking     Seismic phase    
0 引言

相比于被动源探测,主动源探测技术(如深反射地震)能够提供更加精细的地壳结构。利用炸药震源进行深反射地震探测地壳速度结构已成为地质构造演化研究的重要手段(王海燕等,2010侯贺晟等,2012)。但大当量炸药震源对环境的污染大,费用高,可重复性低,在使用中受到越来越多的限制。

寻找更加经济、安全、精确可控的人工震源一直是主动源探测研究的重要内容。目前,气枪震源已经被广泛应用于海洋油气地震勘探。近年来,气枪震源激发技术被逐步移植到陆上区域尺度地下结构的探测中,包括海上气枪激发、陆地接收信号的海陆联测实验和陆地水体气枪激发实验(陈颙等,2007丘学林等,2007林建民等,2008)。

可控震源具有精确可控、重复性好、移动方便的特点,现已成为陆地上常用的非爆炸震源,如精密控制机械震源、液压式可控震源等。精密控制机械震源(Accurately Controlled Routinely Operated Signal System,简称ACROSS)由电机带动2个偏心轮相向转动向地面输出一个垂直向下的合力,通过精确控制电机的转动频率连续向地下发射调频信号,信号振幅随频率成2次方上升,有效频带较窄(王洪体等,2009杨微等,2013)。目前,国内外学者已将精密控制机械震源应用到地壳结构探测(Li et al,2013刘明辉等,2015)、地震断裂带结构动态监测(杨微等,2010)等研究中。液压式可控震源已在陆上油气地震勘探中得到广泛应用。近年来,地震勘探中对低频地震信号的需求不断增大。相对于高频地震信号,低频地震信号具有较强的穿透高速屏蔽层的能力,能够提高深层成像质量(Ziolkowski et al,2003佘德平等,2007),可靠的低频信息还能提高地震反演的准确程度(Baeten et al,2013)。常规可控震源激发频率较高(一般起始频率>6Hz)。高精度低频可控震源是近年来可控震源技术的最新发展成果(陶知非等,2010刘振武等,2013)。东方物探公司依托国家863计划开展了高精度可控震源的研制,推出的LFV3型可控震源是国际上唯一经过野外验证的工业化低频可控震源,能够实现1.5Hz低频信号的激发(陶知非等,2011),现已得到全面推广应用。其激发的低频地震数据已在全波形反演中得到成功应用(陶知非等,2011Baeten et al,2013)。目前,东方公司已完成EV56高精度可控震源的研制工作,其信号控制技术更加成熟可靠,激发信号畸变更低。

低频地震信号在传播过程中不易发生衰减和散射,有效传播距离远,因而低频可控震源在区域尺度地下结构探测中具有很大的潜力。地震信号经远距离传播后,能量非常微弱。极微弱信号的检测与增强是处理远距离传播可控震源数据的关键。针对多次重复激发产生的多道记录,本文评估了不同信号增强方法的处理效果,包括不同的叠加方法和使用SVD提取信号的方法,同时对不同处理顺序的影响也进行了评估。其后,还分析了背景噪声水平和合成地震记录,结果显示,在更远距离上接收有效信号需要合适的场地激发条件、环境噪声水平较低的地震台站以及可控震源重复激发或多台震源车同时激发等来增大激发能量。

1 低频可控震源与观测系统 1.1 实验简介

2015年10月,地学长江计划选取安徽作为先期实验区,实施人工震源地下结构探测科学实验(以下简称安徽实验)。船舶在长江航道马鞍山-安庆段激发气枪震源。通过附近的固定地震台和流动观测地震仪接收信号,以对长江两岸地下结构进行探测。

中国科学技术大学在庐江布置了2条垂直于郯庐断裂的密集测线和覆盖大别山东缘的扇面状台阵。郯庐断裂带南段2条主干断裂(池河-太湖断裂、嘉山-庐江断裂)穿过台阵布设区域(图 1)。

图 1 安徽实验中低频可控震源车激发位置和临时台阵示意图 三角为临时台站;五角星为低频震源车2次激发位置;黑色粗线为池河-太湖断裂、嘉山-庐江断裂

在气枪震源激发间隙,高精度可控震源车EV56在2个选定的激发点(图 1)进行激发实验(表 1)。在可控震源工作过程中,振动平板被紧紧压在地面上,重锤和平板的加速度引起平板对地面的压力变化。这种压力变化被称为振动输出力,可在地表激发弹性波向地下传播。额定振动输出力参数表征了可控震源理论上的最大激发能力,实验所用的EV56型可控震源车额定振动输出力可达56000磅。实际的振动输出力常小于额定振动输出力,震源输出力25%是指实际振动输出力能达到额定振动输出力的25%。

表 1 可控震源扫描参数

第1个激发点位于铜陵,在该点进行了2种重复激发。第1次在2015年10月8~9日进行,用时长20s,1.5~12.0Hz线性升频信号以35%的输出力重复激发了1207次。第2次在2015年10月9~10日进行,用时长20s,1.5~6.0Hz线性升频信号以25%输出力进行了502次重复激发。震源车在铜陵激发点激发效果不佳,有效信号仅在最近的铜陵地震台接收到(相距约5km)。该激发点位于一块由回填土填埋的空地上,浅层介质比较松散,因此,激发效果较差可能与激发点的局部地质条件有关。

第2激发点位于庐江,靠近中国科学技术大学布置的临时台阵(图 1),距最近的台站约6km,在75km内有密集临时台站分布。2015年10月13日,可控震源车用时长20s、1.5~8.0Hz线性升频信号以25%输出力进行了289次激发。本文利用密集分布的临时台站定量评价低频可控震源在庐江激发点的激发效果。

1.2 可控震源激发信号

可控震源在地面激发的是长时间的连续调频信号,通过与地面耦合的振动平板向下传播能量。与炸药震源不同,可控震源激发信号的瞬时能量密度较小,其有效激发能量通过长时间的累积得以增强。

1.2.1 扫描信号与输出力信号

由于机械系统之间、振动底板与地表介质之间的耦合作用,激发的调频信号会发生一定畸变,即传到地面的地面力信号不仅包含基波信号,还包含对应的高阶谐波。谐波的瞬时频率正好为对应基波瞬时频率的整数倍。在可控震源数据处理中,谐波信号通常被视为噪声。谐波的发育强度与激发条件有关(Wei et al,2011)。Stockwell变换是一种有效的时频分析方法,结合了短时傅里叶变换和小波变换的优点,具有良好的多尺度时频分辨能力(Stockwell et al,1996)。本文利用Stockwell变换分析震源激发信号的时频特征。图 2(a)为设计的可控震源线性扫频信号,其对应的参数见表 1图 2(b)为实际输出力信号,与设计的扫描信号(图 2(a))相比,实际输出力信号存在一定畸变。图 2(c)为扫频信号经Stockwell变换得到的时频振幅谱。图 2(d)图 2(b)输出力信号对应的时频振幅谱,由其可观察到除了强能量的基波外,还有2阶、3阶谐波存在,谐波能量较弱。该低频可控震源能有效激发低频信号,频率可低至1.5Hz(图 2(d))。由于低频段时震源车机械系统的耦合性能、振动底板与地面耦合作用等原因,1.5~3.0Hz激发信号相对较弱,3Hz以上信号可以稳定激发(图 2(d))。

图 2 可控震源激发信号 (a)扫描信号;(b)实际输出力信号;(c)扫描信号对应的Stockwell变换频谱分布;(d)实际输出力信号对应的Stockwell变换频谱分布
1.2.2 震源重复性分析

与爆炸震源不同,可控震源在激发过程中不会对场地造成破坏,可以在同一位置进行多次激发。震源车振动过程是精确可控的,每次振动时震源车与场地的耦合状态是一致的,因此在同一位置激发的信号具有可重复性。图 3(a)为震源车289次激发的输出力信号,不同次激发的输出力信号具有极高的一致性。将这些输出力信号的平均值作为参考信号,计算不同次激发输出力信号的互相关系数后发现(图 3(b)),绝大多数都接近1,只有4次小于0.99。这表明可控震源激发具有非常好的可重复性,可以通过对多次重复激发信号进行处理来提取其中的有效信号,同时,震源的可重复性也为动态监测地下介质状态变化提供了条件。

图 3 可控震源信号的可重复性 (a)不同次激发的输出力信号;(b)不同次激发的输出力信号的互相关系数
1.2.3 扫描信号子波压缩

可控震源激发的信号在传播过程中同样遵循地震波传播规律,但长时间的扫描使来自不同路径的信号叠加在一起,仪器接收到的振动记录更加复杂,使我们难以识别其中的震相信息。可控震源数据处理中的重要步骤就是压缩长时间扫描信号,提取其中的有效信号,也就是将可控震源数据转化成类似于爆炸震源形成的地震记录。这常通过将地震记录与震源扫描信号进行互相关来实现,也可以通过反褶积方式实现(Brittle et al,2001杨微等,2013)。反褶积处理后的频谱更宽,波形旁瓣收敛好,但是信噪比对数据处理结果的稳定性影响较大(杨微等,2013),故仅适用于具有一定信噪比的资料。可控震源信号经远距离传播后,有效信号已经被淹没在噪声中,数据信噪比很低。因此,本文采用互相关方式进行处理。

地震响应的褶积模型可以表示为

$ x\left(t \right) = r\left(t \right) * s\left(t \right) $ (1)

式中,x(t)为接收的地震记录;r(t)为脉冲响应;s(t)震源信号;*表示卷积算子。当将震源信号s(t)与地震记录x(t)互相关,则

$ x\left(t \right) \otimes s\left(t \right) = r\left(t \right) * s\left(t \right) \otimes s\left(t \right) $ (2)

将震源扫描信号s(t)的自相关函数称为Klauder子波k(t),它是一个有限频宽0相位子波,其频带宽度取决于扫描信号频谱。因此

$ x\left(t \right) \otimes s\left(t \right) = r\left(t \right) * k\left(t \right) $ (3)

将震源信号s(t)与地震记录x(t)互相关后得到的相关地震记录等价于震源以Klauder子波激发的地震记录。互相关运算不仅将长时间扫描信号压缩成Klauder子波,并且具有滤波作用,相关地震记录只保留了震源信号与地震记录中的共有频率成分。图 4(a)是扫描信号(图 2)对应的Klauder子波,为0相位子波,其频带范围为1.5~8.0Hz(图 4(b)),与扫描信号自身频谱一致。

图 4 庐江激发点1.5~8.0Hz扫描信号对应的klauder子波及其频谱 (a)Klauder子波;(b)Klauder子波对应的频谱
2 实际地震记录与数据处理 2.1 预处理

对记录的连续数据进行预处理以获得单次激发地震记录。首先进行常规预处理,包括去均、去势、去仪器响应、带通滤波(0.5~10.0Hz)、降低采样率等,考虑到可控震源扫描信号频带为1.5~8.0Hz,因而进行0.5~10.0Hz带通滤波与10ms降采样处理;再根据可控震源提供的各次激发的记录时间截取对应单次激发记录。由于最远偏移距约为75km,故我们设置截取记录长度为80s。由于环境噪声随时间的变化较大,我们只选取受环境噪声影响较小的激发记录进行处理。在各单次激发记录中,若某次激发数据方差大于平均值的3倍以上,则认为数据受噪声影响大,其不参与后续处理。

图 5为距震源车最近的6.3km处的流动台站接收到的地震记录。图 5(a)为经过预处理后的289次重复激发对应的记录(每间隔10道显示记录,只显示了前25 s的地震记录)。图 5(b)为与扫描信号互相关形成的相关记录道。对于单次激发而言,有效信号能量弱。图 5中有效信号被强噪声淹没,无法识别。背景噪声水平随着时间而变化,各次激发信号受背景噪声的影响并不相同,因而选取受背景噪声影响较小的记录进行处理。

图 5 重复激发接收记录 (a)接收台站距震源约6.3km原始记录;(b)与扫频信号互相关后得到的互相关记录
2.2 重复激发信号处理方式

单次激发的有效信号能量有限,利用多次重复激发的地震记录可以增强有效信号,处理方式包括多道记录叠加和有效信号分离等2类。由于实际信号处理结果缺乏准确解来评估处理结果的优劣,因此,本文对每种方法都采用含有噪声的合成多道地震记录进行相应处理,以比较不同处理方式的效果。合成多道地震记录中包含相同的15道地震记录,每道都含4个Klauder子波信号(图 6(a)),Klauder子波的最大振幅为1。在数据中加入标准差为0.9的高斯噪声(图 6(b)),有效信号已被强噪声淹没,无法识别。

图 6 含15道记录的合成地震数据 (a)不添加噪声的合成数据;(b)含有强随机噪声的合成数据
2.2.1 叠加方法

地震数据处理过程中,叠加是一种常用的增强有效信号、提高数据信噪比的处理方法(Kennett,2000),可分为线性叠加和非线性加权叠加。非线性加权叠加又可分为时间域加权叠加(Schimmel et al,1997Kennett,2000)和时频域加权叠加(Schimmel et al,2007)。相比于线性叠加,非线性加权叠加充分利用了有效信号的一致性,能够进一步压制噪声。本文对比了不同的叠加方法,包括常规线性叠加、时间域的Semblance加权叠加和相位加权叠加、时频域相位加权叠加等方法。

2.2.1.1 线性叠加

常用的直接线性叠加方法可以表示为

$ {g_{{\rm{ls}}}}\left(t \right) = \frac{1}{N}\sum\limits_{j = 1}^N {{s_j}\left(t \right)} $ (4)

其中,N为记录道数;sj(t)为第j道记录;gls(t)为线性叠加的结果。线性叠加使信号强度随叠加次数N成比例增长,而随机噪声按叠加次数的$ \sqrt N $倍增长,因而线性叠加可以使信噪比增加$ \sqrt N $倍(陆基孟等,2009)。对互相关后的多道记录进行线性叠加处理,图 7为可控震源实际数据线性叠加的结果,其中显示了部分信噪比较高的台站的Z分量记录。由图 7可见,在30km范围内有效信号较清晰,但整体上叠加结果的信噪比较低。

图 7 可控震源重复激发记录线性叠加结果
2.2.1.2 时间域加权叠加

在线性叠加过程中,对每1个时间采样点赋予了同样的权重。地震记录中的有效信号具有一致性,而噪声则互不相关。因此,可以测量记录的一致性,对波形一致性较高的时段赋予较大的权重,这进一步提高了叠加结果的信噪比。加权叠加方法已经在背景噪声提取经验格林函数、地震剖面数据处理中得到成功应用(Schimmel et al,1997Kennett,2000Schimmel et al,2011)。

时间域的加权叠加可以表示为

$ {g_{{\rm{ws}}}}\left(t \right) = C{\left(t \right)^\gamma }{g_{{\rm{ls}}}}\left(t \right) $ (5)

其中,gws(t)为加权叠加结果;gls(t)为线性叠加结果;C(t)为时间域的加权系数;指数γ控制着加权叠加中权重的影响程度,当γ=0时,加权叠加退化为常规线性叠加。

(1) 相位加权叠加

当利用瞬时相位的一致性作为计算权重时(Schimmel et al,1997),则称为相位加权叠加。相位加权系数表示为

$ C\left( t \right) = \frac{1}{N}\left| {\sum\limits_{j = 1}^N {{{\rm{e}}^{i{{\mathit{\Phi}}_j}(t)}}} } \right| $ (6)

其中,Φj(t)为第j道记录的瞬时相位,可以通过记录sj及其Hillbert变换H[sj(t)]获得。记录sj对应的解析信号(复数道信号) Sj

$ {S_j}\left(t \right) = {s_j}\left(t \right) + iH[{s_j}\left(t \right)] $ (7)

同时,解析信号又可由振幅项和瞬时相位来表示

$ {S_j}\left(t \right) = {A_j}\left(t \right){\rm{exp}}\left[ {i{{\mathit{\Phi}} _j}\left(t \right)} \right] $ (8)

其中,Φj(t)为对应的瞬时相位;Aj(t)为对应的振幅项。图 8(a)为可控震源实际数据相位加权叠加结果;图 8(b)为对应的相位加权叠加的权重系数。在数据处理中,加权指数因子γ统一取为2。相比于线性叠加(图 7),相位加权叠加(图 8(a))提高了叠加记录的信噪比,有效信号更加清晰,但同时也改变了叠加记录的波形。

图 8 可控震源重复激发记录相位加权叠加结果 (a)相位加权叠加后的地震记录;(b)对应的相位加权叠加系数

本文通过合成数据评估加权叠加方式对叠加记录的影响。图 9为合成数据相位加权叠加测试结果。图 9(a)为计算出的相位加权系数C(t),归一化后噪声时段的加权系数约为0.4。相位加权系数准确检测出4个有效信号,在有效信号时段内权重系数明显增加,因此应用该加权系数可以进一步增强有效信号。相位加权系数是由逐个点进行计算得到的,并不稳定,需要对加权曲线进行一定的光滑处理。由图 9(b)可见,相位加权叠加压制了叠加数据中的噪声,有效信号更加清晰。图 9(a)中权重系数的峰值恰好对应记录中0相位的klauder子波峰值,加权叠加后波形(图 9(b))的旁瓣发生改变,但波形的峰值位置并未改变。对于可控震源数据得到的互相关记录,震相的走时为对应的klauder子波的峰值位置,因而认为相位加权叠加对震相走时信息的影响较小。

图 9 合成数据相位加权叠加结果 (a)加权叠加系数;(b)相位加权叠加与线性叠加的结果对比

(2) Semblance加权叠加

波形相似性系数(Semblance)是一种常见的衡量多道信号一致性程度的参数(Neidell et al,1971),因而可用Semblance作为非线性叠加的加权系数C(t) (Kennett,2000),该加权叠加方式称为Semblance加权叠加。加权系数C(t)表示为

$ C\left( t \right) = \frac{{\sum\limits_{i = - n/2}^{n/2} {w(i)} {{\left[ {\sum\limits_{j = 1}^N {{s_j}\left( {t + i} \right)} } \right]}^2}}}{{N\sum\limits_{i = - n/2}^{n/2} {w(i)} \sum\limits_{j = 1}^N {s_{_j}^2\left( {t + i} \right)} }} $ (9)

其中,w为加权函数;n为窗函数的大小。加权函数可以取高斯窗函数$ \exp \left({ - \frac{{{{\left({t - {t_0}} \right)}^2}}}{{2{\tau ^2}}}} \right) $,宽度因子τ控制高斯窗的宽度。图 10(a)为可控震源实际数据Semblance加权叠加结果。计算Semblance的高斯窗口的选取依据地震数据的主频确定,一般应大于子波的宽度。此处宽度因子τ为0.3s。图 10(b)为作为加权系数的Semblance值。相比于线性叠加(图 7)和相位加权叠加(图 8(a))的结果,Semblance加权叠加(图 10(a))后有效信号更加清晰。对比相位加权的叠加权重曲线(图 8(b))发现,Semblance系数曲线(图 10(b))更加光滑稳定,因而Semblance加权叠加压制噪声的能力更强。

图 10 可控震源重复激发记录Semblance加权叠加结果 (a)Semblance加权叠加后的地震记录;(b)对应的Semblance加权系数

通过合成数据比较相位加权叠加和Semblance加权叠加2种方式。图 11为合成数据Semblance加权叠加结果,与相位加权系数曲线(图 9(a))相比,Semblance加权叠加是在选定的一定宽度的高斯窗内计算得到的,因而更加光滑稳定。在噪声时段(0~5s),Semblance加权系数(约为0.2) 比相位加权系数(约为0.4) 更小,因而Semblance加权叠加比相位加权叠加更能压制噪声。

图 11 合成数据Semblance加权叠加 (a)加权叠加系数;(b)Semblance加权叠加与线性叠加的结果对比
2.2.1.3 时频域相位加权叠加

实际地震信号往往具有一定的时频分布特征,而且噪声信号在时频域往往呈现随机分散分布的特征(王鹏等,2014)。利用Stockwell变换(Stockwell et al,1996)的可逆性,时间域相位加权叠加可以推广到时频域(Schimmel et al,2007Li et al,2014Thurber et al,2014)。

Stockwell变换可以表述为

$ S\left({\tau, f} \right) = \smallint _{_{ - \infty }}^{^{ + \infty }}s(t)w\left({\tau - t, f} \right){e^{ - i2\pi ft}}{\rm{d}}t $ (10)

其中,f为频率;w(τ-tf)为高斯窗口。

$ w\left({\tau - t, f} \right) = \frac{{|f|}}{{k\sqrt {2\pi } }}{\rm{exp}}\{ [ - {f^2}{\left({\tau - t} \right)^2}]/(2{k^2})\}, \;\;k > 0 $ (11)

在时频域计算相位加权系数Ctf(τf)

$ {C_{tf}}\left({\tau, f} \right) = \left| {\frac{1}{N}\sum\limits_{j = 1}^N {\frac{{{S_j}\left({\tau, f} \right){e^{i2\pi f\tau }}}}{{|{S_j}\left({\tau, f} \right)|}}} } \right| $ (12)

其中,Sj(τf)为第j道地震记录的S变换结果。最终,加权叠加结果由时频域内线性叠加结果乘以对应加权系数,再经Stockwell逆变换到时间域内得到

$ {S_{tf - {\rm{pws}}}}\left(t \right) = {\rm{inv}}\{ {C_{tf}}{\left({\tau, f} \right)^\gamma }{S_{{\rm{ls}}}}\left({\tau, f} \right)\} $ (13)

式中,指数γ为加权因子;Sls(τf)为时频域内线性叠加的结果。图 12为可控震源实际数据时频域相位加权叠加结果。相比于时间域相位加权叠加(图 8(a)),时频域相位加权叠加后有效信号更加清晰。

图 12 可控震源重复激发记录时频域相位加权叠加结果

通过合成数据测试可发现时频域相位加权叠加具有的优势。图 13(a)为合成数据在时频域获得的相位加权系数,加权系数准确地识别出了4个有效信号。加权系数不仅可以确定有效信号在时间上的分布,而且可确定有效信号的频带范围(1.5~8.0Hz)。图 13(b)为合成数据在时频域线性叠加的结果。由图 13(b)可见4个有效信号,噪声呈现随机分散的分布特征。图 13(c)为时频域加权叠加的结果,相比时间域的加权叠加,时频域相位加权叠加不仅能压制噪声时段的噪声,还能有效地压制不在有效信号频带范围内的噪声,因而叠加后有效信号更加清晰。

图 13 合成数据时频域相位加权叠加在时频域的显示 (a)时频域加权系数;(b)时频域线性叠加结果;(c)时频域相位加权叠加结果

图 14为不同加权叠加方式的效果对比。相比于线性叠加,加权叠加能有效提高叠加记录的信噪比。时频域相位加权叠加能有效压制非有效信号频段内的噪声,叠加结果受噪声的干扰更小。相比于时间域的加权叠加方式,时频域相位加权叠加结果与原始信号的相似度更高。

图 14 合成数据加权叠加对比
2.2.2 有效信号直接分离方式

对于重复激发的多道地震记录,依据其中有效信号与噪声的不同特征,利用特定的信号处理方法将地震记录分解成不同的成分,从而达到直接分离有效信号的目的。奇异值分解(Singular Value Decomposition,简称SVD)是一种常用的矩阵分解方式。

对于N道地震记录,每道地震记录有M个采样点,此N×M个数据可用矩阵X来表示

$ \mathit{\boldsymbol{X = }}\left[ {\begin{array}{*{20}{c}} {{x_{11}}} & {{x_{12}}} & \cdots & {{x_{1M}}}\\ {{x_{\, 21}}} & {{x_{22}}} & \cdots & {{x_{2M}}}\\ \vdots & \vdots & \ddots & \vdots \\ {{x_{N1}}} & {{x_{N2}}} & \cdots & {{x_{NM}}} \end{array}} \right] $ (14)

依据奇异值分解理论,矩阵X的奇异值分解可表示为

$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{UE}}{\mathit{\boldsymbol{V}}^T} = \sum\limits_{i = 1}^r {{\sigma _i}} {u_i}v_{_i}^{^T} $ (15)

其中,r为矩阵的秩,上标T表示转置算子;ui为矩阵XXT的第i个特征向量,U =[ u1u2,…,uN]为N×N正交矩阵;vi为矩阵XTX的第i个特征向量,V =[ v1v1,…,vM]为M×M正交矩阵;E =diag(σ1σ2,…,σr)为N×M的矩阵;σi为第i个奇异值。由SVD理论可知,矩阵X的总能量可以表示为

$ {\left\| \mathit{\boldsymbol{X}} \right\|^2} = \sum\limits_{i = 1}^r {\sigma _i^2} $ (16)

通常较大的奇异值主要对应数据中的有效信号的能量,而较小奇异值则对应噪声信号(Freire et al,1988Bekara et al,2007)。因此,可以选取k个较大奇异值来重构地震记录,以达到压制噪声的目的(Bekara et al,2007)。最后,对去噪后的多道地震记录进行叠加处理,获取最终的地震记录。我们将这种处理方式称为SVD提取信号方法。图 15为可控震源实际数据SVD叠加结果,此处我们认为,后一半较小的奇异值主要对应噪声信号,因此选取前一半较大的奇异值来重构地震记录。对比线性叠加结果(图 7),SVD叠加结果并没有显示优越性。

图 15 实际地震数据SVD叠加结果

我们使用了合成数据来验证SVD提取信号的效果。图 16为合成记录SVD叠加测试结果。首先,对合成地震记录进行SVD分解,最大的奇异值明显大于其他14个奇异值(图 16(a));第1个奇异值主要对应有效信号的能量,其余奇异值主要对应噪声能量。然后,用最大的奇异值重构地震记录发现(图 16(b)),有效信号清晰,强噪声得到压制。但重构后的数据中剩余噪声同样具有一致性,这些噪声在后续叠加中并不能被有效消除。含噪数据进行SVD分解后,噪声的存在对所有奇异值均有影响(Cichocki et al,2002),因此,选取部分奇异值重构信号的SVD去噪方法并不能完全去除噪声。

图 16 地震数据SVD去噪 (a)奇异值分布;(b)用最大的奇异值重构的地震数据

将重构后的15个地震记录叠加(图 17),形成最后的叠加记录。由图 17可见,SVD叠加结果与线性叠加结果非常相近。虽然SVD是一种有效的去噪技术,但将标准SVD作为增强重复激发记录有效信号的方法并不显示优越性。

图 17 地震数据奇异值分解(SVD)叠加
2.3 关于数据处理效果的讨论 2.3.1 不同处理顺序的影响

不同的处理顺序常对最终处理效果产生一定的影响。以距震源最近(相距6.3km)台站的Z分量数据比较不同处理顺序的影响。

互相关和线性叠加都属于线性算子,理论上处理顺序对结果没有影响。先互相关后线性叠加与先线性叠加后互相关得到的结果完全相同(图 18)。加权叠加属于非线性算子,不同的处理顺序会对最终结果产生显著影响。相比先时频域相位加权叠加再互相关,先互相关再时频域相位加权叠加能更好地压制噪声,提高了结果的信噪比(图 18)。其原因可能是,先加权叠加一定程度上会损害原始记录中扫频信号,故再进行互相关并不能取得很好的效果。

图 18 处理顺序对处理结果的影响
2.3.2 指数加权因子选取

指数加权因子γ的选取对叠加结果有显著的影响。在同样的噪声水平下,强有效信号相比弱有效信号具有更强的一致性。过小的γ值无法有效压制噪声。当γ=0时,等效于线性叠加,而过大的γ值会压制有效弱信号,不利于信号提取。以先互相关后时频域加权叠加为例,对比γ的不同取值发现,γ取值2~3是较合理的选择,能取得最佳处理效果(图 19)。其他加权叠加方式的指数加权因子γ的选取原则类似。

图 19 时频域加权叠加中加权因子γ的影响
2.3.3 加权叠加方式提高信噪比的能力

线性叠加可使信噪比随叠加次数N的增加而增大$ \sqrt N $倍(陆基孟等,2009),但加权叠加能在多大程度上提高数据的信噪比,目前尚未见相关报道。我们用合成数据比较时频域相位加权叠加及线性叠加时信噪比与叠加次数N之间的关系。此处,将信噪比简单定义为有效信号时段平均振幅与噪声时段平均振幅的比值。图 20为随叠加次数的增加信噪比的变化情况。由图 20可见,随着叠加次数的不断增加,线性叠加提高信噪比的程度不断减缓,而加权叠加并没有相应的减缓趋势。因此,随着叠加次数的不断增加,加权叠加应该具有更大的优势。

图 20 随叠加次数的增加信噪比的变化
2.3.4 对设计扫描信号与实际输出力信号进行数据处理的效果对比

震源车提供的设计扫描信号和实际震源输出力信号可作为可控震源参考信号,通常使用扫描信号与地震记录进行互相关处理。对各次激发的震源输出力信号进行10ms重采样后作为震源参考信号与地震数据进行互相关处理。图 21为经时频域相位加权叠加后的结果,除参考信号外,其他处理过程与图 12相同。对比图 12使用扫描信号作为参考信号的结果可以发现,利用输出力信号作为参考信号的处理效果与其相当。

图 21 震源输出力信号作为参考信号的处理结果
2.3.5 不同处理方式的效果比较

表 2对上述不同处理的结果进行了比较。

表 2 不同处理方式的比较
3 地震信号分析与有效传播距离的探讨 3.1 地震信号分析

根据黄耘等(2011)对郯庐断裂带鲁、苏、皖段及邻区地壳速度结构的研究,本文建立实验区的一维速度结构(图 22)。浅层波速根据表层岩石实验室测定结果(倪红玉等,2015),其他各层速度采用郯庐断裂带下扬子地块的P波速度,S波速度由vS=vP/1.732换算得到。同时,在地表处添加了500m厚的低速沉积层。

图 22 实验区一维速度结构

利用TauP震相走时计算程序(Crotwell et al,1990)计算对应的Pg、Sg震相走时发现,理论走时与有效信号相符(图 23),推断这些信号可能为Pg、Sg波。

图 23 与扫描信号互相关后时频域相位加权叠加得到的地震记录剖面
3.2 有效传播距离的探讨

低频可控震源车激发信号的有效传播距离决定了可控震源在区域尺度地下结构研究中的应用价值。地震波传播距离主要受场地激发条件、地震波在地下介质中的传播以及地震信号接收端背景噪声水平等3方面因素的影响。

为了估计可控震源信号的有效传播距离,我们将可控震源激发视为作用于地表的垂直单力源,利用Fk程序(Zhu et al,2002)合成扫描信号子波(图 2)对应的地震记录。介质模型如图 22,模型近地表设置500m厚的低速沉积层,沉积层具有较强的粘弹性,QS值为40,QP值为80,其余各层为弹性介质。单力源设置为震源车第2次激发的输出力(表 1),震源子波为可控震源扫描信号最大振幅随震中距的增加迅速衰减(图 24(b)),震中距大于30km时,振幅小于10-7 m/s;震中距为100km左右时,振幅约为10-8 m/s。本次实验流动台使用的数据采集器型号为REFTEK130,最小灵敏度约为7.94×10-10 m/s。因此,利用多次重复激发的记录实现区域尺度探测是可能的。

图 24 合成地震记录 (a)地表垂直单力源的合成地震记录;(b)最大振幅随震中距的变化

通常不同台站的背景噪声水平存在很大的差异,差异大小可按数量级衡量。根据地震台站观测环境技术要求的国家标准(杨建思等,2004),环境噪声水平可分为5级,其中最优的Ⅰ级环境噪声水平≤-150 dB。吴建平等(2012)利用华北流动地震台阵研究了华北地区的环境噪声特征,发现对于背景噪声水平较低的台站,1.0~10.0Hz频带内噪声水平小于-150 dB。我们由此估算,对应的噪声幅值可达10-7m/s量级(图 25)。此时噪声幅值与30km外的有效信号幅值接近同一个数量级。因此,震源车需要大量重复激发,才能在远距离处提取到有效信号。

图 25 随机噪声与噪声功率谱(据吴建平等(2012)) (a)0.5~10.0Hz带通滤波后的随机噪声(只显示前15s);(b)对应的功率谱密度函数

为了提高可控震源的有效传播距离,需要通过增加重复激发次数来加大激发能量。在激发前应当选定良好的激发场地条件,对台站噪声水平进行细致检查,选择背景噪声较弱的地震台站,最大限度地减少背景噪声的干扰。这些工作对提高有效信号的传播距离同样非常重要。

4 结论

低频可控震源能有效激发低频地震信号,其激发过程精确可控,对周边环境几乎没有影响,激发效率高,并且激发的信号具有极高的可重复性。

对于多次重复激发的地震记录,本文测试了2类增强有效信号的方式,即叠加方法和有效信号分离方法。结果表明,通过线性叠加处理可以增强有效信号。加权叠加利用有效信号的一致性能够在叠加过程中压制噪声,进一步增强有效信号。时间域的Semblance加权叠加比相位叠加更加稳定,压制噪声的效果更好。时频域相位加权叠加是时间域相位加权叠加的拓展,其充分利用了有效信号的时频分布特征和随机噪声在时频域随机分散分布的特征,在时频域内压制噪声干扰,叠加后的信号更加清晰。利用有效信号与干扰的不同特征可以将多道地震记录分解为不同的信号成分,从而直接分离出有效信号。SVD是一种常见的数据分解方式,对于本次实验,采用标准SVD分离获得的信号同线性叠加的结果类似。

通过对不同处理流程的分析发现,对于预处理后的数据,通过先互相关检测信号再进行时频域相位加权叠加能够取得最佳的处理效果。在本次实验中,经200多次重复激发后,可以在50km范围内提取到有效信号。为了增加有效信号的传播距离,在实验前应当做好场地激发条件和环境噪声水平的调查工作,选择背景噪声水平较低的台站。这些工作对于增加震源激发能量同等重要。

对于低频可控震源用于区域尺度结构的研究,目前进行的工作是初步的,能提取到的有效信号仅限于50km范围内。如何在更大范围内可靠地提取更加丰富的地震信号,如壳内反射折射信号、来自莫霍面的反射信号等,值得进一步研究。另外,研究更加先进的高效激发采集处理技术,提高激发效率,对于增强可控震源在区域尺度研究中的适用性也非常重要。

致谢: 东方物探公司提供了高精度可控震源车用于本次实验研究,中国科学技术大学地震与地球内部物理实验室的众多师生参与了地震台站布置和野外观测,特此感谢。感谢与王未来博士关于背景噪声水平的讨论。
参考文献
陈颙, 张先康, 丘学林, 等. 2007, 陆地人工激发地震波的一种新方法. 科学通报, 52(11): 1317–1321. DOI:10.3321/j.issn:0023-074X.2007.11.017
侯贺晟, 高锐, 贺日政, 等. 2012, 西南天山-塔里木盆地结合带浅深构造关系——深地震反射剖面的初步揭露. 地球物理学报, 55(12): 4116–4125. DOI:10.6038/j.issn.0001-5733.2012.12.024
黄耘, 李清河, 张元生, 等. 2011, 郯庐断裂带鲁苏皖段及邻区地壳速度结构. 地球物理学报, 54(10): 2549–2559. DOI:10.3969/j.issn.0001-5733.2011.10.012
林建民, 王宝善, 葛洪魁, 等. 2008, 大容量气枪震源特征及地震波传播的震相分析. 地球物理学报, 51(1): 206–212.
刘明辉, 周银兴, 彭朝勇, 等. 2015, 用能量累积法检测地震波雷达信号. 地球物理学报, 58(4): 1259–1268. DOI:10.6038/cjg20150414
刘振武, 撒利明, 董世泰, 等. 2013, 地震数据采集核心装备现状及发展方向. 石油地球物理勘探, 48(4): 663–675.
陆基孟, 王永刚. 2009, 地震勘探原理. 东营: 中国石油大学(华东)出版社.
倪红玉, 沈小七, 洪德全, 等. 2015, 2014年金寨ML3.9震群序列特征研究. 地震学报, 37(6): 925–936. DOI:10.11939/jass.2015.06.004
丘学林, 陈颙, 朱日祥, 等. 2007, 大容量气枪震源在海陆联测中的应用:南海北部试验结果分析. 科学通报, 52(4): 463–469.
佘德平, 管路平, 徐颖, 等. 2007, 应用低频信号提高高速玄武岩屏蔽层下的成像质量. 石油地球物理勘探, 42(5): 564–567.
陶知非, 苏振华, 赵永林, 等. 2010, 可控震源低频信号激发技术的最新进展. 物探装备, 20(1): 1–5.
陶知非, 赵永林, 马磊. 2011, 低频地震勘探与低频可控震源. 物探装备, 21(2): 71–76.
王海燕, 高锐, 卢占武, 等. 2010, 深地震反射剖面揭露大陆岩石圈精细结构. 地质学报, 84(6): 818–839.
王洪体, 庄灿涛, 薛兵, 等. 2009, 精密主动地震监测. 地球物理学报, 52(7): 1808–1815.
王鹏, 常旭, 王一博, 等. 2014, 基于时频稀疏性分析法的低信噪比微震事件识别与恢复. 地球物理学报, 57(8): 2656–2665. DOI:10.6038/cjg20140824
吴建平, 欧阳飚, 王未来, 等. 2012, 华北地区地震环境噪声特征研究. 地震学报, 34(6): 818–829.
杨建思, 黄媛, 徐智强, 等. 2004, 地震台站观测环境技术要求(第1部分)测震GB/T19531.2-2004. 北京: 中国标准出版社.
杨微, 葛洪魁, 王宝善, 等. 2010, 由精密控制人工震源观测到的绵竹5.6级地震前后波速变化. 地球物理学报, 53(5): 1149–1157.
杨微, 王宝善, 葛洪魁, 等. 2013, 精密控制机械震源特征及信号检测方法. 中国石油大学学报(自然科学版), 37(1): 50–55、69.
Baeten G, DeMaag J W, Plessix R E, et al. 2013, The use of low frequencies in a full waveform inversion and impedance inversion land seismic case study. Geophysical Prospecting, 61(4): 701–711. DOI:10.1111/gpr.2013.61.issue-4.
Bekara M, Van der Baan M. 2007, Local singular value decomposition for signal enhancement of seismic data. Geophysics, 72(2): V59–V65. DOI:10.1190/1.2435967.
Brittle K F, Lines L R, Dey A K. 2001, Vibroseis deconvolution:a comparison of cross correlation and frequency domain sweep deconvolution. Geophysical Prospecting, 49(6): 675–686. DOI:10.1046/j.1365-2478.2001.00291.x.
Cichocki A, Amari S. 2002, Adaptive blind signal and image processing:learning algorithms and applications. John Wiley & Sons.
Crotwell H P, Owens T J, Ritsema J. 1990, The TauP Toolkit:Flexible seismic travel-time and ray-path utilities. Seism Res Lett, 70(2): 154–160.
Freire S L M, Ulrych T J. 1988, Application of singular value decomposition to vertical seismic profiling. Geophysics, 53(6): 778–785. DOI:10.1190/1.1442513.
Kennett B L N. 2000, Research note:Stacking three-component seismograms. Geophys J Int, 141(1): 263–269. DOI:10.1046/j.1365-246X.2000.00048.x.
Li Q, Gao J. 2014, Application of seismic data stacking in time-frequency domain. IEEE, Geoscience and Remote Sensing Letters, 11(9): 1484–1488. DOI:10.1109/LGRS.2013.2296316.
Li Z, You Q, Ni S, et al. 2013, Waveform retrieval and phase identification for seismic data from the CASS experiment. Pure Appl Geophys, 170(5): 815–830. DOI:10.1007/s00024-012-0585-2.
Neidell N S, Taner M T. 1971, Semblance and other coherency measures for multichannel data. Geophysics, 36(3): 482–497. DOI:10.1190/1.1440186.
Schimmel M, Gallart J. 2007, Frequency-dependent phase coherence for noise suppression in seismic array data. J Geophys Res, 112(B4): B04303.
Schimmel M, Paulssen H. 1997, Noise reduction and detection of weak, coherent signals through phase-weighted stacks. Geophys J Int, 130(2): 497–505. DOI:10.1111/gji.1997.130.issue-2.
Schimmel M, Stutzmann E, Gallart J. 2011, Using instantaneous phase coherence for signal extraction from ambient noise data at a local to a global scale. Geophys J Int, 184(1): 494–506. DOI:10.1111/gji.2010.184.issue-1.
Stockwell R G, Mansinha L, Lowe R P. 1996, Localization of the complex spectrum:the S transform. IEEE Transactions on Signal Processing, 44(4): 998–1001. DOI:10.1109/78.492555.
Thurber C H, Zeng X, Thomas A M, et al. 2014, Phase-weighted stacking applied to low-frequency earthquakes. Bull Seism Soc Am, 104(5): 2567–2572. DOI:10.1785/0120140077.
Wei Z, Hall M A. 2011, Analyses of vibrator and geophone behavior on hard and soft ground. The Leading Edge, 30(2): 132–137. DOI:10.1190/1.3555320.
Zhu L, Rivera L A. 2002, A note on the dynamic and static displacements from a point source in multilayered media. Geophys J Int, 148(3): 619–627. DOI:10.1046/j.1365-246X.2002.01610.x.
Ziolkowski A, Hanssen P, Gatliff R, et al. 2003, Use of low frequencies for sub-basalt imaging. Geophysical Prospecting, 51(3): 169–182. DOI:10.1046/j.1365-2478.2003.00363.x.