偏移是一种地震成像技术，它在空间域或时间域将记录到的同向轴归位到真实的反射界面.传统的地震偏移成像方法基于单程波方程，并将多次散射波当做噪声处理，但对于复杂地质结构及大倾度界面如裂缝，合理的利用多次散射信息至关重要.传统波动方程偏移方法如Kirchhoff积分偏移(Gray，1986; Gray and May，1994; Docherty)难以适应横向速度变化；有限差分波动方程偏移(Vidale，1988; Mufti et al.，1996)近似求解波动方程，存在倾角限制；频率波数域偏移方法(Loewenthal and Mufti，1983)是基于常速度模型，当地下介质速度存在变化时会出现明显的精度问题.有别于传统的基于单程波方程的偏移方法，逆时偏移(Baysal et al.，1983; Levin，1984; Chang and McMechan，1986，1994; Sun and McMechanz，2001; Yoon et al.，2003; 刘欣欣等，2013; 王保利等，2012)基于全波方程，能够精确的处理沿任意方向传播的波场，包括耦合波及横纵波之间的转换波，从而可更精准地将波场归位到其真实的地下空间位置.虽然逆时偏移技术的成功应用实例并不多见，但它在勘探地震中受到越来越多的关注，逆时偏移包括叠前逆时偏移和叠后逆时偏移.

探地雷达电磁波与地震波在运动学和动力学存在一定的相似性，此外，采集到的探地雷达数据和地震数据都是空间位置的时域电压信号，故而在一定程度上勘探地震中的数据处理技术可以运用到探地雷达数据处理中.Fisher(1992)提出用基于张量波动方程的逆时偏移算法对单通道的共偏移距的探地雷达数据进行了偏移处理，处理后的结果在一定程度显示了该算法的优越性.Feng和Sato(2004)利用叠前逆时偏移算法，对地雷探测中采集到的多偏移距数据进行处理，一方面强散射波得到压制，另一方面获得了地下更为精细的雷达剖面图像.Zhou等(2005)应用一种基于遗传算法的偏移速度分析方法，对合成及野外实测的共炮点雷达数据进行分析获得了地下介电常数分布，并以其作为叠前逆时偏移的初始速度模型，取得了较好的结果.

本文开发了电磁波逆时偏移算法，电磁波场的正向及逆向传播是基于时间域有限差分方法，单轴各向异性完全匹配层(UPML)(Fang and Ying，2009; Feng and Dai，2011)吸收边界条件用于吸收边界的反射波，互相关成像条件用于获得空间图像.应用电磁波全波形反演算法，对合成及实测的共偏移距、多偏移距探地雷达数据全波形反演，获得了与实际情况更为接近的地下介质分布情况.通过全波形反演获得的介电常数和电导率模型作为逆时偏移的初始模型，逆时偏移的结果与基于共偏移距雷达数据Kirchhoff偏移成像结果对比发现，逆时偏移成像效果明显好于传统偏移成像结果.

传统的偏移算法基本是基于单程波方程深度外推，所以在对多分量地震资料处理过程中首先要进行纵、横波分离，然后用现有的偏移算法分别处理纵波场和横波场.它的缺陷在于：1)将多分量地震资料当作几个标量波的简单叠加，忽略了地震波的矢量特性，从而影响地震多分量资料的处理与解释精度；2)假如地震波场横波与纵波没有完全分离(目前这种现象广泛存在)，而在实际资料的偏移处理过程中又假设为标量波场，这必然会使偏移结果产生许多人为假象，影响处理精度.逆时偏移与传统的偏移方法有很大的不同，逆时偏移可以看成是波场在时间轴上逆向传播过程，当波场倒退到零时刻，则所有反射波与绕射波的能量都回到最初被反射和绕射的空间位置，应用成像条件即可得到最终的偏移成像剖面.逆时偏移的波动理论方法基于全波方程，波场可以沿各个方向传播，因此该方法不受介质倾角限制，适用于速度任意变化的模型.

