基于伪WVD和CZT的地震信号处理方法

摘要

本发明公开了一种基于伪WVD和CZT的地震信号处理方法，包括以下步骤：S1、对地震信号作Hilbert变换，得到频域信号，将地震信号和频域信号相加，得到解析信号；S2、得到PWVD时频分布结果；S3、得到PWVD的CZT时频分布结果；S4、将PWVD时频分布结果按照频率重排，得到地震信号的高分辨率单频剖面图；对CZT时频分布结果作逆变换，取逆变换后第一行与信号相量相除并取相除结果的实部，得到重构的时域信号。本发明先将原始信号变换到频域，接着在频域中进行高分辨率处理，最后返回到时域，得到高分辨率的时频谱图和重构地震信号。重构信号更为平滑，能较好的反映出薄层的能量分布状况；提升了地震信号对储层边缘细节的刻画能力。

说明书

技术领域

本发明属于地球科学领域，特别涉及一种基于伪WVD和CZT的地震信号处理方法。

背景技术

地震勘探技术是油气勘探的最常用的手段。地震勘探是利用地下介质弹性和密度的差异，通过观测和分析人工地震产生的地震波在地下的传播规律，推断地下岩层的性质和形态的地球物理探勘方法。随着石油勘探程度不断提高，油气勘探的重点转为复杂构造、隐蔽岩性油气藏等类型，所探查的地质结构越发复杂，表现出圈闭规模小、储层厚度薄等特点。整个勘探过程对地震信号的分辨率要求越来越高，而地震波是一种典型的非平稳信号，其能量分布主要存在于低频区域；类似傅里叶变换这样的传统分析方法因其只能在时域或频域内作全局变换的特性，不能对地震信号作时频域的局域化分析，无法在处理非平稳信号中取得良好效果。气藏埋藏深、储层薄、地震能量和频带衰减快等一系列的原因，使得勘探人员难以准确识别勘探目标，很大程度上限制了油气勘探的深化和新目标的发现。地震信号的分辨率是油气勘探工作中获取储层信息的关键因素，对研究薄层地质结构有着重要的意义；因此，针对地震信号的特性，合理有效地提高其分辨率就显得尤为重要。

目前，常用的地震信号处理方法，包括短时傅里叶变换(STFT)、小波变换、Wigner—Ville分布、Cohen类分布等。短时傅里叶变换和小波变换都是线性分布，无法完整的表达出地震信号的时域和频域的分布关系，在分析地震信号时会有一定程度的缺陷。Wigner—Ville分布(WVD)、Cohen类分布都属于双线性分布，Wigner—Ville分布能够直观描述出信号的能量随时域和频域的分布，有着优秀的时频聚焦能力，但其也有着致命的缺点：会产生大量的交叉干扰项，不利于研究人员识别出地下储层构造。伪Wigner—Ville分布(PWVD)属于Cohen类分布的一种，其在Wigner—Ville分布的基础上，通过加入移动的窗函数来改变核函数，从而在很大程度上抑制交叉项的产生。虽然通过抑制交叉项的产生，伪Wigner—Ville分布得到了一个较好的时频分布结果，但其也无法进一步提高地震信号的分辨率。

线性调频Z变换(Chirp-z transform，CZT)是由Lawrence Rabiner提出的一种信号分析方法，最初用于分析语音信号。CZT在Z平面的单位圆上利用螺旋线进行等间隔采样，并且采样点在螺旋线上呈现等角分布，信号的谱分析可以在Z平面上的螺旋线上实现。这种方法计算功率谱的时频分布十分有效，并且有着计算快速、不受数据序列长度限制等优点，非常适合处理复杂的地震信号。通过调整采样间隔和采样点数，可以实现时频谱图的细化，进而提升信号的分辨率。

发明内容

本发明的目的在于克服现有技术的不足，提供一种先将原始信号变换到频域，接着在频域中进行高分辨率处理，最后返回到时域，得到高分辨率的时频谱图和重构地震信号，能够提升地震信号对储层边缘细节的刻画能力的基于伪WVD和CZT的地震信号处理方法。

本发明的目的是通过以下技术方案来实现的：基于伪WVD和CZT的地震信号处理方法，包括以下步骤：

S1、对地震信号作Hilbert变换，得到其频域信号，将地震信号和频域信号相加，得到解析信号；

S2、以时间差和时间序列为自变量，求出解析信号的瞬时自相关函数；对瞬时自相关函数乘上长度为2L-1的矩形窗，得到伪Wigner—Ville分布的核函数；将核函数进行循环位移得到新的核函数，对循环位移后新的核函数作傅里叶变换得到PWVD时频分布结果，其中，2L为地震时域信号的时间采样点；

