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

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

Info

Publication number
CN111965618A
CN111965618A CN202010833608.8A CN202010833608A CN111965618A CN 111965618 A CN111965618 A CN 111965618A CN 202010833608 A CN202010833608 A CN 202010833608A CN 111965618 A CN111965618 A CN 111965618A
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.)
Granted
Application number
CN202010833608.8A
Other languages
English (en)
Other versions
CN111965618B (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

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radar Systems Or Details Thereof (AREA)

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 BDA0002638886750000011
和初始协方差阵
Figure BDA0002638886750000012
以及直角坐标系到视线坐标系的初始转换矩阵M1；
S2：更新视线坐标系
在k＞1时刻，收到雷达各量测后，预测目标状态及其协方差阵，更新视线坐标系；
S3：构造转换量测模型
在步骤S2更新后的视线坐标系内，基于步骤S2中更新后的目标预测及其协方差阵，判断是否能够引入多普勒量测，并根据判断结果构造含或不含多普勒量测的转换量测，估计转换量测的协方差阵以及量测系数矩阵；
S4：估计目标状态
根据步骤S2中更新后的目标预测及其协方差阵，结合步骤S3中得到的转换量测、转换量测的协方差阵以及量测系数矩阵，估计目标状态，根据需要输出目标在直角坐标系的状态估计和估计协方差阵；
S5：下一时刻目标状态估计
进入下一个采样时刻，重复步骤S2～S4，递推估计目标在视线坐标系内的状态。
更进一步地，在所述步骤S1中，在时刻k，目标状态方程和雷达量测方程分别如下：
Xk＝FkXk-1+Gkvk
zk＝HkXk+wk
其中，
Figure BDA0002638886750000021
为k时刻目标状态，xk、
Figure BDA0002638886750000022
是沿视轴方向的位置、速度分量，yk、
Figure BDA0002638886750000023
是垂直视轴方向的位置、速度分量；
Fk为k时刻状态转移矩阵，状态方程采用近常速运动(NCV：Nearly ConstantVelocity)模型时，其表达式如下：
Figure BDA0002638886750000024
Gk为k时刻噪声输入矩阵，目标近常速运动时的表达式如下：
Figure BDA0002638886750000025
系统噪声为vk＝[vx vy]T，vx、vy分别是沿着视轴方向和垂直视轴方向的零均值高斯过程噪声，其协方差阵为Qk，zk为转换量测值，Hk为量测系数矩阵，wk为量测噪声向量，其协方差阵为Rk。
更进一步地，雷达量测分别为斜距量测rm，方位量测θm和多普勒量测
Figure BDA0002638886750000031
量测噪声是零均值高斯白噪声，各量测的标准差为σr、σθ、
Figure BDA0002638886750000032
斜距量测和多普勒量测的相关系数为ρ。
更进一步地，所述步骤S1的具体过程如下：
S11：在视线坐标系内估计目标初始状态
Figure BDA0002638886750000033
目标初始状态
Figure BDA0002638886750000034
为：
Figure BDA0002638886750000035
S12：估计视线坐标系内的初始协方差阵
Figure BDA0002638886750000036
初始协方差阵
Figure BDA0002638886750000037
为：
Figure BDA0002638886750000038
S13：估计初始转换矩阵M1，初始转换矩阵M1为：
Figure BDA0002638886750000039
其中，
Figure BDA00026388867500000310
更进一步地，所述步骤S2的具体过程如下：
S21：估计目标状态一步预测及其协方差阵，分别为：
Figure BDA00026388867500000311
Figure BDA00026388867500000312
其中，
Figure BDA00026388867500000313
是状态一步预测值，xp,1、yp,1是位置一步预测分量，
Figure BDA00026388867500000314
是速度一步预测分量，Fk是k时刻状态转移矩阵，
Figure BDA00026388867500000315
是k-1时刻状态估计；
Figure BDA00026388867500000316
是状态预测协方差阵；
Figure BDA00026388867500000317
是k-1时刻的状态估计协方差阵；Qk是系统噪声vk的协方差阵，Gk为k时刻噪声输入矩阵；
S22：估计视轴旋转矩阵Tk，更新转换矩阵Mk
先根据步骤S21的结果，估计
Figure BDA0002638886750000041
的位置向量相对视轴的转动角度βk，转动角度βk为：
Figure BDA0002638886750000042
再估计视轴旋转矩阵Tk，视轴旋转矩阵Tk为：
Figure BDA0002638886750000043
其中，
Figure BDA0002638886750000044
最后更新转换矩阵Mk，Mk为直角坐标系到视线坐标系的转换矩阵，具体为：
Mk＝TkMk-1；
更新转换矩阵Mk后，视线坐标系的视轴方向改变，k-1时刻的视线坐标系被更新为k时刻的视线坐标系；
S23：在k时刻视线坐标系内更新目标预测及其协方差阵，分别为：
Figure BDA0002638886750000045
Figure BDA0002638886750000046
更新后的目标预测为
Figure BDA0002638886750000047
xp、yp是位置预测分量，
Figure BDA0002638886750000048
是速度预测分量。
更进一步地，所述步骤S3的具体过程如下：
S31：基于步骤S23计算得到的
Figure BDA0002638886750000049
和
Figure BDA00026388867500000410
判断是否能够引入多普勒量测，如果
Figure BDA00026388867500000411
说明多普勒信息可用，执行步骤S32，否则执行步骤S33；
Figure BDA00026388867500000412
的表达式如下所示；
Figure BDA00026388867500000413
其中
Figure BDA00026388867500000414
是
Figure BDA00026388867500000415
第三行、第三列的元素。
Figure BDA00026388867500000416
S32：构造含多普勒量测的转换量测zk，估计其协方差阵Rk及量测系数矩阵Hk，其中：
Figure BDA0002638886750000051
Mk(1:3,1:4)是步骤S22中估计出的Mk的第1～3行、第1～4列的所有元素；
基于步骤S23得到的xp、
Figure BDA0002638886750000052
和
Figure BDA0002638886750000053
估计转换量测zk的协方差阵Rk，协方差阵Rk为：
Figure BDA0002638886750000054
其中：
Figure BDA0002638886750000055
Figure BDA0002638886750000056
Figure BDA0002638886750000057
Figure BDA0002638886750000058
Figure BDA0002638886750000059
是
Figure BDA00026388867500000510
的第一行、第一列的元素；
估计量测系数矩阵Hk，量测系数矩阵Hk为：
Figure BDA00026388867500000511
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 BDA00026388867500000512
估计转换量测zk的协方差阵Rk，协方差阵Rk为：
Figure BDA00026388867500000513
估计量测系数矩阵Hk，量测系数矩阵Hk为：
Figure BDA00026388867500000514
其中：
Figure BDA00026388867500000515
更进一步地，在所述步骤S31中，κ是门限因子，取值不小于5。
更进一步地，所述步骤S4的具体过程如下：
S41：计算转换量测预测残差
Figure BDA0002638886750000061
Figure BDA0002638886750000062
S42：估计转换量测预测残差
Figure BDA0002638886750000063
的协方差阵Sk
Figure BDA0002638886750000064
S43：估计转换量测滤波增益Kk
Figure BDA0002638886750000065
S44：计算转换量测滤波器状态估计
Figure BDA0002638886750000066
Figure BDA0002638886750000067
S45：计算转换量测滤波估计协方差阵
Figure BDA0002638886750000068
Figure BDA0002638886750000069
更进一步地，在所述步骤S5中，直角坐标系内的状态估计
Figure BDA00026388867500000610
和估计协方差阵
Figure BDA00026388867500000611
为：
Figure BDA00026388867500000612
Figure BDA00026388867500000613
本发明还提供了一种融合多普勒量测的转换量测跟踪系统，包括：
初始化模块，用于在k＝0时刻，获得雷达量测；在k＝1时刻，再次获取雷达量测，设雷达对目标均匀采样，采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure BDA00026388867500000614
和初始协方差阵
Figure BDA00026388867500000615
以及直角坐标系到视线坐标系的初始转换矩阵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 BDA0002638886750000071
为k时刻目标状态，xk、
Figure BDA0002638886750000072
是沿视轴方向的位置、速度分量，yk、
Figure BDA0002638886750000081
是垂直视轴方向的位置、速度分量。
Fk为k时刻状态转移矩阵，状态方程采用近常速运动(NCV：Nearly ConstantVelocity)模型，其表达式如下：
Figure BDA0002638886750000082
T为雷达均匀采样的间隔。Gk为k时刻噪声输入矩阵，目标采取近常速运动时的表达式如下：
Figure BDA0002638886750000083
系统噪声为vk＝[vx vy]T，vx、vy分别是沿着视轴方向和垂直视轴方向的零均值高斯过程噪声，其协方差阵为Qk，zk为转换量测值，Hk为量测系数矩阵，wk为量测噪声向量，其协方差阵为Rk。
雷达量测分别为斜距量测rm，方位量测θm和多普勒量测
Figure BDA0002638886750000084
量测噪声是零均值高斯白噪声，标准差为σr、σθ、
Figure BDA0002638886750000085
和
Figure BDA0002638886750000086
的相关系数为ρ。
基于以下步骤实施对目标的跟踪：
步骤1：滤波器初始化
在k＝0时刻，获得量测rm,0、θm,0、
Figure BDA0002638886750000087
在k＝1时刻，获得量测rm,1、θm,1、
Figure BDA0002638886750000088
假设雷达对目标均匀采样，其采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure BDA0002638886750000089
和初始协方差阵
Figure BDA00026388867500000810
以及直角坐标系到视线坐标系的初始转换矩阵M1；
步骤1.1：在视线坐标系内估计初始状态
Figure BDA00026388867500000811
Figure BDA00026388867500000812
步骤1.2：估计视线坐标系内的初始协方差阵
Figure BDA00026388867500000813
Figure BDA0002638886750000091
步骤1.3：估计初始转换矩阵M1：
Figure BDA0002638886750000092
其中
Figure BDA0002638886750000093
步骤2：在k＞1时刻，收到雷达量测后，预测目标状态及其协方差阵，更新视线坐标系；
步骤2.1：估计目标状态一步预测及其协方差阵：
Figure BDA0002638886750000094
Figure BDA0002638886750000095
其中，
Figure BDA0002638886750000096
是状态一步预测值，xp,1、yp,1是位置一步预测分量，
Figure BDA0002638886750000097
是速度一步预测分量，Fk是k时刻状态转移矩阵，
Figure BDA0002638886750000098
是k-1时刻状态估计；
Figure BDA0002638886750000099
是状态一步预测协方差阵；
Figure BDA00026388867500000910
是k-1时刻的状态估计协方差阵；Qk是系统噪声vk的协方差阵，Gk为k时刻噪声输入矩阵。
步骤2.2：估计视轴旋转矩阵Tk，更新转换矩阵Mk；
首先根据步骤2.1的结果，估计
Figure BDA00026388867500000911
的位置向量相对视轴的转动角度βk：
Figure BDA00026388867500000912
接着估计视轴旋转矩阵Tk：
Figure BDA00026388867500000913
其中，
Figure BDA00026388867500000914
最后更新转换矩阵Mk，Mk是直角坐标系到视线坐标系的转换矩阵；
Mk＝TkMk-1
更新转换矩阵Mk后，视线坐标系的视轴方向改变，k-1时刻的视线坐标系被更新为k时刻的视线坐标系；
步骤2.3：在k时刻视线坐标系内更新目标预测及其协方差阵：
Figure BDA0002638886750000101
Figure BDA0002638886750000102
更新后的目标预测为
Figure BDA00026388867500001019
xp、yp是其位置预测分量，
Figure BDA0002638886750000103
是其速度预测分量。
步骤3：构造转换量测模型。
步骤3.1：基于步骤2.3计算得到的
Figure BDA0002638886750000104
和
Figure BDA0002638886750000105
判断是否能够引入多普勒量测，如果
Figure BDA0002638886750000106
说明多普勒信息可用，执行步骤3.2，否则执行步骤3.3；κ是门限因子，取值不小于5；
Figure BDA0002638886750000107
的表达式如下所示；
Figure BDA0002638886750000108
其中
Figure BDA0002638886750000109
是
Figure BDA00026388867500001010
第三行、第三列的元素，
Figure BDA00026388867500001011
是沿着视轴的速度预测，
Figure BDA00026388867500001012
Figure BDA00026388867500001013
需要说明的是，
Figure BDA00026388867500001014
的分子部分是多普勒转换量测的方差，分母部分是沿视轴的速度预测方差。其含义是：多普勒转换量测相对于沿视轴的速度预测是非线性的，为了保证滤波的鲁棒性，多普勒转换量测方差要数倍于视轴速度预测方差。
步骤3.2：构造含多普勒量测的转换量测zk，估计其协方差阵Rk以及量测系数矩阵Hk；
Figure BDA00026388867500001015
其中，Mk(1:3,1:4)是步骤2.2估计出的Mk的第1～3行、第1～4列的所有元素。
基于步骤2.3得到的xp、
Figure BDA00026388867500001016
和
Figure BDA00026388867500001017
估计转换量测zk的协方差阵Rk：
Figure BDA00026388867500001018
其中：
Figure BDA0002638886750000111
Figure BDA0002638886750000112
Figure BDA0002638886750000113
Figure BDA0002638886750000114
Figure BDA0002638886750000115
是
Figure BDA0002638886750000116
的第一行、第一列的元素；
估计量测系数矩阵Hk：
Figure BDA0002638886750000117
步骤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 BDA0002638886750000118
估计转换量测zk的协方差阵Rk：
Figure BDA0002638886750000119
估计量测系数矩阵Hk：
Figure BDA00026388867500001110
其中：
Figure BDA00026388867500001111
Figure BDA00026388867500001112
步骤4：基于步骤2.3得到的
Figure BDA00026388867500001113
以及步骤3.2或步骤3.3得到的zk、Rk和Hk，估计目标状态。
步骤4.1：计算转换量测预测残差
Figure BDA00026388867500001114
Figure BDA00026388867500001115
步骤4.2：估计转换量测预测残差
Figure BDA00026388867500001116
的协方差阵Sk：
Figure BDA00026388867500001117
其中，Hk是量测系数矩阵，
Figure BDA00026388867500001118
是状态预测协方差阵，Rk是转换量测协方差阵。
步骤4.3：估计转换量测滤波增益Kk：
Figure BDA0002638886750000121
步骤4.4：计算转换量测滤波器状态估计：
Figure BDA0002638886750000122
步骤4.5：计算转换量测滤波估计协方差阵：
Figure BDA0002638886750000123
步骤5：下一个采样时刻，重复步骤2～4，递推估计目标在视线坐标系内的状态，用以下公式输出目标在直角坐标系的状态估计
Figure BDA0002638886750000124
和估计协方差阵
Figure BDA0002638886750000125
Figure BDA0002638886750000126
Figure BDA0002638886750000127
实施例二
在本实施例中，考虑对极坐标下的两种目标跟踪场景进行仿真。
设雷达位于原点，量测噪声为零均值高斯分布，标准差为σr＝10m，σθ＝1°，
Figure BDA0002638886750000128
目标初始位置为(0,50)km，全程做近常速运动，速度为(50,0)m/s。雷达采样间隔为1s，蒙特卡洛仿真次数100次，仿真时长200s。
场景1：目标沿直角坐标系各轴的过程噪声标准差为0.1m/s2，斜距和多普勒量测的相关系数ρ＝0.1。
场景2：目标沿直角坐标系各轴的过程噪声标准差为1m/s2，斜距和多普勒量测的相关系数ρ＝-0.1。
图2和图3分别为本发明所提算法在本实施例的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 BDA0002638886750000131
Figure BDA0002638886750000132
其中，
Figure BDA0002638886750000133
和
Figure BDA0002638886750000134
是第i次仿真时x方向和y方向的状态估计误差，N为仿真次数。RMSE越小，说明算法跟踪精度越高。ANEES可以反映滤波器估计的可信程度，ANEES为1时，表明滤波实际误差和估计误差协方差完全一致，置信度最高。
图2和图3是场景1和场景2下，雷达采用不同方法，对目标的融合跟踪性能。其中图2a、图3a是各方法位置精度对比，图2b、图3b是各方法的ANEES对比。通过与其他方法的对比可见：本发明方法的位置精度和其他方法相当，但估计置信度更高(逼近1)，因此综合性能更优。
表1是三种方法在两个场景中仿真100次的耗时对比。由表1可见，静态融合CMKF需要两个递归滤波器和一个静态融合器同时工作，因此耗时最多；本发明方法和序贯融合CMKF方法的计算量在一个量级。
表1、200次蒙特卡洛仿真耗时对比
Figure BDA0002638886750000135
从仿真结果可见，在两种不同场景中，本发明所提方法的置信度均高于静态融合CMKF和序贯融合CMKF方法，位置精度与其他方法相当，计算量小于静态融合CMKF。综合以上，本发明方法在估计精度、计算复杂度和置信度都有不俗表现，有较好的实用意义。
综上所述，上述实施例的融合多普勒量测的转换量测跟踪方法，通过在视线坐标系内对目标跟踪，克服了多普勒量测和目标状态间的强非线性，能够以较小的计算量估计目标状态，在估计精度和置信度指标上表现良好，具有应用推广价值。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (9)

1.一种融合多普勒量测的转换量测跟踪方法，其特征在于，包括以下步骤：
S1：初始化目标状态
在k＝0时刻，获得雷达量测；在k＝1时刻，再次获取雷达量测，设雷达对目标均匀采样，采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure FDA0002638886740000011
和初始协方差阵
Figure FDA0002638886740000012
以及直角坐标系到视线坐标系的初始转换矩阵M1；
S2：更新视线坐标系
在k＞1时刻，收到雷达各量测后，预测目标状态及其协方差阵，更新视线坐标系；
S3：构造转换量测模型
在步骤S2更新后的视线坐标系内，基于步骤S2中更新后的目标预测及其协方差阵，判断是否能够引入多普勒量测，并根据判断结果构造含或不含多普勒量测的转换量测，估计转换量测的协方差阵以及量测系数矩阵；
S4：估计目标状态
根据步骤S2中更新后的目标预测及其协方差阵，结合步骤S3中得到的转换量测、转换量测协方差阵以及量测系数矩阵，估计目标在视线坐标系/直角坐标系的状态估计和协方差阵；
S5：下一时刻目标状态估计
进入下一个采样时刻，重复步骤S2～S4，递推估计目标在视线坐标系/直角坐标系内的状态估计和协方差阵。
2.根据权利要求1所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：在所述步骤S1中，在时刻k，目标状态方程和雷达量测方程分别如下：
Xk＝FkXk-1+Gkvk
zk＝HkXk+wk
其中，
Figure FDA0002638886740000013
为k时刻目标状态，xk、
Figure FDA0002638886740000014
是沿视轴方向的位置、速度分量，yk、
Figure FDA0002638886740000015
是垂直视轴方向的位置、速度分量；
Fk为k时刻状态转移矩阵，状态方程采用近常速模型，其表达式如下：
Figure FDA0002638886740000016
其中，T是雷达均匀采样的间隔，Gk为k时刻噪声输入矩阵，目标近常速运动的表达式如下：
Figure FDA0002638886740000021
系统噪声为vk＝[vx vy]T，vx、vy分别是沿着视轴方向和垂直视轴方向的零均值高斯过程噪声，其协方差阵为Qk，zk为转换量测值，Hk为量测系数矩阵，wk为量测噪声向量，其协方差阵为Rk。
3.根据权利要求2所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：雷达量测分别为斜距量测rm，方位量测θm和多普勒量测
Figure FDA0002638886740000022
量测噪声是零均值高斯噪声，各量测的标准差为σr、σθ、
Figure FDA00026388867400000211
斜距量测和多普勒量测的相关系数为ρ。
4.根据权利要求1所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：所述步骤S1的具体过程如下：
S11：在视线坐标系内估计目标初始状态
Figure FDA0002638886740000023
目标初始状态
Figure FDA0002638886740000024
为：
Figure FDA0002638886740000025
S12：估计视线坐标系内的初始协方差阵
Figure FDA0002638886740000026
初始协方差阵
Figure FDA0002638886740000027
为：
Figure FDA0002638886740000028
S13：估计初始转换矩阵M1，初始转换矩阵M1为：
Figure FDA0002638886740000029
其中，
Figure FDA00026388867400000210
02×2是2行2列的零矩阵。
5.根据权利要求4所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：所述步骤S2的具体过程如下：
S21：估计目标状态一步预测值及其协方差阵，分别为：
Figure FDA0002638886740000031
Figure FDA0002638886740000032
其中，
Figure FDA0002638886740000033
是k时刻的状态一步预测值，xp,1、yp,1是位置一步预测分量，
Figure FDA0002638886740000034
是速度一步预测分量，Fk是k时刻状态转移矩阵，
Figure FDA0002638886740000035
是k-1时刻状态估计；
Figure FDA0002638886740000036
是k时刻状态一步预测协方差阵；
Figure FDA0002638886740000037
是k-1时刻的状态估计协方差阵；Qk是k时刻系统噪声vk的协方差阵，Gk为k时刻噪声输入矩阵；
S22：估计视轴旋转矩阵Tk，更新转换矩阵Mk
先根据步骤S21的结果，估计
Figure FDA0002638886740000038
的位置向量相对视轴的转动角度βk，转动角度βk为：
Figure FDA0002638886740000039
再估计视轴旋转矩阵Tk，视轴旋转矩阵Tk为：
Figure FDA00026388867400000310
其中，
Figure FDA00026388867400000311
最后更新转换矩阵Mk，Mk为k时刻直角坐标系到视线坐标系的转换矩阵，具体为：
Mk＝TkMk-1；
更新转换矩阵Mk后，视线坐标系的视轴改变，k-1时刻视线坐标系更新为k时刻的视线坐标系；
S23：在k时刻视线坐标系内更新目标预测
Figure FDA00026388867400000312
及其协方差阵
Figure FDA00026388867400000313
表达式为：
Figure FDA00026388867400000314
Figure FDA00026388867400000315
更新后的目标预测为
Figure FDA0002638886740000041
xp、yp是位置预测分量，
Figure FDA0002638886740000042
是速度预测分量。
6.根据权利要求5所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：所述步骤S3的具体过程如下：
S31：基于步骤S23计算得到的
Figure FDA0002638886740000043
和
Figure FDA0002638886740000044
判断是否能够引入多普勒量测，如果
Figure FDA0002638886740000045
说明多普勒信息可用，执行步骤S32，否则执行步骤S33，κ是门限因子，取值不小于5；
Figure FDA0002638886740000046
的表达式如下；
Figure FDA0002638886740000047
其中，
Figure FDA0002638886740000048
是
Figure FDA0002638886740000049
第三行、第三列的元素。
Figure FDA00026388867400000410
S32：构造含多普勒量测的转换量测zk，估计其协方差阵Rk以及量测系数矩阵Hk，其中：
Figure FDA00026388867400000411
Mk(1:3,1:4)是步骤S22中估计出的Mk的第1～3行、第1～4列的所有元素；
基于步骤S23得到的xp、
Figure FDA00026388867400000412
和
Figure FDA00026388867400000413
估计转换量测zk的协方差阵Rk，协方差阵Rk为：
Figure FDA00026388867400000414
其中：
Figure FDA00026388867400000415
Figure FDA00026388867400000416
Figure FDA00026388867400000417
Figure FDA00026388867400000418
Figure FDA00026388867400000419
是
Figure FDA00026388867400000420
的第一行、第一列的元素；
估计量测系数矩阵Hk，量测系数矩阵Hk为：
Figure FDA00026388867400000421
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 FDA0002638886740000051
估计转换量测zk的协方差阵Rk，协方差阵Rk为：
Figure FDA0002638886740000052
估计量测系数矩阵Hk，量测系数矩阵Hk为：
Figure FDA0002638886740000053
其中：
Figure FDA0002638886740000054
7.根据权利要求6所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：所述步骤S4的具体过程如下：
S41：计算转换量测预测残差
Figure FDA0002638886740000055
Figure FDA0002638886740000056
S42：估计转换量测预测残差
Figure FDA0002638886740000057
的协方差阵Sk
Figure FDA0002638886740000058
S43：估计转换量测滤波增益Kk
Figure FDA0002638886740000059
S44：计算转换量测滤波器状态估计
Figure FDA00026388867400000510
Figure FDA00026388867400000511
S45：计算转换量测滤波估计协方差阵
Figure FDA00026388867400000512
Figure FDA00026388867400000513
8.根据权利要求7所述的一种融合多普勒量测的转换量测跟踪方法，其特征在于：在所述步骤S5中，直角坐标系内的状态估计
Figure FDA00026388867400000514
和估计协方差阵
Figure FDA00026388867400000515
分别为：
Figure FDA00026388867400000516
Figure FDA0002638886740000061
9.一种融合多普勒量测的转换量测跟踪系统，其特征在于：根据如权利要求1～8任一项所述的跟踪方法对目标进行跟踪工作，包括：
初始化模块，用于在k＝0时刻，获得雷达量测；在k＝1时刻，再次获取雷达量测，设雷达对目标均匀采样，采样间隔为T，估计k＝1时刻视线坐标内的初始状态估计
Figure FDA0002638886740000062
和初始协方差阵
Figure FDA0002638886740000063
以及直角坐标系到视线坐标系的初始转换矩阵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 true CN111965618A (zh) 2020-11-20
CN111965618B 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)

Cited By (4)

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

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106646453A (zh) * 2016-11-17 2017-05-10 电子科技大学 一种基于预测值量测转换的多普勒雷达目标跟踪方法
CN106950562A (zh) * 2017-03-30 2017-07-14 电子科技大学 一种基于预测值量测转换的状态融合目标跟踪方法
CN108303692A (zh) * 2018-01-30 2018-07-20 哈尔滨工业大学 一种解多普勒模糊的多目标跟踪方法
CN108896986A (zh) * 2018-04-23 2018-11-27 电子科技大学 一种基于预测值的量测转换序贯滤波机动目标跟踪方法
CN110501696A (zh) * 2019-06-28 2019-11-26 电子科技大学 一种基于多普勒量测自适应处理的雷达目标跟踪方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106646453A (zh) * 2016-11-17 2017-05-10 电子科技大学 一种基于预测值量测转换的多普勒雷达目标跟踪方法
CN106950562A (zh) * 2017-03-30 2017-07-14 电子科技大学 一种基于预测值量测转换的状态融合目标跟踪方法
CN108303692A (zh) * 2018-01-30 2018-07-20 哈尔滨工业大学 一种解多普勒模糊的多目标跟踪方法
CN108896986A (zh) * 2018-04-23 2018-11-27 电子科技大学 一种基于预测值的量测转换序贯滤波机动目标跟踪方法
CN110501696A (zh) * 2019-06-28 2019-11-26 电子科技大学 一种基于多普勒量测自适应处理的雷达目标跟踪方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
祝依龙 等: "视线坐标系带多普勒观测的雷达目标跟踪", 《雷达科学与技术》 *

