0 引言

地震数据是一种典型的非平稳信号，在震源子波一定的假设条件下，不同时刻反射子波的振幅和频率特征是地下介质物理性质的综合体现^[1]。时频分析技术可方便地提取上述随时间变化的振幅、频率特征，因此该技术已发展成为常规地震解释工具。在地质异常体(如河道、断层等)识别^[2-6]、薄层刻画^[7-8]、油气异常检测^[9-11]、地震衰减参数求取^[12-14]和一些地震数据处理^{[11, 15]}中，各种时频分析方法都得到了广泛应用。

地震信号时频谱的计算方法种类繁多，其中应用最广泛的算法是加窗傅里叶变换。Gabor^[16]首先提出利用高斯滑动时窗截取信号，然后计算该截取信号的傅里叶谱，将其作为该窗口中心处的时频谱。Gabor变换窗函数固定，无法自适应地满足非平稳地震信号频率成分随时间变化的特性。因此，连续小波变换^[2]、S变换^[17]、广义S变换^[18]等窗函数随频率变化的时频分析方法被提出，所得时频谱的时间分辨率在高频端获得显著提升；然而在低频端由于窗口宽，所包含的信息一般来自多个相邻反射子波，造成时频谱上低频时间分辨率极低^[19]。低频信息(如“低频阴影”)通常对油气异常检测具有重要意义，因此上述常规方法尚无法满足当今地震勘探的“高分辨”要求。

研究表明，用反演方式实现加窗傅里叶变换可提高时频谱的时间分辨率。如Puryear等^[20]用约束最小二乘反演策略求解短时傅里叶时频谱，一定程度提高了时频分辨率。借助稀疏正则化策略，Gholami^[11]提出了稀疏短时傅里叶变换方法；在此基础上，Sattari等^[21]以所得时频谱具有最大稀疏度为原则，给出了时变自适应窗口的计算方法；更进一步，Sattari^[22]和Yuan等^[6]分别将稀疏反演时频分析的思想引入S变换和连续小波变换，显著提高了相应方法的时频分辨率。

尽管基于稀疏反演的方法显著提高了时频分辨率，但采用的稀疏约束策略对时间和频率两个方向的能量是无差别压缩，能量通常聚焦成点(即过于稀疏)，与子波真实傅里叶谱的形态相差甚远。此外，因地震子波具有空变特性，目的层聚焦频率在空间存在差异，故难以利用单频信息表征储层。

为了只在时间方向压缩时频能量，并在频率方向尽量保持子波傅里叶谱形态，本文采用组稀疏正则化反演策略，提出了新的Gabor时频分析方法。采用的高斯窗函数也会影响Gabor时频谱的质量^[21]，因此根据地震信号的质心频率给出了高斯窗函数自适应构造方法。合成数据实验表明：本文方法同时提高了Gabor时频谱在低频和高频端的时间分辨率；在实际应用中，低频剖面分辨率与高频剖面相当，为高分辨率地震资料解释提供了便利。

1 方法原理 1.1 基于反演的Gabor变换目标函数

假设u(t)表示一维地震信号，g(t)= $ \mathrm{e}^{\frac{-2 \pi^2 t^2}{(\Delta t \cdot s)^2}}$表示高斯窗函数(其中：Δt为信号采样间隔；s值越大窗口宽度越大，反之亦然)，则常规Gabor时频谱GT(f，τ)可表示为

$ \operatorname{GT}(f, \tau)=\int_{-\infty}^{+\infty} u(t) g(\tau-t) \mathrm{e}^{-\mathrm{i} 2 \pi f t} \mathrm{~d} t $

(1)

式中τ和f对应Gabor时频谱的时间和频率坐标。

对式(1)两端关于频率轴做傅里叶逆变换，可得二维信号

$ u(t) g(\tau-t)=\int_{-\infty}^{+\infty} \mathrm{GT}(f, \tau) \mathrm{e}^{\mathrm{i} \pi 2 f t} \mathrm{~d} f $

