基于时空约束压缩感知的地震数据重建①

石 敏 朱震东 路 昊 朱登明 周 军

(*华北电力大学控制与计算机工程学院 北京102206)

(**中国科学院计算所前瞻研究实验室 北京100190)

(***中国石油集团测井有限公司 西安710077)

0 引言

地震勘探是油气勘探的重要手段。随着油气勘探的不断深入,勘探目标区域的结构和周围环境越来越复杂,使得收集的数据不规则甚至不完整,从而对后续分析解释及油气判断产生了影响。同时,这也意味着更高的勘探成本。如果能够基于低密度采样数据重建出准确的高密度地震数据,便能够改善现有勘探方法,有效地利用获得的数据集可为估算地下石油储量的形成提供更可靠的支持,并从一定程度上降低勘探成本。地震数据在变换域中呈现稀疏分布,并且在时域和空域上也显示出很强的相关性,这为地震数据重建提供了可能。

在这项工作中,针对不完整地震数据的压缩感知重建方法,本文分析了地震勘探和测井勘探的数据特点,研究了一种时空约束模型,对传统压缩感知模型进行了改进。本文通过内核奇异值分解(kernel singular value decomposition,K-SVD)算法来训练超完备字典,并使用改进的稀疏自适应匹配追踪算法(sparsity of adaptive matching pursuit,SAMP)解决相应的优化问题,从而完成压缩感知重建。最终,进行了大量的对比实验,验证了本文算法的重建效果及效率。

总而言之,本文的贡献如下。

(1) 分析了地震数据的特征,改进了经典的压缩感知模型,添加时空相关信息作为压缩感知模型的约束。

(2) 改进了稀疏自适应匹配追踪算法,增加了初始稀疏性估计和可变步长的策略,确保了重建精度的同时,提高了算法效率。

(3) 在真实地震数据和微电阻率成像数据上实现了本文算法,证明了其出色的重建能力和泛化能力。

1 相关工作

基于地震数据的重建方法,通常被划分为以下3 种:第一种方法是基于预测滤波的,通常通过高斯窗口对不规则样本数据进行插值,这种方法导致很多错误出现[1-2]。第二种方法基于波动方程,通过正反演算子解决一个反问题。这种方法结合地下结构的先验信息,并依据波传播的物理特性实现重建,这种方法通常计算量很大[3]。第三种是基于变换的方法,首先对地震数据进行变换,进而在变换域中通过迭代求解等方法实现地震数据的重建。由于具备稳定性和可解释性,该方法得到了较为广泛的发展和应用[4-5]。Abma 等人[6]在地震重建领域引入了压缩感知算法[7-8]。Wen 等人[9]指出,在基于压缩感知方法的重建中,影响地震数据重建效果及效率的3 个主要因素是稀疏变换,迭代算法以及阈值模型。Bora 等人[10]提出了一种基于生成网络的压缩感知算法,该算法实现了基于更少的数据进行重建,但是带来了诸如训练困难和精度不足的问题。

稀疏变换对重建效果有很大的影响,其包括离散余弦变换(DCT),傅立叶变换(Fourier)和超完备字典[11]等。在压缩感知的应用领域中,DCT 和Fourier 较为常见,但这两种变换都难以有效地识别局部特征。短时傅立叶变换改善了这一问题,但是对于诸如地震数据之类的复杂信号,在不同时间的波形变化很大,短时傅立叶变换的时频局部化能力仍然有限。小波分析采用并完善了这种局部化思想[12],但是不能识别方向。直到后来,提出了曲波(Curvelet)变换[13-14],在保证多尺度识别能力情况下优化了多方向识别能力。因此在地震数据重建领域,曲波变换作为最佳稀疏变换方法之一而被经常应用[15]。剪切波变换[16]对多方向识别能力进一步发展,可以使地震信号的稀疏表示更加稀疏。

目前,基于上述变化的压缩感知理论已被用于地震数据重建[17-19],取得了良好的数据重建效果。但是,上述变换方法有一个共同点,变换基都是预先选择的,并不一定能够适合特定场景数据本身的特征。由于实际的地震数据一般非常复杂,通常会涉及多种类型的元素,预设的稀疏基很难对复杂的地震数据进行有效变换。在相关研究中,使用内核奇异值分解算法来学习超完备字典,能够结合场景数据,训练得到更适应数据特征的字典作为稀疏变换基[20-22]。在最近的研究中,基于字典学习的地震数据处理展示了比传统字典更强的稀疏重建能力[23-25]。

迭代算法不仅会影响到整体的重建效果,还极大程度上决定了算法的效率。由于早期迭代算法(如内点算法[26]和梯度投影法[27])的缺点,提出了重建贪婪算法,包括传统的贪婪算法,例如匹配追踪算法和正交匹配追踪算法(OMP)[28-29],以及基于传统算法改进得到的分段正交匹配追踪[30]和规则化正交匹配追踪[31]等,但是这些方法通常只应用于信号具有较低的稀疏度的场景。子空间追踪算法(SP)[32]和压缩采样匹配追踪算法(CoSaMP)[33]是两种具有良好性能的相似算法。Blumensath[34]提出了迭代硬阈值(IHT)及其改进算法,进一步提升了重建效果。但是,这些传统的贪婪算法都需要预先获得信号稀疏性。对此,Thong 等人[35]提出了一种自适应估计信号稀疏性的迭代算法,即稀疏自适应匹配追踪(SAMP)算法。解决了估计信号稀疏性的问题,但是并没有得到最佳的重建精度,并且算法效率较低。

综上所述,研究基于压缩感知的地震数据重建,主要在于研究如何基于地震数据的相关特征,构造稀疏效果更好的变换基以及该场景下更高效快速的迭代算法。另外,传统的压缩感知框架下的重建是针对单帧数据进行的,而地震数据在时域和空域上具有相关性,在重建地震数据时,仅考虑单帧数据的信息,而忽略了连续帧之间的相关信息,这将会从一定程度上影响地震数据的重建效果。因此本文也对此也进行了研究。

2 基于时空约束的压缩感知

2.1 地震数据的时空相关性

地震信号是复杂的,甚至瞬态特征也是不稳定的[36]。由于岩层的密度差异的存在,在不同岩层的交界处会存在地震波反射和地震波折射的现象,因此地震波数据与时间域相关,不同深度的地质信息可以由探测器在不同时间点获取地震波来表示。而震波在同一层的岩石中连续传递,因此地震波数据也与空间域有关,不同空间位置的地质信息可以由波器延测线的水平分布来表示。

目前,传统的压缩感知算法能够基于不完整的地震数据ym∈RM,通过重建得到完整的地震数据fm∈RN(M ＜N)。然而,传统的压缩感知算法仅满足了单帧数据重建的合理性,而忽略了地震数据相邻帧之间的相关性。如图1 所示,利用传统的压缩感知算法,基于单帧信息进行数据重建,重建结果中的高频区域具有较为明显的“失真”现象。因此,在利用压缩感知算法对地震数据进行重建时,有必要利用数据相邻帧之间的相关信息对重建进行优化或约束。

图1 单帧数据重建结果

2.2 基于时空约束的压缩感知算法

压缩感知算法使用观测矩阵用于描述采样数据和完整地震数据之间的关系。

其中,Φm∈RM×N是观测矩阵,而m是地震道号。利用稀疏变换基φ进行稀疏,完整的地震数据fm可以进而表示为

其中,xm是稀疏向量,将两式结合在一起,可以得到压缩感知的基本表示形式:

在传统的压缩感知算法基础上进行了修改,增加时空相关性约束,使其能够更好地重建连续多帧地震数据。设{fm-n,fm-n+1,fm-n+2,…,fm} 为重建的连续n+1 个数据帧,并定义R为连续帧之间变化的能量损失。

可以看出,R可以衡量帧之间的差异,并与相关性呈负相关。将式(2)带入式(6),可以得到:

目标是找到R的最小值,这等同于求解目标。对于稀疏数据,最小L1 范式和最小L0 范式从一定程度上是等效的。因此,可以将目标转换为最小L0 范式的解。

从而可以为式(4)添加时空相关信息R,使重建的数据具有尽可能小的变化损失:

其中,λ代表时空相关损失的比例。还可以在式(9)上进行等效转换,以使加入时空信息的模型与压缩感知的表示相同,便于使用迭代算法来重建数据。本文将λ设为1,然后进行等效转换得到:

3 地震数据重建

本节对地震数据重建算法进行实现。地震数据的特征比较复杂,不同地质环境下的地震数据也有很大差异。非自适应变换基不能很好地表示地震特征。因此,选择利用K-SVD 算法来训练超完备字典,以此用重建所需的稀疏变换基。本文还通过添加初始稀疏性估计和变步长更新策略来改进SAMP算法。基于第3 节中的优化模型,使用改进的SAMP 算法来求解稀疏矩阵。将稀疏矩阵和超完备字典相乘以获得重建数据。完整的流程图如图2 所示。

图2 地震数据重建流程图

3.1 基于K-SVD 的超完备词典学习