Cited By (8)

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

Also Published As

Publication number Publication date
CN111965618B (zh) 2022-09-23

Similar Documents

Publication Publication Date Title
CN111965618B (zh) 一种融合多普勒量测的转换量测跟踪方法及系统
CN111985093B (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
CN106950562B (zh) 一种基于预测值量测转换的状态融合目标跟踪方法
CN108226920B (zh) 一种基于预测值处理多普勒量测的机动目标跟踪系统及方法
CN107045125B (zh) 一种基于预测值量测转换的交互多模型雷达目标跟踪方法
CN107255795B (zh) 基于ekf/efir混合滤波的室内移动机器人定位方法和装置
CN108896986B (zh) 一种基于预测值的量测转换序贯滤波机动目标跟踪方法
CN110849369B (zh) 机器人跟踪方法、装置、设备及计算机可读存储介质
CN110208792B (zh) 同时估计目标状态和轨迹参数的任意直线约束跟踪方法
CN110501696B (zh) 一种基于多普勒量测自适应处理的雷达目标跟踪方法
CN114609912B (zh) 基于伪线性最大相关熵卡尔曼滤波的仅测角目标追踪方法
CN108319570B (zh) 一种异步多传感器空时偏差联合估计与补偿方法及装置
CN111736145B (zh) 一种基于高斯混合概率假设密度滤波的多机动目标多普勒雷达跟踪方法
CN110231620B (zh) 一种噪声相关系统跟踪滤波方法
CN107688179A (zh) 基于多普勒信息辅助的综合概率数据互联方法
CN108871365B (zh) 一种航向约束下的状态估计方法及系统
CN111624594A (zh) 一种基于转换量测重构的组网雷达跟踪方法
CN115204212A (zh) 一种基于stm-pmbm滤波算法的多目标跟踪方法
CN110187337B (zh) 一种基于ls和neu-ecef时空配准的高机动目标跟踪方法及系统
CN116047498A (zh) 基于最大相关熵扩展卡尔曼滤波的机动目标跟踪方法
CN109239704B (zh) 一种基于序贯滤波交互式多模型的自适应采样方法
Pak Gaussian sum FIR filtering for 2D target tracking
CN112034445A (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