CN108896986A - 一种基于预测值的量测转换序贯滤波机动目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于预测值的量测转换序贯滤波机动目标跟踪方法，涉及雷达目标跟踪领域，主要针对机动目标跟踪情况下的多普勒雷达跟踪领域。本发明在交互多模型算法的框架下结合基于预测值的序贯滤波方法实现目标跟踪。其中，模型概率由位置模型概率和多普勒模型概率共同决定，提高了模型概率估计的准确性；同时，基于预测值的量测转换，消除了量测转换误差统计特性与量测值之间的相关性。该算法是一种有效的可同时处理多普勒量测的机动目标跟踪算法。

Description

一种基于预测值的量测转换序贯滤波机动目标跟踪方法

技术领域

本发明涉及雷达目标跟踪领域，主要针对机动目标跟踪情况下的多普勒雷达跟踪领域。

背景技术

将卡尔曼滤波(Kalman Filter,KF)算法运用到目标跟踪时，其状态方程通常用来描述目标的运动特性，往往是在直角坐标系中建立的，其状态向量主要由目标位置、速度及加速度等参量构成；在观测坐标系下，目标的量测信息通常表现为斜距、方位角及俯仰角。在绝大多数跟踪场景中，目标与观测者处于不同的坐标系时，导致了两者的非线性关系。为解决这类非线性滤波估计问题，最简单的方法是使用传统转换量测方式(ConventionalConversion)，即将极(球)坐标系下的量测值变换到直角坐标系下，使得目标状态与量测呈线性关系。然而在互距离(Cross-range)较大时，即目标位于远距离且量测角度误差较大的情况下。传统转换量测的方法会产生较大的偏差。因此，解决转换偏差所带来的跟踪性能下降的问题是非线性滤波中的首要工作。

为了消除非线性变换带来的偏差，许多学者在此提出了一些改进方法。Lerro D等(见文献Lerro D,Bar-Shalom Y.Tracking with debiased consistent convertedmeasurements versus EKF[J].Aerospace&Electronic Systems IEEE Transactions on,1993,29(3):1015-1022.)提出了一种加性去偏的量测转换方式(DCMKF)，即利用传统转换后的量测值去除转换误差，该方法在方位角度量测误差较小时有较好的去偏效果和量测转换一致性，结果表明其跟踪性能优于EKF算法。但由于量测转换误差的均值和方差是采用嵌套方式得到，Mo L.等提出一种无偏转换量测卡尔曼滤波(UCMKF)算法(见文献：Song X,Zhou Y,Bar-Shalom Y.Unbiased converted measurements for tracking[J].Aerospace&Electronic Systems IEEE Transactions on,1998,34(3):1023-1027.)，该方法采用乘性补偿方式去除偏差，弥补了传统转换方式的缺陷，然而进一步研究发现，该方法在量测转换误差均值和方差的推导上存在兼容性问题。修正的无偏转换量测卡尔曼滤波(MUCMKF)算法直接基于量测值对量测转换方差进行估计，解决了UCMKF的兼容性问题。但在推导量测误差协方差矩阵时，修正无偏量测转换卡尔曼滤波算法(MUCMKF)是在量测值已知的条件下进行的，这会导致量测误差协方差和量测值具有一定的相关性，因此其状态估计也是有偏的。对此进一步研究的相关学者提出了基于预测位置的去相关无偏量测转换卡尔曼滤波算法(DUCMKF)解决了以上算法存在的问题。(见文献：Bordonaro S,Willett P,Barshalom Y.Decorrelated unbiased converted measurement Kalman filter[J].IEEETransactions on Aerospace&Electronic Systems,2014,50(2):1431-1444.)

