基于Wigner雙譜對角切片的譜分解技術在油氣檢測中的應用

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

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

2.吉林大學地球探測科學與技術學院，長春 130061

0 引言

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

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

1 方法原理

1.1 地震信號Wigner雙譜對角切片公式推導

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

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

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

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

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

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

1.2 核函數對交叉干擾項的抑制

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

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

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

信號z（t）的局部相關函數、局部三階累積量對角切片分別為

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

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

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

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

第一種是指數核函數。

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

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

對指數核函數表達式及圖1進行分析發現：

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

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

第二種是錐形核函數。

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

與指數核函數對比發現：

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

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

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

通過對指數核函數與錐形核函數優缺點的分析，可以發現，指數核函數對坐標軸外的交叉項有較好的抑制作用，而錐形核函數的突出特點是對時延軸上的交叉項有抑制作用。通過錐形核函數表達式可以發現，它由兩部分組成：一部分是基礎函數（圖3），它是一種類正弦函數；另一部分是一個指數函數（圖4）。可見，指數函數相當于參數變量為時間延遲的一維濾波器，對基礎函數濾波之后就形成了錐形核函數。所以，完全可以利用該指數函數對指數核函數濾波，形成一個新核函數，使其既具有指數核函數對坐標軸外交叉項的抑制能力，也具有對時延軸方向交叉項的抑制能力。該新核函數（圖5）表達式為

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

圖3 錐形核基礎函數模糊函數Fig.3 Basic function in the cone－shaped kernel function

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

2 數值模擬

為了驗證新的核函數對交叉項的抑制作用，本節利用兩分量的模擬信號進行了相應的數值模擬（圖6）。

圖4 指數函數Fig.4 Exponential function in the cone－shaped kernel function

圖5 新核函數模糊函數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所示，加新核函數后的Wigner雙譜對角切片對疊加信號的交叉項有很好的抑制作用。需要強調的是，疊加信號基于Wigner雙譜對角切片的模糊函數的有效信號的中心點與模糊域的坐標原點不一致（圖6c），因此需要通過坐標移動，使其與模糊域坐標原點一致（圖6d），從而使核函數有效抑制交叉項。之后需要將抑制交叉項后的有效信號重新回到原來的位置才能保證最后計算結果的準確（圖6e－f）。

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

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

4 實際應用

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

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

5 結論與討論

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

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

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

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

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

圖8 某一實際數據Fig.8 A field seismic data

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

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

（References）：

［1］徐麗英，徐鳴潔，陳振巖.利用譜分解技術進行薄儲層預測［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］.勘探地球物理進展，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］房文靜，范宜仁，李霞，等.基于測井數據小波變換的準層序自動劃分［J］.吉林大學學報：地球科學版，2007，37（4）：834－836.Fang Wenjing，Fan 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小波變換分頻技術及其在碳酸巖儲層預測中的應用［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］.石油天然氣學報，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］陳學華，賀振華，黃德濟.基于廣義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］陳學華，賀振華，文曉濤，等.基于廣義S變換的裂縫分頻邊緣檢測方法［J］.吉林大學學報：地球科學版，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變換技術預測砂巖儲層［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］馮磊，姜在興.基于匹配追蹤的譜分解方法及其應用［J］.勘探地球物理進展，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］.勘探地球物理進展，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交叉項抑制技術［J］.無線電通信技術，2010，36（4）：39－42.Li Yaodong，Feng 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］張賢達.時間序列分析－高階統計量方法［M］.北京：清華大學出版社，1996.Zhang Xianda.Time Sequence Analysis［M］.Beijing：Tsinghua University Press，1996.

［16］張賢達.現代信號處理［M］.2版.北京：清華大學出版社，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］.北京：電子工業出版社，2006：330.Liu Bo.MATLAB Signal Processing［M］.Beijing：Publishing House of Electronics Industry，2006：330.

［23］王軍，李亞安.多分量信號時頻分析交叉項抑制研究［J］.探測與控制學報，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］.太原理工大學學報，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.