CN105204064A - 一种基于优化系数的混合域傅里叶有限差分偏移方法 - Google Patents

Abstract

为了提高复杂高陡构造区地震偏移成像精度，有力指导油气资源的准确预测及开采，本发明提出一种基于优化系数的混合域傅里叶(Fourier)有限差分叠前深度偏移成像方法，该方法利用padé近似的有理函数对波场外推算子进行展开，然后利用切比雪夫(Chebyshev)多项式优化展开式系数，推导得到新的波场外推算子，降低了与波动方程精确波场外推算子的相对误差，提高了对波场外推算子的逼近程度，且在保证计算效率的同时提高了复杂高陡构造区地震偏移成像的精度。对比改进后的混合域傅立叶有限差分偏移方法与常规的傅立叶有限差分偏移方法(FFD)对Marmousi模型偏移剖面的成像效果，证明了系数优化后的混合域Fourier有限差分叠前深度偏移方法具有更高的偏移成像精度。该方法对于改善复杂高陡构造区地震偏移成像精度具有重要的意义。

Description

一种基于优化系数的混合域傅里叶有限差分偏移方法

技术领域

本发明涉及地震勘探数据资料处理领域，特别是涉及叠前深度偏移成像处理技术。

背景技术

常规地震偏移成像方法是基于均匀各向同性弹性介质的假设条件，但随着我国油气勘探的不断深入，其勘探难度不断增加，复杂高陡构造区成为当今油气勘探的重点，在该区采集的地震原始资料很难同时满足常规地震偏移成像方法的假设条件。波动方程叠前深度偏移作为处理复杂高陡构造区最有效的方法迅速发展起来。

波动方程偏移成像技术始于J.F.Claerbout团队在1971年与1972年提出的基于有限差分近似解的波动方程偏移。波动方程叠前深度偏移的核心是波场延拓算子，常规波动方程叠前深度偏移方法是基于单程波动方程，波动方程偏移算子的优化构造既可以在频率-空间域中实现，也可以在频率-波数域实现，在不同域优化构造的波动方程偏移延拓算子对介质空间变化的适应性不同：在频率-空间域中优化构造的偏移成像延拓算子能够适应介质的空间变化，而在频率-波数域优化构造的波动方程偏移成像算子仅仅能适应空间上的均匀介质。从两种域中的波动方程偏移延拓算子适应性可以看出，单域中的波动方程偏移成像延拓算子具有一定的局限性，受地形起伏、低降速带以及横向剧烈变速的影响，单域中的波动方程偏移成像方法难以满足复杂构造区的精确成像。

由于单域中波动方程偏移成像方法的限制性，自20世界90年代以来，研究人员尝试结合频率-空间域和频率-波数域中波动方程偏移成像延拓算子各自的优点，在混合域(频率-空间域及频率-波数域)中优化构造波动方程偏移成像延拓算子，并且取得了一定的进展。目前应用比较广泛的混合域波动方程偏移成像延拓算子主要包括有分裂步傅立叶算子(SP)、傅立叶有限差分算子(FFD)以及广义屏算子(GSP)等，参考文献有：StoffaPL,FokkemaJT,FreireRM,etal.Split-stepFourierMigration.Geophysics,1990,55(4):410-421；RistowD,RuhlT.Fourierfinite-differencemigration.Geophysics,1994,59(12):1882-1893；LeRousseauJH,HoopMV.Modelingandimagingwiththescalargeneralized-screenalgorithmsinisotropicmedia.Geophysics,2001,66(5):1551-1568.从本质上来说，FFD传播算子和GPS传播算子都是SSF传播算子(也可称为相位屏传播算子)在不同视角上的推广。由于在传播算子的优化构造中应用了Fourier变换，使其在空间域和波数域这两种对偶空间中进行变换，因此混合域波动方程偏移成像延拓算子能够依据介质的空间变化情况自适应地保持频率-波数域波动方程偏移成像延拓算子和频率-空间域波动方程偏移成像延拓算子的特点，尤其是广义屏传播算子还能有效的克服频率-空间域波动方程偏移成像延拓算子对某些叠前道集数据和三维情况存在的一些数值计算方面的麻烦。

