广义S变换与短时傅里叶变换在地震时频分析中的对比研究

黄 斌，张宏兵，王 强，乔 骁

（河海大学地球科学与工程学院，南京 211100）

广义S变换与短时傅里叶变换在地震时频分析中的对比研究

黄 斌，张宏兵，王 强，乔 骁

（河海大学地球科学与工程学院，南京 211100）

短时傅里叶变换（STFT）和广义S变换（GST）都被应用到地震时频分析中，但对两者在信号分析过程中的特点和差异的研究相对较少。通过比较两者的理论公式、窗口函数及地震信号的实际处理效果发现：短时傅里叶变换在地震信号分析过程中整个时频域具有相同的分辨率，整体性较强，缺少时频聚焦能力，不能对信号重点观测区域有针对性的提高时频分辨率；广义S变换对地震高频信号具有较高的时间分辨率，对低频信号具有较高的频率分辨率，且可以通过改变参数p的值对广义S变换窗口函数的形态做出较大的调整，也可以改变λ的值实现窗口形态微调，通过对窗口函数的调整，广义S变换可以对信号特定区域进行时频聚焦。

短时傅里叶变换；S变换；广义S变换；时频窗；地震记录

0 引言

地震勘探数据是一种非平稳的信号，其频率值随时间的变化较快，用传统傅里叶变换进行分析不能将时间值和频率值对应起来，分析效果不佳。采用短时傅里叶变换、S变换和广义S变换等时频分析方法则可以很清晰的反应出地震信号在每个时频点的幅值大小。S变换由R.G.Stockwel[1]于1996年提出，是在短时傅里叶变换和小波变换的基础上发展而来的一种新型时频分析方法。广义S变换提出后即在地震勘探数据处理领域中得到广泛的应用，刘喜武[2]等人用广义S变换来提取地震旋回，结论表明地震数据的振幅和频率的旋回性和沉积旋回具有很好的对应性；高静怀[3]将广义S变换用于研究薄互层地震响应特征，研究得出恰当选用基本小波，广义S变换在薄互层识别中将起到重要作用；陈学华[4]、王长江[5]等人用广义S变换实现了地震资料的高效时频谱分解；姬战怀[6]、齐春艳[7]研究和改进了广义S变换理论，使得广义S变换类型增多，应用范围更广。近年来在雷达信号[8]、图像信息处理[9]等领域也开始使用广义S变换且取得了较好的效果。

虽然短时傅里叶变换和广义S变换已经在地震信号时频分析上得到应用，但是关于两者在信号分析上的优缺点并没有系统的研究，本文首先由短时傅里叶变换理论推导到广义S变换理论，再研究两者窗函数的差别，进而研究广义S变换中两调节参数的具体作用，最后将两者应用到地震记录时频分析中，对比两者在地震信号分析中的特点，为实际数据处理提供参考。

1 广义S变换与短时傅里叶变换理论对比

1.1 从短时傅里叶（SFFT）到广义S变换（GST）

傅里叶变换是将信号的时域和频域联系起来的方法，但是傅里叶变换主要用于平稳信号的时频转换，是信号的一种整体变换，缺乏定位能力。为了解决这个问题，人们提出在时间域加窗的方法，即将时间域分割成不同的片段，对这些片段分别进行傅里叶变换得到对应的频率值，这样就实现了时间域和频率域的局部对应，将一维时间域转换成二维时频域，这也是短时傅里叶变换和广义S变换的基本原理。

其中g(τ -t)为移动窗口函数，τ为时移因子，（1）式是由傅里叶变换公式加了一个窗口函数构成。在时间域用g(t)窗函数去截取信号x(t)，可以认为截取的信号为平稳信号，对其做傅里叶变换得到t时刻附近的信号频域值，移动窗口即可得到时间序列不同t值附近的频域值，这些傅里叶变换的集合就是X(τ ，f)。X(τ ，f)是三维的复合函数，表示信号x(t)随时间和频率变化的相位和幅度值。

根据测不准原理，当窗函数为高斯类型函数时，短时傅里叶变换时频窗口半径最小。高斯窗可用下式表示：

