CN114895296B - 一种基于fda-mimo雷达的多目标跟踪方法 - Google Patents

Abstract

本发明涉及一种基于FDA‑MIMO雷达的多目标跟踪方法，包括：S1、建立目标状态演化模型和多目标叠加量测模型；S2、将上一时刻的多目标后验概率密度由伯努利参数集表征；S3、预测当前时刻的存活伯努利成分，获取当前时刻的新生标签伯努利成分；S4、更新存活伯努利成分中的粒子权重和新生标签伯努利成分中的粒子权重；S5、对更新后的伯努利成分进行融合；S6、对目标个数进行估计，提取目标状态，并返回步骤S3估计下一时刻的目标状态。该方法能够同时估计多个目标的个数和状态，并输出每个目标的航迹，同时根据雷达接收信号的特性，目标跟踪过程中并不考虑检测概率，解决了检测概率未知的多目标跟踪问题。

Description

一种基于FDA-MIMO雷达的多目标跟踪方法

技术领域

本发明属于信号处理技术领域，具体涉及一种基于FDA-MIMO雷达的多目标跟踪方法。

背景技术

多目标跟踪技术(Multi-target tracking,MTT)在军事和民用领域都是一项不可或缺的重要技术。随着技术的快速发展，MTT技术也被广泛的应用于医学、视频跟踪和雷达等各个领域。早期MTT的主要方法为基于数据关联的跟踪方法，但是在虚警目标密集且目标个数未知的情况下，基于数据关联的多目标跟踪方法容易产生“组合爆炸”问题。基于随机有限集(Random finite set,RFS)的多目标跟踪算法则能够有效的避免复杂的数据关联问题。根据近几年学者们的研究,可将RFS框架下的多目标跟踪技术大致分为两种：第一种是通过传递多目标后验概率密度的一阶矩来实现对多目标状态的估计；第二种是将目标状态刻画为多伯努利(Bernoulli)集，通过对多Bernoulli集的后验概率密度进行传递从而实现多目标的状态跟踪。因为第二类滤波器不需要通过聚类的方式提取目标状态，因此它们的跟踪性能不会受到聚类算法的影响，尤其是标签多目标多伯努利滤波，其不仅能够对目标状态进行估计，还可以输出各个目标的航迹。

不同于相控阵雷达(PA)，频率分集阵列(FDA)雷达在各个发射单元之间引入了一个频率差，使得FDA的发射导向矢量不仅包含了目标的角度信息，而且包含了目标的距离信息，因此FDA雷达的发射方向图不仅与角度有关，而且与距离由关。FDA雷达的这一特性给距离-角度联合估计提供了基础依据。为了更好的实现目标的距离和角度的联合估计，FDA常与多输入多输出(Multiple input multiple output,MIMO)系统相互结合。FDA-MIMO雷达作为一种新体制雷达，目前关于它的研究还多集中于参数估计、杂波抑制和抗干扰等方面，而对于FDA-MIMO体制下的单目标跟踪问题研究较少，对于复杂环境下FDA-MIMO雷达多目标跟踪的研究则处于空白阶段。

发明内容

为了解决现有技术中存在的上述问题，本发明提供了一种基于FDA-MIMO雷达的多目标跟踪方法。本发明要解决的技术问题通过以下技术方案实现：

本发明实施例提供了一种基于FDA-MIMO雷达的多目标跟踪方法，包括步骤：

S1、建立目标状态演化模型和多目标叠加量测模型；

S2、将上一时刻的多目标后验概率密度由伯努利参数集表征；

S3、利用所述目标状态演化模型和所述伯努利参数集预测当前时刻的存活伯努利成分，同时获取所述当前时刻的新生标签伯努利成分；

S4、结合所述多目标叠加量测模型，利用Music谱和似然函数更新所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重，得到当前时刻更新后的伯努利成分；

S5、基于欧式距离，对所述更新后的伯努利成分进行融合，得到当前时刻融合后的伯努利成分和融合后的多目标后验概率密度；

S6、根据所述融合后的多目标后验概率密度对目标个数进行估计，根据所述融合后的伯努利成分提取目标状态，得到当前时刻的目标状态估计结果，并返回步骤S3，结合所述融合后的伯努利成分估计下一时刻的目标状态。

在本发明的一个实施例中，所述多目标叠加量测模型为：

Y_k＝[y₁₁,…,y_1N,…,y_N1,…y_NM]^T+N

其中，N为白高斯噪声，y_mn为由第n个阵元发射，经过目标反射后，由第m个阵元接收的经过匹配滤波后的信号；

其中，为第i个目标的多普勒频率，N_k为k时刻的目标个数，ξ_i为第i个目标复散射系数，j为虚数单位，f_n为第n个阵元发射信号的频率，r_i为第i个目标距离，t为时间指标，v_i为第i个目标的速度，λ₀为波长，为由阵元间距引起的发射时延，为由阵元间距引起的接收时延，Δf为频率步进量，n为第n个阵元，f₀为参考频率。

在本发明的一个实施例中，所述伯努利参数集为：

其中，为k-1时刻第i个Bernoulli成分的存在概率，N_k表示k-1时刻的Bernoulli成分总数，k为时间指标，为k-1时刻的空间概率密度；

其中，表示粒子数，为k-1时刻第i个目标状态，为狄拉克脉冲函数，s表示第s个粒子，为刻画第i个目标的第s个粒子，为对应粒子的权重，为一簇加权粒子。

在本发明的一个实施例中，所述存活伯努利成分包括第一存在概率、第一空间分布概率密度和第一粒子权重，其中，

所述第一存在概率为：

