地球物理学报  2015, Vol. 58 Issue (12): 4576-4593   PDF    
金属矿山地地区地震勘探随机噪声的波动方程模拟
李光辉, 李月    
吉林大学信息工程系, 长春 130012
摘要: 消减随机噪声是目前陆地地震勘探数据处理的关键问题之一,分析随机噪声的产生机制及特征是对其进行有效压制的先决条件.本文针对中国南方山地金属矿区的勘探环境,根据随机噪声中包含的自然噪声和人文噪声的发声机理分别确定其噪声源函数,以波动方程作为噪声传播模型对山地地区随机噪声进行建模,将随机噪声作为一个综合波场,并且与实际噪声记录进行比较.随机噪声记录作为时空域的二维随机过程,分别对模拟噪声和实际噪声记录的时间域波形(振动图)特征包括频谱、功率谱密度,相空间轨迹图,统计量特征(能量分布,累积分布,均值,方差,峰度,偏度),和空间域波形(波剖面)特征包括波数谱和统计量特征进行比较,对比结果显示在时空域模拟噪声和实际噪声都有基本相同的性质,证明了本文对随机噪声模拟方法的可行性,为进一步研究随机噪声时空域传播特性以及噪声消除奠定理论基础.
关键词: 金属矿地震勘探     随机噪声     波动方程     波场模拟    
Wave equation modeling of random noise in seismic exploration for metal deposits in mountainous areas
LI Guang-Hui , LI Yue    
Department of Information and Engineering, Jilin University, Changchun 130012, China
Abstract: Random noise attenuation is one of key problems in seismic data processing. The existence of random noise greatly reduces the signal-noise ratio (SNR) of seismic records. Although there are lots of filtering methods available to attain this end, it is inconvenient to select a more appropriate tool for random noise attenuation, of which the characteristics change with the fields of seismic data collection. The understanding of how random noise is generated is the first requirement to solve this problem.
We model seismic random noise on land to analyze the characteristics of noise generated by different sources in seismic records. Taking the noise collected in the mountainous region in Southern China for example, the noise sources include natural sources such as wind friction over the ground surface, tree vibrations and rustles caused by wind loads, and cultural sources including running machines, footsteps of people and animals around the geophones and traffic, factories, people's daily lives in the distance. For convenience of computation, it is assumed that all of the sources contribute as point-sources in their designated areas, the function of each kind of noise source is decided according to the corresponding theory, including wind load theory, effect of mountain on wind speed, transverse vibration of beam, aeroacoustics, pseudoharmonic signal and so forth. The noise propagates by wave equation and random noise record is the superposed wave field. The theoretical model of random noise is built, the factors which influence noise characteristics are analyzed in theory, e.g. wind speed, surface roughness, mountain size etc.
When the source functions are finalized, all kinds of noise wave-fields can be obtained by solving wave equations. The synthetic records, the single channel waveforms and their frequency spectrums of each kind of noise are shown. The results show that the noise caused by branches and leaves of trees rustle in wind is the major high-frequency component. Seismic random noise is a temporal and spatial random process. As a superposed wave-field, it is composed of vibrograms and wave profiles when the distance or time is a constant. Therefore, the characteristics of the simulated noise record are compared with the real noise in the time domain (vibrogram) and space domain (wave profile), respectively, which include frequency spectrum (wave number spectrum in the space domain), power spectral density, phase locus(only in time domain), mean, variance, kurtosis, skewness, frequency distribution, and cumulative distribution function in the time domain. The comparison results both of vibrogram and wave profiles show the similarity between the simulated noise and the real one.
The comparison results demonstrate the feasibility of the proposed method. When a theoretical model of seismic random noise is built, the simulated noise in different data collection regions can be obtained by adjusting the parameters, and noise propagation characters can be analyzed in theory. The simulated noise in the corresponding regions can be used as the background noise instead of white Gaussian noise. A more suitable filtering method and its parameters can be selected and adjusted by analyzing the main component of noise and its mathematical expression.
Key words: Seismic exploration for metal deposits     Random noise     Wave equation     Wave field simulation    
 1 引言

地震勘探方法是资源勘探的主要手段之一,它以人工激发的弹性波在地层中的传播为研究基础.随着国家对矿产资源需求量的增多,以及易探、易采资源的减少,从事勘探理论与应用研究的难度越来越大.在陆地地震勘探数据采集过程中,勘探目标越深、越薄、越不规则,对地震记录质量的提高要求越大.随机噪声是严重影响地震记录信噪比的因素,由于无规律性,给地震资料的处理带来很大的困难,为此人们不断提出、改进消噪方法去除随机噪声,提高地震资料的信噪比(Wu and Yang,2011;李月等,2013Liu et al.,2013; Tian et al.,2014).其中对随机噪声的了解和认识是压制随机噪声的首要问题,人们通过不同环境下的地震噪声记录,对噪声特征值进行分析,了解各种干扰产生的原因及其特性,寻找规律性(Young et al.,1996; Bonnefoy-Claudet et al.,2006; Barajas-Olade and Ramadan,2011; 朱良保和王清东,2011赵盼盼等,2012;鲁来玉等,2014; McNamara and Buland,2004;潘佳铁等,2014),并提出压制方法.至今为止,对随机噪声的进行定性定量地分析还比较少.

