立体层析是将一个局部相干同相轴在炮道集与检波点道集内的射线参数水平分量(p_sx、p_rx, 以下简称地表观测p_x参数)纳入到数据空间之中, 使得数据空间的信息相对传统反射层析更为丰富。除了旅行时t外, 地表观测p_x参数与炮检点坐标也都纳入到立体层析数据空间中, 使得立体层析成为层析成像方法中唯一一种可以同时反演速度、反射点位置、反射张角与构造倾角的方法^[1-4]。

斜率层析属于立体层析的一种特殊情况, 在利用运动学反偏移重建数据空间的过程中, 会出现炮检点位置和旅行时信息已经匹配而地表p_x参数无法匹配的情况, 这时可利用p_x残差完成对速度模型、反射点坐标和散射角的更新^[5-6]。斜率层析算法也因此而得名。本文将上述情况下的各向同性介质斜率层析方法推广到VTI介质。与前人研究工作的不同之处在于, 当VTI介质中的运动学反偏移完成之后, 我们利用p_x残差仅反演三个各向异性参数, 而不反演反射点坐标和散射角这些与地下成像点有关的局部运动学信息。这是由于各向异性参数更新之后, 会立即执行VTI介质中的运动学射线追踪, 以找到更新后介质中的成像点位置(即运动学偏移), 运动学偏移完成后相当于更新了成像点的局部运动学信息。这样做的优点是缓解了算法实际应用中需要联合反演多种物理参数的压力。二维典型理论数据试验结果表明, 仅利用地表观测的p_x残差确实可以完成VTI介质中三个各向异性参数的估计, 将其与VTI介质中的运动学偏移/反偏移相结合, 是一种非常实用的各向异性参数估算策略。本文考虑到qP波依然是大部分实际地震采集中需要面对的主要波动现象, 基于文献[7]提出的VTI介质拟声波程函方程的汉密尔顿形式, 推导了qP波斜率层析所需的Fréchet偏导数公式, 给出了VTI介质中的运动学反偏移与VTI介质qP波斜率层析相结合的处理流程。需要强调的是, 本文所讨论的斜率层析只用到了p_x参数关于各向异性三参数的偏导数, 因此地表观测位置(s_x, r_x)和旅行时t关于地下各向异性介质的偏导数将被略去。

1 二维VTI介质qP波斜率层析偏导数推导

任何一种层析反演算法都必须建立数据空间残差与模型空间残差之间的线性关系, 即:

$ \Delta \mathit{\boldsymbol{d}} = \mathit{\boldsymbol{F}}\Delta \mathit{\boldsymbol{m}} $

(1)

基于方程(1), 利用观测到的数据残差就可以计算出模型空间的更新量, 达到更新初始模型的目的。其关键步骤是建立方程(1)所示的Fréchet偏导数矩阵(以下简称F矩阵)。注意, 在二维VTI介质qP波立体层析中, 数据空间为d=[s_x  r_x  t  p_sx  p_rx], 其中s_x、r_x为炮检点坐标; t为走时; p_sx、p_rx为炮检点处的射线参数水平分量。地下模型空间为m=[x₀  z₀  θ_s  θ_r  v_v  ε  δ], 其中x₀、z₀为地下反射点x的坐标; θ_s、θ_r分别代表从炮点一侧与检波点一侧出射的散射相角; v_v、ε、δ代表VTI介质中的三个各向异性参数。

$ {\mathit{\boldsymbol{F}}_{{\rm{stereo}}}} = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{\kappa _T}}}\frac{{\partial t}}{{\partial x}}}&{\frac{1}{{{\kappa _T}}}\frac{{\partial t}}{{\partial z}}}&{\frac{1}{{{\kappa _T}}}\frac{{\partial t}}{{\partial {\theta _s}}}}&{\frac{1}{{{\kappa _T}}}\frac{{\partial t}}{{\partial {\theta _r}}}}&{\frac{1}{{{\kappa _T}}}\frac{{\partial t}}{{\partial {v_{\rm{v}}}}}}&{\frac{1}{{{\kappa _T}}}\frac{{\partial t}}{{\partial \varepsilon }}}&{\frac{1}{{{\kappa _T}}}\frac{{\partial t}}{{\partial \delta }}}\\ {\frac{1}{{{\kappa _s}}}\frac{{\partial {s_x}}}{{\partial x}}}&{\frac{1}{{{\kappa _s}}}\frac{{\partial {s_x}}}{{\partial z}}}&{\frac{1}{{{\kappa _s}}}\frac{{\partial {s_x}}}{{\partial {\theta _s}}}}&{\frac{1}{{{\kappa _s}}}\frac{{\partial {s_x}}}{{\partial {\theta _r}}}}&{\frac{1}{{{\kappa _s}}}\frac{{\partial {s_x}}}{{\partial {\theta _{\rm{v}}}}}}&{\frac{1}{{{\kappa _s}}}\frac{{\partial {s_x}}}{{\partial \varepsilon }}}&{\frac{1}{{{\kappa _s}}}\frac{{\partial {s_x}}}{{\partial \delta }}}\\ {\frac{1}{{{\kappa _r}}}\frac{{\partial {r_x}}}{{\partial x}}}&{\frac{1}{{{\kappa _r}}}\frac{{\partial {r_x}}}{{\partial z}}}&{\frac{1}{{{\kappa _r}}}\frac{{\partial {r_x}}}{{\partial {\theta _s}}}}&{\frac{1}{{{\kappa _r}}}\frac{{\partial {r_x}}}{{\partial {\theta _r}}}}&{\frac{1}{{{\kappa _r}}}\frac{{\partial {r_x}}}{{\partial {\theta _{\rm{v}}}}}}&{\frac{1}{{{\kappa _r}}}\frac{{\partial {r_x}}}{{\partial \varepsilon }}}&{\frac{1}{{{\kappa _r}}}\frac{{\partial {r_x}}}{{\partial \delta }}}\\ {\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial x}}}&{\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial z}}}&{\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial {\theta _s}}}}&{\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial {\theta _r}}}}&{\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial {\theta _{\rm{v}}}}}}&{\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial \varepsilon }}}&{\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial \delta }}}\\ {\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial x}}}&{\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial z}}}&{\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial {\theta _s}}}}&{\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial {\theta _r}}}}&{\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial {v_{\rm{v}}}}}}&{\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial \varepsilon }}}&{\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial \delta }}} \end{array}} \right] $

(2)

方程(2)展示了二维VTI立体层析Fréchet偏导数矩阵, 其中第一行、第二行、第三行分别为旅行时t、炮点坐标s_x、检波点坐标r_x关于立体层析模型空间各分量的偏导数, 本文不予讨论。考虑到本文讨论的斜率层析仅限于反演各向异性参数, 这里仅关注炮点射线参数水平分量p_sx、检波点射线参数水平分量p_rx关于三个各向异性参数的偏导数, 关于散射相角和反射点位置的偏导数将不予讨论。因此, 本文使用的VTI介质qP波斜率层析Fréchet偏导数矩阵退化为:

$ {\mathit{\boldsymbol{F}}_{{\rm{slope}}}} = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial {v_{\rm{v}}}}}}&{\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial \varepsilon }}}&{\frac{1}{{{\kappa _{psx}}}}\frac{{\partial {p_{sx}}}}{{\partial \delta }}}\\ {\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial {v_{\rm{v}}}}}}&{\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial \varepsilon }}}&{\frac{1}{{{\kappa _{prx}}}}\frac{{\partial {p_{rx}}}}{{\partial \delta }}} \end{array}} \right] $

(3)

注意, 公式(2)、公式(3)中的κ为针对矩阵中每一行的数量级不同设置的均衡系数, 目的是保证每一行的数值范围基本一致。

WANG等^[8]基于重新参数化后的二维VTI介质拟声波程函方程和射线扰动理论^[9], 推导了二维qP波立体层析Fréchet偏导数矩阵。本文采用类似的方式推导出适用于二维qP波的斜率层析Fréchet偏导数矩阵。二维VTI介质拟声波程函方程如下^[7]:

