CN103091711B - 基于时间域一阶弹性波动方程的全波形反演方法及装置 - Google Patents

Abstract

本发明公开一种基于时间域一阶速度-应力弹性波动方程的全波形反演方法及装置，其中方法包括：采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式；根据时间域一阶速度-应力弹性波动方程确定正向传播波场，根据伴随方程确定逆时外推伴随波场；根据正向传播波场、逆时外推伴随波场和梯度表达式确定目标泛函关于模型参数的梯度；利用梯度类迭代算法进行多尺度全波形反演。本发明在逆时传播波场残差前不需要对波场残差进行预处理，不需要考虑质点速度和质点位移的转换问题，可提高反演计算效率。

Description

基于时间域一阶弹性波动方程的全波形反演方法及装置

技术领域

本发明涉及石油地震勘探速度建模技术领域，尤其涉及基于时间域一阶速度-应力弹性波动方程的全波形反演方法及装置。

背景技术

相对于旅行时层析成像方法，全波形反演不仅可以利用地震走时，还可以利用振幅、相位、波形等信息，因此能够得到精确的高分辨率速度模型。随着勘探的不断深入，地震成像需要解决的问题日益精细化和复杂化，单纯依靠纵波信息来推断地下介质模型已远远不够，为了较好地推测出地下介质模型，需要提供弹性多波的成像信息，因而弹性波全波形反演显得十分重要。

全波形反演是根据观测的地震波场记录来推测地下介质参数的一种方法。基于最小二乘原理，当目标泛函取极小时所对应的模型就认为是要求解的模型。在利用最小二乘解求取地下介质模型参数时至关重要的一步是计算目标泛函关于模型参数的梯度即Fréchet导数。直接求解模型Fréchet导数往往比较困难，需要巨大的计算量，为此许多学者引入了其他思想来求解Fréchet导数，其中伴随方法是求解Fréchet导数的一个强有力工具。

目前非线性全波形反演研究较为深入的是基于梯度类的迭代算法，准确的梯度将加快反演的收敛速度，提高反演的计算效率。基于弹性波动方程的传统目标泛函关于模型参数的梯度是通过二阶弹性波动方程系统导出的(Tarantola(1986)，Mora(1987)，Liu和Tromp(2006))。Crase等(1990)引入更加稳健的目标泛函，采用一阶速度-应力弹性波动方程及相应的伴随方程(除右端项之外，其表达形式与正向一阶速度-应力方程一致)来计算正向传播波场和逆时外推伴随波场，利用基于二阶位移弹性波动方程导出的目标泛函关于模型参数的梯度表达式，对海洋实际数据资料进行了弹性波全波形反演研究。Shipp和Singh(2002)利用大角度海洋反射数据集对P波速度进行了反演，在反演过程中，同样采用一阶速度-应力方程对正演波场进行模拟，为了提高计算效率，他们对梯度表达式进行了变形，用正应力来代替位移对空间的导数。Baumstein等(2009)指出在采用上述策略求解梯度时，必然存在着质点位移和质点速度的转换问题，如果想要得到准确的梯度，需要在逆时外推波场残差之前对其进行预处理即求解波场残差对应力或位移的导数。这样，就使得计算效率较为低下。

发明内容

本发明实施例提供一种基于时间域一阶速度-应力弹性波动方程的全波形反演方法，用以提高反演的计算效率，该方法包括：

采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式；

根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，根据所述伴随方程确定逆时外推伴随波场；

根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度；

利用梯度类迭代算法进行多尺度全波形反演；

采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程如下：

$\mathrm{\ρ}\frac{{\∂w}_x}{\∂t}-(\mathrm{\λ}+2\mathrm{\μ})\frac{{\∂w}_{\mathrm{xx}}}{\∂x}-\mathrm{\λ}\frac{{\∂w}_{\mathrm{zz}}}{\∂x}-\mathrm{\μ}\frac{{\∂w}_{\mathrm{xz}}}{\∂z}=g_x;$

$\mathrm{\ρ}\frac{{\∂w}_z}{\∂t}-(\mathrm{\λ}+2\mathrm{\μ})\frac{{\∂w}_{\mathrm{zz}}}{\∂z}-\mathrm{\λ}\frac{{\∂w}_{\mathrm{xx}}}{\∂z}-\mathrm{\μ}\frac{{\∂w}_{\mathrm{xz}}}{\∂x}=g_z;$

$\frac{{\∂w}_{\mathrm{xx}}}{\∂t}-\frac{{\∂w}_x}{\∂x}=g_{\mathrm{xx}};$

$\frac{{\∂w}_{\mathrm{zz}}}{\∂t}-\frac{{\∂w}_z}{\∂z}=g_{\mathrm{zz}};$

$\frac{{\∂w}_{\mathrm{xz}}}{\∂z}-\frac{{\∂w}_x}{\∂z}-\frac{{\∂w}_z}{\∂x}=g_{\mathrm{xz}};$

确定相应的目标泛函关于模型参数的梯度表达式如下：

${\▿}_{\mathrm{\ρ}}\mathrm{\χ}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂v_x}{\∂t}{w}_x+\frac{{\∂v}_z}{\∂t}{w}_z)\mathrm{dt};$

${\▿}_{\mathrm{\λ}}\mathrm{\χ}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{{\∂v}_x}{\∂x}+\frac{{\∂v}_z}{\∂z})(w_{\mathrm{xx}}+w_{\mathrm{zz}})\mathrm{dt};$

${\▿}_{\mathrm{\μ}}\mathrm{\χ}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T[2(\frac{{\∂v}_x}{\∂x}{w}_{\mathrm{xx}}+\frac{{\∂v}_z}{\∂z}{w}_{\mathrm{zz}})+w_{\mathrm{xz}}(\frac{{\∂v}_x}{\∂z}+\frac{{\∂v}_z}{\∂x})]\mathrm{dt};$

其中，ρ为介质密度，w＝(w_x,w_z,w_xx,w_zz,w_xz)^T为逆时外推伴随波场矢量，λ和μ为Lamé常数， $g={(g_x,g_z,g_{\mathrm{xx}},g_{\mathrm{zz}},g_{\mathrm{xz}})}^T=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}\frac{[v(m;x,t)-v^{\mathrm{obs}}(x,t)]\mathrm{\δ}(x-x_r)}{{\mathrm{\σ}}^2(x_s,x,t)},$ v为波场矢量，v(m；x_s,x,t)为正演模拟的波场，v^obs(x_s,x,t)为观测记录的波场，σ²(x_s,x,t)是观测数据的协方差，m为模型向量，ns为震源个数，T为观测记录的时间；

采用基于摄动理论的伴随方法，确定所述梯度表达式如下：

