姜傳金，陳樹民，劉 財，鹿 琪，馮智慧

1.大慶油田有限責(zé)任公司勘探開發(fā)研究院，黑龍江 大慶 163712

2.吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130061

0 引言

譜分解技術(shù)是源于B P Amoco（英國石油阿莫科）公司的一種地震解釋技術(shù)［1］。由于不同的地震頻率對各種地質(zhì)異常體的敏感度不同，所以在刻畫地質(zhì)異常體厚度變化及描述地質(zhì)異常體橫向不連續(xù)性等方面，譜分解技術(shù)已被證明是非常有效的方法。自20世紀(jì)90年代以來，譜分解技術(shù)得到了地質(zhì)、地球物理工作者的廣泛關(guān)注，而且不斷有新的地震譜分解技術(shù)出現(xiàn)［2］。地震信號統(tǒng)計特性隨時間而變，是一種典型的非平穩(wěn)信號，因此時頻分析方法能夠很好地表征其局部時間段頻域特性。目前從各類文獻(xiàn)資料中來看，基于時頻分析方法的譜分解技術(shù)大致分為兩大類：一類是線性時頻分析的方法，包括短時傅里葉變換［3］、小波變換［4－6］、S變換［7］、廣義S變換［8－10］、匹配追蹤［11］等；另一類是非線性時頻分析的方法，主要是指Wigner時頻分布［12－13］。理論上，上述方法各具特色和優(yōu)勢，但也都存在一定的局限。尤其Wigner時頻分布的缺點與其優(yōu)點一樣突出，即存在交叉干擾項問題［14－15］。

筆者利用Wigner高階譜具有高時頻聚集性的特點，提出了一種基于改進(jìn)核函數(shù)的 Wigner雙譜對角切片的譜分解技術(shù)。Wigner高階譜與其他二次型時頻分布一樣同樣存在交叉干擾項的問題，筆者主要圍繞抑制交叉干擾項的問題展開研究，采用模糊域核函數(shù)濾波法抑制Wigner雙譜對角切片的交叉項。分析指數(shù)型核函數(shù)與錐形核函數(shù)優(yōu)缺點后認(rèn)為：指數(shù)核函數(shù)不能抑制信號出現(xiàn)在橫軸和縱軸上的交叉項；與指數(shù)核函數(shù)相比，錐形核函數(shù)在縱軸旁邊存在一定的旁瓣，這樣在進(jìn)行模糊域濾波時，當(dāng)交叉項出現(xiàn)的位置與旁瓣有重合時，交叉項便不能得到很好的抑制。因此，筆者結(jié)合兩者的優(yōu)點給出了改進(jìn)的新核函數(shù)，以解決Wigner雙譜對角切片的模糊函數(shù)中心點校正問題。通過對新的核函數(shù)交叉項抑制能力的數(shù)值模擬，驗證新的核函數(shù)對交叉項有更好的抑制能力。最后對實際地震數(shù)據(jù)進(jìn)行分頻處理，指示油氣有利區(qū)的存在，驗證本文所提出方法的有效性。

1 方法原理

1.1 地震信號Wigner雙譜對角切片公式推導(dǎo)

Wigner高階矩譜與高階譜、高階矩、高階累計量一樣都屬于高階統(tǒng)計量。Wigner雙譜是Wigner高階矩譜的一種，按其變量參數(shù)所在的域劃分，它屬于時頻域的概念［16－17］。

設(shè)x（t）為任一零均值信號，由累積量理論可知，其二階累積量（自相關(guān)）與三階累積量定義為

式中：τ為信號二階累積量中的時間延遲變量；τ1，τ2為信號三階累積量中的時間延遲變量；E代表數(shù)學(xué)期望。

對其二階累積量（自相關(guān)）和三階累積量分別作一維傅里葉變換與二維傅里葉變換，便可以得到其功率譜與雙譜（自雙譜）［18－19］為