在目标跟踪领域内，多普勒雷达除了可以获得位置量测信息，还能够提供多普勒量测信息。理论计算与实践已经证明，处理位置量测信息的同时加入多普勒量测，能够更加有效地提升目标的跟踪精度(见文献：Farina A,Studer F A.Radar data processing I-Introduction and Tracking[J].Memorie Della Societa Astronomica Italiana,1985,56.)。以上量测转换算法都仅仅考虑了雷达的位置量测。先前有关多普勒量测的算法中，通常假设多普勒与距离之间的量测误差是统计独立的。然而Bar-Shalom在研究中发现，对于线性调频情况下的某些波形而言，距离和多普勒的量测误差之间存在一定的负相关性(见文献：Bar-Shalom Y.Negative correlation and optimal tracking with Dopplermeasurements[J].Aerospace&Electronic Systems IEEE Transactions on,2001,37(3):1117-1120.)。为了提升雷达的跟踪精度，又需要减小引入的多普勒量测时导致与位置之间的较强非线性，可将斜距和多普勒量测的乘积构造伪量测，来减小这种相关性。(见文献：ZhanshengDuan,Chongzhao Han and X.Rong Li,Sequential Nonlinear TrackingFilter with Range-rateMeasurements in Spherical Coordinates,7th InternationalConference on Information Fusion,Stockholm,2004,131-138.)段战胜等将多普勒量测信息引入到量测转换卡尔曼滤波算法中，将其推广为序贯滤波跟踪算法，但在求解转换量测误差统计特性时采用了嵌套的条件期望，且误差与量测具有相关性也会引入偏差。

当目标具有机动特性时，以上算法采用固定的系统状态方程无法描述目标的运动状态。交互多模型(IMM，interacting multiple model)算法能够有效的解决机动目标跟踪时存在的问题(见文献：Blom H AP,Barshalom Y.The interacting multiple modelalgorithm for systems with Markovian switching coefficients[J].IEEETransactions on Automatic Control,1988,33(8):780-783.)。该算法假设有限个目标模型存在，每个模型对应于不同的机动输入水平。在计算出各模型为正确的后验概率之后，就可以通过对各模型正确时的状态估计加权求和来给出最终的目标状态估计，其中的加权因子为模型正确的后验概率。针对机动目标情况，将基于预测值的序贯滤波算法与交互多模型算法相结合。本发明首先对量测值进行无偏转换处理，基于预测值推导转换误差统计特性，将序贯滤波算法与交互多模型算法框架相结合，利用基于线性KF算法的位置量测滤波器和基于非线性二阶EKF的伪量测滤波器所获得的模型概率平均值进行状态加权，从而获得一种基于预测值量测转换的序贯滤波机动目标跟踪算法。

发明内容

本发明针对背景技术的不足之处，考虑在机动目标跟踪情况下，利用多普勒雷达所提取的目标位置信息和多普勒信息进行序贯滤波，与交互多模型算法相结合，提出了一种针对机动目标的非线性量测转换序贯滤波跟踪算法。

本发明由混合去偏(位置量测采用相乘去偏方式；伪量测采用相减去偏方式)量测转换模块获得各模型的位置状态信息和多普勒信息的转换量测值，输入到滤波器输入端的交互模块，计算出交互后的各模型输入估计值，输出给各目标运动模型的滤波模块进行滤波：由位置量测滤波器产生的估计结果作为下一时刻的输入交互计算进行迭代，同时输出到伪量测滤波器中进行序贯处理，由伪量测滤波器产生的估计结果作为该模型的滤波估计。所述各目标运动模型滤波模块内的位置状态滤波器和伪量测滤波器，两者的新息计算结果共同输出给各模型概率计算模块，分别得到位置量测概率计算结果和伪量测概率计算结果后，再计算平均概率，作为当前时刻更新的模型概率。输出端的估计组合模块将各模型的滤波估计与平均概率进行加权求和计算得到目标位置状态，实现目标跟踪。

