一种地震数据降噪方法

摘要

本发明公开一种地震数据降噪方法，针对地震数据降噪问题，本发明包括：S1、采用稳健Huber范数替换目标函数的误差项，得到基于Huber范数的目标函数；S2、对经步骤得到的目标函数引入schatter p‑范数，得到基于schatter p‑范数的目标函数；S3、对经步骤S2得到的目标函数进行求解，得到Hankel矩阵的稳健估计；实现在压制随机噪声的同时又能有效去除规则噪声。

权利要求

1.一种地震数据降噪方法，其特征在于，包括：S1、采用稳健Huber范数替换目标函数的误差项，得到基于Huber范数的目标函数；

S2、对经步骤S1得到的目标函数引入schatter p-范数，得到基于schatter p-范数的目标函数；步骤S2基于schatter p-范数的目标函数表达式为：其中，下标Sp表示schatter p-范数，X表示地震数据矩阵，D表示外加噪声，γ为第一正则化系数；

S3、对经步骤S2得到的目标函数进行求解，得到Hankel矩阵的稳健估计。

2.根据权利要求1所述的一种地震数据降噪方法，其特征在于，步骤S1具体包括：S11、地震数据矩阵X的截断SVD的模型表达式如下：s.t.rank(D)＝K

其中，K表示为地震数据的秩，||·||F为Frobenius范数，argmin表示求函数取值最小时的自变量取值，rank(·)表示求矩阵的秩；

S12、将式(1)的表达式转化为正则项目标函数形式：其中，γ为第一正则化系数；

S13、对式(2)的误差项引入Huber范数，得到基于Huber范数的目标函数表达式为：其中，||·||H为Huber范数，X表示地震数据矩阵，D表示外加噪声，arg min(·)表示函数取值最小时的自变量取值，rank(·)表示求矩阵的秩。

3.根据权利要求2所述的一种地震数据降噪方法，其特征在于，步骤S3具体为：S31、将 变换为下式：

E和Z表示随机变量；

S32、采用ADM算法分别求解X、E、Z；

S33、根据步骤S32得到的X、E、Z，得到Hankel矩阵的稳健估计。

4.根据权利要求3所述的一种地震数据降噪方法，其特征在于，S32求解X包括以下分步骤：A1、初始化随机变量Λ，Σ，E和Z；

A2、开始while迭代循环；

A3、固定E和Z求解X的子问题表达式为：其中，λ表示第二正则化系数；

A4、根据步骤A3中表达式的解析解更新X；

X＝N

其中，

A5、使用Λ＝Λ+μ(E-X+D)更新Λ，使用Σ＝Σ+μ(X-Z)；使用μ＝ρμ更新μ；

A6、重复A4至A5，直到收敛。

5.根据权利要求3所述的一种地震数据降噪方法，其特征在于，S32求解E包括以下分步骤：B1、初始化随机变量Λ，Σ，E和Z；

B2、开始while迭代循环；

B3、固定X和Z求解E的子问题表达式为：其中，μ表示第三正则化系数，

B4、使用Λ＝Λ+μ(E-X+D)更新Λ，使用Σ＝Σ+μ(X-Z)；使用μ＝ρμ更新μ；

B5、重复B3至B4，直到收敛。

6.根据权利要求3所述的一种地震数据降噪方法，其特征在于，S32求解E包括以下分步骤：C1、初始化随机变量Λ，Σ，E和Z；

C2、开始while迭代循环；

C3、固定X和E求解Z的子问题表达式为：其中，μ表示第三正则化系数；

C4、使用Λ＝Λ+μ(E-X+D)更新Λ，使用Σ＝Σ+μ(X-Z)；使用μ＝ρμ更新μ；

C5、重复C3至C4，直到收敛。

说明书

一种地震数据降噪方法

技术领域

[0001] 本发明属于地震数据处理领域，特别涉及一种地震数据降噪技术。

背景技术

[0002] 地震勘探是开发矿藏、石油和天然气等资源的重要手段。在勘探的过程当中，数据会不可避免地会受到各种因素的影响，导致收集到的数据中会存在多种类型的噪声：比如地震信号在传播过程中幅度和频率逐渐衰减，收集数据过程中噪声始终存在，这会导致真实信号的信噪比下降，以至于有用信号湮没在噪声中，进而影响我们对地质信息的分析。因此我们需要一种有效地降噪手段来提高地震资料的质量，提高地震勘探的可靠性。

