0 引言

地震信号的时—频表示（TFR）在地震数据处理和解释中得到广泛应用，如时—频域去噪^[1-2]、油气储层识别与储层描述^[3-4]等，因此，TFR方法研发一直受到业界的重视^[5-8]。基于傅里叶变换（FT）导出了一系列时频分析方法，如短时傅里叶变换（STFT）、连续小波变换（CWT）、S变换（ST）以及各种广义S变换（GST）^[9-15]。ST的时间窗是频率的函数，其时—频谱的时间分辨率依赖于频率^[11]。GST能够优化频率依赖性，但没有通用的优化标准分析窗函数，需要根据具体应用调节窗函数^{[4, 15]}。在上述TFR方法中，ST在有效性与计算效率之间具有好的平衡，已被广泛用于储层描述，但存在以下不足：①时—频谱的能量分布中心不在地震波主频处，而向高频方向偏移；②时—频谱的高频时间分辨率较高、低频时间分辨率较低；③在数值计算时奇异频率引起不稳定现象。上述不足导致ST时—频谱不能表示地震属性随深度的变化，而且较低的低频时间分辨率也不适合储层识别，因为低频振幅异常往往是油气储层的标志^[16-19]。由于地震波通过储层时高频能量被吸收^[20]，因此地震记录低频阴影区的能量分布可指示储层分布。

W变换（WT）^[21]是近年来提出的一种新时—频分析技术，它克服了ST的不足，其时—频谱的能量集中在地震波主频附近，而且提高了低频时间分辨率，适合于储层描述。Li等^[3]注意到WT的时—频谱在主频处出现振幅分裂现象，这是由WT的标准偏差函数在主频处的奇异性（不可微）所致；为此，设计了一种新的变频高斯标准偏差函数，具有主频的双参数多元复合指数形式，通过引入趋势因子p，获得了更光滑、更灵活的窗函数，得到广义WT（GWT）。由于新的标准偏差函数在主频处是可微的，消除了WT的时—频谱在主频处的振幅分裂现象，时—频谱的能量集中在主频附近，GWT提供了更好的TFR性能。随后，Chen等^[4]又在GWT中增加了一个标量因子，提出三参数WT（TPWT），可更灵活地调节时—频分辨率。TPWT是分析非平稳信号的有效工具，已成功用于油气储层识别。然而，前人没有详细论述TPWT的可逆性。为此，本文研究了TPWT的可逆性，首先从理论上证明TPWT不是严格可逆的，而是一种近似可逆的变换工具。合成地震数据与实际地震记录的数值计算展示了TPWT逆变换的重建精度，所做工作为TPWT的合理应用提供了参考。

1 TPWT数学理论 1.1 WT

地震记录$ x\left(t\right)\in {L}^{2}\left(R\right) $(平方可积函数空间)的WT定义为^[21]

$ W(\tau ,f)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x\left(t\right){g}_{\sigma }\left(t-\tau ,f;\tau \right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}2\mathrm{\pi }ft\right)\mathrm{d}t $

(1)

式中：$ \tau \mathrm{为} $时移因子；$ f\mathrm{为} $频率；t为时间；$ {g}_{\sigma }\left(t,f;\tau \right)\mathrm{为} $高斯窗函数，即

$ {g}_{\sigma }(t,f;\tau )=\frac{1}{\sqrt[]{2\mathrm{\pi }}\sigma (\tau ,f;k)}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{t}^{2}}{2{\sigma }^{2}(\tau ,f;k)}\right] $

(2)

式中$ \sigma \left(\tau ,f;k\right) $为标准偏差函数，定义为

$ \sigma (\tau ,f;k)=\frac{k}{{f}_{0}\left(\tau \right)+\left|{f}_{0}\left(\tau \right)-f\right|} $

(3)

式中：k为标量因子；$ {f}_{0}\left(\tau \right) $为地震波主频。显然，有

