CN110068865B - 一种低秩矩阵逼近的沙漠地震噪声压制方法 - Google Patents

Abstract

本发明公开了一种基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，包括以下步骤：步骤一、获取含噪地震信号Y，设置参数δ，η，L，τ；步骤二、设置Y⁰＝Y,

步骤三、对k＝1:L做迭代正则化：

步骤四、选择一个主要地震纹理块

步骤五、构建相似块矩阵

步骤六、估计噪声标准差；步骤七、奇异值分解；步骤八、估计奇异值

步骤九、获得奇异值矩阵

步骤十、用阈值τ将信号的较大奇异值和较小奇异值分开；步骤十一、获得估计值

和

步骤十二、利用

和

来获得

和

得出去噪的地震信号

本发明提供的压制方法能够根据含噪地震信号不同纹理块的几何结构直接估计出噪声标准差，对沙漠地震低频噪声压制有明显的效果。

Description

一种低秩矩阵逼近的沙漠地震噪声压制方法

技术领域

本发明涉及地球物理技术领域，尤其涉及一种基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法。

背景技术

在地震信号的收集过程中，不可避免的会受到随机噪声的影响，这会在很大的程度上影响到后续地震资料的分析，如果随机噪声得到了抑制，那么后续对地震资料的分析将会更加准确，因此随机噪声的压制是地震数据应用前的一个非常重要的步骤。沙漠地区随机噪声具有低频，非平稳等特点，而且随机噪声的频带和有效地震信号的频带存在严重混叠，这使得一些适用于抑制白噪声的方法和在频域分出信号和噪声的方法在沙漠地区地震信号噪声压制上不能取得很好的效果。

目前已有大量的方法被提出和应用到实际地震资料中随机噪声的抑制。例如其一，基于稀疏变换的方法如小波变换、curvelet变换，是将地震数据变换到稀疏域中获得相应的稀疏系数，然后对稀疏系数进行阈值处理，再将处理后的稀疏系数变换回时空域中以达到抑制噪声的目的。其二，基于模态分解的方法，如经验模态分解(EMD)、变分模态分解(VMD)，主要是将地震数据分解为许多不同的成分分量，然后选出含有有效信号的分量而舍弃随机噪声的分量以达到抑制噪声的目的。其三，基于预测的方法，如f-x反褶积，主要是利用有效信号的可预测性构建预测滤波器来抑制噪声和增强信号。上述已有方法在应对大套油气层，中浅层勘探随机噪声消减，提高地震记录信噪比方面起到了一定作用，在地震勘探工程中得到较好的效果。但是这些方法由于受某些条件上的约束和自身方法的限制，在沙漠地区地震噪声压制中的应用不能取得很好的实际效果。

近些年来，有大量学者对于信号低秩的方法展开了深入的研究，可以认为地震信号在经过某种信号重组的方法之后是低秩的，如通过堆叠来自地震信号的非局部相似块构建成的矩阵在所给高维空间的低维子空间中是低秩的。如通过EMD经验地将一个不低秩的多倾地震图像分解为多个低秩的单倾地震图像。如利用地震信号构建的hankel矩阵会有一个更低的秩。这些方法去噪的机理也是不同的，考虑到信号奇异值的物理意义，根据不同大小的奇异值分配不同的权重值。通过EMD分解后对每个模态成分构建hankel矩阵，对这个hankel矩阵进行SVD分解后保留与模态成分秩数相同的较大的奇异值，将其它的奇异值置0。将通过地震信号构建的hankel矩阵分解为两个低维子空间中的矩阵的乘积，然后通过加权代价函数来求解这两个矩阵。