现对下文将用到的部分符号进行说明。目标在直角坐标系下的状态可描述为系统在x,y方向上的过程噪声假设为零均值高斯噪声，其协方差Q_k＝diag(q_x,q_y)，q表示各方向上的噪声方差；设多普勒雷达采样周期为T，获取的量测数据：和分别为距离量测、方位角量测和多普勒量测(下标k表示时刻点，上标m量测数据)；其对应的量测噪声和是零均值高斯白噪声序列，方差分别为和且和不相关，和的相关系数为ρ。

本发明的书方案为一种基于预测值的量测转换序贯滤波机动目标跟踪方法，该方法包括：

步骤1：计算量测转换值；

采用无偏量测转换方式计算计算方法如下：

其中，表示k时刻转换后的量测值，表示k时刻转换后的量测值在x方向上的分量，表示k时刻转换后的量测值在y方向上的分量，表示k时刻转换后伪量测值，为乘性偏差补偿因子；和分别为距离量测、方位角量测和多普勒量测，下标k表示时刻点，上标m量测数据；ρ相关系数，σ_r、距离和方位角的量测误差标准差；

步骤2：计算各模型的输入交互；

首先假设目标在k时刻运动模式表示为M_k，若多模型系统中存在N个模型，这N个模型能包含任何可能的目标运动模式，那么：

模型之间由模型i转移到模型j的转移概率为π_ij；

步骤2.1：计算混合概率

k时刻运动模型为j条件下，k-1时刻运动模型为i的混合概率为：

其中，为k-1时刻运动模型i的更新概率，归一化常数：

步骤2.2：计算各模型目标状态估计的交互输入；

通过对上一时刻估计值和协方差的交互计算得到每个模型滤波器输入的状态估计和协方差：

其中：和表示每个模型滤波器输入的状态估计和协方差；

步骤3：模型j的位置量测滤波估计；

基于预测值计算量测转换误差统计特性，并进行序贯滤波，包括以下几个子步骤：现为了公式表达简洁，将原处于每个符号表示当前目标状态为模型j的上标(j)省略，后续中每个符号均隐含有状态模型为j的含义；

步骤3.1:由卡尔曼一步预测方程计算预测状态

其中，分别表示目标在x和y方向上的预测位置，预测速度和预测加速度值；

步骤3.2：计算预测误差协方差P_k,k-1；

步骤3.3：计算极坐标系下的距离r_t、方位角θ_t和多普勒的预测值：

步骤3.4：计算预测值对应的预测误差方差

其中：表示距离的预测方差，0_2×1表示2行1列的零矩阵，表示方位角的预测方差，表示径向速度的预测方差；

步骤3.5：计算量测转换误差协方差R_k；

在上一步计算出预测距离、方位角、多普勒信息r_t、θ_t、及其预测误差方差条件下,计算去偏量测转换误差协方差R_k，下标k表示第k时刻；

其中，协方差矩阵各元素计算方式如下：

步骤3.6：分块量测协方差矩阵，构造位置量测协方差矩阵和伪量测互协方差矩阵；

步骤3.7：计算位置量测新息及新息误差协方差矩阵；

其中，为位置量测观测矩阵，表示k时刻转换后的位置量测值，表示位置的一步预测协方差矩阵；

步骤3.8：计算位置量测滤波器的增益；

步骤3.9：计算位置状态估计和状态估计误差协方差矩阵

步骤4：模型j的伪量测滤波估计；

为了公式表达的简洁性，将原处于每个符号表示当前目标状态为模型j的上标(j)省略，后续中每个符号均隐含有状态模型为j的含义；

步骤4.1：消除相关性；

去除位置量测与伪量测的相关性；将基于预测值求取的量测误差协方差矩阵进行分块：

其中，令

则有：

其中：L_k表示由量测误差协方差构造的分解向量，表示L_k中的第一列元素，表示L_k中的第二列元素，B_k表示由单位阵和L_k构成的分解矩阵，I₂表示2阶单位向量；

由矩阵的Cholesky分解得到：