$ v_{{\rm{nmo}}}^2\left( {1 + 2\eta } \right)p_x^2 + v_v^2p_z^2\left( {1 - 2v_{{\rm{nmo}}}^2\eta p_x^2} \right) = 1 $

(4)

式中, v_v是介质垂向速度; v_nmo是NMO速度。根据文献[7]可知各向异性参数之间有如下关系:η=(ε-δ)/(1+2δ)和v_nmo=v_v(1+2δ)^1/2。为了保证变量之间的独立性, 不妨采用一种等价的程函方程, 其形式如下:

$ v_0^2p_x^2 + v_0^2p_z^2 + 2\varepsilon v_0^2p_x^2\left( {1 - v_0^2p_z^2} \right) + 2\delta v_0^4p_x^2p_z^2 = 1 $

(5)

因此, 二维VTI介质下的qP波汉密尔顿算子可以写为:

$ \begin{array}{l} H = \frac{1}{2}\left( {v_{\rm{v}}^2p_x^2 + v_{\rm{v}}^2p_z^2 - 1} \right) + \varepsilon v_{\rm{v}}^2p_x^2\left( {1 - v_{\rm{v}}^2p_z^2} \right) + \\ \;\;\;\;\;\;\delta v_{\rm{v}}^4p_x^2p_z^2 = 0 \end{array} $

(6)

相应地, 与射线参数水平分量和垂直分量有关的射线追踪一阶微分方程组为:

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}{p_x}}}{{{\rm{d}}\sigma }} = - \frac{{\partial H}}{{\partial x}} = - \left[ {\left( {1 + 2\varepsilon } \right)p_x^2 + p_z^2 - 4\left( {\varepsilon - \delta } \right) \cdot } \right.}\\ {\left. {v_{\rm{v}}^2p_x^2p_z^2} \right]{v_{\rm{v}}}\frac{{\partial v}}{{\partial x}} - v_v^2p_x^2\left( {1 - v_{\rm{v}}^2p_z^2} \right)\frac{{\partial \varepsilon }}{{\partial x}} - v_{\rm{v}}^4p_x^2p_z^2\frac{{\partial \delta }}{{\partial x}}}\\ {\frac{{{\rm{d}}{p_z}}}{{{\rm{d}}\sigma }} = - \frac{{\partial H}}{{\partial z}} = - \left[ {\left( {1 + 2\varepsilon } \right)p_x^2 + p_z^2 - 4\left( {\varepsilon - \delta } \right) \cdot } \right.}\\ {\left. {v_{\rm{v}}^2p_x^2p_z^2} \right]{v_{\rm{v}}}\frac{{\partial v}}{{\partial z}} - v_{\rm{v}}^2p_x^2\left( {1 - v_{\rm{v}}^2p_z^2} \right)\frac{{\partial \varepsilon }}{{\partial z}} - v_{\rm{v}}^4p_x^2p_z^2\frac{{\partial \delta }}{{\partial z}}} \end{array}} \right. $

(7)

式中:σ为沿射线方向的滑动积分变量。

对方程(7)做全微分, 可得到由初始相空间扰动(Δp_x, Δp_z)及各向异性参数扰动(Δv_v, Δε, Δδ)所引起的射线轨迹上的相空间扰动:

$ \frac{{{\rm{d}}\left( {\Delta \mathit{\boldsymbol{\eta }}} \right)}}{{{\rm{d}}\sigma }} = \mathit{\boldsymbol{S}}\Delta \mathit{\boldsymbol{\eta }} + \Delta \mathit{\boldsymbol{w}} + \Delta \mathit{\boldsymbol{e}} + \Delta \mathit{\boldsymbol{d}} $

(8)

其中,

