石油物探  2019, Vol. 58 Issue (2): 199-207  DOI: 10.3969/j.issn.1000-1441.2019.02.005
0
文章快速检索     高级检索

引用本文 

任浩, 李宗杰, 薛姣, 等. 基于稀疏反演的多道匹配追踪地震信号去噪方法及其应用[J]. 石油物探, 2019, 58(2): 199-207. DOI: 10.3969/j.issn.1000-1441.2019.02.005.
REN Hao, LI Zongjie, XUE Jiao, et al. Multichannel matching pursuit based on sparse inversion for seismic data denoising and its application[J]. Geophysical Prospecting for Petroleum, 2019, 58(2): 199-207. DOI: 10.3969/j.issn.1000-1441.2019.02.005.

基金项目

中石化科技部“顺北油田碳酸盐岩断裂体系及缝洞体刻画研究”项目(P17057-9)资助

作者简介

任浩(1990—), 男, 硕士在读, 主要从事地震信号时频分析方法研究。Email:1291009398@qq.com

通信作者

顾汉明(1963—), 男, 教授, 博士生导师, 主要从事油气地震勘探与开发研究。Email:hmgu@cug.edu.cn

文章历史

收稿日期:2018-06-07
改回日期:2018-11-28
基于稀疏反演的多道匹配追踪地震信号去噪方法及其应用
任浩1 , 李宗杰2 , 薛姣1 , 顾汉明1     
1. 中国地质大学(武汉)地球内部多尺度成像湖北省重点实验室, 湖北武汉 430074;
2. 中国石油化工股份有限公司西北油田分公司勘探开发研究院, 新疆乌鲁木齐 830011
摘要:基于超完备子波库的单道匹配追踪方法在自适应分解地震道时未能分离有效信号和噪声, 故难以直接、准确地去除噪声。为此提出了基于稀疏反演的多道匹配追踪去噪方法, 以多个连续地震道为研究对象, 首先采用多方向矢量中值滤波法计算同相轴倾角, 以中间道为标准将同相轴校正后再水平叠加; 然后采用基于基追踪降噪的稀疏反演方法分解叠加道, 选择与叠加道中有效信号匹配的子波作为中间道匹配追踪分解的子波库; 最后根据其它道相对中间道的校正量, 对中间道匹配追踪分解的子波库进行反校正, 分别构建各道匹配追踪分解的子波库。将基于有效信号匹配子波构建的子波库代替超完备子波库, 在搜寻最佳匹配子波时选择性地提取地震道中有效信号, 可直接分离有效信号和噪声。单道信号测试结果验证了该方法去除噪声的可行性; 合成和实际地震数据应用结果表明, 该方法相较于单道匹配追踪方法具有更好的去噪效果, 且去噪后剖面的横向连续性更好。
关键词匹配追踪    连续地震道    稀疏反演    叠加道    超完备子波库    去噪    横向连续性    
Multichannel matching pursuit based on sparse inversion for seismic data denoising and its application
REN Hao1, LI Zongjie2, XUE Jiao1, GU Hanming1     
1. Hubei Subsurface Multi-scale Imaging Key Laboratory, China University of Geosciences, Wuhan 430074, China;
2. Exploration and Development Research Institute of the Northwest Oilfield Branch Corporation, SINOPEC, Urumqi 830011, China
Foundation item: This research is financially supported by the Project from Department of Science and Technology of Sinopec (Grant No.P17057-9)
Abstract: Single-channel matching pursuit algorithm based on the over-complete wavelet library does not distinguish effective signals and noise in the process of adaptively decomposing seismic traces.Hence, it is difficult to remove noise directly and accurately.This paper, therefore, proposes a multi-channel matching pursuit denoising method based on sparse inversion.The method takes multiple continuous seismic traces as the research objects.Firstly, multi-directional vector median filtering method is used to calculate the inclination angle of seismic events.These events are corrected to be horizontal using the central trace as the standard, and then stacked.The sparse inversion method based on the basis pursuit denoising is then used to decompose the stack trace.Next, the wavelets matching the effective signal in the stack trace are selected as the wavelet library of the central trace matching pursuit decomposition.Finally, according to the correction amount of other traces relative to the central trace, the wavelet library of the central trace is inversely corrected to construct the wavelet library of each trace matching pursuit decomposition.The wavelet library constructed by matching wavelets of the effective signal replaces the over-complete wavelet library.In the process of searching for the best matching wavelet, the effective signal components in the seismic trace are selectively extracted, so that the effective signal and noise can be directly distinguished.The single-channel signal test result verified the effectiveness of the proposed method to remove noise.The results of the synthetic and real seismic data showed that the proposed method is better at denoising than the single-channel matching pursuit algorithm and can improve the lateral continuity of the seismic sections after denoising.
Keywords: matching pursuit    continuous seismic traces    sparse inversion    stack trace    over-complete wavelet library    denoising    lateral continuity    

