Elastic parameters waveform inversion method of the second-order time integral wavefield for marine pure P-wave data

WenBo SUN; HongLiang ZHANG; JianHua WANG; GuoChen WU; TieZhao BAO; JiCai DING

doi:10.6038/pg2025HH0540

Progress in Geophysics >

2025 , Vol. 40 >Issue 2: 774 - 786

DOI: https://doi.org/10.6038/pg2025HH0540

Elastic parameters waveform inversion method of the second-order time integral wavefield for marine pure P-wave data

WenBo SUN ^,¹^,² ,
HongLiang ZHANG ^,¹^,²^,* ,
JianHua WANG ¹^,² ,
GuoChen WU ³ ,
TieZhao BAO ¹^,² ,
JiCai DING ¹^,²

Expand

¹ CNOOC Research Institute Ltd., Beijing 100028, China
² National Engineering Research Center of Offshore Oil and Gas Exploration, Beijing 100028, China
³ School of Geosciences, China University of Petroleum, Qingdao 266580, China

Received date: 2024-05-07

Online published: 2025-05-09

Copyright

Fold

Abstract

In marine exploration, when the seismic data in the exploration area only includes P-wave data collected by hydrophones, traditional acoustic full waveform inversion velocity modeling methods cannot invert S-wave velocity information. This paper focuses on ocean fluid-solid coupled media and derives the acoustic-elastic coupled equations using the boundary coupled method of acoustic-elastic wave equations. We validated the potential of using P-wave data for shear wave modeling through numerical simulation results of a three-layer model. An inaccurate velocity model can lead to the failure of full waveform inversion. To address the dependency issue of velocity model in full waveform inversion, we introduce second-order time integration operation and construct a second-order integration objective function. Then, the second-order integral waveform inversion gradient operators suitable for the acoustic-elastic coupled equations were derived, and a P-and S-wave velocity waveform inversion approach based on pure P-wave data in a marine environment was established. Finally, we conduct inversion tests using modified Marmousi2 models of the horizontal and irregular seabed. The inversion results verified the applicability and accuracy of the second-order integrated wave field acoustic-elastic coupled equation inversion method.

Key words： Pure P-wave data; Second-order time integral wavefield; Acoustic-elastic coupled media; Elastic parameters waveform inversion

Cite this article

WenBo SUN , HongLiang ZHANG , JianHua WANG , GuoChen WU , TieZhao BAO , JiCai DING . Elastic parameters waveform inversion method of the second-order time integral wavefield for marine pure P-wave data[J]. Progress in Geophysics, 2025 , 40(2) : 774 -786 . DOI: 10.6038/pg2025HH0540

0 引言

中国海洋油气资源丰富，勘探潜力巨大，在未来具有广阔的勘探开发前景.近年来，随着海洋地震采集技术(OBS)的快速发展，使得获取高精度反演速度可能性增加.在海洋地震勘探工作中，速度是贯穿地震数据采集、处理和解释整个地震勘探过程的重要参数.高精度速度场影响着地震成像、解释和储层预测的可靠性.海洋地震勘探环境可以抽象为上覆海水介质和下伏弹性介质的流-固耦合介质(FSCM，Filuid-Solid Coupled Media)，采用更加符合实际地震波传播的控制方程有助于提高正演模拟精度，进而获得更加精确的反演结果.对此，Hou等(2012)和De Basabe和Sen(2015)将地震波在FSCM下的数值模拟方法分为两种：(1)整体控制方程法，即在流体和固体中采用相同的地震波波动方程进行地震波数值模拟；(2)双区控制方程法，即分别用声波方程和弹性波方程来描述地震波在液体和固体介质中的传播，并在流固边界处利用连续性条件将两者介质联系起来.Stephen(1983, 1988)和Sochacki等(1991)利用纯位移的声学和弹性波方程，结合流固界面的应力、应变连续性条件，采用有限差分法进行数值模拟研究；Yu等(2016)将声波方程代入弹性波方程，采用交错网格法完成了声弹耦合方程的逆时偏移成像；Wu和Wu(2017)提出基于高阶有限差分法的一阶速度-应力声波和弹性波方程耦合算法，并模拟验证了地震波在流-固耦合介质中传播的稳定性；Li等(2019)采用过渡层策略，对不规则海底界面流-固耦合介质进行有限差分地震正演模拟；Zhang等(2023)考虑固相中存在的各向异性，将声-弹耦合波动方程从各向同性介质扩展至VTI介质，并进行交错网格有限差分数值模拟；Yang等(2023)考虑温度效应，建立流-固耦合条件下热弹性介质中地震波的传播方程.

