CN104021519A - 基于gpu架构的密集杂波条件下机动多目标跟踪算法 - Google Patents

Abstract

本发明属于雷达与声纳技术领域，主要涉及结合的联合最大似然-交互式多模型-概率数据关联算法(CJML-IMM-PDA)的实现方法，具体来说是一种基于GPD架构的密集杂波条件下多微弱目标轨迹初始化及维持的实现方法，可在微软提供的软件集成开发平台Visual studio上，针对低信噪比、高杂波条件下的激动多微弱目标进行快速的轨迹初始化并且对已初始化成功的轨迹保持跟踪状态。

Description

基于GPU架构的密集杂波条件下机动多目标跟踪算法

技术领域

本发明属于雷达与声纳技术领域，主要涉及结合的联合最大似然_交互式多模型-概率数据关联算法(CJML-IMM-PDA)的实现方法，具体来说是一种基于GPD架构的密集杂波条件下多微弱目标轨迹初始化及维持的实现方法，可在微软提供的软件集成开发平台Visual studio上，针对低信噪比、高杂波条件下的激动多微弱目标进行快速的轨迹初始化并且对已初始化成功的轨迹保持跟踪状态。

背景技术

密集杂波条件下的多微弱目标跟踪一直以来都是多目标跟踪技术领域的研究热点以及难点，并且该技术在雷达(声纳)系统中有举足轻重的作用。当检测区域内的目标杂波密度较大，信噪比较低时，一般使用检测前跟踪(TBD)算法进行跟踪，TBD算法能够累积多帧观测数据，从而更容易找到目标运动的规律，因此能够获得较好的效果。

TBD算法按结构可分为两大类：批处理和迭代，批处理跟踪器直接从多帧观测数据中进行目标跟踪，由于利用了多帧观测，其效果较好，精度较高，但计算量大。迭代算法主要基于贝叶斯理论，常使用粒子滤波(PF)作为具体的实现手段。在工程应用中，在目标跟踪算法实施之前往往需要对航迹进行初始化，以找到目标的初始状态向量，从而进一步进行跟踪，在批处理TBD算法中，多目标轨迹初始化往往选择的是联合最大似然-概率数据关联(JML-PDA)算法。在轨迹初始化完成后，检测便不是首要任务，因此维持目标轨迹可以交给其它计算量相对较小的算法来实施。由于目标在监测区域内的运动状态未知，传统的单模型滤波算法往往无法得到较好的跟踪效果，交互式多模型-联合概率数据关联算法(IMM-JPDA)基于传统的数学模型，能够与TBD算法较好有机结合起来，实现无缝衔接，因此在密集杂波条件下的多目标跟踪体系中，往往将JML-PDA算法与IMM-JPDA算法结合起来使用，称为CJML-IMM-PDA算法。

在CJML-IMM-PDA算法中，首先基于IMM-JPDA算法基于上一时刻已有的轨迹获得当前时刻下的目标状态估计值，再通过JML-PDA算法搜索出当前时刻下的新生轨迹。其中，JML-PDA算法首先需要获得轨迹(包括新生轨迹与已有轨迹)与收到量测之间的关联事件，并基于对多帧观测数据得到的联合总对数似然比(JLLR)进行最大化，在获得JLLR表达式后基于搜索算法获取JLLR最大值后输出对应的参数向量，常用的搜索算法有：网格搜索法、遗传搜索算法(GA)和基于观测空间反映射到参数空间的直接搜索法。其中，网络搜索法由于计算效率低，精度不理想而未被广泛使用，直接搜索法在通过观测映射到参数空间后能够缩小搜索范围以至提高搜索精度，但在高杂波条件下由于接收到的观测值较多而计算量较大，因此其普遍被应用于小杂波情况，相较之下，在单线程架构的处理器上GA算法的搜索精度与直接搜索法相当，并且其搜索基于固定数量的种子，故其计算复杂度在杂波变化时相对稳定，并且由于在搜索的时候GA算法各种子在可独立计算其Fitness指数，具有先天的并行计算优势，因此我们选择GA算法来搜索LLR的最优值，但GA算法也面临着收敛不稳定的问题，经过GA算法搜索出的目标状态向量在较多情况下只是落入最佳JLLR对应状态向量的周围而非最佳状态向量本身，要想进一步搜索得到最佳的目标状态向量，需要进一步搜索，由于GA算法的输出结果离最佳值不远，采用传统的数学搜索方法也不会耗费过多时间。DFP方法是拟牛顿法中的一种，依靠目标函数对被搜索状态向量的偏导来进一步收敛状态向量到最佳值，可取得令人满意的收敛效果。IMM-JPDA算法的主要计算量体现在目标与量测之间关联事件的获取上，会随着目标数与量测数的增加呈指数增加，即“组合爆炸”，在传统架构的处理器中，若监测区域内的目标较多或杂波密度较大，IMM-JPDA算法根本无法实时实施。因此，CJML-IMM-PDA算法的计算量体现在两个方面：一是获得轨迹与量测之间的关联事件；二是基于搜索算法得到最大JLLR对应的新生轨迹状态向量。

发明内容

本发明的目的在于改善传统基于CPU架构的CJML-IMM-PDA算法在计算实时性上的不足，提出一种基于GPU架构的处理方法，在保持与CPU同等计算精度的前提下实现实时跟踪，可直接应用于工程。

本发明的思路是：在处理目标-量测之间的关联问题时，预先将可能存在的各类关联事件存入存储器中，算法在确定当前目标数及量测数后，通过“查表”的方式获得所有的关联方式，通过“离线”计算关联事件可省略较多时间，此外，本发明将关联事件的记录方式由传统的矩阵方式改为向量方式，可实现一维并行计算；在搜索新生轨迹时，采用对杂波数变化相对不太敏感的GA搜索展开，在GPU架构中独立计算每一种子的Fitness指数，通过多线程一次计算即可完成CPU架构中多次循环才能完成的步骤，从而获得加速比。另外，由于JML-PDA采用多帧观测来进行数据拟合，并且每帧观测往往收到多个观测数据，因此JLLR表达式中往往存在含参数的多项式求和，在为每一种子计算Fitness指数时，在CPU架构中往往需要对求得的JLLR式中的多帧似然值进行数次循环再进行求和，造成极大的时间开销，在GPU架构中，考虑通过多个线程计算每个种子的JLLR表达式，其中，单个线程计算一帧观测数据对应的似然值，再通过多线程递归的方式进行求和，可进一步获得加速比。在对单个向量进行复选的时候，参考多线程对应多项式的方式求JLLR的方式可减少求和的循环次数，加速整个“极大似然”过程。

本发明的目的通过如下步骤实现：

S1、在CPU端对IMM-JPDA算法参数进行初始化，具体如下：

S101、初始化观测环境各项参数，所述参数包括：频率观测方差，角度观测方差，距离观测方差，对应频率、角度和距离的观测方程，虚警概率，检测概率，杂波密度，雷达采样间隔，各模型协方差，确认门的门限γ，k-1时刻轨迹在模型j下的状态向量k-1时刻轨迹在模型j下对应的协方差矩阵目标i对应于模型j的概率

S102、收集观测信息，并将其存入观测矩阵Z中，将每一帧手机到的观测数量存入记录向量Obser_num中；

S103、通过调用CUDA runtime驱动库中的cudaMalloc()与cudaMallocPitch()函数在显存中分别开辟向量与矩阵的存储区域，具体存储区域包括：观测存储矩阵cu_Obser、观测数量记录向量cu_Obser_num、目标在模型条件下于k-1时刻的状态矩阵cu_x_k_1、目标在模型条件下的预测状态矩阵cu_x_pre、目标在模型条件下的更新状态矩阵cu_x_k、目标在模型条件下于k-1时刻的协方差矩阵cu_P_k_1、目标在模型条件下的预测协方差矩阵cu_P_pre、目标在模型条件下的更新协方差矩阵cu_P_k、目标在模型条件下的互协方差矩阵cu_S、观测方程的雅可比矩阵cu_H_k_k_1、目标在各模型下的增益矩阵cu_W；

