地震缺失数据恢复方法、计算机存储介质和计算机系统

摘要

本发明提供了一种地震缺失数据恢复方法、计算机存储介质和计算机系统，该方法确定频率域训练集合与初始字典；利用改进字典学习算法进行字典学习，得到学习后的稀疏基；利用所得稀疏基对频率域训练集合中缺失数据进行稀疏表示，利用凸集投影算法迭代去噪，得到重建的频率域地震数据；将重建的频率域地震数据沿频率方向做傅里叶反变换，得到重建后的时间域地震数据。本发明提出的基于压缩感知的频率域地震缺失数据恢复方法，能够减少计算量，缩短计算时间，减少运行时数据占用空间，同时能够很好保护有效弱信号，使得地震数据重建结果更加可信。

说明书

技术领域

本发明属于地震勘探数据重建技术领域，尤其是一种基于压缩感知的频率域地震缺失数据恢复方法，以及计算机存储介质和计算机系统。

背景技术

随着油气勘探的深入发展，勘探开发逐渐转向复杂的油气藏，高分辨率地震勘探、时移地震等勘探技术对地震勘探数据的规则性及完整性提出了更严格的要求，实现并尽可能提高油气产量的任务，对石油资源开采具有重要的实际意义。随着科学技术的进步，大部分观测系统已经从二维转变成三维，传统的压缩手段难以满足不断增长的数据存储要求，在Nyquist采样定理下，存储野外采集的地震所需消耗很大。克服这些难题对于更加高效采集数据，以及提高深层地震资料品质有着重要意义。还有在地震勘探中，由于地形条件或者经济条件的限制，比如地表存在的障碍物或野外采集出现坏道等原因，从野外采集到的数据常常是不规则和不完整的。但是在后续的地震资料处理中，如地震子波估计、多次波去除、偏移等处理对采样率和地震数据的要求很高，因此这些采集到的不规则或者不完整的地震数据会严重影响最后的处理结果。

Spitz(1991)就提出了采用F-X域预测误差滤波的方法进行重建，但这种方法通常将不规则采样数据当作规则数据来处理，容易造成较大误差，影响插值结果的因素较多，而且对于抗假频方面效果较差，Claerbout(1991)等人则提出了T-X域的预测误差滤波重建方法，但是该种方法计算时花费的时间多，抗噪声能力弱。而后Wang(2002)在Spitz的基础上实现了在f-x-y域中的地震道插值方法，首先从给定的频率切片中估计线性预测算子，然后再预测高频信息，从而对信号进行高质量的重建，该方法可以应用到叠前也可以应到叠后处理，也可以向更高维数据重建发展。Naghizadeh(2009)提出了F-X自适应地震道插值技术，该方法采用指数权重循环最小平方来估计F-X域的自适应预测滤波器，使得能够对空间。变化的同相轴进行波场插值，不过这种方法依靠滤波器长度和遗忘因子两个参数，需要较好的选择参数，Yang(2010)提出一种正则化非稳态自回归插值重构方法，该方法通过整形正则化获取光滑的非稳态预测滤波误差系数，但该类方法计算量大且误差较大；对于基于波动方程的重建算法，Ronen(1987)利用速度信息，使用波动方程求解方法进行了插值，消除了空间假频，Canning(1996)就使用DMO算子和逆DMO算子叠加法进行三维数据的规则化重建，取得了好的效果，Fomel(2003)提出基于偏移距延拓的有限差分滤波算子对三维地震数据重构，但该类方法需要比如速度等地下结构的一些已知条件，且计算量相当大，求解时间过长，对于缺失较大的采样数据重建效果并不理想，所以目前在工业应用上较少；基于压缩感知的重建方法，首先在变换域中对地震信号进行稀疏化，然后依据信号在变换域内的特点，采用不同的数学方法对系数进行处理，再进行逆变换，通过不停的迭代，从而得到重建之后的完整数据，此类方法有很多，最早应用的是基于傅里叶变换的波场重建技术，而后其它变换方法也得到了到了大量的发展，如Radon域恢复技术、聚焦变换重建技术、Curvelet变换方法，Seislet变换方法。

发明内容

