基于贪婪—快速迭代收缩阈值的Radon变换及其在多次波压制中的应用

2022-12-09林柏栎刘书妍

石油地球物理勘探 2022年6期

张全雷芩林柏栎彭博*② 刘书妍②

(①西南石油大学计算机科学学院,四川成都 610500；②油气藏地质及开发工程国家重点实验室(西南石油大学)，四川成都 610500；③电子科技大学信息与通信工程学院，四川成都 611731)

0 引言

1917年，奥地利数学家Radon[1]首次提出Radon变换，经过众多学者的研究，将其引入各领域，并得到了广泛的使用。20世纪70年代，Claerbout等[2]首次将Radon变换引入地震数据处理，为Radon变换在地震勘探中的发展奠定了基础。

Radon变换将地震数据沿特定路径进行曲线积分，目的是将数据中的信号叠加为Radon域的稀疏信号，便于识别与分离。为了提高Radon变换压制多次波的精度，不少学者进行了相关的研究。Thorson等[3]提出了Radon变换最小二乘法和随机反演的理论，将Radon变换建模为一个稀疏逆问题，使Radon变换的分辨率得到一定程度的提高。李元钦等[4]结合平面映射，提出双曲Radon正、反变换，减少变换噪音和能量弥散，提高了地震资料的处理质量。Sacchi等[5-6]在频率域利用柯西稀疏约束提高了Radon变换的分辨率。牛滨华等[7]提出次数为“2”的多项式Radon变换。刘喜武等[8]采用最小二乘反演法实现高分辨率抛物线Radon变换和双曲线Radon变换，提高了精度和分辨率。Ethan等[9]提出加权最小二乘Radon变换，在压制多次波过程中能较好地保持振幅。刘保童等[10]提出一种频率域Radon变换方法，有效地提高了地震数据在Radon域的分辨率。Schonewille等[11]证明了Radon变换在时间域比在频率域更加稀疏，故Radon变换在时间域的分辨率比在频率域高，提高解复用效率。薛亚茹等[12]提出基于Radon变换和正交多项式变换的多方向正交多项式变换的方法，提高了对一次波与多次波的剩余时差的分辨率，在有效压制多次波的同时保留了一次波的AVO特性，增强了Radon变换的保幅特性。Lu[13]提出基于迭代二维模型收缩的混合时频域快速稀疏时不变Radon变换方法，利用混合时频域的优势，该算法既在时间域中提高了Radon模型的稀疏性，又在频率域得到更高的计算效率。薛亚茹等[14]提出高分辨率稀疏Radon变换和正交多项式变换结合的高阶稀疏Radon变换，能有效分离一次波和多次波。Sun等[15]提出基于光滑的 L0范数的高阶、高分辨率频率-曲率域Radon变换，将正交多项式与Radon变换相结合，该算法具有更好的聚焦特性，并在压制多次波的同时能更好地保存一次波的AVO特性，且有更高的计算效率。薛亚茹等[16]提出应用加权迭代软阈值算法的高分辨率Radon变换，将迭代重加权最小二乘法的加权矩阵思想引入迭代软阈值算法中，提高了Radon变换分辨率。孙宁娜等[17]提出基于多窗口联合优化的多次波自适应相减方法，与传统的基于单窗口匹配的多次波自适应相减方法相比，该方法可更有效地压制残余多次波并保护一次波，且计算效率与基于小窗口匹配的传统多次波自适应相减方法相当。马继涛等[18]提出三维高精度保幅Radon变换多次波压制方法，该方法可获得高分辨率的模型域数据，能有效分离一次波与多次波，同时多项式拟合可保护有效波的振幅，高保真地实现多次波的压制。在Radon变换压制多次波的速度提升方面，很多学者也进行了相关的研究。Hampson[19]提出抛物线Radon变换压制多次波，进一步提升了计算效率。熊登等[20]提出一种新的混合域高分辨率抛物线Radon变换，对于多次波的压制既有很高的计算效率又有较高的分辨率。Brahim等[21]提出一种基于奇异值分解的改进的抛物线Radon变换，克服了频域Radon变换实现所需的繁重计算，具有更快的计算速度。张全等[22]提出CPU-GPU异构平台的抛物线Radon变换并行算法，大幅提高了Radon变换压制多次波的计算效率。随着深度学习的发展，不少学者将多次波压制与深度学习相结合。刘小舟等[23]提出数据增广的编解码卷积网络地震层间多次波压制方法，结合去噪卷积神经网络(DnCNN)和U形全卷积神经网络(U-Net)的优势，搭建了适合层间多次波压制的深层编、解码网络，并且利用数据增广提高网络的泛化能力，该网络能够很好地压制层间多次波、保护一次波，提高了多次波的压制效率。