S104、在CPU端将接收到的多帧观测信息矩阵Obser和每帧观测数量Obser_num加载到内存中；

S2、基于马氏转移矩阵P^m在CPU端进行模型交换，计算出目标i对应模型j的预测概率其中，M代表算法在执行时选择的模型集的全部，表示模型从r转移到j的概率，为马氏转移矩阵P^m中的元素；

S3、在CPU端启动多模型预测线程，具体如下：

S301、在模型j条件下的目标状态预测公式为cu_pre＝F_j·cu_x_k_1，其中，F_j代表模型j的转移矩阵；

S302、在模型j条件下的目标协方差矩阵的预测公式为cu_P_pre＝F_j·cu_x_k_1·(F_j)^T+Q_j，其中，(·)^T代表矩阵的转置，Q_j代表模型j的方差；

S303、在模型j条件下的目标互协方差矩阵的计算公式为cu_S＝cu_H_k_k_1·cu_P_pre·(cu_H_k_k_1)^T+R，其中，R为观测的协方差；

S304、在模型j条件下的目标增益矩阵的计算公式为cu_W＝cu_P_pre·cu_H_k_k_1·(cu_P_pre)^-1，其中，(·)^-1表示矩阵求逆；

S4、基于通过离线算得的预设关联事件表获得目标与观测之间的关联事件，同时基于互协方差矩阵较大的模型划定确认门排除一部分关联事件，具体如下：

S401、在S302所述目标协方差矩阵中找出行列式最大的协方差矩阵S_max对应的模型；

S402、找出S303所述行列式最大的协方差矩阵S_max对应的模型；

S403、判断目标i与模型j的关联性，若(z_j)^T·(S_max)^-1·z_j≤γ，则意味着模型j落入目标i确认门中，有关联的可能性，若(z_j)^T·(S_max)^-1·z_j＞γ，则意味着模型j未落入目标i确认门中，没有关联的可能性，其中，确认门的门限γ为经验值；

S5、在CPU端启动模型和目标之间对应概率的计算线程，具体如下：

S501、对第i个关联事件E(i)，在两个模型下计算每个量测与目标关联的似然值，对所述似然值求和；

S502、对所有关联事件，依据事件内量测与目标之间的关联关系加权后归一化得到各量测对于各目标的概率；

S6、在CPU端启动多模型更新线程，具体如下：

S601、基于式cu_x_k＝cu_pre+cu_W·(Z(m)-z_pre(j))算得基于观测m更新得到模型j下目标的状态向量，完成后在模型j下将各观测对应的状态向量同对应的概率加权得到目标在模型j下的状态更新；

S602、基于式cu_P_k＝[cu_x_k-cu_x_k(m)][cu_x_k-cu_x_k(m)]^T算得各量测对应于各模型得到的更新协方差，再将它们同对应的概率加权，最终与预测协方差cu_P_pre相加得到目标在模型j下的更新协方差矩阵；

S7、在CPU端基于目标i对应于模型j的似然值更新模型概率得到

S8、在CPU端计算模型j条件下目标i在k时刻的更新状态向量及协方差矩阵基于将其融合得到最终的目标状态x_k及对应的协方差P_k，其中状态向量融合式为协方差矩阵融合式为 $P^i\left(k\right|k)=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_k^{\mathrm{ij}}\{P_j^i\left(k\right|k)+[x^i\left(k\right|k)-x_j^i\left(k\right|k)\left]\right[x^i\left(k\right|k)-x_j^i\left(k\right|k){]}^T\};$

S9、在CPU端初始化JML-PDA算法参数，具体如下：

S901、初始化GA算法各项参数，包括：随机撒种数量N_p、二进制编码长度、收敛判定阈值、突变概率和交叉概率；

S902、在CPU端随机产生N_p个遗传种子，每个种子含|x|个数值；

S903、通过调用CUDA中的cudaMalloc()函数作为GA算法现存中的向量存储区域，调用CUDA中的cudaMallocPitch()函数作为GA算法显存中的矩阵存储区域，具体为：观测存储矩阵cu_Obser、观测数量记录向量cu_Obser_num、随机父代种子存储矩阵cu_Seed_Old、子代种子存储矩阵cu_Seed_Young、对应种子JLLR存储矩阵cu_Seed_JLLR、计算LLR时需要用到的中间变量存储矩阵cu_JLLR_Per_Observation_Seed、Fitness指数存储向量cu_Fit_Factor、父代配对记录向量cu_Pair_Rec、交叉点记录向量cu_Corss_pos、突变记录矩阵cu_Mut_rec、GA算法收敛结果向量State_Final；

S904、在CPU端将S902所述的N_p个遗传种子通过cudaMemcpy2D()填充显存端的存储区域cu_Seed_Old，并且通过cudaMemcpy2D()将S104所述的Obser拷贝到内存用来填充cu_Obser，通过cudaMemcpy()将S104所述的Obser_num拷贝到内存用来填充cu_Obser_num；

S10、基于k时刻已有轨迹合并新生轨迹与接收到的量测数查表获得关联事件，基于已有轨迹的确认门排除一定数量关联事件；

S11、在CPU端启动JLLR计算线程，具体如下：

S1101、将cu_JLLR_Per_Observation_Seed划分为[max(Obser_num),Np×Nw]的二维矩阵，所述二维矩阵在GPU中以一维向量的方式进行存储，并通过偏移量来进行调用，JLLR值由 $\mathrm{JLLR}=\mathrm{Log}\left\{\underset{i=k}{\overset{k+N_w}{\mathrm{\Π}}}(\underset{E\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Π}}_{t=1}^{K+K^{*}}{\left(P_D^t\right)}^{D(E\left(i\right),t)}{(1-P_D^t)}^{(1-D(E\left(i\right),t))}}{{\mathrm{\λ}}^{d\left[D\left(E\left(i\right)\right)\right]}}.\underset{\underset{e_j\left(i\right)\≠0}{j=1}}{\overset{m\left(i\right)}{\mathrm{\Π}}}{\mathrm{\ρ}}^{e_j\left(i\right)}\left(i\right)p\left(z^j\left(i\right)\right|\left(x^{e_j\left(i\right)}\left(i\right)\right)){|}_{p\left(z^j\left(i\right)\right|\left(x^K\left(i\right)\right))})\right\}$ 算得，其中，D(E(i))＝[d₁(i),d₂(i),...,d_K(i)]为布尔变量构成的向量，所述D(E(i))＝[d₁(i),d₂(i),...,d_K(i)]表示其中任一元素只取0或1，D(E(i),t)表示向量D(E(i))中的第t个元素，表示关联事件E(i)中检测到的目标总数，e_j(i)表示E(i)中与量测j相关联的目标，表示单个目标的似然函数；

S1102、对S1101所述cu_JLLR_Per_Observation_Seed求行和，求得每个遗传种子在对应帧数下的局部联合似然比之和，求得和值之后对每个种子的联合似然比求对数，得到每一帧的局部JLLR；

S1103、根据S1102所述每一帧的局部JLLR求取帧间和得到每个遗传种子的最终LLR值，其中，寻址规则为 $\\begin{array}{c}\mathrm{for\; i}=1:\mathrm{Nw}\\ \mathrm{LLR}={\mathrm{LLR}}_i+{\mathrm{LLR}}_{i+\mathrm{Np}}\\ \mathrm{end}\end{array};$

