几何属性是最重要的地震属性^[1-2]之一, 在地震资料解释中发挥着重要作用。基于地震数据的层位信息, 可获得地层倾角及方位角^[3-4]等几何属性, 估算曲率^[5-9]和相干属性^[10-13]等, 用于地震资料解释。在倾角估计方面, PICOU等^[14]提出了非归一化互相关扫描法; BARNES^[15-16]提出了复数分析法; MARFURT等^[17]发展了相干扫描方法; BAKKER等^[18]提出了梯度结构张量法; LUO等^[19]提出了加权结构法。为减少地震数据振幅横向变化对倾角估算精度的影响, WANG等^[20]深入分析了复偏导数与瞬时振幅及瞬时相位的偏导数间的关系, 提出了基于瞬时相位的梯度结构张量法, 并基于多窗口分析技术提出了基于瞬时相位的多窗口梯度结构张量法^[21]。利用时频分析的多尺度分辨能力对实际地震资料进行高精度时频分析获得高精度的时频瞬时相位谱, 抽取不同频率的多尺度高精度时频瞬时相位谱估算地层倾角和曲率属性可以提高该方法的应用效果。

地层倾角及曲率属性的估算精度依赖于时频瞬时相位谱的精度, 而时频瞬时相位谱的精度受时频分析方法及随机噪声的影响。短时傅里叶变换^[22]、S变换^[23-24]和广义S变换^[25-27]可获得信号的瞬时振幅及初始相位随时间和频率的分布, 即时频瞬时振幅谱和时频初始相位谱; 小波变换^[28-33]可获得信号的瞬时振幅及瞬时相位^[34-38]随时间和尺度(或频率)的分布, 即时频瞬时振幅谱和时频瞬时相位谱。改进的短时傅里叶变换^[39]通过对待分析的实信号解析信号加上对称时窗函数并进行傅里叶变换, 获得待分析实信号的瞬时振幅、初始相位和瞬时相位随时间和频率的分布(分别称为高精度时频瞬时振幅谱、时频初始相位谱和时频瞬时相位谱)。本文利用改进的短时傅里叶变换获得高精度时频瞬时相位谱, 然后抽取不同频率的高精度时频瞬时相位谱估算多尺度曲率属性, 并将该属性应用于实际三维地震数据的不连续性检测。

1 改进的短时傅里叶变换 1.1 基本原理

基于短时傅里叶变换、广义S变换和小波变换, 实信号x(t)改进的短时傅里叶变换^[39]可表示为:

$ \begin{array}{l} H\left( {\tau , f} \right) = \int_{ - \infty }^\infty {z\left( {n\tau + t} \right)p\left[ {\left( {n - 1} \right)\tau + t} \right] \cdot } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}ft} \right){\rm{d}}t \end{array} $

(1)

式中:τ是时间平移因子, f为频率, n={0, 1}, z(t)为利用希尔伯特变换得到的实信号x(t)的解析信号, ${\rm{i}} = \sqrt { - 1} $, p(t)为与频率有关的实函数, 且

$ \left\{ \begin{array}{l} \int_{ - \infty }^\infty {p\left( t \right){\rm{d}}t = 1} \\ p\left( { - t} \right) = p\left( t \right) \end{array} \right. $

(2)

方程(1)可写为:

$ \left\{ \begin{array}{l} {H_0}\left( {\tau , f} \right) = \int_{ - \infty }^\infty {z\left( t \right)p\left( {t - \tau } \right)\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}ft} \right){\rm{d}}t} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;当n = 0时\\ {H_1}\left( {\tau , f} \right) = \int_{ - \infty }^\infty {z\left( {\tau + t} \right)p\left( t \right)\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}ft} \right){\rm{d}}t} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;当n = 1时 \end{array} \right. $

(3)

式中, H₀(τ, f)反映复信号z(t)p(t－τ)中频率成分f的瞬时振幅和初始相位, H₁(τ, f)反映复信号z(t)p(t－τ)中频率成分f在t=τ时刻的瞬时振幅和瞬时相位。若定义

