考虑水层鸣震的海底多分量数据转换波走时反演方法

刘学义, 程玖兵, 王腾飞

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

0 引言

在现代地震观测和高性能计算技术的推动下，多分量地震方法在油气资源勘探领域越来越受到重视.它更全面地记录了地震纵波(P)与横波(S)的矢量振动信息，用于改善“气云区”和复杂地质体(如盐丘、火成岩等)地震成像，也有助于区分岩性，刻画含油气储层及其内部非均质性，检测定向发育的裂缝系统以及分析应力分布特征等(Stewart et al.,2002; 赵邦六,2008).

在多分量数据处理流程中，PS转换波偏移成像处于核心地位，其成败很大程度取决于地震记录P/S分离的可靠性，以及P波和S波偏移速度模型的合理性(芦俊等,2018).通常，在分离出PP波和PS波地震记录之后，采用传统的偏移速度分析或走时层析技术就能得到P波速度模型，而估计S波速度则要遵循PS波入射-反射路径非对称性，采用相适应的速度分析或走时反演方法(Broto et al.,2003;杜启振等,2011;Wang et al.,2019).随着基于波动方程的逆时偏移(RTM)和全波形反演(FWI)等高精度成像与建模技术的推广应用，面向多分量地震数据的RTM与FWI方法也逐渐受到关注(Sears et al.,2008;Ren and Liu,2016).当地震记录缺少可靠的低频成分和大偏移距信号时，FWI很难合理地重构中深层的纵、横波宏观速度模型，影响反射PP波与PS波成像.

这些年来，强化FWI的反射层析功能或发展针对反射波信号的反演方法成为了勘探地球物理界的研究热点，其中有代表性且与本文相关的工作包括Xu 等(2012)提出的反射波形反演(RWI)以及由它退化而成的波动方程反射走时反演(RTI)方法 (Ma and Hale,2013).虽然理论上这些方法经过拓展也能处理多分量地震数据，但都会面临P-S波串扰、多参数(如纵-横波速度)耦合造成的强烈非线性以及由此引起的收敛性问题.为此，除了基于海森矩阵的二阶优化方法(Epanomeritakis et al.,2008)或子空间方法(Baumstein,2014)，近年一些学者提出了P/S模式解耦预条件的弹性波FWI(Wang and Cheng,2017)、弹性波RTI/RWI(Wang et al.,2018)以及弹性波动方程偏移速度分析方法(Wang et al.,2019).在此基础上，Xu 等(2019)基于P/S模式解耦的数据或梯度预条件，建立了串联标量声波方程反射走时反演、弹性波方程RWI和FWI的多级速度建模方法，有效地利用了多分量地震数据的走时与波形信息，构建出高精度的纵、横波速度模型.值得注意的是，上述方法处理实际海底多分量地震资料时要克服一系列挑战.例如，海面的自由界面效应和海底的上-下行波耦合会极大地增加海底多分量地震信号的复杂性.

海底上-下行P/S波分离对多分量地震数据偏移速度分析是非常必要的.当前，海底多分量记录的波场分离方法分为两大类：一类是基于声波方程的分离方法，如PZ叠加(Dragoset and Barr,1994;Osen et al.,1999)，它可以很好地压制鬼波对水听器P波记录的干扰；另一类是基于弹性波方程的分离方法(Schalkwijk et al.,2003;Muijs et al.,2007)，它同时完成上-下行波分离和P/S解耦，压制检波器端的自由表面多次波，但分离的上、下行P波和上行S波记录仍包含源端自由表面多次波.比如，针对由绳系海底地震仪(OBS)在我国东海QY探区采集的多分量地震资料，刘学义等(2021)根据噪声特点对传统海底数据P/S模式分离方法与流程进行了调整，发现分离后的P波和S波记录仍受源端多次波干扰，给后续处理尤其是速度分析和偏移成像带来极大的困扰.

