基于激发主能量优化算法的全波形反演梯度算子的提取方法

摘要

本发明公开了一种基于激发主能量优化算法的全波形反演梯度算子的提取方法，通过建立全波形反演参数模型；采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；将震源置于地表炮点处，利用优化有限差分算法计算该震源产生的正传波场；确定震源波场的正向传播过程中地下各点震源波场最大值及其对应的时刻，并保存地下各点在激发主能量时窗内的波场；再利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量时窗之内，通过基于激发主能量优化算法得到全波形反演梯度算子。本发明在保证精度的前提下，大大减少了I/O操作，提升了计算效率。

说明书

技术领域

本发明涉及石油地震勘探速度建模技术领域，尤其涉及一种基于 激发主能量优化算法的全波形反演梯度算子的提取方法。

背景技术

地震波的传播速度是地震勘探中的一个基本参数，在地震资料 处理和地震解释的诸多环节中扮演着重要的角色。影响地震波传播 速度的因素有很多，诸如地下介质的岩性、密度、孔隙度、埋藏深 度等等，这些因素最终体现在地震波的走时、波形、振幅和相位等 属性上变化。只有充分了解速度及影响因素的相互关系，才有可能 更为准确地估计出地下的结构和构造。我国大多数勘探区域越来越 复杂，呈现出地表构造和地下构造双复杂特征。传统的速度分析方 法已经很难达到所需的精度要求，也会直接影响成像质量。常用的 速度分析方法有：叠加速度分析、偏移速度分析、层析速度反演和 全波形反演。这四种速度分析方法的准确性依次升高，复杂度也依 次升高。

全波形反演技术是目前精度很高的速度分析和反演工具，其利用 叠前地震记录与正演模拟炮记录之间的寻求最佳匹配来获取最优反 演速度场。在反演过程中不仅考虑了地震波传播的旅行时等运动学信 息，而且加入了振幅、相位、波形等动力学信息，能够适应强横向变 速介质和各向异性介质的速度反演。

全波形反演的发展大致可以分为三个阶段，第一个阶段是从20 世纪80年代到90年代，这一时期是全波形反演理论形成和初步验证 阶段，Lailly、Tarantola等人建立了基于波动方程的波形反演的理 论构架，之后一批学者分别利用简单模型对反射和透射波形反演效果 进行了理论验证。Gauthier等人率先对理论模型合成的多偏移距数 据进行了声波和弹性波全波形反演测试，反演结果既展示了全波形反 演方法的应用潜力，也暴露出其根本缺陷——对初始模型的严重依赖 性。

第二个阶段是从上世纪90年代初到21世纪初，为全波形稳定发 展阶段，全波形反演方法形成了多个分支和变种，其中以Pratt等提 出的频率域波形反演方法最具代表性，Pratt提出频率域声波和弹性 波波形反演理论，并利用实际井间地震数据对方法进行验证。该方法 利用稀疏矩阵直接解法求解频率域的离散化声波或弹性波波动方程， 只要频率与物性参数不改变，通过一次稀疏矩阵LU分解和多次高斯 消元法中的前向、后向回代运算，即可完成所有炮点的单频波场响应 计算；应用研究方面以井间地震波形反演和海洋天然气水合物的叠前 波形反演最具影响力。此外为解决全波形反演因初始模型不准而陷入 局部极小问题，Zhou等人先后给出了波动方程走时和波形联合反演 方法，以充分利用走时和波形两类数据对不同尺度模型参数变化的敏 感度的优势，使井间数据全波形反演更为稳定。

自21世纪初至今，全波形反演达到了蓬勃发展的阶段，全波形 反演逐步从二维反演走向三维反演，从声波走向弹性波反演，从单参 数反演到多参数反演，从各向同性到各向异性反演。此外还有 Laplace/Laplace-Fourier域全波形反演、各向异性全波形反演、先 验信息约束的全波形反演和基于超级炮和相位编码的全波形反演。此 外随着计算机水平的进步，全波形反演的计算效率进一步得到了提升， 其实现平台逐步从单机单核走向了多核多节点并行平台、GPU并行平 台。在一定程度上，全波形反演的计算效率完全取决于大数据的I/O 交换。

发明内容

为了提高时间域全波形反演的计算效率，本发明提供了一种基于 激发主能量优化算法的全波形反演梯度算子的提取方法，其利用有限 差分算法计算震源波场和波场误差反传波场，再根据激发主能量优化 算法计算模型更新量的梯度。该方法可以在保障计算精度的同时，减 少大量I/O操作，大大提高全波形反演的计算效率，相比于常规梯度 算法，在保持精度的前提下，提高了近3倍的效率。

本发明解决其技术问题所采用的技术方案是：

一种基于激发主能量优化算法的全波形反演梯度算子的提取方 法，包括：

建立全波形反演参数模型；

利用常规的叠前单炮数据作为观测数据，建立正反演观测系统；

将震源置于地表炮点处，利用优化有限差分算法计算该震源产生 的正传波场；

确定震源波场的正向传播过程中地下各点震源波场最大值及其 对应的时刻，并保存地下各点在激发主能量时窗内的波场；

利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地 下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量 时窗之内，通过基于激发主能量优化算法得到全波形反演梯度算子。

前述的基于激发主能量优化算法的全波形反演梯度算子的提取 方法，具体包括如下步骤：

(1)建立全波形反演参数模型；

(2)采用常规的叠加前单炮数据作为观测数据，建立正反演观 测系统；

(3)将震源置于地表炮点处，利用下式表示的优化有限差分算 法计算该震源产生的正传波场：

$\left(\begin{array}{c}{U}_{i,k}^{n+1}={2U}_{i,k}^n-U_{i,k}^{n-1}+\frac{\Delta t^2{v}^2}{{\mathrm{\Delta x}}^2}[c_0{U}_{i,k}^n+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{i+m,k}^n+U_{i-m,k}^n)]\\ +\frac{\Delta t^2{V}^2}{{\mathrm{\Delta z}}^2}[c_0{U}_{i,k}^n+\underset{m=1}{\overset{N}{\Sigma }}{c}_m(U_{i,k+m}^n+U_{i,k-m}^n)]\end{array}\right)$

其中U为震源产生的正传波场，c_m为优化有限差分系数，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x 和z空间坐标，n为时间坐标；

(4)确定震源波场的正向传播过程中地下各点最大值及其对 应的时刻Ts；

(5)定义激发主能量时窗长度为正传波场中一个子波的时间长 度w，保存地下各点在激发主能量时窗内的波场，即以地下各点的最大值对应时刻Ts为时窗中心的整个主能量时窗之内的波场值；

(6)利用步骤(3)中所述优化有限差分正演模拟算法得到叠前单 炮合成数据RU，并其与观测地震数据d_obs求差得到波场残差RU-d_obs；

(7)利用下式表示的优化有限差分模拟算法计算波场残差逆时 传播波场：

$\left(\begin{array}{c}{L}^{*}\left(R^{*}(\mathrm{RU}-d_{\mathrm{obs}})\right)=U_{\mathrm{bi},k}^{n-1}={2U}_{\mathrm{bi},k}^n-U_{\mathrm{bi},k}^{n+1}+\frac{\Delta t^2{v}^2}{{\mathrm{\Delta x}}^2}[c_0{U}_{\mathrm{bi},k}^n+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{\mathrm{bi}+m,k}^n+U_{\mathrm{bi}-m,k}^n)]\\ +\frac{\Delta t^2{V}^2}{{\mathrm{\Delta z}}^2}[c_0{U}_{\mathrm{bi},k}^n+\underset{m=1}{\overset{N}{\Sigma }}{c}_m(U_{\mathrm{bi},k+m}^n+U_{\mathrm{bi},k-m}^n)]\end{array}\right)$

其中RU-d_obs为波场残差，L^*为反传算子，c_m为优化有限差分系数，U_b为波场残差逆时传播波场，v为介质速度，Δx和Δz为空间采样间隔， Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标；

(8)判断地下任意一点的不同时刻波场残差的逆时传播波场是 否位于该点激发主能量时窗之内，如果位于，则按下式提取基于激发 主能量优化算法的用于全波形反演的梯度算子：

$g\left(v\right)=\frac{\mathrm{\delta E}\left(v\right)}{\mathrm{\delta v}}=\frac{2}{v^3}\underset{s}{\Sigma }\underset{t=\mathrm{Ts}-\frac{w}{2}}{\overset{t=\mathrm{Ts}+\frac{w}{2}}{\Sigma }}\frac{{\partial }^{2}U}{{\partial t}^2}{L}^{*}\left(R^{*}(\mathrm{RU}-d_{\mathrm{obs}})\right)$

其中RU-d_obs为波场残差，L^*为反传算子，为正演波场对时间的二 阶偏导，v为空间速度分布，Ts为激发最大振幅对应时间步，w为激 发主能量时窗长度，g(v)为全波形反演的梯度。

