CN104965223A - 粘声波全波形反演方法及装置 - Google Patents

Abstract

本发明公开一种粘声波全波形反演方法及装置，该方法包括：根据由多个Maxwell体构成的GSLS模型的粘弹性波动方程，确定时间域一阶速度-应力粘声波方程；根据时间域一阶速度-应力粘声波方程和二范数目标函数，确定该方程相应的松弛模量梯度计算公式及纵波速度梯度计算公式；根据时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场；根据正向传播波场、残差逆时反传波场、松弛模量梯度计算公式及纵波速度梯度计算公式，确定松弛模量梯度和纵波速度梯度；根据纵波速度梯度采用局部最优化方法进行二维粘声波全波形反演。可以更好地近似地球介质在地震频带范围内的常Q特征，达到良好的数值频散抑制能力。

Description

粘声波全波形反演方法及装置

技术领域

本发明涉及石油勘探地震速度建模技术领域，尤其涉及粘声波全波形反演方法及装置。

背景技术

理想地球介质为弹性介质，而实际地球介质本身具有粘滞性，导致地震波在传播中的振幅和频率受到严重影响，尤其对于近地表强衰减层等强衰减区域来说，如果不考虑地球介质的吸收衰减作用，正演地震波场与实际观测波场之间的差异较大，造成反演结果产生较大误差。因此，研究地震波在粘弹性介质中传播规律并进行粘弹性介质全波形反演方法研究具有重要意义。

大量学者通过粘弹性介质衰减机理及衰减规律研究，出现了不同的粘弹性模型，同时，研究表明地球介质在地震频带范围内具有近似常Q特征，而广义标准线性固体(generalizedstandard linear solid，GSLS)模型可以很好地近似这种常Q特征。此外，由于GSLS模型相应的波动方程具有在时间域易于模拟的优点，因此基于GSLS模型的粘弹性或粘声波方程被广泛地用于正演模拟、逆时偏移及粘声波全波形反演研究中。Bai等(2014)基于单个标准线性固体的粘弹性模型，采用以位移表示的二阶粘声波方程和中心网格有限差分法，成功地进行了数值模型及实际资料的全波形反演。

Bai等进行粘声波全波形反演时，采用的是基于单个Maxwell(麦克斯韦)体构成的GSLS模型，但研究表明由单个Maxwell体构成的GSLS模型不足以近似地球介质在地震频带范围内的常Q特征。此外，Bai等进行粘声波全波形反演时，采用二阶位移粘声波方程和中心网格有限差分法，要降低数值频散的影响就需要增加有限差分阶数，从而会增加计算量。

发明内容

本发明实施例提供一种粘声波全波形反演方法，用以更好地近似地球介质在地震频带范围内的常Q特征，并达到良好的数值频散抑制能力，该方法包括：

根据GSLS模型的粘弹性波动方程，确定时间域一阶速度-应力粘声波方程，其中，所述GSLS模型由多个Maxwell体构成；

根据所述时间域一阶速度-应力粘声波方程和二范数目标函数，确定所述时间域一阶速度-应力粘声波方程相应的松弛模量梯度计算公式及纵波速度梯度计算公式；

根据所述时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场；

根据所述正向传播波场、残差逆时反传波场、松弛模量梯度计算公式及纵波速度梯度计算公式，确定松弛模量梯度和纵波速度梯度；

根据所述纵波速度梯度，采用局部最优化方法进行二维粘声波全波形反演。

本发明实施例还提供一种粘声波全波形反演装置，用以更好地近似地球介质在地震频带范围内的常Q特征，并达到良好的数值频散抑制能力，该装置包括：

粘声波方程确定模块，用于根据GSLS模型的粘弹性波动方程，确定时间域一阶速度-应力粘声波方程，其中，所述GSLS模型由多个Maxwell体构成；

计算公式确定模块，用于根据所述时间域一阶速度-应力粘声波方程和二范数目标函数，确定所述时间域一阶速度-应力粘声波方程相应的松弛模量梯度计算公式及纵波速度梯度计算公式；

