CN104156984A

CN104156984A - PHD (Probability Hypothesis Density) method for multi-target tracking in uneven clutter environment

Info

Publication number: CN104156984A
Application number: CN201410381497.6A
Authority: CN
Inventors: 杨峰; 史玺; 王永齐; 梁彦; 潘泉; 刘柯利; 陈昊; 史志远; 王碧垚
Original assignee: Northwestern Polytechnical University
Current assignee: Northwestern Polytechnical University
Priority date: 2014-05-30
Filing date: 2014-08-02
Publication date: 2014-11-19
Anticipated expiration: 2034-08-02
Also published as: CN104156984B; CN104019816A

Abstract

An embodiment of the present invention provides a probability hypothesis density method for multi-target tracking in an uneven clutter environment, which relates to the field of intelligent information processing, and can improve the accuracy of calculation results while reducing the amount of calculation. The method includes: if the monitoring area is a sparse clutter area, directly perform standard PHD prediction and update to obtain the estimated number of targets and target status; if the monitoring area is a dense clutter area, then search for a convex hull to determine the clutter area; and When the PHD prediction is performed moment by moment, the echoes not included in the clutter area are selected for calculation. When the PHD is updated, the measurement of the real target is added to the echo that just falls into the clutter area; the posterior intensity The Gaussian components in are pruned and merged to obtain the target strength; then calculate the weight sum of all Gaussian components in the target strength to obtain the estimated number of targets in the monitoring area; extract the Gaussian components with a weight greater than τ in the target strength, as the estimated target state.

Description

A Probabilistic Hypothesis Density Method for Multi-Target Tracking in Inhomogeneous Clutter Environment

技术领域technical field

本发明涉及智能信息处理领域，尤其涉及一种不均匀杂波环境下多目标跟踪的概率假设密度方法。The invention relates to the field of intelligent information processing, in particular to a probability hypothesis density method for multi-target tracking in an uneven clutter environment.

背景技术Background technique

多目标跟踪是研究从目标量测和杂波中估计目标数目与各个目标状态的方法。传统的处理多目标跟踪问题的方法大都是基于数据关联的，需要建立量测和目标的对应关系，如最近邻法(NN)，联合概率数据关联算法(JPDA)，多假设跟踪算法(MHT)。近年来，基于随机有限集(Random Finite Set，RFS)理论的多目标跟踪算法得到很大的关注。Mahler提出的概率假设密度(Probability HypothesisDensity，PHD)滤波算法就是一种在随机有限集框架下，相比完全多目标贝叶斯滤波器，运算复杂度得到有效降低的算法。Multi-target tracking is a method of estimating the number of targets and the state of each target from target measurements and clutter. Most of the traditional methods to deal with multi-target tracking problems are based on data association, and need to establish the corresponding relationship between measurement and target, such as nearest neighbor method (NN), joint probabilistic data association algorithm (JPDA), multiple hypothesis tracking algorithm (MHT) . In recent years, multi-target tracking algorithms based on Random Finite Set (RFS) theory have received great attention. The Probability Hypothesis Density (PHD) filtering algorithm proposed by Mahler is an algorithm that effectively reduces the computational complexity compared with the complete multi-objective Bayesian filter under the framework of random finite sets.

在传统的PHD滤波器中，通常假设在可观测的范围内杂波是均匀分布的，但在实际环境中杂波的空间分布在观测区域是不均匀的，杂波区可分为稀疏杂波区和密集杂波区，在密集杂波区，虚警的数目远远大于真实目标的量测，这种情况下进行多目标跟踪是一个困难的课题。因此，当杂波分布不均匀时，传统的PHD滤波器的跟踪性能会急剧下降。针对此问题，现有技术中有一种方法是在PHD粒子滤波的过程中，通过设定门限，将位于门限之外的观测值滤除，仅利用位于门限内的观测值进行粒子权重的更新。还有一种方法是提出一种应用目标起始和维持规则辅助的SMC-PHD滤波方法，每个粒子除了状态矢量和权值外，还增加了两个辅助变量，一个用来表示粒子的存活年龄，一个用来表示最近N个周期粒子受到观测值更新的情况。In the traditional PHD filter, it is usually assumed that the clutter is uniformly distributed in the observable range, but in the actual environment, the spatial distribution of clutter is not uniform in the observation area, and the clutter area can be divided into sparse clutter In the dense clutter area, the number of false alarms is far greater than the measurement of real targets. In this case, multi-target tracking is a difficult task. Therefore, when the clutter distribution is not uniform, the tracking performance of conventional PHD filters will drop sharply. To solve this problem, there is a method in the prior art that in the process of PHD particle filtering, by setting a threshold, the observations outside the threshold are filtered out, and only the observations inside the threshold are used to update the particle weights. Another method is to propose an SMC-PHD filtering method assisted by the application of target initiation and maintenance rules. In addition to the state vector and weight, each particle also adds two auxiliary variables, one is used to represent the survival age of the particle , one is used to represent the situation that the particle has been updated by the observation value in the last N periods.

上述这些方法，都是针对密集杂波环境下目标跟踪的PHD方法，而在既有稀疏杂波又有密集杂波的复杂的杂波环境下应用这些方法，就会导致在稀疏杂波区和密集杂波区计算量都很大，不能实现快速运算，然而在稀疏杂波区计算量本不用那么大，应用这些方法反而导致计算量大，而且得到的结果也不一定理想。The above-mentioned methods are all PHD methods for target tracking in dense clutter environments, and applying these methods in a complex clutter environment with both sparse clutter and dense clutter will result in the sparse clutter area and The calculation amount in the dense clutter area is very large, and fast calculation cannot be realized. However, the calculation amount in the sparse clutter area is not so large. The application of these methods will lead to a large amount of calculation, and the results obtained may not be ideal.

发明内容Contents of the invention

本发明的实施例提供一种不均匀杂波环境下多目标跟踪的概率假设密度方法，可以在计算量减小的情况下，提高计算结果的准确性。Embodiments of the present invention provide a probability hypothesis density method for multi-target tracking in an uneven clutter environment, which can improve the accuracy of calculation results while reducing the amount of calculation.

为达到上述目的，本发明的实施例采用如下技术方案：In order to achieve the above object, embodiments of the present invention adopt the following technical solutions:

一种不均匀杂波环境下多目标跟踪的概率假设密度方法，包括：A probability hypothesis density method for multi-target tracking in a heterogeneous clutter environment, including:

判断监视区域内n帧累积的回波数是否大于门限值ε，其中，所述n，ε为预设数值；Judging whether the number of echoes accumulated in n frames in the monitoring area is greater than the threshold value ε, wherein the n and ε are preset values;

若否，则直接进行标准PHD预测及PHD更新，获得估计的目标数以及目标状态；If not, perform standard PHD prediction and PHD update directly to obtain the estimated number of targets and target status;

若是，则寻找凸包，确定杂波区；If so, find the convex hull and determine the clutter area;

初始化，获得初始时刻的后验强度，所述初始时刻为0时刻；Initialize to obtain the posterior strength at the initial moment, where the initial moment is 0 moment;

针对k时刻测得的回波值，选择未包含在所述杂波区的回波进行PHD预测计算，根据k-1时刻的后验强度，获得k时刻的预测强度，其中，所述k为大于等于1的整数；For the echo value measured at time k, select echoes not included in the clutter area to perform PHD prediction calculation, and obtain the predicted strength at time k according to the posterior strength at time k-1, where k is an integer greater than or equal to 1;

