唐湘蓉，石战战，彭真明

(1.成都理工大学地球探测与信息技术教育部重点实验室，四川成都610059;2.成都理工大学地球物理学院，四川成都610059;3.成都理工大学工程技术学院，四川乐山614000;4.电子科技大学光电信息学院，四川成都610054)

时频分析技术已成为利用地震资料进行储层预测和烃类检测的重要方法之一。地震波在地下介质中传播时，由于地下介质的复杂性、含油气性以及对地震波吸收衰减程度的差异，地震信号呈现明显的时变和非平稳特征。对于非平稳信号常用的时频分析方法有：短时傅里叶变换(STFT)，Gabor变换(GT)，小波变换(WT)，S变换(ST)，以Wigner-Vile分布(WVD)为代表的Cohen类分布等[1-3]。STFT受窗函数时频分辨率的制约，且时频分辨率是固定的;WT具有多分辨率特性，但本质上是时间-尺度谱;ST变换结合了STFT和WT的特点，具有较高的时频分辨率，但是其频谱图较粗糙，而且受测不准原理的制约;Wigner-Vile分布具有较高的时频分辨率，但由于受交叉项影响严重，其应用也受到了限制。分数域时频分析方法(分数阶傅里叶变换(FrFT))最早由Victor在1980年提出[4]，由于其较高的能量聚集性、可在不同阶次下进行分频计算等优点，已应用于光学信号处理、参数估计、盲信号处理、雷达目标检测等领域[5-8]。该方法在地震信号处理中也逐渐引起重视，Montana等[9]首次利用FrFT进行地震信号的空间预测滤波;刘喜武等[10-11]和刘婉莹等[12]研究了地震信号分数阶域短时Fourier变换、加权伪Wigner分布以及分数阶伪Wigner分布的数值实现方法;陈红等[13]研究了地震信号分数域Gabor变换;彭建亮等[14]提出了基于分数域自适应滤波的地震信号去噪方法;陈颖频等[15-16]基于分数域广义平滑伪Wigner-Vile分布对地震信号进行了时频分析;Tian等[17]针对复杂地震信号提出了基于自适应窗函数的最优分数域Gabor变换;Xu等[18]将S变换推广到分数阶S变换并应用于地震信号处理;Wang等[19]利用FrFT的旋转性进行储层信息提取。Teager能量最早是由Teager等提出的一种能量操作算子[20]，后经Kaiser推导出其离散形式[21]，即Teager-Kaiser(T-K)非线性能量算子，主要用于研究非线性过程中的语音识别、噪声压制、图像增强及边缘检测[22-24],近年来，在地震信号处理中也逐渐引起重视[25-27]。

本文研究分数阶Gabor变换(FrGT)的理论及FrGT最优阶的确定方法，基于最优阶FrGT对地震资料进行谱分解，计算并提取单频分量的T-K能量的瞬时属性。理论信号仿真和实际资料的处理结果表明，FrGT可有效用于地震信号的时频分析，T-K非线性能量可作为储层预测的一种有效属性。

1 方法原理

1.1 FrFT原理

时域信号x(t)的q阶FrFT是一个线性积分运算[13]：

(1)

式中：q为阶数;Xq(u)为x(t)的q阶FrFT;Kq(t,u)为核函数，可以进一步表示为

(2)

式中：j为虚数单位;δ(t)为冲击函数;α=πq/2,q≠2n为时频平面旋转角度;n为正整数。

由(2)式可以看出，当α=0时，x(t)的FrFT为其本身;当α=π/2时，x(t)的FrFT退化为传统傅里叶变换。若将傅里叶变换看作从时间轴t逆时针旋转π/2到频率轴f，则FrFT可看作旋转任意角度α到u轴的算子[6]，FrFT与传统FFT的时频分布关系如图1所示。

图1 FrFT与传统FFT的关系

1.2 FrGT原理

FrGT是对传统Gabor变换应用分数阶化思想的推广[1]，在Gabor核函数中引入了FrFT核函数因子，可以看作是对FrFT的加窗处理，信号x(t)的q阶Gabor展开式如下[16]：