其中，为k-1时刻第i个Bernoulli成分的存在概率，为预测后第i个Bernoulli成分中的粒子个数，lk|k-1,i＝lk-1,i表示预测后对应Bernoulli成分所携带的标签不变，P_S,k表示粒子存活并转移的概率，表示k-1时刻刻画第i个目标的第s个粒子，为对应粒子的权重；

所述第一空间分布概率密度为：

其中，为预测后第i个Bernoulli成分中的粒子数，为归一化的粒子权重，为狄拉克脉冲函数，s指示第s个粒子，S表示存活；

所述第一粒子权重为：

其中，为预测后第i个Bernoulli成分中存活的第s个粒子的权重，为粒子数，为粒子权重，为目标状态的转移核，为预测后的粒子状态，为k-1时刻的粒子状态，P_S,k为粒子存活并转移的概率，S表示存活。

在本发明的一个实施例中，所述新生标签伯努利成分包括第二存在概率、第二空间分布概率密度、第二粒子权重和目标标签，其中，

所述第二存在概率为先验给定的第j个新生伯努利成分的存在概率；

所述第二空间分布概率密度为：

其中，为k时刻新生的第j个Bernoulli成分的空间分布概率密度，为初始化粒子集，为第j个新生Bernoulli成分中第s个粒子的权重，l_k|k-1,j为新生Bernoulli成分的标签，为新生第j个Bernoulli成分中的粒子数，Γ表示新生；

所述第二粒子权重为：

其中，为第j个新生Bernoulli成分中第s个粒子的权重，为新生第j个Bernoulli成分中的粒子数；

所述目标标签为：

l_k|k-1,j＝(k,j)

其中，l_k|k-1,j为新生伯努利成分的标签，k表示时刻，j为出现顺序。

在本发明的一个实施例中，步骤S4包括：

S41、结合所述多目标叠加量测模型，利用当前的目标状态预测值和量测值计算Music谱：

其中，克罗内克乘积kron(b(·),a(·))为联合发射-接收域的权矢量，a(·)和b(·)分别表示发射导向矢量和接收导向矢量，P_Music为Music谱，Y_k为多目标叠加量测值，(·)^H表示共轭转置，U_n,k为k时刻的噪声子空间，为将目标状态映射的角度，为将目标状态映射的径向距离；

S42、利用所述Music谱替换传统似然函数，得到新的似然函数：

g(Y_k|·)＝{P_Music(Y_k|·)}^μ

其中，g(Y_k|·)为新的似然函数，P_Music(Y_k|·)为Music谱，μ为u次方；

S43、保持粒子状态不变，利用所述新的似然函数更新所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重：

其中，为k时刻预测的第i个Bernoulli成分中第s个粒子的权重，对应更新后粒子的权重，为新的似然函数；

S44、根据所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重得到所述更新后的伯努利成分：

其中，为更新后的存在概率，为预测的存在概率，为更新后粒子数，更新前后标签保持不变l_k,i＝l_k|k-1,i，l_k,i为更新后第i个Bernoulli成分的标签，为k时刻更新后粒子的权重，为预测的空间分布概率密度，g(Y_k|x)为似然函数，为根据量测Y_k更新后的空间分布概率密度，为狄拉克脉冲函数，s表示第s个粒子；

S45、根据所述更新后的伯努利成分得到更新后的多目标后验概率密度函数：

其中，π_k为k时刻多目标后验概率密度，为更新后Bernoulli成分的存在概率，为更新后Bernoulli成分的空间分布概率密度，为更新后的Bernoulli成分个数，与预测后的Bernoulli成分个数N_k|k-1保持一致。

在本发明的一个实施例中，步骤S5包括：

S51、根据所述更新后的伯努利成分中的目标状态的空间分布概率密度提取目标状态，得到所述目标状态集；

S52、根据所述目标状态集计算欧式距离矩阵，并根据所述欧式距离矩阵和距离阈值划分目标子群；其中，所述欧式距离矩阵为：

其中，为计算目标状态和之间的欧氏距离；

S53、将同一个所述目标子群内的目标状态对应的伯努利成分进行融合，得到所述融合后的伯努利成分和所述融合后的多目标后验概率密度。

在本发明的一个实施例中，步骤S53包括：

对同一个所述目标子群内的目标状态对应的粒子集进行融合，得到融合后的粒子集：

其中，和分别为刻画第i个和第j个Bernoulli成分的粒子集；

对同一个所述目标子群内的目标状态对应的粒子权重集进行融合，得到融合后的粒子权重集：

其中，和分别为刻画第i个和第j个Bernoulli成分的粒子对应的权重集，和分别为k时刻第i个和第j个Bernoulli成分中的第s个粒子的权重，和为分别为第i个和第j个Bernoulli成分中粒子的个数；

对同一个所述目标子群内的目标状态对应的目标存在概率集进行融合，得到融合后的目标存在概率：

其中，和分别为第i个和第j个Bernoulli成分的存在概率；

对同一个所述目标子群内的目标状态对应的目标标签集进行融合，得到融合后的目标标签：

其中，为融合后的目标标签，和为不同的标签，k₁和k₂为分别表示时间指标k₁和k₂，i和j为出现顺序；

k时刻所述融合后的多目标后验概率密度为：

其中，π_k为融合后的多目标后验概率密度，为融合后的存在概率，为融合后的空间分布概率密度。

在本发明的一个实施例中，根据所述融合后的多目标后验概率密度对目标个数进行估计，根据所述融合后的伯努利成分提取目标状态，得到当前时刻的目标状态估计结果，包括：

S61、根据所述融合后的多目标后验概率密度对目标个数进行估计：

其中，为融合后的前大存在概率的伯努利成分，N_k为k时刻更新后的Bernoulli成分总数，Ω为所有Bernoulli成分的指标集，和分别为更新后第i和第j个Bernoulli成分的存在概率，Ω_e为Ω的子集，Ω_e表示k时刻假设存在的所有目标对应的指标集合；

