广义线性频域Chirplet变换及其应用

Generalized Linear Frequency-Domain Chirplet Transform and Its Application

摘要: 时频分析可以获取地下储层的相关信息，在储层表征中发挥着重要作用。本文提出了广义线性频域Chirplet变换(General linear Frequency-domain chirplet transform, GLFCT)，首先旋转时频基来匹配实际群延迟(group delay, GD)曲线，然后根据每一个时频点的振幅，得到最佳chirp率，最终得到高分辨率的时频表征结果。合成信号的数值验证结果证明了GLFCT处理瞬态信号的有效性，模拟地震记录的测试很好地证实了在储层表征方面的可行性。

Abstract: Time-frequency analysis can obtain the relevant information of underground reservoir, which plays an important role in reservoir characterization. In this paper, the General linear Frequency-domain chirplet transform (GLFCT) is proposed. Firstly, the time-frequency basis is rotated to match the actual group delay (GD) curve. Then, according to the amplitude of each time-frequency point, the optimal chirp rate is obtained, and finally the high-resolution time-frequency characterization results are obtained. The numerical verification results of the synthesized signal demonstrate the effectiveness of GLFCT in processing transient signals, and the test of simulated seismic records well proves the feasibility in reservoir characterization.

1. 引言

时频分析方法能提供时间域和频率域的联合分布信息，是处理地震信号等非平稳信号的常用工具之一。时频分析方法的结果可以刻画地震信号的局部特征，从而反映这些局部特征对应的地质结构和储层信息。传统的时频分析方法主要分为线性时频分析方法和二次型时频分析方法。常用的线性时频分析方法有短时傅里叶变换(Short-time Fourier Transform, STFT) [1] 、连续小波变换(Continuous Wavelet Transform, CWT) [2] 和S变换(Stockwell Transform, ST) [3] ，现已广泛应用于地震资料去噪、储层预测和烃类检测等方面。STFT用移动的窗函数截取信号，并对截取的信号进行傅里叶变换，最终得到了信号的时间–频率谱。Partyka [4] 等于1999年将STFT用于墨西哥湾某工区三维地震数据体，以确定薄层的位置、厚度及刻画地震数据体的不连续性。但是STFT的窗函数是固定的，所以其时频分辨率也是固定的。地震信号作为典型的非平稳信号，其频率成分随着时间不断变化，需要进行多分辨率分析。因此，Morletc [2] 基于可伸缩的小波基提出了CWT，使其具有多分辨率特性，即在低频处时间分辨率低、频率分辨率高，高频处时间分辨率高、频率分辨率低。高静怀等 [5] 于1996年研究了地震资料处理中小波函数的选取问题，提出了一种匹配地震子波的基本小波函数，并用于地震资料去噪及谱成像。Stockwell [3] 将STFT和CWT相结合提出ST，其窗函数随频率变化。由于S变换高斯窗函数的宽度随频率增大而变窄，因而发展出一系列改进的广义S变换 [6] ，在河道形态刻画 [7] 、薄储层厚度检测 [8] 、非均质性储层识别 [9] 等领域取得良好的应用效果。总之，上述方法虽然实现容易，但受Heisenberg测不准原理的约束，其时间分辨率和频率分辨率无法同时得到提高。

WVD (Wigner-Ville distribution, WVD) [10] [11] 作为一种典型的二次型时频分析方法，不受海森伯不确定原理的限制，并具有实值性、能量守恒性和时移频移不变性等数学性质。该类方法被广泛应用在河道形态刻画 [12] 、薄储层预测 [13] 等领域。尽管该类方法在致密河道砂岩形态刻画及其含气性识别中取得了一定的效果，但在处理多分量信号时无法避免交叉项的干扰，极大地限制了其应用 [14] 。为此，提出了一种参数化时频分析方法来提高时频分辨率。Mann在STFT和CWT的基础上引入chirp率，提出了chirplet变换(Chirplet Transform, CT) [15] ，该方法能够为线性调制信号提供高分辨率时频表征结果。然而，当遇到非线性信号时，CT的时频分辨率将有所降低。因此，Yu提出了广义线性Chirplet变换(General linear chirplet transform, GLCT) [16] ，该方法将CT中的参数视为局部时变变量，使中频近似于局部线性，并在每个点重新匹配参数。此外，通过引入广义频率旋转算子和平移算子，提出了多项式CT和样条CT，这些方法适合分析多项式相位信号。然而，目前的参数化时频分析方法的核函数被设计成时间变量的函数，可以有效地表征信号的非线性中频，适用于谐波类信号，但无法为时频曲线几乎平行于频率轴的瞬态信号提供高分辨率的时频表征结果 [17] 。

