目前, 我国油气地震勘探目标逐渐转向小尺度(薄互层)岩性油气藏, 油气地震勘探技术也逐渐向单点高密度、宽带、宽方位(“两宽一高”)地震勘探方向转变。常规的窄带、窄方位地震数据采集只能进行窄带反射系数的估计, 很多情况下只能实现构造的准确成像。“两宽一高”地震数据采集为精确估计背景(偏移)速度和宽带反射系数提供了数据基础。结合来自测井数据的密度建模和合理的深度域反演方法, 得到背景阻抗与宽带反射系数融合的宽带波阻抗成像结果, 是进行精确油气藏描述的重要手段^[1]。

反射波层析反演^[2]是主要的用于背景速度估计宏观速度建模技术。按照实现方式, 反射波层析主要分为成像域和数据域两类反演方法。典型的成像域反演方法如基于射线理论的成像道集层析(ray-based tomographic migration velocity analysis)^[3-5], 通过测量共(偏移距/角度)成像道集(common image gather, 简记为CIG)的剩余时差(residual moveout, RMO)估计反射波时差, 并用射线理论计算反射波路径从而实现速度反演。为考虑波传播的有限频带效应, 也可以用波路径^[6-8]代替射线路径, 如波路径层析^[9]。然而, 成像道集层析中的诸多假设和近似导致了反演的精度较低。高精度的速度反演通常需要额外的模型约束, 如测井约束或者构造约束等^[10-11]。另外, 基于波动方程的偏移速度分析(WEMVA)也是常用的成像域反演方法, 通过极大化成像道集的能量(如叠加能量最大化目标函数^[12])或惩罚成像道集中像的差异性或聚焦性等(如差异相似优化^[13-16], differential semblance optimization, 简记为DSO), 将“像差”向反射波路径反投影, 实现背景速度反演。该方法用波动方程作为正演算子, 克服了射线理论的高频近似假设, 通过采用波动方程产生成像道集, 可以实现自动化反演, 因此适用条件更广。然而, WEMVA方法的计算量巨大, 且存在地下照明不均衡以及泛函梯度假象等问题^[17-20], 影响了反演的收敛速度甚至导致反演不收敛。LUO等^[21]提出了全走时反演方法(full traveltime inversion, FTI), 构建了自动化的成像域走时反演方法^[22]。然而, 该方法仍然需要计算并存储成像道集, 对计算和存储都有较高的需求。

在数据域反演方法中, 叠前地震数据中一次反射/散射波的运动学信息(走时、斜率甚至曲率)的精确测量是一个关键问题。在传统的反射层析中^[23-25], 速度模型由光滑速度和界面构成, 并且需要在叠前数据中追踪并拾取来自目标反射界面的同相轴, 因而仅适用于简单层状模型。在斜率层析中^[26-27], 无需描述模型界面的几何形态及连续追踪同相轴, 因而模型描述和实现方式更加灵活。CHAURIS等^[28]指出, 数据域的斜率层析等价于极小化成像道集的剩余斜率。此外, 基于共反射面元叠加(CRS)得到的波前属性也可以用于速度反演^[29-30]。在立体层析(stereotomography)^[31-32]方法中, 基于扰动射线理论给出了数据(坐标、射线参数、双程走时)对模型(地下反射点位置, 局部反射界面的斜率, 反射射线张角)的Fréchet导数。近年来, 立体层析在模型表达、数据拾取、各向异性、反问题优化等方面取得了长足进步^[33-36]。然而, 受射线近似、多参数反演等影响因素制约, 仍需要引入模型正则化等模型先验信息来帮助反演收敛。

除了射线理论, 也有一些学者提出了基于波动方程传播算子的数据域反射波反演方法。如MA等^[37]用动态时间规整(dynamic time warping, DTW)算法估计两个信号的(时变)时差, 并构造了相应的反射走时反演方法。然而DTW目标函数的定义依赖于地震波振幅, 因而无法消除振幅对时差估计的影响。FENG等^[38]基于波动方程反偏移和互相关时差测量方式, 给出了反射波走时敏感度核函数表达, 但没有给出如何测量反射时差。李辉等^[39]基于地震数据的特征高斯波包分解, 实现了高斯波包框架下的反射走时反演方法。FENG等^[40]提出了基于地震信号稀疏表达方法的反射波层析反演方法, 但该方法仅适用于二维模型。此外, 为避免显式测量反射波时差, 可以利用观测数据和模拟数据互相关函数的某种加权范数直接构造目标泛函^[41]。但BAEK等^[42]指出, 此类泛函存在梯度符号问题, 仅对单个反射层成立, 因而无法处理反射地震数据。

上述分析表明, 对于反射地震信号, 如何能够充分利用一次波的运动学信息并高效地反演高精度的宏观速度模型仍然是一个颇具挑战性的问题。本文将在FENG等^[40]研究成果的基础上, 通过分析反射时差测量方法的局限性, 提出用地下偏移距聚焦代替反射波时差极小化准则, 实现其三维推广。首先, 引入特征波分解方法, 建立特征数据空间; 接着, 分析利用特征波成像准则测量反射时差的局限性, 并给出了一种新的反射波运动学残差测量标准, 即测量炮/检波场的地下偏移距; 随后, 给出了相应的误差泛函及其梯度计算方法。利用二维合成数据验证了本文方法的有效性。最后, 讨论了本文方法的优势和局限性。

1 数据域波动理论走时层析框架 1.1 数据域波动理论走时层析

贝叶斯反演理论是地震波反演成像的基础, 基于最大后验概率密度(MAP), 可以导出代价函数(J(m))^[43]:

$ \begin{array}{l} J(\mathit{\boldsymbol{m}}) = \frac{1}{2}{\rm{\{ }}{\left( {Ru(\mathit{\boldsymbol{m}}) - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right)^{\rm{H}}}\mathit{\boldsymbol{C}}_{\rm{D}}^{ - 1}\left( {Ru(\mathit{\boldsymbol{m}}) - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;{\left( {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}^{{\rm{ prior }}}}} \right)^{\rm{H}}}\mathit{\boldsymbol{C}}_{\rm{M}}^{ - 1}\left( {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}^{{\rm{ prior }}}}} \right){\rm{\} }} \end{array} $

(1)

式中:C_D为数据协方差矩阵; C_M为模型协方差矩阵; d^obs代表实测数据; m为模型参数; u(m)代表正问题; R为限制(采样)算子; m^prior代表参数的先验值。

然而, 基于(1)式进行参数反演仍然存在一些问题。如(1)式假定了正问题可以很好地预测观测数据, 并假定残差符合(广义)高斯分布。在实际应用中, 这些假设通常难以满足, 这也解释了为何全波形反演(FWI)难以在实际数据应用中反演出宽波数的参数场(弹性参数甚至Q值)。为了降低FWI对初始模型的依赖性, 并降低波形拟合的困难, FWI有很多变种, 可以抽象表达为:

$ \mathop {\min }\limits_{\mathit{\boldsymbol{m}} \in \mathit{\boldsymbol{M}}} J(\mathit{\boldsymbol{m}}) = {\left\langle {T\left( {Ru(\mathit{\boldsymbol{m}}), T\left( {{\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right)} \right.} \right\rangle _{\rm{L}}} + \alpha {\left\| {D\mathit{\boldsymbol{m}}} \right\|_p} $

(2)

式中:M为模型参数空间; 〈·〉_L代表数据拟合准则, 如L2范数、Wasserstein距离或者Kantorovich-Rubinstein范数等; 算子T为变换算子(如阻尼Laplace变换、积分变换、包络变换、时窗函数等); D为模型约束算子。然而, 上述波场变换算子无法考虑到实测波场, 因而难以用正问题准确预测。因此, 将正演波场作为系统的状态变量独立出来^[44], 理论上可以发展更有效的反演成像方法。

$ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{m}} \in \mathit{\boldsymbol{M}}, u \in \mathit{\boldsymbol{U}}} J(\mathit{\boldsymbol{m}}) = {\left\langle {T\left( {Ru(\mathit{\boldsymbol{m}}), T\left( {{\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right)} \right.} \right\rangle _{\rm{L}}} + \\ \;\;\;\;\;\alpha {\left\| {D\mathit{\boldsymbol{m}}} \right\|_p}\;\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;A(\mathit{\boldsymbol{m}})\mathit{\boldsymbol{u}} = f \end{array} $

(3)

式中:U为状态变量空间; A为正演算子; f为震源函数。

事实上, 仅将正演波场作为系统的状态变量仍然不够, 只有观测数据中可以被正演算子模拟的那部分信息才可以用于反演。王华忠等^[45]提出了特征波反演成像理论框架, 通过从实测数据中筛选出能被控制方程较准确预测的波场(特征波场^[46]), 构建凸性更好的反问题, 从而实现参数反演。显然, 相比于振幅, 地震波的运动学信息(如走时等)更容易预测, 因而利用走时构造的泛函稳定性更好。同时, 走时对于大尺度模型扰动具有更好的线性特征, 可用于反演大尺度的速度模型。因此, 需要提出高效且稳健的(绝对)走时测量及时差估计方法。

1.2 基于特征波场分解的反射波残差测量方法 1.2.1 地震数据的特征表达与特征数据空间

王华忠等^[46]给出了如下叠前地震数据的特征波场分解(characteristic wavefield decomposition, 简记为CWD)方法

$ \begin{array}{l} c\left( {t, {\mathit{\boldsymbol{p}}_{\rm{s}}}, {\mathit{\boldsymbol{p}}_{\rm{r}}};{\mathit{\boldsymbol{s}}_0}, {\mathit{\boldsymbol{r}}_0}} \right) = \mathop \smallint \limits_{{s_0} - {W_{\rm{s}}}}^{{s_0} + {W_{\rm{s}}}} \mathop \smallint \limits_{r0 - {W_{\rm{r}}}}^{{r_0} + {W_{\rm{r}}}} d\left( {t + {\mathit{\boldsymbol{p}}_{\rm{s}}}\left( {\mathit{\boldsymbol{s}} - {\mathit{\boldsymbol{s}}_0}} \right) + } \right.\\ \;\;\;\;\;\;\;\;{\mathit{\boldsymbol{p}}_{\rm{r}}}\left( {\mathit{\boldsymbol{r}} - {\mathit{\boldsymbol{r}}_0}} \right), \mathit{\boldsymbol{s}}, \mathit{\boldsymbol{r}}){\rm{d}}\mathit{\boldsymbol{s}}{\rm{d}}\mathit{\boldsymbol{r}} \end{array} $

(4)

式中:d(t, s, r)为地震记录; s和r分别为炮点和检波点的坐标; s₀和r₀分别为参考炮点和检波点的坐标; p_s和p_r分别为局部平面波的入射和出射射线参数的水平分量; W_s和W_r分别为炮点和检波点端局部平面波分解的空间窗长度(使得地震波满足局部线性特征); c(t, p_s, p_r; s₀, r₀)为特征波分解后得到的具有特定传播方向的时空局部化的子波。

考虑到(4)式需要预先估计局部平面波的射线参数才能实现特征波场分解, 冯波等^[47]在稀疏反演框架下将特征波场分解问题转化为如下模型参数估计问题(P₀)。

$ {P_0}:\mathop {\min }\limits_c {\left\| \mathit{\boldsymbol{c}} \right\|_0}\quad {\rm{ s}}{\rm{.t}}{\rm{. }}\quad \left\| {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Tc}}} \right\|_2^2 < {\sigma ^2} $

(5)

式中:‖c‖₀表示向量c中非零元素的个数; d=d(t, s, r; s₀, r₀); c=c(t, p_s, p_r; s₀, r₀); 算子T为(4)式中积分过程的伴随算子; σ²为噪声能量阈值。

经过特征波场分解((5)式), 可以直接将叠前地震数据投影到立体空间(炮/检点坐标, 入射/出射射线参数及走时)中。对于带限地震信号的到达时, 有多种定义方式, 如起跳时刻、最大振幅时刻等。本文采用立体层析方法中的走时拾取策略, 用地震道的包络极值定义地震波走时:

$ {T_{i, j, k}} = \mathop {{\mathop{\rm argmax}\nolimits} }\limits_t E\left\{ {\mathit{\boldsymbol{c}}\left( {t;{{\left( {{\mathit{\boldsymbol{p}}_{\rm{s}}}, {\mathit{\boldsymbol{p}}_{\rm{r}}}} \right)}_k}, {\mathit{\boldsymbol{s}}_i}, {\mathit{\boldsymbol{r}}_j}} \right)} \right\} $

(6)

式中:i, j分别是参考炮、检点的下标; k为特征波下标; E为包络算子。

显然, 通过特征波场分解((5)式)及走时定义((6)式), 可以自动生成如下数据空间:

$ \begin{array}{l} \mathit{\boldsymbol{D}} = \left\{ {\left( {T, {\mathit{\boldsymbol{p}}_s}, {\mathit{\boldsymbol{p}}_{\rm{r}}}} \right)_{k = 1}^{{n_p}};{s_i}, {r_j}} \right\}\\ \;\;\;i = 1, 2, \cdots , {n_s}\quad j = 1, 2, \cdots , {n_{\rm{r}}} \end{array} $

(7)

式中:n_s, n_r分别为炮点和检波点个数; n_p为局部同相轴的个数; T, p_s, p_r分别为局部同相轴的走时、入射和出射射线参数。由于(7)式中包含了地震数据中局部相干同相轴的运动学特征, 因而称为特征数据空间。此外, 本文给出的特征波场分解方法并不追求观测数据的精确重构, 仅仅提取了观测数据中符合高维局部平面波特征的波现象, 将其作为观测数据运动学信息的一种抽象与特征表达, 用于后续的背景速度反演。

1.2.2 反射波时差测量方法

为了测量反射波时差, 我们给出如下策略。如图 1所示, 可以在初始模型中从炮点和检波点处分别以观测数据的入射和出射射线参数向地下传播两条射线。若这两条射线相交于P点, 此时有如下走时恒等式:

$ \begin{array}{l} T_{i, j, k}^{{\rm{obs}}} = \Delta {t_{i, j, k}} + {T^{{\rm{cal}}}}\left( {{S_i}, P\left( {{x_i}, {z_i}} \right), p_{\rm{s}}^k;\mathit{\boldsymbol{m}}} \right) + \\ \;\;\;\;\;\;\;\;\;{T^{{\rm{ cal }}}}\left( {{R_j}, P\left( {{x_i}, {z_i}} \right), p_{\rm{r}}^k;\mathit{\boldsymbol{m}}} \right) \end{array} $

(8)

式中:m为速度模型参数化方式, m=1/v; T_{i, j, k}^obs为观测数据走时; T^cal(S_i, P, p^k_s; m)和T^cal(R_j, P, p^k_r; m)分别为初始模型m中从炮点S到交点P以及检波点R到交点P的两条射线的单程走时; Δt_{i, j, k}为反射时差; Q点为真实速度模型中反射点的位置。

此时, 反射时差可以直接写为:

$ \begin{array}{l} \Delta {t_{i, j, k}} = T_{i, j, k}^{{\rm{obs}}} - {T^{{\rm{cal}}}}\left( {{S_i}, P\left( {{x_i}, {z_i}} \right), p_{\rm{s}}^k;\mathit{\boldsymbol{m}}} \right) - \\ \;\;\;\;\;\;\;\;\;\;\;{T^{{\rm{cal}}}}\left( {{R_j}, P\left( {{x_i}, {z_i}} \right), p_{\rm{r}}^k;\mathit{\boldsymbol{m}}} \right) \end{array} $

(9)

显然, 利用本文给出的策略, 可以直接测量立体数据空间中每一个立体数据的反射时差, 既避免了周期跳跃以及串层等可能性, 又消除了振幅因素对时差估计的影响。

1.2.3 地下偏移距测量方法

然而, 在三维空间中, (9)式定义的时差测量方法却难以成立, 原因在于无法保证三维空间中两条射线相交。为了将反射走时反演方法推广至三维, 本文提出用地下某一深度处两条射线的横向距离代替时差测量, 克服三维空间中射线难以相交的局限性。如图 1所示, 若在某个深度z=z_e处, 两条单程射线的走时之和等于观测走时, 即:

$ \begin{array}{l} T_{i, j, k}^{{\rm{obs}}} = {T^{{\rm{cal}}}}\left( {{S_i}, M\left( {{x_{{\rm{se}}}}, {z_{\rm{e}}}} \right), p_{\rm{s}}^k;\mathit{\boldsymbol{m}}} \right) + \\ \;\;\;\;\;\;\;\;\;{T^{{\rm{ cal }}}}\left( {{R_j}, N\left( {{x_{{\rm{ re }}}}, {z_{{\rm{ e }}}}} \right), p_{{\rm{ r }}}^k;\mathit{\boldsymbol{m}}} \right) \end{array} $

(10)

式中:${T^{{\rm{ cal }}}}\left({{S_i}, M\left({{x_{{\rm{ se }}}}, {z_{\rm{e}}}} \right), p_{\rm{s}}^k; \mathit{\boldsymbol{m}}} \right)$和${T^{{\rm{cal}}}}\left({{R_j}, N\left({{x_{{\rm{re}}}}} \right.} \right., {\mathit{z}_{\rm{e}}}), p_{\rm{r}}^k; \mathit{\boldsymbol{m}})$分别为初始模型m中两条射线的单程走时; M(x_se, z_e)和N(x_re, z_e)为两条射线的终点坐标。

显然, 无论两条射线是否相交, 总会存在一个深度使得(10)式成立。定义地下偏移距为:

$ h = \frac{{{x_{{\rm{re}}}} - {x_{{\rm{se}}}}}}{2} $

(11)

若速度模型准确, 则地下偏移距趋近于0。因此, 可以利用地下偏移距构造目标泛函来实现速度反演。

1.3 反演方法

基于特征数据体, 可以直接计算反射时差或地下偏移距, 从而构造相应的目标泛函并实现速度反演。

1.3.1 反射时差目标泛函

利用(9)式的反射时差, 可以构造如下反射时差目标泛函:

$ J(\mathit{\boldsymbol{m}}) = \frac{1}{2}\sum\limits_i {\sum\limits_j {\sum\limits_k \Delta } } t_{i, j, k}^2 $

(12)

目标泛函的梯度为:

$ \nabla J(\mathit{\boldsymbol{m}}) = \sum\limits_i {\sum\limits_j {\sum\limits_k {\frac{{\partial \Delta {t_{i, j, k}}}}{{\partial \mathit{\boldsymbol{m}}}}} } } \Delta {t_{i, j, k}} $

(13)

根据公式(9), 反射时差对模型的Fréchet导数可以表示为:

$ \begin{array}{l} \frac{\partial }{{\partial \mathit{\boldsymbol{m}}}}\Delta {t_{i, j, k}}\;\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right. = - \frac{{\partial {T^{{\rm{cal}}}}}}{{\partial \mathit{\boldsymbol{m}}}}{\left. {\left( {{s_i}, q, p_{\rm{s}}^k;\mathit{\boldsymbol{m}}} \right)} \right|_{{\mathit{\boldsymbol{m}}_0}}} - \\ \;\;\;\;\;\frac{{\partial {T^{{\rm{ cal }}}}}}{{\partial \mathit{\boldsymbol{m}}}}{\left. {\left( {{r_j}, q, p_{\rm{r}}^k;{\mathit{\boldsymbol{m}}_0}} \right)} \right|_{{\mathit{\boldsymbol{m}}_0}}} \end{array} $

(14)

显然, 泛函梯度(13)式取决于走时对模型的Fréchet导数的具体形式, 其计算方法将在1.3.3节中讨论。

1.3.2 地下偏移距目标泛函

利用(11)式测量的地下偏移距, 可以构造如下泛函:

$ J(\mathit{\boldsymbol{m}}) = \frac{1}{2}\sum\limits_i {\sum\limits_j {\sum\limits_k {h_{i, j, k}^2} } } $

(15)

目标泛函的梯度为:

$ \begin{array}{l} \nabla J(\mathit{\boldsymbol{m}}) = \sum\limits_i {\sum\limits_j {\sum\limits_k {\frac{{\partial {h_{i, j, k}}}}{{\partial \mathit{\boldsymbol{m}}}}} } } {h_{i, j, k}} = \\ \;\;\;\;\;\;\;\;\;\sum\limits_i {\sum\limits_j {\sum\limits_k {\frac{1}{2}} } } \frac{{\partial \left( {{x_{{\rm{re}}}} - {x_{{\rm{see}}}}} \right)}}{{\partial \mathit{\boldsymbol{m}}}}{h_{i, j, k}} \end{array} $

(16)

为了计算射线终点横坐标对模型的Fréchet导数, 根据(10)式可得如下隐函数:

$ \begin{array}{l} F\left( {{z_{\rm{e}}}, {x_{{\rm{se}}}}, {x_{{\rm{re}}}}, \mathit{\boldsymbol{m}}} \right) = {T^{{\rm{cal}}}}\left( {{S_i}, M\left( {{x_{{\rm{se}}}}, {z_{\rm{e}}}} \right), p_{\rm{s}}^k;\mathit{\boldsymbol{m}}} \right) + \\ \;\;\;\;\;\;{T^{{\rm{ cal }}}}\left( {{R_j}, N\left( {{x_{{\rm{ re }}}}, {z_{\rm{e}}}} \right), p_r^k;\mathit{\boldsymbol{m}}} \right) - T_{i, j, k}^{{\rm{ obs }}} = 0 \end{array} $

(17)

利用隐函数求导法则, 有:

$ \frac{{\partial {x_{{\rm{re}}}}}}{{\partial \mathit{\boldsymbol{m}}}} = - \frac{{\partial F}}{{\partial \mathit{\boldsymbol{m}}}}/\frac{{\partial F}}{{\partial {x_{{\rm{re}}}}}}, \frac{{\partial {x_{{\rm{se}}}}}}{{\partial \mathit{\boldsymbol{m}}}} = - \frac{{\partial F}}{{\partial \mathit{\boldsymbol{m}}}}/\frac{{\partial F}}{{\partial {x_{{\rm{se}}}}}} $

(18)

其中,

$ \begin{array}{l} \frac{{\partial F}}{{\partial \mathit{\boldsymbol{m}}}} = \frac{\partial }{{\partial \mathit{\boldsymbol{m}}}}{T^{{\rm{cal}}}}\left( {{S_i}, M\left( {{x_{{\rm{se}}}}, {z_{\rm{e}}}} \right), p_s^k;\mathit{\boldsymbol{m}}} \right) + \\ \;\;\;\;\;\;\;\;\frac{\partial }{{\partial \mathit{\boldsymbol{m}}}}{T^{{\rm{cal}}}}\left( {{R_j}, N\left( {{x_{{\rm{re}}}}, {z_{\rm{e}}}} \right), p_{\rm{r}}^k;\mathit{\boldsymbol{m}}} \right) \end{array} $

(19)

$ \begin{array}{l} \frac{{\partial F}}{{\partial {x_{{\rm{se}}}}}} = \frac{\partial }{{\partial {x_{{\rm{se}}}}}}{T^{{\rm{cal}}}}\left( {{S_i}, M\left( {{x_{{\rm{se}}}}, {z_{\rm{e}}}} \right), p_{\rm{s}}^k;\mathit{\boldsymbol{m}}} \right) = \\ \;\;\;\;\;\;\;\;\;{p_x}\left( {{x_{{\rm{se}}}}, {z_{\rm{e}}}} \right) \end{array} $

(20)

$ \begin{array}{l} \frac{{\partial F}}{{\partial {x_{{\rm{re}}}}}} = \frac{\partial }{{\partial {x_{{\rm{re}}}}}}{T^{{\rm{cal}}}}\left( {{R_j}, N\left( {{x_{{\rm{re}}}}, {z_{\rm{e}}}} \right), p_{\rm{r}}^k;\mathit{\boldsymbol{m}}} \right) = \\ \;\;\;\;\;\;\;\;\;{p_x}\left( {{x_{{\rm{re}}}}, {z_{\rm{e}}}} \right) \end{array} $

(21)

泛函梯度((16)式)的核心仍然是走时对模型的Fréchet导数。接下来将讨论如何近似计算走时Fréchet导数。

1.3.3 梯度计算与模型更新

根据走时正问题的具体形式, 可以导出不同的$\partial T/\partial \mathit{\boldsymbol{m}}$的线性近似表达。如射线理论中, $\partial T\left({x, y} \right)/\partial \mathit{\boldsymbol{m}}$为从x到y的射线路径; 而在有限频层析理论中, $\partial T\left({x, y} \right)/\partial \mathit{\boldsymbol{m}}$为从x到y的有限频敏感度核函数^[48-49]; 对于波动方程走时反演, $\partial T\left({x, y} \right)/\partial \mathit{\boldsymbol{m}}$可以在一阶Born近似下用波场对模型的Fréchet导数表达^{[38, 50-51]}。虽然本文在高频近似意义下估计反射时差及地下偏移距, 但走时Fréchet导数的近似计算既可以采用射线近似, 也可以利用波动理论进行。考虑到本文方法定位于反演背景速度, 因而在数值实验中采用了射线近似方法来计算梯度, 计算效率远高于波动理论。

慢度模型可以用高斯-牛顿算法更新:

$ {\mathit{\boldsymbol{m}}^{k + 1}} = {\mathit{\boldsymbol{m}}^k} - H{\left( {{\mathit{\boldsymbol{m}}^k}} \right)^{ - 1}}P\left( {g\left( {{\mathit{\boldsymbol{m}}^k}} \right)} \right) $

(22)

式中:P为梯度预条件算子, 用于消除梯度中的高波数噪声^[39]; $g\left({{\mathit{\boldsymbol{m}}^k}} \right) = \nabla J{\left. {(\mathit{\boldsymbol{m}})} \right|_{\mathit{\boldsymbol{m}} = \mathit{\boldsymbol{m}}k}}$为走时或者地下偏移距泛函梯度; $H\left({{\mathit{\boldsymbol{m}}^k}} \right) = {\left. {\left({\frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{T}}}}}{{\partial \mathit{\boldsymbol{m}}}}\frac{{\partial \mathit{\boldsymbol{r}}}}{{\partial \mathit{\boldsymbol{m}}}}} \right)} \right|_{\mathit{\boldsymbol{m}} = \mathit{\boldsymbol{m}}k}}$; r为反射时差或地下偏移距。