${\▿}_{\mathrm{\ρ}}\mathrm{\χ}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T\frac{\∂\stackrel{\→}{u}}{\∂t}\frac{\∂\stackrel{\←}{u}}{\∂t}+\frac{\∂\stackrel{\→}{w}}{\∂t}\frac{\∂\stackrel{\←}{w}}{\∂t}\mathrm{dt};$

${\▿}_{\mathrm{\λ}}\mathrm{\χ}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂x}+\frac{\∂\stackrel{\→}{w}}{\∂z})(\frac{\∂\stackrel{\←}{u}}{\∂x}+\frac{\∂\stackrel{\←}{w}}{\∂z})\mathrm{dt};$

${\▿}_{\mathrm{\μ}}\mathrm{\χ}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂z}+\frac{\∂\stackrel{\→}{w}}{\∂x})(\frac{\∂\stackrel{\←}{u}}{\∂z}+\frac{\∂\stackrel{\←}{w}}{\∂x})+2(\frac{\∂\stackrel{\→}{u}}{\∂x}\frac{\∂\stackrel{\←}{u}}{\∂x}+\frac{\∂\stackrel{\→}{w}}{\∂z}\frac{\∂\stackrel{\←}{w}}{\∂z})\mathrm{dt};$

进一步包括：对所述梯度表达式进行变形如下：

${\▿}_{\mathrm{\ρ}}\mathrm{\χ}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T\frac{\∂\stackrel{\→}{u}}{\∂t}\frac{\∂{\stackrel{\←}{v}}_x}{\∂t}+\frac{\∂\stackrel{\→}{w}}{\∂t}\frac{\∂{\stackrel{\←}{v}}_z}{\∂t}\mathrm{dt};$

${\▿}_{\mathrm{\λ}}\mathrm{\χ}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂x}+\frac{\∂\stackrel{\→}{w}}{\∂z})(\frac{\∂{\stackrel{\←}{v}}_x}{\∂x}+\frac{\∂{\stackrel{\←}{v}}_z}{\∂z})\mathrm{dt};$

${\▿}_{\mathrm{\μ}}\mathrm{\χ}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂z}+\frac{\∂\stackrel{\→}{w}}{\∂x})(\frac{\∂{\stackrel{\←}{v}}_x}{\∂z}+\frac{\∂{\stackrel{\←}{v}}_z}{\∂x})+2(\frac{\∂\stackrel{\→}{u}}{\∂x}\frac{\∂{\stackrel{\←}{v}}_x}{\∂x}+\frac{\∂\stackrel{\→}{w}}{\∂z}\frac{\∂{\stackrel{\←}{v}}_z}{\∂z})\mathrm{dt};$

其中，u,w表示水平和垂向位移，v_x,v_z表示水平和垂向应力速度，→箭头表示正向传播，←箭头表示逆时传播；ρ为介质密度，λ和μ为Lamé常数，ns为震源个数，T为观测记录的时间；

根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度，包括：

根据所述正向传播波场、逆时外推伴随波场和变形后的所述梯度表达式确定所述目标泛函关于模型参数的梯度；

根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，包括：

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录；

采用高斯平滑法获得初始模型，对初始模型进行正演模拟获得正向传播波场以及人工合成地震记录；

根据所述伴随方程确定逆时外推伴随波场，包括：

对所述伴随方程采用高阶交错网格有限差分法和PML吸收边界条件进行求解得到逆时外推伴随波场；

利用梯度类迭代算法进行多尺度全波形反演，包括：

采用预条件共轭梯度法进行多尺度全波形反演，反演过程中迭代步长采用非精确线性搜索进行求取。

一个实施例中，上述方法还包括：

对所述伴随方程进行量纲分析如下：

[w_xx]＝[w_xz]＝[w_zz]＝M/LT；

对所述梯度表达式进行量纲分析如下：

$\left[{\▿}_{\mathrm{\λ}}\mathrm{\χ}\right]=\left[{\▿}_{\mathrm{\μ}}\mathrm{\χ}\right]=M/{\mathrm{LT}}^2=\left[\mathrm{\σ}\right];$

其中，L,T,M分别表示长度、时间和质量单位，v_xz＝(v_x,v_z)^T,σ＝(σ_xx,σ_zz,σ_xz)^T。

一个实施例中，对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录时，采用的数值方法为交错网格有限差分法，时间是二阶精度，空间是四阶精度，边界采用最佳匹配层PML吸收边界条件。

本发明实施例还提供一种基于时间域一阶速度-应力弹性波动方程的全波形反演装置，用以提高反演的计算效率，该装置包括：

第一确定模块，用于采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式；

第二确定模块，用于根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，根据所述伴随方程确定逆时外推伴随波场；

第三确定模块，用于根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度；

反演处理模块，用于利用梯度类迭代算法进行多尺度全波形反演；

所述第一确定模块具体用于：

采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程如下：

$\mathrm{\ρ}\frac{{\∂w}_x}{\∂t}-(\mathrm{\λ}+2\mathrm{\μ})\frac{{\∂w}_{\mathrm{xx}}}{\∂x}-\mathrm{\λ}\frac{{\∂w}_{\mathrm{zz}}}{\∂x}-\mathrm{\μ}\frac{{\∂w}_{\mathrm{xz}}}{\∂z}=g_x;$

$\mathrm{\ρ}\frac{{\∂w}_z}{\∂t}-(\mathrm{\λ}+2\mathrm{\μ})\frac{{\∂w}_{\mathrm{zz}}}{\∂z}-\mathrm{\λ}\frac{{\∂w}_{\mathrm{xx}}}{\∂z}-\mathrm{\μ}\frac{{\∂w}_{\mathrm{xz}}}{\∂x}=g_z;$

$\frac{{\∂w}_{\mathrm{xx}}}{\∂t}-\frac{{\∂w}_x}{\∂x}=g_{\mathrm{xx}};$

$\frac{{\∂w}_{\mathrm{zz}}}{\∂t}-\frac{{\∂w}_z}{\∂z}=g_{\mathrm{zz}};$

$\frac{{\∂w}_{\mathrm{xz}}}{\∂z}-\frac{{\∂w}_x}{\∂z}-\frac{{\∂w}_z}{\∂x}=g_{\mathrm{xz}};$

确定相应的目标泛函关于模型参数的梯度表达式如下：

${\▿}_{\mathrm{\ρ}}\mathrm{\χ}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂v_x}{\∂t}{w}_x+\frac{{\∂v}_z}{\∂t}{w}_z)\mathrm{dt};$

${\▿}_{\mathrm{\λ}}\mathrm{\χ}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{{\∂v}_x}{\∂x}+\frac{{\∂v}_z}{\∂z})(w_{\mathrm{xx}}+w_{\mathrm{zz}})\mathrm{dt};$

