CN105929391A - 基于后验概率的变结构交互多模型雷达机动目标跟踪方法 - Google Patents

Abstract

本发明属于雷达机动目标跟踪领域，公开了一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，包括：初始化雷达观测值，得到初始状态向量和初始协方差矩阵；选择初始角速度，确定运动模型的状态转移矩阵；由运动模型的个数假定初始先验概率和初始模型概率转移矩阵；进行交互多模型滤波得到状态估计向量和估计协方差矩阵，同时也得到模型的后验概率；利用变结构和图论的思想，通过模型的后验概率对角速度即多模型结构进行调整；通过模型的后验概率对模型概率转移矩阵进行调整，直至跟踪完成。本发明跟踪性能好、计算量小，可用于对高速强机动目标的跟踪。

Description

基于后验概率的变结构交互多模型雷达机动目标跟踪方法

技术领域

本发明涉及雷达机动目标跟踪领域，尤其涉及一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，可用于对高速机动目标的跟踪。

背景技术

由于目标跟踪在国防和民用各领域的重要价值，一直以来很多学者致力于研究目标跟踪算法。随着现代航空的迅速发展，各种飞行器的机动性和运动速度越来越高，但是可靠而精确的跟踪目标始终是目标跟踪系统设计的主要目的与难点。目标跟踪实际上就是对目标状态的跟踪滤波问题，而实现目标跟踪首先要使所建立的目标运动模型与实际的目标运动模型匹配，建立目标的运动模型和自适应滤波是机动目标的两个关键部分。

在当前目标机动能力日益增强的情况下，单模型算法很难满足跟踪速度和精度的需求，而多模型算法可以避免采用单模型时由于目标机动而造成模型的不准确，提高机动目标跟踪性能，从而实现对强机动目标的精确跟踪。经过国内外学者几十年来的不断研究，提出了很多目标模型与算法。目前广泛采用的算法是1984-1989年Blom和Bar-Shalom提出的交互式多模型IMM算法，该算法具有Markov转移概率，且算法中有多个模型并行工作，多个滤波器交互作用得到目标状态估计的结果，由于该模型有较好的自适应能力，因而比较理想。但是IMM算法中由于其模型转移概率是人为设置的，因此引入了人为误差，影响了跟踪性能。除此之外，因为IMM算法中每一个模型都固定不变，在所设的模型都不符合实际运动模型的情况下，IMM算法跟踪性能无疑会很差。现在有为此改进的变结构IMM算法，虽然解决了IMM算法中模型固定不变的问题，但需要提前设好尽可能多的模型集，而模型集增删规则也十分复杂，运算量十分庞大。

发明内容

针对上述现有技术的不足，本发明的目的在于提供一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，能够增强滤波算法的自适应能力，提高雷达在目标强机动下的跟踪精度。

实现本发明的技术思路是：利用变结构和图论的思想，通过交互多模型的后验概率调整滤波器的结构以及转移概率，使多模型采用的机动模型集向着目标真实运动模型逼近，进而使估计信号逼近期望信号，以提高跟踪性能。

为达到上述目的，本发明的实施例采用如下技术方案：

一种基于后验概率反馈的变结构交互多模型机动目标跟踪方法，所述方法包括如下步骤：

步骤1，获取雷达机动目标的位置检测值，对所述位置检测值进行N次采样，得到观测值序列，根据所述观测值序列确定雷达机动目标的第0时刻状态向量和第0时刻协方差矩阵，N为大于1的自然数；

步骤2，雷达机动目标运动模型采用联动式转弯运动模型，对所述联动式转弯运动模型设定r个不同的初始角速度，得到r个对应的目标运动子模型和r个对应的系统状态转移矩阵；分别确定r个目标运动子模型第0时刻先验概率矩阵、第0时刻模型状态转移矩阵；

步骤3，对于第j个目标运动子模型，根据第k-1时刻状态向量和第k-1时刻协方差矩阵，以及其对应的第k-1时刻先验概率矩阵、第k-1时刻模型状态转移矩阵，得到第j个目标运动子模型最新的第k-1时刻状态向量和第k-1时刻协方差矩阵，j∈(1,2,…r)，k≥1，k的初值为1；

步骤4，根据第j个目标运动子模型最新的第k-1时刻状态向量、第k-1时刻协方差矩阵、所述观测值序列和其对应的系统状态转移矩阵，计算第j个目标运动子模型的第k时刻新息值、第k时刻新息协方差矩阵以及第k时刻增益矩阵，从而得到第j个目标运动子模型的第k时刻状态估计向量和第k时刻状态估计协方差矩阵；

步骤5，根据所述第j个目标运动子模型的第k时刻新息值和第k时刻新息协方差矩阵计算第k时刻新息的似然函数；根据所述第k时刻新息的似然函数、所述第j个目标运动子模型的第k-1时刻先验概率矩阵和第k-1时刻模型状态转移矩阵计算所述第j个目标运动子模型的第k时刻后验概率；

