0 引言

谱分解技术是地震属性提取的重要组成部分.它以时频分析方法为基础，将地震数据分解成时间-频率-空间数据体，对该数据体进行计算分析，得到不同的属性，在沉积相划分(Feng et al., 2017)、储层预测(Sinha et al., 2005)、地质体识别(Aminu and Ojo, 2018)等方面得到了广泛应用.

目前，在储层检测中常用的时频分析方法有短时傅里叶变换(Lu and Li, 2013)、S变换(Pinnegar and Mansinha, 2003; Wu and Castagna, 2017)、小波变换(Continuous Wavelet Transform, 简称CWT) (Sinha et al., 2005, 2009; Yomogida, 1994; Daubechies, 1990)以及魏格纳-威尔分布(Wigner-Ville Distribution，简称WVD分布)(Li and Zheng, 2007)等.最初，Castagna等(2003)采用匹配追踪算法对地震信号进行分析，研究了低频阴影与油气储层的联系，而后陈学华等(2009)采用广义S变换实现了三维地震资料的低频阴影检测，减少了流体异常检测的多解性问题.但是，低频阴影仅能够定性确定油气藏的存在，并不能准确指示油气藏的具体方位.Dilay和Eastwood(1995)对地震资料进行分析时发现，在天然气生产期间，产气层的高频衰减非常明显，而其他区域则没有此类高频衰减.基于此，Mitchell等(1996)提出了计算高频衰减分析方法，通过求取地震剖面的高频衰减系数来判断油气特征.张世鑫等(2011)结合Teager-Kaiser非线性能量算子(简称TK能量算子)对高频能量衰减的计算进行了改进，得到了非线性能量衰减方法，能够更加准确指示油气藏的位置和分布范围.但是TK能量算子的适用对象是窄带信号，在应用时需要首先利用小波变换进行多尺度细化分析，增加了计算量和成本.此外，低频信息也是识别油气的一个重要方面.Silin等(2004)在频谱分解的基础上，结合Biot理论与波动方程，推导出了流体流动性属性，被广泛应用于流体识别.Chen等(2012)利用广义S变换进行频谱分解，得到分辨率较高的时频谱，并在Silin的基础上，对流动性属性公式做了补充，提出了低频频段的确定方法.

为了更快速有效地进行油气藏探测，本文提出了以同步挤压小波变换为基础的烃类检测方法，首先探测低频阴影现象，确定油气藏的大致范围，而后，针对重点工区，计算高频能量衰减与流体流动性属性，充分利用高低频信息，自适应选取优势频段，实现油气藏的精确定位.

1 基于同步挤压小波变换的烃类检测方法 1.1 同步挤压小波变换

同步挤压小波变换(Synchrosqueezing Wavelet Transform，简称SWT)是在小波变换的基础上发展起来的新型数学变换工具.小波变换可以看作是信号与经过伸缩平移后的母小波进行褶积的结果，构成了时频谱，完成对信号的稀疏化表达.

对于时间信号f(t)∈L²(R)而言，其小波变换为

(1)

其中，WT_f(a, b)是小波系数，而是母小波的复共轭，a＞0，称为伸缩因子(又被称为尺度因子)，b∈R，称作平移因子.小波基函数为：

(2)

从(2)式可知，小波变换可看作是经过伸缩平移变换后的母小波与信号的褶积，并且把褶积系数排列在时间尺度平面.尺度因子与频率有关，若把尺度映射成频率，那么计算得到的小波系数就分布于时间-频率平面.

Daubechies等(2009)在此基础上，对小波变换进行改进，提出了同步挤压小波变换.

根据帕塞瓦尔定理，公式(1)可以改写为

(3)

其中，ε表示角频率，的傅里叶变换，尺度因子a对复小波的频率进行挤压，e^jbe是平移因子b的傅里叶变换.所以，公式(1)中时间域的褶积通过帕塞瓦尔定理就可以改写成公式(3)中频率域的乘积形式.沿着频率轴进行压缩处理.对于所有的W(a, b)≠0，它的瞬时频率ω_s(a, b)是连续小波变换W(a, b)关于b的偏导，即：

