0 引言

高分辨率地震资料是研究陆相薄层、薄互层储层的基础^[1-2]。时频分析技术现已成为研究薄层的一种重要手段，作为非平稳信号分析的重要工具之一^[3]，已经广泛应用于储层预测和地震资料处理的各个方面。在储层预测方面，李思源等^[4]有效结合Wigner-Ville分布与Chrip-Z变换识别地下微型河道、薄储层等复杂油气储集体；尚平萍等^[5]结合自适应噪声的总体集合经验模态分解和希尔伯特变换进行油气检测。在地震资料处理方面，郭爱华等^[6]通过Shearlet域系数补偿提高叠后地震资料分辨率；江雨濛等^[7]利用同步挤压S变换的时变子波提取方法提高地震资料处理的精度；杨子鹏等^[8]结合广义S变换理论和压缩感知提高地震资料分辨率；纪永祯等^[9]通过贝叶斯学习方法和Wigner-Ville分布提高地震资料的信噪比、分辨率。

连续小波变换(Continuous Wavelet Transform，CWT)通过尺度因子协调时间和频率分辨率，已广泛应用于地震信号的处理和分析^[10]。可选Mexican hat^[11]、Daubechies^[12]、Gabor等小波作为CWT的基函数，其中Gabor小波在表征时频局部特性方面具有如下特性：①在局部频率分析上，Gabor小波具有较高的时频分辨率^[13]；②Gabor小波参数多，具有很强的灵活性和可操作性^[14-15]。然而，Gabor小波的可调因子控制参数是根据经验预设为常数，导致CWT的时频特征的分辨率具有一定的局限性。为提高CWT时频的分辨率，准确识别地震信息，本文介绍了一种基于可调因子Gabor小波(Tunable Factor Gabor Wavelet，TFGW)的连续小波变换(TFGW-CWT)。TFGW-CWT首先采用可调因子的Gabor小波对地震数据进行时频分解，然后基于最小绝对值投影方法优选多个连续小波系数中的极小绝对值，保留小波系数中的最小值，最终形成连续小波系数优化组合的高分辨率时频数据体。在此基础上，通过非负矩阵分解(Non-negative Matrix Factorization，NMF)对TFGW-CWT方法得到的时频数据体进行降维，得到高分辨率的地震数据。模型正演记录和实际数据处理结果证实该方法在薄层识别中具有一定的应用价值。

1 方法原理

本文基于TFGW-CWT进行时频分析，进一步结合NMF对上述时频数据进行降维，得到高分辨率数据，该数据与原始数据具有相似的空间结构。

1.1 TFGW-CWT

L²(R)空间中的任意函数f(t)可用小波基展开，这种展开称为函数f(t)的连续小波变换^[16-17]，其表达式为

$ C(a, b)=\int_{-\infty}^{\infty} f(t) \frac{1}{\sqrt{a}} \bar{\psi}\left(\frac{t-b}{a}\right) \mathrm{d} t $

(1)

式中：a为尺度因子；b为平移因子；ψ为小波基函数(又称母小波)；ψ为ψ的共轭。小波基函数需要满足可容性条件，即

$ {\mathit{\Phi }_\psi } = 2\pi \int_{ - \infty }^\infty {\frac{{|\widetilde \psi (\omega ){|^2}}}{\omega }} {\rm{d}}\omega < \infty $

(2)

式中${\widetilde \psi }$为ψ的频域表示。

可供CWT选择的小波基有很多。为了获取高分辨率的时频变换结果，可改变Gabor小波基的高斯标准差，并且保持复指数频率不变，达到提升Gabor小波基的局部频率分析性能^[18]。改进后的Gabor小波基为

$ \psi_{\mathrm{GW}}(t, b, \sigma)=\exp \left[\frac{-(t-b)^2}{2 \sigma^2}\right] \exp \left[{\rm{i}} \omega_0(t-b)\right] $

(3)

式中：σ为控制高斯包络线的标准差；ω₀为控制中心频率的参数。当σ=1、ω₀≥5，且满足可容许条件时^[19-20]，Gabor小波基即为Morlet小波基，有最佳的时频分辨率，其表达式为

$ \psi_{\mathrm{MW}}(t)=\exp \left(-\frac{t^2}{2}\right) \exp \left({\rm{i}} \omega_0 t\right) $

(4)

为了进一步提高CWT结果的分辨率，用可调因子T代替式(3)中的σ，得到带T的Gabor小波(TFGW)基

$ \psi_{\mathrm{TFGW}}(t)=\exp \left\{-\left[\frac{\omega_0(t-b)}{2 T}\right]^2\right\} \exp \left[\mathrm{i} \omega_0(t-b)\right] $

(5)

T的表达式^[21]为

$ T=\frac{\omega_c}{B_W}=\frac{\omega_0}{\sqrt{2} / \sigma} $

(6)

式中：ω_c为中心频率；B_W为频带宽度。

结合式(1)和式(5)即可得可调因子Gabor小波变换(TFGW-CWT)

$ C_{\mathrm{TFGW}-\mathrm{CWT}}(a, b)=\int_{-\infty}^{\infty} f(t) \frac{1}{\sqrt{a}} \bar{\psi}_{\mathrm{TFGW}}\left(\frac{t-b}{a}\right) \mathrm{d} t $

(7)

需要指出的是，T一般可通过实验取得。进一步利用最小绝对值投影方法^[22]将不同T得到的时频谱进行组合，即每一点都保留绝对值最小的系数作为TFGW-CWT的结果。

为了证实该方法具有高分辨率的时频特性，采用5、10、15和25 Hz的正弦信号合成一个信号

