0 引言

准确的地震波速度在地震勘探中起着至关重要的作用，是地震偏移、岩性识别和储层预测成功的关键。作为目前最先进的高分辨率速度建模方法，全波形反演(FWI)已持续得到学术界和产业界的关注。通常在给定的初始速度模型基础上，FWI以合成地震记录与观测记录差异作为目标函数，计算该目标函数对模型参数的梯度，再结合优化算法迭代更新速度模型参数。FWI迭代时所用梯度通常由伴随状态方法计算得到，即通过正传波场导数与反传波场进行互相关获取^[1-2]。然而，当地震波在衰减介质中传播时，常规的FWI没有考虑地震波衰减，无法恢复可靠的模型速度。而梯度随深度衰减的黏滞声波FWI受到深部梯度变小的影响，收敛速度缓慢。本文致力于开发一种梯度补偿方法以提高黏滞声波FWI的收敛效率和反演结果的精度。

计算效率低一直制约着FWI的发展。目前加速全波形反演的策略有两类。一类加速策略是从优化迭代方法的角度出发，利用高阶优化算法获得更为优越的更新方向，如二阶优化算法中的高斯—牛顿法和截断牛顿法^[3-4]。但这些方法每次迭代至少需要四次波动方程正演计算更新方向，引入不容小觑的额外计算量。因此高阶优化方法务须寻求收敛速度和计算量之间的平衡。此类方法中采用线性搜索技术的L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno)优化算法是性价比最高的二阶优化算法^[5-7]。同时，Castellanos等^[8]指出，L-BFGS方法结合逆Hessian近似矩阵的附加信息预处理可以获得收敛更快的预条件逆Hessian矩阵，即为P-L-BFGS (Preconditioned LBFGS)算法。

另一类加速策略是通过预处理获得收敛更快的梯度，如考虑地震波本身的能量特征^[9-10]或结合其他信息如地层倾信息^[11]进行改进。介质的黏滞效应不仅反映在地震波的能量上^[12]，也会影响全波形反演目标函数的梯度。从这一角度出发，通过补偿衰减为黏滞声波FWI提供更适宜收敛的梯度不失为一种有效的解决策略。该策略的核心在于如何选择合适的衰减补偿方法。Wang等^[13]提出了一种用于衰减补偿的自适应稳定逆时偏移方法，其中稳定因子可以用关于补偿峰值的经验公式显式地表示。与传统的低通滤波稳定方法相比，自适应稳定策略可以高度提升反向时间外推方法的数值稳定性和抗噪性。Zhao等^[14]采取了在完全弹性介质和黏弹性介质中通过外推波场得到的两个算子实现补偿，但该方法受到了特定成像条件的限制。Sun等^[15-16]通过构建稳定的补偿算子，在不同条件下进行波场外推，但计算负担巨大。此外，Chen等^[17]将正则化方法隐式地嵌入波动方程正演模拟中，避免了显式滤波所增加的额外计算量。

描述黏滞效应的模型大体可以分为基于广义标准线性体(SLS)^[18-21]和基于常Q模型两类，后者更常用。基于常Q模型，学者们提出了一系列描述黏滞介质中地震波特性的波动方程^{[11-12, 22-25]}。Zhu等^[26]推导了具有解耦变分数阶拉普拉斯算子的常Q黏滞声波方程。该方程可以用来补偿振幅衰减和校正相位畸变，非常适用于衰减补偿的逆时偏移。Chen等^[27]对该变分数阶解耦方程进行了改进，优化了方程的表示方式，获得了近似常Q固定分数阶拉普拉斯算子黏滞波动方程。该方程保留了原有波方程的解耦性，大大降低了方程的数值计算复杂性，但并没有明显降低黏滞波场的数值模拟精度。在该方程的基础上进行黏滞声波方程的最小二乘逆时偏移可以稳定地补偿介质的黏滞性，获得高分辨率的成像结果^[28]。因此，本文基于该常分数阶拉普拉斯算子常Q黏滞声波方程研究Q补偿稳定化方案，以获得稳定波场补偿，进而获得精度高且收敛更快的全波形反演目标函数的梯度，进而提高黏滞声波方程全波形反演的收敛速度。本文进行速度反演所使用的Q模型是由其他方法获得的，假设已知而不参与反演。

