CN111965618B - 一种融合多普勒量测的转换量测跟踪方法及系统 - Google Patents

一种融合多普勒量测的转换量测跟踪方法及系统 Download PDF

Info

Publication number
CN111965618B
CN111965618B CN202010833608.8A CN202010833608A CN111965618B CN 111965618 B CN111965618 B CN 111965618B CN 202010833608 A CN202010833608 A CN 202010833608A CN 111965618 B CN111965618 B CN 111965618B
Authority
CN
China
Prior art keywords
measurement
matrix
coordinate system
target
covariance matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010833608.8A
Other languages
English (en)
Other versions
CN111965618A (zh
Inventor
盛琥
任清安
杨秀凯
单齐
庞博清
陈磊
唐思圆
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
CETC 38 Research Institute
Original Assignee
CETC 38 Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by CETC 38 Research Institute filed Critical CETC 38 Research Institute
Priority to CN202010833608.8A priority Critical patent/CN111965618B/zh
Publication of CN111965618A publication Critical patent/CN111965618A/zh
Application granted granted Critical
Publication of CN111965618B publication Critical patent/CN111965618B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Abstract

本发明公开了一种融合多普勒量测的转换量测跟踪方法及系统，属于雷达目标跟踪技术领域，包括步骤：S1：初始化目标状态；S2：更新视线坐标系；S3：构造转换量测模型；S4：估计目标状态；S5：下一时刻目标状态估计。为了克服多普勒量测的强非线性，在视线坐标系内对目标跟踪；在每个滤波时刻，追踪雷达到目标视轴的变化并更新视线坐标系，将目标状态、雷达量测转换到更新后的视线坐标系内；视线坐标系内的转换量测和目标状态间满足线性关系，从而提高滤波性能，减小计算量。该融合方法可以扩展到带多普勒的三坐标雷达的跟踪应用中。发明克服了多普勒量测和目标状态间的强非线性，能够以较小的计算量对目标跟踪，在估计精度和置信度都具有较好性能。

Description

一种融合多普勒量测的转换量测跟踪方法及系统
技术领域
本发明涉及雷达目标跟踪技术领域，具体涉及一种融合多普勒量测的转换量测跟踪方法及系统。
背景技术
雷达目标跟踪中，目标状态模型一般在直角坐标系下建立，而量测模型是在极/球坐标系下获得，状态和量测的不兼容产生非线性估计问题。为解决该问题，扩展卡尔曼滤波
(EKF)、无迹滤波(UKF)、容积卡尔曼滤波(CKF)、粒子滤波(PF)、转换量测卡尔曼滤波(CMKF)等多种方法被相继提出。其中的CMKF方法因为精度较高、实时性好、实现简单等特点，在实际中得到广泛应用。CMKF的思想是先将极/球坐标下的非线性量测变换为直角坐标系内的伪线性表达，再用标准卡尔曼滤波器跟踪，具有较好综合性能。
理论计算和实践证明：利用多普勒信息可以有效提高目标跟踪精度。由于多普勒和目标状态的强非线性关系，EKF等方法的处理效果很差。有学者提出用斜距和多普勒量测的乘积构造伪量测，减弱多普勒量测和目标状态的非线性程度，再用CMKF方法序贯滤波跟踪的思想，该方法需要对位置和伪量测去相关，再用位置量测和伪量测分别序贯滤波，计算复杂度较高。另有学者提出基于转换多普勒量测的最优线性无偏滤波器，克服了CMKF方法的固有缺陷，但是该方法的滤波参数的表达式比较复杂。还有学者提出基于静态融合的CMKF方法，用位置量测和伪量测分别独立滤波，再通过静态融合去除二者的相关性，取得满意的结果。但该方法需要两个滤波器和一个静态估计器，计算复杂度高。为克服以上方法缺点，提出一种融合多普勒量测的转换量测跟踪方法及系统。
发明内容
本发明所要解决的技术问题在于：克服多普勒量测和目标状态间的强非线性，以较小的计算量对目标状态精确估计，提供一种融合多普勒量测的转换量测跟踪方法。
本发明是通过以下技术方案解决上述技术问题的，本发明包括以下步骤：
S1：初始化目标状态
在k＝0时刻，获得雷达量测；在k＝1时刻，再次获取雷达量测，设雷达对目标均匀采样，采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure GDA0003800364560000011
和初始协方差阵
Figure GDA0003800364560000012
，以及直角坐标系到视线坐标系的初始转换矩阵M1；
S2：更新视线坐标系
在k＞1时刻，收到雷达各量测后，预测目标状态及其协方差阵，更新视线坐标系；
S3：构造转换量测模型
在步骤S2更新后的视线坐标系内，基于步骤S2中更新后的目标预测及其协方差阵，判断是否能够引入多普勒量测，并根据判断结果构造含或不含多普勒量测的转换量测，估计转换量测的协方差阵以及量测系数矩阵；
S4：估计目标状态
根据步骤S2中更新后的目标预测及其协方差阵，结合步骤S3中得到的转换量测、转换量测的协方差阵以及量测系数矩阵，估计目标状态，根据需要输出目标在直角坐标系的状态估计和估计协方差阵；
S5：下一时刻目标状态估计
进入下一个采样时刻，重复步骤S2～S4，递推估计目标在视线坐标系内的状态。
更进一步地，在所述步骤S1中，在时刻k，目标状态方程和雷达量测方程分别如下：
Xk＝FkXk-1+Gkvk
zk＝HkXk+wk
其中，
Figure GDA0003800364560000021
为k时刻目标状态，xk、
Figure GDA0003800364560000022
是沿视轴方向的位置、速度分量，yk、
Figure GDA0003800364560000023
是垂直视轴方向的位置、速度分量；
Fk为k时刻状态转移矩阵，状态方程采用近常速运动(NCV：Nearly ConstantVelocity)模型时，其表达式如下：
Figure GDA0003800364560000024
Gk为k时刻噪声输入矩阵，目标近常速运动时的表达式如下：
Figure GDA0003800364560000025
系统噪声为vk＝[vx vy]T，vx、vy分别是沿着视轴方向和垂直视轴方向的零均值高斯过程噪声，其协方差阵为Qk，zk为转换量测值，Hk为量测系数矩阵，wk为量测噪声向量，其协方差阵为Rk。
更进一步地，雷达量测分别为斜距量测rm，方位量测θm和多普勒量测
Figure GDA0003800364560000031
量测噪声是零均值高斯白噪声，各量测的标准差为σr、σθ、
Figure GDA0003800364560000032
斜距量测和多普勒量测的相关系数为ρ。
更进一步地，所述步骤S1的具体过程如下：
S11：在视线坐标系内估计目标初始状态
Figure GDA0003800364560000033
目标初始状态
Figure GDA0003800364560000034
为：
Figure GDA0003800364560000035
S12：估计视线坐标系内的初始协方差阵
Figure GDA0003800364560000036
初始协方差阵
Figure GDA0003800364560000037
为：
Figure GDA0003800364560000038
S13：估计初始转换矩阵M1，初始转换矩阵M1为：
Figure GDA0003800364560000039
其中，
Figure GDA00038003645600000310
更进一步地，所述步骤S2的具体过程如下：
S21：估计目标状态一步预测及其协方差阵，分别为：
Figure GDA00038003645600000311
Figure GDA00038003645600000312
其中，
Figure GDA00038003645600000313
是状态一步预测值，xp,1、yp,1是位置一步预测分量，
Figure GDA00038003645600000314
是速度一步预测分量，Fk是k时刻状态转移矩阵，
Figure GDA00038003645600000315
是k-1时刻状态估计；
Figure GDA00038003645600000316
是状态预测协方差阵；
Figure GDA00038003645600000317
是k-1时刻的状态估计协方差阵；Qk是系统噪声vk的协方差阵，Gk为k时刻噪声输入矩阵；
S22：估计视轴旋转矩阵Tk，更新转换矩阵Mk
先根据步骤S21的结果，估计
Figure GDA0003800364560000041
的位置向量相对视轴的转动角度βk，转动角度βk为：
Figure GDA0003800364560000042
再估计视轴旋转矩阵Tk，视轴旋转矩阵Tk为：
Figure GDA0003800364560000043
其中，
Figure GDA0003800364560000044
最后更新转换矩阵Mk，Mk为直角坐标系到视线坐标系的转换矩阵，具体为：
Mk＝TkMk-1；
更新转换矩阵Mk后，视线坐标系的视轴方向改变，k-1时刻的视线坐标系被更新为k时刻的视线坐标系；
S23：在k时刻视线坐标系内更新目标预测及其协方差阵，分别为：
Figure GDA0003800364560000045
Figure GDA0003800364560000046
更新后的目标预测为
Figure GDA0003800364560000047
xp、yp是位置预测分量，
Figure GDA0003800364560000048
是速度预测分量。
更进一步地，所述步骤S3的具体过程如下：
S31：基于步骤S23计算得到的
Figure GDA0003800364560000049
和
Figure GDA00038003645600000410
判断是否能够引入多普勒量测，如果
Figure GDA00038003645600000411
说明多普勒信息可用，执行步骤S32，否则执行步骤S33；
Figure GDA00038003645600000412
的表达式如下所示；
Figure GDA00038003645600000413
其中
Figure GDA00038003645600000414
是
Figure GDA00038003645600000415
第三行、第三列的元素。
Figure GDA00038003645600000416
S32：构造含多普勒量测的转换量测zk，估计其协方差阵Rk及量测系数矩阵Hk，其中：
Figure GDA0003800364560000051
Mk(1:3,1:4)是步骤S22中估计出的Mk的第1～3行、第1～4列的所有元素；
基于步骤S23得到的xp、
Figure GDA0003800364560000052
和
Figure GDA0003800364560000053
估计转换量测zk的协方差阵Rk，协方差阵Rk为：
Figure GDA0003800364560000054
其中：
Figure GDA0003800364560000055
Figure GDA0003800364560000056
Figure GDA0003800364560000057
Figure GDA0003800364560000058
Figure GDA0003800364560000059
是
Figure GDA00038003645600000510
的第一行、第一列的元素；
估计量测系数矩阵Hk，量测系数矩阵Hk为：
Figure GDA00038003645600000511
S33：构造无多普勒量测的转换量测zk，估计其协方差阵Rk及量测系数矩阵Hk，其中：
zk＝Mk(1:2,1:2)[rmcosθm rmsinθm]T
Mk(1:2,1:2)是步骤S22估计出的Mk的第1～2行、第1～2列的所有元素；
基于步骤S23得到的xp和
Figure GDA00038003645600000512
估计转换量测zk的协方差阵Rk，协方差阵Rk为：
Figure GDA00038003645600000513
估计量测系数矩阵Hk，量测系数矩阵Hk为：
Figure GDA00038003645600000514
其中：
Figure GDA00038003645600000515
更进一步地，在所述步骤S31中，κ是门限因子，取值不小于5。
更进一步地，所述步骤S4的具体过程如下：
S41：计算转换量测预测残差
Figure GDA0003800364560000061
Figure GDA0003800364560000062
S42：估计转换量测预测残差
Figure GDA0003800364560000063
的协方差阵Sk
Figure GDA0003800364560000064
S43：估计转换量测滤波增益Kk
Figure GDA0003800364560000065
S44：计算转换量测滤波器状态估计
Figure GDA0003800364560000066
Figure GDA0003800364560000067
S45：计算转换量测滤波估计协方差阵
Figure GDA0003800364560000068
Figure GDA0003800364560000069
更进一步地，在所述步骤S5中，直角坐标系内的状态估计
Figure GDA00038003645600000610
和估计协方差阵
Figure GDA00038003645600000611
为：
Figure GDA00038003645600000612
Figure GDA00038003645600000613
本发明还提供了一种融合多普勒量测的转换量测跟踪系统，包括：
初始化模块，用于在k＝0时刻，获得雷达量测；在k＝1时刻，再次获取雷达量测，设雷达对目标均匀采样，采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure GDA00038003645600000614
和初始协方差阵
Figure GDA00038003645600000615
以及直角坐标系到视线坐标系的初始转换矩阵M1；
更新模块，用于在k＞1时刻，收到雷达各量测后，预测目标状态及其协方差阵，更新视线坐标系；
模型构造模块，用于在更新后的视线坐标系内，基于更新后的目标预测及其协方差阵，判断是否能够引入多普勒量测，并根据判断结果构造含或不含多普勒量测的转换量测，估计转换量测的协方差阵以及量测系数矩阵；
第一估计模块，用于更新后的目标预测及其协方差阵，结合得到的转换量测、转换量测协方差阵以及量测系数矩阵，估计目标在视线坐标系/直角坐标系的状态估计和协方差阵；
第二估计模块，用于在进入下一个采样时刻时重复步骤S2～S4，递推估计目标在视线坐标系/直角坐标系内的状态估计和协方差阵；
控制处理模块，用于向各模块发出指令，完成相关动作；
所述初始化模块、更新模块、模型构造模块、第一估计模块、第二估计模块均与控制处理模块电连接。
本发明相比现有技术具有以下优点：该融合多普勒量测的转换量测跟踪方法，通过在视线坐标系内对目标跟踪，克服了多普勒量测和目标状态间的强非线性，能够以较小的计算量估计目标状态；在估计精度和置信度都具有较好性能，值得被推广使用。
附图说明
图1是本发明实施例一中转换量测跟踪方法的流程示意图；
图2a是本发明实施例二中场景1下各方法位置误差对比图；
图2b是本发明实施例二中场景1下各方法ANEES对比图；
图3a是本发明实施例二中场景2下各方法位置误差对比图；
图3b是本发明实施例二中场景2下各方法ANEES对比图。
具体实施方式
下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
实施例一
本实施例提供一种技术方案：一种融合多普勒量测的转换量测跟踪方法，为了克服多普勒量测的强非线性，在视线坐标系内对目标跟踪。在每个滤波时刻，追踪雷达到目标视轴的变化，将目标状态、雷达量测转换到视线坐标系内。在视线坐标系内，转换量测和目标状态之间满足线性关系，能够对目标高精度跟踪。该方法可以扩展到带多普勒的三坐标雷达的跟踪应用中，本方法主要包括四部分：第一部分，初始化目标状态；第二部分，预测目标状态，追踪并更新视线坐标系；第三部分，在当前视线坐标系内构造转换量测；第四部分，在视线坐标系内，融合多普勒量测，估计目标状态。
具体包括如下技术内容：
假设一部具备斜距、方位和多普勒量测的雷达，对目标进行融合跟踪，雷达位于视线坐标系的坐标原点。在k时刻，目标状态方程和雷达量测方程分别如下：
Xk＝FkXk-1+Gkvk
zk＝HkXk+wk
其中
Figure GDA0003800364560000081
为k时刻目标状态，xk、
Figure GDA0003800364560000082
是沿视轴方向的位置、速度分量，yk、
Figure GDA0003800364560000083
是垂直视轴方向的位置、速度分量。
Fk为k时刻状态转移矩阵，状态方程采用近常速运动(NCV：Nearly ConstantVelocity)模型，其表达式如下：
Figure GDA0003800364560000084
T为雷达均匀采样的间隔。Gk为k时刻噪声输入矩阵，目标采取近常速运动时的表达式如下：
Figure GDA0003800364560000085
系统噪声为vk＝[vx vy]T，vx、vy分别是沿着视轴方向和垂直视轴方向的零均值高斯过程噪声，其协方差阵为Qk，zk为转换量测值，Hk为量测系数矩阵，wk为量测噪声向量，其协方差阵为Rk。
雷达量测分别为斜距量测rm，方位量测θm和多普勒量测
Figure GDA0003800364560000086
量测噪声是零均值高斯白噪声，标准差为σr、σθ、
Figure GDA0003800364560000087
和
Figure GDA0003800364560000088
的相关系数为ρ。
基于以下步骤实施对目标的跟踪：
步骤1：滤波器初始化
在k＝0时刻，获得量测rm,0、θm,0、
Figure GDA0003800364560000089
在k＝1时刻，获得量测rm,1、θm,1、
Figure GDA00038003645600000810
假设雷达对目标均匀采样，其采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure GDA00038003645600000811
和初始协方差阵
Figure GDA00038003645600000812
以及直角坐标系到视线坐标系的初始转换矩阵M1；
步骤1.1：在视线坐标系内估计初始状态
Figure GDA00038003645600000813
Figure GDA00038003645600000814
步骤1.2：估计视线坐标系内的初始协方差阵
Figure GDA0003800364560000091
Figure GDA0003800364560000092
步骤1.3：估计初始转换矩阵M1：
Figure GDA0003800364560000093
其中
Figure GDA0003800364560000094
步骤2：在k＞1时刻，收到雷达量测后，预测目标状态及其协方差阵，更新视线坐标系；
步骤2.1：估计目标状态一步预测及其协方差阵：
Figure GDA0003800364560000095
Figure GDA0003800364560000096
其中，
Figure GDA0003800364560000097
是状态一步预测值，xp,1、yp,1是位置一步预测分量，
Figure GDA0003800364560000098
是速度一步预测分量，Fk是k时刻状态转移矩阵，
Figure GDA0003800364560000099
是k-1时刻状态估计；
Figure GDA00038003645600000910
是状态一步预测协方差阵；
Figure GDA00038003645600000911
是k-1时刻的状态估计协方差阵；Qk是系统噪声vk的协方差阵，Gk为k时刻噪声输入矩阵。
步骤2.2：估计视轴旋转矩阵Tk，更新转换矩阵Mk；
首先根据步骤2.1的结果，估计
Figure GDA00038003645600000912
的位置向量相对视轴的转动角度βk：
Figure GDA00038003645600000913
接着估计视轴旋转矩阵Tk：
Figure GDA00038003645600000914
其中，
Figure GDA00038003645600000915
最后更新转换矩阵Mk，Mk是直角坐标系到视线坐标系的转换矩阵；
Mk＝TkMk-1
更新转换矩阵Mk后，视线坐标系的视轴方向改变，k-1时刻的视线坐标系被更新为k时刻的视线坐标系；
步骤2.3：在k时刻视线坐标系内更新目标预测及其协方差阵：
Figure GDA0003800364560000101
Figure GDA0003800364560000102
更新后的目标预测为
Figure GDA0003800364560000103
xp、yp是其位置预测分量，
Figure GDA0003800364560000104
是其速度预测分量。
步骤3：构造转换量测模型。
步骤3.1：基于步骤2.3计算得到的
Figure GDA0003800364560000105
和
Figure GDA0003800364560000106
判断是否能够引入多普勒量测，如果
Figure GDA0003800364560000107
说明多普勒信息可用，执行步骤3.2，否则执行步骤3.3；κ是门限因子，取值不小于5；
Figure GDA0003800364560000108
的表达式如下所示；
Figure GDA0003800364560000109
其中
Figure GDA00038003645600001010
是
Figure GDA00038003645600001011
第三行、第三列的元素，
Figure GDA00038003645600001012
是沿着视轴的速度预测，
Figure GDA00038003645600001013
Figure GDA00038003645600001014
需要说明的是，
Figure GDA00038003645600001015
的分子部分是多普勒转换量测的方差，分母部分是沿视轴的速度预测方差。其含义是：多普勒转换量测相对于沿视轴的速度预测是非线性的，为了保证滤波的鲁棒性，多普勒转换量测方差要数倍于视轴速度预测方差。
步骤3.2：构造含多普勒量测的转换量测zk，估计其协方差阵Rk以及量测系数矩阵Hk；
Figure GDA00038003645600001016
其中，Mk(1:3,1:4)是步骤2.2估计出的Mk的第1～3行、第1～4列的所有元素。
基于步骤2.3得到的xp、
Figure GDA00038003645600001017
和
Figure GDA00038003645600001018
估计转换量测zk的协方差阵Rk：
Figure GDA00038003645600001019
其中：
Figure GDA0003800364560000111
Figure GDA0003800364560000112
Figure GDA0003800364560000113
Figure GDA0003800364560000114
Figure GDA0003800364560000115
是
Figure GDA0003800364560000116
的第一行、第一列的元素；
估计量测系数矩阵Hk：
Figure GDA0003800364560000117
步骤3.3：构造无多普勒量测的转换量测zk，估计其协方差阵Rk，以及量测系数矩阵Hk：
zk＝Mk(1:2,1:2)[rmcosθm rmsinθm]T
其中，Mk(1:2,1:2)是步骤2.2估计出的Mk的第1～2行、第1～2列的所有元素。
基于步骤2.3得到的xp和
Figure GDA0003800364560000118
估计转换量测zk的协方差阵Rk：
Figure GDA0003800364560000119
估计量测系数矩阵Hk：
Figure GDA00038003645600001110
其中：
Figure GDA00038003645600001111
Figure GDA00038003645600001112
步骤4：基于步骤2.3得到的
Figure GDA00038003645600001113
以及步骤3.2或步骤3.3得到的zk、Rk和Hk，估计目标状态。
步骤4.1：计算转换量测预测残差
Figure GDA00038003645600001114
Figure GDA00038003645600001115
步骤4.2：估计转换量测预测残差
Figure GDA00038003645600001116
的协方差阵Sk：
Figure GDA00038003645600001117
其中，Hk是量测系数矩阵，
Figure GDA00038003645600001118
是状态预测协方差阵，Rk是转换量测协方差阵。
步骤4.3：估计转换量测滤波增益Kk：
Figure GDA0003800364560000121
步骤4.4：计算转换量测滤波器状态估计：
Figure GDA0003800364560000122
步骤4.5：计算转换量测滤波估计协方差阵：
Figure GDA0003800364560000123
步骤5：下一个采样时刻，重复步骤2～4，递推估计目标在视线坐标系内的状态，用以下公式输出目标在直角坐标系的状态估计
Figure GDA0003800364560000124
和估计协方差阵
Figure GDA0003800364560000125
Figure GDA0003800364560000126
Figure GDA0003800364560000127
实施例二
在本实施例中，考虑对极坐标下的两种目标跟踪场景进行仿真。
设雷达位于原点，量测噪声为零均值高斯分布，标准差为σr＝10m，σθ＝1°，
Figure GDA0003800364560000128
目标初始位置为(0,50)km，全程做近常速运动，速度为(50,0)m/s。雷达采样间隔为1s，蒙特卡洛仿真次数100次，仿真时长200s。
场景1：目标沿直角坐标系各轴的过程噪声标准差为0.1m/s2，斜距和多普勒量测的相关系数ρ＝0.1。
场景2：目标沿直角坐标系各轴的过程噪声标准差为1m/s2，斜距和多普勒量测的相关系数ρ＝-0.1。
图2a、图2b和图3a、图3b分别为本发明所提算法在本实施例的2种场景下，和静态融合CMKF方法及序贯融合CMKF方法对目标跟踪的精度对比图。
通过对转换量测跟踪算法的分析，我们选取静态融合CMKF方法(见文献:Statically Fused Converted Position and Doppler Measurement Kalman Filters[J].IEEE Transactions on AES,2014,50(1):300-318)和序贯融合CMKF方法(见文献：极坐标系中带多普勒量测的雷达目标跟踪[J].系统仿真学报,2004,16(12):2860-2863)与本发明所提方法进行比较。
所有算法在均使用相同的目标初始状态，采用两点差分法起始，目标跟踪性能指标包括平均归一化估计误差平方(ANEES)和位置均方根误差(RMSE)，具体定义如下：
Figure GDA0003800364560000131
Figure GDA0003800364560000132
其中，
Figure GDA0003800364560000133
和
Figure GDA0003800364560000134
是第i次仿真时x方向和y方向的状态估计误差，N为仿真次数。RMSE越小，说明算法跟踪精度越高。ANEES可以反映滤波器估计的可信程度，ANEES为1时，表明滤波实际误差和估计误差协方差完全一致，置信度最高。
图2a、图2b和图3a、图3b是场景1和场景2下，雷达采用不同方法，对目标的融合跟踪性能。其中图2a、图3a是各方法位置精度对比，图2b、图3b是各方法的ANEES对比。通过与其他方法的对比可见：本发明方法的位置精度和其他方法相当，但估计置信度更高(逼近1)，因此综合性能更优。
表1是三种方法在两个场景中仿真100次的耗时对比。由表1可见，静态融合CMKF需要两个递归滤波器和一个静态融合器同时工作，因此耗时最多；本发明方法和序贯融合CMKF方法的计算量在一个量级。
表1、200次蒙特卡洛仿真耗时对比
Figure GDA0003800364560000135
从仿真结果可见，在两种不同场景中，本发明所提方法的置信度均高于静态融合CMKF和序贯融合CMKF方法，位置精度与其他方法相当，计算量小于静态融合CMKF。综合以上，本发明方法在估计精度、计算复杂度和置信度都有不俗表现，有较好的实用意义。
综上所述，上述实施例的融合多普勒量测的转换量测跟踪方法，通过在视线坐标系内对目标跟踪，克服了多普勒量测和目标状态间的强非线性，能够以较小的计算量估计目标状态，在估计精度和置信度指标上表现良好，具有应用推广价值。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (4)

1.一种融合多普勒量测的转换量测跟踪方法，其特征在于，包括以下步骤：
S1：初始化目标状态
在k＝0时刻，获得雷达量测；在k＝1时刻，再次获取雷达量测，设雷达对目标均匀采样，采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure FDA0003782282440000011
和初始协方差阵
Figure FDA0003782282440000012
以及直角坐标系到视线坐标系的初始转换矩阵M1；
S2：更新视线坐标系
在k＞1时刻，收到雷达各量测后，预测目标状态及其协方差阵，更新视线坐标系；
S3：构造转换量测模型
在步骤S2更新后的视线坐标系内，基于步骤S2中更新后的目标预测及其协方差阵，判断是否能够引入多普勒量测，并根据判断结果构造含或不含多普勒量测的转换量测，估计转换量测的协方差阵以及量测系数矩阵；
S4：估计目标状态
根据步骤S2中更新后的目标预测及其协方差阵，结合步骤S3中得到的转换量测、转换量测协方差阵以及量测系数矩阵，估计目标在视线坐标系/直角坐标系的状态估计和协方差阵；
S5：下一时刻目标状态估计
进入下一个采样时刻，重复步骤S2～S4，递推估计目标在视线坐标系/直角坐标系内的状态估计和协方差阵；
在所述步骤S1中，在时刻k，目标状态方程和雷达量测方程分别如下：
Xk＝FkXk-1+Gkvk
zk＝HkXk+wk
其中，
Figure FDA0003782282440000013
为k时刻目标状态，xk、
Figure FDA0003782282440000014
是沿视轴方向的位置、速度分量，yk、
Figure FDA0003782282440000015
是垂直视轴方向的位置、速度分量；
Fk为k时刻状态转移矩阵，状态方程采用近常速模型，其表达式如下：
Figure FDA0003782282440000016
其中，T是雷达均匀采样的间隔，Gk为k时刻噪声输入矩阵，目标近常速运动的表达式如下：
Figure FDA0003782282440000021
系统噪声为vk＝[vx vy]T，vx、vy分别是沿着视轴方向和垂直视轴方向的零均值高斯过程噪声，其协方差阵为Qk，zk为转换量测值，Hk为量测系数矩阵，wk为量测噪声向量，其协方差阵为Rk；
雷达量测分别为斜距量测rm，方位量测θm和多普勒量测
Figure FDA0003782282440000022
量测噪声是零均值高斯噪声，各量测的标准差为σr、σθ、σr，斜距量测和多普勒量测的相关系数为ρ；
所述步骤S1的具体过程如下：
S11：在视线坐标系内估计目标初始状态
Figure FDA0003782282440000023
目标初始状态
Figure FDA0003782282440000024
为：
Figure FDA0003782282440000025
S12：估计视线坐标系内的初始协方差阵
Figure FDA0003782282440000026
初始协方差阵
Figure FDA0003782282440000027
为：
Figure FDA0003782282440000028
S13：估计初始转换矩阵M1，初始转换矩阵M1为：
Figure FDA0003782282440000029
其中，
Figure FDA00037822824400000210
02×2是2行2列的零矩阵；
所述步骤S2的具体过程如下：
S21：估计目标状态一步预测值及其协方差阵，分别为：
Figure FDA00037822824400000211
Figure FDA0003782282440000031
其中，
Figure FDA0003782282440000032
是k时刻的状态一步预测值，xp,1、yp,1是位置一步预测分量，
Figure FDA0003782282440000033
是速度一步预测分量，Fk是k时刻状态转移矩阵，
Figure FDA0003782282440000034
是k-1时刻状态估计；
Figure FDA0003782282440000035
是k时刻状态一步预测协方差阵；
Figure FDA0003782282440000036
是k-1时刻的状态估计协方差阵；Qk是k时刻系统噪声vk的协方差阵，Gk为k时刻噪声输入矩阵；
S22：估计视轴旋转矩阵Tk，更新转换矩阵Mk
先根据步骤S21的结果，估计
Figure FDA0003782282440000037
的位置向量相对视轴的转动角度βk，转动角度βk为：
Figure FDA0003782282440000038
再估计视轴旋转矩阵Tk，视轴旋转矩阵Tk为：
Figure FDA0003782282440000039
其中，
Figure FDA00037822824400000310
最后更新转换矩阵Mk，Mk为k时刻直角坐标系到视线坐标系的转换矩阵，具体为：
Mk＝TkMk-1；
更新转换矩阵Mk后，视线坐标系的视轴改变，k-1时刻视线坐标系更新为k时刻的视线坐标系；
S23：在k时刻视线坐标系内更新目标预测
Figure FDA00037822824400000311
及其协方差阵
Figure FDA00037822824400000312
表达式为：
Figure FDA00037822824400000313
Figure FDA00037822824400000314
更新后的目标预测为
Figure FDA00037822824400000315
xp、yp是位置预测分量，
Figure FDA00037822824400000316
是速度预测分量；
所述步骤S3的具体过程如下：
S31：基于步骤S23计算得到的
Figure FDA00037822824400000317
和
Figure FDA00037822824400000318
判断是否能够引入多普勒量测，如果
Figure FDA00037822824400000319
说明多普勒信息可用，执行步骤S32，否则执行步骤S33，κ是门限因子，取值不小于5；
Figure FDA0003782282440000041
的表达式如下；
Figure FDA0003782282440000042
其中，
Figure FDA0003782282440000043
是
Figure FDA0003782282440000044
第三行、第三列的元素；
Figure FDA0003782282440000045
S32：构造含多普勒量测的转换量测zk，估计其协方差阵Rk以及量测系数矩阵Hk，其中：
Figure FDA0003782282440000046
Mk(1:3,1:4)是步骤S22中估计出的Mk的第1～3行、第1～4列的所有元素；
基于步骤S23得到的xp、
Figure FDA0003782282440000047
和
Figure FDA0003782282440000048
估计转换量测zk的协方差阵Rk，协方差阵Rk为：
Figure FDA0003782282440000049
其中：
Figure FDA00037822824400000410
Figure FDA00037822824400000411
Figure FDA00037822824400000412
Figure FDA00037822824400000413
Figure FDA00037822824400000414
是
Figure FDA00037822824400000415
的第一行、第一列的元素；
估计量测系数矩阵Hk，量测系数矩阵Hk为：
Figure FDA00037822824400000416
S33：构造无多普勒量测的转换量测zk，估计其协方差阵Rk和量测系数矩阵Hk，其中：
zk＝Mk(1:2,1:2)[rmcosθm rmsinθm]T
Mk(1:2,1:2)是步骤S22估计出的Mk的第1～2行、第1～2列的所有元素；
基于步骤S23得到的xp和
Figure FDA00037822824400000417
估计转换量测zk的协方差阵Rk，协方差阵Rk为：
Figure FDA00037822824400000418
估计量测系数矩阵Hk，量测系数矩阵Hk为：
Figure FDA0003782282440000051
其中：
Figure FDA0003782282440000052
2.根据权利要求1所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：所述步骤S4的具体过程如下：
S41：计算转换量测预测残差
Figure FDA0003782282440000053
Figure FDA0003782282440000054
S42：估计转换量测预测残差
Figure FDA0003782282440000055
的协方差阵Sk
Figure FDA0003782282440000056
S43：估计转换量测滤波增益Kk
Figure FDA0003782282440000057
S44：计算转换量测滤波器状态估计
Figure FDA0003782282440000058
Figure FDA0003782282440000059
S45：计算转换量测滤波估计协方差阵
Figure FDA00037822824400000510
Figure FDA00037822824400000511
3.根据权利要求2所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：在所述步骤S5中，直角坐标系内的状态估计
Figure FDA00037822824400000512
和估计协方差阵
Figure FDA00037822824400000513
分别为：
Figure FDA00037822824400000514
Figure FDA00037822824400000515
4.一种融合多普勒量测的转换量测跟踪系统，其特征在于：根据如权利要求1～3任一项所述的跟踪方法对目标进行跟踪工作，包括：
初始化模块，用于在k＝0时刻，获得雷达量测；在k＝1时刻，再次获取雷达量测，设雷达对目标均匀采样，采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure FDA00037822824400000516
和初始协方差阵
Figure FDA0003782282440000061
以及直角坐标系到视线坐标系的初始转换矩阵M1；
更新模块，用于在k＞1时刻，收到雷达各量测后，预测目标状态及其协方差阵，更新视线坐标系；
模型构造模块，用于在更新后的视线坐标系内，基于更新后的目标预测及其协方差阵，判断是否能够引入多普勒量测，并根据判断结果构造含或不含多普勒量测的转换量测，估计转换量测的协方差阵以及量测系数矩阵；
第一估计模块，用于更新后的目标预测及其协方差阵，结合得到的转换量测、转换量测协方差阵以及量测系数矩阵，估计目标在视线坐标系/直角坐标系的状态估计和协方差阵；
第二估计模块，用于在进入下一个采样时刻时重复步骤S2～S4，递推估计目标在视线坐标系/直角坐标系内的状态估计和协方差阵；
控制处理模块，用于向各模块发出指令，完成相关动作；
所述初始化模块、更新模块、模型构造模块、第一估计模块、第二估计模块均与控制处理模块电连接。
CN202010833608.8A 2020-08-18 2020-08-18 一种融合多普勒量测的转换量测跟踪方法及系统 Active CN111965618B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010833608.8A CN111965618B (zh) 2020-08-18 2020-08-18 一种融合多普勒量测的转换量测跟踪方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010833608.8A CN111965618B (zh) 2020-08-18 2020-08-18 一种融合多普勒量测的转换量测跟踪方法及系统

Publications (2)

Publication Number Publication Date
CN111965618A CN111965618A (zh) 2020-11-20
CN111965618B true CN111965618B (zh) 2022-09-23

Family

ID=73389123

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010833608.8A Active CN111965618B (zh) 2020-08-18 2020-08-18 一种融合多普勒量测的转换量测跟踪方法及系统

Country Status (1)

Country Link
CN (1) CN111965618B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112966720B (zh) * 2021-02-05 2023-05-23 中国电子科技集团公司第三十八研究所 一种基于blue的雷达与红外观测数据融合方法、系统
CN113009468B (zh) * 2021-02-26 2023-05-26 中国电子科技集团公司第三十八研究所 一种视线坐标系内的解耦cmkf跟踪方法及系统
CN116224320B (zh) * 2023-02-17 2023-09-22 昆明理工大学 一种极坐标系下处理多普勒量测的雷达目标跟踪方法
CN116047495B (zh) * 2023-03-31 2023-06-02 昆明理工大学 一种用于三坐标雷达的状态变换融合滤波跟踪方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106646453B (zh) * 2016-11-17 2019-04-05 电子科技大学 一种基于预测值量测转换的多普勒雷达目标跟踪方法
CN106950562B (zh) * 2017-03-30 2020-02-18 电子科技大学 一种基于预测值量测转换的状态融合目标跟踪方法
CN108303692B (zh) * 2018-01-30 2021-02-05 哈尔滨工业大学 一种解多普勒模糊的多目标跟踪方法
CN108896986B (zh) * 2018-04-23 2022-06-03 电子科技大学 一种基于预测值的量测转换序贯滤波机动目标跟踪方法
CN110501696B (zh) * 2019-06-28 2022-05-31 电子科技大学 一种基于多普勒量测自适应处理的雷达目标跟踪方法

Also Published As

Publication number Publication date
CN111965618A (zh) 2020-11-20

Similar Documents

Publication Publication Date Title
CN111965618B (zh) 一种融合多普勒量测的转换量测跟踪方法及系统
CN111985093B (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
CN108226920B (zh) 一种基于预测值处理多普勒量测的机动目标跟踪系统及方法
CN106950562B (zh) 一种基于预测值量测转换的状态融合目标跟踪方法
CN107255795B (zh) 基于ekf/efir混合滤波的室内移动机器人定位方法和装置
CN108896986B (zh) 一种基于预测值的量测转换序贯滤波机动目标跟踪方法
CN110208792B (zh) 同时估计目标状态和轨迹参数的任意直线约束跟踪方法
CN110501696B (zh) 一种基于多普勒量测自适应处理的雷达目标跟踪方法
CN111736145B (zh) 一种基于高斯混合概率假设密度滤波的多机动目标多普勒雷达跟踪方法
CN114609912B (zh) 基于伪线性最大相关熵卡尔曼滤波的仅测角目标追踪方法
CN108319570B (zh) 一种异步多传感器空时偏差联合估计与补偿方法及装置
CN110231620B (zh) 一种噪声相关系统跟踪滤波方法
CN108871365B (zh) 一种航向约束下的状态估计方法及系统
CN111624594A (zh) 一种基于转换量测重构的组网雷达跟踪方法
CN111693984B (zh) 一种改进的ekf-ukf动目标跟踪方法
CN110187337B (zh) 一种基于ls和neu-ecef时空配准的高机动目标跟踪方法及系统
CN115204212A (zh) 一种基于stm-pmbm滤波算法的多目标跟踪方法
CN116047498A (zh) 基于最大相关熵扩展卡尔曼滤波的机动目标跟踪方法
US7928898B2 (en) Method for determining the kinematic state of an object, by evaluating sensor measured values
CN106871905B (zh) 一种非理想条件下高斯滤波替代框架组合导航方法
CN110989341B (zh) 一种约束辅助粒子滤波方法及目标跟踪方法
CN112034445A (zh) 基于毫米波雷达的车辆运动轨迹跟踪方法和系统
CN116929343A (zh) 位姿估计方法、相关设备及存储介质
CN111190173B (zh) 一种基于预测值量测转换的相控阵雷达目标跟踪方法
CN113030945A (zh) 一种基于线性序贯滤波的相控阵雷达目标跟踪方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant