CN110161565B - 一种地震数据重建方法 - Google Patents

Abstract

本发明公开了一种地震数据重建方法，包括以下步骤：获取待重建的地震数据及其采集矩阵；对地震数据和采集矩阵进行分块处理；建立每个小块的地震数据及其对应的采集矩阵的双曲正切函数正则化的稀疏反演模型；通过每个小块数据求取该稀疏反演模型的最优解，得到重建后的地震数据；对所有分块数据进行数据重建，并将所有重建的数据拼接成完整的重建数据。本发明公开的地震数据重建方法计算速率快、精度高，不仅适用于二维数据，同样适用于三维数据和五维数据；不仅适用于叠前地震数据，也适用于叠后地震数据；不仅适用于陆上地震勘探数据，也适用于海洋地震勘探数据。

Description

一种地震数据重建方法

技术领域

本发明涉及地震数据处理技术领域，具体涉及一种高时效、高精度的地震数据重建方法。

背景技术

地震勘探是人类获取地下信息的重要手段，在油气勘探、煤田勘探、环境与工程勘探等领域具有重要作用。在地震勘探的数据采集过程中，地震道缺失、采样数据不规则、不满足采集设计方案与采样定理的情形时有发生。其产生原因包括坏道、地表障碍物影响或在预处理过程中剔除含强烈噪声的道记录等，这些均会导致数据缺失，数据缺失又会导致采样假频，使得地震道中出现假象，从而影响地震数据处理过程中的多次波消除、振幅随偏移距变化分析、偏移成像和全波形反演等关键处理环节的效果，进而最终影响地震勘探的精度。在此背景下，地震数据重建的研究逐渐兴起。对采样不完备的数据进行地震数据重建，获得满足采样定理的数据，对于地震数据处理和解释的精准非常重要。

地震数据重建方法主要有基于预测误差滤波的方法和基于矩阵和张量完备的方法等。其中，基于预测误差滤波的方法利用不含假频的低频数据预测高频数据，在傅里叶变换域实现数据重建，但是这类方法需要假设地震数据的同相轴是线性的，并且该方法适用于规则采样不足的数据。基于矩阵和张量完备的方法假设完整采样的地震数据经过某种重排后的矩阵的奇异值分解得到奇异值很少，而采样数据经过同样方式重排得到的矩阵的奇异值不是稀疏的，这类方法同样假设同相轴是线性的。目前，地震数据的处理工作具有数据量大和计算量大的特点，目前的地震数据重建过程面临着计算速度慢、重建数据精度低等问题。

发明内容

本发明的目的在于提供一种地震数据重建方法，用以解决目前的地震数据重建计算速度慢、重建数据精度低的问题。

本发明提供一种地震数据重建方法，包括以下步骤：

步骤A：获取待重建的地震数据及其采集矩阵；

步骤B：对地震数据和采集矩阵进行分块处理；

步骤C：建立每个小块的地震数据及其对应的采集矩阵的双曲正切函数正则化的稀疏反演模型，其中，所述双曲正切函数正则化的稀疏反演模型为：

式中，min表示最小化符号，J_σ(x)为目标函数，F^-1为傅里叶变换的逆，x为地震数据在傅里叶变换域的系数，:＝表示将minJ_σ(x)定义为

s.t.表示前面目标函数最小化的约束条件为SF^-1x＝d，σ为超参数，d为待重建的地震数据，S为采集矩阵，

表示地震数据在傅里叶变换域系数的第i个分量的一种双曲正切函数； N为地震数据在傅里叶变换域的系数x的分量个数，xi表示地震数据在傅里叶变换域的系数x的第i个分量；

步骤D：通过每个小块数据求取该稀疏反演模型的最优解，得到重建后的地震数据；

步骤E：对所有分块数据进行数据重建，并将所有重建的数据拼接成完整的重建数据。

作为一种具体的技术方案，所述步骤D包括以下步骤：

步骤D11：计算该稀疏反演模型当前迭代计算中迭代解的梯度；

步骤D12：判断该梯度的修正牛顿方向；