本文首先给出常分数阶拉普拉斯算子黏滞波动方程，阐明该方程的数值解法；再回顾FWI理论，分析地层黏滞性对黏滞声波FWI的影响，说明衰减补偿对反演的重要性；然后阐述基于该波动方程的补偿方案，获得含补偿项的补偿黏滞声波方程，用于伴随波场的计算；最后应用数值实例验证基于该补偿方案计算的梯度能够加快黏滞声波FWI的收敛。

1 方法原理 1.1 常分数阶拉普拉斯算子黏滞声波方程

Zhu等^[26]提出的解耦分数阶拉普拉斯算子黏滞声波方程为

$ \left[\frac{1}{{c}^{2}}\frac{{\partial }^{2}}{\partial {t}^{2}}+\eta {\left(-{\nabla }^{2}\right)}^{\gamma +1}+\tau \frac{\partial }{\partial t}{\left(-{\nabla }^{2}\right)}^{\gamma +\frac{1}{2}}\right]p(\boldsymbol{x}, t)=0 $

(1)

式中：$ p(\boldsymbol{x}, t) $为 t时刻、空间 x处的声波波场；$ {\nabla }^{2} $为拉普拉斯算子；γ、c、τ、η定义为

$ \left\{\begin{array}{l}\gamma =\frac{1}{\mathrm{\pi }Q}\\ \tau ={c}_{0}^{2\gamma -1}{\omega }_{0}^{-2\gamma }\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }\gamma \right)\\ c={c}_{0}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{\pi }\gamma }{2}\\ \eta ={c}_{0}^{2\gamma }{\omega }_{0}^{-2\gamma }\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\pi }\gamma \right)\end{array}\right. $

(2)

其中：c为波速；γ为无量纲的分数阶数，与品质因子Q有关；c₀为参考频率$ {\omega }_{0} $对应的速度。式(1)将地下介质的黏滞效应对地震波相位畸变和振幅衰减分别通过$ \eta {\left(-{\nabla }^{2}\right)}^{\gamma +1} $和$ \tau \frac{\partial }{\partial t}{\left(-{\nabla }^{2}\right)}^{\gamma +1/2} $控制。在均匀介质条件下，将式(1)转换到波数域，对应的衰减声波波场$ P\left(k, t\right) $的解析解可表示为

$ P\left(k, t\right)={\mathrm{e}}^{-\alpha t}\left[{A}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta t\right)+{A}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta t\right)\right] $

(3)

式中：波数$ k=\sqrt{{k}_{x}^{2}+{k}_{z}^{2}} $，$ {k}_{x} $和$ {k}_{z} $分别为波数在x和z方向的分量；A₁、A₂为待定系数，由初始条件和边界条件确定；α、β为