$ \left\{ \begin{array}{l} {H_0}\left( {\tau , f} \right) = {A_0}\left( {\tau , f} \right)\exp\left[ {{\rm{i}}{\theta _0}\left( {\tau , f} \right)} \right]\\ {H_1}\left( {\tau , f} \right) = {A_1}\left( {\tau , f} \right)\exp\left[ {{\rm{i}}{\theta _1}\left( {\tau , f} \right)} \right] \end{array} \right. $

(4)

且

$ \left\{ \begin{array}{l} {A_0}(\tau , {f_\tau }) = \max \left[ {{A_0}\left( {\tau , f} \right)} \right]\\ {A_1}(\tau , {f_\tau }) = \max \left[ {{A_1}\left( {\tau , f} \right)} \right]\\ {\theta _0}(\tau , {f_\tau }) = {\theta _0}(\tau , f = {f_\tau })\\ {\theta _1}(\tau , {f_\tau }) = {\theta _1}(\tau , f = {f_\tau }) \end{array} \right. $

(5)

则称A₀(τ, f), θ₀(τ, f), A₀(τ, f_τ), θ₀(τ, f_τ)和f_τ分别为实信号x(t)在t=τ时刻的高精度时频瞬时振幅谱、时频初始相位谱、时频峰值瞬时振幅谱、时频峰值初始相位谱和时频峰值频率; A₁(τ, f), θ₁(τ, f), A₁(τ, f_τ)和θ₁(τ, f_τ)分别为实信号x(t)在t=τ时刻的高精度时频瞬时振幅谱、时频瞬时相位谱、时频峰值瞬时振幅谱和时频峰值瞬时相位谱。由于

$ {H_0}\left( {\tau , f} \right) = {H_1}\left( {\tau , f} \right)\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}f\tau } \right) $

(6)

因此, 有

$ \left\{ \begin{array}{l} {\theta _0}\left( {\tau , f} \right) = {\theta _1}\left( {\tau , f} \right) - 2{\rm{ \mathsf{ π} }}f\tau \\ {A_0}\left( {\tau , f} \right) = {A_1}\left( {\tau , f} \right) \end{array} \right. $

(7)

1.2 合成记录分析

对改进的短时傅里叶变换、短时傅里叶变换及小波变换进行了时频分析对比, 如图 1、图 2所示。这里假设短时傅里叶变换可表示为:

$ S\left( {\tau , f} \right) = \int_{ - \infty }^\infty {x\left( t \right)p\left( {t - \tau } \right)\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}ft} \right){\rm{d}}t} $

(8)

且

$ p\left( t \right) = \frac{{\lambda {{\left| f \right|}^\beta }}}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}\exp \left( { - \frac{{{\lambda ^2}{{\left| f \right|}^{2\beta }}}}{2}{t^2}} \right) $

(9)

式中, f为频率, λ>0且β≥0。假设小波变换可表示为:

$ W\left( {\tau , a} \right) = \frac{1}{{\sqrt a }}\int_{ - \infty }^\infty {x(t){q^*}\left( {\frac{{t - \tau }}{a}} \right)} {\rm{d}}t $

(10)

式中, a为尺度算子, q(t)为小波函数, *表示取共轭, 且

$ q(t) = p(t)\exp\left( {{\rm{i}}2{\rm{ \mathsf{ π} }}ft} \right) $

(11)

方程(10)的小波变换可表示为:

$ \begin{array}{l} W\left( {\tau , a = 1} \right) = \int_{ - \infty }^\infty {x\left( {t + \tau } \right)p\left( t \right)\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}ft} \right){\rm{d}}t} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ = }}\exp \left( {{\rm{i}}2{\rm{ \mathsf{ π} }}f\tau } \right)S\left( {\tau , f} \right) \end{array} $

(12)

图 1a为合成记录, 图 1b和图 1c分别为改进的短时傅里叶变换得到的高精度时频瞬时振幅谱A₁(τ, f)和A₀(τ, f), 图 1d和图 1e分别为小波变换和短时傅里叶变换得到的时频瞬时振幅谱。对比图 1b~图 1e可知, 各种时频瞬时振幅谱的变化趋势完全相同, 改进的短时傅里叶变换得到的高精度时频瞬时振幅谱变化范围为0~3, 与合成记录的瞬时振幅变化区间相同, 而小波变换和短时傅里叶变换得到的时频瞬时振幅谱在0~1.5之间变化。图 1f和图 1g分别为改进的短时傅里叶变换得到的高精度时频瞬时振幅谱与小波变换得到的时频瞬时振幅谱的差值与比值。由图 1f可见, 两者的差值全部大于0;由图 1g可见, 当频率从4Hz变化到26Hz时, 比值系数为2.0, 较大比值主要分布在低频和高频端。

