0 问题的提出

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

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

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

$ \min \limits_{\boldsymbol{x} \in \boldsymbol{Q}}F(\boldsymbol{x})=f(\boldsymbol{x})+\varphi(\boldsymbol{x}) $

(1)

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

$ p_{L}(\boldsymbol{y})=\mathop {\arg \;\;{\text{min}}}\limits_{x \in \boldsymbol{R}} \left[\frac{1}{2}\left\|\boldsymbol{x}-\boldsymbol{y}+\frac{1}{L} \nabla f(\boldsymbol{y})\right\|_{2}^{2}+\varphi(\boldsymbol{x})\right] $

(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的初始估计值L₀大于实际值时，整个迭代过程的L偏大，无法快速收敛^[11]；②当函数f(x)的局部曲率在初始迭代点附近较大，但在最优点附近减小时，如果L一直偏大，很可能导致“周波跳跃”而无法收敛。

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

1 FISTA方法原理

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

$ \|\boldsymbol{W x}-\boldsymbol{S}\|_{2}^{2}+\lambda\|\boldsymbol{x}\|_{1} $

(3)

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

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

2 自适应步长FISTA方法原理

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

$ t_{k+1}^{2}-t_{k+1}=t_{k}^{2} $

(4)

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

$ t_{k+1}^{2}-t_{k+1}=\frac{L_{k+1}}{L_{k}} t_{k}^{2} $

(5)

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

3 理论模型分析

为了验证改进算法的稀疏脉冲反褶积效果，首先合成一个理论模型(图 3)，分别利用FISTA算法和改进算法对理论模型进行稀疏脉冲反褶积。分两种情况进行计算：①初始输入常数L值过大，L=10L₀(图 4)；②初始输入常数L值过小，L=0.3L₀(图 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算法具有更好的抗噪能力。

4 实际数据处理

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

5 结束语

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