一种基于块重组期望对数似然的地震勘探噪声压制方法与流程

本发明属于随机噪声的压制方法，具体涉及的是地震勘探数据处理领域，以地震勘探随机噪声的统计特征作为内部先验，通过分类重组控制块期望对数似然方法的高斯混合子模型匹配和参数选取，实现非平稳地震勘探随机噪声压制。

背景技术：

油气资源作为国家的战略资源，对国计民生有着重要的影响。地震勘探是油气资源的主要勘探手段，高质量的地震勘探记录对于地下地质结构成像、解释及对地震勘探数据进一步处理起到至关重要的作用。然而，地震勘探记录质量受强随机噪声的影响导致信噪比降低。在实际地震勘探测量中，需要在测区内布置多个检波器且记录时间较长，不同空间位置的检波器在同一时刻接收到的信号往往差别较大，甚至相邻记录道的数据变化规律也不尽相同，所采集的反射波、折射波等信号所经过的传递介质也不尽相同，获取到能够反映地下地质结构的信号微弱且复杂。同时，地震勘探数据也伴随着大量的强随机噪声，随机噪声强度及性质亦随接收时间和检波器空间位置变化而变化。强地震勘探随机噪声的存在阻碍甚至淹没有效信号，对地震勘探信号识别及进一步的处理产生不利的影响。压制低信噪比地震勘探记录中的随机噪声，恢复微弱的有效信号，对提高地震勘探记录质量和后续处理精度，实现油气资源的高精度地震勘探有着重要的基础意义。

针对地震勘探随机噪声压制以及微弱有效信号提取等一系列问题，地震勘探学者经过多年来的研究提出了多种去噪方法,包括时空反卷积、非线性扩散方法、基于小波变换的去噪方法、基于shearlet去噪方法、非局部均值等。这些方法取得较好的去噪效果，有的已经得到广泛应用。时空反卷积去噪方法是一种基于地震勘探同相轴线性预测的去噪方法，它对于线性同相轴有着很好的恢复效果，但是该方法不能直接应用于复杂分布的弯曲同相轴。非线性扩散方法将非线性扩散方程与能量扩散过程联系起来，该方法能够在压制地震勘探随机噪声的同时保持边缘细节信息。但是该方法在处理实际地震勘探数据时容易产生震荡，甚至增强噪声。小波变换方法是一种多尺度的几何分析方法，具有时频域上表达信号特征的能力，这一优点使得它在地震勘探去噪领域应用广泛。但它的问题在于没有充分地考虑地震勘探信号的方向性，因而它不能有效地表达信号的边缘细节。shearlet方法是对二维信号进行稀疏表示的一种方法，它弥补了小波变换方向性表达不足的问题，相对于小波变换，基于shearlet去噪方法能更好地保留信号中的细节部分。但是它对高频噪声抑制能力不足，容易受到伪吉布斯现象的影响。非局部均值(non-localmeans)方法是近年来发展起来的基于块的去噪方法，该方法利用地震勘探信号结构相似性以及随机噪声冗余性作为内部先验知识，利用非局部邻域数据的加权平均实现噪声压制。非局部均值方法在处理强噪声干扰下的地震勘探信号时，由于对相似结构信号赋予较大的权值，在滤波过程中对细节有较大程度的保留。但该方法仅利用信号相似结构这一内部先验，缺少外部先验知识的描述，在低信噪比情况下易产生信号失真，导致噪声残留，甚至会增强噪声。

期望块对数似然算法(expectedpatchloglikelihood)是一种基于先验知识的去噪方法，所谓的先验知识指的是待处理数据所属的某类数理统计模型，它是对该类数据的共性进行总结却又不失特性上的描述。基于期望块对数似然的去噪方法通常采用高斯混合模型，采用与待处理数据无关的训练样本，以最大期望算法(expectationmaximum)对高斯混合模型进行训练获得先验模型。利用训练好的高斯混合模型先验，期望块对数似然方法以迭代的方式，在最大似然概率前提条件下，对待处理数据块建模，结合局部的数据块信息和整体结构信息对待处理数据进行去噪处理。由于期望块对数似然方法在去噪过程中综合考虑局部先验和整体信号结构，因此既能保证局部建模的高效性，又能避免局部滤波的块效应问题，该方法能够在压制噪声的同时实现最优的信号保持效果。然而，期望块对数似然方法假设随机噪声为平稳噪声，当随机噪声分布在时间和空间上差异较大时，非平稳的噪声分布影响期望块对数似然方法的高斯混合子模型匹配以及滤波参数的准确性，继而导致强随机噪声难以去除的问题。