图 2a和图 2d分别为改进的短时傅里叶变换得到的高精度时频瞬时相位谱和高精度时频初始相位谱, 两者变化趋势相同; 图 2b和图 2e分别为小波变换和短时傅里叶变换得到的时频瞬时相位谱和时频初始相位谱, 两者变化趋势也相同。图 2c为改进的短时傅里叶变换的高精度时频瞬时相位谱与小波变换的时频瞬时相位谱之差值, 图 2f为改进的短时傅里叶变换的高精度时频初始相位谱与短时傅里叶变换的时频初始相位谱之差值, 差值的变化区间为[-45°, 45°]。当频率从4Hz变化到26Hz时, 差值几乎为零, 较大的差值主要分布在低频和高频端。

综合分析可知, 在相同时窗函数情况下, 改进的短时傅里叶变换具有更好的时频精度。

2 基于高精度时频瞬时相位谱的多尺度曲率属性 2.1 基于瞬时相位的曲率属性

曲率属性估算的核心是利用地震数据准确估算地层倾角, 而利用地震信号的振幅信息估算地层倾角是常用的方法。由于振幅不仅反映地层倾角信息, 还反映波阻抗横向变化, 因此基于振幅估算的曲率属性受地层倾角及波阻抗横向变化的综合影响, 在波阻抗横向变化剧烈的地区应用效果不佳。瞬时相位仅反映地层倾角信息, 与波阻抗横向变化无关, 因此可基于地震信号的瞬时相位估算地层倾角, 最终获得较高精度的曲率属性。

以下通过二维合成地震记录及相关属性进行说明。图 3a为无噪合成地震记录, 共三套层位:第一套(0.05s处)和第三套(0.15s处)均为水平层; 第二套层位(0.10s附近)共有7个起伏(中心分别在30, 90, 150, 210, 270, 330和390道处), 且起伏程度随道号增大而增加, 另外, 当道号取值为440~470时, 第二套层位仅有波阻抗变化而无倾角变化。对图 3a合成地震记录数据进行希尔伯特变换, 得到的瞬时相位如图 3b所示:当道号为440~470时, 第二套层位对应的地震信号瞬时相位无任何变化。图 3c为利用图 3a所示的地震记录数据计算得到的地层倾角:当道号为440~470时, 第二套层位的倾角信息受地震信号振幅的影响严重, 与理论模型不匹配。图 3d为利用图 3b所示的瞬时相位计算得到的地层倾角, 与理论模型的倾角信息完全匹配。

综上所述, 基于瞬时相位估算的地层倾角信息可以避免地震信号振幅的空间变化带来的假象。基于时频分析中低频时频谱具有较高频率分辨率、高频时频谱具有较高时间分辨率的多尺度思想, 可利用改进的短时傅里叶变换方法对地震资料进行时频谱分析, 获得高精度的时频瞬时相位谱, 并抽取高精度的时频瞬时相位谱进行曲率属性计算, 以提高曲率属性的估算精度。

2.2 基于高精度时频瞬时相位谱的多尺度曲率属性

为了减小随机噪声对时频瞬时相位谱的影响, 首先利用改进的短时傅里叶变换对地震信号进行时频分析, 并利用高精度的时频峰值瞬时振幅谱和高精度的时频峰值瞬时相位谱对输入的地震信号进行重构; 然后利用改进的短时傅里叶变换对重构信号进行时频分析得到高精度的时频瞬时相位谱; 最后抽取共频率的高精度时频瞬时相位谱进行多尺度曲率属性计算。

基于方程(4)和方程(5), 输入实信号x(t)的重构信号可表示为:

$ \begin{array}{l} r\left( \tau \right) = {\mathop{\rm Re}\nolimits} \left[ {{H_1}\left( {\tau , {f_\tau }} \right)} \right] = {A_1}\left( {\tau , {f_\tau }} \right) \cdot \\ \;\;\;\;\;\;\;\;\;\;\cos \left[ {{\theta _1}\left( {\tau , {f_\tau }} \right)} \right] \end{array} $

(13)

图 4验证了高精度时频峰值瞬时振幅谱和高精度时频峰值瞬时相位谱重构地震信号的可行性。图 4a为20Hz余弦信号; 图 4b和图 4c分别为该余弦信号的瞬时振幅和瞬时相位; 图 4d为图 4a所示余弦信号加入强随机噪声的结果; 图 4e和图 4f分别为图 4d所示20Hz含噪余弦信号的瞬时振幅和瞬时相位, 其均被随机噪声污染。

