基于罚函数的全波形反演优化模型及其算法研究

殷雅倩 叶沐芊 张灵溪

摘 要：地震波反演技术在地质勘探中具有重要意义，该文在总结归纳拉格朗日乘子法和约化空间法的基础上，提出了基于罚函数的求解地震波反演方法的优化模型及其算法，利用共轭梯度法求取地震波反演中最小二乘问题的最优解。在实际模型中的数值实验表明，该模型和算法是行之有效的。

关键词：全波形反演 最小二乘法 罚函数 模型

中图分类号：P631.4 文献标识码：A 文章编号：1674-098X（2015）08（c）-0036-02

全波形反演技术在近年来不断受到很多学者的关注，其研究在石油、工程等诸多勘探问题中具有重要意义。该文对全波形反演问题求解中的拉格朗日乘子法[1]和约化空间法[2]做了简要介绍，建立了基于罚函数的反演优化模型，并且提出利用交替算法求解全波形反演的最小二乘问题的方法，最后，本文用Matlab语言编写出算法，在实际模型中验证了模型和算法的有效性。

20世纪60年代末，Backus等[3]提出地震波反演理论。受当时较为发达的医学层析成像技术的影响，地质学家从中受到启发，利用地震波对地下不同介质的反映，准确描绘地下地层结构。此后，Claerbout等[4]陆续提出了建立在波动方程基础上的地球反演理论框架，地球反演技术得以发展起来。

全波形反演是利用完整波场信息进行地震波反演的方法[5]，其优势在于可以精确刻画模型细节，使得反演效果更加出色。Tarantola等[6]首先在1982年提出了基于最小二乘法进行反演的时间域全波形反演，随后Pratt[7]又将其扩展至频率域。频率域全波形反演较之时间域全波形反演，拥有更快的计算速度和更高的计算精度。在求解频率域全波形反演时，该文在研究拉格朗日乘子法和约化空间法的基础上，建立了基于罚函数的反演优化模型，利用交替方向算法对其最小二乘问题进行求解。在数值实验中得到体现。

1 全波形反演模型及主要方法

频率域的全波形反演主要是基于Helmhotz波动方程的模型求解，通过对反射地震波的数据收集，进行反演的模拟，在这个过程中，要求使得模拟值与真实值之间的差值最小。也可简要描述为如下约束最小二乘问题：

（1）

其中，为波场与观测数据之间的转化算子，为实验次数，表示第次实验，为波场，代表实际观测数据，为地层模型，能够反映地震波在不同介质中传播的差距；为该公式的约束条件，是联系波场、波源信息和地层模型的模型方程，其中，是离散之后的Helmhotz算子，是圆频，是作用在波场上的离散的Laplacian算子，代表波源信息。

求解上述约束最小二乘问题的常用方法有两种，一种是Lagrange乘子法，其优化目标函数更改如下：

。

（2）

其中是Lagrange乘子，表示共轭转置。

该方法的优越性在于以Lagrange乘子的形式将约束条件增加到目标函数中，使得搜索区域扩大，减少局部最小值对优化问题的影响。同时，Lagrange乘子法能够有效避免每次迭代都要对波动方程进行精确求解的弊端，减少计算量，提高了运算效率。但是在实际问题中，往往对较大范围的区域进行反演，如果每次迭代更新都对的波场信息进行储存，可能会造成存储量太大。因此，Lagarange乘子法对于全波形反演并不完全适用。

另外一种方法则是先求解Helmhotz方程，得到，再将求得的代入目标函数得到关于的单变量函数：

。 （3）

这种方法称为约化空间法，将代入到目标函数并对求导，得到梯度：

。 （4）

满足共轭方程。

约化空间法与Lagrange乘子法相比，在每次更新时无需对波场信息进行存储，即存储量较小，但是在随后的梯度计算时，则需要对波动方程以及共轭方程进行精确求解。由于求解这两步的难度较大，当在处理大型问题时，计算代价就会较高。