${\▿}_{\mathrm{\μ}}\mathrm{\χ}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T[2(\frac{{\∂v}_x}{\∂x}{w}_{\mathrm{xx}}+\frac{{\∂v}_z}{\∂z}{w}_{\mathrm{zz}})+w_{\mathrm{xz}}(\frac{{\∂v}_x}{\∂z}+\frac{{\∂v}_z}{\∂x})]\mathrm{dt};$

其中，ρ为介质密度，w＝(w_x,w_z,w_xx,w_zz,w_xz)^T为逆时外推伴随波场矢量，λ和μ为Lamé常数， $g={(g_x,g_z,g_{\mathrm{xx}},g_{\mathrm{zz}},g_{\mathrm{xz}})}^T=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}\frac{[v(m;x,t)-v^{\mathrm{obs}}(x,t)]\mathrm{\δ}(x-x_r)}{{\mathrm{\σ}}^2(x_s,x,t)},$ v为波场矢量，v(m；x_s,x,t)为正演模拟的波场，v^obs(x_s,x,t)为观测记录的波场，σ²(x_s,x,t)是观测数据的协方差，m为模型向量，ns为震源个数，T为观测记录的时间；

所述第一确定模块具体用于：

采用基于摄动理论的伴随方法，确定所述梯度表达式如下：

${\▿}_{\mathrm{\ρ}}\mathrm{\χ}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T\frac{\∂\stackrel{\→}{u}}{\∂t}\frac{\∂\stackrel{\←}{u}}{\∂t}+\frac{\∂\stackrel{\→}{w}}{\∂t}\frac{\∂\stackrel{\←}{w}}{\∂t}\mathrm{dt};$

${\▿}_{\mathrm{\λ}}\mathrm{\χ}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂x}+\frac{\∂\stackrel{\→}{w}}{\∂z})(\frac{\∂\stackrel{\←}{u}}{\∂x}+\frac{\∂\stackrel{\←}{w}}{\∂z})\mathrm{dt};$

${\▿}_{\mathrm{\μ}}\mathrm{\χ}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂z}+\frac{\∂\stackrel{\→}{w}}{\∂x})(\frac{\∂\stackrel{\←}{u}}{\∂z}+\frac{\∂\stackrel{\←}{w}}{\∂x})+2(\frac{\∂\stackrel{\→}{u}}{\∂x}\frac{\∂\stackrel{\←}{u}}{\∂x}+\frac{\∂\stackrel{\→}{w}}{\∂z}\frac{\∂\stackrel{\←}{w}}{\∂z})\mathrm{dt};$

在采用基于摄动理论的伴随方法，确定所述梯度表达式后，进一步对所述梯度表达式进行变形如下：

${\▿}_{\mathrm{\λ}}\mathrm{\χ}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂x}+\frac{\∂{\stackrel{\←}{v}}_x}{\∂z})(\frac{\∂{\stackrel{\→}{v}}_x}{\∂x}+\frac{\∂{\stackrel{\←}{v}}_z}{\∂z})\mathrm{dt};$

${\▿}_{\mathrm{\μ}}\mathrm{\χ}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂x}+\frac{\∂\stackrel{\→}{w}}{\∂z})(\frac{\∂{\stackrel{\←}{v}}_x}{\∂x}+\frac{\∂{\stackrel{\←}{v}}_z}{\∂z})\mathrm{dt};$

${\▿}_{\mathrm{\μ}}\mathrm{\χ}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂\stackrel{\→}{u}}{\∂z}+\frac{\∂\stackrel{\→}{w}}{\∂x})(\frac{\∂{\stackrel{\←}{v}}_x}{\∂z}+\frac{\∂{\stackrel{\←}{v}}_z}{\∂x})+2(\frac{\∂\stackrel{\→}{u}}{\∂x}\frac{\∂{\stackrel{\←}{v}}_x}{\∂x}+\frac{\∂\stackrel{\→}{w}}{\∂z}\frac{\∂{\stackrel{\←}{v}}_z}{\∂z})\mathrm{dt};$

其中，u,w表示水平和垂向位移，v_x,v_z表示水平和垂向应力速度，→箭头表示正向传播，←箭头表示逆时传播；ρ为介质密度，λ和μ为Lamé常数，ns为震源个数，T为观测记录的时间；

所述第二确定模块具体用于：

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录；

采用高斯平滑法获得初始模型，对初始模型进行正演模拟获得正向传播波场以及人工合成地震记录；

所述第二确定模块具体用于：

对所述伴随方程采用高阶交错网格有限差分法和PML吸收边界条件进行求解得到逆时外推伴随波场；

所述第三确定模块具体用于：

根据所述正向传播波场、逆时外推伴随波场和变形后的所述梯度表达式确定所述目标泛函关于模型参数的梯度；

所述反演处理模块具体用于：

采用预条件共轭梯度法进行多尺度全波形反演，反演过程中迭代步长采用非精确线性搜索进行求取。

一个实施例中，上述装置还包括：

量纲分析模块，用于：

对所述伴随方程进行量纲分析如下：

[w_xx]＝[w_xz]＝[w_zz]＝M/LT；

对所述梯度表达式进行量纲分析如下：

$\left[{\▿}_{\mathrm{\λ}}\mathrm{\χ}\right]=\left[{\▿}_{\mathrm{\μ}}\mathrm{\χ}\right]=M/{\mathrm{LT}}^2=\left[\mathrm{\σ}\right];$

其中，L,T,M分别表示长度、时间和质量单位，v_xz＝(v_x,v_z)^T,σ＝(σ_xx,σ_zz,σ_xz)^T。

一个实施例中，所述第二确定模块具体用于：

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录时，采用的数值方法为交错网格有限差分法，时间是二阶精度，空间是四阶精度，边界采用PML吸收边界条件。

本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演方法及装置，直接推导时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式，采用其进行多尺度全波形反演时，在逆时传播波场残差之前不需要对波场残差进行预处理即可得到准确的梯度，在求取目标泛函关于模型参数梯度的过程中不需要考虑质点速度和质点位移的转换问题，而且减少了伴随波场对空间导数的计算，提高了反演的计算效率。

附图说明

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

图1为本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演方法的处理流程图；

图2为本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演装置的结构示意图。

具体实施方式

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

图1为本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演方法的处理流程图。如图1所示，本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演方法可以包括：

步骤101、采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式；

步骤102、根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，根据所述伴随方程确定逆时外推伴随波场；

步骤103、根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度；

步骤104、利用梯度类迭代算法进行多尺度全波形反演。

由图1所示流程可以得知，本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演方法，直接推导时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式，采用其进行多尺度全波形反演时避免了质点位移和质点速度之间的转换问题而且在求取梯度过程中不需要对波场残差进行预处理，提高了计算效率。

具体实施时，采用基于摄动理论的伴随方法，直接推导时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式。具体的，采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程可以如下：

