声波测井作为三大测井方法之一, 是油气勘探与开发的重要技术。阵列声波测井记录的全波列主要包含纵波、横波、伪瑞利波及斯通利波, 其中伪瑞利波由于受较多因素影响而难以实际应用^[1-2]。模式波的速度、幅度衰减和频率变化分别与地层运动学、动力学、吸收滤波性质关系密切^[3]。因此, 通过研究全波列中各种模式波的到时、幅度衰减、频谱、相位等信息可以直接或间接地获取地层的岩石类型、储集性、流体性质等信息。全波列数据蕴含的地层信息十分丰富^[4-5], 且能够反映井眼周围和距离井眼较远地层的力学和流体特性^[6], 因此阵列声波测井数据的处理和信息提取显得尤为重要。

阵列声波测井数据是典型的非线性、非平稳信号, 用传统的傅里叶变换进行处理存在很多局限^[7], 且仅能获得信号在时域或频域的全局特征。时频分析方法可解决时域和频域的局部化矛盾, 通过时频分析方法可以获知信号的某个频率成分出现在什么时间。该方法能将信号的时间、频率和幅度(能量)这3个与地层性质关系密切的信息联系起来, 得到幅度(能量)在时-频坐标系中的分布。

目前时频分析方法已被广泛应用于处理阵列声波测井数据。田鑫等^[8]、陶钧^[9]、刘娇等^[10]用短时傅里叶变换提取了阵列声波测井信号中各模式波的时差; 王祝文等^[11]用重排谱图方法处理阵列声波测井资料来识别断裂、构造; 王祝文等^[12]分析了不同岩性构造破碎带的Choi-Williams能量分布特征; 张学涛等^[13]根据阵列声波信号的Choi-Williams分布特征来区分油、水层; 高云娇^[14]研究了不同性质储集层的Page分布特征; 杨闯等^[15]通过Choi-Williams分布判断裂缝地层的破碎程度。随着时频分析方法在阵列声波测井信号信息提取方面的深入应用, 发展出联合两种时频分布处理信号的方法, 主要是联合经验模态分解(EMD)或分数阶傅里叶变换和其它时频分布, 取得了良好的效果。王祝文等^[16]联合使用EMD和平滑伪Wigner-Ville分布提取阵列声波信号所蕴含的裂缝特性; WANG等^[17]、王晓丽^[18]联合分数阶傅里叶变换和平滑伪Wigner-Ville分布处理阵列声波测井信号; XIANG等^[19]与向旻等^[20]联合分数阶傅里叶变换和Choi-Williams分布、Born-Jordan分布提取阵列声波测井信号信息; 徐方慧等^[21]联合EMD和傅里叶变换分析了火成岩裂缝地层的特征。

各种时频分析方法特性不同, 处理效果也存在差异。EMD可自适应地将复杂信号分解为有限个固有模态函数(IMF), 分析IMF分量的时频分布可获得信号不同分辨尺度的信息。本文对比了基于EMD的短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布及重排谱图这4种时频分析方法, 首先对单点实测阵列声波测井全波列进行处理, 然后利用时频分布的边缘特性处理井段数据, 主要从时间分辨率、频率分辨率和交叉项干扰几方面比较分析4种方法优缺点, 同时研究干层、水层、油水同层、油层阵列声波测井信号的时频特征。

1 方法原理 1.1 经验模态分解(EMD)

经验模态分解(EMD)由HUANG^[22]于1998年提出, 可以自适应地将非平稳信号分解为若干个具有实际意义的固有模态函数(IMF)^[23]。固有模态函数需满足以下2点^[24]:①在整个信号中, 零值点与极值点的个数相等或最多相差一个; ②信号上任意一个数据点, 由所有的局部极大值点组成的上包络线、所有局部极小值点组成的下包络线的均值等于零, 即关于时间轴对称。IMF的优点在于:任何复杂的信号都可以表示为若干个IMF的叠加, 而每阶IMF都代表了信号中一种简单频率的振荡。

信号x(t)经EMD可以表示为^[25]:

$ x\left( t \right) = \sum\limits_{q = 1}^N {{c_q}} \left( t \right) + {r_N}\left( t \right) $

(1)

式中:N为分解次数; c_q(t)为第q个IMF分量; r_N(t)为残余分量, 表示原始信号整体的振动趋势。

1.2 短时傅里叶变换

信号x(t)∈L²(R)的短时傅里叶变换(STFT)与其逆变换可表示为^[26]:

