專利名稱：一種魯棒的地震信號(hào)瞬時(shí)頻率的估計(jì)方法
技術(shù)領(lǐng)域：
本發(fā)明涉及信號(hào)處理領(lǐng)域，具體地說涉及到地震信號(hào)處理領(lǐng)域。

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

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

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

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


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

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

為

利用加權(quán)振幅幾何級(jí)數(shù)的方法計(jì)算新解析信號(hào)的頻率域函數(shù)


其中N為加權(quán)階數(shù)，max(·)為取最大值函數(shù)； 步驟(3.4)利用快速傅里葉逆變換計(jì)算xi(t)對(duì)應(yīng)的新的解析信號(hào)


IFFT為快速傅立葉逆變換； 步驟(3.5)對(duì)于時(shí)間點(diǎn)i，得到多個(gè)解析信號(hào)的估計(jì)

k＝i-L+1，i-L+2，…，i，將該點(diǎn)對(duì)應(yīng)多個(gè)估計(jì)的均值作為該點(diǎn)對(duì)應(yīng)的新的解析信號(hào)

其中n為時(shí)間點(diǎn)i所估計(jì)出的解析信號(hào)的個(gè)數(shù)，

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

該解析信號(hào)表示為

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

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

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



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

具體實(shí)施例方式 普通的希爾伯特變換是將一個(gè)信號(hào)進(jìn)行-90°的相移?？梢詫⑵淅斫鉃閷⒁粋€(gè)信號(hào)通過一個(gè)全通濾波器

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

綜合起來即為

對(duì)于一個(gè)觀測(cè)信號(hào)對(duì)應(yīng)的解析信號(hào)的頻率域響應(yīng)，可以通過如下方法來實(shí)現(xiàn)給定一個(gè)信號(hào)x(t)，其頻率域響應(yīng)為

它的解析信號(hào)表示為

其中

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

而

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

計(jì)算yi(t)的頻率域響應(yīng)

利用(1)式求出該信號(hào)對(duì)應(yīng)的解析信號(hào)的頻率域響應(yīng)

變換該頻域函數(shù)

其中

是振幅譜，

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

為信號(hào)xi(t)對(duì)應(yīng)的新的解析信號(hào)的頻率域響應(yīng)，對(duì)其進(jìn)行傅里葉逆變換即可得到xi(t)對(duì)應(yīng)的新的解析信號(hào)


其中IFFT為快速傅里葉逆變換。由于是逐點(diǎn)選取長(zhǎng)度為L(zhǎng)的觀測(cè)信號(hào)的，因此對(duì)于每個(gè)時(shí)間點(diǎn)i上會(huì)出現(xiàn)多個(gè)解析信號(hào)的估計(jì)

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

之后，可以將其寫成

其中

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