一种天波雷达多路径量测聚类方法与流程

1.本发明属于雷达数据处理领域中的天波雷达多路径量测聚类子系统技术领域，特别涉及一种天波雷达多路径量测聚类方法。


背景技术：

2.天波雷达利用电磁波在电离层与地/海面之间的反射作用传输高频能量，可实现对高价值目标的远程预警。由于目标与雷达站之间存在多条电离层传播路径，单目标在雷达射线坐标系下会产生多路径量测，进而可能导致多路径航迹，使得态势感知和情报出现误判。因此，多路径量测为天波雷达数据处理带来严峻挑战。
3.现有的天波雷达数据处理主要集中在目标跟踪，其按照多路径量测融合阶段的不同可以分为航迹级跟踪融合和量测级跟踪融合。航迹级跟踪融合将问题分解为单路径目标跟踪和多路径航迹融合两个子问题。量测级跟踪融合将单路径数据关联扩展到多路径数据关联，通过目标-量测-路径关联，在量测级融合目标信息，跟踪性能较航迹级跟踪算法具有明显提升。量测级跟踪算法在预处理阶段需要进行航迹起始，即将来自同一目标的不同路径量测进行模式辨识，并转化到目标坐标系进行聚类融合，得到目标的统一态势，并初始化航迹状态。错误的航迹起始可能导致航迹冗余或目标漏跟。然而，目前的量测级跟踪算法主要关注跟踪阶段，很少研究航迹起始。
4.天波雷达航迹起始本质上需要解决多路径量测聚类问题，即对同一目标的多路径量测进行模式辨识和聚类融合。因此，如何发展出一种天波雷达多路径量测聚类算法，在天波雷达数据处理领域具有一定的理论及实际工程意义。


技术实现要素：