S12、将cu_Seed_JLLR矩阵复制到内存填充Seed_JLLR，通过遍历方式寻找其最大、最小与均值，计算出Fitness指数线性算式中两个参数

$a=\frac{1}{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\×\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}$

if(a＜0)

$b=\frac{\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\×\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\×\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)};$

else

$b=\frac{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\×\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\×\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}$

S13、在CPU断创建计算Fitness指数的工作线程，通过公式cu_Fit_Factor＝a×JLLR+b计算出每个遗传种子的Fitness指数；

S14、将cu_Fit_Factor拷贝到内存填充Fit_Factor，计算出遗传种子的平均Fitness指数，将所述遗传种子的Fitness指数基于平均数归一化；

S15、基于S14所述归一化后的Fitness指数选择出N_p个父代种子，将所述父代种子随机配对，配对后的结果所对应的遗传种子标签存储于向量Pair_Rec中，将每一个遗传种子交叉点的位置存储于Corss_pos中，将每一个遗传种子的突变记录存储于cu_Mut_rec中；

S16、将S14所述选择出的N_p个父代种子与S15所述标签存储向量Pair_Rec拷贝至GPU端的cu_Seed_Old与cu_Pair_Rec中；

S17、在CPU端启动种子交叉线程，初始化个线程，所述线程寻址为父代种子进行交叉或克隆，得到子代，即，将个线程用于填充向量Corss_pos，若Corss_pos元素等于0，则即不进行交叉进行克隆动作，若Corss_pos元素位于上，则在相应位置进行交叉动作；

S18、在CPU端启动子代突变线程，用突变完成后的子代填充父代存储矩阵，启动N_p线程，若cu_Mut_rec第一行对应值小于阈值，则按照第二行对应位置记录的突变点进行突变，将随机位置上的二进制数取反，若大于阈值，则不突变；

S19、在CPU端启动JLLR计算线程，将得到的新一代种子JLLR值拷贝回内存填充Seed_JLLR，同时得到新一代种子LLR均值mean_JLLR并且进行收敛测试，若Seed_JLLR大于或者等于mean_JLLR的遗传种子个数大于阈值，则收敛，若Seed_JLLR大于或者等于mean_JLLR的遗传种子个数小于或者等于阈值，则不收敛；

S20、若S19所述收敛测试结果为收敛，则GA算法结束，若S19所述收敛测试结果为不收敛，则返回步骤S10；

S21、将S19所述收敛对应的遗传种子复制回内存端，计算所述遗传种子的每一个参数的均值，得到收敛后的向量State_GA；

S22、令k＝1，在CPU端初始化校正矩阵，得到单位矩阵H，其中，H与S21所述向量State_GA维度相同；

S23、在CPU端启动梯度计算线程，计算得到梯度cu_Grad；

S24、在CPU端将S23所述cu_Grad拷贝回内存端用来填充Grad_k向量，循环计算方向向量d＝-H×Grad_k；

S25、在CPU端通过一维搜索方法确定搜索步长，具体为：

S2501、初始化搜索步长λ＝0，设定搜索步长迭代步进λ_step＝1；

S2502、令x_k＝State_GA，在CPU端启动x_k的JLLR计算线程；

S2503、在CPU端将在S2502所述得到的JLLR值拷贝回内存用来填充JLLR_x_k，令JLLR_x_k＝-JLLR_x_k；

S2504、令λ＝λ+λ_step，x_k+1＝State_GA+λ×d，在CPU端启动x_k+1的JLLR计算线程；

S2505、在CPU端将S2504所得到的JLLR值拷贝回内存用来填充JLLR_x_k+1，令JLLR_x_k+1＝-JLLR_x_k+1；

S2506、在CPU端进行判定，若JLLR_x_k+1＞JLLR_x_k则停止，若JLLR_x_k＝JLLR_x_k+1则返回步骤S22；

S26、在CPU端进行收敛判断，若当前||Grad_k+1||小于阈值，则令State_Final＝x_k+1，若当前||Grad_k+1||大于等于阈值，则移至下一个||Grad_k+1||继续进行收敛判断，其中，所述Grad_k+1根据步骤S24中的Grad_k计算得到；

S27、若k＝|x|，则令x₁＝x_k+1，返回步骤S22，否则，更新校正矩阵H，返回步骤S23，其中，校正矩阵H更新式为p＝x_k+1-x_k，q＝Grad_k+1-Grad_k；

S28、输出最优参数向量State_Final，释放CPU和GPU端存储空间；

S29、计算似然比检测值ΛH₁/H₀与航迹校验门限T_vali，验证新生轨迹是否存在，若存在则计算新生轨迹的协方差矩阵，将所述新生轨迹加入已存在轨迹，若不存在则返回步骤S10，具体如下：

S2901、将S28所述最优参数向量State_Final映射到量测空间，并进行量测检验，收集落入关联门限的量测，更新量测矩阵与量测数量记录向量；

S2902、计算对应于State_Final的JLLR值JLLR_State；

S2903、基于算得ΛH₁/H₀值，其中μ₁为航迹存在条件下单帧观测的JLLR期望值，σ₁为航迹存在条件下单帧观测的JLLR标准差，所述μ₁和σ₁在一个跟踪区域内只需计算一次即可，都可以通过在线计算或者都可以通过离线计算；

S2904、根据公式计算门限T_vali；

S2905、判断新生轨迹是否存在，若存在则计算新生轨迹的协方差矩阵，将所述新生轨迹加入已存在轨迹，若不存在则返回步骤S10。

进一步地，S1102所述求得每个遗传种子在对应帧数下的局部联合似然比之和时，需要对max(Obser_num)进行判定：若max(Obser_num)≥64，则通过多线程递归的方式进行求和；若max(Obser_num)＜64，则使用一个线程求和。

进一步地，S24所述循环计算的次数<9。

进一步地，S2903所述计算方式为离线计算。

本发明的有益效果是：

本发明在保持与CPU同等计算精度的前提下实现实时跟踪，可直接应用于工程。

在处理目标-量测之间的关联问题时，预先将可能存在的各类关联事件存入存储器中，算法在确定当前目标数及量测数后，通过“查表”的方式获得所有的关联方式，通过“离线”计算关联事件可省略较多时间，此外，本发明将关联事件的记录方式由传统的矩阵方式改为向量方式，可实现一维并行计算；在搜索新生轨迹时，采用对杂波数变化相对不太敏感的GA搜索展开，在GPU架构中独立计算每一种子的Fitness指数，通过多线程一次计算即可完成CPU架构中多次循环才能完成的步骤，从而获得加速比。另外，由于JML-PDA采用多帧观测来进行数据拟合，并且每帧观测往往收到多个观测数据，因此JLLR表达式中往往存在含参数的多项式求和，在为每一种子计算Fitness指数时，在CPU架构中往往需要对求得的JLLR式中的多帧似然值进行数次循环再进行求和，造成极大的时间开销，在GPU架构中，考虑通过多个线程计算每个种子的JLLR表达式，其中，单个线程计算一帧观测数据对应的似然值，再通过多线程递归的方式进行求和，可进一步获得加速比。在对单个向量进行复选的时候，参考多线程对应多项式的方式求JLLR的方式可减少求和的循环次数，加速整个“极大似然”过程。

附图说明

图1为改变信噪比SNR的值至3dB、6dB、9dB、12dB时算法初始化目标轨迹所用的采样时间以及算法的跟踪精度。

图2为不同杂波密度条件下JML-PDA算法在CPU与GPU平台执行时间对比。

图3为不同种子数量条件下JML-PDA算法在CPU与GPU平台执行时间对比。

具体实施方式

本发明提出一种基于GPU架构的处理方法，具体实现步骤如下：

S1、在CPU端对IMM-JPDA算法参数进行初始化，具体如下：

