基于L－BFGS算法和同时激发震源的频率多尺度全波形反演

张生强，刘春成，韩立国，杨小椿

1.吉林大学地球探测科学与技术学院，长春 130026

2.中国海洋石油研究总院，北京 100027

0 引言

在地震勘探的每一个环节中，速度均起着至关重要的作用。目前没有任何直接手段可以获得地下介质的速度分布，通常只能利用地表观测到的地震资料来确定地层速度。用于确定地层速度的基本方法有速度分析、射线层析成像和叠前（逆时）偏移等，但得到的结果均比较粗糙。全波形反演可以利用叠前地震记录，迭代反演出高精度的地下介质速度。它不仅能够显示地下介质的构造信息，还能显示地下介质的速度变化。因此，研究全波形反演方法具有极其广泛的意义和价值。

全波形反演技术经历了漫长的发展过程。20世纪80年代，Tarantola［1－2］提出了基于广义最小二乘反演理论的时间域全波形反演方法，对近20多年多维地震反演理论的发展产生了深远的影响。Marfurt［3］明确指出，对于含有大量炮点的波动方程模拟问题，使用频率域有限元或有限差分法是最有效的手段。为提高计算效率，20世纪80年代末90年代初，Pratt等［4－7］将这一思想与 Tarantola的全波形反演方法相结合，发展了频率域声波和弹性波全波形反演方法。频率域全波形反演中各个频率数据是相互独立的，可以并行计算，相对于时域全波形反演具有计算高效、数据选择灵活等优势。频率采样间隔选择具有一定的灵活性，早期的频率域反演通常使用等间隔频率采样。Song等［8］首先将低频数据波形反演结果作为相邻高频数据波形反演的初始模型以提高稳定性。该方法随后被广泛应用于全波形反演研究中，并被称为串行反演策略，如Sourbier等［9］应用该方法，有效地防止了求解过程陷入目标函数的局部极小值。Pratt［10］随后提出将有效频段内频率域数据由低频到高频分组，每组内各频率成分同时参与反演（即组内并行反演策略），将较低频分组反演结果作为相邻较高频率分组数据反演的初始模型（即组间串行反演策略），有利于压制噪音和反演假象。

全波形反演计算成本高、效率低、内存需求大，这些都是限制其发展的主要原因。笔者提出在频率多尺度全波形反演中将L－BFGS（limited memory BFGS）算法与同时激发震源技术相结合的方法来改善这一现状。同时激发震源技术可以大大地提高波形反演过程中正演模拟的计算效率。求解全波形反演问题的优化算法有很多，如最速下降法、共轭梯度法和拟牛顿法等，Nocedal［11］提出的 L－BFGS数值优化方法是求解优化问题最有效的拟牛顿方法之一。2007年Qin Zhiwei（Tony）［12］对共轭梯度法、BFGS算法以及2种有限内存优化算法进行了对比，证明了L－BFGS算法是一种占用内存小、收敛速度快、计算精度高的算法，适用于求解大规模地球物理优化问题。但是由于全波形反演是一个强非线性反问题，针对这个问题L－BFGS算法可能会收敛到局部极值，得不到全局最优解。为了克服这个缺点，笔者使用频率多尺度反演方法，即将有效频段内频率域数据由低频到高频分组，采取组内并行反演、组间串行的反演策略。这样可以有效地避免收敛于局部极值，通常都能够寻找到全局最优值。最后通过数值模拟试验对这种频率多尺度全波形反演方法及其抗噪能力和反演速度扰动能力进行了测试，收到了很好的成效。

1 方法原理

频率域波动方程全波形反演是利用波形，即波场的振幅、频率以及相位信息，将数据残差反传波场与初始模型的入射波场互相关，更新地下介质参数，进而对地下介质物性参数成像的方法。其具有揭示复杂地质背景下构造与岩性细节信息的潜力。基于波动方程的全波形反演本身是个强非线性反问题，解这类问题的经典方法是首先将问题线性化，再通过逐步迭代最终确定一个使理论计算值与观测记录较好拟合的反演结果。从数学的角度来讲，全波形反演可以对应地分为2个问题：大型稀疏矩阵线性方程组的求解以及非线性局部最优化问题。

