CN105652321A - 一种粘声各向异性最小二乘逆时偏移成像方法 - Google Patents

Abstract

本发明公开了一种粘声各向异性介质最小二乘逆时偏移成像方法，属于石油勘探领域，本发明通过一种粘声各向异性拟微分方程实现对地震波在粘声各向异性介质中传播的准确模拟，通过添加稳定的规则化算子实现地震波稳定的逆时传播，在最小二乘反演的框架下，构建了新的粘声各向异性偏移算子和反偏移算子，通过求取梯度方向，更新原始成像剖面，达到提高成像精度的目的。本发明能够对同时存在粘滞性和各向异性的复杂地下构造进行高精度成像，既能克服传统成像方法在同时处理同时存在粘滞性和各向异性介质的不足，又通过最小二乘方式消除了成像噪音，得到真振幅的成像剖面。

Description

一种粘声各向异性最小二乘逆时偏移成像方法

技术领域

本发明属于石油勘探领域，具体涉及一种粘声各向异性最小二乘逆时偏移成像方法。

背景技术

地球介质广泛存在着粘滞性和各向异性，而且这两种性质常常是同时存在的，介质的粘滞性影响地震波传播的能量，而介质的各向异性影响地震波传播的形态。因此，在进行地震波偏移成像方法研究时，需要同时考虑地下介质的粘滞性和各向异性，消除两者对地震波传播的影响，传统基于理想的弹性介质及其声学假设都是不成立的，需要对同时存在粘滞性和各向异性介质的成像方法进行研究。

传统的逆时偏移成像方法存在低频成像噪音、成像能量不均衡、假象严重及分辨率低等缺点，最小二乘逆时偏移方法是一种基于线性化反演理论的保幅成像方法，其基本思想是在宏观背景模型的基础上估计出一个最优化的扰动部分。很多模型及实际资料试验表明，它可以提高偏移剖面的分辨率，压制偏移假象，更重要的是，它可以给出真实地下反射界面的相对振幅。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种粘声各向异性最小二乘逆时偏移成像方法，本发明能够对同时存在粘滞性和各向异性的复杂地下构造进行高精度成像，既能克服传统成像方法在同时处理同时存在粘滞性和各向异性介质的不足，又通过最小二乘方式消除了成像噪音，得到真振幅的成像剖面。

为了实现上述目的，本发明采用如下技术方案：

所述的粘声各向异性最小二乘逆时偏移成像方法，具体包括如下步骤：

步骤1：输入速度模型、各向异性Thomsen参数模型、品质因子Q模型及慢度扰动模型，并建立观测系统；

步骤2：采用有限差分方法进行粘声各向异性介质正演模拟：

正演模拟采用如下的波动方程：

$\left\{\begin{array}{c}\frac{{\∂}^{2}p}{\∂t^2}+\frac{{\mathrm{\τv}}_{px}}{2}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p}{\∂t}-v_{px}^2\frac{{\∂}^{2}p}{\∂x^2}-v_{pn}{v}_{pz}\frac{{\∂}^{2}q}{\∂z^2}=0\\ \frac{{\∂}^{2}q}{\∂t^2}+\frac{{\mathrm{\τv}}_{pz}}{2}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q}{\∂t}-v_{pn}{v}_{pz}\frac{{\∂}^{2}p}{\∂x^2}-v_{pz}^2\frac{{\∂}^{2}q}{\∂z^2}=0\end{array}\right.---\left(1\right)$

将衰减项和在波数域进行差分离散：

$\begin{array}{c}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂U}{\∂t}{|}_{t=n}\≈\frac{1}{\mathrm{\Δ}t}{F}^{-1}\[\left|k_x\right|F(U^n-U^{n-1})\]\\ \sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂U}{\∂t}{|}_{t=n}\≈\frac{1}{\mathrm{\Δ}t}{F}^{-1}\[\left|k_z\right|F(U^n-U^{n-1})\]\end{array}---\left(2\right)$

时间差分离散形式：

$\frac{{\∂}^{2}U}{\∂t^2}{|}_{t=n}\≈\frac{1}{{\mathrm{\Δt}}^2}(U^{n+1}+U^{n-1}-2U^n)---\left(3\right)$

空间差分离散形式：

$\begin{array}{c}\frac{{\∂}^{2}U}{\∂x^2}{|}_{x=i}\≈D_x^{2}U=\frac{1}{{\mathrm{\Δx}}^2}\[c_{0}U+\underset{m=1}{\overset{M}{\Σ}}{c}_m(U_{i+m}+u_{i-m})\]\\ \frac{{\∂}^{2}U}{\∂z^2}{|}_{z=j}\≈D_z^{2}U=\frac{1}{{\mathrm{\Δz}}^2}\[c_{0}U+\underset{m=1}{\overset{M}{\Σ}}{c}_m(U_{j+m}+u_{j-m})\]\end{array}---\left(4\right)$

带入波动方程，得到粘声各向异性介质正演模拟差分递推公式：

$\left\{\begin{array}{c}{p}^{n+1}=2p^n-p^{n-1}\±\frac{{\mathrm{\τv}}_{px}}{2}{\mathrm{\ΔtF}}^{-1}\[\left|k_x\right|F(p^n-p^{n-1})\]+{\mathrm{\Δt}}^2{v}_{px}^2{D}_x^2{p}^n+{\mathrm{\Δt}}^2{v}_{pn}{v}_{pz}{D}_z^2{q}^n)\\ q^{n+1}=2q^n-q^{n-1}\±\frac{{\mathrm{\τv}}_{pz}}{2}{\mathrm{\ΔtF}}^{-1}\[\left|k_z\right|F(q^n-q^{n-1})\]+{\mathrm{\Δt}}^2{v}_{pn}{v}_{pz}{D}_x^2{p}^n+{\mathrm{\Δt}}^2{v}_{pz}^2{D}_z^2{q}^n)\end{array}\right.---\left(5\right)$

其中，p和q分别表示水平和垂直方向的应力分量，v_pz和分别为对称轴和对称面方向的相速度，为动校正速度，ε和δ为Thomsen各向异性参数，n为时间坐标，τ＝τ_ε/τ_σ-1是一个无量纲的变量，τ_σ和τ_ε为松弛时间，可由品质因子Q求得ω为角速度，Δt为时间采样间隔，k_x和k_z分别是水平和垂直方向的波数，F和F^-1分别表示傅里叶变换和傅里叶反变换，和分别表示水平和垂直方向的空间二阶偏导数。

用U代表波场p或者q，则有限差分公式可以写为：

$\left\{\begin{array}{c}{D}_x^{2}U=\frac{1}{{\mathrm{\Δx}}^2}\[c_{0}U+\underset{m=1}{\overset{M}{\Σ}}{c}_m(U_{i+m}+u_{i-m})\]\\ D_z^{2}U=\frac{1}{{\mathrm{\Δz}}^2}\[c_{0}U+\underset{m=1}{\overset{M}{\Σ}}{c}_m(U_{j+m}+u_{j-m})\]\end{array}\right.---\left(6\right)$

其中Δx和Δz为空间采样间隔，i和j为x和z方向的空间坐标，c表示有限差分系数，M表示差分阶数；

步骤3：输入粘声各向异性探区的实际叠前炮记录；

步骤4：将炮记录由检波器处反传到整个波场，反传波场满足如下的波动方程：

$\left\{\begin{array}{c}\frac{{\∂}^{2}p}{\∂t^2}+\frac{{\mathrm{\τv}}_{px}}{2}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p}{\∂t}-v_{px}^2\frac{{\∂}^{2}p}{\∂x^2}-v_{pn}{v}_{pz}\frac{{\∂}^{2}q}{\∂z^2}=0\\ \frac{{\∂}^{2}q}{\∂t^2}+\frac{{\mathrm{\τv}}_{pz}}{2}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q}{\∂t}-v_{pn}{v}_{pz}\frac{{\∂}^{2}p}{\∂x^2}-v_{pz}^2\frac{{\∂}^{2}q}{\∂z^2}=0\end{array}\right.---\left(7\right)$

