探地雷达(Ground Penetrating Radar, GPR)偏移作为一种高分辨率成像技术，是获取地下复杂结构分布的有效手段之一(Leuschen and Plumb, 2001; Streich et al., 2007; 冯德山和戴前伟, 2008; Bradford, 2015).传统的GPR单程波偏移算法如克希霍夫积分法(Moran et al., 2000; Zhuge et al., 2010)、有限差分波动方程法(冯德山等, 2011)、频率波数域偏移法(Bitri and Grandjean, 1998; Leuschen and Plumb, 2000)、相位移波动方程偏移法(Lopera et al., 2007)、分裂式傅里叶变换偏移法(Sena et al., 2006)大都将多次散射波作为噪声处理，且具有不适应横向速度变化，存在倾角限制的缺点，对复杂地质结构的偏移精度差.逆时偏移成像(Chattopadhyay and McMechan, 2008)基于波动理论的双程波算子，能够精确地处理沿任意方向传播的电磁波，从而精准地将电磁波场归位到其真实的地下空间位置，具有精度高、相位准确和不受横向变速和高陡倾角影响的优点，克服了上述单程波偏移算法所面临的成像缺陷，能刻画正确的运动学和动力学信息.因此，深入研究GPR逆时偏移成像方法，全面提高复杂目标体的成像质量，有助于提高实测资料的解释精度.

电磁波与地震波在运动学和动力学存在一定的相似性，且GPR信号和地震信号都是空间位置的时域电压信号，故GPR逆时偏移成像研究在一定程度上可借鉴地震勘探的偏移算法和理论(Chattopadhyay and McMechan, 2008; Zhu et al., 2014).近几年来，逆时偏移被广泛应用于地震资料处理的同时，许多学者将逆时偏移算法引入到GPR数据处理中：Fisher等(1992)提出了一种基于矢量波动方程的逆时算法并应用到共偏移距GPR数据处理；Leuschen和Plumb(2001)在传统逆时偏移算法的基础上，采用线性逆散射理论提出了一种基于匹配滤波技术的逆时偏移算法，改善了GPR逆时偏移成像的精度；傅磊等(2014)采用零延时互相关成像条件的逆时偏移算法实现了合成及实测的机载雷达数据逆时偏移，与传统偏移算法的结果进行对比，证明了逆时偏移算法的优越性；Liu等(2014)将全波形反演的结果作为偏移速度模型，改善了GPR逆时偏移的精度；Bradford(2015)等针对起伏地形的逆时偏移问题，提出了一种考虑电磁波振幅补偿的逆时偏移算法并应用于合成及实测GPR数据的处理中；Lu等(2016)开发了一种适合起伏地表的GPR逆时偏移成像算法，有效处理了起伏地表问题；朱尉强和黄清华(2016)考虑到电导率对电磁波衰减的影响，通过在电磁波的正演模拟及逆时传播过程中改变衰减项的正负号，保证了逆时传播的时间对称性，提出了一种GPR衰减补偿逆时偏移成像方法，数值实验验证了该方法可提升高电导率区域的成像效果.

上述GPR逆时偏移算法大都采用互相关成像条件，在互相关成像条件计算过程中，源点正传电磁波场与接收点反传电磁波场的计算顺序相反，需要存储所有时刻的源点正传电磁波场，及巨大的计算机存储空间.现有的源点电磁波场存储策略如保存边界处的电磁波场值、检波点法和随机边界条件法，需要在接收点电磁波场逆推的同时，对源点正传电磁波场进行重建，严重降低了成像效率(刘红伟等, 2010；陈婷和何兵寿, 2014).此外，由于直达波和激励源的影响，基于互相关成像条件的成像剖面存在严重的低频噪声，大大降低了成像精度(刘玉金和李振春, 2015).为此，Nguyen和McMechan(2012)在声波逆时偏移成像中提出了激发振幅成像条件，它是一种特殊的互相关成像条件，等效于在成像时刻用检波点电磁波场与单位脉冲的源点电磁波场进行互相关并进行归一化.相比于互相关成像条件，该成像条件在源点正传过程中只需保存最大能量密度时刻的源点电磁波场及相应的走时，而不需要存储任何源点正传电磁波场，且在接收点逆推过程中无需进行源点电磁波场重建，在保证具有互相关成像精度的同时，极大地提高了计算效率.目前，激发振幅成像条件被成功应用于声波、弹性波(张智等, 2013; Gu et al., 2014; Nguyen and McMechan, 2014)逆时偏移成像问题中，并取得了良好的效果.