利用改进的短时傅里叶变换对图 4d所示含噪信号进行时频分析, 得到高精度时频瞬时振幅谱和高精度时频瞬时相位谱, 并基于高精度时频峰值瞬时振幅谱及高精度时频峰值瞬时相位谱进行信号重构。图 5a为图 4a所示无噪信号的高精度时频瞬时振幅谱; 图 5b为图 4a所示无噪信号的高精度时频瞬时相位谱; 图 5c为图 4d所示含噪信号的高精度时频瞬时振幅谱, 其受随机噪声影响较小; 图 5d为图 4d所示含噪信号的高精度时频瞬时相位谱, 其受随机噪声干扰严重。图 6a为图 5c中各时刻的高精度时频峰值瞬时振幅谱, 图 6b为图 5d中各时刻的高精度时频峰值瞬时相位谱; 图 6c为利用图 6a和图 6b重构的信号, 其与图 4a所示的20Hz余弦信号完全一致, 说明可基于高精度时频峰值瞬时振幅谱和高精度时频峰值瞬时相位谱重构信号。

3 应用实例

利用某工区三维地震资料检验了基于高精度时频瞬时相位谱的多尺度曲率属性估算方法的应用效果。处理流程如下:首先利用方程(9)所示的时窗函数对原始地震数据进行改进的短时傅里叶变换, 得到高精度时频瞬时振幅谱和高精度时频瞬时相位谱; 然后利用高精度时频峰值瞬时振幅谱及高精度时频峰值瞬时相位谱进行地震数据重构; 再利用方程(9)所示的时窗函数对重构地震数据进行改进的短时傅里叶变换, 得到重构数据的高精度时频瞬时振幅谱和高精度时频瞬时相位谱, 并抽取共频率的高精度时频瞬时振幅谱和高精度时频瞬时相位谱, 分别估算最负曲率; 最后沿地震解释层位做切片对比(图 7~图 10)。图 7a所示原始数据几乎不含随机噪声, 主要检验是否可以利用高精度时频峰值瞬时振幅谱及高精度时频峰值瞬时相位谱重构原始地震数据。图 7b为重构数据振幅切片, 对比可见, 除蓝色箭头标注处外, 重构数据和原始数据几乎没有差异。

图 8a和图 8b分别为重构数据30Hz和50Hz的高精度时频瞬时振幅谱切片; 图 8c和图 8d分别为基于重构数据30Hz和50Hz的高精度时频瞬时振幅谱估算的最负曲率切片。对比图 8c和图 8d可见:图 8c主要反映大尺度地质体结构, 图 8d主要反映小尺度地质体结构, 两者互补。对比图 8a和图 8c可见:高精度时频瞬时振幅谱相对稳定区域的曲率精度较高, 高精度时频瞬时振幅谱变化剧烈区域的曲率精度较低。对比图 8b和图 8d得到的结论相同。

图 9a和图 9b分别为重构数据30Hz和50Hz的高精度时频瞬时相位谱切片, 图 10a和图 10b分别为基于重构数据30Hz和50Hz的高精度时频瞬时相位谱估算的最负曲率切片。对比图 10a和图 10b可见:相对于基于低频成分的高精度时频瞬时相位谱估算的曲率属性而言, 基于高频成分的高精度时频瞬时相位谱估算的曲率属性具有更高的分辨率。图 10a比图 8c具有更高的分辨率, 图 10b比图 8d具有更高的分辨率, 说明基于高精度时频瞬时振幅谱估算的曲率属性容易受地震数据振幅空间变化的影响, 不能完全反映地层的真实结构; 基于高精度时频瞬时相位谱估算的曲率属性较基于高精度时频瞬时振幅谱估算的曲率属性具有更高的精度。

4 结论

1) 改进的短时傅里叶变换可获得信号的高精度时频瞬时振幅谱、时频瞬时相位谱和时频初始相位谱。

2) 瞬时相位反映了地下构造的变化, 基于瞬时相位可精确估算地层倾角, 提高曲率属性的估算精度。但瞬时相位对噪声极其敏感, 需充分压制随机噪声并提高瞬时相位的精度。

3) 基于改进的短时傅里叶变换获得的高精度时频瞬时相位谱估算的多尺度曲率属性可提高宏观及微观地下构造的识别精度。为降低随机噪声对时频瞬时相位谱的影响, 进一步提高基于高精度时频瞬时相位谱的多尺度曲率属性的估算精度, 可基于改进的短时傅里叶变换获得的高精度时频峰值瞬时振幅谱和高精度时频峰值瞬时相位谱对原始信号进行重构, 再利用重构信号的高精度时频瞬时相位谱进行多尺度的曲率分析。