專利名稱：一種魯棒的地震信號瞬時頻率的估計方法
技術領域：
本發(fā)明涉及信號處理領域，具體地說涉及到地震信號處理領域。

背景技術：
地震勘探是一種利用人工地震技術激發(fā)地震波，然后通過地面布置的檢波器得到由地下地層分界面反射回來的聲波，通過對數(shù)據(jù)進行處理，提取有用的地震屬性以查明地下的地質構造或巖性特征，為尋找油氣藏或其他勘探目的提供重要的工具。
地震信號的瞬時頻率估計是地震信號的屬性提取中的關鍵問題之一，是其他一些屬性提取的基礎。地震信號的瞬時頻率估計在識別巖性和儲層分析方面具有重要的作用，因此準確估計地震信號的瞬時頻率具有重要的應用價值。
目前，常用的瞬時頻率估計方法是基于希爾伯特變換的解析信號方法(B.Boashash，“Estimating and interpreting the instantaneous frequency of asignal-PartlFundamentals，”Proc.IEEE，vol.80，no.4，pp.520-538，1992.4)。該方法介紹如下解析信號(又稱復信號)由實部和虛部組成。它的實部是原信號，虛部是原信號的希爾伯特變換。對于一個普通的信號x(t)＝A(t)cos(θ(t))，其中A(t)是振幅，θ(t)是相位，則其希爾伯特變換為

根據(jù)解析信號的定義，解析信號f(t)可以表示為

利用歐拉公式可以得到f(t)＝A(t)exp(iθ(t))，其中exp(·)表示自然指數(shù)函數(shù)(下同)，i為虛數(shù)單位。瞬時頻率S(t)定義為相位的時間變化率

(J.Ville，“Theorie etapplication de la notion de signal analytic”，Cables et Transmissions，vol.2A，pp.61-74，1948)。至此就是利用普通的希爾伯特變換的方法估計信號的瞬時頻率的全部過程。
該方法在信噪比較高的環(huán)境下可以得到可信度較高的結果，但是在信噪比較低的環(huán)境下它的結果受噪聲的影響很嚴重。在努力降低噪聲對結果帶來的影響方面，Luo Yi曾經(jīng)用一種稱為廣義希爾伯特變換的方法來估計信號的瞬時相位，從而更好的顯示河道和斷裂帶(Luo Yi.Seismic data processing method toenhance fault and channel displayUS006963804[P].2005-11-08)，這種技術對噪聲環(huán)境下瞬時相位的估計有一定的效果。
針對上述現(xiàn)有技術存在的問題，本發(fā)明的目的是提出一種魯棒的地震信號瞬時頻率的估計方法。使用該方法能夠信噪比低的情況下能夠準確的估計出信號的瞬時頻率。


發(fā)明內容
本發(fā)明采用的技術方案是一種魯棒的地震信號瞬時頻率的估計方法，該方法包括以下步驟 步驟(1)利用地震勘探設備，采集地震數(shù)據(jù) 在反射波地震勘探中采用一點放炮多點接收的常規(guī)方法，在預勘探區(qū)域設置多個接收點，所述多個接收點散布在包括炮點在內的地面上一個勘探目標區(qū)域范圍的矩形網(wǎng)格點上，使用地震檢波設備獲取爆炸波及地震波的波形； 步驟(2)對步驟(1)中獲取的波形數(shù)據(jù)進行地震資料的處理，經(jīng)過反褶積、抽取共中心點道集、動校正、疊加和偏移過程之后獲得疊后的每道數(shù)據(jù)x(t)，t＝1，2，…，H，采樣時間間隔T，T為1ms、2ms或4ms，設定 窗長L和加權階數(shù)N，其中L取20～30，N取2～5； 步驟(3)計算機按以下步驟計算新的解析信號
步驟(3.1)沿時間點i＝1，2，…，H-L+1逐點取出一段窗長為L的觀測信號zi＝{x(i+j-1)，j＝1，2，…，L}，進行加窗處理yi＝zi·hi，其中hi為長度L的高斯窗函數(shù)，求出該信號的頻率域響應

其中FFT為快速傅立葉變換； 步驟(3.2)按下述公式計算該信號對應的解析信號的頻率域響應
步驟(3.3)變換

為

利用加權振幅幾何級數(shù)的方法計算新解析信號的頻率域函數(shù)


其中N為加權階數(shù)，max(·)為取最大值函數(shù)； 步驟(3.4)利用快速傅里葉逆變換計算xi(t)對應的新的解析信號


IFFT為快速傅立葉逆變換； 步驟(3.5)對于時間點i，得到多個解析信號的估計

k＝i-L+1，i-L+2，…，i，將該點對應多個估計的均值作為該點對應的新的解析信號

其中n為時間點i所估計出的解析信號的個數(shù)，

表示求和； 步驟(4)利用新的解析信號

該解析信號表示為

按下面的公式求出地震信號的瞬時相位θ(t)，

tan-1(·)表示反正切函數(shù)； 步驟(5)相位展開，如果|θ(i+1)-θ(i)|＞σ，θ(i+1)＝θ(i+1)+2πk，使得相鄰兩個相位之間的差值不超過σ，其中σ是閾值，取為π，k為任意整數(shù)； 步驟(6)利用估計出來的瞬時相位，利用相位差分方法來計算地震信號的瞬時頻率S(t)，采用下面的這個公式

其中T為采樣 時間間隔，單位是秒。
由上述本發(fā)明采用的技術方案可以看出，本發(fā)明是在地震觀測信號中逐點選取一定長度的信號進行處理，然后再估計整個觀測信號對應的解析信號；在解析信號的頻率域對振幅進行了幾何級數(shù)的加權，這種方法很好的壓制了噪聲；對于在每個時間點上得到的多個估計，利用多個估計的均值作為該點的解析信號的值，使得估計出來的新的解析信號更加的準確。這樣通過解析信號估計地震信號的瞬時相位進而估計地震信號的瞬時頻率將更加的準確。
圖1是本發(fā)明方法的流程示意圖。
在我們的仿真實驗中，分別產(chǎn)生兩種類型的信號，第一種是普通的正弦信號，分別在有噪聲和沒有噪聲的情況下；第二種是普通的chirp(調頻)信號，分別在有噪聲和沒有噪聲的情況下。在比較中，分別給出了原信號，原信號真實的瞬時頻率，用普通的希爾伯特變換得到的瞬時頻率，以及用本發(fā)明方法得到的瞬時頻率。
圖2、圖3、圖4、圖5分別給出了不同情況下兩種方法的結果?？梢钥闯鲈跊]有噪聲的情況下，本方法和普通的方法一樣能準確的估計出信號的瞬時頻率；在有噪聲的情況下，普通的方法無法有效的估計出信號的瞬時頻率，但是本發(fā)明方法卻能夠準確的估計出信號的瞬時頻率。



圖1、本發(fā)明所屬方法的計算機程序流程圖。
圖2、正弦信號在沒噪聲情況下的瞬時頻率估計 a，瞬時頻率為100Hz的正弦信號； b，正弦信號的真實瞬時頻率； c，利用普通希爾伯特變換得到的瞬時頻率； d，利用本發(fā)明方法得到的瞬時頻率。
圖3、正弦信號在有噪聲情況下的瞬時頻率估計 a，添加了15dB高斯白噪聲的瞬時頻率為100Hz的正弦信號； b，正弦信號的真實瞬時頻率； c，利用普通希爾伯特變換得到的瞬時頻率； d，利用本發(fā)明方法得到的瞬時頻率。
圖4、調頻信號在無噪聲情況下的瞬時頻率估計 a，瞬時頻率從100Hz到500Hz的二次調頻信號； b，調頻信號的真實瞬時頻率； c，利用普通希爾伯特變換得到的瞬時頻率； d，利用本發(fā)明方法得到的瞬時頻率。
圖5、調頻信號在有噪聲情況下的瞬時頻率估計 a，添加了15dB高斯白噪聲的瞬時頻率從100Hz到500Hz的二次調頻信號； b，正弦信號的真實瞬時頻率； c，利用普通希爾伯特變換得到的瞬時頻率； d，利用本發(fā)明方法得到的瞬時頻率。

具體實施例方式 普通的希爾伯特變換是將一個信號進行-90°的相移?？梢詫⑵淅斫鉃閷⒁粋€信號通過一個全通濾波器

得到該信號的希爾伯特變換。該濾波器具有如下特性

綜合起來即為

