基于高精度时频瞬时相位谱的多尺度曲率及其应用

李 勇,张固澜,何承杰,王晓凯,王婷一,殷 瑛,张建军

(1.西南石油大学地球科学与技术学院,四川成都610500;2.西安交通大学电子与信息工程学院,陕西西安710049;3.中国石油集团东方地球物理公司新兴物探开发处,河北涿州072751)

几何属性是最重要的地震属性[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]可表示为:

(1)

(2)

方程(1)可写为:

(3)

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

(4)

且

(5)

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

(6)

因此,有

(7)

1.2 合成记录分析

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

(8)

且

(9)

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

(10)

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

(11)

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

(12)

图1a为合成记录,图1b和图1c分别为改进的短时傅里叶变换得到的高精度时频瞬时振幅谱A1(τ,f)和A0(τ,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时,差值几乎为零,较大的差值主要分布在低频和高频端。

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

图1 合成记录时频分析a 合成记录; b 改进的短时傅里叶变换得到的高精度时频瞬时振幅谱(n=1); c 改进的短时傅里叶变换得到的高精度时频瞬时振幅谱(n=0); d 小波变换得到的时频瞬时振幅谱; e 短时傅里叶变换得到的时频瞬时振幅谱; f 图1b与图1d的差值; g 图1b与图1d的比值

图2 合成记录时频分析a 改进的短时傅里叶变换得到的高精度时频瞬时相位谱; b 小波变换得到的时频瞬时相位谱; c 图2a与图2b的差值; d 改进的短时傅里叶变换得到的高精度时频初始相位谱; e 短时傅里叶变换得到的时频初始相位谱; f 图2d与图2e的差值

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所示的瞬时相位计算得到的地层倾角,与理论模型的倾角信息完全匹配。

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

图3 合成地震记录及相关属性分析a 合成地震记录; b 瞬时相位; c 基于振幅的地层倾角信息; d 基于瞬时相位的地层倾角信息

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

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

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

(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余弦信号完全一致,说明可基于高精度时频峰值瞬时振幅谱和高精度时频峰值瞬时相位谱重构信号。

图4 合成信号瞬时属性分析a 无噪信号; b 无噪信号的瞬时振幅; c 无噪信号的瞬时相位; d 含噪信号; e 含噪信号瞬时振幅; f 含噪信号瞬时相位

图5 无噪信号的高精度时频瞬时振幅谱(a)和相位谱(b)与含噪信号的高精度时频瞬时振幅谱(c)和相位谱(d)

图6 含噪信号的高精度时频峰值瞬时振幅谱(a)和瞬时相位谱(b)及其重构信号(c)

3 应用实例

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

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

图7 原始数据振幅切片(a)与重构数据振幅切片(b)

高精度时频瞬时振幅谱相对稳定区域的曲率精度较高,高精度时频瞬时振幅谱变化剧烈区域的曲率精度较低。对比图8b和图8d得到的结论相同。

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

图8 沿层切片a 重构数据30Hz的高精度时频瞬时振幅谱; b 重构数据50Hz的高精度时频瞬时振幅谱; c 重构数据30Hz高精度时频瞬时振幅谱估算的最负曲率; d 重构数据50Hz高精度时频瞬时振幅谱估算的最负曲率

图9 沿层切片a 重构数据30Hz的高精度时频瞬时相位谱; b 重构数据50Hz的高精度时频瞬时相位谱

图10 沿层切片a 重构数据30Hz高精度时频瞬时相位谱估算的最负曲率; b 重构数据50Hz高精度时频瞬时相位谱估算的最负曲率

4 结论

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

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

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