步骤6，根据步骤3至步骤5，依次得到r个目标运动子模型的第k时刻状态估计向量、第k时刻状态估计协方差矩阵以及第k时刻后验概率，根据所述r个目标运动子模型的第k时刻状态估计向量、第k时刻状态估计协方差矩阵以及第k时刻后验概率得到r个目标运动子模型第k时刻总体状态估计向量和第k时刻总体状态估计协方差矩阵；

步骤7，根据r个目标运动子模型第k时刻后验概率对r个目标运动子模型的角速度进行修正，并将修正后的角速度作为r个目标运动子模型的最新角速度；

步骤8，根据r个目标运动子模型第k时刻后验概率和第k-1时刻后验概率对r个目标运动子模型第k-1时刻模型状态转移矩阵进行修正，得到r个目标运动子模型第k时刻模型状态转移矩阵，并将其作为r个目标运动子模型最新的模型状态转移矩阵；

步骤9，令k的值加1，并依次重复执行步骤3至步骤8，直到k＞N-1，从而得到r个目标运动子模型对雷达机动目标跟踪的最终模型状态转移矩阵，根据r个目标运动子模型的最终模型状态转移矩阵对雷达机动目标进行跟踪。

本发明与现有技术相比的优点在于：1.跟踪性能高：本发明中引入了变结构和图论的思想，与传统IMM算法中模型集和转移概率固定不同的是，本发明中滤波器可以利用多模型的后验概率自适应调节模型集和转移概率，从而提高了跟踪精度，减小了误差，同时提高了对机动目标的跟踪性能；2.计算简单：现有变结构交互多模型算法都经过一系列复杂算法来达到增删模型集的目的，虽然可能对更复杂的机动情况达到全面的自适应，更有可能使结果背道而驰。本发明仅使用多模型的后验概率对模型的某一参数进行修正就达到修改模型集的目的，大幅度降低了计算复杂度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法的流程示意图；

图2为本发明对雷达机动目标跟踪的效果示意图；

图3为本发明与传统IMM方法对雷达机动目标跟踪三维位置误差比较示意图；

图4为本发明与传统IMM方法对雷达机动目标跟踪三维速度误差比较示意图。

具体实施方式

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

一种基于后验概率反馈的变结构交互多模型机动目标跟踪方法，参照图1，所述方法包括如下步骤：

对机动目标运动状态进行初始化：

步骤1，步骤1，获取雷达机动目标的位置检测值，对所述位置检测值进行N次采样，得到观测值序列，根据所述观测值序列确定雷达机动目标的第0时刻状态向量和第0时刻协方差矩阵，N为大于1的自然数。

步骤1具体包括：

(1a)对雷达机动目标的位置进行检测，得到雷达机动目标的位置检测值，并对所述位置检测值进行N次采样，得到观测值序列Z(k)，k＝1,2,…,N；

从观测值序列Z(k)中取出前三个采样点的观测值Z(1)，Z(2)，Z(3)，其中：

Z(1)＝[xx(1) yy(1) zz(1)]^T

Z(2)＝[xx(2) yy(2) zz(2)]^T

Z(3)＝[xx(3) yy(3) zz(3)]^T

式中xx(·)，yy(·)，zz(·)分别表示雷达机动目标在x轴、y轴、z轴的位置，[·]^T表示矩阵的转置；

(1b)根据所述观测值序列Z(k)确定雷达机动目标的第0时刻状态向量根据雷达机动目标的第0时刻状态向量确定雷达机动目标的第0时刻协方差矩阵P₀：

根据观测值Z(1)，Z(2)，Z(3)，获得雷达机动目标的第0时刻状态向量如下：

${\hat{X}}_0=\left[\begin{array}{c}xx\left(3\right)\\ \left(xx\right(3)-xx(2\left)\right)/T\\ \left(\right(xx\left(3\right)-xx\left(2\right))/T-(xx\left(2\right)-xx\left(1\right))/T)/T\\ yy\left(3\right)\\ \left(yy\right(3)-yy(2\left)\right)/T\\ \left(\right(yy\left(3\right)-yy\left(2\right))/T-(yy\left(2\right)-yy\left(1\right))/T)/T\\ zz\left(3\right)\\ \left(zz\right(3)-zz(2\left)\right)/T\\ \left(\right(zz\left(3\right)-zz\left(2\right)/T)-(zz\left(2\right)-zz\left(1\right)/T\left)\right)/T\end{array}\right],$

其中，(xx(3)-xx(2))/T、(yy(3)-yy(2))/T、(zz(3)-zz(2))/T分别表示雷达机动目标在第3采样点处的x轴方向速度、y轴方向速度、z轴方向速度；((xx(3)-xx(2))/T-(xx(2)-xx(1))/T)/T、((yy(3)-yy(2))/T-(yy(2)-yy(1))/T)/T、((zz(3)-zz(2)/T)-(zz(2)-zz(1)/T))/T分别表示雷达机动目标在第3采样点的x轴方向加速度、y轴方向加速度、z轴方向加速度，T为采样间隔；