技术实现要素：

本发明提供一种基于块重组期望对数似然的地震勘探噪声压制方法，以解决强随机噪声难以去除的问题。

本发明采取的技术方案是，包括下列步骤：

s1.地震勘探数据分类重组，首先，对地震勘探数据分块处理，得到尺寸相同的地震勘探数据块；然后，根据地震勘探数据块的内部先验对地震勘探数据块分类，计算地震勘探数据块的统计分布，判定各地震勘探数据块统计分布的相似程度，将具有相似统计分布的地震勘探数据块判为同一类别，得到地震勘探数据类，最后，依次对地震勘探数据类扫描，排列同类的地震勘探数据块，将它们重组为地震勘探数据组；

s2.分别对待处理地震勘探数据组进行期望对数似然去噪，任取一地震勘探数据组,首先进行重叠分块，得到地震勘探数据重叠块，然后将预先训练好的、包含多个子模型的高斯混合模型作为先验模型，分别计算出每个重叠块在不同高斯子模型下的似然概率，以最大似然概率为选取准则，将具有最大似然概率的高斯子模型作为重叠块的先验，利用该先验构造滤波器对重叠块滤波，最后将滤波后的重叠块进行重构，得到去噪后的地震勘探数据组，重复上述步骤对地震勘探数据组迭代处理，直到到达预定的迭代次数，完成对当前待处理数据组的去噪，按上述步骤实现所有地震勘探数据组的去噪处理；

s3.地震勘探数据重构，首先按s1的扫描重组方法，从去噪后的地震数据组分离出地震勘探数据块，然后，将去噪后的地震勘探数据块按s1的分块矩阵重构，得到整个地震勘探数据的去噪结果。

本发明所述步骤s1地震勘探数据分类重组，具体表述为：

s11.输入的地震勘探观测数据为检波器记录的地下反射信号和环境噪声、人文噪声等地震勘探随机噪声的叠加，表示为y＝x+v，其中地震观测数据为y，x为有效地震勘探信号，v为随机噪声，对地震勘探数据y进行分块处理，子块分割算子为di，得到的第i个地震勘探数据块表示为：

yi＝diy＝xi+vi,i＝1,2,…,m。

其中m为地震勘探数据块数目；

s12.地震勘探数据块多阈值分类，首先计算地震勘探数据块yi的噪声标准差σi，求得所有块噪声标准差分布区间[σmin,σmax]，其中σmin为块噪声标准差最小值，σmax为块噪声标准差最大值，采用多阈值方法对地震勘探数据块分类：将块噪声标准差分布区间[σmin,σmax]等间隔分割成k个子区间，区间的间隔为第k类的判决阈值为tk＝σmin+q(k-1),k＝1,2,…,k+1。将地震勘探数据块yi归为第k类，表示为：

ci＝kiftk≤σi<tk+1,i＝1,2,…,m；k＝1,2,…,k,

其中k为分类数，ci为地震勘探数据块yi的分类索引；

s13.构建地震勘探数据组。分别对每类地震勘探数据按从左到右、从上到下的顺序依次扫描，将相同类别的地震勘探数据块重组，得到具有平稳分布随机噪声的地震勘探数据组，记第k类待处理地震勘探数据组：

yk＝{yi|ci＝k}

本发明所述步骤s2基于块期望对数似然方法以迭代的方式对处理地震勘探数据组进行去噪，具体表述如下：

s21.重叠分块，在第s次迭代中，对待处理数据组的去噪结果进行重叠分块，得到大小为l×l的地震勘探数据重叠块，第k类地震勘探数据组的第j个地震勘探数据块表示为：

