基于VMD的高分辨率频率衰减分析算法

薛雅娟，陈辉，刘哲哿，杨佳，钟剑丹，周娟，李东芳

1.气象信息与信号处理四川省高校重点实验室(成都信息工程大学)，四川 成都 610225

2.数学地质四川省重点实验室(成都理工大学)，四川 成都 610059

地震波的能量衰减是地震波在地下介质中传播期间发生的能量损耗，是介质的一种内在属性。除了岩层的不均匀性和几何扩散等造成的弹性损耗会导致地震波传播过程中的能量衰减以外，岩石的黏弹性性质也会导致地震波能量衰减，使地震波的弹性能转化为热能然后耗散掉，造成地震波的总能量衰减[1-6]。由地震波弹性势能转化为热能形成的地震波本征衰减与各种介质和流体特性(如流体含量，渗透率，黏度)有关[7]，使其成为潜在的储层监测工具[8-9]。在相对稳定的垂直和水平岩性几乎没有变化的地层结构中，地震波的衰减主要是由流体特性引起的[10-13]。实验室和实际数据测量均表明，在大多数频率带宽上黏性流体饱和岩石的地震波衰减比干燥岩石的地震波衰减更为明显，且地震波的高频能量衰减大于低频能量衰减[10-14]。

能量吸收分析方法[15]是一种通过计算吸收系数来估计介质高频能量吸收的有效的衰减估计方法，已被广泛用于油气探测中[13-21]。能量吸收分析方法的关键是准确计算吸收系数或者是衰减梯度。常规能量吸收分析方法使用两点斜率法计算吸收系数，仅适用于旁瓣很少或几乎没有旁瓣的平滑频谱或具有高信噪比(SNR)的地震数据[17]。近年来为了提高衰减估计的准确性，具有更高时频分辨率的时频分析方法，如小波变换等被引入能量吸收分析方法中，同时，为了提高吸收系数计算的准确性，最小二乘法也被引入替代传统的二点斜率法[15-16]，为避免不同频率地震波衰减不同造成的相互影响而采用自适应经验模态分解及其衍生方法[13,17-19]、量子信息处理算法[20]的高精度衰减估计方法也引入到能量吸收分析中。

对于地震波反射数据，能量吸收分析方法通常是通过地震波的谱信息来计算。然而，目前依赖于地震波谱信息分析的能量吸收分析方法中很少考虑反射系数的影响。实际上，在许多沉积地层中，地震波反射谱的主要成分来自相互干扰的近距离反射，即众所周知的调谐效果[22]，地震数据谱的形状(包括幅度、频率和相位)不可避免地受到调谐效应的影响[23-25]。为了抑制反射系数的影响，校正地震波频谱形状，给出更为准确可靠的衰减估计结果，本文引入变分模态分解(VMD)方法，结合平均对数谱距离(ALSD)发展一种自适应鲁棒的能量吸收分析方法。

1 基于VMD的抑制反射系数影响的频率衰减分析算法

能量吸收分析方法主要利用指数衰减函数exp(-af)在时频域中估计高频能量吸收的数量[14,17]：

X(t，f)=ke(-af)

(1)

式中：X(t，f)为地震波x(t)的时频分布；k为常数；f为频率；a为吸收系数或称作衰减梯度，被广泛应用于烃类指示。对式(1)两边取以e为底的对数，得：

lnX(t，f)=lnk-af

(2)

目前衰减梯度估计中，通常对式(2)直接利用两点斜率法或者最小二乘法估计衰减梯度。这里，为了抑制反射系数影响，校正地震波频谱，首先考虑地震波的平稳卷积模型：

x(t)=w(t)*r(t)

(3)

式中：w(t)为地震子波函数；r(t)为反射系数序列。

式(3)在频域中变为：

X(f)=W(f)R(f)

(4)

对式(4)两边取对数得：

lnX(f)=lnW(f)+lnR(f)

(5)

通常情况下，反射系数序列r(t)具有类似随机噪声特性[26]，相对于子波函数，它的频谱变化非常快，而子波函数相对于反射系数序列通常变化很慢，通常是一个低频信号。地震信号在对数域可以看作是低频子波函数序列和相对高频反射系数序列的和。

接下来，引入VMD方法进行分析。由于VMD分解是在时频域中进行处理的，因此首先利用具有较高分辨率的广义S变换将地震信号转换到时频域。此处VMD方法可以看作是一个自适应、多频带滤波器，VMD是一个约束变分问题，表示为[27]：

(6)