针对上述问题，本发明提出了一种新的基于压缩感知的频率域地震缺失数据恢复方法，以及计算机存储介质和计算机系统。

根据本发明的第一个方面，本发明提出的地震缺失数据恢复方法主要包括以下步骤：

S100，将待处理的地震数据从时间域变换到频率域，构成频率域训练矩阵；

S200，对频率域训练矩阵进行字典学习，以获得学习后的稀疏基，其中，所述稀疏基用于对频率域训练矩阵中的缺失数据的稀疏进行表示；

S300，根据稀疏基，重建频率域训练矩阵中的缺失数据；

S400，将重建后的频率域训练矩阵中的数据从频率域变换回时间域。

根据本发明的实施例，上述步骤S100包括以下步骤：

S110，对待处理的地震数据沿时间域进行傅里叶变换，以获得频率域地震数据；

S120，确定地震波有效信号的频带范围，根据频带范围确定相应的频率切片，根据频率切片对频率域地震数据进行分割；

S130，将分割后的频率域地震数据向量化并依此排列，以形成频率域训练矩阵。

根据本发明的实施例，上述步骤S200包括以下步骤：

S210，对频率域训练矩阵设置初始字典；

S220，利用K-SVD学习算法对初始字典进行优化处理，以获得学习后的稀疏基。

根据本发明的实施例，上述步骤S220中，所述K-SVD学习算法为基于交替最小二乘法进行最优化求解的K-SVD学习算法。

根据本发明的实施例，上述基于交替最小二乘算法进行最优化求解的K-SVD学习算法，包括以下步骤：

将下式(1)分解成式(2)和式(3)，分别对式(2)和式(3)进行求解，

式(1)中，M是频率域训练矩阵中样本的个数；n是样本的长度；D字典；x

式(2)代表稀疏表示，寻找系数矩阵X，使得DX尽可能接近Y；

式(3)代表字典的求解过程，寻找适应当前数据的字典D，将Y变换为相应的稀疏表示形式；

其中，在求解式(3)的过程中，首先将式(3)中与d

式(4)中，d

然后通过定义矩阵P

式(6)中，

将式(6)所表示的最优化问题改写成式(7)所表示的最优化问题

最后根据交替最小二乘法，将式(7)的求解转化为式(8)的两个计算式的交替迭代计算，以获得最优解，

根据本发明的实施例，在求解式(2)的过程中，采用批量正交匹配追踪算法完成矩阵逆运算。

根据本发明的实施例，上述待处理的地震数据为二维或三维含噪地震数据。

根据本发明的实施例，上述步骤S300中，根据稀疏基，通过凸集投影法迭代去噪重建频率域训练矩阵中的缺失数据。

根据本发明的另一方面，本发明还提供一种计算机存储介质，其中存储有用于实现上述方法的计算机程序。

根据本发明的另一方面，本发明还提供一种计算机系统，包括存储器和处理器，所述处理器用于执行所述存储器中存储的计算机程序，所述计算机程序用于实现上述方法。

与现有技术相比，本发明提供的地震缺失数据恢复方法具有如下优点或有益效果：

1)本发明提出的地震缺失数据恢复方法在稀疏表示过程中利用用批量正交匹配追踪算法算法来进行矩阵求逆，以减少在矩阵求逆的过程中所消耗的时间；

2)本发明提出的地震缺失数据恢复方法采用交替最小二乘算法代替奇异值分解，求解最优化问题，减少计算量；

3)本发明提出的地震缺失数据恢复方法采用频率域的训练集合代替时间域的训练集合，使最后重建地震数据更加平滑连续，从而能够避免在时间域重建结果中由于相邻时间点振幅不一致产生高频噪声以及有效信号的损失；

4)本发明提出的地震缺失数据恢复方法不仅能减少计算量、减少存储空间，还能够避免在时间域重建三维地震数据造成有效信号的损失，以及相邻时间点数据不一致而产生的高频噪声，从而获得更高质量的三维地震数据重建结果；

5)与之前的插值方法相比，本发明提出的地震缺失数据恢复方法计算效率更高、参数设置简单、不需要速度等先验信息的约束。在字典学习的基础上，能够获得适应不同类型地震数据的稀疏基，具有更强的适用性并且采集效率更高，尤其适用于地震数据重建。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