(1c)根据第0时刻状态向量获得第0时刻协方差矩阵P₀：其中，E[·]表示数学期望。

对机动目标跟踪模型参数进行初始化：

步骤2，雷达机动目标运动模型采用联动式转弯运动模型，对所述联动式转弯运动模型设定r个不同的初始角速度，得到r个对应的目标运动子模型和r个对应的系统状态转移矩阵；分别确定r个目标运动子模型第0时刻先验概率矩阵、第0时刻模型状态转移矩阵。

根据机动目标的运动特性，可知协CT(联动式转弯运动)模型因角速度不同几乎可以包含所有运动可能状态。

步骤2具体包括：

(2a)对雷达机动目标运动模型设定r个不同的初始角速度ω_i，从而得到r个对应的系统状态转移矩阵Φ_i：

${\mathrm{\Φ}}_{i=1,2,...,r}=\left[\begin{array}{ccc}1& \frac{\sin \left({\mathrm{\ω}}_{i}T\right)}{{\mathrm{\ω}}_i}& \frac{1-\cos \left({\mathrm{\ω}}_{i}T\right)}{{{\mathrm{\ω}}_i}^2}\\ 0& \cos \left({\mathrm{\ω}}_{i}T\right)& \frac{\sin \left({\mathrm{\ω}}_{i}T\right)}{{\mathrm{\ω}}_i}\\ 0& -{\mathrm{\ω}}_i\sin \left({\mathrm{\ω}}_{i}T\right)& \cos \left({\mathrm{\ω}}_{i}T\right)\end{array}\right]$

根据目标运动子模型的个数确定r个目标运动子模型的第0时刻先验概率μ矩阵和第0时刻模型状态转移矩阵p：

μ＝[1/r 1/r … 1/r]_1×r

$p={\left[\begin{array}{cccc}0.8& 0.2/r-1& ...& 0.2/r-1\\ 0.2/r-1& 0.8& ...& 0.2/r-1\\ ...& ...& M& ...\\ 0.2/r-1& 0.2/r-1& ...& 0.8\end{array}\right]}_{r\×r}$

其中，i＝1,2,…r，[·]_1×r表示该矩阵为1行r列的矩阵，[·]_r×r表示该矩阵为r行r列的矩阵。

输入交互，即以模型j为例j∈(1,2,…r)进行如下步骤：

步骤3，对于第j个目标运动子模型，根据第k-1时刻状态向量和第k-1时刻协方差矩阵，以及其对应的第k-1时刻先验概率矩阵、第k-1时刻模型状态转移矩阵，得到第j个目标运动子模型最新的第k-1时刻状态向量和第k-1时刻协方差矩阵，j∈(1,2,…r)，k≥1，k的初值为1。

步骤3具体包括：

对于第j个目标运动子模型，j∈(1,2,…r)，根据第k-1时刻先验概率矩阵、第k-1时刻模型状态转移矩阵、第k-1时刻状态向量和第k-1时刻协方差矩阵，得到第k-1时刻第j个目标运动子模型的最新的第k-1时刻状态向量和第k-1时刻协方差矩阵

${\mathrm{\μ}}_{i|j}(k-1|k-1)=\frac{p_{ij}(k-1){\mathrm{\μ}}_i(k-1)}{\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{p}_{ij}(k-1){\mathrm{\μ}}_i(k-1)}$

${\hat{X}}_j^0(k-1|k-1)=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\hat{X}}_i(k-1|k-1){\mathrm{\μ}}_{i|j}(k-1|k-1)$

$\begin{array}{c}{P}_j^0(k-1|k-1)=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i|j}(k-1|k-1)\\ \·\{P_i(k-1|k-1)+\[{\hat{X}}_i(k-1|k-1)-\\ {\hat{X}}_j^0(k-1|k-1)\]\·{\[{\hat{X}}_i(k-1|k-1)-{\hat{X}}_j^0(k-1|k-1)\]}^T\}\end{array}$

其中，μ_i(k-1)为第k-1时刻目标运动子模型i发生的概率，p_ij(k-1)为第k-1时刻模型i转换到模型j的概率，r代表目标运动子模型的总数，μ_i|j(k-1|k-1)为第k-1时刻在转换结果为目标运动子模型j的条件下目标运动子模型i发生的概率。

对于模型j，以步骤3求得的及观测值Z(k)作为输入按照如下步骤进行滤波：

步骤4，根据第j个目标运动子模型最新的第k-1时刻状态向量、第k-1时刻协方差矩阵、所述观测值序列和其对应的系统状态转移矩阵，计算第j个目标运动子模型的第k时刻新息值、第k时刻新息协方差矩阵以及第k时刻增益矩阵，从而得到第j个目标运动子模型的第k时刻状态估计向量和第k时刻状态估计协方差矩阵。

步骤4具体包括：

