基于全变分低秩逼近的三维张量随机加足迹噪声去噪方法

1.本发明属于地震数据处理领域，特别涉及一种随机加足迹噪声衰减技术。


背景技术：

2.采集足迹被明确地描述为发生在三维地震时间或水平切片上的周期性振幅伪影，这些伪影实质上是由与测量几何形状不可分割的振幅变化引起的，而不是地下的地质。这些足迹伪影显著地模糊了地震剖面的地质构造解释，尤其是通过地震属性分析或avo(振幅随偏移量变化)等方法进行的解释。一个值得注意的事实是，随机噪声也不可避免地伴随着足迹，出现在真实的三维野外数据中。
3.对于足迹噪声的去除，关键是遇到一个合适的变换域，将有价值的信号和足迹噪声彻底分离。最近，地震文献中记录了大量的变换域，它们大致分为三大类:频域、时频域和字典学习(dictionary learning,dl)。第一类将足迹作为空间周期噪声，在频域使用各种滤波技术去除足迹，包括经典波数滤波、自适应波数滤波、谱陷波滤波、双通频率波数滤波和波数陷波滤波。第二类使用时频变换来洞察时间和频率域的足迹，从而显著减少信号和足迹重叠，便于滤波。典型的时间-频率变换包括小波变换、三维复小波变换、平稳小波变换(stationarywavelettransform,swt)、基追踪和曲波变换。然而，滤波技术既能消除这些不良痕迹，又能使原始地震信号不受影响，这是一个矛盾。
4.与解析傅里叶变换和时频变换完全不同的是，最后一个类恰当地分配所学字典的原子(而不是之前的稀疏系数)来表示清晰的图像和足迹模式。一旦获得了表示足迹模式的原子，基于dl的去噪方法就可以对它们进行操作，对所有学习到的原子进行二进制模式分类或过滤任务重构足迹组件，然后从输入中减去足迹得到去噪结果。然而，尽管目前基于dl的方法在足迹去除方面取得了显著的进展，但该去除过程存在一个关键缺陷:作为典型的矩阵dl方法，它们必须将3-d变平，为了便于处理，将三维地震数据中隐藏的复杂高维图像结构自然湮灭，从而明显影响了基于dl方法的去噪性能。
5.对于随机噪声的衰减，地球物理学界进行了大量的研究工作，大致可以分为三个研究方向：空域滤波、稀疏表示(sparse representation,sr)和低秩近似(low rank approximation,lra)。第一个方向利用空间掩模中滤波系数和相关图像像素的乘积之和来产生去噪结果。典型的空域滤波包括非局部均值滤波、中值滤波、形态滤波和非线性扩散滤波。第二个方向是在给定的稀疏变换中，利用sr表示系数尽可能少的无噪声地震数据；但是高斯随机噪声不能用这种方法表示。有几种稀疏变换可以用于图像变换，如小波变换，curvelet变换，小波变换，shearlet变换，学习词典。最后一个方向对一个数据矩阵或张量采用矩阵lra等多种降阶方法实现地震图像去噪，基于hankel矩阵的lra方法，张量奇异值分解(tensor singular value decomposition,tsvd)，candecomp/parafac(cp)分解。当然，目前的随机噪声抑制技术也经常被用于消除足迹，如非线性扩散滤波和奇异值分解(singular value decomposition,svd)。但是，很明显，这些去噪方法并不能完全去除混合噪声，重建出高质量的地震图像，因为之前的方法都是强调随机噪声的抑制而不是足迹的
去除。


技术实现要素：