观察（2）式可知短时傅里叶变换时移因子τ一旦确定，窗函数也就确定，其窗宽和幅值不会随着频率值改变。为了让窗函数和频率值联系起来，R.G. Stockwe[2]提出的S变换（ST）将（2）式中的σ与频率f相关联，令可得S变换的窗口函数为:

将（3）式带入（1）式就可以得到S变换的公式：

广义S变换又将S变换的窗函数加以改造，通过添加两个调节因子λ和p来达到调节窗函数值的目的，其窗函数为：

将（5）式带入（1）中即可得广义S变换的公式：

通过上述公式推导可知S变换是在继承短时傅里叶变换原理的基础上，通过采用高斯类窗函数来提高时频分辨率并且在窗函数中加入频率变量，使得窗口的宽度和幅值能随着频率值做线性变化，具有与频率相关的的分辨率，这就大大改善了短时傅里叶变换窗口在整个频率轴上固定不变的缺点。广义S变换是在S变换的窗口函数中引入调节参数，使得窗口函数随频率变化的规律可以调节，解决了S变换的时窗宽度和幅值随频率的变化值固定的这一缺点，使得窗口函数可以根据实际需要更加灵活的调控。

1.2 窗函数形态特征

图1a、图1b、图1c分别为短时傅里叶变换高斯窗，S变换窗，广义S变换窗。从图中可以明显看出：①短时傅里叶窗口函数窗宽和幅值沿频率轴没有发生改变。②S变换窗口函数幅值随频率的增大而增大，窗宽随频率的增大而降低，在高频处有高幅值和窄窗口，可以判定其高频处的时间分辨率很高，在低频处具有低幅值和宽窗口，可以得出其低频时间分辨率较低，根据测不准原理又可得其高频处的频率分辨率很低而低频处的频率分辨率则较高。③广义S变换窗函数形态大致和S变换窗函数相同，但是通过改变λ和p值可以发现其在高频处的幅值较S变换窗幅值大，窗宽则变得更窄，表明广义S变换通过参数调节使得其在高频处的时间分辨率较S变换提高了。

为了更具体的了解广义S变换窗口函数中调节参数λ和p值的调节特点，本文采用控制变量的方法，分别对单变量λ和p值做图分析。

首先分析λ和p值的改变对窗函数的幅值和宽的影响。固定p=1，λ分别取0.8、0.9、1.0、1.1、1.2，如图2a所示，观察可知随着λ的增大，窗的宽度逐渐减小，幅值逐渐增大，但是减小及增大的幅度不大。另外λ=1时为S变换的时窗曲线，可知若要在S变换的基础上减小窗宽和加大幅值可以适当增加λ的值，反之则适当减小λ值。

固定λ=1，p值分别取0.8、0.9、1.0、1.1、1.2，如图2b所示，观察可知随着p值的增大时窗幅值急剧增大，时窗宽度变窄的幅度也比较快。虽然效果和λ值变化时大致相同，但是p值的改变导致的时窗变化要比λ大的多。

通过对图1c的观察发现，改变λ、p的值对时窗的脊（时窗在各个频率值剖面最大值的连线）的形态也会发生改变。为了研究其变化特性，分别取p= 0.8，λ=0.8、1.0、1.2；p=1.0，λ=0.8、1.0、1.2；p=1.2， λ=0.8、1.0、1.2作窗脊幅值随频率f的变化曲线如图3所示。

图1 时频窗口Figure 1 Time-frequency window

图2 λ、ρ值的改变对窗口幅值和宽度的影响Figure 2 Impacting on window amplitude and width from variation of λ and ρ values

图3 窗脊的变化曲线Figure 3 Window ridge variation curve

观察图3可知p值固定时λ值的变化对曲线弯曲情况影响不大，对窗脊的幅值大小有一定影响，即λ越大，脊线的整体幅值越大。λ固定时，p值的变化对曲线的弯曲程度有较大影响，p＞1时，曲线向幅值减小方向凹陷，表明窗脊沿频率减小的方向幅值减小的速率加快；p=1时，窗脊曲线随频率的变化大致呈线关系性，即幅值沿频率轴变化的速率不变，表明S变换的时窗函数的幅值与频率值呈线性变化关系；p＜1时，窗脊曲线大体往幅值增大方向凸起，表明窗脊沿频率减小的方向幅值减小的速率变慢。p值的增大也会导致窗脊曲线的幅值整体增大，且增幅很大，比λ增大时带来的增幅更大。

