基于Adam优化的二维自适应TFPF地震去噪算法

孟繁磊,范秦寅,穆丽红

(1.长春大学 电子信息工程学院,长春 130022；2.大阪大学 机械科学部门,吹田 564-0053， 日本；3.吉林省计量科学研究院 几何量室,长春 130103)

0 引 言

地震勘探是油气田矿等资源勘探的主要手段。随着地震勘探开发工作的深入,易探易采资源的减少,勘探目标向更深、更薄、更不规则的复杂地区发展的同时,对勘探的精度及保真度要求越来越高,从事勘探基础理论研究与应用的难度越来越大。一方面对于浅薄储集层(如新疆准噶尔盆地春风油田浅薄储集层)[1]、地层岩性油气藏(如塔里木盆地碳酸盐岩油气),现有的地震资料不能满足精细刻画缝洞体系、准确预测油气富集区域的开发需求[2],因此需要更高精度的新数据处理方法。另一方面,无论是为了更好地理解形成和控制陆内成矿的深部动力学过程,还是深部找寻资源的现实需求,勘查走向深部已成为发展趋势[3]。如塔里木盆地、青藏高原腹地、四川盆地、华北地区及松辽盆地的深层勘探,都需要地震勘探的新理论及交叉学科的新研究成果解决更多的工程技术难题。上述全部勘探的基础是从品质好的资料出发开展后续解释工作。压制地震勘探强随机噪声提高信噪比和分辨率是改善地震记录品质尤为重要的环节,也是地震信号研究与处理的难点与热点之一。

近年来,国内外学者就消减强随机噪声提高地震勘探资料信噪比问题进行了不断探索,并已取得了实际意义的成果。如可控方向滤波器[4]、稀疏去噪[5]、粒子滤波[6]、变分模态分解[7]、Shearlet[8]以及深度卷积网络[9]等方法。基于加窗Wigner-Ville分布(WVD：Wigner-Ville Distribution)的时频峰值滤波技术[10]TFPF(Time-Frequency Peak Filter)能处理极低信噪比情况下(-7 dB)的非平稳信号,因此非常适合处理地震数据。对于所有频率成分,传统 TFPF仅使用固定窗长进行滤波。短窗长不能有效压制随机噪声,长窗长会造成有效信号的衰减。基于这个问题,文献[11-12]提出了径向道以及自适应径向道TFPF,选择了更靠近反射同相轴的径向道方向,而非一维算法沿地震道的方向进行TFPF；Tian等[13]提出了变离心率双曲Trace变换的TFPF；2015年,Zeng等[14]证明TFPF算法可以近似等效于一个时不变低通滤波器,并结合凸优化的基本思想,提出了基于滤波器凸集的自适应滤波算法,最优化求解采用维特比算法。然而,该文献只是在一维情况下的讨论,二维信号的时空关联性有可能被忽视甚至被破坏,同时,维特比算法的速度较慢。随后,曾谦[15]提出了阵列信号时空自适应滤波算法。在一维自适应滤波框架的基础上选取方向导数作为待优化函数的惩罚项,可以兼顾信号的方向性和时空相关性。优化函数的求解使用了传统的梯度下降法,在每次迭代运算中,梯度下降根据自变量当前位置,沿着当前位置的梯度更新自变量。然而,自变量的迭代方向仅仅取决于自变量当前位置,并且目标函数自变量每个元素的学习率在迭代过程中保持不变,这对收敛速度与精度都会造成影响。

笔者在基于凸集构造的二维自适应滤波思想基础上,结合线性等效后的时频峰值滤波,提出了二维自适应TFPF算法。使用Adam算法求取自适应框架中目标函数的全局最优解,Adam算法在传统梯度下降的基础上使用了动量变量和按元素平方的指数加权移动平均变量,同时,目标函数自变量中每个元素都分别拥有自己的学习率,这可以较梯度下降法更快速稳定地逼近最优解。

1 二维自适应TFPF

1.1 TFPF原理

TFPF是一种基于瞬时频率(IF：Instantaneous Frequency)估计的非线性信号平滑算法。利用加窗的WVD能量聚集性好的优势,估计出经过调频的含噪信号的瞬时频率进行去噪。

假设一组理想观测地震信号

s=x+g

(1)

其中s∈n×1是含噪观测信号,x∈n×1是不含噪的有效信号,g∈n×1是随机噪声。值得注意的是,TFPF是一维去噪算法,需要对二维含噪地震记录进行逐道滤波。

首先,将s调制成解析信号z

(2)

其中μ>0是缩放系数。显然,s是z的瞬时频率。然后,计算关于z的加窗WVD

(3)

其中w(wτ向量形式)为窗函数,且长度为2L+1。最后,从Wi,f中得到峰值处的瞬时频率即为x估计值

(4)

(5)

其中*表示卷积运算,h为TFPF的冲激响应,且和窗函数w具有相同的长度。经证明,h是一个时不变低通滤波器。由h与w之间的映射关系易知TFPF的滤波效果被窗函数唯一决定。