(3)

其中，gm,k,q为展开系数;h(t)为综合窗函数;hm,k,q(t)=h(t-mT)Wq,k(t)，是由分数阶综合窗函数进行平移和调制生成的，Wq,k=exp{j[-1/2(t2+(kΩsinα)2)cotα+kΩt]}，T和Ω分别为时间域和频率域采样间隔，m和k分别为时间域和频率域采样点数。

由(3)式可知，当TΩ>2π时，时频域采样栅格过稀，缺少足够信息恢复原信号x(t)，称为欠采样;而当TΩ<2π时，时频域采样栅格过密，类似一维信号采样率过高，会产生冗余，称为过采样;当TΩ=2π时，称为临界采样。

由(1)式得到信号x(t)的q阶FrGT为

(4)

其中，γm,k,q(t)=γ(t-mT)Wq,k(t)，γ(t)为正交分析窗，“*”为复共轭算子。将(4)式代入(3)式得到FrGT完备性条件：

(5)

当给定综合窗h(t)时，由(5)式可以求出分析窗函数γ(t)，代入(4)式可以求出信号的q阶FrGT。实际应用时总是取TΩ<2π，即过采样，这样能够保证FrGT的稳定。

1.3 FrGT最优阶确定

非平稳信号时频分析中，时间分辨率和频率分辨率是一对矛盾体，既具有任意高的时间分辨率又具有任意高的频率分辨率的时频分析方法是不存在的，任何时频表示方法都只能在一定程度上近似信号在时频点(t,f)的能量密度。要获得分数域时频聚集性最好的一个关键问题就是确定分数域变换的最优阶次。根据不确定原理，信号分析时频分辨率由时宽-带宽积(time-bandwidth product,TBP)决定，TBP越小，它的时频分辨率越高。

连续信号x(t)的时宽-带宽积定义为[17]

(6)

式中：Tx为信号x(t)的时间分辨率，即有限时间宽度(时宽);Bx为信号x(t)的频率分辨率，即有限频率宽度(带宽)。

对(6)式进行拓展得到分数阶的广义时宽-带宽积的表达形式：

(7)

式中：xq(t)为信号x(t)的q阶FrGT;Tq和Bq分别为分数阶时宽和带宽，表达式如下：

1.4 T-K非线性能量计算原理

Teager能量算子最早是由Teager等在研究语音信号非线性建模时提出，能够跟踪信号瞬时能量[20]。连续信号的Teager能量算子方程为

(10)

式中：x(t)为连续单频信号;ΨTx(t)为连续单频信号的Teager能量算子。

Kaiser推导了方程(10)的离散形式，也称为Teager-Kaiser(T-K)能量方程：

ΨTx(n)=x2(n)-x(n-1)x(n+1)

(11)

式中：x(n)为单频离散时间序列;ΨTx(n)为离散序列的T-K能量离散算子。

大量研究表明[21-27]，单频离散信号某时刻的T-K能量与振幅平方、频率的平方成正比。T-K能量是刻画波场能量密度、表征波场能量分布的一种有效算子，对于单频信号是严格成立的。但实际地震信号是复杂的，会产生交叉因子的干扰，因此需利用优良的时频分析方法将复杂的地震信号分解成一系列单频信号的线性组合，对每一单频信号分别计算T-K能量。

1.5 计算流程

FrGT算法中的关键是确定最优阶，阶数的选取直接影响算法对信号的分析性能。因为FrGT具有周期性，为了得到最优阶，算法实现过程中可在一个周期内对阶数进行遍历，实现流程为：

1) 输入信号x(n)，设定窗函数h(n)，计算和确定正交分析函数γ(n);

2) 根据正交分析函数γ(n)对输入信号x(n)进行加窗处理;

3) 确定分数域变换小的最优阶次,设q,Δq,qmin,qmax,qopt分别表示FrGT的任意分数阶、递增步长、最小阶、最大阶数及最优阶，令q=qmin,∀q∈[qmin,qmax]，给定步长Δq，令q=q+Δq，计算分数域下广义时宽-带宽积GTBP，优选出使GTBP最小的最优阶qopt;

