地震勘探所利用的数据主要是反射波。反射波场由背景速度扰动和高波数速度扰动共同构成^[1-2], 因此基于反射波的成像分为叠前深度偏移成像(确定高波数扰动)和反射波层析成像(确定背景速度扰动)两个步骤。偏移成像用于定位地下反射位置, 其实现了从数据空间到成像空间的映射; 层析成像用于反演传播路径上的速度场, 其实现了从成像空间到模型空间的映射。当然, 这两个步骤也可以耦合在一起进行数据域的迭代反演成像(如全波形反演)。相比于常规射线走时层析, 结合偏移成像在成像域进行波动方程线性化走时层析速度反演是当前比较实用有效且精度较高的技术组合。

基于傍轴近似的常规射线类成像方法是目前应用于地震数据层析及偏移的重要方法。传统的射线追踪方法一般局限于对射线路径及走时的描述, 其实现灵活高效、没有倾角限制且容易拓展到起伏地表。但是, 高频近似下的常规射线追踪认为中心射线代表着地震波的主能量, 在实现上仅仅利用中心射线来描述地震波传播, 这样的近似处理只能反映地震波的高频运动学特征。且在数值计算上对于复杂介质可能存在焦散及多到达时问题, 因此常规射线类成像方法(包括偏移和层析)在复杂构造情况下的应用效果并不十分理想。

高斯束传播算子是射线传播算子的发展, HILL^[3-4], HALE^[5]以及POPOV等^[6-7]将高斯束方法应用到偏移处理中并取得了较好的成像效果。Hill于1990年首先提出高斯束叠后深度偏移方法, 并于2001年将该方法推广到叠前深度偏移, 利用共偏移距道集中的对称性解决了高斯束叠前深度偏移中的执行效率问题, 非常适合海上拖缆地震数据的成像处理。对于陆上地震数据的成像处理, GRAY^[8]于2005年提出了快速精确的、适合于陆地起伏地表数据的炮域高斯束叠前深度偏移方法, 解决了高斯束偏移方法中的局部平面波分解因地表高程变化及地表速度横向变化带来的问题。自HILL奠定了高斯束偏移方法的理论基础以来, 衍生出了一系列实用化的束偏移技术(如控制束偏移和快速束偏移及自适应束偏移等)。高斯束叠前深度偏移技术的实现主要包括单个独立的高斯束传播的求解及所有高斯束叠加成像两个步骤。单个独立的高斯束传播分两步求得, 即通过运动学射线追踪求取中心射线的路径及走时, 通过动力学射线追踪获取中心射线附近的走时及高频能量分布。利用相互独立的高斯束描述波传播, 既保持了射线方法的高效性和灵活性, 又考虑了波场的动力学特征。高斯束偏移利用相互独立的高斯束叠加并成像, 解决了射线类方法中的多路径问题, 兼具了初至波到达时Kirchhoff积分偏移的灵活性和波动方程偏移的精确性, 是一种精致、精确且实现上灵活高效的深度偏移方法。

与成像方法类似, 在层析速度反演中也有介于常规射线理论和波动理论之间的方法, 如胖射线层析^[9]、菲涅尔体层析^[10-12]、高斯束层析^[13]、波路径层析^[14]等。胖射线层析从射线层析向波动层析走近了一步, 但其理论基础不够严格, 仅仅是对波动理论感性认识的结果。基于菲涅尔体的波形层析方法是近年来发展较快的重要方法。它利用射线周围的波动性质来同时拟合走时和波场的振幅, 其本质是有限频层析的一种。SEMTCHENOK等^[13]最早提出了高斯束层析方法, 利用高斯束展布范围内的波场分布建立层析方程, 对于一个炮检对通过建立多个层析方程来减小层析反演的病态性, 在一定程度上提高了反演稳定性。但是SEMTCHENOK在建立层析敏感度核函数(以下简称为核函数)时未给出详细的推导过程, 仅仅是根据直观认识利用高斯束的局部波场振幅作为加权系数进行求解计算, 其理论基础不严格。上面所述的层析速度反演方法大多在数据域实现, 在面对复杂构造低信噪比数据时，受叠前数据复杂波场及低信噪比的影响, 实际应用较为困难。邵荣峰等^[15]将高斯束层析方法应用到成像域, 并通过自动拾取等技术实现了偏移速度分析, 万弘等^[16]将构造约束等实用化方法加入到高斯束层析中以提高反演的稳健性, 但两者的理论基础仍然与文献[13]相似, 同样是简单采用了高斯束振幅作为加权系数构建层析核函数。李辉等^[17]从高斯束作为正演工具的角度重新推导了成像域高斯束层析核函数, 推导得到的层析核函数仍然是以高斯束振幅作为加权系数进行计算。为了更好地将偏移成像与层析反演结合起来, 本文将直接从偏移成像条件出发, 推导建立成像走时扰动与速度扰动的线性关系, 更直观地体现层析与偏移的关系, 给出新的成像域走时层析核函数。鉴于高斯束传播算子在叠前深度偏移领域中具有较强的实用性^[18-22], 且能够高效提取三维共方位-反射角度成像道集^[23], 本文从高斯束偏移角度道集出发, 在波动方程的一阶Born近似和Rytov近似下, 推导成像域反射波走时层析方程及其核函数, 从而发展一种新的成像域的波动方程线性化走时层析反演方法, 并利用高斯束传播算子计算该走时层析核函数(简称基于高斯束传播算子的成像域波动方程线性化走时层析方法为高斯束层析), 从而形成结合高斯束偏移及高斯束层析的迭代反演建模及成像方法和流程。