地震数据不可避免地包含噪声, 噪声会影响到地震属性提取和地震解释等后续工作, 因此压制噪声是地震数据处理过程中的关键环节[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], 间隔为一个采样点, 依次计算不同倾角时窗口中连续地震道与中间道的距离, 计算距离的函数定义如下。

$ D\left[ {\mathit{\boldsymbol{X}}\left( {{p_{i,k}}} \right)} \right] = \sum\limits_{j = 1}^{2{W_1} + 1} {{{\left\| {{\mathit{\boldsymbol{X}}_{{W_{1 + 1}}}}\left( {{p_{i,k}}} \right) - {\mathit{\boldsymbol{X}}_j}\left( {{p_{i,k}}} \right)} \right\|}_2}} $ (1)
$ {p_{i,m}} = \mathop {\arg \min }\limits_{{p_{i,k}} \in \left[ {{p_{\min }},{p_{\max }}} \right]} D\left[ {\mathit{\boldsymbol{X}}\left( {{p_{i,k}}} \right)} \right] $ (2)

式中: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]:

$ s\left( t \right) = w\left( t \right) * r\left( t \right) + n\left( t \right) $ (3)

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

$ s\left( t \right) = \sum\limits_{i = 1}^N {\left[ {{w_i}\left( t \right) * {r_i}\left( t \right)} \right]} + n\left( t \right) $ (4)

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

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{S}}_{M \times 1}} = \left( {{\mathit{\boldsymbol{W}}_1},{\mathit{\boldsymbol{W}}_2}, \cdots ,{\mathit{\boldsymbol{W}}_i}, \cdots ,{\mathit{\boldsymbol{W}}_{M \times N}}} \right)\left( {\begin{array}{*{20}{c}} {{R_1}}\\ {{R_2}}\\ \vdots \\ {{R_i}}\\ \vdots \\ {{R_{M \times N}}} \end{array}} \right) + {\mathit{\boldsymbol{n}}_{M \times 1}} = }\\ {{\mathit{\boldsymbol{D}}_{M \times \left( {M \times N} \right)}}{\mathit{\boldsymbol{m}}_{\left( {M \times N} \right) \times 1}} + {\mathit{\boldsymbol{n}}_{M \times 1}}} \end{array} $ (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)式可表示为(略去各个量的下角标):

$ \mathit{\boldsymbol{S}} = \mathit{\boldsymbol{Dm}} $ (6)

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

$ \mathop {\min }\limits_m {\left\| \mathit{\boldsymbol{m}} \right\|_1}\;\;\;{\rm{s}}.\;{\rm{t}}.\;\;\;\;\mathit{\boldsymbol{S}} = \mathit{\boldsymbol{Dm}} $ (7)

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

