0 引言

在地震数据采集过程中，为了得到更多的地下构造和岩性信息，通常采用多炮激发多检波器接收的观测系统。但是，实际数据处理的时候，有时不需要如此稠密的震源采样^[1-2]。所以，应该寻求合适的方法，在保证计算精度的情况下减少采用的震源数，即减少每次迭代的正演次数，减少计算量，提高计算效率。

众所周知，全波形反演计算量巨大。因为一次模型迭代更新至少需要三次正演，而正演计算量不容小觑^[3]。所以，全波形反演方法一度因为计算量问题而停滞不前。计算机技术的发展缓解了计算量问题，给全波形反演带来了曙光，使其进入了快速发展时期。至今30多年，全波形反演取得了许多突破性成果^[4-5]。但是，计算量大的问题仍然没有得到彻底解决，至今仍是研究人员需要突破的关键性问题。国内研究人员主要从算法的收敛性和GPU (众核处理器)技术来加快全波形反演的速度。例如，邓武斌等^[6]提出将L-BFGS (limited-memory Broyden-Fletcher-Goldfarb-Shanno)算法和SPGL1(spectral projected gradient for L1 minimization)算法联合应用于全波形反演中；苏超等^[7]使用GPU并行计算技术对数据进行处理；龙桂华等^[8]在反演过程中用Hessian矩阵的对角线元素来做梯度类方法的预条件算子，能够吸收高斯牛顿法的二次收敛优势, 使得算法具有较快的收敛速度；刘璐等^[9]提出利用一种新的拟牛顿公式对DFP (Davidon-Fletcher-Powell)和BFGS算法进行修正，改进后的BFGS算法在近似Hessian矩阵逆矩阵时，不仅考虑了梯度和模型信息，还加入了目标函数本身的信息，在保证反演精度的同时，明显提高了反演效率；成景旺等^[10]提出伪Hessian矩阵的预处理方法，该方法吸收了高斯牛顿法的二次收敛优势，在不增加计算量的前提下，加快收敛速度。国外学者为减少计算量也提出了很多有效的方法。Li等^[11]提出使用压缩感知方法来加快全波形反演的速度。震源编码技术通过同时激发震源成为典型的减少计算量的方法之一^[12]；但是，其在加快运算速度的同时引进了震源之间的串扰噪声，从而降低了反演的精度^[13]。为此，研究人员为了发挥编码技术的优点，发现编码技术的运用效果取决于超级炮中编码的震源数^[14]。所以，在应用此技术的时候要注意选择合适的震源数参与编码。为了克服串扰噪声，对合适的震源编码步骤，如随机编码方案^[15]等进行研究是必然的。Schuster等^[16]提出相编码技术来压制串扰噪声。全部震源参与编码较部分震源参与编码有效率，但它也是对串扰噪声最敏感的方法^[17]。值得注意的是，应用L2范数时，震源编码技术不适于拖拽检波器数据。Choi和Alkhalifah^[18-19]引进了卷积波场并发展了应用互相关目标函数的震源编码方法来反演海洋检波器数据。

本文是关于应用逐减随机震源采样法减少全波形反演计算量问题的研究。减少每次迭代中参与运算的震源数量就会增加全波形反演的效率。随机震源采样方法，即每次迭代时从全部震源中随机抽取一部分震源参与运算，每次迭代所用震源不同，获得的信息几乎包含全部震源参与运算时获得的信息，这样就在保证了计算精度的同时提高了计算效率。许多学者已成功将随机震源法应用于全波形反演中并取得了良好效果。2010年，Moldovean^[20]在3D模型中实现了时间域规则和随机震源采样方法；2012年，Ha等^[21]提出循环震源采样方法并应用于Laplace域声波方程；2014年，Wang等^[22]提出不规则震源采样方法并应用于基于VTI (vertical transverse isotropy)介质的声波波动方程。相较于震源编码技术，逐减随机震源采样法不会产生串扰噪声。但是，若是采用周期采样方法会产生信号混叠问题。所以，在选取采样方法的时候一般偏向于选择逐减随机震源采样法。因为它可以将混叠噪声转化为低级随机噪声^[23]。

本文提出了逐减随机震源采样法，并将其应用于频率域二维黏滞声波波动方程中；还将其与无记忆拟牛顿算法、L-BFGS算法结合，以证明无记忆拟牛顿算法和L-BFGS进行全波形反演过程中应用逐减随机震源采样法的可行性。

1 方法原理 1.1 逐减随机震源采样方法