其中rj为第j个地震勘探数据重叠块的分块算子；

s22.地震勘探数据重叠块匹配先验模型，计算地震勘探数据重叠块在第g个高斯混合子模型下的似然概率g＝1,2,...,g，求得具有最大对数似然概率的高斯混合子模型：

其中gmax为具有最大似然概率的高斯混合子模型索引；

s23.地震勘探数据重叠块滤波，利用gmax高斯混合的子模型作为先验，构造滤波器对地震勘探数据重叠块进行滤波，表示为:

其中为：

σgmax为gmax高斯混合子模型的协方差矩阵，分别为第k类待处理数据组的噪声方差；

s24.重构得到地震勘探数据组的滤波结果，对滤波后的地震勘探数据重叠块进行重构，得到地震勘探数据组的去噪结果：

其中λ和β均为正则化参数，重复步骤s21到s24对地震勘探数据组迭代处理，直到达到预定的迭代次数，得到待处理数据组yk的去噪结果k＝1,2,…,k。

本发明所述步骤s3地震勘探数据重构的具体表述如下：

s31.按s1的扫描重组方法，从去噪后的地震数据组中分离出地震勘探数据块，其中分离出的去噪后的地震勘探数据块i＝1,2,…,m；

s32.将去噪后的地震勘探数据块按s1的分块矩阵重构：

得到整个地震勘探数据的去噪结果

本发明提出一种块重组期望对数似然去噪算法并应用于地震勘探随机噪声压制，基于地震勘探数据的内部先验，提出利用噪声统计特征对数据重组分类，将噪声水平相近的数据子块聚集在一类。然后，以内部先验为依据，控制期望块对数似然的高斯混合子模型匹配、滤波器构造以及滤波参数的选取，进而实现对地震勘探数据的去噪。本发明将噪声统计特征(内部先验)与高斯混合模型(外部先验)相结合，提出一种基于块重组期望对数似然的地震勘探随机噪声压制方法，有效降低了非平稳随机噪声引起的数据块与高斯混合模型匹配精度的影响，拓展了块对数似然方法对非平稳分布随机噪声的处理能力，为具有非平稳分布随机噪声的地震勘探数据处理提供了一种有效、可靠的去噪方法。

本发明主要解决了如下技术问题：一、解决了地震勘探随机噪声分布非平稳导致的噪声残留问题；本发明步骤s1利用多阈值方法，基于噪声统计特征对待处理数据分类重组，得到噪声分布平稳的待处理数据组，达到了期望对数似然方法对噪声平稳分布的要求，有效地扩展了期望对数似然方法对地震勘探数据降噪的适用性。二、解决了低信噪比地震记录去噪过程中地震勘探信号衰减问题；本发明步骤s2利用训练好的高斯混合模型建立地震勘探数据统计模型，结合局部的数据块模型和整体结构信息对待处理数据进行去噪，能够在压制噪声的同时实现信号保持，为低信噪比地震记录去噪过程中地震勘探信号衰减问题提供有效解决方案。三、减少块效应，本发明步骤s3对去噪后的地震数据块重构过程中，结合整体信号保持与局部滤波的优势，避免了地震勘探数据分块处理的块效应问题。

本发明的优点：本发明提出基于块重组期望对数似然地震勘探数据去噪方法，利用训练好的高斯混合模型建立地震勘探数据统计模型，结合局部的数据块模型和整体结构信息对待处理数据进行去噪，能够在压制噪声的同时实现信号保持。此外，本发明以噪声统计特征作为内部先验，对含噪数据进行重组分类，使得重叠块与所匹配的先验之间拟合程度更高，避免地震勘探噪声非平稳性对子模型选取、滤波器构造以及正则参数选取的影响，这样得到的重构数据在噪声压制以及细节保持方面均有提升。本发明扩展了块对数似然方法对非平稳分布随机噪声的处理能力，尤其对具有非平稳分布随机噪声的地震勘探数据去噪效果明显，为低信噪比地震勘探数据处理提供了一种有效、可靠的去噪方法。

附图说明