$ \left\{\begin{array}{l}\alpha =\frac{1}{2}{\tau }_{\mathrm{c}}{\left|k\right|}^{2\gamma +1}\begin{array}{cc}& {\tau }_{\mathrm{c}}=\frac{1}{Q}{c}_{0}^{2\gamma +1}{\omega }_{0}^{-2\gamma +1}\end{array}\\ \beta ={\eta }_{\mathrm{c}}{\left|k\right|}^{\gamma +1}\begin{array}{cc}& \end{array}{\eta }_{\mathrm{c}}={c}_{0}^{\gamma +1}{\omega }_{0}^{-\gamma }\end{array}\right. $

(4)

由式(3)可以看出，α和β分别控制着波场的振幅和相位，且振幅随时间呈指数衰减。在实际情况中，Q和c随空间变化，故$ \gamma $也是空间变化的，振幅衰减和相位畸变的程度也会发生对应的变化。应用式(1)进行波场模拟，需要计算变分数阶拉普拉斯算子，这将大大降低计算效率。若采用分数阶的平均值转化为常分数阶拉普拉斯方程，又将降低模拟的精度。为此采用Chen等^[27]的简化方程进行数值模拟和全波形反演。该简化方程兼具可操作性和较高数值精度的优点，具体为

$ \begin{array}{l}\frac{1}{\lambda {c}_{0}^{2}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\pi }\gamma \right)}\frac{{\partial }^{2}p}{\partial {t}^{2}}=\left(1-\frac{2}{\epsilon \mathrm{\pi }Q}\right){\nabla }^{2}p-\\ \frac{2}{\epsilon \mathrm{\pi }Q}{\left(\frac{{c}_{0}}{{\omega }_{\mathrm{d}}}\right)}^{\epsilon }{\left(-{\nabla }^{2}\right)}^{1+0.5\epsilon }p○\frac{1}{Q{c}_{0}}\frac{\partial }{\partial t}G\left(p\right)\end{array} $

(5)

式中：ω_d是主角频率；λ=$ {\left({\omega }_{\mathrm{d}}/{\omega }_{0}\right)}^{2\gamma } $；“$ ○ $”取负号；$ G\left(p\right)={\left(-{\nabla }^{2}\right)}^{0.505}p $；ε是应用泰勒展开进行近似简化时引入的参量，经过测算一般选为1/16。

1.2 基于伴随状态法的目标函数梯度计算

全波形反演通常选用$ {\mathrm{L}}_{2} $范数作为目标函数，其梯度通过伴随状态法计算获得^{[10, 29]}。若先将式(5)简单表示为

$ F\left({c}_{0}, p\right)=0 $

(6)

式中F表示为黏滞声波波动方程系统，则全波形反演的梯度可以写为

$ g\left(\boldsymbol{x}\right)=\int \frac{\partial F\left({c}_{0}, p\right)}{\partial {c}_{0}}{p}_{\mathrm{a}\mathrm{d}\mathrm{j}}\left(\boldsymbol{x}, t\right)\mathrm{d}t $

(7)

伴随波场$ {p}_{\mathrm{a}\mathrm{d}\mathrm{j}} $由式(5)对应的伴随方程获得。

根据式(5)，式(7)可以写为

$ \begin{array}{l}g=\int {p}_{\mathrm{a}\mathrm{d}\mathrm{j}}\left[-\frac{2}{{\left(\frac{{\omega }_{\mathrm{d}}}{{\omega }_{0}}\right)}^{2/\left(\mathrm{\pi }Q\right)}\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\left(\frac{1}{2Q}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{1}{Q}\right){c}_{0}^{3}}\times \right.\\ \begin{array}{cc}& \frac{{\partial }^{2}p}{\partial {t}^{2}}+\frac{2}{\mathrm{\pi }Q}{\left(\frac{1}{{\omega }_{\mathrm{d}}}\right)}^{1/16}{\left({c}_{0}\right)}^{-15/16}{\left(-{\nabla }^{2}\right)}^{33/32}p\end{array}-\\ \begin{array}{cc}& \left.\frac{1}{Q{c}_{0}^{2}}\frac{\partial }{\begin{array}{c}\partial t\\ \end{array}}{\left(-{\nabla }^{2}\right)}^{0.505}p\right]\mathrm{d}t\end{array}\end{array} $

(8)

由式（7）和式（8）可知，梯度受震源波场和伴随波场的黏滞效应影响。若对波场实施衰减补偿，则由此计算的梯度能量也会增强。梯度迭代更新方法中，作为二阶优化算法的L-BFGS算法非常有效，因为它能直接计算梯度与逆Hessian的乘积，从而通过无矩阵递归算法直接确定模型更新量。加入预条件算子的P-L-BFGS方法相比L-BFGS方法能更快收敛，故本文采用P-L-BFGS算法迭代更新模型。

1.3 基于稳定化策略的稳定衰减补偿方法

由于补偿是衰减的逆过程，则式(3)对应的补偿形式应为

$ \widetilde{P}\left(k, t\right)={\mathrm{e}}^{\alpha t}\left[{A}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta t\right)+{A}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta t\right)\right] $

(9)

由上式可知，在补偿过程中地震波场的振幅呈现指数型增长，容易导致数值不稳定。为确保波场的稳定补偿，本文提出一种含稳定因子的Q补偿方法。引入稳定因子$ \sigma $的波场补偿解析表达式为

$ \widetilde{P}\left(k, t\right)=\frac{{\mathrm{e}}^{-\alpha t}}{{\mathrm{e}}^{-2\alpha t}+{\sigma }^{2}}\left[{A}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta t\right)+{A}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta t\right)\right] $

(10)

当$ \sigma $取0时，式(10)退化为式(9)，稳定化程度随着$ \sigma $的增大而增强。利用三角函数的周期性，式(10)可等价为

$ {\widetilde{P}}^{n+1}-{B}_{1}{\widetilde{P}}^{n}+{B}_{2}{\widetilde{P}}^{n-1}=0 $

(11)

式中：$ {\widetilde{P}}^{n} $表示$ n\mathrm{\Delta }t $时刻的波数域波场，其中$ \mathrm{\Delta }t $为时间步长；系数B₁、B₂为

