1.   引言

天波雷达(Over-The-Horizon Radar, OTHR)利用电磁波在电离层与地/海面之间的反射作用传输高频能量，可实现对高价值目标的远程预警^[1-3]。由于目标与雷达之间存在多条电离层传播路径，单目标在雷达极坐标系下会产生多路径量测，进而可能导致多路径航迹，使得态势感知和情报出现误判。因此，多路径量测为天波雷达数据处理带来严峻挑战。现有的天波雷达数据处理主要集中在目标跟踪，其按照多路径量测融合阶段的不同可以分为航迹级跟踪融合和量测级跟踪融合。航迹级跟踪融合将问题分解为单路径目标跟踪和多路径航迹融合两个子问题^[4,5]。量测级跟踪融合将单路径数据关联扩展到多路径数据关联，通过目标-量测-路径关联，在量测级融合目标信息^[6,7]，性能较航迹级跟踪具有明显提升。量测级跟踪在预处理阶段需要进行航迹起始^[6,7]，即融合同一个目标的多路径量测得到航迹的初始状态。错误的航迹起始可能导致航迹冗余或目标漏跟。然而，目前的量测级跟踪算法主要关注跟踪阶段，很少研究航迹起始。

天波雷达航迹起始本质上需要解决多路径量测聚类问题。由于同一个目标经不同传播路径在雷达极坐标系得到的多路径量测不具备空间相邻特性，无法直接进行聚类。而多路径量测通过正确的传播路径反变换到目标坐标系后可以表示目标的状态估计，因此具有空间相邻特性，可以利用该特征进行聚类。故天波雷达量测聚类需要对同一个目标的多路径量测进行传播路径辨识和聚类融合，两者相互耦合，需要统一处理。同时，天波雷达量测模型假设1个目标通过1种传播路径至多产生1个量测^[4,5]，因此1个量测只能选择1种传播路径，且来自同一种传播路径的量测不能聚类，这两个约束统称为多路径聚类约束。

在数据挖掘领域，针对一大类数据模型已经提出了多种聚类算法，例如基于密度的聚类算法(Density-Based Spatial Clustering of Applications with Noise, DBSCAN)^[8]、K均值^[9]、相似性传播(Affinity Propagation, AP)^[10-13]等，但是针对多路径约束聚类问题的研究却很少。在传感器数据处理领域，聚类算法也存在广泛应用。针对多检测传感器系统，文献[14,15]提出一种两阶段多目标跟踪算法，在第1阶段将量测值自动分组为确定数量的聚类团，在第2阶段对聚类团进行目标跟踪。对于大型传感器网络，跟踪前进行量测聚类融合也是一种有效的解决方案^[16-18]。文献[19,20]在先验信息较差甚至没有先验信息的多传感器融合场景，采用量测聚类进行目标检测与状态估计。上述传感器量测聚类算法都没有考虑多路径约束，无法直接应用于天波雷达多路径量测聚类。对于天波雷达多路径量测聚类问题，专门研究较少，可直接采用多路径多假设(Multiple Hypothesis, MH)航迹融合^[4,5]的静态假设树构造算法，通过枚举量测可能来源的目标和传播路径递归构建假设树，同时评估每个聚类假设的概率及目标状态。该算法为最优解算法，其聚类假设数目是量测和传播路径数目的指数函数^[5]，可以通过剪枝提升计算效率，但是本质上仍具有指数时间复杂度，无法满足杂波环境下目标跟踪实时性要求。

AP聚类是文献[10-13]提出的一种聚类算法，其通过识别代表性样本，即聚类中心，及相应的子集对数据进行聚类。算法输入为成对数据之间的相似性，通过数据之间的消息迭代交换，可得到一组高质量的代表性样本和相应的子集。相比于K均值等聚类技术，AP聚类不需要提前指定聚类团的个数，聚类精度较高，算法复杂度较低。