在粘声各向异性波场反向传播补偿能量衰减的过程中，解的高频成分呈现指数增长形态，会导致波场反传的不稳定，常规的处理方法是在波场延拓过程中，进行傅里叶变换，在频率域设计高频滤波器进行压制，这要求每个时刻进行一次傅里叶变换，会增加大量的计算量，而且傅里叶变换与高频滤波是在全局进行的，对于复杂模型，将会消耗掉部分有效信号。在粘声各向异性波场反向传播补偿能量衰减的过程中，在粘声各向异性拟微分方程中添加下式表示的两个方向的规则化项压制反向传播的不稳定：

$\begin{array}{c}\mathrm{\ϵ}\frac{{\mathrm{\τv}}_{px}^2}{2}\frac{\∂}{\∂t}\frac{{\∂}^{2}p}{\∂x^2}\\ \mathrm{\ϵ}\frac{{\mathrm{\τv}}_{pz}^2}{2}\frac{\∂}{\∂t}\frac{{\∂}^{2}q}{\∂z^2}\end{array}---\left(8\right)$

其中，ε表示规则化参数，为一个小正数。

得到稳定的反传过程中的波动方程：

$\left\{\begin{array}{c}\frac{{\∂}^{2}p}{\∂t^2}+\frac{{\mathrm{\τv}}_{px}}{2}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p}{\∂t}-v_{px}^2\frac{{\∂}^{2}p}{\∂x^2}-v_{pn}{v}_{pz}\frac{{\∂}^{2}q}{\∂z^2}-\mathrm{\ϵ}\frac{{\mathrm{\τv}}_{px}^2}{2}\frac{\∂}{\∂t}\frac{{\∂}^{2}p}{\∂x^2}=0\\ \frac{{\∂}^{2}q}{\∂t^2}+\frac{{\mathrm{\τv}}_{pz}}{2}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q}{\∂t}-v_{pn}{v}_{pz}\frac{{\∂}^{2}p}{\∂x^2}-v_{pz}^2\frac{{\∂}^{2}q}{\∂z^2}-\mathrm{\ϵ}\frac{{\mathrm{\τv}}_{pz}^2}{2}\frac{\∂}{\∂t}\frac{{\∂}^{2}q}{\∂z^2}=0\end{array}\right.---\left(9\right)$

步骤5：采用互相关成像条件进行成像，得到常规逆时偏移成像结果；

步骤6：对步骤5所得的成像结果应用下式表示的粘声各向异性反偏移算子进行反偏移，得到基于伯恩近似的炮记录：

$\left\{\begin{array}{c}\frac{1}{v_{pz0}^2}\frac{{\∂}^{2}dp(x,t)}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz0}}\sqrt{1+2\mathrm{\ϵ}}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂dp(x,t)}{\∂t}-(1+2\mathrm{\ϵ})\frac{{\∂}^{2}dp(x,t)}{\∂x^2}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^{2}dq(x,t)}{\∂z^2}\\ \≈m\left(x\right)(\frac{{\∂}^2{p}_0(x,t)}{\∂t^2}+\frac{{\mathrm{\τv}}_{pz0}\sqrt{1+2\mathrm{\ϵ}}}{4}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p_0(x,t)}{\∂t})\\ \frac{1}{v_{pz0}^2}\frac{{\∂}^{2}dp(x,t)}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz0}}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂p_0(x,t)}{\∂t}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^{2}dp(x,t)}{\∂x^2}-\frac{{\∂}^{2}dq(x,t)}{\∂z^2}\\ \≈m\left(x\right)(\frac{{\∂}^2{q}_0(x,t)}{\∂t^2}+\frac{{\mathrm{\τv}}_{pz0}}{4}-\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q_0(x,t)}{\∂t})\end{array}\right.---\left(10\right)$

其中，为定义的反射系数，dv_pz为速度扰动，dp(x,t)和dq(x,t)为扰动波场，x和t分别表示时间和空间坐标；

步骤6：对成像结果应用粘声各向异性反偏移算子进行反偏移；粘声各向异性波动方程可以改写为：

$\{\begin{array}{c}\frac{1}{v_{pz}^2}\frac{{\∂}^{2}p}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz}}\sqrt{1+2\mathrm{\ϵ}}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p}{\∂t}-(1+2\mathrm{\ϵ})\frac{{\∂}^{2}p}{\∂x^2}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^{2}q}{\∂z^2}=f_x\\ \frac{1}{v_{pz}^2}\frac{{\∂}^{2}w}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz}}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q}{\∂t}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^{2}p}{\∂x^2}-\frac{{\∂}^{2}q}{\∂z^2}=f_z\end{array}---\left(11\right)$

其中f_x和f_z为x方向和z方向的震源项；相速度v_pz由背景速度v_pz0及其扰动量dv_pz组成：

v_pz＝v_pz0+dv_pz(12)

波场p(x,t)和q(x,t)也由背景波场p₀(x,t)和q₀(x,t)及其扰动量dp(x,t)和dq(x,t)组成：

p(x,t)＝p₀(x,t)+dp(x,t)

(13)

q(x,t)＝q₀(x,t)+dq(x,t)

将背景速度和背景波场带入波动方程(7)中，可得:

$\{\begin{array}{c}\frac{1}{v_{pz0}^2}\frac{{\∂}^2{p}_0(x,t)}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz0}}\sqrt{1+2\mathrm{\ϵ}}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p_0(x,t)}{\∂t}-(1+2\mathrm{\ϵ})\frac{{\∂}^2{p}_0(x,t)}{\∂x^2}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^2{q}_0(x,t)}{\∂z^2}=f_x\\ \frac{1}{v_{pz0}^2}\frac{{\∂}^2{q}_0(x,t)}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz0}}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q_0(x,t)}{\∂t}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^2{p}_0(x,t)}{\∂x^2}-\frac{{\∂}^2{q}_0(x,t)}{\∂z^2}=f_z\end{array}---\left(14\right)$

应用Taylor级数展开，略去高阶项，得到如下的近似公式：

$\begin{array}{c}\frac{1}{v_{pz}^2}\≈\frac{1}{v_{pz0}^2}-\frac{2{\mathrm{dv}}_{pz}}{v_{pz0}^3}\\ \frac{1}{v_{pz}}\≈\frac{1}{v_{pz0}}-\frac{{\mathrm{dv}}_{zp}}{v_{pz0}^2}\end{array}---\left(15\right)$

将(12)(13)(15)式代入波动方程(7)中，与(14)式相减，并略去高阶项，得到粘声各向异性LSRTM的反偏移算子：

$\{\begin{array}{c}\frac{1}{v_{pz0}^2}\frac{{\∂}^{2}dp(x,t)}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz0}}\sqrt{1+2\mathrm{\ϵ}}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂dp(x,t)}{\∂t}-(1+2\mathrm{\ϵ})\frac{{\∂}^{2}dp(x,t)}{\∂x^2}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^{2}dq(x,t)}{\∂z^2}\\ =\frac{2{\mathrm{dv}}_{pz}}{v_{pz0}^3}\frac{{\∂}^2{p}_0(x,t)}{\∂t^2}+\frac{{\mathrm{dv}}_{pz}}{v_{pz0}^2}\frac{\mathrm{\τ}\sqrt{1+2\mathrm{\ϵ}}}{2}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p_0(x,t)}{\∂t}+O\left({\mathrm{dv}}_{pz}^2\right)\\ \frac{1}{v_{pz0}^2}\frac{{\∂}^{2}dq(x,t)}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz0}}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂dq(x,t)}{\∂t}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^{2}dp(d,t)}{\∂x^2}-\frac{{\∂}^{2}dq(x,t)}{\∂z^2}\\ =\frac{2{\mathrm{dv}}_{pz}}{v_{pz0}^3}\frac{{\∂}^2{q}_0(x,t)}{\∂t^2}+\frac{{\mathrm{dv}}_{pz}}{v_{pz0}^2}\frac{\mathrm{\τ}}{2}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q_0(x,t)}{\∂t}+O\left({\mathrm{dv}}_{pz}^2\right)\end{array}---\left(16\right)$

步骤7：将步骤6所得的炮记录与实际炮记录相减，判断是否满足误差条件；若：判断结果是满足误差条件，则执行步骤10；

或判断结果是不满足误差条件，则执行步骤8；

步骤8：利用下式表示的粘声各向异性梯度公式求取梯度更新方向，对步骤5所得的成像剖面进行更新；