全波形反演(FWI)是以波动方程为基础，通过最小化模拟和观测数据之间的残差，从而获得高分辨率的地下物理参数(如速度、密度等)的方法.常规全波形反演对初始速度模型的依赖性强，不准确的初始速度模型会导致基于局部寻优算法的收敛陷入局部极值.为了降低这种依赖性，Bunks等(1995)提出了多尺度反演策略；Pratt等(Pratt等，1998；Pratt，1999)提出在频率域中频率由低到高的分频反演策略.在上述反演策略思想的影响下，许多学者提出了一些基于全波形反演框架下的低频背景速度的反演方法.如Laplace域、Laplace-Fourier域、对数域全波形反演方法(Shin and Cha, 2008, 2009; Shin and Min, 2006)，反射波形反演(Xu et al., 2012；孙思宇等，2021)，积分类全波形反演方法(Donno et al., 2013; 陈生昌和陈国新, 2016), 基于包络信号的全波形反演方法(Wu et al., 2014; Luo and Wu, 2015; 梁展源等, 2019; 李冀蜀等, 2023; Xiong et al., 2023)，基于自适应滤波器的波形反演方法(Warner and Guasch, 2016)，波动方程全走时反演(吴彦等, 2017)等.在低频背景模型恢复的基础上在进行全波形反演，一定程度上降低了目标函数落入局部极值的风险.

本文首先分析海上纵波资料存在的转换波信息，以说明纯纵波资料可以指示海底以下弹性层介质的潜力.基于适用于流-固耦合介质的声-弹耦合波动方程，构建全波形反演二阶积分波场目标函数，进而计算梯度算子，通过迭代反演获取低频大尺度速度模型，并以反演出来的低频速度模型作为初始速度模型进行常规全波形反演，更新模型的中高频信息.最后采用更加符合实际的重采样的水平海底界面的Marmousi2模型和修改后的起伏海底界面的Marmousi2模型进行反演测试，验证了本文方法利用海上纯纵波资料恢复弹性参数波形反演方法的适用性.

1 方法原理

1.1 基于上覆海水介质的声-弹耦合方程

海洋环境中上层流体中剪切应力为零，地震波在此区域的传播可以使用一阶速度压力声波方程进行描述：

(1)

$\begin{equation*}\boldsymbol{L}_{\text {acoustic }} \boldsymbol{T}=\boldsymbol{f}_{\mathrm{p}}, \end{equation*}$

其中$\boldsymbol{T}=\left[p, v_{x}, v_{z}\right]^{\mathrm{T}}$，p是流体压力，v_x和v_z为纵、横向速度矢量，$\boldsymbol{f}_{\mathrm{p}}=\left[f_{\mathrm{p}}, 0, 0\right]^{\mathrm{T}}$是震源项，L是声学偏微分算子矩阵：

(2)

$\boldsymbol{L}_{\text {acoustic }}=\left[\begin{array}{ccc}\partial / \partial t & -\rho v_{\mathrm{P}}^{2} \partial / \partial x & -\rho v_{\mathrm{P}}^{2} \partial / \partial z \\-\partial / \partial x & \rho \partial / \partial t & 0 \\0 & -\partial / \partial z & \rho \partial / \partial t\end{array}\right] .$

地震波在海洋中激发，穿越海底流固分界面后以弹性波的形式传播，故在海洋固相介质中地震波是无源的，地震波在该区域的传播可由一阶速度-应力弹性波方程来描述，

(3)

$\boldsymbol{L}_{\text {elastic }} \boldsymbol{U}=0, $

其中$\boldsymbol{U}=\left(v_x, v_z, \tau_{x x}, \tau_{z z}, \tau_{x z}\right)^{\mathrm{T}}$，L_elastic是弹性偏微分算子矩阵:

(4)

$\boldsymbol{L}_{\text {elastic }}=\left[\begin{array}{ccccc}\rho \partial / \partial t & 0 & -\partial / \partial x & 0 & -\partial / \partial z \\0 & \rho \partial / \partial t & 0 & -\partial / \partial z & -\partial / \partial x \\-\rho v_{\mathrm{p}}^2 \partial / \partial x & -\rho\left(v_{\mathrm{P}}^2-2 v_{\mathrm{s}}^2\right) \partial / \partial z & \partial / \partial t & 0 & 0 \\-\rho\left(v_{\mathrm{P}}^2-2 v_{\mathrm{S}}^2\right) \partial / \partial x & -\rho v_{\mathrm{P}}^2 \partial / \partial z & 0 & \partial / \partial t & 0 \\-\rho v_{\mathrm{S}}^2 \partial / \partial z & -\rho v_{\mathrm{S}}^2 \partial / \partial x & 0 & 0 & \partial / \partial t\end{array}\right] \text {, }$

