基于正则化理论的时频分析方法及应用

doi:10.11720/wtyht.2023.1369

基于正则化理论的时频分析方法及应用

张金强^,¹^,²^,³

1.中国石化页岩油气勘探开发重点实验室,北京 102206

2.中国石油化工股份有限公司石油勘探开发研究院,北京 102206

3.页岩油气富集机理与有效开发国家重点实验室,北京 100026

A regularization theory-based method for time-frequency analysis and its applications

ZHANG Jin-Qiang^,¹^,²^,³

1. Key Laboratory of Shale Oil and Gas Exploration & Production,Sinopec,Beijing 102206,China

2. Sinopec Petroleum Exploration and Production Research Institute,Beijing 102206,China

3. State Key Laboratory of Shale Oil and Gas Enrichment Mechanisms and Effective Development,Beijing 100026,China

第一作者: 张金强(1970-),男,博士,高级工程师,1991年获中国石油大学(华东)勘查地球物理专业学士学位,1994年、1998年分别获中国石油大学(北京)应用地球物理专业硕士学位、应用地球物理专业博士学位;现就职于中国石化股份公司石油勘探开发研究院,主要从事岩石物理、地震储层预测等科研工作。Email: zhangjq.syky@sinopec.com

责任编辑: 叶佩

收稿日期: 2022-07-22 修回日期: 2023-02-1

基金资助:

中石化科技攻关项目“常压页岩气地球物理评价技术研究”(P21087-3)
“麻黄山西区延长—延安组有效开发关键技术”(P22184)

Received: 2022-07-22 Revised: 2023-02-1

摘要

时频分析方法在地震勘探中有广泛的应用,因而获得具有良好时频分辨率的时频分析算法至关重要。传统的时频分析方法存在着一定的局限性,为克服这些局限性,提出了基于正则化理论的时频分析方法。该方法认为,短时窗信号是不同频率谐波的叠加,应从求解反问题的角度考察时频分析问题。在此视角下,时频分析问题具有不适定性,为得到有意义的时频谱,需要在正则化理论框架下进行时频分析。考察了正则化理论中常用的L₁范数约束、L₂范数约束以及最小支撑约束条件下的求解方法,并将3种约束函数的求解方法统一到同一个求解框架中。通过数值分析表明,最小支撑约束的时频分析方法具有较高的时频分辨率。将方法系统应用于一个特定研究区的实际资料,获得了具有较高时频分辨率的时频数据体,并利用单频数据体清晰刻画了储层的平面展布范围,展示了方法良好的应用前景。

关键词： 时频分析; 正则化理论; L₁范数约束; L₂范数约束; 最小支撑约束; 时频谱

Abstract

Time-frequency analysis (TFA) has been widely used in seismic exploration,thus it is crucial to develop a TFA algorithm with high time-frequency resolution.Given the limitations of conventional TFA methods,this study proposed a TFA method based on the regularization theory.The proposed method considers the signal in a short-time window as a superposition of harmonics with different frequencies and takes the TFA problem as an inverse problem.From this perspective,the TFA problem is ill-posed and needs to be solved based on the regularization theory to get a significant time-frequency spectrum.The solution methods under the conditions of L₁ and L₂ norm constraints and the minimum support constraint are commonly used in the regularization theory.This study investigated these solution methods and unified them into the same solution framework.Numerical analysis shows that the TFA method under the condition of the minimum support constraint yielded high time-frequency resolution.This method was systematically applied to the actual data of a specific study area,producing a time-frequency data volume with high time-frequency resolution.Moreover,the planar reservoir distribution was clearly characterized using a single-frequency data volume,demonstrating the promising application prospect of the method.

Keywords： time-frequency analysis; regularization theory; L₁ norm constraint; L₂ norm constraint; minimium support constraint; time-frequency spectrum