$\begin{array}{c}g=\underset{t}{\∫}\frac{2}{v_{vp0}^2}\left(\frac{{\∂}^2{p}_0(x,t)}{\∂t^2}{p}^R\right(x,t)+\frac{{\∂}^2{q}_0(x,t)}{\∂t^2}{q}^R(x,t\left)\right)\\ +\frac{1}{v_{vp0}}\left(\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p_0(x,t)}{\∂t}{p}^R\right(x,t)+\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q_0(x,t)}{\∂t}{q}^R(x,t\left)\right)dt\end{array}$

其中，g为梯度，p₀(x,t)和q₀(x,t)为各向异性背景波场，p^R(x,t)和q^R(x,t)分别为p₀(x,t)和q₀(x,t)的伴随变量。

步骤9：将步骤8所得的成像剖面返回步骤6进行反偏移；

步骤10：输出最终的成像剖面。

本发明所带来的有益技术效果：

本发明提出了一种粘声各向异性介质最小二乘逆时偏移成像方法，与现有技术相比，本发明能够同时考虑地下介质的粘滞性和各向异性，通过提供一种稳定的粘声各向异性最小二乘偏移算子和反偏移算子，既能克服传统成像方法在同时处理同时存在粘滞性和各向异性介质的不足，又通过最小二乘方式消除了成像噪音，得到真振幅的成像剖面，开发基于粘声各向异性介质的最小二乘逆时偏移成像技术，为同时存在粘滞性和各向异性探区的采集的地震数据进行后续的解释工作提供高精度成像基础。

附图说明

图1为本发明粘声各向异性介质最小二乘逆时偏移成像方法的流程框图。

图2为偏移速度模型。

图3为品质因子Q模型。

图4为各向异性Thomsen参数Epsilon模型。

图5为各向异性Thomsen参数Delta模型。

图6为相对慢度扰动模型。

图7为粘声各向异性常规逆时偏移成像结果。

图8为声波各向同性最小二乘逆时偏移成像结果。

图9为粘声各向同性最小二乘逆时偏移成像结果。

图10为声波各向异性最小二乘逆时偏移成像结果。

图11为本发明的偏移成像结果。

图12为不同成像方法模型的误差曲线。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

附图1为本发明的粘声各向异性介质最小二乘逆时偏移成像方法的流程图，粘声各向异性介质最小二乘逆时偏移成像方法，包括如下步骤：

输入速度模型、各向异性Thomsen参数模型、品质因子Q模型及慢度扰动模型，并建立观测系统；

采用有限差分方法进行粘声各向异性介质正演模拟；

输入粘声各向异性探区的实际叠前炮记录；

将炮记录由检波器处反传到整个波场，在粘声各向异性波场反向传播补偿能量衰减的过程中，在粘声各向异性拟微分方程中添加两个规则化项压制反向传播的不稳定；

采用互相关成像条件进行成像，得到常规逆时偏移成像结果；

对成像结果应用粘声各向异性反偏移算子进行反偏移，得到基于伯恩近似的炮记录，并实际炮记录相减，判断是否满足误差条件，如果满足条件输出最终的成像剖面，如果不满足条件求取梯度更新方向，对成像剖面进行更新，再次进行反偏移，直至满足误差条件。

将本发明粘声各向异性介质的最小二乘成像方法，应用于国际标准的Marmousi模型数据，取得了理想的计算效果。输入速度模型(如图2所示)、品质因子Q模型(如图3所示)、各向异性Thomsen参数Epsilon模型(如图4所示)、各向异性Thomsen参数Delta模型(如图5所示)以及慢度扰动模型(如图6所示)，并建立观测系统；波场正演模拟和反向传播，并采用互相关成像条件进行成像，得到常规逆时偏移成像结果(如图7所示)；对常规逆时偏移成像结果应用粘声各向异性反偏移算子进行反偏移，得到基于伯恩近似的炮记录，并实际炮记录相减，判断是否满足误差条件，如果不满足条件求取梯度更新方向，对成像剖面进行更新，再次进行反偏移，直至满足误差条件，如果满足条件输出最终的偏移成像结果(如图11所示)。本发明的粘声各向异性最小二乘逆时偏移方法所得到的成像剖面，相比于声波各向同性最小二乘偏移结果(如图8所示)和粘声各向同性最小二乘偏移结果(如图9所示)，其反射同相轴能够正确归位，成像位置较为清晰，由相位畸变等造成的同相轴紊乱现象基本被消除，剖面中震源效应被抑制，断层、背斜、高陡构造以及不整合面等地下复杂构造准确成像，相比于声波各向异性最小二乘逆时偏移成像结果(如图10所示)，中深层反射能量也得到补偿，相对于粘声各向异性常规逆时偏移成像结果(如图7所示)，深部能量强，成像剖面整体均衡化，得到了更为精确可靠的保幅成像剖面。