(2)

将式(2)中的二维时频谱GT(f，τ)重排为列向量，即表示为如下矩阵形式

$ \bar{\boldsymbol{W}} \boldsymbol{u}=\boldsymbol{F x} $

(3)

式中：列向量x为GT(f，τ)重排所得；列向量u由时域地震信号u(t)构成；F为分块对角阵，其中每个子阵都为N点傅里叶逆变换矩阵F_N^－1(N为地震信号采样点数)

$ \boldsymbol{F}=\left[\begin{array}{lll} \boldsymbol{F}_N^{-1} & & \\ & \ddots & \\ & & \boldsymbol{F}_N^{-1} \end{array}\right]_{N^2 \times N^2} $

(4)

W为列分块矩阵

$ \overline{\boldsymbol{W}}=\left[\begin{array}{c} \overline{\boldsymbol{W}}_1 \\ \vdots \\ \overline{\boldsymbol{W}}_i \\ \vdots \\ \overline{\boldsymbol{W}}_N \end{array}\right]_{N^2 \times N} $

(5)

其子阵$\overline{\boldsymbol{W}}_i$为高斯窗函数$g\left(\tau_i-t\right)$构成的对角阵。不难发现, 欲获得地震信号$\boldsymbol{u}$的表达式, 需在式(3) 等号两端同时乘以$\overline{\boldsymbol{W}}$的广义逆矩阵$\boldsymbol{W}$。由$\overline{\boldsymbol{W}}$的结构特点可知对应的$\boldsymbol{W}$为行分块矩阵

$ \boldsymbol{W}=\left[\begin{array}{llll} \boldsymbol{W}_1 & \cdots & \boldsymbol{W}_i & \cdots \boldsymbol{W}_N \end{array}\right]_{N \times N^2} $

(6)

式中子阵$\boldsymbol{W}_i$为函数$g\left(\tau_i-t\right) / h(t)$构成的对角矩阵。为使$\boldsymbol{W} \overline{\boldsymbol{W}}$为单位矩阵, 需使参数$h(t)=\sum\limits_{i=1}^N\left[g\left(\tau_i\right.\right.$ $-t)]^2 。$

式(3)两端同时乘以W，可得

$ \boldsymbol{u}=\boldsymbol{W F x} $

(7)

以式(7)为正演公式，在最小二乘意义下求解Gabor时频谱所用的目标函数可表示为

$ J(\boldsymbol{x})=\frac{1}{2}\|\boldsymbol{W Fx} -\boldsymbol{u}\|_2^2 $

(8)

式中‖●‖₂表示向量的L₂范数。

1.2 组稀疏正则化约束

组稀疏正则化由L_2，1范数实现，其表达式为

$ \begin{array}{l} R(\boldsymbol{x})=\mu\|\boldsymbol{x}\|_{2.1}=\mu \sum\limits_{i=1}^N\left\|\boldsymbol{x}^i\right\|_2 \\ \quad=\mu \sum\limits_{i=1}^N \sqrt{\left(x^{i , 1}\right)^2+\left(x^{i, 2}\right)^2+\cdots+\left(x^{i, N}\right)^2} \end{array} $

(9)

式中：μ为正则化参数；$\boldsymbol{x}^i=\left[x^{i, 1}, x^{i, 2}, \cdots, x^{i, N}\right]^{\mathrm{T}} $, 表示第i个采样点处频谱构成的列向量。

图 1展示了相应的分组策略，其中黑色方框内相同时刻的频谱被分配到同组。由式(9)可知，组与组之间的关系是绝对值之和，类比于L₁范数可促进组间稀疏性，组内各元素之间的关系是平方和开根号，类比于L₂范数可获得平滑解，因此获得的Gabor时频谱只在时间方向具有稀疏性，满足上述仅在时间方向压缩时频能量的要求。

式(8)与式(9)相加可得最终的反演目标函数

