分布式穿墙成像雷达由多个雷达节点与融合中心组成，多个雷达节点部署于多个视角，对建筑环境中多个人体目标进行协同探测，典型工作场景如图1所示。

假设雷达节点配备$M$个发射单元与$ N $个接收单元，$M$个发射阵元依次发射超宽带脉冲信号$s(t)$，发射信号经过目标与环境反射后被$ N $个接收阵元接收。当节点$i$的阵元${T_m}$发射信号时，接收阵元${R_n}$接收到的回波信号可表示为：

式中，$\sigma _{mn}^i$为目标复散射系数；$\tau _{mn}^i$为信号传播时延；$w_{mn}^i(t)$为墙体、家具等静止强散射体回波；$\beta _{mn}^i(t)$为多径杂波；$\varphi _{mn}^i(t)$为随机噪声。

图  1  分布式穿墙成像雷达目标探测示意图

由于人体目标散射截面积远远小于墙体等建筑结构，同时电磁波在墙体内部的衰减导致信号穿透墙体后能量大幅损失，式(1)中静止强散射体回波$ w_{mn}^i(t) $强度远大于目标回波。考虑到不同探测周期中静止强散射体回波基本保持不变，而人体目标存在轻微晃动或大幅运动，本文采用脉冲对消器^[11] 抑制墙体回波，经过对消后回波表示为${\bar r^i}_{mn}\left( t \right)$。

反投影成像算法^[12] 因其实现简单、无需探测场景先验信息等优点广泛应用于穿墙透视成像中。反投影成像算法的基本思想为将$MN$个通道回波分别映射到二维成像空间，并相干叠加得到二维雷达图像。将探测区域离散为$Q$个像素点，根据反投影成像算法，节点$i$成像结果中像素点${{\boldsymbol{x}}_q}$的像素值计算如下：

式中，${\tilde \tau _{mn}}({{\boldsymbol{x}}_q})$为计算时延。本文采用固定时延补偿法^[13] 补偿电磁波穿透墙体引起的时延，并基于能量相干因子算法^[14] 抑制旁瓣。

分布式穿墙成像雷达系统中，所有雷达节点基于自身坐标系进行目标定位跟踪。假设每个雷达节点阵列中心位于自身坐标系的原点处，且其他雷达节点位置与角度未知。为了融合多节点信息，需要建立统一的坐标系。

本文将其中一个雷达节点作为参考节点，将其自身坐标系作为参考坐标系，其他节点坐标系可以表示为参考坐标系的旋转平移坐标变换，如图2所示，进而通过求解各个节点的坐标变换参数实现节点配准。

图  2  雷达节点配准示意图

假设目标$p$在参考坐标系中位置为${{\boldsymbol{y}}_p}$，则雷达节点$i$检测结果中目标$p$的量测可表示为：

式中，${{\boldsymbol{\beta}} _i} = {\left[ {\Delta {x_i},\Delta {y_i}} \right]^{\rm{T}}}$为偏移量；H为观测矩阵；${\boldsymbol{v}}_p^i$为量测误差；${{\boldsymbol{\varPhi }}_i}$为旋转矩阵，定义为：

式中，${\varphi _i}$为旋转角度。

由于雷达节点分辨率的限制，目标图像方位向扩展与目标距离成正比，因此不同雷达节点成像结果中目标幅度与形状均不相同^[15] ，因此需要根据成像结果估计目标位置与量测误差。二值图像${I_{{\rm{bin}}}}(x,y)$的$\left( {i + j} \right)$阶矩定义为^[16] ：

本文采用目标图像形心作为雷达量测，假设量测误差服从高斯分布，并正比于图像扩展。则节点$i$获取的量测值$ {\boldsymbol{z}}_p^i $与量测误差协方差矩阵${\boldsymbol{R}}_p^i$可根据图像矩计算为：

式中，$\rho $为尺度因子。

假设雷达节点$i$获取的量测点集合为${{\boldsymbol{M}}_i} = \left( {\boldsymbol{z}}_1^i, {\boldsymbol{z}}_2^i, \cdots , {\boldsymbol{z}}_{{m_i}}^i \right)^{\rm{T}}$，其中${m_i}$为节点$i$获取的量测点数量。参考节点获取的量测点集合表示为${{\boldsymbol{M}}_0}$。

雷达节点配准本质为将各个雷达节点的坐标系经过旋转与平移变换后对齐到参考坐标系，该问题等效为寻找一组旋转平移变换参数${{\boldsymbol{\theta }}_i}$，使坐标变换后的点集$T\left( {{{\boldsymbol{M}}_i},{{\boldsymbol{\theta }}_i}} \right)$与参考点集${{\boldsymbol{M}}_0}$之间距离最小。因此雷达节点配准可表示为：

