核范数与广义全变差联合约束的地震随机噪声压制方法与流程

本发明属于石油地球物理勘探技术领域，具体技术为核范数与广义全变差联合约束的地震随机噪声压制方法，可以有效地压制地震数据中的随机噪声，提高地震资料信噪比，同时保持地震数据的局部光滑性，提高地震资料保真度。

背景技术：

石油勘探发展到今天，勘探开发逐渐转向复杂型、隐蔽型、深层和非常规油气藏，对地震资料质量要求越来越高。受复杂采集环境影响，野外采集的地震资料常含有很强的随机噪声，严重的甚至将有效信息完全淹没。地震随机噪声压制是地震资料处理中需要面对的重要问题。

地震随机噪声压制方法有很多，这些方法通常不直接在时间域内对地震资料进行去噪处理，而是在某种变换域内对其变换系数进行相应处理，之后进行反变换来达到去噪目的。

用于地震数据随机噪声压制的变换有很多种，如fourier变换、radon变换、wavelet变换、curvelet变换、seislet变换、ridgelet变换、shearlet变换等。fourier变换是一种全局变换，基于fourier变换的阈值去噪方法可能会产生吉布斯现象；radon变换虽然能较好地识别直线特征，但地震波前通常并不是直线；wavelet变换在地震去噪方面表现出一定的优越性，但它主要用于表示各向同性的点奇异特征，对于地震波前这种曲线奇异特征并不能进行有效的稀疏表示，因此基于wavelet变换的地震数据压缩和去噪等方法，都不可避免地在边界和波前奇异附近存在不光滑和模糊现象。近年来发展起来的curvelet变换、seislet变换、ridgelet变换、shearlet变换等则弥补了传统小波变换的不足，具有良好的多尺度性、多方向性和各向异性，更加适合表示地震波前特征，该变换一经提出，便在地震数据处理领域获得广泛应用。

基于变换域的地震随机噪声压制方法都可统一在信号稀疏性表示理论中，并在地震数据处理领域引起广泛关注。例如curvelet变换得到广泛应用正是由于其具有很好的稀疏表示能力。但在许多情况下，地震数据的稀疏性未必能有效地表现地震数据的内部结构特征，特别是相邻地震道间的波形相似性。而基于相邻地震道波形相似性的噪声压制方法，如核范数约束的噪声压制等，很少考虑到地震数据的局部光滑性，使得随机噪声压制地震数据在地质体边缘部位产生强烈的“阶梯抖动”现象。

技术实现要素：

本发明的目的在于提供核范数与广义全变差联合约束的地震数据随机噪声压制方法，综合地震数据处理、图像处理、最优化理论、矩阵分析等多个学科，将图像处理中的广义全变差约束、矩阵分析中的核范数约束、最优化理论的分裂bergman迭代有机结合，适用于地震数据随机噪声压制处理等多个环节，有效提高地震资料的信噪比与保真度。

为达到以上技术目的，本发明提供以下技术方案。

核范数与广义全变差联合约束的地震随机噪声压制方法，其核心是地震数据矩阵核范数约束、广义全变差约束、分裂bergman迭代求解。该方法依次包括以下步骤：

(1)输入三维地震数据体，将地震数据按inline线进行输入，再对每条inline线地震数据按时间采样点与cmp号排列成矩阵，其中时间采样点按行排列，cmp按列排列；

(2)构建目标函数；

(3)(3)采用分裂bergman迭代策略，将目标函数分解为两个子问题进行迭代，针对具体子问题采用相应算法进行求解；

(4)采用奇异值软阈值化算法求解子问题一，求取辅助变量矩阵；

(5)采用梯度投影算法求解子问题二，求取地震数据矩阵；

(6)判断bergman迭代是否停止，若停止迭代则输出随机噪声压制后的地震数据矩阵；

(7)对每条inline线地震数据都进行随机噪声压制处理，得到整个随机噪声压制后的地震数据体。

其中，所述的步骤(2)具体过程为：构建核范数与广义全变差联合约束随机噪声压制方法的目标函数；

所述的步骤(4)具体过程为：将地震数据矩阵视为已知量，输入核范数约束系数与辅助参数，采用奇异值软阈值化算法求解子问题一，求取辅助变量矩阵；

所述的步骤(5)具体过程为：输入广义全变差约束系数、上一次bergman迭代得到的地震数据矩阵，以及步骤(4)得到的辅助变量矩阵，采用梯度投影算法求解子问题二，求取更新后的地震数据矩阵。

本发明提供的核范数与广义全变差联合约束的地震随机噪声压制方法，优点在于：有机结合核范数和广义全变差，对地震随机噪声的压制处理进行联合约束；通过多道地震数据矩阵的核范数约束有效地表征相邻地震道之间的波形相似性，通过广义全变差约束约束地震数据的局部光滑性，从而更有效地压制地震数据中的随机噪声，提高地震资料信噪比，并保持地震数据的局部光滑性，提高地震资料保真度。

附图说明

图1为本发明的流程图。

具体实施方式

结合实施例说明本发明的具体技术方案。

核范数与广义全变差联合约束的地震随机噪声压制方法，如图1的流程：

步骤一：输入含躁地震数据体。将地震数据按inline线进行输入，再对每条inline线地震数据按时间采样点与cmp号排列成矩阵d，其中时间采样点按行排列，cmp按列排列。

