地震勘探是透视地球内部、探查油气资源的重要手段.受复杂地表条件影响，地震勘探资料广泛存在强随机噪声，极大地制约了探测目标的实现.压制随机噪声，提高地震勘探资料信噪比是地震资料处理的关键环节.但随机噪声不具有确定性，仅存在统计特性，在低信噪比情况下压制随机噪声会损失一部分有效信号，严重影响进一步解释和成像.因此，压制噪声且保留有效信号是国内外学者一直致力于研究的热点问题(Herrmann et al.，2008; Lu Y H and Lu W K，2009; 邓小英等，2010；Deng et al.，2011；Liu et al.，2012; Naghizadeh and Sacchi，2012; 袁艳华等，2013; 李月等，2013).

时频峰值滤波(TFPF)是在时频分析基础上发展的随机噪声压制方法，在解决低信噪比、非平稳信号估计方面优势明显，目前TFPF已成功应用于地震勘探随机噪声消减(Barkat and Boashash，1999；Boashash and Mesbah，2004；Lin et al.，2007; Rosh and el Kahoo and Siahkoohi，2009).TFPF是一维随机噪声压制方法，无偏估计信号的前提条件是信号近似线性和噪声为高斯白噪声.在处理实际观测信号时，利用伪魏格纳威力分布(PWVD)的窗函数对信号局部线性化处理，从而实现对信号的无失真恢复(Swiercz，2006).但TFPF压制噪声和保持信号对窗长的要求是矛盾的，恢复有效信号需要小窗长以满足线性条件，而压制随机噪声需要长窗长(李月等，2009; Li et al.，2013).时变窗长时频峰值滤波基于信号和噪声的特征调节窗函数长度，有效解决了时频峰值滤波方法保持信号和压制噪声矛盾的问题(林红波等，2011; Lin et al.，2014).但这种方法对于有效信号和噪声的判别精度要求较高，对复杂的实际地震勘探记录滤波效果有待改善.因此，将一维TFPF推广到二维是TFPF方法研究的重要内容.二维TFPF方法结合地震同相轴的连续性，通过构建与同相轴时距曲线接近的径向或二次滤波轨线拓展时频峰值滤波(Wu et al.，2011; Tian and Li，2014).二维时频峰值滤波具有更好的保幅性能和噪声压制效果.但该方法假设地震勘探随机噪声是平稳的，其均值为常数且二阶矩不随时间变化，没有考虑实际地震勘探随机噪声空间分布特点.而地震勘探记录受环境因素和地表地质条件复杂变化影响，存在各道地震记录噪声不相关，且各道噪声功率不等的随机噪声，与阵列信号处理中的空间非平稳噪声类似(Moghaddamjoo，1991; 吴云韬等，2003).在空间非平稳地震勘探随机噪声情况下，二维时频峰值滤波方法会受噪声空间非平稳特性的影响，导致压制空间非平稳地震勘探随机噪声效果不理想，需要根据噪声空间特性改进滤波算法.

本文着重研究空间非平稳地震勘探随机噪声情况下径向时频峰值滤波随机噪声压制方法.首先介绍空间非平稳随机噪声，研究空间非平稳地震勘探随机噪声空间分布特点，并进一步分析空间非平稳地震勘探随机噪声在径向时频峰值滤波过程中的影响.在此基础上，提出基于绝对级差统计量(ROAD)的径向时频峰值滤波方法，抑制空间非平稳地震勘探随机噪声.

在地震信号处理研究中，通常将地震勘探随机噪声看作各阶统计量不随时间变化的平稳随机噪声处理.在空间分布上，各道地震记录的随机噪声功率通常比较接近，二阶统计量不随空间变化，对二维滤波方法影响较小.但在复杂地表条件和环境因素的影响下，有些地震记录中随机噪声的二阶统计特性随空间剧烈变化，不满足随空间变化宽平稳的条件.我们借用阵列信号处理中空间非平稳噪声来描述这种地震勘探随机噪声.

空间非平稳噪声是阵列信号处理研究中广泛关 注的问题(Moghaddamjoo，1991; 吴云韬等，2003). 空间非平稳噪声是指各天线阵元之间噪声不相关，且各个阵元噪声功率不等的随机噪声.对于这种空间非平稳噪声，假设噪声为零均值、时间和空间上均不相关的高斯过程，且噪声功率不相等，则噪声协方差矩阵为：

地震道相当于阵列信号处理中阵元接收到的信号，地震勘探记录中也存在各道地震勘探随机噪声不相关，但各道随机噪声功率随空间变化较大的情况，与阵列信号处理中空间非平稳随机噪声近似.以中国某地区实际地震勘探随机噪声记录为例，分析 各道地震勘探随机噪声空间分布特征和非平稳特性.

(1)各地震道噪声功率空间分布特征.统计空间非平稳随机噪声记录(图 1a)和平稳随机噪声记录(图 1c)各地震道噪声方差分析噪声功率的变化，图 1b为每一道空间非平稳噪声记录的方差，可以看出各地震道随机噪声方差随机波动，在40道和100道附近变化最为剧烈，该记录噪声功率随空间变化剧烈.与之相对应，图 1d中空间平稳地震勘探随机噪声的方差随空间变化较小，没有剧烈波动，在空间滤波过程中可近似为平稳随机噪声处理.

(2)空间非平稳特性度量.采用衡量非平稳特性的相关积分来分析地震勘探随机噪声沿空间采样的平稳特性(关惠玲等，2003).从地震记录中任选一点，任取三个方向构建采样路径，如图 1a中直线D1，D2和D3所示，分别统计图 1a和图 1c沿空间轨线D1，D2和D3采样随机噪声的相关积分，以度量空间采样噪声的非平稳性，统计结果记录于表 1.

从表 1中数据可以看出，图 1a记录中每条轨线上噪声的相关积分在0.77左右，图 1c记录中每条 轨线上噪声的相关积分为0.81.图 1a记录的相关积分值小，说明沿空间采样得到的噪声记录产生较大波动，在空间方向非平稳.这种空间非平稳随机噪声特性不满足常规地震勘探噪声压制方法对噪声性质的假设，易影响二维地震勘探噪声压制方法的去噪效果.本文将基于空间非平稳地震勘探随机噪声的特点与表征，研究其解决方案.

图 1

Fig. 1

图 1 地震勘探随机噪声记录 (a)空间非平稳地震勘探噪声记录；(b)记录(a)各道方差；(c)空间平稳地震勘探随机噪声记录；(d)记录(c)的各道方差. Fig. 1 Seismic noise record (a)Spatial nonstationary seismic noise；(b)The variance of seismic traces in(a)；(c)Spatial stationary seismic noise；(d)The variance of seismic traces in(c).

表 1 (Table 1)

表 1 中文标题 Table 1 英文标题 相关积分 D1 D2 D3 空间非平稳噪声 0.77 0.78 0.78 空间平稳噪声 0.81 0.81 0.81

表 1 径向轨线噪声相关积分 Table 1 The correlation integral of seismic noise on radial traces

径向时频峰值滤波沿径向滤波轨线执行时频峰值滤波.径向时频峰值滤波的滤波轨线可以看作是一系列斜率为k的平行直线，如图 2中平行虚线所示.对含有加性随机噪声的地震记录x(t，d)，沿平行径向滤波轨线重采样数据表示为

径向时频峰值滤波的滤波轨线和同相轴平行时，如图 2加粗轨线Trace1，径向轨线上的数据样本为同相轴同相位的数据，则径向轨线上数据幅度变化缓慢(图 2b).与沿时间方向的第一道地震信号(图 2c)相比，径向轨线上的信号更接近于线性，更好地满足时频峰值滤波获得无偏估计的线性条件.因此沿平行径向轨线时频峰值滤波有利于改善时频峰值滤波方法对信号的保持能力.然而，时频峰值滤波方法假设噪声是高斯白噪声，在空间非平稳地震勘探随机噪声情形下，滤波轨线上噪声的性质会受到噪声空间非平稳特性的影响，因此我们需要进一步分析实际地震勘探随机噪声在平行径向轨线上的特点及解决方案.

图 2

Fig. 2

图 2 平行径向滤波轨线示意图 (a)合成记录的平行径向滤波轨线；(b)径向道Trace1上的地震记录；(c)第一道地震记录. Fig. 2 Parallel radial traces (a)Parallel radial traces on synthetic seismic data;(b)Samples on the radial trace Trace1;(c)The 1st seismic trace in(a).

在利用径向时频峰值滤波方法滤波时，需针对这种空间非平稳噪声在径向轨线上噪声的特性改进滤波方法.以中国某地区40道实际地震勘探随机噪声为例分析地震勘探随机噪声空间分布表征.图 3a所示随机噪声记录各道功率波动较大，地震记录中16道至22道记录功率变化最为剧烈，噪声功率在空间上呈现非平稳分布.

当沿空间径向轨线(图 3a中灰色虚线所示)滤 波时，在径向轨线上的噪声性质会受到噪声空间非平稳性的影响.任取一道径向轨线Tr1，其上的信号如图 3b所示，采样点幅值在0.92 s附近突然增大，其波形类似脉冲噪声，这主要是由于地震勘探随机噪声在16道至22道噪声记录功率存在突变，在径向轨线Tr1相应区域(0.9至0.95 s之间)表现为近似脉冲的噪声.在采用径向时频峰值滤波压制噪声时，这类空间非平稳地震勘探随机噪声在径向轨线上的噪声性质不是平稳的高斯白噪声，不满足时频峰值滤波无偏估计的假设，会影响滤波效果.

图 3

Fig. 3

图 3 实际地震勘探随机噪声在径向滤波轨迹上的表征 (a)实际地震勘探随机噪声；(b)径向道上的随机噪声；(c)Tr1上随机噪声的绝对级差统计量(ROAD). Fig. 3 Field seismic r and om noise and its property on radial trace (a)Field seismic r and om noise;(b)R and om noise on radial traces;(c)ROAD of the seismic noise on Tr1.

对实际地震勘探随机噪声(图 4a)径向时频峰值滤波处理，滤波结果如图 4b所示.可见在噪声强度变化剧烈区域(方框所示)存在大量残留噪声.因为滤波轨线上脉冲噪声在滤波过程中影响附近点滤波结果幅值，0幅值变为非零值，甚至引起幅值符号的变化.径向轨线上脉冲噪声的相邻点对应于时空域相邻地震道上的点，继而引起时空域记录相邻道上对应点的幅值变化，产生残留噪声.剔除径向轨线上的脉冲噪声，然后再沿径向轨线时频峰值滤波处理，滤波结果(图 4c)在噪声强度变化剧烈区域几乎没有残留的随机噪声.从上述实验结果我们可以推断，残留噪声产生原因是随机噪声空间非平稳，可以通过抑制径向道上的脉冲噪声来改善径向时频峰值滤波压制空间非平稳随机噪声的效果.为此，我们结合脉冲噪声检测方法改进径向时频峰值滤波使之能够更好地压制这种空间非平稳地震勘探随机噪声.

图 4

Fig. 4

图 4 实际地震勘探随机噪声在径向滤波结果 (a)实际地震勘探随机噪声；(b)径向时频峰值滤波结果;(c)剔除脉冲噪声后RTFPF结果. Fig. 4 Field seismic r and om noise and filtered result by RTFPF (a)Field seismic r and om noise;(b)Filtered result by RTFPF;(c)Result from RTFPF on data removed pulse noise.

这里用绝对级差统计量(ROAD)检测脉冲噪声(Garnett et al.，2005)，脉冲噪声的ROAD统计量在强度上要高于非脉冲噪声的ROAD统计量，由此可以判断出脉冲噪声.

令Ω_M表示以0为中心的局部邻域窗的(2M+1)坐标，即Ω_M={u|－M≤u≤M}.对于径向轨线上的数据x′(t)，在以t时刻数据x′(t)为中心的局部邻域内，R_t定义为x′(t)与邻域内其他点的绝对差，即

图 3c为径向轨线Tr1上随机噪声的ROAD统计量值的变化.对照图 3b分析可见，对应于滤波轨线上强度较高噪声(0.9至0.95 s之间)，ROAD统计量值大于2，明显大于其他时刻噪声的ROAD值.较大的ROAD值说明径向轨线上这种能量突然变大的噪声近似脉冲噪声.这些近似脉冲噪声点对应地震勘探噪声记录中功率变化剧烈的地震道，这种近似脉冲的噪声是噪声能量空间不平稳地震勘探随机噪声在径向滤波轨线上的表现.

基于ROAD统计量，我们提出二次迭代径向时频峰值滤波算法.该方法首先计算径向滤波轨线上数据的ROAD统计量，将大于阈值的点判为脉冲噪声.然后对以脉冲噪声为中心的邻域内的数据局部 时频峰值滤波处理，滤波结果替代ROAD大于阈值T的脉冲噪声，我们将压制脉冲噪声后的记录表示为

通过脉冲噪声检测和迭代径向时频峰值滤波，ROAD径向时频峰值滤波方法能够有效压制空间非平稳地震勘探随机噪声.

为了验证所提出方法的有效性，我们将ROAD径向时频峰值滤波方法应用于合成地震记录.利用主频为30 Hz的Ricker子波构建40道地震信号模型，包含2个同相轴.Ricker子波的表达式为：

