CN102322861A - 一种航迹融合方法 - Google Patents

Abstract

本发明属于多源信息融合技术领域，公开了一种航迹融合方法。包括：利用多传感器的观测信息建立数据间的相对距离矩阵；计算支持度函数，获得支持度矩阵，建立方程组，求解加权因子；将权因子与其相应的观测值相乘，分别滤波获取对应的滤波值，将获取的全部滤波值相加得到观测系数化后的滤波融合值；将滤波融合值作为状态更新输入值，利用Kalman滤波，分步对目标状态估计值进行估计更新。本发明通过对多传感器的观测信息进行观测系数化的滤波融合，在不增加数据处理复杂度的条件下，降低了观测信息不确定性对航迹融合的影响，同时在观测系数化滤波融合过程，考虑了观测信息的相关性，进而提高了观测精度，获得了对目标的可靠跟踪。

Description

一种航迹融合方法

技术领域

本发明属于多源信息融合技术领域，尤其涉及一种航迹融合方法。

背景技术

面对日益复杂的信息环境及信息表现形式的多样性、信息数量的巨大性、信息之间关系的复杂性，以及要求对信息处理的实时性，不仅超出了人脑的综合处理能力，而且单个传感器的量测信息存在着不精确性、不可靠性、不一致性及不完整性，因此必须依靠多传感器协同工作，利用多传感器信息融合技术综合优化处理信息，以最大限度地获取观测目标的状态、特征等完整信息，实现目标的精确定位及跟踪。信息融合是一门信息综合处理技术，即利用多传感器或多源数据的互补性、冗余性，使各传感器提供的各种信息能接近无损失和无误地加以融合，从而提高数据的质量。目标跟踪作为信息融合的一个经典应用及各个领域中的关键部分，引起了广泛的兴趣。针对此问题的研究，相继提出了一系列的航迹融合方法，不断地改善航迹融合质量，提高了目标的跟踪性能。

目前，针对航迹融合问题的方法主要分为两类，一类是集中式航迹融合方法，另一类是分布式航迹融合方法。前者由于信息损失最小，使其融合精度较高，然而当其中某个传感器失效时会对融合性能产生较大的影响，稳定性不高，且计算量大；后者虽融合精度相对较低，但结构简单，通信量小，同时具有局部跟踪能力。通过对分布式融合系统的研究，已提出的融合方法有简单方差凸组合法、加权融合法、基于分步式滤波的航迹融合法等。其中简单方差凸组合法是最早提出的融合估计方法，该方法在假设各局部状态估计误差统计独立的情况下实现的，计算比较简单；加权融合法是利用权因子简单组合多源信息，易实现，在工程上应用广泛，此外，针对加权因子的制定，现在已提出了许多分配权值的方法，如极大似然法、最小二乘法等，通过对制定权因子算法的改进，易达到对加权融合法性能的提高；相对于以上方法，基于分步式滤波的航迹融合法，主要是利用多传感器的观测值，结合系统先前信息对目标状态的预测估计值，通过Kalman滤波依次对目标状态估计值进行更新，从而得到基于全局系统的融合估计值。该方法虽然结构相对复杂，计算量较大，但因充分考虑了各个传感器所提供的观测信息，使其与其他算法相比，具有相对较高的信息融合精度。然而由于受到传感器本身性能及外界随机干扰的影响，使各传感器提供的观测信息具有不确定性及相关性，从而造成该融合算法稳健性及融合性能的降低。

发明内容

本发明的目的是为了解决现有的分步式滤波的航迹融合方法由于各传感器提供的观测信息的不确定性及相关性造成融合性能降低的问题，提出了一种航迹融合方法。

本发明的技术方案是：一种航迹融合方法，包括如下步骤：

S1.利用预处理后的多传感器的观测信息建立数据间的相对距离矩阵；

S2.依据步骤S1得到的相对距离矩阵，计算支持度函数，获得支持度矩阵，建立方程组，求解加权因子；