式中：f為信號功率譜中的頻率變量；f1，f2為信號雙譜中的頻率變量。

當(dāng)f1＝f2＝f時，稱W2x（t，f）為 Wigner雙譜對角切片［22］，其表達(dá)式為

1.2 核函數(shù)對交叉干擾項的抑制

在研究模糊函數(shù)時，人們發(fā)現(xiàn)一個重要的事實［23］：在模糊域，交叉項傾向于遠(yuǎn)離原點，而信號項則聚集在原點附近；因此，減小交叉項的一種很自然的方法是在模糊域用核函數(shù)對模糊函數(shù)進(jìn)行濾波，濾去交叉項，然后，再由模糊函數(shù)的傅里葉變換求相應(yīng)的時頻分布。這種抑制交叉項的方法稱為模糊域核函數(shù)濾波法。應(yīng)當(dāng)指出，交叉項的抑制與信號項的維持是一對矛盾。因為交叉項的減小必然會對信號項產(chǎn)生拉平的負(fù)面作用，從而減低時頻分布的聚集性。

例如信號z（t）＝x（t）＋y（t）。由褶積原理知，z（t）的Wigner分布與Wigner雙譜對角切片分別為

其中：2（Wxy（t，f））為Wigner分布的交叉項；與為 Wigner雙譜對角切片的交叉項。

信號z（t）的局部相關(guān)函數(shù)、局部三階累積量對角切片分別為

則信號z（t）的Wigner分布模糊函數(shù)定義為

同理，信號z（t）的 Wigner雙譜對角切片模糊函數(shù)定義為

其中：τ表示時間延遲；v表示頻偏（多分量信號間的頻率差值）。二者構(gòu)成時延－頻偏平面，即模糊域。利用模糊域核函數(shù)濾波法，建立 Wigner分布與Wigner雙譜對角切片計算公式：

在模糊域核函數(shù)濾波法中有2種比較常用的核函數(shù)。

第一種是指數(shù)核函數(shù)。

Choi與 Williams［24］在Cohen類分布中引入指數(shù)核函數(shù)（圖1），用該核與Wigner分布構(gòu)成的時頻分布稱為Choi－Williams分布，其模糊域表達(dá)式為

式中：τ為時間延遲；v為頻偏；α為控制核函數(shù)形狀的常數(shù)。

對指數(shù)核函數(shù)表達(dá)式及圖1進(jìn)行分析發(fā)現(xiàn)：

1）由式（16）易驗證φcw（0，0）＝1，φcw（0，v）＝1，φcw（τ，0）＝1。這表明，指數(shù)核函數(shù)對原點（0，0）以及橫軸（τ軸）和縱軸（v軸）上信號的模糊函數(shù)沒有任何影響。因此，若信號模糊函數(shù)的交叉項出現(xiàn)在橫軸和縱軸上，則它們將不能被抑制，從而時頻分布中相對應(yīng)的交叉項也不能被抑制。

2）當(dāng)τ≠0和v≠0時，信號的模糊函數(shù)在坐標(biāo)軸以外的交叉項都能夠得到一定程度的抑制，從而可以減少與這些模糊交叉項相對應(yīng)的時頻分布交叉項。

第二種是錐形核函數(shù)。

為了克服指數(shù)核函數(shù)不能抑制坐標(biāo)軸上存在交叉項的問題，Zhao等［25］提出了錐形核函數(shù)，其模糊域表達(dá)式為

與指數(shù)核函數(shù)對比發(fā)現(xiàn)：

1）由式（17）易驗證，錐形核函數(shù)除了對坐標(biāo)以外的信號模糊函數(shù)有影響外，對橫軸（τ軸）上的信號的模糊函數(shù)也有抑制作用。從而時頻分布中相對應(yīng)的交叉項也能被抑制。

2）錐形核函數(shù)也有其缺點。如圖2所示，與指數(shù)核函數(shù)相比，其模糊函數(shù)在縱軸旁邊存在一定的旁瓣。這樣在進(jìn)行模糊域濾波時，當(dāng)交叉項出現(xiàn)的位置與旁瓣有重合時，交叉項便不能得到很好的抑制。

圖1 指數(shù)核函數(shù)模糊函數(shù)圖Fig.1 Ambiguity function of the exponential kernel function

通過對指數(shù)核函數(shù)與錐形核函數(shù)優(yōu)缺點的分析，可以發(fā)現(xiàn)，指數(shù)核函數(shù)對坐標(biāo)軸外的交叉項有較好的抑制作用，而錐形核函數(shù)的突出特點是對時延軸上的交叉項有抑制作用。通過錐形核函數(shù)表達(dá)式可以發(fā)現(xiàn)，它由兩部分組成：一部分是基礎(chǔ)函數(shù)（圖3），它是一種類正弦函數(shù)；另一部分是一個指數(shù)函數(shù)（圖4）?？梢姡笖?shù)函數(shù)相當(dāng)于參數(shù)變量為時間延遲的一維濾波器，對基礎(chǔ)函數(shù)濾波之后就形成了錐形核函數(shù)。所以，完全可以利用該指數(shù)函數(shù)對指數(shù)核函數(shù)濾波，形成一個新核函數(shù)，使其既具有指數(shù)核函數(shù)對坐標(biāo)軸外交叉項的抑制能力，也具有對時延軸方向交叉項的抑制能力。該新核函數(shù)（圖5）表達(dá)式為

圖2 錐形核函數(shù)模糊函數(shù)Fig.2 Ambiguity function of the cone－shaped kernel function

圖3 錐形核基礎(chǔ)函數(shù)模糊函數(shù)Fig.3 Basic function in the cone－shaped kernel function

式中：τ為時間延遲；v為頻偏；α、β為控制核函數(shù)形狀的常數(shù)。

2 數(shù)值模擬

為了驗證新的核函數(shù)對交叉項的抑制作用，本節(jié)利用兩分量的模擬信號進(jìn)行了相應(yīng)的數(shù)值模擬（圖6）。

圖4 指數(shù)函數(shù)Fig.4 Exponential function in the cone－shaped kernel function

圖5 新核函數(shù)模糊函數(shù)Fig.5 Ambiguity function of the new kernel function

z（t）為主頻分別為30Hz、65Hz的零均值信號x（t）和y（t）的疊加信號（圖6a），即z（t）＝x（t）＋y（t）。疊加信號z（t）的 Wigner雙譜對角切片存在交叉項（圖6b）。如圖6d－f所示，加新核函數(shù)后的Wigner雙譜對角切片對疊加信號的交叉項有很好的抑制作用。需要強調(diào)的是，疊加信號基于Wigner雙譜對角切片的模糊函數(shù)的有效信號的中心點與模糊域的坐標(biāo)原點不一致（圖6c），因此需要通過坐標(biāo)移動，使其與模糊域坐標(biāo)原點一致（圖6d），從而使核函數(shù)有效抑制交叉項。之后需要將抑制交叉項后的有效信號重新回到原來的位置才能保證最后計算結(jié)果的準(zhǔn)確（圖6e－f）。

3 基于Wigner雙譜對角切片的譜分解技術(shù)流程

“基于Wigner雙譜對角切片的譜分解技術(shù)”的常規(guī)流程在實際應(yīng)用時有其局限性。即，方法流程中間生成Wigner雙譜對角切片三維數(shù)據(jù)體往往數(shù)據(jù)量很大，會使matlab程序報錯：內(nèi)存溢出?？梢酝ㄟ^由二維多道地震數(shù)據(jù)直接生成單頻剖面的方法避免上述問題。具體方法（圖7）為：首先初始化多個與地震數(shù)據(jù)行列數(shù)相等的零矩陣，即對第一道地震信號做Wigner雙譜對角切片后，得到一個二維時頻剖面，對此剖面做頻率切片，得到一個一維的單頻數(shù)據(jù)；然后將該數(shù)據(jù)存儲到零矩陣對應(yīng)的第一列中；依此方法類推，處理完整個地震記錄，便可以得到單頻剖面。基于Wigner分布的譜分解技術(shù)也可以采用上述方法。