逆时偏移算法的实现概括起来主要包括三个步骤：第一是震源波场正向延拓；第二是接收波场的逆时延拓；第三是应用成像条件成像.前两个步骤的实现可以基于相同的数值方法，本文应用时间域有限差分法(刘四新和曾昭发，2007; 刘四新等，2006; Taflove and Hagness，2005; Liu and Sato，2005; Liu et al.，2007; Liu et al.，2011)实现电磁波场在空间的正向和逆向传播.影响叠前逆时偏移成像效果的一个非常重要因素是成像条件(Chattopadhyay and McMechan，2008)的应用，Claerbout(1971)提出的互相关成像条件是最基础的成像条件，它采用零延迟互相关成像条件，该成像条件的不足之处在于对首波及背景散射波做互相关会产生低频噪声，其次应用该成像条件得到的成像结果与地下介质的反射系数没有对应关系.鉴于零延迟互相关成像条件的弊端，后续逐步发展了反褶积成像条件(Guitton et al.，2007)、震源补偿成像条件(Guitton et al.，2006)、扩展成像条件(Sava and Fomel，2006)以及角度域成像条件(Yoon and Marfurt，2006)等.标准的逆时偏移采用零延迟互相关成像条件，其计算公式为：

本文开发的电磁波逆时偏移成像算法如图 1所示，对于合成数据，首先给定真实模型，通过FDTD模拟电磁波在模型中的传播过程，同时记录接收波场.逆时偏移成像的第一步是通过FDTD进行激励源波场的正向传播模拟，使用雷克子波作为激励源函数，电磁波最大传播时间序列为n_max，为了节约存储空间，将每个时刻模型空间的外围波场存储；第二步是电磁波波场的逆向传播模拟，该步骤包括两个方面的反向传播过程：(一)利用包围模型空间外围全时序波场恢复电磁波场的正向传播过程，这样做的目的其实就是以时间换空间；(二)利用接收波场进行电磁波场的逆向传播模拟；第三步是成像条件的应用，因为互相关条件容易实现且无稳定性问题，本文选取互相关作为成像条件.波场逆向传播的过程中，利用成像条件获得每个时刻的偏移图像，当波场传播到时间零点，将各个时刻的图像叠加即可获得单源逆时偏移成像结果.此外，通过去噪技术去除低频影响，再将所有单源偏移结果叠加得到最终的逆时偏移成像结果.

图 1

Fig. 1

图 1 电磁波逆时偏移成像算法流程Fig. 1 RTM algorithm flow for electromagnetic wave

文中构建了一个复杂地质结构模型如图 2所示，X方向长度为1.2 m，Z方向长度为0.7 m.图中蓝色区域代表空气，暗红色区域代表土壤层，其相对节电常数为6，电导率为5 mS·m^-1，平均厚度约为0.12 m，该土壤层的底界面起伏变化；黄色区域代表干沙层，其相对介电常数为4，电导率为1 mS·m^-1，平均厚度约为0.48 m.由于土壤层与干沙层之间存在大倾角起伏界面，故而该模型中存在变化剧烈的横向速度.干沙中埋有2个异常体，其中蓝色方形区域代表一个高速异常体，其相对介电常数为1；圆心位于(0.8，0.4)m处圆形暗红色区域代表一个低速异常体，其相对介电常数为6，电导率为4 mS·m^-1.

图 2

Fig. 2

图 2 相对介电常数模型Fig. 2 Relative dielectric constant model

典型的雷达反射测量中，电场与测线往往是垂直的，因此我们应用二维TM模FDTD求解Maxwell方程组，同时采用UPML吸收边界条件来减少边界处反射.为确保稳定性及减少数值频散，设置网格大小为0.004 m，发射源为主频1600 MHz的雷克子波，时间取样间隔0.008 ns，采样点为1024.首先合成了共偏移距剖面探地雷达数据，第一道信号发射天线的水平距离位于0.05 m，接收天线与发射天线间距为0.05 m，发射天线的移动步长为0.012 m，沿地面总该采集了88道雷达信号，水平方向上覆盖了1.05 m的水平距离.图 3a是合成的共偏移距雷达剖面，横坐标代表道号，纵坐标代表雷达波的双程走时，从图可知在2.5~4.5 ns之间的双曲线同向轴是来自土壤层和沙层界面的反射波，值得注意的是出现了来自不规则界面拐点处的若干散射波，这些散射杂波会对地下目标的识别造成困扰.水平位置位于第20道，顶点约在6 ns出现的双曲线同向轴来自方形高速异常体的反射波，而水平位置位于第60道，顶点约为4.8 ns处出现的双曲线同向轴则是来自圆形低速异常体的反射波.合成的共炮点雷达数据如图 3b所示，水平轴代表道号，纵轴是雷达波的双程走时.共炮点雷达数据集与共偏移距雷达剖面在某种程度上存在一定的相似性，来自界面及异常的反射波出现的双程走时相当，此外反射波的曲线形态也存在一定的相似性.