随机噪声根据其噪声源类型分为自然噪声和人文噪声(Bonnefoy-Claudet et al.,2006b):自然噪声是指自然因素产生的噪声,例如风吹过时,风对地面的作用,植被枝叶的拍打声,植被随风摆动、高大 建筑物抖动引起的地面微震等(Ward and Crawford,1966). 这里将由风引起的噪声统称为风成噪声.若采集地点临近海、河流时,还有海浪拍打海岸,河水流动引起的地面震动;另外还有大规模气候变化、火山爆发等引起的远源微震.人文噪声是指人类活动引起的噪声,例如人类居住活动、工厂运作、交通噪声,以及检波器附近机器开动,人、动物走动引起的震动等等.随机噪声是一个时空域二维的随机过程,山地地区地形复杂,植被覆盖率高,给地震数据的采集带来很大的干扰,本文结合山地地区实际数据采集情况对噪声源进行模拟,用波动方程描述传播过程,将随机噪声作为各个噪声源的叠加波场,对山地地区随机噪声进行建模,并与实际噪声的时空域特征进行对比,时域波形特征包括频谱,功率谱密度,相轨迹图,统计量特征(均值,方差,能量分布,累积分布,峰度,偏度);空域波形特征包括波数谱及统计量特征.通过建立金属矿山地地区地震勘探随机噪声的理论模型,为进一步研究随机噪声的时空域传播特性以及选择合适的噪声压制方法提供理论指导作用.

2 山地地区随机噪声模拟

随机噪声由风吹草动、海浪、雷电、开动的机器、人车行走等外力产生的,其特性取决于当地的地质、地理、气象等因素以及接收点周边的环境.根据前人研究总结随机噪声中各类噪声频率分布范围如图 1所示.自然噪声中,由气候变化,季风引起的地壳运动,雷电、海洋、火山爆发等引起的微震幅度微小,因此只考虑由风作用引起的各类干扰.

图 1 随机噪声频率范围分布 Fig. 1 Frequency range of r and om noise

与非山地地区相比,山地地区地形复杂,植被覆盖情况高,尽管自然噪声源的种类大致相同,但是由于地形对风的影响导致山地地区风噪声的性质与非山地地区不尽相同.下面的章节中分别对风对地表和树木的作用力这两类噪声的源函数进行模拟,并通过波动方程对随机噪声进行建模.随机噪声是一系列噪声的叠加结果,假设噪声源以点源的形式分布在接收装置周围,如图 2.

图 2 各类噪声点源分布示意图 Fig. 2 Schematic of point source distribution of various noise
2.1 自然噪声

假设地表是均匀、各向同性的半无限大弹性介 质,风吹过地表时风荷载对地表的作用力引起地表发生形变(Sorrells et al.,1971),风荷载对树木的作用力导致树木发生振动引起的地面抖动以及风吹过时枝叶摆动发出的声音耦合到地面被接收装置检测到(Young et al.,1996; McNamara and Buland,2004).

2.1.1 风对地表的作用力引起的噪声

在风振理论中,风速可以看成是平均风速和脉动风速的和(王之宏,1994),表达式为

式中v为平均风速,vf为脉动风速.

平均风速随高度和地表粗糙度呈指数变化(王之宏,1994),表达式为

v为10 m高度的平均风速,作为参考风速;vzz(m)高度处的风速;α为地面粗糙度指数,根据建筑工程标准对地表粗糙度的分类(中华人民共和国住房和城乡建设部,2012),山地属于B类,α=0.16.

脉动风速谱采用风振理论中广泛使用的Davenport风速谱(Davenport,1961),表达式为

S(ω)为脉动风速谱,ω为脉动风角频率,$\chi $为湍流积分尺度系数且$\chi =\frac{600\omega }{\pi {{{\bar{v}}}_{10}}}$,K为地表拖拽力系数,B类地区K=0.0022.

由式(3),根据Shinozuka声振动理论(Shinozaka and Jan,1972),脉动风速时域波形为

N为分量数量,ωi为各个分量的频率,εi为各个分量的初始相位,是(0,2π)内均匀分布的随机数,Δω为频率间距,假设把频率范围(ωLωH)分为N等分,则有Δω=ωi-ωi-1=$\frac{{{\omega }_{\text{H}}}-\omega \text{L}}{N}$(阎石和郑伟,2005).

式(4)是平地脉动风速的表达式,适用于非山地地区风干扰模拟,而山地地区风速还受到山体的影响.山体对风速的影响很大,通常采用一个无量纲参数——加速效应来具体描述,加速效应(speed-up effect)指在山地地形中,某高地平均风速比平均相应高度平均风速有所增加的效应,一般来说在山顶的近地面加速效应最为明显.山体对风速影响的水平距离,一般在向风面为山高的5~10倍,背风面为15倍.山脊越高,坡度越缓,在背风面影响的距离越远.对于山体形状的模拟,主要有三角形、钟形、高斯型、余弦平方形和余弦形等(李鑫,2010).不同的山体形状,对于加速效应的影响不同,考虑到山地的复杂性可取各种模型结果的平均值或者用某一种山体模型进行研究.

本文采用余弦山体模型(图 3),模型函数为

zm为山地表面高度处,H为山顶高度,L1是山顶到H/2高度处的水平距离,s=H/L1表示山的陡缓程度,r在二维山体中等于横坐标x,在三维山体等于$\sqrt{{{x}^{2}}+C{{y}^{2}}}$,C是表示山体三维形状的常数.