$ \begin{aligned} f= & 10 \sin (10 \pi t)+8 \sin (20 \pi t)+ \\ & 5 \sin (30 \pi t)+\sin (50 \pi t) \end{aligned} $

(8)

进行分析。由TFGW-CWT方法的时频分析结果(图 1) 可以看出：当$T=3 \sqrt{2}$时，合成信号在频率10 Hz与15 Hz之间存在能量交叉干扰(图 1a黑色椭圆处)，且25 Hz所对应的能量较低，特征不明显(图 1a红色椭圆处)；随着T的不断增大，合成信号的时频谱逐渐变窄，10 Hz与15 Hz之间能量交叉干扰也逐渐减弱，其中25 Hz处的特征更清晰(图 1b、图 1c红色椭圆处)。通过T的两个特定值($10 \sqrt{2}$和$50 \sqrt{2}$)得到两个不同的时频谱，然后利用最小绝对值投影方法进行比较，得到一个绝对值最小的时频谱(图 1d)。由图 1d可见，TFGW-CWT计算的时频谱更聚焦、特征更明显，分辨率更高。经典的连续小波变换和广义S变换的时频谱(图 1e和图 1f)在10 Hz和25 Hz处能量较弱(黑色和红色椭圆处)。

1.2 NMF算法原理

NMF是一个最优化问题，可以用迭代方法求解。NMF算法定义为：对于任意给定的一个非负矩阵S，其能够寻找到两个非负矩阵W和H，满足S=WH，从而将一个非负的矩阵分解为两个非负矩阵的乘积^[23]。其中，S是待分解矩阵；W为基矩阵；H为系数矩阵或权重矩阵。将式(7)得到的小波系数C_TFGW-CWT作为待分解矩阵S。据NMF算法求取一个系数矩阵H，并用H代替S，即可实现小波系数的降维，从而得到高精度的数据结构体。

为求取最优H结果，首先必须定义一个目标函数。针对NMF问题的目标函数有多种，可以利用S和WH的欧氏距离^[24]衡量，即

$ E=\|\boldsymbol{S}-\boldsymbol{W H}\|_2^2 $

(9)

求解上述优化问题迭代公式为

$ H_{d j} \leftarrow H_{d j} \frac{\left(\boldsymbol{W}^{\mathrm{T}} \boldsymbol{S}\right)_{d j}}{\left(\boldsymbol{W}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{H}\right)_{d j}} $

(10)

$ W_{i d} \leftarrow W_{i d} \frac{\left(\boldsymbol{S H}^{\mathrm{T}}\right)_{i d}}{\left(\boldsymbol{W H} \boldsymbol{H}^{\mathrm{T}}\right)_{i d}} $

(11)

1.3 基于TFGW-CWT和NMF的高分辨率处理方法

以TFGW-CWT和NMF两种算法为基础，充分利用上述方法在局部频率、空间结构信息的基础上形成一种基于NMF挖掘高分辨率时频特征的数据处理方法，其具体实现过程(图 2)为：首先，基于TFGW-CWT对地震数据进行时频分析，获得T取不同值时的连续小波系数；使用最小绝对值投影方法组合不同T的连续小波系数，作为时频分析结果；进一步利用NMF算法不断迭代优化(式(9)~式(11))对S进行分解，得到W(基矩阵)和H(特征数据结构关系体)；最后，用H代替TFGW-CWT变换的频率数据体，即可得到高精度的数据结构体。

2 理论模型试验

为了验证TFGW-CWT和NMF算法结合的有效性，设计如图 3所示的时间域楔状体模型。采用主频为35 Hz的Ricker子波与反射系数褶积方法得到模型合成记录(图 4a)。

用本文方法对合成记录进行处理，结果如图 4b所示。可以看出，本文方法提高了正演记录的分辨率，尖灭点(椭圆所示)更接近真实位置，能更精确识别楔状体。经本文方法处理后，带宽从10~70 Hz拓展到10~120 Hz，高频部分信息增强，低频部分基本保持与原始合成记录一致(图 5)。为了验证本文方法的抗噪性，对正演记录加入信噪比为10 dB的高斯白噪声(图 4c)，处理结果如图 4d所示。对比图 4d与图 4c和图 4b可以看出，本文方法不仅提高了分辨率，还具有一定的抗噪性。

3 实际应用

采用X工区三维数据验证本文方法的实用性。图 6a是原始地震剖面，仅AW1井可用于标定。图 6b是本文方法处理后的剖面。从地震层位解释的角度看，在图 6a与图 6b上均可以拾取T1、T11、T2等层位，表明处理前、后大的地震层位解释比较一致，地震剖面的构造特征相似。进一步，从T1和T11层间的反射特征看，处理后的剖面层间细节更丰富，且与AW1井测井曲线更吻合(图 6圆圈所示)。经本文方法处理后，分辨率更高，实际地震数据的频带从7~80 Hz拓展到7~115 Hz(图 7)。

4 结束语

本文基于可调因子的Gabor小波变换，利用最小绝对值投影方法降低连续CWT中邻近频率的交叉干扰，得到一个高分辨率的时频数据体。进一步结合NMF算法挖掘时频结果的关键特征，获得聚焦性强的关键特征频率数据，最终达到提高地震资料分辨率之目的。本文方法为提高地震资料分辨率处理提供了一个新的途径，在薄层识别中有一定的应用前景。

本文方法虽然可以达到提高地震资料分辨率的目的，但难免会产生假频和假的同相轴。借助钻、测井及地质资料构建基于模型和数据驱动的处理方法、降低多解性是今后的研究重点。