0 引言

随着油气勘探技术的发展，如何进行复杂岩性油气藏成像已成为一个亟待解决的问题。叠前深度偏移在复杂构造成像中起着至关重要的作用。逆时偏移(reverse time migration，简称RTM)是一种基于双程波动方程的先进叠前深度偏移方法，能利用回转波和多次反射波对复杂构造进行成像。RTM的主要特点是能够解决任意横向变速和高陡倾角等复杂地质构造成像问题，特别是盐丘和潜山侧翼成像。因此，该方法克服了常规地震偏移(Kirchhoff积分法偏移、有限差分法偏移、频率波数域波动方程偏移等)对高陡倾角和横向变速大等复杂地质构造成像效果不好的缺陷，但同时会产生低频噪声污染和大量的计算量、存储量(Fletcher et al，2006；Guitton et al，2007；Symes，2007；Chattopadhyay et al，2008；Jones，2008；Costa et al，2009；Xu et al，2011；Nguyen et al，2013；Weibull et al，2013)。

在实际地震数据采集中，由于受采集成本或地形的约束，一般情况下道间距和时间采样间隔较大，这会影响RTM成像波形的光滑性。对于RTM成像构造解释，波形的光滑性通常可以忽略不计。然而对于储层成像精细解释，因其需要更高的分辨率和更精确的成像，故波形不光滑性的影响不能忽视。如果采用很小的空间采样间隔，可能有助于减少波形的不光滑性，但计算、存储成本随空间采样间隔呈4次方增长(三维)或3次方增长(二维)。因此，这在实践中是不可行的。

在本文中，我们首先分析了空间、时间采样间隔对RTM成像波形光滑性的影响；然后采用三次样条插值对不光滑的成像波形进行空间插值，以解决RTM波形不光滑的问题。最后利用复杂模型进行验证，并对RTM波形校正前后的地震属性(瞬时振幅、瞬时频率、瞬时相位等)进行了分析。

1 波形光滑性分析

RTM通常采用源波场和接收波场互相关作为成像条件，其公式可以写为(Symes，2007；Chattopadhyay et al，2008；Claerbout，1971)

$ \boldsymbol{I}(\overrightarrow{\boldsymbol{X}})=\int_{0}^{t_{\max }} P_{\mathrm{s}}(\overrightarrow{\boldsymbol{X}}, t) P_{\mathrm{r}}(\overrightarrow{\boldsymbol{X}}, t) \mathrm{d} t $

(1)

其中，$\boldsymbol{I}(\overrightarrow{\boldsymbol{X}}) $为RTM成像；P_s$ (\overrightarrow{\boldsymbol{X}}, t)$、P_r$(\overrightarrow{\boldsymbol{X}}, t) $分别为源波场和接收波场；$\overrightarrow{\boldsymbol{X}} $为空间坐标矢量；2-D：$ \overrightarrow{\boldsymbol{X}}$=(x，z)，3-D：$\overrightarrow{\boldsymbol{X}} $=(x，y，z)；t∈(0，t_max)，为时间。互相关成像条件(式(1))在源为脉冲时等价于准确的反射率成像公式，而在源为短区间子波时是反射率成像公式的一个合理近似(Claerbout，1971)。

为了计算源波场和接收波场，这里采用一阶应力-速度声波方程，其2-D公式为