$ \sigma (\tau ,f;k)=\left\{\begin{array}{c}\begin{array}{cc}\frac{k}{f}& f > {f}_{0}\end{array}\\ \begin{array}{cc}\frac{k}{2{f}_{0}-f}& \begin{array}{c}\end{array}f\le {f}_{0}\end{array}\end{array}\right. $

(4)

$ \sigma (\tau ,f;k) $对$ f $的偏导数为

$ \frac{\partial \sigma(\tau, f ; k)}{\partial f}=\left\{\begin{array}{cl} -\frac{\sigma^2(\tau, f ; k)}{k} & f>f_0 \\ \text { 不存在 } & f=f_0 \\ \frac{\sigma^2(\tau, f ; k)}{k} & f<f_0 \end{array}\right. $

(5)

可见$ \sigma (\tau ,f;k) $在$ {f}_{0} $处不可微。

1.2 TPWT

把$ {g}_{\sigma }\left(t,f;\tau \right) $改为多变量组合椭圆函数形式，即

$ g_\sigma(t, f ; \tau)=\frac{1}{\sqrt{2 \pi} \sigma(\tau, f ; \boldsymbol{\varPhi})} \exp \left[-\frac{t^2}{2 \sigma^2(\tau, f ; \boldsymbol{\varPhi})}\right] $

(6)

则标准偏差函数定义为

$ \sigma(\tau, f ; \boldsymbol{\varPhi})=\frac{1}{\sqrt[p]{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p}} $

(7)

式中$ \boldsymbol{\varPhi}=\left\{{k}_{1},{k}_{2},p\right\} $是一组调节参数，标量因子$ {k}_{1} $调节三参数高斯窗的宽度，标量因子$ {k}_{2} $控制$ f $和$ {f}_{0} $之间高斯窗尺度的份额比例，p为高斯窗的趋势因子，用于调节标准偏差的变化率。

结合式(6)和式(7)，得到TPWT^[4]

$ \begin{array}{l}\mathrm{T}\mathrm{P}\mathrm{W}\mathrm{T}(\tau ,f)\mathrm{ }=\frac{\sqrt[p]{{\left[\frac{{f}_{0}\left(\tau \right)}{{k}_{1}}\right]}^{p}+{\left[\frac{\left|{f}_{0}\left(\tau \right)\mathrm{ }-f\right|}{{k}_{2}}\right]}^{p}}}{\sqrt[]{2\mathrm{\pi }}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x\left(t\right)\times \\ \\ \mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{(t-\tau )}^{2}}{2}{\left\{\sqrt[p]{{\left[\frac{{f}_{0}\left(\tau \right)}{{k}_{1}}\right]}^{p}+{\left[\frac{\left|{f}_{0}\left(\tau \right)\mathrm{ }-f\right|}{{k}_{2}}\right]}^{p}}\right\}}^{2}\right)\times \\ \mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}2\mathrm{\pi }f\right)\mathrm{d}t\end{array} $

(8)

不同参数的TPWT的表现形式不同：

（1）当$ {k}_{1}={k}_{2}=1 $且$ f={f}_{0}\left(\tau \right) $或$ {f}_{0}\left(\tau \right)=0 $时，TPWT与ST等价；

（2）当$ {k}_{2}\ne 1 $且$ {f}_{0}\left(\tau \right)=0 $时，TPWT与GST等价；

（3）当$ {k}_{1}={k}_{2} $且$ p=1 $时，TPWT退化为WT；

（4）当$ p\to \mathrm{\infty } $时，TPWT退化为STFT。

式(7)可改写为

$ \sigma(\tau, f ; \boldsymbol{\varPhi})= \begin{cases}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{f-f_0(\tau)}{k_2}\right]^p\right\}^{-\frac{1}{p}} & f>f_0 \\ f_0(\tau) & f=f_0 \\ \left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{f_0(\tau)-f}{k_2}\right]^p\right\}^{-\frac{1}{p}} & f<f_0\end{cases} $

(9)

$ \sigma (\tau ,f;\boldsymbol{\varPhi}) $对$ f $的偏导数为

$ \begin{array}{l}\frac{\partial \sigma (\tau ,f;k)}{\partial f}=\\ \begin{cases}-\frac{\left[f-f_0(\tau)\right]^{p-1}}{k_2^p\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{f-f_0(\tau)}{k_2}\right]^p\right\}^{\frac{1+p}{p}}} & f>f_0 \\ 0 & f=f_0 \\ \frac{\left[f_0(\tau)-f\right]^{p-1}}{k_2^p\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{f_0(\tau)-f}{k_2}\right]^p\right\}^{\frac{1+p}{p}}} & f<f_0 \\ \end{cases}\end{array} $

(10)