从下面描述的实施例并参考附图，本发明的其它优点和细节将变得显而易见。以下是示意图并示出：

图1为本发明实施例二的地震缺失数据恢复方法的工作流程图；

图2为本发明实施例二中的模型试验的合成地震记录：(a)原始数据(b)含噪数据(c)掩码矩阵(d)缺失数据；

图3为本发明实施例三中的基于K-SVD经典算法和K-SVD改进典算法所得字典的对比：图3(a)为基于K-SVD经典算法所得时间域字典，图3(b)、图3(c)为基于K-SVD改进算法所得频率域字典实部和虚部；

图4为本发明实施例三中模型试验得到的基于K-SVD经典算法(a)和基于K-SVD改进典算法(b)重建地震数据结果；

图5本发明实施例三中主测线20处基于K-SVD经典算法和基于K-SVD改进典算法重建地震结果和差剖面的对比：图5(a)、图5(c)为基于K-SVD经典算法时间域地震数据重建结果及误差，图5(b)、图5(d)为基于K-SVD改进典算法重建结果及误差；

图6本发明实施例三中初始字典和本发明方法所得字典，图6(a)为初始字典实部，图6(b)为初始字典虚部，图6(c)为改进算法所得字典实部，图6(d)为改进算法所得字典虚部；

图7本发明实施例三中实际数据处理结果，图7(a)为实际数据，图7(b)为50％随机采样加噪数据，图7(c)为时间域基于K-SVD经典学习算法所得的地震数据重建结果，图7(d)为频率域基于K-SVD改进学习算法所得的地震数据重建结果。

具体实施方式

实施例一

针对地震数据重建问题，本发明提出对频率域训练集合进行字典学习，在压缩感知的框架下对缺失数据进行稀疏表示，优选利用凸集投影算法对频率域三维地震数据进行重建，通过多次迭代获得完整的地震数据。

具体地，本发明提出地震缺失数据恢复方法主要包括以下步骤：

S100，将待处理的地震数据从时间域变换到频率域，构成频率域训练矩阵；

S200，对频率域训练矩阵进行字典学习，以获得学习后的稀疏基，其中，所述稀疏基用于对频率域训练矩阵中的缺失数据的稀疏进行表示；

S300，根据稀疏基，重建频率域训练矩阵中的缺失数据；

S400，将重建后的频率域训练矩阵中的数据从频率域变换回时间域。

其中，对频率域训练矩阵进行字典学习，以获得学习后的稀疏基的步骤，优选采用基于交替最小二乘法进行最优化求解的K-SVD学习算法来完成。

相比于传统的在时间域利用K-SVD学习算法(以下简称K-SVD经典算法)来实现数据重建的方式，本发明的在频率域利用基于交替最小二乘法进行最优化求解的K-SVD学习算法(以下简称K-SVD改进算法)所获得的结果会更加精确，且计算效率更高。

下面通过几个实施例来说明二种方式的不同之处以及各自所带来的优势和有益效果。

实施例二

在本实施例中，本发明提出的地震缺失数据恢复方法主要包括以下步骤：

1)选择待处理的地震数据构建训练集合，其中沿时间域对训练集合中的数据进行傅里叶变换，将这些数据从时间域变化到频率域，获得频率域地震数据；

2)确定有效信号的频带范围，根据频带范围选择对应的频率切片，利用频率切片对频率域地震数据进行分割，并将切割后的数据向量化，依次排列，构建成频率域训练矩阵；

3)为处理好的频率域训练矩阵给定一个初始字典，利用改进的K-SVD字典学习算法(即K-SVD改进算法)对初始字典进行优化，得到优化后的字典以及学习后的稀疏基；

4)学习所得稀疏基对频率域训练矩阵中的缺失数据的稀疏进行表示，利用凸集投影算法迭代去噪，重建频率域训练矩阵中的缺失数据；

5)对重建后的频率域训练矩阵进行傅里叶反变换，将其中的数据从频率域变化到时间域，以获得最终重建后的时间域的三维地震数据。

下面对上述方法的工作原理进行详细的介绍。

(一)改进的K-SVD字典学习算法

