一种类RNN的改进ISTA稀疏脉冲反褶积

2019-08-12潘树林杨海飞蒋从元秦子雨

石油物探 2019年4期

潘树林,闫柯,杨海飞,蒋从元,秦子雨

(1.西南石油大学地球科学与技术学院,四川成都610500;2.长庆油田长北作业分公司,陕西榆林719000;3.四川职业技术学院电子电气工程系,四川遂宁629000)

稀疏脉冲反褶积可利用有限带宽地震记录反演得到地下反射系数序列,是一种常用的地震资料高分辨率处理方法[1]。该方法在传统的最小二乘反褶积的基础上,加一个稀疏先验约束项,通过求解带约束的目标函数,进而得到高分辨率的处理结果。求解该目标函数的方法主要有以下两类。第一类为匹配追踪类贪婪迭代算法[2-5],直接使用L0范数作为稀疏约束项,通过枚举匹配迭代的方法求解未知参数。该类方法绝对收敛,但运算量巨大。第二类为CHEN等[6]提出的基追踪类凸优化算法,此类方法利用L1范数代替L0范数,以便使用线性规划来求解。该类方法仅在保证期望足够稀疏时,才能获得理想的反演结果。为了进一步改善基追踪类方法对数据稀疏度的要求,张繁昌等[7]、曹静杰等[8]通过引入地质模型约束和非凸Lp范数正则化方法对其进行了改进。基追踪类方法结构简单,便于实现,计算效率高,是当前较主流的求解方法。应用该类方法时,要得到好的反褶积效果,地震子波必须准确,但在实际生产中,难以从地震数据中提取出非常准确的地震子波,因而制约了稀疏脉冲反褶积方法的应用效果。

针对上述问题,国内外学者进行了很多研究。第一类方法为DONOHO等[9]采用过完备字典的方法,通过将地震子波扩充为包含多个频率成分的过完备字典,消除了单个频率成分子波带来的缺点。这类算法最大的特点在于利用多个频率成分的过完备子波字典代替单个子波字典,求解结果的精度取决于子波字典的准确性。第二类方法为字典学习的方法[10-12],通过学习修改的方式来更新地震子波,但是该类方法只是利用误差来交替修改反射系数与地震子波,而不是直接通过误差反传来修改,因此计算误差较大。

迭代阈值收缩算法(iterative shrinkage-thresholding algorithm,ISTA)是基追踪类方法中一种重要的求解算法,其收敛性受输入子波的影响。循环神经网络(recurrent neural network,RNN)是一种递归神经网络,可以根据输入数据进行自主学习,改善处理结果,在很多领域取得了较好的应用效果。本文通过对比ISTA结构和RNN网络结构,结合两种方法提出了一种类RNN的改进ISTA稀疏脉冲反褶积方法,该算法是将ISTA转换为类RNN网络结构,并引入RNN中BPTT的误差反传的思想修改地震子波,克服了子波不准确造成的反演不收敛问题。最后用理论模型数据和吐哈盆地雁木西地区的实际资料对本文方法进行了应用测试。

1 类RNN的改进ISTA原理

1.1 ISTA

稀疏脉冲反褶积可以表达为:

(1)

其中,y表示地震记录向量,W表示由地震子波组成的矩阵,x表示待求解反射系数向量,λ为正则化参数。

BECK等[13]提出了ISTA,用于求解具有L2-L1结构的复合函数。利用迭代软阈值的方法求解(1)式,可得:

(2)

(3)

ISTA迭代流程如图1所示。

图1 ISTA迭代流程示意

显然,ISTA求解反射系数的关键在于地震子波矩阵We的准确确定,但是在实际应用中,很难准确确定地震子波,所以传统的ISTA在反褶积应用过程中因为地震子波的问题而存在较大的局限性。

1.2 类RNN的改进ISTA原理

RNN是一种节点定向连接成环的人工神经网络,其本质特征是在处理单元之间既有内部的反馈连接又有前馈连接,可以在网络运行中对内部参数不断进行优化,这是一类适用于处理序列数据的神经网络,它与基础的神经网络最大的不同就是在层之间的神经元也建立了权值连接[14-20]。

图2为RNN基本结构及展开形式。左侧为RNN基本结构,垂向黑色粗箭头表示“循环”,这意味着RNN的循环体现在隐藏层;右侧为RNN展开结构,隐藏层各神经元之间通过不同权值(见小箭头上的权值符号)连接,即上一步处理的隐藏层(Nt-1)将会影响到下一步处理的隐藏层(Nt)。RNN结构包含了RNN正向计算(黑色细箭头方向)和反向误差修正(红色细箭头方向)两个核心内容,其中M和U为计算权值,V为激活函数。本节提出的类RNN的改进ISTA同样包含了正向计算和反向误差修正两个重要环节。

