变分模态分解与包络导数算子结合的时频分析方法及溶洞储层预测

2021-05-15宋维琪陈礼豪陈俊安杨子鹏

石油地球物理勘探 2021年2期

武迪宋维琪刘军陈礼豪陈俊安杨子鹏

(①中国石油大学(华东)地球科学与技术学院，山东青岛 266580；②中国石化西北油田分公司勘探开发研究院，新疆乌鲁木齐 830011)

0 引言

深层缝洞型碳酸盐岩储层中缝洞空间分布散乱，不同缝洞的规模、形态和内部结构等特征差异明显，并且由于储层埋深较大，造成地震信号信噪比低。因此由缝洞引起的“串珠状”地震响应与背景分离度差，储层识别困难[1-2]。碳酸盐岩地层中的裂缝、溶洞往往是良好的油气聚集空间[3-4]，因此缝洞体储层预测具有十分重要的意义。

油气的存在导致地震信号频率和能量异常，通过谱分解技术，容易发现隐藏的特定频段的异常信息，目前已经被广泛用于储层预测、烃类检测等方面[5-8]。近年来，在经典谱分析方法的基础上，涌现出经验小波变换[9]、稀疏时频分析[10]、同步压缩变换[11]等多种有效算法。作为一种时频分析手段，联合经验模态分解(Empirical Mode Decomposition，EMD)与希尔伯特变换的希尔伯特—黄变换(HHT)方法由Huang等[12]提出，能够自适应地将信号分解为一系列的固有模态函数(Intrinsic Mode Function，IMF)，继而分析各IMF的瞬时频率、瞬时振幅等[13]。EMD具有更高的时频分辨能力，但也存在模态混叠、端点效应等问题[14]。为了解决模态混叠问题，人们相继提出引入噪声辅助的集合经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)、完备集合经验模态分解(Complete Ensemble Empirical Mode Decomposition，CEEMD)等方法[15-16]，在一定程度上消除了模态混叠，已被用于地震资料处理[17]及烃类检测[18-20]，但由于其本质上缺乏完备的理论支撑，在实际应用中存在一定局限性。Dragomiretskiy等[21]提出了变分模态分解(Variational Mode Decomposition，VMD)方法，具有完善的理论基础[22]，能够将多分量信号自适应、非递归地分解为一系列具有带限性质的IMF，有效控制了模态混叠问题。由于VMD的本质是维纳滤波的推广，因而具有良好的抗噪性，在地震资料去噪[23-24]及时频分析领域[25]取得了良好的效果。

Teager等[26]提出了Teager能量算子，Kaiser[27]给出了其离散形式——Teager-Kaiser(TK)算子，因具有良好的局部分析能力及考虑了频率特性，de Matos等[28]最早将其引入地震勘探邻域，用于分析碳酸盐岩储层。陈学华等[29]、唐湘蓉等[30]在时频域进行储层预测及流体识别。Xue等[31]在EMD的基础上，利用TK算子代替希尔伯特变换获得瞬时属性进行联合时频分析，解决了希尔伯特变换的端点效应问题，然而TK算子对噪声敏感的特性限制了其应用。O′Toole等[32]提出了具有良好非负特性的包络导数算子(Envelope Derivative Operator，EDO)，与TK算子相比，同为具有良好的追踪信号瞬时变化能力的非线性能量算子，但EDO抗噪性和稳定性更好[33]，已经被用于轴承故障检测等方面[34]。在本文中，为了满足该类非线性能量算子的应用条件，利用VMD技术将地震信号分解为包含低频、高频有效信息的IMF，再引入EDO获取瞬时振幅、瞬时频率，克服了希尔伯特变换的端点效应且具有更高的时频分辨率，最终形成一种高精度时频分析方法，并应用于中国西部N区缝洞型碳酸盐岩储层预测，取得了良好的效果。

1 方法原理

1.1 VMD

VMD的目标是将输入信号分解为一系列子信号，分解后得到的模态分量具有稀疏特性并且可以重建原始信号。与EMD相比，VMD具有完备的理论基础，从本质上看，VMD是维纳滤波的推广。

VMD非递归地将多分量信号自适应地分解为具有特定稀疏属性的带限模态分量，首先定义具有中心频率的有限带宽的IMF

u(t)=A(t)cosφ(t)

(1)

