最小平方約束反演譜分析方法的應(yīng)用效果分析

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

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

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

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

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

1 基本原理

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

(1)

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

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

(2)

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

(3)

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

(4)

這里，α的取值為

(5)

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

(6)

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

(7)

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

2 模型試算分析

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

2.1 單子波試算分析

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

圖1 子波頻譜分析a 子波及計(jì)算時(shí)窗; b 不同方法頻譜對(duì)比

2.2 模擬道試算分析

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

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

圖2 合成余弦信號(hào)時(shí)頻分析a 合成余弦信號(hào); b STFT時(shí)頻切片; c S變換時(shí)頻切片; d CLSSA時(shí)頻切片

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

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

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

3.1 單頻剖面分析

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

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

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

3.2 沿層屬性分析

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

參 考 文 獻(xiàn)

[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í)變分頻技術(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ì),等.時(shí)頻域油氣儲(chǔ)層低頻陰影檢測[J].地球物理學(xué)報(bào),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ù)在含氣儲(chǔ)層預(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ù)在油氣識(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