基于三参数小波的瞬时频率分析方法

基于三参数小波的瞬时频率分析方法

地震信号的瞬时性特征包括瞬时振幅、瞬时频率、瞬时相位和瞬时振幅。瞬时属性可以用于分析地层岩性变化、地下构造及反演地层中的岩性参数。分析瞬时属性的方法有时间-频率域方法和复信号方法等。时间-频率域方法有短时Fourier变换、小波变换、S变换及广义S变换、Wigner分布、匹配追踪等。复信号计算方法也有多种,但地震信号处理领域常用的是Hilbert变换法。该方法首先计算待分析地震信号的Hilbert变换,并将其结果作为虚部;然后将待分析的地震信号作为实部构成复信号,并利用该复信号计算瞬时频率等瞬时属性。我们将这类方法称为Hilbert变换属性提取法。Hilbert变换属性提取法已经在地层结构和地下岩性变化研究及参数反演中得到了广泛应用,并取得了很好的效果。但该方法也存在着局限性:计算瞬时频率等属性时对噪声敏感,难以进行低信噪比资料分析;不能提供多分辨率瞬时属性,即无法对地下地质体进行最佳分辨率解释等。为了提高地震资料解释的可靠性,Grossman等提出了基于小波变换的地震资料多分辨率解释方法。高静怀等讨论了小波变换用于地震资料多分辨率分析时小波函数的选择问题,并构造出了适合多分辨率解释的小波——匹配地震子波(或待分析有效信号)的小波。为了突破Hilbert变换属性提取法存在的局限,高静怀等研究了待分析信号的解析小波变换与Hilbert变换之间的关系,给出了在小波变换域计算地震信号对应的解析信号,进而进行瞬时属性分析的方法。该方法在采用匹配地震子波的小波或Morlet小波进行瞬时频率分析时,如果基本小波中的调制频率较小,则小波函数的时间局部化性质变差,从而影响分析效果。为此,高静怀等借鉴Harrop等人的研究成果,提出了三参数小波。本文基于小波变换域瞬时属性提取的原理,分析了Morlet小波、匹配地震子波的小波及三参数小波用于计算瞬时频率的抗噪性能及分辨率。1在小波转换区中提取瞬时属性原则的原则1.1有效信号能量分布的确定若小波函数g(t)和它的Fourier变换g(ω)满足对于一个给定信号s(t)∈L2(R,dt),s(t)关于小波函数g(t)的小波变换定义为式中:t,b∈R;a>0;g(t)是g(t)的复共轭;L1和L2分别代表绝对可积和平方可积函数空间。小波变换的反变换公式为式中:s(t)为s(t)对应的复信号;Ω为待分析信号在时间-尺度域的能量分布区域。令gR(t)表示g(t)的实部,R(ω)为gR(t)的傅里叶变换。gR(t)和R(ω)满足对于任意的能量有限信号s(t),即s(t)∈L2(R,dt),有式中:H[s(t)]代表s(t)的Hilbert变换。(8)式表明,当基本小波满足(6)、(7a)和(7b)式时,由(5)式得到的复信号就是s(t)对应的解析信号。图la为待分析的含噪信号,图lb为其小波变换结果(取模)。可见,在小波变换域中,有效信号的能量分布比较集中,随机噪声等干扰能量分布比较分散,噪声在有效信号能量分布区仅有一小部分。因此,在有效信号能量分布区内计算地震信号对应的复信号,可以使得噪声对计算结果的影响较小。对于有效信号能量分布范围的确定,可以通过有效信号与噪声在统计特性上的差别来进行。图lc是分离出的有效信号能量分布区域,图1d是在有效信号能量分布区利用(4)式重构的有效信号。比较图la与图1d可见,待分析信号的噪声得到了很好的压制。1.2瞬时属性的确定将待分析信号s(t)对应的复信号记为(t),(t)的虚部记为s1(t)。s(t)对应的瞬时属性为式中:e(t),θ(t)和f(t)分别代表s(t)的瞬时幅度、瞬时相位和瞬时频率。本文在数值计算中采用阻尼瞬时频率,即式中:=max[e2(t)];ε为阻尼因子,0<ε<1,当信噪比较高的时候,应该取较小的ε值。2小波函数选择2.1小波的局部化性质匹配地震子波的小波定义为:式中:A,σ,τ和β分别为基本小波的振幅、调制频率、能量衰减因子和能量延迟因子。为了书写简便,用矢量Г=(A,σ,τ,β)表示参变量A,σ,τ,β的一组取值,则(13)式可表示为式中:当采用Morlet小波和匹配地震子波的小波分析振幅慢变的信号时,σ2/(8τ)取足够大的值即可。在这种情况下,(14)式中的修正项R(t;Γ)可以略去。然而,如果用这类小波分析振幅快速变化的信号,为了得到较高的时间分辨率和频率分辨率,就要求σ取较小的值。在这种情况下,修正项R(t;Γ)不能忽略。虽然匹配地震子波的小波在σ值较小时仍满足容许条件(Morlet小波与此类似),但其时间局部化性质不能令人满意,即|Ψ(t)|由单峰变为多峰。使用具有多峰的分析小波对待分析信号作小波变换,会导致信号关于多个不同位置局部化,在瞬时频率和瞬时振幅分析中产生假象。2.2地震波和改进的molret小波的时间域特征分析三参数小波的定义为式中:σ为三参数小波的调制频率;σ,τ,β∈R,且σ,τ≥0。为了书写方便,用矢量Λ=(σ,τ,β)表示参数σ,τ,β的集合,则Ψ(t;σ,τ,β)可记为Ψ(t;Λ),其它量类似,(16a)式可简写为式中:当τ=0.5,β=0,Λ1=(σ,0.5,0)时将(20)式～(22)式代入(16b)式得此时,三参数小波相当于改进的Morlet小波。当σ>5.33时,令A2=(σ>5.33,0.5,0),p(Λ2)和q(Λ2)都近似等于,于是此即Morlet小波。因此,三参数小波包含了Morlet小波和改进的Morlet小波。当σ取值很小,例如σ=1时,基于三参数小波的瞬时属性分析具有较高的时间分辨率。图2给出了Λ1=(1,0.5,0)时三参数小波的时间域波形及其傅里叶振幅谱,可见三参数小波在时间域具有单峰,而且其傅里叶振幅谱的负频率成分幅度很小。根据解析信号的定义,三参数小波可近似地看作解析信号。以匹配地震子波的小波、Morlet小波等为基本小波,利用(5)式得到的复信号,在何种条件下为解析信号?文献对此已有详细讨论。定性地说,以三参数小波为基本小波,利用(5)式计算复信号时,σ越大复信号越接近于解析信号,σ较小时复信号可近似看成解析信号。在有些情况下,由近似解析的复信号得到的瞬时频率等属性,更能刻画信号中所含的特殊信息。3基于三个参数的小波时间频率计算3.1复信号属性的确定图3为主频fc为50Hz的Ricker子波,子波周期为20ms,表达式为为了观察三参数小波中参数对瞬时频率等属性分析结果的影响,分别以Λ1=(1,0.5,0),Λ2=(2,0.5,0),Λ3=(3,0.5,0)和Λ6=(6,0.5,0)等4组参数构造基本小波,计算Ricker子波的瞬时频率,结果如图4所示。由图4可以清晰地看出,当参数为Λ1=(1,0.5,0)时,瞬时频率在时间轴中间有一个尖峰且正好对应Ricker子波的尖峰位置;随着σ的增大,此尖脉冲逐渐展宽,到Λ6=(6,0.5,0)(Morlet小波)的时候,瞬时频率曲线的形状与Hilbert变换属性提取方法的计算结果几乎相同,不出现尖峰。上述结果启发我们,采用与实信号对应的近似解析的复信号形式上定义瞬时频率等属性,比用解析信号定义的同类属性时间分辨率更高,更有益于薄互层结构的地震响应分析。选取50Hz和30Hz两个Ricker子波并加入随机噪声,峰-峰信噪比为25:1,如图5所示。当参数向量分别为Λ1=(1,0.5,0)和Λ6=(6,0.5,0),也就是所选取的小波分别为时间分辨率较高的三参数小波和传统的Morlet小波时,计算出的瞬时频率分别如图6所示。由图6可见,采用Λ1=(1,0.5,0)得到的瞬时频率噪声幅度相对较低,且存在尖的峰值,有利于确定反射界面的位置。若采用基于Hilbert变换的瞬时频率估计方法,对于图5给出的含噪信号,则无法得到有效的瞬时频率。3.2不同参数构造的差互层模型图7为薄互层模型的反射系数序列及其与50Hz的Ricker子波合成的地震记录。反射系数幅度相同(0.5),正负相间,时间为(200,201,203,206,210,215,221,228,236,245,255,264,272,279,285,290,294,297,299,300)ms。分别以Λ1=(1,0.5,0)和Λ6=(6,0.5,0)为参数构造基本小波,计算出薄互层模型的瞬时频率,如图8所示。由图8b可见,瞬时频率的极小值点准确地确定了(200,215,221,228,236,245,255,264,272,279,285,300)ms处的反射系数位置;瞬时频率的整体变化趋势与图8a一致。也就是说,用三参数小波计算出的瞬时频率属性,既可以反映薄互层的变化趋势,也可以刻画薄互层的内部结构。3.3波变换域时效分析方法对比图9是对大庆油田某地区三维数据体采用Hilbert变换属性提取方法和小波变换域瞬时频率分析方法计算的瞬时频率切片对比图。通过比较可以看出,采用小波变换域方法计算的瞬时频率属性更清楚地反映了火山岩主体的展布范围及储层内部的结构变化。图中火山岩主体上的xs301井和xs3井是产气井。4小波的选择应用1)在三参数小波中,参数取不同的值可得到不同的小波。用不同的小波得到的瞬时频率对地震信号的刻画