步骤D13：根据该修正牛顿方向更新迭代解；

步骤D14：计算该更新后的解在时间空间域的投影，得到下一次迭代计算的新迭代解，直至当满足预设迭代要求或迭代次数阈值时，停止迭代并以该新迭代解为最优解。

进一步地，所述步骤D11包括以下步骤：

步骤D111：计算当前迭代计算中迭代解在变换域的系数；

步骤D112：根据该系数计算该迭代解的梯度。

作为上述步骤D的一种等价描述，所述步骤D包括以下步骤：

步骤D21：设定循环次数IterN、初始超参数σ₀和终止迭代超参数σ_IterN，并将观测得到的地震数据作为初始解d₀；

步骤D22：计算当前第k次迭代的迭代解在变换域的系数x_k，其表达式为：

x_k＝Fd_k (2)

式中，k表示迭代序号；F表示傅里叶变换；当k＝0时，表示初始迭代，d_k表示第k次迭代时重建所得的地震数据；

计算当前第k次迭代的修正牛顿方向g_k，其表达式为：

式中，σ_k为第k次迭代的超参数；

并求出更新解：

式中，λ为步长；

步骤D23：求更新解

在时间空间域的投影

进行k的自增1运算，如果k＜IterN，重复所述步骤D22到所述步骤D23的迭代过程；如果k≥IterN，则迭代结束；

步骤D24：输出小块数据的最后一次迭代时重建所得的地震数据d_IterN。

本发明的有益效果是：

本发明公开的地震数据重建方法，充分考虑到缺失道对一些多道处理环节的影响，采用快速傅里叶变换实现采样数据的稀疏变换，能够快速实现地震数据重建；建立的双曲正切函数正则化的稀疏反演模型可以提高地震数据重建的精度；采用的修正牛顿- 投影法在求解建立稀疏反演模型时，比求解基于一范数正则化的反演模型需要更少的迭代次数，能够更加快速的得到重建的结果；对地震数据进行分块处理，对每个小块进行地震数据重建，可以提高地震数据重建的计算效率。本发明公开的地震数据重建方法计算速率快、计算精度高，不仅适用于二维数据，同样适用于三维数据和五维数据；不仅适用于叠前地震数据，也适用于叠后地震数据；不仅适用于陆上地震勘探数据，也适用于海洋地震勘探数据。

附图说明

图1为实施例1提供的原始Marmousi速度模型正演得到的炮记录和采样的结果；

图2(a)为实施例1提供的采用凸集投影法50次迭代后的重建结果；

图2(b)为实施例1提供的采用本发明涉及的地震数据重建法50次迭代后的重建结果；

图3(a)为实施例1提供的原始数据第三道和凸集投影法重建数据的第三道的比较；

图3(b)为实施例1提供的原始数据第三道和本发明涉及的地震数据重建法重建数据的第三道的比较；

图4(a)为实施例1提供的一个真实的地震数据；

图4(b)为实施例1提供的采样二分之一的地震数据；

图5(a)为实施例1提供的采用凸集投影法20次迭代的重建结果；

图5(b )为实施例1提供的采用本发明涉及的地震数据重建法20次迭代的重建结果；

图6为实施例1提供的原始数据、凸集投影法重建数据和本发明涉及的地震数据重建法重建数据的第三道的比较。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

实施例1提供了一种地震数据重建方法，该地震数据重建方法包括以下步骤：

步骤A：获取待重建的地震数据及其采集矩阵；

这里，待重建的地震数据可以是采集的原始地震数据经预处理后得到的地震数据，也可以是经过叠加、偏移等处理后的数据。首先，需要获取该待重建的地震数据，以及地震数据相对应的采集矩阵。

步骤B：对地震数据和采集矩阵进行分块处理；

这里，采用的分块方式可以是有重叠的分块，也可以是无重叠的分块。对待重建的地震数据和采集矩阵采用相同的分块方式。

步骤C：建立每个小块的地震数据及其对应的采集矩阵的双曲正切函数正则化的稀疏反演模型；

其中，该双曲正切函数正则化的稀疏反演模型为：