波场确定模块，用于根据所述时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场；

梯度确定模块，用于根据所述正向传播波场、残差逆时反传波场、松弛模量梯度计算公式及纵波速度梯度计算公式，确定松弛模量梯度和纵波速度梯度；

全波形反演模块，用于根据所述纵波速度梯度，采用局部最优化方法进行二维粘声波全波形反演。

本发明实施例中，采用由多个Maxwell体构成的GSLS模型进行粘声波全波形反演，相对于现有技术采用由单个Maxwell体构成的GSLS模型进行粘声波全波形反演而言，可以更好地近似地球介质在地震频带范围内的常Q特征；并且，在进行粘声波全波形反演时，采用时间域一阶速度-应力粘声波方程和高阶交错网格有限差分法，这样在差分阶数相同时，能够比现有技术的中心网格有限差分法有更好的数值频散抑制能力。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：

图1为本发明实施例中粘声波全波形反演方法的示意图；

图2为本发明实施例中粘声波全波形反演装置的示意图；

图3为本发明实施例中未考虑粘弹性影响的声波全波形反演结果示意图；

图4为本发明实施例中粘声波全波形反演结果示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

实际地球介质本身具有粘滞性，导致地震波在传播中的振幅和频率受到严重影响，尤其对于近地表强衰减层等强衰减区域来说，如果不考虑地球介质的吸收衰减作用，会造成全波形反演结果产生较大误差。本发明实施例选取GSLS粘弹性模型，实现基于二维一阶速度-应力粘声波方程的粘声波全波形反演方法。该方法通过采用由多个Maxwell体构成的GSLS模型，可以更好地近似地球介质的常Q特征；通过采用高阶交错网格有限差分法，具有更好地数值频散抑制能力；由于考虑了地球介质的粘弹性特征，能够得到更为准确的纵波速度反演结果。

图1为本发明实施例中粘声波全波形反演方法的示意图。如图1所示，本发明实施例的粘声波全波形反演方法可以包括：

步骤101、根据GSLS模型的粘弹性波动方程，确定时间域一阶速度-应力粘声波方程，其中，所述GSLS模型由多个Maxwell体构成；

步骤102、根据所述时间域一阶速度-应力粘声波方程和二范数目标函数，确定所述时间域一阶速度-应力粘声波方程相应的松弛模量梯度计算公式及纵波速度梯度计算公式；

步骤103、根据所述时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场；

步骤104、根据所述正向传播波场、残差逆时反传波场、松弛模量梯度计算公式及纵波速度梯度计算公式，确定松弛模量梯度和纵波速度梯度；

步骤105、根据所述纵波速度梯度，采用局部最优化方法进行二维粘声波全波形反演。

具体实施时，选取由多个Maxwell体构成的GSLS模型来描述地球介质的粘弹性特征，从基于GSLS模型的粘弹性波动方程出发，得到时间域一阶速度-应力粘声波方程。实施例中，可以根据理想声介质中不存在剪切应力的特点，确定时间域一阶速度-应力粘声波方程如下：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_i}{\∂t}+\frac{\∂p}{\∂x_i}=f_p\\ \frac{\∂p}{\∂t}+M(1+{\mathrm{\τ}}^p)\frac{\∂v_k}{\∂x_k}+{\mathrm{\Σ}}_{l=1}^L{r}_l=0\\ \frac{\∂r_l}{\∂t}+\frac{1}{{\mathrm{\τ}}_{\mathrm{\σ}l}}({\mathrm{M\τ}}^p\frac{\∂v_k}{\∂x_k}+r_l)=0\end{array}\right.$

其中，p为流体平均压力，v_i为波场速度分量，ρ为介质密度，f_p为震源项，L为GSLS模型中Maxwell体的个数，r_l和τ_σl分别为第l个Maxwell体对应的记忆变量和应力松弛时间，τ^p为纵波松弛时间，M为松弛模量，t为时间，x_i和x_k表示波场分量方向，v_k为x_k方向上的波场速度分量，

