地球物理学报  2019, Vol. 62 Issue (2): 680-696   PDF    
利用EMD和子波振幅谱与相位谱关系的时变子波提取方法
张鹏, 戴永寿, 谭永成, 张红倩, 王春娴     
中国石油大学(华东)信息与控制工程学院, 山东青岛 266580
摘要:本文首先分析了地震波在黏弹介质的传播规律,基于黏弹介质地震波动方程总结了时变子波振幅谱和相位谱的关系,从而得出结论,准确估计子波相位谱初值和不同时刻的子波振幅谱是实现时变子波准确提取的必要条件.在此基础上,针对传统方法限制子波振幅谱形态且受限于分段平稳假设的问题,提出了一种利用EMD(Empirical Mode Decomposition)和子波振幅谱与相位谱关系的时变子波提取方法,根据子波对数振幅谱光滑连续而反射系数对数振幅谱振荡剧烈的特点,采用EMD方法将不同时刻地震记录的对数振幅谱分解为一组具有不同振荡尺度的模态分量,通过滤除振荡剧烈分量、重构光滑连续分量提取时变子波振幅谱;再应用子波振幅谱和相位谱的关系提取时变子波相位谱,将分别提取的振幅谱和相位谱逐点进行合成,最终实现时变子波的准确提取.本文方法不需要求取Q值,适用于变Q值的情况,具有良好的抗噪性能.数值仿真和叠后实际资料处理结果表明,相比传统的分段提取方法,利用本文方法提取的时变子波准确度更高,研究成果对提高地震资料分辨率具有重要意义.
关键词: 时变子波提取      EMD      振幅谱      相位谱      地震波波动方程     
A time-varying wavelet extraction method using EMD and the relationship between wavelet amplitude and phase spectra
ZHANG Peng, DAI YongShou, TAN YongCheng, ZHANG HongQian, WANG ChunXian     
College of Information and Control Engineering, China University of Petroleum(East China), Qingdao Shandong 266580, China
Abstract: The propagation and attenuation law of seismic wavelets is analysed, and the relationship between the amplitude and phase spectra of time-varying wavelets is obtained by deriving the wave equation for a viscoelastic medium. We find that accurately estimating the amplitude spectrum at different time points and the initial phase spectrum is necessary to extract accurate time-varying wavelets. Based on this conclusion, to overcome the defects of the classical methods, we propose a time-varying wavelet extraction method that utilizes EMD (Empirical Mode Decomposition) and the relationship between the wavelet amplitude and phase spectra. According to the differences that the logarithm amplitude spectra of wavelets and reflection coefficient are smooth and oscillating respectively, the logarithm amplitude spectra of the seismogram at different time points are decomposed into multi-layer components with different oscillation scales by EMD, and the amplitude spectra of time-varying wavelets can be estimated by filtering the oscillating components and reconstructing the smooth components. Then the phase spectra of time-varying wavelets are estimated by applying the obtained relationship, so that time-varying wavelet extraction is achieved by combining the estimated amplitude and phase spectra of wavelets correspondingly. This method does not need to calculate the Q value and can be applied to the case of variable Q, exhibiting good anti-noise performance. The numerical simulation and real seismic data processing results demonstrate that the proposed method can improve the accuracy of time-varying wavelet extraction compared to the classical method.
Keywords: Time-varying wavelet extraction    EMD    Amplitude spectrum    Phase spectrum    Seismic wave equation    
0 引言

地震子波是影响地震记录分辨率的重要因素,传统的地震子波提取方法基于子波时不变的假设,然而由于地下介质对地震波能量的吸收衰减作用,实际地震资料中的地震子波具有时变特征,即地震资料具有非平稳性.时变子波的准确提取是实现高精度地震资料处理和解释的必要条件,也是目前迫切需要解决的重要问题,同时时变子波提取的方法研究也将成为推动正反演、偏移、反褶积等地震资料处理和解释技术进一步发展的“助推剂”.

近年来,为了进一步提高地震资料处理精度,国内外学者针对非平稳地震记录处理、时变子波提取开展了一系列研究工作.Clarke(1968)首先考虑了地震记录的非平稳特性,提出了非平稳褶积模型;Margrave(1998)王守君等(2015)在此基础上分别论证了地层的吸收和频散效应是地震子波振幅和相位具有时变特征的根本原因;王晓华等(1998)揭示了时变子波的有效提取对于谱反演高分辨率反射系数的重要意义.Stockwell(2007)Sinha等(2009)采用时频分析方法补偿地震记录的能量衰减和频率损失;Wang等(2014)提出了基于贝叶斯反演的地层衰减补偿方法,以实际地震资料为先验信息实现了对非平稳地震记录的补偿.但这类方法在实际应用中易出现欠补偿或过补偿等问题,影响了子波提取精度.Rosa和Ulrych(1991)最先提出谱模拟反褶积方法,Zhou等(2014)将该方法进一步推广,利用时频谱模拟方法提取时变子波振幅谱,有效反映了时变子波能量衰减、频率损失等特征.Wang(2002, 2006)提出了稳定的反Q滤波方法,可以同时实现子波振幅的补偿和相位的校正.该方法必须已知地层品质因子Q,但在实际资料处理中往往难以准确提取Q值.van der Baan(2008, 2012)将非平稳地震记录分解成若干等长的重叠片段,在子波常相位的假设下,在每段内提取子波平均相位谱,再延拓到每一采样点,实现时变子波相位的提取和校正;Economou和Vafidis(2010, 2012)分段提取时变子波并基于提取的子波进行反褶积处理以实现衰减补偿和相位校正;Wang等(2010, 2012)提出了一种自适应分子分解的方法,无须人为指定窗函数的数据长度和分段数目,根据地震记录特性自动进行分段,进一步优化了分段子波提取方法.但现有时变子波提取方法对子波衰减规律的假设仍然存在局限性,时频谱模拟方法限制了子波振幅谱的形态,分段方法受限于分段平稳假设,导致提取的子波难以真实反映地震资料的时变特征.

