CN111025385A - 一种基于低秩和稀疏约束的地震数据重建方法 - Google Patents

Abstract

本发明提供一种基于低秩和稀疏约束的地震数据重建方法，包括：S1、对观测地震数据进行预变换得到变换结果矩阵，作为迭代初始值；S2、在每一次更新过程中，对变换结果矩阵进行奇异值分解，利用分解得到的左、右奇异矩阵构建正交矩阵；S3、利用正交矩阵建立模型，并采用迭代算法求解所述模型，得到更新后的变换结果矩阵；S4、判断更新后的变换结果矩阵是否收敛，若收敛则执行步骤S5，否则将更新后的变换结果矩阵用作下一次更新，回到步骤S2；S5、对收敛的变换结果矩阵进行反变换得到重建的地震数据。本发明的有益效果：充分考虑了地震数据的冗余结构，兼顾细节和全局结构特点，显著提高了重建数据的质量；具有较高的峰值信噪比。

Description

一种基于低秩和稀疏约束的地震数据重建方法

技术领域

本发明涉及地震信号处理领域，尤其涉及一种基于低秩和稀疏约束的地震数据重建方法。

背景技术

地震勘探是研究地下地质构造的重要方法，然而受限于地震数据采集设备和地质环境等其他因素，获得的地震数据通常并不完整。缺失的地震数据会直接影响后续的偏移成像、反演以及地质构造的解释，因此，对地震数据的重建尤为重要。

经典的地震数据重建方法有：(1)基于空间预测滤波的重建方法，将待重建的地震数据与滤波器进行卷积，包括反假频的f-x域地震道插值方法和i-x域的预测误差滤波技术；(2)基于波动方程的重建方法，利用波的传播进行地震数据重建，所述方法需要波的速度等先验信息；(3)基于稀疏变换的重建方法，利用稀疏性先验信息，所述方法分为固定基的方法和基于学习的方法，其中，固定基的方法包括Fourier变换、Radon变换和Curvelet变换等，基于学习的方法通常自适应学习使数据表示更为稀疏的字典；(4)基于低秩的地震数据重建方法，利用低秩性先验信息，所述方法基于假设：完整的地震数据具有低秩结构，缺失地震道和随机噪声会增加目标矩阵或张量的秩，由此，地震数据的重建问题可以转化为矩阵或张量的降秩问题。

上述现有的地震数据重建方法中，稀疏性先验和低秩性先验由于能够获得较好的重建效果而成为目前较为流行的先验条件，受到越来越多研究者的关注。稀疏性先验、低秩性先验是从不同角度考虑数据的冗余结构，将地震数据的稀疏性和低秩性联合考虑可以充分挖掘并利用地震数据中的这种冗余结构，兼顾细节和全局结构的特点，更加有助于提高重建的质量。

发明内容

有鉴于此，本发明提供了一种结合低秩性和稀疏性的地震数据重建算法，充分挖掘并利用了地震数据中的冗余结构，兼顾细节和全局结构的特点，显著提高了地震数据重建的质量。

本发明提供一种基于低秩和稀疏约束的地震数据重建方法，包括以下步骤：

S1、对观测地震数据M进行预变换得到初始变换结果矩阵，利用所述初始变换结果矩阵进行迭代更新；

S2、在第l次更新过程中，对上一次更新得到的变换结果矩阵P_l进行奇异值分解，利用分解得到的左、右奇异矩阵构建两个正交矩阵A_l、B_l；

S3、利用所述正交矩阵A_l、B_l建立变换结果矩阵P的更新模型，所述模型为：

目标函数：

约束条件：

采用迭代算法求解所述模型描述的最小值问题，得到更新后的变换结果矩阵P_l+1；其中，l表示更新次数，

表示所述步骤S1中预变换的反运算；λ表示稀疏项的权重，λ＞0，

表示使地震数据稀疏的变换，N＝P，A_l、B_l表示步骤S2中构建的正交矩阵；M表示观测地震数据，Ω表示未缺失数据的索引，

表示线性算子，用于使X的未缺失部分与M一致，||·||_*表示核范数运算，Tr(·)表示求矩阵的迹，||·||₁表示L₁范数运算；

S4、判断更新后的变换结果矩阵是否收敛，若收敛，则执行步骤S5，否则将所述更新后的变换结果矩阵用于下一次更新，回到步骤S2；

S5、对收敛的变换结果矩阵进行反变换得到重建的地震数据。