式中，$f{\mkern 1mu} ( \cdot )$为两个点集之间距离度量函数；${{T}}\left( {{{\boldsymbol{M}}_i},{{\boldsymbol{\theta }}_i}} \right)$表示点集${{\boldsymbol{M}}_i}$的旋转平移变换，变换参数为${{\boldsymbol{\theta}} _i}$。

对于二维配准问题，需要估计${{\boldsymbol{\theta}} _i} = \left[ {\Delta {x_i},\Delta {y_i},{\varphi _i}} \right]$中3个参数，其中$\Delta {x_i}$与$\Delta {y_i}$分别为$x$与$y$方向平移参数，${\varphi _i}$为旋转参数。则点集${{\boldsymbol{M}}_i}$的旋转平移变换可表示为：

式中，$ \otimes $为Kronecker积；${{\boldsymbol{1}}_{{m_i} \times 1}}$表示${m_i}$行1列的全1矩阵。

在节点$i$自身坐标系下，目标$p$位置${\boldsymbol{x}}_p^i$的概率密度函数可表示为：

式中，$\det \left( \cdot \right)$为矩阵行列式。节点$i$观测到的所有目标位置的概率密度函数可表示为${m_i}$个高斯模型的叠加，得到高斯混合模型：

式中，${\boldsymbol{z}}_j^i$为第$j$个量测值；${{\boldsymbol{R}}_j}$为对应量测误差协方差矩阵。

采用${L_2}$距离作为上式中两个点集之间的距离度量函数，两个高斯混合模型的${L_2}$距离定义为：

式(12)中${L_2}$距离共包含3项，第一项与坐标系变换无关；对于刚性变换(对称、平移与旋转变换)，第二项为固定值。因此${L_2}$距离可化简为两个高斯混合模型乘积的积分：

两个高斯分量乘积的积分表达式为：

则式(8)中$f{\mkern 1mu} ( \cdot )$可表示为：

式中，$\tilde {\boldsymbol{z}}_j^i$表示经过坐标变换后的量测点；$\bar {\boldsymbol{R}} = {\boldsymbol{R}}_l^0 + {\boldsymbol{\tilde R}}_j^i$，${\boldsymbol{\tilde R}}_j^i$为变换后量测误差协方差矩阵。

${\boldsymbol{\tilde R}}_j^i$可通过奇异值分解算法计算，首先对原始误差协方差矩阵${\boldsymbol{\tilde R}}_j^i$做奇异值分解，得到：

式中，${\lambda _1}$、${\lambda _2}$为奇异值；向量${({u_{11}},{u_{21}})^{\rm{T}}}$与${({u_{12}},{u_{22}})^{\rm{T}}}$为两个主方向的方向向量。

由于旋转与平移变换只改变两个主方向，则变换后的量测误差协方差矩阵${\boldsymbol{\tilde R}}_j^i$计算为：

在雷达节点配准问题中，变换参数${{\boldsymbol{\theta }}_i}$的估计精度主要取决于点集${{\boldsymbol{M}}_0}$与${{\boldsymbol{M}}_i}$中量测点的数量、量测误差大小、量测点的空间分布与优化算法。当量测点过少时，配准问题可能产生多解。为了在完全无先验信息的条件下实现雷达节点配准，需要点集${{\boldsymbol{M}}_0}$与${{\boldsymbol{M}}_i}$中至少存在3个源于共同观测目标的量测点。

当探测场景中目标数量不足3个时，可采用多周期量测点进行配准，此时式(8)中$f{\mkern 1mu} ( \cdot )$可表示为：

式中，${{\boldsymbol{M}}_i} = \bigcup {\boldsymbol{M}} {\mkern 1mu} _i^k$；${\boldsymbol{M}}{\mkern 1mu} _i^k$为第$k$周期节点$i$获取的量测点集合。

本文采用粒子群优化(particle swarm optimization, PSO)算法^[17] 求解式(8)中的优化问题，流程如下。

1) 随机初始化$J$个三维粒子，3个维度分别对应${{\boldsymbol{\theta }}_i}$的3个参数。每个粒子具有位置与速度两个属性，粒子位置表示为${{\boldsymbol{P}}_j} = \left[ {{p_{j1}},{p_{j2}},{p_{j3}}} \right]$，粒子速度表示为${{\boldsymbol{V}}_j} = \left[ {{v_{j1}},{v_{j2}},{v_{j3}}} \right]$，其中${p_{jd}} \in \left[ {{P_{d\min }},{P_{d\max }}} \right]$，且${p_{jd}} + {v_{jd}} \in \left[ {{P_{d\min }},{P_{d\max }}} \right]$，${P_{d\min }}$与${P_{d\max }}$分别为第$d$个维度参数上下界；