因此,为进一步提高时变子波提取精度,本文深入分析了地震波在黏弹介质中的传播规律,针对叠后地震记录从黏弹介质地震波波动方程中总结推导得出了时变子波振幅谱和相位谱的相互关系,为时变子波间接提取的方法研究提供理论依据和指导.在此基础上,针对现有方法存在的问题,引入时变信号处理中的EMD方法,提出一种利用EMD和子波振幅谱与相位谱关系的时变子波提取方法,在不求取Q值的情况下就能实现时变子波的提取.理论分析、仿真实验和实际资料处理结果表明,相比于传统的分段提取方法,利用本文方法提取的时变子波更符合地震波的衰减规律,因此提取精度更高.同时通过采用时频滤波技术有效提高了方法的抗噪性能,即在噪声干扰下依然能够准确地提取时变地震子波.

1 黏弹介质中子波振幅谱和相位谱的关系

黏弹介质中地震波的波动方程在频率域可以表示为

(1)

该方程基于品质因子Q与频率无关的假设.其中,U表示地震波的谱,z表示地震波的传播距离,ω表示波的角频率,k(ω)为z方向的波数,它可以由角频率ω、地震波相速度v(ω)和品质因子Q表示:

(2)

其中,j表示虚数符号.(1)式经过微分算子分解,可得单程波方程(Wang, 2006):

(3)

(4)

(3) 式代表前进方向的波,(4)式代表后退方向的波.由于本文研究的是非平稳地震记录的形成机理及时变子波的衰减规律,因此主要对(3)式进行研究.对(3)式进行积分,可得方程的解析解为:

(5)

其中,Δz表示深度增量.将(2)式代入(5)式,可得到地震子波在z方向的频域表达式:

(6)

则地震子波在深度域的振幅谱A和相位谱φ可分别如(7)式和(8)式所示:

(7)

(8)

接下来,将(7)和(8)式分别从深度域变换到时间域.首先将(7)式等号两边同时取对数,可得

(9)

参考Rickett进行深-时转换时的假设(Rickett, 2006),定义时间间隔Δtz/v(ω),则(9)式在时间域即可表示为

(10)

其中,A(tt, ω)和A(t, ω)分别表示不同时刻地震子波的振幅谱.将(10)式中的-ωΔt/Q定义为地层衰减因子,则(10)式描述了时变子波振幅谱和地层衰减因子的关系.

然后,对(8)式进行深-时转换.参考Rickett对时间间隔Δt的定义(Rickett, 2007; Braga and Moraes, 2013),使得Δtz/v(ωr),v(ωr)为参考频率ωr的相速度,而ωr为地震记录频带的最高频率.(8)式在时间域可以表示为

(11)

Aki和Richards(1980)给出了v(ω)/v(ωr)的麦克劳林展开式:

(12)

将(12)式代入(11)式可得时变子波相位谱和地层衰减因子的关系如(13)式所示:

(13)

由(10)式和(13)式可知,时变子波振幅谱和相位谱分别与地层衰减因子具有相关关系.因此,以地层衰减因子为中间量,即可得到黏弹介质中时变子波振幅谱和相位谱之间的关系如(14)式所示:

(14)

式中,A(tt, ω)、A(t, ω)、φ(t, ω)和φ(tt, ω)为四个未知量.从中可以得出结论,估计出不同时刻子波的振幅谱A(tt, ω)和A(t, ω)以及某一时刻子波的相位谱φ(t, ω),即可在不求取Q值的情况下得到各个时刻时变子波的振幅谱和相位谱,从而实现时变子波提取.

在此基础上,本文提出了一种利用EMD和子波振幅谱与相位谱关系的时变子波提取方法,首先采用EMD方法提取不同时刻时变子波振幅谱,然后利用(14)式所示的子波振幅谱和相位谱的关系,提取出不同时刻时变子波的相位谱,从而实现时变子波的提取.

2 时变子波提取方法 2.1 时频分析和时频滤波

非平稳褶积模型可表示为(Margrave et al., 2011)

(15)

*表示卷积,x(t)为非平稳地震记录,w(t, v)表示时变子波,r(t)为反射系数序列.

非平稳地震记录具有非平稳信号的特征,而时频分析是表征信号在每一时间点处的频率成分,是研究非平稳信号的有力工具.故对地震记录进行时频分析能够充分描述其中的非平稳特性.采用改进的广义S变换将非平稳地震记录变换到时频域(陈学华等, 2009),其窗函数的表达式为

(16)

其中q, p为大于零的调节因子.地震记录x(t)的时频谱为

(17)

野外采集的地震资料信噪比较低,虽然叠前进行了各种噪声压制,但对于一些能量相对较弱的噪声难以识别和彻底压制,因此叠后地震记录中仍存在一定量的噪声,且主要表现为随机噪声,可认为幅值较小,在频带内分布均匀.由于地震记录的振幅谱表现为带限信号,所以随机噪声将主要影响地震记录振幅谱的高频部分,将对时变子波提取结果产生严重的影响.为了提高时变子波提取方法的抗噪性能,对地震记录的时频谱进行时频滤波,即在不同时刻地震记录的频谱中分别加入巴特沃斯低通滤波器H(t, f),巴特沃斯低通滤波器的幅频特性如图 1所示,其中fp, fs分别为通带截止频率和阻带截止频率.根据巴特沃斯低通滤波理论,fp对应信号峰值的0.707倍(对应-3 dB)处的频率,fs对应信号峰值的0.001倍(对应-60 dB)处的频率,因此fp, fs的取值能够根据不同时刻地震记录的频谱特征和信噪比自动确定,即可以实现自适应滤波.经时频滤波后地震记录的时频谱可表示为

(18)

图 1 低通滤波器的幅频特性 Fig. 1 Amplitude frequency characteristic of the low pass filter
2.2 基于EMD方法提取时变子波振幅谱 2.2.1 EMD方法的基本原理

