0 引言 逆时成像技术在现行高精度地震成像技术系列中，是理论较为成熟、成像最为精确的技术之一。与克希霍夫叠前深度成像技术和单程波地震深度成像技术相比，逆时成像技术具有以下诸多优点[1-4]：①对地震波动方程的近似较少；②能够解决多值走时问题；③适合于任意陡倾角及速度在横纵向变化较为剧烈等情况下的成像问题；④可实现多次波、回转反射波等通常被认为是干扰波类型波场的准确成像。然而，该技术自20世纪80年代被提出以来，并没有得到充分重视和广泛应用，主要原因是逆时成像技术存在庞大的计算量和巨大的存储量等瓶颈技术问题。近几年来，随着CPU/GPU高性能协同并行计算和大容量磁盘的并行快速存储等技术的飞速发展，较大程度上缓解了上述问题，从而改善了逆时成像技术的工业化应用现状。虽然克希霍夫叠前深度成像技术仍是当前工业界主流的地震成像技术，但其正逐渐被高精度全波场逆时成像技术所代替[5-7]。 目前，业内对地震波叠前逆时成像技术已开展了诸多研究[8-10]，主要包括叠前道集保真预处理、高精度速度建模方法（复杂构造层析反演、全波形反演等）、波动方程数值计算理论（规则网格有限差分法、交错网格有限差分法等数值离散算法，旁轴近似吸收边界、阻尼吸收边界、PML吸收边界条件等边界条件）、改进的逆时成像条件（互相关法成像条件、归一化互相关成像条件、时空移相关成像条件等）优化、逆时成像后数据体的去噪处理（拉普拉斯去噪、低通滤波去噪、扩散滤波去噪等）、GPU加速计算策略等[11-13]。这些理论方法的研究在一定程度上推动了高精度逆时成像技术的实用化进程，因此，逆时偏移技术能够在滨里海、墨西哥湾和北海等地区的盐下或盐丘侧翼等成像方面取得重大突破。虽然该技术在国内受陆地地震资料近地表等因素的综合影响[14-17]，应用进程相对缓慢，但目前也已取得了较为丰富的研究成果，如高陡构造区域成像[18-19]、复杂山前带成像[20]、起伏地表成像[21-22]和碳酸盐岩缝洞储层成像等[23-25]。 目前常用的叠前逆时成像技术普遍采用基于常规相关法的逆时成像条件，但因该成像条件在偏移成像路径上将不相关的炮、检点波场进行了互相关处理，从而引入了能量较强的背景噪声，且这种噪声在后续处理中较难有效地实现保幅压制[26]；此外，常规逆时成像技术没有考虑震源子波对不同成像结果的影响，也没有对逆时成像的关键参数进行优化。本次研究在前人研究和实践的基础上，从逆时成像技术的基本原理出发，分析了影响成像精度的关键因素，以期改进逆时成像核心算法，提高成像精度和效率。 1 逆时偏移基本原理 以各向同性介质三维地震波动方程[8]为例，其数学表达式为 $ \frac{{{\partial ^2}P}}{{\partial {t^2}}} = {V^2}\left\{ {\frac{{{\partial ^2}P}}{{\partial {x^2}}} + \frac{{{\partial ^2}P}}{{\partial {y^2}}} + \frac{{{\partial ^2}P}}{{\partial {z^2}}}} \right\} $ (1) 式中：P为地震波场，无量纲；V为声波在介质中的传播速度，m/s；t为时间，s；x，y，z分别代表空间的3个方向坐标。 在地震波场数值模拟中，受数值计算稳定性条件限制，时间步长通常较小，因此，由时间离散引起的数值误差也较小，时间导数项可采用2阶精度有限差分法进行数值离散。空间导数项主要受计算机内存等因素的限制，只能用有限空间解决无限空间的问题，因此需要采用高阶有限差分来最大程度减小空间导数与解析解之间的误差。考虑到在保证计算精度的前提下高阶空间导数可以允许更大的空间步长，减小计算量，因此，可以建立时间二阶、空间十六阶精度的有限差分算法。根据差分近似的各向异性分析[27]，十六阶差分可以有效提高数值模拟精度，满足应用需求，并能实现GPU/CPU协同加速处理。 逆时偏移包含2个主要步骤：用小时间样值推算大时间样值的波场正向传播过程，用大时间样值推算小时间样值的波场逆时延拓过程，具体公式为 $ {P^{n + 1}} = 2{P^n} - {P^{n - 1}} + Q + {f_{{\rm{source}}}}({x_0},{y_0},{z_0},t) $ (2) $ {P^{n - 1}} = 2{P^n} - {P^{n + 1}} + Q + {f_{{\rm{record}}}}({x_1},{y_1},{z_1},t) $ (3) 式中：x0，y0，z0为震源的三维空间坐标；x1，y1，z1为检波器的三维空间坐标；t为时间，s；n = 0，1，2，…，N；N为总采样点数；Pn+1，Pn，Pn-1分别为下一时刻、当前时刻和上一时刻的波场；fsource为震源函数；frecord为地震记录。 $ \begin{array}{l} Q = {\left( {\frac{{V\Delta t}}{{\Delta x}}} \right)^2}[ - 3.054\;844\;1 \times P_0^n + 1.777\;778 \times (P_{x,1}^n + P_{x, - 1}^n) - \\ \quad 0.311\;111\;11 \times (P_{x,2}^n + P_{x, - 2}^n) + 0.075\;420\;875 \times \left( {P_{x,3}^n + P_{x, - 3}^n} \right) - \\ \quad 0.017\;676\;768\left( {P_{x,4}^n + P_{x, - 4}^n} \right) + 0.003\;480\;963\;5 \times \left( {P_{x,5}^n + P_{x, - 5}^n} \right) - \\ \quad 0.000\;518\;000\;52 \times \left( {P_{x,6}^n + P_{x, - 6}^n} \right) - 0.000\;050\;742\;908\left( {P_{x,7}^n + P_{x, - 7}^n} \right) - \\ \quad 0.000\;002\;428\;127 \times \left( {P_{x,8}^n + P_{x, - 8}^n} \right)] + \\ \quad {\left( {\frac{{V\Delta t}}{{\Delta y}}} \right)^2}[ - 3.054\;844\;1 \times P_0^n + 1.777\;778 \times \left( {P_{y,1}^n + P_{y, - 1}^n} \right) - \\ \quad 0.311\;111\;11 \times \left( {P_{y,2}^n + P_{y, - 2}^n} \right) + 0.075\;420\;875 \times \left( {P_{y,3}^n + P_{y, - 3}^n} \right) - \\ \quad 0.017\;676\;768\left( {P_{y,4}^n + P_{y, - 4}^n} \right) + 0.003\;480\;963\;5 \times \left( {P_{y,5}^n + P_{y, - 5}^n} \right) - \\ \quad 0.000\;518\;000\;52 \times \left( {P_{y,6}^n + P_{y, - 6}^n} \right) - 0.000\;050\;742\;908\left( {P_{y,7}^n + P_{y, - 7}^n} \right)\\ \quad - 0.000\;002\;428\;127 \times \left( {P_{y,8}^n + P_{y, - 8}^n} \right)] + \\ \quad {\left( {\frac{{V\Delta t}}{{\Delta z}}} \right)^2}[ - 3.054\;844\;1 \times P_0^n + 1.777\;778 \times \left( {P_{z,1}^n + P_{z, - 1}^n} \right) - \\ \quad 0.311\;111\;11 \times \left( {P_{z,2}^n + P_{z, - 2}^n} \right) + 0.075\;420\;875 \times \left( {P_{z,3}^n + P_{z, - 3}^n} \right) - \\ \quad 0.017\;676\;768\left( {P_{z,4}^n + P_{z, - 4}^n} \right) + 0.003\;480\;963\;5 \times \left( {P_{z,5}^n + P_{z, - 5}^n} \right) - \\ \quad 0.000\;518\;000\;52 \times \left( {P_{z,6}^n + P_{z, - 6}^n} \right) - 0.000\;050\;742\;908\left( {P_{z,7}^n + P_{z, - 7}^n} \right) - \\ \quad 0.000\;002\;428\;127 \times \left( {P_{z,8}^n + P_{z, - 8}^n} \right)] \end{array} $ (4) 式中：Δt为时间步长，s；Δx，Δy，Δz分别为空间x，y，z等3个方向的网格大小，m；Pi中的下标i（i =-8，-7，…，7，8）为空间样点偏离差分中心点的距离，m；P0为差分中心点的波场值。 用有限的数值空间来模拟无限的空间问题时，必然引入人工截断边界的问题。本次研究选用褶积PML吸收边界[28]（Convolutional PML）来解决人工截断边界的难题。与常规一阶声波偏微分方程不同，各向同性介质三维地震波动方程[参见式（1）]中的二阶导数需要2次应用褶积PML吸收边界条件，即通过2次复数域下的拉伸处理，将公式 $\bar x = x + \frac{1}{{j\omega + \alpha }}\int_0^x {\sigma \left(s \right){\rm{d}}s} $ 变换为 $\frac{\partial }{{\partial x}} \Rightarrow \frac{\partial }{{\partial \tilde x}} = \left({1 + \frac{{\sigma \left(x \right)}}{{j\omega + \alpha }}} \right)\frac{\partial }{{\partial x}}$ ，再变换回时间域。以x方向二阶导数为例，其具体推导过程如下： $ \begin{array}{l} \frac{{{\partial ^2}}}{{\partial {{\tilde x}^2}}}P = \frac{\partial }{{\partial \tilde x}}\left( {\frac{\partial }{{\partial \tilde x}}P} \right) = \frac{\partial }{{\partial \tilde x}}\left( {\frac{\partial }{{\partial x}}P + \psi } \right) = \frac{\partial }{{\partial x}}\left( {\frac{\partial }{{\partial x}}P + \psi } \right) + \\ \quad \quad \frac{{{\partial ^2}}}{{\partial {x^2}}}P + \frac{\partial }{{\partial x}} + \psi + \xi , \end{array} $ $ \begin{array}{c} {\psi ^n} = a{\psi ^{n - 1}} + b{\left( {\frac{\partial }{{\partial x}}P} \right)^n},{\xi ^n} = a{\xi ^{n - 1}} + b{\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}}P + \frac{\partial }{{\partial x}} + \psi } \right)^n},\\ a = \exp \left\{ { - \left[ {\sigma \left( x \right) + \alpha } \right]\Delta t} \right\},b\frac{{\sigma \left( x \right)}}{{\sigma \left( x \right) + \alpha }}\left( {a - 1} \right),\\ \sigma \left( x \right) = {c_1}c\;{d^3}\left( x \right) \end{array} $ (5) 式中：d(x)为PML层厚度（网格点数）；α，c为常数；c1为用户可调参数（经测试，c1一般取1~8的实数，其数值越大，吸收效果越好；其次，α = 0时，也可以取得较好的边界吸收效果）。 2 关键因素分析 为综合提高逆时成像精度，满足大庆探区实际生产需求[29]，结合数值模拟及实例应用，系统分析有限差分阶数、吸收边界条件、逆时成像条件、子波相位、背景噪声、逆时成像参数等关键因素对逆时成像效果的影响。根据分析结果，对这些影响因素逐一进行压制、校正或参数优化。 2.1 有限差分阶数 以均匀介质波场快照为例，用差分近似的各向异性分析理论[27]研究有限差分阶数对波场数值模拟精度的影响。图 1(a)~(d)的差分阶数依次为二阶、四阶、八阶和十六阶。从图 1可以看出，差分阶数为二时，数值频散现象最为严重；当差分阶数提高到四阶时，数值频散现象有所减弱；当差分阶数提高到八阶时，数值频散现象得到有效压制；当差分阶数提高到十六阶时，波场快照的数值计算精度最高，不同频率成分的相速度基本一致，初至波能量更加集中。由此可见，在计算网格、频率等参数一定情况下，差分阶数越高，成像精度也越高。 图 1(e)和图 1(f)分别为二阶差分和十六阶差分情况下Marmousi模型[4]的逆时偏移成像结果。采用二阶差分精度时，由于不同频率成分波场传播的相速度不一致，导致波场正演和逆时延拓过程中均有较强的数值频散能量被引入，这些能量导致地震资料信噪比降低、复杂断块区域构造特征模糊；采用十六阶差分精度时，引入的数值频散噪声能量得到了有效压制，信噪比得以提高。由此可见，通过提高有限差分阶数，可以有效提高逆时成像结果的可信度和准确性。具体差分阶数可根据逆时成像处理的网格、速度以及GPU程序设计的难易程度等因素综合确定。 下载eps/tif图 图 1 不同差分阶数情况下的波场快照和逆时成像结果 Fig. 1 Snapshots and reverse-time migration results with different difference orders 2.2 吸收边界条件 以层状介质波场快照为例，分析PML厚度对波场数值模拟精度的影响。图 2(a)和图 2(b)的PML吸收层厚度分别为3个和10个网格点。从图中可以看出，PML厚度较小时，其边界反射能量较强，与有效波场能量混叠；PML厚度增加到一定程度时，可以有效削弱或消除边界反射，提高计算结果的精度。 下载eps/tif图 图 2 不同PML厚度情况下波场快照和逆时成像结果 Fig. 2 Snapshots and reverse-time migration results with different PML thickness 图 2(c)和图 2(d)是PML吸收层厚度分别为3个网格点和10个网格点情况下Marmousi模型[4]的逆时偏移结果。从图中可以看出，PML厚度较小时，由于波场正演和逆时延拓过程中均增加了由人为截断边界而引入的虚假反射波场，因此有效模拟区域的逆时成像效果受到严重干扰；PML厚度较大时，人为截断边界引入的噪声能量可以得到有效压制，从而提高信噪比。 由此可见，通过增加PML厚度，可以有效提高逆时成像结果的可信度和准确性。PML厚度所占网格点数越多，所占用的计算时间就越长，因此，PML的具体厚度可根据逆时成像精度的要求、GPU程序设计的难易程度、生产任务的周期等因素综合确定。 2.3 逆时成像条件 以层状介质逆时成像脉冲响应为例，分析不同逆时成像条件对逆时成像精度的影响[26]。图 3(a)和图 3(b)分别为采用常规炮检点相关法逆时成像条件和基于行波分离炮检点相关法逆时成像条件的脉冲响应。从图 3(a)可以看出，经过速度界面的脉冲响应存在4条成像曲线，自上而下依次为：以炮点和检波点为焦点的椭圆曲线1，以炮点投影和检波点投影为焦点的椭圆曲线2，以炮点和检波点投影为焦点的椭圆曲线3以及以炮点投影和检波点为焦点的椭圆曲线4。其中，曲线1是实现逆时成像的主要能量曲线，其他3条曲线对逆时成像的贡献很小，而曲线3以及曲线4是引入逆时成像背景噪声的主要能量曲线。从图 3(b)可以看出，引入基于行波分离的逆时成像条件后，曲线1以及曲线2得到有效保留，曲线3以及曲线4得到有效压制。 下载eps/tif图 图 3 不同逆时成像条件下脉冲响应和成像结果 Fig. 3 Impulse response and reverse-time migration results with different reverse-time imaging conditions 图 3(c)和图 3(d)分别为采用常规炮检点相关法逆时成像条件下和基于行波分离炮检点相关法逆时成像条件下的逆时成像结果。从图 3(c)可以看出，逆时成像结果中低频背景噪声的能量较强，同相轴较为粗大，同时部分同相轴能量被掩盖；从图 3(d)可以看出，采用基于行波分离炮检点相关法逆时成像条件时，成像结果中的低频背景噪声得到有效压制，信噪比得到提高，地层细节被刻画得更加清晰。由此可见，通过分析逆时成像过程中背景噪声的成因机理，并通过优化逆时成像条件对背景噪声进行有针对性的压制，可以有效提高逆时成像的精度。 2.4 子波相位 以某地区实际地震资料为例，分析子波相位校正前后逆时成像结果与克希霍夫叠前深度成像结果之间存在的差异。图 4(a)~(c)依次为常规逆时成像结果、克希霍夫叠前深度成像结果和子波相位校正后的逆时偏移成像结果。从图 4(a)和图 4(b)可以看出，常规逆时成像结果与克希霍夫叠前深度成像结果在波形、能量和深度上均存在一定的差异，这是因为常规逆时成像波动方程通常采用一阶或二阶偏微分进行运算，这2种方程在震源加载时引入了90°的相位差（这里采用的是二阶偏微分方程）。从图 4(b)和图 4(c)可以看出，通过对常规逆时成像结果[图 4(a)]进行90°相移（子波相位校正）便可以有效消除常规逆时成像结果与克希霍夫叠前深度成像结果的差异性，从而保证逆时成像结果的有效性。 下载eps/tif图 图 4 子波相位校正前后逆时成像结果与常规克希霍夫叠前深度成像结果对比 Fig. 4 Results of conventional reverse-time imaging, Kirchhoff prestack depth imaging and reverse-time imaging in wavelet phase correction 2.5 背景噪声 成像过程中产生低频噪声是地震波叠前逆时成像方法区别于其他成像方法的显著特征之一，这些低频噪声由震源波场和检波点波场中不相关的波场沿着传播路径通过相关成像而形成，并以较强背景噪声的形式叠加在逆时成像剖面上[26]。本次研究采用一种非线性扩散滤波迭代去噪方法[12]来实现逆时成像中背景噪声与有效地层成像信号的相对分离。图 5(a)为含有背景噪声的Marmousi模型[4]逆时成像结果，图 5(b)和图 5(c)分别为压制掉的背景噪声和压制背景噪声后的逆时成像结果。从图 5可以看出，逆时成像的背景噪声具有波数低、全局平滑等特点，有效的地层反射信号被掩盖在这种背景噪声之下，通过扩散滤波，这种噪声得到有效压制[图 5(d)]，掩盖在背景噪声能量之下的有效反射信号（地层细节）被恢复出来，逆时成像精度得到提高[图 5(c)]。 下载eps/tif图 图 5 背景噪声压制前后逆时成像结果对比 Fig. 5 Reverse-time migration results before and after background noise suppression 2.6 逆时成像参数 以多层介质逆时成像脉冲响应[25]为例，分析逆时成像参数优化策略问题。图 6(a)为二维层状介质速度模型；图 6(b)为纵横向均满足频散关系情况下的脉冲响应，其横向偏移网格和纵向偏移步长均为10 m；图 6(c)为横向不满足频散关系情况下的脉冲响应，其横向偏移网格为20 m，纵向偏移步长为10 m；图 6(d)为纵向不满足频散关系情况下的脉冲响应，其横向偏移网格为10 m，纵向偏移步长为20 m。根据数值频散关系，逆时偏移网格大小为10 m时，允许的最大偏移频率为80 Hz，大于实际偏移频率（60 Hz）。从图中可看出，图 6(b)不存在由数值频散引入的偏移干扰能量，计算结果准确可靠；当横向偏移网格不满足数值频散关系时[图 6(c)]，由数值频散引起的偏移干扰仅出现在浅层局部位置，而不往纵向深层传递；当纵向偏移步长不满足数值频散关系时[图 6(d)]，由数值频散引起的偏移干扰在横向及纵向的深层均进行了传递。由此可见，逆时成像的纵向偏移步长必须满足数值频散关系，而横向偏移网格可适当放宽数值频散关系的约束，从而保证在提高逆时成像频率的同时提高逆时成像精度。 下载eps/tif图 图 6 基于脉冲响应分析的逆时成像参数优化 Fig. 6 Reverse-time migration parameter optimization based on impulse response analysis 3 实际资料成像效果分析 为了进一步说明逆时成像技术关键问题经过优化后的成像效果，以大庆探区2个实际陆地地震资料为例进行逆时成像的应用效果分析。图 7(a)和图 7(b)为大庆探区S工区分别采用常规相关法逆时成像技术和改进的高精度逆时成像技术的应用结果。从图 7可以看出，虽然2种逆时成像技术均能实现火山岩地层的有效成像，但改进的高精度逆时成像技术的应用结果在火山岩内部地层反射特征、边缘刻画以及地层连续性等细节方面成像效果更好。 下载eps/tif图 图 7 大庆探区S工区不同逆时偏移方法成像效果对比 Fig. 7 Different reverse-time migration results comparison in S project in Daqing exploration area 图 8(a)和图 8(b)为大庆探区H工区分别采用克希霍夫叠前深度成像技术和改进的高精度逆时成像技术的成像结果。从图 8可以看出，虽然2种地震成像技术均能实现复杂构造及地层的有效成像，但克希霍夫叠前深度成像技术的成像结果在复杂构造区存在较强的划弧噪声，地层细节被掩盖，而改进的高精度逆时成像技术的成像结果在相应位置处的地层成像清晰，边界刻画清楚。 下载eps/tif图 图 8 大庆探区H工区不同成像方法应用效果对比 Fig. 8 Different migration results comparison in H project in Daqing exploration area 大庆探区2个不同工区的应用效果表明，改进的高精度逆时成像技术在复杂构造和复杂波场成像条件下在复杂构造细节刻画方面展示出比较明显的应用优势，虽然其计算效率受制于现有的GPU集群、内存及并行存储能力等因素，但在构造刻画理论研究方面，该技术已经能够超过现有克希霍夫叠前深度成像的计算效果，而对于更精细的岩性刻画方面，由于逆时成像技术效率受限于具体的计算网格大小和GPU协同计算能力等因素，因此仍有较大的提升空间。 4 结论 （1）从逆时成像技术的基本原理出发，推导构建了十六阶有限差分算法及褶积PML吸收边界条件，形成了改进的高精度逆时成像技术的核心。 （2）基于脉冲响应和计算模型，系统研究了影响逆时成像技术成像精度和效率的7项因素，并分别给出了具体的优化策略，综合提高了逆时成像技术的算法精度和整体性能。 （3）改进的高精度逆时成像技术在大庆探区2个典型工区的实际地震资料复杂构造和火成岩成像处理中实现了有效应用，取得了较为明显的应用效果，地层细节刻画优于常规相关法逆时成像结果和克希霍夫叠前深度成像结果。