其中，v_p0为未考虑吸收衰减效应的纵波速度，表示取实部，i为虚数单位(i的平方等于-1)，ω₀为参考角频率。

再根据时间域一阶速度-应力粘声波方程，基于二范数目标函数，推导出该方程相应的松弛模量梯度及纵波速度梯度计算公式；实施例中，松弛模量梯度计算公式可以如下：

${\mathrm{\δM}}^{\′}=-\underset{source}{\Σ}{\∫}_0^{T}dt\frac{p^{foward}{p}^{back}}{{\left(M\right(1+{\mathrm{\τ}}^p\left)\right)}^2}((1+{\mathrm{\τ}}^p)+{\mathrm{\Σ}}_{l=1}^L\frac{1}{{\mathrm{\τ}}_{\mathrm{\σ}l}}{\mathrm{\τ}}^p)$

其中，p^foward为正向传播波场相应的应力，p^back为残差逆时反传波场相应的应力，T为波场传播时间，表示对所有炮求和，

纵波速度梯度计算公式可以如下：

其中，v_p为纵波速度。

具体实施时，根据所述时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场，可以包括：利用真实纵波速度模型和Q模型，采用高阶交错网格有限差分法数值求解所述时间域一阶速度-应力粘声波方程，得到观测炮记录；利用初始纵波速度模型和Q模型，采用高阶交错网格有限差分法数值求解所述时间域一阶速度-应力粘声波方程，得到正向传播波场各时刻的快照及合成炮记录；利用所述观测炮记录和合成炮记录，得到残差记录；利用初始纵波速度模型和Q模型，采用高阶交错网格有限差分法逆时外推所述时间域一阶速度-应力粘声波方程，得到残差逆时反传波场各时刻的快照。

实施例中，在确定正向传播波场及残差逆时反传波场后，根据正向传播波场、残差逆时反传波场、松弛模量梯度计算公式及纵波速度梯度计算公式，确定松弛模量梯度和纵波速度梯度；再根据纵波速度梯度，采用局部最优化方法进行二维粘声波全波形反演。其中，根据纵波速度梯度采用局部最优化方法进行二维粘声波全波形反演，可以包括：在每次迭代中，采用共轭梯度法更新松弛模量和纵波速度模型，共轭方向采用弗莱彻和李维斯提出的FR(Fletcher-Reeves)方法来计算。

由上述实施例可以得知，本发明实施例中粘声波全波形反演方法，由于采用由多个Maxwell体构成的GSLS模型，可以更好地近似地球介质的常Q特征；由于是基于时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法，比基于二阶位移方程的粘声波反演方法在采用相同差分阶数时具有更好地数值频散抑制能力；由于在时间域实现，因此波场正演和残差逆时反传比较直接快速；由于考虑了地球介质本身的粘滞性，能够得到更为准确的纵波速度反演结果。

基于同一发明构思，本发明实施例中还提供了一种粘声波全波形反演装置，如下面的实施例所述。由于该装置解决问题的原理与粘声波全波形反演方法相似，因此该装置的实施可以参见粘声波全波形反演方法的实施，重复之处不再赘述。

图2为本发明实施例中粘声波全波形反演装置的示意图。如图2所示，本发明实施例中粘声波全波形反演装置可以包括：

粘声波方程确定模块201，用于根据GSLS模型的粘弹性波动方程，确定时间域一阶速度-应力粘声波方程，其中，所述GSLS模型由多个Maxwell体构成；

计算公式确定模块202，用于根据所述时间域一阶速度-应力粘声波方程和二范数目标函数，确定所述时间域一阶速度-应力粘声波方程相应的松弛模量梯度计算公式及纵波速度梯度计算公式；

波场确定模块203，用于根据所述时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场；

梯度确定模块204，用于根据所述正向传播波场、残差逆时反传波场、松弛模量梯度计算公式及纵波速度梯度计算公式，确定松弛模量梯度和纵波速度梯度；

