0 引言

叠前深度偏移是目前油气地震勘探资料处理中的关键技术。获取叠前深度偏移地震记录要经过炮记录采集或正演、常规处理和偏移归位处理等诸多环节。而地震正演与偏移成像存在计算量大、运行周期长等问题，且常规处理和偏移流程较为繁琐。为提高正演和偏移的效率，研究人员做了许多工作：王德利等^[1]对黏弹介质中三维波动方程进行了有限差分MPI并行正演及验证；Micikevicius^[2]针对GPU的有限差分实现给出了快速算法；龙桂华等^[3]应用GPU实现了三维地震波交错网格有限差分模拟的加速计算；李博等^[4]、刘守伟等^[5]、Foltinek等^[6]研究了基于GPU加速的高阶有限差分叠前逆时偏移，取得了显著的成果；李振华等^[7]利用无记忆拟牛顿—模拟退火法对偏移算子方程进行求解，提出了一种正则化偏移成像的全局优化快速算法；屠宁^[8]对基于压缩感知的快速最小二乘逆时偏移进行了研究，降低了最小二乘逆时偏移中波动方程的计算量。

尽管如此，获取叠前深度偏移剖面流程繁琐的根本问题依旧未能得到解决。为简化这一复杂流程、快速获得偏移剖面，人们往往应用一维褶积进行偏移剖面的模拟，但对于二维、三维复杂速度模型模拟效果不甚理想^[9]。Hamran等^[10-11]、Gelius等^[12-14]、Wu等^[15]、Zhu等^[16]、朱小三等^[17]在广义散射层析成像的研究中详细介绍了波场的成像条件及成像分辨率算子，实现了高精度的散射层析成像反演；马在田^[18]研究了地震偏移成像分辨率与观测系统、成像点、震源等的关系，给出地震成像分辨率的定量计算公式，阐述了观测地震道的分辨率和成像分辨率的空间变化的特征。

Lecomte等^[19-22]分析了叠前深度偏移成像分辨率算子，提出用散射波数算子表征地震成像分辨率函数，并以此实现局部小区域的叠前深度直接模拟成像。在此基础上，本文发展了一种基于深度域广义卷积技术的叠前深度偏移地震记录直接模拟方法，在已知速度模型和野外采集参数的前提下，无需经过波场正演和偏移处理，便可以直接快速生成偏移剖面，得到理论成像结果。根据理论成像结果，可以检测观测系统的优劣和地震数据的质量，为优化观测系统参数和评价成像质量提供参考，以满足数据采集和处理中质量控制的需求^[23-25]。

1 方法原理

对地下介质正演和成像过程可以用算子的形式分别表示为

$ U\left( {\omega , {\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}} \right) = F\left( {\omega , {\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}\mid \mathit{\boldsymbol{r}}} \right)S(\mathit{\boldsymbol{r}}) $

(1)

$ \begin{array}{l} I\left( {{\mathit{\boldsymbol{r}}^\prime }} \right) = B\left( {{\mathit{\boldsymbol{r}}^\prime }\mid \omega , {\mathit{\boldsymbol{r}}_{\rm{s}}}, {\mathit{\boldsymbol{r}}_{\rm{g}}}} \right)U\left( {\omega , {\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}} \right)\\ \;\;\;\;\;\;\;\; = B\left( {{\mathit{\boldsymbol{r}}^\prime }\mid \omega , {\mathit{\boldsymbol{r}}_{\rm{s}}}, {\mathit{\boldsymbol{r}}_{\rm{g}}}} \right)F\left( {\omega , {\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}\mid \mathit{\boldsymbol{r}}} \right)S(\mathit{\boldsymbol{r}}) \end{array} $

(2)

式中：F是模型的正演算子，用于生成模型S的正演数据U；B是成像算子，用于生成U的成像结果I；r表示模型中散射点位置；r′表示模型中散射点成像位置；r_g为检波点位置；r_s为震源位置(图 1)。由式(2)可知，通过算子B(r′|ω, r_s, r_g)F(ω, r_g, r_s|r)，可实现由模型S直接获取成像结果I。

标量声波方程可表示为

$ {\nabla ^2}P(\mathit{\boldsymbol{r}}, t) = \frac{1}{{C(\mathit{\boldsymbol{r}})}}\frac{{{\partial ^2}P(\mathit{\boldsymbol{r}}, t)}}{{\partial {t^2}}} $

(3)

式中：▽²为拉普拉斯算子；P(r，t)为地震波场；C(r)为速度场。将式(3)变换到频率域，并引入速度扰动函数O(r)=1-C₀(r)/C(r)，在Born近似下可得