图 3

Fig. 3

图 3 (a)共偏移距雷达剖面，88道信号在X方向覆盖了约1.05 m，(b)雷达数据集包含17条共炮点雷达剖面，所有图像进行了幅度归一化.Fig. 3 (a) The common-offset radar profile, 88 traces signal in total covering about 1.05 m in X direction, (b) The radar gather consisting of 17 common-shot profiles, both figures are normalized by the maximum amplitude.

逆时偏移成像算法中，影响最终成像效果的一个最重要因素是初始模型.为了获得一个较高精确度的初始模型，本文采用时间域全波形反演技术，反演地下介质的相对介电常数分布，从而为逆时偏移成像提供更为精准的初始模型.全波形反演技术试图利用包括走时、幅度、相位等全波形信息来反演地下参数分布，它比常规的基于射线的成像方法能提供更加可靠的参数信息，本文应用的时间域全波形反演技术细节可参考(Meles et al.，2010，2011)，其主要目的是建立包含观测数据和模拟数据的目标函数(2)式：

时间域全波形反演算法中，不管是单参数分别反演还是介电常数和电导率联合反演，初始介电常数和电导率模型的选取至关重要.考虑到如图 2所示的真实相电常数模型中存在起伏较大的界面，地下介电常数模型不易给定，故而这里选取的初始介电常数模型如图 4a所示，蓝色区域代表空气；黄色区域代表地下介质，其相对介电常数为4.5.图 4b是时间域全波形反演25次迭代反演后得到的相对介电常数结果，首先土壤与沙层之间的起伏界面得到了清晰的刻画，起伏界面的形状及位置与真实模型匹配较好，另外在反演的界面上方相对介电常数要高于界面下方的相对介电常数；反演图像中埋于沙层中的两个异常体得到了较好的刻画，两个异常体的水平位置与实际模型相对应.右侧低速异常体的位置与真实模型一致，而且其内部相对介电常数值高于周围背景值；左侧高速异常体的位置与真实模型一致，其内部相对介电常数值低于背景值.

图 4

Fig. 4

图 4 (a)初始介电常数模型； (b)全波形反演25次迭代后的介电常数反演结果Fig. 4 (a) The initial permittivity model for the FWI inversion； (b) the inversion permittivity results after 25 iterations

图 5给出了真实模型、初始模型及反演结果之间的详细数值比较，黑色实线代表真实模型的相对介电常数值，蓝色虚线代表初始模型相对介电常数值，红色实线是时间域全波形反演迭代25次后的反演结果.图 5a是经过左侧圆形异常体A1-A2线(见图 2)的相对介电常数比较，横轴为相对介电常数值，纵轴对应的空间位置.小于0.1 m时，反应在真实模型中是相对介电常数为1空气，可以发现三条曲线完全重合.0.1~0.3 m之间，给出的初始模型与真实模型之间存在1.5的差异，而全波形反演结果基本介于初始模型与真实模型之间.在0.3 m处是界面突变点，反演结果中在同样位置出现了拐点.0.5 m附近是异常体所在的位置，反演结果的数值大小及拐点位置与真实模型基本一致.图 5b是经过右侧圆形异常体B1-B2线(见图 2)的相对介电常数比较，从中可以得到与左图相似的结论：全波反演的结果较为准确的刻画了地下相对介电常数分布情况，可以为逆时偏移提供更为精准的初始模型.

图 5

Fig. 5

图 5 真实模型，初始模型及时间域全波反演结果在不同位置的介电常数比较 (a) A1-A2线(见图2)介电常数比较; (b) B1-B2线(见图2)介电常数比较.Fig. 5 The permittivity comparison of true model, initial model and FWI inversion results at different position (a) Permittivity comparison of A1-A2 line (the left black dash line indicated in Fig.2); (b) Permittivity comparison of B1-B2 line (the right black dash line indicated).