所以，在完成上-下行P/S分离之后，海底多分量地震数据处理还需克服源端自由表面多次波的影响(Brown and Yan,1999).对于P波记录，通行的做法是利用预测滤波(Wiggins,1988)或基于波动方程的自由表面多次波预测-相减法(SRME: Verschuur et al.,1992)进一步滤除源端多次波，然后进行速度分析与偏移成像.近年来，人们也提出了一些视多次波为有效信号的偏移成像方法(刘伊克等,2018;Liu et al.,2020).不过，所有这些方法都仅关注反射PP波.在PS波速度分析与偏移成像中考虑源端P波水层鸣震，需克服海底介质的复杂影响和PS波传播路径的非对称性，挑战很大，但相关研究比较少(Dawson et al.,2010).

本文基于海洋声学与地球物理文献资料，分析地震波在典型海底界面的能量分配规律，调查软海底情况下源端多次波的类型，明确源端水层鸣震对P/S分离之后获得的上行S波数据(即本文所说的转换波数据)的主导性干扰(见图1至图3).若已知纵波速度模型和PP波成像结果，针对横波速度建模问题，我们提出一种考虑源端水层鸣震的波动方程转换波反射走时反演方法,以摆脱常规方法因鸣震干扰带来的数据拟合困难和收敛性问题.

1 上行S波数据多次波干扰分析

海底多分量地震记录受到海面和海底这两个特殊界面的双重影响.大量海洋声学调查和地质、地球物理资料表明，除了基岩或老地层直接出露的局部区域，海底大都覆盖一层未固结软性沉积物，具有速度随深度逐渐增加、泊松比或纵-横波速度比较大的特性(Carlson et al.,1984; Berge et al.,1991).海底介质的纵波速度一般低于两千米每秒,横波速度在几十到几百米每秒之间(Hamilton,1979;Sidler and Holliger,2010).海底表层沉积物孔隙度较大，随着深度增加介质孔隙度逐渐变小，密度逐步增大，以孔隙度为35%的砂质沉积为例，密度在2100 kg·m-3左右(Buckingham,2005).针对较普遍的软质海底环境，下面分析地震波在海底界面的能量分配规律，推断在海底观测的上行S波中的源端多次波的主导成分.在给定海水P波速度为1.5 km·s-1、密度为1.02 g·cm-3的情况下，依据海底介质纵、横波速度和密度取值范围，各四次均匀采样构成64种参数组合，然后用Zoeppritz方程的线性近似(Aki and Richards,2002)计算海底界面处P波、转换S波的透射与反射系数，最后统计并绘制其临界角之前的平均值(如图1所示，四个大图板代表四种横波速度，每个图板中的四列代表四种纵波速度，每列四组代表四种密度取值).图2展示了其中两种典型软海底界面的透射/反射系数.不难发现，尽管大部分P波能量经海底透射到地层中，但其中的透射转换S波并不强；即便对于横波速度较大、泊松比较小的海底地层，海底界面的PP波反射也强于PS透射波.所以在软海底情况下，在海底界面发生的模式转换对最终分离出的S波记录贡献可以忽略(如图2).这种能量分配特点表明，软海底多分量地震记录的源端多次波干扰主要是一些P波水层鸣震或自由表面多次波.这与前人一些研究相吻合(Caldwell,1999; Suarez,2000).因此，在海底界面实施上-下行P/S分离之后，S波数据中除了源自于海底之下一些阻抗界面的反射PS波，也包含一些源端多次波，其中由海面与海底各反射一次的PPS波(图3粗实线)路径最短、能量最强.而当海底以下不存在盐丘等强波阻抗界面时，地下界面反射波场相关的PPS波能量比较弱(图3细虚线),本文不予考虑.后文数值试验部分基于一阶速度-应力方程模拟的海底多分量地震记录也与上述分析相吻合.

图1 海底界面反射/透射系数平均值统计图板 纵轴代表透射/反射系数平均值；横轴为纵/横波速度四次采样；每一组含4次密度采样， 取值依次为1.5 g·cm-3、1.7 g·cm-3、1.9 g·cm-3和2.1 g·cm-3.Fig.1 Average statistics of transmission/reflection coefficient at seabed interface The vertical axis represents the average of transmission/reflection coefficient; the horizontal axis represents four samples of P-/S- wave velocities, and each group has 4 samples (1.5 g·cm-3、1.7 g·cm-3、1.9 g·cm-3 and 2.1 g·cm-3) of density.

