CN113406618B - 一种tws雷达多目标连续跟踪方法 - Google Patents

一种tws雷达多目标连续跟踪方法 Download PDF

Info

Publication number
CN113406618B
CN113406618B CN202110691909.6A CN202110691909A CN113406618B CN 113406618 B CN113406618 B CN 113406618B CN 202110691909 A CN202110691909 A CN 202110691909A CN 113406618 B CN113406618 B CN 113406618B
Authority
CN
China
Prior art keywords
track
observation
flight path
time interval
current
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110691909.6A
Other languages
English (en)
Other versions
CN113406618A (zh
Inventor
周共健
桑航
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202110691909.6A priority Critical patent/CN113406618B/zh
Publication of CN113406618A publication Critical patent/CN113406618A/zh
Application granted granted Critical
Publication of CN113406618B publication Critical patent/CN113406618B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/66—Radar-tracking systems; Analogous systems
    • G01S13/72—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
    • G01S13/723—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar by using numerical data
    • G01S13/726—Multiple target tracking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/006—Theoretical aspects

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明涉及一种TWS雷达多目标连续跟踪方法，该方法包括对TWS雷达接收到的每一个观测，选取当前时刻之前获取的观测组成第一预选观测集；根据第一预选观测集建立航迹，并结合航迹的状态进行滤波，计算下一次更新航迹的预测估计及对应的波门时间区间；对各个已有的航迹，每次更新，在当前航迹对应的波门时间区间结束时刻，根据波门时间区间确定采集时间区间，获取采集时间区间内的观测，组成第二预选观测集；在采集时间区间结束时刻，利用第二预选观测集对当前航迹的状态进行更新，并结合航迹的状态进行滤波，计算下一次更新航迹的预测估计及对应的波门时间区间。本发明能够实现TWS雷达低延迟跟踪，并解决扫描边界引起的航迹断裂问题。

Description

一种TWS雷达多目标连续跟踪方法
技术领域
本发明涉及目标检测技术领域，尤其涉及一种TWS雷达多目标连续跟踪方法、计算机设备及计算机可读存储介质。
背景技术
大多数边扫描边跟踪雷达(TWS雷达)为机械扫描雷达，靠天线的机械旋转来改变天线波束的方向，从而对监视区域进行扫描。在扫描到目标后，天线将对目标的观测返回到雷达，雷达利用接收到的观测对目标进行跟踪。对于圆周扫描TWS雷达，当目标处于监视区域内的不同方位角的位置时，天线扫描到它的时间也不同，但是在同一扫描周期内，这些目标的状态更新却是同时进行，这导致从接收观测到处理观测之间有一定的延迟，即出现延迟问题。此外，在对目标进行跟踪的时候，传统的跟踪方法都假设在一个扫描周期内天线会探测到有且只有一个源于目标的观测，在一个扫描周期内只更新一次航迹状态。但是在实际运用中，这种假设并不是总成立的。对于圆周扫描TWS雷达，当目标处在扫描边界(即扫描开始或结束的位置)附近时，其在一个扫描周期内被探测到的次数可能不等于1，即出现边界问题，进而影响到跟踪的准确性、可靠性。
发明内容
(一)要解决的技术问题
本发明要解决的技术问题是解决现有技术无法实现TWS雷达低延迟跟踪、无法保证目标越过扫描边界时的跟踪稳定性的问题。
(二)技术方案
为了解决上述技术问题，本发明提供了一种TWS雷达多目标连续跟踪方法，包括：
步骤(1)、对TWS雷达接收到的每一个观测，选取当前时刻之前获取的观测组成第一预选观测集，所述第一预选观测集内的观测的接收时刻与当前时刻间隔0.8T～1.2T且未被已有的航迹利用，T为TWS雷达的扫描周期；根据所述第一预选观测集建立航迹，并结合航迹的状态进行滤波，计算下一次更新航迹的预测估计及对应的波门时间区间；
步骤(2)、对各个已有的航迹，均采用如下方式在连续跟踪过程中迭代更新航迹的估计及对应的波门时间区间：
每次更新，在当前航迹对应的波门时间区间结束时刻，根据波门时间区间确定采集时间区间，获取采集时间区间内的观测，组成第二预选观测集；在所述采集时间区间结束时刻，利用所述第二预选观测集对当前航迹的状态进行更新，并结合航迹的状态进行滤波，计算下一次更新航迹的预测估计及对应的波门时间区间。
优选地，步骤(1)中，根据所述第一预选观测集建立航迹时，对所述第一预选观测集中的每一个观测，均判断是否满足如下初始化条件：
Figure GDA0003812235820000021
其中，zx re和zy re分别表示当前时刻接收到的观测的横坐标和纵坐标位置，zx i和zy i分别表示所述第一预选观测集内第i个观测的横坐标和纵坐标位置，Tr为表示当前时刻与接收到所述第一预选观测集内第i个观测的时刻的时间差；vmax表示目标最大速度，verror表示目标速度计算的误差；
若所述第一预选观测集内第i个观测满足初始化条件，则根据当前时刻的观测及第一预选观测集内第i个观测进行航迹初始化，建立航迹，得到航迹的初始状态。
优选地，步骤(2)中，根据波门时间区间确定采集时间区间，包括：
若当前航迹与其余已有的航迹对应的波门时间区间均没有交叉，则令采集时间区间与当前航迹对应的波门时间区间相同；
若当前航迹与其余已有的航迹对应的波门时间区间有交叉，且有交叉的航迹对应的最晚波门时间区间结束时刻等于或早于当前航迹对应的波门时间区间结束时刻，则令有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，当前航迹对应的波门时间区间结束时刻为采集时间区间结束时刻；
若当前航迹与其余已有的航迹对应的波门时间区间有交叉，且有交叉的航迹对应的最晚波门时间区间结束时刻晚于当前航迹对应的波门时间区间结束时刻，则令有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，有交叉的航迹对应的最晚波门时间区间结束时刻为采集时间区间结束时刻。
优选地，步骤(2)中，计算下一次更新航迹的预测估计及对应的波门时间区间，包括：
计算所述第二预选观测集中各观测对应的杂波密度；
根据波门对所述第二预选观测集中的观测进行选择，得到用于更新的观测集；
计算所述第二预选观测集中各个观测与当前航迹的似然函数；
结合似然函数，对用于更新的观测集中的观测对应的杂波密度进行调制；
将用于更新的观测集与当前航迹进行数据关联，得到目标存在的后验概率和数据关联后验概率；
进行航迹分量控制；
进行航迹状态更新，得到当前航迹对应的目标轨迹状态；
计算每个航迹分量下一次更新航迹的预测估计及对应的波门时间区间。
优选地，步骤(2)中，计算所述第二预选观测集中各观测对应的杂波密度，包括：
寻找观测zk(i)与所述第二预选观测集Zτ(k)中其余观测之间的第n小的距离rn(i)，n为大于0的整数，zk(i)∈Zτ(k)，k表示更新次数；若所述第二预选观测集Zτ(k)中的观测数小于(n+1)，则令从第k次更新扫描周期的开始到当前时刻所接收到的观测合集为
Figure GDA0003812235820000041
寻找观测zk(i)在观测合集
Figure GDA0003812235820000042
中第n小的距离rn(i)；
计算观测zk(i)对应的稀疏性，表达式为：
γ(zk(i))＝V(rn(i))/n；
其中，
Figure GDA0003812235820000043
Γ(·)为伽马函数，l为空间的维数；
根据稀疏性计算杂波密度，表达式为：
Figure GDA0003812235820000044
优选地，步骤(2)中，第k次更新时，根据波门对所述第二预选观测集Zτ(k)中的观测进行选择，包括：
对于当前航迹τ的航迹分量ck-1的模型σ，选择的观测集yk(ck-1,σ)满足：
Figure GDA0003812235820000045
其中，y为观测集yk(ck-1,σ)中的观测，
Figure GDA0003812235820000046
表示当前航迹τ第k次更新时航迹分量ck-1的模型σ的预测观测，Sk(ck-1,σ)表示当前航迹τ第k次更新时航迹分量ck-1的模型σ的新息矩阵，
Figure GDA0003812235820000047
为Sk(ck-1,σ)的逆矩阵，g为波门大小，表达式为：
Figure GDA0003812235820000048
chi2inv(·)为χ2分布函数的反函数，PG表示波门概率，l表示观测y的维数；
对于当前航迹τ不同的模型，将所有选择的观测集观测合并成用于更新的观测集
Figure GDA0003812235820000051
优选地，步骤(2)中，进行航迹分量控制，包括：
将用于更新的观测集yk中的每个观测yk(i)与每个航迹分量ck-1关联，组成暂时航迹分量；暂时航迹分量的个数为
Figure GDA0003812235820000052
其中Ck-1表示第k-1次更新后航迹分量的个数，mk表示用于更新的观测集yk中的观测个数；
计算每个暂时航迹分量的相对概率，表达式为：
Figure GDA0003812235820000053
其中，
Figure GDA0003812235820000054
表示暂时航迹分量
Figure GDA0003812235820000055
的相对概率，βk(i)表示用来建立暂时航迹分量
Figure GDA0003812235820000056
的观测yk(i)与航迹τ的数据关联后验概率，p(ck-1)表示航迹τ第k-1次更新后航迹分量ck-1的相对概率，
Figure GDA0003812235820000057
表示观测zk(i)与当前航迹τ的航迹分量ck-1的似然函数，
Figure GDA0003812235820000058
表示观测zk(i)与当前航迹τ的似然函数；
进行航迹分量树叶裁剪，若一个暂时航迹分量
Figure GDA0003812235820000059
的相对概率
Figure GDA00038122358200000510
小于门限值，则去掉暂时航迹分量
Figure GDA00038122358200000511
进行航迹分量子树裁剪，从当前第k次更新向前追溯以前更新的航迹分量，对于第κ次更新，保留相对概率p(cκ)最大的航迹分量；其中κ＝k-N，N为子树裁剪深度；
将航迹分量树叶裁剪和航迹分量子树裁剪后，有一个暂时航迹分量子集Θ得以保留，而Θ的补集
Figure GDA00038122358200000512
则被删除；
计算补集
Figure GDA00038122358200000513
的相对概率Δp，表达式为：
Figure GDA00038122358200000514
Θ中的暂时航迹分量为确认航迹分量，用ck表示，航迹分量ck相对概率的表达式为：
Figure GDA00038122358200000515
调整目标存在的后验概率为：
p(χk|Yk)＝p(χk|Yk)(1-Δp)。
优选地，步骤(2)中，计算每个航迹分量下一次更新航迹对应的波门时间区间，包括：
根据计算航迹分量ck下一次更新航迹的预测估计得到的新息矩阵，对当前航迹τ的航迹分量ck的模型σ，计算相应的波门方位角的范围[θmin,θmax]；
波门方位角的范围的表达式为：
Figure GDA0003812235820000061
Figure GDA0003812235820000062
其中，g为波门大小，θ0为波门中心方位角，Sij,i,j∈{1,2}为航迹分量ck的模型σ的新息矩阵Sk(ck,σ)的各个元素；
根据波门方位角的范围，计算当前航迹τ的航迹分量ck的模型σ对应的波门时间范围；
计算出当前航迹τ的每个航迹分量的每个模型对应的波门时间范围后，将所有航迹分量的所有模型对应的最早波门时间区间起始时刻作为当前航迹的波门时间区间起始时刻，所有航迹分量的所有模型对应的最晚波门时间区间结束时刻作为当前航迹的波门时间区间结束时刻。
本发明还提供了一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现上述任一项所述TWS雷达多目标连续跟踪方法的步骤。
本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项所述TWS雷达多目标连续跟踪方法的步骤。
(三)有益效果
本发明的上述技术方案具有如下优点：本发明提供了一种TWS雷达多目标连续跟踪方法、计算机设备及计算机可读存储介质，本发明在多目标连续跟踪过程中，通过目标的航迹建立对应的波门时间区间，根据波门时间区间及时对航迹的状态进行更新，从而最大程度地降低跟踪延迟，实现TWS雷达低延迟跟踪，提高跟踪的准确性，并且本发明改变了原来按扫描周期更新的思路，改为按时间更新，每探测一次就更新一次，从而解决一个扫描周期内对目标的探测次数和更新次数不匹配造成的航迹断裂的问题，提高航迹维持能力。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是TWS雷达跟踪延迟示意图；
图2(a)是TWS雷达一种边界问题示意图；
图2(b)是TWS雷达另一种边界问题示意图；
图3是本发明实施例中一种TWS雷达多目标连续跟踪方法步骤示意图；
图4是目标运动轨迹图；
图5(a)是IMM-LMITS方法对图4的目标跟踪结果图；
图5(b)是本发明实施例中一种TWS雷达多目标连续跟踪方法对图4的目标跟踪结果图；
图6(a)示出了IMM-LMITS方法对目标1的跟踪延迟；
图6(b)示出了本发明实施例中一种TWS雷达多目标连续跟踪方法对目标1的跟踪延迟；
图7(a)示出了IMM-LMITS方法对目标2的跟踪延迟；
图7(b)示出了本发明实施例中一种TWS雷达多目标连续跟踪方法对目标2的跟踪延迟；
图8(a)示出了IMM-LMITS方法和本发明实施例中一种TWS雷达多目标连续跟踪方法跟踪目标1的位置RMSE；
图8(b)示出了IMM-LMITS方法和本发明实施例中一种TWS雷达多目标连续跟踪方法跟踪目标2的位置RMSE；
图9(a)示出了IMM-LMITS方法和本发明实施例中一种TWS雷达多目标连续跟踪方法跟踪目标1的速度RMSE；
图9(b)示出了IMM-LMITS方法和本发明实施例中一种TWS雷达多目标连续跟踪方法跟踪目标2的速度RMSE；
图10示出了IMM-LMITS方法和本发明实施例中一种TWS雷达多目标连续跟踪方法的真实航迹率比较。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
如前所述，请参阅图1，使用TWS雷达进行跟踪，设雷达天线按照逆时针旋转，空间中有两个目标，跟踪它们的航迹为T1(k1)和T2(k2)，若两航迹分别在t1和t2时刻被天线扫描到。随后，天线继续扫描，在tend时刻，天线结束这一扫描周期的扫描。目标的更新在tend时刻进行，则两航迹T1(k1)，T2(k2)的跟踪延迟分别为tdelay1＝tend-t1,tdelay2＝tend-t2。目标越接近扫描周期开始的位置，即目标位置对应的天线方位角越小，跟踪延迟越大。
此外，传统的跟踪方法在目标跨过边界时，也会产生边界问题。如图2(a)所示，天线按照逆时针旋转，若航迹T1(k1)所跟踪的目标第k1次扫描发生在t1时刻，第k1+1次扫描发生在t2时刻，而目标在t1和t2时刻之间顺时针穿过了边界，t2时刻早于该扫描周期的结束时刻tend，如此一来，在同一扫描周期，目标被天线探测到两次，返回了两个观测值。如图2(b)所示，天线按照逆时针旋转，在第k-1个扫描周期，若天线在t1时刻扫描到航迹状态T1(k1)所跟踪的目标，此后天线结束第k-1个扫描周期，进入第k个扫描周期。在此期间，目标逆时针穿过扫描边界，导致天线在第k个扫描周期没有扫描到目标。在第k+1个扫描周期t2时刻天线才重新扫描到目标，此时相应的航迹状态为T1(k1+1)。如此一来，天线在第k个扫描周期返回零个观测，即在第k个扫描周期天线是扫描不到目标的。对于传统方法来说，第k-1个扫描周期对第k个扫描周期目标观测zk的预测观测实际上出现在第k+1个扫描周期。这样一来，跟踪方法的预测观测与实际接收到的观测就出现了错位的现象。在第k个扫描周期结束时，传统方法利用预测观测对第k个扫描周期收到的所有观测进行观测选择以选出来自目标的观测。但实际上它所期望的来自目标的观测zk存在于第k+1个扫描周期的观测集。这样一来传统方法会认为第k个扫描周期出现了漏警，它在缺少观测的前提下对航迹状态进行更新估计，并再次预测下一次观测。在第k+1扫描周期结束时，传统方法继续进行更新。这一次传统方法所期待的观测zk+1实际上在第k+2个扫描周期时才被扫描到。如此一来，在边界问题发生后，传统方法一直都无法正确得到来自目标的观测。这个时候可能有两种后果：(1)航迹波门足够大，在第k+1扫描周期关联到了目标观测zk。这是一种错误关联，会导致跟踪系统的准确度降低，还可能导致航迹断裂的发生。(2)航迹波门较小，一直无法关联到来自目标的观测。这会导致航迹断裂的发生。但是这一航迹断裂的现象不会持续太久。几个扫描周期后一个新的航迹将被起始来跟踪当前目标。
想要最大程度降低跟踪延迟，最佳方案是在天线扫到航迹的波门末端时就对该航迹的状态进行更新。这种方案的具体实现方法依托于基于波门的观测选择技术。在观测选择技术中，波门是观测的一种范围，它以接近1的概率包含了来自当前航迹所跟踪的目标的观测。因此，当天线扫过波门时，只有极小的概率没有扫描到来自目标的观测。在这个时候更新就能实现理论上跟踪延迟的最小值。本发明对波门的形状进行研究，并利用几何关系求出天线扫过波门末端时的方位角，从而得到对应的时刻，由此，更新的时刻在波门建立时就可以确定。当TWS雷达的时间达到波门末端对应时刻时，就对相应的航迹进行更新，从而实现最大程度地降低跟踪延迟，解决延迟问题。对于边界问题，其根源在于在一个扫描周期内对目标的探测次数和更新次数不匹配。因此，本发明取缔了原来按扫描周期更新的思路，改为按时间更新，每探测一次就更新一次，这样一来也彻底解决了边界问题。本发明结合线性多目标方法(LM)和交互作用多模型方法(IMM)实现对多个目标的跟踪，本发明提供的连续跟踪方法可简称为IMM-LMCITS方法。
下面描述以上构思的具体实现方式。
如图3所示，本发明实施例提供的一种TWS雷达多目标连续跟踪方法包括：
步骤301、对TWS雷达接收到的每一个观测zre，选取获取观测zre的当前时刻之前获取的观测，组成第一预选观测集
Figure GDA0003812235820000101
第一预选观测集
Figure GDA0003812235820000102
内的观测相应的TWS雷达接收时刻与当前时刻间隔0.8T～1.2T，第一预选观测集
Figure GDA0003812235820000103
内的观测且未被已有的航迹利用，T为TWS雷达旋转一周的扫描周期，设第一预选观测集
Figure GDA0003812235820000104
包括的观测数为
Figure GDA0003812235820000105
k表示更新次数，k大于等于2，k＝1对应第一个扫描周期；
根据第一预选观测集
Figure GDA0003812235820000106
建立航迹，并结合航迹的状态进行滤波，计算下一次更新(即第k+1次更新)航迹的预测估计及对应的波门时间区间。
本发明提供的方法中，每接收到一个观测，就进行航迹初始化的操作，尝试建立航迹。为提高新航迹的可靠性，可采用速度波门作为建立新航迹的条件，若满足条件则建立新的航迹，若不满足则保存接收到的观测，以备后续使用。显然，若第一预选观测集
Figure GDA0003812235820000111
为空集，认为不能建立航迹。
优选地，步骤301中，根据第一预选观测集建立航迹时，以第一预选观测集中的每个观测为中心建立速度波门，判断是否能建立新的航迹，即，对第一预选观测集中的每一个观测，均判断该观测是否满足如下初始化条件：
Figure GDA0003812235820000112
其中，zx re和zy re分别表示当前时刻TWS雷达接收到的观测zre的横坐标和纵坐标位置，zx i和zy i分别表示第一预选观测集
Figure GDA0003812235820000113
内第i个观测的横坐标和纵坐标位置，
Figure GDA0003812235820000114
Tr为表示当前时刻与TWS雷达接收到第一预选观测集
Figure GDA0003812235820000115
内第i个观测的时刻的时间差；vmax表示目标最大速度，verror表示目标速度计算的误差，可根据实际情况设置；
若第一预选观测集
Figure GDA0003812235820000116
内第i个观测满足上述初始化条件，则根据当前时刻TWS雷达接收到的观测及第一预选观测集
Figure GDA0003812235820000117
内第i个观测进行航迹初始化，建立航迹，得到航迹的初始状态。
建立新的航迹时，根据当前观测zre的位置可确定新航迹的位置，再结合当前观测zre与航迹中前一个观测之间的时间戳之差可计算出新航迹的速度，进而确定该航迹的初始状态。航迹初始状态的协方差为：
Figure GDA0003812235820000118
其中，Rc表示观测在直角坐标系下的观测噪声协方差矩阵。
步骤302、对各个已有的航迹，均采用如下方式在连续跟踪过程中迭代更新航迹的估计及对应的波门时间区间：
每次更新，在当前航迹对应的波门时间区间的结束时刻，根据波门时间区间确定采集时间区间，获取采集时间区间内TWS雷达接收到的观测，组成第二预选观测集；
在采集时间区间的结束时刻，利用第二预选观测集对当前航迹的状态进行更新，并结合航迹的状态进行滤波，计算下一次更新航迹的预测估计及对应的波门时间区间。
本发明的方法在连续跟踪过程中，利用第二预选观测集Zτ(k)对当前航迹τ的状态进行更新，更新完成后对当前航迹τ下一次更新的状态进行预测，得到下一次更新航迹τ的估计，以便实现滤波。通过本发明的方法，当TWS雷达的时间达到波门末端对应时刻(即波门时间区间的结束时刻)就对相应的航迹进行更新，从而有效缩减延迟。
考虑到在连续跟踪的过程中可能有多个航迹，而各个航迹的波门可能有交叉，即，对应的波门时间区间有交叉，因此，为进一步优化航迹更新的时刻，步骤302中，根据波门时间区间确定采集时间区间，包括：
若当前航迹τ与其余已有的航迹对应的波门时间区间均没有交叉，则令采集时间区间与当前航迹对应的波门时间区间相同；
若当前航迹τ与其余已有的航迹对应的波门时间区间有交叉，且所有与当前航迹τ有交叉的航迹对应的最晚的波门时间区间结束时刻等于或早于当前航迹对应的波门时间区间结束时刻，则令所有与当前航迹τ有交叉的航迹及当前航迹对应的最早的波门时间区间起始时刻为采集时间区间起始时刻，当前航迹对应的波门时间区间结束时刻为采集时间区间结束时刻；
若当前航迹τ与其余已有的航迹对应的波门时间区间有交叉，且所有与当前航迹τ有交叉的航迹对应的最晚的波门时间区间结束时刻晚于当前航迹对应的波门时间区间结束时刻，则令所有与当前航迹τ有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，所有与当前航迹τ有交叉的航迹对应的最晚的波门时间区间结束时刻为采集时间区间结束时刻。其中，“早”指的是时间上靠前，相应的数值较小，“晚”指的是时间上靠后，相应的数值也更大。
进一步地，上述方式在实现时，可通过判断当前航迹τ是否被推迟更新进行区分，若航迹τ没有被标记为推迟更新的航迹，则按照航迹τ对应的波门时间区间结束时刻进行更新，若航迹τ被标记为推迟更新的航迹，则需要调整进行更新的时刻，具体地，步骤302包括：
步骤302-1、第k次更新时，在当前航迹τ对应的波门时间区间[ts(k),td(k)]的结束时刻td(k)，判断当前航迹τ是否被标记为推迟更新的航迹，是则跳转执行步骤302-7，否则执行步骤302-2；
步骤302-2、判断当前航迹τ是否和其余航迹有波门交叉，是则执行步骤302-4，否则执行步骤302-3；
步骤302-3、令采集时间区间与当前航迹对应的波门时间区间相同，即，取当前航迹τ的波门时间区间[ts(k),td(k)]内的观测，组成第二预选观测集Zτ(k)；
利用第二预选观测集Zτ(k)对当前航迹τ进行更新和预测，计算当前航迹τ下一次更新的波门时间区间[ts(k+1),td(k+1)]；
步骤302-4、若当前航迹τ与其余航迹波门交叉，但其余航迹最晚的波门结束时刻td_max不晚于航迹τ的波门结束时刻td(k)，即td_max≤td(k)，则执行步骤302-5；
若当前航迹τ与其余波门交叉，但其余航迹最晚的波门结束时刻td_max晚于航迹τ的波门结束时刻td(k)，即td_max＞td(k)，则执行步骤302-6；td_max为当前航迹以及与其有交叉的各航迹对应的最晚的波门时间区间结束时刻；
步骤302-5、令有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，当前航迹对应的波门时间区间结束时刻为采集时间区间结束时刻，即，取出时间区间[ts_min,td(k)]内的观测组成第二预选观测集Zτ(k)，ts_min为当前航迹以及与其有交叉的各航迹对应的最早的波门时间区间起始时刻；
利用第二预选观测集Zτ(k)对当前航迹τ进行更新和预测，计算当前航迹τ下一次更新的波门时间区间[ts(k+1),td(k+1)]；
步骤302-6、令有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，有交叉的航迹对应的最晚波门时间区间结束时刻为采集时间区间结束时刻，即，td(k)'＝td_max＝max(td1,...,tdn,td(k))，其中td1,...,tdn为n个与航迹τ交叉的波门结束时刻，ts(k)'＝ts_min＝min(ts1,...,tsn,ts(k))，其中ts1,...,tsn为n个与航迹τ交叉的波门起始时刻，采集时间区间为[ts(k)',td(k)']，将航迹τ标记为推迟更新的航迹，并执行步骤302-7；
步骤302-7、在采集时间区间结束时刻td(k)'，取出时间区间[ts(k)',td(k)']内的观测，组成第二预选观测集Zτ(k)；
利用第二预选观测集Zτ(k)对当前航迹τ进行更新和预测，计算当前航迹τ下一次更新的波门时间区间[ts(k+1),td(k+1)]。
这一个优选的实施方式中，步骤302可通过如下表1所示的伪代码实现多目标更新过程：
表1多目标更新过程伪代码
Figure GDA0003812235820000141
Figure GDA0003812235820000151
优选地，步骤302中，计算下一次更新航迹的预测估计及对应的波门时间区间，包括：
步骤302-A、计算第二预选观测集Zτ(k)中各观测对应的杂波密度。
进一步地，步骤302-A包括：
首先，寻找观测zk(i)与第二预选观测集Zτ(k)中其余观测之间的第n小的距离rn(i)，n为大于0的整数，优选为n＝2，zk(i)∈Zτ(k)，i的取值为从1到第二预选观测集Zτ(k)的观测数，k表示更新次数；若第二预选观测集Zτ(k)中的观测数小于(n+1)，则令从第k次更新相应扫描周期的开始到当前时刻所接收到的观测合集为
Figure GDA0003812235820000152
寻找观测zk(i)在观测合集
Figure GDA0003812235820000153
中第n小的距离rn(i)以便计算杂波密度。
其次，计算观测zk(i)对应的稀疏性，表达式为：
γ(zk(i))＝V(rn(i))/n；
其中，
Figure GDA0003812235820000154
Γ(·)为伽马函数，l为空间的维数，优选为l＝2；
最后，根据稀疏性计算杂波密度，表达式为：
Figure GDA0003812235820000155
ρk(i)表示第二预选观测集Zτ(k)中观测zk(i)对应的杂波密度，对第二预选观测集Zτ(k)中每个观测均执行如上计算过程。
步骤302-B、根据波门对第二预选观测集Zτ(k)中的观测进行选择，得到用于更新的观测集yk。
一个航迹可能对应多个航迹分量，每个航迹分量可能对应多个模型，优选地，步骤302-B中，第k次更新时，根据波门对第二预选观测集Zτ(k)中的观测进行选择，包括：
对于当前航迹τ的航迹分量ck-1的模型σ，选择的观测集yk(ck-1,σ)满足：
Figure GDA0003812235820000161
其中，y为观测集yk(ck-1,σ)中的观测，
Figure GDA0003812235820000162
表示当前航迹τ第k次更新时航迹分量ck-1的模型σ的预测观测，Sk(ck-1,σ)表示当前航迹τ第k次更新时航迹分量ck-1的模型σ的新息矩阵，
Figure GDA0003812235820000163
为Sk(ck-1,σ)的逆矩阵，上标“T”表示转置，g为波门大小，表达式为：
Figure GDA0003812235820000164
chi2inv(·)为χ2分布函数的反函数，PG表示波门概率，l表示观测向量y的维数，一般情况下，观测向量y的维数为2；
对于当前航迹τ不同的模型，将所有选择的观测集观测合并成一个用于更新的观测集yk，表达式为：
Figure GDA0003812235820000165
步骤302-C、计算第二预选观测集Zτ(k)中各个观测与当前航迹的似然函数。
优选地，第二预选观测集Zτ(k)中的观测zk(i)与当前航迹τ的航迹分量ck-1的模型σ的似然函数表达式为：
Figure GDA0003812235820000166
其中，zk(i)∈Vk(ck-1,σ)表示观测zk(i)被当前航迹τ的航迹分量ck-1的模型σ选中，Vk(ck-1,σ)表示当前航迹τ的航迹分量ck-1的模型σ的波门范围，i的取值为从1到第二预选观测集Zτ(k)的观测数，N(·)为高斯分布的概率密度函数；
得到观测zk(i)与当前航迹τ的航迹分量ck-1的似然函数表达式为：
Figure GDA0003812235820000171
得到观测zk(i)与当前航迹τ的似然函数表达式为：
Figure GDA0003812235820000172
其中，μk|k-1(ck-1,σ)表示航迹τ第k次更新时航迹分量ck-1的模型σ的模型预测概率，p(ck-1)表示航迹τ第k-1次更新后航迹分量ck-1的相对概率。
步骤302-D、结合似然函数，对用于更新的观测集yk中的观测yk(i)对应的杂波密度进行调制。
优选地，步骤302-D中，第k次更新时，包括：
结合似然函数，计算观测yk(i)是航迹τ对应的目标在第k次更新的观测结果的先验概率Pτ(i)，表达式为：
Figure GDA0003812235820000173
其中，θk(i),i＞0表示观测yk(i)是源自于目标检测的事件，θk(0)表示没有一个观测是目标检测的事件，Yk-1表示从初始到第k-1次更新TWS雷达所接收到的所有观测的合集,
Figure GDA0003812235820000174
表示航迹τ第k次更新时目标存在的先验概率，
Figure GDA0003812235820000175
表示航迹τ对应的目标存在的事件，PD表示TWS雷达对目标的发现概率，PG表示波门概率，mk表示用于更新的观测集yk中的观测数，i的取值为从1到mk；
得到对于航迹τ，观测yk(i)对应的、调制后的杂波密度为：
Figure GDA0003812235820000181
其中，航迹η为除当前航迹航迹τ之外的其余已有的航迹。
若yk(i)不被其他航迹选中，则它对应的杂波密度不变。
步骤302-E、将用于更新的观测集yk与当前航迹τ进行数据关联，得到目标存在的后验概率
Figure GDA0003812235820000182
和数据关联后验概率
Figure GDA0003812235820000183
Yk表示从初始时刻到第k次更新TWS雷达所接收到的所有观测的合集。
优选地，步骤302-E包括：
计算第k次更新，观测yk(i)的似然比Λk，表达式为：
Figure GDA0003812235820000184
得到航迹τ第k次更新后目标存在的后验概率
Figure GDA0003812235820000185
表达式为：
Figure GDA0003812235820000186
其中，
Figure GDA0003812235820000187
表示航迹τ第k次更新时目标存在的先验概率；
得到观测yk(i)与航迹τ的数据关联后验概率表达式为：
Figure GDA0003812235820000188
步骤302-F、进行航迹分量控制。
优选地，步骤302-F包括：
将用于更新的观测集yk中的每个观测yk(i)与每个航迹分量ck-1关联，组成暂时航迹分量；暂时航迹分量的个数为
Figure GDA0003812235820000189
其中Ck-1表示第k-1次更新后航迹分量的个数，mk表示用于更新的观测集yk中的观测个数；
计算每个暂时航迹分量的相对概率，表达式为：
Figure GDA0003812235820000191
其中，
Figure GDA0003812235820000192
表示暂时航迹分量
Figure GDA0003812235820000193
的相对概率，βk(i)表示用来建立暂时航迹分量
Figure GDA0003812235820000194
的观测yk(i)与航迹τ的数据关联后验概率，p(ck-1)表示航迹τ第k-1次更新后航迹分量ck-1的相对概率，
Figure GDA0003812235820000195
表示观测zk(i)与当前航迹τ的航迹分量ck-1的似然函数，
Figure GDA0003812235820000196
表示观测zk(i)与当前航迹τ的似然函数；
进行航迹分量树叶裁剪，若一个暂时航迹分量
Figure GDA0003812235820000197
的相对概率
Figure GDA0003812235820000198
小于设置的门限值，则去掉该暂时航迹分量
Figure GDA0003812235820000199
进行航迹分量子树裁剪，从当前第k次更新向前追溯以前更新的航迹分量，对于第κ次更新，将每一个航迹分量cκ进行比较，保留相对概率p(cκ)最大的航迹分量，κ＝k-N，N为子树裁剪深度，一般取N＝3；即，若航迹分量cκ在κ层对应的相对概率p(cκ)最大，则保留；记该保留下来的航迹分量为
Figure GDA00038122358200001910
由此，在第k次更新所有不是由
Figure GDA00038122358200001911
衍生的航迹分量都将被删除，将航迹分量树叶裁剪和航迹分量子树裁剪后，有一个暂时航迹分量子集Θ得以保留，而暂时航迹分量子集Θ的补集
Figure GDA00038122358200001915
则被删除；
计算补集
Figure GDA00038122358200001912
的相对概率Δp，表达式为：
Figure GDA00038122358200001913
子集Θ中的暂时航迹分量为确认航迹分量，用ck表示，航迹分量ck相对概率的表达式为：
Figure GDA00038122358200001914
调整目标存在的后验概率为：
p(χk|Yk)＝p(χk|Yk)(1-Δp)；
上式等式左侧为调整后的目标存在的后验概率，等式右侧为未调整的目标存在的后验概率。
步骤302-G、进行航迹状态更新，得到当前航迹τ对应的目标轨迹状态。
优选地，步骤302-G包括：
计算当前航迹τ每个航迹分量的状态估计值和协方差，表达式为：
Figure GDA0003812235820000201
其中，ck-1是确认航迹分量ck对应的上一次更新的航迹分量，
Figure GDA0003812235820000202
表示当前航迹τ的航迹分量ck的模型σ的状态估计值，Pk|k(ck,σ)表示当前航迹τ的航迹分量ck的模型σ的状态协方差，
Figure GDA0003812235820000203
表示航迹τ第k次更新时航迹分量ck-1的模型σ的状态预测值，Pk|k-1(ck-1,σ)表示航迹τ第k次更新时航迹分量ck-1的模型σ的预测协方差，
Figure GDA0003812235820000204
表示卡尔曼滤波估计，H表示观测矩阵，R表示观测噪声协方差矩阵。
各模型的后验模型概率为：
Figure GDA0003812235820000205
其中，μk|k-1(ck-1,σ)表示航迹τ第k次更新时航迹分量ck-1的模型σ的模型预测概率；
当前航迹τ对应的目标轨迹状态由各航迹分量组合得到，表达式为：
Figure GDA0003812235820000206
Figure GDA0003812235820000207
其中，
Figure GDA0003812235820000208
表示航迹τ第k次更新后的航迹状态，Pk|k表示航迹τ第k次更新后的航迹状态协方差。
步骤302-H、计算每个航迹分量下一次更新航迹的预测估计及对应的波门时间区间。
优选地，步骤302-H中，计算航迹每个航迹分量每个模型下一次被更新时的状态预测、预测观测、目标存在的先验概率和航迹波门时间，具体地：
假设航迹分量ck共有M个模型，第σ个模型对应的航迹状态均值为
Figure GDA0003812235820000211
航迹状态协方差为Pk|k(σ)。
首先计算每个模型的模型预测概率，表达式为：
Figure GDA0003812235820000212
其中，Γσm表示目标在第k次更新遵循第σ个模型的条件下，在第k+1次更新切换到第m个模型的概率，μk+1|k(m)表示模型m的模型预测概率。
接下来求出混合模型概率和混合状态，表达式为：
Figure GDA0003812235820000213
Figure GDA0003812235820000214
Figure GDA0003812235820000215
其中，
Figure GDA00038122358200002111
表示第k+1次更新是m模型的条件下第k次更新是σ模型的先验概率，简称混合模型概率，
Figure GDA0003812235820000216
表示模型m的混合模型状态，
Figure GDA0003812235820000217
表示模型m的混合模型状态协方差。
对于航迹分量ck的模型m，状态预测的表达式为：
Figure GDA0003812235820000218
Figure GDA0003812235820000219
其中，
Figure GDA00038122358200002110
表示航迹τ第k+1次更新时航迹分量ck的模型m状态预测值，Pk+1|k(ck,m)表示航迹τ第k+1次更新时航迹分量ck的模型m的预测协方差矩阵。F表示状态转移矩阵，FT为F的转置，Q表示过程噪声矩阵。
得到预测观测表达式为：
Figure GDA0003812235820000221
Sk(ck,m)＝HPk+1|k(ck,m)HT+Q；
其中，
Figure GDA0003812235820000222
表示当前航迹第k+1次更新时航迹分量ck的模型m的预测观测，Sk(ck,m)表示当前航迹第k次更新时航迹分量ck的模型m的新息矩阵。
目标存在的先验概率表达式为：
Figure GDA0003812235820000223
其中，
Figure GDA0003812235820000224
表示航迹τ第k+1次更新时目标存在的先验概率，γ为目标存在性的转换概率，可根据实际情况自行设置。
优选地，步骤302中，计算每个航迹分量下一次更新航迹对应的波门时间区间，包括：
根据计算航迹分量ck下一次更新航迹的预测估计得到的新息矩阵，对当前航迹τ的航迹分量ck的模型σ，计算相应的波门方位角的范围[θmin,θmax]；
波门方位角的范围的表达式为：
Figure GDA0003812235820000225
Figure GDA0003812235820000226
其中，g为波门大小，θ0为波门中心方位角，Sij,i,j∈{1,2}为航迹分量ck的模型σ的新息矩阵Sk(ck,σ)的各个元素；
根据波门方位角的范围，计算当前航迹τ的航迹分量ck的模型σ对应的波门时间范围；
对航迹τ的航迹分量ck的模型σ，波门时间范围表达式为：
Figure GDA0003812235820000231
其中，
Figure GDA0003812235820000232
为天线扫描速度，t0为波门中心时间，将航迹分量ck的模型σ的波门时间范围记为
Figure GDA0003812235820000233
其中a为编号，假设航迹τ的航迹分量数为C，模型个数为M，则a的取值范围为[1,CM]。
在这一步需要对每个航迹分量的每一个模型的时间范围进行计算。计算出当前航迹τ的每个航迹分量的每个模型对应的波门时间范围后，将所有航迹分量的所有模型对应的最早波门时间区间起始时刻作为当前航迹的波门时间区间起始时刻ts＝min{ts 1,...,ts CM}，所有航迹分量的所有模型对应的最晚波门时间区间结束时刻作为当前航迹的波门时间区间结束时刻td＝max{td 1,...,td CM}。
步骤301中，根据第一预选观测集建立航迹，并结合航迹的状态进行滤波，计算下一次更新航迹的预测估计及对应的波门时间区间，具体实现方式可参考上述步骤302，在此不再进一步赘述。
请参阅图4至图9(b)本发明还对提出的TWS雷达多目标连续跟踪方法进行了仿真验证，雷达位于笛卡尔坐标系的原点(0m,0m)。有两个目标(即目标1、目标2)，均经历匀速-匀加速-匀速-转弯-匀速-匀减速-转弯-匀速的运动过程。两个目标的初始状态如表2所示，x表示目标位置的横坐标，
Figure GDA0003812235820000234
表示目标在横坐标上的速度，
Figure GDA0003812235820000235
表示目标在横坐标上的加速度，y表示目标位置的纵坐标，
Figure GDA0003812235820000236
表示目标在纵坐标上的速度，
Figure GDA0003812235820000237
表示目标在纵坐标上的加速度。
表2两个目标的初始状态
Figure GDA0003812235820000238
两个目标加速阶段的加速度分别为[-0.1m/s2,-0.2m/s2]和[-0.3m/s2,0.15m/s2]，减速阶段的加速度分别为[0.1m/s2,0.2m/s2]和[0.3m/s2,-0.15m/s2]。在转弯阶段，它们的角速度均为π/20。实验进行100次蒙特卡洛仿真，每次仿真模拟80次更新过程。监视区域中的杂波密度为2.0×10-6/scan/m2，这里的杂波密度也是均匀分布，目标的发现概率PD＝0.9。天线的扫描周期T为1s，以逆时针方向扫描。过程噪声设置为ds＝0.01m/s2，距离、方位角的观测噪声标准差分别为σr＝5m，σθ＝0.01deg。
得到两个目标的运动轨迹如图4所示，本发明利用现有技术中传统的集成航迹分裂跟踪方法(英文缩写为IMM-LMITS，简称ITS)与本发明提供的TWS雷达多目标连续跟踪方法(英文缩写为IMM-LMCITS，简称CITS)进行对比。请参阅图5(a)至图10，从图5(a)和图5(b)可以看出两种方法都成功地跟踪上了目标，但是对图5(a)和图5(b)中确认航迹进行局部放大可知，当目标1穿过扫描边界时，传统的IMM-LMITS方法跟踪航迹出现了断裂的情况，而本发明的连续跟踪IMM-LMCITS方法跟踪航迹没有出现断裂。这说明传统的IMM-LMITS方法在穿过边界时受到边界问题的影响，导致航迹断裂，而本发明提供的方法很好地解决了边界问题。根据图6(a)至图7(b)可知本发明提供的IMM-LMCITS方法的跟踪延迟远小于IMM-LMITS方法。根据图8(a)和图8(b)可知，对于目标1，IMM-LMITS方法的位置RMSE在目标跨过边界后(第24到58次更新期间)远大于IMM-LMITS方法，对于目标2，两者的RMSE十分接近，这是因为它没有跨过边界，没有出现边界问题。由图9(a)和图9(b)可知，两种方法对目标速度的跟踪准确度十分接近。从图10可以看出，在跟踪稳定后，IMM-LMCITS方法一直保持着比较高的真实航迹率，而IMM-LMCITS方法在第24和58次更新附近的真实航迹率有一个大幅度的下降，这是因为目标受边界问题影响，航迹发生了断裂。本发明提出的TWS雷达多目标连续跟踪方法大幅降低了跟踪延迟，并且跟踪延迟的降低并没有导致跟踪准确率下降，且当目标跨过边界时，本发明所提出的连续跟踪方法依然能够持续对目标进行跟踪，不会出现航迹断裂的情况。
特别地，在本发明一些优选的实施方式中，还提供了一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现上述任一实施方式中TWS雷达多目标连续跟踪方法的步骤。
在本发明另一些优选的实施方式中，还提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述任一实施方式中TWS雷达多目标连续跟踪方法的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述TWS雷达多目标连续跟踪方法实施例的流程，在此不再重复说明。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (9)

1.一种TWS雷达多目标连续跟踪方法，其特征在于，包括如下步骤：
步骤(1)、对TWS雷达接收到的每一个观测，选取当前时刻之前获取的观测组成第一预选观测集，所述第一预选观测集内的观测的接收时刻与当前时刻间隔0.8T～1.2T且未被已有的航迹利用，T为TWS雷达的扫描周期；根据所述第一预选观测集建立航迹，并结合航迹的状态进行滤波，计算下一次更新航迹的预测估计及对应的波门时间区间；
步骤(2)、对各个已有的航迹，均采用如下方式在连续跟踪过程中迭代更新航迹的估计及对应的波门时间区间：
每次更新，在当前航迹对应的波门时间区间结束时刻，根据波门时间区间确定采集时间区间，获取采集时间区间内的观测，组成第二预选观测集；在所述采集时间区间结束时刻，利用所述第二预选观测集对当前航迹的状态进行更新，并结合航迹的状态进行滤波，计算下一次更新航迹的预测估计及对应的波门时间区间；
步骤(2)中，根据波门时间区间确定采集时间区间，包括：
若当前航迹与其余已有的航迹对应的波门时间区间均没有交叉，则令采集时间区间与当前航迹对应的波门时间区间相同；
若当前航迹与其余已有的航迹对应的波门时间区间有交叉，且有交叉的航迹对应的最晚波门时间区间结束时刻等于或早于当前航迹对应的波门时间区间结束时刻，则令有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，当前航迹对应的波门时间区间结束时刻为采集时间区间结束时刻；
若当前航迹与其余已有的航迹对应的波门时间区间有交叉，且有交叉的航迹对应的最晚波门时间区间结束时刻晚于当前航迹对应的波门时间区间结束时刻，则令有交叉的航迹及当前航迹对应的最早波门时间区间起始时刻为采集时间区间起始时刻，有交叉的航迹对应的最晚波门时间区间结束时刻为采集时间区间结束时刻。
2.根据权利要求1所述的TWS雷达多目标连续跟踪方法，其特征在于，
步骤(1)中，根据所述第一预选观测集建立航迹时，对所述第一预选观测集中的每一个观测，均判断是否满足如下初始化条件：
Figure FDA0003812235810000021
其中，zx re和zy re分别表示当前时刻接收到的观测的横坐标和纵坐标位置，zx i和zy i分别表示所述第一预选观测集内第i个观测的横坐标和纵坐标位置，Tr为表示当前时刻与接收到所述第一预选观测集内第i个观测的时刻的时间差；vmax表示目标最大速度，verror表示目标速度计算的误差；
若所述第一预选观测集内第i个观测满足初始化条件，则根据当前时刻的观测及第一预选观测集内第i个观测进行航迹初始化，建立航迹，得到航迹的初始状态。
3.根据权利要求1所述的TWS雷达多目标连续跟踪方法，其特征在于，
步骤(2)中，计算下一次更新航迹的预测估计及对应的波门时间区间，包括：
计算所述第二预选观测集中各观测对应的杂波密度；
根据波门对所述第二预选观测集中的观测进行选择，得到用于更新的观测集；
计算所述第二预选观测集中各个观测与当前航迹的似然函数；
结合似然函数，对用于更新的观测集中的观测对应的杂波密度进行调制；
将用于更新的观测集与当前航迹进行数据关联，得到目标存在的后验概率和数据关联后验概率；
进行航迹分量控制；
进行航迹状态更新，得到当前航迹对应的目标轨迹状态；
计算每个航迹分量下一次更新航迹的预测估计及对应的波门时间区间。
4.根据权利要求3所述的TWS雷达多目标连续跟踪方法，其特征在于，
步骤(2)中，计算所述第二预选观测集中各观测对应的杂波密度，包括：
寻找观测zk(i)与所述第二预选观测集Zτ(k)中其余观测之间的第n小的距离rn(i)，n为大于0的整数，zk(i)∈Zτ(k)，k表示更新次数；若所述第二预选观测集Zτ(k)中的观测数小于n+1，则令从第k次更新扫描周期的开始到当前时刻所接收到的观测合集为
Figure FDA0003812235810000031
寻找观测zk(i)在观测合集
Figure FDA0003812235810000032
中第n小的距离rn(i)；
计算观测zk(i)对应的稀疏性，表达式为：
γ(zk(i))＝V(rn(i))/n；
其中，
Figure FDA0003812235810000033
Γ(·)为伽马函数，
Figure FDA0003812235810000035
为空间的维数；
根据稀疏性计算杂波密度，表达式为：
Figure FDA0003812235810000034
5.根据权利要求4所述的TWS雷达多目标连续跟踪方法，其特征在于，
步骤(2)中，第k次更新时，根据波门对所述第二预选观测集Zτ(k)中的观测进行选择，包括：
对于当前航迹τ的航迹分量ck-1的模型σ，选择的观测集yk(ck-1,σ)满足：
Figure FDA0003812235810000041
其中，y为观测集yk(ck-1,σ)中的观测，
Figure FDA0003812235810000042
表示当前航迹τ第k次更新时航迹分量ck-1的模型σ的预测观测，Sk(ck-1,σ)表示当前航迹τ第k次更新时航迹分量ck-1的模型σ的新息矩阵，
Figure FDA0003812235810000043
为Sk(ck-1,σ)的逆矩阵，g为波门大小，表达式为：
Figure FDA0003812235810000044
chi2inv(·)为χ2分布函数的反函数，PG表示波门概率，l表示观测y的维数；
对于当前航迹τ不同的模型，将所有选择的观测集观测合并成用于更新的观测集
Figure FDA0003812235810000045
6.根据权利要求5所述的TWS雷达多目标连续跟踪方法，其特征在于，
步骤(2)中，进行航迹分量控制，包括：
将用于更新的观测集yk中的每个观测yk(i)与每个航迹分量ck-1关联，组成暂时航迹分量；暂时航迹分量的个数为
Figure FDA0003812235810000046
其中Ck-1表示第k-1次更新后航迹分量的个数，mk表示用于更新的观测集yk中的观测个数；
计算每个暂时航迹分量的相对概率，表达式为：
Figure FDA0003812235810000047
其中，
Figure FDA0003812235810000048
表示暂时航迹分量
Figure FDA0003812235810000049
的相对概率，βk(i)表示用来建立暂时航迹分量
Figure FDA00038122358100000410
的观测yk(i)与航迹τ的数据关联后验概率，p(ck-1)表示航迹τ第k-1次更新后航迹分量ck-1的相对概率，
Figure FDA0003812235810000051
表示观测zk(i)与当前航迹τ的航迹分量ck-1的似然函数，
Figure FDA0003812235810000052
表示观测zk(i)与当前航迹τ的似然函数；
进行航迹分量树叶裁剪，若一个暂时航迹分量
Figure FDA0003812235810000053
的相对概率
Figure FDA0003812235810000054
小于门限值，则去掉暂时航迹分量
Figure FDA0003812235810000055
进行航迹分量子树裁剪，从当前第k次更新向前追溯以前更新的航迹分量，对于第κ次更新，保留相对概率p(cκ)最大的航迹分量；其中κ＝k-N，N为子树裁剪深度；
将航迹分量树叶裁剪和航迹分量子树裁剪后，有一个暂时航迹分量子集Θ得以保留，而Θ的补集
Figure FDA0003812235810000056
则被删除；
计算补集
Figure FDA0003812235810000057
的相对概率Δp，表达式为：
Figure FDA0003812235810000058
Θ中的暂时航迹分量为确认航迹分量，用ck表示，航迹分量ck相对概率的表达式为：
Figure FDA0003812235810000059
调整目标存在的后验概率为：
p(χk|Yk)＝p(χk|Yk)(1-Δp)。
7.根据权利要求6所述的TWS雷达多目标连续跟踪方法，其特征在于，
步骤(2)中，计算每个航迹分量下一次更新航迹对应的波门时间区间，包括：
根据计算航迹分量ck下一次更新航迹的预测估计得到的新息矩阵，对当前航迹τ的航迹分量ck的模型σ，计算相应的波门方位角的范围[θmin,θmax]；
波门方位角的范围的表达式为：
Figure FDA0003812235810000061
Figure FDA0003812235810000062
其中，g为波门大小，θ0为波门中心方位角，Sij,i,j∈{1,2}为航迹分量ck的模型σ的新息矩阵Sk(ck,σ)的各个元素；
根据波门方位角的范围，计算当前航迹τ的航迹分量ck的模型σ对应的波门时间范围；
计算出当前航迹τ的每个航迹分量的每个模型对应的波门时间范围后，将所有航迹分量的所有模型对应的最早波门时间区间起始时刻作为当前航迹的波门时间区间起始时刻，所有航迹分量的所有模型对应的最晚波门时间区间结束时刻作为当前航迹的波门时间区间结束时刻。
8.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至7中任一项所述TWS雷达多目标连续跟踪方法的步骤。
9.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至7中任一项所述TWS雷达多目标连续跟踪方法的步骤。
CN202110691909.6A 2021-06-22 2021-06-22 一种tws雷达多目标连续跟踪方法 Active CN113406618B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110691909.6A CN113406618B (zh) 2021-06-22 2021-06-22 一种tws雷达多目标连续跟踪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110691909.6A CN113406618B (zh) 2021-06-22 2021-06-22 一种tws雷达多目标连续跟踪方法

Publications (2)

Publication Number Publication Date
CN113406618A CN113406618A (zh) 2021-09-17
CN113406618B true CN113406618B (zh) 2022-10-25

Family

ID=77682398

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110691909.6A Active CN113406618B (zh) 2021-06-22 2021-06-22 一种tws雷达多目标连续跟踪方法

Country Status (1)

Country Link
CN (1) CN113406618B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106199584A (zh) * 2016-07-05 2016-12-07 电子科技大学 一种基于量测存储的航迹起始方法
CN106597411A (zh) * 2016-12-30 2017-04-26 无锡市雷华科技有限公司 雷达信号处理方法
CN108469604A (zh) * 2018-02-07 2018-08-31 哈尔滨工业大学 Tws雷达空时联合数据模拟方法及系统
CN108490429A (zh) * 2018-02-07 2018-09-04 哈尔滨工业大学 Tws雷达多目标跟踪方法及系统
CN112114308A (zh) * 2019-06-20 2020-12-22 哈尔滨工业大学 一种扇扫雷达空时联合目标跟踪方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106199584A (zh) * 2016-07-05 2016-12-07 电子科技大学 一种基于量测存储的航迹起始方法
CN106597411A (zh) * 2016-12-30 2017-04-26 无锡市雷华科技有限公司 雷达信号处理方法
CN108469604A (zh) * 2018-02-07 2018-08-31 哈尔滨工业大学 Tws雷达空时联合数据模拟方法及系统
CN108490429A (zh) * 2018-02-07 2018-09-04 哈尔滨工业大学 Tws雷达多目标跟踪方法及系统
CN112114308A (zh) * 2019-06-20 2020-12-22 哈尔滨工业大学 一种扇扫雷达空时联合目标跟踪方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
A Track-Before-Detect Method for Rotating Radars with Revisit Interval Uncertainty;Peiyuan Li,et al;《2019 International Conference on Control, Automation and Information Sciences (ICCAIS)》;20191231;p1-6 *
Fast data association approaches for multi-target tracking;Yaotian Zhang,et al;《2016 CIE International Conference on Radar (RADAR)》;20161231;p1-5 *
TWS雷达空时联合多目标跟踪方法;王倩倩;《中国优秀硕士学位论文全文数据库 信息科技辑》;20190115;全文 *
交互作用多模型集成航迹分裂目标跟踪方法;桑航等;《指挥控制与仿真》;20200430;第38-45页 *
空时联合建模与目标跟踪方法;赵颖超;《中国优秀硕士学位论文全文数据库 信息科技辑》;20200215;全文 *

Also Published As

Publication number Publication date
CN113406618A (zh) 2021-09-17

Similar Documents

Publication Publication Date Title
CN101770024B (zh) 多目标跟踪方法
CN111289965B (zh) 一种多目标雷达快速跟踪方法及系统
CN103942536A (zh) 一种迭代更新轨迹模型的多目标跟踪方法
CN107436434B (zh) 基于双向多普勒估计的航迹起始方法
Liu et al. Deinterleaving pulse trains in unconventional circumstances using multiple hypothesis tracking algorithm
CN113406618B (zh) 一种tws雷达多目标连续跟踪方法
CN113406617B (zh) 一种扇扫雷达多目标连续跟踪方法
CN108490429B (zh) Tws雷达多目标跟踪方法及系统
Zhu et al. Particle labeling PHD filter for multi-target track-valued estimates
Musicki et al. Multi-target tracking in clutter without measurement assignment
Chen et al. Group tracking of flock targets in low-altitude airspace
García-Fernández et al. Multitarget tracking using the joint multitrack probability density
Radosavljević et al. Limits of target tracking in heavy clutter
Coraluppi et al. MCMC and MHT Approaches to Multi-INT surveillance
Aoki et al. Suboptimal JPDA for tracking in the presence of clutter and missed detections
CN108646235B (zh) 基于点迹位置聚合识别空间散布固定杂波的方法
JP6400251B1 (ja) 電波探知装置
CN110031797A (zh) 用于被动传感系统对具有非连续特性目标的检测跟踪方法
Li et al. Graph-based Track Stitching Method for Low Altitude Marine Multi-target Tracking
Hanusa et al. Contact clustering and fusion for preprocessing multistatic active sonar data
Barlow et al. New Bayesian track-before-detect design and performance study
Gan et al. Variational Tracking and Redetection for Closely-spaced Objects in Heavy Clutter
Memon Smoothing data association for linear multi-target tracking in clutter
Halimeh et al. Evolutionary resampling for multi-target tracking using probability hypothesis density filter
CN111175738B (zh) 一种基于多模型隶属控制的相控阵雷达目标快速建航方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant