逆时偏移(RTM)技术相对其它偏移方法具有成像精度高的优势^[1-3], 能够处理速度横向变化剧烈、陡倾角构造等复杂地质情况^[4]。近年来, 随着计算机技术的快速发展, RTM方法逐渐走向实用化^[5-7]。基于弹性介质假设的逆时偏移方法在理论上考虑了实际介质中波传播多波多分量问题, 能得到较为准确的成像结果, 因而备受关注。

弹性波逆时偏移问题的研究主要包括震源、检波点波场重建和成像条件应用。不少学者在弹性波RTM过程中先进行纵横波波场分离然后再成像^[8-10]。采用的主要方法是基于Helmholtz分解的波场分离^[11], 可分为直接波场分离法^[12-14]和间接波场分离法。直接波场分离法是在空间域对波场进行散度和旋度运算分离得到纵、横波波场, 但存在振幅相位校正问题。AKI等^[15]选择空间域矢量分解方法解决了成像时极性反转问题, 但其波场的物理意义发生了改变; ZHU^[16]尝试在波数域进行波场分离, 但该方法需要进行数据域转换以及波数域归一化处理等才可保证其物理意义。间接波场分离技术则是改造波动方程, 在波场延拓过程中分解得到纵横波分量。马德堂等^[17]利用各向同性介质纵横波分离方程实现波场完全分解。XIAO等^[18]通过纵波方程得到纵波波场再与弹性波波场作差得到分离的横波波场, 此类方法能够得到与原波场一致的纵横波波场。WANG等^[19]将基于纵横波矢量解耦的波动方程应用于弹性波逆时偏移成像, 得到独立的纵横波偏移结果。吴潇等^[20]对比了5种波场分离方法在弹性波逆时偏移中的应用效果。与直接波场分离法相比, 间接波场分离法成像质量更好, 但这类方法需要改造波动方程, 增加了计算量和内存存储。

应用准确成像条件是弹性波逆时偏移成像的关键环节。激发时间成像条件是利用提取震源与检波点波场对应成像时间波场值, 得到最终偏移结果^[21-22]; 还有通过提取震源、检波点波场作比值进行成像^[23-24]。该类方法计算效率高, 且有一定的保幅特性, 但转换波成像物理意义有待研究。较为常用的互相关成像条件能够利用全波场信息, 有学者通过选择性成像^[25]和滤波的方法压制背向散射噪声干扰^[26], 但该压制技术同时也会对有效信号造成损伤。ROCHA等^[26-27]利用弹性波波动方程及能量守恒原则, 提出新的成像条件, 成像结果表示反射系数能量, 无需对整个波场进行纵横波分离, 也避免了极性反转等传统弹性波偏移成像的问题。张晓语等^[28-29]在ROCHA等提出方法的基础上发展了基于能量密度自解耦成像条件, 此成像条件能够有效压制低频散射成像干扰, 主要应用了二阶方程, 对于具有较高效率的交错网格一阶方程尚无具体实现方式。本文发展了一阶速度应力弹性波能量范数成像条件的逆时偏移方法, 通过构造中间应变量给出了一阶方程情况下各分量转换规则, 提出了具体波场更新方式。最后通过数值算例测试了该方法对弹性波逆时偏移中低频噪声的压制效果。

1 方法原理 1.1 弹性波逆时偏移理论

逆时偏移理论框架可按其实现步骤来描述:①震源波场正向外推; ②检波点波场逆时延拓; ③成像条件应用。震源波场正向外推可简化为如下边值问题^[1]:

$ \left\{ \begin{array}{l} {U_{\rm{s}}}\left( {x,z,t} \right) = 0\;\;\;\;\;\;\;\;\;\;t \le 0\\ {U_{\rm{s}}}\left( {x,z,t} \right) = f\left( t \right)\;\;\;\;\;t > 0 \end{array} \right. $

(1)

式中:U_s(x, z, t)为震源波场值; f(t)为加载震源函数。波场值更新采用数值离散化的波动方程求取。检波点波场逆时延拓则与正向外推相反^[1], 可表达为:

$ \left\{ {\begin{array}{*{20}{l}} {{U_{\rm{r}}}\left( {x,z,t} \right) = 0}&{t > {T_{\max }}}\\ {{U_{\rm{r}}}\left( {x,z,t} \right) = g\left( t \right)}&{t \le {T_{\max }}} \end{array}} \right. $

(2)

式中:g(t)为地震记录值; U_r(x, z, t)为检波点波场值; T_max为最大记录时间。公式(2)中将检波点地震记录作为初始条件, 利用波场延拓算子求取进行逆时延拓。逆时偏移是记录每一计算递推时刻的震源波场和检波点波场, 然后选取相应成像条件得到最终偏移剖面。

本文方法基于完全弹性假设, 波场延拓过程采用各向同性介质一阶速度应力弹性波方程^[3]:

$ \left\{ \begin{array}{l} \frac{{\partial {\sigma _{xx}}}}{{\partial t}} = \rho v_{\rm{P}}^2\left( {\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}} \right) - 2\rho v_{\rm{S}}^2\frac{{\partial {v_z}}}{{\partial z}}\\ \frac{{\partial {\sigma _{zz}}}}{{\partial t}} = \rho v_{\rm{P}}^2\left( {\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_z}}}{{\partial z}}} \right) - 2\rho v_{\rm{S}}^2\frac{{\partial {v_x}}}{{\partial x}}\\ \frac{{\partial {\sigma _{xz}}}}{{\partial t}} = 2\rho v_{\rm{S}}^2\frac{1}{2}\left( {\frac{{\partial {v_z}}}{{\partial x}} + \frac{{\partial {v_x}}}{{\partial z}}} \right)\\ \frac{{\partial {v_x}}}{{\partial t}} = \frac{1}{\rho }\left( {\frac{{\partial {\sigma _{xx}}}}{{\partial x}} + \frac{{\partial {\sigma _{xz}}}}{{\partial z}}} \right)\\ \frac{{\partial {v_z}}}{{\partial t}} = \frac{1}{\rho }\left( {\frac{{\partial {\sigma _{xz}}}}{{\partial x}} + \frac{{\partial {\sigma _{zz}}}}{{\partial z}}} \right) \end{array} \right. $

(3)

式中:σ_xx、σ_zz、σ_xz为应力分量; v_x、v_z为速度分量; v_P、v_S分别为介质中纵、横波速度; ρ为密度。本文采用伪谱法对此方程进行数值离散。

1.2 能量范数成像条件

逆时偏移中常用的互相关成像条件^[3]表达为:

$ I = \sum\limits_{e,t} {{U_{\rm{s}}}} {U_{\rm{r}}} $

(4)

式中:I为成像结果; U_s为震源波场; U_r为检波点波场; e表示炮索引号; t为时间。在应用该成像条件时, 存在检波点波场与震源波场相同路径传播的互相关噪声。

能量范数成像条件首先从能量角度以一维方程声波定义一个空间域(Ω)内能量表达^[26]:

$ E\left( t \right) = \frac{1}{2}\int_\mathit{\Omega } {\left[ {{{\left( {\frac{{\partial U}}{{\partial t}}} \right)}^2} + {v^2}{{\left| {\nabla U} \right|}^2}} \right]{\rm{d}}x} $

(5)

式中:U为波场值; v为介质中波速。定义波场能量范数为:

$ \left\| U \right\|_E^2 = \sum\limits_{x,t} {\left[ {\frac{{{\partial ^2}U}}{{\partial {t^2}}} + {v^2}{{\left| {\nabla U} \right|}^2}} \right]} $

(6)

类似两个波场U、V的内积能量范数可表达为:

$ {\left\langle {U,V} \right\rangle _E} = \sum\limits_t {\left( {\frac{{\partial U}}{{\partial t}}\frac{{\partial V}}{{\partial t}} + v\nabla U \cdot v\nabla V} \right)} $

(7)

借助(4)式, U、V为正反传波场时, 能量范数成像条件^[26]为:

$ {I_E} = \sum\limits_{e,t} {\left( {\frac{{\partial U}}{{\partial t}}\frac{{\partial V}}{{\partial t}} + v\nabla U \cdot v\nabla V} \right)} $

(8)

在弹性介质假设条件下, 公式(8)中两个波场对应为弹性波场。

1.3 能量范数成像条件噪声压制分析

下面进一步分析上述成像条件噪声压制原理。将检波点波场与震源波场看作是多方向的平面波合成, 分别将震源和检波点波场变换至频率域, 互相关成像条件可以表示为:

$ I\left( \omega \right) = \sum\limits_{e,\omega } {{{\tilde U}_{\rm{s}}}\tilde U_{\rm{r}}^ * } $

(9)

式中:I(ω)为频率域互相关成像结果; $\widetilde{U}_{s} \widetilde{U}_{\mathrm{r}}^{*} $为频率域震源波场和检波点共轭波场的乘积。同理能量范数成像条件公式(8)可以表示为:

$ \begin{array}{l} {I_E}\left( k \right) = \sum\limits_{e,\omega } {\left[ { - {v^2}\left( {{\mathit{\boldsymbol{k}}_{\rm{s}}} \cdot {\mathit{\boldsymbol{k}}_{\rm{r}}}} \right) + {\omega ^2}} \right]{{\tilde U}_{\rm{s}}}\tilde U_{\rm{r}}^ * } \\ \;\;\;\;\;\;\;\;\; = \sum\limits_{e,\omega } {\left[ { - {v^2}\left\| {{k_{\rm{s}}}} \right\|\left\| {{k_{\rm{r}}}} \right\|\cos 2\theta + {\omega ^2}} \right]{{\tilde U}_{\rm{s}}}\tilde U_{\rm{r}}^ * } \end{array} $

(10)

式中: k_s、k_r分别为震源波场和检波点波场的波数向量; θ为波场波数向量夹角。利用频散关系:

$ \left\| {{\mathit{\boldsymbol{k}}_{\rm{s}}}} \right\| = \left\| {{\mathit{\boldsymbol{k}}_{\rm{r}}}} \right\| = \frac{\omega }{v} $

(11)

得到能量范数成像条件^[26]为:

$ {I_E}\left( k \right) = \sum\limits_{e,\omega } {\left[ { - {\omega ^2}\cos 2\theta + {\omega ^2}} \right]{{\tilde U}_{\rm{s}}}\tilde U_{\rm{r}}^ * } $

(12)

在ω²项前加入一个cos2θ_c项, 将其近似表示为:

$ {I_E}\left( k \right) \approx \sum\limits_{e,\omega } {\left[ { - {\omega ^2}\cos 2\theta + {\omega ^2}\cos 2{\theta _c}} \right]{{\tilde U}_{\rm{s}}}\tilde U_{\rm{r}}^ * } $

(13)

当选择θ_c与θ相等时, 则该夹角的能量范数为0, 将其反变换回空间域可以得到:

$ {I_E} = \sum\limits_{e,t} {\left[ {\cos \left( {2{\theta _c}} \right)\frac{{\partial {U_{\rm{s}}}}}{{\partial t}}\frac{{\partial {U_{\rm{r}}}}}{{\partial t}} + v\nabla {U_{\rm{s}}} \cdot v\nabla {U_{\rm{r}}}} \right]} $

(14)

因此, 可以选择θ_c来压制震源波场与检波点波场任意角度互相关造成的干扰噪声。互相关成像条件的后散射干扰可以近似看作两个波场成180°互相关干扰, 这里设置θ_c为90°入射, 即得到最终压制反向散射干扰成像条件为:

$ {I_E} = \sum\limits_{e,t} {\left( { - \frac{{\partial {U_s}}}{{\partial t}}\frac{{\partial {U_r}}}{{\partial t}} + v\nabla {U_{\rm{s}}} \cdot v\nabla {U_{\rm{r}}}} \right)} $

(15)

将其应用于二维弹性波逆时偏移成像中得到垂直分量和水平分量成像条件:

$ \begin{array}{*{20}{c}} {{{\left( {{I_E}} \right)}_x} = \sum\limits_{e,t} {\left\{ { - \frac{{\partial {U_{{\rm{s}}x}}}}{{\partial t}}\frac{{\partial {U_{{\rm{r}}x}}}}{{\partial t}} + \left( {v_{\rm{P}}^2 - v_{\rm{S}}^2} \right) \cdot } \right.} }\\ {\left( {\frac{{\partial {U_{{\rm{s}}x}}}}{{\partial x}} + \frac{{\partial {U_{{\rm{s}}z}}}}{{\partial z}}} \right)\left( {\frac{{\partial {U_{{\rm{r}}x}}}}{{\partial x}} + \frac{{\partial {U_{{\rm{r}}z}}}}{{\partial z}}} \right) + }\\ {\left. {v_{\rm{S}}^2\left( {\frac{{\partial {U_{{\rm{s}}x}}}}{{\partial x}} + \frac{{\partial {U_{{\rm{s}}z}}}}{{\partial x}}} \right)\left( {\frac{{\partial {U_{{\rm{r}}x}}}}{{\partial x}} + \frac{{\partial {U_{{\rm{r}}z}}}}{{\partial x}}} \right)} \right\}} \end{array} $

(16)

$ \begin{array}{*{20}{c}} {{{\left( {{I_E}} \right)}_z} = \sum\limits_{e,t} {\left\{ { - \frac{{\partial {U_{{\rm{s}}z}}}}{{\partial t}}\frac{{\partial {U_{{\rm{r}}z}}}}{{\partial t}} + \left( {v_{\rm{P}}^2 - v_{\rm{S}}^2} \right) \cdot } \right.} }\\ {\left( {\frac{{\partial {U_{{\rm{s}}x}}}}{{\partial x}} + \frac{{\partial {U_{{\rm{s}}z}}}}{{\partial z}}} \right)\left( {\frac{{\partial {U_{{\rm{r}}x}}}}{{\partial x}} + \frac{{\partial {U_{{\rm{r}}z}}}}{{\partial z}}} \right) + }\\ {\left. {\left. {v_{\rm{S}}^2\left( {\frac{{\partial {U_{{\rm{s}}x}}}}{{\partial z}} + \frac{{\partial {U_{{\rm{s}}z}}}}{{\partial z}}} \right)\frac{{\partial {U_{{\rm{r}}x}}}}{{\partial z}} + \frac{{\partial {U_{{\rm{r}}z}}}}{{\partial z}}} \right)} \right\}} \end{array} $

(17)

式中:(I_E)_x、(I_E)_z分别为能量范数成像条件偏移结果的水平分量和垂直分量; U_sx、U_sz分别为震源波场x分量和z分量; U_rx、U_rz分别为检波点波场x分量和z分量。

1.4 能量范数成像条件弹性波逆时偏移

公式(16)和公式(17)利用一阶速度应力弹性波方程求取时需引入中间变量, 首先假设:

$ \begin{array}{*{20}{c}} {{\varepsilon _{xx}} = \frac{{\partial {U_x}}}{{\partial x}}}&{{\varepsilon _{zz}} = \frac{{\partial {U_z}}}{{\partial z}}}\\ {{\varepsilon _{xz}} = \frac{{\partial {U_x}}}{{\partial z}}}&{{\varepsilon _{zx}} = \frac{{\partial {U_z}}}{{\partial x}}} \end{array} $

(18)

对(18)式取时间一阶导数, 变换得到:

$ \begin{array}{*{20}{c}} {\frac{{\partial {\varepsilon _{xx}}}}{{\partial t}} = \frac{\partial }{{\partial t}}\frac{{\partial {U_x}}}{{\partial x}} = \frac{{\partial {v_x}}}{{\partial x}}}\\ {\frac{{\partial {\varepsilon _{zz}}}}{{\partial t}} = \frac{\partial }{{\partial t}}\frac{{\partial {U_z}}}{{\partial z}} = \frac{{\partial {v_z}}}{{\partial z}}}\\ {\frac{{\partial {\varepsilon _{xz}}}}{{\partial t}} = \frac{\partial }{{\partial t}}\frac{{\partial {U_x}}}{{\partial z}} = \frac{{\partial {v_x}}}{{\partial z}}}\\ {\frac{{\partial {\varepsilon _{zx}}}}{{\partial t}} = \frac{\partial }{{\partial t}}\frac{{\partial {U_z}}}{{\partial x}} = \frac{{\partial {v_z}}}{{\partial x}}} \end{array} $

(19)

由(19)式假设的中间分量可以通过速度分量空间偏导求取, 对(19)式利用时间二阶差分和空间伪谱法进行数值离散, 有:

$ \begin{array}{*{20}{c}} {\varepsilon _{xx}^{n + 1} = \varepsilon _{xx}^n + \Delta t\frac{{\partial {v_x}}}{{\partial x}}}\\ {\varepsilon _{zz}^{n + 1} = \varepsilon _{zz}^n + \Delta t\frac{{\partial {v_z}}}{{\partial z}}}\\ {\varepsilon _{xz}^{n + 1} = \varepsilon _{xz}^n + \Delta t\frac{{\partial {v_x}}}{{\partial z}}}\\ {\varepsilon _{zx}^{n + 1} = \varepsilon _{zx}^n + \Delta t\frac{{\partial {v_z}}}{{\partial x}}} \end{array} $

(20)

式中: $ \varepsilon_{x x}^{n+1}、\varepsilon_{z z}^{n+1}、\varepsilon_{x z}^{n+1}、\varepsilon_{z x}^{n+1}$为n+1时刻应变波场值; $ \varepsilon_{x x}^{n}、\varepsilon_{z z}^{n}、\varepsilon_{x z}^{n}、\varepsilon_{z x}^{n}$为n时刻应变波场值; v_x、v_z分别为震源波场速度分量的水平分量和垂直分量。按照(20)式可求取任意时刻$ \partial U_{x} / \partial x、\partial U_{z} / \partial x、\partial U_{z} / \partial x 、\partial U_{z} / \partial z$而得到最终成像结果。