图1是本发明基于块重组期望对数似然的地震勘探数据去噪方法的流程图；

图2a是不含噪声的人工合成地震数据图；

图2b是将图2a所示数据添加高斯白噪声的含噪数据图；

图3a是图2b所示数据记录噪声标准差分布情况图；

图3b是将图2b所示数据分类重组结果图；

图4a是采用现有的非局部均值技术对图2b所示数据的降噪结果图；

图4b是采用本发明方法对图2b所示数据的降噪结果图；

图5是实际共炮点地震勘探数据图；

图6是采用非局部均值方法对图5所示数据去噪的效果图；

图7是采用本发明方法对图5所示数据去噪的效果图。

具体实施方式

包括下列步骤：

s1.地震勘探数据分类重组，首先，对地震勘探数据分块处理，得到尺寸相同的地震勘探数据块。然后，根据地震勘探数据块的内部先验对地震勘探数据块分类。计算地震勘探数据块的统计分布，判定各地震勘探数据块统计分布的相似程度，将具有相似统计分布的地震勘探数据块判为同一类别，得到地震勘探数据类。最后，依次对地震勘探数据类扫描，排列同类的地震勘探数据块，将它们重组为地震勘探数据组；

得到的各个待处理地震勘探子块的噪声水平相近，为下一步匹配先验，构造滤波器打下了基础；

s2.分别对待处理地震勘探数据组进行期望对数似然去噪；任取一地震勘探数据组,首先进行重叠分块，得到地震勘探数据重叠块。然后将预先训练好的、包含多个子模型的高斯混合模型作为先验模型，分别计算出每个重叠块在不同高斯子模型下的似然概率，以最大似然概率为选取准则，将具有最大似然概率的高斯子模型作为重叠块的先验，利用该先验构造滤波器对重叠块滤波；最后将滤波后的重叠块进行重构，得到去噪后的地震勘探数据组；重复上述步骤对地震勘探数据组迭代处理，直到到达预定的迭代次数，完成对当前待处理数据组的去噪；按上述步骤实现所有地震勘探数据组的去噪处理；

s3.地震勘探数据重构，首先按s1的扫描重组方法，从去噪后的地震数据组分离出地震勘探数据块，然后，将去噪后的地震勘探数据块按s1的分块矩阵重构，得到整个地震勘探数据的去噪结果；

所述步骤s1基于内部先验对地震勘探数据分类重组，其具体表述为：

s11.输入的地震勘探观测数据为检波器记录的地下反射信号和和环境噪声、人文噪声等地震勘探随机噪声的叠加，表示为y＝x+v，其中地震观测数据为y，x为有效地震勘探信号，v为随机噪声；由于地震勘探检波器分布距离远，由风吹草动、人机运动引起的地震勘探随机噪声在各道地震勘探记录中差异较大，噪声的统计分布在时间、空间上不平稳，而现有滤波技术通常假设地震勘探随机噪声平稳，导致随机噪声压制不彻底，期望块对数似然滤波方法采用与噪声相关的全局参数进行滤波，在处理这种非平稳地震勘探随机噪声时，由于地震勘探数据的噪声非平稳主要表现为噪声强度在时间和空间上差异，在采用基于块的滤波方法时，这种差异造成块的内部先验差异较大，在块先验模型选取时产生误差，造成滤波结果中的局部噪声残留问题，因此，需要在块期望对数似然滤波过程中采用不同的滤波参数来选择先验模型才能获得好的滤波效果，本发明基于地震勘探数据的噪声统计量作为内部先验，对地震勘探数据进行分类重组，使获得的地震勘探数据组的随机噪声平稳，本发明对地震勘探数据y进行分块处理，为避免块效应的产生，每个地震勘探数据块边界处均有重叠部分，分块算子记为di，得到的第i个地震勘探数据块表示为：

yi＝diy＝xi+vi,i＝1,2,…,m

其中m为地震勘探数据块数目；

