CN110287537B - 用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法 - Google Patents

用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法 Download PDF

Info

Publication number
CN110287537B
CN110287537B CN201910447795.3A CN201910447795A CN110287537B CN 110287537 B CN110287537 B CN 110287537B CN 201910447795 A CN201910447795 A CN 201910447795A CN 110287537 B CN110287537 B CN 110287537B
Authority
CN
China
Prior art keywords
frequency
value
matrix
frequency deviation
data
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.)
Active
Application number
CN201910447795.3A
Other languages
English (en)
Other versions
CN110287537A (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.)
Northwest University
Original Assignee
Northwest 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 Northwest University filed Critical Northwest University
Priority to CN201910447795.3A priority Critical patent/CN110287537B/zh
Publication of CN110287537A publication Critical patent/CN110287537A/zh
Application granted granted Critical
Publication of CN110287537B publication Critical patent/CN110287537B/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/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00—Adaptive networks
    • H03H21/0012—Digital adaptive filters
    • H03H21/0025—Particular filtering methods
    • H03H21/0029—Particular filtering methods based on statistics
    • H03H21/003—KALMAN filters
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00—Adaptive networks
    • H03H21/0012—Digital adaptive filters
    • H03H21/0043—Adaptive algorithms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Computing Systems (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Complex Calculations (AREA)
  • Filters That Use Time-Delay Elements (AREA)

Abstract

本发明属于频率标准故障检测技术领域，公开了一种用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法，所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法包括：野值自适应卡尔曼滤波利用观测数据进行递推滤波的同时，根据事实的得到的量测新息的实际方差与理论方差的比值；根据先验频率偏差矩阵得到后验频率偏差的矩阵，据先验频率偏差的先验均方误差得到频率偏差的增益矩阵，根据增益矩阵得到后验频率偏差均方误差的矩阵。带有频率跳变或者野值的样本数据常常会使卡尔曼滤波器对系统的状态预报内进行错误的修正，使得滤波的结果发生偏移，甚至发散，在新息进行加权后可以使滤波后的数据更加接近真实状态。

Description

用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法
技术领域
本发明属于频率标准故障检测技术领域，尤其涉及一种用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法。
背景技术
针对频率标准跳变的检测方法现如今使用最多的就是最小二乘法、块平均、连续平均、最大似然估计、卡尔曼滤波算法及动态阿伦方差算法。在使用最小二乘法时，需要自动忽略变量误差后进行曲线上相应点的拟合，并且使用这种方法进行频率跳变检测时会出现检测精度不够高的问题，并且检测较多的数据时会使用较长的检测时间；块平均是在相邻窗口中的平均值与选定阈值之间进行对比的方法，这种检测方法只适用于对少量的数据进行处理，处理较多的数据时需要使用较长的时间；连续平均是一种使用搜索跳转序列的检测方法，这种方法在进行频率跳变检测时会出现假跳的现象，误检率相对较高；最大似然估计可以较为快速的检测到频率跳变，但是确是以牺牲频率跳变检测概率才能达到；动态阿伦方差方法是近几年的一种新方法，这种方法在检测频率跳变时需要使用过去的数据进行检测，所以检测速率较为缓慢，因此通常都被用来检测频率标准的稳定性，可以当做检测频率标准是否发生异常的一种方法。
卡尔曼滤波器，就目前来讲，是一个研究大热门，因为卡尔曼滤波器是一种无偏线性最小方差估计算法，属于预测性滤波算法，具有递归的优点，计算成本低，对于频率跳变具有相对较高的检测概率；同时在进行滤波时可以基于内部的分析来更新数据，无需附加方程来检测异常数据，但是卡尔曼滤波器如果在无法确定被研究对象的精确数学模型及精确噪声统计特性时，实际测量数目会不断的增加，进而导致滤波的估计误差的均值和估计误差协方差可能越来越大，最终使得滤波逐渐失去准确估计的作用，无法直接真实的反映物理过程，使得模型与获得的量测量值不匹配，滤波精度就会被削弱，严重时甚至会造成滤波发散，由于模型建立过于粗糙或失真所引起的发散称为滤波发散。并且随着卡尔曼滤波递推算法的过程，滤波的步数的增加，会导致舍入误差逐渐累积，一旦计算机字长不足，就有一定的风险使得估计误差方差真失去其非负定性甚至失去对称性，导致滤波增益矩阵逐渐失去合适的加权作用而发生发散。同时，测量仪器的不精确可能会导致单个野值的出现，野值的存在会使得滤波器误判为频标故障，导致实验结果出现错误。一般在实践中，噪声的统计特性无法被精确估计，同时伴随着科技的发展，对滤波检测概率有了更高的要求，卡尔曼滤波无法满足所需的精确度以及检测概率，同时在容错能力上，卡尔曼滤波算法确实比较不足。在实际的工程应用中，由于工作环境、测量环境的影响，或者测量仪器的测量问题导致个别数据出现突发性的误差，这些值我们称它为野值。在进行频差数据测量时，需要高精度的测量，测量方法复杂，进行测量时由于不可控的因素导致出现异常的离群值，对于这些离群值只有剔除影响才能得到更精确的数据。但是卡尔曼滤波器对于野值缺少一定的容错能力。所以需要有新的技术来提高滤波检测概率以及容错能力。
综上所述，现有技术存在的问题是：卡尔曼滤波器在无法确定被研究对象的精确数学模型及精确噪声统计特性时，滤波精度就会被大大削弱，严重时甚至会造成滤波发散。解决上述技术问题的难度：现实情况中，卡尔曼滤波器模型无法做到十分精确，同时噪声统计特性大多以估计计算为准，所以使用卡尔曼滤波器时，可能会出现滤波精度不足的问题，并且可能会出现错判的情况，所以如何能够将数据进行自动调整，识别出野值的存在，以及自适应的调整数据精度就是现有的一个问题。
解决上述技术问题的意义：如果能够添加一个抗野值自适应因子根据实际数据进行自适应调整滤波后数据，那么将能够自动识别野值并且进行自适应的调整，同时能够根据测量数据以及先验数据调整提高现有滤波精度。
发明内容
针对现有技术存在的问题，本发明提供了一种用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法。
本发明是这样实现的，一种用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法，所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法包括：
第一步，抗野值自适应卡尔曼滤波器利用观测数据进行递推滤波的同时，根据事实得到的量测新息的实际方差与理论方差的比值，得到平均频率偏差的滤波方程；
第二步，根据先验频率偏差矩阵得到后验频率偏差的矩阵，据先验频率偏差的先验均方误差得到频率偏差的增益矩阵，根据增益矩阵得到后验频率偏差均方误差的矩阵。
进一步，所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法的测量模型为yk＝hT·Xk+vk，yk为第k个采样点的观测向量的平均频率偏差，vk表示协方差矩阵Pk的测量噪声；
进一步，野值自适应卡尔曼滤波利用观测数据进行递推滤波的同时，根据事实的得到的量测新息的实际方差与理论方差的比值；得到平均频率偏差的滤波方程，从而得到频率偏差的新息为
Figure BDA0002074182600000031
其中Zk为k时刻的测量值，h为测量系数矩阵，根据先验频率偏差矩阵得到后验频率偏差的矩阵
Figure BDA0002074182600000032
据先验频率偏差的先验均方误差得到频率偏差的增益矩阵，Kk＝Pk·hT(hPkhT+R)-1，根据增益矩阵得到后验频率偏差均方误差的矩阵
Figure BDA0002074182600000033
进一步，所述测量模型为卡尔曼滤波器时钟模型，考虑白色频率噪声WFN以及随机游走噪声RWFN两种时钟噪声，y(k)＝yW(k)+yRW(k)+d(k)，yW(k)为白色频率噪声，为均值为0，方差为
Figure BDA0002074182600000034
的统计独立高斯随机噪声，yRW(k)为随机游走噪声，yRW(k)由yRW(k-1)及ΔyRW(k-1)组成，其中yRW(0)＝0；
另一部分ΔyRW(k-1)为均值为0，方差为
Figure BDA0002074182600000041
的统计独立的高斯随机变量，d(k)为确定性慢频率漂移，是一个递推函数，d(k)＝d(k-1)+T·d，T为测量时间间隔，为1s,d(0)代表标准频率偏移；d为一个常量代表线性漂移；Xk为频率标准在时间tk的平均频率表示，状态预测向量为
Figure BDA0002074182600000042
状态估计向量为
Figure BDA0002074182600000043
Xk＝A·Xk-1+B·Uk，
Figure BDA0002074182600000044
其中，A代表状态转移系数矩阵，为
Figure BDA0002074182600000045
是表示状态量从Xk-1表达Xk的系数矩阵，B表示模型误差的系数转移矩阵，为
Figure BDA0002074182600000046
第k-1个采样点的模型误差协方差用Qk-1表示。
进一步，所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法为后验频率偏差
Figure BDA0002074182600000047
其中Φk(τk)为抗野值修正因子，
Figure BDA0002074182600000048
Φk为适当的选取的按光滑的函数，称为压缩影响函数；Dk为权矩阵；
当选取Φk＝1时，滤波为自适应卡尔曼滤波；Φk≠1时，则为提到的抗野值自适应滤波；若出现野值而τk增大时，Φk减小，减小野值的影响权重；其中Ck为适当选取的门限常数序列，取值为
Figure BDA0002074182600000049
λk为
Figure BDA00020741826000000410
矩阵的最大特征值。
Figure BDA00020741826000000411
为χ2(m)分布的置信度为(1-α)×100％的上分为点，通常α取0.05或0.025；m为量测数据的维数。
进一步，所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法具体包括：先验频率偏差数据和后验频率偏差数据相互影响，在通过当前时刻的先验概率偏差数据得到当前的时刻后验频率偏差需要加入抗野值自适应因子，当在k时刻出现由于测量仪器故障导致的野值时，抗野值因子根据k时刻之前的数据通过τk和
Figure BDA0002074182600000051
进行自适应的调整，当τk小于
Figure BDA0002074182600000052
时，则此时抗野值因子取值为1，如果τk大于
Figure BDA0002074182600000053
时，则抗野值因子取值为
Figure BDA0002074182600000054
也就是说当此时频差数据单个数据突变，则通过该式进行自适应的修正，最终将无异常的数据保持在合理的范围内。
进一步，抗野值自适应卡尔曼滤波跳变检测需要测得观测频差数据，经过卡尔曼滤波得到先验频率偏差，通过计算得到频率偏差的新息，然后通过加权平均后和检测阈值进行比较；具体步骤如下：
步骤一，选择最优累计数N，将铷原子钟测量的频差数据划分窗口，构建提出的检验统计量的预测残留；
步骤二，通过抗野值自适应卡尔曼滤波器计算出在时间tk时的更新矩阵
Figure BDA0002074182600000055
根据最优累计数N进行N步预测，得出
Figure BDA0002074182600000056
步骤三，根据N步预测通过式
Figure BDA0002074182600000057
得到相应的N步更新矩阵，得出
Figure BDA0002074182600000058
步骤四，在第k个采样周期，通过测量频率偏差和预测频率偏差的差值，得到n步预测残差，得到矩阵ε＝[εk+1 εk+2 … εk+N]T，其中，
Figure BDA0002074182600000059
步骤五，将得到N个采样周期的预测残差进行平均，得到检测统计量
Figure BDA00020741826000000510
步骤六，将每个检测统计量的值和检测阈值T进行比较，大于阈值则判定为铷频率标准发生频率跳变，小于阈值T则铷频率标准正常运行。
本发明的另一目的在于提供一种应用所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法的卡尔曼滤波器。
本发明的另一目的在于提供一种应用所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法的频率跳变检测系统。
本发明的另一目的在于提供一种应用所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法的物联网信息处理系统。
综上所述，本发明的优点及积极效果为：本发明的频率标准的频率跳变检测算法抗野值自适应卡尔曼滤波器在进行滤波过程中可以有效防止滤波发散，在滤波过程中当前测量值在滤波中占有较大比重，需要对相应的量测噪声有相对精确的估计。为了卡尔曼滤波能够更好的抗野值，本发明的方法使用将新息序列进行加权的方法来消除野值的影响，也就是公式中的抗野值因子的部分。带有频率跳变或者野值的样本数据常常会使卡尔曼滤波器对系统的状态预报内进行错误的修正，使得滤波的结果发生偏移，甚至发散，在新息进行加权后可以使滤波后的数据更加接近真实状态。
附图说明
图1是本发明实施例提供的用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法流程图。
图2是本发明实施例提供的用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法实现流程图。
图3是本发明实施例提供的跳变检测流程图。
图4是本发明实施例提供的故障模拟系统在t＝7390s时刻产生野值的频率偏差数据显示示意图。
图5是本发明实施例提供的抗野值自适应卡尔曼滤波器对野值的抑制效果图。
图6是本发明实施例提供的卡尔曼滤波器对野值的抑制效果图。
图7是本发明实施例提供的卡尔曼滤波在频率标准上的状态估计示意图。
图8是本发明实施例提供的抗野值自适应卡尔曼滤波在频率标准上的状态估计示意图。
图9是本发明实施例提供的抗野值自适应卡尔曼滤波与卡尔曼滤波频率跳变的累积数和检测概率之间的关系示意图。
图10是本发明实施例提供的铷频率标准PRS10频率跳变大小和检测概率(PD)的关系图曲线示意图。
图11是本发明实施例提供的两种不同滤波针对不同频率跳变对应的检测概率(PD)示意图。
图12是本发明实施例提供的使用原子钟测得一组数据示意图。
图13是本发明实施例提供的通过卡尔曼滤波器检测结果示意图。
图14是本发明实施例提供的抗野值自适应卡尔曼滤波器的检测结果示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
针对卡尔曼滤波器在无法确定被研究对象的精确数学模型及精确噪声统计特性时，滤波精度就会被大大削弱，严重时甚至会造成滤波发散；本发明提供了一种抗野值自适应卡尔曼滤波方法，不仅能够保持传统卡尔曼滤波器递归的优点，同时具有低的计算成本，并且在传统卡尔曼滤波的基础上提高了频率跳变的检测概率，有效的防止滤波扩散以及对野值的抑制。
下面结合附图对本发明的应用原理作详细的描述。
如图1所示，本发明实施例提供的用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法包括以下步骤：
S101：野值自适应卡尔曼滤波利用观测数据进行递推滤波的同时，根据事实的得到的量测新息的实际方差与理论方差的比值；
S102：根据先验频率偏差矩阵得到后验频率偏差的矩阵，据先验频率偏差的先验均方误差得到频率偏差的增益矩阵，根据增益矩阵得到后验频率偏差均方误差的矩阵。
所述算法模型为卡尔曼滤波器时钟模型，并且只考虑白色频率噪声(WFN)以及随机游走噪声(RWFN)两种时钟噪声，y(k)＝yW(k)+yRW(k)+d(k)，yW(k)为白色频率噪声，为均值为0，方差为
Figure BDA0002074182600000081
的统计独立高斯随机噪声，yRW(k)为随机游走噪声，yRW(k)由yRW(k-1)及ΔyRW(k-1)组成，其中yRW(0)=0，另一部分ΔyRW(k-1)为均值为0，方差为
Figure BDA0002074182600000082
的统计独立的高斯随机变量，d(k)为确定性慢频率漂移，是一个递推函数，d(k)＝d(k-1)+T·d，T为测量时间间隔，为1s,d(0)代表标准频率偏移；d为一个常量代表线性漂移。Xk为频率标准在时间tk的平均频率表示，状态预测向量为
Figure BDA0002074182600000083
状态估计向量为
Figure BDA0002074182600000084
Xk＝A·Xk-1+B·Uk，
Figure BDA0002074182600000085
其中，A代表状态转移系数矩阵，为
Figure BDA0002074182600000086
是表示状态量从Xk-1表达Xk的系数矩阵，B表示模型误差的系数转移矩阵，为
Figure BDA0002074182600000087
第k-1个采样点的模型误差协方差用Qk-1表示。
所述算法的测量模型为yk＝hT·Xk+vk，yk为第k个采样点的观测向量的平均频率偏差，vk表示协方差矩阵Pk的测量噪声。野值自适应卡尔曼滤波利用观测数据进行递推滤波的同时，根据事实的得到的量测新息的实际方差与理论方差的比值。如此得到平均频率偏差的滤波方程，以此得到该频率偏差的信息为
Figure BDA0002074182600000088
其中Zk为k时刻的测量值，h为测量系数矩阵，根据先验频率偏差矩阵得到后验频率偏差的矩阵
Figure BDA0002074182600000089
据先验频率偏差的先验均方误差得到频率偏差的增益矩阵，Kk＝Pk·hT(hPkhT+R)-1，根据增益矩阵得到后验频率偏差均方误差的矩阵
Figure BDA00020741826000000810
所述算法的核心算法为后验频率偏差
Figure BDA0002074182600000091
其中Φk(τk)为抗野值修正因子，
Figure BDA0002074182600000092
Φk为适当的选取的按光滑的函数，称为压缩影响函数；Dk为权矩阵。当选取Φk＝1时，上述滤波即为自适应卡尔曼滤波；而Φk≠1时，则为提到的抗野值自适应滤波。若出现野值而τk增大时，Φk减小，从而减小野值的影响权重。其中Ck为适当选取的门限常数序列，取值为
Figure BDA0002074182600000093
λk为
Figure BDA0002074182600000099
矩阵的最大特征值。
Figure BDA0002074182600000094
为χ2(m)分布的置信度为(1-α)×100％的上分为点，通常α取0.05或0.025；m为量测数据的维数。
下面结合附图对本发明的应用原理作进一步的描述。
以铷原子钟输出频率跳变频差数据为参考。
如图2所示，后一时刻的先验频率偏差通过前一时刻的后验频率偏差递推而来，直接影响下一时刻的后验频率偏差。在进行检测时，先验频率偏差数据和后验频率偏差数据相互影响，所以改变后验频率偏差则会改变先验频率偏差，所以在通过当前时刻的先验概率偏差数据得到当前的时刻后验频率偏差需要加入抗野值自适应因子，当在k时刻时出现了由于测量仪器故障导致的野值时，抗野值因子是一个压缩影响函数，会根据k时刻之前的数据通过τk和
Figure BDA0002074182600000095
进行自适应的调整，当τk小于
Figure BDA0002074182600000096
时，则此时抗野值因子取值为1，如果τk大于
Figure BDA0002074182600000097
时，则抗野值因子取值为
Figure BDA0002074182600000098
也就是说当此时频差数据单个数据突变，则通过该式进行自适应的修正，最终将无异常的数据保持在合理的范围内。
如图3所示，抗野值自适应卡尔曼滤波跳变检测需要测得观测频差数据，经过卡尔曼滤波得到先验频率偏差，通过计算得到频率偏差的新息，然后通过加权平均后得到的后验频率偏差，将该值和检测阈值进行比较。具体步骤如下：
步骤一：选择最优累计数N，将铷原子钟测量的频差数据划分窗口，构建提出的检验统计量的预测残留。
步骤二：通过抗野值自适应卡尔曼滤波器计算出在时间tk时的更新矩阵
Figure BDA0002074182600000101
根据最优累计数N进行N步预测，得出
Figure BDA0002074182600000102
步骤三：根据N步预测通过式
Figure BDA0002074182600000103
得到相应的N步更新矩阵，得出
Figure BDA0002074182600000104
步骤四：在第k个采样周期，通过测量频率偏差和预测频率偏差的差值，得到n步预测残差，得到矩阵ε＝[εk+1 εk+2 … εk+N]T，其中，
Figure BDA0002074182600000105
步骤五：将得到N个采样周期的预测残差进行平均，得到检测统计量
Figure BDA0002074182600000106
步骤六：将每个检测统计量的值和检测阈值T进行比较，大于阈值则判定为铷频率标准发生频率跳变，小于阈值T则铷频率标准正常运行。
下面结合实验对本发明的技术效果作详细的描述。
以PRS10铷原子钟为基准源，搭建一套故障模拟系统进行算法验证，该系统采样时间统为t＝1。图4为故障模拟的野值频差数据显示。在当t＝7390s时产生了一个1.5×10-10的野值，使用抗野值自适应卡尔曼滤波器和卡尔曼滤波器的抑制效果图分别如图5和图6。野值的出现使抗野值因子减小，从而使野值的影响减少，提高滤波精度。
当故障模拟系统出现2.5×10-11跳变时，使用卡尔曼滤波器的滤波效果以及抗野值自适应卡尔曼滤波器的滤波效果如图7和8所示。
检测概率、误检率和累积数是验证滤波效果的几种常用的数据，系统频率跳变设置如图7和图8，通过频率跳变滤波检测的两种检测方法的累积数和检测概率的关系验证新方法的性能优势，如图9所示。抗野值自适应卡尔曼滤波器在达到同样的检测概率时累积数为16远远小于卡尔曼滤波器的25。
设定累积数为11，误检率为10-3，通过检测概率和误检率之间的关系来进行仿真验证。
如图10所示，相同频率跳变，抗野值自适应卡尔曼滤波的检测概率高于卡尔曼滤波，当抗野值自适应卡尔曼滤波检测概率到达0.99时，卡尔曼滤波的检测概率只有0.2左右。
如图11所示，当累积数为11，频率跳变为2.5×10-11时，两种方法的检测概率都会随着误检率的增加而增加。在两种滤波方法针对频率跳变的误检率都为1×10-3时，卡尔曼滤波的检测概率约为0.6，但是抗野值自适应卡尔曼滤波的检测概率在0.95左右，比卡尔曼滤波的检测概率提高了0.35左右。
如图12所示，使用原子钟测得一组数据。原子钟在t＝3000s处发生了2.5×10-11的频率跳变，通过卡尔曼滤波器和抗野值自适应卡尔曼滤波器的检测结果分别如图13和14所示。卡尔曼滤波器的不精确性可能会导致滤波器无法正确检测出跳变故障的发生。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法，其特征在于，所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法包括：
第一步，野值自适应卡尔曼滤波利用观测数据进行递推滤波的同时，根据事实的得到的量测新息的实际方差与理论方差的比值，得到平均频率偏差的滤波方程；
第二步，根据先验频率偏差矩阵得到后验频率偏差的矩阵，据先验频率偏差的先验均方误差得到频率偏差的增益矩阵，根据增益矩阵得到后验频率偏差均方误差的矩阵；
所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法的测量模型为yk＝hT·Xk+vk，yk为第k个采样点的观测向量的平均频率偏差，vk表示协方差矩阵Pk的测量噪声；
野值自适应卡尔曼滤波利用观测数据进行递推滤波的同时，根据事实的得到的量测新息的实际方差与理论方差的比值；得到平均频率偏差的滤波方程，从而得到频率偏差的信息为
Figure FDA0004134977000000011
其中Zk为k时刻的测量值，h为测量系数矩阵，根据先验频率偏差矩阵得到后验频率偏差的矩阵
Figure FDA0004134977000000012
据先验频率偏差的先验均方误差得到频率偏差的增益矩阵，Kk＝Pk·hT(hPkhT+R)-1，根据增益矩阵得到后验频率偏差均方误差的矩阵
Figure FDA0004134977000000013
所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法具体包括：先验频率偏差数据和后验频率偏差数据相互影响，在通过当前时刻的先验概率偏差数据得到当前的时刻后验频率偏差需要加入抗野值自适应因子，当在k出现由于测量仪器故障导致的野值时，抗野值因子根据k时刻之前的数据通过τk和
Figure FDA0004134977000000014
进行自适应的调整，当τk小于等于
Figure FDA0004134977000000015
时，则此时抗野值因子取值为1，如果τk大于
Figure FDA0004134977000000021
时，则抗野值因子取值为
Figure FDA0004134977000000022
也就是说当此时频差数据单个数据突变，则通过
Figure FDA0004134977000000023
进行自适应的修正，最终将无故障异常的数据保持在合理的范围内；
抗野值自适应卡尔曼滤波跳变检测需要测得观测频差数据，经过卡尔曼滤波得到先验频率偏差，通过计算得到频率偏差的新息，然后通过加权平均后得到后验频率偏差，将该值和检测阈值进行比较；具体步骤如下：
步骤一，选择最优累计数N，将铷原子钟测量的频差数据划分窗口，构建提出的检验统计量的预测残留；
步骤二，通过抗野值自适应卡尔曼滤波器计算出在时间tk时的更新矩阵
Figure FDA0004134977000000024
根据最优累计数N进行N步预测，得出
Figure FDA0004134977000000025
步骤三，根据N步预测通过式
Figure FDA0004134977000000026
得到相应的N步更新矩阵，得出
Figure FDA0004134977000000027
步骤四，在第k个采样周期，通过测量频率偏差和预测频率偏差的差值，得到n步预测残差，得到矩阵ε＝[εk+1 εk+2…εk+N]T，其中，
Figure FDA0004134977000000028
步骤五，将得到N个采样周期的预测残差进行平均，得到检测统计量
Figure FDA0004134977000000029
步骤六，将每个检测统计量的值和检测阈值T进行比较，大于阈值则判定为铷频率标准发生频率跳变，小于阈值T则铷频率标准正常运行。
2.如权利要求1所述的用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法，其特征在于，所述测量模型为卡尔曼滤波器时钟模型，考虑白色频率噪声WFN以及随机游走噪声RWFN两种时钟噪声，y(k)＝yW(k)+yRW(k)+d(k)，yW(k)为白色频率噪声，均值为0，方差为
Figure FDA00041349770000000210
的统计独立高斯随机噪声，yRW(k)为随机游走噪声，yRW(k)由yRW(k-1)及ΔyRW(k-1)组成，其中yRW(0)＝0；
另一部分ΔyRW(k-1)为均值为0，方差为
Figure FDA0004134977000000031
的统计独立的高斯随机变量，d(k)为确定性慢频率漂移，是一个递推函数，d(k)＝d(k-1)+T·d，T为测量时间间隔，为1s,d(0)代表标准频率偏移；d为一个常量代表线性漂移；Xk为频率标准在时间tk的平均频率表示，状态预测向量为
Figure FDA00041349770000000311
状态估计向量为
Figure FDA0004134977000000032
Xk＝A·Xk-1+B·Uk，
Figure FDA0004134977000000033
其中，A代表状态转移系数矩阵，为
Figure FDA0004134977000000034
是表示状态量从Xk-1表达Xk的系数矩阵，B表示模型误差的系数转移矩阵，为
Figure FDA0004134977000000035
第k-1个采样点的模型误差协方差用Qk-1表示。
3.如权利要求1所述的用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法，其特征在于，所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法为后验频率偏差
Figure FDA0004134977000000036
其中Φk(τk)为抗野值修正因子，
Figure FDA0004134977000000037
Φk为适当的选取的按光滑的函数，称为压缩影响函数；Dk为权矩阵；
当选取Φk＝1时，滤波为自适应卡尔曼滤波；Φk≠1时，则为提到的抗野值自适应滤波；若出现野值而τk增大时，Φk减小，减小野值的影响权重；其中Ck为适当选取的门限常数序列，取值为
Figure FDA0004134977000000038
λk为
Figure FDA0004134977000000039
矩阵的最大特征值；
Figure FDA00041349770000000310
为χ2(m)分布的置信度为(1-α)×100％的上分为点，α取0.05或0.025；m为量测数据的维数。
4.一种应用权利要求1～3任意一项所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法的频率跳变检测系统。
5.一种应用权利要求1～3任意一项所述用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法的物联网信息处理系统。
CN201910447795.3A 2019-05-27 2019-05-27 用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法 Active CN110287537B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910447795.3A CN110287537B (zh) 2019-05-27 2019-05-27 用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910447795.3A CN110287537B (zh) 2019-05-27 2019-05-27 用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法

Publications (2)

Publication Number Publication Date
CN110287537A CN110287537A (zh) 2019-09-27
CN110287537B true CN110287537B (zh) 2023-05-05

Family

ID=68002571

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910447795.3A Active CN110287537B (zh) 2019-05-27 2019-05-27 用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法

Country Status (1)

Country Link
CN (1) CN110287537B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110991015B (zh) * 2019-11-19 2023-05-02 三峡大学 融合卡尔曼和小波的mems陀螺仪自适应抗野值去噪方法
CN113176725B (zh) * 2021-03-05 2023-01-24 北京大学 基于卡尔曼滤波和/或dfb的激光芯片原子钟及实现方法
CN114710252B (zh) * 2022-03-17 2023-05-16 陕西国防工业职业技术学院 一种用于精密时钟同步的滤波方法及系统
CN117590838B (zh) * 2024-01-19 2024-05-03 深圳市凯威达电子有限公司 一种微处理器老化智能检测方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102323602A (zh) * 2011-05-30 2012-01-18 哈尔滨工程大学 一种基于自适应二阶卡尔曼滤波器的载波跟踪环路及其滤波方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6898582B2 (en) * 1998-12-30 2005-05-24 Algodyne, Ltd. Method and apparatus for extracting low SNR transient signals from noise
CN101793522B (zh) * 2010-03-31 2011-04-13 上海交通大学 基于抗差估计的稳健滤波方法
CN106324521B (zh) * 2016-09-05 2018-09-11 北京理工大学 一种联合估计动力电池系统参数与荷电状态的方法
CN107576330A (zh) * 2017-09-07 2018-01-12 西北大学 一种基于wlan指纹的室内动态感知策略的定位方法
CN107832575B (zh) * 2017-10-10 2021-07-16 中国航空无线电电子研究所 基于伪测量的带反馈机动目标异步航迹融合算法
CN107807278A (zh) * 2017-12-06 2018-03-16 河海大学 基于h∞扩展卡尔曼滤波的低频振荡信号参数辨识方法
CN109459705B (zh) * 2018-10-24 2021-04-27 江苏理工学院 一种抗野值鲁棒无迹卡尔曼滤波的动力电池soc估计方法
CN109724599B (zh) * 2019-03-12 2023-08-01 哈尔滨工程大学 一种抗野值的鲁棒卡尔曼滤波sins/dvl组合导航方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102323602A (zh) * 2011-05-30 2012-01-18 哈尔滨工程大学 一种基于自适应二阶卡尔曼滤波器的载波跟踪环路及其滤波方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
北斗导航系统频率校准技术研究与应用;曾亮 等;《第八届中国卫星导航学术年会论文集——S06原子钟与时频技术》;1-5 *

Also Published As

Publication number Publication date
CN110287537A (zh) 2019-09-27

Similar Documents

Publication Publication Date Title
CN110287537B (zh) 用于频标输出跳变检测的抗野值自适应卡尔曼滤波方法
US10496515B2 (en) Abnormality detection apparatus, abnormality detection method, and non-transitory computer readable medium
CN112148557B (zh) 一种性能指标实时预测方法、计算机设备及存储介质
CN109813342B (zh) 一种惯导-卫星组合导航系统的故障检测方法和系统
CN110687555B (zh) 一种导航卫星原子钟弱频率跳变在轨自主快速检测方法
JP2011522269A (ja) 異常な擬似距離測定値から無線ナビゲーション受信機ユーザを保護するための方法
CN113326616A (zh) 抗慢变量测粗差的容错卡尔曼滤波方法
CN114567288B (zh) 基于变分贝叶斯的分布协同非线性系统状态估计方法
CN108469609B (zh) 一种用于雷达目标跟踪的检测信息滤波方法
Wibowo et al. Sensor array fault detection technique using kalman filter
CN106153046B (zh) 一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法
CN115932913B (zh) 一种卫星定位伪距修正方法及装置
CN100450047C (zh) 基于模式识别的自适应检测时钟重置的方法
JP2012053040A (ja) 局所的航法衛星システム又は全地球的航法衛星システムにおけるロバストで改善された空間信号の正確性パラメータの演算
CN112965964B (zh) 一种实测飞参数据的野值检测方法、系统及计算机相关产品
CN112710306A (zh) 列车用bds与ins组合导航的自定位方法
CN110852397A (zh) 一种基于相对波动的自适应信号融合方法及系统
Ma et al. A new enabling variational inference model for approximating measurement likelihood in filtering nonlinear system
CN117798744B (zh) 一种数控机床运行状态监测方法
Bai et al. A novel progressive Gaussian approximate filter with variable step size based on a variational Bayesian approach
CN113805204B (zh) 接收机dcb的更新方法及其装置
CN117250582A (zh) 一种tdoa异常值探测方法
CN114167375A (zh) 一种绝对中位差结合最小描述长度的cfar检测方法
Gandhi Robust Kalman filters using generalized maximum likelihood-type estimators
CN117388889A (zh) 一种基于dbscan的gnss定位方法和装置

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