1. 引言
地震信号分析的核心目标之一是准确提取地震信号的时频特征,而由于地下介质对地震波会有一定的衰减和散射效应,导致其振幅和相位会随时间变化而变化,所以地震信号往往是时变的、非平稳的[1]。因此,需要对地震信号进行分析。时频分析能很好的描述信号的频率成分是如何随时间变化而变化的,也能很好地刻画信号的局部时频特征[2]。
常见的时频分析方法有短时傅里叶变换(STFT) [3]、Wigner-Ville分布(WVD) [4]、连续小波变换(CWT) [5]、S变换(ST) [6]等。短时傅里叶变换首次实现了信号的时频局部化表征,而它受限于海森堡不确定性原理,导致时间和频率分辨率无法同时优化;Wigner-Ville分布作为最早的非线性时频分析工具,虽具有理想的时频聚集性,但存在严重的交叉项干扰问题;连续小波变换具有多尺度分析能力,不过小波基的选择需依赖先验知识;S变换结合了STFT和CWT的优点,时频分辨率自适应变化,且无需窗函数选择。但目前S变换仍然存在以下问题:一是S变换存在时频谱中主频向高频方向偏移的现象,在对地震信号做时频分析后会出现假主频,导致地层定位不准确;二是S变换虽然能实现自适应的调整时频谱在不同频率的时频分辨率,但它的窗函数过于单一,形状变化固定,不能进一步的提升时频效果。
目前为了解决主频偏移现象,最大的改进就是提出了去尺度化变换。在改进S变换的过程中,Wu和Castagna移除了S变换中有关频率的线性项,从而来构建去尺度化S变换算法[7];Zhang等设计了一种去尺度化的改进S变换的频域形式,并将之用于地震储层定位[8];刘乃豪等利用S变换的多尺度多分辨率特性,对不同的实际地震数据进行了处理和分析[9]。然而,去除频率相关的线性项会导致低频的时频分辨率降低,因此还需要考虑如何去优化时频分辨率。
针对这第二点,现在很多学者在S变换的基础上对其进行改进,提出了一系列新方法。高静怀等针对地震信号的特点,修改了基本小波,并提出了对应的广义S变换,有效地改善了时频分辨率[10];陈学华等通过引入两个可变参数来改进高斯窗函数,提升复杂信号表征能力[11];Liu等在此基础上添加了一个常数项,以此来优化频率分辨率,能更好区分密集频段成分[12];Li等采用四参数实现窗函数的自适应调整,可以同步优化时频分辨率,有利于更精准表征快速瞬变信号[13]。
在这篇文章中,将去尺度化与广义S变换结合,提出了一种四参数调控的去尺度广义S变换(QCDGST)。一方面,文中所提出的变换通过引入了四个可变参数来更灵活的控制窗函数,进一步提高时频分辨率;另一方面,该变换仍然具有去尺度化S变换(UST)的频率特性。利用该方法对单个子波、合成地震记录以及Marmousi2模型进行了分析,结果显示在解决主频偏移的基础上能有效地提高时频分辨率。
2. 基本理论
2.1. S变换
Stockwell在1996年发展的S变换,创新性地引入频率相关窗函数[6],其数学表达式为:
其中,
为时频谱矩阵,
表示待分析信号,t和f分别表示时间和频率向量,
控制窗口函数在时间轴上的位置。该变换结合了STFT的时间–频率局部化特性和CWT的多尺度、多分辨率分析能力,从短时傅里叶变换来看,S变换利用形状可变的窗口截取信号进行分析;从连续小波变换来看,S变换实际上就是将小波的基函数用高斯窗函数来代替,而且高斯窗口可以克服了短时傅里叶变换窗口高度和宽度固定的缺陷,其会随频率而变化。
高斯窗函数表达式为:
可以看出,S变换中的窗函数引入了一个与频率成反比关系的尺度因子,当频率比较小时,窗口宽度大,频率分辨率较好;当频率较大时,窗口宽度小,时间分辨率好。
标准S变换在一定程度上改进了传统时频分析中窗口固定、分辨率单一的问题,其窗口函数与频率呈线性关系而变化。但由于这种关系的限制,当处理结构更为复杂的信号时,ST在时频分辨率方面仍难以达到理想效果。
2.2. 广义S变换
为了优化S变换的时频分辨率,高静怀团队在2003年提出了广义S变换,这是一种S变换的改进方法,它通过引入自适应调节因子,让窗口不再局限于随频率线性变化,而是可以依据信号的不同来改变[10],其改进的时频分布表示为:
其中,
为频率自适应调节因子。
本文在该研究的基础上,提出对调节因子进行修改,通过增加四个可变参数进一步改进高斯窗函数,以让其对不同变化的信号具有更一般的适应性。定义改进的窗函数为:
其中,改进的标准差函数为:
由此可以得到该窗口下的S变换:
在该变换中,选择的窗口宽度不宜过大或过小,因此我们首先需要设定四个参数的大概范围,窗宽需要满足:
其中,最小为nTs,最大为lTs,Ts是采样周期,理论上,n和l为整数,且n要大于2,根据经验和大量实验,这里n选择为10,l为1000,因此对于采样频率为1000Hz的信号,窗口最大为1秒。
因此可以给出四个参数满足的范围,k,p,m,n的下限和上限分别设置为0和2。
此外,新提出的窗口满足S变换的窗口归一化条件:
这一条件可以保证该广义S变换也是可逆的。
对高斯窗函数的改进是在动态标准差的分子上添加了两个分量
和
,分母上将
修改为
和
,它们共同决定时频分辨率的缩放。当
时,该变换与标准的S变换等价,为了说明该改进的有效性,图1展示了不同参数的取值对窗口大小的影响。
四个参数
可以根据信号自适应地调整窗口大小,在本文中,这里的四个参数需要手动确定,我们可以根据实际信号的需求,合理地去调整参数的取值,图2(a)~(d)分别分析了当频率f = 10 Hz时不同参数对高斯窗口形状的影响,从中可以知道它们对于改善时频分辨率的作用及意义:
(1) 参数k通过调控线性项改变窗口大小,随着k值的增大,窗函数的宽度会减小,它的变化对低频段的影响比较大,因此改变它的值可以确保低频时窗口不会过宽;
(2) 参数p以
的形式位于分母上,可以有效地降低高频段的影响,随着p值的增大,高频段的窗函数的宽度会减小,它对低频段的几乎没有影响;
(3) 参数m引入了指数项
,主要影响的是高频段的分辨率;
(4) 参数n作为分子上的线性项参数,在高频段,它的影响虽然比不过
,但它可以确保当频率较小时,窗口不会过宽。
Figure 1. Effect of parameters on window size
图1. 参数对窗口大小的影响
Figure 2. Effects of different parameters on window shape. (a) Effect of parameter k (p = 1, m = 1, n = 0); (b) Effect of parameter p (k = 0, m = 1, n = 0); (c) Effect of parameter m (k = 0, p = 1, n = 0); (d) Effect of parameter n (k = 0, p = 1, m = 1)
图2. 不同参数对窗口形状的影响。(a) 参数k的影响(p = 1, m = 1, n = 0);(b) 参数p的影响(k = 0, m = 1, n = 0);(c) 参数m的影响(k = 0, p = 1, n = 0);(d) 参数n的影响(k = 0, p = 1,m = 1)
2.3. 去尺度广义S变换
S变换通过与频率绝对值成正比的因子对频率相关高斯窗口进行归一化,这种归一化会使时频谱的主频偏向更高的频率,Wu和Castagna在此基础上提出去尺度化S变换(UST),删除了窗函数中的一个尺度因子[7],该UST可以表示为:
虽然它能有效地改善主频偏移的现象,但是它的窗口变化单一、时频分辨率不足,不能像广义S变换那么灵活的选取窗口大小,这就会导致时频效果不好,那么在确定主频时也会受到其影响,会出现主频仍有细微的偏移。
本文充分考虑到标准S变换与去尺度S变换的优缺点,将去尺度化与2.2节提出的广义S变换相结合,利用二者的优点进行互补,改善两者的不足,在保证主频不再偏移的情况下,尽可能地去改善时频分辨率,将这种变换称之为四参数调控的去尺度广义S变换(QCDGST),它的数学表达式为:
由2.2节中提到的窗函数归一化条件,我们知道,QCDGST也满足该条件,因此这确保了QCDGST也是可逆的,它的逆变换为:
3. 模拟数据及结果分析
3.1. 单个雷克子波
为了说明传统时频分析会存在主频向高频偏移的现象,本小节将比较不同时频分析方法对单个雷克子波的处理结果,采用的子波主频为30 Hz,中心位置为0.5 s,截取的部分为0.3 s~0.7 s,如图3(a),分
Figure 3. Ricker wavelet with a dominant frequency of 30 Hz and a center position of 0.5 s and its time-spectrum. (a) Ricker wavelet; (b) STFT time-spectrum; (c) ST time-spectrum; (d) UST time-spectrum; (e) AGST time-spectrum; (f) QCDGST time-spectrum
图3. 主频为30 Hz,中心位置0.5 s雷克子波及其时频谱。(a) 雷克子波;(b) STFT时频谱;(c) ST时频谱;(d) UST时频谱;(e) AGST时频谱;(f) QCDGST时频谱
别对其做短时傅里叶变换(STFT)、标准S变换(ST)、去尺度S变换(UST)、自适应广义S变换(AGST)以及本文的变换,得到的结果如图3(b)~(f)所示,其中红色虚线代表频率为30 Hz的位置。从该结果中可以看出,STFT虽然没有发生偏移,但由于窗口形状固定,它的时频分辨率在低频段和高频段呈现一样的分辨率,无法区分地震信号的频率特性;ST和AGST都很明显地出现了主频向高频段偏移的现象;UST和QCDGST很好地避免了这个问题,而且本文方法通过引入可变的四参数修改高斯窗函数,可以获得更高的时频分辨率。
另外,对该雷克子波做傅里叶变换,得到如图4所示的振幅谱,可以看出,雷克子波的主频在30 Hz然后在每个变换的时频结果中,固定时间t = 0.5 s,提取该时间点处的振幅谱,得到的对比结果见表1。从表中数据可以看出,对于主频为30 Hz的雷克子波来说,ST和AGST有明显的主频偏移现象;UST和QCDGST能很好地改善这一问题,以便能更精确地定位主频位置。
Figure 4. Ricker wavelet amplitude spectrum
图4. 雷克子波振幅谱
Table 1. Comparison of the main frequencies estimated by different time-frequency analysis methods
表1. 不同时频分析方法估计的主频对比(单位:Hz)
时频分析方法 |
STFT |
ST |
UST |
AGST |
QCDGST |
估计主频(Hz) |
31.06 |
37.08 |
30.06 |
37.07 |
30.06 |
3.2. 合成地震记录
现在以3.1节中的雷克子波为基准,合成多层地震记录,分别设置雷克子波主频为15 Hz、30 Hz、45 Hz,时间为0到1 s,由于实际地震会有振幅衰减,所以该实验中设置了子波的振幅随时间增加而减小。其中,主频为15 Hz的雷克子波分布在时间为0.8 s处,记为信号
,主频为30 Hz的雷克子波分布在时间为0.4 s、0.6 s处,记为信号
,主频为45 Hz的雷克子波分布在时间为0.1 s、0.25 s,记为信号
,然后将三者相加构建一个地震记录
。
图5(a)给出了合成地震记录的时域表现形式,图5(b)~(f)分别是STFT、ST、UST、AGST、QCDGST的时频结果,为了能有最优的时频效果,本文多次尝试不同的参数取值,并依据大量的实验,选取这里的参数[k, p, m, n] = [1.5, 1.5, 0.6, 0.1]。首先从时频分辨率的角度来分析,不管频率怎么变化,STFT的结果是不变的,它对高低频段的分辨率都是一样;对于标准S变换,由于添加了尺度因子,高频时窗宽明显变小,分辨率更高了,但低频时窗口过宽,导致时间分辨率比较差;UST在S变换的基础上将窗函数的一个尺度因子去除了,但较S变换,它的结果在低频处分辨率会更差;与传统的时频分析方法相比,AGST不管高频还是低频,它的分辨率都要更好一些;从QCDGST的结果中可以发现,低频段几乎没有信号与信号之间的交叉,而且高频段的拖尾也较短,说明时频定位更准确。表2给出了各个时频分析结果的CM值,这个值越大,说明效果越好,CM值的数学定义如下:
该变换归一化后的结果为:
Figure 5. Synthetic seismic record and its time-frequency analysis results. (a) Synthetic seismic record; (b) STFT time-frequency spectrum; (c) ST time-frequency spectrum; (d) UST time-frequency spectrum;(e) AGST time-frequency spectrum; (f) QCDGST time-frequency spectrum
图5. 合成地震记录及其时频分析结果。(a) 合成的地震记录;(b) STFT时频谱;(c) ST时频谱;(d) UST时频谱;(e) AGST时频谱;(f) QCDGST时频谱
Table 2. Comparison of CM values corresponding to different time-frequency analysis methods
表2. 不同时频分析方法对应的CM值比较
时频分析方法 |
STFT |
ST |
UST |
AGST |
QCDGST |
CM值 |
0.0082 |
0.0084 |
0.0080 |
0.0090 |
0.0098 |
其次从主频是否偏移的角度来看,现在已经有了时频谱,因此可以固定时间t = 0.1 s,取该时间点处的振幅谱进行比较,振幅谱如图6所示,其中红色虚线代表的是频率为45 Hz的位置。
Figure 6. Comparison of amplitude spectra of different methods at t = 0.1 s (45 Hz is marked as a red dashed line)
图6. 不同方法在t = 0.1 s处的振幅谱对比(45 Hz标记为红色虚线)
表3给出了该位置的估计主频,ST和AGST明显地发生了主频偏移,它们的估计主频分别为56.23 Hz、55.22 Hz;从STFT和UST的振幅谱来看,它们的主频都不是在45 Hz这个地方,而是有移动迹象;对于本文提出的方法QCDGST,它的估计主频为45.18 Hz,较其他方法最为接近真实主频。
Table 3. Comparison of the main frequencies estimated by different time-frequency analysis methods
表3. 不同时频分析方法估计的主频对比(单位:Hz)
时频分析方法 |
STFT |
ST |
UST |
AGST |
QCDGST |
估计主频(Hz) |
46.18 |
56.23 |
44.18 |
55.22 |
45.18 |
该模拟数据说明本文提出的方法不管是在改善时频分辨率上,还是优化主频偏移问题上,都有比较好的效果,进而说明了该方法的可行性。
4. Marmousi2模型
在本节中,我们采用Marmousi2模型进行对比分析,采样率为0.02 s,如图7(a)所示,现在从该模型中提取某一道数据(比如第450道),其中黑色虚线就代表第450道,从图中明显可以看出,在时间为2 s、4.5 s、9 s、11 s附近都有反射层,信号在这些地方波动很大,图7(b)是第450道地震数据的时域形式,现在对第450道地震数据进行时频分析,截取8 s与10 s之间的信号,并对比其时频分辨率和主频位置变化。
图8是每个时频分析方法对该部分地震数据处理的结果,这里的参数取值和3.2节类似,通过取不同的参数值来确定时频分析的效果,选取效果最好的那组参数,此时参数[k, p, m, n] = [1.5, 1.7, 1.1, 0.05]。我们注意到,该道的地震信号在8 s与10 s之间存在两个反射层,因此在截取的信号中,分别在8.8 s和9.4 s附近的能量就会大一些,因此时频分析的好坏就取决于是否能够清晰地分开这两个能量。从实验中发现,五种方法中,STFT、ST、UST的时间分辨率都比较差,AGST和QCDGST的分辨率比较高,但是AGST的时频分布仍比较分散些,它的拖尾现象较严重,所以本文提出的方法在平衡分辨率这方面更有优势。
Figure 7. Marmousi2 model. (a) Marmousi2 model cross-section; (b) Selected seismic signal No. 450
图7. Marmousi2模型。(a) Marmousi2模型剖面图;(b) 选取的第450道地震信号
Figure 8. Time-frequency analysis results of 8 s~10 s signals using different methods. (a) Cropped signal; (b) STFT time-frequency spectrum; (c) ST time-frequency spectrum; (d) UST time-frequency spectrum;(e) AGST time-frequency spectrum; (f) QCDGST time-frequency spectrum
图8. 不同方法对8 s~10 s信号的时频分析结果。(a) 截取的信号;(b) STFT时频谱;(c) ST时频谱;(d) UST时频谱;(e) AGST时频谱;(f) QCDGST时频谱
另外,通过固定时间为8.8 s,对比各个方法估计的主频发现,ST和AGST仍然存在主频向高频偏移的现象,而UST和QCDGST估计的主频都是在真实主频附近,虽然STFT的结果也很接近,但分辨率太差了,它的时频谱并不能用来解决实际问题。具体的估计频率见表4。
Table 4. Comparison of the main frequencies estimated by different time-frequency analysis methods
表4. 不同时频分析方法估计的主频对比(单位:Hz)
时频分析方法 |
STFT |
ST |
UST |
AGST |
QCDGST |
估计主频(Hz) |
4.58 |
5.95 |
4.44 |
5.80 |
4.51 |
5. 结论与展望
5.1. 结论
本文针对传统时频分析方法的窗口形状固定、变换单一和S变换中主频偏移等问题,提出了一种新的时频分析方法QCDGST。该方法通过在高斯窗函数中添加四个可变参数,可以让窗口的形状更多变,以此来适应不同的地震信号。通过对地震数据的分析可以证明本文方法是有效的。实验表明,本文方法在时频聚焦性、分辨率具有更好的性能,对高斯窗函数的修改可以让时频能量更集中,减少信号的能量发散、拖尾等现象。然后在此基础上引入去尺度化,可以很好地解决S变换及其改进方法中的主频向高频偏移问题。通过分析合成地震数据和Marmousi2模型,我们知道,去除高斯窗函数的一个尺度因子可以保证主频不会偏移,从而能从地震信号中提取到更为精确的时频特征。另外,通过计算信号重建的均方根误差,与其他方法进行对比,QCDGST的重构误差是最小的,说明它的逆变换能很好地保留信号原始特征,进一步证明了该变换要优于其他方法。
5.2. 展望
本文所提出的方法虽然能够改善时频分辨率,但由于引入了四个可变的参数,这些参数会极大地影响最终分析结果,因此对于不同的信号来说,选取什么样的参数才能有最优的分辨率,所以根据信号的不同自适应的选择参数是一个待研究的问题。