自适应步长FISTA算法稀疏脉冲反褶积

潘树林 闫 柯 李凌云 蒋从元 石林光

(①西南石油大学地球科学与技术学院，四川成都 610500；②中国石油西南油气田分公司,四川成都 610051；③胜利油田分公司物探研究院，山东东营 257000；④四川职业技术学院电子电气工程系,四川遂宁 629000)

0 问题的提出

稀疏脉冲反褶积是目前较为常用的提高地震资料分辨率的方法，其目的是利用带限的地震记录反演地下反射系数，从而得到高分辨率的地震数据[1-2]。稀疏脉冲反褶积的原理是在传统的最小二乘反褶积的基础上，加上一个稀疏先验约束项，通过求解L2范数与L1范数之和形式的目标函数得到高分辨率的反演结果[3-4]。

基追踪[5]方法利用Nesterov[6-8]提出的近端梯度的概念，使用二次函数近似表示目标函数，最终利用梯度法求解。该方法具有结构简单、实现容易、计算效率高等特点，是主流求解方法。

FISTA算法是Beck等[9]提出的基追踪方法，用于求解具有复合结构函数的最值

(1)

式中：F(x)为具有复合结构的目标函数，x为反射系数向量；f(x)为具有Lipschitz连续梯度的凸函数；φ(x)为在可行域Q上不可微的凸函数。根据Nesterov[6-8]提出的近端梯度法，式(1)对应的近端算子可以写作

(2)

式中x,y∈R，y为x的近似值。由式(2)可以看出，Nesterov[6-8]采用二次函数近似式(1)的f(x)，通过求解式(2)迭代求解式(1)。式(2)的常数L即为内部梯度的步长，它决定了FISTA算法的收敛度。通常无法准确给出L，FISTA算法采用线性搜索方法寻找最佳L[10]。然而，线性搜索只能使L朝着增大的方向搜索，严重影响了FISTA算法的收敛性。线性搜索算法的主要缺点为：①当L的初始估计值L0大于实际值时，整个迭代过程的L偏大，无法快速收敛[11]；②当函数f(x)的局部曲率在初始迭代点附近较大，但在最优点附近减小时，如果L一直偏大，很可能导致“周波跳跃”而无法收敛。

针对上述问题前人进行了研究[12-16]，但效果不好。为此，本文提出了一种改进的自适应步长FISTA方法，在没有外部参数调整的情况下仍能够求解，并利用实验验证了方法的可靠性。

1 FISTA方法原理

时间域矩阵形式的稀疏脉冲反褶积的目标函数为[17]

(3)

式中：S为地震记录向量；W为由地震子波组成的矩阵；λ为正则化参数变量。式(3)前半部分为最小二乘约束项，代表通过褶积模型生成的地震记录与实际地震记录的相似程度，后半部分为数据的稀疏先验项。

根据前文所述，式(3)可以利用FISTA算法求解，其求解过程如图1所示。由于采用线性搜索的方法，常数L只能沿着增大的方向变化，因此L的初始值决定了算法的收敛度，从而造成FISTA算法的不可控性。

图1 FISTA计算流程图

2 自适应步长FISTA方法原理

当常数L不变时，辅助序列tk+1与L无关，可以直接由下式确定[18]

(4)

为了能够满足理论上的收敛性，构造一个新的辅助序列，该辅助序列不仅只由tk和tk+1决定，还与前、后两次迭代的常数Lk和Lk+1有关，即

(5)

由此，结合FISTA算法和式(5)可以得到一种自适应步长的改进算法(下文简称改进算法)，算法流程如图2所示。改进算法在不同模型下均可取得满意的计算效果。通过对比FISTA算法流程(图1)和改进算法流程(图2)可知，在每次迭代中，当L不满足F[pLk(yk)]≤QLk[pLk(yk),yk]时： ①图1通过Lk=nLk-1线性增大常数L，而在下一次迭代中直接使用Lk。②图2则在每次迭代中，仍通过Lk=bLk-1线性增大常数L，但是在下一次迭代前，先通过Lk=aLk-1收缩L，然后利用收缩后的Lk进行后续迭代；与此同时，在每一次迭代过程中，利用前、后两次迭代常数L的信息更新辅助序列tk，使其在理论上可收敛。

图2 自适应步长FISTA方法流程图

其中01。a为收缩因子，在每一次迭代前收缩常数L，b为放大因子，在每一次迭代中放大常数L。参数a，b的取值对计算结果影响较小，经实验得a=0.7、b=2.0

3 理论模型分析

为了验证改进算法的稀疏脉冲反褶积效果，首先合成一个理论模型(图3)，分别利用FISTA算法和改进算法对理论模型进行稀疏脉冲反褶积。分两种情况进行计算：①初始输入常数L值过大，L=10L0(图4)；②初始输入常数L值过小，L=0.3L0(图5)。整个过程中取a=0.7、b=2.0。

由图4可见：当初始L过大时，改进算法反褶积效果优于传统基于FISTA算法(图4a)；当初始L过大时，传统方法无法调整L的大小，改进算法能够自适应地调整L(图4b)；在相同迭代次数的情况下，由于常数L的调整使改进算法收敛度远远高于传统方法(图4c)。由图5可见：当初始L过小时，改进算法与传统方法计算过程大致相同，因此改进算法的计算效果与传统算法的差异不大(图5a)；改进算法在相同迭代次数下，收敛度略高于传统算法(图5b、图5c)。

采用带随机噪声的数据进一步测试算法效果。在理论合成地震记录(图3c)中加入随机噪声，得到了一个低信噪比理论模型(图6)，对该资料分别使用过大和过小L值进行反演，获得不同方法反褶积结果(图7)。可见，基于传统FISTA的稀疏脉冲反褶积方法在低信噪比情况下无法得到准确的反褶积结果，由改进算法能得到较准确的反褶积结果。上述分析进一步证明了改进方法具有更好的收敛性，能在不同信噪比下得到理想的反演结果，较常规FISTA算法具有更好的抗噪能力。

图3 理论模型

图4 初始L过大时反褶积结果

图5 初始L过小时反褶积结果

图6 含噪声理论模型

图7 不同方法反褶积结果(SNR=3.85dB)

4 实际数据处理

选取吐哈盆地雁木西北部地震资料测试改进算法的适用性[19-20]。图8为雁木西T3井区Line 299测线地震剖面，可见地震分辨率较低，无法满足精细解释的需求。为此，分别采用基于传统FISTA算法的稀疏脉冲反褶积以及改进算法对图8数据进行高分辨率处理。由处理结果可见： 在相同的时间段内，改进算法反褶积结果(图9b)的同相轴数目较FISTA算法(图9a)多；经过基于FISTA算法的稀疏脉冲反褶积以及改进算法处理的地震资料的分辨率都得到提升，但是后者的频带更宽(图10)、细节刻画更清楚(图11c)，利于高精度构造解释。

图8 雁木西T3井区Line 299测线地震剖面

图9 反褶积剖面

图10 图9数据的频谱

5 结束语

本文提出了一种基于自适应步长FISTA算法的稀疏脉冲反褶积方法，该方法在FISTA算法的基础上，通过在每一次迭代之前适当减小常数L，然后利用线性搜索的方式寻找最优的常数L，以达到自适应调整L的目的。为了使算法达到理论收敛，通过结合前、后两次的L，对传统FISTA算法的辅助序列进行修改，最终使整套算法在理论上得以收敛。理论模型与实际地震资料的处理、分析结果表明，改进方法具有更好的收敛性，能在不同信噪比下得到理想的反演结果，较常规FISTA算法具有更强的抗噪能力。

图11 图9的局部放大