为了评价逆时偏移算法应用到探地雷达数据处理中的效果，选取Kirchhoff偏移算法与之比较.对于逆时偏移成像算法，我们可以将其处理共偏移距数据，亦可以处理共炮点数据.

对于如图 3a所示的共偏移距雷达数据，基于真速度模型，用Kirchhoff偏移成像算法获得的偏移成像结果图 6a所示.从图可知，来自不规则界面的反射波得到一定程度的抑制，可以大致识别出土壤与沙层的界面，然而土壤与沙层界面的散射波没有完全收敛，仍然存在较严重的杂波干扰.偏移剖面中，来自两个异常体的反射波能量基本上得到了压制.左侧高速异常体水平位置与真实模型一致，但是垂直位置与实际模型存在一定的偏差；右侧低速异常体的位置与实际模型基本一致.此外需要注意的是成像结果中出现了一些干扰杂波，这些杂波的存在将使得雷达数据的解释工作更加复杂.因此，对于存在横向速度变化较大的模型，即便使用真实速度模型进行Kirchhoff偏移，也得不到良好的偏移结果.

图 6

Fig. 6

图 6 不同偏移方法应用不同初始模型得到的成像结果 (a) 基于共偏移距数据的Kirchhoff偏移结果(真介电常数模型作为初始模型); (b)基于共偏移距数据的RTM结果(真介电常数模型作为初始模型); (c) 基于共偏移距数据的RTM结果(地下相对介电常数取4作为初始模型); (d) 基于共偏移距数据的RTM结果(地下相对介电常数取6作为初始模型); (e) 基于共偏移距数据的RTM结果(FWI反演结果作为初始模型); (f) 基于共炮点雷达数据集的RTM结果(真介电常数模型作为初始模型). 所有图像都进行了振幅归一化.Fig. 6 Imaging results from different migration methods with different initial model (a) The Kirchhoff migration result based on common offset radar data (true dielectric constant model as the initial model); (b) RTM result based on common offset radar data (true dielectric constant model as the initial model); (c) RTM result based on common offset radar data (dielectric constant of 4 in the underground as the initial model); (d) RTM result based on common offset radar data (dielectric constant of 6 in the underground as the initial model); (e) RTM result based on common offset radar data (FWI result as the initial model); (f) RTM result based on common shot radar data (true dielectric constant model as the initial model). All figures are normalized by their maximum amplitude.

同样，对于如图 3a所示的共偏移距雷达数据，将真介电常数模型作为初始模型参数，忽略电导率的影响，用逆时偏移成像获得的成像结果如图 6b所示，可以发现来自起伏界面的散射波已被完全压制住了，土壤与沙层间的界面可以被清晰的识别出来，成像后的界面形状及位置与真实模型完全一致.来自两个异常体的反射波聚焦较好，成像结果显示两个异常体的形态大小，空间位置与真实模型完全一致，此外偏移剖面中存在的杂波干扰较小.总体来说，逆时偏移成像效果明显好于Kirchhoff偏移结果.

逆时偏移成像算法中，初始模型的选择至关重要，初始模型的建立可以结合已知的先验信息、CMP速度分析、或者反射层析成像结果.为了研究初始模型对成像精度的影响，文中分别选取地下相对介电常数为4的均匀半空间，地下相对介电常数为6的均匀半空间以及选取全波形反演的结果作为初始速度模型.

对于共偏移距探地雷达数据，当选取地下相对介电常数为4的均匀半空间作为初始介电常数模型时，采用逆时偏移技术获得的成像结果如图 6c所示，从图可知来自起伏界面的绕射波没有完全归位，杂波干扰严重，来自地下异常体的散射波能量聚焦效果较差.当选取地下相对介电常数为6的均匀半空间作为初始介电常数模型时，逆时偏移成像结果如图 6d所示，来自起伏界面的反射波能量聚焦较好，能清晰的识别起伏界面的位置及形态，但是来自地下异常体的散射波能量没有完全聚焦.当选取全波形反演的结果作为初始介电常数模型时，逆时偏移成像的结果如图 6e所示，从图可知，来自界面的反射波能量聚焦较好，能够清晰的识别起伏界面的形态和位置；来自地下的两个异常体的散射波能量聚焦较好，其空间位置与模型基本一致.由此可知，当初始模型越接近于真实模型时，逆时偏移成像结果越准确；全波形反演能够为逆时偏移成像提供更为准确的初始模型，将全波形反演的结果作为逆时偏移成像的初始模型是行之有效的.