图 3 余弦山体形状示意图 Fig. 3 Sketch of the cosine mountain shape

加速效应通常用一个无量纲的参数来定量描述,即加速比为

式中U(z)表示山地地面以上z高度处的平均风速,U0(z)表示平面地面以上z高度处的平均风速.

Taylor和Lee(Taylor and Lee,1984)提出的“初始算法”(“original Guidelines”)是通过ΔSmax计算不同高度处计算不同高度处的ΔS,公式为

其中AB是两个经验常数,根据山体的几何条件而有所不同,但未考虑地面粗糙度等其他因素的不同,当山体为二维连续时,A=3.5,B=1.55(李鑫,2010).

由式(2)可知,平均风速随高度和地表粗糙度呈指数变化,脉动风速均方根值在山体不同位置的分布规律与平均风速分布规律基本一致.定义脉动风速均方根的增大比ΔSσ为(李鑫,2010)

σ(z)为山地风场脉动风速均方根值,σ0(z)为平均风场脉动风速均方根值.

由于不同坡度山体背风面山脚的脉动风速均方根增大最大值均出现在0.8H处,坡度越大,脉动风速均方根增大值越大,在0.8H上下均大体呈直线分布.根据以上特点,建立了三折线型模型为(孙毅等,2011)

ΔSσmax为0.8H处的脉动风均方根最大增大值,zσ为山体影响高度,表达式为

山地地区脉动风速谱在垂向坐标上有较大差异,根据Kolmogrov理论(冯宏等,2013),频域内脉动风速功率谱为

f表示脉动风频率,σ表示摩擦速度,F表示无量纲频率,AB是两个常数,α,β,γ分别表示谱的幂指数,且满足γ-αβ=-$\frac{2}{3}$.

根据不同的实测数据对A,B,α,β,γ进行假设得到不同形式的水平脉动风速谱.对于山地地区,其湍流强度和阵风因子明显大于现用规范推荐值,而湍流积分尺度比规范建议值小很多,式(13)中参数取值分别为(田红旗,2009):A=11.4,B=14.44,α=1,β=5/3,γ=1,σ(z)取脉动风速均方根值,F=$\frac{f{{z}_{\text{ref}}}}{{{{\bar{v}}}_{10}}}$,其中zref 为参考高度,其余参数定义同上.

S(z,f)为山地脉动风速功率谱表达式,由脉动风压谱的概念及Weiner-Khintchine定理得到脉动风压功率谱为

Sp为脉动风压谱,ρ为空气密度,U(z)z高度处的平均风速,S(z,f)为脉动风速谱.

根据式(4)的理论得到山地地区脉动风压的时域波形为

各参数定义同上.

风压是平均风压和脉动风压的和,公式为

其中P是风压,Pf是脉动风压,P是平均风压,当空气密度根据统一标准取值即ρ=1.239 kg·m-3时,平均风压为P=U(z)2/1600.

当风从左向右刮时,对地面的作用力点源分布如图 2黑点所示,考虑到风的影响范围,假设风对检 波器周围地表的作用点随机均匀分布在半径为200 m 的圆形区域内,式(16)为方程(29)的震源函数,当平地参考风速 10=10 m·s-1时,通过有限差分法对方程(29)求解并叠加所有波场分量得到风噪声模拟 记录如图 4.抽取任意单道,其时域波形及谱特征如图 5所示,可以看出风吹地表引起的噪声主要能量都集中在5 Hz以内,频率较低,正如图 1中所示风对地表作用力引起的干扰小于10 Hz.

图 4 山地地区风吹过地表时引起的噪声记录 Fig. 4 Synthetic record of noise generated by force of wind on the ground surface in mountainous area

图 5 风吹过地表时引起的噪声单道(No.9)波形
(a)时域波形;(b)频谱;(c)功率谱密度.
Fig. 5 Single channel noise(No.9)waveforms
(a)Time domain waveform;(b)Frequency spectra;(c)Power spectral density.

2.1.2 树木振动引起的噪声

假设山地地区树木相对接收装置的分布如图 2星号所示.风吹过时,树木在风荷载的作用下发生振动引起地面微震.为方便计算,首先对树的模型进行简化,如图 6a,树冠面积采用抛物面,CH为树冠高 度,树冠中心为树冠高度的60%处,也是风对 树冠作用力的中心,树总高度为TH,CW为树冠宽度 CW=0.22×TH,树冠面积表达式为CA= $\frac{2}{3}$×CW×CH(Coder,2000).在顺风向,整个风力由拖曳力和升力两个风量构成,本文只考虑拖曳力的作用,根据风速风压关系公式,拖曳力为

式中CD分别为拖曳力系数,ρ为空气密度,v为风速,S为受力面积.

式中CD=0.013,由于树冠受力面积远大于树干的受力面积,因此只考虑树冠的受力作用,风对树冠的总作用力为:WP=CA×Fd.下面对树木进行一些条件假设(Bourg,2002; 王琳,2006):

(1)将树干简化为直径随高度按幂函数变化的 弹性杆,弹性杆长度为l,为地面到树冠中心的高度.

(2)弹性杆根部固定,根部直径为d0,任意高度h处直径为d(h)=d0${{\left( 1-\alpha \frac{l}{h} \right)}^{m}}$,其中α为横截面变化系数,m为树干随高度呈m次幂变化.

