CN111708013A - 一种距离坐标系目标跟踪滤波方法 - Google Patents

一种距离坐标系目标跟踪滤波方法 Download PDF

Info

Publication number
CN111708013A
CN111708013A CN202010620953.3A CN202010620953A CN111708013A CN 111708013 A CN111708013 A CN 111708013A CN 202010620953 A CN202010620953 A CN 202010620953A CN 111708013 A CN111708013 A CN 111708013A
Authority
CN
China
Prior art keywords
distance
state
measurement
coordinate system
expression
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
CN202010620953.3A
Other languages
English (en)
Other versions
CN111708013B (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.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
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 Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202010620953.3A priority Critical patent/CN111708013B/zh
Publication of CN111708013A publication Critical patent/CN111708013A/zh
Application granted granted Critical
Publication of CN111708013B publication Critical patent/CN111708013B/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
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/66—Radar-tracking systems; Analogous systems
    • G01S13/70—Radar-tracking systems; Analogous systems for range tracking only

Landscapes

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

Abstract

本发明涉及一种距离坐标系目标跟踪滤波方法、计算机设备及计算机可读存储介质，该方法包括：在距离‑多普勒子空间对目标匀速运动建模，获得距离坐标系状态方程及对应距离量测的量测方程；从雷达处获取距离量测，若当前跟踪的周期k≤3，则利用基于目标匀速运动模型的初始化方法进行滤波初始化，若当前跟踪的周期k>3则跳过滤波初始化执行下一步骤；利用距离量测进行非线性滤波，根据k‑1周期的状态估计和状态估计协方差，更新k周期的状态估计和状态估计协方差。本发明能够针对不含角度量测和多普勒量测的距离量测序列进行滤波，从而估计出目标运动的距离‑多普勒状态。

Description

一种距离坐标系目标跟踪滤波方法
技术领域
本发明涉及空间目标跟踪技术领域，尤其涉及一种距离坐标系目标跟踪滤波方法、计算机设备及计算机可读存储介质。
背景技术
传统的固定单站雷达需要同时观测目标距离和方位角才能解算目标位置，如果只有距离量测和多普勒量测，或只有距离量测，固定单站雷达无法独立完成对目标的定位和跟踪，需要多站雷达组网探测实现测距交叉定位，从而获得目标状态的可观测性。这一问题类似于无源定位中仅用角度量测的多站时差定位。随着具有抗反辐射导弹能力的被动雷达的发展和应用，无源定位跟踪方法得到了广泛研究。相对于无源定位，无角度量测的目标跟踪是一个比较新的课题，最近几年才受到关注。
目前，无角度量测的目标跟踪研究主要涉及以下两个方面：(1)基于数据关联的集中式融合方法：此方法是处理仅用距离量测和多普勒量测进行目标跟踪最为直接的方法。然而为数众多的鬼影点将会对这种方法形成极大的困扰，由此造成多维分配问题，当目标数目较多时，多维分配的相关算法处理起来就非常复杂，运算量也会很大，实用性和时效性难以控制。(2)基于数据关联的分布式跟踪方法：在描述距离的动态模型时，传感器所提供的多普勒观测能够提供重要信息，因此许多学者提出了分级处理的思想，先对单传感器的距离量测和多普勒量测进行相关，建立某种意义上的局部航迹并剔除杂波，再对局部航迹进行关联去鬼影。然而单传感器距离量测和多普勒量测相关处理中一般采用的都是匀速经验模型或匀加速经验模型，这种模型往往比较粗糙，与真实的距离和多普勒演化规律不符，在进行状态估计相关处理时，难以获得比较理想的性能。当分布式系统中的部分传感器或全部传感器只能提供目标距离量测时，采用交叉定位或分级处理的思想进行跟踪，就需要每个传感器具有根据距离量测序列解算出目标精确距离的能力，因此，需要提供一种能够只依赖距离量测进行目标跟踪的方法。
发明内容
本发明的目的是提供一种针对不含角度量测和多普勒量测的距离量测序列进行滤波，从而估计出目标运动的距离-多普勒状态的目标跟踪滤波方法。
为了实现上述目的，本发明提供了一种距离坐标系目标跟踪滤波方法，该方法包括如下步骤：
S1、在距离-多普勒子空间对目标匀速运动建模，获得距离坐标系状态方程及对应距离量测的量测方程；
S2、从雷达处获取距离量测，若当前跟踪的周期k≤3，则利用基于目标匀速运动模型的初始化方法进行滤波初始化，若当前跟踪的周期k>3则跳过滤波初始化执行步骤S3；其中k为正整数，进行滤波初始化时，获取k＝1、2、3周期的距离量测，利用k＝3周期的状态向量与k＝1、2、3周期的距离真值之间的关系，以距离量测替代距离真值，得到k＝3周期的状态估计，进而利用不敏变换计算k＝3周期的状态估计协方差；
S3、利用距离量测进行非线性滤波，根据k-1周期的状态估计和状态估计协方差，更新k周期的状态估计和状态估计协方差；
S4、判断是否结束跟踪过程，若不结束，则返回执行步骤S2。
优选地，所述步骤S1中，在距离-多普勒子空间对目标匀速运动建模时，在只有距离量测的情况下，量测zk表示为：
Figure BDA0002565131610000021
其中，
Figure BDA0002565131610000022
为目标的距离量测，rk为目标的距离真值，
Figure BDA0002565131610000023
为距离量测误差，距离量测误差为零均值高斯白噪声，方差为
Figure BDA0002565131610000031
获得距离坐标系下状态方程表示为：
xk+1＝f(xk)+vk
对于匀速运动，状态方程表示为：
Figure BDA0002565131610000032
其中，
Figure BDA0002565131610000033
表示状态向量，f为描述状态向量随时间演化规律的非线性函数，
Figure BDA0002565131610000034
表示从k周期的距离-多普勒状态演化而来的、没有过程噪声污染的k+1周期的距离，
Figure BDA0002565131610000035
为多普勒，
Figure BDA0002565131610000036
为转换多普勒的一阶导数，T为雷达采样间隔，vk为过程噪声，q为笛卡尔坐标系中沿x轴和y轴方向过程噪声的标准差，vk的方差表示为：
Figure BDA0002565131610000037
其中，
Figure BDA0002565131610000038
Figure BDA0002565131610000039
Figure BDA00025651316100000310
Figure BDA00025651316100000311
Figure BDA00025651316100000312
Figure BDA00025651316100000313
对应距离量测的量测方程表示为：
Figure BDA0002565131610000041
其中，H为量测矩阵，wk为量测噪声，对应的量测噪声协方差矩阵为Rk。
优选地，所述步骤S2中，进行滤波初始化时，在不考虑随机扰动的情况下，得到状态方程表示为：
Figure BDA0002565131610000042
将包含k、k-1、k-2周期的状态方程联立组成方程组，解方程组，以k、k-1、k-2周期对应的距离真值表示k周期的状态向量的各元素，表达式为：
Figure BDA0002565131610000043
利用距离量测替代距离真值，得到k周期的状态向量表达式为：
Figure BDA0002565131610000044
其中，
Figure BDA0002565131610000045
是由k-2、k-1、k周期对应的距离量测组成的向量，g是表征状态向量和k-2、k-1、k周期对应的距离量测组成的向量之间非线性关系的向量值函数；带入k＝3，计算相应的状态向量作为k＝3周期的状态估计；
采用不敏变换计算k＝3周期的状态估计协方差。
优选地，所述步骤S2中，采用不敏变换计算k＝3周期的状态估计协方差时，包括如下步骤：
首先计算向量rk的2nx+1个采样点
Figure BDA0002565131610000051
及其相应的权值Wi，表达式为：
Figure BDA0002565131610000052
其中
Figure BDA0002565131610000053
nx是向量rk的维数，λ是满足nx+λ≠0的标量参数，
Figure BDA0002565131610000054
是矩阵
Figure BDA0002565131610000055
均方根的第i行或第i列；
然后计算各采样点的映射值和相应的转移状态，映射值表达式为：
Figure BDA0002565131610000056
转移状态表达式为：
Figure BDA0002565131610000057
最后计算状态估计协方差，表达式为：
Figure BDA0002565131610000058
优选地，所述步骤S3中，根据距离量测进行非线性滤波时，采用转换量测卡尔曼滤波方法、无迹卡尔曼滤波方法、扩展卡尔曼滤波方法或粒子滤波方法中的一种。
优选地，所述步骤S3中，根据距离量测进行非线性滤波时，采用无迹卡尔曼滤波方法，从k＝4周期开始滤波，包括如下步骤：
S3-1、通过不敏变换计算2nx+1个采样点
Figure BDA0002565131610000059
及相应的权重Wi，表达式为：
Figure BDA00025651316100000510
其中，nx是状态向量
Figure BDA00025651316100000511
的维度，λ是满足nx+λ≠0的标量参数，
Figure BDA0002565131610000061
是矩阵(nx+λ)Pk-1|k-1均方根的第i行或第i列；
S3-2、计算状态一步预测
Figure BDA0002565131610000062
表达式为：
Figure BDA0002565131610000063
S3-3、计算一步预测协方差Pk|k-1，表达式为：
Figure BDA0002565131610000064
S3-4、计算滤波增益Kk，表达式为：
Figure BDA0002565131610000065
Figure BDA0002565131610000066
Kk＝Pxz(Pzz)-1
其中
Figure BDA0002565131610000067
为采样点对应的量测预测，
Figure BDA0002565131610000068
为量测预测，Pzz为量测预测协方差矩阵，Pxz为状态和量测之间的互协方差矩阵；
S3-5、更新状态估计
Figure BDA0002565131610000069
表达式为：
Figure BDA00025651316100000610
S3-6、更新状态估计协方差Pk|k，表达式为：
Pk|k＝Pk|k-1-KkPzz(Kk)′。
本发明还提供了一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现上述任一项所述的距离坐标系目标跟踪滤波方法的步骤。
本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项所述的距离坐标系目标跟踪滤波方法的步骤。
本发明的上述技术方案具有如下优点：本发明提供了一种距离坐标系目标跟踪滤波方法、计算机设备及计算机可读存储介质，本发明对目标匀速运动在距离-多普勒子空间建模，针对不含角度量测和多普勒量测的距离量测序列进行滤波，从而估计出目标运动的距离-多普勒状态。本发明提供了一种只依赖距离量测进行目标跟踪的方法，并且在初始化时利用目标匀速运动模型进行滤波初始化，相比常用的两点差分法，由于引入了真实目标运动模型信息，可以获得更好的滤波初始化精度。
附图说明
图1示出了本发明实施例提供的一种距离坐标系目标跟踪滤波方法步骤示意图；
图2示出了采用两点差分初始化的距离滤波方法和本发明实施例提供的跟踪滤波方法分别得到的距离均方根误差对比结果；
图3示出了采用两点差分初始化的距离滤波方法和本发明实施例提供的跟踪滤波方法分别得到的多普勒均方根误差对比结果；
图4示出了采用两点差分初始化的距离滤波方法和本发明实施例提供的跟踪滤波方法分别得到的距离-多普勒状态向量第三个元素均方根误差对比结果；
图5示出了采用两点差分初始化的距离滤波方法和本发明实施例提供的跟踪滤波方法分别得到的平均归一化误差平方对比结果。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明实施例提供的一种距离坐标系目标跟踪滤波方法，包括如下步骤：
S1、在距离-多普勒子空间对目标匀速运动建模，获得距离坐标系状态方程及对应距离量测的量测方程。
S2、从雷达处获取距离量测，若当前跟踪的周期k≤3，则利用基于目标匀速运动模型的初始化方法进行滤波初始化，若当前跟踪的周期k>3则跳过滤波初始化执行步骤S3。其中k为雷达的扫描周期数(也即跟踪周期数)，是正整数。
进行滤波初始化时，获取k＝1、2、3周期的距离量测，利用k＝1、2、3周期的距离量测得到k＝3周期的状态估计和状态估计协方差，以k＝3周期对应的状态估计和状态估计协方差为初始状态和初始协方差，确定初始状态和初始协方差后，即完成滤波初始化。需要说明的是，在k≤3的情况下，即在未完成滤波初始化的情况下，不执行步骤S3。具体地，进行滤波初始化时，利用k＝3周期的状态向量与k＝1、2、3周期的距离真值之间的关系，以距离量测替代距离真值，得到k＝3周期的状态估计，进而利用不敏变换(UnscentedTransformation，UT)计算k＝3周期的状态估计协方差。
S3、利用距离量测进行非线性滤波，根据k-1周期的状态估计和状态估计协方差，更新k周期的状态估计和状态估计协方差。
S4、判断是否结束跟踪过程，若不结束，则返回执行步骤S2。
优选地，步骤S1中，在距离-多普勒子空间对目标匀速运动建模时，在只有距离量测的情况下，量测zk表示为：
Figure BDA0002565131610000081
其中，
Figure BDA0002565131610000082
为(k周期)目标的距离量测，对于只有距离量测的情况，量测zk等于距离量测
Figure BDA0002565131610000083
rk为(k周期)目标的距离真值，
Figure BDA0002565131610000084
为(k周期)距离量测误差，距离量测误差为零均值高斯白噪声，方差为
Figure BDA0002565131610000085
距离坐标系状态方程可表示为：
xk+1＝f(xk)+vk
其中，xk为(k周期)距离-多普勒子空间的状态向量(简称状态向量，或距离-多普勒状态向量)，f为描述状态向量随时间演化规律的非线性函数，vk为过程噪声。
对于匀速运动，获得对应的距离坐标系的状态向量和状态方程分别表示为：
Figure BDA0002565131610000091
Figure BDA0002565131610000092
其中，
Figure BDA0002565131610000093
表示从k周期的距离-多普勒状态演化而来的、没有过程噪声污染的k+1周期的距离，
Figure BDA0002565131610000094
为多普勒(真值)，
Figure BDA0002565131610000095
为转换多普勒的一阶导数(真值)，T为雷达采样间隔，q为笛卡尔坐标系中沿x轴和y轴方向过程噪声的标准差，vk为过程噪声，是距离-多普勒子空间的零均值高斯噪声，过程噪声vk的方差为：
Figure BDA0002565131610000096
其中，
Figure BDA0002565131610000097
Figure BDA0002565131610000098
Figure BDA0002565131610000099
Figure BDA00025651316100000910
Figure BDA00025651316100000911
Figure BDA00025651316100000912
由过程噪声vk的方差的各矩阵元素可以看出，距离-多普勒子空间的过程噪声只与距离、多普勒和转换多普勒的一阶导数有关，而与目标在笛卡尔空间的状态无关。考虑到在实际应用中，由于真值rk、
Figure BDA00025651316100000916
和
Figure BDA00025651316100000913
未知，可以用相应的估计值
Figure BDA00025651316100000914
和
Figure BDA00025651316100000915
替代。
距离量测是距离-多普勒状态向量的线性函数，根据匀速运动对应的状态向量形式，对应距离量测的量测方程表示为：
Figure BDA0002565131610000101
其中，H为量测矩阵，wk为量测噪声，对应的量测噪声协方差矩阵为Rk。
本发明根据仅有的距离量测，利用上面的距离坐标系状态方程和相应的量测方程，针对匀速运动，能够估计出对应的距离-多普勒状态，实现目标跟踪。
优选地，步骤S2中，利用基于目标匀速运动模型的初始化方法进行滤波初始化时，在不考虑随机扰动的情况下，可得到状态方程表示为：
Figure BDA0002565131610000102
将包含k、k-1、k-2周期的状态方程联立组成方程组，解方程组，以k、k-1、k-2周期对应的距离真值表示k周期的状态向量的各元素，表达式为：
Figure BDA0002565131610000103
利用k、k-1、k-2周期对应的距离量测替代距离真值，得到k周期的状态向量表达式为：
Figure BDA0002565131610000104
其中，
Figure BDA0002565131610000111
是由k-2、k-1、k周期对应的距离量测组成的向量，g是表征状态向量和k-2、k-1、k周期对应的距离量测组成的向量之间非线性关系的向量值函数。对于k＝3周期，相应状态向量可由k＝1、k＝2、k＝3周期对应的距离量测表示，即初始状态向量可由过去连续三个扫描间隔的距离量测计算。带入k＝3，计算相应的状态向量作为k＝3周期的状态估计，即确定初始状态。
通过上述k周期的状态向量表达式计算出来的初始状态向量与真值之间的误差项中，分母部分含有距离量测误差
Figure BDA00025651316100001110
的表述，在计算协方差求数学期望的积分过程中，会出现不可积的情况。为了便于处理，本发明采用不敏变换计算k＝3周期的状态估计协方差。
进一步地，步骤S2中，采用不敏变换计算k＝3周期的状态估计协方差时，包括如下步骤：
首先计算向量rk的2nx+1个采样点
Figure BDA0002565131610000112
及其相应的权值Wi，表达式为：
Figure BDA0002565131610000113
其中
Figure BDA0002565131610000114
nx是向量rk的维数，λ是满足nx+λ≠0的标量参数，
Figure BDA0002565131610000115
是矩阵
Figure BDA0002565131610000116
均方根的第i行或第i列；
然后计算各采样点的映射值和相应的转移状态，其中映射值表达式为：
Figure BDA0002565131610000117
转移状态表达式为：
Figure BDA0002565131610000118
最后计算状态估计协方差，表达式为：
Figure BDA0002565131610000119
带入k＝3计算k＝3周期的状态估计协方差，即可得到初始协方差。
由于量测与目标状态之间是非线性关系，因此在滤波过程中需要采用非线性滤波方法，常用的非线性滤波方法包括转换量测卡尔曼滤波方法、无迹卡尔曼滤波方法、扩展卡尔曼滤波方法以及粒子滤波方法等，可采用其中一种实现非线性滤波。
进一步地，步骤S3中，根据距离量测进行非线性滤波时，采用无迹卡尔曼滤波方法，从k＝4周期开始滤波，包括如下步骤：
S3-1、通过不敏变换计算(在k-1周期的状态估计
Figure BDA0002565131610000121
附近选取的)2nx+1个采样点
Figure BDA0002565131610000122
及相应的权重Wi，表达式为：
Figure BDA0002565131610000123
其中，nx是状态向量
Figure BDA0002565131610000124
的维度，λ是满足nx+λ≠0的标量参数，
Figure BDA0002565131610000125
是矩阵(nx+λ)Pk-1|k-1均方根的第i行或第i列；
S3-2、计算状态一步预测
Figure BDA0002565131610000126
表达式为：
Figure BDA0002565131610000127
S3-3、计算一步预测协方差Pk|k-1，表达式为：
Figure BDA0002565131610000128
S3-4、计算滤波增益Kk，表达式为：
Figure BDA0002565131610000129
Figure BDA00025651316100001210
Kk＝Pxz(Pzz)-1
其中，
Figure BDA00025651316100001211
为采样点对应的量测预测，
Figure BDA00025651316100001212
为量测预测，Pzz为量测预测协方差矩阵，Pxz为状态和量测之间的互协方差矩阵；
S3-5、更新(k周期的)状态估计
Figure BDA00025651316100001213
表达式为：
Figure BDA00025651316100001214
S3-6、更新(k周期的)状态估计协方差Pk|k，表达式为：
Pk|k＝Pk|k-1-KkPzz(Kk)′。
在完成滤波初始化后，从k＝4周期开始迭代，根据步骤S2得到的、k＝3周期的状态估计(初始状态)和状态估计协方差(初始协方差)更新k＝4周期的状态估计和状态估计协方差，下一次计算根据距离量测基k＝4周期的状态估计和状态估计协方差更新k＝5周期的状态估计和状态估计协方差，以此类推进行非线性滤波。
当目标距离随时间线性变化时，利用两点差分法初始化滤波器比较精确。但实际情况并不如此，在初始化阶次比较低的情况下，两点差分法也只是近似满足精度要求。在本发明中，距离量测被用来初始化转换多普勒的高阶导数，两点差分法会带来比较大的近似误差。因此，本发明基于目标匀速运动模型的状态方程提出了一种新的初始化方法，首先利用确定性系统的状态演化方程，推导当前周期的目标状态和过去几个连续扫描间隔的距离真值之间的函数关系，然后用过去几个连续扫描间隔的距离量测替代距离真值，表示当前周期的目标状态；同时根据它们之间的函数关系，利用UT变换计算初始协方差。因为状态模型(即目标匀速运动模型)是精确的，新的初始化方法也相对精确，特别是在距离随时间变化非线性的高阶场景下。
为了验证本发明提出的仅用距离观测进行状态估计方法的有效性和基于目标匀速运动模型的初始化方法的优越性，本发明还进行了相应的数值仿真及性能比较。如图2至图5所示，基于1000次蒙特卡洛数值仿真，采用均方根误差(Root Mean Squared Error，RMSE)评价估计性能，并利用后验克拉美罗下限(Posterior Cramer-Rao Lower Bound，PCRLB)作为可能达到最优性能的参考，同时，采用平均归一化误差平方(AverageNormalized Error Squared，ANES)评价估计的一致性。对比的方法为采用两点差分初始化的距离滤波方法与本发明提供的距离坐标系目标跟踪滤波方法。
仿真场景设定雷达位于坐标原点，雷达采样间隔T＝5s，即每隔5s雷达传回目标距离量测。目标做匀速运动，在笛卡尔坐标系下的过程噪声设定为零均值高斯白噪声，其标准差设定为典型值q＝0.01m/s2。考虑两种典型场景：距离雷达较远的低速运动，初始位置为(50km,50km)，初始速度为10m/s，航向-45度，量测噪声标准差σr＝300m；距离雷达较近的高速运动，初始位置为(5km,5km)，初始速度为500m/s，航向-45度，量测噪声标准差σr＝50m。雷达扫描次数(即k的最大值)设为100。
图2-图5分别示出了初始状态为(5km,5km),(500m/s,-45度),σr＝50m场景下，采用两点差分初始化的距离滤波方法(简称两点差分法)和本发明实施例提供的跟踪滤波方法(简称本发明所提方法)分别得到的距离RMSE、多普勒RMSE、距离-多普勒状态向量第三个元素RMSE，以及ANES对比结果，可以看到本发明所提方法随着迭代次数的增加很快收敛，并且接近于PCRLB，ANES也落在98％置信区间之内。这说明对于匀速运动在距离-多普勒子空间建立的模型是精确的，本发明在只有距离量测的情况下采用该模型可有效估计距离-多普勒状态。
同时，由图2至图5可以明显看到采用本发明所提方法进行初始化，距离-多普勒状态向量各分量的RMSE相比两点差分法进行初始化收敛更快，并且能始终保持良好的一致性，而两点差分法滤波的一致性明显恶化了。这是因为在距离雷达较近量测噪声较小的情况下，距离-多普勒状态向量各分量伴随着高速匀速运动的非线性变化效应更突出，对此场景下的运动目标跟踪滤波初始化，两点差分法相比本发明所提方法初始化处理更为粗糙，如图3和图4所示，两点差分法对状态向量后两个分量的初始化误差明显比本发明所提方法误差更大，而这种误差一直伴随着整个滤波过程。
特别地，在本发明一些优选的实施方式中，还提供了一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现上述任一实施方式中所述的距离坐标系目标跟踪滤波方法的步骤。
在本发明另一些优选的实施方式中，还提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述任一实施方式中所述的距离坐标系目标跟踪滤波方法的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程，在此不再重复说明。
综上所述，本发明对目标匀速运动在距离-多普勒子空间建模，推导了准确的距离坐标系状态方程，提出了一种新的距离坐标系目标跟踪滤波方法，能够针对不含角度量测和多普勒量测的距离量测序列进行滤波，从而估计出目标运动的距离-多普勒状态。本发明还提出了基于状态模型的初始化方法，根据状态转移矩阵和距离量测计算滤波器初值(即初始状态)，利用UT方法计算初始协方差。相比常用的两点差分法，由于引入了真实目标运动模型信息，可以获得更好的滤波初始化精度。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (8)

1.一种距离坐标系目标跟踪滤波方法，其特征在于，包括如下步骤：
S1、在距离-多普勒子空间对目标匀速运动建模，获得距离坐标系状态方程及对应距离量测的量测方程；
S2、从雷达处获取距离量测，若当前跟踪的周期k≤3，则利用基于目标匀速运动模型的初始化方法进行滤波初始化，若当前跟踪的周期k>3则跳过滤波初始化执行步骤S3；其中k为正整数，进行滤波初始化时，获取k＝1、2、3周期的距离量测，利用k＝3周期的状态向量与k＝1、2、3周期的距离真值之间的关系，以距离量测替代距离真值，得到k＝3周期的状态估计，进而利用不敏变换计算k＝3周期的状态估计协方差；
S3、利用距离量测进行非线性滤波，根据k-1周期的状态估计和状态估计协方差，更新k周期的状态估计和状态估计协方差；
S4、判断是否结束跟踪过程，若不结束，则返回执行步骤S2。
2.根据权利要求1所述的距离坐标系目标跟踪滤波方法，其特征在于，
所述步骤S1中，在距离-多普勒子空间对目标匀速运动建模时，在只有距离量测的情况下，量测zk表示为：
Figure FDA0002565131600000011
其中，
Figure FDA0002565131600000012
为目标的距离量测，rk为目标的距离真值，
Figure FDA0002565131600000013
为距离量测误差，距离量测误差为零均值高斯白噪声，方差为
Figure FDA0002565131600000014
获得距离坐标系下状态方程表示为：
xk+1＝f(xk)+vk
对于匀速运动，状态方程表示为：
Figure FDA0002565131600000021
其中，
Figure FDA0002565131600000022
表示状态向量，f为描述状态向量随时间演化规律的非线性函数，
Figure FDA0002565131600000023
表示从k周期的距离-多普勒状态演化而来的、没有过程噪声污染的k+1周期的距离，
Figure FDA0002565131600000024
为多普勒，
Figure FDA0002565131600000025
为转换多普勒的一阶导数，T为雷达采样间隔，vk为过程噪声，q为笛卡尔坐标系中沿x轴和y轴方向过程噪声的标准差，vk的方差表示为：
Figure FDA0002565131600000026
其中，
Figure FDA0002565131600000027
Figure FDA0002565131600000028
Figure FDA0002565131600000029
Figure FDA00025651316000000210
Figure FDA00025651316000000211
Figure FDA00025651316000000212
对应距离量测的量测方程表示为：
Figure FDA00025651316000000213
其中，H为量测矩阵，wk为量测噪声，对应的量测噪声协方差矩阵为Rk。
3.根据权利要求2所述的距离坐标系目标跟踪滤波方法，其特征在于，
所述步骤S2中，进行滤波初始化时，在不考虑随机扰动的情况下，得到状态方程表示为：
Figure FDA0002565131600000031
将包含k、k-1、k-2周期的状态方程联立组成方程组，解方程组，以k、k-1、k-2周期对应的距离真值表示k周期的状态向量的各元素，表达式为：
Figure FDA0002565131600000032
利用距离量测替代距离真值，得到k周期的状态向量表达式为：
Figure FDA0002565131600000033
其中，
Figure FDA0002565131600000034
是由k-2、k-1、k周期对应的距离量测组成的向量，g是表征状态向量和k-2、k-1、k周期对应的距离量测组成的向量之间非线性关系的向量值函数；带入k＝3，计算相应的状态向量作为k＝3周期的状态估计；
采用不敏变换计算k＝3周期的状态估计协方差。
4.根据权利要求3所述的距离坐标系目标跟踪滤波方法，其特征在于，
所述步骤S2中，采用不敏变换计算k＝3周期的状态估计协方差时，包括如下步骤：
首先计算向量rk的2nx+1个采样点
Figure FDA0002565131600000041
及其相应的权值Wi，表达式为：
Figure FDA0002565131600000042
其中
Figure FDA0002565131600000043
nx是向量rk的维数，λ是满足nx+λ≠0的标量参数，
Figure FDA0002565131600000044
是矩阵
Figure FDA0002565131600000045
均方根的第i行或第i列；
然后计算各采样点的映射值和相应的转移状态，映射值表达式为：
Figure FDA0002565131600000046
转移状态表达式为：
Figure FDA0002565131600000047
最后计算状态估计协方差，表达式为：
Figure FDA0002565131600000048
5.根据权利要求4所述的距离坐标系目标跟踪滤波方法，其特征在于：所述步骤S3中，根据距离量测进行非线性滤波时，采用转换量测卡尔曼滤波方法、无迹卡尔曼滤波方法、扩展卡尔曼滤波方法或粒子滤波方法中的一种。
6.根据权利要求5所述的距离坐标系目标跟踪滤波方法，其特征在于，
所述步骤S3中，根据距离量测进行非线性滤波时，采用无迹卡尔曼滤波方法，从k＝4周期开始滤波，包括如下步骤：
S3-1、通过不敏变换计算2nx+1个采样点
Figure FDA0002565131600000049
及相应的权重Wi，表达式为：
Figure FDA0002565131600000051
其中，nx是状态向量
Figure FDA0002565131600000052
的维度，λ是满足nx+λ≠0的标量参数，
Figure FDA0002565131600000053
是矩阵(nx+λ)Pk-1|k-1均方根的第i行或第i列；
S3-2、计算状态一步预测
Figure FDA0002565131600000054
表达式为：
Figure FDA0002565131600000055
S3-3、计算一步预测协方差Pk|k-1，表达式为：
Figure FDA0002565131600000056
S3-4、计算滤波增益Kk，表达式为：
Figure FDA00025651316000000512
Figure FDA0002565131600000057
Kk＝Pxz(Pzz)-1
其中
Figure FDA0002565131600000058
为采样点对应的量测预测，
Figure FDA0002565131600000059
为量测预测，Pzz为量测预测协方差矩阵，Pxz为状态和量测之间的互协方差矩阵；
S3-5、更新状态估计
Figure FDA00025651316000000510
表达式为：
Figure FDA00025651316000000511
S3-6、更新状态估计协方差Pk|k，表达式为：
Pk|k＝Pk|k-1-KkPzz(Kk)′。
7.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至6中任一项所述的距离坐标系目标跟踪滤波方法的步骤。
8.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至6中任一项所述的距离坐标系目标跟踪滤波方法的步骤。
CN202010620953.3A 2020-07-01 2020-07-01 一种距离坐标系目标跟踪滤波方法 Active CN111708013B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010620953.3A CN111708013B (zh) 2020-07-01 2020-07-01 一种距离坐标系目标跟踪滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010620953.3A CN111708013B (zh) 2020-07-01 2020-07-01 一种距离坐标系目标跟踪滤波方法

Publications (2)

Publication Number Publication Date
CN111708013A true CN111708013A (zh) 2020-09-25
CN111708013B CN111708013B (zh) 2022-06-07

Family

ID=72544930

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010620953.3A Active CN111708013B (zh) 2020-07-01 2020-07-01 一种距离坐标系目标跟踪滤波方法

Country Status (1)

Country Link
CN (1) CN111708013B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113219406A (zh) * 2021-04-29 2021-08-06 电子科技大学 一种基于扩展卡尔曼滤波的直接跟踪方法
CN114942429A (zh) * 2022-05-31 2022-08-26 哈尔滨工业大学 仅用距离-多普勒观测的双站雷达协同跟踪方法及装置
CN115032626A (zh) * 2022-06-07 2022-09-09 哈尔滨工业大学 距离平方域的目标跟踪方法、装置、电子设备及存储介质
CN115451962A (zh) * 2022-08-09 2022-12-09 中国人民解放军63629部队 一种基于五变量卡诺图的目标跟踪策略规划方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2009103242A (ru) * 2009-02-02 2010-08-10 Общество с ограниченной ответственностью "ОКБ Траверз" (RU) Способ сопровождения пилотируемой воздушной цели
CN102508238A (zh) * 2011-10-14 2012-06-20 北京理工大学 一种基于坐标旋转变换的雷达跟踪方法
CN103048655A (zh) * 2013-01-11 2013-04-17 中国人民解放军空军预警学院 天波超视距雷达频域超分辨微多径测高方法
CN105954742A (zh) * 2016-05-19 2016-09-21 哈尔滨工业大学 一种球坐标系下带多普勒观测的雷达目标跟踪方法
US9465108B1 (en) * 2014-12-03 2016-10-11 The United States Of America As Represented By The Secretary Of The Navy System and method for target doppler estimation and range bias compensation using high duty cycle linear frequency modulated signals

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2009103242A (ru) * 2009-02-02 2010-08-10 Общество с ограниченной ответственностью "ОКБ Траверз" (RU) Способ сопровождения пилотируемой воздушной цели
CN102508238A (zh) * 2011-10-14 2012-06-20 北京理工大学 一种基于坐标旋转变换的雷达跟踪方法
CN103048655A (zh) * 2013-01-11 2013-04-17 中国人民解放军空军预警学院 天波超视距雷达频域超分辨微多径测高方法
US9465108B1 (en) * 2014-12-03 2016-10-11 The United States Of America As Represented By The Secretary Of The Navy System and method for target doppler estimation and range bias compensation using high duty cycle linear frequency modulated signals
CN105954742A (zh) * 2016-05-19 2016-09-21 哈尔滨工业大学 一种球坐标系下带多普勒观测的雷达目标跟踪方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
CHENWANG 等: "Real-time FPGA-based Kalman filter for constant and non-constant velocity periodic error correction", 《PRECISION ENGINEERING》 *
GONGJIAN ZHOU 等: "Statically Fused Converted Position and Doppler Measurement Kalman Filters", 《IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS》 *
王瑞昕: "基于毫米波雷达的多目标检测与跟踪技术研究", 《中国优秀博硕士学位论文全文数据库(硕士)工程科技Ⅱ辑》 *
荣里 等: "基于极坐标多普勒伪状态的最佳线性无偏估计算法", 《火力与指挥控制》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113219406A (zh) * 2021-04-29 2021-08-06 电子科技大学 一种基于扩展卡尔曼滤波的直接跟踪方法
CN113219406B (zh) * 2021-04-29 2022-08-05 电子科技大学 一种基于扩展卡尔曼滤波的直接跟踪方法
CN114942429A (zh) * 2022-05-31 2022-08-26 哈尔滨工业大学 仅用距离-多普勒观测的双站雷达协同跟踪方法及装置
CN115032626A (zh) * 2022-06-07 2022-09-09 哈尔滨工业大学 距离平方域的目标跟踪方法、装置、电子设备及存储介质
CN115032626B (zh) * 2022-06-07 2024-07-12 哈尔滨工业大学 距离平方域的目标跟踪方法、装置、电子设备及存储介质
CN115451962A (zh) * 2022-08-09 2022-12-09 中国人民解放军63629部队 一种基于五变量卡诺图的目标跟踪策略规划方法
CN115451962B (zh) * 2022-08-09 2024-04-30 中国人民解放军63629部队 一种基于五变量卡诺图的目标跟踪策略规划方法

Also Published As

Publication number Publication date
CN111708013B (zh) 2022-06-07

Similar Documents

Publication Publication Date Title
CN111708013B (zh) 一种距离坐标系目标跟踪滤波方法
CN111722214A (zh) 雷达多目标跟踪phd实现方法
CN111722213B (zh) 一种机动目标运动参数的纯距离提取方法
CN111711432B (zh) 一种基于ukf和pf混合滤波的目标跟踪算法
CN111736144B (zh) 一种仅用距离观测的机动转弯目标状态估计方法
Cao et al. Extended object tracking using automotive radar
Jahan et al. Implementation of underwater target tracking techniques for Gaussian and non-Gaussian environments
Campbell et al. An algorithm for large-scale multitarget tracking and parameter estimation
Sun et al. Modeling and tracking of maneuvering extended object with random hypersurface
CN116609776A (zh) 复杂环境下的基于人工势场法的星凸形扩展目标跟踪方法
CN116224320B (zh) 一种极坐标系下处理多普勒量测的雷达目标跟踪方法
CN116047495B (zh) 一种用于三坐标雷达的状态变换融合滤波跟踪方法
Zhou et al. State estimation with destination constraints
Hou et al. Bearing‐only underwater uncooperative target tracking for non‐Gaussian environment using fast particle filter
Bunch et al. Dynamical models for tracking with the variable rate particle filter
Li et al. Motion modeling and state estimation in range-squared coordinate
Xiong et al. The linear fitting Kalman filter for nonlinear tracking
CN115327503A (zh) 基于高斯粒子滤波的固定单站无源定位方法和相关装置
Voronina et al. Algorithm for constructing trajectories of maneuvering object based on bearing-only information using the Basis Pursuit method
Chen et al. Improved Unscented Kalman Filtering Algorithm Applied to On-vehicle Tracking System
Ma et al. Variational Bayesian cubature Kalman filter for bearing-only tracking in glint noise environment
Sönmez et al. Analysis of performance criteria for optimization based bearing only target tracking algorithms
Sun et al. A new maneuvering target tracking method using adaptive cubature Kalman filter
Qian et al. Improvement of bearings only target tracking using smoothing
Moawad et al. Study a bearing-only moving ground target tracking problem using single seismic sensor

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