$ {P_{{\rm{sc}}}}\left( {{\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}} \right) = - \int_v {{P_{{\rm{in}}}}} \left( {\mathit{\boldsymbol{r}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}} \right)k_0^2{\rm{O}}(\mathit{\boldsymbol{r}})G\left( {{\mathit{\boldsymbol{r}}_{\rm{g}}}, \mathit{\boldsymbol{r}}} \right){\rm{d}}\mathit{\boldsymbol{r}} $

(4)

式中：速度扰动函数O(r)一定程度上可近似看作反射率；C₀(r)为背景速度场；k₀=ω/C₀为背景波数；P_sc、P_in分别代表散射波场和入射波场；v为散射体范围。对于一个非均匀且无衰减的背景场，可以使用动态射线追踪得到其格林函数G(r_g，r)^[26]，结合Stamnes^[27]给出的格林函数近似解，并引入一个偏移算子l(r_g, r_s, r′)，则此时散射点的成像结果可表示为

$ \begin{array}{l} \hat O\left( {{\mathit{\boldsymbol{r}}^\prime }} \right) = \int_v O (\mathit{\boldsymbol{r}})\left\{ {l\left( {{\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}, {\mathit{\boldsymbol{r}}^\prime }} \right) \times } \right.\\ \left. {\;\;\;\;\;\;\;\exp \left\{ {{\rm{i}}\omega \left[ {\tau \left( {{\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}, \mathit{\boldsymbol{r}}} \right) - \tau \left( {{\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}, {\mathit{\boldsymbol{r}}^\prime }} \right)} \right]} \right\}{\rm{d}}{\mathit{\boldsymbol{r}}_{\rm{g}}}{\rm{d}}{\mathit{\boldsymbol{r}}_{\rm{s}}}} \right\}{\rm{d}}\mathit{\boldsymbol{r}}\\ \;\;\;\;\;\;\; = \int_v O (\mathit{\boldsymbol{r}})D\left( {\mathit{\boldsymbol{r}}, {\mathit{\boldsymbol{r}}^\prime }} \right){\rm{d}}\mathit{\boldsymbol{r}} \end{array} $

(5)

式中：D(r, r′)即是直接成像算子，本文将其称为点扩散算子；a和τ为