$ \left\{\begin{array}{l}{B}_{1}=\frac{2{\mathrm{e}}^{-\alpha \mathrm{\Delta }t}\left({\mathrm{e}}^{-2\alpha t}+{\sigma }^{2}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta \mathrm{\Delta }t\right)}{{\mathrm{e}}^{-2\alpha \left(t+\mathrm{\Delta }t\right)}+{\sigma }^{2}}\\ {B}_{2}=\frac{{\mathrm{e}}^{-2\alpha t}+{\sigma }^{2}{\mathrm{e}}^{-2\alpha \mathrm{\Delta }t}}{{\mathrm{e}}^{-2\alpha \left(t+\mathrm{\Delta }t\right)}+{\sigma }^{2}}\end{array}\right. $

(12)

将式(11)整理成与震源波场传播类似的表达形式

$ \frac{{\tilde{P}}^{n+1}-2{\widetilde{P}}^{n}+{\widetilde{P}}^{n-1}}{\mathrm{\Delta }{t}^{2}}={C}_{1}{\widetilde{P}}^{n}+{C}_{2}\frac{{\widetilde{P}}^{n}-{\widetilde{P}}^{n-1}}{\mathrm{\Delta }t} $

(13)

根据式（12），可以获得C₁和C₂的表达式

$ \left\{\begin{array}{l}{C}_{1}=\frac{2{\mathrm{e}}^{-\alpha \mathrm{\Delta }t}\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta \mathrm{\Delta }t\right)-{\mathrm{e}}^{-2\alpha t}-1}{\mathrm{\Delta }{t}^{2}}\\ {C}_{2}=\left(1-{\mathrm{e}}^{-2\alpha \mathrm{\Delta }t}\right)\frac{{\mathrm{e}}^{-2\alpha t}-{\sigma }^{2}}{{\mathrm{e}}^{-2\alpha t}+{\mathrm{e}}^{-2\alpha \mathrm{\Delta }t}}\end{array}\right. $

(14)

对C₁利用泰勒近似$ \mathrm{c}\mathrm{o}\mathrm{s}\left(\beta \mathrm{\Delta }t\right)\approx 1-\frac{1}{2}{\beta }^{2}\mathrm{\Delta }{t}^{2} $，$ {\mathrm{e}}^{-\alpha \mathrm{\Delta }t}\approx 1 $，$ {\mathrm{e}}^{-2\alpha \mathrm{\Delta }t}\approx 1 $，对C₂利用泰勒近似$ {\mathrm{e}}^{-2\alpha \mathrm{\Delta }t}-1\approx -2\alpha \mathrm{\Delta }t $分别进行简化，可得

$ \left\{\begin{array}{l}{C}_{1}\approx -{\beta }^{2}\\ {C}_{2}\approx \frac{2\alpha \left({\mathrm{e}}^{-2\alpha t}-{\sigma }^{2}\right)}{{\mathrm{e}}^{-2\alpha \left(t+\mathrm{\Delta }t\right)}+{\sigma }^{2}}\end{array}\right. $

(15)

$ {C}_{2} $可进一步简化为

$ {C}_{2}=2\alpha \frac{{\mathrm{e}}^{-2\alpha t}-{\sigma }^{2}}{{\mathrm{e}}^{-2\alpha \left(t+\mathrm{\Delta }t\right)}+{\sigma }^{2}}={\tau }_{\mathrm{c}}C $

(16)

式中

$ C={\left|k\right|}^{2\gamma +1}\left(\frac{{\mathrm{e}}^{-{\tau }_{\mathrm{c}}{\left|k\right|}^{2\gamma +1}t}-{\sigma }^{2}}{{\mathrm{e}}^{-{\tau }_{\mathrm{c}}{\left|k\right|}^{2\gamma +1}\left(t+\mathrm{\Delta }t\right)}+{\sigma }^{2}}\right) $

(17)

将式(13)由波数域转换到空间域，即获得时空域的波动方程。同时利用

$ {\tau }_{\mathrm{c}}{k}^{2\gamma +1}=\frac{1}{Q}{c}_{0}^{2\gamma +1}{\omega }_{0}^{-2\gamma }{k}^{2\gamma +1}\approx \frac{\lambda {c}_{0}}{Q}{k}^{1.01} $

(18)

可进一步简化。由于 $ {\tau }_{\mathrm{c}} $随空间变化，为简便计算式(17)，取最大值$ {\tau }_{\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{x}} $替代在空间域变化的$ {\tau }_{\mathrm{c}} $。将最终简化的$ {C}_{1} $和$ {C}_{2} $代入式(13)，得到最终的含稳定因子的波场传播方程。同时参照Chen等^[27]文献中的第一种思路，可得到最终时空域补偿方程。该补偿方程除最后一项外，其他项与式(5)相同，即右端$ ○ $取正号且$ G\left(p\right)={C}_{\mathrm{s}}p $，其中$ {C}_{\mathrm{s}} $为C在空间域的对应响应。

补偿稳定效果取决于稳定因子$ {\sigma }^{2} $的选取。当稳定因子取0时，该方程退化为无稳定策略的补偿方程。稳定因子越大，对补偿过程中波场能量就有越强的约束作用，但如果选取过大，可能会损失波场信息。为适应不同黏滞性的衰减模型，通常选用$ {\sigma }^{2}=r\mathrm{e}\mathrm{x}\mathrm{p}\left(-\tau \mathrm{c}\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{\Delta }t\right) $，其中$ r $为常数，一般令$ r={10}^{-2} $。

2 数值模拟算例 2.1 稳定衰减补偿策略的有效性测试

为了验证本文衰减补偿方法的有效性，选用具有强衰减区的二维BP气囱模型(图 1)进行测试。该模型网格数为398×161，空间采样间隔为10 m。Q模型由速度模型根据经验公式Q $ =3.516{c}^{2.2} $生成，其中速度c的单位为km/s。该模型在上部的气囱区吸收非常强，当地震波穿过该区域时会快速衰减，导致下层的反射信号非常弱，应用常规方法在气囱区下方无法获得可靠的反演结果。

初始速度模型(图 1c)的顶层60 m设为准确速度，其下速度由准确模型光滑而得，只保留了速度变化趋势。以200 m间隔在模型顶部布设20炮，首炮位于模型左边界右侧90 m处。接收方式为双边可变观测系统，检波器的位置随震源位置移动而移动，检波点间距为10 m。设置最大炮检距为1800 m，以期望接收到一定角度范围内的深层反射信息。震源子波选用主频为10 Hz的Ricker子波，参考角频率为200 rad/s。为了充分记录波场信息及满足稳定性条件，设置记录的时间样点数为2501，采样间隔设置为0.8 ms。

选用$ {\sigma }^{2}={10}^{-2}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\tau }_{\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{\Delta }t\right) $为补偿时的稳定因子，比较补偿前、后的反传波场及其表征补偿效果。图 2为第10炮不同时刻补偿前、后反传波场快照对比，可见未经衰减补偿的波场振幅相对较弱，而经补偿后恢复了振幅信息，在高衰减区域明显增强，表明该衰减补偿算法能实现有效的振幅补偿，且补偿后的反传波场没有不稳定现象。

2.2 基于稳定梯度衰减补偿的反演测试

为了验证稳定Q补偿对黏滞声波FWI的有效性，利用BP气囱模型进行反演测试。图 3为第一次迭代的梯度。由梯度趋势可以看出，Q补偿后，梯度的整体趋势大体不变，但模型中浅层尤其是高衰减区域的梯度能量明显增强，反映的结构信息也更清晰。

采用炮点平方预处理方式，P-L-BFGS法迭代20、40和100次反演结果如图 4所示。可以看出，基于梯度Q补偿的全波形反演在初始阶段的迭代就可以恢复大体的速度模型，而利用未补偿的梯度无法在反演初期获得大体的速度趋势。虽然随着迭代次数的增加，梯度未补偿黏滞声波全波形反演效果也得到改善，但即便迭代至100次，仍然未能反演出准确的速度。对比x=1000、1500和3200 m的最终反演的速度纵向曲线(图 5)可以看出，基于梯度Q补偿的全波形反演不仅能获得更精确的层速度，而且较清晰地揭示了深部的大倾角倾斜地层，验证了本文方法在加快黏滞声波全波形反演的良好效果。

2.3 不同补偿稳定因子选取测试

稳定因子$ {\sigma }^{2}=r\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\tau }_{\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{\Delta }t\right) $中r分别取不同值时，第10炮0.4s时刻的反传波场如图 6所示。对比图 2c与图 6可以看出，选用的稳定因子都可以使补偿的反传波场保持稳定，同时补偿后的反传波场均较未补偿的反传波场能量有明显增强，大稳定因子比小稳定因子补偿后的波场能量弱。