进一步地，所述步骤S1中的预变换为纹理块变换；所述步骤S3中，

表示反纹理块变换算子。

进一步地，所述步骤S2的具体过程为：对上一次更新得到的变换结果矩阵P_l进行奇异值分解，得到左奇异矩阵U_l和右奇异矩阵V_l：

[U_l，∑_l，V_l]＝svd(P_l)，

其中，

m、n分别表示矩阵P_l的行数和列数，u₁，u₂，…，u_m表示左奇异矩阵U_l的列向量，v₁，v₂，…，v_n表示右奇异矩阵V_l的列向量，

表示对应的奇异值矩阵；利用所述左奇异矩阵U_l和右奇异矩阵V_l的前r列向量，构造两个正交矩阵A_l、B_l，r＜min(m，n)且为正整数，其中：

A_l＝(u₁，u₂，…，u_r)^T B_l＝(v₁，v₂，…，v_r)^T。

进一步地，所述步骤S3中，通过构造增广的拉格朗日函数来迭代求解所述模型描述的最小值问题，具体过程为：引入三个拉格朗日乘子Y、Z、F和惩罚项系数，当迭代系数k＝1时，给定迭代初值X¹、W¹、N¹、Y¹、Z¹、F¹，其中N¹＝P_l，X¹、W¹为维数大小与M相同的随机矩阵，Z¹为维数大小与M相同的零矩阵，Y¹、F¹均为零矩阵，且矩阵维度与P_l相同；在进行第k次迭代时，依次更新P^k+1、X^k+1、W^k+1、N^k+1、Y^k+1、Z^k+1、F^k+1七个参数；迭代过程中，若N^k+1满足收敛条件，则迭代结束，得到P_l+1＝P^k+1，否则递增k，继续进行迭代。

进一步地，第k次迭代的具体过程如下：

使用奇异值收缩方法更新P^k+1：

上述运算式的计算方式为：

其中，x＝U∑V^T表示x的奇异值分解，且S_τ(∑)＝diag(max{σ_i-τ，0})，σ_i表示奇异值矩阵∑中的元素，

表示纹理块变换算子，β表示增广的拉格朗日函数中的惩罚项系数；

更新X^k+1：

其中，Ω^c表示缺失数据的索引，

表示线性算子；

更新W^k+1：

其中，上述运算式的计算方式为：

sgn(·)表示符号函数，λ表示目标函数中稀疏项的权重；

更新N^k+1：

判断N^k+1是否收敛：

||N^k+1-N^k||_F≤ε₁，

其中，||·||_F表示Frobenius范数，ε₁为给定的第一精度阈值；当满足上述条件时，停止迭代，输出P_l+1＝P^k+1；

更新Y^k+1：

Y^k+1＝Y^k+β(N^k+1-P^k+1)；

更新Z^k+1：

更新F^k+1：

进一步地，所述步骤S4中，收敛条件为：

||P_l+1-P_l||_F≤ε₂，

式中，ε₂为给定的第二精度阈值；若满足所述收敛条件，则迭代结束，执行步骤S5；否则递增l，回到步骤S2。

本发明提供的技术方案带来的有益效果是：本发明同时考虑了地震数据具有的低秩属性和稀疏属性，并从不同角度考虑了地震数据的冗余结构，充分挖掘并利用地震数据中的冗余特性，兼顾了细节和全局结构特点，显著提高了地震数据的重建效果。

附图说明

图1是本发明实施例一提供的基于低秩和稀疏约束的地震数据重建方法的流程图；

图2是本发明实施例一重建得到的地震数据与完整地震数据的PSNR对比图；

图3是本发明实施例一重建得到的地震数据与完整地震数据的效果对比图；

图4是本发明实施例二重建得到的地震数据与完整地震数据的效果对比图；

图5是本发明实施例二重建得到的地震数据与完整地震数据的频谱分析图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。

实施例一

请参考图1，本实施例提供了一种基于低秩和稀疏约束的地震数据重建方法，包括以下步骤：

S1、对观测地震数据进行预变换得到变换结果矩阵；本实施例采用仿真地震数据，共256道，每道数据包含256个时间采样点，所述仿真地震数据的缺失比例从10％到70％，所述不同缺失比例的仿真地震数据构成256×256维观测地震数据M；优选地，本实施例的预变换采用纹理块变换(texture patch-based transform)，所述预变换的目的在于使得观测地震数据具有低秩性，也可选择其他变换方式，比如Hankel变换等。