根据所述k时刻的预测强度，用更新的回波进行PHD更新，获得k时刻的后验强度，其中，所述更新的回波包括k时刻测得的回波值中所述未包含在所述杂波区的回波，以及k时刻真实目标的量测刚好落入所述杂波区的回波；所述后验强度为高斯混合形式，包括各个高斯成分的权值、均值和协方差；According to the predicted strength at time k, the updated echo is used to perform PHD update to obtain the posterior strength at time k, wherein the updated echo includes the echo value measured at time k that is not included in the The echo in the clutter area, and the measurement of the real target at time k just falls into the echo in the clutter area; the posterior intensity is in the form of a Gaussian mixture, including the weight, mean and covariance of each Gaussian component ;

将所述后验强度中权值小于T的高斯成分滤除，将高斯成分间的距离小于U的高斯成分进行合并，获得目标强度；其中，U≤10，T≤10^-5；Filter out the Gaussian components whose weights are less than T in the posterior intensity, and combine the Gaussian components whose distance between the Gaussian components is less than U to obtain the target intensity; wherein, U≤10, T≤10 ^-5 ;

计算所述目标强度中所有高斯成分的权值和，得到所述监视区域内的估计目标数；提取所述目标强度中权值大于τ的高斯成分，作为估计的目标状态，其中，所述τ≥0.5。Calculate the weight sum of all Gaussian components in the target intensity to obtain the estimated number of targets in the monitoring area; extract the Gaussian components with a weight greater than τ in the target intensity as the estimated target state, where the τ ≥0.5.

可选的，所述寻找凸包，确定杂波区，包括：Optionally, the searching for the convex hull and determining the clutter area include:

应用近似传播(AP)聚类算法对所述监视区域内的n帧累积的回波进行聚类，获得聚类结果；Applying an approximate propagation (AP) clustering algorithm to cluster the accumulated echoes of n frames in the monitoring area to obtain a clustering result;

对于所述聚类结果中聚类点数大于η的区域，用Graham scan方法寻找凸包，确定杂波区，所述η根据仿真场景而定。For the region where the number of clustering points is greater than η in the clustering result, the Graham scan method is used to find the convex hull to determine the clutter area, and the η is determined according to the simulation scene.

上述技术方案提供的不均匀杂波环境下多目标跟踪的概率假设密度方法，利用AP聚类算法，首先对监视区域内多帧累积的回波进行聚类，用寻找凸包的方法确定杂波区，然后再逐观测时刻进行PHD预测和更新，在PHD预测时，不用杂波区内的回波，但当进行PHD更新时，需要将发现的真实目标的量测刚好落入杂波区中的回波取出进行PHD更新，这样可以保证算法的高精度，使估计结果更准确。另外，当监视区域内的回波数大于某个门限时才进行聚类和确定杂波区的运算，否则不作此运算直接进行标准的PHD预测和更新，这样计算过程简单，对于不均匀杂波区可以实现快速运算，从而实现不均匀杂波环境下PHD方法跟踪性能的提升。The probability assumption density method of multi-target tracking in the uneven clutter environment provided by the above technical solution uses the AP clustering algorithm to cluster the echoes accumulated in multiple frames in the monitoring area first, and then determines the clutter by finding the convex hull area, and then perform PHD prediction and update time by observation time. In the PHD prediction, the echo in the clutter area is not used, but when the PHD update is performed, the measurement of the real target found should just fall into the clutter area. The echo is taken out for PHD update, which can ensure the high precision of the algorithm and make the estimation result more accurate. In addition, when the number of echoes in the monitoring area is greater than a certain threshold, the calculation of clustering and determining the clutter area is performed; otherwise, the standard PHD prediction and update is directly performed without this operation, so that the calculation process is simple. For the uneven clutter area Fast calculation can be realized, so as to improve the tracking performance of the PHD method in the uneven clutter environment.

附图说明Description of drawings

图1为本发明实施例提供的一种不均匀杂波环境下多目标跟踪的概率假设密度方法的流程示意图；FIG. 1 is a schematic flowchart of a probability hypothesis density method for multi-target tracking in an uneven clutter environment provided by an embodiment of the present invention;

图2为本发明实施例提供的一种AP聚类算法数据点之间消失传递示意图；Fig. 2 is a schematic diagram of vanishing transfer between data points of an AP clustering algorithm provided by an embodiment of the present invention;

图3为本发明实施例提供的一种应用Graham scan算法寻找凸包的示意图；Fig. 3 is a kind of application Graham scan algorithm provided by the embodiment of the present invention to find the schematic diagram of convex hull;

图4为本发明实施例提供的一种不均匀杂波环境下进行多目标跟踪的仿真场景图；FIG. 4 is a simulation scene diagram of multi-target tracking in an uneven clutter environment provided by an embodiment of the present invention;

图5为本发明实施例提供的一种应用本发明提供的方法仿真出来的状态估计图；FIG. 5 is a state estimation diagram simulated by applying the method provided by the present invention provided by the embodiment of the present invention;

图6为本发明实施例提供的一种应用本发明提供的方法仿真出来的目标数估计图；Fig. 6 is a target number estimation diagram simulated by applying the method provided by the present invention provided by the embodiment of the present invention;

图7为本发明实施例提供的一种应用高斯混合PHD方法仿真出来的状态估计图；FIG. 7 is a state estimation diagram simulated by using a Gaussian mixture PHD method provided by an embodiment of the present invention;

图8为本发明实施例提供的一种应用高斯混合PHD方法仿真出来的目标数估计图。FIG. 8 is an estimation diagram of the number of targets simulated by using the Gaussian mixture PHD method provided by the embodiment of the present invention.

具体实施方式Detailed ways

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only some, not all, embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without making creative efforts belong to the protection scope of the present invention.

本发明实施例提供了一种不均匀杂波环境下多目标跟踪的概率假设密度方法，如图1所示，所述方法包括以下步骤：An embodiment of the present invention provides a probability hypothesis density method for multi-target tracking in an uneven clutter environment. As shown in FIG. 1 , the method includes the following steps:

101、判断监视区域内n帧累积的回波数是否大于ε。101. Determine whether the number of echoes accumulated in n frames in the monitoring area is greater than ε.

其中，所述n，ε为预设数值；所述ε是依据监视区域内的回波数而由使用者预设的，所述n也是由使用者自己预设的。Wherein, the n and ε are preset values; the ε is preset by the user according to the number of echoes in the monitoring area, and the n is also preset by the user.

若监视区域内n帧累积的回波数大于ε，则表明所述监视区域为密集杂波区域，进行步骤103-107；若监视区域内n帧累积的回波数小于ε，则表明所述监视区域为稀疏杂波区域，进行步骤102。If the number of echoes accumulated in n frames in the monitoring area is greater than ε, it indicates that the monitoring area is a dense clutter area, and steps 103-107 are performed; if the number of echoes accumulated in n frames in the monitoring area is less than ε, it indicates that the monitoring area is To sparse the clutter area, go to step 102 .

102、直接进行标准PHD预测及PHD更新，获得估计的目标数以及目标状态。102. Directly perform standard PHD prediction and PHD update, and obtain estimated target number and target status.

在这里所述的标准PHD预测及PHD更新，即为现有技术中的传统的PHD预测及PHD更新，本领域技术人员都清楚了解其计算过程，故在此不再进行赘述。The standard PHD prediction and PHD update mentioned here are the traditional PHD prediction and PHD update in the prior art, and those skilled in the art are well aware of the calculation process, so it will not be repeated here.