1.5 伪谱法数值离散求解方式

前文给出了能量范数成像条件求解方式, 利用(3)式、(16)式、(17)式、(18)式、(19)式和(20)式可实现基于能量范数成像条件逆时偏移。图 1给出了该成像条件对应的完整一阶速度应力弹性波方程波场延拓各分量更新过程, 其中$\sigma_{x x}^{n}、\sigma_{z z}^{n}、\sigma_{x z}^{n}、v_{x}^{n}、v_{z}^{n} $和$ \sigma_{x x}^{n+1}、\sigma_{z z}^{n+1}、\sigma_{x z}^{n+1}、v_{x}^{n+1}、v_{z}^{n+1}$分别为各分量n时刻和n+1时刻值。

2 模型试算分析

为验证算法正确性, 设计两层介质模型进行试算分析, 模型参数如图 2所示, 密度由一般经验公式计算得到。对单炮正演得到地震记录进行偏移, 震源采用30 Hz雷克子波, 位置为(1 000 m, 0);201道接收, 道间距10 m, 检波点深度为0, 0.5 ms采样, 1.0 s记录长度。

图 3为利用伪谱法弹性波互相关成像条件和能量范数成像条件对单炮记录进行偏移的结果。从图 3a和图 3b可以看出, 由于直达波近似水平方向传播与检波点波场绕射产生的P波和S波产生的干扰被较好压制(箭头处), 在界面上方的大部分绕射产生的干扰很大程度上得到了消除; 从图 3c和图 3d可以看出, 椭圆部分低频噪声部分被压制。多炮偏移叠加结果如图 4所示。从图 4a和图 4b可以看出, 界面上方的低频散射干扰明显降低; 图 4c和图 4d验证了能量范数成像条件的有效去噪效果。

设计凹陷模型(图 5)进行多分量逆时偏移试算, 分别测试互相关成像条件的偏移、纵横波解耦的弹性波偏移、能量范数成像条件的弹性波逆时偏移和多炮偏移叠加结果, 如图 6所示。与直接利用v_x、v_z分量进行互相关成像(图 6a, 图 6b)相比, 基于纵横波分离的成像PP波与SS波偏移结果更为清晰, 第二界面上方的低频干扰明显降低(图 6c, 图 6f)。PS波与SP波的偏移效果则较差(图 6d, 图 6e), 说明不同类型波场的直接互相关成像方法带来了更多噪声, 甚至有效信息被覆盖, 而关于转换波的成像方法则还需要利用其它技术来完成。基于能量范数成像条件的偏移结果x分量与z分量存在的低频噪声明显降低(图 6g, 图 6h), 正反传波场的非界面处互相关噪声得到了较好的压制。

3 复杂模型试算

截取Marmousi模型的一部分进行本文方法的噪声压制测试, 模型大小为2 000 m×700 m, 网格大小5 m×5 m, 速度参数如图 7所示, 密度由一般经验公式计算得到。正演模拟共41炮, 炮点位置从0到2 000 m, 炮间距50 m。单炮401道接收, 道间距5 m, 分布于整个模型地表附近, 0.5 ms采样, 2 s记录长度。炮点、检波点深度均为0, 采用30 Hz雷克子波点震源(同时激发纵波和横波)。从最终偏移叠加剖面(图 8)可以看出, 与互相关成像条件偏移结果(图 8a, 图 8b)相比, 基于能量范数成像条件的结果(图 8g, 图 8h)中低频噪声得到了较好的压制; 与纵横波分离后的成像结果(图 8c, 图 8d, 图 8e, 图 8f)相比, 图 8g和图 8h的信噪比较高。

4 结论

本文从能量范数的角度详细分析了能量范数成像条件及其噪声压制原理, 将这一成像条件发展到了一阶速度应力弹性波逆时偏移方法中。通过引入中间应变分量进行所需分量更新, 实现了一阶速度应力弹性波方程逆时偏移方法能量范数成像条件的应用。通过复杂模型试算, 对比v_x分量、v_z分量互相关成像条件偏移结果, 本文方法能够压制低频噪声; 对比纵横波波场分离偏移, 本文方法无需考虑转换波极性反转问题, 避免进行解耦波动方程延拓。