基于折射波静校正的未知参数穿墙成像合成延迟时算法*

李家强**1a,，刘 然，徐才秀2，陈金立1a,，陈焱博，朱艳萍

(1.南京信息工程大学 a.气象灾害预报预警与评估协同创新中心;b.电子与信息工程学院，南京210044；2.江苏北斗卫星导航检测中心公司，南京210044)

0 引 言

穿墙成像雷达(Though-the-Wall Radar Imaging，TWRI)作为一种无损检测技术，能够探测视觉非透明障碍物后的目标，在反恐、救灾等方面发挥了重要作用，广泛应用于军事、民用等领域[1]。

近年来，在穿墙成像雷达方面有许多研究致力于对未知参数的单层墙体后目标进行成像。其中文献[2]讨论了相对介电常数和墙壁厚度对目标点成像的影响，通过设置两组与墙壁不同距离的天线阵列，对预先假设的几组墙壁厚度和介电常数的场景进行成像，两组阵列所成图像都会随着墙体厚度和相对介电常数的误差而产生偏移，其中两者交叉的点即为真实目标位置。由于该方法需要对设置的不同墙体参数进行成像，此过程计算复杂，且墙体参数的估计结果误差较大。另外，支持向量机方法(Support Vector Machine,SVM)[3]根据散射信号提取出信号作为输入、墙体参数作为输出得到的函数关系，对墙后目标进行成像。此方法所估计出的墙体参数值易受目标位置、大小以及形状的变化影响。文献[4]运用两列垂直天线阵列对穿墙场景进行成像，利用中心关联目标的位置关系对穿墙场景进行重聚焦，但该方法易受目标的种类和位置的影响。对于多层墙体成像的研究，文献[5]主要针对成像聚焦时延的计算问题，其每层墙体的参数都是假设已知的；文献[6]通过构造滤波器去除回波域中墙体产生的影响以获得清晰的成像结果，该方法每次更新值时都要进行重新计算并成像，故计算效率低，且成像效果易受实际墙体参数的影响;文献[7]通过对墙体厚度和相对介电常数的二维搜索，得到延时补偿值的取值范围，进而对该取值范围进行遍历，找到最佳聚焦时延近似补偿由于墙体存在产生的异常传播时延。在对墙体参数进行搜索时，该方法不能精确地计算聚焦时延的取值范围，导致在遍历时迭代次数过多。

在实际应用中，墙体一般并非单层均匀介质，且墙体参数未知。由于电磁波在墙体内的传播速度比在空气中传播速度慢，导致其回波时延增长，成像过程中目标点位置会产生偏移或出现散焦现象，这给成像工作带来巨大挑战。聚焦时延的计算是影响穿墙成像质量至关重要的因素，因此，解决未知参数多层墙体的时延计算问题是目前穿墙成像面临的挑战之一。为解决未知参数的多层墙体介质，如常见的双层墙体成像问题，本文提出了基于折射波静校正的合成延迟时法对墙后目标进行精准地聚焦成像，通过计算出每一个通道内的聚焦时延，进而依次对每个通道中由于墙体存在产生的异常传播时延进行一一校正。

1 穿墙成像雷达场景模型

本文运用多发多收的天线阵列探测墙后目标。如图1所示，墙体由两层不同介质组成，内层墙体为石膏材料，外层墙体为混凝土材料，两者厚度分别为d1和d2，相对介电常数分别为ε1和ε2。天线阵列紧贴墙体依次摆放，发射信号为超宽带窄脉冲信号[7]，信号经过墙体的透射和反射至目标点处，再由接收天线所接收。

图1 电磁波发射接收传播时延

本文运用基于时域有限差分法(Finite Difference Time Domain,FDTD)的GprMax2D/3D[8]软件对穿墙场景进行仿真并获取回波数据。设发射信号为s(t)为负一阶高斯脉冲信号，其表达式为

(1)

设第m(m=1,2,…,M)根发射天线坐标为(xtm,0)，第n(n=1,2,…,N)根接收天线坐标为(xrn,0)，目标点q坐标为(xq,yq)，在不考虑噪声的影响下，则接收天线接收的回波信号为aqs(t-τmqn)，其中aq为目标的回波幅值，τmqn为电磁波从第m根发射天线到第n根接收天线的传播时延。由图1可知,

(2)