$\mathrm{\ρ}\frac{{\∂w}_x}{\∂t}-(\mathrm{\λ}+2\mathrm{\μ})\frac{{\∂w}_{\mathrm{xx}}}{\∂x}-\mathrm{\λ}\frac{{\∂w}_{\mathrm{zz}}}{\∂x}-\mathrm{\μ}\frac{{\∂w}_{\mathrm{xz}}}{\∂z}=g_x;$

$\mathrm{\ρ}\frac{{\∂w}_z}{\∂t}-(\mathrm{\λ}+2\mathrm{\μ})\frac{{\∂w}_{\mathrm{zz}}}{\∂z}-\mathrm{\λ}\frac{{\∂w}_{\mathrm{xx}}}{\∂z}-\mathrm{\μ}\frac{{\∂w}_{\mathrm{xz}}}{\∂x}=g_z;$

$\frac{{\∂w}_{\mathrm{xx}}}{\∂t}-\frac{{\∂w}_x}{\∂x}=g_{\mathrm{xx}};$

$\frac{{\∂w}_{\mathrm{zz}}}{\∂t}-\frac{{\∂w}_z}{\∂z}=g_{\mathrm{zz}};$

$\frac{{\∂w}_{\mathrm{xz}}}{\∂z}-\frac{{\∂w}_x}{\∂z}-\frac{{\∂w}_z}{\∂x}=g_{\mathrm{xz}};$

确定相应的目标泛函关于模型参数的梯度表达式可以如下：

${\▿}_{\mathrm{\ρ}}\mathrm{\χ}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{\∂v_x}{\∂t}{w}_x+\frac{{\∂v}_z}{\∂t}{w}_z)\mathrm{dt};$

${\▿}_{\mathrm{\λ}}\mathrm{\χ}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T(\frac{{\∂v}_x}{\∂x}+\frac{{\∂v}_z}{\∂z})(w_{\mathrm{xx}}+w_{\mathrm{zz}})\mathrm{dt};$

${\▿}_{\mathrm{\μ}}\mathrm{\χ}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T[2(\frac{{\∂v}_x}{\∂x}{w}_{\mathrm{xx}}+\frac{{\∂v}_z}{\∂z}{w}_{\mathrm{zz}})+w_{\mathrm{xz}}(\frac{{\∂v}_x}{\∂z}+\frac{{\∂v}_z}{\∂x})]\mathrm{dt};$

其中，ρ为介质密度，w＝(w_x,w_z,w_xx,w_zz,w_xz)^T为逆时外推伴随波场矢量，λ和μ为Lamé(拉梅)常数， $g={(g_x,g_z,g_{\mathrm{xx}},g_{\mathrm{zz}},g_{\mathrm{xz}})}^T=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}\frac{[v(m;x,t)-v^{\mathrm{obs}}(x,t)]\mathrm{\δ}(x-x_r)}{{\mathrm{\σ}}^2(x_s,x,t)},$ v为波场矢量，v(m；x_s,x,t)为正演模拟的波场，v^obs(x_s,x,t)为观测记录的波场，σ²(x_s,x,t)是观测数据的协方差，m为模型向量，ns为震源个数，T为观测记录的时间。

下面详细说明具体的推导过程：

时间域二维各向同性均匀介质中的一阶速度-应力弹性波动方程写成算子形式可以表示为：

Lv＝f                 (1)

其中， $L=A{\∂}_x+B{\∂}_z+C{\∂}_t,v={(v_x,v_z,{\mathrm{\σ}}_{\mathrm{xx}},{\mathrm{\σ}}_{\mathrm{zz}},{\mathrm{\σ}}_{\mathrm{xz}})}^T,f={(f_x,f_z,f_{\mathrm{xx}},f_{\mathrm{zz}},f_{\mathrm{xz}})}^T---\left(2\right)$

$A=\left(\begin{array}{ccccc}0& 0& -1& 0& 0\\ 0& 0& 0& 0& -1\\ -(\mathrm{\λ}+2\mathrm{\μ})& 0& 0& 0& 0\\ -\mathrm{\λ}& 0& 0& 0& 0\\ 0& -\mathrm{\μ}& 0& 0& 0\end{array}\right),B=\left(\begin{array}{ccccc}0& 0& 0& 0& -1\\ 0& 0& 0& -1& 0\\ 0& -\mathrm{\λ}& 0& 0& 0\\ 0& -(\mathrm{\λ}+2\mathrm{\μ})& 0& 0& 0\\ -\mathrm{\μ}& 0& 0& 0& 0\end{array}\right),C=\left(\begin{array}{ccccc}\mathrm{\ρ}& & & & \\ & \mathrm{\ρ}& & & \\ & & 1& & \\ & & & 1& \\ & & & & 1\end{array}\right)---\left(3\right)$

ρ为介质密度，λ和μ为Lamé常数，L为正演差分算子，v为波场矢量，f为震源矢量。(1)式满足的初始条件和边界条件分别为：

v(x,0)＝0,σ_xz|_z＝0＝σ_zz|_z＝0＝0           (4)

全波形反演是利用观测的地震波场记录来求取地下介质的模型参数的一种方法。基于最小二乘原理，目标泛函取极小值所对应的模型参数就是要求解的模型。目标泛函可表示为：

$\mathrm{\χ}\left(m\right)=\frac{1}{2}\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T{\∫}_{\mathrm{\Ω}}\frac{{\left|\right|v(m;x_s,x,t)-v^{\mathrm{obs}}(x_s,x,t)\left|\right|}^2}{{\mathrm{\σ}}^2(x_s,x,t)}\mathrm{\δ}(x-x_r)\mathrm{d\Ωdt}---\left(5\right)$

其中，v^obs(x_s,x,t)为观测记录的波场，v(m；x_s,x,t)为正演模拟的波场，σ²(x_s,x,t)是观测数据的协方差，m为模型向量，ns为震源个数，T为观测记录的时间，Ω为地下介质模型的区域。根据摄动理论，模型的扰动会产生相应的扰动波场。扰动模型和扰动波场满足如下关系：

$\mathrm{L\δv}=\hat{f}---\left(6\right)$

其中， $\mathrm{\δv}={({\mathrm{\δv}}_x,{\mathrm{\δv}}_z,{\mathrm{\δ\σ}}_{\mathrm{xx}},{\mathrm{\δ\σ}}_{\mathrm{zz}},{\mathrm{\δ\σ}}_{\mathrm{xz}})}^T,\hat{f}=(-{\mathrm{\δC}\∂}_t-{\mathrm{\δA}\∂}_x-{\mathrm{\δB}\∂}_z)\·v---\left(7\right)$

Tarantola(1986)定义伴随算子：

<L^*w,δv>＝<w,Lδv>              (8)

其中，L^*为伴随算子，w＝(w_x,w_z,w_xx,w_zz,w_xz)^T为伴随波场矢量，<,>代表对时间和空间的积分。把(6)式带入(8)式并利用分部积分法可得：