$ {\rm{STF}}{{\rm{T}}_x}\left( {t,f} \right) = \int\limits_{ - \infty }^{ + \infty } {x(\tau )h(t - \tau ){{\rm{e}}^{ - {\rm{i}}2{\rm{ \mathsf{ π} }}ft}}{\rm{d}}t} $

(2)

$ x(t) = \frac{1}{{2{\rm{ \mathsf{ π} }}h(0)}}\int\limits_{ - \infty }^{ + \infty } {{\rm{STF}}{{\rm{T}}_x}(t,f){{\rm{e}}^{{\rm{i}}2{\rm{ \mathsf{ π} }}ft}}{\rm{d}}f} $

(3)

式中:t为时间; f为频率; τ为时移; i为虚数单位; h(t)为对称的时间宽度很短的窗函数。

将(2)式离散化可得^[27]:

$ {\rm{STF}}{{\rm{T}}_x}\left( {n,f} \right) = \sum\limits_{m = - \infty }^{ + \infty } x (m)h(n - m){{\rm{e}}^{ - {\rm{i}}2{\rm{ \mathsf{ π} }}fm}} $

(4)

式中:h(m)为实数窗序列。

对于多分量信号$x(t) = \sum\limits_{k = 1}^M {{x_k}} (t)$, 其短时傅里叶变换为:

$ {\rm{STF}}{{\rm{T}}_x}\left( {t,f} \right) = \sum\limits_{k = 1}^M {STF{T_{xk}}\left( {t,f} \right)} $

(5)

式中:M为信号分量总数; x_k(t)为信号x(t)的第k个分量; STFT_xk(t, f)为第k个分量的短时傅里叶变换。(5)式表明短时傅里叶变换处理结果无交叉项。

根据(4)式可知, 短时傅里叶变换的本质是加窗的傅里叶变换:在时域用一个固定的时间窗来截取一小段信号, 假定时窗内的信号为平稳信号并对其做傅里叶变换, 通过连续不断移动窗函数并做傅里叶变换得到该信号的短时傅里叶变换^[28]。因此, 短时傅里叶变换的性能在很大程度上取决于窗函数的类型和窗口的宽度。不同的窗函数其功能特性有很大的差别。窗函数类型选择的总原则是使其频谱中主瓣宽度尽量窄, 旁瓣高度尽量小, 旁瓣峰值衰减尽量快^[29]。另外, 要得到高时间分辨率就要求加短时窗, 而要得到高频率分辨率则要求加长时窗。根据Heisenberg不确定性原理, 短时傅里叶变换无法在时域和频域同时取得最佳分辨率。

短时傅里叶变换对整个信号加不变时窗, 更适用于处理局部平稳长度大的非平稳信号, 其主要优点是易于实现、计算快捷^[30]、无交叉项; 但选定窗函数及窗口长度就意味着确定了时间分辨率和频率分辨率, 且时间分辨率、频率分辨率受Heisenberg不确定性原理制约, 难以适应信号频率变化的不同要求, 这是其最大不足。

1.3 Choi-Williams分布

20世纪60年代中期, COHEN^[31]提出用统一的表达式来表示众多的双线性时频分布, 后被称为Cohen类时频分布。信号x(t)的Cohen类时频分布为:

$ \begin{array}{*{20}{c}} {P(t,f) = \int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {x\left( {u + \frac{\tau }{2}} \right){x^*}\left( {u - \frac{\tau }{2}} \right)g(\tau ,v) \cdot } } } }\\ {{{\rm{e}}^{ - {\rm{i}}2{\rm{ \mathsf{ π} }}\left( {tv + f\tau - uv} \right)}}{\rm{d}}u{\rm{d}}v{\rm{d}}\tau } \end{array} $

(6)

式中:v为频移; u为时间积分变量; *表示取共轭; g(τ, v)为核函数, 不同的核函数对应着不同Cohen类时频分布。Cohen类时频分布的本质是以核函数加权的模糊函数的二维傅里叶变换。双线性时频分布的时间分辨率、频率分辨率不再受Heisenberg不确定性原理约束, 但缺点是计算相对线性时频分布复杂, 而且固有的不足就是存在交叉干扰项。多分量信号$x(t) = \sum\limits_{k = 1}^M {{x_k}} (t)$的Cohen类时频分布为:

$ \begin{array}{*{20}{c}} {{P_x}\left( {t,f} \right) = \sum\limits_{k = 1}^M {{P_{xk}}} \left( {t,f} \right) + 2\sum\limits_{k = 1,k \ne {k_1}}^M {\mathop \sum \limits_{{k_1} = 1}^M } \cdot }\\ {{\mathop{\rm Re}\nolimits} \left[ {{P_{{x_k}{x_{k1}}}}\left( {t,f} \right)} \right]} \end{array} $

