CN103955892A - 一种目标跟踪方法及扩展截断无迹卡尔曼滤波方法、装置 - Google Patents

一种目标跟踪方法及扩展截断无迹卡尔曼滤波方法、装置 Download PDF

Info

Publication number
CN103955892A
CN103955892A CN201410134331.4A CN201410134331A CN103955892A CN 103955892 A CN103955892 A CN 103955892A CN 201410134331 A CN201410134331 A CN 201410134331A CN 103955892 A CN103955892 A CN 103955892A
Authority
CN
China
Prior art keywords
probability density
density function
sigma
constantly
posterior probability
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
CN201410134331.4A
Other languages
English (en)
Other versions
CN103955892B (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.)
Shenzhen University
Original Assignee
Shenzhen University
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 Shenzhen University filed Critical Shenzhen University
Priority to CN201410134331.4A priority Critical patent/CN103955892B/zh
Publication of CN103955892A publication Critical patent/CN103955892A/zh
Application granted granted Critical
Publication of CN103955892B publication Critical patent/CN103955892B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Image Analysis (AREA)

Abstract

本发明公开了一种目标跟踪方法、系统及扩展截断无迹卡尔曼滤波方法、装置,该扩展截断无迹卡尔曼滤波方法包括根据无迹变换获取原始先验概率密度函数;根据原始先验概率密度函数获取第一后验概率密度函数;根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数;根据修正先验概率密度函数获取第二后验概率密度函数;根据原始先验概率密度函数以及第二后验概率密度函数获取联合后验概率密度函数。本发明能够解决观测函数不具有唯一反函数的问题,有效减少目标状态先验分布方差,自适应地根据观测信息的精度进行状态更新,有效提高滤波精度且实用性较高。

Description