5.本发明解决的技术问题是：为了解决现有技术的不足之处，提高天波雷达多路径量测聚类性能，同时满足研究所实际需求，本发明提出一种基于相似性传播的多路径量测聚类算法。本发明的优点是能够自动估计聚类团数目；单次消息传播的计算复杂度为量测与路径数乘积的平方。
6.本发明的技术方案是：步骤1：建立模型并对模型初始化，包括以下子步骤：
7.步骤1.1：多路径量测建模：通过所有电离层传播路径将量测坐标配准得到路径相关量测
8.步骤1.2：聚类建模；
9.步骤1.3：计算先验概率：根据局部传感器的检测概率以及杂波密度，计算聚类变量的初始置信p(c)，并根据多路径聚类约束得到聚类变量约束；
10.步骤1.4：计算似然函数：根据路径相关量测对聚类变量c的统计依赖，计算全局似然函数
11.步骤2：采用最大和置信传播算法进行迭代优化，包括以下子步骤：
12.步骤2.1：因子图构建：根据后验概率的因子分解式构建航迹关联因子图；
13.步骤2.2：消息初始：根据步骤1计算的先验概率p(c)和全局似然函数对因子图中的证据因子节点进行消息初始；
14.步骤2.3：消息迭代：在聚类变量和约束1因子节点和约束2因子节点之间进行消息迭代传播，当消息收敛或达到最大迭代次数时终止迭代；
15.步骤2.4：聚类团提取：根据聚类变量的最大后验概率估计提取量测的传播路径以及同一个目标的量测聚类团；
16.步骤3：航迹融合：将聚类团内多路径量测通过简单凸组合算法，得到目标的量测融合状态估计。
17.本发明进一步的技术方案是：所述步骤1.1中，在天波雷达目标探测场景的某一时刻下，k＝1,...,t，表示责任区内第k个目标的状态，t为目标个数。i＝1,...,m，表示第i个量测，m为量测个数；电离层传播路径数为p，假设电离层状态已知，对于传播路径τ，雷达极坐标系到麦卡托坐标系的坐标配准可以通过非线性变换实现t
τ
:量测集合定义为量测yi的路径相关量测集合定义为全体路径相关量测集合定义为假设目标检测概率为杂波在雷达探测区域内服从均匀分布，每次扫描的杂波个数服从泊松分布。
18.本发明进一步的技术方案是：所述步骤1.2中，定义多路径聚类变量集合满足以下约束：
19.约束1：一个量测只能通过一种电离层传播路径得到，一个目标通过一种电离层传播路径至多产生一个量测，表示为i,m＝1,...,m，t,τ＝1,...,p；
20.约束2：一个量测只能选择一个聚类中心，表示为i＝1,...,m；
21.约束3：如果存在量测选择某个量测作为聚类中心，则该量测必须选择自己作为聚类中心，表示为m＝1,...,m。
22.满足约束1-3的聚类矩阵集合记为
23.路径相关量测对聚类变量c的统计依赖可以表示为全局似然函数根据一般假设，对不同的量测和路径独立，全局似然函数可以因子分解为
[0024][0025]
本发明进一步的技术方案是：所述先验概率模型为
[0026]
本发明进一步的技术方案是：所述步骤2.1中的关联因子图为：
[0027][0028]
本发明进一步的技术方案是：所述步骤2中的聚类团提取包括以下子步骤：
[0029]
子步骤1：提取量测的传播路径：为了估计量测yi的传播路径，将所有传向量测yi的消息相加，并取使之最大化的路径，即
[0030][0031]
量测yi的传播路径为
[0032]
子步骤2：在正确路径相关量测i＝1,...,m，中提取同一个目标的量测聚类团；为了估计聚类变量的值，将所有传入聚类变量的消息相加，并取使之最大化的值，即
[0033][0034]
然后通过聚类变量的估计值提取出聚类中心及其子集；
[0035]
子步骤3：将聚类团内多路径量测通过简单凸组合算法得到目标的量测融合状态估计。
[0036]
发明效果
[0037]
本发明的技术效果在于：本发明涉及一种基于相似性传播的天波雷达多路径量测聚类方法，可以同时实现多路径量测的模式辨识以及量测聚类。首先，将相似性传播聚类模型扩展到多路径聚类，并构建多路径聚类约束。然后，构建聚类变量联合概率密度函数及其对应的因子图。该模型的求解是np-hard问题，采用最大和置信传播算法求取近似解。多路径相似性传播聚类算法能够自动估计聚类团数目；使用最大和置信传播算法近似计算聚类变量的最大后验概率，避免了多路径量测聚类假设的枚举，单次消息传播的计算复杂度为量测与路径数乘积的平方。仿真结果表明，多路径相似性传播聚类算法较多假设聚类算法相比，在检测正确率、漏检率、均方根误差、ospa和运行时间性能指标相当情况下可以大大降低运算时间。
附图说明
[0038]
图1为多路径量测及量测聚类问题图，子图(a)为多路径量测与多假设坐标配准示意图，子图(b)为多路径量测与杂波图，子图(c)为路径相关量测图。
[0039]
图2为多路径量测聚类因子图。
[0040]
图3为目标平行运动融合结果图。
[0041]
图4为麦卡托坐标系下目标位置图。
[0042]
图5为多路径量测和杂波图，其中
○
(红色)表示多路径量测，
·
(黑色)表示杂波
[0043]
图6为融合结果图，子图(a)为相似性传播聚类算法融合结果图，子图(b)为多假设算法融合结果图。
·
(黑色)、
○
(红色)、*(蓝色)分别表示真实目标、正确检测和错误检测
[0044]
图7为性能指标随目标个数变化图，子图(a)-(e)分别为检测正确率、漏检率、均方根误差、ospa和运行时间随目标个数变化图。
具体实施方式
[0045]
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
[0046]
参见图1-图7，本发明的技术方案是：一种天波雷达多路径量测聚类方法，包括以下步骤：
[0047]
步骤一：系统建模与问题描述，包括以下子步骤：
[0048]
子步骤一：多路径量测模型：
[0049]
考虑天波雷达目标探测场景，在某一时刻，k＝1,...,t，表示责任区内第k个目标的状态，t为目标个数。i＝1,...,m，表示第i个量测，m为量测个数。对于多路径量测系统，量测yi可能通过特定的传播路径由目标反射得到，也有可能来自杂波。电离层传播路径数为p，假设电离层状态已知，对于传播路径τ，雷达极坐标系到麦卡托坐标系的坐标配准可以通过非线性变换实现t
τ
:量测集合定义为量测yi的路径相关量测集合定义为全体路径相关量测集合定义为假设目标检测概率为杂波在雷达探测区域内服从均匀分布，每次扫描的杂波个数服从泊松分布。
[0050]
子步骤二：多路径量测聚类模型：
[0051]
对于路径相关量测和考虑二值聚类变量若是的聚类中心，则否则令为联合聚类变量。为了方便描述，定义多路径聚类变量集合根据多路径量测模型及聚类模型，假设合理的聚类矩阵满足下面约束：
[0052]
约束1：一个量测只能通过一种电离层传播路径得到，一个目标通过一种电离层传播路径至多产生一个量测，表示为i,m＝1,...,m，t,τ＝1,...,p；
[0053]
约束2：一个量测只能选择一个聚类中心，表示为i＝1,...,m；
[0054]
约束3：如果存在量测选择某个量测作为聚类中心，则该量测必须选择自己作为聚类中心，表示为m＝1,...,m。
[0055]
满足约束1-3的聚类矩阵集合记为
[0056]
路径相关量测对聚类变量c的统计依赖可以表示为全局似然函数根据一般假设，对不同的量测和路径独立，全局似然函数可以因子分解为
[0057][0058]
子步骤三：聚类变量的先验概率模型：
[0059]
假设给定目标，聚类变量对不同量测和不同路径是独立的，聚类变量联合先验概率密度函数可以分解为
[0060][0061]
子步骤四：问题描述：
[0062]
多路径量测聚类的目的是根据路径相关量测集合估计聚类变量c，在贝叶斯框架下该问题可以转化为最大后验概率估计问题，即
[0063][0064]
步骤二：基于相似性传播的多路径量测聚类，包括以下子步骤：
[0065]
子步骤一：根据后验概率f(a|z)的因子分解式构建航迹关联因子图：
[0066]
使用贝叶斯公式可以得到以下公式
[0067][0068]
上式表示的联合后验概率密度函数可以表示为因子图，变量节点为聚类变量，因子节点包括证据因子节点约束2因子节点和约束3因子节点通过在因子图上迭代运行最大和置信传播算法可以得到聚类变量的最大后验概率近似解。
[0069]
子步骤二：消息传播
[0070]
根据最大和置信传播算法的消息传递公式，变量节点和因子节点之间传递五种类型的消息，即
[0071][0072][0073][0074][0075][0076]
通过应用一些数学技巧，我们可以得到如下简化消息：
[0077][0078][0079]
[0080]
消息和都初始化为0，通过消息传播公式不断迭代更新，直到相邻两次迭代消息变化小于终止阈值δ
l
或达到最大迭代次数ι
max
。
[0081]
步骤三：聚类团提取
[0082]
可以根据聚类变量的最大后验概率估计提取量测的传播路径以及同一个目标的量测聚类团。第一步，为了估计量测yi的传播路径，将所有传向量测yi的消息相加，并取使之最大化的路径。第二步，在正确路径相关量测i＝1,...,m，中提取同一个目标的量测聚类团。为了估计聚类变量的值，将所有传入聚类变量的消息相加，并取使之最大化的值。通过聚类变量的估计值提取出聚类中心及其子集。将聚类团内多路径量测通过简单凸组合算法得到目标的量测融合状态估计。
[0083]
本发明解决天波雷达多路径量测聚类问题。该算法通过构建多路径量测聚类的概率图模型，将聚类问题转化为概率图模型隐变量的推断问题，采用最大和置信传播算法近似求解聚类矩阵的最大后验概率。算法优点包括可以自动识别聚类数目，单次消息迭代的计算复杂度为量测个数和传播路径个数乘积的平方。
[0084]
本发明中基于相似性传播的天波雷达多路径量测聚类方法的步骤如下：
[0085]
步骤1、初始化：
[0086]
步骤a、多路径量测建模：对于多路径量测系统，量测yi可能通过特定的传播路径由目标反射得到，也有可能来自杂波。通过所有电离层传播路径将量测坐标配准得到路径相关量测
[0087]
步骤b、聚类建模：对于路径相关量测和考虑二值聚类变量若是的聚类中心，则否则令为联合聚类变量；
[0088]
步骤c、先验概率：根据局部传感器的检测概率以及杂波密度，计算聚类变量的初始置信p(c)，并根据多路径聚类约束得到聚类变量约束；
[0089]
步骤d、似然函数：根据路径相关量测对聚类变量c的统计依赖，计算全局似然函数
[0090]
步骤2、采用最大和置信传播算法进行迭代优化：
[0091]
步骤a、因子图构建：根据后验概率的因子分解式构建航迹关联因子图；
[0092]
步骤b、消息初始：根据步骤1计算的先验概率p(c)和对因子图中的证据因子节点进行消息初始；
[0093]
步骤c、消息迭代：在聚类变量和约束1因子节点和约束2因子节点之间进行消息迭代传播，当消息收敛或达到最大迭代次数时终止迭代；
[0094]
步骤d、聚类团提取：根据聚类变量的最大后验概率估计提取量测的传播路径以及同一个目标的量测聚类团。
[0095]
步骤3：航迹融合：将聚类团内多路径量测通过简单凸组合算法得到目标的量测融合状态估计。
[0096]
现结合实施例、附图对本发明作进一步描述。本研究各部分具体实施细节如下：
[0097]
1.系统建模与问题描述
[0098]
1.1多路径量测模型
[0099]
考虑天波雷达目标探测场景，在某一时刻，k＝1,...,t，表示责任区内第k个目标的状态，t为目标个数。目标状态空间采用麦卡托坐标系。i＝1,...,m，表示第i个量测，m为量测个数。对于多路径量测系统，量测yi可能通过特定的传播路径由目标反射得到，也有可能来自杂波。多路径量测模型可以表示为
[0100][0101]
其中，h
τ
(
·
)，τ＝1,2,...,p，为传播路径τ的量测方程，p为已知的传播路径数，w
τ
为量测噪声，服从协方差为r
τ
的零均值高斯分布，因此，量测也服从高斯分布。假设电离层状态已知，对于传播路径τ，雷达极坐标系到麦卡托坐标系的坐标配准可以通过非线性变换实现t
τ
:为了同时满足精度和计算要求，采用无迹变换对t
τ
进行线性化。因此，坐标配准到麦卡托坐标系的路径相关量测的状态估计近似服从高斯分布。量测集合定义为量测yi的路径相关量测集合定义为全体路径相关量测集合定义为假设目标检测概率为且对不同的传播路径相互独立。假设杂波在雷达探测区域内服从均匀分布，每次扫描的杂波个数服从泊松分布，其中雷达探测区域体积记为v，杂波密度为f
fa
，则泊松分布的均值为
[0102]
1.2多路径量测聚类问题
[0103]
天波雷达多路径量测聚类需要同时识别量测的电离层传播路径，并将来自同一目标的不同路径量测转化到目标坐标系进行聚类。下面结合一个具体实例，阐述该问题：假设雷达探测场景存在2个目标，电离层存在4种传播路径。图1-(a)左侧部分展示了天波雷达多路径量测模型，即单个目标可以通过多种电离层传播路径在雷达系产生多个量测。图1-(a)右侧部分展示了雷达系量测多假设坐标配准，即由于量测的传播路径未知，将其通过全部传播路径配准到目标坐标系得到路径相关量测。图1-(b)为雷达系多路径量测与杂波图，可以看到目标1被传播路径1、2、4检测到，未被传播路径3检测到，目标2被全部传播路径检测到，另外还检测到6个杂波。图1-(c)为目标坐标系下的路径相关量测图，红色虚线框表示目标所有正确路径相关量测，即算法的正确结果，黑色虚线框表示干扰聚类团。
[0104]
1.3多路径量测聚类模型
[0105]
对于路径相关量测和考虑二值聚类变量若是的聚类中心，则否则令为联合聚类变量。为了方便描述，定义多路径聚类变量集合根据多路径量测模型及聚类模型，假设合理的聚类矩阵满足下面约束：
[0106]
约束1：一个量测只能通过一种电离层传播路径得到，一个目标通过一种电离层传播路径至多产生一个量测，表示为i,m＝1,...,m，t,τ＝1,...,p；
[0107]
约束2：一个量测只能选择一个聚类中心，表示为i＝1,...,m；
[0108]
约束3：如果存在量测选择某个量测作为聚类中心，则该量测必须选择自己作为聚类中心，表示为m＝1,...,m。
[0109]
根据上述3个约束，分别得到如下约束方程组：
[0110][0111][0112][0113]
其中，&和|分别表示逻辑与和逻辑或。满足约束1-3的聚类矩阵集合记为
[0114]
路径相关量测对聚类变量c的统计依赖可以表示为全局似然函数根据一般假设，对不同的量测和路径独立，全局似然函数可以因子分解为
[0115][0116]
其中，表示除源自量测yi和传播路径τ以外的路径相关量测，因为根据约束1，这些路径相关量测不可能与来自同一个目标。假设对不同的量测和路径独立，可以进一步因子分解为
[0117][0118]
为量测因此是确定的，上式可以写为
[0119][0120]
其中，定义为
[0121]
[0122]
因此，可得
[0123][0124]
1.4聚类变量的先验概率
[0125]
假设给定目标，聚类变量对不同量测和不同路径是独立的，聚类变量联合先验概率密度函数可以分解为
[0126][0127]
其中，可以定义为
[0128][0129]
因此，可得
[0130][0131]
其中，z为归一化常数，定义为
[0132][0133]
1.5问题描述
[0134]
多路径量测聚类的目的是根据路径相关量测集合估计聚类变量c，在贝叶斯框架下该问题可以转化为最大后验概率估计问题，即
[0135][0136]
通过精确搜索最优聚类矩阵使目标函数最大化是np难问题，需要采用近似算法求解。
[0137]
2相似性量测聚类
[0138]
2.1联合后验概率与因子图
[0139]
使用贝叶斯公式可以得到以下公式
[0140][0141]
其中，为证据因子，定义如下
[0142][0143]
上式表示的联合后验概率密度函数可以表示为如图2所示的因子图，其中
需要注意的是，如消息传播图3所示，每一个变量节点都与对应的证据因子节点连接，为了图示简洁，图2的证据因子节点未画出。
[0144]
根据最大和置信传播算法的消息传递公式，变量节点和因子节点之间传递五种类型的消息，如图3所示，其中即
[0145][0146][0147][0148][0149][0150]
通过应用一些数学技巧，我们可以得到如下简化消息：
[0151][0152][0153][0154]
其中，ψ
m,t
的定义为
[0155][0156][0157][0158][0159]
β和η被省略了。消息和都初始化为0，通过上式不断迭代更新，直到相邻两次迭代消息变化小于终止阈值δ
l
或达到最大迭代次数ι
max
。
[0160]
2.4聚类团提取
[0161]
可以根据聚类变量的最大后验概率估计提取量测的传播路径以及同一个目标的量测聚类团。第一步，提取量测的传播路径。为了估计量测yi的传播路径，将所有传向量测yi的消息相加，并取使之最大化的路径，即
[0162][0163]
量测yi的传播路径为第二步，在正确路径相关量测i＝1,...,m，中提取同一个目标的量测聚类团。为了估计聚类变量的值，将所有传入聚类变量的消息相加，并取使之最大化的值，即
[0164][0165]
然后通过聚类变量的估计值提取出聚类中心及其子集。最后将聚类团内多路径量测通过简单凸组合算法得到目标的量测融合状态估计。
[0166]
3仿真实验
[0167]
3.1仿真场景
[0168]
考虑天波雷达多目标探测场景，如图4所示。仿真场景参数设置如下：电离层为两层结构，即e层和f层，高度分别为he＝100km、hf＝260km，传播路径数为p＝4；目标个数t＝8；雷达接收机经纬度分别为143.20
°
和24.29
°
，瞄准角为325
°
，距离发射机100km；径向距、方位角和径向速度的范围分别为[1100,1800]km，[0.1,0.7]rad和[-0.3,0.3]km/s；量测噪声协方差为r
τ
＝diag(25km2,10-6
km2/s2,9
×
10-6
rad2)。消息传递迭代终止阈值为δ
l
＝10-6
，最大迭代次数为ι
max
＝1000。
[0169]
天波雷达量测聚类是为了得到航迹起始后目标统一的态势，因此通过目标检测与估计测试聚类算法的性能。性能评估之前需要找到检测目标和真实目标一一对应的匹配关系，x轴和y轴的匹配阈值都为10km。评估指标包括：(1)检测正确率(
↑
)：可以与真实目标匹配的检测目标数与检测目标总数之比；(2)漏检率(
↓
)：未与检测目标匹配的真实目标数与真实目标总数之比；(3)均方根误差(
↓
)：可以匹配上的检测目标与真实目标位置的均方根误差；(4)最优次子模式分配(optimal sub pattern assignment,ospa
↓
)[28]；(5)运行时间(
↓
)：算法运行时间。
↑
表示指标的值越大性能越好，
↓
反之。检测正确率越高、漏检率越低表示算法目标检测能力越强，均方根误差越低表示定位精度越高，ospa越低表示综合性能越优，算法运行时间越短，表示算法计算效率越高。蒙特卡洛仿真次数为4000次。蒙特卡洛仿真时，对指标在蒙特卡洛次数上求平均。对比算法为多路径多假设聚类算法，最佳假设树的保留个数分别取为k＝200,400,600,。
[0170]
3.2仿真结果与分析
[0171]
图5展示了pd＝0.75，pf＝2
×
10-5
时，8个目标的多路径量测以及杂波位置分布。图6-(a)和图6-(b)分别表示了相似性传播算法和多假设算法在k＝800时的聚类结果，可以看出相似性传播算法成功检测出全部目标，且无错误检测，多假设算法漏检了目标t3、t5、t7，且产生了两个错误检测。
[0172]
图7(a)-图7-(e)展示了算法性能指标随目标个数的变化。随着目标个数的增加，相似性传播算法的检测正确率呈现上升趋势，这是由于杂波密度确定时误检测数基本确定，正确检测数会随着目标个数的增加而增加，因此检测正确率会上升；相似性传播算法的
漏检率基本保持稳定，这是因为漏检率主要由目标检测概率和杂波密度决定；与检测正确率类似，相似性传播算法的ospa指标略有下降，这是由于误检测数不变而正确检测数增加。而随着目标个数的增加，多假设算法的检测正确率呈现降低趋势，漏检率和ospa呈现上升趋势，这是由于随着目标个数的增加，多假设算法需要枚举的假设个数呈指数式增加，且远远大于算法中设定的假设树保留个数，由于假设树保留的个数不变，多假设算法对假设树空间的搜索能力降低，导致检测正确率、漏检率和ospa性能指标下降。两种算法的均方根误差指标随着目标个数的增加呈现上升趋势，这是由于随着目标个数的增加，杂波被当作目标通过一种传播路径产生的量测与其他真实多路径量测聚类，导致定位精度的降低。图7-(e)为两种算法运行时间随目标个数变化图，实线的坐标轴在左侧，虚线的坐标轴在右侧，可以看出，两种算法运行时间均随着目标个数的增加而增加，这是由于随着目标个数的增加，导致了更大的计算消耗。对于多假设算法，随着假设树的保留数目增加，检测正确率、漏检率、均方根误差和ospa指标性能会提升，而运行时间明显增加，这是由于多假设聚类算法假设树的保留数目越多，可以枚举的聚类假设越多，算法达到较优解的可能越大，同时会导致计算量大幅增加。对于相似性传播和多假设两种聚类算法，相似性传播聚类算法的检测正确率、漏检率、均方根误差、ospa和运行时间指标性能明显优于多假设算法，这是由于相似性传播算法通过消息传递得到的近似最大后验概率可以达到邻域最优，同时计算复杂度较低；而多假设算法通过枚举搜索最优的聚类解，理论上算法性能随着假设树的保留数目增加会不断提升，直到达到最优解，但其时间和空间复杂度都会指数式增长，实际上无法满足实时性和存储要求。