EMD方法由美国Huang(1998)提出,其主要的优点在于基函数可以从信号自身获得,克服了小波变换中选择基函数的困难.根据信号时间尺度的不同,EMD可以将复杂的信号分解成若干个按振荡程度由高到低排列的固有模态函数(Intrinsic Mode Function,IMF),可以利用这个性质对信号进行分解和分离.每个IMF必须满足两个条件:一是在整个数据区内,局部极大值和局部极小值的数目之和与过零点数目必须相等或至多相差一个;二是由极大值与极小值确定的包络线均值为零.当仅剩一个趋势分量时,分解停止(Han and van der Baan, 2013).对信号x(t)进行EMD分解的基本步骤为:

(1) 求取原始信号x(t)的局部极大值和局部极小值,拟合生成上、下包络线;

(2) 计算上、下包络线的均值;

(3) 用原始信号减去步骤(2)中的均值;

(4) 判断是否满足IMF的两个筛选条件,若满足,进入下一步;若不满足,用步骤(3)中结果返回步骤(1)进一步计算;

(5) 分解结果作为第i个IMF分量;

(6) 判断是否满足EMD分解终止条件,若满足转入步骤(8),若不满足转入步骤(7);

(7) 从原始信号中减去该层IMF分量作为新的原始信号返回步骤(2),计算第i+1层IMF分量;

(8) 得到n层IMF分量和剩余分量.

经分解得到一系列IMF分量后可表示为

(19)

其中,IMFi为第i个IMF分量,c(t)为分解剩余分量.

针对EMD存在的模态混叠问题,本文在分解过程中叠加高斯白噪声,利用白噪声具有频率均匀分布的统计特性来消除模态混叠效应(Wu and Huang, 2009).具体步骤如下:

(1) 给原始信号x(t)添加一组白噪声(标准差取被分析信号标准差的0.1~0.4倍),得到信号s(t);

(2) 对s(t)进行EMD分解,得到IMFs;

(3) 重复以上(1)、(2)两步,每次添加不同的噪声序列,分解后得到各自的IMF分量组;

(4) 相应IMF的均值为分解的最终结果.

2.2.2 时变子波振幅谱的提取

根据(15)式所示的非平稳褶积模型,非平稳地震记录可以表示为时变子波和反射系数序列的褶积,根据信号处理理论中时域褶积则频域相乘的关系(吴大正, 1998),对(15)式等式两边同时进行时频变换,则地震记录的时频谱可以表示为时变子波与反射系数序列时频谱的乘积.同时,Margrave等(2011)经过研究,得到如下表达式:

(20)

其中,表示非平稳地震记录经时频变换得到的时频谱,表示反射系数序列经时频变换得到的时频谱,表示震源子波的傅里叶谱,α(t, f)表示衰减因子,则表示传播过程中的时变子波在时刻t的傅里叶谱,即时变子波的时频谱.本文用W(τ, f)表示式(20)中的,即时变子波的时频谱,用X(τ, f)和R(τ, f)分别表示非平稳地震记录和反射系数序列的时频谱,则有

(21)

X(τ, f),W(τ, f)和R(τ, f)分别反映了地震记录、子波和反射系数的频谱随时间的变化情况.根据信号处理理论对非周期信号频谱的定义,信号的频谱可以由它的振幅谱和相位谱表示(吴大正, 1998; Maly and Mowlaee, 2016)

(22)

式中,|F(jω)|和φ(ω)分别是频谱函数F(jω)的模和相位,ω=2πf.因此在时频域利用振幅谱和相位谱可将(21)式表示为(Zhou et al., 2014; Dai et al., 2016)

(23)

其中,j代表虚数符号,|X(τ, f)|,|W(τ, f)|和|R(τ, f)|分别表示地震记录、时变子波和反射系数的振幅谱,而φx(τ, f),φw(τ, f)和φr(τ, f)分别表示地震记录、时变子波和反射系数的相位谱.对(23)式等式两边同时取对数,则地震记录、时变子波和反射系数的对数振幅谱可以表示为

(24)

其中,ln|X(τ, f)|表示地震记录的对数振幅谱.ln|R(τ, f)|为反射系数的对数振幅谱,具有振荡剧烈的特点,表征为高频信号.ln|W(τ, f)|为时变子波的对数振幅谱,具有光滑连续的特点,表征为低频信号.因此可以将地震记录的对数振幅谱看作含噪声信号,反射系数的对数振幅谱看作加性噪声,而将时变子波的对数振幅谱看作待估计的有效信号.这一特性符合如(25)式所示的含有噪声信号模型.

(25)

其中,s(t),x(t)和n(t)分别表示含噪声信号、有效信号和噪声信号.因此,时变子波振幅谱的提取问题即可转化为信号去噪问题.

信号去噪的目标是从含噪声数据中恢复有效信号,滤除噪声干扰.维纳滤波、基于傅里叶变换的低通滤波等方法仅能应用于平稳信号,而针对非平稳信号的去噪问题主要采用小波分析方法,但小波分析须预先定义好小波基,缺乏自适应性,影响了实际应用效果.而与小波分析相比,EMD方法的自适应性体现在分解过程中,它的基函数直接从信号本身产生,不同的信号会产生不同的基函数,EMD方法还具有自适应滤波特性和自适应分辨率特性(Gan et al., 2014).EMD方法的完备性体现在整个分解过程没有能量的损失,可以用分解得到的分量进行信号恢复和重构(王婷, 2010).基于上述优势,EMD近年来被广泛采用以解决信号去噪问题,取得了良好的应用效果,相比传统方法去噪精度更高.