(4)

接下来，对其进行挤压操作.利用尺度与频率之间的关系，把计算得到的小波系数W(a, b)从时间尺度平面映射到时间-频率平面(b, a)→(b, ω(a, b))，同步挤压小波变换的系数只在以ω_l为中心的频率区间[ω_l-Δω/2, ω_l+Δω/2]内，同时Δω=ω_l-ω_l-1，同步挤压小波变换为

(5)

(5) 式表明，同步挤压小波变换是将小波变换的时频谱在频率方向上进行挤压处理得到的.同步挤压小波变换的逆变换为

(6)

其中，C_ψ是与所选小波有关的常数.图 1是复合信号s(t)的示意图，由三个分信号组成，分别为

其中，t∈[0, 10]，共有2048个采样点.分别采用小波变换和同步挤压小波变换进行时频分析得到图 2.从图 2中可以看出，小波变换时频谱的能量聚焦度远低于同步挤压小波变换，并且，SWT的时频谱上三个信号的频率成分非常清晰明了，而CWT时频谱上，频率成分并不清晰.因此，SWT得到的时频谱的分辨率更高，有利于后续对信号的深入分析.

1.2 流体流动性估计

Goloshubin等(2006)提出将流体活动性属性表示为地层含流体时地震振幅与频率的比值.当储层中含气时，振幅与频率的比值较不含气时在低频端表现为正异常，而在高频端表现为负异常.基于此，Silin等(2004)将Biot理论与波动方程结合，推导出既可描述地震波的动力学特征，又能将储层密度、渗透率和流体黏度等油藏参数联系起来的弹性波方程，被广泛应用于流体识别中.

考虑到动力学原理和孔隙弹性介质，地震波的低频近似表达可认为是弹性波在含流体的孔隙介质中传播，反射系数可表示为

(7)

其中，R₀和R₁是反射系数，无量纲，，k是渗透率，η是储层流体黏度，ρ_b是储层的体积密度.流体流动性因子M则定义为岩石渗透率与孔隙中流体黏度的比值，即：

(8)

将公式(8)代入公式(7)，并对角频率求偏导，即：

(9)

由(9)式可得：

(10)

定义，代入(10)式得：

(11)

该式表明，流体流动性因子M与的平方呈正比关系.Chen等(2012)提出F可以通过测井数据计算出来，可看作是一个比例因子，因此，可忽略该比例因子，通过计算来计算相对关系.在频谱分解后，瞬时振幅谱可以代表在该频率处的反射系数，因此

采用SWT计算地震记录的时频谱，则在每个时间点上，均可得到一个瞬时振幅谱A(t, ω)，则为

(12)

为了更好地表达相对关系，对其进行归一化处理，即：

(13)

基于同步挤压小波变换分析高精度流体流动性属性的常规流程为：首先对数据进行SWT时频分析，并计算得到每一道地震数据的时频谱.其次，计算比较不同频段的数据特征，选取合适的低频频段.随后对选定的频率段内的瞬时振幅谱进行斜率计算，最后进行归一化处理，寻找含油气区域.

在此过程中，关键步骤之一是如何确定合适的低频区间.常规算法根据地震数据的不同确定不同的频段，使得选取过程计算量大，且耗时较长.通过实际数据测试，本文首先计算数据的重心频率和重心频率所对应的频谱能量，再选取占重心能量所对应频谱能量的10%~90%所对应的频率段作为计算频段，即频段位置由其积分能量决定.

1.3 高频能量衰减梯度

频率吸收衰减在频谱分析中占据重要地位.地震波的吸收衰减反映了地下介质的非弹性属性.当地层中含有油气时，地震波在传播过程中会产生较为强烈的吸收衰减.这种吸收衰减主要发生在信号的高频能量部分，因此，Mitchell等(1996)提出了能量衰减梯度这一属性来衡量地层的吸收衰减能力.这一属性对烃类反应较为敏感，在油气藏检测中得到广泛应用和推广.

能量衰减梯度技术的核心是在地震信号的时频谱中检测高频能量的衰减增加，因此该项技术又被称为高频能量衰减梯度.黏滞弹性理论指出，在地震波的传播过程中，由于地层的吸收衰减作用，地震波的能量衰减与传播距离之间的关系可由(14)式描述：

