CN112115419A - 系统状态估计方法、系统状态估计装置 - Google Patents

系统状态估计方法、系统状态估计装置 Download PDF

Info

Publication number
CN112115419A
CN112115419A CN202010963665.8A CN202010963665A CN112115419A CN 112115419 A CN112115419 A CN 112115419A CN 202010963665 A CN202010963665 A CN 202010963665A CN 112115419 A CN112115419 A CN 112115419A
Authority
CN
China
Prior art keywords
state
estimation
observation
covariance
equation
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
CN202010963665.8A
Other languages
English (en)
Other versions
CN112115419B (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 CN202010963665.8A priority Critical patent/CN112115419B/zh
Publication of CN112115419A publication Critical patent/CN112115419A/zh
Application granted granted Critical
Publication of CN112115419B publication Critical patent/CN112115419B/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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F41—WEAPONS
    • F41G—WEAPON SIGHTS; AIMING
    • F41G3/00—Aiming or laying means
    • G—PHYSICS
    • G08—SIGNALLING
    • G08G—TRAFFIC CONTROL SYSTEMS
    • G08G5/00—Traffic control systems for aircraft, e.g. air-traffic control [ATC]

Landscapes

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

Abstract

本申请公开了一种系统状态估计方法、装置，用于对目标系统的状态进行估计，该方法包括：获取目标系统的状态方程和观测方程；根据目标系统的上一时刻的协方差，对所述状态方程和观测方程进行无迹变换，得到状态先验估计、协方差先验估计、观测先验估计，以及状态和观测之间的协方差；根据所述状态先验估计、协方差先验估计、观测先验估计和状态和观测之间的协方差，构建线性回归方程；根据所述线性回归方程计算估计误差；采用基于模糊交叉熵构造的损失函数对所述估计误差进行优化，得到所述目标系统的状态估计，该方法提高了无迹卡尔曼滤波算法在非高斯噪声环境中的鲁棒性，提高无迹卡尔曼滤波算法的估算精度。

Description

系统状态估计方法、系统状态估计装置
技术领域
本申请涉及卡尔曼滤波算法领域，尤其涉及一种系统状态估计方法、系 统状态估计装置。
背景技术
估计问题在精确制导、预警系统、航空交通以及智能监控等领域发挥着 关键作用。卡尔曼滤波算法(Kalman Filter,KF)是一种经典的状态估计方法， 可以对线性系统进行最小均方误差估计。然而，卡尔曼滤波算法并不适用于 非线性系统，因此研究人员相继提出了一些优化卡尔曼滤波算法的方法，例 如无迹卡尔曼滤波(Unscented KalmanFilter,UKF)算法，无迹卡尔曼滤波算法 主要是利用一组精确的点集逼近状态的概率分布，并通过非线性方程进行传 播。无迹卡尔曼滤波算法是近年来滤波的主要方法之一，但是，由于无迹卡 尔曼滤波算法是基于最小均方误差估计的准则进行优化的，使得系统在在非 高斯噪声环境中的鲁棒性下降，从而导致当系统受到非高斯噪声的干扰时， 无迹卡尔曼滤波算法的性能急剧下降。
申请内容
本申请实施例提供一种系统状态估计方法、系统状态估计装置，以提高 无迹卡尔曼滤波算法在非高斯噪声环境中的鲁棒性，提高估算的精准度。
一种系统状态估计方法，用于对目标系统的状态进行估计，包括：
获取目标系统的状态方程和观测方程；
根据目标系统的上一时刻的协方差，对所述状态方程和观测方程进行无 迹变换，得到状态先验估计、协方差先验估计、观测先验估计，以及状态和 观测之间的协方差；
根据所述状态先验估计、协方差先验估计、观测先验估计和状态和观测 之间的协方差，构建线性回归方程；
根据所述线性回归方程计算估计误差；
采用基于模糊交叉熵构造的损失函数对所述估计误差进行优化，得到所 述目标系统的状态估计。
优选地，所述状态方程包括状态函数，且所述状态方程为n维状态向量；
所述方法通过以下步骤计算所述状态先验估计和协方差先验估计：
根据所述上一时刻的协方差和所述状态方程获取k-1时刻的多个第一 sigma点集，且所述k-1时刻的第一sigma点集为：
Figure BDA0002681456900000021
Figure BDA0002681456900000022
其中，
Figure BDA0002681456900000023
是矩阵(n+λ)P(k-1|k-1)平方根的第i列， P(k-1|k-1)为所述上一时刻的协方差，n为状态维度，λ为复合比例因子，且 λ＝α2(n+φ)-n，α和φ均为预设值；
根据所述状态函数，将所述k-1时刻的第一sigma点集转换为k时刻的多 个第二sigma点集：
χi*(k|k-1)＝f(k-1,χi(k-1|k-1)),for i＝0...2n
其中，χi*(k|k-1)表示第二sigma点集，f表示状态函数；
根据所述k时刻的多个第二sigma点集，计算所述状态先验估计和协方 差先验估计：
Figure BDA0002681456900000031
Figure BDA0002681456900000032
其中，所述
Figure BDA0002681456900000033
为状态先验估计，P(k|k-1)为协方差先验估计，
Figure BDA0002681456900000034
Q(k-1)为预设的状态协方差矩阵。
优选地，所述观测方程包括观测函数；所述方法通过以下步骤计算所述 观测先验估计，以及状态和观测之间的协方差：
根据所述状态先验估计和协方差先验估计，获取所述k时刻的多个第三 sigma点集，且所述第三sigma点集为：
Figure BDA0002681456900000035
Figure BDA0002681456900000036
根据所述观测函数，将所述k时刻的第三sigma点集转换为k时刻的第 四sigma点集：
γi(k)＝h(k,χi(k|k-1)),for i＝0...2n
其中，γi(k)为第四sigma点集，h()为观测函数；
根据所述k时刻的第四sigma点集，计算所述观测先验估计：
Figure BDA0002681456900000041
根据所述第二sigma点集、状态先验估计、第四sigma点集和观测先验估 计计算所述状态和观测之间的协方差：
Figure BDA0002681456900000042
优选地，所述观测方程包括观测噪声；
通过以下步骤构建所述线性回归方程，包括：
获取观测斜率矩阵，所述观测斜率矩阵为：
H(k)＝(P-1(k|k-1)Pxy(k))T；
根据所述观测斜率矩阵计算线性回归方程，所述线性回归方程为：
Figure BDA0002681456900000043
其中，所述I为n×n的单位矩阵，
Figure BDA0002681456900000044
且r(k)为观测噪 声，线性回归方程的协方差为
Figure BDA0002681456900000045
R(k)为预设观测协方差 矩阵。
优选地，所述估计误差为：
e(k)＝D(k)-W(k)x(k)
其中，所述e(k)为估计误差，
Figure BDA0002681456900000046
Figure BDA0002681456900000047
x(k)为状态方程，S(k)通过对所述线性回归方程的协方差 进行Cholesky分解得到。
优选地，所述损失函数为：
Figure BDA0002681456900000051
其中，a为预设值，σ为高斯核函数的核宽，μik为k时刻第i维的模糊隶属 度，ei(k)＝di(k)-wi(k)x(k)，ei(k)是e(k)的第i个元素，di(k)是D(k)的第i个元素， wi(k)是W(k)的第i行元素，L是e(k)的维数，
Figure BDA0002681456900000052
根据所述损失函数计算得到所述目标系统的状态估计为：
Figure BDA0002681456900000053
优选地，所述σ的取值为：
Figure BDA0002681456900000054
其中，σi是第i维误差的核宽，ei是第i维的误差，σ为核宽的预设值。
8.如权利要求6所述的系统状态估计方法，其特征在于，在得到所述目 标系统的状态，计算所述目标系统的协方差后验估计：
Figure BDA0002681456900000055
其中，
Figure BDA0002681456900000056
Figure BDA0002681456900000057
Figure BDA0002681456900000058
Sp(k|k-1)为对协方差先验估计P(k|k-1)进行Cholesky分解得到，Sr(k)为 对R(k)进行Cholesky分解得到。
优选地，所述μik的取值为：
Figure BDA0002681456900000061
一种系统状态估计装置，用于对目标系统的状态进行估计，包括：
获取单元，用于获取目标系统的状态方程和观测方程；
无迹变换单元，用于根据上一时刻的协方差，对所述状态方程和观测方 程进行无迹变换，得到状态先验估计、协方差先验估计、观测先验估计，以 及状态和观测之间的协方差；
线性方程构建单元，用于根据所述状态先验估计、协方差先验估计、观 测先验估计和状态和观测之间的协方差，构建线性回归方程；
误差估计单元，用于根据所述线性回归方程计算估计误差；
优化单元，用于采用基于模糊交叉熵构造的损失函数对所述估计误差进 行优化，得到所述目标系统的状态。
上述系统状态估计方法、装置通过对目标系统进行无迹变换，以获得状 态先验估计、协方差的先验估计、观测先验估计，以及状态和观测之间的协 方差，并根据得到的数据构建线性回归方程并获得估计误差，最后采用最大 模糊交叉熵准则优化无迹卡尔曼滤波算法得到的估计误差，使得该状态估计 方法可以较好地处理非线性非高斯系统，提高了无迹卡尔曼滤波算法在非高 斯噪声环境中的鲁棒性，提高无迹卡尔曼滤波算法的估算精度。
附图说明
为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例的 描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅 仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性 劳动性的前提下，还可以根据这些附图获得其他的附图。
图1是本申请一实施例中系统状态估计方法的流程图；
图2是本申请实施例中的目标真实运动轨迹图；
图3是各个算法的跟踪轨迹图；
图4是各个算法的均方根误差对比图；
图5是本申请一实施例中系统状态估计装置的原理框图。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行 清楚、完整地描述，显然，所描述的实施例是本申请一部分实施例，而不是 全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创 造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
在本申请的描述中，需要说明的是，术语“第一”、“第二”、“第三”等仅用 于区分描述，而不能理解为指示或暗示相对重要性。
本申请提出了一种系统状态估计方法，该方法首先基于模糊信息理论构 建了一种模糊交叉熵，利用模糊隶属度较好地表示不同样本对状态估计的不 同影响，然后采用最大模糊交叉熵准则(MFCC)优化无迹卡尔曼滤波算法计算 得到的误差，以获得最优的状态估计。该方法可以用于对目标系统进行状态 估计，其中，目标系统可以是精确制导、预警系统、航空交通以及智能监控 等领域的系统，本实施例对目标系统不做具体限定。
首先介绍本申请实施例所构建的模糊交叉熵：
交叉熵是两个随机变量之间的广义相似性度量，定义如下：
V(X,Y)＝E[κ(X,Y)]＝∫κ(x,y)dFXY(x,y) (1)
其中，
Figure BDA0002681456900000081
是两个随机变量，FXY(x,y)是这两个变量的联合分布函数， E是数学期望，κ(·,·)是一个具有平移不变性的Mercer核函数。
本文所用的核函数为高斯核函数：
Figure BDA0002681456900000082
其中，e＝x-y，核宽σ>0。
在实际情况中，我们通常只能得到一组有限的数据而不知道他们的联合 分布FXY。为此，我们用样本均值来估计交叉熵：
Figure BDA0002681456900000083
其中，e(i)＝x(i)-y(i)，
Figure BDA0002681456900000084
是从FXY中抽取的N个样本。
从交叉熵的定义可以看出，其对于所有样本都具有相同的权值1/N。而在 实际情况中，不同样本对于状态估计的作用应该是不尽相同的，也就是说， 不同的样本应该具有不同的权值。对此，基于模糊信息处理理论，定义如下 模糊交叉熵：
Figure BDA0002681456900000085
其中，a为加权指数，μi表示变量x(i)和y(i)之间的模糊隶属度，且满足如 下条件：
Figure BDA0002681456900000086
在了解了模糊交叉熵之后，本实施将从获得目标系统的状态方程和观测 方程开始介绍本申请所提出的系统状态估计方法。如图1所示，该方法包括 如下步骤：
S10：获取目标系统的状态方程和观测方程。其中，这里所说的目标系统 可以是线性系统，也可以是非线性系统，目标系统的状态方程和观测可以根 据系统的运行参数或者其他方式获得，由于如何获得目标系统的状态方程和 观测方程并非本实施例的重点，因此这里便不展开论述。本实施例将以目标 系统为非线性系统为例进行解释。可以理解地，本实施例所提出的系统状态 估计方法可以用来对非线性系统的状态进行估计，自然而然地，也能对线性 系统的状态进行估计。示例性地，目标系统的状态方程x(k)和观测方程y(k)分别为：
x(k)＝f(k-1,x(k-1))+q(k-1) (1)
y(k)＝h(k,x(k))+r(k) (2)
其中，
Figure BDA0002681456900000091
是k时刻的n维状态向量，
Figure BDA0002681456900000092
是k时刻的m维观测向量， f是状态函数，h是观测函数。过程噪声q(k-1)和观测噪声r(k)的均值为零， 噪声协方差矩阵为Q(k-1)和观测协方差矩阵为R(k)，且：
E[q(k-1)qT(k-1)]＝Q(k-1),E[r(k)rT(k)]＝R(k) (3)
S20：根据目标系统的上一时刻的协方差，对状态方程和观测方程进行无 迹变换，得到状态先验估计、协方差先验估计、观测先验估计，以及状态和 观测之间的协方差。
其中，目标系统的上一时刻的协方差为P(k-1|k-1)，在一些情况下，上 一时刻的协方差为初始时刻的协方差时，此时通常为为预设值，可以设置为 单位矩阵，并且会附加一些噪声，当然，本实施例的上一时刻的协方差矩阵 并不限定于单位矩阵，具体可根据需求设置。对状态方程和观测方程进行无 迹变换，主要包括时间更新和状态更新。
时间更新：
从k-1时刻的状态方程和上一时刻的协方差中生成2n+1个第一sigma点 集：
Figure BDA0002681456900000101
其中
Figure BDA0002681456900000102
是矩阵(n+λ)P(k-1|k-1)平方根的第i列，n为状 态维度，λ为复合比例因子，定义如下：
λ＝α2(n+φ)-n (10)
其中，α为预设值，其确定了sigma点的分布，通常选择为一个小的正数； φ也为预设值，其设置可以是为3-n。
通过状态函数f将第一sigma点集转换为第二sigma点集：
χi*(k|k-1)＝f(k-1,χi(k-1|k-1)),for i＝0...2n (11)
则状态先验估计
Figure BDA0002681456900000103
和协方差先验估计P(k|k-1)分别为：
Figure BDA0002681456900000104
Figure BDA0002681456900000105
其中
Figure BDA0002681456900000106
观测更新：
从状态和协方差的先验估计中生成2n+1个第三sigma点集：
Figure BDA0002681456900000111
通过观测函数h将第三sigma点集转换为第四转换sigma点集：
γi(k)＝h(k,χi(k|k-1)),for i＝0...2n (16)
则观测先验估计为：
Figure BDA0002681456900000112
状态和观测之间的协方差为：
Figure BDA0002681456900000113
S30：根据状态先验估计、协方差先验估计、观测先验估计和状态和观测 之间的协方差，构建线性回归方程。
获取观测斜率矩阵，所述观测斜率矩阵可定义为：
H(k)＝(P-1(k|k-1)Pxy(k))T (19)
则观测方程可以近似为
Figure BDA0002681456900000114
构建如下线性回归方程：
Figure BDA0002681456900000115
其中，所述I为n×n的单位矩阵，
Figure BDA0002681456900000116
且r(k)为观测噪 声，线性回归方程的协方差为
Figure BDA0002681456900000117
R(k)为预设观测协方差 矩阵。
S40：根据线性回归方程计算估计误差。
为了更好地表示出估计误差，以便进行后续的误差优化，在计算估计误 差之前，需要将现行长城进行转换，具体转换过程为：首先对线性回归方程 的协方差进行Cholesky：
Figure BDA0002681456900000121
其中，对协方差先验估计P(k|k-1)进行Cholesky分解得到Sp(k|k-1)，对 R(k)进行Cholesky分解得到Sr(k)。
将式(21)两边左乘以S-1(k)，得到
D(k)＝W(k)x(k)+e(k) (23)
其中，
Figure BDA0002681456900000122
e(k)＝S-1(k)ξ(k)。此时，E[e(k)eT(k)]＝I。由此得到目标系统的估计误差为： e(k)＝D(k)-W(k)x(k)。
上述步骤通过对目标系统的重构得到线性回归方程，并得到目标系统的 估计误差e(k)，该误差表示了状态估计与先验估计之间的差异以及实际观测值 与预测观测值之间的差异。接下来对其进行优化以获得最优估计，也即是目 标系统的状态。
S50：采用基于模糊交叉熵构造的损失函数对估计误差进行优化，得到目 标系统的状态估计。
关于模糊交叉熵本实施例的开头部分已经做过介绍，再次不再赘述。基 于模糊交叉熵构造的损失函数为：
Figure BDA0002681456900000131
其中，a为预设值，σ为高斯核函数的核宽，μik为k时刻第i维的模糊隶属 度，ei(k)＝di(k)-wi(k)x(k)，ei(k)是e(k)的第i个元素，di(k)是D(k)的第i个元素， wi(k)是W(k)的第i行元素，L是e(k)的维数，
Figure BDA0002681456900000132
通过对损失函数进行计算得到最优解，该最优解为目标系统的状态估计：
Figure BDA0002681456900000133
上述极值的约束条件为
Figure BDA0002681456900000134
可以用拉格朗日乘数法来求解：
Figure BDA0002681456900000135
最优化的一阶必要条件为
Figure BDA0002681456900000136
Figure BDA0002681456900000137
由式(28)可得
Figure BDA0002681456900000138
将式(29)代入式(27)得
Figure BDA0002681456900000139
因而
Figure BDA0002681456900000141
将上式代入式(29)得
Figure BDA0002681456900000142
用类似的方法可以得到x(k)的值。
令
Figure BDA0002681456900000143
可得
Figure BDA0002681456900000144
因此可以利用定点迭代算法进行求解：
Figure BDA0002681456900000145
将式(32)转化为矩阵形式：
x(k)＝(WT(k)U(k)W(k))-1WT(k)U(k)D(k) (35)
其中
Figure BDA0002681456900000146
且
Figure BDA0002681456900000147
Figure BDA0002681456900000148
将上式(35)进行转化：
Figure BDA0002681456900000149
Figure BDA00026814569000001410
Figure BDA0002681456900000151
由(36)(37)得
Figure BDA0002681456900000152
由矩阵求逆公式(A-BD-1C)-1＝A-1+A-1B(D-CA-1B)-1CA-1得
Figure BDA0002681456900000153
由(36)-(38)可得
Figure BDA0002681456900000154
结合(35)(40)(41)得：
Figure BDA0002681456900000155
其中
Figure BDA0002681456900000156
协方差矩阵为
Figure BDA0002681456900000157
上述实施例通过对目标系统进行无迹变换，以适用于非线性系统，获得 获得状态先验估计、协方差的先验估计、观测先验估计，以及状态和观测之 间的协方差，并根据得到的数据构建线性回归方程，最后通过模糊交叉熵对 误差进行优化，以获得目标系统的状态估计。
上述实施例基于模糊信息理论构建了一种模糊交叉熵，利用模糊隶属度 较好地表示估计误差的不同维度对状态估计的不同影响，然后，采用最大模 糊交叉熵准则优化无迹卡尔曼滤波算法得到的估计误差，使得该状态估计方 法可以较好地处理非线性非高斯系统，提高了无迹卡尔曼滤波算法在非高斯 噪声环境中的鲁棒性，提高无迹卡尔曼滤波算法的估算精度。
此外，模糊交叉熵的性能主要取决于核宽的选择。核宽过小，虽然鲁棒 性会提高，但是收敛速度过慢，耗时长，有时甚至会陷入发散状态；核宽过 大，虽然收敛速度快，但是性能会有所减弱。因此，如何选取合适的核宽是 亟待解决的问题。在实施例中，主要是对两方面的误差进行优化：一方面是 状态估计与先验估计之间的差异；另一方面是实际观测值与预测观测值之间 的差异。由于现实情况中过程噪声和观测噪声是完全不同的，因此应该针对 每一维选取合适的核宽而不是直接统一选取固定的核宽。
对此，本实施例采用自适应调整核宽的方法，将核宽设置为估计误差的 绝对值除以2的平方根和预设的核宽之间的最大值，即：
Figure BDA0002681456900000161
其中，σi是第i维误差的核宽，ei是第i维的误差，σ0为核宽的预设值。
本实施例采用自适应的方法设置核宽，提高了系统状态估计方法的估计 性能。
为了验证本实施例所提出的系统状态估计方法的估计性能，本实施例采 用了均方根误差(Root-mean Square Error,RMSE)作为性能指数，以验证该方法 的估计性能。均方根误差定义如下：
Figure BDA0002681456900000162
其中，M代表蒙特卡洛运行次数，K表示每次蒙特卡洛运行的总步长。
采用本实施例所提出的系统状态估计方法对一下系统进行验证：
实验一：
单变量非平稳增长模型，其状态方程和观测方程如下：
Figure BDA0002681456900000171
Figure BDA0002681456900000172
其中，过程噪声和观测噪声均为混合高斯分布：
qk～0.8N(0,0.1)+0.2N(0,10)
rk～0.8N(0,1)+0.2N(0,400)
同时采用无迹卡尔曼滤波算法、最大交叉熵滤波算法和本实施例所提出 的系统状态估计方法进行仿真，得到的系统状态估计的均方根误差如下表1 所示：
表1
算法 均方根误差
无迹卡尔曼滤波算法 11.8561
最大交叉熵滤波算法 5.2197
系统状态估计方法 5.1885
在本次仿真实验中，设置K＝500，并进行100次蒙特卡洛运行，即M＝100。
通过表1可以知，无迹卡尔曼滤波算法的性能最差，这是因为无迹卡尔 曼滤波算法的抗非高斯噪声能力较弱。最大交叉熵滤波算法与本实施例的基 于模糊交叉熵的系统状态估计方法效果较好，这是由于交叉熵准则能更好的 处理非高斯噪声信号。且系统状态估计方法的性能优于最大交叉熵滤波算法， 说明模糊交叉熵较好地处理了普通交叉熵中存在的问题，提高了估计性能。
实验二：
本实验利用纯方位目标跟踪仿真数据进行验证，采用如下目标跟踪模型：
Figure BDA0002681456900000181
yk＝h(xk)+ek (54)
Figure BDA0002681456900000182
其中，状态向量为
Figure BDA0002681456900000183
xk、yk、zk分别表示目标在k时 刻的位置，
Figure BDA0002681456900000184
分别表示目标k时刻在xk、yk、zk方向上的速度。
过程噪声，vk-1～N(0,Q)，Q＝diag([0.012km2s4 0.012km2s4])。观测噪声ek～(0,R)， R＝diag([0.0012rad2 0.0012rad2])。(Si,x,Si,y,Si,z)，i＝1,2分别表示两个被传感器的位 置。被动传感器观测站1的位置为(0,5km,0)，被动传感器观测站2的位置为 (0,-5km,0)。
在本次仿真实验中，K＝80，并进行100次蒙特卡洛运行，即M＝100。
图2给出了目标的真实运动轨迹，图3给出了各个算法的跟踪结果，其 中，TrueTrajcectory表示真实运动轨迹，MFC-UF表示本实施例所提出的系 统状态估计方法(即基于模糊交叉熵的滤波算法)，MCUF表示最大交叉熵 滤波算法，UKF为无迹卡尔曼滤波算法。图4是各个算法的均方根误差对比 图，其中图4(a)为X方向均方根误差，图4(b)为Y方向均方根误差，图4(c) 为Z方向均方根误差，图4(d)为位置均方根误差。
图4可以看出，本实施例提出的系统状态估计方法的性能优于最大交叉 熵滤波算法和无迹卡尔曼滤波算法。主要原因在于本实施例的系统状态估计 方法采用模糊交叉熵准则优化无迹卡尔曼滤波算法，引入了模糊隶属度处理 不同样本对状态估计的不同影响，提高了状态估计的准确性。
表2不同观测噪声下的均方根误差
观测噪声 UKF MCUF MFC-UF
0.001<sup>2</sup> 0.2816 0.2662 0.2543
0.005<sup>2</sup> 0.9274 1.0433 0.8214
0.01<sup>2</sup> 1.8155 1.3061 1.3609
为了更好地分析所提算法的性能，我们固定过程噪声为0.012，针对不同的 观测噪声进行实验，结果如表2所示。随着观测噪声的提高，所有算法的性 能都会下降。观测噪声较小时，本实施例的系统状态估计方法的性能优于无 迹卡尔曼滤波和最大交叉熵无迹滤波算法。
应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后， 各过程的执行顺序应以其功能和内在逻辑确定，而不应对本申请实施例的实 施过程构成任何限定。
在一实施例中，如图5所示，提供一种系统状态估计装置，用于对目标 系统的状态进行估计，该装置包括：
获取单元10，用于获取目标系统的状态方程和观测方程；
无迹变换单元20，用于根据上一时刻的协方差，对所述状态方程和观测 方程进行无迹变换，得到状态先验估计、协方差先验估计、观测先验估计， 以及状态和观测之间的协方差；
线性方程构建单元30，用于根据所述状态先验估计、协方差先验估计、 观测先验估计和状态和观测之间的协方差，构建线性回归方程；
误差估计单元40，用于根据所述线性回归方程计算估计误差；
优化单元50，用于采用基于模糊交叉熵构造的损失函数对所述估计误差 进行优化，得到所述目标系统的状态。
关于系统状态估计装置的具体限定可以参见上文中对于系统状态估计方 法的限定，在此不再赘述。
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以 上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而 将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划 分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。
以上所述实施例仅用以说明本申请的技术方案，而非对其限制；尽管参 照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解： 其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技 术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱 离本申请各实施例技术方案的精神和范围，均应包含在本申请的保护范围之 内。

Claims (10)

1.一种系统状态估计方法，用于对目标系统的状态进行估计，其特征在于，包括：
获取目标系统的状态方程和观测方程；
根据目标系统的上一时刻的协方差，对所述状态方程和观测方程进行无迹变换，得到状态先验估计、协方差先验估计、观测先验估计，以及状态和观测之间的协方差；
根据所述状态先验估计、协方差先验估计、观测先验估计和状态和观测之间的协方差，构建线性回归方程；
根据所述线性回归方程计算估计误差；
采用基于模糊交叉熵构造的损失函数对所述估计误差进行优化，得到所述目标系统的状态估计。
2.如权利要求1所述的系统状态估计方法，其特征在于，所述状态方程包括状态函数，且所述状态方程为n维状态向量；
所述方法通过以下步骤计算所述状态先验估计和协方差先验估计：
根据所述上一时刻的协方差和所述状态方程获取k-1时刻的多个第一sigma点集，且所述k-1时刻的第一sigma点集为：
Figure FDA0002681456890000011
Figure FDA0002681456890000012
其中，
Figure FDA0002681456890000013
是矩阵(n+λ)P(k-1|k-1)平方根的第i列，P(k-1|k-1)为所述上一时刻的协方差，n为状态维度，λ为复合比例因子，且λ＝α2(n+φ)-n，α和φ均为预设值；
根据所述状态函数，将所述k-1时刻的第一sigma点集转换为k时刻的多个第二sigma点集：
χi*(k|k-1)＝f(k-1,χi(k-1|k-1)),for i＝0...2n
其中，χi*(k|k-1)表示第二sigma点集，f表示状态函数；
根据所述k时刻的多个第二sigma点集，计算所述状态先验估计和协方差先验估计：
Figure FDA0002681456890000021
Figure FDA0002681456890000022
其中，所述
Figure FDA0002681456890000023
为状态先验估计，P(k|k-1)为协方差先验估计，
Figure FDA0002681456890000024
for i＝1...2n，Q(k-1)为预设的状态协方差矩阵。
3.如权利要求2所述的系统状态估计方法，其特征在于，所述观测方程包括观测函数；所述方法通过以下步骤计算所述观测先验估计，以及状态和观测之间的协方差：
根据所述状态先验估计和协方差先验估计，获取所述k时刻的多个第三sigma点集，且所述第三sigma点集为：
Figure FDA0002681456890000025
Figure FDA0002681456890000026
根据所述观测函数，将所述k时刻的第三sigma点集转换为k时刻的第四sigma点集：
γi(k)＝h(k,χi(k|k-1)),for i＝0...2n
其中，γi(k)为第四sigma点集，h()为观测函数；
根据所述k时刻的第四sigma点集，计算所述观测先验估计：
Figure FDA0002681456890000031
根据所述第二sigma点集、状态先验估计、第四sigma点集和观测先验估计计算所述状态和观测之间的协方差：
Figure FDA0002681456890000032
4.如权利要求3所述的系统状态估计方法，其特征在于，所述观测方程包括观测噪声；
通过以下步骤构建所述线性回归方程，包括：
获取观测斜率矩阵，所述观测斜率矩阵为：
H(k)＝(P-1(k|k-1)Pxy(k))T；
根据所述观测斜率矩阵计算线性回归方程，所述线性回归方程为：
Figure FDA0002681456890000033
其中，所述I为n×n的单位矩阵，
Figure FDA0002681456890000034
且r(k)为观测噪声，线性回归方程的协方差为
Figure FDA0002681456890000035
R(k)为预设观测协方差矩阵。
5.如权利要求4所述的系统状态估计方法，其特征在于，所述估计误差为：
e(k)＝D(k)-W(k)x(k)
其中，所述e(k)为估计误差，
Figure FDA0002681456890000041
Figure FDA0002681456890000042
x(k)为状态方程，S(k)通过对所述线性回归方程的协方差进行Cholesky分解得到。
6.如权利要求5所述的系统状态估计方法，其特征在于，所述损失函数为：
Figure FDA0002681456890000043
其中，a为预设值，σ为高斯核函数的核宽，μik为k时刻第i维的模糊隶属度，ei(k)＝di(k)-wi(k)x(k)，ei(k)是e(k)的第i个元素，di(k)是D(k)的第i个元素，wi(k)是W(k)的第i行元素，L是e(k)的维数，
Figure FDA0002681456890000044
根据所述损失函数计算得到所述目标系统的状态估计为：
Figure FDA0002681456890000045
7.如权利要求6所述的系统状态估计方法，其特征在于，所述σ的取值为：
Figure FDA0002681456890000046
其中，σi是第i维误差的核宽，ei是第i维的误差，σ为核宽的预设值。
8.如权利要求6所述的系统状态估计方法，其特征在于，在得到所述目标系统的状态，计算所述目标系统的协方差后验估计：
Figure FDA0002681456890000047
其中，
Figure FDA0002681456890000051
Figure FDA0002681456890000052
Figure FDA0002681456890000053
Sp(k|k-1)为对协方差先验估计P(k|k-1)进行Cholesky分解得到，Sr(k)为对R(k)进行Cholesky分解得到。
9.如权利要求6所述的系统状态估计方法，其特征在于，所述μik的取值为：
Figure FDA0002681456890000054
10.一种系统状态估计装置，用于对目标系统的状态进行估计，其特征在于，包括：
获取单元，用于获取目标系统的状态方程和观测方程；
无迹变换单元，用于根据上一时刻的协方差，对所述状态方程和观测方程进行无迹变换，得到状态先验估计、协方差先验估计、观测先验估计，以及状态和观测之间的协方差；
线性方程构建单元，用于根据所述状态先验估计、协方差先验估计、观测先验估计和状态和观测之间的协方差，构建线性回归方程；
误差估计单元，用于根据所述线性回归方程计算估计误差；
优化单元，用于采用基于模糊交叉熵构造的损失函数对所述估计误差进行优化，得到所述目标系统的状态。
CN202010963665.8A 2020-09-14 2020-09-14 系统状态估计方法、系统状态估计装置 Active CN112115419B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010963665.8A CN112115419B (zh) 2020-09-14 2020-09-14 系统状态估计方法、系统状态估计装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010963665.8A CN112115419B (zh) 2020-09-14 2020-09-14 系统状态估计方法、系统状态估计装置

Publications (2)

Publication Number Publication Date
CN112115419A true CN112115419A (zh) 2020-12-22
CN112115419B CN112115419B (zh) 2024-07-12

Family

ID=73802831

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010963665.8A Active CN112115419B (zh) 2020-09-14 2020-09-14 系统状态估计方法、系统状态估计装置

Country Status (1)

Country Link
CN (1) CN112115419B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112486134A (zh) * 2020-12-14 2021-03-12 中国科学技术大学 一种多对象的采集控制方法、装置及控制设备
CN113625552A (zh) * 2021-08-16 2021-11-09 西南大学 对状态受限非线性系统进行鲁棒状态估计的方法及装置
RU2767463C2 (ru) * 2020-04-09 2022-03-17 Федеральное государственное автономное учреждение "Военный инновационный технополис "ЭРА" Устройство для расчета функций распределения потоков сообщений на основе оценочных данных параметров систем информационного обмена
CN114371232A (zh) * 2021-12-22 2022-04-19 天津国科医工科技发展有限公司 基于卡尔曼滤波算法的色谱滤波方法、装置、介质、系统
CN115587629A (zh) * 2022-12-07 2023-01-10 中国科学院上海高等研究院 协方差膨胀系数估计方法、模型训练方法、存储介质终端
CN115792796A (zh) * 2023-02-13 2023-03-14 鹏城实验室 基于相对观测等效模型的协同定位方法、装置及终端
CN115823951A (zh) * 2023-01-09 2023-03-21 中国兵器装备集团自动化研究所有限公司 一种搜索与跟踪航迹融合方法、装置、设备及存储介质
CN115859039A (zh) * 2023-03-01 2023-03-28 南京信息工程大学 一种车辆状态估计方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050251328A1 (en) * 2004-04-05 2005-11-10 Merwe Rudolph V D Navigation system applications of sigma-point Kalman filters for nonlinear estimation and sensor fusion
CN104202019A (zh) * 2014-08-25 2014-12-10 北京理工大学 带有未知过程噪声协方差阵递推估计的卡尔曼滤波方法
CN105378496A (zh) * 2013-09-05 2016-03-02 日本康奈可株式会社 估计装置以及估计方法
US20200132775A1 (en) * 2017-06-14 2020-04-30 Mitsubishi Electric Corporation State estimation device

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050251328A1 (en) * 2004-04-05 2005-11-10 Merwe Rudolph V D Navigation system applications of sigma-point Kalman filters for nonlinear estimation and sensor fusion
CN105378496A (zh) * 2013-09-05 2016-03-02 日本康奈可株式会社 估计装置以及估计方法
CN104202019A (zh) * 2014-08-25 2014-12-10 北京理工大学 带有未知过程噪声协方差阵递推估计的卡尔曼滤波方法
US20200132775A1 (en) * 2017-06-14 2020-04-30 Mitsubishi Electric Corporation State estimation device

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
张开元;吉兴全;于永进;: "采用ARUKF算法的电力系统动态状态估计", 中国科技论文, no. 11, pages 23 - 28 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2767463C2 (ru) * 2020-04-09 2022-03-17 Федеральное государственное автономное учреждение "Военный инновационный технополис "ЭРА" Устройство для расчета функций распределения потоков сообщений на основе оценочных данных параметров систем информационного обмена
CN112486134A (zh) * 2020-12-14 2021-03-12 中国科学技术大学 一种多对象的采集控制方法、装置及控制设备
CN112486134B (zh) * 2020-12-14 2022-04-19 中国科学技术大学 一种多对象的采集控制方法、装置及控制设备
CN113625552A (zh) * 2021-08-16 2021-11-09 西南大学 对状态受限非线性系统进行鲁棒状态估计的方法及装置
CN114371232A (zh) * 2021-12-22 2022-04-19 天津国科医工科技发展有限公司 基于卡尔曼滤波算法的色谱滤波方法、装置、介质、系统
CN114371232B (zh) * 2021-12-22 2024-03-22 天津国科医工科技发展有限公司 基于卡尔曼滤波算法的色谱滤波方法、装置、介质、系统
CN115587629A (zh) * 2022-12-07 2023-01-10 中国科学院上海高等研究院 协方差膨胀系数估计方法、模型训练方法、存储介质终端
CN115587629B (zh) * 2022-12-07 2023-04-07 中国科学院上海高等研究院 协方差膨胀系数估计方法、模型训练方法、存储介质终端
CN115823951A (zh) * 2023-01-09 2023-03-21 中国兵器装备集团自动化研究所有限公司 一种搜索与跟踪航迹融合方法、装置、设备及存储介质
CN115823951B (zh) * 2023-01-09 2023-04-18 中国兵器装备集团自动化研究所有限公司 一种搜索与跟踪航迹融合方法、装置、设备及存储介质
CN115792796A (zh) * 2023-02-13 2023-03-14 鹏城实验室 基于相对观测等效模型的协同定位方法、装置及终端
CN115859039A (zh) * 2023-03-01 2023-03-28 南京信息工程大学 一种车辆状态估计方法

Also Published As

Publication number Publication date
CN112115419B (zh) 2024-07-12

Similar Documents

Publication Publication Date Title
CN112115419B (zh) 系统状态估计方法、系统状态估计装置
CN108897286B (zh) 一种基于分散式非线性动态关系模型的故障检测方法
CN106599368B (zh) 基于改进粒子提议分布和自适应粒子重采样的FastSLAM方法
Dovžan et al. Recursive fuzzy c-means clustering for recursive fuzzy identification of time-varying processes
CN104376581B (zh) 一种采用自适应重采样的高斯混合无迹粒子滤波算法
Tzeng Design of fuzzy wavelet neural networks using the GA approach for function approximation and system identification
CN106487358A (zh) 一种基于统计线性回归的最大相关熵容积卡尔曼滤波方法
CN103902812B (zh) 一种粒子滤波方法、装置及目标跟踪方法、装置
CN110377942B (zh) 一种基于有限高斯混合模型的多模型时空建模方法
CN106933106A (zh) 一种基于模糊控制多模型算法的目标跟踪方法
CN105205313A (zh) 模糊高斯和粒子滤波方法、装置及目标跟踪方法、装置
CN109284677A (zh) 一种贝叶斯滤波目标跟踪算法
CN103902819A (zh) 基于变分滤波的粒子优化概率假设密度多目标跟踪方法
CN105354860A (zh) 基于箱粒子滤波的扩展目标CBMeMBer跟踪方法
CN106067034B (zh) 一种基于高维矩阵特征根的配电网负荷曲线聚类方法
CN107861492A (zh) 一种基于裕度统计量的广义非负矩阵分解故障监测方法
CN114819054B (zh) 一种基于物理信息神经网络的电力电子系统状态监测方法
CN114449452A (zh) 一种异构设备室内定位算法
CN108871365B (zh) 一种航向约束下的状态估计方法及系统
CN113452349B (zh) 一种基于贝叶斯序贯重要性积分的卡尔曼滤波方法
CN113221311A (zh) 一种大气边界层风速的不确定性量化方法
Du et al. A novel locally regularized automatic construction method for RBF neural models
CN113190960B (zh) 一种基于非等维状态混合估计的并行imm机动目标跟踪方法
CN112865748A (zh) 基于递归最小二乘的在线分布式多任务图滤波器构建方法
El Farissi et al. Application of neuro-fuzzy in the recognition of control chart patterns

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