一种基于动态目标泛函的全波形反演多解性评估方法

黄建平, 崔 超

(1.中国石油大学地球科学与技术学院,山东青岛 266580;2.海洋国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛 266071)

一种基于动态目标泛函的全波形反演多解性评估方法

黄建平1,2, 崔 超1,2

(1.中国石油大学地球科学与技术学院,山东青岛 266580;2.海洋国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛 266071)

全波形反演是一种精度较高的地下速度建模方法,但其存在计算量巨大、反演多解性、初始模型依赖性强、适用条件苛刻等困难,生产实用化仍然较为困难。将基于动态目标泛函的全波形反演方法应用于反演解的讨论中,利用已有的反演解作为先验信息,通过动态调整目标泛函,在保证地下介质模型能够准确解释地震数据的前提下得到尽可能不同的反演结果。进一步根据数据约束强弱差异选取不同的反演策略,新方法能够对数据覆盖较弱区域的反演结果进行有效评估。在实现算法和反演流程的基础上,通过典型海底模型进行试算。结果表明，新方法能够在保证解的准确性前提下对全波形反演的特征给出合理分析,相对于传统方法计算效率也有了明显提高。

全波形反演; 多解性; 目标泛函; 可靠性评估

全波形反演作为目前最有发展潜力的速度场建模手段之一[1-2],在理论发展和实际应用中受到了广泛关注[1,3-9]。然而由于地球物理理论不完善[10-14]、野外数据质量差[15-17]、初始模型构建困难[18-20]等多方面因素的影响,全波形反演结果往往不够理想。对于全波形反演解的讨论和评估成为近年来的研究重点之一。前人提出了全局寻优法[21-26]、概率分析法[27]、“Bootstraping”法[28]等反演解的分析方法。由于地球物理问题具有数据量大、计算量大等特征[29],以及波形反演问题本身所具有的复杂性[30-31],传统方法往往需要极大的计算、存储量,并且适用条件苛刻[30],目前缺少一种快速高效的全波形反演问题解的评估方法。为此,笔者引入一种基于动态目标泛函的全波形反演解的评估方法[29],根据数据对模型敏感程度不同给出不同的求解策略,并对该方法的适用性及优缺点进行讨论。

1 方法原理

1.1 传统的全波形反演理论

全波形反演通过对目标泛函的优化达到求解地下最优模型的目的。本文中将目标泛函表示为χ(m),目标泛函可以表示为空间目标泛函核函数的积分形式,即

χ(m)=∬L(m)dxdt=〈L(m)〉.

(1)

式中,L为χ(m)的核函数;m为模型参数;x为空间坐标;t为时间采样点。利用链式法则,目标泛函关于模型的一阶导数可以转化为其对于数据的一阶导数与数据关于模型的一阶导数的乘积,即

(2)

式中,u为正演数据;δm为模型扰动量;δu=muδm为数据对于模型的敏感核函数。对于该核函数的求解,可利用有限差分法对模型的不同参数进行扰动,然而由于全波形反演问题中的参数数目众多,并且正演的计算量巨大,现阶段上式的有限差分解法应用较为困难。利用伴随状态法[32]求解该核函数是一种较为有效的策略,将波动方程表示为

B(u)=s.

(3)

式中,B为正演算子;s为震源项。对式(3)关于模型参数m求导有

(4)

(5)

(6)

(7)

(8)

(9)

(10)

同时,正演算子的初值条件转化为伴随算子的终值条件。式(10)中震源波场的伴随场为观测数据与正演数据残差的反传波场。利用上述结论,公式(2)可以简化为

(11)

(12)

由式(12)可见,目标泛函关于模型参数的导数项可以表示为震源正传波场与地震数据残差反传波场的零偏移互相关。因此目标泛函关于模型参数导数项的求解可以简化为两次正演过程,一次为震源波场的正传,另一次为检波器波场残差的反传。

全波形反演通过对模型参数的迭代更新求得最优解,其迭代公式可以表示为

mk+1=mk-αkHkmχ.

(13)

式中,k为迭代次数;Hk为对梯度的修正项,将模型参数m表示为列向量，Hk为一大规模矩阵,可采用共轭梯度法、L-BFGS方法等;αk为步长,可由抛物插值方法获得[7]。

1.2 基于动态目标泛函的全波形反演

