一种深度学习与预测滤波结合的地震数据反假频插值方法与流程

本发明涉及地震信号处理技术领域，尤其涉及一种深度学习与预测滤波结合的地震数据反假频插值方法。

背景技术：

物理或经济上的限制使采集的地震信号通常存在地震道缺失或空间采样间隔过大的问题，这将会影响后续的偏移成像、反演和解释。地震差值问题分为不规则缺失插值和规则缺失插值，其中，规则缺失插值对数据进行上采样，是昂贵的密集空间采样的替代方法，可以去除过大的空间采样间隔导致的频率混叠现象。

预测滤波是一类重要的地震数据反假频插值方法，例如f-x域反假频插值方法将地震记录视为线性同相轴的叠加，从低频f估计的预测误差滤波器(pef)可用于对高频2f的数据进行插值，后来该方法被推广至t-x域、f-x-y域等。然而该方法基于线性同相轴假设，而地震数据同相轴大多是非线性的，因此实际应用时要考虑加窗技术，以保证分析窗口中的数据近似线性。非平稳的滤波器是另一种处理非线性同相轴数据的插值方法，非平稳滤波器将滤波器系数的求解建模为欠定问题，通常需要引入正则化项来约束滤波器平滑变化。

近年来，机器学习和深度学习技术为插值问题提供了新的思路，利用大数据的优势可突破先验假设的限制，自动从数据中挖掘特征，极大地减少了人工工作量。例如，利用支持向量回归(svr)从训练数据中学习连续回归超平面，用以表示缺失数据与输出完整数据之间的关系。深度学习是机器学习的子类，深层神经网络相比于svr可以非线性地提取数据更高层的特征，已有的研究包括使用残差网络对规则缺失数据插值、生成对抗网络对叠后gap缺失数据插值、u-net网络对随机缺失数据插值等，这些方法构建不同的网络结构，从大量的训练数据中学习退化数据(缺失或预插值)到完整数据的映射，训练好的网络能够直接用于相似数据的插值，但难以学到一个对不同数据具有普适性的映射，同时由于深度网络本身的不可解释性，这些方法的插值原理难以解释，插值效果也难以预估。

技术实现要素：

有鉴于此，本发明提供了一种深度学习与预测滤波结合的地震数据反假频插值方法，考虑t-x域中的非平稳预测误差滤波器(predictionerrorfilter，pef)，设计网络学习退化数据和滤波器之间的关系，输入为缺失数据，输出为滤波器系数，然后使用该滤波器重建缺失的地震道，本发明同时集成了深度学习和预测滤波方法的优点，能够有效地提取地震数据的局部倾角特征，训练好的网络对不同类型的数据具有较好的适应性，计算效率更高且测试阶段不需要人工调参。

本发明提供一种深度学习和预测滤波结合的地震数据反假频插值方法，包括以下步骤：

s1、对仿真地震数据进行不同尺度的空间下采样以扩充样本集，并从扩充后的样本集中随机截取若干m×n大小的数据块，然后将所述数据块分为训练样本集和测试样本集，m表示数据块的行数，n表示数据块的列数；

s2、构建预测误差滤波器网络，所述预测误差滤波器网络包括深度学习以及预测滤波，其中，所述深度学习由残差网络实现，所述残差网络利用训练数据学习缺失的地震数据与非平稳预测滤波器之间的映射关系；所述预测滤波利用学习得到的非平稳预测滤波器对缺失的地震数据进行预测；

s3、利用步骤s1中得到的训练样本集对预测误差滤波器网络进行训练，并利用测试样本集评估所述预测误差滤波器网络的性能；利用训练完成后的预测误差滤波器网络对地震数据进行反假频插值。

进一步地，所述步骤s1中不同尺度的空间下采样包括对仿真地震数据进行2倍下采样和4倍下采样，且所述训练样本集和所述测试样本集截取自仿真地震数据的不同炮记录。

进一步地，所述步骤s2中，所述残差网络的输入为规则缺失的地震数据p(t,x)表示时间-空间域内完整的地震数据，t表示反射时间，x表示空间数据，mask表示用于构造规则缺失的地震数据的掩模算子，表示哈达玛积；所述残差网络的输出为非平稳预测滤波器的系数a，系数a的矩阵大小为(2r+1)×m×n，r表示非平稳预测滤波器的半径。