(3)任意高度处横截面积为A(h),大小为A(h)=A0${{\left( 1-\alpha \frac{l}{h} \right)}^{2m}}$,A0是根部面积大小,A0=$\frac{\pi }{4}$$d_{0}^{2}$.

(4)树干为均质,密度大小为ρt.

(5)任何高度处惯性矩表达式为:I(h)=I0A0${{\left( 1-\alpha \frac{l}{h} \right)}^{4m}}$,根部惯性矩大小为I0= $\frac{\pi }{64}d_{0}^{4}$.

(6)此弹性杆初始倾斜角度为θ.

根据上面给出的6条假设,把树干简化为一端固定,自由端有集中质量团的变截面弹性杆,得到树干的力学模型如图 6b所示.对树木的力学模型进行分析,考虑到一端固定一端有集中质量的弹性杆边界条件,可以得到弹性杆振动方程,用y(x,t)表示弹性杆的振动响应,由经典的梁的横向振动理论,对于一般的变截面梁,横向自由振动方程为

E为树干弹性模量,t为时间.

方程(18)的边界条件为

初始条件为
其中${{\left. \frac{\partial y}{\partial x} \right|}_{x=0}}$=θ是树干的倾斜程度,为简化计算,我们这里假设θ=0,即树是垂直于地面的.

在风荷载作用下,树的振动响应主要是受迫振 动响应.梁的受迫振动方程为

边界条件,初始条件同式(19)、(20),本文均为零.

对方程(21)进行求解得到树在风荷载作用下的振动响应作为波动方程(29)的震源函数,解在自由边界,无应力条件下的波动方程得到树振动引起的噪声.由方程(21)可知树的高度,根部直径,树冠底部直径,树冠质量,树干的弹性模量等都是影响噪声 的因素,由于文中计算的是整个植被产生的噪声,单个因素影响可以忽略,云南山地植被以阔叶林为主,树木大小取平均水平.当平地参考风速 v 10=10 m·s-1时,方程(21)的解作为方程(29)的源函数并通过有限差分法进行求解并叠加所有解得到的由于树木振动引起的噪声记录如图 7.抽取任意单道,其时域波 形及谱特征如图 8,图中显示树木风振动引起的噪声记录会产生类似于有效信号的短轴,且主要频率集中在5 Hz以内.

图 6 (a)树木的简化模型;(b)树木的力学模型 Fig. 6 (a)Simplified model of tree;(b)Mechanical model of tree

图 7 树木风振引起的噪声记录 Fig. 7 Synthetic records of noise generated by tree vibration in wind

图 8 树木振动引起的噪声单道(No.5)波形
(a)时域波形;(b)频谱;(c)功率谱密度.
Fig. 8 Single channel noise(No.5)
(a)Waveforms in time domain;(b)Frequency spectra;(c)Power spectral density.
2.1.3 风吹树木发出的声音引起的噪声

在风荷载作用下,除了引起树木振动导致的噪声,还有枝叶在风中摆动发出的声音导致的噪声.在气动声学理论中,气流流过障碍物时,由于空气的黏滞性,在障碍物背面产生气流涡流诱发的气动噪声,例如生活中常听到的风吹过电线时发出的声音,风吹过树枝、树叶时发出的声音等,如图 9所示.

图 9 风吹过障碍物时气动噪声产生示意图 Fig. 9 Schematic diagram of aerodynamic noise

Lighthill 声相似原理就是通过模拟发声源类 比气流噪声,其原理图如图 10.风速是平均风速和脉动风速的和,平均风对障碍物产生静态荷载,并不会产生气动噪声,脉动风作为随时间变化的随机变量,对障碍物产生脉动风荷载,Lighthill 类比就是将这部分脉动压力用一系列声源代替,如同障碍物表面放置了一系列的噪声源,能发出与脉动风压诱发的大小一样的气动噪声(张强,2012).物体在流体中的发声问题满足Ffowcs Williams-Hawkings方程(简称FW-H方程)(Williams and Hawkings,1969):

式中p′=p-p0为声压,p为总压力,p0为未受扰动时流场压力,pij为压应力,c为声速,vi为流体在xi向速度分量,Tij为Lighthill应力张量,ρ0为流体密度,ζ用于描述物面所在的位置,且满足无量纲条件| Δ ζ|=1,δ(ζ)为Dirac-δ函数.

方程(22)中,右边第一项是Lighthill声源项,为四极子声源项,第二项表示由表面脉动压力引起的声源(力分布),是偶极子声源项,第三项表示由加速度引起的声源(流体位移分布),是单极子声源项.通过式(17)将气流运动与声压联系起来.本文风吹树木发出的声音是由偶极子声源引起的,方程(22)右边其余两项均为零.单极子声源是构成其他两种声源的基础,是指媒质中流入的质量或热量不均匀时形成的声源(或叫简单声源),单极子声源如调制气流声源(语气声,气流扬声器),脉冲喷气等.单极子和脉动球体一样,声能量均匀的向各个方向辐射,如图 11a所示(Russell et al.,1999).单极子声源强度定义为

式中a为单极子声源的半径,U为表面振动速度,这里取脉动风速.

图 10 Lighthill 声相似原理示意图 Fig. 10 Diagram of Lighthill acoustic analogy

通过求解与方程(29)类似的无限大空间内的声波波动方程,可以得到声强为Q时的声压为