本文介绍的逐减随机震源采样方法如图 1所示。即每次迭代时，只从所有震源中选取一定数量的震源，每次选取时，每个震源被选中的概率是一样的；但是，这是概率事件，或许会发生使用某些震源次数比另一些多的情况。随机震源采样方法加入了随机因素，使得到的地下构造和岩性信息具有随机性，即每次迭代会得到一些不同的信息，直至迭代完毕使用了几乎所有得到的地下信息，使得在降低计算量的同时亦没有失去太多地下信息，保证了反演精度几乎不变。

随着迭代次数的增加，逐减随机震源采样方法参与计算的震源数逐渐减少。此方法在模型迭代开始时参与计算的震源数目较多，这样可以获取更多信息，得到好的反演模型，为下一次迭代奠定好的初始模型基础。

1.2 频率域全波形反演

之所以选择在频率域进行基于逐减随机震源采样法的全波形反演，是因为频率域反演具有计算高效、频率选择灵活、易引入衰减因子、适于多炮计算等优点^[24]。

1.2.1 频率域黏滞声波方程正演

虽然对二维声波方程模拟波在介质中传播的研究已经很多，但是实践证明，二维黏滞声波方程能更好地模拟波的传播；因为其加入了品质因子，考虑了地层的黏滞吸收作用，更好地贴合了实际介质。所以本文研究的波动方程为频率域二维黏滞声波波动方程：

式中：(x, z)=ρ(x, z)ṽ²(x, z, ω)，是体积模量；，是复速度，v(x, z)是声波方程中的实速度；Q(x, z)是品质因子；ρ(x, z)是密度；P(x, z)是波场值；ω是圆频率；ζ(x, z)是震源项。

所用边界条件为PML (完全匹配层)吸收边界条件，正演模拟方法为交错网格有限差分法。对方程(1)进行差分离散后可得到简化后的形式：

式中：A是阻抗矩阵，取决于模拟频率和介质性质；P是波场向量；F是震源向量。求解矩阵方程(2)，可得到对应于圆频率ω的波场值。

1.2.2 基于逐减随机震源采样法的频率域全波形反演流程

本文在频率域进行了随机震源采样类方法的实现。在此给出频率域逐减随机震源采样法的全波形反演流程(图 2)。其中，逐减随机震源采样法应用于“计算目标函数对模型参数的梯度和海森矩阵”这一步。每个震源计算出来的梯度的加和即为最终得到的梯度值。每次迭代过程中，通过随机震源采样方法选取一部分震源参与运算，求得梯度值的一个近似，然后用此近似梯度值参与到后面的运算。

1.2.3 全波形反演原理

本文将全波形反演中的目标函数设置为

式中：C(m)是反演的目标函数；δP_i是第i炮模拟波场δP_i^cal和理论波场P_i^obs的差；S_{d, i}是第i炮与偏移距相关的权系数；T代表矩阵转置；*代表矩阵的复共轭；m为待反演参数。

目标函数关于反演参数的一阶导数如下：

本文使用拟牛顿类局部优化算法求取搜索方向：

式中：H是海森矩阵；Δm是模型扰动量。由于使用的是拟牛顿类优化算法，那么H^-1是通过最近几次迭代得到的模型参数值和梯度值构造而得，满足拟牛顿条件即可。本文所用拟牛顿方法为无记忆拟牛顿算法，Liu等^[25]首次将其应用于全波形反演中，其原理及公式可以参考文献[25]，在此不再赘述。

得到模型扰动量后，更新模型可通过下式得到：

式中，α_k是迭代步长，可通过多种方式得到，本文是采用抛物线拟合方式求取得到的。

2 模型试算

本文分别截取了Overthrust模型和Marmousi模型具有典型构造特征的部分作为试算模型，如图 3所示。截取的模型大小均为75×300。其中：Overthrust模型的纵横网格间距均20 m；Marmousi模型的纵横网格间距均为24 m。观测系统为：在z=60 m深度处布设300炮震源，Overthrust模型和Marmousi模型炮间隔分别为20 m和24 m；检波器布设在地表，共300个，两模型各检波器之间的距离分别为20 m和24 m；震源位置与检波器位置一一对应，分布在网格点上；两模型最大偏移距分别为6.0 km和7.2 km。本文所用随机方法为：在Overthrust模型中，参与计算的震源最多95个，最少20个，中间随着迭代次数增加而从95递减；在Marmousi模型中，参与计算的震源最多85个，最少10个，中间随着迭代次数增加而从85递减。经过大量实验发现，两模型各自选择震源的方式达到的效果最好。

本文在频率域进行反演。震源主频为15 Hz。反演频率从小到大依次为3.17，5.37，9.77，15.14，17.58，21.24，27.10，31.01 Hz，共8个频率；每个频率最大迭代次数为40次。迭代终止条件为反演速度参数值变化率小于10^-5。初始模型由理论模型通过高斯平滑得到，如图 4所示。