2 基于罚函数的反演优化模型提出

由于Lagrange乘子法和约化空间法在频率域的全波形反演上均有优势，但也有各自的不足，本文将两者的优势综合，将约束条件引入到罚函数中[8]，建立了基于罚函数的频率域全波形反演优化模型。

罚函数的模型摒弃了伴随波场，减少了计算量和存储量。同时扩大了原有目标函数的搜索空间，有效的限制了出现局部最小点的情况。罚函数模型将波动方程作为惩罚项，令物理条件更加松弛。在对模型实际求解时，可以通过适当调节惩罚因子，从而平衡约束条件误差以及目标误差，有利于得到更为便捷有效的求解方法。引入罚函数之后，可将目标函数改写为：

。

（5）

首先设为固定值，将上式展开，目标函数是关于的函数，对求导，最小化，对每次实验，均有：

。

（6）

整理之后，波场可由下列公式求解得到：

。 （7）

其中。将上式改写成矩阵形式，则有：

。 （8）

接着将设为固定值，对求导最小化，则有：

。（9）

将代入上式，整理可得：

。（10）

由于上式中的的系数矩阵是一个对角阵，我们可以直接求出。同时，令，可用牛顿法求解的迭代格式：

。 （11）

由此可以写出相应的算法，并利用已知模型进行验证。

算法如下：

（1）设置的初始值；

（2）进行下列交替方向迭代；

（3）根据公式（7）求出；

（4）再根据公式（11）更新；

（5）不断迭代直到很小，退出循环。

3 数值实验

通过数值实验将算法应用到求解地震波反演问题中，从而驗证算法的有效性。如图1所示，建立一个均匀速度的初始地层模型，并在其中心放置一个正方形块体。在区域左端放置51个震源，在同层右端的相应位置放置51个接收器，多个结果叠加可以使模拟结果更加接近于真实值。将深度y方向和区域长度x方向分别划分出51条网格线，形成51*51的网格区域，设定震源频率为10 Hz，网格间距为20 m，网格边界条件取吸收边界条件。

所有计算均在一台CPU为2.39Hz，内存为4GB的计算机上运行；编程语言为MATLAB R2012b，其中机器精度为1.1×1016。

图2是10次迭代之后的反演图像，已经可以比较清楚的反映地层结构。

4 结语

该文在拉格朗日乘子法和约化空间法的基础上，提出了基于罚函数的反演优化模型，并且利用交替方向算法进行迭代；数值实验中，该模型和算法能够反演出较好的结果。但是如何改善算法的速度，以及如何提高算法的实用性，还有待进一步研究。

参考文献

[1] Haber E， Ascher U M， Oldenburg D. On optimization techniques for solving nonlinear inverse problems[J].Inverse problems， 2000， 16（5）：1263-1280.

[2] Pratt R G. Frequency-domain elastic wave modeling by finite differences： A tool for crosshole seismic imaging[J].Geophysics， 2012，55（5）：626-632.

[3] Backus G E，Gilbert F.The Resolving Power of Gross Earth Data[J].Geophysical Journal of the Royal Astronomical Society，1968，16（2）：169-205.

[4] Claerbout，J.F.Toward a unified theory of reflector mapping[J].Geophysics，1971（3）：467-481.

[5] 卞愛飞，於文辉，周华伟.频率域全波形反演方法研究进展[J].地球物理学进展，2010（3）：982-993.

[6] Tarantola A，Valette B.Generalized non-linear inverse problems solved using the least squares criterion[J]. Review of geophysics and space physics，1982，20（2）：219-232.

[7] Pratt G R， Shin C S， Hicks G J. Gaussnewton And Full Newton Methods In Frequencyspace Seismic Waveform. Inversion[J]. Geophysical Journal International，1998，133（2）：341-362.

[8] Tristan Van Leeuwen， Herrmann F J. Mitigating local minima in full-waveform inversion by expanding the search space[J]. Geophysical Journal International，2013，195（1）：661-667.