图 7展示了由图 6所示的反传波场计算的第一次迭代的梯度。对比图 7与图 3可见：与小稳定因子相比，大稳定因子补偿的梯度数值较小，但依然大于未补偿波场计算的梯度。

反演的最大迭代次数设为100，如果目标函数停止下降也可提前终止迭代。r取不同值的反演结果如图 8所示。对比最终反演结果(图 4c与图 8)可以看出，由于衰减效应，梯度未补偿的反演结果存在明显的错误。相对而言，Q补偿的反演结果更优越。基于不同的稳定因子进行梯度 Q补偿的反演结果不尽相同。当补偿稳定因子过大时，如$ r={10}^{-1} $，补偿过程被抑制，获得的反演结果（图 8a）反而出现偏差，只迭代17次就陷入局部极小，无法继续迭代。当补偿稳定因子过小时，如$ r={10}^{-32} $，反传过程出现不稳定现象，终止于第66次迭代，获得极不准确的反演结果(图 8d)。当补偿稳定因子较小但处于合理范围时，既能保证反演顺利进行同时又保证补偿过程稳定，如$ r={10}^{-2} $、$ r={10}^{-4} $及$ r={10}^{-16} $，迭代100次的反演效果（图 4c右、图 8b和图 8c）均比未补偿时的结果更准确。

应用平均相对误差定量评价反演结果的精度

$ E=\frac{1}{N}\sum\limits_{i=1}^{N}\frac{\left|{c}_{i}-{c}_{i}^{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}\right|}{{c}_{i}^{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}} $

(19)