式中：信号的相位φ(t)需满足其对时间t的一阶偏导数φ′(t)≥0；A(t)为包络。

定义了式(1)的IMF之后，VMD过程包括变分问题的构造和求解，若将分析信号分解成K个IMF，则对应的约束变分模型为

(2)

为了求解式(2)的最优解，引入二次惩罚因子α和Lagrange乘法算子λ(t)，则约束变分问题转换为非约束变分问题。增广Lagrange公式为

(3)

(4)

式中n为迭代次数。对应的中心角频率为

(5)

通过

(6)

(7)

则迭代结束，得到K个IMF。

由VMD分解信号时，需要预先设定IMF的个数K，在测试模拟信号时发现K取l+1(l为合成信号分量的个数)时效果最好。由于实际信号复杂多变，吴文轩等[35]、刘尚坤[36]分别利用峭度及互信息准则等方法确定K值，文中利用经验方法，确定当K=3时可取得最佳的分解效果。

1.2 包络导数算子

1.2.1 TK算子

对于连续信号x(t)，TK算子被定义为二阶微分方程的形式[26]

(8)

x(m)=a(m)cosφ(m)

(9)

式中a(m)为调幅信号，m为采样点号。x(m)的离散形式的TK算子为[27]

ψd[x(m)]=x2(m)-x(m-1)x(m+1)

(10)

TK算子是仅利用差分运算计算瞬时能量的非线性算子，考虑了信号的频率特性，最明显的优势是良好的局域特性及计算的简洁、高效性。

1.2.2 EDO

对于信号瞬时能量的求取，典型的方法是以振幅值的平方进行量化，或者以包络的形式

(11)

表征。式中H[·]为希尔伯特变换。为了引入频率域的信息，可以在式(11)中以导函数的形式引入加权滤波器。根据傅里叶变换(FT)的性质，有

(12)

则包络导数算子为[32]

(13)

式(13)结合了信号的频率变化信息与随时间变化的包络，称为包络导数算子(EDO)。EDO的定义形式与TK算子较相似，只有第二项存在区别。O′Toole等[32]详细对比了TK算子、EDO用于x(t)=Acos(ω0t+φ)(A为振幅，ω0为初始角频率，φ为相位)或x(t)=Aertcos(ω0t+φ)(r为调节系数，用于描述振幅的变化)形式的简单信号时的特性，均取得较好效果。对于线性组合信号y(t)=x1(t)+x2(t)

其中

(14)

将TK算子、EDO应用于y(t)，分别得到

ψ[y(t)]=ψ[x1(t)]+ψ[x2(t)]+

cos[(ω1-ω2)t+φ1-φ2]}

(15)

Γ[y(t)]=Γ[x1(t)]+Γ[x2(t)]+

a[sin(ω1t+φ1)sin(ω2t+φ2)+

cos(ω1t+φ1)cos(ω2t+φ2)]

=Γ[x1(t)]+Γ[x2(t)]+

acos[(ω1-ω2)t+φ1-φ2]

(16)

式中a=2A1A2ω1ω2。由式(15)、式(16)可见：Γ[y(t)]包含一个等于信号y(t)两个分量频率之差的项；ψ[y(t)]则包含了额外的附加项，这将在TK能量计算时出现负值。图1为信号y(t)及TK算子、EDO能量计算结果。由图可见：对于能量的变化特征，TK算子、EDO均取得了良好的追踪结果，但在部分极值点处TK算子出现负值，这将在进一步分离能量时出现奇异值；EDO则具有良好的非负特性。

图1 信号y(t)(上)及TK算子、EDO能量(下)计算结果

TK算子采用向前差分法获得离散形式，为了得到EDO的离散形式，定义以下中心差分形式

(17)

因此，离散形式的EDO为

h2(m+1)+h2(m-1)]+

h(m+1)+h(m-1)]

(18)

式中h(·)=Η[x(·)]。

1.2.3 基于VMD-EDO的时频分析方法

基于VMD-EDO的时频分析方法如图2所示。

图2 基于VMD-EDO的谱分解算法

首先，将地震信号进行VMD，原始数据中的有效信息将主要由一个或几个IMF所反映，选取包含储层信息最多的IMF进行分析。

