0 引言

为了处理油气勘探中日益复杂的地下介质情况以及产生的新问题，新的偏移理论不断被提出，更先进的成像技术也相应发展并应用到各个领域，不仅提高了计算效率，成像精度也在不断提高。但是获取叠前深度偏移剖面时流程繁琐的根本问题依旧未能得到解决。为简化这一复杂流程、快速获得高精度偏移剖面，研究人员也做了很多工作：Fletcher等^[1]首先提出了深度域等效子波的概念，利用等效子波完成了快速成像；周东勇等^[2]在总结了已有的研究基础上推导了深度域等效子波的表达式，并基于这一理论完成了深度域合成地震记录；向晨等^[3]、郝洪鉴等^[4]相继完成了基于点扩散函数的深度域快速成像技术以及深度域反演；段伟国等^[5]、张庆臣^[6]及姚振岸等^[7]分别完成了基于点扩散函数的迭代加权最小二乘反演、高精度逆时偏移和预条件黏声最小二乘偏移；岳玉波等^[8]利用射线理论计算点扩散函数；袁江洋等^[9]对深度域合成地震记录的精度问题进行了分析。

快速成像技术虽然在不断发展和完善，但是成像结果仍需考虑很多方面，成像理论也有待完善。Gelius等^[10-11]提出广义散射层析方法，率先引入成像算子；Hamran等^[12-13]也实现了基于局部平面波的散射层析成像；徐海等^[14]系统分析了现有地震数据高分辨率处理方法的适用性及不足，探索了一种基于高频噪声定量约束函数的高分辨率处理方法，有效且保真地提高地震资料信噪比与分辨率；Lecomte等^[15-18]分析了叠前深度偏移成像分辨率算子，提出用散射波数算子表征地震成像分辨率函数，并以此实现局部小区域的叠前深度地震记录直接模拟成像；吴涵等^[19-20]提出了一种基于深度域广义卷积技术的叠前深度偏移地震记录直接模拟方法；焦峻峰等^[21]对模拟成像方法进行了改善，考虑了地震子波在深度域中的空变。

在已有基础上，本文发展了一种基于粗网格的多算子叠前深度直接模拟成像方法，此方法对已知的速度模型先进行粗网格离散，构建粗成像算子，测试并选取最优插值方法对算子进行插值，使用插值后成像算子进行直接模拟成像，以提高计算效率。后续又在粗网格的基础上计算了多算子并进行模拟成像，使用多个算子控制模拟成像，得到更符合理论的成像结果，对观测系统的建立与优化有一定的指导意义。

1 模拟成像方法原理

模拟叠前局部成像方法是一种卷积技术，不需要首先生成合成数据，然后进行偏移以计算储层模型的偏移地震响应。根据测量、激发子波、波模式和局部储层结构，快速获得模拟叠前偏移剖面。此方法源于对地震分辨率^{[12, 14-15]}和地震成像新方法^[16-18]进行的研究，最初是为了更好地理解偏移剖面而开发的。

研究此方法，就必须首先定义模拟叠前局部成像方法的关键参数，即散射波数向量$ ( $图 1$ ) $，可使用散射点处的入射向量与反射向量进行表示

$ \boldsymbol{K}(\omega ,{r}_{0},{r}_{\mathsf{s}},{r}_{\mathsf{g}})={\boldsymbol{k}}_{\mathsf{g}}(\omega ,{r}_{0},{r}_{\mathsf{g}})-{\boldsymbol{k}}_{\mathsf{s}}(\omega ,{r}_{\mathsf{s}},{r}_{0}) $

(1)

式中：$ {\boldsymbol{k}}_{\mathsf{s}} $为沿着激发点到散射点的射线路径的波数向量；$ {\boldsymbol{k}}_{\mathsf{g}} $为沿着成像点到接收点的射线路径的波数向量；$ \omega $表示角频率；r₀表示模型中散射点位置；$ {r}_{\mathsf{s}} $和$ {r}_{\mathsf{g}} $分别表示激发点位置与接收点位置。入射波数向量与反射波数向量共同决定了散射向量$ \boldsymbol{K} $，而$ \boldsymbol{K} $又控制着模拟叠前成像的结果。

如果分别用算子的形式描述模型的正演过程和成像过程，则有

$ \boldsymbol{U}(\omega ,{r}_{\mathsf{g}},{r}_{\mathsf{s}})=\boldsymbol{D}(\omega ,{r}_{\mathsf{g}},{r}_{\mathsf{s}})\boldsymbol{V}\left(r\right) $