在以上研究的基础上，为了进一步提高混合域中波动方程偏移成像方法的成像精度，本发明提出一种基于优化系数的混合域傅立叶(Fourier)有限差分叠前深度偏移成像方法，该方法利用padé近似的有理函数对波场外推算子进行展开，然后利用切比雪夫(Chebyshev)多项式优化展开式系数，推导得到新的波场外推算子，降低了与波动方程精确波场外推算子的相对误差，提高了对波场外推算子的逼近程度。利用本发明获得的偏移成像剖面，在保证计算效率的同时提高了复杂高陡构造区地震偏移成像的精度，能较好地适应介质的空间变化。

发明内容

为了解决单域(频率-空间域或频率-波数域)中波动方程偏移成像方法受限制性问题，本发明目的是在混合域(频率-空间域及频率-波数域)中通过数值计算的方法优化构造波动方程偏移成像延拓算子，使得优化构造的混合域波动方程偏移成像延拓算子与精确偏移成像延拓算子的误差减小，从而提高复杂构造区偏移成像剖面的成像精度，有力指导复杂构造区油气的勘探与开采。

本发明的创新之处在于利用padé近似的有理函数逼近混合域中的波动方程偏移成像延拓算子，并利用切比雪夫多项式函数辅以优化padé近似的有理函数逼近展开式系数，从而推导出新的波动方程偏移成像延拓算子，使优化构造的波动方程偏移成像延拓算子与精确的单平方根算子的误差减小，从而提高复杂高陡构造区的地震偏移成像精度。

本发明是一种基于优化系数的混合域傅立叶(Fourier)有限差分叠前深度偏移成像方法，其特征在于：将波动方程在频率-空间域通过傅立叶变换到频率-波数域，然后在频率-波数域利用padé近似的有理函数逼近波动方程偏移成像延拓算子，最后利用Chebyshev多项式对padé有理函数逼近式系数进行优化，得到新的混合域波动方程偏移成像延拓算子，利用计算得到的新的混合域波动方程偏移成像延拓算子对模型进行偏移成像。

本发明方法的特征在于它包含下列步骤：

1)利用傅里叶变换的微分性质，将二维全波列方程从时间—空间域变换到频率—波数域，表达式如下：

$\{k_z+{\[\frac{{\mathrm{\ω}}^2}{v^2(x,z)}-k_x^2\]}^{\frac{1}{2}}\}\·\{k_z-{\[\frac{{\mathrm{\ω}}^2}{v^2(x,z)}-k_x^2\]}^{\frac{1}{2}}\}\·u(k_x,k_z;\mathrm{\ω})=0---\left(1\right)$

其中：k_x和k_z分别为x、z方向上的波数，单位为：/m；v为介质速度，单位为：m/s；ω为角频率，单位为：rad/s；u(ω；k_x,k_z)为频率—波数域地震波场值；

2)将(1)式进行因式分解，得到z方向上的波数k_z，表达式如下：

$k_z={\[\frac{{\mathrm{\ω}}^2}{v^2(x,z)}-k_x^2\]}^{\frac{1}{2}}=\frac{\mathrm{\ω}}{v(x,z)}{\[1-\frac{v^2(x,z)}{{\mathrm{\ω}}^2}\·k_x^2\]}^{\frac{1}{2}}---\left(2\right)$

在实际地震资料中，介质速度分为常数背景速度和扰动速度，设背景速度为v₀(x,z)，式(2)可化为：

$k_z={\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}-k_x^2\]}^{\frac{1}{2}}+{\[\frac{{\mathrm{\ω}}^2}{v^2(x,z)}-k_x^2\]}^{\frac{1}{2}}-{\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}-k_x^2\]}^{\frac{1}{2}}---\left(3\right)$