2) 根据式(15)或式(18)计算每个粒子处目标函数值${f_j}$；

3) 对于每个粒子，将其当前位置处目标函数值${f_j}$与其移动过程中计算得到的最小的目标函数值$f_j^l$作比较，若${f_j} < f_j^l$，则将当前粒子位置作为局部最优位置${\boldsymbol{P}}{\mkern 1mu} _j^l = \left[ {p_{j1}^l,p_{j2}^l,p_{j3}^l} \right]$；

4) 对于每个粒子，将其当前位置处目标函数值${f_j}$与所有粒子移动过程中计算得到的最小的目标函数值${f_n}$作比较，若${f_j} < {f_n}$，则将当前粒子位置作为全局最优位置${{\boldsymbol{P}}_n} = \left[ {{\mkern 1mu} {p_{n1}},{p_{n2}},{p_{n3}}} \right]$；

5) 更新所有粒子的位置与速度：

式中，$\omega $、${c_1}$、${c_2}$为调节因子；${r_{1d}}$、${r_{2d}}$为0～1之间均匀分布的随机数。

6) 若目标函数值${f_n}$小于预设阈值或达到最大迭代次数，将全局最优位置作为优化结果，即${{\boldsymbol{\theta }}_i} = {{\boldsymbol{P}}_n}$；反之，返回步骤2)。

为了验证基于高斯混合模型的雷达节点配准算法的有效性，本文采用数值计算软件模拟分布式穿墙成像雷达多目标探测。仿真中采用两个雷达节点探测两个运动目标，全局坐标系下，仿真场景如图3所示。

图  3  仿真场景

两个节点均采用两发四收天线阵列，全局坐标系下节点1阵列中心设为(5 m, 0 m)，发射天线坐标分别为(4.85 m，0 m)与(5.15 m，0 m)，接收天线位置分别为(4.8875 m，0 m)、(4.9625 m，0 m)、(5.0375 m，0 m)与(5.1125 m，0 m)，阵列角度0°；节点2阵列中心设为(0 m, 5 m)，阵列布局与节点1相同，阵列角度90°。两个雷达节点均发射步进频连续波信号，频率范围为1.6～2.2 GHz。为简化问题分析，仿真中没有加入墙体。两个目标均作匀速直线运动，目标1起点为(2 m, 3 m)，终点为(8 m, 8 m)；目标2起点为(8 m, 2 m)，终点为(3 m, 7 m)。两个目标同时开始运动，运动过程共持续8个周期。

产生多周期回波信号后，采用逆傅立叶变换算法进行脉冲压缩。脉冲压缩结果中加入了高斯白噪声，以最大值作为信号强度，信噪比设置为8 dB。部分周期成像结果如图4所示。

图  4  雷达节点成像结果

每周期成像结果经过阈值检测后得到目标区域，进而根据式(6)与式(7)分别计算目标形心与量测误差协方差矩阵，得到量测点集。

以节点1自身坐标系为参考坐标系，采用本文提出的基于高斯混合模型的雷达节点配准算法计算节点2旋转与平移变换参数${\boldsymbol{\theta }} = \left[ {\Delta x,\Delta y,\varphi } \right]$，实现节点配准。粒子群优化算法中，$\omega = 0.6$，${c_1} = 0.9$，${c_2} = 0.9$，$\Delta x$与$\Delta y$搜索空间均设置为−20～20 m，$\varphi $搜索空间为−180°～180°，粒子数$J = 500$。然而，由于搜索空间过大时粒子群优化算法不易收敛，且配准时阵列角度参数对配准结果影响较大，因此实现中将角度搜索空间等分为36个区间，每次对10°范围内进行优化求解，最后将36个区间内最优解作为配准结果。

雷达节点配准结果如图5所示。图5a为配准前两个节点在自身坐标系下获取到的量测，实心点为节点1量测，空心点为节点2量测；图5b为配准结果。由于目标图像存在扩展与轻微形变，目标图像形心不能完全反映目标准确位置，因此节点2量测点集经过坐标变换后无法与节点1量测点集完全重合。

图  5  雷达节点配准仿真结果

