基于迭代收缩阈值网络的地震数据重构研究

范 帅，邢 磊,2，李倩倩,2

(1.中国海洋大学 海底科学与探测技术教育部重点实验室，山东 青岛 266100；2.海洋国家实验室 海洋矿产资源评价与探测技术功能实验室，山东 青岛 266071)

1 引 言

地震数据重建是勘探地震学中一个重要的研究方向，近几十年来学者们已经提出了多种地震数据重构的方法[1,2]，如插值法。地震数据插值，也可以称为地震数据恢复，是地震资料处理中极其重要的一环。传统的地震数据插值方法可大致包含三大类：基于预测滤波器的插值算法，基于波动方程的插值算法以及基于信号分析处理的插值算法。Donoho等人(2006年)[3]提出压缩感知理论，即信号的采样频率可以低于Nyquist频率，并在采样后通过合适的恢复算法能较为完整的恢复原信号。压缩感知理论可以对数据同时进行采样和压缩，这样可以去除大量冗余数据并且释放有限的储存资源，压缩感知理论的提出推动了地震数据重构方法的发展。

Wang等人(2011年)[4]，将地震波场重构问题看作一个压缩感知问题，提出一种基于小波变换的分段随机子采样方法，通过求解L1范数来解决压缩感知问题，提高了波场恢复的质量。曹静杰等人(2012年)[5]在Curvelet域利用梯度投影的算法求解L0范数约束下的压缩感知问题。Gou等人(2014年)[6]通过Seislet变换对复杂地震数据进行稀疏后，应用Bregman迭代算法求解压缩感知问题中的混合范数优化问题。作为一种多尺度稀疏变换，Seislet变换可以将地震信号进行有效的压缩，因此对变换域的选择也尤为重要。随着字典学习的火热，周亚同等人(2014年)[7]将K-SVD字典学习方法应用于地震数据重建上，这个方法的应用是在压缩感知框架下的。该方法通过用K-SVD字典训练大量的地震样本数据，并在训练完成后得到超完备字典，并采用正则化正交匹配追踪(ROMP,Regularization Orthogonal Matching Pursuit)来对地震数据实现重建，可以实现地震数据的自适应稀疏扩展。但基于压缩感知的地震数据重构方法在不同变换域中的重构效果有很大区别，重构效果不稳定。

近几年深度学习在图像处理和图像修复领域体现出了非常好的效果， Cui等人(2014年)[8]通过对图像使用降采样后，将高分辨率图像降采样成低分辨率图像，将降采样后的低分辨率图像作为训练数据集输入进深度卷积神经网络中进行训练，之后使用训练过的网络对图像进行重建。Wang等人(2016年)[9]基于深度学习在图像处理方面的优势，使用高质量的核磁共振图像对神经网络进行训练来加速磁共振成象(MRI,Magnetic Resonance Imaging)的方法。同年，Sun等人[10]在压缩感知可以快速磁共振成象的基础上，通过将ADMM(ADMM,Alternating Direction Method of Multipliers)算法与深度神经网络结合，将网络的结构变为ADMM算法的展开形式来进行磁共振成象的压缩感知重建。Zhang等人(2018年)[11]将压缩感知凸优化算法中的迭代软阈值算法的每一步迭代过程用深度神经网络的方法展开，在自然图像重建和磁共振成像(MRI)中获得了良好的效果。Mandelli等人(2018)[12]构建了一种卷积自动编码器的特定类型的卷积神经网络来对损坏的地震数据进行插值和去噪。

受到图像处理和修复的启发，以及传统的地震数据重构方法存在的运算量大、重构成本大的问题，本文通过将迭代收缩软阈值算法和深度神经网络结合起来，将其应用在地震数据的重构方法中，通过使用大数据集对网络进行训练，实现端到端的映射，网络中的参数都可以通过数据集自动学习并更新，无需人为手动设置每步参数。实验结果证明，该方法可以高效、精确地恢复地震信号。

2 迭代收缩阈值算法