$ \left\{\begin{aligned} \dot{v}_{\mathrm{x}} &=\frac{1}{\rho} \cdot \frac{\partial \sigma}{\partial x} \\ \dot{v}_\rm{z} &=\frac{1}{\rho} \cdot \frac{\partial \sigma}{\partial z} \\ \dot{\sigma} &=K\left(\frac{\partial v_{\mathrm{x}}}{\partial x}+\frac{\partial v_\rm{z}}{\partial z}\right) \end{aligned}\right. $

(2)

其中，v_x、v_z分别为水平、垂直质点振动速度；σ为声波应力场，这里代表源波场P_s$(\overrightarrow{\boldsymbol{X}}, t) $或接收波场P_r$(\overrightarrow{\boldsymbol{X}}, t) $；K=ρv_p²为杨氏模量，v_p为地震纵波速度；ρ为密度；字母上面的点表示对时间求导。

这里采用交错网格有限差分方法(时间二阶，空间四阶)进行求解，求解公式如下(Yang et al，2015)

$ \left\{ \begin{matrix} \frac{vx_{i, k}^{n+1}-vx_{i, k}^{n}}{\Delta t}=\frac{1}{\rho }\cdot \frac{-{{a}_{1}}\sigma _{i+3/2, k}^{n+1/2}+{{a}_{2}}\sigma _{i+1/2, k}^{n+1/2}-{{a}_{2}}\sigma _{i-1/2, k}^{n+1/2}+{{a}_{1}}\sigma _{i-3/2, k}^{n+1/2}}{\Delta x} \\ \frac{vz_{i, k}^{n+1}-vz_{i, k}^{n}}{\Delta t}=\frac{1}{\rho }\cdot \frac{-{{a}_{1}}\sigma _{i, k+3/2}^{n+1/2}+{{a}_{2}}\sigma _{i, k+1/2}^{n+1/2}-{{a}_{2}}\sigma _{i, k-1/2}^{n+1/2}+{{a}_{1}}\sigma _{i, k-3/2}^{n+1/2}}{\Delta z} \\ \frac{\sigma _{i, k}^{n+1/2}-\sigma _{i, k}^{n-1/2}}{\Delta t}\text{=}K\left(\begin{matrix} \frac{-{{a}_{1}}vx_{i+3/2, k}^{n}+{{a}_{2}}vx_{i+1/2, k}^{n}-{{a}_{2}}vx_{i-1/2, k}^{n}+{{a}_{1}}vx_{i-3/2, k}^{n}}{\Delta x} \\ +\frac{-{{a}_{1}}vz_{i, k+3/2}^{n}+{{a}_{2}}vz_{i, k+1/2}^{n}-{{a}_{2}}vz_{i, k-1/2}^{n}+{{a}_{1}}vz_{i, k-3/2}^{n}}{\Delta z} \\ \end{matrix} \right). \\ \end{matrix} \right. $

(3)

此处，$a_{1}=\frac{1}{24}, a_{2}=\frac{9}{8} $。

本文采用一个简单层状模型分析空间、时间采样间隔对RTM波形光滑性的影响。模型大小为800m×400m，上层速度为2000m/s，下层速度为3000m/s。震源采用主频45Hz的Ricker子波，震源位置在(400m，0m)，检波器位置在(200m，0m)。我们首先计算和分析了在相同时间采样间隔(0.25ms)的情况下，不同空间采样间隔(1m×1m、2m×2m、4m×4m)时的RTM成像和波形(图 1)。从图 1可以看出，空间采样间隔对波形光滑性的影响较大，随着空间采样间隔的增大，波形的不光滑性增加，即使在空间网格为1m×1m时，在中心峰值处仍然看到波形的不光滑性。然后，计算和分析了在相同空间采样间隔(4m×4m)的情况下，不同时间采样间隔(0.25、1.00、4.00ms)时的RTM成像和波形(图 2)。从图 2可以看出，几种波形不光滑程度相当，也就是说时间采样间隔对波形的光滑性影响很小，当然其前提是采样间隔不能过大而引起频散。对于基于互相关成像条件的深度偏移，波场采样间隔应为采样定理所给值的50%(Zhang et al，2003)。

2 波形校正方法

在实际数据采集中，常常采用较大的空间采样间隔来减少RTM计算、存储成本，但这会导致RTM成像波形不光滑。如果采用先插值后成像的方法，势必造成RTM计算、存储量过大，成本太高。因此，本文采用先成像后插值的方法，利用三次样条函数对不光滑的RTM波形进行空间插值，以得到光滑的波形。为了方便起见，只考虑垂向深度插值，因为在实际数据采集中主要考虑垂向波形。RTM波形校正公式如下

$ \begin{aligned} \tilde{I}(z) &=M_{i} \frac{\left(z_{i+1}-z\right)^{3}}{6 h_{i}}+M_{i+1} \frac{\left(z-z_{i}\right)^{3}}{6 h_{i}}+\left[I\left(z_{i}\right)-\frac{M_{i} h_{i}^{2}}{6}\right] \frac{z_{i+1}-z}{h_{i}}+\\ &\left[I\left(z_{i+1}\right)-\frac{M_{i+1} h_{i}^{2}}{6}\right] \frac{z-z_{i}}{h_{i}}, (i=0, 1, \cdots, n) \end{aligned} $

(4)

其中，I(z_i)为插值前的RTM成像；$\tilde{I}(z) $为插值后的RTM成像；z∈ [z_i，z_i+1]为垂向插值区间；h_i=z_i+1-z_i为垂向空间插值间隔。参数M_i计算公式为

$ \left(\begin{array}{ccccc}{2} & {\lambda_{1}} & {} & {} & {} \\ {\mu_{2}} & {2} & {\lambda_{2}} & {} & {} \\ {} & {\ddots} & {\ddots} & {\ddots} & {} \\ {} & {} & {\mu_{n-2}} & {2} & {\lambda_{n-2}} \\ {} & {} & {} & {\mu_{n-1}} & {2}\end{array}\right)\left(\begin{array}{c}{M_{1}} \\ {M_{2}} \\ {\vdots} \\ {M_{n-2}} \\ {M_{n-1}}\end{array}\right)=\left(\begin{array}{c}{d_{1}-\mu_{1} M_{0}} \\ {d_{2}} \\ {\vdots} \\ {d_{n-2}} \\ {d_{n-1}-\lambda_{n-1} M_{n}}\end{array}\right). $

(5)

这里采用自然边界条件

$ \left\{\begin{array}{l}{I^{\prime \prime}\left(z_{0}\right)=M_{0}=0} \\ {I^{\prime \prime}\left(z_{n}\right)=M_{n}=0}\end{array}\right. $

(6)

此处，M₀、M_n为RTM成像在2个端点处关于深度的二阶导数值。

其他参数计算公式如下

$ \begin{align} & \left\{ \begin{array}{*{35}{l}} {{\mu }_{i}}=\frac{{{h}_{i-1}}}{{{h}_{i-1}}+{{h}_{i}}} \\ {{\lambda }_{i}}=1-{{\mu }_{i}} \\ {{d}_{i}}=6I\left[ {{z}_{i-1}}, {{z}_{i}}, {{z}_{i+1}} \right] \\ I\left[ {{z}_{i-1}}, {{z}_{i}}, {{z}_{i+1}} \right]=\frac{I\left[ {{z}_{i}}, {{z}_{i+1}} \right]-I\left[ {{z}_{i-1}}, {{z}_{i}} \right]}{{{z}_{i+1}}-{{z}_{i-1}}} \\ I\left[ {{z}_{i}}, {{z}_{i+1}} \right]=\frac{I\left({{z}_{i+1}} \right)-I\left({{z}_{i}} \right)}{{{h}_{i}}} \cdot \\ \end{array} \right. \\ & (i=1, 2, \cdots n-1) \\ \end{align} $

(7)

图 3为简单模型RTM成像波形校正前后的波形对比图。从图 3可以看出，波形校正方法可以很好地解决波形不光滑问题。

3 波形校正应用

在进行储层解释时，通常选取储层顶底地震属性横向空变稳定性较好的反射层作为标准层位，然后进行储层沉积演化解释，其中波形光滑性程度直接影响反射层位地震属性空变稳定性。这里根据实际地质构造情况，构建了一个带有3个散射点的单斜模型(图 4(a))。假定散射点代表储层，其上下的单斜构造为标准层位，其可当作储层的顶底，这样就可以通过地震属性技术进行储层沉积演化解释。如果标准层位地震属性具有稳定的空变特性，将会得到更好的储层精细解释精度，以减少地震解释误差。单斜模型大小为21975m×5000m，空间采样间隔12.50m×6.25m。3个散射点P₁、P₂、P₃速度分别为3500、4000、2400m/s。模型中单斜构造速度从上到下依次为2200、2600、3000、3400、3800、5000m/s。震源子波采用主频45Hz的Ricker子波，第1炮震源位于(9975m，0m)，从左往右移动，炮间距为50m，共有241炮，每炮799道，道间距为12.5m。图 4(b)、4(c)分别为波形校正前后的RTM成像对比，从图 4(b)、4(c)可以看出，构造成像差别不大。图 5为波形校正前后的RTM成像局部放大和散射点P₂波形对比。从图 5可以看出，波形差别较大，原始成像波形不光滑，经过波形校正后，波形变得光滑。图 6为波形校正前后的L2层拉平剖面和沿层地震属性(瞬时振幅、瞬时频率、瞬时相位)空变稳定性对比，因为是二维剖面，故沿层地震属性只反映一个方向(横向)的空间变化特征。从图 6可以看出，波形校正后，瞬时振幅空变稳定性变化不大，但瞬时频率和瞬时相位横向空变稳定性明显变好，特别是瞬时相位更明显一些。

4 结论与讨论

本文对RTM波形光滑性问题进行了研究，针对稀疏网格导致的RTM波形不光滑问题，选取三次样条插值方法进行垂向空间插值来进行波形校正。利用简单、复杂模型进行验证，结果均表明本文给出的波形校正方法不仅提高了波形的光滑度，而且提高了沿层地震属性(特别是瞬时频率和瞬时相位)横向空变稳定性，从而可以提高储层精细解释的精度。本文虽然只考虑了空间垂向深度波形插值问题，但该方法可以推广到二维、三维空间波形插值情形。同样，本文探讨的二维RTM波形问题，同样可以推广到三维RTM，这样就可以提取三维数据体的沿等时面的地震属性，更好地看出波形校正对地震属性空变稳定性的影响，当然这需要更大的计算量和存储量及进一步的研究。

致谢: 在本文研究过程中得到中国石油集团东方地球物理勘探有限责任公司王永军博士的支持和帮助，在此表示感谢。同时审稿人的认真审阅和提出的宝贵意见，以及编辑的认真负责，使得本文更加严谨和科学，在此一并表示感谢。