从图 5b可以看出，径向时频峰值滤波结果中大部分随机噪声得到很好地抑制，同相轴恢复良好，但对于随机噪声强度波动较大的16道至22道地震记录附近，虽然大部分随机噪声能够被压制，但在径向时频峰值滤波结果中残留许多随机干扰，这些随机干扰噪声将对后续处理产生不利的影响.相比之下，图 5c所示的ROAD径向时频峰值滤波结果背景更为干净，同相轴保持良好，不存在残留的随机干扰.可见我们提出的方法能够更好地抑制空间非平稳地震勘探随机噪声，保留有效信号.

为了进一步认识空间非平稳随机噪声对径向时频峰值滤波的影响，我们对比分析经ROAD处理前后所有径向轨线上的地震数据.图 6a为ROAD处理前的径向轨线上地震数据，在每条径向轨线上均出现类似脉冲的随机噪声.经ROAD处理后(图 6b)，径向轨线上的脉冲噪声被抑制，而有效信号没有受到影响.分析上述结果我们得到以下结论，随机噪声在空间非平稳会在径向滤波轨线上表现为脉冲噪声，若直接进行径向时频峰值滤波处理，压制噪声效果不理想，有大量残留的随机干扰.ROAD径向时频峰值滤波方法能够有效抑制脉冲噪声且不损伤有效信号，更好地压制空间非平稳地震勘探随机噪声.

图 5

Fig. 5

图 5 实际地震勘探随机噪声合成记录 (a)加入实际地震勘探随机噪声的合成数据；(b)径向时频峰值滤波处理结果；(c)本文方法处理结果. Fig. 5 Synthetic seismic data with field seismic r and om noise (a)Synthetic data;(b)Result by using radial TFPF;(c)Result by using our method.

图 6

Fig. 6

图 6 实际地震勘探随机噪声ROAD处理前后对比 (a)径向轨迹上的实际地震勘探随机噪声；(b)经ROAD处理的径向道上随机噪声. Fig. 6 Comparison between the seismic data on radial traces before and after ROAD processing (a)Field seismic r and om noise on radial traces;(b)R and om noise on radial traces after ROAD processing.

为了验证本文方法的实际应用效果，我们将本文方法应用于中国某地区实际共炮点记录，并与径向道时频峰值滤波进行对比.图 7给出了161道共炮点记录，每道时长5.5 s，采样时间为1 ms.该共炮点记录中含有较强的随机噪声，且部分区域随机 噪声能量较强，从空间上看各道噪声功率变化较明显，如框图所示.分别采用径向时频峰值滤波和本文方法对该共炮点记录进行处理，本文方法以炮点为中心向两侧取径向滤波轨线，结果分别显示于图 7b和图 7c.

图 7

Fig. 7

图 7 实际地震资料滤波结果对比 (a)共炮点地震勘探记录；(b)径向时频峰值滤波处理结果；(c)本文方法处理结果. Fig. 7 Comparison of filtering results of field seismic data (a)Field common-shot-point data;(b)Result from radial TFPF;(c)Result from our method.

比较两种方法的滤波结果可见，两种方法滤波结果中的随机噪声均得到有效抑制，同相轴能够从随机噪声中显现出来.但对原始记录中随机噪声突然增强部分(如矩形框所示)，径向时频峰值滤波结果在该区域依然残留部分随机干扰.而在ROAD径向时频峰值滤波结果中，该区域几乎没有残留的随机干扰，与径向时频峰值滤波结果相比，背景更干净，同相轴保持良好，没有畸变.对比上述滤波结果可见，本文的方法和已有方法相比，能够在保留有效信号的同时，较好地压制低信噪比地震记录中的随机噪声，对空间非平稳地震勘探随机噪声也能够有效压制，拓展了径向时频峰值滤波抑制噪声的能力.

本文针对径向时频峰值滤波不能有效去除空间非平稳地震勘探随机噪声的问题，提出了一种基于绝对级差统计量的径向时频峰值滤波算法.这类空间非平稳地震勘探随机噪声在径向轨线局部区域存在类似脉冲的噪声，与各道功率接近的地震勘探随机噪声相比具有较低的相关积分值，非平稳度更强.径向时频峰值滤波压制空间非平稳地震勘探随机噪声效果不理想，滤波结果在强度变化剧烈的地震道相邻区域存在大量残留噪声.理论分析和实验结果表明，结合绝对级差统计量和迭代径向时频峰值滤波方法能够检测并抑制空间非平稳地震勘探随机噪 声，实现对地震记录保幅压噪.实际地震记录处理结果也表明，ROAD径向时频峰值滤波改进了径向时频峰值滤波压制空间非平稳随机噪声能力和同相轴保持效果，将其应用于实际地震勘探处理是可行且有效的.