对于给定的训练集合和初始字典，求解下列最优化问题得到学习后的字典：

式中，M是训练集合中样本的个数；n是样本的长度；D字典；x

利用式(1)直接求解字典D和系数x

式(2)代表稀疏表示，寻找系数矩阵X，使得DX尽可能接近Y；式(3)代表字典的求解过程，寻找适应当前数据的字典D，将Y变换为相应的稀疏表示形式，地震数据模型复杂度随着字典更新的迭代次数的变化而减小。

通常可以利用正交匹配追踪算法对式(2)进行求解。然而，对于式(3)的求解，本实施例优选采用的是改进的K-SVD字典学习算法，按照顺序依次更新字典中的原子。具体地，考虑对第j个原子进行更新，将式(3)中与d

式中，d

在此，对原子和系数同时更新，即解决下面推导出来的优化问题：

对该优化问题的解法，常用的方法是利用秩1的逼近，通过奇异值分解(SingularValue Decomposition)，但是这会使得系数中非零个数增加不再稀疏。为了同时保持稀疏性和减少目标函数的值，本实施例优选仅使用

式中，

在上述过程中，对于最优化问题的求解是用奇异值分解，取出他们的一行和一列作为最终结果。其中，除了第一行和第一列，其他的元素都没涉及使用，因此这中间有很多的计算是不必要的。针对于此，本实施例提出不用奇异值分解进行最优化问题的求解，而是采用交替最小二乘方法代替奇异值分解求解最优化问题。换言之，对于最优化问题

的求解，采用下面这两个式子交替迭代计算：

如此，可以有效减少在奇异值分解上多余的运算，提高K-SVD算法的计算效率。

进一步地，由于利用传统的正交匹配追踪算法求解式(2)时，求解

以此避免对矩阵直接求逆。即，将矩阵求逆过程通过Cholesky分解，转化为求解上三角、下三角方程的过程，从而减少计算量，提高计算效率。

从下表表1中可以看出，在对较多的信号进行处理时，利用批量正交匹配追踪算法能够更加快速地处理数据，减少了计算所消耗的时间，提高了整体效率。

表1用K-SVD经典算法和K-SVD改进算法计算所需的时间的对比

(2)频率域地震数据重建算法

应当说明的是，本发明提出的地震缺失数据恢复方法也适用于二维和三维的含噪地震数据。

当地震数据既含有噪声，同时也存在缺失地震数据时，本发明提出的地震缺失数据恢复方法优选采用凸集投影算法迭代去噪，重建频率域训练矩阵中的缺失数据。

具体得，令信号为X，随机噪声为ε，采样矩阵为S，则含噪缺失地震数据Y：

Y＝S·×(X+ε) (10)

同样，使用上标f和下标i来表示时间切片(f为频率)和第i个数据块。

对于一个数据块的去噪重建最优化问题可以写成如下形式：

式(11)求解可以分解成：

如果利用正交匹配追踪算法求解式(12)，对于λ值的选取比较困难，其最终解相当于含噪数据与“去噪”数据对的混合，λ越小，含噪数据所占比重越小。因此，本实施例优选使用凸集投影法(POCS)迭代去噪重建地震数据。

Y

随着迭代次数的增加，λ应该从0逐渐增加到1，ρ是引进补偿能量。并且每一次迭代之后，剩余噪声能量越来越小，ε的值应该越来越小。

实施例三

下面结合具体应用实例进一步说明本发明的地震缺失数据恢复方法的有效性和进步性。