S3.将步骤S2得到的各个传感器的权因子与其相应的观测值相乘，分别滤波获取对应的滤波值，将获取的全部滤波值相加得到观测系数化后的滤波融合值；

S4.将步骤S3得到的滤波融合值作为状态更新输入值，利用Kalman滤波，分步对目标状态估计值进行估计更新。

所述步骤S1具体过程如下：计算数据间相对距离，建立相对距离矩阵，在k时刻，利用多传感器对同一目标的观测值，计算相对距离d_ij(k)：

d_ij(k)＝|z_i(k)-z_j(k)|i，j＝1，2，L，N，其中，z_i(k)表示预处理后的第i个传感器的观测信息，N是传感器的个数；

在k时刻N个传感器的相对距离矩阵D(k)记为：

$D\left(k\right)={\left[\begin{array}{cccc}{d}_{11}\left(k\right)& d_{12}\left(k\right)& L& d_{1N}\left(k\right)\\ d_{21}\left(k\right)& d_{22}\left(k\right)& L& d_{2N}\left(k\right)\\ L& L& L& L\\ d_{N1}\left(k\right)& d_{N2}\left(k\right)& L& d_{\mathrm{NN}}\left(k\right)\end{array}\right]}_{N\×N}.$

所述步骤S2具体包括如下分步骤：

S21.计算相对距离矩阵D(k)中元素的最大值D_Max(k)；

S22.计算支持度函数r_ij(k)，i，j＝1，2，L，N；

S23.建立支持度矩阵，N个传感器数据之间的支持度矩阵Sup_Matrix(k)为：

$\mathrm{Sup}\_\mathrm{Matrix}\left(k\right)={\left[\begin{array}{cccc}{r}_{11}\left(k\right)& r_{12}\left(k\right)& L& r_{1N}\left(k\right)\\ r_{21}\left(k\right)& r_{22}\left(k\right)& L& r_{2N}\left(k\right)\\ L& L& L& L\\ r_{N1}\left(k\right)& r_{N2}\left(k\right)& L& r_{\mathrm{NN}}\left(k\right)\end{array}\right]}_{N\×N};$

S24.计算权因子，具体过程如下：存在一组非负数：l₁，l₂，…，l_N，使得权因子为：w_i＝l₁r_i1+l₂r_i2+L+l_Nr_iN，i＝1，2，L，N，并且满足

，其矩阵形式为：W＝Sup_Matrix(k)*L，其中，W＝[w₁ w₂ L w_N]^T，L＝[l₁ l₂ L l_N]^T，建立矩阵方程λL＝Sup_Matrix(k)*L，计算特征值λ及相应的特征向量，取矩阵最大模特征值max(|λ|)相对应的特征向量，得到一组非负系数l₁，l₂，…，l_N，进而求得一个观测值对所有数据的总体支持度，即获得第i个传感器观测值的权因子可为：

i＝1，2，L，N。

步骤S4具体包括如下步骤：

S41.状态初始化，当k＝1时，利用等权因子加权思想求得全局状态估计值

及估计误差协方差P(1|1)分别为：

$\hat{X}\left(1\right|1)=\frac{{\hat{X}}_1\left(1\right|1)+{\hat{X}}_2\left(1\right|1)+L+{\hat{X}}_N\left(1\right|1)}{N}$

$P\left(1\right|1)=\frac{P_1\left(1\right|1)+P_2\left(1\right|1)+L+P_N\left(1\right|1)}{N^2}$

其中，

及P_i(1|1)分别为第i个传感器的局部估计和相应的估计误差协方差初始值；

S42.预测值求解，当k＝2，3，L，t/T，其中，t为观测时间，T为融合中心采样周期，基于全局系统的状态估计值

和估计误差协方差P(k-1|k-1)的一步预测估计值

及相应的预测误差协方差P(k|k-1)分别为：P(k|k-1)＝FP(k-1|k-1)F^T+GQG^T，其中，F为状态转移矩阵、G为过程噪声输入矩阵、Q为系统噪声协方差矩阵，*^T表示矩阵*的转置运算；