S62、从所述融合后的伯努利成分中提取前个目标存在概率的成分，得到所述目标状态：

其中，为粒子状态，为粒子个数，为更新后刻画第i个目标的第s个粒子的状态，为对应粒子的权重；

S63、根据所述目标个数和所述目标状态得到所述目标状态估计结果。

与现有技术相比，本发明的有益效果：

1、本发明的多目标跟踪方法首先建立了雷达多目标叠加量测模型，其次利用Music谱替代传统似然函数，结合新的似然函数和随机集框架下的多伯努利滤波实现对假目标存在情况下的多目标跟踪，最后提出一种基于欧氏距离的航迹融合算法，解决了同一目标航迹过估计问题、多目标叠加量测的强非线性问题以及无论假目标存在主瓣区域或者旁瓣区域的假目标存在问题，能够同时估计多个目标的个数和状态，并输出每个目标的航迹，同时根据雷达接收信号的特性，目标跟踪过程中并不考虑检测概率，解决了检测概率未知的多目标跟踪问题。

2、本发明的多目标跟踪方法引入标签技术，不仅可以实现对目标个数和状态的估计，还能够输出各个目标的航迹。

附图说明

图1为本发明实施例提供的一种基于FDA-MIMO雷达的多目标跟踪方法的流程示意图；

图2为本发明实施例提供的一种FDA-MIMO雷达的结构示意图；

图3为本发明实施例提供的一种真实目标轨迹的示意图；

图4为本发明实施例提供的一种不同μ情况下OSPA距离(SNR＝-15dB)的示意图；

图5为本发明实施例提供的一种不同μ情况下平均OSPA距离的示意图；

图6为本发明实施例提供的SNR＝-15dB时IL-MeMBer滤波器的多目标跟踪效果图；

图7为本发明实施例提供的SNR＝-18dB时IL-MeMBer滤波器的多目标跟踪效果图；

图8为本发明实施例提供的SNR＝-15dB时I-MeMBer滤波器的多目标跟踪效果图；

图9为本发明实施例提供的一种不同时刻的Music空间谱和IL-MeMBer滤波器的跟踪结果示意图；

图10本发明实施例提供的一种不同SNR下IL-MeMBer滤波器目标个数估计图；

图11为本发明实施例提供的一种不同SNR下IL-MeMBer滤波器多目标状态OSPA距离估计图；

图12为本发明实施例提供的一种不同算法目标个数估计图；

图13为本发明实施例提供的一种不同算法OSPA距离估计图。

具体实施方式

下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。

实施例一

请参见图1，图1为本发明实施例提供的一种基于FDA-MIMO雷达的多目标跟踪方法的流程示意图。

本实施例在FDA-MIMO雷达系统下，采用改进的带标签的多目标多伯努利(Improved labeled multi-target multi-Beroulli,IL-MeMBer)滤波器进行多目标跟踪，其进行多目标跟踪的思路为：首先，在虚假目标存在的情况下，给出了FDA-MIMO雷达多目标叠加量测模型。第二，利用传统点量测似然函数和Music空间谱的相似特征，利用Music空间谱替代传统点量测似然函数，推导了FDA-MIMO多目标叠加量测模型下的似然函数。第三，针对目标个数时变、目标个数未知、未经恒虚警率(Constant False-Alarm Rate,CFAR)处理的雷达回波信号不存在检测概率以及传统多目标多伯努利(Multi-Target Multi-Bernoulli,MeMBer)滤波器无法输出目标航迹等问题，通过引入标签的策略推导了不考虑检测概率情况下的IL-MeMBer滤波器的解析解。最后，针对传统MeMBer滤波器在目标跟踪过程中可能多估计目标的问题，提出了基于欧式距离的融合策略。除此之外，给出了IL-MeMBer滤波器的序贯蒙特卡洛实现。

具体的，该多目标跟踪方法包括步骤：

S1、建立目标状态演化模型和多目标叠加量测模型。

具体的，目标状态演化模型为：

x_k＝f(x_k-1,v_k-1) (1)

其中，f为状态转移函数(线性或非线性)，k为时间指标，v为过程噪声，[x,y]为目标位置信息，为目标速度信息。

请参见图2，图2为本发明实施例提供的一种FDA-MIMO雷达的结构示意图。图2中，FDA-MIMO雷达包含N个发射阵元和M个接收阵元，其中，N＝M，各个阵元之间的间距d为半波长，发射端相邻阵元发射的信号之间存在一个恒定的频率差Δf，第一个阵元的发射频率为参考频率f₁，第n个发射阵元发射的信号为发射信号经过N_k＝N_k,T+N_k,Γ(N_k,T,N_k,Γ分别表示真、假目标个数)个目标反射，接收端匹配滤波后的信号建模为多目标叠加量测模型：

Y_k＝[y₁₁,…,y_1N,…,y_N1,…y_NM]^T+N (2)

其中，N为白高斯噪声，y_mn为由第n个阵元发射，经过目标反射后，由第m个阵元接收的经过匹配滤波后的信号。，

其中，为第i个目标的多普勒频率，N_k为k时刻的目标个数，ξ_i为第i个目标复散射系数，j为虚数单位，f_n为第n个阵元发射信号的频率，r_i为第i个目标距离，t为时间指标，v_i为第i个目标的速度，λ₀为波长，为由阵元间距引起的发射时延，为由阵元间距引起的接收时延，Δf为频率步进量，n为第n个阵元，f₀为参考频率。

假设在每个时刻进行K次快拍的采样，则k时刻采样得到的数据矩阵为Y_k中包含了k时刻所有目标的量测信息，同时Y_k也是定义的多目标叠加量测模型。

S2、将上一时刻的多目标后验概率密度由伯努利参数集表征。