进一步地，所述步骤s2中，所述残差网络由7个残差块堆积而成，每个残差块均包含两个卷积层，所述卷积层的卷积核大小为5×5，非线性激活函数采用relu函数，并使用批标准化层加速收敛；除第一层和最后一层的卷积层外，其余卷积层的输入和输出通道数均为64；所述残差网络为一次下采样和一次上采样的编码解码结构，其中，所述下采样由步长为2的卷积层实现，所述上采样由转置卷积实现。

进一步地，所述步骤s2中，所述预测滤波中利用学习得到的非平稳预测滤波器对缺失的地震数据进行预测的具体过程为：

利用第x道地震数据对第x+1道地震数据进行预测：

其中，表示时间t处第x+1道地震数据的预测结果，p(t+k,x)表示时间t+k处第x道地震数据，a(k,t,x)表示(t,x)处的非平稳预测滤波器的第k个系数，k＝-r,-r+1,…,r。

进一步地，所述步骤s3中，利用步骤s1中得到的训练样本集对预测误差滤波器网络进行训练的具体过程为：

将训练样本集中原始完整的数据块作为标签数据y，利用掩模算子mask构造规则缺失的地震数据其中，x0为m×n的二维矩阵，然后将所述规则缺失的地震数据x0归一化到[-1,1]之间得到数据x，从而构造训练样本对(x,y)；所述归一化的方式为x＝x0/max(|xi,j|)，x表示归一化后的地震数据，x0表示输入的规则缺失的地震数据，xi,j表示x0中第i行第j列的元素，0＜i≤m，0＜j≤n，i、j为正整数。

进一步地，利用所述构造训练样本对(x,y)对所述预测误差滤波器网络进行训练，具体过程为：

向残差网络中输入数据x，并根据残差网络输出的非平稳预测滤波器系数a对标签数据y进行一次预测，得到第1步预测的结果y⁽¹⁾：

其中，0＜i≤m，0＜j≤n，i、j为正整数，y⁽¹⁾中的第1至n-1列数据是对标签数据y中第2至n列数据的预测，a(k,i,j)表示(i,j)处的非平稳预测滤波器的第k个系数，k＝-r,-r+1,…,r；

根据非平稳预测滤波器的系数a对标签数据y进行多步预测：

其中，y^(l+1)表示第l+1步预测结果，y^(l)表示第l步预测结果；

使用l1范数约束预测结果和标签数据y的距离，由此得到损失函数：

其中，s表示最多预测的地震道的道数；使用adam梯度下降法优化损失函数，完成对预测误差滤波器网络的一轮训练。所述adam梯度下降法的批尺寸为48，初始学习率为0.001，且每5轮学习率乘以0.1，共训练15轮，完成对预测误差滤波器网络的训练。

进一步地，所述步骤s3中，利用测试样本集评估所述预测误差滤波器网络的性能时，采用的评价指标为重建信噪比：

其中，表示利用所述预测误差滤波器网络对测试样本集的中测试样本进行地震数据插值的结果，yt表示测试样本集中的标签数据，表示f范数的平方。

本发明提供的技术方案带来的有益效果是：本发明同时集成了深度学习和预测滤波方法的优点，相较于其他深度学习方法，具有更好的泛化性且滤波重建过程具有可解释性，对不同数据自适应性更好；相较于传统预测滤波方法，可使用gpu并行计算，计算效率更高，无需加窗策略和人工调参。

附图说明

图1是本发明实施例提供的训练样本集中的四个示例图；

图2是本发明实施例提供的预测误差滤波器网络的结构示意图；

图3是本发明实施例提供的残差网络的结构示意图；

图4是本发明实施例提供的规则缺失的地震数据的示意图；

图5是本发明实施例提供的pefnet方法与resnet、f-x、pwd方法的插值结果对比图；

图6是本发明实施例提供的pefnet方法对真实数据的插值结果的示意图。

具体实施方式

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

发明的实施例提供了一种深度学习与预测滤波结合的地震数据反假频插值方法，采用二维仿真地震数据，包括以下步骤：