4) 根据公式(4)计算最优阶次下的FrGT变换，得到信号的分数域时频分布，分数域频率可通过反旋转换算为实际信号x(n)的真实频率;

5) 对谱分解后每一单频数据计算T-K瞬时能量，构造T-K瞬时能量谱数据体;

6) 对T-K瞬时能量谱数据体求取瞬时属性用于储层预测。

2 理论信号仿真和实际信号测试

为了验证上述方法的有效性，以线性调频信号(linear frequency modulation signal,LFM)作为测试信号。LFM是公认检验时频分析方法优劣的模型之一。图2a中合成信号由两个起始归一化频率为0.01和0.10的LFM叠加而成;计算中，时间t取0～250ms;图2b为图2a合成信号的STFT，可以看出信号各频率成分，但信号时频分辨率较低;图2c为合成信号的WVD分布，虽然时频能量局部化有较大的改善，具有很高的时频分辨率，但交叉项影响严重，其时频谱解释困难，影响其在非平稳信号时频分析中的应用;图2d到图2f为合成信号不同阶次的FrGT，当q=1.850(图2e)时，其时频分辨率可以和WVD分布相比较，且不受交叉项影响，检测信号局部特征能力明显提高，而q=1.570(图2d)和q=2.198(图2f)时，分辨率较低，FrGT优越性不能体现。不同阶数下对应的GTBP如表1所示，当q=1.850时，GTBP最小，其FrGT时频分布具有更高分辨率。

表1 合成信号不同阶次q下的GTBP

图3是新疆TH油田某区块地震信号的T-K能量谱瞬时分析图。其中，图3a为过T1井单道地震记录，含气储层位置如图中虚框所示;图3b为最优阶FrGT时频谱;图3c是基于最优阶(q=1.06)FrGT的T-K能量谱;对比图3b和图3c可见，在已知的含油气层段(经测井与开发证实的)，均明显可见高能量分布特征(图中红色箭头所示位置)，相对于最优阶(q=1.06)FrGT时频谱，T-K能量谱时间延续度较小，具有更高的时频聚焦性和时间分辨率，其时间局部化能力更强，更有利于精确地指示储层位置。

图2 合成信号的FrGT仿真结果a 合成LFM信号; b STFT时频谱; c WVD时频谱; d FrGT时频谱(q=1.570); e FrGT时频谱(q=1.850); f FrGT时频谱(q=2.198)

图3 TH油田某区块单道地震信号的T-K能量谱分析a 过井单道地震记录; b 最优阶(q=1.06)FrGT谱; c 基于最优阶(q=1.06)FrGT的T-K能量谱

3 实际资料处理结果分析

为了进一步评价基于分数阶时频域T-K能量储层预测方法的有效性，选取了新疆TH油田某区块经过保幅处理的三维叠后偏移地震数据进行处理分析。该区三叠系主力储层是中油组，以辫状河河道砂体沉积为主。图4是实际地震资料的剖面分析图。图4a为过井剖面，目的层在2970～3050ms，含气储层位于粉色椭圆虚框内;图4b是均方根振幅剖面，储层发育带显示高振幅异常，但相对背景差异较小，难以精确圈定储层范围;图4c是q=1.00时FrGT的T-K非线性能量剖面，含气储层也显示高能量的“亮点”反射，岩性边界较清晰，但储层内部非均匀性模糊;图4d是最优阶q=1.06时FrGT的T-K非线性能量剖面，储层高能量异常更为突出，对储层范围和岩性边界的检测更精确，且储层内部非均匀性刻画明显变好。图5是实际地震资料的平面分析图。图5a，图5b，图5c和图5d分别是沿目的层抽取的原始振幅、均方根振幅、q=1时的T-K非线性能量和q=1.06时的T-K非线性能量切片。从图5可以看出,在q=1.06时，瞬时T-K非线性能量切片中不仅油气储层的高能量异常明显,而且储层分布范围、岩性边界及储层内部非均匀性刻画得也是最清晰，该范围内已有多口高产油气井证实。