上述方案更进一步包括：

步骤(1)中全波形反演的模型参数求取采用迭代方法求解，迭 代公式为

m^k+1＝m^k+Δm(1)

其中m^k+1，m^k分别表示第k次迭代模型参数和第k+1次迭代模型参数， Δm是参数模型的更新量，其表达式为：

Δm＝-αg(2)

其中g是参数模型更新量的梯度算子；α是计算步长。所以利用先验 信息建立全波形反演参数的初始模型；

步骤(3)中将震源置于地表炮点处，采用优化有限差分算法计 算震源激发的正传波场；

正传波场U满足如下的波动方程：

$\frac{{\partial }^{2}U}{\partial x^2}+\frac{{\partial }^{2}U}{{\partial z}^2}=\frac{1}{v^2(x,z)}\frac{{\partial }^{2}U}{{\partial t}^2}---\left(3\right)$

其中U(x,z,t)为震源激发的正传波场；v(x,z)为介质速度；

利用Taylor展开进行时间和空间差分，将时间离散形式：

$U(t+\mathrm{\Delta t})=U\left(t\right)-U(t-\mathrm{\Delta t})+{\left(\mathrm{v\Delta t}\right)}^2(\frac{{\partial }^{2}U}{{\partial x}^2}+\frac{{\partial }^{2}U}{{\partial z}^2})+O\left({\mathrm{\Delta t}}^2\right)---\left(4\right)$

和空间差分的形式

$\left(\begin{array}{c}\frac{{\partial }^{2}U}{{\partial x}^2}=\frac{1}{{\mathrm{\Delta x}}^2}[c_0{U}_{i,k}+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{i+m,k}+U_{i-m,k})]\\ \frac{{\partial }^{2}U}{{\partial z}^2}=\frac{1}{\Delta z^2}[c_0{U}_{i,k}+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{i,k+m}+U_{i,k-m})]\end{array}\right)---\left(5\right)$

带入波动方程，得到各向同性介质正传波场数值计算公式

$\left(\begin{array}{c}{U}_{i,k}^{n+1}={2U}_{i,k}^n-U_{i,k}^{n-1}+\frac{{\mathrm{\Delta t}}^2{v}^2}{{\mathrm{\Delta x}}^2}[c_0{U}_{i,k}^n+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{i+m,k}^n+U_{i-m,k}^n)]\\ +\frac{\Delta t^2{V}^2}{{\mathrm{\Delta z}}^2}[c_0{U}_{i,k}^n+\underset{m=1}{\overset{N}{\Sigma }}{c}_m(U_{i,k+m}^n+U_{i,k-m}^n)]\end{array}\right)---\left(6\right)$

其中c_m为优化有限差分系数，v为介质速度，Δx和Δz为 空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间 坐标；

步骤(4)是在震源波场的正向传播过程中，根据正传波场U计算 地下各点的最大值及其对应的时刻Ts；

步骤(7)中利用优化有限差分模拟算法计算波场残差逆时传播 波场，其中逆时传播波场满足如下方程：

$\frac{{\partial }^2{U}_b}{{\partial x}^2}+\frac{{\partial }^2{U}_b}{{\partial z}^2}=\frac{1}{v^2(x,z)}\frac{{\partial }^2{U}_b}{{\partial t}^2}---\left(7\right)$

其中U_b(x,z,t)为波场误差的逆时传播波场；v(x,z)为介质速度，

利用Taylor展开进行时间和空间差分，即将时间离散形式：

$U_b(t-\mathrm{\Delta t})=U_b\left(t\right)-U_b(t+\mathrm{\Delta t})+{\left(\mathrm{v\Delta t}\right)}^2(\frac{{\partial }^2{U}_b}{{\partial x}^2}+\frac{{\partial }^2{U}_b}{{\partial z}^2})+O\left({\mathrm{\Delta t}}^2\right)---\left(8\right)$

和空间差分的形式

$\left(\begin{array}{c}\frac{{\partial }^2{U}_b}{{\partial x}^2}=\frac{1}{{\mathrm{\Delta x}}^2}[c_0{U}_{\mathrm{bi},k}+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{\mathrm{bi}+m,k}+U_{\mathrm{bi}-m,k})]\\ \frac{{\partial }^2{U}_b}{{\partial z}^2}=\frac{1}{{\mathrm{\Delta z}}^2}[c_0{U}_{\mathrm{bi},k}+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{\mathrm{bi},k+m}+U_{\mathrm{bi},k-m})]\end{array}\right)---\left(9\right)$