图2 两种典型软海底界面的反射/透射系数曲线 (a) vP=1.7 km·s-1,vS=0.1 km·s-1,ρ=1.9 g·cm-3; (b) vP=1.7 km·s-1,vS=0.5 km·s-1,ρ=1.9 g·cm-3.Fig.2 Transmission/reflection coefficient curves of two typical soft seabed interfaces

图3 海底上行S波传播路径示意图 黑线表示P波传播路径，红线表示S波传播路径.Fig.3 Schematic diagram of propagation paths of up-going S-waves The black lines are the paths of P-waves, and the red lines are the paths of S-waves.

2 考虑源端水层鸣震的转换波走时反演方法

国内外迄今采集了很多海底多分量地震数据.为了利用纵波数据，最常用的前处理手段是联合水检和陆检信号压制检波器端海水鸣震(李维新等,2018).为了充分利用记录到的纵波与横波信号，通常要在海底界面上实施上-下行P/S 波分离，压制检波器端的水层鸣震或海面多次波，残余的源端多次波还需要其他多次波消去法进行压制(Edme and Sing,2009).近年来，陆续出现了一些可同时处理一次反射波、表面多次波(甚至层间多次波)的偏移成像和速度分析方法(Lu et al.,2015;李志娜和李振春,2016;叶月明等,2019;Liu et al.,2020),但这些研究仍局限于纵波数据.

根据前文分析，在上-下行P/S 波分离得到的上行S波数据中，除了人们需要的PS波，还包括PPS波等源端水层鸣震和其他噪声.我们在前期研究中发现，这些源端多次波会给基于转换波数据(上行S波数据)的横波速度分析带来极大的困扰.为此，下文在波动方程反射走时反演理论框架下，首先给出反射PS波及含源端水层鸣震路径的S波(即PPS波)的Born模拟方程，然后建立转换反射波走时拟合目标泛函，分析PS波和PPS波这两种成分对泛函梯度或反射波敏感核的影响，进而提出合适的梯度预条件方法.这里假设前期处理已获得可靠的纵波速度模型、海底界面位置以及PP波成像结果，聚焦宏观横波速度模型构建问题.

2.1 正问题

考虑到走时反演对弹性效应不敏感，故采用声波方程来相对经济地描述纵波和转换横波的传播过程.常密度标量声波方程在频率域可写成：

L(mp(x),ω)u(x,xs,ω)=-F(ω)δ(x-xs),

(1)

L(mp0(x),ω)up0(x,xs,ω)=-F(ω)δ(x-xs),

(2)

L(mp0(x),ω)up1(x,xs,ω)=Kp(x,ω)up0(x,xs,ω),

(3)

式中Kp(x,ω)=ω2mp1(x)是高波数纵波慢度扰动形成的二次虚源，一阶散射场up1通常称为PP波.对于一次转换/反射PS波，替换(3)式中的模型参数，也可用一阶Born模拟合成运动学可靠的数据，表示为

L(ms0(x),ω)us0(x,xs,ω)=Ks(x,ω)up0(x,xs,ω),

(4)

式中Ks(x,ω)=ω2ms1(x)，ms0和ms1分别代表横波慢度平方模型的背景和扰动分量，us0表示PS波，其传播路径如图3所示.

对于源端水层鸣震中占主导的PPS波，可认为是经海底反射到达海面的一次反射波作为新的震源，经水层传播到地层中在一些阻抗界面发生P-S模式转换与反射，最后上传到海底被检波器接收.假定已通过海洋声学调查或地球物理方法获得海底界面信息(Yang et al.,2015;Perdomo et al.,2018)，则可由公式(2)和(3)模拟仅仅由海底界面反射到达海面的PP波；然后借助下面两个方程模拟出运动学可靠的PPS波，即

L(mp0(x),ω)up2(x,xs,ω)=-up1(x,xs,ω)δ(x-x0),

(5)

L(ms0(x),ω)us1(x,xs,ω)=Ks(x,ω)up2(x,xs,ω),

(6)

(7)

其中加权系数w0和w1用于调节两种波场成分的相对幅值，其大小主要与两种转换波传播路程引起的几何扩散效应差异(在浅水区可忽略不计)以及PPS波发生在海底、海面的两次反射有关.加权系数的设置方法参见附录A.