(4a)根据第k时刻目标运动子模型j的状态转移矩阵Φ_j(k)对第k-1时刻目标运动子模型j的状态向量进行一步预测，得到对应目标运动子模型的一步状态预测值

(4b)根据状态转移矩阵Φ_j(k)对第k-1时刻目标运动子模型j的状态协方差矩阵进行一步预测，得到对应目标运动子模型的一步协方差预测值P_j(k|k-1)：其中，Q_j表示模型j的过程噪声协方差；

(4c)获取第k时刻目标运动子模型j的量测矩阵H_j(k)，利用第k时刻目标运动子模型j的量测矩阵H_j(k)对目标运动子模型j的一步状态预测值进行量测转换，得到对应目标运动子模型的量测预测值其中，H_j(k)表示雷达机动目标第k时刻目标运动子模型j的量测矩阵；

(4d)利用当前时刻观测值Z(k)和目标运动子模型j的量测预测值相减得到对应目标运动子模型的新息值V_j(k)：

(4e)由目标运动子模型j的一步协方差预测值P_j(k|k-1)和量测矩阵H_j(k)确定对应目标运动子模型的新息协方差矩阵S_j(k)：其中，R(k)是观测噪声在第k时刻的协方差矩阵；

(4f)根据目标运动子模型j的一步协方差预测值P_j(k|k-1)、量测矩阵H_j(k)以及新息协方差矩阵S_j(k)计算增益矩阵K_j(k)：(·)^H表示求矩阵的共轭转置，(·)^-1表示求矩阵的逆；

(4g)利用增益矩阵K_j(k)对目标运动子模型j的一步状态预测值和一步协方差预测值P_j(k|k-1)进行滤波，得到目标运动子模型j的第k时刻状态估计向量和第k时刻状态估计协方差矩阵P_j(k|k)：

${\hat{X}}_j\left(k\right|k)={\hat{X}}_j\left(k\right|k-1)+K_j\left(k\right)V_j\left(k\right)$

$P_j\left(k\right|k)=P_j\left(k\right|k-1)-K_j\left(k\right)S_j\left(k\right){K_j}^T\left(k\right).$

对机动目标跟踪模型进行更新：

步骤5，根据所述第j个目标运动子模型的第k时刻新息值和第k时刻新息协方差矩阵计算第k时刻新息的似然函数；根据所述第k时刻新息的似然函数、所述第j个目标运动子模型的第k-1时刻先验概率矩阵和第k-1时刻模型状态转移矩阵计算所述第j个目标运动子模型的第k时刻后验概率。

步骤5具体包括：

(5a)由目标运动子模型j的新息V_j(k)及其协方差矩阵S_j(k)计算新息的似然函数Λ_j(k)：

(5b)利用模型状态转移矩阵和新息的似然函数Λ_j(k)对目标运动子模型j的先验概率进行更新，得到目标运动子模型j的后验概率μ_j(k)：

${\mathrm{\μ}}_j\left(k\right)=\frac{{\mathrm{\Λ}}_j\left(k\right)\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{p}_{ij}(k-1){\mathrm{\μ}}_i(k-1)}{\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\Λ}}_j\left(k\right)\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{p}_{ij}(k-1){\mathrm{\μ}}_i(k-1)}.$

输出交互：

步骤6，根据步骤3至步骤5，依次得到r个目标运动子模型的第k时刻状态估计向量、第k时刻状态估计协方差矩阵以及第k时刻后验概率，根据所述r个目标运动子模型的第k时刻状态估计向量、第k时刻状态估计协方差矩阵以及第k时刻后验概率得到r个目标运动子模型第k时刻总体状态估计向量和第k时刻总体状态估计协方差矩阵。

步骤6具体包括：

根据所述r个目标运动子模型第k时刻状态估计向量、第k时刻状态估计协方差矩阵以及第k时刻后验概率得到r个目标运动子模型交互后第k时刻状态估计向量和第k时刻状态估计协方差矩阵P(k|k)：

$\hat{X}\left(k\right|k)=\underset{j=1}{\overset{r}{\mathrm{\Σ}}}{\hat{X}}_j\left(k\right|k){\mathrm{\μ}}_j\left(k\right)$

$P\left(k\right|k)=\underset{j=1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\μ}}_j\left(k\right)\·\{P_j\left(k\right|k)+\[{\hat{X}}_j\left(k\right|k)-\hat{X}\left(k\right|k)\]\·{\[{\hat{X}}_j\left(k\right|k)-\hat{X}\left(k\right|k)\]}^T\}.$

按照如下步骤根据后验概率建立变结构模型：

步骤7，根据r个目标运动子模型第k时刻后验概率对r个目标运动子模型的角速度进行修正，并将修正后的角速度作为r个目标运动子模型的最新角速度。

步骤7具体包括：

若目标运动子模型数r为奇数即r＝2n+1，则在第k时刻的模型集为其特征参数分别对应r个目标运动子模型的角速度且在模型集中从小到大排列，称为中心模型；对应模型集的后验概率为同理若目标运动子模型数r为偶数即r＝2n，则在第k时刻的模型集为为中心模型；对应模型集的后验概率为