迭代收缩阈值算法(ISTA,Iterative Shrinkage-Thresholding Algorithm)最早可以追溯到2004年，Daubechies等人(2004年)[13]提出了一种迭代算法，迭代收缩软阈值算法由于其简单性、迭代速度快、迭代所需的算法复杂度小，因此被广泛的研究。

首先，地震数据重构问题可以同构构建一个线性方程组解决

y=Φx

(1)

其中，x为地震信号；Φ为观测矩阵；y为观测信号。通过使用最小二乘法可以求解上述问题：

(2)

在实际情况下，x一般都是不稀疏的，无论x是自然信号或是地震信号，它需要通过某种变换实现稀疏，例如稀疏基，x可以通过稀疏基来实现稀疏，稀疏基矩阵用Ψ来表示，s为稀疏变换后的稀疏系数。因此，公式(1)也就变成了：

y=ΦΨs

(3)

在式(3)中，y和Φ以及Ψ都是已知的值，因此只要求解出s就可以实现对地震数据的重建。一般而言，地震信号都是在时空域中进行表现，但由于Φ和Ψ都是已知的，因此引入感知矩阵的概念来简化公式，感知矩阵A就等于A=ΦΨ。式中观测矩阵Φ∈RM×N，Ψ∈RN×N，其中M≪N，即将稀疏信号(K-Sparse)从N维空间通过线性投影到M维空间当中。

y=Φx=ΦΨs=As

(4)

求解上述问题的方法是通过对f(x)=‖Ax-b‖2求导，得到求导之后的函数：

f′(x)=2AT(Ax-b)

(5)

要获得x的值只要让f′(x)=0就可以求得f(x)的最小值，前提是当f(x)为凸函数，且它的局部最小值就是它的最小值时。如果矩阵A为非奇异矩阵，即矩阵的行列式不为0，有逆矩阵时，该问题有精确解，x=A-1b。若当矩阵A为奇异矩阵时，由于A没有逆矩阵，导致方程等于0时无法求解，方程没有精确解。因此这时需要引入惩罚项来协助进行求解。

由于无法得到精确解，因此引入惩罚项的概念，使得x可以尽量稀疏：

min‖x‖1subject to‖Ax-b‖2≤ε

(6)

minx{F(x)≡‖Ax-b‖2+λ‖x‖1}

(7)

由于式(7)中函数并不是标准的凸函数，因此需要对其进行转换，将式(7)转换为

=f(α)+λg(α)

(8)

这个L1范数的正则化问题通常非常难解决，通过使用迭代收缩阈值算法可以解决上述问题，具体流程如图1所示。

图1 迭代收缩阈值算法流程图Fig.1 Iterative shrinking threshold algorithm flowchart

3 迭代收缩阈值网络

拟解决以下正则化问题：

(9)

式中,ψx表示对x进行某种稀疏变换后的变换系数;λ为正则化参数，一般来说是预定义的。根据迭代收缩阈值算法，式(9)可以转换为式(10)和式(11)。解决式(10)问题的方法可以是在对式(11)进行迭代：

r(k)=x(k-1)-ρΦT(Φx(k-1)-y)

(10)

(11)

在式(10)和式(11)中，k表示ISTA迭代指数;ρ是步长。式(10)是琐碎的，得到的结果是输入数据，而式(11)实际上是所谓近端映射的特例，即proxλφ(r(k))，此时φ(x)=‖Ψx‖1。形式上，由proxλφ(r)表示的正则化φ的近映射定义为：

(12)

由于地震数据往往是非线性的，无法通过简单的线性变换对其进行处理。因此通过设置卷积层以及应用ReLU(整流线性单元)函数来实现非线性功能，将两个卷积层与ReLU函数连接起来，实现对地震数据的非线性处理。

