油气勘探

地震数据时频分析的W变换及其改进方法

  • 王仰华 , 1, 2 ,
  • 饶莹 , 3 ,
  • 赵振聪 1, 3
展开
  • 1 英国帝国理工大学资源地球物理学院,英国伦敦 SW7 2AZ
  • 2 中国地质大学(北京),北京 100083
  • 3 中国石油大学(北京)地球物理学院,北京 102249
饶莹(1979-),女,江西南昌人,中国石油大学(北京)教授,主要从事地震波场模拟、地震波形反演与油气藏特征描述等研究工作。地址:北京市昌平区,中国石油大学(北京)地球物理学院,邮政编码:102249。E-mail:

王仰华(1963-),男,江苏盐城人,英国帝国理工大学终身教授,英国皇家工程院院士,中国工程院外籍院士,主要从事应用地球物理研究。地址:英国伦敦,帝国理工大学资源地球物理学院,邮政编码:SW7 2AZ。E-mail:

收稿日期: 2024-04-19

  修回日期: 2024-05-30

  网络出版日期: 2024-08-02

基金资助

国家自然科学基金(42055402)

The W transform and its improved methods for time-frequency analysis of seismic data

  • WANG Yanghua , 1, 2 ,
  • RAO Ying , 3 ,
  • ZHAO Zhencong 1, 3
Expand
  • 1 Resource Geophysics Academy, Imperial College, London SW7 2AZ, UK
  • 2 China University of Geosciences (Beijing), Beijing 100083, China
  • 3 College of Geophysics, China University of Petroleum (Beijing), Beijing 102249, China

Received date: 2024-04-19

  Revised date: 2024-05-30

  Online published: 2024-08-02

摘要

针对常规线性时频分析方法无法同时达到时间域和频率域上的高分辨率和能量聚焦,尤其在低频区域分辨率较差的问题,为提高线性时频分析方法在低频区域的分辨率,提出了W变换方法,在线性变换中引入瞬时频率参数,构建与地震数据瞬时频率相匹配的分析时窗。将W变换方法与典型的非线性时频分析方法——魏格纳-维利分布(WVD)进行了对比。WVD方法展示的是时间-频率域中的能量分布,明确指示子波的时间重心和频率重心,而W变换中由于引入了任意时间位置所对应的瞬时频率作为变换参数,所展示的时间-频率谱因而也具有明确的能量聚焦重心,因此,W变换与WVD方法具有直接对标意义。分析了近年来W变换的3种改进方法的发展状况,详细阐述从常规W变换、谐变时窗W变换、分数阶W变换、直到线性正则W变换的发展和演进。通过W变换在河道砂体识别和溶洞检测方面的3个应用实例,验证W变换可以提高时间-频率谱的分辨率与能量聚焦。

本文引用格式

王仰华 , 饶莹 , 赵振聪 . 地震数据时频分析的W变换及其改进方法[J]. 石油勘探与开发, 2024 , 51(4) : 774 -782 . DOI: 10.11698/PED.20240260

Abstract

The conventional linear time-frequency analysis method cannot achieve high resolution and energy focusing in the time and frequency domains at the same time, especially in the low frequency region. In order to improve the resolution of the linear time-frequency analysis method in the low-frequency region, we have proposed a W transform method, in which the instantaneous frequency is introduced as a parameter into the linear transformation, and the analysis time window is constructed which matches the instantaneous frequency of the seismic data. In this paper, the W transform method is compared with the Wigner-Ville distribution (WVD), a typical nonlinear time-frequency analysis method. The WVD method that shows the energy distribution in the time-frequency domain clearly indicates the gravitational center of time and the gravitational center of frequency of a wavelet, while the time-frequency spectrum of the W transform also has a clear gravitational center of energy focusing, because the instantaneous frequency corresponding to any time position is introduced as the transformation parameter. Therefore, the W transform can be benchmarked directly by the WVD method. We summarize the development of the W transform and three improved methods in recent years, and elaborate on the evolution of the standard W transform, the chirp-modulated W transform, the fractional-order W transform, and the linear canonical W transform. Through three application examples of W-transform in fluvial sand body identification and reservoir prediction, it is verified that W-transform can improve the resolution and energy focusing of time-frequency spectra.

0 引言

