基于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.