S101、初始化观测环境各项参数，所述参数包括：频率观测方差，角度观测方差，距离观测方差，对应频率、角度和距离的观测方程，虚警概率，检测概率，杂波密度，雷达采样间隔，各模型协方差，确认门的门限γ，k-1时刻轨迹在模型j下的状态向量k-1时刻轨迹在模型j下对应的协方差矩阵目标i对应于模型j的概率

S102、收集观测信息，并将其存入观测矩阵Z中，将每一帧手机到的观测数量存入记录向量Obser_num中；

S103、通过调用CUDA runtime驱动库中的cudaMalloc()与cudaMallocPitch()函数在显存中分别开辟向量与矩阵的存储区域，具体存储区域包括：观测存储矩阵cu_Obser、观测数量记录向量cu_Obser_num、目标在模型条件下于k-1时刻的状态矩阵cu_x_k_1、目标在模型条件下的预测状态矩阵cu_x_pre、目标在模型条件下的更新状态矩阵cu_x_k、目标在模型条件下于k-1时刻的协方差矩阵cu_P_k_1、目标在模型条件下的预测协方差矩阵cu_P_pre、目标在模型条件下的更新协方差矩阵cu_P_k、目标在模型条件下的互协方差矩阵cu_S、观测方程的雅可比矩阵cu_H_k_k_1、目标在各模型下的增益矩阵cu_W；

S104、在CPU端将接收到的多帧观测信息矩阵Obser和每帧观测数量Obser_num加载到内存中；

S2、基于马氏转移矩阵P^m在CPU端进行模型交换，计算出目标i对应模型j的预测概率其中，M代表算法在执行时选择的模型集的全部，表示模型从r转移到j的概率，为马氏转移矩阵P^m中的元素；

S3、在CPU端启动多模型预测线程，具体如下：

S301、在模型j条件下的目标状态预测公式为cu_pre＝F_j·cu_x_k_1，其中，F_j代表模型j的转移矩阵；

S302、在模型j条件下的目标协方差矩阵的预测公式为cu_P_pre＝F_j·cu_x_k_1·(F_j)^T+Q_j，其中，(·)^T代表矩阵的转置，Q_j代表模型j的方差；

S303、在模型j条件下的目标互协方差矩阵的计算公式为cu_S＝cu_H_k_k_1·cu_P_pre·(cu_H_k_k_1)^T+R，其中，R为观测的协方差；

S304、在模型j条件下的目标增益矩阵的计算公式为cu_W＝cu_P_pre·cu_H_k_k_1·(cu_P_pre)^-1，其中，(·)^-1表示矩阵求逆；

S4、基于通过离线算得的预设关联事件表获得目标与观测之间的关联事件，同时基于互协方差矩阵较大的模型划定确认门排除一部分关联事件，具体如下：

S401、在S302所述目标协方差矩阵中找出行列式最大的协方差矩阵S_max对应的模型；

S402、找出S303所述行列式最大的协方差矩阵S_max对应的模型；

S403、判断目标i与模型j的关联性，若(z_j)^T·(S_max)^-1·z_j≤γ，则意味着模型j落入目标i确认门中，有关联的可能性，若(z_j)^T·(S_max)^-1·z_j＞γ，则意味着模型j未落入目标i确认门中，没有关联的可能性，其中，确认门的门限γ为经验值；

S5、在CPU端启动模型和目标之间对应概率的计算线程，具体如下：

S501、对第i个关联事件E(i)，在两个模型下计算每个量测与目标关联的似然值，对所述似然值求和；

S502、对所有关联事件，依据事件内量测与目标之间的关联关系加权后归一化得到各量测对于各目标的概率；

S6、在CPU端启动多模型更新线程，具体如下：

S601、基于式cu_x_k＝cu_pre+cu_W·(Z(m)-z_pre(j))算得基于观测m更新得到模型j下目标的状态向量，完成后在模型j下将各观测对应的状态向量同对应的概率加权得到目标在模型j下的状态更新；

S602、基于式cu_P_k＝[cu_x_k-cu_x_k(m)][cu_x_k-cu_x_k(m)]^T算得各量测对应于各模型得到的更新协方差，再将它们同对应的概率加权，最终与预测协方差cu_P_pre相加得到目标在模型j下的更新协方差矩阵；

S7、在CPU端基于目标i对应于模型j的似然值更新模型概率得到具体为：

S701、根据公式计算模型j对应目标i的似然值其中， $b=\frac{\mathrm{\&lambda;}(1-P_D{P}_G)}{P_D{P}_G},$ $P_j^t\left(k\right)=\left\{\begin{array}{cc}0,& {\mathrm{\&omega;}}_{\mathrm{it}}\&NotEqual;1\\ N({\tilde{Z}}_j^t;0,S^t)P_D^t& \mathrm{others}\end{array}\right.,$ V_k代表监测区域面积，P_G表示正确量测落入确认门的概率，N(·)表示高斯分布概率密度表达式；

S702、根据公式对进行更新，在得到后，通过式归一化得到各模型对应各目标的更新概率，其中，

S8、在CPU端计算模型j条件下目标i在k时刻的更新状态向量及协方差矩阵基于将其融合得到最终的目标状态x_k及对应的协方差P_k，其中状态向量融合式为协方差矩阵融合式为 $P^i\left(k\right|k)=\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_k^{\mathrm{ij}}\{P_j^i\left(k\right|k)+[x^i\left(k\right|k)-x_j^i\left(k\right|k)\left]\right[x^i\left(k\right|k)-x_j^i\left(k\right|k){]}^T\};$

S9、在CPU端初始化JML-PDA算法参数，具体如下：

S901、初始化GA算法各项参数，包括：随机撒种数量N_p、二进制编码长度、收敛判定阈值、突变概率和交叉概率；

S902、在CPU端随机产生N_p个遗传种子，每个种子含|x|个数值；

S903、通过调用CUDA中的cudaMalloc()函数作为GA算法现存中的向量存储区域，调用CUDA中的cudaMallocPitch()函数作为GA算法显存中的矩阵存储区域，具体为：观测存储矩阵cu_Obser、观测数量记录向量cu_Obser_num、随机父代种子存储矩阵cu_Seed_Old、子代种子存储矩阵cu_Seed_Young、对应种子JLLR存储矩阵cu_Seed_JLLR、计算LLR时需要用到的中间变量存储矩阵cu_JLLR_Per_Observation_Seed、Fitness指数存储向量cu_Fit_Factor、父代配对记录向量cu_Pair_Rec、交叉点记录向量cu_Corss_pos、突变记录矩阵cu_Mut_rec、GA算法收敛结果向量State_Final；

S904、在CPU端将S902所述的N_p个遗传种子通过cudaMemcpy2D()填充显存端的存储区域cu_Seed_Old，并且通过cudaMemcpy2D()将S104所述的Obser拷贝到内存用来填充cu_Obser，通过cudaMemcpy()将S104所述的Obser_num拷贝到内存用来填充cu_Obser_num；

S10、基于k时刻已有轨迹合并新生轨迹与接收到的量测数查表获得关联事件，基于已有轨迹的确认门排除一定数量关联事件；

S11、在CPU端启动JLLR计算线程，具体如下：

S1101、将cu_JLLR_Per_Observation_Seed划分为[max(Obser_num),Np×Nw]的二维矩阵，所述二维矩阵在GPU中以一维向量的方式进行存储，并通过偏移量来进行调用，JLLR值由 $\mathrm{JLLR}=\mathrm{Log}\left\{\underset{i=k}{\overset{k+N_w}{\mathrm{\&Pi;}}}(\underset{E\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Pi;}}_{t=1}^{K+K^{*}}{\left(P_D^t\right)}^{D(E\left(i\right),t)}{(1-P_D^t)}^{(1-D(E\left(i\right),t))}}{{\mathrm{\&lambda;}}^{d\left[D\left(E\left(i\right)\right)\right]}}.\underset{\underset{e_j\left(i\right)\&NotEqual;0}{j=1}}{\overset{m\left(i\right)}{\mathrm{\&Pi;}}}{\mathrm{\&rho;}}^{e_j\left(i\right)}\left(i\right)p\left(z^j\left(i\right)\right|\left(x^{e_j\left(i\right)}\left(i\right)\right)){|}_{p\left(z^j\left(i\right)\right|\left(x^K\left(i\right)\right))})\right\}$ 算得，其中，D(E(i))＝[d₁(i),d₂(i),...,d_K(i)]为布尔变量构成的向量，所述D(E(i))＝[d₁(i),d₂(i),...,d_K(i)]表示其中任一元素只取0或1，例如：d_K(i)＝1，则表示在关联事件E(i)中的第K个目标被检测到，相应地，d_K(i)＝0即表示未被检测，D(E(i),t)表示向量D(E(i))中的第t个元素，表示关联事件E(i)中检测到的目标总数，e_j(i)表示E(i)中与量测j相关联的目标， $p\left(z^j\left(i\right)\right|\left(x^{e_j\left(i\right)}\left(i\right)\right))=N(z;H^k(x^{e_j\left(i\right)}\left(i\right),x^s\left(i\right)),R)$ 表示单个目标的似然函数；

S1102、对S1101所述cu_JLLR_Per_Observation_Seed求行和，求得每个遗传种子在对应帧数下的局部联合似然比之和，求得和值之后对每个种子的联合似然比求对数，得到每一帧的局部JLLR，所述求得每个遗传种子在对应帧数下的局部联合似然比之和时，需要对max(Obser_num)进行判定：若max(Obser_num)≥64，则通过多线程递归的方式进行求和；若max(Obser_num)＜64，则使用一个线程求和。例如：初始化4个线程，每个线程负责16个数相加，则只需要循环16次得到4个和值，最终再将四个和值相加即可求得总和，一共需要17个步骤，可实现求和的加速，在max(Obser_num)数量不多时无需采用这种方式，使用一个线程求和即可。

S1103、根据S1102所述每一帧的局部JLLR求取帧间和得到每个遗传种子的最终LLR值，其中，寻址规则为 $\\begin{array}{c}\mathrm{for\; i}=1:\mathrm{Nw}\\ \mathrm{LLR}={\mathrm{LLR}}_i+{\mathrm{LLR}}_{i+\mathrm{Np}}\\ \mathrm{end}\end{array};$

S12、将cu_Seed_JLLR矩阵复制到内存填充Seed_JLLR，通过遍历方式寻找其最大、最小与均值，计算出Fitness指数线性算式中两个参数

$a=\frac{1}{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}$

if(a＜0)

$b=\frac{\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)};$

else

$b=\frac{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}$

S13、在CPU断创建计算Fitness指数的工作线程，通过公式cu_Fit_Factor＝a×JLLR+b计算出每个遗传种子的Fitness指数；

S14、将cu_Fit_Factor拷贝到内存填充Fit_Factor，计算出遗传种子的平均Fitness指数，将所述遗传种子的Fitness指数基于平均数归一化；

S15、基于S14所述归一化后的Fitness指数选择出N_p个父代种子，将所述父代种子随机配对，配对后的结果所对应的遗传种子标签存储于向量Pair_Rec中，将每一个遗传种子交叉点的位置存储于Corss_pos中，将每一个遗传种子的突变记录存储于cu_Mut_rec中，具体如下：

S1501、基于S14所述归一化后的Fitness指数选择出N_p个父代种子，将所述父代种子随机配对，配对后的结果所对应的遗传种子标签存储于向量Pair_Rec中，将Pair_Rec填充为1:N_p的N_p个元素，即，循环N_p次，每一次随机生成一个[1,N_p]之间的整数Mov_step，将Pair_Rec右移Mov_step次；

S1502、S1501所述的生成N_p个[0,l]之间的随机整数并将所述随机整数用于填充向量Corss_pos，即，若Corss_pos元素等于0，则即不进行交叉进行克隆动作，若Corss_pos元素位于[1,l]上，则在相应位置进行交叉动作，其中，l为经验值；

S1503、生成N_p个位于[0,1]的均匀分布的随机数，并且在每一个位置进行判断，若其大于突变概率P_mut，则在第二行对应位置生成位于[1,l]间的随机整数；

S16、将S14所述选择出的N_p个父代种子与S15所述标签存储向量Pair_Rec拷贝至GPU端的cu_Seed_Old与cu_Pair_Rec中；

S17、在CPU端启动种子交叉线程，初始化个线程，所述线程寻址为父代种子进行交叉或克隆，得到子代，即，将个线程用于填充向量Corss_pos，若Corss_pos元素等于0，则即不进行交叉进行克隆动作，若Corss_pos元素位于上，则在相应位置进行交叉动作；

S18、在CPU端启动子代突变线程，用突变完成后的子代填充父代存储矩阵，启动N_p线程，若cu_Mut_rec第一行对应值小于阈值，则按照第二行对应位置记录的突变点进行突变，将随机位置上的二进制数取反，若大于阈值，则不突变；

S19、在CPU端启动JLLR计算线程，将得到的新一代种子JLLR值拷贝回内存填充Seed_JLLR，同时得到新一代种子LLR均值mean_JLLR并且进行收敛测试，若Seed_JLLR大于或者等于mean_JLLR的遗传种子个数大于阈值，则收敛，若Seed_JLLR大于或者等于mean_JLLR的遗传种子个数小于或者等于阈值，则不收敛；

S20、若S19所述收敛测试结果为收敛，则GA算法结束，若S19所述收敛测试结果为不收敛，则返回步骤S10；

S21、将S19所述收敛对应的遗传种子复制回内存端，计算所述遗传种子的每一个参数的均值，得到收敛后的向量State_GA；

S22、令k＝1，在CPU端初始化校正矩阵，得到单位矩阵H，其中，H与S21所述向量State_GA维度相同；

S23、在CPU端启动梯度计算线程，计算得到梯度cu_Grad；

S24、在CPU端将S23所述cu_Grad拷贝回内存端用来填充Grad_k向量，循环计算方向向量d＝-H×Grad_k，所述循环计算的次数<9；

S25、在CPU端通过一维搜索方法确定搜索步长，具体为：

S2501、初始化搜索步长λ＝0，设定搜索步长迭代步进λ_step＝1；

S2502、令x_k＝State_GA，在CPU端启动x_k的JLLR计算线程；

S2503、在CPU端将在S2502所述得到的JLLR值拷贝回内存用来填充JLLR_x_k，令JLLR_x_k＝-JLLR_x_k；

S2504、令λ＝λ+λ_step，x_k+1＝State_GA+λ×d，在CPU端启动x_k+1的JLLR计算线程；

S2505、在CPU端将S2504所得到的JLLR值拷贝回内存用来填充JLLR_x_k+1，令JLLR_x_k+1＝-JLLR_x_k+1；

S2506、在CPU端进行判定，若JLLR_x_k+1＞JLLR_x_k则停止，若JLLR_x_k＝JLLR_x_k+1则返回步骤S22；

S26、在CPU端进行收敛判断，若当前||Grad_k+1||小于阈值，则令State_Final＝x_k+1，若当前||Grad_k+1||大于等于阈值，则移至下一个||Grad_k+1||继续进行收敛判断，其中，所述Grad_k+1根据步骤S24中的Grad_k计算得到；

