石油物探  2019, Vol. 58 Issue (4): 533-540  DOI: 10.3969/j.issn.1000-1441.2019.04.007
0
文章快速检索     高级检索

引用本文 

潘树林, 闫柯, 杨海飞, 等. 一种类RNN的改进ISTA稀疏脉冲反褶积[J]. 石油物探, 2019, 58(4): 533-540. DOI: 10.3969/j.issn.1000-1441.2019.04.007.
PAN Shulin, YAN Ke, YANG Haifei, et al. A sparse spike deconvolution method based on Recurrent Neural Network like improved Iterative Shrinkage Thresholding Algorithm[J]. Geophysical Prospecting for Petroleum, 2019, 58(4): 533-540. DOI: 10.3969/j.issn.1000-1441.2019.04.007.

基金项目

国家自然基金项目(NSFC 41674095)、“油气藏地质及开发工程”国家重点实验室开放基金项目(PLN201733)和天然气地质四川省重点实验室开放基金项目(2015trqdz03)共同资助

作者简介

潘树林(1979—), 男, 博士, 副教授, 主要从事观测系统优化设计、地震数据处理新方法和微地震监测技术等方面的研究工作。Email:149769166@qq.com

文章历史

收稿日期:2019-01-18
改回日期:2019-04-22
一种类RNN的改进ISTA稀疏脉冲反褶积
潘树林1 , 闫柯1 , 杨海飞2 , 蒋从元3 , 秦子雨1     
1. 西南石油大学地球科学与技术学院, 四川成都 610500;
2. 长庆油田长北作业分公司, 陕西榆林 719000;
3. 四川职业技术学院电子电气工程系, 四川遂宁 629000
摘要:稀疏脉冲反褶积方法对提高地震资料分辨率有着重要作用, 迭代阈值收缩算法(ISTA)是其核心算法, 首先利用地震数据提取子波, 再利用ISTA求解反射系数。当地震子波提取不准确时, 反褶积效果不理想。为此, 在ISTA基础上, 结合循环神经网络(RNN)中反向传播(BPTT)的思想, 研究形成了一种类RNN的改进ISTA稀疏脉冲反褶积方法。该算法首先使用常规手段从实际地震数据中提取地震子波, 构建反褶积的子波字典; 然后将构建的地震子波字典作为已知的初始条件, 结合ISTA求取的反射系数; 再根据BPTT算法思想, 将求取的反射系数与子波褶积并与实际数据进行比较, 反向修改地震子波; 最终, 经过多次迭代修改获得合理的地震子波字典, 并利用该地震子波字典求解实际地震数据的反射系数序列。为验证算法的有效性, 采用不同信噪比的理论地震记录, 给定存在较大误差的初始子波, 进行了反褶积计算。采用传统的ISTA和类RNN的改进ISTA进行对比处理, 结果表明, 改进ISTA具有较好的抗噪能力和子波自适应能力, 可使实测地震资料的有效频带拓展约1.5倍, 能够较好地适应实际地震资料的反褶积处理。
关键词稀疏脉冲反褶积    分辨率    ISTA    地震子波    信噪比    循环神经网络    反向传播    
A sparse spike deconvolution method based on Recurrent Neural Network like improved Iterative Shrinkage Thresholding Algorithm
PAN Shulin1, YAN Ke1, YANG Haifei2, JIANG Congyuan3, QIN Ziyu1     
1. School of Earth Science and Technology, Southwest Petroleum University, Chengdu 610500, China;
2. North Operation Branch of Changqing Oilfield, PetroChina, Yulin 719000, China;
3. Department of Electrical and Electronic Engineering, Sichuan Vocational and Technical College, Suining 629000, China
Foundation item: This research is financially supported by the National Natural Science Foundation of China (Grant No.NSFC 41674095), the Open Projects Fund of the State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation (Grant No.PLN201733) and the Open Projects Fund of the Natural Gas and Geology Key Laboratory of Sichuan Province (Grant No.2015trqdz03)
Abstract: The sparse spike deconvolution plays an important role in improving the resolution of seismic data.A successful deconvolution requires accurate wavelet data for the calculation of the reflection coefficient, which is performed using the core algorithm of the Iterative Shrinkage Thresholding Algorithm (ISTA).An inaccurate extraction of wavelet data can affect the result of the deconvolution.In this study, an RNN-like ISTA algorithm is proposed, which combines the concept of Back-Propagation Through Time (BPTT) in Recurrent Neural Network (RNN) with the traditional ISTA algorithm.First, seismic wavelets are extracted from the actual seismic data to construct a wavelet dictionary, which is taken as the initial condition, and the reflection coefficient is calculated using ISTA algorithm.Subsequently, the reflection coefficient is convoluted with the wavelets to construct seismic data by means of the BPTT algorithm, and a seismic wavelet correction is performed by comparing the obtained seismic data with the actual data.Finally, a reasonable seismic wavelet dictionary is obtained after several iterations, which can be used to calculate the actual reflection coefficient series.Tests on theoretical seismic records with different signal-to-noise ratio showed that the improved algorithm has a better anti-noise and wavelet adaptive abilities than the conventional ISTA algorithm.Moreover, the frequency band of the actual seismic data can be expanded by about 1.5 times.These results demonstrate that the proposed method can be successfully applied to the deconvolution of actual seismic data.
Keywords: sparse spike deconvolution    resolution    Iterative Shrinkage Thresholding Algorithm (ISTA)    seismic wavelet    SNR    Recurrent Neural Network (RNN)    Back-Propagation Through Time (BPTT)    

