基于截断牛顿法的频率域全波形反演方法

引用本文

周斯琛, 李振春, 张敏, 张凯. 基于截断牛顿法的频率域全波形反演方法[J]. 物探与化探, 2017,41(1): 147-152.
ZHOU Si-Chen, LI Zhen-Chun, ZHANG Min, ZHANG Kai. Full waveform inversion in frequency domain using the truncated Newton method[J]. Geophysical and Geochemical Exploration, 2017,41(1): 147-152. 复制到剪切板

Doi:10.11720/wtyht.2017.1.23
Permissions

《物探与化探》编辑部所有

基于截断牛顿法的频率域全波形反演方法

周斯琛, 李振春, 张敏, 张凯

中国石油大学(华东) 地球科学与技术学院,山东青岛 266580

作者简介: 周斯琛(1993-),男,硕士在读,就读于中国石油大学(华东),研究方向为地球物理数据处理。

收稿日期: 2016-05-06

基金: 山东省自然科学基金项目(ZR2013DL012); 高等学校博士学科点专项科研基金项目(20120133120012); 国家自然科学基金项目 (41504100)

摘要

全波形反演方法可以视为大型非线性最小化问题。其中Hessian算子对反演结果有着重要的影响,传统的优化方法只能近似地表示Hessian算子,反演精度较低,收敛速度较慢,且对于反演目标照明不足的深部区域,往往出现参数无法聚焦的情况。而一种新的优化方法截断牛顿法,通过计算Hessian矩阵与已知向量乘积的形式,能够获得更精确的Hessian算子信息,从而解决以上问题。本文基于截断牛顿法在频率域实现全波形反演,通过模型试算表明,截断牛顿法相对于有限内存BFGS(Limited-memory Broyden-Fletcher-Goldfarb-Shanno,L-BFGS)法,能够得到更精确的反演结果,同时能提高收敛速度,尤其对于照明不足的深部区域,截断牛顿法有更明显的优势。

关键词: 全波形反演; 截断牛顿法; Hessian算子; 频率域

中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2017)01-0147-06

Full waveform inversion in frequency domain using the truncated Newton method

ZHOU Si-Chen, LI Zhen-Chun, ZHANG Min, ZHANG Kai

School of Geoscience and Technology,China University of Petroleum (East China),Qingdao 266580,China

Abstract

The full waveform inversion method can be regarded as a large nonlinear minimization problem.In this method,Hessian operator exerts significant influence on the inversion result.Nevertheless,traditional optimization methods can only express the Hessian operator approximately,which will lead to low inversion accuracy,slow convergence speed,and the result that the parameter cannot be focused,especially for deep inversion region target with poor illumination.In contrast,the truncated Newton method,a new optimization method,can obtain the information of Hessian operator more accurately by the computation of the product of Hessian matrix and a known vector,and thus can solve the problem mentioned above.Therefore,this paper achieves full waveform inversion in the frequency domain based on truncated Newton method.And the model test shows that the truncated Newton method has more precise inversion results and improves the efficiency of inversion,especially for deep area with insufficient illumination compared with the limited memory BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno,L-BFGS) method.

Keyword: full waveform inversion; truncated Newton method; Hessian operator; frequency domain

Show Figures

0 引言

全波形反演(full waveform inversion, FWI)能够利用地震记录精确估计地下物性参数, FWI过程是正演波场与观测波场达到最小化的过程^{[1, 2, 3, 4, 5, 6]}。FWI的思想最早由Lailly^[7]和Tarantala^[8]提出的, 并首先应用在时间域中。由于时间域FWI公式的复杂性, 其计算效率相对较低。Pratt在1990年将FWI的思想推广至频率域中^[9], 并通过分析证明与时间域FWI相比, 频率域FWI具有天然多尺度优势^[10], 且在模型较小的情况下使用多源反演具有较高的计算效率。