4 實際應(yīng)用

為了進(jìn)一步驗證基于 Wigner雙譜對角切片分頻技術(shù)的實用性，筆者對圖8的實際數(shù)據(jù)按照頻率由低到高的原則進(jìn)行了分頻處理，并提取了各單頻率剖面，尋找頻率異常，驗證油氣顯示區(qū)域是否符合低頻能量較強、高頻能量較弱的現(xiàn)象（圖9）。如果油氣顯示區(qū)域滿足上述現(xiàn)象，就可以認(rèn)為利用筆者提出的譜分解技術(shù)成功地驗證了目的層段中油氣的存在。實際數(shù)據(jù)中的目的層為400～500ms，在橢圓標(biāo)記內(nèi)遇到油氣顯示。

從圖9中標(biāo)記內(nèi)的能量顯示可以看出：在20 Hz、31Hz這樣的低頻區(qū)域內(nèi)，油氣顯示區(qū)域內(nèi)有能量顯示；而隨著頻率的增加，在45Hz、50Hz這樣的高頻區(qū)域油氣顯示區(qū)域內(nèi)的能量團(tuán)消失。由此表明，筆者提出的基于改進(jìn)核函數(shù)的 Wigner分布譜分解技術(shù)成功地驗證了目的層段中油氣的存在。

5 結(jié)論與討論

1）筆者利用Wigner高階矩譜具有高時頻聚集性的特點，提出了一種基于改進(jìn)核函數(shù)的 Wigner雙譜對角切片的譜分解技術(shù)。與以往基于線性時頻分析理論的譜分解技術(shù)相比，不存在時窗的限制，因此具有更高的時頻分辨率。

2）在提出了改進(jìn)核函數(shù)的基礎(chǔ)上，解決了Wigner雙譜對角切片模糊函數(shù)中心點與核函數(shù)中心點不一致導(dǎo)致核函數(shù)抑制能力下降的問題。無論數(shù)值模擬與實際資料應(yīng)用都證明了其適用性。因此該譜分解技術(shù)可以作為油氣檢測的一種有效的輔助手段。

圖6 新核函數(shù)對模擬信號交叉項的抑制Fig.6 Suppression of new kernel function for the cross－terms with analog signals

3）在實際生產(chǎn)中，為了避免多解性，還需要與其他地球物理方法進(jìn)行綜合解釋才能確認(rèn)油氣的存在。此外受地層速度及薄層組合復(fù)雜性的影響，油氣檢測結(jié)果難免出現(xiàn)多解性，使得本次預(yù)測的精度受到限制。

圖7 改進(jìn)的分頻流程圖Fig.7 Modified flowchart of the frequency－divided technique based on the Wigner bispectrum diagonal slice

圖8 某一實際數(shù)據(jù)Fig.8 A field seismic data

圖9 實際數(shù)據(jù)的分頻顯示Fig.9 Single－frequency profiles of the seismic data shown in Fig.8

4）該項技術(shù)在實際應(yīng)用中以全區(qū)或多個構(gòu)造單元作為研究對象時，應(yīng)該以地質(zhì)先驗信息及經(jīng)保幅處理的常規(guī)屬性為基礎(chǔ)，從而可以使預(yù)測結(jié)果有據(jù)可依。如果簡單地單獨應(yīng)用難免以偏蓋全，摻雜過多的多解性。如果能夠分區(qū)域、分相帶展開研究，預(yù)測的結(jié)果會更加精確，會最大限度地降低多解性。

（References）：

