一種地震信號的分數(shù)域局部功率譜計算方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于地震信號處理領(lǐng)域。方法基于最優(yōu)分數(shù)階傅里葉變換，將現(xiàn)代功率譜 估計與分數(shù)域時頻分析技術(shù)相結(jié)合，提出了新的局部功率譜計算方法。通過該方法可以提 高信號功率譜圖的時頻分辨率，可以在地震信號處理中獲得高分辨率的局部功率譜特征， 為后續(xù)儲層流體預(yù)測和識別創(chuàng)造有利條件。 技術(shù)背景
[0002] 分數(shù)階傅里葉變換起源于1929年Wiener的研宄，興起于1980年Namias的研宄， 當時他利用特征值的任意次冪運算明確的提出了階數(shù)為分數(shù)的傅里葉變換的概念。20世 紀90年代以后，分數(shù)傅里葉變換作為一種新型的信號分析工具受到越來越多的關(guān)注，并在 光學領(lǐng)域最先得到廣泛應(yīng)用。1994年Almeida揭示了分數(shù)傅里葉變換和傳統(tǒng)時頻分析工具 的關(guān)系，得出了分數(shù)傅里葉變換可以解釋為信號在時頻面上坐標軸繞原點逆時針旋轉(zhuǎn)的重 要結(jié)論。至此分數(shù)傅里葉變換被賦予明確的物理意義，吸引著越來越多的研宄者參與相關(guān) 研宄。分數(shù)階傅里葉變換較傳統(tǒng)傅里葉變換的一個顯著優(yōu)點就是能為時頻分析提供更大的 選擇余地，以獲得某些性能上的改進。
[0003] 功率譜描述了隨機過程中各頻率成分平均功率的大小，在工程上廣泛應(yīng)用于隨機 信號處理，是隨機信號的重要表征形式之一。它源于1898年Schuster提出的周期圖概念。 20世紀50年代后，功率譜估計方法發(fā)展迅速，間接法、AR模型法、Burg譜估計法、特征空間 分解法等相繼出現(xiàn)。到現(xiàn)在信號功率譜密度估計方面已有許多新的算法，這些新方法連同 演變出來的各種算法不下幾十種。這些方法作為一種功率譜估計手段在生物醫(yī)學、語音、雷 達、地震勘探等領(lǐng)域得到了廣泛應(yīng)用。
[0004] 然而功率譜估計給出的是一種基于信號全局特征的頻譜信息，這只對平穩(wěn)信號產(chǎn) 生較好的分析效果。對于非平穩(wěn)信號，在不同的時間段包含的頻譜信息是不一樣的，若對 信號直接做功率譜分析，就無法得到信號頻譜隨時間變化的特征。1946年Gabor提出了 一種對信號局部加窗后再做傅里葉變換的方法即Gabor變換。這種變換能夠很好的反映 出信號局部的頻譜信息，于是得到了快速發(fā)展。隨后短時傅里葉變換、Wigner-Vile分布、 Cohen分布、小波變換、S變換等時頻分析技術(shù)相繼出現(xiàn)。功率譜圖與短時傅里葉變換誕 生于同一時期，它可由短時傅里葉變換幅值的平方求得，它是一種不同于短時傅里葉變換、 Wigner-Vile分布、小波變換等，且能夠反映信號在時間和頻率上的能量分布特性的一種優(yōu) 秀的時頻分析技術(shù)。
[0005] 在地震勘探領(lǐng)域，利用地震資料直接進行儲層流體識別，已成為地震儲層預(yù)測技 術(shù)的主流發(fā)展方向之一。頻譜分解技術(shù)是其中一種行之有效的技術(shù)，它主要利用不同性質(zhì) 的流體對地震波的吸收、衰減不同這一特點，通過地震信號不同層位下的頻率屬性來反映 地下傳輸介質(zhì)的物理特性，從而為判斷該層位的油氣特性提供依據(jù)。地震信號作為一種非 平穩(wěn)信號，頻譜分解的核心是時頻分析技術(shù)?，F(xiàn)有的時頻分析技術(shù)因為時頻分辨率、交叉項 等問題制約著其在地震信號處理中的應(yīng)用效果，同時也激發(fā)了人們尋求新的時頻分析方法 的熱情。
[0006] 2002年Lufiye等人將分數(shù)域傅里葉變換的思想引入到時頻分析中，有效地提高 了信號的時頻分析效果。分數(shù)域時頻分析在傳統(tǒng)時頻分析的基礎(chǔ)上利用時頻旋轉(zhuǎn)特性，通 過判定某一旋轉(zhuǎn)角度下信號具有的最小時間帶寬積來達到最優(yōu)的時頻分辨率。2013年陳穎 頻等將基于分數(shù)階的自適應(yīng)最優(yōu)階Gabor變換應(yīng)用到地震信號頻譜分解中，取得了較好的 效果。2015年田琳等利用反卷積方法求取了分數(shù)階Gabor譜圖并應(yīng)用到地震儲層頻譜分解 中。該方法仍可以看作是在直接法計算功率譜的概念下做的進一步處理，即現(xiàn)階段的分數(shù) 域時頻分析技術(shù)尚未與功率譜估計技術(shù)結(jié)合。

