一种基于形态成分分析的采集脚印压制方法

摘要

本发明公开一种基于形态成分分析的采集脚印压制方法，包括：步骤101：根据三维地震数据体的沿层切片或者等时切片的地震记录中采集脚印波形形态特征构造二维局部离散余弦变换，并与二维平稳小波变换联合构成超完备字典；步骤102：通过对地震记录信号的原始数据逐层执行的方式，在每个沿层切片或者等时切片内利用基于形态成分分析的方法进行采集脚印的初步压制；步骤103：确定二维局部离散余弦变换的低截频；步骤104：重复步骤102‑103直到所有切片数据处理完成，利用基于形态成分分析的方法进行信噪分离，压制采集脚印噪声，最终实现对三维地震数据体中采集脚印噪声的压制。利用本发明对三维地震数据体进行采集脚印压制，达到了提高了地震资料信噪比的目的。

说明书

技术领域

本发明属于地震勘探中数据处理技术领域，特别涉及一种地震采集数据中采集脚印压制的方法。

背景技术

采集脚印一般是由于地震观测系统的滚动排列方式以及震源和检波器侧线间隔造成的不完全采样而引起地震成像中产生的一种周期性振幅假象，通常可以在三维地震数据的时间或者深度切片中观察到。在地层-岩性等隐蔽型油气藏的勘探过程中，为了满足复杂储层预测的需求，需要对地震资料进行精细的结构解释和岩性圈闭解释。然而，采集脚印噪声会严重干扰地震解释过程，甚至会被误判为地质结构或者岩性异常，严重影响地震资料解释的精度和可靠性。因此研究压制和衰减采集脚印的方法十分重要，具有重大的理论意义与市场价值。

随着信号稀疏性理论的发展，Starck等人提出了形态成分分析的混合信号分解方法。形态成分分析是指，根据复杂信号的组成成分波形形态特征，将两种具有不同原子特征的变换字典构成超完备字典，实现对复杂信号更稀疏的表示方式以及更有效的信息识别能力，实现两种成分的分离。字典通常是根据经验从已知的数学变换中选择或者构造，而选出的字典是否充分满足形态成分分析的假设，是形态成分分析方法能否成功的关键。使用形态成分分析的方法对信号与谐波噪声实现分离来达到压制谐波噪声的目的，需要根据有效信号与谐波噪声的波形形态特征构造合适的稀疏变换。已有文献研究选择使用连续小波变换来稀疏表示有效信号比较合适。但谐波噪声较为复杂，需要构造较为合适匹配的变换来稀疏表示，并选择确定变换的具体参数，目前相关研究较少。

现有技术：

优化采集方法。采用宽方位角的观测系统，可降低数据记录方位角分布的不均匀性，在地面进行全方位观测，使得地震资料记录了地下反射点各个方向上的波场特征，获得完整的地震波场信息；宽方位角的观测系统还有助于减少炮点的反向散射噪声，改善了静校正效果，这些都有助于减弱采集脚印。

现有技术的缺点：

只能针对新采集数据，无法处理之前采集的数据；需要针对不同采集系统做不同的设计，可移植性差。

发明内容

本发明的目的在于提供了一种基于形态成分分析的采集脚印压制方法，以解决上述技术问题。针对三维地震数据体沿层切片或者等时切片内的采集脚印和有效波信号的波形形态特征不同，构造合适的稀疏表示变换，即选取二维局部离散余弦变换表示采集脚印，选取二维平稳小波变换表示有效波信号，并将这两类波形形态字典联合构成超完备字典，并根据数据特征选择确定二维平稳小波变换变换参数，以及使用坐标分块松弛算法方法实现信噪分离以达到有效压制三维地震数据体的沿层切片或者等时切片的地震数据中的采集脚印的目的。

为了实现上述目的，本发明采用如下技术方案：

基于形态成分分析的采集脚印压制方法，包括以下步骤：

步骤101：根据三维地震数据体的沿层切片或者等时切片的地震记录中采集脚印波形形态特征构造二维局部离散余弦变换，并与二维平稳小波变换联合构成超完备字典，同时确定二维平稳小波变换的分解级数J和二维局部离散余弦变换的分析窗尺寸W；

步骤102：通过对地震记录信号的原始数据逐层执行的方式，在每个沿层切片或者等时切片内利用基于形态成分分析的方法进行采集脚印的初步压制；

步骤103：确定二维局部离散余弦变换的低截频；

步骤104：重复步骤02-03直到所有切片数据处理完成，利用基于形态成分分析的方法进行信噪分离，压制采集脚印噪声，最终实现对三维地震数据体中采集脚印噪声的压制。

