CN107561576B - 基于字典学习正则化稀疏表示的地震信号恢复方法 - Google Patents

Abstract

本发明公开了一种基于字典学习正则化稀疏表示的地震信号恢复方法，包括以下步骤：S1、将张量乘积方法应用于张量字典学习过程中，构建目标函数；S2、通过交替迭代算法求解张量稀疏系数；S3、通过拉格朗日对偶的方法训练张量字典；S4、迭代更新张量字典和张量稀疏系数，实现缺失地震信号的重构。本发明提出的基于字典学习正则化稀疏表示的地震信号恢复方法，首先在张量的分解中引入t‑product算子，在此基础上构造新的目标函数，然后分别通过ADM算法和拉格朗日对偶算法求解稀疏系数和张量字典，从而实现地震信号的稀疏编码，最终实现张量数据的稀疏表示和缺失地震信号的恢复，提高了地震数据重构的效果。

Description

基于字典学习正则化稀疏表示的地震信号恢复方法

技术领域

本发明属于地震信号处理技术领域，特别涉及一种基于字典学习正则化稀疏表示的地震信号恢复方法。

背景技术

随着煤层气、页岩气等非常规能源勘探开发的不断深入，地震勘探对数据的规则性和完整性提出了更高要求。然而在野外地震采集过程中，高山的阻挡、河流的阻碍或是湖泊的分布使得检波器的放置变得十分困难；在城市和村庄附近，建筑物的存在也让我们很难在相应的位置放置接收装置；同时由于发射接收设备的损耗使一些地方采集不到地下的信息，造成部分数据的缺失。叠前地震信号的部分缺失，使得形成的叠后数据层位中出现断开、出现不连续的现象。混入噪声的叠前数据同相轴不清晰，波形走势不明显，形成的叠后数据出现模糊，所以对缺失数据重构同时压制混入的噪声变得十分必要。

针对地震数据的缺失或是混入噪声的情况，通过地震数据重构可以填补出相应的缺失部分，进一步在插值的同时有效的压制混入的各种噪声。缺失地震数据插值和去噪是地震数据预处理的重要方面，为后续的反演问题提供了完整的数据支持，可以获取更加完善的地下地质结构信息，为复杂地区勘探开发提供合理有效的指导和强有力的技术支持。

多维地震数据缺失的严重程度在很多情况下会有所不同，即有的地方只有一道数据缺失，而有的地方是多道地震数据缺失，有的地方随机缺失，所以需要采取不同的方法有效且合理的重构多维地震信号。根据地震数据的特点，对于低秩性相对较好的情况，用低秩张量分解的方法重构缺失的且混入噪声的地震数据；对于低秩性不好的情况，用冗余字典学习的方法处理少量缺失的地震数据。

压缩感知和稀疏表示已经成为信号处理中一个比较活跃的话题，该方法有效实现数据的压缩，数据的传输和数据的恢复。基于稀疏表示的字典学习是地震信号重构的另一种有效方法，该方法充分利用地震数据自身的特点学习其相应的特征，能够更好的表示地震数据的有效信息。最初始的稀疏表示方法采用固定基的线性组合，但是这样并不能达到最优的稀疏性，为此构造自适应的冗余字典成为一个发展的方向。冗余字典不再强制不同原子间的正交性，因为正交性限制其自适应性的范围。冗余字典中的每一个列向量仍然可以看成一个基，不过此时的基是不固定的，可以随着数据的改变而改变，新输入的数据可以通过这些自适应基的线性组合来表示，每个基的系数形成了稀疏系数分布于稀疏矩阵中。稀疏表示过程分为冗余字典的训练和稀疏系数的求解，通过依次迭代处理来得到理想的结果。

在稀疏字典学习问题中，仅仅知道输入数据来获取数据对/> 使得目标函数/>最小；并且满足/>具有稀疏性(/>中仅有少量的非0元素)，其中A可以认为是基矩阵，如果没有噪声，可以认为/>传统的字典学习方法无论在图像处理还是在地震数据去噪、缺失补全中，大部分的输入数据都是二维的，而对于多维数据也是先将其转化为低维数据进行处理，这样就忽略了多维数据间固有的相关性和空间上的相似性。更进一步来说，二维数据的处理效果会在一定程度上劣于直接选取多维数据进行运算。因此，常规的冗余字典学习方法往往过于粗糙，得到的地震数据不够准确，难以满足实际生产需要。