s12.基于内部先验的地震勘探数据块多阈值分类，计算地震勘探数据块yi的噪声标准差σi，求得所有块噪声标准差分布区间[σmin,σmax]，其中σmin为块噪声标准差最小值，σmax为块噪声标准差最大值。采用多阈值方法对地震勘探数据块分类：将块噪声标准差分布区间[σmin,σmax]等间隔分割成k个子区间，区间的间隔为第k类的判决阈值为：

tk＝σmin+q(k-1),k＝1,2,…,k+1

其中k为分类数；

将地震勘探数据块yi归为第k类，表示为：

ci＝kiftk≤σi<tk+1,i＝1,2,…,m；k＝1,2,…,k，

其中ci为地震勘探数据块yi的分类索引。

本发明通过对地震勘探数据分类重组，使得每类中的块内部先验-噪声标准差相似，提高地震勘探数据块匹配高斯混合子模型精度，从而获得好的随机噪声压制效果。

s13.构建地震勘探数据组，分别对每类地震勘探数据按从左到右、从上到下的顺序依次扫描，将相同类别的地震勘探数据块重组，得到具有平稳分布随机噪声的地震勘探数据组。记第k类待处理地震勘探数据组：

yk＝{yi|ci＝k}

待处理地震勘探数据组中的随机噪声分布平稳，降低了地震勘探随机噪声非平稳性对去噪结果所带来的影响；

所述步骤s2基于块期望对数似然方法以迭代的方式对处理地震勘探数据组进行去噪，具体表述为：

s21.重叠分块，在第s次迭代中，对待处理数据组的去噪结果进行重叠分块，得到大小为l×l的地震勘探数据重叠块，第k类地震勘探数据组的第j个地震勘探数据块表示为：

其中rj为第j个地震勘探数据重叠块的分块算子；

s22.地震勘探数据重叠块匹配先验模型，计算地震勘探数据重叠块在第g个高斯混合子模型下的似然概率g＝1,2,...,g，n(·)为高斯子模型概率密度函数,g为高斯混合子模型数，σk是第k组噪声标准差，计算得到具有最大对数似然概率的高斯混合子模型：

其中gmax为具有最大似然概率的高斯混合子模型索引；

可以看出对重叠块建模的高斯子模型受噪声标准差这一内部先验影响较大。由于重叠分块数量非常大，块期望对数似然滤波方法采用整体数据的标准差代替重叠块方差进行子模型匹配，在处理非平稳随机噪声时，重叠块噪声标准差与整体数据标准差的差异会在子模型匹配过程中产生偏差，造成强噪声压制不彻底，因此，本发明提出基于块重组期望对数似然去噪方法，通过对地震勘探数据重组，使得每组数据的随机噪声平稳，此时采用数据组的噪声标准差建模有助于提高重叠块与先验子模型的匹配度，这对于处理这种非平稳随机噪声是非常必要的；

s23..地震勘探数据重叠块滤波，利用gmax高斯混合的子模型作为先验，构造滤波器对地震勘探数据重叠块进行滤波，表示为:

其中为：

σgmax为gmax高斯混合子模型的协方差矩阵，分别为第k类待处理数据组的噪声方差，在获取期望重叠块统计先验的情况下，减弱输入重叠块与随机噪声的相关性，使得到的结果向期望重叠块逼近；

s24.重构得到地震勘探数据组的滤波结果，对滤波后的地震勘探数据重叠块进行重构，得到地震勘探数据组的去噪结果：

其中λ和β均为正则化参数，用于平衡全局滤波与局部的关系，λ控制细节保留，β则控制噪声压制，它们的取值与当前处理的地震勘探数据组噪声方差有关，在处理地震勘探数据时，由于随机噪声分布非平稳，本发明根据各个地震勘探数据组的噪声水平的不同设定不同的λ和β，从而在处理地震勘探数据时能够有效压制地震勘探随机噪声且不损伤细节；

重复步骤s21到s24对地震勘探数据组迭代处理，直到达到预定的迭代次数，得到待处理数据组yk的去噪结果k＝1,2,…,k。

所述步骤s3地震勘探数据重构的具体表述为：

s31.按s1的扫描重组方法，从去噪后的地震数据组中分离出地震勘探数据块，其中分离出的去噪后的地震勘探数据块i＝1,2,…,m；