通过以上对p，λ的分析可以得出，p值的改变对时窗的窗宽、幅值、窗脊的变化具有决定性的作用，λ对窗的形态改变远不及p值，在实际应用过程中可以起到对窗口微调的作用。

2 广义S变换与短时傅里叶变换在地震记录上的应用对比

2.1 地震记录模型的建立

在忽略各种能量衰减和干扰波的作用下，地震反射波记录x（t）可以看成是地震子波b（t）和地震反射系数w（t）的褶积加上噪声n（t）的结果，本文中的噪声选用正态分布的随机信号，其公式：

地震反射系数是上下地层波阻抗之间的函数，若上层波阻抗大于下层波阻抗，反射系数为负值，反之，反射系数为正值。本文采用的是Ricker子波模型[10]，模型公式：

式（8）中参数g＞0是子波的主频率[11]，本文g取频率为60 Hz，接近实际地震资料反射波主频。图4是本文模拟的地震信号图形，图4a为设置的10个反射系数，每一个反射系数值代表一个地层界面；图4b为ricker子波模型，图4c为其频谱图。可以看出ricker子波的主频为60 Hz，最高频约为150 Hz。图4d为合成的单道地震记录，采样间隔为2 ms，采样点数为300，时间从0～0.598 s。图4e为其频谱图，可知合成的地震记录在50～65 Hz以及80～100 Hz有较大的振幅值，最大频率约为250 Hz。

2.2 地震记录时频分析结果对比

对图4d所示的合成地震记录利用matlab编程软件分别对其做STFT、ST、GST时频分析，结果如图5所示。

图4 合成单道地震记录Figure 4 Synthetic single channel seismic records

图5a为短时傅里叶时频分析的结果，上文对其窗函数研究时提到，其窗函数形态与率值的大小无关，无法对特定的区域进行聚焦，所以其分析效果比较差，沿时间轴只能大概看出反射层的位置，对于相距很近且反射系数较小的地层如图4a所示的前三个反射系数值所代表的地层，在图5a中基本无法识别出来。图4e中的50～65 Hz及80～100 Hz两个高振幅频率区间对应STFT时频图中0.25～0.3 s，30～100 Hz的能量团，但是沿频率轴方向观察，能量团粘联在一起，无法区分两个高振幅频率区间，故其频率分辨率也不足。

由于ST法是GST法取λ=1，p=1时的特例，故也属于GST法，结果分析时统一看做GST法。图5b为ST法时频图，有效信号在0.1～0.4 s，30～120 Hz内，也不能得到特别好的时频分辨率，但是相对图5a来说已有所提升，尤其在低频频率分辨率上。有效信号区间外的一些浅色能量团为噪声信号，特点是噪声能量团比较分散不具连续性，但信噪比很低的地震信号，噪声会对有效信号的识别带来较大的影响。图5c为调节参数λ=1.5，p=1.1的GST时频图，由上文窗函数的讨论可知，这种条件下的窗函数较ST变换的窗函数在高频处的幅值要大，窗宽要窄，故其在信号高频处的时间分辨率比ST法要高。图5c中的每个条带状能量团代表一个地层，各条带在时间尺度上间隔比较明显，可见其时间分辨率比前面两图都要高。为了得到较好的频率分辨率将GST窗的调节参数适当降低，如图5d所示的时频图采用λ=0.8，p=0.8的窗函数，可知其频率分辨率较前面3个图的要高。前面3个图都无法分辨的50～65 Hz及80～100 Hz两个高振幅频率，在图5d中0.25～0.3 s区域沿频率轴方向被清楚的分成了2个高能量团，表明在这种参数条件下的GST频率分辨率高。图5b、图5c、图5d联合起来分析就可以比较细致的刻画出地震信号在时频域的能量分布，为后期滤波等处理提供技术支持。

图5 单道地震记录时频分析Figure 5 Single channel seismic records time-frequency analysis

3 结论