[0003] 现有研究工作主要集成在随机噪声的降噪上，这方面有大量的研究工作，总结起来可以分成两类：(1)时域滤波方法：比如经典结构导向滤波，中值滤波、保边滤波等；(1)变换域方法：包括常用的tau-p变换、拉东变换、f-x域滤波、eigenimage滤波器、EMD、基于小波变换分解，radial trace域和curvelet变换进行去除，变换域方法由于在变换域系数上操作，逆变换回去之后会产生明显的副作用，比如滤波后的道集记录上会出现明显的蚯蚓化现象等。更重要的是以上两个框架无法很好去除规则噪声，所以在实际地震资料处理时，需要建立一个处理流程来进行高斯和非高斯噪声压制：先用以上方法去除高斯白噪声，然后再用专门方法除规则噪声，这样做一方面非常费时费力；另外一方面随机噪声的平滑效果会影响规则噪声的去除。

发明内容

[0004] 为了解决上述技术问题，本发明提出一种地震数据降噪方法，引入稳健Huber范数，提升非高斯去噪能力，同时保留对高斯噪声的降噪能力。

[0005] 本发明采用技术方案为：一种地震数据降噪方法，包括：

[0006] S1、采用稳健Huber范数替换目标函数的误差项，得到基于Huber范数的目标函数；

[0007] S2、对经步骤得到的目标函数引入schatter p-范数，得到基于schatter p-范数的目标函数；

[0008] S3、对经步骤S2得到的目标函数进行求解，得到Hankel矩阵的稳健估计。

[0009] 进一步地，步骤S1具体包括：

[0010] S11、地震数据矩阵X的截断SVD的模型表达式如下：

[0011]

[0012] s.t.rank(D)＝K

[0013] 其中，K表示为地震数据的秩，||·||F为Frobenius范数，argmin表示求函数取值最小时的自变量取值，rank(·)表示求矩阵的秩；

[0014] S12、将式(1)的表达式转化为正则项目标函数形式：

[0015]

[0016] 其中，γ为第一正则化系数；

[0017] S13、对式(2)的误差项引入Huber范数，得到基于Huber范数的目标函数表达式为：

[0018]

[0019] 其中，||·||H为lp-范数，X表示地震数据矩阵，D表示外加噪声，arg min(·)表示函数取值最小时的自变量取值，rank(·)表示求矩阵的秩。

[0020] 进一步地，步骤S2基于schatter p-范数的目标函数表达式为：

[0021]

[0022] 其中，下标Sp表示schatter p-范数。

[0023] 进一步地，步骤S3具体为：

[0024] S31、将 变换为下式：

[0025]

[0026] S32、采用ADM算法分别求解X、E、Z；

[0027] S33、根据步骤S32得到的X、E、Z，得到Hankel矩阵的稳健估计。

[0028] 进一步的，S32求解X包括以下分步骤：

[0029] A1、初始化随机变量Λ，Σ，E和Z；

[0030] A2、开始while迭代循环；

[0031] A3、固定E和Z求解X的子问题表达式为：

[0032]

[0033] 其中，λ表示第二正则化系数；

[0034] A4、根据步骤A3中表达式的解析解更新X；

[0035]

[0036] 其中，

[0037] A5、使用Λ＝Λ+μ(E-X+D)更新Λ，使用Σ＝Σ+μ(X-Z)；使用μ＝ρμ更新μ；

[0038] A6、重复A4至A5，直到收敛。

[0039] 进一步的，S32求解E包括以下分步骤：

[0040] B1、初始化随机变量Λ，Σ，E和Z；

[0041] B2、开始while迭代循环；

[0042] B3、固定X和Z求解E的子问题表达式为：

[0043]

[0044] 其中，μ表示第三正则化系数，

[0045] B4、使用Λ＝Λ+μ(E-X+D)更新Λ，使用Σ＝Σ+μ(X-Z)；使用μ＝ρμ更新μ；

[0046] B5、重复B3至B4，直到收敛。

[0047] 进一步的，S32求解E包括以下分步骤：

[0048] C1、初始化随机变量Λ，Σ，E和Z；

[0049] C2、开始while迭代循环；

[0050] C3、固定X和E求解Z的子问题表达式为：

[0051]

[0052] 其中，μ表示第三正则化系数；

[0053] C4、使用Λ＝Λ+μ(E-X+D)更新Λ，使用Σ＝Σ+μ(X-Z)；使用μ＝ρμ更新μ；

[0054] C5、重复C3至C4，直到收敛。

[0055] 本发明的有益效果：本发明的方法，通过基于Huber范数的目标函数代替原来截断SVD所采用的F范数目标函数的误差项，然后schattern p-范数来对目标函数的秩最小化项进行正则化，最后迭代ADM估计法用于近似求解基于schattern p-范数的目标函数，从而获得Hankel矩阵的稳健估计；实现在压制随机噪声的同时又能有效去除规则噪声。

附图说明

[0056] 图1为本发明实施例提供的地震数据低秩特性示意图；

[0057] 图2为本发明实施例提供的Huber范数示意图；

[0058] 图3为本发明实施例提供的方案流程图；