对于如图 3b所示的共炮点探地雷达数据，采用逆时偏移成像算法对其进行偏移成像，其结果如图 6d所示，图中来自起伏界面的反射波能量得到聚焦，起伏界面的空间位置及形态与实际模型基本一致；来自地下异常体的散射波能量聚焦较好，其空间位置与模型一致.偏移剖面中，在地面处出现了若干个和源位置一致的强能量点；其反射波子波略宽于如图 6b中的反射波子波宽度.

本文在实验室沙槽进行了物理模拟实验，如图 7a所示，沙槽长8 m，宽1.6 m，深0.8 m，沙槽四壁为混泥土，一个PVC盒以45°倾斜角埋于沙中，其埋深约0.24 m.图 7c展示的是PVC盒埋于干沙中的示意图，PVC盒厚度为0.04 m，其中间为空气，顶端位于(0.75，0.24)m处，干沙相对介电常数是4，PVC的相对介电常数是3.1.图 7b展示了本实验中采用的 MALA ProEx雷达控制单元及两套中心频率为1600 MHz的高频天线.当采集共炮点雷达数据时，雷达主机外接两个高频模块，同时将两套高频天线分别接于高频模块；当采集共偏移距雷达数据时，雷达主机外接一个高频模块，同时只需要一个高频天线.

图 7

Fig. 7

图 7 物理模拟实验 (a) PVC盒子倾斜掩埋在沙中; (b) 实验用到的仪器：Mala ProEx控制单元及 两套1600 MHz的高频天线; (c) PVC盒子埋于干沙中的示意图.Fig. 7 Laboratory experiment setup (a) The PVC box tilted buried in sand; (b) Instruments used to carry out radar measurements consisting of Mala ProEx control unit and two sets of 1.6 GHz antenna; (c) A sketch of the PVC box embedded in sand.

首先用1600 MHz高频天线进行了共偏移距测量，测量过程中天线紧贴干沙表面，图 8a所示是实测的共偏移距探地雷达剖面，横轴为道号，纵轴为雷达波双程走时.图中0~1 ns之间是天线耦合直达波，顶点位于3.2 ns的双曲线是来自PVC盒子顶端的散射波.同时可以发现剖面中出现了一些微弱的散射杂波，杂波是由于干沙局部不均匀性造成的.

图 8

Fig. 8

图 8 物理实验实测数据 (a) 实测的共偏移距剖面; (b) 实测的包含36条共炮点 雷达剖面的雷达数据集，所有图像均进行了归一化.Fig. 8 Physical experiments radar data (a) The measured common-offset profile; (b) The measured radar gather consisting of 36 common-shot radar profiles, all figures are normalized by their maximum amplitude.

其次，采集了共炮点雷达数据集，第一个发射天线位于起始位置，第一个接收位置距离发射天线的间距为0.1 m，每一个发射位置对应21个接收位置，即单炮雷达数据包含了21道信号，接收器间的间隔是0.01 m，发射天线的移动间距为0.05 m.图 8b是实测的共炮点探地雷达数据集，它包含36个共炮点探地雷达剖面数据.这些数据将用于时间域全波形反演，以此来反演地下介质分布，从而为叠前逆时偏移成像提供初始速度模型.

实测数据必须经过必要的数据预处理，包括对共炮点探地雷达数据进行直流偏移和带通滤波，同时考虑到本文所用的时间域全波形反演算法是二维的，而实测雷达数据是三维的，故采用(Bleistein，1986)提出的三维转二维技术，将三维实测数据转为二维情况下的数据，从而消除它们在幅度和相位上的差异.由于在沙中仅埋有一个PVC盒子，故而在时间域全波形反演中，采用相对介电常数为4均匀模型作为的初始相对介电常数模型(图 9a).以实测雷达数据和计算雷达数据的均方根作为反演迭代次数的目标函数，如图 10所示，经过15次迭代后，目标函数的误差改变小于1%.图 9b是迭代15次后的相对介电常数反演结果，图中可清晰辨别出倾斜掩埋的PVC盒子，其形状和位置与实际模型基本一致，倾角也非常接近45°，此外反演结果中小散射点清晰可见.