时间-频率域谱分析是非稳态信号处理与分析的主要技术手段与工具,被广泛应用于地震数据的分析与处理中[1]。时间-频率域谱分析方法通过将一维时间域地震数据拓展到二维的时间-频率域,在高维空间中分析和处理地震数据,进而更有效地表征地震数据中的信息,指导油气藏的勘探与开发。
根据基本方法不同,常规时频分析方法可分为线性时频分析方法和非线性时频分析方法两种。常用的线性时频分析方法包括:短时傅里叶变换[2]、小波变换[3]、S变换[4]、W变换[5]等。常用的非线性时频分析方法包括魏格纳-维利分布(Wigner-Ville Distribution,WVD)[6]、匹配追踪技术[7]、基于频率重排的时频分析方法[8]等。
关于非线性时频分析方法,地震数据处理与分析中应用最为广泛的方法之一为WVD[6]。WVD可视为复数域信号与其共轭信号之间的瞬时互相关,获得的时间-频率谱分辨率与能量聚集度较高。但是不同信号成分之间的互相关会在时间-频率谱中造成大量信号分量之间的串扰。多种改进方法目标都是为了压制信号分量之间的串扰。比如,Choi和Williams提出的平滑伪WVD(SP-WVD)[9] 可以有效压制串扰;Wang等提出的多道最大熵WVD方法同时利用多个WVD核函数消除多种成分之间的串扰[10];Rao等提出利用多道最大熵WVD估算地震数据主频,用于识别地层中的岩溶空洞,提高地震异常体识别的精度[11];Alsalmi和Wang提出的掩码WVD方法,构建稳定的掩码滤波器对WVD时频谱进行后处理[12]。上述这些改进方法都可视为对时间-频率谱进行滤波处理,无疑会滤除一些有效能量。
WVD方法中的信号分量串扰问题也可以通过匹配追踪技术解决[13]。Wang将匹配追踪引入到地震数据处理和分析领域,并提出利用Morlet子波匹配地震道时的解析时频谱[7]。为了提高匹配追踪的可靠性,Wang进一步提出了多道匹配追踪技术[14],利用相邻道之间的平均地震道作为标准道,在提高时间-频率谱横向连续性以及压制随机噪声的同时,提高了匹配追踪的收敛效率。Zhao等提出了结构自适应的多道匹配追踪算法[15],进一步提高了匹配追踪技术在复杂构造地区地震数据处理中的适应性。以匹配追踪为基础的WVD方法(MP-WVD),利用匹配追踪算法将地震数据中的每个信号分量提取出来,再单独生成时间-频率谱,这样就避免了不同分量之间的串扰。但是,匹配追踪通过构建超完备的子波字典,将待分析信号与超完备字典中的子波进行匹配,并不断迭代直到满足终止条件。因此,该计算方法非常耗时。
综上所述,WVD方法是典型的非线性时频分析方法,因为其非线性特性,导致了地震信号不同分量之间的串扰;不同改进方法试图解决这个问题,一来会导致分辨率的降低,二来会大大增加运算耗时。
线性时频分析方法主要利用待分析地震信号与特定窗函数的褶积求取时间-频率谱,表达形式简单,计算效率较高,因此在地震数据处理与分析中具有较为广泛的应用。线性时频分析变换的研究主要集中于如何选取最佳的时窗分析函数,以达到最佳时频分析效果。
线性时频分析方法因其线性特性,不存在非线性时频分析方法中的串扰问题。但是,传统的线性时频分析方法无法同时达到时间域和频率域上的高分辨率和能量聚焦,尤其在低频区域其分辨率较差。
为了解决线性时频分析方法中的传统变换方法在低频区域分辨率较低导致油气藏识别精度不高等问题,Wang提出了W变换(即王氏变换)[5]。该方法利用地震数据的瞬时主频改进传统线性变换中的窗核函数,并提出了改进地震数据主频估算的表达形式,大大提高了线性时频分析方法在油藏描述等方面的应用。自W变换方法提出以来,相关研究人员开展了一系列针对W变换窗核函数的改进方案研究。Li等针对常规W变换在主频附近的导数奇异性问题,将窗核函数表达成多项式形式,提出了广义W变换[16],提高了W变换在主频附近的分辨率。同时,Chen等提出3参数W变换以提高地震数据时间-频率谱分辨率[17],该方法与广义W变换的思想相同,通过引入3个辅助参数,控制平滑Gaussian窗的时间与频率分辨率。Luo和Zong将同步提取变换思想引入W变换中,提出了同步提取W变换[18],改善了W变换在薄层与河道检测中的应用效果。
以上方法均通过改进W变换窗核函数的时间或者频率尺度,提高时间-频率谱的分辨率,但对于强非稳态性地震数据的适用性仍然较差。基于该限制,笔者团队通过对强非稳态地震数据的分析,提出了谐变调制窗的W变换方法,改进时频分析窗核函数时间-频率域中的旋转方向,进一步提高时频分析的分辨率与能量聚焦能力[19]。另外笔者团队还进一步对W变换进行了广义化,结合分数阶傅里叶变换与W变换,提出分数阶W变换,通过旋转时间-频率平面选取地震数据的最优时间-频率谱表征[20]。将上述时间-频率域中的旋转关系进一步广义化,笔者团队提出线性正则W变换[21]。实际地震数据的应用表明,这些W变换的改进方法可以有效进行基于地震数据属性的油藏描述。
本文主要阐述W变换的原理及其改进方法的发展状况,并与典型的非线性时频分析方法——WVD方法对比,对比W变换与WVD的应用效果差异。常规WVD方法展示的是时间-频率域中的能量分布,明确指示子波的时间重心和频率重心,因此对W变换具有直接对标意义。W变换中引入任意时间位置所对应的瞬时频率作为变换参数,所展示的时间-频率谱因而也具有明确的能量聚焦重心。如何实现能量聚焦,保证时间-频率谱的高分辨率则是针对W变换的不同改进方法的思想初衷。另外,本文还给出了W变换在河道砂体检测和溶洞检测方面的3个应用实例。

1 W变换原理

