地震反射波反演二阶优化方法及其应用

程玖兵, 王腾飞, 徐文才

海洋地质国家重点实验室，同济大学，上海 200092

0 引言

地球内部存在复杂的多尺度的非均匀性.它决定了弹性参数(如纵、横波速度与密度等)的空间分布与变化特征，影响地震波的传播与散射，甚至导致波场衰减与各向异性效应.在主动源地震探测过程中，人们利用地震记录(包括初至与反射波等震相)去构建地下弹性参数模型，定位和描述介质内部的不连续界面，进而解析地质构造与沉积模式，预测岩性、孔隙流体以及应力场的空间分布，服务于油气及其他矿产资源勘探开发、区域地质或岩石圈构造演化研究和一些工程与环境地质任务.

以地震波的主波长为参照，弹性非均匀性可以分解到不同尺度，包括长波长、中波长和短波长分量.一般认为，前两个分量构成光滑的“背景模型”，控制波场运动学(如走时与相位)，而短波长分量构成震荡的“扰动模型”，主要影响波场动力学(如振幅).因此，在实践中，背景模型和扰动模型的建立通常是分开进行的(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 正问题

地震波传播过程通常用波动方程来描述，在频率域可简要表示为如下矩阵-矢量形式：

L(x,ω)u(x,ω)=f(x,ω)，

(1)

其中u为波场矢量，f代表震源矢量,x代表空间坐标；L称为复“阻抗”矩阵，是由频率ω和地下介质弹性参数m(x)定义的传播算子，描述波场与参数之间的非线性关系.在声学近似意义下，u通常指声压场；在弹性介质中，u一般包含质点的水平、垂直位移(或速度)分量.基于声压波场理论的地震成像与反演最关心纵波速度这个模型参数，有时也考虑介质的密度变化和吸收衰减.欲更贴切地描述地震波传播的物理过程，需采用更复杂的波动理论，引入横波速度、各向异性参数或品质因子等介质参数.

大量观测与实验表明，地下介质的弹性非均匀性具有复杂的多尺度特征.在地震波成像与反演应用中，弹性参数模型一般被分解为光滑的背景模型m0和震荡的扰动模型m1，即

m(x)=m0(x)+m1(x).

(2)

这种尺度分解一般以地震主波长为参照，比如背景模型包含长波长与中波长参数结构，控制波的传播，而扰动模型由刻画短波长特征的参数分量构成，导致波场散射(Tarantola,1984a;Jannane et al.,1989).假设扰动分量远小于背景分量，总波场可分解为入射场u0与一阶散射场u1之和，即u=u0+u1.于是(1)式演化成一阶Born模拟方程：

L(x,ω)u0(x,ω)=f(xs,ω)，

(3a)

L(x,ω)u1(x,ω)=Q(x,ω)u0(x,ω),

(3b)

1.2 反射波形反演

波动方程反射波形反演(RWI)以模拟与观测反射波数据波形残差最小化为目标，除了重构弹性模型的宏观(或长波长)结构,还可适当恢复部分中波长分量(Xu et al.,2012;徐文才,2020).由于反射波u1受背景模型m0与扰动模型m1共同影响，同时估计这两个模型分量面临极强的非线性与复杂的耦合问题.本文将这个反问题分解成嵌套在一起的两个子反演，在迭代过程中交替更新模型的背景与扰动分量，结合频率延拓策略实现宽谱建模目标.

于是，估计背景模型的反问题可定义为

(4)

(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)

1.2.1 一阶优化

(8)

(9)

于是，对方程(3a)和(3b)关于背景模型求一阶偏导数，得到相应的Fréchet微商(推导见附录A)：

(10)

或写成Jm0=L-1F0，其中F0是维数为q×l的虚源矩阵，其列向量满足：

(11)

同理，方程(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)

根据附录B中的推导，扰动模型的泛函梯度可表示成：

(15)

图1 背景模型与扰动模型泛函梯度计算示意图

在交替反演过程中，一旦扰动模型获得更新，就由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，其对角元素与孔径、照明有关，非对角元素反映带限数据条件下空间结点之间的关联性或耦合性.

(18)

引入海森算子可较好地克服因孔径、照明与频带限制对泛函梯度的模糊化作用，使模型更新更均衡、更合理.这将在后文数值实验中得到验证.

一般来讲，海森矩阵可拆解成两项：

(19)

其中第一项只涉及Fréchet微商，定义了所谓的近似海森矩阵：

(20)

第二项涉及波场关于模型参数的二阶偏导数，即

Rmj=

(21)

当数据残差较大或非线性波场效应(如多次散射)较强时，海森算子第二项的作用不可忽视.对于FWI或LSM问题，考虑这一项有助于处理多次散射效应(Pratt et al.,1998; Métivier et al.,2013).不过，第二项可能会破坏矩阵正定性，在求解牛顿方程时要考虑这一因素.在局部线性假设条件下，可求解如下高斯-牛顿方程：

(22)

获得模型更新方向.因为近似海森矩阵具有正定性，高斯-牛顿法总能给出合理的下降方向，并避免海森算子非线性项引入的计算复杂性.

1.2.3 海森矩阵数值分析

虽然本文二阶优化RWI方法不需要显式计算海森矩阵及其逆，但为了观察海森矩阵的特点，本节以常密度声介质单参数(即速度)反演为例，先给出反射波敏感核与海森矩阵的计算式，然后基于小模型演示背景与扰动速度海森矩阵的差异.给定震源s和检波器r的单次观测，背景慢度的反射波敏感核按(A2)可表示成：

Jm0(x,r,s,ω)=ω2[u0(x,s,ω)G1(x,r,ω)

+u1(x,s,ω)G0(x,r,ω)]

=ω2f(ω)[G0(x,s,ω)G1(x,r,ω)

+G1(x,s,ω)G0(x,r,ω)].

(23)

它涉及入射场与一阶散射场的频率域乘积或时间域褶积，揭示了反射波路径上震源端、检波器端的二次散射对背景慢度变化的敏感性.代(23)入(20)式，近似海森矩阵的元素由模型参数结点x与x′对应的Fréchet微商时间域互相关构成，在频率域写成：

×[G0(x′,s,ω)G1(x′,r,ω)]

×[G1(x′,s,ω)G0(x′,r,ω)],

(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推导扰动模型的反射波敏感核和近似海森矩阵的计算式.

图2 背景模型与扰动模型近似海森矩阵计算示意图

在高频近似意义下，基于射线追踪的正演方法假定地震波仅沿着所谓的费马射线传播.相应地，射线层析利用走时数据更新射线路径上的慢度参数，其中走时敏感核对应模型每个离散单元内射线段的长度(Nolet，2008).波动方程走时层析或全波形反演则依据地震波能量沿着“波路径”的传播，揭示波路径上模型参数对波场传播与散射的作用(Woodward,1992;Tromp et al.,2005)，对应的波形敏感核沿着连接震源、检波器的费马射线向四周展开；敏感核的宽度可用来衡量走时层析、偏移成像以及波形反演的分辨率(Woodward,1992).当地震子波频带无限宽时，带限波路径就蜕变成费马射线.

本文理论推导采用平方慢度模型参数化方式，但数值实验展示更直观的速度模型.如图3，在均匀背景含有一个高斯异常体的模型底部设一反射界面.首先合成反射地震数据，然后按前文公式计算背景、扰动模型的泛函梯度和近似海森矩阵.为了方便观察，图中仅显示单次观测的泛函梯度(或波路径)，而近似海森矩阵则考虑了所有观测的贡献.一方面，背景速度模型的波路径是由大角度(接近180°)相遇的正向入射场(或散射场)与伴随散射场(或入射场)零延迟互相关得到的，波路径的宽度要远小于其长度(如图3b).基于反射波形数据反演背景速度模型时，沿射线方向分辨率最低，垂直射线方向最高.另一方面，扰动速度模型的波路径是由小角度(小于90°)相遇的正向与伴随入射场零延迟互相关得到的，对应于叠前偏移的“微笑型椭圆”，侧重更新反射界面或模型短波长分量(图3c).

图3 简单模型RWI泛函梯度与近似海森矩阵

如图2所示，近似海森矩阵是由反射波敏感核的互相关得到的.上面二维模型参数结点为100×100，其海森矩阵的维数为10000×10000.如图3d，背景速度模型的近似海森矩阵对角元素幅值最大，非对角元素幅值也不小，故而很稠密，这与相应的反射波敏感核(或波路径)展布有一定宽度有关.相反，在图3e中，除了左上角之外，扰动速度模型的近似海森矩阵几乎是对角绝对占优的，且沿对角方向幅值逐渐变小，这与几何扩散效应有关.浅层扰动模型的敏感核涉及到大角度相遇的正向与反向入射场的相互作用(类似FWI中直达波和透射波对应的敏感核)，因此海森矩阵左上角较稠密.在射线走时层析中，海森算子对泛函梯度的预条件相当于对模型网格内射线总长度(或照明)进行归一化校正；而在本文反射波形反演中，背景模型的海森矩阵很大程度可消除孔径限制、照明不均以及数据带限对泛函梯度的模糊化，使模型更新方向更合理、更均衡，尤其是明显提升垂向分辨率(见图4).可见，海森矩阵与泛函梯度共同影响模型更新方向和反演分辨能力.

图4 背景速度模型更新方向对比

2 数值实验

这部分基于理论模型合成数据和实际地震资料检验本文波动方程反射波形反演二阶优化方法的有效性.针对反问题的病态性和零空间问题，我们借助构造倾角约束的正则化处理(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)进行多尺度反演建模.

图5 SEG/EAGE推覆体模型合成数据

图6a显示了RWI第一次迭代获得的背景速度更新方向.受观测孔径、照明和数据带限效应影响，负梯度方向在模型空间很不均衡，且受一些假象污染.如图6b所示，高斯-牛顿方向引入了基于海森矩阵的去模糊化处理，提供了比较均衡、合理的更新方向(除了在模型底部两侧照明严重不足的区域).由于每一轮迭代都朝着比较可靠的方向前进，高斯-牛顿法最终重构出了更精确、含更多中波长成分的速度模型(如图6c与6d).如图7，共轭梯度法(CG)经过100次迭代，误差函数仅下降到40%左右，而截断高斯-牛顿法(TGN)外循环20次(内部循环最大迭代次数为5，统计平均将近4次)，误差函数下降到了20%以下.把交替反演得到的背景速度和扰动速度叠加起来，重构出宽谱的层速度模型(图8).连同伪井曲线(图9)可以看出，高斯-牛顿RWI基本实现高分辨率速度建模目标.

图6 背景速度模型反演

图7 反射波形反演残差下降曲线

图8 共轭梯度法(左)与高斯-牛顿法(右)反演结果

图9 伪井速度曲线对比

2.2 东海实际地震资料

接下来基于东海长江坳陷某二维拖缆地震资料检验方法实用性.该数据共有851炮，炮间距为37.5 m，最大炮检距为4 km，检波器间隔为12.5 m，记录长度为4 s.用于实验的地震数据已完成滤除涌浪噪声、压制海面多次波等常规处理(图10).为了克服平面外传播效应的影响，对所有炮记录采用了三维到二维的振幅和相位校正(Wang and Rao,2009).初始层速度模型由常规叠前时间偏移处理获得的均方根速度转换而成.图11显示了基于初始模型的逆时偏移(RTM)剖面及沿测线均匀分布的8个共成像点道集(CIG).浅层CIG同相轴向上弯曲，表明速度偏小，而中深层CIG同相轴明显向下弯曲，表明速度偏大，导致反射波没有准确归位，一些断层绕射波没有收敛.

图10 典型共炮道集

图11 基于初始速度模型的逆时偏移

作为对比，先用商业软件基于叠前深度偏移的反射走时层析技术(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)方法更新速度模型的逆时偏移结果.为了便于区分二者在模型更新上的差异，图中偏移剖面叠合了总的速度更新量.可见，按本文二阶优化思路改善了速度更新，提高了地质一致性，使反射波与绕射波归位、聚焦效果进一步提升，更精细地刻画了一些中深层的沉积序列和左侧基底上覆的断层结构.

图12 商业软件反射走时层析速度建模基础上的逆时偏移

图13 本文反射波形反演速度建模基础上的逆时偏移

图14 逆时偏移局部对比

图15 不同反射波形反演算法基础上的逆时偏移(左为速度更新量与偏移剖面叠合图；右为共成像点道集)

3 讨论

本文数值实验只涉及声压场反射波数据纵波速度建模与偏移成像问题.若要逼近实际地下波场传播的物理过程，有时需要考虑纵-横波模式转换与耦合、各向异性乃至吸收衰减等复杂波现象，因此要在建模与成像过程中引入横波速度、密度、品质因子以及各向异性等参数(如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)

该式可写成：

Jm0=L-1F0,

(A4)

其中F0是维数为q×l的虚源矩阵，其列向量满足：

(A5)

因此，就背景模型某个特定网格结点，Fréchet导数可按下式计算：

假定在检波器r处设一点脉冲，则有背景格林函数G0=L-1δ(x-r)以及对应的扰动格林函数LG1=QG0.于是，与检波器单次观测对应的Fréchet导数可表示成：

(A7)

把(A7)代入(8)式，则可构建背景模型对应的泛函梯度：

由于L-1每一列与每个参数结点处脉冲源的格林函数波场相对应，且根据格林函数的互易性，有(L-1)†=L-1(Virieux and Operto,2009)，于是(A8)式化简得到

(A9)

(A10a)

(A10b)

附录B 扰动模型的敏感核函数、泛函梯度与近似海森矩阵

由(12)与(13)式可知，扰动模型对应的Fréchet导数满足：

(B1)

引入与某检波器r处点脉冲对应的背景格林函数G0，则有：

(B2)

(B2)代入(8)式，则把扰动模型对应的泛函梯度表示成：

(B3)

同样，利用格林函数互易性，(B3)式进一步可写成：

(B4)

(B5)

给定震源s和检波器r的单次观测，扰动速度模型的反射波敏感核按(B2)式可表示成：

Jm1(x,r,s,ω)=ω2f(ω)[G0(x,r,ω)G0(x,s,ω)].

(B6)

它仅涉及入射场的乘积(时间域褶积)，揭示了一次散射对扰动速度变化的敏感性.注意，敏感核是频率或时间的函数，只有与数据残差相结合才体现相应数据成分对模型更新的作用.

根据(20)式，对应的近似海森矩阵的元素可按如下公式计算：

×G0(x,r,ω)]*×[G0(x′,s,ω)G0(x′,r,ω)].

(B7)

可见，扰动模型的近似海森矩阵是相应Fréchet微商的互相关构建出来的(图2c).注意，上式已在前人文献中给出(如Tang,2009;任浩然等,2013).