FWI可以看作大型非线性最小化问题, 其中一个关键环节便是优化求解问题。对于FWI来说, 通常使用局部优化方法求解, 主要分为梯度类方法和牛顿类方法。梯度类方法主要包括最速下降法与共轭梯度法; 牛顿类方法主要有高斯牛顿法, 拟牛顿法等。Tarantala^[8]与Pratt^{[11, 12]}分别在时间域和频率域中使用了最速下降法作为FWI的优化方法, 该方法虽然计算效率较高, 但收敛速度慢, 且存在折叠效应。为了进一步提高收敛效率, 利用当前梯度方向与上一次模型更新量的组合, 提出使用预条件的共轭梯度法进行优化求解。总体上梯度类方法在计算效率与存储消耗上虽具有优势, 但收敛速度都较慢。因此, 同时利用梯度信息与Hessian算子信息的牛顿类方法得到了发展, 牛顿类方法的难度在于对Hessian算子的直接求解, Operto提出高斯— 牛顿法^{[13, 14]}, 近似Hessian算子进行求解, 但该方法对计算量与存储要求很高。另一种优化方法为BFGS法, 避免了直接求解Hessian矩阵, 在迭代过程中利用梯度构建Hessian矩阵^[15]。但此方法消耗内存巨大, 不利于处理大规模数据。为改善BFGS方法, L-BFGS法通过最近几次的迭代信息估算得到曲率信息, 再直接求解近似的Hessian逆算子, 节省了内存消耗^[16]。

截断牛顿(truncated newton, TN)法与上述介绍的优化方法不同, TN法是利用精确的Hessian算子信息进行模型更新。但只需要计算Hessian矩阵与已知向量乘积的总体^{[17, 18]}, 不需要直接表示出Hessian算子, 该方法能够节省内存, 同时充分利用Hessian矩阵的信息提高收敛速度和反演精度。

文中结合频率域FWI与TN法的优势, 提出了基于TN法的频率域FWI, 提高了反演精度与收敛速度。文中首先讨论了频率域全波形反演的基本理论, 然后介绍了TN法基本原理与流程。最后将两者结合起来, 对简单异常体模型和Marmousi模型分别进行了试算。其中对异常体模型的试算, 来验证TN法频率域FWI的有效性; 对Marmousi模型的试算, 分别采用L-BFGS法与TN法两种不同的优化方法进行反演, 并对相应的反演结果与收敛速度进行分析比较, 得出TN法在频率域FWI中反演精度与收敛速度上具有明显优势。

1 频率域全波形反演基本理论

FWI是正演模拟得到的波场与观测波场之间差的最小化过程, 频率域的正演过程可以通过

$\begin{matrix} A (m) P = S \end{matrix}$ （1）

实现。式中:A(m)表示阻抗矩阵; m表示地下介质模型参数; P表示波场向量; S表示震源向量。

频率域FWI的目标函数C(m)可表示为最小二乘的形式:

$\begin{matrix} C (m) = \frac{1}{2} δ d^{T} δ d^{*}, \end{matrix}$ （2）

式中δ d=d_cal-d_obs, d_cal是模型正演得到的波场, d_obs表示为观测波场; T表示转置; * 表示复共轭。从数值计算的角度上考虑, FWI是一个大型非线性最小化问题, 由于其非线性复杂程度高, 因此采取局部优化方法求解。模型更新定义为

$\begin{matrix} m_{k + 1} = m_{k} + α_{k} p_{k}, \end{matrix}$ （3）

式中:m_k, α _k和p_k分别表示第k次迭代的模型参数、步长以及更新方向。其中的步长可以通过抛物线拟合法或者线性搜索法得到。目标函数C(m)用二阶泰勒级数近似为

$\begin{matrix} \begin{matrix} C (m_{k} + α_{k} p_{k}) = C (m_{k}) + α_{k} {p^{T}}_{k} \nabla C (m_{k}) + \\ \frac{1}{2} α_{k}^{2} {p^{T}}_{k} \nabla^{2} C (m_{k}) p_{k} + o (‖ α_{k} p_{k} ‖^{3}) 。 \end{matrix} \end{matrix}$ （4）

求取目标函数最小值需要依照式(3)的迭代过程, 更新方向的求取是个重要环节。为了获得牛顿更新方向, 对式(4)进行二阶近似得到

