（十）联合概率数据互联原理及应用（JPDA）

前言

一、JPDA原理及算法步骤

（一）算法步骤

1.确认矩阵计算

2.确认矩阵拆分

3.互联概率计算

4.状态及协方差更新

二、仿真验证

（一）模型构建

（二）仿真结果

总结

引用文献

前言

本文主要针对数据关联算法中的联合概率数据互联算法（JPDA）进行了分析与介绍，介绍了其相应的原理与算法步骤，并且针对两个杂波下的匀速运动目标进行了仿真实验，验证了算法的准确性和有效性。后续会逐步介绍其他数据关联算法，读者可以根据自身需求对代码进行修改。

一、JPDA原理及算法步骤

数据互联算法主要用于解决航迹起始之后点迹与航迹的正确互联问题，也就是要解决量测的不确定性问题，为各航迹下一时刻滤波提供正确的量测值。量测的不确定性主要来源于虚假航迹、目标虚警以及多目标点迹判断等原因。点击与航迹的数据互联算法主要包括两种，即极大似然类数据互联算法和贝叶斯类数据互联算法，我们主要研究贝叶斯类数据互联算法，包括以下几种算法：

最近邻域算法（NNSF）

概率数据互联算法（PDA）

综合概率数据互联算法（IPDA）

联合概率数据互联算法（JPDA）

全邻模糊聚类数据互联算法（ANFC）

……

本文主要介绍联合概率数据互联（JPDA）算法，该算法是在联合概率数据互联（PDA）的基础上发展而来的，但是PDA主要用于解决杂波条件下的单雷达单目标问题，JPDA主要用于适用于杂波条件下的多目标数据互联。杂波条件下的多目标数据互联是多目标跟踪中重要但是难以处理的问题，当回波落入多个目标的跟踪波门中时，数据关联过程就要考虑到各个目标航迹对于量测值的竞争问题。

（一）算法步骤

1.确认矩阵计算

确认矩阵定义如下：

$\Omega =[w_{jt}]=\begin{bmatrix} w_{10}&\cdots &w_{1T} \\ \vdots & \cdots & \vdots \\ w_{m_{k}0} & \cdots & w_{m_{k}T} \end{bmatrix}$

其中， $w_{jt}=1$ 表示量测j(j = 1…… $m_{k}$ )落入目标t(t = 0……T)的跟踪波门中， $w_{jt}=0$ 表示量测j未落入目标t的跟踪波门中，当t=0时表示量测为杂波或者虚警，不存在目标，且 $w_{j0}=1$ 。如下图所示：

对应部分代码程序为：

m1 = 0;
[~,n_z_total] = size(Z);
Q1 = zeros(10000,n_target+1); %每个时刻所有的量测点的数量不超过10000，预先准备矩阵
index_cal = [];
z_mid = [];
for k = 1:n_z_total
    flag = 0;
    for q = 1:n_target
        d= Z(:,k) - Z_mid_pre(:,q);
        D(k,q) = d'*inv(S_mid_pre{q})*d;
        if  D(k,q)<= g_sigma
            flag = 1;
            Q1(m1+1,1)=1; %杂波虚警那一列,只有回波落入门内，才有可能是虚警
            Q1(m1+1,q+1)=1; %如果落入回波内，就进行对应的滤波
        end
    end
    if flag == 1
        m1 = m1+1;
        z_mid = [z_mid,Z(:,k)];
        index_cal = [index_cal,k];
    end
end
Q2 = Q1(1:m1,1:n_target+1);

2.确认矩阵拆分

确认矩阵拆分的原则有两个：

基本假设	拆分原则
每个量测有唯一的源，不考虑不可分辨的情况	确认矩阵每行中仅选出一个1作为互联矩阵在该行的唯一非零元素
对于一个给定目标最多一个量测以其为源	在互联矩阵中，除第一列外，其余各列最多只能有一个非零元素

则上述矩阵可以拆分为以下：

对应部分代码为：

%计算各个库
for i = 1 : n1
    k = 0;
   for j = 1:n2
       d = zeros(1,n2);      
       if Q2(i,j) == 1 
           k = k+1;
          d(1,j) = 1;
          c_lib{i}{k} = d;
       end
   end
end
%%
A = zeros(n1,n2,100000);
%拼接各个库并进行判断
num = 0;
for i  = 1:n_total
    c_judge = [];
    for j =  1 : n1
     c_judge = [c_judge;c_lib{j}{index(i,j)}];
    end
    f = sum(c_judge(:,2:end));
    if max(f)<= 1
       num = num+1;
       A(:,:,num) = c_judge;
    end   
end
A1 = A(:,:,1:num);

3.互联概率计算

