深地震反射剖面探测是地震学调查工具的首选，能够得到剖面下的精细结构，便于准确剖析复杂的构造地质问题(Clowes et al.，1987; 王海燕等，2010). 因此，深地震反射剖面一般穿越构造复杂的结合带，而在这些构造复杂地区(如地形起伏较大，断裂带发育)中，激发和接受条件都很差，导致了深地震反射数据一般反射信号弱和信噪比低，干扰类型多而复杂且差别大(Bois et al.，1988; 吴宣志等，1995; 杨文采等，1999; 崔兴宝，2003; Ross et al.，2004; 杨卓欣等，2006; 刘金凯等，2010; Karplus et al.，2011; Kumar et al.，2011). 此外，为了兼顾获得剖面下的浅、中、深层反射信息，爆炸震源激发时采取不同药量和井深组合以及变观测系统接收，从而导致激发、接收条件变化较大，同一深地震反射剖面的地震信号在频带、振幅、能量、相位等方面都存在着较大的差异(朱小三等，2013). 深地震反射数据的采集成本很高，而且重复采集几乎不可能，如果这些炮集位于关键目标构造带内，就会对整个区域的构造解释产生很大的影响. 因此，对于如何有效地从深地震反射数据中剔除干扰波，提高信噪比，前人曾做过大量尝试(朱小三等，2013)，虽然取得了一定的效果，但仍有不足.

S变换结合了短时傅里叶变换与小波变换的优点，克服了短时傅里叶变换时频窗口在时频平面中不可变的问题，其时窗宽度随频率呈反向变化，又无须满足小波变换中小波在时间域均值为零的容许性条件(Stockwell et al.，1996; 赵淑红和朱光明，2007; 陈学华等，2008). 因此，S变换是非平稳信号时频分析的良好工具，其高质量的时频表示为具体应用中设计有效的时频滤波器提供了保证. S变换与傅里叶变换的直接联系使得其算法实现简洁高效(Stockwell et al.，1996). 本文利用S变换的这些优点重新处理了西南天山深地震反射剖面中的废炮数据(高锐等，2008; He et al.，2008; Gao et al.，2013)，通过综合对比分析，效果明显，这些废炮数据得以重新使用.

S变换是由美国地球物理学家Stockwell于1996年提出的一种类似于短时傅里叶变换(Gabor，1946)的时频谱表示方法，函数h(t)的S变换定义为(Stockwell et al.，1996)：

图 1

Fig. 1

图 1 不同频率f所对应的时窗函数ω(t,f,τ)(其中τ=0)Fig. 1 Gaussian time window function ω(t,f,τ) for various values of f (τ=0)

高斯窗函数的缩放特性类似于连续小波(Mallat，1999)，而式(1)中用到的小波并不满足小波变换的小波在时间域均值为零的容许性条件(Stockwell et al.，1996). 因此，S变换并不是严格意义上的小波变换. 而且，S变换小波的振荡部分e^－i2πft是傅里叶复正弦曲线(complex Fourier sinusoid)，当τ变化时，它不随高斯窗形状的改变而改变(Pinnegar and Mansinha，2003). 因此，S变换小波的实部和虚部的形状随着高斯窗在时间上的移动(即τ的改变)而改变. 真正的小波变换的高斯小波在时间轴移动时(即τ改变时)，形状并不改变(Parolai，2009). 小波变换方法有较高的时间分辨率，低频区分辨率较高，高频区分辨率较低，但其相位局部化，造成各频率相位基准不一，从而造成解释的困难，而S变换对其相位进行了校正(Stockwell et al.，1996； Chakraborty and Okaya，1995; 陈雨红等，2006).

S变换也可以由h(t)的傅里叶频谱H(f)表示(Stockwell，1999)：

S变换表示的是局部的频谱特征，将其在时间方向积分可以容易地得到傅里叶频谱(Stockwell et al.，1996)：