发明内容

本发明的目的在于克服现有技术的不足，提供一种首先在张量的分解中引入t-product算子，在此基础上构造新的目标函数，然后分别通过ADM算法和拉格朗日对偶算法求解稀疏系数和张量字典，从而实现地震信号的稀疏编码，最终实现张量数据的稀疏表示和缺失地震信号的恢复基于字典学习正则化稀疏表示的地震信号恢复方法。

本发明的目的是通过以下技术方案来实现的：基于字典学习正则化稀疏表示的地震信号恢复方法，包括以下步骤：

S1、将张量乘积方法应用于张量字典学习过程中，构建目标函数；

S2、通过交替迭代算法求解张量稀疏系数；

S3、通过拉格朗日对偶的方法训练张量字典；

S4、迭代更新张量字典和张量稀疏系数，实现缺失地震信号的重构。

进一步地，所述步骤S1具体实现方法为：将每一个输入向量y∈Rⁿ利用稀疏向量a₁,a₂,...,a_p∈Rⁿ表示，稀疏向量的系数x∈Rⁿ，所以y≈∑_ja_jx_j；输入向量与其稀疏向量表达式之间的误差y-∑_ja_jx_j符合0均值、协方差为σ的高斯分布；

对于稀疏系数，每个系数x_j的先验分布表示为P(x_j)∝exp(-αφ(x_j))；其中α为自定义的常数；φ(x_j)为稀疏函数：

因此，字典学习稀疏表示模型表示为：

其中，λ为调节参数，用于平衡误差项和稀疏正则项之间的权重；Y、A、X分别表示二维数据、二维字典和二维稀疏系数；||·||_F表示F范数；||·||₁表示张量的l₁范数；A∈C表示二维字典符合DCT或者正交结构；

引入t-product张量乘积，将字典学习拓展到多维数据处理过程中，对于三阶张量数据其目标函数的形式为，

其中，表示学习的张量字典；/>为张量稀疏系数；*表示循环卷积运算符；/>表示张量字典/>符合DCT或者正交结构；通过学习数据对/>从而实现张量的稀疏编码。

进一步地，所述步骤S2具体实现方法为：对于多维地震数据，选取三维地震数据作为张量的输入数据，三阶张量输入数据张量稀疏系数/>张量字典则目标函数为：

其中，*表示循环卷积运算符，式(5)表明张量字典的侧切面满足数据归一性约束；

通过ADM的方法迭代求解稀疏系数的学习，其中需要时域和频域的转换，通过DFT变换将其转化到频域对每个频率的数据进行处理；最终通过迭代处理不断更新张量稀疏系数；具体实现过程为：

S21、初始化调节参数λ、平衡参数ρ，初始化最大迭代次数Miter；

S22、对稀疏系数进行迭代求解：当求解张量稀疏系数时，将张量字典作为已知变量，因此目标函数表示为：

将公式(6)变形为：

添加变量将式(7)转化为增广拉格朗日的形式：

在式(8)中，平衡参数ρ＞0；

采用ADM的方法对式子(8)求解，得到：

针对变量求解：

对于 表示/>的第k个前向切面；式(12)变换到频域在每个频率切片上求解：

函数是可导的，令/>得到：

针对变量求解：

式(15)中的||·||₁范数的求解采用近似操作符P_τ(·)，即用软阈值操作符来逐点处理每个元素；因此，式子(15)的解为其中软阈值操作符为S_τ(·)＝sign(·)max(|·|-τ,0)；如果是l₀范数的情况，最优近似采用硬阈值的方法，即/>其中硬阈值操作符为H_τ(·)＝max(|·|-τ,0)；

针对变量求解：

当满足收敛条件或是达到最大迭代次数时停止迭代，得到稀疏系数。