带入逆时传播波场满足的方程，则得出波场误差的逆时传播波场数值 计算公式：

$\left(\begin{array}{c}{L}^{*}\left(R^{*}(\mathrm{RU}-d_{\mathrm{obs}})\right)=U_{\mathrm{bi},k}^{n-1}={2U}_{\mathrm{bi},k}^n-U_{\mathrm{bbi},k}^{n+1}+\frac{\Delta t^2{v}^2}{{\mathrm{\Delta x}}^2}[c_0{U}_{\mathrm{bi},k}^n+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{\mathrm{bi}+m,k}^n+U_{\mathrm{bi}-m,k}^n)]\\ +\frac{{\mathrm{\Delta t}}^2{V}^2}{{\mathrm{\Delta z}}^2}[c_0{U}_{\mathrm{bi},k}^n+\underset{m=1}{\overset{N}{\Sigma }}{c}_m(U_{\mathrm{bi},k+m}^n+U_{\mathrm{bi},k-m}^n)]\end{array}\right)---\left(10\right)$

其中RU-d_obs为波场残差，L^*为反传算子，c_m为优化有限差分系数，U_b为波场残差逆时传播波场，v为介质速度，Δx和Δz为空间采样间隔， Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标。

本发明的有益效果是：传统的全波形算法在构造梯度时需要保存 地震震源波场全部时刻波场信息，基于激发主能量优化梯度算法定义 了构成梯度的主能量，只需保存构成激发主能量时窗的正传波场，并 与相应反传波场计算即可得到模型更新量的梯度，大大减少了I/O操 作，在保持精度的前提下，可以提高3倍的效率。

附图说明

图1是基于激发主能量优化算法的全波形反演梯度算子提取方 法的工作流程图。

图2是真实参数模型图。

图3是初始参数模型图。

图4是模型各点激发最大振幅图。

图5是模型各点激发最大振幅所对应时刻。

图6是模型数据常规梯度算法计算的梯度。

图7是模型数据激发主能量优化算法计算的梯度。

图8是实际资料常规梯度算法计算的梯度。

图9是实际资料激发主能量优化算法计算的梯度。

图10是实际资料常规梯度算法反演结果。

图11是实际资料激发主能量优化算法反演结果。

图12是实际资料单道反演结果对比图。

图13是实际资料反演效率对比图。

具体实施方式

为使本发明的目的、特征和优点能更明显易懂，下文特举出较佳 实施例，并配合所附图式，作详细说明如下：

实施例1

结合附图1，一种基于激发主能量优化算法的全波形反演梯度算 子的提取方法，包括如下步骤：

建立全波形反演参数模型；

采用常规的叠前单炮数据作为观测数据，建立正反演观测系统；

将震源置于地表炮点处，利用优化有限差分算法计算该震源产生 的正传波场；

确定震源波场的正向传播过程中地下各点震源波场最大值及其 对应的时刻，并保存地下各点在激发主能量时窗内的波场；

利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地 下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量 时窗之内，通过基于激发主能量优化算法得到全波形反演梯度算子。

具体步骤为：

步骤1：全波形反演的模型参数求取采用迭代方法求解，迭代公 式为

m^k+1＝m^k+Δm(1)

其中m^k+1，m^k分别表示第k次迭代模型参数和第k+1次迭代模型参数， Δm是参数模型的更新量，其表达式为：

Δm＝-αg(2)

其中g是参数模型更新量的梯度算子；α是计算步长。所以利用先验 信息建立全波形反演参数的初始模型。

步骤2：因为本发明是针对叠前单炮数据，所以采用常规的叠前 单炮数据作为观测数据，建立正反演观测系统。

步骤3：将震源置于地表炮点处，计算正传波场U。正传波场U满 足如下的波动方程：

$\frac{{\partial }^{2}U}{\partial x^2}+\frac{{\partial }^{2}U}{{\partial z}^2}=\frac{1}{v^2(x,z)}\frac{{\partial }^{2}U}{{\partial t}^2}---\left(3\right)$

其中U(x,z,t)为震源激发的正传波场；v(x,z)为介质速度。

本发明采用优化有限差分算法计算震源激发的正传波场。利用 Taylor展开进行时间和空间差分，将时间离散形式：