噪声估计是低秩矩阵逼近去噪框架中十分重要的一个环节，估计出来的噪声标准差会用于计算纯净信号的奇异值和相应的权重阈值，这些权重阈值是去噪的决定性参数，所以噪声估计的准确性会直接关系到去噪结果的好坏。现有的噪声估计方法是通过计算含噪地震记录堆叠相似块矩阵与当前迭代去噪的对应矩阵的欧式距离，直观的将它认为是噪声，这种方法的去噪结果受噪声水平初值设置的影响比较大，而且当前迭代去噪的结果会影响到下次迭代去噪的结果，如果噪声水平初值设置的不准确，或者是在某次迭代中滤除的不仅仅只有噪声，那么去噪的结果将会不是很好。同时，考虑到奇异值的物理意义，有效信息的奇异值大小和信号的能量有关，一般能量较高的信号的奇异值较大，能量较低的信号的奇异值较小，现有的噪声估计方法中，较小的奇异值在经过权重阈值处理后会变成0，这就会导致低能量信号的损失。尤其是在实际地震资料处理中，由于实际地震资料的情况较为复杂，如果还是用原来的噪声估计方式会不太准确。

发明内容

本发明为解决目前的技术不足之处，提供了一种基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，能够通过根据含噪地震信号不同纹理块的几何结构直接估计出噪声标准差，无需手动设置噪声水平初值，也不用考虑前次迭代去噪的效果，低秩矩阵逼近不需要考虑信号的频率，其在沙漠地震噪声压制上具有显著效果。

本发明提供的技术方案为：一种基于几何文理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，包括以下步骤：

步骤一、获取含噪地震信号Y，并且对所述含噪地震信号Y做k＝1:L的正则化处理：

其中，Y^k表示第k次迭代处理后的含噪地震信号，

表示第k次迭代处理后的去噪地震信号，

表示第k次迭代处理后的较大奇异值的去噪地震信号，

表示第k次迭代处理后的较小奇异值的去噪地震信号，δ，η，τ表示参数，并且设置初始值

步骤二、在含噪地震信号中设置搜索范围，选择主要地震纹理块矩阵，再选择与所述主要地震纹理块矩阵相似的地震纹理块矩阵构成含噪地震信号矩阵；

步骤三、所述地震纹理块矩阵的噪声标准差满足：

其中，

表示估计的噪声标准差，ω_min(·)表示矩阵的最小特征值，W表示挑选出来弱纹理块，cov(·)表示协方差矩阵；

步骤四、将所述含噪地震信号矩阵中的所有地震纹理块矩阵获得的噪声标准差初值，并且设置初始阈值，当所述地震纹理块的梯度协方差矩阵的最大特征值小于初始阈值时，所述地震纹理块为弱纹理块；

其中，所述地震纹理块的梯度协方差矩阵满足：

其中，

表示地震纹理块的梯度协方差矩阵，p_j表示含噪地震信号矩阵

中的纹理块矩阵，D_h和D_v分别表示水平和垂直微分算子；

初始阈值为人为设置，在第k次迭代处理后，阈值满足：

ρ＝F^-1(v,α,β)；

其中，F^-1表示逆伽马累积分布函数，v表示显著性水平，α表示伽马分布的形状参数，β表示伽马分布的尺度参数；

步骤五、根据所述弱纹理块的噪声标准差得到纯净地震信号的奇异值：

其中，

表示纯净地震信号矩阵的奇异值，

表示含噪地震信号矩阵的奇异值，

表示所述弱纹理块的噪声标准差，m是构成所述含噪地震信号矩阵的块的数目；

步骤六、所述奇异值分解得到奇异值矩阵，根据阈值将较大奇异值和较小奇异值分开并保留；

步骤七、得出去噪的地震信号。

优选的是，所述含噪地震信号矩阵中包括含噪平坦地震块，其分解为完美平坦地震块和噪声块；

其中，所述完美平坦地震块的梯度是0，所述含噪平坦地震块的梯度协方差矩阵的期望值为

其中，两个对角线元素

的统计特性是相同的，其中Ω＝h或v。

优选的是，所述伽马分布的形状参数α和尺度参数β分别为：

其中，α表示伽马分布的形状参数，β表示伽马分布的尺度参数，N表示噪声块中元素的个数，

表示噪声标准差，

代表矩阵

的迹。

优选的是，在所述步骤六中，所述奇异值的分解过程为：