综上所述,EMD方法能够用于提取时变子波振幅谱,主要理论依据是:将反射系数的对数振幅谱视为叠加在时变子波对数振幅谱上的噪声干扰,利用EMD的自适应分解和分析能力,将地震记录的对数振幅谱分解为一组具有不同振荡尺度的分量,振荡剧烈的高频分量对应着反射系数的对数振幅谱(噪声分量),光滑连续的低频分量对应着时变子波的对数振幅谱(有效信号分量),通过将光滑连续的分量进行重构即可实现时变子波振幅谱的准确提取.因此,本文利用EMD方法将不同时刻地震记录的对数振幅谱分解为一组具有不同振荡尺度的固有模态函数(IMF)分量IMF1, IMF2, …, IMFn,反射系数的振幅谱对IMF分量的支配作用逐渐降低,子波振幅谱对IMF分量的支配作用逐渐增强,即IMF分量中反射系数的振幅能量逐渐减弱而子波的振幅能量逐渐增强,然后用光滑连续的几个IMF对子波对数振幅谱进行重构,即:

(26)

为时变子波对数振幅谱的估计值.因此,必然存在某个IMF分量,使得该分量之后子波支配的振幅能量超过反射系数支配的振幅能量.时变子波振幅谱提取的目标就是找到这个分界IMF分量的索引值k=js,使得利用该分量以后的IMF分量进行重构的结果误差最小.为了实现这个目标,Boudraa和Cexus(2007)提出了连续均方误差(CMSE)的准则,即:

(27)

其中,N为信号的长度.基于该准则,索引值js可以由(28)式给出.

(28)

由于反射系数的振幅能量主要集中在振荡剧烈的分量中,随着分解的进行,反射系数的能量逐渐减小,于是可以将IMF能量首次发生转折的位置作为反射系数起主导作用的模态与子波起主导作用的模态的分界分量.

综上所述,利用EMD提取时变子波振幅谱的步骤可以归结如下:

(1) 将时频变换和滤波处理后不同时刻地震记录的对数振幅谱进行EMD分解得到一组具有不同振荡尺度的IMFk, k=1, …, n-1以及剩余分量c(t);

(2) 利用(26)式计算, k=1, …, n-1;

(3) 利用(27)式计算, k=1, …, n-1;

(4) 利用(28)式计算js的值;

(5) 利用(26)式对信号进行重构,得到时变子波的对数振幅谱,再经过指数变换即可提取时变子波振幅谱.

传统的基于谱模拟的时变子波振幅谱提取方法须预先对子波振幅谱进行定义,只有在振幅谱的形态满足定义式时,提取结果才准确,因此应用范围受限.而本文方法充分利用了EMD自适应性强的优势,对不同时刻地震记录的对数振幅谱分别进行分解和分析,无须预先定义和假设,自动地滤除高频分量,重构低频的光滑连续分量实现子波振幅谱的准确提取,克服了谱模拟方法的缺陷,对任意形态的子波振幅谱均可有效估计,从而提高了时变子波振幅谱的提取精度.

2.3 利用子波振幅谱和相位谱的关系提取时变子波相位谱

利用时频谱模拟方法提取得到了不同时刻子波的振幅谱A(tt, ω)、A(t, ω),在此基础上继续估计子波相位谱初值φ(t, ω),即可利用(14)式实现各个时刻时变子波振幅谱和相位谱的准确提取.

高阶累积量和高阶谱由于包含了信号的完整信息被地球物理工作者引入到混合相位子波提取中.根据高阶统计量的性质,通过地震记录的双谱可以重构出子波的相位谱初值.现有的相位重构方法有很多种,其中Matsuora-Ulrych算法(Matsuora and Ulrych, 1984)不存在累积误差,且运算简单,故本文选用此算法进行子波相位谱初值的估计.

双谱是信号三阶累积量的傅里叶变换.地震记录的三阶累积量可以表示为

(29)

通过傅里叶变换得到地震记录的三阶谱,即双谱,记为Bx(ω1, ω2),如(30)式所示:

(30)

根据地震记录、反射系数和地震子波的关系,Bx(ω1, ω2)又可以表示为

(31)

其中,Br(ω1, ω2)和Bω(ω1, ω2)分别为反射系数和地震子波的双谱.根据文献(Matsuora and Ulrych, 1984),Br(ω1, ω2)=E[r3(t)]=C≠0,(31)式可以表示为

(32)

其中,|Bx(ω1, ω2)|为地震记录双谱的振幅谱,ψx(ω1, ω2)为地震记录双谱的相位谱,W(ω)=|W(ω)|exp(iφ(ω)),W(ω), |W(ω)|和φ(ω)分别表示地震子波的频谱、振幅谱和相位谱,因此地震记录双谱的相位谱和地震子波的相位谱可以如(33)式所示:

(33)

利用(33)式即可从地震记录的双谱中提取子波的相位谱初值.

将估计的不同时刻子波的振幅谱和相位谱初值代入式(14)即可得到所有时刻子波的振幅谱和相位谱,将振幅谱和相位谱逐点在频域进行合成,再通过逆傅里叶变换即可得到时变子波的时域表达.各个时刻提取的时变子波可表示为

(34)

其中,为提取的时变子波,F-1表示逆傅里叶变换.综上所述,利用子波振幅谱和相位谱的关系提取时变子波的主要流程如图 2所示.

图 2 利用EMD和子波振幅谱与相位谱关系提取时变子波的方法流程 Fig. 2 Flowchart of the time-varying wavelet extraction method using EMD and the relationship between wavelet amplitude and phase spectra

传统的分段子波提取方法是在分段平稳的假设条件下分段提取子波,再延拓到各个时刻实现时变子波的提取,该方法受限于分段平稳的假设,同时延拓结果难以符合实际的衰减特征,影响了提取结果的准确性.而本文方法无需分段平稳的假设,并依据地震波在黏弹介质中的衰减规律精细化地提取时变子波,克服了传统方法的缺陷,提高了时变子波提取精度.

3 数值仿真验证