考虑到参数化时频分析方法在瞬态信号中的局限性，本文对原有的Chirplet变换进行修改，提出了一种广义线性频域Chirplet变换(General linear Frequency-domain chirplet transform, GLFCT)方法。该方法通过旋转时频基来匹配瞬态信号的实际GD曲线，从而为瞬态信号提供高分辨率的时频表征结果。在接下来的章节中，首先阐述了GLFCT的方法原理，并用合成信号和模拟地震道，对该方法的表征效果进行阐释，实验结果表明该方法在地震信号处理方面具有很大潜力。

2. 原理

对于信号，其短时傅里叶变换为：

$\text{STFT}\left(t,\omega \right)={\displaystyle {\int }_{-\infty }^{+\infty }x\left(u\right)g^{*}\left(u-t\right){\text{e}}^{-i\omega u}\text{d}u}$ (1)

其中t为时间， $\omega$ 为频率，u为时移变量，i为虚数单位， $g\in L\left(R^2\right)$ 表示一个非负对称的归一化时窗，即高斯窗函数为：

$g\left(t\right)=\frac{1}{\sqrt{2\pi }\sigma }{\text{e}}^{-\frac{1}{2}{\left(\frac{t}{\sigma }\right)}^2}$ (2)

其中， $\sigma$ 为标准差。式(1)在频域的对偶可写为：

$\text{SFFT}\left(t,\omega \right)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{+\infty }X\left(\eta \right)G^{*}\left(\eta -\omega \right){\text{e}}^{i\eta t}\text{d}\eta }$ (3)

其中， $\eta$ 为频移变量。式(3)称为短频傅里叶变换(short-frequency Fourier transform, SFFT)。 $X\left(\omega \right)$ 是 $x\left(t\right)$ 的傅里叶变换， $G\left(\omega \right)$ 是窗函数 $g\left(t\right)$ 的傅里叶变换。利用傅里叶变换，STFT和SFFT之间的关系可以推导为：

$\text{STFT}\left(t,\omega \right)={\text{e}}^{-i\omega t}\text{SFFT}\left(t,\omega \right)$ (4)

这表明除了因子 ${\text{e}}^{-i\omega t}$ ，STFT与SFFT是相同的。因此，二者有以下关系：

${\left|\text{STFT}\left(t,\omega \right)\right|}^2={\left|\text{SFFT}\left(t,\omega \right)\right|}^2$ (5)

与STFT类似，SFFT适用于分析具有恒定GD的信号。类似于SFFT的推导方式，CT的对偶可以表示为：

$V_X^G\left(t,\omega ,c\right)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{+\infty }X\left(\eta \right)G^{\text{*}}\left(\eta -\omega \right){\text{e}}^{i\left(\left(\eta -\omega \right)t+c\frac{{\left(\eta -\omega \right)}^2}{2}\right)}\text{d}\eta }$ (6)

式(6)为GLFCT的表达式。其中，c表示chirp率， ${\text{e}}^{ic\frac{{\left(\eta -\omega \right)}^2}{2}}$ 表示解调算子。由式(6)可知信号的相位函数为：

$\tau \left(t,\omega ,\eta ,c\right)=\left(\eta -\omega \right)t+c\frac{{\left(\eta -\omega \right)}^2}{2}$ (7)

由于信号的GD定义为相位函数对频率的导数的负值，由式(7)可得：

$\frac{\partial \tau \left(t,\omega ,\eta ,c\right)}{\partial \eta }=t+c\left(\eta -\omega \right)=-\text{GD}\left(\eta \right)$ (8)

