具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

如前所述，请参阅图1，使用扇扫雷达进行跟踪，对于目标A，若天线在T₁时刻顺时针扫描到它，当前扫描周期的结束时刻为T_end1，跟踪延迟为delay_a＝T_end1-T₁；对于目标B，若天线在T₂时刻逆时针扫描到它，当前扫描周期的结束时刻为T_end2，跟踪延迟为delay_b＝T_end2-T₂。可以看出，在一个扫描周期中，目标越早被扫到，对该目标的跟踪延迟越大。此外，对于同一个目标，天线扫描的方向不同，跟踪延迟也不同。

想要最大程度降低跟踪延迟，最佳方案是在天线扫到航迹的波门末端时就对该航迹的状态进行更新。这种方案的具体实现方法依托于基于波门的观测选择技术。在观测选择技术中，波门是观测的一种范围，它以接近1的概率包含了来自当前航迹所跟踪的目标的观测。因此，当天线扫过波门时，只有极小的概率没有扫描到来自目标的观测。在这个时候更新就能实现理论上跟踪延迟的最小值。本发明对波门的形状进行研究，并利用几何关系求出天线扫过波门末端时的方位角，从而得到对应的时刻，由此，更新的时刻在波门建立时就可以确定。当扇扫雷达的时间达到波门末端对应时刻时，就对相应的航迹进行更新，从而实现最大程度地降低跟踪延迟。本发明提供的连续跟踪方法可简称为S-CIPDA(SectorCountinuous Intergrated Probabilistic Data Association)跟踪方法。

下面描述以上构思的具体实现方式。

扇扫雷达的扫描周期定义为雷达天线从监视区域一端扫描到另一端的时间。如图2(a)和图2(b)所示，在扫描周期开始时为加速阶段，天线进行匀加速运动，到达指定扫描角速度后进行匀速运动。在扫描周期快结束时为减速阶段，天线进行匀减速运动，速度变为零时，天线刚好扫描到监视区域的另一端。随后，天线又进行反向扫描。在每一个扫描周期中，雷达的监视区域分为加速区、匀速区和减速区。假设天线的扫描角度范围为[-β,β]，β为天线扫描的扇区的半宽度，天线匀速扫描时的扫描角速度为ω₀，天线加、减速阶段扫描时的加速度绝对值为α，天线在加、减速阶段扫过的方位角范围为天线单次扫描的扫描周期为T＝2ω₀/α+2(β-δ)/ω₀。天线的扫描方向用f表示，当f＝-1时，天线为顺时针扫描，当f＝1时，天线为逆时针扫描。

如图3所示，本发明实施例提供的一种扇扫雷达多目标连续跟踪方法包括：

步骤301、对扇扫雷达接收到的每一个观测z_re，取该观测z_re所在扫描周期的前一个扫描周期内获取的、未被已有航迹利用的观测，组成第一预选观测集设第一预选观测集内包括的观测数为k表示当前扇扫雷达接收观测的扫描周期数，k大于等于2；

根据第一预选观测集建立航迹，并结合航迹的状态进行滤波，计算下一个扫描周期(即第k+1个扫描周期)航迹的预测估计及对应的、用于下一次更新(即第k+1次更新)航迹的波门时间区间。

本发明提供的方法中，每接收到一个观测，就进行航迹初始化的操作，尝试建立航迹。为提高新航迹的可靠性，可采用速度波门作为建立新航迹的条件，若满足条件则建立新的航迹，若不满足则保存接收到的观测，以备后续使用。显然，若第一预选观测集为空集，认为不能建立航迹。

优选地，步骤301中，根据第一预选观测集建立航迹时，以第一预选观测集中的每个观测为中心建立速度波门，判断是否能建立新的航迹，即，对第一预选观测集中的每一个观测，均判断该观测是否满足如下初始化条件：