$ \begin{array}{l} \Delta \mathit{\boldsymbol{\eta }} = \left[ {\begin{array}{*{20}{c}} {\Delta {p_x}}\\ {\Delta {p_z}} \end{array}} \right]\\ \mathit{\boldsymbol{S}} = \left[ {\begin{array}{*{20}{c}} { - \frac{{{\partial ^2}H}}{{\partial {x^2}}}}&{ - \frac{{{\partial ^2}H}}{{\partial x\partial z}}}&{ - \frac{{{\partial ^2}H}}{{\partial x\partial {p_x}}}}&{ - \frac{{{\partial ^2}H}}{{\partial x\partial {p_z}}}}\\ { - \frac{{{\partial ^2}H}}{{\partial z\partial x}}}&{ - \frac{{{\partial ^2}H}}{{\partial {z^2}}}}&{ - \frac{{{\partial ^2}H}}{{\partial z\partial {p_x}}}}&{ - \frac{{{\partial ^2}H}}{{\partial z\partial {p_z}}}} \end{array}} \right] \end{array} $

(9a)

根据射线扰动理论^[9], 关于模型空间3个各向异性参数的扰动汉密尔顿形式可以写为:

$ \begin{array}{*{20}{c}} {\Delta H = \Delta {H_{{v_{\rm{v}}}}} + \Delta {H_\varepsilon } + \Delta {H_\delta } = \frac{{\partial H}}{{\partial {v_{\rm{v}}}}}\Delta {v_{\rm{v}}} + }\\ {\frac{{\partial H}}{{\partial \varepsilon }}\Delta \varepsilon + \frac{{\partial H}}{{\partial \delta }}\Delta \delta } \end{array} $

(9b)

(8) 式中后三项形式如下:

$ \begin{array}{l} \Delta \mathit{\boldsymbol{w}} = \left[ {\begin{array}{*{20}{c}} { - \frac{{{\partial ^2}H}}{{\partial x\partial {v_{\rm{v}}}}}\Delta {v_{\rm{v}}}}\\ { - \frac{{{\partial ^2}H}}{{\partial z\partial {v_{\rm{v}}}}}\Delta {v_{\rm{v}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \frac{{\partial \Delta {H_{{v_{\rm{v}}}}}}}{{\partial x}}}\\ { - \frac{{\partial \Delta {H_{{v_{\rm{v}}}}}}}{{\partial z}}} \end{array}} \right]\\ \Delta \mathit{\boldsymbol{e}} = \left[ {\begin{array}{*{20}{c}} { - \frac{{{\partial ^2}H}}{{\partial x\partial \varepsilon }}\Delta \varepsilon }\\ { - \frac{{{\partial ^2}H}}{{\partial z\partial \varepsilon }}\Delta \varepsilon } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \frac{{\partial \Delta {H_\varepsilon }}}{{\partial x}}}\\ { - \frac{{\partial \Delta {H_\varepsilon }}}{{\partial z}}} \end{array}} \right]\\ \Delta \mathit{\boldsymbol{d}} = \left[ {\begin{array}{*{20}{c}} { - \frac{{{\partial ^2}H}}{{\partial x\partial \delta }}\Delta \delta }\\ { - \frac{{{\partial ^2}H}}{{\partial z\partial \delta }}\Delta \delta } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \frac{{{\partial ^2}H}}{{\partial z\partial \delta }}\Delta \delta }\\ { - \frac{{\partial \Delta {H_\delta }}}{{\partial z}}} \end{array}} \right] \end{array} $

(9c)

结合(9a)、(9b)、(9c)式, 利用射线传播矩阵理论^[10]即可求得方程(8)的传播矩阵解:

$ \begin{array}{l} \Delta \mathit{\boldsymbol{\eta }}\left( {{\sigma _1}} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPi} }}\left( {{\sigma _1},{\sigma _0}} \right)\left[ {\Delta \mathit{\boldsymbol{\eta }}\left( {{\sigma _0}} \right) + \int_{\sigma 0}^{{\sigma _1}} {{\mathit{\boldsymbol{ \boldsymbol{\varPi} }}^{ - 1}}} \left( {{\sigma _1},{\sigma _0}} \right) \cdot } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {\left( {\Delta \mathit{\boldsymbol{w}} + \Delta \mathit{\boldsymbol{e}} + \Delta \mathit{\boldsymbol{d}}} \right){\rm{d}}\sigma } \right] \end{array} $

(10)