式中 |p(r,θ,t)|为单极子声压幅值且| p(r,θ,t)|= $\frac{Q\rho ck}{4\pi r}$,ρ 为空气密度,c为声传播速度,k为波数,r为 接收点到声源的距离,通常情况下满足 ka$\ll $1,kr$\gg $1,ω为角频率,t为时间.各参数定义同上,可以看出声源均匀向各个方向辐射声波,不存在指向性,如图 11b所示.

图 11 (a)单极子声源辐射,(b)单极子声源指向性 Fig. 11 (a)Acoustic radiation of monopole,(b)theoretical directivity pattern of monopole

偶极子声源是指当流体中有障碍物存在时,流体与物体产生的不稳定的反作用力形成的.偶极子声源是强度相同,方向相反的两个单极子声源组成,满足kd$\ll $1,k是波数,d是两个单极子声源之间的距离,偶极子声源是力声源,如图 12a所示.

图 12 (a)偶极子声源,(b)偶极子声源指向性 Fig. 12 (a)Acoustic radiation of dipole,(b)Theoretical directivity pattern of dipole

通过求解无限大空间内的声波波动方程,得到偶极子声源辐射声压为

式中| pd(r,θ,t)|偶极子声压幅值且| pd(r,θ,t)|= $\frac{Q\rho ck}{4\pi r}$×kd×cosθ,各参数定义同上,θ为源点接收点之间的距离向量与Z轴之间的夹角.偶极子声源辐射的声能量随cosθ变化,存在指向性(Russell et al.,1999),如图 12b所示.

采用树的简化模型,当有风吹过树冠时,由于树枝树叶受力面积的大小不同,产生了不同频率范围的风吹声,满足

式中fw为风吹声的频率,V为物体在流体中的运动速度,这里指风速,L为物体的特征长度,可以看出,风速越大,树枝、树叶越细,越小,风吹声的频率就越大.

在数据采集过程中接收到的风吹树木发出的声音是所有频率范围风吹声的叠加.由式(25)可以看出,除了风速、源点接收点之间的距离等常规参数的影响,由于偶极子声源的指向性,源点接收点之间的 距离矢量与Z轴之间的夹角也是十分重要的影响因素,本文偶极子声源指向性cosθ等于树的高度与树与接收点之间距离的比值,当参考风速v10=10 m·s-1时, 由式(25)得到各个点源引起的声波波场并将所有波场叠加得到山地地区风吹树木引起的气流噪声 记录(图 13),单道波形及谱特征(图 14),可以看出风吹树木引起的气流噪声是风成噪声高频部分的主要来源,频率范围在0~200 Hz左右,也是图 1中所示高频部分的风干扰.式(26)中L的取值导致了这部分噪声在风速确定的情况下,在某一频率能量比较集中.

图 13 风吹树木引起的气流噪声记录 Fig. 13 Synthetic record of noise generated by leaves and braches of trees rustled from wind

图 14 风吹树木引起的气流噪声单道(No.19)波形
(a)时域波形;(b)频谱;(c)功率谱密度.
Fig. 14 Single channel noise(No.19)
(a)Waveforms in time domain;(b)Frequency spectra;(c)Power spectral density.
2.2 人文噪声

人文噪声主要是由人类活动引起的噪声,根据噪声源与接收点的距离大致上分为近场噪声和远场噪声.

2.2.1 近场噪声

根据云南山地地区的实际勘探环境,近场人文噪声主要是地震测线附近机器开动,人车走动等引起的噪声,频带较窄(图 1).假设其噪声点源分布如图 2中方块所示,考虑到机器开动的连续性,其震源函数采用高斯调幅余弦时间函数为(Bonnefoy-Claudet et al.,2006b)

式中A为幅值,bw为相对带宽,t为时间变量,τ为时延,fc为中心频率,ζ为初始相位.

式(27)为波动方程(29)的噪声源函数,通过求解波动方程后将所有的解叠加得到近场人文噪声记录(图 15),单道时域波形及其谱特征(图 16),可以看出近场人文噪声频谱范围较窄,主要能量集中在20 Hz以内,正如图 1中描述的人车走动等干扰的分布范围.

图 15 近场人文噪声记录 Fig. 15 Synthetic record of near-field cultural noise

图 16 近场人文噪声单道(No.8)波形
(a)时域波形;(b)频谱;(c)功率谱密度.
Fig. 16 Single channel noise(No.8)
(a)Waveforms in time domain;(b)Frequency spectra;(c)Power spectral density.
2.2.2 远场噪声

远场噪声是采集区远处由于人类居住活动、交 通、工厂等产生的噪声,远场噪声是宽频带(0~200 Hz以上)的平稳随机噪声,可由多个(理论上是无穷多个)不同周期不同随机初相位的简谐波叠加而成(McNamara and Buland,2004),点源分布如图 2中加号所示.

远场噪声点源震源函数表达式为

j为第j个分量,a为幅值,k为角频率,ω为角频率,ε为初始相位,x为距离,t为时间变量.

远场噪声源与接收点装置的相对位置比较稳定,可近似为不随距离变化,只随时间变化的源函数,因此xl=0.将式(28)作为波动方程的源函数,得出波动方程的解并将一系列不同频率不同初始相位的解叠加得到远场人文噪声记录如图 17,单道波形及其谱特征如图 18,远场噪声幅值较小,频谱范围较宽,符合图 1中所示的人文噪声频率范围.