6.为解决上述技术问题，本发明提出一种基于全变分低秩逼近的三维张量随机加足迹噪声去噪方法，对三维地震数据进行随机噪声和足迹去除。
7.本发明采用的技术方案为：一种基于全变分低秩逼近的三维张量随机加足迹噪声去噪方法，包括：
8.s1、待处理的噪声图像张量记为表达式为：
[0009][0010]
其中，表示噪声图像张量，表示清晰图像张量，表示足迹，表示随机噪声；
[0011]
s2、采用以下去噪模型，在和同时存在的情况下，从估计清晰图像张量估计清晰图像张量
[0012]
其中，为核范数最小化的清晰地震数据低秩张量模型，表示足迹噪声的utv模型；
[0013]
s3、引入三个辅助变量将中上的无约束最小化问题转化为有约束的最小化问题，然后调用sbi来求解有约束最小化问题，如下所示：
[0014][0015][0016]
s4、将进一步转化为后续的无约束最小化：
[0017][0018]
其中α、β、γ表示bregman惩罚参数，是bregman张量；
[0019]
s5、使用交替最小化方案，求解得到清晰图像张量
[0020]
本发明的有益效果：本发明提出了一种新的utv-lrta三维地震随机加足迹噪声衰减方法，本发明的方法将去噪问题建模为低秩稀疏张量分解模型；具体地：在本发明的utv-lrta模型中，通过有效地将足迹污染的地震图像分解为utv的图像分量和足迹分量来去除足迹噪声，而通过使用张量核范数正则化三阶地震张量的低秩来稀释随机噪声；然后，采用基于sbi的算法对提出的utv-lrta进行了有效求解。通过人工合成实验和现场实验验证了本发明方法的性能；实验结果表明，本发明算法对随机噪声加足迹的去噪能力较好，计算复
杂度较好。
附图说明
[0021]
图1为本实施例提供的地震图像；
[0022]
其在，(a)为带足迹噪声的地震图像原始切片，(b)为原始切片在x轴上的导数，(c)为原始切片在y轴上的导数；
[0023]
图2为本发明的模型示意图；
[0024]
图3为本实施例提供的本发明方法与现有技术在合成数据集上的单道对比示意图；
[0025]
其中，(a)为本发明的utv-lrta算法去噪结果与干净的合成数据在单道上的比较结果，(b)为fr-net算法去噪结果与干净的合成数据在单道上的比较结果，(c)为sgrdl算法去噪结果与干净的合成数据在单道上的比较结果，(d)为drr算法去噪结果与干净的合成数据在单道上的比较结果；
[0026]
图4为本发明的sbi算法在合成数据集上的收敛曲线；
[0027]
图5为本发明的方法与现有技术在合成数据集上的去噪效果对比；
[0028]
其中，(a)为噪声数据，(b)为本发明提出的utv-lrta方法的足迹去除结果，(c)为本发明提出的utv-lrta得到相应的残差，(d)为真实数据，(e)为fr-net方法的足迹去除结果，(f)为fr-net方法得到相应的残差，(g)为sgrdl方法的足迹去除结果，(h)为sgrdl方法得到相应的残差，(i)为drr方法的足迹去除结果，(j)为drr方法得到相应的残差；
[0029]
图6为本发明的方法与现有技术在合成数据集time切片上的去噪效果对比；
[0030]
其中，(a)为原第163个time的切片，(b)为本发明提出的utv-lrta在切片上的足迹去除结果，(c)为本发明提出的utv-lrta方法得到相应的残差，(d)为地面真实数据，(e)为fr-net方法在切片上的足迹去除结果，(f)为fr-net方法得到相应的残差，(g)为sgrdl方法在切片上的足迹去除结果，(h)为sgrdl方法得到相应的残差，(i)为drr方法在切片上的足迹去除结果，(j)为drr方法得到相应的残差；
[0031]
图7为本发明的方法与现有技术在合成数据集inline切片上的去噪效果对比；
[0032]
其中，(a)为原第30个inline的切片，(b)为本发明提出的utv-lrta在切片上的足迹去除结果，(c)为本发明提出的utv-lrta方法得到相应的残差，(d)为地面真实数据，(e)为fr-net方法在切片上的足迹去除结果，(f)为fr-net方法得到相应的残差，(g)为sgrdl方法在切片上的足迹去除结果，(h)为sgrdl方法得到相应的残差，(i)为drr方法在切片上的足迹去除结果，(j)为drr方法得到相应的残差；
[0033]
图8为本发明的方法与现有技术在合成数据集xline切片上的去噪效果对比；
[0034]
其中，(a)为原第70个xline的切片，(b)为本发明提出的utv-lrta在切片上的足迹去除结果，(c)为本发明提出的utv-lrta方法得到相应的残差，(d)为地面真实数据，(e)为fr-net方法在切片上的足迹去除结果，(f)为fr-net方法得到相应的残差，(g)为sgrdl方法在切片上的足迹去除结果，(h)为sgrdl方法得到相应的残差，(i)为drr方法在切片上的足迹去除结果，(j)为drr方法得到相应的残差；
[0035]
图9为本发明的方法与现有技术在penobscot-3d数据集上的去噪效果对比；
[0036]
其中，(a)为噪声，(b)为本发明提出的utv-lrta的足迹去除结果，(c)为本发明提
出的utv-lrta方法得到相应的残差，(d)为fr-net方法的足迹去除结果，(e)为fr-net方法得到相应的残差，(f)为sgrdl方法的足迹去除结果，(g)为sgrdl方法得到相应的残差，(h)为drr方法的足迹去除结果，(i)为drr方法得到相应的残差；
[0037]
图10为本发明的方法与现有技术在penobscot-3d数据集常数时间片上的去噪效果对比；
[0038]
其中，(a)为原第163个time的切片，(b)为本发明提出的utv-lrta在切片上的足迹去除结果，(c)为本发明提出的utv-lrta方法得到相应的残差，(d)为fr-net方法在切片上的足迹去除结果，(e)为fr-net方法得到相应的残差，(f)为sgrdl方法在切片上的足迹去除结果，(g)为sgrdl方法得到相应的残差，(h)为drr方法在切片上的足迹去除结果，(i)为drr方法得到相应的残差；
[0039]
图11为本发明的方法与现有技术在penobscot-3d数据集inline切片上的去噪效果对比；
[0040]
其中，(a)为原第10个inline的切片，(b)为本发明提出的utv-lrta在切片上的足迹去除结果，(c)为本发明提出的utv-lrta方法得到相应的残差，(d)为fr-net方法在切片上的足迹去除结果，(e)为fr-net方法得到相应的残差，(f)为sgrdl方法在切片上的足迹去除结果，(g)为sgrdl方法得到相应的残差，(h)为drr方法在切片上的足迹去除结果，(i)为drr方法得到相应的残差；
[0041]
图12为本发明的方法与现有技术在penobscot-3d数据集xline切片上的去噪效果对比；
[0042]
其中，(a)为原第70个xline的切片，(b)为本发明提出的utv-lrta在切片上的足迹去除结果，(c)为本发明提出的utv-lrta方法得到相应的残差，(d)为fr-net方法在切片上的足迹去除结果，(e)为fr-net方法得到相应的残差，(f)为sgrdl方法在切片上的足迹去除结果，(g)为sgrdl方法得到相应的残差，(h)为drr方法在切片上的足迹去除结果，(i)为drr方法得到相应的残差。
具体实施方式
[0043]
为便于本领域技术人员理解本发明的内容，首先对以下技术术语进行说明：
[0044]
定义1.给定张量和是一个大小为n1×
n4×
n3的张量，定义为：
[0045][0046]
其中，代表全体实数，n1、n2、n4、n3为整数，本实施例中采用n1×
n2×
n3、n2×
n4×
n3、n1×
n4×
n3来表示张量在三个维度上的尺寸大小；
[0047]
一般情况下，fold()和unfold()操作定义如下：
[0048][0049]
x
(1)
表示x在第三个维度n3上的第一个二维切片；
[0050]
定义circ(
·
)为一个n1n3×
n2n3的矩阵：
[0051][0052]
定义2(t-svd).对于其张量svd分解为
[0053][0054]
其中，*表示卷积运算；和是正交张量；是对角张量，θ正面切片(即n1 x n2方向上的切片)是对角矩阵。t-svd和矩阵svd具有相同的特征，即坚持和的正交性质，以及θ的对角线结构。在傅里叶域中，t-svd的每个正面切片等价于一个矩阵svd。t-svd的这个定义引出了张量核范数的概念，在下面进一步定义。
[0055]
定义3.张量核范数(tnn).一个张量的张量核范数，记为定义为的所有正面切片的核范数的平均值，并进一步转化为傅里叶域块对角矩阵的核范数，即
[0056][0057]
其中是在傅里叶域中的块对角矩阵：
[0058][0059]
本发明中相关参数符号的含义说明如表1所示：
[0060]
表1参数符号含义说明
[0061]
[0062][0063]
以下结合附图1-12对本发明的内容进行详细说明：
[0064]
本发明方法的处理过程如图2所示，图2中split bregman iteration optimization表示分裂bregman迭代优化算法，noisy image表示噪声图像，footprint表示足迹噪声，clear image表示干净图像，random noise表示随机噪声，denoising image表示去噪后的图像，unidirectional total variation表示单向总变差，tensor nuclear norm表示张量核范数；本发明的方法具体包括以下步骤：
[0065]
1、首先分析现有技术存在的问题
[0066]
足迹和随机噪声效应一般假设为可加性效应，而不是可乘性效应，混合噪声退化模型可以写成
[0067][0068]
其中，和分别表示噪声图像张量和真实图像张量。形式上，本发明的任务是在足迹和随机噪声同时存在的情况下，从给定的图像张量估计潜在的清晰图像张量具体地说，通过最小化以下方程：
[0069][0070]
由式(13)可知，式(7)中的任务被表述为病态反问题。关键是构造适当的正则化项来约束解，从而去除足迹和随机噪声。基于正则化理论，去噪模型可以表示为数据保真项和正则化项的结合：
[0071][0072]
式(9)中的第一项是数据保真度项，它提供了期望的清晰图像张量和噪声图像张量之间的相似性。第二项是对图像施加约束的正则化项，目的是惩罚图像中不可取的性质。最后一项也是一个正则化项，与项不同，项表示足迹噪声的一些先验知识。τ》0和λ》0都是正则化参数，它控制着数据保真度和两个不同正则化项之间的权衡。
[0073]
2、问题建模
[0074]
本发明针对上述问题的建模过程，包括：干净地震数据和足迹的低秩张量模型和
utv模型；提出了拟议的utv-lrta；提出了能够同时去除高斯噪声和足迹的utv-lrta模型。具体为：
[0075]
1)清晰图像张量的低管秩张量模型：地震去噪问题将获取地面真实数据视为一个低秩矩阵或张量，而噪声会增加该数据矩阵或张量的秩。因此，当数据受到高斯噪声污染时，降阶方法是一种有效的随机噪声衰减方法。在这项工作中，可以考虑基于张量核范数(tnn)最小化的清晰地震数据低秩张量模型，其定义为：
[0076][0077]
众所周知，基于张量核范数最小化的降阶方法在地震图像去噪应用中取得了一定的降噪效果。
[0078]
2)足迹的utv模型：可以通过utv模型充分挖掘足迹的方向性特征。在图1中展示了张量的前片在两个方向上的导数。原始切片如图1(a)所示，切片在x轴和y轴上的导数分别如图1(b)-1(c)所示。由于足迹的方向性，如图1(b)所示，可以观察到条纹严重破坏空间水平梯度，而如图1(c)所示的空间垂直梯度完全不受影响。因此，本发明提出了utv，利用足迹的结构和方向特性，具有更好的细节保持性能，定义为：
[0079][0080]
其中和表示对应方向上的一阶导数算子。
[0081]
3)提出的utv-lrta模型：用低管阶张量模型和(9)式中的utv正则化方法将和分别代入公式(11)，得到关于矩阵x的最终能量泛函：
[0082][0083]
其中，λ1，λ2，和τ表示正则化参数，用于调整每一项的权重，以平衡模型。
[0084]
3、求解公式(12)的方法
[0085]
a.目标函数
[0086]
在使用欠定方程(12)求解时所面临的挑战是，1-范数项本身是不可微的，即在区间的某些点上不可达的函数。分裂bregman迭代(sbi)是解决“基于1范数的正则化最小化问题(当然还有tv模型)的有效工具，由于其快速求解的特性，最近引起了相当大的关注。因此，sbi自然被扩展到式(12)中提出的utv-lrta模型，以求解明确的三阶数据张量通过引入三个辅助变量将式(12)中上的无约束最小化问题转化为有约束的最小化问题，然后调用sbi来求解有约束最小化问题，如下所示：
[0087][0088]
随后，通过严格应用bregman迭代，上述约束问题(13)可以进一步转化为后续的无约束最小化：
[0089][0090]
其中α，β和γ表示bregman惩罚参数，而是bregman张量。
[0091]
b.求解算法
[0092]
使用交替最小化方案来解决式(14)中的问题，即迭代优化一个变量，同时固定其他变量。因此，问题式(14)可以分解为4个更简单的最小化子问题，每个子问题都可以独立求解。
[0093]
1)更新x
l+1
：与相关的子问题为
[0094][0095]
为了求解式(15)，可以解出的每一个正面切片
[0096][0097]
请注意，问题式(16)是一个实变量的凸函数，因此，它更直接地找到最优解，即
[0098][0099]
式中对张量的二维输入额片进行正则矩阵转置，l为迭代数。利用二维傅里叶变换找到上述问题的闭形式解，其解如下：
[0100][0101]
其中
[0102]
[0103][0104]
式中，“*”表示复共轭，“o”表示按分量相乘，和分别表示快速傅里叶变换矩阵及其逆矩阵。
[0105]
2)更新：与相关的子问题为
[0106][0107]
显然，(21)中的子问题必须使用一个张量svd来解决，它类似于矩阵的svd：
[0108][0109]
然后利用软收缩算子得到的闭合形式解如下：
[0110][0111]
其中
[0112][0113]
显然，这种收缩操作的实现是非常直接的。在每个迭代步骤中，只需要计算一个tsvd分解并执行基本的张量操作。
[0114]
3)更新与相关的子问题为
[0115][0116]
(25)中的子问题可以用一个收缩算子来求解：
[0117][0118]
4)更新与相关的子问题为
[0119][0120]
与相关子问题一样，子问题(27)的解析解为：
[0121][0122]
表2在合成数据集上定量比较现有技术和本发明方法的足迹去除结果
[0123]
方法合成数据集drr13.74dbfr-net20.26dbsgrdl21.68db本发明的utv-lrta37.57db
[0124]
表3与现有方法的计算时间比较
[0125]
方法合成数据集penobscot-3ddrr1450.11s410.84sfr-net1173.51s358.11ssgrdl161.79s87.12s本发明的utv-lrta40.71s21.26s
[0126]
5)更新bregman张量和最后，和是bregman张量，更新如下：
[0127][0128]
本发明提出的utv-lrta的求解过程的收敛条件当为满足时，或者迭代次数超过预设阈值时，外层循环停止。在本实施例中，设置∈＝10-5
。另外，sbi方法可以保证分裂bregman迭代算法对式(14)的收敛性，如图4所示。图4中rse为相对标准误差，the convergence area为收敛区域，the number of iterations为迭代次数。
[0129]
通过对合成数据和实际现场数据进行一系列对比实验来确定本发明的utv-lrta模型的有效性。首先，实验配置包括性能指标、使用的对比算法和参数设置，分为两部分：一部分采用合成数据，另一部分采用实际现场数据。
[0130]
a.性能指标
[0131]
在合成数据实验中，地面真实数据是可用的，可以与足迹去除结果进行比较。除了视觉比较之外，很自然地使用定量评估指标来比较fr-net和sota方法的去噪性能。本实施例还使用信噪比(snr)作为定量评估指标。
[0132][0133]
式中分别表示干净图像和估计图像。
[0134]
b.对比算法
[0135]
由于本发明提出的utv-lrta同时具有经典的dcae和utv组件，我们将我们提出的方法与经典的电视正则化dcae模型进行比较，以验证使用utv所带来的好处。根据第一部分的地震文献综述，我们还将我们的方法与最近的足迹去除模型和基于svd的滤波进行了比较，后者学习了如何有效地从污染数据中消除噪声。这些sota方法详细如下。
[0136]
1)fr-net:通过svd将噪声信号的hankel矩阵向量空间分解为信号子空间和噪声子空间；
[0137]
2)统计引导剩余字典学习(statistics-guided residual dictionary learning,sgrdl)：
[0138]
sgrdl生成一个学习信号字典原子的算法框架和从学到的原子分离足迹特征。在这个实例中，sgrdl利用一种基于统计分析的方法将字典原子分类为受足迹影响的或不受足迹影响的。
[0139]
3)阻尼秩减少(damped rank reduction,drr):通过在传统的截断svd模型中引入阻尼算子，drr倾向于将噪声信号的hankel矩阵向量空间划分为有用的信号子空间和噪声子空间。
[0140]
实验中，drr和sgrdl在matlab中实现，cbdnet和fr-net在谷歌colab平台上的tensorflow3；请记住，colab在平台上随机分配gpu给用户。
[0141]
c.参数设置
[0142]
选择合适的参数是许多算法的共同难点，经验调整是确定参数的常用方法。本实施例方法包含α,β,γ,τ,λ1,λ26个正则化参数，参数取决于具体的条纹噪声水平增强的随机噪声电平。虽然本发明的方法涉及很多参数，但参数具有较好的鲁棒性。由于本实施例涉及许多条纹级别，根据实验中测试图像的不同退化级别，这里本实施例设置了α＝0.1，β＝50，γ＝50，τ＝1，λ1＝0.1，λ2＝0.1的人工合成数据和α＝20，β＝10，γ＝10，τ＝1，λ1＝0.1，λ2＝0.1的penobscot工区数据。
[0143]
d.对合成数据的验证
[0144]
首先，将本发明所提出的utv-lrta应用于一个合成物理模型数据集，该数据集包含785
×
790条轨迹，每条轨迹包含400个时间样本。故意在不含地面真实值的数据中加入随机噪声和足迹，从而使带噪声的地震图像成为utv-lrta的输入，如图5(a)、6(a)、7(a)和8(a)所示。注意，utv-lrta根本看不到原始地震图像，因为它只用于计算信噪比值。本发明希望utv lrta能够准确地从相应的有噪声的3d图像中重建出无噪声的原始图像。
[0145]
最后将得到的去噪图像可视化，图5为混合噪声环境下不同方法的三维合成数据，信噪比为6.02db。图5显示了utv-lrta和sota方法对噪声数据量的结果，这些结果由马达加斯加软件包2显示。在信噪比方面，fr-net、drr和sgrdl方法的值分别为21.68、13.74和20.26db，而本发明的utv-lrta模型的值更高，为37.57db，如表2所示。重要的是，我们发现utv-lrta在图像内容方面提供了非常接近地面真实的发现，而sota方法导致了对地面真实的可见偏差，如图5、6、7和8。此外，sota和fr-net方法的足迹残差可以通过单迹比较更直观地证明，如图3所示。正如观察到的，很明显，本发明提出的方法有效地执行足迹去除，而不损害有用的信号。图3中time表示时间，amplitude表示振幅。
[0146]
e.真实现场数据的验证
[0147]
在本实施例中，使用真实的penobscot-3d数据集来测试所提出的方法，该数据集是由dgb地球科学3公开发布的真实转台数据。penobscot-3d数据被浅层时间片上的随机噪声和足迹污染，如图9(a)、10(a)、11(a)和12(a)所示，每个足迹图案由一系列不同长度的条纹组成，这些条纹沿交叉线方向排列。
[0148]
图9、图10、图11、图12给出了所有测试方法的去噪效果。其中，图9对比了本发明提出的utv-lrta和sota方法在整个数据体中的降噪性能，图10、11、12分别为恒定时间切片、xline切片和inline切片的降噪结果。从图9、10、11、12中混合了随机噪声和足迹的噪声可以看出，本发明提出的utv-lrta方法实现了竞争性去噪结果。fr-net在去除足迹方面具有良好的性能，但对随机噪声的抑制效果较差。drr无法抑制足迹噪声，如图10(h)、(i)所示。对于sgrdl，随机和占用空间都被删除；然而，很明显，在图10(f)、(g)所示的噪声降低结果中仍然存在相当大的足迹。在图9、10、11、12中可以观察到，utv-lrta恢复图像的视觉质量有所提高，成功地抑制了随机噪声和足迹，保留了地震图像的局部细节和结构信息。图5-12
中time表示时间，crossline表示联络线地震解释剖面，inline表示主测线地震解释剖面。
[0149]
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。