为此，本文在上述前人工作基础上，针对电磁波逆时偏移问题将激发振幅成像条件应用于GPR逆时偏移成像中.电磁波场的正向及逆向传播采用基于三角形剖分的有限单元法(Finite element Time Domain method, FETD)进行计算，单轴各向异性完全匹配层吸收边界条件用于吸收模型截断边界处的超强反射波.应用基于激发振幅成像条件的GPR逆时偏移算法对倾斜界面模型和起伏界面模型的多偏移距数据进行逆时偏移成像，并与归一化互相关成像条件的逆时偏移结果进行对比，结果表明基于激发振幅成像条件的逆时偏移结果明显优于归一化互相关成像条件的逆时偏移结果.

GPR逆时偏移的基本原理是将地表记录到的接收点电磁波场在时间轴上进行逆向传播，当电磁波场逆推至零时刻，则所有反射波与绕射波的能量都回到最初被反射和绕射的空间位置，然后应用成像条件可获取最终的偏移成像剖面(Claerbout, 1971; 底青云等, 2000).因此，GPR逆时偏移实现过程主要包括3个环节：(1)源点电磁波的正向传播，计算从源点到各成像空间点的时间，为初值问题；(2)接收点电磁波的逆向传播，将接收点电磁波按照时间逆序插入到对应的空间位置，为边值问题；(3)成像条件的应用(有关成像条件的讨论见下节).前两个过程，可利用相同的算法进行计算，本文采用基于三角形剖分的FETD算法实现电磁波场在时空域中的正向和逆向传播.

根据电磁波理论，考虑二维地电模型的走向为y轴，则TM模式下电场y分量满足的波动方程为(戴前伟等, 2012)

(1)

式中，E_y为y方向电场强度(V·m^-1)；ε, σ, μ₀分别表示为介质的介电常数(F·m^-1)，电导率(S·m^-1)和真空的磁导率(H·m^-1)；S为激励源函数.利用伽辽金有限元法可推导式(1)的有限元方程(王洪华和戴前伟, 2014)为

(2)

式中，M，K′，K分别表示质量矩阵，阻尼矩阵和刚度矩阵，它们只与介质的物性参数和几何分布有关；S为源向量.采用三角形网格对计算区域进行剖分如图 1所示，每个矩形网格被两条对角线再剖分成四个三角形网格，以提高复杂GPR模型的正向与逆向电磁波场的计算精度(徐世浙, 1994).同时，在计算区域的外边界采用完全匹配层吸收边界条件处理截断边界处的超强反射波(Drossaert and Giannopoulos, 2007; 李静等, 2010; 王洪华和戴前伟, 2013).

图 1

Fig. 1

图 1 有限元网格剖分及节点编号示意图 Fig. 1 Sketch map of mesh division and node numbering in finite element method

逆时偏移的成像剖面受直达波和激励源等的影响，在剖面顶部会形成弧形强干扰；同时反射波在电性差异的界面会形成呈低波数特征的串扰现象.这些低频噪声(Yoon and Marfurt, 2006; Symes, 2007; Chattopadhyay and McMechan, 2008; 刘红伟等, 2010)会严重干扰GPR逆时偏移成像的精度.为高效地去除上述低频噪声，本文采用空间高通滤波(Mulder and Plessix, 2004)进行去除，其滤波系数如下：

(3)

