基于稀疏反演的多道匹配追踪地震信号去噪方法及其应用

任 浩,李宗杰,薛 姣,顾汉明

(1.中国地质大学(武汉)地球内部多尺度成像湖北省重点实验室,湖北武汉430074;2.中国石油化工股份有限公司西北油田分公司勘探开发研究院,新疆乌鲁木齐830011)

地震数据不可避免地包含噪声,噪声会影响到地震属性提取和地震解释等后续工作,因此压制噪声是地震数据处理过程中的关键环节[1]。

增强有效信号或减弱噪声通常包括两个步骤:①定义一个可区分有效信号与噪声的特征;②研发基于该定义的分离算法[2]。考虑到相邻地震道中有效信号是相关的,噪声是非相关的,基于该特征进行同相位叠加可增强有效信号,并减弱噪声。地震信号中有效信号和噪声在时间域是混叠的,不易分离,谱分解技术将时域地震信号映射到其它域后再进行处理,取得了良好的效果。目前主流的谱分解技术包括傅里叶变换、经验模态分解(empirical mode decomposition,EMD)、同步压缩小波变换(synchrosqueezing wavelet transform,SWT)、小波阈值法(wavelet thresholding approach,WTA)、匹配追踪算法(matching pursuit,MP)和反演谱分解(inverse spectral decomposition,ISD)等。傅里叶变换是将地震信号转换到频率域,根据噪声频率与有效信号频率范围的差异特征,设置合适的滤波器压制噪声,但这种方法对非平稳地震信号的处理效果不佳[3]。经验模态分解[4]是将地震信号表示成多阶固有模态函数分量之和,低阶分量频率高,高阶分量频率低,去除包含高频噪声的低阶分量即可实现去噪,但该方法造成了部分有效信号损失,针对这一问题,相继有学者提出了改进的方法[1,5-8]。同步压缩小波变换[9]具有高时频分辨率特性和信号重构能力,该方法根据有效信号的时频分布范围,将时频谱中相应时频成分重构有效信号以达到去除噪声的目的[10-12]。小波阈值法[13]是将地震信号变换到时间-尺度域,利用噪声能量相对较弱的特点,将小于阈值的系数设置为零,由新的小波系数经小波反变换得到去除噪声后的信号。匹配追踪算法[14]是将地震信号分解为一系列子波的线性组合,赵天姿等[15]人工挑选出与有效信号匹配的子波重构信号,后来有学者依据多个连续地震道中有效信号的横向相关性,提出了多道匹配追踪地震信号分解方法[16-17]。反演谱分解[18]与匹配追踪算法类似,也是一种基于子波库的分解方法,该方法预先设定误差水平值,将地震信号表示成稀疏时频子波的线性相加,将残差信号当作噪声[19-20]。

匹配追踪分解是采用逐步搜寻子波库中与残差信号最佳匹配子波的方式分解地震信号,并非直接且准确地去除噪声。搜寻最佳匹配子波是分解中最为关键的步骤,本文方法以子波库为着手点来解决这一问题,采用与地震信号中有效信号匹配的子波构建的子波库代替超完备子波库,将地震信号表示为有效信号匹配子波的线性叠加与残差信号之和,将残差信号当作噪声,即达到直接去除噪声的目的。

1 基于稀疏反演的多道匹配追踪分解原理

如何实现直接分离有效信号和噪声？首先,以连续多道的中间道为标准,校正同相轴为水平后同相叠加求取叠加道;接着,稀疏反演分解叠加道,选取叠加道中有效信号的匹配子波作为连续多道的中间道匹配追踪分解的子波库;最后,根据过中间道上每一点的同相轴倾角大小,反校正中间道子波库的子波中心时间,分别构建各地震道的子波库,约束匹配追踪分解的扫描范围,选择性地分离有效信号。

1.1 连续多道同相轴校正

叠后数据剖面中局部范围内的同相轴可认为是线性的,但其方向通常并非水平,如果直接叠加会造成波形畸变,因此在叠加前需进行水平校正。HUO等[21]提出了计算连续地震道同相轴倾角的方法。该方法根据倾角的大小,以中间道为基准,校正同相轴为水平后叠加地震道,叠加道可保持连续地震道的波形特征,具体如下。

选取2W1+1个连续地震道作为一组,在中间道上坐标为(i,W1+1)各点处定义长度为2W1+1,宽度为2W2+1、倾角为p的窗口。倾角改变对应着右侧采样点的变化,p=0时窗口水平,p>0时窗口相对p=0为逆时针旋转,反之,p<0时窗口相对p=0为顺时针旋转。设定p的取值范围为[pmin,pmax],间隔为一个采样点,依次计算不同倾角时窗口中连续地震道与中间道的距离,计算距离的函数定义如下。