一种目标跟踪方法及扩展截断无迹卡尔曼滤波方法、装置
技术领域
本发明涉及非线性滤波领域,特别是涉及一种目标跟踪方法、系统及扩展截断无迹卡尔曼滤波方法、装置。
背景技术
在导航和制导系统、跟踪快速变化的无线信道的信道状态信息、跟踪飞机的实时位置等科学领域中经常需要使用滤波技术以实现对目标的实时跟踪。现有技术中常采用如下几种滤波方法:第一种为由卡尔曼于1960年提出的卡尔曼滤波(KF),该方法可以处理具有高斯分布噪声的线性系统以得到系统状态的最小均方差估计RMMSE;第二种为使用序贯蒙特卡洛方法(SMC),例如粒子滤波(PF);第三种为卡尔曼滤波类型的算法,例如扩展卡尔曼滤波、无迹卡尔曼滤波、截断无迹卡尔曼滤波、积分卡尔曼滤波等。
本申请发明人在长期研发中发现,现有技术中第一种KF方法只能处理线性系统,而在现实情况中,一个动态系统的真实的过程和观测模型通常都是非线性的,且过程和观测噪声都是非高斯的;第二种例如粒子滤波方法其计算量随着采样粒子数的增加而剧增,一般很难进行实际应用;第三种方法中的扩展卡尔曼滤波会引入较大的近似误差以及引起发散现象,对于截断无迹卡尔曼滤波方法则由于要求当前时刻的观测函数是双向单射的,即要求观测函数具有唯一的反函数,而在实际情况中观测函数通常都不具有唯一的反函数,因此难以进行实际应用;此外扩展卡尔曼滤波、无迹卡尔曼滤波、积分卡尔曼滤波等方法在先验分布方差增大时,非线性的影响增加,观测函数的非线性更加明显,此时滤波性能受到较大影响而使得滤波精度较低。
发明内容
本发明主要解决的技术问题是提供一种目标跟踪方法、系统及扩展截断无迹卡尔曼滤波方法、装置,能够解决观测函数不具有唯一反函数的问题,有效减少目标状态先验分布方差,自适应地根据观测信息的精度进行状态更新,有效提高滤波精度且实用性较高。
为解决上述技术问题,本发明的第一方面是:提供一种目标跟踪方法,包括:对目标进行感测以获得当前目标观测时刻的目标观测向量;根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数;根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;根据第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数;利用目标状态的联合后验概率密度函数对目标状态进行估计,以获得当前目标观测时刻的目标状态估计值;输出当前目标观测时刻的目标状态估计值,以实现对目标的跟踪。
为解决上述技术问题,本发明的第二方面是:提供一种目标跟踪系统,包括:传感器以及扩展截断无迹卡尔曼滤波装置,传感器连接至扩展截断无迹卡尔曼滤波装置,传感器用于对目标进行感测以获得当前目标观测时刻的目标观测向量;扩展截断无迹卡尔曼滤波装置包括:原始先验概率密度函数获取模块,用于根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;第一后验概率密度函数获取模块,用于根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;修正先验概率密度函数获取模块,用于根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数;第二后验概率密度函数获取模块,用于根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;联合后验概率密度函数获取模块,用于根据第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数;目标状态估计模块,用于利用目标状态的联合后验概率密度函数对目标状态进行估计,以获得当前目标观测时刻的目标状态估计值;目标状态估计值输出模块,用于输出当前目标观测时刻的目标状态估计值。
为解决上述技术问题,本发明的第三方面是:提供一种扩展截断无迹卡尔曼滤波方法,包括:根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数;根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;根据原始先验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数,完成扩展截断无迹卡尔曼滤波过程。
其中,根据无迹变换获取当前目标观测目标状态的原始先验概率密度函数的步骤具体包括:根据无迹变换获取N个Sigma点其中i=1,2,...,N;
将Sigma点代入非线性状态转移函数,以获取Sigma点对应的N个预测值具体如下式所示:
χ 0 , k | k - 1 i = f ( χ 0 i ) - - - ( 1 )
其中,f(·)为非线性状态转移函数;
根据预测值获取k时刻目标状态的原始先验概率密度函数对应的均值和协方差,具体如下所示:
x ^ p , 0 , k | k - 1 = Σ i = 1 N w i χ 0 , k | k - 1 i - - - ( 2 )
P p , 0 , k | k - 1 = Q k + Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) T - - - ( 3 )
其中,为Sigma点对应的权值,Qk为过程噪声向量的协方差,k时刻表示当前目标观测时刻,为k时刻目标状态的原始先验概率密度函数对应的均值,Pp,0,k|k-1为k时刻目标状态的原始先验概率密度函数对应的协方差。
其中,根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数的步骤具体包括:将预测值代入非线性观测函数,以获取预测值对应的N个观测值具体如下式所示:
z 0 , k | k - 1 i = h ( χ 0 , k | k - 1 i ) - - - ( 4 )
其中,h(·)为非线性观测函数;
获取N个观测值的均值具体如下式所示:
z ^ 0 , k | k - 1 = Σ i = 1 N w i z 0 , k | k - 1 i - - - ( 5 )
根据均值获取对应的新息协方差Pzz,0,k|k-1以及交叉协方差Pxz,0,k|k-1,具体如下所示:
P zz , 0 , k | k - 1 = R k + Σ i = 1 N w i ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) T - - - ( 6 )
P xz , 0 , k | k - 1 = Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) T - - - ( 7 )
其中,Rk为观测噪声协方差;
根据新息协方差Pzz,0,k|k-1、交叉协方差Pxz,0,k|k-1以及k时刻目标状态的原始先验概率密度函数获取k时刻目标状态的第一后验概率密度函数对应的均值和协方差,具体如下所示:
x ^ u , 0 , k | k = x ^ p , 0 , k | k - 1 + P xz , 0 , k | k - 1 P zz , 0 , k | k - 1 - 1 ( z k - z ^ 0 , k | k - 1 ) - - - ( 8 )
P u , 0 , k | k = P p , 0 , k | k - 1 - P xz , 0 , k | k - 1 P zz , 0 , k | k - 1 - 1 P xz , 0 , k | k - 1 T - - - ( 9 )
其中,为k时刻目标状态的第一后验概率密度函数对应的均值,Pu,0,k|k为k时刻目标状态的第一后验概率密度函数对应的协方差,zk为k时刻的目标观测向量。
其中,根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数的步骤具体包括:根据k时刻目标状态的原始先验概率密度函数获取误差协方差,具体如下式所示:
P xx , 0 , k | k - 1 = Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) T - - - ( 10 )
其中,Pxx,0,k|k-1为误差协方差;
根据统计线性回归理论以及误差协方差Pxx,0,k|k-1获取非线性观测函数的线性回归参数Hk、dk,具体如下所示:
H k = P xz , 0 , k | k - 1 T P xx , 0 , k | k - 1 - 1 - - - ( 11 )
d k = z ^ 0 , k | k - 1 - H k x ^ p , 0 , k | k - 1 - - - ( 12 )
根据线性回归参数Hk、dk以及k时刻的目标观测向量zk获取修正先验概率密度函数的均值和协方差,具体如下所示:
μa,1=Hk*(zk-dk)    (13)
x ^ p , 1 , k | k - 1 = μ a , 1 μ b , 1 - - - ( 14 )
P p , 1 , k | k - 1 = Σ a , 1 Σ ab , 1 ( Σ ab , 1 ) T Σ b , 1 - - - ( 15 )
Σ a , 1 = ∫ ( a k - μ 1 ) ( a k - μ 1 ) T p 1 ( a k ; z k ) d a k = H ~ - 1 R ( H ~ - 1 ) T - - - ( 16 )
μ b , 1 = μ b , 0 + Σ ab , 0 T Σ a , 0 - 1 ( μ a , 1 - μ a , 0 ) - - - ( 17 )
Σ ab , 1 = Σ a , 1 ( Σ a , 0 - 1 ) T Σ ab , 0 - - - ( 18 )
Σ b , 1 = Γ - ( μ b , 1 - μ b , 0 ) ( μ b , 1 - μ b , 0 ) T + Σ ab , 0 T Σ a , 0 - 1 × [ Σ a , 1 + ( μ a , 1 - μ a , 0 ) ( μ a , 1 - μ a , 0 ) T ] ( Σ a , 0 - 1 ) T Σ ab , 0 - - - ( 19 )
Γ = Σ b , 0 - Σ ab , 0 T Σ a , 0 - 1 Σ ab , 0 - - - ( 20 )
x ^ p , 0 , k | k - 1 = Σ i = 1 N w i χ 0 , k | k - 1 i = μ a , 0 μ b , 0 - - - ( 21 )
P p , 0 , k | k - 1 = Q k + Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i x ^ p , 0 , k | k - 1 ) T = Σ a , 0 Σ ab , 0 ( Σ ab , 0 ) T Σ b , 0 - - - ( 22 )
其中,为修正先验概率密度函数的均值,Pp,1,k|k-1为修正先验概率密度函数的协方差,a为目标的状态分量。
其中,根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数的步骤具体包括:
根据无迹变换以及修正先验概率密度函数的均值和协方差Pp,1,k|k-1获取N个Sigma点
将Sigma点代入非线性观测函数h(·),以获取Sigma点对应的N个观测值具体如下式所示:
z 1 , k | k - 1 i = h ( χ 1 , k | k - 1 i ) - - - ( 23 )
获取N个观测值的均值具体如下式所示:
z ^ 1 , k | k - 1 = Σ i = 1 N w i z 1 , k | k - 1 i - - - ( 24 )
根据均值获取对应的新息协方差Pzz,1,k|k-1以及交叉协方差Pxz,1,k|k-1,具体如下所示:
P zz , 1 , k | k - 1 = R k + Σ i = 1 N w i ( z 1 k | k - 1 i - z ^ 1 , k | k - 1 ) ( z 1 , k | k - 1 i - z ^ 1 , k | k - 1 ) T - - - ( 25 )
P xz , 1 , k | k - 1 = Σ i = 1 N w i ( χ 1 , k | k - 1 i - x ^ p , 1 , k | k - 1 ) ( z 1 , k | k - 1 i - z ^ 1 , k | k - 1 ) T - - - ( 26 )
根据新息协方差Pzz,1,k|k-1、交叉协方差Pxz,1,k|k-1以及k时刻的目标观测向量zk获取k时刻目标状态的第二后验概率密度函数对应的均值和协方差,具体如下所示:
x ^ u , 1 , k | k = x ^ p , 1 , k | k - 1 + P xz , 1 , k | k - 1 P zz , 1 , k | k - 1 - 1 ( z k - z ^ 1 , k | k - 1 ) - - - ( 27 )
P u , 1 , k | k = P p , 1 , k | k - 1 - P xz , 1 , k | k - 1 P zz , 1 , k | k - 1 - 1 P xz , 1 , k | k - 1 T - - - ( 28 )
其中,为k时刻目标状态的第二后验概率密度函数对应的均值,Pu,1,k|k为k时刻目标状态的第二后验概率密度函数对应的协方差。
其中,根据原始先验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数的步骤具体包括:
根据第一后验概率密度函数以及第二后验概率密度函数各自对应的均值获取目标状态估计权值,具体如下所示:
μ 0 ( x ^ u , 0 , k | k ) = 1 | P zz , 0 , k | k - 1 | · exp ( - ( z k - h ( x ^ u , 0 , k | k ) ) 2 2 ) - - - ( 29 )
μ 1 ( x ^ u , 1 , k | k ) = 1 | P zz , 1 , k | k - 1 | · exp ( - ( z k - h ( x ^ u , 1 , k | k ) ) 2 2 ) - - - ( 30 )
α = μ 1 ( x ^ u , 1 , k | k ) μ 0 ( x ^ u , 0 , k | k ) + μ 1 ( x ^ u , 1 , k | k ) - - - ( 31 )
其中,α为目标状态估计权值;
根据第一后验概率密度函数、第二后验概率密度函数各自对应的均值、协方差以及目标状态估计权值α获取k时刻目标状态的联合后验概率密度函数,具体如下所示:
p ( x k | z 1 : k ) = N ( x ^ k | k , P k | k ) - - - ( 32 )
x ^ k | k = α k · x ^ u , 1 , k | k + ( 1 - α k ) · x ^ u , 0 , k | k - - - ( 33 )
P k | k = α k · [ P u , 1 , k | k + ( x ^ u , 1 , k | k - x ^ k | k ) ( x ^ u , 1 , k | k - x ^ k | k ) T ] + ( 1 - α k ) · [ P u , 0 , k | k + ( x ^ u , 0 , k | k - x ^ k | k ) ( x ^ u , 0 , k | k - x ^ k | k ) T ] - - - ( 34 )
其中,p(xk|z1:k)为k时刻目标状态的联合后验概率密度函数,为k时刻目标状态的联合后验概率密度函数的均值,Pk|k为k时刻目标状态的联合后验概率密度函数的协方差。
为解决上述技术问题,本发明的第四方面是:提供一种扩展截断无迹卡尔曼滤波装置,包括:原始先验概率密度函数获取模块,用于根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;第一后验概率密度函数获取模块,用于根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;修正先验概率密度函数获取模块,用于根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数;第二后验概率密度函数获取模块,用于根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;联合后验概率密度函数获取模块,用于根据第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数。
本发明的有益效果是:区别于现有技术的情况,本发明根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数,能够解决观测函数不具有唯一反函数的问题,有效减少目标状态先验分布方差;根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;根据原始先验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数,实现自适应地根据观测信息的精度进行状态更新,有效提高滤波精度且实用性较高。
附图说明
图1是本发明扩展截断无迹卡尔曼滤波方法一实施方式的流程图;
图2是本发明扩展截断无迹卡尔曼滤波方法一实施方式中根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数的流程图;
图3是本发明扩展截断无迹卡尔曼滤波方法一实施方式中根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数的流程图;
图4是本发明扩展截断无迹卡尔曼滤波方法一实施方式中根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数的流程图;
图5是本发明扩展截断无迹卡尔曼滤波方法一实施方式中根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数的流程图;
图6是本发明扩展截断无迹卡尔曼滤波方法一实施方式中根据第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数的流程图;
图7是本发明扩展截断无迹卡尔曼滤波方法一实施方式中EKF、UKF、QKF、PF和ETUKF的均方根误差对比图;
图8是本发明扩展截断无迹卡尔曼滤波方法一实施方式中EKF、UKF、QKF、PF和ETUKF在不同观测噪声方差下的均方根误差对比图;
图9是本发明扩展截断无迹卡尔曼滤波装置一实施方式的原理框图;
图10是本发明目标跟踪系统一实施方式的结构示意图。
具体实施方式
下面将结合本发明实施方式中的附图,对本发明实施方式中的技术方案进行清楚、完整地描述,显然,所描述的实施方式仅仅是本发明一部分实施方式,而不是全部的实施方式。基于本发明中的实施方式,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式,均属于本发明保护的范围。
本发明的扩展截断无迹卡尔曼滤波ETUKF为从一系列的不完全及包含噪声的目标观测向量中,对目标在某一时刻的状态进行估计。本发明扩展截断无迹卡尔曼滤波方法的基本理论具体如下所述:
假设k时刻目标的状态向量为xk=[ak T,bk T]T并且nx=na+nb,ak、bk为目标的状态分量,其中ak具体可为目标的位置信息,bk具体可为目标的速度信息,对应的观测方程可为:
zk=h(ak)+vk
进一步作如下两个假设:(1)观测函数为连续的、单射的函数。(2)噪声的概率密度函数为有界的、连通的,即
pη(v)=0
η ∉ I n ⋐ R n z
其中,Iη表示一个nx维的连通区域。基于上述假设(2),在目标状态已知的条件下观测似然函数为:
p ( z k | x k ) = p ( z k | a k ) = p η ( z k - h ( a k ) ) χ I η ( z k - h ( a k ) ) - - - ( a )
其中,表示子集Iη的指示函数。
根据上述假设(1),式(a)可以进一步表示为:
p ( z k | x k ) = p v ( z k - h ( a k ) ) χ I x ( z ) ( x k )
其中,
I x ( z k ) = { x k | x k = [ ( h - 1 ( z k - v k ) ) T , b k T ] T , v k ∈ I v , b k ∈ R n b } = I a ( z k ) × R n b - - - ( b )
根据贝叶斯准则,xk的后验概率密度函数对应为:
p ( x k | z k ) ∝ p ( z k | x k ) · p 0 ( x k ) = p v ( z k - h ( a k ) ) χ I x ( z ) ( x k ) p 0 ( x k ) ∝ p ( z k | x k ) · p 1 ( x k ; z k )
其中:
p 1 ( x k ; z k ) = 1 ϵ 1 p 0 ( x k ) χ I x ( z ) ( x k ) - - - ( c )
其中,ε1为标准化常数,从式(c)可以看出,基于当前观测时刻的目标观测向量zk,式(c)构建了一种修正的先验概率密度函数p1(·),可以有效减少原始先验概率密度函数p0(·)的方差。
基于上述基本理论,请参阅图1,本发明扩展截断无迹卡尔曼滤波方法一实施方式包括以下步骤:
步骤S11:根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;
无迹变换(Unscented Transformation,UT)是用于计算经过非线性变换的随机变量统计的一种方法,无迹变换不需要对非线性状态和测量模型进行线性化,而是对状态向量的概率密度函数进行近似化,近似化后的概率密度函数仍然是高斯的,但它表现为一系列选取好的Sigma采样点,其获取Sigma点及Sigma点对应权值的具体方法为:
假设xk|k为一个nx维随机向量,xk|k的均值和协方差分别为Pk|k,根据下述公式获取(2na+1)个Sigma点χi和χi对应的权值wi
χ 0 = x ^ k | k , i = 0 χ i = x ^ k | k + ( ( n a + κ ) P k | k ) i , i = 1 , . . . , n a χ i + n a = x ^ k | k - ( ( n a + κ ) P k | k ) i , i = 1 , . . . , n a
w 0 = κ / ( n a + κ ) i = 0 w i = 1 / [ 2 ( n a + κ ) ] i = 1 , . . . , n a w i + n x = 1 / [ 2 ( n a + κ ) ] i = 1 , . . . , n a
其中,κ是一个尺度参数,κ为满足(na+κ)≠0的任何数值,是均方根矩阵(na+κ)Pk|k的第i行或第j列,na为状态向量的维数。
本步骤根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数,请参阅图2,具体包括以下子步骤:
子步骤S111:根据无迹变换获取N个Sigma点;
假设在上一目标观测时刻(k-1时刻)系统的初始状态均值以及协方差Pk-1|k-1为已知,即k-1时刻目标状态的联合后验概率密度函数的均值和协方差Pk-1|k-1为已知,基于上述均值和协方差Pk-1|k-1根据无迹变换获取N个Sigma点其中i=1,2,...,N,N=2na+1,na=nx+nv+nw,nx表示目标状态x的维数,nv表示观测噪声v的维数,nw表示过程噪声w的维数。
子步骤S112:将Sigma点代入非线性状态转移函数,以获取Sigma点对应的N个预测值;
将上述N个Sigma点分别代入非线性状态转移函数,以获取N个Sigma点对应的N个预测值具体如下式(1)所示:
χ 0 , k | k - 1 i = f ( χ 0 i ) - - - ( 1 )
其中,f(·)为非线性状态转移函数。
子步骤S113:根据预测值获取k时刻目标状态的原始先验概率密度函数对应的均值和协方差。
根据N个预测值获取k时刻目标状态的原始先验概率密度函数对应的均值和协方差,具体如下式(2)、(3)所示:
x ^ p , 0 , k | k - 1 = Σ i = 1 N w i χ 0 , k | k - 1 i - - - ( 2 )
P p , 0 , k | k - 1 = Q k + Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) T - - - ( 3 )
其中,为N个Sigma点各自对应的权值,Qk为过程噪声向量的协方差,k时刻表示当前目标观测时刻,为k时刻目标状态的原始先验概率密度函数对应的均值,Pp,0,k|k-1为k时刻目标状态的原始先验概率密度函数对应的协方差。
步骤S12:根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;
本步骤进一步基于上述原始先验概率密度函数进行观测更新以获取当前目标观测时刻(k时刻)目标状态的第一后验概率密度函数,请参阅图3,具体包括以下子步骤:
子步骤S121:将预测值代入非线性观测函数,以获取预测值对应的N个观测值;
将上述子步骤S112获取的预测值代入非线性观测函数,以获取N个预测值对应的N个观测值具体如下式(4)所示:
z 0 , k | k - 1 i = h ( χ 0 , k | k - 1 i ) - - - ( 4 )
其中,h(·)为非线性观测函数。
子步骤S122:获取N个观测值的均值;
获取N个观测值的均值具体如下式(5)所示:
z ^ 0 , k | k - 1 = Σ i = 1 N w i z 0 , k | k - 1 i - - - ( 5 )
子步骤S123:根据均值获取对应的新息协方差以及交叉协方差;
根据观测值的均值获取对应的新息协方差Pzz,0,k|k-1以及交叉协方差Pxz,0,k|k-1,具体如下式(6)、(7)所示:
P zz , 0 , k | k - 1 = R k + Σ i = 1 N w i ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) T - - - ( 6 )
P xz , 0 , k | k - 1 = Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) T - - - ( 7 )
其中,Rk为观测噪声协方差。
子步骤S124:根据新息协方差、交叉协方差以及k时刻目标状态的原始先验概率密度函数获取k时刻目标状态的第一后验概率密度函数对应的均值和协方差。
根据新息协方差Pzz,0,k|k-1、交叉协方差Pxz,0,k|k-1以及k时刻(当前目标观测时刻)目标状态的原始先验概率密度函数获取k时刻目标状态的第一后验概率密度函数对应的均值和协方差,具体如下所示:
x ^ u , 0 , k | k = x ^ p , 0 , k | k - 1 + P xz , 0 , k | k - 1 P zz , 0 , k | k - 1 - 1 ( z k - z ^ 0 , k | k - 1 ) - - - ( 8 )
P u , 0 , k | k = P p , 0 , k | k - 1 - P xz , 0 , k | k - 1 P zz , 0 , k | k - 1 - 1 P xz , 0 , k | k - 1 T - - - ( 9 )
其中,为k时刻目标状态的第一后验概率密度函数对应的均值,Pu,0,k|k为k时刻目标状态的第一后验概率密度函数对应的协方差,zk为k时刻的目标观测向量。
步骤S13:根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数;
为了获取修正先验概率密度函数p1(·),对应作如下三个近似假设:(1)观测函数h(·)是局部线性的;(2)目标的状态分量a的边缘先验概率密度p0(a)在区域Ia(z)是常数;(3)观测噪声w满足修改后的噪声和真实噪声有相同的一二阶矩,E[w]=0,cov[w]=R。
根据假设1,由于观测函数h(·)是局部线性的,在点对观测函数h(·)进行一阶泰勒级数展开,定义为ak的最大似然估计,并令这是因为根据当前目标观测,是最可能的目标状态,最大化似然函数可得到a,因此有:
h ( a ) = h ( a ~ ( z ) ) + H ~ ( a - a ~ ( z ) )
其中 为雅克比矩阵在的取值。
代入上述基本理论中的式(b)可得:
I a ( z k ) = { a k | a k = a ~ k ( z k ) - H ~ - 1 w k , w k ∈ I η }
接着上述式(c)转化为:
p 1 ( x k ; z k ) = p 1 ( a k , b k ; z k ) = 1 ϵ 1 χ I a ( z ) ( a k ) p 0 ( b k | a k ) p 0 ( a k ) - - - ( d )
由于上述假设p0(ak)在区域Ia(z)中为常数,因此式(d)可写为:
p 1 ( a k , b k ; z k ) = 1 ϵ 2 χ I a ( z ) ( a k ) p 0 ( b k | a k ) - - - ( e )
其中,ε2为归一化常数。根据式(e),修正先验概率密度函数p1(·)中目标的状态分量a对应的均值μa,1和协方差Σa,1可以近似如下:
μ a , 1 = ∫ a k p 1 ( a k ; z k ) da k = a ~ k ( z k )
Σ a , 1 = ∫ ( a k - μ 1 ) ( a k - μ 1 ) T p 1 ( a k ; z k ) d a k = H ~ - 1 R ( H ~ - 1 ) T
其中,a为目标的状态分量。
另外,原始先验概率密度函数p0(xk)=p0(ak,bk)的均值和协方差Pp,0,k|k-1可以分解为:
x ^ p , 0 , k | k - 1 = Σ i = 1 N w i χ 0 , k | k - 1 i = μ a , 0 μ b , 0
P p , 0 , k | k - 1 = Q k + Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i x ^ p , 0 , k | k - 1 ) T = Σ a , 0 Σ ab , 0 ( Σ ab , 0 ) T Σ b , 0
因此修正先验概率密度函数p1(·)的均值和协方差Pp,1,k|k-1分别如下所示:
x ^ p , 1 , k | k - 1 = μ a , 1 μ b , 1
P p , 1 , k | k - 1 = Σ a , 1 Σ ab , 1 ( Σ ab , 1 ) T Σ b , 1
由上述推导过程可以看出,要获取修正先验概率密度函数p1(·),观测函数h(·)必须为双向单射的,即要求观测函数h(·)具有唯一的反函数,而在实际情况中观测函数h(·)通常不具有唯一的反函数。针对上述观测函数为非双向单射的问题,本发明根据统计线性回归理论获取修正先验概率密度函数,请参阅图4,其具体包括以下子步骤:
子步骤S131:根据k时刻目标状态的原始先验概率密度函数获取误差协方差;
根据k时刻目标状态的原始先验概率密度函数获取误差协方差,具体如下式(10)所示:
P xx , 0 , k | k - 1 = Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) T - - - ( 10 )
其中,Pxx,0,k|k-1为误差协方差。
子步骤S132:根据统计线性回归理论以及误差协方差获取非线性观测函数的线性回归参数;
为了获取上述反函数h-1(zk),目标为将非线性函数近似成如下线性估计的形式其中Hk、dk分别为由最小均方根误差准则确定的线性回归参数。本实施方式根据统计线性回归理论以及误差协方差Pxx,0,k|k-1获取非线性观测函数的线性回归参数Hk、dk的过程具体如下所述。
对于线性回归参数Hk、dk有:
{Hk,dk}=argminE(ek Tek)    (f)
其中,ek为线性化误差,将ek代入式(f),对dk求导并令其为0,得到:
( - 2 ) E ( z k - H k x k - d k ) = 0 ⇔ d k = z ‾ k - H k x ‾ k
其中, x ‾ k = E ( x k ) = x ^ p , 0 , k | k - 1 , z ‾ k = E ( z k ) = z ^ 0 , k | k - 1 . 将dk代入式(f),得到:
e T e = [ ( z k - z ‾ k ) - H k ( x k - x ‾ k ) ] T [ ( z k - z ‾ k ) - H k ( x k - x ‾ k ) ]
对Hk求偏导并令其为0,得到:
( - 2 ) E { [ ( z k - z ‾ k ) - H k ( x k - x ‾ k ) ] [ z k - z ‾ k ] ) } = 0
进一步可得线性回归参数Hk、dk,具体如下式(11)、(12)所示:
H k = P xz , 0 , k | k - 1 T P xx , 0 , k | k - 1 - 1 - - - ( 11 )
d k = z ^ 0 , k | k - 1 - H k x ^ p , 0 , k | k - 1 - - - ( 12 )
其中,Pxz,0,k|k-1即为上述公式(7)中的交叉协方差,Pxx,0,k|k-1即为上述公式(10)中的误差协方差,即为上述公式(5)中的均值,即为上述公式(2)中的均值。
子步骤S133:根据线性回归参数以及k时刻的目标观测向量获取修正先验概率密度函数的均值和协方差。
由上述论述可知:
μa,1=∫akp1(ak;zk)dak=h-1(zk)
进一步根据线性回归参数Hk、dk以及k时刻的目标观测向量zk获取修正先验概率密度函数的均值和协方差,具体如下式(13-22)所示:
μa,1=Hk*(zk-dk)    (13)
x ^ p , 1 , k | k - 1 = μ a , 1 μ b , 1 - - - ( 14 )
P p , 1 , k | k - 1 = Σ a , 1 Σ ab , 1 ( Σ ab , 1 ) T Σ b , 1 - - - ( 15 )
Σ a , 1 = ∫ ( a k - μ 1 ) ( a k - μ 1 ) T p 1 ( a k ; z k ) d a k = H ~ - 1 R ( H ~ - 1 ) T - - - ( 16 )
μ b , 1 = μ b , 0 + Σ ab , 0 T Σ a , 0 - 1 ( μ a , 1 - μ a , 0 ) - - - ( 17 )
Σ ab , 1 = Σ a , 1 ( Σ a , 0 - 1 ) T Σ ab , 0 - - - ( 18 )
Σ b , 1 = Γ - ( μ b , 1 - μ b , 0 ) ( μ b , 1 - μ b , 0 ) T + Σ ab , 0 T Σ a , 0 - 1 × [ Σ a , 1 + ( μ a , 1 - μ a , 0 ) ( μ a , 1 - μ a , 0 ) T ] ( Σ a , 0 - 1 ) T Σ ab , 0 - - - ( 19 )
Γ = Σ b , 0 - Σ ab , 0 T Σ a , 0 - 1 Σ ab , 0 - - - ( 20 )
x ^ p , 0 , k | k - 1 = Σ i = 1 N w i χ 0 , k | k - 1 i = μ a , 0 μ b , 0 - - - ( 21 )
P p , 0 , k | k - 1 = Q k + Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i x ^ p , 0 , k | k - 1 ) T = Σ a , 0 Σ ab , 0 ( Σ ab , 0 ) T Σ b , 0 - - - ( 22 )
其中,为修正先验概率密度函数的均值,Pp,1,k|k-1为修正先验概率密度函数的协方差,a为目标的状态分量。
步骤S14:根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;
本步骤进一步基于上述修正先验概率密度函数进行观测更新以获取当前目标观测时刻(k时刻)目标状态的第二后验概率密度函数,请参阅图5,本步骤具体包括以下子步骤:
子步骤S141:根据无迹变换以及修正先验概率密度函数的均值和协方差获取N个Sigma点;
根据无迹变换以及修正先验概率密度函数的均值和协方差Pp,1,k|k-1获取N个Sigma点
子步骤S142:将Sigma点代入非线性观测函数,以获取Sigma点对应的N个观测值;
将Sigma点代入非线性观测函数h(·),以获取Sigma点对应的N个观测值具体如下式(23)所示:
z 1 , k | k - 1 i = h ( χ 1 , k | k - 1 i ) - - - ( 23 )
子步骤S143:获取N个观测值的均值;
获取N个观测值的均值具体如下式(24)所示:
z ^ 1 , k | k - 1 = Σ i = 1 N w i z 1 , k | k - 1 i - - - ( 24 )
其中,公式(24)中的wi同样指代上述公式(2)、(3)等所示的权值wi
子步骤S144:根据均值获取对应的新息协方差以及交叉协方差;
根据均值获取对应的新息协方差Pzz,1,k|k-1以及交叉协方差Pxz,1,k|k-1,具体如下式(25)、(26)所示:
P zz , 1 , k | k - 1 = R k + Σ i = 1 N w i ( z 1 k | k - 1 i - z ^ 1 , k | k - 1 ) ( z 1 , k | k - 1 i - z ^ 1 , k | k - 1 ) T - - - ( 25 )
P xz , 1 , k | k - 1 = Σ i = 1 N w i ( χ 1 , k | k - 1 i - x ^ p , 1 , k | k - 1 ) ( z 1 , k | k - 1 i - z ^ 1 , k | k - 1 ) T - - - ( 26 )
子步骤S145:根据新息协方差、交叉协方差以及k时刻的目标观测向量获取k时刻目标状态的第二后验概率密度函数对应的均值和协方差。
根据新息协方差Pzz,1,k|k-1、交叉协方差Pxz,1,k|k-1以及k时刻的目标观测向量zk获取k时刻目标状态的第二后验概率密度函数对应的均值和协方差,具体如下式(27)、(28)所示:
x ^ u , 1 , k | k = x ^ p , 1 , k | k - 1 + P xz , 1 , k | k - 1 P zz , 1 , k | k - 1 - 1 ( z k - z ^ 1 , k | k - 1 ) - - - ( 27 )
P u , 1 , k | k = P p , 1 , k | k - 1 - P xz , 1 , k | k - 1 P zz , 1 , k | k - 1 - 1 P xz , 1 , k | k - 1 T - - - ( 28 )
其中,为k时刻目标状态的第二后验概率密度函数对应的均值,Pu,1,k|k为k时刻目标状态的第二后验概率密度函数对应的协方差。
步骤S15:根据原始先验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数。
本步骤根据步骤S12获取的第一后验概率密度函数以及步骤S14获取的第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数,请参阅图6,本步骤具体包括以下子步骤:
子步骤S151:根据第一后验概率密度函数以及第二后验概率密度函数各自对应的均值获取目标状态估计权值;
当目标观测向量zk精度较高时,基于修正先验概率密度函数得到的目标状态的估计结果更可信;当目标观测向量zk精度较低时,基于原始先验概率密度函数得到的目标状态的估计结果更可信。基于此,本发明通过一目标状态估计权值体现目标观测精度高低对目标状态的估计结果的影响:当目标观测向量zk精度较高时,目标状态估计权值更接近于1;当目标观测zk精度较低时,目标状态估计权值更接近于0。根据这一原则,考虑当前的估计结果和目标观测,本发明定义如下的高斯模糊隶属函数:
μ ( x ) = K · exp ( - ( z k - h ( x ) ) 2 2 σ 2 )
其中,Κ为一已知参数,σ2表示新息方差。
根据第一后验概率密度函数对应的均值以及第二后验概率密度函数对应的均值获取目标状态估计权值,具体如下式(29-31)所示:
μ 0 ( x ^ u , 0 , k | k ) = 1 | P zz , 0 , k | k - 1 | · exp ( - ( z k - h ( x ^ u , 0 , k | k ) ) 2 2 ) - - - ( 29 )
μ 1 ( x ^ u , 1 , k | k ) = 1 | P zz , 1 , k | k - 1 | · exp ( - ( z k - h ( x ^ u , 1 , k | k ) ) 2 2 ) - - - ( 30 )
α = μ 1 ( x ^ u , 1 , k | k ) μ 0 ( x ^ u , 0 , k | k ) + μ 1 ( x ^ u , 1 , k | k ) - - - ( 31 )
其中,α为目标状态估计权值。由式(29-31)可以看出,当zk的差别越小时,表示估计结果越接近真实的目标状态,此时目标状态估计权值αk相应变小;反之,目标状态估计权值αk相应变大。即目标状态估计权值αk自适应地随着目标观测向量的精度进行调整。
子步骤S152:根据第一后验概率密度函数、第二后验概率密度函数各自对应的均值、协方差以及目标状态估计权值获取k时刻目标状态的联合后验概率密度函数。
根据第一后验概率密度函数对应的均值协方差Pu,0,k|k、第二后验概率密度函数对应的均值协方差Pu,1,k|k以及目标状态估计权值α获取k时刻目标状态的联合后验概率密度函数,具体如下式(32-34)所示:
p ( x k | z 1 : k ) = N ( x ^ k | k , P k | k ) - - - ( 32 )
x ^ k | k = α k · x ^ u , 1 , k | k + ( 1 - α k ) · x ^ u , 0 , k | k - - - ( 33 )
P k | k = α k · [ P u , 1 , k | k + ( x ^ u , 1 , k | k - x ^ k | k ) ( x ^ u , 1 , k | k - x ^ k | k ) T ] + ( 1 - α k ) · [ P u , 0 , k | k + ( x ^ u , 0 , k | k - x ^ k | k ) ( x ^ u , 0 , k | k - x ^ k | k ) T ] - - - ( 34 )
其中,p(xk|z1:k)为k时刻目标状态的联合后验概率密度函数,为k时刻目标状态的联合后验概率密度函数的均值,Pk|k为k时刻目标状态的联合后验概率密度函数的协方差。
上述步骤S11对应为本发明扩展截断无迹卡尔曼滤波方法的时间更新阶段,步骤S12对应为本发明扩展截断无迹卡尔曼滤波方法基于原始先验概率密度函数的观测更新阶段,步骤S13-14为本发明扩展截断无迹卡尔曼滤波方法基于修正先验概率密度函数的观测更新阶段,步骤S15为本发明扩展截断无迹卡尔曼滤波方法的联合状态更新阶段。
以下将以一个例子对本发明的ETUKF方法的性能进行评估以及与现有的扩展卡尔曼滤波(EKF)、无迹卡尔曼滤波(UKF)、积分卡尔曼滤波(QKF)和粒子滤波(PF)的性能进行对比,具体如下所述。
该例子考虑单变量非平稳增长模型(UNGM),该模型的过程模型是高度非线性的,观测模型是非平稳的。非线性过程模型和观测模型可以写为:
x k = α x k - 1 + β x k - 1 1 + x k - 1 2 + γ cos ( 1.2 ( k - 1 ) ) + n k
z k = φ 2 x k 2 + v k k ≤ 30 φ 1 x k 3 - 2 + v k k > 30
其中nk是均值为零方差为1的高斯噪声,vk是均值为零方差为0.01的高斯噪声,β=25、α=0.5、γ=8、φ1=0.2、φ2=0.05。在每次蒙特卡洛仿真中,假设状态x0的初始分布式是[0 1]上的均匀分布。
图7给出了本发明的ETUKF与现有的EKF、UKF、QKF、PF的均方根误差对比结果。从图7中可以看出ETUKF的性能明显优于EKF、UKF和QKF,ETUKF的性能与PF接近。导致EKF、UKF和QKF的性能较差的一个原因是高度非线性的过程模型和非平稳的观测模型所产生的近似误差的增加。
图8给出了不同观测噪声方差下本发明的ETUKF与现有的EKF、UKF、QKF、PF的均方根误差对比结果。由图8中可以看出,当观测信息量非常丰富或不丰富时,本发明的ETUKF在所有情况下都具有鲁棒性。ETUKF的性能明显优于EKF、UKF和QKF,ETUKF的性能与PF接近。
此外,下表1给出了上述各个滤波方法的计算时间。
表1 不同滤波方法计算时间的对比
上述各个方法各执行了100次蒙特卡洛仿真。由表1可以看出,对于一维UNGM,PF的计算时间比EKF、UKF、QKF和ETUKF高得多,ETUKF和QKF所用时间接近。然而,在四维BOT的情况下,QKF的计算时间(117.847s)比ETUKF高得多(30.913s),造成该情况的一个主要原因是:QKF遭受维数灾难,其高斯-厄米特积分点的数目随着状态向量的维数的增加呈指数增长。在各个方法中,EKF计算最快。但是,对照上表1和图7或图8,可以发现在计算效率这方面本发明的ETUKF优于现有的EKF、UKF、QKF、PF。
可以理解,本发明扩展截断无迹卡尔曼滤波方法一实施方式根据当前目标观测向量对原始先验概率密度函数进行修正,能够减少目标状态先验分布方差,有效解决观测函数的非线性对滤波精度的影响,提高滤波精度;根据统计线性回归理论而获取修正先验概率密度函数,能够解决传统滤波算法进行观测更新时要求观测函数必须具有唯一反函数的问题,能够进行实际应用;在联合后验概率密度函数中引入目标状态估计权值,实现自适应地根据观测信息的精度进行状态更新,能够有效提高滤波精度,提高对目标的跟踪性能;此外,本发明的扩展截断无迹卡尔曼滤波方法并不基于粒子滤波,计算量较简单,实用性较高。
请参阅图9,本发明实施方式还提供一种扩展截断无迹卡尔曼滤波装置,该扩展截断无迹卡尔曼滤波装置具体包括以下各个模块:
原始先验概率密度函数获取模块21,用于根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数,其中该原始先验概率密度函数获取模块21获取原始先验概率密度函数所采用的具体方法可以参阅上一实施方式中的步骤S11。
第一后验概率密度函数获取模块22,用于根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数,其中该第一后验概率密度函数获取模块22获取第一后验概率密度函数所采用的具体方法可以参阅上一实施方式中的步骤S12。
修正先验概率密度函数获取模块23,用于根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数,其中该修正先验概率密度函数获取模块23获取修正先验概率密度函数所采用的具体方法可以参阅上一实施方式中的步骤S13。
第二后验概率密度函数获取模块24,用于根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数,其中该第二后验概率密度函数获取模块24获取第二后验概率密度函数所采用的具体方法可以参阅上一实施方式中的步骤S14。
联合后验概率密度函数获取模块25,用于根据第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数,其中该联合后验概率密度函数获取模块25获取联合后验概率密度函数所采用的具体方法可以参阅上一实施方式中的步骤S15。
此外,在其他实施方式中也可采用其他模块框架结构实现本发明的扩展截断无迹卡尔曼滤波方法而不局限于本实施方式所提供的扩展截断无迹卡尔曼滤波装置对应的模块框架结构,此处不作过多限制。
本发明实施方式还提供一种目标跟踪方法,该目标跟踪方法包括:
对目标进行感测以获得当前目标观测时刻的目标观测向量。
根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数。
根据原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数。
根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正原始先验概率密度函数,以获取修正先验概率密度函数。
根据修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数。
根据第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数。
利用目标状态的联合后验概率密度函数对目标状态进行估计,以获得当前目标观测时刻的目标状态估计值。
输出当前目标观测时刻的目标状态估计值,以实现对飞机、航空飞行器、车辆等目标的跟踪。
请参阅图10,本发明实施方式还提供一种目标跟踪系统,该目标跟踪系统包括传感器31以及扩展截断无迹卡尔曼滤波装置32,传感器31连接至扩展截断无迹卡尔曼滤波装置32。
传感器31用于对目标进行感测以获得当前目标观测时刻的目标观测向量,其中传感器31具体可为红外、雷达等。
扩展截断无迹卡尔曼滤波装置32处理来自传感器31的观测数据,其处理过程可参阅上述的目标跟踪方法或扩展截断无迹卡尔曼滤波方法实施方式,在此不再赘述。
以上所述仅为本发明的实施方式,并非因此限制本发明的专利范围,凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换,或直接或间接运用在其他相关的技术领域,均同理包括在本发明的专利保护范围内。