进一步的，步骤101中根据有效信号波形形态特征与采集脚印波形形态特征确定形态成分分析的两种变换字典，有效信号选择二维平稳小波变换，采集脚印选择局部离散余弦变换，并构成超完备字典。

应用形态成分分析的对象是含有两种具有不同形态特征的成分：

s＝s₁+s₂，

式中：s表示待分析信号；s₁、s₂表示信号中两种具有不同形态特征的成分；分别提取出s₁、s₂这两种成分是形态成分分析的目标；假设s₁和s₂可以分别由字典Φ₁和Φ₂有效的稀疏表示，但是用Φ₂稀疏表示s₁和用Φ₁稀疏表示s₂时稀疏性差；

选择二维平稳小波变换作为稀疏表示有效信号成分的变换字典，其中二维平稳小波正变换为：

式中H_j和G_j分别代表第j层分解的滤波器组。

二维平稳小波变换的反变换为：

式中和分别代表H_j和G_j的对偶滤波器组。

选择二维局部离散余弦变换(IV型)作为稀疏表示采集脚印的变换字典，其中二维局部离散余弦变换的正变换为：

式中，f(i,j)表示待分析信号，表示待分析信号的二维局部离散余弦变换系数，k₁,k₂＝0,1,...N-1。

二维局部离散余弦变换的反变换为：

根据形态成分分析理论，用上面选定的字典Φ₁即二维平稳小波变换和Φ₂即二维局部离散余弦变换，构成超完备字典，稀疏表示信号s，计算稀疏表示系数：

式中：x₁为重构系数中与Φ₁对应的部分；x₂为重构系数中与Φ₂对应的部分。为拉格朗日乘子。

进一步的，步骤102中使用分块坐标松弛算法对对原始地震记录数据进行沿层切片或者等时切片处理，实现采集脚印的压制，该优化问题可以通过坐标分块松弛算法求解。坐标分块松弛算法的基本思想是交替迭代的计算x₁和x₂。其主要内容步骤为：

初始化：初始迭代步数k＝0，初始解

用来表示信号成分s₁的系数初始解，用来表示信号成分s₂的系数初始解；

迭代：每步迭代时迭代步数k增加1，并计算：

式中，T_λ为硬阈值函数；与Φ₁构成一对正反变换，与Φ₂构成一对正反变换；

终止条件：当小于预设的值时，继续迭代对结果的影响足够小时，迭代终止；

输出：

为分离的信号成分s₁的最终变换系数，为分离的信号成分s₂的最终变换系数；

分块坐标松弛算法中，硬阈值函数公式如下：

式中：为硬阈值函数，λ为硬阈值，为系数矩阵的元素，k＝1,2,...,N，N为系数矩阵的尺寸。

进一步的，步骤103即确定二维局部离散余弦变换的低截频实现中，具体包括：

步骤301：按照间隔时间ΔT对三维地震数据体取M个等时切片用于试验局部余弦变换的重构变换的低截频参数，它们的位置分别为t₀,t_ΔT,t_2ΔT,…,t_(M-1)ΔT；

步骤302：对于位置t_kΔT,k＝0,1,...,M-1，给出二维局部离散余弦变换的重构变换的一系列低截频参数，并针对每一个低截频参数，求解获得有效波信号和采集脚印噪声的分离；通过对比各个参数下的信噪分离效果确定在当前时间位置取得最佳信噪分离效果的低截频参数f_kΔT；

步骤303：利用线性插值方法，根据时间位置t_kΔT，k＝0,1,...,M-1以及其对应的低截频参数f_kΔT，k＝0,1,...,M-1，得到整个三维地震数据体其它时间位置的二维局部离散余弦变换的重构变换的低截频参数；

进一步的，步骤04中按照从浅层到深层逐层处理三维地震数据体的各个时间切片，通过求解式

得到有效波信号和采集脚印噪声的表示向量x_e和x_f，从而得到从当前时间切片去除的采集脚印噪声s_f＝Φ_fx_f以及对应的有效波信号为s_e＝s-s_f。

进一步的，预设的值取10^-6。

进一步的，硬阈值λ取中按照从大到小排列的第个值。

相对于现有技术，本发明具有以下有益效果：

本发明利用形态成分分析理论，将采集脚印噪声也看作叠后地震波场中的一种信号分量，构造合适的稀疏表示变换即二维局部离散余弦变换，与二维平稳小波变换构成超完备字典，并根据数据特征选择确定二维平稳小波变换参数，以及使用坐标分块松弛算法方法实现信噪分离达到压制谐波噪声的目的。本发明不仅能够有效的压制采集脚印，去除部分随机噪声和偏移处理假象，而且有效信号具有较高保真性。

附图说明

图1为局部余弦基时域波形图；