为了进一步验证雷达配准算法的准确性与鲁棒性，分别采用500次蒙特卡洛仿真测试无漏检无虚警、有漏检无虚警、有漏检有虚警条件下的配准算法性能。在有漏检条件下，每个目标检测概率设置为80%；在有虚警条件下，每个周期有10%的概率在探测区域中随机位置出现假目标。此外，为了测试粒子群算法中粒子数量对参数估计精度的影响，分别计算了粒子数$J = 100$与$J = 500$时，3种条件下节点2的变换参数，并计算变换参数与真实值的误差，6种情况下误差累积分布曲线如图6所示。

图  6  不同条件下误差累积分布曲线

由结果可知，粒子数为500时，3种条件下位置误差小于0.05 m且角度误差小于1°的概率大于90%，粒子数减少至100时，配准精度下降。当粒子数相同时，无漏检无虚警条件下配准误差明显小于有漏检无虚警及有漏检有虚警条件。对于雷达配准问题，漏检与虚警均会导致另一个节点获取到的量测点集中没有量测点与之对应，从而增大配准难度，导致配准误差增大。

在无漏检无虚警条件下，采用500次蒙特卡洛仿真测试了文献[9] 中提出的基于图模型的雷达节点配准算法性能，并与本文算法进行对比，得到两种算法误差累积分布曲线如图7所示。

图  7  两种算法误差累积分布曲线

基于图模型的雷达节点配准算法利用多节点成像结果中多目标空间分布结构的相似性，构建完全图模型求解两个节点探测结果中对应目标的关联关系，最后计算出坐标变换参数。然而，该算法需要已知两个节点共同观测到的目标数量，因此无法处理漏检与虚警情况。此外，该算法未考虑量测误差的影响，量测误差较大时配准精度降低。粒子数为500时，本算法配准性能优于基于图模型的雷达节点配准算法；当粒子数降至100时，本算法性能略低于图模型算法。由于本算法无需已知两个节点共同观测到的目标个数，且配准过程中考虑到了不同量测点量测误差的差异，在实际应用中具有更强的鲁棒性。

为了验证实际应用中本文提出的雷达节点配准算法的有效性，采用两个两发四收穿墙成像雷达节点协同探测房间内两个运动人体目标，雷达节点参数与仿真相同。实测场景如图8a所示，房间尺寸约为8 m×8 m，墙体材质为实心砖，厚度约为0.22 m。以墙角为原点，建立全局坐标系如图8b所示。在全局坐标系中，节点1阵列中心坐标为(4 m, 0 m)，节点2阵列中心坐标为(0 m, 5 m)。两个目标均做匀速直线运动，在全局坐标系下，目标1起点为(2.5 m, 4.5 m)，终点为(5.0 m, 1.5 m)；目标2起点为(6.0 m, 5.5 m)，终点为(4.5 m, 3.0 m)。

图  8  实验场景

目标探测共持续18个周期，雷达节点信号处理方法与仿真相同。两个节点成像运算均在自身坐标系下进行，成像范围为10 m×10 m，像素点尺寸设置为0.03 m，部分周期成像结果如图9所示。

图  9  运动目标成像结果

由于墙体材质的非均匀性，目标图像出现扩展与形变，且不同周期目标图像形状差异明显，甚至部分周期出现图像分裂。此时目标图像形心明显偏离实际位置，雷达获取到的量测点存在较大量测误差。同时由于运动过程中目标姿态实时变化，不同周期目标相对强度存在显著差异，部分周期出现图像缺失。此外，两个节点成像结果中墙体位置附近均出现遮蔽鬼影，带来了大量虚假量测。

经过目标检测后，两个节点获取到的量测点如图10a所示。进而采用本文提出的基于高斯混合模型的雷达节点配准算法计算节点2的旋转平移变换参数，粒子群优化器参数与仿真相同。雷达节点配准结果如图10b所示。

图  10  雷达节点配准结果

由于本实验中目标图像形变与扩展与仿真相比更为严重，雷达节点量测误差更大，经过配准后两个节点的量测点集仍然存在一定偏差。与雷达节点真实位置相比，20次实验中本算法输出配准结果平均位置误差为0.19 m，平均角度误差为3.15°。此外，由于该实验中漏检与虚警问题严重，无法获取两个雷达节点共同观测到的目标数量，文献[9] 中基于图模型的雷达节点配准算法无法使用。

本文提出了一种基于高斯混合模型的分布式穿墙成像雷达空间配准方法，将雷达节点获取的量测点集建模为高斯混合模型，并采用${L_2}$距离计算各个节点量测点集与参考节点量测点集之间的距离，最终基于粒子群优化算法求解出旋转平移变换参数。该算法无需额外辅助定位设备，仅根据各个节点探测结果实现雷达节点空间配准。仿真与实测结果表明，与现有算法相比，该方法无需求解多节点多目标之间的关联问题，且对于量测误差具有较好的鲁棒性。