1.2 二维自适应滤波框架

下面将构造滤波器输出的凸集,并在其内设计一个含有时空相关性的二次凸函数,在该函数关于某个指标取得最小值时所对应的滤波器输出即为最优解。这样,二维自适应TFPF等效于带有箱型约束的凸优化问题,并存在全局唯一的最优估计。

现选取由不同窗长得到的一组冲激响应集合,H={h1,h2,…,hK(K≥2)},则H的凸包为

(6)

其中conv(H)是包含H的最小凸集,hc∈conv(H)为不同窗长参数下TFPF冲激响应的凸组合,合理地选取不同的窗长参数使conv(H)可以完全而连续地覆盖信号的频谱范围。根据TFPF输出与冲激响应的线性卷积关系能方便地得到一组滤波器输出的凸集{Y1,Y2,…,YK},注意Yi∈n×m是经线性TFPF逐道滤波后的二维数据。这里使用矩阵向量化(vectorization,vec(·))的表示方式,将n×m的二维地震数据的列向量首尾相接转换为nm×1的列向量,所以输出的凸集可写为Y={vec(Y1),vec(Y2),…,vec(YK)}∈nm×K,同时由式(6)可知,利用hc逐道卷积二维含噪地震数据S并做向量化处理可得

y=vec{hc*S}∈conv(Y)

(7)

这样,在每个新的采样点i=1,2,…,nm上,子集yi是连续有界的,即

yi∈Ωi=[li,ui]=[minyk,i,maxyk,i]=conv(Yi)

(8)

设置Ω=Ω1×Ω2×…×Ωnm=conv(Y),便得到经向量化后的凸包。进而,二维自适应TFPF后的所有可能估计值y均在可行域Ω={y:l≤y≤u}中,其中u=sup{Y}和l=inf{Y}分别为Y的上界和下界。至此,选择最佳的估计值y∈conv(Y)等价于选择hc∈conv(H),根据希尔伯特投影理论,在可行域Ω上设计一个关于滤波输出的凸函数J(y),必有全局最小值。可行域Ω的建立过程如图1所示。

图1 最优解建立凸集过程

1.3 目标函数

在建立了滤波器输出的某个凸集Ω后,根据信号特点,需建立自适应滤波的目标函数(凸函数)。假设不含噪信号足够光滑,且至少存在一阶绝对连续的导数。因此,将地震勘探随机噪声消减问题表述为凸优化问题,也就是寻找一个定义在Ω上的关于y的目标函数并使之最小化,拟建立带惩罚项的最小二乘准则(PLS：Penalized Least Squares)作为该问题的目标函数[11]

(9)

注意到式(9)采用了矩阵向量化(vec)表示。式(9)第1项是最简单的范数逼近问题----加权的最小二乘逼近,这项表示滤波结果y对含噪信号s的拟合度,强调保幅。其中A∈nm×nm是非负对角权矩阵,表示内积运算。另外,之前假设不含噪信号足够光滑,至少存在一阶绝对连续的导数,则在最小二乘的基础上加入了罚函数(式(9)中第2项与第3项),其中λ称作平滑因子；D为差分算子,Dr和Dc分别执行对Y的行和列的差分操作,相当于对采样前的连续二维信号Y在两个坐标方向上的偏导数。也就是说：Dcy=vec(Dc0Y),Dry=vec(Dr0Y),其中

(10)

参数|p|=1,可将其理解为cosθ或sinθ,所以第2项与第3项表示了沿着arccosp或arcsinp方向的方向导数的均方值,强调这一方向的光滑度,也就是强调去噪。这样,该目标函数引入了时空关联性。

采用加权的PLS为信号在每个采样点选择符合约束的滤波器。如果权值是“自适应”的,则达成了在每个采样点自适应滤波(选择滤波器)的初衷。因此,乘在最小二乘上的权值A和乘在惩罚项上的平滑因子λ,共同决定了滤波器输出的跟踪性能与平滑性能的折衷。对A和λ的赋值方法如下

Aa,a=|ua-la-ψ|,ψ=σ(max‖Hk‖F-min‖Hk‖F),λ=(κψ)2,κ>0

(11)

(12)

当某时刻权值很大时,滤波器会倾向于选择较多的拟合度,反之则选择较多的平滑度,使滤波估计可以在时空域某方向上达到最佳平滑效果。

1.4 最优化问题的求解

目标函数J(y)：nm→是一个实值函数。笔者讨论的优化问题的含义是指能使函数J达到最小时,从约束凸集Ω中找出与之相对应的向量y,即arg minJ(y)。由前述可知,在凸集Ω上定义的凸函数有唯一最优解。

笔者采用投影Adam算法求解,Adam是基于梯度下降的优化算法[16]

y(k+1)=y(k)-αJ(y(k))

(13)

图2 以几何方式演示利用梯度法求二元单值函数J：2→的极小值点

但给定了α,由于初始点y(0)中的各个元素位置不同,则在梯度下降迭代y中的nm×1个元素时,可能导致某些元素越过目标函数最优解,这就需要一个较小的α。然而,这也会造成自变量朝最优解移动变慢[18]。