S27、若k＝|x|，则令x₁＝x_k+1，返回步骤S22，否则，更新校正矩阵H，返回步骤S23，其中，校正矩阵H更新式为p＝x_k+1-x_k，q＝Grad_k+1-Grad_k；

S28、输出最优参数向量State_Final，释放CPU和GPU端存储空间；

S29、计算似然比检测值ΛH₁/H₀与航迹校验门限T_vali，验证新生轨迹是否存在，若存在则计算新生轨迹的协方差矩阵，将所述新生轨迹加入已存在轨迹，若不存在则返回步骤S10，具体如下：

S2901、将S28所述最优参数向量State_Final映射到量测空间，并进行量测检验，收集落入关联门限的量测，更新量测矩阵与量测数量记录向量；

S2902、计算对应于State_Final的JLLR值JLLR_State；

S2903、基于算得ΛH₁/H₀值，其中μ₁为航迹存在条件下单帧观测的JLLR期望值，σ₁为航迹存在条件下单帧观测的JLLR标准差，所述μ₁和σ₁在一个跟踪区域内通过离线计算，计算一次即可；

S2904、根据公式计算门限T_vali；

S2905、判断新生轨迹是否存在，若存在则计算新生轨迹的协方差矩阵，将所述新生轨迹加入已存在轨迹，若不存在则返回步骤S10。

本发明的效果可以通过以下仿真进一步说明：

(1)实验条件

CPU架构工作条件：

硬件：

CPU：Core^TMi3-2120CPU3.3GHz内存

内存：3.00GB

软件：Windows7SP132位操作系统，Matlab2012b

GPU架构工作条件：

硬件：

GPU：Nvidia GeForce GT610、48Cores1.62GHz

显存：1.00GB

软件：

Visual Studio2010,CUDA Driver Version:5.0

(2)实验内容及结果：

在一个二维坐标的声纳应用场景中，传感器一共观测了70个采样时刻，其中，传感器在第1时刻由[9000m,0m]位置开始运动，其速度为[-6m/s,1m/s]，在第51个采样时刻改变运动方向，其速度向量变为[6,1]。在采样过程中，有两个目标先后出现监测区域中，其中，目标T1在在第1时刻出现，其初始状态的运动向量为：[2000m,10m/s,11000m,10m/s]^T，当出现后便沿着初始速度做匀速直线运动，直到第40时刻发生机动，开始沿逆时针方向作角速度为0.9°/s的圆周运动，至50时刻恢复做匀速直线运动至70时刻；目标T2在第20时刻出现，其初始状态向量为：[9000m,-10m/s,14000m,10m/s]^T，第50时刻发生机动，沿逆时针方向做角速度为0.9°/s的圆周运动，至60时刻恢复匀速直线运动。两个目标的状态转移方程噪声功率均为0.01，目标的存活概率P_s,k＝0.99。

声纳传感器接收目标的角度、多谱勒频率及回波的幅度信息，其采样周期T＝20s。使用IMM-JPDA算法来维持目标的运动轨迹，一共包括两个模型：匀速直线运动模型F_CV和匀转弯运动模型F_CT，它们分别为：

$F_{\mathrm{CV}}=\left(\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right)$

$F_{\mathrm{CT}}=\left[\begin{array}{ccccc}1& \frac{\sin \left(\mathrm{\&omega;T}\right)}{\sin \left(\mathrm{\&omega;}\right)}& 0& -\frac{(1-\cos \left(\mathrm{\&omega;T}\right))}{\sin \left(\mathrm{\&omega;}\right)}& 0\\ 0& \cos \left(\mathrm{\&omega;T}\right)& 0& -\sin \left(\mathrm{\&omega;T}\right)& 0\\ 0& \frac{(1-\cos \left(\mathrm{\&omega;T}\right))}{\sin \left(\mathrm{\&omega;}\right)}& 1& \frac{\sin \left(\mathrm{\&omega;T}\right)}{\sin \left(\mathrm{\&omega;}\right)}& 0\\ 0& \sin \left(\mathrm{\&omega;T}\right)& 0& \cos \left(\mathrm{\&omega;T}\right)& 0\\ 0& 0& 0& 0& 1\end{array}\right]$

其中转弯模型的角速度ω＝0.6°/s。目标的观测方程为：

$z\left(k\right)=\left[\begin{array}{c}a\tan \left(\frac{x^t-x^s}{y^t-y^s}\right)\\ \mathrm{\&gamma;}[1-\frac{{\stackrel{.}{x}}^t\sin \mathrm{\&theta;}+{\stackrel{.}{y}}^t\cos \mathrm{\&theta;}}{c}]\end{array}\right]$

式中，γ代表声波频率，c代表声波在介质中的传输速度，角度观测的噪声标准差σ_θ＝0.5π/180deg，频率观测的噪声标准差σ_γ＝1Hz。传感器接收到的回波幅度a呈瑞利分布，当回波来源于目标时，其概率密度方程表示为p₁(a)，若回波是杂波，其概率密度分布被表示为p₀(a)，它们的具体形式如下：

$p_1\left(a\right)=\frac{a}{1+d}\exp (-\frac{a^2}{2(1+d)}),a\&GreaterEqual;0$

$p_0\left(a\right)=\mathrm{aexp}(-\frac{a^2}{2}),a\&GreaterEqual;0$

式中，d代表监测环境的信噪比，因此，若声纳对信号的检测门限为τ，那么传感器对目标的检测概率P_D及虚警概率P_FA计算公式为：

$P_D=\exp (-\frac{{\mathrm{\&tau;}}^2}{2(1+d)})$

$P_{\mathrm{FA}}=\exp (-\frac{{\mathrm{\&tau;}}^2}{2})$

通过下式算得通过检测门限后信号的幅度分布及幅度似然比ρ：

$p_0^{\mathrm{\&tau;}}\left(a\right)=\frac{1}{P_{\mathrm{FA}}}\mathrm{aexp}(-\frac{a^2}{2}),a>\mathrm{\&tau;}$

$p_1^{\mathrm{\&tau;}}\left(a\right)=\frac{1}{P_D}\frac{a}{1+d}\exp (-\frac{a^2}{2(1+d)}),a>\mathrm{\&tau;}$

$\mathrm{\&rho;}=\frac{p_1^{\mathrm{\&tau;}}\left(a\right)}{p_0^{\mathrm{\&tau;}}\left(a\right)}$

在本仿真条件中，τ＝2.64，传感器的角度观测范围为[0,π]，声波频率的接收范围为[290,310]，发射频率为300Hz。

在使用JML-PDA算法负责目标航迹起始时，观测滑窗内使用10帧观测数据构建目标JLLR，使用GA搜索算法对目标的最大JLLR及对应的状态向量进行粗略搜索，GA算法一共生成100个种子，种子间交叉概率P_corss＝0.8，突变概率P_mut＝0.01，对种子使用32位二进制编码，DFP算法收敛条件设为0.001。

为了更加全面地说明本发明的各项指标，现基于以上仿真条件做如下实验：

实验一：跟踪精度测试。分别改变信噪比SNR的值至3dB、6dB、9dB、12dB时算法初始化目标轨迹所用的采样时间以及算法的跟踪精度分别如表1和图1所示。

信噪比(dB)	3	6	9	12
					目标1初始化时间(采样时刻)	1.31	1.16	1.00	1.00
目标2初始化时间(采样时刻)	1.28	1.17	1.00	1.00
					目标1初始化误差(米)	97.28	69.44	42.63	28.55
目标2初始化误差(米)	100.87	72.68	44.01	36.69

表1

实验二：杂波密度与计算量测试。分别改变监测区域内的杂波密度至0.1/deg·Hz、1/deg·Hz、2/deg·Hz、5/deg·Hz以及10/deg·Hz时在GPU平台和CPU平台上对应的算法运行时间如图2所示。

