CN111310110A - 一种高维耦合不确定***混合状态估计方法 - Google Patents

一种高维耦合不确定***混合状态估计方法 Download PDF

Info

Publication number
CN111310110A
CN111310110A CN202010205361.5A CN202010205361A CN111310110A CN 111310110 A CN111310110 A CN 111310110A CN 202010205361 A CN202010205361 A CN 202010205361A CN 111310110 A CN111310110 A CN 111310110A
Authority
CN
China
Prior art keywords
state
model
representing
parameter
dimensional
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
CN202010205361.5A
Other languages
English (en)
Other versions
CN111310110B (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.)
Zhengzhou University of Light Industry
Original Assignee
Zhengzhou University of Light Industry
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 Zhengzhou University of Light Industry filed Critical Zhengzhou University of Light Industry
Priority to CN202010205361.5A priority Critical patent/CN111310110B/zh
Publication of CN111310110A publication Critical patent/CN111310110A/zh
Application granted granted Critical
Publication of CN111310110B publication Critical patent/CN111310110B/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
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/02—Total factory control, e.g. smart factories, flexible manufacturing systems [FMS] or integrated manufacturing systems [IMS]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提出了一种高维耦合不确定***混合状态估计方法，其步骤为：首先，构建高维深耦合不确定***状态、参数和测量的模型，并设计对应的观测器模型以得到状态和参数的估计值；其次，对***离散化将其分解为低维离散化混合模型，进而得到低维离散化混合模型和参数模型；最后，以观测器输出估计值为辅助信号，利用容积卡尔曼滤波算法分别对低维离散化混合模型的状态估计值进行滤波处理，输出低维离散化混合模型的状态值。本发明通过观测器输出的估计值修正***模型，能够在保障***稳定性的前提下，有效提高***状态估计精度，同时低维容积卡尔曼滤波算法降低滤波计算过程的计算维数，适用于具有不确定性的高维、耦合、非线性***状态估计。

Description