對于一個觀測信號對應的解析信號的頻率域響應，可以通過如下方法來實現(xiàn)給定一個信號x(t)，其頻率域響應為

它的解析信號表示為

其中

為它的希爾伯特變換，i為虛數(shù)單位；它的頻域表達式為

而

綜合起來就是
即我們對該觀測信號先求解它的快速傅里葉變換，然后將正頻率的幅值乘以2，負頻率的幅值置0即可得到該觀測信號對應的解析信號的頻率域響應。
具體來說，本發(fā)明的方法是在一個范圍內逐點依次對整個信號中的每一小段地震信號進行解析信號的估計，然后再估計整段信號的解析信號。對于獲得的地震觀測信號x(t)，t＝1，2，…，H，沿時間點i＝1，2，…，H-L+1依次取窗長L＝20的觀測信號zi(t)＝{x(i)，x(i+1)，…，x(i+L-1)}，進行加窗處理yi＝zi·hi，其中hi為長度L的高斯窗函數(shù)，

計算yi(t)的頻率域響應

利用(1)式求出該信號對應的解析信號的頻率域響應

變換該頻域函數(shù)

其中

是振幅譜，

是相位譜。對獲得的頻域函數(shù)進行振幅加權處理
其中加權階數(shù)N＝3。新獲得的頻率域函數(shù)

為信號xi(t)對應的新的解析信號的頻率域響應，對其進行傅里葉逆變換即可得到xi(t)對應的新的解析信號


其中IFFT為快速傅里葉逆變換。由于是逐點選取長度為L的觀測信號的，因此對于每個時間點i上會出現(xiàn)多個解析信號的估計

其中k一般為i-L+1，i-L+2，…，i，對于多個估計，利用它們的均值作為該點上的解析信號的值 n為每個點上獲得的估計信號的個數(shù)。獲得了新的解析信號

之后，可以將其寫成

其中