每层卷积层的卷积核不相同，如第一个卷积层(无偏置项)的卷积核的设为符合正态分布的(1×3×3)的三阶张量，共有32个卷积核；第二个卷积核设为符合正态分布的(32×3×3)的三阶张量，同样具有32个，这一步骤对应着对输入数据的稀疏变换。每一层步骤都对应着ISTA算法的一次迭代，可以看到F(x)=BReLU(Ax)。因此，将F(x)带入到式(12)中，用F(x)替换Ψx，可以得到以下式子：

(13)

利用ISTA算法很容易就能得到式(14)的解：

(14)

图2 F(·)函数结构Fig.2F(·) function structure

(15)

4 数值实验

4.1 模拟数据

使用规则采样以及随机采样对模拟数据进行采样。模拟数据大小总共为1 101×20 000，也就是将100炮数据沿x轴排列，每炮200道接收，总共的道数为20 000道，时间为1 101 ms。在这将100炮数据体重新抽样排列，每炮抽取100道以及100道的前800 ms传播时间作为输入数据，并将这100炮按照放炮次序进行排列，得到100×800×100的模拟数据体。

本次实验使用的地质模型为12层地层形成的地质体，震源为50 Hz的Ricker子波，模拟计算地震波在12层地质模型中得到的地震波波场，计算方法为声波方程。由于深度神经网络需要大量数据集，因此100炮地震数据显然不足以得到一个优秀的训练模型，数据量过小容易造成过拟合。因此，在这里将100×800×100的数据体中随机选取95炮，数据体大小为95×800×100，使用随机测量矩阵对这95炮数据进行数据集的扩充，即对这95炮数据体中每一炮都进行10次随机的采样，使得95炮数据体变为950×800×100的数据体，将该数据体作为ISTA-Net中的训练集；同理，在抽取95炮数据体后，遗留的5炮可以作为测试集，在950×800×100的地震数据体训练完毕后，将剩余5炮的数据导入，作为测试集进行测试，并将最终重构的结果进行保存。

将数据导入ISTA-Net中，预设总共的对整个ISTA-Net进行200次的迭代运算，在每次ISTA-Net中设置6个阶段层，每一批次导入20个训练块，总共950炮数据，因此设置950个训练块。在每一个训练周期完成后，使用学习率为0.000 1 的Adam算法对ISTA-Net的参数进行优化，训练200次后终止。图3表示的是模拟地震记录第5炮，其色标柱是振幅值。图3～图14的色标柱相同。图4表示的为对其进行随机采样50 % 后得到的采样矩阵。

图3 模拟地震记录第4炮数据Fig.3 The 4th shot data of simulated seismic records

图4 随机采样50 %的第4炮数据Fig.4 Randomly sample 50 % of the 4th shot data

此次模拟数据的实验的结果使用信噪比(SNR，Signal to Noise Ratio)进行衡量，衡量标准SNR的公式为

(16)

在式(15)中，参数p代表最终训练得到的重构结果，r为输入的原始地震数据，通过计算SNR可以得到重构效果的优劣，SNR越大，则表示重构后的地震数据的质量越高，效果越好。

图5是在使用迭代收缩阈值算法在使用Curvelet变换后的图像，图6表示的是Curvelet重构数据误差，Curvelet重构数据的SNR=12.691 2。从图5中可以看出，基于Curvelet变换后的重构数据同相轴较为连续， SNR较高，但重构误差相对较大，且在多道有连续缺失的位置，重构的效果并不好。

图5 Curvelet重构数据Fig.5 Curvelet reconstruction data

图6 Curvelet重构误差Fig.6 Curvelet reconstruction error

图7中表示的是利用迭代收缩阈值网络方法重构的地震数据，图8中是利用该方法的重构误差，使用迭代收缩阈值网络方法重构的数据的SNR=20.898 1。从图中可以明显看出，重构后的地震数据地震同相轴连续性好，无论是直达波还是反射波，其重构效果较好且绝对误差较小，连续缺失地震道的位置重构效果依然很好。

图7 迭代收缩阈值网络重构数据Fig.7 Reconstruction data of iterative shrinkage threshold network