式中：c为电磁波在空气中的传播速度,τtm1与τtm2分别为电磁波在外层和内层墙体中发射天线延迟时，τrn1与τrn2分别为电磁波在两层墙体中接收天线的延迟时,l1与l2为电磁波在空气中的传播路径。由后向投影(Back Projection)算法知，成像场景区域某一像素点p(xp,yp)的像素值为

(3)

式中：τmpn为成像区域任意点p从第m根发射天线到第n根接收天线的传播时延。

若两层墙体参数已知，则可通过折射定理计算折射点坐标，进而计算出电磁波在墙体内的传播时延，完成对墙后目标的成像与精准定位。而在实际应用中，墙体厚度以及相对介电常数这些先验信息通常是未知的，这给后续的成像工作带来影响，导致目标位置偏移和成像散焦。本文提出基于折射波静校正的合成延迟时法，通过计算每一个通道的聚焦时延，使得目标更接近真实位置，图像聚焦程度更高。

2 基于折射波静校正的合成延迟时算法

2.1 折射波静校正原理

静校正技术属于地震勘探领域，是一项处理地震资料的技术。假设存在两层不同介质，电磁波在下层介质中的传播速度大于在上层介质的传播速度时，则透射波会沿水平界面滑动，滑行波在上层介质中激发新的波动即折射波。基于基本折射方程[9]，折射波静校正通过读取折射波初至时间[10]计算介质的厚度以及速度信息，从而得到地震波在介质中的传播时间，完成对由于地表条件起伏引起的时差的校正。故应用折射波静校正必须具备两个基本条件：一是要存在水平的折射面；二是下层速度要大于上层速度[11]。在本文设置的穿墙成像场景模型中，由于墙体设定为平行，故其拥有相对稳定的折射面；另外，由于外层墙体的相对介电常数比内层墙体相对介电常数大，电磁波在介质内的传播速度与其相对介电常数成反比，故电磁波在外层墙体内传播速度比在内层传播速度小。因此，折射波静校正理论运用于穿墙成像领域能够计算电磁波在墙体内的传播时延，以校正每一通道中由于墙体存在引起的异常传播时延误差。

在本文设置的穿墙场景中，墙体由两层不同介质组成。当电磁波从发射天线发出后，一部分透射波会沿着两层墙体的界面滑动，另一部分透射波在内层墙体与空气的界面上滑动，引起墙体振动，产生折射波。折射波的折射路径如图2所示，在本文模型中折射面有两个：内层墙体与空气接触的界面设为折射面1，两层相对介电常数不同的墙体接触的界面设为折射面2。由基本折射方程知：

图2 折射波延迟时路径

(4)

(5)

2.2 外层墙体电磁波传播延迟时

(6)

令

(7)

可得到一系列延迟时超定方程组，则式(5)可写为

Gx=b。

(8)

式中：G=[g1,g2,…,gm]T为该超定方程组的系数矩阵，x为发射天线以及接收天线延迟时构成的未知数向量。对式(8)作差可得

(9)

对式(9)构建最小二乘：

(10)

对x求导得

(11)

式中：tr[xxTGTG]可写为tr[xTGTGx]，整理可得

(12)

其中：

(13)

将式(13)代入式(12)中得

(14)

所以有

(15)

则

(16)

则有

x=(GTG)-1b。

(17)

求得的向量x为该超定方程组的解，即电磁波在外层墙体内的发射点传播延迟时τm,2和接收点传播延迟时τn,2。

2.3 内层墙体电磁波传播时延

(18)

(19)

对于相邻接收天线阵元n与n-1，其截距时间差[17]可表示为

(20)

则求得的相邻发射天线截距时间差即相邻发射天线延迟时间增量，而相邻接收天线截距时间差即为相邻接收天线延迟时间增量。

τ=τm,1+τn,1+τm,2+τn,2+τ″ 。

(21)

其中：

(22)

3 实验与分析

采用GprMax2D/3D建立图1所示的穿墙仿真场景。天线阵列紧贴外侧墙面摆放，采用3发13收天线装置，目标为球形理想导体，位于坐标(1.1,1.2)处。若不考虑墙体存在，使用后向投影算法对目标场景进行直接成像，其成像结果如图3所示，可见由于墙体存在引起的异常传播时延，成像目标位置会严重偏离真实坐标。

图3 原始成像效果图

图4 接收天线1的回波信号波形