2 数值试验

数值试验分为两部分。首先用偏大及偏小的初始模型测试CRI泛函梯度方向, 并分析梯度性态。然后用二维合成数据测试反演方法。

2.1 梯度方向测试

前文中提到DSO泛函梯度存在垂向条带状噪声^{[17, 20]}, 因而影响DSO泛函的收敛速度甚至导致其不收敛。本文采用文献[20]中的测试模型(包含数个截断的水平反射界面, 如图 2a所示)来测试CRI泛函梯度。速度模型参数为:X和Z方向采样点数分别为401和101, 网格间距均为5m。地表炮点和检波点的位置相同, 从250m开始至1750m结束, 间隔10m。因此地震数据共151炮, 每炮151道。考虑到测试模型为二维速度场, 我们采用反射走时目标泛函并用射线理论近似计算泛函梯度, 其优点在于计算效率极高(仅需要做2×N_c次初值射线追踪即可, N_c为特征波数据的个数)。梯度计算采用了10m×10m的网格。分别用扰动-10%(1700m/s)和+10%(2300m/s)的初始速度模型计算泛函梯度, 结果如图 2b和图 2c所示。由于模型离散化导致的稀疏射线路径使得CRI梯度中存在一些高波数噪声, 因此, 采用梯度滤波方法^[39]获得梯度的光滑部分, 如图 2d和图 2e所示。无论初始速度偏高或者偏低, CRI梯度均能获得正确的速度更新方向。结合梯度滤波方法, 可以更好地保留梯度中的低波数成分, 因而有利于反演收敛。另一方面, 虽然采用有限频理论或者波动理论计算的CRI梯度更加平滑, 但计算量增加了一个甚至几个数量级。考虑到CRI反演方法定位于提供背景速度模型, 因此在后续的反演测试中采用了射线理论近似计算梯度, 并结合梯度去噪方法兼顾反演收敛速度和计算效率。