式中:pi,k是过点(i,W1+1)的第k个倾角,当k=m时距离最小,可认为pi,m是经过该点的同相轴倾角;Xj(pi,k)表示在点(i,W1+1)处倾角为pi,k的窗口中第j道的数据向量,向量长度为2W2+1。如图1所示,中间道第57个采样点处的倾角为p57,m。采用该方法计算过中间道上每一个采样点的同相轴倾角,再依据倾角大小校正同相轴,得到同相轴水平的地震道集。如图2所示,图2a为合成的连续地震道;经校正后得到图2b;将图2a和图2b中的连续地震道叠加求平均,分别得到图2c和图2d中的蓝色虚线。对比可知:校正后的叠加道保持了波形特征,直接叠加地震道则会造成波形畸变。

图1 计算连续地震道同相轴倾角

图2 合成的连续地震道同相轴校正a 合成的连续地震道; b 校正后的连续地震道; c 中间道(红色实线)和未经校正的叠加道(蓝色虚线)对比; d 中间道(红色实线)和经校正后的叠加道(蓝色虚线)对比

1.2 稀疏反演分解提取匹配子波

地震褶积模型可以描述为:一个地震信号s(t)由地震子波w(t)和反射系数r(t)的褶积加上随机噪声[22]:

(3)

式中:“*”表示褶积运算。(3)式无法表示信号频率成分随时间的变化,考虑到地震信号包含有多种频率成分,可认为地震信号由不同频率的子波经地下界面反射叠加而成,(3)式可表示为:

(4)

式中:wi(t)为主频为下角标对应频率的子波;ri(t)是与wi(t)对应的频率依赖的伪反射系数;N为地震信号的频率最大值。将(4)式写成矩阵与向量乘积的形式,则地震褶积模型如下。

(5)

式中:SM×1是长度为M的地震信号向量;Wi(i=1,2,…,M×N)是长度为M的时频子波向量;Ri(i=1,2,…,M×N)为与之对应的系数;DM×(M×N)是M行、M×N列的矩阵,即超完备子波库;m(M×N)×1表示时频子波的系数向量;nM×1是噪声向量。根据地震信号傅里叶变换得到的振幅谱确定信号的频率范围,在该范围内以1Hz为间隔进行采样作为时频子波的主频,将这些子波沿时间轴平移建立超完备子波库。求解(5)式可得系数向量m(M×N)×1,将其重排成M×N的矩阵,即为地震信号的时频主频谱。时频主频谱平面上的每一个坐标点均对应着一个时频子波,所有坐标点的集合对应着超完备子波库。不考虑噪声的情况下(5)式可表示为(略去各个量的下角标):

(6)

(6)式是一个欠定方程,有无数个解,为了减小多解性,引入稀疏约束,即用尽量少的时频子波表示地震信号,从而得到高分辨率的时频主频谱。向量的稀疏度常用L0范数表示,但含L0范数的方程不易求解,将其转化为L1范数方程:

(7)

(7)式称为基追踪(basis pursuit,BP)问题。当噪声存在时,基追踪问题转变为基追踪降噪(basis pursuit de-noising,BPDN)问题:

(8)

式中:正参数σ是对数据中噪声水平的估计。当σ=0时,即为BP问题,因此BP问题是BPDN问题的特殊情况。本文运用L1范数最小化的谱投影梯度(spectral projected-gradient for L1 minimization,SPGL1)算法求解(8)式,该算法的基本思想是将BPDN问题转化为最小绝对值收敛和选择算法(least absolute shrinkage and selection operator,LASSO)问题,通过迭代求得符合BPDN条件的解。LASSO问题表示为:

(9)

式中:τ为m的L1范数约束。地震信号稀疏反演分解得到的时频主频谱平面上有能量的坐标点对应着构成地震信号的子波,称为匹配子波,因此通过稀疏反演可从超完备子波库中提取与地震信号匹配的子波。

1.3 基于稀疏反演的多道匹配追踪分解

地震信号匹配追踪分解受到噪声的强烈影响[16],结果具有非唯一性,本文提出的基于稀疏反演的多道匹配追踪算法能够减小噪声干扰,直接分离有效信号,具体实现步骤如下。

选出2W1+1个连续地震道{S1,…,SW1+1,…,S2W1+1}作为一组,经校正后叠加表示为:

(10)

(11)

(12)

(13)

因此,每次迭代更新后的残差信号为:

(14)

(15)

2 模型试算

2.1 含噪单道信号试算

图3a为合成的单道信号,由不同主频、时移和相位的Ricker子波组成,其子波主频从上到下依次为60,40,20,35,35,18,56Hz,加入高斯白噪声后如图3b所示。图3c为有效信号的理论时频主频谱,图中每一个具有能量的坐标点均代表着一个时频子波,表明信号中含有该时频子波成分,能量越强,与信号越匹配。建立超完备复Ricker子波库,设定σ值,图3d为合成的含噪信号稀疏反演分解得到的时频主频谱,与图3c对比发现,虽然能量分布范围有所扩大,但准确地定位出了时频主频谱的实际分布位置。稀疏反演分解方法具有较强的抗噪能力,能够精确地从超完备子波库中选取出与有效信号匹配的子波。

图3 合成的含噪单道信号时频主频谱a 合成的单道信号; b 加噪单道信号; c 有效信号的理论时频主频谱; d 合成的含噪信号稀疏反演分解的时频主频谱

