0 引言

在地震勘探中，受经济因素或自然条件的限制，现场采集的地震信号会出现空道或坏道，造成空间假频、成像分辨率低等问题，进而影响后续的地震资料处理和解释。因此，将低空间采样率的地震数据以较高精度恢复成高空间采样率数据成为一项重要课题，众多学者陆续提出一系列地震数据重建法，可分为模型驱动和数据驱动两类。

基于模型驱动的地震数据重建法主要包括：波动方程重建法^[1-3]、稀疏变换重建法^[4-6]和预测滤波重建法^[7-8]及低秩矩阵完备重建法^[9-11]。这类方法往往依靠地震数据的先验信息构建模型，因此重建精度也会受其影响。如波动方程法的重建精度会十分依赖于所获取的地下介质速度信息。同时，此类方法所构建模型会随数据特征而改变，欠普适性。

随着人工智能机器学习和深度学习的飞速发展，基于数据驱动的地震数据重建法也随之兴起。相较于模型驱动重建法，数据驱动方法通过对地震数据自身的学习获取特征分布规律从而实现数据重建。典型的人工智能机器学习类方法包括字典学习重建算法^[12-13]及以支持向量回归(Support Vector Regression，SVR)算法^[14]为代表的浅层机器学习算法。人工智能深度学习类方法则是通过多层神经网络提取地震数据特征进而完成重建，包括深度卷积网络^[15]、残差网络^[16]及U-net网络^[17]等方法。

基于数据驱动重建法的优势在于不需传统方法中包括稀疏性、低秩性等先验假设，对地震数据具有普适性。但深度学习重建法对数据量的要求极高，训练时间过长，且对硬件配置有较高要求，有一定的局限性；相对而言，机器学习重建法的使用条件更灵活，但在特征提取方面仍有待提高。综合考虑这些因素，本文选择机器学习中具有良好泛化能力和鲁棒性的SVR算法，同时引入Gabor变换提取地震数据的纹理特征以弥补机器学习自身特征提取能力的不足。

Gabor滤波器作为一种纹理特征提取方法，具有可提取目标的局部空间和频率域信息的良好特性，有助于地震数据物理信息的提取与利用。在地震勘探领域，Gabor变换主要用于地震数据处理^[18]，如利用Gabor算法重建地震数据^[19]、地震信号去噪^[20]及提高地震数据分辨率^[21]，充分表明Gabor变换能有效挖掘地震数据的物理信息。

综上所述，本文引入Gabor变换作为纹理特征提取方法，提取地震数据的物理信息，并加入SVR算法的特征向量中，期待增强机器学习方法的特征提取能力，从而提高经训练得到的回归模型的预测精度。

1 支持向量回归重建算法

SVR算法是支持向量机为解决回归问题而衍生出的机器学习方法。在SVR中，假定一个训练样本

$ \left\{\boldsymbol{x}_i, y_i\right\}_i^N $

(1)

式中：x_i∈R^n_d为输入向量, n_d为维数；y_i∈R为输出标签；N为样本总数。构造最优线性函数

$ f(\boldsymbol{x})=\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x}+b $

(2)

式中：w为权重；b为偏置。为了最小化‖w‖²，引入不敏感损失函数ε，可得

$ \begin{aligned} & \min _{\boldsymbol{w}, b} \frac{1}{2}\|\boldsymbol{w}\|^2+C \sum\limits_{i=1}^N l_{\varepsilon}\left[f\left(\boldsymbol{x}_i\right)-y_i\right] \\ & \;\;\quad\text { s. t. } \quad\quad\left|y_i-f(\boldsymbol{x})\right| \leqslant \varepsilon \\ & \end{aligned} $

(3)

式中：C为正则化常数；l_ε为ε-不敏感损失函数。

为处理函数在ε精度不能估计的样本，引入非负的松弛变量ξ和ξ^*，式(3)可改写成

$ \min _\limits{w, b} \frac{1}{2}\|\boldsymbol{w}\|^2+C \sum\limits_{i=1}^N\left(\xi_i+\xi_i^*\right) $