Claims (9)

1.一种目标跟踪方法,其特征在于,包括:
对目标进行感测以获得当前目标观测时刻的目标观测向量;
根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;
根据所述原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;
根据统计线性回归理论以及所述当前目标观测时刻的目标观测向量修正所述原始先验概率密度函数,以获取修正先验概率密度函数;
根据所述修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;
根据所述第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数;
利用所述目标状态的联合后验概率密度函数对目标状态进行估计,以获得当前目标观测时刻的目标状态估计值;
输出所述当前目标观测时刻的目标状态估计值,以实现对目标的跟踪。
2.一种目标跟踪系统,其特征在于,包括:
传感器以及扩展截断无迹卡尔曼滤波装置,所述传感器连接至扩展截断无迹卡尔曼滤波装置,所述传感器用于对目标进行感测以获得当前目标观测时刻的目标观测向量;
所述扩展截断无迹卡尔曼滤波装置包括:
原始先验概率密度函数获取模块,用于根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;
第一后验概率密度函数获取模块,用于根据所述原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;
修正先验概率密度函数获取模块,用于根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正所述原始先验概率密度函数,以获取修正先验概率密度函数;
第二后验概率密度函数获取模块,用于根据所述修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;
联合后验概率密度函数获取模块,用于根据所述第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数;
目标状态估计模块,用于利用所述目标状态的联合后验概率密度函数对目标状态进行估计,以获得当前目标观测时刻的目标状态估计值;
目标状态估计值输出模块,用于输出所述当前目标观测时刻的目标状态估计值。
3.一种扩展截断无迹卡尔曼滤波方法,其特征在于,包括:
根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;
根据所述原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;
根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正所述原始先验概率密度函数,以获取修正先验概率密度函数;
根据所述修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;
根据所述原始先验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数,完成所述扩展截断无迹卡尔曼滤波过程。
4.根据权利要求3所述的方法,其特征在于,所述根据无迹变换获取当前目标观测目标状态的原始先验概率密度函数的步骤具体包括:
根据无迹变换获取N个Sigma点其中i=1,2,...,N;
将所述Sigma点代入非线性状态转移函数,以获取Sigma点对应的N个预测值具体如下式所示:
χ 0 , k | k - 1 i = f ( χ 0 i ) - - - ( 1 )
其中,f(·)为所述非线性状态转移函数;
根据所述预测值获取k时刻目标状态的原始先验概率密度函数对应的均值和协方差,具体如下所示:
x ^ p , 0 , k | k - 1 = Σ i = 1 N w i χ 0 , k | k - 1 i - - - ( 2 )
P p , 0 , k | k - 1 = Q k + Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) T - - - ( 3 )
其中,为Sigma点对应的权值,Qk为过程噪声向量的协方差,k时刻表示当前目标观测时刻,为k时刻目标状态的原始先验概率密度函数对应的均值,Pp,0,k|k-1为k时刻目标状态的原始先验概率密度函数对应的协方差。
5.根据权利要求4所述的方法,其特征在于,所述根据所述原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数的步骤具体包括:
将所述预测值代入非线性观测函数,以获取预测值对应的N个观测值具体如下式所示:
z 0 , k | k - 1 i = h ( χ 0 , k | k - 1 i ) - - - ( 4 )
其中,h(·)为所述非线性观测函数;
获取所述N个观测值的均值具体如下式所示:
z ^ 0 , k | k - 1 = Σ i = 1 N w i z 0 , k | k - 1 i - - - ( 5 )
根据所述均值获取对应的新息协方差Pzz,0,k|k-1以及交叉协方差Pxz,0,k|k-1,具体如下所示:
P zz , 0 , k | k - 1 = R k + Σ i = 1 N w i ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) T - - - ( 6 )
P xz , 0 , k | k - 1 = Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( z 0 , k | k - 1 i - z ^ 0 , k | k - 1 ) T - - - ( 7 )
其中,Rk为观测噪声协方差;
根据所述新息协方差Pzz,0,k|k-1、交叉协方差Pxz,0,k|k-1以及k时刻目标状态的原始先验概率密度函数获取k时刻目标状态的第一后验概率密度函数对应的均值和协方差,具体如下所示:
x ^ u , 0 , k | k = x ^ p , 0 , k | k - 1 + P xz , 0 , k | k - 1 P zz , 0 , k | k - 1 - 1 ( z k - z ^ 0 , k | k - 1 ) - - - ( 8 )
P u , 0 , k | k = P p , 0 , k | k - 1 - P xz , 0 , k | k - 1 P zz , 0 , k | k - 1 - 1 P xz , 0 , k | k - 1 T - - - ( 9 )
其中,为所述k时刻目标状态的第一后验概率密度函数对应的均值,Pu,0,k|k为所述k时刻目标状态的第一后验概率密度函数对应的协方差,zk为k时刻的目标观测向量。
6.根据权利要求5所述的方法,其特征在于,所述根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正所述原始先验概率密度函数,以获取修正先验概率密度函数的步骤具体包括:
根据所述k时刻目标状态的原始先验概率密度函数获取误差协方差,具体如下式所示:
P xx , 0 , k | k - 1 = Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) T - - - ( 10 )
其中,Pxx,0,k|k-1为所述误差协方差;
根据统计线性回归理论以及误差协方差Pxx,0,k|k-1获取非线性观测函数的线性回归参数Hk、dk,具体如下所示:
H k = P xz , 0 , k | k - 1 T P xx , 0 , k | k - 1 - 1 - - - ( 11 )
d k = z ^ 0 , k | k - 1 - H k x ^ p , 0 , k | k - 1 - - - ( 12 )
根据所述线性回归参数Hk、dk以及k时刻的目标观测向量zk获取修正先验概率密度函数的均值和协方差,具体如下所示:
μa,1=Hk*(zk-dk)    (13)
x ^ p , 1 , k | k - 1 = μ a , 1 μ b , 1 - - - ( 14 )
P p , 1 , k | k - 1 = Σ a , 1 Σ ab , 1 ( Σ ab , 1 ) T Σ b , 1 - - - ( 15 )
Σ a , 1 = ∫ ( a k - μ 1 ) ( a k - μ 1 ) T p 1 ( a k ; z k ) d a k = H ~ - 1 R ( H ~ - 1 ) T - - - ( 16 )
μ b , 1 = μ b , 0 + Σ ab , 0 T Σ a , 0 - 1 ( μ a , 1 - μ a , 0 ) - - - ( 17 )
Σ ab , 1 = Σ a , 1 ( Σ a , 0 - 1 ) T Σ ab , 0 - - - ( 18 )
Σ b , 1 = Γ - ( μ b , 1 - μ b , 0 ) ( μ b , 1 - μ b , 0 ) T + Σ ab , 0 T Σ a , 0 - 1 × [ Σ a , 1 + ( μ a , 1 - μ a , 0 ) ( μ a , 1 - μ a , 0 ) T ] ( Σ a , 0 - 1 ) T Σ ab , 0 - - - ( 19 )
Γ = Σ b , 0 - Σ ab , 0 T Σ a , 0 - 1 Σ ab , 0 - - - ( 20 )
x ^ p , 0 , k | k - 1 = Σ i = 1 N w i χ 0 , k | k - 1 i = μ a , 0 μ b , 0 - - - ( 21 )
P p , 0 , k | k - 1 = Q k + Σ i = 1 N w i ( χ 0 , k | k - 1 i - x ^ p , 0 , k | k - 1 ) ( χ 0 , k | k - 1 i x ^ p , 0 , k | k - 1 ) T = Σ a , 0 Σ ab , 0 ( Σ ab , 0 ) T Σ b , 0 - - - ( 22 )
其中,为所述修正先验概率密度函数的均值,Pp,1,k|k-1为所述修正先验概率密度函数的协方差,a为目标的状态分量。
7.根据权利要求6所述的方法,其特征在于,所述根据所述修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数的步骤具体包括:
根据无迹变换以及修正先验概率密度函数的均值和协方差Pp,1,k|k-1获取N个Sigma点
将所述Sigma点代入非线性观测函数h(·),以获取Sigma点对应的N个观测值具体如下式所示:
z 1 , k | k - 1 i = h ( χ 1 , k | k - 1 i ) - - - ( 23 )
获取所述N个观测值的均值具体如下式所示:
z ^ 1 , k | k - 1 = Σ i = 1 N w i z 1 , k | k - 1 i - - - ( 24 )
根据所述均值获取对应的新息协方差Pzz,1,k|k-1以及交叉协方差Pxz,1,k|k-1,具体如下所示:
P zz , 1 , k | k - 1 = R k + Σ i = 1 N w i ( z 1 k | k - 1 i - z ^ 1 , k | k - 1 ) ( z 1 , k | k - 1 i - z ^ 1 , k | k - 1 ) T - - - ( 25 )
P xz , 1 , k | k - 1 = Σ i = 1 N w i ( χ 1 , k | k - 1 i - x ^ p , 1 , k | k - 1 ) ( z 1 , k | k - 1 i - z ^ 1 , k | k - 1 ) T - - - ( 26 )
根据所述新息协方差Pzz,1,k|k-1、交叉协方差Pxz,1,k|k-1以及k时刻的目标观测向量zk获取k时刻目标状态的第二后验概率密度函数对应的均值和协方差,具体如下所示:
x ^ u , 1 , k | k = x ^ p , 1 , k | k - 1 + P xz , 1 , k | k - 1 P zz , 1 , k | k - 1 - 1 ( z k - z ^ 1 , k | k - 1 ) - - - ( 27 )
P u , 1 , k | k = P p , 1 , k | k - 1 - P xz , 1 , k | k - 1 P zz , 1 , k | k - 1 - 1 P xz , 1 , k | k - 1 T - - - ( 28 )
其中,为所述k时刻目标状态的第二后验概率密度函数对应的均值,Pu,1,k|k为所述k时刻目标状态的第二后验概率密度函数对应的协方差。
8.根据权利要求7所述的方法,其特征在于,所述根据所述原始先验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数的步骤具体包括:
根据所述第一后验概率密度函数以及第二后验概率密度函数各自对应的均值获取目标状态估计权值,具体如下所示:
μ 0 ( x ^ u , 0 , k | k ) = 1 | P zz , 0 , k | k - 1 | · exp ( - ( z k - h ( x ^ u , 0 , k | k ) ) 2 2 ) - - - ( 29 )
μ 1 ( x ^ u , 1 , k | k ) = 1 | P zz , 1 , k | k - 1 | · exp ( - ( z k - h ( x ^ u , 1 , k | k ) ) 2 2 ) - - - ( 30 )
α = μ 1 ( x ^ u , 1 , k | k ) μ 0 ( x ^ u , 0 , k | k ) + μ 1 ( x ^ u , 1 , k | k ) - - - ( 31 )
其中,α为所述目标状态估计权值;
根据所述第一后验概率密度函数、第二后验概率密度函数各自对应的均值、协方差以及目标状态估计权值α获取k时刻目标状态的联合后验概率密度函数,具体如下所示:
p ( x k | z 1 : k ) = N ( x ^ k | k , P k | k ) - - - ( 32 )
x ^ k | k = α k · x ^ u , 1 , k | k + ( 1 - α k ) · x ^ u , 0 , k | k - - - ( 33 )
其中,p(xk|z1:k)为所述k时刻目标状态的联合后验概率密度函数,为所述k时刻目标状态的联合后验概率密度函数的均值,Pk|k为所述k时刻目标状态的联合后验概率密度函数的协方差。
9.一种扩展截断无迹卡尔曼滤波装置,其特征在于,包括:
原始先验概率密度函数获取模块,用于根据无迹变换获取当前目标观测时刻目标状态的原始先验概率密度函数;
第一后验概率密度函数获取模块,用于根据所述原始先验概率密度函数获取当前目标观测时刻目标状态的第一后验概率密度函数;
修正先验概率密度函数获取模块,用于根据统计线性回归理论以及当前目标观测时刻的目标观测向量修正所述原始先验概率密度函数,以获取修正先验概率密度函数;
第二后验概率密度函数获取模块,用于根据所述修正先验概率密度函数获取当前目标观测时刻目标状态的第二后验概率密度函数;
联合后验概率密度函数获取模块,用于根据所述第一后验概率密度函数以及第二后验概率密度函数获取当前目标观测时刻目标状态的联合后验概率密度函数。
CN201410134331.4A 2014-04-03 2014-04-03 一种目标跟踪方法及扩展截断无迹卡尔曼滤波方法、装置 Expired - Fee Related CN103955892B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201410134331.4A CN103955892B (zh) 2014-04-03 2014-04-03 一种目标跟踪方法及扩展截断无迹卡尔曼滤波方法、装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410134331.4A CN103955892B (zh) 2014-04-03 2014-04-03 一种目标跟踪方法及扩展截断无迹卡尔曼滤波方法、装置