全波形反演模块205，用于根据所述纵波速度梯度，采用局部最优化方法进行二维粘声波全波形反演。

具体实施时，波场确定模块203具体可以用于：

利用真实纵波速度模型和Q模型，采用高阶交错网格有限差分法数值求解所述时间域一阶速度-应力粘声波方程，得到观测炮记录；

利用初始纵波速度模型和Q模型，采用高阶交错网格有限差分法数值求解所述时间域一阶速度-应力粘声波方程，得到正向传播波场各时刻的快照及合成炮记录；

利用所述观测炮记录和合成炮记录，得到残差记录；

利用初始纵波速度模型和Q模型，采用高阶交错网格有限差分法逆时外推所述时间域一阶速度-应力粘声波方程，得到残差逆时反传波场各时刻的快照。

具体实施时，全波形反演模块205具体可以用于：

在每次迭代中，采用共轭梯度法更新松弛模量和纵波速度模型，共轭方向采用FR方法来计算。

具体实施时，所述时间域一阶速度-应力粘声波方程可以如下：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_i}{\∂t}+\frac{\∂p}{\∂x_i}=f_p\\ \frac{\∂p}{\∂t}+M(1+{\mathrm{\τ}}^p)\frac{\∂v_k}{\∂x_k}+{\mathrm{\Σ}}_{l=1}^L{r}_l=0\\ \frac{\∂r_l}{\∂t}+\frac{1}{{\mathrm{\τ}}_{\mathrm{\σ}l}}({\mathrm{M\τ}}^p\frac{\∂v_k}{\∂x_k}+r_l)=0\end{array}\right.$

其中，p为流体平均压力，v_i为波场速度分量，ρ为介质密度，f_p为震源项，L为GSLS模型中Maxwell体的个数，r_l和τ_σl分别为第l个Maxwell体对应的记忆变量和应力松弛时间，τ^p为纵波松弛时间，M为松弛模量，t为时间，x_i和x_k表示波场分量方向，v_k为x_k方向上的波场速度分量，

其中，v_p0表示未考虑吸收衰减效应的纵波速度，表示取实部，i为虚数单位(i的平方等于-1)，ω₀为参考角频率，

所述松弛模量梯度计算公式可以如下：

${\mathrm{\δM}}^{\′}=-\underset{source}{\Σ}{\∫}_0^{T}dt\frac{p^{foward}{p}^{back}}{{\left(M\right(1+{\mathrm{\τ}}^p\left)\right)}^2}((1+{\mathrm{\τ}}^p)+{\mathrm{\Σ}}_{l=1}^L\frac{1}{{\mathrm{\τ}}_{\mathrm{\σ}l}}{\mathrm{\τ}}^p)$

其中，p^foward为正向传播波场相应的应力，p^back为残差逆时反传波场相应的应力，T为波场传播时间，表示对所有炮求和，

所述纵波速度梯度计算公式可以如下：

其中，v_p为纵波速度。

综上所述，本发明实施例中基于时间域一阶速度-应力粘声波方程的粘声波全波形反演装置方法及装置，由于采用由多个Maxwell体构成的GSLS模型，可以更好地近似地球介质的常Q特征；由于采用高阶交错网格有限差分法，比基于二阶位移方程的粘声波反演方法在采用相同差分阶数时具有更好地数值频散抑制能力；其中基于广义标准线性固体模型的粘声波全波形反演，直接推导时间域松弛模量梯度和纵波速度梯度计算公式，由于在时间域实现，因此波场正演和残差逆时反传比较直接快速；由于考虑了地球介质本身的粘滞性，能够得到更为准确的纵波速度反演结果。