式中:σ₀代表射线出射位置, σ₁代表射线传播的最终位置, Π(σ₁, σ₀)表示从射线出射位置到最终位置的传播矩阵。由于本文仅利用p_x作为数据残差反演各向异性参数, 而不反演出射位置处的散射角信息, 因此(10)式可进一步简化为:

$ \begin{array}{l} \Delta \mathit{\boldsymbol{\eta }}\left( {{\sigma _1}} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPi} }}\left( {{\sigma _1},{\sigma _0}} \right)\left[ {\int_{\sigma 0}^{{\sigma _1}} {{\mathit{\boldsymbol{ \boldsymbol{\varPi} }}^{ - 1}}} \left( {{\sigma _1},{\sigma _0}} \right) \cdot } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {\left( {\Delta \mathit{\boldsymbol{w}} + \Delta \mathit{\boldsymbol{e}} + \Delta \mathit{\boldsymbol{d}}} \right){\rm{d}}\sigma } \right] \end{array} $

(11)

利用(11)式即可获得二维VTI各向异性介质qP波斜率层析所需的数据分量与模型空间各分量之间的一阶扰动关系, 也意味着斜率层析Fréchet偏导数矩阵F得到建立, 后续通过求解线性方程组(1), 就可以反演各向异性参数。

2 二维VTI介质qP波斜率层析中运动学反偏移及实现流程

首先介绍文献[5]和文献[6]提出的运动学反偏移概念以及由此导出的各向同性介质中的斜率层析方法。如图 1所示, v是地下成像点(x, z)处的速度, t_s、t_r分别是地下成像点(x, z)到炮点s、检波点r的单程走时, h是半偏移距, m是炮点s与检波点r的中点。图 1a表示地下模型信息与地表数据信息的几何关系; 图 1b为共偏移距偏移成像剖面, 在该剖面上拾取构造倾角ξ; 图 1c为偏移距域共成像点道集(CIG), 在该道集上拾取剩余曲率(RMO)的斜率信息tanφ。基于偏移后共偏移距成像剖面上搜索到的构造倾角ξ, 向地表发射线找到对应的炮检点坐标, 对应的射线水平参数分量设为p_sx和p_rx。

CHAURIS等^[5]首先发现了地表共炮点道集、共检波点道集内的同相轴斜率信息p_sx^*、p_rx^*(即射线出射到地表的水平射线参数信息)与成像域共成像点道集内的运动学信息(图 1b中的构造倾角ξ与图 1c中的剩余曲率的斜率tanφ)以及射线出射到地表时的射线水平参数分量p_sx、p_rx之间的定量关系:

$ \left( {p_{sx}^ * - {p_{sx}}} \right) + \left( {p_{rx}^ * - {p_{rx}}} \right) = 0 $

(12a)

$ \left( {p_{sx}^ * - {p_{sx}}} \right) + \left( {p_{rx}^ * - {p_{rx}}} \right) = \alpha \sin \varphi $

(12b)

式中, α=2vcosβcosξ, β是反射角, 它与两个散射角θ_s、θ_r以及构造倾角ξ的关系如下:

$ \cos {\theta _s} + \cos {\theta _r} = 2\cos \beta \cos \xi $

(12c)

由公式(12)可以发展一种特殊的层析反演方法:先分别计算炮点和检波点的旅行时表, 找到地下的等旅行时面, 再从检波点开始沿p_rx进行射线追踪, 找到地下反射点, 反追到炮点。根据公式(12), 由于背景速度是错误的, 通过射线正演获得的p_sx与实际观测到的p_sx^*之间必然有残差δp_sx。利用δp_sx以及传统的各向同性立体层析偏导数^[1]即可更新背景速度、反射点坐标和两个散射角。注意, 由于使用了检波点端的p_rx参数和旅行时的联合成像条件, 这时的数据残差仅有共检波点道集内观测到的δp_sx。这就是最早发展的各向同性介质中的斜率层析算法。CHAURIS等^[5]将这种基于共偏移距偏移成像剖面内错误的构造倾角ξ以及在共成像道集CIG上拾取的剩余曲率的斜率tanφ、向地表进行射线追踪并找到对应的炮检点位置的过程, 命名为运动学反偏移。杨锴等^[11]将结构张量算法引入运动学反偏移, 实现了更为精确高效的运动学信息提取; 同时进一步指出公式(12)的运动学反偏移虽然在理论上是正确的, 但是在实际应用中通过p_sx、p_rx换算真实的p_sx^*、p_rx^*时往往会出现一定的数值误差。熊凯等^[12]建议, 在射线追踪到地表找到正确的炮检点位置后, 直接利用结构张量实时搜索p_sx^*、p_rx^*。由于结构张量计算效率高, 只需很小的计算量即可获得稳健可靠的地表p_x参数信息, 因此这是一种更为实用的准确获取地表p_x参数信息的方式。

