CN108614804B - 基于信噪比检验的正则化卡尔曼滤波方法 - Google Patents

基于信噪比检验的正则化卡尔曼滤波方法 Download PDF

Info

Publication number
CN108614804B
CN108614804B CN201810379450.4A CN201810379450A CN108614804B CN 108614804 B CN108614804 B CN 108614804B CN 201810379450 A CN201810379450 A CN 201810379450A CN 108614804 B CN108614804 B CN 108614804B
Authority
CN
China
Prior art keywords
state
parameters
parameter
signal
estimated
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.)
Expired - Fee Related
Application number
CN201810379450.4A
Other languages
English (en)
Other versions
CN108614804A (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.)
Information Engineering University of PLA Strategic Support Force
Original Assignee
Information Engineering University of PLA Strategic Support Force
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 Information Engineering University of PLA Strategic Support Force filed Critical Information Engineering University of PLA Strategic Support Force
Priority to CN201810379450.4A priority Critical patent/CN108614804B/zh
Publication of CN108614804A publication Critical patent/CN108614804A/zh
Application granted granted Critical
Publication of CN108614804B publication Critical patent/CN108614804B/zh
Expired - Fee Related 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Algebra (AREA)
  • Evolutionary Biology (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Complex Calculations (AREA)
  • Filters That Use Time-Delay Elements (AREA)

Abstract

本发明提供一种基于信噪比检验的正则化卡尔曼滤波方法。该方法包括：步骤1、获取tk+1时刻待估状态参数Xk+1的初始状态估计
Figure DDA0001640687990000011
和典则参数θk+1的最小二乘估计
Figure DDA0001640687990000012
步骤2、根据待估状态参数Xk+1的第i个参数
Figure DDA0001640687990000013
的状态估计
Figure DDA0001640687990000014
确定
Figure DDA0001640687990000015
的信噪比统计量，i＝1,2…，t‑1，t，t为待估状态参数的个数；步骤3、根据所述
Figure DDA0001640687990000016
的信噪比统计量，将待估状态参数Xk+1分为涉扰状态参数和非涉扰状态参数；步骤4、根据最小二乘估计
Figure DDA0001640687990000017
确定涉扰状态参数的岭参数
Figure DDA0001640687990000018
和非涉扰状态参数的岭参数
Figure DDA0001640687990000019
步骤5、根据所述岭参数
Figure DDA00016406879900000110
和
Figure DDA00016406879900000111
确定待估状态参数Xk+1的修正矩阵以对所述初始状态估计
Figure DDA00016406879900000112
进行修正。本发明在降低状态参数估计方差的同时有效减少了偏差的引入，使得状态参数估计结果在均方误差意义下更优。

Description

基于信噪比检验的正则化卡尔曼滤波方法
技术领域
本发明涉及动态数据处理技术领域，尤其涉及一种基于信噪比检验的正则化卡尔曼滤波方法。
背景技术
卡尔曼Kalman滤波是目前动态数据处理最常用的方法之一，在大地测量、卫星导航和卫星定轨等方面得到了深入的研究和广泛的应用。离散动态***中观测矩阵的病态性会对Kalman滤波状态估计产生很大的的影响，为了克服观测矩阵病态性的影响，提高参数估计的精度，许多学者给出了改进算法。
目前，已有学者从有偏估计的角度提出了解决离散动态***中的病态性问题的一些方法：(1)将有偏估计与Kalman滤波相结合，提出了Biased Kalman Filter。(2)将岭回归与Kalman滤波相结合，通过对增益阵进行修正，以克服观测矩阵病态性对滤波值的不良影响。(3)将Stein压缩估计和岭回归与Kalman滤波相结合，提出了相应的压缩型Kalman滤波和岭型Kalman滤波以及它们的算法，给出了压缩系数和岭参数的选取方法。
实际上，观测矩阵的病态性对每个状态参数估计的危害大小是不同的，这种危害的大小既与状态参数本身数量级大小有关，也与参数所对应观测矩阵数据列参与复共线性的程度有关，但是以上改进方法均没有考虑各个参数受到病态性影响大小的差异，而是对所有参数进行完全一致的修正。
发明内容
本发明提供一种基于信噪比度量的双参数岭型卡尔曼滤波方法，利用信噪比统计量对参数状态估计受病态性影响大小进行度量，根据度量结果对卡尔曼滤波递推过程采取针对性措施，改进了岭型卡尔曼滤波算法，进一步降低了观测矩阵的病态性对状态估计的影响。
本发明提供了一种基于信噪比检验的正则化卡尔曼滤波方法，该方法包括：
步骤1、获取tk+1时刻待估状态参数Xk+1的初始状态估计
Figure GDA0003745529980000021
和典则参数θk+1的最小二乘估计
Figure GDA0003745529980000022
步骤2、根据待估状态参数Xk+1的第i个参数
Figure GDA0003745529980000023
的状态估计
Figure GDA0003745529980000024
确定
Figure GDA0003745529980000025
的信噪比统计量，i＝1,2…,t-1,t，t为待估状态参数的个数；
步骤3、根据所述
Figure GDA0003745529980000026
的信噪比统计量，将待估状态参数Xk+1分为涉扰状态参数和非涉扰状态参数；
步骤4、根据最小二乘估计
Figure GDA0003745529980000027
确定涉扰状态参数的岭参数
Figure GDA0003745529980000028
和非涉扰状态参数的岭参数
Figure GDA0003745529980000029
步骤5、根据所述岭参数
Figure GDA00037455299800000210
和
Figure GDA00037455299800000211
确定待估状态参数Xk+1的修正矩阵以对所述初始状态估计
Figure GDA00037455299800000212
进行修正。
进一步地，该方法还包括：
获取tk+1时刻的法矩阵Nk+1；
若判断获知法矩阵Nk+1的矩阵条件数大于预设条件数阈值，则对tk+1时刻每个状态参数进行信噪比统计。
进一步地，所述步骤2具体为：
根据下式
Figure GDA00037455299800000213
确定
Figure GDA00037455299800000214
的信噪比统计量
Figure GDA00037455299800000215
其中，
Figure GDA00037455299800000216
为
Figure GDA00037455299800000217
的方差。
进一步地，所述步骤3具体为：
若
Figure GDA00037455299800000218
则
Figure GDA00037455299800000219
为涉扰状态参数；
若
Figure GDA00037455299800000220
则
Figure GDA00037455299800000221
为非涉扰状态参数；
其中
Figure GDA00037455299800000222
ω为显著性水平，
Figure GDA00037455299800000223
是中心χ2分布
Figure GDA00037455299800000227
的上侧ω分位点。
进一步地，所述步骤4具体为：
所述岭参数
Figure GDA00037455299800000224
和
Figure GDA00037455299800000225
分别根据下式
Figure GDA00037455299800000226
Figure GDA0003745529980000031
确定；其中，
Figure GDA0003745529980000032
表示
Figure GDA0003745529980000033
中的最大值，
Figure GDA0003745529980000034
s为涉扰状态参数的个数。
进一步地，所述待估状态参数Xk+1的修正矩阵为：
Figure GDA0003745529980000035
本发明的有益效果：
本发明提供的基于信噪比检验的正则化卡尔曼滤波方法，根据每个状态参数的信噪比统计量的大小，将待估状态参数分为涉扰参数和非涉参数两部分；然后根据最小二乘估计
Figure GDA0003745529980000036
确定两个岭参数
Figure GDA0003745529980000037
和
Figure GDA0003745529980000038
从而对这两部分参数的状态估计进行不同强度岭型修正。对于涉扰参数，使其修正岭参数
Figure GDA0003745529980000039
相对较大，对于非涉参数，使其修正岭参数
Figure GDA00037455299800000310
相对较小。这种精细化的处理在降低状态参数估计方差的同时有效减少了岭型Kalman滤波中偏差的引入，使得状态参数估计结果在均方误差意义下更优。
附图说明
图1为本发明实施例提供的基于信噪比检验的正则化卡尔曼滤波方法的流程示意图；
图2为本发明实施例提供的各滤波方法中法矩阵的条件数的变化示意图；
图3为本发明实施例提供的各滤波方法中均方误差的变化示意图；
图4为本发明实施例提供的各滤波方法中状态估计误差的变化示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在对本发明提供的基于信噪比检验的正则化卡尔曼滤波方法进行描述之前，下面对离散动态***、卡尔曼滤波基本方程、岭型卡尔曼滤波状态估计以及离散动态***中的观测矩阵的病态性对卡尔曼滤波状态估计的影响作相关介绍。
本发明实施例采用如下状态空间模型描述离散动态***：
Xk+1＝Φk+1,kXk+Wk (1)
Yk＝HkXk+Vk (2)
式中，Xk∈Rt为离散动态***在时刻tk的状态；Yk∈Rn为对应的观测信号；Φk+1,k为t×t非奇异矩阵，是tk时刻至tk+1时刻的一步状态转移矩阵；Hk为n×t观测矩阵；Wk∈Rt为服从正态分布的输入噪声；Vk∈Rn为服从正态分布的观测噪声。称式(1)为状态方程，式(2)为观测方程。
同时，假设Wk和Vk满足
E(Wk)＝0,Cov(Wk,Wj)＝E(WkWj T)＝Qkδkj
E(Vk)＝0,Cov(Vk,Vj)＝E(VkVj T)＝Rkδkj
Cov(Wk,Vj)＝E(WkVj T)＝0 (3)
式中，Qk为输入噪声的协方差阵，假设Qk为非负定阵，Rk为观测噪声的协方差阵，假设Rk为正定阵；δkj为Kronecher-δ函数，其定义为
Figure GDA0003745529980000041
由式(3)可以看出，Wk和Vk是均值为零、协方差阵为Qk和Rk的不相关白噪声。
下面对卡尔曼滤波算法的基本方程作相关介绍：
状态一步预测：
Figure GDA0003745529980000051
一步预测协方差阵：
Figure GDA0003745529980000052
滤波增益矩阵：
Figure GDA0003745529980000053
状态更新：
Figure GDA0003745529980000054
估计误差方差阵：
Figure GDA0003745529980000055
式(4)～(8)递推卡尔曼滤波的基本方程。给定初始值
Figure GDA0003745529980000056
和P0，根据tk+1时刻的观测量Yk+1，就可递推计算得tk+1时刻的状态估计
Figure GDA0003745529980000057
由卡尔曼滤波基本方程可知，tk+1时刻的状态估计又可表示为：
Figure GDA0003745529980000058
即
Figure GDA0003745529980000059
是式(10)的解：
Figure GDA00037455299800000510
设
Figure GDA00037455299800000511
Nk+1称为卡尔曼滤波的法矩阵。Nk+1为非负定阵，对Nk+1进行特征值分解，
Figure GDA00037455299800000512
于是
Figure GDA00037455299800000513
其中，
Figure GDA00037455299800000514
Figure GDA0003745529980000061
如果观测矩阵Hk+1存在病态性，则Hk+1和
Figure GDA0003745529980000062
的联合作用很可能使得Nk+1存在病态性，实际工作表明，
Figure GDA0003745529980000063
对观测矩阵病态性的控制作用是微弱的，并不能消除观测矩阵病态性对状态估计的不良影响。
由式(12)可以看到Nk+1存在一个或多个小特征值
Figure GDA0003745529980000064
此时若Yk+1和
Figure GDA0003745529980000065
存在很小的观测误差或偏差，上式中的小特征值的倒数会对误差或偏差产生放大作用，从而使得估值偏离真值很远。
为了减小小特征值对观测结果的影响，在上述算法的基础上，提出了岭型卡尔曼滤波方法，tk+1时刻的岭型卡尔曼滤波状态估计为：
Figure GDA0003745529980000066
岭型卡尔曼滤波利用岭参数αk+1对小特征值
Figure GDA0003745529980000067
进行镇压，降低状态估计的方差，从而减弱小特征值对观测误差的放大作用。
但是，岭型卡尔曼滤波存在两个缺陷，一是没有利用病态信息，从而导致岭参数的修正具有盲目性，是对所有参数进行完全一致的修正；二是岭型卡尔曼滤波引入了偏差，岭型卡尔曼滤波引入的偏差在不断的递推过程中有可能被放大，进而影响估计精度。
由上面分析可知，tk+1时刻的状态估计
Figure GDA0003745529980000068
与法矩阵相关(见式(10))，若法矩阵Nk+1存在病态性，则tk+1时刻的状态估计将变得极不稳定，这是观测矩阵的病态性对状态估计施加影响的地方。法矩阵Nk+1存在病态性的原因是它的数据列之间存在复共线性关系，从而导致tk+1时刻的状态估计不好，但是并非所有的状态估计都不好，法矩阵的病态性只对参与复共线性数据列所对应的状态参数估计影响较大，而对未参与复共线性数据列所对应的状态参数估计影响较小。
本发明在岭型卡尔曼滤波方法的基础上，提出一种基于信噪比检验的正则化卡尔曼滤波方法(简称DPRTKF)，对受观测矩阵病态性较小和较大的这两部分状态参数进行不同强度的修正。图1为本发明实施例提供的基于信噪比检验的正则化卡尔曼滤波方法的流程示意图。如图1所示，该方法包括以下步骤：
S101、获取tk+1时刻待估状态参数Xk+1的初始状态估计
Figure GDA0003745529980000071
和典则参数θk+1的最小二乘估计
Figure GDA0003745529980000072
具体地，本发明提供的方法需预先获取tk+1时刻待估状态参数Xk+1的初始状态估计
Figure GDA0003745529980000073
和典则参数θk+1的最小二乘估计
Figure GDA0003745529980000074
由于本发明旨在对现有卡尔曼滤波方法得到的状态估计
Figure GDA0003745529980000075
进行修正，因此需要预先获取利用现有卡尔曼滤波方法得到的初始状态估计
Figure GDA0003745529980000076
以及典则参数θk+1的最小二乘估计
Figure GDA0003745529980000077
例如，给定状态参数的初始值
Figure GDA0003745529980000078
及其均方误差
Figure GDA0003745529980000079
利用卡尔曼滤波基本方程：
状态一步预测方程：
Figure GDA00037455299800000710
一步预测协方差阵方程：
Figure GDA00037455299800000711
滤波增益矩阵方程：
Figure GDA00037455299800000712
可得：
Figure GDA00037455299800000713
Figure GDA00037455299800000714
而典则参数
Figure GDA00037455299800000715
其最小二乘估计为
Figure GDA00037455299800000716
S102、根据待估状态参数Xk+1的第i个参数
Figure GDA00037455299800000717
的状态估计
Figure GDA00037455299800000718
确定
Figure GDA00037455299800000719
的信噪比统计量，i＝1,2…,t-1,t，t为待估状态参数的个数；
S103、根据所述
Figure GDA00037455299800000720
的信噪比统计量，将待估状态参数Xk+1分为涉扰状态参数和非涉扰状态参数；
S104、根据最小二乘估计
Figure GDA00037455299800000721
确定涉扰状态参数的岭参数
Figure GDA00037455299800000722
和非涉扰状态参数的岭参数
Figure GDA0003745529980000081
S105、根据所述岭参数
Figure GDA0003745529980000082
和
Figure GDA0003745529980000083
确定待估状态参数Xk+1的修正矩阵以对所述初始状态估计
Figure GDA0003745529980000084
进行修正。
本发明提供的基于信噪比检验的正则化卡尔曼滤波方法，根据每个状态参数的信噪比统计量的大小，将待估状态参数分为涉扰参数和非涉参数两部分；然后根据最小二乘估计
Figure GDA0003745529980000085
确定两个岭参数
Figure GDA0003745529980000086
和
Figure GDA0003745529980000087
从而对这两部分参数的状态估计进行不同强度岭型修正。对于涉扰参数，使其修正岭参数
Figure GDA0003745529980000088
相对较大，对于非涉参数，使其修正岭参数
Figure GDA0003745529980000089
相对较小。这种精细化的处理在降低状态参数估计方差的同时有效减少了岭型Kalman滤波中偏差的引入，使得状态参数估计结果在均方误差意义下更优。
在上述实施例的基础上，该方法还包括以下步骤：
获取tk+1时刻的法矩阵Nk+1；
若判断获知法矩阵Nk+1的矩阵条件数大于预设条件数阈值，则对tk+1时刻每个状态参数进行信噪比统计。
具体地，若法矩阵Nk+1的矩阵条件数大于预设条件数阈值，则认为观测矩阵为病态矩阵，而观测矩阵的病态性会对待估状态参数的状态估计结果产生不良影响，因此需要对对tk+1时刻每个状态参数进行信噪比统计。如此，可以为后续基于每个状态参数的信噪比统计量对待估状态参数进行分组，将待估状态参数分为受影响较小的非涉扰状态参数和受影响较大的涉扰状态参数。
在上述各实施例的基础上，该方法中的所述步骤S102具体为：根据下式
Figure GDA00037455299800000810
确定
Figure GDA00037455299800000811
的信噪比统计量
Figure GDA00037455299800000812
其中，
Figure GDA00037455299800000813
为
Figure GDA00037455299800000814
的方差。
在上述各实施例的基础上，该方法中的所述步骤S103具体为：
若
Figure GDA00037455299800000815
则
Figure GDA00037455299800000816
为非涉扰状态参数；
若
Figure GDA0003745529980000091
则
Figure GDA0003745529980000092
为涉扰状态参数；
其中
Figure GDA0003745529980000093
ω为显著性水平，
Figure GDA0003745529980000094
是中心χ2分布
Figure GDA0003745529980000095
的上侧ω分位点。
具体地，当
Figure GDA0003745529980000096
时，认为复共线性对相应状态估计
Figure GDA0003745529980000097
的危害比较严重，其估计效果不好；当
Figure GDA0003745529980000098
时，认为复共线性对相应状态估计
Figure GDA0003745529980000099
的危害比较小，其估计效果较好。通过计算信噪比统计量，可将tk+1时刻的状态参数分为两部分
Figure GDA00037455299800000910
其中信噪比统计量较小的s个参数为
Figure GDA00037455299800000911
受复共线性的危害较大，称为涉扰状态参数；信噪比统计量较大的t-s个参数为
Figure GDA00037455299800000912
它们受复共线性的危害较小，称为非涉扰状态参数。对于
Figure GDA00037455299800000913
由于其估计较差，对这部分参数作较大的修正；对于
Figure GDA00037455299800000914
对其作较小的修正。实际应用中，阈值
Figure GDA00037455299800000915
可根据实际情况灵活确定。
在上述各实施例的基础上，该方法中的步骤S104具体为：
所述岭参数
Figure GDA00037455299800000916
和
Figure GDA00037455299800000917
分别根据下式
Figure GDA00037455299800000918
Figure GDA00037455299800000919
确定；其中，
Figure GDA00037455299800000920
表示
Figure GDA00037455299800000921
中的最大值，
Figure GDA00037455299800000922
s为涉扰状态参数的个数。
具体地，将所有信噪比统计量
Figure GDA00037455299800000923
按从小到大顺序重新排列并编号为
Figure GDA00037455299800000924
令
Figure GDA00037455299800000925
Figure GDA00037455299800000926
为信号噪声比，则其倒数
Figure GDA00037455299800000927
为噪声信号比，噪声信号比越大，则对应状态参数的估计受复共线性危害程度就越大。用
Figure GDA0003745529980000101
噪声信号比均值
Figure GDA0003745529980000102
和
Figure GDA0003745529980000103
噪声信号比均值
Figure GDA0003745529980000104
的比值来定量地反应
Figure GDA0003745529980000105
和
Figure GDA0003745529980000106
参数估计受复共线性危害程度之间的比例关系，从而确定ck+1。
在上述各实施例的基础上，该方法中的所述待估状态参数Xk+1的修正矩阵为：
Figure GDA0003745529980000107
具体地，修正矩阵Zk+1为t行t列矩阵，矩阵的前s行是主对角线均为
Figure GDA0003745529980000108
的对角矩阵，用来修正涉扰状态参数
Figure GDA0003745529980000109
矩阵的前后t-s行是主对角线均为
Figure GDA00037455299800001010
的对角矩阵，用来修正非涉扰状态参数
Figure GDA00037455299800001011
下面通过又一具体实施例对本发明进行进一步解释说明。
步1.初始化：给定状态参数的初始值
Figure GDA00037455299800001012
及其均方误差
Figure GDA00037455299800001013
步2.时间更新：
Figure GDA00037455299800001014
Figure GDA00037455299800001015
步3.状态更新：
Figure GDA00037455299800001016
Figure GDA00037455299800001017
Figure GDA00037455299800001018
步4.判断法矩阵
Figure GDA0003745529980000111
的条件数是否大于500，若大于500，则进行步5，否则返回步2。
步5.进行信噪比度量，确定涉扰状态参数和非涉扰状态参数。
步6.确定两个岭参数
Figure GDA0003745529980000112
及参数ck+1。
步7.利用双参数岭型卡尔曼滤波对状态估计进行修正，并计算均方误差矩阵
Figure GDA0003745529980000113
Figure GDA0003745529980000114
其中，
Figure GDA0003745529980000115
步8.令
Figure GDA0003745529980000116
并返回步2，利用卡尔曼滤波对下一时
刻状态参数进行估计。
由上述内容可知：
Figure GDA0003745529980000117
所以
Figure GDA0003745529980000118
是
Figure GDA0003745529980000119
的一个线性变换；
又因为
Figure GDA00037455299800001110
所以
Figure GDA00037455299800001111
是状态参数Xk的一个有偏估计；
同时
Figure GDA00037455299800001112
所以
Figure GDA00037455299800001113
是状态参数Xk的一种压缩型有偏估计。
此外，
Figure GDA0003745529980000121
Figure GDA0003745529980000122
Figure GDA0003745529980000123
由于
Figure GDA0003745529980000124
所以
Figure GDA0003745529980000125
又因为，当选择合适的初始值
Figure GDA0003745529980000126
及其均方误差
Figure GDA0003745529980000127
时，Kalman滤波
Figure GDA0003745529980000128
可以是无偏估计，而普通岭型Kalman滤波
Figure GDA0003745529980000129
和本发明提供的正则化滤波
Figure GDA00037455299800001210
均为Kalman滤波
Figure GDA00037455299800001211
的压缩估计，且
Figure GDA00037455299800001212
对Kalman滤波
Figure GDA00037455299800001213
的压缩强度大于
Figure GDA00037455299800001214
对Kalman滤波
Figure GDA00037455299800001215
的压缩强度，所以从偏差角度看，
Figure GDA00037455299800001216
的偏差大于
Figure GDA00037455299800001217
而小于
Figure GDA00037455299800001218
由于
Figure GDA00037455299800001219
所以从条件数意义下
Figure GDA00037455299800001220
算法条件数小于Kalman滤波估计
Figure GDA00037455299800001221
大于岭型Kalman滤波
Figure GDA00037455299800001222
算法条件数。
综上分析，本发明提供的基于信噪比检验的正则化卡尔曼滤波方法既有普通Kalman滤波所不具备的抗病态性，又克服了一般岭型Kalman滤波盲目而过度的抗病态性，从而减少了过度抗病态性而引起偏差过大的损害，达到了更加精准恰当的效果。
下面通过再一具体实施例对本发明进一步解释说明。
为了验证本文病态性分析的正确性以及新算法的有效性，采用如下算例来进行模拟试验。设一步状态转移矩阵Φk+1,k为
Figure GDA0003745529980000131
设计观测矩阵Hk的5个列向量为
ak1＝[15.57 44.02 20.42 18.74 49.20 44.92 55.48 59.28 94.39 128.0296.00 131.42 127.21 252.90 409.20 463.70 510.22]T
ak2＝[2643 2048 3940 6505 5723 11520 5779 5969 8461 20106 13313 1077145543 36194 34703 39204 86533]T
ak4＝[18.0 9.5 12.8 36.7 35.7 24.0 43.3 46.7 76.7 180.5 60.9 103.7126.8 157.7 169.4 331.4 371.6]T
ak5＝[4.45 6.92 4.28 3.90 5.50 4.60 5.62 5.15 6.18 6.15 5.88 4.68 4.885.57 10.78 7.05 6.35]T
ak3＝2ak1+0.5ak4+ek，ek～N17(0,0.052I)
Hk＝[ak1 ak2 ak3 ak4 ak5]
Qk＝I5，Rk＝0.52×I17分别为状态噪声序列和观测噪声序列的协方差阵，初始值
Figure GDA0003745529980000132
为
Figure GDA0003745529980000133
其中X0＝[200,15,35,16,-2.8,6]Τ为真值，
Figure GDA0003745529980000134
为初始误差矩阵。
本发明实施例中设定观测矩阵数据列中存在近似线性关系，导致法矩阵具有严重的病态性。运用三种滤波算法进行计算，并比较它们的条件数、均方误差和偏差，效果对比如图2至图4所示。
(1)由图2可以看出，普通Kalman滤波算法的法矩阵条件数一直在104的数量级，超过了500的阈值，说明观测矩阵存在复共线性确实会导致法矩阵病态性，而RTKF算法(即图中的普通岭型卡尔曼滤波)和DPRTKF算法(即图中的双参数岭型卡尔曼滤波)的法矩阵的条件数得到了有效地控制。
(2)由图3可以看出RTKF算法和DPRTKF算法有效削弱了Kalman滤波中的病态性，在均方误差意义下，DPRTKF算法优于普通Kalman滤波算法和RTKF算法。
(3)由图4可以看出，相对于Kalman滤波算法，RTKF算法适当牺牲了状态估计的无偏性换取了方差的极大减小，使得估计结果较为稳定，而DPRTKF算法在降低状态估计方差的同时减少了偏差的引入，使得解与真值之间的欧式距离小于Kalman滤波算法及RTKF算法。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (4)

1.基于信噪比检验的正则化卡尔曼滤波方法，应用于离散动态***中病态性观测矩阵的状态参数估计过程，采用状态空间模型描述离散动态***：
Xk+1＝Φk+1,kXk+Wk (1)
Yk＝HkXk+Vk (2)
式中，Xk∈Rt为离散动态***在时刻tk的状态；Yk∈Rn为对应的观测信号；Φk+1,k为t×t非奇异矩阵，是tk时刻至tk+1时刻的一步状态转移矩阵；Hk为n×t观测矩阵；Wk∈Rt为服从正态分布的输入噪声；Vk∈Rn为服从正态分布的观测噪声；其特征在于，所述离散动态***指大地测量***、卫星导航***和卫星定轨***中任意一种，所述方法包括：
步骤1、获取tk+1时刻待估状态参数Xk+1的初始状态估计
Figure FDA0003740170520000011
和典则参数θk+1的最小二乘估计
Figure FDA0003740170520000012
步骤2、根据待估状态参数Xk+1的第i个参数
Figure FDA0003740170520000013
的状态估计
Figure FDA0003740170520000014
确定
Figure FDA0003740170520000015
的信噪比统计量，i＝1,2…，t-1，t，t为待估状态参数的个数；所述步骤2具体为：
根据下式
Figure FDA0003740170520000016
确定
Figure FDA0003740170520000017
的信噪比统计量
Figure FDA0003740170520000018
其中，
Figure FDA0003740170520000019
为
Figure FDA00037401705200000110
的方差；
步骤3、根据所述
Figure FDA00037401705200000111
的信噪比统计量，将待估状态参数Xk+1分为涉扰状态参数和非涉扰状态参数；具体包括：
若
Figure FDA00037401705200000112
则
Figure FDA00037401705200000113
为涉扰状态参数；
若
Figure FDA00037401705200000114
则
Figure FDA00037401705200000115
为非涉扰状态参数；
其中
Figure FDA00037401705200000116
ω为显著性水平，
Figure FDA00037401705200000117
是中心χ2分布
Figure FDA00037401705200000118
的上侧ω分位点；
步骤4、根据最小二乘估计
Figure FDA0003740170520000021
确定涉扰状态参数的岭参数
Figure FDA0003740170520000022
和非涉扰状态参数的岭参数
Figure FDA0003740170520000023
步骤5、根据所述岭参数
Figure FDA0003740170520000024
和
Figure FDA0003740170520000025
确定待估状态参数Xk+1的修正矩阵以对所述初始状态估计
Figure FDA0003740170520000026
进行修正。
2.根据权利要求1所述的方法，其特征在于，在所述步骤1之前还包括：
获取tk+1时刻的法矩阵Nk+1；
若判断获知法矩阵Nk+1的矩阵条件数大于预设条件数阈值，则对tk+1时刻每个状态参数进行信噪比统计。
3.根据权利要求1所述的方法，其特征在于，所述步骤4具体为：
所述岭参数
Figure FDA0003740170520000027
和
Figure FDA0003740170520000028
分别根据下式
Figure FDA0003740170520000029
Figure FDA00037401705200000210
确定；其中，
Figure FDA00037401705200000211
表示
Figure FDA00037401705200000212
中的最大值，
Figure FDA00037401705200000213
s为涉扰状态参数的个数，将所有信噪比统计量
Figure FDA00037401705200000214
按从小到大顺序重新排列并编号为
Figure FDA00037401705200000215
4.根据权利要求3所述的方法，其特征在于，所述待估状态参数Xk+1的修正矩阵为：
Figure FDA0003740170520000031
CN201810379450.4A 2018-04-25 2018-04-25 基于信噪比检验的正则化卡尔曼滤波方法 Expired - Fee Related CN108614804B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810379450.4A CN108614804B (zh) 2018-04-25 2018-04-25 基于信噪比检验的正则化卡尔曼滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810379450.4A CN108614804B (zh) 2018-04-25 2018-04-25 基于信噪比检验的正则化卡尔曼滤波方法

Publications (2)

Publication Number Publication Date
CN108614804A CN108614804A (zh) 2018-10-02
CN108614804B true CN108614804B (zh) 2022-09-02

Family

ID=63660835

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810379450.4A Expired - Fee Related CN108614804B (zh) 2018-04-25 2018-04-25 基于信噪比检验的正则化卡尔曼滤波方法

Country Status (1)

Country Link
CN (1) CN108614804B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE102020115919A1 (de) * 2020-06-17 2021-12-23 Fette Compacting Gmbh Verfahren zum Betreiben einer Mischeinrichtung einer Anlage

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102064798A (zh) * 2010-12-17 2011-05-18 北京大学 一种负反馈自适应在线实时滤波方法及***
WO2012177335A1 (en) * 2011-06-21 2012-12-27 Exxonmobil Upstream Research Company Improved dispersion estimation by nonlinear optimization of beam-formed fields
CN103684349A (zh) * 2013-10-28 2014-03-26 北京理工大学 一种基于递推协方差阵估计的卡尔曼滤波方法
CN103902834A (zh) * 2014-04-14 2014-07-02 重庆大学 一种基于岭估计和l曲线法的结构损伤识别方法
CN104376383A (zh) * 2014-11-27 2015-02-25 东北大学 一种基于地理信息***的电网电压监视与预测***及方法
CN107425548A (zh) * 2017-09-11 2017-12-01 河海大学 一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法
CN107729845A (zh) * 2017-10-20 2018-02-23 开沃新能源汽车集团有限公司 一种基于子空间特征值分解的实测频响函数降噪方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
TW200803262A (en) * 2006-06-23 2008-01-01 Pixart Imaging Inc Bumping algorithm used to estimate signal frequency and displacement-estimating method and device

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102064798A (zh) * 2010-12-17 2011-05-18 北京大学 一种负反馈自适应在线实时滤波方法及***
WO2012177335A1 (en) * 2011-06-21 2012-12-27 Exxonmobil Upstream Research Company Improved dispersion estimation by nonlinear optimization of beam-formed fields
CN103684349A (zh) * 2013-10-28 2014-03-26 北京理工大学 一种基于递推协方差阵估计的卡尔曼滤波方法
CN103902834A (zh) * 2014-04-14 2014-07-02 重庆大学 一种基于岭估计和l曲线法的结构损伤识别方法
CN104376383A (zh) * 2014-11-27 2015-02-25 东北大学 一种基于地理信息***的电网电压监视与预测***及方法
CN107425548A (zh) * 2017-09-11 2017-12-01 河海大学 一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法
CN107729845A (zh) * 2017-10-20 2018-02-23 开沃新能源汽车集团有限公司 一种基于子空间特征值分解的实测频响函数降噪方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A novel GLM-based method for the Automatic IDentification of functional Events (AIDE) in fNIRS data recorded in naturalistic environments;Pinti Paola 等;《Neuroimage》;20170731;第155卷;291-304 *
Nonparametric hemodynamic deconvolution of fMRI using homomorphic filtering;Sreenivasan K. R. 等;《IEEE Transactions on Medical Imaging》;20141212;第34卷(第5期);1155-1163 *
一种改进的高斯近似滤波方法;黄玉龙 等;《自动化学报》;20151231;第42卷(第3期);385-401 *
基于信噪比检验的双参数岭型Kalman滤波及其在BDS星地联合定轨中的应用;李豪 等;《大地测量与地球动力学》;20181115;第38卷(第11期);1143-1148 *

Also Published As

Publication number Publication date
CN108614804A (zh) 2018-10-02

Similar Documents

Publication Publication Date Title
US20210159964A1 (en) Direction-of-arrival estimation and mutual coupling calibration method and system with arbitrary sensor geometry and unknown mutual coupling
CN111985093B (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
Brüggemann et al. Inference in VARs with conditional heteroskedasticity of unknown form
US10754922B2 (en) Method and apparatus for sensor fusion
CN110061716B (zh) 一种基于最小二乘和多重渐消因子的改进kalman滤波方法
EP2784537A1 (en) Inversion method and apparatus based on polarimetric interferometric synthetic aperture radar
Felus Application of total least squares for spatial point process analysis
CN107870001A (zh) 一种基于椭球拟合的磁力计校正方法
CN109239653B (zh) 一种基于子空间分解的多辐射源被动直接时差定位方法
CN103293517B (zh) 基于脊参数估计的对角加载稳健自适应雷达波束形成方法
CN110826021B (zh) 一种非线性工业过程鲁棒辨识和输出估计方法
CN110196410B (zh) 一种阵列天线主瓣干扰抑制方法及***
Metref et al. A non-Gaussian analysis scheme using rank histograms for ensemble data assimilation
CN117236084B (zh) 一种木工机械加工生产智能管理方法及***
CN110032709B (zh) 一种用于地理坐标转换中异常点的定位与估值方法
CN108614804B (zh) 基于信噪比检验的正则化卡尔曼滤波方法
CN113238227B (zh) 一种结合深度学习的改进最小二乘相位解缠方法及***
CN111398523B (zh) 基于分布的传感器数据校准方法及***
Lederer et al. Estimating the Lasso's effective noise
CN111142062B (zh) 一种利用Toeplitz特性的无栅格目标波达方向估计方法
CN113569393B (zh) 一种稳健的绝对定向算法
CN114048431A (zh) 一种基于协方差矩阵重构和admm的波束形成方法
CN103684349B (zh) 一种基于递推协方差阵估计的卡尔曼滤波方法
CN114266223B (zh) 机台的故障确定方法、装置、设备和计算机可读存储介质
CN111175692A (zh) 基于分层合成Lasso先验模型的离格稀疏贝叶斯DOA估计方法

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20220902

CF01 Termination of patent right due to non-payment of annual fee