一种基于变分模态分解的低秩张量地震数据去噪方法

本发明涉及图像处理技术领域，具体涉及一种基于变分模态分解的低秩张量地震数据去噪方法。

背景技术：

地震数据收集会受到多种因素的影响，例如：振动干扰；收集设备故障；以及传输过程中的数据丢失。这些影响以噪声的形式反映在地震数据中，降低了地震数据的质量，并增加了后期地震数据处理和分析的难度。为了减少噪声对地震数据的干扰，学者们提出了许多方法。例如，初始随机噪声降低，通过fx经验模式分解的随机噪声衰减，广义总变异(tv)模型，基于规则邻域块的局部自适应小波阈值降噪等。勘探地球物理学文献中通常存在三种不同类别的噪声衰减方法。第一种是基于空间域的特征。非局部均值算法(nlm)使用图像在空间域中的自相似性来抑制高斯噪声。自适应动态加权中值滤波算法(adwmf)使用基于高斯表面的加权算子对邻域中的像素进行加权，以提高地震数据的质量。第二种基于频率域特性。小波阈值降噪对图像信号进行小波分解以获得小波系数，然后对小波系数进行阈值校正，以达到降噪的目的。curvelet变换去噪算法通过使用频域中信号特征的稀疏表示来消除随机噪声的目的。第三种是基于低秩近似策略，例如高阶奇异值分解(hosvd)用于自然图像去噪。wtr1方法将二维地震资料转化为三维张量，并使用固有的低秩张量近似模型来得到最终的去噪图像。

技术实现要素：

针对现有技术中的上述不足，本发明提供了一种基于变分模态分解的低秩张量地震信号去噪方法。

为了达到上述发明目的，本发明采用的技术方案为：

一种基于变分模态分解的低秩张量地震数据去噪方法，包括如下步骤：

s1、获取地震数据信号，利用变分模态分解将地震数据分解为多个子模态，并求解地震数据信号的三维张量；

s2、构造低秩张量最小化去噪模型，对步骤s1求解的地震张量进行去噪，得到去噪后的地震张量；

s3、将去噪后的地震张量重构为地震剖面，达到去噪后的地震数据。

上述方案的有益效果是，利用vmd将地震资料在频率域中进行分解，通过构造地震资料张量来突出地震资料的频率域信息，并采用低秩张量最小化去噪模型，综合地震资料空间域、频率域信息对其进行去噪。本发明不仅考虑了地震资料的空间域相似性，还考虑了地震资料的频率域相关性。实现了将地震资料由二维转化到三维进行处理，采用低秩张量最小化去噪模型，能够更好的实现对地震资料的去噪，去除噪声的同时又能够更好的保留原地震资料的结构纹理特征。

进一步的，所述步骤s1具体为：

s11、将地震数据信号分解为多个有限带宽模态，其中每个模态的带宽之和最小且各模态之和等于地震数据信号；

s12、利用交替方向乘子法求解地震数据的三维张量。

上述进一步方案的有益效果是：将二维地震数据进一步转化为三维张量，实现了维度的提升。

进一步的，所述步骤s11中将地震数据信号分解为多个有限带宽模态的公式表示为：

其中，t为时间，uk表示第k个有限带宽模态，{uk}为所有有限带宽模态的集合，{ωk}为模态{uk}的中心频率，λ(t)为拉格朗日乘子，α为二次惩罚因子，为求导函数，k为模态索引，δ(t)为狄拉克函数，*为卷积运算符，f(t)表示原始地震信号，j为复数单位，表示为二范数的平方。

上述进一步方案的有益效果是：给出了vmd计算的相关公式，进一步将地震数据信号分解为多个有限带宽模态。

进一步的，所述步骤s12中利用交替方向乘子法求解地震张量的具体步骤为：

s121、初始化{u1},{ω1},λ(t),n＝0；

s122、循环n＝n+1，k＝1:k,利用要求3所述公式更新uk,ωk,λ(t)，其中ω＞0；

s123、重复s122步骤，满足迭代分解精度：时，输出分量uk(k＝1,2,…,k)，ε为分解精度阈值。

上述进一步方案的有益效果是：详细阐述利用交替方向乘子法求解地震张量的具体步骤。

进一步的，所述步骤s2具体包括：

s21、将步骤s1得到的三维张量分解成多个三维全层子块；

s22、任意选取一个三维全层子块，我们以欧几里得距离作为相似性度量，计算步骤s21中每个三维全层子块与它的空间域相似性；

s23、选择空间域相似性小于相似度阈值的三维全层子块进行频率域堆叠，得到新的地震数据的三维张量ti；

s24、利用低秩张量对步骤s23生成的新的地震数据的三维张量ti进行去噪，得到三维全层子块去噪后的张量χi；

s25、固定不同的三维全层子块，重复s23以及s24，可以得到不同的χi，将得到的所有噪后的张量χi重构还原得到去噪后的地震数据的三维张量。

