APP下载

基于贪婪—快速迭代收缩阈值的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)。

图8 SRTIS 法不同迭代次数下的多次波压制结果对比

图9 SRTFIS 法不同迭代次数下的多次波压制结果对比

图10 SRTGFIS 法不同迭代次数下的多次波压制结果对比

在相同的计算环境和实验数据下,对程序运行所需时间进行了测试。当达到相同的压制精度时,SRTIS、SRTFIS及SRTGFIS算法耗时如表2所示(各算法分别运行10次求得的平均值),SRTGFIS算法优于前两种,计算效率更高。

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

4 结论和认识

本文将Greedy FISTA算法引入到Radon变换压制多次波的逆问题求解中,构建了SRTGFIS算法。模拟数据和实际数据的处理结果均表明:该算法对多次波的压制效果优于LSRT算法;与SRTIS和SRTFIS算法相比,当压制精度相同时,效率更高。

随着深度学习方法的不断发展,在其他学科领域均已展现出了巨大的优势,将其引入到多次波压制,避免传统多次波压制算法中逆问题求解的时间消耗,是下一步提高多次波压制效率的研究方向。

猜你喜欢

压制步长抛物线
基于Armijo搜索步长的BFGS与DFP拟牛顿法的比较研究
巧求抛物线解析式
赏析抛物线中的定比分点问题
基于随机森林回归的智能手机用步长估计模型
基于Armijo搜索步长的几种共轭梯度法的分析对比
空射诱饵在防空压制电子战中的应用
抛物线变换出来的精彩
玩转抛物线
基于动态步长的无人机三维实时航迹规划
大型压制模用高分子粉体履平机的研制