成像条件决定着逆时偏移的成像质量和计算效率.借鉴弹性波逆时偏移成像理论，常用的GPR逆时偏移的成像条件主要有零时间成像条件(Chang and McMechan, 1986; 底青云等, 2000; 冯德山等, 2011)，零延时互相关成像条件(Chattopadhyay and McMechan, 2008; Fletcher et al., 2009)，归一化互相关成像条件(Schleicher et al., 2008; Chen and He, 2014).其中，零时间成像条件根据爆炸反射界面成像原理，在计算模型速度减半的条件下将采集的GPR信号逆推至零时刻从而得到成像剖面，目前主要用于GPR数据的叠后逆时偏移处理.本文重点讨论用于GPR叠前逆时偏移的归一化互相关成像条件和激发振幅成像条件.

Claerbout(1971)根据成像点的入射波到时与反射波到时相同提出了零延时互相关成像条件，它通过计算源点正传电磁波场和接收点反传电磁波场的零延时互相关，从而获得成像剖面.标准的零延迟互相关成像条件公式为

(4)

式中I(x, y)为地下空间的成像结果，S(x, y, t)为源点电磁波场，R(x, y, t)为接收点电磁波场，，分别表示对源点和时间进行叠加.当接收点电磁波传播至t时刻时，将t时刻的源点电磁波场和接收点电磁波场相乘，最后将所有源点和所有时刻的对应网格点相加即可得到最终成像结果.然而应用该成像条件时，GPR剖面中的直达波和背景散射波做互相关会产生低频噪声，且成像结果与地下介质的反射系数没有对应关系.鉴于零延时互相关成像条件上述不足，Schleicher(2008)等在零延迟互相关成像条件的基础上，提出了归一化互相关成像条件，其表达式为

(5)

由式(5)可知：归一化互相关成像条件是在零延时互相关成像条件的基础上，用源点电磁波场进行归一化，对成像噪声具有一定地压制效果.然而，在应用归一化互相关成像条件时，要求在逆时外推接收点电磁波场的每个时刻做成像运算时的源点正传电磁波为已知.因此，GPR数据进行逆时偏移计算之前，需要计算并保存所有时刻的源点电磁波场，这将导致计算机存储量都非常大.

激发振幅成像条件(Nguyen and McMechan, 2012)通过在源点电磁波场的正传过程中，计算出每个时刻所有网格点的能量密度，并保存相应的最大能量密度对应的走时和电磁波场值；而在接收点电磁波场的逆传过程中，利用保存的源点电磁波场做归一化，从而得到与反射系数相关的成像剖面.激发振幅成像条件表达式为

(6)

(7)

式(6)和式(7)中，R(x, z, t_k)分别表示最大能量密度对应的接收点电磁波场；t_k为最大能量密度时刻源点电磁波场的走时；S_max(x, z, t_k)为最大能量密度时刻对应的源点电磁波场；E_y(x, z, t)为接收点电磁波场.由式(6—7)可知：激发振幅成像条件只需保存所有网格点处最大能量密度时刻的源点电磁波场及相应的走时；同时在检波点电磁逆传过程中，在每个网格点提取最大能量密度时刻的检波点电磁波场，并用最大能量密度的时刻源点电磁波场进行归一化，避免了归一化互相关成像条件中源点电磁波场的存储和重建问题，在降低计算机存储量的同时，提高了计算效率.相比于归一化成像条件的成像结果为反射系数(源点电磁波场与接收点电磁波场的比值)的平方，激发振幅成像条件的成像结果为反射系数，可获得依赖反射系数的成像剖面.不同于声波传播，电磁波传播过程中某些网格点处场值始终为零或是很小，式(6)的计算会出现数值不稳定现象.为此，借鉴张智等(2013)和Gu等(2014)在弹性波逆时偏移成像中的改进方法，可得到如下形式的激发振幅成像条件：

(8)

(9)

其中，S表示最大能量密度时刻对应的源点电磁波场的绝对值后再求平均值，ND为网格点总数，sign[S(x, y)]表示对S(x, y)取符号.