图12为不同成像方法模型的误差曲线，本发明的成像方法在迭代过程中的误差曲线相比于其他成像方法收敛速度更快。

地震勘探技术向复杂构造地区不断发展，因此，对存在粘滞性和各向异性介质的波场传播性质进行研究，提高存在粘滞性和各向异性介质的成像精度对适应国内地震勘探的发展趋势是非常必要的。

为此本发明提供了一种粘声各向异性最小二乘逆时偏移成像技术，开发基于同时存在粘滞性和各向异性介质的高精度成像技术，为后续的复杂地质构造解释工作提供成像基础。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (3)

1.一种粘声各向异性介质最小二乘逆时偏移成像方法，其特征在于：包括如下步骤：

步骤1：输入速度模型、各向异性Thomsen参数模型、品质因子Q模型及慢度扰动模型，并建立观测系统；

步骤2：采用有限差分方法进行粘声各向异性介质正演模拟，得到下式表示的粘声各向异性介质正演模拟差分递推公式；

$\left\{\begin{array}{c}{p}^{n+1}=2p^n-p^{n-1}\±\frac{{\mathrm{\τv}}_{px}}{2}\mathrm{\Δ}t{F}^{-1}\[\left|k_x\right|F(p^n-p^{n-1})\]+{\mathrm{\Δt}}^2{v}_{px}^2{D}_x^2{p}^n+{\mathrm{\Δt}}^2{v}_{pn}{v}_{pz}{D}_z^2{q}^n)\\ q^{n+1}=2q^n-q^{n-1}\±\frac{{\mathrm{\τv}}_{pz}}{2}\mathrm{\Δ}t{F}^{-1}\[\left|k_z\right|F(q^n-q^{n-1})\]+{\mathrm{\Δt}}^2{v}_{pn}{D}_x^2{p}^n+{\mathrm{\Δt}}^2{v}_{pz}^2{D}_z^2{q}^n)\end{array}\right.$

其中，p和q分别表示水平和垂直方向的应力分量，v_pz和分别为对称轴和对称面方向的相速度，为动校正速度，ε和δ为Thomsen各向异性参数，n为时间坐标，τ＝τ_ε/τ_σ-1是一个无量纲的变量，τ_σ和τ_ε为松弛时间，可由品质因子Q求得：ω为角速度，Δt为时间采样间隔，k_x和k_z分别是水平和垂直方向的波数，F和F^-1分别表示傅里叶变换和傅里叶反变换，和分别表示水平和垂直方向的空间二阶偏导数；

用U代表p或者q：