实验三：搜索算法种子数与计算量测试。分别改变GA算法搜索种子数至100、200、300、400、500及600时算法在GPU平台和CPU平台上对应的运行时间如图3所示。

Claims (4)

1.基于GPU架构的密集杂波条件下机动多目标跟踪算法，其特征在于，包括以下步骤：

S1、在CPU端对IMM-JPDA算法参数进行初始化，具体如下：

S101、初始化观测环境各项参数，所述参数包括：频率观测方差，角度观测方差，距离观测方差，对应频率、角度和距离的观测方程，虚警概率，检测概率，杂波密度，雷达采样间隔，各模型协方差，确认门的门限γ，k-1时刻轨迹在模型j下的状态向量k-1时刻轨迹在模型j下对应的协方差矩阵目标i对应于模型j的概率

S102、收集观测信息，并将其存入观测矩阵Z中，将每一帧手机到的观测数量存入记录向量Obser_num中；

S103、通过调用CUDA runtime驱动库中的cudaMalloc()与cudaMallocPitch()函数在显存中分别开辟向量与矩阵的存储区域，具体存储区域包括：观测存储矩阵cu_Obser、观测数量记录向量cu_Obser_num、目标在模型条件下于k-1时刻的状态矩阵cu_x_k_1、目标在模型条件下的预测状态矩阵cu_x_pre、目标在模型条件下的更新状态矩阵cu_x_k、目标在模型条件下于k-1时刻的协方差矩阵cu_P_k_1、目标在模型条件下的预测协方差矩阵cu_P_pre、目标在模型条件下的更新协方差矩阵cu_P_k、目标在模型条件下的互协方差矩阵cu_S、观测方程的雅可比矩阵cu_H_k_k_1、目标在各模型下的增益矩阵cu_W；

S104、在CPU端将接收到的多帧观测信息矩阵Obser和每帧观测数量Obser_num加载到内存中；

S2、基于马氏转移矩阵P^m在CPU端进行模型交换，计算出目标i对应模型j的预测概率其中，M代表算法在执行时选择的模型集的全部，表示模型从r转移到j的概率，为马氏转移矩阵P^m中的元素；

S3、在CPU端启动多模型预测线程，具体如下：

S301、在模型j条件下的目标状态预测公式为cu_pre＝F_j·cu_x_k_1，其中，F_j代表模型j的转移矩阵；

S302、在模型j条件下的目标协方差矩阵的预测公式为cu_P_pre＝F_j·cu_x_k_1·(F_j)^T+Q_j，其中，(·)^T代表矩阵的转置，Q_j代表模型j的方差；

S303、在模型j条件下的目标互协方差矩阵的计算公式为cu_S＝cu_H_k_k_1·cu_P_pre·(cu_H_k_k_1)^T+R，其中，R为观测的协方差；

S304、在模型j条件下的目标增益矩阵的计算公式为cu_W＝cu_P_pre·cu_H_k_k_1·(cu_P_pre)^-1，其中，(·)^-1表示矩阵求逆；

S4、基于通过离线算得的预设关联事件表获得目标与观测之间的关联事件，同时基于互协方差矩阵较大的模型划定确认门排除一部分关联事件，具体如下：

S401、在S302所述目标协方差矩阵中找出行列式最大的协方差矩阵S_max对应的模型；

S402、找出S303所述行列式最大的协方差矩阵S_max对应的模型；

S403、判断目标i与模型j的关联性，若(z_j)^T·(S_max)^-1·z_j≤γ，则意味着模型j落入目标i确认门中，有关联的可能性，若(z_j)^T·(S_max)^-1·z_j＞γ，则意味着模型j未落入目标i确认门中，没有关联的可能性，其中，确认门的门限γ为经验值；

S5、在CPU端启动模型和目标之间对应概率的计算线程，具体如下：

S501、对第i个关联事件E(i)，在两个模型下计算每个量测与目标关联的似然值，对所述似然值求和；

S502、对所有关联事件，依据事件内量测与目标之间的关联关系加权后归一化得到各量测对于各目标的概率；

S6、在CPU端启动多模型更新线程，具体如下：

S601、基于式cu_x_k＝cu_pre+cu_W·(Z(m)-z_pre(j))算得基于观测m更新得到模型j下目标的状态向量，完成后在模型j下将各观测对应的状态向量同对应的概率加权得到目标在模型j下的状态更新；

S602、基于式cu_P_k＝[cu_x_k-cu_x_k(m)][cu_x_k-cu_x_k(m)]^T算得各量测对应于各模型得到的更新协方差，再将它们同对应的概率加权，最终与预测协方差cu_P_pre相加得到目标在模型j下的更新协方差矩阵；

S7、在CPU端基于目标i对应于模型j的似然值更新模型概率得到

S8、在CPU端计算模型j条件下目标i在k时刻的更新状态向量及协方差矩阵基于将其融合得到最终的目标状态x_k及对应的协方差P_k，其中状态向量融合式为协方差矩阵融合式为 $P^i\left(k\right|k)=\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_k^{\mathrm{ij}}\{P_j^i\left(k\right|k)+[x^i\left(k\right|k)-x_j^i\left(k\right|k)\left]\right[x^i\left(k\right|k)-x_j^i\left(k\right|k){]}^T\};$

S9、在CPU端初始化JML-PDA算法参数，具体如下：

S901、初始化GA算法各项参数，包括：随机撒种数量N_p、二进制编码长度、收敛判定阈值、突变概率和交叉概率；

S902、在CPU端随机产生N_p个遗传种子，每个种子含|x|个数值；

S903、通过调用CUDA中的cudaMalloc()函数作为GA算法现存中的向量存储区域，调用CUDA中的cudaMallocPitch()函数作为GA算法显存中的矩阵存储区域，具体为：观测存储矩阵cu_Obser、观测数量记录向量cu_Obser_num、随机父代种子存储矩阵cu_Seed_Old、子代种子存储矩阵cu_Seed_Young、对应种子JLLR存储矩阵cu_Seed_JLLR、计算LLR时需要用到的中间变量存储矩阵cu_JLLR_Per_Observation_Seed、Fitness指数存储向量cu_Fit_Factor、父代配对记录向量cu_Pair_Rec、交叉点记录向量cu_Corss_pos、突变记录矩阵cu_Mut_rec、GA算法收敛结果向量State_Final；

S904、在CPU端将S902所述的N_p个遗传种子通过cudaMemcpy2D()填充显存端的存储区域cu_Seed_Old，并且通过cudaMemcpy2D()将S104所述的Obser拷贝到内存用来填充cu_Obser，通过cudaMemcpy()将S104所述的Obser_num拷贝到内存用来填充cu_Obser_num；

S10、基于k时刻已有轨迹合并新生轨迹与接收到的量测数查表获得关联事件，基于已有轨迹的确认门排除一定数量关联事件；

S11、在CPU端启动JLLR计算线程，具体如下：

S1101、将cu_JLLR_Per_Observation_Seed划分为[max(Obser_num),Np×Nw]的二维矩阵，所述二维矩阵在GPU中以一维向量的方式进行存储，并通过偏移量来进行调用，JLLR值由 $\mathrm{JLLR}=\mathrm{Log}\left\{\underset{i=k}{\overset{k+N_w}{\mathrm{\&Pi;}}}(\underset{E\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Pi;}}_{t=1}^{K+K^{*}}{\left(P_D^t\right)}^{D(E\left(i\right),t)}{(1-P_D^t)}^{(1-D(E\left(i\right),t))}}{{\mathrm{\&lambda;}}^{d\left[D\left(E\left(i\right)\right)\right]}}.\underset{\underset{e_j\left(i\right)\&NotEqual;0}{j=1}}{\overset{m\left(i\right)}{\mathrm{\&Pi;}}}{\mathrm{\&rho;}}^{e_j\left(i\right)}\left(i\right)p\left(z^j\left(i\right)\right|\left(x^{e_j\left(i\right)}\left(i\right)\right)){|}_{p\left(z^j\left(i\right)\right|\left(x^K\left(i\right)\right))})\right\}$ 算得，其中，D(E(i))＝[d₁(i),d₂(i),...,d_K(i)]为布尔变量构成的向量，所述D(E(i))＝[d₁(i),d₂(i),...,d_K(i)]表示其中任一元素只取0或1，D(E(i),t)表示向量D(E(i))中的第t个元素，表示关联事件E(i)中检测到的目标总数，e_j(i)表示E(i)中与量测j相关联的目标，表示单个目标的似然函数；