1.1 频率域正演

频率域正演是频率域全波形反演的基础，正演的精度和效率在很大程度上决定了反演的精度和效率。本文研究了二阶常密度声波方程频域正演。目前，频率域全波形反演中使用的正演方法包括有限差分法、有限元法和有限体积法三大类，笔者采用有限差分法（9点有限差分格式）求解频域声波方程，即将时空域的微分方程转化为频率－空间域的大型稀疏矩阵线性方程组问题，并采用LU分解法求解。该方法的优点在于可以重复使用已分解的矩阵因子，提高了求解相同频率不同震源的正演计算效率。

1.1.1 Helmholtz方程

对时域二阶标量声波方程作关于时间的Fourier变换可得到频率域标量声波方程，即Helmholtz方程［13］：

其中：▽2为 Laplace算子；A＝－▽2－k2为Helmholtz算子；u为地震波场；k＝ω/v为波数，ω为角频率，v为介质速度；S为震源项。

采用LU分解法求解方程（1）即可求得地震波场u。地震波场u和模型参数m之间的关系是非线性的，可以通过算子D（D为m和u之间的非线性函数）来定义：

1.1.2 同时激发震源

全波形反演中，正演模拟最佳震源选择的一个重要标准是多炮模拟的有效性。2008年Berkhout［14］提出了混合激发采集，即不同空间位置的多个震源按照随机线性编码方式激发组成超级炮。笔者在全波形反演过程中选择了混合激发采集策略（同时激发震源技术，如图1），主要原因有二：1）提高声波有限差分算法多炮正演模拟的计算效率，使用同时激发震源技术来代替每一炮独立的正演，这样可以在一个炮记录里同时模拟多炮正演；2）可以用更多的炮来对地下进行照明。

图1 混合激发采集（同时激发震源）方法示意图Fig.1 Schematic diagram of blended acquisition（simultaneous sources）method

超炮道集的构建可由下式表达：

式中，P指包含多个单炮记录共频道集Pj的道集集合。混合采集采用随机线性编码，不同震源按照一定的随机时间延迟激发。而对于同时激发震源，也采用随机线性编码，只是不同震源按照相同的时间激发。

当使用单炮正演时，每炮正演耗时为tseq，共Ns炮，则正演总用时为Nstseq。若对这些震源采用随机线性编码，按照相同的时间激发，组成同时激发震源，此时正演耗时为tsim。在这2种情况下，由于tsim≪Nstseq，所以同时激发震源技术大大地提高了正演模拟的计算效率［15］。

1.2 全波形反演

对于地震勘探中的每一个震源－检波器，定义波场残差Δd为检波器位置处的观测波场dobs与正演模拟计算的波场dcal（m）之间的差值，即

式中：模型m表示地下离散计算域中的一些物理参数，本文中这些模型参数是指定义在每个数值网格节点上的P波速度，下文将用v来代替m表示速度模型。这时需要定义一个目标函数C（v）用来衡量波场观测值与计算值之差。全波形反演的目的就是通过迭代更新使C（v）达到最小，最终求出v。根据最小二乘原理，目标函数为

其中，＊表示共轭转置算子。

1.2.1 波恩近似及反演问题的线性化

根据波恩近似理论，更新后的速度模型v可以表示为速度模型初值v0与速度模型扰动Δv之和，即

在初始速度模型v0点附近对目标函数进行二阶Taylor－Lagrange展开：

其中：整数M表示速度模型v中的元素个数；O（v3）表示误差项。式（7）两端同时对速度模型参数vl求导并写成矩阵形式：

当目标函数的一阶导数为0时，则目标函数在v0附近达到极小值，同时得到速度模型扰动的表达式：