GLFCT的核心思想是通过旋转时频基来匹配实际GD曲线，从而产生能量集中的时频表征。在长度为L、频率为 ${\omega }_c$ 的窗口中，GLFCT与目标GD曲线的匹配过程如图1所示。黑线是信号 $X\left(\eta \right)$ 的真实GD曲线，红线段定义为时频基。当时频基从平行于频率轴的位置绕点 $\left(X\left({\omega }_c\right),{\omega }_c\right)$ 旋转一个角度 $\alpha$ 时，它与整个窗长内的实际GD曲线匹配。

Figure 1. Schematic diagram of using GLFCT to match GD curves

图1. 用GLFCT匹配GD曲线示意图

对于每一个时频点 $\left(t,\omega \right)$ ，如果时频基与实际GD曲线完全匹配，那么GD曲线周围的时频表征结果能量较为集中，且振幅 $\left|V_X^G\left(t,\omega ,c\right)\right|$ 在所有值中达到最大值。然后对于每一个时频点 $\left(t,\omega \right)$ ，根据 $\left|V_X^G\left(t,\omega ,c\right)\right|$ 的振幅，我们可以得到最佳chirp率 $\hat{c}$ ：

$\hat{c}=\arg \underset{c}{\max }\left\{\left|V_X^G\left(t,\omega ,c\right)\right|\right\}$ (9)

下一个需要解决的问题是如何确定解调算子 ${\text{e}}^{ic\frac{{\left(\eta -\omega \right)}^2}{2}}$ 。根据Refs的研究结果，通过引入 ${\text{e}}^{ic\frac{{\left(\eta -\omega \right)}^2}{2}}$ 会对时频表征结果产生旋转效应。对于一个信号，采样时间 $T_s$ 和采样频率 $F_s$ 形成一个TF平面 $t\in \left(0,T_s\right)$ ， $f\in \left(0,\frac{F_s}{2}\right)$ ，所以可以在这个TF面内引入一个旋转参数 $\beta$ ：

$c=\tan \left(\beta \right)\cdot \frac{F_s}{2T_s},\beta \in \left(-\frac{\pi }{2},\frac{\pi }{2}\right)$ (10)

假设参数 $\beta$ 有 $N_c$ 个值，那么可以将TF平面按 $N_c$ 值平均分成 $N_c+1$ 个截面：

$\beta =-\frac{\pi }{2}+\frac{\pi }{N_c+1},-\frac{\pi }{2}+\frac{2\pi }{N_c+1},\cdots ,-\frac{\pi }{2}+\frac{N_c\pi }{N_c+1}$ (11)

根据式(10)，信号 $X\left(\eta \right)$ 的GLFCT公式可改写为：

$V_X^G\left(t,\omega ,\beta \right)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{+\infty }X\left(\eta \right)G^{*}\left(\eta -\omega \right){\text{e}}^{i\left(\left(\eta -\omega \right)t+\tan \left(\beta \right)\cdot \frac{F_s}{2T_s}\cdot \frac{{\left(\eta -\omega \right)}^2}{2}\right)}\text{d}\eta }$ (12)

3. 测试与分析

为验证GLFCT方法的有效性，将该方法应用于合成信号和模拟地震道的分析与处理。

3.1. 合成信号

首先本节将GLFCT应用到以下瞬态信号中，信号表达式如下：

$X\left(\omega \right)={\text{e}}^{\frac{0.01}{2\pi }\omega }\cdot {\text{e}}^{i2\pi \cdot \left(5\omega +2\pi \times \sin \left(\frac{0.0004}{\pi }{\omega }^2\right)\right)}$ (13)

该信号采样频率为100 Hz，信号持续时间为10 s。

图2给出了STFT、GLCT和GLFCT方法的时频表示结果。如图2(a)所示，STFT在GD曲线周围存在能量扩散现象，这是由于受到Heisenberg不确定性原理的限制，即时间分辨率和频率分辨率中的任何一个都只能以牺牲另一个为代价来增强。GLCT得到的时频表征结果如图2(b)所示。与STFT相比，GLCT产生的时频表征效果更差，这是因为GLCT将核函数设计成时间变量函数，通过估计信号的中频将能量集中在中频附近。因此，GLCT对具有慢时变瞬时频率的非线性信号具有良好的时频表征能力。然而，对于在时频域中GD曲线快速变化的瞬态信号，GLCT将提供模糊的时频表征结果。相较于STFT和GLCT，GLFCT通过旋转时频基匹配瞬态信号的GD曲线，从而产生了一个分辨率更高的时频表征结果，如图2(c)所示。综上所示，对于瞬态信号，GLFCT通过将能量集中在GD曲线附近，提供了高分辨率的时频表征结果。