式中，min表示最小化符号，J_σ(x)为目标函数，F^-1为傅里叶变换的逆，x为地震数据在傅里叶变换域的系数，:＝表示将minJ_σ(x)定义为

s.t.表示前面目标函数最小化的约束条件为SF^-1x＝d，σ为超参数，d为待重建的地震数据，S为采集矩阵，

表示地震数据在傅里叶变换域系数的第i个分量的一种双曲正切函数；N为地震数据在傅里叶变换域的系数x的分量个数，x_i表示地震数据在傅里叶变换域的系数x的第i个分量；

步骤D：通过每个小块数据求取该稀疏反演模型的最优解，得到重建后的地震数据：

作为一种具体地实施方式，包括以下步骤：

步骤D11：计算该稀疏反演模型当前迭代计算中迭代解的梯度；

步骤D12：判断该梯度的修正牛顿方向；

步骤D13：根据该修正牛顿方向更新迭代解；

步骤D14：计算该更新后的解在时间空间域的投影，得到下一次迭代计算的新迭代解，直至当满足预设迭代要求时，停止迭代并以该新迭代解为最优解，或者迭代计算的次数达到预设的迭代次数阈值时，停止迭代并以该新迭代解为最优解。

在一种可能的实施方式中，当迭代计算的次数达到预设的迭代次数阈值时，停止迭代并以该新迭代解为最优解。在另一种方式中，也可以设置最小超参数，当迭代计算中得到的超参数减小到与该最小超参数相等时，停止迭代计算，并以此次迭代得到的新迭代解为最优解。即当迭代完成后，最后一次迭代计算中得到的新迭代解即为最优解，作为最终的重建地震数据。

进一步地，步骤D11还包括以下步骤：

步骤D111：计算当前迭代计算中迭代解在变换域的系数；

步骤D112：根据该系数计算该迭代解的梯度。

步骤E：对所有分块数据进行数据重建，并将所有重建的数据拼接成完整的重建数据。

这样，即实现了本发明涉及的的高效、高精度地震数据重建方法。该方法通过双曲正切函数建立稀疏反演模型并求解，在时间空间域中求取基于修正牛顿方向的预测解，然后利用投影运算获得迭代解，并通过多次迭代获得最终地震数据重建结果。

为验证本发明方法的可行性和有效性，以下举两个应用实例。

应用实例1：

由图1(a)可知，经Marmousi模型正演得到的地震炮记录，共含有96道，时间采样为600，采样率为4毫秒，对该数据进行60％的随机采样得到的数据见图1(b)。

这里，分别采用经典的凸集投影法和本发明公开的地震数据重建方法对上述地震数据进行重建，得到结果分别见图2(a)和图2(b)，其中，图2(a)为凸集投影法50次迭代的结果，凸集投影法选择最小阈值为0.001，其计算时间为0.177675秒，图2(a) 的信噪比为14.7561dB。图2(b)为本发明涉及的地震数据重建方法的结果，参数σ_min＝0.001，内部迭代次数1，外部迭代次数为50次，计算时间为0.199839秒，图2(b) 的信噪比为16.8293dB。图3(a)和图3(b)为第三道数据的重建结果比较。可见在相同的迭代次数、相似的计算时间条件下，本发明方法能够取得更高的计算精度。

表1两种方法经不同迭代次数所得重建结果信噪比的比较

通过对比可见，在取得相当的重建结果的前提下，本发明涉及的地震数据重建法的迭代次数和计算时间大约为基于一范数的凸集投影法的一半。

应用实例2：

如图4(a)所示，某地区的真实地震数据的一部分的道数为100，时间采样个数为100，时间采样率为4毫秒。对该数据采样50％的结果见图4(b)。基于一范数正则化的凸集投影法计算20次迭代，得到的结果见图5(a)，其信噪比为16.5436dB，计算时间为0.018395s。本发明涉及的地震数据重建法计算20次迭代得到的结果为图5(b)，信噪比为24.0718dB，计算时间为0.019996s。