(14)

式中，A₀是地震波的初始振幅，α是吸收系数，f是频率，A是地震波传播一定距离后的振幅.在计算过程中，利用时频变换将地震数据从时间域转换到时间频率域，在时频域计算每个时间点上的衰减系数，建立时间与衰减系数的关系.地震波的固有衰减与地层结果、裂缝及流体等有关，若地层中含有油气，则油气对高频部分的吸收将明显高于其他因素带来的吸收衰减作用，显示为衰减异常.

常规衰减梯度的计算流程为，对每一道地震记录进行时频分析处理，在每个时间点上抽取瞬时频谱，计算占地震总能量的65%及85%的两点频率，将过这两点的斜线斜率作为衰减梯度a，即：

(15)

式中y_85%和f_85%及y_65%和f_65%分别代表占总能量85%和65%时所对应的能量和频率.

在实际数据计算时，每个时间点的瞬时频谱的高频部分能量衰减形态并不是线性的，而呈现指数型的衰减趋势，因此，本文采用非线性拟合来对高频部分能量趋势进行拟合，则有：

(16)

其中，b表示高频吸收衰减梯度，B为初始振幅.本文采用Nelder-Mead算法对公式(16)进行非线性求解，得到高频吸收衰减梯度.

为了更稳定地计算能量衰减梯度，利用重心频率来选取频段.信号的重心可以反映出信号总体频率及能量成分的变化，利用多级重心的概念对频率进行分段分析可更加稳定地讨论信号的能量变化.以信号的重心为分界线，在重心两侧分别按照重心的定义计算各自重心，称之为二级重心，以此类推，可计算信号的三级重心及多级重心.

假设信号x(t)，其频谱为U(f)，如图 3所示，信号的一级重心如图中圆点所示，其计算公式为

(17)

信号的二级重心如图中方点所示，计算公式为

(18)

图中三角形为信号的三级重心，计算公式为

(19)

本文仅用到三级重心，选取一级重心与三级重心〈f_c〉_{3, 4}之间的频段为待计算频段，即图 3中两条黑线之间.一般而言，三级重心频率〈f_c〉_{3, 1}和〈f_c〉_{3, 4}之间的频段为信号频谱的主频带，二级重心频率〈f_c〉_{2, 2}和三级重心频率〈f_c〉_{3, 4}为其高频部分.由于实际信号通常较为复杂，为确保结果稳定，选定一级重心频率与三级重心频率〈f_c〉_{3, 4}之间的频段为待处理的高频频段.选定频段后，对该频段内的能量变化进行非线性拟合，得到能量衰减梯度.

1.4 工作流程

在实际处理中，首先对地震剖面进行低频阴影检测，定性确定油气藏的位置，再对数据进行流体流动性属性和高频能量衰减梯度计算，精确描述油气藏位置，整体流程如图 4所示，即：

(1) 对地震数据应用同步挤压小波变换计算其时频谱；

(2) 以(1)中的结果为基础，提取单频剖面；

(3) 寻找低频阴影，定性确定油气藏范围；

(4) 对工区数据进行流体流动性属性和高频能量衰减梯度计算；

(5) 分析(4)中的结果，确定烃类的分布范围和位置.

2 模型算例

为了验证本文提出的计算流体流动性属性的可行性，首先采用Marmousi Ⅱ模型部分数据进行试验.图 5为Marmousi Ⅱ部分模型数据，在CDP420~550、时间深度1.35~1.45 s的范围里存在一砂岩透镜体，此透镜体内含有油气，顶底反射振幅较大.