s1、根据仿真地震数据构造训练样本集和测试样本集：对仿真地震数据进行不同尺度的空间下采样以扩充样本集，并从中随机截取若干m×n大小的数据块，将其分为训练样本集以及测试样本集，m表示数据块的行数，n表示数据块的列数。

由于地震勘探成本昂贵，缺少大量的公开实际数据，且实际数据往往存在缺失，因此使用仿真数据构建训练集，本实施例步骤s2中构建的预测误差滤波器网络在仿真数据上训练完成后可直接用于实际地震数据的重建，而在有实际完整地震数据的条件下也可以进一步利用实际地震数据对所述预测误差滤波器网络进行微调。本实施例利用500炮仿真地震数据构造训练样本集以及测试样本集，为增加数据局部斜率的多样性，首先使用不同尺度的空间下采样方式对数据进行扩充，具体地，在空间上分别对地震数据进行2倍下采样和4倍下采样，得到的斜率分别为原地震数据的2倍和4倍；对于每种采样方式，从1-400炮地震数据下采样的结果中随机截取8000个m×n大小的数据块形成训练样本集，从401-500炮地震数据下采样的结果中随机截取2000个m×n大小的数据块形成测试样本集，m表示数据块的行数(即地震数据的时间)，n表示数据块的列数(即地震道的数量)。优选地，m＝n＝128。需要说明的是，本实施例还对样本集中的数据块进行50％概率的左右方向翻转，以满足训练数据的多样性。请参考图1，其为训练样本集中的四个示例，横坐标(tracenumber)表示地震道数，纵坐标(timesamplenumber)表示采样时间点。

s2、构建预测误差滤波器网络(pefnet)：所述pefnet包括深度学习部分和预测滤波部分，请参考图2，所述深度学习由残差网络f实现，残差网络f从大量的训练数据中学习缺失的地震数据与非平稳预测滤波器之间的映射关系。

具体地，残差网络f的输入为规则缺失的地震数据p(t,x)表示时间-空间(t-x)域内完整的地震数据，t表示反射时间，x表示空间数据，mask表示用于构造规则缺失的地震数据的掩模算子，本实施例中采用偶数道缺失来构造地震数据，因此mask中奇数列元素为1、偶数列元素为0，表示哈达玛积。残差网络f的输出为t-x域非平稳预测滤波器的系数a，其大小为(2r+1)×m×n，r表示非平稳预测滤波器的半径，本实施例中r＝3。

请参考图3，残差网络f由7个残差块block1、block2、……、block7堆积而成，残差块block1的输入为残差网络f的输入，此后每个残差块的输出作为下一个残差块的输入，残差块block7的输出即为残差网络f的输出；每个残差块均包含两个卷积层，卷积核大小为5×5，非线性激活函数使用relu函数，并使用标准化层(bn)加速收敛；除第一层和最后一层卷积层外，其余卷积层的输入和输出通道数均为64；残差网络f为一次下采样和一次上采样的编码解码结构，其中，残差块block1后进行下采样编码，下采样由步长为2的卷积层实现，残差块block4后进行上采样解码，上采样由转置卷积实现。

所述预测滤波使用残差网络f学习得到的非平稳预测滤波器，由第x道地震数据对第x+1道地震数据进行预测：

其中，表示时间t处第x+1道地震数据预测结果，p(t+k,x)表示时间t+k处第x道地震数据，a(k,t,x)表示(t,x)处的滤波器的第k个系数，k＝-r,-r+1,…,r。

s3、利用步骤s1中得到的训练样本集对pefnet进行训练，利用测试样本集评估网络性能；利用训练完成后的pefnet对地震数据进行反假频插值。

具体地，将训练样本集中原始完整的数据块作为标签数据y，人为将地震数据的偶数道置为0得到规则缺失的地震数据x0，其中，x0为m×n的二维矩阵，然后将所述规则缺失的地震数据x0归一化到[-1,1]之间得到数据x，从而构造训练样本对(x,y)。归一化方式为x＝x0/max(|xi,j|)，x表示归一化后的地震数据，x0表示输入的规则缺失的地震数据，xi,j表示x0中第i行第j列的元素，0＜i≤m，0＜j≤n，i、j为正整数。

