高维地震数据Wiener中心滤波方法

张力起,张 猛,2,王华忠,秦广胜

(1.波现象与智能反演成像研究组(WPI),同济大学海洋与地球科学学院,上海200092;2.中国石油化工股份有限公司胜利油田分公司物探研究院,山东东营257000;3.中国石油化工股份有限公司中原油田分公司物探研究院,河南濮阳457001)

近十余年,油气地震勘探技术的显著进展主要体现在“两宽一高”地震数据采集技术的普遍使用和地震波反演成像方法技术的逐渐实用化。噪声在地震波反演成像中影响巨大,它决定了反演成像方法的收敛性及成像质量。“两宽一高”地震数据采集得到了(尽可能)无假频且高维的叠前地震数据,为在高维数据(信号)空间中进行彻底的噪声压制提供了良好的数据基础。

勘探地震中的地震波场,可以表达成信号、相干噪声和随机噪声3部分的叠加,可记为D=S+Nc+Nr,其中,D为观测数据,S为地震信号,Nc为相干噪声,Nr为随机噪声。噪声压制的思想框架通常是将含噪信号视为Hilbert空间中的一个平方可积函数,选取合适的基函数对信号进行稀疏、近似的表达(信号预测的正问题);对上述信号预测正问题,在一定的估计理论基础上,求解最佳的信号预测器,从而实现信号的最佳估计。在这样的函数逼近思想下,地震数据的去噪、插值、压缩和特征提取等正问题是统一的。

在Bayes估计理论框架下,假定信号可以用既定的预测器进行预测,噪声满足某种统计分布,信号自身满足一些先验信息(例如具有稀疏性),可以通过后验概率最大化估计出最佳预测信号。若假定信号是线性可预测的且噪声满足高斯分布,那么可以在均方误差或L2范数误差最小准则下对信号实现最佳预测。常见的方法包括f-x滤波[1]、奇异谱分析(依据观测数据的自相关函数分解)[2]。信号在一些变换空间如小波基、曲波基、物理小波基和字典基[3-5]中,具有更好的稀疏性,将其与L0/L1范数下稀疏反演理论[6]结合可以更好地实现信号预测,达到信号分析的目的。这是当前现代信号分析的重点发展方向。此外,按照地震数据(随机信号)的统计特征设计滤波器是另一种重要的滤波方法,例如高斯滤波器、中值滤波器。将此类统计滤波方法扩展到高维数据的情况是必然的,例如非局部均值滤波器[7]、三维块匹配滤波器[8-9],此时数据或图像中结构信息在构建统计滤波器时是非常关键的。据此,基于信号逼近理论的滤波方法和基于信号结构信息的统计滤波方法逐渐有融合的趋势。

基于AR模型的预测滤波方法通常是在频率空间域(f-x)采用前向或后向预测滤波方法,或利用前向(后向)预测算子的共轭对称性同时进行后向(前向)预测滤波,从而达到压制噪声的目的。GÜLÜNAY[10]将非因果滤波(类似高维Wiener中心滤波)应用于三维叠后数据的去噪;WANG[11]将f-x域地震数据插值推广至三维的f-x-y域进行数据插值;LIU等[12]利用ARMA模型和非因果空间预测算子进行f-x-y滤波。

“两宽一高”地震勘探中采集的巨量地震数据(一块上百平方千米的工区叠前数据达到十几到几十TB),对高维空间、高效且高质量的噪声压制技术的研发有着迫切的需求。本文未采用更复杂的信号预测器(譬如高维小波变换),而是将AR模型的f-x滤波器(前向或后向预测滤波器)修改为Wiener中心预测滤波器,同时将高维数据排成Hankel矩阵的形式使其与Wiener中心预测滤波器对应,在预测误差的L2范数定义下求解对应的法方程获得最佳Wiener中心滤波器系数来构建高维线性预测滤波器,实现最佳滤波。对于含有复杂波场的地震数据,为满足局部线性假设,采用对数据局部取窗的方式进行去噪。数据窗长度以及滤波算子长度的选取会明显地影响去噪效果和执行效率。为了提高本文方法的计算效率,利用MPI进程级并行和OpenMP线程级并行来加速计算,将本文方法应用于合成数据和实际数据,结果表明其有良好的去噪效果。