按照下述方法调整模型间距：

第k时刻模型间距调整参数为 其中δ⁽ⁱ⁾表示模型区分距离，一般取

(7a)模型概率重心无跳变：此时则可根据其他模型概率判定模型集不变或向中心模型靠拢；

其中，t₁表示不可能模型的门限值，一般取t₁<0.1；

(7b)模型概率重心向左跳变：此时则可根据其他模型概率判定模型集不变或向左侧模型靠拢；

${\mathrm{\&omega;}}_{k+1}^{\max \left(i\right)}={\mathrm{\&omega;}}_{k+1}^{mid}+{\mathrm{\&lambda;}}_k^{r\left(i\right)}$

其中，t₂表示有效模型的门限值，一般取t₂＞0.9；

(7c)模型概率重心向右跳变：此时则可根据其他模型概率判定模型集不变或向右侧模型靠拢；

${\mathrm{\&omega;}}_{k+1}^{\min \left(i\right)}={\mathrm{\&omega;}}_{k+1}^{mid}-{\mathrm{\&lambda;}}_k^{l\left(i\right)}.$

按照如下步骤基于后验概率改变转移概率矩阵：

步骤8，根据r个目标运动子模型第k时刻后验概率和第k-1时刻后验概率对r个目标运动子模型第k-1时刻模型状态转移矩阵进行修正，得到r个目标运动子模型第k时刻模型状态转移矩阵，并将其作为r个目标运动子模型最新的模型状态转移矩阵。

步骤8具体包括：

(8a)对于目标运动子模型j，其后验概率在第k时刻与第k-1时刻的比值λ_j(k)为：

${\mathrm{\&lambda;}}_j\left(k\right)=\frac{{\mathrm{\&mu;}}_j\left(k\right)}{{\mathrm{\&mu;}}_j(k-1)}$

(8b)利用前后时刻后验概率比值λ_j(k)对第k-1时刻模型状态转移矩阵进行修正，第k时刻模型状态转移矩阵：

$p_{ij}\left(k\right)=\frac{{\mathrm{\&lambda;}}_j\left(k\right)p_{ij}(k-1)}{\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j\left(k\right)p_{ij}(k-1)}.$

判断跟踪是否完成：

步骤9，令k的值加1，并依次重复执行步骤3至步骤8，直到k＞N-1，从而得到r个目标运动子模型对雷达机动目标跟踪的最终模型状态转移矩阵，根据r个目标运动子模型的最终模型状态转移矩阵对雷达机动目标进行跟踪。

本发明的效果可以通过以下仿真进一步说明：

1.仿真条件：

假设目标的起始位置为[60 40 20]km,起始速度为[-184 184 150]m/s，转弯速度300m/s，在1-60s，156-180s，246-300s目标在x-z平面做左转45°圆周运动；在64-155s目标在三维斜面做左转1.87°圆周运动；在181-245s目标在三维斜面做右转2.8°圆周运动。距离误差的标准差为Δr＝127m，方位角误差的标准差为Δθ＝0.1°,仰角测量误差的标准差为对目标进行300s的观测，采样间隔T为1s，Monte Carlo仿真100次。

选择的运动模型是CT模型。初始模型集选取M₀＝{-π/3,-π/18,-π/3}，模型区分距离δ＝π/18。

2.仿真内容

仿真1，采用本发明方法对仿真条件中所述的目标进行跟踪仿真，结果如图2。从图2可以看出，本发明方法的跟踪轨迹很接近机动目标运动的真实轨迹。

仿真2，对用本发明方法与传统IMM方法的跟踪误差进行比较，图3表示本发明方法和传统IMM算法在距离上的跟踪误差比较；图4表示本发明方法和传统IMM算法在速度上的跟踪误差比较。

图3和图4表明本发明方法的跟踪误差在整体趋势上小于IMM算法的跟踪误差，且跟踪稳定，更适合对高速强机动目标的跟踪。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (9)

1.一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，所述方法包括如下步骤：

步骤1，获取雷达机动目标的位置检测值，对所述位置检测值进行N次采样，得到观测值序列，根据所述观测值序列确定雷达机动目标的第0时刻状态向量和第0时刻协方差矩阵，N为大于1的自然数；

步骤2，雷达机动目标运动模型采用联动式转弯运动模型，对所述联动式转弯运动模型设定r个不同的初始角速度，得到r个对应的目标运动子模型和r个对应的系统状态转移矩阵；分别确定r个目标运动子模型第0时刻先验概率矩阵、第0时刻模型状态转移矩阵；

步骤3，对于第j个目标运动子模型，根据第k-1时刻状态向量和第k-1时刻协方差矩阵，以及其对应的第k-1时刻先验概率矩阵、第k-1时刻模型状态转移矩阵，得到第j个目标运动子模型最新的第k-1时刻状态向量和第k-1时刻协方差矩阵，j∈(1,2,…r)，k≥1，k的初值为1；

