0 引言

在地震数据现场采集过程中，受水库、大堤、村庄、矿场等障碍物或禁采区的影响，地震数据往往是不完整或不规则的。如果不对缺失的数据进行恢复，会影响偏移成像和油藏描述的精度。因此，缺失数据的重建在地震勘探中具有重要意义。

目前，缺失数据重建的方法主要包括基于预测滤波的重建方法、基于波动方程的重建方法、基于降秩正则化的重建方法以及基于稀疏变换的压缩感知重建方法四类。基于预测滤波的重建方法是利用分频策略，以低频信息预测并构建高频信息，经典方法包含F-X域预测滤波^[1]、F-K域预测滤波^[2]等，此类方法通常用于规则采样数据的插值加密。基于波动方程的重建方法是通过正、反算子实现地震数据恢复^[3-4]，此类方法虽可以对非规则采样的数据进行重建，但仍受限于速度场的准确程度以及全波场计算的巨大运算量。基于降秩正则化的重建方法是假设地震数据具有低秩特征，数据的缺失和含噪会破坏该特征，由此通过降秩算法实现地震数据的重建^[5-6]。基于稀疏变换的压缩感知重建方法是将数据空间投影到模型空间，在保证模型空间数据稀疏的前提下，通过求解反问题实现缺失地震道的重建。常用的稀疏变换包含傅里叶变换^[7-8]、小波变换^[9]、Radon变换^[10-12]、Curvelet变换^[13-17]、Seislet变换^[18]、Shearlet变换^[19]等，上述变换的基函数都是固定的，往往不能根据地震数据的特征进行自适应调整。随着压缩感知技术的发展，提出了自适应学习字典的新型稀疏变换^[20-22]。与固定基不同，此类变换的基函数是从数据自身中学习而来，因其在学习过程中能够揭示数据内部的相似性，故相较于固定基具有更强的稀疏表示能力，更有利于地震数据重建^[23]。K-奇异值分解(K-SVD)是最具代表性的字典学习算法，周亚同等^[24]和李慧等^[25]将其应用到了地震数据重建方面，取得了不错的应用效果。但K-SVD算法具有较高的计算复杂度，难以训练维数较高的字典，不利于高维数据的处理。为解决上述问题，Li等^[26]提出了适合高维数据处理的快速结构字典学习算法，并将其应用于三维医学图像的去噪。

本文将快速结构字典学习引入三维地震数据重建，提出了一种基于快速结构字典学习的稀疏域三维地震数据重建方法。该方法通过快速结构字典学习算法训练训练集产生三维自适应字典，并将三维自适应字典作为稀疏基对地震数据进行稀疏表示，之后根据三维地震数据缺失道的分布情况构建观测矩阵，在重建阶段使用正则化正交匹配追踪重构算法对数据进行保幅重建。该方法的优点在于采用的是三维自适应字典，具有更强的稀疏表示能力，可以充分利用三维数据空间上的连续性约束，减少重建误差。

1 方法原理 1.1 压缩感知地震数据重建模型

地震数据重建模型可以描述为

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{Rf}} $

(1)

式中：y为观测数据，即含缺失道的地震数据；R是元素为0和1组成的观测矩阵，采样点处为1，缺失点处为0；f为完整数据。根据压缩感知理论，式(1)可以改写为

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{D}}^{\rm{H}}}\mathit{\boldsymbol{x}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} x}} $

(2)

式中：Φ=RD^H为压缩感知信息算子；D^H表示稀疏基D的共轭转置；x=Df为稀疏域中的系数。由于x是稀疏的，重建问题可以转化为最小L₀范数的优化问题

$ \mathit{\boldsymbol{\tilde x}} = \min {\left\| \mathit{\boldsymbol{x}} \right\|_0}\;\;\;满足\;\;\;\mathit{\boldsymbol{y}} = \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{D}}^{\rm{H}}}\mathit{\boldsymbol{x}} $

(3)

本文采用正则化正交匹配追踪算法^[27]对式(3)进行稀疏迭代求解，该算法相较于其他追踪算法具有运算效率高及重建精度高的优势。

在估算出系数$\mathit{\boldsymbol{\tilde x}}$后，由下式进行重建

$ \mathit{\boldsymbol{\tilde f}} = {\mathit{\boldsymbol{D}}^{\rm{H}}}\mathit{\boldsymbol{\tilde x}} $

(4)

式中$\mathit{\boldsymbol{\tilde f}}$为重建后的数据。

1.2 快速结构字典学习

