1.   引言

探地雷达（Ground Penetrating Radar，GPR）是一种用于探测地下场景的雷达系统，其利用电磁波在介质电磁特性变化处的反射和散射现象来获取地下场景的信息^［1］。GPR系统相比于传统对地探测方式，具有无损性、分辨率高、易于操作等优势，因此在武器探测、冰层冻土探测、管线定位以及考古探测等工程及环境领域得到了广泛应用^［2-5］。

近年来，随着探地雷达应用场景的不断增多，传统GPR系统收发一体、多点采集的数据获取方式越来越不能满足应用中实时数据采集和处理的需求，因此多输入多输出（Multiple Input Multiple Output， MIMO）雷达体制被引入GPR系统中，其利用空间中分布的多个发射正交波形的收发天线，从多角度照射目标^［6］，既明显提高了数据采集速度，也丰富了接收的目标散射数据，从而提升了最终目标的分辨率^［7］。但随之而来的是MIMO-GPR下回波数据处理算法的改进工作，其中偏移成像技术作为GPR数据解释的基础步骤，其处理质量将直接影响后续对地下场景中目标的探测及分类等的准确性。然而近年来对MIMO-GPR的研究多是着眼于MIMO阵列拓扑设计^［8］或利用MIMO雷达体制进行高精度地质探测^［9］及目标定位^［10］等，而针对MIMO-GPR系统的成像算法研究相对较少，因此探索基于MIMO-GPR的成像技术很有必要。

传统偏移成像算法的研究可追溯到20世纪60年代末，Schneider等人首次引入了衍射叠加偏移技术来重建成像区域的反射率图^［11］，即通过反向传播每个采集到的时域距离剖面来将观测数据从采集位置移回其真实的空间位置，这种方法需要对所有接收信号进行后向传播并积分，从而完成偏移区域的衍射叠加计算。由于这种方法实现简单，灵活性强，因此被广泛用于不同场景下的偏移成像处理。然而这种算法也存在一些缺点，例如当回波中只有一个尖峰时，这种算法也会使得其出现在所有可能的反射点的位置，从而带来较大的旁瓣^［12］，这主要是因为该方法并不是基于严格的波动方程理论。为改善衍射叠加方法的缺陷，提高成像质量，Schneider提出了严格基于波动方程的基尔霍夫偏移方法，而后X.Liu等人对其做了进一步完善并应用于诸多成像问题^［13-14］，这种方法代表了惠更斯原理的定量描述，根据Kirchhoff积分定理，利用接收波场值及其导数表示任意点的场值，表现出了比衍射叠加方法更高的成像质量和更低的旁瓣伪影。

由于最初提出该传统Kirchhoff成像方法时主要用于地震数据处理，所以其充分利用了地震数据采集中逐点收发的特点基于爆炸源模型得到了该算法，因而其也仅适用于收发一体的雷达体制，而MIMO探地雷达的收发天线分布于空间中不同位置，这也就直接导致传统成像算法无法直接应用于处理MIMO-GPR采集的收发天线位置不同时的回波数据。为解决这一问题，Zhuge X等人在传统SISO的成像方法的基础上根据MIMO场景对算法进行了直接拓展，提出了适用于MIMO情况下的Kirchhoff成像方法^［15］，但未给出严格证明，也没有考虑远点衰减以及介质分层的情况，因此并不能真正应用于实际探地雷达成像中；而后Zhongmin Wang等人在基于有效相位中心原理和标量衍射理论提出了一种用于近场MIMO成像的距离偏移方法^［16］，实现了频域MIMO成像，但此方法也仅适用于阵元间隔相同的情况。为了弥补之前工作中的不足，本文提出了一种基于MIMO-GPR的时域偏移成像算法，并给出了频域的严格证明，同时引入远点补偿和折射点近似计算方法，实现了对地下目标的精确成像。最后，本文对所提算法进行了仿真验证，仿真结果显示，利用本文方法所得的成像结果相较于利用传统成像方法的结果有更低的旁瓣伪影和更高的图像分辨率。

2.   MIMO-GPR成像算法推导

2.1.   衍射叠加原理