$U(t+\mathrm{\Delta t})=U\left(t\right)-U(t-\mathrm{\Delta t})+{\left(\mathrm{v\Delta t}\right)}^2(\frac{{\partial }^{2}U}{{\partial x}^2}+\frac{{\partial }^{2}U}{{\partial z}^2})+O\left({\mathrm{\Delta t}}^2\right)---\left(4\right)$

和空间差分的形式

$\left(\begin{array}{c}\frac{{\partial }^{2}U}{{\partial x}^2}=\frac{1}{{\mathrm{\Delta x}}^2}[c_0{U}_{i,k}+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{i+m,k}+U_{i-m,k})]\\ \frac{{\partial }^{2}U}{{\partial z}^2}=\frac{1}{\Delta z^2}[c_0{U}_{i,k}+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{i,k+m}+U_{i,k-m})]\end{array}\right)---\left(5\right)$

带入波动方程，得到各向同性介质正传波场数值计算公式

$\left(\begin{array}{c}{U}_{i,k}^{n+1}={2U}_{i,k}^n-U_{i,k}^{n-1}+\frac{{\mathrm{\Delta t}}^2{v}^2}{{\mathrm{\Delta x}}^2}[c_0{U}_{i,k}^n+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{i+m,k}^n+U_{i-m,k}^n)]\\ +\frac{\Delta t^2{V}^2}{{\mathrm{\Delta z}}^2}[c_0{U}_{i,k}^n+\underset{m=1}{\overset{N}{\Sigma }}{c}_m(U_{i,k+m}^n+U_{i,k-m}^n)]\end{array}\right)---\left(6\right)$

其中c_m为优化有限差分系数，v为介质速度，Δx和Δz为 空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间 坐标；

步骤4：在震源波场的正向传播过程中，根据正传波场U计算地 下各点的最大值及其对应的时刻Ts。

步骤5：定义激发主能量时窗长度为正传波场中一个子波的时间 长度w，保存以地下各点的最大值对应时刻Ts为时窗中心的整个 主能量时窗之内的波场值，即保存地下各点在激发主能量时窗内的波 场，为计算梯度奠定基础。

步骤6：利用优化有限差分正演模拟算法得到叠前单炮合成数据 RU，并将其与观测地震数据d_obs求差得到波场残差RU-d_obs，

步骤7：利用优化有限差分模拟算法计算波场残差逆时传播波场。 逆时传播波场满足如下方程：

$\frac{{\partial }^2{U}_b}{{\partial x}^2}+\frac{{\partial }^2{U}_b}{{\partial z}^2}=\frac{1}{v^2(x,z)}\frac{{\partial }^2{U}_b}{{\partial t}^2}---\left(7\right)$

其中U_b(x,z,t)为波场误差的逆时传播波场；v(x,z)为介质速度。

本发明采用优化有限差分算法计算波场误差的逆时传播波场。利 用Taylor展开进行时间和空间差分，即将时间离散形式：

$U_b(t-\mathrm{\Delta t})=U_b\left(t\right)-U_b(t+\mathrm{\Delta t})+{\left(\mathrm{v\Delta t}\right)}^2(\frac{{\partial }^2{U}_b}{{\partial x}^2}+\frac{{\partial }^2{U}_b}{{\partial z}^2})+O\left({\mathrm{\Delta t}}^2\right)---\left(8\right)$

和空间差分的形式

$\left(\begin{array}{c}\frac{{\partial }^2{U}_b}{{\partial x}^2}=\frac{1}{{\mathrm{\Delta x}}^2}[c_0{U}_{\mathrm{bi},k}+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{\mathrm{bi}+m,k}+U_{\mathrm{bi}-m,k})]\\ \frac{{\partial }^2{U}_b}{{\partial z}^2}=\frac{1}{{\mathrm{\Delta z}}^2}[c_0{U}_{\mathrm{bi},k}+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{\mathrm{bi},k+m}+U_{\mathrm{bi},k-m})]\end{array}\right)---\left(9\right)$

带入逆时传播波场满足的方程，则可以得出波场误差的逆时传播波场 数值计算公式：

$\left(\begin{array}{c}{L}^{*}\left(R^{*}(\mathrm{RU}-d_{\mathrm{obs}})\right)=U_{\mathrm{bi},k}^{n-1}={2U}_{\mathrm{bi},k}^n-U_{\mathrm{bbi},k}^{n+1}+\frac{\Delta t^2{v}^2}{{\mathrm{\Delta x}}^2}[c_0{U}_{\mathrm{bi},k}^n+\underset{m=1}{\overset{M}{\Sigma }}{c}_m(U_{\mathrm{bi}+m,k}^n+U_{\mathrm{bi}-m,k}^n)]\\ +\frac{{\mathrm{\Delta t}}^2{V}^2}{{\mathrm{\Delta z}}^2}[c_0{U}_{\mathrm{bi},k}^n+\underset{m=1}{\overset{N}{\Sigma }}{c}_m(U_{\mathrm{bi},k+m}^n+U_{\mathrm{bi},k-m}^n)]\end{array}\right)---\left(10\right)$