S43.分步式滤波更新，利用步骤S3得到的观测系数化后的滤波融合值，记为Z_fil(k|k)，对一步预测估计值

及相应预测误差协方差P(k|k-1)依次更新，具体实现过程如下：

S431.令

P₁(k|k-1)＝P(k|k-1)，根据Kalman滤波公式有：

${\hat{X}}_1\left(k\right|k)=F\hat{X}(k-1|k-1)+K_1\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_1\left(k\right){\hat{X}}_1\left(k\right|k-1)]$

P₁(k|k)＝[I-K₁(k)H₁(k)]P₁(k|k-1)

其中，H₁(k)、R₁(k)分别为第一个传感器的量测矩阵及观测噪声协方差，I为单位矩阵，K₁(k)为第一个传感器的滤波增益矩阵，即有：

$K_1\left(k\right)=P_1\left(k\right|k-1)H_1^T\left(k\right){[H_1\left(k\right)P_1\left(k\right|k-1)H_1^T\left(k\right)+R_1\left(k\right)]}^{-1};$

S432.令

P₂(k|k-1)＝P₁(k|k)，则根据Kalman滤波公式有：

${\hat{X}}_2\left(k\right|k)=\widehat{X_1}\left(k\right|k)+K_2\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_2\left(k\right){\hat{X}}_2\left(k\right|k-1)]$

$=F\hat{X}(k-1|k-1)+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{K}_i\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_i\left(k\right){\hat{X}}_i\left(k\right|k-1)]$

$P_2\left(k\right|k)=[I-K_2\left(k\right)H_2\left(k\right)]P_2\left(k\right|k-1)=\underset{i=1}{\overset{2}{\mathrm{\Π}}}[I-K_i\left(k\right)H_i\left(k\right)]P_1\left(k\right|k-1)$

其中，H₂(k)、R₂(k)为第二个传感器的量测矩阵及观测噪声协方差，K₂(k)为第二个传感器的滤波增益矩阵，即有：

S433.依次可以得到第i个传感器的状态更新估计值

及相应的误差协方差P_i(k|k)，i＝3，……，N，第N个传感器的状态更新估计值

及相应的误差协方差P_N(k|k)为：

${\hat{X}}_N\left(k\right|k)=\widehat{X_{N-1}}\left(k\right|k)+K_N\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_N\left(k\right){\hat{X}}_N\left(k\right|k-1)]$

$=F\hat{X}(k-1|k-1)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{K}_i\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_i\left(k\right){\hat{X}}_i\left(k\right|k-1)]$

$P_N\left(k\right|k)=[I-K_N\left(k\right)H_N\left(k\right)]P_N\left(k\right|k-1)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}[I-K_i\left(k\right)H_i\left(k\right)]P_N\left(k\right|k-1)$

其中，H_N(k)、R_N(k)为第N个传感器的量测矩阵及观测噪声协方差，K_N(k)为第N个传感器的滤波增益矩阵，即有：

S44.计算全局融合估计值，可求得k时刻基于全局系统的融合估计值

及相应的估计误差协方差P(k|k)分别为：

本发明的有益效果：本发明的方法通过对多传感器的观测信息进行观测系数化的滤波融合，在不增加数据处理复杂度的条件下，降低了观测信息不确定性对航迹融合的影响，同时在观测系数化滤波融合过程，考虑了观测信息的相关性，进而提高了观测精度，获得了对目标的可靠跟踪。

附图说明

图1为本发明的航迹融合方法流程图。

图2为本发明的观测系数化滤波融合模型示意图。

图3为本发明的观测系数化滤波融合过程基本流程图。

图4为本发明的分步式滤波融合基本流程图。

图5为本发明的分步式滤波更新流程图。

图6为实施例的多个传感器的融合跟踪与单个传感器下的局部跟踪的误差曲线图。

图7为本发明的方法与基于分步式滤波的航迹融合算法的误差曲线比较示意图。

图8为本发明实施例的观测系数化的滤波融合与单个传感器下观测值的误差曲线图。

具体实施方式

下面结合附图和具体的实施例对本发明做进一步的说明。

如图1所示，本发明的航迹融合方法包括如下步骤：

S1.利用预处理后的多传感器的观测信息建立数据间的相对距离矩阵；

S2.依据步骤S1得到的相对距离矩阵，计算支持度函数，获得支持度矩阵，建立方程组，求解加权因子；

S3.将步骤S2得到的各个传感器的权因子与其相应的观测值相乘，分别滤波获取对应的滤波值，将获取的全部滤波值相加得到观测系数化后的滤波融合值；

S4.将步骤S3得到的滤波融合值作为状态更新输入值，利用Kalman滤波，分步对目标状态估计值进行估计更新。

本发明的方法主要有两部分组成，一是观测系数化的滤波融合，二是分步式滤波融合，其中，步骤S2和S3构成了观测系数化的滤波融合，步骤S4实现了分步式滤波融合。步骤S1主要是利用预处理后的多传感器的观测信息获得相对距离矩阵，具体过程如下：计算数据间相对距离，建立相对距离矩阵，在k时刻，利用多传感器对同一目标的观测值，计算相对距离d_ij(k)：d_ij(k)＝|z_i(k)-z_j(k)|i，j＝1，2，L，N，其中，z_i(k)表示预处理后第i个传感器的观测信息，N是传感器的个数。可以看出，d_ij(k)越大，则两个传感器观测数据间的差别越大，即知它们的相互支持度越小。

进而可以得到k时刻N个传感器的相对距离矩阵D(k)：

$D\left(k\right)={\left[\begin{array}{cccc}{d}_{11}\left(k\right)& d_{12}\left(k\right)& L& d_{1N}\left(k\right)\\ d_{21}\left(k\right)& d_{22}\left(k\right)& L& d_{2N}\left(k\right)\\ L& L& L& L\\ d_{N1}\left(k\right)& d_{N2}\left(k\right)& L& d_{\mathrm{NN}}\left(k\right)\end{array}\right]}_{N\×N}.$

步骤S2和S3是观测系数化的滤波融合的实现，首先需要根据观测值选择支持度函数，建立数据间支持度矩阵，并利用数据间支持度矩阵给每个传感器制定相应的加权因子w_i，(i＝1，2，L，N)，然后将各个传感器的权因子乘以其相应的观测值，分别滤波获取对应的滤波值，最后将它们全部相加，从而获得观测系数化后的滤波融合值Z_fil(k|k)。这里的滤波采用的是交互式多模型滤波算法。

图2为观测系数化的滤波融合模型，获得观测系数化的滤波融合值的具体过程如下：

S21.计算相对距离矩阵D(k)中元素的最大值D_Max(k)；

S22.计算支持度函数r_ij(k)，i，j＝1，2，L，N；

可以看出，数据间的相对距离越大，支持度值就越小，即当数据间的相对距离为最大值时，支持度函数值为零；当数据对自身的相对距离为零时，数据对自身的支持度函数值为1。

S23.建立支持度矩阵，N个传感器数据之间的支持度矩阵Sup_Matrix(k)为：

$\mathrm{Sup}\_\mathrm{Matrix}\left(k\right)={\left[\begin{array}{cccc}{r}_{11}\left(k\right)& r_{12}\left(k\right)& L& r_{1N}\left(k\right)\\ r_{21}\left(k\right)& r_{22}\left(k\right)& L& r_{2N}\left(k\right)\\ L& L& L& L\\ r_{N1}\left(k\right)& r_{N2}\left(k\right)& L& r_{\mathrm{NN}}\left(k\right)\end{array}\right]}_{N\×N};$