如图6为凸集投影法和本发明涉及的地震数据重建法重建数据的第三道(缺失道)的原始数据的比较。由图6可知，在相同的迭代次数、相似的计算时间条件下，本发明涉及的地震数据重建法能够取得更高的计算精度。

实施例2

实施例2提供了一种地震数据重建方法，在实施例1的基础上，实施例2对实施例1的步骤D做了另一种等价表述：

步骤D采用修正牛顿-投影法求取该稀疏反演模型的最优解，具体的步骤包括：

步骤D21：设定循环次数IterN、、初始超参数σ₀和终止迭代超参数σ_IterN，并将观测得到的地震数据作为初始解d₀；

步骤D22：计算当前第k次迭代的迭代解在变换域的系数x_k＝Fd_k，式中，k表示迭代序号；F表示傅里叶变换；当k＝0时，表示初始迭代，d_k表示第k次迭代时重建所得的地震数据；

计算当前第k次迭代的修正牛顿方向：

式中，σ_k为第k次迭代的超参数；

并求出更新解：

λ为步长；

步骤D23：求更新解

在时间空间域的投影

进行k的自增1运算，如果k＜IterN，重复所述步骤D22到所述步骤D23的迭代过程；如果k≥IterN，则迭代结束；

步骤D24：输出小块数据的最后一次迭代时重建所得的地震数据d_IterN。

本发明公开的地震数据重建方法，充分考虑到缺失道对一些多道处理环节(如波动方程偏移、表面多次波消除和谱估计等)的影响，在偏移之前加密空间方向的采样，防止因采样过疏而在偏移过程中产生空间假频，以双曲正切函数作为正则化函数建立稀疏反演模型，通过修正牛顿-投影法求解稀疏反演模型实现地震数据重建，提高了地震数据重建的精度。本发明公开的地震数据重建方法计算速率快、精度高，将会在地震勘探中发挥重要作用，可广泛应用于油气勘探、煤田勘探、环境与工程勘探等领域中。

虽然，上文中已经用一般性说明及具体实施例对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (1)

1.一种地震数据重建方法，其特征在于，包括以下步骤：

步骤A：获取待重建的地震数据及其采集矩阵；

步骤B：对地震数据和采集矩阵进行分块处理；

步骤C：建立每个小块的地震数据及其对应的采集矩阵的双曲正切函数正则化的稀疏反演模型，其中，所述双曲正切函数正则化的稀疏反演模型为：

式中，min表示最小化符号，J_σ(x)为目标函数，F^-1为傅里叶变换的逆，x为地震数据在傅里叶变换域的系数，:＝表示将minJ_σ(x)定义为

s.t.表示前面目标函数最小化的约束条件为SF^-1x＝d，其中σ为超参数，d为待重建的地震数据，S为采集矩阵，

表示地震数据在傅里叶变换域系数的第i个分量的一种双曲正切函数；N为地震数据在傅里叶变换域的系数x的分量个数，x_i表示地震数据在傅里叶变换域的系数x的第i个分量；

步骤D：通过每个小块数据求取该稀疏反演模型的最优解，得到重建后的地震数据；

步骤D21：设定循环次数IterN、初始超参数σ₀和终止迭代超参数σ_IterN，并将观测得到的地震数据作为初始解d₀；

步骤D22：计算当前第k次迭代的迭代解在变换域的系数x_k，其表达式为：

x_k＝Fd_k (2)

式中，k表示迭代序号；F表示傅里叶变换；当k＝0时，表示初始迭代，d_k表示第k次迭代时重建所得的地震数据；

计算当前第k次迭代的修正牛顿方向g_k，其表达式为：

式中，σ_k为第k次迭代的超参数；

求出更新解：

式中，λ为步长；

步骤D23：更新解

在时间空间域的投影：

进行k的自增1运算，如果k＜IterN，重复所述步骤D22到所述步骤D23的迭代过程；如果k≥IterN，则迭代结束；

步骤D24：输出小块数据的最后一次迭代时重建所得的地震数据d_IterN；

步骤E：对所有分块数据进行数据重建，并将所有重建的数据拼接成完整的重建数据。

Images