S3、将循环位移后新的核函数乘上加权系数得到序列f

S4、将PWVD时频分布结果按照频率重排，得到地震信号的高分辨率单频剖面图；再对CZT时频分布结果作逆变换，然后取逆变换后数据的第一行与信号相量相除并取相除结果的实部，得到重构的时域信号。

进一步地，所述步骤S1包括以下子步骤：

S11、获取地震剖面，利用该地震剖面得到单道地震信号x(n)；

S12、将单道地震信号x(n)作Hilbert变换得到频域信号H[x(n)]，然后同时域地震信号相加，得到解析信号：

z(n)＝x(n)+jH[x(n)] (2)。

进一步地，所述步骤S2包括以下子步骤：

S21、以时间差l和时间序列n为自变量，其中l＝-L+1,…0,1,…L；n＝0,1,…2L-1；求出解析信号z(n)的瞬时自相关函数r(n,l)：

r(n,l)＝z(n+l)z(n-l)

式中，z(n+l)为n+l时刻的解析信号，z(n-l)

S22、将瞬时自相关函数r(n,l)乘上长度为2L-1的矩形窗h(l)，且当|l|＞L时，h(l)＝0；得到加窗后的离散瞬时自相关函数为R

R

式中，h(-l)

S23、令Z(n,l)＝z(n+l)h(l)，得到伪Wigner—Ville分布的核函数K

式中，Z(n,-l)

S24、对核函数K

式(6)中：Z(n,-l+2L)

S25、对核函数K

式(7)中：e

进一步地，所述步骤S3包括以下子步骤：

S31、将步骤S2中得到的PWVD的核函数K

式中，A

S32、对序列f

式(9)和式(10)中，Q是为求CZT所构造的序列长度，取值为最接近N+M-1的二次幂，N为离散信号的时间序列长度，N＝2L；M为CZT的采样点数；

S33、分别对g(l),f(l)作傅里叶变换，再将傅里叶变换得到的G(r)、F(r)在频域内相乘得到Y(r)：

G(r)＝fft[g(l)],F(r)＝fft[f(l)] (11)

Y(r)＝F(r)G(r) (12)；