其中，表示位置量测噪声，表示位置量测观测矩阵；

其中：分别表示k时刻的x,y方向的转换量测值和转换后的伪量测值， 分别表示k时刻的x,y方向的转换量测误差均值和伪量测转换误差均值；

步骤4.2：计算与伪量测信息相关的量测转换误差协方差

步骤4.3：这里使用二阶EKF对伪量测进行序贯滤波估计，计算伪量测新息及其协方差方法如下；

其中，表示与伪量测转换误差的均值，表示由位置估计结果变换得到的伪量测真实值，表示与伪量测相关的量测转换误差协方差矩阵，表示位置滤波的估计协方差，为在处的Jacobian矩阵；用来表示位置估计误差协方差矩阵的第i行第j列的元素，则

步骤4.4：计算伪量测滤波增益

步骤4.5：计算伪量测估计结果及估计误差协方差:

步骤4.6：最终模型j的滤波估计及误差自相关矩阵为：

步骤5：模型概率计算模块；

步骤5.1:计算位置量测滤波器似然函数；

其中：表示位置新息自相关矩阵，表示位置新息，上标j表示运动模型，p表示与位置量测相关的信息，下标k表示当前时刻；

步骤5.2:计算位置模型概率；

其中，C^(j),p表示模型j的位置归一化常数；

步骤5.3:计算伪量测滤波器似然函数；

其中，表示伪量测新息自相关矩阵，表示伪量测新息，上标j表示运动模型，ε表示与伪量测相关的信息，下标k表示当前时刻；

步骤5.4:计算伪量测模型概率；

其中，C^(j),ε表示模型j的伪量测归一化常数；

步骤5.5:计算运动模型j平均概率

步骤6：计算加权的状态估计和估计误差自相关矩阵P_k,k；

步骤7：重复步骤1-6进行递推滤波算法，实现目标跟踪。

从上述步骤可见，本发明在交互多模型算法的框架下结合基于预测值的序贯滤波方法实现目标跟踪。其中，模型概率由位置模型概率和多普勒模型概率共同决定，提高了模型概率估计的准确性；同时，基于预测值的量测转换，消除了量测转换误差统计特性与量测值之间的相关性。该算法是一种有效的可同时处理多普勒量测的机动目标跟踪算法。

在本发明中每个模型可看作是一个序贯滤波器，都包含有一个位置量测滤波器和伪量测滤波器。每个序贯滤波都有两个计算模型概率的部分，分别是位置滤波估计新息及其协方差矩阵计算位置信息的模型概率和多普勒量测滤波估计新息及其协方差矩阵计算多普勒信息的模型概率估计。为了提高模型概率的准确性，本发明分别计算位置信息和多普勒信息的模型概率估计后计算它们的均值。

附图说明

图1为本发明一种基于预测值的序贯滤波机动目标跟踪流程图；

图2为本发明交互多模型算法示意图；

图3为本发明具体实施方式中场景1 RMSE曲线；

图4为本发明具体实施方式中场景2 RMSE曲线；

图5为本发明具体实施方式中场景3 RMSE曲线。

其中，图3a，图4a，图5a是SQU-PRE和SQU-PRE-IMM的RMSE对比结果；图3b，图4b，图5b是SQ-COVN-IMM和SQU-PRE-IMM的RMSE对比结果。

具体实施方式

在此次实施方案中，考虑对二维平面内机动场景下的某一运动目标进行跟踪，目标初始值位置为(x＝5000m,y＝1000m)，初始速度为(v_x＝25m/s,v_y＝2m/s)。0～60s目标做匀速运动，60s～85s做匀加速运动，加速度为(-2m/s²,2m/s²)，85s～135s做匀速运动，135s～150s做匀加速运动，加速度为(3m/s²,-0.5m/s²)，150s～200s做匀速运动。径向速度和距离的相关系数为0.5。针对距离、方位角、多普勒量测噪声标准差取不同的值，做如下场景设置：