一种高维耦合不确定***混合状态估计方法
技术领域
本发明涉及混合状态估计技术领域，特别是指一种高维耦合不确定***混合状态估计方法。
背景技术
针对线性精确模型，线性卡尔曼滤波算法提供了***的最优递推解。但对于实际复杂工程***而言，***模型通常是非线性、高维、耦合的，且***模型存在参数不确定性。此时，经典的基于模型的估计方法会出现计算量过大，估计精度降低，甚至发散的问题。目前还没有一种有效的针对高维、非线性、不确定***的状态估计方法。因而开展高维、非线性、不确定***模型下的状态估计方法具有一定的前沿性和实际应用价值。
发明内容
针对现有的状态估计方法存在计算量大、估计精度低且结果发散的技术问题，本发明提出了一种高维耦合不确定***混合状态估计方法，能够实现高维、非线性、不确定非线性***的精确、有效、稳定的状态和参数估计，降低估计过程计算维数，并且在确保估计过程稳定性的前提下，提高***状态估计的精度。
本发明的技术方案是这样实现的：
一种高维耦合不确定***混合状态估计方法，其步骤如下：
S1、构建高维深耦合不确定***状态模型、参数模型和测量模型；
S2、根据***状态模型、参数模型和测量模型分别构建状态和参数的观测器模型，并利用状态和参数的观测器模型估计***状态和参数的估计值；
S3、利用欧拉法对步骤S1中的***状态模型、参数模型和测量模型进行离散化处理，得到低维离散化混合模型，并对低维离散化混合模型中的不确定参数进行建模，得到参数离散模型；
S4、基于步骤S2中的***状态和参数的估计值，利用容积卡尔曼滤波算法分别对步骤S3中的低维离散化混合模型和参数离散模型进行滤波处理，输出低维离散化混合模型的状态值。
所述步骤S1中的高维深耦合不确定***状态模型、参数模型和测量模型分别为：
Figure BDA0002419320900000011
Figure BDA0002419320900000012
y(t)＝h(x(t),t)+υ(t) (3)，
其中，
Figure BDA0002419320900000021
表示***状态的一阶导数，f(·)为非线性状态转移矩阵，x(t)表示随时间t变化的状态变量，θ(t)表示不确定参数，u(t)表示***输入变量，ω1(t)为状态过程噪声，
Figure BDA0002419320900000022
表示***不确定参数的一阶导数，ω2(t)为参数过程噪声，y(t)为测量值，h(·)为测量转移矩阵，υ(t)为测量噪声。
所述根据***状态模型、参数模型和测量模型分别构建状态和参数的观测器模型分别为：
Figure BDA0002419320900000023
Figure BDA0002419320900000024
其中，
Figure BDA0002419320900000025
表示***状态估计值的一阶导数，
Figure BDA0002419320900000026
表示线性状态估计值，
Figure BDA0002419320900000027
表示参数估计值的一阶导数，
Figure BDA0002419320900000028
表示输入变量估计值，K1和K2表示对应的观测器增益矩阵，
Figure BDA0002419320900000029
表示***实际测量值与观测器的输出测量值之间的偏差，
Figure BDA00024193209000000210
表示观测器的输出测量值。
所述低维离散化混合模型为：
Figure BDA00024193209000000211
Figure BDA00024193209000000212
Figure BDA00024193209000000213
其中，
Figure BDA00024193209000000214
表示k时刻的非线性状态分量，
Figure BDA00024193209000000215
表示k-1时刻的非线性状态分量，
Figure BDA00024193209000000216
表示k时刻的线性状态分量，
Figure BDA00024193209000000217
表示k-1时刻的线性状态分量，
Figure BDA00024193209000000218
表示k-1时刻的非线性状态不确定参数，
Figure BDA00024193209000000219
表示k-1时刻的线性状态不确定参数，
Figure BDA00024193209000000220
表示k时刻与非线性***状态相关的输入变量，
Figure BDA00024193209000000221
表示k时刻与线性***状态相关的输入变量，
Figure BDA00024193209000000222
表示与非线性***状态相关的过程噪声，
Figure BDA00024193209000000223
表示与线性***状态相关的过程噪声，
Figure BDA00024193209000000224
表示非线性状态转移矩阵，
Figure BDA00024193209000000225
表示与非线性状态相关的线性状态转移矩阵，
Figure BDA00024193209000000226
表示线性状态转移矩阵，
Figure BDA00024193209000000227
表示与线性状态相关的非线性状态转移矩阵，hk(·)表示与非线性状态相关的测量转移矩阵，，Bk(·)表示与线性状态相关的测量转移矩阵。
所述非线性状态不确定参数和线性状态不确定参数均为动态变参数时，参数离散模型为：
Figure BDA00024193209000000228
Figure BDA00024193209000000229
其中，
Figure BDA00024193209000000230
表示k时刻与非线性状态有关的不确定参数，
Figure BDA00024193209000000231
表示k-1时刻与非线性状态有关的不确定参数，
Figure BDA00024193209000000232
表示k时刻与线性状态有关的不确定参数，
Figure BDA00024193209000000233
表示k-1时刻与线性状态有关的不确定参数，
Figure BDA00024193209000000234
表示不确定参数
Figure BDA00024193209000000235
的转移矩阵，
Figure BDA00024193209000000236
表示不确定参数
Figure BDA00024193209000000237
的转移矩阵，ξk-1和εk-1均表示零均值高斯白噪声；
所述非线性状态不确定参数和线性状态不确定参数均为静态常参数时，参数离散模型为：
Figure BDA0002419320900000031
Figure BDA0002419320900000032
所述低维离散化混合模型和参数离散模型进行滤波处理，输出低维离散化混合模型中非线性状态的状态值的方法为：
S41、根据初始滤波误差协方差
Figure BDA0002419320900000033
和k-1时刻的状态观测器估计值
Figure BDA0002419320900000034
计算下一时刻的非线性状态滤波值和滤波误差协方差，具体方法为：
S41.1、生成初始容积点：
Figure BDA0002419320900000035
其中，[1]i表示一组n维状态空间单位矢量点，i＝1,2...M表示权值，且M＝2n，δi表示单位初始容积点，
Figure BDA0002419320900000036
表示非线性状态估计值，Sk-1表示估计误差协方差的平方根，χi,k-1表示***状态的初始容积点；
S41.2、计算非线性状态预测值
Figure BDA0002419320900000037
Figure BDA0002419320900000038
其中，
Figure BDA0002419320900000039
表示状态传递容积点，
Figure BDA00024193209000000310
表示k-1时刻的线性状态分量，
Figure BDA00024193209000000311
表示与非线性状态相关的不确定参数，
Figure BDA00024193209000000312
表示与非线性状态相关的输入变量；
S41.3、根据步骤S41.2中的非线性状态预测值
Figure BDA00024193209000000313
计算预测误差协方差
Figure BDA00024193209000000314
Figure BDA00024193209000000315
其中，Qn表示过程噪声
Figure BDA00024193209000000316
的方差；
S42、根据预测误差协方差
Figure BDA00024193209000000317
和非线性状态预测值
Figure BDA00024193209000000318
计算测量预测值
Figure BDA00024193209000000319
和测量估计误差协方差
Figure BDA00024193209000000320
具体方法为：
S42.1、生成预测容积点χi,k|k-1：
Figure BDA00024193209000000321
其中，Sk|k-1表示估计误差协方差的平方根因子，δi表示单位初始容积点；
S42.2、计算测量预测值
Figure BDA00024193209000000322
Figure BDA00024193209000000323
其中，
Figure BDA0002419320900000041
表示测量传递容积点，
Figure BDA0002419320900000042
表示k时刻的线性状态分量，hk(·)表示与非线性状态相关的测量转移矩阵，Bk(·)表示与线性状态相关的测量转移矩阵；
S42.3、根据测量预测值
Figure BDA0002419320900000043
计算测量预测误差协方差
Figure BDA0002419320900000044
Figure BDA0002419320900000045
其中，R为测量噪声υk的方差；
S43、根据步骤S41.2中的非线性状态预测值
Figure BDA0002419320900000046
和步骤S42.2中的测量预测值
Figure BDA0002419320900000047
计算预测互协方差
Figure BDA0002419320900000048
Figure BDA0002419320900000049
S44、利用预测互协方差
Figure BDA00024193209000000410
和测量误差协方差
Figure BDA00024193209000000411
更新k时刻的滤波误差协方差
Figure BDA00024193209000000412
并利用非线性状态预测值
Figure BDA00024193209000000413
测量值yk和测量预测值
Figure BDA00024193209000000414
计算低维离散化混合模型的非线性状态的状态值
Figure BDA00024193209000000415
Figure BDA00024193209000000416
Figure BDA00024193209000000417
其中，
Figure BDA00024193209000000418
本技术方案能产生的有益效果：
(1)本发明基于模型参数解耦，可以将不确定参数对***具体状态的影响进行直观的建模，可解决***模型参数不确定下的参数建模问题；
(2)本发明通过构建低维离散化混合模型，可将高维***分解为若干低维***的混合形式，通过这种混合，达到降低***维数的目的，减少后续估计过程的计算维数；
(3)本发明通过设计参数观测器对混合模型中的参数进行估计，克服了经典估计理论在***模型参数不确定下的稳定性降低甚至发散的问题；
(4)本发明基于混合模型和参数观测器的容积卡尔曼滤波算法，克服了传统容积卡尔曼滤波算法对精确***模型参数的依赖，且与无迹卡尔曼滤波算法和蒙特卡洛粒子滤波算法相比，基于容积规则的容积卡尔曼滤波算法所需的样本点数更少，计算量进一步减少。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明的流程图；
图2为本发明的参数观测器流程图；
图3为本发明基于参数观测器输出的容积卡尔曼滤波算法流程图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明实施例提供了一种高维耦合不确定***混合状态估计方法，具体步骤如下：
S1、构建高维深耦合不确定***状态模型、参数模型和测量模型，通用形式为：
Figure BDA0002419320900000051
Figure BDA0002419320900000052
y(t)＝h(x(t),t)+υ(t) (3)，
其中，
Figure BDA0002419320900000053
表示***状态的一阶导数，f(·)为非线性状态转移矩阵，x(t)表示随时间t变化的状态变量，θ(t)表示不确定参数，u(t)表示***输入变量，ω1(t)为状态过程噪声，
Figure BDA0002419320900000054
表示表示***不确定参数的一阶导数，ω2(t)为参数过程噪声，y(t)为测量值，h(·)为测量转移矩阵，υ(t)为测量噪声。
S2、根据***状态模型、参数模型和测量模型分别构建状态和参数的观测器模型，并利用状态和参数的观测器模型估计***状态和参数的估计值；假设***是完全能控和能观的。那么可基于观测器理论分别设计如图2所示的观测器模型。其中，观测器增益分别与不同的观测器相对应。构建的观测器模型为：
Figure BDA0002419320900000055
Figure BDA0002419320900000056
其中，
Figure BDA0002419320900000057
表示***状态估计值的一阶导数，
Figure BDA0002419320900000058
表示非线性状态估计值，
Figure BDA0002419320900000059
表示参数估计值的一阶导数，
Figure BDA00024193209000000510
表示输入变量估计值，K1和K2表示对应的观测器增益矩阵，
Figure BDA00024193209000000511
表示***实际测量值与观测器的输出测量值之间的偏差，
Figure BDA00024193209000000512
表示观测器的输出测量值。
S3、利用欧拉法对步骤S1中的***状态模型、参数模型和测量模型进行离散化处理，得到低维离散化混合模型：
Figure BDA0002419320900000061
Figure BDA0002419320900000062
Figure BDA0002419320900000063
其中，非线性状态和线性状态是相互影响的，
Figure BDA0002419320900000064
表示k时刻的非线性状态分量，
Figure BDA0002419320900000065
表示k-1时刻的非线性状态分量，
Figure BDA0002419320900000066
表示k时刻的线性状态分量，
Figure BDA0002419320900000067
表示k-1时刻的线性状态分量，
Figure BDA0002419320900000068
表示非线性状态不确定参数，
Figure BDA0002419320900000069
表示线性状态不确定参数，
Figure BDA00024193209000000610
表示k时刻与非线性***状态相关的输入变量，
Figure BDA00024193209000000611
表示k时刻与线性***状态相关的输入变量，
Figure BDA00024193209000000612
表示均值为0，方差为Qn的非线性过程噪声，
Figure BDA00024193209000000613
表示均值为0，方差为Ql的线性过程噪声，
Figure BDA00024193209000000614
表示非线性状态转移矩阵，
Figure BDA00024193209000000615
表示与非线性状态相关的线性状态转移矩阵，
Figure BDA00024193209000000616
表示线性状态转移矩阵，
Figure BDA00024193209000000617
表示与线性状态相关的非线性状态转移矩阵，hk(·)表示与非线性状态相关的测量转移矩阵，Bk(·)表示与线性状态相关的测量转移矩阵。
所述非线性状态不确定参数和线性状态不确定参数均为动态变参数时，对低维离散化混合模型中的不确定参数进行建模，得到的参数离散模型为：
Figure BDA00024193209000000618
Figure BDA00024193209000000619
其中，
Figure BDA00024193209000000620
表示k时刻与非线性状态相关的不确定参数，
Figure BDA00024193209000000621
表示k时刻与线性状态相关的不确定参数，
Figure BDA00024193209000000622
表示不确定参数
Figure BDA00024193209000000623
的转移矩阵，
Figure BDA00024193209000000624
表示不确定参数
Figure BDA00024193209000000625
的转移矩阵，ξk-1和εk-1均表示零均值高斯白噪声；
所述非线性状态不确定参数和线性状态不确定参数均为静态常参数时，对低维离散化混合模型中的不确定参数进行建模，得到的参数离散模型为：
Figure BDA00024193209000000626
Figure BDA00024193209000000627
S4、基于步骤S2中的***状态和参数的估计值，利用容积卡尔曼滤波算法分别对步骤S3中的低维离散化混合模型和参数离散模型进行滤波处理，根据公式(4)得到k时刻的状态估计值
Figure BDA00024193209000000628
和
Figure BDA00024193209000000629
根据公式(5)得到参数估计值
Figure BDA00024193209000000630
和
Figure BDA00024193209000000631
并将
Figure BDA00024193209000000632
和
Figure BDA00024193209000000633
作为k时刻的低维离散化混合模型的修正信号来修正容积卡尔曼滤波算法，最终输出低维离散化混合模型的状态值，以非线性状态变量
Figure BDA00024193209000000634
为例，如图3所示，具体方法为：
S41、根据初始估计误差协方差
Figure BDA00024193209000000635
和k-1时刻的非线性状态估计值
Figure BDA00024193209000000636
计算非线性状态预测值和估计误差协方差，具体方法为：
S41.1、生成初始容积点：
Figure BDA0002419320900000071
其中，[1]i表示一组n维状态空间单位矢量点，i＝1,2...M＝2n表示权值，δi表示单位初始容积点，
Figure BDA0002419320900000072
表示非线性状态估计值，Sk-1表示估计误差协方差的平方根，经由
Figure BDA0002419320900000073
计算得到，χi,k-1表示***状态的初始容积点；
S41.2、计算非线性状态预测值
Figure BDA0002419320900000074
Figure BDA0002419320900000075
其中，
Figure BDA0002419320900000076
表示状态传递容积点，
Figure BDA0002419320900000077
表示k-1时刻的线性状态分量，
Figure BDA0002419320900000078
表示与非线性状态相关的不确定参数，
Figure BDA0002419320900000079
表示与非线性状态相关的输入变量估计值，
Figure BDA00024193209000000710
步骤S3中对应的与非线性状态相关的参观测器的输出值，
Figure BDA00024193209000000711
表示k-1时刻步骤S3中线性状态观测器输出值；
S41.3、根据步骤S41.2中的非线性状态预测值
Figure BDA00024193209000000712
计算预测误差协方差
Figure BDA00024193209000000713
Figure BDA00024193209000000714
其中，Qn表示过程噪声
Figure BDA00024193209000000715
的方差；
S42、根据预测误差协方差
Figure BDA00024193209000000716
和非线性状态预测值
Figure BDA00024193209000000717
计算测量预测值
Figure BDA00024193209000000718
和测量估计误差协方差
Figure BDA00024193209000000719
具体方法为：
S42.1、生成预测容积点χi,k|k-1：
Figure BDA00024193209000000720
其中，Sk|k-1表示估计误差协方差的平方根因子，经由
Figure BDA00024193209000000721
计算得到，δi表示单位初始容积点；
S42.2、计算测量预测值
Figure BDA00024193209000000722
Figure BDA00024193209000000723
其中，
Figure BDA00024193209000000724
表示测量传递容积点，
Figure BDA00024193209000000725
表示k时刻步骤S3中线性状态观测器输出值，hk(·)表示与非线性状态相关的测量转移矩阵，Bk(·)表示与线性状态相关的测量转移矩阵；
S42.3、根据测量预测值
Figure BDA0002419320900000081
计算测量误差协方差
Figure BDA0002419320900000082
Figure BDA0002419320900000083
其中，R为测量噪声υk的方差；
S43、根据步骤S41.2中的非线性状态预测值
Figure BDA0002419320900000084
和步骤S42.2中的测量预测值
Figure BDA0002419320900000085
计算预测互协方差
Figure BDA0002419320900000086
Figure BDA0002419320900000087
S44、利用预测互协方差
Figure BDA0002419320900000088
和测量误差协方差
Figure BDA0002419320900000089
更新k时刻的估计误差协方差
Figure BDA00024193209000000810
并利用非线性状态预测值
Figure BDA00024193209000000811
测量值yk和测量预测值
Figure BDA00024193209000000812
计算低维离散化混合模型的非线性状态估计值
Figure BDA00024193209000000813
Figure BDA00024193209000000814
Figure BDA00024193209000000815
其中，
Figure BDA00024193209000000816
将步骤S44中得到的低维离散化混合模型的状态值
Figure BDA00024193209000000817
进行存储并输出，将估计误差协方差
Figure BDA00024193209000000818
进行存储并反馈到步骤S41执行下一时刻的状态估计，循环运行，实现高维深耦合不确定***混合非线性状态
Figure BDA00024193209000000819
的估计。
同理，采用步骤S41至步骤S44的处理方法可得到高维深耦合不确定***混合线性状态
Figure BDA00024193209000000820
的估计。
基于“当前”统计模型的机动目标模型在目标跟踪领域有着广泛的应用，在三维空间中目标的位置、速度和加速度构成了一组9维空间状态模型。
为了说明发明的实施过程，选取任意一维空间的模型进行具体说明：
1.一维机动目标连续状态方程如下：
Figure BDA00024193209000000821
其中，
Figure BDA00024193209000000822
x(t)表示当前目标的位置，
Figure BDA00024193209000000823
表示当前目标的速度，
Figure BDA00024193209000000824
表示当前目标的加速度，ω1(t)为状态噪声干扰，在“当前”统计模型条件下，当目标以某一加速度机动时，通常可以采用非零均值的时间相关模型来表示加速度的机动，具体如下：
Figure BDA00024193209000000825
a(t)＝-αa(t)+ω2(t) (24)，
其中，
Figure BDA0002419320900000091
为加速度均值，a(t)为加速度的零均值有色噪声干扰，α为时间常数的倒数，ω2(t)为零均值白噪声。
2.一维空间测量方程为：
y(t)＝h(X(t))+υ1(t) (25)，
其中，y(t)表示测量的位置，υ1(t)表示测量噪声干扰，h(·)表示测量状态转移矩阵。
3.针对***模型(22)和(24)，设计位置、速度、加速度以及加速度有色噪声的观测器如下：
Figure BDA0002419320900000092
Figure BDA0002419320900000093
通过选取合适的观测器增益，保证观测器的全局稳定性，进而由观测器估计出目标的位置、速度、加速度以及加速度干扰的状态估计值。
4.然后对公式(26)进行离散化处理，采样时间为△t，得到离散状态***方程如下：
Figure BDA0002419320900000094
通过变换可转化为以下混合模型状态方程：
Figure BDA0002419320900000095
Figure BDA0002419320900000096
加速度干有色噪声干扰的离散化模型为：
ak+1＝ak+ω6,k (31)；
测量离散化模型为：
yk＝[1 0 0]Xk+υk (32)；
因此，根据公式(30)、(31)和(32)可知低维离散化混合模型中的变量分别为：
Figure BDA0002419320900000101
Figure BDA0002419320900000102
Figure BDA0002419320900000103
Figure BDA0002419320900000104
Figure BDA0002419320900000105
Figure BDA0002419320900000106
Figure BDA0002419320900000107
Figure BDA0002419320900000108
Figure BDA0002419320900000109
最后，利用容积卡尔曼滤波算法对以上参数进行滤波，即可得到低维离散化混合模型的状态值。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种高维耦合不确定***混合状态估计方法，其特征在于，其步骤如下：
S1、构建高维深耦合不确定***状态模型、参数模型和测量模型；
S2、根据***状态模型、参数模型和测量模型分别构建状态和参数的观测器模型，并利用状态和参数的观测器模型估计***状态和参数的估计值；
S3、利用欧拉法对步骤S1中的***状态模型、参数模型和测量模型进行离散化处理，得到低维离散化混合模型，并对低维离散化混合模型中的不确定参数进行建模，得到参数离散模型；
S4、基于步骤S2中的***状态和参数的估计值，利用容积卡尔曼滤波算法分别对步骤S3中的低维离散化混合模型和参数离散模型进行滤波处理，输出低维离散化混合模型的状态值。
2.根据权利要求1所述的高维耦合不确定***混合状态估计方法，其特征在于，所述步骤S1中的高维深耦合不确定***状态模型、参数模型和测量模型分别为：
Figure FDA0002419320890000011
Figure FDA0002419320890000012
y(t)＝h(x(t),t)+υ(t) (3)，
其中，
Figure FDA0002419320890000013
表示***状态的一阶导数，f(·)为非线性状态转移矩阵，x(t)表示随时间t变化的状态变量，θ(t)表示不确定参数，u(t)表示***输入变量，ω1(t)为状态过程噪声，
Figure FDA0002419320890000014
表示***不确定参数的一阶导数，ω2(t)为参数过程噪声，y(t)为测量值，h(·)为测量转移矩阵，υ(t)为测量噪声。
3.根据权利要求1或2所述的高维耦合不确定***混合状态估计方法，其特征在于，所述根据***状态模型、参数模型和测量模型分别构建状态和参数的观测器模型分别为：
Figure FDA0002419320890000015
Figure FDA0002419320890000016
其中，
Figure FDA0002419320890000017
表示***状态估计值的一阶导数，
Figure FDA0002419320890000018
表示线性状态估计值，
Figure FDA0002419320890000019
表示参数估计值的一阶导数，
Figure FDA00024193208900000110
表示输入变量估计值，K1和K2表示对应的观测器增益矩阵，
Figure FDA00024193208900000111
表示***实际测量值与观测器的输出测量值之间的偏差，
Figure FDA00024193208900000112
表示观测器的输出测量值。
4.根据权利要求1所述的高维耦合不确定***混合状态估计方法，其特征在于，所述低维离散化混合模型为：
Figure FDA00024193208900000113
Figure FDA00024193208900000114
Figure FDA0002419320890000021
其中，
Figure FDA0002419320890000022
表示k时刻的非线性状态分量，
Figure FDA0002419320890000023
表示k-1时刻的非线性状态分量，
Figure FDA0002419320890000024
表示k时刻的线性状态分量，
Figure FDA0002419320890000025
表示k-1时刻的线性状态分量，
Figure FDA0002419320890000026
表示k-1时刻的非线性状态不确定参数，
Figure FDA0002419320890000027
表示k-1时刻的线性状态不确定参数，
Figure FDA0002419320890000028
表示k时刻与非线性***状态相关的输入变量，
Figure FDA0002419320890000029
表示k时刻与线性***状态相关的输入变量，
Figure FDA00024193208900000210
表示与非线性***状态相关的过程噪声，
Figure FDA00024193208900000211
表示与线性***状态相关的过程噪声，
Figure FDA00024193208900000212
表示非线性状态转移矩阵，
Figure FDA00024193208900000213
表示与非线性状态相关的线性状态转移矩阵，
Figure FDA00024193208900000214
表示线性状态转移矩阵，
Figure FDA00024193208900000215
表示与线性状态相关的非线性状态转移矩阵，hk(·)表示与非线性状态相关的测量转移矩阵，，Bk(·)表示与线性状态相关的测量转移矩阵。
5.根据权利要求4所述的高维耦合不确定***混合状态估计方法，其特征在于，所述非线性状态不确定参数和线性状态不确定参数均为动态变参数时，参数离散模型为：
Figure FDA00024193208900000216
Figure FDA00024193208900000217
其中，
Figure FDA00024193208900000218
表示k时刻与非线性状态有关的不确定参数，
Figure FDA00024193208900000219
表示k-1时刻与非线性状态有关的不确定参数，
Figure FDA00024193208900000220
表示k时刻与线性状态有关的不确定参数，
Figure FDA00024193208900000221
表示k-1时刻与线性状态有关的不确定参数，
Figure FDA00024193208900000222
表示不确定参数
Figure FDA00024193208900000223
的转移矩阵，
Figure FDA00024193208900000224
表示不确定参数
Figure FDA00024193208900000225
的转移矩阵，ξk-1和εk-1均表示零均值高斯白噪声；
所述非线性状态不确定参数和线性状态不确定参数均为静态常参数时，参数离散模型为：
Figure FDA00024193208900000226
Figure FDA00024193208900000227
6.根据权利要求1所述的高维耦合不确定***混合状态估计方法，其特征在于，所述低维离散化混合模型和参数离散模型进行滤波处理，输出低维离散化混合模型中非线性状态的状态值的方法为：
S41、根据初始滤波误差协方差
Figure FDA00024193208900000228
和k-1时刻的状态观测器估计值
Figure FDA00024193208900000229
计算下一时刻的非线性状态滤波值和滤波误差协方差，具体方法为：
S41.1、生成初始容积点：
Figure FDA00024193208900000230
其中，[1]i表示一组n维状态空间单位矢量点，i＝1,2...M表示权值，且M＝2n，δi表示单位初始容积点，
Figure FDA00024193208900000231
表示非线性状态估计值，Sk-1表示估计误差协方差的平方根，χi,k-1表示***状态的初始容积点；
S41.2、计算非线性状态预测值
Figure FDA0002419320890000031
Figure FDA0002419320890000032
其中，
Figure FDA0002419320890000033
表示状态传递容积点，
Figure FDA0002419320890000034
表示k-1时刻的线性状态分量，
Figure FDA0002419320890000035
表示与非线性状态相关的不确定参数，
Figure FDA0002419320890000036
表示与非线性状态相关的输入变量；
S41.3、根据步骤S41.2中的非线性状态预测值
Figure FDA0002419320890000037
计算预测误差协方差
Figure FDA0002419320890000038
Figure FDA0002419320890000039
其中，Qn表示过程噪声
Figure FDA00024193208900000310
的方差；
S42、根据预测误差协方差
Figure FDA00024193208900000311
和非线性状态预测值
Figure FDA00024193208900000312
计算测量预测值
Figure FDA00024193208900000313
和测量估计误差协方差
Figure FDA00024193208900000314
具体方法为：
S42.1、生成预测容积点χi,k|k-1：
Figure FDA00024193208900000315
其中，Sk|k-1表示估计误差协方差的平方根因子，δi表示单位初始容积点；
S42.2、计算测量预测值
Figure FDA00024193208900000316
Figure FDA00024193208900000317
其中，
Figure FDA00024193208900000318
表示测量传递容积点，
Figure FDA00024193208900000319
表示k时刻的线性状态分量，hk(·)表示与非线性状态相关的测量转移矩阵，Bk(·)表示与线性状态相关的测量转移矩阵；
S42.3、根据测量预测值
Figure FDA00024193208900000320
计算测量预测误差协方差
Figure FDA00024193208900000321
Figure FDA00024193208900000322
其中，R为测量噪声υk的方差；
S43、根据步骤S41.2中的非线性状态预测值
Figure FDA00024193208900000323
和步骤S42.2中的测量预测值
Figure FDA00024193208900000324
计算预测互协方差
Figure FDA00024193208900000325
Figure FDA00024193208900000326
S44、利用预测互协方差
Figure FDA00024193208900000327
和测量误差协方差
Figure FDA00024193208900000328
更新k时刻的滤波误差协方差
Figure FDA00024193208900000329
并利用非线性状态预测值
Figure FDA00024193208900000330
测量值yk和测量预测值
Figure FDA00024193208900000331
计算低维离散化混合模型的非线性状态的状态值
Figure FDA0002419320890000041
Figure FDA0002419320890000042
Figure FDA0002419320890000043
其中，
Figure FDA0002419320890000044
CN202010205361.5A 2020-03-20 2020-03-20 一种高维耦合不确定***混合状态估计方法 Active CN111310110B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010205361.5A CN111310110B (zh) 2020-03-20 2020-03-20 一种高维耦合不确定***混合状态估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010205361.5A CN111310110B (zh) 2020-03-20 2020-03-20 一种高维耦合不确定***混合状态估计方法