$ \left\{ {\begin{array}{*{20}{l}} {a\left( {{\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}, \mathit{\boldsymbol{r}}} \right) = - \frac{{\gamma {k_0}}}{2}\sqrt {\frac{{{\rm{i}}{k_0}}}{{2\pi }}} \frac{1}{{\sqrt {\left| {\mathit{\boldsymbol{r}} - {\mathit{\boldsymbol{r}}_{\rm{s}}}} \right|} \cdot \sqrt {\left| {{\mathit{\boldsymbol{r}}_{\rm{g}}} - \mathit{\boldsymbol{r}}} \right|} }}}\\ {\tau \left( {{\mathit{\boldsymbol{r}}_{\rm{g}}}, {\mathit{\boldsymbol{r}}_{\rm{s}}}, \mathit{\boldsymbol{r}}} \right) = \frac{{\left| {\mathit{\boldsymbol{r}} - {\mathit{\boldsymbol{r}}_{\rm{s}}}} \right| + \left| {{\mathit{\boldsymbol{r}}_{\rm{g}}} - \mathit{\boldsymbol{r}}} \right|}}{{{{\rm{C}}_0}}}} \end{array}} \right. $

(6)

其中γ表征振幅项。

结合Beylkin等^[28-29]提出的理论，在成像点的邻域内将振幅项和相位项Talor展开，结合程函方程并使点扩散算子趋近于δ函数，同时为了得到携带地震波响应的叠前深度偏移剖面，需将地震子波作用在点扩散算子上，则有

$ D\left( {\mathit{\boldsymbol{r}}, {\mathit{\boldsymbol{r}}^\prime }} \right) = \int {\int {s(\omega )} } \exp \left[ {{\rm{i}}\mathit{\boldsymbol{K}}\left( {\mathit{\boldsymbol{r}} - {\mathit{\boldsymbol{r}}^\prime }} \right)} \right]{{\rm{d}}^2}\mathit{\boldsymbol{K}} $

(7)

式中：s(ω)为地震子波的频谱；K=-k_g+k_s=ωτ(r_g, r_s, r′)，其中k_s为沿着震源点到成像点的射线路径上的波数向量，k_g为沿着成像点到检波器的射线路径上的波数向量，且有|k_s|=|k_g|=k₀。由此，可以使用地震子波信息、观测系统参数和射线路径信息在波数域构建点扩散算子，与地下速度模型的速度扰动或反射率信息相互作用，即可得到理论成像结果。同时，还可以将扰动点近似看成是一个小区域内的照明点，即扰动点处的点扩散算子可以对扰动点临近区域进行局部模拟成像。

2 点扩散算子分析

由上文可知，将点扩散算子与速度扰动函数褶积即可得到模拟成像结果。速度扰动函数可以近似看作反射率，由地下速度模型直接得出。空间域的褶积运算相当于波数域的乘积运算，点扩散算子可以在波数域由震源信息、观测系统信息及其射线路径信息构建。为说明具体算法和效果，本文以主频为25Hz的Ricker子波作为震源，起始炮点位置为50m，终止炮点位置为4600m，炮点距为50m，每炮共461道接收，道间距为10m，计算了Marmousi模型(图 2)中四个点处的点扩散算子(图 3)，以此观察点扩散算子的形态。由图 3可以看出，点扩散算子随着深度的增加，地震波频谱带宽变窄、波长变长，表明成像分辨率逐渐降低。图 4、图 5展示了(1500m，2000m)附近局部区域(1000m×1000m)的模拟成像过程，可以看出，在给定速度模型信息、震源信息、观测系统信息的情况下, 构建的点扩散算子可以对地下构造进行准确模拟成像。

3 模型试算

通过简单凹槽模型和Marmousi模型，验证叠前深度偏移地震记录直接模拟方法的准确性及其在指导观测系统设计方面的高效性。

3.1 简单凹槽模型

建立图 6所示的凹槽模型，尺寸为200×300个网格，网格间距Δx=Δz=5m。采用主频为30Hz的Ricker子波作为震源进行试算，起始炮点位置为5m，终止炮点位置为1000m，炮点距为50m(共20炮)，每炮共200道接收，道间距5m。叠前深度偏移地震记录直接模拟结果(图 7)能够准确反映地下层位位置、形态以及对应的波形信息。

3.2 Marmousi模型

为了检验算法的适用性和准确性，使用图 2所示的Marmousi速度模型进行试算。该模型尺寸为461×284个网格，网格间距Δx=Δz=10m。采用主频为25Hz的Ricker子波作为震源，起始炮点位置为50m，终止炮点位置为4600m，炮点距为50m (共92炮)，每炮共461道接收，道间距为10m。图 8为叠前深度偏移地震记录直接模拟结果，能够准确反映地下地质构造信息。与地下反射系数深度域褶积波形相似度很高(图 9)，相关系数为0.9494，验证了叠前深度偏移地震记录直接模拟方法的准确性。

图 10为相同观测参数下Marmousi模型不同方法的叠前深度偏移结果(正演数据由空间八阶、时间二阶的波动方程有限差分正演方法计算)。与图 8对比可以看出：逆时偏移成像结果在深部细节上，尤其是深部低速目标区以及中部层状背斜构造处，成像效果较差(图 10a深部方框所示)，且浅部存在明显的噪声干扰(图 10a浅部方框所示)；最小二乘逆时偏移能够对深部细节进行准确成像，但深部与浅部反射率信息对应性较差(图 10b方框所示)；Kirchhoff偏移在中部层状背斜构造和深部低速层附近的成像效果明显较差，不能够完整准确描述地下真实构造(图 10c方框所示)。由此可见，叠前深度偏移地震记录直接模拟方法能够为成像质量提供评价的参考标准。

采用主频为25Hz的Ricker子波作为震源，起始炮点位置为0，炮点距分别为100、50、30、20、10m，炮点右侧接收，最小炮检距为10m，最大炮检距为1000m，道间距为10m，共100道接收，计算不同覆盖次数(分别为5、10、17、25、50)的直接模拟成像结果。图 11展示了不同观测系统参数下Marmousi模型深度2500m附近低速目标区域叠前深度偏移地震记录直接模拟结果，覆盖次数越高的观测系统成像效果越好。炮点距为20m(覆盖次数为25)时成像准确(图 11d)，成像剖面第340道波形与对应反射系数深度域褶积波形相关系数为0.8190(图 12)；炮点距为30m(覆盖次数为17)时成像结果(图 11c)也可识别出目标层及其附近构造，第340道波形与对应反射系数深度域褶积波形相关系数为0.7878(图 12)。兼顾经济、效率两方面要求，建议采用20~30m的炮点距进行地震采集。

4 结束语

本文研究了一种在已知速度模型和野外采集参数的前提下直接模拟叠前深度偏移记录的方法。该方法可以简单高效地获得叠前深度偏移的理论成像结果，省去了繁琐的炮记录正演、常规处理及偏移等一系列过程。其模拟成像精度能够准确反映地质构造信息及反射率信息，可为评价观测系统和成像质量提供参考标准，以满足采集观测系统设计、成像处理和储层分析等领域的研究需求。