步骤4，根据第j个目标运动子模型最新的第k-1时刻状态向量、第k-1时刻协方差矩阵、所述观测值序列和其对应的系统状态转移矩阵，计算第j个目标运动子模型的第k时刻新息值、第k时刻新息协方差矩阵以及第k时刻增益矩阵，从而得到第j个目标运动子模型的第k时刻状态估计向量和第k时刻状态估计协方差矩阵；

步骤5，根据所述第j个目标运动子模型的第k时刻新息值和第k时刻新息协方差矩阵计算第k时刻新息的似然函数；根据所述第k时刻新息的似然函数、所述第j个目标运动子模型的第k-1时刻先验概率矩阵和第k-1时刻模型状态转移矩阵计算所述第j个目标运动子模型的第k时刻后验概率；

步骤6，根据步骤3至步骤5，依次得到r个目标运动子模型的第k时刻状态估计向量、第k时刻状态估计协方差矩阵以及第k时刻后验概率，根据所述r个目标运动子模型的第k时刻状态估计向量、第k时刻状态估计协方差矩阵以及第k时刻后验概率得到r个目标运动子模型第k时刻总体状态估计向量和第k时刻总体状态估计协方差矩阵；

步骤7，根据r个目标运动子模型第k时刻后验概率对r个目标运动子模型的角速度进行修正，并将修正后的角速度作为r个目标运动子模型的最新角速度；

步骤8，根据r个目标运动子模型第k时刻后验概率和第k-1时刻后验概率对r个目标运动子模型第k-1时刻模型状态转移矩阵进行修正，得到r个目标运动子模型第k时刻模型状态转移矩阵，并将其作为r个目标运动子模型最新的模型状态转移矩阵；

步骤9，令k的值加1，并依次重复执行步骤3至步骤8，直到k>N-1，从而得到r个目标运动子模型对雷达机动目标跟踪的最终模型状态转移矩阵，根据r个目标运动子模型的最终模型状态转移矩阵对雷达机动目标进行跟踪。

2.根据权利要求1所述的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，步骤1具体包括：

(1a)对雷达机动目标的位置进行检测，得到雷达机动目标的位置检测值，并对所述位置检测值进行N次采样，得到观测值序列Z(k)，k＝1,2,…,N；

从观测值序列Z(k)中取出前三个采样点的观测值Z(1)，Z(2)，Z(3)，其中：

Z(1)＝[xx(1) yy(1) zz(1)]^T

Z(2)＝[xx(2) yy(2) zz(2)]^T

Z(3)＝[xx(3) yy(3) zz(3)]^T

式中xx(·)，yy(·)，zz(·)分别表示雷达机动目标在x轴、y轴、z轴的位置，[·]^T表示矩阵的转置；

(1b)根据所述观测值序列Z(k)确定雷达机动目标的第0时刻状态向量根据雷达机动目标的第0时刻状态向量确定雷达机动目标的第0时刻协方差矩阵P₀：

根据观测值Z(1)，Z(2)，Z(3)，获得雷达机动目标的第0时刻状态向量如下：

${\hat{X}}_0=\left[\begin{array}{c}xx\left(3\right)\\ \left(xx\right(3)-xx(2\left)\right)/T\\ \left(\right(xx\left(3\right)-xx\left(2\right))/T-(xx\left(2\right)-xx\left(1\right))/T)/T\\ yy\left(3\right)\\ \left(yy\right(3)-yy(2\left)\right)/T\\ \left(\right(yy\left(3\right)-yy\left(2\right))/T-(yy\left(2\right)-yy\left(1\right))/T)/T\\ zz\left(3\right)\\ \left(zz\right(3)-zz(2\left)\right)/T\\ \left(\right(zz\left(3\right)-zz\left(2\right)/T)-(zz\left(2\right)-zz\left(1\right)/T\left)\right)/T\end{array}\right],$

其中，(xx(3)-xx(2))/T、(yy(3)-yy(2))/T、(zz(3)-zz(2))/T分别表示雷达机动目标在第3采样点处的x轴方向速度、y轴方向速度、z轴方向速度；((xx(3)-xx(2))/T-(xx(2)-xx(1))/T)/T、((yy(3)-yy(2))/T-(yy(2)-yy(1))/T)/T、((zz(3)-zz(2)/T)-(zz(2)-zz(1)/T))/T分别表示雷达机动目标在第3采样点的x轴方向加速度、y轴方向加速度、z轴方向加速度，T为采样间隔；

(1c)根据第0时刻状态向量获得第0时刻协方差矩阵P0：其中，E[·]表示数学期望。

3.根据权利要求1所述的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，步骤2具体包括：

(2a)对雷达机动目标运动模型设定r个不同的初始角速度ω_i，从而得到r个对应的系统状态转移矩阵Φ_i：