1 Bayes框架下的地震数据分析理论

地震数据预处理(地震信号处理)中包含的大多数问题,譬如去噪、插值、多次波压制等正问题的构建都可以在Bayes反演理论框架下统一。实际采集的地震数据被认为是随机信号,满足一定的概率分布。地震数据去噪的正问题是建立对数据中所蕴含信号的预测理论,通常利用函数逼近理论建立包含参数的信号预测模型。预测结果与实测数据的残差是随机的,满足一定的先验概率分布。因此,随机地震观测数据可表示为:

(1)

(2)

式中:P(m|dobs)=P(dobs|m)P(m)/P(dobs),P(m|dobs)为后验概率密度函数,其中,P(dobs|m)为实测数据dobs对模型参数m的先验概率密度函数,P(m)为模型参数m的先验概率密度函数,P(dobs)为实测数据dobs的先验概率密度函数。P(m|dobs)可以理解为当观测数据为实测数据dobs时,对应不同的参数m时的概率。在高维情况下,后验概率密度函数P(m|dobs)几乎无法由实际计算得到。因此,常取后验概率密度函数P(m|dobs)最大值时对应的参数m为最终的估计结果:

(3)

一般地,在假设P(dobs|m)和P(m)这两个概率密度函数都是高斯函数且正问题是线性的情况下,有:

(4)

式中:C为常数;S(m)为代价函数。

其中,代价函数S(m)为:

(5)

式中:mprior是模型参数m的先验值;CD是数据自相关矩阵;CM是模型参数自相关矩阵。采用各种数值优化算法求解公式(5),均可得到估计的模型参数m。

信号预测算子G(·)可记为G,它由各种选定的基函数族中的基函数线性叠加产生。基函数族可分为两类,第1类是预先选定的基函数簇,包括Fourier基函数、余弦基函数、Radon基函数、小波基函数和曲波基函数等;第2类是由数据驱动获得的基函数簇,包括K-L变换、奇异谱分析和字典学习获得的基函数等。模型参数m为线性信号预测器的系数。基函数实质上是信号的特征成分。基函数选择恰当,则可以用很少的特征成分组合很好地表达信号,这是信号(与图像)分析中最根本的基础。在上述Bayes框架下,对模型参数进行最佳估计后,即可实现对数据中所蕴含信号的最佳表达。

上述抽象的理论框架建立了数据分析的理论基础,原则上适用于任何信号和图像分析。本文中提出的高维Wiener中心滤波方法则是基于信号是线性可预测且噪声满足高斯分布的假设,在L2范数误差最小准则下对信号进行建模和最佳预测。

2 f-x域中线性信号模型的建立与信号预测

Wiener滤波可实现对线性信号或相干噪声的最佳预测,达到压制随机噪声(也包括相干噪声)和提高信噪比的目的。在频率空间域,这些线性信号可以表示成具有线性相位移的谐波叠加。时间空间域中的二维信号s(t,x)在f-x域表示为:

(6)

式中:p为水平射线参数;Δx为道间距;m为道数;Δt为相邻两道间的时移量;A(f)为子波的频谱;φ=2πfpΔx,φ代表由p确定的线性同相轴的线性相位移。

对于含有一个线性同相轴的任意一个单频,相邻道之间存在如下关系:

(7)

式中:a=ei2πfpΔx;Si-1为频率域第i-1道的值;Si为频率域第i道的值。(7)式说明各道信号之间存在可预测性。类比AR模型,当含有p个线性同相轴时,可以得到前向线性预测滤波器{gi},i∈[1,p]:

(8)

(9)

(10)

图1 算子长度为3时前向预测滤波(圆圈)和后向预测滤波(方框)示意

(11)

(12)

为求解滤波器系数g,需要建立如下误差泛函J(g):

(13)

(13)式的法方程为:

(14)

(15)

式中:Df为前向预测滤波器构建的矩阵;gf为前向预测滤波器系数;Db为后向预测滤波器构建的矩阵;gb为后向预测滤波器系数。

3 高维Wiener中心滤波

根据实际数据各自构建法方程,然后求解得到的前向预测滤波器和后向预测滤波器(基于AR模型构建),只适用于噪声满足高斯分布、振幅不随空间变化的局部线性信号预测,若不满足该条件,则在压制噪声的同时会破坏信号的振幅和相位信息。Wiener中心预测方法适用于振幅缓变的局部线性信号。利用数据的多级Hankel矩阵排布来构建Wiener中心预测滤波器的法方程,进而实现高维数据的预测滤波,该方法称为高维Wiener中心滤波方法。

3.1 Hankel矩阵的构建

对于一维数据(向量),其构建的Hankel矩阵是指每一条交叉对角线上的元素都相等的矩阵,即Hankel矩阵的项hi,j=hi-1,j+1。对于高维数据,也同样能够构建出Hankel矩阵,即块状Hankel矩阵(block Hankel matrix)。多维信号以二维数据为例,D∈Rm×n所代表的矩阵为:

(16)

(17)

(18)

式中:l2+k2-1=m。图2为将3×3的二维数据排成块状Hankel矩阵结构。

图2 由二维数据构建的块状Hankel矩阵结构

3.2 Wiener滤波中心预测滤波

将二维数据变换到f-x域后,单个频率片数据为一维数据。采用图3中的一维Wiener中心预测滤波方法进行预测滤波。假定数据长度为6,Wiener中心滤波算子长度为2,一维Wiener中心滤波器的构建过程如图4所示,可以看出,只需要对一维数据进行Hankel排布就可以得到公式(14)中的D和S。

构建多级Hankel矩阵可以将一维Wiener中心滤波方法推广到高维。以二维Wiener中心预测滤波

s

s

k,l

(19)

图5 二维Wiener中心预测滤波示意(算子大小为5×5)

图6 根据单频空间二维数据矩阵构建D和向量d

若是按行优先的方式进行数据的Hankel矩阵排布,首先对窗内(红色框)数据按行优先方式排成一行;然后将窗沿行方向滑动一个单位长度;再将窗内数据按行优先方式排成一行,依照上述方式滑动窗并取数据进行排布;直到窗滑动到数据末端;再将窗按列滑动一个单位长度,重复上述按行方向滑动窗取数据并排布的操作;重复以上过程直到窗滑动到数据末端。

高维Wiener中心滤波的计算步骤与f-x滤波计算步骤类似,具体如下:

1) 将数据从时间空间域变换到频率空间域;

2) 对每个频率切片数据进行Hankel矩阵排布,构造法方程,估计滤波器系数;

3) 利用滤波器系数对变换后的复数域数据进行褶积运算;

4) 将滤波后的结果变换回时间空间域。

为了满足信号的空间局部线性假设条件,对数据进行取窗处理。因三维地震数据的单个频率片数据为二维复数域(如图7蓝色框所示),故对二维复数域数据取窗数据(图7中黑色框所示),对缺少的数据进行补零后,将其排成Hankel矩阵,最终滤波后输出数据的区域如图7红色框所示,具体算法流程如图8所示。

图7 二维单个频率片数据的I/O示意

图8 实际数据处理流程

