CN101706583B - 一种多偏移距vsp成像的局域化相空间方法 - Google Patents

Abstract

本发明公开了一种多偏移距VSP成像的局域化相空间方法，采用基于Gauss-Daubechies紧标架(G-D紧标架)波场分解的局域化相空间单程波延拓算子对VSP上行波波场和震源下行波场进行延拓，并借助延拓算子的渐近展开解析式提高计算效率；为了能够有效的减弱成像剖面上的偏移假象，采用基于局部平面假设的成像条件：结合延拓波场的波数(对应传播方向)信息，在成像点处，基于镜面反射确定主入射角对应的主反射角，对局域化平面波进行相关成像并进行叠加。本方法考虑了局域化相空间延拓算子的渐近展开形式，计算效率较高；采用基于局部平面假设的成像条件，在不增加计算量的同时，也能够有效减弱偏移假象和干扰噪声。该成像方法有助于提高地震处理和解释的有效性和可靠性，可用于油气勘探中的地震信号处理，对井旁地下构造形态进行精细成像，消除假异常，提高油气田勘探开发的精度。

Description

一种多偏移距VSP成像的局域化相空间方法

技术领域：

本发明属于地震勘探领域，涉及一种地震勘探中的地震信号处理方法，尤其是一种多偏移距VSP成像的局域化相空间方法。

背景技术：

地震勘探的垂直地震剖面(VSP)观测系统中，反射波仅经过地表一次，接收到的地震记录分辨率较高，使得利用多偏移距VSP资料进行成像能够对地下结构进行精细刻画。

对多偏移距VSP资料上行波场进行叠前深度偏移成像常规的方法是采用和地面资料相类似的处理方法，例如射线类方法，基于单程波波动方程的偏移成像，波动方程逆时偏移等。

基于单程波波动方程的多偏移距VSP成像方法一般包括两个方面：一是对震源和接收波场进行深度延拓；二是利用成像条件对延拓得到的波场进行成像。VSP接收上行波和震源下行波波场的延拓可以采用相位移加插值以及局域化相空间(例如标架分解，局部余弦变换等)的波场递推方法等，即将波场分解到变换域(频域或相空间域)，得到分解系数，利用相应的变换域波场延拓算子，对波场的分解系数在深度方向进行递推。常规成像条件采用的是基于时间一致性的零延时相关成像条件，即在空间成像点处，将上行波场和下行波场的零延时相关系数作为成像值，包括时间域，频率域方法，它只考虑了上下行波的传播时间，没有考虑成像点处波场的方向性。

现有技术：(1)基于局域化相空间波场分解进行波场延拓：对深度z处的频率域VSP上行波场u_U(x，z，ω)和震源下行波场u_D(x，z，ω)在空间x方向做标架分解，得到标架分解系数，这里ω表示频率；然后利用相应的局域化相空间中的波场延拓算子将标架分解系数递推到深度z+Δz，将递推得到的深度z+Δz处的标架分解系数进行标架重构，就得到对应的上下行波波场分别为u_U(x，z+Δz，ω)，u_D(x，z+Δz，ω)，对频域波场u_U(x，z+Δz，ω)，u_D(x，z+Δz，ω)做反傅里叶变换可得到时间域波场u_U(x，z+Δz，t)，u_D(x，z+Δz，t)，这里t为时间；

(2)采用基于时间一致性的零延时相关成像条件进行成像，具体的频率域公式为：

$R(x,z+\mathrm{\Δz})={\∫u}_D^{*}(x,z+\mathrm{\Δz},\mathrm{\ω})u_U(x,z+\mathrm{\Δz},\mathrm{\ω})\mathrm{d\ω},$

式中，右上角的星号表示取复共轭；与频域公式等价的时间域成像条件表达式为：

$R(x,z+\mathrm{\Δz})={\∫u}_D^{*}(x,z+\mathrm{\Δz},t)u_U(x,z+\mathrm{\Δz},t)\mathrm{dt},$

式中，R(x，z+Δz)为空间点(x，z+Δz)处的成像值。

现有技术的缺点：

(1)基于局域化相空间的波场延拓算子具有空间和波数的双重局域化特性，能够较好地适应介质的横向非均匀，但是在计算机模拟中，由于采用的对偶标架没有解析的数学表达式，波场延拓算子的积分式采用数值离散求和的计算方法，运算量较大；且波场延拓算子没有明确的解析表达式，不易详细分析其局域化特性；

(2)基于时间一致性的零延时相关成像条件(包括时间域，频率域方法)，它只考虑了上下行波的传播时间，没有考虑成像点处波场的方向性；对于VSP资料，由于其接收孔径小，覆盖次数少，使得叠加后的成像剖面上仍然有较多的偏移绕射假象。

发明内容：

本发明针对常规局域化相空间波场延拓算子数值计算效率较低，没有解析表达式的问题，在满足渐近展开的条件下，采用延拓算子的渐近展开解析形式；针对常规单程波成像条件易引起成像剖面上偏移假象的问题，提供一种采用基于局部平面假设的成像条件进行多偏移距VSP成像的方法，这种成像条件考虑了成像点处波场传播的方向性。该成像方法能够对VSP进行准确成像，可以减少偏移剖面上偏移假象，可用于油气勘探中的地震信号处理，以查明地下构造形态，消除假异常，提高地震处理和油气解释的有效性和可靠性。

本发明的目的是通过以下技术方案来解决的，一种多偏移距VSP成像的局域化相空间方法：

步骤1：对VSP的上行波场u_U(x，z，t)和震源下行波场u_D(x，z，t)做时间t方向的傅里叶变换，得到频率域波场u_U(x，z，ω)和u_D(x，z，ω)；将频率ω离散成K个频率分量ω_k，k＝1，…，K，对VSP上行波场和震源下行波场的每一个频率分量u_U(x，z，ω_k)，u_D(x，z，ω_k)进行步骤2-5；

步骤2：对频域上下行波波场在空间x方向进行Gauss-Daubechies紧标架(以下简称G-D紧标架)分解；

步骤3：基于G-D紧标架的局域化相空间上下行波波场延拓算子

的计算和对应四维数组的计算机存储；

步骤4：采用步骤3中得到的局域化相空间上下行波波场延拓算子，将深度z+(j-1)Δz处的标架分解系数 $U_U({\stackrel{\‾}{x}}_{m_0},{\stackrel{\‾}{\mathrm{\ξ}}}_{n_0},z+(j-1)\mathrm{\Δz},{\mathrm{\ω}}_k),U_D({\stackrel{\‾}{x}}_{m_0},{\stackrel{\‾}{\mathrm{\ξ}}}_{n_0},z+(j-1)\mathrm{\Δz},{\mathrm{\ω}}_k)$ 采用下面的公式在深度方向进行延拓：

$U_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m_0}{\mathrm{\&Sigma;}}\underset{n_0}{\mathrm{\&Sigma;}}{P}_U<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>U_U({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z+(j-1)\mathrm{\&Delta;z},{\mathrm{\&omega;}}_k),$

$U_D({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m_0}{\mathrm{\&Sigma;}}\underset{n_0}{\mathrm{\&Sigma;}}{P}_D<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>U_D({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z+(j-1)\mathrm{\&Delta;z},{\mathrm{\&omega;}}_k),$

得到深度z+jΔz处的标架分解系数U_U(x_m，ξ_n，z+jΔz，ω_k)，U_D(x_m，ξ_n，z+jΔz，ω_k)；

步骤5：对步骤4得到的深度z+jΔz处的标架分解系数利用下面的公式进行空间重构：

$u_U^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m}{\mathrm{\&Sigma;}}{g}_{\mathrm{mn}}\left(x\right)U_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k);$

$u_D^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m}{\mathrm{\&Sigma;}}{g}_{\mathrm{mn}}\left(x\right)U_D({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k);$

式中，u_U ^a(x，ξ_n，z+jΔz，ω_k)，u_D ^a(x，ξ_n，z+jΔz，ω_k)表示重构得到的(x，z+jΔz)处的(上下行波)局域化平面波分量，标架函数的具体表达式为 $g\left(x\right)={\mathrm{\&pi;}}^{-\frac{1}{4}}{e}^{-\frac{x^2}{2}};$

步骤6：采用基于局部平面假设的相关成像条件对步骤5得到的局域化平面波分量进行成像，得到深度z+jΔz处的成像值R(x，z+jΔz)。

所述步骤2中对上行波波场u_U(x，z，ω_k)的G-D紧标架分解表达式为：

$U_U({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z,{\mathrm{\&omega;}}_k)=<u_U,{\tilde{g}}_{m_0,n_0}>={\&Integral;}_{-\&infin;}^{\&infin;}{u}_U(x,z,{\mathrm{\&omega;}}_k){\tilde{g}}_{m_0,n_0}^{*}\left(x\right)\mathrm{dx},$

式中，星号表示取复共轭， ${\tilde{g}}_{m_0,n_0}\left(x\right)=\frac{1}{\sqrt{s}}\tilde{g}\left(\frac{x-{\stackrel{\&OverBar;}{x}}_{m_0}}{s}\right)e^{i{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}x}$ 是对偶标架基

经过伸缩s，平移 ${\stackrel{\&OverBar;}{x}}_{m_0}=m_0\mathrm{\&chi;}$ 和调制 ${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}=n_0$ 得到的标架函数组，χ和

分别表示空间-波数域(相空间)平面上空间和波数方向的采样间隔，m₀，n₀为对应的采样点标号，取值范围分别为[-M/2，M/2-1]和[-N/2，N/2-1]，M，N为空间和波数采样点总数， $i=\sqrt{-1};$ 对偶标架的具体表达式为 $\tilde{g}\left(x\right)=\frac{1}{A}{\mathrm{\&pi;}}^{-\frac{1}{4}}{e}^{-\frac{x^2}{2}},$ 其中A为标架冗余度，一般取为4；同理，采用同样的计算公式可得震源下行波的G-D紧标架分解系数为

所述步骤3按照如下步骤：

(1)基于G-D紧标架的局域化相空间上行波波场延拓算子的具体表达式为：

$P_U<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0,}\mathrm{\&Delta;z}>=\frac{s}{A\sqrt{\mathrm{\&pi;}}}\exp (-\frac{s^2}{4}{({\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0})}^2)$

$\&times;\exp \left[i({\stackrel{\&OverBar;}{x}}_{m_0}{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}-{\stackrel{\&OverBar;}{x}}_m{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n)\right]I({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}),$

式中，x_m，为空间平移因子，ξ_n，为波数调制因子，

的具体表达式为：

$I({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z})={\&Integral;}_{-\&infin;}^{\&infin;}\exp [-s^2{(\mathrm{\&xi;}-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c)}^2-\mathrm{i\&xi;}({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\mathrm{i\&zeta;}\left(\mathrm{\&xi;}\right)\mathrm{\&Delta;z}]\mathrm{d\&xi;},$

式中， $\mathrm{\&zeta;}\left(\mathrm{\&xi;}\right)=\sqrt{\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}-{\mathrm{\&xi;}}^2},$ v为波传播速度，ξ表示波数变量，Re(ζ)≥0，Im(ζ)≤0， ${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c=\frac{{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n+{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}}{2};$ 将

表达式中被积函数的指数项表示为 $\mathrm{\&psi;}\left(\mathrm{\&xi;}\right)=-s^2{(\mathrm{\&xi;}-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c)}^2-\mathrm{i\&xi;}({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\mathrm{i\&zeta;}\left(\mathrm{\&xi;}\right)\mathrm{\&Delta;z},$ 其在ξ_c点对ξ的一阶，二阶导数分别为： ${\mathrm{\&psi;}}_1=\frac{\&PartialD;\mathrm{\&psi;}}{\&PartialD;\mathrm{\&xi;}}{|}_{\mathrm{\&xi;}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}=-i[({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\frac{{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c\mathrm{\&Delta;z}}{{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c}],$ ${\mathrm{\&psi;}}_2=\frac{{\&PartialD;}^2\mathrm{\&psi;}}{\&PartialD;{\mathrm{\&xi;}}^2}{|}_{\mathrm{\&xi;}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}=-2s^2+\mathrm{i\&Delta;z}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2{{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c}^3},$ 定义鞍点 ${\mathrm{\&xi;}}_s={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c-\frac{{\mathrm{\&psi;}}_1}{{\mathrm{\&psi;}}_2},$ 则渐近展开的条件可以表示为：|Im(ξ_s)|□ω_k/v和 $|{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c\&PlusMinus;\frac{\mathrm{\&omega;}}{v}|>\frac{1}{\sqrt{\left|{\mathrm{\&psi;}}_2\right|}};$ 当满足渐近展开的条件时，延拓算子采用其渐近展开式

具体表达式为：

${\stackrel{\&OverBar;}{P}}_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z})=\frac{\sqrt{2}s}{A\sqrt{\left|{\mathrm{\&psi;}}_2\right|}}\exp [-\frac{{({\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0})}^2}{2{\mathrm{\&Delta;}}_{\mathrm{\&xi;}}^2}-\frac{{({\stackrel{\&OverBar;}{x}}_m-{\stackrel{\&OverBar;}{x}}_{m_0}-\mathrm{\&Delta;}z\tan {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_c)}^2}{2{\mathrm{\&Delta;}}_x^2}]$

式中，各变量的表达式为： ${\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c=\sqrt{\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}-{{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}^2},$ ${\mathrm{\&Delta;}}_x=\sqrt{2}s\sqrt{[1+{\left(\frac{\mathrm{\&Delta;z}}{Z_R}\right)}^2},$ ${\mathrm{\&Delta;}}_{\mathrm{\&xi;}}=\sqrt{2}/s,\cos {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_c={\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_{c}v/{\mathrm{\&omega;}}_k,$ Z_R＝2s²ω_kcos³θ_c/v，α＝arg(ψ₂)。

从 $P_U<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>,{\stackrel{\&OverBar;}{P}}_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z})$ 的表达式可以看出，局域化相空间延拓算子在波数和空间方向均具有较好的局域化特性；鉴于此，在满足计算精度的前提下，可以根据m，n的取值限定m₀，n₀的取值范围：下标m的取值范围为[-M/2，M/2-1]，对应的m₀取值范围为[m-r_x，m+r_x]，r_x为空间采样的半径；下标n的取值范围为[-N/2，N/2-1]，对应n₀的取值范围为[n-r_ξ，n+r_ξ]，r_ξ为波数采样的半径；

(2)将

由

被积函数指数项指数部分-iζ(ξ)Δz改为iζ(ξ)Δz，即可得到震源下行波延拓算子

经过与(1)同样的推导，即可得下行波延拓算子的渐近展开条件以及对应的渐近展开解析式

将深度坐标离散化为z+jΔz，j＝1，2，…，J，重复步骤4和步骤5，步骤2得到的上下行波波场标架分解系数

为步骤4的波场延拓初值。

所述步骤6按照如下步骤：

根据波数与传播角度的对应关系，得到主入射角和反射角角度，分别记为 ${\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l=-{\sin }^{-1}\left(\frac{v}{{\mathrm{\&omega;}}_k}{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_l\right),$ ${\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_n={\sin }^{-1}\left(\frac{v}{{\mathrm{\&omega;}}_k}{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n\right);$ 这里角度选取为与延拓方向z轴夹角，逆时针旋转为正，图2所示主入射角θ_l＜0，对应的主反射角θ_n＞0；

利用基于局部平面假设的相关成像条件进行成像主要分三种情况：

(1)反射界面1是水平界面，波场的局域化平面波分量镜面反射时有θ_n＝-θ_l，相应的成像条件记为：

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}{u}_D^{a^{*}}(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

上式表示先对具有相同波数值ξ_l的上下行波波场分量计算相关，然后把所有的相关分量累加，作为点(x，z+jΔz)处的成像值；利用角度作为变量表示的成像条件记为：

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}{\cos }^2{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l{u}_D^{a*}(x,{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

它的含义是指对具有相同主反射角和入射角的波场分量进行相关成像；

(2)若(x，z+jΔz)处为倾斜地层，与水平方向夹角为θ₀，逆时针旋转为正，反射界面2，则成像条件应为

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}\cos {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l\cos ({\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l-2{\mathrm{\&theta;}}_0)u_D^{a*}(x,{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l+2{\mathrm{\&theta;}}_0,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

(3)另外，如果为弯曲界面，应取成像点处界面的切线方向作为局部地层倾角θ₀，然后利用情况(2)给出的公式进行成像。

本方法考虑了局域化相空间延拓算子的渐近展开形式，计算效率较高；采用基于局部平面假设的成像条件，在不增加计算量的同时，也能够有效减弱偏移假象和干扰噪声。该成像方法能够增强地震处理和油气解释的有效性和可靠性，提高石油勘探开发的精度。

附图说明：

图1为本发明的多偏移距VSP成像的局域化相空间方法流程图；

图2为本发明的地震波镜面反射示意图：(a)水平反射界面，(b)倾斜及弯曲界面。图中，θ_l，θ_n分别为入射波和反射波的主传播方向与z轴的夹角；(b)中的地层倾角θ₀定义为反射界面与水平方向x轴的夹角。

图3为本发明的模拟VSP成像结果图：(a)速度模型；(b)VSP上行波记录；(c)采用现有技术一得到的成像结果；(d)采用本发明提供的方法得到的成像结果。

图4多炮VSP成像叠加结果图：(a)采用现有技术一得到的成像结果；(b)采用本发明提供的方法得到的成像结果。

具体实施方式：

下面结合附图对本发明做进一步详细描述：

图1所示为本发明的流程图。本发明采用基于Gauss-Daubechies紧标架(以下简记为G-D紧标架)波场分解的局域化相空间单程波延拓算子对VSP上行波波场和震源下行波场进行延拓，并借助延拓算子的渐近展开解析式提高计算效率；为了能够有效的减弱成像剖面上的偏移假象，采用基于局部平面假设的成像条件：结合延拓波场的波数(对应传播方向)信息，在成像点处，基于镜面反射确定主入射角对应的主反射角，对局域化平面波进行相关成像并进行叠加。该成像方法能够对VSP进行准确成像，可以减少偏移剖面上偏移假象，可用于油气勘探中的地震信号处理，以查明井旁地下构造形态，消除假异常，提高油气勘探开发的精度。

如图1所示，本发明所提供的局域化相空间多偏移距VSP成像方法，具体通过如下步骤实施：

步骤1：对VSP的上行波场u_U(x，z，t)和震源下行波场u_D(x，z，t)做时间t方向的傅里叶变换，得到频率域波场u_U(x，z，ω)和u_D(x，z，ω)；

将频率ω离散成K个频率分量ω_k，k＝1，…，K，对VSP上行波场和震源下行波场的每一个频率分量u_U(x，z，ω_k)，u_D(x，z，ω_k)进行步骤2-5：

步骤2：对频域上下行波波场在空间x方向进行G-D紧标架分解。具体为，对上行波波场u_U(x，z，ω_k)的G-D紧标架分解表达式为：

$U_U({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z,{\mathrm{\&omega;}}_k)=<u_U,{\tilde{g}}_{m_0,n_0}>={\&Integral;}_{-\&infin;}^{\&infin;}{u}_U(x,z,{\mathrm{\&omega;}}_k){\stackrel{~*}{g}}_{m_0,n_0}\left(x\right)\mathrm{dx},$

式中，星号表示取复共轭， ${\tilde{g}}_{m_0,n_0}\left(x\right)=\frac{1}{\sqrt{s}}\tilde{g}\left(\frac{x-{\stackrel{\&OverBar;}{x}}_{m_0}}{s}\right)e^{i{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}x}$ 是对偶标架基经过伸缩s，平移 ${\stackrel{\&OverBar;}{x}}_{m_0}=m_0\mathrm{\&chi;}$ 和调制 ${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}=n_0$

得到的标架函数组，χ和分别表示空间-波数域(相空间)平面上空间和波数方向的采样间隔，m₀，n₀为对应的采样点标号，取值范围分别为[-M/2，M/2-1]和[-N/2，N/2-1]，M，N为空间和波数采样点总数，i为虚数单位；对偶标架的具体表达式为 $\tilde{g}\left(x\right)=\frac{1}{A}{\mathrm{\&pi;}}^{-\frac{1}{4}}{e}^{-\frac{x^2}{2}},$ 其中A为标架冗余度，一般可以取为4。

同理，采用同样的计算公式可得震源下行波的G-D紧标架分解系数为

步骤3：基于G-D紧标架的局域化相空间上下行波波场延拓算子

的计算和对应四维数组的计算机存储，四个变量为x_m，ξ_n，

其具体含义同步骤2，这里深度延拓步长Δz一般取固定值。

(1)基于G-D紧标架的局域化相空间上行波波场延拓算子的具体表达式为：

$P_U<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0,}\mathrm{\&Delta;z}>=\frac{s}{A\sqrt{\mathrm{\&pi;}}}\exp (-\frac{s^2}{4}{({\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0})}^2)$

$\&times;\exp \left[i({\stackrel{\&OverBar;}{x}}_{m_0}{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}-{\stackrel{\&OverBar;}{x}}_m{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n)\right]I({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}),$

式中，x_m，

为空间平移因子，ξ_n，

为波数调制因子，

的表达式为：

$I({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z})={\&Integral;}_{-\&infin;}^{\&infin;}\exp [-s^2{(\mathrm{\&xi;}-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c)}^2-\mathrm{i\&xi;}({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\mathrm{i\&zeta;}\left(\mathrm{\&xi;}\right)\mathrm{\&Delta;z}]\mathrm{d\&xi;},$

式中， $\mathrm{\&zeta;}\left(\mathrm{\&xi;}\right)=\sqrt{\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}-{\mathrm{\&xi;}}^2},$ v为波传播速度，ξ表示波数变量，Re(ζ)≥0，Im(ζ)≤0， ${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c=\frac{{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n+{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}}{2};$ 将

表达式中被积函数的指数项表示为 $\mathrm{\&psi;}\left(\mathrm{\&xi;}\right)=-s^2{(\mathrm{\&xi;}-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c)}^2-\mathrm{i\&xi;}({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\mathrm{i\&zeta;}\left(\mathrm{\&xi;}\right)\mathrm{\&Delta;z},$ 其在ξ_c点对ξ的一阶，二阶导数分别为： ${\mathrm{\&psi;}}_1=\frac{\&PartialD;\mathrm{\&psi;}}{\&PartialD;\mathrm{\&xi;}}{|}_{\mathrm{\&xi;}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}=-i[({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\frac{{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c\mathrm{\&Delta;z}}{{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c}],$ ${\mathrm{\&psi;}}_2=\frac{{\&PartialD;}^2\mathrm{\&psi;}}{\&PartialD;{\mathrm{\&xi;}}^2}{|}_{\mathrm{\&xi;}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}=-2s^2+\mathrm{i\&Delta;z}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2{{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c}^3},$ 定义鞍点 ${\mathrm{\&xi;}}_s={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c-\frac{{\mathrm{\&psi;}}_1}{{\mathrm{\&psi;}}_2},$ 则渐近展开的条件可以表示为：|Im(ξ_s)|□ω_k/v和 $|{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c\&PlusMinus;\frac{\mathrm{\&omega;}}{v}|>\frac{1}{\sqrt{\left|{\mathrm{\&psi;}}_2\right|}}.$ 当满足渐近展开的条件时，采用其渐近展开式具体表达式为：

${\stackrel{\&OverBar;}{P}}_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z})=\frac{\sqrt{2}s}{A\sqrt{\left|{\mathrm{\&psi;}}_2\right|}}\exp [-\frac{{({\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0})}^2}{2{\mathrm{\&Delta;}}_{\mathrm{\&xi;}}^2}-\frac{{({\stackrel{\&OverBar;}{x}}_m-{\stackrel{\&OverBar;}{x}}_{m_0}-\mathrm{\&Delta;}z\tan {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_c)}^2}{2{\mathrm{\&Delta;}}_x^2}]$

式中， ${\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c=\sqrt{\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}-{{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}^2},$ ${\mathrm{\&Delta;}}_x=\sqrt{2}s\sqrt{[1+{\left(\frac{\mathrm{\&Delta;z}}{Z_R}\right)}^2},$ ${\mathrm{\&Delta;}}_{\mathrm{\&xi;}}=\sqrt{2}/s,$ cosθ_c＝ζ_cv/ω_k，Z_R＝2s²ω_kcos³θ_c/v，α＝arg(ψ₂)。

从 $P_U<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>,{\stackrel{\&OverBar;}{P}}_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z})$ 的表达式可以看出，局域化相空间延拓算子在波数和空间方向均具有较好的局域化特性；鉴于此，在满足计算精度的前提下，可以根据m，n的取值限定m₀，n₀的取值范围：下标m的取值范围为[-M/2，M/2-1]，对应的m₀取值范围为[m-r_x，m+r_x]，r_x为空间采样的半径；下标n的取值范围为[-N/2，N/2-1]，对应n₀的取值范围为[n-r_ξ，n+r_ξ]，r_ξ为波数采样的半径。

(2)将中

被积函数指数项指数部分-iζ(ξ)Δz改为iζ(ξ)Δz，即可得到震源下行波延拓算子

经过与(1)同样的推导，即可得下行波延拓算子的渐近展开条件以及对应的渐近展开解析式

将深度坐标离散化为z+jΔz，j＝1，2，…，J，重复步骤4，5，步骤2得到的上下行波波场标架分解系数

为步骤4的波场延拓初值。

步骤4：采用步骤3中得到的局域化相空间上下行波波场延拓算子，将深度z+(j-1)Δz处的标架分解系数 $U_U({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z+(j-1)\mathrm{\&Delta;z},{\mathrm{\&omega;}}_k),U_D({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z+(j-1)\mathrm{\&Delta;z},{\mathrm{\&omega;}}_k)$ 采用下面的公式在深度方向进行延拓：

$U_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m_0}{\mathrm{\&Sigma;}}\underset{n_0}{\mathrm{\&Sigma;}}{P}_U<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>U_U({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z+(j-1)\mathrm{\&Delta;z},{\mathrm{\&omega;}}_k),$

$U_D({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m_0}{\mathrm{\&Sigma;}}\underset{n_0}{\mathrm{\&Sigma;}}{P}_D<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>U_D({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z+(j-1)\mathrm{\&Delta;z},{\mathrm{\&omega;}}_k),$

得到深度z+jΔz处的标架分解系数U_U(x_m，ξ_n，z+jΔz，ω_k)，U_D(x_m，ξ_n，z+jΔz，ω_k)。

步骤5：对步骤4得到的深度z+jΔz处的标架分解系数利用下面的公式进行空间重构：

$u_U^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m}{\mathrm{\&Sigma;}}{g}_{\mathrm{mn}}\left(x\right)U_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k);$

$u_D^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m}{\mathrm{\&Sigma;}}{g}_{\mathrm{mn}}\left(x\right)U_D({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k);$

式中，u_U ^a(x，ξ_n，z+jΔz，ω_k)，u_D ^a(x，ξ_n，z+jΔz，ω_k)表示重构得到的(x，z+jΔz)处的(上下行波)局域化平面波分量，标架函数的具体表达式为 $g\left(x\right)={\mathrm{\&pi;}}^{-\frac{1}{4}}{e}^{-\frac{x^2}{2}}.$

步骤6：采用基于局部平面假设的相关成像条件对步骤5得到的局域化平面波分量进行成像，得到深度z+jΔz处的成像值R(x，z+jΔz)。

参考图2的地震波平面反射示意图，根据波数与传播角度的对应关系，得到主入射角和反射角角度，分别记为 ${\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l=-{\sin }^{-1}\left(\frac{v}{{\mathrm{\&omega;}}_k}{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_l\right),{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_n={\sin }^{-1}\left(\frac{v}{{\mathrm{\&omega;}}_k}{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n\right);$ 这里角度选取为与延拓方向z轴夹角，逆时针旋转为正，则图2中示意的主入射角θ_l＜0，对应的主反射角θ_n＞0。利用基于局部平面假设的相关成像条件进行成像主要分三种情况：

(1)如图2(a)反射界面1所示的水平界面，波场的局域化平面波分量镜面反射时有θ_n＝-θ_l，相应的成像条件记为：

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}{u}_D^{a*}(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

上式表示先对具有相同波数值ξ_l的上下行波波场分量计算相关，然后把所有的相关分量累加，作为点(x，z+jΔz)处的成像值；利用角度作为变量表示的成像条件记为：

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}{\cos }^2{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l{u}_D^{a*}(x,{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

它的含义是指对具有相同主反射角和入射角(即主入射角和主反射角大小相等，方向相反)的波场分量进行相关成像。

(2)若(x，z+jΔz)处为倾斜地层，与水平方向夹角为θ₀(逆时针旋转为正)，如图2(b)所示的反 射界面2，则成像条件应为

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}\cos {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l\cos ({\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l-2{\mathrm{\&theta;}}_0)u_D^{a*}(x,{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l+2{\mathrm{\&theta;}}_0,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

(3)另外，如果为弯曲界面(如图2(b)所示的反射界面3)，应取成像点处界面的切线方向作为局部地层倾角θ₀，然后利用情况(2)给出的公式进行成像。对实际的VSP记录，可以结合该地区已知的地震资料，给出大致的地层倾角。

下面利用本发明提供的基于局域化相空间的VSP成像方法，对(1)模型VSP和(2)某油田实际多偏移距VSP资料进行成像，并和现有技术一的成像结果进行了对比讨论：

(1)如图3所示，图3(a)是速度模型，其中的黑竖线表示井中的检波器排列，共110个检波器，深度从135m到1225m，间隔为10m；选取主频50Hz的Ricker子波做为震源，并放置于地表，偏移距为0.7km。图3(b)为模拟得到的上行波记录，时间采样间隔为1ms。分别利用现有技术一和本发明的步骤1-6对震源下行波场和图3(b)所示的上行波场进行延拓，成像；图3(c)，3(d)是对应的成像结果，其中黑虚线表示反射界面，黑虚线的长短大致表示理论上的界面成像范围。对比3(c)和3(d)的成像剖面，可以看出，图3(c)的整个空间剖面(包括有效成像范围内)内有较多的偏移假象，例如A₀，B₀所指空间位置的绕射假象，从而较严重地影响了成像质量，而图3(d)上的反射界面较为清晰，与理论上的反射界面位置吻合较好，偏移绕射假象较弱。

(2)图4为实际多偏移距VSP的成像叠加结果，4(a)，4(b)分别为利用现有技术一和本发明的技术方案得到的成像结果；对比图4(a)，(b)可得，图4(b)的同相轴连续性较好，偏移假象和噪声干扰较图4(a)要弱。

另外，由于本发明技术方案的步骤3在满足渐近展开的条件下，采用延拓算子的渐近展开式，所以图3(d)，4(b)较对应的图3(c)，4(d)的计算效率要高。

综上所述，本发明提供的基于局域化相空间波场分解的多偏移距VSP成像方法，采用的延拓算子在空间和波数方向都有较好的局部化特性，并且在满足渐近展开的条件下，可以利用延拓算子的高频渐近展开式，从而提高计算效率；结合相空间中延拓波场的特点，采用了基于局部平面假设的相关成像条件，在不增加计算量的同时，也能够有效减弱偏移假象和噪声干扰。该方法主要用于油气勘探中的地震信号处理，以查明井旁地下构造形态，消除假异常。

最后需要说明的是，以上模型和实际资料算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围，在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。

Claims (1)

1.一种多偏移距VSP成像的局域化相空间方法，其特征在于：

步骤1：对VSP的上行波场u_U(x，z，t)和震源下行波场u_D(x，z，t)做时间t方向的傅里叶变换，得到频率域波场u_U(x，z，ω)和u_D(x，z，ω)，其中x表示空间x方向，z表示深度；将频率ω离散成K个频率分量ω_k，k＝1，…，K，对VSP上行波场和震源下行波场的每一个频率分量进行步骤2-5；

步骤2：对频域上下行波波场在空间x方向进行G-D紧标架分解；

所述步骤2中对上行波波场u_U(x，z，ω_k)的G-D紧标架分解表达式为：

$U_U({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z,{\mathrm{\&omega;}}_k)=<u_U,{\tilde{g}}_{m_0,n_0}>={\&Integral;}_{-\&infin;}^{\&infin;}{u}_U(x,z,{\mathrm{\&omega;}}_k){\stackrel{~*}{g}}_{m_0,n_0}\left(x\right)\mathrm{dx},$

式中，星号表示取复共轭，

是对偶标架基

经过伸缩s，平移

和调制得到的标架函数组，χ和Ξ分别表示空间-波数域平面上空间和波数方向的采样间隔，m₀，n₀为对应的采样点标号，取值范围分别为[-M/2，M/2-1]和[-N/2，N/2-1]，M，N为空间和波数采样点总数，i为虚数单位；对偶标架的具体表达式为

其中A为标架冗余度，取为4；同理，采用同样的计算公式可得震源下行波的G-D紧标架分解系数为

步骤3：基于G-D紧标架的局域化相空间上下行波波场延拓算子

的计算和对应四维数组的计算机存储，其中，

为空间平移因子，

为波数调制因子，Δz为深度延拓步长；

所述步骤3按照如下步骤：

(1)基于G-D紧标架的局域化相空间上行波波场延拓算子的具体表达式为：

$P_U<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>=\frac{s}{A\sqrt{\mathrm{\&pi;}}}\exp (-\frac{s^2}{4}{({\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0})}^2)$

$\&times;\exp \left[i({\stackrel{\&OverBar;}{x}}_{m_0}{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}-{\stackrel{\&OverBar;}{x}}_m{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n)\right]I({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}),$

式中，

为空间平移因子，为波数调制因子，s为伸缩s；

的表达式为：

$I({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z})={\&Integral;}_{-\&infin;}^{\&infin;}\exp [-s^2{(\mathrm{\&xi;}-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c)}^2-\mathrm{i\&xi;}({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\mathrm{i\&zeta;}\left(\mathrm{\&xi;}\right)\mathrm{\&Delta;z}]\mathrm{d\&xi;},$

式中，Re(ζ)≥0，Im(ζ)≤0，

v是速度；将

表达式中被积函数的指数项表示为 $\mathrm{\&psi;}\left(\mathrm{\&xi;}\right)=-s^2{(\mathrm{\&xi;}-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c)}^2-\mathrm{i\&xi;}({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\mathrm{i\&zeta;}\left(\mathrm{\&xi;}\right)\mathrm{\&Delta;z},$ 其在

点对ξ的一阶，二阶导数分别为： ${\mathrm{\&psi;}}_1=\frac{\&PartialD;\mathrm{\&psi;}}{\&PartialD;\mathrm{\&xi;}}{|}_{\mathrm{\&xi;}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}=-i[({\stackrel{\&OverBar;}{x}}_{m_0}-{\stackrel{\&OverBar;}{x}}_m)-\frac{{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c\mathrm{\&Delta;z}}{{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c}],$ ${\mathrm{\&psi;}}_2=\frac{{\&PartialD;}^2\mathrm{\&psi;}}{\&PartialD;{\mathrm{\&xi;}}^2}{|}_{\mathrm{\&xi;}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}=-{2s}^2+\mathrm{i\&Delta;z}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2{{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c}^3},$ 定义鞍点

则渐近展开的条件表示为：当满足渐近展开的条件时，采用其渐近展开式

具体表达式为：

${\stackrel{\&OverBar;}{P}}_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z})=\frac{\sqrt{2}s}{A\sqrt{\left|{\mathrm{\&psi;}}_2\right|}}\exp [-\frac{{({\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0})}^2}{{2\mathrm{\&Delta;}}_{\mathrm{\&xi;}}^2}-\frac{{({\stackrel{\&OverBar;}{x}}_m-{\stackrel{\&OverBar;}{x}}_{m_0}-\mathrm{\&Delta;}z\tan {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_c)}^2}{2{\mathrm{\&Delta;}}_x^2}]$

$\&times;\exp \left\{i\right[\frac{{({\stackrel{\&OverBar;}{x}}_m-{\stackrel{\&OverBar;}{x}}_{m_0}-\mathrm{\&Delta;}z\tan {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_c)}^2}{{2\mathrm{\&Delta;}}_x^2}-\frac{\mathrm{\&Delta;z}}{Z_R}+{\stackrel{\&OverBar;}{x}}_{m_0}{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0}-{\stackrel{\&OverBar;}{x}}_m{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n+{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c({\stackrel{\&OverBar;}{x}}_m-{\stackrel{\&OverBar;}{x}}_{m_0})-{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c\mathrm{\&Delta;z}+(\mathrm{\&pi;}-\mathrm{\&alpha;})/2\left]\right\},$

式中， ${\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_c=\sqrt{\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}-{{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_c}^2},$ ${\mathrm{\&Delta;}}_x=\sqrt{2}s+\sqrt{[1+{\left(\frac{\mathrm{\&Delta;z}}{Z_R}\right)}^2},$ ${\mathrm{\&Delta;}}_{\mathrm{\&xi;}}=\sqrt{2}/s,$ $\cos {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_c={\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_{c}v/{\mathrm{\&omega;}}_k,$ $Z_R={2s}^2{\mathrm{\&omega;}}_k{\cos }^3{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_c/v,$ α＝arg(ψ₂)；

从

的表达式可以看出，局域化相空间延拓算子在波数和空间方向均具有较好的局域化特性；鉴于此，在满足计算精度的前提下，能够根据m，n的取值限定m₀，n₀的取值范围：下标m的取值范围为[-M/2，M/2-1]，对应的m₀取值范围为[m-r_x，m+r_x]，r_x为空间采样的半径；下标n的取值范围为[-N/2，N/2-1]，对应n₀的取值范围为[n-r_ξ，n+r_ξ]，r_ξ为波数采样的半径；

(2)将中被积函数指数项指数部分-iζ(ξ)Δz改为iζ(ξ)Δz，即可得到震源下行波延拓算子

经过与(1)同样的推导，即可得下行波延拓算子的渐近展开条件以及对应的渐近展开解析式

将深度坐标离散化为z+jΔz，j＝1，2，…，J，重复步骤4和步骤5，步骤2得到的上下行波波场标架分解系数

为步骤4的波场延拓初值；

步骤4：采用步骤3中得到的局域化相空间上下行波波场延拓算子，将深度z+(j-1)Δz处的标架分解系数

采用下面的公式在深度方向进行延拓：

$U_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m_0}{\mathrm{\&Sigma;}}\underset{n_0}{\mathrm{\&Sigma;}}{P}_U{<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>U}_U({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z+(j-1)\mathrm{\&Delta;z},{\mathrm{\&omega;}}_k),$

$U_D({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m_0}{\mathrm{\&Sigma;}}\underset{n_0}{\mathrm{\&Sigma;}}{P}_D{<{\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n;{\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},\mathrm{\&Delta;z}>U}_D({\stackrel{\&OverBar;}{x}}_{m_0},{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_{n_0},z+(j-1)\mathrm{\&Delta;z},{\mathrm{\&omega;}}_k),$

得到深度z+jΔz处的标架分解系数

其中，Δz为深度延拓步长，j为正整数；

步骤5：对步骤4得到的深度z+jΔz处的标架分解系数利用下面的公式进行空间重构：

$u_U^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m}{\mathrm{\&Sigma;}}{g}_{\mathrm{mn}}\left(x\right)U_U({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k);$

$u_D^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)=\underset{m}{\mathrm{\&Sigma;}}{g}_{\mathrm{mn}}\left(x\right)U_D({\stackrel{\&OverBar;}{x}}_m,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k);$

式中，

表示重构得到的(x，z+jΔz)处的局域化平面波分量，标架函数的具体表达式为

其中，空间采样点标号为m，波数采样点标号为n，g_mn(x)为m、n处的标架函数；

步骤6：采用基于局部平面假设的相关成像条件对步骤5得到的局域化平面波分量进行成像，得到深度z+jΔz处的成像值R(x，z+jΔz)；

所述步骤6按照如下步骤：

根据波数与传播角度的对应关系，得到主入射角和反射角角度，分别记为

这里角度选取为与延拓方向z轴夹角，逆时针旋转为正，主入射角

对应的主反射角

利用基于局部平面假设的相关成像条件进行成像主要分三种情况：

(1)反射界面1是水平界面，波场的局域化平面波分量镜面反射时有

相应的成像条件记为：

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}{u}_D^{a*}(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

上式表示先对具有相同波数值的上下行波波场分量计算相关，然后把所有的相关分量累加，作为点(x，z+jΔz)处的成像值；利用角度作为变量表示的成像条件记为：

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}{\cos }^2{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l{u}_D^{a*}(x,{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

它的含义是指对具有相同主反射角和入射角的波场分量进行相关成像；

(2)若(x，z+jΔz)处为倾斜地层，与水平方向夹角为θ₀，逆时针旋转为正，则成像条件应为

$R(x,z+\mathrm{j\&Delta;z})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l}{\mathrm{\&Sigma;}}\frac{{{\mathrm{\&omega;}}_k}^2}{v^2}\cos {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l\cos ({\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l-{2\mathrm{\&theta;}}_0)u_D^{a*}(x,{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k)u_U^a(x,{-\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_l+2{\mathrm{\&theta;}}_0,z+\mathrm{j\&Delta;z},{\mathrm{\&omega;}}_k),$

(3)另外，如果为弯曲界面，应取成像点处界面的切线方向作为局部地层倾角θ₀，然后利用情况(2)给出的公式进行成像。

Images