其中，z_x ^re和z_y ^re分别表示当前时刻扇扫雷达接收到的观测z_re的横坐标和纵坐标位置，z_x ⁱ和z_y ⁱ分别表示第一预选观测集内第i个观测的横坐标和纵坐标位置，T_r为表示当前时刻与扇扫雷达接收到第一预选观测集内第i个观测的时刻的时间差；v_max表示目标最大速度，v_error表示目标速度计算的误差，可根据实际情况设置；

若第一预选观测集内第i个观测满足上述初始化条件，则根据当前时刻扇扫雷达接收到的观测及第一预选观测集内第i个观测进行航迹初始化，建立航迹，得到航迹的初始状态。

建立新的航迹时，根据当前观测z_re的位置可确定新航迹的位置，再结合当前观测z_re与航迹中前一个观测之间的时间戳之差可计算出新航迹的速度，进而确定该航迹的初始状态。航迹初始状态的协方差为：

其中，R_c表示观测在直角坐标系下的观测噪声协方差矩阵。

步骤302、对各个已有的航迹，均采用如下方式在连续跟踪过程中迭代更新航迹的估计及对应的波门时间区间：

每个扫描周期内，在当前航迹对应的波门时间区间的结束时刻，根据波门时间区间确定采集时间区间，获取采集时间区间内扇扫雷达接收到的观测，组成第二预选观测集；

在采集时间区间的结束时刻，利用第二预选观测集对当前航迹的状态进行更新，并结合航迹的状态进行滤波，计算下一个扫描周期航迹的预测估计及对应的、用于下一次更新航迹的波门时间区间。

本发明的方法在连续跟踪过程中，利用第二预选观测集Z_τ(k)对当前航迹τ的状态进行更新，更新完成后对当前航迹τ下一个扫描周期内的状态进行预测，得到下一个扫描周期航迹τ的估计，以便实现滤波。通过本发明的方法，当扇扫雷达的时间达到波门末端对应时刻(即波门时间区间的结束时刻)就对相应的航迹进行更新，从而有效缩减延迟。

考虑到在连续跟踪的过程中可能有多个航迹，而各个航迹的波门可能有交叉，即，对应的波门时间区间有交叉，因此，为进一步优化航迹更新的时刻，步骤302中，根据波门时间区间确定采集时间区间，包括：

若当前航迹τ与其余已有的航迹对应的波门时间区间均没有交叉，则令采集时间区间与当前航迹对应的波门时间区间相同；

若当前航迹τ与其余已有的航迹对应的波门时间区间有交叉，且所有与当前航迹τ有交叉的航迹对应的最晚的波门时间区间结束时刻等于或早于当前航迹对应的波门时间区间结束时刻，则令所有与当前航迹τ有交叉的航迹及当前航迹对应的最早的波门时间区间起始时刻为采集时间区间起始时刻，当前航迹对应的波门时间区间结束时刻为采集时间区间结束时刻；

若当前航迹τ与其余已有的航迹对应的波门时间区间有交叉，且所有与当前航迹τ有交叉的航迹对应的最晚的波门时间区间结束时刻晚于当前航迹对应的波门时间区间结束时刻，则令所有与当前航迹τ有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，所有与当前航迹τ有交叉的航迹对应的最晚的波门时间区间结束时刻为采集时间区间结束时刻。其中，“早”指的是时间上靠前，相应的数值较小，“晚”指的是时间上靠后，相应的数值也更大。

进一步地，上述方式在实现时，可通过判断当前航迹τ是否被推迟更新进行区分，若航迹τ没有被标记为推迟更新的航迹，则按照航迹τ对应的波门时间区间结束时刻进行更新，若航迹τ被标记为推迟更新的航迹，则需要调整进行更新的时刻，具体地，步骤302包括：

步骤302-1、第k个扫描周期内(也即第k次更新时)，在当前航迹τ对应的波门时间区间[t_s(k),t_d(k)]的结束时刻t_d(k)，判断当前航迹τ是否被标记为推迟更新的航迹，是则跳转执行步骤302-7，否则执行步骤302-2；