则单平方根算子的误差为:

$\mathrm{\ξ}={\[\frac{{\mathrm{\ω}}^2}{v^2(x,z)}-k_x^2\]}^{\frac{1}{2}}-{\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}-k_x^2\]}^{\frac{1}{2}}---\left(4\right)$

3)利用padé有理函数对式(2)逼近，设R_(m,n)(d)是f(d)的最佳逼近式，则有：

$R_{(m,n)}\left(d\right)=\frac{a_0{T}_0\left(d\right)+a_1{T}_1\left(d\right)+...+a_m{T}_m\left(d\right)}{T_0\left(d\right)+b_1{T}_1\left(d\right)+...+b_n{T}_n\left(d\right)}---\left(5\right)$

式中，a₀,a₁,…,a_m；b₁,…,b_n为常系数，T₀(d),T₁(d),…为第一类切比雪夫多项式，且在式中d＝v(x,z)·k_x/ω_，d∈(-1,1)。考虑到算子的计算效率，在保证成像精度的基础上，算子阶数不能太高，因此在下面算子计算过程中都取m＝n＝1，则有：

$R_{(1,1)}\left(d\right)=\frac{a_0{T}_0\left(d\right)+a_1{T}_1\left(d\right)}{T_0\left(d\right)+b_1{T}_1\left(d\right)}---\left(6\right)$

4)利用Chebyshev多项式对padé有理函数逼近式(6)中的系数a₀,a₁,b₁进行优化，取Chebyshev多项式前三项，与padé有理函数逼近式R_(1,1)(d)相等，可得到二阶近似，进而求得a₀,a₁,b₁满足的线性代数方程组为：

$\left\{\begin{array}{c}{a}_0=c_0+\frac{1}{2}{b}_1{c}_1\\ a_1=c_1+b_1(c_0+\frac{1}{2}{c}_2)\\ 0=\frac{1}{2}{b}_1(c_1+c_3)+c_2\end{array}\right.---\left(3\right)$

式(6)中c₀,c₁,c₂为Chebyshev多项式系数。求解式(7)可得a₀,a₁,b₁，将其带入(6)式得：

$R_{(1,1)}\left(d\right)=\frac{(c_0-c_1\frac{c_2}{c_1+c_3})T_0\left(d\right)+(c_1-c_2\frac{(2c_0+c_2)}{c_1+c_3})T_1\left(d\right)}{T_0\left(d\right)-\left(2\frac{2c_0}{c_1+c_3}\right)T_1\left(d\right)}---\left(8\right)$

由于f(d)为偶函数，在求解f(d)的Chebyshev多项式的系数时，奇数项系数c₁＝c₃＝…＝c_2M-1＝0，此时求出的系数a₀,a₁,b₁都等于零，使得式(6)没有意义，因此需对f(d)的Chebyshev多项式中的奇数项进行剔除，只保留偶数项的Chebyshev多项式，则式(8)可化为：

$R_{(1,1)}\left(d\right)=\frac{a_0+a_1(2d^2-1)}{1+b_1(2d^2-1)}=1+\frac{a_0-1+(a_1-b_1)(2d^2-1)}{1+b_1(2d^2-1)}=1-\frac{1-a_0+a_1-b_1-2(a_1-b_1)d^2}{1-b_1+2b_1{d}^2}---\left(9\right)$

对式(9)进行简化_，令A＝1-a₀+a₁-b₁，B＝2(a₁-b₁)，C＝1-b₁，D＝2b₁，则式(9)可简化为：

$R_{(1,1)}\left(d\right)=1-\frac{A-{\mathrm{Bd}}^2}{C+{\mathrm{Dd}}^2}---\left(10\right)$

5)计算混合域中单平方根偏移成像延拓算子，将(10)式代入(4)，表达式如下：