2.2 反射波走时拟合目标泛函及其梯度计算

波动方程反射走时反演的关键在于定量估计模拟反射波与观测反射波的走时拟合误差.主流方法包括动态图像匹配(DIW)(Hale,2013)和局部互相关(Van Leeuwen & Mulder,2010).Xu 等(2019)将Ma 和Hale(2013)提出的基于DIW的反射走时反演方法拓展到了多分量地震数据，但他们没有考虑水层鸣震等多次波干扰.在将分离出来的上行S波记录作为观测数据与模拟数据进行走时匹配时，必然受到一些较强的源端水层鸣震(如PPS波)的影响(可见数值试验部分的例证).为此，我们认为既要在模拟转换反射波时同时考虑PS波和PPS波，又要在反演过程中避免二者相互串扰对模型更新方向的畸变.下面介绍本文的方法和处理步骤.

τ=argminlD(l(xr,xs,t)),

(8)

式中l(xr,xs,t)为时空变化的时移量，D是l的函数，用来描述观测数据与模拟数据的L2范数残差，即：

(9)

于是将转换波走时反演的目标函数定义为

(10)

在局部优化算法基础上，逐步更新横波速度模型使目标函数最小化.

借鉴Ma 和Hale(2013)的PP反射波走时反演方法，可推导出转换波(含PS波和PPS波)反射走时反演的泛函梯度表达式：

(11)

(12)

注意，公式(11)表示泛函梯度是由正向扰动场与伴随入射场零延迟互相关获得，代表地下反射界面与海底检波器之间上行波路径对横波速度的敏感性，这与PS、PPS转换波运动学特征一致.

假设背景横波速度扰动较小，真实反射数据可由时移后的模拟反射数据替代(Luo and Schuster,1991)，即：

(14)

则伴随源可简化为DIW估计的走时残差τ(xr,xs,t)对模拟记录偏导数波场的加权，进而把公式(12)改写为

2.3 转换反射波走时敏感核分析

这部分基于水平层状介质模型(图4)合成OBC数据进行泛函梯度或敏感核分析.震源位于(360 m,15 m)处，采用主频为20 Hz的雷克子波，在海底(深度200 m)以5 m的间隔接收模拟的多分量地震信号(图4a).图5显示按本文Born模拟方法合成的含PS波和PPS波的转换波记录.

图4 水平层状介质模型 (a) 纵波速度； (b) 横波速度.五角星表示炮点位置，倒三角表示检波器位置.Fig.4 Horizontal layered model (a) P-wave velocity； (b) S-wave velocity. The star is for source, and the triangles are for receivers.

图5 共炮点转换波数据 右侧子图中黑线表示P波路径，红线表示S波路径.Fig.5 Converted waves of common-shot gather In the rightsubgraphs, black lines represent P-waves path, and red lines represent S-waves path.

图6显示按前文公式计算的单炮检对(炮点、检波器位置已在图4b中标识)对应的转换波走时关于横波速度的敏感核，包含低波数的转换反射波路径和高波数的偏移脉冲响应.图6a最为复杂，这是由PS波、PPS波的偏移脉冲响应以及转换反射波路径相互混叠造成.对于模拟数据，可在梯度计算过程中对参与互相关的伴随入射波场和正向扰动波场进行PS波与PPS波分离，进而获得四种分解后的敏感核(图6b-e).可见，通过对转换波记录进行波模式分离，可以摒弃串模式互相关产生的高波数干扰(图6d和6e)，有助于压制对低波数转换反射波路径的污染(见图6b和6c).注意，图6b和6c中的低波数敏感核形状虽然相似，但倾斜度有差异(后者稍陡一些)，这与该炮检对PS和PPS波在检波器端反射路径的倾角差异相吻合.图中箭头指示的残余高波数干扰是单一模式转换波记录的偏移脉冲响应，对更新宏观速度模型没有帮助，经多炮检对敏感核叠加后明显减弱，也可在波场互相关时滤除小角度相遇的正向与伴随波场进一步加以压制.

2.4 梯度预条件处理

前文敏感核数值分析表明，对参与泛函梯度计算的正向和伴随波场进行PS波与PPS波模式分离，可以明显改善泛函梯度的合理性.基于声波方程Born模拟算子可以单独模拟正向传播的这两种震相，但由于直接在观测记录中区分PS与PPS很困难，因此获得相应的伴随波场并不容易.

由于公式(14)引入了对观测记录的运动学替代处理，因此只要在Born模拟时区分这两种震相，就可以按下面方程分别实现PS和PPS伴随波场的逆时传播，即

(16)

(17)

其中αi取值范围为0～1，且满足α0+α1=1.当α0=1且α1=0时，表示只用PS波走时信息；当α0=0且α1=1时，表示只用PPS波走时信息.考虑到波场中PS波能量占主导，α0取值通常不小于α1.上式表示将类似图6b和图6c中的单震相转换反射波敏感核加权叠加，获得合理的横波速度更新方向.

图6 单炮检对情况下横波速度对应的泛函梯度及其分解 (a) 原始泛函梯度； (b) PS波泛函梯度； (c) PPS波泛函梯度； (d) 正向扰动PPS波与伴随入射PS波互相关产生的梯度假象； (e) 正向扰动PS波与伴随入射PPS波互相关产生的梯度假象.Fig.6 The gradient with respect to S-wave velocity and its decomposition for a single pair of shot and receiver (a) Original functional gradient； (b) Functional gradient of PS waves； (c) Functional gradient of PPS waves； (d) Gradient artifacts generated by the cross-correlation between forward perturbed PPS waves and adjoint incident PS waves； (e) Gradient artifacts generated by the cross-correlation between forward perturbed PS waves and adjoint incident PPS waves.

一旦确定更新方向，则按局部寻优方法进行背景横波速度模型的更新：

(18)

3 数值实验

转换波走时反演流程图如图7所示.本文聚焦于利用转换波走时信息进行宏观横波速度建模，所以假设前期处理已获得可靠的纵波速度模型、海底界面位置以及PP波成像结果，并且已经完成海底多分量数据的P/S分离工作.

图7 转换波走时反演流程图Fig.7 Flowchart of converted wave travel time inversion

3.1 水平层状模型

首先基于水平层状介质模型检验方法有效性.该模型水深为200 m，海底界面以下纵、横波速度随深度线性增加，在海底附近泊松比较大(图8).海水表面采用自由表面边界条件，基于弹性波速度-应力方程有限差分算法合成多分量观测数据，震源为主频25 Hz的雷克子波，炮间距为50 m，共40炮均匀分布于水深5 m处，检波器置于海底(如图8a中三角形所示)，间隔为5 m，最大炮检距为700 m，记录时长1.15 s.图9a至9c显示了炮点位于(1250 m，5 m)处的三分量记录.基于弹性波理论进行上-下行P/S波分离(刘学义等,2021)，获得的海底上行S波数据如图9d所示.可见，它仍然受到源端自由表面多次波的强烈干扰，也存在一些较弱的层间多次波噪声.总体看来，经海底、海面各反射一次的水层鸣震形成了在源端多次波中占主导的PPS波信号.

图8 水平层状介质模型 (a) 纵波速度； (b) 横波速度； (c) 密度.Fig.8 Horizontal layered medium model (a) P-wave velocity； (b) S-wave velocity； (c) Density.

图9 数值模拟的三分量记录和分离后的转换波数据 (a) 声压分量； (b) vz分量； (c) vx分量； (d) 海底分离出的上行S波.Fig.9 Synthetic 3C common shot gathers and the decomposed up-going S-waves (a) Pressure； (b) Vertical particle velocity； (c) Horizontal particle velocity； (d) Decomposed up-going S-waves.