线性时频分析方法的核心思路为:构建一个时频分析的时窗核函数,利用待分析信号与时窗核函数之间的褶积关系来表征待分析信号的时间-频率谱。W变换是一个典型的线性时频分析方法,可表征为[5]
$WT\left( t,f \right)=\frac{\left| {{f}_{0}}\left( \tau \right)+\Delta f \right|}{k\sqrt{2\text{ }\!\!\pi\!\!\text{ }}}\int_{-\infty }^{\infty }{x(t){{\text{e}}^{-\frac{{{\left( t-\tau \right)}^{2}}{{\left[ {{f}_{0}}\left( \tau \right)+\Delta f \right]}^{2}}}{2{{k}^{2}}}}}{{\text{e}}^{-\text{i}2\text{ }\!\!\pi\!\!\text{ }f\tau }}\text{d}\tau }$
式中${{f}_{0}}(t)$为时间$t$处的瞬时频率。有关瞬时频率的估算及其改进公式请参考文献[5]。
(1)式中时窗核函数为:
$w\left( t \right)={{\text{e}}^{-\frac{{{t}^{2}}{{\left[ {{f}_{0}}\left( t \right)+\Delta f \right]}^{2}}}{2{{k}^{2}}}}}$
通过调节$k$可以控制时间-频率谱分辨率。
分析(2)式中W变换的时窗核函数可以看出,在时间位置$t$,不同频率所对应的时窗核函数以${{f}_{0}}\left( t \right)$为重心,当$f<{{f}_{0}}\left( t \right)$$f>{{f}_{0}}\left( t \right)$时,时窗尺度逐渐变窄,因此频率域的频带逐渐变宽,表明其在低频和高频部分都具有较高的时间分辨率。
如果不考虑待分析信号的瞬时频率,即${{f}_{0}}\left( \tau \right)=0$,W变换可以蜕变为S变换[4]
$ST\left( t,f \right)=\int_{-\infty }^{\infty }{x\left( t \right){{\text{e}}^{-\frac{{{\left( t-\tau \right)}^{2}}{{f}^{2}}}{2{{k}^{2}}}}}{{\text{e}}^{-\text{i}2\text{ }\!\!\pi\!\!\text{ }f\tau }}\text{d}\tau }$
其中时窗核函数为:
$w\left( t \right)={{\text{e}}^{-\frac{{{t}^{2}}{{f}^{2}}}{2{{k}^{2}}}}}$
从(4)式可以看出,在时间位置$t$,S变换的时窗长度随着频率而变化。当频率较高时,时窗尺度较短,但是,当频率较低时,时窗变长。因此S变换在低频区域的时间分辨率较低,影响其在油藏描述中的应用。
图1展示的是S变换与W变换时窗对比图,待分析信号的主频为40 Hz。正如(4)式所示,S变换的时窗尺度跟频率成反比(见图1a)。W变换时窗在低于主频40 Hz时(见图1b),明显比S变换的时窗窄,表明其在低频部分具有更高的时间方向上的分辨率。
图1 S变换与W变换时窗对比
图2a构建了含有不同延迟、不同频率的模拟信号,该模拟信号由4个不同主频的地震子波组成。图2b显示,S变换的时间-频率谱在低于子波主频的区域,其时间分辨率逐步降低。而W变换的提出摆脱了线性时频分析方法在低频区域分辨率低的困境(见图2c),可以大大提高复杂油藏的描述能力,并在工程地球物理中具有较大的应用前景和潜力。
图2 S变换与W变换时间-频率谱对比
上述W变换由于引入任意时间位置所对应的瞬时频率作为变换参数,所展示的时间-频率谱具有明确的能量聚焦重心。WVD方法[6]展示的是时间-频率域中能量分布,也明确指示了子波的时间重心和频率重心。虽然WVD方法是一种非线性时频分析方法,但是它对W变换方法具有直接的对标意义。
地震数据$x(t)$的WVD可以表征为:
${{S}_{\text{WVD}}}(t,f)=\int_{-\infty }^{\infty }{z\left( t+\frac{\tau }{2} \right){{z}^{*}}\left( t-\frac{\tau }{2} \right){{\text{e}}^{-\text{i}2\text{ }\!\!\pi\!\!\text{ }f\tau }}\text{d}\tau }$
上式中$z(t)$是地震数据$x(t)$的解析信号,不是原始的地震数据$x(t)$。因此,WVD可以视为任意时刻地震解析信号的瞬时自相关的频谱。解析信号$z(t)$经过傅里叶变换之后,其频率没有负值[1],这样就避免了瞬时自相关时可能出现的正频率与负频率分量之间的相互串扰。但是,瞬时自相关即为“二次幂”,WVD为非线性变换方法,其结果不可避免地存在着不同信号分量之间的相互串扰。
平滑伪WVD方法可以有效压制串扰[9]。另有一种相对复杂的WVD办法,即以匹配追踪为基础的WVD方法。通过应用匹配追踪技术[7],把地震数据中的信号分量逐一分解出来,再用(5)式WVD算法逐一生成时间-频率谱,就可以避免不同信号分量之间的串扰。
图3展示了利用不同WVD方法计算得到的时间-频率谱。图3a图2a是同一模拟地震信号。图3b显示的是常规WVD生成的时间-频率谱,不同地震子波之间的相互串扰在时间-频率谱中形成了大量假象谱。图3c则是采用SP-WVD方法生成的时间-频率谱,有效滤除了不同信号分量之间的串扰[22]图3d中MP-WVD方法计算得到的时间-频率谱分辨率最高,其时间跨度和频率跨度与子波的主频和特征宽度有严格的对应关系。Wang给出了其时间-频率谱中的时间长度和频率宽度的展布表达式[7]。但是,应用匹配追踪技术不断将待分析信号与字典中的子波进行匹配,最终获取最佳的稀疏表征,所采用的是贪婪算法,计算效率极低。
图3 不同WVD方法计算得到的时间-频率谱
综合图2图3的分析,MP-WVD方法所得到的时间-频率谱分辨率最高,可是计算效率极低。W变换作为线性变换方法的代表之一,其计算效率远高于MP-WVD等非线性变换方法。但是,W变换时间-频率谱的分辨率没有MP-WVD时间-频率谱分辨率高,常规W变换方法需要进一步改进。