受AP聚类的启发，本文将AP聚类扩展到多路径聚类约束模型，并提出一种新的多路径AP聚类算法。首先，将量测通过所有传播路径变换到目标坐标系形成多个路径相关量测，其中，通过正确传播路径得到的路径相关量测为真实路径相关量测，否则为虚假路径相关量测。其次，聚类团被建模为可以自动识别数目的代表性路径相关量测，并将每个真实路径相关量测分配给最合适的代表性路径相关量测，多路径聚类约束通过两个方面加到模型中，一是同一个量测和同一种传播路径得到的路径相关量测不能互为代表性路径相关量测，二是同一种传播路径产生的路径相关量测不能选择同一个代表性路径相关量测。最后，构建聚类变量联合概率密度函数及其对应的因子图，本模型的求解是非确定性多项式难题(Non-deterministic Polynomial-hard, NP-hard)^[10]，采用最大和置信传播算法求取近似解，并提取量测的传播路径以及量测聚类团。本文主要贡献如下：

(1) 将原始AP聚类算法拓展到多路径数据模型，并提出一种新的多路径AP聚类算法；

(2) 多路径AP聚类算法继承了原始AP聚类算法的优势，能够自动估计聚类团数目；单次消息传播的时间复杂度为量测与路径个数乘积的平方。

本文其余部分组织如下：第2节描述多路径量测模型及聚类问题；第3节介绍多路径AP聚类算法及其实现流程；第4节和第5节分别给出仿真实验分析和结束语。

2.   系统建模与问题描述

2.1   多路径量测模型

考虑天波雷达目标探测场景，在某一时刻，${{\boldsymbol{x}}_k} \triangleq {[{x_k},{y_k}]^{\rm T}} \in \mathcal{G}$表示第$k$个目标的状态，其中，$\mathcal{G}$为目标坐标系，类似文献[7]，$\mathcal{G}$选择麦卡托坐标系^[21]，${x_k}$和${y_k}$分别为该目标在$\mathcal{G}$上的横轴和纵轴坐标值，$k = 1,2, \cdots ,K$，$K$为目标个数。目前，天波雷达的目标高度测量误差大，因此本文假设目标高度信息缺失。${\boldsymbol{{z}}_i} \triangleq {[{r_i},{\zeta _i}]^{\rm T}} \in \mathcal{R}$，表示雷达的第$i$个量测，其中，$\mathcal{R}$为雷达极坐标系，${r_i}$和${\zeta _i}$分别为该量测在$\mathcal{R}$上的径向距离和方位角^[7]，$i \in \mathcal{M} \triangleq \{ 1, 2, \cdots ,M\}$，$M$为量测个数。对于多路径量测系统，量测${z_i}$可能通过特定的传播路径由目标反射得到，也有可能来自杂波，可以表示为

其中，${h_\tau }( \cdot )$,$\tau \in \mathcal{P} \triangleq \{ 1,2, \cdots ,P\}$，为传播路径$\tau$的状态量测转换关系，$P$为已知的传播路径数，${{\boldsymbol{w}}_\tau }$为量测噪声，服从协方差为${{\boldsymbol{R}}_\tau }$的0均值高斯分布。目标检测概率为$P_{\text{d}}^\tau$，且对不同的传播路径相互独立。杂波在雷达探测区域内服从均匀分布，每次扫描的杂波个数服从泊松分布，其中，雷达探测区域体积记为$V$，杂波密度为${f_{{\text{FA}}}}$，泊松分布的均值为${\lambda _{{\text{FA}}}} = {f_{{\text{FA}}}}V$。

假设电离层状态$\mathcal{I}$已知，对于传播路径$\tau$，雷达极坐标系到麦卡托坐标系的坐标配准可以通过非线性变换实现${T_\tau }:\mathcal{R} \times \mathcal{I} \to \mathcal{G}$^[7]。量测${{\boldsymbol{z}}_i}$通过${T_\tau }$坐标配准到麦卡托坐标系得到路径相关量测$\tilde z_i^\tau \in \mathcal{G}$，其中，通过正确传播路径得到的路径相关量测为真实路径相关量测，否则为虚假路径相关量测。为了同时满足精度和计算要求，采用无迹变换^[22]对${T_\tau }$进行线性化。因此，路径相关量测$\tilde z_i^\tau$的状态估计近似服从高斯分布。量测集合定义为${\boldsymbol{Z}} \triangleq {[{\boldsymbol{z}}_1^{\rm T}{\boldsymbol{z}}_2^{\rm T} \cdots {\boldsymbol{z}}_M^{\rm T}]^{\rm T}}$，量测${{\boldsymbol{z}}_i}$的路径相关量测集合定义为${\tilde {\boldsymbol{Z}}_i} \triangleq [{(\tilde {\boldsymbol{z}}_i^1)^{\rm T}}{(\tilde {\boldsymbol{z}}_i^2)^{\rm T}} \cdots {(\tilde {\boldsymbol{z}}_i^P)^{\rm T}}]^{\rm T}$，全体路径相关量测集合定义为$\tilde {\boldsymbol{Z}} \triangleq {[\tilde {\boldsymbol{Z}}_1^{\rm T}\tilde {\boldsymbol{Z}}_2^{\rm T} \cdots \tilde {\boldsymbol{Z}}_M^{\rm T}]^{\rm T}}$。