可以看出，支持度矩阵Sup_Matrix(k)仅表示两数据间的相互支持度，并不能反映一个观测值对所有数据的总体支持度。

S24.计算权因子，具体过程如下：存在一组非负数：l₁，l₂，…，l_N，使得权因子为：w_i＝l₁r_i1+l₂r_i2+L+l_Nr_iN i＝1，2，L，N，满足

其矩阵形式为：W＝Sup_Matrix(k)*L，其中，W＝[w₁ w₂ L w_N]^T，L＝[l₁ l₂ L l_N]^T，建立矩阵方程λL＝Sup_Matrix(k)*L，计算特征值λ及相应的特征向量，取矩阵最大模特征值max(|λ|)相对应的特征向量，得到一组非负系数l₁，l₂，…，l_N，进而求得一个观测值对所有数据的总体支持度，即获得第i个传感器观测值的权因子可为：

i＝1，2，L，N；

至此将步骤S24得到的各个传感器观测值的权因子w_i(i＝1，2，L，N)乘以其相应的观测值，分别利用交互式多模型算法获取对应的滤波值，分别记为：wz₁_fil(k|k)，wz₂_fil(k|k)，L，wz_N_fil(k|k)，将它们全部相加获得观测系数化后的滤波融合值Z_fil(k|k)，记为：Z_fil(k|k)＝wz₁(k|k)+wz₂(k|k)+L+wz_N(k|k)。交互式多模型算法为本领域的公知技术常识，不再进行详细描述。

步骤S4分步式滤波融合的实现。分步式滤波融合作为整个融合方法的第二步，是提高航迹融合精度的关键，主要过程是在已获得k-1时刻基于全局系统的目标状态估计值和估计误差协方差P(k-1|k-1)下，求得相应的预测值，然后利用Kalman滤波器和观测系数化后的滤波融合值Z_fil(k|k)，依次更新k时刻的目标状态估计值，最后得到基于全局的融合估计值

和相应的估计误差协方差矩阵P(k|k)。图4为分步式滤波融合的基本实现流程，具体过程如下：

S41.状态初始化，当k＝1时，利用等权因子加权思想求得全局状态估计值

及估计误差协方差P(1|1)分别为：

$\hat{X}\left(1\right|1)=\frac{{\hat{X}}_1\left(1\right|1)+{\hat{X}}_2\left(1\right|1)+L+{\hat{X}}_N\left(1\right|1)}{N}$

$P\left(1\right|1)=\frac{P_1\left(1\right|1)+P_2\left(1\right|1)+L+P_N\left(1\right|1)}{N^2}$

其中，

及P_i(1|1)分别为第i个传感器的局部估计和相应的估计误差协方差初始值；

S42.预测值求解，当k＝2，3，L，t/T，其中，t为观测时间，T为融合中心采样周期，基于全局系统的状态估计值

和估计误差协方差P(k-1|k-1)的一步预测估计值及相应的预测误差协方差P(k|k-1)分别为：

P(k|k-1)＝FP(k-1|k-1)F^T+GQG^T，其中，F为状态转移矩阵、G为过程噪声输入矩阵、Q为系统噪声协方差矩阵，*^T表示矩阵*的转置运算。这里F、G、Q根据系统状态方程和目标观测方程确定，为本领域的公知技术常识。

S43.分步式滤波更新，利用步骤S3得到的观测系数化后的滤波融合值，记为Z_fil(k|k)，对一步预测估计值及相应预测误差协方差P(k|k-1)依次更新，该步骤基本流程如图5所示，具体实现过程如下：

S431.令

P₁(k|k-1)＝P(k|k-1)，根据Kalman滤波公式有：

${\hat{X}}_1\left(k\right|k)=F\hat{X}(k-1|k-1)+K_1\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_1\left(k\right){\hat{X}}_1\left(k\right|k-1)]$

P₁(k|k)=[I-K₁(k)H₁(k)]P₁(k|k-1)