令伴随算子满足如下伴随方程：

$L^{*}w=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}\frac{[v(m;x,t)-v^{\mathrm{obs}}(x,t)]\mathrm{\&delta;}(x-x_r)}{{\mathrm{\&sigma;}}^2(x_s,x,t)}---\left(10\right)$

对(10)式展开：

$\begin{array}{c}\mathrm{\&rho;}\frac{{\&PartialD;w}_x}{\&PartialD;t}-(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;x}-\mathrm{\&lambda;}\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;x}-\mathrm{\&mu;}\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;z}=g_x;\\ \mathrm{\&rho;}\frac{{\&PartialD;w}_z}{\&PartialD;t}-(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;z}-\mathrm{\&lambda;}\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;z}-\mathrm{\&mu;}\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;x}=g_z;\\ \frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;t}-\frac{{\&PartialD;w}_x}{\&PartialD;x}=g_{\mathrm{xx}};\\ \frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;t}-\frac{{\&PartialD;w}_z}{\&PartialD;z}=g_{\mathrm{zz}};\end{array}---\left(11\right)$

$\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;t}-\frac{{\&PartialD;w}_x}{\&PartialD;z}-\frac{{\&PartialD;w}_z}{\&PartialD;x}=g_{\mathrm{xz}};$

其中 $g={(g_x,g_z,g_{\mathrm{xx}},g_{\mathrm{zz}},g_{\mathrm{xz}})}^T=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}\frac{[v(m;x,t)-v^{\mathrm{obs}}(x,t)]\mathrm{\&delta;}(x-x_r)}{{\mathrm{\&sigma;}}^2(x_s,x,t)}---\left(12\right)$

由(11)式可以看出，基于时间域一阶速度-应力方程推导出来的伴随方程和正向的一阶速度-应力方程不再一致，这区别于传统的伴随方程。对(11)式进一步分析，定义：

${\&PartialD;}_t\stackrel{\&LeftArrow;}{u}=w_x,{\&PartialD;}_t\stackrel{\&LeftArrow;}{w}=w_z---\left(13\right)$

其中，←箭头表示逆时传播，把(13)式带入(11)式，并略去右端项可得：

$\begin{array}{c}{\mathrm{\&rho;}\&PartialD;}_{\mathrm{tt}}^2\stackrel{\&LeftArrow;}{u}-(\mathrm{\&lambda;}+2\mathrm{\&mu;}){\&PartialD;}_{\mathrm{xx}}^2\stackrel{\&LeftArrow;}{u}-(\mathrm{\&lambda;}+\mathrm{\&mu;}){\&PartialD;}_{\mathrm{xz}}^2\stackrel{\&LeftArrow;}{w}-{\mathrm{\&mu;}\&PartialD;}_{\mathrm{zz}}^2\stackrel{\&LeftArrow;}{u}=0\\ {\mathrm{\&rho;}\&PartialD;}_{\mathrm{tt}}^2\stackrel{\&LeftArrow;}{w}-(\mathrm{\&lambda;}+2\mathrm{\&mu;}){\&PartialD;}_{\mathrm{zz}}^2\stackrel{\&LeftArrow;}{w}-(\mathrm{\&lambda;}+\mathrm{\&mu;}){\&PartialD;}_{\mathrm{xz}}^2\stackrel{\&LeftArrow;}{u}-{\mathrm{\&mu;}\&PartialD;}_{\mathrm{xx}}^2\stackrel{\&LeftArrow;}{w}=0\end{array}---\left(14\right)$

(14)式与去除震源项的二阶位移表示的弹性波动方程是一致的，不同的是它满足终止条件和新的边界条件((9)式)。由此可得，(11)式是二阶位移表示的弹性波动方程的一种新的分裂形式。

目标泛函(5)式的变分为：

$\mathrm{\&delta;\&chi;}\left(m\right)=\frac{1}{2}\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T{\&Integral;}_{\mathrm{\&Omega;}}\frac{{\left|\right|v(m;x_s,x,t)-v^{\mathrm{obs}}(x_s,x,t)\left|\right|}^2}{{\mathrm{\&sigma;}}^2(x_s,x,t)}\mathrm{\&delta;}(x-x_r)\mathrm{\&delta;vd\&Omega;dt}---\left(15\right)$

把(6)式和(10)式带入(8)式并展开与(15)式联立可得目标泛函关于密度和Lamé常数的梯度表达式：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;v_x}{\&PartialD;t}{w}_x+\frac{{\&PartialD;v}_z}{\&PartialD;t}{w}_z)\mathrm{dt}---\left(16\right)$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{{\&PartialD;v}_x}{\&PartialD;x}+\frac{{\&PartialD;v}_z}{\&PartialD;z})(w_{\mathrm{xx}}+w_{\mathrm{zz}})\mathrm{dt}---\left(17\right)$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T[2(\frac{{\&PartialD;v}_x}{\&PartialD;x}{w}_{\mathrm{xx}}+\frac{{\&PartialD;v}_z}{\&PartialD;z}{w}_{\mathrm{zz}})+w_{\mathrm{xz}}(\frac{{\&PartialD;v}_x}{\&PartialD;z}+\frac{{\&PartialD;v}_z}{\&PartialD;x})]\mathrm{dt}---\left(18\right)$

当目标泛函是应力的函数时即：

$\mathrm{\&chi;}\left(m\right)=\frac{1}{2}\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T{\&Integral;}_{\mathrm{\&Omega;}}{\left|\right|\mathrm{\&sigma;}(m;x_s,x,t)-{\mathrm{\&sigma;}}^{\mathrm{obs}}(x_s,x,t)\left|\right|}^2\mathrm{\&delta;}(x-x_r)\mathrm{d\&Omega;dt}---\left(19\right)$

那么目标泛函关于弹性张量的梯度的单位应与应力的单位一致。

具体实施时，还可以对所述伴随方程进行量纲分析如下：

[w_xx]＝[w_xz]＝[w_zz]＝M/LT；

对所述梯度表达式进行量纲分析如下：

$\left[{\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}\right]=\left[{\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}\right]=M/{\mathrm{LT}}^2=\left[\mathrm{\&sigma;}\right];$

其中，L,T,M分别表示长度、时间和质量单位，v_xz＝(v_x,v_z)^T,σ＝(σ_xx,σ_zz,σ_xz)^T。

具体的，对(1)式进行量纲分析可得：

[v_xz]＝L/T,[σ]＝M/LT²             (20)

其中，v_xz＝(v_x,v_z)^T,σ＝(σ_xx,σ_zz,σ_xz)^T，L,T,M分别表示长度、时间和质量单位。在逆时外推之前，不对应力波场残差做任何预处理，直接把应力波场残差加载到伴随方程的右端，并对(10)式进行量纲分析有：

[w_xx]＝[w_xz]＝[w_zz]＝M/LT             (21)