利用构造的训练样本对(x,y)对pefnet进行训练，具体地，向残差网络f中输入数据x，并根据残差网络f输出的非平稳预测滤波器系数a对标签数据y进行一次预测，得到第1步预测的结果y⁽¹⁾：

其中，y⁽¹⁾中的第1至n-1列数据是对标签数据y中第2至n列数据的预测，a(k,i,j)表示(i,j)处的滤波器的第k个系数，k＝-r,-r+1,…,r；再次使用滤波器系数a对第1步预测的结果y⁽¹⁾进行预测，得到第2步预测的结果y⁽²⁾：

利用l1范数约束预测结果与标签数据的距离，由此得到残差网络f的损失函数：

最后使用adam梯度下降法优化损失函数，由此利用训练样本集对pefnet进行训练，其中，本实施例中adam梯度下降法的批尺寸为48，初始学习率为0.001，且每5轮学习率乘以0.1，共训练15轮。

利用测试样本集对pefnet的性能进行评价，评价指标采用重建信噪比：

其中，表示利用pefnet对测试样本集的中测试样本进行地震数据插值的结果，yt表示测试样本集中的标签数据，表示f范数的平方。

需要说明的是，利用训练完成后的pefnet对地震数据进行反假频插值时，首先将地震数据构造成m×n大小的数据块，并作为输入数据输入到训练完成后的pefnet中，然后利用pefnet中的残差网络f得到非平稳预测滤波器的系数，最后利用所述非平稳预测滤波器的系数通过预测滤波得到插值结果，具体地，对本实施例中偶数道缺失的地震数据进行插值重建：

其中，0＜2n≤128，n为正整数。

请参考图4，图4中(b)部分为完整地震数据，所述完整地震数据的时间和空间采样间隔分别为2ms和12.5m，共121个地震道，每个地震道包含450个时间采样点；图4中(a)部分为间隔25m的规则采样数据，在斜率较大的地方可以观察到较明显的混叠现象。

请参考图5，本实施例对比了resnet、f-x、pwd以及本发明提出的pefnet方法对图4中的仿真地震数据的插值结果，信噪比分别为17.34db、15.94db、19.38db、以及19.50db，图5中(a)、(b)、(c)部分分别为pefnet方法的重建结果、残差以及滤波器脉冲响应，图5中(d)、(g)部分分别为resnet方法的重建结果及残差，图5中(e)、(h)部分分别为f-x方法的重建结果及残差，图5中(f)、(i)部分分别为pwd方法的重建结果及残差。由图5可知，相较于其他三种方法，pefnet方法可实现较好的插值结果，且整体信噪比最高；从插值结果和残差图可以看出，f-x方法对弯曲同相轴数据的重建误差较大，且存在少量伪影；resnet方法的整体误差比pefnet大，且图5中(g)部分的箭头处可以看出，resnet方法对斜率较大的同相轴误差明显增大；图5中(i)部分的箭头处可以看出，pwd方法由于只能估计一个局部倾角，因此在两个同相轴交叉处存在明显的误差。图5中(b)部分所示的残差图可知，pefnet方法与pwd方法由于使用了类似格式的滤波器，导致pefnet方法同样在同相轴交叉处性能下降，但由于pefnet方法的滤波器具有一定自适应性，所以误差较pwd方法小。另外，由于pefnet方法的插值重建由预测滤波实现，可在一定程度上得到解释，图5中(c)部分为pefnet方法输出的滤波器在四个不同位置处的脉冲响应，即滤波器对脉冲信号的连续15道预测结果，可以看出在数据符合一个局部平面波模型时，pefnet输出的滤波器类似相位延迟算子；而在两个局部平面波叠加的情况下，pefnet在交叉点处能一定程度上的自适应提取到两个不同的局部倾角特征，表明训练好的pefnet具有一定推广性。

请参考图6，其为pefnet方法对真实数据进行处理的结果，所述数据为墨西哥湾深水区调查中获取的海洋二维炮点，包含由盐体引起的复杂衍射事件。其中，图6中(a)部分为两倍采样间隔数据，(b)部分为原始完整数据，(c)、(d)部分分别为pefnet的插值结果和残差，信噪比为17.02db，(d)部分显示的残差图表明本发明对有交叉同相轴的真实数据具有良好的自适应性。

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

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

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