由式(8)、(9)可知，在应用激发振幅成像条件进行逆时偏移之前，需要事先计算并保存所有网格点处最大能量密度的走时及相应的源点电磁波场.由于本文只利用电场分量E_y进行逆时偏移成像，因此在源点电磁波场正传过程中，所有网格点的能量密度w的计算公式可定义如下：

(10)

使用基于三角形网格剖分的FETD算法计算源点及检波点电磁波场，在计算能量密度时，需要在时空域进行插值.为提高计算效率，本文采用线性插值.在源点电磁波场正传的每个时刻，计算所有网格点的最大能量密度，并与前一时刻保存的网格点最大能量密度w(i－1)进行比较，若当前时刻的能量密度w(i)大于w(i－1)，则更新w(i－1)、并保存相应的走时t^i－1及电磁波场E_y^i－1，直到最大时刻，最终可获得每个网格点相应的最大能量密度对应的走时和源点电磁波场.

基于激发振幅成像条件的GPR逆时偏移算法流程如图 2所示，在给定偏移速度模型的基础上，第一步是采用FETD计算源点电磁波场的正传过程，激励源函数采用雷克子波，在源点电磁波场正传的每个时刻，计算各个网格点的能量密度，并与前一时刻保存的最大能量密度w(i－1)进行比较，若当前时刻的能量密度w(i)大于w(i－1)，则更新w(i－1)，并保存相应的走时t(i－1)及电磁波场E_y(i－1)直到最大时刻T_max，则可获得每个网格点的最大能量密度对应的走时和源点电磁波场值；第二步是计算接收电磁波场的逆传过程，利用保存的最大能量密度对应的走时t，提取出相应时刻的接收点电磁波场值，直到0时刻；第三步是激发振幅成像条件的应用，利用保存的源点电磁波场对提取出的接收点电磁波场做归一化，即可得到单个源点的逆时偏移剖面.此外，通过空间高通滤波技术去除低频噪声的影响，并将所有源点的逆时偏移结果进行叠加，最终得到逆时偏移成像剖面.

图 2

Fig. 2

图 2 基于激发振幅成像条件的GPR逆时偏移成像流程图 Fig. 2 RTM algorithm flow for GPR based on excitation amplitude imaging condition

依据上述方法原理编制了基于激发振幅成像条件的GPR逆时偏移成像程序，利用该程序分别对倾斜界面模型和起伏界面模型的多偏移距GPR数据进行逆时偏移处理，并与归一化互相关逆时偏移成像结果进行对比分析.

图 3为2.0 m×1.0 m的倾斜界面模型示意图，模型上部存在一个倾角9°的倾斜界面，倾斜界面上方介质的相对介电常数为ε₁=5，电导率为σ₁=0.001 S·m^-1；下方介质的相对介电常数为ε₂=10，电导率为σ₂=0.002 S·m^-1；埋深为0.65 m的位置处存在一条水平界面，其下方介质的相对介电常数为ε₃=15，电导率为σ₃=0.005 S·m^-1.利用二维FETD正演算法计算电磁波在该模型中的正向及逆向传播时，计算区域被剖分成200×400的正方形网格单元，网格间距为0.005 m×0.005 m，每个正方形网格被两条对角线再剖分成四个三角形.为消除人工截断边界处产生的强反射波，采用完全匹配层边界条件以吸收模型截断边界处的强反射波，完全匹配层厚度为0.1 m.采用共中心点观测方式，发射源信号采用中心频率为500 MHz的雷克子波，采样时间间隔为0.01 ns，时窗长度为16 ns.在地表布设11个发射源(倒三角形表示)，发射源的间距为0.1 m；第一个反射源放置在水平距离0.5 m；每个发射源两侧各布置50个接收源，接收点的间距为0.01 m，水平方向上覆盖了0.5 m.

图 3

Fig. 3

图 3 倾斜模型示意图 倒三角表示11个发射源的位置，虚线为第6个发射源的偏移孔径. Fig. 3 Schematic map of the inclined model The inverted triangles denotes the position of 11 transmitting source and the dash lines denote the migration aperture of the 6th transmitting source.