对新的目标泛函梯度(16-18)式进行量纲分析(假设对时间的积分对梯度的单位没有影响)：

$\left[{\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}\right]=\left[{\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}\right]=M/{\mathrm{LT}}^2=\left[\mathrm{\&sigma;}\right]---\left(22\right)$

由此可以验证，本发明实施例不需要对波场残差进行预处理就可以得到准确的梯度。

具体实施时，若采用传统方案，Baumstein等(2009)指出若想要得到准确的梯度，需要在逆时外推波场残差之前对其进行预处理即求解波场残差对应力或位移的导数。利用一阶速度-应力弹性波动方程和相应的伴随方程(除右端项外，和一阶速度应力方程一样)来求解二阶位移弹性波动方程推导出的目标泛函关于模型梯度的表达式为传统策略。其中梯度表达形式为：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt}.$

当观测记录为质点速度时，本发明实施例给出了一种不需要对波场残差进行预处理就可以得到准确梯度的技巧，此时目标泛函关于模型的梯度表达式可以相应地变形为：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;t}\mathrm{dt}---\left(23\right)$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt}---\left(24\right)$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;z}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt}---\left(25\right)$

其中，u,w表示水平和垂向位移，v_x,v_z表示水平和垂向应力速度，→箭头表示正向传播，←箭头表示逆时传播；ρ为介质密度，λ和μ为Lamé常数，ns为震源个数，T为观测记录的时间。

即，实施中，可以先用传统方法采用基于摄动理论的伴随方法确定梯度表达式如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

在确定梯度表达式后，进一步包括：

对梯度表达式进行变形如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;z}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt};$

则后续根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度，可以包括：根据所述正向传播波场、逆时外推伴随波场和变形后的所述梯度表达式确定所述目标泛函关于模型参数的梯度。

具体实施时，根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，可以包括：

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录；

采用高斯平滑法获得初始模型，对初始模型进行正演模拟获得正向传播波场以及人工合成地震记录。

具体实施时，对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录时，采用的数值方法可以为交错网格有限差分法，时间是二阶精度，空间是四阶精度，边界采用PML(Perfect Matched Layer，最佳匹配层)吸收边界条件。

具体实施时，根据所述伴随方程确定逆时外推伴随波场，可以包括：

对所述伴随方程采用高阶交错网格有限差分法和PML吸收边界条件进行求解得到逆时外推伴随波场。

上述采用基于多重网格法的多尺度全波形反演方法是考虑到：采用复杂模型记录的地震观测数据进行反演时，目标函数往往存在着大量的局部极小值。为了减小局部极小值，降低全波形反演非线性的程度，实施中可以采用基于多重网格法的多尺度全波形反演方法。

具体实施时，利用梯度类迭代算法进行多尺度全波形反演，可以包括：

采用预条件共轭梯度法进行多尺度全波形反演，反演过程中迭代步长采用非精确线性搜索进行求取。

如此实施是考虑到：当采用共轭梯度法进行全波形反演时，一个不可避免的问题是求解出的目标泛函梯度在震源位置附近会产生较大的异常值。另外震源的振幅会随着传播距离的增大逐渐减小，因而位于深部的散射体较浅部的散射体对梯度的贡献相对较小。为了加速收敛，消除震源位置的异常幅值，平衡浅部散射体和深部散射体对梯度的贡献，在较佳的实施例中可以采用预条件共轭梯度法。

基于同一发明构思，本发明实施例中还提供了一种基于时间域一阶速度-应力弹性波动方程的全波形反演装置，如下面的实施例所述。由于该装置解决问题的原理与基于时间域一阶速度-应力弹性波动方程的全波形反演方法相似，因此该装置的实施可以参见基于时间域一阶速度-应力弹性波动方程的全波形反演方法的实施，重复之处不再赘述。

图2为本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演装置的结构示意图。如图2所示，本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演装置可以包括：

第一确定模块201，用于采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式；

第二确定模块202，用于根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，根据所述伴随方程确定逆时外推伴随波场；

第三确定模块203，用于根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度；

反演处理模块204，用于利用梯度类迭代算法进行多尺度全波形反演。

具体实施时，所述第一确定模块201具体可以用于：

采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程如下：

$\mathrm{\&rho;}\frac{{\&PartialD;w}_x}{\&PartialD;t}-(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;x}-\mathrm{\&lambda;}\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;x}-\mathrm{\&mu;}\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;z}=g_x;$

$\mathrm{\&rho;}\frac{{\&PartialD;w}_z}{\&PartialD;t}-(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;z}-\mathrm{\&lambda;}\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;z}-\mathrm{\&mu;}\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;x}=g_z;$

$\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;t}-\frac{{\&PartialD;w}_x}{\&PartialD;x}=g_{\mathrm{xx}};$

$\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;t}-\frac{{\&PartialD;w}_z}{\&PartialD;z}=g_{\mathrm{zz}};$

$\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;t}-\frac{{\&PartialD;w}_x}{\&PartialD;z}-\frac{{\&PartialD;w}_z}{\&PartialD;x}=g_{\mathrm{xz}};$

确定相应的目标泛函关于模型参数的梯度表达式如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;v_x}{\&PartialD;t}{w}_x+\frac{{\&PartialD;v}_z}{\&PartialD;t}{w}_z)\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{{\&PartialD;v}_x}{\&PartialD;x}+\frac{{\&PartialD;v}_z}{\&PartialD;z})(w_{\mathrm{xx}}+w_{\mathrm{zz}})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T[2(\frac{{\&PartialD;v}_x}{\&PartialD;x}{w}_{\mathrm{xx}}+\frac{{\&PartialD;v}_z}{\&PartialD;z}{w}_{\mathrm{zz}})+w_{\mathrm{xz}}(\frac{{\&PartialD;v}_x}{\&PartialD;z}+\frac{{\&PartialD;v}_z}{\&PartialD;x})]\mathrm{dt};$

其中，ρ为介质密度，w＝(w_x,w_z,w_xx,w_zz,w_xz)^T为逆时外推伴随波场矢量，λ和μ为Lamé常数， $g={(g_x,g_z,g_{\mathrm{xx}},g_{\mathrm{zz}},g_{\mathrm{xz}})}^T=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}\frac{[v(m;x,t)-v^{\mathrm{obs}}(x,t)]\mathrm{\&delta;}(x-x_r)}{{\mathrm{\&sigma;}}^2(x_s,x,t)},$ v为波场矢量，v(m；x_s,x,t)为正演模拟的波场，v^obs(x_s,x,t)为观测记录的波场，σ²(x_s,x,t)是观测数据的协方差，m为模型向量，ns为震源个数，T为观测记录的时间。

具体实施时，上述装置还可以包括：

量纲分析模块，用于：

对所述伴随方程进行量纲分析如下：