场景一：传感器的距离量测误差为方位角的量测误差为度，多普勒径向速度的量测误差为0.1m/s；

场景二：传感器的距离量测误差为方位角的量测误差为度，多普勒径向速度的量测误差为0.3m/s；

场景三：传感器的距离量测误差为方位角的量测误差为度，多普勒径向速度的量测误差为0.5m/s；

在上述3个场景下，采用基于传统量测转换序贯滤波的交互多模型方法(SQ-COVN-IMM)和基于预测值的序贯滤波方法(SQU-PRE)与本发明的基于预测值的量测转换序贯滤波机动目标跟踪方法(SQU-PRE-IMM)进行500次蒙特卡罗仿真。

图3、图4和图5分别采取上述三种机动目标跟踪方法的仿真结果的位置估计均方误差(RMSE)曲线。其中，图(3a)，(4a)，(5a)是SQU-PRE和SQU-PRE-IMM的RMSE对比结果，图(3b)，(4b)，(5b)是SQ-COVN-IMM和SQU-PRE-IMM的RMSE对比结果。

实施结果分析：一种基于预测值的序贯滤波机动目标跟踪方法可以实现带有多普勒量测的机动场景下的目标跟踪，并且跟踪精度比基于传统量测转换序贯滤波的交互多模型方法高；由于考虑了目标的机动特性，其跟踪精度也比没使用交互多模型的基于预测值的序贯滤波算法高。

三个场景下，根据SQU-PRE和SQU-PRE-IMM的对比结果可以看出，由于目标进行机动运动，因此SQU-PRE-IMM在交互多模型框架下进行目标跟踪，能获得较好的跟踪效果，而SQU-PRE在匀速运动阶段能够保持较好的跟踪效果，但从目标开始做匀加速运动时，滤波器中的模型与目标真实运动状态不相匹配，滤波误差大幅度增大，直至85秒时目标恢复匀速运动，滤波误差才开始下降，在135秒进行第二次匀加速运动时，曲线又出现发散情况。另一方面，对比SQ-CONV-IMM可以看出，SQU-PRE-IMM具有较好的跟踪性能和更高的跟踪精度，三个场景下随着量测误差的增加更加明显。

综上所述，基于预测值的序贯滤波机动目标跟踪算法SQU-PRE-IMM可以有效的实现机动目标跟踪，且具有较高的跟踪精度。

Claims (1)

1.一种基于预测值的量测转换序贯滤波机动目标跟踪方法，该方法包括：

步骤1：计算量测转换值；

采用无偏量测转换方式计算计算方法如下：

其中，表示k时刻转换后的量测值，表示k时刻转换后的量测值在x方向上的分量，表示k时刻转换后的量测值在y方向上的分量，表示k时刻转换后伪量测值，为乘性偏差补偿因子；和分别为距离量测、方位角量测和多普勒量测，下标k表示时刻点，上标m量测数据；ρ相关系数，σ_r、距离和方位角的量测误差标准差；

步骤2：计算各模型的输入交互；

首先假设目标在k时刻运动模式表示为M_k，若多模型系统中存在N个模型，这N个模型能包含任何可能的目标运动模式，那么：

模型之间由模型i转移到模型j的转移概率为π_ij；

步骤2.1：计算混合概率

k时刻运动模型为j条件下，k-1时刻运动模型为i的混合概率为：

其中，为k-1时刻运动模型i的更新概率，归一化常数：

步骤2.2：计算各模型目标状态估计的交互输入；

通过对上一时刻估计值和协方差的交互计算得到每个模型滤波器输入的状态估计和协方差：

其中：和表示每个模型滤波器输入的状态估计和协方差；

步骤3：模型j的位置量测滤波估计；