[0059] 图4为本发明实施例提供的原始叠后地震剖面图一；

[0060] 图5为本发明实施例提供的降噪后的叠后剖面图一；

[0061] 图6为本发明实施例提供的经本发明方法剖面一去除的干扰；

[0062] 图7为本发明实施例提供的原始叠后地震剖面图二；

[0063] 图8为本发明实施例提供的降噪后的叠后剖面二；

[0064] 图9为本发明实施例提供的经本发明方法剖面二去除的干扰。

具体实施方式

[0065] 为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

[0066] 本发明目的在于提出一种地震数据降噪方法能够在压制随机噪声的同时又能有效去除规则噪声。由于地震数据具有两个特性：(1)如图1所示，干净的地震数据具有低秩特性，而任何噪声都将会破坏地震数据的相关性，从而破坏地震数据的低秩结构；(2)地震随机噪声为带限的高斯白噪声，也可近似为高斯分布的噪声，规则噪声可笼统视为非高斯分布的噪声。基于以上两个特性，本发明将随机和规则噪声同步去除方法问题创新转变为从含噪的地震数据中提取低秩结构的技术问题。

[0067] 本发明方法解决以上技术问题的具体思路如下：(1)schatten p-范数正则化提升低秩结构稳定性：由于规则噪声等非高斯噪声表现为集体离群值，如果采用常规核范数的方法无法克服离群值影响，所以本发明创新地使用schatten p-范数代替核范数来消除离群值影响；(2)Huber范数从波形角度稳健去除非高斯噪声：如图2所示Huber范数为L1范数和L2范数融合的范数，兼具L1范数和L2范数的优点，在迭代过程中从波形上消除高斯和非高斯噪声影响；(3)Hankel矩阵提升地震数据低秩结构：Hankel矩阵等结构化矩阵常常用来提升矩阵的秩，矩阵的秩越低则恢复越好，所以Hankel矩阵在本发明中是数据的表达形式。

[0068] 本实施例中简要介绍基于Hankel矩阵的地震数据降噪基本原理，基于Hankel矩阵的地震噪声衰减和地震数据重建的详细内容可见文献“Oropeza,V.,and M.D.Sacchi,2011,Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis:Geophysics,76,no.3,V25–V32,doi:10.1190/1.3552706”。
基于Hankel矩阵的地震数据降噪具体主要可以分成两个步骤：

[0069] (1)基于滑动时窗的地震数据Hankel变换

[0070] 本实施例主要以降秩滤波的二维(t-x)地震数据中的实现为例进行说明，但是降秩滤波也可以在广泛应用到三维和五维的地震数据中。降秩滤波主要利用地震数据的低秩性，通过Hankel变换可以进一步降低理想地震信号的低秩性，同时升高含噪地震资料的秩。Hankel变换主要涉及两个步骤：

[0071] 步骤一：对地震数据进行傅里叶变换：小窗口内中的地震数据可以通过频率-空间域内平面波的叠加获得：

[0072]

[0073] 其中， j＝1,2,...,N是空间轴上的地震道集索引值，ω表示时间频率。假设地震数据由具有不同射线参数Pk的K个线性事件组成。这里，Ak(ω)表示第k平面波的复振幅，Δx表示地震剖面之间的空间间隔。

[0074] 步骤二：Hankel变换：降秩滤波将单频D(ω)＝(D1(ω)，D2(ω)，...，DN(ω))T嵌入信号来构建如下Hankel矩阵：

[0075]

[0076] 其中 表示Hankel算子。

[0077] 为了便于计算，本实施例中选择 使Hankel矩阵逼近为方阵，省略符号ω，并且对所有频率进行分析。对于K平面波的叠加，可得rank
(MX)＝K。

[0078] (2)基于Hankel域的低秩逼近

[0079] 由于外加噪声D将会增加矩阵X的秩，降秩滤波器通过降低秩，从而达到降噪的目的。降秩滤波可以表示成如下：

[0080]

[0081] 其中， 是反对角平均算子， 是降秩操作符，它通过秩K矩阵逼近M， 是Hankel算子。运算符 通过平均对角线将Hankel形式转换成向量。对于多维信号则采用块Hankel矩阵和阻尼对角平均算子。降秩操作 可以通过截断的SVD，随机SVD或采用Lanczos双对角化方法和高效矩阵向量乘法的快速算法来实现，快速算法采用快速傅里叶变换。

[0082] 基于上述原理，本发明的技术方案为：一种地震数据降噪方法，如图1所示，包括以下步骤：

[0083] S1、采用稳健Huber范数替换目标函数的误差项，得到基于Huber范数的目标函数；具体包括以下分步骤：

[0084] S11、地震数据降噪的截断低秩分解模型