进一步地，所述步骤S3具体实现方法为：

固定稀疏系数从目标函数式(3)中学习张量字典/>将优化问题转化为：

[r]表示1到r的集合，r表示切面的个数；中每个侧切面的约束使式(17)变为非凸问题，同时因为循环卷积*的存在造成求解上的困难，通过循环卷积*所具有的性质将式(17)转化为一系列子问题，通过求解每个子问题最终得到/>的解的形式，具体包括以下子步骤：

S31、将优化问题变换到频域，得到：

其中，表示将/>变换到频域；/>分别表示/>的频域形式；表示/>的第k个前向切面；式(20)表示k个切面的F范数的和小于或等于k；

基于拉格朗日对偶算法，将式(20)表示为：

其中，对偶变量β_j≥0，j∈[r]；设Σ＝diag(λ)，得到：

式子(22)为原问题，对求导，得到：

令得到：

将带入原问题，求得对偶问题：

对偶函数关于λ_i的梯度表示为：

其中e_i表示在i处为1的单位向量；用牛顿法或是共轭梯度法求解对偶问题，通过最大化将求解的结果Σ带入式子(28)中：

根据式(28)求得的求得张量字典/>

本发明的有益效果是：

1、本发明提出的基于字典学习正则化稀疏表示的地震信号恢复方法，首先在张量的分解中引入t-product算子，在此基础上构造新的目标函数，然后分别通过ADM算法和拉格朗日对偶算法求解稀疏系数和张量字典，从而实现地震信号的稀疏编码，最终实现张量数据的稀疏表示和缺失地震信号的恢复，提高了地震数据重构的效果。

2、本发明根据二维字典学习的方法采用新的张量乘积将其应用于高维张量稀疏表示中，有效的解决了部分从二维数据扩展到多维数据所带来的问题，相比二维稀疏编码的方法，张量稀疏编码更加有效的利用了多维数据的空间信息和数据间的紧密相关性。

附图说明

图1为本发明的地震信号恢复方法流程图；

图2为本实施例原始切面的数据；

图3为本实施例的整体数据随机列缺失15％后该切面的效果图；

图4为本实施例重构后的切面图；

图5为本实施例重构后的残差效果图。

具体实施方式

下面结合附图进一步说明本发明的技术方案。

如图1所示，基于字典学习正则化稀疏表示的地震信号恢复方法，包括以下步骤：

S1、将张量乘积方法应用于张量字典学习过程中，构建目标函数；具体实现方法为：稀疏编码就是通过部分基向量的线性组合来近似趋近输入的向量，基向量的组合可以有效的提取输入数据的主要特征。将每一个输入向量y∈Rⁿ利用稀疏向量a₁,a₂,...,a_p∈Rⁿ表示，稀疏向量的系数x∈Rⁿ，所以y≈∑_ja_jx_j；输入向量与其稀疏向量表达式之间的误差y-∑_ja_jx_j符合0均值、协方差为σ的高斯分布；

对于稀疏系数，每个系数x_j的先验分布表示为P(x_j)∝exp(-αφ(x_j))；其中α为自定义的常数；φ(x_j)为稀疏函数：

因此，字典学习稀疏表示模型表示为：

其中，λ为调节参数，用于平衡误差项和稀疏正则项之间的权重；Y、A、X分别表示二维数据、二维字典和二维稀疏系数；||·||_F表示F范数；||·||₁表示张量的l₁范数；A∈C表示二维字典符合DCT或者正交结构；

引入t-product张量乘积，将字典学习拓展到多维数据处理过程中，对于三阶张量数据其目标函数的形式为，

其中，表示学习的张量字典；/>为张量稀疏系数；*表示循环卷积运算符；/>表示张量字典/>符合DCT或者正交结构；通过学习数据对/>从而实现张量的稀疏编码。

S2、通过交替迭代算法求解张量稀疏系数；具体实现方法为：对于多维地震数据，选取三维地震数据作为张量的输入数据，三阶张量输入数据张量稀疏系数张量字典/>则目标函数为：