阈值滤波主要分为硬阈值法和软阈值法，最早是由斯坦福大学的Donoho D L和Johnstone I M在估计小波系数时提出的(Donoho and Johnstone，1994a，1994b; Donoho et al.，1995; Donoho，1995).硬阈值法是指: 信号的系数大于一个确定的阈值λ时，保持不变，系数小于阈值λ时置零. 硬阈值法运算量小，操作简单，但函数在阈值λ处不连续，估计的方差较大，对数据的微小变化敏感(Donoho and Johnstone，1994a，1994b; 侯建华，2003).软阈值法是指: 信号的系数大于一个确定的阈值λ时，系数逐渐变小使λ趋近于0，系数小于阈值λ时置零. 软阈值法估计的偏差较大，对代表有用信息的数据进行了压缩，易出现过平滑现象(Donoho et al.，1994a，1994b; 侯建华，2003). 但对大多数的信号而言，软阈值滤波的效果要好于硬阈值滤波(Yoon and Vaidyanathan，2004). 因此本文采用软阈值滤波方法对数据进行处理，

在本文的信号去噪处理过程中，x代表S变换系数的幅值，λ取初至波到达之前的噪声每一频点的平均值，即

本文所用到的数据来自盆山结合带深部结构与油气远景研究项目(He et al.，2008; 侯贺晟等，2012; Gao et al.，2013). 地表条件复杂，经过平原带、山前带、高山区等各种类型的地形和构造结合带，原始单炮资料属典型的山地低信噪比资料(高锐等，2008; He et al.，2008; Gao et al.，2013)，这类资料在深地震反射剖面(Alsdorf et al.，1998; Langin et al.，2003; Karplus et al.，2011; Kumar et al.，2011)中普遍存在. 在数据采集过程中，如图 2所示的炮集遇到不明震源干扰. 当时无法很好地排除其影响，只能弃之不用. 考虑到数据缺失对整个区域构造解释的影响，本文选取其中一个废炮记录进行处理，以期最大可能的获得有用信号. 该单炮的激发药量为40 kg，炮点埋深为24 m，炮点桩号为1541.5，检波器排列的起止桩号为1042-2041，道间距为50 m，共1000道，采样率为4 ms. 文中选取桩号1042-1441共400道30 s的数据进行尝试性处理(高锐等，2008; He et al.，2008; Gao et al.，2013).

图 2

Fig. 2

图 2 原始单炮记录(方框A1、B1、C1为图8中对比分析区域)Fig. 2 The raw shot set ( boxes A1, B1, C1 are used to compare in Fig.8)

图 2为原始单炮记录，显示该炮集中的干扰波能量较强，信噪比较低，5 s以上反射波组较弱，连续性差，很难识别有效信号. 图 3为所有道集的频谱叠加图，该单炮记录的能量主要集中在5~50 Hz范围内，其中5~25 Hz的频带范围(图中方框a区域)包含了该单炮深部和浅部有效信号的主要能量; 25~50 Hz的频带范围(图中方框b区域)包含少量的浅部有效信号，但以强干扰噪声为主. 这与处理报告分析结果相似(高锐等，2008). 由于深地震反射主要目的是为了探测深部结构信息，目标信号以低频为主，且25~50 Hz频带范围内的强干扰噪声对突出显示有用信号不利，因此本文只针对5~25 Hz范围内的信号进行精细化处理.

图 3

Fig. 3

图 3 原始单炮记录频谱叠加分析图(方框a反映深浅部主要有效信息，方框b反映强干扰噪声和部分浅部有效信息)Fig. 3 The stacked Fourier amplitude spectra of raw shot set (The box of a reveals seismic reflection signals, the box of b reveals strong disturbing waves and a little shallow layers reflection signals)

利用巴特沃斯时不变带通滤波器，保留有效频带5~25 Hz范围之内的信号. 地震波在地下介质传播，由于大地吸收衰减使得地震波的能量随着传播距离的增加而逐渐的衰减，尤其对于大偏移距(超过数公里，甚至十多公里)的深地震反射而言，这种现象尤为明显. 为了突出显示深部的有效能量信号，对经过带通滤波之后的单炮记录进行瞬时自动增益. 在实际处理中选取的瞬时自动增益时间门值通常在256~1024 ms范围内(渥·伊尔马滋，1994)，本文选取512 ms的瞬时自动增益时间门对带通滤波之后的单炮记录进行处理，处理结果如图 5所示. 带通滤波的效果明显，有效滤除了大部分的干扰噪声，信噪比得到显著的提高，深部和浅部的有效能量信号均得以突显，浅部10 s以上的壳内反射信息丰富，深部20~25 s范围内Moho界面反射清晰可见. 但是在5~25 Hz的通带范围内依然存在混频干扰噪声，这部分噪声严重影响了相邻的有效波组之间的区分，同时掩盖了部分深部弱反射的有效信号. 因此，为了更加清晰地区分有效波组和突出显示弱反射有效信号，本文利用基于S变换的软阈值滤波方法对瞬时自动增益后的带通数据(图 4)进行处理，逐频点压制有效频带范围(5~25 Hz)内的混频噪声.