图3为本发明实施例中未考虑粘弹性影响的声波全波形反演结果示意图；图4为本发明实施例中粘声波全波形反演结果示意图。图3示出了未考虑粘弹性影响的声波全波形反演的纵波速度，图4示出了应用本发明实施例获得的粘声波全波形反演的纵波速度，如图3、图4所示范，对于粘声波观测记录，由于声波全波形反演方法没有考虑地球介质的吸收衰减效应，反演得到的纵波速度是不准确的，而本发明实施例的粘声波反演方法得到了较为准确的反演结果。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种粘声波全波形反演方法，其特征在于，包括：

根据GSLS模型的粘弹性波动方程，确定时间域一阶速度-应力粘声波方程，其中，所述GSLS模型由多个麦克斯韦Maxwell体构成；

根据所述时间域一阶速度-应力粘声波方程和二范数目标函数，确定所述时间域一阶速度-应力粘声波方程相应的松弛模量梯度计算公式及纵波速度梯度计算公式；

根据所述时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场；

根据所述正向传播波场、残差逆时反传波场、松弛模量梯度计算公式及纵波速度梯度计算公式，确定松弛模量梯度和纵波速度梯度；

根据所述纵波速度梯度，采用局部最优化方法进行二维粘声波全波形反演。

2.如权利要求1所述的方法，其特征在于，根据所述时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场，包括：

利用真实纵波速度模型和Q模型，采用高阶交错网格有限差分法数值求解所述时间域一阶速度-应力粘声波方程，得到观测炮记录；

利用初始纵波速度模型和Q模型，采用高阶交错网格有限差分法数值求解所述时间域一阶速度-应力粘声波方程，得到正向传播波场各时刻的快照及合成炮记录；

利用所述观测炮记录和合成炮记录，得到残差记录；

利用初始纵波速度模型和Q模型，采用高阶交错网格有限差分法逆时外推所述时间域一阶速度-应力粘声波方程，得到残差逆时反传波场各时刻的快照。

3.如权利要求1所述的方法，其特征在于，根据所述纵波速度梯度，采用局部最优化方法进行二维粘声波全波形反演，包括：

在每次迭代中，采用共轭梯度法更新松弛模量和纵波速度模型，共轭方向采用弗莱彻-李维斯FR方法来计算。

4.如权利要求1至3任一项所述的方法，其特征在于，所述时间域一阶速度-应力粘声波方程如下：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_i}{\∂t}+\frac{\∂p}{\∂x_i}=f_p\\ \frac{\∂p}{\∂t}+M(1+{\mathrm{\τ}}^p)\frac{\∂v_k}{\∂x_k}+{\mathrm{\Σ}}_{l=1}^L{r}_l=0\\ \frac{\∂r_l}{\∂t}+\frac{1}{{\mathrm{\τ}}_{\mathrm{\σ}l}}(M{\mathrm{\τ}}^p\frac{\∂v_k}{\∂x_k}+r_l)=0\end{array}\right.$

其中，p为流体平均压力，v_i为波场速度分量，ρ为介质密度，f_p为震源项，L为GSLS模型中Maxwell体的个数，r_l和τ_σl分别为第l个Maxwell体对应的记忆变量和应力松弛时间，τ^p为纵波松弛时间，M为松弛模量，t为时间，x_i和x_k表示波场分量方向，v_k为x_k方向上的波场速度分量，

其中，v_p0为未考虑吸收衰减效应的纵波速度，表示取实部，i为虚数单位，ω₀为参考角频率，

所述松弛模量梯度计算公式如下：

${\mathrm{\δM}}^{\′}=-\underset{source}{\Σ}{\∫}_0^{T}dt\frac{p^{foward}{P}^{back}}{{\left(M\right(1+{\mathrm{\τ}}^p\left)\right)}^2}\left(\right(1+{\mathrm{\τ}}^p)+{\mathrm{\Σ}}_{l=1}^L\frac{1}{{\mathrm{\τ}}_{\mathrm{\σ}l}}{\mathrm{\τ}}^p)$

其中，p^foward为正向传播波场相应的应力，p^back为残差逆时反传波场相应的应力，T为波场传播时间，表示对所有炮进行求和，