图2为钟形窗函数(k＝0)；

图3为二维8×8离散余弦变换原子；

图4为二维平稳小波变换原子；

图5为本发明方法处理数据的流程图；

图6为确定局部余弦变换重构变换的低截频参数的流程图；

图7A为实际数据一条主测线剖面；图7B为本发明方法从图7A中分离的有效波信号；图7C为本发明方法从图7A中去除的采集脚印噪声；

图8A为实际数据体中1.7s的时间切片；图8B本发明方法从图8A中分离的有效波信号；图8C本发明方法从图8A中去除的采集脚印噪声。

具体实施方式

为了使本发明的目的、技术方案更加清楚明白，下面结合附图和具体实施方式对本发明做进一步详细的说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为本发明的限定。

在本发明实施例中，提出了一种基于形态成分分析的采集脚印压制方法，包括以下步骤：

步骤101：根据三维地震数据体的沿层切片或者等时切片的地震记录中采集脚印波形形态特征构造二维局部离散余弦变换，并与二维平稳小波变换联合构成超完备字典，同时确定二维平稳小波变换的分解级数J和二维局部离散余弦变换的分析窗尺寸W；

步骤102：通过对三维地震数据体的原始数据逐层执行的方式，在每个沿层切片或者等时切片内利用基于形态成分分析的方法进行采集脚印的初步压制；

步骤103：确定二维局部离散余弦变换的低截频；

步骤104：重复步骤102-103直到所有切片数据处理完成，利用基于形态成分分析的方法进行信噪分离，压制采集脚印噪声，最终实现对三维地震数据体中采集脚印噪声的压制。

步骤01中根据三维地震记录中采集脚印波形形态特征构造二维局部离散余弦变换，并与二维平稳小波变换构成超完备字典，具体包括：

选择二维局部离散余弦变换作为稀疏表示采集脚印的变换字典，离散余弦变换中常用的局部余弦基有I型和IV两种，本发明主要采用IV型余弦基，局部IV型余弦基函数(LCB-IV)的定义如下：

上式中，b_mn(x)表示具有波数指标m和位置指标n的局部余弦基，如图1所示，其中为区间起始位置，为区间终止位置，为区间长度。钟形窗函数B_n(x)是定义在闭区间的光滑函数，且ε和ε′分别为区间左、右两侧的重叠半径，

其中β(x)为轮廓函数(递增或递减变换)，其表达式通过递归定义为

这里k≥0，且

轮廓函数β_k+1(x)的光滑程度将随k的增大而增加，例如k＝0时的轮廓函数为

如图2所示为此时对应的钟形窗函数。

对于函数f∈C^N，离散IV型余弦变换(DCT-IV)定义如下:

而其逆变换(IDCT-IV)定义为:

选择二维平稳小波变换作为稀疏表示有效信号的变换字典，对于给定一维的尺度函数φ(t)与小波函数ψ(t)，二维平稳小波变换的尺度函数和小波函数通过φ(t)和ψ(t)以张量积的方式得到：

二维尺度函数φ(x,y)：φ(x,y)＝φ(x)φ(y)；

二维水平方向小波函数Ψ^H(x,y)：ψ^H(x,y)＝φ(x)ψ(y)；

二维垂直方向小波函数ψ^V(x,y)：ψ^V(x,y)＝ψ(x)φ(y)；

二维对角方向小波函数ψ^D(x,y)：ψ^D(x,y)＝ψ(x)ψ(y)。

二维平稳小波变换将第j层的信号低频部分分解为第j+1层的低频部分和垂直、水平和对角方向的高频部分，其中信号的低频部分对应行为低频、列为低频的信号；信号的水平高频部分对应行为低频、列为高频的信号；信号的垂直高频部分对应行为高频、列为低频的信号；信号的对角高频部分对应行为高频、列为高频的信号。利用多孔算法来实现平稳小波变换，定义滤波器组H和G，那么H_j和G_j分别代表第j层分解的滤波器组，通过对H和G的各个系数之间插入2^j-1个零得到，对于任何j≥0，第j层的平稳小波变换的正变换为：

如果H_j和G_j的对偶滤波器组分别为和那么可以得到第j层的平稳小波变换的反变换为：

如图3与图4所示，分别为二维8×8的离散余弦变换原子和二维平稳小波变换原子,其中，图3中的二维离散余弦变换原子从左边到右边其水平方向的频率逐渐升高，而从上边到下边其垂直方向的频率逐渐升高；图4为三个不同分解尺度的二维平稳小波变换原子，在每个尺度内均有水平、垂直和对角三个方向，完全符合形态成分分析理论的要求，取得预期的对不同信号分量的稀疏表示和期望分离效果。