其中：目标函数的一阶导矩阵为梯度，表示在v0点处目标函数最速上升的方向；目标函数的二阶导矩阵为Hessian矩阵，表示目标函数在v0点附近的弯曲度。求解方程（9）的方法通常称为牛顿方法。地震波场和速度模型之间的关系是非线性的（u＝D（v）），在全波形反演中，反演需要迭代多次才能收敛到目标函数的极小值。

1.2.2 L－BFGS算法

考虑如下优化问题：

其中：函数C∈R，是实值函数，称C（v）为问题（10）的目标函数；集合Ω为问题的定义域。

求解最优化问题（10），通常采用迭代算法，其迭代格式如下：

其中：步长α（k）由线搜索策略确定；d（k）为搜索方向。在众多的优化方法中，Nocedal［11］提出的BFGS方法因其稳定的数值效果和快速收敛性而被公认为是求解优化问题（10）最有效的拟牛顿算法。该方法在构造搜索方向时只需要利用目标函数及其一阶导数的信息，避免了Hessian矩阵的计算，减少了计算量，并且保持超线性收敛的优点。它的基本思想是用某个给定的对称正定矩阵H（k）近似于Hessian矩阵的逆，即通过校正公式逐步修正矩阵H（k），使其越来越近似于［▽2C（v（k））］－1。该算法中Hessian矩阵的逆近似H（k）的修正公式［11］为

但BFGS方法的一个显著缺点是需要在内存中存储Hessian矩阵的逆近似，当需要求解的问题规模很大时，需要消耗很多的计算机内存，因此不适合求解大规模优化问题（例如地球物理问题）。针对该问题，Nocedal［11］又提出了一种计算高效的 L－BFGS方法，其主要思想是：在当前迭代点v（k）处，对给定的初始对称正定矩阵和非负整数m，利用前面m步的信息，对进行修正m次产生H（k）。与标准的BFGS方法不同，L－BFGS算法无需存储矩阵H（k），仅存储m＋1个向量组就能计算H（k＋1）。

利用式（12）导出L－BFGS算法第k＋1次迭代的修正矩阵H（k＋1）的表达式：

综上所述，求解优化问题（10）的L－BFGS算法步骤如下。

第1步，给定初始速度模型v0∈Ω、初始对称正定矩阵H0、非负整数m、误差限ε＞0，令k＝0。

第3步，使用反向传播线搜索方法，确定步长α（k）＞0满足下面的 Wolfe－Powell条件：

第4步，令m＝min｛k＋1，m｝，取由公式（13）确定

第5步，若k＝k＋1，转第2步。

在 L－BFGS算法步骤中，计算H（k）g（k）的一个高效公式见文献［11］。Liu和 Nocedal［16］证明，在适当的条件下，L－BFGS方法在求解一致凸函数极小化问题时具有全局收敛性，大量的数值试验结果表明，此算法非常适合求解大规模优化问题。在实际计算中，可以根据问题规模的大小以及计算机容许的内存，通过选择适当的m来控制所需的存储量。通常情况下，m的取值为3～20［16］。

1.2.3 频率多尺度串行反演方法

频率域全波形反演开始前需要确定选用数据的最低和最高频率。低频数据对于恢复大尺度宏观速度模型极为重要，因此反演最低频率通常都设置在数据有效频段内能使用的最低频率附近。高频数据则对恢复速度模型细节信息极为重要，最大频率通常由反演目标、有效数据频段、频散误差优化和硬件条件等因素综合确定。笔者使用的频率多尺度串行反演方法是将有效频段内频率域数据由低频到高频分组，每组内各频率成分同时参与反演，将较低频分组反演结果作为相邻较高频率分组数据反演的初始模型，有利于压制噪音和反演假象。