[w_xx]＝[w_xz]＝[w_zz]＝M/LT；

对所述梯度表达式进行量纲分析如下：

$\left[{\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}\right]=\left[{\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}\right]=M/{\mathrm{LT}}^2=\left[\mathrm{\&sigma;}\right];$

其中，L,T,M分别表示长度、时间和质量单位，v_xz＝(v_x,v_z)^T,σ＝(σ_xx,σ_zz,σ_xz)^T。

具体实施时，所述第一确定模块201具体可以用于：

采用基于摄动理论的伴随方法，确定所述梯度表达式如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

在采用基于摄动理论的伴随方法，确定所述梯度表达式后，进一步对所述梯度表达式进行变形如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;z}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt};$

其中，u,w表示水平和垂向位移，v_x,v_z表示水平和垂向应力速度，→箭头表示正向传播，←箭头表示逆时传播；ρ为介质密度，λ和μ为Lamé常数，ns为震源个数，T为观测记录的时间；

所述第三确定模块203具体可以用于：根据所述正向传播波场、逆时外推伴随波场和变形后的所述梯度表达式确定所述目标泛函关于模型参数的梯度。

具体实施时，所述第二确定模块202具体可以用于：

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录；

采用高斯平滑法获得初始模型，对初始模型进行正演模拟获得正向传播波场以及人工合成地震记录。

具体实施时，所述第二确定模块202具体可以用于：

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录时，采用的数值方法为交错网格有限差分法，时间是二阶精度，空间是四阶精度，边界采用PML吸收边界条件。

具体实施时，所述第二确定模块202具体可以用于：

对所述伴随方程采用高阶交错网格有限差分法和PML吸收边界条件进行求解得到逆时外推伴随波场。

具体实施时，所述反演处理模块204具体可以用于：

采用预条件共轭梯度法进行多尺度全波形反演，反演过程中迭代步长采用非精确线性搜索进行求取。

综上所述，本发明实施例中基于时间域一阶速度-应力弹性波动方程的全波形反演方法及装置，直接推导时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式，采用其进行多尺度全波形反演时，在逆时传播波场残差之前不需要对波场残差进行预处理即可得到准确的梯度，在求取目标泛函关于模型参数梯度的过程中不需要考虑质点速度和质点位移的转换问题，而且减少了伴随波场对空间导数的计算，提高了反演的计算效率。

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

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

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

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

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

Claims (6)

1.一种基于时间域一阶速度-应力弹性波动方程的全波形反演方法，其特征在于，该方法包括：

采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式；

根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，根据所述伴随方程确定逆时外推伴随波场；

根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度；

利用梯度类迭代算法进行多尺度全波形反演；

采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程如下：

$\mathrm{\&rho;}\frac{{\&PartialD;w}_x}{\&PartialD;t}-(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;x}-\mathrm{\&lambda;}\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;x}-\mathrm{\&mu;}\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;z}=g_x;$

$\mathrm{\&rho;}\frac{{\&PartialD;w}_z}{\&PartialD;t}-(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;z}-\mathrm{\&lambda;}\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;z}-\mathrm{\&mu;}\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;x}=g_z;$

$\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;t}-\frac{{\&PartialD;w}_x}{\&PartialD;x}=g_{\mathrm{xx}};$

$\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;t}-\frac{{\&PartialD;w}_z}{\&PartialD;z}=g_{\mathrm{zz}};$

$\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;t}-\frac{{\&PartialD;w}_x}{\&PartialD;z}-\frac{{\&PartialD;w}_z}{\&PartialD;x}=g_{\mathrm{xz};}$

确定相应的目标泛函关于模型参数的梯度表达式如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{{\&PartialD;v}_x}{\&PartialD;t}{w}_x+\frac{{\&PartialD;v}_z}{\&PartialD;t}{w}_z)\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{{\&PartialD;v}_x}{\&PartialD;x}+\frac{{\&PartialD;v}_z}{\&PartialD;z})(w_{\mathrm{xx}}+w_{\mathrm{zz}})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T[2(\frac{{\&PartialD;v}_x}{\&PartialD;x}{w}_{\mathrm{xx}}+\frac{{\&PartialD;v}_z}{\&PartialD;z}{w}_{\mathrm{zz}})+w_{\mathrm{xz}}(\frac{{\&PartialD;v}_x}{\&PartialD;z}+\frac{{\&PartialD;v}_z}{\&PartialD;x})]\mathrm{dt};$

其中，ρ为介质密度，w＝(w_x,w_z,w_xx,w_zz,w_xz)^T为逆时外推伴随波场矢量，λ和μ为Lamé常数， $g={(g_x,g_z,g_{\mathrm{xx}},g_{\mathrm{zz}},g_{\mathrm{xz}})}^T=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}\frac{[v(m;x,t)-v^{\mathrm{obs}}(x,t)]\mathrm{\&delta;}(x-x_r)}{{\mathrm{\&sigma;}}^2(x_s,x,t)},$ v为波场矢量，v(m；x_s,x,t)为正演模拟的波场，v^obs(x_s,x,t)为观测记录的波场，σ²(x_s,x,t)是观测数据的协方差，m为模型向量，ns为震源个数，T为观测记录的时间；

采用基于摄动理论的伴随方法，确定所述梯度表达式如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

进一步包括：对所述梯度表达式进行变形如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=\underset{\mathrm{is}=10}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;z}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\text{\&PartialD;x}}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt};$

其中，u,w表示水平和垂向位移，v_x,v_z表示水平和垂向应力速度，→箭头表示正向传播，←箭头表示逆时传播；ρ为介质密度，λ和μ为Lamé常数，ns为震源个数，T为观测记录的时间；

根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，包括：

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录；

采用高斯平滑法获得初始模型，对初始模型进行正演模拟获得正向传播波场以及人工合成地震记录；

根据所述伴随方程确定逆时外推伴随波场，包括：

对所述伴随方程采用高阶交错网格有限差分法和PML吸收边界条件进行求解得到逆时外推伴随波场；

根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度，包括：

根据所述正向传播波场、逆时外推伴随波场和变形后的所述梯度表达式确定所述目标泛函关于模型参数的梯度；

利用梯度类迭代算法进行多尺度全波形反演，包括：

采用预条件共轭梯度法进行多尺度全波形反演，反演过程中迭代步长采用非精确线性搜索进行求取。

2.如权利要求1所述的方法，其特征在于，还包括：

对所述伴随方程进行量纲分析如下：

[w_xx]＝[w_xz]＝[w_zz]＝M/LT；

对所述梯度表达式进行量纲分析如下：

$\left[{\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}\right]=\left[{\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}\right]=M/{\mathrm{LT}}^2=\left[\mathrm{\&sigma;}\right];$

其中，L,T,M分别表示长度、时间和质量单位，v_xz＝(v_x,v_z)^T,σ＝(σ_xx,σ_zz,σ_xz)^T。

