[发明专利]一种粘声各向异性最小二乘逆时偏移成像方法

摘要

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

主权项

一种粘声各向异性介质最小二乘逆时偏移成像方法，其特征在于：包括如下步骤：步骤1：输入速度模型、各向异性Thomsen参数模型、品质因子Q模型及慢度扰动模型，并建立观测系统；步骤2：采用有限差分方法进行粘声各向异性介质正演模拟，得到下式表示的粘声各向异性介质正演模拟差分递推公式；$\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.$其中，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：将实际叠前炮记录由检波器处反传到整个波场，在粘声各向异性波场反向传播补偿能量衰减的过程中，在粘声各向异性拟微分方程中添加两个规则化项压制反向传播的不稳定；所述在粘声各向异性拟微分方程中添加的两个规则化项为：$\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}$其中，ε表示规则化参数；步骤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{\∂dq(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}{\∂z^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\left(x,t\right)+\sqrt{-\frac{{\∂}^2}{\∂z^2}}\frac{\∂q_0(x,t)}{\∂t}{q}^R\left(x,t\right)\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)的伴随变量；步骤8：输出最终的成像结果。