为了验证利用子波振幅谱和相位谱相关关系提取时变子波的有效性,本文采用ARMA模型描述的混合相位子波,与满足独立同分布及伯努利-高斯分布,长度1000 ms,采样间隔1 ms的反射系数序列,基于如(15)式所示的非平稳褶积模型构造非平稳地震记录,如图 3所示.初始混合相位子波的振幅谱和相位谱如图 4所示.实验中采用的初始子波其ARMA描述下的差分形式为

(35)

图 3 合成地震记录 (a)混合相位子波;(b)反射系数序列;(c)非平稳地震记录
(t=1~1000 ms, Q=150).
Fig. 3 Synthetic seismogram (a) Mixed-phase wavelet; (b) Reflection coefficient sequence; (c) Non-stationary seismogram (t=1~1000 ms, Q=150).
图 4 初始混合相位子波的振幅谱和相位谱 (a)子波振幅谱;(b)子波相位谱. Fig. 4 Amplitude and phase spectra of the initial mixed-phase seismic wavelet (a) Amplitude spectrum of the wavelet; (b) Phase spectrum of the wavelet.

其Z域系统函数为

(36)

3.1 时变子波提取测试

采用EMD方法从图 3c所示的非平稳地震记录中提取时变地震子波振幅谱,将提取结果与时频谱模拟方法提取的振幅谱、理论子波的振幅谱及地震记录的振幅谱进行对比如图 5所示.从中可以看出,相比于时频谱模拟方法,利用EMD方法提取的时变子波振幅谱在低频部分的拟合精度得到了提高,从而进一步提高了时变子波振幅谱的提取精度.

图 5 时变子波振幅谱提取结果 虚线表示理论值,点划线表示谱模拟方法的提取结果,实线表示EMD方法的提取结果,点线表示地震记录的振幅谱. (a) 111 ms处的提取结果;(b) 381 ms处的提取结果;(c) 619 ms处的提取结果;(d) 892 ms处的提取结果. Fig. 5 Results of time-varying wavelet amplitude spectra extraction Theoretical values are expressed by the dashed line, and the extraction results using spectral modeling method are expressed by the dot dash line, and the extraction results using EMD method are expressed by the solid line, and the amplitude spectra of the seismogram are expressed by the dotted line. (a) Extraction results at 111 ms; (b) Extraction results at 381 ms; (c) Extraction results at 619 ms; (d) Extraction results at 892 ms.

图 3c所示非平稳地震记录的前200 ms划分为一个时窗,然后采用双谱重构方法提取该时窗内的子波相位谱并将其作为(14)式中时变子波相位谱的初值φ(t, ω).将φ(t, ω)以及提取的不同时刻子波的振幅谱代入(14)式可以提取不同时刻的时变子波相位谱,时变子波相位谱提取结果及提取误差分别如图 67所示.从中可以看出,提取的时变子波相位谱仅在少数频率点处存在误差,整体吻合度较高.

(37)

图 6 时变子波相位谱提取结果(111 ms, 381 ms, 619 ms, 892 ms) (a)理论子波的相位谱;(b)提取的子波相位谱. Fig. 6 Results of time-varying wavelet phase spectra extraction (111 ms, 381 ms, 619 ms, 892 ms) (a) Theoretical wavelet phase spectra; (b) Extracted wavelet phase spectra.
图 7 时变子波相位谱提取误差 (a) 111 ms处提取的子波相位谱与理论值的误差;(b) 381 ms处提取的子波相位谱与理论值的误差;(c) 619 ms处提取的子波相位谱与理论值的误差;(d) 892 ms处提取的子波相位谱与理论值的误差. Fig. 7 Errors of time-varying wavelet phase spectra extraction (a) Errors between the extracted and theoretical phase spectra at 111 ms; (b) Errors between the extracted and theoretical phase spectra at 381 ms; (c) Errors between the extracted and theoretical phase spectra at 619 ms; (d) Errors between the extracted and theoretical phase spectra at 892 ms.

n表示数据长度,observedt表示观测值,theoreticalt表示理论值.

将分别提取的时变子波振幅谱和相位谱按照(34)式逐点合并,实现时变子波的提取,将提取结果分别与理论子波和传统的分段方法提取的子波进行对比如图 8所示,提取误差如图 9所示.从中可以看出,相比于传统的分段提取方法,利用本文方法提取的时变子波与理论子波的误差更小,准确度更高.利用提取的时变子波对地震记录进行反褶积处理,反演反射系数,反演结果和反演误差分别如图 10图 11所示.利用(37)式所示的均方误差量化评价反演误差,传统方法的反演均方误差为0.0604,本文方法的反演均方误差为0.0293,结果表明相比于分段方法,利用本文方法提取的子波反演得到的反射系数更接近理论值.

图 8 时变子波提取结果(111 ms, 381 ms, 619 ms, 892 ms) (a)理论子波;(b)传统的分段方法提取的时变子波;(c)本文方法提取的时变子波. Fig. 8 Time-varying wavelet extraction results (111 ms, 381 ms, 619 ms, 892 ms) (a) Theoretical wavelets; (b) Extracted wavelets using the conventional segmentation method; (c) Extracted wavelet using the proposed method.
图 9 时变子波提取误差(111 ms, 381 ms, 619 ms, 892 ms) (a)分段方法的提取误差;(b)本文方法的提取误差. Fig. 9 Errors of extracted time-varying wavelets (111 ms, 381 ms, 619 ms, 892 ms) (a) Extraction errors of the conventional segmentation method; (b) Extraction errors of the proposed method.
图 10 反射系数反演结果 (a)理论反射系数;(b)分段方法反演的反射系数;(c)本文方法反演的反射系数. Fig. 10 Results of reflection coefficient inversion (a) Theoretical reflection coefficient; (b) Inverted reflection coefficient from the conventional segmentation method; (c) Inverted reflection coefficient from the proposed method.
图 11 反射系数反演误差 (a)分段方法的反演误差;(b)本文方法的反演误差. Fig. 11 Errors of the reflection coefficient inversion (a) Inversion errors of the conventional segmentation method; (b) Inversion errors of the proposed method.