其中，*表示循环卷积运算符，式(5)表明张量字典的侧切面满足数据归一性约束；

通过ADM的方法迭代求解稀疏系数的学习，其中需要时域和频域的转换，通过DFT变换将其转化到频域对每个频率的数据进行处理；最终通过迭代处理不断更新张量稀疏系数；具体实现过程为：

S21、初始化调节参数λ、平衡参数ρ，初始化最大迭代次数Miter；

S22、对稀疏系数进行迭代求解：当求解张量稀疏系数时，将张量字典作为已知变量，因此目标函数表示为：

将公式(6)变形为：

添加变量将式(7)转化为增广拉格朗日的形式：

在式(8)中，平衡参数ρ＞0；

采用ADM的方法对式子(8)求解，得到：

针对变量求解：

对于 表示/>的第k个前向切面；式(12)变换到频域在每个频率切片上求解：

函数是可导的，令/>得到：

针对变量求解：

式(15)中的||·||₁范数的求解采用近似操作符P_τ(·)，即用软阈值操作符来逐点处理每个元素；因此，式子(15)的解为其中软阈值操作符为S_τ(·)＝sign(·)max(|·|-τ,0)；如果是l₀范数的情况，最优近似采用硬阈值的方法，即/>其中硬阈值操作符为H_τ(·)＝max(|·|-τ,0)；

针对变量求解：

当满足收敛条件或是达到最大迭代次数时停止迭代，得到稀疏系数。

S3、通过拉格朗日对偶的方法训练张量字典；具体实现方法为：

固定稀疏系数从目标函数式(3)中学习张量字典/>将优化问题转化为：

[r]表示1到r的集合，r表示切面的个数；中每个侧切面的约束使式(17)变为非凸问题，同时因为循环卷积*的存在造成求解上的困难，通过循环卷积*所具有的性质将式(17)转化为一系列子问题，通过求解每个子问题最终得到/>的解的形式，具体包括以下子步骤：

S31、将优化问题变换到频域，得到：

其中，表示将/>变换到频域；/>分别表示/>的频域形式；表示/>的第k个前向切面；式(20)表示k个切面的F范数的和小于或等于k；

基于拉格朗日对偶算法，将式(20)表示为：

其中，对偶变量β_j≥0，j∈[r]；设Σ＝diag(λ)，得到：

式子(22)为原问题，对求导，得到：

令得到：

将带入原问题，求得对偶问题：

对偶函数关于λ_i的梯度表示为：

其中e_i表示在i处为1的单位向量；用牛顿法或是共轭梯度法求解对偶问题，通过最大化将求解的结果Σ带入式子(28)中：

根据式(28)求得的求得张量字典/>

S4、迭代更新张量字典和张量稀疏系数，实现缺失地震信号的重构。

下面通过具体实施例进一步验证本发明的地震信号恢复方法的有效性。本实施例采用F3工区的数据，从中选取了大小为100*100*20的三维数据A0，对基于字典学习正则化稀疏表示的地震信号恢复方法进行验证。对F3工区中的A0数据整体缺失15％，然后进行数据重构。为了方便，只显示xline号为15的切面来查看恢复前后的效果。图2为原始切面的数据，从图中可知原始数据有的地方同相轴比较清晰、连续性强，而有的地方同相轴不清晰、连续性弱，总的来说包含多种信息特征。图3为整体数据随机列缺失15％后该切面的效果，其中有的相邻道全部缺失。图4为利用本发明的方法重构后的切面图，图5为利用本发明的方法重构后的残差效果图。从图中可以看出，利用本发明提出的算法能将整道缺失的数据恢复出来，算法残差较小，缺失部分恢复后波动比较小，由此可知本发明提出的算法能够更好的利用多维数据的信息。

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

Claims (1)

1.基于字典学习正则化稀疏表示的地震信号恢复方法，其特征在于，包括以下步骤：