本文在Lu[13]的研究基础上，将Liang等[24]提出的贪婪—快速迭代收缩阈值算法(Greedy Fast Ite-rative Shrinkage Thresholding Algorithm，Greedy FISTA)应用于混合域抛物线Radon变换，构建了一种快速稀疏时不变Radon变换进行地震多次波压制。实验结果表明，本文方法有效地提高了Radon变换压制多次波的精度和效率。

1 Radon变换多次波压制基本原理

在地震数据处理中，常采用抛物线Radon变换实现一次波与多次波分离。动校正后，道集中的一次波被拉平，多次波呈“抛物线”形态，根据二者在Radon域的形态差异实现多次波压制。时域与Radon域地震数据可用数学模型表示为

d=Lm

(1)

式中：d为时空域地震数据；m为地震数据在Radon域的表示；L为Radon变换算子矩阵。通常情况下，d为已知数据，m为未知数据，L通常不可逆，常用共轭转置矩阵LH近似L的逆矩阵，故式(1)在数学上为一逆问题。求解逆问题的经典方法为最小二乘(LS)法

(2)

L通常是病态的，且LS法不适用于求解病态方程，常添加正则化项解决此问题。基于L1范数的稀疏正则化是目前常用的提高Radon变换分辨率的方法

(3)

式中λ>0为正则化系数。在时域和频域构成的混合域(简称混合域)实现Radon变换

d=F-1[LF(m)]

(4)

式中：F(·)为傅里叶变换算子；F-1(·)为傅里叶逆变换算子。频率域的Radon变换算子为[26]

(5)

式中：i=0，1，…，Nq，其中Nq为曲率q的个数；j=0，1，…，Nx，其中Nx为地震数据道数；f为频率；xj为第j道的炮检距。

Radon变换压制多次波在混合域的目标方程为

(6)

该目标方程的求解是影响多次波压制效率的重要因素。

2 Radon变换算法及其改进

对式(6)的求解，传统的迭代重加权最小二乘(Iterative Reweighted Least Squares，IRLS)算法涉及大规模的求逆运算，耗时长。而基于梯度的迭代算法，计算复杂度小，且算法结构简单。本文主要将当前应用较广的迭代收缩阈值算法(Iterative Shrinkage-Thresholding Algorithm，ISTA)、快速迭代收缩阈值算法(Fast Iterative Shrinkage-Thresho-lding Algorithm，FISTA)以及Greety FISTA与抛物线Radon变换结合对多次波进行压制。

2.1 基于迭代收缩阈值的Radon变换算法

ISTA是经典梯度算法的扩展，计算简单，是目前最常用的方法之一。Lu[13]将二维ISTA引入Radon变换，提出了基于迭代二维模型收缩的混合域快速稀疏时不变Radon变换(an accelerated sparse time-invariant Radon transform in the mixed frequency-time domain based on iterative 2D model shrinkage，SRTIS)。SRTIS算法流程如下。

(1)输入时空域地震数据d，最大迭代次数K；

(2)初始化模型m0，当前迭代次数k=0，迭代步长t>0；

(3)若k≤K，迭代更新