接下来,我们将上述理论进行推广,测试本文方法在Q分段变化时的有效性, t=1~500 ms,Q1=150,t=501~1000 ms,Q2=120.同样采用EMD方法提取不同时刻子波的振幅谱,采用双谱重构方法估计子波相位谱初值,然后利用(14)式提取时变子波振幅谱和相位谱,将分别提取的振幅谱和相位谱按照(34)式逐点合成,实现时变子波的提取,提取结果和提取误差如图 1213所示.利用提取的子波对非平稳地震记录进行反褶积,反演反射系数,反演结果和反演误差如图 1415所示.利用(37)式所示的均方误差量化评价反演误差,传统方法的反演均方误差为0.0422,本文方法的反演均方误差为0.0236.结果表明,利用本文方法在Q不变以及Q分段变化的情况下均是有效可行的,且子波提取精度和反射系数反演精度均优于传统的分段提取方法.

图 12 Q1=150, Q2=120条件下的时变子波提取结果 (111 ms, 381 ms, 619 ms, 892 ms)
(a)理论子波;(b)传统的分段方法提取的时变子波;(c)本文方法提取的时变子波.
Fig. 12 Time-varying wavelet extraction results when Q1=150 and Q2=120 (111 ms, 381 ms, 619 ms, 892 ms)
(a) Theoretical wavelets; (b) Extracted wavelets when using the conventional segmentation method; (c) Extracted wavelets when using the proposed method.
图 13 Q1=150, Q2=120条件下的时变子波提取误差 (111 ms, 381 ms, 619 ms, 892 ms)
(a)分段方法的提取误差;(b)本文方法的提取误差.
Fig. 13 Errors of the extracted time-varying wavelets when Q1=150 and Q2=120 (111 ms, 381 ms, 619 ms, 892 ms) (a) Extraction errors of the conventional segmentation method; (b) Extraction errors of the proposed method.
图 14 Q1=150, Q2=120条件下的反射系数反演结果 (a)理论反射系数;(b)分段方法反演的反射系数;(c)本文方法反演的反射系数. Fig. 14 Results of the reflection coefficient inversion when Q1=150 and Q2=120 (a) Theoretical reflection coefficient; (b) Inverted reflection coefficient from the conventional segmentation method; (c) Inverted reflection coefficient from the proposed method.
图 15 Q1=150, Q2=120条件下的反射系数反演误差 (a)分段方法的反演误差;(b)本文方法的反演误差. Fig. 15 Errors of the reflection coefficient inversion when Q1=150 and Q2=120 (a) Inversion errors of the conventional segmentation method; (b) Inversion errors of the proposed method.
3.2 抗噪性能测试

为了验证本文方法的抗噪性能,向非平稳地震记录中分别加入不同强度的高斯白噪声和高斯色噪声,使得信噪比分别为10 dB、7 dB和3 dB,然后利用本文方法提取时变地震子波.图 16表示加入不同强度高斯白噪声的地震记录.图 1718分别为不同强度高斯白噪声干扰下的时变子波提取结果及其提取误差.图 19表示加入不同强度高斯色噪声的地震记录.图 2021分别为不同强度高斯色噪声干扰下的时变子波提取结果及其提取误差.

图 16 加入不同高斯白噪声的地震记录 (a) SNR=10 dB;(b) SNR=7 dB;(c) SNR=3 dB. Fig. 16 Seismograms with Gaussian white noise
图 17 高斯白噪声干扰下的时变子波提取结果(111 ms, 381 ms, 619 ms, 892 ms) (a)理论子波;(b) SNR=10 dB的子波提取结果;(c) SNR=7 dB的子波提取结果;(d) SNR=3 dB的子波提取结果. Fig. 17 Time-varying wavelet extraction results with different intensities of Gaussian white noise (111 ms, 381 ms, 619 ms, 892 ms) (a) Theoretical wavelets; (b) Extracted wavelets when SNR=10 dB; (c) Extracted wavelets when SNR=7 dB; (d) Extracted wavelets when SNR=3 dB.
图 18 高斯白噪声干扰下的时变子波提取误差(111 ms, 381 ms, 619 ms, 892 ms) (a) SNR=10 dB的子波提取误差;(b) SNR=7 dB的子波提取误差;(c) SNR=3 dB的子波提取误差. Fig. 18 Errors of the extracted time-varying wavelets with different intensities of Gaussian white noise (111 ms, 381 ms, 619 ms, 892 ms) (a) Extraction errors when SNR=10 dB; (b) Extraction errors when SNR=7 dB; (c) Extraction errors when SNR=3 dB.
图 19 加入不同高斯色噪声的地震记录 (a) SNR=10 dB;(b) SNR=7 dB;(c) SNR=3 dB. Fig. 19 Seismograms with Gaussian colored noise
图 20 高斯色噪声干扰下的时变子波提取结果(111 ms, 381 ms, 619 ms, 892 ms) (a)理论子波;(b) SNR=10 dB的子波提取结果;(c) SNR=7 dB的子波提取结果;(d) SNR=3 dB的子波提取结果. Fig. 20 Time-varying wavelet extraction results with different intensities of Gaussian colored noise (111 ms, 381 ms, 619 ms, 892 ms) (a) Theoretical wavelets; (b) Extracted wavelets when SNR=10 dB; (c) Extracted wavelets when SNR=7 dB; (d) Extracted wavelets when SNR=3 dB.
图 21 高斯色噪声干扰下的时变子波提取误差(111 ms, 381 ms, 619 ms, 892 ms) (a) SNR=10 dB的子波提取误差;(b) SNR=7 dB的子波提取误差;(c) SNR=3 dB的子波提取误差. Fig. 21 Errors of the extracted time-varying wavelets with different intensities of Gaussian colored noise (111 ms, 381 ms, 619 ms, 892 ms) (a) Extraction errors when SNR=10 dB; (b) Extraction errors when SNR=7 dB; (c) Extraction errors when SNR=3 dB.