(7)

式中:$\sum\limits_{k = 1}^M {{P_{{x_k}}}} (t, f)$为自主项; $2\sum\limits_{k = 1, k \ne {k_1}}^M {\sum\limits_{{k_1} = 1}^M {{\mathop{\rm Re}\nolimits} } } \left[ {{P_{{x_k}{x_k}_1}}(t, f)} \right]$为交叉项; Re[·]表示取实部。(7)式表明交叉项的幅度可达自主项的两倍, 会严重影响有用信号时频特征的识别。交叉项主要通过设计核函数来抑制。

在(6)式中g(τ, v)=1, 取核函数, 即得到Wigner-Ville分布。Wigner-Ville分布的优点是时频分辨率很高, 但交叉项干扰十分严重。CHOI等^[32]于1989年提出Choi-Williams分布, 提高时间和频率分辨率, 同时抑制交叉项。(6)式中取核函数g(τ, v)=${{\rm{e}}^{ - \frac{{v2\tau 2}}{\sigma }}}$即为Choi-Williams分布:

$ \begin{array}{*{20}{c}} {{\rm{CW}}{{\rm{D}}_x}(t,f) = \iint{\sqrt{\frac{\text{ }\!\!{\rm{\pi}}\!\!\text{ }\sigma }{{{\tau }^{2}}}}}x\left( {u + \frac{\tau }{2}} \right) \cdot }\\ {{x^*}\left( {u - \frac{\tau }{2}} \right){{\rm{e}}^{ - \frac{{{{\rm{ \mathsf{ π} }}^2}\sigma {{\left( {u - t} \right)}^2}}}{{4{\tau ^2} - {\rm{i}}2{\rm{ \mathsf{ π} }}f\tau }}}}{\rm{d}}u{\rm{d}}\tau } \end{array} $

(8)

式中:σ(σ>0)是缩放因子, 其降低自相关分辨率用于抑制交叉项。为得到高自主项分辨率, σ值应该大; 而为减少交叉项的影响, σ值应该小^[33]。σ的选择应针对每个具体应用进行优化, 针对幅度和频率变化快的信号应该取较大的σ, LIN^[34]认为σ选择为0.1~10.0时应用效果较好。

1.4 平滑伪Wigner-Ville分布

平滑伪Wigner-Ville分布是在Wigner-Ville分布的基础上再加入一个时间窗函数和频率窗函数, 以此来抑制交叉项。信号x(t)的平滑伪Wigner-Ville分布(SPWD)为:

$ \begin{array}{*{20}{c}} {{\rm{SPWD}}\left( {t,f} \right) = \int\limits_{ - \infty }^{ + \infty } {{h_2}(\tau )} \int\limits_{ - \infty }^{ + \infty } {{h_1}} (u - t)x\left( {u + \frac{\tau }{2}} \right) \cdot } \\ {{x^*}\left( {u - \frac{\tau }{2}} \right){\text{d}}u{{\text{e}}^{ - {\text{i}}2{\rm{\pi }}f\tau }}{\text{d}}\tau } \end{array} $

(9)

式中:h₁(u)为时间窗函数; h₂(τ)为频率窗函数。

1.5 重排谱图

谱图的本质是信号x(t)的短时傅里叶变换模的平方^[35]:

$ \begin{gathered} {S_x}\left( {t,f} \right) = {\left| {{\rm{STF}}{{\rm{T}}_x}(t,f)} \right|^2} = \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left| {\int\limits_{ - \infty }^{ + \infty } x (\tau )h(\tau - t){{\text{e}}^{ - {\text{i}}2{\text{\pi }}ft}}{\text{d}}t} \right|^2} \hfill \\ \end{gathered} $

(10)

谱图是一种实值、非负的双线性分布, 其时频分辨率和短时傅里叶变换一样受到约束, 且存在干扰项。最初提出重排的目的是抑制交叉项, 从而改良谱图的效果^[36]。谱图可看作是信号x(t)的Wigner-Ville分布和分析窗的Wigner-Ville分布的二维卷积:

$ {S_x}\left( {t,f} \right) = \int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {{W_x}(u,\xi ){W_h}(t - u,f - \xi ){\text{d}}u{\text{d}}\xi } } $

(11)

式中:ξ为频率积分变量; W_h是窗函数h(t)的Wigner-Ville分布。