单道计算示意图如图 6所示，(a)为地震道，在1.3~1.4 s范围内含有油气，(b)为相应的时频谱，(c)为某一时间样点的瞬时谱，通过计算，重心频率为40 Hz，计算其积分能量，确定计算频段为19~38 Hz，(d)为该地震道所对应的流体流动性曲线，(e)为采用小波变换得到的流体流动性曲线，在含油气区域内，流体流动性明显增强，与模型情况相符.对该地震剖面进行流体流动性计算，得到最终结果如图 7所示，(a)和(b)分别是采用SWT和CWT所得到的结果，图中透镜体的流动性明显高于其他区域，与模型一致.对比图 6d和6e可知，在含油气部分，小波变换得到的流体流动性的范围较大，而同步挤压小波变换得到的区域范围较为精确.对比图 7a和7b可知，小波变换得到的流体流动性结果的分辨率低于同步挤压小波变换得到的结果，这一现象与两种变换的分辨率有关.因此，基于同步挤压小波变换的流体流动性结果能够较为精确地指示模型中油气藏的位置.

图 8为Marmousi Ⅱ模型的一个含气地震道以及该道对应的时频谱，与图 5中采用相同的地震道，图 8c为某一时刻对应的瞬时频谱，占能量65%和85%所对应的频率分别为45 Hz和60 Hz，连接这两处频率能量的直线为应用公式(15)计算得到的斜率，而一级重心频率与三级重心频率〈f_c〉_{3, 4}分别为40 Hz和65 Hz，连接这两处频率能量的粗曲线代表非线性拟合方法计算的衰减梯度.由图可知，采用本文提出的频段选择与非线性拟合的方法能够更加稳定准确地描述能量衰减规律.图 8d为采用本文提出的方法得到的该地震道的衰减梯度，图 8e为传统算法得到的地震道的衰减梯度.两者对比可知，在含气部分，图 8d中曲线的值明显增加，图 8e曲线对应区域的值也有所增加，但是其规律性比图 8d明显较差.采用本文提出的新算法对图 5的地震剖面进行计算，结果如图 9所示，图中透镜体的值明显大于其他区域，表明该区域内的衰减梯度较大，则可能含有油气，与模型情况相符.

3 实际资料应用

采用一套海上叠后数据进行处理.图 10为该工区数据的示意图，该套数据共有232条主测线，249条联络测线，井位于图中所示位置处，目前已探明在井位置处深度约4 s的地方含有油气.含气储层为第三系浊积砂岩，储层埋深约3291~3779 m，在叠加地震剖面上显示为亮点响应，横向上反射同相轴连续性较好.图 11为过井剖面，该数据的采样率为4 ms，时间范围为2~5 s，共249道，井资料显示，该工区内的含气层大约位于3.9 s附近.采用低频阴影进行定性检测，再利用流体流动性属性与高频能量衰减梯度进行精细分析，三种结果联合进行含气储层位置的确定.

首先，采用SWT对该数据进行处理，计算其低频阴影.为了进行对比，传统的小波变换谱分解也用于处理该数据.图 12(a)、(c)、(e)为CWT的分频结果，而图 12(b)、(d)、(f)为SWT的分频结果，两者均展示了20 Hz、40 Hz和60 Hz的分频结果.图中椭圆代表已知油气藏的位置.

从图中可知，CWT与SWT在频率切片的4 s附近均表现出了低频能量强而高频能量弱的特征，能量减弱的位置恰好位于含气储层的下方，符合传统“低频阴影”的认知，成功进行了低频阴影的检测，与实际情况相符.需要注意的是，在进行时频分析时，与CWT相比，SWT的时频分辨率更高，能量聚焦性更好，并且SWT的分频剖面上，能量的衰减表现得比CWT的结果更强烈.

为了进一步验证低频阴影的结果，对同一过井剖面进行流体流动性属性计算.与低频阴影检测的试验方式类似，将其与以小波变换为基础的计算结果进行比较.

图 13是实际地震道的流体流动性计算图，4 s左右振幅增强，时频谱上对应位置能量较大，选取某一时刻的瞬时谱，重心频率为32 Hz，积分能量为重心频率能量的10%和90%所对应的频率分别为10 Hz和29 Hz，计算该道的流体流动性结果如图 13d所示.在4 s附近，流体流动性明显增强，高于其他时间段，证明该时间段对应的地层含有油气.图 14是整个地震剖面的流体流动性因子计算结果，(a)与(b)分别是利用CWT与SWT的计算结果.从图中可以看出，椭圆区域内的流体流动性属性明显增大，呈现亮点，说明此处含有油气，结合之前低频阴影的识别结果，可以进一步确定此处含有油气.虽然图 14a与14b均指示出了油气的位置，但是采用SWT计算得到的结果分辨率更高，油气藏的位置指示更加明显，整个剖面比小波变换得到的剖面更加干净明了.