图 4

Fig. 4

图 4 时不变带通滤波(5~25 Hz)后瞬时自动增益的单炮记录(方框A2、B2、C2为在图8中效果对比区域, 与图2中A1、B1、C1分别代表相同时窗)Fig. 4 The shot set after time-invariant bandpass filtering (5~25 Hz) and AGC (boxes A2, B2, C2 are used to compare in Fig.8, which are the same time windows with A1，B1，C1 shown in Fig.2，separately)

图 5

Fig. 5

图 5 阈值λ随道号和频率的变换Fig. 5 The variation of threshold λ with trace and frequency

对带通滤波后瞬时自动增益的单道数据进行S变换，根据公式(9)求取初至波到达之前实际噪声各频点的阈值λ. 图 5直观地描述了阈值λ随道号和频率的变化，如1041-1091道整体噪声能量相对较弱，阈值λ整体相对较小；1241-1291道中15~25 Hz的区域和1310-1390道中8~18 Hz的区域噪声能量较强，阈值λ较大. 依据图 5中的λ值，逐频点对各道|S(τ，f)|的系数进行软阈值滤波处理. 处理结果如图 6 所示，与时不变带通滤波后瞬时自动增益的单炮记录(图 4)相比，图 6显示的单炮记录更加干净，初至明确，同相轴连续，Moho反射界面清晰，且深部的反射波组更加丰富，有利于进一步提高地层分辨率.

图 6

Fig. 6

图 6 基于S变换的软阈值滤波后的单炮记录 (方框A3、B3、C3为图8中的对比区域, 与图2中A1、B1、C1和图4中A2、B2、C2分别代表相同时窗)Fig. 6 The shot set after S-transform-based soft threshold filtering (boxes A3, B3, C3 are used to compare in Fig.8, which are the same time windows with A1，B1，C1 shown in Fig.2 and A2，B2，C2 shown in Fig.4，separately)

为了更好地说明滤波效果，图 7给出了第1292道数据在不同处理阶段的地震波形态及其相对应的S变换时频分析图. 初至波大致在2.2 s左右到达(如图 7中红线所示). 原始数据(图 7a)的S变换时频分析图(图 7b)可以看出地震波的有效信号主要集中在5~25 Hz的频带范围内(图 7b黄色方框所示)，再次印证了图 3频谱分析图所确定的时不变带通滤波的通带范围(5~25 Hz)的合理性. 同时也明显地看出深部地震信号能量较弱，不利于地震信号的分析处理. 因此在带通滤波之后，对地震信号进行瞬时自动增益处理，以便突出显示深部的有效能量，处理结果如图 7c所示. 为了压制5~25 Hz范围内的混频噪声干扰，本文运用基于S变换的软阈值滤波方法对图 7c数据进行处理，之后的处理结果(图 7e)清晰地显示该方法能够明显压制有效频带范围(5~25 Hz)内的干扰波，初至之前的噪声被显著地压制，而在初至之后，干扰波被压制的同时，有效信号的能量基本被保留，显著地提高了信噪比.

图 7

Fig. 7

图 7 第1292道数据在不同数据处理阶段的波形及其对应的S变换时频分析图 (a) 原始单道记录; (c) 时不变带通滤波(5~25 Hz)后瞬时自动增益的单道记录; (e)基于S变换的软阈值滤波后的单道记录; (b), (d), (f)分别对应图(a), (c), (e)的S变换时频分析图.Fig. 7 The seismic record of trace 1292 and its S transform amplitude spectra on different data processing stages (a) The raw record; (c) The record after time-invariant bandpass filtering and AGC; (e) The record after soft threshold filtering; (b), (d) and (f) are the S transform amplitude spectra of the Fig.(a), (c) and (e).