图4 实际地震资料处理剖面分析a 过井原始地震振幅剖面; b 均方根振幅剖面; c FrGT(q=1.00)T-K能量剖面; d FrGT(q=1.06)T-K能量剖面

图5 实际地震资料处理切片分析a 原始振幅切片; b 均方根振幅切片; c FrGT(q=1.00)T-K能量切片; d 最优阶FrGT(q=1.06)T-K能量切片

4 结论与讨论

FrGT是针对非平稳信号提出的一种高精度时频分析方法，通过优选出最优阶次能够获得最佳的时频分辨率。T-K能量是一种非线性能量，能反映地震信号能量密度随时间变化的主要信息，有效地突出强能量，使“亮点”更亮，能有效反映含油气储层的空间展布及内部的非均匀性。本文经最优阶的FrGT谱分解后每一单频分量计算T-K能量，构造出的T-K能量谱具有更高的时频聚焦性和时间分辨率，能有效突出信号中高能量部分。理论信号仿真和实际资料试处理结果表明，基于分数阶时频域Teager-Kaiser非线性能量的储层预测方法，能够较精确地检测储层空间展布、岩性边界及内部非均匀性，提高储层预测的精度和可信度。

参 考 文 献

[1] Portnoff M R.Time-frequency representation of digital signals and systems based on short-time Fourier analysis[J].IEEE Transactions on Acoustics,Speech,and Signal Processing,1980,28(1):55-69

[2] Zibulski M,Zeevi Y Y.Oversampling in the Gabor scheme[J].IEEE Transactions on Signal Processing,1993,40(3):2679-2687

[3] Simon C,Ventosa S,Schimmel M.The S transform and its inverses:side effects of discrediting and filtering[J].IEEE Transaction on Signal Processing,2007,55(10):4928-4936

[4] Victor N.The fractional order Fourier transform and its application to quantum mechanics[J].IMA Journal of Applied Mathematics:Institute of Mathematics and its Applications,1980,25(3):241-265

[5] Ozaktas H M,Ankan O,Kutay M A,et al.Digital computation of the fractional Fourier transform[J].IEEE Transactions on Signal Processing,1996,44(9):2141-2150

[6] Ozaktas H M,Kutay M A.Introduction to the fractional Fourier transform and its applications[J].Advances in Imaging and Electron Physics,1999,106:239-291

[7] Aydin A,Yalcin C.A fractional Gabor expansion[J].Journal of the Franklin Institute,2003,340(5):391-397

[8] Ervin S,Igor D,Jubisa S.Fractional Fourier transform as a signal processing tool:an overview of recent developments[J].Signal Processing,2011,91(6):1351-1369

[9] Montana C A,Margrave G F.Spatial prediction filtering in fractional Fourier domains[J].Expanded Abstracts of 72ndAnnual Internat SEG Mtg,2004,241-244

[10] 刘喜武,刘婉莹,刘洪,等.地震信号广义时频分析及其数值实现[J].物探化探计算技术,2007,29(5):386-390

Liu X W,Liu W Y,Liu H,et al.Generalized seismic signal time frequency analysis and numerical algorithms[J].Computing Techniques for Geophysical and Geochemical Exploration,2007,29(5):386-390

[11] 刘喜武,张宁,勾永峰，等.地震勘探信号时频分析方法对比与应用分析[J].地球物理学进展,2008,23(3):743-753

Liu X W,Zhang N,Gou Y F,et al.The comparison and application of time frequency analysis methods to seismic signal[J].Progress in Geophysics,2008,23(3):743-753

[12] 刘婉莹,孟凡顺,刘喜武,等.地震数据分数阶域空间预测滤波方法研究[J].中国海洋大学学报,2008,38(6):1008-1012

Liu W Y,Meng F S,Liu X W,et al.Study on spatial prediction filtering in fractional fourier domains for seismic data[J].Periodical of Ocean University of China,2008,38(6):1008-1012