【發(fā)明內(nèi)容】

[0007] 本發(fā)明是針對目前譜圖計算提出的新的計算方法，它結(jié)合了現(xiàn)代功率譜估計方法 并引入分數(shù)域最優(yōu)階思想提高局部功率譜的分辨率，以提高整體信號功率譜圖的時頻特 性。
[0008] 為了解決上述技術(shù)問題，達到上述目的，本發(fā)明采用如下技術(shù)方案：
[0009] 本發(fā)明提供了一種地震信號的分數(shù)域局部功率譜計算方法，其特征在于包括以下 步驟：
[0010] 步驟（1)、對輸入信號做P階分數(shù)傅里葉變換，求取分數(shù)域信號的時寬和頻寬及時 間帶寬積；
[0011] 步驟（2)、對變換階次P遍歷，即求得各個階次下分數(shù)域信號的時寬和頻寬及時 間帶寬積；
[0012] 步驟（3)、對各個階次下的時間帶寬積進行遍歷搜索，找到時間帶寬積最小值以及 對應(yīng)的階次、對應(yīng)的時寬、頻寬；
[0013] 步驟（4)、根據(jù)求得的時寬和頻寬構(gòu)造高斯窗函數(shù)，對高斯窗函數(shù)進行相應(yīng)階次的 分數(shù)階傅里葉變換求得最終的最優(yōu)窗函數(shù)，該窗函數(shù)一般是一個復(fù)函數(shù)；
[0014] 步驟（5)、將信號左端和右端進行周期性延拓，用最優(yōu)窗函數(shù)從信號第一個點為中 心開始截取固定長度，即信號本身的長度，截取的信號稱為原信號的局部，其功率譜反映的 是原信號的局部功率譜特征，對獲得局部應(yīng)用功率譜計算方法，如間接法、AR模型法、Burg 譜估計等，即可獲得有別于傳統(tǒng)譜圖的局部功率譜時頻分布。
[0015] 上述技術(shù)方案中，對輸入信號x(n)做p階分數(shù)傅里葉變換如下：
[0017] 其中p為變換階次，a為變換角度，u為分數(shù)域變量，e為常數(shù)，j為虛數(shù)單位，t為 時間變量，xp (u)為p階分數(shù)域的變換系數(shù)，當初始值p=0時，xp (u)即為x (n)，計算xp (u)的時寬Txp，對Xp(u)再做傅里葉變換得Xp+1 (v)，可計算出頻寬Bxp:

[0020] 其中I I *| |2為信號二范數(shù)的平方即信號能量，n u為信號P階分數(shù)域的時間中心， nv為信號p階分數(shù)域的頻率中心，即p+i階分數(shù)域的時間中心，求取階次p下的時間帶寬 積：
[0021] TBPp=Txp.Bxp (4)
[0022] 上述技術(shù)方案中，鑒于分數(shù)階傅里葉變換關(guān)于p = 2的對稱性，對p在0~2之間 遍歷，取步長為〇. 01，即P = P+0. 01，按步驟（1)計算新階次下的xp(u)、時寬Txp、Xp+1(v)、 頻寬Bx p、時間帶寬積TBPP，遍歷結(jié)束可獲得一個時寬向量Tx，頻寬向量Bx，時間帶寬積向 量 TBP。
[0023] 上述技術(shù)方案中，對時間帶寬積進行遍歷搜索，找到時間帶寬積最小值以及對應(yīng) 的階次及對應(yīng)的時寬、頻寬，
[0025] 即階次P(l滿足使信號在該階次下的時間帶寬積是最小的。
[0026] 上述技術(shù)方案中，由時寬TXd，頻寬8乂(|構(gòu)造分數(shù)階高斯窗函數(shù)gwin(u)(如式6)， 再對gwin(u)做階分數(shù)傅里葉變換，可得對應(yīng)的時域局部功率譜最優(yōu)窗函數(shù)g(n)，n為 時間序列：
[0029] 上述技術(shù)方案中，用窗函數(shù)g(n)截取輸入信號x(n)，獲得信號在m時刻的局部 y (n, m) = x(n) ? g(n-m)，m為窗的中心（其初始值為1)，y (n, m)是一個局部觀測信號，可 以結(jié)合功率譜估計的多種方法對y(n，m)求取功率譜，當窗函數(shù)的變量m逐漸滑動時即獲得 了信號x(n)的局部功率譜在不同時間下的時頻分布。
[0030] 上述技術(shù)方案中，在m時刻的局部功率譜可用直接法求取如下：
[0032] N為信號y (n，m)的長度，k為對應(yīng)的頻率點，STPSD (m，k)表示信號在m時刻在頻 率分量k處的功率譜密度。
[0033] 上述技術(shù)方案中，在m時刻的局部功率譜可用間接法求取如下：
[0034] ①計算y(n，m)的自相關(guān)函數(shù)r(l，m):
[0036] 其中1為時延；
[0037]②計算局部功率譜：
[0039] 上述技術(shù)方案中，在m時刻的局部功率譜可用AR模型法求取如下：
[0040] ①此方法是將局部觀測信號y (n，m)看作是一方差為的高斯白噪聲信號通過一 線性時不變最小相位系統(tǒng)產(chǎn)生的，若用M階自回歸模型來描述系統(tǒng)，系統(tǒng)表述為：
[0042] 其中 < 為m時刻系統(tǒng)待求的第k個模型系數(shù)，v (n，m)為m時刻觀測信號對應(yīng)的白 噪聲信號；
[0043] ②將式（11)兩邊同時乘以/(n_l，m)(y平移1個單位，并取共軛）并取期望可得 y的自相關(guān)函數(shù)構(gòu)成的方程組式（12)如下：
[0045] 式（12)寫成矩陣形式即轉(zhuǎn)化為Yule-Walker方程如式（13):
[0047]其中白噪聲方差可由下式（14)求得：
[0049] ③求解方程獲得系數(shù)<丨，則可求得局部功率譜式（15)如下：
[0051] 與現(xiàn)有技術(shù)相比，本發(fā)明具有如下優(yōu)勢：
[0052] 功率譜估計中現(xiàn)代譜估計較經(jīng)典譜估計有明顯較高的分辨率，本發(fā)明將其與分數(shù) 域時頻分析技術(shù)相結(jié)合，通過引入分數(shù)域最優(yōu)階思想，獲得較傳統(tǒng)窗函數(shù)性能更優(yōu)秀的復(fù) 窗函數(shù)，在截斷信號時能夠提高局部功率譜估計效果，從而最終提高地震信號功率譜圖的