基于三參數(shù)小波的瞬時頻率分析方法

基于三參數(shù)小波的瞬時頻率分析方法

地震信號的瞬時性特征包括瞬時振幅、瞬時頻率、瞬時相位和瞬時振幅。瞬時屬性可以用于分析地層巖性變化、地下構(gòu)造及反演地層中的巖性參數(shù)。分析瞬時屬性的方法有時間-頻率域方法和復(fù)信號方法等。時間-頻率域方法有短時Fourier變換、小波變換、S變換及廣義S變換、Wigner分布、匹配追蹤等。復(fù)信號計算方法也有多種,但地震信號處理領(lǐng)域常用的是Hilbert變換法。該方法首先計算待分析地震信號的Hilbert變換,并將其結(jié)果作為虛部;然后將待分析的地震信號作為實部構(gòu)成復(fù)信號,并利用該復(fù)信號計算瞬時頻率等瞬時屬性。我們將這類方法稱為Hilbert變換屬性提取法。Hilbert變換屬性提取法已經(jīng)在地層結(jié)構(gòu)和地下巖性變化研究及參數(shù)反演中得到了廣泛應(yīng)用,并取得了很好的效果。但該方法也存在著局限性:計算瞬時頻率等屬性時對噪聲敏感,難以進(jìn)行低信噪比資料分析;不能提供多分辨率瞬時屬性,即無法對地下地質(zhì)體進(jìn)行最佳分辨率解釋等。為了提高地震資料解釋的可靠性,Grossman等提出了基于小波變換的地震資料多分辨率解釋方法。高靜懷等討論了小波變換用于地震資料多分辨率分析時小波函數(shù)的選擇問題,并構(gòu)造出了適合多分辨率解釋的小波——匹配地震子波(或待分析有效信號)的小波。為了突破Hilbert變換屬性提取法存在的局限,高靜懷等研究了待分析信號的解析小波變換與Hilbert變換之間的關(guān)系,給出了在小波變換域計算地震信號對應(yīng)的解析信號,進(jìn)而進(jìn)行瞬時屬性分析的方法。該方法在采用匹配地震子波的小波或Morlet小波進(jìn)行瞬時頻率分析時,如果基本小波中的調(diào)制頻率較小,則小波函數(shù)的時間局部化性質(zhì)變差,從而影響分析效果。為此,高靜懷等借鑒Harrop等人的研究成果,提出了三參數(shù)小波。本文基于小波變換域瞬時屬性提取的原理,分析了Morlet小波、匹配地震子波的小波及三參數(shù)小波用于計算瞬時頻率的抗噪性能及分辨率。1在小波轉(zhuǎn)換區(qū)中提取瞬時屬性原則的原則1.1有效信號能量分布的確定若小波函數(shù)g(t)和它的Fourier變換g(ω)滿足對于一個給定信號s(t)∈L2(R,dt),s(t)關(guān)于小波函數(shù)g(t)的小波變換定義為式中:t,b∈R;a>0;g(t)是g(t)的復(fù)共軛;L1和L2分別代表絕對可積和平方可積函數(shù)空間。小波變換的反變換公式為式中:s(t)為s(t)對應(yīng)的復(fù)信號;Ω為待分析信號在時間-尺度域的能量分布區(qū)域。令gR(t)表示g(t)的實部,R(ω)為gR(t)的傅里葉變換。gR(t)和R(ω)滿足對于任意的能量有限信號s(t),即s(t)∈L2(R,dt),有式中:H[s(t)]代表s(t)的Hilbert變換。(8)式表明,當(dāng)基本小波滿足(6)、(7a)和(7b)式時,由(5)式得到的復(fù)信號就是s(t)對應(yīng)的解析信號。圖la為待分析的含噪信號,圖lb為其小波變換結(jié)果(取模)。可見,在小波變換域中,有效信號的能量分布比較集中,隨機(jī)噪聲等干擾能量分布比較分散,噪聲在有效信號能量分布區(qū)僅有一小部分。因此,在有效信號能量分布區(qū)內(nèi)計算地震信號對應(yīng)的復(fù)信號,可以使得噪聲對計算結(jié)果的影響較小。對于有效信號能量分布范圍的確定,可以通過有效信號與噪聲在統(tǒng)計特性上的差別來進(jìn)行。圖lc是分離出的有效信號能量分布區(qū)域,圖1d是在有效信號能量分布區(qū)利用(4)式重構(gòu)的有效信號。比較圖la與圖1d可見,待分析信號的噪聲得到了很好的壓制。1.2瞬時屬性的確定將待分析信號s(t)對應(yīng)的復(fù)信號記為(t),(t)的虛部記為s1(t)。s(t)對應(yīng)的瞬時屬性為式中:e(t),θ(t)和f(t)分別代表s(t)的瞬時幅度、瞬時相位和瞬時頻率。本文在數(shù)值計算中采用阻尼瞬時頻率,即式中:=max[e2(t)];ε為阻尼因子,0<ε<1,當(dāng)信噪比較高的時候,應(yīng)該取較小的ε值。2小波函數(shù)選擇2.1小波的局部化性質(zhì)匹配地震子波的小波定義為:式中:A,σ,τ和β分別為基本小波的振幅、調(diào)制頻率、能量衰減因子和能量延遲因子。為了書寫簡便,用矢量Г=(A,σ,τ,β)表示參變量A,σ,τ,β的一組取值,則(13)式可表示為式中:當(dāng)采用Morlet小波和匹配地震子波的小波分析振幅慢變的信號時,σ2/(8τ)取足夠大的值即可。在這種情況下,(14)式中的修正項R(t;Γ)可以略去。然而,如果用這類小波分析振幅快速變化的信號,為了得到較高的時間分辨率和頻率分辨率,就要求σ取較小的值。在這種情況下,修正項R(t;Γ)不能忽略。雖然匹配地震子波的小波在σ值較小時仍滿足容許條件(Morlet小波與此類似),但其時間局部化性質(zhì)不能令人滿意,即|Ψ(t)|由單峰變?yōu)槎喾?。使用具有多峰的分析小波對待分析信號作小波變換,會導(dǎo)致信號關(guān)于多個不同位置局部化,在瞬時頻率和瞬時振幅分析中產(chǎn)生假象。2.2地震波和改進(jìn)的molret小波的時間域特征分析三參數(shù)小波的定義為式中:σ為三參數(shù)小波的調(diào)制頻率;σ,τ,β∈R,且σ,τ≥0。為了書寫方便,用矢量Λ=(σ,τ,β)表示參數(shù)σ,τ,β的集合,則Ψ(t;σ,τ,β)可記為Ψ(t;Λ),其它量類似,(16a)式可簡寫為式中:當(dāng)τ=0.5,β=0,Λ1=(σ,0.5,0)時將(20)式～(22)式代入(16b)式得此時,三參數(shù)小波相當(dāng)于改進(jìn)的Morlet小波。當(dāng)σ>5.33時,令A(yù)2=(σ>5.33,0.5,0),p(Λ2)和q(Λ2)都近似等于,于是此即Morlet小波。因此,三參數(shù)小波包含了Morlet小波和改進(jìn)的Morlet小波。當(dāng)σ取值很小,例如σ=1時,基于三參數(shù)小波的瞬時屬性分析具有較高的時間分辨率。圖2給出了Λ1=(1,0.5,0)時三參數(shù)小波的時間域波形及其傅里葉振幅譜,可見三參數(shù)小波在時間域具有單峰,而且其傅里葉振幅譜的負(fù)頻率成分幅度很小。根據(jù)解析信號的定義,三參數(shù)小波可近似地看作解析信號。以匹配地震子波的小波、Morlet小波等為基本小波,利用(5)式得到的復(fù)信號,在何種條件下為解析信號?文獻(xiàn)對此已有詳細(xì)討論。定性地說,以三參數(shù)小波為基本小波,利用(5)式計算復(fù)信號時,σ越大復(fù)信號越接近于解析信號,σ較小時復(fù)信號可近似看成解析信號。在有些情況下,由近似解析的復(fù)信號得到的瞬時頻率等屬性,更能刻畫信號中所含的特殊信息。3基于三個參數(shù)的小波時間頻率計算3.1復(fù)信號屬性的確定圖3為主頻fc為50Hz的Ricker子波,子波周期為20ms,表達(dá)式為為了觀察三參數(shù)小波中參數(shù)對瞬時頻率等屬性分析結(jié)果的影響,分別以Λ1=(1,0.5,0),Λ2=(2,0.5,0),Λ3=(3,0.5,0)和Λ6=(6,0.5,0)等4組參數(shù)構(gòu)造基本小波,計算Ricker子波的瞬時頻率,結(jié)果如圖4所示。由圖4可以清晰地看出,當(dāng)參數(shù)為Λ1=(1,0.5,0)時,瞬時頻率在時間軸中間有一個尖峰且正好對應(yīng)Ricker子波的尖峰位置;隨著σ的增大,此尖脈沖逐漸展寬,到Λ6=(6,0.5,0)(Morlet小波)的時候,瞬時頻率曲線的形狀與Hilbert變換屬性提取方法的計算結(jié)果幾乎相同,不出現(xiàn)尖峰。上述結(jié)果啟發(fā)我們,采用與實信號對應(yīng)的近似解析的復(fù)信號形式上定義瞬時頻率等屬性,比用解析信號定義的同類屬性時間分辨率更高,更有益于薄互層結(jié)構(gòu)的地震響應(yīng)分析。選取50Hz和30Hz兩個Ricker子波并加入隨機(jī)噪聲,峰-峰信噪比為25:1,如圖5所示。當(dāng)參數(shù)向量分別為Λ1=(1,0.5,0)和Λ6=(6,0.5,0),也就是所選取的小波分別為時間分辨率較高的三參數(shù)小波和傳統(tǒng)的Morlet小波時,計算出的瞬時頻率分別如圖6所示。由圖6可見,采用Λ1=(1,0.5,0)得到的瞬時頻率噪聲幅度相對較低,且存在尖的峰值,有利于確定反射界面的位置。若采用基于Hilbert變換的瞬時頻率估計方法,對于圖5給出的含噪信號,則無法得到有效的瞬時頻率。3.2不同參數(shù)構(gòu)造的差互層模型圖7為薄互層模型的反射系數(shù)序列及其與50Hz的Ricker子波合成的地震記錄。反射系數(shù)幅度相同(0.5),正負(fù)相間,時間為(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)為參數(shù)構(gòu)造基本小波,計算出薄互層模型的瞬時頻率,如圖8所示。由圖8b可見,瞬時頻率的極小值點準(zhǔn)確地確定了(200,215,221,228,236,245,255,264,272,279,285,300)ms處的反射系數(shù)位置;瞬時頻率的整體變化趨勢與圖8a一致。也就是說,用三參數(shù)小波計算出的瞬時頻率屬性,既可以反映薄互層的變化趨勢,也可以刻畫薄互層的內(nèi)部結(jié)構(gòu)。3.3波變換域時效分析方法對比圖9是對大慶油田某地區(qū)三維數(shù)據(jù)體采用Hilbert變換屬性提取方法和小波變換域瞬時頻率分析方法計算的瞬時頻率切片對比圖。通過比較可以看出,采用小波變換域方法計算的瞬時頻率屬性更清楚地反映了火山巖主體的展布范圍及儲層內(nèi)部的結(jié)構(gòu)變化。圖中火山巖主體上的xs301井和xs3井是產(chǎn)氣井。4小波的選擇應(yīng)用1)在三參數(shù)小波中,參數(shù)取不同的值可得到不同的小波。用不同的小波得到的瞬時頻率對地震信號的刻畫