所述纵波速度梯度计算公式如下：

其中，v_p为纵波速度。

5.一种粘声波全波形反演装置，其特征在于，包括：

粘声波方程确定模块，用于根据GSLS模型的粘弹性波动方程，确定时间域一阶速度-应力粘声波方程，其中，所述GSLS模型由多个Maxwell体构成；

计算公式确定模块，用于根据所述时间域一阶速度-应力粘声波方程和二范数目标函数，确定所述时间域一阶速度-应力粘声波方程相应的松弛模量梯度计算公式及纵波速度梯度计算公式；

波场确定模块，用于根据所述时间域一阶速度-应力粘声波方程，采用高阶交错网格有限差分法确定正向传播波场及残差逆时反传波场；

梯度确定模块，用于根据所述正向传播波场、残差逆时反传波场、松弛模量梯度计算公式及纵波速度梯度计算公式，确定松弛模量梯度和纵波速度梯度；

全波形反演模块，用于根据所述纵波速度梯度，采用局部最优化方法进行二维粘声波全波形反演。

6.如权利要求5所述的装置，其特征在于，所述波场确定模块具体用于：

利用真实纵波速度模型和Q模型，采用高阶交错网格有限差分法数值求解所述时间域一阶速度-应力粘声波方程，得到观测炮记录；

利用初始纵波速度模型和Q模型，采用高阶交错网格有限差分法数值求解所述时间域一阶速度-应力粘声波方程，得到正向传播波场各时刻的快照及合成炮记录；

利用所述观测炮记录和合成炮记录，得到残差记录；

利用初始纵波速度模型和Q模型，采用高阶交错网格有限差分法逆时外推所述时间域一阶速度-应力粘声波方程，得到残差逆时反传波场各时刻的快照。

7.如权利要求5所述的装置，其特征在于，所述全波形反演模块具体用于：

在每次迭代中，采用共轭梯度法更新松弛模量和纵波速度模型，共轭方向采用FR方法来计算。

8.如权利要求5至7任一项所述的装置，其特征在于，所述时间域一阶速度-应力粘声波方程如下：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_i}{\∂t}+\frac{\∂p}{\∂x_i}=f_p\\ \frac{\∂p}{\∂t}+M(1+{\mathrm{\τ}}^p)\frac{\∂v_k}{\∂x_k}+{\mathrm{\Σ}}_{l=1}^L{r}_l=0\\ \frac{\∂r_l}{\∂t}+\frac{1}{{\mathrm{\τ}}_{\mathrm{\σ}l}}(M{\mathrm{\τ}}^p\frac{\∂v_k}{\∂x_k}+r_l)=0\end{array}\right.$

其中，p为流体平均压力，v_i为波场速度分量，ρ为介质密度，f_p为震源项，L为GSLS模型中Maxwell体的个数，r_l和τ_σl分别为第l个Maxwell体对应的记忆变量和应力松弛时间，τ^p为纵波松弛时间，M为松弛模量，t为时间，x_i和x_k表示波场分量方向，v_k为x_k方向上的波场速度分量，

其中，v_p0表示未考虑吸收衰减效应的纵波速度，表示取实部，i为虚数单位，ω₀为参考角频率，

所述松弛模量梯度计算公式如下：

${\mathrm{\δM}}^{\′}=-\underset{source}{\Σ}{\∫}_0^{T}dt\frac{p^{foward}{P}^{back}}{{\left(M\right(1+{\mathrm{\τ}}^p\left)\right)}^2}\left(\right(1+{\mathrm{\τ}}^p)+{\mathrm{\Σ}}_{l=1}^L\frac{1}{{\mathrm{\τ}}_{\mathrm{\σ}l}}{\mathrm{\τ}}^p)$

其中，p^foward为正向传播波场相应的应力，p^back为残差逆时反传波场相应的应力，T为波场传播时间，表示对所有炮求和，

所述纵波速度梯度计算公式如下：

其中，v_p为纵波速度。