Publications (2)

Publication Number Publication Date
CN103955892A true CN103955892A (zh) 2014-07-30
CN103955892B CN103955892B (zh) 2015-10-28

Family

ID=51333161

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410134331.4A Expired - Fee Related CN103955892B (zh) 2014-04-03 2014-04-03 一种目标跟踪方法及扩展截断无迹卡尔曼滤波方法、装置

Country Status (1)

Country Link
CN (1) CN103955892B (zh)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105447574A (zh) * 2015-11-10 2016-03-30 深圳大学 一种辅助截断粒子滤波方法、装置及目标跟踪方法及装置
CN107300697A (zh) * 2017-06-07 2017-10-27 南京航空航天大学 基于无人机的运动目标ukf滤波方法
WO2017185688A1 (zh) * 2016-04-26 2017-11-02 深圳大学 一种在线目标跟踪方法及装置
CN107659989A (zh) * 2017-10-24 2018-02-02 东南大学 无线传感器网络节点分布式测量休眠和目标跟踪方法
CN109117965A (zh) * 2017-06-22 2019-01-01 长城汽车股份有限公司 基于卡尔曼滤波器的系统状态预测装置和方法
CN109903266A (zh) * 2019-01-21 2019-06-18 深圳市华成工业控制有限公司 一种基于样本窗的双核密度估计实时背景建模方法及装置
WO2020024243A1 (zh) * 2018-08-03 2020-02-06 深圳大学 增量核密度估计器的生成方法、装置和计算机可读存储介质
CN111181529A (zh) * 2020-01-17 2020-05-19 中山大学 一种应用于非线性高斯模型的平滑约束扩展卡尔曼滤波方法
CN111954158A (zh) * 2020-07-01 2020-11-17 珠海高凌信息科技股份有限公司 基于rss地图的联合滤波室内单目标跟踪方法、装置及介质
CN113625552A (zh) * 2021-08-16 2021-11-09 西南大学 对状态受限非线性系统进行鲁棒状态估计的方法及装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080231509A1 (en) * 2005-08-18 2008-09-25 Mitsubishi Electric Corporation Gps Positioning Method and Gps Position Device
CN102608595B (zh) * 2012-03-14 2013-06-12 西安电子科技大学 基于分布式相干处理mimo米波雷达的目标定位方法
CN103618326A (zh) * 2013-11-13 2014-03-05 清华大学 基于卡尔曼滤波的风电场中储能系统充放电控制方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080231509A1 (en) * 2005-08-18 2008-09-25 Mitsubishi Electric Corporation Gps Positioning Method and Gps Position Device
CN102608595B (zh) * 2012-03-14 2013-06-12 西安电子科技大学 基于分布式相干处理mimo米波雷达的目标定位方法
CN103618326A (zh) * 2013-11-13 2014-03-05 清华大学 基于卡尔曼滤波的风电场中储能系统充放电控制方法

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105447574A (zh) * 2015-11-10 2016-03-30 深圳大学 一种辅助截断粒子滤波方法、装置及目标跟踪方法及装置
CN105447574B (zh) * 2015-11-10 2018-07-03 深圳大学 一种辅助截断粒子滤波方法、装置及目标跟踪方法及装置
WO2017185688A1 (zh) * 2016-04-26 2017-11-02 深圳大学 一种在线目标跟踪方法及装置
CN107300697A (zh) * 2017-06-07 2017-10-27 南京航空航天大学 基于无人机的运动目标ukf滤波方法
CN109117965A (zh) * 2017-06-22 2019-01-01 长城汽车股份有限公司 基于卡尔曼滤波器的系统状态预测装置和方法
CN107659989B (zh) * 2017-10-24 2020-08-04 东南大学 无线传感器网络节点分布式测量休眠和目标跟踪方法
CN107659989A (zh) * 2017-10-24 2018-02-02 东南大学 无线传感器网络节点分布式测量休眠和目标跟踪方法
WO2020024243A1 (zh) * 2018-08-03 2020-02-06 深圳大学 增量核密度估计器的生成方法、装置和计算机可读存储介质
CN109903266A (zh) * 2019-01-21 2019-06-18 深圳市华成工业控制有限公司 一种基于样本窗的双核密度估计实时背景建模方法及装置
CN109903266B (zh) * 2019-01-21 2022-10-28 深圳市华成工业控制股份有限公司 一种基于样本窗的双核密度估计实时背景建模方法及装置
CN111181529A (zh) * 2020-01-17 2020-05-19 中山大学 一种应用于非线性高斯模型的平滑约束扩展卡尔曼滤波方法
CN111181529B (zh) * 2020-01-17 2022-02-08 中山大学 一种应用于非线性高斯模型的平滑约束扩展卡尔曼滤波方法
CN111954158A (zh) * 2020-07-01 2020-11-17 珠海高凌信息科技股份有限公司 基于rss地图的联合滤波室内单目标跟踪方法、装置及介质
CN113625552A (zh) * 2021-08-16 2021-11-09 西南大学 对状态受限非线性系统进行鲁棒状态估计的方法及装置

