基于梯度投影法的全变差正则化全波形反演

doi:10.11720/wtyht.2022.1320

基于梯度投影法的全变差正则化全波形反演

姚含^,, 徐海

贵州省地质环境监测院,贵州贵阳 550001

Total variation regularized full waveform inversion based on gradient projection method

YAO Han^,, XU Hai

Guizhou Geological Environment Monitoring Institute,Guiyang 550001,China

责任编辑: 叶佩

收稿日期: 2021-06-22 修回日期: 2022-04-13

基金资助:

贵州高层次人才科研启动专项资金资助项目(0203001018040)

Received: 2021-06-22 Revised: 2022-04-13

作者简介 About authors

姚含(1972-),男,2006年毕业于贵州工业大学,主要研究方向包括水文地质、工程地质及环境地质。Email: 124121080@qq.com

摘要

为降低地震全波形反演的不适定性,常用方法是引入未知模型的先验信息,从而将反演问题正则化。但是,传统正则化方法在包含多个先验信息的情况下,仍然面临挑战。本文提出一种扩展的全波形反演公式,其中包含对模型的凸集约束。本文以慢度平方作为反演的模型参数,展示了如何在施加全变差约束的同时,施加边界约束令其保持在一个物理意义上的可行范围内。为验证本文所提算法的适用性,分别开展简单模型及国际标准地质模型数值实验研究,结果表明,全变差正则化的引入可以提高光滑背景模型下高速扰动体的重构效果。

关键词： 梯度投影法; 全变差正则化; 全波形反演

Abstract

To reduce the ill-posedness of seismic full waveform inversion,a common method is to introduce prior information to regularize the inversion problem.Traditional regularization methods still face challenges even when they contain multiple prior information.This study proposed an extended full waveform inversion formula,which includes the convex set constraints on models.Specifically,this study showed how to constrain the total variation of the slowness square while forcing the constraint to keep it within a physical reality range.To verify the applicability of the algorithm proposed in this study,numerical experiments on simple models and international standard geological models were carried out.The results show that the introduction of total variation regularization can improve the reconstruction of high-speed disturbances under smooth background models.

Keywords： gradient projection method; total variation regularization; full waveform inversion