103、寻找凸包，确定杂波区。103. Find the convex hull and determine the clutter area.

可选的，步骤103包括以下步骤：Optionally, step 103 includes the following steps:

A、应用AP聚类算法对所述监视区域内的n帧累积的回波进行聚类，获得聚类结果。A. Applying the AP clustering algorithm to cluster the echoes accumulated in n frames in the monitoring area to obtain a clustering result.

AP(Affinity Propagation，近似传播)聚类算法是根据N个数据点之间的相似度进行聚类，这些相似度可以是对称的，即两个数据点互相之间的相似度一样(如欧式距离)；也可以是不对称的，即两个数据点互相之间的相似度不等。这些相似度组成N×N的相似度矩阵S(其中N为有N个数据点)。AP算法不需要事先指定聚类数目，相反它将所有的数据点都作为潜在的聚类中心，称之为exemplar。The AP (Affinity Propagation, Approximate Propagation) clustering algorithm is based on the similarity between N data points, which can be symmetrical, that is, the similarity between two data points is the same (such as the Euclidean distance ); it can also be asymmetric, that is, the similarity between two data points is not equal to each other. These similarities form an N×N similarity matrix S (wherein N means that there are N data points). The AP algorithm does not need to specify the number of clusters in advance, instead it regards all data points as potential cluster centers, called exemplar.

AP聚类算法的主要思想是：以S矩阵的对角线上的数值S(k′，k′)作为k′点能否成为聚类中心的评判标准，这意味着该值越大，这个点成为聚类中心的可能性也就越大，这个值便是参考度P。聚类的数量受到参考度P的影响，如果认为每个数据点都有可能作为聚类中心，那么P就应取相同的值。如果取输入的相似度的均值作为P的值，得到聚类数量是中等的。如果取最小值，得到类数较少的聚类。AP聚类算法中传递两种类型的消息，即Responsibility和Availability。R(i′，k′)表示从点i′发送到候选聚类中心k′的数值消息，反映k′点是否适合作为i′点的聚类中心。A(i′，k′)则表示从候选聚类中心k′发送到i′的数值消息，反映i′点是否选择k′作为其聚类中心。R(i′，k′)与A(i′，k′)越强，则k′点作为聚类中心的可能性就越大，并且i′点隶属于以k′点为聚类中心的聚类的可能性也越大。AP聚类算法通过迭代过程不断更新每一个点的吸引度和归属度值，直到产生M个高质量的exemplar，同时将其余的数据点分配到相应的聚类中。The main idea of the AP clustering algorithm is: take the value S(k′, k′) on the diagonal of the S matrix as the criterion for judging whether point k′ can become a cluster center, which means that the larger the value, the The more likely the point is to become the cluster center, this value is the reference degree P. The number of clusters is affected by the reference degree P. If each data point is considered to be a cluster center, then P should take the same value. If the mean value of the input similarity is taken as the value of P, the number of clusters obtained is moderate. If the minimum value is taken, clusters with fewer classes are obtained. Two types of messages are passed in the AP clustering algorithm, namely Responsibility and Availability. R(i', k') represents the numerical message sent from point i' to candidate cluster center k', reflecting whether point k' is suitable as the cluster center of point i'. A(i', k') represents the numerical message sent from the candidate cluster center k' to i', reflecting whether point i' selects k' as its cluster center. The stronger R(i', k') and A(i', k') are, the more likely point k' is to be the cluster center, and point i' belongs to the cluster with point k' as the cluster center The probability of clustering is also greater. The AP clustering algorithm continuously updates the attractiveness and belonging value of each point through an iterative process until M high-quality exemplars are generated, and at the same time, the remaining data points are assigned to the corresponding clusters.

在这里需要说明的是：What needs to be explained here is:

1、exemplar：指的是聚类中心。1. exemplar: refers to the cluster center.

2、similarity：数据点i′和点j′的相似度记为S(i′，j′)。是指点j′作为点i′的聚类中心的相似度。一般使用欧氏距离来计算。在该算法中，所有点与点的相似度全部取为负值。这样，相似度值越大说明点与点的距离越近，便于后续计算。2. Similarity: The similarity between data point i' and point j' is recorded as S(i', j'). is the similarity of point j' as the cluster center of point i'. Generally, Euclidean distance is used to calculate. In this algorithm, all point-to-point similarities are taken as negative values. In this way, the larger the similarity value, the closer the distance between points is, which is convenient for subsequent calculations.

3、preference：数据点i′的参考度称为P(i′)或S(i′，i′)。是指点i′作为聚类中心的参考度。一般取S相似度值的中值。3. Preference: The reference degree of data point i' is called P(i') or S(i', i'). is the reference degree of point i′ as the cluster center. Generally, the median value of the S similarity value is taken.

4、Responsibility：R(i′，k′)用来描述点k′适合作为数据点i′的聚类中心的程度。4. Responsibility: R(i', k') is used to describe the degree to which point k' is suitable as the cluster center of data point i'.

Availability：A(i′，k′)用来描述点i′选择点k′作为其聚类中心的适合程度。Availability: A(i', k') is used to describe the suitability of point i' to select point k' as its cluster center.

如图2所示，为Responsibility(图2中记为R)和Availability(图2中记为A)两者之间的关系。As shown in Figure 2, it is the relationship between Responsibility (marked as R in Figure 2) and Availability (marked as A in Figure 2).

R(i′，k′)、A(i′，k′)与R(k′，k′)的计算公式为：The calculation formulas of R(i', k'), A(i', k') and R(k', k') are:

R(i′，k′)＝S(i′，k′)-max{A(i′，j′)+S(i′，j′)}(j′∈{1，2，…，N}&j′≠k′) (1)R(i',k')=S(i',k')-max{A(i',j')+S(i',j')}(j'∈{1,2,...,N }&j′≠k′) (1)

R(k′，k′)＝P(k′)-max{A(k′，j′)+S(k′，j′)}(j′∈{1，2，…，N}&j′≠k′) (3)R(k', k')=P(k')-max{A(k', j')+S(k', j')}(j'∈{1,2,...,N}&j' ≠k′) (3)

上面可以看出，当P(k′)较大使得R(k′，k′)较大时，A(i′，k′)也较大，从而类代表k′作为最终聚类中心的可能性较大；同样，当越多的P(i′)较大时，越多的类代表倾向于成为中心的聚类中心。因此，增大或减小P可以增加或减少AP输出的聚类数目。It can be seen from the above that when P(k′) is large so that R(k′, k′) is large, A(i′, k′) is also large, so that the class represents the possibility of k′ as the final clustering center similarly, when more P(i′) is larger, more classes represent cluster centers that tend to be centers. Therefore, increasing or decreasing P can increase or decrease the number of clusters output by AP.

5、Damping factor(阻尼系数)λ：主要是起收敛作用的。每次迭代，吸引度R_i′和归属度A_i′要与上一次的R_i′-1和A_i′-1进行加权更新。公式如下：5. Damping factor (damping factor) λ: mainly for convergence. At each iteration, the attractiveness R _i' and the belonging A _i' are weighted and updated with the last R _i'-1 and A _i'-1 . The formula is as follows:

R_i′＝(1-λ)*R_i′+λ*R_i′-1 (4)R _i' = (1-λ)*R _i' +λ*R _i'-1 (4)

A_i′＝(1-λ)*A_i′+λ*A_i′-1 (5)A _i' = (1-λ)*A _i' +λ*A _i'-1 (5)

其中，λ∈[0.5，1)。where, λ∈[0.5, 1).

AP聚类算法的具体步骤如下：The specific steps of the AP clustering algorithm are as follows:

(1)计算N个点之间的相似度值，将值放在S矩阵中，再选取P值(一般取S的中值)。(1) Calculate the similarity value between N points, put the value in the S matrix, and then select the P value (generally the median value of S).

(2)设置一个最大迭代次数，迭代过程开始后，计算每一次的R值和A值，根据R(k′，k′)+A(k′，k′)的值来判断是否为聚类中心，当迭代次数超过最大值或者当聚类中心连续多少次迭代不发生改变时终止计算，获得聚类结果。(2) Set a maximum number of iterations. After the iteration process starts, calculate the R value and A value each time, and judge whether it is a cluster according to the value of R(k', k')+A(k', k') Center, when the number of iterations exceeds the maximum value or when the number of consecutive iterations of the cluster center does not change, the calculation is terminated and the clustering result is obtained.

B、对于所述聚类结果中聚类点数大于η的区域，用Graham scan方法寻找凸包，确定杂波区。B. For the region where the number of clustering points is greater than n in the clustering result, use the Graham scan method to find the convex hull and determine the clutter area.

所述η为预设值，可以根据仿真场景而定。Said n is a preset value, which can be determined according to the simulation scenario.

其中，Graham scan算法是一种计算一组平面有限点的凸包的算法。其求解凸包的具体步骤如下：Among them, the Graham scan algorithm is an algorithm for calculating the convex hull of a set of finite points on a plane. The specific steps to solve the convex hull are as follows:

(1)在所有点中选取y坐标最小的一点H，当作基点。如果存在多个点的y坐标都为最小值，则选取x坐标最小的一点。坐标相同的点应排除。然后按照其它各点P和基点构成的向量＜H，P＞与x轴的夹角进行排序，夹角由大至小进行顺时针扫描，反之则进行逆时针扫描。实现中无需求得夹角，只需根据向量的内积公式求出向量的模即可。以图3为例，基点为H，根据夹角由小至大排序后依次为H，K，C，D，L，F，G，E，I，B，A，J。下面进行逆时针扫描。(1) Select a point H with the smallest y coordinate among all points as the base point. If there are multiple points whose y-coordinates are all minimum, select the point with the minimum x-coordinate. Points with the same coordinates should be excluded. Then sort according to the angle between the vector <H, P> formed by other points P and the base point and the x-axis. The angle is scanned clockwise from large to small, otherwise, it is scanned counterclockwise. In the implementation, there is no need to obtain the included angle, and it is only necessary to find the modulus of the vector according to the inner product formula of the vector. Take Figure 3 as an example, the base point is H, and the order of the included angles is H, K, C, D, L, F, G, E, I, B, A, J after sorting from small to large. Scan counterclockwise below.

(2)线段＜H，K＞一定在凸包上，接着加入C。假设线段＜H，C＞也在凸包上，因为就H，K，C三点而言，它们的凸包就是由此三点所组成。但是接下来加入D时会发现，线段＜K，D＞才会在凸包上，所以将线段＜K，C＞排除，C点不可能是凸包。即当加入一点时，必须考虑到前面的线段是否会出现在凸包上。(2) The line segment <H, K> must be on the convex hull, and then add C. Assume that the line segment <H, C> is also on the convex hull, because as far as the three points H, K, and C are concerned, their convex hull is composed of these three points. But when adding D next, it will be found that the line segment <K, D> will be on the convex hull, so the line segment <K, C> is excluded, and point C cannot be a convex hull. That is, when adding a point, we must consider whether the previous line segment will appear on the convex hull.

(3)从基点开始，凸包上每条相临的线段的旋转方向应该一致，并与扫描的方向相反。如果发现新加的点使得新线段与上线段的旋转方向发生变化，则可判定上一点必然不在凸包上。实现时可用向量叉积进行判断，设新加入的点为P_n+1，上一点为P_n，再上一点为P_n-1。顺时针扫描时，如果向量＜P_n-1，P_n＞与＜P_n，P_n+1＞的叉积为正(逆时针扫描判断是否为负)，则将上一点删除。删除过程需要回溯，将之前所有叉积符号相反的点都删除，然后将新点加入凸包。在图3中，加入K点时，由于线段＜H，K＞相对于＜H，C＞为顺时针旋转，所以C点不在凸包上，应该删除，保留K点。接着加入D点，由于线段＜K，D＞相对＜H，K＞为逆时针旋转，故D点保留。按照上述步骤进行扫描，直到点集中所有的点都遍历完成，即得到如图3所示的由点H、J、B、G、D以及K连线形成的凸包，确定出杂波区。(3) Starting from the base point, the rotation direction of each adjacent line segment on the convex hull should be consistent and opposite to the scanning direction. If it is found that the newly added point makes the rotation direction of the new line segment and the previous line segment change, it can be determined that the previous point must not be on the convex hull. The vector cross product can be used for judgment during implementation. Let the newly added point be P _n+1 , the previous point be P _n , and the next point be P _n-1 . When scanning clockwise, if the cross product of the vectors <P _n-1 , P _n > and <P _n , P _n+1 > is positive (scan counterclockwise to judge whether it is negative), then delete the previous point. The deletion process requires backtracking, deleting all previous points with opposite signs of the cross product, and then adding new points to the convex hull. In Figure 3, when point K is added, since the line segment <H, K> rotates clockwise relative to <H, C>, point C is not on the convex hull, so point C should be deleted and point K should be retained. Then add point D, because the line segment <K, D> rotates counterclockwise relative to <H, K>, so point D remains. Scan according to the above steps until all the points in the point set are traversed, that is, the convex hull formed by the connecting lines of points H, J, B, G, D and K is obtained as shown in Figure 3, and the clutter area is determined.

当然，寻找凸包，确定杂波区还可以有现有技术中的其他方法，在这里不再详述。Of course, there may be other methods in the prior art to find the convex hull and determine the clutter area, which will not be described in detail here.

1041、初始化，获得初始时刻的后验强度，所述初始时刻为0时刻。1041. Initialize. Obtain the posterior strength at the initial time, where the initial time is time 0.

在进入仿真场景进行PHD预测时，应当假设目标运动模型和观测模型都满足线性高斯条件；假设目标的存活概率和探测概率相互独立，且与目标状态无关；假设新生目标随机集和衍生目标随机集的强度均为高斯混合形式；用J₀个高斯成分的加权和进行初始化。When entering the simulation scene for PHD prediction, it should be assumed that both the target motion model and the observation model satisfy the linear Gaussian condition; the survival probability and detection probability of the target are assumed to be independent of each other and have nothing to do with the target state; the new target random set and the derived target random set are assumed The intensities of are all in Gaussian mixture form; initialized with the weighted sum of J ₀ Gaussian components.

假设目标运动模型和观测模型都满足线性高斯条件：Assume that both the target motion model and the observation model satisfy the linear Gaussian condition:

目标运动模型： f_k|k-1(x|ζ)＝N(x；F_k-1ζ，Q_k-1) (6)Target motion model: f _k|k-1 (x|ζ)=N(x; F _k-1 ζ, Q _k-1 ) (6)

观测模型： g_k(z|x)＝N(z；H_kx，R_k) (7)Observation model: g _k (z|x)=N(z; H _k x, R _k ) (7)

其中，N(·；m，P)表示其密度均值为m，方差为P的高斯分布；F_k-1为状态转移矩阵，Q_k-1为系统噪声协方差矩阵；H_k为观测矩阵，R_k为测量噪声协方差矩阵。Among them, N(·;m, P) represents the Gaussian distribution whose density mean is m and variance is P; F _k-1 is the state transition matrix, Q _k-1 is the system noise covariance matrix; H _k is the observation matrix, R _k is the measurement noise covariance matrix.

假设目标的存活概率和探测概率相互独立，且与目标状态无关：Assume that the survival probability and detection probability of the target are independent of each other and independent of the target state:

存活概率： p_S，k(x)＝p_S，k (8)Survival probability: p _{S, k} (x) = p _{S, k} (8)

探测概率： p_D，k(x)＝p_D，k (9)Probability of detection: p _{D, k} (x) = p _{D, k} (9)

假设新生目标随机集和衍生目标随机集的强度均为高斯混合形式：Assume that the intensities of both the nascent target random set and the derived target random set are Gaussian mixtures:

新生目标随机集的强度： $γ_{k} (x) = Σ_{i = 1}^{J_{γ, k}} ω_{γ, k}^{(i)} N (x; m_{γ, k}^{(i)}, P_{γ, k}^{(i)}) - - - (10)$ The strength of the freshman target random set: $γ_{k} (x) = Σ_{i = 1}^{J_{γ, k}} ω_{γ, k}^{(i)} N (x; m_{γ, k}^{(i)}, P_{γ, k}^{(i)}) - - - (10)$

衍生目标随机集的强度：The strength of the random set of derived targets:

${β β}_{k k | | k k - - 11} ((x x | | ζ ζ)) = = {Σ Σ}_{j j = = 11}^{{J J}_{β β,, k k}} {ω ω}_{β β,, k k}^{((j j))} N N ((x x;; {F f}_{β β,, k k - - 11}^{((j j))} ζ ζ + + {d d}_{β β,, k k - - 11}^{((j j))},, {Q Q}_{β β,, k k - - 11}^{((j j))})) - - - - - - ((1111))$

其中，J_γ，k，i＝1，2，...，J_γ，k等参数决定了新生目标随机集的强度，为新生目标随机集强度的第i个高斯成分的权值、均值和协方差，J_γ，k是k时刻新生目标高斯成分的个数。同理：J_β，k，j＝1，2，...，J_β，k等参数决定了由目标ζ衍生的目标随机集的强度，为新生目标随机集强度的第i个高斯成分的权值、均值和协方差，J_β，k是k时刻新生目标高斯成分的个数。在这里所述的各个参数的在初始时刻即0时刻时的值都是已知的。在后续时刻的值可以根据初始时刻即0时刻时的值推导计算出来。Among them, J _{γ, k} , i=1, 2, ..., J _{γ, k} and other parameters determine the strength of the random set of newborn targets, is the weight, mean and covariance of the i-th Gaussian component of the random set strength of the newborn target, J _{γ, and k} is the number of Gaussian components of the newborn target at time k. Similarly: J _{β, k} , j = 1, 2, ..., J _{β, k} and other parameters determine the strength of the target random set derived from the target ζ, is the weight, mean and covariance of the i-th Gaussian component of the random set strength of the newborn target, J _{β, k} is the number of Gaussian components of the newborn target at time k. The values of the various parameters described here at the initial time ie time 0 are all known. The value at the subsequent time can be derived and calculated based on the value at the initial time, ie, time 0.

用J₀个高斯成分的加权和进行初始化，获得初始时刻的后验强度：Initialize with a weighted sum of J ₀ Gaussian components to obtain the posterior strength at the initial moment:

${v v}_{00} = = {Σ Σ}_{i i = = 11}^{{J J}_{00}} {ω ω}_{00}^{((i i))} N N ((x x;; {m m}_{00}^{((i i))},, {P P}_{00}^{((i i))})) - - - - - - ((1212))$

其中N(·；m，P)是均值为m，方差为P的高斯分布。其权值的和是期望的初始目标数：公式(12)中的其他相关参数J₀，在用户应用本发明进行仿真时，由用户根据具体仿真环境设定自己设定。where N(·;m,P) is a Gaussian distribution with mean m and variance P. The sum of their weights is the desired number of initial targets: Other related parameters J ₀ in formula (12), When the user applies the present invention for simulation, the user sets his own settings according to the specific simulation environment.

1042、针对k时刻得到的回波值，选择未包含在所述杂波区的回波进行PHD预测计算，根据k-1时刻的后验强度，获得k时刻的预测强度。1042. For the echo value obtained at time k, select echoes not included in the clutter area to perform PHD prediction calculation, and obtain the predicted intensity at time k according to the posterior intensity at time k-1.

本发明实施例中的k时刻、k-1时刻为观测时刻，其中k-1时刻为k时刻的前一观测时刻。k-1＝0时刻为初始时刻，此时，k＝1。Time k and time k-1 in the embodiment of the present invention are observation times, wherein time k-1 is an observation time before time k. The time k-1=0 is the initial time, and at this time, k=1.

在逐时刻进行PHD预测时，对于每个观测时刻得到的回波值，若其中有回波在步骤103中的杂波区，则去掉那些回波，再进行PHD预测。When performing PHD prediction time by time, for the echo values obtained at each observation time, if there are echoes in the clutter area in step 103, those echoes are removed, and then PHD prediction is performed.

假定k-1时刻的后验强度具有如下混合高斯形式：Assume that the posterior strength at time k-1 has the following mixed Gaussian form:

${v v}_{k k - - 11} ((x x)) = = {Σ Σ}_{i i = = 11}^{{J J}_{k k - - 11}} {ω ω}_{k k - - 11}^{((i i))} N N ((x x;; {m m}_{k k - - 11}^{((i i))},, {P P}_{k k - - 11}^{((i i))})) - - - - - - ((1313))$

在k-1＝0时，公式(13)可以用公式(12)表示：When k-1=0, formula (13) can be represented by formula (12):

${v v}_{k k - - 11} ((x x)) = = {v v}_{00} = = {Σ Σ}_{i i = = 11}^{{J J}_{00}} {ω ω}_{00}^{((i i))} N N ((x x;; {m m}_{00}^{((i i))},, {P P}_{00}^{((i i))})) - - - - - - ((1313 a a))$

则针对k时刻得到的回波值，选择未包含在所述杂波区的回波进行PHD预测计算，根据k-1时刻的后验强度，获得k时刻的预测强度，k时刻的预测强度也为高斯混合形式，如下给出：Then, for the echo value obtained at time k, select echoes not included in the clutter area for PHD prediction calculation, and obtain the predicted strength at time k according to the posterior strength at time k-1, and the predicted strength at time k is also In the form of a Gaussian mixture, it is given as:

v_k|k-1(x)＝v_S，k|k-1(x)+v_β，k|k-1(x)+γ_k(x) (14)v _k|k-1 (x)=v _{S, k|k-1} (x)+v _{β, k|k-1} (x)+γ _k (x) (14)

其中，γ_k(x)在公式(10)中已给出。Among them, γ _k (x) has been given in formula (10).

${v v}_{S S,, k k | | k k - - 11} ((x x)) = = {p p}_{S S,, k k} {Σ Σ}_{j j = = 11}^{{J J}_{k k - - 11}} {ω ω}_{k k - - 11}^{((j j))} N N ((x x;; {m m}_{S S,, k k | | k k - - 11}^{((j j))},, {P P}_{S S,, k k | | k k - - 11}^{((j j))})) - - - - - - ((1515))$

${m m}_{S S,, k k | | k k - - 11}^{((j j))} = = {F f}_{k k - - 11} {m m}_{k k - - 11}^{((j j))} - - - - - - ((1616))$

${P P}_{S S,, k k | | k k - - 11}^{((j j))} = = {Q Q}_{k k - - 11} + + {F f}_{k k - - 11} {P P}_{k k - - 11}^{((j j))} {F f}_{k k - - 11}^{T T} - - - - - - ((1717))$

${v v}_{β β,, k k | | k k - - 11} ((x x)) = = {Σ Σ}_{j j = = 11}^{{J J}_{k k - - 11}} {Σ Σ}_{l l = = 11}^{{J J}_{β β,, k k}} {ω ω}_{k k - - 11}^{((j j))} {ω ω}_{β β,, k k}^{((l l))} N N ((x x;; {m m}_{β β,, k k | | k k - - 11}^{((j j,, l l))},, {P P}_{β β,, k k | | k k - - 11}^{((j j,, l l))})) - - - - - - ((1818))$

${m m}_{β β,, k k | | k k - - 11}^{((j j,, l l))} = = {F f}_{β β,, k k - - 11}^{((l l))} {m m}_{k k | | k k - - 11}^{((j j))} + + {d d}_{β β,, k k - - 11}^{((l l))} - - - - - - ((1919))$

${P P}_{β β,, k k | | k k - - 11}^{((j j,, l l))} = = {Q Q}_{β β,, k k - - 11}^{((l l))} + + {F f}_{β β,, k k - - 11}^{((l l))} {P P}_{β β,, k k - - 11}^{((j j))} {(({F f}_{β β,, k k - - 11}^{((l l))}))}^{T T} - - - - - - ((2020))$

由公式(13)中含有的参数以及J_k-1可以计算出公式(10)以及公式(15)-(20)，然后通过公式(10)以及公式(15)-(20)就可以计算出公式(14)中k时刻的预测强度v_k|k-1(x)。By the parameters contained in formula (13) And J _k-1 can calculate the formula (10) and the formula (15)-(20), and then can calculate the k time forecast in the formula (14) through the formula (10) and the formula (15)-(20) Intensity v _k|k-1 (x).

在这里，可以将由公式(10)，以及(15)-(20)，计算得到k时刻的预测强度记为如下高斯混合形式：Here, the predicted strength at time k calculated by formula (10) and (15)-(20) can be recorded as the following Gaussian mixture form:

${v v}_{k k | | k k - - 11} ((x x)) = = {Σ Σ}_{i i = = 11}^{{J J}_{k k | | k k - - 11}} {ω ω}_{k k | | k k - - 11}^{((i i))} N N ((x x;; {m m}_{k k | | k k - - 11}^{((i i))},, {P P}_{k k | | k k - - 11}^{((i i))})) - - - - - - ((21 twenty one))$

105、根据所述k时刻的预测强度，用更新的回波进行PHD更新，获得k时刻的后验强度，其中，所述更新的回波包括k时刻所得的回波值中所述未包含在所述杂波区的回波，以及k时刻真实目标的量测刚好落入所述杂波区的回波。105. According to the predicted intensity at time k, use the updated echo to perform PHD update to obtain the posterior intensity at time k, wherein the updated echo includes the echo values not included in the echo value obtained at time k. The echo of the clutter area and the measurement of the real target at time k just fall into the echo of the clutter area.

所述后验强度为高斯混合形式，包括各个高斯成分的权值、均值和协方差。The posterior strength is in the form of a Gaussian mixture, including the weight, mean and covariance of each Gaussian component.

在进行PHD更新时，需要对杂波区中的回波进行自适应判定，若确定杂波区里的都是杂波，而无真实目标的量测，则杂波区内的回波均不用；但若发现真实目标的量测刚好落入杂波区中，则将其取出放入PHD更新要用的回波中，改变杂波区的大小。判断是否落入杂波区的方法为：用目标状态计算量测，将得出的量测与杂波区的回波比对，若发现量测落入杂波区，则将量测取出，放入PHD更新需用的回波中；反之，则进行PHD更新时用不包括杂波区的回波。When performing PHD update, it is necessary to make an adaptive judgment on the echoes in the clutter area. If it is determined that all the clutter areas are clutter and there is no real target measurement, the echoes in the clutter area will not be used. ; But if it is found that the measurement of the real target just falls into the clutter area, then take it out and put it into the echo used for PHD update, and change the size of the clutter area. The method of judging whether it falls into the clutter area is: use the target state to calculate the measurement, and compare the obtained measurement with the echo in the clutter area. If it is found that the measurement falls into the clutter area, take the measurement out. Put it into the echoes needed for PHD update; otherwise, use the echoes that do not include the clutter area when performing PHD update.

由步骤1042可以获得k时刻的预测强度为如下高斯混合形式：From step 1042, the predicted strength at time k can be obtained as the following Gaussian mixture form:

则k时刻的后验强度也为高斯混合形式：Then the posterior intensity at time k is also in the form of Gaussian mixture:

${v v}_{k k} ((x x)) = = ((11 - - {p p}_{D D.,, k k})) {v v}_{k k | | k k - - 11} ((x x)) + + \underset{z z &Element; &Element; {Z Z}_{k k}}{Σ Σ} {v v}_{D D.,, k k} ((x x;; z z)) - - - - - - ((22 twenty two))$

其中，in,

${v v}_{D D.,, k k} ((x x;; z z)) = = {Σ Σ}_{j j = = 11}^{{J J}_{k k | | k k - - 11}} {ω ω}_{k k}^{((j j))} ((z z)) N N ((x x;; {m m}_{k k | | k k}^{((j j))} ((z z)),, {P P}_{k k | | k k}^{((j j))})) - - - - - - ((23 twenty three))$

${ω ω}_{k k}^{((j j))} ((z z)) = = \frac{{p p}_{D D.,, k k} {ω ω}_{k k | | k k - - 11}^{((j j))} {q q}_{k k}^{((j j))} ((z z))}{{k k}_{k k} ((z z)) + + {p p}_{D D.,, k k} {Σ Σ}_{l l = = 11}^{{J J}_{k k | | k k - - 11}} {{ω ω}_{k k | | k k - - 11}^{((l l))} q q}_{k k}^{((l l))} ((z z))} - - - - - - ((24 twenty four))$

${m m}_{k k | | k k}^{((j j))} ((z z)) = = {m m}_{k k | | k k - - 11}^{((j j))} + + {K K}_{k k}^{((j j))} ((z z - - {H h}_{k k} {m m}_{k k | | k k - - 11}^{((j j))})) - - - - - - ((2525))$

${P P}_{k k | | k k}^{((j j))} = = [[I I - - {K K}_{k k}^{((j j))} {H h}_{k k}]] {P P}_{k k | | k k - - 11}^{((j j))} - - - - - - ((2626))$

${K K}_{k k}^{((j j))} = = {P P}_{k k | | k k - - 11}^{((j j))} {H h}_{k k}^{T T} {(({H h}_{k k} {P P}_{k k | | k k - - 11}^{((j j))} {H h}_{k k}^{T T} + + {R R}_{k k}))}^{- - 11} - - - - - - ((2727))$

由公式(21)中含有的参数以及J_k|k-1可以计算出公式(23)-(27)，然后通过公式(23)-(27)就可以计算出公式(22)中k时刻的后验强度v_k(x)。By the parameters contained in formula (21) And J _k|k-1 can calculate the formula (23)-(27), and then the posterior strength v _k (x) at time k in the formula (22) can be calculated through the formula (23)-(27).

106、将所述后验强度中权值小于T的高斯成分滤除，将高斯成分间的距离小于U的高斯成分进行合并，获得目标强度；其中，U≤10，T≤10^-5。106. Filter out Gaussian components whose weights are smaller than T in the posterior strength, and combine Gaussian components whose distance between Gaussian components is smaller than U to obtain a target strength; wherein, U≤10, T≤10 ^-5 .

设定一个阈值T(T≤10^-5)，当一些高斯成分的权值小于T时，将这些高斯成分滤除掉，具体公式可以表示如下：Set a threshold T (T≤10 ^-5 ), when the weight of some Gaussian components is smaller than T, filter out these Gaussian components, the specific formula can be expressed as follows:

${\overset{- -}{v v}}_{k k} ((x x)) = = \frac{{Σ Σ}_{l l = = 11}^{{J J}_{k k}} {ω ω}_{k k}^{((l l))}}{{Σ Σ}_{j j = = {N N}_{p p} + + 11}^{{J J}_{k k}} {ω ω}_{k k}^{((j j))}} {Σ Σ}_{i i = = {N N}_{p p} + + 11}^{{J J}_{k k}} {ω ω}_{k k}^{((i i))} N N ((x x;; {m m}_{k k}^{((i i))},, {P P}_{k k}^{((i i))})) - - - - - - ((2828))$

其中，i＝N_p+1，...，J_k是大于阈值的权值。in, i=N _p +1, . . . , J _k are weights greater than the threshold.

设定一个阈值U(U＜10)，当一些高斯成分间的距离小于U时，将这些高斯成分进行合并。给出一个最大允许的高斯成分数目J_max。A threshold U (U<10) is set, and when the distance between some Gaussian components is smaller than U, these Gaussian components are combined. Gives a maximum allowed number of Gaussian components J _max .

令l＝0， $I = {i = 1,2, . . ., J_{k} | ω_{k}^{(i)} > T} .$ let l=0, $I = {i = 1,2, . . ., J_{k} | ω_{k}^{(i)} > T} .$

重复repeat

l：＝l+1 (29)l:=l+1 (29)

$j j = = arg arg \underset{i i &Element; &Element; I I}{max max} {ω ω}_{k k}^{((i i))} - - - - - - ((3030))$

$L L : : = = {{i i &Element; &Element; I I | | {(({m m}_{k k}^{((i i))} - - {m m}_{k k}^{((j j))}))}^{T T} {(({P P}_{k k}^{((i i))}))}^{- - 11} (({m m}_{k k}^{((i i))} - - {m m}_{k k}^{((j j))})) \leq \leq U u}} - - - - - - ((3131))$

${\overset{~ ~}{ω ω}}_{k k}^{((l l))} = = \underset{i i &Element; &Element; L L}{Σ Σ} {\overset{~ ~}{ω ω}}_{k k}^{((i i))} - - - - - - ((3232))$

${\overset{~ ~}{m m}}_{k k}^{((l l))} = = \frac{11}{{\overset{~ ~}{ω ω}}_{k k}^{((l l))}} \underset{i i &Element; &Element; L L}{Σ Σ} {ω ω}_{k k}^{((i i))} {x x}_{k k}^{((i i))} - - - - - - ((3333))$

${\overset{~ ~}{P P}}_{k k}^{((l l))} = = \frac{11}{{\overset{~ ~}{ω ω}}_{k k}^{((l l))}} \underset{i i &Element; &Element; L L}{Σ Σ} {ω ω}_{k k}^{((i i))} (({P P}_{k k}^{((i i))} + + (({\overset{~ ~}{m m}}_{k k}^{((l l))} - - {m m}_{k k}^{((i i))})) {(({\overset{~ ~}{m m}}_{k k}^{((l l))} - - {m m}_{k k}^{((i i))}))}^{T T})) - - - - - - ((3434))$

I：＝I\L (35)I:＝I\L (35)

直到 until

如果l＞J_max，用J_max个权值较大的高斯成分代替得到的即为合并后的高斯成分。If l＞J _max , use J _max Gaussian components with larger weights instead owned is the combined Gaussian component.

其具体合并过程本领域技术人员都清楚了解，在此不再详述。The specific merging process is well understood by those skilled in the art and will not be described in detail here.

进行高斯合并后的所述目标强度也是高斯混合形式，可以如下所示：The target intensity after Gaussian combination is also a Gaussian mixture form, which can be shown as follows:

$\overset{~ ~}{v v} ((x x)) = = {Σ Σ}_{i i = = 11}^{{J J}_{max max}} {\overset{~ ~}{ω ω}}_{k k}^{((i i))} N N ((x x;; {\overset{~ ~}{m m}}_{k k}^{((i i))},, {\overset{~ ~}{P P}}_{k k}^{((i i))})) - - - - - - ((3636))$

107、计算所述目标强度中所有高斯成分的权值和，得到所述监视区域内的估计目标数；提取所述目标强度中权值大于τ的高斯成分，作为估计的目标状态，其中，所述τ≥0.5。107. Calculate the sum of the weights of all Gaussian components in the target intensity to obtain the estimated number of targets in the monitoring area; extract the Gaussian components with a weight greater than τ in the target intensity as the estimated target state, wherein the Said τ≥0.5.

获取估计目标数：通过混合高斯项的权值和获取监视区域内的目标数，即Obtain the estimated number of targets: by mixing the weight of the Gaussian term and obtaining the number of targets in the monitoring area, that is

${\overset{^^}{n no}}_{k k} = = {Σ Σ}_{i i = = 11}^{{J J}_{max max}} {\overset{~ ~}{ω ω}}_{k k}^{((i i))} - - - - - - ((3737))$

获取目标状态：提取权值大于τ(τ≥0.5)的高斯成分，即Obtain the target state: extract the Gaussian components whose weight is greater than τ (τ≥0.5), namely

${\overset{^^}{X x}}_{k k} = = {{{\overset{~ ~}{m m}}_{k k}^{((i i))} : : {\overset{~ ~}{ω ω}}_{k k}^{((i i))} > > τ τ}} - - - - - - ((3838))$

为了验证本发明方案，下面给出面向不均匀杂波环境下多目标跟踪的自适应概率假设密度方法仿真结果。In order to verify the solution of the present invention, the simulation results of the adaptive probability hypothesis density method for multi-target tracking in an uneven clutter environment are given below.

仿真场景：考虑不均匀杂波环境下监视区域内有三个目标的场景。Simulation scenario: Consider a scenario where there are three targets in the surveillance area in an uneven clutter environment.

目标1在k＝1时刻出现，初始状态为[-25m，4m/s，10m，2m/s]^T，在k＝38时刻消失；目标2在k＝10时刻出现，初始状态为[550m，-15m/s，100m，4m/s]^T，在k＝60时刻消失；目标3在k＝22时刻出现，初始状态为[-270m，18m/s，-190m，6m/s]^T，在k＝51时刻消失。Target 1 appears at time k=1, with an initial state of [-25m, 4m/s, 10m, 2m/s] ^T , and disappears at time k=38; target 2 appears at time k=10, with an initial state of [550m, -15m/s, 100m, 4m/s] ^T , disappear at time k=60; target 3 appears at time k=22, the initial state is [-270m, 18m/s, -190m, 6m/s] ^T , at k=51 time disappears.

目标状态方程为：The target state equation is:

X_k+1＝F_kX_k+v_k (39)X _k+1 = F _k X _k +v _k (39)

其中，状态向量X_k＝[x_1，k，x_2，k，x_3，k，x_4，k]^T，[x_1，k，x_3，k]^T是k时刻的目标位置，[x_2，k，x_4，k]^T是k时刻的目标速度。采样间隔T＝1s。状态转移矩阵F_k、过程噪声的协方差分别是 $F_{k} = [\begin{matrix} 1 & T & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & T \\ 0 & 0 & 0 & 1 \end{matrix}], Q = [\begin{matrix} \frac{T^{4}}{4} & 0 & \frac{T^{3}}{2} & 0 \\ 0 & T^{2} & 0 & \frac{T^{3}}{2} \\ \frac{T^{3}}{2} & 0 & \frac{T^{4}}{4} & 0 \\ 0 & \frac{T^{3}}{2} & 0 & T^{2} \end{matrix}] .$ 式中，过程噪声是零均值的高斯白噪声，且和量测噪声序列相互独立。目标存活概率为p_S＝0.975。Among them, the state vector X _k =[x _{1, k} , x _{2, k} , x _{3, k} , x _{4, k} ] ^T , [x _{1, k} , x _{3, k} ] ^T is the target position at time k, [ x _{2, k} , x _{4, k} ] ^T is the target velocity at time k. Sampling interval T=1s. The state transition matrix F _k and the covariance of the process noise are respectively $f_{k} = [\begin{matrix} 1 & T & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & T \\ 0 & 0 & 0 & 1 \end{matrix}], Q = [\begin{matrix} \frac{T^{4}}{4} & 0 & \frac{T^{3}}{2} & 0 \\ 0 & T^{2} & 0 & \frac{T^{3}}{2} \\ \frac{T^{3}}{2} & 0 & \frac{T^{4}}{4} & 0 \\ 0 & \frac{T^{3}}{2} & 0 & T^{2} \end{matrix}] .$ In the formula, the process noise is Gaussian white noise with zero mean value, and it is independent of the measurement noise sequence. The target probability of survival is p _S =0.975.

量测方程为：The measurement equation is:

Z_k＝H_kX_k+ω_k (40)Z _k ＝H _k X _k +ω _k (40)

式中，量测噪声是零均值的白噪声。量测矩阵和量测噪声协方差分别是 $H_{k} = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{matrix}], R = [\begin{matrix} 100 & 0 \\ 0 & 100 \end{matrix}] .$ 目标探测概率为p_D＝0.96。where the measurement noise is white noise with zero mean. The measurement matrix and the measurement noise covariance are respectively $h_{k} = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{matrix}], R = [\begin{matrix} 100 & 0 \\ 0 & 100 \end{matrix}] .$ The target detection probability is p _D =0.96.