式中：u表示非线性和非平稳信号；ck表示第k个本征模态函数(IMF)；ωk是角频率。

ck在ωk周围几乎是紧凑的。通常，使用二次惩罚和拉格朗日乘数来解决问题式(6)。 通过采用乘数的交替方向方法(ADMM)，可以获得每个筛选过程中的每个模式和中心频率。VMD方法的更多细节及主要步骤详见文献[27]。经过VMD分解后，原始信号u被分解为一系列窄带本征模态函数ck的和：

(7)

式中：n为IMF的个数。

对于对数域地震信号lnX(f)，经过VMD分解后，可以得到一系列从低频到高频的IMF。从而，可以找到一个分离点，在分离点之前的信号主要反映低频的子波函数序列，分离点后的信号主要反映高频的反射系数序列。为此，引入ALSD作为有用的相似性度量以有效地区分子波函数主导的IMF和反射系数序列主导的IMF。ALSD是两个对数频谱之间的距离度量，被广泛应用于语音编码领域中用于测量量化或内插性能[28-30]。通过最大化ALSD来确定子波函数主导的IMF和反射系数序列lnR(ω)主导的IMF之间的分离点：

(8)

其中，对于两个IMF之间的ALSD定义为：

(9)

式中：m是ck的长度。然后，t1时刻地震子波函数W(ω)主导的子信号可以通过下式重构：

(10)

提取到地震子波函数主导的子信号后，结合式(2)进一步可以获得衰减梯度估计式：

(11)

式中：f1为地震谱峰值频率。考虑到处理后的窄频谱特性，f2按下式进行动态选取：

(12)

式中：fm为频谱主要范围截止点。

为了给出更为准确的衰减梯度估计值，此处利用最小二乘法拟合获得衰减梯度a。设拟合多项式φ(f)=a0+a1f，则最小二乘法是对于给定的数据y=s(f)的数据点(fi，lnW(t，fi)；i=1，2，…，m)求近似曲线φ(f)，且使φ(f)与原曲线y=s(f)偏差最小，可以表示为：

(13)