$\left\{\begin{array}{c}{D}_x^{2}U=\frac{1}{{\mathrm{\Δx}}^2}\[c_{0}U+\underset{m=1}{\overset{M}{\Σ}}{c}_m(U_{i+m}+u_{i-m})\]\\ D_z^{2}U=\frac{1}{{\mathrm{\Δz}}^2}\[c_{0}U+\underset{m=1}{\overset{M}{\Σ}}{c}_m(U_{j+m}+U_{j-m})\]\end{array}\right.$

其中，Δx和Δz为空间采样间隔，i和j为x和z方向的空间坐标，c表示有限差分系数，M表示差分阶数；

步骤3：输入粘声各向异性探区的实际叠前炮记录；

步骤4：将实际叠前炮记录由检波器处反传到整个波场，在粘声各向异性波场反向传播补偿能量衰减的过程中，在粘声各向异性拟微分方程中添加两个规则化项压制反向传播的不稳定；

步骤5：采用互相关成像条件进行成像，得到常规逆时偏移成像结果；

步骤6：对步骤5所得的常规逆时偏移成像结果，应用下式表示的粘声各向异性反偏移算子进行反偏移，得到基于伯恩近似的炮记录；

$\left\{\begin{array}{c}\frac{1}{v_{pz0}^2}\frac{{\∂}^{2}dp(x,t)}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz0}}\sqrt{1+2\mathrm{\ϵ}}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂dp(x,t)}{\∂t}-(1+2\mathrm{\ϵ})\frac{{\∂}^{2}dp(x,t)}{\∂x^2}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^{2}dq(x,t)}{\∂z^2}\\ =m\left(x\right)\left(\frac{{\∂}^2{p}_0(x,t)}{\∂t^2}+\frac{{\mathrm{\τv}}_{pz0}\sqrt{1+2\mathrm{\ϵ}}}{4}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p_0(x,t)}{\∂t}\right)\\ \frac{1}{v_{pz0}^2}\frac{{\∂}^{2}dq(x,t)}{\∂t^2}+\frac{\mathrm{\τ}}{2}\frac{1}{v_{pz0}}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂dp(x,t)}{\∂t}-\sqrt{1+2\mathrm{\δ}}\frac{{\∂}^{2}dp(x,t)}{\∂x^2}-\frac{{\∂}^{2}dq(x,t)}{\∂z^2}\\ =m\left(x\right)\left(\frac{{\∂}^2{q}_0(x,t)}{\∂t^2}+\frac{{\mathrm{\τv}}_{pz0}}{4}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂q_0(x,t)}{\∂t}\right)\end{array}\right.$

其中，为定义的反射系数，dv_pz为速度扰动，dp(x,t)和dq(x,t)为扰动波场，x和t分别表示时间和空间坐标；

步骤7：将步骤6所得的基于伯恩近似的炮记录与实际叠前炮记录相减，判断是否满足误差条件；

若：判断结果是满足误差条件，则执行步骤8；

或判断结果是不满足误差条件，利用下式表示的粘声各向异性梯度公式求取梯度更新方向，对步骤5所得的成像结果进行更新，然后执行步骤6；

$\begin{array}{c}g=\underset{t}{\∫}\frac{2}{v_{vp0}^2}\left(\frac{{\∂}^2{p}_0(x,t)}{\∂t^2}{p}^R\right(x,t)+\frac{{\∂}^2{q}_0(x,t)}{\∂t^2}{q}^R(x,t\left)\right)\\ +\frac{1}{v_{vp0}}\left(\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p_0(x,t)}{\∂t}{p}^R\right(x,t)+\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q_0(x,t)}{\∂t}{q}^R(x,t))dt\end{array}$

其中，g为梯度，p₀(x,t)和q₀(x,t)为各向异性背景波场，p^R(x,t)和q^R(x,t)分别为p₀(x,t)和q₀(x,t)的伴随变量；

步骤8：输出最终的成像结果。

2.根据权利要求1所述的粘声各向异性介质最小二乘逆时偏移成像方法，其特征在于：在步骤4中，所述反传波场满足如下的波动方程：

$\left\{\begin{array}{c}\frac{{\∂}^{2}p}{\∂t^2}-\frac{{\mathrm{\τv}}_{px}}{2}\sqrt{-\frac{{\∂}^2}{\∂x^2}}\frac{\∂p}{\∂t}-v_{px}^2\frac{{\∂}^{2}p}{\∂x^2}-v_{pn}{v}_{pz}\frac{{\∂}^{2}q}{\∂z^2}=0\\ \frac{{\∂}^{2}q}{\∂t^2}-\frac{{\mathrm{\τv}}_{pz}}{2}\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q}{\∂t}-v_{pn}{v}_{pz}\frac{{\∂}^{2}p}{\∂x^2}-v_{pz}^2\frac{{\∂}^{2}q}{\∂z^2}=0\end{array}\right..$

3.根据权利要求1所述的粘声各向异性介质最小二乘逆时偏移成像方法，其特征在于：在步骤4中，所述在粘声各向异性拟微分方程中添加的两个规则化项为：

$\mathrm{\ϵ}\frac{{\mathrm{\τv}}_{px}^2}{2}\frac{\∂}{\∂t}\frac{{\∂}^{2}p}{\∂x^2},$

$\mathrm{\ϵ}\frac{{\mathrm{\τv}}_{pz}^2}{2}\frac{\∂}{\∂t}\frac{{\∂}^{2}q}{\∂z^2}$

其中，ε表示规则化参数。