[0085] 干净的地震数据经过SVD分解后发现：地震信号的大部分能量集中在少量几个特征值上面，大部分能量集中在前面几个特征值上，这样的特性就称为地震数据的低秩特性。更进一步推广，有效的地震数据的能量分布在低秩的部分成分上，噪声则分布在除了低秩以外的成分上。因此采用截断SVD方法就可对地震数据进行降噪。

[0086] 更具体来说，地震数据矩阵X的截断SVD的模型可以表达如下：

[0087]

[0088] s.t.rank(D)＝K

[0089] 其中，K表示为地震数据的秩，||·||F为Frobenius范数， 为矩阵 公式(1)中最优化问题本质上是基于Frobenius范数的凸优化。

[0090] S12、将公式(1)转化为正则项目标函数形式：

[0091]

[0092] 其中γ为正则化系数。

[0093] S13、为了进一步提升对非高斯噪声的压噪能力，本发明对式(2)误差项引入稳健范数以提升对非高斯的去噪能力，同时也需要保留对高斯噪声的降噪能力。则将将公式(2)转化为基于稳健Huber范数的矩阵分解方式：

[0094]

[0095] 其中，||·||H为lp-范数，由于Huber范数也是凸函数，所以Huber范数的引入不会改变目标函数(3)的性质，其依旧是凸函数。

[0096] S2、对经步骤得到的目标函数引入schatter p-范数，得到基于schatter p-范数的目标函数；具体过程为：

[0097] 由于实际地震数据存在高斯噪声和非高斯噪声，因此除了在目标函数误差项需要考虑非高斯和高斯噪声的情况，也需要在低秩逼近时候考虑非高斯的情况，schatten 1-范数能够比普通秩最小化方法能够更好地低秩逼近，而且对离群值稳健。矩阵的schatten 1-范数即为公式(4)所描述的核范数优化模型，将其扩展到schatten p-范数(0＜p＜∞)，则有如下优化模型：

[0098]

[0099] 当p＝1时，schatter p-范数即为核范数。如果定义σ0＝0，则可以将schattern p-范数扩展到p＝0的情况，即为秩函数。引入了schatter p-范数之后，公式(3)目标函数可变为：

[0100]

[0101] S3、对经步骤S2得到的目标函数进行求解，得到Hankel矩阵的稳健估计。具体过程为：

[0102] 为了求解公式(3)的问题，本发明选择使用一种新的Alternating Direction Method(ADM)(具体可参考：Lin Z,Liu R,Su Z.Linearized alternating direction method with adaptive penalty for low-rank representation[C]//Advances in neural information processing systems.2011:612-620.)。为了ADM求解方便，将公式(5)表示为：

[0103]

[0104] 根据ADM算法，可以将公式(5)和(6)分解成以下几个子问题求解。

[0105] a、当固定E和Z则求解以下子问题：

[0106]

[0107] 其中 和 对公式(7)的子问题，容易得到以下解析解：

[0108]

[0109] b、当固定X,Z我们可以交替求解如下子问题：

[0110]

[0111] 其中，

[0112] c、当固定X,E我们可以交替求解另外一个子问题：

[0113]

[0114] 其中 G为X的近似估计值。公式(9)和(10)属于经典最优化问题求解，可以采用现有成熟算法求解即可，比如共轭梯度算法。

[0115] 综上，以求解X为例，步骤S3的求解流程为：

[0116] A1、初始化随机变量Λ，Σ，E和Z；

[0117] A2、开始while迭代循环；

[0118] A3、使用公式(8)更新X；

[0119] A4、使用Λ＝Λ+μ(E-X+D)更新Λ，使用Σ＝Σ+μ(X-Z)；

[0120] A5、使用μ＝ρμ更新μ；本领域技术人员应知这里的ρ的取值根据实际中的精度要求设置。

[0121] A6、重复A3至A5，直到收敛。收敛条件为前后两次迭代目标函数值小于预设门限时，表示达到收敛条件。设定门限根据实际应用中的精度要求计算，一般为了减少迭代次数，设置为0.001也是可行的。

[0122] 将上述求解流程中的步骤S33替换为使用公式(9)更新E；则得到了求解E的流程；将上述求解流程中的步骤S33替换为使用公式(10)更新Z；则得到了求解Z的流程。

[0123] 图4至图6为本实施例提供的原始叠后地震剖面图一的处理过程；图4为原始叠后地震剖面图一，图5为降噪后的叠后剖面一，图6为经本发明方法剖面一去除的干扰；图7只图9为本实施例提供的原始叠后地震剖面图二的处理过程；图7为原始叠后地震剖面图二，图8为降噪后的叠后剖面二，图9为经本发明方法剖面二去除的干扰；可以看出，实际去噪结果符合理论假设，并且实际效果非常好，降噪效果很明显。图4至图9中的Rows表示行，Columns表示列。

[0124] 本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。