CN105607125A - 基于块匹配算法和奇异值分解的地震资料噪声压制方法 - Google Patents

Abstract

本发明公开了一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，解决现有方法由于没有考虑到地震资料同相轴固有的局部和非局部相似性，造成的噪声压制不充分和有效信息衰减的问题，该方法包括相似性分组、三维变换奇异值滤波去噪、根据基础估计重新分组并计算维纳收缩系数和维纳滤波四个步骤。该方法在理论上采用块匹配算法，考虑到地震资料同相轴的局部和非局部相似性的固有特点，采用二维奇异值分解和一维奇异值分解结合，充分利用地震资料的结构信息，实现地震资料的稀疏表示，能够较大程度地减小有效信号衰减，并提高噪声压制效果。

Description

基于块匹配算法和奇异值分解的地震资料噪声压制方法

技术领域

本发明属于地震资料处理领域，具体涉及一种基于块匹配算法和奇异值分解的地震资料噪声压制方法。

背景技术

有关于地震资料噪声压制，并已经成熟应用到实际地震资料处理中的方法主要包括：自适应滤波，f-k滤波，多项式拟合，独立成分分析，时频峰值滤波等。M.Naghizadeh在Seismicdatainterpolationanddenoisinginthefrequency-wavenumberdomain中提出的f-k滤波，根据信号与噪声的视速度差异完成噪声压制，主要应用于去除频率较低的面波噪声。多项式拟合的拟合中心点随正交多项式系数而改变，导致处理后数据主频降低、断点变形。

随着基于多尺度分析的小波变换方法被引入到地震信号去噪中来，多分辨率分析方法为地震资料处理开辟了新的方向。ShucongL和ErgenG等在SeismicDataDenoisingSimulationResearchBasedonWaveletTransform所使用的小波变换，在处理二维及更高维的信号方面存在方向局限性，并不能很好的描述有效信号的方向信息，并且二维滤波方法时常会引起信号畸变，导致出现虚假同相轴。

近年来如Ridgelet变换、Curvelet变换、Contourlet变换、Sharlet变换等具有多分辨率、多方向性、多尺度化的几何分析法成为热门方法，具体可以参考：王梅和侯振杰等的基于Ridgelet变换的多尺度去噪算法研究；孟阁阁，王德利等的基于2DCurvelet变换的多炮地震数据去噪方法研究；刘成明等的基于Shearlet变换的地震随机噪声压制。它们通过少数非零变换系数重构有效信号，在各向异性上实现了良好的稀疏表示。但上述方法采用固定基函数对信号进行分解，忽略了地震同相轴固有的时空连续性，导致噪声压制效果不佳。

发明内容

本发明的目的在于解决现有方法由于没有考虑到地震资料同相轴固有的局部和非局部相似性，造成的噪声压制不充分和有效信息衰减的问题，提出一种基于地震资料同相轴固有的局部和非局部相似性的，可以明显提高地震资料信噪比，并能较大程度保留有效信息的噪声压制方法。

本发明的目的是通过以下方案实现的：

一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，包括以下步骤：

步骤一、相似性分组：将整个地震资料切割成过完备的子块，遍历地选择参考子块，计算参考子块与一定邻域内的其他子块的欧氏距离，判断两子块是否相似，并把相似二维子块按照三维数据的方式存放起来成为相似组；

步骤二、三维奇异值滤波去噪：对相似组中每个子块先进行二维离散余弦变换，选用频域二维奇异值分解滤波，再将相似块向量化，进行频域一维奇异值分解，并进行奇异值保留去噪；二维离散余弦反变换后，将过完备子块加权聚合成原始尺寸地震资料，得到基础估计；

步骤三、计算维纳收缩系数：将基础估计重新分块分组，并且通过新的相似组计算变换域的维纳收缩系数；

步骤四、维纳滤波：按照新的相似分组，对原始信号进行三维离散余弦变换，采用维纳滤波实现噪声压制，三维离散余弦反变换并聚合重构得到最终噪声压制后的地震资料。

所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，其中，步骤一相似性分组包括以下具体步骤：

1)预处理：先将整个地震资料按照一定的步进长度过完备的分割成N₁×N₁大小的时-空域子块Y_x；然后对各个子块的二维频谱进行阈值r′处理，将大于阈值的频谱值保留下来，小于阈值的频谱值置0；

2)计算块间距离：将预处理后的子块频谱r′(T(Y_x))作为计算块间距离的元素，选定一个参考块Y_R，选定搜索窗大小为N_h×N_h，通过2-范数的计算参考块与其相应搜索窗内其他子块的块间距离：

$d(Y_R,Y_x)=\frac{\left|\right|r^{\′}\left(T\left(Y_R\right)\right)-r^{\′}\left(T\left(Y_x\right)\right)|{|}_2^2}{{\left(N_1\right)}^2}$

3)堆叠成三维数组：当两个子块的距离d(Y_R,Y_x)小于阈值τ时，判定两块相似，然后将所有与同一参考块相似的二维子块，连同参考块一起存放成一个三维数组，即相似组Z_R。

所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，其中，步骤二三维奇异值滤波去噪包括以下具体步骤：

1)二维离散余弦变换：进行二维离散余弦变换，将相似组中各个子块变换到频域；

2)奇异值处理：

(1)二维奇异值分解：对于一个相似组中的任一二维子块A，A是N₁×N₁的，它的特征矩阵U_2D和V_2D分别根据它的行-行协方差矩阵和列-列协方差矩阵求得，将A_i,i∈[0,N₁-1]投影到U_2D和V_2D上，得到投影系数Λ_i，满足：

$A_i=U_{2D}{\mathrm{\Λ}}_{2D}{V}_{2D}^T$

Λ_i即奇异值矩阵，将奇异值矩阵中大于阈值T₁的奇异值保留，小于T₁的奇异值置0，然后重构；

(2)一维奇异值分解：将相似组中所有子块向量化，进行一维奇异值分解，获得奇异值矩阵并进行重构；

3)离散余弦反变换、加权聚合：经所述步骤2)处理后的奇异值矩阵重构得到的是噪声压制后的频域数据块，将其进行二维离散余弦反变换，然后对二维离散余弦反变换结果进行加权重构，得到实现大部分噪声压制的该相似组所有子块对应位置的时域的的基础估计

${\hat{Y}}^{basic}=\frac{\underset{Y_R}{\Σ}\underset{Z_R}{\Σ}{\mathrm{\ω}}_R^{ht}\×{\hat{Y}}_R}{\underset{Y_R}{\Σ}\underset{Z_R}{\Σ}{\mathrm{\ω}}_R^{ht}}$

其中，N_R是未置0的系数个数

所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，其中，步骤三计算维纳收缩系数具体包括以下步骤：

1)重新分块：将经过所述步骤二获得的基础估计过完备地分成N₁×N₁的子块选择参考块重新计算子块间欧氏距离：

$d({\hat{Y}}_R^{basic},{\hat{Y}}_x^{basic})=\frac{\left|\right|{\hat{Y}}_R^{basic}-{\hat{Y}}_x^{basic}|{|}_2^2}{{\left(N_1\right)}^2}$

其中，是经过所述步骤二得到的基础估计；

2)堆叠成三维数组：当两个子块的距离小于阈值τ时，判定两块相似，并将分组方法传递给原始含噪地震资料，将判定为相似的子块作为一个数组存放起来形成一个新的相似组

3)三维离散余弦变换：首先对相似组中每个二维块独立地进行二维离散余弦变换，然后对所有子块的同一位置的数据点进行一维离散余弦变换，三维离散余弦变换后，地震资料被变换到频域；

4)计算维纳收缩系数：

根据基础估计计算维纳收缩系数W_R：

$W_R=\frac{\left|T_{3D}{\left({\hat{Y}}_R^{basic}\right)}^2\right|}{{\left|T_{3D}\left({\hat{Y}}_R^{basic}\right)\right|}^2+{\mathrm{\σ}}^2}$

其中，σ²是噪声方差，T_3D为三维离散余弦变换，是的分组块。

所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，其中，步骤四维纳滤波具体包括以下步骤：

1)维纳滤波：在所述步骤三中新的相似组已经经过三维离散余弦变换，在此将频域块与对应的维纳收缩系数相乘实现滤波；维纳滤波的计算方法是：

其中，σ²是噪声方差，T_3D为三维离散余弦变换，是经过所述步骤二获得的基础估计的分组块；

2)三维离散余弦反变换：按照先一维离散余弦反变换，再二维离散余弦反变换的倒序方式，对每个子块依次实现反变换；

3)加权聚合：对二维离散余弦反变换结果进行加权聚合，得到最终噪声压制后的地震资料

${\hat{Y}}^{wie}=\frac{\underset{Y_R}{\Σ}\underset{Z_R^{wie}}{\Σ}{\mathrm{\ω}}_R^{wie}\×{\hat{Y}}_R^{wie}}{\underset{Y_R}{\Σ}\underset{Z_R^{wie}}{\Σ}{\mathrm{\ω}}_R^{wie}}$