图 17 远场人文噪声记录 Fig. 17 Synthetic record of far-field cultural noise

图 18 远场人文噪声单道(No.23)波形
(a)时域波形;(b)频谱;(c)功率谱密度.
Fig. 18 Single channel noise(No.23)
(a)Waveforms in time domain;(b)Frequency spectra;(c)Power spectral density.
3 结果对比

上述模拟了山地地区各类噪声源引起的噪声进行模拟,这里主要将模拟噪声和实际噪声的特征进行对比,实际噪声与模拟噪声记录如图 19ab所示,可以看出模拟噪声和实际噪声记录基本相似.随机噪声作为时空域的二维随机变量,是噪声源通过波动方程传播的综合波场,因此将时间域(振动图)和空间域(剖面图)的特征逐一进行对比,以证明对随机噪声建模方法的有效性.

图 19 山地地区随机噪声记录对比
(a)实际噪声;(b)模拟噪声.
Fig. 19 Contrast of r and om noise in Yunnan mountainous area
(a)Real noise record;(b)Synthetic noise record.
3.1 振动图特征对比

振动图特征除了频谱和功率谱密度,还包括混沌性和统计量特征(包括均值、方差、峰度、偏度、频率分布以及累积分布).分别从图 19a、b中抽取任意相同道,文中抽取第6道,其时域波形及其频谱、功率谱密度对比分别如图 20abc,可以看出单道时域波形和谱特征都较为相似,频谱较宽,主要频率带在0~50 Hz左右.通过对比各种频率范围的噪声在随机噪声中所占的比重可以看到风对树木的作用力是造成山地地区高频干扰的主要来源.

图 20 单道(No.6)波形对比
(a)时域波形;(b)频谱;(c)功率谱密度.
Fig. 20 Single channel noise(No.6)
(a)Waveform comparison in time domain waveform;(b)Frequency spectra;(c)Power spectral density.

随机噪声是否存在混沌性通过其相轨迹中是否存在奇怪吸引子来判断.周期序列的相轨迹图是一个周而复始的曲线,不存在奇怪吸引子;而混沌序列的相轨迹图则是在一定区域内不规则的混乱曲线,存在奇怪吸引子,实际噪声和模拟单道时域波形的相轨迹分别如图 21ab所示,可以看出模拟噪声和实际噪声有相似的相轨迹,并且都不是简单的周期性的往复曲线,存在奇怪吸引子,具有混沌特性.

图 21 相轨迹图
(a)实际噪声;(b)模拟噪声.
Fig. 21 Phase locus comparison of the channel vibrogram
(a)Real noise;(b)Simulated noise.

模拟噪声与实际噪声振动图的频率分布如图 22ab,累积分布对比如图 22c.统计量特征对比如表 1,可以看出两者趋于一致.通过对模拟噪声和实际噪声振动图的各项性质进行对比,不难发现模拟噪声与实际噪声在时域内基本相似.

图 22 单道振动图频率与累积分布
(a)实际噪声频率分布;(b)模拟噪声频率分布;(c)累积分布对比.
Fig. 22 Frequency distribution and cumulative distribution function(CDF)of single-channel vibrogram
(a)Real noise;(b)Simulated noise;(c)CDF comparison between real and simulated noise.

表 1 模拟噪声与实际噪声振动波形统计性质对比 Table 1 Comparison of statistics of waveforms for real and simulated noise
3.2 波剖面特征对比

在地震勘探记录中,波剖面是在某一时刻由空间阵列检波器接收到的多道记录,即随机噪声的空域波形,选取任意时刻,文中取t=644 ms时,实际噪声和模拟噪声的波形对比及其波数谱如图 23ab.统计特性如图 24表 2所示.通过对模拟噪声和实际噪声波剖面的各项性质进行对比,不难发现模拟噪声与实际噪声在空域内基本相似.

图 23 t=644 ms 时刻波剖面波形对比
(a)空间域波形对比;(b)波数谱对比.
Fig. 23 Comparison of wave profiles at t=644 ms
(a)Waveforms in space domain;(b)Wave number spectrum.

图 24 t=644 ms波剖面噪声频率与累积分布
(a)实际噪声频率分布;(b)模拟噪声频率分布;(c)累积分布对比.
Fig. 24 Frequency distribution and cumulative distribution function(CDF)in profile t=644 ms
(a)Frequency distribution of real noise;(b)Frequency distribution of simulated noise;(c)CDF comparison between real and simulated noise.

表 2 模拟噪声与实际噪声波剖面统计性质对比 Table 2 Comparison of statistics of wave profiles for real and simulated noise

通过比较模拟噪声场和实际噪声场任意位置(振动图)和任意时刻(剖面图)的各项特性,结果显示了模拟噪声与实际噪声基本一致,说明了地震勘探随机噪声模拟方法的可行性.