基于预测值计算量测转换误差统计特性，并进行序贯滤波，包括以下几个子步骤：现为了公式表达简洁，将原处于每个符号表示当前目标状态为模型j的上标(j)省略，后续中每个符号均隐含有状态模型为j的含义；

步骤3.1:由卡尔曼一步预测方程计算预测状态

其中，分别表示目标在x和y方向上的预测位置，预测速度和预测加速度值；

步骤3.2：计算预测误差协方差P_k,k-1；

步骤3.3：计算极坐标系下的距离r_t、方位角θ_t和多普勒的预测值：

步骤3.4：计算预测值对应的预测误差方差

其中：表示距离的预测方差，0_2×1表示2行1列的零矩阵，表示方位角的预测方差，表示径向速度的预测方差；

步骤3.5：计算量测转换误差协方差R_k；

在上一步计算出预测距离、方位角、多普勒信息r_t、θ_t、及其预测误差方差条件下,计算去偏量测转换误差协方差R_k，下标k表示第k时刻；

其中，协方差矩阵各元素计算方式如下：

步骤3.6：分块量测协方差矩阵，构造位置量测协方差矩阵和伪量测互协方差矩阵；

步骤3.7：计算位置量测新息及新息误差协方差矩阵；

其中，为位置量测观测矩阵，表示k时刻转换后的位置量测值，表示位置的一步预测协方差矩阵；

步骤3.8：计算位置量测滤波器的增益；

步骤3.9：计算位置状态估计和状态估计误差协方差矩阵

步骤4：模型j的伪量测滤波估计；

为了公式表达的简洁性，将原处于每个符号表示当前目标状态为模型j的上标(j)省略，后续中每个符号均隐含有状态模型为j的含义；

步骤4.1：消除相关性；

去除位置量测与伪量测的相关性；将基于预测值求取的量测误差协方差矩阵进行分块：

其中，令

则有：

其中：L_k表示由量测误差协方差构造的分解向量，表示L_k中的第一列元素，表示L_k中的第二列元素，B_k表示由单位阵和L_k构成的分解矩阵，I₂表示2阶单位向量；

由矩阵的Cholesky分解得到：

其中，表示位置量测噪声，表示位置量测观测矩阵；

其中：分别表示k时刻的x,y方向的转换量测值和转换后的伪量测值 分别表示k时刻的x,y方向的转换量测误差均值和伪量测转换误差均值；

步骤4.2：计算与伪量测信息相关的量测转换误差协方差

步骤4.3：这里使用二阶EKF对伪量测进行序贯滤波估计，计算伪量测新息及其协方差方法如下；

其中，表示与伪量测转换误差的均值，表示由位置估计结果变换得到的伪量测真实值，表示与伪量测相关的量测转换误差协方差矩阵，表示位置滤波的估计协方差，为在处的Jacobian矩阵；用来表示位置估计误差协方差矩阵的第i行第j列的元素，则

步骤4.4：计算伪量测滤波增益

步骤4.5：计算伪量测估计结果及估计误差协方差:

步骤4.6：最终模型j的滤波估计及误差自相关矩阵为：

步骤5：模型概率计算模块；

步骤5.1:计算位置量测滤波器似然函数；

其中：表示位置新息自相关矩阵，表示位置新息，上标j表示运动模型，p表示与位置量测相关的信息，下标k表示当前时刻；

步骤5.2:计算位置模型概率；

其中，C^(j),p表示模型j的位置归一化常数；

步骤5.3:计算伪量测滤波器似然函数；

其中，表示伪量测新息自相关矩阵，表示伪量测新息，上标j表示运动模型，ε表示与伪量测相关的信息，下标k表示当前时刻；

步骤5.4:计算伪量测模型概率；

其中，C^(j),ε表示模型j的伪量测归一化常数；

步骤5.5:计算运动模型j平均概率

步骤6：计算加权的状态估计和估计误差自相关矩阵P_k,k；

步骤7：重复步骤1-6进行递推滤波算法，实现目标跟踪。