为了反演横波速度，假定已知准确的纵波背景速度、海底界面模型和反射PP波成像结果.在这里，纵波背景速度由真实速度垂向平滑获得(图10a).考虑走时反演对波场模拟振幅精度要求不高，而且可以由权系数调节合成转换波记录中的PS波和PPS波相对振幅关系，因此要求海底界面模型深度准确、能大致反映从水层到海底的速度扰动即可，文中不再展示.PP成像结果由反射PP波记录逆时偏移获得(图10b).初始横波速度模型是真实层速度垂向大尺度平滑后并减去一个常量得到(图11a)，密度设为常数.公式(17)中两种震相权重相等，即α0=α1=0.5.接下来进行两组反演试验：第一组的正演引擎只能模拟PS波，则出现模拟记录与含有较强PPS波信号的观测记录很难匹配的问题，获取不到比较合理的走时残差，导致反演结果不合理(图11b和图12蓝色实线).第二组的正演引擎同时模拟出PS波和PPS波信号，转换波走时残差估计变得相对容易，最终得到了合理的反演结果(图11c和图12红色实线).注意，模型底部速度无法更新是因为没有反射波场照明这一层.

图10 先验模型 (a) 背景纵波速度； (b) 反射PP波逆时偏移结果.Fig.10 The prior model (a) Background P-wave velocity； (b) Reverse time migration (RTM) of the reflected PP waves.

图11 含源端水层鸣震的转换波数据横波速度反演结果 (a) 初始背景横波速度； (b) 正演引擎无法模拟PPS波时的反演结果； (c) 本文方法走时反演结果.Fig.11 The inversion S-wave velocity with converted wave recordings including water-layer source-side reverberations (a) The initial background S-wave velocity； (b) The result of travel-time inversion with a forward modeling engine that can not simulates PS waves； (c) The proposed method.

图12 伪井横波速度曲线Fig.12 Pseudo-well curves of the S-wave velocity models

图13把两组实验由初始和反演的横波速度模型模拟的转换波记录与原始的观测记录进行对比.当正演引擎不能模拟PPS波时，会导致走时匹配紊乱，第二个PS波同相轴错误地匹配了观测记录的PPS波同相轴.一旦正演引擎可以同时模拟PS波和PPS波，即便初始模型对应的模拟数据与观测数据存在较大运动学误差，通过多轮反演速度模型也能朝着正确方向进行更新，最终使模拟和观测记录基本匹配.注意，由于观测记录由弹性波方程合成，而反演中采用的正演引擎是声波方程Born模拟算子，模拟和观测数据在振幅、波形方面存在一定偏差是正常的.

图13 第20炮模拟与观测数据匹配情况 (a) 正演引擎无法模拟PPS波的方法； (b) 本文方法.Fig.13 Data matching of 20th common-shot gathers (a) The forward modeling engine can not model PPS waves； (b) The proposed method.

3.2 断块模型

通过断块模型合成数据检验方法对复杂介质情况的适用性.图14显示了相应的纵、横波速度与密度模型，其中水深为245 m.这里基于弹性波速度-应力方程Born模拟算子合成含PS波与PPS波的转换波数据，视为观测记录.震源为主频25 Hz的雷克子波，共100炮，炮点以40 m的间隔均匀分布于海面，检波器以5 m的均匀间隔布置于海底，最大炮检距为1.0 km.第50炮转换波记录如图15所示，观测数据中PS波受到PPS波的严重干扰(如15c箭头指示).

图14 断块模型 (a) 纵波速度； (b) 横波速度； (c) 密度.Fig.14 The fault model (a) P-wave velocity; (b) S-wave velocity; (c) Density.

图15 第50炮转换波记录 (a) PS波； (b) 源端水层鸣震相关的PPS波； (c) 合成的转换波观测记录(含PS和PPS波).Fig.15 Converted waves of No.50 common-shot gathers (a) PS waves； (b) PPS waves related to source-side water layer reverberations； (c) The synthetic recordings containing PS and PPS waves.

