基于期望最大迭代算法的雷达多目标联合检测及跟踪方法与流程

1.本发明属于雷达技术领域，更进一步涉及雷达目标检测技术领域中的一种基于期望最大迭代算法的雷达多目标联合检测及跟踪方法。本发明可用于雷达回波信号图像中多个目标进行联合检测及跟踪。


背景技术：

2.多径环境下雷达多目标联合检测及跟踪是拓宽联合检测及跟踪范围的一个重要方向，在海上侦察、缉毒、交通管制等领域有着重要的作用。由于天波超视距雷达othr(over the horizon radar)中的多径传播现象、多基地配置的目标群的多次延迟回波以及高分辨率雷达的目标扩展特性等原因，雷达每次扫描一个目标可能产生多个量测值，且目标的数量是未知的，所以如何对每个目标产生的多个量测值进行处理，实现联合检测及跟踪是目前雷达领域面临的挑战。
3.b k habtemariam，r tharmarasa，t thayaparan，m mallick，and t kirubarajan在其发表的论文“a multiple
‑
detection joint probabilistic data association filter”(ieee journal of selected topics in signal processing，2013，7(3):461
‑
471)中提出了一种基于多检测概率数据互联滤波器md
‑
pdaf(multiple detection joint probabilistic data association filter)算法的雷达多目标联合检测及跟踪方法。该方法用组合关联事件处理单个目标产生多个量测的概率，通过多目标多路径的联合检测跟踪mtmp
‑
jdt(joint detection and tracking of multitarget and multipath)，检测在杂波和漏检情况下每次扫描的多个目标。该算法计算了观测数据与每个目标之间的关联概率，且认为所有有效回波都可能源于每个特定目标，只是它们源于特定目标的概率不同。该方法能够处理每次扫描单个目标产生的多个量测，且不需要任何关于目标和杂波的先验信息，但是，该方法仍然存在的不足之处是，md
‑
pdaf计算观测数据与每个目标之间的关联概率，当目标和量测数目增多时，需要处理大量组合关联事件，算法的计算量将过于复杂，在实际情况中应用时计算耗时久，耗费大量时间和精力。
4.西安电子科技大学在其申请的专利文献“多雷达系统中目标检测跟踪联合处理方法”(专利申请号：cn201610817863.7，申请公布号：cn106526584a)中公开了一种多雷达系统中目标检测跟踪联合处理方法。该方法建立目标的运动模型以及多雷达系统对目标的量测模型，采用贝叶斯检测器，根据从跟踪器反馈得到目标的跟踪状态，设置恒虚警检测门限，在波门内平均虚警率恒定的前提下，提升目标的平均检测概率和系统的跟踪性能。但是，该方法仍然存在的不足之处是，在复杂环境下目标存在状态和目标跟踪状态之间会存在耦合，导致雷达处理数据时直接给出目标存在状态或目标运动状态存在会带来性能损失，影响多雷达系统中目标检测跟踪结果的准确性。


技术实现要素：

