刘笛， 胡英*， 陈辉， 李军， 方玉霞
1 地球勘探与信息技术教育部重点实验室(成都理工大学)， 成都 610059 2 地质灾害防治与地质环境保护国家重点实验室(成都理工大学)， 成都 610059 3 数学地质四川省重点实验室(成都理工大学)， 成都 610059 4 成都理工大学计算机与网络安全学院(牛津布鲁克斯学院)， 成都 610059
瞬时谱分解技术对于描述非平稳信号的时变特征具有重要的实用价值(Marfurt and Kirlin，2001；
潘晓等，2020；
潘辉等，2021；
田琳和胡津健，2021)，其实质是利用数学手段对地震信号进行时频分析(徐天吉和程冰洁，2008)，以提取地震信号的特征参数，进而直接检测油气藏的存在，是地球物理学的研究热点之一(Huang et al.，2016；
江雨濛等，2021).用于地质勘探的传统时频分析方法主要包括短时傅里叶变换(Short-Time Fourier Transform，STFT)(邓攻等，2015)、小波变换(Continue Wavelet Transform，CWT)(Sinha et al.，2009；
李伟等，2017)、S变换(S Transform，ST)(Wu and Castagna，2017)、保幅S变换(Amplitude Preserving S-Transform，APST)(Wang，2016；
Wang and Lu，2018)以及广义S变换(General S Transform，GST)(陈学华等，2009，2011；
Fang et al.，2021)等，其中，由于小波变换多尺度多分辨率的特性，使其在地震信号处理中有着区别于短时傅里叶变换和S变换的独特优势，当小波基的波形与地震子波相似时，CWT的时频谱将具有更高的时频聚焦性以及良好的噪声鲁棒性(Tian et al.，2022).高静怀等(2006)在改进Morlet小波的基础上，提出了一种三参数小波变换(Three-Parameter Wavelet Transform，TPWT)，即使在中心频率较小时也可以满足小波的容许性条件，通过其中三个参数的调节，可以灵活匹配各种地震信号.但因为受到海森堡不确定性原理的限制，致使传统时频分析方法不能在时间和频率上以任意精度定位信号(Auger et al.，2013；
韩利等，2016；
Yuan et al.，2019)，难以准确地描述非线性平稳信号的非线性特征(Nikoo et al.，2016；
严海滔等，2019).随着油气藏勘探领域向深层复杂储层转变，迫切需要发展高精度的时频分析方法用于描述复杂多变的地震信号(乐友喜等，2018).
Daubechies等(2011)基于小波变换提出了一种基于相位的时频分析技术，即同步挤压变换算法(Synchrosqueezing Transform，SST)，该算法将小波变换后的结果在时频域内进行挤压和重排，使得信号的瞬时频率能够更加接近真实的频率，提高了信号时频分析结果的时频分辨率；
受到同步挤压变换算法的启发，Yu等(2017)提出了基于STFT的同步提取变换(SET)，与同步挤压变换不同的是其仅提取与信号瞬时频率高度相关的时频信息，进而产生比SST能量更聚焦的时频结果，并且提高了噪声鲁棒性(Li et al.，2020a,b).有学者(高静怀等，2018；
Wang et al.，2020；
Chen et al.，2021；
曾爱平等，2022)将此类后处理方法应用在地震资料处理中，进行非常规油气藏的烃类检测和储层识别，取得了良好的应用效果(石战战等，2018).
SET虽然已经成为提高时频结果分辨率的有效时频分析工具，然而受到STFT时域窗口或频域窗口固定的限制，在分析变化较快或较慢的信号时，效果依旧不理想(Shi et al.，2021).因此，选择合适的时频原子至关重要(Tian et al.，2022).诸多学者在这一领域进行了许多尝试，如时间方向同步提取(Li et al.，2020a,b)、同步提取S变换(冯永鑫等，2020)及同步提取广义S变换(Hu et al.，2020)等，在同步提取变换的基础上进一步提升时频分析的效果，其中S变换及广义S变换时频谱的频率偏高(Li et al.，2022；
Chen et al.，2022；
Wang et al.，2021)，易造成时频表征假象，保幅保频S变换的提出(Wang et al.，2021)为S变换的改进奠定新的思路，改善了频率上移的现象.但是S变换类方法时窗固定不易匹配不同类型地震信号(Wang，2021)，同时还存在低频段处时间分辨较差的问题(Li et al.，2022).为了保证后续进行同步提取变换的时频分辨率和准确性，具有多尺度多分辨率特性的三参数小波变换根据频率-尺度进行伸缩和平移可以得到比STFT、ST和GST等方法更稳定的时频分辨率，并且三个可调参数使其有效匹配包含快变分量的信号，提高了匹配地震信号的灵活度.因此，本文聚焦致密砂岩气藏，提出了基于同步提取三参数小波变换(SETPWT)的高分辨率后处理谱分解方法，该方法结合三参数小波灵活多变的优势，通过同步提取算法，进一步提高了三参数小波变换的分辨能力.在接下来的章节中，首先阐述了SETPWT的方法原理，并用合成数据体现其能量聚焦性和抗噪性，然后通过实际地震资料处理的结果，说明该方法在致密砂岩储层含气性检测方面的有效性和实用性.
1.1 三参数小波变换回顾
对于任意待分析信号s(t)∈L2(R)，t为时间，R为实数集合，L2表示可积函数空间，则三参数小波变换定义为：
(1)
ψ(t;Γ)=e-τ(t-β)2{p(Γ)×[cos(σt)-k(Γ)]
+iq(Γ)sin(σt)}，
(2)
其中：
(3)
(4)
(5)
式中，Γ=(σ,τ,β)，σ为基小波的调制频率，控制公式中三角函数的频率，影响小波的震荡程度，σ越大，小波震荡越剧烈；
τ作为能量衰减因子，控制衰减函数的衰减速度，从而控制窗口大小；
β为能量延迟因子.其中β对基小波形态的影响比较复杂，当β取值为基小波周期的整数倍时，小波仅在时间轴上发生时移，否则，β还会使小波产生形变.β的取值与待分析信号的相位有关，β=0时，可匹配零相位子波；
β≠0时，可匹配非零相位子波.在地震数据实际处理过程中，β非零虽然可以更好的匹配非零地震子波，取得更好的时频结果，但是由于地层结构复杂，β取值使基小波发生时移所产生的地质体埋深的解释差异尚未进行讨论，因此，目前学者们使用三参数小波进行地震数据处理时，一般设定β=0.
不同参数下的小波形态如图1所示.固定参数τ的值为1，调节σ的值，观察小波基形态变化，从图1a中可以看到，σ越大，小波震荡程度越大；
当固定σ的值不变，调节τ的值，如图1b所示，τ越大，小波衰减速度越快，窗口越窄.
1.2 同步提取三参数小波变换
1.2.1 基于单频信号的SETPWT理论
首先考虑一个单频解析信号s(t)=A(t)e2πiφ(t)，φ(t)表示该信号的瞬时相位，其理想时频(Ideal Time-Frequency Analysis，ITFA)定义为：
ITFA=A(t)·δ(ω-φ′(t))，
(6)
(7)
(8)
利用狄拉克函数，取与上述瞬时频率估计值高度相关的时频信息，将同步提取算子SEO写作如下形式：
(9)
(10)
将公式(7)代入(10)可得：
(11)
同步提取三参数小波变换含3个可调参数Γ=(σ,τ,β)，其中σ≥0,τ>0，且β≥0.不同的参数组合将产生不同形态的小波基，给同步提取三参数小波变换的结果带来不同的计算精度.
1.2.2 基于多频信号的SETPWT理论及重构
公式(6)—(11)的推导基于单频信号s(t)，但是在实际应用处理中，待分析信号往往是多频率成分的.为了增强算法的实用性，给出如下多频信号的同步提取三参数小波变换推导过程.
对于一个多频信号ms(t)，其公式为：
(12)
其中，msk(t)=Ak(t)ei2πφk(t)，Ak(t)和φk(t)是第k个频率分量的瞬时振幅和相位函数，满足|φ′k(t)-φ′m(t)|>2Δ，∀m=k-1，k,m∈{1,2,…,n}，其中φ′k(t)和φ′m(t)是瞬时频率，Δ表示频率紧支集.
假设，|A′k(t)|≤ε且φ′k(t)≤ε′，其中ε和ε′都是很小的值，那么一个单频信号msk(t)可以被认为是一个纯谐波信号.因此，根据公式(7)，msk(t)的三参数小波变换结果为：
(13)
此时，使第k个频率分量对时间求偏导，可以得到：
(14)
信号ms(t)的三参数小波变换结果为：
(15)
所以，根据公式(8)，结合公式(14)和(15)，该多频信号的瞬时频率估计为：
(16)
(17)
则公式(16)可以近似为：
(18)
因此，多频信号msk(t)的同步提取三参数小波变换(SETPWT)的表达式为：
(19)
由于每个频率分量信号的小波变换系数分布在不同时间-尺度区域中各自的中心频率轨迹周围，因此SETPWT对信号ms(t)的时频结果可以表示为：
(20)
公式(20)表示信号ms(t)的时频结果可以被看做是每个频率分量的所有时间-尺度子区域的叠加.由于同步提取变换类方法仅提取与瞬时频率高度相关的信息，其余信息完全被剔除，因此同步提取三参数小波变换只能近似重构，并不能精确重构.信号ms(t)近似重构公式为：
(21)
作为一种小波变换的后处理技术，SETPWT具有明显的多分辨率特性，且由于是对小波变换结果进行了提取处理，因此其时频谱能量相对于传统时频分析方法更聚焦.
2.1 单频信号测试与分析
首先本节将同步提取三参数小波变换应用到以下单频信号中，该信号表达式为：
f(t)=(1-0.2cos(πt))·cos(250πt-10sin(8πt)
(22)
其中，t∈[0,1]，采样频率为2048 Hz，其波形如图2a所示，随着时间增加，信号振幅增强.分别使用CWT、TPWT、SET和SETPWT对该合成信号进行处理，分别得到如图2b—e的时频谱，并将其中t∈[0.45,0.55]，频率[45,100]范围内的信息进行局部放大，并将细节图排列在右侧列对应位置.
图2 (a) 原始信号图像；
分别由 (b) CWT；
(c) TPWT；
(d) SET；
(e) SETPWT 计算产生的频谱图Fig.2 (a) Waveform of the synthetic signal；
The spectrum calculated by (b) CWT；
(c) TPWT；
(d) SET；
(e) SETPWT
从时频结果可以看出，CWT和TPWT可以较准确表征信号的时频轨迹及能量变化，并且表现出低频处频率分辨率高，高频处时间分辨率高的特点，但是时频能量在脊线附近相对分散.SET和SETPWT都有相较于CWT结果更高的可读性，充分体现了SEO的提取能力.然而，SET由于窗函数固定，时频分辨率受到限制，在信号频率调制部分时频脊线附近的能量有些许发散，且出现了错误的能量团提取，如图2c及其细节图所示；
SETPWT通过调节参数，最终当Γ=(3,1,0)时取得较好的时频结果，利用同步提取因子将时频脊线附近的能量去除，得到一条清晰且能量聚焦性强的时频能量脊线(如图2d所示)，细节图中的局部细节相较于另外的方法也更加精确，未出现错误能量的提取现象.
表1 不同方法处理该信号的任意熵Table 1 Rényi entropy calculated by different methods of the signal
2.2 多频信号测试与分析
将本文的方法用于如下加噪变幅调频信号，其表达式为：
(23)
其中，t∈[0,1]，采样频率为2048 Hz.n(t)表示信噪比(Signal-to-Noise Ratio，SNR)为8 dB的白高斯噪声，模型x1(t)的瞬时频率变化缓慢；
模型x2(t)在t∈[0,0.6]范围内具有很强的调频特性和广泛的瞬时频率范围，t∈[0.6,1]时为一个频率为400 Hz的谐波；
模型x3(t)为频率为100 Hz的谐波，以扩展多分量信号在频率域上的范围.
无噪声信号模型X(t)的示意图如3a所示，分别使用CWT、TPWT、SET和SETPWT进行时频分析得到图3b—e，为了便于观察，如图3b—e的矩形框所示，提取时频谱中时间[0,0.2]，频率[270,570]范围内的局部信息放大，进行对比分析，对应显示在图3b—e的右侧.从整体来看，四种方法都可以表征信号的每个分量，但是CWT时频谱的能量聚焦程度较低，并且在t∈[0,0.2]范围内，未能将模型x1(t)和模型x2(t)清晰地分离开，模型x1(t)和模型x3(t)彼此之间也相互影响；
TPWT通过调节其参数，获得的时频曲线相对CWT更聚焦，并将其中模型x1(t)和模型x3(t)完整分离；
SET结果的能量聚焦性大幅提高，但是在矩形框内信号强频率调制位置能量提取产生有误能量团.相比之下，SETPWT的能量提取最为准确，多个分量之间完整分离，且能量脊线聚焦，可读性高，对于频率急剧变化的位置也可以很好的进行瞬时频率估计，体现其对非平稳信号的适用性.
图3 (a) 未加噪信号X(t)的原始信号的图像；
分别由 (b) CWT；
(c) TPWT；
(d) SET；
(e) SETPWT 计算产生的频谱图Fig.3 (a) Waveform of the synthetic noiseless signal X(t)；
The spectrum calculated by (b) CWT；
(c) TPWT；
(d) SET；
(e) SETPWT
为进一步验证同步提取三参数小波对噪声的鲁棒性，在原始信号基础上加SNR=8 dB的白高斯噪声.加入噪声之后的信号模型Y(t)如图4a所示，CWT受到噪声影响较大，TPWT相对于CWT具有较好的噪声鲁棒性(如图4b、c)；
对比观察图4d、e，SET的结果明显受到噪声影响，矩形块内瞬时频率曲线上的能量严重发散.而SETPWT的时频谱即使在噪声影响的情况下，仍然可以表示出清晰的时频曲线，并且与下方的chirp信号分离开来，信号的瞬时特征得到了很好的刻画，而且时频谱的分辨率在时间和频率域上均呈现了最佳的状态.
图4 (a) 加噪信号Y(t)的原始信号图像；
分别由(b) CWT；
(c) TPWT；
(d) SET；
(e) SETPWT计算产生的频谱图Fig.4 (a) Waveform of the synthetic signal Y(t) with noise；
The spectrum calculated by (b) CWT；
(c) TPWT；
(d) SET；
(e) SETPWT
为了表征SETPWT处理该信号在不同SNR噪声下的性能，图5中分别计算了具有不同SNR(0 dB、7 dB、14 dB、21 dB、28 dB)的加噪信号的Rényi熵.可以看出SETPWT时频结果的Rényi熵一直处于最低的状态，也就是说，SETPWT可以为加噪信号提供最好的时频浓度.
图5 不同SNR加噪信号的Rényi熵Fig.5 The Rényi entropy of adding-noise signal with different SNRs
通过SETPWT在该加噪强调制信号的应用，表现该方法有效对抗噪声影响的优越性，同时由于三参数小波变换可调节参数丰富，可以适应多种类型信号的特点，使得信号模型中频率突变的地方也得到了精细刻画，综合表现该方法的抗噪性以及对非平稳突变信号刻画的稳定性.
当地震信号穿过含气储层时，流体流动会导致地震信号能量和频率的异常衰减，该现象主要表现为高频能量相对损失和低频能量相对增强(Castagna et al.，2003；
Yin et al.，2015).通过分析不同地质体的频率地震响应特征，利用谱分解技术识别油气藏.因此，采用中国四川盆地西部的中江气田二维地震数据进行全频带瞬时谱分解，研究该方法在地震信号时频处理中的适用性.该实际资料包含286道地震记录，每道记录有551个时间采样点，采样间隔为2 ms.图6a为该实际资料的地震剖面图，从左至右两条竖线分别表示井A和井B.其中，已知井A为发育良好的含气井，井B是一口干井，分别由矩形框和椭圆框标注.
首先提取井A和井B单道地震信号，分别如图6b、d所示，并对其进行三参数小波变换分析，分别如图6c、e所示，与图6a对应观察，将地震反射层信息用黑色矩形框在图6b—e中标定显示.从井A和井B的单道地震信号的时频结果图中可以看出，两道地震信号在反射层处均有明显的能量增强现象.然后提取地震反射层的频率信息，如图7所示，井A振幅在28 Hz处达到峰值，然后随频率的增加而减小，而井B振幅在此频率段的变化为在28 Hz附近先增强后减弱，后在40 Hz处再增强，且高频处幅值高于低频处(如图7中黑色矩形框所示).根据两个过井道地震信号的频率信息，分别选取28 Hz和40 Hz作为低频剖面和高频剖面，用于识别致密砂岩气藏.
图6 (a) 中江气田原始地震剖面图；
(b) 井A单道地震信号谱；
(c) 井A时频谱；
(d) 井B单道地震信号谱；
(e) 井B时频谱Fig.6 (a) Seismic section from the Zhongjiang Gas Field；
(b) Seismic signal of well A；
(c) The TPWT spectrum of well A；
(d) Seismic signal of well B；
(e) The TPWT spectrum of well B
图7 过井A和井B单道数据振幅谱Fig.7 The amplitude spectrum of well A and B
为了验证SETPWT相对于CWT、TPWT和SET的优异性，分别采用这四种时频分析方法对实际资料进行处理，根据上述单道地震频谱信息得出的结论，从处理结果中提取28 Hz和40 Hz两个常频率剖面，用于储层含气性检测，结果如图8所示.并将矩形框处的结果进行放大，便于细节观察，如图9所示.从图8可以看出，四种方法共有的特点是低频剖面(如图8a、c、e、g)上与井A相交的矩形框中能量较高，而高频剖面(如图8b、d、f、h)上的相应位置有明显的能量衰减现象，根据地震波穿过含气储层时，流体流动会导致地震信号低频能量增强，高频能量衰减这一特性可以判定，四种方法均可验证井A是一口含气井，井B是一口干井.由此可见，四种方法均能有效对致密砂岩含气性进行识别，但在描述致密砂岩含气性的精度上存在一定差异.众所周知，精确圈定储层对识别储层具有重要意义，因此进一步分析低频剖面中矩形框内致密砂岩含气性特征的准确性.CWT大致能够表现地震信号的异常响应(如图9a、b)，但是难以对致密砂岩气藏产生准确响应；
与CWT相比，TPWT(如图9c、d)和SET(如图9e、f)更清晰的描述了该气藏的低频特征，但是能量相对发散.经过SEO对能量谱进行提取后，SETPWT的处理结果(如图9e、f)相较于CWT、TPWT和SET的刻画精度显著提高，时频分辨率和能量聚焦性进一步提升.对比四种方法的高频剖面(如图9b、d、f、h)，SETPWT的能量衰减表现得更加强烈，致密砂岩含气储层产生的低频强能量，高频弱能量特征也更加明显.
图8 由(a)(b) CWT；
(c)(d) TPWT；
(e)(f) SET；
(g)(h) SETPWT 计算产生的低频(28 Hz)剖面和高频(40 Hz)剖面的对比显示Fig.8 Dominant-frequency (28 Hz) and high-frequency (40 Hz) sections generated by (a)(b) CWT；
(c)(d) TPWT；
(e)(f) SET；
(g)(h) SETPWT
图9 由 (a)(b) CWT；
(c)(d) TPWT；
(e)(f) SET；
(g)(h) SETPWT 计算产生的低频(28 Hz)剖面和高频(40 Hz)剖面中矩形框局部细节图对比显示Fig.9 Local detail sections in dominant-frequency (28 Hz) and high-frequency (40 Hz) sections generated by (a)(b) CWT；
(c)(d) TPWT；
(e)(f) SET；
(g)(h) SETPWT
总之，SETPWT能够有效提取地震信号经过含气储层的异常响应特征，高、低频剖面对比证实了SETPWT在致密砂岩储层含气性识别的优势.
本文提出了同步提取三参数小波变换，并应用于中江地区致密砂岩含气性检测中，该变换利用地震波经过含气储层会产生“低频能量强，高频能量弱”这一特性，可以检测油气藏的存在.通过与CWT、TPWT以及SET对比发现，本文方法具有更高的时频分辨率，能够进行准确精细的频率脊线提取，对于强调制变幅信号，不仅可以准确揭示待分析信号的幅度信息，并且在加噪的情况下也依然可以反映原始信号的能量曲线，体现其抗噪性以及对非平稳信号的适应性.此外，实际地震资料算例表明，本文方法能够提供能量更集中的时频表示，证明其在非常规油气藏含气性识别领域中具有一定潜力.
猜你喜欢时频小波分辨率构造Daubechies小波的一些注记科技风(2021年19期)2021-09-07基于MATLAB的小波降噪研究电子制作(2019年13期)2020-01-14EM算法的参数分辨率数学物理学报(2019年3期)2019-07-23原生VS最大那些混淆视听的“分辨率”概念家庭影院技术(2018年9期)2018-11-02基于改进的G-SVS LMS 与冗余提升小波的滚动轴承故障诊断制造技术与机床(2017年10期)2017-11-28基于深度特征学习的图像超分辨率重建自动化学报(2017年5期)2017-05-14一种改进的基于边缘加强超分辨率算法成都信息工程大学学报(2017年6期)2017-03-16基于时频分析的逆合成孔径雷达成像技术舰船科学技术(2015年8期)2015-02-27对采样数据序列进行时频分解法的改进电测与仪表(2014年17期)2014-04-04基于FPGA小波变换核的设计电测与仪表(2014年8期)2014-04-04
