基于分数阶Wigner-Ville分布的地震信号谱分解

基于分数阶Wigner-Ville分布的地震信号谱分解
张晓燕;彭真明;张萍;何艳敏;田琳
【摘要】本文将经典Wigner-Ville分布(WVD)拓展到分数阶Wigner-Ville分布(FrWVD),并将其引入到复杂地震信号的时频分析和谱分解中.利用分数阶傅里叶变换的旋转特性和广义时频带宽积理论可以获得比经典WVD更高的时频分辨率,提出的最优分数阶伪WVD能更好地解决交叉项抑制问题.理论信号仿真及实际地震资料的处理结果验证了该方法的有效性.同时,分数域频谱成像的引入为油气储层流体识别提供了一条新的途径.
【期刊名称】《石油地球物理勘探》
【年(卷),期】2014(049)005
【总页数】7页(P839-845)
【关键词】分数阶Wigner-Ville分布;最优分数阶伪Wigner-Ville分布;交叉项抑制;时频谱分解
【作者】张晓燕;彭真明;张萍;何艳敏;田琳
【作者单位】电子科技大学光电信息学院,四川成都610054;电子科技大学光电信息学院,四川成都610054;电子科技大学光电信息学院,四川成都610054;电子科技大学光电信息学院,四川成都610054;电子科技大学光电信息学院,四川成都610054
【正文语种】中文
【中图分类】P631
1 引言
地震信号是一种典型的非平稳信号，时频分析方法则是处理非平稳信号的有力工具［1］。

传统的时频分析方法主要包括两类，即以短时傅里叶变换（Short Time Fourier Transform，STFT）、Gabor变换［2，3］、小波变换、S变换为代表
的线性时频分析方法和以 Wigner－Ville分布（Wigner－Ville Distribution，WVD）、Cohen类时频分布为代表的二次型时频分布［4，5］。

线性时频分析的时频聚集性往往不高，分辨率较低。

而WVD具有非常高的时频分辨率，但由于它的非线性性质，处理多分量信号必然存在交叉项干扰，因此在很大程度上限制了其应用范围。

Cohen类时频分布即是针对WVD的交叉项问题提出的一类核函数去
除交叉项方法，主要分为固定核函数方法和自适应核函数方法。

由于固定核函数方法的核函数形状固定，只能较好地去除某一类信号的交叉项，不具有普适性。

自适应核函数方法，如径向高斯核分布、自适应最优核分布等，虽能自适应地调整核函数的形状，但其迭代过程较复杂，运算效率太低，因此在大数据处理中也受到了限制［4～6］。

分数域时频分析技术是近年来出现的一种新的时频分析技术［7～13］，最基本的方法就是分数阶傅里叶变换（Fractional Fourier Transform，FrFT）。

分数域时频分析可在传统时频分析的基础上利用Fr－FT特有的旋转特性提高时频分辨率，能更好地检测到信号的突变及异常信息，尤其适合复杂油气藏储层预测及流体性质判别。

分数阶 Wigner－Ville分布（Fractional Wigner－Ville Distribution，FrWVD）最初由Dragonman［9］提出，在光学领域既可用于远场衍射也可用于近场衍射。

但因交叉项的存在，FrWVD在信号处理领域并未得到广泛应用，只被简单地用于检测线性调频信号［10］，尤其在地震信号处理中还是一个空白。

此外，现有
FrWVD的定义形式有所差异，并未形成统一的理论体系。

分析前人的研究成果，笔者认为Dragonman提出的FrWVD更适合于地震信号的处理。

在此基础上，文中结合分数阶伪WVD与分数阶STFT之间的关系推导了最优分数阶伪WVD的计算公式。

理论信号仿真及实际地震资料处理结果表明，最优分数阶伪WVD能在保证高分辨率的前提下较好地去除交叉项，且不会增加额外的计算成本。

2 FrWVD理论
2.1 基本原理
关于FrWVD，目前还没有形成统一理论，大体上有以下3种定义。

2.1.1 定义一
由 Dragoman［9］提出，表达式为
式中：α＝pπ／2为信号时频图顺时针旋转的角度，其中p为FrFT的阶次；xp （u）为信号x（t）在阶次p下的FrFT；“＊”表示取共轭；Wα（u，v）为最终得到的FrWVD。