为了连接海底界面的两个控制方程，应考虑流-固界面上应力和速度的法向分量的连续性，在只考虑水平界面的情况下，流-固界面处的应力和质点速度(位移)连续条件如下：

(5)

$\left\{\begin{array}{l}\tau_{x x}=\tau_{z z}=-P \\\tau_{x z}=0\end{array}\right.,$

(6)

$\boldsymbol{v}_{\text {acoustic }}(\boldsymbol{x}, t)=\boldsymbol{v}_{\text {elastic }}(\boldsymbol{x}, t) \text {. }$

结合上述公式(1)—(5)，得声-弹耦合方程：

(7)

$\boldsymbol{L}_{\text {couple }} \boldsymbol{\varPhi}=\boldsymbol{F}，$

其中$\boldsymbol{\varPhi}=\left(P, v_x, v_z, \tau_{x x}, \tau_{z z}, \tau_{x z}\right)^{\mathrm{T}}$，震源F =[ f_p 0 0 0 0 0 ]^T在流体中激发，其形式为：

(8)

$\boldsymbol{L}_{\text {couple }}=\left[\begin{array}{cccccc}\partial / \partial t & \rho\left(v_{\mathrm{P}}^2-v_{\mathrm{S}}^2\right) \partial / \partial x & \rho\left(v_{\mathrm{P}}^2-v_{\mathrm{S}}^2\right) \partial / \partial z & 0 & 0 & 0 \\\partial / \partial x & \rho \partial / \partial t & 0 & -\partial / \partial x & 0 & -\partial / \partial z \\\partial / \partial z & 0 & \rho \partial / \partial t & 0 & -\partial / \partial z & -\partial / \partial x \\0 & -\rho v_{\mathrm{S}}^2 \partial / \partial x & \rho v_{\mathrm{S}}^2 \partial / \partial z & \partial / \partial t & 0 & 0 \\0 & \rho v_{\mathrm{S}}^2 \partial / \partial x & -\rho v_{\mathrm{S}}^2 \partial / \partial z & 0 & \partial / \partial t & 0 \\0 & -\rho v_{\mathrm{S}}^2 \partial / \partial z & -\rho v_{\mathrm{S}}^2 \partial / \partial x & 0 & 0 & \partial / \partial t\end{array}\right] .$

1.2 常规声-弹耦合方程全波形反演方法

常规全波形反演通过观测和模拟数据间的差异构建最小二乘目标函数：

(9)

$\begin{equation*}E(\boldsymbol{m})=\frac{1}{2} \sum\limits_{s, r} \int_{t}\left[d_{\text {obs }}\left(\boldsymbol{x}_{\mathrm{s}}, \boldsymbol{x}_{\mathrm{r}}, t\right)-d_{\text {cal }}\left(\boldsymbol{x}_{\mathrm{s}}, \boldsymbol{x}_{\mathrm{r}}, t\right)\right]^{2} \mathrm{~d} t ,\end{equation*}$

其中，m是模型参数，s表示炮点位置，r表示检波点位置，d_obs为纵波观测数据，d_cal为纵波模拟数据.

当采取梯度类算法进行求解时，梯度为目标函数对模型参数的一阶导数，即：

(10)

$\begin{align*}\nabla E(\boldsymbol{m})= & \frac{1}{2} \sum\limits_{s, r} \int_{t} \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{m}} *\left[d_{\text {obs }}\left(\boldsymbol{x}_{\mathrm{s}}, \boldsymbol{x}_{\mathrm{r}}, t\right)-\right. \\& \left.d_{\text {cal }}\left(\boldsymbol{x}_{\mathrm{s}}, \boldsymbol{x}_{\mathrm{r}}, t\right)\right] \mathrm{d} t, \end{align*}$

其中，$\frac{\partial \boldsymbol{u}}{\partial \boldsymbol{m}}$为雅可比矩阵，Li和Wu(2020)通过一阶Born近似推导得到基于声-弹耦合方程全波形反演目标函数对纵波速度、横波速度，密度梯度算子表达式为：

(11)