步骤302-2、判断当前航迹τ是否和其余航迹有波门交叉，是则执行步骤302-4，否则执行步骤302-3；

步骤302-3、令采集时间区间与当前航迹对应的波门时间区间相同，即，取当前航迹τ的波门时间区间[t_s(k),t_d(k)]内的观测，组成第二预选观测集Z_τ(k)；

利用第二预选观测集Z_τ(k)对当前航迹τ进行更新和预测，计算当前航迹τ下一个扫描周期更新的波门时间区间[t_s(k+1),t_d(k+1)]，调整调整扫描方向参数f；其中扫描方向参数f用于调整扇扫雷达天线的扫描方向；

步骤302-4、若当前航迹τ与其余航迹波门交叉，但其余航迹最晚的波门结束时刻t_{d_max}不晚于航迹τ的波门结束时刻t_d(k)，即t_{d_max}≤t_d(k)，则执行步骤302-5；

若当前航迹τ与其余波门交叉，但其余航迹最晚的波门结束时刻t_{d_max}晚于航迹τ的波门结束时刻t_d(k)，即t_{d_max}＞t_d(k)，则执行步骤302-6；t_{d_max}为当前航迹以及与其有交叉的各航迹对应的最晚的波门时间区间结束时刻；

步骤302-5、令有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，当前航迹对应的波门时间区间结束时刻为采集时间区间结束时刻，即，取出时间区间[t_{s_min},t_d(k)]内的观测组成第二预选观测集Z_τ(k)，t_{s_min}为当前航迹以及与其有交叉的各航迹对应的最早的波门时间区间起始时刻；

利用第二预选观测集Z_τ(k)对当前航迹τ进行更新和预测，计算当前航迹τ下一个扫描周期更新的波门时间区间[t_s(k+1),t_d(k+1)]，调整调整扫描方向参数f；

步骤302-6、令有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，有交叉的航迹对应的最晚波门时间区间结束时刻为采集时间区间结束时刻，即，t_d(k)'＝t_{d_max}＝max(t_d1,...,t_dn,t_d(k))，其中t_d1,...,t_dn为n个与航迹τ交叉的波门结束时刻，t_s(k)'＝t_{s_min}＝min(t_s1,...,t_sn,t_s(k))，其中t_s1,...,t_sn为n个与航迹τ交叉的波门起始时刻，采集时间区间为[t_s(k)',t_d(k)']，将航迹τ标记为推迟更新的航迹，并执行步骤302-7；

步骤302-7、在采集时间区间结束时刻t_d(k)'，取出时间区间[t_s(k)',t_d(k)']内的观测，组成第二预选观测集Z_τ(k)；

利用第二预选观测集Z_τ(k)对当前航迹τ进行更新和预测，计算当前航迹τ下一个扫描周期更新的波门时间区间[t_s(k+1),t_d(k+1)]，调整调整扫描方向参数f。

这一个优选的实施方式中，步骤302可通过如下表1所示的伪代码实现多目标更新过程：

表1多目标更新过程伪代码

优选地，步骤302中，计算下一个扫描周期航迹的预测估计及对应的波门时间区间，包括：

步骤302-A、计算第二预选观测集Z_τ(k)中各观测对应的杂波密度。

进一步地，步骤302-A包括：

首先，寻找观测z_k(i)与第二预选观测集Z_τ(k)中其余观测之间的第n小的距离rⁿ(i)，n为大于0的整数，优选为n＝2，z_k(i)∈Z_τ(k)，i的取值为从1到第二预选观测集Z_τ(k)的观测数，k表示扫描周期数；若第二预选观测集Z_τ(k)中的观测数小于(n+1)，则令从当前扫描周期(即第k个扫描周期)的开始到当前时刻所接收到的观测合集为寻找观测z_k(i)在观测合集中第n小的距离rⁿ(i)以便计算杂波密度。

其次，计算观测z_k(i)对应的稀疏性，表达式为：

γ(z_k(i))＝V(rⁿ(i))/n；