Adam使用了如下的动量变量v(k)和按梯度元素平方的指数加权移动平均变量s(k),并将它们中每个元素在初始时刻时设置为0[19]。

v(k)=β1v(k-1)+(1-β1)g(k)

(14)

s(k)=β2v(k-1)+(1-β2)g(k)⊗g(k)

(15)

(16)

(17)

y(k+1)=y(k)-g′(k)

(18)

对于式(18)的迭代算法,当梯度为零时可作为在理论上的停止规则。但在实际应用中,很难得到g′(k)=0的一步迭代,因此可以计算相邻两个迭代点差值的范数‖y(k+1)-y(k)‖,如果小于预设的阈值(如10-5),则停止迭代。Adam优化算法的基本流程如图3所示。

图3 Adam优化算法的基本流程

(19)

点PΩ[y]称为y到Ω的投影,PΩ为投影算子。可以看出,PΩ是Ω中最接近y的点。

2 仿真结果与分析

2.1 人工合成地震记录

首先给出一个复杂人工合成地震记录测试二维TFPF的去噪能力。如图4a所示,纯净信号包含5个由雷克子波形成的反射同相轴,其形状有直线轴、双曲轴、交叉轴与断轴,其频率分别为45 Hz,35 Hz,30 Hz,25 Hz和20 Hz。共56道,每道信号有1 500个采样点,采样频率为 1 000 Hz。为了更好地模拟真实情况,在纯净信号中添加了实际的地震初至前噪声,使信噪比为-4.766 1 dB,如图4b所示。从图4b中可以看到,有效信号几乎被完全淹没。分别采用传统TFPF、一维和二维自适应TFPF处理含噪信号,传统TFPF的窗长设为9点,后两种方法采用的TFPF滤波器集合参数确保得到一个合理的对比：H={h1,h1 499},对应于矩形窗分别为1和1 499点数据长度,该组滤波器的通带足以覆盖全部信号频谱。其中二维TFPF的方向导数参数p从0～1取值时,输出信噪比的变化情况如图4c所示,可以看出在p=0.1时取得最大输出信噪比。滤波结果分别如图4d～图4f所示。对比3种算法可以看出,TFPF在强实际地震噪声背景下的去噪效果并不明显,仍然有大量的随机噪声残留；两个自适应TFPF均较好地恢复出有效同相轴,但二维算法去噪效果更明显,背景更干净。

图4 复杂人工合成地震记录例子

另外二维方法使用Adam优化,图5a给出了学习率为0.005时梯度下降和Adam算法的迭代次数与损失函数值的关系,从中可以看出Adam比梯度下降算法有更快的收敛速度。在图5b中,当学习率设置为0.1时,梯度下降算法由于学习率过大导致梯度弥散,而Adam只需20次迭代即可稳定。

a 学习率为0.005 b 学习率为0.1

接着,为了更加深入地对比,对上述纯净、含噪、传统TFPF、一维与二维自适应TFPF 5个记录分别比较了F-K谱图,如图6a～图6e所示,可以明显地看出,二维自适应TFPF 的F-K谱比一维的方法更接近纯净信号,这说明二维自适应方法在复杂地形与强噪声背景下依然可以有效恢复信号。

图6 各记录F-K谱图

最后利用

(20)

分别计算不同记录的输出信噪比定量分析实验结果,式(20)中x0是无噪信号,x是滤波后的信号,M和N是矩阵的维度。实验分析结果列入表1中,由表1可看出,二维自适应TFPF的输出信噪比最高。

表1 各记录的SNR

2.2 实际地震记录

图7 实际地震记录对比

笔者通过处理实际地震资料验证了二维自适应TFPF方法压制地震勘探随机噪声的适用性。选取两幅实际地震记录。图7a给出了一个168道的共炮点地震记录,每道包含2 025个采样点,记录的采样频率为250 Hz。从图7a中可以发现,记录中同相轴的形状十分不规则,很多有效轴的清晰度和连续性都被随机噪声所干扰。采用一维和二维自适应TFPF分别处理这两幅地震资料,其滤波器集合参数均为H={h1,h2 025},二维方法使用Adam优化,初始学习率为0.1。结果分别如图 7b、图7c所示。其中一维方法处理后的记录,随机噪声并没有完全去除；而二维方法的噪声压制得更加彻底,几乎没有残留,同相轴清晰可见。有代表性的区域已用白色矩形框标出。

3 结 语

笔者针对一维自适应TFPF没用利用空间相关信息的问题,研究了基于滤波器凸集的二维自适应TFPF算法。在一维自适应目标函数框架下加入了方向导数参数,同时使用投影的Adam算法对目标函数进行优化,从而实现了快速稳定的二维滤波。通过人工合成实验和实际资料处理结果表明,笔者提出的算法充分利用了地震同相轴的横向相关性,能更清晰地恢复出有效同相轴,并且对噪声有较好的压制作用。