笔者在详细介绍了同时激发震源方法、L－BFGS算法、频率多尺度串行反演方法以及它们优点的基础上，提出了在频率多尺度全波形反演中将LBFGS数值优化算法与同时激发震源技术相结合的方法，以期能在计算效率和存储量方面有所改进。其实现流程是：首先需要在选定反演频率的初始频带下给定一个初始速度模型，使用同时激发震源作正演模拟，利用观测数据和正演模拟数据的差值构造目标函数；然后使用L－BFGS优化算法对目标函数进行优化，并求出目标函数取得极小值时对应的速度模型；接下来将该速度模型作为相邻较高频带数据反演的初始速度模型，使用上述相同的方法求出该频带对应的速度模型；这样反复进行频带循环计算，直到达到最高频带循环次数为止，最后得出的速度模型即为最终反演得到的速度模型。

2 数值模拟试验

为了测试本文提出的全波形反演方法的性能，分别对复杂Marmousi模型、高速楔形体模型及逆冲断层模型进行了速度反演、抗噪能力研究和反演速度扰动能力研究。在各模型的试算中选用了以下相同的参数：震源子波为雷克子波，主频12Hz；接收时间长度为2.4s；频率采样间隔为0.42Hz，反演了2.92～25.41Hz间的10个频带，共55个频率采样点，其中每个频带包含10个频率，每个频带最多迭代计算30次；单炮震源数等于相应模型横向网格点数，并利用这些单炮组合成20个同时激发震源；检波器个数为相应模型横向网格点数；检波器和震源均位于模型的最顶部。

引入拟合误差来定量刻画速度反演的正确程度，拟合误差W定义为

其中：N为模型总道数；M为每道采样点数；v′ij为速度反演结果；vij为理论速度值；i、j则分别为采样点号和道号。

2.1 Marmousi模型

通过对 Marmousi纵波速度模型（图2a）的反演，来验证并分析笔者提出的全波形反演方法。模型网格为128×384，网格间距24m。原始Marmousi速度模型大尺度平滑后的速度设定为初始速度模型（图2b）。对初始速度模型从低频到高频逐次反演10个频带，得到最终速度模型（图2c）。

图2 Marmousi模型全波形反演Fig.2 Full waveform inversion of Marmousi model

对于该模型，采用10个频带的单炮正演一次所需要的时间大约为85s，同时激发震源（384炮）正演一次所需时间大约为700s；因此，单炮震源正演完384炮所需的时间32 640s，要远远大于同时激发震源正演所需要的时间。所以，相对于常规单炮全波形反演，笔者的全波形反演方法在计算效率上得到了很大提高，并且在计算过程中可以明显地发现在对计算机内存的占用方面也得到了很大的改善。对比图2a和图2c可以看出，用笔者的全波形反演方法通过有限频带反演得到的速度模型与实际模型非常接近，反演效果非常好。利用式（14）计算最终反演结果与实际Marmousi模型的拟合误差较小（W＝0.095 9），定量地说明了反演得到的速度精度较高。另外，在研究过程中发现，初始频段的选取和初始速度模型的选取是一种矛盾的关系：即若选取的初始频段频率较高，反演方法对初始速度模型的要求就比较高，初始速度模型比较准确时才能得到比较好的反演结果；相反，若选取的初始频段频率较低，反演方法对初始速度模型的要求降低，比较差一点的初始速度模型也能得到比较好的反演结果。

2.2 高速楔形体模型

利用信噪比（SNR）对压制随机噪声的质量进行评价，定义式［17］如下：

式中：S0为原始模型正演地震记录；S1为含噪声地震记录或压制噪声后的地震记录。

为了检验全波形反演方法对噪声干扰的适应程度，同时也为了问题的说明方便，笔者利用简单高速楔形体对含噪声模型进行试算。实际模型是在均匀介质中包含一个高速的楔形体，均匀介质和高速楔形体的纵波速度分别为3 500m/s和4 000m/s（图3a），模型网格为160×320，网格间距20m。以3 500m/s的均匀背景速度场作为初始速度模型，将加噪声的地震记录数据用于笔者的全波形反演方法中，图3b为最终反演结果。图4a是对截取的第10个同时激发震源作用于原始模型的时间域正演地震记录加入高斯随机白噪生成的含噪地震记录，SNR为11.147 3dB。图4b为加噪数据最终反演结果的同时激发震源时间域正演地震记录。从图3可以看到，反演的楔形体很清楚，与实际速度模型非常接近。利用式（14）计算最终反演结果与实际速度模型的拟合误差较小（W＝0.024 1）；从而定量地说明了该反演方法反演得到的速度精度很高。由图4可以明显看出，随机噪声得到了压制。根据式（15）计算可知，信噪比提高为22.251 8dB。所以笔者提出的全波形反演方法具有一定的抗噪能力，对于提高地震资料信噪比也是比较有效的。