3.如权利要求1所述的方法，其特征在于，对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录时，采用的数值方法为交错网格有限差分法，时间是二阶精度，空间是四阶精度，边界采用最佳匹配层PML吸收边界条件。

4.一种基于时间域一阶速度-应力弹性波动方程的全波形反演装置，其特征在于，该装置包括：

第一确定模块，用于采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程和相应的目标泛函关于模型参数的梯度表达式；

第二确定模块，用于根据所述时间域一阶速度-应力弹性波动方程确定正向传播波场，根据所述伴随方程确定逆时外推伴随波场；

第三确定模块，用于根据所述正向传播波场、逆时外推伴随波场和所述梯度表达式确定所述目标泛函关于模型参数的梯度；

反演处理模块，用于利用梯度类迭代算法进行多尺度全波形反演；

所述第一确定模块具体用于：

采用基于摄动理论的伴随方法，确定时间域一阶速度-应力弹性波动方程的伴随方程如下：

$\mathrm{\&rho;}\frac{{\&PartialD;w}_x}{\&PartialD;t}-(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;x}-\mathrm{\&lambda;}\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;x}-\mathrm{\&mu;}\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;z}=g_x;$

$\mathrm{\&rho;}\frac{{\&PartialD;w}_z}{\&PartialD;t}-(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;z}-\mathrm{\&lambda;}\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;z}-\mathrm{\&mu;}\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;x}=g_z;$

$\frac{{\&PartialD;w}_{\mathrm{xx}}}{\&PartialD;t}-\frac{{\&PartialD;w}_x}{\&PartialD;x}=g_{\mathrm{xx}};$

$\frac{{\&PartialD;w}_{\mathrm{zz}}}{\&PartialD;t}-\frac{{\&PartialD;w}_z}{\&PartialD;z}=g_{\mathrm{zz}};$

$\frac{{\&PartialD;w}_{\mathrm{xz}}}{\&PartialD;t}-\frac{{\&PartialD;w}_x}{\&PartialD;z}-\frac{{\&PartialD;w}_z}{\&PartialD;x}=g_{\mathrm{xz};}$

确定相应的目标泛函关于模型参数的梯度表达式如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{{\&PartialD;v}_x}{\&PartialD;t}{w}_x+\frac{{\&PartialD;v}_z}{\&PartialD;t}{w}_z)\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{{\&PartialD;v}_x}{\&PartialD;x}+\frac{{\&PartialD;v}_z}{\&PartialD;z})(w_{\mathrm{xx}}+w_{\mathrm{zz}})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=+\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T[2(\frac{{\&PartialD;v}_x}{\&PartialD;x}{w}_{\mathrm{xx}}+\frac{{\&PartialD;v}_z}{\&PartialD;z}{w}_{\mathrm{zz}})+w_{\mathrm{xz}}(\frac{{\&PartialD;v}_x}{\&PartialD;z}+\frac{{\&PartialD;v}_z}{\&PartialD;x})]\mathrm{dt};$

其中，ρ为介质密度，w＝(w_x,w_z,w_xx,w_zz,w_xz)^T为逆时外推伴随波场矢量，λ和μ为Lamé常数， $g={(g_x,g_z,g_{\mathrm{xx}},g_{\mathrm{zz}},g_{\mathrm{xz}})}^T=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}\frac{[v(m;x,t)-v^{\mathrm{obs}}(x,t)]\mathrm{\&delta;}(x-x_r)}{{\mathrm{\&sigma;}}^2(x_s,x,t)},$ v为波场矢量，v(m；x_s,x,t)为正演模拟的波场，v^obs(x_s,x,t)为观测记录的波场，σ²(x_s,x,t)是观测数据的协方差，m为模型向量，ns为震源个数，T为观测记录的时间；

所述第一确定模块具体用于：

采用基于摄动理论的伴随方法，确定所述梯度表达式如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}\frac{\&PartialD;\stackrel{\&LeftArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;\stackrel{\&LeftArrow;}{w}}{\&PartialD;z})\mathrm{dt};$

在采用基于摄动理论的伴随方法，确定所述梯度表达式后，进一步对所述梯度表达式进行变形如下：

${\&dtri;}_{\mathrm{\&rho;}}\mathrm{\&chi;}=-\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;t}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;t}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;t}\mathrm{dt};$

${\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt};$

${\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}=\underset{\mathrm{is}=10}{\overset{\mathrm{ns}}{\mathrm{\&Sigma;}}}{\&Integral;}_0^T(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\&PartialD;z}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;x})(\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;z}+\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;x})+2(\frac{\&PartialD;\stackrel{\&RightArrow;}{u}}{\text{\&PartialD;x}}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_x}{\&PartialD;x}+\frac{\&PartialD;\stackrel{\&RightArrow;}{w}}{\&PartialD;z}\frac{\&PartialD;{\stackrel{\&LeftArrow;}{v}}_z}{\&PartialD;z})\mathrm{dt};$

其中，u,w表示水平和垂向位移，v_x,v_z表示水平和垂向应力速度，→箭头表示正向传播，←箭头表示逆时传播；ρ为介质密度，λ和μ为Lamé常数，ns为震源个数，T为观测记录的时间；

所述第二确定模块具体用于：

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录；

采用高斯平滑法获得初始模型，对初始模型进行正演模拟获得正向传播波场以及人工合成地震记录；

所述第二确定模块具体用于：

对所述伴随方程采用高阶交错网格有限差分法和PML吸收边界条件进行求解得到逆时外推伴随波场；

所述第三确定模块具体用于：

根据所述正向传播波场、逆时外推伴随波场和变形后的所述梯度表达式确定所述目标泛函关于模型参数的梯度；

所述反演处理模块具体用于：

采用预条件共轭梯度法进行多尺度全波形反演，反演过程中迭代步长采用非精确线性搜索进行求取。

5.如权利要求4所述的装置，其特征在于，还包括：

量纲分析模块，用于：

对所述伴随方程进行量纲分析如下：

[w_xx]＝[w_xz]＝[w_zz]＝M/LT；

对所述梯度表达式进行量纲分析如下：

$\left[{\&dtri;}_{\mathrm{\&lambda;}}\mathrm{\&chi;}\right]=\left[{\&dtri;}_{\mathrm{\&mu;}}\mathrm{\&chi;}\right]=M/{\mathrm{LT}}^2=\left[\mathrm{\&sigma;}\right];$

其中，L,T,M分别表示长度、时间和质量单位，v_xz＝(v_x,v_z)^T,σ＝(σ_xx,σ_zz,σ_xz)^T。

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

对所述时间域一阶速度-应力弹性波动方程进行离散求解，得到观测地震记录时，采用的数值方法为交错网格有限差分法，时间是二阶精度，空间是四阶精度，边界采用PML吸收边界条件。