其中，H₁(k)、R₁(k)分别为第一个传感器的量测矩阵及观测噪声协方差，I为单位矩阵，K₁(k)为第一个传感器的滤波增益矩阵，即有：

$K_1\left(k\right)=P_1\left(k\right|k-1)H_1^T\left(k\right){[H_1\left(k\right)P_1\left(k\right|k-1)H_1^T\left(k\right)+R_1\left(k\right)]}^{-1};$

S432.令

P₂(k|k-1)＝P₁(k|k)，则根据Kalman滤波公式有：

${\hat{X}}_2\left(k\right|k)=\widehat{X_1}\left(k\right|k)+K_2\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_2\left(k\right){\hat{X}}_2\left(k\right|k-1)]$

$=F\hat{X}(k-1|k-1)+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{K}_i\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_i\left(k\right){\hat{X}}_i\left(k\right|k-1)]$

$P_2\left(k\right|k)=[I-K_2\left(k\right)H_2\left(k\right)]P_2\left(k\right|k-1)=\underset{i=1}{\overset{2}{\mathrm{\Π}}}[I-K_i\left(k\right)H_i\left(k\right)]P_1\left(k\right|k-1)$

其中，H₂(k)、R₂(k)为第二个传感器的量测矩阵及观测噪声协方差，K₂(k)为第二个传感器的滤波增益矩阵，即有：

S433.依次可以得到第i个传感器的状态更新估计值

及相应的误差协方差P_i(k|k)，i＝3，……，N，第N个传感器的状态更新估计值

及相应的误差协方差P_N(k|k)为：

${\hat{X}}_N\left(k\right|k)=\widehat{X_{N-1}}\left(k\right|k)+K_N\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_N\left(k\right){\hat{X}}_N\left(k\right|k-1)]$

$=F\hat{X}(k-1|k-1)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{K}_i\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_i\left(k\right){\hat{X}}_i\left(k\right|k-1)]$

$P_N\left(k\right|k)=[I-K_N\left(k\right)H_N\left(k\right)]P_N\left(k\right|k-1)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}[I-K_i\left(k\right)H_i\left(k\right)]P_N\left(k\right|k-1)$

其中，H_N(k)、R_N(k)为第N个传感器的量测矩阵及观测噪声协方差，K_N(k)为第N个传感器的滤波增益矩阵，即有：

S44.计算全局融合估计值，可求得k时刻基于全局系统的融合估计值及相应的估计误差协方差P(k|k)分别为：

P(k|k)＝P_N(k|k)。

在典型航迹蛇形航路条件下，以均方根误差作为目标跟踪性能指标对本发明的方法进行仿真验证，仿真结果如图6、7、8所示。其中图6表示多个传感器的融合跟踪与单个传感器下的局部跟踪的误差曲线图，图7表示本发明的方法与现有的基于分步式滤波的航迹融合方法的误差曲线图，图8表示观测系数化的滤波融合与单个传感器下观测值的误差曲线图。由图6、7、8中可以看出，在目标机动情况下，随着运动速度的增加，每当目标发生机动时，相较于任何单个传感器的局部跟踪，多传感器对目标的融合跟踪精度得到了明显的提高，且由于利用了观测系数化后的滤波融合值作为分步式融合的状态更新输入值，使得该方法在提高观测精度的同时，融合跟踪性能比现有的基于分步式滤波的航迹融合方法得到了明显的改善，从而更好地保证了对目标的可靠跟踪。

可以看出，本发明的方法通过对多传感器的观测信息进行观测系数化的滤波融合，在不增加数据处理复杂度的条件下，降低了观测信息不确定性对航迹融合的影响，同时在观测系数化滤波融合过程，考虑了观测信息的相关性，进而提高了观测精度，获得了对目标的可靠跟踪。本发明方法的主要特点是对多源信息进行多层次融合，并使用观测系数化后的滤波融合值作为状态更新输入值，同时本发明适用于对多种目标运动形式的信息融合，在工程应用上具有可实施性。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (5)