根据(11)式, 信号的Wigner-Ville分布的加权平均值被分配到点(t, f)的邻域。重排原理的关键思想是:重心更能代表信号的局部能量分布, 应该将这些值分配到时频域的重心而不是其几何中心。

重排指的是将任意一点(t, f)处计算得到的谱图值分配到它附近信号能量的重心(, ):

$ \hat t = \frac{{\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } u } {W_x}(u,\xi ){W_h}(t - u,f - \xi ){\rm{d}}u{\rm{d}}\xi }}{{\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{W_x}} } (\tau ,\xi ){W_h}(t - u,f - \xi ){\rm{d}}\tau {\rm{d}}\xi }} $

(12)

$ \hat f = \frac{{\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } v } {W_x}(u,v){W_h}(t - u,f - \xi ){\rm{d}}u{\rm{d}}\xi }}{{\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{W_x}} } (u,v){W_h}(t - u,f - \xi ){\rm{d}}u{\rm{d}}\xi }} $

(13)

重排后的谱图在任何一点(t′, f′)处的值等于所有重排至该点的谱图值之和。重排后的谱图为:

$ S_x^{(r)}\left( {{t^\prime },{f^\prime }} \right) = \int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {{S_x}(t,f)\delta \left( {{t^\prime } - \hat t} \right)\delta \left( {{f^\prime } - \hat f} \right){\rm{d}}t{\rm{d}}f} } $

(14)

这种新分布的一个非常有用的特性是它不仅利用了短时傅里叶变换的幅度信息, 而且利用了短时傅里叶变换的相位信息^[37]。这一点可以从下面的重排算子中看出:

$ \hat t = - \frac{{{\rm{d}}{\varphi _x}(t,f)}}{{{\rm{d}}f}} $

(15)

$ \hat f = f + \frac{{{\rm{d}}{\varphi _x}(t,f)}}{{{\rm{d}}t}} $

(16)

式中:φ_x(t, f)是信号x(t)的短时傅里叶变换的相位, 即φ_x(t, f)=arg[STFT_x(t, f)]。实际应用中, 这两个公式由于计算复杂, 须用下式代替:

$ \hat t = t -{\rm{ Re}}\left\{ {\frac{{{\rm{STF}}{{\rm{T}}_x}\left( {t,f;{T_h}(t)} \right){\rm{STFT}}_x^*(t,f;h)}}{{{{\left| {{\rm{STF}}{{\rm{T}}_x}(t,f;h)} \right|}^2}}}} \right\} $

(17)

$ \hat f = f - {\rm{Im}}\left\{ {\frac{{{\rm{STF}}{{\rm{T}}_x}\left[ {t,f;{D_h}(t)} \right]{\rm{STFT}}_x^*(t,f;h)}}{{{{\left| {{\rm{STF}}{{\rm{T}}_x}(t,f;h)} \right|}^2}}}} \right\} $

(18)

式中:Im[·]表示取虚部; STFT_x(t, f; h)表示信号x(t)加窗函数h(t)得到的短时傅里叶变换结果; T_h(t)=t×h(t)、D_h(t)=dh/dt(t)是平滑窗函数。

重排谱图不满足双线性性质, 但仍保留时移和频移不变性、非负性、能量守恒性的特点^[38]。重排后的谱图不仅能减小交叉项影响, 也能增强时频聚焦性^[39]。重排谱图与短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布方法的本质、适用条件、优缺点总结见表 1。

1.6 时频分布的边缘特性

一个时频分布的边缘特性能提供信号的重要信息。时间边缘的定义是:

$ {m_f}(t) = \int\limits_{ - \infty }^{ + \infty } {{\rm{TFR}}(t,f){\rm{d}}f} $

(19)

式中:TFR(t, f)为信号的时频分布。

频率边缘的定义是:

$ {m_t}(f) = \int\limits_{ - \infty }^{ + \infty } {{\rm{TFR}}(t,f){\rm{d}}t} $

(20)

一个时频分布的时间边缘的物理意义是信号的瞬时能量, 频率边缘的物理意义是信号的能量谱密度。

2 时频分析方法对比

图 1所示为一个原始的阵列声波测井全波列、经EMD得到前3个IMF分量和对应的频谱图。本文所有图对应的幅度(能量)均为相对刻度。全波列(图 1a)主要由纵波、横波、伪瑞利波、斯通利波组成; 其频谱(图 1b)较复杂, 高频成分的频率主要介于7.0~11.0kHz、低频成分的频率小于5.0kHz。EMD分解后得到的各IMF分量的频谱成分则比较简单:从IMF1到IMF3其频率是逐渐降低的, IMF1(图 1d)主要是大于7.0kHz的成分, IMF2(图 1f)主要是2.5~5.0kHz的成分, IMF3(图 1h)主要是小于3.0kHz的成分。经过对比可以发现, 各IMF分量的频谱分别与原始波列的频谱的不同区域对应性良好。

为了对比短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布、重排谱图的处理效果, 利用这4种时频分析方法制作原始波列的二维等值线图(图 2)。短时傅里叶变换表征信号的幅度在时频面上的分布情况, 其它3种时频分析方法则表征信号的能量在时频面上的分布情况, 为了方便对比, 将短时傅里叶变换结果取平方转换为能量分布。从图 2a短时傅里叶变换处理结果可以看出, 纵波出现在0.98~1.21ms, 主频约为9.4kHz; 横波出现在1.64~2.15ms, 主频约为8.8kHz, 波峰时间为1.90ms; 斯通利波出现在2.61~3.95ms, 主频为2.6kHz, 波峰时间为2.90ms。从图 2b Choi-Williams分布处理结果可以看出, 纵波出现在1.03~1.16ms, 主频约为9.3kHz; 横波出现在1.92~2.11ms, 主频约为8.9kHz, 波峰时间为1.88ms; 斯通利波出现在2.60~3.75ms, 主频为2.6kHz, 波峰时间为2.90ms。从图 2c平滑伪Wigner-Ville分布处理结果可以看出, 纵波出现在0.98~1.33ms, 主频约为9.3kHz; 横波出现在1.64~2.17ms, 主频约为8.8kHz, 波峰时间为1.88ms; 斯通利波出现在2.50~3.80ms, 主频为2.6kHz, 波峰时间为2.90ms。图 2d重排谱图显示纵波出现在1.04~1.17ms, 主频约为9.4kHz; 横波出现在1.71~2.05ms, 主频约为8.8kHz, 波峰时间为1.88ms; 斯通利波出现在2.65~3.73ms, 主频为2.9kHz, 波峰时间为2.87ms。短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布处理原始波列获得纵波、横波、斯通利波相应的出现时间、主频、波峰时间参数基本一致; 重排谱图的纵波、横波的出现时间、主频、波峰时间参数与其它3种方法结果基本相同, 但斯通利波的主频、波峰时间参数与其它3种方法结果相差较大。这是由于短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布是将时频分析值分配到时频域的几何中心, 而重排谱图则是将其分配到时频域的重心。对比图 2a、图 2b、图 2c、图 2d可以看出, Choi-Williams分布、平滑伪Wigner-Ville分布的时频分辨率较短时傅里叶变换高一些, 但相差不大; 重排谱图的时频聚焦性、时频分辨率在4种方法中最高, 有利于精准地确定波峰时间和主频。

由图 1可知, 原始波列经EMD分解得到的IMF分量具有不同的分辨尺度。为明确各IMF分量的模式波成分, 同时从不同分辨尺度上对比这几种时频分析方法, 对前3个IMF分量进行4种时频分析, 结果见图 3。分析IMF1分量的时频分析结果(图 3a、图 3d、图 3g、图 3j)可知, IMF1分量主要包含纵波、横波的信息; 分析IMF2分量的时频分析结果(图 3b、图 3e、图 3h、图 3k)可知, IMF2分量主要包含高频斯通利波、伪瑞利波的信息; 分析IMF3分量的时频分析结果(图 3c、图 3f、图 3i、图 3l)可知, IMF3分量主要包含低频斯通利波的信息。因此, 原始波列通过EMD分解可以自适应地分解为具有一定物理意义的IMF分量, 对不同的IMF分量分别进行时频分析处理可以得到不同尺度的信息。

分别对比IMF1分量(图 3a、图 3d、图 3g、图 3j)、IMF2分量(图 3b、图 3e、图 3h、图 3k)和IMF3分量(图 3c、图 3f、图 3i、图 3l)的4种时频分析图可以看出, 短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布的结果比较接近, 意味着时频分辨率相差不大。Choi-Williams分布的频率分辨率比短时傅里叶变换稍高一些, 平滑伪Wigner-Ville分布的时间分辨率和频率分辨率都比较高。在IMF2和IMF3分量的短时傅里叶变换(图 3b、图 3c)、Choi-Williams分布(图 3e、图 3f)、平滑伪Wigner-Ville分布(图 3h、图 3i)结果图中均只能看到主要有一个波峰, 重排谱图中却有几个时间、频率都很接近的波峰, 对应频谱图(图 1f、图 1h)中几个频率接近的峰, 说明重排谱图时频分辨率在4种方法中最高, 有利于更精准地定位波峰时间、主频, 也能更细致地刻画信号成分的变化。

从交叉项方面来看, 图 3显示4种方法均无交叉项影响纵波、横波、斯通利波的识别问题, 这是因为:短时傅里叶变换是线性时频分析方法, 本身无交叉项问题; Choi-Williams分布、平滑伪Wigner-Ville分布均为双线性时频分析方法, 分别通过核函数、加窗平滑取得较好的抑制交叉项效果; 重排谱图将时频分析值重排至时频域的重心, 减小交叉项同时改善了时频聚集性。

3 井段数据处理结果对比

图 4为不同流体性质储集层的常规测井曲线和阵列声波波形。4个深度段测井解释岩性均为杂色粗面岩, 其中4529~4534m常规测井解释为干层, 4568~4573m常规测井解释为水层, 4559~4563m常规测井解释为油水同层, 4388~4393m常规测井解释为油层。从图 4第5列可以看出, 干层的单极子波形中纵波、横波、斯通利波清晰可辨; 水层的单极子波形中纵波和横波衰减严重、斯通利波也有一定的衰减; 油水同层的单极子波形显示纵波衰减较严重、横波和斯通利波也有衰减; 油层的单极子波形显示斯通利波衰减十分严重。IMF1频谱图显示IMF1分量主要是大于8.0kHz的成分, 结合IMF1波形中IMF1分量主要分布在0~2.50ms判断, IMF1分量主要为频率高且到时靠前的纵波和横波; IMF2频谱图显示IMF2分量主要是2.0~6.0kHz的成分, 结合IMF2波形中IMF2分量主要分布在2.50ms之后, 但在0~2.50ms区域内也有分布判断, IMF2分量主要为伪瑞利波和高频斯通利波; IMF3频谱图显示IMF3分量主要是小于3.0kHz的成分, 结合IMF3波形中IMF3分量主要分布在2.50ms之后判断, IMF3分量主要为低频斯通利波, 验证了前文对单点声波全波列分析得到的结果。由此可见, 单极子波形通过EMD方法被分解为不同分辨尺度的IMF分量, 而IMF分量的波形图和频谱图可以分别从时域、频域表征信号, 但却无法指示某种频率成分对应的时间。由上文讨论可知, 时频分析方法可以解决这一问题, 而且通过提取时频分布的时间边缘和频率边缘特性可以对井段数据进行处理, 能够更直观地看到不同方法之间的差异和处理效果。

图 5为不同流体性质储集层声波全波列的IMF1分量边缘特征图。对于干层声波全波列的IMF1分量而言, 其纵波、横波幅度均较高, 因此IMF1的短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布、重排谱图 4种方法的时间边缘都能较清晰地显示两条明显的条带, 波峰时间和能量稳定, 分别对应纵波的波峰时间0.97ms和横波的波峰时间1.69ms。其中重排谱图的时间边缘中条带更窄、边界更清晰, 说明能量更集中、时间聚焦性更好, 因而波峰时间的定位更精准。而在IMF1的短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布3种方法的频率边缘图中, 主要表现为一条宽的条带, 且能量分布较均匀, 波峰能量不突出, 很难从中识别出纵波、横波, 因此无法获知纵波、横波的主频信息。而在IMF1的重排谱图的频率边缘中可以看到两条较明显的条带, 分别对应纵波的主频12kHz和横波的主频10kHz, 说明重排谱图的频率分辨率高于其它3种方法。下面将以干层的IMF1的边缘特性为基础来分析水层、油水同层及油层的IMF1分量的边缘特性。

水层IMF1分量时间边缘显示纵波和横波能量衰减明显, 由于纵波幅度偏低, 在短时傅里叶变换和平滑伪Wigner-Ville分布的时间边缘图中基本显示不出纵波, Choi-Williams分布、重排谱图的时间边缘图中均能看到代表纵波的条带, 但重排谱图的时间边缘显示更加清晰。纵波波峰时间延后至1.08ms, 横波波峰时间延后至1.93ms, 频率边缘则显示纵波和横波主频变化不明显、能量衰减严重。

纵波、横波的波峰时间变化与地层弹性常数和密度有关, 纵波的速度公式为v_P=$\sqrt {\left( {K + 4\mu /3} \right)/\rho } $, 横波的速度公式为v_S=$\sqrt {\mu /\rho } $, 其中K为体积模量, μ为剪切模量, ρ为地层密度。水层的体积模量、剪切模量均小于干层, 从补偿密度曲线可以看出, 水层密度小于干层, 但密度减小带来的影响远小于体积模量和剪切模量, 所以水层的纵波、横波速度低于干层, 表现为波峰时间延后。

纵波、横波能量的衰减主要是地层颗粒之间相对滑动、流体的吸收、声耦合率的改变几个因素的影响。水层密度降低, 地层颗粒相对滑动所需克服的摩擦阻力增大, 热损耗也增大。对于水及与水相近的液体介质, 声波的吸收系数为α=[2η(2πf)²]/(3ρ_fv_f³), 式中η, ρ_f, v_f分别为流体介质的粘滞系数、密度和声速, f为声波的频率。由此可见, 吸收系数与介质密度成反比关系, 与介质声波传播速度的三次方成反比。水层密度降低, 导致纵波、横波速度减小, 使吸收系数增大。同时, 地层含水会导致声阻抗发生变化, 声耦合率下降, 也会阻碍能量的传递。因此, 与干层相比, 水层的纵波、横波能量衰减严重。纵波和横波的主频均无明显变化则说明纵波、横波各频率成分能量的相对大小未显著改变。

油水同层的IMF1分量时间边缘显示纵波和横波能量均有衰减, 衰减的程度介于水层、油层之间。在短时傅里叶变换和平滑伪Wigner-Ville分布的时间边缘图中纵波显示不明显, Choi-Williams分布的时间边缘图中能看到代表纵波的条带, 而重排谱图的时间边缘图中纵波条带显示最清晰。纵波波峰时间滞后, 横波波峰时间延后, 频率边缘则显示纵波和横波主频无明显变化、能量发生衰减。

油层IMF1分量4种方法的时间边缘图差别较大, 其中Choi-Williams分布和平滑伪Wigner-Ville分布的结果比较一致, 可见代表纵波的条带, 但不够清晰, 同时在1.86ms处能量很高; 短时傅里叶变换与重排谱图的结果比较相近, 其中重排谱图的分辨率更高。如果只看Choi-Williams分布和平滑伪Wigner-Ville分布的结果会认为横波的波峰时间延后、能量升高, 而从补偿密度曲线可以看出油层密度略小于水层和干层, 克服地层阻力的热损耗应该更大, 加上流体的吸收和声耦合率的下降, 能量升高是不合理的。造成这种假象的原因可能是这2种方法的时间分辨率不够高导致出现在横波波至之后的高频伪瑞利波的能量被叠加到一起。重排谱图的时间边缘不仅能清晰地显示代表纵波的条带, 而且能将高频伪瑞利波与横波区分开。与干层相比, 油层的纵波、横波能量出现衰减, 但衰减程度没有水层大, 说明虽然油层的密度略小于水层, 地层颗粒相对滑动的热损耗大一些, 但流体对能量的吸收、声耦合率的改变对能量传递的阻碍程度没有水层高。虽然油层的剪切模量、体积模量、密度均小于干层, 但(K+4μ/3)在数值上减小的程度远小于密度, 而剪切模量和密度在数值上减小的程度相差无几, 使纵波速度减小, 横波速度无明显变化。因此, 油层的纵波波峰时间延后至1.07ms, 横波波峰时间无明显变化。

在油层IMF1的频率边缘图中, 短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布的显示为一条宽的条带, 且波峰能量不突出, 无法从中识别出纵波、横波。而在重排谱图的频率边缘图中可以看到2条明显的条带, 分别对应纵波的主频12.0kHz和横波的主频9.0kHz。与干层相比, 纵波主频无变化, 说明纵波各频率成分能量的相对大小未显著改变, 而横波主频有所减小则说明横波中的高频成分发生较大的衰减。

图 6为不同流体性质储集层声波全波列的IMF2分量边缘特征图。由前文分析可知IMF2分量主要包含的是伪瑞利波和高频斯通利波。干层的IMF2时间边缘中高频斯通利波能量最高, 波峰时间为2.89ms; 频率边缘图中Choi-Williams分布、平滑伪Wigner-Ville分布的能量较低, 短时傅里叶变换的能量较高, 而重排谱图能量最高, 显示十分明显, 高频斯通利波的主频为3.0kHz, 说明重排谱图通过将时频分析值重排使得信号能量得到聚焦, 与其它3种方法相比具有很大的优势。与干层相比, 水层的高频斯通利波的能量衰减较多, 但波峰时间延后不明显, 主频无明显变化。油水同层的高频斯通利波的能量衰减程度和水层差别不大, 波峰时间延后不明显, 主频无明显变化, 与水层的差别在于伪瑞利波条带的出现。油层的高频斯通利波的能量衰减十分严重, 图中的条带为伪瑞利波成分; 频率边缘图中也显示其能量衰减严重, 图中有显示的部分应为伪瑞利波成分。

图 7为不同流体性质储集层声波全波列的IMF3分量边缘特征图, 干层的IMF3时间边缘中低频斯通利波能量较高, 波峰时间为2.87ms; 频率边缘图中短时傅里叶变换、Choi-Williams分布、平滑伪Wigner-Ville分布的能量较低, 而重排谱图能量很高, 低频斯通利波的主频约2.0kHz, 再次展现了重排谱图的优势。对比干层、水层的低频斯通利波的能量衰减较多, 但波峰时间延后不明显, 主频无明显变化。油水同层的低频斯通利波的能量衰减程度和水层差别不大, 波峰时间延后不明显, 主频无明显变化。油层的低频斯通利波的能量衰减十分严重。总的来看, 油水同层的IMF2、IMF3分量时频边缘特征与水层区别不明显, 用IMF1分量的时频边缘特征区别水层、油水层的效果更好。

斯通利波是沿井壁传播的一种面波, 由射向井壁的声波与井壁界面接触时产生, 同时受到井壁两侧介质性质的影响。与干层地层相比, 水层和油层虽然改变了地层的剪切模量、体积模量、密度等参数, 但高频斯通利波、低频斯通利波的波峰时间均无明显变化, 说明地层参数的综合影响未造成高频斯通利波、低频斯通利波波峰时间的明显延后。斯通利波对流体十分敏感, 地层中有流体存在时, 斯通利波的能量会发生衰减。但高、低频斯通利波的主频均无明显变化, 且干层、水层及油层的高频斯通利波的能量均大于低频斯通利波的能量, 说明斯通利波各频率成分能量的相对大小未显著改变。不同流体性质储集层的声波全波列的边缘特征总结见表 2。

4 结束语

本文对比了4种基于EMD的时频分析方法处理阵列声波测井全波列的效果, 通过处理一个单点数据和提取时频分布的边缘特性进行井段数据处理比较了它们的时间分辨率、频率分辨率、抑制交叉项性能, 可得到以下结论。

1) 阵列声波测井信号通过EMD得到的IMF1分量主要包含纵波、横波信息; IMF2分量主要包含高频的斯通利波和伪瑞利波信息; IMF3分量主要包含低频的斯通利波信息。与直接对阵列声波测井信号进行时频分析处理相比, 对EMD得到的IMF分量做时频分析处理可获得不同分辨尺度的信息, 能更大程度地挖掘阵列声波测井信号中蕴含的信息。

2) 通过提取时频分布的边缘特性可以突破时频分布单点处理的局限, 达到处理井段数据的目的, 扩大了时频分布在阵列声波数据处理中的应用。

3) 重排谱图通过将时频分析值重排至时频域的重心, 在减小交叉项的同时改善时频聚焦性。在文中的4种方法中其时频分辨率最高、聚焦性最强, 能更好地刻画所处理信号中包含的频率、时间都比较接近的组分的特征。由时间边缘图可知, 该方法一方面使能量很低的纵波也能较清晰地显示, 另一方面对波峰时间的定位更精准, 有利于分析地层含油性不同时模式波能量峰的时间变化; 由频率边缘图也可看出, 重排谱图法使能量聚焦明显, 便于分析主频的变化。

4) 在4种方法中, 基于EMD的重排谱图区分干层、水层、油层的效果最好。以干层波列的纵波、横波、高频斯通利波、低频斯通利波对应的波峰时间、主频、能量为基准, 水层的纵波、横波能量衰减明显, 高频斯通利波、低频斯通利波能量均有一定的衰减; 油层的纵波、横波能量也出现衰减, 但衰减程度没有水层大, 高、低频斯通利波能量衰减则十分严重。时间边缘图显示水层的纵波、横波波峰时间有一定的滞后, 而油层纵波波峰时间滞后、横波波峰时间则无明显变化, 水层和油层的高、低频斯通利波波峰时间变化均不明显。频率边缘图显示水层纵波、横波主频无显著变化, 而油层的纵波主频无明显变化、横波主频降低, 水层、油层的高、低频斯通利波波峰主频变化均不明显。由此可见, 与干层相比, 水层和油层具有不同的时频边缘特征, 可以利用基于EMD的重排谱图方法达到区分水层、油层的目的。

本文中实际储层的岩性为杂色粗面岩, 但所用分析方法及效果可推广应用于其它岩性的储层。