联合概率数据互联就是计算每个量测与多个目标互联的概率，假设 $\theta _{jt}\left ( k \right )$ 表示量测j源自目标t的事件，其中 $0\leq t\leq T$ ， $\theta _{j0}\left ( k \right )$ 表示量测来源于虚警或者杂波。

按照单目标概率数据互联中条件概率的定义为： $\beta _{jt}\left ( k \right )=P_{r}\left ( \theta _{jt}\left ( k \right )\mid Z^{k} \right )$

则其满足 $\sum_{j=0}^{m_{k}}\beta _{jt}\left ( k \right )=1$

假设确认矩阵按照拆分原则最终拆分为 $n_{k}$ 个矩阵，则表示有 $n_{k}$ 个联合事件，则第j个量测与目标互联的概率为：

$\beta _{jt}\left ( k \right )=P_{r}\left ( \theta _{jt}\left ( k \right )\mid Z^{k} \right )=P_{r}\left ( \bigcup_{i=1}^{n_{k}}\theta _{jt}^{i}\left ( k \right ) \mid Z^{k}\right )=\sum_{i=1}^{m_{k}}\omega _{jt}^{i}\left [ \theta _{k} \right ]P_{r}\left ( \theta _{i}\left ( k \right ) \mid Z^{k}\right )$

其中， $\theta _{jt}^{i}\left ( k \right )$ 表示量测j在第i个联合事件中源于目标t的事件， $\theta _{i}\left ( k \right )$ 表示第i个联合事件。

