CN111444475A - 一种应用在飞行测试数据分析的容错ckf滤波融合方法 - Google Patents

一种应用在飞行测试数据分析的容错ckf滤波融合方法 Download PDF

Info

Publication number
CN111444475A
CN111444475A CN202010212538.4A CN202010212538A CN111444475A CN 111444475 A CN111444475 A CN 111444475A CN 202010212538 A CN202010212538 A CN 202010212538A CN 111444475 A CN111444475 A CN 111444475A
Authority
CN
China
Prior art keywords
filter
sub
noise
fault
measurement
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
CN202010212538.4A
Other languages
English (en)
Other versions
CN111444475B (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 Research Institute of Guangdong Ocean University
Original Assignee
Shenzhen Research Institute of Guangdong Ocean 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 Research Institute of Guangdong Ocean University filed Critical Shenzhen Research Institute of Guangdong Ocean University
Priority to CN202010212538.4A priority Critical patent/CN111444475B/zh
Publication of CN111444475A publication Critical patent/CN111444475A/zh
Application granted granted Critical
Publication of CN111444475B publication Critical patent/CN111444475B/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
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64F—GROUND OR AIRCRAFT-CARRIER-DECK INSTALLATIONS SPECIALLY ADAPTED FOR USE IN CONNECTION WITH AIRCRAFT; DESIGNING, MANUFACTURING, ASSEMBLING, CLEANING, MAINTAINING OR REPAIRING AIRCRAFT, NOT OTHERWISE PROVIDED FOR; HANDLING, TRANSPORTING, TESTING OR INSPECTING AIRCRAFT COMPONENTS, NOT OTHERWISE PROVIDED FOR
    • B64F5/00—Designing, manufacturing, assembling, cleaning, maintaining or repairing aircraft, not otherwise provided for; Handling, transporting, testing or inspecting aircraft components, not otherwise provided for
    • B64F5/60—Testing or inspecting aircraft components or systems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M99/00—Subject matter not provided for in other groups of this subclass
    • G01M99/005—Testing of complete machines, e.g. washing-machines or mobile phones

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Probability & Statistics with Applications (AREA)
  • Software Systems (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Algebra (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Databases & Information Systems (AREA)
  • Operations Research (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Manufacturing & Machinery (AREA)
  • Transportation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Feedback Control In General (AREA)

Abstract

本发明涉及一种应用在飞行测试数据分析的容错CKF滤波融合方法。本发明大体包括三部分内容：第一部分，根据实际运动目标进行系统建模；第二部分，根据复杂工程环境下量测系统中存在乘性噪声相关以及故障问题设计了乘性噪声相关的容错CKF滤波器；第三部分，根据无重置式联邦滤波器的不足，提出了改进的滤波融合方法。本发明既能处理乘性噪声相关的滤波估计问题，又能应对系统发生故障的情形，同时也能处理所有子滤波器均发生故障的极端情况，极大地提高了系统的容错性，保证了飞行测试数据的精度。

Description

一种应用在飞行测试数据分析的容错CKF滤波融合方法
技术领域
本发明涉及一种应用在飞行测试数据分析的容错CKF滤波融合方法，属于目标跟踪领域。
背景技术
飞行实验测试在评估飞行器质量和性能中有着非常重要的意义，利用测量过程中采集到的实际飞行数据对飞行器的状态跟踪轨迹分析是评估飞行器性能的重要手段，高精度的目标跟踪数据对于评估和分析飞行器运行过程中的质量和稳定性极为重要。
飞行器和测量设备在运行过程中会受到复杂的环境、设备作用距离和通信电磁波的干扰，使得获得的实际飞行测试数据被噪声污染而不能直接被使用。为了解决噪声干扰问题，1960年提出的卡尔曼滤波理论，在飞行器飞行测试数据分析方面得到了重要应用。因此，开展高精度的飞行器飞行测试的数据分析研究，对提高飞行试验飞行器性能评估的准确性和稳定性具有重要意义。
由于目标跟踪系统多为非线性，卡尔曼滤波只能适用于线性系统。因此，大量非线性滤波方法得以提出，其中容积卡尔曼滤波(CKF)因具有较高的滤波估计精度而被广泛使用。随着实际工程环境日趋复杂，传感器量测会受到乘性噪声和野值的干扰，导致滤波精度降低。
随着科学技术的迅猛发展，人们对滤波性能的要求越来越高，多传感器信息融合技术受到广大科研人员和工程技术人员的欢迎。融合结构常分为两种：集中式融合和分布式融合。其中，分布式融合中的无重置式联邦滤波器容错性较好，且具有较快的计算速度。但是传统的无重置式联邦滤波器也面临着两个问题：一，对故障子系统隔离之后，由于该故障子系统无法获得融合估计的结果，因此不能再进行故障检测，将一直被隔离造成传感器资源浪费，同时也会对融合稳定性造成影响。二，对于出现所有子系统均发生故障的情况，融合估计将无法继续进行。
发明内容
为了应对上面提到的系统出现故障的信息融合问题，本发明设计了容错处理模块代替故障隔离模块，并对故障检测方法做出了改进，且对子滤波器提出乘性噪声相关的容错CKF滤波算法，有效实现了对飞行器状态的实时估计。
本发明大体包括三部分内容。第一部分根据实际运动目标进行系统建模；第二部分，根据量测子系统中出现故障以及出现乘性噪声相关的情况设计了乘性噪声相关的容错CKF滤波算法；第三部分，设计出改进的滤波融合方法，用于对飞行器飞行测试数据进行分析，得到飞行器实时运动状态。
本发明包括以下步骤：
步骤1.系统建模，假设系统有N个传感器，每个传感器构成一个子滤波器，且子滤波器与主滤波器的状态变量相同，考虑如下具有乘性噪声相关的离散时间非线性系统，其状态方程和第i(i＝1,2,…,N)个传感器的量测方程分别如下：
Figure BDA0002423308600000021
式中，
Figure BDA0002423308600000022
是k时刻的系统状态向量，其是由x方向位移和速度以及y方向位移和速度构成，f为已知非线性过程函数，过程噪声wk-1是零均值方差为Qk-1的高斯白噪声向量，
Figure BDA0002423308600000023
是k时刻第i个子滤波器的量测向量，hi为第i个子滤波器非线性量测函数，Ai,k＝diag{1+ui,k,…,1+ui,k}＝(1+ui,k)I，ui,k和vi,k分别是第i个子滤波器k时刻乘性和加性高斯白噪声，且具有相关性，vi,k＝[v1,k,…,vm,k]T，且E{vi,k}＝μi,v1＝μi,v[1,…,1]T，wk-1和vi,k互不相关，σi为随机向量，用来描述野值大小，当ρk＝0时表示无故障，当ρk＝1时，表示发生故障。
Figure BDA0002423308600000024
式中，δkj为Kronecher-δ函数，μi,u和
Figure BDA0002423308600000025
分别为第i个子滤波器乘性量测噪声均值和方差，μi,v和
Figure BDA0002423308600000026
分别为第i个子滤波器加性量测噪声均值和方差，di,k为第i个子滤波器噪声相关系数。
步骤2.模型转换，将量测方程表示为非线性量测和虚拟量测噪声的总和如下：
Figure BDA0002423308600000031
式中，
Figure BDA0002423308600000032
为第i个子滤波器k时刻的虚拟量测噪声，其相对应的均值为
Figure BDA0002423308600000033
方差为
Figure BDA0002423308600000034
步骤3.给出子滤波器的乘性噪声相关的容错CKF滤波算法，具体如下：
步骤3.1时间更新阶段，已知上一时刻子滤波器的状态估计值
Figure BDA0002423308600000035
和估计误差协方差矩阵Pi,k-1，该滤波算法的时间更新如下：
Figure BDA0002423308600000036
Figure BDA0002423308600000037
Figure BDA0002423308600000038
Figure BDA0002423308600000039
Figure BDA00024233086000000310
式中，Sk-1为由估计误差协方差矩阵Pi,k-1经过cholesky分解得到，ξj为提前确定的cubature点，
Figure BDA00024233086000000311
和
Figure BDA00024233086000000312
都为经计算和传播后的第j个cubature点，
Figure BDA00024233086000000313
为状态预测值，Pk|k-1为预测误差协方差矩阵。
步骤3.2给出E{hi(xk)}、Var{hi(xk)}和Cov{xk,hi(xk)}的计算方法：
Figure BDA00024233086000000314
Figure BDA00024233086000000315
Figure BDA00024233086000000316
Figure BDA00024233086000000317
Figure BDA0002423308600000041
Figure BDA0002423308600000042
式中，Sk|k-1为由预测误差协方差矩阵Pk+1|k经过cholesky分解得到，
Figure BDA0002423308600000043
和
Figure BDA0002423308600000044
都为k时刻第j个cubature点，E{hi(xk)}和Var{hi(xk)}分别为hi(xk)的均值和方差，Cov{xk,hi(xk)}为xk和hi(xk)的互协方差。
步骤3.3给出虚拟量测噪声均值
Figure BDA0002423308600000045
和虚拟量测噪声方差
Figure BDA0002423308600000046
的计算方法：
Figure BDA0002423308600000047
Figure BDA0002423308600000048
步骤3.4给出量测预测值
Figure BDA0002423308600000049
新息协方差矩阵Pzz,k|k-1和互协方差矩阵Pxz,k|k-1的计算方法：
Figure BDA00024233086000000410
Figure BDA00024233086000000411
Pxz,k|k-1＝(1+μi,u)Cov{xk,hi(xk)} (19)
步骤3.5给出故障检测方法以及容错策略，根据子滤波器当前滤波新息与理论新息协方差的不一致程度判断子滤波器是否发生故障，其步骤如下：
Figure BDA00024233086000000412
Figure BDA00024233086000000413
Figure BDA00024233086000000414
Figure BDA00024233086000000415
式中，εi,k为滤波新息，αi,k为故障检测函数，Ti,D为故障检测阈值，可提前根据预警率设定，λi,k为调节因子。
将调节因子λi,k引入到滤波增益矩阵Kk中，其公式如下：
Kk＝λi,kPxz,k|k-1(Pzz,k|k-1)-1 (24)
步骤3.6给出状态估计值
Figure BDA0002423308600000051
和估计误差协方差Pi,k的计算方法：
Figure BDA0002423308600000052
Pi,k＝Pk|k-1-KkPzz,k|k-1(Kk)T (26)
步骤4.给出滤波融合算法，具体算法如下：
步骤4.1信息分配，只在初始时刻进行一次信息分配，具体如下：
Figure BDA0002423308600000053
Figure BDA0002423308600000054
Figure BDA0002423308600000055
式中，
Figure BDA0002423308600000056
为全局状态估计值，Pg,0为其对应的估计误差协方差矩阵，Qg,0为第i个子滤波器的噪声方差矩阵，信息分配系数βi为：
Figure BDA0002423308600000057
步骤4.2给出时间更新和量测更新，按照步骤3对每一个子滤波器各自独立地进行时间更新和量测更新的各个环节，得到各个子滤波器的状态估计值
Figure BDA0002423308600000058
和估计误差协方差矩阵Pi,k。
步骤4.3给出主滤波器的融合算法：
Pg,k＝[(P1,k)-1+(P2,k)-1+…+(PN,k)-1]-1 (31)
Figure BDA0002423308600000059
步骤4.4给出信息反馈策略：
主滤波器完成当前时刻的信息融合后，将全局融合信息反馈给发生故障的子滤波器，直到该子滤波器恢复正常以后，则不再进行信息反馈；
Figure BDA00024233086000000510
Pi,k＝Pg,k (34)
对于没发生故障的子滤波器，主滤波器不对其进行信息反馈。
本发明的有益效果：本发明能够处理系统中出现故障以及乘性量测噪声与加性量测噪声相关情形，且有效解决所有传感器均发生故障的极端情况，同时也避免了将故障子滤波器隔离造成的资源浪费问题，将其用于飞行测试数据分析中，可以实时估计出飞行器状态，有效得到高精度的飞行测试数据。
附图说明：
图1：本发明的乘性噪声相关的容错CKF滤波算法流程图。
图2：本发明的容错CKF滤波融合算法结构图。
具体实施方法
本发明提出一种应用在飞行测试数据分析的容错CKF滤波融合方法。本发明首先根据实际目标的运动状态建立模型，其次给出子滤波器的乘性噪声相关的容错CKF滤波算法的步骤，最后给出主滤波器融合算法和信息反馈策略。乘性噪声相关的容错CKF滤波算法流程图如图1所示，容错CKF滤波融合算法结构图如图2所示，包括以下几个步骤：
步骤1.系统建模，假设系统有N个传感器，每个传感器构成一个滤波子滤波器，且子滤波器与主滤波器的状态变量相同，考虑如下具有乘性噪声相关的离散时间非线性系统，其状态方程和第i(i＝1,2,…,N)个子滤波器的量测方程分别如下：
Figure BDA0002423308600000061
式中，
Figure BDA0002423308600000062
是k时刻的系统状态向量，其是由x方向位移和速度以及y方向位移和速度构成，f为已知非线性过程函数，过程噪声wk-1是零均值方差为Qk-1的高斯白噪声向量，
Figure BDA0002423308600000063
是k时刻第i个子滤波器的量测向量，hi为第i个子滤波器非线性量测函数，Ai,k＝diag{1+ui,k,…,1+ui,k}＝(1+ui,k)I，ui,k和vi,k分别是第i个子滤波器k时刻乘性和加性高斯白噪声，且具有相关性，vi,k＝[v1,k,…,vm,k]T，且E{vi,k}＝μi,v1＝μi,v[1,…,1]T，wk-1和vi,k互不相关，σi为随机向量，用来描述野值大小，当ρk＝0时表示无故障，当ρk＝1时，表示发生故障。
Figure BDA0002423308600000071
式中，δkj为Kronecher-δ函数，μi,u和
Figure BDA0002423308600000072
分别为第i个子滤波器乘性量测噪声均值和方差，μi,v和
Figure BDA0002423308600000073
分别为第i个子滤波器加性量测噪声均值和方差，di,k为第i个子滤波器噪声相关系数。
步骤2.模型转换，将量测方程表示为非线性量测和虚拟量测噪声的总和如下：
Figure BDA0002423308600000074
式中，
Figure BDA0002423308600000075
为第i个子滤波器k时刻的虚拟量测噪声，其相对应的均值为
Figure BDA0002423308600000076
方差为
Figure BDA0002423308600000077
步骤3.给出子滤波器的乘性噪声相关的容错CKF滤波算法，每一个子滤波器各自独立地进行时间更新和量测更新的各个环节，得到各个子滤波器的状态估计值
Figure BDA0002423308600000078
和估计误差协方差矩阵Pi,k，具体如下：
步骤3.1时间更新阶段，已知上一时刻子滤波器的状态估计值
Figure BDA0002423308600000079
和估计误差协方差矩阵Pi,k-1，该滤波算法的时间更新如下：
Figure BDA00024233086000000710
Figure BDA00024233086000000711
Figure BDA00024233086000000712
Figure BDA00024233086000000713
Figure BDA00024233086000000714
式中，Sk-1为由估计误差协方差矩阵Pi,k-1经过cholesky分解得到，ξj为提前确定的cubature点，
Figure BDA0002423308600000081
和
Figure BDA0002423308600000082
都为经计算和传播后的第j个cubature点，
Figure BDA0002423308600000083
为状态预测值，Pk|k-1为预测误差协方差矩阵。
步骤3.2给出E{hi(xk)}、Var{hi(xk)}和Cov{xk,hi(xk)}的计算方法：
Figure BDA0002423308600000084
Figure BDA0002423308600000085
Figure BDA0002423308600000086
Figure BDA0002423308600000087
Figure BDA0002423308600000088
Figure BDA0002423308600000089
式中，Sk|k-1为由预测误差协方差矩阵Pk+1|k经过cholesky分解得到，
Figure BDA00024233086000000810
和
Figure BDA00024233086000000811
都为k时刻第j个cubature点，E{hi(xk)}和Var{hi(xk)}分别为hi(xk)的均值和方差，Cov{xk,hi(xk)}为xk和hi(xk)的互协方差。
步骤3.3给出虚拟量测噪声均值
Figure BDA00024233086000000812
和虚拟量测噪声方差
Figure BDA00024233086000000813
的计算方法：
Figure BDA00024233086000000814
Figure BDA00024233086000000815
步骤3.4给出量测预测值
Figure BDA00024233086000000816
新息协方差矩阵Pzz,k|k-1和互协方差矩阵Pxz,k|k-1的计算方法：
Figure BDA00024233086000000817
Figure BDA00024233086000000818
Pxz,k|k-1＝(1+μi,u)Cov{xk,hi(xk)} (19)
步骤3.5给出故障检测方法以及容错策略，根据子滤波器当前滤波新息与理论新息协方差的不一致程度判断子滤波器是否发生故障，其步骤如下：
Figure BDA0002423308600000091
Figure BDA0002423308600000092
Figure BDA0002423308600000093
Figure BDA0002423308600000094
式中，εi,k为滤波新息，ai,k为故障检测函数，Ti,D为故障检测阈值，可提前根据预警率设定，λi,k为调节因子。
将调节因子λi,k引入到滤波增益矩阵Kk中，其公式如下：
Kk＝λi,kPxz,k|k-1(Pzz,k|k-1)-1 (24)
步骤3.6给出状态估计值
Figure BDA0002423308600000095
和估计误差协方差Pi,k的计算方法：
Figure BDA0002423308600000096
Pi,k＝Pk|k-1-KkPzz,k|k-1(Kk)T (26)
步骤4.给出滤波融合算法，具体算法如下：
步骤4.1.给出初始化信息分配过程，只在初始时刻进行一次信息分配，具体如下：
Figure BDA0002423308600000097
Figure BDA0002423308600000098
Figure BDA0002423308600000099
式中，
Figure BDA00024233086000000910
为全局状态估计值，Pg,0为其对应的协方差矩阵，Qg,0为主滤波器的噪声方差矩阵，信息分配系数βi为：
Figure BDA0002423308600000101
步骤4.2给出时间更新和量测更新，按照步骤3对每一个子滤波器各自独立地进行时间更新和量测更新的各个环节，得到各个子滤波器的状态估计值
Figure BDA0002423308600000102
和估计误差协方差矩阵Pi,k。
步骤4.3给出主滤波器的融合算法，主滤波器将各子滤波器的信息进行融合，得到全局状态估计值和对应的估计误差协方差矩阵：
Pg,k＝[(P1,k)-1+(P2,k)-1+…+(PN,k)-1]-1 (31)
Figure BDA0002423308600000103
步骤4.4给出信息反馈策略：
主滤波器完成当前时刻的信息融合后，将全局融合信息反馈给发生故障的子滤波器，直到该子滤波器恢复正常以后，则不再进行信息反馈；
Figure BDA0002423308600000104
Pi,k＝Pg,k (34)
对于没发生故障的子滤波器，主滤波器不对其进行信息反馈。
按照步骤4的要求进行循环迭代，得到飞行器的实时状态估计值。
本发明所述容错CKF滤波融合方法可以处理乘性噪声相关的滤波估计问题，当子滤波器发生故障时，通过对子滤波器进行容错处理，不但保证了子滤波器的鲁棒性，而且提高了全局融合估计精度，即使所有子滤波器均发生故障也能保证滤波估计精度，有效提高了复杂工程环境下飞行测试数据的精度。

Claims (1)

1.一种应用在飞行测试数据分析的容错CKF滤波融合方法，其特征在于该方法包括以下步骤：
步骤1.系统建模
假设系统有N个传感器，每个传感器构成一个子滤波器，且子滤波器与主滤波器的状态变量相同，考虑如下具有乘性噪声相关的离散时间非线性系统，其状态方程和第i个传感器的量测方程分别如下：
Figure FDA0002423308590000011
式中，
Figure FDA0002423308590000012
是k时刻的系统状态向量，其是由x方向位移和速度以及y方向位移和速度构成，f为已知非线性过程函数，过程噪声wk-1是零均值方差为Qk-1的高斯白噪声向量；
Figure FDA0002423308590000013
是k时刻第i个子滤波器的量测向量，i＝1,2,…,N，hi为第i个子滤波器非线性量测函数，Ai,k＝diag{1+ui,k,…,1+ui,k}＝(1+ui,k)I，ui,k和vi,k分别是第i个子滤波器k时刻乘性和加性高斯白噪声，且具有相关性，vi,k＝[v1,k,…,vm,k]T，且E{vi,k}＝μi,v1＝μi,v[1,…,1]T，wk-1和vi,k互不相关，σi为随机向量，用来描述野值大小，当ρk＝0时表示无故障，当ρk＝1时，表示发生故障；
Figure FDA0002423308590000014
式中，δkj为Kronecher-δ函数，μi,u和
Figure FDA0002423308590000015
分别为第i个子滤波器乘性量测噪声均值和方差，μi,v和
Figure FDA0002423308590000016
分别为第i个子滤波器加性量测噪声均值和方差，di,k为第i个子系统噪声相关系数；
步骤2.模型转换
将量测方程表示为非线性量测和虚拟量测噪声的总和：
Figure FDA0002423308590000021
式中，
Figure FDA0002423308590000022
Figure FDA0002423308590000023
为第i个子滤波器k时刻的虚拟量测噪声，其相对应的均值为
Figure FDA0002423308590000024
方差为
Figure FDA0002423308590000025
步骤3.给出子滤波器的乘性噪声相关的容错CKF滤波算法，具体如下：
步骤3.1时间更新阶段，已知上一时刻子滤波器的状态估计值
Figure FDA0002423308590000026
和估计误差协方差矩阵Pi,k-1，则时间更新如下：
Figure FDA0002423308590000027
Figure FDA0002423308590000028
Figure FDA0002423308590000029
Figure FDA00024233085900000210
Figure FDA00024233085900000211
式中，Sk-1为由估计误差协方差矩阵Pi,k-1经过cholesky分解得到，ξj为提前确定的cubature点，
Figure FDA00024233085900000212
和
Figure FDA00024233085900000213
为经计算和传播后的第j个cubature点，
Figure FDA00024233085900000214
为状态预测值，Pk|k-1为预测误差协方差矩阵；
步骤3.2计算E{hi(xk)}、Var{hi(xk)}和Cov{xk,hi(xk)}：
Figure FDA00024233085900000215
Figure FDA00024233085900000216
Figure FDA00024233085900000217
Figure FDA00024233085900000218
Figure FDA00024233085900000219
Figure FDA0002423308590000031
式中，Sk|k-1为由预测误差协方差矩阵Pk+1|k经过cholesky分解得到，
Figure FDA0002423308590000032
和
Figure FDA0002423308590000033
都为k时刻第j个cubature点，E{hi(xk)}和Var{hi(xk)}分别为hi(xk)的均值和方差，Cov{xk,hi(xk)}为xk和hi(xk)的互协方差；
步骤3.3计算虚拟量测噪声均值
Figure FDA0002423308590000034
和虚拟量测噪声方差
Figure FDA0002423308590000035
Figure FDA0002423308590000036
Figure FDA0002423308590000037
步骤3.4计算量测预测值
Figure FDA0002423308590000038
新息协方差矩阵Pzz,k|k-1和互协方差矩阵Pxz,k|k-1的：
Figure FDA0002423308590000039
Figure FDA00024233085900000310
Pxz,k|k-1＝(1+μi,u)Cov{xk,hi(xk)} (19)
步骤3.5给出故障检测以及容错策略，根据子滤波器当前滤波新息与理论新息协方差的不一致程度判断子滤波器是否发生故障：
Figure FDA00024233085900000311
Figure FDA00024233085900000312
Figure FDA00024233085900000313
Figure FDA00024233085900000314
式中，εi,k为滤波新息，αi,k为故障检测函数，Ti,D为故障检测阈值，λi,k为调节因子；
将调节因子λi,k引入到滤波增益矩阵Kk中，其公式如下：
Kk＝λi,kPxz,k|k-1(Pzz,k|k-1)-1 (24)
步骤3.6计算状态估计值
Figure FDA0002423308590000041
和估计误差协方差Pi,k：
Figure FDA0002423308590000042
Pi,k＝Pk|k-1-KkPzz,k|k-1(Kk)T (26)
步骤4.给出滤波融合算法：
步骤4.1信息分配，只在初始时刻进行一次信息分配，具体如下：
Figure FDA0002423308590000043
Figure FDA0002423308590000044
Figure FDA0002423308590000045
式中，
Figure FDA0002423308590000046
为全局状态估计值，Pg,0为其对应的估计误差协方差矩阵，Qi,0为第i个子滤波器的噪声方差矩阵，Qg,0为主滤波器的噪声方差矩阵，信息分配系数βi为：
Figure FDA0002423308590000047
步骤4.2给出时间更新和量测更新，按照步骤3对每一个子滤波器各自独立地进行时间更新和量测更新的各个环节，得到各个子滤波器的状态估计值
Figure FDA00024233085900000410
和估计误差协方差矩阵Pi,k；
步骤4.3给出主滤波器的融合算法：
Pg,k＝[(P1,k)-1+(P2,k)-1+…+(PN,k)-1]-1 (31)
Figure FDA0002423308590000048
步骤4.4给出信息反馈策略：
主滤波器完成当前时刻的信息融合后，将全局融合信息反馈给发生故障的子滤波器，直到该子滤波器恢复正常以后，则不再进行信息反馈；
Figure FDA0002423308590000049
Pi,k＝Pg,k (34)
对于没发生故障的子滤波器，主滤波器不对其进行信息反馈。
CN202010212538.4A 2020-03-24 2020-03-24 一种应用在飞行测试数据分析的容错ckf滤波融合方法 Active CN111444475B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010212538.4A CN111444475B (zh) 2020-03-24 2020-03-24 一种应用在飞行测试数据分析的容错ckf滤波融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010212538.4A CN111444475B (zh) 2020-03-24 2020-03-24 一种应用在飞行测试数据分析的容错ckf滤波融合方法

Publications (2)

Publication Number Publication Date
CN111444475A true CN111444475A (zh) 2020-07-24
CN111444475B CN111444475B (zh) 2023-07-14

Family

ID=71629444

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010212538.4A Active CN111444475B (zh) 2020-03-24 2020-03-24 一种应用在飞行测试数据分析的容错ckf滤波融合方法

Country Status (1)

Country Link
CN (1) CN111444475B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113691237A (zh) * 2021-07-27 2021-11-23 浙江工商大学 一种加权融合鲁棒滤波方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4816913A (en) * 1987-11-16 1989-03-28 Technology, Inc., 64 Pixel interpolation circuitry as for a video signal processor
CN109470267A (zh) * 2018-11-02 2019-03-15 佛山科学技术学院 一种卫星姿态滤波方法
CN109829938A (zh) * 2019-01-28 2019-05-31 杭州电子科技大学 一种应用在目标跟踪的自适应容错容积卡尔曼滤波方法
CN110032812A (zh) * 2019-04-18 2019-07-19 河海大学 一种基于自适应容积卡尔曼滤波的动态状态估计方法
CN110196443A (zh) * 2019-06-06 2019-09-03 中国人民解放军战略支援部队信息工程大学 一种飞行器的容错组合导航方法及系统

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4816913A (en) * 1987-11-16 1989-03-28 Technology, Inc., 64 Pixel interpolation circuitry as for a video signal processor
CN109470267A (zh) * 2018-11-02 2019-03-15 佛山科学技术学院 一种卫星姿态滤波方法
CN109829938A (zh) * 2019-01-28 2019-05-31 杭州电子科技大学 一种应用在目标跟踪的自适应容错容积卡尔曼滤波方法
CN110032812A (zh) * 2019-04-18 2019-07-19 河海大学 一种基于自适应容积卡尔曼滤波的动态状态估计方法
CN110196443A (zh) * 2019-06-06 2019-09-03 中国人民解放军战略支援部队信息工程大学 一种飞行器的容错组合导航方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
QUANBO GE: "Cubature information filters with correlated noises and their applications in decentralized fusion", 《SIGNAL PROCESSING》 *
XIAOQING HU等: "Generalised Kalman filter tracking with multiplicative measurement noise in a wireless sensor network", 《IET SIGNAL PROCESSING》 *
雷浩然: "弹载SINS/GNSS组合导航系统研究", 《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113691237A (zh) * 2021-07-27 2021-11-23 浙江工商大学 一种加权融合鲁棒滤波方法
CN113691237B (zh) * 2021-07-27 2024-01-02 浙江工商大学 一种加权融合鲁棒滤波方法

Also Published As

Publication number Publication date
CN111444475B (zh) 2023-07-14

Similar Documents

Publication Publication Date Title
CN109829938B (zh) 一种应用在目标跟踪的自适应容错容积卡尔曼滤波方法
Xiong et al. Detection of satellite attitude sensor faults using the UKF
Brown et al. Efficient semiparametric estimation of expectations
US7180443B1 (en) Reduced state estimator for systems with physically bounded parameters
CN111178385A (zh) 一种鲁棒在线多传感器融合的目标跟踪方法
CN109709934B (zh) 一种飞行控制系统故障诊断冗余设计方法
CN110823217A (zh) 一种基于自适应联邦强跟踪滤波的组合导航容错方法
Baghernezhad et al. Computationally intelligent strategies for robust fault detection, isolation, and identification of mobile robots
Zio et al. Failure and reliability predictions by infinite impulse response locally recurrent neural networks
CN108871365B (zh) 一种航向约束下的状态估计方法及系统
Tehrani et al. Dynamic neural network-based estimator for fault diagnosis in reaction wheel actuator of satellite attitude control system
Wang et al. Uncertainty evaluation and optimization of INS installation measurement using Monte Carlo Method
CN111444475A (zh) 一种应用在飞行测试数据分析的容错ckf滤波融合方法
Nie et al. A robust unscented Kalman filter for intermittent and featureless aircraft sensor faults
Chen et al. A New Particle Predictor for Fault Prediction of Nonlinear Time‐varying Systems
CN109324591B (zh) 基于鲁棒估计器的随机混杂系统的故障诊断装置及方法
CN112034713B (zh) 非理想网络环境下移动目标的最优状态估计方法及系统
Chen et al. Multivariate simulation assessment for virtual metrology
Dubois et al. Performance evaluation of a moving horizon estimator for multi-rate sensor fusion with time-delayed measurements
CN116558406A (zh) 基于状态域的gnss-加速度计集成桥梁变形监测突变故障探测方法
CN115438728A (zh) 一种基于鲁棒滤波的目标融合方法及系统
CN105549003A (zh) 一种汽车雷达目标跟踪方法
Xiong et al. Unscented Kalman filter for fault detection
Chen et al. Sequential inverse covariance intersection fusion kalman filter for multi-sensor systems with packet dropouts and multiplicative noise
Wang et al. Constrained two-stage Kalman filter for real-time state estimation of systems with time-varying measurement noise covariance

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