PDF (3053KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

姚含, 徐海. 基于梯度投影法的全变差正则化全波形反演[J]. 物探与化探, 2022, 46(4): 977-981 doi:10.11720/wtyht.2022.1320

YAO Han, XU Hai. Total variation regularized full waveform inversion based on gradient projection method[J]. Geophysical and Geochemical Exploration, 2022, 46(4): 977-981 doi:10.11720/wtyht.2022.1320

0 引言

我国油气资源勘探从构造勘探阶段逐步走向岩性勘探阶段,勘探难度逐渐加大。地震勘探具有勘探深度大、精度高、数据量大等特点,是目前油气勘探最有效的方法,应用广泛。为进一步提高地震勘探精度,地震全波形反演方法迅速发展起来,成为国内外学者研究的热点。

全波形反演(full waveform inversion,FWI)是一种高精度、高分辨率的速度反演方法^[1-2]。FWI利用地震资料中的运动学和动力学信息,通过最小化模拟数据与实际观测数据之间的误差来重构地下介质模型^[3-4],能得到地下介质中更多细节的参数变化特征^[5]。由于实际观测中,子波未知,数据频带有限,且通常存在一定程度的噪声干扰,波场与介质参数之间呈现较强的非线性关系,FWI体现为高度欠定问题,对初始模型的精度要求较高,导致FWI在实际应用中达不到预期的效果。

为了减少FWI的不适定性,加入先验模型约束的正则化项,能使反演具有更好的稳定性^[6-7]。全变差(total variation,TV)约束方法是一种非光滑约束,由Rudin等^[8]在1992年提出,主要用于图像的去噪处理,能有效处理不连续性解的求解问题,重构图像中的一些间断点,保留边缘信息。近几年来,全变差约束已经应用在电法成像^[9]、图像去噪^[10-11]、偏移速度分析^[12-13]、混合数据偏移^[14]等领域。地下地层通常存在较多复杂的构造变化带及岩性突变带,导致地层物理参数在这些地方产生不连续变化,应用基于TV约束的FWI可以有效保留这些不连续特征^[15]。

本文采用梯度投影法^[16],构建基于全变差约束的全波形反演目标函数,实现多约束集下的声波方程全波形反演,逐步搜索到全局极值点,避免陷入局部极值,提高全波形速度反演的精度与速度。

1 方法原理

FWI是一个最优化问题,通过对观测波场与理论波场的残差求取极小值来构造目标函数进行迭代更新,最终重构地下介质的参数模型。当前常用的目标函数主要是通过数据残差的l₂范数来构建,即:

(1)

G (m) = \frac{1}{2} ‖ d_{c a l} (m) - d_{o b s} ‖^{2},

我们考虑模型参数m,它对应于速度平方的倒数,是一个实向量m∈R^N,其中N是空间离散化中的点数。对于每个源和频率,波场和观测数据分别用u_sv∈C^N和d_obs∈ $C^{N_{r}}$ 表示,其中,s=1,…,N_s指源,v=1,…,N_v指频率,N_r为接收器的数量。d_cal=Pu_sv,P是将波场投射到接收器位置的算子。

目标函数对m求一阶导数可以得到梯度矩阵,公式为:

(2)

\begin{array}{r} \nabla G (m) & = \frac{\partial G (m)}{\partial m} \\ = R e [{(\frac{\partial d_{c a l} (m)}{\partial m})}^{T} (d_{c a l} (m) - d_{o b s})^{H}] \\ = R e [J^{T} Δ d^{H}], \end{array}

其中,J为Frechet导数矩阵,上标T和H分别表示矩阵的转置和共轭。对目标函数求二阶导数可以得到近似Hessian矩阵,公式为:

(3)

H_{a} = R e [J^{T} J^{H}] 。

根据高斯牛顿法可以得到模型的更新量为:

(4)

Δ m = - H_{a}^{- 1} \nabla G {(m)}_{m_{0}} 。

利用梯度下降法^[17]最小化目标函数G,模型更新量可以写成:

(5)

Δ m = a r g \underset{Δ m}{m i n} Δ m^{T} \nabla G (m) + \frac{1}{2} Δ m^{T} H_{a} Δ m 。

2 全变差正则化

如果把m表示成N₁×N₂的速度模型,网格间距为h,用(i, j)代表各个网格点坐标,其中i=1,2,…, N₁,j=1,2,…, N₂,则在纵向和横向两个方向的梯度分别为:

(6)

\begin{array}{l} {(\nabla m)}_{i, j}^{z} = \{\begin{array}{l} \frac{m_{i + 1, j} - m_{i, j}}{h} & i f i < N_{1} \\ 0 & i f i = N_{1} \end{array} \\ {(\nabla m)}_{i, j}^{x} = \{\begin{array}{l} \frac{m_{i, j + 1} - m_{i, j}}{h} & i f j < N_{2} \\ 0 & i f j = N_{2} \end{array} \end{array}

全变差定义式如下:

(7)

\begin{array}{l} ‖ m ‖_{T V} = \sum_{i \leq N_{1}, j \leq N_{2}} ‖ {(\nabla m)}_{i, j} ‖ \\ = \frac{1}{h} \sum_{i j} \sqrt{(m_{i + 1, j} - m_{i, j})^{2} + (m_{i, j + 1} - m_{i, j})^{2}}, \end{array}

式(7)代表离散模型中每一点的离散梯度的l₂范数之和。假设Neumann边界条件使这些差值在边界处为零。我们可以通过定义一个差分算子D来更紧凑地表示‖m‖_TV,这样Dm就是离散梯度的一个连接,(Dm)_n表示在n索引位置对应的离散梯度向量,n=1,…,N,N=N₁×N₂,可以定义:

(8)

‖ m ‖_{T V} = \overset{N}{\sum_{n = 1}} ‖ {(D m)}_{n} ‖ 。

将全变差约束项应用于常规全波形反演算法中,可以得到基于全变差约束的目标函数,公式为:

(9)

G (m) = \frac{1}{2} ‖ d_{c a l} (m) - d_{o b s} ‖^{2} + λ ‖ m ‖_{T V},

对于式(9),如果添加约束m∈[B₁,B₂]和‖m‖_TV≤τ,则式(9)具有等价的优化形式:

(10)

\begin{array}{l} m i n G (m) = \frac{1}{2} ‖ d_{c a l} (m) - d_{o b s} ‖^{2} 约 束 于 \\ m^{n} + Δ m \in [B_{1}, B_{2}], ‖ m ‖_{T V} \leq τ, \end{array}

其中,τ是与λ有关的一个正常数,τ值越小,全变差约束项在目标函数中的权重作用越大。此时反问题可以表示为:

(11a)

m^{n + 1} = m^{n} + Δ m

(11b)

\begin{array}{l} Δ m = a r g \underset{Δ m}{m i n} Δ m^{T} \nabla G (m^{n}) + \frac{1}{2} Δ m^{T} H_{a} Δ m m^{n} + \\ Δ m \in [B_{1}, B_{2}], ‖ m ‖_{T V} \leq τ 。 \end{array}

根据改进的原始对偶混合梯度(PDHG)方法,求解式(11)的全变差约束问题相当于求解以下问题的鞍点(Δm,p):

(12)

\begin{array}{l} \underset{Δ m}{m i n} \underset{p}{m a x} [Δ m^{T} \nabla G (m^{n}) + \frac{1}{2} Δ m^{T} H_{a} Δ m + \\ p^{T} D (m^{n} + Δ m) - τ ‖ p ‖_{\infty, 2}], \end{array}

其中,p是拉格朗日乘子,mⁿ+Δm∈[B₁,B₂],‖·‖_∞,2为‖·‖_1,2的对偶模,相当于求最大值,而不是l₂模的和。

修改后的PDHG方法需要迭代:

(13)

\begin{array}{l} p^{k + 1} = a r g \underset{p}{m i n} τ ‖ p ‖_{\infty, 2} - p^{T} D (m^{n} + Δ m^{k}) + \\ \frac{1}{2 δ} ‖ p - p^{k} ‖^{2} \\ Δ m^{k + 1} = a r g \underset{Δ m}{m i n} Δ m^{T} \nabla G (m^{n}) + \frac{1}{2} Δ m^{T} H_{a} Δ m + \\ Δ m^{T} D^{T} (2 p^{k + 1} - p^{k}) + \frac{1}{2 α} ‖ Δ m - Δ m^{k} ‖^{2}, \\ m^{n} + Δ m \in [B_{1}, B_{2}] \end{array}

迭代公式可以更明确地写成:

(14)

\begin{array}{l} p^{k + 1} = p^{k} + δ D (m^{n} + Δ m^{k}) - \\ \prod_{‖ \cdot ‖_{1,2} \leq τ δ} [p^{k} + δ D (m^{n} + Δ m^{k})] \\ Δ m^{k + 1} = {(H_{a} + \frac{1}{α} I)}^{- 1} \cdot \\ [- \nabla G (m^{n}) + \frac{Δ m^{k}}{α} - D^{T} (2 p^{k + 1} - p^{k})], \end{array}

其中,δ和α为迭代步长,满足0<αδ<1/‖D^TD‖, $\prod_{‖ \cdot ‖_{1,2} \leq τ δ}$ (z)表示z在半径为τδ的球上‖·‖_1,2范数的正交投影。终止条件为:

(15)

m a x (\frac{‖ p^{k + 1} - p^{k} ‖}{‖ p^{k + 1} ‖}, \frac{‖ Δ m^{k + 1} - Δ m^{k} ‖}{‖ Δ m^{k + 1} ‖}) \leq 10^{- 5} 。

3 数值模拟

考虑一个2D合成实验,模拟区域网格点大小为200×200,网格宽度h为10 m。图1所示的合成速度模型有一个由较慢的平滑背景包围的恒定的高速区域。将该平滑背景作为反演的初始模型m⁰。类似于van Leeuwen和Herrmann^[18]中的例1,将N_s=34个源置于矩形高速异常体左侧,N_r=81个接收器置于矩形高速异常体右侧。源q_sv对应于主频为30 Hz的Ricker小波。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 真实速度模型

Fig.1 True velocity model

在5~35 Hz共31个不同频率上开展反演实验。本文考虑无噪声和微噪声两类数据,在有噪声的情况下,对每个频率指标v分别加入随机高斯噪声d_v,标准差为0.05‖d_v‖/ $\sqrt{N_{s} N_{r}}$ 。这可能不是一个真实的噪声模型,但它至少可以表明该方法对数据中的少量噪声是鲁棒的。

正则化参数τ有3种不同的选择:① τ_opt,即慢度平方的总变化量;② τ_large=1 000τ_opt, 此时τ_large足够大,全变差约束没有影响;③τ_small=0.875τ_opt,比理想的选择稍微小一点。注意,通过使用从式(5)开始的高斯牛顿步骤作为初始模型,式(11)中的凸子问题在τ_large情况下立即收敛。对于所有实验,惩罚参数λ固定为1。从5 Hz数据和6 Hz数据开始,使用计算出的m作为反演6 Hz和7 Hz数据的初始模型,按照由低频到高频的顺序,以此类推。每一个频带都进行50次外部迭代,以确保凸子问题收敛,达到式(15)所示终止条件时,迭代停止。

图2为开展的6次数值实验结果。在无噪声和有噪声的情况下,TV正则化均减少了反演模型中的振荡现象,并使得高速异常体得到更好的恢复。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 不同正则化参数τ在有噪声和无噪声情况下反演的速度模型

Fig.2 Inversed velocity model from noise free data and noisy data under different regularization parameter τ

从噪声数据中反演得到的不连续点的位置要比从无噪声数据中反演的更准确。这是因为噪声在接收器上造成了更大的不连续,从而加强了其他地方TV正则化的效果。τ_small实验证明了这一点,增加TV正则化的权重会导致较大不连续点的分辨率更高。

选取国际标准复杂地质模型开展全波形反演实验。真实速度模型如图3a所示,取自墨西哥湾典型复杂岩体的合成数据——SEAM模型。模拟区域网格点大小为267×384,网格宽度h为20 m。在靠近地面处放置总共72个源和191个接收器。初始模型采用图3b所示平滑模型。图3c及图3d为无噪声数据下的反演结果,图3e及图3f为噪声数据下的反演结果。图3d及图3f为选取τ_small正则化参数下的反演结果,图3c及图3e为选取τ_large正则化参数下的反演结果。由图可知,全变差正则化有助于消除伪影,主要是减少盐下的噪声。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 Seam模型的真实速度模型、初始速度模型,以及不同正则化参数τ在有噪声和无噪声情况下反演的速度模型

Fig.3 True and initial velocity model of Seam model, and inversed velocity model from noise free data and noisy data under different regularization parameter τ

4 结论及讨论

本文提出了一种在计算上可行的梯度投影算法,该算法可将附加凸约束条件下的全波形反演二次惩罚公式最小化。本文展示了在模型上添加全变差约束和边界约束时,全波形反演中凸子问题的求解过程。数值实验表明,在数据有限或含有噪声的情况下,TV正则化可以通过消除伪影和精确识别不连续点的位置,提高模型反演效果。

本文围绕的梯度投影法,其本质依赖于目标函数的一阶偏微分信息,这使得收敛速度受到一定制约。因此,理论上,我们可以采用二阶方法进行非线性优化,以探索收敛速度是否有提升。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Virieux

, Operto

An overview of full-waveform inversion in exploration geophysics

[J]. Geophysics, 2009, 74(6):WCC1-WCC26.

DOI:10.1190/1.3238367 URL [本文引用: 1]

[2]

Andrew

P V

, Malcolm

Elastic versus acoustic inversion for marine surveys

[J]. Geophysical Journal International, 2018, 215(1):1003-1021.

DOI:10.1093/gji/ggy303 URL [本文引用: 1]

[3]

Lailly

. The seismic inverse problem as a sequence of before-stack migrations[C]// Conference on Inverse Scattering: Theory and Applications,Expanded Abstracts, 1983:206-220.

[本文引用: 1]

[4]

Tarantola

Inversion of seismic reflection data in the acoustic approximation

[J]. Geophysics, 1984, 49(8):1259-1266.

DOI:10.1190/1.1441754 URL [本文引用: 1]

[5]

Mora

, Wu

Elastic versus acoustic inversion for marine surveys

[J]. Geophysical Journal International, 2018, 214(1):596-622.

DOI:10.1093/gji/ggy166 URL [本文引用: 1]

[6]

Lin

Y Z

, Huang

L J

Acoustic and elastic-waveform inversion using a modified total-variation regularization scheme

[J]. Geophysical Journal International, 2015, 200:489-502.

DOI:10.1093/gji/ggu393 URL [本文引用: 1]

[7]

, Liu

, Wu

, et al.

A high-order total-variation regularisation method for full-waveform inversion

[J]. Journal of Geophysics and Engineering, 2021, 18(2): 241-252.

DOI:10.1093/jge/gxab010 URL [本文引用: 1]

[8]

Rudin

, Osher

, Fatemi

Nonlinear total variation based noise removal algorithms

[J]. Physica D:Nonlinear Phenomena, 1992, 60:259-268.

DOI:10.1016/0167-2789(92)90242-F URL [本文引用: 1]

[9]

Chung

E T

, Chan

T F

, Tai

X C

Electrical impedance tomography using level set representation and total variational regularization

[J]. Journal of Computational Physics, 2005, 205:357-372.

DOI:10.1016/j.jcp.2004.11.022 URL [本文引用: 1]

[10]

屈勇, 曹俊兴, 朱海东, 等.

一种改进的全变分地震图像去噪技术

[J]. 石油学报, 2011, 32(5):815-819.

DOI:10.7623/syxb201105011 [本文引用: 1]

消除地震记录中的噪声是地震资料处理的重要环节,它对地震资料的后续处理和解释都有重要意义。在对传统全变分图像去噪模型原理分析的基础上,提出了一种加入扩散张量的改进全变分图像去噪方法。这种方法能够加强全变分模型中扩散项沿地层纹理方向的滤波作用,使得地层沿层理方向的连续性得到一致性增强；同时控制垂直于地层层理方向上的扩散以保护边缘信息。对川西某实际地震资料的处理结果表明,该方法能够有效地消除地震图像中的噪声,增强同相轴的连续性,提高地震资料的信噪比。

, Cao

J X

, Zhu

H D

, et al.

An improved total variation technique for seismic image denoising

[J]. Acta Petrolei Sinica, 2011, 32(5):815-819.

[本文引用: 1]

[11]

公成敏

全变分原理在地震数据去噪中的应用

[J]. 计算机与数字工程, 2014, 42(7):1271-1274.

[本文引用: 1]

Gong

C M

Application of variational principle in seismic data denoising

[J]. Computer & Digital Engineering, 2014, 42(7):1271-1274.

[本文引用: 1]

[12]

Anagaw

A Y

. Total variation and adjoint stat methods for seismic wavefield imaging[M]. Alberta: Physics Department of University of Alberta, 2009.

[本文引用: 1]

[13]

Anagaw

A Y

, Sacchi

M D

Edge-preserving seismic imaging using the total variation method

[J]. Joutnal of Geophysics and Engineering, 2012, 9(2):138-146.

[本文引用: 1]

[14]

卢昕婷, 韩立国, 张盼, 等.

基于全变分原理的多震源混合数据直接偏移方法

[J]. 地球物理学报, 2015, 58(9):3335-3345.

[本文引用: 1]

X T

, Han

L G

, Zhang

, et al.

Direct migration method of multi-source blended data based on total variation

[J]. Chinese Journal of Geophysics, 2015, 58(9):3335-3345.

[本文引用: 1]

[15]

Esser

, Guasch

, Leeuwen

T V

, et al.

Total variation regularization strategies in full-waveform inversion

[J]. SIAM Journal on Imaging sciences, 2018, 11(1):376-406.

DOI:10.1137/17M111328X URL [本文引用: 1]

[16]

Birgin

E G

, Martínez

J M

, Raydan

Nonmonotone spectral projected gradient methods on convex sets

[J]. SIAM Journal on Control and Optimization, 2000, 10(4):1196-1211.

[本文引用: 1]

[17]

Bertsekas

D P

. Nonlinear programming:2^nd edition[M]. Belmont: Athena Scientific, 1999.

[本文引用: 1]

[18]

van Leeuwen

, Herrmann

F J

Mitigating local minima in full-waveform inversion by expanding the search space

[J]. Geophysical Journal International, 2013, 195:661-667.

DOI:10.1093/gji/ggt258 URL [本文引用: 1]

An overview of full-waveform inversion in exploration geophysics

2009

... 全波形反演(full waveform inversion,FWI)是一种高精度、高分辨率的速度反演方法^[1-2].FWI利用地震资料中的运动学和动力学信息,通过最小化模拟数据与实际观测数据之间的误差来重构地下介质模型^[3-4],能得到地下介质中更多细节的参数变化特征^[5].由于实际观测中,子波未知,数据频带有限,且通常存在一定程度的噪声干扰,波场与介质参数之间呈现较强的非线性关系,FWI体现为高度欠定问题,对初始模型的精度要求较高,导致FWI在实际应用中达不到预期的效果. ...

Elastic versus acoustic inversion for marine surveys

2018

1983

Inversion of seismic reflection data in the acoustic approximation

1984

Elastic versus acoustic inversion for marine surveys

2018

Acoustic and elastic-waveform inversion using a modified total-variation regularization scheme

2015

... 为了减少FWI的不适定性,加入先验模型约束的正则化项,能使反演具有更好的稳定性^[6-7].全变差(total variation,TV)约束方法是一种非光滑约束,由Rudin等^[8]在1992年提出,主要用于图像的去噪处理,能有效处理不连续性解的求解问题,重构图像中的一些间断点,保留边缘信息.近几年来,全变差约束已经应用在电法成像^[9]、图像去噪^[10-11]、偏移速度分析^[12-13]、混合数据偏移^[14]等领域.地下地层通常存在较多复杂的构造变化带及岩性突变带,导致地层物理参数在这些地方产生不连续变化,应用基于TV约束的FWI可以有效保留这些不连续特征^[15]. ...

A high-order total-variation regularisation method for full-waveform inversion

2021

Nonlinear total variation based noise removal algorithms

1992

Electrical impedance tomography using level set representation and total variational regularization

2005

一种改进的全变分地震图像去噪技术

2011

一种改进的全变分地震图像去噪技术

2011

全变分原理在地震数据去噪中的应用

2014

全变分原理在地震数据去噪中的应用

2014

2009

Edge-preserving seismic imaging using the total variation method

2012

基于全变分原理的多震源混合数据直接偏移方法

2015

基于全变分原理的多震源混合数据直接偏移方法

2015

Total variation regularization strategies in full-waveform inversion

2018

Nonmonotone spectral projected gradient methods on convex sets

2000

... 本文采用梯度投影法^[16],构建基于全变差约束的全波形反演目标函数,实现多约束集下的声波方程全波形反演,逐步搜索到全局极值点,避免陷入局部极值,提高全波形速度反演的精度与速度. ...

1999

... 利用梯度下降法^[17]最小化目标函数G,模型更新量可以写成: ...

Mitigating local minima in full-waveform inversion by expanding the search space

2013

... 考虑一个2D合成实验,模拟区域网格点大小为200×200,网格宽度h为10 m.图1所示的合成速度模型有一个由较慢的平滑背景包围的恒定的高速区域.将该平滑背景作为反演的初始模型m⁰.类似于van Leeuwen和Herrmann^[18]中的例1,将N_s=34个源置于矩形高速异常体左侧,N_r=81个接收器置于矩形高速异常体右侧.源q_sv对应于主频为30 Hz的Ricker小波. ...

〈

〉