PDF (6166KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

张金强. 基于正则化理论的时频分析方法及应用[J]. 物探与化探, 2023, 47(4): 965-974 doi:10.11720/wtyht.2023.1369

ZHANG Jin-Qiang. A regularization theory-based method for time-frequency analysis and its applications[J]. Geophysical and Geochemical Exploration, 2023, 47(4): 965-974 doi:10.11720/wtyht.2023.1369

0 引言

时频分析将一维时间信号变换到时频域,可以突出时间信号在不同时间局部的频域特征,在地震资料处理与地震储层预测中获得了广泛的应用。在地震资料处理领域,将地震信号变换到时频域有利于区分噪声和信号,因此,时频分析方法被广泛地应用于地震噪声压制^{[1⇓⇓⇓-5]}。在地震储层预测中,由于不同频率的地震信号对不同厚度储层的响应不同,时频分析方法被广泛地应用于河道、三角洲等特殊储层预测中^[6⇓⇓-9]。

地球物理领域常用的时频分析方法主要有短时傅里叶变换、连续小波变换、S变换及广义S变换、WVD、Hilbert-Huang变换以及匹配追踪等,这些方法各有其特点和局限。短时傅里叶变换、连续小波变换以及S变换都是线性时频分析方法^[10-11],都需要对时间信号加窗截断得到局部信号,从而获得不同时间局部的频域特征^[12-13]。这类方法尽管时频聚集特性不同(小波变换及S变换较短视窗傅里叶变换具有更好的时频聚集特性),但是都会受到海森伯格测不准原理的约束。WVD(Wigner-Ville 分布)具有较高的时频分辨率,但是要受到交叉项的影响^[14-15],导致其应用局限。Hilbert-Huang变换(HHT)是基于经验模态分解的方法,具有较高的时频聚焦特性,但是也存在模态混叠等问题^[16]。匹配追踪方法(MP)通过建立完备的时频信号字典,利用反复迭代寻求以最少时频原子获得信号的最佳匹配^[17-18]。这种方法由于采用稀疏表示技术,通常可以获得比连续小波变换和Gabor变换更好的时频分辨率,但是由于其本身遍历性的迭代算法,计算成本较高。

稀疏时频分析方法是当前较受关注的一类时频分析方法。Puryear等^[19]提出了约束最小二乘时频分析方法,获得了较高的时频分辨率。李海英等^[20]基于正则化方法提出了最小二乘约束下的时频分析方法,二者都未对不同正则化约束函数的求解方法进行充分讨论。田琳等^[21]讨论了基于L₁范数约束的稀疏短时傅里叶变换时频分析方法。本文在前人研究的基础上,结合有关反演理论,基于正则化理论框架提出一种更具普遍性的稀疏时频分析方法,并对不同正则化约束函数条件下的求解方法进行了讨论,得到了具有较广适用范围的算法流程。

1 方法原理

在传统的离散傅里叶变换理论中,有限离散信号与其离散谱之间存在以下关系:

(1)

X_{m} = \overset{N - 1}{\sum_{n = 0}} x_{n} e^{- i 2 π n Δ t m Δ f}

(

m = 0,1, \dots, N - 1

(2)

x_{n} = \frac{1}{N} \overset{N - 1}{\sum_{m = 0}} X_{m} e^{i 2 π n Δ t m Δ f}

(

n = 0,1, \dots, N - 1

式中: $x_{n}$ 为以Δ $t$ 为采样间隔,在[0,(N-1)Δ $t$ ]区间内取值的有限离散信号; $X_{m}$ 为以Δ $f$ 为基频的有限离散信号的离散谱。Δ $t$ 、Δ $f$ 之间满足:

(3)Δf=

\frac{1}{N Δ t}

。

式(1)和式(2)可表达为矩阵形式:

X=Ax,
x=

\frac{1}{N}

A^*X,

式中:X= ${(X_{0}, X_{1}, \dots, X_{N - 1})}^{T}$ ;x= ${(x_{0}, x_{1}, \dots, x_{N - 1})}^{T}$ ;A为矩阵;A^*为矩阵A的共轭转置矩阵,A矩阵的n行m列元素为:a_nm=e^-i2π^nm/N。在上述理论体系中,A为对称正定矩阵,组成A矩阵的行向量或列向量是相互正交的。按式(1)定义的地震信号的离散谱包含了整个地震道频谱特征,在应用中要识别和凸显信号频率随时间而变化的特征就要进行时频分解。一个自然延伸的想法是对地震信号用短时窗截断,研究短时窗内地震信号的频谱特征,连续滑动时窗进而得到整个地震道的时频谱。这是短时窗傅里叶变换的基本思想,但是这种方法得到的频谱等于信号频谱与窗函数频谱的褶积,导致信号频谱的频域分辨率降低。为提高短时窗信号的频谱分辨率,以反演的思路考察这一问题。

设离散短时窗信号为d,根据傅里叶变换的基本理论,d可以由不同频率的谐波叠加而得到,参照式(2):

(4)d=Fm,

式中: $d = {(d_{0}, d_{1}, \dots, d_{L - 1})}^{T}$ 是代表离散短时窗信号的数据向量,其维度为 $L \times 1$ ;F为核矩阵,其第n行m列的元素为: $F_{n m} = e^{i 2 π n Δ t m Δ f}$ ,按此定义的核矩阵其第n行所对应的行向量为:

F_{n} = (1, e^{i 2 π n Δ t Δ f}, e^{i 4 π n Δ t Δ f}, \dots, e^{i 2 (M - 1) π n Δ t Δ f}),

所表征的为同一时刻不同频率谐波的函数值。核矩阵的第m列对应的列向量为:

F_{m} = {(1, e^{i 2 π Δ t m Δ f}, e^{i 4 π Δ t m Δ f}, \dots, e^{i 2 π (L - 1) Δ t m Δ f})}^{T},

所表征的为同一频率谐波不同时刻的函数值。F矩阵的维度为L×M,L为短时窗信号的样点数,M为不同频率谐波的个数,根据具体问题确定有意义的频带范围和频率采样间隔,最终确定M。m= ${(m_{0}, m_{1}, \dots, m_{M - 1})}^{T}$ 对应于不同频率谐波的加权系数向量,即离散信号的频谱。与式(2)不同的是Δt、Δf之间不再满足式(3),因此矩阵F的行向量和列向量不再是相互正交的。在上述定义下,问题可以描述为,已知离散短时窗信号d和核矩阵F,设法获得m的一个估计 $\overset{︿}{m}$ ,使得F $\overset{︿}{m}$ ≈d。这一问题通常是不适定的,需要在正则化理论框架下求解,根据正则化理论,上述问题可以转化为求解m,使得

(5)

E (m) = {(F m - d)}^{2} + α S (m)

取得最小值。其中

{(F m - d)}^{2} = {(F m - d)}^{*} (F m - d),

表征的是利用估计频谱得到离散信号与实测信号的偏差。S(m)为约束函数,α为决定约束函数贡献大小的平衡参数,α越大约束函数的贡献越大。选择不同的约束函数,得到的解也不同。常用的约束函数有

(6)

S_{L_{1}} (m) = \overset{M - 1}{\sum_{i = 0}} |m_{i}|

(7)

S_{L_{2}} (m) = \overset{M - 1}{\sum_{i = 0}} {|m_{i}|}^{2}

(8)

S_{β} (m) = \overset{M - 1}{\sum_{i = 0}} \frac{{|m_{i}|}^{2}}{{|m_{i}|}^{2} + β^{2}}

式中: $m_{i}$ 为m的第i个元素;β为一个任意小的数。式(6)、(7)所表示的约束函数在反演中经常使用,Zhang等^[22]认为L₁范数约束较之L₂范数约束的反演结果具有更好的稀疏性;而Portniaguine等^[23-24]研究证明以式(8)作为约束函数可以得到最佳聚焦的反演图像。

尽管约束函数有多种形式,但多数约束函数可以从形式上表示为某个函数的L₂范数形式,即 $S (m)$ 可以表示为 $S (m) = (f (m), f (m))$ 。对于式(6)有

(9)

S_{L_{1}} (m) = \overset{M - 1}{\sum_{i = 0}} |m_{i}| = \overset{M - 1}{\sum_{i = 0}} {({|m_{i}|}^{1 / 2})}^{2}

。

对于式(8)有

(10)

S_{β} (m) = \overset{M - 1}{\sum_{i = 0}} \frac{{|m_{i}|}^{2}}{{|m_{i}|}^{2} + β^{2}} = \overset{M - 1}{\sum_{i = 0}} {[\frac{|m_{i}|}{{({|m_{i}|}^{2} + β^{2})}^{\frac{1}{2}}}]}^{2}

。

按照以上表述,更进一步,可以将约束函数表示为

(11)

S (m) = (f (m), f (m)) = (W_{e} (m) m, W_{e} (m) m)

对照式(11)和式(9)有

(12)

W_{e} (m) = d i a g ({|m_{i}|}^{- 1 / 2})

对照式(11)和式(10)有

(13)

W_{e} (m) = d i a g [{({|m_{i}|}^{2} + β^{2})}^{- \frac{1}{2}}]

这样,式(5)可以表示为

(14)

E (m) = {(F m - d)}^{2} + α (W_{e} (m) m, W_{e} (m) m)

式(14)用矩阵的形式表示为

(15)

E (m) = {(F m - d)}^{2} + α {(W_{e} m)}^{2}

对应于式(6)、(7)、(8)的约束函数的矩阵形式分别为

(16)

W_{e} = d i a g ({|m_{i}|}^{- \frac{1}{2}})

(17)

W_{e} = I

(18)

W_{e} = d i a g [{({|m_{i}|}^{2} + β^{2})}^{- \frac{1}{2}}]

式(15)可以进一步变换为如下形式:

(19)

E (m) = {(F {W_{e}}^{- 1} W_{e} m - d)}^{2} + α {(W_{e} m)}^{2}

令:

F_{w} = F {W_{e}}^{- 1}

m_{w}

W_{e} m

最终将式(5)变换为如下形式:

(20)

E (m) = {(F_{w} m_{w} - d)}^{2} + α {(m_{w})}^{2}

式(20)与经典的带正则化项的最小二乘问题形式上完全一致,所不同是F_w是与m有关的矩阵,因此需要进行迭代求解。以使用式(8)作为约束函数的情况说明迭代过程。

首先,注意到采用常规L₂范数(即式7)作为约束函数时,W_e=I,此时可以用常规最小二乘得到问题的解:

(21)m=

{(F^{*} F + α I)}^{- 1}

F^*d 。

可以用式(21)作为问题的初始解 $m_{0}$ ,在以后的迭代过程中,首先利用式(17)得到 $W_{e}$ ,并利用式(19)计算 $F_{w}$ 与 $m_{w 0}$ ,然后利用下式求解m_w及m:

(22)m_w=

{(F_{w}^{*} F_{w} + α I)}^{- 1} F_{w}^{*}

(23)m=

{W_{e}}^{- 1}

m_w 。

利用更新的m后,进入到下一个迭代过程,重复计算W_e,直至再次更新m。上述循环的终止原则可以选择使式(5)值趋于稳定,也可以根据经验选择合适的迭代次数。

上述过程得到一个短时窗信号的高分辨率频谱,实际应用中,我们应用一个窗函数对地震道进行逐点滑动,这样就可以得到整个地震道的时频谱。窗函数可以选择汉宁窗函数,由于汉宁窗函数具有由中心点向两端迅速衰减的特征,所以我们估算的其实是中心点附近数据的频谱。算法的这种特征,要求将时窗中心点作为时间零点,因此核矩阵应按此进行调整。

对于地震道时频分析,可以先将地震道进行希尔伯特变换,得到对应的复地震道,在此基础上进行相应的时频分析。

2 数值算例

设计一个由两个单频信号叠加而成的信号: $x (t) = s i n (2 π f_{1} t) + s i n (2 π f_{2} t)$ ,取 $f_{1}$ 、 $f_{2}$ 分别等于20 Hz和50 Hz,两个单频信号和合成后的信号如图1所示。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 两个单频信号的合成信号

a—20 Hz正弦信号;b—50 Hz正弦信号;c—图1a、b中两个信号的合成信号

Fig.1 Synthetic signal of the two sinuasoidal signals

a—sinuasoidal signal of 20 Hz;b—sinuasoidal signal of 50 Hz;c—summary signal of fig.1a and fig.1b

分别应用式(6)、(7)、(8)作为约束函数,按照本文所提出的方法求取信号的时频谱。为方便标记,将应用式(6)、(7)、(8)作为约束函数的算法标记为RL1TFS、RL2TFS、RMSTFS,同时应用连续小波变换求取信号的时频谱(CWT)作为对比和参照。图2展示了几种算法的对比结果。由图2可以看出,利用式(8)作为约束函数的算法具有最好的时频分辨率,非常清晰地分辨出了不断时间的两个单频信号,较之常用的连续小波变换具有更高频率分辨率。而L₂范数作为约束函数的时频谱图像光滑,时频分辨率较低。以L₁范数作为约束函数项时频谱大致与小波变换的时频谱具有相同的时频分辨能力。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 不同时频分析方法得到时频谱对比

a—CWT时频谱;b—RL2TFS时频谱;c—RL1TFS时频谱;d—RMSTFS时频谱

Fig.2 Comparisons of time-frequency spectrum of different time frequency analysis method

a—time-frequency spectrum of CWT;b—time-frequency spectrum of RL2TFS;c—time-frequency spectrum of RL1TFS;d—time-frequency spectrum of RMSTFS

通过实验分析可以看出,最小支撑约束稀疏谱算法具有较L₁范数约束和L₂范数约束的算法更好的时频分辨率。因此重点对此方法进行进一步研究。为进一步说明最小支撑约束稀疏谱算法对短时窗信号的频谱的分辨能力,研究不同长度时窗条件下算法估算频谱的精度,设计如下算例。雷克子波是地震勘探中经常用到的子波,该子波有理论频谱,便于研究。因此我们使用不同长度的汉宁窗对雷克子波进行截断,比较短时傅里叶变换与最小支撑约束稀疏谱算法得到的频谱。首先,我们用长度为40 ms的汉宁窗截断雷克子波,研究两种算法得到的频谱(图3),然后用长度为20 ms的汉宁窗截断雷克子波,研究其频谱(图4)。在进行短时傅里叶变换时,为得到高精度频谱,同样以汉宁窗截断信号,将短时窗信号外数据充零并延长数据长度至1 024 ms。对比图3与图4可以看出,首先最小支撑约束稀疏谱较之传统的短时傅里叶变换的频谱更为接近雷克子波的理论频谱,其次,当汉宁窗的长度为40 ms时,短时傅里叶变换的频谱形态与理论谱虽有较大差异,但在高频段仍有一定的符合度,当汉宁窗的长度变为20 ms时,短时傅里叶变换的频谱已经完全不能反映雷克子波的频谱特征了,最小支撑约束稀疏谱虽然与理论谱的差异增大,但仍可以较好地拟合雷克子波的理论频谱。另外,图3c及图4c中都可以看到利用复地震道为输入的最小支撑约束频谱(红色曲线)较以原始地震道为输入的最小支撑约束频谱(藏青色曲线)更接近于雷克子波的理论曲线。在汉宁窗长度为40 ms时,二者整体比较接近,当汉宁窗长度为20 ms时,利用复地震道作为输入的最小支撑约束频谱明显地优于利用原始地震道作为输入的结果。这主要是因为,复地震道虚部相当于对原始地震道作了90°相移,增加了数据冗余,同时复谐波的实部与虚部也呈90°相移,这样利用复地震道反演频谱增加了反演的稳定性和精度。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 40 ms汉宁窗作用于雷克子波得到的信号及其对应的频谱

a—40 ms汉宁窗函数及雷克子波叠合;b—汉宁窗作用于雷克子波得到信号;c—不同时频分析方法计算图3b信号的频谱曲线

Fig.3 40 ms Hanning function applied on the ricker wavelet and the correpsonding spectrum

a—40 ms Hanning function and Riker wavelet;b—signal obtained by applying Hanning on ricker wavelet;c—frequency spectrum of the signal (fig. 3b) obtained by different time-frequency method

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 20 ms汉宁窗作用于雷克子波得到的信号及其对应的频谱

a—20 ms汉宁窗函数及雷克子波叠合;b—汉宁窗作用于雷克子波得到信号;c—不同时频分析方法计算的图4b信号的频谱曲线

Fig.4 20 ms Hanning function applied on the ricker wavelet and the correpsonding spectrum

a—20 ms Hanning function and Riker Wavelet;b—signal obtained by applying Hanning on ricker wavelet;c—frequency spectrum of the signal(fig. 4b) obtained by different time-frequency method

为使模型更加贴近实际情况,进一步,设计如下两个算例。在第1个算例中,令地震道由若干对不同厚度的反射系数对组成,每对反射系数大小相同,符号相同。在第2个算例中,地震反射系数的时间厚度与第一个算例相同,但每对反射系数大小相同,符号相反。两个算例中,反射系数对的时间厚度都是由4 ms增加到32 ms,每个反射系数对较上一个系数对厚度增加4 ms。两个算例均采用反射系数与主频为30 Hz的雷克子波进行褶积得到合成地震道,并分别应用连续小波变换以及最小支撑约束的稀疏谱算法计算地震道的频谱。结果如图5、图6所示。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 顶底反射系数相同不同时间厚度地层的合成地震记录及不同算法计算的时频谱

a—地震反射系数序列;b—合成地震道;c—CWT时频谱;d—RMSTFS时频谱

Fig.5 Synthetic record of different thikness layers with same recflction coefficient on top and bottom and the corresponding time-frequency spectrum

a—seismic reflection series;b—synthetic record;c—spectrum of CWT;d—spectrum of RMSTFS

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 顶底反射系数相反不同时间厚度地层的合成地震记录及不同算法计算的时频谱

a—地震反射系数序列;b—合成地震道;c—CWT时频谱;d—RMSTFS时频谱

Fig.6 Synthetic record of different thikness layers with revese recflction coefficients on top and bottom and the corresponding time-frequency spectrum

a—seismic reflection series;b—synthetic record;c—spectrum of CWT;d—spectrum of RMSTFS

对比两个算例的结果,可以看出:连续小波变换的频谱在低频段的时间分辨率低,有明显的自上而下的条纹,是明显的假象。最小支撑约束的时频谱得到的结果显示时频分辨率更高,时频谱能量团聚焦清晰,主能量团时间跨度普遍小于连续小波变换时频谱上对应主能量团的时间跨度,说明方法具有较好的时间分辨率。在图5d中,反射系数对的时间厚度为12 ms时反射的顶底清晰可分,并且在此之后,反射顶底之间频谱显示出调谐作用造成的陷波点随着时间厚度的增大有规律地向低频偏移的现象。连续小波变换的频谱也表现出类似的现象,但其能量团聚焦差,与地层反射的顶底对应关系差。

对比两个算例还可以发现,当反射系数对大小相同符号相同时,时频谱的分辨率更高,其调谐作用的陷波点偏向低频段,因而也更为明显。为更清楚地解释此现象,计算了两种反射系数对在不同厚度条件下调谐作用滤波器的频率响应曲线,如图7所示。当反射系数对大小相同、符号相反时,厚度为16 ms时第一个陷波点为62.5 Hz,由于地震子波本身的频带范围限制,不能清晰观察。当厚度增加到24 ms时,第一个陷波点为41.67 Hz,在时频谱上可以清晰观察。对于大小相等、符号相反的反射系数对,当厚度为16 ms时,陷波点即在地震频带内,可以清晰观察。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 不同厚度反射系数对调谐滤波器的响应曲线

a—图5a中部分反射系数对的调谐滤波器响应曲线;b—图6a中部分反射系数对的调谐滤波器响应曲线

Fig.7 The frequency response curves of reflection pairs with different time thickness

a—the frequency response curves for reflection paris of fig. 5a;b—the frequency response curves for reflection paris of fig. 6a

将时间信号变换到时频域这个过程必然是具有多解性的,对于一些简单的算例,如上文中所设计的算例,我们可以通过理论分析的方法清楚地获知其理论时频谱的特征,因而易于识别哪些现象是地层本身的反映,哪些现象是算法的假象,而对于现实中更为复杂的时间信号,则往往难以判定,因此也更要求算法能够捕获信号中地层的响应特征。通过数据实验算例可以看出,在正则化理论框架下,利用最小支撑约束条件求解的最小支撑稀疏时频谱具有良好的时频谱特征,在时频域都具有良好的分辨能力,同时可以清楚地展现地层的调谐、陷波作用的基本特征,是一种非常有应用潜力的时频分析方法。通过实验算例分析,也启示我们:由于地层反射波的时频谱存在谐振区(振幅极大值)和陷频区(振幅极小值),因此,在时频域振幅分析时,需根据地震数据的频率范围和地层厚度变化范围合理选取时频分析的频率范围;从另一方面说,选择了某一频率得到的时频剖面,对其振幅的变化的分析和解释也应注意谐振区和陷频区的影响。

3 实际资料应用

方法在鄂尔多斯盆地西缘某工区进行了系统应用。该区以侏罗系延安组为主要目的层,储层为河流—三角洲沉积,河道砂体横向变化快,预测困难。应用Partyka等^[25]提出的做法,首先对地震道进行时频分解,然后将抽取不同频率对应的能谱值形成单频数据体。图8为利用连续小波变换和最小支撑约束稀疏时频分析得到的不同频率的能谱剖面。从这些剖面可以看出,在低频段,如10 Hz和20 Hz的能谱剖面,连续小波变换得到能谱剖面时间分辨率低,在剖面上部缺失有意义的低频信号,10 Hz的能谱剖面尤为明显,而最小支撑约束稀疏时频分析得到的能谱剖面则在剖面的不同部位都呈现出有意义的信号。在高频信号段,如60 Hz,两种方法得到的能谱剖面相似度较高,最小支撑约束稀疏时频分析得到能谱剖面的时间分辨率仍优于连续小波变换。对10 Hz的低频能谱剖面进行局部放大,如图9所示。在图9b的能谱剖面的红色椭圆内的高能谱值反映了工区内主要目的层的储层响应,而在连续小波变换的对应剖面上则没有类似的异常响应。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 连续小波变换和最小支撑约束稀疏谱得到不同频率的能谱剖面

a—连续小波变换10 Hz能谱剖面;b—RMSTFS 10 Hz能谱剖面;c—连续小波变换20 Hz能谱剖面;d—RMSTFS 20 Hz能谱剖面;e—连续小波变换60 Hz能谱剖面;f—RMSTFS 60 Hz能谱剖面

Fig.8 Sections of different frequency spectrum obtained by CWT and RMSTFS

a—section of 10 Hz spectral spectrum of CWT;b—section of 10 Hz spectral spectrum of RMSTFS;c—section of 20 Hz spectral spectrum of CWT;d—section of 20 Hz spectral spectrum of RMSTFS;e—section of 60 Hz spectral spectrum of CWT;f—section of 60 Hz spectral spectrum of RMSTFS

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 连续小波变换及最小支撑约束稀疏时频分析得到10 Hz能谱剖面对比

a—连续小波变换10 Hz能谱剖面(局部放大);b—最小支撑约束时频分析10 Hz能谱剖面(局部放大)

Fig.9 Comparisons of 10 Hz spectrum section obtained by CWT and RMSTFS

a—section of 10 Hz spectral spectrum of CWT(zoomed);b—section of 10 Hz spectral spectrum of RMSTFS(zoomed)

进一步在单频数据体上,提取目的层的沿层切片。图10为在连续小波变换和最小支撑约束时频分析得到的10 Hz单频能谱数据体上得到的沿层切片。根据钻井及地质分析认为,在图10的中部区域有一近SN向展布的河流沉积储层展布。X3、X6、X8井位于河道主体部位,具有较厚的砂体厚度,X1井位于河道边部,砂体厚度薄,X2、X4、X5及X7井均为泛滥平原沉积,砂体缺失。在最小支撑约束时频分析得到10 Hz单频数据体切片上可以清楚地凸显中部的河道沉积异常(X1井以东,X4井及X5井以西),而在对应的连续小波变换数据体切片上则难以得到这种清楚的地质认识。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 不同算法计算的能谱沿层切片

a—连续小波变换10 Hz能谱沿层切片;b—最小支撑约束时频分析10 Hz能谱沿层切片

Fig.10 Horizon slices of 10 Hz spectrum by different algorithms

a—Horizon slice of 10 Hz CWT spectrum;b—Horizon slice of 10 Hz RMSTFS spectrum

4 结论

利用求解反问题的思路考察时频分析问题,时频分析问题可以描述为正则化理论框架下的反问题。不同的正则化约束函数约束下的时频谱求解可以统一到同一个求解框架。数值计算结果表明,最小支撑约束的时频分析方法具有最佳的时频分辨率。方法在特定研究区的应用显示了良好的时频分析方法有助于凸显储层的地震响应,在地震储层预测中具有较好的应用前景。

方法通过迭代求解信号的时频谱,运算量较大。因此,算法在运算效率的提高上有待进一步研究。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

曹思远, 袁殿.

高分辨率地震资料处理技术综述

[J]. 新疆石油地质, 2016, 37(1):112-119.