以上定义形式可理解为先对信号做某一阶次的FrFT，再对变换后的信号进行WVD计算，最终得到信号的FrWVD，且算法复杂度为O（N2logN）（N为信号点数）。

该定义的提出，使得原本只能用于远场衍射的WVD通过阶次p的选择也可用于近场衍射。

2.1.2 定义二
由Olcay等［10］提出，基本思想是将传统的Hermite时间和频率算子推广到分数阶变量u的Hermite算子，并找到该算子的酉算子，然后利用特征函数法先找到分数阶模糊函数，再对其做二维傅里叶变换得到FrWVD。

最终表达式为
式中：α的含义同式（1）；x（t）为分析信号；Wα（t，u）为最终得到的FrWVD。

定义二与定义一对信号的处理效果是一样的，该定义的出发点是对包括定义一在内的联合分数阶时频表示（FrWVD、FrGT、STFrFT）在数学上利用Hermite算子
做出严格的、正规的推导和解释，但式（2）的复杂度（O（N3logN）较式（1）高，因此迄今未得到广泛应用。

2.1.3 定义三
由陈喆［11］提出，表达式为
式中Kα（t，u）为FrFT的核函数。

该算法的本质是由FrFT的核函数取代FFT的核函数ei2πft。

文献［11］主要利用分数阶模糊函数对三次相位信息的敏感性质，将该算法用于检测线性调频信号三次相位，但未见进一步的研究进展报道。

定义一描述的FrWVD很好地体现了FrFT的旋转特性，在阶次由0～4的增长过
程中时频图表现为相对0阶时频图（经典WVD的时频图）的顺时针旋转。

在某
些阶次下时频图较常规时频图具有更小的时频带宽积，且在所有阶次下都能以较高的时频聚集性和分辨率反映出信号的能量分布特征。

定义三描述的FrWVD几乎没有FrFT的旋转特征，在阶次由0～4的增长过程中，时频图表现出一种对能量团的拉伸和压缩过程，且当阶次取2的倍数或邻近阶次时，信号的时频图聚集在原点附近，甚至完全看不出信号的真实能量分布特征。

综上所述，定义一所描述的FrWVD更合理、更适合于非平稳地震信号处理，故文中将采用定义一的FrWVD相关知识进行理论信号仿真和实际资料处理。

2.2 最优分数阶伪WVD
本文定义分数阶伪WVD的表达式为
式中h（t）表示一维窗函数。

该式也可由分数阶STFT表示，即
式中：xp（u）的含义与式（1）相同；STFTxp（u，v）表示xp（u）的分数阶STFT，即
因此可以通过计算最优的STFT计算最优的分数阶伪 WVD。

文献［12，13］分别讨论了最优STFT的求取方法，且所求出的最优窗都能使旋转了的时频图反旋转回具有实际意义的时频图，其中文献［12］给出的分数域STFT表达式与式（6）相同。

文中利用广义时频带宽积（Generalized Time－Bandwidth Product，GTBP）准则以及FrFT的旋转特性，最终得出了基于GTBP准则的最优窗函数hGTBP（t），即
式中：γ＝Bxp／Txp，Bxp和Txp分别表示信号xp的带宽和时宽；hGTBP（t）表示最终得到的基于GTBP准则的最优窗函数。

由式（7）可见，将hGTBP（t）用于计算必然是过于复杂的。

文献［13］给出的分数域STFT表达式为
式中p、α以及Kα（t，u）的含义与前文相同。

该式类似于FrWVD的第三种定义，也是将FFT的核函数替换为了FrFT的核函数。

文献［13］巧妙地利用FrFT 的旋转特性、时移和频移特性求出了基于GTBP准则的最优窗函数，即
式中：表示信号xp的在最小时频带宽积（Time－Bandwidth Product，TBP）准则下的高斯窗函数；ropt表示作用于窗函数hTBP（xp）的最优阶次，通过计算
获取该最优阶次。