为更准确地重建地震数据，本实施例同时考虑地震数据的低秩性和稀疏性，并结合截断核范数以及L₀范数，构建一个新的地震数据重建模型：

其中，

表示反纹理块变换算子，

表示地震数据进行纹理块变换后得到的变换结果矩阵，m、n分别表示矩阵P的行数和列数，本实施例中M进行纹理块变换后m＝64、n＝1024；||·||_r表示截断核范数运算，

σ_i为P的第i个奇异值，||·||₀表示L₀范数运算；

表示使地震数据稀疏的变换，比如DCT变换、小波变换、曲波变换等，λ表示稀疏项的权重，λ＞0；Ω表示未缺失数据的索引，

表示线性算子，用于使X的未缺失部分与M一致。

需要说明的是，本实施例对纹理块变换结果矩阵P进行迭代更新来求解所述地震数据重建模型(1)，l表示更新次数，步骤S1中由观测地震数据M进行纹理块变换得到的变换结果矩阵作为迭代更新的初始值。模型中截断核范数以及L₀范数均为非凸，难以求解，因此本实施例将其转换成其他易于求解的形式，所述L₀范数采用L₁范数代替，所述截断核范数转化为核范数减去构造的矩阵的迹的形式，具体求解过程包括步骤S2-S5。

S2、在进行第l次更新时，对矩阵P_l进行奇异值分解，得到左奇异矩阵U_l和右奇异矩阵V_l：

[U_l，∑_l，V_l]＝svd(P_l)，

其中，

m、n分别表示矩阵P_l的行数和列数，u₁，u₂，…，u_m表示左奇异矩阵U_l的列向量，v₁，v₂，…，v_n表示右奇异矩阵V_l的列向量；

表示对应的奇异值矩阵。利用所述左奇异矩阵U_l和右奇异矩阵V_l的前r个列向量，构造两个正交矩阵A_l、B_l：

A_l＝(u₁，u₂，…，u_r)^T B_l＝(v₁，v₂，…，v_r)^T，

其中，r＜min(m，n)，r即为地震数据重建模型(1)中截断核范数的r参数，本实施例中r＜64。

S3、利用所述正交矩阵A_l、B_l建立模型对第l次更新中的变换结果矩阵P_l进行更新，所述模型如下：

目标函数：

约束条件：

其中，P_l+1表示更新后的变换结果矩阵，

N＝P_l，λ表示稀疏项的权重，λ＞0；||·||_*表示核范数运算，Tr(·)表示求矩阵的迹，||·||₁表示L₁范数运算；所述模型(2)即为地震数据重建模型(1)的具体求解形式。

本实施例通过构造增广的拉格朗日函数来迭代求解所述模型(2)描述的最小值问题，得到更新后的变换结果矩阵P_l+1。具体地，引入三个拉格朗日乘子Y、Z、F和惩罚项系数，当迭代系数k＝1时，给定迭代初值X¹、W¹、N¹、Y¹、Z¹、F¹，其中N¹＝P_l，X¹、W¹为维数大小与M相同的随机矩阵，Z¹为维数大小与M相同的零矩阵，Y¹、F¹均为零矩阵，且矩阵维度与P_l相同；在进行第k次迭代时，依次更新P^k+1、X^k+1、W^k+1、N^k+1、Y^k+1、Z^k+1、F^k+1七个参数；迭代过程中，若N^k+1满足收敛条件，则迭代结束，得到P_l+1＝P^k+1，否则递增k，即k＝k+1，继续进行更新。所述求解过程如下：

S301、使用奇异值收缩方法更新P^k+1：

上述运算式的计算方式为：

其中，x＝U∑V^T表示x的奇异值分解，且S_τ(∑)＝diag(max{σ_i-τ，0})，σ_i表示奇异值矩阵∑中的元素，β表示增广的拉格朗日函数中的惩罚项系数。

S302、更新X^k+1：

其中，Ω^c表示缺失数据的索引，

表示线性算子。

S303、更新W^k+1：

其中，上述运算式的计算方式为：

sgn(·)表示符号函数，λ表示目标函数中稀疏项的权重。

S304、更新N^k+1：

判断N^k+1是否收敛：

||N^k+1-N^k||_F≤ε₁，

其中，||·||_F表示Frobenius范数，ε₁为给定的第一精度阈值；当满足上述条件时，停止迭代，输出P_l+1＝P^k+1；