传统的全波形反演理论往往只能给出一个反演解,即通过对目标泛函的逐步寻优获得能够解释数据的最优模型,但是由于地球物理反演问题本身所具有的特征,能有效解释同一观测数据的介质模型往往有多种,主要原因包括:

(1)地球物理反演问题本身是一个欠定问题,由于观测数据的有限性和对地下介质先验知识的局限性,能够解释同一观测数据的模型往往有多个。

(2)全波形反演问题是一个强的非线性问题,目标泛函的性态较为复杂[20],存在较多的局部极小,对于全局极小的查找极为困难。因此全波形反演的反演结果往往是众多局部极值之一。

(3)由于观测系统照明的限制,地下介质存在对于观测数据没有影响或者影响极小的区域,因此目标泛函关于模型的导数矩阵在这些位置为0或者接近于0,即模型参数在该区域取任意值都可以对数据做较好地解释,全波形反演方法对此模型参数不具有反演能力或者反演能力很微弱。

Rawlinson等[29]通过修改目标泛函,将已有的反演结果作为先验信息,在传统目标泛函的基础上添加动态正则化项,在保证模型能够有效解释观测数据的基础上,尽量保证每次反演所得到的反演解不同。通过求解最后多组最优解的均值与均方差,从而对全波形反演多解性问题进行分析,并对反演结果的可靠性进行定量评估。修改后的目标泛函可以表示为

χ1=χ+S(m).

(14)

其中

式中,S(m)为动态正则化项;Mj为第j个已有的反演解;λ、p、ε为正则化项中人为控制参数。对于上述目标泛函关于模型参数求导可得

(15)

其中

由公式(15)可见,修改后的目标泛函在已有的反演结果位置处生成局部极大值,能够保证后续的反演结果与已有反演解不同。

1.3 目标泛函的性态分析与反演策略

图1给出了一个一维不同参数情况下的动态正则化项性态,令M=50,m在1～100变化。分析可见,λ主要决定动态正则化项的分布范围,取值过大则动态正则化项只能在已有解很小的邻域内进行约束,导致不同反演解之间十分接近;而取值过小会使动态正则化项的作用范围过大,影响目标泛函整体性态,甚至导致反演失败。p的作用与λ相似,在实际应用中往往取值为1。ε主要决定动态正则化项在已有解位置处的峰值大小,其取值越大对应的峰值越小。本文中每次反演均从初始模型开始,并令ε=0。总的来说,选取参数的原则为:正则化项的作用既要具有一定区域性,又能够有效改变目标泛函在已有解位置处的性态。

图1 不同参数取值情况下的动态正则化项的分布性态Fig.1 Regularization character of dynamic misfit function with different parameters

需要着重指出的是:在数据对模型变化不敏感的区域,即反演前后模型变化极小的区域,正则化项在初始模型位置生成一极大值,修正后的目标泛函在该位置处关于模型参数的导数始终为0,因此不能生成多个最优解。然而实际上模型参数在该区域取任意值对观测数据都可以作较好地解释,为了能够有效解决这一问题,在生成第一个最优解时查找出反演前后模型参数变化小于某一阈值的区域,认为数据在该区域对于模型的敏感性极小。为保证修正后的目标泛函在该区域关于模型参数的导数项不为0,可在反演过程中对Mj在该区域添加一定量的随机扰动,以保证每次所得到的最优解不同。

1.4 基于动态目标泛函的全波形反演流程

本文中利用已有的反演解作为先验信息，通过动态地修改目标泛函以获得尽可能不同的反演解。其反演的流程如图2所示。其中k表示内部循环的迭代次数(式(13))，j表示不同的外部循环次数。将终止条件(1)表示为目标泛函小于给定阈值或反演迭代次数k达到给定上限。实际上，每次内部循环都是一次常规意义的全波形反演过程，在每次内部循环结束时，返回一个反演结果Mj，将其作为已知反演解代入动态目标泛函，然后根据判断条件(2)决定是否继续生成反演解。由于每次反演都利用初始模型m0作为输入，因此内部循环开始时令k=0。将终止条件(2)表示为所得解的个数达到给定上限。

图2 基于动态目标泛函的全波形反演方法流程Fig.2 Flowchart of full waveform inversion by using dynamic function

2 模型测试

为了测试本文中方法的适用性,采用如图3(a)所示的海底地下模型。其中模型纵横向尺寸为4 km×8 km, 每100 m做一次正演模拟,共80炮,在海水表面从0开始每20 m放置一个检波器。地震子波采用主频为15 Hz的雷克子波,反演过程中假定震源子波和水层速度已知。对水层以下采用三角平滑方式得到反演的初始模型(图3(b))。