其中，

表示含噪地震信号矩阵，U和V分别表示左奇异值向量矩阵和右奇异值向量矩阵，

表示第k次迭代处理后的纯净信号矩阵，w_i表示对应与奇异值的非负权重值，

Σ表示奇异值对角矩阵，

表示奇异值矩阵，

表示奇异值矩阵的元素，Σ_ii表示奇异值对角矩阵中的元素，并且所述奇异值矩阵中的元素满足：

其中，Σ_ii表示矩阵的元素，

表示含噪地震信号矩阵的奇异值。

优选的是，所述奇异值分为较大奇异值矩阵和较小奇异值矩阵，如下所示：

其中，

表示仅保留

中较大奇异值的矩阵，

表示仅保留

中较小奇异值的矩阵，并且设置阈值τ，奇异值小于所述阈值为较小奇异值的矩阵，奇异值大于所述阈值为较大奇异值的矩阵。

优选的是，通过所述较大奇异值矩阵和所述较小奇异值矩阵得到所述纯净信号矩阵：

其中，

表示第k次迭代处理后的纯净信号矩阵，

表示保留纯净信号矩阵中较高奇异值获得的矩阵，

表示保留纯净信号矩阵中较低奇异值获得的矩阵。

优选的是，在所述步骤七中，通过合并所述保留纯净信号矩阵中较低奇异值获得的矩阵和所述保留纯净信号矩阵中较高奇异值获得的矩阵得到保留较大奇异值和保留较小奇异值的两幅地震信号，将所述两幅地震信号合并得到第k次迭代处理后的纯净地震信号。

优选的是，在所述步骤七中，将所述两幅地震信号和所述第k次迭代处理后的纯净地震信号经过所述步骤一中的正则化处理，得到去噪的地震信号。

本发明与现有技术相比较所具有的有益效果：

本发明基于几何纹理的噪声估计方法，根据含噪地震信号不同纹理块的几何结构直接估计出噪声标准差，无需手动设置噪声水平初值，也不用考虑前次迭代去噪的效果，而且将奇异值矩阵分解为较大奇异值和较小奇异值，防止有效信号的丢失。

附图说明

图1为本发明所述合成地震记录的去噪理想信号示意图。

图2为本发明所述合成地震记录的含噪信号示意图。

图3为本发明所述合成地震记录的小波变换去噪结果示意图。

图4为本发明所述合成地震记录的f-x反褶积去噪结果示意图。

图5为本发明所述合成地震记录的WNNM去噪结果示意图。

图6为本发明所述合成地震记录的本发明方法去噪结果示意图。

图7为本发明所述合成地震记录的去噪理想结果的第三道单道示意图。

图8为本发明所述合成地震记录的含噪信号的第三道单道示意图。

图9为本发明所述合成地震记录的小波变换去噪结果的第三道单道示意图。

图10为本发明所述合成地震记录的f-x反褶积去噪结果的第三道单道示意图。

图11为本发明所述合成地震记录的WNNM去噪结果的第三道单道示意图。

图12为本发明所述合成地震记录的本发明方法去噪结果的第三道单道示意图。

图13为本发明所述合成记录的理想信号去噪结果的FK谱示意图。

图14为本发明所述合成记录的含噪信号的FK谱示意图。

图15为本发明所述合成记录的小波变换去噪结果的FK谱示意图。

图16为本发明所述合成记录的f-x反褶积去噪结果的FK谱示意图。

图17为本发明所述合成记录的WNNM去噪结果的FK谱示意图。

图18为本发明所述合成记录的本发明方法去噪结果的FK谱示意图。

图19为本发明所述实际沙漠地震的实际共炮点记录示意图。

图20为本发明所述实际沙漠地震的小波变换处理结果示意图。

图21为本发明所述实际沙漠地震的f-x反褶积处理结果示意图。

图22为本发明所述实际沙漠地震的WNNM处理结果示意图。

图23为本发明所述实际沙漠地震的本发明方法处理结果示意图。

图24为本发明所述的原始地震数据局部放大示意图。