［1］徐麗英，徐鳴潔，陳振巖.利用譜分解技術(shù)進(jìn)行薄儲層預(yù)測［J］.石油地球物理勘探，2006，41（3）：299－302.Xu Liying，Xu Mingjie，Chen Zhenyan.Using Spectrum Decomposition Technique for Prediction of Thin Reservoir［J］.Oil Geophysical Prospecting，2006，41（3）：299－302.

［2］張進(jìn)鐸，楊平，王云雷.地震信息的譜分解技術(shù)及其應(yīng)用［J］.勘探地球物理進(jìn)展，2006，29（4）：235－238.Zhang Jinduo，Yang Ping，Wang Yunlei.Spectral Decomposition of Seismic Data and Its Application in Oil and Gas Exploration［J］.Progress in Exploration Geophysics，2006，29（4）：235－238.

［3］Marfurt K J，Kirlin R L.Narrow－Band Spectral Analysis and Thin－Bed Tuning［J］.Geophysics，2001，66（4）：1274－1283.

［4］Sinha S，Routh P S，Anno P D，et al.Spectral Decomposition of Seismic Data with Continuous－Wavelet Transforms［J］.Geophysics，2005，70（6）：19－25.

［5］房文靜，范宜仁，李霞，等.基于測井?dāng)?shù)據(jù)小波變換的準(zhǔn)層序自動劃分［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，2007，37（4）：834－836.Fang Wenjing，F(xiàn)an Yiren，Li Xia，et al.Parasequence Automatical Partition Based on Wavelet Transform of Logging Data［J］.Journal of Jilin University：Earth Science Edition，2007，37（4）：834－836.

［6］馬志霞，孫贊東.Gabor－Morlet小波變換分頻技術(shù)及其在碳酸巖儲層預(yù)測中的應(yīng)用［J］.石油物探，2010，49（1）：42－45.Ma Zhixia，Sun Zandong.Spectral Decomposition Technique Based on Gabor－Morlet Wavelet Transform and Its Application to Carbonate Reservoir Prediction［J］.Geophysical Prospecting for Petroleum，2010，49（1）：42－45.

［7］楊海濤，朱仕軍，文中平，等.基于S－變換的譜分解效果分析［J］.石油天然氣學(xué)報，2009，31（5）：267－270.Yang Haitao，Zhu Shijun，Wen Zhongping，et al.Effect Analysis of Spectral Decomposition Based on S Transform［J］.Journal of Oil and Gas Technology，2009，31（5）：267－270.

［8］陳學(xué)華，賀振華，黃德濟(jì).基于廣義S變換的地震資料高效時頻譜分解［J］.石油地球物理勘探，2008，43（5）：530－534.Chen Xuehua，He Zhenhua，Huang Deji.High－Efficient Time－Frequency Spectrum Decomposition of Seismic Data Based on Generalized S Transform［J］.Oil Geophysical Prospecting，2008，43（5）：530－534.

［9］陳學(xué)華，賀振華，文曉濤，等.基于廣義S變換的裂縫分頻邊緣檢測方法［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，2011，41（5）：1605－1609.Chen Xuehua，He Zhenhua，Wen Xiaotao，et al.Fracture Multi－Frequency Edge Detection Based Generalized S Transform［J］.Journal of Jilin University：Earth Science Edition，2011，41（5）：1605－1609.

［10］王寶江，王大興，于波，等.利用廣義S變換技術(shù)預(yù)測砂巖儲層［J］.石油地球物理勘探，2006，41（4）：423－426.Wang Baojiang，Wang Daxing，Yu Bo，et al.Using Generalized S Transform Technique to Predict Sandstone Reservoir［J］.Oil Geophysical Prospecting，2006，41（4）：423－426.

［11］馮磊，姜在興.基于匹配追蹤的譜分解方法及其應(yīng)用［J］.勘探地球物理進(jìn)展，2009，32（1）：33－36.Feng Lei，Jiang Zaixing.Spectral Decomposition Based on Matching Pursuit and Its Application［J］.Progress in Exploration Geophysics，2009，32（1）：33－36.