$\begin{matrix} \begin{matrix} C (m_{k} + α_{k} p_{k}) \approx C (m_{k}) + α_{k} {p^{T}}_{k} \nabla C (m_{k}) + \\ \frac{1}{2} α_{k}^{2} {p^{T}}_{k} \nabla^{2} C (m_{k}) p_{k} = F_{k} (p_{k}), \end{matrix} \end{matrix}$ （5）

其中函数F_k(p_k)表示C(m_k+α _kp_k)的二阶近似。而对于F_k(p_k)最小值点对应的变量p_k称为牛顿更新方向, 可以表示为

$\begin{matrix} p_{k} = - [\nabla^{2} C (m_{k} {)]}^{- 1} \nabla C (m_{k}) 。 \end{matrix}$ （6）

当目标函数与二阶近似方程F_k(p_k)十分接近的情况下, 牛顿更新方向是可靠的。式(6)中∇²C(m_k)即为Hessian算子, 可以用H(m_k)表示, 将式(6)转换形式为

$\begin{matrix} H (m_{k}) p_{k} = - \nabla C (m_{k}) 。 \end{matrix}$ （7）

使用牛顿法对目标函数优化的重点在于对Hessian算子或Hessian逆算子的求解, 由于计算量和内存消耗上的考虑, 直接求解Hessian算子或Hessian逆算子计算成本较高。这也是牛顿法不能直接应用于全波形反演中的主要原因。所以发展了高斯牛顿法, BFGS法以及L-BFGS法等, 都是通过近似Hessian算子或Hessian逆算子的途径, 进而求解牛顿更新方向。

2 截断牛顿法基本原理

上述牛顿类方法只是求取近似的Hessian算子或Hessian逆算子, 并没有充分利用Hessian算子的信息, 而TN法能够改善此问题。TN法包含内外两部分循环, 其中外部循环为模型更新迭代, 依照式(3)进行; 内部循环设置为一个线性反演问题, 一般采用共轭梯度法。在内部循环中采用共轭梯度法获取牛顿更新方向, 只需要在过程中计算Hessian矩阵与已知向量的乘积, 不需要存储和直接计算Hessian算子, 因此TN法大大节省了内存的消耗并提高了收敛速度。Hessian算子与已知向量乘积的形式可以表示为

$\begin{matrix} H (m_{k}) v = \frac{\nabla C (m_{k} + h v) - \nabla C (m_{k})}{h}, \end{matrix}$ （8）

式中:v为已知向量, h表示为某一趋于0的正实数。从H(m_k)v的数学表示可以看出, 在求取一次 H(m_k)v的过程中, 需要进行四次正演的计算, 无论是在计算量与计算效率上都是可以接受的。TN法基本流程可以用图1简单表示出来。

该截断牛顿法的内部循环采用的是共轭梯度形式, 其中需要设立一个循环终止条件, 根据Eisenstat和Walker的思想^[19], 终止条件设立是与目标函数二阶近似精确程度有关, 对于共轭梯度内部循环的终止条件可以设为

$\begin{matrix} ‖ H (m_{k}) Δ m_{k} + \nabla C (m_{k}) ‖ \leq η_{k} ‖\nabla C (m_{k}) ‖, \end{matrix}$ （9）

式中:η _k即为强迫项, η _k值的大小反映了式(5)二阶近似的精确程度。精确程度越高, η _k的值对应较小, 反之, 精确程度越低, η _k值相应较大。而η _k的取值可以定义为

$\begin{matrix} η_{k} = \frac{‖\nabla C (m_{k}) - \nabla C (m_{k - 1}) - α_{k - 1} H (m_{k - 1}) Δ m_{k - 1} ‖}{‖\nabla C (m_{k - 1}) ‖} 。 \end{matrix}$ （10）

	Figure Option View Download New Window
	图1 TN法算法流程

3 模型试算