对比图1和图2可以发现,ISTA与RNN在结构上有着相似之处。将ISTA中的每一次迭代看作一个时间层,输入向量中的每一个值看作是一个神经元,软阈值函数hθ等价于激活函数V,那么ISTA便可以看作一个类似的RNN网络,如图3所示。

图2 RNN基本结构及展开形式

图3 类RNN的改进ISTA

由图3可知,矩阵S为循环层参数,即为时间方向参数,We为各层神经元之间的参数,即为层方向参数。

2 改进算法正反向计算方法

2.1 正向算法

由图2可知,RNN的正向计算可表示为:

(4)

式中:Et为中间变量。

公式(3)可以改写为如下形式:

(5)

其中,公式(5)中的y相当于公式(4)中的输入xt,S相当于隐藏层权值M,xk相当于公式(4)中的隐藏层参数Nt,hθ相当于激活函数。

由此可以给出改进ISTA的正向计算流程(图4),其计算步骤与常规的ISTA一致。

2.2 反向误差修正算法

图4 正向算法流程

图5中,xr为训练时给定的反射系数,xk为训练过程中计算获得的反射系数,h′为软阈值求导,Hk,Hk+1为梯度下降过程中的控制函数。由图5可知,需要先通过迭代逐步反向求解各层误差,然后求和,得到关于初始输入的误差δB、δθ、δWe,进一步利用图5流程最后一步的梯度下降方法更新这3个参数。

图5 改进算法的反向误差修正流程

结合正向算法与反向算法,便可以根据训练数据多次迭代修改与地震子波有关的参数,得到更好的收敛效果,更加有利于后续的反褶积计算。

3 理论模型数据试算

首先合成了一个不含噪声的地震单道数据,使用50Hz零相位Ricker子波与第1和第2个反射系数进行褶积,40Hz零相位Ricker子波与第3和第4个反射系数进行褶积,30Hz零相位Ricker子波与第5和第6个反射系数进行褶积,最终得到如图6所示的理论地震道,该数据包含512个采样点,6个反射位置。

图6 理论模型a 地震子波; b 反射系数; c 合成地震记录

利用图6中的模型数据,令初始输入的子波矩阵为零相位40Hz的Ricker子波,结合前文所述的正向算法与反向算法,对初始子波进行修正,经过本文算法优化后的子波如图7所示,从图7可以看出,虽然输入仅有40Hz的子波,但是经过算法自我修正后,可以得到与图6a多个子波接近的结果,且经过算法修正后获得的子波集合与期望子波集合基本吻合。使用优化后的子波矩阵进行反褶积处理,得到如图8所示的结果。给定同样的初始子波(40Hz的Ricker子波),进行传统ISTA反褶积方法处理,结果如图8中蓝线所示。

由图8可见,由于传统ISTA不具备子波修正能力,当采用40Hz子波作为输入子波时,算法收敛后计算结果与期望输出存在明显差异,反演结果出现了许多虚假信息。而类RNN的改进ISTA在计算过程中不断优化子波集,最终在不同的位置获得不同的最优子波。在此基础上计算的结果明显优于传统ISTA。这说明类RNN的改进ISTA能够通过反向修正子波集合,消除初始输入子波不准确对反褶积带来的影响。

为了说明本文方法的抗噪性,在图6c合成地震记录的基础上加入随机噪声,得到一个如图9所示的低信噪比地震记录,传统ISTA与本文方法的反褶积结果如图10所示。

图7 原始子波与优化后的子波

图8 传统ISTA和改进ISTA反褶积处理的结果

图9 含噪声地震记录(信噪比=-4.55dB)

图10 含噪声地震记录反褶积结果

由图10可以看出,传统的ISTA在地震资料信噪比较低时,难以获得准确的反褶积结果,而采用改进ISTA可获得较为可靠的反褶积结果,说明该方法具有良好的抗噪能力。

4 实际地震资料应用

图11为吐哈盆地雁木西地区T3井区Line 299线主要目的层段局部地震剖面。分别采用传统ISTA和改进ISTA对其进行反褶积处理,结果如图12所示。

对比图11和图12可以看出,采用传统ISTA进行反褶积处理后的剖面整体同相轴变多,变细,且以井旁地震道分辨率提高最为明显,但是距离井位置越远处理效果越不明显(图12a)。相比之下,采用改进ISTA的反褶积剖面中,不论是井旁道还是距离井位置较远的地震道,分辨率均得到了显著提升,且与测井合成记录的吻合程度大大提高,说明采用改进ISTA的反褶积处理结果可靠。从图11和图12中井旁红色方框中局部放大的地震剖面(图13)可以看出,采用传统ISTA处理后,原剖面中部分无法分辨的同相轴得到了有效分离,但是分辨率仍然偏低,而采用改进ISTA处理后的剖面,不仅使原来无法分离的同相轴得到了有效分离,而且同相轴连续性增强,与井资料对比良好,分辨率明显提高。