劉炳楊，李勝軍，高建虎，劉軍迎

(1.中國石油天然氣股份有限公司勘探開發(fā)研究院西北分院，甘肅蘭州730020;2.中國石油天然氣集團公司油藏描述重點實驗室，甘肅蘭州730020)

解析信號的時頻特性分析研究最早可追溯到1946年Gabor[1]提出的Gabor展開。現(xiàn)今常用的時頻分析方法大體可以分為3類。第1類是線性變換方法，包括短時傅里葉變換[2]、小波變換[3]、S變換[4]和廣義S變換[5-6]等，這類方法采用線性變換形式描述信號的頻譜隨時間的變化規(guī)律，算法較為簡單且計算效率較高，但也存在一定問題，如短時傅里葉變換方法用固定的時窗取一段信號進行分析，無法調(diào)節(jié)時間和頻率分辨率；小波變換方法用小波庫對信號進行分解，具有多分辨率特點，但其時間尺度的概念與頻率關(guān)系不直接；S變換以及廣義S變換引入了滑動高斯窗函數(shù)，其時窗大小隨頻率而變，能夠在低頻端獲得較高的頻率分辨率而在高頻端獲得較高的時間分辨率，但其窗函數(shù)形式固定，一定程度上限制了它的發(fā)展。第2類常用的時頻分析方法是二次型時頻變換法，如Wigner分布。這種方法具有較高的時頻分辨率，但美中不足是信號的時頻分量間存在交叉項，干擾時頻圖中的有效信號。第3類時頻分析方法是近年來比較熱門的貪婪算法，包括匹配追蹤(MP)[7]以及Hilbert-Huang變換(HHT)[8]方法，這兩種方法都是通過搜索從信號中減去一個最優(yōu)的基信號，不斷循環(huán)直到滿足停止條件，這類方法能夠獲取稀疏的高分辨率時頻譜，但計算量較大，且都是在單個循環(huán)中尋求最優(yōu)解而不是全局最優(yōu)解?？梢哉f，每一種方法都有其優(yōu)點和缺點，難以滿足所有需求。近年來，隨著分頻地震解釋[9]、分頻AVO[10-11]以及吸收衰減[12-13]等其它基于地震數(shù)據(jù)體頻率變化的儲層分析檢測方法的發(fā)展，對高精度時頻分析算法的需求越來越迫切。

傳統(tǒng)反演問題的求解方法是構(gòu)造誤差函數(shù)的某一范數(shù)(通常是2范數(shù))并使其為最小，在這一條件下對未知量進行迭代求解(目前常用Tikhonov回歸方法)。此外，為了使求解過程更加準(zhǔn)確、穩(wěn)定，求解過程中一般還要引入其它約束條件，如L1范數(shù)、正則約束等。

最近，發(fā)展了一類新興的、基于反演的譜分析方法。這類方法將譜分解問題描述成一個線性反問題，并通過加入其它約束條件來對頻譜進行迭代求解，可以獲得高時頻分辨率的譜。韓利[14]介紹了一種以L1范數(shù)為約束條件的稀疏反演譜分解方法，能夠獲得信號的稀疏解，得到高分辨率頻譜，但其穩(wěn)定性方面則略有欠缺。Puryear等[15]借鑒傳統(tǒng)反演的思想，提出了一種最小平方約束譜分析(CLSSA)方法，這種方法實際上可以說是Portniaguine等[16]所提出的反演譜分解(Inverse Spectra Decomposition)方法中的一種。我們對最小平方約束反演譜分析(Constrained Least-Square Spectra Analysis，CLSSA)[15]這一高精度時頻分析方法的應(yīng)用效果進行了分析研究，發(fā)現(xiàn)該方法能夠準(zhǔn)確地拾取小時窗范圍內(nèi)地震信號的頻率。應(yīng)用最小平方約束反演譜分析方法對國內(nèi)某油田實際地震數(shù)據(jù)進行頻譜分析的結(jié)果表明，其應(yīng)用效果較傳統(tǒng)分頻方法有明顯改善，更為準(zhǔn)確、清晰地刻畫了研究區(qū)內(nèi)目標(biāo)層位發(fā)育的辮狀河道。

1 基本原理

假設(shè)時窗地震信號為d，對應(yīng)的頻譜為m，F(xiàn)為傅氏反變換矩陣，可得到形如(1)式的正問題：

(1)