${\mathrm{\&Phi;}}_{i=1,2,\mathrm{...},r}=\left[\begin{array}{ccc}1& \frac{\sin \left({\mathrm{\&omega;}}_{i}T\right)}{{\mathrm{\&omega;}}_i}& \frac{1-\cos \left({\mathrm{\&omega;}}_{i}T\right)}{{{\mathrm{\&omega;}}_i}^2}\\ 0& \cos \left({\mathrm{\&omega;}}_{i}T\right)& \frac{\sin \left({\mathrm{\&omega;}}_{i}T\right)}{{\mathrm{\&omega;}}_i}\\ 0& -{\mathrm{\&omega;}}_i\sin \left({\mathrm{\&omega;}}_{i}T\right)& \cos \left({\mathrm{\&omega;}}_{i}T\right)\end{array}\right]$

根据目标运动子模型的个数确定r个目标运动子模型的第0时刻先验概率μ矩阵和第0时刻模型状态转移矩阵p：

μ＝[1/r 1/r …1/r]_1×r

$p={\left[\begin{array}{cccc}0.8& 0.2/r-1& ...& 0.2/r-1\\ 0.2/r-1& 0.8& ...& 0.2/r-1\\ ...& ...& M& ...\\ 0.2/r-1& 0.2/r-1& ...& 0.8\end{array}\right]}_{r\&times;r}$

其中，i＝1,2,…r，[·]_1×r表示该矩阵为1行r列的矩阵，[·]_r×r表示该矩阵为r行r列的矩阵。

4.根据权利要求1所述的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，步骤3具体包括：

对于第j个目标运动子模型，j∈(1,2,…r)，根据第k-1时刻先验概率矩阵、第k-1时刻模型状态转移矩阵、第k-1时刻状态向量和第k-1时刻协方差矩阵，得到第k-1时刻第j个目标运动子模型的最新的第k-1时刻状态向量和第k-1时刻协方差矩阵

${\mathrm{\&mu;}}_{i|j}(k-1|k-1)=\frac{p_{ij}(k-1){\mathrm{\&mu;}}_i(k-1)}{\underset{i=1}{\overset{r}{\&Sigma;}}{p}_{ij}(k-1){\mathrm{\&mu;}}_i(k-1)}$

${\hat{X}}_j^0(k-1|k-1)=\underset{i=1}{\overset{r}{\&Sigma;}}{\hat{X}}_i(k-1|k-1){\mathrm{\&mu;}}_{i|j}(k-1|k-1)$

$\begin{array}{c}{P}_j^0(k-1|k-1)=\underset{i=1}{\overset{r}{\&Sigma;}}{\mathrm{\&mu;}}_{i|j}(k-1|k-1)\\ \&CenterDot;\{P_i(k-1|k-1)+\&lsqb;{\hat{X}}_i(k-1|k-1)-\\ {\hat{X}}_j^0(k-1|k-1)\&rsqb;\&CenterDot;{\&lsqb;{\hat{X}}_i(k-1|k-1)-{\hat{X}}_j^0(k-1|k-1)\&rsqb;}^T\}\end{array}$

其中，μ_i(k-1)为第k-1时刻目标运动子模型i发生的概率，p_ij(k-1)为第k-1时刻模型i转换到模型j的概率，r代表目标运动子模型的总数，μ_i|j(k-1|k-1)为第k-1时刻在转换结果为目标运动子模型j的条件下目标运动子模型i发生的概率。

5.根据权利要求1所述的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，步骤4具体包括：

(4a)根据第k时刻目标运动子模型j的状态转移矩阵Φ_j(k)对第k-1时刻目标运动子模型j的状态向量进行一步预测，得到对应目标运动子模型的一步状态预测值

(4b)根据状态转移矩阵Φ_j(k)对第k-1时刻目标运动子模型j的状态协方差矩阵进行一步预测，得到对应目标运动子模型的一步协方差预测值P_j(k|k-1)：其中，Q_j表示模型j的过程噪声协方差；

(4c)获取第k时刻目标运动子模型j的量测矩阵H_j(k)，利用第k时刻目标运动子模型j的量测矩阵H_j(k)对目标运动子模型j的一步状态预测值进行量测转换，得到对应目标运动子模型的量测预测值 其中，H_j(k)表示雷达机动目标第k时刻目标运动子模型j的量测矩阵；

(4d)利用当前时刻观测值Z(k)和目标运动子模型j的量测预测值相减得到对应目标运动子模型的新息值V_j(k)：

(4e)由目标运动子模型j的一步协方差预测值P_j(k|k-1)和量测矩阵H_j(k)确定对应目标运动子模型的新息协方差矩阵S_j(k)：其中，R(k)是观测噪声在第k时刻的协方差矩阵；

(4f)根据目标运动子模型j的一步协方差预测值P_j(k|k-1)、量测矩阵H_j(k)以及新息协方差矩阵S_j(k)计算增益矩阵K_j(k)：(·)^H表示求矩阵的共轭转置，(·)^-1表示求矩阵的逆；