$ \begin{aligned} \varXi(\boldsymbol{x}) & =J(\boldsymbol{x})+R(\boldsymbol{x}) \\ & =\frac{1}{2}\|\boldsymbol{W F} \boldsymbol{x}-\boldsymbol{u}\|_2^2+\mu\|\boldsymbol{x}\|_{2, 1} \end{aligned} $

(10)

式中J(x)和R(x)分别表征数据拟合度和组稀疏度。

1.3 投影快速软阈值迭代(PFIST)求解

求解式(10)过程中，为了简化函数J(x)的求导运算，首先对其进行如下等价变形

$ \begin{aligned} \varXi(\hat{\boldsymbol{x}}) & =J(\hat{\boldsymbol{x}})+R(\hat{\boldsymbol{x}}) \\ & =\frac{1}{2}\|\boldsymbol{W} \hat{\boldsymbol{x}}-\boldsymbol{u}\|_2^2+\mu\left\|\boldsymbol{F}^{-1} \hat{\boldsymbol{x}}\right\|_{2, 1} \end{aligned} $

(11)

式中$ \hat{\boldsymbol{x}}=\boldsymbol{F} \boldsymbol{x}$。PFIST算法可对上式快速优化求解^[23]，第k步核心迭代方程为

$ \hat{\boldsymbol{x}}_{k+1 / 2}=\boldsymbol{F}\left\{\operatorname{prox}_{\alpha \mu}^{2, 1}\left[\boldsymbol{F}^{-1}\left(\hat{\boldsymbol{x}}_k-\alpha \boldsymbol{d}_k\right)\right]\right\} $

(12)

式中：d_k=$ \nabla J=\boldsymbol{W}^{\mathrm{T}}\left(\boldsymbol{W} \hat{\boldsymbol{x}}_k-\boldsymbol{u}\right)$为J($\hat{\boldsymbol{x}} $)的梯度(上标T为共轭转置)；α为梯度下降的最优步长，它与梯度和正演算子的关系^[24]为

$ \alpha=\frac{\boldsymbol{d}_k^{\mathrm{T}} \boldsymbol{d}_k}{\boldsymbol{d}_k^{\mathrm{T}} \boldsymbol{W}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{d}_k} $

(13)

prox_αμ^{2, 1}(●)为近端投影算子，对组内元素的作用为^[25]

$ \operatorname{prox}_{\alpha \mu }^{2, 1}\left(\boldsymbol{x}^i\right)=\boldsymbol{x}^i \cdot \max \left\{1-\frac{\alpha \mu}{\left\|\boldsymbol{x}^i\right\|_2}, 0\right\} $

(14)

式中“max”表示取最大值。在PFIST理论框架下，目标函数的详尽优化步骤如表 1所示。

1.4 高斯窗函数自适应构造方法

实际应用中，总是希望高斯窗函数截取的地震信号与地震反射子波尽可能接近。若将高斯窗函数视为滤波器，则其通带宽度H应与子波延续时间T₀相等。图 2为高斯窗函数通带宽度的确定方式，能量衰减3 dB(即最大振幅的0.7071倍)以内的部分视为通带，其余部分为阻带。令高斯函数g(t)等于0.7071时即可获得通带边界的坐标±v，进而得到通带宽度的计算公式

$ H=2 v=2 \Delta t s \sqrt{\frac{-\ln 0.7071}{2 \pi^2}} $

(15)

实际地震道中T₀难以精确计算，但可通过地震道的瞬时质心频率f_c对其进行估算

$ T_0=\frac{1}{f_{\mathrm{c}}} $

(16)

据上分析，令H=T₀，可得s与质心频率的关系

$ s=\frac{\frac{1}{2 \Delta t f_{\mathrm{c}}}}{\sqrt{\frac{-\ln 0.7071}{2 \pi^2}}} $

(17)

因此，构造自适应高斯窗函数的步骤为：①利用常规Gabor变换获得地震信号的时频谱；②计算τ_i(i=1, 2, ⋯，N)时刻频谱的质心频率f_c(τ_i)；③对质心频率曲线进行平滑处理；④利用式(17)计算各时刻的窗函数参数s(τ_i)，并代入高斯窗函数公式即获得各个时刻的窗函数。