根据形态成分分析理论，用上面选定的字典Φ₁即二维平稳小波变换和Φ₂即二维局部离散余弦变换构成超完备字典稀疏表示信号s，计算稀疏表示系数：

式中：x₁为重构系数中与Φ₁对应的部分；x₂为重构系数中与Φ₂对应的部分。为拉格朗日乘子。该优化问题通过坐标分块松弛算法求解。

步骤02中使用分块坐标松弛算法对对原始地震记录数据进行沿层切片或者等时切片处理，实现采集脚印的压制，该优化问题通过坐标分块松弛算法求解。坐标分块松弛算法的基本思想是交替迭代的计算x₁和x₂。其主要内容步骤为：

初始化：初始迭代步数k＝0，初始解

用来表示信号成分1的系数初始解，用来表示信号成分2的系数初始解；

迭代：每步迭代时迭代步数k增加1，并计算：

式中，T_λ为硬阈值函数；与Φ₁构成一对正反变换，与Φ₂构成一对正反变换；

终止条件：当小于预设的值时，继续迭代对结果的影响足够小时，迭代终止；

输出：

为分离的信号成分1的最终变换系数，为分离的信号成分2的最终变换系数；

分块坐标松弛算法中，硬阈值函数公式如下：

式中：为硬阈值函数，λ为硬阈值，通常取中按照从大到小排列的第个值，为系数矩阵的元素，k＝1,2,...,N，N为系数矩阵的尺寸，Φ为变换字典，s为时域信号。

如图6所示，步骤103中确定二维局部离散余弦变换的重构变换的低截频，主要包括如下步骤：

步骤301：按照间隔时间ΔT对三维地震数据体取M个等时切片用于试验局部余弦变换的重构变换的低截频参数，它们的位置分别为t₀,t_ΔT,t_2ΔT,…,t_(M-1)ΔT；

步骤302：对于位置t_kΔT,k＝0,1,...,M-1，给出二维局部离散余弦变换的重构变换的一系列低截频参数，并针对每一个低截频参数，求解上述最优化问题获得有效波信号和采集脚印噪声的分离；通过对比各个参数下的信噪分离效果确定在当前时间位置取得最佳信噪分离效果的低截频参数；

步骤303：利用线性插值方法，根据时间位置t_kΔT，k＝0,1,...,M-1以及其对应的低截频参数f_kΔT，k＝0,1,...,M-1，得到整个三维地震数据体其它时间位置的二维局部离散余弦变换的重构变换的低截频参数。

在本实施例中，提出了一种基于形态成分分析的采集脚印压制方法，从而达到了有效压制三维地震数据中采集脚印目的。本发明具有以下有益效果：

利用本发明的方法对三维地震数据体进行采集脚印压制，不仅能够有效的压制采集脚印，去除部分随机噪声和偏移处理假象，而且有效信号具有较高保真性。

下面结合一个具体实施例对上述变换构造及参数确定方法进行具体说明，且该实施例仅是为了更好说明本发明，并不构成对本发明的不当限定。

图7A为某油田的海上三维地震资料的一条主测线剖面，该资料的采集时间较早，由于受到当时采集技术和采集仪器的限制，不同测线之间的采集质量差异较大，远缆采集质量不好，使得偏移数据受采集脚印的影响非常严重；另外，在海上地震采集过程，往往由于季节更替和洋流的影响，海水的温度和盐度都会有较大的变化，容易造成地震波的传播速度发生改变从而导致海上地震资料中产生明显的采集脚印。

如图7B和图7C所示，通过本发明方法对该含有采集脚印的地震记录(图7A)进行处理，得到对应的有效信号和去除的采集脚印噪声。本发明方法可以有效的对采集脚印噪声进行压制，还去除了部分随机噪声和偏移处理假象，处理后的资料变得更加清晰、同相轴的连续性更强，从而印证了本发明方法有较好的有效波信号保真度。

图8A该三维地震数据体中1.7s的时间切片，可以看出，该三维数据体中采集脚印主要沿主测线方向分布，并且采集脚印现象非常明显，导致很难从原始时间切片中确定这种振幅异常是由地下介质横向变化造成还是采集脚印的影响。

如图8B-图8C所示，应用本发明的方法对该三维数据体进行处理的结果，本发明方法很好地压制了脚印噪声，且较好地保持了剖面内的有效波信号的各种细微结构特征。

以上的模型和实际资料算例中，利用本发明的方法构造采集脚印的稀疏表示变换，对地震资料数据进行采集脚印压制，不仅能够有效的压制采集脚印噪声，而且有效信号具有较高保真性，为后续资料的分析奠定基础。

最后需要说明的是，以上模型和实际资料算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围。凡在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。