图4a是根据图3b所示的信号单道匹配追踪分解得到的重构信号(红色虚线),由于分解过程中无法区分有效信号(绿色实线)和噪声,重构信号时引入了部分噪声。将理论时频主频谱(图3c)中有能量的坐标点对应的时频子波(有效信号的匹配子波)设定为子波库,采用与单道匹配追踪分解相同的迭代次数,得到图4b所示的重构信号(红色虚线),可以发现重构信号和有效信号基本重合,说明将地震信号中有效信号的最佳匹配子波作为子波库,可分离有效信号和噪声,准确重构有效信号。稀疏反演方法能较精确地选出有效信号的匹配子波,由此与匹配追踪分解结合也可实现有效信号的准确重构。

2.2 含噪连续多道信号试算

图5a为合成的地震数据剖面,包括两个斜线型和1个抛物线型同相轴;加入噪声后的剖面如图5b所示,一共有50道,分成10组,每组5道;水平校正每1组地震道后的剖面如图5c所示;图5d为所有组的叠加道,叠加道信噪比相较于中间道信噪比更高。第1组中间道的信噪比RSN=0.77,叠加道信噪比RSN=1.67(RSN=Jsignal/Jnoise,其中Jsignal为有效信号能量,Jnoise为噪声能量),高信噪比有利于稀疏反演提取出更准确匹配子波。每道有500个时间采样点,间隔为1ms,经傅里叶变换分析,所有道的频率范围为0～100Hz,以1Hz为主频率间隔采样,即有100个主频采样点,可建立矩阵大小为500×50000的超完备复Ricker子波库。图6a和图6b分别为采用本文方法去噪后的叠加剖面和重构误差剖面;图6c和图6d 分别为采用单道匹配追踪方法去噪后的叠加剖面和重构误差剖面,迭代次数均为20次。对比发现,两种方法均可有效压制数据中噪声,但采用本文方法得到的重构结果与原地震数据误差小,去噪效果佳。

图4 有效信号(绿色实线)的重构a 单道匹配追踪分解的重构信号(红色虚线); b 有效信号匹配子波作为子波库的匹配追踪分解的重构信号(红色虚线)

图5 合成的地震数据剖面及同相轴校正a 合成的地震数据剖面; b 加入噪声后的地震数据剖面; c 同相轴校正后剖面; d 所有组的叠加道(共10组)

图6 不同方法的去噪效果a 本文方法去噪后的叠加剖面; b 本文方法去噪后的重构误差剖面; c 单道匹配追踪分解去噪后的叠加剖面; d 单道匹配追踪分解的重构误差剖面

3 实际数据处理

图7a为某探区实际叠后数据剖面,共有50个地震道,每道400个采样点,采样间隔为4ms。采用本文方法对实际数据进行处理:设置7个连续地震道为1组,为了提高去噪后剖面的横向连续性,相邻两组共享3个地震道,除最后1组对该组内全部道进行分解外,其它组只分解前4道。采用单道匹配追踪方法和本文方法分解每一个地震道的迭代次数均为100次,图7b和图7c分别为采用本文方法去噪后得到的叠加剖面和残差剖面,处理结果表明,采用本文方法可较好地压制随机噪声,同时保持地震数据的横向连续性。这是由于子波库中的子波与地震道中有效信号是匹配的,在匹配追踪分解的过程中可有选择性地提取有效信号成分,避免了将噪声引入重构信号,同时,叠加道很好地保留了连续地震道中横向连续的有效信号成分,利用叠加道稀疏反演分解提取的子波建立子波库,匹配追踪分解识别了横向连续的有效信号,因此重构信号具有清晰的同相轴。图7d和图7e分别为采用单道匹配追踪方法去噪后得到的叠加剖面和残差剖面,相较于采用本文方法得到的处理结果,该方法得到结果中仍有部分噪声未去除,且同相轴不清晰。

图7 实际地震数据去噪效果a 实际叠后数据剖面; b 采用本文方法去噪后得到的叠加剖面; c 采用本文方法去噪后得到的残差剖面; d 单道匹配追踪分解去噪后得到的叠加剖面; e 单道匹配追踪分解去噪后得到的残差剖面

4 结论

本文提出的基于稀疏反演的多道匹配追踪去噪方法,利用叠加道的高信噪比、稀疏反演分解地震信号的抗噪性和解的稀疏表示特征,提取与有效信号匹配的子波作为子波库,在匹配追踪分解过程中分离有效信号和噪声,实现了信号的精确重构。拟合算例和实际地震数据处理结果表明,本文方法相较于单道匹配追踪方法具有以下优点:①去除噪声的效果佳;②去噪后数据剖面上同相轴横向连续性好;③合理地增加每一组中地震道数,运算效率高。但本文方法也存在着不足之处,当地震剖面中存在断层或者较多跳点时,连续地震道同相轴校正不准确,造成叠加道失真,构建的子波库中子波与有效信号不匹配,导致匹配追踪分解无法分离噪声和有效信号,这些问题有待于今后进一步的研究。