对于单个地震剖面，基于同步挤压小波变换的流体流动性属性计算结果符合实际勘探结果.因此，对整个工区数据进行流体流动性计算，结果如图 15a和15b所示.在时间深度为4 s左右，空间范围在如图 10中方框所示的范围内，流体流动性较大，该区域结果如图 15c和15d所示.由图可知，油气藏的深度范围在3.95~4.05 s之间，平面范围大约在160~197号联络测线和192~220号主测线所围起来的平面内，整体呈现一个扁平圆柱状.流体流动性属性的结果与实际钻井结果相符，可行性较高.分析对比图 15a与15b，通过同步挤压小波变换得到的结果更清晰，含油气部分与不含油气部分的颜色对比更加明显，边缘更为锐化，能够更为清晰地表达油气藏的范围与位置.因此，基于同步挤压小波变换的流体流动性属性计算在实际数据勘测中是可行的，其结果与实际勘探结果符合度较高.

本文对该套数据计算其频率衰减梯度.首先，计算过井剖面的能量衰减属性.单道结果如图 16所示，图中地震道与图 5中一致，在大约4 s左右含有油气，该地震道的能量衰减梯度与流体流动性的计算结果一致，均在4 s左右，衰减梯度增加，证明该处的衰减增强，可能含有流体.

图 17是过井剖面的能量衰减梯度计算结果，在CDP范围150~200，时间深度3.9 s的区域内，衰减梯度增强，与该处含有流体对应，与测井结果一致.图 18a是整个工区的计算结果，平面范围上，与图 10中方框部分对应区域在4 s左右深度范围内的能量衰减梯度明显增加，颜色较深，为了更好地了解该区域的高频能量衰减梯度变化的情况，将该区域的计算结果单独分析，如图 18b所示，在深度大约为4 s左右，平面范围大约在160~197号联络测线和192~220号主测线所围起来的平面内，衰减梯度明显增加，并且其形状与流体流动性估计得到的结果一致，与已知油气藏的位置相符，证明基于SWT的高频能量衰减梯度算法具有较好的实际应用价值，其实际资料处理结果可信.

对比图 18与图 16，两者在相同的地方均出现较高的值，能够较好地指示油气的位置.纵观整个剖面，流体流动性属性的结果更为干净，说明高频能量衰减梯度的算法存在较大的改进空间.

此外，由于流体流动性估计和高频能量衰减梯度需要用到低频段和高频段的信息，所采用的地震资料应尽量保留完整频带信息，以提高本方法的准确率.同时需要注意的是，由于薄层调谐作用，决定薄互层反射特征的主要因素不是岩性或油气，而是储层结构，为了保持振幅的相对关系，在进行研究时，地震资料应符合保幅要求.

4 结论

频谱分解是对地震资料分析的常用手段.同步挤压小波变换得到的时频谱比常规小波变换得到的时频谱分辨率更高，在后续的处理中，其结果的分辨率亦优于小波变换的结果.

本文提出了基于同步挤压小波变换的频谱分解烃类检测技术，包括基于同步挤压小波变换的低频阴影检测、流体流动性属性计算和高频能量衰减梯度等，低频阴影可定性指示出油气藏的方位，流体流动性属性与高频能量衰减梯度可较为准确地定位油气藏，综合利用地震信息的高低频信息，实现了多属性联合检测油气藏.通过对实际地震资料的处理表明：基于同步挤压小波变换对地震数据进行低频阴影检测、流体流动性估计以及高频能量衰减梯度计算结果均优于常规小波变换.本文以重心频率为基础，定义了流体流动性属性与高频能量衰减梯度的频率计算范围，与常规算法相比，大大降低了计算量.此外，采用非线性拟合对高频能量衰减梯度进行重定义，更为准确地描述了高频信息的衰减.