s32.将去噪后的地震勘探数据块按s1的分块矩阵重构：

得到整个地震勘探数据的去噪结果

本发明结合局部去噪和整体重构，对所有滤波的地震勘探数据组重构，从而有效解决基于块的滤波方法的块效应问题。

为验证本发明方法的有效性，将本发明分别应用于合成地震勘探记录和实际共炮点地震勘探记录。所采用的合成地震勘探记录由雷克子波生成，如图2a所示。其中3条同相轴主频分别为100hz、90hz、70hz，传播速度分别为6000m/s、5000m/s、3000m/s，雷克子波幅值均为1，其抽样频率为1000hz。在合成记录中加入高斯白噪声，且在11-20、21-30、41-50地震道分别添加标准差为0.245、0.196、0.294的随机噪声以仿真实际地震勘探随机噪声，得到的含噪合成记录的平均信噪比为-3.06db，如图2b所示。

首先按本发明s11将含噪合成地震勘探记录分块处理，估计地震勘探数据块的噪声标准差：

s12.统计噪声标准差的分布情况(图3a所示)，可知标准差分布不均匀。计算得到块噪声标准差分布区间为[0,0.35]。设子区间的间隔为q＝0.1，对地震勘探数据块分类。

s13.将同类地震勘探数据块依次组合在一起，完成数据的分类重组，如图3b所示，可见每个待处理数据组噪声分布都是均匀的。然后按s2和s3实现对合成地震勘探数据的块重组期望对数似然去噪，并与非局部均值方法的去噪结果进行了比较。附图4a和附图4b分别为非局部均值方法和本发明的去噪结果。经比较可见，两种方法均有较好的去噪效果，大部分地震勘探背景噪声得到压制。但非局部均值处理结果中，在噪声较强区域(红框区域)仍有噪声残留。与之相比，本发明方法不仅有效地去除了红框区域内的强随机噪声，同相轴也被更好地恢复，连续性好。

在信噪比提高方面，我们改变合成记录的噪声强度，分别采用非局部均值和本发明方法处理。经非局部均值方法和本发明方法处理后的地震勘探记录信噪比记录于表1。由表1可知本方面方法在信噪比提高方面明显优于非局部均值方法。

表1

为验证本发明的应用效果，将本发明应用于实际共炮点地震勘探记录。所采用的实际地震勘探记录为中国某矿区的36道共炮点记录，如附图5所示。该记录采样间隔为1ms，道间距为30m。该记录中存在大量随机噪声，且在第3、10、25和30道区域噪声水平明显增强。强噪声的存在掩盖了有效信号，部分区域甚至无法识别同相轴。

分别采用非局部均值方法和本发明提出方法对共炮点记录进行处理，得到的结果如附图6和附图7所示。从非局部均值处理结果可见，大部分地震勘探随机噪声被压制，但由于地震勘探随机噪声的非平稳性影响滤波效果，在非局部均值的处理结果中，得到的信号过于平滑且产生失真。相对于非局部均值方法，本发明的处理结果在压制随机噪声的同时，能够很好地保持有效信号的幅度和结构，同相轴得到更好地恢复。

合成地震勘探记录和共炮点地震勘探记录的处理结果表明，本发明以噪声统计特征作为内部先验进行块分类重组，能够将具有相同噪声水平的地震勘探数据块归为一类，分别采用块期望对数似然方法处理地震勘探数据组，这样就避免了由非平稳噪声引起的数据块建模误差，进而提高了块期望对数似然去噪方法重构数据的准确性，达到更好的噪声压制效果。本发明提出的块重组期望对数似然去噪方法实现了基于内部先验自适应调节数据块建模参数，能够有效压制非平稳地震勘探随机噪声，更好地恢复地震勘探信号。与非局部均值方法相比，本发明在强随机噪声压制、信号保持和信噪比提高方面更具优势。本发明可广泛地应用于复杂分布的地震勘探随机噪声压制。

本发明内容并不局限于以上实例，熟悉本领域的技术人员可以在本发明提出的框架之下，进行等同变换或者变形，这些均包含在本发明的权利要求之内。