可见$ \sigma (\tau ,f;k) $在$ {f}_{0} $处可微。

TPWT引入趋势因子$ p $克服WT的振幅分裂现象，通过$ {k}_{1} $和k₂可更灵活地调节分辨率。

2 TPWT的可逆性

在理论上，TPWT的逆变换具有如下形式

$ x\left(t\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{T}\mathrm{P}\mathrm{W}\mathrm{T}(\tau ,f)\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}2\mathrm{\pi }ft\right)\mathrm{d}\tau \mathrm{d}f $

(11)

可分解为

$ X\left(f\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{T}\mathrm{P}\mathrm{W}\mathrm{T}(\tau ,f)\mathrm{d}\tau $

(12)

$ x\left(t\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}X\left(f\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}2\mathrm{\pi }ft\right)\mathrm{d}f $

(13)

式(12)为无损可逆条件^{[13, 21]}，式(13)为傅里叶逆变换（IFT）。当且仅当式(12)和式(13)同时成立时，才能保证式(11)成立。由式（8）可知

$ \begin{array}{l}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{T}\mathrm{P}\mathrm{W}\mathrm{T}\left(\tau ,f\right)\mathrm{d}\tau ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\sqrt[p]{{\left[\frac{{f}_{0}\left(\tau \right)}{{k}_{1}}\right]}^{p}+{\left[\frac{\left|{f}_{0}\left(\tau \right)\mathrm{ }-f\right|}{{k}_{2}}\right]}^{p}}}{\sqrt[]{2\mathrm{\pi }}}\times \\ {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x\left(t\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{(t-\tau )}^{2}}{2}{\left\{\sqrt[p]{{\left[\frac{{f}_{0}\left(\tau \right)}{{k}_{1}}\right]}^{p}+{\left[\frac{\left|{f}_{0}\left(\tau \right)\mathrm{ }-f\right|}{{k}_{2}}\right]}^{p}}\right\}}^{2}\right)\times \\ \mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}2\mathrm{\pi }ft\right)\mathrm{d}t\mathrm{d}\tau \\ ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left\{{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{1}{\sqrt[]{2\mathrm{\pi }}}{\left\{{\left[\frac{{f}_{0}\left(\tau \right)}{{k}_{1}}\right]}^{p}+{\left[\frac{\left|{f}_{0}\left(\tau \right)\mathrm{ }-f\right|}{{k}_{2}}\right]}^{p}\right\}}^{\frac{1}{p}}\right.\times \\ \left.\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{(\tau -t)}^{2}}{2}{\left\{{\left[\frac{{f}_{0}\left(\tau \right)}{{k}_{1}}\right]}^{p}+{\left[\frac{\left|{f}_{0}\left(\tau \right)-f\right|}{{k}_{2}}\right]}^{p}\right\}}^{\frac{2}{p}}\right)\mathrm{d}\tau \right\}\times \\ x\left(t\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}2\mathrm{\pi }ft\right)\mathrm{d}t\end{array} $

(14)

对比式(14)与式(12)可知

$ \begin{aligned} & \int_{-\infty}^{+\infty} \frac{1}{\sqrt{2 \pi}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p\right\}^{\frac{1}{p}} \times \\ & \exp \left(-\frac{(\tau-t)^2}{2}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p\right\}^{\frac{2}{p}} \mathrm{~d} \tau=1\right. \end{aligned} $

(15)

令$ u=\frac{1}{\sqrt[]{2}}{\left\{{\left[\frac{{f}_{0}\left(\tau \right)}{{k}_{1}}\right]}^{p}+{\left[\frac{\left|{f}_{0}\left(\tau \right)-f\right|}{{k}_{2}}\right]}^{p}\right\}}^{\frac{1}{p}}(\tau -t) $，则