$\mathrm{\ξ}=(\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)})+\frac{\mathrm{\ω}}{v(x,z)}\[-\frac{A-B\frac{v^2(x,z)}{{\mathrm{\ω}}^2}{k}_x^2}{C+D\frac{v^2(x,z)}{{\mathrm{\ω}}^2}{k}_x^2}\]-\frac{\mathrm{\ω}}{v_0(x,z)}\[-\frac{A-B\frac{v_0^2(x,z)}{{\mathrm{\ω}}^2}{k}_x^2}{C+D\frac{v_0^2(x,z)}{{\mathrm{\ω}}^2}{k}_x^2}\]---\left(11\right)$

将式(11)代入式(3)，其表达式如下：

$k_z={\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}-k_x^2\]}^{\frac{1}{2}}+\frac{\mathrm{\ω}}{v(x,z)}\[1-\frac{A-B\frac{v^2(x,z)}{{\mathrm{\ω}}^2}{k}_x^2}{C+D\frac{v^2(x,z)}{{\mathrm{\ω}}^2}{k}_x^2}\]-\frac{\mathrm{\ω}}{v_0(x,z)}\[1-\frac{A-B\frac{v_0^2(x,z)}{{\mathrm{\ω}}^2}{k}_x^2}{C+D\frac{v_0^2(x,z)}{{\mathrm{\ω}}^2}{k}_x^2}\]---\left(12\right)$

根据傅里叶变换的微分性质，利用傅里叶反变换将式(12)变换到频率—空间域，得单程波动方程为：

$\begin{array}{c}\frac{\∂u(x,z;\mathrm{\ω})}{\∂z}=-i\[{\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}-k_x^2\]}^{\frac{1}{2}}+(\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)})+\frac{\mathrm{\ω}}{v(x,z)}\[-\frac{A+B\frac{v^2(x,z)}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2}}{C-D\frac{v^2(x,z)}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2}}\]\\ -\frac{\mathrm{\ω}}{v_0(x,z)}\[-\frac{A+B\frac{v_0^2(x,z)}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2}}{C-D\frac{v_0^2(x,z)}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2}}\]u(x,z;\mathrm{\ω})\end{array}---\left(13\right)$

则混合域中偏移成像延拓算子Λ可表示为：

$\begin{array}{c}\mathrm{\Λ}={\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}+\frac{{\∂}^2}{\∂x^2}\]}^{\frac{1}{2}}+(\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)})+\[A\frac{\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)}}{C+D\frac{v^2(x,z)+v_0^2(x,z)}{{\mathrm{\ω}}^2}\·\frac{{\∂}^2}{\∂x^2})}\\ +B\frac{\frac{v(x,z)-v_0(x,z)}{\mathrm{\ω}}\·\frac{{\∂}^2}{\∂x^2}}{C+D\frac{v^2(x,z)+v_0^2(x,z)}{{\mathrm{\ω}}^2}\·\frac{{\∂}^2}{\∂x^2})}\]=A_1+A_2+A_3\end{array}---\left(14\right)$

由式(14)可知，混合域中偏移成像延拓算子可分解为三个算子A₁,A₂,A₃，其中是频率—波数域关于常数背景速度的相移偏移算子；表示频率—空间域的时移校正算子，相当于分步傅里叶法(SSF)里的一阶校正项； $A_3=\[A\frac{\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)}}{C+D\frac{v^2(x,z)+v_0^2(x,z)}{{\mathrm{\ω}}^2}\·\frac{{\∂}^2}{\∂x^2})}+B\frac{\frac{v(x,z)-v_0(x,z)}{\mathrm{\ω}}\·\frac{{\∂}^2}{\∂x^2}}{C+D\frac{v^2(x,z)+v_0^2(x,z)}{{\mathrm{\ω}}^2}\·\frac{{\∂}^2}{\∂x^2})}\]$ 表示频率—空间域扰动波场的差分算子，则混合域中波动方程外推过程为：