上述进一步方案的有益效果是：对三维张量进行去噪，得到较为干净的地震数据三维张量。

进一步的，所述步骤s23具体为：

s231、将所选择的三维全层子块的每一层按照字典顺序排列为列向量：

s232、对所选择的三维全层子快的所有层按照层数的顺序依次排列成为二维矩阵；

s233、对所有满足空间域相似性小于相似度阈值的三维全层子块执行步骤s231-s232，得到二维矩阵组；

s234、按照三维全层子块的排列顺序将二维矩阵组堆叠，得到新的地震数据的三维张量。

上述进一步方案的有益效果是：得到由三维全层子块组成的新的地震数据的三维张量。

进一步的，所述步骤s24中去噪方式表示为：

其中，χi为去噪后的张量，h(t)为低秩张量，γ是折中参数，t为引入的中间变量，为frobenius范数的平方。

上述进一步方案的有益效果是：利用上述公式，对s234得到的ti进行去噪。

进一步的，所述低秩张量h(t)表示为：

其中，st是张量t经过tucker分解的核心张量；我们将张量t在第i维上进行分解，记为t(i)，则rank(t(i))表示为t(i)的秩。

上述进一步方案的有益效果是：给出低秩张量h(t)的表达式，易于实验的进行。

进一步的，所述步骤s3具体为：对三维张量的每一层进行堆叠求和运算得到去噪后的地震数据信号。

上述进一步方案的有益效果是：将三维张量通过叠加求和转化到二维地震数据中来。

附图说明

图1为本发明基于变分模态分解的低秩张量地震数据去噪方法流程示意图。

图2为本发明实施例中所用合成资料和实际资料展示图。

图3为本发明实施例中合成资料的噪声图和去噪对比图。

图4为本发明实施例中实际资料的去噪对比图。

图5为本发明实施例中实际资料的残差图像。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

s1、获取地震数据信号，利用变分模态分解将地震数据分解为多个子模态，并求解地震数据信号的三维张量；

在本实施例里，变分模态分解(vmd)通过设置合理的预设模态数k，将原始信号分解为一系列具有中心频率ωk(k＝1,2,…,k)的有限带宽模态uk(k＝1,2,…,k)，δ(t)为狄拉克函数，*为卷积运算符。每个uk的估计带宽之和最小，同时各分量信号之和等于原始信号。

在将原始信号分解为uk(k＝1,2,…,k)时，我们根据的是以下公式：

为了将上述约束变分问题转化为无约束变分问题，我们引入拉格朗日乘法算子λ(t)和二次惩罚因子α，将上述公式转化为：

通过交替方向乘子法，我们对其进行求解。

具体步骤为：①初始化{u1},{ω1},λ(t),n＝0；

②循环n＝n+1，k＝1:k,利用步骤s11的公式更新uk,ωk,λ(t)，其中ω＞0；

③重复中间环节②，满足迭代分解精度：时，输出分量uk(k＝1,2,…,k)

s2、构造低秩张量最小化去噪模型，对步骤s1求解的地震张量进行去噪，得到去噪后的地震张量；

本实施例中，具体包括如下步骤：

s21、将步骤s1得到的三维张量分解成多个三维全层子块，经变分模态分解过后，地震资料分解为三维张量其中q是迭代次数，设初始q＝0。dh,dw,ds表示地震张量的空间高度，空间宽度和层数。每层包含有限频带内的地震数据。

s22、计算步骤s21中每个三维全层子块的空间域相似性，本实施例中，为了提升空间域相似性，我们对p^(q)进行块匹配。将p^(q)分成一系列小块pi，(i代表小块数量)。固定一小块pi，本实施例以欧几里得距离作为相似性度量来确定相似性，并找到一组相似的三维全层小块pi(k)，(k＝1,…,n；n是与pi相似的全层小块的数量)。每个pi(k)也是具有两个空间模式和一个光谱模式的三维张量。

s23、选择空间域相似性小于相似度阈值的三维全层子块进行频率域堆叠，得到新的地震数据的三维张量；

本实施例中，根据不同频段的地震资料具有不同的特征，直接将每个频率层的小块进行堆叠是不方便的。因此，本步骤包括如下：

s231、将所选择的三维全层子块的每一层按照字典顺序排列为列向量；对于pi(k)的第k个小块，我们将其第t层记为并将按字典顺序排列为列向量

s232、对所选择的三维全层子块的所有层按照层数的顺序依次排列成为二维矩阵；对于每个pi(k)得到的所有本实施例根据层数t的顺序，将其依次排列成二维矩阵lk，并记为2dflp:lk(k＝1,…,n)。

lk＝unfolding(pi(k))

对所有满足空间域相似性小于相似度阈值的三维全层子块执行步骤s231-s232，得到二维矩阵组，表示为