S1102、对S1101所述cu_JLLR_Per_Observation_Seed求行和，求得每个遗传种子在对应帧数下的局部联合似然比之和，求得和值之后对每个种子的联合似然比求对数，得到每一帧的局部JLLR；

S1103、根据S1102所述每一帧的局部JLLR求取帧间和得到每个遗传种子的最终LLR值，其中，寻址规则为 $\\begin{array}{c}\mathrm{for\; i}=1:\mathrm{Nw}\\ \mathrm{LLR}={\mathrm{LLR}}_i+{\mathrm{LLR}}_{i+\mathrm{Np}}\\ \mathrm{end}\end{array};$

S12、将cu_Seed_JLLR矩阵复制到内存填充Seed_JLLR，通过遍历方式寻找其最大、最小与均值，计算出Fitness指数线性算式中两个参数

$a=\frac{1}{\max \left(\mathrm{JLLE}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLE}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}$

if(a＜0)

$b=\frac{\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)};$

else

$b=\frac{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}{\max \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)+\min \left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)-2\&times;\mathrm{mean}\left(\mathrm{JLLR}\right[1:\mathrm{Np}\left]\right)}$

S13、在CPU断创建计算Fitness指数的工作线程，通过公式cu_Fit_Factor＝a×JLLR+b计算出每个遗传种子的Fitness指数；

S14、将cu_Fit_Factor拷贝到内存填充Fit_Factor，计算出遗传种子的平均Fitness指数，将所述遗传种子的Fitness指数基于平均数归一化；

S15、基于S14所述归一化后的Fitness指数选择出N_p个父代种子，将所述父代种子随机配对，配对后的结果所对应的遗传种子标签存储于向量Pair_Rec中，将每一个遗传种子交叉点的位置存储于Corss_pos中，将每一个遗传种子的突变记录存储于cu_Mut_rec中；

S16、将S14所述选择出的N_p个父代种子与S15所述标签存储向量Pair_Rec拷贝至GPU端的cu_Seed_Old与cu_Pair_Rec中；

S17、在CPU端启动种子交叉线程，初始化个线程，所述线程寻址为父代种子进行交叉或克隆，得到子代，即，将个线程用于填充向量Corss_pos，若Corss_pos元素等于0，则即不进行交叉进行克隆动作，若Corss_pos元素位于上，则在相应位置进行交叉动作；

S18、在CPU端启动子代突变线程，用突变完成后的子代填充父代存储矩阵，启动N_p线程，若cu_Mut_rec第一行对应值小于阈值，则按照第二行对应位置记录的突变点进行突变，将随机位置上的二进制数取反，若大于阈值，则不突变；

S19、在CPU端启动JLLR计算线程，将得到的新一代种子JLLR值拷贝回内存填充Seed_JLLR，同时得到新一代种子LLR均值mean_JLLR并且进行收敛测试，若Seed_JLLR大于或者等于mean_JLLR的遗传种子个数大于阈值，则收敛，若Seed_JLLR大于或者等于mean_JLLR的遗传种子个数小于或者等于阈值，则不收敛；

S20、若S19所述收敛测试结果为收敛，则GA算法结束，若S19所述收敛测试结果为不收敛，则返回步骤S10；

S21、将S19所述收敛对应的遗传种子复制回内存端，计算所述遗传种子的每一个参数的均值，得到收敛后的向量State_GA；

S22、令k＝1，在CPU端初始化校正矩阵，得到单位矩阵H，其中，H与S21所述向量State_GA维度相同；

S23、在CPU端启动梯度计算线程，计算得到梯度cu_Grad；

S24、在CPU端将S23所述cu_Grad拷贝回内存端用来填充Grad_k向量，循环计算方向向量d＝-H×Grad_k；

S25、在CPU端通过一维搜索方法确定搜索步长，具体为：

S2501、初始化搜索步长λ＝0，设定搜索步长迭代步进λ_step＝1；

S2502、令x_k＝State_GA，在CPU端启动x_k的JLLR计算线程；

S2503、在CPU端将在S2502所述得到的JLLR值拷贝回内存用来填充JLLR_x_k，令JLLR_x_k＝-JLLR_x_k；

S2504、令λ＝λ+λ_step，x_k+1＝State_GA+λ×d，在CPU端启动x_k+1的JLLR计算线程；

S2505、在CPU端将S2504所得到的JLLR值拷贝回内存用来填充JLLR_x_k+1，令JLLR_x_k+1＝-JLLR_x_k+1；

S2506、在CPU端进行判定，若JLLR_x_k+1＞JLLR_x_k则停止，若JLLR_x_k＝JLLR_x_k+1则返回步骤S22；

S26、在CPU端进行收敛判断，若当前||Grad_k+1||小于阈值，则令State_Final＝x_k+1，若当前||Grad_k+1||大于等于阈值，则移至下一个||Grad_k+1||继续进行收敛判断，其中，所述Grad_k+1根据步骤S24中的Grad_k计算得到；

S27、若k＝|x|，则令x₁＝x_k+1，返回步骤S22，否则，更新校正矩阵H，返回步骤S23，其中，校正矩阵H更新式为p＝x_k+1-x_k，q＝Grad_k+1-Grad_k；

S28、输出最优参数向量State_Final，释放CPU和GPU端存储空间；

S29、计算似然比检测值ΛH₁/H₀与航迹校验门限T_vali，验证新生轨迹是否存在，若存在则计算新生轨迹的协方差矩阵，将所述新生轨迹加入已存在轨迹，若不存在则返回步骤S10，具体如下：

S2901、将S28所述最优参数向量State_Final映射到量测空间，并进行量测检验，收集落入关联门限的量测，更新量测矩阵与量测数量记录向量；

S2902、计算对应于State_Final的JLLR值JLLR_State；

S2903、基于算得ΛH₁/H₀值，其中μ₁为航迹存在条件下单帧观测的JLLR期望值，σ₁为航迹存在条件下单帧观测的JLLR标准差，所述μ₁和σ₁在一个跟踪区域内只需计算一次即可，都可以通过在线计算或者都可以通过离线计算；

S2904、根据公式计算门限T_vali；

S2905、判断新生轨迹是否存在，若存在则计算新生轨迹的协方差矩阵，将所述新生轨迹加入已存在轨迹，若不存在则返回步骤S10。

2.根据权利要求1所述基于GPU架构的密集杂波条件下机动多目标跟踪算法，其特征在于：S1102所述求得每个遗传种子在对应帧数下的局部联合似然比之和时，需要对max(Obser_num)进行判定：若max(Obser_num)≥64，则通过多线程递归的方式进行求和；若max(Obser_num)＜64，则使用一个线程求和。

3.根据权利要求1所述基于GPU架构的密集杂波条件下机动多目标跟踪算法，其特征在于：S24所述循环计算的次数<9。

4.根据权利要求1所述基于GPU架构的密集杂波条件下机动多目标跟踪算法，其特征在于：S2903所述计算方式为离线计算。