其中， ${\mathrm{\ω}}_R^{wie}=\frac{1}{{\mathrm{\σ}}^2\×\left|\right|W_R|{|}_2^2}$

与现有技术方案相比，本发明的方法在理论上采用块匹配算法，考虑到地震资料的固有特点——同相轴的局部和非局部相似性；采用二维奇异值分解和一维奇异值分解结合，充分利用地震资料的结构信息，实现地震资料的稀疏表示，能够较大程度地减小有效信号衰减，并提高噪声压制效果。从数据方面来看，通过仿真模型可以看出，表征信号保留程度的均方根误差减小了0.002，表征噪声压制效果的信噪比增加了大约3.4db，在250～500Hz的高频段噪声压制效果尤其明显。

附图说明

图1为本发明方法流程图。

图2为块匹配图示。

图3为纯净地震仿真信号

图4为图3所示仿真信号的加入-5db噪声地震仿真信号

图5为原始方法处理后的地震仿真信号

图6为本发明方法处理后的地震仿真信号

图7为第30道仿真信号处理对比图

图8为第30道仿真信号放大对比图

图9为第30道信号幅度谱对比图

图10为实际带噪地震资料

图11为原始方法处理后的实际地震资料

图12为本发明方法处理后的实际地震资料

图13为实际信号放大对比图

具体实施方式

本发明提出的地震资料噪声压制方法，是一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，该算法结合了三维块匹配算法和奇异值分解算法。

本发明中块匹配是指将待处理信号过完备的分成特定大小的子块，遍历性地在一定搜索邻域内寻找满足某种相似规范的相似块，然后将所有与同一参考块匹配的二维块与参考块一起，按照三维数组的形式存放起来。过完备是指分成的子块是相互重叠的，其中心点数据的序号差为分块时选择的步进长度，如图2所示。遍历性地搜索是指按顺序将每一个子块作为参考块，考虑后续的计算过程。相似性分组的处理方法充分的考虑到同相轴的非局部相似性。

在频率域，不同道的地震信号具有相同的有效频带，并且在地震道的方向上这些频率具有可预测性、横相干性；而噪声是随机的，频带较宽，也不存在可预测性和相关性。信号在频率域能得到比时域更稀疏的表示方法，由此，本发明采用在频率域进行奇异值分解去噪。

考虑到从属于不同块的同相轴间具有非局部相似性，本发明采用一维奇异值分解实现同相轴非局部冗余的去除。一维奇异值分解充分保留了地震资料的非局部结构信息，在实现噪声压制的同时，也能较大程度地保留有效信息。

考虑到从属于同一块的同相轴间具有局部相似性，本发明采用二维奇异值分解实现同相轴局部冗余的去除。二维奇异值分解充分保留了地震资料的局部结构信息，在实现噪声压制的同时，也较大程度地保留了有效信息。

三维块匹配算法和奇异值分解相结合充分考虑了同相轴局部和非局部的相似性，避免了现有方法中没有考虑同相轴局部和非局部的相似性而导致的噪声压制不充分和信号畸变。

本发明提出的基于块匹配算法和奇异值分解的地震资料噪声压制方法主要分为4步，包括相似性分组、三维变换奇异值滤波去噪、根据基础估计重新分组并计算维纳收缩系数和维纳滤波四个步骤。如图1所示。

步骤一、相似性分组：将整个地震资料切割成过完备的子块，遍历地选择参考子块，计算参考子块与一定邻域内的其他子块的欧氏距离，判断两子块是否相似，并把相似二维子块按照三维数据的方式存放起来成为三维数组。

步骤二、三维奇异值滤波去噪：对相似组中每个子块先进行二维离散余弦变换，选用频域二维奇异值分解滤波。再将相似块向量化，进行频域一维奇异值分解，进行奇异值保留去噪。二维离散余弦反变换后，将过完备子块加权聚合成原始尺寸地震资料，得到基础估计。

步骤三、计算维纳收缩系数：将基础估计重新分块分组，并且通过新的相似分组计算变换域的维纳收缩系数。

步骤四、维纳滤波：按照新的相似分组，对原始信号进行三维离散余弦变换，采用维纳滤波实现噪声压制，三维离散余弦反变换并聚合重构得到最终估计，即最终噪声压制后的地震资料。

下面对各步骤进一步介绍如下：

步骤一、相似性分组：

1.预处理

先将整个地震资料按照一定的步进长度过完备的分割成N₁×N₁大小的时-空域子块Y_x，准备计算其相似性差值。但是，当地震资料中包含噪声较多时，可能由于噪声方差过大而引起相似规范的计算误差，从而引起分组错误。因此，在计算块间距离之前，首先对各个子块的二维频谱进行阈值r′处理。将大于阈值的频谱值保留下来，小于阈值的频谱值置0。预处理过程能够减小下一步计算块间距离过程中噪声对相似判断的错误影响。