$ \begin{aligned} & \mathrm{d} u=\frac{1}{\sqrt{2}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p\right\}^{\frac{1}{p}} \mathrm{~d} \tau+ \\ & (\tau-t) \mathrm{d}\left(\frac{1}{\sqrt{2}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p\right\}^{\frac{1}{p}}\right) \\ & =\left\{\begin{array}{l} \frac{1}{\sqrt{2}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p\right\}^{\frac{1}{p}} \mathrm{~d} \tau+ \\ (\tau-t) \frac{1}{\sqrt{2}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p\right\}^{\frac{1}{p}} \times \\ \left\{\frac{\left[f_0(\tau)\right]^{p-1}}{k_1^p}+\frac{\left[f_0(\tau)-f\right]^{p-1}}{k_2^p}\right\} f_0^{\prime}(\tau) \mathrm{d} \tau \quad f<f_0 \\ \frac{1}{\sqrt{2}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{f-f_0(\tau)}{k_2}\right]^p\right\}^{\frac{1}{p}} \mathrm{~d} \tau+ \\ (\tau-t) \frac{1}{\sqrt{2}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{f-f_0(\tau)}{k_2}\right]^p\right\}^{\frac{1}{p}} \times \\ \left\{\left[\frac{\left[f_0(\tau)\right]^{p-1}}{k_1^p}+\frac{\left[f-f_0(\tau)\right]^{p-1}}{k_2^p}\right\} f_0^{\prime}(\tau) \mathrm{d} \tau \quad f>f_0\right. \end{array}\right. \end{aligned}$

(16)

只有当

$ \mathrm{d}u\approx \frac{1}{\sqrt[]{2}}{\left\{{\left[\frac{{f}_{0}\left(\tau \right)}{{k}_{1}}\right]}^{p}+{\left[\frac{\left|{f}_{0}\left(\tau \right)-f\right|}{{k}_{2}}\right]}^{p}\right\}}^{\frac{1}{p}}\mathrm{d}\tau $

(17)

时，才能得到

$ \begin{aligned} & \int_{-\infty}^{+\infty} \frac{1}{\sqrt{2 \pi}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p\right\}^{\frac{1}{\rho}} \times \\ & \exp \left(-\frac{(\tau-t)^2}{2}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{\left|f_0(\tau)-f\right|}{k_2}\right]^p\right\}^{\frac{2}{p}} \mathrm{~d} \tau\right. \\ & \approx \frac{1}{\sqrt{\pi}} \int_{-\infty}^{+\infty} \exp \left(-u^2\right) \mathrm{d} u=1 \end{aligned} $

(18)

此时，由式(14)可得

$ \begin{aligned} & \int_{-\infty}^{+\infty} \operatorname{TPWT}(\tau, f) \mathrm{d} \tau \approx \int_{-\infty}^{+\infty} x(t) \exp (-\mathrm{i} 2 \pi f t) \mathrm{d} t \\ & \quad=X(f) \end{aligned} $

(19)

因此

$ {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{T}\mathrm{P}\mathrm{W}\mathrm{T}(\tau ,f)\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}2\mathrm{\pi }ft\right)\mathrm{d}\tau \mathrm{d}f\approx x\left(t\right) $

(20)

式（20）在理论上证明了TPWT并不是严格可逆的，而是一种近似可逆的变换工具，这是由于在定义高斯窗函数时使用了非平稳的频率权重所致。

比较式(16)和式(17)可知，在推导近似逆变换公式（式(20)）的过程中引入了近似表达式（式(17)），而忽略了$ \mathrm{d}u $的以下项

$ \begin{gathered} (\tau-t) \frac{1}{\sqrt{2}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{f_0(\tau)-f}{k_2}\right]^p\right\}^{\frac{1}{p}} \times \\ \left\{\frac{\left[f_0(\tau)\right]^{p-1}}{k_1^p}+\frac{\left[f_0(\tau)-f\right]^{p-1}}{k_2^p}\right\} f_0^{\prime}(\tau) \mathrm{d} \tau \quad f<f_0 \end{gathered} $

(21)

或

$ \begin{gathered} (\tau-t) \frac{1}{\sqrt{2}}\left\{\left[\frac{f_0(\tau)}{k_1}\right]^p+\left[\frac{f-f_0(\tau)}{k_2}\right]^\rho\right\}^{\frac{1}{p}} \times \\ \left\{\frac{\left[f_0(\tau)\right]^{p-1}}{k_1^\rho}+\frac{\left[f-f_0(\tau)\right]^{p-1}}{k_2^\rho}\right\} f_0^{\prime}(\tau) \mathrm{d} \tau \quad f>f_0 \end{gathered} $

(22)