本实施例中，上一时刻是指k时刻的上一时刻即k-1时刻。

具体的，对感兴趣的目标区域进行目标状态初始化，根据初始化的滤波器参数，假设IL-MeMBer滤波器k-1时刻产生的带标签的伯努利成分为则k-1时刻的多目标后验概率密度可以由伯努利(Bernoulli)参数集表示：

其中，为k-1时刻第i个Bernoulli成分的存在概率，N_k表示k-1时刻的Bernoulli成分总数，k为时间指标，为k-1时刻的空间概率密度。

进一步的，空间概率密度可以由一簇加权粒子表示：

其中，表示粒子数，为k-1时刻第i个目标状态，为狄拉克脉冲函数，s表示第s个粒子，为刻画第i个目标的第s个粒子，为对应粒子的权重，为一簇加权粒子。

S3、利用所述目标状态演化模型和所述伯努利参数集预测当前时刻的存活伯努利成分，同时获取所述当前时刻的新生标签伯努利成分。

在一个具体实施例中，Bernoulli成分的预测分为存活部分的预测和新生部分的预测。

具体的，存活伯努利成分包括第一存在概率、第一空间分布概率密度和第一粒子权重。本实施例中，对各个存活的伯努利成分中的粒子，利用目标状态演化模型进行状态预测，即对第一粒子权重进行预测、对存活伯努利成分的第一存在概率和第一空间分布概率密度进行预测。

其中，第一存在概率为：

其中，为k-1时刻第i个Bernoulli成分的存在概率，为预测后第i个Bernoulli成分中的粒子个数，l_k|k-1,i＝l_k-1,i表示预测后对应Bernoulli成分所携带的标签不变，P_S,k表示粒子存活并转移的概率，表示k-1时刻刻画第i个目标的第s个粒子，为对应粒子的权重。

所述第一空间分布概率密度为：

其中，为预测后第i个Bernoulli成分中的粒子数，为归一化的粒子权重，为狄拉克脉冲函数，s指示第s个粒子，S表示存活。

所述第一粒子权重为：

其中，为预测后第i个Bernoulli成分中存活的第s个粒子的权重，为粒子数，为粒子权重，为目标状态的转移核，为预测后的粒子状态，为k-1时刻的粒子状态，P_S,k为粒子存活并转移的概率，S表示存活。

具体的，根据预设参数新生带标签的Bernoulli成分用以捕获新生目标。新生标签伯努利成分包括第二存在概率、第二空间分布概率密度、第二粒子权重和目标标签。新生的Bernoulli成分可以参数化为其中，

第二存在概率(例如0.02)为先验给定的第j个新生伯努利成分的存在概率。

第二空间分布概率密度可以由粒子集表征：

其中，初始化的粒子集由Gaussian分布采样得到，和分别为先验给定的均值和协方差。

第二粒子权重为粒子的初始化权重，其表示为：

其中，为第j个新生Bernoulli成分中第s个粒子的权重，为新生第j个Bernoulli成分中的粒子数。

新生Bernoulli成分的标签即目标标签由新生的时刻k以及新生的顺序j确定：

l_k|k-1,j＝(k,j) (12)

其中，l_k|k-1,j为新生伯努利成分的标签，k表示时刻，j为出现顺序。

所以，k时刻预测后的多目标后验概率密度可以表示为：

S4、结合所述多目标叠加量测模型，利用Music谱和似然函数更新所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重，得到当前时刻更新后的伯努利成分。

S41、结合所述多目标叠加量测模型，利用当前的目标状态预测值和量测值计算Music谱。

具体的，给定k时刻多目标叠加量测Y_k，利用当前的目标状态预测值和量测值计算Music谱：

其中，克罗内克乘积kron(b(·),a(·))为联合发射-接收域的权矢量，a(·)和b(·)分别表示发射导向矢量和接收导向矢量，P_Music为Music谱，Y_k为多目标叠加量测值，(·)^H表示共轭转置，U_n,k为k时刻的噪声子空间，由R_k的前个小特征值对应的特征向量构成，R_k为k时刻接受信号的协方差矩阵，N为发射阵元数，M为接收阵元数，为通过利用最小描述长度(Minimum description length,MDL)算法估计的到的第k个时刻的目标个数，为将目标状态映射的角度，为将目标状态映射的径向距离。

其中，a(θ_i,r_i)为依赖于目标距离r_i和角度θ_i的发射导向矢量，b(θ_i)表示仅依赖于目标角度的接收导向矢量，为第n个阵元和参考阵元之间的相位差。

S42、利用所述Music谱替换传统似然函数，得到新的似然函数：

g(Y_k|·)＝{P_Music(Y_k|·)}^μ (17)

其中，g(Y_k|·)为新的似然函数，P_Music(Y_k|·)为Music谱，μ为μ次方。

本实施例利用Music谱替代传统似然函数，解决了FDA-MIMO雷达多目标叠加量测的强非线性问题。

S43、保持粒子状态不变，利用所述新的似然函数更新所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重。

具体的，在IL-MeMBer滤波器的序贯蒙特卡洛实现的更新过程中，粒子的状态保持不变粒子的权重根据新的似然函数计算：

其中，为k时刻预测的第i个Bernoulli成分中第s个粒子的权重，对应更新后粒子的权重，为新的似然函数；

S44、根据所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重得到所述更新后的伯努利成分：

其中，为更新后的存在概率，为预测的存在概率，为更新后粒子数，更新前后标签保持不变l_k,i＝l_k|k-1,i，l_k,i为更新后第i个Bernoulli成分的标签，为k时刻更新后粒子的权重，为预测的空间分布概率密度，g(Y_k|x)为似然函数，为根据量测Y_k更新后的空间分布概率密度，为狄拉克脉冲函数，s表示第s个粒子；