Publications (2)

Publication Number Publication Date
CN111310110A true CN111310110A (zh) 2020-06-19
CN111310110B CN111310110B (zh) 2023-05-05

Family

ID=71149858

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010205361.5A Active CN111310110B (zh) 2020-03-20 2020-03-20 一种高维耦合不确定***混合状态估计方法

Country Status (1)

Country Link
CN (1) CN111310110B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113222268A (zh) * 2021-05-24 2021-08-06 郑州轻工业大学 一种基于多模式推理的烟草烘烤质量预测模型建立方法
CN113345532A (zh) * 2021-06-01 2021-09-03 江南大学 聚苯醚生产过程的快速状态估计方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5991525A (en) * 1997-08-22 1999-11-23 Voyan Technology Method for real-time nonlinear system state estimation and control
CN106646356A (zh) * 2016-11-23 2017-05-10 西安电子科技大学 一种基于卡尔曼滤波定位的非线性***状态估计方法
WO2018014602A1 (zh) * 2016-07-19 2018-01-25 东南大学 适于高维gnss/ins深耦合的容积卡尔曼滤波方法
CN110824363A (zh) * 2019-10-21 2020-02-21 江苏大学 一种基于改进ckf的锂电池soc和soe联合估算方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5991525A (en) * 1997-08-22 1999-11-23 Voyan Technology Method for real-time nonlinear system state estimation and control
WO2018014602A1 (zh) * 2016-07-19 2018-01-25 东南大学 适于高维gnss/ins深耦合的容积卡尔曼滤波方法
CN106646356A (zh) * 2016-11-23 2017-05-10 西安电子科技大学 一种基于卡尔曼滤波定位的非线性***状态估计方法
CN110824363A (zh) * 2019-10-21 2020-02-21 江苏大学 一种基于改进ckf的锂电池soc和soe联合估算方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
许大星等: "未知状态模型下基于高阶容积卡尔曼滤波和神经网络的状态估计算法", 《计算机应用与软件》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113222268A (zh) * 2021-05-24 2021-08-06 郑州轻工业大学 一种基于多模式推理的烟草烘烤质量预测模型建立方法
CN113222268B (zh) * 2021-05-24 2023-04-07 郑州轻工业大学 一种基于多模式推理的烟草烘烤质量预测模型建立方法
CN113345532A (zh) * 2021-06-01 2021-09-03 江南大学 聚苯醚生产过程的快速状态估计方法
CN113345532B (zh) * 2021-06-01 2022-07-15 江南大学 聚苯醚生产过程的快速状态估计方法