式中：d可以是實地震信號，也可以是由原時窗信號及其Hilbert變換得到的虛地震信號共同形成的復(fù)地震信號。試驗證明，用復(fù)信號進行計算效果更好一些。

為了保證解的唯一性以及穩(wěn)定性，Puryear等[15]引入了模型權(quán)重Wm以及數(shù)據(jù)權(quán)重Wd。Wm初始為單位矩陣，隨迭代求解過程而改變；而對于窗函數(shù)Wd，我們選用的是Hanning窗。

(2)

可以推導(dǎo)得到形如(3)式的以L2范數(shù)進行約束的成本函數(shù)

(3)

(3)式的最小范數(shù)解為

(4)

這里，α的取值為

(5)

式中：αF為一常系數(shù)，一般為0.010～0.001，用以保證算法的穩(wěn)定性。那么很容易得到所求頻譜為

(6)

為了得到更加收斂的頻譜，可以進行多次迭代求解，只需更新模型權(quán)重Wm即可，其更新公式為

(7)

采用本方法對頻譜進行求解，可以直接利用(4)式和(6)式一次計算求得，還可以引入(7)式進行迭代求解。迭代次數(shù)越多，獲取的頻譜越發(fā)的收斂，甚至可以獲得脈沖狀的頻譜。但這樣也會帶來更大的計算量，需要依實際需要進行均衡選擇。原則上，當(dāng)選取信號的時窗長度較小時，可以適當(dāng)增加迭代次數(shù)，以獲取更高的頻率分辨率。

2 模型試算分析

為了分析CLSSA方法進行譜分解的效果，分別對單個子波、混疊余弦信號進行了試算，并同短時傅里葉變換(STFT)以及S變換方法進行了對比。

2.1 單子波試算分析

取一30Hz主頻的Ricker理論子波進行試算，并與傅里葉方法進行對比(子波信號采樣間隔為1ms)。以子波峰值點為中心，用Hanning窗分別取40，20ms時窗數(shù)據(jù)(圖1a)進行頻譜試算，用來比較不同方法在小時窗情況下的頻率分辨率，結(jié)果如圖1b所示。其中，用CLSSA方法進行頻譜分析時均采用了一次迭代。為了便于對比，對各頻譜做了歸一化處理，以便疊合在一起顯示。圖1b 中黑線是對整個子波采用傅氏變換得到的頻譜，假定為真實頻譜；綠線為STFT方法的分析結(jié)果，可明顯看到其頻率分辨率較差，40ms時窗情況下(綠色實線)勉強可以識別出子波主頻，但當(dāng)時窗長度縮小為20ms時(綠色虛線)，主頻完全無法識別；藍線為僅采用原始信號的CLSSA方法分析結(jié)果，可見40ms時窗情況下(藍色實線)所得頻譜的頻率分辨率相比于STFT明顯得到改善，但當(dāng)時窗縮小為20ms時(藍色虛線)，依然無法正確識別子波主頻；紅線為對原始信號進行Hilbert變換得到虛部道，然后與原信號一起合成復(fù)數(shù)道，再采用CLSSA方法對復(fù)數(shù)道進行譜分析的結(jié)果，可以看到40ms時窗情況下(紅色實線)可以得到高頻率分辨率的頻譜，當(dāng)時窗縮小為20ms時，依舊能夠準(zhǔn)確地識別子波的主頻，且頻率分辨率較40ms時窗的STFT計算結(jié)果高得多?？梢?，相較于STFT方法，采用復(fù)數(shù)道的CLSSA方法可以用較小的分析時窗得到更為精確的分頻結(jié)果，在提高時間分辨率的同時還能改善頻率分辨率。

圖1 子波頻譜分析a 子波及計算時窗; b 不同方法頻譜對比

2.2 模擬道試算分析

為了進一步對CLSSA方法進行評價分析，設(shè)計了一個用不同頻率余弦信號混疊而成的模擬記錄道。模擬記錄道長300ms，1ms間隔采樣，其中，0～100ms時窗內(nèi)為20Hz余弦信號；100～200ms時窗內(nèi)為50Hz的余弦信號；200～300ms時窗內(nèi)為(20+50)Hz混疊的余弦信號(圖2a)。分別采用STFT，S變換和CLSSA方法進行時頻分析，結(jié)果分別如圖2b，圖2c和圖2d所示(其中STFT方法和CLSSA方法采用的分析時窗為40ms)。對比可以看出，STFT方法頻率分辨率最差，尤其是對于兩頻率混合的余弦信號，幾乎無法分辨；S變換方法獲得的時頻譜相較于STFT方法有所改善，基本可以分辨混疊信號，但效果不太理想；而用CLSSA方法經(jīng)過15次迭代后得到的時頻譜效果好得多，可清晰、準(zhǔn)確地識別出包括混合信號在內(nèi)的各時窗內(nèi)信號的頻率，分頻效果理想。