K-SVD 算法的本质是迭代思想,每次迭代使用范数稀疏约束跟踪计算稀疏系数,并利用奇异值分解算法更新字典原子,使得稀疏系数与字典能够得到同步更新,最终可以根据稀疏约束来自适应地训练出超完备字典。

假设,矩阵Y=是N个完整的地震数据切片的集合,代表了字典学习过程中的训练数据,矩阵X=是N个与训练数据Y=所对应的稀疏表示系数向量集合,矩阵φ∈Rn×K表示超完备字典。那么超完备字典训练的过程就可以表示为一个优化问题:

其中,T0是稀疏表示系数中非零元素的最大数量。超完备字典训练的整体步骤如下。

步骤1利用某种变换基进行字典初始化。

步骤2根据已知字典φ,使用迭代算法求解每个样本yi的稀疏系数向量xi,即

步骤3更新字典φ。如果满足收敛条件抑或达到预先设置的迭代次数上限,则获得最终字典φ,否则转向步骤2。设向量dk为要更新的字典φ的第k列原子,此时样本集的分解形式可表示为

其中,向量是与dk相对应的X中第k行向量;矩阵Ek代表提取dk之后的误差矩阵。

3.2 基于自适应动态步长SAMP 算法的重建

在地震数据重建领域通常情况下待重建数据的稀疏性是未知的,而在这种情况下,SAMP 算法可以通过设置固定步长s动态更新稀疏度来达到最优。通常设定稀疏度为最小值,并在每次迭代中,根据残值判断是否需要增加稀疏度,逐渐逼近信号实际的稀疏度,从而得到稀疏度的最佳估计,最终实现重建结果。

由于SAMP 算法中步长是固定的,如果将步长s设置太大,则信号的真实稀疏性可能会被跳过,得到的结果可能并非全局最优解,这将导致重建精度降低;如果步长s设置太小,迭代次数将会显著增加,严重影响了算法的效率。因此,本文对SAMP 算法进行了两方面的改进,即初始稀疏度K0估计和动态步长更新策略。

其中λ∈(0,1) 为步长变化率。本文实验设定为η=0.1,λ=0.5。改进后的SAMP 算法完整步骤如下。

输入传感矩阵θ,相关度n,观测向量{y1,y2,…,yn},迭代次数M,稀疏系数变化率阈值η,初始步长s,步长变化率λ,邻接矩阵A。

输出信号稀疏表示系数估计^x。

步骤1根据式(10)构造时空传感矩阵A,y=(y1,y2,…,yn)T。

步骤2初始化:g=ATy,F0=Ø,K0=1。

步骤3取g中K0个最大值的索引组成F0。

步骤4如果,则K0=K0+1,重复步骤4;否则,到步骤5。

步骤5初始化:^x=0,r0=y,I=K0,k=1。

4 实验结果和分析

4.1 量化指标

本文在真实地震数据进行量化实验,来验证本文算法的可行性和有效性。衡量数据重建效果的量化指标是信噪比(SNR)以及峰值信噪比(PSNR)。

其中,f是原始数据,是重建的数据,MSE是原始数据和重建数据的均方误差。

4.2 稀疏变换矩阵对地震数据重建的影响

本文中的数据集来自实际中收集的完整地震数据。为了进行实验分析和验证,将其设为601 ×626切片数据,用作原始数据集,如图3(左)所示为原始地震切片数据。为了验证重建算法的有效性,对其进行了50%高斯随机抽样,如图3(右)所示。为了更好地显示重建结果,计算重建前后的SNR 和PSNR,计算结果如表1 所示。

图3 原始地震数据和50%采样地震数据

表1 重建前后SNR 和PSNR 对比

为了验证本文中使用的K-SVD 词典学习在地震数据上的稀疏重建能力,设置了一种生成式压缩感知(GCS)[10]算法进行比较。在数据恢复之前,使用K-SVD 算法训练得到601 ×1052 的超完备字典,然后将学习到的超完备字典用作变换基,利用本文提出的迭代算法对50%采样数据进行重建,重建结果如图4(右)所示。使用生成式压缩感知算法获得的重建结果如图4(左)所示。

图4 使用GCS(左)或使用超完备字典(右)作为转换基础的重建结果

从图5 的直观重建结果可以看出,尽管生成式压缩感知的重建结果可以获得具有良好图像质量的地震数据,但是重建的地震数据与原始数据中的高频信息并不十分吻合。这反映了基于神经网络的生成式压缩感知重建方法的不确定性,不适用于地震数据重建。另一方面,不同地貌的地震数据差异很大,基于网络训练的模型无法普遍使用,导致该方法在地震数据重建领域的局限性。