2.2 二维合成数据测试

理论分析可知, 地震波走时对背景速度扰动有更好的线性特征。我们采用二维合成数据验证CRI方法。速度模型取自Sigsbee速度模型的一部分(x方向0~6700m, z方向2300~6100m), 如图 3a所示。速度模型以沉积岩为主, 并有数组断层发育。x和z方向的采样点数分别为801和601, 网格间距为7.5m。单边观测系统, 最小偏移距300m, 最大偏移距3900m, 道间距为15m。第一炮在x=0处激发, 炮间距30m, 共201炮。震源子波采用20Hz主频的Ricker子波。记录时间5s, 采样间隔1ms。显然, 这样的观测系统接收到的地震数据主要为反射波。

初始模型采用1500m/s的常速模型(与真实模型的地表速度相同), 采用反射波走时反演方法进行速度反演。反演终止准则为走时残差的L2范数小于初始走时残差的L2范数的1%。同时, 考虑到计算效率, 泛函梯度用射线近似计算, 并在迭代过程中动态调整梯度的平滑半径。

图 3b为CRI走时反演得到的速度模型, 与真实模型的光滑背景部分吻合程度较好。为进一步对比反演结果, 我们抽取不同地表位置(x=2000, 4000, 6000m)的纵向速度, 如图 4所示。图中红色直线、蓝色曲线和黑色曲线分别为初始速度、反演得到的速度和真实速度。显然, CRI反演的速度与真实速度的变化趋势较为吻合, 较好地恢复了中-大尺度的背景速度结构。为对比成像结果和道集, 我们进行了高斯束偏移(GBM)并输出角度域共成像点道集(ADCIG)。图 5为利用初始模型、反演模型和真实模型计算得到的GBM角度道集。可以看出, 经过反射走时反演之后, 角度道集基本拉平。图 6a, 图 6b和图 6c分别表示利用初始模型、反演模型和真实模型计算得到的GBM偏移剖面。对比速度模型、成像结果和角度道集可以看出, 由于深度速度横向变化较为剧烈(深度大于3500m)且缺乏足够的照明角度, 导致反演结果欠佳。同时, 最深部的水平反射界面的成像深度及形态与真实模型存在一定偏差。这种小尺度的速度结构变化需要用更高精度的反演方法(如全波形反演等)进一步提高速度模型的分辨率。