如图2所示，图2(a)为原始数据，图2(b)则是加上标准差为0.5的高斯随机噪声得到的含噪数据，图2(c)是掩码矩阵，其中黑色部分为采样部分，白色部分为缺失部分，图2(d)是含噪缺失后的地震数据。选取一个三维数据作为训练数据体进行字典学习，图3是利用K-SVD经典算法和利用K-SVD改进算法训练出来的字典，其中上部图3(a)是基于K-SVD经典算法获得的字典，从图中可以看出该算法所得字典有较多杂乱的原子，若这些原子参与到地震数据重建中，则会导致部分地震数据结果有误差；下部则是基于K-SVD改进算法获得的复数字典，其中左边图3(b)是字典的实部，右边图3(c)是字典的虚部，从图中可以看出该算法训练结果只有极少数杂乱字典。此外，在训练过程中，在同样的限制条件下，K-SVD改进算法所需要的稀疏系数的个数要比K-SVD经典算法要少。利用上述两个字典对图2中的缺失含噪地震数据进行重建，其结果如图4所示，其中左边图4(a)是基于K-SVD经典算法重建的地震数据结果，右边图4(b)是基于K-SVD改进算法重建的地震数据结果。从图中可以看出缺失数据、基于K-SVD经典算法的重建结果、基于K-SVD改进算法的重建结果的信噪比分别是-2.44dB、11.13dB、17.12dB。从信噪比的角度分析，K-SVD改进算法能够更加有效处理缺失含噪地震数据；从视觉效果的角度分析，K-SVD经典算法所得结果在无信号处依然有一定的较强的能量，并且从表面的数据来看，其结果不太连续，有较多的毛刺，而K-SVD改进算法能够很好克服这样的缺点，其结果更加清晰、平滑、连续，与原始数据最接近。

如图5所示，提取第20个主测线的剖面，从图5(a)、图5(c)中能够明显看出，基于K-SVD经典算法重建的地震数据会产生很多毛刺状数据，这可能是因为在相邻的时间切片上，同样的位置计算所得结果受噪声影响，会产生有较大的差异，而且从差剖面上看，残差的能量依然残留一些有效信号，因此虽然K-SVD经典算法也能够较好重建地震数据，但是其结果并不能让人满意。而基于K-SVD改进算法重建的地震数据图5(b)不但信噪比高，而且其残差图5(d)很小，基本看不出有效信号的损失，结果平滑连续，相较于图5(c)部分，有明显的改善。

具体地，如图6所示，图6(a)、图6(b)是初始字典的实部和虚部，图6(c)、图6(d)是基于K-SVD改进算法所得的字典的实部和虚部。从图中可以看出，改进后的字典明显比初始的字典更加平滑连续有规律。使用学习后的字典进行实际地震数据重建，其结果如图7所示，图7(a)是实际三维地震数据的二维显示，图7(b)是50％随机采样含噪数据，分别用时间域重建方法和本发明的方法进行地震数据重建，这两种方法都能够有效恢复缺失数据的能量。同样地，基于K-SVD经典算法重建的地震数据结果依然存在同样的问题：沿时间轴的数据不连续平滑，有较强的毛刺感，而频率域地震数据重建结果则平滑连续，能够更好保护有效信号能量。相应地，如图7(c)所示，时间域的地震数据重建结果依然可见地震同相轴有大量毛刺状数据，相邻的时间点数据不一致相差较大，不能很好的恢复且能量较弱的信号。相反，基于K-SVD改进算法获得的结果较为符合要求(如图7(d)所示)，不仅能够避免数据的不一致性，而且还能减少有效信号能量的损失，能够更好的恢复地震数据。

实施例四

此外，本发明还提供一种计算机存储介质，其中存储有用于实现上述方法的计算机程序。

实施例五

此外，本发明还提供一种计算机系统，其包括存储器和处理器，所述处理器用于执行所述存储器中存储的计算机程序，其中所述计算机程序用于实现上述方法。

本发明方法将提出基于压缩感知频率域数据重建，能够减少计算量，缩短计算时间，减少运行时数据占用空间，同时能够很好保护有效弱信号，使得地震数据重建结果更加可信。

应该理解的是，本发明所公开的实施例不限于这里所公开的特定处理步骤或材料，而应当延伸到相关领域的普通技术人员所理解的这些特征的等同替代。还应当理解的是，在此使用的术语仅用于描述特定实施例的目的，而并不意味着限制。

说明书中提到的“实施例”意指结合实施例描述的特定特征、或特性包括在本发明的至少一个实施例中。因此，说明书通篇各个地方出现的短语“实施例”并不一定均指同一个实施例。

本领域的技术人员应该明白，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。本领域的技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

结合本文中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储介质(RAM)、内存、只读存储介质(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM或技术领域内所公知的任意其它形式的存储介质中。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的保护范围，仍须以所附的权利要求书所界定的范围为准。