$u(x,z+\mathrm{\Δ}z\·,\mathrm{\ω})=u(x,z;\mathrm{\ω})\·e^{-{\mathrm{ik}}_z\·\mathrm{\Δ}z}=u(x,z;\mathrm{\ω})\·e^{-{\mathrm{iA}}_1\·\mathrm{\Δ}z}\·e^{-{\mathrm{iA}}_2\·\mathrm{\Δ}z}\·e^{-{\mathrm{iA}}_3\·\mathrm{\Δ}z}---\left(15\right)$

6)按照(15)式并以相关型成像原理对地震数据进行偏移延拓成像，获得最终的偏移成像剖面。

本发明方法的依据是：

1)叠前深度偏移是解决复杂构造区偏移成像最有效的方法，能够对速度横向剧烈变化和构造高陡区地震资料精确成像。

2)混合域中波动方程叠前深度偏移成像方法结合了单域中波动方程叠前深度偏移方法的优点，能够适应介质的空间变化，可有效克服单域中波动方程叠前深度偏移成像方法的限制性问题。

本发明具有以下优点：

1)利用padé近似的有理函数逼近混合域中的波动方程偏移成像延拓算子，并利用切比雪夫多项式函数辅以优化padé近似的有理函数逼近展开式系数，从而推导出新的波动方程偏移成像延拓算子，其与精确的波动方程偏移成像延拓算子的误差相对减小。

2)利用新的优化偏移成像延拓算子进行偏移成像，提高了偏移成像精度，偏移剖面上构造形态清晰，提高了构造解释的准确性，利于指导山地高陡构造区油气藏的勘探开发。

3)混合域中波动方程叠前深度偏移成像方法结合了单域中波动方程叠前深度偏移方法的优点，能够适应介质的空间变化，可有效克服单域中波动方程叠前深度偏移成像方法的限制性问题。

附图说明

图1为算法实现流程图；

图2优化构造的新算子与精确算子的相对误差期望曲线图；横坐标为角度(θ)，纵坐标为不同优化算子误差值Er(θ)；在算子相对误差计算过程中，计算角度个数N＝86.0，图中红色曲线表示的是连分式式展开的傅立叶有限差分(FFD)算子误差，误差期望值为E＝0.78547，黑色曲线表示的是新算子误差，误差期望值为E＝0.04483；

图3为Marmousi模型速度场；该图横坐标表示偏移距，纵坐标表示深度；Marmousi模型构造复杂，尤其是在偏移距为4000m-7000m之间是成像的难点，横向速度变化剧烈，该模型的速度场参数为：横向采样点数为737，纵向采样点数为750，横向采样间隔为12.5m，纵向采样间隔为4m；

图4基于连分式展开的傅立叶有限差分算子的叠前深度偏移成像剖面；该图纵坐标表示深度，横坐标表示偏移距；整体上从图中可看出构造的基本形态，但近偏移距的成像效果差，同相轴连续性差，且断层构造的断面不连续，断点模糊不清(如图中白色框中黑色箭头所示)，另外在叫深层的断层断面一下同相轴模糊不清，无法识别地层的形态(如图中黑色框中所示)。

图5为本发明方法获得的混合域叠前深度偏移成像剖面；该图横坐标表示偏移距，纵坐标表示深度；对比图4可以看出，整体上剖面的信噪比和分辨率得到了提高，构造形态清晰，断层的连通性较好，断层的断面连续，断点清晰(如图中白色框中黑色箭头所示)，深层地层形态也能较清晰地识别。

具体实施方式：

下面根据附图对本发明的技术方案的主要实现原理、具体实施方式等进行详细描述。

假设频散方程单平方根算子的最佳逼近式为：

$R_{(1,1)}\left(d\right)=\frac{a_0+a_1(2d^2-1)}{1+b_1(2d^2-1)}=1+\frac{a_0-1+(a_1-b_1)(2d^2-1)}{1+b_1(2d^2-1)}=1-\frac{1-a_0+a_1-b_1-2(a_1-b_1)d^2}{1-b_1+2b_1{d}^2}---\left(16\right)$