$\begin{align*}\frac{\partial E}{\partial v_{\mathrm{p}}}= & \int\limits_{t}-2 \rho v_{\mathrm{P}} P^{*}\left(\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}\right) \mathrm{d} t \\\frac{\partial E}{\partial v_{\mathrm{s}}}= & \int\limits_{t} 2 \rho v_{\mathrm{S}}\left[P^{*}\left(\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}\right)+\tau_{x z}^{*}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)+\right. \\& \left.\tau_{x x}^{*}\left(\frac{\partial u}{\partial x}-\frac{\partial w}{\partial z}\right)+\tau_{z z}^{*}\left(\frac{\partial w}{\partial z}-\frac{\partial u}{\partial x}\right)\right] \\\frac{\partial E}{\partial \rho}= & \int\limits_{t}-\left(u^{*} \frac{\partial^{2} u}{\partial t^{2}}+w^{*} \frac{\partial^{2} w}{\partial t^{2}}\right) \mathrm{d} t ，\end{align*}$

引入算子J表示雅可比矩阵，向量η表示观测数据与模拟数据间的残差，则公式(11)可简化为：

(12)

$\begin{equation*}\nabla E=\boldsymbol{J}^{\mathrm{T}} \boldsymbol{\eta} .\end{equation*}$

1.3 基于声-弹耦合方程二阶积分波场全波形反演方法

声-弹耦合波动方程反演使用L₂模建立正演模拟数据和观测数据之间的积分波场目标函数为：

(13)

$\begin{array}{l}{E_Q}({\boldsymbol{m}}) = \frac{1}{2}\sum\limits_{s,r} {} {\int\limits_t}|Q\left[ {{d_{{\rm{obs }}}}\left( {{{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{x}}_{\rm{r}}},t} \right)} \right] - \\~~~~~~~~~~~~~~~~~~Q\left[ {{d_{{\rm{cal}}}}\left( {{{\boldsymbol{x}}_s},{{\boldsymbol{x}}_r},t} \right)} \right]{|^2}{\rm{d}}t，\end{array}$

其中Q为二阶时间积分运算符，其满足：

(14)

$\begin{equation*}\mathrm{Q}[\boldsymbol{u}(\boldsymbol{x}, t)]=F^{-1}\left(\frac{\tilde{\boldsymbol{u}}(\boldsymbol{x}, \omega)}{(-\mathrm{i} \omega)^{2}}\right), \end{equation*}$

其中F和F^－1分别表示傅里叶正、逆变换，ω为角频率.

本文中密度ρ均为常密度，根据一阶Born近似假设，声弹耦合介质可分为背景介质(v_P, v_S)和扰动介质(Δv_P, Δv_S).将背景和扰动参数代入式(7)可得：

(15)

$\begin{equation*}\boldsymbol{L}_{\text {couple }}\left(v_{\mathrm{P}}+\delta v_{\mathrm{P}}, v_{\mathrm{S}}+\delta v_{\mathrm{S}}\right)(\boldsymbol{\varPhi}+\Delta \boldsymbol{\varPhi})=0, \end{equation*}$

其中$\Delta \boldsymbol{\varPhi}=\left(\delta P, \delta v_{x}, \delta v_{z}, \delta \tau_{x x}, \delta \tau_{z z}, \delta \tau_{x z}\right)^{\mathrm{T}}$为波场扰动.式(7)减去背景介质产生得背景波场得：

(16)

$\begin{equation*}\boldsymbol{L}_{\text {couple }} \Delta \boldsymbol{\varPhi}=-\boldsymbol{L}_{v_{\mathrm{P}}}^{\prime} \boldsymbol{\varPhi} \delta v_{\mathrm{P}}-L_{v_{\mathrm{S}}}^{\prime} \boldsymbol{\varPhi} \delta v_{\mathrm{S}} ,\end{equation*}$

其中$\boldsymbol{L}_{v_{\mathrm{P}}, ~}^{\prime}, ~ \boldsymbol{L}_{v_{\mathrm{S}}}^{\prime}$表示L相对于弹性参数(v_P, v_S)的偏导数，并表示为：

(17)

$\boldsymbol{L}_{v_{\mathrm{P}}}^{\prime}=\left[\begin{array}{cccccc}0 & \partial / \partial x & \partial / \partial z & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0\end{array}\right],$

(18)

$\boldsymbol{L}_{v_{\mathrm{S}}}^{\prime}=\left[\begin{array}{cccccc}0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & -\partial / \partial x & \partial / \partial z & 0 & 0 & 0 \\0 & \partial / \partial x & -\partial / \partial z & 0 & 0 & 0 \\0 & -\partial / \partial z & -\partial / \partial x & 0 & 0 & 0\end{array}\right] .$

根据式(16)，可进一步推导出：

(19)

$\begin{equation*}\frac{\partial \boldsymbol{\varPhi}}{\partial \boldsymbol{m}}=-\boldsymbol{L}_{\text {couple }}{}^{-1} ~\boldsymbol{L}_{\mathrm{m}}^{\prime} \boldsymbol{\varPhi} ，\end{equation*}$