$ \mathop {\min }\limits_m {\left\| \mathit{\boldsymbol{m}} \right\|_1}\;\;\;{\rm{s}}.\;{\rm{t}}.\;\;\;{\left\| {\mathit{\boldsymbol{Dm}} - \mathit{\boldsymbol{S}}} \right\|_2} \le \sigma $ (8)

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

$ \mathop {\min }\limits_m {\left\| {\mathit{\boldsymbol{Dm}} - \mathit{\boldsymbol{S}}} \right\|_2}\;\;{\rm{s}}.\;{\rm{t}}.\;\;\;{\left\| \mathit{\boldsymbol{m}} \right\|_1} \le \tau $ (9)

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

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

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

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

$ \mathit{\boldsymbol{\bar S}} = \mathit{\boldsymbol{Dm}} + \mathit{\boldsymbol{\bar n}} $ (10)

式中: $ \mathit{\boldsymbol{\bar S}}$为叠加道; D为超完备子波库; m为子波的系数向量; $ {\mathit{\boldsymbol{\bar n}}}$为残差值。建立超完备子波库, 采用SPGL1算法求解(10)式, 得到向量m的值, 选取向量中不为零的值对应的子波作为中间地震道匹配追踪分解(第W1+1道)的子波库DW1+1。反校正子波库DW1+1, 分别建立各地震道匹配追踪分解的子波库Dj(j=1, 2, …, 2W1+1)。

匹配追踪分解原理是基于最佳匹配原则, 不断迭代并从子波库中选取与残差信号最匹配的子波, 将原信号分解为一系列匹配子波的线性组合[14]。具体如下:对于组内某一地震道Sj, 将其匹配追踪分解的子波库Dj中任意一个时频子波记作Wγ, 初始输入的信号为Sj, 初始残差信号为Rj(0)=Sj, 每一次迭代提取一个与当前的残差信号最匹配的子波Wjn(以第n次迭代为例), 即:

$ \mathit{\boldsymbol{W}}_j^n = \mathop {\arg \max }\limits_{{\mathit{\boldsymbol{W}}_\gamma } \in {\mathit{\boldsymbol{D}}_j}} \frac{{\left| {\left\langle {\mathit{\boldsymbol{R}}_j^{\left( {n - 1} \right)},{\mathit{\boldsymbol{W}}_\gamma }} \right\rangle } \right|}}{{\left\| {{\mathit{\boldsymbol{W}}_\gamma }} \right\|}} $ (11)

式中:〈·, ·〉为两个信号的内积运算符号; ‖Wγ‖= $\sqrt {\left\langle {{\mathit{\boldsymbol{W}}_\gamma }, {\mathit{\boldsymbol{W}}_\gamma }} \right\rangle } $为子波信号的归一化算子; R(n-1)为第n-1次迭代的残差信号。

对于每一次选取的最佳匹配子波, 其对应的振幅系数ajn为:

$ a_j^n = \frac{{\left| {\left\langle {\mathit{\boldsymbol{R}}_j^{\left( {n - 1} \right)},\mathit{\boldsymbol{W}}_j^n} \right\rangle } \right|}}{{{{\left\| {\mathit{\boldsymbol{W}}_j^n} \right\|}^2}}} $ (12)

此次迭代分解的子成分信号为ajnWjn。用xjn表示每次迭代产生的子成分信号集合为:

$ \mathit{\boldsymbol{x}}_j^n = a_j^n\mathit{\boldsymbol{W}}_j^n,n = 1,2,3, \cdots $ (13)

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

$ \mathit{\boldsymbol{R}}_j^n = \mathit{\boldsymbol{R}}_j^{\left( {n - 1} \right)} - \mathit{\boldsymbol{x}}_j^n $ (14)

当残差信号Rjn的能量低于预设阈值或迭代次数达到最大迭代次数时终止迭代, 地震道Sj表示为:

$ {\mathit{\boldsymbol{S}}_j} = \sum\limits_{n = 1}^N {a_j^n\mathit{\boldsymbol{W}}_j^n} + \mathit{\boldsymbol{R}}_j^n $ (15)