为了满足局部线性假设条件,需对地震数据进行取窗处理。不同的数据窗长度和不同的算子长度会造成去噪结果和效率的明显差异。长数据窗会破坏局部线性假设,无益于提高去噪能力,同时还会增加计算量;而短数据窗则导致该算法去噪能力下降,无法适应复杂波场的情况。因此,需要根据数据波场的复杂程度选取合适的数据窗长度和算子长度。一般情况下,算子长度选取为3或5,空间窗长度选取范围为15～30,时间窗长度选取范围为100～200个采样点。由于采用了分块处理数据的方式并且预测滤波是在单个频率片上进行的,故可以采用并行策略加速运算。并行主要利用MPI和OpenMP进行加速。利用MPI多进程分块读取高维窗数据体,各个进程对各自的数据体独立进行高维去噪;在每个进程均利用OpenMP多线程分别处理单个频率片数据。

4 测试应用及分析

对三维数据体取窗后获得的数据体同样是三维的,而频率空间域中滤波算子的维度则是二维的。公式(20)为信噪比(SNR)的定义。

(20)

式中:{xi},i∈[1,N]为无噪声数据;{ni},i∈[1,N]为含噪声数据;N为数据长度。

4.1 三维合成地震数据

图9中合成的三维数据大小为600×60×60,时间采样间隔为4ms,空间采样间隔为1m,子波主频为30Hz。图9a为三维合成数据中某条测线的地震记录,图9b为三维合成数据加入-6dB高斯白噪声后该测线的地震记录(对应图9a同一条测线)。采用不同长度的一维前向和后向预测滤波器沿两个方向各做一次一维预测滤波,其滤波结果分别如图9c、图9e 所示,数据残差剖面中均出现了强同相轴(图9d、图9f),说明这种滤波方法破坏了信号结构,压制噪声的能力差。分别采用不同大小的数据窗和二维Wiener中心滤波算子进行预测滤波,结果如图10 所示。可以看出,当数据窗(100×10×10)和Wiener滤波算子(3×3)较小时,滤波效果较差;当数据窗(100×30×30)和Wiener滤波算子(5×5)较大时,其数据残差剖面包含的有效信号较小且信噪比较大,滤波效果较好。

图11中合成的三维合成单炮炮集数据大小为1500×120×120,时间采样间隔为4ms,空间采样间隔为10m,子波主频为30Hz。图11a为三维合成炮集数据中的某条测线的地震记录,图11b为加入-10dB高斯白噪声后三维合成炮集数据的地震记录(对应图11a同一条测线)。采用二维Wiener中心滤波方法进行预测滤波,图11c显示了数据时窗大小为300×30×30,滤波算子大小为3×3时的二维Wiener中心滤波后的地震记录,图11d为采用二维Wiener

图9 三维合成数据及一维前向和后向预测滤波结果a 三维合成数据中某条测线的地震记录; b 加入-6dB高斯白噪声后三维合成数据中某条测线的地震记录(对应图9a同一条测线); c 数据窗大小为100×10×10,算子长度均为3时前向预测和后向预测滤波后的结果; d 图9c与图9a的数据残差剖面(SNR为15.69dB); e 数据窗大小为100×10×10,算子长度均为5时前向预测和后向预测滤波后的地震记录; f 图9e与图9a的数据残差剖面(SNR为12.98dB)

中心滤波方法得到的图11c与图11a的数据残差剖面,图中无明显的连续的同相轴,并且噪声幅值小。可以看出Wiener中心滤波方法在压制噪声的同时较好地保留了有效信号。

4.2 三维实际地震资料

图12a为某地区实际采集的四维数据(叠前CDP道集)中某个CDP道集(三维数据体),其数据大小为1501×120×34。当数据窗大小为200×30×20,滤波算子为3×3时,采用二维Wiener中心滤波方法对数据进行预测滤波,得到的CDP道集和残差剖面分别如图12b和图12c所示。可以看出二维Wiener中心滤波方法能够压制随机噪声和部分相干噪声,在振幅存在局部变化时(图12a和图12b中的红圈)能较好地保留信号。

