CN109490957B - 一种基于空间约束压缩感知的地震数据重建方法 - Google Patents

Abstract

本发明属于油田地震大数据重建技术领域，尤其涉及一种基于空间约束压缩感知的地震数据重建方法，包括：使用一部分数据作为训练数据，使用K‑SVD字典学***滑，说明本发明所设计的算法利用到了空间的相关信息，重建结果更加真实。

Description

一种基于空间约束压缩感知的地震数据重建方法

技术领域

本发明属于油田地震大数据重建技术领域，尤其涉及一种基于空间约束压缩感知的地震数据重建方法。

背景技术

数据重建是数据处理的重要部分。在信号领域，由于环境、设备以及人为等因素采集到的信号数据并不一定是完整的。如果使用不完整的数据进行数据解释和分析的话，则分析结果会存在较大的偏差，所以在数据解释分析之前需要对数据进行重建。另外对于地震勘探这种数据量较大的采集工作中，大量的数据会在采集、存储和运输等各个环节产生巨大的成本。因此一方面希望尽可能减少采集到的数据，另一方面希望重建出来的数据尽可能地精确。

常规的地震数据采样方法是基于Nyquist采样定理，对地震信号的采样间距具有一定的要求，如果采样频率过低则会出现假频现象，影响数据的重建。而压缩感知理论表明：基于信号的稀疏性，在低于Nyquist欠采样的情况下，对少数的采样点通过合适的重建方法仍能准确地重构信号。通常地震信号在某个变换域内是稀疏的，为利用压缩感知理论重建地震数据提供了可能。

地震数据重建即对不完整的人造地震采样数据进行插值处理，恢复出完整或者采样率更高的数据。早在1981年，Larner就对不完整的地震道恢复和野外地震数据采集设计进行了深入的讨论和研究。传统的地震数据重建方法分三类：第一类方法是基于预测滤波的方法，即采用分频预测思路，由低频信息预测高频信息。这类方法通常将非规则采样数据当作规则数据处理，并通过高斯窗进行插值，较易引入误差。第二类方法是基于波动方程的方法，即通过DMO或AMO正、反演算子迭代求解一个反问题，这类方法利用波传播的物理性质重建地震波场，但需要地下结构的先验信息，且计算量很大。第三类为基于某种变换的方法，即先对地震数据进行某种变换，然后在变换域重建。这类方法由于原理直观、计算结果稳健而得到广泛应用。

传统的方法面对低采样率和不均匀采样数据重建都是比较棘手的问题。而传统的规则均匀采样受到Nyquist采样定理的限制。而新发展起来的压缩感知理论认为即使采样频率低于Nyquist极限，也有可能恢复出满足一定精度要求的完整数据。目前压缩感知算法已经应用到了很多领域。压缩感知技术首先要求信号是稀疏的或者可压缩的，但大部分信号本身并不稀疏。不过，如果其在某个变换域内满足此条件，同样适用于压缩感知理论。自从该理论框架被提出以来，常用的变换方法主要有离散余弦变换、傅里叶变换、小波变换和曲波变换，以及逐渐采用的学习型超完备冗余字典等。

离散余弦变换(DCT)是信号处理领域最常用的变换之一，但是DCT变换是一种全局变换，并不能对图像的局部特征进行有效的识别。大部分压缩感知的应用领域都选用傅里叶变换作为稀疏变换基，但是傅里叶变换是在整个时间域内的积分，是一种全局的变换，不能很好地刻画某个局部时间的频谱特征，所以，在处理地震数据这种有明显突变现象的特征时，傅里叶变换并不是最理想的选择。Gabor提出的短时傅里叶变换，可以更好地刻画信号的局部特征，提取信号在局部时间间隔内的频谱信息。它的基本思想是通过加窗的方式将信号划分成许多小的时间间隔，然后在每一个时间窗内做傅里叶分析，以达到识别该时间间隔内局部频率的目的。