S305、更新Y^k+1：

Y^k+1＝Y^k+β(N^k+1-P^k+1)；

S306、更新Z^k+1：

S307、更新F^k+1：

上述迭代过程结束后，得到P_l+1。

S4、判断经步骤S3更新后的变换结果矩阵P_l+1是否收敛；所述收敛条件为：

||P_l+1-P_l||_F≤ε₂，

其中，ε₂为给定的第二精度阈值；若满足上述收敛条件，则迭代结束，执行步骤S5；否则递增l，即l＝l+1，回到步骤S2。

S5、对收敛的变换结果矩阵P_l+1进行反纹理块变换

得到重建的地震数据X。

本实施例采用峰值信噪比(peak signal-to-noise ratio，PSNR)对重建效果进行评价：

其中，peakval表示地震数据的最大数值，MSE表示重建的地震数据和原始未缺失地震数据的均方误差；所述峰值信噪比值越大，表明重建的地震数据与完整的地震数据越接近，重建效果越好。

请参考图2，分别利用本实施例的基于低秩和稀疏约束的重建方法(TNNSR)，以及基于截断核范数的重建方法(TNNR)，对10％-70％缺失比例的仿真地震数据进行重建，得到PSNR的变化趋势图。从图2可以看出，随着缺失比例的增加，两种方法重建的地震数据的PSNR均会下降，但基于低秩和稀疏约束的重建方法明显在各种缺失比例下均要优于TNNR方法。

两种方法重建得到的地震图像也有明显差异，请参考图3，图3(a)表示完整的地震数据，图3(b)表示随机整道缺失70％的地震数据，图3(c)表示本实施例采用基于低秩和稀疏约束的方法重建得到的地震数据，图3(d)表示采用TNNR方法重建得到的地震数据。对比图3(c)、图3(d)中圆圈部分数据，TNNR方法重建得到的数据明显没有采用本实施例重建得到的数据光滑，本实施例重建得到的图像更加接近图3(a)中的完整地震数据。

实施例二

本实施例采用真实地震数据，共232道，每道数据包含512个时间采样点，缺失比例为50％的观测地震数据M，即M为512×232维矩阵，进行纹理块变换后的变换结果矩阵

m＝64，n＝1856，其余过程与实施例一相同；同样采用PSNR对重建效果进行评价，本实施例重建得到的地震数据的PSNR为33.1841dB，采用TNNR方法进行重建的地震数据的PSNR为30.1977dB；请参考图4，其为本实施例重建得到的地震数据与完整地震数据的效果对比图，图4(a)表示表示完整的地震数据，图4(b)表示随机整道缺失50％的地震数据，图4(c)表示本实施例重建的地震数据，图4(d)表示采用TNNR方法进行重建的地震数据，其中，图4(d)中仍存在未重建完整的“竖道”，且图4(c)、(d)中的方框部分数据也证明了本实施例重建得到的地震数据优于TNNR方法重建得到的地震数据。

请参考图5，其为本实施例重建得到的地震数据与完整地震数据的频谱分析图，图5(a)表示完整地震数据的频谱图，图5(b)表示缺失比例为50％的地震数据的频谱图，图5(c)表示本实施例重建得到的地震数据的频谱图，图5(d)表示采用TNNR方法重建得到的地震数据的频谱图；图5(c)相较于图5(d)，聚焦性更好、能量更集中，也更接近原始图像，即图5(a)，说明对于真实地震数据，本发明提供的基于低秩和稀疏约束的地震数据重建方法优于传统的TNNR方法。

在本文中，所涉及的前、后、上、下等方位词是以附图中零部件位于图中以及零部件相互之间的位置来定义的，只是为了表达技术方案的清楚及方便。应当理解，所述方位词的使用不应限制本申请请求保护的范围。

在不冲突的情况下，本文中上述实施例及实施例中的特征可以相互结合。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于低秩和稀疏约束的地震数据重建方法，其特征在于，包括以下步骤：

S1、对观测地震数据M进行预变换得到初始变换结果矩阵；

S2、在第l次更新过程中，对上一次更新得到的变换结果矩阵P_l进行奇异值分解，利用分解得到的左、右奇异矩阵构建两个正交矩阵A_l、B_l；

S3、利用所述正交矩阵A_l、B_l建立变换结果矩阵P的更新模型，所述模型为：

目标函数：