式中:Wjn是与地震道有效信号匹配的子波, 是构成有效信号的子波成分。(15)式等号右边第1项不包含噪声, 第2项残差Rjn可认为是噪声, 如此便可区分有效信号和噪声。对每一组地震道均按以上步骤处理, 即可去除地震数据中的噪声。

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所示的重构信号(红色虚线), 可以发现重构信号和有效信号基本重合, 说明将地震信号中有效信号的最佳匹配子波作为子波库, 可分离有效信号和噪声, 准确重构有效信号。稀疏反演方法能较精确地选出有效信号的匹配子波, 由此与匹配追踪分解结合也可实现有效信号的准确重构。

图 4 有效信号(绿色实线)的重构 a单道匹配追踪分解的重构信号(红色虚线); b有效信号匹配子波作为子波库的匹配追踪分解的重构信号(红色虚线)
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次。对比发现, 两种方法均可有效压制数据中噪声, 但采用本文方法得到的重构结果与原地震数据误差小, 去噪效果佳。

图 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 结论

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

参考文献
[1]
刘霞, 黄阳, 黄敬, 等. 基于经验模态分解(EMD)的小波熵阈值地震信号去噪[J]. 吉林大学学报(地球科学版), 2016, 46(1): 262-269.
LIU X, HUANG Y, HUANG J, et al. Wavelet entropy threshold seismic signal denoising based on empirical mode decomposition(EMD)[J]. Journal of Jilin University(Earth Science Edition), 2016, 46(1): 262-269.
[2]
HORNBOSTEL S. Spatial prediction filtering in the t-xand t-xdomains[J]. Geophysics, 1991, 56(12): 2019-2026. DOI:10.1190/1.1443014
[3]
GACI S. The use of wavelet-based denoising techniques to enhance the first-arrival picking on seismic traces[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(8): 4558-4563. DOI:10.1109/TGRS.2013.2282422
[4]
HUANG N E, SHEN Z, LONG S R. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society of London, Series A:Mathematical, Physical and Engineering Sciences, 1998, 45(4): 903-995.
[5]
蔡建华, 王先春, 胡惟文. 基于经验模态分解与小波阈值的MT信号去噪方法[J]. 石油地球物理勘探, 2013, 48(2): 303-307.
CAI J H, WANG X C, HU W W. MT signal denoising method based on empirical mode decomposition and wavelet threshold[J]. Oil Geophysical Prospecting, 2013, 48(2): 303-307.
[6]
李月, 彭蛟龙, 马海涛, 等. 过渡内蕴模态函数对经验模态分解去噪结果的影响研究及改进算法[J]. 地球物理学报, 2013, 56(2): 626-634.
LI Y, PENG J L, MA H T, et al. Study of the influence of transition IMF on EMD do-noising and the improved algorithm[J]. Chinese Journal of Geophysics, 2013, 56(2): 626-634.
[7]
秦晅, 蔡建超, 刘少勇, 等. 基于经验模态分解互信息熵与同步压缩变换的微地震信号去噪方法研究[J]. 石油物探, 2017, 56(5): 658-666.
QIN X, CAI J C, LIU S Y, et al. Microseismic data denoising method based on EMD mutual information entropy and synchrosqueezing transform[J]. Geophysical Prospecting for Petroleum, 2017, 56(5): 658-666. DOI:10.3969/j.issn.1000-1441.2017.05.006
[8]
王姣, 李振春, 王德营. 基于CEEMD的地震数据小波阈值去噪方法研究[J]. 石油物探, 2014, 53(2): 164-172.
WANG J, LI Z C, WANG D Y. A method for wavelet threshold denoising of seismic data based on CEEMD[J]. Geophysical Prospecting for Petroleum, 2014, 53(2): 164-172. DOI:10.3969/j.issn.1000-1441.2014.02.006
[9]
DAUBECHIES I, LU J F, WU H T. Synchrosqueezed wavelet transforms:An empirical mode decomposition-like tool[J]. Applied and Computational Harmonic Analysis, 2011, 30(2): 243-261. DOI:10.1016/j.acha.2010.08.002
[10]
秦晅, 宋维琪. 基于同步压缩变换微地震弱信号提取方法研究[J]. 石油物探, 2016, 55(1): 60-66.
QIN X, SONG W Q. Weak signal extraction method of microseismic data based on synchrosqueezing transform[J]. Geophysical Prospecting for Petroleum, 2016, 55(1): 60-66. DOI:10.3969/j.issn.1000-1441.2016.01.008
[11]
SIAHSAR M A N, GHOLTASHI S, KAHOO A R, et al. Sparse time-frequency representation for seismic noise reduction using low-rank and sparse decomposition[J]. Geophysics, 2016, 81(2): V117-V124. DOI:10.1190/geo2015-0341.1
[12]
SHANG S, HAN L G, HU W. Seismic data analysis using synchrosqueezing wavelet transform[J]. Expanded Abstracts of 83rd Annual Internat SEG Mtg, 2013, 4330-4334.
[13]
DONOHO D L. De-noising by soft-thresholding[J]. IEEE Trans Inform Theory, 1995, 41(3): 613-627. DOI:10.1109/18.382009
[14]
MARLLAT S, ZHANG Z. Matching persuits with time-ferquency dictionaries[J]. IEEE Transactions on Antennas and Propagation, 1993, 41(12): 3397-3415.
[15]
赵天姿, 宋炜, 王尚旭. 基于匹配追踪算法的时频滤波去噪方法[J]. 石油物探, 2008, 47(4): 367-371.
ZHAO T Z, SONG W, WANG S X. Time-frequency filtering denoising method based on matching pursuit algorithm[J]. Geophysical Prospecting for Petroleum, 2008, 47(4): 367-371. DOI:10.3969/j.issn.1000-1441.2008.04.009
[16]
WANG Y H. Multichannel matching pursuit for seismic trace decomposition[J]. Geophysics, 2010, 75(4): V61-V66. DOI:10.1190/1.3462015
[17]
FENG X, ZHANG X B, LIU C, et al. Single-channel and multi-channel orthogonal matching pursuit for seismic trace decomposition[J]. Journal of Geophysics and Engineering, 2017, 14(1): 90-99. DOI:10.1088/1742-2140/14/1/90
[18]
PORTNIAGUINE O, CASTAGNA J. Inverse spectral decomposition[J]. Expanded Abstracts of 74th Annual Internat SEG Mtg, 2004, 1786-1789.
[19]
HAN L, HAN L G, LI Z. Inverse spectral decomposition with the SPGL1 algorithm[J]. Journal of Geophysics and Engineering, 2012, 9(4): 423-427. DOI:10.1088/1742-2132/9/4/423
[20]
张博伦, 王林飞, 蔡正辉, 等. SPGL1算法在地震资料随机噪声压制中的应用[J]. 中国海洋大学学报, 2016, 46(12): 81-86.
ZHANG B L, WANG L F, CAI Z H, et al. Application of SPGL1 in random noise suppression of seimic data[J]. Periodical of Ocean University of China, 2016, 46(12): 81-86.
[21]
HUO S, LUO Y, KELAMIS P G. Simultaneous sources separation via multidirectional vector-median filtering[J]. Geophysics, 2012, 77(4): V123-V131. DOI:10.1190/geo2011-0254.1
[22]
韩利, 刘春成, 张益明, 等. 地震复谱分解技术及其在烃类检测中的应用[J]. 地球物理学报, 2016, 59(3): 1095-1101.
HAN L, LIU C C, ZHANG Y M, et al. Seismic complexspectal decomposition and its application on hydrocarbon detection[J]. Chinese Journal of Geophysics, 2016, 59(3): 1095-1101.