图25为本发明所述的小波变换滤波结果局部放大示意图。

图26为本发明所述的f-x反褶积滤波结果局部放大示意图。

图27为本发明所述的WNNM滤波结果局部放大示意图。

图28为本发明所述的本发明方法滤波结果局部放大示意图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

本发明提供了一种基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，包括以下：

沙漠地区含噪地震数据可以用下面模型表示：

Y＝X+E， (1)

其中Y是含噪地震数据，X是潜在的纯净信号，E是方差为σ²的随机噪声。本发明的目的就是要从Y中恢复出X。由于在地震信号较为复杂的情况下X本身并不是低秩的，利用地震信号的非局部自相似性，在一特定大小的搜索窗内，选择一个主块

然后再选择与主块相似的块将它们堆叠成一个矩阵

同时可以找到对应的纯净信号矩阵

和随机噪声矩阵

来表示(1)式。

通过这种数据重组之后，可以认为矩阵

是低秩的，因此就可以利用低秩矩阵逼近的方法从

中恢复出

达到去噪的目的了。

求解(2)式最原始的方法定义为：

但是式(3)中的目标函数是非凸的，求解它是一个NP难问题，但是可以通过凸松弛转换为核范数最小化(NNM)，如下所示：

其中γ是固定正常数，||·||_F表示Frobenius范数，

表示矩阵

的核范数，||·||_*＝∑_i|λ_i(·)|，λ_i(·)表示矩阵的第i个奇异值。可以通过奇异值分解的方法求出(4)式的最优解。

其中Σ_ii表示矩阵Σ的对角线元素，

可以观察到NNM对每一个奇异值都做相同的阈值处理，这没有考虑到奇异值的实际物理意义，地震信号的有效信息主要是由大的奇异值表示，因此较大的奇异值与较小的奇异值相比应该要分配更小的阈值。因此使用加权核范数最小化(WNNM)，根据不同的奇异值分配不同的权重阈值，定义为：

其中||·||_w,*＝∑_i|w_iλ_i(·)|，w＝[w₁,w₂,…,w_n]，w_i是对应与奇异值λ_i(·)的非负权重值。w_i与的

关系定义为

其中c＞0是常数，m是构成矩阵

的块的数目，避免分母为0，让ε＝10^-16，从(8)中可以看出奇异值越大，其对应的w_i越小。虽然当权重矩阵w是非降序排列时(7)式是非凸的，存在一个局部极小值点是它的最优解。(7)式的求解方式与(4)式相似，将(5)式中的γ换成w_i即可，即

其它的求解步骤与求解(4)式相同。

在求解w_i时存在一个问题，就是无法事先知道纯净地震信号的奇异值

但是可以通过含噪地震信号的奇异值和含噪地震信号中的噪声方差估计出来，定义为：

其中

通过计算是可以得到的，所以现在的问题是如何获得σ，我们的工作也是围绕σ展开的，从公式(8)、(9)、(10)可以看出σ²的准确性会直接影响到

估计的准确性，同时就会影响到w_i的值，会影响到信号去噪的效果，因此我们引入了基于几何纹理的噪声估计方法，根据含噪地震信号不同纹理块的几何结构直接估计出噪声标准差，无需手动设置噪声水平初值，也不用考虑前次迭代去噪的效果。

基于几何纹理的噪声估计最先是应用到图像中，这种估计方法在图像中有较好的实际应用效果，估计出来的噪声标准差与实际值非常接近。于是本发明将这种估计方法应用到含噪地震信号中。噪声水平可以根据地震纹理块中的弱纹理块估计为

其中ω_min(·)表示矩阵的最小特征值，W表示挑选出来弱纹理块，cov(·)表示协方差矩阵。因此将基于几何纹理的噪声估计方法应用到地震信号中的关键问题是如何从地震纹理块中挑选出弱纹理块。首先地震纹理块的梯度协方差矩阵

可以定义为：

其中p_j表示矩阵

中的纹理块矩阵，D_h和D_v分别表示水平和垂直微分算子。