式（9）显然较式（7）精简很多，且从文献［13］的仿真结果
可知，利用式（9）的窗函数STFT可获得具有极高分辨率的时频图。

但式（9）的获取过程涉及到p和r两个不同的最优阶，文献［13］并没有很清晰地给出p和r 之间的关系。

文中将综合文献［12，13］求取最优窗的过程，推导出STFT的基
于GTBP准则的最优窗函数，即
式中：Rα 为旋转算子，定义为Rα｛D（t，f）｝＝D（tcosα＋fsinα，－tsinα＋fcosα），其中α为信号时频图顺时针旋转角度；xp为信号x的p阶FrFT；其余符号的含义与上文相同。

从式（10）可以看出，通过巧妙应用FrFT时频旋转算子，原本同时在信号和窗函数都起作用的FrFT最终完全转换到了窗函数中，且该窗函数解决了时频轴在旋转
中失去物理意义的问题。

最终的hGTBP（t）表达式为
则最优窗函数可以描述为：最优阶下信号在最小TBP准则下窗函数的－p阶FrFT。

可将以式（11）为最优窗函数的STFT应用于式（5），以得到有实际意义的最优分数阶伪WVD，即
3 理论信号仿真与分析
为了验证最优分数阶伪 WVD（式（12））去除交叉项的有效性，同时采用Choi
－William窗函数、平滑伪WVD两种经典方法对信号做处理。

仿真所用到的两组多分量信号（图1）为
式中：x1（t）为余弦信号与抛物线信号的叠加；x2（t）为余弦信号与线性调频信
号的叠加。

利用Choi－William窗函数法和平滑伪WVD去除交叉项是两种应用较广泛的方法：前者是固定核函数中满足最多WVD数学特性、且对交叉项去除具有一定作用；后者在时间和频率方向分别对原图像做平滑处理，且可以独立选择时域和频域窗函数，应用也较广泛。

图2为对x1（t）、x2（t）未去交叉项时频分布采用Choi－William窗函数、平滑伪 WVD、最优分数阶伪WVD去除交叉项效果对比，由图中可见：利用Choi－William窗函数法（图2a右上、图2b右上）和平滑伪WVD（图2a左下、图2b左下）虽然都在一定程度上抑制了交叉项，但是在信号分量交叉部分及边缘部分仍然存在明显的干扰成分；利用最优分数阶伪WVD（式（12））则可在保持较高分辨率的同时，更彻底地去除交叉干扰项（图2a右下、图2b右下）。

该例证明了文中提出的最优分数阶伪WVD在去除交叉项方面的有效性。

图1 两组多分量信号（a）x1（t）＝cosπt＋ei0.045πt3；（b）x2（t）＝cosπt ＋ei0.5πt2－0.02πt2两种信号的采样周期均为0.0625s，采样频率为16Hz，t∈［－8s，8s］
图2 对x1（t）（a）、x2（t）（b）未去交叉项时频分布（左上）采用Choi－William窗函数（右上）、平滑伪WVD（左下）、最优分数阶伪WVD（右下）
去除交叉项效果对比
4 实际地震资料处理及分析
图3为川东北地区A测线CDP叠加剖面，测井解释及钻井资料显示碳酸盐岩鲕滩储层位于1600～2100ms附近（图中椭圆标记范围），测井资料及测试结果显示该区域储层是一含气层。

笔者利用不同方法对包含目的层数据范围的地震信号进行了处理、分析和对比。

图4为基于传统Choi－William窗函数的WVD进行谱分
解得到的单频剖面，其频率范围为15～65Hz。

图5为基于平滑伪WVD的谱分解
效果。

图3 川东北地区A测线叠加剖面
图4 基于传统Choi－William窗函数的WVD进行谱分解得到的单频剖面（a）
15Hz；（b）25Hz；（c）35Hz；（d）45Hz；（e）55Hz；（f）65Hz
图5 基于平滑伪WVD的谱分解效果（a）15Hz；（b）25Hz；（c）35Hz；（d）45Hz；（e）55Hz；（f）65Hz
由图4、图5可见：当频率为15～35Hz时（图4a～图4c、图5a～图5c），产
气层段（1900～2000ms附近）的能量团逐渐增强；当频率为45～65Hz时（图
4d～图4f、图5d～图5f），强能量团逐渐消失，表明高频能量被气层吸收，含
气储层的主频向低频移动，出现了低频异常现象，依此可以识别储层。