2.计算块间距离

相似性的规范一般采用计算信号两个块之间的一些距离参数，例如欧几里得距离。本发明将预处理后的子块频谱r′(T(Y_x))作为计算块间距离的元素，选定一个参考块Y_R，选定搜索窗大小为N_h×N_h，通过2-范数的计算参考块与其相应搜索窗内其他子块的块间距离：

$d(Y_R,Y_x)=\frac{\left|\right|r^{\′}\left(T\left(Y_R\right)\right)-r^{\′}\left(T\left(Y_x\right)\right)|{|}_2^2}{{\left(N_1\right)}^2}$

通过计算块间距离能够得知参考块邻域中的其他子块与这一参考块的相似程度，为下一步相似性分组做好了准备。

3.堆叠成三维数组

通过计算块间距离，已经得知参考块的一定邻域内其他块与当前参考块的相似程度，为了后续滤波过程的顺利实施，需要将判定为相似的时-空域子块作为一个数组存放起来。当两个子块的距离d(Y_R,Y_x)小于阈值τ时，判定两块相似，然后将所有与同一参考块相似的二维子块，连同参考块一起存放成一个三维数组，即相似组Z_R。经过分组后，具有非局部相似性的子块被存放在一起，为一维奇异值分解做好了准备。

步骤二、三维奇异值滤波去噪

1.二维离散余弦变换

对相似组Z_R进行三维奇异值滤波之前，首先进行二维离散余弦变换，将相似组Z_R中各个子块变换到频域。上面提到过在频率域，不同道的地震信号具有相同的有效频带，并且在地震道的方向上这些频率具有可预测性、横相干性；而噪声是随机的，频带较宽，也不存在可预测性和相关性。也就是说在频域中，有效信号的能量能更集中在某几个少数不为0的点上，即稀疏表示。这样做扩大了有效信息与噪声的差异，能够增强奇异值滤波的滤波效果。

2.奇异值处理

(1)二维奇异值分解

经过二维离散余弦变换，子块的地震数据变换到频域，实现初步的稀疏表示。对于一个相似组中的任一二维子块A，A是N₁×N₁的，它的特征矩阵U_2D和V_2D分别根据它的行-行协方差矩阵和列-列协方差矩阵求得，将A_i,i∈[0,N₁-1]投影到U_2D和V_2D上，得到投影系数Λ_i，满足

$A_i=U_{2D}{\mathrm{\Λ}}_i{V}_{2D}^T$

Λ_i即奇异值矩阵。奇异值分解使得地震资料能用更少的非零系数表示完全，实现了进一步稀疏表示。而且，奇异值矩阵中较大的奇异值代表地震资料中的有效信息，较小的奇异值代表噪声信息。因此，将大于阈值T₁的奇异值保留，小于T₁的奇异值置0，然后重构，这一步考虑了局部相似性，实现了块内冗余的去除。

(2)一维奇异值分解

二维奇异值分解只考虑到地震资料的局部相似性，而没有考虑到其非局部相似性。因此，将相似组中所有子块向量化，即将二维数据子块变成一个列向量，进行一维奇异值分解。

每个子块是N₁×N₁的，假设此参考块有M-1个相似块，那么，构造出来的相似矩阵B的尺寸为N₁×N₁行，M列。对于此矩阵必然存在两个正交矩阵 $U_{1D}\∈R^{N_1\×N_1}$ 和 $V_{1D}\∈R^{N_1\×N_1},$ 满足

$B=U_{1D}{\mathrm{\Λ}}_{1D}{V}_{1D}^T$

Λ_1D即奇异值矩阵，同样将Λ_1D中将大于阈值T₂的奇异值保留，小于T₂的奇异值置0，然后重构，这一步骤考虑了同相轴的非局部相似性，实现了块间冗余的去除。

3.离散余弦反变换、加权聚合

处理后的奇异值矩阵重构得到的是噪声压制后的频域数据块，将其进行二维离散余弦反变换，得到该相似组所有子块对应位置的时域的基础估计但块是相互重叠的，因此需要加权重构。对应权重计算方法是：

${\mathrm{\ω}}_R^{ht}=\frac{1}{{\mathrm{\σ}}^2\×N_R}$

由于噪声方差σ²和未置0的系数个数N_R的存在，权重代表了该参考块中有效信息的占比，有效信息越多，权重越大。相互重叠的过完备块的基础估计按照下式方法聚合。