2 方法验证 2.1 一维衰减信号试算

利用图 3a所示一维衰减地震信号验证自适应高斯函数构造方法。分析该一维信号的常规Gabor变换时频谱(图 3b)，为了使获得的质心频率曲线(图 3b黄色实线)尽量光滑，窗函数参数s可取较大值(150)，进一步平滑后的质心频率曲线为该图中的红色虚线。由式(17)求取的参数s如图 3c所示，参数s随时间增大而增大表明高斯窗函数的通带宽度也随时间增大，符合衰减地震信号中反射子波延续时间越来越长的客观规律。

图 3d中红色曲线为若干时刻处截取的地震信号(对应窗函数为黑色曲线)，对比图中蓝色曲线所示的原始信号可知，截取的地震信号与实际反射子波一致性较好，证明利用图 3c的s值构造的高斯窗函数通带宽度适用于该衰减地震信号。

图 4对比了不同方法获得的时频谱，其中图 4a、图 4b分别为采用自适应高斯窗函数和固定高斯窗函数(s=30)的Gabor变换结果。可见自适应Gabor变换(图 4a)在深层的高斯窗函数宽度大，时间分辨率低。尽管s=30的固定窗函数Gabor变换(图 4b)分辨率高，但由于窗函数通带宽度过小而破坏了反射子波，进而在低频端出现严重的“分叉”假象。图 4c中广义S变换时频谱在高频端时间分辨率较高，但低频端混叠严重，分辨率低。本文方法所得的时频谱如图 4d所示，其低频端与高频端的分辨率相当，且都明显优于其他方法。

2.2 正则化参数μ对时频谱的影响

正则化参数μ是决定本文方法时频谱效果的关键参数，图 5对比了采用不同μ时的计算结果。当μ取0时，式(11)退化为无约束最小二乘问题，此时获得的时频振幅谱(图 5a)与常规自适应Gabor变换结果(图 4a)完全一致。随着μ增大(图 5b和图 4d)，其时频谱的时间分辨率提高，但并不意味着μ值越大越好，过大的μ值会压制有效弱信号(如图 5c中0.6~1.1 s层段)。实际中一般对典型井旁地震道进行多次试验，选择既能获得高分辨率又能使时频能量团较好对应反射子波(特别是感兴趣的弱反射)的μ值。

2.3 薄层模型“低频阴影”检测

为了进一步说明本文方法在高分辨率时频分析方面的优势，设计了图 6a所示的“泥包砂”薄层模型。该砂岩储层时间厚度为20 ms，模型参数如表 2所示。在CDP30~CDP70之间砂岩含气，并引起低Q值异常，其他部分假设为弹性无衰减介质。以40 Hz雷克子波为源子波，采用非稳态褶积模型^[24]进行正演模拟，所得地震剖面如图 6b所示。由于含气砂岩对地震波的衰减和频散作用，可观察到砂岩底界反射波的延续时间明显变长，表明地震波低频分量占比增加。

分别用常规Gabor变换、广义S变换和本文方法计算图 6b中第50道的时频谱，如图 7所示。常规Gabor变换时频谱中含气砂岩顶、底界面反射子波对应的能量团混叠严重、时间分辨率低(图 7a)；广义S变换时频谱中高频端顶、底界面能量团可分辨，但低频端时间分辨率低(图 7b)；本文方法在低频端和高频端都可清晰分辨含气砂岩顶、底界面能量团，且底界面反射子波低频占比升高的现象更直观(图 7c)。图 8为30 Hz和85 Hz对应的三种方法得到的单频剖面，可见尽管常规Gabor变换(图 8a)和广义S变换(图 8b)能够获得“低频阴影”现象，但二者的低频剖面时间分辨率低，无法揭示薄储层的顶、底界面，而高频剖面揭示的顶、底界面由于能量混叠使时间位置发生移动。本文方法所得低频剖面和高频剖面分辨率大致相当，薄砂岩储层顶、底界面时间位置精确，进一步体现了本文方法的“高分辨”优势(图 8c)。