$\omega _{jt}^{i}\left ( \theta _{i}\left ( k \right ) \right )=\left\{\begin{matrix} 1,\theta _{jt}^{i}\left ( k \right )\subset \theta _{i}\left ( k \right )\\0, other \end{matrix}\right.$

公式表示在第i个联合事件中，如果量测j源于目标q，则为1，否则为0。

$\theta _{i}\left ( k \right )=\bigcap_{j=1}^{m_{k}}\theta _{jt}^{i}\left ( k \right )$

公式表示计算联合事件i的发生概率。

该部分的核心思想为：先计算每个联合事件发生的概率（即将联合事件中为1的时间概率相乘，取交集）；之后进行归一化（即使得所有联合事件概率相加为1）；之后计算量测j源于目标t的概率（即将各个联合事件中涉及到量测j源于目标q的概率相加，取并集）。

量测互联指示：

$\tau _{j}\left [ \theta _{i}\left ( k \right ) \right ]=\sum_{i=1}^{T}\omega _{jt}^{i}\left ( \theta _{i}\left ( k \right ) \right )=\left\{\begin{matrix} 1\\ 0 \end{matrix}\right.$

公式表示量测j在 $\theta _{i}\left ( k \right )$ 中是否与一个真实目标相联。

目标检测指示：

$\sigma _{j}\left [ \theta _{i}\left ( k \right ) \right ]=\sum_{i=1}^{m_{k}}\omega _{jt}^{i}\left ( \theta _{i}\left ( k \right ) \right )=\left\{\begin{matrix} 1\\ 0 \end{matrix}\right.$

公式表示任意量测在联合事件 $\theta _{i}\left ( k \right )$ 中是否与目标t互联，亦是否被检测到。

假量测数量：

$\phi _{j}\left [ \theta _{i}\left ( k \right ) \right ]=\sum_{j=1}^{m_{k}}\left ( 1- \tau _{j}\left [ \theta _{i}\left ( k \right ) \right ]\right )$

$\phi _{j}\left [ \theta _{i}\left ( k \right ) \right ]$ 表示在联合事件 $\theta _{i}\left ( k \right )$ 中假量测的数量。

互联概率计算：

JPDA联合概率数据计算具有两种形式，根据假量测数量 $\phi _{j}\left [ \theta _{i}\left ( k \right ) \right ]$ 的概率质量函数不同而不同（即杂波的分布方式不同而不同）[1]：

        泊松分布：

$P_{r}\left ( \theta _{i}\left ( k \right ) \mid Z^{k}\right )=\frac{\lambda^{\phi _{j}\left [ \theta _{i}\left ( k \right ) \right ]} }{c}\prod_{j=1}^{m_{k}}N_{tj}\left [ z_{j}\left ( k \right ) \right ]^{\tau _{j}\left [ \theta _{i}\left ( k \right ) \right ]}\prod_{t=1}^{T}\left ( P_{D}^{t}\right )^{\sigma _{j}\left [ \theta _{i}\left ( k \right ) \right ]}\left ( 1-P_{D}^{t}\right )^{1-\sigma _{j}\left [ \theta _{i}\left ( k \right ) \right ]}$

        其中， $\lambda$ 为泊松分布参数，c为归一化常数

        均匀分布：

$P_{r}\left ( \theta _{i}\left ( k \right ) \mid Z^{k}\right )=\frac{\phi _{j}\left [ \theta _{i}\left ( k \right ) \right ]\textsc{!}}{V^{\phi _{j}\left [ \theta _{i}\left ( k \right ) \right ] }c}\prod_{j=1}^{m_{k}}N_{tj}\left [ z_{j}\left ( k \right ) \right ]^{\tau _{j}\left [ \theta _{i}\left ( k \right ) \right ]}\prod_{t=1}^{T}\left ( P_{D}^{t}\right )^{\sigma _{j}\left [ \theta _{i}\left ( k \right ) \right ]}\left ( 1-P_{D}^{t}\right )^{1-\sigma _{j}\left [ \theta _{i}\left ( k \right ) \right ]}$

        其中，c为归一化常数，V为确认区域的体积

其中， $n_{z}$ 为量测向量的维数， $N[z_{j}(k)]$ 的具体形式为:

$f(z_{j}(k))=\frac{1}{\left ( 2\pi \right )^{n_{z}}\left | S\left (k \right ) \right |^{1/2}}exp\left \{ -\frac{1}{2} \left ( z_{j}(k)-z_{j}(k\mid k-1) \right ){}'S^{-1}\left ( k \right )\left ( z_{j}(k)-z_{j}(k\mid k-1) \right )\right \}$

对应部分代码为：

% target_indicator =  %目标指示器
Pr = zeros(1,num);
for i1 = 1:num %对每个关联矩阵的概率进行计算
    False_num = m1;
    N=1;
    for j = 1:m1
        mea_indicator =  sum(A_matrix(j,2:end,i1));%观测指示器
        if mea_indicator ==1
            False_num = False_num-1;
            m0 = find(A_matrix(j,2:end,i1));
            b=(z_mid(:,j)-Z_mid_pre(:,m0))'*inv(S_mid_pre{m0})*(z_mid(:,j)-Z_mid_pre(:,m0));
            N=N/sqrt(det((2*pi)^(nz/2)*S_mid_pre{m0}))*exp(-1/2*b);
        end
    end
    if Pd==1
        a=1;
    else
        a=1;
        for j=1:c
            target_indicator=sum(A_matrix(:,j+1,i));                               %参考文献中式4-49
            a=a*Pd^target_indicator*(1-Pd)^(1-target_indicator);                   %计算检测概率
        end
    end
    %       V = sum(ellipse_Volume)*clutter_density;                                         %表示整个空域的体积

    Pr(i1)=N*a*(0.00000000000000000000000001)^False_num;
end
Pr=Pr/sum(Pr);

4.状态及协方差更新

在得到最终互联概率矩阵之后，利用卡尔曼滤波对目标的状态进行更新，公式如下：

状态向量更新：

$x^{t}\left ( k\mid k \right )=x^{t}\left ( k\mid k-1 \right )+K^{t}\left ( k \right ) V^{t}\left ( k \right )$

其中， $V^{t}(k)=\sum_{j=0}^{m_{k}}\beta _{jt}\left ( k \right )V_{j}^{t}(k)$ 为k时刻对目标t用j个量测组合而成的新息向量。

协方差更新：

$P^{t}\left ( k\mid k \right )=P^{t}\left ( k\mid k-1 \right )-(1-\beta _{0t})K^{t}\left ( k \right )S^{t}\left ( k \right )\left [ K^{t}\left ( k \right ) \right ]^{\top }+K^{t}\left ( k \right )\left [ \sum_{j=0}^{m_{k}}\beta _{jt}(k)V_{j}^{t}\left ( k \right )\left (V_{j}^{t}\left ( k \right )\right )^{\top } -V^{t}\left ( k \right )\left (V^{t}\left ( k \right )\right )^{\top } \right ]\left ( K^{t}\left ( k \right ) \right )^{\top }$

二、仿真验证

（一）模型构建

建立两个匀速目标的CV模型，模型具体构建方法见之前的博客：

（四）交互式多模型（IMM）中二维平面CV/CA/CT模型构建https://blog.csdn.net/h320130/article/details/135141783 设置杂波密度，跟踪门限等参数，具体代码如下：


sim_time=150;        %仿真时间
T=1;                %T为采样间隔
n_sim = sim_time/T; %仿真次数
MC_number=5;        %Monte Carlo仿真次数
n_target=2;         %目标个数
initial_v=[1500 30 500 40; 
           500 40 1500 30];%初始目标状态
nf = 0.1; %过程噪声强度
P0 = 1*eye(4);     %各目标的初始跟踪误差
clutter_density = 0.000001; %杂波密度
R = 1*eye(2); %目标量测误差
%生成目标的运动
X_location = cell(1,n_target);
X_location(:) = {zeros(4,n_sim)};
for q = 1:n_target
    [X_location{q},F,Q] = track(initial_v(q,:)',n_sim,T,nf,P0);
end