${\hat{Y}}^{basic}=\frac{\underset{Y_R}{\Σ}\underset{Z_R}{\Σ}{\mathrm{\ω}}_R^{ht}\×{\hat{Y}}_R}{\underset{Y_R}{\Σ}\underset{Z_R}{\Σ}{\mathrm{\ω}}_R^{ht}}$

由此，得到已经实现大部分噪声压制的基础估计

步骤三、计算维纳收缩系数

1.重新分块

由于基础估计已经实现大部分的噪声压制，能够实现更精确的相似性分组。将基础估计过完备地分成N₁×N₁的子块选择参考块按照同样的方法重新计算子块间欧氏距离：

$d({\hat{Y}}_R^{basic},{\hat{Y}}_x^{basic})=\frac{\left|\right|{\hat{Y}}_R^{basic}-{\hat{Y}}_x^{basic}|{|}_2^2}{{\left(N_1\right)}^2}$

这里是上一步得到的基础估计，得到的是基础估计中每一参考块与其周围邻域中各个子块的相似性差值，为下一步相似性分组做好准备。

2.堆叠成三维数组

通过计算块间距离，已经得知参考块的一定邻域内其他块与当前参考块的相似程度，为了后续滤波过程的顺利实施，需要将判定为相似的子块作为一个数组存放起来。当两个子块的距离小于阈值τ时，判定两块相似，并将分组方法传递给原始含噪地震资料，形成一个新的三维数组重新分组后得到更为精确的非局部相似性，能更为准确地保留地震资料的非局部结构信息。

3.三维离散余弦变换

在进行滤波处理之前，首先将地震资料由时域变换到频域，本发明采用的是三维离散余弦变换。三维离散余弦变换是首先对相似组中每个二维块独立地进行二维离散余弦变换，然后对所有块的同一位置的数据点进行一维离散余弦变换。三维离散余弦变换后，地震资料被变换到频域，实现了初步的稀疏表示。

4.计算维纳收缩系数

另外，根据基础估计计算下一步滤波所需要的维纳收缩系数W_R。

$W_R=\frac{\left|T_{3D}{\left({\hat{Y}}_R^{basic}\right)}^2\right|}{{\left|T_{3D}\left({\hat{Y}}_R^{basic}\right)\right|}^2+{\mathrm{\σ}}^2}$

这里σ²是噪声方差，T_3D为三维离散余弦变换，是基础估计的分组块。通过维纳系数的计算方式可以看出，维纳系数代表了该子块的频域能量与估计的含噪信息的频域能量的比值，即有效信息越多，维纳收缩系数越大，这也是它的去噪原理。

步骤四、维纳滤波

1.维纳滤波

在第三个大步骤中，通过各个小步的计算，已经完成了所有维纳滤波的准备过程，在这一部分中将实现维纳滤波去除地震资料中的噪声。维纳滤波的计算方法是：

在上述步骤三中分组块已经经过三维离散余弦变换，在此需要将频域块与对应的维纳收缩系数相乘实现滤波。上面已经提到包含有效信息的数据点对应的维纳收缩系数较大，也就是说通过维纳滤波子块中代表有效信息的数据点得以保留，而代表噪声的数据点被压制，实现了地震资料子块的噪声去除。

2.三维离散余弦反变换

进行三维离散余弦反变换将频域子块变回到时域，得到过完备的块的估计。三维离散余弦反变换过程按照先一维离散余弦反变换，再二维离散余弦反变换的倒序方式，对每个子块依次实现反变换。但此时块的估计是相互重叠的，需要下一步的加权聚合得到最终的噪声压制结果。

3.加权聚合

同样，对过完备块进行加权求和，权重是根据维纳收缩系数计算得来的。

${\mathrm{\ω}}_R^{wie}=\frac{1}{{\mathrm{\σ}}^2\×\left|\right|W_R|{|}_2^2}$

权重同样代表该块中有效信息的占比。聚合方法与之前类似：

${\hat{Y}}^{wie}=\frac{\underset{Y_R}{\Σ}\underset{Z_R^{wie}}{\Σ}{\mathrm{\ω}}_R^{wie}\×{\hat{Y}}_R^{wie}}{\underset{Y_R}{\Σ}\underset{Z_R^{wie}}{\Σ}{\mathrm{\ω}}_R^{wie}}$

这样就得到了最终噪声压制后的地震资料

实施例

一.仿真信号模型实施例

仿真条件如下：在Matlab实验平台，采用100道，每道2000个采样点的合成地震信号仿真模型来进行算法的性能验证。该记录包含三条同相轴，主频分别为45Hz、35Hz和25Hz，视速度分别为1000m/s、1400m/s和1500m/s，采样频率是1000Hz，炮检距为10m。不包含噪声的纯净地震信号仿真模型如图3所示。