3 讨论与结论

基于特征波场分解和包络走时定义准则, 可以提取叠前地震数据的运动学信息, 从而构造特征数据体。在此基础之上, 通过分析特征波在地下的聚焦特性, 可以自动测量反射波残差(时差或者地下偏移距), 避免了地震波振幅对反射波运动学信息残差的影响。利用极小化反射波时差或地下偏移距两种准则均可实现背景速度反演。理论分析可知:反射波时差测量方法仅适用于二维模型, 而地下偏移距测量方法可以适用于二维和三维模型, 因而适用性更广。数值实验表明:本文方法无需长偏移距观测数据或低频信息, 对初始模型依赖性低、计算效率高, 且整个反演流程可以实现全自动化, 可以为后续的宽波数带速度建模提供较为可靠的低波数背景速度模型。

对于三维观测系统, 本文中采用的特征波分解方法要求炮点和检波点在主测线(In-line)和联络测线(Cross-line)的采样满足地震数据的高维局部平面波假设。因此, 若联络测线方向上炮点或检波点采样间隔过大, 会影响射线参数的反演精度, 进而影响速度反演的可信度。所以, 后续的研究需要在地震数据稀疏分解的同时, 给出数据可信度(或精度)的估计, 作为数据加权矩阵融入目标泛函, 提高反问题的稳健性。

在计算效率方面, 本文在数值实验中采用射线近似计算泛函梯度, 结合梯度去噪方法实现背景速度更新。虽然提升了计算效率, 但牺牲了反演精度, 因而只能反演低波数的速度结构。下一步将研究如何利用波动理论快速计算泛函梯度, 并结合构造特征等先验信息的约束尽可能恢复中波数的速度成分。此外, 为达到宽带波阻抗成像的要求, 背景各向异性参数及Q值估计也是下一步的研究方向。

致谢: 感谢中石油勘探开发研究院及西北分院、中海油研究院和湛江分公司、中国石油化工股份有限公司石油物探技术研究院和胜利油田分公司对波现象与智能反演成像研究组(WPI)研究工作的资助与支持。