地震記錄通常用地震子波與反射系數(shù)的褶積來描述，為了分析CLSSA方法對地震信號的時頻分辨效果，對一個由單子波構(gòu)成的信號道進行了試算分析。子波選用的是30Hz主頻的Ricker子波，記錄道時長200ms，1ms間隔采樣，如圖3a所示。分別用STFT，S變換和CLSSA方法進行時頻分析(STFT方法和CLSSA方法的分析時窗長度為20ms)，結(jié)果分別如圖3b，圖3c和圖3d所示。分析可見，STFT方法得到的分頻結(jié)果時間分辨率還可以，但頻率分辨率較差；S變換方法得到的分頻結(jié)果的頻率分辨率較STFT方法有所提高，但時間分辨率變得很差，尤其是低頻段；而采用CLSSA方法所得的分頻結(jié)果，無論是時間分辨率還是頻率分辨率要遠(yuǎn)遠(yuǎn)優(yōu)于STFT方法和S變換方法，分頻效果較為理想。

圖2 合成余弦信號時頻分析a 合成余弦信號; b STFT時頻切片; c S變換時頻切片; d CLSSA時頻切片

圖3 單子波道時頻分析a 單子波道; b STFT方法時頻切片; c S變換時頻切片; d CLSSA時頻切片

3 實際資料應(yīng)用效果分析

通過模型數(shù)據(jù)分析結(jié)果可知，CLSSA方法具有十分理想的時頻分辨率，因此該方法獲得了廣泛的應(yīng)用[17-18]。我們分別采用STFT，S變換和CLSSA方法對國內(nèi)某油田實際地震資料進行分頻分析，提取最大振幅屬性以及主頻屬性。

3.1 單頻剖面分析

從研究區(qū)內(nèi)切出201線×201道的地震數(shù)據(jù)(道間距為20m)進行研究(數(shù)據(jù)采樣間隔為2ms)。圖4a為原始剖面，圖中黃線標(biāo)示的是本次研究的目標(biāo)層位。通過頻譜分析，目標(biāo)層段的主頻約為15Hz。因此，分別采用STFT，S變換和CLSSA方法抽取15Hz單頻數(shù)據(jù)體作對比分析，結(jié)果分別如圖4b，圖4c和圖4d所示(STFT和CLSSA方法的分析時窗為20ms)。對比可見，STFT方法和CLSSA方法提取的單頻切片都具有較好的時間分辨率，相比之下S變換方法提取的單頻剖面的時間分辨率顯得很差；STFT方法抽取的單頻剖面，其能量軸的強弱關(guān)系與原始地震記錄剖面完全一致，而CLSSA方法抽取的單頻剖面則更加突出了高能量異常體。

另外，圖4a箭頭所示為河道發(fā)育處，通過觀察我們發(fā)現(xiàn)，沿目標(biāo)層位發(fā)育的河道沒有表現(xiàn)出明顯的強反射軸，反而在河道邊緣表現(xiàn)出弱反射特征。

圖4 實際地震數(shù)據(jù)單頻剖面a 原始地震剖面; b STFT方法15Hz振幅剖面; c S變換方法15Hz振幅剖面; d CLSSA方法15Hz振幅剖面

3.2 沿層屬性分析

CLSSA方法的主要優(yōu)勢在于能用較小的時窗準(zhǔn)確識別地震道某時刻的頻率，用來識別地震數(shù)據(jù)體空間范圍內(nèi)的頻率變化有望取得很好的效果。

首先沿圖4所示研究區(qū)目的層開40ms時窗，分別采用STFT，S變換和CLSSA方法對時窗數(shù)據(jù)進行計算得到頻譜，提取最大振幅值及其對應(yīng)的峰值頻率(主頻)，結(jié)果如圖5和圖6所示。

圖5 沿層最大振幅屬性切片(40ms時窗)a 原始數(shù)據(jù)切片; b STFT方法最大振幅切片; c S變換方法最大振幅切片; d CLSSA方法最大振幅切片

圖6 沿層主頻屬性切片(40ms時窗)a STFT方法主頻切片; b S變換方法主頻切片; c CLSSA方法主頻切片