建立随机噪声的理论模型使得对随机噪声的研究不再局限于对其性质,特征参数的总结,而是能够依据不同测区的环境,通过改变相关参数得到相应的模拟噪声,以山地地区的随机噪声为例,改变风速条件,地表地质参数,树木特征参数等,就能从非山地地区的噪声场得到山地地区噪声场,为下一步研究其传播特性和压制打下基础.首先,随机噪声的性质不是一成不变的,而是随传播介质等因素的变化发生改变,通过随机噪声建模,可以从理论上解释引起噪声性质发生变化的原因.其次,通过对随机噪声进行建模得到不同测区的模拟随机噪声后,在合成地震记录时,不再任意测区的合成记录都使用高斯白噪声作为其背景噪声,而是将该测区的模拟随机噪声当作背景噪声,为后续处理实际数据时方法选择,参数调整减小工作量.再次,能够从理论上分析不同测区随机噪声的组成部分及其特性,从而在大量的噪声压制方法中选择出适当的滤波方法.总之,对随机噪声进行理论建模,对以后的噪声压制和野外数据采集都有一定的理论指导作用.

4 结论

本文将随机噪声分为自然噪声和人文噪声两类,根据云南山地地区地震勘探数据的采集环境,自然噪声主要是风吹地表引起的地表发生形变,吹过树木时,树木发生振动以及枝叶摆动发出的声音引起的噪声,人文噪声主要是机器开动,人走动等引起的近场噪声和工厂、交通等引起的远场噪声.分别对不同的噪声源函数进行模拟,以波动方程为传播方程,通过求解非齐次波动方程得到各个噪声源引起的噪声,叠加后得到的综合波场就是时空域的随机噪声.分别将模拟噪声和实际噪声的振动图特征和波剖面特征进行对比,可以看出模拟造声与实际噪声高度相似,证明了文中模拟方法的有效性,为之后对压制随机噪声选择合适的方法打下基础.

附录

本文使用弹性波波动方程和声波波动方程.假设激发波以球面形式向外传播被检波器接收.波动方程描述了波传播的基本规律,通过波动方程解表征随机噪声的波场特征.根据波传播理论,波在理想介质中传播时,满足方程:

式中${{\nabla }^{2}}=\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{z}^{2}}}$,$\phi $为位移势,c为波传播速度,t为时间,r 为距离矢量,M(r,t)为震源函数.

根据弹性动力学理论,弹性波在均匀、各向同性、半无限大的理想介质中传播时,假设震源点位于 地表,Z轴垂直向下,方程(29)满足无应力边界条件为

式中λ、μ为拉梅常数,ux、uz分别为X方向、Z方向的位移分量且可以表示为
式中φ、ψ分别为横波和纵波的标量位.

声波在无限大理想介质中传播时满足自由边界条件.得到各个噪声源函数后选择合适的方法(本文选取有限差分法)求解波动方程,得到不同噪声源下的噪声波形,叠加后得到的综合波场即为随机噪声.