为了进一步突出显示整体滤波效果，选取了三个子区域(A、B和C)在不同处理阶段的结果进行对比，如图 8所示. A1、B1和C1区域为原始信号(图 2); A2、B2和C2区域为经过时不变带通滤波(有效频带范围5~25 Hz)并进行瞬时自动增益后的信号(图 4); A3、B3和C3区域是在A2、B2和C2处理基础之上执行S变换软阈值滤波的结果(图 6).

图 8

Fig. 8

图 8 处理效果局部对比图(A1、B1、C1为原始记录, 见图2; A2、B2、C2为时不变带通滤波(5~25 Hz)后瞬时自动增益的记录, 见图4; A3、B3、C3为基于S变换软阈值滤波之后的记录, 见图6)Fig. 8 Comparison of local shot set (A1, B1, C1 is the raw record in Fig.2; A2, B2, C2 is the record after time-invariant bandpass filtering (5~25 Hz) and AGC in Fig.4; A3, B3, C3 is the record after soft threshold filtering based on S transform in Fig.6)

原始单炮局部记录(图 8的A1)中可以看到有效波组，但是波组毛刺较多，经过带通滤波后瞬时自动增益的单炮局部记录(图 8的A2)中反射波组更加圆滑，反射波同相轴更加的清晰，再经过基于S变换软阈值滤波后的单炮局部记录(图 8的A3)部分反射波同相轴更加连续，同时使得部分弱反射波组更加突出. 原始单炮局部记录(图 8的B1)中1215-1240道有较强的高频噪声，有效波混杂其中，经过带通滤波后瞬时自动增益的单炮局部记录(图 8的B2)高频噪声被有效去除，再通过基于S变换软阈值滤波后(图 8的B3)不仅反射波同相轴更加连续，而且反射波组相当丰富，相邻反射波组也被有效地分开，有效弱反射信号能量进一步增强，有利于薄地层的识别和同相轴的追踪. 原始单炮局部记录(图 8的C1)中有较强的高频噪声，基本看不出有效反射波信息，经过带通滤波后瞬时自动增益的单炮局部记录(图 8的C2)高频噪声被有效去除，可以识别出几个连续的有效反射波组，再经过基于S变换软阈值滤波后(图 8的C3)背景噪声被进一步压制，较弱有效反射波得到进一步的突显，信噪比显著提高. 从图 8的对比可以看出，基于S变换的软阈值滤波在深地震反射单炮数据滤波处理中取得了显著的效果.

虽然经过本文精细化处理之后，原来的废炮数据现得以重新使用，但如图 6所示的单炮记录中不明震源所形成的规则强干扰仍未完全消除，为减少强干扰数据对后续共深度点叠加数据的影响，需要剔除强干扰信号. 如果坏道或强干扰波出现在关键构造目标区内，剔除后所造成的数据缺失将不利于揭示该区域复杂的地质构造，重建剔除区域的地震波场(吴蔚等，2014)十分重要. 因此，关于深地震反射数据的处理仍需进一步研究.

本文运用了基于S变换的软阈值滤波方法对西南天山的废炮数据进行实验性的去噪处理，取得了较好的效果. 与传统的滤波方法(朱小三等，2013)相比，基于S变换的软阈值滤波方法去噪效果明显，压制了有效频带范围内的混频干扰噪声，突出了弱反射信号，提高了信噪比，更便于有效反射波组的连续追踪和薄地层的识别. 另外，S变换可以在频率域实现，有较高的计算效率. 因此，基于S变换的软阈值滤波方法不仅可以用于地震反射数据去噪处理，也适用于其他地震学数据去噪处理.

致谢 感谢中国石化集团华东石油局第六物探大队采集了大量的野外资料，感谢侯贺晟副研究员从数据库中提取了本文处理所需的原始数据，感谢北京派特森科技发展有限公司在数据处理中提供的宝贵经验，感谢中国地质科学院地质研究所岩石圈中心的所有成员对作者写作过程中的指导与帮助. 文中的数据处理和成图是通过Matlab、Mathematica和CorelDRAW软件来完成的.