的较大特征值往往反映了p_j更多的纹理信息。理论上来说，重组地震信号

中的含噪平坦地震块

可以分解为完美平坦的地震块

和噪声块

如下所示：

完美平坦的地震数据块的梯度是0，所以含噪平坦地震块

的梯度协方差矩阵的期望值为

其中两个对角线元素的统计特性是相同的，两个对角线元素

(其中Ω＝h或v)的分布可以由伽马分布来近似，该伽马分布的形状参数α和尺度参数β分别为

其中N表示矩阵

中元素的个数，σ表示噪声标准差，

代表矩阵

的迹。

为了挑选出地震纹理块中的弱纹理块，采用零假设检验的方法，零假设内容为“所给的块为含有随机噪声的平坦块”，当所给地震纹理块的梯度协方差矩阵的最大特征值小于某个阈值ρ时，认为假设成立，反之不成立。阈值ρ定义为

ρ＝F^-1(v,α,β)， (16)

其中F^-1表示逆伽马累积分布函数，v表示显著性水平，必须人为设置，一般设为0.95～0.99。从(15)、(16)式可以看出噪声标准差σ也是计算阈值ρ必不可少的参数。

本发明采用迭代机制来估计最后的噪声水平，将其与地震信号结合，归纳如下：首先由

中的所有纹理块产生的协方差矩阵通过(11)式估计出噪声标准差初值

并根据

通过(16)式计算出初始阈值ρ⁽⁰⁾，根据阈值ρ⁽⁰⁾通过零假设检验即可挑出弱纹理块W⁽¹⁾，再根据(11)式可以估计出

重复迭代，直到迭代出的σ_n与前一次迭代的结果不超过10^-6，认为噪声标准差σ_n趋于稳定，最后的σ_n就是最真实的噪声水平。

低秩矩阵逼近是一个迭代去噪的过程，因此在迭代过程中不可避免地会导致一些地震信号的丢失。考虑到地震信号奇异值的物理意义，一般高能量有效信号的奇异值都比较大，但是还有一些能量比较低的信号的奇异值相对较小，迭代过程中通过(9)式的阈值处理后，有可能会将这部分的奇异值变成0或非常接近0，就会导致这部分信号的丢失。为了减少这种丢失，在每次迭代前，用阈值τ将信号的较大奇异值和较小奇异值分开，即奇异值小于阈值τ的为较小奇异值，奇异值大于阈值τ的为较大奇异值，如下所示：

其中

和

分别是仅保留

中较大奇异值和较小奇异值的矩阵。(7)式的最优解就可以写为

然后，通过合并所有的块矩阵

和

获得两幅地震信号记录

和

和

的差值可以反映同时包含在

中但是具有不同能量水平的边缘和纹理地震信号，在迭代正则项中加上高、低能量地震信号的差值来减少迭代去噪过程中一些不必要的信号丢失。

为了更清楚的描述，我们将上述低秩矩阵逼近的算法总结为算法1，具体包括：

1)、输入含噪地震信号Y，设置参数δ，η，L，τ；

2)、初始化

3)、设置k＝1:L；

4)、对k＝1:L做迭代正则化：

5)、在第k次迭代处理后的含噪地震信号Y^k中设置搜索范围，选择主要地震纹理块

后，选择与主要地震纹理块相似的块构成含噪地震信号矩阵

6)、用(11)式估计噪声标准差

7)、奇异值分解

8)、用(10)式估计奇异值

9)、用(8)式计算权重向量w_i；

10)、用(9)式获得奇异值矩阵

11)、用阈值τ和(17)式分解奇异值

12)、用(18)式获得估计值

和

13)、利用

和

来获得

和

14)、得出去噪的地震信号

15)、将上述步骤迭代，直到获得的值趋于平稳，即为最大去噪结果的地震信号。