图3 高速楔形体模型全波形反演Fig.3 Full waveform inversion of high－speed wedge model

图4 高速楔形体模型同时激发震源时间域正演地震记录Fig.4 Simultaneous source time domain forward seismic records of high－speed wedge model

2.3 逆冲断层模型

在经典的波场传播研究中，通常将实际介质简化为均匀介质或层状均匀介质，常规的地震勘探工作即建立在该假设的基础上。然而实际上，小尺度范围内的速度和密度的扰动广泛存在于地壳之中，并通过向各个方向上散射入射波的能量，共同影响人工地震数据的质量。所以，通常建立的均匀介质模型并不符合实际地下情况，需要对这些小尺度的异常加以考虑。在地震勘探中，介质的随机扰动可以理解为产生随机噪声的一种来源，所以它像随机噪声一样，对勘探精度有很大影响。为了考虑这些细小但是无法忽略的异常，笔者使用全波形反演方法对只具有速度扰动的随机介质逆冲断层模型进行试算。自相关函数、尺度、方差和粗糙系数是随机介质的4个要素［18］。介质随机扰动可以分为3种情况：1）以随机介质为目标体；2）目标体被随机介质遮挡；3）上述两种情况的复合。笔者选用的是第三种情况。

逆冲断层纵波速度模型如图5a所示，模型网格为187×801，网格间距25m。在图5a模型上加入相关长度为200m×20m、方差为3%、粗糙系数为2的Von Karman型随机介质，生成只具有速度扰动的随机介质模型（如图5b所示，为方便对比，浅部450m以上没有加入速度扰动）。对图5b模型做大尺度平滑后得到初始速度模型（图5c），并利用笔者提出的全波形反演方法进行反演，得到最终速度模型（图5d）。从图5d可以看到，反演得到的最终速度模型很清楚，与具有速度扰动特性的实际速度模型（图5b）非常接近。利用式（14）计算最终反演结果与实际速度模型的拟合误差W＝0.036 0，也定量说明了反演得到的速度精度较高，同时验证了笔者所提出的全波形反演方法具有反演速度扰动（随机介质）能力。

3 结论与展望

1）模型试算的结果表明：基于L－BFGS算法的频率多尺度全波形反演可以通过有限的频带反演出精度很高的速度；同时也证明了笔者提出的将LBFGS数值优化算法与同时激发震源技术相结合的方法有效地减少了内存开销，提高了计算效率，适用于求解全波形反演问题。

2）分析了随机噪声对笔者提出的全波形反演方法的影响，通过模型试算可知，该方法具有一定的抗噪能力，表现出较好的可靠性、有效性和稳定性。

图5 逆冲断层模型全波形反演Fig.5 Full waveform inversion of thrust fault model

3）介质中的随机扰动对高分辨率地震勘探有较大影响，笔者尝试了用全波形反演方法反演随机介质，将随机场理论引入到全波形反演中，充分考虑了地下介质的随机特性，使反演结果更符合实际地质情况。反演算法成功地将介质的速度扰动特征反演出来，且分辨率较高。

4）全波形反演是构造高分辨率地层速度模型的有效手段，随着计算机软硬件技术的推进、地震模拟技术的进步和大规模优化计算技术的发展，全波形反演必将在不远的将来被广泛地应用于生产实际。

（References）：