(4g)利用增益矩阵K_j(k)对目标运动子模型j的一步状态预测值和一步协方差预测值P_j(k|k-1)进行滤波，得到目标运动子模型j的第k时刻状态估计向量和第k时刻状态估计协方差矩阵P_j(k|k)：

${\hat{X}}_j\left(k\right|k)={\hat{X}}_j\left(k\right|k-1)+K_j\left(k\right)V_j\left(k\right)$

$P_j\left(k\right|k)=P_j\left(k\right|k-1)-K_j\left(k\right)S_j\left(k\right){K_j}^T\left(k\right).$

6.根据权利要求1所述的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，步骤5具体包括：

(5a)由目标运动子模型j的新息V_j(k)及其协方差矩阵S_j(k)计算新息的似然函数Λ_j(k)：

(5b)利用模型状态转移矩阵和新息的似然函数Λ_j(k)对目标运动子模型j的先验概率进行更新，得到目标运动子模型j的后验概率μ_j(k)：

${\mathrm{\&mu;}}_j\left(k\right)=\frac{{\mathrm{\&Lambda;}}_j\left(k\right)\underset{i=1}{\overset{r}{\&Sigma;}}{p}_{ij}(k-1){\mathrm{\&mu;}}_i(k-1)}{\underset{i=1}{\overset{r}{\&Sigma;}}{\mathrm{\&Lambda;}}_j\left(k\right)\underset{i=1}{\overset{r}{\&Sigma;}}{p}_{ij}(k-1){\mathrm{\&mu;}}_i(k-1)}.$

7.根据权利要求1所述的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，步骤6具体包括：

根据所述r个目标运动子模型第k时刻状态估计向量、第k时刻状态估计协方差矩阵以及第k时刻后验概率得到r个目标运动子模型交互后第k时刻状态估计向量和第k时刻状态估计协方差矩阵P(k|k)：

$\hat{X}\left(k\right|k)=\underset{j=1}{\overset{r}{\&Sigma;}}{\hat{X}}_j\left(k\right|k){\mathrm{\&mu;}}_j\left(k\right)$

$P\left(k\right|k)=\underset{j=1}{\overset{r}{\&Sigma;}}{\mathrm{\&mu;}}_j\left(k\right)\&CenterDot;\{P_j\left(k\right|k)+\&lsqb;{\hat{X}}_j\left(k\right|k)-\hat{X}\left(k\right|k)\&rsqb;\&CenterDot;{\&lsqb;{\hat{X}}_j\left(k\right|k)-\hat{X}\left(k\right|k)\&rsqb;}^T\}.$

8.根据权利要求1所述的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，步骤7具体包括：

若目标运动子模型数r为奇数即r＝2n+1，则在第k时刻的模型集为其特征参数分别对应r个目标运动子模型的角速度且在模型集中从小到大排列，称为中心模型；对应模型集的后验概率为同理若目标运动子模型数r为偶数即r＝2n，则在第k时刻的模型集为 为中心模型；对应模型集的后验概率为

按照下述方法调整模型间距：

第k时刻模型间距调整参数为 其中δ⁽ⁱ⁾表示模型区分距离，一般取

(7a)模型概率重心无跳变：此时则可根据其他模型概率判定模型集不变或向中心模型靠拢；

其中，t₁表示不可能模型的门限值，一般取t₁<0.1；

(7b)模型概率重心向左跳变：此时则可根据其他模型概率判定模型集不变或向左侧模型靠拢；

${\mathrm{\&omega;}}_{k+1}^{\max \left(i\right)}={\mathrm{\&omega;}}_{k+1}^{mid}+{\mathrm{\&lambda;}}_k^{r\left(i\right)}$

其中，t₂表示有效模型的门限值，一般取t₂>0.9；

(7c)模型概率重心向右跳变：此时则可根据其他模型概率判定模型集不变或向右侧模型靠拢；

${\mathrm{\&omega;}}_{k+1}^{\min \left(i\right)}={\mathrm{\&omega;}}_{k+1}^{mid}-{\mathrm{\&lambda;}}_k^{l\left(i\right)}.$

9.根据权利要求1所述的一种基于后验概率的变结构交互多模型雷达机动目标跟踪方法，其特征在于，步骤8具体包括：

(8a)对于目标运动子模型j，其后验概率在第k时刻与第k-1时刻的比值λ_j(k)为：

${\mathrm{\&lambda;}}_j\left(k\right)=\frac{{\mathrm{\&mu;}}_j\left(k\right)}{{\mathrm{\&mu;}}_j(k-1)}$

(8b)利用前后时刻后验概率比值λ_j(k)对第k-1时刻模型状态转移矩阵进行修正，第k时刻模型状态转移矩阵：

$p_{ij}\left(k\right)=\frac{{\mathrm{\&lambda;}}_j\left(k\right)p_{ij}(k-1)}{\underset{i=1}{\overset{r}{\&Sigma;}}{\mathrm{\&lambda;}}_j\left(k\right)p_{ij}(k-1)}.$