1.一种航迹融合方法，其特征在于，包括如下步骤：

S1.利用预处理后的多传感器的观测信息建立数据间的相对距离矩阵；

S2.依据步骤S1得到的相对距离矩阵，计算支持度函数，获得支持度矩阵，建立方程组，求解加权因子；

S3.将步骤S2得到的各个传感器的权因子与其相应的观测值相乘，分别滤波获取对应的滤波值，将获取的全部滤波值相加得到观测系数化后的滤波融合值；

S4.将步骤S3得到的滤波融合值作为状态更新输入值，利用Kalman滤波，分步对目标状态估计值进行估计更新。

2.根据权利要求1所述的航迹融合方法，其特征在于，所述步骤S1具体过程如下：计算数据间相对距离，建立相对距离矩阵，在k时刻，利用多传感器对同一目标的观测值，计算相对距离d_ij(k)：

d_ij(k)＝|z_i(k)-z_j(k)|i，j＝1，2，L，N，其中，z_i(k)表示预处理后的第i个传感器的观测信息，N是传感器的个数；

在k时刻N个传感器的相对距离矩阵D(k)记为：

$D\left(k\right)={\left[\begin{array}{cccc}{d}_{11}\left(k\right)& d_{12}\left(k\right)& L& d_{1N}\left(k\right)\\ d_{21}\left(k\right)& d_{22}\left(k\right)& L& d_{2N}\left(k\right)\\ L& L& L& L\\ d_{N1}\left(k\right)& d_{N2}\left(k\right)& L& d_{\mathrm{NN}}\left(k\right)\end{array}\right]}_{N\×N}.$

3.根据权利要求2所述的航迹融合方法，其特征在于，所述步骤S2具体包括如下分步骤：

S21.计算相对距离矩阵D(k)中元素的最大值D_Max(k)；

S22.计算支持度函数r_ij(k)，

i，j＝1，2，L，N；

S23.建立支持度矩阵，N个传感器数据之间的支持度矩阵Sup_Matrix(k)为：

$\mathrm{Sup}\_\mathrm{Matrix}\left(k\right)={\left[\begin{array}{cccc}{r}_{11}\left(k\right)& r_{12}\left(k\right)& L& r_{1N}\left(k\right)\\ r_{21}\left(k\right)& r_{22}\left(k\right)& L& r_{2N}\left(k\right)\\ L& L& L& L\\ r_{N1}\left(k\right)& r_{N2}\left(k\right)& L& r_{\mathrm{NN}}\left(k\right)\end{array}\right]}_{N\×N};$

S24.计算权因子，具体过程如下：存在一组非负数：l₁，l₂，…，l_N，使得权因子为：w_i＝l₁r_i1+l₂r_i2+L+l_Nr_iN，i＝1，2，L，N，并且满足

其矩阵形式为：W＝Sup_Matrix(k)*L，J，其中，W＝[w₁ w₂ L w_N]^T，L＝[l₁ l₂ L l_N]^T，建立矩阵方程λL＝Sup_Matrix(k)*L，计算特征值λ及相应的特征向量，取矩阵最大模特征值max(|λ|)相对应的特征向量，得到一组非负系数，l₁，l₂，…，l_N，进而求得一个观测值对所有数据的总体支持度，即获得第i个传感器观测值的权因子可为：

i＝l，2，L，N。

4.根据权利要求2或3所述的航迹融合方法，其特征在于，所述步骤S3中的滤波采用的是交互式多模型滤波算法。

5.根据权利要求4所述的航迹融合方法，其特征在于，所述步骤S4具体包括如下步骤：

S41.状态初始化，当k＝l时，利用等权因子加权思想求得全局状态估计值

及估计误差协方差P(1|1)分别为：

$\hat{X}\left(1\right|1)=\frac{{\hat{X}}_1\left(1\right|1)+{\hat{X}}_2\left(1\right|1)+L+{\hat{X}}_N\left(1\right|1)}{N}$