这种方法虽然在一定程度上实现了局部化，但本质上是具有单一分辨率的分析方法。但对于像地震数据这样的复杂信号，在不同时刻的波形变化较大，短时傅里叶变换的时频局部化能力还是有限的。小波分析则继承和发展了短时傅里叶变换的局部化思想，其窗口大小固定，但形状可以随着频率的变化而变换，根据频率的不同来调整时间分辨率，弥补了短时傅里叶变换的窗口大小和形状不能随频率变化的缺点。不过，小波变换不具备方向识别能力，只能捕捉点奇异特性。后来，一种被称作Curvelet(曲波)的变换被发展起来，其变换基由不同尺寸和方向的曲线状元素组成，具有多尺度和多方向识别能力，被认为是地震数据稀疏表达的最优方法之一。近期发展的Shearlet变换具有更敏感的方向性，相比Curvelet变换，能对地震信号进行更稀疏的表示，使基于压缩感知的地震数据重建效果更好，但仍然存在不能根据待处理数据自适应选择的问题。

发明内容

为了解决传统地震数据重建算法需要满足Nyquist采样定理限制的问题、使用压缩感知算法重建地震数据稀疏基难以选择和缺少帧连续性信息的问题、以及重建算法需要知道重建数据稀疏度和重建效率低的问题，本发明提出了一种基于空间约束压缩感知的地震数据重建方法，包括：

步骤1：使用一部分数据作为训练数据，使用K-SVD字典学习训练超完备字典来重建原始的地震数据；

步骤2：使用联合稀疏分解的方法，提取共有的空间信息，并改造压缩感知算法中的感知矩阵；

步骤3：对稀疏度自适应匹配追踪算法(SAMP)进行改进，引入初始稀疏度估计的方法，采用变步长的策略对数据进行重建。

所述重建原始的地震数据采用如下方法：

为原始的地震数据，

为稀疏解x的估计，

为超完备字典，稀疏解x中非零的个数K要远远小于N，通过观测矩阵Φ得到采集到的不完整的地震数据y，

y＝θx

传感矩阵

θ满足：

σ是常数且取值范围为(0，1)。

所述联合稀疏分解的方法通过将数据向量分割为共同部分和特殊部分后处理得到目标函数和约束条件为：

s为稀疏向量，δ_t+n为第n列数据的稀疏向量，y_t+n为第n列不完整的地震数据，θ_t+n为传感矩阵的第n列数据，Φ_t+n为观测矩阵的第n列数据，

为超完备字典。

所述步骤3具体包括：

设定稀疏解x中非零的个数K的初始值并判断如果

则依次增加K初始值步长K₀直到不等式不成立，同时得到绝对值最大的索引F的初始估计值F₀，其中F₀为θ^Ty中元素绝对值最大的前K₀个索引，θ为传感矩阵，y为不完整的地震数据，

为传感矩阵的第F₀列数据的转置；

引入阈值η和步长变化率λ∈(0，1)来控制步长的动态变换，当

时，

最终得到的K值为：

为稀疏估计，σ是常数且取值范围为(0，1)，v为初始步长，入为步长变化率；a_t为第t次迭代步长变化率、N为最大迭代次数、t为当前迭代轮数、

为第t轮重构的稀疏系数、

为到第t轮为止最佳的重构的稀疏系数。

5、根据权利要求4所述方法，其特征在于，所述对稀疏度自适应匹配追踪算法进行改进包括：

输入：传感矩阵θ，相关度n，观测向量{y₁，y₂，...，y_n}，迭代次数M，阈值η，初始步长v，步长变化率入；

输出：信号稀疏表示系数估计

1)构造空间传感矩阵

y＝(y₁，y₂，...，y_n)^T；

2)初始化：定义K₀为初始稀疏度，F₀是重构结果的非零项索引集合，g是中间变量，g＝A^Ty，

K₀＝1；

3)取g中K₀个最大值的索引组成F₀；

4)如果

则K₀＝K₀+1，重复(2)；

5)初始化：

初始残差r₀＝y，中间变量I＝K₀，当前轮数k＝1；

6)计算并选取|A^Tr_k-1|中I个最大值，将这些值对应A的序列号j构成集合S_k；

7)C_k＝F_k-1∪S_k，中间变量C_k中所有的序号所对应的A的列向量构成A_k；