研究區(qū)內(nèi)發(fā)育河道砂體，從圖5a原始地震數(shù)據(jù)沿層切片中隱約可見。圖5b為40ms時窗長度下，采用STFT方法提取的最大振幅值沿層切片，河道隱約可見，但展布特征不清晰；圖5c為采用S變換方法提取的最大振幅屬性，可見，對河道的識別效果較STFT方法要好得多；圖5d為采用CLSSA方法提取的沿層最大振幅切片，同樣可以清晰地展現(xiàn)研究區(qū)內(nèi)的河道特征，并且同S變換提取的切片相比，CLSSA方法提取的切片中河道的低振幅特性更加突顯。

當(dāng)河道砂體中飽含流體時，會導(dǎo)致地震波主頻向低頻端衰減，特別是含氣情況下，這種衰減現(xiàn)象變得尤為明顯。這是利用地震數(shù)據(jù)頻率變化屬性進行儲層及含油氣性檢測的理論基礎(chǔ)。我們提取了主頻屬性對研究區(qū)進行分析。圖6為圖5中用不同方法求取的最大振幅值處所對應(yīng)的峰值頻率，即主頻切片。從圖6中可見，河道發(fā)育處地震波主頻明顯降低。圖6a為STFT方法提取的主頻切片，河道特征并不明顯，很大一部分淹沒在背景信息里；圖6b是S變換方法提取的主頻切片，相較于STFT方法的主頻切片，下面的“人”字形河道以及左上部的曲流河特征都得到了較為清晰的展現(xiàn)，但同時也存在局部特征不連續(xù)現(xiàn)象；圖6c為CLSSA方法提取的主頻切片，河道特征更加清晰連續(xù)，尤其是紅框所示區(qū)域，相較于前兩種方法，雕刻效果有較大改善。

進一步縮小時窗，對比STFT方法和CLSSA方法的效果。圖7為20ms時窗分析結(jié)果。STFT方法的最大振幅切片(圖7a)中，僅隱約可見河道，而主頻切片(圖7b)出現(xiàn)大片零頻率區(qū)域，說明在20ms時窗下，STFT方法已經(jīng)無法準(zhǔn)確地進行頻率識別。而相比之下，CLSSA方法拾取的最大振幅切片(圖7c)中河道特征依然很顯著，同40ms時窗(圖5d)相比沒有大的變化；主頻切片(圖7d)中，對河道的刻畫效果比圖6c更好，分辨率有所提高，雕刻更加細(xì)微。

將分析時窗長度縮小到10ms，STFT方法以及CLSSA方法所得分頻結(jié)果如圖8所示?？梢钥吹剑琒TFT方法已經(jīng)完全不能進行頻率識別，主頻切片(圖8b)中絕大部分區(qū)域都成了零頻率，而少數(shù)非零頻率也拾取得完全不準(zhǔn)確，看起來更像是奇異值點。而CLSSA方法提取的最大振幅(圖8c) 和主頻(圖8d)依然能夠準(zhǔn)確拾取頻率特征，清晰刻畫了河道特征。

實際資料應(yīng)用結(jié)果表明，與STFT和S變換相比，CLSSA方法能夠更為準(zhǔn)確地拾取地震記錄中的頻率變化，對研究區(qū)內(nèi)目標(biāo)層位引起局部主頻下降的河道砂體進行精細(xì)雕刻。此外，后續(xù)減小時窗后STFT和CLSSA方法的應(yīng)用效果對比也表明，CLSSA方法對時窗長度要求較低，能夠用一個很小的時窗準(zhǔn)確拾取信號頻率，也就是說，CLSSA方法在用小時窗保證時間分辨率的同時，還能夠獲得很好的頻率分辨率。

圖7 20 ms時窗沿層屬性a STFT最大振幅切片; b STFT主頻切片; c CLSSA最大振幅切片; d CLSSA主頻切片

圖8 10ms時窗沿層屬性a STFT最大振幅切片; b STFT主頻切片; c CLSSA最大振幅切片; d CLSSA主頻切片

4 結(jié)束語