①通过对STFT及GST窗函数的对比研究可得：STFT的窗函数一旦固定，其时频分辨率在整个时频域上相同且无法改变；GST的窗函数在高频处具有幅值高，窗窄，在低频处具有幅值低，窗宽等特点，所以其在信号高频处具有高频率聚焦性能，在低频处具有高时间聚焦性能。通过对GST窗函数调节参数的研究表明：p值的改变对时窗的窗宽、幅值、窗脊的变化具有决定性的作用，λ对窗的形态改变远不及p值，在实际应用过程中可以起到对窗口微调的作用。

②通过对合成地震记录采用多种时频方法对比分析表明：STFT法在地震数据频谱分析过程中虽然具有一定整体性但缺乏时频聚焦能力，既不能得到高的频率分辨率也不能得到好的时间分辨率，在实际地震数据时频分析过程中可以作为一种参照方法。GST法虽然不能在一张时频图中同时得到高的时间分辨率和频率分辨率，但是可以改变参数值，分别得出地震信号的高频率分辨率的时频图和高时间分辨率的时频图，联合起来分析就可以较为清楚的了解地震信号在时频域的能量分布，在对时频分辨率要求较高的地震时频分析过程中，可以利用其刻画特定区域的时频能量分布。

[1]Stockwell R G，Mansinha L，Lowe R P.Localization of the complex spectrum:the Stransform[J].IEEE on Signal Processing，1996，44（4）: 998-1001.

[2]刘喜武，刘洪，李幼铭，等.基于广义S变换研究地震地层特征[J].地球物理学进展，2006，21（2）:440-451.

[3]高静怀，陈文超，李幼铭，等.广义S变换与薄互层地震响应分析[J].地球物理学报，2003，46（4）:526-532.

[4]陈学华，贺振华，黄德济.基于广义S变换的地震资料高效时频谱分解[J].石油地球物理勘探，2008，43（5）:530-534.

[5]王长江，杨培杰，罗红梅，等.基于广义S变换的时变分频技术[J].石油物探，2013，52（5）:489-494.

[6]姬战怀，严胜刚.关于“一种改进广义S变换”的讨论[J].石油地球物理勘探，2015，50（6）:1224-1230.

[7]齐春艳，李彦鹏，彭继新，等.一种改进的广义S变换[J].石油地球物理勘探，2010，45（2）:215-218.

[8]刘传武，张智军，毕笃彦.S变换在雷达目标识别中的应用[J].系统仿真学报，2008，20（12）:3290-3292.

[9]甄莉，彭真明.基于广义S变换的图像局部时频分析[J].航空学报，2008，29（4）:1013-1019.

[10]梅金顺，王润秋.Ricker类地震子波[J].石油地球物理勘探，2012（s1）:8-14.

[11]韩海英，王志章，王宗俊，等.基于Ricker类地震子波的匹配追踪[J].石油物探，2014，53（1）:17-25.

Contrastive Study of Generalized S Transform with Short-time Fourier Transform in Seismic Time-frequency Analysis

Huang Bin,Zhang Hongbing,Wang Qiang and Qiao Xiao
(School of Earth Science and Engineering,Hohai University,Nanjing,Jiangsu 211100)

The short-time Fourier transform(STFT)and generalized S transform(GST)are all applied to the seismic time-frequency analysis procedure,but without systematic study of difference and peculiarity between the two in signal analysis procedure.The paper through correlation of theoretical formulae,window functions and real processed effects of the two has found that the STFT during the seismic signal analysis procedure has identical resolution in whole time-frequency domain with better integrity,but lack of time-frequency focusing capability and can not enhance time-frequency resolution in key signal observation area.The GST has higher time resolution on seismic high frequency signals,higher frequency resolution on low frequency signals;and can make major adjustment on GST window functions through change of parameter p;also can change λ value to realize fine adjustment of window configuration; through adjustment of window functions,GST can carry out time-frequency focusing in signal specific areas.

STFT;S transform;GST;time-frequency window;seismic record

P631.4

A

10.3969/j.issn.1674-1803.2017.01.13

1674-1803(2017)01-0059-05

国家自然科学基金（41674113)（41374116）

黄斌（1990—），男，江西抚州人，河海大学，地质资源与地质工程专业，硕士，研究方向为地震勘探数据处理。

2016-12-05

责任编辑：孙常长