衍射叠加原理是时域偏移成像方法的基础，在连续均匀介质中，未经过聚焦处理的点状反射物的回波数据所表现出的形状称为衍射双曲线，因此在偏移成像时，对于回波中的每一个波峰数据，都要在成像平面中构造一条曲线，然后取每条曲线在该点的强度值并叠加，从而得到该点的像素值。如果所用GPR雷达为收发一体式，且接收到的回波数据为$u(k,t)$，则由爆炸源模型，根据逆时偏移算法原理^［11］，可将待测场景中某一点$(x,z)$的像素值表示为：

其中，$v$代表介质中电磁波的传播速度，$f(x,z)$表示$(x,z)$位置处的像素值，$u(k,t)$表示第$k$道数据中$t$时刻的回波场值，$R_{\mathrm{1}}$和$R_{\mathrm{2}}$分别代表电磁波从发射天线传播到成像点和从成像点传播回接收天线的距离，$w(x,z)$表示加权系数，主要用于补偿收发机的空间分布不均匀所带来的成像误差。

2.2.   Kirchhoff算法的时域推广

从传统的SISO-GPR体制到MIMO体制，所需解决的问题主要是如何充分有效利用多个接收机所采集的全部数据。对于Kirchhoff这类衍射叠加算法而言，可以从算法原理出发进行考虑，衍射叠加算法的本质就是对接收信号进行反向传播并积分，而从这个角度，MIMO体制下的Kirchhoff成像可以认为是对每一个收发对的回波信号进行反向传播并叠加，故对MIMO-Kirchhoff算法成像表达式的推导可以通过对常规Kirchhoff算法进行一般化推广而得到。因此，本文先对常规Kirchhoff算法进行推导，再根据MIMO场景的不同特点对成像公式进行一般化处理。

在均匀各向同性介质中，矢量波动方程可以简化为标量波动方程，利用其Kirchhoff积分计算式并进行时间反向后，可以将波场表示为 ：

此时，可以选定$Q$面为半球面和水平地面形成的闭合曲面。由于在半球面上积分项对场值的贡献为0，可将上式化简，而后令$\frac{\partial r}{\partial z}=-\mathrm{c}\mathrm{o}\mathrm{s}\theta$，$\theta$表示阵列与成像点的连线和竖直方向的夹角，可得常规Kirchhoff成像表达式为：

其中，$S$表示传统SISO-GPR系统数据采集平面上的所有采集点，$v$表示电磁波在介质中的传播速度，$r$表示收发阵元到成像点$(x,z)$的距离，若用$(x_{\mathrm{0}},z_{\mathrm{0}})$表示收发阵元的位置，则距离$r=\sqrt[]{{\left(x-x_{\mathrm{0}}\right)}^{\mathrm{2}}+{\left(z-z_{\mathrm{0}}\right)}^{\mathrm{2}}}$，$\tau$表示接收信号的反向传播时延，即$\tau =\frac{r}{v}$。

但为使得Kirchhoff偏移成像算法可以满足MIMO阵列配置的要求，需要对上式进行修正。首先由于此时雷达系统的天线不再收发一体，所以电磁波发射到目标的路径和从目标反射到接收机的路径不再相同，因此需要用${\theta }_{\mathrm{1}}$和$r_{\mathrm{1}}$来分别表示成像点和发射天线之间的夹角和距离，用${\theta }_{\mathrm{2}}$和$r_{\mathrm{2}}$来分别表示成像点和接收天线之间的夹角和距离，具体几何关系见下图1。

图  1  收发天线和成像点的位置示意图

Fig.  1.  The Schematic diagram of position of receiving and transmitting antenna and imaging point

下载: 全尺寸图片 幻灯片

此时上述公式可以修正为

其中$S_{\mathrm{0}}$和$S_{\mathrm{1}}$分别表示雷达系统中的发射孔径和接收孔径。

2.3.   MIMO-Kirchhoff算法的频域证明

上述表达式是在常规Kirchhoff算法的基础上针对MIMO雷达的不同之处而做出的普适性外推，本节中，我们将在频率-波数域中对上述成像表达式展开推导证明。

首先假设所有发射机的发射信号强度都相同，并通过傅里叶变换将发射信号和目标散射到接收天线的场变换到波数-频率域，而后根据波场延拓关系^［9］，可以将发射波场从发射天线延拓至散射点处，则可得到散射前的波场为：