式(21)和式(22)均与$ {f}_{0}\left(\tau \right) $、$ {f}_{0}^{\text{'}}\left(\tau \right) $以及$ \left\{{k}_{1},{k}_{2},p\right\} $有关。所以，逆TPWT的重建误差与原始输入数据和变换参数$ \left\{{k}_{1},{k}_{2},p\right\} $有关。下文通过数值计算显示不同输入数据和变换参数对重建误差的影响。

3 数值计算

图 1为不同方法计算的主频为60 Hz的Morlet小波$ \psi \left(t\right) $（图 1a）的时—频谱。由图可见，ST的时—频谱的能量分布中心向高频偏移（图 1b），WT的时—频谱在主频60 Hz处存在明显的奇异性，振幅发生分裂（图 1c~图 1e），TPWT的时—频谱的能量分布中心在主频处（图 1f）。图 1e在0.4 s处的归一化振幅（图 2）更清楚地展示了振幅分裂现象。图 3为图 1b、图 1d、图 1f在0.4 s处的归一化振幅与图 1d、图 1f的标准偏差函数σ随频率的变化。由图可见：①ST的时—频谱的能量分布中心向高频偏移（图 3a）；TPWT的时—频谱的能量分布中心在主频处（图 3a），而且消除了振幅分裂现象。②由于WT的σ在主频处不可微，所以其时—频谱振幅在主频处存在分裂现象，TPWT的σ在主频处可微，所以其时—频谱振幅在主频处不存在分裂现象（图 3b）。

为了检验TPWT逆变换的近似程度，分别采用合成地震记录和实际地震数据进行分解与重构计算。首先，合成了一个地震记录$ x\left(t\right) $(图 4a)，由9个雷克子波组成，主频分别为70、56、60、52、40、48、30、15、10 Hz，用复数道分析方法^[24-25]计算$ x\left(t\right) $的主频(图 4b), 然后分别得到ST(图 4c)、WT(图 4d)、TPWT（图 4e)的时—频谱。图 5为不同逆变换重建的振幅谱与标准傅里叶谱|X(f)|的对比。可见：由逆ST重建的振幅谱只在低频端与|X(f)|有微小偏差(图 5a、图 5d），这是由数值计算噪声引起^[22-23]；由逆WT(图 5b、图 5d）、逆TPWT(图 5c、图 5d)重建的振幅谱则与|X(f)|的偏差较大，这是由式(18)的计算误差所致。图 6为不同逆变换重建的地震记录与图 4a的对比。可见：由于ST在理论上是严格可逆的，数据重建误差很小(图 6a、图 6d)，且由数值计算噪声引起；由逆TPWT重建的地震记录则与图 4a差别较大(图 6c、图 6d），主要由理论近似造成，如图 6d的最大振幅误差绝对值约为0.4，即相对误差约为20%。为了观察变换参数$ \boldsymbol{\varPhi }=\left\{{k}_{1},{k}_{2},p\right\} $对重建数据的影响，用三组不同的参数分别计算(图 7)，得到最大振幅误差绝对值分别为0.4270、0.2320、0.2883(图 7d），即相对误差分别为21.35%、11.6%、14.42%(图 7e)。

图 8为实际资料叠加剖面，随机抽取第95道（图 9a），利用复数道分析方法计算瞬时频率，再用一个低通滤波器通过加权平均消除数值噪声^[21]，从而获得主频（图 9b），最终得到不同变换参数的TPWT时—频谱（图 9c~图 9e）。可见，重建数据的最大振幅误差的绝对值分别为401.45、426.92、381.55，即相对误差分别为12.07%、12.83%、11.47%（图 9f、图 9g）。在上述数值计算中，相对误差最小值为11.47%，最大值为21.35%，即误差是明显的。因此，TPWT的不可逆导致其应用受到局限，在去噪、高分辨率处理等需要重建数据的处理领域不宜使用。

4 结论

（1）TPWT既解决了ST的低频时间分辨率低与时—频谱能量分布中心向高频偏移的问题，也克服了WT时—频谱在主频处的振幅分裂现象，能更准确地描述油气储层，更有利于地震解释。

（2）在理论上TPWT不是严格可逆的，而是一种近似可逆的变换工具，这与FT、ST不同。

（3）合成地震数据与实际地震记录的数值计算结果显示，与原始地震数据相比，利用逆TPWT重建地震道的相对误差为11.47%~21.35%，即理论上的不可逆导致较大的重建误差，严重影响TPWT的应用范围，在去噪、高分辨率处理等需要重建数据的处理领域不宜使用。