杂波的空间分布是非均匀的，在某几个区域杂波很密集。整个监视空域是[-400m，600m]×[-400m，800m]，在整个区域每个观测时刻的平均杂波数是5。密集杂波区域为：在[-200m，-150m]×[-200m，-150m]区域，每个观测时刻的平均杂波数是20；在[50m，100m]×[0m，100m]区域，每个观测时刻的平均杂波数是20；在[300m，350m]×[100m，150m]区域，每个观测时刻的平均杂波数是20。The spatial distribution of clutter is non-uniform, and clutter is very dense in some areas. The entire monitoring airspace is [-400m, 600m]×[-400m, 800m], and the average clutter number at each observation moment in the entire area is 5. The dense clutter area is: in the [-200m, -150m]×[-200m,-150m] area, the average clutter number at each observation moment is 20; in the [50m, 100m]×[0m, 100m] area, each The average clutter number of each observation moment is 20; in the [300m, 350m]×[100m, 150m] area, the average clutter number of each observation moment is 20.

仿真场景见图4。The simulation scene is shown in Figure 4.

仿真时，n＝5，ε＝200，η＝80。During simulation, n=5, ε=200, η=80.

在仿真实验中，本发明方法与传统的高斯混合PHD(简称为：GM-PHD)方法(标准PHD的一种实现方式)进行对比。本发明仿真结果见图5及图6，高斯混合PHD的仿真结果见图7及图8。In the simulation experiment, the method of the present invention is compared with the traditional Gaussian mixture PHD (referred to as: GM-PHD) method (an implementation of standard PHD). The simulation results of the present invention are shown in Fig. 5 and Fig. 6, and the simulation results of the Gaussian mixture PHD are shown in Fig. 7 and Fig. 8 .