基于可靠的背景纵波速度(图16a)、海底界面模型以及反射PP波逆时偏移剖面(图16b)进行横波速度反演建模.如图17a所示，初始横波速度模型是真实模型大尺度光滑而成，密度则设为常数.同样进行两组试验，第一组反演算法采用的正演引擎仅模拟PS波，同样发现因为无法与含有PPS波的观测记录匹配，导致反演的横波速度很不合理，尤其在中间断块附近存在明显异常(图17b).本文方法由于同时模拟PS与PPS波，使得转换反射波走时残差估计比较可靠，最终获得了合理的反演结果(图17c).图18显示的目标函数下降曲线表明，本文方法收敛性更好，最终走时残差更小.如图19所示，在水平位置分别为1 km、2 km和3 km处抽取的伪井速度曲线也表明，本文方法合理地考虑了源端水层鸣震对转换反射波走时反演的影响，使横波速度反演过程更稳健，更准确.相反，如果模拟的转换波记录不包含PPS波，反演过程中走时残差最小化过程发生模拟PS波与观测PPS波串相匹配，导致错误的模型更新方向，使得横波速度沿深度交替出现非常异常的剧烈更新(见图19蓝色曲线).图20展示了不同横波速度模型情况下转换PS波逆时偏移结果.由于初始模型误差大，反射界面形态和深度存在明显误差，深层几个散射点对应的绕射波没有收敛.在反演算法中仅模拟PS波导致不合理的走时残差估计进而影响速度更新，使得断块相关的界面以及深层散射点成像存在明显失真.本文方法反演的横波速度模型基本保证所有界面形态和深度的准确性，取得与真实速度模型偏移结果相当的成像效果.

图16 (a) 背景纵波速度； (b) PP波逆时偏移结果Fig.16 (a) Background P-wave velocity； (b) RTM of reflected PP waves

图17 (a) 初始横波速度； (b) 正演引擎无法模拟PPS波时的反演结果； (c) 本文方法反演结果Fig.17 (a) The initial background S-wave velocity； (b) The result of travel-time inversion with a forward modeling engine that can not simulates PS waves； (c) The proposed method

图18 目标泛函下降曲线Fig.18 Normalized descent curve of objective function

图19 伪井横波速度曲线 (a) x=1 km; (b) x=2 km; (c) x=3 km.Fig.19 Pseudo-well curves of the S-wave velocity models at the lateral locations of (a) x=1 km, (b) x=2 km, (c) x=3 km

图20 反射PS波逆时偏移结果 (a) 准确速度模型； (b) 初始横波速度模型； (c) 基于图17b速度模型； (d) 基于图17c速度模型.Fig.20 Reverse time migration of the reflected PS waves The PS image built by (a) true velocity, (b) the initial background S-wave velocity, (c) the S-wave velocity shown in Fig.17b, (d) the S-wave velocity shown in Fig.17c.

4 讨论

本文方法建立在已获取准确海底位置、较好完成多分量数据P/S分离以及获得较准确纵波背景速度和PP波成像结果基础上，这些数据处理结果的准确性直接影响最终横波速度建模精度.海底位置影响源端水层鸣震P波旅行时，在纵波背景速度固定不变情况下，由海底位置不准确引起的鸣震P波走时扰动在反演迭代过程中会逐渐投影到横波背景速度中.此外，纵波背景速度误差同样会迭代投影到横波速度反演结果中，导致横波速度反演结果误差增大.在震源子波频带不变的情况下，PP波成像剖面分辨率直接影响转换波反偏移数据的频带.当PP波成像结果分辨率较低时，相应反偏移转换波数据整体频带会偏低(表现为记录中同相轴变粗)，从而影响DIW走时残差拾取精度，最终影响反演结果.数值实验表明，当横波初始模型离真实值偏差较大时，会影响一次波与多次波的相对走时关系，给基于DIW算法的走时反演带来挑战，导致最终无法获得准确的反演结果.

从正问题描述中发现，我们利用公式(2)和(4)模拟PS波数据，利用公式(2)、(3)、(5)和(6)模拟PPS波数据.其中公式(3)用来描述水层散射场(PP波)，仅需包含海底和水层的速度模型，而不需要全空间模型，所以本文方法数值实验中正演数据模拟阶段计算成本大约为只考虑PS波时的两倍.本文提出的梯度预处理方法需要单独求取PS波、PPS波的梯度.相比于只利用PS波，本文方法求取泛函梯度时计算成本加倍.