［1］Tarantola A.Inversion of Seismic Reflection Data in the Acoustic Approximation［J］.Geophysics，1984，49（8）：1259－1266.

［2］Tarantola A.A Strategy for Nonlinear Elastic Inversion of Seismic Reflection Data［J］.Geophysics，1986，51（10）：1893－1903.

［3］Marfurt K J.Accuracy of Finite－Difference and Finite－Element Modeling of the Scalar and Elastic Wave E－quation［J］.Geophysics，1984，49（5）：533－549.

［4］Pratt R G.Frequency－Domain Elastic Wave Modeling by Finite Differences：A Tool for Crosshole Seismic Imaging［J］.Geophysics，1990，55（5）：626－632.

［5］Pratt R G，Worthington M H.Inverse Theory Applied to Multi－Source Cross－Hole Tomography：Part 1：Acoustic Wave－Equation Method［J］.Geophysical Prospecting，1990，38（3）：287－310.

［6］Pratt R G.Inverse Theory Applied to Multi－Source Cross－Hole Tomography：Part 2：Elastic Wave－Equation Method［J］.Geophysical Prospecting，1990，38（3）：311－329.

［7］Pratt R G，Shin C，Hicks G J.Gauss－Newton and Full Newton Methods in Frequency－Space Seismic Waveform Inversion［J］.Geophysical Journal International，1998，133（2）：341－362.

［8］Song Z M，Williamson P R，Pratt R G.Frequency－Domain Acoustic－Wave Modeling and Inversion of Crosshole Data：Part 2：Inversion Method，Synthetic Experiments and Real－Data Results［J］.Geophysics，1995，60（3）：796－809.

［9］Sourbier F，Operto S，Virieux J，et al.FWT2D：A Massively Parallel Program for Frequency－Domain Full－Waveform Tomography of Wide－Aperture Seismic Data：Part 1：Algorithm［J］.Computers and Geosciences，2009，35（3）：487－495.

［10］Pratt R G.Seismic Waveform Inversion in the Frequency Domain：Part 1：Theory and Verification in a Physical Scale Model［J］.Geophysics，1999，64（3）：888－901.

［11］Nocedal J.Updating Quasi－Newton Matrices with Limited Storage［J］.Mathematics of Computation，1980，35：773－782.

［12］Qin Z W（Tony）.The Relationships Between CG，BFGS，and Two Limited－Memory Algorithms［J］.Electronic Journal of Undergraduate Mathematics，2007，12：5－20.

［13］韩利，韩立国，李翔，等.二阶声波方程频域PML边界条件及频域变网格步长并行计算［J］.吉林大学学报：地球科学版，2011，41（4）：1226－1232.Han Li，Han Liguo，Li Xiang，et al.PML Boundary Conditions for Second－Order Acoustic Wave Equations and Variable Grid Parallel Computation in Frequency－Domain Modeling［J］.Journal of Jilin University：Earth Science Edition，2011，41（4）：1226－1232.

［14］Berkhout A J.Changing the Mindset in Seismic Data Acquisition［J］.The Leading Edge，2008，27（7）：924－938.

［15］Neelamani R，Krohn C E，Krebs J R，et al.Efficient Seismic Forward Modeling Using Simultaneous Random Sources and Sparsity［J］.Geophysics，2010，75（6）：WB15－WB27.

［16］Liu D C，Nocedal J.On the Limited Memory BFGS Method for Large Scale Optimization［J］.Mathematical Programming，1989，45（1）：503－528.

［17］张亚红.Focal变换及其在地震数据去噪和插值中的应用研究［D］.长春：吉林大学，2010.Zhang Yahong.Focal Transformation Research and Application in Seismic Data Denoising and Interpolation［D］.Changchun：Jilin University，2010.

［18］王恩利.基于各向异性的复合介质弹性波场模拟与特征分析［D］.长春：吉林大学，2008.Wang Enli.The Forward Modeling and Character A－nalysis for Complex Media Based on Anisotropic Media［D］.Changchun：Jilin University，2008.