为验证本发明方法对低信噪比地震资料的噪声压制效果，在纯净地震信号仿真模型中加入-5db的高斯白噪声，如图4所示，可以看到有效信息大量淹没在噪声之中。

1、相似性分组

将2000×100的带噪地震仿真模型按照块步进长度为2，分成1000×50个8×8的子块。首先将每个子块通过二维傅里叶变换转换到频域，经过阈值为0.9的预处理。频域值大于0.9的数据点保留，小于0.9的数据点置零。然后以每个子块为中心，在其周围19×19的邻域内搜索其相似块，相似规范为2-范数距离。把2-范数距离小于3的判为两块相似，与其对应的参考块存放在一起作为一个相似组，而距离大于3的则不做处理。

2、三维奇异值滤波去噪

将时域同一相似组中的每个子块进行二维离散余弦变换转换到频域，实现初步的稀疏表示。然后对每个子块独立的进行二维奇异值分解，以1.2作为阈值处理二维奇异值，大于1.2的奇异值保留下来，小于1.2的奇异值作为噪声分量置0，然后重构回去，得到该相似组中每个子块的估计，是去除块内冗余的结果。然后将所有子块向量化，即将二维数据块转化为一个列向量，然后依次按列存放成一个二维数据块，做一维奇异值分解。这一步以1作为阈值，重构得到的是去除这一相似组中块间冗余的结果。由于子块是相互重叠的，按照上文提到的公式计算权重，加权求和即得到整个地震仿真模型的基础估计。

3、计算维纳收缩系数

将2000×100的地震资料基础估计按照块步进长度为2，分成1000×50个8×8的子块。然后以每个子块为中心，在其周围19×19的邻域内搜索其相似块，相似规范为2-范数距离。把2-范数距离小于1的判为两块相似，与其对应的参考块存放在一起作为一个相似组，而距离大于1的则不做处理。同时将这种分组方法，通过子块序号的方式传递给原始带噪地震资料，将原始带噪地震资料重新分组。

将原始带噪地震资料新的相似组中的每个二维块进行二维离散余弦变换，然后再对每个二维块相同位置的数据点进行一维离散余弦变换，将时域信号变换到频域，实现地震资料的初步稀疏表示。并通过各数据块的频谱计算对应的维纳收缩系数，为下一步滤波做好准备工作。

4、维纳滤波

将每一子块与其对应的维纳收缩系数相乘实现滤波。这一步不需要阈值处理，因为维纳收缩系数的计算方法使得包含有效信息的点系数较大，噪声点系数较小。因此子块与其相乘就能实现有效信息的保留和噪声压制效果。滤波后，将频域子块按照先一维离散余弦反变换，再二维离散余弦反变换的顺序变换回时域。然后将重叠的子块估计加权聚合，就能得到2000×100的地震资料仿真模型最终的噪声压制效果。

经过原始三维块匹配算法处理后的噪声压制效果如图5所示，可以看到与图3的纯净信号相比，图5仍包含少量噪声。而图6，本发明算法处理后的地震信号明显剩余噪声更少，有更强的噪声压制效果。表1从均方根误差和信噪比两个方面，量化地证明本发明算法取得的更好的效果。

表1去噪效果对比表

	均方根误差	信噪比
			原始三维块匹配	0.025	13.43db
本发明算法	0.023	16.86db

图7在100道地震信号仿真模型中选取了第30道地震信号，图8放大其中一条同相轴进行观察。可以清晰地看到本发明算法处理后的信号曲线在不含信号部分更为平滑，包含信号的部分曲线与纯净信号曲线更为贴合。这就证明本发明算法有更好的噪声压制效果，同时也能更好地保留有效信号的结构信息。

图9从第30道地震信号的幅度谱的角度证明这一结论。从幅度谱来看，本发明算法在250Hz～500Hz的频率段更为贴合纯净信号的幅度谱，而原始三维块匹配在这一频段仍包含噪声。

二.实际地震资料实施例

实际地震资料是一张168道，每道6000点的共炮点记录，如图10。采样频率为1000Hz，炮检距为30m。该资料中包含随机噪声、较为混乱的同相轴和设备带来的其他噪声。

1、相似性分组

将6000×168的带噪地震仿真模型按照块步进长度为4，分成1500×42个8×8的子块。首先将每个子块通过二维傅里叶变换转换到频域，经过阈值为1×10^-5的预处理。频域值大于1×10^-5的数据点保留，小于1×10^-5的数据点置零。然后以每个子块为中心，在其周围39×39的邻域内搜索其相似块，相似规范为2-范数距离。把2-范数距离小于2×10^-9的判为两块相似，与其对应的参考块存放在一起作为一个相似组，而距离大于2×10^-9的则不做处理。