图8 迭代收缩阈值网络重构误差Fig.8 Reconstruction error of iterative shrinkage threshold network

4.2 实际数据

本次实验使用的数据来源于神狐海域的一条拖缆测线的数据，拖缆数据中存在着多种噪声，在对数据进行采样之前，通过时频分析发现数据中存在低频干扰，因此在采样前对数据使用了一个带通滤波将低频干扰滤掉。整个数据体大小超过5个G，因此效仿模拟数据的截取，将拖缆数据截取部分用于此次实验。拖缆数据单炮由360道接收，将其中100炮的数据提取出来，并且每一炮中提取51道到150道的数据，时间方向上选择1 500 ms到2 300 ms处地震同相轴较为明显的部分，将其截取并导出，数据体大小为100×800×100。其中第4炮的地震数据如图9所示，图10为随机采样50 %后的实际拖缆数据。数据产生流程与模拟数据产生流程相似，随机抽取100炮中的95炮数据，对它进行重采样10倍后，可以得到950×800×100的数据体，这个数据体作为训练集进行训练，而剩余的5炮作为测试集，测试集数据体大小为5×800×100。

图9 实际拖缆数据Fig.9 Actual streamer data

图10 随机采样50 %的拖缆数据Fig.10 Randomly sample 50 % of the streamer data

图11是在迭代收缩阈值算法使用Curvelet变换后的重构数据，图12为使用Curvelet变换后的重构误差，重构后的数据SNR=5.286 0。从图11中可以看出，基于Curvelet变换后的实际重构数据同相轴不连续，SNR较低，且在多道有连续缺失的位置以及噪声较大的位置，重构的效果并不好。

图11 Curvelet重构数据Fig.11 Curvelet reconstruction data

图12 Curvelet重构误差Fig.12 Curvelet reconstruction error

图13中表示的是使用迭代收缩阈值方法重构的实际拖缆数据，图14表示的是利用该方法重构的误差，重构后的数据的SNR=11.073 6。从图中可以明显看出，重构后的地震数据的地震同相轴连续性较好，数据上半部分直达波的恢复效果较好，绝对误差相对较小，连续缺失地震道的位置重构效果较好，重构时去除了一部分噪声，突出了有效同相轴。

图13 迭代收缩阈值网络重构数据Fig.13 Reconstruction data of iterative shrinkage threshold network

图14 迭代收缩阈值网络重构误差Fig.14 Reconstruction error of iterative shrinkage threshold network

5 结 论

基于复杂地质条件下的地震数据缺失的问题，结合传统的凸优化算法以及近年来在图像处理效果较好的深度神经网络方法得到迭代收缩阈值方法，并将其应用在模拟数据和实际拖缆数据的地震数据重构上，通过大量的实验，主要得到以下几点结论：

1)通过将迭代收缩阈值算法与深度神经网络相结合，分析发现迭代收缩阈值算法的每一步迭代过程可以用深度神经网络的每一层进行替换，因此可以将迭代收缩阈值算法映射到神经网络中去，结合深度神经网络的特点，克服了迭代阈值收缩算子需要预设步长、迭代次数以及变换域不符的问题，通过将迭代收缩阈值算法映射到深度神经网络中，各个参数可以实现端到端的学习，不需要对各个参数进行预设，其可以通过神经网络进行自我更新学习，并将稀疏变换矩阵用稀疏变换算子代替，采用非线性可学习的稀疏变换，可以对数据进行更好的稀疏操作。

2)通过将迭代收缩阈值网络重构方法应用于模拟数据和实际拖缆数据上，发现该方法在各个方面的重构效果都要好于传统的迭代收缩阈值算法的表现，尤其是该方法在随机采样时对连续缺失道的恢复效果远好于其他方法。在实际地震数据中迭代收缩阈值网络重构方法可以用于对含噪数据进行一些去噪操作，使得同相轴可以更好的凸显出来。该方法的重构后的数据分辨率较好，SNR值较高，且同相轴连续。