参考文献
[1] Barajas-Olade C, Ramadan A. 2011. Is it possible to conduct seismic wind noise experiments in a wind tunnel?. //Presented at the 73rd EAGE Conference and Exhibition Incorporating SPE EUROPEC. Vienna, Austria: European Association of Geoscientists and Engineers, 238.
[2] Bonnefoy-Claudet S, Cotton F, Brad P Y. 2006a. The nature of noise wavefield and its applications for site effects studies: A literature review. Earth-Science Reviews, 79(3-4): 205-227.
[3] Bonnefoy-Claudet S, Cotton F, Brad P Y. 2006b. H/V ratio: a tool for site effects evaluation. Results from 1-D noise simulations. Geophys. J. Int., 167(2): 827-837.
[4] Bourg D M. 2002. Physics for Game Developers. Beijing: Publishing House of Electronics Industry, 135-137.
[5] Coder K D. 2000. Estimating wind forces on tree crowns. Athens, Georgia: The University of Georgia, http://warnell.forestry.uga.edu/service/library/for00-016/for00-016.pdf.
[6] Davenport A G. 1961. The spectrum of horizontal gustiness near ground in high winds. Roy. Meteor. Soc., 87(372): 194-211.
[7] Feng H, Xiao Z Z, Li Z L, et al. 2013. Research on the wind speed spectrum of complex mountainous environment. Journal of Hunan University (Natural Sciences) (in Chinese), 40(1): 27-32.
[8] Li X. 2010. Study on characteristics of wind field in ground layer on hilly terrain (in Chinese)[Ph. D. thesis]. Chongqing: Chongqing University.
[9] Li Y, Peng J L, Ma H T, et al. 2013. Study of the influence of transition IMF on EMD do-noising and the improved algorithm. Chinese J. Geophys. (in Chinese), 56(2): 626-634, doi: 10.6038/cjg20130226.
[10] Liu Y P, Li Y, Nie P F, et al. 2013. Spatiotemporal time-frequency peak filtering method for seismic random noise reduction. IEEE Geoscience and Remote Sensing Letters, 10(4): 756-760.
[11] Lu L Y, He Z Q, Ding Z F, et al. 2014. Investigation of ambient seismic noise sources in the North China array. Chinese J. Geophys. (in Chinese), 57(3): 822-836, doi: 10.6038/cjg20140312.
[12] McNamara D E, Buland R P. 2004. Ambient noise levels in the Continental United States. Bulletin of the Seismological Society of America, 94(4): 1517-1527.
[13] Ministry of Housing and Urban-Rural Development of the People's Republic of China (MOHURD). 2012. GB 50009-2012 Load Code for the Design of Building Structures (in Chinese). Beijing: China Building Industry Press.
[14] Pan J T, Wu Q J, Li Y H, et al. 2014. Ambient noise tomography in northeast China. Chinese J. Geophys. (in Chinese with), 57(3): 812-821, doi: 10.6038/cjg20140311.
[15] Russell D A, Titlow J P, Bemmen Y J. 1999. Acoustic monopoles, dipoles, and quadrupoles: An experiment revisited. Am. J. Phys., 67(8): 660-664.
[16] Shinozaka M, Jan C M. 1972. Digital simulation of random processes and its applications. Journal of Sound and Vibration, 25(1): 111-128.
[17] Sorrells G G, McDonald J A, Der Z A, et al. 1971. Earth motion caused by local atmospheric pressure changes. Geophys. J. Roy. Astr. Soc., 26(1-4): 83-98.
[18] Sun Y, Li Z L, Huang H J, et al. 2011. Experimental research on mean and fluctuating wind velocity in hilly terrain wind field. Acta Aerodynamica Sinica (in Chinese), 29(5): 593-599.
[19] Taylor P A, Lee R J. 1984. Simple guidelines for estimating wind speed variations due to small scale topographic features. Climatol. Bull,18(2): 3-22.
[20] Tian H Q. 2009. Advances in Industrial Aerodynamics (in Chinese). Changsha: Central South University Press, 20-23.
[21] Tian Y N, Li Y, Yang B J. 2014. Variable-eccentricity hyperbolic-trace TFPF for seismic random noise attenuation. IEEE Transactions on Geoscience and Remote Sensing, 52(10): 6449-6458.
[22] Wang L. 2006. Building and analysis of dynamic model for windthrow of spruce (in Chinese)[Ph. D. thesis]. Harbin: Harbin Institute of Technology.
[23] Wang Z H. 1994. Simulation of wind loading. Journal of Building Structures (in Chinese), 15(1): 44-52.
[24] Ward H S, Crawford R. 1966. Wind-induced vibrations and building modes. Bull. Seismol. Soc. Am., 56(4): 793-813.
[25] Williams J E F, Hawkings D L. 1969. Sound generation by turbulence and surfaces in arbitrary motion. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 264(1151): 321-342.
[26] 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 letter, 8(5): 874-878.
[27] Yan S, Zheng W. 2005. Wind load simulation by superposition of harmonic). Journal of Shenyang Arch. and Civ. Eng. Univ. (in Chinese), 21(1): 1-4.
[28] Young C J, Chael E P, Withers M M, et al. 1996. A comparison of the high-frequency ( > 1 Hz) surface and subsurface noise environment at three sites in the United States. Bull. Seismol. Soc. Am., 86(5): 1516-1528.
[29] Zhang Q. 2012. Aeroacoustics Foundation (in Chinese). Beijing: National Defence Industry Press, 96-142.
[30] 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 J. Geophys. (in Chinese), 55(1): 137-145, doi: 10.6038/j.issn.0001-5733.2012.01.013.
[31] Zhu L B, Wang Q D. 2011. An expression of the cross-correlation of ambient seismic noise: A derivation based on the surface-wave theory. Chinese J. Geophys. (in Chinese), 54(7): 1835-1841, doi: 10.3969/j.issn.0001-5733.2011.07.017.
[32] 冯宏,肖正直,李正良等. 2013. 复杂山地环境下脉动风速谱研究. 湖南大学学报(自然科学版), 40(1): 27-32.
[33] 李鑫. 2010. 山地地形的近地风场特性研究[硕士论文]. 重庆: 重庆大学.
[34] 李月, 彭蛟龙, 马海涛等. 2013. 过渡内蕴模态函数对经验模态分解去噪结果的影响研究及改进算法. 地球物理学报, 56(2): 626-634, doi: 10.6038/cjg20130226.
[35] 鲁来玉, 何正勤, 丁志峰等. 2014. 基于背景噪声研究云南地区面波速度非均匀性和方位各向异性. 地球物理学报, 57(3): 822-836, doi: 10.6038/cjg20140312.
[36] 潘佳铁, 吴庆举, 李永华等. 2014. 中国东北地区噪声层析成像. 地球物理学报, 57(3): 812-821, doi: 10.6038/cjg20140311.
[37] 孙毅, 李正良, 黄汉杰等. 2011. 山地风场平均及脉动风速特性试验研究. 空气动力学学报, 29(5): 593-599.
[38] 田红旗. 2009. 工业空气动力学研究进展. 长沙: 中南大学出版社, 20-23.
[39] 王琳. 2006. 云杉风倒动力学模型的建立与分析[硕士论文]. 哈尔滨: 哈尔滨工业大学.
[40] 王之宏. 1994. 风荷载的模拟研究. 建筑结构学报, 15(1): 44-52.
[41] 阎石, 郑伟. 2005. 简谐波叠加法模拟风谱. 沈阳建筑大学学报(自然科学版), 21(1): 1-4.
[42] 张强. 2012. 气动声学基础. 北京: 国防工业出版社, 96-142.
[43] 赵盼盼, 陈九辉, Campillo M等. 2012. 汶川地震区地壳速度相对变化的环境噪声自相关研究. 地球物理学报, 55(1): 137-145, doi: 10.6038/j.issn.0001-5733.2012.01.013.
[44] 中华人民共和国住房和城乡建设部. GB 50009-2012建筑结构荷载规范. 北京: 中国建筑工业出版社.
[45] 朱良保, 王清东. 2011. 地震背景噪声互相关函数的面波理论表达形式. 地球物理学报, 54(7): 1835-1841, doi: 10.3969/j.issn.0001-5733.2011.07.017.