本發(fā)明屬于地球物理勘探中的信號處理領(lǐng)域，具體涉及一種基于同步擠壓變換的地震加權(quán)平均瞬時頻率提取方法。
背景技術(shù)：
：地震屬性是指由地震數(shù)據(jù)經(jīng)過數(shù)學(xué)變換導(dǎo)出的有關(guān)地震波的幾何形態(tài)、運動學(xué)特征、動力學(xué)特征和統(tǒng)計學(xué)特征，這些特征可以用來預(yù)測地下巖性，描述油藏內(nèi)部的儲集特征等，在地震資料處理過程中發(fā)揮了很大的作用。20世紀(jì)70年代，地震屬性分析開始被引入地震解釋中。20世紀(jì)90年代以來，由于儲層描述和三維數(shù)據(jù)體解釋的需要，地震屬性分析技術(shù)迅速發(fā)展。但目前還沒有一個公認(rèn)的地震屬性分類，按照物理意義，地震屬性可以分為時間、振幅、頻率、相干、衰減等幾大類，按照屬性拾取方法，每類又可以分為瞬時屬性、單道時窗屬性、多道時窗屬性。其中瞬時屬性是在地震波到達(dá)的位置上拾取的屬性，包括瞬時振幅、瞬時相位、瞬時頻率、瞬時帶寬等等。近年來，地震屬性的數(shù)量突增，然而，瞬時屬性仍然是地震數(shù)據(jù)地質(zhì)解釋的支柱，也是目前商業(yè)處理軟件必備的技術(shù)模塊。在信號處理領(lǐng)域，以瞬時振幅、瞬時相位、瞬時頻率為代表的信號瞬時屬性的概念由來已久，在Gabor和Ville的工作之后，大量學(xué)者在相關(guān)領(lǐng)域做了研究，Boashash對這些工作做了綜述。Taner等人于1979年引入復(fù)地震道分析，提出通過復(fù)地震道求取瞬時屬性的方法，并給出了這些屬性的物理意義以及在地震解釋中的應(yīng)用。瞬時振幅與相鄰層的巖性變化及油氣聚集有關(guān)。瞬時相位反映界面的不連續(xù)性、斷層、不整合面和層序邊界等。瞬時頻率，作為最常用的瞬時屬性之一，可以有效刻畫地層的厚度和巖性變化，指示油氣的分布等。Zeng利用瞬時頻率異常來指示薄層。Gao等人利用瞬時頻率進(jìn)行地震資料Q值估計。本發(fā)明主要研究地震資料瞬時頻率的提取。瞬時頻率屬性是一種時變參數(shù)，定義為輸入非平穩(wěn)信號隨著時間變換的頻譜峰值位置。根據(jù)其原始定義，瞬時頻率提供了一種對信號在頻率域能量聚集性的度量方式。常用的提取地震瞬時頻率的方法可以分為以下三類：(1)Hilbert變換方法。對地震道做Hilbert變換，轉(zhuǎn)化為復(fù)地震道，其中實部為原地震道數(shù)據(jù)，虛部為其Hilbert變換。得到復(fù)地震道之后，就可以在每個采樣點計算瞬時頻率。在Taner的工作之后，Hilbert變換法被廣泛用于計算地震瞬時頻率，至今大部分商業(yè)軟件仍采用該方法。很多學(xué)者對該方法進(jìn)行了發(fā)展和改進(jìn)，Barnes于1996年提出二維復(fù)地震道分析的概念。Luo等人于2003年提出廣義Hilbert變換并給出其在地球物理方面的應(yīng)用。Barnes于2007年提出加權(quán)瞬時頻率的概念。Lu和Zhang將Hilbert變換方法進(jìn)行推廣，提出基于自適應(yīng)濾波器計算瞬時頻率的方法。(2)基于時頻分析的方法。時頻分析方法利用信號的時頻分布來求取瞬時頻率。Boashash等人提出基于時頻分析的自適應(yīng)瞬時頻率估計方法。Stankovic等人利用自適應(yīng)窗時頻分布計算瞬時頻率。高靜懷等人提出了在相空間計算地震瞬時頻率的方法。Marfurt等人提出基于窄帶譜分析的地震資料瞬時頻率提取方法。Huang等人提出基于經(jīng)驗?zāi)Ｊ椒纸?EMD)計算瞬時頻率的方法。Han等人將完全總體經(jīng)驗?zāi)Ｊ椒纸?CEEMD)方法用于提取地震資料的瞬時頻率。(3)基于反演的方法。Fomel等人提出利用反演的方法來獲得局部屬性，相比于瞬時頻率，局部頻率物理意義更加明確，且應(yīng)用效果明顯。Liu等人提出基于反演的方法，計算地震資料的瞬時頻率。技術(shù)實現(xiàn)要素：本發(fā)明的目的在于克服上述不足，提供一種基于同步擠壓變換的地震加權(quán)平均瞬時頻率提取方法，具有良好的抗噪性能和精度，將該方法應(yīng)用于實際資料，能夠更加清晰地刻畫儲層的特征。為了達(dá)到上述目的，本發(fā)明包括以下步驟：步驟一，根據(jù)實地震道進(jìn)行三參數(shù)小波變換；步驟二，計算重排準(zhǔn)則，通過該準(zhǔn)則進(jìn)行能量重排可得到更集中的時頻表示；步驟三，將能量重排，得到了擠壓后的時頻分布；步驟四，進(jìn)行有效信號能量分布空間的確定，得到有效信號對應(yīng)的變換系數(shù)；步驟五，進(jìn)行加權(quán)平均瞬時頻率提取，完成平均瞬時頻率。所述步驟一中，實地震道s(t)的小波變換定義為Ws(a,b)=a-1∫s(t)ψ‾(t-ba)dt=12π∫-∞∞S(ω)ψ^‾(aω)eiωbdω,---(1)]]>其中，為ψ(t)的復(fù)共軛，a為尺度因子，b為時移因子，是基本小波ψ(t)的Fourier變換，是s(t)的Fourier變換，基本小波平方可積，無直流分量，滿足以下容許性條件：Cψ=∫-∞∞|ψ(ω)|2|ω|dω<∞.---(2)]]>三參數(shù)小波的定義為：ψ(t;σ,τ,β)=e-τ(t-β)2{p(σ,τ,β)[cos(σt)-k(σ,τ,β)]+iq(σ,τ,β)sin(σt)},---(3)]]>其中，τ為能量衰減因子，β為能量延遲時間，σ為分析小波調(diào)制頻率，(σ,τ,β∈R且σ,τ＞0)，用向量Λ＝(σ,τ,β)記參數(shù)集合，則ψ(t；σ,τ,β)可記為ψ(t；Λ)，其他量類似，以向量的形式，上式可簡寫為：ψ(t;Λ)=e-τ(t-β)2{p(Λ)[cos(σt)-k(Λ)]+iq(Λ)sin(σt)},---(4)]]>其中，k(Λ)=e-σ24τ[cos(βσ)+iq(Λ)p(Λ)sin(βσ)],p(Λ)=(2τπ)14[4(e-σ22τ-e-3σ28τ)×cos2(βσ)+1-e-σ22τ]-12,q(Λ)=(2τπ)14[4(e-σ22τ-e-3σ28τ)×sin2(βσ)+1-e-σ22τ]-12.---(5)]]>對(4)式作Fourier變換，得到其頻域形式為ψ^(ω;Λ)=∫-∞∞ψ(t;Λ)e-iωtdt=πτp(Λ)+q(Λ)2e-iβ(ω-σ)-(ω-σ)24τ+πτp(Λ)-q(Λ)2e-iβ(ω+σ)-(ω+σ)24τ-πτp(Λ)k(Λ)e-iβω-ω24τ.---(6)]]>所述步驟二中，假設(shè)一個單頻余弦信號，即s(t)＝Acos(ω0t)，s(t)的Fourier變換為根據(jù)(1)式，信號s(t)關(guān)于三參數(shù)小波ψ(t；Λ)的小波變換結(jié)果為Ws(a,b)=12π∫-∞+∞s^(ω)ψ^‾(aω;Λ)eiωbdω=A2π[ψ^‾(aω0;Λ)eiω0b+ψ^‾(-aω0;Λ)e-iω0b].---(7)]]>因為三參數(shù)小波沒有負(fù)頻率分量，即ω＜0，且只考慮正尺度，則(7)式可以簡化為Ws(a,b)=A2πψ^‾(aω0;Λ)eiω0b.---(8)]]>從(8)式可以看到，假如三參數(shù)小波ψ(t；Λ)的峰值頻率為ξM，則小波變換的結(jié)果將在尺度a＝ξM/ω0處取到最大值，并以這個能量最大的尺度為中心形成一個尺度帶，造成能量的擴(kuò)散，為了得到更集中的時頻分布，計算重排準(zhǔn)則——瞬時頻率：ωs(a,b)=∂bWs(a,b)2πiWs(a,b),---(9)]]>其中，Ws(a,b)≠0，當(dāng)Ws(a,b)＝0時，信號在此處無能量，(9)式將無窮大，計算瞬時頻率ωs(a,b)沒有意義，所以計算時要求當(dāng)噪聲較小時，可以選擇閾值如果噪聲較大，需要抑制噪聲，選取自適應(yīng)閾值作為閾值其中ση是和噪聲水平相關(guān)，它是由小波變換系數(shù)的前nv個尺度來估計的，median是matlab中的庫函數(shù)，代表取中值，0.6745是針對高斯噪聲的正則化因子，假設(shè)有N個時刻，所以可以對N個時刻求得的ση取均值，能得到最優(yōu)的閾值。當(dāng)參考信號是單頻余弦信號時，其小波變換會形成一個尺度帶，但是可以由(9)計算其瞬時頻率為ωs(a,b)＝ω0(10)從(10)式可以看到，(9)式定義的瞬時頻率就是單頻余弦信號的頻率，即對于一個余弦信號，它的三參數(shù)小波變換得到的時間-尺度域結(jié)果會在某個能量最大的尺度附近形成一個尺度帶，但是這些尺度對應(yīng)小波系數(shù)通過(10)式計算出來的瞬時頻率都為余弦信號的頻率ω0，因此想象將這些尺度的能量都集中到ω0上，從而得到更集中的時頻表示。所述步驟三中，(9)式計算得到瞬時頻率后，得知能量往這個頻率上集中，也就是將小波變換系數(shù)累加到這個頻率成分上，即(a,b)→(ωs(a,b),b)，進(jìn)行能量重排。對于給定的信號s(t)∈L2(R)的三參數(shù)小波變換為Ws(a,b)，有如下表達(dá)式：∫0+∞Ws(a,b)a-1da=12π∫-∞+∞∫0+∞s^(ω)ψ^‾(aω;Λ)eiωba-1dadω=12π∫0+∞∫0+∞s^(ω)ψ^‾(aω;Λ)eiωba-1dadω=∫0+∞ψ^‾(ω;Λ)ωdω12π∫0+∞s^(ω)eiωbdω.---(11)]]>記而是s(t)的解析信號，從而(11)式可以寫為：∫0+∞Ws(a,b)a-1da=Cψsa(b).---(12)]]>因為sa(b)是s(b)的解析信號，故有從(12)中可以看到，通過小波變換系數(shù)再乘以因子a-1得到的結(jié)果和原信號的解析信號只差一個常數(shù)因子，這很容易得到其反變換，如(13)式所示；如果將Ws(a,b)a-1都分配給上面提到相應(yīng)的瞬時頻率成分上，則得到擠壓后的時頻分布，由此得出時間-尺度域到時間-頻率域的映射如下：Ts(ω,b)=∫{a:a>0,ωs(a,b)=ω,Ws(a,b)≠0}Ws(a,b)a-1da.---(14)]]>通過(14)式，在某個固定的時刻b，計算其瞬時頻率ωs(a,b)，將所有瞬時頻率都為某一頻率ω的小波系數(shù)通過(14)累加到一起，這樣就完成了能量的重分配，得到了擠壓后的時頻分布，(14)的等價形式為：Ts(ω,b)=∫{a:a>0,Ws(a,b)≠0}Ws(a,b)a-1δ(ωs(a,b)-ω)da.---(15)]]>所述步驟四中，為了更有效地確定有效信號的能量分布空間，提出了對重排后的時間頻率域系數(shù)運用百分比閾值策略進(jìn)行噪聲壓制，即軟閾值操作，定義為S(Ts)=Ts-δTs|Ts|,|Ts|≥δ0,|Ts|<δ,---(16)]]>其中閾值參數(shù)δ由下式計算：δ＝p·max{|Ts|}，(17)其中p為百分比，Ts為變換后的時間頻率域系數(shù)。所述步驟五中，確定有效信號對應(yīng)的能量分布空間后，就可以利用穩(wěn)定性更強(qiáng)、更為稀疏的時頻域變換系數(shù)，計算加權(quán)瞬時頻率：f(t)=∫-∞∞f*Ts(2πf,b)ξ(t-b)df∫-∞∞Ts(2πf,b)ξ(t-b)df,---(18)]]>其中f為頻率，b和t為時間，ξ(t)為漢明窗。與現(xiàn)有技術(shù)相比，本發(fā)明在三參數(shù)小波變換的基礎(chǔ)上，通過估計重排準(zhǔn)則(瞬時頻率)后進(jìn)行能量重排，得到更加準(zhǔn)確、更加稀疏的時頻表示，以更加準(zhǔn)確地確定有效信號能量分布空間；在此基礎(chǔ)上引入閾值去噪，進(jìn)行噪聲壓制；最后，給出了含噪信號加權(quán)瞬時頻率的估計方法。本發(fā)明通過計算無噪和含噪合成地震記錄的瞬時頻率，對比不同方法的抗噪性能，本文方法的計算結(jié)果具有良好的抗噪性能和精度，將該方法應(yīng)用于實際資料，能夠更加清晰地刻畫儲層的特征。附圖說明圖1為本發(fā)明40Hz的Ricker子波的時頻域譜圖；其中，(a)無噪的Ricker子波，(b)、(c)無噪Ricker子波的小波變換和本發(fā)明提出變換的時頻圖；(d)加入信噪比為5dB的高斯白噪聲后的Ricker子波，(e)、(f)含噪Ricker子波的小波變換和本發(fā)明提出變換的時頻圖；圖2為本發(fā)明無噪合成地震記錄和采用不同方法計算出來的瞬時頻率；其中，(a)合成的真記錄，(b)基于HT計算的瞬時頻率，(c)、(d)分別為基于Fourier變換和同步擠壓變換計算的加權(quán)平均瞬時頻率；圖3為本發(fā)明含噪合成地震記錄和采用不同方法計算出來的瞬時頻率；其中，(a)含噪的合成地震記錄，(b)基于HT計算的瞬時頻率，(c)、(d)分別為基于Fourier變換和同步擠壓變換計算的加權(quán)平均瞬時頻率；圖4為本發(fā)明2D實際地震記錄剖面；其中，箭頭標(biāo)注的位置為三角洲砂體疊置，因地震數(shù)據(jù)分辨率較低和噪聲影響，并不能反映疊置特征；圖5為本發(fā)明用不同方法計算得到的2D實際地震數(shù)據(jù)的瞬時頻率；其中，(a)基于HT計算的瞬時頻率，(b)、(c)、(d)分別為基于Fourier變換、小波變換和同步擠壓變換計算的加權(quán)平均瞬時頻率。具體實施方式下面結(jié)合附圖和實施例對本發(fā)明做進(jìn)一步說明。(1)小波變換定義實地震道s(t)的小波變換定義為Ws(a,b)=a-1∫s(t)ψ‾(t-ba)dt=12π∫-∞∞S(ω)ψ^‾(aω)eiωbdω,---(1)]]>其中，為ψ(t)的復(fù)共軛，a為尺度因子，b為時移因子，是基本小波ψ(t)的Fourier變換，是s(t)的Fourier變換。基本小波平方可積，無直流分量，滿足以下容許性條件：Cψ=∫-∞∞|ψ(ω)|2|ω|dω<∞.---(2)]]>針對薄互層地震信號含有快速變化的振幅和頻率分量的特點，以高靜懷提出的最佳匹配地震子波(BMSW)的小波為基礎(chǔ)，借鑒Harrop等人的工作對傳統(tǒng)的Morlet小波加以改進(jìn)，高靜懷等人提出了一種小波函數(shù)——三參數(shù)小波(TPW)。由于具有三個參數(shù)，通過調(diào)節(jié)，它可以靈活的適應(yīng)我們的需要，將它應(yīng)用于層序檢測和高精度地震資料分辨中都取得了很明顯的效果。本發(fā)明選取三參數(shù)小波作為小波變換的母小波。三參數(shù)小波的定義為：ψ(t;σ,τ,β)=e-τ(t-β)2{p(σ,τ,β)[cos(σt)-k(σ,τ,β)]+iq(σ,τ,β)sin(σt)},---(3)]]>其中，τ為能量衰減因子，β為能量延遲時間，σ為分析小波調(diào)制頻率，(σ,τ,β∈R且σ,τ＞0)。為了書寫方便，用向量Λ＝(σ,τ,β)記參數(shù)集合，則ψ(t；σ,τ,β)可記為ψ(t；Λ)，其他量類似。以向量的形式，上式可簡寫為：ψ(t;Λ)=e-τ(t-β)2{p(Λ)[cos(σt)-k(Λ)]+iq(Λ)sin(σt)},---(4)]]>其中，k(Λ)=e-σ24τ[cos(βσ)+iq(Λ)p(Λ)sin(βσ)],p(Λ)=(2τπ)14[4(e-σ22τ-e-3σ28τ)×cos2(βσ)+1-e-σ22τ]-12,q(Λ)=(2τπ)14[4(e-σ22τ-e-3σ28τ)×sin2(βσ)+1-e-σ22τ]-12.---(5)]]>對(4)式作Fourier變換，得到其頻域形式為ψ^(ω;Λ)=∫-∞∞ψ(t;Λ)e-iωtdt=πτp(Λ)+q(Λ)2e-iβ(ω-σ)-(ω-σ)24τ+πτp(Λ)-q(Λ)2e-iβ(ω+σ)-(ω+σ)24τ-πτp(Λ)k(Λ)e-iβω-ω24τ.---(6)]]>(2)計算重排準(zhǔn)則(瞬時頻率)假設(shè)一個單頻余弦信號，即s(t)＝Acωo0s(，ts)(t)的Fourier變換為根據(jù)(1)式，信號s(t)關(guān)于三參數(shù)小波ψ(t；Λ)的小波變換結(jié)果為Ws(a,b)=12π∫-∞+∞s^(ω)ψ^‾(aω;Λ)eiωbdω=A2π[ψ^‾(aω0;Λ)eiω0b+ψ^‾(-aω0;Λ)e-iω0b].---(7)]]>因為我們選擇的三參數(shù)小波幾乎沒有負(fù)頻率分量，即且只考慮正尺度，則(7)式可以簡化為Ws(a,b)=A2πψ^‾(aω0;Λ)eiω0b.---(8)]]>從(8)式可以看到，假如三參數(shù)小波ψ(t；Λ)的峰值頻率為ξM，則小波變換的結(jié)果將在尺度a＝ξM/ω0處取到最大值，并以這個能量最大的尺度為中心形成一個尺度帶，造成能量的擴(kuò)散，為了得到更集中的時頻分布，計算重排準(zhǔn)則——瞬時頻率：ωs(a,b)=∂bWs(a,b)2πiWs(a,b),---(9)]]>其中，Ws(a,b)≠0。當(dāng)Ws(a,b)0＝時，信號在此處無能量，(9)式將無窮大，計算瞬時頻率ωs(a,b)沒有意義，所以計算時要求當(dāng)噪聲較小時，可以選擇閾值如果噪聲較大，需要抑制噪聲，我們選取自適應(yīng)閾值作為閾值其中ση是和噪聲水平相關(guān)的，它是由小波變換系數(shù)的前nv個尺度來估計的，median是matlab中的庫函數(shù)，代表取中值，0.6745是針對高斯噪聲的正則化因子。假設(shè)有N個時刻，所以可以對N個時刻求得的ση取均值，能得到最優(yōu)的閾值。當(dāng)參考信號是單頻余弦信號時，其小波變換會形成一個尺度帶，但是可以由(9)計算其瞬時頻率為ωs(a,b)＝ω0(10)從(10)式可以看到，(9)式定義的瞬時頻率就是單頻余弦信號的頻率，即對于一個余弦信號，它的三參數(shù)小波變換得到的時間-尺度域結(jié)果會在某個能量最大的尺度附近形成一個尺度帶，但是這些尺度對應(yīng)小波系數(shù)通過(10)式計算出來的瞬時頻率都為余弦信號的頻率ω0，因此我們可以想象將這些尺度的能量都集中到ω0上，從而得到更集中的時頻表示。(3)能量重排——時間尺度域能量映射到時間頻率域(9)式計算得到瞬時頻率后，我們知道能量應(yīng)該往這個頻率上集中，也就是將小波變換系數(shù)累加到這個頻率成分上，即(a,b)→(ωs(a,b),b)，進(jìn)行能量重排。對于給定的信號s(t)∈L2(R)的三參數(shù)小波變換為Ws(a,b)，有如下表達(dá)式：∫0+∞Ws(a,b)a-1da=12π∫-∞+∞∫0+∞s^(ω)ψ^‾(aω;Λ)eiωba-1dadω=12π∫0+∞∫0+∞s^(ω)ψ^‾(aω;Λ)eiωba-1dadω=∫0+∞ψ^‾(ω;Λ)ωdω12π∫0+∞s^(ω)eiωbdω.---(11)]]>記而是s(t)的解析信號，從而(11)式可以寫為：∫0+∞Ws(a,b)a-1da=Cψsa(b).---(12)]]>因為sa(b)是s(b)的解析信號，故有從(12)中可以看到，通過小波變換系數(shù)再乘以因子a-1得到的結(jié)果和原信號的解析信號只差一個常數(shù)因子，這很容易得到其反變換，如(13)式所示。如果我們將Ws(a,b)a-1都分配給我們上面提到相應(yīng)的瞬時頻率成分上，則就可以得到擠壓后的時頻分布，由此我們得出時間-尺度域到時間-頻率域的映射如下：Ts(ω,b)=∫{a:a>0,ωs(a,b)=ω,Ws(a,b)≠0}Ws(a,b)a-1da.---(14)]]>通過(14)式，在某個固定的時刻b，計算其瞬時頻率ωs(a,b)，將所有瞬時頻率都為某一頻率ω的小波系數(shù)通過(14)累加到一起，這樣就完成了能量的重分配，得到了擠壓后的時頻分布，(14)的等價形式為：Ts(ω,b)=∫{a:a>0,Ws(a,b)≠0}Ws(a,b)a-1δ(ωs(a,b)-ω)da.---(15)]]>(4)有效信號能量分布空間的確定當(dāng)選用合適的小波函數(shù)(本發(fā)明選取三參數(shù)小波)對含噪信號進(jìn)行小波變換，將其投影到時間尺度域的時候，有效信號的能量將會分布在較小的子空間V，噪聲的能量會擴(kuò)散到比較大的子空間V′，甚至整個時間-尺度域。換句話說，有效信號和少數(shù)的系數(shù)相對應(yīng)，而噪聲幾乎分布在全部的系數(shù)上。如果我們確定了有效信號對應(yīng)的子空間V，得到有效信號對應(yīng)的系數(shù)，將其它系數(shù)置零，則噪聲就會在變換域被壓制，信噪比得到提高。本發(fā)明通過同步擠壓變換之后，時間尺度域系數(shù)經(jīng)過重排，變換到時間頻率域，得到更為稀疏的時頻結(jié)果，有效信號能量分布空間變小，更有利于噪聲壓制。為了更有效地確定有效信號的能量分布空間，我們提出了對重排后的時間頻率域系數(shù)運用百分比閾值策略進(jìn)行噪聲壓制，即軟閾值操作，定義為S(Ts)=Ts-δTs|Ts|,|Ts|≥δ0,|Ts|<δ,---(16)]]>其中閾值參數(shù)δ由下式計算：δ＝p·max{|Ts|}，(17)其中p為百分比，Ts為本發(fā)明變換后的時間頻率域系數(shù)。通過該濾波策略，得到有效信號對應(yīng)的時頻域系數(shù)，這些系數(shù)構(gòu)成噪聲被壓制以后的信號對應(yīng)的時間-頻率域譜圖。圖1(a)畫出了不含噪的40Hz的Ricker子波的譜圖，圖1(b)和(c)為小波變換和本發(fā)明變換后的時頻譜圖，可見本發(fā)明變換能夠得到更加集中、更加稀疏的時頻表示。圖1(d)為含噪的40Hz的Ricker子波(信噪比SNR＝5dB)的譜圖，可見它已經(jīng)被噪聲污染，圖1(e)和(f)為運用濾波策略后得到的小波變換和本發(fā)明變換的時頻譜圖，可以看出由于本發(fā)明變換更為稀疏，噪聲得到明顯壓制，有效信號被清晰地表現(xiàn)出來。(5)加權(quán)平均瞬時頻率提取確定有效信號對應(yīng)的能量分布空間后，就可以利用本發(fā)明得到的穩(wěn)定性更強(qiáng)、更為稀疏的時頻域變換系數(shù)，計算加權(quán)瞬時頻率：f(t)=∫-∞∞f*Ts(2πf,b)ξ(t-b)df∫-∞∞Ts(2πf,b)ξ(t-b)df,---(18)]]>其中f為頻率，b和t為時間，ξ(t)為漢明窗。實施例：(1)合成信號算例我們對30Hz的Ricker子波同反射系數(shù)序列卷積生成的合成地震記錄進(jìn)行研究。圖2(a)為不含噪的合成記錄，(b)、(c)、(d)分別為Hilbert變換計算的瞬時頻率、基于Fourier變換和本發(fā)明計算得到的加權(quán)瞬時頻率?？梢?b)和(c)的瞬時頻率估計結(jié)果起伏較大，不夠穩(wěn)定；本發(fā)明的估計結(jié)果更為穩(wěn)定，瞬時頻率主要分布在Ricker子波主頻30Hz附近、20Hz-40Hz范圍內(nèi)。圖3(a)為加上高斯白噪聲的含噪記錄(信噪比SNR＝5dB)，(b)、(c)、(d)分別為Hilbert變換計算的瞬時頻率、基于Fourier變換和本發(fā)明計算得到的加權(quán)瞬時頻率。從圖3(b)可以看出，通過傳統(tǒng)的Hilbert變換方法得到的瞬時頻率結(jié)果中，有效信號的瞬時頻率因為噪聲的存在被完全掩蓋。基于Fourier加權(quán)平均的瞬時頻率估計結(jié)果，雖然具有一定的抗噪性，但是估計結(jié)果仍然不夠穩(wěn)定。通過本發(fā)明方法計算的結(jié)果(圖3(d))仍然能夠清晰地、穩(wěn)定地顯示有效信號的瞬時頻率。(2)實際地震資料算例為了進(jìn)一步檢驗本發(fā)明方法的有效性，我們將其用于實際地震資料。圖4是三維數(shù)據(jù)體的一個剖面，該資料隨機(jī)噪聲十分嚴(yán)重，嚴(yán)重影響儲層預(yù)測。根據(jù)井?dāng)?shù)據(jù)信息揭示，層位H6和H7之間為三角洲相，紅色箭頭所指處為三角洲砂體疊置，但受地震分辨率和噪聲限制，無法識別三角洲砂體邊界和疊置。圖5(a-d)分別是采用Hilbert變換計算的瞬時頻率和基于Fourier變換、小波變換、本發(fā)明方法計算的加權(quán)瞬時頻率。圖5(a)中可以看出三角洲砂體疊置，但由于信噪比較低導(dǎo)致砂體不夠連續(xù)。圖5(b)中由于引入窗函數(shù)，在一定程度上降低了噪聲影響，但同時分辨率也隨之受影響，并不能很好地刻畫三角洲砂體疊置現(xiàn)象。圖5(c)在一定程度上壓制了噪聲，檢測精度相對較好，但是由于受限于時頻分辨率不高，噪聲壓制不夠明顯，剖面信噪比仍舊較低。圖5(d)結(jié)果抗噪性明顯增強(qiáng)，估計得到的瞬時頻率信噪比大大改善，三角洲砂體的疊置現(xiàn)象更加清晰，可見本方法是檢測巖性的有效手段。當(dāng)前第1頁1 2 3