式中：$ {c}_{i} $、$ {c}_{i}^{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}} $分别为反演结果和真实模型第i个样点的速度值；N为速度模型的离散点数。

图 9为不同r迭代17次和66次的反演结果误差对比。从图 9a可以看到，当$ r={10}^{-1} $时，补偿稳定因子过大，抑制了波场振幅，获得的反演结果偏差大。相对于传统未补偿的反演结果，使用较小稳定因子的梯度补偿方法的反演结果的误差更小。由图 9b可以看出，当补偿稳定因子过小时，如$ r={10}^{-32} $，即便最终的反演误差低于传统未补偿的误差，但由于未能保证补偿过程的稳定，出现明显的数值不稳定现象，影响反演结果。对不同稳定因子的反演结果进行误差分析可以看出，随着稳定因子的增大，反演误差呈现出先降后升的现象。在全波形反演过程中进行梯度补偿时，应选用不抑制补偿、同时保证数值计算稳定的补偿稳定因子。

以图 9的误差曲线推算出最优的补偿稳定因子为$ r={10}^{-2} $，故进一步测试其对二维BP盐丘模型(图 10)的适用性。该模型网格数为625×120，空间采样间隔为10 m，其Q值同样由速度模型根据经验公式生成。模型的空间采样间隔为10 m。初始模型(图 10c)保留了顶部60 m以上的准确速度，其余由准确速度模型平滑获得。以100 m间隔在模型表面布设60炮，首个炮点距模型左边界170 m；接收方式为双边可变观测系统，为接收大入射角的反射信息，设置最大炮检距为1200 m，道间距为10 m。选用主频为3、5 Hz的Ricker子波进行正、反演，参考角频率为200 rad/s，记录的样点数为2001，时间间隔为0.8 ms。主频为3 Hz子波的正演数据的反演结果作为主频为5 Hz子波的正演数据反演的初始模型，分别反演15和10次。对比反演结果(图 11、图 12)，梯度Q补偿可以有效加速反演收敛。即便选用较低频的数据，经过梯度补偿后，FWI结果的盐丘速度精度有所提高，断层附近层间微弱的变化得到凸显，对于盐丘下方高衰减区域的表征也更加清晰，证实了本文方法对不同模型的适应性。

3 结论

本文基于常分数阶拉普拉斯算子黏滞声波方程提出了一种基于Q补偿的梯度预处理方法，能加速黏滞声波全波形反演，改善反演效果。该梯度通过Q补偿后的反传波场与正传波场互相关获得。针对特定常分数阶拉普拉斯算子黏滞波动方程的反传波场Q补偿，本文提出了一种有效的稳定补偿策略，可以确保计算过程的稳定，且提升反演效果。得到如下结论。

(1) 基于常分数阶拉普拉斯算子黏滞声波方程，通过在时间—波数域的解析解中引入稳定因子，以确保Q补偿反传波场计算的稳定，保证了梯度计算过程中的数值稳定。数值算例表明，未经稳定因子作用的直接补偿会使反传波场不稳定，所提方法能有效解决补偿过程中因高频成分指数放大而导致的数值不稳定问题。

(2) BP气囱模型和盐丘模型全波形反演算例表明，实施所提出的稳定的Q补偿梯度预处理方法能为黏滞声波全波形反演提供更有利于收敛的目标函数梯度，相比于未补偿的黏滞声波全波形反演的梯度，补偿的梯度能量整体有所提升，不仅能提高反演的收敛速度，还能提供更精确的反演结果，尤其是在高衰减区域。

(3) 在本文提出的Q补偿全波形反演中，稳定因子需要选在一定范围内才有提升效果。在满足波场稳定延拓的前提下，应该依据所反演区域的黏滞性($ {\tau }_{\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{x}} $)尽量选择合适的稳定因子以确保相位信息准确。本文试算结果显示，稳定因子可选用$ {\sigma }^{2}={10}^{-2}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\tau }_{\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{\Delta }t\right) $。