2.2   多路径量测聚类实例

下面结合一个具体实例，阐述多路径量测聚类问题。假设雷达探测场景内存在2个目标，电离层存在4种传播路径。图1(a)展示了天波雷达多路径量测模型，即单个目标可以通过多种传播路径在雷达极坐标系下产生多个量测。图1(b)展示了雷达极坐标系量测坐标配准过程，由于量测的传播路径未知，将其通过全部传播路径配准到目标坐标系得到路径相关量测，其中填充图形表示真实路径相关量测，未填充图形表示虚假路径相关量测。

图 1  多路径量测与坐标配准示意图

下载: 全尺寸图片 幻灯片

图2展示了雷达极坐标系下的多路径量测与杂波，目标1被传播路径1,2,4检测到，未被传播路径3检测到，目标2被全部传播路径检测到，另外还检测到6个杂波，可以看到同一个目标的多路径量测在雷达极坐标系下不具备空间相邻特性。图3展示了目标坐标系下的路径相关量测图，红色虚线框内为同一个目标的全部真实路径相关量测，红色虚线框外为虚假路径相关量测和杂波产生的路径相关量测，其中，黑色虚线框内表示可能的干扰聚类团。从图3可以看到同一个目标的真实路径相关量测具有空间相邻特性，其可作为传播路径辨识和聚类的特征。

图 2  多路径量测与杂波图

下载: 全尺寸图片 幻灯片

图 3  路径相关量测图

下载: 全尺寸图片 幻灯片

若采用MH思想枚举所有可能的聚类假设^[4,5]，当量测个数为$M$、传播路径为$P$时，场景中可能存在的目标数目为^[23]

其中，$C_M^i$表示从$M$个量测中取出$i$个元素的所有组合的个数，$A_P^i$表示从$P$个传播路径中取出$i$个元素的所有排列的个数。可能存在的聚类假设数目具有如式(3)的递归形式^[5]

其中，$H(0,P) = 1$。当量测个数$M = 10$，传播路径$P = 9$时，可能存在的目标个数为$N(10,9) \approx 6 \times {10^7}$，可能存在的聚类假设数目为$H(10,9) \approx 4.89 \times {10^{12}}$。显然，枚举算法在实际使用中无法满足实时性和存储要求。

2.3   多路径聚类模型

多路径聚类变量定义为${\boldsymbol{C}} = \{ c_{i,\tau }^{m,t}:i \in \mathcal{M}, \tau \in \mathcal{P},m \in \mathcal{M}\backslash \{ i\} ,t \in \mathcal{P}\backslash \{ \tau \} \} \cup \{ c_{m,t}^{m,t}:m \in \mathcal{M},t \in \mathcal{P}\}$, $c_{i,\tau }^{m,t} \in \{ 0,1\}$, $i \ne m,\tau \ne t$，当$\tilde z_m^t$是$\tilde z_i^\tau$的聚类中心时，$c_{i,\tau }^{m,t} = 1$，否则，$c_{i,\tau }^{m,t} = 0$，即

以及$c_{m,t}^{m,t} \in \{ 0,1\}$，当$\tilde z_m^t$为聚类中心时，$c_{m,t}^{m,t} = 1$，否则$c_{m,t}^{m,t} = 0$，即

全体路径相关量测$\tilde {\boldsymbol{Z}}$对聚类变量${\boldsymbol{C}}$的统计依赖可以表示为全局似然函数$p(\tilde {\boldsymbol{Z}}|{\boldsymbol{C}})$。根据一般假设^[6,7]，$\tilde {\boldsymbol{Z}}$对不同的量测和路径独立，全局似然函数$p(\tilde {\boldsymbol{Z}}|{\boldsymbol{C}})$可以因子分解为