稀疏脉冲反褶积可利用有限带宽地震记录反演得到地下反射系数序列, 是一种常用的地震资料高分辨率处理方法[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

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

$ \underset{x}{\operatorname{argmin}}\left(\frac{1}{2}\|\boldsymbol{W x}-\boldsymbol{y}\|_{2}^{2}+\lambda\|\boldsymbol{x}\|_{1}\right) $ (1)

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

显然, (1)式由L2范数和L1范数两部分组成, 前者$\frac{1}{2}\|\boldsymbol{W x}-\boldsymbol{y}\|_{2}^{2}$表示估计反射系数对应的地震记录与实际地震数据的最小二乘误差, 后者λx1表示估计反射系数中含有的有效反射系数的个数, 这意味着估计的反射系数不仅需要满足对应地震记录与原始地震数据误差达到最小, 而且应满足反射系数达到个数最少的条件。

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

$ x_{k+1}=h_{\theta}\left(x_{k}-\frac{1}{L} \boldsymbol{W}^{\mathrm{T}}\left(\boldsymbol{W} x_{k}-\boldsymbol{y}\right)\right) $ (2)

式中:WTW的转置矩阵; θ为阈值; 令p=xk$\frac{1}{L} \boldsymbol{W}^{\mathrm{T}}\left(\boldsymbol{W} x_{k}-\boldsymbol{y}\right)$, 则$h_{\theta}(p)=\left\{\begin{array}{l}{p+\theta, p \leqslant-\theta} \\ {0, |p| \leqslant \theta} \\ {p-\theta, p \geqslant \theta}\end{array}\right.$为软阈值函数; L为Lipchitz常数(L>0);xk+1xk为第k+1、k次迭代的最佳反射系数。令$\boldsymbol{W}_{\mathrm{e}}=\frac{1}{L} \boldsymbol{W}^{\mathrm{T}}, \boldsymbol{S}=$ $\boldsymbol{I}-\frac{1}{L} \boldsymbol{W}^{\mathrm{T}} \boldsymbol{W}$, 公式(2)可化简为:

$ x_{k+1}=h_{\theta}\left(\boldsymbol{W}_{e} \boldsymbol{y}+\boldsymbol{S} x_{k}\right) $ (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正向计算(黑色细箭头方向)和反向误差修正(红色细箭头方向)两个核心内容, 其中MU为计算权值, V为激活函数。本节提出的类RNN的改进ISTA同样包含了正向计算和反向误差修正两个重要环节。

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

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

图 3 类RNN的改进ISTA

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

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

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

$ \begin{aligned} \boldsymbol{E}_{t} &=\boldsymbol{U} x_{t}+\boldsymbol{M} N_{t} \\ O_{t} &=V\left(\boldsymbol{E}_{t}\right) \end{aligned} $ (4)

式中:Et为中间变量。

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

$ \begin{array}{l}{\boldsymbol{E}_{k}=\boldsymbol{W}_{\mathrm{e}} \boldsymbol{y}+\boldsymbol{S} x_{k}} \\ {x_{k+1}=h_{\theta}\left(\boldsymbol{E}_{k}\right)}\end{array} $ (5)

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

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

图 4 正向算法流程
2.2 反向误差修正算法

这里采用BPTT算法[20], 包含3个步骤:①正向计算每一个神经元的输出值; ②反向计算出每一个神经元的误差值, 它可以表示为误差函数对神经元加权输入的偏导数; ③计算权值参数的梯度, 并利用梯度下降法更新权重。设第lt时刻的误差值为δtl, 由BPTT算法可知, δtl沿着两个方向反向传播, 一个方向是沿着层反方向传递到上一层网络, 可得δtl-1, 这部分误差只和权重矩阵U有关; 另一个是沿着时间反方向传递到初始时刻t-1, 可得δ1l, 该部分误差只和权重矩阵M有关。参考BPTT算法, 结合图 3和公式(5), 建立如图 5所示的改进算法反向误差修正流程。

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

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

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

3 理论模型数据试算

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

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

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

图 7 原始子波与优化后的子波
图 8 传统ISTA和改进ISTA反褶积处理的结果

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

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

图 9 含噪声地震记录(信噪比=-4.55 dB)
图 10 含噪声地震记录反褶积结果

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

4 实际地震资料应用

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

图 11 雁木西地区T3井Line 299线主要目的层段局部地震剖面
图 12 采用传统ISTA(a)和改进ISTA(b)进行反褶积处理的结果(Line299)

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

图 13 采用不同方法反褶积结果剖面细节对比(Line299) a原始剖面; b传统ISTA; c改进ISTA

分别对图 13中的3个剖面进行频谱分析, 结果如图 14所示。显然, 基于改进ISTA的反褶积处理结果有效频带宽度从原始剖面的10~40 Hz拓展为5~50 Hz左右, 明显优于传统ISTA。

图 14 频谱对比
5 结论

1) 与传统的基于ISTA的稀疏脉冲反褶积方法相比, 采用类RNN的改进ISTA实施反褶积处理, 可以获得更加稳定、可靠的反褶积处理结果, 且方法的抗噪性更佳。

2) 基于类RNN的改进ISTA的稀疏脉冲反褶积方法可极大提高地震资料分辨率。实际资料的处理结果表明, 采用该方法可使原始地震资料的有效频带拓展约1.5倍, 而且除拓展高频外, 还可以极大地保护低频有效信息, 处理结果更有利于进一步储层精细描述或反演处理。

参考文献
[1]
余江奇, 曹思远, 陈红灵, 等. 改进阈值的Curvelet变换稀疏反褶积[J]. 石油地球物理勘探, 2017, 52(3): 426-433.
YU J Q, CAO S Y, CHEN H L, et al. Sparse deconvolution based on Curvelet transform of improved threshold[J]. Oil Geophysical Prospecting, 2017, 52(3): 426-433.
[2]
梁东辉, 陈生昌. 基于L0范数稀疏约束的地震数据反褶积[J]. 石油物探, 2014, 53(4): 397-403.
LIANG D H, CHEN S C. Deconvolution of seismic data based on L0 norm sparse constraint[J]. Geophysical Prospecting for Petroleum, 2014, 53(4): 397-403. DOI:10.3969/j.issn.1000-1441.2014.04.004
[3]
MALLAT S G, ZHANG Z. Matching pursuits with time-frequency dictionaries[J]. IEEE Transactions on Signal Processing, 1993, 41(12): 3397-3415. DOI:10.1109/78.258082
[4]
CAI T T, WANG L. Orthogonal matching pursuit for sparse signal recovery with noise[J]. IEEE Transactions on Information Theory, 2011, 57(7): 4680-4688. DOI:10.1109/TIT.2011.2146090
[5]
李海山, 杨午阳, 田军, 等. 匹配追踪煤层强反射分离方法[J]. 石油地球物理勘探, 2014, 49(5): 866-870.
LI H S, YANG W Y, TIAN J, et al. Coal seam strong reflection separation with matching pursuit[J]. Oil Geophysical Prospecting, 2014, 49(5): 866-870.
[6]
CHEN S, DONOHO D L, SAUNDERS M A. Atomic decomposition by basis pursuit[J]. SIAM Journal on Scientific Computing, 1998, 20(1): 33-61.
[7]
张繁昌, 彭德木, 张营革, 等. 基于对偶对数障碍规划算法的基追踪反演[J]. 石油物探, 2017, 56(2): 273-279.
ZHANG F C, PENG D M, ZHANG Y G, et al. A basis pursuit inversion method based on the primal-dual log-barrier programming algorithm[J]. Geophysical Prospecting for Petroleum, 2017, 56(2): 273-279. DOI:10.3969/j.issn.1000-1441.2017.02.014
[8]
曹静杰. 基于广义高斯分布和非凸Lp范数正则化的地震稀疏盲反褶积[J]. 石油地球物理勘探, 2016, 51(3): 428-433.
CAO J J. Seismic sparse blind deconvolution based on generalized Gaussian distribution and non-convex Lp norm regularization[J]. Oil Geophysical Prospecting, 2016, 51(3): 428-433.
[9]
DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582
[10]
孔德辉, 彭真明. 利用改进的在线字典学习估计时变子波[J]. 石油地球物理勘探, 2016, 51(5): 901-908.
KONG D H, PENG Z M. Time-varying wavelet estimation based on improved online dictionary learning[J]. Oil Geophysical Prospecting, 2016, 51(5): 901-908.
[11]
练秋生, 王小娜, 石保顺, 等. 基于多重解析字典学习和观测矩阵优化的压缩感知[J]. 计算机学报, 2015, 38(6): 1162-1171.
LIAN Q S, WANG X N, SHI B S, et al. Compressed sensing based on multiple analytical dictionary learning and observation matrix optimization[J]. Journal of Computer Science, 2015, 38(6): 1162-1171.
[12]
MAIRAL J, BACH F, PONCE J, et al."Online dictionary learning for sparse coding." International Conference on Machine Learning[C]//Proceedings of the 26th Annual International Conference on Machine Learning.Montreal: ICML, 2009: 689-696
[13]
BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM Journal on Imaging Sciences, 2009, 2(1): 183-202. DOI:10.1137/080716542
[14]
LECUN Y, BENGIO Y, HINTON G. Deep learning[J]. Nature, 2015, 521(7553): 436-444. DOI:10.1038/nature14539
[15]
SCHMIDHUBER J. Deep Learning in neural networks:An overview[J]. Neural Networks, 2015, 61(1): 85-117.
[16]
GRAVES A, MOHAMED A R, HINTON G.Speech recognition with deep recurrent neural networks[C]//2013 IEEE International Conference on Acoustics, Speech and Signal Processing.Vancouver: IEEE, 2013: 6645-6649 http://www.cs.toronto.edu/~fritz/absps/RNN13.pdf
[17]
GRAVES A.Long short-term memory[M].Supervised sequence labelling with recurrent neural networks.Berlin Heidelberg: Springer, 2012: 1735-1780
[18]
PINEDA F J. Generalization of back-propagation to rcurrent neural networks[J]. Physical Review Letters, 1987, 59(19): 2229-2232. DOI:10.1103/PhysRevLett.59.2229
[19]
WERBOS P J. Backpropagation through time:what it does and how to do it[J]. Proceedings of the IEEE, 1990, 78(10): 1550-1560. DOI:10.1109/5.58337
[20]
KASAC J, DEUR J, NOVAKOVIC B, et al. A conjugate gradient-based BPTT-like optimal control algorithm[J]. 2009 IEEE Multi-conference on Systems and Control, 2009, 861-866.