8)求y＝A_kx_k的最小二乘解，

9)从

中选取出绝对值最大的I项记为

对应的A_k中对应的的I项记为A_kI，对应的A的序列号记为F是计算的中间结果；

10)更新残差

11)重构所得

在F处的非零项，其值为

12)如果满足终止条件，即r＝0或者达到最大迭代次数，则

输出

13)如果

那么转至(14)，否则转至(15)；

14)

重复(6)，并更新步长；

15)如果||r||₂≥||r_k-1||₂，那么I＝I+v，否则定义计算中间结果F_k＝F，r_k＝r，k＝k+1，

重复(6)，并更新k。

本发明的有益效果：

解决了传统地震数据重建算法需要满足Nyquist采样定理限制的问题；

解决了使用压缩感知算法重建地震数据稀疏基难以选择的问题；

解决了压缩感知算法重建数据缺少帧连续性信息的问题；

解决了重建算法需要知道重建数据稀疏度和重建效率低的问题。

本发明的重建的结果不但细节比较清晰，运算时间相较于IRLS和SAMP大幅地降低，而且横向的过度更加的平滑，说明本发明所设计的算法利用到了空间的相关信息，重建结果更加真实。

附图说明

图1为本发明的方法流程图

图2a为原始地震数据

图2b为50％的随机抽道后的原始地震数据

图3为使用本发明所设计的方法对50％抽道数据进行重建后的结果；

图4多种算法与本发明的方法的对照试验重建结果

图5为对图4局部放大图。

具体实施方式

下面结合附图，对实施例作详细说明。

A.具有空间相关性的压缩感知算法

基于压缩感知的地震数据重建可以表示为：

y＝Φf (1)

式中：y∈R^M为采集到的不完整的地震数据，f∈R^N为原始的完整的地震数据(M＜N)。Φ∈R^M×N为观测矩阵。使用超完备字典

对完整的地震数据f进行稀疏表示，则可以表示为：

式中，稀疏解x中非零的个数K要远远小于N，再通过观测矩阵Φ得到采集到的不完整的地震数据y，表示为

y＝θx (3)

式中，传感矩阵

Φ和

不相关，最后进行地震数据的重建，即

是x的估计。最终通过以下公式来重建原始的地震数据

为保证K稀疏解x的准确性，θ需要满足RIP，即

式中σ是常数且取值范围为(0，1)。

分析以上公式可知，传统的压缩感知算法是对一维向量f进行的恢复。但是地震切片数据是二维的，二维空间属性之间存在一定的相关性，空间上越相近的两点属性值越接近。假设

与

为重建后相邻的两列数据，则可以表示为：

我们采用联合稀疏分解的方式将

与

分割为一个共同部分和两个特殊部分，即：

其中μ是

和

的共同部分，为基数

和稀疏向量s的积。两个特殊部分分别表示为v_t和v_t+1，可以表示为基数

和稀疏向量δ_t和δ_t+1的积。则式(7)的目标函数可以修改为：

min||δ_t-δ_t+1||₁ (9)

由于δ_t与δ_t+1中没有公共部分，则式(9)可以转变为：

min||δ_t，δ_t+1||₁ (10)

将式(10)与式(4)合并可得：

上式可以近似表达为：

可以看出，上式满足压缩感知的表示形式，可以使用压缩感知的重建算法来进行求解。我们可以将上式推广到n个相关向量形式，即：

根据式(13)，我们在求

的过程中保证非零项尽可能地少，那么就会导致非零项尽可能地集中在共有部分s中，这样就可以保证空间上相邻的向量(m₁，m₂，...，m_n)尽可能地接近，满足了空间上相邻的向量尽可能相似的条件。

B.自适应动态步长SAMP算法

实际勘探到的地震数据稀疏度往往是未知的，SAMP算法在稀疏的未知的情况下通过设定固定步长s，每轮迭代中根据残差值的比对情况，决定信号运算过程当中稀疏度的估计值是否增大，从而近似求出信号稀疏度的最优估计，最终实现压缩结果。该算法流程如附图1所示。