图 4为第6个发射源偏移孔径内，采用(8)式计算得到的最大能量密度时刻的电场y分量的源点电磁波场及相应的走时.第6发射源的偏移孔径如图 2中的虚线所示.其中图 4a为计算的走时等值线图，图 4b为最大能量密度时刻的电场y分量的源点电磁波等值线图；图中虚线表示两个真实界面的位置.由图 4a可知：电场y分量以标准的半同心圆向外扩散传播，经时深转换，与理论上的电磁波传播时间相一致.图 4b中电磁波能量在深度为0.2 m范围内以标准的半同心圆向外扩散传播，在0.2~0.4 m范围内0.6 m位置处电磁波能量发生变化，出现两条明显的能量突变界面，这是由于电磁波在遇到电性差异的物性分界面时会发生反射及绕射所致，这两条能量突变界面与真实模型的电性参数分界面的位置刚好吻合，验证了本文中最大能量密度时刻源点电磁波场及相应走时计算方法的正确性和有效性.

图 4

Fig. 4

图 4 倾斜模型第6个发射源偏移孔径内的最大能量密度走时等值线(a)和电场y分量的源点电磁波场(b) Fig. 4 Maximum energy traveltime contour (a) and source wavefiled of electric y component (b) of the 6th transmitting source in incline model

在获得第6个发射源偏移孔径内的最大能量密度时刻及其源点电磁波场的基础上，利用GPR逆时偏移成像算法对第6个发射源偏移孔径内接收的GPR数据进行计算，获得的成像剖面如图 5所示.其中图 5a和5b分别为基于归一化互相关成像条件和激发振幅成像条件的逆时偏移成像结果，经过高通滤波之后的成像剖面.由图中可知，成像剖面都能看到两条非常清晰的反射界面，这两条反射界面与真实模型的界面位置非常相符.对比图 5a和图 5b可知：激发振幅逆时偏移成像剖面中，反射界面的串扰更小，特别是第二条反射界面.因此，激发振幅成像条件的分辨率要优于归一化互相关成像条件，这是因为前者的偏移剖面在量上为反射系数，而后者的偏移剖面在量上为反射系数的平方.此外，空间高通滤波技术将成像剖面中直达波和地滚波切除不完全造成的低频噪声进行了有效压制.但是对比图 5a和图 5b可知：基于激发振幅成像条件的逆时偏移成像剖面中的低频噪声非常微弱，几乎不可见；而基于归一化互相关成像条件的逆时偏移成像剖面中的低频噪声还清晰可见.

图 5

Fig. 5

图 5 倾斜模型的第6个发射源偏移孔径内的成像剖面 (a)归一化互相关成像条件; (b)激发振幅成像条件. Fig. 5 Imaging profiles of the migration aperture of the 6th transmitting source in the incline model (a) Imaging profile from the normalized correlation imaging condition; (b) Imaging profile from excitation amplitude imaging condition.

图 6为倾斜模型中11个发射源偏移孔径内接收到GPR数据的逆时偏移叠加成像剖面，其中图 6a和6b分别为基于激发振幅成像条件和归一化互相关成像条件的逆时偏移成像结果.由图可见：两种成像条件获得的偏移剖面的界面形态都非常清晰，由于模型两端在偏移成像过程中叠加次数为1，而中间区域进行了多次叠加，因此界面两端的成像不如中间界面清楚，且成像剖面中都存在明显的串扰噪声.对比图 6a和图 6b可知：图 6a中成像结果的分辨率明显优于图 6b.为了更好地对比两种不同成像条件的逆时偏移成像结果的分辨率，将图 6中水平位置处1.0 m的叠加能量进行对比如图 7所示.由图可知：基于激发振幅成像条件的成像结果能量更汇聚，串扰现象更小.对比结果表明：相比于归一化互相关成像条件，激发振幅成像条件的成像结果的分辨率更高，成像质量更好.

图 6

Fig. 6