其中，

Γ(·)为伽马函数，为空间的维数，优选为

最后，根据稀疏性计算杂波密度，表达式为：

ρ_k(i)表示第二预选观测集Z_τ(k)中观测z_k(i)对应的杂波密度，对第二预选观测集Z_τ(k)中每个观测均执行如上计算过程。

步骤302-B、根据波门对第二预选观测集Z_τ(k)中的观测进行选择，得到用于更新的观测集。

一个航迹可能对应多个运动模型，例如一个航迹可能对应匀速模型、加速模型、转弯模型等，航迹的每个模型有相应的模型概率，表示该航迹的符合该运动模型的概率。

考虑到航迹的模型可能并不唯一，优选地，步骤302-B中，第k次更新时，根据波门对第二预选观测集Z_τ(k)中的观测进行选择，包括：

对于当前航迹τ的模型σ，选择的观测集y_k(σ)满足：

其中，y为观测集y_k(σ)中的观测，表示当前航迹τ第k次更新时模型σ的预测观测，上标“T”表示转置，S_k(σ)当前航迹τ第k次更新时模型σ的新息矩阵，为S_k(σ)的逆矩阵，S_k(σ)属于第k-1次更新时航迹的估计，g为波门大小，表达式为：

chi2inv(·)为χ²分布函数的反函数，P_G表示波门概率，l表示观测向量y的维数，一般情况下，观测向量y的维数为2；

对于当前航迹τ不同的模型，将所有选择的观测集观测合并成一个用于更新的观测集y_k，表达式为：

步骤302-C、计算第二预选观测集Z_τ(k)中各个观测与当前航迹的似然函数。

优选地，第二预选观测集Z_τ(k)中的观测z_k(i)与当前航迹τ的模型σ的似然函数表达式为：

其中，z_k(i)∈V_k(σ)表示观测z_k(i)被当前航迹τ的模型σ选中，V_k(σ)表示当前航迹τ的模型σ的波门范围，i的取值为从1到第二预选观测集Z_τ(k)的观测数，N(·)为高斯分布的概率密度函数；

得到观测z_k(i)与当前航迹τ的似然函数表达式为：

其中，μ_k|k-1(σ)表示航迹τ第k次更新时模型σ的模型预测概率，表示观测z_k(i)与当前航迹τ的模型σ的似然函数。

步骤302-D、结合似然函数，对用于更新的观测集y_k中的观测y_k(i)对应的杂波密度进行调制。

优选地，步骤302-D包括：

结合似然函数，计算观测y_k(i)是航迹τ对应的目标在第k次更新的观测结果的先验概率P^τ(i)，表达式为：

其中，θ_k(i),i＞0表示观测y_k(i)是源自于目标检测的事件，θ_k(0)表示没有一个观测是目标检测的事件，Y^k-1表示从初始扫描周期到第k-1个扫描周期扇扫雷达所接收到的所有观测的合集,表示航迹τ第k次更新时目标存在的先验概率，表示航迹τ对应的目标存在的事件，P_D表示扇扫雷达对目标的发现概率，P_G表示波门概率，m_k表示用于更新的观测集y_k中的观测数，i的取值为从1到m_k；

得到对于航迹τ，观测y_k(i)对应的、调制后的杂波密度为：

其中，航迹η为除当前航迹航迹τ之外的其余已有的航迹。

若y_k(i)不被其他航迹选中，则它对应的杂波密度不变。

步骤302-E、将用于更新的观测集y_k与当前航迹τ进行数据关联，得到目标存在的后验概率和数据关联后验概率Y^k表示从初始扫描周期到第k个扫描周期扇扫雷达所接收到的所有观测的合集。

优选地，步骤302-E包括：

计算第k次更新，观测y_k(i)的似然比Λ_k，表达式为：

得到航迹τ第k次更新后目标存在的后验概率表达式为：

其中，表示航迹τ第k次更新时目标存在的先验概率；

得到数据关联后验概率表达式为：