分析SAMP算法可知，如果步长s初始值过小，则迭代次数明显会增大，从而算法执行时间过长；而步长s如果设置过大有可能会跳过信号真实稀疏度的值，这样会导致重构精度大打折扣。本发明结合地震数据的具体特点，提出了针对地震数据具有更好效果的基于空间序列的SAMP改进算法。

根据相关理论可得，当观测矩阵0以参数(K，σ)满足RIP性质时，如果K₀≥K，则

其中F₀为θ^Ty中元素绝对值最大的前K₀个索引。

根据以上结论的逆否命题可得，当

时，K₀＜K。从这个命题可以得到对K初始估计方法：K₀取初始值1，如果

则依次增加K₀直到不等式不成立，同时得到F的初始估计值F₀。

本发明引入阈值η和步长变化率λ∈(0，1)来控制步长的动态变换，当

时，

最终得到的K值为：

结合上述内容，本发明所涉及的基于空间序列的SAMP改进算法步骤如下：

输入：传感矩阵θ，相关度n，观测向量{y₁，y₂，...，y_n}，迭代次数M，阈值η，初始步长v，步长变化率λ；

输出：信号稀疏表示系数估计

(1)根据公式(13)构造空间传感矩阵

y＝(y₁，y₂，...，y_n)^T。

(2)初始化：g＝A^Ty，

K₀＝1。

(3)取g中K₀个最大值的索引组成F₀。

(4)如果

则K₀＝K₀+1，重复(2)

(5)初始化：

r₀＝y，I＝K₀，k＝1。

(6)计算并选取|A^Tr_k-1|中I个最大值，将这些值对应A的序列号j构成集合S_k。

(7)C_k＝F_k-1∪S_k，C_k中所有的序号所对应的A的列向量构成A_k。

(8)求y＝A_kx_k的最小二乘解，

(9)从

中选取出绝对值最大的I项记为

对应的A_k中对应的的I项记为A_kI，对应的A的序列号记为F。

(10)更新残差

(11)重构所得

在F处的非零项，其值为

(12)如果满足终止条件(r＝0或者达到最大迭代次数)，

输出

(13)如果

那么转至(14)，否则转至(15)。

(14)

(15)如果||r||₂≥||r_k-1||₂，那么I＝I+v，否则F_k＝F，r_k＝r，k＝k+1，

重复(6)。

为了验证本发明设计的重建算法的可行性以及有效性，对真实的地震数据进行数值实验。衡量数据重建效果的指标分别为信噪比(SNR)和峰值信噪比(PSNR)，即：

式中：y为原始数据；

为重建后的数据；MSE为原始数据与重建数据的均方误差。由公式(11，12)可知，SNR和PSNR越大代表误差越小，重建的效果越好。

附图2(a)为原始地震数据，数据大小为601*626。为验证本发明设计的重构算法的效果，对它进行50％的随机抽道，附图2(b)所示。

在数据恢复之前，首先使用K-SVD字典学习得到一个601*1052的超完备字典，然后使用学习到的超完备字典作为稀疏基矩阵，使用本发明所设计的算法对50％抽道数据进行重建，重建结果如附图3所示。

计算本发明设计的方法重建结果的SNR和PSNR：SNR为10.514，PSNR为49.044。本发明设计的算法取得了良好重建精度效果。

与地震数据残缺部分非常均匀不同，测井的电成像数据残缺部分具有残缺部分比较宽，而且不规则的特点。为了证明本发明所设计的算法可以应对此类问题，对一段真实的测井的电成像数据进行重建，原始数据大小为360*1000，分别使用SP、OMP、CoSaMP、IHT、IRLS、SAMP算法与本发明所设计的算法做对照试验，其中SP、OMP、CoSaMp、IHT算法的K值设定为50，SAMP的步长设定为5。重建结果如附图4所示。对局部进行放大，显示结果附图5所示。