图 6 倾斜模型的11个源点叠加成像剖面 (a)激发振幅成像条件; (b)归一化互相关成像条件. Fig. 6 Stacked imaging profiles of 11 sources in the incline model (a) The normalized correlation imaging condition; (b) Excitation amplitude imaging condition.

图 7

Fig. 7

图 7 倾斜界面模型不同成像条件的逆时偏移结果在水平位置1.0 m处的叠加能量对比 Fig. 7 Comparison of stacked energy with different imaging conditions at position of 1.0 m in the incline model

图 8为1.6 m×1.2 m的起伏界面模型示意图，埋深约0.4 m处存在一条起伏界面，包含两个山峰和一个山谷地形，最大的起伏变化深度为0.2 m.起伏界面上方介质的相对介电常数为ε₁=5，电导率为σ₁=0.001 S·m^-1；下方介质的相对介电常数为ε₂=12，电导率为σ₂=0.005 S·m^-1.起伏界面的下方设置有一个正方形空洞和一个圆形空洞异常体，其中正方形空洞的中心位置为(0.68 m, 0.35 m)，边长为0.5 m；圆形空洞的圆形位置为(1.15 m, 0.68 m)，直径大小为0.5 m.利用二维FETD正演算法计算在该模型中电磁波正向及逆向传播时，计算区域被剖分成240×320的正方形网格单元，网格间距0.005 m×0.005 m，每个正方形网格内用两条对角线再剖分成四个三角形.为消除人工截断边界处产生的强反射波，采用完全匹配层边界条件以吸收模型截断边界处的强反射波减少模型边界处的超强反射，完全匹配层厚度为0.1 m.采用等偏移距和多偏移距观测方式，发射源信号采用中心频率为500 MHz的雷克子波，采样时间间隔为0.01 ns，时窗长度为20 ns.在地表布设15个发射天线，发射天线的间距为0.1 m；第一个反射天线放置在源点位置处；每个发射天线右侧布置80道接收天线，接收点的间距为0.01 m，水平方向上覆盖了0.8 m.

图 8

Fig. 8

图 8 起伏界面模型示意图 Fig. 8 Sketch map of fluctuation interface model

图 9是利用FETD法对起伏界面模型进行计算得到的共偏移距雷达剖面.由图可知在6~12 ns之间出现的双曲线同相轴是来自于起伏界面的反射波，其中起伏界面的拐点，特别是山谷的最低点和山峰的最高点产生了若干散射波，这些散射杂波会对地下目标的成像造成干扰.水平位置位于0.35 m，顶点约在13 ns出现的双曲线反射波来自正方形异常体的反射波，而水平位置位于1.15 m，顶点约为12 ns处出现的双曲线反射波则是来自圆形异常体的反射波.圆形空洞和正方形异常体产生的反射波呈双曲线，不能准确定位异常体的真实空间位置.由此可见，利用逆时偏移成像算法对GPR数据进行偏移成像，准确定位异常体的真实空间位置显的非常有必要.

图 9

Fig. 9

图 9 起伏界面模型的共偏移距GPR正演剖面 Fig. 9 Common offset simulated GPR profile of fluctuation interface model

图 10为起伏界面模型中15个发射源的逆时偏移成像叠加剖面，图 10a和10b分别为互相关成像条件和基于激发振幅成像条件的逆时偏移成像结果.由图可见，两种成像条件获得的偏移剖面的起伏界面形态都非常清晰，绕射波收敛且与真实位置非常相符；由于模型两端在偏移成像过程中叠加次数为1，而中间区域进行了多次叠加，因此界面两端的成像不如中间界面清楚，且成像剖面下方存在明显的串扰噪声.偏移剖面中圆状异常体和正方形异常体产生的绕射波收敛，并归位到其真实位置(如图中实线所示).对比图 10(a，b)可知：图 10b中的起伏界面成像结果宽度更小，界面更清晰；圆形和正方形空洞异常体的偏移结果更接近真实位置，垂直分辨率更好，偏移结果更清晰.此外，图 10a中起伏界面上方出现明显的低频噪声，而图 10b中地表附近的低频噪声更微弱.对比结果表明：相比于归一化互相关成像条件，激发振幅成像条件的成像结果的分辨率更高，成像质量更好，低频噪声更弱.