S34、对Y(r)作逆变换，得到y(l)；最终对y(l)乘上系数后得到PWVD的CZT分布结果CZT(z

y(l)＝ifft[Y(r)] (13)

式中，

进一步地，所述步骤S4包括以下子步骤：

S41、将步骤S2中得到的PWVD时频分布结果PWVD

S42、将步骤S3中得到的CZT时频分布结果CZT(z

ifK

式中，ifft为逆傅里叶变换；

S43、将逆变换结果ifK

X

式中，x(n)为原始时域信号，re表示信号的实部。

本发明的有益效果是：本发明基于伪Wigner—Ville分布和Chrip-Z变换对地震信号进行处理，先将原始信号变换到频域，接着在频域中进行高分辨率处理，最后返回到时域，得到高分辨率的时频谱图和重构地震信号。本发明首次将CZT与PWVD相结合，并用于地震信号的分析当中；在PWVD的基础上，实现了地震信号时频谱图细化，进一步提升了信号分辨率。由CZT逆变换得到的重构时域信号与原始时域信号相比，重构信号更为平滑，能较好的反映出薄层的能量分布状况；同时，使得地层的层间关系更加明显，提升了地震信号对储层边缘细节的刻画能力，更利于识别出油气矿藏，为地下薄互层介质的高精度识别提供了更为精确的地震资料。

附图说明

图1为本发明的基于伪WVD和CZT的地震信号处理方法的流程图；

图2为CZT的变换路径示意图；

图3是CZT的计算过程示意图；

图4是仿真信号的时域波形图；

图5(a)、(b)、(c)、(d)、(e)、(f)分别为图4的仿真信号的短时傅里叶变换结果图、短时傅里叶变换结果的三维图、伪Wigner—Ville分布结果图、伪Wigner—Ville分布结果的三维图、Chrip-Z变换结果图、Chrip-Z变换结果的三维图；

图6是仿真信号经过基于PWVD的CZT算法处理后再由CZT作逆变换得到的重构信号；

图7为数值模拟地震信号的时域波形、STFT、PWVD以及CZT分布图；

图8(a)、(b)分别为数值模拟信号PWVD的剖面图和CZT的剖面图；

图9(a)、(b)分别为数值模拟信号的原始时域信号图和经过基于PWVD的CZT算法处理后再由CZT逆变换得到的重构信号图；

图10是Crossline数据的短时傅里叶变换剖面图；

图11为经过基于PWVD的CZT算法处理的Crossline数据单频剖面图；

图12为Inline数据的短时傅里叶变换剖面图；

图13为经过基于PWVD的CZT算法处理的Inline数据单频剖面图；

图14(a)、(b)、(c)、(d)分别是原始Crossline数据的时域信号、CZT逆变换后的Crossline时域信号、原始Inline数据的时域信号和CZT逆变换后的Inline时域信号。

具体实施方式

下面结合附图进一步说明本发明的技术方案。

如图1所示，本发明的一种基于伪WVD和CZT的地震信号处理方法，包括以下步骤：

S1、对地震信号x(n)作Hilbert变换，得到其频域信号H[x(n)]，将地震信号x(n)和频域信号H[x(n)]相加，得到解析信号z(n)；具体包括以下子步骤：

S11、获取地震剖面，利用该地震剖面得到单道地震信号x(n)；

S12、将单道地震信号x(n)作Hilbert变换得到频域信号H[x(n)]，然后同时域地震信号相加，得到含有时域和频域信息的解析信号：

z(n)＝x(n)+jH[x(n)] (2)。

S2、以时间差l和时间序列n为自变量，求出解析信号z(n)的瞬时自相关函数r(n,l)；对瞬时自相关函数r(n,l)乘上长度为2L-1的矩形窗h(l)，得到伪Wigner—Ville分布的核函数K

具体包括以下子步骤：

S21、以时间差l和时间序列n为自变量，其中l＝-L+1,…0,1,…L；n＝0,1,…2L-1；求出解析信号z(n)的瞬时自相关函数r(n,l)：

r(n,l)＝z(n+l)z(n-l)

式中，z(n+l)为n+l时刻的解析信号，z(n-l)

S22、将瞬时自相关函数r(n,l)乘上长度为2L-1的矩形窗h(l)，且当|l|＞L时，h(l)＝0；得到加窗后的离散瞬时自相关函数为R

R

式中，h(-l)

S23、令Z(n,l)＝z(n+l)h(l)，得到伪Wigner—Ville分布的核函数K

式中，Z(n,-l)

S24、由于傅里叶变换(FFT)的计算域为非负数(l＝0,1,2·····2L-1)，因此需要对核函数K

式(6)中：Z(n,-l+2L)

S25、对核函数K

式(7)中：e

S3、将循环位移后新的核函数K

具体包括以下子步骤：

S31、将步骤S2中得到的PWVD的核函数K

式中，A

S32、由PWVD的推导公式可知，K

式(9)和式(10)中，Q是为求CZT所构造的序列长度，取值为最接近N+M-1的二次幂，N为离散信号的时间序列长度，N＝2L；M为CZT的采样点数；

S33、分别对g(l),f(l)作傅里叶变换，再将傅里叶变换得到的G(r)、F(r)在频域内相乘得到Y(r)：

G(r)＝fft[g(l)],F(r)＝fft[f(l)] (11)

Y(r)＝F(r)G(r) (12)；

S34、对Y(r)作逆变换，得到y(l)；最终对y(l)乘上系数后得到PWVD的CZT分布结果CZT(z

y(l)＝ifft[Y(r)] (13)

式中，

本步骤所涉及到的CZT的原理如下：

已知有限长度的离散信号序列x(n)的Z变换为X(z)，有

令Z

式中A

由式(16)可知，当r＝0时，

CZT变换路径如图2所示。

对信号的进行频谱分析时，应在单位圆上实现CZT。A

由Bluestein公式：

令

式(20)即为CZT的定义，可以看出CZT是通过先对离散信号乘上系数

S4、将PWVD时频分布结果PWVD

具体包括以下子步骤：

S41、将步骤S2中得到的PWVD时频分布结果PWVD

S42、将步骤S3中得到的CZT时频分布结果CZT(z

ifK

式中，ifft为逆傅里叶变换；

S43、逆变换结果ifK

X

式中，x(n)为原始时域信号，re表示信号的实部。

下面通过具体信号对本发明的信号处理效果进行验证。

实例1仿真信号

为了验证明基于伪Wigner—Ville分布和Chrip-Z变换的地震信号处理提高分率算法的可行性和有效性。首先采用如下仿真信号，信号的峰值在0.2s、0.6s、1s处，分别对应的频率分量为40Hz、20Hz、10Hz，这是模拟了单道地震信号的结果；信号的采样频率为1000Hz，采样点数为1200。信号的时域波形如图4所示。

图5(a)为仿真信号的短时傅里叶变换结果；图5(b)为短时傅里叶变换结果的三维图；图5(c)为仿真信号的伪Wigner—Ville分布结果；图5(d)为伪Wigner—Ville分布结果的三维图；图5(e)为仿真信号的Chrip-Z变换结果；图5(f)为Chrip-Z变换结果的三维图。从仿真信号处理后的时频谱图可以看出，加入了窗函数的PWVD不仅有着优良的时频聚焦性还保留了WVD的数学特性，较好的抑制了干扰项的生成，得到了一个较为不错的时频分析结果。通过将CZT与PWVD结合，利用CZT设置不同的幅角参数的方式，不仅实现了细化仿真信号时频谱图，还进一步提升仿真信号的时频聚焦能力，从而提升了信号的时频分辨率。

图6为由CZT作逆变换得到的重构信号。可以看出，通过CZT的逆变换基本还原了原始信号的时域波形；由于仿真信号较为简单平滑和不含噪声，并且信号经过WVD、CZT后其幅度和相位信息也会发生变化，因此重构信号在0.2s、0.6s附近的峰值与原始信号有些许不一致。

实例2正演数值模拟信号

在仿真信号验证了算法的有效性后，接着将算法应用于数值模拟地震信号。所采用的理论地震信号采用二维CDP的方式存放，总共包含101个信道，采样周期为1ms，每一列信道有231个时间采样点，采样频率为1000Hz。

图7(a)～(c)分别为第1列信道的时域波形、第60列信道的时域波形和第101列信道的时域波形；图7(d)～(f)分别为第1列信道的STFT、第60列信道的STFT和第101列信道的STFT；(g)～(i)分别为第1列信道的PWVD、第60列信道的PWVD和第101列信道的PWVD；(j)～(l)分别为第1列信道的CZT分布、第60列信道的CZT分布和第101列信道的CZT分布。通过对正演数值模拟信号处理得到的时频谱图，可以看出在单道地震信号的PWVD的基础上，CZT通过设置采样频段和细化点数，在Z平面上螺旋采样，以采样的方式对时频分布进行插值处理，并且通过设置相应的采样间隔，从而进一步提升原始的PWVD的分辨率，达到了频谱细化的效果。

图8(a)为数值模拟信号PWVD的剖面图；图8(b)为数值模拟信号CZT的剖面图。对比可以看出通过本发明中步骤S2、S3的技术方案，实现了剖面图的频谱细化，提升了信号的时频分辨率。

图9(a)为数值模拟信号的原始时域信号，含有101道信号；图9(b)为CZT逆变换得到的重构信号。对比可以看出通过本发明实例中步骤S4的技术方案，由CZT得到的重构信号更加平滑，不仅有着更好的分辨率和能量聚焦性，且储层走向的局部边缘细节更加清晰。

实例3实际地震资料

基于以上理论模型研究，将本发明方法用于四川盆地西部的川西坳陷的试验区马井1井地下6142.0～6152.4米的上储层和位于6171.8～6228.0米的下储层的实际地震资料。该资料含有水平走向的Crossline数据和垂直走向的Inline数据，共1011道，时间采样点数为1001。

图10为Crossline的短时傅里叶变换剖面图；图11(a)～(f)为基于PWVD的CZT算法时频处理的Crossline数据单频剖面图，分别对应频率35Hz、45Hz、55Hz、65Hz、75Hz、85Hz；图12为Inline的短时傅里叶变换剖面图。图13(a)～(f)为经过基于PWVD的CZT算法处理的Inline数据单频剖面图,分别对应频率35Hz、45Hz、55Hz、65Hz、75Hz、85Hz。对比可以看出，在经过PWVD以及CZT后，不仅信号的能量分布有着更为明显的间隔，不同高度的频谱带明显变细，提升了时频分辨率，从而实现了Crossline、Inline数据的频谱细化。细化后的时频谱图能更好的分辨出储层的水平走向的地质构造，区分出储层的边缘细节，更利于识别出地下薄互层介质。

图14(a)为原始Crossline数据的时域信号；图14(b)为CZT逆变换后的Crossline时域信号；图14(c)为原始Inline数据的时域信号；图14(d)为CZT逆变换后的Inline时域信号。通过对比可以看出，相较于原始信号，逆变换的时域信号不仅增加了储层的边缘细节，使得不同水平高度的地层之间的间隔更加明显；同时也对储层能量分布相近的区域作出了区分，更利于识别储层。通过对细化后的时频图作逆变换，从而使得时域信号的精度有了提升。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。