(4)

引入拉格朗日函数(对应算子为α_i、α_i^*)和KKT(Karush-Kuhn-Tucke)最优化条件，得到回归函数

$ f(\boldsymbol{x})=\sum\limits_{i=1}^N\left(\alpha_i^*-\alpha_i\right) \boldsymbol{x}_i^{\mathrm{T}} \boldsymbol{x}+b $

(5)

对于非线性回归，基本思想是通过非线性变换ϕ(·)将原始数据映射到足够高维特征空间，然后在高维空间中基于线性回归解决低维空间的非线性回归问题。核函数K(x, x′)可将该非线性映射转换为隐性映射，使运算更简捷、灵活。常用的有多项式核函数、高斯径向(RBF)核函数和Sigmoid核函数，以及根据Ricker子波推导得到的更适用于地震数据的Ricker子波核函数^[22-23]。核函数的选择会对SVR的性能产生极大影响。

SVR地震数据重建算法^[24]的基本思想是通过预处理将矩阵形式的地震数据转化为SVR所需的向量点对格式，即特征向量—标签形式，经过训练得到回归模型，进而预测缺失的地震数据。该算法的特征向量仅包含原始地震数据，未考虑加入更多有效的地震数据物理信息，因此本文将引入特征提取方法以提高重建精度。

2 Gabor变换

选取Gabor滤波器提取地震数据中的纹理特征，作为地震数据的物理信息。Gabor滤波器的本质是在短时傅里叶变换中加入高斯函数作为窗函数，通过窗函数对信号进行时频分析，故具有良好的局部化特性，调整Gabor滤波器的方向和尺度更灵活方便。利用Gabor滤波器提取特征的步骤为：首先设计滤波器的参数和滤波器组的布局；然后从输出结果中提取有效特征。实际中常用的是二维Gabor滤波器，也称空间滤波器，是一个高斯函数与正弦平面波的乘积，可在给定区域提取频域特征，若像素坐标设为(x，y)，其表达式为

$ \begin{aligned} G(x, y ; \lambda, \theta, \psi, \sigma, \gamma)= & \exp \left(-\frac{x^{\prime 2}+\gamma^2 y^{\prime 2}}{2 \sigma^2}\right) \times \\ & \exp \left[\mathrm{i}\left(2 \pi \frac{x^{\prime}}{\lambda}+\psi\right)\right] \end{aligned} $

(6)

实部和虚部分别为

$ \begin{aligned} g(x, y ; \lambda, \theta, \psi, \sigma, \gamma)= & \exp \left(-\frac{x^{\prime 2}+\gamma^2 y^{\prime 2}}{2 \sigma^2}\right) \times \\ & \cos \left(2 \pi \frac{x^{\prime}}{\lambda}+\psi\right) \end{aligned} $

(7)

$ \begin{aligned} g(x, y ; \lambda, \theta, \psi, \sigma, \gamma)= & \exp \left(-\frac{x^{\prime 2}+\gamma^2 y^{\prime 2}}{2 \sigma^2}\right) \times \\ & \sin \left(2 \pi \frac{x^{\prime}}{\lambda}+\psi\right) \end{aligned} $

(8)

其中

$ \begin{gathered} x^{\prime}=x \cos \theta+y \sin \theta \\ y^{\prime}=-x \sin \theta+y \cos \theta \end{gathered} $

(9)

式中：λ为Gabor滤波器波长；θ为核函数方向；ψ为相位偏移；σ为高斯标准差；γ为空间纵横比。

二维Gabor滤波器的实部可进行平滑滤波，虚部可用作边缘检测。一个经典的空间Gabor函数h(x，y)及其傅里叶变换H(u，v)有以下形式