5.本发明的目的在于针对上述现有技术的不足，提出了一种基于期望最大迭代算法
的雷达多目标联合检测及跟踪方法，旨在解决在目标运动状态复杂多变、高杂波率、低信噪比环境下多目标检测及跟踪问题，联合估计目标存在状态和目标运动状态，实现多目标的联合检测及跟踪。
6.实现本发明目的的思想是：基于期望最大迭代算法，对量测数据迭代循环进行预处理操作，降低量测数据的杂波率、提高信噪比，以适应不理想的目标检测及跟踪环境。同时根据hybird forwardandbackward算法计算待检测及跟踪目标中的目标存在状态概率，以上次迭代的目标存在状态与目标运动状态的估计结果作为当前迭代的输入，使得在复杂多变的目标运动状态环境下的目标检测性能随着迭代次数逐步提高。再利用目标存在状态的概率和目标运动状态的互耦关系，由当前迭代计算出的目标存在状态概率得到当前迭代的目标运动状态，避免目标存在状态与目标运动状态独立估计带来的估计性能损失。在多目标跟踪场景中目标存在状态是未知信息的情况下，实现多目标联合检测及跟踪。
7.本发明的具体步骤如下：
8.步骤1，对目标存在状态和运动状态进行初始化：
9.(1a)将k
‑
l
‑
1时刻目标运动状态代入目标动力学方程中，得到k
‑
l到k时刻目标运动状态的初始状态
10.(1b)将k
‑
l
‑
1时刻目标存在状态代入目标存在方程，得到k
‑
l到k时刻目标存在状态的初始状态其中，k表示当前迭代时平滑操作的截止时刻，k的取值范围为5到100，l表示当前迭代时平滑操作的窗口时间，取值为5；
11.步骤2，对量测数据进行预处理：
12.(2a)从[k
‑
l,k]中选择一个未被选过的时刻；
[0013]
(2b)从[1,t
j
]中选择一个未被选过的目标，其中，t
j
表示当前时刻待检测及跟踪目标的总数；
[0014]
(2c)从[1,m]中选择一个未被选过的路径，其中，m表示待检测及跟踪目标的量测路径的总数；
[0015]
(2d)按照下式，计算所选定时刻所选目标在所选定路径的预测量测值：
[0016][0017]
其中，表示所选择的第j个时刻第t个待检测及跟踪目标在第m个量测路径的预测量测值，h
m
(
·
)表示路径量测方程，表示上一次迭代得到的目标运动状态，若当前迭代为第一次迭代时，则为初始化时的目标运动状态；
[0018]
(2e)计算所选时刻量测数据中的每个量测值与预测量测值之间的马氏距离，若该马氏距离小于或等于η，则将该量测值存储至有效量测集v
t,m,j
后执行步骤(2f)；否则，直接执行步骤(2f)；
[0019]
(2f)判断是否选完[1,m]中所有路径，若是，则执行步骤(2g)，否则，执行步骤(2c)；
[0020]
(2g)判断是否选完[1,t
j
]中所有目标，若是，则执行步骤(2h)，否则，执行步骤
(2b)；
[0021]
(2h)计算所选时刻预处理后的总量测集；
[0022]
(2i)判断是否选完[k
‑
l,k]中所有时刻，若是，则执行步骤(3)；否则，执行步骤(2a)；
[0023]
步骤3，根据hybirdforwardandbackward算法计算目标存在状态的概率：
[0024][0025][0026][0027][0028]
其中，p(
·
)表示概率函数，α
j
(
·
)表示前项目标存在状态概率，β
j
(
·
)表示后项目标存在状态概率，s
t,j
表示中第j时刻第t个待检测及跟踪目标的目标存在状态，σ表示求和操作，γ用来判断待检测及跟踪的目标是否存在，当γ＝0时表示待检测及跟踪的目标不存在，当γ＝1时表示待检测及跟踪的目标存在，α
j
‑1(
·
)表示当前跟踪前一时刻的目标存在状态概率，表示已知第j
‑
1时刻第t个待检测及跟踪目标的目标存在状态为γ时，第j时刻第t个待检测及跟踪目标的目标存在状态的概率，s
t,j
‑1表示中第j
‑
1时刻第t个待检测及跟踪目标的目标存在状态，exp表示以自然常数e为底的指数操作，ξ
t,j
‑1(
·
)表示当前跟踪上一时刻关联事件与目标存在状态的关系，β
j 1
(
·
)表示当前跟踪下一时刻目标存在状态概率，表示已知第j时刻第t个待检测及跟踪目标的目标存在状态为γ时，第j 1时刻第t个待检测及跟踪目标的目标存在状态的概率，ξ
t,j
(
·
)表示当前跟踪时刻关联事件与目标存在状态的关系，ε[
·
]表示数学期望操作，表示第j时刻所有关联事件集中第t个待检测及跟踪目标在第m种量测路径下的关联标志，θ
j
表示第j时刻每个待检测及跟踪目标与每种量测路径、每个量测数据中量测值的所有关联事件集，log(
·
)表示取对数操作，表示在第m种量测路径下的检测概率，表示上次迭代估计得到的目标存在状态，若为第一次迭代，为初始化时的目标存在状态；
[0029]
步骤4，按照下式，计算当前迭代时待检测及跟踪的目标与跟踪目标的量测路径、量测数据中的量测值三者之间的关联概率：
[0030][0031]
其中，表示当前迭代时第j时刻第t个待检测及跟踪目标与第m种量测
路径、量测数据中第n个量测值之间的关联概率，p(
·
|
·
)表示“|”的右侧为已知数据时左侧未知关联事件发生的概率，表示第j时刻第t个待检测及跟踪目标、第m种跟踪目标的量测路径和量测数据中第n个量测值三者之间的一个关联事件，表示从k
‑
l到k时刻的所有量测数据，表示上次迭代估计得到的目标运动状态，若为第一次迭代，则为初始化时的目标运动状态；
[0032]
步骤5，根据目标存在状态与目标运动状态的互耦关系，按照下式计算当前迭代时的目标运动状态及其协方差矩阵：
[0033][0034][0035]
其中，表示当前迭代时第j时刻第t个待检测及跟踪目标的目标运动状态，σ(
·
)表示求和操作，表示当前迭代时第j时刻第t个待检测及跟踪目标与第m种跟踪目标的量测路径的关联概率，表示第j时刻第t个待检测及跟踪目标在第m种跟踪目标的量测路径的目标运动状态平滑值的协方差矩阵，表示第j时刻第t个待检测及跟踪目标的目标运动状态平滑值，φ
t,j
(r 1)表示当前迭代时第j时刻第t个待检测及跟踪目标的目标运动状态的协方差矩阵；
[0036]
步骤6，判断迭代误差θ是否小于或等于ε，若是，则执行步骤7，否则，将迭代次数加1作为当前迭代次数后执行步骤2；
[0037]
步骤7，输出当前迭代的目标运动状态和目标存在状态。
[0038]
本发明与现有技术相比具有以下优点：
[0039]
第一，由于本发明基于期望最大迭代算法，对量测数据迭代循环进行预处理操作，克服现有技术无法精确去除杂波以及未保留足够的有效量测的问题，使得本发明保留有效量测的同时具有降低计算复杂度，提高迭代收敛速度，提高系统稳定性，更符合雷达应用真实情况的优点。
[0040]
第二，由于本发明基于期望最大迭代算法，在估计目标存在状态时，以上次迭代的目标存在状态与目标运动状态的估计结果作为当前迭代的输入，解决了现有方法中单次估计精度不高的问题，再根据hybird forward and backward算法计算待检测及跟踪目标中的目标存在状态概率，解决了实际检测及跟踪场景中目标时隐时现、速度捷变等复杂运动状态的问题，使得本发明提高了在复杂目标运动场景下的目标检测概率。
[0041]
第三，由于本发明利用了目标存在状态的概率和目标运动状态的互耦关系，在估计目标运动状态时，由当前迭代计算出的目标存在状态概率得到当前迭代的目标运动状态，解决了现有方法中未利用互耦关系而直接对目标存在状态与目标运动状态进行独立估
计的问题，使得本发明提高了在目标存在状态未知的情况下的目标跟踪精度。
附图说明
[0042]
图1为本发明的流程图；
[0043]
图2为目标运动轨迹图；
[0044]
图3为4种量测路径下量测数据和杂波的分布图；
[0045]
图4是两组算法检测及跟踪目标的轨迹图。
具体实施方式
[0046]
下面结合附图和实施例对本发明做进一步的描述。
[0047]
参照图1，本发明实现的具体实施步骤做进一步的描述。
[0048]
步骤1，对目标存在状态和运动状态进行初始化。
[0049]
将k
‑
l
‑
1时刻目标运动状态代入目标动力学方程中，得到k
‑
l到k时刻目标运动状态的初始状态
[0050]
将k
‑
l
‑
1时刻目标存在状态代入目标存在方程，得到k
‑
l到k时刻目标存在状态的初始状态其中，k表示当前迭代时平滑操作的截止时刻，k的取值范围为5到100，l表示当前迭代时平滑操作的窗口时间，取值为5。
[0051]
步骤2，对量测数据进行预处理。
[0052]
第一步，从[k
‑
l,k]中选择一个未被选过的时刻。
[0053]
第二步，从[1,t
j
]中选择一个未被选过的目标，其中，t
j
表示当前时刻待检测及跟踪目标的总数。
[0054]
第三步，从[1,m]中选择一个未被选过的路径，其中，m表示待检测及跟踪目标的量测路径的总数。
[0055]
第四步，按照下式，计算所选定时刻所选目标在所选定路径的预测量测值：
[0056][0057]
其中，表示所选择的第j个时刻第t个待检测及跟踪目标在第m个量测路径的预测量测值，h
m
(
·
)表示路径量测方程，是代表目标运动状态与量测值的数学函数关系，该关系是由实际检测及跟踪的场景确定，不同检测及跟踪的场景对应h
m
(
·
)的不同表达式，表示上一次迭代得到的目标运动状态，若当前迭代为第一次迭代时，则为初始化时的目标运动状态。
[0058]
第五步，计算所选时刻量测数据中的每个量测值与预测量测值之间的马氏距离，若该马氏距离小于或等于η，则将该量测值存储至有效量测集v
t,m,j
后执行本步骤的第六步；否则，直接执行本步骤的第六步。
[0059]
第六步，判断是否选完[1,m]中所有路径，若是，则执行本步骤的第七步，否则，执行本步骤的第三步。
[0060]
第七步，判断是否选完[1,t
j
]中所有目标，若是，则执行本步骤的第八步，否则，执
行本步骤的第二步。
[0061]
第八步，计算所选时刻预处理后的总量测集：
[0062][0063]
其中，v
j
表示所选时刻预处理后的总量测集，∪表示取并集操作，max{
·
}表示取最大值操作，v
t,1,j
表示第j个时刻第t个待检测及跟踪目标在第1个量测路径下的有效量测集，v
t,m,j
表示第j个时刻第t个待检测及跟踪目标在第m个量测路径下的有效量测集。
[0064]
第九步，判断是否选完[k
‑
l,k]中所有时刻，若是，则执行步骤3；否则，执行本步骤的第一步。
[0065]
步骤3，根据hybird forward and backward算法计算目标存在状态的概率。
[0066][0067][0068][0069][0070]
其中，p(
·
)表示概率函数，α
j
(
·
)表示前项目标存在状态概率，β
j
(
·
)表示后项目标存在状态概率，s
t,j
表示中第j时刻第t个待检测及跟踪目标的目标存在状态，σ表示求和操作，γ用来判断待检测及跟踪的目标是否存在，当γ＝0时表示待检测及跟踪的目标不存在，当γ＝1时表示待检测及跟踪的目标存在，α
j
‑1(
·
)表示当前跟踪前一时刻的目标存在状态概率，表示已知第j
‑
1时刻第t个待检测及跟踪目标的目标存在状态为γ时，第j时刻第t个待检测及跟踪目标的目标存在状态的概率，s
t,j
‑1表示中第j
‑
1时刻第t个待检测及跟踪目标的目标存在状态，exp表示以自然常数e为底的指数操作，ξ
t,j
‑1(
·
)表示当前跟踪上一时刻关联事件与目标存在状态的关系，β
j 1
(
·
)表示当前跟踪下一时刻目标存在状态概率，表示已知第j时刻第t个待检测及跟踪目标的目标存在状态为γ时，第j 1时刻第t个待检测及跟踪目标的目标存在状态的概率，ξ
t,j
(
·
)表示当前跟踪时刻关联事件与目标存在状态的关系，ε[
·
]表示数学期望操作，表示第j时刻所有关联事件集中第t个待检测及跟踪目标在第m种量测路径下的关联标志，θ
j
表示第j时刻每个待检测及跟踪目标与每种量测路径、每个量测数据中量测值的所有关联事件集，log(
·
)表示取对数操作，表示在第m种量测路径下的检测概率，表示上次迭代估计得到的目标存在状态，若为第一次迭代，为初始化时的目标存在状态；
[0071]
步骤4，按照下式，计算当前迭代时待检测及跟踪的目标与跟踪目标的量测路径、量测数据中的量测值三者之间的关联概率：
[0072][0073]
其中，
[0074][0075][0076][0077]
表示当前迭代时第j时刻第t个待检测及跟踪目标与第m种量测路径、量测数据中第n个量测值之间的关联概率，p(
·
|
·
)表示“|”的右侧为已知数据时左侧未知关联事件发生的概率，表示第j时刻第t个待检测及跟踪目标、第m种跟踪目标的量测路径和量测数据中第n个量测值三者之间的一个关联事件，表示从k
‑
l到k时刻的所有量测数据，表示上次迭代估计得到的目标运动状态，若为第一次迭代，则为初始化时的目标运动状态，表示一个关联事件的先验概率，表示关联事件与量测数据之间的条件概率，λ表示量测数据中所含杂波的杂波密度，n
j
表示第j时刻量测数据中量测值总数，！表示阶乘，n(
·
)表示高斯概率函数，z
n,j
表示第j时刻量测数据中第n个量测值，x
t,j
表示第j时刻第t个待检测及跟踪目标，r
m,j
表示第j时刻量测数据中在第m种量测路径下噪声的协方差矩阵。
[0078]
步骤5，根据目标存在状态与目标运动状态的互耦关系，按照下式计算当前迭代时的目标运动状态及其协方差矩阵：
[0079][0080][0081]
其中，
[0082][0083]
[0084][0085]
其中，表示当前迭代时第j时刻第t个待检测及跟踪目标的目标运动状态，σ(
·
)表示求和操作，表示当前迭代时第j时刻第t个待检测及跟踪目标与第m种跟踪目标的量测路径的关联概率，表示第j时刻第t个待检测及跟踪目标在第m种跟踪目标的量测路径的目标运动状态平滑值的协方差矩阵，表示第j时刻第t个待检测及跟踪目标的目标运动状态平滑值，φ
t,j
(r 1)表示当前迭代时第j时刻第t个待检测及跟踪目标的目标运动状态的协方差矩阵；表示当前迭代时第j时刻第t个待检测及跟踪目标在第m种量测路径下的前项估计的目标运动状态，表示当前迭代时第j时刻第t个待检测及跟踪目标在第m种量测路径下的后项估计的目标运动状态，t个待检测及跟踪目标在第m种量测路径下的后项估计的目标运动状态，与皆可通过扩展卡尔曼滤波器或无损卡尔曼滤波器获得，滤波器的输入为当前迭代时第j时刻第t个目标与第m种跟踪目标的量测路径相关联的虚拟量测值与当前迭代时第j时刻第t个目标与第m种跟踪目标的量测路径相关联的虚拟量测值的协方差矩阵分别由下式得到：
[0086][0087][0088]
其中，代表当前迭代时第j时刻第t个目标与第m种跟踪目标的量测路径相关联的量测值的协方差矩阵。
[0089]
步骤6，计算迭代误差：
[0090][0091]
其中，表示当前迭代估计出的目标运动状态，表示上次迭代估计出的目标运动状态；判断迭代误差θ是否小于或等于ε，若是，则执行步骤7，否则，将迭代次数加1作为当前迭代次数后执行步骤2。
[0092]
步骤7，输出当前迭代的目标运动状态和目标存在状态。
[0093]
下面结合附图对本发明的效果做进一步的描述。
[0094]
1.仿真条件
[0095]
本发明的仿真运行系统为cpu主频3.3ghz、内存8g、显卡gtx960、windows操作系统，仿真软件采用matlab。
[0096]
本发明仿真实验的电离层模型采用e层和f层双层电离层模型，待检测及跟踪目标的量测路径总数为4，雷达收发距离为100km，噪声协方差矩阵为diag(25,1e
‑
6,9e
‑
6)，门概率p
g
＝0.971。
[0097]
2.仿真内容与结果分析：
[0098]
本发明的仿真实验有两个。
[0099]
仿真实验1是采用期望最大迭代算法对一台天波超视距雷达othr检测及跟踪区域内的6个运动目标及杂波进行联合检测及跟踪，检测概率内的6个运动目标及杂波进行联合检测及跟踪，检测概率杂波数量为125，通过在matlab上进行仿真，得到如图2和图3所示的仿真实验结果。
[0100]
图2表示采用期望最大迭代算法对一台天波超视距雷达othr检测及跟踪区域内的6个运动目标进行联合检测及跟踪，通过在matlab上进行仿真，得到的目标运动轨迹图；
[0101]
图3表示相同仿真条件下得到的4种量测路径下量测数据和杂波的分布图。
[0102]
图2中的横坐标表示方位，单位为度，纵坐标表示距离，单位为km，用标示目标运动的起点，
“×”
标示目标运动的终点，t1[1,20]表示第1个目标在时刻1出现时刻20消失；t2[1,20]表示第2个目标在时刻1出现时刻20消失；t3[1,30]表示第3个目标在时刻1出现时刻30消失；t4[1,20]&[20,35]表示第4个目标在时刻1出现时刻20消失，时刻20出现时刻35消失；t5[5,20]&[30,40]表示第5个目标在时刻5出现时刻20消失，时刻30出现时刻40消失；t6[1,10]&[15,25]&[30,40]表示第6个目标在时刻1出现时刻10消失，时刻15出现时刻25消失，时刻30出现时刻40消失。
[0103]
图2由和
“×”
标示的6条直线段或虚线段表示了6个目标的运动轨迹，说明采用期望最大迭代算法能够检测到目标。
[0104]
图3中的横坐标表示方位，单位为度，纵坐标表示距离，单位为km，以标示的所有点表示雷达通过f层电离层发射信号以及通过f层电离层接收回波的此条量测路径内的量测数据，以“*”标示的所有点表示雷达通过f层电离层发射信号以及通过e层电离层接收回波的此条量测路径内的量测数据，以“ ”标示的所有点表示雷达通过e层电离层发射信号以及通过f层电离层接收回波此条量测路径内的量测数据，以
“☆”
标示的所有点表示雷达通过e层电离层发射信号以及通过e层电离层接收回波的此条量测路径内的量测数据，“·”表示天波超视距雷达othr检测及跟踪区域内4种量测路径下量测数据中所含的杂波。
[0105]
图3由“*”、“ ”、
“☆”
在检测区域中的分布表示在4种量测路径下都存在量测数据，说明期望最大迭代算法在杂波密度较高时能够检测到4种量测路径下的量测数据。
[0106]
仿真实验2是分别采用期望最大迭代算法以及现有技术中的多检测概率数据互联滤波器md
‑
pdaf算法对一台天波超视距雷达othr检测及跟踪区域内的6个运动目标及杂波进行联合检测及跟踪，通过在matlab上进行仿真，得到如图4所示的仿真实验结果。
[0107]
所述的现有技术中的多检测概率数据互联滤波器md
‑
pdaf算法是指：
[0108]
b k habtemariam，r tharmarasa，t thayaparan，m mallick，and t kirubarajan在其发表的论文“a multiple
‑
detection joint probabilistic data association filter”(ieee journal of selected topics in signal processing，2013，7(3):461
‑
471)中提出的算法。
[0109]
图4(a)表示采用现有技术中的多检测概率数据互联滤波器md
‑
pdaf算法对一台天
波超视距雷达othr检测及跟踪区域内的6个运动目标及杂波进行联合检测及跟踪，通过在matlab上进行仿真，得到的仿真实验结果，图4(b)表示采用期望最大迭代算法对一台天波超视距雷达othr检测及跟踪区域内的6个运动目标及杂波进行联合检测及跟踪，通过在matlab上进行仿真，得到的仿真实验结果。
[0110]
图4中的横坐标表示方位角，单位为度，纵坐标表示距离范围，单位为km，图4中的以标示的所有点表示雷达通过f层电离层发射信号以及通过f层电离层接收回波的此条量测路径内的量测数据，以“*”标示的所有点表示雷达通过f层电离层发射信号以及通过e层电离层接收回波的此条量测路径内的量测数据。
[0111]
图4(a)由多个“*”近似组成的6段直线分布表示6个被检测及跟踪到的目标，检测及跟踪区域中其它单个“*”表示虚假航迹，可见，6个运动目标在两种量测路径下都被成功跟踪，但6段直线中和“*”重合度不够好，且虚假航迹较多，说明md
‑
pdaf算法在两种量测路径下检测及跟踪同一个目标会存在偏差，检测及跟踪结果不够准确；图4(b)由多个“*”近似组成的6段直线分布表示6个被检测及跟踪到的目标，可见，6个运动目标在两种量测路径下都被成功跟踪，且6段直线中和“*”重合度较高，虚假航迹很少，说明期望最大迭代算法在两种量测路径下检测及跟踪同一个目标基本没有偏差，检测及跟踪结果更准确。