3 实际应用

图 9a为中国西北地区某深层砂岩油气储层地震剖面，井旁地震道的CDP号为14，粉色曲线为自然电位测井曲线，在时间为3.0 s处(红色箭头位置)钻遇工业油层。图 9b~图 9d分别展示利用自适应Gabor变换(窗函数参数见图 10b)、广义S变换和本文方法获得的井旁地震道时频谱，可见由于含油砂岩对地震波的强衰减作用，3.00~3.15 s的时频能量团(红色虚线框)向低频方向移动。对比三种方法可知，本文方法所得时频谱在低频端和高频端的时间分辨率最高(图 9d)。

图 10为图 9a的单频分析结果。在平滑质心频率剖面(图 10a)上，由于黏弹性衰减作用使时间大于3.0 s层段的质心频率出现低异常特征；s整体趋势依然符合地震子波延续时间随旅行时增大而变长的规律(图 10b)。分别应用自适应Gabor变换(图 10c、图 10d)、广义S变换(图 10e、图 10f)、本文方法(图 10g、图 10h)得到21 Hz和38 Hz单频剖面。可见三种方法的低频剖面上含油储层下部都出现强能量异常，但Gabor变换和广义S变换高频剖面上依然存在强能量，若无井数据对照，易对解释人员造成困扰；高频剖面上的强能量实际上缘于储层顶界面尚未经历衰减的反射子波，由于分辨率不足难以区分。而本文方法具有高分辨率优势，能够区分储层顶、底界面(图 10g和图 10h局部放大图中的红色箭头所示)，顶界面能量在低频和高频剖面上都很稳定，而底界面能量在高频剖面上显著降低，与薄层模型分析中的规律一致。

4 讨论

Gabor时频分析的本质是对信号进行分解，因此有必要从信号分解的角度考查本文方法与常规Gabor分析的差异。为得到各时间采样点处的时间域分解信号，将Gabor分解结果x代入下式

$ \boldsymbol{X}=\tilde{\boldsymbol{W}} \circ\operatorname{matrix}(\boldsymbol{Fx} , N, N) $

(18)

式中：matrix(●, N, N)表示将向量“●”重排为N行×N列的矩阵；N行×N列的矩阵$ \widetilde{\boldsymbol{W}} $中各列由g(τ_i－t)/h(t)构成；符号“$ \circ$”表示两个矩阵点乘。N行×N列的矩阵X即表示Gabor分析在各采样点处分解出的地震信号。图 11a和图 11b分别给出了常规Gabor分解和本文Gabor分解方法对图 3a的时间域分解结果及其叠加信号(红色虚线)，可见常规方法和本文方法的分解结果叠加后与原始信号(蓝色实线)都很吻合，但本文方法由于采用了组稀疏约束，相当于高斯窗函数仅作用于少数采样点后截取信号的合成结果就可拟合原始信号，因此分解结果在τ轴方向是稀疏的，对应时频谱上就会有更高的时间分辨率。常规Gabor分析在反射子波之间的采样点也会截取信号，所得分解结果的时间分辨率必然较低。

5 结论

从Gabor反变换角度出发，本文将Gabor时频谱的求解归结为反演问题，并在反演目标函数中加入组稀疏正则化约束项，实现了时频谱仅在时间方向的压缩，进而同步提高低频端和高频端的时间分辨率。通过推导高斯窗函数参数与反射子波质心频率之间的关系，避免了Gabor变换中高斯窗函数设置的盲目性。模型数据测试和实际地震资料处理结果表明，本文方法的时间分辨率明显优于常规Gabor变换和广义S变换。在实际地震数据的“低频阴影”分析中，本文方法所得低频和高频剖面能同时识别薄储层的顶、底界面。该方法的高分辨优势降低了解释不确定性，体现出较强应用潜力。