$ \left\{\begin{array}{l} g(x, y)=\frac{1}{2 \pi \sigma_x \sigma_y} \exp \left[-\frac{1}{2}\left(\frac{x^2}{\sigma_x^2}+\frac{y^2}{\sigma_y^2}\right)\right] \\ h(x, y)=g(x, y) \exp (2 \pi \mathrm{j} W x) \\ H(u, v)=\exp \left\{-\frac{1}{2}\left[\frac{(u-W)^2}{\sigma_u^2}+\frac{v^2}{\sigma_v^2}\right]\right\} \end{array}\right. $

(10)

式中：g(x，y)为高斯函数；σ_x和σ_y是高斯函数在两坐标轴上的标准差；W是沿横轴的复正弦函数的频率。

计算中需将Gabor函数分解成实部h_R(x，y)和虚部h_I(x，y)，进而对图像I进行滤波，即

$ S(x, y)=\sqrt{\left(h_{\mathrm{R}} * I\right)^2(x, y)+\left(h_{\mathrm{I}} * I\right)^2(x, y)} $

(11)

式中：(h*I)表示图像I与滤波器h的卷积；S(x, y)为提取的特征图像。以S(x，y)为母小波，对其参数进行调整，得到一系列自相似滤波器，称为Gabor小波

$ h_{m n}=a^{-m} h\left(x^{\prime}, y^{\prime}\right) \quad a>1 $

(12)

其中

$ \begin{gathered} x^{\prime}=a^{-m}(x \cos \theta+y \sin \theta) \\ y^{\prime}=a^{-m}(-x \cos \theta+y \sin \theta) \\ \theta=n \pi / K \end{gathered} $

(13)

式中：a^－m为尺度因子；m和n的取值范围分别为m=0, 1, …，S－1和n=0, 1, …，K－1，其中S和K为尺度和方向的数目。通过改变m和n，就能得到一组具有不同方向和尺度的Gabor滤波器，提取图像多个方向的空间频率(尺度)和局部性纹理特征。

如图 1所示，构建一个波长分别为3、7、15，方向为0°、45°、90°的Gabor滤波器组。

3 基于Gabor纹理学习的重建算法

基于Gabor变换的支持向量回归重建算法(Gabor SVR)主要分为模型训练阶段和预测重建阶段。其中模型训练分为预处理阶段、训练特征提取及向量对构建阶段、回归模型建立阶段。

预处理阶段：为了确保训练模型适用于预测阶段中的缺失地震数据，应对训练数据进行预处理。实验使用与缺失的地震数据相同的掩码矩阵对完整的训练数据进行下采样。然后，对下采样数据进行上采样或预插值。

训练特征提取及向量对构建阶段：首先提取训练特征，包括Gabor变换提取的纹理特征及原始数据特征，还要将其转换成符合SVR框架输入的向量对格式，即特征向量—标签形式。

原始数据特征的向量对构建的具体过程如下：在预插值后的地震数据(i，j)处，以该处为中心取大小为m×m的局部矩阵小块，并转换为长度为m²的行向量，即得到数据特征a_B；然后从地震数据中提取Gabor纹理特征。据上文所述，提取Gabor特征即是对地震数据进行实数部分的Gabor变换，首先需设计Gabor滤波器的参数，确定正弦函数的λ和θ两个参数，对于不同种类的地震数据，应选择适合的参数值，保证达到最佳重建效果。提取完成后，将所得数据分割成众多子数据块，目的是降低特征维度，最后统一转化成行向量(过程与原始数据特征转换步骤一致)，即可得到原始地震数据在该方向和尺度变换后的特征向量a_G。将a_B与a_G合并作为特征向量X，标签$\hat{Y}$的设计则是选取特征向量构建时居中的(i，j)处对应的像素值，整个流程如图 2所示。

回归模型建立阶段：将构建好的特征向量和标签输入SVR框架进行训练，得到非线性回归模型

$ f(\boldsymbol{x})=\sum\limits_{i=1}^N\left(\alpha_i-\alpha_i^*\right) K\left(\boldsymbol{x}_i, \boldsymbol{x}\right)+b $

(14)

用于预测缺失的地震数据。

预测重建阶段：对测试数据进行与训练数据同样的预处理；然后通过Gabor变换提取缺失数据的纹理特征，与缺失数据相结合，一并构成测试数据的特征向量；最后将这些特征向量X输入训练好的回归模型f(x)，从而预测缺失的像素值，完成地震数据的重建。

4 数值实验

为了验证本文所提Gabor SVR方法的有效性，分别测试了合成地震数据和两块不同大小的实际地震数据的重建效果，并与SVR算法对比。数据重建效果用信噪比(SNR)衡量，具体表达式为

$ \mathrm{SNR}=10 \lg \frac{\|I\|_{\mathrm{F}}^2}{\left\|I_{\mathrm{n}}-I\right\|_{\mathrm{F}}^2} $

(15)

式中：I_n代表重建数据；I为原始数据。

4.1 合成数据实验

首先选用大小为128×128、地震道空间采样率为25%的合成地震数据集验证重建效果。图 3a为完整合成地震数据集，图 3b为地震道空间采样率占全部的25%的数据集。训练集选择4张完整的地震数据集，将对其进行预处理，选取3×3的局部图像块及Gabor纹理特征提取后，经合并得到63504个训练点对，用于构建回归模型。

Gabor SVR重建结果(图 3c)的信噪比为42.29 dB；SVR算法重建结果(图 3d)的信噪比为41.52 dB，低于Gabor SVR方法的信噪比，证明Gabor变换提取的纹理信息能提高地震数据信息的利用率，进而提高算法的重建精度。

4.2 实际数据实验

针对大小为128×128的实际地震数据集(图 4a)，按地震道空间采样率为50%得到对应的缺失数据(图 4b)，以该实际地震数据测试算法的重建效果。选取8张128×128完整地震数据集作为训练集，将训练集进行预处理和Gabor变换提取特征后，可提取127008个向量对，经训练得到回归模型，进行测试集的预测，得到实际地震数据的重建结果。其中Gabor SVR方法重建结果(图 4c)的信噪比为35.00 dB，SVR方法重建结果(图 4d)的信噪比为32.20 dB。从实验数值看，Gabor SVR重建后的信噪比高于SVR重建后的信噪比，证明Gabor SVR针对实际地震数据具有更好的重建性能，且Gabor变换提取实际地震数据的纹理信息对提高SVR算法的重建质量效果显著。

为了更直观地对比两种方法的重建结果，分别用两种方法抽取一些地震道进行比较。可看出Gabor SVR方法重建的地震道波形更接近于原始地震道，表明重建效果更好(图 5)。实验还测试了仅有25%的地震道空间采样率的重建效果，Gabor SVR方法重建结果的信噪比为22.36 dB，SVR算法结果的信噪比为20.83 dB，表明Gabor SVR方法在地震道缺失极严重的情形下同样优于SVR方法。

第二组实际数据实验中，选用大小为512×128的实际数据(图 6a)，仍采用地震道空间采样率为50%的缺失数据(图 6b)。训练集选取6张地震数据进行训练，经过同样的步骤处理，可提取训练集中的385560个数据点对，作为特征向量训练回归模型。Gabor SVR方法重建结果(图 6c)的信噪比为59.65 dB；SVR方法重建结果(图 6d)的信噪比为59.14 dB，低于Gabor SVR方法的信噪比，表明Gabor SVR重建效果更好。

同时也测试了两种方法对该实际地震数据以地震道空间采样率为25%时的重建效果，Gabor SVR重建后的信噪比为46.34 dB，而SVR方法的信噪比为44.46 dB，证明针对此数据集Gabor SVR方法在地震道缺失多的情况下重建效果仍优于SVR方法。

4.3 算法讨论部分

SVR算法是基于数据驱动的学习类算法，将地震数据的物理信息与隐含在数据中的统计规律结合，是提高重建算法精度的重要途经。本文数值实验证明，Gabor SVR算法将地震数据的纹理特征融入数据特征向量，提高了物理信息的利用率，进而提高重建算法的精度。SVR算法中核函数是支持向量回归机特征提取的关键因素，与地震数据匹配的核函数更有助于物理特征的挖掘和学习，是提高重建精度的另一重要途径。Ricker子波核函数是由可模拟地震信号的Ricker子波函数经数学变形推导而来，即

$ K\left(\boldsymbol{x}, \boldsymbol{x}^{\prime}\right)=\left(N-2 \pi^2 f^2\left\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\|^2\right) \mathrm{e}^{-\pi^2 f^2\left\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\|^2} $

(16)

式中：N是输入空间特征向量x的数量；f为频率。采用该核函数可提高算法的重建精度^[22-23]。Gabor SVR算法采用RBF核函数，故本文提出Gabor Ricker SVR算法，即将Gabor SVR算法中的RBF核函数替换为Ricker子波核函数。

针对第一组实际数据重建问题，采用Gabor Ricker SVR算法重建结果(图 7)的信噪比为35.44 dB，略高于Gabor SVR方法。再次证明物理信息的充分挖掘和利用可提高数据驱动算法重建精度。

除重建精度外，重建效率也是重要的评价指标之一。基于数据驱动的学习类算法通常依赖大量的训练数据，但随着数据量的增加，计算成本也会随着增加。众所周知，训练数据多到一定程度，必然会在样本中包含冗余信息，且带来高昂的计算成本。从统计学角度出发，随机抽样是解决该问题的有效方法。由于随机性，在降低样本数量的同时，可包含绝大部分的特征信息。因此基于随机抽样选取不同数量的训练样本，测试三种方法(SVR、Gabor SVR、Gabor Ricker SVR)在50%地震道空间采样率时的地震数据重建效果。

表 1记录的是随机抽取不同百分比的训练数据时，三种方法重建结果的信噪比，图 8更清晰、细致展示了三种方法信噪比的变化过程。从纵向对比可见，无论是哪种训练数据随机抽样的百分比，Gabor SVR算法的信噪比均高于SVR算法，Gabor Ricker SVR又高于Gabor SVR(RBF)。此外，Gabor SVR算法和Gabor Ricker SVR算法采用5%训练数据时的重建信噪比均高于SVR算法采用100%训练数据时的信噪比。该结果充分说明：训练过程中物理信息利用率的提高，可不断提高算法的重建精度；也印证物理特征与数据特征的融合，是提高数据驱动类算法重建质量的有效途径之一。

表 2记录的是基于随机抽样训练数据，不同百分比训练数据量时，三种算法所需训练时间。首先纵向对比，在相同训练数据量下，由于Gabor SVR、Gabor Ricker SVR在训练过程中加入了物理特征，计算的复杂度增大。在100%(127008)训练样本下，SVR、Gabor SVR及Gabor Ricker SVR三种方法的计算时间分别为2908.44 s，4589.00 s和9067.42 s。但达到相同的重建精度，Gabor SVR和Gabor Ricker SVR算法仅需极少训练数据。以SVR算法的最高重建信噪比32.20 dB为例，耗时为2908.44 s，其他两种方法仅需5%的训练数据即可达到该信噪比，而耗时仅需约4 s。因此，证明物理信息的融入可降低对巨大数据量的依赖性，也充分印证了提高地震数据物理信息利用率的重要性。

综上所述，随机抽样部分训练样本可包含地震数据绝大部分特征，减少信息冗余，在信噪比损失不大的情况下提高重建效率，也证明Gabor变换提取的纹理信息，即地震数据的隐含特征更有助于回归模型的构建，从而提高精度，且在随机抽样少量训练数据的条件下表现更突出，既节省计算成本，缩短训练时间，还提高了重建质量。

5 结束语

本文提出一种基于Gabor纹理特征的数据驱动类地震数据重建方法。利用Gabor滤波器局部空间信息和频域信息提取能力获取地震数据的物理信息并加入到特征向量中，再通过SVR算法构建回归学习模型，用于缺失数据重建。大量的数值实验结果证明Gabor变换提取的地震数据物理信息有助于提高SVR方法的重建质量，具有更高重建信噪比。