CN111985093A - 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法 - Google Patents

一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法 Download PDF

Info

Publication number
CN111985093A
CN111985093A CN202010766995.8A CN202010766995A CN111985093A CN 111985093 A CN111985093 A CN 111985093A CN 202010766995 A CN202010766995 A CN 202010766995A CN 111985093 A CN111985093 A CN 111985093A
Authority
CN
China
Prior art keywords
covariance
noise
state
updating
state estimation
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
CN202010766995.8A
Other languages
English (en)
Other versions
CN111985093B (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 Engineering University
Original Assignee
Harbin Engineering 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 Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN202010766995.8A priority Critical patent/CN111985093B/zh
Publication of CN111985093A publication Critical patent/CN111985093A/zh
Application granted granted Critical
Publication of CN111985093B publication Critical patent/CN111985093B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • 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
    • 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/68—Radar-tracking systems; Analogous systems for angle tracking only
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Filters That Use Time-Delay Elements (AREA)

Abstract

本发明提供一种带噪声估计器的自适应UKF状态估计方法，包括：建立目标跟踪非线性离散状态空间模型并确立UT采样点和权值。将采样点经非线性函数传递，加权处理得到一步预测和预测协方差矩阵。完成预测协方差与观测噪声协方差更新。采用UT方法更新滤波相关参数。联合后验概率密度函数的变分近似的参数解算并完成步骤四中的参数更新。联合后验概率密度函数的变分近似的一步预测协方差与观测噪声协方差解算并更新。输出后验状态估计和估计协方差矩阵。本发明的方法在带有非高斯噪声并且噪声统计特性未知的情况下完成目标跟踪过程中的状态估计任务，其精度和时间消耗低于标准UKF算法，具有良好的实际应用价值。

Description

一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
技术领域
本发明涉及一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，属于涉及被动雷达对目标跟踪过程中噪声未知情况的非线性状态估计领域，特别涉及强非线性、噪声非高斯且未知时变的非理想条件下。具体是指对目标状态在强非线性条件下过程噪声和量测噪声处于非高斯分布且统计特性先验未知的信号处理的方法。
背景技术
动态系统的状态估计问题在目标跟踪领域起着广泛的关注，而针对非线性系统一直是状态估计的难题，滤波器则是处理非线性动态系统状态估计问题的有力手段，因此滤波器的设计则是在非线性动态系统中起着尤为重要的作用。雷达对目标的跟踪过程中，雷达通过自身位置信息通过测得目标角度信息来获取观测目标的位置信息，进一步对目标轨迹进行估计和预测实现跟踪及精确打击。但是由于物质的本质是非线性的，采用经典线性卡尔曼滤波(KF)算法是难以实现精确估计。专家对非线性的状态估计做了大量的研究，其中包括最典型的方法(Extended Kalman Filter,EKF)，此方法是基于非线性函数近似的思想，并要求保证非线性函数连续可微或可导，对于此方法在技术层面要进行雅可比(Jacobian)矩阵的求解。但是在工程应用中，EKF对强非线性系统的处理往往存在很大的截断误差，导致跟踪精度降低，严重时导致滤波发散。另外有学者依靠概率密度函数近似的方法提出解决非线性状态估计问题的方法，例如无迹卡尔曼滤波(Unscented KalmanFilter,UKF),这些非线性变换得到的结果可以满足系统状态分布的后验均值和协方差，完成状态估计，并可以达到三阶泰勒精度逼近任意非线性，效果要优与EKF。但是由于这一类滤波器要求对噪声统计特性先验已知才能对非线性系统精准的估计，而在现实环境中由于受到外界环境干扰以及建模误差，很难采用这一思想的滤波器达到最优状态估计。目前人们在目标跟踪领域内应用的UKF都是基于假设的标准正态分布的高斯噪声，但噪声对滤波器带来的影响也不可忽视，如果对噪声设计一种在线的估计方法并提高滤波算法在不同环境下的自适应能力，方可提升状态估计器的精度。
发明内容
本发明的目的是克服上述滤波算法的局限性，解决现代目标跟踪未知噪声先验信息和建模不精确的问题，降低噪声对系统的干扰，减少动力学模型误差的影响，提高滤波器在非理想环境下的自适应能力和鲁棒性，进一步提高滤波器的精度，因此本发明提出一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法。本发明与标准UKF在强非线性目标跟踪模型非高斯噪声的非理想环境下进行比较，通过蒙特卡洛验证展示出本发明具有更高的状态估计精度和更好的鲁棒性，通过时间消耗分析出本专利可提高状态估计的效率，充分体现出本专利的优越性和有效性，在实际工程实践中具有良好的应用价值。
本发明的目的是这样实现的：步骤如下：
步骤一：建立目标跟踪非线性离散的状态模型和量测模型；
步骤二：根据状态空间模型确立无迹变换采样点和权值；
步骤三：将采样点经非线性函数传递，并进行加权处理得到状态一步预测和状态一步预测协方差矩阵；
步骤四：一步预测协方差与观测噪声协方差更新；
步骤五：采用UT方法更新滤波相关参数；
步骤六：联合后验概率密度函数的变分近似的参数解算并完成步骤四中的参数更新；
步骤七：联合后验概率密度函数的变分近似的一步预测协方差与观测噪声协方差解算并更新；
步骤八：状态滤波更新，输出状态估计和估计协方差矩阵，进入下一次迭代。
本发明还包括这样一些结构特征：
1.步骤一所述的建立描述目标跟踪非线性离散的状态模型和量测模型为：
Figure BDA0002615081120000021
其中:状态方程为xk＝fk-1(xk-1)+wk-1，量测方程为zk＝hk(xk)+vk，k表示第k步，k-1表示第k-1步；xk表示第k步的n维跟踪目标参数状态向量，zk为第k步的m维跟踪目标的量测向量，f(·)和h(·)均为已知的非线性函数，wk-1和vk分别代表k-1步的n维随机系统噪声和第k步m维的量测噪声。
2.步骤二所述的根据状态空间模型确立UT采样点和权值具体过程如下：
Figure BDA0002615081120000022
Figure BDA0002615081120000031
其中，式(2)和式(3)分别为UT采样点ξ和对应的权值ω，
Figure BDA0002615081120000032
和
Figure BDA0002615081120000033
分别表示对应的初始均值权值和初始协方差权值，
Figure BDA0002615081120000034
表示第i步骤所用的权值，n表示系统维数，κ用于确定样本点与均值
Figure BDA0002615081120000035
在k-1时刻之间的距离，Pk-1为k-1时刻的协方差矩阵，α是采样点分布的状态控制，β≥0为待选参数，为非负权系数。
3.步骤三所述将采样点经非线性函数传递，并进行加权处理得到状态一步预测和状态一步预测协方差矩阵的过程为：
Figure BDA0002615081120000036
Figure BDA0002615081120000037
Figure BDA0002615081120000038
其中，f(·)表示非线性函数，∑表示求和运算。
4.步骤四所述一步预测协方差与观测噪声协方差更新的过程为：
Figure BDA0002615081120000039
其中，
Figure BDA00026150811200000310
Figure BDA00026150811200000311
τ和ρ称为协调系数和遗忘因子，Pk|k-1表示在k-1时刻的预测。
5.步骤五所述采用UT方法更新滤波相关参数的过程为：
Figure BDA00026150811200000312
Figure BDA00026150811200000313
其中，等式左侧上标表示迭代次数，下标表示运行时刻，等式右侧N(·；μ,P)表示为均值为μ方差为P的高斯分布。
6.步骤六所述的联合后验概率密度函数的变分近似的参数解算并完成步骤四中的参数更新的过程为：
Figure BDA0002615081120000041
Figure BDA0002615081120000042
其中，G(·)表示伽马分布函数，q(·)表示p(·)近似概率密度函数，通过近似后验概率密度函数q(·)进行最小化Kullback-Leibler散度，q(·)和p(·)之间的Kullback-Leibler散度表示为：
Figure BDA0002615081120000043
将一步预测协方差和量测噪声协方差建模为Inverse-Wishart分布：
Figure BDA0002615081120000044
Figure BDA0002615081120000045
其中，IW(·；a,b)表示Inverse-Wishart分布，a和b分别表示自由度参数和逆尺度矩阵；通过构建联合后验概率密度函数：
Figure BDA0002615081120000046
其中式(15)中的最优解满足：
Figure BDA0002615081120000047
其中，E(·)表示求期望运算，log(·)表示自然对数运算，θ是Θ中的任意元素，
Figure BDA0002615081120000048
表示Θ中的非θ元素，cθ表示随机变量θ的常数；联合概率密度表示为：
Figure BDA0002615081120000051
令θ＝ξk则有：
Figure BDA0002615081120000052
Figure BDA0002615081120000053
令θ＝λk则有：
Figure BDA0002615081120000054
Figure BDA0002615081120000055
令θ＝Pk|k-1则有：
Figure BDA0002615081120000056
Figure BDA0002615081120000057
令θ＝Rk则有：
Figure BDA0002615081120000061
Figure BDA0002615081120000062
7.步骤七所述联合后验概率密度函数的变分近似的一步预测协方差与观测噪声协方差解算并更新的过程为：
Figure BDA0002615081120000063
Figure BDA0002615081120000064
8.步骤八所述状态滤波更新，输出状态估计和估计协方差矩阵，进入下一次迭代的过程为：
令θ＝xk则有：
Figure BDA0002615081120000065
Figure BDA0002615081120000066
Figure BDA0002615081120000067
其中，采用UT方法计算
Figure BDA0002615081120000068
的表达式如下所示：
Figure BDA0002615081120000069
与现有技术相比，本发明的有益效果是：1.在目标跟踪中，标准UKF在非高斯噪声的影响下导致精度变差，本专利采用变分贝推理的方法联合估计状态一步预测和噪声协方差矩阵，从而在非高斯噪声的情况下提高目标跟踪。2.在现有的技术中，计算时间都要高于标准UKF算法，这是以时间为代价换取的精度，本专利的优点在非高斯状态估计过程中，单步运行时间低于标准UKF而精度高于标准UKF算法，提高了算法的效率。
附图说明
图1是本发明提供的方法与基于二阶UKF的目标跟踪应用中对目标x坐标轴方向位置估计的均方误差曲线；图中A表示本发明提供的方法对目标x坐标轴方向位置估计的均方误差曲线；图中B表示提供的二阶UKF方法对目标x坐标轴方向位置估计的均方误差曲线。
图2是本发明提供的方法与基于二阶UKF的目标跟踪应用中对目标y坐标轴方向位置估计的均方误差曲线。
图3是本发明提供的方法与基于二阶UKF的目标跟踪应用中单步时间消耗曲线。
图4是本发明的算法流程图。
具体实施方式
下面结合附图与具体实施方式对本发明作进一步详细描述。
结合图1至图4，
实施方式一、本发明的步骤如下：
步骤一：建立目标跟踪非线性离散的状态模型和量测模型；
步骤二：根据状态空间模型确立UT采样点和权值；
步骤三：将采样点经非线性函数传递，并进行加权处理得到状态一步预测和状态一步预测协方差矩阵；
步骤四：一步预测协方差与观测噪声协方差更新；
步骤五：采用UT方法更新滤波相关参数；
步骤六：联合后验概率密度函数的变分近似的参数解算并完成步骤四中的参数更新；
步骤七：联合后验概率密度函数的变分近似的一步预测协方差与观测噪声协方差解算并更新；
步骤八：状态滤波更新，输出状态估计和估计协方差矩阵，进入下一次迭代。
实施方式二、本实施方法与具体实施方法一所述的基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法的区别在于步骤一所建立目标跟踪非线性离散的状态模型和量测模型为：
Figure BDA0002615081120000071
其中状态方程为xk＝fk-1(xk-1)+wk-1，量测方程为zk＝hk(xk)+vk，k表示第k步，k-1表示第k-1步。xk表示第k步的n维跟踪目标参数状态向量，zk为第k步的m维跟踪目标的量测向量，f(·)和h(·)均为已知的非线性函数，wk-1和vk分别代表k-1步的n维随机系统噪声和第k步m维的量测噪声。满足wk-1～N(0,Qk-1)和vk～N(0,Rk-1)，N(·)表示正态分布，并且随机系统噪声服从均值为0方差为Qk-1的高斯分布，Qk-1表示第k-1步系统噪声的方差阵。随机量测噪声向量服从均值为0方差为Rk的高斯分布，Rk表示第k步量测噪声的方差阵，且满足wk-1与vk不相关。
实施方式三、结合图4具体说明本专利实施方式，本实施方式与具体实施方式一或二所述的基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法的区别在于，步骤二所述的根据状态空间模型确立UT采样点和权值：
Figure BDA0002615081120000081
Figure BDA0002615081120000082
其中，式(2)和式(3)分别为UT采样点ξ和对应的权值ω，
Figure BDA0002615081120000083
和
Figure BDA0002615081120000084
分别表示对应的初始均值权值和初始协方差权值，
Figure BDA0002615081120000085
表示第i步骤所用的权值，n表示系统维数，κ用于确定样本点与均值
Figure BDA0002615081120000086
在k-1时刻之间的距离，Pk-1为k-1时刻的协方差矩阵，α是采样点分布的状态控制，β≥0为待选参数，为非负权系数。
实施方式四、本实施方式与具体实施方式三所述的基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法的区别在于，步骤三所述的将采样点经非线性函数传递，并进行加权处理得到状态一步预测和状态一步预测协方差矩阵的过程为：
Figure BDA0002615081120000087
Figure BDA0002615081120000091
Figure BDA0002615081120000092
其中，f(·)表示非线性函数，∑表示求和运算。
实施方式五、本实施方式与具体实施方式四所述的基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法的区别在于，步骤四所述的将采样点经非线性函数传递，并进行加权处理得到一步预测协方差与观测噪声协方差更新的过程为：
Figure BDA0002615081120000093
其中，
Figure BDA0002615081120000094
Figure BDA0002615081120000095
τ和ρ称为协调系数和遗忘因子，Pk|k-1表示在k-1时刻的预测。
实施方式六、本实施方式与具体实施方式五所述的基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法的区别在于，步骤五所述的采用UT方法更新滤波相关参数的过程为：
Figure BDA0002615081120000096
Figure BDA0002615081120000097
其中，等式左侧上标表示迭代次数，下标表示运行时刻，等式右侧N(·；μ,P)表示为均值为μ方差为P的高斯分布。
实施方式七、本实施方式与具体实施方式六所述的基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法的区别在于，步骤六所述的联合后验概率密度函数的变分近似的参数解算并完成步骤四中的参数更新的过程为：
Figure BDA0002615081120000098
Figure BDA0002615081120000099
其中，G(·)表示伽马分布函数，q(·)表示p(·)近似概率密度函数，通过近似后验概率密度函数q(·)进行最小化Kullback-Leibler散度，q(·)和p(·)之间的Kullback-Leibler散度表示为：
Figure BDA0002615081120000101
将一步预测协方差和量测噪声协方差建模为Inverse-Wishart分布：
Figure BDA0002615081120000102
Figure BDA0002615081120000103
其中，IW(·；a,b)表示Inverse-Wishart分布，a和b分别表示自由度参数和逆尺度矩阵。通过构建联合后验概率密度函数：
Figure BDA0002615081120000104
其中式(15)中的最优解满足：
Figure BDA0002615081120000105
其中，E(·)表示求期望运算，log(·)表示自然对数运算，θ是Θ中的任意元素，
Figure BDA0002615081120000106
表示Θ中的非θ元素，cθ表示随机变量θ的常数。联合概率密度表示为：
Figure BDA0002615081120000107
令θ＝ξk则有：
Figure BDA0002615081120000108
Figure BDA0002615081120000111
令θ＝λk则有：
Figure BDA0002615081120000112
Figure BDA0002615081120000113
令θ＝Pk|k-1则有：
Figure BDA0002615081120000114
Figure BDA0002615081120000115
令θ＝Rk则有：
Figure BDA0002615081120000116
Figure BDA0002615081120000117
实施方式八、基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法的区别在于，步骤七联合后验概率密度函数的变分近似的一步预测协方差与观测噪声协方差解算并更新的过程为：
Figure BDA0002615081120000121
Figure BDA0002615081120000122
实施方式九、基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法的区别在于，步骤八状态滤波更新，输出状态估计和估计协方差矩阵，进入下一次迭代的过程为：令θ＝xk则有：
Figure BDA0002615081120000123
Figure BDA0002615081120000124
Figure BDA0002615081120000125
其中，采用UT方法计算
Figure BDA0002615081120000126
的表达式如下所示：
Figure BDA0002615081120000127
在目标行进过程中数据信息不断的更新，雷达站能够获得目标方位信息，但此方位信息含有非高斯噪声分量，具体表现在过程非高斯和量测非高斯，但二者处于不同的自由度。通常应用EKF，UKF，CKF等方法来获得目标的状态信息，但由于受到非高斯噪声的影响，对于对目标定位跟踪精度要求较高的场合,上述滤波方法并不能满足要求。本发明提供的方法比现有的方法具备更高的估计精度和鲁棒性，下面以具体实例来说明本发明的优越性。具体如下：
依据纯方位目标跟踪问题，建立如下描述目标跟踪系统的状态方程和量测方程：
Figure BDA0002615081120000131
其中，第k步的状态向量
Figure BDA0002615081120000132
表示笛卡尔坐标系内目标的位置，zk表示观测向量，wk和vk表示系统噪声和量测噪声且对应的噪声方差为Qk和
Figure BDA0002615081120000133
矩阵
Figure BDA0002615081120000134
表征了目标x方向的速度为0.9m/s，y方向的速度为1m/s；随机系统噪声wk方差为
Figure BDA0002615081120000135
随机量测噪声vk的方差为Rk＝0.025rad2；初始真实状态设定值x0＝[20m 5m]，初始协方差设定为P0＝[0.1m2 0m2；0m20.1m2]表征了目标初始位置的不确定性。针对非高斯噪声本文设定wk和vk表达如下：
Figure BDA0002615081120000136
Figure BDA0002615081120000137
其中，w.p.表示百分比。
实施过程：仿真过程中根据式(35)均方误差MSE做为衡量滤波器的性能指标：
Figure BDA0002615081120000138
其中，N表示Monte Carlo次数。对目标位置估计的均方误差越小代表跟踪精度越高，效果越好。仿真时间为100s，在相同的仿真条件下，采用二阶UKF和本发明提出的带噪声估计器的自适应UKF进行对比，所有方法均进行500次Monte Carlo仿真。
实施效果：图1给出了纯方位目标跟踪过程中，本发明提供的方法与二阶UKF、对x方向位置估计的均方误差曲线，图2给出了纯方位目标跟踪过程中，本发明提供的方法与二阶UKF、对y方向位置估计的均方误差曲线。在图1和图2中,曲线A代表本发明提供方法状态估计值的均方误差，曲线B代表二阶UKF提供方法状态估计值的均方误差。状态估计值的均方误差越小代表估计精度越高，性能越好。从图1和图2可以看出，本专利提出新型带噪声估计器的自适应UKF状态估计方法优于二阶UKF。图3给出单步运行的时间消耗，不难看出本专利提出新型带噪声估计器的自适应UKF状态估计方法时间要小于二阶UKF，提高状态估计的效率，节省运行成本。从以上实施例不难看出，相对于二阶UKF，本发明提供方法能够在消耗时间更少的同时，获得更高的精度以及获取对目标更准确的跟踪。
综上，本发明设计的是基于一种新型带噪声估计器的自适应无迹卡尔曼滤波状态估计方法。涉及被动雷达对目标跟踪过程中噪声未知情况的非线性状态估计领域，特别涉及强非线性、噪声非高斯且未知时变的非理想条件下。包括1、建立目标跟踪非线性离散的状态模型和量测模型。2、根据状态空间模型确立UT采样点和权值。3、将采样点经非线性函数传递，并进行加权处理得到状态一步预测和状态一步预测协方差矩阵。4、一步预测协方差与观测噪声协方差更新。5、采用UT方法更新滤波相关参数。6、联合后验概率密度函数的变分近似的参数解算并完成步骤四中的参数更新。7、联合后验概率密度函数的变分近似的一步预测协方差与观测噪声协方差解算并更新。8、状态滤波更新，输出状态估计和估计协方差矩阵。本发明的方法在带有非高斯噪声并且噪声统计特性未知的情况下完成目标跟踪过程中的状态估计任务，其精度高于标准UKF但时间低于标准UKF算法，具有良好的实际应用价值。

Claims (9)

1.一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤如下：
步骤一：建立目标跟踪非线性离散的状态模型和量测模型；
步骤二：根据状态空间模型确立无迹变换采样点和权值；
步骤三：将采样点经非线性函数传递，并进行加权处理得到状态一步预测和状态一步预测协方差矩阵；
步骤四：一步预测协方差与观测噪声协方差更新；
步骤五：采用UT方法更新滤波相关参数；
步骤六：联合后验概率密度函数的变分近似的参数解算并完成步骤四中的参数更新；
步骤七：联合后验概率密度函数的变分近似的一步预测协方差与观测噪声协方差解算并更新；
步骤八：状态滤波更新，输出状态估计和估计协方差矩阵，进入下一次迭代。
2.根据权利要求1所述的一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤一所述的建立描述目标跟踪非线性离散的状态模型和量测模型为：
Figure FDA0002615081110000011
其中:状态方程为xk＝fk-1(xk-1)+wk-1，量测方程为zk＝hk(xk)+vk，k表示第k步，k-1表示第k-1步；xk表示第k步的n维跟踪目标参数状态向量，zk为第k步的m维跟踪目标的量测向量，f(·)和h(·)均为已知的非线性函数，wk-1和vk分别代表k-1步的n维随机系统噪声和第k步m维的量测噪声。
3.根据权利要求1所述的一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤二所述的根据状态空间模型确立UT采样点和权值具体过程如下：
Figure FDA0002615081110000012
Figure FDA0002615081110000021
其中，式(2)和式(3)分别为UT采样点ξ和对应的权值ω，
Figure FDA0002615081110000022
和
Figure FDA0002615081110000023
分别表示对应的初始均值权值和初始协方差权值，
Figure FDA0002615081110000024
表示第i步骤所用的权值，n表示系统维数，κ用于确定样本点与均值
Figure FDA0002615081110000025
在k-1时刻之间的距离，Pk-1为k-1时刻的协方差矩阵，α是采样点分布的状态控制，β≥0为待选参数，为非负权系数。
4.根据权利要求1所述的一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤三所述将采样点经非线性函数传递，并进行加权处理得到状态一步预测和状态一步预测协方差矩阵的过程为：
Figure FDA0002615081110000026
Figure FDA0002615081110000027
Figure FDA0002615081110000028
其中，f(·)表示非线性函数，∑表示求和运算。
5.根据权利要求1所述的一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤四所述一步预测协方差与观测噪声协方差更新的过程为：
Figure FDA0002615081110000029
其中，
Figure FDA00026150811100000210
Figure FDA00026150811100000211
τ和ρ称为协调系数和遗忘因子，Pk|k-1表示在k-1时刻的预测。
6.根据权利要求1所述的一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤五所述采用UT方法更新滤波相关参数的过程为：
Figure FDA00026150811100000212
Figure FDA0002615081110000031
其中，等式左侧上标表示迭代次数，下标表示运行时刻，等式右侧N(·；μ,P)表示为均值为μ方差为P的高斯分布。
7.根据权利要求1所述的一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤六所述的联合后验概率密度函数的变分近似的参数解算并完成步骤四中的参数更新的过程为：
Figure FDA0002615081110000032
Figure FDA0002615081110000033
其中，G(·)表示伽马分布函数，q(·)表示p(·)近似概率密度函数，通过近似后验概率密度函数q(·)进行最小化Kullback-Leibler散度，q(·)和p(·)之间的Kullback-Leibler散度表示为：
Figure FDA0002615081110000034
将一步预测协方差和量测噪声协方差建模为Inverse-Wishart分布：
Figure FDA0002615081110000035
Figure FDA0002615081110000036
其中，IW(·；a,b)表示Inverse-Wishart分布，a和b分别表示自由度参数和逆尺度矩阵；通过构建联合后验概率密度函数：
Figure FDA0002615081110000037
其中式(15)中的最优解满足：
Figure FDA0002615081110000038
其中，E(·)表示求期望运算，log(·)表示自然对数运算，θ是Θ中的任意元素，
Figure FDA0002615081110000039
表示Θ中的非θ元素，cθ表示随机变量θ的常数；联合概率密度表示为：
Figure FDA0002615081110000041
令θ＝ξk则有：
Figure FDA0002615081110000042
Figure FDA0002615081110000043
令θ＝λk则有：
Figure FDA0002615081110000044
Figure FDA0002615081110000045
令θ＝Pk|k-1则有：
Figure FDA0002615081110000046
Figure FDA0002615081110000047
令θ＝Rk则有：
Figure FDA0002615081110000051
Figure FDA0002615081110000052
8.根据权利要求1所述的一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤七所述联合后验概率密度函数的变分近似的一步预测协方差与观测噪声协方差解算并更新的过程为：
Figure FDA0002615081110000053
Figure FDA0002615081110000054
9.根据权利要求1所述的一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法，其特征在于：步骤八所述状态滤波更新，输出状态估计和估计协方差矩阵，进入下一次迭代的过程为：
令θ＝xk则有：
Figure FDA0002615081110000055
Figure FDA0002615081110000056
Figure FDA0002615081110000057
其中，采用UT方法计算
Figure FDA0002615081110000058
的表达式如下所示：
Figure FDA0002615081110000061
CN202010766995.8A 2020-08-03 2020-08-03 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法 Active CN111985093B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010766995.8A CN111985093B (zh) 2020-08-03 2020-08-03 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010766995.8A CN111985093B (zh) 2020-08-03 2020-08-03 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法

Publications (2)

Publication Number Publication Date
CN111985093A true CN111985093A (zh) 2020-11-24
CN111985093B CN111985093B (zh) 2022-06-21

Family

ID=73445031

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010766995.8A Active CN111985093B (zh) 2020-08-03 2020-08-03 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法

Country Status (1)

Country Link
CN (1) CN111985093B (zh)

Cited By (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112468116A (zh) * 2020-12-01 2021-03-09 哈尔滨工程大学 一种基于Gibbs采样器的自适应平滑方法
CN112836354A (zh) * 2021-01-12 2021-05-25 中南大学 一种目标跟踪定位方法、系统、装置及可读存储介质
CN113011475A (zh) * 2021-01-29 2021-06-22 深圳信息职业技术学院 考虑相关噪声和随机参数矩阵的分布式融合算法
CN113407909A (zh) * 2021-07-15 2021-09-17 东南大学 一种用于非解析复数非线性系统的无味算法
CN113419278A (zh) * 2021-06-21 2021-09-21 大庆油田有限责任公司 一种基于状态空间模型与支持向量回归的井震联合多目标同时反演方法
CN113452349A (zh) * 2021-06-28 2021-09-28 中山大学 一种基于贝叶斯序贯重要性积分的卡尔曼滤波方法
CN113591020A (zh) * 2021-07-23 2021-11-02 江南大学 一种基于轴对称盒空间滤波的非线性系统状态估计方法
CN113659962A (zh) * 2021-08-03 2021-11-16 青岛迈金智能科技有限公司 一种盘爪踏频计及用于盘爪踏频计的参数优化方法
CN113670315A (zh) * 2021-08-25 2021-11-19 江南大学 一种基于变分迭代卡尔曼滤波的李群重尾干扰噪声动态飞行器姿态估计方法
CN113759309A (zh) * 2021-08-31 2021-12-07 河海大学 一种室内移动目标定位方法、装置和计算机设备
CN113779497A (zh) * 2021-10-13 2021-12-10 东南大学 一种解决量测信息存在随机时延和丢包的目标跟踪方法
CN113821952A (zh) * 2021-09-18 2021-12-21 国网浙江省电力有限公司舟山供电公司 一种基于卡尔曼滤波算法的数字孪生无砟轨道优化方法
CN113971752A (zh) * 2021-09-18 2022-01-25 江苏大学 一种抗观测数据干扰的多车协同状态估计方法
CN114609912A (zh) * 2022-03-18 2022-06-10 电子科技大学 基于伪线性最大相关熵卡尔曼滤波的仅测角目标追踪方法
CN114614797A (zh) * 2022-05-12 2022-06-10 之江实验室 基于广义最大非对称相关熵准则的自适应滤波方法和系统
CN114897111A (zh) * 2022-07-15 2022-08-12 深圳市协和传动器材有限公司 一种用于监控间歇式凸轮分割器运行状况的方法
CN115201799A (zh) * 2022-09-09 2022-10-18 昆明理工大学 一种用于声呐的时变卡尔曼滤波跟踪方法
CN116595897A (zh) * 2023-07-17 2023-08-15 广东工业大学 一种基于消息传递的非线性动态系统状态估计方法和装置
CN116884516A (zh) * 2023-09-05 2023-10-13 河南科技学院 一种基于svm-ukf数据融合的pm2.5浓度预测方法
CN117908034A (zh) * 2024-03-20 2024-04-19 西北工业大学 基于自适应波束跟踪水下目标高稳健模基doa估计方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015033503A1 (ja) * 2013-09-05 2015-03-12 カルソニックカンセイ株式会社 推定装置及び推定方法
CN109459033A (zh) * 2018-12-21 2019-03-12 哈尔滨工程大学 一种多重渐消因子的机器人无迹快速同步定位与建图方法
CN110455287A (zh) * 2019-07-24 2019-11-15 南京理工大学 自适应无迹卡尔曼粒子滤波方法
CN111047627A (zh) * 2019-11-14 2020-04-21 中山大学 一种平滑约束无迹卡尔曼滤波方法及目标跟踪方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015033503A1 (ja) * 2013-09-05 2015-03-12 カルソニックカンセイ株式会社 推定装置及び推定方法
CN109459033A (zh) * 2018-12-21 2019-03-12 哈尔滨工程大学 一种多重渐消因子的机器人无迹快速同步定位与建图方法
CN110455287A (zh) * 2019-07-24 2019-11-15 南京理工大学 自适应无迹卡尔曼粒子滤波方法
CN111047627A (zh) * 2019-11-14 2020-04-21 中山大学 一种平滑约束无迹卡尔曼滤波方法及目标跟踪方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
F. DENG 等: "Adaptive unscented Kalman filter for parameter and state estimation of nonlinear high-speed objects", 《JOURNAL OF SYSTEMS ENGINEERING AND ELECTRONICS》 *
赵洪山 等: "基于自适应无迹卡尔曼滤波的电力系统动态状态估计", 《电网技术》 *

Cited By (32)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112468116A (zh) * 2020-12-01 2021-03-09 哈尔滨工程大学 一种基于Gibbs采样器的自适应平滑方法
CN112468116B (zh) * 2020-12-01 2023-06-16 哈尔滨工程大学 一种基于Gibbs采样器的自适应平滑方法
CN112836354A (zh) * 2021-01-12 2021-05-25 中南大学 一种目标跟踪定位方法、系统、装置及可读存储介质
CN113011475A (zh) * 2021-01-29 2021-06-22 深圳信息职业技术学院 考虑相关噪声和随机参数矩阵的分布式融合算法
CN113011475B (zh) * 2021-01-29 2022-12-02 深圳信息职业技术学院 考虑相关噪声和随机参数矩阵的分布式融合方法
CN113419278A (zh) * 2021-06-21 2021-09-21 大庆油田有限责任公司 一种基于状态空间模型与支持向量回归的井震联合多目标同时反演方法
CN113419278B (zh) * 2021-06-21 2022-03-08 大庆油田有限责任公司 一种基于状态空间模型与支持向量回归的井震联合多目标同时反演方法
CN113452349A (zh) * 2021-06-28 2021-09-28 中山大学 一种基于贝叶斯序贯重要性积分的卡尔曼滤波方法
CN113407909B (zh) * 2021-07-15 2024-01-09 东南大学 一种用于非解析复数非线性系统的无味算法的计算方法
CN113407909A (zh) * 2021-07-15 2021-09-17 东南大学 一种用于非解析复数非线性系统的无味算法
CN113591020B (zh) * 2021-07-23 2024-03-01 江南大学 一种基于轴对称盒空间滤波的非线性系统状态估计方法
CN113591020A (zh) * 2021-07-23 2021-11-02 江南大学 一种基于轴对称盒空间滤波的非线性系统状态估计方法
CN113659962A (zh) * 2021-08-03 2021-11-16 青岛迈金智能科技有限公司 一种盘爪踏频计及用于盘爪踏频计的参数优化方法
CN113670315B (zh) * 2021-08-25 2024-03-15 无锡北微传感科技有限公司 一种基于变分迭代卡尔曼滤波的李群重尾干扰噪声动态飞行器姿态估计方法
CN113670315A (zh) * 2021-08-25 2021-11-19 江南大学 一种基于变分迭代卡尔曼滤波的李群重尾干扰噪声动态飞行器姿态估计方法
CN113759309A (zh) * 2021-08-31 2021-12-07 河海大学 一种室内移动目标定位方法、装置和计算机设备
CN113971752A (zh) * 2021-09-18 2022-01-25 江苏大学 一种抗观测数据干扰的多车协同状态估计方法
CN113971752B (zh) * 2021-09-18 2024-03-19 江苏大学 一种抗观测数据干扰的多车协同状态估计方法
CN113821952A (zh) * 2021-09-18 2021-12-21 国网浙江省电力有限公司舟山供电公司 一种基于卡尔曼滤波算法的数字孪生无砟轨道优化方法
CN113821952B (zh) * 2021-09-18 2024-02-20 国网浙江省电力有限公司舟山供电公司 一种基于卡尔曼滤波算法的数字孪生无砟轨道优化方法
CN113779497B (zh) * 2021-10-13 2022-11-18 东南大学 一种解决量测信息存在随机时延和丢包的目标跟踪方法
CN113779497A (zh) * 2021-10-13 2021-12-10 东南大学 一种解决量测信息存在随机时延和丢包的目标跟踪方法
CN114609912A (zh) * 2022-03-18 2022-06-10 电子科技大学 基于伪线性最大相关熵卡尔曼滤波的仅测角目标追踪方法
CN114609912B (zh) * 2022-03-18 2023-10-03 电子科技大学 基于伪线性最大相关熵卡尔曼滤波的仅测角目标追踪方法
CN114614797B (zh) * 2022-05-12 2022-09-30 之江实验室 基于广义最大非对称相关熵准则的自适应滤波方法和系统
CN114614797A (zh) * 2022-05-12 2022-06-10 之江实验室 基于广义最大非对称相关熵准则的自适应滤波方法和系统
CN114897111A (zh) * 2022-07-15 2022-08-12 深圳市协和传动器材有限公司 一种用于监控间歇式凸轮分割器运行状况的方法
CN115201799A (zh) * 2022-09-09 2022-10-18 昆明理工大学 一种用于声呐的时变卡尔曼滤波跟踪方法
CN116595897B (zh) * 2023-07-17 2023-10-10 广东工业大学 一种基于消息传递的非线性动态系统状态估计方法和装置
CN116595897A (zh) * 2023-07-17 2023-08-15 广东工业大学 一种基于消息传递的非线性动态系统状态估计方法和装置
CN116884516A (zh) * 2023-09-05 2023-10-13 河南科技学院 一种基于svm-ukf数据融合的pm2.5浓度预测方法
CN117908034A (zh) * 2024-03-20 2024-04-19 西北工业大学 基于自适应波束跟踪水下目标高稳健模基doa估计方法

Also Published As

Publication number Publication date
CN111985093B (zh) 2022-06-21

Similar Documents

Publication Publication Date Title
CN111985093B (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
CN111178385B (zh) 一种鲁棒在线多传感器融合的目标跟踪方法
CN110503071B (zh) 基于变分贝叶斯标签多伯努利叠加模型的多目标跟踪方法
CN108896986B (zh) 一种基于预测值的量测转换序贯滤波机动目标跟踪方法
CN108319570B (zh) 一种异步多传感器空时偏差联合估计与补偿方法及装置
CN104199022B (zh) 一种基于目标模态估计的临近空间高超声速目标跟踪方法
CN111693984B (zh) 一种改进的ekf-ukf动目标跟踪方法
CN110231620B (zh) 一种噪声相关系统跟踪滤波方法
CN111125885A (zh) 一种基于改进克里金插值算法的asf修正表构建方法
CN108871365B (zh) 一种航向约束下的状态估计方法及系统
CN111913175A (zh) 一种传感器短暂失效下带补偿机制的水面目标跟踪方法
CN107290742A (zh) 一种非线性目标跟踪系统中平方根容积卡尔曼滤波方法
Bai et al. A Robust Generalized $ t $ Distribution-Based Kalman Filter
CN111711432A (zh) 一种基于ukf和pf混合滤波的目标跟踪算法
CN115204212A (zh) 一种基于stm-pmbm滤波算法的多目标跟踪方法
CN114236480A (zh) 一种机载平台传感器系统误差配准算法
CN110989341B (zh) 一种约束辅助粒子滤波方法及目标跟踪方法
CN114445459B (zh) 基于变分贝叶斯理论的连续-离散最大相关熵目标跟踪方法
CN105549003A (zh) 一种汽车雷达目标跟踪方法
CN115685128A (zh) 一种机动目标场景下的雷达目标跟踪算法及电子设备
CN115469314A (zh) 一种均匀圆环阵稳健水下目标方位跟踪方法及系统
CN114611068A (zh) 一种高机动目标跟踪方法
CN112241583A (zh) 一种最小化后验距离的传感器路径优化方法
Singh et al. Cubature and quadrature based continuous-discrete filters for maneuvering target tracking
Jian et al. Research on Performance of Terrain Matching Algorithm for Underwater Autonomous Vehicle Based on Particle Filter

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