Figure 2. Time spectrum calculated by (a) STFT, (b) GLCT and (c) GLFCT respectively

图2. 分别由(a) STFT，(b) GLCT和(c) GLFCT计算得到的时频谱

为了定量评价不同方法的性能，计算了这些时频表征结果的Rényi熵，并将其列在表1中。Rényi熵是常用的定量评价指标，其值越小表示时频能量聚焦性越好。从表中可以看出，GLCT因为产生了干扰现象，Rényi熵值最大；GLFCT方法的Rényi熵值是最小的，这说明在这些时频分析方法中，GLFCT能生成能量最聚焦的时频表征结果。

Table 1. Rényi entropy of time-frequency analysis method

表1. 时频分析方法的Rényi熵

3.2. 模拟地震道

Ricker小波被广泛用作地震信号模型，以提供地震波形特征的近似描述。通常，将其与反射系数序列进行褶积以生成模拟地震道。为了验证GLFCT在处理地震信号时的有效性和可行性，使用模拟地震道进行测试。该模拟地震道包含两组反射对(四个反射系数)，采样频率为1000 Hz。第一组反射系数的对应振幅分别为0.5和−0.5，第二组反射系数的对应振幅同样为0.5和−0.5。第一和第二组反射系数对之间的时间间隔分别为4和5 ms。将反射系数与主频分别为40、40、30、30的零相位Ricker小波进行褶积，从而得到了模拟地震道，如图3(a)所示。

对该模拟地震道分别使用STFT、GLCT和GLFCT来获得相应的时频表征结果，如图3(b)~(d)所示。从图3(b)中可以看出，由STFT表征的时频谱能量分布相对分散，这是由于Heisenberg不确定性原理导致的。GLCT时频结果中的两组反射对都出现了能量干扰现象，导致其不能有效地分离两组子波，时频谱出现模糊现象，如图3(c)所示。与STFT和GLCT的表征结果相比，GLFCT能清晰的区分出两组反射对，并且时频能量分布更为集中，如图3(d)所示。上述分析表明，在处理地震信号时GLFCT在一定程度上优于STFT和GLCT。

地震波在流体储层中传播时，其振幅会由于岩石骨架或孔隙流体分布不均等原因产生衰减，导致地震信号的能量减少，这主要表现在高频时地震信号能量显著降低，而低频时地震信号能量保持良好。时频分析方法能充分利用地震信号中所包含的有效信息，提高地震资料的利用效率，挖掘地震勘探资料中的有效地质信息，从而揭示不同尺度地质体下的不同频率响应特征，最终实现储层的精细表征。GLFCT通过提高地震信号的分辨率，从而挖掘地震资料中的有效信息，以期揭示地震信号异常时频响应规律，最终实现储层表征。

Figure 3. (a) Simulated seismic signal, time-frequency calculated by (b) STFT, (c) GLCT and (d) GLFCT respectively

图3. (a) 模拟地震信号，分别由(b) STFT，(c) GLCT和(d) GLFCT计算得到的时频

4. 结论

本文提出了一种处理瞬态信号的新时频分析方法——GLFCT，该方法通过旋转时频基来匹配实际GD曲线，并根据振幅值选取最佳chirp率，从而产生高分辨率的时频表征结果。合成信号分析结果表明，该方法能匹配瞬态信号的真实GD曲线，解决了GLCT无法为瞬态信号提供高分辨率时频表征结果的问题。对模拟地震记录的测试证实了GLFCT能够为地震信号提供能量高分辨率的时频表征结果，进而表明该方法在储层表征中的可行性。然而，GLFCT在地震资料处理中的应用需要更深入的研究，并且该方法的时频分辨率需进一步提高。