其中$∆z^{\text{'}}=z_{\mathrm{1}}-z_{\mathrm{0}}^{\text{'}}$表示波场延拓距离，$U_{S_{\mathrm{0}}}^{\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}}$表示发射天线辐射出的源波场，$H(k_x,∆z,w)$函数表示电磁波场在介质中的外推延拓运算，$w$表示电磁波信号的角频率。同样的可以将接收波场逆向延拓至散射点处，即得散射后的波场为：

其中，波场延拓距离$∆z=z_{\mathrm{1}}-z_{\mathrm{0}}$，$U_{S_{\mathrm{0}}}^r$表示接收天线收到的波场。由于当前散射场中介质分界面相对明显，可以进一步假设Born近似在本散射过程中仍然适用，则可以得到散射后的波场和散射前波场的关系式为：

式中，$f_r$表示散射体的散射系数。

一般情况下，成像处理对分辨率要求较高且待成像区域普遍较浅，因此GPR所用发射信号通常为高频脉冲信号，有时还会采用超宽带（UWB）信号^［17］，所以此处可以将发射信号理想化为一冲激信号，则发射波场$U_{S_{\mathrm{0}}}^{\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}}$为一常数，同时根据上述$H(k_x,∆z,w)$函数的特点可以得到：

而后将上式逆傅里叶变换回空间-时间域，可得：

根据参考文献［13］可得式中，

而后将式（10）中的卷积运算写成二维积分形式，并对其进一步离散、化简，即可得到：

其中，$r_{\mathrm{1}}$和$r_{\mathrm{2}}$分别表示发射天线和接收天线到成像点的距离，$\tau$和${\tau }^{\text{'}}$分别表示发射天线和接收天线到成像点的时延。此时，反射系数在空间中的分布反映的就是最终GPR逆时偏移处理所得的成像结果。另外，按照图1中所示的几何关系，可得上式和式（4）所表示的MIMO-Kirchhoff算法成像表达式是一致的。

2.4.   远点补偿处理

根据电磁波在介质中的传播规律可知，当介质为非线性时，电磁波在介质中传播时除了相位会以传播常数随距离变化外，其幅值也会以衰减常数随距离指数衰减。因此，距离收发天线越远的目标，其回波的波场幅值越低。而上述成像表达式（12）中，由于被积分项的系数亦会随着成像点离收发天线的距离增大而减小，因此利用上式对雷达回波数据进行成像处理，必定会使得本就微弱的深处目标的场值变得更小，从而带来较为严重的有限孔径问题，故需要对上式进行修正以补偿远点的场值。此处，为了给远点目标以更大的权值，可以将成像表达式中的系数修正为$r_{\mathrm{1}}{r}_{\mathrm{2}}$，即可得到最终的成像表达式：

3.   空耦下的MIMO-GPR成像处理

3.1.   折射点计算

对于实际应用中的一般情况，GPR天线通常都无法布设于待探测场景所在的介质中，而往往是位于空气中，且距待探测场景的界面有一定距离，此时的成像场景示意图如下图2。

图  2  空耦场景下收发阵列和成像点的几何关系示意图

Fig.  2.  The Geometric relationship between transmitter receiver array and imaging point in space coupled scene

下载: 全尺寸图片 幻灯片

从上图可以看出，由于空气和介质分界面的作用，电磁波的传输路径不再是一条直线，而是一条被折射点$\left(x_r,\mathrm{0}\right)$分为两段的折线，此处为了几何关系表示清晰，只绘制了发射天线发出的电磁波到达成像点的路径，而根据光路可逆原理，电磁波从成像点到接收天线的路径与此相似，故此处只需研究发射路径即可。根据Snell折射定律以及图中所示的几何关系，可得

其中，${\theta }_i$和${\theta }_r$分别表示入射角和折射角。结合上式（14）~（15），可得

可以发现，这是一个关于$x_r$的一元四次方程，而由一元多次方程的一般理论可知此方程没有一般通用的求解方法，在实际中一般采用牛顿法等数值解法进行求解。但由于牛顿迭代法需要较多的迭代次数，计算量很大，会使得算法的运算速度大大下降，所以Mast等人提出了一种折射点近似求解方法^［18］，可将复杂的求解过程简单化，其给出的计算公式为：

其中$x_l=\frac{z_A}{z_A+h}{x}_k+\frac{h}{z_A+h}{x}_A$，而后Zhou对其进行了进一步优化^［19］，得到了当成像点和收发天线方位向距离较大时仍适用的近似求解折射点的方法：

3.2.   空耦下的MIMO-Kirchhoff算法

在空耦合的情况下，利用波场延拓关系，重新表示散射前和散射后的波场。由于此时的介质不连续，且波场在介质分界面处折射前后只有幅值的变化，而不发生相位的移动，因此可以将上述二式重新表示为：

其中$H_{\mathrm{1}}(k_x,∆z,w)$表示电磁波场在空气层中的外推延拓运算。

而后将上述两式代入式（8）中，并对其进行化简、逆傅里叶变换和离散后，可得空耦合情况下MIMO-Kirchhoff的成像表达式为：

其中$r_a$和$r_b$分别表示发射天线到折射点和折射点到散射点的距离，$r_a^{\text{'}}$和$r_b^{\text{'}}$分别表示散射点到折射点和折射点到接收天线的距离，${\epsilon }_{r\mathrm{1}}$和${\epsilon }_{r\mathrm{2}}$分别表示上下两层介质的介电常数。

至此，本文成功利用波场延拓关系，将MIMO-Kirchhoff算法拓展到空气耦合的应用场景中，并得到了空耦适用的成像表达式。

4.   数值仿真

为了充分验证算法的性能，本文分别建立了地耦合多点目标模型和空耦合W形目标模型，并利用基于FDTD的仿真工具gprMax^［20］对上述模型进行仿真，从而获得了雷达回波数据。在下文中，将对上述两种模型及相应的仿真结果分别进行介绍和分析。

4.1.   地耦合多点目标模型

为了说明文中MIMO-Kirchhoff算法的有效性以及远点补偿处理的必要性，建立了如下图3所示的地耦合多点目标模型。

图  3  地耦合多点目标模型场景示意图

Fig.  3.  Scene sketch of ground coupled multi point target model

下载: 全尺寸图片 幻灯片

如下图所示，本仿真场景中$Z=\mathrm{0}$处为介质分界面，收发天线阵均为线阵，且紧贴分界面，其中天线均为赫兹偶极子天线，发射和接收阵元均从$x=\mathrm{0.2}\mathrm{m}$处开始布置，发射阵元间隔为0.16 m，接收阵元间隔为0.032 m；介质的介电常数${\epsilon }_r$=6，电导率为0，其中有3个半径均为0.03 m的小球，具体位置如图3所示。

由于小球半径相对较小，同时也为获得较高的距离向分辨率，本模型所用的发射信号频率为6 GHz，波形为Ricker子波。利用gprMax对模型进行仿真，可以得到该模型的B-Scan回波如下图4（a）所示。

图  4  多点目标模型B-Scan图

Fig.  4.  B-scan graph of multi-point target model

下载: 全尺寸图片 幻灯片

从图4（a）中可以发现，回波数据中最强波为界面直耦波，利用空场景对消法将其去除后，得到的回波如图4（b）所示。从图中可以较明显的看到三条双曲线，这也就代表地下的三个点目标。而后用与MIMO阵列中接收阵元数相同道数的传统SISO体制的Kirchhoff算法^［21］作为对比，分别绘制传统算法的成像结果和本文中未经远点补偿的MIMO-Kirchhoff算法的成像结果如下图5所示。

图  5  基于不同算法的多点目标成像结果

Fig.  5.  Multi point target imaging results based on different algorithms

下载: 全尺寸图片 幻灯片

从图5（a）中可以看到，传统算法得到的像中包含很多旁瓣伪影，且由于采集到的目标散射信息不够丰富，距离维分辨率也略差一些，而由图5（b）可以看出未经补偿的MIMO-Kirchhoff算法虽然可以实现对较近目标的精确成像，但深层目标在成像结果中基本完全不可见。

而后绘制补偿后的MIMO-Kirchhoff算法成像结果，如下图6所示，可以看出成像结果中不仅没有明显的旁瓣伪影，还实现了深层目标的精确成像。但对于浅层目标，补偿后的成像质量不如补偿前好，这主要是因为远点补偿处理在增强深层目标的同时也增强了近目标附近的噪声，略微降低了成像质量。

图  6  补偿后的MIMO-Kirchhoff算法成像结果

Fig.  6.  Imaging results of MIMO-Kirchhoff algorithm after compensation

下载: 全尺寸图片 幻灯片

为定量描述算法的成像质量，本文引入积分旁瓣比（Integrates Side Lobe Ratio， ISLR）来衡量几种成像算法对图像中旁瓣能量的抑制程度。ISLR的定义为目标所有旁瓣能量与主瓣能量的比值：

其中$E_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$和$E_{\mathrm{m}\mathrm{a}\mathrm{i}\mathrm{n}}$分别表示成像结果的总能量和目标主瓣能量。从上式可知，图像的ISLR值越低，说明算法对图像中旁瓣和杂波的抑制能力越强，成像质量越好。

而后本文分别测量上述三种算法所成图像中目标点的位置坐标，并计算其各自在方位维和距离维的ISLR值，最终将所得结果整理如下表1所示。

表  1  三种算法的成像指标对比

Tab.  1.  Comparison of imaging indexes of three algorithms

成像算法	成像目标位置	方位维ISLR	距离维ISLR
传统Kirchhoff成像算法	(0.600,0.093)(0.600,0.245)(1.206,0.246)	4.540	1.775
补偿前MIMO-Kirchhoff算法	(0.600,0.095)	0.617	-2.071
补偿后MIMO-Kirchhoff算法	(0.600,0.095)(0.600,0.248)(1.203,0.247)	0.940	-0.546

下载: 导出CSV 

| 显示表格

由上表1可知，根据成像结果中的目标位置，传统Kirchhoff算法所成像的位置偏差最大，平均偏差距离为0.0064 m，而补偿后的MIMO-Kirchhoff算法成像目标位置偏差最小，平均偏差距离为0.0037 m；同时，根据成像结果的ISLR对比，传统Kirchhoff成像算法对图像旁瓣和杂波的抑制能力最差，而补偿前MIMO-Kirchhoff算法对杂波的抑制能力最强，因此当目标处于较近距离时使用此种算法最佳，当目标处于较远距离时则需要对算法进行补偿，以显示出深层目标。

4.2.   空耦合W形目标模型

为了说明空耦合情况下MIMO-Kirchhoff算法的有效性，本文建立了如下图7（a）所示的空耦合W形目标模型。其中，天线阵列配置和介质参数设置均与图3中相同，但天线阵列距介质分界面的竖直距离为0.1 m，介质中目标的形状为W形，其各顶点的位置如图7（a）所示，发射信号仍采用6 GHz的Ricker波形。

图  7  空耦合W形目标模型的真实场景及成像结果

Fig.  7.  Real scene and imaging results of space coupled W-shaped target model

下载: 全尺寸图片 幻灯片

利用gprMax对上述目标模型进行建模仿真，并对回波数据进行空场景对消处理，而后利用式（21）所示的成像表达式对处理后的数据进行偏移成像，最终的成像结果如图7（b），从中可以看出本文所提的空耦MIMO-Kirchhoff算法在地下场景相对复杂时也可以达到不错的成像效果，充分说明了本算法的有效性。

5.   结论

本文通过对传统探地雷达偏移方法进行推广，提出了一种基于Kirchhoff的MIMO探地雷达成像方法，解决了传统探地雷达不能实时采集数据并成像的缺陷。通过远点补偿修正成像表达式中的加权系数，实现了对较深层目标的探测和成像，而后引入折射点近似计算方法并对成像表达式进行修正，得到了空耦场景下的MIMO-Kirchhoff算法，最后本文利用gprMax对地耦合多点目标模型和空耦合W形目标模型进行建模和仿真，结果表明该方法在地耦和空耦场景下都能达到较好的成像效果，具有较高的分辨率，可以实现对地下场景或目标的精确成像。

由于算法权值的作用，在增强远点的目标的同时也增强了部分噪声，使得成像结果的ISLR降低，在后续工作中可以考虑利用Sigmoid等非线性加权函数对成像表达式的权值进行修正。