在一种实施例中，在分块时，设置地震纹理块的大小为2×2，也就是一个地震纹理块的矩阵大小为2×2的，总共分了50589个地震纹理块，然后将每个地震纹理块变成一列数，即每个地震纹理块都变成4×1的，然后将所有的地震纹理块合成一个矩阵，变成4×50589个的大矩阵，由于矩阵数据太多无法全列出来，以四个地震纹理块为例

地震纹理块1：

地震纹理块2：

地震纹理块3：

地震纹理块4：

合并后构成含噪地震信号矩阵：

本发明为运用迭代正则化方法确定噪声水平，第一次运用迭代正则化时将所有的地震纹理块均认为是弱纹理块，因此不需要计算每个地震纹理块的梯度协方差矩阵，设置阈值ρ＝INF，含噪地震信号矩阵的协方差矩阵，为：

含噪地震信号矩阵的协方差矩阵的特征值为：

λ₁＝3.99648244928769；

λ₂＝4.00178249846769；

λ₃＝4.04493252716955；

λ₄＝4.93125003670187；

噪声方差就为最小特征值，即噪声方差为λ₁＝3.99648244928769，噪声标准差为

继续进行迭代正则化，直到迭代出的噪声标准差与上一次迭代出的噪声标准差波动范围不超过10^-6，噪声标准差趋于稳定，即为真实的噪声水平。

为了验证本发明所提方法的有效性，本发明首先对合成地震数据进行测试。如图1所示是一个合成的地震模型，总共有7个有效事件，这些有效事件是由主频为30HZ和25HZ的雷克子波产生。然后在合成的地震数据中加入实际沙漠地区随机噪声，让合成记录的信噪比为-5.1402dB，如图2所示。在处理这幅记录的时候，设置搜索窗的大小为40×40，地震纹理块的大小为13×13，这个设置只适合于处理这个噪声水平下的记录，在处理不同噪声水平的合成记录时，搜索窗和地震块的大小都不一样，噪声水平较低时，将它们设置的小一些，噪声水平较高时，设置的大一些。同时设置参数δ＝0.1，η＝0.4，τ＝0.5来处理所有信噪比下的合成记录。

在对比实验中，小波变换，f-x反褶积和WNNM被用来处理相同的合成记录来和本发明所提的方法进行比较，将每种对比方法的参数调到使去噪结果最好。如图3、图4、图5所示分别是这三种方法对应的处理结果，图6是本发明方法处理的结果。为了更全面地观察和比较这四种方法的去噪结果，本发明又做了单道的分析。为此抽取了合成记录的第3道进行对比，如图7、图8、图9、图10、图11、图12所示分别是纯净信号，含噪信号，小波变换，f-x反褶积，WNNM和本发明方法的第三道单道去噪效果图。然后本发明又利用FK谱图对这四种方法的去噪结果在频域里进行了分析，如图13～图18分别代表纯净信号，含噪信号，小波变换去噪结果，f-x反褶积去噪结果，WNNM去噪结果和本发明方法去噪结果的FK谱图。

经过以上全方位多角度的分析后，了解到小波变换的噪声压制不彻底，有效信号也没有很好的恢复；f-x反褶积的噪声压制不是很彻底，与有效信号同频带的噪声都保留下来，也存在一部分信号的丢失；WNNM虽然噪声压制的比较好，但是有效信号的恢复情况并不理想；本发明方法压制噪声最为彻底，有效信号的恢复情况也是最好的，处理结果与理想信号最为接近。

接下来，用这四种方法处理了加入不同水平的实际沙漠地区随机噪声的合成记录，通过计算去噪结果的信噪比和均方误差来定量分析本发明方法的有效性，信噪比(SNR)和均方误差(MSE)的计算公式如下：

其中，X是理想信号，

是处理后的信号。

这四种方法处理结果的信噪比和均方误差分别在表.1，表.2中给出。信噪比越高，均方误差越低表示处理结果最好，从两个表中可以看出，本发明方法处理后的信噪比一直是最高的，均方误差一直是最小的。通过定量分析也可以得出本发明方法处理结果是最好的。

表一不同方法去噪结果的信噪比

表二不同方法去噪结果的均方误差