然而图4、图5的单频剖面中虚假频率成分太多，影响高精度油气层判别精度。

为此，文中
采用本文提出的最优分数阶伪 WVD对图3数据进行处理，图6为基于最优分数
阶伪WVD的谱分解效果。

对比最优分数阶伪WVD（图6）与常规WVD（图4、图5）的谱分解效果可见，前者较后者更能清晰地反映出油气层的位置，充分说明了最优分数阶伪WVD用于地震信号处理的有效性。

图6 基于最优分数阶伪WVD的谱分解效果（a）15Hz；（b）25Hz；（c）
35Hz；（d）45Hz；（e）55Hz；（f）65Hz
5 结束语
文中总结了现有三种定义形式的分数阶Wigner－Ville分布（FrWVD）各自的优
势与局限，并将定义一的FrWVD引入地震信号频谱分解。

此外，结合最优STFT
与FrWVD之间的关系，推导了最优分数阶伪WVD的计算公式。

理论信号仿真及实际地震资料处理结果表明，最优分数阶伪WVD在一定程度上解决了限制FrWVD在信号处理中广泛应用的交叉项问题，且时频分析结果较传统WVD具有更高的时频聚集性和分辨率。

文中方法为经典时频分析理论的拓展，且能获得较好
的时频分析效果，为储层流体识别提供了一种新思路。

参考文献
［1］ Hardy H H，Richard A B，Jonathan D G.Frequency estimates of seismic traces.Geophysics，2003，68（1）：370－378.
［2］陈红，彭真明，王峻等.地震信号分数阶Gabor变换谱分解方法及应用.地球物理学报，2011，54（3）：867－873.
Chen Hong，Peng Zhenming，Wang Jun et al.Spectral decomposition of seismic signal based on fractional Gabor transform and its application.Chinese J Geophys，2011，54（3）：867－873.
［3］ Chen Y，Peng Z，He Z et al.The optimal fractional Gabor transform based on the adaptive window function and its application.Applied Geophysics，2013，10（3）：305－313.
［4］ Sejdic E，Djurovic I，Jiang J.Time－frequency feature representation using energy concentration：An overview of recent advances.Digital Signal Processing，2009，19（1）：153－183.
［5］ Zou H，Lu X，Dai Q et al.Nonexistence of cross－term free time－frequency distribution with concentration of Wigner－Ville distribution.Science in China Series F：Information Sciences，2002，45（3）：174－180.
［6］陈颖频，彭真明.基于分数域广义平滑伪Wigner－Ville分布的地震信号时频分析.中国地球物理，2012，504－505.
Chen Yingpin，Peng Zhenming.Seismic signal timefrequency analysis based on generalized smoothing pseudo Wigner－Vile distribution in fractional domain.GCS，2012，504－505.
［7］ Alieva T，Bastiaans M J.Wigner distribution and fractional Fourier transform.IEEE ISSPA，2001，168－169.
［8］ Wang Z W，Wang X L，Xiang M et al.Reservoir information extraction using a fractional Fourier transform and a smooth pseudo Wigner－Ville distribution.Applied Geophysics，2012，9（4）：391－398. ［9］ Dragonman D.Fractional Wigner distribution function.J Opt Soc Amer，1996，13（3）：474－478.
［10］ Olcay A et al.Joint fractional signal representations.Journal of the Franklin Institute，2000，337（4）：365－378.
［11］陈喆.自适应与分数阶非平稳信号处理的研究［博士学位论文］.辽宁大连：大连理工大学，2003，69－83.
［12］ Lutfiye D，Orhan A.Short－time Fourier transform：Two fundamental properties and an optimal implementation.IEEE Transactions on Signal Processing，2003，51（5）：1231－1242.
［13］ Chen Y，Peng Z.A novel method of the optimal STFr－FT and its application in seismic signal processing.International Conference on Computational Problem－Solving，2012，328－331.。