S45、根据所述更新后的伯努利成分得到更新后的多目标后验概率密度函数：

其中，π_k为k时刻多目标后验概率密度，为更新后Bernoulli成分的存在概率，为更新后Bernoulli成分的空间分布概率密度，为更新后的Bernoulli成分个数，与预测后的Bernoulli成分个数N_k|k-1保持一致。

S5、基于欧式距离，对所述更新后的伯努利成分进行融合，得到当前时刻融合后的伯努利成分和融合后的多目标后验概率密度。

具体的，根据预先设置的欧式距离阈值，对更新后的伯努利成分进行融合。融合后的粒子集合和权重集合为分别为融合前粒子集的并集和权重集的并集。更新后伯努利成分的标签为最先出现的标签。

步骤S5具体包括步骤：

S51、根据所述更新后的伯努利成分中的目标状态的空间分布概率密度提取目标状态，得到所述目标状态集

S52、根据所述目标状态集计算欧式距离矩阵，并根据所述欧式距离矩阵和距离阈值划分目标子群。

首先，根据提取得到的目标状态集计算欧式距离矩阵

其中，为计算目标状态和之间的欧氏距离。

然后，根据所述欧式距离矩阵和距离阈值划分目标子群。具体的，将k时刻估计到的多个目标看作多个群目标，根据D和距离阈值(例如200m)将所有目标划分为N_k个子群，则任意一个子群内的目标都认为是同一个目标估计得到的状态，划分后的目标子群可以表示为：

其中，p表示第p个子群，N_k表示共划分为N_k个子群。

S53、将同一个所述目标子群内的目标状态对应的伯努利成分进行融合，得到所述融合后的伯努利成分和所述融合后的多目标后验概率密度。

具体的，同一个目标子群内的目标状态对应的伯努利成分包括粒子集、粒子权重集、目标存在概率和目标标签，则对同一个目标子群内任意两个伯努利成分进行融合包括：

对同一个所述目标子群内的目标状态对应的粒子集进行融合，得到融合后的粒子集：

其中，和分别为刻画第i个和第j个Bernoulli成分的粒子集；

对同一个所述目标子群内的目标状态对应的粒子权重集进行融合，得到融合后的粒子权重集：

其中，和分别为刻画第i个和第j个Bernoulli成分的粒子对应的权重集，和分别为k时刻第i个和第j个Bernoulli成分中的第s个粒子的权重，和为分别为第i个和第j个Bernoulli成分中粒子的个数。

对同一个所述目标子群内的目标状态对应的目标存在概率集进行融合，得到融合后的目标存在概率：

其中，和分别为第i个和第j个Bernoulli成分的存在概率。

对同一个所述目标子群内的目标状态对应的目标标签集进行融合，得到融合后的目标标签。具体的，因为目标标签是根据出现的时刻和先后顺序来创造的，因此对于不同的目标标签和的融合可以公式化为：

其中，为融合后的目标标签，和为不同的标签，k₁和k₂为分别表示时间指标k₁和k₂，i和j为出现顺序。

进一步的，k时刻融合后的多目标后验概率密度为：

其中，π_k为融合后的多目标后验概率密度，为融合后的存在概率，为融合后的空间分布概率密度。

S6、根据所述融合后的多目标后验概率密度对目标个数进行估计，根据所述融合后的伯努利成分提取目标状态，得到当前时刻的目标状态估计结果，并返回步骤S1，结合所述融合后的伯努利成分估计下一时刻的目标状态。

具体的，根据所述融合后的多目标后验概率密度对目标个数进行估计，根据所述融合后的伯努利成分提取目标状态，得到当前时刻的目标状态估计结果包括步骤：

S61、根据所述融合后的多目标后验概率密度对目标个数进行估计：

其中，为融合后的前大存在概率的伯努利成分，N_k为k时刻更新后的Bernoulli成分总数，Ω为所有Bernoulli成分的指标集，和分别为更新后第i和第j个Bernoulli成分的存在概率，Ω_e为Ω的子集，Ω_e表示k时刻假设存在的所有目标对应的指标集合。

S62、从所述融合后的伯努利成分中提取前个目标存在概率的成分，得到所述目标状态。具体的，目标状态从融合后的前个大存在概率的伯努利成分中提取，具体可描述为对应粒子集的加权求和。

其中，为粒子状态，为粒子个数，为更新后刻画第i个目标的第s个粒子的状态，为对应粒子的权重；

S63、根据所述目标个数和所述目标状态得到所述目标状态估计结果。

具体的，目标状态估计结果包括步骤S61得到的目标个数和步骤S62得到的目标状态。

进一步的，在对当前目标状态进行估计的同时，将融合后k时刻的伯努利成分带入步骤S1中进行迭代，对k+1时刻的多目标状态进行估计。

本实施例的多目标跟踪方法首先建立了雷达多目标叠加量测模型，其次利用Music谱替代传统似然函数，结合新的似然函数和随机集框架下的多伯努利滤波实现对假目标存在情况下的多目标跟踪，最后提出一种基于欧氏距离的航迹融合算法，解决同一目标航迹过估计问题，解决了多目标叠加量测的强非线性问题，能够有效的解决无论假目标存在主瓣区域或者旁瓣区域的假目标存在问题，能够同时估计多个目标的个数和状态，并输出每个目标的航迹，并且根据雷达接收信号的特性，目标跟踪过程中并不考虑检测概率，解决了检测概率未知的多目标跟踪问题。另外，该方法引入标签技术，不仅可以实现对目标个数和状态的估计，还能够输出各个目标的航迹。

实施例二

在实施例一的基础上，本实施例通过仿真对基于FDA-MIMO雷达的多目标跟踪方法进行说明。

1、仿真参数设置：