2、三维奇异值滤波去噪

将时域同一相似组中的每个子块进行二维离散余弦变换转换到频域，实现初步的稀疏表示。然后对每个子块独立的进行二维奇异值分解，以1×10^-4作为阈值处理二维奇异值，大于1×10^-4的奇异值保留下来，小于1×10^-4的奇异值作为噪声分量置0，然后重构回去，得到该相似组中每个子块的估计，是去除块内冗余的结果。然后将所有子块向量化，即将二维数据块转化为一个列向量，然后依次按列存放成一个二维数据块，做一维奇异值分解。这一步以0.0012作为阈值，重构得到的是去除这一相似组中块间冗余的结果。由于子块是相互重叠的，按照上文提到的公式计算权重，加权求和即得到整个地震仿真模型的基础估计。

3、计算维纳收缩系数

将6000×168的地震资料基础估计按照块步进长度为4，分成1500×42个8×8的子块。然后以每个子块为中心，在其周围39×39的邻域内搜索其相似块，相似规范为2-范数距离。把2-范数距离小于1×10^-15的判为两块相似，与其对应的参考块存放在一起作为一个相似组，而距离大于1×10^-15的则不做处理。同时将这种分组方法，通过子块序号的方式传递给原始带噪地震资料，将原始带噪地震资料重新分组。

将原始带噪地震资料新的相似组中的每个二维块进行二维离散余弦变换，然后再对每个二维块相同位置的数据点进行一维离散余弦变换，将时域信号变换到频域，实现地震资料的初步稀疏表示。并通过各数据块的频谱计算对应的维纳收缩系数，为下一步滤波做好准备工作。

4、维纳滤波

将每一子块与其对应的维纳收缩系数相乘实现滤波。这一步不需要阈值处理，因为维纳收缩系数的计算方法使得包含有效信息的点系数较大，噪声点系数较小。因此子块与其相乘就能实现有效信息的保留和噪声压制效果。滤波后，将频域子块按照先一维离散余弦反变换，再二维离散余弦反变换的顺序变换回时域。然后将重叠的子块估计加权聚合，就能得到6000×168的地震资料仿真模型最终的噪声压制效果。

原始三维块匹配和本发明算法处理后的效果分别如图11和图12所示。可以清晰看到本发明算法处理后的地震资料有更为清晰的同相轴。放大后观察，也能看出图13(b)本算法处理后含的噪声比图13(a)原始三维块匹配算法处理后的噪声更少。

Claims (5)

1.一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，其特征在于，包括以下步骤：

步骤一、相似性分组：将整个地震资料切割成过完备的子块，遍历地选择参考子块，计算参考子块与一定邻域内的其他子块的欧氏距离，判断两子块是否相似，并把相似二维子块按照三维数据的方式存放起来成为相似组；

步骤二、三维奇异值滤波去噪：对相似组中每个子块先进行二维离散余弦变换，选用频域二维奇异值分解滤波，再将相似块向量化，进行频域一维奇异值分解，并进行奇异值保留去噪；二维离散余弦反变换后，将过完备子块加权聚合成原始尺寸地震资料，得到基础估计；

步骤三、计算维纳收缩系数：将基础估计重新分块分组，并且通过新的相似组计算变换域的维纳收缩系数；

步骤四、维纳滤波：按照新的相似分组，对原始信号进行三维离散余弦变换，采用维纳滤波实现噪声压制，三维离散余弦反变换并聚合重构得到最终噪声压制后的地震资料。

2.如权利要求1所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，其特征在于，所述步骤一相似性分组包括以下具体步骤：

1)预处理：先将整个地震资料按照一定的步进长度过完备的分割成N₁×N₁大小的时-空域子块Y_x；然后对各个子块的二维频谱进行阈值r′处理，将大于阈值的频谱值保留下来，小于阈值的频谱值置0；

2)计算块间距离：将预处理后的子块频谱r′(T(Y_x))作为计算块间距离的元素，选定一个参考块Y_R，选定搜索窗大小为N_h×N_h，通过2-范数的计算参考块与其相应搜索窗内其他子块的块间距离：

$d(Y_R,Y_x)=\frac{\left|\right|r^{\′}\left(T\right(Y_R\left)\right)-r^{\′}\left(T\right(Y_x\left)\right)|{|}_2^2}{{\left(N_1\right)}^2}$