在獲得信號的瞬時相位之后，由于求出的相位的范圍在(-π，π]之間，相鄰的兩個相位會出現(xiàn)跳變，使得相鄰的相位不連續(xù)，因此需要對相位進行展開，使得相鄰的兩個相位之間的差值不超過一個固定的閾值σ(J.M.Tribolet，a new phase unwrapping algorithm，IEEE Transactionson acoustics，speech，and signal processing，vol.ASSP-23，No.2，1977.4)，這里閾值σ＝π。
目前在利用解析信號計算其瞬時頻率方面，出現(xiàn)了多種計算方法(E.Barnes，The calculation of instantaneous frequency and instantaneous bandwidthGeophysics，57，1992.11)，有基于兩點的有限脈沖響應的微分器(Scheuer.T.E，Oldenburg.D.W，Local phase velocity from complex seismic data，Geophysics，1988，53pp.1503)，有基于三點的有限脈沖響應的微分器(Boashash.B，O’Shea.P，Ristic.B，Statisticalcomputation comparison of some estimators for instantaneous frequency，Proc IEEEICASSP-91，1991，3193～3196)，相位差分法包括前向有限差分法、后向有限差分法、中心有限差分法(Boashash.B，Estimating and interpreting the instantaneousfrequency of a signal-Part2algorithms and applications.Proc.IEEE，vol.80，no.4，pp.540-568，1992.4)。在本發(fā)明里利用相位后向差分法計算信號的瞬時頻率S(t)，即利用求出的相位然后進行后向差分
普通的希爾伯特變換是對整個信號進行變換的，在頻域沒有做任何處理；而本發(fā)明的方法是先逐點選擇一小段信號，然后對每一段信號進行如下處理先進行加窗處理，然后進行頻率域的振幅幾何級數(shù)加權，接著利用傅里葉逆變換獲得新的解析信號。再進行求均值獲得魯棒的估計結果。該方法的好處是由于噪聲在頻率域表現(xiàn)為一些很小的值，這樣利用振幅幾何級數(shù)加權時可以大大的壓制噪聲；同時反變換后對同一個點上多個值進行求平均，因此這樣得到的結果將會更加的魯棒，受噪聲的干擾小。
總之，本發(fā)明適用于地震信號的瞬時頻率的估計，不僅可以在高信噪比的環(huán)境中準確的估計出信號的瞬時頻率，在低信噪比的環(huán)境下仍然能準確的估計出信號的瞬時頻率。
權利要求
1.一種魯棒的地震信號瞬時頻率的估計方法，其特征在于，該方法依次含有以下步驟
步驟(1)利用地震勘探設備，采集地震數(shù)據(jù)
在反射波地震勘探中采用一點放炮多點接收的方法，在預勘探區(qū)域設置多個接收點，所述多個接收點散布在包括炮點在內的地面上一個勘探目標區(qū)域范圍的矩形網(wǎng)格點上，使用地震檢波設備獲取爆炸波及地震波的波形；
步驟(2)對步驟(1)中獲取的波形數(shù)據(jù)進行地震資料的常規(guī)處理，經(jīng)過反褶積、抽取共中心點道集、動校正、疊加和偏移過程之后獲得疊后的每道數(shù)據(jù)x(t)，t＝1，2，…，H，采樣時間間隔T，T為1ms、2ms或4ms，設定窗長L和加權階數(shù)N；
步驟(3)計算機按以下步驟計算新的解析信號
步驟(3.1)沿時間點i＝1，2，…，H-L+1逐點取出一段窗長為L的觀測信號zi＝{x(i+j-1)，j＝1，2，…，L}，進行加窗處理yi＝zi·hi，其中hi為長度L的高斯窗函數(shù)，求出該信號的頻率域響應
其中FFT為快速傅立葉變換；
步驟(3.2)按下述公式計算該信號對應的解析信號的頻率域響應
步驟(3.3)變換
為
利用加權振幅幾何級數(shù)的方法計算新解析信號的頻率域函數(shù)
其中N為加權階數(shù)，max(·)為取最大值函數(shù)；
步驟(3.4)利用快速傅里葉逆變換計算xi(t)對應的新的解析信號
IFFT為快速傅立葉逆變換；
步驟(3.5)對于時間點i，得到多個解析信號的估計
k＝i-L+1，i-L+2，…，i，將該點對應多個估計的均值作為該點對應的新的解析信號
其中n為時間點i所估計出的解析信號的個數(shù)，
表示求和；
步驟(4)利用新的解析信號
該解析信號表示為
按下面的公式求出地震信號的瞬時相位θ(t)，
tan-1(·)表示反正切函數(shù)；
步驟(5)相位展開，如|θ(i+1)-θ(i)|＞σ，θ(i+1)＝θ(i+1)+2πk，使得相鄰兩個相位之間的差值不超過σ，其中σ是閾值，取為π，k為任意整數(shù)；
步驟(6)利用估計出來的瞬時相位，利用相位差分方法來計算地震信號的瞬時頻率S(t)，采用下面的這個公式
其中T為采樣時間間隔，單位是秒。
2.根據(jù)權利要求1所述的一種魯棒的瞬時頻率的估計方法，其特征在于，加權階數(shù)N為2～5之間的整數(shù)。
3.根據(jù)權利要求1所述的一種魯棒的瞬時頻率的估計方法，其特征在于，濾波窗長L根據(jù)采樣時間間隔來選取，L為20～30。
4.根據(jù)權利要求1所述的一種魯棒的瞬時頻率的估計方法，其特征在于，所選擇的窗函數(shù)是矩形窗、高斯窗、漢寧窗、漢明窗、三角窗之一。
全文摘要
一種魯棒的地震信號瞬時頻率的估計方法屬于地震信號處理技術領域。其特征在于該方法包括下列步驟利用地震勘探設備，采集地震數(shù)據(jù)；獲得疊后地震數(shù)據(jù)；根據(jù)每道地震數(shù)據(jù)的長度，設定窗長和加權階數(shù)；逐點取一段窗長的觀測信號，進行加窗處理，計算該信號對應的解析信號的頻率域響應；再對其進行振幅的幾何級數(shù)加權，然后利用快速傅里葉逆變換得到新的解析信號；每個時間點會得到多段解析信號的估計值，利用它們的均值作為該點對應的新的解析信號；利用新的解析信號計算信號的瞬時相位，再計算瞬時頻率。本發(fā)明不僅可以準確地計算地震信號的瞬時頻率，而且信噪比低的情況下依然能夠準確的估計瞬時頻率。
文檔編號G01V1/28GK101825722SQ20101013335
公開日2010年9月8日 申請日期2010年3月25日 優(yōu)先權日2010年3月25日
發(fā)明者陸文凱, 張長開 申請人:清華大學