地球内部存在复杂的多尺度的非均匀性.它决定了弹性参数(如纵、横波速度与密度等)的空间分布与变化特征,影响地震波的传播与散射,甚至导致波场衰减与各向异性效应.在主动源地震探测过程中,人们利用地震记录(包括初至与反射波等震相)去构建地下弹性参数模型,定位和描述介质内部的不连续界面,进而解析地质构造与沉积模式,预测岩性、孔隙流体以及应力场的空间分布,服务于油气及其他矿产资源勘探开发、区域地质或岩石圈构造演化研究和一些工程与环境地质任务.
以地震波的主波长为参照,弹性非均匀性可以分解到不同尺度,包括长波长、中波长和短波长分量.一般认为,前两个分量构成光滑的“背景模型”,控制波场运动学(如走时与相位),而短波长分量构成震荡的“扰动模型”,主要影响波场动力学(如振幅).因此,在实践中,背景模型和扰动模型的建立通常是分开进行的(Claerbout, 1985).建立背景模型尤其是反演纵、横波宏观速度主要依靠走时层析或偏移速度分析.在反射地震勘探领域,基于叠前深度偏移共成像点道集剩余曲率信息的反射走时层析至今仍是速度建模的主要手段(Stork,1992;Woodward et al., 2008).背景模型一旦确定,基于反射或单次散射信号估计弹性参数扰动可视为一个线性反问题.要么基于线性正演算子的伴随算子,由射线或波动理论基础上的地震偏移生成反射界面图像(Claerbout, 1971; Schneider, 1978;McMechan, 1983);要么利用正演算子的近似逆,借助线性逆散射(Miller et al., 1987;Beylkin and Burridge, 1990; 毛伟建等,2021)或真振幅偏移(Bleistein et al., 2001; 张宇, 2006),甚至基于迭代的线性波形反演或最小二乘偏移,即LSM(Clayton and Stolt, 1981; Tarantola, 1984a;Guy and Rene-Edouard, 1999; 任浩然等, 2013; 陈生昌, 2018)定量估计阻抗界面的反射率或弹性参数扰动.
将背景模型与扰动模型反演过程完全分开,不便于考虑两者的耦合关系,也影响反演的收敛性和分辨率.于是,在高性能计算平台的支撑下,无需进行模型尺度分解,以模拟与观测数据所有震相拟合误差最小化为目标的非线性全波形反演,即FWI(Lailly, 1983;Tarantola, 1984b)在近十年重新引起了关注.它具有“走时层析”和“偏移成像”双重功能,可捕捉不同角度散射场信号对弹性参数长、中、短波长分量的敏感性,因此具备高分辨率建模潜力(Mora, 1989;Neves and Singh, 1996).由于人工源地震数据缺少有效的低频成分,FWI的成功案例多限于透射波、回转波或折射波等大角度波场信号充足的井间地震、长炮检距地面地震数据(Pratt, 1999; Brenders and Pratt, 2007).自2005年前后基于双程波方程的逆时偏移(RTM)投入实际生产以来,FWI作为配套的速度建模技术在一些探区尤其是针对海洋长缆观测地震数据取得了很好的实用效果(如Ratcliffe et al., 2014).然而,大量常规地震数据的炮检距与方位角分布范围有限,大角度波场很难穿透到中、深层,FWI极易受非线性和“周波跳跃”问题的困扰(Virieux and Operto, 2009; 董良国等, 2013).即便采用从低频到高频(Bunks et al., 1995)、层剥离(Wang and Rao, 2009)或散射角滤波(Alkhalifah, 2016)等多尺度策略,或者采用更凸的目标泛函(如Luo et al., 2016; Yang and Engquist, 2018)或者扩展模型(如Huang et al., 2018),都只能部分克服FWI在处理反射波数据时面临的挑战.
面向反射地震数据的非线性反演可以不采用经典FWI的理论框架,除了依据反射波数据拟合误差设计目标函数,还可以依据反射波的成像效果(如聚焦性)设计目标函数.数据域的非线性反射波形反演,即RWI(Clément et al., 2001; Xu et al., 2012)和成像域的波动方程偏移速度分析, 即WEMVA(Sava and Biondi, 2004; Shen and Symes, 2008)都属于这一范畴.它们均采用模型尺度分解思想,只不过把重构背景与扰动模型置于一个联合的反问题之中,由嵌套在一起的两个“子反演”交替进行求解.成像域反射走时层析通过叠前深度偏移估计扰动模型或反射率,针对背景模型的目标函数凸性好,对初始模型要求较低,对应的反演分辨率也低一些(Díaz et al., 2013).数据域反射波形拟合目标函数凸性差,仍受“周波跳跃”影响,故通常在初始迭代阶段采用走时拟合标准,当长波长分量得以重构之后再切换到波形拟合标准,进一步弥补一些模型的中波长分量(如Ma and Hale, 2013; Wang et al., 2018).
由于反射数据对大尺度的密度变化不敏感,而与宏观速度之间存在很强的非线性关系,曾有学者尝试用全局优化方法去反演背景速度模型(如Snieder et al., 1989).考虑到计算可行性,人们在处理实际问题时几乎都采用局部优化方法.在迭代的局部寻优过程中,目标泛函关于模型参数的一阶偏导数和二阶偏导数(即泛函梯度与海森算子)共同定义了模型的更新方向.为了降低算法复杂度,最速下降与共轭梯度等一阶优化方法被广泛用于FWI、LSM和RWI等反问题(Virieux and Operto, 2009; 黄建平等, 2014; 姚刚和吴迪, 2017).不过,泛函梯度易受观测孔径与数据频带限制引起的模糊化作用的影响,甚至受多参数串扰脚印与多次散射效应的污染,故导致反演过程收敛较慢,精度也有较大提升空间.二阶优化思想就是引入海森算子对泛函梯度的预条件处理,从而改善收敛性、提高分辨率、压制多参数耦合(Pratt et al., 1998; Métivier et al., 2013; Wang et al., 2016; Pan et al., 2017).根据采用精确还是近似海森算子,出现了拟牛顿法(如L-BFGS)、高斯-牛顿法以及(全)牛顿法等不同的二阶优化算法实现(Pratt et al., 1998; Nocedal and Wright, 2006).近10年来,已有不少文献讨论针对FWI、LSM等反问题的二阶优化方法(如高凤霞等, 2013; 刘璐等, 2013; 任浩然等, 2013; Métivier et al., 2013; Wang et al., 2016; Pan et al., 2017; Rao and Wang, 2017).但是,针对RWI的二阶优化方法研究才刚开始,关于海森算子的特征、二阶优化RWI的实际效果以及所涉及的算法复杂性都值得去调查(Cheng et al., 2020; Wang et al., 2021; 徐文才, 2020).
本文面向观测孔径有限的常规反射地震数据,在二阶优化反演理论框架下,推导背景和扰动模型的反射波敏感核、泛函梯度以及海森算子,提出并验证基于近似海森矩阵的高斯-牛顿反射波形反演方法.论文结构如下:首先采用波动方程的Born近似推导反射波形拟合目标函数关于背景模型和扰动模型的泛函梯度与海森算子,借助简单模型分析海森矩阵的特征及对泛函梯度的预条件作用.然后建立基于高斯-牛顿算法的二阶优化反射波形反演方法.最后利用理论模型合成数据和东海实际地震资料,检验该方法的有效性与实用性.
1 二阶优化反射波形反演方法 1.1 正问题地震波传播过程通常用波动方程来描述,在频率域可简要表示为如下矩阵-矢量形式:
(1) |
其中u为波场矢量,f代表震源矢量, x代表空间坐标;L称为复“阻抗”矩阵,是由频率ω和地下介质弹性参数m(x)定义的传播算子,描述波场与参数之间的非线性关系.在声学近似意义下,u通常指声压场;在弹性介质中,u一般包含质点的水平、垂直位移(或速度)分量.基于声压波场理论的地震成像与反演最关心纵波速度这个模型参数,有时也考虑介质的密度变化和吸收衰减.欲更贴切地描述地震波传播的物理过程,需采用更复杂的波动理论,引入横波速度、各向异性参数或品质因子等介质参数.
大量观测与实验表明,地下介质的弹性非均匀性具有复杂的多尺度特征.在地震波成像与反演应用中,弹性参数模型一般被分解为光滑的背景模型m0和震荡的扰动模型m1,即
(2) |
这种尺度分解一般以地震主波长为参照,比如背景模型包含长波长与中波长参数结构,控制波的传播,而扰动模型由刻画短波长特征的参数分量构成,导致波场散射(Tarantola, 1984a; Jannane et al., 1989).假设扰动分量远小于背景分量,总波场可分解为入射场u0与一阶散射场u1之和,即u=u0+u1.于是(1)式演化成一阶Born模拟方程:
(3a) |
(3b) |
这里L已退化为背景介质中的传播算子(矩阵),入射场u0经扰动算子(矩阵)Q作用构成二次虚源或散射源.若将式(3a)中的震源改为点脉冲,则由Born模拟方程(组)正演得到背景格林函数G0和扰动格林函数G1,即有u0(x, ω)=f(ω)G0(x, ω)且u1(x, ω)=f(ω)G1(x, ω).在常密度声介质中,标量声压场传播算子可写成L(x, ω)=▽2+ω2M0(x),其中▽2为拉普拉斯算子(矩阵),对角阵M0(x)的非零元素由模型空间的背景慢度平方m0(x)构成,扰动算子(矩阵)Q(x, ω)=ω2M1(x)中的对角阵M1(x)由模型空间的扰动慢度平方m1(x)构成.
1.2 反射波形反演波动方程反射波形反演(RWI)以模拟与观测反射波数据波形残差最小化为目标,除了重构弹性模型的宏观(或长波长)结构, 还可适当恢复部分中波长分量(Xu et al., 2012; 徐文才, 2020).由于反射波u1受背景模型m0与扰动模型m1共同影响,同时估计这两个模型分量面临极强的非线性与复杂的耦合问题.本文将这个反问题分解成嵌套在一起的两个子反演,在迭代过程中交替更新模型的背景与扰动分量,结合频率延拓策略实现宽谱建模目标.
于是,估计背景模型的反问题可定义为
(4) |
且χ(m0)为二范数意义下的误差函数:
(5) |
其中△dm0=u1(m0)|m1-dobs代表给定扰动模型条件下因背景模型误差导致的反射数据拟合误差,“†”代表共轭转置.由于背景模型与反射数据之间存在复杂的非线性关系,由(4)和(5)式定义的反问题的非线性很强.
类似地, 估计扰动模型的反问题可定义为
(6) |
其中△dm1=u1(m1)|m0-dobs为基于更新后的背景模型,由扰动模型误差导致的反射数据拟合误差.由式(3b)可知,扰动模型与反射波响应成线性关系,因此(6)式本质上对应一个线性反射波形反演问题, 即所谓的最小二乘逆时偏移问题(Tarantola, 1984a; Guy and Rene-Edouard, 1999).假定有n次观测(相当于有n个炮检对),有p类不同的物理参数(比如考虑纵、横波速度与密度时,p=3),以q代表每个物理参数离散结点数,则模型参数l=q×p,于是m0和m1是维数为l的向量,传播矩阵L维数为q×q;对单个频率成分,dobs和△dm1是维数为n的向量,u0和u1是维数为q的向量.
在每个频段交替求解上述两个反问题,均可在局部优化理论框架下采用如下迭代格式:
(7) |
式中mj(0)为该频段的初始模型,α(i)为更新步长(通常用线搜索法估计),△mj(i)为第i次迭代的更新方向,由目标泛函关于模型参数的一阶和二阶偏导数来确定.
1.2.1 一阶优化最速下降(SD)法和共轭梯度(CG)法是最常用的一阶优化方法(Nocedal and Wright, 2006).前者沿负梯度方向进行模型更新,即△mj(i)=-gmj(i),其中泛函梯度是误差函数关于模型参数的一阶偏导数,满足:
(8) |
其中
(9) |
于是,对方程(3a)和(3b)关于背景模型求一阶偏导数,得到相应的Fréchet微商(推导见附录A):
(10) |
或写成Jm0=L-1F0,其中F0是维数为q×l的虚源矩阵,其列向量满足:
(11) |
可见,按(10)式计算Fréchet微商是由模型空间各点的虚源f0(l)驱动的.由于算子L、波场u0和u1都是背景模型m0的函数,因此(10)式表示的Fréchet微商是依赖于模型m0本身的.
同理,方程(3a)和(3b)关于扰动模型求一阶偏导数,得到相应的Fréchet微商:
(12) |
或写成Jm1=L-1F1,其中虚源矩阵F1的列向量满足:
(13) |
在特定参数化方式下,关于扰动模型的Fréchet微商是不依赖于模型m1的.例如,对于常密度声介质,若用慢度平方来表征模型,则(13)式右端偏微分项仅与频率有关.因此,在Born近似意义下,针对扰动模型的波形反演可视为线性问题.
如果先计算模型空间每个虚源点的Fréchet微商,后按(8)式构建泛函梯度,会涉及太多波场模拟,在实际应用中不可取.人们通常采用伴随状态法,借助正向波场和逆向(或伴随)波场的零延迟互相关来计算泛函梯度(Tarantola, 1984b; Tromp et al., 2005; Plessix, 2006).根据附录A中的推导,背景模型的泛函梯度满足:
(14) |
式中正向入射场u0和散射场u1由式(3a)和(3b)给定的状态方程所描述,伴随入射场
根据附录B中的推导,扰动模型的泛函梯度可表示成:
(15) |
式中伴随入射场
在交替反演过程中,一旦扰动模型获得更新,就由Born模拟计算正向和伴随散射场,进而构建关于背景模型的泛函梯度,并按(9)式更新背景模型;接下来重新模拟正向和伴随入射场,构建扰动模型的泛函梯度并按(9)式更新扰动模型.如此交替进行直至达到终止条件.当初始模型误差较小时,共轭梯度方法可以稳健地逼近真实模型.
1.2.2 二阶优化牛顿法认为误差函数在初始模型m附近满足二次型,故可将其做二阶泰勒近似:
(16) |
其中海森算子H由误差函数在初始模型m处的二阶偏导数构成,即:
(17) |
假定模型空间包含p类参数,每类含有q个离散结点(即l=q×p),则海森算子可表示成维数为l×l的方阵.若将其排布成p×p个矩阵块,则对角块蕴含与观测孔径、波场照明以及数据带限相关的信息,非对角块主要反映非同类参数之间的耦合效应,也携带了带限数据条件下p类参数空间结点之间的关联信息(Pratt et al., 1998; Wang et al., 2016).对于单参数问题,海森矩阵维度降为q×q,其对角元素与孔径、照明有关,非对角元素反映带限数据条件下空间结点之间的关联性或耦合性.
欲在二次型意义下逼近误差函数的极小值,则引入海森矩阵对泛函梯度的预条件处理,即△mj=-Hmj-1▽mjχ.就地震成像与速度建模而言,即便是二维单参数问题,模型维数的数量级至少是106,海森矩阵非常庞大,不适合显式计算和求逆.因此,更新方向一般通过求解牛顿方程获得,即:
(18) |
引入海森算子可较好地克服因孔径、照明与频带限制对泛函梯度的模糊化作用,使模型更新更均衡、更合理.这将在后文数值实验中得到验证.
一般来讲,海森矩阵可拆解成两项:
(19) |
其中第一项只涉及Fréchet微商,定义了所谓的近似海森矩阵:
(20) |
第二项涉及波场关于模型参数的二阶偏导数,即
(21) |
当数据残差较大或非线性波场效应(如多次散射)较强时,海森算子第二项的作用不可忽视.对于FWI或LSM问题,考虑这一项有助于处理多次散射效应(Pratt et al., 1998; Métivier et al., 2013).不过,第二项可能会破坏矩阵正定性,在求解牛顿方程时要考虑这一因素.在局部线性假设条件下,可求解如下高斯-牛顿方程:
(22) |
获得模型更新方向.因为近似海森矩阵具有正定性,高斯-牛顿法总能给出合理的下降方向,并避免海森算子非线性项引入的计算复杂性.
牛顿方程或高斯-牛顿方程常由共轭梯度(CG)法迭代求解.它们涉及的海森-向量积要么采用二阶伴随状态法(Akcelik et al., 2003),要么采用散射积分法(Chen et al., 2007)来构建.这两种方法有内在联系(Tromp et al., 2005), 具体哪种效率更高取决于反问题的几何规模,尤其是震源-检波器数量比值,以及计算资源情况,例如考虑计算与存储开销的平衡等(Chen et al., 2007;Liu et al., 2015).本文所有数值实验中高斯-牛顿方程的求解都采用散射积分法.此时,近似海森矩阵与模型更新矢量的积可写成
虽然本文二阶优化RWI方法不需要显式计算海森矩阵及其逆,但为了观察海森矩阵的特点,本节以常密度声介质单参数(即速度)反演为例,先给出反射波敏感核与海森矩阵的计算式,然后基于小模型演示背景与扰动速度海森矩阵的差异.给定震源s和检波器r的单次观测,背景慢度的反射波敏感核按(A2)可表示成:
(23) |
它涉及入射场与一阶散射场的频率域乘积或时间域褶积,揭示了反射波路径上震源端、检波器端的二次散射对背景慢度变化的敏感性.代(23)入(20)式,近似海森矩阵的元素由模型参数结点x与x′对应的Fréchet微商时间域互相关构成,在频率域写成:
(24) |
其中“*”代表复共轭,右边第一项代表震源端Fréchet微商零延迟互相关,第二项代表检波器端Fréchet微商零延迟互相关.因为Fréchet微商由震源端和检波器端两项组成,所以两个参数结点对应Fréchet微商的互相关共包含四项,其中不同端的Fréchet微商在空间上几乎不重叠,其互相关贡献很小,故(24)式仅保留了相同端Fréchet微商的互相关,见图 2a和2b.如前文所讲,扰动模型的更新是由嵌套在一起的最小二乘偏移来实现的.已有文献讨论这个过程中如何构建和应用海森算子(如Symes, 2008; Tang, 2009; Chen and Sacchi, 2017),因此仅在附录B推导扰动模型的反射波敏感核和近似海森矩阵的计算式.
在高频近似意义下,基于射线追踪的正演方法假定地震波仅沿着所谓的费马射线传播.相应地,射线层析利用走时数据更新射线路径上的慢度参数,其中走时敏感核对应模型每个离散单元内射线段的长度(Nolet,2008).波动方程走时层析或全波形反演则依据地震波能量沿着“波路径”的传播,揭示波路径上模型参数对波场传播与散射的作用(Woodward, 1992; Tromp et al., 2005),对应的波形敏感核沿着连接震源、检波器的费马射线向四周展开;敏感核的宽度可用来衡量走时层析、偏移成像以及波形反演的分辨率(Woodward, 1992).当地震子波频带无限宽时,带限波路径就蜕变成费马射线.
本文理论推导采用平方慢度模型参数化方式,但数值实验展示更直观的速度模型.如图 3,在均匀背景含有一个高斯异常体的模型底部设一反射界面.首先合成反射地震数据,然后按前文公式计算背景、扰动模型的泛函梯度和近似海森矩阵.为了方便观察,图中仅显示单次观测的泛函梯度(或波路径),而近似海森矩阵则考虑了所有观测的贡献.一方面,背景速度模型的波路径是由大角度(接近180°)相遇的正向入射场(或散射场)与伴随散射场(或入射场)零延迟互相关得到的,波路径的宽度要远小于其长度(如图 3b).基于反射波形数据反演背景速度模型时,沿射线方向分辨率最低,垂直射线方向最高.另一方面,扰动速度模型的波路径是由小角度(小于90°)相遇的正向与伴随入射场零延迟互相关得到的,对应于叠前偏移的“微笑型椭圆”,侧重更新反射界面或模型短波长分量(图 3c).
如图 2所示,近似海森矩阵是由反射波敏感核的互相关得到的.上面二维模型参数结点为100×100,其海森矩阵的维数为10000×10000.如图 3d,背景速度模型的近似海森矩阵对角元素幅值最大,非对角元素幅值也不小,故而很稠密,这与相应的反射波敏感核(或波路径)展布有一定宽度有关.相反,在图 3e中,除了左上角之外,扰动速度模型的近似海森矩阵几乎是对角绝对占优的,且沿对角方向幅值逐渐变小,这与几何扩散效应有关.浅层扰动模型的敏感核涉及到大角度相遇的正向与反向入射场的相互作用(类似FWI中直达波和透射波对应的敏感核),因此海森矩阵左上角较稠密.在射线走时层析中,海森算子对泛函梯度的预条件相当于对模型网格内射线总长度(或照明)进行归一化校正;而在本文反射波形反演中,背景模型的海森矩阵很大程度可消除孔径限制、照明不均以及数据带限对泛函梯度的模糊化,使模型更新方向更合理、更均衡,尤其是明显提升垂向分辨率(见图 4).可见,海森矩阵与泛函梯度共同影响模型更新方向和反演分辨能力.
这部分基于理论模型合成数据和实际地震资料检验本文波动方程反射波形反演二阶优化方法的有效性.针对反问题的病态性和零空间问题,我们借助构造倾角约束的正则化处理(Yu et al., 2020)压制泛函梯度中不合理的高波数假象,规范速度模型的更新,提升反演结果的地质一致性.
2.1 SEG/EAGE二维推覆体模型本实验从SEG/EAGE推覆体模型截取其主体部分,包含逆冲断层、古河道等典型构造与沉积特征(图 5a).模型空间上有501×187个采样点,采样间隔为25×25 m,以主频为15 Hz的雷克子波以100 m为间隔合成120炮道集作为观测记录,最大炮检距为4 km,时间采样为2 ms,记录时长为3 s(见图 5b).为了检验RWI方法,将炮记录初至波切除,仅保留反射波数据.图 5c和5d为初始速度模型及相应的逆时偏移剖面.尽管初始速度长波长分量大致合理,但其误差还是导致一些构造成像失真,分辨率较低.接下来从4 Hz到22 Hz分三个频段(4~8 Hz、10~16 Hz和18~22 Hz)进行多尺度反演建模.
图 6a显示了RWI第一次迭代获得的背景速度更新方向.受观测孔径、照明和数据带限效应影响,负梯度方向在模型空间很不均衡,且受一些假象污染.如图 6b所示,高斯-牛顿方向引入了基于海森矩阵的去模糊化处理,提供了比较均衡、合理的更新方向(除了在模型底部两侧照明严重不足的区域).由于每一轮迭代都朝着比较可靠的方向前进,高斯-牛顿法最终重构出了更精确、含更多中波长成分的速度模型(如图 6c与6d).如图 7,共轭梯度法(CG)经过100次迭代,误差函数仅下降到40%左右,而截断高斯-牛顿法(TGN)外循环20次(内部循环最大迭代次数为5,统计平均将近4次),误差函数下降到了20%以下.把交替反演得到的背景速度和扰动速度叠加起来,重构出宽谱的层速度模型(图 8).连同伪井曲线(图 9)可以看出,高斯-牛顿RWI基本实现高分辨率速度建模目标.
接下来基于东海长江坳陷某二维拖缆地震资料检验方法实用性.该数据共有851炮,炮间距为37.5 m,最大炮检距为4 km,检波器间隔为12.5 m,记录长度为4 s.用于实验的地震数据已完成滤除涌浪噪声、压制海面多次波等常规处理(图 10).为了克服平面外传播效应的影响,对所有炮记录采用了三维到二维的振幅和相位校正(Wang and Rao, 2009).初始层速度模型由常规叠前时间偏移处理获得的均方根速度转换而成.图 11显示了基于初始模型的逆时偏移(RTM)剖面及沿测线均匀分布的8个共成像点道集(CIG).浅层CIG同相轴向上弯曲,表明速度偏小,而中深层CIG同相轴明显向下弯曲,表明速度偏大,导致反射波没有准确归位,一些断层绕射波没有收敛.
作为对比,先用商业软件基于叠前深度偏移的反射走时层析技术(Woodward et al., 2008),从初始模型和相应的共成像点道集出发,经数轮迭代构建层速度模型,并由逆时偏移得到新的成像结果(图 12).由于层速度模型得到修正,2 km以上平缓地层界面成像连续性有一定改善,坳陷内部复杂构造、基底界面成像也有明显改观,CIG同相轴剩余曲率得到大幅度消减.鉴于射线走时层析的分辨率瓶颈,宏观速度模型还不够准确,还能观察到CIG中、深层一些同相轴仍有一定的剩余时差.当初始模型偏离真实情况太远时,反射波形反演仍然受到周波跳跃问题的困扰.为了从相同的初始模型出发稳健地进行速度更新,我们对波动方程RWI算法略微做了调整,即在反演前期阶段采用凸性较好的走时拟合标准(Ma and Hale, 2013;Xu et al., 2019),在稳健地估计出速度模型长波长分量之后,再切换到波形拟合标准,进一步提取一些中波长分量并构建出更准确的层速度模型(图 13a).相应的逆时偏移效果由浅到深都得到了明显提升,CIG同相轴基本都被拉平(图 13b和13c).由图 12b和13b及其局部放大(图 14)对比可知,本文RWI方法提高了速度建模精度,使逆时偏移除了使2 km以上的水平层状地层界面更清晰,坳陷内部和穿刺基底的一些断层绕射波收敛也更彻底,从而高分辨率地刻画了一些复杂的小断块,使深部复杂基底的解释变得更容易.进一步地,图 15对比了基于共轭梯度RWI(CG-RWI)方法和高斯-牛顿RWI(GN-RWI)方法更新速度模型的逆时偏移结果.为了便于区分二者在模型更新上的差异,图中偏移剖面叠合了总的速度更新量.可见,按本文二阶优化思路改善了速度更新,提高了地质一致性,使反射波与绕射波归位、聚焦效果进一步提升,更精细地刻画了一些中深层的沉积序列和左侧基底上覆的断层结构.
本文数值实验只涉及声压场反射波数据纵波速度建模与偏移成像问题.若要逼近实际地下波场传播的物理过程,有时需要考虑纵-横波模式转换与耦合、各向异性乃至吸收衰减等复杂波现象,因此要在建模与成像过程中引入横波速度、密度、品质因子以及各向异性等参数(如Tarantola, 1986; 石玉梅等, 2010).近年来,随着声波方程单参数FWI、LSM和RWI技术走向工业应用,基于更复杂波动方程的多参数反演研究正逐年增多(如Choi et al., 2008; Operto et al., 2013; Alkhalifah and Plessix, 2014; Wang and Cheng, 2017; Yang et al., 2016; 王腾飞, 2017;邹鹏和程玖兵, 2020).为了提高收敛效率、压制参数耦合,考虑海森算子的多参数FWI、LSM方法受到了一些学者的关注(如Wang et al., 2016; Pan et al., 2017; Chen and Sacchi, 2017).近来,徐文才(2020)在弹性波RWI理论框架下研究了纵、横波速度和密度三参数情况下的海森矩阵,提出了弹性各向同性介质反射波形反演高斯-牛顿方法.
本文方法原理和主要理论推导适合声介质也适合(黏)弹性、各向异性介质.一旦考虑更实际的波场物理和多参数问题,就要采用更复杂的波动方程正演引擎,泛函梯度、海森-向量积以及模型更新的计算成本会大幅攀升,这是多参数反射波形反演二阶优化方法走向实际应用必须应对的挑战.
4 结论针对常规地面地震数据中、深层速度建模面临的挑战,本文在波动方程反射波形反演理论框架下,推导了背景与扰动模型的反射波敏感核、泛函梯度以及海森算子,建立了基于近似海森矩阵的反射波形反演高斯-牛顿方法.基于简单模型揭示了两个模型分量近似海森矩阵的特征与差异,证实了海森矩阵对泛函梯度的去模糊化作用可优化速度模型更新方向.从SEG/EAGE推覆体模型合成数据实验结果可知,相比于常用的共轭梯度法,利用海森矩阵的高斯-牛顿法提升了反射波形反演的收敛性与宽谱建模能力.东海长江坳陷拖缆地震数据处理实例表明,商业软件射线走时层析技术的分辨率瓶颈制约了波动方程逆时偏移发挥其潜力;而本文二阶优化反射波形反演方法改善了速度长波长与中波长分量的重构,提高了偏移速度模型的可靠性,从而发挥出逆时偏移的成像分辨率优势,精细刻画了坳陷内部复杂断裂系统以及深部基底形态.本文方法在处理三维实际地震资料和多参数问题时,还面临一些计算层面的挑战,今后可结合数据压缩、下采样与随机优化等技术提升它的实用化水平.
附录A 背景模型的敏感核函数和泛函梯度对方程(3a)和(3b)关于背景模型求一阶偏导数,得到
(A1) |
和
(A2) |
联合(A1)和(A2)式,可得到背景模型对应的Fréchet导数(矩阵):
(A3) |
该式可写成:
(A4) |
其中F0是维数为q×l的虚源矩阵,其列向量满足:
(A5) |
因此,就背景模型某个特定网格结点,Fréchet导数可按下式计算:
(A6) |
假定在检波器r处设一点脉冲,则有背景格林函数G0=L-1δ(x-r)以及对应的扰动格林函数LG1=QG0.于是,与检波器单次观测对应的Fréchet导数可表示成:
(A7) |
把(A7)代入(8)式,则可构建背景模型对应的泛函梯度:
(A8) |
由于L-1每一列与每个参数结点处脉冲源的格林函数波场相对应,且根据格林函数的互易性,有(L-1)†=L-1(Virieux and Operto, 2009),于是(A8)式化简得到
(A9) |
式中逆向(或伴随)入射场
(A10a) |
(A10b) |
由(12)与(13)式可知,扰动模型对应的Fréchet导数满足:
(B1) |
引入与某检波器r处点脉冲对应的背景格林函数G0,则有:
(B2) |
(B2) 代入(8)式,则把扰动模型对应的泛函梯度表示成:
(B3) |
同样,利用格林函数互易性,(B3)式进一步可写成:
(B4) |
其中伴随入射场
(B5) |
给定震源s和检波器r的单次观测,扰动速度模型的反射波敏感核按(B2)式可表示成:
(B6) |
它仅涉及入射场的乘积(时间域褶积),揭示了一次散射对扰动速度变化的敏感性.注意,敏感核是频率或时间的函数,只有与数据残差相结合才体现相应数据成分对模型更新的作用.
根据(20)式,对应的近似海森矩阵的元素可按如下公式计算:
(B7) |
可见,扰动模型的近似海森矩阵是相应Fréchet微商的互相关构建出来的(图 2c).注意,上式已在前人文献中给出(如Tang, 2009; 任浩然等, 2013).
致谢 本文实例是依托国家科技重大专项子课题“长江坳陷油气资源潜力评价”研究靶区的地震数据完成的.感谢中海石油有限公司上海分公司提供的支持.
Akcelik V, Bielak J, Biros G, et al. 2003. High resolution forward and inverse earthquake modeling on terascale computers. //Proceedings of the 2003 ACM/IEEE Conference on Supercomputing. New York, NY, USA: ACM.
|
Alkhalifah T, Plessix R É. 2014. A recipe for practical full-waveform inversion in anisotropic media: an analytical parameter resolution study. Geophysics, 79(3): R91-R101. DOI:10.1190/geo2013-0366.1 |
Alkhalifah T. 2016. Full-model wavenumber inversion: An emphasis on the appropriate wavenumber continuation. Geophysics, 81(3): R89-R98. DOI:10.1190/geo2015-0537.1 |
Beylkin G, Burridge R. 1990. Linearized inverse scattering problems in acoustics and elasticity. Wave Motion, 12(1): 15-52. DOI:10.1016/0165-2125(90)90017-X |
Bleistein N, Stockwell J W Jr, Cohen J K. 2001. Mathematics of Multidimensional Seismic Imaging, Migration, and Inversion. New York: Springer-Verlag.
|
Brenders A J, Pratt R G. 2007. Full waveform tomography for lithospheric imaging: results from a blind test in a realistic crustal model. Geophysical Journal International, 168(1): 133-151. DOI:10.1111/j.1365-246X.2006.03156.x |
Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion. Geophysics, 60(5): 1457-1473. DOI:10.1190/1.1443880 |
Chen K, Sacchi M D. 2017. Elastic least-squares reverse time migration via linearized elastic full-waveform inversion with pseudo-Hessian preconditioning. Geophysics, 82(5): S341-S358. DOI:10.1190/geo2016-0613.1 |
Chen P, Jordan T H, Zhao L. 2007. Full three-dimensional tomography: a comparison between the scattering-integral and adjoint-wavefield methods. Geophysical Journal International, 170(1): 175-181. DOI:10.1111/j.1365-246X.2007.03429.x |
Chen S C. 2018. Seismic reflection data waveform imaging based on reflection wave equation. Oil Geophysical Prospecting (in Chinese), 53(3): 487-501. |
Cheng J B, Wang T F, Xu W C. 2020. Hessian based reflection waveform inversion. //82nd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, doi: 10.3997/2214-4609.202011480.
|
Choi Y, Min D J, Shin C. 2008. Frequency-domain elastic full waveform inversion using the new pseudo-Hessian matrix: experience of elastic Marmousi-2 synthetic data. Bulletin of the Seismological Society of America, 98(5): 2402-2415. DOI:10.1785/0120070179 |
Clément F, Chavent G, Gómez S. 2001. Migration based traveltime waveform inversion of 2-D simple structures: A synthetic example. Geophysics, 66(3): 845-860. DOI:10.1190/1.1444974 |
Claerbout J F. 1971. Toward a unified theory of reflector mapping. Geophysics, 36(3): 467-481. DOI:10.1190/1.1440185 |
Claerbout J F. 1985. Imaging the Earth's Interior. New York: Blackwell Science Inc.
|
Clayton R W, Stolt R H. 1981. A Born-WKBJ inversion method for acoustic reflection data. Geophysics, 46(11): 1559-1567. DOI:10.1190/1.1441162 |
Díaz E, Sava P, Yang T N. 2013. Data-domain and image-domain wavefield tomography. The Leading Edge, 32(9): 1064-1072. DOI:10.1190/tle32091064.1 |
Dong L G, Chi B X, Tao J X, et al. 2013. Objective function behavior in acousticfull-waveform inversion. Chinese Journal of Geophysics (in Chinese), 56(10): 3445-3460. DOI:10.6038/cjg20131020 |
Gao F X, Liu C, Feng X, et al. 2013. Comparisons and analyses of several optimization methods in the application of frequency-domain full waveform inversion. Progress in Geophysics (in Chinese), 28(4): 2060-2068. DOI:10.6038/pg20130450 |
Guy C, Rene-Edouard P. 1999. An optimal true-amplitude least-squares prestack depth-migration operator. Geophysics, 64(2): 508-515. DOI:10.1190/1.1444557 |
Huang G H, Nammour R, Symes W W. 2018. Volume source-based extended waveform inversion. Geophysics, 83(5): R369-R387. DOI:10.1190/geo2017-0330.1 |
Huang J P, Sun Y S, Li Z C, et al. 2014. Least-squares split-step migration based on frequency-division encoding. Oil Geophysical Prospecting (in Chinese), 49(4): 702-707. |
Jannane M, Beydoun W, Crase E, et al. 1989. Wavelengths of earth structures that can be resolved from seismic reflection data. Geophysics, 54(7): 906-910. DOI:10.1190/1.1442719 |
Lailly P. 1983. The seismic inverse problem as a sequence of before stack migrations: Conference on Inverse Scattering, Theory and Application, Society for Industrial and Applied Mathematics, Expanded Abstracts, 206-220.
|
Liu L, Liu H, Zhang H, et al. 2013. Full waveform inversion based on modified quasi-Newton equation. Chinese Journal of Geophysics (in Chinese), 56(7): 2447-2451. DOI:10.6038/cjg20130730 |
Liu Y Z, Yang J Z, Chi B X, et al. 2015. An improved scattering-integral approach for frequency-domain full waveform inversion. Geophysical Journal International, 202(3): 1827-1842. DOI:10.1093/gji/ggv254 |
Luo Y, Ma Y, Wu Y, et al. 2016. Full-traveltime inversion. Geophysics, 81(5): R261-R274. DOI:10.1190/geo2015-0353.1 |
Métivier L, Brossier R, Virieux J, et al. 2013. Full waveform inversion and the truncated Newton method. SIAM Journal on Scientific Computing, 35(2): B401-B437. DOI:10.1137/120877854 |
Ma Y, Hale D. 2013. Wave-equation reflection traveltime inversion with dynamic warping and full-waveform inversion. Geophysics, 78(6): R223-R233. DOI:10.1190/geo2013-0004.1 |
Mao W J, Li W Q, Ouyan W. 2021. Review of seismic inverse scattering migration and inversion. Reviews of Geophysics and Planetary Physics (in Chinese), 52(1): 27-44. |
McMechan G A. 1983. Migration by extrapolation of time-dependent boundary values. Geophysical Prospecting, 31(3): 413-420. DOI:10.1111/j.1365-2478.1983.tb01060.x |
Miller D, Oristaglio M, Beylkin G. 1987. A new slant on seismic imaging: Migration and integral geometry. Geophysics, 52(7): 943-964. DOI:10.1190/1.1442364 |
Mora P. 1989. Inversion=migration+tomography. Geophysics, 54(12): 1575-1586. DOI:10.1190/1.1442625 |
Neves F A, Singh S C. 1996. Sensitivity study of seismic reflection/refraction data. Geophysical Journal International, 126(2): 470-476. DOI:10.1111/j.1365-246X.1996.tb05303.x |
Nocedal J, Wright S J. 2006. Numerical Optimization. New York: Springer.
|
Nolet G. 2008. A Breviary of Seismic Tomography. Cambridge: Cambridge University Press.
|
Operto S, Gholami Y, Prieux V, et al. 2013. A guided tour of multiparameter full-waveform inversion with multicomponent data: from theory to practice. The Leading Edge, 32(9): 1040-1054. DOI:10.1190/tle32091040.1 |
Pan W Y, Innanen K A, Liao W Y. 2017. Accelerating hessian-free gauss-newton full-waveform inversion via l-BFGS preconditioned conjugate-gradient algorithm. Geophysics, 82(2): R49-R64. DOI:10.1190/geo2015-0595.1 |
Plessix R E. 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysical Journal International, 167(2): 495-503. DOI:10.1111/j.1365-246X.2006.02978.x |
Pratt R G, Shin C, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x |
Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part 1:Theory and verification in a physical scale model. Geophysics, 64(3): 888-901. DOI:10.1190/1.1444597 |
Rao Y, Wang Y H. 2017. Seismic waveform tomography with shot-encoding using a restarted L-BFGS algorithm. Scientific Reports, 7(1): 8494. DOI:10.1038/s41598-017-09294-y |
Ratcliffe A, Privitera A, Conroy G, et al. 2014. Enhanced imaging with high-resolution full-waveform inversion and reverse time migration: A North Sea OBC case study. The Leading Edge, 33(9): 986-992. DOI:10.1190/tle33090986.1 |
Ren H R, Huang G H, Wang H Z, et al. 2013. A research on the Hessian operator in seismic inversion imaging. Chinese Journal of Geophysics (in Chinese), 56(7): 2429-2436. DOI:10.6038/cjg20130728 |
Sava P, Biondi B. 2004. Wave-equation migration velocity analysis.I. Theory. Geophysical Prospecting, 52(6): 593-606. DOI:10.1111/j.1365-2478.2004.00447.x |
Schneider W A. 1978. Integral formulation for migration in two and three dimensions. Geophysics, 43(1): 49-76. DOI:10.1190/1.1440828 |
Shen P, Symes W W. 2008. Automatic velocity analysis via shot profile migration. Geophysics, 73(5): VE49-VE59. DOI:10.1190/1.2972021 |
Shi Y M, Yao P C, Sun H S, et al. 2010. Density inversion and porosity estimation using seismic data. Chinese Journal of Geophysics (in Chinese), 53(1): 197-204. DOI:10.3969/j.issn.0001-5733.2010.01.022 |
Snieder R, Xie M Y, Pica A, et al. 1989. Retrieving both the impedance contrast and background velocity: A global strategy for the seismic reflection problem. Geophysics, 54(8): 991-1000. DOI:10.1190/1.1442742 |
Stork C. 1992. Reflection tomography in the postmigrated domain. Geophysics, 57(5): 680-692. DOI:10.1190/1.1443282 |
Symes W W. 2008. Approximate linearized inversion by optimal scaling of prestack depth migration. Geophysics, 73(2): R23-R35. DOI:10.1190/1.2836323 |
Tang Y X. 2009. Target-oriented wave-equation least-squares migration/inversion with phase-encoded Hessian. Geophysics, 74(6): WCA95-WCA107. DOI:10.1190/1.3204768 |
Tarantola A. 1984a. Linearized inversion of seismic reflection data. Geophysical Prospecting, 32(6): 998-1015. DOI:10.1111/j.1365-2478.1984.tb00751.x |
Tarantola A. 1984b. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
Tarantola A. 1986. A strategy for nonlinear elastic inversion of seismic reflection data. Geophysics, 51(5): 1893-1903. |
Tromp J, Tape C, Liu Q Y. 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels. Geophysical Journal International, 160(1): 195-216. |
Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367 |
Wang T F, Cheng J B. 2017. Elastic full-waveform inversion based on mode decomposition: the approach and mechanism. Geophysical Journal International, 209(2): 606-622. DOI:10.1093/gji/ggx038 |
Wang T F, Cheng J B, Guo Q, et al. 2018. Elastic wave-equation-based reflection kernel analysis and traveltime inversion using wave mode decomposition. Geophysical Journal International, 215(1): 450-470. DOI:10.1093/gji/ggy291 |
Wang T F, Cheng J B, Geng J H. 2021. Reflection full waveform inversion with second-order optimization using the adjoint-state method. Journal of Geophysical Research: Solid Earth, 126(8): e2021JB022135. |
Wang Y, Rao Y. 2009. Reflection seismic waveform tomography. Journal of Geophysical Research: Solid Earth, 114(B3): B03304. DOI:10.1029/2008JB005916 |
Wang Y W, Dong L G, Liu Y Z, et al. 2016. 2D frequency-domain elastic full-waveform inversion using the block-diagonal pseudo-Hessian approximation. Geophysics, 81(5): R247-R259. DOI:10.1190/geo2015-0678.1 |
Woodward M J. 1992. Wave-equation tomography. Geophysics, 57(1): 15-26. DOI:10.1190/1.1443179 |
Woodward M J, Nichols D, Zdraveva O, et al. 2008. A decade of tomography. Geophysics, 73(5): VE5-VE11. DOI:10.1190/1.2969907 |
Xu W C, Wang T F, Cheng J B. 2019. Elastic model low-to-intermediate wavenumber inversion using reflection traveltime and waveform of multicomponent seismic data. Geophysics, 84(1): R109-123. DOI:10.1190/geo2018-0306.1 |
Xu S, Wang D L, Chen F, et al. 2012. Inversion on reflected seismic wave. //82nd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, doi: 10.1190/segam2012-1473.1.
|
Yang J Z, Liu Y Z, Dong L G. 2016. Simultaneous estimation of velocity and density in acoustic multiparameter full-waveform inversion using an improved scattering-integral approach. Geophysics, 81(6): R399-R415. DOI:10.1190/geo2015-0707.1 |
Yang Y N, Engquist B. 2018. Analysis of optimal transport and related misfit functions in full-waveform inversion. Geophysics, 83(1): A7-A12. DOI:10.1190/geo2017-0264.1 |
Yao G, Wu D. 2017. Reflection full waveform inversion. Science China Earth Sciences, 60(10): 1783-1794. DOI:10.1007/s11430-016-9091-9 |
Yu Y, Cheng J B, Yang T, et al. 2020. Angle-domain differential reflection traveltime tomography. IEEE Geoscience and Remote Sensing Letters, 17(2): 202-206. DOI:10.1109/LGRS.2019.2920385 |
Zhang Y. 2006. The theory of true-amplitude one-way wave equation migration. Chinese Journal of Geophysics (in Chinese), 49(5): 1410-1430. |
Zou P, Cheng J B. 2020. Visco-acoustic wave equation reflection inversion for the Q model. Chinese Journal of Geophysics (in Chinese), 63(1): 287-297. DOI:10.6038/cjg2020M0574 |
陈生昌. 2018. 基于反射波动方程的地震反射数据波形成像. 石油地球物理勘探, 53(3): 487-501. |
董良国, 迟本鑫, 陶纪霞, 等. 2013. 声波全波形反演目标函数性态. 地球物理学报, 56(10): 3445-3460. DOI:10.6038/cjg20131020 |
高凤霞, 刘财, 冯晅, 等. 2013. 几种优化方法在频率域全波形反演中的应用效果及对比分析研究. 地球物理学进展, 28(4): 2060-2068. DOI:10.6038/pg20130450 |
黄建平, 孙郧松, 李振春, 等. 2014. 一种基于分频编码的最小二乘裂步偏移方法. 石油地球物理勘探, 49(4): 702-707. |
刘璐, 刘洪, 张衡, 等. 2013. 基于修正拟牛顿公式的全波形反演. 地球物理学报, 56(7): 2447-2451. DOI:10.6038/cjg20130730 |
毛伟建, 李武群, 欧阳威. 2021. 地震逆散射偏移与反演综述. 地球与行星物理论评, 52(1): 27-44. |
任浩然, 黄光辉, 王华忠, 等. 2013. 地震反演成像中的Hessian算子研究. 地球物理学报, 56(7): 2429-2436. DOI:10.6038/cjg20130728 |
石玉梅, 姚逢昌, 孙虎生, 等. 2010. 地震密度反演及地层孔隙度估计. 地球物理学报, 53(1): 197-204. DOI:10.3969/j.issn.0001-5733.2010.01.022 |
王腾飞. 2020. 弹性波模式解耦全波形反演方法[博士论文]. 上海, 同济大学.
|
徐文才. 2020. 基于Hessian与P/S解耦的反射地震波波形反演方法[博士论文]. 上海: 同济大学.
|
姚刚, 吴迪. 2017. 反射波全波形反演. 中国科学: 地球科学, 47(10): 1220-1232. |
张宇. 2006. 振幅保真的单程波方程偏移理论. 地球物理学报, 49(5): 1410-1430. DOI:10.6038/cjg2020M0574 |
邹鹏, 程玖兵. 2020. 黏声方程Q值反射波反演. 地球物理学报, 63(1): 287-297. |