3)堆叠成三维数组：当两个子块的距离d(Y_R,Y_x)小于阈值τ时，判定两块相似，然后将所有与同一参考块相似的二维子块，连同参考块一起存放成一个三维数组，即相似组Z_R。

3.如权利要求1所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，其特征在于，所述步骤二三维奇异值滤波去噪包括以下具体步骤：

1)二维离散余弦变换：进行二维离散余弦变换，将相似组中各个子块变换到频域；

2)奇异值处理：

(1)二维奇异值分解：对于一个相似组中的任一二维子块A，A是N₁×N₁的，它的特征矩阵U_2D和V_2D分别根据它的行-行协方差矩阵和列-列协方差矩阵求得，将A_i,i∈[0,N₁-1]投影到U_2D和V_2D上，得到投影系数Λ_i，满足：

$A_i=U_{2D}{\mathrm{\Λ}}_i{V}_{2D}^T$

Λ_i即奇异值矩阵，将奇异值矩阵中大于阈值T₁的奇异值保留，小于T₁的奇异值置0，然后重构；

(2)一维奇异值分解：将相似组中所有子块向量化，进行一维奇异值分解，获得奇异值矩阵并进行重构；

3)离散余弦反变换、加权聚合：经所述步骤2)处理后的奇异值矩阵重构得到的是噪声压制后的频域数据块，将其进行二维离散余弦反变换，然后对二维离散余弦反变换结果进行加权重构，得到实现大部分噪声压制的该相似组所有子块对应位置的时域的的基础估计

${\hat{Y}}^{basic}=\frac{\underset{Y_R}{\Σ}\underset{Z_R}{\Σ}{\mathrm{\ω}}_R^{ht}\×{\hat{Y}}_R}{\underset{Y_R}{\Σ}\underset{Z_R}{\Σ}{\mathrm{\ω}}_R^{ht}}$

其中，N_R是未置0的系数个数。

4.如权利要求1所述的一种基于三维块匹配算法和奇异值分解的地震资料噪声压制方法，其特征在于，所述步骤三计算维纳收缩系数具体包括以下步骤：

1)重新分块：将经过所述步骤二获得的基础估计过完备地分成N₁×N₁的子块选择参考块重新计算子块间欧氏距离：

$d({\hat{Y}}_R^{basic},{\hat{Y}}_x^{basic})=\frac{\left|\right|{\hat{Y}}_R^{basic}-{\hat{Y}}_x^{basic}|{|}_2^2}{{\left(N_1\right)}^2}$

其中，是经过所述步骤二得到的基础估计；

2)堆叠成三维数组：当两个子块的距离小于阈值τ时，判定两块相似，并将分组方法传递给原始含噪地震资料，将判定为相似的子块作为一个数组存放起来形成一个新的相似组

3)三维离散余弦变换：首先对相似组中每个二维块独立地进行二维离散余弦变换，然后对所有子块的同一位置的数据点进行一维离散余弦变换，三维离散余弦变换后，地震资料被变换到频域；

4)计算维纳收缩系数：

根据基础估计计算维纳收缩系数W_R：

$W_R=\frac{\left|T_{3D}{\left({\hat{Y}}_R^{basic}\right)}^2\right|}{|T_{3D}\left({\hat{Y}}_R^{basic}\right){|}^2+{\mathrm{\σ}}^2}$

其中，σ²是噪声方差，T_3D为三维离散余弦变换，是的分组块。

5.如权利要求1所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法，其特征在于，所述步骤四维纳滤波具体包括以下步骤：

1)维纳滤波：在所述步骤三中新的相似组已经经过三维离散余弦变换，在此将频域块与对应的维纳收缩系数相乘实现滤波；维纳滤波的计算方法是：

其中，σ²是噪声方差，T_3D为三维离散余弦变换，是经过所述步骤二获得的基础估计的分组块；

2)三维离散余弦反变换：按照先一维离散余弦反变换，再二维离散余弦反变换的倒序方式，对每个子块依次实现反变换；

3)加权聚合：对二维离散余弦反变换结果进行加权聚合，得到最终噪声压制后的地震资料

${\hat{Y}}^{wie}=\frac{\underset{Y_R}{\Σ}\underset{Z_R^{wie}}{\Σ}{\mathrm{\ω}}_R^{wie}\×{\hat{Y}}_R^{wie}}{\underset{Y_R}{\Σ}\underset{Z_R^{wie}}{\Σ}{\mathrm{\ω}}_R^{wie}}$

其中， ${\mathrm{\ω}}_R^{wie}=\frac{1}{{\mathrm{\σ}}^2\×\left|\right|W_R|{|}_2^2}$