图4为采用传统的全波形反演方法得到的第50、100、250次迭代的反演,结果可见:由于给出的初始模型较为准确,随着迭代次数的增加,海底模型细节的刻画逐渐趋于明显,高波数成分逐渐被恢复。图5抽取了位于2 km位置处的纵向速度剖面对比,可见随着迭代次数的增加,地下介质分界面准确归位,尤其是位于2～2.2 km位置处的低速体被有效反演。但是由于观测系统的限制,全波形反演对模型的右下侧区域反演能力不足,反演结果较为模糊。而在模型的左侧,由于深层倾斜高速体的强反射作用受到较强照明,因而反演结果较好。在深层高速基地区域由于不能接收到有效的反射信息,模型更新较弱。

图3 反演所用的准确速度场与初始速度场Fig.3 True model and initial model used by inversion

图4 常规全波形反演不同迭代次数后的反演结果Fig.4 Inversion results after traditional full waveform inversion for different iteration times

图5 由图3和图4中模型中抽取2 km位置处的单道对比Fig.5 Single vertical velocity profile comparision from Fig.3 and Fig.4 at 2 km

为了能够定量地衡量全波形反演对不同模型区域的反演能力,并对反演模型可靠性进行评估,本文中采用基于动态目标泛函的全波形反演方法进行反演。共得到10个最优解,其中参数的设置为:λ=7×10-7,p=1,ε=0。由于动态目标泛函反演的目的在于对反演的多解性和结果的可靠性做评估,并非得到十分准确的反演结果,且基于较少次数的迭代就可以有效地反映出全波形反演对模型的更新能力。为了有效的降低计算量,每次反演过程只迭代20次。对反演结果进行加权平均得到的平均模型如图6所示,可见加权平均后的结果对模型仍具有一定的恢复能力,浅层层位更为清晰。对于10次基于动态目标泛函所得到的反演结果做均方差分析,如图7所示。浅层数据约束较强的中间区域(区域A)多次反演结果的均方差较小。表明全波形反演在数据覆盖强的区域具有稳定的反演能力,所得到的反演结果可靠性也比较高,而在深层基底(区域B)及模型右下侧数据照明较弱区域(区域C)均方差较大,说明不同反演次数得到的反演结果相差较大。观察图8可见不同反演结果(Mj)对应的整体数据残差变化较小,说明深层基底(B)以及数据照明较弱区域(C)模型参数的变化对数据影响不大,全波形反演在该区域反演能力不足且反演结果不可靠。对于这一现象的具体解释为:对于深层基底(B),由于不能接收穿过该区域的反射波,因此该区域的速度场变化对于数据几乎不产生影响;对于模型右下侧照明较弱区域(C),由于缺少较强的反射界面以及观测系统的限制,该区域速度的变化对于目标泛函的影响较小。

图6 利用动态的目标泛函得到的10个反演结果的平均值Fig.6 Average of 10 inversion results using full waveform inversion with dynamic misfit function

图7 利用动态的目标泛函得到的10个反演结果的归一化后的均方差Fig.7 Mean square deviation of 10 inversion results using full waveform inversion with dynamic misfit function

图8 利用动态的目标泛函得到的10个反演结果对应的归一化后的目标泛函Fig.8 Normalized value of misfit function corresponding to 10 inversion results using dynamic misfit function

3 结 论

(1)在数据约束较强的区域,通过修改目标泛函在已知解位置处的性态能够保证多次反演所得到的反演结果不同。

(2)在数据约束较弱的区域,通过对已有反演解的强制扰动,保证每次所得到的解均不同,从而对该区域的反演结果可信度做出评估。

(3)该方法不需要对后验协方差矩阵求解,也避免了“Bootstraping”方法只能适用于超定问题的局限性,只需要在不同的求解过程中修改目标泛函,因此本文方法具有较强的适用性。

[1] TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation [J]. Geophysics, 1984,49(8):1259-1266.

[2] PLESSIX R E. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications [J]. Geophys J Int, 2006,167(2):495-503.

[3] PRATT R, SHIN C, HICKS G. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion [J]. Geophys J Int, 1998,133(2):341-362.

[4] SHIN C, YOUNG H C. Waveform inversion in the Laplace domain [J]. Geophys J Int, 2008,173(3):922-931.