字典学习的本质是从原始数据中提取一个训练集，根据稀疏约束条件从训练集中训练出与原始数据特征相似的超完备字典，再使用该字典对原始数据进行稀疏表示。由于训练的超完备字典是完全自适应的，故其相较于固定字典具有更强的稀疏表示能力。快速结构字典学习结合了多簇追踪(Multiple Cluster Pursuit，MCP)稀疏分解算法和最小二乘快速逼近SVD算法的优点，具有计算效率高、适合高维数据稀疏表示的特点。快速结构字典学习目标函数定义为

$ \mathop {\min }\limits_{\mathit{\boldsymbol{D}},\mathit{\boldsymbol{X}}} \left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{DX}}} \right\|_{\rm{F}}^2\;\;\;满足\;\;\;{\left\| \mathit{\boldsymbol{X}} \right\|_p} \le T $

(5)

式中：D为训练的超完备自适应字典；Y是训练集向量化后的矩阵；X是稀疏系数矩阵；‖·‖_F表示Frobenius范数；‖·‖_p表示L_p范数(0≤p≤1)；T为稀疏系数中非零系数的数量上限值。由于稀疏系数X和训练集Y为矩阵形式，故对于三维地震数据的处理只需将其向量化后组成矩阵便可实现。

快速结构化字典学习方法的实施步骤如下。

(1) 稀疏编码：固定字典D后，使用MCP稀疏分解算法逐列更新系数矩阵X

$ \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_i}} \left\{ {\left\| {{\mathit{\boldsymbol{y}}_i} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{x}}_i}} \right\|_2^2} \right\}\;\;\;\;满足\;\;\;{\left\| {{\mathit{\boldsymbol{x}}_i}} \right\|_p} \le T $

(6)

式中：y_i为训练集Y中的第i列向量；x_i为稀疏系数矩阵X中的第i列向量。

(2) 字典更新：固定X后采用逼近奇异值分解(SVD)的策略更新字典D

$ \mathop {\min }\limits_\mathit{\boldsymbol{D}} \left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{DX}}} \right\|_{\rm{F}}^2 $

(7)

为了求解方便，式(7)可以改写为

$ \begin{array}{l} \left\{ {{\mathit{\boldsymbol{d}}_h},\mathit{\boldsymbol{x}}_h^{\rm{T}}} \right\} = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{d}}_h},\mathit{\boldsymbol{x}}_h^{\rm{T}}} \left\| {\mathit{\boldsymbol{Y}} - \sum\limits_{j \ne h} {{\mathit{\boldsymbol{d}}_j}\mathit{\boldsymbol{x}}_j^{\rm{T}}} - {\mathit{\boldsymbol{d}}_h}\mathit{\boldsymbol{x}}_h^{\rm{T}}} \right\|_{\rm{F}}^2\\ \;\;\;\;\; = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{d}}_h},\mathit{\boldsymbol{x}}_h^{\rm{T}}} \left\| {{\mathit{\boldsymbol{E}}_h} - {\mathit{\boldsymbol{d}}_h}\mathit{\boldsymbol{x}}_h^{\rm{T}}} \right\|_{\rm{F}}^2\;\;\;满足\;\;{\left\| {{\mathit{\boldsymbol{d}}_h}} \right\|_2} = 1 \end{array} $

(8)

式中：d_h为当前字典D的第h列向量；x_h^T为稀疏系数矩阵X中的第h行向量；E_h为残差矩阵。一般通过对残差矩阵E_h进行SVD更新d_h和x_h^T，但残差矩阵E_h的SVD往往需要巨大的运算成本，并不适应于三维数据的处理。因此，本文使用交替优化策略^[28]逼近SVD，以实现式(8)的快速求解。具体地，令函数$g\left( {{\mathit{\boldsymbol{d}}_h}, \mathit{\boldsymbol{x}}_h^{\rm{T}}} \right) = \left\| {{\mathit{\boldsymbol{E}}_h} - {\mathit{\boldsymbol{d}}_h}\mathit{\boldsymbol{x}}_h^{\rm{T}}} \right\|_{\rm{F}}^2$对d_h和x_h^T的偏导数都为零，可得

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{d}}_h} = \frac{{{\mathit{\boldsymbol{E}}_h} - \mathit{\boldsymbol{x}}_h^{\rm{T}}}}{{{{\left\| {{\mathit{\boldsymbol{E}}_h} - \mathit{\boldsymbol{x}}_h^{\rm{T}}} \right\|}_2}}}\\ \mathit{\boldsymbol{x}}_h^{\rm{T}} = \mathit{\boldsymbol{E}}_h^{\rm{T}}{\mathit{\boldsymbol{d}}_h} \end{array} \right. $

(9)

仅需一次迭代式(9)便可收敛到最优，极大地提高了字典更新的效率。

MCP稀疏分解算法是在K均值聚类对字典进行结构化处理(图 1)的基础上，利用多簇原子选择机制，使每次迭代都能实现多个最优原子的优选，提高了稀疏编码的效率。MCP稀疏分解算法描述如下。