mk+1=Tα{mk+2tF-1{(LHL)-1×

LH[F(d)-LF(mk)]}}

k+1→k

(4)输出 Radon域地震数据mk。

算法中，一般要求迭代步长t∈(0,1/λmax(LHL))，其中λmax(·)表示求最大特征值，Tα为近似算子，定义为

(7)

(8)

相较于一些传统方法，SRTIS提高了Radon模型的稀疏性，从而提高了多次波的压制精度，同时还有更快的收敛速度。

2.2 基于快速迭代收缩阈值的Radon变换算法

ISTA每次迭代只需考虑当前点的信息更新迭代起点，导致算法迭代速度较慢。FISTA在ISTA的基础上考虑当前点和前一点更新迭代起点，比ISTA具有更快的收敛速度。Li等[26]将该算法引入到Radon变换，提出了基于快速迭代收缩的混合域快速稀疏时不变Radon变换(an accelerated sparse time-invariant Radon transform in the mixed frequency-time domain based on fast iterative shrinkage-thresholding algorithm，SRTFIS)。SRTFIS算法流程如下。

(1)输入时空域地震数据d，最大迭代次数K；

(2)初始化模型m0，当前迭代次数k=0，迭代步长t>0，p0=1,y=m0；

(3)若k≤K，迭代更新

(4)输出Radon域地震数据mk。

算法中一般要求t=1/λmax(LHL)。

与SRTIS相比，SRTFIS使用yk代替mk，即在迭代步骤中，SRTFIS算法的迭代点mk不仅仅依赖于前一个迭代点mk-1，还在yk上使用了前两点{mk-1,mk-2}的线性组合，提高了收敛速度。惯性参数ak用于控制mk-mk-1的增长速度。

2.3 基于贪婪—快速迭代收缩阈值的Radon变换算法

ISTA计算简单，但在收敛速度上较慢；FISTA的收敛速度虽然快于ISTA，但对于FISTA，当ak→1时，会引起振荡现象，对收敛速度造成一定的减缓。海量地震数据的处理，对算法的处理效率提出了更高的要求。Liang等[24]提出的Greedy FISTA在保持算法计算简单的同时，全局收敛速度比ISTA、FISTA更快，有更好的收敛性能，本文将该算法应用于稀疏Radon变换算法中，提出一种基于贪婪的快速迭代收缩的混合域快速稀疏时不变Radon变换(an accelerated sparse time-invariant Radon transform in the mixed frequency-time domain based on greedy fast iterative shrinkage-thresholding algorithm，SRTGFIS)。SRTGFIS算法流程如下。

(1)输入时空域地震数据d，最大迭代次数K。

(2)初始化模型m0，当前迭代次数k=0，迭代步长t>0,y1=m0,步长t的收缩因子ξ<1,算法收敛的控制因子S>1。

(3)若k≤K，迭代更新

yk=(mk-mk-1)

mk+1=Tα{yk+2tF-1{(LHL)-1×

LH[F(d)-LF(yk)]}}

k+1→k

重启条件：若(yk-mk+1)T(mk+1-mk)≥0，则

yk=mk；

收敛条件：若‖mk+1-mk‖≥S‖m1-m0‖，则

t=max{ξt,t}

k+1→k

(4)输出Radon域地震数据mk。

该算法中，一般要求t∈[1/λmax(LHL),2/λmax(LHL)]。

为了解决FISTA中ak→1引起的振荡问题，Greedy FISTA对此引入了重启动方法，将pk重置为1，迫使ak从0开始增加，当ak→1引起下一次振荡时，进行重启，如此循环。为了缩短重启动的间隔，减小振荡周期，以此提高收敛速率，将ak取为常数1，但由于ak控制迭代步长t的大小，ak为常数将导致初始迭代步长过大，容易造成算法发散，因此Greedy FISTA算法还有一个保证收敛的条件。

整个重启动框架在保证收敛的情况下，缩短了重启间隔和振荡周期，相较于ISTA、FISTA算法，该算法提高了收敛速度。