然后用积分目标函数(式(13))关于弹性参数的偏导数来计算梯度：

(20)

$\begin{align*}\nabla_{\mathrm{m}} E_{Q}= & \int\limits_{0}^{t_{\max }}-Q\left(\boldsymbol{\varPhi}^{\mathrm{T}}\right)\left(\boldsymbol{L}_{\mathrm{m}}^{\prime}\right)^{\mathrm{T}}\left(\boldsymbol{L}_{\text {couple }}{}^{-1}\right)^{\mathrm{T}} \\& {\left[Q\left(P_{\text {obs }}\right)-Q\left(P_{\text {cal }}\right)\right] \mathrm{d} t } ，\end{align*}$

给出伴随波场ψ计算方程如下：

(21)

$\begin{equation*}\left(\boldsymbol{L}_{\text {couple }}{}^{-1}\right)^{\mathrm{T}} \boldsymbol{\psi}=Q\left(P_{\text {obs }}\right)-Q\left(P_{\text {cal }}\right), \end{equation*}$

引入伴随波场将梯度改写为：

(22)

$\left\{\begin{array}{l}\nabla_{v_{\mathrm{P}}} E_{Q}=-\int_{0}^{t_{\max }} Q\left(\boldsymbol{\varPhi}^{\mathrm{T}}\right)\left(\boldsymbol{L}_{v_{\mathrm{P}}}^{\prime}\right)^{\mathrm{T}} \boldsymbol{\psi} \mathrm{~d} t \\\nabla_{v_{\mathrm{S}}} E_{Q}=-\int_{0}^{t_{\max }} Q\left(\boldsymbol{\varPhi}^{\mathrm{T}}\right)\left(\boldsymbol{L}_{v_{\mathrm{S}}}^{\prime}\right)^{\mathrm{T}} \boldsymbol{\psi} \mathrm{~d} t\end{array}\right.，$

其中$\boldsymbol{\psi}=\left(\psi_{\mathrm{P}}, \psi_{x}, \psi_{z}, \psi_{x x}, \psi_{z z}, \psi_{x z}\right)^{\mathrm{T}}$表示向后伴随波场.

2 数值测试

2.1 Ricker子波时间积分

我们假设波场u(t)一个主频为10 Hz的Ricker子波，如图 1黑线所示，对Ricker子波进行归一化的时间一阶积分，如图 1蓝线所示，其对应的频谱如图 2中蓝线所示；对Ricker子波进行归一化的时间二阶积分，如图 1红线所示，其对应的频谱如图 2中红线所示.综合图 1的时间信号和图 2的频谱可以看出，随着对原始Ricker进行归一化积分运算，积分运算后的Ricker子波相较于原始子波出现较为明显的低频化，且随着积分阶数的增加，低频分量占比逐渐提高.

显示原图|下载原图ZIP|生成PPT

图1 Ricker子波及其时间积分(主频为10 Hz)

黑线为Ricker子波；蓝线为Ricker子波的时间一阶积分；红线为Ricker子波的时间二阶积分.

Fig 1 Ricker wavelet with a dominant frequency of 10 Hz and its time integral

The black line represents the Ricker wavelet, the blue line denotes the first-order time integral of Ricker wavelet and the red line denotes the second-order time integral of the Ricker wavelet.

显示原图|下载原图ZIP|生成PPT

图2 Ricker子波及其时间积分的频谱(主频为10 Hz)

黑线为Ricker子波的频谱；蓝线为Ricker子波时间一阶积分的频谱；红线为Ricker子波时间二阶积分的频谱.

Fig 2 The frequency spectrum of a Ricker wavelet with a dominant frequency of 10 Hz and its time integral

The black line represents the spectrum of the Ricker wavelet, the blue line denotes the spectrum of first-order time integral of ricker wavelet, the red line denotes the spectrum of second-order time integral of ricker wavelet.

2.2 流-固耦合介质三层模型

为验证分析海上声-弹耦合方程正演模拟的准确性和稳定性，本节将进行一个流-固边界耦合的三层模型的数值模拟.我们设计一个301×151网格大小的层状模型，如图 3所示，炮点位于坐标(600 m，10 m)处，空间步长均为10 m，时间采样间隔为0.8 ms，震源选择主频是18 Hz雷克子波.介质为三层介质模型，第一层纵波速度1500 m/s，横波速度0 m/s；第二层纵波速度2500 m/s，横波速度1470 m/s；第三层纵波速度3000 m/s，横波速度1764 m/s，密度均为常数.

显示原图|下载原图ZIP|生成PPT

图3 三层弹性介质模型示意图

(a)纵波速度模型；(b)横波速度模型.

Fig 3 Elastic three-layered model

(a) Longitudinal wave velocity model; (b) Shear wave velocity model.

图 4a和图 4b分别为用声波方程和声-弹耦合方程正演模拟得到的接收记录，图 4a为纯声波记录，主要观测到的波有①P₀^T直达波、②P₀^TP₁^RP₀^T反射波、③P₀^TP₁^TP₂^RP₁^TP₀^T反射波和④P₀^TP₁^TP₀^T第一层界面处的滑行波四种波，此外层间也会产生多次波，由于振幅弱，因此在记录中显示不明显.图 4b为声-弹耦合波动方程正演模拟得到的记录，波形依次为①P₀^T直达波、②P₀^TP₁^RP₀^T反射波、③P₀^TP₁^TP₂^RP₁^TP₀^T反射波、④P₀^TP₁^TP₀^T第一层界面处的滑行波、⑤P₀^TP₁^TS₂^RP₁^TP₀^T/P₀^TS₁^TP₂^RP₁^TP₀^T反射波和⑥P₀^TS₁^TS₂^RP₁^TP₀^T反射波，此外，记录中还存在诸多层间多次反射的转换波，因其振幅能量弱，故其在观测记录中并不清晰. 上述波中P和S表示波的类型(纵波和横波)，T表示透射波，R表示反射波，下标表示波在第几层传播.在上述模型测试中，声-弹耦合方程接收到的单分量记录中包含转换横波信息，验证了基于声-弹耦合方程的纯纵波资料反演横波速度的潜力.

显示原图|下载原图ZIP|生成PPT

图4 (a) 声波方程模拟记录；(b)声-弹耦合方程模拟记录

Fig 4 (a) Seismic records simulated by acoustic equation; (b) Seismic records simulated by acoustic-elastic coupled equation

2.3 水平海底界面的Marmousi2模型

为了验证基于声-弹耦合波动方程的时间二阶积分波场波形反演的正确性和有效性.本节重采样后修改Marmousi2模型，在模型顶端设置一层海水层，海水层纵波速度设置为1500 m/s，横波速度设置为0 m/s.Marmousi2模型的网格化参数为nx=465，nz=148，dx=10 m，dz=10 m.图 5a为Marmousi2纵波速度模型，图 5b为Marmousi2横波速度模型，其中横波速度是由纵波速度计算得到(V_s=V_p/1.732).反演测试在海水层中震源共激发65次，炮间距70 m，每炮465道，道间距10 m.总时间采样点数为2500，时间采样间隔为0.8 ms.震源采用主频为8 Hz的雷克子波.图 6a和图 6b分别为线性化初始纵波速度模型和横波速度模型，即采用真实模型的最大速度和最小速度，设置梯度进行线性平滑.

显示原图|下载原图ZIP|生成PPT

图5 Marmousi2模型

(a)纵波速度；(b)横波速度.

Fig 5 Marmousi2

(a) P-wave velocity; (b) S-wave velocity.

显示原图|下载原图ZIP|生成PPT

图6 Marmousi2模型

(a)线性初始纵波速度；(b)线性初始横波速度.

Fig 6 Marmousi2 model

(a)Linear initial P-wave velocity; (b) Linear initial S-wave velocity.

图 7a和图 7b是声-弹耦合波方程时间二阶积分波场弹性参数波形反演(简称IFWI)迭代60次结果.由结果可见，全波形反演纵横波速度反演结果大致构造轮廓基本正确，由于反演过程中的积分运算相当于低通滤波，所以反演的分辨率不高.我们将图 7所示模型作为初始速度，进行全波形反演(简称FWI).图 8是以图 7模型为初始速度的常规全波形反演结果.图 9a是真实纵波接收记录，图 9b—d为不同迭代次数模型的数值模拟记录结果，对比可以看出反演逐步向真实目标记录贴近.作为比较，我们采用图 6的线性初始速度模型作为常规全波形反演的初始速度模型，图 10为以图 6所示线性初始速度模型为初始速度的常规反演速度的结果.对比图 8和图 10的反演结果，基于声-弹耦合波方程时间二阶积分波场的弹性参数波形反演和常规全波形反演的逐步反演结果，构造轮廓清晰，速度更加接近真实模型.

显示原图|下载原图ZIP|生成PPT

图7 时间二阶积分波场全波形反演结果

(a)反演纵波速度；(b)反演横波速度.

Fig 7 The inversion results by full waveform inversion of the second-order time integral wavefield

(a) Inverted S-wave velocity; (b) Inverted S-wave velocity.

显示原图|下载原图ZIP|生成PPT

图8 时间二阶积分波场+全波形反演(IntI₂+FWI)结果

(a)反演纵波速度；(b)反演横波速度.

Fig 8 The second-order time integral wavefield+full waveform inversion (IntI₂+FWI) results

(a) Inverted P-wave velocity; (b) Inverted S-wave velocity.

显示原图|下载原图ZIP|生成PPT

图9 地震记录

(a)观测记录；(b)初始模型模拟记录；(c)积分波场反演记录；(d)积分波场+全波形反演记录.

Fig 9 Seismic records

(a) Observed data; (b) Record calculated using initial model; (c) Record calculated using integral wavefield; (d) Record calculated final inverted model.

显示原图|下载原图ZIP|生成PPT

图10 以线性模型为初始速度的常规全波形反演结果

(a)反演纵波速度；(b)反演横波速度.

Fig 10 The results of a single full waveform inversion with a linear model as the initial velocity

(a) Inverted P-wave velocity; (b) Inverted S-wave velocity.

为了更好的对比两种反演结果的细节，图 11分别为Marmousi模型两种反演结果的第1500 m、2450 m、4000 m处深度-速度(纵波速度、横波速度)曲线对比图.图中黑线为真实速度，蓝线为线性初始速度，绿线为单纯的常规全波形反演速度，红线为二阶积分波场全波形反演加常规全波形反演结果.由图 11可以明显看出，二阶积分波场全波形反演加常规全波形反演结果反演曲线更加贴合实际速度曲线，验证了在只输入纯纵波资料的情况下，基于声-弹耦合方程的二阶积分波场全波形反演加常规全波形反演反法能够很好的重建出真实纵、横波速度模型.

显示原图|下载原图ZIP|生成PPT

图11 第1500 m、2450 m、4000 m道(左，中，右)的两种反演结果的(a)深度-纵波速度和(b)深度-横波速度曲线对比

图中黑线为真实速度，绿线为初始速度，蓝线为常规全波形反演结果，红线为时间二阶积分+全波形反演的结果.

Fig 11 Comparison of (a) depth-P-wave velocity and (b) depth-S-wave velocity curves for two inversion results at traces of 1500 m, 2450 m, and 4000 m (left, middle, right)

The black line in the figure represents the true velocity, the green line represents the initial velocity, the blue line represents the conventional full waveform inversion result, and the red line represents the time second-order integration+full waveform inversion result.

2.4 起伏海底界面的Marmousi2模型

实际海底界面并非理想的水平界面，因此我们选择不规则海底界面的Marmousi2模型进行反演用以验证基于声-弹耦合波动方程的时间二阶积分波场波形反演的正确性和有效性.顶端海水层纵波速度依然设置为1500 m/s，横波速度设置为0 m/s，其他反演参数与2.3节相同.图 12a和图 12b分别为真实不规则海底界面Marmousi2纵波速度模型和横波速度模型，图 13a和图 13b分别为真实模型线性化后的初始纵波速度模型和横波速度模型.

显示原图|下载原图ZIP|生成PPT

图12 Marmousi2模型

(a)纵波速度；(b)横波速度.

Fig 12 Marmousi2

(a) P-wave velocity; (b) S-wave velocity.

显示原图|下载原图ZIP|生成PPT

图13 Marmousi2模型

(a)线性初始纵波速度；(b)线性初始横波速度.

Fig 13 Marmousi2 model

(a) Linear initial P-wave velocity; (b) Linear initial S-wave velocity.

图 14a和图 14b是声-弹耦合波方程时间二阶积分波场弹性参数波形反演(简称IFWI)迭代60次结果.由结果可见，全波形反演纵横波速度反演结果大致构造轮廓基本正确，由于本文正演所用网格为矩形网格，在起伏海底界面处会产生“阶梯状”非物理散射噪声，致使反演结果界面处存在相应的噪声.为消除这种虚假散射，许多学者基于有限差分法给出解决方案，如贴体网格法(蒋丽丽, 2010)、垂向网格映射法(Qu et al., 2020)、广义有限差分法(贾宗峰等, 2022；Li et al., 2023)等，本文不作过多阐述.图 15是以图 14模型为初始速度的常规全波形反演结果，由图可得，出了起伏界面附近的噪声，其他区域轮廓清晰，构造清楚，基本验证了本文反演方法对起伏海底界面的适用性.

显示原图|下载原图ZIP|生成PPT

图14 时间二阶积分波场全波形反演结果

(a)反演纵波速度；(b)反演横波速度.

Fig 14 The inversion results by full waveform inversion of the second-order time integral wavefield

(a) Inverted S-wave velocity; (b) Inverted S-wave velocity.

显示原图|下载原图ZIP|生成PPT

图15 时间二阶积分波场+全波形反演(IntI₂+FWI)结果

(a)反演纵波速度；(b)反演横波速度.

Fig 15 The second-order time integral wavefield+full waveform inversion (IntI₂+FWI) results (a) Inverted P-wave velocity; (b) Inverted S-wave velocity.

为了更好的对比两种反演结果的细节，图 16分别为起伏海底界面Marmousi模型的IFWI反演结果和IFWI+FWI反演结果的第1500 m、2450 m、4000 m处深度-速度(纵波速度、横波速度)曲线对比图.图中黑线为真实速度，蓝线为线性初始速度，绿线为二阶积分波场全波形反演速度，红线为二阶积分波场全波形反演加常规全波形反演结果.由图 16可以明显看出，二阶积分波场全波形反演加常规全波形反演结果反演曲线基本贴近真实模型曲线，验证了在只输入海上纯纵波资料的情况下，本文方法对起伏海底界面的适用性.

显示原图|下载原图ZIP|生成PPT

图16 第1500 m、2450 m、4000 m道(左，中，右)的起伏海底界面反演结果的(a)深度-纵波速度和(b)深度-横波速度曲线对比

图中黑线为真实速度，绿线为初始速度，蓝线为常规全波形反演结果，红线为时间二阶积分+全波形反演的结果.

Fig 16 Comparison of the irregular seabed interface inversion results for the 1500 m, 2450 m, and 4000 m traces (left, middle, right) in terms of (a) depth-P-wave velocity, and (b) depth-S-wave velocity curves

3 结论

本文针对海洋勘探环境，提出了声波和弹性波方程边界耦合方法用于模拟地震波在海水中激发向下穿越流体和弹性固体后被水听器接收的传播过程，并建立海上纯纵波资料的横波初始速度建模方法.声-弹耦合方程模拟的地震记录包含丰富转换波信息，为只采用纯纵波资料反演横波速度提供支撑.针对全波形反演对初始速度模型的强依赖性的问题，本文构建了基于声-弹耦合方程的二阶积分目标函数，并验证了积分计算有效的增强了低频成分，复杂Marmousi2模型测试结果表明基于声-弹耦合方程的二阶积分弹性波形反演能够有效的恢复纵横波速度模型的低频成分，以低频恢复的反演结果作为下一步常规波形反演的初始速度，最终良好的反演结果验证了声-弹耦合方程的二阶积分弹性波形反演方法的适用性和准确性.

感谢评审专家提出的宝贵意见.

References

Publishing order | Descend order by publishing year | Descend order by cited within

Bunks

, Saleck

F M

, Zaleski

, et al. Multiscale seismic waveform inversion. Geophysics, 1995, 60 (5): 1457- 1473.

模态框（Modal）标题

Abstract

Cite this article

0 引言

1 方法原理

1.1 基于上覆海水介质的声-弹耦合方程

1.2 常规声-弹耦合方程全波形反演方法

1.3 基于声-弹耦合方程二阶积分波场全波形反演方法

2 数值测试

2.1 Ricker子波时间积分

图1 Ricker子波及其时间积分(主频为10 Hz)

图2 Ricker子波及其时间积分的频谱(主频为10 Hz)

2.2 流-固耦合介质三层模型

图3 三层弹性介质模型示意图

图4 (a) 声波方程模拟记录；(b)声-弹耦合方程模拟记录

2.3 水平海底界面的Marmousi2模型

图5 Marmousi2模型

图6 Marmousi2模型

图7 时间二阶积分波场全波形反演结果

图8 时间二阶积分波场+全波形反演(IntI2+FWI)结果

图9 地震记录

图10 以线性模型为初始速度的常规全波形反演结果

图11 第1500 m、2450 m、4000 m道(左，中，右)的两种反演结果的(a)深度-纵波速度和(b)深度-横波速度曲线对比

2.4 起伏海底界面的Marmousi2模型

图12 Marmousi2模型

图13 Marmousi2模型

图14 时间二阶积分波场全波形反演结果

图15 时间二阶积分波场+全波形反演(IntI2+FWI)结果

图16 第1500 m、2450 m、4000 m道(左，中，右)的起伏海底界面反演结果的(a)深度-纵波速度和(b)深度-横波速度曲线对比

3 结论

References

图8 时间二阶积分波场+全波形反演(IntI₂+FWI)结果

图15 时间二阶积分波场+全波形反演(IntI₂+FWI)结果