上述实验结果表明,在不同强度的高斯白噪声和色噪声的干扰下,利用本文方法提取的时变子波依然是准确的,验证了该方法具有良好的抗噪性能.

4 实际资料处理

本文针对胜利油田某区块叠后地震剖面,利用EMD和子波振幅谱与相位谱的关系提取时变子波.以第81道为例,子波提取结果如图 22所示.从中可以看出,提取的时变子波符合地层吸收作用引起的动态衰减特性.为进一步验证提取结果的准确性和有效性,利用提取的子波对实际地震资料进行多道反褶积处理,处理结果如图 23所示,反褶积前后地震资料的振幅谱如图 24所示,结果表明,相比于传统的分段方法,利用本文方法进行反褶积处理后,地震资料分辨率更高,频带更宽,振幅能量更强,从而可以得出结论,利用EMD和子波振幅谱与相位谱的关系提取时变子波对提高地震资料分辨率具有重要的意义.

图 22 实际资料中提取时变子波 (a) 1.2 s处提取的子波;(b) 1.5 s处提取的子波;(c) 1.8 s处提取的子波. Fig. 22 Time-varying wavelet extraction in actual seismic data (a) Extracted wavelet at 1.2 s; (b) Extracted wavelet at 1.5 s; (c) Extracted wavelet at 1.8 s.
图 23 实际地震资料的时变子波反褶积处理结果 (a)原始叠后地震资料;(b)利用分段方法进行反褶积后的地震资料;(c)利用本文方法进行反褶积后的地震资料. Fig. 23 Time-varying wavelet deconvolution results of actual seismic data (a) Original pro-stack seismic data; (b) Seismic data after deconvolution with the conventional segmentation method; (c) Seismic data after deconvolution with the proposed method.
图 24 实际资料反褶积前后的振幅谱对比 (a)原始地震资料的振幅谱;(b)分段反褶积后地震资料的振幅谱;(c)本文方法反褶积后地震资料的振幅谱. Fig. 24 Comparison of the amplitude spectra before and after time-varying wavelet deconvolution (a) Amplitude spectrum of the original seismic data; (b) Amplitude spectrum of the seismic data after deconvolution with the conventional segmentation method; (c) Amplitude spectrum of the seismic data after deconvolution with the proposed method.

为了进一步针对实际资料处理评估本文方法的可信度和有效性,我们将实际资料的反褶积结果与测井得到的反射系数序列进行对比,结果如图 25所示,观察图中虚线框部分可以看出:分段方法提取的子波准确度不高导致反演的反射系数存在极性相反现象,且难以分辨相邻反射系数;而利用本文方法提取的时变子波精度更高,使得反射系数反演的效果更佳.利用(37)式量化评价反演误差,传统方法反演结果的均方误差为0.3940,本文方法反演结果的均方误差为0.1327.综合比较上述定性和定量分析结果可以看出,相比于传统的分段提取方法,利用本文方法进行反射系数反演的结果更接近测井反射系数,精度更高,进一步表明本文方法在实际资料处理中具有重要的应用价值.

图 25 实际资料的反射系数反演结果 (a)测井反射系数;(b)利用分段方法反演的反射系数;(c)利用本文方法反演的反射系数. Fig. 25 Reflection coefficient sequence inversion results of actual seismic data (a) Reflection coefficient sequence from the well log; (b) Inverted reflection coefficient sequence using the conventional segmentation method; (c) Inverted reflection coefficient sequence using the proposed method.
5 结论

本文深入分析了黏弹介质中地震波的传播规律,从地震波波动方程中总结了时变子波振幅谱和相位谱之间的关系并将其应用于时变子波提取中,针对现有方法中存在的问题,引入时变信号处理中的EMD方法,提出了一种利用EMD和子波振幅谱与相位谱关系的时变子波提取方法,通过理论分析、数值仿真实验和实际资料处理可以得到如下结论:

(1) 准确估计子波相位谱初值和不同时刻的振幅谱即可利用子波振幅谱与相位谱的关系提取时变子波的振幅谱和相位谱,从而实现时变地震子波的准确提取.

(2) 本文方法能够根据子波对数振幅谱光滑连续而反射系数对数振幅谱振荡剧烈的特点,从非平稳地震记录的对数振幅谱中提取子波的振幅信息,突破了传统的谱模拟方法对子波振幅谱形态的限制.同时利用子波振幅谱和相位谱的关系提取时变子波相位谱,突破了分段平稳假设的局限性.因此相比于传统方法,利用本文方法提取的时变子波精度更高,对提高地震资料分辨率具有重要的意义.

(3) 本文方法在Q不变、Q分段变化的情况下均是有效可行的,且在不同强度高斯白噪声和高斯色噪声的干扰下,利用本文方法提取的时变子波均是准确的,表明该方法具有一定的实际应用价值.