将上述各向同性介质中的运动学反偏移推广到VTI介质非常简单, 只需要将各向同性介质的运动学射线追踪更换为VTI介质的运动学射线追踪, 将各向同性克希霍夫叠前深度偏移更换为VTI介质克希霍夫叠前深度偏移即可。

图 2为VTI介质的运动学反偏移及斜率层析流程。简述如下:根据初始VTI介质三参数模型实施克希霍夫叠前深度偏移; 在小偏移距成像剖面上人工选取斜率层析所需的数据点位置, 在这些位置利用搜索到的剩余曲率信息追踪大偏移距剖面内的数据点位置, 基于这些位置即可利用结构张量进行搜索, 获得该位置处的构造倾角ξ与RMO曲线斜率信息^[12]; 基于构造倾角和偏移距信息, 即可向地表发射线找到对应的炮点和检波点, 获得p_sx、p_rx; 按照旅行时和炮检点坐标信息分别在共炮点道集和共检波点道集实施搜索, 可获得正确的慢度水平分量(即同相轴斜率)p_sx^*、p_rx^*。获得地表观测p_x参数的数据残差后, 利用公式(11)即可反演各向异性参数的更新量。在各向异性参数获得更新之后, 重新从地表向地下发射线并根据成像条件找到新的成像点位置(即运动学偏移), 在找到新的成像点位置的同时得到了新的地表p_x参数残差和RMO斜率信息。如果斜率残差或者RMO斜率信息达到精度门槛值, 则反演结束。如果斜率残差或者RMO斜率信息没有达到精度门槛值, 则重复上述过程。

3 理论数据测试

利用二维VTI多层背斜模型数据, 验证了本文计算公式的正确性和实现流程的可靠性。图 3为6层背斜VTI参数模型, 基于该模型进行VTI波动方程正演模拟。正演数据共601炮, 每炮801道, 炮间距20m, 道间距10m, 最大偏移距4000m, 记录时间为4s。图 4展示了正演数据的3个单炮记录。大量实际资料分析表明, δ参数通常是地震响应较弱的一个参数, 目前工业界主要通过利用井资料对各向同性深度偏移成像剖面进行标定后, 再利用真实深度与成像深度的差转换得到δ剖面, 这是一个比较成熟的技术^[13]。因此, 本文重点讨论如何利用提取出的地表p_x参数残差反演垂直速度v_v和ε, 将δ剖面设置为正确值, 不对其进行反演。

图 5展示了初始v_v剖面和初始ε剖面(常梯度)。为了检验本文算法联合反演这两个参数的能力, 这里使用的初始v_v各向同性模型为经过倾角时差校正(DMO)速度分析后转到深度域、并做了大尺度平滑后的速度模型。

基于上述两个剖面以及正确的δ剖面, 实施第1轮VTI介质克希霍夫深度偏移, 获得初始小偏移距叠加成像剖面(偏移距范围0~500m)和若干初始共成像点道集(图 6)。基于图 6a所示的初始成像剖面, 拾取了700多个数据点(图 6b), 然后从每个数据点出发, 利用结构张量算法从零偏移距到最大偏移距每隔5个偏移距向地表发出一对射线, 这样得到的数据空间一共有1.4×10⁴个射线对。图 6d展示了1500m偏移距的偏移成像剖面局部放大结果, 以及CDP1430处的CIG, 图中水平方向红线对应深度为820m。两条红线的交点为A, 利用结构张量算法可以快速地在1500m偏移距成像剖面内搜索出关于A点的构造倾角ξ, 同时在CDP1430处的CIG内搜索出RMO曲线的斜率信息tanφ。