［12］吳小羊，劉天佑.基于時頻重排的地震信號 Wigner－Ville分布時頻分析［J］.石油地球物理勘探，2009，44（2）：201－205.Wu Xiaoyang，Liu Tianyou.Time－Frequency Analysis on Wigner－Ville Distribution of Seismic Signal Based on Time－Frequency Rearrangement［J］.Oil Geophysical Prospecting，2009，44（2）：201－205.

［13］楊海濤，朱仕軍，文中平.基于Wigner－Ville的譜分解效果分析［J］.勘探地球物理進(jìn)展，2009，32（1）：37－39.Yang Haitao，Zhu Shijun，Wen Zhongping，et al.Effect Analysis of Spectral Decomposition Based on Wigner－Ville Distribution［J］.Progress in Exploration Geophysics，2009，32（1）：37－39.

［14］李耀東，馮濤，袁超偉.基于自項窗的WVD交叉項抑制技術(shù)［J］.無線電通信技術(shù)，2010，36（4）：39－42.Li Yaodong，F(xiàn)eng Tao，Yuan Chaowei.Auto－Term Window Method for Cross－Term Reduction of Wigner－Ville Distribution［J］.Radio Communications Technology，2010，36（4）：39－42.

［15］張賢達(dá).時間序列分析－高階統(tǒng)計量方法［M］.北京：清華大學(xué)出版社，1996.Zhang Xianda.Time Sequence Analysis［M］.Beijing：Tsinghua University Press，1996.

［16］張賢達(dá).現(xiàn)代信號處理［M］.2版.北京：清華大學(xué)出版社，2002.Zhang Xianda.Model Signal Processing［M］.2nd ed.Beijing：Tsinghua University Press，2002.

［17］Tugnait J K.On Time Delay Estimation with Unknown Spatially Correlated Gaussian Noise Using Fourth－Ordercumulants and Cross Cumulants［J］.IEEE Trans Signal Processing，1991，39（6）：1259－1267.

［18］Mendel J M.Tutorial on Higher－Order Statistics Insignal Orocessing and System Theory［J］.Theoretical Results and Some Applications Proc IEEE，1991，79（3）：278－305.

［19］Nikias L C，Mendel M J.Signal Processing with Higher－Order Spectra［J］.IEEE Trans Signal Processing，1993，10（3）：10－37.

［20］Gerr N L.Introducing a Third－Order Wigner Distribution［J］.Proc IEEE，1988，3：290－292.

［21］Lee S K，White P R.Higer－Order Time－Fpequecncy Analysis and Its Application to Fault Detection in Rotating Machinery［J］.Mechanical Systems and Signal Processing，1997，11（4）：637－650.

［22］劉波.MATLAB信號處理［M］.北京：電子工業(yè)出版社，2006：330.Liu Bo.MATLAB Signal Processing［M］.Beijing：Publishing House of Electronics Industry，2006：330.

［23］王軍，李亞安.多分量信號時頻分析交叉項抑制研究［J］.探測與控制學(xué)報，2004，26（4）：13－17.Wang Jun，Li Yaan.The Research of Cross－Component Suppression of Multi－Component Signals Using Time Frequency Analysis［J］.Journal of Detection ＆Control，2004，26（4）：13－17.

［24］秦太龍，程珩，薛松.抑制Wigner－Ville分布交叉項方法的比較［J］.太原理工大學(xué)學(xué)報，2008，36（9）：548－550.Qin Tailong，Cheng Heng，Xue Song.The Comparison of Methods That Inhibit Wigner－Ville Distribution’s Cross－Term［J］.Journal of Taiyuan University of Technology，2008，36（9）：548－550.

［25］Zhao Yunxin，Atlas E L，Marks J R.The Use of Cone－Shaped Kernels for Generalized Time－Frequency Distributions of Nonstationary Signals［J］.IEEE Transaction on Acoustics，Speech，and Signal Processing，1990，38（7）：1084－1091.