其中，$c_{i,\tau }^{m,t} = 1$表示路径相关量测$\tilde {\boldsymbol{z}}_i^\tau$和$\tilde {\boldsymbol{z}}_m^t$来自同一个目标，由于$\tilde {\boldsymbol{z}}_i^\tau$和$\tilde {\boldsymbol{z}}_m^t$分别服从高斯分布$N({\widehat{{\boldsymbol{z}}}}_{i}^{\tau },{{\boldsymbol{P}}}_{i}^{\tau })$以及$N(\hat {\boldsymbol{z}}_m^t,{\boldsymbol{P}}_m^t)$，可得$\tilde {\boldsymbol{z}}_i^\tau - \tilde {\boldsymbol{z}}_m^t$也服从高斯分布$N(0, {\boldsymbol{P}}_i^\tau + {\boldsymbol{P}}_m^t)$，可得

$c_{i,\tau }^{m,t} = 0$表示路径相关量测$\tilde z_i^\tau$和$\tilde z_m^t$不是来自同一个目标，因此$p(\tilde z_i^\tau ,\tilde z_m^t|c_{i,\tau }^{m,t} = 0)$可以定义为${f_{{\text{FA}}}}$。

为了方便描述，定义如下多路径量测聚类变量集合${\mathcal{C}_i} \triangleq \{ c_{i,\tau }^{m,t}:\;\tau \in \mathcal{P},m \in \mathcal{M}\backslash \{ i\} ,\;t \in \mathcal{P}\backslash \{ \tau \} \} \cup$$\{ c_{i,\tau }^{i,\tau }: \tau \in \mathcal{P}\}$，${\mathcal{C}^m} \triangleq \{ c_{i,\tau }^{m,t}:t \in \mathcal{P},i \in \mathcal{M}\backslash \{ m\} ,\tau \in \mathcal{P}\backslash \{ t\} \} \cup \{ c_{m,t}^{m,t}:t \in \mathcal{P}\}$。根据多路径量测模型及聚类模型，假设合理的聚类矩阵满足下面约束：

(1) 约束1 1个量测只能选择1种传播路径，且只能选择1个路径相关量测作为其聚类中心，表示为${I_i}({\mathcal{C}_i})$，$i \in \mathcal{M}$；

(2) 约束2 如果存在量测1来源于路径1并选择某个来源于路径2的量测2作为其聚类中心，则量测2只能来源于路径2且必须选择自己作为聚类中心，并且量测2不能再作为路径1下其他量测的聚类中心，表示为${E_m}({\mathcal{C}^m})$，$m \in \mathcal{M}$。

由上述两个约束可以得到约束方程组

其中，$|$表示逻辑或。聚类变量${\boldsymbol{C}}$的全体约束可以定义为

满足约束$\varOmega ({{\boldsymbol{C}}})$的聚类变量集合记为$\mathcal{C}$。

多路径量测聚类的先验概率可以定义为

其中，$p(c_{i,\tau }^{m,t}) = 1,i \ne m,\tau \ne t$，以及

其中，$\lambda$为偏向参数，值越大产生的聚类团个数越多，反之亦然。

根据贝叶斯公式，聚类变量${\boldsymbol{C}}$的后验概率分布正比于联合概率密度函数，即可以表示为

将式(6)、式(10)和式(11)代入式(13)可得

由于$\tilde {\boldsymbol{z}}_m^t$和$\tilde {\boldsymbol{z}}_i^\tau$已知，$p({\tilde {\boldsymbol{z}}}_i^\tau ,{\tilde {\boldsymbol{z}}}_m^t|c_{i,\tau }^{m,t})p(c_{i,\tau }^{m,t})$和$p(c_{m,t}^{m,t})$可以分别表示为证据因子$s(c_{i,\tau }^{m,t})$和$s(c_{m,t}^{m,t})$，式(14)可以表示为

对于$s(c_{i,\tau }^{m,t})$和$s(c_{m,t}^{m,t})$，定义对数似然比

其中，${\text{const}}$为与$c_{i,\tau }^{m,t}$无关的常数，以及