2 W变换的改进方法

由于常规线性时频分析方法的变换运算基于正交时频平面进行,因此应用范围仅限于能量在谐波空间具有较高聚焦性的信号。针对强非稳态地震数据瞬时频率变化较快的特点,常规时频分析方法无法实现较高的能量聚焦度。近年来,笔者团队针对强非稳态地震数据开展了基于W变换的各种改进方法研究,其中包括谐变时窗W变换、分数阶W变换以及线性正则W变换。与常规W变换拉伸或者压缩时窗核函数不同,该类W变换通过改变时窗核函数的角度或者时频平面角度提高时间-频率谱的分辨率与能量聚焦,进而提高地震数据对古河道、薄层等目标的提取能力。

2.1 谐变时窗W变换

为了非稳态地震数据时频分析能够达到较高的能量聚焦,笔者团队提出在常规W变换中引入旋转角度,即谐变时窗W变换。其数学表达为[19]
$\begin{align} & W{{T}_{\text{C}}}\left[ \tau,f,c\left( \tau \right) \right]=\frac{\left| {{f}_{0}}\left( \tau \right)+\Delta f \right|}{\sqrt{2\text{ }\!\!\pi\!\!\text{ }}k}\int_{-\infty }^{\infty }{x\left( t \right)w\left( t-\tau,f;\tau \right)}\times \\ & \ \ {{\text{e}}^{\text{i }\!\!\pi\!\!\text{ }c\left( t-\tau \right){{\left( t-\tau \right)}^{2}}}}{{\text{e}}^{-\text{i}2\text{ }\!\!\pi\!\!\text{ }ft}}\text{d}t \end{align}$
其中时窗函数被修正为调谐时窗:
${{w}_{\text{c}}}(t,f;\tau )=w(t,f;\tau ){{\text{e}}^{\text{i }\!\!\pi\!\!\text{ }c\left( t \right){{t}^{2}}}}$
调谐率c(t)的引入可以改变时频分析时窗的角度,进而匹配地震数据瞬时频率。选取最佳调谐率,首先将时频平面分解成($N+1$)个区域,角度取值为:
$\theta =\left\{ -\frac{\text{ }\!\!\pi\!\!\text{ }}{2}+\frac{\text{ }\!\!\pi\!\!\text{ }}{N+1},-\frac{\text{ }\!\!\pi\!\!\text{ }}{2}+\frac{\text{2 }\!\!\pi\!\!\text{ }}{N+1},\cdots,-\frac{\text{ }\!\!\pi\!\!\text{ }}{2}+\frac{N\text{ }\!\!\pi\!\!\text{ }}{N+1} \right\}$
调谐率与旋转角度之间的对应关系为:
$c=-\frac{{{F}_{\text{s}}}}{2{{T}_{\text{s}}}}\tan \theta $
对应的谐变时窗W变换则可以表示为:
$\begin{align} & W{{T}_{\text{C}}}\left( \tau,f,\theta \right)=\frac{\left| {{f}_{0}}\left( \tau \right)+\Delta f \right|}{\sqrt{2\text{ }\!\!\pi\!\!\text{ }}}\int_{-\infty }^{\infty }{x\left( t \right)w\left( t-\tau,f \right)}\times \\ & \ \ {{\text{e}}^{-\text{i }\!\!\pi\!\!\text{ }\left( {{{F}_{\text{s}}}}/{2{{T}_{\text{s}}}}\; \right){{\left( t-\tau \right)}^{2}}\tan \theta }}{{\text{e}}^{-\text{i}2\text{ }\!\!\pi\!\!\text{ }ft}}\text{d}t \end{align}$
由于该过程引入N个连续小波变换(Continuous Wavelet Transform,CW变换),时频计算的效率必须提高。另外,该过程不可避免地引入了大量人工噪声,有可能降低时频分析的分辨率。笔者团队提出谐变率可以表示为瞬时频率对时间的一阶导数,这样就消除了多个CW变换引起的人工噪声,同时提高了CW变换的效率[19]
为了提高W变换的计算效率,以满足大规模地震数据处理的需求,笔者团队借鉴离散傅里叶变换矩阵的方法,将W变换表征为矩阵与向量的乘积形式:
${{W}_{\text{c}}}\left[ :,\text{ }:,\text{ }c\left( t \right) \right]=KX$
矩阵K可以在图形处理器(GPU)上进行计算,大大节省了计算时间。

