2. 吉林大学信息工程系, 长春 130012;
3. 吉林大学地球物理系, 长春 130026
2. Department of Information and Engineering, Jilin University, Changchun 130012, China;
3. Department of Geophysics, Jilin University, Changchun 130026, China
地震勘探是油气勘探的重要手段之一.由于数据采集环境复杂,地震勘探记录中常常混杂有大量的随机噪声.随机噪声的存在成为了获取高精度地震记录的主要障碍之一.要有效地压制随机噪声,首先要对地震勘探随机噪声的性质进行科学的研究.这里所说的随机噪声是指地震记录中的非相干噪声,它主要是由风、环境噪声、记录仪器噪声和检波器与地面耦合不好引起的(Yilmaz and Doherty, 2001).随机噪声没有固定的视频率和视速度,通过统计方法对随机噪声的特性进行研究是相对可行的途径.目前对于地震噪声的研究是地球物理学界的热点之一,研究主要集中在微震和环境地球物理学领域(Persson, 2003; 鲁来玉等, 2009; Groos and Ritter, 2009; Akhouari et al., 2011; 赵盼盼等, 2012; 潘佳铁等, 2014).在这些领域内,地球物理学家关心的频段远低于地震勘探中有效信号分布的频段(McNamara and Buland, 2004),因此,上述领域中的研究成果不能直接应用于地震勘探领域.与此同时,目前地震勘探随机噪声的特性还没有得到较系统的研究,只有少量论文发表.基于以上原因,本文对陆地地震勘探随机噪声的统计特性进行了系统的分析.
本文中主要研究的随机噪声特性包括:平稳性、高斯性和线性.这里所说的平稳性是指广义平稳性(wide-sense stationary),即随机过程的一阶统计量和二阶统计量是时不变的(Koopmans, 1995; Chatfield, 2003).高斯性是指地震勘探随机噪声幅度是服从高斯分布的.线性是指随机噪声是由线性系统产生的(Theiler et al., 1992).从物理学角度看,如果随机噪声是线性的,意味着在噪声源同检波器接收到的数据之间是线性关系,也就是说可以通过线性模型来对随机噪声的传播过程进行建模(Bath, 1974).在地震勘探领域,随机噪声普遍认为是平稳、高斯随机过程.陆基孟(1991)认为地震勘探中的随机干扰是具有各态历经性质的平稳随机过程.李庆忠(1993)认为地震勘探记录中的随机噪声具有(空间与时间)平稳随机性.同时,在理想条件下,即存在大量独立噪声源时,根据中心极限定理,地震勘探随机噪声应该是服从高斯分布的(Bendat and Piersol, 1994).基于此认识,地球物理学家提出了一系列随机噪声压制方法,例如维纳滤波(Mendel, 1977)和F-X预测滤波(Canales, 1984).然而在实际应用中,上述方法在某些情况下处理结果并不理想,有效信号容易发生畸变,去噪结果信噪比提高有限(Gulunay, 2008; Wu et al., 2011).因此,针对噪声特性的研究对于改进和设计噪声压制方法具有重要的理论和实际意义.
在本文中,噪声数据是按照实际地震勘探要求,在不同典型地理环境(沙漠、草原和山地地区)中测得的.在此基础上,应用一系列基于双谱(bispectrum)的假设检验方法(Hinich, 1982; Hinich and Messer, 1995)对噪声的性质进行研究.上述假设检验方法主要根据检测序列双谱值的分布特点,对检测序列的性质进行判断.以上方法已成功应用于医学信号检测(Shirazi and Moussavi, 2011)、海洋噪声特性分析(Brockett et al., 1987)和震相分析(Persson, 2003)等领域.
2 实验数据及方法 2.1 噪声数据文中应用到的数据是按照实际地震勘探要求采集得到的.使用到的检波器型号为JF-20DX-10,它是一种地震勘探中常用的垂直分量检波器.分别在沙漠、草原和山地地区记录地震勘探随机噪声.同时,记录环境比较偏僻,人文噪声较弱.沙漠地区每幅噪声记录由512道记录组成,道间距为50 m;草原地区每幅噪声记录也包含512道记录,道间距为30 m;山地地区每幅噪声记录由256道记录构成,道间距为40 m.三个地区记录的采样频率均为1 ms.
Cooper和Cook (1984)指出在无人工震源激发的情况下,地震记录中的相干噪声源主要包括:机动车噪声、机械噪声和输电线路带来的工频干扰(50 Hz或60 Hz的固定频率干扰)等.由于记录噪声区域比较偏僻,人文噪声相对较弱,因此可以近似认为采集到的噪声记录中只包含随机噪声.同时,基于相干噪声的相关特性,计算了每幅地震记录中噪声的互相关函数,只有没有明显相干噪声存在的记录才能被用以统计随机噪声的特性.图 1为1幅3 s长草原地区实际噪声记录,观察发现记录中并没有明显的相干噪声存在,说明假定噪声记录中只含有随机噪声的假设是合理的.
在实际情况下,测得的噪声序列都是有限长的,只能基于有限长序列对产生噪声的随机过程的特性进行估计.一般情况下,双谱估计结果的估计方差较大(Brillinger, 1975), 减小估计方差常用的方法是将长噪声序列分割成等长度的子序列,通过对子序列双谱进行平均或对得到结果在双频域进行平滑从而减小估计的方差.在本文中,将长度为N的序列x(n)分为P段长为L的子序列xp(l)(p=1, 2, …, P, l=1, 2, …, L).子序列的功率谱估计表示为P(k)=
(1) |
B(k1, k2)并不是序列x(n)双谱的一致估计.为了获得待测序列双谱的一致估计
(2) |
式中
本文采用的平稳性检验方法(Hinich and Messer, 1995)通过检验待测序列的双谱值在双频域外三角形区域(简称OT域)内的分布情况,从而对待测序列的平稳性进行判断.OT域是由三条直线围成的范围:2k1+k2=L,k1+k2=
(3) |
式中参数取值范围满足:2k1+k2 < L,k1+k2>
(4) |
可以证明当L <
高斯性和线性检验方法(Hinich, 1982)主要通过检验待测序列的双谱值在主体域(PD域)内的分布情况来对待测序列的高斯性和线性进行判断.PD域是由三条直线围成的范围:2k1+k2=L,k1=k2和k2=0.定义变量:
(5) |
如果待测序列满足高斯分布,Γ(k1, k2)为一单位方差复高斯变量.在此基础上构建检验统计量D:
(6) |
对高斯随机过程而言,假设PD域存在Q个不相关的双谱值,D应该近似服从自由度为2Q的中心χ2分布χ2Q2(0).进而由χ2分布特性可知,当自由度2Q足够大(2Q>50)时,D可以近似认为服从均值为2Q、方差为4Q的正态分布,因此可以构建高斯性检验统计量:
(7) |
DG满足标准正态分布,在显著水平设定为5%的情况下,选定高斯性判决门限THRG=1.96.当检验统计量|DG|≤THRG时,待测序列被判定为高斯序列;当检验统计量|DG| >THRG时,待测序列被判定为非高斯序列.
高斯序列一般情况下也是线性的(Persson, 2003),这里我们只需要对非高斯序列的高斯性进行判断.如果待测序列是线性非高斯序列,检验统计量D应该近似满足自由度为2Q非中心χ2分布χ2Q2(
(8) |
在此基础上,求取χ2Q2(
(9) |
在显著水平为5%的情况下,采用单侧置信区间,线性假设判定门限THRL=1.64.当检验统计量|DL |≤THRL时,待测序列被判定为线性时间序列;当检验统计量|DL |≤THRL时,待测序列被判定为非线性时间序列.
3 实验结果对于噪声特性的研究主要分三步进行:首先,对噪声的平稳性进行判断.因为高斯性和线性检验方法需要待测序列为平稳或广义平稳时间序列,这里需要通过平稳性检验剔除非平稳噪声;其次,应用高斯性检验方法,对随机噪声的高斯性进行判断;最后,对随机噪声的线性特性进行研究.服从高斯分布的随机噪声序列一般是线性的,故只需要对非高斯序列的随机噪声的线性特性进行判断即可.在实际情况下,地震勘探记录的时长都很短,本文只关注随机噪声在短时段内的统计特性.
3.1 随机噪声的平稳性应用平稳性检验方法对随机噪声的平稳性进行研究.图 2a和2b分别给出采自沙漠和山地地区的3 s长随机噪声的时域波形图,图 2c和2d分别为两道随机噪声记录对应的时频谱图.如图所示,采自沙漠的随机噪声记录时域波动更为平稳、时频域能量分布更为均匀,直观判断其平稳性相对较好.平稳性检验结果也证明了这一点,沙漠地区噪声记录的检验统计量DST=0.76小于判决门限值1.64,被判定为平稳时间序列,而山地地区噪声记录的检验统计量DST=2.77大于判决门限值,被判定为非平稳时间序列.图 3所示为不同地区不同时长随机噪声序列中平稳序列所占的百分比.由图可知,随机噪声的平稳性随记录时长的增加而变差,同时采集环境对随机噪声的平稳性也有着重要影响,对于相同时长随机噪声而言,沙漠地区随机噪声平稳性最好,山地地区随机噪声平稳性最差.非平稳随机噪声记录的出现,说明地震勘探随机噪声并不是传统意义上认为的平稳随机过程,而应归为局部平稳随机过程.也就是说,非平稳随机噪声序列在某些时段内的子序列可能是平稳的.同时,对于小于3 s随机噪声而言,超过65%的随机噪声被判定为平稳时间序列,因此时长小于3 s的随机噪声可以近似认为是平稳随机过程.
在此基础上,在频域内比较平稳和非平稳随机噪声记录在能量分布上的差异.对随机噪声的功率谱进行归一化处理后,求取了其在频域的分布情况,结果如图 4所示.图 4a为草原地区平稳随机噪声能量的分布情况,图 4b为草原地区非平稳随机噪声对应的结果,从中发现平稳随机噪声和非平稳随机噪声在低频段( < 50 Hz)具有相近的能量分布,平稳随机噪声较非平稳随机噪声而言其能量在中高频段([50~250 Hz])分布较为集中,而非平稳随机噪声在中高频段具有更多的能量且分布更为无序.图 4c到4f为沙漠和山地地区随机噪声能量在频域的分布情况,可以得到与草原地区相近的结论,平稳随机噪声能量分布相对集中,而非平稳随机噪声的能量分布在中高频段相对无序.以上结果表明,地震勘探随机噪声不是平稳随机过程,其平稳性受到噪声时长和采集环境的影响,同时平稳随机噪声较非平稳随机噪声,其能量在频域的分布情况更为规则.
采用前述高斯性检验法.图 5a和5b分别为两道采自沙漠地区和草原地区3 s长随机噪声记录的时域波形,沙漠地区随机噪声序列的高斯检验统计量为DG=3.51>1.96,被判定为非高斯序列,草原地区随机噪声序列的高斯统计量DG=0.37 < 1.96,被判定为高斯序列.图 5c和5d中分别比较了随机噪声的概率分布同与随机噪声具有相同均值和方差的高斯函数之间的差异,从图中可以看出,草原地区随机噪声序列的概率分布同高斯函数吻合较好,而沙漠地区随机噪声序列的概率分布同高斯函数之间存在较大差异,这说明沙漠地区随机噪声序列不服从高斯分布.对不同地区时长小于6 s的随机噪声的高斯性进行研究,结果如图 6所示.不同地区非高斯噪声序列所占百分比均高于60%,故短时长随机噪声序列应该是不服从高斯分布的,同时采集环境对随机噪声的高斯性也有影响,环境相对复杂的山地地区随机噪声的高斯性相对较好,而环境最为简单的沙漠地区随机噪声高斯性最弱.
采用线性检验法对随机噪声线性进行研究.图 7a和7b分别为两道采自沙漠地区和山地地区的随机噪声记录的时域波形,沙漠地区随机噪声记录的线性检验统计量DL=0.93 < 1.64被判定为线性时间序列,山地地区随机噪声记录的线性统计量DL=3.51>1.64被判定为非线性时间序列.DVV (Gautama et al., 2004)是一种常用的时间序列线性检验算法,DVV散点图是DVV算法的一种可视化表达.图 7c和7d为沙漠地区随机噪声序列和山地地区随机噪声序列的DVV散点图,沙漠地区随机噪声序列的特征曲线几乎与对角线重合,说明待测序列线性度很高,山地地区噪声序列的结果偏离了对角线,说明该噪声序列中存在很强的非线性成分.对不同时长随机噪声线性度进行比较,结果如图 8所示.非线性道记录百分比随噪声记录时长的增加而逐渐减小,同时,采集环境对于随机噪声线性影响作用突出,以2 s长随机噪声为例,沙漠地区91%的噪声记录被判定为线性的,而山地地区只有68%.但是总体来看,就短时长随机噪声而言,非线性道记录所占百分比较小,短时长随机噪声可以近似认为是线性系统产生的.
设计和改进地震勘探随机噪声消减算法的基础是对地震勘探随机噪声的特性有一个准确、科学的认知.在信号处理领域中,平稳性是非常重要的性质,其用来反映随机过程的时不变特性.平稳随机噪声由于噪声特性不随时间发生明显变化,在压制复杂度上较非平稳随机噪声要简单许多,同时压制算法也相对成熟.上述研究结果表明,非平稳随机噪声较平稳随机噪声而言,在中高频段具有更为丰富的能量.可以推知,随机噪声的非平稳性主要是由偶发的中高频干扰造成的,压制随机噪声的中高频能量可以改善随机噪声的平稳性.同时,地震勘探有效反射信息其频率普遍较低.因此,可以通过低通滤波器,设置恰当的通带频段,在有效保护地震反射有效信息的前提下,对随机噪声的中高频能量进行压制,达到改善随机噪声平稳性的目的.
基于上述思想,本文设计了如下的实验.选取非平稳噪声记录,应用低通滤波器对噪声记录进行处理,低通滤波器的拐点频率设置为100 Hz.图 9a为由100道非平稳实际噪声序列组成的记录,将非平稳噪声记录经过低通滤波处理后,得到处理后噪声记录如图 9b所示.通过直观判断可以发现,非平稳噪声记录中,噪声的特性时变特征明显,经过低通滤波处理后,噪声的中高频能量被压制,噪声序列特性变得相对稳定.由此,我们可以推断低通滤波压制高频干扰可以有效地改善随机噪声的平稳性.分别在两幅噪声记录中加入模拟地震反射信号,结果如图 9c和9d所示.两幅含噪记录的三条同相轴主频由上至下分别为40 Hz,30 Hz和25 Hz,偏移速度由上至下分别为1800 m·s-1,3000 m·s-1和5000 m·s-1,两幅含噪记录的信噪比均为0 dB.采用维纳滤波分别对两幅含噪记录中的噪声进行压制,结果如图 9e和9f所示.同时,分别计算了两幅去噪结果的信噪比,结果表明图 9e所示记录处理后信噪比约为4.81 dB,图 9f所示记录处理后信噪比约为7.32 dB.由维纳滤波的特性可知,当待处理噪声平稳性较好时,维纳滤波可以很好地压制噪声,取得较好的噪声消减结果;当噪声特性复杂多变时,维纳滤波的处理能力受到了限制,噪声消减结果在某些情况下可能不甚理想.两幅含噪记录在处理结果方面的差异充分说明低通滤波后噪声在平稳性上要优于原噪声记录,这同由图 9a和9b得到的关于噪声记录特性的初步判断是相一致的.同时,实验结果也说明,压制中高频干扰可以有效地改善随机噪声的平稳特性,这为后续设计随机噪声消减算法提供了一种新的思路.
风及环境噪声一直被认为是地震勘探随机噪声的主要噪声源.在短时间内,采集环境可以认为是稳定的,因此短时长的随机噪声平稳性较好;随着记录时间的增加,采集环境的稳定性难以维持,随机噪声的平稳性遭到破坏.就随机噪声产生机理而言,McNamara和Buland (2004)认为物体因为风引起震动,当这种震动通过耦合的方式传导到地表,引起地表震动,这被认为是随机噪声的主要成因,同时树通过根将震动传导入地将会产生高频噪声.山地地区的环境条件和噪声源分布情况最为复杂,随机噪声的平稳性最差,而环境最为简单的沙漠地区随机噪声平稳性最好.可以推知复杂的植被和地形条件共同作用影响随机噪声的平稳性.就高斯性而言,在理想条件下,随机噪声应该是服从高斯分布的.实际上,用于分析的随机噪声序列是由有限个噪声源生成的,这不能满足理想条件,从这方面看随机噪声的非高斯特性是合理的.山地地区随机噪声的高斯性最好,可能是由复杂的植被和地形条件造成的,这使得山地地区噪声源的分布情况更为接近理想假设.同时,由于短时长随机噪声平稳性较好,意味着短时长随机噪声具有相似的震动形式,不同时刻点之间具有较强的相关性,由此看来,它也不太可能服从高斯分布.就线性特性而言,结果表明沙漠地区随机噪声线性最好,山地地区随机噪声线性最差,可以推知随机噪声的线性特性同采集环境的复杂程度密切相关.具体的说,山地地区噪声分布情况最为多变,因此随机噪声的产生和传播机制在三个地区中最为复杂,所以在某些情况下通过简单的线性模型对随机噪声进行建模是比较困难的.
6 结论本文主要研究地震勘探随机噪声的平稳性、高斯性和线性.结果表明随机噪声不是严格平稳的,其平稳性随着噪声时长的增加而变弱,同时采集环境的植被和地形复杂程度也对随机噪声的平稳性有影响.但是在噪声时长很短的情况下,采集环境的影响就被弱化了,例如,对于3 s记录而言,至少65%的记录被判定为平稳时间序列,可以被认为是平稳的.此外,短时长地震勘探随机噪声不能认为是高斯随机过程,例如,3 s长的噪声序列超过70%被判定为非高斯序列.就线性特性而言,采集环境的复杂度对噪声的线性有着重要的影响,沙漠和草原地区随机噪声线性特性较好,而山地地区随机噪声线性特性稍差,但总体来说,短时长随机噪声还是可以近似认为是线性随机过程的产物.本文的研究成果对地震勘探随机噪声建模、改进和设计随机噪声压制方法等领域提供一定的基础条件和理论支持.
Akhouari E S, El Hassan, Laasri A, et al. 2011. Signal stationarity testing and detecting of its abrupt change.//2011 Saudi International Electronics, Communications and Photonics Conference. Riyadh:IEEE, 1-5. | |
Bath M. 1974. Spectral Analysis in Geophysics (Development in Solid Earth Geophysics). Amsterdam:Elsevier: 563. | |
Bendat J S, Piersol A G. Random Data:Analysis and Measurement Procedures.New York: Wiley Interscience Publication, 1994. | |
Brillinger D R. Time Series:Data Analysis and Theory.New York: Holt, Rinehart and Winston, 1975: 500. | |
Brockett P L, Hinich M, Wilson G R. 1987. Nonlinear and non-Gaussian ocean noise. The Journal Acoustical Society of America, 82(4): 1386-1394. DOI:10.1121/1.395273 | |
Canales L L. 1984. Random noise reduction.//Society of Exploration Geophysicists. Atlanta, USA:SEG, 525-527. | |
Chatfield C. The Analysis of Time Series:An Introduction.(6th ed).Washington, D. C: Chapman & Hall/CRC, 2003. | |
Cooper H W, Cook R E. 1984. Seismic data gathering. Proceedings of the IEEE, 72(10): 1266-1275. DOI:10.1109/PROC.1984.13016 | |
Gautama T, Mandic D P, van Hulle M M. 2004. A novel method for determining the nature of time series. IEEE Transactions on Biomedical Engineering, 51(5): 728-735. DOI:10.1109/TBME.2004.824122 | |
Groos J C, Ritter J R. 2009. Time domain classification and quantification of seismic noise in an urban environment. Geophysical Journal International, 179(2): 1213-1231. DOI:10.1111/gji.2009.179.issue-2 | |
Gulunay N. 2008. Two different algorithms for seismic interference noise attenuation. The Leading Edge, 27(2): 176-181. DOI:10.1190/1.2840364 | |
Hinich M J. 1982. Testing for Gaussianity and Linearity of a stationary time series. Journal of Time Series Analysis, 3(3): 169-176. DOI:10.1111/j.1467-9892.1982.tb00339.x | |
Hinich M J, Messer H. 1995. On the principal domain of the discrete bispectrum of a stationary signal. IEEE Transactions on Signal Processing, 43(9): 2130-2134. DOI:10.1109/78.414775 | |
Koopmans L H. The Spectral Analysis of Time Series.San Diego: Academic Press, 1995: 37. | |
Li Q Z. The Way to Obtain a Better Resolution in Seismic Prospecting(in Chinese).Beijing: Petroleum Industry Press, 1993: 112-114. | |
Lu J M. The Principles of Seismic Prospecting and Data Interpretation(in Chinese).Beijing: China University of Petroleum Press, 1991: 140-142. | |
Lu L Y, He Z Q, Ding Z F, et al. 2009. Investigation of ambient noise source in North China array. Chinese Journal of Geophysics(in Chinese), 52(10): 2566-2572. DOI:10.3969/j.issn.0001-5733.2009.10.015 | |
McNamara D E, Buland R P. 2004. Ambient noise levels in the Continental United States. Bull. Seismol. Soc. Amer., 94(4): 1517-1527. DOI:10.1785/012003001 | |
Mendel J M. 1977. White-noise estimators for seismic data processing in oil exploration. IEEE Transactions on Automatic Control, 22(5): 694-706. DOI:10.1109/TAC.1977.1101597 | |
Pan J T, Wu Q J, Li Y H, et al. 2014. Ambient noise tomography in northeast China. Chinese Journal of Geophysics(in Chinese), 57(3): 812-821. DOI:10.6038/cjg20140311 | |
Persson L. 2003. Statistical tests for regional seismic phase characterizations. Journal of Seismology, 7(1): 19-33. DOI:10.1023/A:1021216313892 | |
Shirazi S S, Moussavi Z. 2011. Investigating the statistical properties of the swallowing sounds.//Proceedings of the 33rd Annual International Conference of the IEEE Engineering in Medicine and Biology Society. Boston, MA:IEEE, 6021-6024. | |
Theiler J, Eubank S, Longtin A, et al. 1992. Testing for nonlinearity in time series:the method of surrogate data. Physica D, 58(1-4): 77-94. DOI:10.1016/0167-2789(92)90102-S | |
Wu N, Li Y, Yang B J. 2011. Noise attenuation for 2-D seismic data by radial-trace time-frequency peak filtering. IEEE Geoscience and Remote Sensing Letters, 8(5): 874-878. DOI:10.1109/LGRS.2011.2129552 | |
Yilmaz Ö, Doherty S M. Seismic Data Analysis:Processing, Inversion, and Interpretation of Seismic Data.(2nd ed).Tulsa: Society of Exploration Geophysicists, 2001. | |
Zhao P P, Chen J H, Campillo M, et al. 2012. Crustal velocity changes associated with the Wenchuan M8.0 earthquake by auto-correlation function analysis of seismic ambient noise. Chinese Journal of Geophysics(in Chinese), 55(1): 137-145. DOI:10.6038/j.issn.0001-5733.2012.01.013 | |
李庆忠. 走向精确勘探的道路.北京: 石油工业出版社, 1993. | |
陆基孟. 地震勘探原理及资料解释.北京: 中国石油大学出版社, 1991. | |
鲁来玉, 何正勤, 丁志峰, 等. 2009. 华北科学探测台阵背景噪声特征分析. 地球物理学报, 52(10): 2566–2572. DOI:10.3969/j.issn.0001-5733.2009.10.015 | |
潘佳铁, 吴庆举, 李永华, 等. 2014. 中国东北地区噪声层析成像. 地球物理学报, 57(3): 812–821. DOI:10.6038/cjg20140311 | |
赵盼盼, 陈九辉, CampilloM, 等. 2012. 汶川地震区地壳速度相对变化的环境噪声自相关研究. 地球物理学报, 55(1): 137–145. DOI:10.6038/j.issn.0001-5733.2012.01.013 | |