(1) 输入训练数据集Y、初始字典D、聚类数k、最大迭代次数S和阈值ε。

(2) 初始化：稀疏解X=0；残差E₀=Y；设置每类原子索引集U_j、中心原子索引集F₀、多原子索引集J₀、所选原子索引集I₀均为空集；迭代次数s=1。

(3) 应用K均值聚类方法对字典进行结构化处理，其具体步骤为：

步骤一：将字典D划分为k组，并计算子字典D_Uj(j=1, …, k)的左奇异值向量来设置中心原子c_j^*；步骤二：重新计算所属类，即计算原子d_i的新所属类指标j^*(d_j)=arg maxj|〈c_j^*, d_i〉|²；并更新原子索引集U_j={i:j^*(d_i)=j}；

步骤三：计算新类D_Uj(j=1, 2，…, k)的左奇异值向量来更新中心原子c_j^*；

步骤四：重复步骤二和步骤三，直至D_Uj(j=1, 2，…, k)不再变化。

(4) 采用多原子选择机制优选多个最优原子，其具体步骤为：

步骤一：计算残差E_s和中心原子c_j^*的内积P_j=|〈E_s, c_j^*〉|，并通过阈值λ_s，选择多个原子${\mathit{\boldsymbol{F}}_s} = \left\{ {j:\left| {{\mathit{\boldsymbol{P}}_j}} \right| > {\lambda _s}} \right\}, {\lambda _s} = \mu {\lambda _{s - 1}}$，μ为阈值下降步长；

步骤二：在每个所选组中搜索最优原子J_s= $\left\{ {i:\arg {\max\limits_{i \in {{\bf{U}}_l}}}\left| {\left\langle {{\mathit{\boldsymbol{E}}_s}, {\mathit{\boldsymbol{d}}_i}} \right\rangle } \right|, l \in {\mathit{\boldsymbol{F}}_s}} \right\}$，并将新搜索的原子集J_s与前一步所选原子集I_s-1合并，即I_s=I_s-1∪J_s；

步骤三：将信号投影到子字典D_Is更新系数${{\bf{x}}_s} = {\left( {{\bf{D}}_{{\mathit{\boldsymbol{I}}_s}}^ \top {{\bf{D}}_{{{\bf{I}}_s}}}} \right)^{ - 1}}{\bf{D}}_{{I_s}}^ \top {\bf{y}}$，同时更新残差E_s=Y-Dx_s；

步骤四：判断是否停止迭代。若|E_s|²＜ε或s＞S，则停止迭代，否则令s=s+1重复步骤一~步骤三，开始下次迭代。

(5) 输出稀疏系数X和最优原子集合I。

1.3 训练集来源

数据分块是数据处理中的常用手段^[29]，即将整个数据分成大小相等的数据子集后再进行处理，在字典学习中其往往是训练集的样本来源。图 2展示的是三维数据向量化的过程，需要将三维数据分成大小一致的小正方体后向量化，依次存储到新矩阵中。实际处理中往往希望得到更多数据子集，故要求数据子集之间的重叠度最大，并且在数据边界处进行补零延拓。

1.4 初始字典

离散余弦变换(Discrete Cosine Transform，DCT)因其构建方式简单、稀疏表示性强被广泛应用于信号处理领域。本文以三维DCT字典作为初始字典对三维地震数据进行重建，DCT字典可由DCT获得，三维DCT的定义^[30]为

$ \begin{array}{l} F\left( {u,v,w} \right) = H\left( {u,v,w} \right)\frac{1}{{\sqrt {l \times m \times n} }} \times \\ \;\;\;\;\;\;\;\sum\limits_{i = 0}^{l - 1} {\sum\limits_{j = 0}^{m - 1} {\sum\limits_{z = 0}^{n - 1} {y\left( {i,j,z} \right)G\left( {u,v,w,i,j,z} \right)} } } \end{array} $

(10)

$ H\left( {u,v,w} \right) = c\left( u \right)c\left( v \right)c\left( w \right) $

(11)