约束条件：

采用迭代算法求解所述模型描述的最小值问题，得到更新后的变换结果矩阵P_l+1；其中，l表示更新次数，

表示所述步骤S1中预变换的反运算；λ表示稀疏项的权重，λ＞0，

表示使地震数据稀疏的变换，N＝P，A_l、B_l表示步骤S2中构建的正交矩阵；M表示观测地震数据，Ω表示未缺失数据的索引，

表示线性算子，用于使X的未缺失部分与M一致，||·||_*表示核范数运算，Tr(·)表示求矩阵的迹，||·||₁表示L₁范数运算；

S4、判断更新后的变换结果矩阵是否收敛，若收敛，则执行步骤S5，否则将所述更新后的变换结果矩阵用于下一次更新，回到步骤S2；

S5、对收敛的变换结果矩阵进行反变换得到重建的地震数据。

2.根据权利要求1所述的基于低秩和稀疏约束的地震数据重建方法，其特征在于，所述步骤S1中的预变换为纹理块变换；所述步骤S3中，

表示反纹理块变换算子。

3.根据权利要求1所述的基于低秩和稀疏约束的地震数据重建方法，其特征在于，所述步骤S2的具体过程为：对上一次更新得到的变换结果矩阵P_l进行奇异值分解，得到左奇异矩阵U_l和右奇异矩阵V_l：

[U_l，∑_l，V_l]＝svd(P_l)，

其中，

m、n分别表示矩阵P_l的行数和列数，u₁，u₂，…，u_m表示左奇异矩阵U_l的列向量，v₁，v₂，…，v_n表示右奇异矩阵V_l的列向量，

表示对应的奇异值矩阵；利用所述左奇异矩阵U_l和右奇异矩阵V_l的前r列向量，构造两个正交矩阵A_l、B_l，r＜min(m，n)且为正整数，其中：

A_l＝(u₁，u₂，…，u_r)^TB_l＝(v₁，v₂，…，v_r)^T。

4.根据权利要求1所述的基于低秩和稀疏约束的地震数据重建方法，其特征在于，所述步骤S3中，通过构造增广的拉格朗日函数来迭代求解所述模型描述的最小值问题，具体过程为：引入三个拉格朗日乘子Y、Z、F和惩罚项系数，当迭代系数k＝1时，给定迭代初值X¹、W¹、N¹、Y¹、Z¹、F¹，其中N¹＝P_l，X¹、W¹为维数大小与M相同的随机矩阵，Z¹为维数大小与M相同的零矩阵，Y¹、F¹均为零矩阵，且矩阵维度与P_l相同；在进行第k次迭代时，依次更新P^k+1、X^k+1、W^k+1、N^k+1、Y^k+1、Z^k+1、F^k+1七个参数；迭代过程中，若N^k+1满足收敛条件，则迭代结束，得到P_l+1＝P^k+1，否则递增k，继续进行迭代。

5.根据权利要求1或4所述的基于低秩和稀疏约束的地震数据重建方法，其特征在于，第k次迭代的具体过程如下：

使用奇异值收缩方法更新P^k+1：

上述运算式的计算方式为：

其中，x＝U∑V^T表示x的奇异值分解，且S_τ(∑)＝diag(max{σ_i-τ，0})，σ_i表示奇异值矩阵∑中的元素，

表示纹理块变换算子，β表示增广的拉格朗日函数中的惩罚项系数；

更新X^k+1：

其中，Ω^c表示缺失数据的索引，

表示线性算子；

更新W^k+1：

其中，上述运算式的计算方式为：

sgn(·)表示符号函数，λ表示目标函数中稀疏项的权重；

更新N^k+1：

判断N^k+1是否收敛：

||N^k+1-N^k||_F≤ε₁，

其中，||·||_F表示Frobenius范数，ε₁为给定的第一精度阈值；当满足上述条件时，停止迭代，输出P_l+1＝P^k+1；

更新Y^k+1：

Y^k+1＝Y^k+β(N^k+1-P^k+1)；

更新Z^k+1：

更新F^k+1：

6.根据权利要求1所述的基于低秩和稀疏约束的地震数据重建方法，其特征在于，所述步骤S4中，收敛条件为：

||P_l+1-P_l||_F≤ε₂，

式中，ε₂为给定的第二精度阈值；若满足所述收敛条件，则迭代结束，执行步骤S5；否则递增l，回到步骤S2。

Images