图 9

Fig. 9

图 9 FWI的初始介电常数模型 (a) 相对介电常数定为4; (b) 15次反演迭代后的介电常数结果.Fig. 9 The initial permittivity model for the FWI (a) The dielectric constant is set to 4; (b) The inversion permittivity results after 15 iterations.

图 10

Fig. 10

图 10 以实测雷达数据和计算雷达数据间均方根值作为迭代次数的函数(均方根值经过以做归一化，而且15次迭代后均方根改变值小于1%)Fig. 10 RMS values as a function of the iteration number for the observed radar data and the calculated data (The RMS is normalized to its maximum and after 15 iterations; the RMS misfit changes less than 1%)

通过时间域全波形反演获得了地下介电常数分布情况，基于反演结果给定的初始介电常数模型，分别进行了Kirchhoff偏移成像和逆时偏移成像，其结果见图 11.图 11a是对实测的共偏移距雷达数据进行Kirchhoff偏移成像得到的结果，显然来自PVC盒顶的双曲线同向轴能量收敛到一个点，该聚焦点几乎与PVC盒顶端位置吻合，然而PVC盒身在偏移图像中并没有得到成像，此外，成像结果中出现了许多散射杂点，这些散射杂点不利于雷达数据的解释.图 11b是对实测共炮点探地雷达数据进行逆时偏移成像处理获得的结果，从图可知，成像结果中倾斜PVC盒子整体(白色椭圆框所示)形态及倾角与实际模型基本一致，其偏移结果明显好于Kirchhoff偏移成像结果.

图 11

Fig. 11

图 11 (a)实测共偏移距剖面Kirchhoff偏移结果; (b)实测共炮点雷达数据得到的逆时偏移成像结果，两图均已归一化处理Fig. 11 (a) Kirchhoff migration result from measured common-offset profile, (b) RTM result from measured common-shot radar gather; all figures are normalized by their maximum amplitude

本文开发了电磁波二维逆时偏移成像算法，其中电磁波场的正向传播和反向传播均采用时间域有限差分方法，互相关成像条件用于获得成像结果.考虑到全空间波场存储内存开销较大，文中将成像空间外边界做为存储空间，记录其各个时刻的场值，记录的场值完全能够恢复电磁波场的正向传播过程.

数值模拟实验中，构建了横向速度变化较大的模型，合成了共偏移距和共炮点雷达数据.对于合成的共偏移距雷达数据，同时采用Kirchhoff偏移成像和逆时偏移成像算法，对比成像结果发现，逆时偏移成像剖面中，来自起伏界面及地下异常体的空间位置和几何形态与实际模型一致，而Kirchhoff偏移剖面中来自起伏界面和异常体的反射波(散射波)能量聚焦不完全，剖面中仍然存在明显的杂波，可见逆时偏移成像算法明显好于Kirchhoff偏移算法.然而需要注意的是，逆时偏移成像算法是基于正演计算，每完成一次偏移需要进行两次以上正演计算，而正演计算非常耗时，故而其计算效率远低于Kirchhoff偏移算法.采用逆时偏移成像处理共炮点雷达合成数据，得到的成像剖面与该算法处理共偏移距雷达数据获得的成像剖面基本一致，不同之处在于前者成像剖面中的子波宽度略宽于后者.

影响逆时偏移成像效果的关键因素之一是包括介电常数和电导率在内的初始模型的选择，数值实验表明更精准的初始模型能获得更高精度的偏移成像结果.为了获得更精确的初始模型，本文通过时间域全波形反演算法对探地雷达数据进行全波形反演，反演后的地下相对介电常数模型作为逆时偏移成像的初始速度模型.数值模拟实验及实际物理模拟实验表明，通过全波形反演获取的信息作为逆时偏移成像的初始模型是行之有效的，同时逆时偏移成像的效果在一定程度上检验了全波反演结果的准确性.