基于上述拾取的数据点位置、构造倾角及道集剩余曲率信息, 在初始三参数的模型中进行VTI运动学射线追踪及二分法搜索^[11], 即可找到对A点成像有贡献的炮检点坐标。这里直接列出其追踪结果:炮点坐标为7100m, 检波点坐标为8600m。射线追踪到炮点坐标7100m时的p_sx参数为-4.5833×10^-4, 射线追踪到检波点坐标8600m时的p_rx参数为-1.8541×10^-4, 真实的p_sx^*为-4.0894×10^-4, 真实的p_rx^*为-1.6207×10^-4。图 6e、图 6f显示了共炮点记录和共检波点记录以及与这个炮检对对应的正确的地表p_x参数。需要强调的是, 这两个真实的p_sx^*、p_rx^*都是在地下发射线到地表找到正确的炮检点位置之后, 根据其旅行时与炮检点坐标信息实时搜索得到的, 因此其精度完全有保障。

获得地表观测p_x参数的数据残差后, 利用公式(11)即可求得各向异性参数的更新量。在各向异性参数获得更新之后, 重新从地表炮检点位置向地下发射线, 根据本文采用的等旅行时面共偏移距成像条件p_sx^*+p_rx^*=p_sx+p_rx, 可以搜索到新的成像点位置。在找到新的成像点位置的同时, 得到了新的地表p_x残差信息, 当p_x残差达到精度门槛值时, 反演结束。本次实验共迭代12轮, 误差泛函基本上不再下降。图 7a、图 7b分别展示了最终的v_v和ε参数反演结果, 图 7c展示了泛函下降曲线; 图 8对比了水平位置5km和10km处的v_v和ε(其中红线为真实模型、黑线为初始模型、蓝线为最终联合反演结果), 可以看到最终模型和真实模型的变化趋势非常接近。

图 9、图 10分别为基于初始模型和最终反演模型的CIG道集和VTI叠前深度偏移结果。图 9a对比了5km处基于初始模型和最终反演模型得到的CIG道集, 图 9b对比了10km处基于初始模型和最终反演模型得到的CIG道集, 可见基于最终反演模型获得的CIG均得到了有效拉平。图 10对比了基于初始模型和最终模型的PSDM偏移剖面, 可以看出中深部的成像品质大为提升, 说明本文方法的反演结果能够满足叠前深度偏移成像的要求。

4 结论和讨论

本文将各向同性介质中的斜率层析方法推广到二维VTI介质qP波的情况。基于射线扰动理论推导了地表p_x参数关于地下三个各向异性参数的偏导数, 建立了适合于二维VTI介质qP波斜率层析的线性方程组, 并将各向同性介质中的运动学偏移/反偏移策略推广到VTI介质中, 最终实现了二维VTI介质qP波斜率层析方法。研究结果表明, VTI介质中的运动学偏移/反偏移与二维VTI介质qP波斜率层析线性方程组相结合, 可以很好地实现二维VTI介质中的各向异性参数重建。该方法也可以方便地推广到三维或者TTI介质。

需要指出的是, 无论是各向同性介质还是各向异性介质中的斜率层析方法, 对初始模型都有一定要求:初始成像体的品质能够保证CIG上提取的RMO斜率tanφ和共偏移距成像剖面上提取的构造倾角ξ都是可靠的。如果这个要求无法满足, 那么后续通过运动学反偏移获取的地表p_x参数精度将难以保证。

其次, 在斜率层析中忽略旅行时对反射点位置和散射角的导数项不会降低算法的适用性。这是因为, 如果在立体层析反演中能够将数据残差通过合理的方式转移到某一个数据分量上, 则会使得数据空间得到合理的压缩, 这对提高反演的稳定性是有益的。这一点已经在前人发展的斜率层析和控制定向层析方法中得到了证实。同时, 本文方法没有将反射点位置和散射角纳入到模型空间内, 实际上又压缩了传统斜率层析的模型空间, 进一步提高了反演的稳定性, 有利于本文方法的推广应用。