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

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

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

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

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

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

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

式中:W^T为W的转置矩阵; θ为阈值; 令p=x_k－ $\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);x_k+1、x_k为第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)可化简为:

ISTA迭代流程如图 1所示。

图 1(Fig. 1)

图 1 ISTA迭代流程示意

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

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

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

图 2(Fig. 2)

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

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

图 3(Fig. 3)

图 3 类RNN的改进ISTA

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

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

式中:E_t为中间变量。

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

其中, 公式(5)中的y相当于公式(4)中的输入x_t, S相当于隐藏层权值M, x_k相当于公式(4)中的隐藏层参数N_t, h_θ相当于激活函数。

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

图 4(Fig. 4)

图 4 正向算法流程

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

图 5(Fig. 5)

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

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

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

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

图 6(Fig. 6)

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

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

图 7(Fig. 7)

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

图 8(Fig. 8)

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

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

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

图 9(Fig. 9)

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

图 10(Fig. 10)

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

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

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

图 11(Fig. 11)

图 11 雁木西地区T3井Line 299线主要目的层段局部地震剖面

图 12(Fig. 12)

图 12 采用传统ISTA(a)和改进ISTA(b)进行反褶积处理的结果(Line299)

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

图 13(Fig. 13)

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

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

图 14(Fig. 14)

图 14 频谱对比

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

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