然后，采用能量分离(ESA)算法获得各分量信号EDO能量的瞬时频率和瞬时振幅[37]。由于ESA算法只适用于单分量信号，不能直接用于地震信号，而VMD方法将地震信号分解为一系列窄带信号，近似满足非线性能量算子的计算条件，可获得EDO能量。通过三点对称差分运算，得

(19)

(20)

这里，瞬时频率ω(m)和瞬时振幅A(m)都是时间的函数，因此可以定义一个三维空间[t，A(t)，ω(t)][38]

(21)

最终可以获得IMF的时频分布，进一步可以分离得到不同频段的高精度瞬时谱。

2 模型验证

为了验证VMD-EDO时频分析方法的效果，建立了地质模型(图3a、表1)。采用主频30Hz的雷克

图3 基于模型的VMD-EDO谱分解

子波作为震源进行模拟，并添加7%的随机噪声。基于模型的VMD-EDO谱分解结果表明：在合成地震记录上缝洞体呈“串珠状”响应(图3b)；在VMD-EDO低频剖面上缝洞体呈强能量体特征(图3c)；相对于地层信息，在VMD-EDO高频剖面上缝洞体能量衰减明显(图3d)。因此VMD-EDO谱分解结果总体体现了“低频能量加强、高频能量衰减”的特征，证明了方法的有效性。

表1 模型参数

3 实际资料应用分析

为进一步验证VMD-EDO时频分析方法的效果，选取中国西部N区的叠后地震资料进行测试。N区位于顺托果勒低隆的北部，处于阿瓦提、满加尔坳陷与沙雅隆起的结合部，紧邻南部海相烃源岩灶，同时发育本地寒武系烃源岩，构造位置有利，油源充足，是油气长期运移、聚集的有利区。目前研究表明，N区奥陶系油气成藏条件优越，储层地震响应主要以同相轴错断、异常弯曲和“串珠状”为主。由于缝洞储层自身的规模、形态和内部结构差异明显以及埋深较大，对于时频分析精度提出了更高的要求。

图4为A井井旁地震道CDP105及VMD结果，图5为IMF2振幅谱及由TK算子、EDO得到的瞬时振幅。由图可见：①IMF1为超低频信息(图4b上)，而中低频及高频有效信息主要分别包含在IMF2(图4b中)与IMF3(图4b下)中。②整体来看，对于4600～4700ms层段的油气异常特征，分别采用TK算子(图5b上)与EDO(图5b下)得到的瞬时振幅相近，并且均避免了希尔伯特变换的端点问题。③由于TK能量计算时出现负值，在部分极值点附近瞬时振幅出现部分奇异值，需要采用取绝对值才能避免此问题(图5b上)；EDO计算结果更平滑，且EDO良好的非负特性有效避免了奇异值(图5b下)。

图6为A井井旁道 CDP109时频谱。由图可见，有效信息主要分布在15～40Hz频段，在4650ms附近出现强能量异常，VMD-EDO时频谱(图6a)的时频分辨率高于连续小波变换(CWT)时频谱(图6b)。

撒利明等[39]采用“低频共振、高频衰减”特征检测油气，薛雅娟等[18]分别讨论了低频(14～18Hz)、高频(26～30Hz)的时频特征，并预测了鄂尔多斯盆地海相碳酸盐岩储层(频带范围为12～30Hz)。图7为过A井地震剖面及VMD结果。由图可见：在过A井地震剖面的CDP105、4650ms附近钻遇油气，呈小“串珠状”地震响应(图7a)；IMF2(图7c)、IMF3(图7d)分别反映了中低频、高频有效信息。图8、图9分别为VMD-EDO、CWT-EDO的低频、高频剖面。由图可见：①在含油气区域(红圈位置)，两种方法都取得了良好的检测结果，具体表现为：在低频剖面上呈不连续的强亮点及周边存在弱能量团(图8a、图9a)；随着频率升高，在高频剖面上含油气区域(红圈位置)的能量衰减明显(图8b、图9b)。②VMD-EDO方法检测结果的分辨率更高，能更好地区分不同深度的能量异常。图10为对应图7c、图7d的VMD-EDO低频、高频目的层沿层切片。由图可见，A井、B井均为良好的产油井，两口井所在区域均明显表现为低频段强能量异常(图10a)、高频能量衰减(图10b)特征，进一步验证了方法的有效性。因此，在缝洞型碳酸盐岩储层中，基于VMD-EDO的高精度时频分析算法的储层预测效果较好。