3 实验数据

3.1 模拟数据

本文实验的模拟数据来自SeismicLab工具包合成的多次波模拟数据(图1)。该模拟数据只有一个地震道集，共49道，每道有1001个样点，采样间隔为4ms。分别将频率域最小二乘法Radon变换(LSRT)、SRTIS、SRTFIS、SRTGFIS四种算法用于模拟数据进行比较(图2)。

图1 合成的多次波全波场模拟数据

比较四种算法对模拟数据多次波的压制效果(图2)，其中SRTIS、SRTFIS和SRTGFIS算法的K均为20。由图可见，LSRT的压制效果一般，处理之后还有多次波的残余(图2a红色箭头所指)，而SRTIS(图2b)、SRTFIS(图2c)和SRTGFIS(图2d)均能完全压制多次波。

图2 不同算法对模拟数据的多次波压制效果对比

图3 三种算法的收敛速度比较

在多次波压制精度相当的情况下，比较三种算法不同迭代次数的多次波压制结果(图4)，对比三种算法对模拟数据的计算速度(表1)。图4a为SRTIS在不同迭代次数下压制多次波的结果，当迭代到第7次时，多次波残余明显(红色箭头所示，下同)，迭代到第16次时，多次波基本被压制，但依然有残余；图4b为SRTFIS在不同迭代次数下压制多次波的结果，迭代到第7次时，多次波有少量残余，迭代到第8次，多次波基本被压制；图4c为SRTGFIS在不同迭代次数下压制多次波的结果，迭代到第6次时，压制效果好于SRTFIS，迭代到第7次时，多次波已基本被压制。

图4 不同算法不同迭代次数多次波压制结果对比

同时，在相同的计算环境和实验数据下，对程序运行时间进行了测试。当达到相同的压制精度时，SRTIS、SRTFIS和SRTGFIS三种算法耗时(由各算法分别运行10次求得的平均值)见表1，可见，SRTGFIS较其他两种算法计算效率更高。

表1 模拟数据三种算法达到相同压制精度所需迭代次数和时间对比

3.2 实际数据

实际数据来自SeismicLab工具包，只有1个地震道集，包含92道，每道401个样点，采样间隔为4ms。

图5为实际包含多次波与一次波的全波场数据。对于SRTIS、SRTFIS和SRTGFIS算法，当t=1/λmax(LHL)时，步长太小，多次波的压制效率均太低。考虑模拟数据的测试情况，对于实际数据，将步长设置为0.10时，SRTIS、SRTFIS和SRTGFIS三种算法分别需要迭代20、14、10次才能达到较好的多次波压制效果，故初始化Radon模型相关参数为：m0=0，t=0.10，K=20。LSRT、SRTIS、SRTFIS、SRTGFIS四种算法多次波的压制结果如图6所示。LSRT算法的压制结果中还有残留的多次波(图6a红色箭头所指)；SRTIS(图6b)、SRTFIS(图6c)和SRTGFIS(图4d)算法的多次波压制效果均优于LSRT算法。

图5 实际数据

图6 四种方法对实际数据的多次波压制结果对比

实际数据三种算法的收敛曲线(图7)显示，SRTIS的收敛速度最慢，SRTGFIS算法的收敛速度最快，由于ak=1导致初始迭代步长过大，引起收敛曲线的波动。

图7 三种算法的收敛速度比较

在多次波压制精度相当的情况下，分别比较三种算法不同迭代次数下的多次波压制效果(图8～图10)。SRTIS算法在迭代到第18次时，多次波有少量残余(图8中红色箭头所指，下同)，迭代到第20次时，多次波得到基本压制；SRTFIS算法在迭代到第12次时，多次波有少量残余，迭代到第14次时，多次波基本压制(图9)；SRTGFIS算法在迭代到第8次时，多次波有少量残余，迭代到第10次时，多次波得到很好压制(图10)。