2.1 Overthrust模型

图 5为利用Overthrust模型、基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法反演得到的速度模型。比较图 5a、b可以看出，针对Overthust模型，两种方法反演得到的速度模型几乎一致，没有差别。比较模型中具有典型构造特征的部分(圈出来的部分)可以看出，两方法均反演出了断裂带和推覆构造这样的构造分布信息。并且，由图 5b可以看出，逐减随机震源采样法反演得到的速度模型没有引进随机串扰噪声。

为了更直观地比较两种方法的反演精度，从两种方法反演得到的Overthrust模型和理论Overthrust模型中抽取地表第40个点和第140个点处、即位于水平位置x=800 m处和x=2 800 m处的纵向反演曲线进行对比，结果如图 6所示。由图 6中圈出来的部分可以看出，两种方法对于模型上部的反演都很精确，对于底部速度突变点处，利用全部震源反演得出的速度值较精确些，但是相差不大，可以接受。

现在从数值上精确确定基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法的计算精度。本文利用反演得到的模型参数值与理论模型参数值的相对误差M_err来判定反演模型的精确度，其计算公式为

式中：m_cal是反演模型参数；m_true是理论模型参数；N是模型的总网格点数。

本文所用计算机运行内存为16 G，CPU为运行频率为2.40 GHz的Intel (R) Xeon (R)，所有计算均使用同一台计算机。通过对Overthrust模型进行试算，并结合式(7)，得到了基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法的总计算用时和相对误差(表 1)。

从表 1可以看出，逐减随机震源采样法频率域全波形反演在保证反演精度的情况下，能快速减少计算量，较之前的计算效率提高了2.54倍。

本文将无记忆拟牛顿算法与随机震源采样法结合并将其应用于Overthrust模型进行试算，取得了不错的反演结果。

2.2 Marmousi模型

图 7为利用Marmousi模型、基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法反演得到的速度模型。比较图 7a、b可以看出，针对Marmousi模型，两种方法反演得到的速度模型几乎一致。比较模型中具有典型构造特征的部分(圈出来的部分)可以看出，两方法均反演出了构造断裂带和速度突变点等信息。并且，由图 7b可以看出，逐减随机震源采样法反演得到的速度模型没有引进随机串扰噪声。

为了更直观地比较两方法的反演精度，从两种方法反演得到的Marmousi模型和理论Marmousi模型中抽取地表第40个点和第140个点处、即位于水平位置x=1 000 m处和x=3 360 m处的纵向反演曲线进行对比，结果如图 8所示。在图 8中，z=1 400~1 500 m处反演速度拟合的误差很大。这是因为，观测系统布设在地表，能量随深度逐减，使得处于地表的检波器记录到的深部信息量减少，反演时信息不足导致不能准确得出此处精确速度值；并且，此处速度强烈突变，不利于反演。除此之外，两种方法对于模型上部和下部的反演都很精确。

通过Marmousi模型试算，并结合式(7)，得到了基于全部震源的全波形反演方法和基于逐减随机震源采样的全波形反演方法的总计算用时和相对误差(表 2)。

从表 2可以看出，逐减随机震源采样法频率域全波形反演在保证反演精度的情况下，能快速减少计算量，较之前的计算效率提高了4.44倍。

在Marmousi模型试算中，本文使用的优化算法是L-BFGS方法，也取得了满意的反演结果。

综上，从两模型的实验结果可以看出，随机震源采样方法具有不可替代的计算优势。

3 结论和展望

1)模型试算的结果表明：基于逐减随机震源采样法的频率域全波形反演可以通过有限的频带反演出精度很高的速度值；此方法合理地减少了参与计算的震源数，即减少了正演次数，提高了计算效率，缓解了全波形反演计算量大的问题。

2)通过分析反演得到的速度图，可以看出该方法没有引进随机串扰噪声，即该方法具有一定的抗噪能力，表现出较好的可靠性、有效性和稳定性。

3)逐减随机震源采样法与无记忆拟牛顿算法、L-BFGS算法结合完成了全波形反演工作，证明该方法可以与各种反演优化算法结合，说明其稳定可靠；利用Overthrust模型和Marmousi模型进行试算，说明该方法适用模型不单一，再次证明了该方法的稳定可靠。

基于随机震源采样法的全波形反演用于2D模型所能提高的效率有限，预计若将其应用于3D模型计算效率会大幅度提高。本文所做研究只是为3D模型提供借鉴及经验。在将来的发展中，研究如何利用最少的震源获取最多的地下地质信息是今后努力的方向。