处理前	小波变换	f-x反褶积	WNNM	本发明
					SNR(dB)	SNR(dB)	SNR(dB)	SNR(dB)	SNR(dB)
1.0935	0.0479	0.0412	0.0332	0.0131
					-1.0050	0.0699	0.0318	0.0222	0.0100
-3.0694	0.1046	0.0412	0.0332	0.0131
					-5.1402	0.1604	0.0495	0.0557	0.0211
-7.0256	0.2402	0.0754	0.1156	0.0499
					-9.0900	0.3781	0.1082	0.1880	0.0877

为了验证本发明所提方法的实际使用性，将其用于处理实际沙漠地震数据，图19为某沙漠地区的地震共炮点记录，该记录总共有127道，采样频率为500Hz，可以看到该记录中两侧的同相轴几乎被随机噪声淹没，不易辨别。同样，小波变换、f-x反褶积、WNNM和本发明方法被应用于处理该记录，这四种方法处理结果分别对应于图20～图23。在处理该记录时，设置本发明方法中搜索窗大小为30×30，地震纹理块的大小为7×7，同时将WNNM中的参数与本发明方法中的参数设为一样。

从图20～图23中可以看出小波变换的低频噪声压制效果并不太好，而且还有效信号发生了畸变。f-x虽然整体对低频噪声有一定的压制效果，但是可以看到在该记录两侧噪声水平较高的部分并没有将有效信号很好的恢复出来。WNNM对记录深部地区的低频噪声的压制并不是太好，本发明方法整体上对低频噪声的压制较为理想，而且可以看出在记录两测随机噪声较强的部分，本发明方法仍然可以将有效信号很好的恢复出来。

为了更清晰的评价这四种方法去噪结果的质量，将该记录的的第77-110道，第1000-1900个采样点之间的记录局部放大如图24所示，图24～图28分别对应图19～图23的局部放大图。通过比较局部放大图，可以更明显的看出本发明方法与其他三种方法相比可以更好的压制沙漠地区低频噪声，同时可以很好的恢复出有效信号。

低秩矩阵逼近主要利用地震信号在时空域的非局部相似性，而不需要考虑信号的频率，因此将其应用到沙漠地震低频噪声压制中会有一定的效果。考虑到权重阈值与噪声标准差有着直接的关系，本发明利用基于几何纹理的噪声估计方法可以获得比较接近真实的噪声标准差，从而获得更加精确的权重阈值，同时考虑到不同大小的奇异值也能代表不同能量的有效信号，因此利用截断奇异值来减少迭代过程中有效信号的丢失。通过与小波变换、f-x反褶积、WNNM的对比说明了本发明方法无论是对沙漠随机噪声的压制还是有效信号的恢复都有最好的效果。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (8)

1.一种基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，其特征在于，包括以下步骤：

步骤一、获取含噪地震信号Y，并且对所述含噪地震信号Y做k＝1:L的正则化处理：

其中，Y^k表示第k次迭代处理后的含噪地震信号，

表示第k次迭代处理后的去噪地震信号，

表示第k次迭代处理后的较大奇异值的去噪地震信号，

表示第k次迭代处理后的较小奇异值的去噪地震信号，δ，η表示参数，τ表示奇异值的阈值，并且设置初始值Y^k＝Y，

步骤二、在含噪地震信号中设置搜索范围，选择主要地震纹理块矩阵，再选择与所述主要地震纹理块矩阵相似的地震纹理块矩阵构成含噪地震信号矩阵；

步骤三、所述地震纹理块矩阵的噪声标准差满足：

其中，

表示估计的噪声标准差，ω_min(·)表示矩阵的最小特征值，W表示挑选出来弱纹理块，cov(·)表示协方差矩阵；

步骤四、将所述含噪地震信号矩阵中的所有地震纹理块矩阵获得噪声标准差初值，并且设置初始阈值，当所述地震纹理块的梯度协方差矩阵的最大特征值小于初始阈值时，所述地震纹理块为弱纹理块；