2.2 分数阶W变换

上述谐变时窗W变换将W变换分析时窗在正交时频平面中进行旋转,并通过匹配地震数据瞬时频率获取最佳分辨率以及能量聚焦度。但是,对于频率成分相同、且时间间隔较短的混叠信号,该方法效果仍然不够理想,无法将不同地震信号分量进行有效分离,因而达不到时频分析所需要的理想分辨率。
分数阶傅里叶变换作为傅里叶变换的广义化方法,通过旋转时频平面,可以更有效表征时变非平稳信号,如非稳态地震数据。为了提高W变换对于非稳态地震数据的能量聚焦度与分辨率,笔者团队将分数阶傅里叶变换与W变换相结合,提出了分数阶W变换(FrWT)方法[20]。信号$x(t)$$\alpha $阶分数阶W变换可以表示为:
$WT_{\text{Fr,}\alpha }^{{}}(\tau,u)=\int_{-\infty }^{\infty }{x(t)w\left( \tau -t,u \right)}{{K}_{\alpha }}(t,u)\text{d}t$
其中,${{K}_{\alpha }}\left( t,u \right)$$\alpha $阶分数阶傅里叶变换算子,可以表示为[1]
${{K}_{\alpha }}\left( t,u \right)=\left\{ \begin{matrix} {{A}_{\alpha }}{{\text{e}}^{\text{i }\!\!\pi\!\!\text{ }\left( {{u}^{2}}\cot \alpha -2ut\csc \alpha +{{t}^{2}}\cot \alpha \right)}} \\ \delta \left( u-t \right)\text{ } \\ \delta \left( u+t \right)\text{ } \\\end{matrix} \right.\begin{matrix} \alpha \ne k\text{ }\!\!\pi\!\!\text{ } \\ \alpha =2k\text{ }\!\!\pi\!\!\text{ } \\ \alpha =\left( 2k+1 \right)\text{ }\!\!\pi\!\!\text{ } \\\end{matrix}$
${{A}_{\alpha }}=\sqrt{{\left( 1-\text{i}\cot \alpha \right)}/{2\text{ }\!\!\pi\!\!\text{ }}\;}$
$w\left( t,u \right)=\frac{{{\left| {{f}_{\text{inst}}}+\Delta u \right|}^{p}}}{\sqrt{2\text{ }\!\!\pi\!\!\text{ }}q}\exp \left[ -\frac{{{t}^{2}}{{\left( {{f}_{\text{inst}}}+\Delta u \right)}^{2p}}}{2{{q}^{2}}} \right]$
其中 $\Delta u=\left| u\csc \alpha -{{f}_{\text{inst}}} \right|$
从(12)式可以看出,FrWT可以视为谐变时窗W变换,时频平面旋转角度为$\theta =\left( {\alpha -\text{ }\!\!\pi\!\!\text{ }}/{2}\; \right)$。当$\alpha ={\text{ }\!\!\pi\!\!\text{ }}/{2}\;$以及$p=q=1$时,FrWT可以蜕化为常规的W变换。
FrWT在W变换的基础上多引入了3个辅助参数αpq,因此,FrWT具有更高的灵活度。辅助参数的优化可以分为以下几个步骤:①确定最佳旋转角度$\alpha $,与谐变时窗W变换相同,当且仅当最佳旋转角度与待分析信号调谐率一致时,能量聚焦性最好;②针对不同$p$$q$,求取信号的分数阶W时间-频率谱;③针对不同阶数分数阶的W时间-频率谱,求取能量集中度:
$CM\left( p,q \right)=\frac{1}{\sqrt{\int_{-\infty }^{\infty }{\int_{-\infty }^{\infty }{{{\left| WT_{\text{Fr,}\ \alpha,\text{ }(p,q)}^{{}}\left( \tau,u \right) \right|}^{2}}\text{d}\tau \text{d}u}}}}$
④定义最优$p$$q$为:${{\left( p,q \right)}_{\text{opt}}}=\underset{\left( p,q \right)}{\mathop{\max }}\,\left| CM\left( p,q \right) \right|$
由于采用了(11)式的计算方法,尽管FrWT需要进行多次W变换,但时间成本较低,可以用于大规模地震数据的时频分析。分数阶W变换兼具常规W变换与分数阶傅里叶变换的特点,可以在保证低频区域高分辨率不变的同时,通过改变分数阶W变换的阶数,提高非稳态地震数据时间-频率谱的分辨率。

2.3 线性正则W变换

将上述时间-频率域中的旋转关系进一步广义化,即将W变换与线性正则变换相结合,笔者团队提出了线性正则W变换方法[21]。线性正则变换可以表征比旋转关系更为广义的变换,是仿射变换的一种[23]。仿射变换常用于图像处理、光学传播等相关应用。由于引入了3个自由参数,线性正则变换也具有比分数阶傅里叶变换更加广义的表征。
线性正则可以数学表示为:
${{X}_{\text{M}}}\left( u \right)={{L}_{\text{M}}}\left[ x(t) \right]\left( u \right)=\left\{ \begin{matrix} \int_{-\infty }^{\infty }{x(t)}{{K}_{\text{M}}}\left( t,u \right)\text{d}t\text{ }B\ne 0 \\ \sqrt{D}{{\text{e}}^{\text{i}{CD{{u}^{2}}}/{2}\;}}x\left( Du \right)\text{ }B=0 \\\end{matrix} \right.$
其中,${{K}_{\text{M}}}\left( t,u \right)$为线性正则变换核,可以数学表达为:
${{K}_{\text{M}}}\left( t,u \right)=\frac{1}{\sqrt{2\text{i }\!\!\pi\!\!\text{ }B}}{{\text{e}}^{\text{i}{\left( A{{t}^{2}}-2ut+D{{u}^{2}} \right)}/{\left( 2B \right)}\;}}$
由(17)式和(18)式可以看出,相较于常规傅里叶变换,线性正则变换添加了4个参数$A,\text{ }B,\text{ }C$$D$,这4个参数可以组成矩阵:
$M=\left[ \begin{matrix} A & B \\ C & D \\\end{matrix} \right]$
该矩阵满足辛条件,即$AD-BC=1$。该条件是为了确保信号在时频面的总支撑区域不会改变,且线性正则变换具备完美的重构能力。由(17)式和(18)式可以看出,实际计算过程中,仅当B=0时需考虑参数C的存在。而B=0时,线性正则变换退化为简单的调谐变换。
通过调整矩阵的参数,线性正则变换可以蜕化为包括菲涅尔变换、傅里叶变换以及分数阶傅里叶变换在内的诸多变换,因此,线性正则变换具有比分数阶傅里叶变换更高的自由度,应用范围也更加广泛。
将线性正则变换与W变换相结合,线性正则W变换可表示为:
$\begin{aligned}W T_{\mathrm{LC}, \mathrm{M}}(\tau, u)= & \int_{-\infty}^{\infty} x(t) w(\tau-t, u) \mathrm{e}^{\mathrm{i} \pi\left(A t^{2}-D \tau^{2}\right) /(2 B)} \mathrm{e}^{-\mathrm{i} 2 \pi(u / B) t} \mathrm{~d} t= \\& \int_{-\infty}^{\infty} x(t) w_{\mathrm{M}}(\tau-t, u) \mathrm{e}^{-\mathrm{i} 2 \pi(u / B) t} \mathrm{~d} t\end{aligned}$
其中:
$w\left( t,u,\tau \right)=\frac{{{\left| {{f}_{0}}\left( \tau \right)+\Delta u \right|}^{p}}}{\sqrt{2\text{ }\!\!\pi\!\!\text{ }}q}{{\operatorname{e}}^{-{{{t}^{2}}{{\left[ {{f}_{0}}\left( \tau \right)+\Delta u \right]}^{2p}}}/{\left( 2{{q}^{2}} \right)}\;}}$
$\Delta u=\left| {u}/{B}\;-{{f}_{0}}\left( \tau \right) \right|$
线性正则W变换可以对常规时频算法的支撑域进行包括平移、伸缩以及剪切在内的各种线性变换。由(21)式可以看出,当矩阵M取值为$M=[\cos \alpha,$ $\sin \alpha,\text{ }-\sin \alpha,\text{ }\cos \alpha ]$时,线性正则W变换蜕化为$\alpha $阶的分数阶W变换;当$M=[\text{ }0]$时,线性正则W变换蜕化为常规W变换。不同W变换之间的支撑域如图4所示。可以看出,从常规W变换到分数阶W变换,再到线性正则W变换,是不断泛化的过程,提高了W变换的自由度。
图4 不同W变换时窗在时间-频率平面中的显示图
线性正则W变换的参数寻优与分数阶W变换大致相同,区别在于在线性正则变换中,当且仅当${A}/{B}\;$的值与信号调谐率的相反数相同时,时间-频率谱的能量聚焦达到最高。线性正则W变换的提出,不仅在非稳态地震数据中薄层油气藏识别中具有较高的应用前景,而且还开拓了W变换在光学、图像处理学中的相关应用。
图5对比了不同W变换方法计算得到的时间-频率谱。所有W变换方法都在低于主频的区域得到了时间方向的很好的分辨率,而改进的W方法尤其在主频为低频的子波分离方面取得了明显效果。
图5 不同W变换方法计算得到的时间-频率谱

3 W变换应用实例

本文展示3个W变换改进方法的实际应用效果。首先,提取实际地震数据,验证W变换在河道砂体地震检测中的应用效果。第1个实例的地震沿层切片上(见图6a)可以识别大量古河道,但是河道边界不清晰,连续性也需改进。图6b为沿层切片时间信息,图6c图6d为对应的SP-WVD时间-频率谱沿层切片和常规W变换得到的时间-频率谱沿层切片,频率为40 Hz。W变换的时间-频率谱切片显示较为清晰的河道边缘,完整地刻画了河道走向。因此,W变换时间-频率谱对于地下储层边界具有很好的刻画效果。
图6 常规W变换与SP-WVD方法进行河道识别实例
第2个实例验证W变换在溶洞检测中的应用效果。图7a是从实际三维地震数据体中抽取的地震剖面,该段地震剖面具有明显的“串珠状”地震反射。但是,时间域地震剖面不能很好地表征“串珠状”地震反射是否由溶洞影响造成。本文利用SP-WVD方法、MP-WVD方法和常规W变换方法对该剖面进行时频分析,并抽取40 Hz对应的频率切片。
图7 SP-WVD变换、MP-WVD变换与常规W变换溶洞检测效果对比实例
上述3种方法可以有效将原始地震剖面分解为不同频率的地震频谱切片,有利于地震数据的溶洞检测识别。但是,与常规基于正交时窗变换的时频分析方法不同,W变换在低频中也具有较高的分辨率,有助于地震数据“低频阴影”的探明,进而指导钻探过程。
第3个实例验证不同W变换方法在多期河道检测中的应用效果。图8a为三维地震数据体的等时切片,可以看出该区域发育多条古河道,但是在地震数据等时剖面上分辨率较差,河道边界不清晰,河道展布不明确。图8b图8e则是常规W变换、调谐时窗W变换、分数阶W变换和线性正则W变换计算得到的时间-频率谱RGB融合对比图。RGB三色分别对应20,30,40 Hz的单频数据。
图8 常规W变换、调谐时窗W变换、分数阶W变换与线性正则W变换河道检测对比实例
从RGB融合切片图上可以看出,该区域发育有多条河道,且在融合切片上比地震数据等时切片上的边界清晰。分数阶W变换和线性正则W变换得到的时间-频率谱中河道更为清晰,且干扰较少,有效提高了地震古河道探测的分辨率。

4 结论

对W变换时频分析方法及其改进方法进行了综述性介绍。近年提出的W变换是典型的线性时频分析方法。线性时频分析通过地震信号与特定窗函数的褶积而实现,理论清晰,运算简单,但是通常无法同时达到时间域和频率域上的高分辨率和能量聚焦,尤其在低频区域其分辨率较差。针对线性时频分析方法在低频区域分辨率较差的问题,笔者近年提出了W变换,引入地震数据瞬时频率作为参数之一,构建关于地震数据瞬时频率对称的分析时窗,提高了线性时频分析方法在低频区域的分辨率。
将W变换与WVD方法进行对标。WVD方法是非线性时频分析方法的代表,因其非线性特性导致的地震数据不同信号分量之间的串扰问题,严重降低了时频分析的效果。但是,WVD方法展示的时频域中的能量分布,明确指示地震子波的时间重心和频率重心,对W变换具有直接的对标意义。W变换中引入了任意时间位置所对应的瞬时频率作为变换参数,因此所展示的时间-频率谱也具有明确的能量聚焦重心。
详细阐述了W变换的不同改进方法,其目的是对W变换方法不断广义化,利用包括尺度变换、旋转、剪切等在内的线性变换,进一步提高W变换时间-频率谱的分辨率与能量聚焦。这些W变换方法具有对于非稳态地震数据的适用性,尤其是对于富含油气藏区域,通过对时频轴进行不同程度的旋转变换可以提高时间-频率谱的能量聚焦度。这些方法在油藏描述、古河道识别等方面均取得了良好的实际应用效果。
符号注释:
ABCD——线性正则变换中的参变量,无因次;$c(t)$——信号谐变率,Hz/s;${{A}_{\alpha }}$——分数阶傅里叶变换归一化系数,无因次;CM——时间-频率谱能量聚焦度,无因次; f——频率,Hz;${{f}_{0}}\left( t \right)$——时间t处地震数据瞬时频率,Hz;${{f}_{0}}\left( \tau \right)$——分析时窗位置地震数据瞬时频率,Hz;${{f}_{\text{inst}}}$——瞬时频率,Hz;Fs——频率采样间隔,Hz;i——虚数单位;k——时频分析时窗尺度因子,无因次;K——不同待分析频率下离散调谐W变换矩阵;${{L}_{\text{M}}}$——线性正则变换算子;${{K}_{\alpha }}\left( t,u \right)$——α阶W变换算子;${{K}_{\text{M}}}\left( t,u \right)$——线性正则变换核函数;M——线性正则变换中的参变量ABCD组成的矩阵;N——广义线性调谐变换中时间-频率平面划分个数;pq——控制时窗长度的辅助参变量,无因次;${{S}_{\text{SP-WVD}}}\left( t,f \right)$——SP-WVD时间-频率谱,无因次;${{S}_{\text{WVD}}}\left( t,f \right)$——WVD时间-频率谱,无因次;$ST\left( t,f \right)$——S变换时间-频率谱,无因次;t——时间,s;Ts——时间采样间隔,s;u——辅助参变量,无因次;$w\left( t \right)$——W变换的时窗核函数;$w\left( t,f;\tau \right)$——W变换时窗函数;${{w}_{\text{c}}}\left( t,f;\tau \right)$——调谐时窗W变换时窗函数;${{w}_{\text{M}}}\left( \tau,u \right)$——线性正则W变换时窗函数;${{W}_{\mathbf{c}}}\left[ :,:,c\left( t \right) \right]$——调谐W变换时间-频率谱矩阵;$WT\left( t,f \right)$——W变换时间-频率谱,无因次;$W{{T}_{\text{C}}}\left[ \tau,f,c\left( \tau \right) \right]$——调谐时窗W变换时间-频率谱,无因次;$WT_{\text{Fr,}\ \alpha }^{{}}\left( \tau,u \right)$——分数阶W变换时间-频率谱,无因次;$WT_{\text{LC,}\ \text{M}}^{{}}\left( \tau,u \right)$——线性正则W变换时间-频率谱,无因次;$c\left( t \right)$$x\left( t \right)$——地震信号;X——待分析地震信号矩阵;${{X}_{\text{M}}}$——线性正则变换频谱,dB;$z\left( t \right)$——地震信号解析形式;α——分数阶W变换阶数;τ——分析时窗位置,s;π——圆周率,3.14;*——褶积;ϕ——子波相位项;θ——时频平面角度;$c(t)$δ——单位阶跃函数;Δf——地震数据主频与待分析频率之间差的绝对值,Hz;Δu——不同变换下待分析频率域地震数据主频之间差的绝对值,Hz。
[1]
WANG Y H. Time-frequency analysis of seismic signals[M]. Hoboken, NJ: John Wiley & Sons, Ltd, 2023.

[2]
PORTNOFF M R. Time-frequency representation of digital signals and systems based on short-time Fourier analysis[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1980, 28(1): 55-69.

[3]
DAUBECHIES I. The wavelet transform, time-frequency localization and signal analysis[J]. IEEE Transactions on Information Theory, 1990, 36(5): 961-1005.

[4]
STOCKWELL R G, MANSINHA L, LOWE R P. Localization of the complex spectrum: The S transform[J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001.

[5]
WANG Y H. The W transform[J]. Geophysics, 2021, 86(1): V31-V39.

[6]
WIGNER E. On the quantum correction for thermodynamic equilibrium[J]. Physical Review, 1932, 40(5): 749-759.

[7]
WANG Y H. Seismic time-frequency spectral decomposition by matching pursuit[J]. Geophysics, 2007, 72(1): V13-V20.

[8]
OBERLIN T, MEIGNEN S, PERRIER V. The Fourier-based synchrosqueezing transform[R]. Florence:2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2014: 315-319.

[9]
CHOI H I, WILLIAMS W J. Improved time-frequency representation of multicomponent signals using exponential kernels[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1989, 37(6): 862-871.

[10]
WANG Y H, RAO Y, XU D. Multichannel maximum-entropy method for the Wigner-Ville distribution[J]. Geophysics, 2020, 85(1): V25-V31.

[11]
RAO Y, GUO Y X, XU D. Detecting karst voids based on dominant frequencies of seismic profiles[J]. Pure and Applied Geophysics, 2021, 178(8): 3057-3067.

[12]
ALSALMI H, WANG Y H. Mask filtering to the Wigner-Ville distribution[J]. Geophysics, 2021, 86(6): V489-V496.

[13]
MALLAT S G, ZHANG Z F. Matching pursuits with time-frequency dictionaries[J]. IEEE Transactions on Signal Processing, 1993, 41(12): 3397-3415.

[14]
WANG Y H. Multichannel matching pursuit for seismic trace decomposition[J]. Geophysics, 2010, 75(4): V61-V66.

[15]
ZHAO Z C, RAO Y, WANG Y H. Structure-adapted multichannel matching pursuit for seismic trace decomposition[J]. Pure and Applied Geophysics, 2023, 180(3): 851-861.

[16]
LI R, ZHU X, ZHOU Y Z, et al. Generalized W transform and its application in gas-bearing reservoir characterization[J]. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 1-5.

[17]
CHEN H, PENG L, CHEN X P, et al. A three-parameter W transform and its application to gas reservoir identification[J]. Geophysics, 2022, 87(5): V521-V532.

[18]
LUO C P, ZONG Z Y. The synchroextracting algorithm based on W transform and its application in channel characterization[J]. IEEE Geoscience and Remote Sensing Letters, 2023, 20: 1-5.

[19]
ZHAO Z C, RAO Y, WANG Y H. The W transform with a chirp-modulated window[J]. Geophysics, 2023, 88(2): V93-V100.

[20]
ZHAO Z C, RAO Y, WANG Y H. The fractional-order W-transform[J]. IEEE Geoscience and Remote Sensing Letters, 2024, 21: 1-5.

[21]
ZHAO Z C, RAO Y, WANG Y H. The linear canonical W transform[J]. Geophysics, 2024, 89: in press.

[22]
JEONG J, WILLIAMS W J. Alias-free generalized discrete-time time-frequency distributions[J]. IEEE Transactions on Signal Processing, 1992, 40(11): 2757-2765.

[23]
MOSHINSKY M, QUESNE C. Linear canonical transformations and their unitary representations[J]. Journal of Mathematical Physics, 1971, 12(8): 1772-1780.

文章导航

/