图 10

Fig. 10

图 10 起伏界面模型的逆时偏移成像剖面 (a)归一化互相关成像条件；(b)激发振幅成像条件. Fig. 10 Reverse time imaging profiles of fluctuation interface model (a) The normalized correlation imaging condition; (b) Excitation amplitude imaging condition.

为对比上述两种偏移成像条件的计算效率，分别统计了基于归一化互相关成像条件和激发振幅成像条件的逆时偏移成像算法对倾斜界面模型和起伏界面模型进行成像的计算时间和对硬盘的需求量，如表 1和表 2所示.逆时偏移计算过程在Lenovo ThinkPad T450s(Corel i5, 4.0 GHz)上测试运行，源点电磁波场和检波点电磁波场的计算采用FETD串行程序.归一化互相关成像条件由于需要存储电磁波场快照，其硬盘存储空间和I/O吞吐操作是巨大的，基于归一化互相关成像条件的逆时偏移成像算法对倾斜界面模型的硬盘需求达到10.6 GB，CPU耗时为55 min，平均每个反射源的偏移计算时间为5 min；而激发振幅成像条件仅需在每个时间步上计算所有网格点的能量密度，并保存最大能量密度走时的源点电磁波场及相应的走时，其硬盘需求可忽略不计，CPU耗时为29 min，平均每个发射源的偏移计算时间为2.6 min.此外，对于起伏界面模型的逆时偏移成像，基于归一化互相关成像条件的逆时偏移硬盘需求达到15.8 GB，CPU耗时为85 min，平均每个反射源的偏移计算时间为5.3 min；基于激发振幅成像条件的逆时偏移成像的硬盘需求忽略不计，CPU耗时为44 min，平均每个发射源的偏移计算时间为3.0 min.因此，从计算效率上看，基于激发振幅成像条件的逆时偏移的计算效率比归一化互相关成像条件的逆时偏移成像提高了2倍多，且硬盘存储量可忽略不计.

表 1 Table 1

表 1 逆时偏移成像的CPU耗时统计 Table 1 The statistics for CPU time consuming of reverse time migration imaging 成像条件 归一化互相关成像条件 激发振幅成像条件 倾斜界面模型 CPU耗时(min) 55 29 发射源点数 11 11 平均每个发射源耗时(min) 5.0 2.6 起伏界面模型 CPU耗时(min) 85 44 发射源点数 15 15 平均每个发射源耗时(min) 7.0 3.0

表 1 逆时偏移成像的CPU耗时统计 Table 1 The statistics for CPU time consuming of reverse time migration imaging

表 2 Table 2

表 2 逆时偏移成像的硬盘需求情况统计(单位：GB) Table 2 The statistics for hard disk occupation of reverse time migration imaging 成像条件 倾斜界面模型 起伏界面模型 归一化互相关成像条件 10.6 15.8 激发振幅成像条件 0 0

表 2 逆时偏移成像的硬盘需求情况统计(单位：GB) Table 2 The statistics for hard disk occupation of reverse time migration imaging

从逆时偏移原理出发，成功将激发振幅成像条件引入到GPR逆时偏移成像中，避免了源点电磁波场的存储和重建，为探地雷达逆时偏移提供了一种新的偏移手段.为提高复杂GPR模型的逆时偏移计算精度和效率，采用基于三角形剖分的FETD模拟电磁波场的正向及逆向传播，激发振幅成像条件用于获取逆时偏移成像剖面.两个典型模型的逆时偏移成像结果表明：相比于基于归一化互相关成像条件的逆时偏移，基于激发振幅成像条件的逆时偏移在成像精度相当的前提下大大节省了存储量，计算效率提高了近2倍，同时偏移剖面的空间分辨率更高，低频噪声更弱.