最小平方約束譜分析(CLSSA)方法采用反演的思路對信號頻譜進行迭代求解，能夠獲得信號的更為精準(zhǔn)的頻譜。模型數(shù)據(jù)試算結(jié)果表明，CLSSA方法相較于傳統(tǒng)時頻分析方法，能夠同時獲得較高的時間分辨率和頻率分辨率。實際資料的應(yīng)用效果也表明，CLSSA方法用很小的時窗就能夠十分準(zhǔn)確地進行主頻識別，能夠?qū)ρ芯繀^(qū)內(nèi)引起頻率變化的地質(zhì)異常體(河道)進行準(zhǔn)確識別和細(xì)致刻畫。鑒于該方法對分析時窗長度要求較低，建議應(yīng)用時選取較小的時窗進行分析，這樣在保證高頻率分辨率的同時，還能夠獲得更高的時間分辨率。

參 考 文 獻

[1] Gabor D.Theory of communication,part 1:the analysis of information[J].Electrical Engineers-Part III:Radio and Communication Engineering,1946,93(26):429-441

[2] Nawab S H,Quatieri T F.Short-time fourier transform[J].Advanced Topics in Signal Processing,1988,6(2):289-337

[3] Daubechies I.The wavelet transform,time-frequency localization and signal analysis[J].Information Theory,IEEE Transactions on,1990,36(5):961-1005

[4] Stockwell R G,Mansinha L,Lowe R P.Localization of the complex spectrum:the S transform[J].Signal Processing,IEEE Transactions on,1996,44(4):998-1001

[5] Adams M D,Kossentini F,Ward R K.Generalized S transform[J].Signal Processing,IEEE Transactions on,2002,50(11):2831-2842

[6] 王長江,楊培杰,羅紅梅,等.基于廣義 S 變換的時變分頻技術(shù)[J].石油物探,2013,52(5):489-494

Wang C J,Yang P J,Luo H M,et al.Time variable frequency division based on generalized S transform[J].Geophysical Prospecting for Petroleum,2013,52(5):489-494

[7] Wang Y H.Seismic time-frequency spectral decomposition by matching pursuit[J].Geophysics,2006,72(1):V13-V20

[8] Huang N E,Shen S S.Hilbert-Huang transform and its applications[M].Singapore:World Scientific Publishing Company,2005:1-311

[9] 陳學(xué)華,賀振華,黃德濟,等.時頻域油氣儲層低頻陰影檢測[J].地球物理學(xué)報,2009,52(1):215-221

Chen X H,He Z H,Huang D J,et a1.Low frequency shadow detection of gas reservoirs in time-frequency domain[J].Chinese Journal of Geophysics,2009,52(1):215-221

[10] 寧媛麗.頻散介質(zhì)基于反演譜分解的AVO方法研究[D].長春：吉林大學(xué),2012

Ning Y L.Study on AVO dispersive medium using inverse spectral decomposition[D].Changchun:Jilin University,2012

[11] 何兵紅,吳國忱,郭念民.基于疊前地震道集的 FVO 分析方法[J].石油地球物理勘探,2013，48(1):94-102

He B H,Wu G C，Guo N M.FVO analysis using pre-stack seismic data[J].Oil Geophysical Prospecting,2013,48(1):94-102

[12] 王小杰,印興耀,吳國忱.基于 S 變換的吸收衰減技術(shù)在含氣儲層預(yù)測中的應(yīng)用研究[J].石油物探,2012,51(1):37-42

Wang X J,Yin X Y,Wu G C.The application of an S transform based absorption and attenuation technique for prediction of gas-bearing reservoir[J].Geophysical Prospecting for Petroleum,2012,51(1):37-42

[13] 刁瑞,李振春,韓文功,等.基于廣義 S 變換的吸收衰減分析技術(shù)在油氣識別中的應(yīng)用[J].石油物探,2011,50(3):260-265

Diao R,Li Z C,Han W G,et al.Application of absorption and attenuation analysis technique based on generalized S transform for hydrocarbon identification[J].Geophysical Prospecting for Petroleum,2011,50(3):260-265

[14] 韓利.高分辨率全譜分解方法研究[D].長春:吉林大學(xué),2013

Han L.Research on the methods of high-resolution full spectrum decomposition[D].Changchun:Jilin University,2013

[15] Puryear Charles I,Oleg N Portniaguine,Carlos M Cobos,et al.Constrained least-squares spectral analysis:application to seismic data[J].Geophysics,2012,77(5):V143-V167

[16] Portniaguine O,Castagna J.Inverse spectral decomposition[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,1786-1789

[17] Qi J,Castagna J.Application of a PCA fault-attribute and spectral decomposition in Barnett Shale fault detection[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013,1421-1425

[18] Oyem A,Castagna J.Layer thickness estimation from the frequency spectrum of seismic reflection data[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013,1451-1455