References
Aki K, Richards P G. 1980. Quantitative Seismology. San Francisco:W. H. Freeman: 5-43.
Boudraa A O, Cexus J C. 2007. EMD-based signal filtering. IEEE Transactions on Instrumentation and Measurement, 56(6): 2196-2202. DOI:10.1109/TIM.2007.907967
Braga I L S, Moraes F S. 2013. High-resolution gathers by inverse Q filtering in the wavelet domain. Geophysics, 78(2): V53-V61. DOI:10.1190/geo2011-0508.1
Chen X H, He Z H, Huang D J, et al. 2009. Low frequency shadow detection of gas reservoirs in time-frequency domain. Chinese Journal of Geophysics (in Chinese), 52(1): 215-221.
Clarke G K C. 1968. Time-varying deconvolution filters. Geophysics, 33(6): 936-944. DOI:10.1190/1.1439987
Dai Y S, Wang R R, Li C, et al. 2016. A time-varying wavelet extraction using local similarity. Geophysics, 81(1): V55-V68.
Economou N, Vafidis A. 2010. Spectral balancing GPR data using time-variant bandwidth in the t-f domain. Geophysics, 75(3): J19-J27. DOI:10.1190/1.3374464
Economou N, Vafidis A. 2012. GPR data time varying deconvolution by kurtosis maximization. Journal of Applied Geophysics, 81: 117-121. DOI:10.1016/j.jappgeo.2011.09.004
Gan Y, Sui L F, Wu J F, et al. 2014. An EMD threshold de-noising method for inertial sensors. Measurement, 49: 34-41. DOI:10.1016/j.measurement.2013.11.030
Han J J, van der Baan M. 2013. Empirical mode decomposition for seismic time-frequency analysis. Geophysics, 78(2): O9-O19. DOI:10.1190/geo2012-0199.1
Huang N E, Shen Z, Long S R, et al. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society of London. Series A:Mathematical, Physical and Engineering Sciences, 454(1971): 903-995. DOI:10.1098/rspa.1998.0193
Maly A, Mowlaee P. 2016. On the importance of harmonic phase modification for improved speech signal reconstruction.//2016 IEEE International Conference on Acoustics, Speech and Signal Processing. Shanghai: IEEE, 584-588.
Margrave G F. 1998. Theory of nonstationary linear filtering in the Fourier domain with application to time-variant filtering. Geophysics, 63(1): 244-259. DOI:10.1190/1.1444318
Margrave G F, Lamoureux M P, Henley D C. 2011. Gabor deconvolution:Estimating reflectivity by nonstationary deconvolution of seismic data. Geophysics, 76(3): W15-W30. DOI:10.1190/1.3560167
Matsuoka T, Ulrych T J. 1984. Phase estimation using the bispectrum. Proceedings of the IEEE, 72(10): 1403-1411. DOI:10.1109/PROC.1984.13027
Rickett J. 2006. Integrated estimation of interval-attenuation profiles. Geophysics, 71(4): A19-A23. DOI:10.1190/1.2209722
Rickett J. 2007. Estimating attenuation and the relative information content of amplitude and phase spectra. Geophysics, 72(1): R19-R27. DOI:10.1190/1.2399451
Rosa A L R, Ulrych T J. 1991. Processing via spectral modeling. Geophysics, 56(8): 1244-1251. DOI:10.1190/1.1443144
Sinha S, Routh P, Anno P. 2009. Instantaneous spectral attributes using scales in continuous-wavelet transform. Geophysics, 74(2): WA137-WA142. DOI:10.1190/1.3054145
Stockwell R G. 2007. A basis for efficient representation of the S-transform. Digital Signal Processing, 17(1): 371-393. DOI:10.1016/j.dsp.2006.04.006
van der Baan M. 2008. Time-varying wavelet estimation and deconvolution by kurtosis maximization. Geophysics, 73(2): V11-V18.
van der Baan M. 2012. Bandwidth enhancement:Inverse Q filtering or time-varying Wiener deconvolution?. Geophysics, 77(4): V133-V142. DOI:10.1190/geo2011-0500.1
Wang L L, Gao J H, Zhang M. 2010. A method for absorption compensation based on adaptive molecular decomposition. Applied Geophysics, 7(1): 74-87. DOI:10.1007/s11770-010-0010-5
Wang L L, Gao J H, Zhao W, et al. 2012. Enhancing resolution of nonstationary seismic data by molecular-Gabor transform. Geophysics, 78(1): V31-V41.
Wang S D, Song H J, Yang D F. 2014. Seismic attenuation compensation by Bayesian inversion. Journal of Applied Geophysics, 111: 356-363. DOI:10.1016/j.jappgeo.2014.10.012
Wang S J, Fang Z Y, Shi W Y, et al. 2015. Marine seismic wavelet zero-phasing technology and its application. Geophysical Prospecting for Petroleum (in Chinese), 54(5): 551-559.
Wang T. 2010. Research on EMD algorithm and its application in signal denoising[Ph. D. thesis] (in Chinese). Harbin: Harbin Engineering University.
Wang X H, Qin Y L, Huang Z P, et al. 1998. Time-variant wavelet deconvolution based on spectrum analysis. Geophysical Prospecting for Petroleum (in Chinese), 37(S1): 109-112.
Wang Y H. 2002. A stable and efficient approach of inverse Q filtering. Geophysics, 67(2): 657-663. DOI:10.1190/1.1468627
Wang Y H. 2006. Inverse Q-filter for seismic resolution enhancement. Geophysics, 71(3): V51-V60. DOI:10.1190/1.2192912
Wu D Z. 1998. Analysis of Signals and Linear Systems (in Chinese). 3rd ed. Beijing: Higher Education Press.
Wu Z H, Huang N E. 2009. Ensemble empirical mode decomposition:a noise-assisted data analysis method. Advances in Adaptive Data Analysis, 1(1): 1-44.
Zhou H L, Wang J, Wang M C, et al. 2014. Amplitude spectrum compensation and phase spectrum correction of seismic data based on the generalized S transform. Applied Geophysics, 11(4): 468-478. DOI:10.1007/s11770-014-0456-y
陈学华, 贺振华, 黄德济, 等. 2009. 时频域油气储层低频阴影检测. 地球物理学报, 52(1): 215-221.
王守君, 方中于, 史文英, 等. 2015. 海洋地震资料子波零相位化技术研究与应用. 石油物探, 54(5): 551-559. DOI:10.3969/j.issn.1000-1441.2015.05.008
王晓华, 秦义龙, 黄真萍, 等. 1998. 基于频谱分析的时变子波反褶积. 石油物探, 37(S1): 109-112.
王婷. 2010. EMD算法研究及其在信号去噪中的应用[博士论文].哈尔滨: 哈尔滨工程大学. http://cdmd.cnki.com.cn/Article/CDMD-10217-1011020271.htm
吴大正. 1998. 信号与线性系统分析. 3版. 北京: 高等教育出版社.