式(16)中：d＝v(x,z)·k_x/ω_，d∈(-1,1)_，a₀,a₁,b₁为待求的系数。

对(16)式进行有理切比雪夫逼近，可得系数a₀＝1.358，1a₁＝1.2866，b₁＝0.3684，则逼近式(16)表示为：

$R_{(1,1)}\left(x\right)=\frac{1.3581-1.2866{\left(\frac{v(x,z)k_x}{\mathrm{\ω}}\right)}^2}{1.3684-0.7368{\left(\frac{v(x,z)k_x}{\mathrm{\ω}}\right)}^2}---\left(17\right)$

设波在介质中的传播角度为θ，则则式(17)与精确频散方程的误差为：

$E_r\left(\mathrm{\θ}\right)=\frac{\frac{1.3581-1.2866{\sin }^2\mathrm{\θ}}{1.3684-0.7368{\sin }^2\mathrm{\θ}}-cos\mathrm{\θ}}{\cos \mathrm{\θ}}---\left(18\right)$

则频率-空间域的优化系数的频散方程单平方根算子为：

$\begin{array}{c}\mathrm{\Λ}={\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}+\frac{{\∂}^2}{\∂x^2}\]}^{\frac{1}{2}}+(\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)})\\ \frac{0.5601\×\[\frac{\mathrm{\ω}}{v_0(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)}\]+1.8364\×\frac{v(x,z)-v_0(x,z)}{\mathrm{\ω}}\·\frac{{\∂}^2}{\∂x^2}}{0.6316+0.7368\×\frac{v^2(x,z)+v_0^2(x,z)}{{\mathrm{\ω}}^2}\·\frac{{\∂}^2}{\∂x^2})}\end{array}---\left(19\right)$

则优化系数的混合域叠前深度偏移的下行波的外推过程可表示为：

$\begin{array}{c}u(x,z+\mathrm{\Δ}z;\mathrm{\ω})=u(x,z;\mathrm{\ω})\·e^{-{\mathrm{ik}}_z\·\mathrm{\Δ}z}\\ =u(x,z;\mathrm{\ω})\·e^{-i{\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}+\frac{{\∂}^2}{\∂x^2}\]}^{\frac{1}{2}}\·\mathrm{\Δ}z}\·e^{-i(\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)})\·\mathrm{\Δ}z}\·e^{-i\frac{0.5601\×\[\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)}\]+1.8364\×\frac{v(x,z)-v_0(x,z)}{\mathrm{\ω}}\·\frac{{\∂}^2}{\∂x^2}}{0.6316+0.7368\×\frac{v^2(x,z)+v_0^2(x,z)}{{\mathrm{\ω}}^2}\·\frac{{\∂}^2}{\∂x^2})}\·\mathrm{\Δ}z}\end{array}---\left(20\right)$

波场从深度z到z+Δz的下延拓可通过三步来完成，第一步是在频率—波数域中计算相移算子 ${\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}+\frac{{\∂}^2}{\∂x^2}\]}^{\frac{1}{2}},$ 余下两步是分别在频率空间域中用有限差分法计算偏移算子和

为了获得偏移成像剖面，采用本发明的偏移成像方法对二维Marmousi模型进行叠前深度偏移成像。横向采样点数为737，纵向采样点数为750，横向采样间隔为12.5m，纵向采样间隔为4m。图3为二维盐丘模型的速度模型。对比图4、图5，本发明提出偏移方法取得很好的成像结果。图5成像效果相对于图4有了明显的改善，图5的断层构造形态更清晰，构造边界更明显，断面连续，断点可以清晰识别。从整体上看，本发明提出的基于优化系数的混合域傅立叶有限差分叠前深度偏移成像方法对强横向速度变化的介质具有更高的描述精度。

Claims (2)