步骤302-F、进行航迹状态更新，得到当前航迹τ对应的目标轨迹状态。

优选地，步骤302-F包括：

首先计算给定第i个观测y_k(i)的后验模型概率，表达式为：

其中，μ_k|k-1(σ)表示航迹τ第k次更新时模型σ的模型预测概率；

接着计算当前航迹τ各模型的模型概率，表达式为：

然后计算后验数据关联概率，表达式为：

其中，β_k(i)表示数据关联后验概率，表示观测z_k(i)与当前航迹τ的模型σ的似然函数，表示观测z_k(i)与当前航迹τ的似然函数。

计算给定观测y_k(i)时当前航迹τ每个模型的状态估计值和协方差，表达式为：

其中，表示给定观测y_k(i)时当前航迹τ模型σ的状态估计值，P_k|k(i，σ)表示给定观测y_k(i)时当前航迹τ模型σ的状态协方差，表示航迹τ第k次更新时模型σ的状态预测值，P_k|k-1(σ)表示航迹τ第k次更新时模型σ的预测协方差，表示卡尔曼滤波估计，H表示观测矩阵，R表示观测噪声协方差矩阵。

得到当前航迹τ对应的目标轨迹状态表达式为：

其中，表示航迹τ第k次更新后的航迹状态，P_k|k表示航迹τ第k次更新后的航迹状态协方差。

步骤302-G、计算下一个扫描周期航迹的预测估计及对应的波门时间区间。

优选地，步骤302-G中，计算航迹的每个模型下一次被更新时的状态预测、预测观测、目标存在的先验概率和航迹波门时间，具体地：

假设共有M个模型，第σ个模型对应的航迹状态均值为航迹状态协方差为P_k|k(σ)。

首先计算每个模型的模型预测概率，表达式为：

其中，Γ_σm表示目标在第k次更新遵循第σ个模型的条件下，在第k+1次更新切换到第m个模型的概率，μ_k+1|k(m)表示模型预测概率。

接下来求出混合模型概率和混合状态，表达式为：

其中，表示第k+1次更新是m模型的条件下第k次更新是σ模型的先验概率，简称混合模型概率，表示模型m的混合模型状态，表示模型m的混合模型状态协方差。

对于模型m，状态预测的表达式为：

其中，表示航迹τ第k+1次更新时模型m状态预测值，P_k+1|k(m)表示航迹τ第k+1次更新时模型m的预测协方差矩阵。F表示状态转移矩阵，F^T为F的转置，Q表示过程噪声矩阵。

得到预测观测表达式为：

S_k(m)＝HP_k+1|k(m)H^T+Q；

其中，表示当前航迹第k+1次更新时模型m的预测观测，S_k(m)表示当前航迹第k次更新时模型m的新息矩阵。

目标存在的先验概率表达式为：

其中，表示航迹τ第k+1次更新时目标存在的先验概率，γ为目标存在性的转换概率，可根据实际情况自行设置。

优选地，步骤302中，计算下一个扫描周期航迹对应的波门时间区间，包括：

根据计算下一个扫描周期航迹的预测估计得到的新息矩阵S_k(m)，对当前航迹的模型m计算波门方位角的范围[θ_min,θ_max]；

进一步地，波门方位角的范围的表达式为：

其中，g为波门大小，θ₀为波门中心方位角，即模型m相应的预测观测的方位角，S_ij，i，j∈{1，2}为模型m相应的新息矩阵S_k(m)的各个元素；对当前航迹的各模型均采用上述方式对应计算波门方位角的范围；

根据波门方位角的范围，计算当前航迹的各模型的波门时间范围，如模型m的波门时间范围为

将各模型对应的最早波门时间区间起始时刻作为当前航迹的波门时间区间起始时刻，各模型对应的最晚波门时间区间结束时刻作为当前航迹的波门时间区间结束时刻，即，对第k个扫描周期，航迹τ的波门时间区间起始时刻t_s(k)＝min{t_s ¹，...，t_s ^M}，波门时间区间结束时刻为t_d(k)＝max{t_d ¹，...，t_d ^M}，