s234、按照三维全层子块的排列顺序将二维矩阵组堆叠，得到新的地震数据的三维张量，表示为：

根据所有三维全层小块pi(k)(k＝1,…,n)得到的二维矩阵组ωp(i,k)，本实施例按照顺序直接将其堆叠形成三维张量因此，地震张量ti的频率域相关性和空间域相似性可以很好地反映出来。

s24、利用低秩张量对步骤s23生成的新的地震数据的三维张量进行去噪，得到三维全层子块去噪后的张量，

本实施例中，通过构造以下优化问题：

本实施例对地震张量ti进行去噪，得到去噪后的张量χi。公式中，后者的frobenius范数为保真项，γ是折衷参数，t为引入的中间变量。h(t)为引入的张量秩：

式中，是张量t经过tucker分解的核心张量，0＜t＜1。在公式中，它的第一项物理解释为tucker分解中核心张量的大小。第二项约束了目标张量的kronecker基的数量，这符合cp分解的内在机制。因此，tucker分解和cp分解的稀疏性都考虑到了。我们将张量t在第i维上进行分解，记为t(i)。

选取不同的小块重复上述步骤，得到所有去噪后的χi，按照展开方式的逆过程，将所有得到的χi还原回三维张量，即得到去噪后的三维张量

s3、将去噪后的地震张量重构为地震剖面，达到去噪后的地震数据。

即，对三维张量的每一层进行堆叠求和运算，则可以得到去噪后的地震资料。叠加求和公式表示为：p为三维张量的层数；

本实施例中，为了验证本发明提出的地震资料去噪方法的有效性，本发明将所提出的基于变分模态分解(vmd)的低秩张量地震资料去噪方法与nlm方法、ksvd方法、tv方法、ncsr方法进行去噪对比实验。本发明使用的数据集包括模拟和实际地震资料，通过理想原图和去噪后图像计算snr、ssim以及psnr来衡量算法的有效性。

其中，m、n分别表示图像长和宽上的像素点的个数，g(i,j)和f(i,j)分别表示原图像与去噪后的图像在点(i,j)处的灰度值。

其中，μx为x的均值，μy为y的均值，为x的方差，为y的方差，σxy为x和y的协方差，c1＝(k1l)²，c2＝(k2l)²为两个常数，用来避免分母为0，l为像素值的范围，实验取l＝255，k1＝0.01，k2＝0.03。

其中，maxi是表示图像点灰度值的最大数值，如果每个采样点用8位表示，那么maxi就为255，m,n表示二维图像的大小。

本实施例中，为了展示去噪效果优劣，选取信噪比为-1的高斯噪声剖面图，如图3(a)所示。分别使用nlm方法、ksvd方法、tv方法、ncsr方法和本文方法进行去噪处理，结果如图3(b)—图3(f)所示。由图可以看出，加入高斯噪声后模型的断层边界模糊不清，断点处的反射基本淹没在噪声中。对比方法在一定程度上能去除部分噪声，但去噪后断层边界以及断点的恢复效果不明显。从图3(f)所示的去噪剖面可以看出，相比较其他方法，本方法在很好地压制噪声同时，保留了有效信息，同相轴更加光滑，去噪后的断层边界最为清晰，断点位置也最为明确。

如表1所示，用本文方法与对比方法进行去噪处理，测试结果如表1所示。通过表1所列结果可以看出，在不同噪声水平下，本文方法在整体去噪性能(psnr)基本优于其他几种对比方法，同时本方法在图像质量的提升上(ssim)也明显优于其他几种对比方法，信噪比(snr)相较于其他几种对比方法有显著提高，说明经本文方法处理后的图像质量更好。

表1

图4为应用nlm、ksvd、tv、ncsr以及本文方法来处理图4(a)中的实际地震剖面的去噪效果对比。从图4(a)可以看出，整个信号都受到强噪声的干扰，很难识别几乎被完全掩盖的有效信息。仔细对比可发现，对比实验(图4b-e)噪声没有得到有效地压制，恢复边界能力较差，断层边界模糊；本文方法在地震剖面处理结果中，噪声得到了较有效的压制，同相轴更为清晰、连续；断层边界解释更为准确。

图5a-e为五种方法去噪的残差图像。在ncsr处理的方法噪声中，能看到较为明显的断层轮廓，表明经过ncsr算法去噪后，地震剖面的断层信息丢失比较严重。nlm，ksvd，tv和本文方法可以很好地保护地震剖面的基本结构和细节信息，并且本文方法在保留断层结构的情况下消除了更多的随机噪声。综合去噪后剖面和残差剖面可以看出，本文方法去噪效果最好，更有利于后续相关实验。

综上所述，通过模拟和真实数据集实验表明，本发明所提出的一种基于变分模态分解(vmd)的低秩张量地震资料去噪方法，不仅能对地震资料的噪声进行有效去除，同时可以更好地保留地震剖面断层边界信息。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。