在仿真实验中共有9个目标随机出现在观测区域内，所有目标均服从匀速直线运动模型，实验中目标的出现和死亡时刻如图3所示，图3为本发明实施例提供的一种真实目标轨迹的示意图。实验中持续跟踪的时间长度为100S，观测区域大小为在预测过程中目标的存活概率为Ps＝0.98，目标存活并转移的概率为f_k|k-1(x|·)＝1。光速3×10⁸m/s，参考频率1×10¹⁰Hz，频率步进量3×10⁴Hz，发射阵元和接收阵元个数15，采样快拍数2000。匀速直线运动模型：

其中，v_k服从均值为0，协方差为Q的白Gaussian过程噪声，

参数q＝5。观测区域假目标个数建模为参数λ_c＝3的Poisson分布，假目标位置随机分布于观测区域内。新生目标概率为r^(1,2)＝0.02，r^(3,4)＝0.03新生目标分布为服从的高斯分布，其中：

m⁽¹⁾＝[5200m；0m·s^-1；3000m；0m·s^-1]；

m⁽²⁾＝[4000m；0m·s^-1；6000m；0m·s^-1]；

m⁽³⁾＝[0m；0m·s^-1；6900m；0m·s^-1]；

m⁽⁴⁾＝[4000m；0m·s^-1；7450m；0m·s^-1]；

P＝[40m²；40m²·s^-2；40m²；40m²·s^-2]。

在SMC实现过程中，新生Bernoulli成分中的粒子数目为300，存活Bernoulli成分中的粒子数目的最大值伪500。对于算法的跟踪性能，本实施例中采用著名的采用最优子模式分配距离(Optimal Subpattern Assignment，OSPA)来评估。因为传统的RMSE误差仅能看出目标位置估计误差，而OSPA距离则能够联合估计目标的个数和位置误差。

在上述仿真参数下，对本实施例方法与传统方法在FDA-MIMO体制下的多目标跟踪性能进行仿真，并对实验结果进行分析。

2、仿真内容：

仿真1：对式g(Y_k|·)＝{P_Music(Y_k|·)}^μ中不同μ值下的多目标跟踪性能进行分析，并确定不同SNR下的近似最优μ值。

请参见图4，图4为本发明实施例提供的一种不同μ情况下OSPA距离(SNR＝-15dB)的示意图。图4中，“mean”表示所有时刻OSPA的均值。由图4可以看出在SNR＝-15dB时，μ＝6.4的情况下多目标跟踪性能最优。当μ＝6.5和μ＝6.6时，多目标跟踪性能稍微有所下降，但性能下降有限，说明μ值小幅度的变化对IL-MeMBer滤波器的多目标跟踪性能的影响有限。当μ＝4和μ＝6时，IL-MeMBer滤波器的多目标跟踪性能明显下降，尤其是在目标出生(1s,15s,28s,40s,55s,66s)和目标死亡(55s,65s,70s,80s,95s)时刻附近的多目标的状态估计会出现较大的误差。

请参见图5，图5为本发明实施例提供的一种不同μ情况下平均OSPA距离的示意图。图5中统计了不同SNR下100个时刻的平均OSPA距离估计，由图5可以看出一个合适的μ值是IL-MeMBer滤波器所必须的。结合图4和图5可知，如果μ值的选取过大或者过小，都会造成较大的估计误差，图5中每条曲线的谷值被认为是相应SNR下的近似最优的μ值。统计不同SNR下的近似最优的μ值：SNR＝-10dB,μ取7.1；SNR＝-15dB,μ取6.4；SNR＝-18dB,μ取5.0。

仿真2：对IL-MeMBer滤波器的跟踪性能进行实验验证。

请参见图6、图7和图8，图6为本发明实施例提供的SNR＝-15dB时IL-MeMBer滤波器的多目标跟踪效果图，图7为本发明实施例提供的SNR＝-18dB时IL-MeMBer滤波器的多目标跟踪效果图，图8为本发明实施例提供的SNR＝-15dB时I-MeMBer滤波器的多目标跟踪效果图。由图6和图7可以看出SNR＝-15dB和SNR＝-18dB时IL-MeMBer滤波器都能很好的估计多个目标的状态。对比图6和图7可以看出，SNR减小时，目标未知估计的精度会有所下降。图8中I-MeMBer滤波器指的是增加本文中所提的航迹融合步骤后的MeMBer滤波器。由图8可以看出，I-MeMBer滤波器也能够在SNR＝-15dB时实现对多目标状态的估计，但是因为I-MeMBer滤波器没有标签，所以它输出的是目标的点迹，而不是如图6和图7所示的目标的航迹。

请参见图9，图9为本发明实施例提供的一种不同时刻的Music空间谱和IL-MeMBer滤波器的跟踪结果示意图。如图9所示，假目标的个数、位置在不同时刻是具备随机性的，除此之外，真实目标的个数和状态是时变的。在不同时刻的Music空间谱表明无论是在假目标处或者在真实目标位置处，都会形成一个略高于背景强度的谱峰(SNR＝-15dB)；显然单纯的采用Music算法是无法处理假目标的情况的；相反，IL-MeMBer滤波器则能够利用目标状态在不同时刻的关联性，将假目标剔除掉，而且无需考虑假目标是在如图9(a)和图9(f)所示的主瓣或者旁瓣。

请参见图10和图11，图10本发明实施例提供的一种不同SNR下IL-MeMBer滤波器目标个数估计图，图11为本发明实施例提供的一种不同SNR下IL-MeMBer滤波器多目标状态OSPA距离估计图。图11(a)、图11(b)和图11(c)分别为多目标位置估计OSPA距离、估计得到的OSPA距离中的位置分量和估计得到的OSPA距离中的夺目目标的势分量。结合图10和图11(c)可以看出，SNR不同时，IL-MeMBer滤波器都能较为准确的估计得到目标的个数。图11(b)则反映出随着SNR的增大，IL-MeMBer滤波器在多目标位置估计方面会更准确。图11(a)则从全局上反映出IL-MeMBer滤波器的多目标跟踪性能会随着SNR的提高而提升。