(2)

式中：$ \boldsymbol{V} $表示速度模型；$ \boldsymbol{D} $表示正演算子，用于生成$ \boldsymbol{V} $的正演数据$ \boldsymbol{U} $；r表示模型中散射点位置。

引入成像算子$ \boldsymbol{T} $并作用在$ \boldsymbol{U} $上，则可以得到成像结果$ \boldsymbol{I} $

$ \begin{array}{l}\boldsymbol{I}\left({r}^{\text{'}}\right)=\boldsymbol{T}\left({r}^{\text{'}}\right|\omega ,{r}_{\mathsf{g}},{r}_{\mathsf{s}}\left)\boldsymbol{U}\right(\omega ,{r}_{\mathsf{g}},{r}_{\mathsf{s}})\\ =\boldsymbol{T}\left({r}^{\text{'}}\right|\omega ,{r}_{\mathsf{g}},{r}_{\mathsf{s}}\left)\boldsymbol{D}\right(\omega ,{r}_{\mathsf{g}},{r}_{\mathsf{s}}\left)V\right(r)\\ =\boldsymbol{F}(\omega ,{r}_{\mathsf{g}},{r}_{\mathsf{s}})\boldsymbol{V}\left(r\right)\end{array} $

(3)

式中$ \boldsymbol{F} $即直接成像算子，本文称其为点扩散算子，可直接对$ \boldsymbol{V} $进行模拟成像。标量声波方程可表示为

$ {\nabla }^{2}\boldsymbol{P}(r,t)=\frac{1}{{v}^{2}\left(r\right)}\frac{{\partial }^{2}\boldsymbol{P}(r,t)}{\partial {t}^{2}} $

(4)

式中：$ \boldsymbol{P}(r,t) $为地震波场；$ v\left(r\right) $为速度场。将式(4)变换到频率域，引入速度扰动函数$ O\left(r\right)=1-{v}_{0}\left(r\right)/v\left(r\right) $，其中$ {v}_{0}\left(r\right) $为背景速度场，在Born近似下得

$ \boldsymbol{P}_{\mathsf{s}\mathsf{c}}({r}_{\mathsf{g}},{r}_{\mathsf{s}})=-{\int }_{\boldsymbol{\varOmega }}\boldsymbol{P}_{in}(r,{r}_{s})\boldsymbol{k}_{0}^{2}O\left(r\right)\boldsymbol{G}({r}_{\mathsf{g}},r)\mathsf{d}r $

(5)

式中：$ \boldsymbol{G} $(r_g, r)为格林函数；$ \boldsymbol{k}_{0}=\omega /{v}_{0} $为背景波数；$ \boldsymbol{P}_{\mathsf{s}\mathsf{c}} $、$ \boldsymbol{P}_{\mathsf{i}\mathsf{n}} $分别代表散射波场和入射波场；Ω为散射体范围。反射率在一定程度上可近似为速度扰动函数$ O\left(r\right) $并参与运算。

对于非均匀且无衰减的背景场，使用动态射线追踪得到其格林函数$ \boldsymbol{G}({r}_{\mathsf{g}},r) $，结合Stamnes^[22]给出的格林函数近似解，并引入一个偏移算子$ \boldsymbol{m}({r}_{\mathsf{g}},{r}_{\mathsf{s}},{r}^{\text{'}}) $，则此时散射点的成像结果可表示为

$ \begin{array}{l}\stackrel{}{\widehat{\boldsymbol{O}}}={\int }_{\boldsymbol{\varOmega }}O\left(r\right)\iiint \frac{\boldsymbol{a}\left({r}_{\mathsf{g}},{r}_{\mathsf{s}},r\right)}{\boldsymbol{a}\left({r}_{\mathsf{g}},{r}_{\mathsf{s}},{r}^{\mathsf{\text{'}}}\right)}\boldsymbol{m}\left[\boldsymbol{a}\left({r}_{\mathsf{g}},{r}_{\mathsf{s}},{r}^{\mathsf{\text{'}}}\right)\right]\times \\ \mathsf{e}\mathsf{x}\mathsf{p}\left\{\mathsf{i}\omega \left[\boldsymbol{\tau }\right({r}_{g},{r}_{s},r)-\boldsymbol{\tau }({r}_{g},{r}_{s},r\mathsf{\text{'}}\left)\right]\right\}\mathsf{d}{r}_{\mathsf{g}}\mathsf{d}{r}_{\mathsf{s}}\mathsf{d}r\\ ={\int }_{\boldsymbol{\varOmega }}\boldsymbol{O}\left(r\right)\boldsymbol{F}(r,r\mathsf{\text{'}})\mathsf{d}r\end{array} $

(6)

式中$ \boldsymbol{\alpha } $和$ \boldsymbol{\tau } $分别为

$ \boldsymbol{\alpha }({r}_{\mathsf{g}},{r}_{\mathsf{s}},r)=-\frac{\gamma {\boldsymbol{k}}_{0}}{2}\sqrt[]{\frac{\mathsf{i}{\boldsymbol{k}}_{0}}{2\mathsf{\pi }}}\frac{1}{\sqrt[]{|r-{r}_{\mathsf{s}}|}\sqrt[]{|{r}_{\mathsf{g}}-r|}} $

(7)

$ \boldsymbol{\tau }({r}_{\mathsf{g}},{r}_{\mathsf{s}},r)=\frac{|r-{r}_{\mathsf{s}}|+|{r}_{\mathsf{g}}-r|}{{v}_{0}} $

(8)

式中$ \gamma $表征振幅项。

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

$ \boldsymbol{M}(r,r\mathsf{\text{'}})=\iint \boldsymbol{S}\left(\omega \right)\mathsf{e}\mathsf{x}\mathsf{p}\left[\mathsf{i}\boldsymbol{K}\right(r-r\mathsf{\text{'}}\left)\right]{\mathsf{d}}^{2}\boldsymbol{K} $

(9)

式中：$ \boldsymbol{S}\left(\omega \right) $为地震子波的频谱；$ \boldsymbol{K}=-{\boldsymbol{k}}_{\mathsf{s}}+{\boldsymbol{k}}_{\mathsf{g}}=\omega \nabla \boldsymbol{\tau }({r}_{\mathsf{g}},{r}_{\mathsf{s}},{r}^{\text{'}}) $，且$ \left|{\boldsymbol{k}}_{\mathsf{s}}\right|=\left|{\boldsymbol{k}}_{\mathsf{g}}\right|={\boldsymbol{k}}_{0} $。由此，结合地震子波信息、观测系统参数，利用射线追踪求取射线路径信息得到散射波数向量，在波数域建立点扩散算子。再作用到地下速度模型的速度扰动，即与反射系数作用，并通过傅氏逆变换即可得到深度域理论成像结果。

此处的成像算子实质上就是由地下速度模型及观测系统共同构建的滤波器，模拟成像结果则是地下真实反射率经过此滤波器得到的理论成像结果。

多算子协同模拟成像则是在模拟成像的基础上，设置合适的偏移成像窗，在相同的速度模型和观测系统的前提下构造多个控制成像的点扩散算子$ {F}_{1}\mathsf{、}{F}_{2}\mathsf{、}\cdots 、{F}_{n} $，与对应窗内的速度扰动函数$ {O}_{1}\mathsf{、}{O}_{2}\mathsf{、}\cdots 、{O}_{n} $进行空间域褶积，最终组合完成多算子模拟成像

$ \left[\begin{array}{c}{I}_{1}\\ {I}_{2}\\ ⋮\\ {I}_{n}\end{array}\right]=\left[\begin{array}{c}{F}_{1}\\ {F}_{2}\\ ⋮\\ {F}_{n}\end{array}\right]\left[\begin{array}{cccc}{R}_{1}& {R}_{2}& \cdots & {R}_{n}\end{array}\right] $

(10)

式中R为反射系数。

2 算子插值与模拟成像模型试算 2.1 不同插值方法

求取散射点处的散射波数向量需要预先求取散射点处的射线路径，本文选用算法稳定的多模板快速行进法计算走时，再使用最短路径法求取射线路径信息^[25-28]。

为提高计算效率，先对模型进行粗网格离散，建立粗点扩散算子并进行不同方法插值测试。此方法区别于常见的对地震射线追踪的走时线性插值^[29-30]，插值域也从空间域变为波数域。

2.2 简单模型模拟成像

先使用一个501×501大小的网格且夹带低速层的简单速度模型(图 2a)，使用粗网格对模型进行剖分以提高计算效率。使用1∶5比例网格抽稀原速度模型(图 2b)，求取射线路径(图 2c)、反射系数(图 2d)。观测系统参数：均匀激发20炮，每炮100道接收，地震子波为主频30 Hz的Ricker子波。点扩散算子布置在(3500 m，3500 m)处。

散射波数向量的不完备导致粗算子很多方向存在信息缺失，这将会对模拟成像结果造成影响。由图 3可见，不同的插值方法得到的算子在形态上与未插值算子保持一致，但也有明显区别。相比于粗点扩散算子(图 3a)，最近邻插值法(图 3b)得到的算子仅相当于等比例放大，依旧缺失信息。双线性插值(图 3c)与双三次插值(图 3d)不仅将原有的小波数图放大，以便与反射系数波数图在波数域相乘的大小，而且通过线性计算补充了部分新的值。双三次插值算子产生的新数值数量明显多于双线性插值算子，新值虽然可以填补算子缺失信息，但也可能对成像结果造成影响。

使用反射系数与地震子波在深度域做非稳态褶积，可以得到未经点扩散算子作用的理想情况的成像结果，以此作为标准对比不同插值方法的点扩散算子的成像效果。相比于双线性插值法(图 4c)与双三次插值法(图 4d)，最近邻插值法模拟成像结果干扰明显较多，信噪比相对较低(图 4b)。使用双线性插值法和双三次插值法得到的结果都能较好地减少噪声，对反射系数的模糊程度较低。但是由于双三次插值法使点扩散算子插值回原来大小的同时，在原本没有散射向量的区域产生了更多的数值，这些数值不一定有益于模拟成像结果，故其与双线性插值法模拟结果还是有微弱的区别。

用反射系数与单位脉冲作褶积，得到理想成像结果，选用三种插值法模拟成像结果与理想成像结果的第100、第300道地震记录进行对比(表 1)。由表可见，同一速度模型相同粗网格的情况下，三种插值方法提升的计算效率几乎相等，综合考虑了成像效果，采用双线性插值方法得到的模拟成像结果受粗网格因素的影响也最小，因此，此插值方法是较好的选择，后续对算子插值均采用双线性插值方法。

2.3 复杂速度模型粗网格模拟成像分析

使用1361×561网格的Marmousi模型进行复杂速度模型的粗网格模拟成像。地震子波为主频30 Hz的Richer子波，地面激发20炮，道间距20 m，每炮227道接收，计算(1800 m, 1800 m)处算子。分别使用横向1∶3、纵向1∶2和横向1∶6、纵向1∶4的两种粗网格速度模型求取射线路径构建粗点扩散算子。

可以由射线路径图(图 5)看出，(1800 m, 1800 m)处射线路径在原速度模型基础上求取的射线覆盖面积要大于粗网格速度模型，表明了原模型射线路径携带的照明信息更多。即原速度模型射线在理论照明程度上要优于粗网格模型。反映到点扩散算子上，即使用粗网格速度模型建立的点扩散算子(图 6b、图 6c)即使通过插值得到与原点扩散算子大小相同的算子也会缺失不少信息，而且网格越粗，信息缺失越严重。

在相同观测系统下，使用粗网格速度模型进行射线追踪虽然可以节省计算时间，但是对于散射点附近的成像效果会造成影响。计算粗网格算子(图 7c、图 7d)与原模型算子(图 7b)成像结果的相关系数分别为0.7902和0.7713，总计算效率提高了约84%和89%。综合考虑计算效率、模拟成像结果和可能存在的不确定因素，对于复杂模型也可以使用粗网格构建点扩散算子进行模拟成像。

2.4 算子插值模拟成像方法分析

对波数域粗算子进行插值，测试方法中较优的为双线性插值法。用此方法对粗点扩散算子进行插值并与反射系数进行褶积，得到的模拟成像结果干扰最少。基于粗网格构建的粗算子插值模拟成像可很大程度上提高计算效率，但是对于不同复杂程度的模型要具体问题具体分析，选择合适的粗网格在很大程度上减少计算时间的同时也不会偏离获取理论成像结果的初衷。

3 多点布置算子模拟成像 3.1 多点算子布置原则与理论模型测试

传统的模拟成像结果是通过在目的层放置散射点，根据观测系统与速度模型构建点扩散算子进行模拟成像。对于散射点附近的局部成像较为准确，却不适合控制整个地下空间的模拟成像。本文提出多算子模拟成像方法，设置多个偏移成像窗并布置多个算子，计算多点处散射向量，构建多点处成像算子还原观测系统对于地下介质成像的全局控制效果。

根据速度模型与成像精度要求设置偏移成像窗的尺寸，成像窗的长、宽比应控制在1~2范围内，尽量保证成像窗中心位置的算子对整个成像窗成像效果的控制。又因为增加新的算子会增加计算量，所以布置算子是应综合考虑计算时间与成像精度要求。

使用1361×561网格的Marmousi模型进行粗网格多算子模拟成像测试。地面均匀激发40炮，每炮227道接收，地震子波为主频30 Hz的Ricker子波，使用横向1∶3、纵向1∶2的粗网格得到一个网格为454×281的新模型。布置3×2个算子，算子位于(1125 m, 700 m)、(1125 m, 2100 m)、(2250 m, 700 m)、(2250 m, 2100 m)、(3375 m, 700 m)和(3375 m, 2100 m)。每个算子控制一个近似正方形的成像区域。为了方便表述，给这些位置赋予1号到6号的位置编号(图 8)。分别构建1~6号位置的点扩散算子(图 9)，完成多算子以及单算子模拟成像(图 10)。

通过对比不同算子情况的模拟成像结果，发现使用单点扩散算子控制全局的模拟成像得到的仅为单一散射点处照明信息，不适宜用作全局的成像算子(图 10b、图 10c、图 10d)。例如使用4号算子进行插值后模拟成像，其方位信息仅能保证对倾角较小的地层有较好的成像效果，不符合全局理论成像结果(图 10c)。同理，2号算子的方位信息较为完备，也不适宜作用到照明不好的区域进行模拟成像(图 10b)。所以最合适的是使用多散射点进行综合模拟成像。

3.2 实际数据测试

观测系统信息：地面均匀激发20炮，每炮125道接收，地震子波为30 Hz的Ricker子波。具体成像流程参考图 11。

选取网格为2001×1201的纵向速度变化剧烈且夹有一层低速异常的实际速度模型。使用射线类偏移成像中平滑处理方法^[31-33]得到不同平滑次数和平滑系数的三种不同平滑程度的速度模型。如图 12b、图 12c为两种简单平滑速度模型，图 12d则使用了纵向更大的平滑系数检测存在速度误差情况的模拟成像结果。

对四种模型使用横向1:5、纵向1:4的粗网格进行离散，算子均匀布置在4条平行于x轴的直线z=750、2250、3750、5250 m和5条平行于z轴的直线x=1000、3000、5000、7000、9000 m的20个交点的位置。

根据每个散射点处的射线路径信息可以构建控制每个成像窗的波数域点扩散算子，并在波数域进行插值(图 13)。因为速度平滑的原因，不同平滑模型的反射系数相较于原模型的反射系数也更模糊，符合实际情况中宏观的深度域地层速度模型，如图 14所示。

多算子模拟成像需要使用成像窗内反射系数与算子在波数域相乘；平滑速度模型也符合只能得到宏观的深度域层速度模型；得不到准确的速度反射系数这一实际情况。即使使用适合射线类偏移成像的平滑方法，平滑后的速度模型在射线照明区域上有优势，但是平滑后的反射系数也影响多算子成像效果。

在原模型多算子成像清晰的区域在平滑速度模型多算子成像中也能很好地模拟成像，成像结果也符合平滑后的速度模型。如图 15所示，观测系统无论是对2000 m处的低速异常体还是3500~4500 m处的速度层，其成像控制效果均较好；同样对于成像结果干扰较多的近地表浅层，四个模型的多算子成像也都能体现出观测系统对此区域成像效果的控制不佳。所以对于所提实际数据，如果针对浅层进行地震勘探，需要改进观测系统；若针对深层进行地震勘探，则可以使用此观测系统。因此，可根据目的层具体问题具体分析地设计和优化观测系统。

4 结论

本文在已有的研究基础上改进了模拟叠前深度域偏移成像方法，提出了粗网格多算子模拟成像，在粗网格的基础上使用多散射点算子综合模拟成像。理论模型和实际数据测试表明：对于不同复杂程度的模型使用不同程度粗网格，并设置合适的局部偏移成像窗，在不影响计算速度与成像效果的前提下可以得到观测系统对全局理论模拟成像结果的控制；此方法使用多算子代替单算子的模拟成像效果，在保持高计算效率的同时，更直观地展示了观测系统对不同区域理论成像结果的控制效果。