由式(16)可以看出，$\tilde s_{i,\tau }^{m,t}$正比于$\tilde {\boldsymbol{z}}_m^t$和$\tilde {\boldsymbol{z}}_i^\tau$距离的负值，表示$\tilde {\boldsymbol{z}}_m^t$和$\tilde {\boldsymbol{z}}_i^\tau$的相似度，可以看出聚类特征为同一个目标的多路径量测在目标坐标系中的位置相邻特性。由式(17)可以看出，$\tilde s_{m,t}^{m,t}$等于偏向因子$\lambda$，通过设置合理的偏向参数可以实现自动识别聚类团数目。

2.4   问题描述

天波雷达多路径量测聚类的目的是根据路径相关量测集合$\tilde {\boldsymbol{Z}}$估计聚类变量${\boldsymbol{C}}$，在贝叶斯框架下该问题可以转化为最大后验概率估计问题，即

通过精确搜索最优聚类变量使得目标函数式(18)最大化是NP-hard^[10]，需要采用近似算法求解。

3.   多路径约束AP聚类算法

3.1   消息传播方法

式(15)表示的联合后验概率密度函数可以表示为如图4所示的因子图，其中黑色实线圆圈表示变量节点，蓝色方块${I_i} \triangleq {I_i}({\mathcal{C}_i})$和橙色方块${E_m} \triangleq {E_m}({\mathcal{C}^m})$均为因子节点。需要注意的是，每一个变量节点都与对应的证据因子节点连接，为了图示简洁，证据因子未画出。为体现因子图结构，图4还绘制了不满足多路径约束的聚类变量，如红色虚线圆圈所示。

图 4  因子图和不满足约束的聚类变量

下载: 全尺寸图片 幻灯片

联合后验概率密度$p({\boldsymbol{C}}|\tilde {\boldsymbol{Z}})$的近似最大后验概率估计$\tilde f({\boldsymbol{C}})$可以通过在因子图上迭代运行最大和置信传播算法得到。最大和置信传播算法在因子图上进行两类消息传播，即从变量节点到每个相邻因子节点的消息，以及从因子节点到每个相邻变量节点的消息。变量节点到因子节点的消息为其他相邻因子节点到该变量节点的消息之和，因子节点到变量节点的消息为因子节点的对数与其他变量节点到因子节点的消息之和对其他变量节点取最大化，通过这两类消息传递可以保证消息收敛到邻域最大值^[24]。根据最大和置信传播算法，变量节点和因子节点之间传递5种类型的消息，如图5所示，其中$s_{i,\tau }^{m,t} \triangleq s(c_{i,\tau }^{m,t})$

图 5  消息传播因子图

下载: 全尺寸图片 幻灯片

通过应用文献[10,11]中使用的一些数学技巧，可以得到简化消息为

其中，$\zeta _{m,t}^{i,\tau }$, $\phi _{m,t}^{i,\tau }$, $\xi _{m,t}^\tau$, ${\psi _{m,t}}$的定义为