$ c\left( u \right) = \left\{ {\begin{array}{*{20}{l}} 1&{u = 0}\\ {\sqrt 2 }&{u \ne 0} \end{array}} \right. $

(12)

$ \begin{array}{*{20}{c}} {G\left( {u,v,w,i,j,z} \right) = \cos \frac{{\left( {2i + 1} \right)u{\rm{ \mathsf{ π} }}}}{{2l}} \times }\\ {\cos \frac{{\left( {2j + 1} \right)v{\rm{ \mathsf{ π} }}}}{{2m}}\cos \frac{{\left( {2z + 1} \right)w{\rm{ \mathsf{ π} }}}}{{2n}}} \end{array} $

(13)

式中：y(i, j, z)为三维地震数据；i、j、z为三维数据三个方向的序号；l、m、n为三个方向的样点数。

1.5 数据重建

快速结构字典学习三维地震数据重建方法是在压缩感知技术框架下，首先对观测数据进行向量化处理以构造出用于字典学习的训练集，之后利用快速结构字典学习算法对训练集进行训练，并将训练产生的自适应字典作为地震数据稀疏表示的稀疏基，然后根据观测数据缺失道分布情况构建观测矩阵，并将其与自适应字典及正则化正交匹配追踪算法结合对目标函数进行稀疏促进求解，最后将求解结果进行逆向量化处理便可实现缺失数据的恢复重建，具体的重建流程如图 3所示。

2 数据处理效果分析

为了对重建结果进行定量评价，引入信噪比

$ {\rm{SNR}} = 10 \times \lg \frac{{\left\| \mathit{\boldsymbol{f}} \right\|_{\rm{F}}^2}}{{\left\| {\mathit{\boldsymbol{f}} - \mathit{\boldsymbol{\tilde f}}} \right\|_{\rm{F}}^2}} $

(14)

式中：f为原始地震数据；$\mathit{\boldsymbol{\tilde f}}$为重建后的地震数据。

2.1 模型数据测试

用图 4a所示SEG盐丘模型的炮集数据(大小为471×201×201)测试本文三维数据重建方法的可行性。图 4b是图 4a数据随机缺失70%的结果。图 4c是使用本方法的重建结果(SNR=21.3dB，计算耗时为9104.4s)，可见缺失的地震道得到了很好的恢复，重建结果也与原始数据比较吻合。图 4d为重建结果与原始数据的误差，可见重建误差很小，且基本不包含有效信号，说明本方法可以对三维缺失数据实现有效重建。图 5为通过快速结构字典学习算法训练用于重建图 4c的部分三维自适应字典。

2.2 实际数据重建

为说明本文方法相较于二维字典学习逐切片重建方式具有更高重建精度的优点，将两种方法应用于实际地震资料(数据量为320×160×160)的重建。图 6a是实际地震资料的切片显示，图 6b是图 6a数据随机缺失50%的结果。图 6c、图 6d分别为二维字典学习逐Crossline切片重建的结果及其与原始数据的残差。图 6e、图 6f分别是本文方法重建结果及残差。对比两种方法的重建结果及残差可以看出，两者虽都能实现对缺失数据的重建，但二维字典学习逐切片重建方法的残差相对较大，原因在于其仅从Crossline单个方向上进行重建，未考虑Inline方向及时间切片上的相互关系。而本文方法考虑到了切片之间的关联性，重建结果基本与原始数据无异，重建残差很小。

表 1为对图 6两种方法重建结果的信噪比和计算用时对比。与图 6结论一致，本文方法重建结果的SNR要高于二维字典学习逐切片重建结果。但从耗时方面分析，本文方法耗时要长于二维字典学习逐切片重建方法，原因在于本文方法需要训练三维自适应字典，大量的训练样本导致了训练耗时增加。

图 7为图 6数据Inline50剖面的F-K谱对比，可见，随机欠采样会引起频谱泄露，造成随机缺失数据全频带出现低幅值随机噪声(图 7b)。对比二维字典学习逐切片重建数据(图 7c)、本文方法重建数据(图 7e)和原始数据的F-K谱(图 7a)，可见两种方法重建结果的F-K谱均与原始数据的F-K谱差别不大，说明两种方法都可恢复缺失数据。对比图 7d与图 7f可见，两种重建方法均能消除随机欠采样所产生的全频带随机噪声，但就消除效果而言，图 7d消除的噪声幅值相较于图 7f较大，且主要分布在有效信号频谱附近，这表明二维字典学习重建残差中包含更多的有效信号。

上述结果表明，本文方法可以实现三维地震数据的高精度重建，且具有保幅性良好的特点。

3 结束语

本文提出一种快速结构字典学习三维数据重建方法。该方法在压缩感知理论框架下，利用快速结构字典学习算法对训练集进行训练，产生三维自适应字典，之后使用三维自适应字典、观测矩阵以及正则化正交匹配追踪算法对数据进行高精度重建。理论模型和实际资料的重建结果表明：相比二维字典逐切片重建方法，三维字典由于利用了地震数据三个维度上的完整结构和相互关系，能够更好地提高重建精度和保持重建结果的振幅变化特性。需要指出的是，本文虽在训练字典方面采用了快速求解算法，但需要训练三维自适应字典，运算效率仍低于二维字典学习逐切片重建方法，如何进一步提高效率将是下一步的研究重点。