Also Published As

Publication number Publication date
CN103955892B (zh) 2015-10-28

Similar Documents

Publication Publication Date Title
CN103955892A (zh) 一种目标跟踪方法及扩展截断无迹卡尔曼滤波方法、装置
CN103902812B (zh) 一种粒子滤波方法、装置及目标跟踪方法、装置
Wang et al. Spherical simplex-radial cubature Kalman filter
Yan et al. State estimation for asynchronous multirate multisensor dynamic systems with missing measurements
Hoteit et al. Particle Kalman filtering: A nonlinear Bayesian framework for ensemble Kalman filters
CN105205313A (zh) 模糊高斯和粒子滤波方法、装置及目标跟踪方法、装置
CN102082560A (zh) 一种基于集合卡尔曼滤波的粒子滤波方法
Andrieu et al. Self-triggered continuous–discrete observer with updated sampling period
CN101819682A (zh) 基于马尔科夫链蒙特卡洛粒子滤波的目标跟踪方法
Sterk et al. Predictability of extreme values in geophysical models
Li et al. Variance change-point detection in panel data models
CN103607772A (zh) 一种基于lmbp神经网络的泰勒定位算法
Bae et al. Multirate moving horizon estimation combined with parameter subset selection
Zjavka Multi-site post-processing of numerical forecasts using a polynomial network substitution for the general differential equation based on operational calculus
Iratni et al. On-line robust nonlinear state estimators for nonlinear bioprocess systems
Ruslan et al. Parameters effect in Sampling Importance Resampling (SIR) particle filter prediction and tracking of flood water level performance
Petrie et al. Ensemble-based data assimilation and the localisation problem
CN103793614A (zh) 突变滤波算法
CN104270119B (zh) 基于非线性未知随机偏差的两阶段容积卡尔曼滤波方法
Berera et al. Information production in homogeneous isotropic turbulence
CN104316905A (zh) 处理飞行时间测距数据的自适应卡尔曼滤波的方法
CN104202019A (zh) 带有未知过程噪声协方差阵递推估计的卡尔曼滤波方法
O'Kane et al. Application of statistical dynamical turbulence closures to data assimilation
Xu et al. Data assimilation with a barotropically unstable shallow water system using representer algorithms
Wen et al. A defensive marginal particle filtering method for data assimilation

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20151028

Termination date: 20170403