[13] 陈红,彭真明,王峻,等.地震信号分数阶Gabor变换谱分解方法及应用[J].地球物理学报,2011,54(3):867-873

Chen H,Peng Z M,Wang J,et al.Spectral decomposition of seismic signal based on fractional Gabor transform and its application[J].Chinese Journal of Geophysics,2011,54(3):867-873

[14] 彭建亮,彭真明,张杰,等.基于分数域自适应滤波的地震信号去噪方法[J].地球物理学进展,2012,27(4):1730-1737

Peng J L,Peng Z M,Zhang J,et al.De-noising method of seismic signal based on adaptive filtering in fractional domain[J].Progress in Geophysics,2012,27(4):1730-1737

[15] 陈颖频,彭真明.基于分数域广义平滑伪Wigner-Vile分布的地震信号时频分析[C]∥第二十八届中国地球物理年会论文集.北京：中国地球物理学会，2012：504-505

Chen Y P,Peng Z M.Seismic signal time-frequency analysis based on generalized smoothing pseudo Wigner-Vile distribution in fractional domain[C]∥28thChinese Geophysical Annual Conference Expanded Abstracts.Beijing:Chinese Geophysical Society,2012:504-505

[16] Chen Y P,Peng Z M,He Z H,et al.The optimal fractional Gabor transform based on the adaptive window function and its application[J].Applied Geophysics,2013,10(3):305-313

[17] Tian L,Peng Z M.Determining the optimal order of fractional Gabor transform based on kurtosis maximization and its application[J].Journal of Applied Geophysics,2014,108:152-158

[18] Xu D P,Guo K.Fractional S transform,part 1:theory[J].Applied Geophysics,2012,9(1):73-79

[19] Wang Z W,Wang X L,Xiang W,et al.Reservoir information extraction using a fractional Fourier transform and a smooth pseudo Wigner-Ville distribution[J].Applied Geophysics,2012,9(4):391-400

[20] Teager H,Teager S.Evidence for nonlinear production mechanisms in the vocal tract[J].Speech Production & Speech Modeling,1990,55:241-261

[21] Kaiser J F.On a simple algorithm to calculate the energy of a signal[J].IEEE CASSP,1990(2):381-384

[22] Jabloun F,Cetin A E,Erzin E.Teager energy based feature parameters for speech recognition in car noise[J].IEEE Signal Processing Letters,1999,6(10):259-261

[23] 高亚召,赵霞.基于Teager能量算子的自适应小波语音增强[J].电声技术,2009,33(1):58-62

Gao Y Z,Zhao X.Adaptive thresholding wavelet speech enhancement based on Teager energy operator[J].Audio Engineering,2009,33(1):58-62

[24] 李辉,郑海起,杨绍普.基于EMD和Teager能量算子的轴承故障诊断研究[J].振动与冲击,2008,27(10):16-18

Li H,Zheng H Q,Yang S P.Bearings fault diagnosis based on EMD and Teager Kaiser energy operator[J].Journal of Vibration and Shock,2008,27(10):15-17

[25] Matos D,Marfurt K J.Wavelet transform Teager-Kaiser energy applied to a carbonate field in Brazil[J].The Leading Edge,2009,28(6):708-713

[26] 熊晓军,王小品,贺振华.基于非线性能量分析的鲕滩储层预测方法[J].西南石油大学学报(自然科学版),2010,32(4):1-5

Xiong X J,Wang X P,He Z H.Oolitic reservoir prediction method based on nonlinear energy analysis[J].Journal of Southwest Petroleum University(Science&Technology Edition),2010,32(4):1-5

[27] 陈学华,钟文丽,贺振华.基于时频域Teager主能量的储层检测[J].石油地球物理勘探,2011,46(3):434-437

Chen X H,Zhong W L,He Z H.Reservoir characterization based on Teager dominant energy in time-frequency domain[J].Oil Geophysical Prospecting,2011,46(3):434-437