1.本发明涉及一种基于优化系数的混合域傅立叶有限差分叠前深度偏移方法，其特征在于包括如下步骤：

步骤1：将二维声波方程变换到频率-波数域，并用波场的背景速度v₀(x,z)代替介质的实际速度v(x,z)造成的频散方程的单平方根的误差，其误差如(1)式所示：

$\mathrm{\ξ}={\[\frac{{\mathrm{\ω}}^2}{v^2(x,z)}-k_x^2\]}^{\frac{1}{2}}-{\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}-k_x^2\]}^{\frac{1}{2}}---\left(1\right)$

(1)式中：x、z为方向坐标，k_x，k_z分别为x、z方向上的波数，单位为：/m，v₀(x,z)为波场的背景速度，单位为：m/s，v(x,z)为波场的实际速度，单位为：m/s，ω为角频率，单位为：rad/s；

设d＝v(x,z)·k_x/ω，d∈(-1,1)，则(1)式中根式可以写为：

$f\left(d\right)={(1-d^2)}^{\frac{1}{2}}---\left(2\right)$

步骤2：利用padé有理函数逼近(2)式，表达式如(3)式所示：

(3)式中：d＝v(x,z)·k_x/ω，d∈(-1,1)，a₀，a₁，b₁为逼近式系数，T₀(d)，T₁(d)为第一类切比雪夫多项式；

步骤3利用切比雪夫多项式函数对(3)式中的系数a₀，a₁，b₁进行优化，取切比雪夫多项式函数的前三项与(3)式相等，系数a₀＝1.358，1a₁＝1.2866，b₁＝0.3684，则(3)式可表示为：

$R_{(1,1)}\left(x\right)=\frac{1.3581-1.2866{\left(\frac{v(x,z)k_x}{\mathrm{\ω}}\right)}^2}{1.3684-0.7368{\left(\frac{v(x,z)k_x}{\mathrm{\ω}}\right)}^2}---\left(4\right)$

步骤4将(4)式代入到(1)式并变换到频率-空间域，代入到优化系数的混合域叠前深度偏移下行波外推方程中，如(5)式所示：

$\begin{array}{c}u(x,z+\mathrm{\Δ}z;\mathrm{\ω})=u(x,z;\mathrm{\ω})\·e^{-{\mathrm{ik}}_z\·\mathrm{\Δ}z}\\ =u(x,z;\mathrm{\ω})\·e^{-i{\[\frac{{\mathrm{\ω}}^2}{v_0^2(x,z)}+\frac{{\∂}^2}{\∂x^2}\]}^{\frac{1}{2}}\·\mathrm{\Δ}z}\·e^{-i(\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)})\·\mathrm{\Δ}z}\·e^{-i\frac{0.5601\×\[\frac{\mathrm{\ω}}{v(x,z)}-\frac{\mathrm{\ω}}{v_0(x,z)}\]+1.8364\×\frac{v(x,z)-v_0(x,z)}{\mathrm{\ω}}\·\frac{{\∂}^2}{\∂x^2}}{0.6316+0.7368\×\frac{v^2(x,z)+v_0^2(x,z)}{{\mathrm{\ω}}^2}\·\frac{{\∂}^2}{\∂x^2})}\·\mathrm{\Δ}z}\end{array}---\left(15\right)$

步骤5按照(5)式并以相关型成像原理对地震数据进行偏移延拓成像，获得最终的偏移成像叠加剖面。

2.根据权利要求1所述的一种基于优化系数的混合域的傅立叶有限差分叠前深度偏移方法，其特征在于：利用padé近似的有理函数逼近混合域中的波动方程偏移成像延拓算子，并利用切比雪夫多项式函数辅以优化padé近似的有理函数逼近展开式系数，从而推导出新的波动方程偏移成像延拓算子，使优化构造的波动方程偏移成像延拓算子与精确的单平方根算子的误差减小，从而提高复杂高陡构造区的地震偏移成像精度。