以及，消息$\theta _{i,\tau }^{m,t}(c_{i,\tau }^{m,t})$,$\beta _{t,\tau }^{m,t}(c_{t,\tau }^{m,t})$和$\eta _{i,\tau }^{m,t}(c_{i,\tau }^{m,t})$被隐含在上述简化消息中，因此被省略了。从式(24)—式(26)的消息传播公式可以看出，聚类变量之间存在两类消息，分别为$\tilde {\boldsymbol z}_m^t$对$\tilde {\boldsymbol z}_i^\tau$的吸引度$\tilde \rho _{i,\tau }^{m,t}$和$\tilde {\boldsymbol z}_i^\tau$对$\tilde {\boldsymbol z}_m^t$的归属度$\tilde \alpha _{i,\tau }^{m,t}$。吸引度$\tilde \rho _{i,\tau }^{m,t}$如式(24)，$\tilde s_{i,\tau }^{m,t}$表示$\tilde {\boldsymbol z}_m^t$与$\tilde {\boldsymbol z}_i^\tau$的相似度，$\max [\tilde s_{i,\tau '}^{m',t'} + \tilde \alpha _{i,\tau '}^{m',t'}]$表示其他路径相关量测$\tilde {\boldsymbol z}_{m'}^{t'}$对$\tilde {\boldsymbol z}_i^\tau$的相似度和归属度之和的最大值，因此$\tilde s_{i,\tau }^{m,t} - \max [\tilde s_{i,\tau '}^{m',t'} + \tilde \alpha _{i,\tau '}^{m',t'}]$可以表示$\tilde {\boldsymbol z}_m^t$对$\tilde {\boldsymbol z}_i^\tau$的吸引度。$\tilde \alpha _{m,t}^{m,t}$表示$\tilde {\boldsymbol z}_m^t$对自身的归属度，如式(25)，${\psi _{m,t}}$为$\tilde {\boldsymbol z}_m^t$接收到的正吸引度之和，$\max [\tilde \rho _{m,t'}^{m,t'} + {\psi _{m,t'}}]$表示${z_m}$的其他路径相关量测$\tilde {\boldsymbol z}_m^{t'}$接收到的正吸引度之和的最大值。$\tilde \alpha _{i,\tau }^{m,t}$表示$\tilde {\boldsymbol z}_i^\tau$对$\tilde {\boldsymbol z}_m^t$的归属度，如式(26)，$\tilde \rho _{m,t}^{m,t}$表示$\tilde {\boldsymbol z}_m^t$对自身的吸引度，$\phi _{m,t}^{i,\tau }$表示$\tilde {\boldsymbol z}_m^t$接收到除量测${z_i}$和传播路径$\tau$以外的正吸引度之和，$\tilde \rho _{m,t}^{m,t} + \zeta _{m,t}^{i,\tau } + \xi _{m,t}^\tau$表示$\tilde {\boldsymbol z}_m^t$为聚类中心且$c_{i,\tau }^{m,t} = 0$时$\tilde {\boldsymbol z}_m^t$接受到的正吸引度之和，$\max [\tilde \rho _{m,t'}^{m,t'} + {\psi _{m,t'}}]$表示$\tilde {\boldsymbol z}_m^t$不为聚类中心时${z_m}$的其他路径相关量测$\tilde {\boldsymbol z}_m^{t'}$接受到了正吸引度之和的最大值。因此，$\tilde {\boldsymbol z}_m^t$在$c_{i,\tau }^{m,t} = 1$时接受到的正吸引度之和越大，$\tilde {\boldsymbol z}_m^t$在$c_{i,\tau }^{m,t} = 0$时或$\tilde {\boldsymbol z}_m^{t'},t' \ne t$接收到的正吸引度之和越小，$\tilde {\boldsymbol z}_i^\tau$对$\tilde {\boldsymbol z}_m^t$的归属度越大。吸引度$\tilde \rho _{i,\tau }^{m,t}$和归属度$\tilde \alpha _{i,\tau }^{m,t}$都初始化为0，迭代终止条件为相邻两次迭代消息变化小于终止阈值${\delta _l}$或达到最大迭代次数${\iota _{{\text{max}}}}$。

3.2   聚类团提取与量测融合

可以根据消息更新结果提取量测的传播路径以及同一个目标的量测聚类团。第1步，提取量测的传播路径。为了估计量测${z_i}$的传播路径，将所有传向量测${{\boldsymbol{z}}_i}$的消息相加，并取使之最大化的传播路径${\hat \tau ^i}$作为其传播路径，即

第2步，在真实路径相关量测$\tilde {\boldsymbol{z}}_i^{{{\hat \tau }^i}},i = 1,2, \cdots ,M$，中提取同一个目标的量测聚类团。为了估计聚类变量的值，将所有传入聚类变量的消息相加，并取使之最大化的值，即

然后，通过聚类变量的估计值提取出聚类中心及其对应的多路径量测。最后，将聚类团内多路径量测通过简单凸组合算法^[25]得到目标的量测融合状态估计。

3.3   算法流程

综上所述，多路径AP聚类算法流程如下：

步骤1　多假设坐标配准。将量测通过所有传播路径配准到目标坐标系得到路径相关量测。

步骤2　证据因子初始化。根据式(16)和式(17)初始化证据因子。

步骤3　消息迭代更新。根据式(24)—式(26)进行消息迭代更新，直至达到迭代终止条件。

步骤4　聚类团提取。根据消息更新结果提取量测的传播路径以及同一个目标的量测聚类团。

步骤5　量测融合。根据聚类团内多路径量测融合得到目标的融合状态估计。

3.4   算法特性分析

步骤1需要$MP$次坐标变换，时间复杂度为$O(MP)$。步骤2初始化证据因子消息最少需要$MP(MP - 1)/2$次距离计算，时间复杂度为$O({M^2}{P^2})$。步骤3的计算消耗主要来自式(24)—式(26)的迭代消息更新，由于各个消息中存在大量重复部分，这些重复部分只运算1次可以进一步降低计算量，总的时间复杂度为$O({\iota _{{\text{max}}}}{M^2}{P^2})$，${\iota _{{\text{max}}}}$为最大迭代次数。步骤4聚类团提取时需要${M^2}{P^2}$次求和以及$2MP$次比较，时间复杂度为$O({M^2}{P^2})$。步骤5至多进行$M$次量测融合，每次融合至多有$P$个量测，因此时间复杂度至多为$O(MP)$。综上所述，多路径AP聚类算法的时间复杂度为$O({\iota _{{\text{max}}}}{M^2}{P^2})$。

综上可得，多路径AP聚类算法可以看作AP聚类算法在多路径数据模型下的推广。两个算法都有3个消息传递公式，多路径AP聚类的消息传递过程中同时考虑了多路径约束。与AP聚类算法一致，通过合理的偏向参数$\lambda$，多路径AP聚类算法能够自动估计聚类团数目。

4.   仿真实验

4.1   仿真场景

考虑天波雷达多目标探测场景，如图6所示。仿真场景参数设置如下：电离层为两层结构，即E层和F层，高度分别为${h_{\rm{E}}} = 100$ km, ${h_{\rm{F}}} = 260$ km，传播路径数为$P = 4$；目标个数$T = 8$；雷达接收机经纬度分别为东经143.20°和南纬24.29°，法向角为$325{{\text{ }}^ \circ }$，距离发射机100 km；径向距离、方位角和径向速度的范围分别为$[1100,1800]$ km, $[0.1,0.7]$ rad和$[ - 0.3,0.3]$ km/s；量测噪声协方差为${R_\tau } = {\text{diag}}(25{\text{ k}}{{\text{m}}^2}$, $9 \times {10^{ - 6}}{\text{ ra}}{{\text{d}}^2})$。偏向参数设置为$\lambda = {10^{ - 7}}$，可以实现聚类团数目的自动识别。消息传递迭代终止阈值为${\delta _l}{ = 10^{ - 6}}$，最大迭代次数为${\iota _{{\text{max}}}} = 1000$。

图 6  麦卡托坐标系下目标位置图

下载: 全尺寸图片 幻灯片

天波雷达量测聚类是为了得到航迹起始后目标统一的态势，因此通过目标检测与定位评估聚类算法的性能。性能评估之前需要找到检测目标和真实目标一一对应的匹配关系，$x$轴和$y$轴的匹配阈值都为$10$ km。评估指标包括：(1)检测正确率($\uparrow$)：可以与真实目标匹配的检测目标数与检测目标总数之比；(2)漏检率($\downarrow$)：未与检测目标匹配的真实目标数与真实目标总数之比；(3)均方根误差($\downarrow$)：可以匹配上的检测目标与真实目标位置的均方根误差；(4)最优次子模式分配(Optimal Sub Pattern Assignment, OSPA $\downarrow$)^[26]；(5)运行时间($\downarrow$)：算法运行时间。$\uparrow$表示指标的值越大性能越好，$\downarrow$反之。OSPA参数设置如下，$p = 2$和$c = 19.08{\text{ km}}$。蒙特卡罗仿真次数为4000次。蒙特卡罗仿真时，对指标在蒙特卡罗次数上求平均。对比算法为多路径MH聚类算法，最佳假设树的保留个数分别取为$N = 200,400,600,800$。所有算法的仿真实验都是在同一台个人计算机(联想，Intel(R) Core(TM) i7-9700 3.00 GHz处理器，16 GB内存和Windows 10操作系统)上利用MATLAB R2020a仿真环境完成的。

4.2   仿真结果与分析

图7展示了$P_{\text{d}}^\tau = 0.75$, $\tau = 1,2,3,4$, ${f_{{\text{FA}}}} = 2 \times {10^{ - 5}}$时，8个目标的多路径量测以及杂波位置分布。

图 7  多路径量测和杂波图

下载: 全尺寸图片 幻灯片

图8(a)和图8(b)分别表示多路径AP聚类算法和MH算法在$K = 800$时的聚类结果，可以看出多路径AP聚类算法成功检测出全部目标，且无错误检测，MH算法漏检了目标${T_1}$, ${T_3}$, ${T_6}$，产生了3个误检。

图 8  融合结果

下载: 全尺寸图片 幻灯片

图9(a)—图9(e)展示了算法性能指标随目标个数的变化，其中图9(e)中实线(对应图例中包含Left)的坐标轴在左侧，虚线(对应图例中包含Right)的坐标轴在右侧，曲线AP(Left)和AP(Right)分别表示大刻度和小刻度坐标系下多路径AP聚类算法的运行时间。总体来看，随着目标个数的增加，多路径AP聚类算法的各项性能指标变化较小，且都明显优于MH算法，这是由于该算法通过局部量测之间的消息传递可以得到邻域最优聚类解，且受目标个数影响较小；而MH聚类算法性能明显变差，这是由于随着目标个数增加，如表1所示，MH聚类算法需要枚举的假设个数呈指数式增加，由于采用了剪枝策略，使得MH聚类算法对假设树空间的搜索能力降低，通过增加保留假设树的数目，虽然聚类精度有所改善，但运行时间急剧上升。需要注意的是，多路径AP聚类算法的检测正确率和均方根误差随目标个数的增加都呈上升趋势，这是由于杂波密度确定时，误检测数基本确定，随着目标个数的增加，正确检测数会增加，但是杂波与目标真实量测聚类的情况也会增加，导致融合后的定位均方根误差增加。图9(c)和图9(d)还展示了量测通过真实聚类结果融合时的最优均方根误差和OSPA随目标个数变化，可以看出，多路径AP聚类算法与最优结果相差较小，可以说明所提算法具有较高的近似精度。在算法运行时间和内存消耗方面，分别如图9(e)和表2所示，可以看出多路径AP聚类与MH聚类算法的运行时间和内存消耗都随目标个数而增加。同时，MH聚类算法的运行时间和内存消耗随保留假设树数目的增加而增加。对比两种算法，MH聚类算法所需的运行时间远高于多路径AP聚类算法，内存消耗略高于多路径AP聚类算法，这是由于在编程实现中MH算法对每一个量测序贯扩展假设树并进行剪枝操作，剪枝操作后程序会自动释放被剪枝假设树的内存，使得MH算法的运行时间较长，而内存消耗在有限范围内。此外，当目标个数为8时，多路径AP聚类算法单次运行时间为0.24 s、内存消耗在100 kB左右，满足实时性和计算资源要求。

图 9  性能指标随目标个数变化图

下载: 全尺寸图片 幻灯片

表 1  假设个数随着目标个数变化表

目标数	1	2	3	4	5	6	7	8
假设数	3e6	2e7	9e7	3e8	7e8	2e9	3e9	6e9

下载: 导出CSV 

| 显示表格

表 2  内存消耗(kB)随目标个数变化表

算法	目标数
	1	2	3	4	5	6	7	8
AP	70.8	71.50	70.5	77.2	92.2	96.8	98.8	99.7
MH(N=200)	77.1	72.2	78.1	103.5	102.2	97.6	100.0	102.7
MH(N=400)	83.3	80.9	78.4	106.9	116.0	110.5	104.3	105.3
MH(N=600)	81.1	102.1	125.3	94.7	145.7	95.5	133.3	128.5
MH(N=800)	92.3	107.5	119.3	117.0	130.4	139.1	131.1	136.3

下载: 导出CSV 

| 显示表格

5.   结束语

本文针对天波雷达多路径量测聚类问题，将AP聚类扩展到多路径约束模型，并提出一种新的多路径AP聚类算法。该算法可以同时对多路径量测进行电离层传播路径辨识和聚类，能够自动估计聚类团数目，单次消息传播的时间复杂度为量测与传播路径个数乘积的平方，可用于天波雷达量测级目标跟踪的航迹快速起始。仿真实验分析表明，所提算法较多路径多假设聚类算法具有更好的聚类性能。对于后续研究工作，一方面，如果未来天波雷达能够提供更加准确的高度信息，从而降低定位误差，可进一步改善本文提出的多路径AP聚类算法性能；另一方面，可考虑将AP聚类扩展到天波雷达组网量测聚类，解决天波雷达组网量测级跟踪航迹快速起始问题。