对比两种方法，可以看出，在不均匀杂波环境下，即在既有稀疏杂波区又有密集杂波区的情况下，本发明(图5和图6显示)可以更好的估计每个时刻的目标状态及目标数，从图中可看出，相比高斯混合PHD，本发明估计每个时刻目标状态及目标数的精度有显著的提高，达到了提升目标跟踪性能的目的。Comparing the two methods, it can be seen that the present invention (shown in Figures 5 and 6) can better estimate The target state and the number of targets at each moment can be seen from the figure. Compared with the Gaussian mixture PHD, the accuracy of the present invention to estimate the target state and the number of targets at each moment is significantly improved, and the purpose of improving the target tracking performance is achieved.

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。The above is only a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Anyone skilled in the art can easily think of changes or substitutions within the technical scope disclosed in the present invention. Should be covered within the protection scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. a probability hypothesis density method for multiple target tracking under inhomogeneous clutter environment, is characterized in that, comprising:

Whether the number of echoes that judges n frame accumulation in monitor area is greater than threshold value ε, wherein, described n, ε is default value;

If not, directly carry out standard P HD prediction and PHD and upgrade, obtain number of targets and the dbjective state estimated;

If so, find convex closure, determine clutter district;

Initialization, the posteriority intensity of acquisition initial time, described initial time was 0 moment;