(a)回波时刻离散点

如上文所述，首先通过解式(8)超定方程组求得外层墙体的发射天线与接收天线延迟时，再计算每个通道内的发射天线与接收天线延迟时增量来确定内层墙体的延迟时，最后令其两者相加则可得电磁波在两层墙体内的传播延迟时。图6给出了本文算法的流程图。

图6 合成延迟时算法流程图

为验证本文算法在不同墙体两层介质参数下的有效性，设置以下三组不同情况墙体参数值，每组墙体由介质1和介质2组成，相对介电常数和厚度分别为(ε1,d1)和(ε2,d2)，如表1所示。运用基于折射波静校正的合成延迟时算法对三种不同墙体参数情况时穿墙场景的成像结果如图7所示，同时对比了运用回波域滤波方法[6]以及自聚焦方法[7]在不同墙体参数值下的成像情况。与回波域滤波算法对比，本文所提算法对墙后目标位置的定位更准确，且更不易受墙体参数的影响。自聚焦算法成像遍历电磁波的传播时延的取值范围，通过图像评价标准选取成像结果最佳图像确定电磁波传播时延的值，故当墙体参数改变时，其成像目标位置不易受墙体参数影响，而本文算法与自聚焦算法对比，其成像目标位置更精准。

表1 不同墙体参数取值

图7 本文算法与两种算法成像结果对比

为进一步分析算法的性能，对本文算法成像结果的图像熵进行对比。图像熵能够衡量墙后目标物体成像的聚焦程度，其值越小则聚焦程度越好。若像素点个数为X×Y，p(i,j)表示第i行第j列的像素值，则该图像的图像熵ξ可表示为[14]

(23)

根据式(23) 求得在不同墙体参数情况下本文所提算法与回波域滤波算法和自聚焦算法成像结果的图像熵值，其对比如表2所示。

表2 三种算法成像结果图像熵对比

由表2可知，在不同墙体参数环境下，本文算法与两种算法相比，其图像熵值较小，即本文算法成像结果在聚焦程度方面比上述两种算法效果更优。

分别对回波域滤波算法、自聚焦算法以及本文所提算法进行计算复杂度对比分析。对于回波域滤波算法，其计算复杂度可表示为

C0=mnlgp+mN+3N。

(24)

式中：P为采样点数，m、n分别为穿墙成像区域网格横向和纵向的划分个数。对于自聚焦算法，其计算复杂度可表示为

C1=LmnP+mN。

(25)

式中：L为自聚焦算法延时补偿值取值范围的遍历次数。而本文所提基于折射波静校正的合成延迟时算法，其计算复杂度如式(26)所示：

C2=mnlg 2P。

(26)

图8给出了三种算法在不同采样点数情况下计算复杂度的曲线图。由图8可知，在相同采样点数条件下，本文所提算法相较于其他两种算法，其计算复杂度相对较低。

图8 三种算法的计算复杂度对比曲线

为验证在噪声影响下本文所提算法的鲁棒性，在设置同一墙体参数下运用本文算法以及上述两种算法对穿墙场景进行成像，并给出了不同信噪比条件下该算法的输入输出目标杂波比(Target Clutter Ratio,TCR)。由图9可知，相较于自聚焦算法，基于折射波静校正的合成延迟时算法在不同信噪比的噪声情况下能够输出更高的目标杂波比；相较于回波域滤波方法，本文所提算法在噪声信噪比为-3 dB之后能输出更高的目标杂波比。

图9 噪声影响下TCR的变化

4 结 论

针对未知参数的多层墙体穿墙雷达成像目标散焦或位置偏移问题，本文提出基于折射波合成延迟时静校正算法对电磁波的异常传播时延进行校正并成像。该算法首先读取回波信息建立超定方程组计算外层墙体内电磁波的传播时延，再通过合成延迟时增量估计内层墙体每一通道的电磁波传播时延，最后对每一通道电磁波的异常传播时延进行一一校正，使得成像目标被校正到准确位置。仿真结果表明，本文所提算法能够有效校正墙体存在聚焦时延的偏差问题。与传统方法相比，本文所提方法能够更加精确地将目标点校正至正确位置，且图像聚焦效果优于传统算法。理论推导以及程序仿真结果表明，本文方法在未知墙体参数穿墙成像的实际应用中具有一定的参考和指导意义。