一种改进的最小二乘时频分析方法

杨培杰,隋风贵

(1.中国石油化工股份有限公司胜利油田分公司勘探开发研究院,山东东营257015;2.中国石油化工股份有限公司胜利油田分公司,山东东营257001)

时频分析[1]在地球物理勘探中的应用十分广泛。目前比较常用的时频分析方法主要有短时傅里叶变换[2]、小波变换[3]、S变换[4]、希尔伯特-黄变换[5]、二次型时频分布[6]和匹配追踪分解[7]等。短时傅里叶变换相当于地震道在一个固定的时间窗上与正弦基的互相关,一旦窗函数选定,时频分辨率便可确定下来,这使得短时傅里叶变换对非平稳地震信号的分析存在局限性;小波变换是地震道与小波字典的互相关,是对短时傅里叶变换的改进,小波变换受测不准原理的影响,需要在时间分辨率和频率分辨率之间进行权衡,因此很难满足精细储层预测的要求;S变换是短时傅里叶变换和小波变换的继承和发展,S变换具有良好的时频特性,在时频分辨率上有所改进;希尔伯特-黄变换将信号分解为有限的几个固有模态函数,进而得到信号的时频谱,其不足是对噪声较敏感;二次型时频分布在时频平面存在非常严重的交叉项,导致其在地震勘探中的应用并不广泛;匹配追踪分解是将地震信号看作由许多被称作为原子信号的加权总和,将每个原子信号的频谱进行叠加就可得到地震道的时频谱。

最小二乘法在地震勘探中的应用十分广泛,约束最小二乘地震反演[8]使用最小二乘方法构建目标函数,极大提髙了解的收敛速度,其优点是反演过程稳定、计算效率高;最小二乘偏移[9-10]是一种基于线性化反演理论的真振幅成像方法,它将成像问题当作一个反问题来处理,以寻求更接近于真实的地下反射系数;罗水亮等[11]提出了一种基于变形Pride模型和Lee模型(P-L模型)的矩阵方程迭代精细横波预测方法,其核心是通过最小二乘将横波估计目标函数的最优化问题转化为求解线性矩阵方程组问题,从而极大地提高了横波估计的效率和稳定性。

本文在前人研究成果[12-16]的基础上,提出了改进的约束最小二乘时频分析方法(improved constrained least squares time frequency analysis,ICLS)。基于反演的思路将频谱分析转换成一个求解反问题的过程,通过正则化约束最小二乘方法求解该问题,进而得到地震信号的时频分析结果,理论模型和实际应用验证了本文方法的有效性。

1 方法原理

定义如下由待求的频谱和不同频率的三角函数计算信号的正演过程:

d=Sh

(1)

式中:d表示待分析的信号,在实际应用中,d往往是实数,为了提高时频估计结果的稳定性和分辨率,可以通过Hilbert变换将d转换成复数的形式;S表示不同频率的三角函数构成的矩阵;h表示待求的频谱向量。已知信号d和矩阵S求频谱h,这是一个典型的反问题,这里用正则化最小二乘的思路进行求解。

S的具体形式如下:

S=

(2)

式中:Δt表示时间采样间隔;β表示时间采样点数;Δf表示频率采样间隔;α表示频率采样点数。

对于精确的数据,待求的频谱向量为:

h=S-1d

(3)

但在实际计算过程中,正演方程((1)式)往往是超定的或是欠定的,所以没有精确解,即:

d=She+e

(4)

式中:he为计算的频谱结果;e表示实际数据中存在的误差。

通过普通最小二乘法即可以得到he的解:

he=(S′S)-1S′d

(5)

式中:S′表示S的共轭转置矩阵。

进一步,为了得到信号的时频谱,需要以待分析点为参考,对待分析信号加窗函数,目的是提取待分析点处前后的若干点,当待分析的信号加上窗函数后,矩阵S的元素变得相关,因此,为了获得稳定的解,需要加约束,引入对角矩阵Z:

(6)

式中:Nw=ceil(β/2),ceil表示将实数向无穷远的方向舍入到最接近的整数;Diag表示求对角矩阵。此时,(1)式可以进一步写成:

dw=Swh

(7)

式中:dw=Zd;Sw=ZS。

利用最小二乘法求解(7)式。为了进一步提高求解的稳定性和唯一性,需要对求解过程加约束[17],由于L2范数正则化可以有解析形式的解,更易于求解,因此选用L2范数正则化约束最小二乘法(附录A),最终可以得到ICLS时频分析结果:

he=S′w(SwS′w+μΩ)-1dw

(8)

式中:Ω表示约束对角矩阵[13];μ为权重系数,用于调节约束的大小,实际求解过程中,可不断调节μ的大小,直到获得一个时间分辨率和频率分辨率折中的时频分析结果。加入μΩ约束分量后,对于信噪比较低、时间段较短的信号的频谱分析会有更好的效果,可以有效地提高频谱分析过程的稳定性。逐点加窗计算(8)式即可得到待求信号段的时频谱。

进一步,可以通过公式(9)方便地实现ICLS反变换:

dinv=She

(9)

式中:dinv为经过反变换后重构的信号。

2 模型试算

设计一个由3个信号叠加而成的复合信号,如图1 所示,3个信号分别为平稳单频信号①、频率随时间变化的信号②和加高斯窗函数的单频高频信号③。

图1 复合信号

采用本文方法(ICLS)对该复合信号进行时频分析,如图2所示,从图2中可以看出,时频分析结果客观、准确地描述了该复合信号频率成分的变化,分辨率高。图3为复合信号的瞬时振幅及本文方法时频谱的局部放大显示,可以看出,时频谱为零的地方对应着复合信号瞬时振幅零点的位置,时频谱非零值对应着复合信号的瞬时振幅非零值的位置,也就是说,本文方法时频谱不仅可以描述频谱随时间的变化,同时也描述了原始信号瞬时振幅的大小。图4为本文方法反变换复合信号重构结果,可以看出,重构后复合信号与原始复合信号几乎重合,说明了本文方法的准确性。

图2 复合信号的ICLS时频谱

图3 复合信号瞬时振幅及其ICLS时频谱

图4 ICLS反变换

3 实际应用

实际资料来自胜利油田某工区,原始剖面如图5所示,在1.7s左右、CDP100～CDP200处有一背斜油气藏。分别采用本文方法(ICLS)和广义S变换(STFT)[18]对该剖面进行时频分析,结果如图6至图8 所示。从图6至图8可以看出,ICLS和STFT的时频分析结果具有较好的一致性,STFT的时频分析结果在频率较高的频谱分量上分辨率也较高,但在低频的频谱分量上分辨率较低,这已经无法满足目前隐蔽油气藏精细勘探的要求,相比而言,ICLS时频分析结果的分辨率更高、准确性更好,如对10Hz频谱分量依然具有很好的聚焦能力。

图8 STFT和ICLS的时频谱分析结果(50Hz)

图6 STFT和ICLS的时频谱分析结果(10Hz)

图5 原始地震剖面

图9为ICLS时频分析结果的10Hz频谱分量除以50Hz频谱分量的频谱比剖面,箭头处数值较大,可以明显看出该背斜存在低频增加、高频衰减的现象,指示了一套气藏的存在。

图7 STFT和ICLS的时频谱分析结果(30Hz)

图9 频谱比剖面(10Hz/50Hz)

4 结论

基于正则化约束最小二乘估计的思路,提出了一种改进的最小二乘时频分析方法。模型试算结果表明,采用本文方法可以得到准确的、高分辨率的时频分析结果;实际应用结果表明,本文方法对于地震信号的高中低频部分都有很高的分辨率和稳定性,可以有效地提高储层预测的准确度和可靠性,在隐蔽油气藏精细勘探中具有广阔的应用前景。