步骤301中，根据第一预选观测集建立航迹，并结合航迹的状态进行滤波，计算下一个扫描周期航迹的预测估计及对应的波门时间区间，具体实现方式可参考上述步骤302，在此不再进一步赘述。

优选地，根据波门方位角的范围，计算当前航迹的各模型的波门时间范围时，设对于方位角θ，天线从扫描开始到扫描到方位角θ花费的时间为ΔT_s，当天线处于匀速扫描区域，|θ|＜β-δ，ΔT_s的表达式为：

其中，ω₀为天线匀速扫描时的扫描角速度，α为天线加、减速阶段的加速度绝对值，天线在加减速阶段扫过的方位角范围为天线的扫描范围为[-β，β]，天线的扫描方向用f表示，当f＝-1时，天线为顺时针扫描，当f＝1时，天线为逆时针扫描；

当天线处于加速阶段，θ＞β-δ且f＝-1，或者θ＜-β+δ且f＝1，ΔT_s的表达式为：

当天线处于减速阶段，θ＞β-δ且f＝1，或者θ＜-β+δ且f＝-1，ΔT_s的表达式为：

请参阅图4至图9(b)本发明还对提出的扇扫雷达多目标连续跟踪方法进行了仿真验证，雷达位于笛卡尔坐标系的原点(0m，0m)。有两个目标，均经历匀速-匀加速-匀速-转弯-匀速-转弯-匀速的运动过程，两个目标(即目标1、目标2)的初始状态如表2所示，x表示目标位置的横坐标，表示目标在横坐标上的速度，表示目标在横坐标上的加速度，y表示目标位置的纵坐标，表示目标在纵坐标上的速度，表示目标在纵坐标上的加速度。

表2两个目标的初始状态

两个目标加速阶段的加速度分别为[-0.2m/s²,0.1m/s²]和[0.2m/s²,0.1m/s²]。在转弯阶段，它们的角速度均为π/30，进行100次蒙特卡洛仿真，每次仿真对70次更新过程进行模拟。监视区域中的杂波密度为1.0×10^-7/scan/m²且杂波满足均匀分布。雷达天线扫描的扇区的半宽度为β＝π/3，匀速扫描时的扫描速度为(π/6)rad/s，角加速度为(π/3)rad/s²。根据这些数据可以进一步得到加、减速区域的方位角范围为δ＝π/24，单次扫描周期为T＝4.5s。目标运动的过程噪声设置为ds＝0.01m/s²，距离、方位角的观测噪声标准差分别为σ_r＝5m，σ_θ＝0.01deg。

得到两个目标的运动轨迹如图4所示，本发明利用现有技术中的IPDA方法(即传统的概率互联方法)进行对比，如图5(a)至图9(b)，从图5(a)和图5(b)可以看出IPDA方法和本发明的S-CIPDA方法都成功地跟踪上了目标，但根据图6(a)至图7(b)可知S-CIPDA方法的跟踪延迟远低于IPDA方法，IPDA方法对目标的跟踪延迟与目标位置和扫描方向有关，跟踪延迟也随着更新而不断起伏，且根据图8(a)至图9(b)可知，在大多数时候IPDA方法和本发明的S-CIPDA方法的跟踪准确度十分接近，有时IPDA方法的跟踪准确度相比S-CIPDA方法要低一点。本发明提出的扇扫雷达多目标连续跟踪方法大幅降低了跟踪延迟，并且跟踪延迟的降低并没有导致跟踪准确率下降。

特别地，在本发明一些优选的实施方式中，还提供了一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现上述任一实施方式中扇扫雷达多目标连续跟踪方法的步骤。

在本发明另一些优选的实施方式中，还提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述任一实施方式中扇扫雷达多目标连续跟踪方法的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述扇扫雷达多目标连续跟踪方法实施例的流程，在此不再重复说明。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。