其中RU-d_obs为波场残差，L^*为反传算子，c_m为优化有限差分系数，U_b为波场残差逆时传播波场，v为介质速度，Δx和Δz为空间采样间隔， Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标。

步骤8：因为常规的时间域内基于时间相关条件的梯度算子表达 式为：

$g\left(v\right)=\frac{\mathrm{\delta E}\left(v\right)}{\mathrm{\delta v}}=\frac{2}{v^3}\underset{s}{\Sigma }\underset{t}{\Sigma }\frac{{\partial }^{2}U}{{\partial t}^2}{L}^{*}\left(R^{*}(\mathrm{RU}-d_{\mathrm{obs}})\right)---\left(11\right)$

其中RU-d_obs为波场残差，L^*为反传算子，为正演波场对时间的二 阶偏导，v为空间速度分布。

而本发明的激发主能量优化算法是在激发主能量时窗内利用正 传波场和波场误差逆时波场来计算梯度，其表达式为：

$g\left(v\right)=\frac{\mathrm{\delta E}\left(v\right)}{\mathrm{\delta v}}=\frac{2}{v^3}\underset{s}{\Sigma }\underset{t=\mathrm{Ts}-\frac{w}{2}}{\overset{t=\mathrm{Ts}+\frac{w}{2}}{\Sigma }}\frac{{\partial }^{2}U}{{\partial t}^2}{L}^{*}\left(R^{*}(\mathrm{RU}-d_{\mathrm{obs}})\right)---\left(12\right)$

其中RU-d_obs为波场残差，L^*为反传算子，为正演波场对时间的二 阶偏导，v为空间速度分布，Ts为激发最大振幅对应时间步，w为激 发主能量时窗长度。

所以在步骤7之后，判断地下各点不同时刻的波场残差逆时传 播波场是否位于该点激发主能量时窗之内，如果位于，则按公式(12) 提取基于激发主能量优化算法的用于全波形反演的梯度算子。

激发主能量优化梯度算法在波场的正向传播过程中只需保存构 成激发主能量时窗的正传波场，并与相应反传波场计算即可得到模型 更新量的梯度，大大减少了I/O操作，在保持精度的前提下，可以提 高3倍的效率。

应用实验例1

本发明基于激发主能量优化算法计算的梯度算子的提取方法， 应用于模型数据，取得了理想的计算效果。图2为真实参数模型； 图3为初始参数模型；首先将震源置于地表炮点处，利用优化有限 差分算法计算该震源产生的正传波场；确定震源波场的正向传播过 程中地下各点震源波场最大值(图4)及其对应的时刻(图5)，并保存 地下各点在激发主能量时窗内的波场；再利用优化有限差分模拟算 法计算波场残差逆时传播波场，判断地下各点不同时刻的波场残差 逆时传播波场是否位于该点激发主能量时窗之内，通过基于激发主 能量优化算法得到全波形反演梯度算子(图7)。与常规算法计算梯度 结果(图6)相比，基于激发主能量优化算法计算的梯度结果精度与其 相当，但是计算速度更快，计算效率大大提高。

应用实验例2

应用本发明基于激发主能量优化算法计算的梯度算子的提取方 法，应用于油田某块区的地震数据，取得了精确的反演效果和明显的 效率提高。图8为常规梯度算法计算的梯度，图9为基于激发主能量 优化算法计算的梯度，二者计算精度相当。图10为常规梯度算法的 反演结果，图11为基于激发主能量优化算法的反演结果，图12为单 道反演结果对比图，反演结果表面主能量梯度算法和常规梯度算法误 差在5％以下，即在全波形反演过程中，基于激发主能量优化算法的 梯度提取方法满足精度的需求。图13为反演效率对比图。常规梯度 算法进行一次迭代的时间为320min，激发主能量优化梯度算法迭代 一次的时间为90min。总结，在多炮的全波形反演算法中，在相同的 参数情况下，激发主能量优化梯度算法相比于常规梯度算法，在保持 精度的前提下，提高了近3倍的效率。