For the Echo Rating that the k moment records, select the echo that is not included in described clutter district to carry out PHD prediction and calculation, according to the posteriority intensity in k-1 moment, obtain the predicted intensity in k moment, wherein, described k is more than or equal to 1 integer;

According to the predicted intensity in described k moment, with upgrade echo carry out PHD renewal, obtain the posteriority intensity in k moment, wherein, the echo of described renewal comprises the echo that is not included in described clutter district described in the Echo Rating that the k moment records, and the measurement of k moment real goal just falls into the echo in described clutter district; Described posteriority intensity is Gaussian Mixture form, comprises weights, average and the covariance of each gauss component;

Weights in described posteriority intensity are less than to the gauss component filtering of T, the gauss component that the distance between gauss component is less than to U merges, and obtains target strength; Wherein, U≤10, T≤10 ^-5;

Calculate all gauss components in described target strength weights and, obtain the estimating target number in described monitor area; Extract the gauss component that weights in described target strength are greater than τ, as the dbjective state of estimating, wherein, described τ >=0.5.

2. method according to claim 1, is characterized in that, described searching convex closure, determines clutter district, comprising:

Approximate (AP) clustering algorithm of propagating of application carries out cluster to the echo of the n frame accumulation in described monitor area, obtains cluster result;

Count and be greater than the region of η for cluster in described cluster result, find convex closure by Graham scan method, determine clutter district, described η determines according to simulating scenes.