[5] SHIN C, YOUNG H C. Waveform inversion in the Laplace-Fourier domain [J]. Geophys J Int, 2009,177(3):1067-1079

[6] 杨积忠,刘玉柱,董良国.变密度声波方程多参数全波形反演策略[J].地球物理学报,2014,57(2):628-643.

YANG Jizhong, LIU Yuzhu, DONG Liangguo. A multi-parameter full waveform inversion strategy for acoustic media with variable density [J]. Chinese Journal Geophysics, 2014,57(2):628-643.

[7] KÖHN D. Time domain 2D elastic full waveform tomography [D]. Kiel: Kiel University, 2011.

[8] MORA P. Nonlinear two-dimensional elastic inversion of multi-offset seismic data [J]. Geophysics, 1987,52(9):1211-1228.

[9] VIRIEUX J, OPERTO S. An overview of full-waveform inversion in exploration geophysics [J]. Geophysics, 2009,74(6):WCC1-WCC26.

[10] ALAN R L. Fourth-order finite-difference P-SV seismograms [J]. Geophysics, 1988,53(11):1425-1436.

[11] 李庆洋,李振春,黄建平,等. 基于贴体全交错网格的起伏地表正演模拟影响因素 [J]. 吉林大学学报(地球科学版),2016,46(3):920-929.

LI Qingyang, LI Zhenchun, HUANG Jianping, et al. Factor analysis of seismic modeling with topography based on a fully staggered body-fitted grids [J]. Journal of Jilin University (Earth Science Edition),2016,46(3):920-929.

[12] 李娜,李振春,黄建平,等.Lebedev网格与标准交错网格耦合机制下的复杂各向异性正演模拟[J].石油地球物理勘探,2014,49(1):121-131.

LI Na, LI Zhenchun, HUANG Jianping, et al. Numerical simulation with coupling Lebedev and standard staggered grid schemes for complex anisotropic media [J]. Oil Geophysical Prospecting, 2014,49(1):121-131.

[13] 李振春,李庆洋,黄建平,等.一种稳定的高精度双变网格正演模拟与逆时偏移方法[J]. 石油物探, 2014,53(2):127-136.

LI Zhenchun, LI Qingyang, HUANG Jianping, et al. A stable and high-precision dual-variable grid forward modeling and reverse time migration method [J]. Geophysical Prospecting for Petroleum, 2014,53(2):127-136.

[14] 黄建平,杨宇,李振春,等. 基于 M-PML 边界的 Lebedev 网格起伏地表正演模拟方法及稳定性分析[J].中国石油大学学报(自然科学版), 2016,40(4):47-56.

HUANG Jianping, YANG Yu, LI Zhenchun, et al. Lebedev grid finite difference modeling for irregular free surface and stability analysis based on M-PML boundary condition [J]. Journal of China University of Petroleum (Edition of Natural Science), 2016,40(4):47-56.

[15] BIONDI B, ALMOMIN A. Tomographic full-waveform inversion (TFWI) by combining FWI and wave-equation migration velocity analysis [J]. The Leading Edge, 2013,32(9):1074-1080.

[16] BIONDI B, ALMOMIN A. Simultaneous inversion of full data bandwidth by tomographic full-waveform inversion [J]. Geophysics, 2014,79(3):WA129-WA140.

[17] WANG H, SINGH S C, AUDEBERT F, et al. Inversion of seismic refraction and reflection data for building long-wavelength velocity models [J]. Geophysics, 2015,80(2):R81-R93.

[18] FICHTNER A, KENNETT B L N, IGEL H, et al. Theoretical background for continental-and global-scale full-waveform inversion in the time-frequency domain[J]. Geophys J Int, 2008,175(2):665-685.

[19] MORA P. Inversion=migration+tomography[J]. Geophysics, 1989,54(12):1575-1586.

[20] 董良国,迟本鑫,陶纪霞,等.声波全波形反演目标函数性态[J].地球物理学报,2013,56(10):3445-3460.

DONG Liangguo, CHI Benxin, TAO Jixia, et al. Objective function behavior in acoustic full-waveform inversion [J]. Chinese Journal Geophysics, 2013,56(10):3445-3460.

[21] LOMAX A, SNIEDER R. Finding sets of acceptable solutions with a genetic algorithm with application to surface wave group dispersion in Europe [J]. Geophys Res Lett, 1994,21(24):2617-2620.

[22] SAMBRIDGE M. Exploring multidimensional landscapes without a map [J]. Inverse Problems, 1998,14(3):427-440.

[23] SAMBRIDGE M. Geophysical inversion with a neighborhood algorithm-I. searching parameter space [J]. Geophys J Int, 1999a,138(2):479-494.

[24] SAMBRIDGE M. Geophysical inversion with a neighborhood algorithm-II. appraising the ensemble [J]. Geophys J Int, 1999b,138(3):727-746.

[25] KOPER K D, WYSESSION M E, WIENS D A. Multimodal function optimization with a niching genetic algorithm: a seismological example [J].Bull Seism Soc Am, 1999,89(4):978-988.

[26] MOSEGAARD K, SAMBRIDGE M. Monte Carlo analysis of inverse problems [J]. Inverse Problems, 2002,18(3):R29-R54.

[27] TARANTOLA A. Inverse problem theory [M]. Amsterdam:Elsevier, 1987:44-57.

[28] EFRON B, TIBSHIRANI, ROBERT J. An introduction to the bootstrap [M]. New York: Chapman and Hall/CRC, 1993.

[29] RAWLINSON N, M SAMBRIDGE, SAYGIN E. A dynamic objective function technique for generating multiple solution models in seismic tomography [J]. Geophys J Int, 2008,174(1):295-308.

[31] 李振春,雍鹏,黄建平,等. 基于矢量波场分离弹性波逆时偏移成像[J]. 中国石油大学学报(自然科学版),2016,40(1):42-48.

LI Zhenchun, YONG Peng, HUANG Jianping, et al. Elastic wave reverse time migration based on vector wavefield seperation [J]. Journal of China University of Petroleum (Edition of Natural Science), 2016,40(1):42-48.

[32] FICHTNER A, TRAMPERT J. Hessian kernels of seismic data functionals based upon adjoint techniques [J]. Geophys J Int, 2011,185(2):775-798.

Multi-solutionanalysisoffullwaveforminversionbasedondynamicmisfitfunction

HUANG Jianping1,2， CUI Chao1,2

(1.SchoolofGeosciencesinChinaUniversityofPetroleum,Qingdao266580,China;2.LaboratoryforMarineMineralResources,QingdaoNationalLaboratoryforMarineScienceandTechnology,Qingdao266071,China)

Full wave-form inversion (FWI) is one of the most promising velocity model building tools for its high precision and resolution which theoretically can utilize all the information carried by seismic record. However, because of the multi-solution property of the geophysical inversion problem and the limitation of the quality of seismic data, FWI can only generate one of the local optimize solutions which partly explain the seismic data. The analysis of the inversion result and the multi-solution property is one of the key problems in geophysical inversion. In order to mitigate the computational burden and strict application conditions in conventional result analysis methods, we introduce a dynamic function to analyze the multi-solution property of FWI. The technique designed for generating multiple solutions is based on the modification of objective function using the information of the previous results. The searching for a new model is achieved by adding a feedback term which creates a local maximum at each point in parameter space filled with previously computed models. For the model where the objective function defines a local maximum, the gradient is randomly perturbed to create a family of distinct solutions which helps evaluate the model where poorly illuminated by seismic data. The application on the sea model verifies that the method is efficient in terms of analyzing the inversion results of conventional FWI while keeping the accuracy. The merit of the proposed method is that it only needs to produce a relatively small ensemble of solutions, since each model will substantially differ from all others to the extent permitted by the data and the computational burden compared with the traditional method is significantly reduced.

full waveform inversion; multi-solution; objective function; reliability analysis

2017-02-18

山东省重大创新工程(2017CXGC1608，2017CXGC1602)；中国科学院战略性先导科技专项(XDA14010303)；国际合作重点项目(41720104006)；国家油气重大专项课题(2016ZX05006-004，2016ZX05014001，2016ZX05002)

黄建平(1982-)，男，教授，博士，研究方向为地震波传播与成像技术。 E-mail:jphuang@upc.edu.cn。

1673-5005(2017)06-0064-07

10.3969/j.issn.1673-5005.2017.06.007

P 631.4

A

黄建平,崔超. 一种基于动态目标泛函的全波形反演多解性评估方法[J].中国石油大学学报(自然科学版),2017,41(6):64-70.

HUANG Jianping, CUI Chao. Multi-solution analysis of full waveform inversion based on dynamic misfit function [J].Journal of China University of Petroleum (Edition of Natural Science), 2017,41(6):64-70.

(编辑 修荣荣)