结合叠前深度偏移成像提取的角度域共成像点道集, 成像域的反射波层析成像可以将反投影路径分解成两个透射波分支, 类似于高斯束或单程波偏移成像的处理方式。其基本思想是认为成像点的走时扰动由炮点到成像点的走时扰动以及成像点到检波点的走时扰动之和决定。从而, 成像域反射波走时层析成像可以建立起包含两个核函数分支的走时扰动对速度扰动的线性关系。这两个核函数分别代表从炮点到成像点以及从成像点到检波点的透射波路径。考虑到陆上地震数据采样时的地表非水平、空间采样不规则和数据信噪比低, 为了更方便高效地提取方位-反射角成像道集, 将偏移成像与层析成像更好地配套起来, 本文选择高斯束算子作为波传播算子用于计算格林函数, 发展起高斯束偏移和高斯束层析配套方法。相对于常规射线层析, 由于高斯束层析核函数具有更严格的理论推导、更逼近波传播实际而反演精度更高, 同时由于构建的层析矩阵稀疏性较低而使得反演稳定性更强。基于高斯束传播算子的成像域走时层析方法与高斯束偏移相结合, 可形成具有典型特征波成像特点、适应复杂构造低信噪比数据的成像与建模工具。

1 高斯束叠前深度偏移及方位-反射角度道集提取方法

高斯束偏移从实现上讲, 主要包括运动学和动力学射线追踪及高斯束波场叠加两个部分。前者给出了单个独立高斯束的计算, 后者则表达了空间中任意一点的格林函数可由其邻域内所有高斯束的叠加得到, 如(1)式所示。该格林函数可用于实现高斯束叠前深度偏移, 亦可用于求解高斯束层析核函数。

$ G\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x'}};\omega } \right) = \frac{{{\rm{i}}\omega }}{{2{\rm{ \mathsf{ π} }}}}\int {\frac{{{\rm{d}}{p_x}{\rm{d}}{p_y}}}{{{p_z}}}{u_{{\rm{GB}}}}\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x'}},\mathit{\boldsymbol{p}};\omega } \right)} $

(1)

(1)式表示从震源点x =(x, y, z)出发, 在目标点x′=(x′, y′, z′)处的格林函数由每个相互独立的高斯束在该点的积分得到。其中, ω为圆频率; p为慢度矢量, 代表传播方向; u_GB表示从震源点x =(x, y, z)出发, 在目标点x′=(x′, y′, z′)接收的单个高斯束波场, 其表达为(2)式, 表示高斯束的构建通过运动学射线追踪确定其中心射线的轨迹和走时τ(s), 通过动力学射线追踪获取中心射线附近的动力学参数Q和M, 从而得到单个高斯束波场。该表达式在文献[3]中已经给出并讨论了其具体计算方法, 本文仅给出表达式, 不再赘述。

$ \begin{array}{l} {u_{{\rm{GB}}}}\left( {s,{q_1},{q_2}} \right) = \sqrt {\frac{{v\left( s \right)\det \mathit{\boldsymbol{Q}}\left( {{s_0}} \right)}}{{v\left( {{s_0}} \right)\det \mathit{\boldsymbol{Q}}\left( s \right)}}} \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left\{ {{\rm{i}}\omega \left[ {\tau \left( s \right) + \frac{1}{2}{\mathit{\boldsymbol{q}}^{\rm{T}}}\mathit{\boldsymbol{M}}\left( s \right)\mathit{\boldsymbol{q}}} \right]} \right\} \end{array} $

(2)

给出了高斯束计算格林函数的表达式之后, 通过频率域波场相关的成像条件即可实现高斯束叠前深度偏移。进一步, 根据高斯束叠前深度偏移成像条件, 可以推导成像域高斯束走时层析核函数, 并通过公式(1)利用高斯束积分表达格林函数实现该核函数的计算, 从而将高斯束偏移与高斯束层析有机结合起来, 实现复杂介质的偏移与层析迭代。

对于炮域实现的高斯束叠前深度偏移, 需要从炮点和接收点分别进行高斯束正演传播, 以求得从炮点出发的下行高斯束波场和从反射点出发到达检波器的上行高斯束波场。然后类似于单程波波动方程偏移, 使用上行波场和下行波场的互相关成像条件:

$ \begin{array}{l} {I_{{\rm{cs}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \frac{{ - 1}}{{2{\rm{ \mathsf{ π} }}}}\int {{\rm{d}}\omega \int {{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}\frac{{\partial {G^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{r}}},\omega } \right)}}{{\partial {z_{\rm{r}}}}}} } \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{G^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega } \right){D_{\rm{s}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega } \right) \end{array} $

(3)

(3)式给出的是单炮成像结果。其中, G(x, x_s, ω)与G(x, x_r, ω)分别表示从炮点x_s和从检波点x_r到成像点x的格林函数, “*”代表共轭; D_s(x_s, x_r, p_r, ω)表示基于炮道集的加高斯窗局部倾斜叠加。将(1)式代入(3)式, 进一步表示成:

$ \begin{array}{*{20}{c}} {{I_{{\rm{cs}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega } \right) = {\omega ^3}\int {{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}\int {{\rm{d}}{\mathit{\boldsymbol{p}}_{\rm{s}}}\int {{\rm{d}}{\mathit{\boldsymbol{p}}_{\rm{r}}}U_{{\rm{GB}}}^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega } \right)} } } \cdot }\\ {U_{{\rm{GB}}}^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{r}}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega } \right){D_{\rm{s}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega } \right)} \end{array} $

(4)

(4)式说明了任一单炮的成像结果是由从炮点出发的所有下行方向高斯束波场与从检波点出发的所有上行方向高斯束波场的相关得到。当速度存在误差导致成像存在误差时, 可以从(4)式出发推导建立速度误差与成像误差的线性关系, 给出基于高斯束偏移成像道集的成像域高斯束走时层析核函数, 将在下一节详细介绍。

在叠前深度偏移过程中提取方位-反射角度道集是进行后续成像域走时层析速度反演的必备过程。计算地下成像点的方位-反射角的一种有效且高效的方法是分别估算从炮点和检波点出发到达成像点处的波场传播方向。一旦获得两者的传播方向, 即可通过简单的向量代数运算得到方位角和反射角(张角)。对于射线类偏移方法而言, 这个过程相对容易, 只需通过旅行时场的空间导数分别计算震源波场的射线慢度p_s和检波点波场射线慢度p_r即可。对高斯束函数式((2)式)解析求解, 即可方便高效地计算旅行时场的空间导数。

假设高斯束邻域内任意一点Q, 其对应的中心射线上的点R, 根据(2)式所示的高斯束函数, 则点Q的旅行时可以由点R的旅行时及动力学射线追踪参量表示为:

$ \begin{array}{l} T\left( Q \right) = T\left( R \right) + \frac{1}{2}{\mathit{\boldsymbol{q}}^{\rm{T}}}{\mathop{\rm Re}\nolimits} \left( \mathit{\boldsymbol{M}} \right)\mathit{\boldsymbol{q}} = \\ \;\;\;\;\;\;\;\;\;\;\;T\left( R \right) + \frac{1}{2}\sum\limits_{i,j = 1}^2 {{q_i}{q_j}{M_{i,j}}} \end{array} $

(5)

得到高斯束的旅行时可进一步通过旅行时场的空间导数得到点Q的高斯束波场慢度矢量P(Q)=∂T(Q)/∂x, 因此可以得到慢度矢量的解析表达式:

$ \begin{array}{l} {P_k}\left( Q \right) = \frac{{\partial T\left( Q \right)}}{{\partial {x_k}}} = \frac{{\partial T\left( R \right)}}{{\partial {x_k}}} + \frac{1}{2}\sum\limits_{i,j = 1}^2 {\left( {{q_i}\frac{{\partial {q_j}}}{{\partial {x_k}}} + } \right.} \\ \left. {\;\;\;\;\;\;\;{q_j}\frac{{\partial {q_i}}}{{\partial {x_k}}}} \right){M_{i,j}} = {P_k}\left( R \right) + \frac{1}{2}\sum\limits_{i,j = 1}^2 {\left( {{q_i}\frac{{\partial {q_j}}}{{\partial {x_k}}} + } \right.} \\ \left. {\;\;\;\;\;\;\;{q_j}\frac{{\partial {q_i}}}{{\partial {x_k}}}} \right){M_{i,j}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;k = 1,2,3 \end{array} $

(6)

分别得到炮点波场和检波点波场的慢度矢量后, 根据反射角和方位角的定义, 即可给出具体求解方法, 详见参考文献[21], 此处不再赘述。

2 基于高斯束成像道集的成像域高斯束层析方法

层析反演方法可以统一表示成正问题的一阶近似:

$ \mathit{\boldsymbol{K}}_\mathit{\boldsymbol{m}}^\mathit{\boldsymbol{d}}\Delta \mathit{\boldsymbol{m = }}\Delta \mathit{\boldsymbol{d}} $

(7)

式中: K_m^d为核函数, 取决于正问题的理论基础, 包括基于射线理论、有限频理论、波动理论等方法; Δd为反演所采用的数据属性, 包括利用走时、相位、振幅或波形等属性; Δm为待反演的参数, 包括速度、密度、各向异性参数等。根据不同的反演需求和组合, 进而衍生出射线(走时)层析、有限频(走时)层析、菲涅尔体(走时、振幅、波形)层析、波动方程(走时、波形)层析、相位层析、速度层析、密度层析、各向异性参数层析等方法。

仍然基于层析方程式(7), 本文从高斯束偏移成像条件出发, 推导出基于高斯束积分计算格林函数的成像域波动方程线性化走时层析(简称为高斯束层析)核函数。从而实现基于高斯束偏移成像道集的成像域高斯束层析方法。

从公式(4)出发, 去掉所有的积分项, 仅考虑任意一个炮检对、任一传播方向的成像结果(所有炮检对、所有方向高斯束的成像结果叠加即为单炮成像结果), 将角度域高斯束偏移(GBM)成像条件重写为(8)式。为了后续表达易于区分, 将从震源出发的下行高斯束波场表示为S(x, p_s, ω; x_s, x_r), 将从检波点出发的上行高斯束波场表示为R(x, p_r, ω; x_s, x_r), 有:

$ \begin{array}{*{20}{c}} {{I_{{\rm{GBM}}}}\left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right) = S\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) \cdot }\\ {{R^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)} \end{array} $

(8)

(8)式是对应频率域任一炮检对某一传播方向p的高斯束成像结果(共方位-反射角成像道集)。其中, x =(x, y, z)表示成像点坐标; θ, φ分别表示成像点的反射张角及方位角。

在波动方程的一阶Born近似下, 波场U可以分解为背景波场U₀和扰动波场ΔU:

$ U = {U_0} + \Delta U $

(9)

因此从成像条件(8)式可以近似得到扰动像:

$ \begin{array}{l} {I_{{\rm{GBM}}}}\left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right) \approx \Delta S\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) \cdot \\ \;\;\;\;R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) + {S_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) \cdot \\ \;\;\;\;\Delta {R^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) \end{array} $

(10)

式中:ΔS, ΔR分别为一阶Born近似散射场; S₀, R₀分别为从炮点和检波点出发到成像点的背景波场。该成像条件表明, 成像点x处像的扰动来自炮点端和检波点端两个分支的影响。

进一步地, 根据文献[24]和文献[25], 一阶Born近似散射场ΔS与ΔR可以表示为:

$ \begin{array}{l} \Delta S\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) = \int\limits_{{V_{\rm{S}}}} {2k_0^2\frac{{\Delta v\left( \mathit{\boldsymbol{y}} \right)}}{{{v_0}\left( \mathit{\boldsymbol{y}} \right)}}{G_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;\mathit{\boldsymbol{y}}} \right)} \cdot \\ \;\;\;\;\;{S_0}\left( {\mathit{\boldsymbol{y}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right){\rm{d}}\mathit{\boldsymbol{y}} \end{array} $

(11a)

$ \begin{array}{l} \Delta {R^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) = \int\limits_{{V_{\rm{S}}}} {2k_0^2\frac{{\Delta v\left( \mathit{\boldsymbol{y}} \right)}}{{{v_0}\left( \mathit{\boldsymbol{y}} \right)}}{G_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},} \right.} \\ \;\;\;\;\left. {\omega ;\mathit{\boldsymbol{y}}} \right)R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right){\rm{d}}\mathit{\boldsymbol{y}} \end{array} $

(11b)

式中:V_S和V_R分别表示从炮点和从检波点到成像点的Born波路径; k₀=ω/v₀表示背景模型波数; Δv为待反演的速度扰动; G(x, p, ω; y)表示由高斯束积分计算的从点y到点x的格林函数; S₀(y, p_s, ω; x_s, x_r)和R₀(y, p_r, ω; x_s, x_r)分别表示在背景速度模型中从炮点和检波点传播到空间任意一点y的波场。

如果将(11)式代入(10)式, 形式上可以给出像的扰动(左端项)与速度扰动(右端项)的关系式。但实际操作时, 显然不能直接利用该式进行层析反演。对比数据域层析反演可以分析:数据域反演利用正演数据与实测数据的差在某种范数下(一般是二范数)最小作为误差泛函, 其实测数据是客观的, 可直接用来做逼近标准; 如果直接利用(10)式进行成像域反演, 则由于客观上无法得到真实像I_GBM(x, θ, φ, ω)从而无法得到扰动像ΔI_GBM(x, θ, φ, ω), 因此其本质是利用一个未知的中间量来估计另一未知量Δv。直接利用像的扰动这个概念来进行反演缺乏严格的判断标准。从另一个角度分析, 像的扰动是一个综合概念, 其实质包括了走时(位置)扰动和振幅扰动。实际计算时需要将像的扰动退化为走时扰动(深度扰动)或振幅扰动, 从而分别建立与速度扰动的关系式。由于像域的振幅扰动影响因素太复杂, 实际层析反演一般退化为仅利用走时扰动。

考虑退化到仅利用走时扰动进行层析反演, 进一步地, 将(10)式重写为:

$ \begin{array}{l} \Delta {I_{{\rm{GBM}}}}\left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right) \approx \frac{{\Delta S\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{{S_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}} \cdot \\ {S_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) + \\ \frac{{\Delta {R^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{S_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) \cdot \\ R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)\\ = {I_0}\left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right) \cdot \left[ {\frac{{\Delta S\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{{S_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}} + } \right.\\ \left. {\frac{{\Delta {R^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}} \right] \end{array} $

(12)

整理得:

$ \begin{array}{*{20}{c}} {\frac{{I\left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right)}}{{{I_0}\left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right)}} - 1 \approx \frac{{\Delta S\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{{S_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}} + }\\ {\frac{{\Delta {R^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}} \end{array} $

(13)

在波动方程的Rytov近似下, 波场可以表示为u=(A₀+ΔA)e^{i(φ₀+Δφ)}。成像值是两个波场相关得到, 因此同样可以表示I=(A₀+ΔA)e^{i(φ₀+Δφ)}, I₀(x, θ, φ, w)=A₀e^iφ₀, 由于ΔA≪A₀, 则(13)式可以进一步表示成:

$ \begin{array}{l} {{\rm{e}}^{{\rm{i}}\Delta \varphi }} - 1 \approx \frac{{\Delta S\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{{S_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\frac{{\Delta {R^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}} \end{array} $

(14)

由于Δφ趋于零, 则(14)式左端项近似等于iΔφ, 两边取虚部得:

$ \begin{array}{*{20}{c}} {\Delta \varphi \left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right) = {\mathop{\rm Im}\nolimits} \left[ {\frac{{\Delta S\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{{S_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}} \right] + }\\ {{\mathop{\rm Im}\nolimits} \left[ {\frac{{\Delta {R^ * }\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}{{R_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}} \right]} \end{array} $

(15)

进一步地, 根据文献[26], 单频相位扰动与单频走时扰动有近似关系:

$ \Delta t\left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right) = \frac{{\Delta \varphi \left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right)}}{\omega } $

(16)

同时将(11)式及(16)式代入(15)式整理, 最终得到成像域单频走时扰动与速度扰动的关系式:

$ \begin{array}{*{20}{c}} {\Delta {t_{{\rm{GBM}}}}\left( {\mathit{\boldsymbol{x}},\theta ,\varphi .\omega } \right) = \int\limits_{{V_{\rm{S}}}} {K_S^F\left( {\mathit{\boldsymbol{y}},\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)\Delta v\left( \mathit{\boldsymbol{y}} \right)} \cdot }\\ {{\rm{d}}\mathit{\boldsymbol{y}} + \int\limits_{{V_{\rm{R}}}} {K_R^F\left( {\mathit{\boldsymbol{y}},\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)\Delta v\left( \mathit{\boldsymbol{y}} \right){\rm{d}}\mathit{\boldsymbol{y}}} } \end{array} $

(17)

其中, K^F为频率域走时层析核函数, 其两个分支分别表示为:

$ \begin{array}{l} K_S^F\left( {\mathit{\boldsymbol{y}},\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) = \\ \;\;{\mathop{\rm Im}\nolimits} \left[ {\frac{{2\omega }}{{v_0^3\left( \mathit{\boldsymbol{y}} \right)}}\frac{{{G_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;\mathit{\boldsymbol{y}}} \right){G_0}\left( {\mathit{\boldsymbol{y}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{{G_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}} \right] \end{array} $

(18a)

$ \begin{array}{l} K_R^F\left( {\mathit{\boldsymbol{y}},\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) = \\ \;\;{\mathop{\rm Im}\nolimits} \left[ {\frac{{2\omega }}{{v_0^3\left( \mathit{\boldsymbol{y}} \right)}}\frac{{{G_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;\mathit{\boldsymbol{y}}} \right)G_0^ * \left( {\mathit{\boldsymbol{y}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{G_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}} \right] \end{array} $

(18b)

从(18)式可以看出, 成像域走时层析核函数的表现形式与LIU等^[12]给出的数据域菲涅尔体走时层析核函数的表现形式类似, 其本质都是有限频核函数, 关键是背景速度下格林函数的计算。

需要说明的是, 在背景模型中的波场可以表示成子波与格林函数的乘积, 因此在推导得到(18)式的过程中约去了子波项(假设子波不变), 仅剩下格林函数项。

由于实际操作时, 走时扰动Δt在时空域(角度域成像道集)测量得到, 与频率无关。因此, 定义带限地震信号的走时扰动可以用单频走时扰动加权叠加得到:

$ \Delta t = \int_{{\omega _0} - \Delta \omega }^{{\omega _0} + \Delta \omega } {W\left( \omega \right)\Delta t\left( \omega \right){\rm{d}}\omega } $

(19)

式中:W(ω)表示归一化的加权函数, 本文采用高斯函数$ \mathit{g}\left( \mathit{\omega } \right)=\frac{1}{\sqrt{2\pi \mathit{\delta }}}\text{exp}\left[-\frac{{{\left( \mathit{\omega }\text{-}{{\mathit{\omega }}_{\text{0}}} \right)}^{2}}}{2{{\mathit{\delta }}^{\text{2}}}} \right]$, 以ω₀为高斯函数的期望值(对称中心), Δω为高斯函数的标准差(半宽度), 则

$ W\left( \omega \right)\frac{{g\left( \omega \right)}}{{\int_{{\omega _0} - \Delta \omega }^{{\omega _0} + \Delta \omega } {g\left( \omega \right){\rm{d}}\omega } }} $

(20)

最终得到成像域带限走时扰动与速度扰动关系式:

$ \begin{array}{*{20}{c}} {\Delta {t_{{\rm{GBM}}}}\left( {\mathit{\boldsymbol{x}},\theta ,\varphi } \right) = \int\limits_{{V_{\rm{S}}}} {K_S^T\left( {\mathit{\boldsymbol{y}},\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}};{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)\Delta v\left( y \right){\rm{d}}y} + }\\ {\int\limits_{{V_{\rm{R}}}} {K_R^T\left( {\mathit{\boldsymbol{y}},\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}};{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)\Delta v\left( y \right){\rm{d}}\mathit{\boldsymbol{y}}} } \end{array} $

(21)

其中, K^T为带限走时层析核函数, 其两个分支分别表示为:

$ \begin{array}{l} K_S^T\left( {\mathit{\boldsymbol{y}},\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}};{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) = \int {{\rm{d}}\omega W\left( \omega \right)} \cdot \\ {\mathop{\rm Im}\nolimits} \left[ {\frac{{2\omega }}{{v_0^3\left( \mathit{\boldsymbol{y}} \right)}}\frac{{{G_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;\mathit{\boldsymbol{y}}} \right){G_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{{G_0}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{s}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}} \right] \end{array} $

(22a)

$ \begin{array}{l} K_R^T\left( {\mathit{\boldsymbol{y}},\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}};{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right) = \int {{\rm{d}}\omega W\left( \omega \right)} \cdot \\ {\mathop{\rm Im}\nolimits} \left[ {\frac{{2\omega }}{{v_0^3\left( \mathit{\boldsymbol{y}} \right)}}\frac{{{G_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;\mathit{\boldsymbol{y}}} \right)G_0^ * \left( {\mathit{\boldsymbol{y}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{G_0^ * \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{p}}_{\rm{r}}},\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}},{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right)}}} \right] \end{array} $

(22b)

至此, 基于高斯束角度道集、在波动方程的一阶Born近似和Rytov近似下导出了成像域带限走时层析方程(21)式及其对应的核函数表达式(22)式。从(22)式可以看出, 该核函数的本质是有限频核函数, 其求解主要是背景波场中格林函数的计算。利用(1)式的高斯束积分计算格林函数是一种精度较高且计算量较小的实用化计算方式。

需要说明的是, 尽管最终得到的有限频核函数在推导过程引入了多次线性化处理, 但有限频核函数相对于高频近似的射线层析核函数仍然更接近于波动传播的实际情况, 是对波动层析的一种折中近似。用此方法替换常规射线层析能提高精度, 且能借鉴常规射线层析流程实现技术的实用化。为进一步提高对复杂地质体(复杂构造、强横向变速等)的成像效果, 结合高斯束成像的优势, 有必要在后续研究中加强非线性项的合理近似和算法优化的方法研究。

下面通过数值试验分析单频和带限走时层析核函数的特征。设计背景模型参数:网格个数N_x=N_y=N_z=201, 网格间距dx=dy=dz=10 m, 速度为3 000 m/s, 震源和检波点坐标分别为(500, 1 000, 1 000)和(1 500, 1 000, 1 000), 带限走时层析核函数的频率范围取0~40 Hz, 单频核函数计算取中心频率20 Hz。图 1给出了对应的单频及带限走时层析核函数。从图 1g和图 1h上可以看出, 核函数在炮检连线的中心线上的敏感度为零, 这与常规射线层析的核函数仅分布在射线上的假设完全相反。刘玉柱^[10]详细分析了这种差异, 并指出了常规射线层析之所以仍然能取得较好效果的原因。

图 1反映的是单频和带限走时层析核函数。根据(22)式, 该核函数仅仅是成像域走时层析核函数的一个分支, 两个分支(炮点到成像点及成像点到检波点)则构成了成像域走时层析反演核函数的完整形态, 如图 2所示。用此带有一定宽度的高斯束层析核函数替换常规射线层析核函数能更准确逼近波实际传播方式, 提高层析反演精度和稳定性。

成像域高斯束层析反演的实际操作流程及具体计算方法可以完全借鉴常规射线层析反演的框架和算法, 包括成像剖面层位及成像道集RMO的自动拾取, 矩阵求解, 层位约束正则化方法等等。两者的区别仅仅是在核函数的表达与计算上, 因此该方法的实用性较强。

需要进一步说明的是, 成像域走时层析的具体实现与其走时误差的具体计算形式有关。XIE等^[27-29]给出了基于炮偏移后按照炮索引排列的成像道集的剩余时差(RMO)所对应的核函数的表现形式, 该核函数的两个分支并不对称, 这与成像道集上每一个RMO代表一炮的成像走时误差有关。本文所推导的核函数与XIE等给出的核函数表达形式一致, 但由于是针对高斯束偏移所提取的方位角度域成像道集, 所提取的RMO仅与方位角及反射角有关^[30], 核函数则表现为关于剖面上反射轴法方向对称的两个分支, 这一点并不会影响理论上的精度, 但会使得层析的实现更加自然方便。

3 理论模型及实际数据应用 3.1 理论模型数据

根据某地区实际地质构造设计的速度模型如图 3a所示, 断层发育, 地层高陡, 速度横向变化较大。横纵向采样点为731×550, 横向采样间隔为10 m, 纵向采样间隔为5 m。利用声波正演得到叠前炮记录(图 3b), 炮间距和道间距都是10 m, 中心放炮, 每炮361道。利用等梯度速度模型作为层析初始模型(图 4a)。对初始模型进行高斯束叠前深度偏移, 输出偏移剖面(图 4b)及角度道集(图 5), 可以看出初始模型对应的偏移剖面上的各层位成像深度并不准确, 绕射波也没有完全收敛。图 5对比了CDP400位置处提取的初始模型角度道集和层析后角度道集, 由于速度偏低, 初始角度道集同相轴上翘; 经过层析反演更新速度模型, 偏移后角度道集得到了拉平。图 6是高斯束层析经过5次迭代后得到的速度模型及其相应的高斯束偏移剖面。图 7是常规射线层析经过8次迭代得到的速度模型及其高斯束偏移剖面。与常规射线层析对比可知, 本文发展的高斯束层析方法能够提供更加丰富的速度信息, 迭代次数也小于常规射线层析。从抽取的单道速度曲线对比(图 8)可以看出, 高斯束层析反演的速度值分辨率更高, 更逼近真实速度。对比图 6b和图 7b的偏移深度可以看出, 高斯束层析偏移结果更加逼近真实深度, 反射界面归位到正确的深度位置, 绕射波收敛, 断层位置聚焦更好, 更新后偏移提取的角度道集更加接近真实速度模型下偏移提取的角度道集(图 5), 说明高斯束层析对复杂构造模型的速度反演精度优于常规射线层析方法。

3.2 实际地震资料

实际资料来自中国东北某复杂构造工区。该区基底之下是上古生界变质岩, 变质岩之上发育火山岩。目的层波场复杂, 存在多组断裂系统, 速度变化大, 深度域建模困难, 精确成像难度较大。处理要求纵向上能刻画火山机构特征和喷发期次, 横向上能分辨接触关系。合理突出与火山岩接触的地层反射。前期进行了时间域速度建模, 以此作为深度域初始模型进行后续高斯束层析反演建模。

本文的研究工作是在前期处理人员进行了预处理、叠前时间偏移及深度域初始建模及偏移成像的基础上开展的, 主要通过高斯束层析反演后的速度模型进行高斯束叠前深度偏移来体现本文方法的实际应用效果。从图 9所示的三维体模型上可以看出, 高斯束层析反演的速度模型分辨率明显高于常规射线层析反演结果。进一步对比过井速度剖面及偏移剖面可以更清楚地反映该效果。

从图 10所示的层析反演的速度模型上看, 高斯束层析结果体现了速度模型的更多细节信息, 进一步, 通过与井曲线的对比可以看出, 高斯束层析反演速度模型的精度更高(图 11)。

图 12为利用射线层析模型及高斯束层析模型对实际数据进行高斯束偏移的结果, 图 13为图 12中蓝框的放大显示。从图 12及图 13可以看出, 同样利用高斯束偏移方法, 高斯束层析反演的速度模型对应的偏移结果断裂成像质量更高, 断点更干脆, 同相轴更连续, 整体成像质量更高。从该实际资料处理可以认识到, 结合高斯束层析与高斯束偏移技术的处理方式, 可以更好地适应复杂构造资料的深度域速度建模, 提供更高质量的成像结果。

4 结论

本文从高斯束积分表达格林函数出发, 引出了高斯束偏移及提取方位反射角成像道集方法。基于角度域高斯束偏移成像条件, 在波动方程的一阶Born近似和Rytov近似下, 推导给出了成像域波动方程线性化走时层析成像方程及其核函数表达式, 并利用高斯束传播算子计算该核函数。利用高斯束层析核函数替代常规射线层析核函数(该核函数为常数1), 可以改进层析反演精度, 加快反演收敛, 从而形成了基于高斯束算子的偏移成像与层析成像的联合迭代速度建模与偏移处理方法, 相对于常规方法, 该方法具有更高的精度和更强的实用性。

本文发展的高斯束层析反演方法利用高斯束传播算子计算成像域走时层析核函数, 提供了一种新的成像域波动方程线性化近似走时层析反演方法。该方法主要利用反射数据的走时信息, 对初始模型的依赖性较低。但该方法仅能反演速度模型的低波数成分, 因而主要用于背景速度建模, 为偏移成像服务及为更高精度的反演方法(如FWI)提供初始模型。高斯束层析方法与高斯束偏移技术相结合, 可形成具有典型特征波成像特点的、适应低信噪比数据的成像与建模工具, 真正体现偏移成像与速度建模一体化的处理思想。