图11 三维合成炮集数据二维Wiener中心滤波结果a 三维合成炮集数据中某条测线的地震记录; b 加入-10dB高斯白噪声后三维合成炮集数据中某条测线的地震记录(对应图11a同一条测线); c 数据窗大小为300×30×30,滤波算子大小为3×3时二维Wiener中心滤波后的地震记录; d 图11c与图11a的数据残差剖面(SNR为12.24dB)

图12 实际数据二维Wiener中心滤波结果a 野外采集得到的某个CDP道集; b 数据窗大小为200×30×20,滤波算子大小为3×3时二维Wiener中心滤波后的CDP道集; c 图12b与图12a的数据残差剖面

某地区的三维实际数据叠后剖面如图13a所示,数据大小为1501×567×10。当数据窗大小为300×20×12,滤波算子为3×3时的采用二维Wiener中心滤波方法对数据进行预测滤波,得到的叠后剖面和残差剖面分别如图13b和图13c所示。可以看出,二维Wiener中心滤波方法能够压制中、深层杂乱的噪声,并且很好地保留浅、中层连续的同相轴以及断层信息。

图13 实际数据二维Wiener中心滤波结果a 三维野外实际数据叠后剖面; b 数据窗大小为300×20×12,滤波算子大小为3×3时二维Wiener中心滤波后的叠后剖面; c 图13b与图13a的数据残差剖面

5 讨论与结论

“两宽一高”地震勘探针对的是小尺度的精细油藏和复杂油藏,相应的成像方法也逐渐进入以FWI+LS_RTM为标志的反演成像阶段。噪声对成像的影响需要被重新审视。偏移成像中的随机噪声可通过高覆盖的成像叠加被压制,对成像结果影响较小;而反演成像过程中噪声则会严重影响反演成像方法的收敛性及反演结果的精度(分辨率)。另外,针对当前“两宽一高”地震数据采集得到的巨量数据体,在实际生产中实用的、高维空间中可实现的且效率和效果达到均衡的去噪方法还是相当缺乏的。

本文所讨论的高维Wiener中心滤波方法,是由基于AR模型中的前向预测滤波器和后向预测滤波器发展得到的,其要求信号是线性可预测的并且噪声满足高斯分布的。对于非线性的信号,可采用对数据局部取窗的方式满足局部线性可预测的假设条件;对于噪声不满足高斯分布假设的信号,需要修正Bayes理论框架下的信号预测理论和方法以达到更好的去噪效果。地震数据基本都是高维数据体,采用低维空间中的预测算子很难利用其高维空间信息,也无法达到满意的去噪效果。本文通过对数据体进行Hankel矩阵排布解决了这一问题,结合Wiener中心预测滤波方法,实现了三维或者更高维数据空间中的噪声压制。实验结果表明,高维Wiener中心滤波方法能够在有效压制噪声的同时较好地保留信号,但存在计算效率低的问题。本文采用MPI和OpenMP并行计算可以部分地提高计算效率。但是本文方法也存在一些问题,如对于非规则数据、复杂波场(振幅变化较为强烈)以及在信噪比相对较低的情况下则很难最佳地预测信号并压制噪声。

到目前为止,地震数据去噪(插值)真正的问题依然是所建立的地震信号预测算子不能很好地预测实测数据中的信号,以及噪声的概率分布是未知的。当前,信号建模的合理性、参数估计的精度、算法计算效率等方面远未达到理想的程度,这使得去噪结果达不到反演成像的要求。“两宽一高”的数据采集方式提供了更好的数据基础,但(尤其在复杂近地表探区)地震数据在空间时间上的复杂变化使得当前已有的去噪技术并没有达到令人满意的、满足反演成像要求的程度,未来仍需要继续深入研究。

致谢:感谢中石油勘探开发研究院及西北分院、中海油研究院和湛江分公司、中国石油化工股份有限公司石油物探技术研究院和胜利油田分公司对波现象与智能反演成像研究组(WPI)研究工作的资助与支持。