针对以上研究内容, 首先使用一个简单的异常体模型验证TN法频率域全波形反演方法的有效性。然后使用Marmousi模型进一步试算, 分别使用TN法与L-BFGS法实现频率域FWI, 分别从反演精度与收敛速度两方面进行对比说明, 验证TN法在优化求解上的优越性。

3.1 异常体模型试算

真实模型网格个数为101× 101, 网格大小设为20 m× 20 m, 其中包含的异常体为一个网格个数为9× 9的正方形, 异常体速度设为1 800 m/s, 背景速度为1 500 m/s, 如图2a所示。反演时选用的初始模型即为没有异常体的背景模型, 震源函数为10 Hz主频的雷克子波, 设置25个炮点位于模型的顶端, 模型边界上各有50个接收点, 选用两组频率进行反演, 第一组:2、4、6、8、10 Hz; 第二组:10、12、14、16、18 Hz, 其中每个频率组最大迭代次数为20次。

通过两组频率的反演, 最终得到反演结果如图2b所示, 形态大小基本上与真实异常体模型完全一致。为了更准确地反映反演的结果, 抽取了x=1 km 处的垂直速度曲线, 就同位置处真实模型与反演结果作对比, 如图2c左所示, 真实模型与反演结果在异常体位置处的速度大小基本吻合, 且精度较高。证明基于TN法的频率域FWI方法的有效性。

	Figure Option View Download New Window
	图2 异常体模型试算 a— 真实模型; b— 反演结果; c— x=1 km处垂直速度曲线

3.2 Marmousi模型试算

文中使用的Marmousi模型和初始模型如图2a、b所示, 模型网格个数为681× 141, 网格大小设为25 m× 25 m, 初始模型通过真实模型平滑所得到, 其中有165个震源与660个接收点设置在模型的顶部。震源函数设置为主频10 Hz的雷克子波, 反演过程中根据经验选取了4组频率进行反演:第一组1、1.5、2.5、3、4 Hz; 第二组4、4.5、5.5、6.5、7 Hz; 第三组7、8.5、9、9.5、10 Hz; 第四组10、12、14、16、18 Hz每组频率最大迭代次数为20次。在反演中, 除使用的优化方法不同之外, 其他参数都一致, 所得反演结果如图3c、d所示。

	Figure Option View Download New Window
	图3 Marmousi模型试算 a— 真实模型; b— 初始模型; c— TN法反演结果; d— L-BFGS法反演结果

从两种方法反演的结果来看, 无论是TN法还是L-BFGS法, 对于模型基本构造形态的反演都比较精确, 在中浅层两种方法反演的效果比较好; 但是位于深层的构造, 尤其是在照明较差的区域, 可以明显看出TN法相对于L-BFGS法具有明显优势。为了更好地说明这个问题, 从两个反演结果中抽取x=5 km处和x=10 km处的垂直速度曲线, 并与同位置处真实模型的速度曲线进行对比, 如图4所示。从图中可以清晰地看出, 在深度2 km以内的区域, 无论是TN法还是L-BFGS法反演出的速度都十分准确。但随着深度的不断增大, 深层的照明程度较低, 此时TN法由于对牛顿更新方向更加精确的求解, 在照明程度较低的中深层也能反演出较好的结果, 而L-BFGS方法在中深层反演的结果相对较差。

	Figure Option View Download New Window
	图4 速度曲线对比 a— x=5 km处速度曲线; b— x=10 km处速度曲线

此外, 文中又针对TN法与L-BFGS法在收敛速度上进行对比, 在相同的迭代次数情况下, 对比其目标函数收敛的快慢, 如图5所示。

从归一化目标函数的收敛曲线可以清晰地观察两种方法的收敛快慢, TN法能利用较少的迭代次数达到更精确的结果, 也表示TN法在单次迭代更新中能够获取更准确的牛顿更新方向, 体现了TN法在收敛速度上的优势。

	Figure Option View Download New Window
	图5 归一化目标函数的收敛曲线

4 结论

文中详细介绍了基于TN法的频率域FWI基本原理, 并对模型进行了试算, 得出如下结论:

对于频率域FWI, 在相同反演条件下, 基于TN法的频率域FWI比L-BFGS法在速度反演的精度和收敛速度上存在明显优势。尤其对于深层照明不足的区域, 由于TN法利用的是更加精确的Hessian算子信息, 相对于L-BFGS法在反演结果上有一定程度上的改善。但在迭代的过程中, 由于TN法存在一个内循环, 导致其在计算量上相对于L-BFGS法有一定程度上的增加, 合理的分配内循环次数以及采用混合迭代的策略, 可以提高其计算效率, 这部分内容也值得进一步研究。总体而言, TN法在全波形反演中具有一定的应用前景。

致谢:感谢中国石油大学(华东)地震波传播与成像课题组的老师和师兄、师姐的大力指导, 感谢审稿人的指导意见。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Virieux J, Operto S. An overview of full-waveform inversion in exploration geophysics[J]. Geophysics, 2009, 74(6): WCC1-WCC26. [本文引用:1]
[2]	高凤霞. 频率域波动方程多参数全波形反演方法研究[D]. 长春: 吉林大学, 2014. [本文引用:1]
[3]	Plessix R E, Perkins C. Thematic Set: Full waveform inversion of a deep water ocean bottom seismometer dataset[J]. First Break, 2010, 28(4): 71-78. [本文引用:1]
[4]	Xue Z, Zhu H. Full waveform inversion with sparsity constraint in seislet domain[C]//85^th Annual International Meeting, SEG, Expand ed Abstracts, 2015. [本文引用:1]
[5]	Prieux V, Brossier R, Operto S, et al. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field, Part 1: imaging compressional wave speed, density and attenuation[J]. Geophysical Journal International, 2013, 194(3): 1640-1664. [本文引用:1]
[6]	李媛媛, 李振春, 张凯. 频率域多尺度弹性波全波形反演[J]. 石油物探, 2015, 54(3): 317-323. [本文引用:1]
[7]	Lailly P. The seismic inverse problem as a sequence of before stack migrations[C]//Philadelphia: Conference on inverse scattering: theory and application, Society for Industrial and Applied Mathematics, 1983: 206-220. [本文引用:1]
[8]	Tarantola A. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 1984, 49(8): 1259-1266. [本文引用:2]
[9]	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. [本文引用:1]
[10]	Vigh D, Starr E W. Comparisons for waveform inversion, time domain or frequency domain?[C]//78^th Annual International Meeting, SEG, Expand ed Abstracts, 2008. [本文引用:1]
[11]	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. [本文引用:1]
[12]	Pratt R G, Shipp R M. Seismic waveform inversion in the frequency domain, Part 2: Fault delineation in sediments using crosshole data[J]. Geophysics, 1999, 64(3): 902-914. [本文引用:1]
[13]	Operto S, Ravaut C, Improta L, et al. Quantitative imaging of complex structures from dense wide-aperture seismic data by multiscale traveltime and waveform inversions: a case study[J]. Geophysical prospecting, 2004, 52(6): 625-651. [本文引用:1]
[14]	Operto S, Gholami Y, Prieux V, et al. A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice[J]. The Leading Edge, 2013, 32(9): 1040-1054. [本文引用:1]
[15]	刘璐, 刘洪, 张衡, 等. 基于修正拟牛顿公式的全波形反演[J]. 地球物理学报, 2013, 56(7): 2447-2451. [本文引用:1]
[16]	Brossier R, Operto S, Virieux J. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[J]. Geophysics, 2009, 74(6): WCC105-WCC118. [本文引用:1]
[17]	Nash S G. A survey of truncated-Newton methods[J]. Journal of Computational and Applied Mathematics, 2000, 124(1): 45-59. [本文引用:1]
[18]	王义, 董良国. 基于截断牛顿法的 VTI 介质声波多参数全波形反演[J]. 地球物理学报, 2015, 58(8): 2873-2885. [本文引用:1]
[19]	Eisenstat S C, Walker H F. Choosing the forcing terms in an inexact Newton method[J]. SIAM Journal on Scientific Computing, 1996, 17(1): 16-32. [本文引用:1]

2009

0.0