步骤二：构建核范数与广义全变差联合约束的地震随机噪声压制方法的目标函数，即

其中，f(x)即为目标函数，d为原始含躁地震数据矩阵，x为噪声压制后地震数据矩阵，为矩阵x-d的frobenius范数，||x||*为x的核范数，||x||gtv为x的广义全变差，λ与μ分别为相应约束系数。

目标函数(1)中，矩阵x的核范数为其奇异值之和，即

其中λ(x)是矩阵x的奇异值组成的矢量，σi是矩阵x的第i个奇异值，||λ(x)||1为矢量λ(x)的1范数，而矩阵x的广义全变差为

其中为导数算子，为二阶导数算子，vec为矩阵矢量化转换算子。

因此，目标函数(1)可等价地写为

步骤三：采用分裂bergman迭代策略，将目标函数(4)分解为两个子问题迭代，针对具体子问题采用相应算法进行求解。

分裂bergman迭代策略的基本思想是引进一个辅助变量矩阵θ，将目标函数(4)改写为

其中f(x,θ)为改写后的目标函数，β为bergman迭代辅助参数。

由分裂bergman迭代原理可得，随着辅助参数β逐渐增大，目标函数(5)会渐近收敛于目标函数(4)。由此，将目标函数(4)分解为以下两个子问题进行迭代：

子问题一：将地震数据矩阵x视为已知量，求解下式(6)以求取辅助变量矩阵θ，

其中f(θ)为子问题一目标函数。

子问题二：将辅助变量矩阵θ视为已知量，求解下式(7)以求取迭代更新后的地震数据矩阵x，

其中上式中的f(x)为子问题二目标函数。

不失一般性，假设已进行到第l次bergman迭代。

步骤四：将地震数据矩阵x视为已知量，输入核范数约束系数λ与辅助参数β，采用奇异值软阈值化算法求解子问题一，求取辅助变量矩阵θ。一般地，在第一次bergman迭代时，可将地震数据矩阵x取为输入的原始地震数据矩阵d。算法的具体计算过程为：

首先，对地震数据矩阵x进行奇异值分解，得到

x＝uλv^t(8)

其中u与v分别为矩阵x的左右奇异矩阵，λ为奇异值矩阵，上标t为矩阵转置。

然后，对x的奇异值进行软阈值化，即

λτ＝sτ(λ)(9)

其中τ＝λ/β，sτ(λ)即为对λ进行软阈值化，λτ为奇异值软阈值化结果，计算公式为

最后，按下式(11)计算辅助变量矩阵θ，即

θ＝uλτv^t(11)

步骤五：输入广义全变差约束系数μ以及上一次bergman迭代得到的地震数据矩阵x，再将步骤四得到的辅助变量θ代入式(7)中，采用梯度投影算法求解子问题二，求取更新后的地震数据矩阵x。同样，在第一次bergman迭代时，可将矩阵x取为输入的原始地震数据矩阵d。算法的具体计算过程为：

1)令

i为单位矩阵；

2)计算矩阵a^ta的最大特征值ζ；

3)计算

k＝x+(1/ζ)a^t(k-ax)(14)

4)计算(w,q)＝l^*(x)，其中l^*是作用在x上的线性算子，它将实数矩阵集合r^n×m中的元素映射为实数矩阵对集合(r^(n-1)×m,r^n×(m-1))的元素，(w,q)为计算结果，计算公式为

wi,j＝xi,j-xi+1,j,1≤i≤n-1,1≤j≤m(15)

qi,j＝xi,j-xi,j+1,1≤i≤n,1≤j≤m-1(16)

其中wi,j、qi,j与xi,j分别为矩阵w、q、x的第i行第j列的分量。

5)计算

(p,t)＝pro{(w,q)+(ζ/8μ)l^*[k-(μ/ζ)l(w,q)]}(17)

其中l为作用在(w,q)上的线性算子，将实数矩阵对集合(r^(n-1)×m,r^n×(m-1))中的元素映射为实数矩阵集合r^n×m中的元素，计算公式为

l(w,q)i,j＝wi,j-wi-1,j+qi,j-qi,j-1,1≤i≤n,1≤j≤m(18)

pro为作用在一组实数矩阵对上的梯度投影算子，它将实数矩阵对(r^(n-1)×m,r^n×(m-1))中的元素映射为实数矩阵对(r^(n-1)×m,r^n×(m-1))中的元素，即

(p,t)＝pro(w,q)(19)

(p,t)为梯度投影的结果，计算公式为

同样，pi,j、ti,j分别为矩阵p、t的第i行第j列的分量。

6)计算

x＝k-(μ/ζ)l(p,t)(22)

步骤六：若满足下面的bergman迭代容限条件，则停止bergman迭代，输出地震数据矩阵x，

其中l为当前迭代次数，ε是bergman迭代容限参数，xl-1为上一次迭代得到的地震数据矩阵，xl为当前第l次迭代得到的地震数据矩阵；否则，令

β＝1.2β(24)

并转入步骤四，进行下一次迭代。

步骤七：对输入的每条inline线地震数据都进行上述随机噪声压制处理运算，得到整个随机噪声压制后的地震数据体。