从恢复效果可以看出，由于OMP、SP、CoSaMP、IHT算法需要预先设定稀疏度，在稀疏度设定不合适的时候重建结果略显模糊。而实际的电成像问题中稀疏度并不知道，因此这些算法并不适合此类问题。IRLS、SAMP算法重建结果细节比较清晰，但是需要的运算时间过长。而且这六种重建方法共有的问题是重建的图像在残缺部分都有明显的重建痕迹。而发明所设计的算法重建的结果不但细节比较清晰，运算时间相较于IRLS和SAMP大幅地降低，而且横向的过度更加的平滑，说明本发明所设计的算法利用到了空间的相关信息，重建结果更加真实。

此实施例仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (1)

1.一种基于空间约束压缩感知的地震数据重建方法，其特征在于，包括：

步骤1：使用一部分数据作为训练数据，使用K-SVD字典学习训练超完备字典来重建原始的地震数据；

步骤2：使用联合稀疏分解的方法，提取共有的空间信息，并改造压缩感知算法中的感知矩阵；

步骤3：对稀疏度自适应匹配追踪算法(SAMP)进行改进，引入初始稀疏度估计的方法，采用变步长的策略对数据进行重建；

所述重建原始的地震数据采用如下方法：

为原始的地震数据，

为稀疏解x的估计，

为超完备字典，稀疏解x中非零的个数K要远远小于N，通过观测矩阵Φ得到采集到的不完整的地震数据y，

y＝θx

传感矩阵

θ满足：

σ是常数且取值范围为(0，1)；

所述联合稀疏分解的方法通过将数据向量分割为共同部分和特殊部分后处理得到目标函数和约束条件为：

s为稀疏向量，δ_t+n为第n列数据的稀疏向量，y_t+n为第n列不完整的地震数据，θ_t+n为传感矩阵的第n列数据，Φ_t+n为观测矩阵的第n列数据，

为超完备字典；

所述步骤3具体包括：

设定稀疏解x中非零的个数K的初始值并判断如果

则依次增加K初始值步长K₀直到不等式不成立，同时得到绝对值最大的索引F的初始估计值F₀，其中F₀为θ^Ty中元素绝对值最大的前K₀个索引，θ为传感矩阵，y为不完整的地震数据，

为传感矩阵的第F₀列数据的转置；

引入阈值η和步长变化率λ∈(0，1)来控制步长的动态变换，当

时，

最终得到的K值为：

为稀疏估计，σ是常数且取值范围为(0，1)，v为初始步长，λ为步长变化率；a_t为第t次迭代步长变化率、N为最大迭代次数、t为当前迭代轮数、

为第t轮重构的稀疏系数、

为到第t轮为止最佳的重构的稀疏系数；

所述对稀疏度自适应匹配追踪算法进行改进包括：

输入：传感矩阵θ，相关度n，观测向量{y₁，y₂，…，y_n}，迭代次数M，阈值η，初始步长v，步长变化率λ；

输出：信号稀疏表示系数估计

1)构造空间传感矩阵

y＝(y₁，y₂，...，y_n)^T；

2)初始化：定义K₀为初始稀疏度，F₀是重构结果的非零项索引集合，g是中间变量，g＝A^Ty，

K₀＝1；

3)取g中K₀个最大值的索引组成F₀；

4)如果

则K₀＝K₀+1，重复2)；

5)初始化：

初始残差r₀＝y，中间变量I＝K₀，当前轮数k＝1；

6)计算并选取|A^Tr_k-1|中I个最大值，将这些值对应A的序列号j构成集合S_k；

7)C_k＝F_k-1∪S_k，中间变量C_k中所有的序号所对应的A的列向量构成A_k；

8)求y＝A_kx_k的最小二乘解，

9)从

中选取出绝对值最大的I项记为

对应的A_k中对应的的I项记为A_kI，对应的A的序列号记为F是计算的中间结果；

10)更新残差

11)重构所得

在F处的非零项，其值为

12)如果满足终止条件，即r＝0或者达到最大迭代次数，则

输出

13)如果

那么转至14)，否则转至15)；

14)

重复6)，并更新步长；

15)如果||r||₂≥||r_k-1||₂，那么I＝I+v，否则定义计算中间结果F_k＝F，r_k＝r，k=k+1，

重复6)，并更新k。

Images