从图3中我们发现，当海底较浅时PPS波与PS波路径基本一致，两者照明范围基本相同.随着水深增加PPS波照明范围越来越大，对横波背景速度反演的帮助可能会逐步增大，这一表现后续需做更详细研究.另外，不同水深和海底介质条件源端多次波对横波速度分析的影响不尽相同.结合本文工作，今后还应当从如下几方面开展深入的研究：首先，将本文方法应用到一些浅水探区的实际多分量地震资料，在实践中进行检验与完善，且与对转换波记录进行源端多次波压制，再实施反射走时反演的方法流程进行对比，明确不同技术路线的适用性.其次，虽然不对转换波数据做源端多次波压制时，利用本文方法能反演获得较精确的横波速度模型.但基于现有转换波成像方法，数据中如果含有PPS波等多次波会造成偏移剖面假象严重.所以后续研究面向转换波数据中包含多次波的相关成像方法具有重要实用意义.再次，在硬海底环境下，在海底界面直接透射转换形成的PSS波会使转换波走时反演问题更复杂.尽管可参照本文思路引入针对PSS波的Born模拟并探究相应的梯度预条件技术，但实际可行性有待验证.此外，在深水(水深超过500 m)或超深水(水深超过1500 m)情况下，源端自由表面多次波在水层中传播时间长，主要干扰一些较深层的一次反射与转换波信号.针对深水区海底多分量地震横波速度建模与转换波成像问题开展相关研究也是很有意义的.

5 结论

海底多分量地震数据经过上-下行P/S波分离之后，获得的P波和S波记录仍然受到源端多次波尤其是水层鸣震的干扰，给基于反射PP波和PS波的速度分析和偏移成像造成极大的困扰.近年来，在全波形反演或反射波形反演理论框架下，出现了同时考虑一次反射波与自由表面多次波的反演成像方法，但主要还是面向P波数据，缺乏针对S波数据的相关研究.针对数据域P/S分离之后的转换反射波走时反演问题，本文聚焦浅水、软海底情况，由海底界面地震波能量分配规律推断海底波场分离后的上行S波数据主要受源端水层鸣震(尤其是PPS波)干扰.在反演过程中，如果正演引擎仅考虑PS波，很难与同时包含PS波、PPS波的观测记录进行走时拟合与残差估算，导致不合理的模型更新.为此，本文构建了描述PS波和PPS波的Born模拟方程，利用伴随状态法提出了同时考虑这两种转换波震相的反射走时反演方法，并引入震相分离梯度预条件方法，使横波速度模型的更新方向合理化.水平层状模型和复杂断块模型合成数据实验结果表明，该方法为浅水、软海底环境基于转换波数据的横波速度建模提供了一种有潜力的技术手段.

致谢感谢多位审稿人对本文提出的宝贵建议.

附录A 基于PS波和PPS波合成转换S波的加权系数

源端不同阶次水层鸣震P波产生的下行P波场如附图1 所示.图中红线表示一次透射P波，在S波速度扰动处产生PS反射波.假设源端直达波场能量为1，则其能量为T；蓝线表示源端一阶水层鸣震P波，在S波速度扰动处产生PPS反射波，其能量为-TR；绿线表示源端二阶水层鸣震P波，在S波速度扰动处产生PPPS反射波，其能量为TR2.

附图1 反射S波能量关系示意图 T为下行P波在海底处的透射系数；R为反射系数； -1为自由表面反射系数.Fig.A1 Energy relationship of the reflected S- waves T represents the transmission coefficient of the down-going P- waves at the seabed；R is reflection coefficient；-1 is the reflection coefficient of free surface.

S波速度扰动不变情况下，如果不考虑几何扩散效应，反射波能量大小取决于入射波场能量，所以转换波中任意阶自由表面多次波与一次波能量关系表示为

(A1)

式中，n=1,2,3,…，E表示波场能量求取算子，R表示下行P波在海底界面处的反射系数，可由多分量数据P/S分离(刘学义等,2021)过程中估计出的海底介质纵横波速度和密度计算获得，或利用其他海底参数估计方法(Amundsen and Reitan,1995;Muijs et al.,2004)计算获得.公式(7)中的不同波场成分之间相对幅值需满足公式(A1)，所以将其分别代入公式(A1)求得的wn即为待求加权系数.由于走时反演对模拟数据的振幅精度的要求并不高，根据介质参数估计出海底反射系数之后，结合实际数据情况给出相对合理的加权系数即可.