$P\left(1\right|1)=\frac{P_1\left(1\right|1)+P_2\left(1\right|1)+L+P_N\left(1\right|1)}{N^2}$

其中，

及P_i(1|1)分别为第i个传感器的局部估计和相应的估计误差协方差初始值；

S42.预测值求解，当k＝2，3，L，t/T，其中，t为观测时间，T为融合中心采样周期，基于全局系统的状态估计值

和估计误差协方差P(k-l|k-1)的一步预测估计值

及相应的预测误差协方差P(k|k-1)分别为：

P(k|k-1)＝FP(k-l|k-1)F^T+GQG^T，其中，F为状态转移矩阵、G为过程噪声输入矩阵、Q为系统噪声协方差矩阵，*^T表示矩阵*的转置运算；

S43.分步式滤波更新，利用步骤S3得到的观测系数化后的滤波融合值，记为Z_fil(k|k)，对一步预测估计值

及相应预测误差协方差P(k|k-1)依次更新，具体实现过程如下：

S431.令

P₁(k|k-1)＝P(k|k-1)，根据Kalman滤波公式有：

${\hat{X}}_1\left(k\right|k)=F\hat{X}(k-1|k-1)+K_1\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_1\left(k\right){\hat{X}}_1\left(k\right|k-1)]$

P₁(k|k)＝[I-K₁(k)H₁(k)]P₁(k|k-1)

其中，H₁(k)、R₁(k)分别为第一个传感器的量测矩阵及观测噪声协方差，I为单位矩阵，K₁(k)为第一个传感器的滤波增益矩阵，即有：

$K_1\left(k\right)=P_1\left(k\right|k-1)H_1^T\left(k\right){[H_1\left(k\right)P_1\left(k\right|k-1)H_1^T\left(k\right)+R_1\left(k\right)]}^{-1};$

S432.令

P₂(k|k-1)＝P₁(k|k)，则根据Kalman滤波公式有：

${\hat{X}}_2\left(k\right|k)=\widehat{X_1}\left(k\right|k)+K_2\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_2\left(k\right){\hat{X}}_2\left(k\right|k-1)]$

$=F\hat{X}(k-1|k-1)+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{K}_i\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_i\left(k\right){\hat{X}}_i\left(k\right|k-1)]$

$P_2\left(k\right|k)=[I-K_2\left(k\right)H_2\left(k\right)]P_2\left(k\right|k-1)=\underset{i=1}{\overset{2}{\mathrm{\Π}}}[I-K_i\left(k\right)H_i\left(k\right)]P_1\left(k\right|k-1)$

其中，H₂(k)、R₂(k)为第二个传感器的量测矩阵及观测噪声协方差，K₂(k)为第二个传感器的滤波增益矩阵，即有：

S433.依次可以得到第i个传感器的状态更新估计值

及相应的误差协方差P_i(k|k)，i＝3，……，N，第N个传感器的状态更新估计值

及相应的误差协方差P_N(k|k)为：

${\hat{X}}_N\left(k\right|k)=\widehat{X_{N-1}}\left(k\right|k)+K_N\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_N\left(k\right){\hat{X}}_N\left(k\right|k-1)]$

$=F\hat{X}(k-1|k-1)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{K}_i\left(k\right)[Z\_\mathrm{fil}\left(k\right|k)-H_i\left(k\right){\hat{X}}_i\left(k\right|k-1)]$

$P_N\left(k\right|k)=[I-K_N\left(k\right)H_N\left(k\right)]P_N\left(k\right|k-1)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}[I-K_i\left(k\right)H_i\left(k\right)]P_N\left(k\right|k-1)$

其中，H_N(k)、R_N(k)为第N个传感器的量测矩阵及观测噪声协方差，K_N(k)为第N个传感器的滤波增益矩阵，即有：

S44.计算全局融合估计值，可求得k时刻基于全局系统的融合估计值

及相应的估计误差协方差P(k|k)分别为：

Images