Also Published As

Publication number Publication date
CN111310110B (zh) 2023-05-05

Similar Documents

Publication Publication Date Title
Särkkä et al. Bayesian filtering and smoothing
CN105785999B (zh) 无人艇航向运动控制方法
CN103217175B (zh) 一种自适应容积卡尔曼滤波方法
Singer Parameter estimation of nonlinear stochastic differential equations: simulated maximum likelihood versus extended Kalman filter and Itô-Taylor expansion
CN104112079A (zh) 一种模糊自适应变分贝叶斯无迹卡尔曼滤波方法
KR100816269B1 (ko) 언센티드 필터를 적용한 강인한 동시 위치 추정 및 지도작성 방법
CN107179693B (zh) 基于Huber估计的鲁棒自适应滤波和状态估计方法
CN111310110A (zh) 一种高维耦合不确定***混合状态估计方法
CN103940433A (zh) 一种基于改进的自适应平方根ukf算法的卫星姿态确定方法
CN109827579B (zh) 一种组合定位中滤波模型实时校正的方法和***
CN104019817A (zh) 一种用于卫星姿态估计的范数约束强跟踪容积卡尔曼滤波方法
CN104283529B (zh) 未知测量噪声方差的平方根高阶容积卡尔曼滤波方法
CN113587926A (zh) 一种航天器空间自主交会对接相对导航方法
CN109341690B (zh) 一种鲁棒高效的组合导航自适应数据融合方法
Qian et al. An INS/DVL integrated navigation filtering method against complex underwater environment
CN109582915B (zh) 应用于纯方位跟踪的改进非线性可观测度自适应滤波方法
CN114139109A (zh) 一种目标跟踪方法、***、设备、介质及数据处理终端
CN112986977A (zh) 一种克服雷达扩展卡尔曼航迹滤波发散的方法
CN110912535B (zh) 一种新型无先导卡尔曼滤波方法
CN109655057B (zh) 一种六推无人机加速器测量值的滤波优化方法及其***
CN109188422B (zh) 一种基于lu分解的卡尔曼滤波目标跟踪方法
CN115580261A (zh) 基于鲁棒策略和动态增强ckf结合的状态估计方法
CN111736144A (zh) 一种仅用距离观测的机动转弯目标状态估计方法
CN115655285A (zh) 一种权值及参考四元数修正的无迹四元数姿态估计方法
CN115438728A (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