仿真3：对比所提IL-MeMBer滤波器、增加了轨迹融合的传统的MeMBer(I-MeMBer)滤波器和不进行轨迹融合的IL-MeMBer滤波器的多目标跟踪性能。仿真3中SNR＝-15dB,μ＝6.4。

请参见图12，图12为本发明实施例提供的一种不同算法目标个数估计图，图12中MDL指的是IL-MeMBer滤波器中在计算每个粒子所在位置处的Music空间谱之前采用MDL算法估计得到的目标个数。从图12可以看出在假目标存在情况下，MDL算法估计得到的目标个数与真实目标个数是存在明显差异的，很明显，IL-MeMBer滤波器能够通过利用真实目标在时间维度上的连续性消除假目标；未经轨迹融合的IL-MeMBer滤波器则会过估计目标个数，其目标个数估计曲线中的谷值是因为目标死亡造成的目标个数估计值骤然下降。IL-MeMBer滤波器在目标个数估计方面同I-MeMBer滤波器几乎一致，但是I-MeMBer滤波器是无法输出每一个目标的轨迹的。

请参见图13，图13为本发明实施例提供的一种不同算法OSPA距离估计图。由图13(b)可以看出在位置估计上未经轨迹融合的IL-MeMBer滤波器在位置估计方面要优于经过轨迹融合后的IL-MeMBer滤波器，造成这一现象的原因是航迹融合技术把属于同一个目标的多条估计得到的轨迹融合为一个目标，显然融合后的目标位置估计误差取决于所有估计值的误差，所提融合后的位置估计误差一定是大于融合前同一个目标的所有估计值的最小误差的。由图13(c)可以看出未经航迹融合的IL-MeMBer滤波器是完全无法对目标的个数进行有效地估计的。由图13(a)则可以看出，IL-MeMBer滤波器和I-MeMBer滤波器的多目标跟踪性能几乎一致，它们的多目标跟踪性能都要优于未经轨迹融合的IL-MeMBer滤波器。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (9)

1.一种基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，包括步骤：

S1、建立目标状态演化模型和多目标叠加量测模型；

S2、将上一时刻的多目标后验概率密度由伯努利参数集表征；

S3、利用所述目标状态演化模型和所述伯努利参数集预测当前时刻的存活伯努利成分，同时获取所述当前时刻的新生标签伯努利成分；

S4、结合所述多目标叠加量测模型，利用Music谱和似然函数更新所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重，得到当前时刻更新后的伯努利成分；

S5、基于欧式距离，对所述更新后的伯努利成分进行融合，得到当前时刻融合后的伯努利成分和融合后的多目标后验概率密度；

S6、根据所述融合后的多目标后验概率密度对目标个数进行估计，根据所述融合后的伯努利成分提取目标状态，得到当前时刻的目标状态估计结果，并返回步骤S3，结合所述融合后的伯努利成分估计下一时刻的目标状态。

2.根据权利要求1所述的基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，所述多目标叠加量测模型为：

Y_k＝[y₁₁,…,y_1N,…,y_N1,…y_NM]^T+N

其中，N为白高斯噪声，y_mn为由第n个阵元发射，经过目标反射后，由第m个阵元接收的经过匹配滤波后的信号；

其中，为第i个目标的多普勒频率，N_k为k时刻的目标个数，ξ_i为第i个目标复散射系数，j为虚数单位，f_n为第n个阵元发射信号的频率，r_i为第i个目标距离，t为时间指标，v_i为第i个目标的速度，λ₀为波长，为由阵元间距引起的发射时延，为由阵元间距引起的接收时延，Δf为频率步进量，n为第n个阵元，f₀为参考频率，d表示各个阵元之间的间距，θ_i表示表示第i个目标相对于雷达的角度，c表示光速。

3.根据权利要求1所述的基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，所述伯努利参数集为：

其中，为k-1时刻第i个Bernoulli成分的存在概率，N_k-1表示k-1时刻的Bernoulli成分总数，k为时间指标，为k-1时刻的空间概率密度，l_k-1,i表示k-1时刻第i个Bernoulli成分的标签；

其中，表示粒子数，为k-1时刻第i个目标状态，为狄拉克脉冲函数，s表示第s个粒子，为刻画第i个目标的第s个粒子，为对应粒子的权重，为一簇加权粒子。

4.根据权利要求1所述的基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，所述存活伯努利成分包括第一存在概率、第一空间分布概率密度和第一粒子权重，其中，

所述第一存在概率为：

其中，为k-1时刻第i个Bernoulli成分的存在概率，为预测后第i个Bernoulli成分中的粒子个数，lk|k-1,i＝lk-1,i表示预测后对应Bernoulli成分所携带的标签不变，P_S,k表示粒子存活并转移的概率，表示k-1时刻刻画第i个目标的第s个粒子，为对应粒子的权重；

所述第一空间分布概率密度为：

其中，为预测后第i个Bernoulli成分中的粒子数，为归一化的粒子权重，为狄拉克脉冲函数，s指示第s个粒子，S表示存活；

所述第一粒子权重为：

其中，为预测后第i个Bernoulli成分中存活的第s个粒子的权重，为粒子数，为粒子权重，为目标状态的转移核，为预测后的粒子状态，为k-1时刻的粒子状态，P_S,k为粒子存活并转移的概率，S表示存活。

5.根据权利要求1所述的基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，所述新生标签伯努利成分包括第二存在概率、第二空间分布概率密度、第二粒子权重和目标标签，其中，

所述第二存在概率为先验给定的第j个新生伯努利成分的存在概率；

所述第二空间分布概率密度为：

其中，为k时刻新生的第j个Bernoulli成分的空间分布概率密度，为初始化粒子集，为第j个新生Bernoulli成分中第s个粒子的权重，l_k|k-1,j为新生Bernoulli成分的标签，为新生第j个Bernoulli成分中的粒子数，Γ表示新生；

所述第二粒子权重为：

其中，为第j个新生Bernoulli成分中第s个粒子的权重，为新生第j个Bernoulli成分中的粒子数；

所述目标标签为：

l_k|k-1,j＝(k,j)

其中，l_k|k-1,j为新生伯努利成分的标签，k表示时刻，j为出现顺序。

6.根据权利要求1所述的基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，步骤S4包括：

S41、结合所述多目标叠加量测模型，利用当前的目标状态预测值和量测值计算Music谱：

其中，克罗内克乘积kron(b(·),a(·))为联合发射-接收域的权矢量，a(·)和b(·)分别表示发射导向矢量和接收导向矢量，P_Music为Music谱，Y_k为多目标叠加量测值，(·)^H表示共轭转置，U_n,k为k时刻的噪声子空间，为将目标状态映射的角度，为将目标状态映射的径向距离；

S42、利用所述Music谱替换传统似然函数，得到新的似然函数：

g(Y_k|·)＝{P_Music(Y_k|·)}^μ

其中，g(Y_k|·)为新的似然函数，P_Music(Y_k|·)为Music谱，μ为u次方；

S43、保持粒子状态不变，利用所述新的似然函数更新所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重：

其中，为k时刻预测的第i个Bernoulli成分中第s个粒子的权重，对应更新后粒子的权重，为新的似然函数，表示第i个Bernoulli成分中第s个粒子从k-1时刻到k时刻的预测状态；

S44、根据所述存活伯努利成分中的粒子权重和所述新生标签伯努利成分中的粒子权重得到所述更新后的伯努利成分：

其中，为更新后的存在概率，为预测的存在概率，为更新后粒子数，更新前后标签保持不变l_k,i＝l_k|k-1,i，l_k,i为更新后第i个Bernoulli成分的标签，为k时刻更新后粒子的权重，为预测的空间分布概率密度，g(Y_k|x)为似然函数，为根据量测Y_k更新后的空间分布概率密度，为狄拉克脉冲函数，s表示第s个粒子；

S45、根据所述更新后的伯努利成分得到更新后的多目标后验概率密度函数：

其中，π_k为k时刻多目标后验概率密度，为更新后Bernoulli成分的存在概率，为更新后Bernoulli成分的空间分布概率密度，为更新后的Bernoulli成分个数，与预测后的Bernoulli成分个数N_k|k-1保持一致。

7.根据权利要求1所述的基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，步骤S5包括：

S51、根据所述更新后的伯努利成分中的目标状态的空间分布概率密度提取目标状态，得到所述目标状态集；

S52、根据所述目标状态集计算欧式距离矩阵，并根据所述欧式距离矩阵和距离阈值划分目标子群；其中，所述欧式距离矩阵为：

其中， 为计算目标状态和之间的欧氏距离，i表示第i个目标，j表示第j个目标，l_k,i表示第k个时刻i个目标的标签，l_k,j表示第k个时刻第j个目标的标签；

S53、将同一个所述目标子群内的目标状态对应的伯努利成分进行融合，得到所述融合后的伯努利成分和所述融合后的多目标后验概率密度。

8.根据权利要求7所述的基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，步骤S53包括：

对同一个所述目标子群内的目标状态对应的粒子集进行融合，得到融合后的粒子集：

其中，和分别为刻画第i个和第j个Bernoulli成分的粒子集；

对同一个所述目标子群内的目标状态对应的粒子权重集进行融合，得到融合后的粒子权重集：

其中，和分别为刻画第i个和第j个Bernoulli成分的粒子对应的权重集，和分别为k时刻第i个和第j个Bernoulli成分中的第s个粒子的权重，和为分别为第i个和第j个Bernoulli成分中粒子的个数；

对同一个所述目标子群内的目标状态对应的目标存在概率集进行融合，得到融合后的目标存在概率：

其中，和分别为第i个和第j个Bernoulli成分的存在概率；

对同一个所述目标子群内的目标状态对应的目标标签集进行融合，得到融合后的目标标签：

其中，为融合后的目标标签，和为不同的标签，k₁和k₂为分别表示时间指标k₁和k₂，i和j为出现顺序；

k时刻所述融合后的多目标后验概率密度为：

其中，π_k为融合后的多目标后验概率密度，为融合后的存在概率，为融合后的空间分布概率密度，N_k表示融合后的目标个数。

9.根据权利要求1所述的基于FDA-MIMO雷达的多目标跟踪方法，其特征在于，根据所述融合后的多目标后验概率密度对目标个数进行估计，根据所述融合后的伯努利成分提取目标状态，得到当前时刻的目标状态估计结果，包括：

S61、根据所述融合后的多目标后验概率密度对目标个数进行估计：

其中，为融合后的前个大存在概率的伯努利成分，N_k为k时刻更新后的Bernoulli成分总数，Ω为所有Bernoulli成分的指标集，和分别为更新后第i和第j个Bernoulli成分的存在概率，Ω_e为Ω的子集，Ω_e表示k时刻假设存在的所有目标对应的指标集合；

S62、从所述融合后的伯努利成分中提取前个目标存在概率的成分，得到所述目标状态：

其中，为粒子状态，为粒子个数，为更新后刻画第i个目标的第s个粒子的状态，为对应粒子的权重；

S63、根据所述目标个数和所述目标状态得到所述目标状态估计结果。