为了验证K-SVD 字典学习方法在地震数据上优于传统的稀疏变换方法,本文使用地震数据重建中常用的4 种常见稀疏变换矩阵进行数据重建作为对比实验,它们分别是傅立叶变换、离散余弦变换、小波变换和曲波变换。表2 列出了每个稀疏基的重建结果的平均SNR 和PSNR。从表中可以看出,傅立叶变换矩阵的重建精度最差,而超完备字典是其中最佳的稀疏基。

表2 5 种稀疏基重建地震数据的SNR 和PSNR 对比

通过以上实验证明,本文选择的K-SVD 字典学习方法对地震数据的稀疏重建能力较强。

4.3 迭代算法对地震数据重建的影响

为了比较不同迭代算法之间的性能差异,使用超完备字典作为数据重建的稀疏变换基,并使用不同算法重建50%的采样数据10 次并计算平均值。重建效果和运行时间如表3 所示。

表3 不同迭代算法对比

从实验结果的分析可以看出,OMP 算法在重建精度及执行效率方面均不佳。SP、CoSaMP 和IHT算法的运行时间较短,但重建精度相对较低,且需要设置稀疏性。IRLS 和SAMP 不需要设置稀疏度,并且具有更好的重建效果,但是运行时间相对较长。本文方法具有最佳的重建效果。并且与IRLS 和SAMP 算法相比,大大缩短了运行时间。

本文增设不同采样率实验,基于10%～70%样本数据进行了重建。其中OMP、SP、CoSaMP、IHT 设置固定稀疏度K=50,SAMP 算法设置步长s=5。不同算法的重建效果,见图5 和图6。

图5 不同采样率下的SNR 对照图

图6 不同采样率下的PSNR 对照图

从实验结果的分析来看,需要设置稀疏度的重建算法(OMP、IHT、SP、CoSaMP)无法很好地处理不同采样率情况。而无需设置稀疏度的算法(IRLS,SAMP 和本文方法)在多种采样率情况下效果相对较好,且高采样率情况重建效果更佳。但是,SAMP算法需要设置步长,如果采样率较低且步长设置得太高,则会导致跳过最准确的稀疏度。当采样率较高时,如果步长较小,则将导致过多的迭代。在实验中,本文方法和IRLS 可以随着采样率的增加而保持稳定的增长,但是本文方法所需的时间比IRLS 所需的时间要短得多。

4.4 基于微电阻率成像数据的重建实验

为了更好地验证算法的泛化能力,增加了基于微电阻率成像数据的重建实验。微电阻率成像数据也具有空间相关性的特征,但不同于地震数据的不规则特征。微电阻率成像数据的不完整部分相对较宽,具有一定的规律性。为了证明本文中的算法也可以解决此类问题,本文使用多种算法基于真实的微电阻率成像数据进行重建,原始数据大小为360×1000。重建结果显示在图7 中。对局部进行放大,显示结果如图8 所示。

图7 原始电成像数据及不同算法重建结果

图8 局部放大后的结果

从恢复效果可以看出,除去本文算法,其他算法的重建结果都有不同程度的“失真”。由于OMP、SP、CoSaMP 和IHT 算法需要预设稀疏度,稀疏度设置不佳,重建结果会稍微模糊。在实际生产过程中,具体的稀疏度通常未知,因此这些算法不适用于该类问题。IRLS 和SAMP 算法重建结果细节较好,但是所需的计算时间太长。而且,这6 种重建方法重建结果的缺陷部位均有明显的重建痕迹。与其他算法相比,本文算法利用了空间相关信息,使重建结果细节更清晰,水平过渡更平滑,并且相比IRLS 和SAMP 大大减少了操作时间。证明了本文算法具有较强的稀疏重建能力和泛化能力。

5 结论

本文提出了一种基于时空约束压缩感知的地震数据重构算法。该算法对传统的压缩感知理论模型进行了修改,增加了时空约束,并使用K-SVD 算法基于现有地震数据训练并获得超完备字典来代替传统变换。然后,对SAMP 算法进行了改进,提出初始稀疏估计和可变步长更新策略,从而大幅减少了迭代次数,提高了重构精度。将生成式压缩感知和传统压缩感知方法作为本文算法的对照实验,并在真实地震数据及微电阻率成像数据上进行实验,结果表明了本文提出的重建方法相比其他方法具有较强的重建能力。此外,本文方法与其他方法一样,在低采样率情况下无法很好地重建。生成式压缩感知方法在本文实验中没有能够展示出很好的结果,这是由于该方法约束不足而具有很强的不确定性,仍然将基于神经网络的重建方法作为未来的研究点,这可能将会是实现基于更低采样率重建的新突破口。