其中，所述地震纹理块的梯度协方差矩阵满足：

其中，

表示地震纹理块的梯度协方差矩阵，p_j表示含噪地震信号矩阵

中的纹理块矩阵，D_h和D_v分别表示水平和垂直微分算子；

初始阈值为人为设置，在第k次迭代处理后，阈值满足：

ρ＝F^-1(v,α,β)；

其中，F^-1表示逆伽马累积分布函数，v表示显著性水平，α表示伽马分布的形状参数，β表示伽马分布的尺度参数；

步骤五、根据所述弱纹理块的噪声标准差得到纯净地震信号的奇异值：

其中，

表示纯净地震信号矩阵的奇异值，

表示含噪地震信号矩阵的奇异值，

表示所述弱纹理块的噪声标准差，m是构成所述含噪地震信号矩阵的块的数目；

步骤六、所述奇异值分解得到奇异值矩阵，根据阈值将较大奇异值和较小奇异值分开并保留；

步骤七、得出去噪的地震信号。

2.根据权利要求1所述的基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，其特征在于，在所述步骤二中，所述含噪地震信号矩阵中包括含噪平坦地震块，其分解为完美平坦地震块和噪声块；

其中，所述完美平坦地震块的梯度是0，所述含噪平坦地震块的梯度协方差矩阵的期望值为

其中，两个对角线元素

的统计特性是相同的，其中Ω＝h或v，

表示含噪平坦地震块的梯度协方差矩阵，

表示完美平坦地震块的梯度协方差矩阵，

表示噪声块的梯度协方差矩阵，

表示噪声块。

3.根据权利要求2所述的基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，其特征在于，在所述步骤四中，所述伽马分布的形状参数α和尺度参数β分别为：

其中，α表示伽马分布的形状参数，β表示伽马分布的尺度参数，N表示噪声块中元素的个数，

表示噪声标准差，

代表矩阵

的迹。

4.根据权利要求1所述的基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，其特征在于，在所述步骤六中，所述奇异值的分解过程为：

其中，

表示含噪地震信号矩阵，U和V分别表示左奇异值向量矩阵和右奇异值向量矩阵，

表示第k次迭代处理后的纯净信号矩阵，w_i表示对应与奇异值的非负权重值，

Σ表示奇异值对角矩阵，

表示奇异值矩阵，

表示奇异值矩阵的元素，Σ_ii表示奇异值对角矩阵中的元素，并且所述奇异值矩阵中的元素满足：

其中，Σ_ii表示矩阵的元素，

表示含噪地震信号矩阵的奇异值。

5.根据权利要求4所述的基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，其特征在于，所述奇异值分为较大奇异值矩阵和较小奇异值矩阵，如下所示：

其中，

表示仅保留

中较大奇异值的矩阵，

表示仅保留

中较小奇异值的矩阵，并且设置阈值τ，奇异值小于所述阈值为较小奇异值的矩阵，奇异值大于所述阈值为较大奇异值的矩阵。

6.根据权利要求5所述的基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，其特征在于，通过所述较大奇异值矩阵和所述较小奇异值矩阵得到所述纯净信号矩阵：

其中，

表示第k次迭代处理后的纯净信号矩阵，

表示保留纯净信号矩阵中较高奇异值获得的矩阵，

表示保留纯净信号矩阵中较低奇异值获得的矩阵。

7.根据权利要求6所述的基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，其特征在于，在所述步骤七中，通过合并所述保留纯净信号矩阵中较低奇异值获得的矩阵和所述保留纯净信号矩阵中较高奇异值获得的矩阵得到保留较大奇异值和保留较小奇异值的两幅地震信号，将所述两幅地震信号合并得到第k次迭代处理后的纯净地震信号。

8.根据权利要求7所述的基于几何纹理噪声估计的低秩矩阵逼近的沙漠地震噪声压制方法，其特征在于，在所述步骤七中，将所述两幅地震信号和所述第k次迭代处理后的纯净地震信号经过所述步骤一中的正则化处理，得到去噪的地震信号。

Images