式中：δi为近似曲线在点(fi，ln(W(t，fi))处的偏差。

对式(13)右边求对于ai的偏导数，列出系数矩阵求解，确定拟合系数a1。则该时刻的衰减梯度a=|a1|。基于VMD的抑制反射系数影响的能量吸收分析算法的处理流程如图1所示。

图1 基于VMD的抑制反射系数影响的频率衰减分析算法的流程图

2 实际地震数据应用

利用四川盆地川西坳陷某气田的叠后偏移地震数据进行分析。该数据前期处理中经过了地表一致性振幅补偿提高振幅保真性的处理，但不包括自适应反褶积或时变频谱白化、几何扩散补偿、自动增益AGC等幅度补偿处理。该气田研究区河道砂发育，储层厚度较薄，物性较差且致密，平均孔隙度8.66%，平均渗透率0.21 mD。受非均质性的影响，不同砂组及砂体含气性差异较大，气水分布十分复杂，在构造高、低部位均有天然气及地层水分布。该研究区存在的一类储层从地球物理特征角度来看具有中等波阻抗、弱地震响应的特征，笔者认为其属于一类隐蔽储层，在叠后数据处理中传统方法较难识别。

2.1 过井道基于VMD的抑制反射系数影响的衰减特征及算法分析

表1 IMF分量的ALSD

图2 含气水井A的测井解释结果

图3 过井道衰减分析

2.2 二维地震剖面处理

图4为研究区一条过含气井B的叠后偏移地震剖面。含气井B为斜井，储层位置如图4中椭圆区域所示，目标层地震响应较弱，椭圆区域地震信号与周围地震信号差异不明显，含气性检测困难，属于隐蔽储层。含气井B的测井解释结果如图5所示。该剖面基于VMD的抑制反射系数影响的衰减梯度剖面如图6(a)所示，可以看出，椭圆区域有强振幅异常响应，由于该段内岩性稳定，排除岩性、地层等因素影响，认为该强振幅异常特征表明该椭圆区域存在烃类，符合测井解释结果的含气性检测结果。而在基于小波变换的衰减梯度剖面(见图6(b))中，在相同椭圆区域位置并没有强振幅异常响应，说明该算法不能很好地给出符合测井解释结果的含气性检测结果，表明针对隐蔽储层的含气性检测问题，传统方法无法有效检测到微弱地震响应信号的强振幅异常特征。

图4 过含气井B的叠后偏移地震剖面

图5 含气井B的测井解释结果

图6 基于VMD抑制反射系数影响的衰减梯度剖面与基于小波变换的衰减梯度剖面对比

3 讨论

3.1 VMD分解次数选取

基于VMD的抑制反射系数影响的频率衰减分析方法中，重构抑制反射系数影响的频谱是该算法的一个关键步骤。合适的VMD分解次数影响着重构抑制反射系数影响的频谱的准确性。尽管通常情况下，选择IMF的数量(即分解级数)是VMD的关键问题，VMD分解次数设置过大，虽然可以尽量地避免子波信息和反射系数信息的重叠，但是会导致分解产生的IMF分量过分解，捕获额外的噪声或者导致模态重复；而如果VMD分解次数设置过小，会导致数据分割不足，导致子波信息和反射系数信息重叠较多，产生模态混叠。目前虽然关于VMD分解次数选取有很多研究并发展了很多方法[31-34]，但在本次研究的算法中选择并不困难。在本次研究的算法中，主要通过计算各个IMF分量的ALSD来确定子波函数主导的IMF和反射系数序列主导的IMF之间的分离点，而不是获得每个IMF的精确表达，因此，由过多模式引起的模式重复的传统难题在本次研究中不需要考虑，只需要避免生成太少的模式导致的模态混叠问题，即子波信息和反射系数信息重叠太多的情况。这种情况很容易解决，对于实际地震数据，通常使用过井道结合测井解释结果来获得适当的VMD分解级数。此外，需要说明的是，尽管本次研究使用VMD算法来分离子波信息和反射系数信息，但是实际上子波信息和反射系数信息在对数谱中并不是完全分开的，由于筛选出的信息主要是子波信息，且反射系数信息相对子波信息频率更高，体现在频谱上主要是高频细节波动信息，而子波信息决定了总体频谱曲线的趋势，个别细节微小波动信息不会影响到总体频谱曲线趋势，从而对衰减梯度计算结果的影响很小。图6(a)的结果也表明了该方法可以给出准确度很高的衰减梯度估计值。

此外，由于VMD的强局域分解能力和抗噪声性能[27]，噪声信息往往会被提取为单独的IMF信号，因此本文的方法也可以适用于低信噪比地震数据的处理。

3.2 衰减梯度估计的高频优势频段选取

基于VMD的抑制反射系数影响的频率衰减分析方法中，影响结果精度的另一个因素是高频优势频段的选取。传统衰减梯度估计算法中，通常选取总能量65%到总能量85%的频段进行衰减梯度计算。本次研究中，针对采用了广义S变换后时频谱分辨率较高、频带较窄，且经过IMF选择重构后的抑制反射系数频谱频带又变窄的特性，采用了动态频段选取方式，频段选取方式如式(12)所示。该式中采取了高频优势频段最大频段计算范围为40 Hz的处理方式。在实际地震数据处理中，也可以根据实际地震数据特征，适当调整该参数。

3.3 含气性检测结果分析

含烃类区域通常在衰减梯度剖面中表现为强振幅异常特征，排除岩性、地层等其他因素影响，通常认为衰减梯度剖面中的强振幅异常特征可以给出含烃类信息的解释结果。本次研究中针对四川盆地川西坳陷某气田的叠后偏移地震数据进行分析，选取了隐蔽储层段进行研究。从过含气水井A的地震道分析及过含气井B的地震剖面分析中可以看到，常规基于小波变换的衰减梯度计算方法检测不到隐蔽储层中的烃类信息，不能给出很好的烃类解释结果，而基于VMD的抑制反射系数影响的衰减分析方法可以检测到这类弱含气响应信号的强振幅异常信息，给出符合测井解释结果的含烃类检测结果。

4 结束语

本文提出了一种基于VMD和ALSD的抑制反射系数影响的频率衰减分析算法，该算法避免了常规频率衰减分析算法中反射系数的影响，能实现更准确地估计频率衰减信息。为了抑制反射系数影响，该算法在时频域对地震道的对数谱进行VMD分解，结合ALSD判断子波函数主导的IMF和反射系数序列主导的IMF之间的分离点，通过重构获得抑制了反射系数影响的频谱。另一方面，为了保证提取的频率衰减信息的准确性，对重构的抑制了反射系数影响的频谱采用了动态选取频段的方式计算频率衰减信息，实现窄频带情况下高频衰减信息计算频段的优势选取。四川盆地川西坳陷某气田的叠后偏移地震数据处理结果表明，该算法可以有效识别出一些储层特征不明显的隐蔽储层，检测出弱烃类信号引起的能量异常，给出符合测井解释结果的烃类解释。