S1、将张量乘积方法应用于张量字典学习过程中，构建目标函数；具体实现方法为：将每一个输入向量y∈Rⁿ利用稀疏向量a₁,a₂,...,a_p∈Rⁿ表示，稀疏向量的系数x∈Rⁿ，所以y≈∑_ja_jx_j；输入向量与其稀疏向量表达式之间的误差y-∑_ja_jx_j符合0均值、协方差为σ的高斯分布；

对于稀疏系数，每个系数x_j的先验分布表示为P(x_j)∝exp(-αφ(x_j))；其中α为自定义的常数；φ(x_j)为稀疏函数：

因此，字典学习稀疏表示模型表示为：

其中，λ为调节参数，用于平衡误差项和稀疏正则项之间的权重；Y、A、X分别表示二维数据、二维字典和二维稀疏系数；||·||_F表示F范数；||·||₁表示张量的l₁范数；A∈C表示二维字典符合DCT或者正交结构；

引入t-product张量乘积，将字典学习拓展到多维数据处理过程中，对于三阶张量数据其目标函数的形式为，

其中，表示学习的张量字典；/>为张量稀疏系数；*表示循环卷积运算符；/>表示张量字典/>符合DCT或者正交结构；通过学习数据对/>从而实现张量的稀疏编码；

S2、通过交替迭代算法求解张量稀疏系数；具体实现方法为：对于多维地震数据，选取三维地震数据作为张量的输入数据，三阶张量输入数据张量稀疏系数/>张量字典/>则目标函数为：

其中，*表示循环卷积运算符，式(5)表明张量字典的侧切面满足数据归一性约束；

通过ADM的方法迭代求解稀疏系数的学习，其中需要时域和频域的转换，通过DFT变换将其转化到频域对每个频率的数据进行处理；最终通过迭代处理不断更新张量稀疏系数；具体实现过程为：

S21、初始化调节参数λ、平衡参数ρ，初始化最大迭代次数Miter；

S22、对稀疏系数进行迭代求解：当求解张量稀疏系数时，将张量字典作为已知变量，因此目标函数表示为：

将公式(6)变形为：

添加变量将式(7)转化为增广拉格朗日的形式：

在式(8)中，平衡参数ρ＞0；

采用ADM的方法对式子(8)求解，得到：

针对变量求解：

对于 表示/>的第k个前向切面；式(12)变换到频域在每个频率切片上求解：

函数是可导的，令/>得到：

针对变量求解：

式(15)中的||·||₁范数的求解采用近似操作符P_τ(·)，即用软阈值操作符来逐点处理每个元素；因此，式子(15)的解为其中软阈值操作符为S_τ(·)＝sign(·)max(|·|-τ,0)；如果是l₀范数的情况，最优近似采用硬阈值的方法，即/>其中硬阈值操作符为H_τ(·)＝max(|·|-τ,0)；

针对变量求解：

当满足收敛条件或是达到最大迭代次数时停止迭代，得到稀疏系数；

S3、通过拉格朗日对偶的方法训练张量字典；具体实现方法为：

固定稀疏系数从目标函数式(3)中学习张量字典/>将优化问题转化为：

[r]表示1到r的集合，r表示切面的个数；中每个侧切面的约束使式(17)变为非凸问题，同时因为循环卷积*的存在造成求解上的困难，通过循环卷积*所具有的性质将式(17)转化为一系列子问题，通过求解每个子问题最终得到/>的解的形式，具体包括以下子步骤：

S31、将优化问题变换到频域，得到：

其中，表示将/>变换到频域；/>分别表示/>的频域形式；/>表示/>的第k个前向切面；式(20)表示k个切面的F范数的和小于或等于k；

基于拉格朗日对偶算法，将式(20)表示为：

其中，对偶变量β_j≥0，j∈[r]；设Σ＝diag(λ)，得到：

式子(22)为原问题，对求导，得到：

令得到：

将带入原问题，求得对偶问题：

对偶函数关于λ_i的梯度表示为：

其中e_i表示在i处为1的单位向量；用牛顿法或是共轭梯度法求解对偶问题，通过最大化将求解的结果Σ带入式子(28)中：

根据式(28)求得的求得张量字典/>

S4、迭代更新张量字典和张量稀疏系数，实现缺失地震信号的重构。