CN111650617B - 基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质 - Google Patents

基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质 Download PDF

Info

Publication number
CN111650617B
CN111650617B CN202010524159.9A CN202010524159A CN111650617B CN 111650617 B CN111650617 B CN 111650617B CN 202010524159 A CN202010524159 A CN 202010524159A CN 111650617 B CN111650617 B CN 111650617B
Authority
CN
China
Prior art keywords
vector
time
frequency
representing
crystal oscillator
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
CN202010524159.9A
Other languages
English (en)
Other versions
CN111650617A (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.)
State Grid Corp of China SGCC
Electric Power Research Institute of State Grid Hunan Electric Power Co Ltd
State Grid Hunan Electric Power Co Ltd
Original Assignee
State Grid Corp of China SGCC
Electric Power Research Institute of State Grid Hunan Electric Power Co Ltd
State Grid Hunan Electric Power Co Ltd
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 State Grid Corp of China SGCC, Electric Power Research Institute of State Grid Hunan Electric Power Co Ltd, State Grid Hunan Electric Power Co Ltd filed Critical State Grid Corp of China SGCC
Priority to CN202010524159.9A priority Critical patent/CN111650617B/zh
Publication of CN111650617A publication Critical patent/CN111650617A/zh
Application granted granted Critical
Publication of CN111650617B publication Critical patent/CN111650617B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/02—Frequency selective networks
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/02—Frequency selective networks
    • H03H17/0202—Two or more dimensional filters; Filters for complex signals
    • H03H2017/0205—Kalman filters

Landscapes

  • Engineering & Computer Science (AREA)
  • Signal Processing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Physics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质，本发明方法包括根据GPS秒信号、晶振信号的误差特性建立两者的频差数学模型；针对频差数学模型进行sigma采样得到k时刻的sigma点集；针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程；整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中以减少GPS较大野值对校频精度的影响。本发明基于新息加权自适应不敏卡尔曼滤波技术，将GPS信号的跳变野值和随机误差，以及晶振信号的累积误差和频率漂移进行滤波处理，最终提高GPS锁定二级频标的最终校频精度。

Description

基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、 ***及介质
技术领域
本发明涉及GPS驯服晶振信号的方法，具体涉及一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质，用于实现GPS信号对晶振频率的高精度驯服。
背景技术
随着网络技术的快速发展，***对时间基准的要求也越来越高。电力***中时间基准对事故分析、提高电能质量、电网调度操作以及电能分时计费有着重要的意义。
目前时间(频率)基准按照应用领域和性能指标可划分为一级频率标准(铯原子频标)、二级频率标准(包括铷原子频标和高稳晶体振荡器)和其它频率标准(高稳晶体振荡器意外的其他晶体振荡器)。虽然一级频标长期稳定度和准确度较高，但是其价格十分昂贵，难以应用在对成本有严格要求的民用场合；二级频标的价格相对较低，其长期稳定度和准确度却大不如前者，故难以应用于高精度场合。
现有技术利用GPS锁定二级频标(高稳晶振)原理如附图1所示：首先测量GPS接收机输出的秒脉冲信号PPS与晶振分频后输出的脉冲信号之间的相对频差，然后对频差信号进行滤波处理，最后生成控制信号，不断调节晶振频率从而输出高准确度和高稳定度的频率输出。该技术在一定程度上解决了上述问题，即利用较低的成本实现高长期稳定度和平均准确度的频率标准。但是在实际应用中由于GPS的1PPS信号存在随机误差和跳变野值、晶振因老化和温度等特性产生的频率漂移为频率的校准造成了一定的困难。
为了解决以上问题，有的学者提出采用多次采样取平均的方法，对测量的频差数据进行处理，计算简单但是校频时间较长。有学者通过建立晶振信号和GPS信号之间的频差数学模型，利用最小二乘法对其误差进行计算和补偿，但该方法计算较为复杂且没有考虑GPS信号的野值对校频精度的影响。还有学者提出利用卡尔曼滤波算法来消除频差信号中的噪声，但是现有利用卡尔曼滤波算法来消除频差信号中的噪声的方式没有考虑晶振频率漂移的影响，导致校频精度有限。
发明内容
本发明要解决的技术问题：针对现有技术的上述问题，提供一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质，本发明基于新息加权自适应不敏卡尔曼滤波技术，将GPS信号的跳变野值和随机误差，以及晶振信号的累积误差和频率漂移进行滤波处理，最终提高GPS锁定二级频标的最终校频精度。
为了解决上述技术问题，本发明采用的技术方案为：
一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法，包括：
1)根据GPS秒信号、晶振信号的误差特性建立两者的频差数学模型；
2)针对频差数学模型进行sigma采样得到k时刻的sigma点集；
3)针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程；
4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中以减少GPS较大野值对校频精度的影响。
可选地，步骤1)的详细步骤包括：
1.1)获取GPS秒信号序列Y，GPS秒信号序列Y中的任意k时刻的GPS秒信号的世界统一时间UTC的函数表达式如式(1)所示；
Yk＝k-εk (1)
上式中，Yk表示GPS秒信号序列Y中的k时刻的GPS秒信号的世界统一时间UTC，εk为k时刻的GPS秒信号的误差，且满足ε～N(0,σ2)；
获取晶振信号序列Y′，晶振信号序列Y′中的任意k时刻的晶振信号的世界统一时间UTC的函数表达式如式(2)所示；
Y′k＝k+a+bk+ck2 (2)
上式中，Y′k表示晶振信号序列Y′中的k时刻的晶振信号的世界统一时间UTC，a为初始相位误差，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数；
1.2)根据GPS秒信号序列Y、晶振信号序列Y′做差得到频差序列，且频差序列中任意k时刻的频差信号的函数表达式如下式所示；
Xk＝Y′k-Yk＝a+bk+ck2+εk (3)
上式中，Xk表示频差序列中k时刻的频差信号。
可选地，步骤2)的详细步骤包括：
2.1)基于频差数学模型建立频差的预测方程：
X'k＝X'k-1+2ck+b-c+ω(k) (4)
上式中，X'k表示k时刻的频差信号的预测值，X'k-1表示k-1时刻的频差信号的预测值，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数，ω(k)表示k时刻的过程噪声；
2.2)令u(k)＝2ck+b-c表示输入向量，建立频差的状态方程和量测方程：
Figure GDA0003976378670000031
上式中，x(k+1)表示k+1时刻的频差状态向量，x(k)表示k时刻的频差状态向量，u(k)表示输入向量，z(k+1)表示k+1时刻的量测值；ω(k)表示k时刻的过程噪声向量，其协方差矩阵为Q；γ(k)表示k时刻的量测噪声向量，其协方差矩阵为R；
2.3)对频差状态向量进行扩维处理xa＝[xT ωT γT]，其中xa表示扩维后的状态向量，xT表示频差状态向量的转置，ωT表示过程噪声向量的转置，γT表示量测噪声向量的转置；确定频差状态向量的初始状态：
Figure GDA0003976378670000032
上式中，
Figure GDA0003976378670000033
表示频差状态向量初值，E(x0)表示频差向量初值的均值，x0表示频差向量初值，P0表示频差向量的协方差向量初值，
Figure GDA0003976378670000034
表示频差向量初值的协方差；
确定频差状态向量的扩维处理后的状态向量初始状态：
Figure GDA0003976378670000035
上式中，
Figure GDA0003976378670000036
表示扩维后的状态向量初值，
Figure GDA0003976378670000037
表示扩维后的状态向量均值，
Figure GDA0003976378670000038
表示表示扩维后的状态向量，
Figure GDA0003976378670000039
表示频差状态向量初值，
Figure GDA00039763786700000310
表示扩维后的状态向量协方差初值，
Figure GDA00039763786700000311
表示求取扩维后状态向量初值的协方差，P0表示频差向量的协方差向量初值，Q表示k时刻的过程噪声向量ω(k)的协方差矩阵，R表示k时刻的量测噪声向量γ(k)的协方差矩阵；
2.4)针对频差状态向量采用对称sigma采样，得到k时刻状态的sigma点集
Figure GDA00039763786700000312
i＝1,…,L，其中
Figure GDA00039763786700000313
表示第i个维数为a的sigma点，L为sigma点采样个数，L＝2a+1，此时状态维数a＝n+q+m，其中n表示频差状态向量的维数，q表示过程噪声向量的维数，m表示量测噪声向量的维数；且sigma点集
Figure GDA0003976378670000041
的函数表达式如下式所示：
Figure GDA0003976378670000042
上式中，
Figure GDA0003976378670000043
表示状态向量的均值，Pxx表示状态向量的协方差，γ为系数，
Figure GDA0003976378670000044
a为状态维数，λ＝α2(a+κ)-a，κ为比例系数，γPxx的平方根矩阵的第i行或列记为
Figure GDA0003976378670000045
2.5)确定k时刻状态的sigma点集中各sigma点对应的权重：
Figure GDA0003976378670000046
上式中，Wi m表示均值权值，a为状态维数，λ＝α2(a+κ)-a，κ为比例系数；
Figure GDA0003976378670000047
上式中，Wi c表示协方差权值，Wi m表示均值权值，a为状态维数，λ＝α2(a+κ)-a，κ为比例系数，α为用于确定状态向量均值
Figure GDA0003976378670000048
周围sigma点的分布的系数，1e-4≤α<1；β为常数系数。
可选地，步骤3)中针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程时，以
Figure GDA0003976378670000049
表示粒子
Figure GDA00039763786700000410
的前n维组成的列向量；
Figure GDA00039763786700000411
表示粒子
Figure GDA00039763786700000412
的n+1维到n+q维组成的列向量；
Figure GDA00039763786700000413
表示粒子
Figure GDA00039763786700000414
的n+q+1维到n+q+m维组成的列向量，建立任意k时刻的预测过程的状态方程如式(12)～(14)所示，建立任意k时刻的更新方程如式(15)～(21)所示；
Figure GDA00039763786700000415
Figure GDA00039763786700000416
Figure GDA00039763786700000417
式(12)～(14)中，
Figure GDA0003976378670000051
表示k+1时刻的sigma点预测值，
Figure GDA0003976378670000052
表示k时刻状态向量sigma点的估计值，u(k)表示输入向量，
Figure GDA0003976378670000053
表示k时刻过程噪声向量的sigma点估计值；
Figure GDA0003976378670000054
表示k+1时刻状态向量预测值，Wi m表示均值权值，L为sigma点采样个数；P(k+1,k)表示k+1时刻的协方差预测值；
Figure GDA0003976378670000055
Figure GDA0003976378670000056
Figure GDA0003976378670000057
Figure GDA0003976378670000058
Figure GDA0003976378670000059
Figure GDA00039763786700000510
Figure GDA00039763786700000511
式(15)～(21)中，zi(k+1,k)表示k+1时刻的量测向量估计值，
Figure GDA00039763786700000512
表示k时刻状态向量sigma点的估计值，u(k)为输入向量，
Figure GDA00039763786700000513
表示k+1时刻量测噪声向量的sigma点估计值；
Figure GDA00039763786700000514
表示k+1时刻量测向量的预测值，Wi m表示均值权值，L为sigma点采样个数；Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差，αk+1为自适应调整因子，Wi c表示协方差权值；Pxz(k+1,k)表示k+1时刻状态向量和量测向量估计偏差的互协方差，
Figure GDA00039763786700000515
表示k+1时刻的sigma点预测值，
Figure GDA00039763786700000516
表示k+1时刻的状态向量预测值；K(k+1)表示k+1时刻的滤波增益；
Figure GDA00039763786700000517
表示k+1时刻状态向量的最优估计值，z(k+1)表示k+1时刻的量测值；P(k+1,k+1)表示k+1时刻协方差最优估计值。
可选地，自适应调整因子αk+1的计算函数表达式如下式所示：
Figure GDA0003976378670000061
上式中，tr(Pzz(k+1,k))表示矩阵Pzz(k+1,k)的迹，tr(ΔZk+1ΔZk+1 T)表示矩阵ΔZk+1ΔZk+1 T的迹，Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差，ΔZk+1为量测时间更新残差值，
Figure GDA0003976378670000062
可选地，步骤4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中的函数表达式如式(23)～(27)所示：
Figure GDA0003976378670000063
Figure GDA0003976378670000064
Figure GDA0003976378670000065
Figure GDA0003976378670000066
Dk+1＝P(k+1,k)+Rk+1 (27)
式(23)～(27)中，
Figure GDA0003976378670000067
表示k+1时刻状态向量的最优估计值，
Figure GDA0003976378670000068
表示k+1时刻状态向量的预测值，K(k+1)表示k+1时刻的滤波增益值，φk+1表示修正函数，τk+1表示衡量k+1时刻量测误差大小的变量，ek+1表示k+1时刻的量测误差；z(k+1)表示k+1时刻的量测值；
Figure GDA0003976378670000069
表示k+1时刻量测向量预测值；τk+1服从自由度为m的χ2分布，当出现野值时τk+1增大、φk+1(·)减小；Ck+1为选取的门限值或常数序列，λk+1为
Figure GDA00039763786700000610
矩阵的最大特征值，Kk+1表示k+1时刻的卡尔曼滤波增益，Dk+1表示k+1时刻衡量量测误差大小的权矩阵，P(k+1,k)表示k+1时刻的协方差预测值，Rk+1表示k+1时刻量测噪声方差。
此外，本发明还提供一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服***，包括计算机设备，该计算机设备被编程或配置以执行所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的步骤。
此外，本发明还提供一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服***，包括计算机设备，该计算机设备的存储器上存储有被编程或配置以执行所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
此外，本发明还提供一种带有GPS的智能终端，包括GPS模块、微处理器和存储器，该微处理器被编程或配置以执行所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的步骤，或该存储器上存储有被编程或配置以执行所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
此外，本发明还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
和现有技术相比，本发明具有下述优点：本发明提供了一种对GPS锁定晶振频率技术中的频差信号进行处理的新息序列加权的自适应不敏卡尔曼滤波方法，该方法主要包括四个部分：一是根据GPS的1PPS信号和晶振信号的误差特性建立两者的频差数学模型；二是根据频差数学模型确定***的量测和预测方程，进行sigma采样；三是确定自适应不敏卡尔曼算法的预测方程及更新方程；四是整合新息序列加权到上述算法之中，以减少GPS较大野值对校频精度的影响。通过以上四步，可以实现新息序列加权的自适应不敏卡尔曼滤波，本具有以下技术效果：1、采用不敏卡尔曼滤波算法对GPS锁定晶振频率过程中的频差信号进行处理，在对GPS信号的随机误差进行处理的同时也考虑了晶振由于老化和温度影响而产生的非线性误差；2、将新息序列加权和自适应滤波技术引入不敏卡尔曼滤波过程，可以减少GPS信号中较大跳变野值以及***噪声和量测噪声的协方差矩阵不断变化对于滤波结果的影响。
附图说明
图1为现有技术中GPS锁定晶振频率原理图。
图2为本发明实施例方法的基本流程示意图。
具体实施方式
如图2所示，本实施例基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法包括：
1)根据GPS秒信号(1PPS信号)、晶振信号的误差特性建立两者的频差数学模型；
2)针对频差数学模型进行sigma采样得到k时刻的sigma点集；
3)针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程；
4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中以减少GPS较大野值对校频精度的影响。
获取GPS秒信号序列Y时，GPS秒信号序列Y(大小为n)中的GPS秒信号的世界统一时间UTC可分别表示为：
1-ε1，2-ε2，…，k-εk，…，n-εn
获取晶振信号序列Y′时，晶振信号序列Y′(大小为n)中的晶振信号的世界统一时间UTC可分别表示为：
1+a+b+c，1+a+2b+4c，…，k+a+bk+ck2，…，n+a+bn+cn2
其中，a为初始相位误差，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数；因此，本实施例中步骤1)的详细步骤包括：
1.1)获取GPS秒信号序列Y，GPS秒信号序列Y中的任意k时刻的GPS秒信号的世界统一时间UTC的函数表达式如式(1)所示；
Yk＝k-εk (1)
上式中，Yk表示GPS秒信号序列Y中的k时刻的GPS秒信号的世界统一时间UTC，εk为k时刻的GPS秒信号的误差，且满足ε～N(0,σ2)；
获取晶振信号序列Y′，晶振信号序列Y′中的任意k时刻的晶振信号的世界统一时间UTC的函数表达式如式(2)所示；
Y′k＝k+a+bk+ck2 (2)
上式中，Y′k表示晶振信号序列Y′中的k时刻的晶振信号的世界统一时间UTC，a为初始相位误差，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数；
1.2)根据GPS秒信号序列Y、晶振信号序列Y′做差得到频差序列，且频差序列中任意k时刻的频差信号的函数表达式如下式所示；
Xk＝Y′k-Yk＝a+bk+ck2+εk (3)
上式中，Xk表示频差序列中k时刻的频差信号。
本实施例中，步骤2)的详细步骤包括：
2.1)基于频差数学模型建立频差的预测方程：
X'k＝X'k-1+2ck+b-c+ω(k) (4)
上式中，X'k表示k时刻的频差信号的预测值，X'k-1表示k-1时刻的频差信号的预测值，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数，ω(k)表示k时刻的过程噪声；
2.2)令u(k)＝2ck+b-c表示输入向量，建立频差的状态方程和量测方程：
Figure GDA0003976378670000081
上式中，x(k+1)表示k+1时刻的频差状态向量，x(k)表示k时刻的频差状态向量，u(k)表示输入向量，z(k+1)表示k+1时刻的量测值；ω(k)表示k时刻的过程噪声向量，其协方差矩阵为Q；γ(k)表示k时刻的量测噪声向量，其协方差矩阵为R；
2.3)由于***具有噪声项，因此需要对状态向量进行扩维处理，对频差状态向量进行扩维处理xa＝[xT ωT γT]，其中xa表示扩维后的状态向量，xT表示频差状态向量的转置，ωT表示过程噪声向量的转置，γT表示量测噪声向量的转置；确定频差状态向量的初始状态：
Figure GDA0003976378670000091
上式中，
Figure GDA0003976378670000092
表示频差状态向量初值，E(x0)表示频差向量初值的均值，x0表示频差向量初值，P0表示频差向量的协方差向量初值，
Figure GDA0003976378670000093
表示频差向量初值的协方差；
确定频差状态向量的扩维处理后的状态向量初始状态：
Figure GDA0003976378670000094
上式中，
Figure GDA0003976378670000095
表示扩维后的状态向量初值，
Figure GDA0003976378670000096
表示扩维后的状态向量均值，
Figure GDA0003976378670000097
表示表示扩维后的状态向量，
Figure GDA0003976378670000098
表示频差状态向量初值，
Figure GDA0003976378670000099
表示扩维后的状态向量协方差初值，
Figure GDA00039763786700000910
表示求取扩维后状态向量初值的协方差，P0表示频差向量的协方差向量初值，Q表示k时刻的过程噪声向量ω(k)的协方差矩阵，R表示k时刻的量测噪声向量γ(k)的协方差矩阵；
2.4)针对频差状态向量采用对称sigma采样，得到k时刻状态的sigma点集
Figure GDA00039763786700000911
i＝1,…,L，其中
Figure GDA00039763786700000912
表示第i个维数为a的sigma点，L为sigma点采样个数，L＝2a+1，此时状态维数a＝n+q+m，其中n表示频差状态向量的维数，q表示过程噪声向量的维数，m表示量测噪声向量的维数；且sigma点集
Figure GDA00039763786700000913
的函数表达式如下式所示：
Figure GDA0003976378670000101
上式中，
Figure GDA0003976378670000102
表示状态向量的均值，Pxx表示状态向量的协方差，γ为系数，
Figure GDA0003976378670000103
a为状态维数，λ＝α2(a+κ)-a，κ为比例系数(比例系数κ通常可取值为0或者3-a，可用于调整sigma点和x均值
Figure GDA0003976378670000104
之间的距离)，γPxx的平方根矩阵的第i行或列记为
Figure GDA0003976378670000105
2.5)确定k时刻状态的sigma点集中各sigma点对应的权重：
Figure GDA0003976378670000106
上式中，Wi m表示均值权值，a为状态维数，λ＝α2(a+κ)-a，κ为比例系数；
Figure GDA0003976378670000107
上式中，Wi c表示协方差权值，Wi m表示均值权值，a为状态维数，λ＝α2(a+κ)-a，κ为比例系数，α为用于确定状态向量均值
Figure GDA0003976378670000108
周围sigma点的分布的系数，1e-4≤α<1；β为常数系数。
本实施例步骤3)中针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程时，以
Figure GDA0003976378670000109
表示粒子
Figure GDA00039763786700001010
的前n维组成的列向量；
Figure GDA00039763786700001011
表示粒子
Figure GDA00039763786700001012
的n+1维到n+q维组成的列向量；
Figure GDA00039763786700001013
表示粒子
Figure GDA00039763786700001014
的n+q+1维到n+q+m维组成的列向量，建立任意k时刻的预测过程的状态方程如式(12)～(14)所示，建立任意k时刻的更新方程如式(15)～(21)所示；
Figure GDA00039763786700001015
Figure GDA00039763786700001016
Figure GDA00039763786700001017
式(12)～(14)中，
Figure GDA0003976378670000111
表示k+1时刻的sigma点预测值，
Figure GDA0003976378670000112
表示k时刻状态向量sigma点的估计值，u(k)表示输入向量，
Figure GDA0003976378670000113
表示k时刻过程噪声向量的sigma点估计值；
Figure GDA0003976378670000114
表示k+1时刻状态向量预测值，Wi m表示均值权值，L为sigma点采样个数；P(k+1,k)表示k+1时刻的协方差预测值；
Figure GDA0003976378670000115
Figure GDA0003976378670000116
Figure GDA0003976378670000117
Figure GDA0003976378670000118
Figure GDA0003976378670000119
Figure GDA00039763786700001110
Figure GDA00039763786700001111
式(15)～(21)中，zi(k+1,k)表示k+1时刻的量测向量估计值，
Figure GDA00039763786700001112
表示k时刻状态向量sigma点的估计值，u(k)为输入向量，
Figure GDA00039763786700001113
表示k+1时刻量测噪声向量的sigma点估计值；
Figure GDA00039763786700001114
表示k+1时刻量测向量的预测值，Wi m表示均值权值，L为sigma点采样个数；Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差，αk+1为自适应调整因子，Wi c表示协方差权值；Pxz(k+1,k)表示k+1时刻状态向量和量测向量估计偏差的互协方差，
Figure GDA00039763786700001115
表示k+1时刻的sigma点预测值，
Figure GDA00039763786700001116
表示k+1时刻的状态向量预测值；K(k+1)表示k+1时刻的滤波增益；
Figure GDA00039763786700001117
表示k+1时刻状态向量的最优估计值，z(k+1)表示k+1时刻的量测值；P(k+1,k+1)表示k+1时刻协方差最优估计值。
本实施例中，自适应调整因子αk+1的计算函数表达式如下式所示：
Figure GDA0003976378670000121
上式中，tr(Pzz(k+1,k))表示矩阵Pzz(k+1,k)的迹，tr(ΔZk+1ΔZk+1 T)表示矩阵ΔZk+1ΔZk+1 T的迹，Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差，ΔZk+1为量测时间更新残差值，
Figure GDA0003976378670000122
在实际应用中，由于GPS信号受噪声干扰还会产生较大的跳变野值，带有野值的频差信号将会使滤波结果发生偏移甚至发散。为了解决上述不敏卡尔曼算法的抗野值问题，有学者将新息序列加权的方法引入卡尔曼滤波来消除野值的影响。当观测数据不包含野值时，滤波算法能够正常运行，对观测噪声进行滤除；当观测数据包含野值时，能够克服或减少野值对于滤波结果的影响，使滤波结果更加准确。本实施例中，步骤4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中的函数表达式如式(23)～(27)所示：
Figure GDA0003976378670000123
Figure GDA0003976378670000124
Figure GDA0003976378670000125
Figure GDA0003976378670000126
Dk+1＝P(k+1,k)+Rk+1 (27)
式(23)～(27)中，
Figure GDA0003976378670000127
表示k+1时刻状态向量的最优估计值，
Figure GDA0003976378670000128
表示k+1时刻状态向量的预测值，K(k+1)表示k+1时刻的滤波增益值，φk+1表示修正函数，τk+1表示衡量k+1时刻量测误差大小的变量，ek+1表示k+1时刻的量测误差；z(k+1)表示k+1时刻的量测值；
Figure GDA0003976378670000129
表示k+1时刻量测向量预测值；τk+1服从自由度为m的χ2分布，当出现野值时τk+1增大、φk+1(·)减小；Ck+1为选取的门限值或常数序列，λk+1为
Figure GDA00039763786700001210
矩阵的最大特征值，Kk+1表示k+1时刻的卡尔曼滤波增益，Dk+1表示k+1时刻衡量量测误差大小的权矩阵，P(k+1,k)表示k+1时刻的协方差预测值，Rk+1表示k+1时刻量测噪声方差。
综上所述，目前应用于低成本高精度频率标准的主要技术为GPS锁定二级频标(晶振)技术，为减少GPS锁定二级频标技术中GPS信号的随机误差和跳变野值以及晶振信号的非线性频率偏移对于频差信号的影响，本实施例方法提供了一种对GPS锁定晶振频率技术中的频差信号进行处理的新息序列加权的自适应不敏卡尔曼滤波方法，该方法主要包括四个部分：一是根据GPS的1PPS信号和晶振信号的误差特性建立两者的频差数学模型；二是根据频差数学模型确定***的量测和预测方程，进行sigma采样；三是确定自适应不敏卡尔曼算法的预测方程及更新方程；四是整合新息序列加权到上述算法之中，以减少GPS较大野值对校频精度的影响。通过以上四步，可以实现新息序列加权的自适应不敏卡尔曼滤波，本具有以下技术效果：1、采用不敏卡尔曼滤波算法对GPS锁定晶振频率过程中的频差信号进行处理，在对GPS信号的随机误差进行处理的同时也考虑了晶振由于老化和温度影响而产生的非线性误差；2、将新息序列加权和自适应滤波技术引入不敏卡尔曼滤波过程，可以减少GPS信号中较大跳变野值以及***噪声和量测噪声的协方差矩阵不断变化对于滤波结果的影响。
此外，本实施例还提供一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服***，包括计算机设备，其特征在于，该计算机设备被编程或配置以执行权利要求1～5中任意一项所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的步骤。
此外，本实施例还提供一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服***，包括计算机设备，该计算机设备的存储器上存储有被编程或配置以执行前述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
此外，本实施例还提供一种带有GPS的智能终端，包括GPS模块、微处理器和存储器，该微处理器被编程或配置以执行前述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的步骤，或该存储器上存储有被编程或配置以执行前述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
此外，本实施例还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行前述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
本领域内的技术人员应明白，本申请的实施例可提供为方法、***、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可读存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。本申请是参照根据本申请实施例的方法、设备(***)、和计算机程序产品的流程图和/的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (9)

1.一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法，其特征在于包括：
1)根据GPS秒信号、晶振信号的误差特性建立两者的频差数学模型；
2)针对频差数学模型进行sigma采样得到k时刻的sigma点集；
3)针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程；
4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中以减少GPS较大野值对校频精度的影响；
步骤3)中针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程时，以
Figure FDA0003976378660000011
表示粒子
Figure FDA0003976378660000012
的前n维组成的列向量；
Figure FDA0003976378660000013
表示粒子
Figure FDA0003976378660000014
的n+1维到n+q维组成的列向量；
Figure FDA0003976378660000015
表示粒子
Figure FDA0003976378660000016
的n+q+1维到n+q+m维组成的列向量，建立任意k时刻的预测过程的状态方程如式(12)～(14)所示，建立任意k时刻的更新方程如式(15)～(21)所示；
Figure FDA0003976378660000017
Figure FDA0003976378660000018
Figure FDA0003976378660000019
式(12)～(14)中，
Figure FDA00039763786600000110
表示k+1时刻的sigma点预测值，
Figure FDA00039763786600000111
表示k时刻状态向量sigma点的估计值，u(k)表示输入向量，
Figure FDA00039763786600000112
表示k时刻过程噪声向量的sigma点估计值；
Figure FDA00039763786600000113
表示k+1时刻状态向量预测值，Wi m表示均值权值，L为sigma点采样个数；P(k+1,k)表示k+1时刻的协方差预测值；
Figure FDA00039763786600000114
Figure FDA00039763786600000115
Figure FDA00039763786600000116
Figure FDA00039763786600000117
Figure FDA00039763786600000118
Figure FDA0003976378660000021
Figure FDA0003976378660000022
式(15)～(21)中，zi(k+1,k)表示k+1时刻的量测向量估计值，
Figure FDA0003976378660000023
表示k时刻状态向量sigma点的估计值，u(k)为输入向量，
Figure FDA0003976378660000024
表示k+1时刻量测噪声向量的sigma点估计值；
Figure FDA0003976378660000025
表示k+1时刻量测向量的预测值，Wi m表示均值权值，L为sigma点采样个数；Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差，αk+1为自适应调整因子，Wi c表示协方差权值；Pxz(k+1,k)表示k+1时刻状态向量和量测向量估计偏差的互协方差，
Figure FDA0003976378660000026
表示k+1时刻的sigma点预测值，
Figure FDA0003976378660000027
表示k+1时刻的状态向量预测值；K(k+1)表示k+1时刻的滤波增益；
Figure FDA0003976378660000028
表示k+1时刻状态向量的最优估计值，z(k+1)表示k+1时刻的量测值；P(k+1,k+1)表示k+1时刻协方差最优估计值。
2.根据权利要求1所述的基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法，其特征在于，步骤1)的详细步骤包括：
1.1)获取GPS秒信号序列Y，GPS秒信号序列Y中的任意k时刻的GPS秒信号的世界统一时间UTC的函数表达式如式(1)所示；
Yk＝k-εk (1)
上式中，Yk表示GPS秒信号序列Y中的k时刻的GPS秒信号的世界统一时间UTC，εk为k时刻的GPS秒信号的误差，且满足ε～N(0,σ2)；
获取晶振信号序列Y′，晶振信号序列Y′中的任意k时刻的晶振信号的世界统一时间UTC的函数表达式如式(2)所示；
Y′k＝k+a+bk+ck2 (2)
上式中，Y′k表示晶振信号序列Y′中的k时刻的晶振信号的世界统一时间UTC，a为初始相位误差，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数；
1.2)根据GPS秒信号序列Y、晶振信号序列Y′做差得到频差序列，且频差序列中任意k时刻的频差信号的函数表达式如下式所示；
Xk＝Y′k-Yk＝a+bk+ck2+εk (3)
上式中，Xk表示频差序列中k时刻的频差信号。
3.根据权利要求2所述的基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法，其特征在于，步骤2)的详细步骤包括：
2.1)基于频差数学模型建立频差的预测方程：
X'k＝X'k-1+2ck+b-c+ω(k) (4)
上式中，X'k表示k时刻的频差信号的预测值，X'k-1表示k-1时刻的频差信号的预测值，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数，ω(k)表示k时刻的过程噪声；
2.2)令u(k)＝2ck+b-c表示输入向量，建立频差的状态方程和量测方程：
Figure FDA0003976378660000031
上式中，x(k+1)表示k+1时刻的频差状态向量，x(k)表示k时刻的频差状态向量，u(k)表示输入向量，z(k+1)表示k+1时刻的量测值；ω(k)表示k时刻的过程噪声向量，其协方差矩阵为Q；γ(k)表示k时刻的量测噪声向量，其协方差矩阵为R；
2.3)对频差状态向量进行扩维处理xa＝[xT ωT γT]，其中xa表示扩维后的状态向量，xT表示频差状态向量的转置，ωT表示过程噪声向量的转置，γT表示量测噪声向量的转置；确定频差状态向量的初始状态：
Figure FDA0003976378660000032
上式中，
Figure FDA0003976378660000033
表示频差状态向量初值，E(x0)表示频差向量初值的均值，x0表示频差向量初值，P0表示频差向量的协方差向量初值，
Figure FDA0003976378660000034
表示频差向量初值的协方差；
确定频差状态向量的扩维处理后的状态向量初始状态：
Figure FDA0003976378660000035
上式中，
Figure FDA0003976378660000041
表示扩维后的状态向量初值，
Figure FDA0003976378660000042
表示扩维后的状态向量均值，
Figure FDA0003976378660000043
表示表示扩维后的状态向量，
Figure FDA0003976378660000044
表示频差状态向量初值，
Figure FDA0003976378660000045
表示扩维后的状态向量协方差初值，
Figure FDA0003976378660000046
表示求取扩维后状态向量初值的协方差，P0表示频差向量的协方差向量初值，Q表示k时刻的过程噪声向量ω(k)的协方差矩阵，R表示k时刻的量测噪声向量γ(k)的协方差矩阵；
2.4)针对频差状态向量采用对称sigma采样，得到k时刻状态的sigma点集
Figure FDA0003976378660000047
i＝1,…,L，其中
Figure FDA0003976378660000048
表示第i个维数为a的sigma点，L为sigma点采样个数，L＝2a+1，此时状态维数a＝n+q+m，其中n表示频差状态向量的维数，q表示过程噪声向量的维数，m表示量测噪声向量的维数；且sigma点集
Figure FDA0003976378660000049
的函数表达式如下式所示：
Figure FDA00039763786600000410
上式中，
Figure FDA00039763786600000411
表示状态向量的均值，Pxx表示状态向量的协方差，γ为系数，
Figure FDA00039763786600000412
a为状态维数，λ＝α2(a+κ)-a，κ为比例系数，γPxx的平方根矩阵的第i行或列记为
Figure FDA00039763786600000413
2.5)确定k时刻状态的sigma点集中各sigma点对应的权重：
Figure FDA00039763786600000414
上式中，Wi m表示均值权值，a为状态维数，λ＝α2(a+κ)-a，κ为比例系数；
Figure FDA00039763786600000415
上式中，Wi c表示协方差权值，Wi m表示均值权值，a为状态维数，λ＝α2(a+κ)-a，κ为比例系数，α为用于确定状态向量均值
Figure FDA00039763786600000416
周围sigma点的分布的系数，1e-4≤α<1；β为常数系数。
4.根据权利要求3所述的基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法，其特征在于，自适应调整因子αk+1的计算函数表达式如下式所示：
Figure FDA0003976378660000051
上式中，tr(Pzz(k+1,k))表示矩阵Pzz(k+1,k)的迹，tr(ΔZk+1ΔZk+1 T)表示矩阵ΔZk+1ΔZk+1 T的迹，Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差，ΔZk+1为量测时间更新残差值，
Figure FDA0003976378660000052
5.根据权利要求3所述的基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法，其特征在于，步骤4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中的函数表达式如式(23)～(27)所示：
Figure FDA0003976378660000053
Figure FDA0003976378660000054
Figure FDA0003976378660000055
Figure FDA0003976378660000056
Dk+1＝P(k+1,k)+Rk+1 (27)
式(23)～(27)中，
Figure FDA0003976378660000057
表示k+1时刻状态向量的最优估计值，
Figure FDA0003976378660000058
表示k+1时刻状态向量的预测值，K(k+1)表示k+1时刻的滤波增益值，φk+1表示修正函数，τk+1表示衡量k+1时刻量测误差大小的变量，ek+1表示k+1时刻的量测误差；z(k+1)表示k+1时刻的量测值；
Figure FDA0003976378660000059
表示k+1时刻量测向量预测值；τk+1服从自由度为m的χ2分布，当出现野值时τk+1增大、φk+1(·)减小；Ck+1为选取的门限值或常数序列，λk+1为
Figure FDA00039763786600000510
矩阵的最大特征值，Kk+1表示k+1时刻的卡尔曼滤波增益，Dk+1表示k+1时刻衡量量测误差大小的权矩阵，P(k+1,k)表示k+1时刻的协方差预测值，Rk+1表示k+1时刻量测噪声方差。
6.一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服***，包括计算机设备，其特征在于，该计算机设备被编程或配置以执行权利要求1～5中任意一项所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的步骤。
7.一种基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服***，包括计算机设备，其特征在于，该计算机设备的存储器上存储有被编程或配置以执行权利要求1～5中任意一项所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
8.一种带有GPS的智能终端，包括GPS模块、微处理器和存储器，其特征在于，该微处理器被编程或配置以执行权利要求1～5中任意一项所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的步骤，或该存储器上存储有被编程或配置以执行权利要求1～5中任意一项所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
9.一种计算机可读存储介质，其特征在于，该计算机可读存储介质上存储有被编程或配置以执行权利要求1～5中任意一项所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
CN202010524159.9A 2020-06-10 2020-06-10 基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质 Active CN111650617B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010524159.9A CN111650617B (zh) 2020-06-10 2020-06-10 基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010524159.9A CN111650617B (zh) 2020-06-10 2020-06-10 基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质

Publications (2)

Publication Number Publication Date
CN111650617A CN111650617A (zh) 2020-09-11
CN111650617B true CN111650617B (zh) 2023-03-03

Family

ID=72343261

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010524159.9A Active CN111650617B (zh) 2020-06-10 2020-06-10 基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质

Country Status (1)

Country Link
CN (1) CN111650617B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112578415B (zh) * 2020-11-06 2023-10-13 中国科学院国家空间科学中心 一种基于自适应滤波器的数字频率锁定方法及环路
CN112380774B (zh) * 2020-11-23 2022-04-15 青岛柯锐思德电子科技有限公司 一种基于残差回声状态网络的动态建模方法及***
CN112713881B (zh) * 2020-12-10 2022-11-01 国网四川省电力公司电力科学研究院 一种基于边缘计算的同步时钟维持***与方法
CN113419414B (zh) * 2021-07-13 2023-03-21 贵州省计量测试院 一种具有gnss驯服和间隔时间保持能力的标准计时***
CN113959433B (zh) * 2021-09-16 2023-12-08 南方电网数字平台科技(广东)有限公司 一种组合导航方法及装置
CN118100914B (zh) * 2024-04-18 2024-07-05 西安乾景防务技术有限公司 基于人工智能的双频率合成器控制方法及***

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6711230B1 (en) * 2002-09-27 2004-03-23 Nortel Networks Limited Reference timing signal oscillator with frequency stability
US8643444B2 (en) * 2012-06-04 2014-02-04 Broadcom Corporation Common reference crystal systems
KR20140070722A (ko) * 2012-11-26 2014-06-11 한국전자통신연구원 다중속도시스템 결합 장치 및 운용 방법
JP2018165618A (ja) * 2017-03-28 2018-10-25 セイコーエプソン株式会社 信号処理装置、検出装置、物理量測定装置、電子機器及び移動体
CN110336279B (zh) * 2019-07-17 2020-11-20 国网湖南省电力有限公司 基于并网变换器的电力***振荡自适应抑制方法、***及介质

Also Published As

Publication number Publication date
CN111650617A (zh) 2020-09-11

Similar Documents

Publication Publication Date Title
CN111650617B (zh) 基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质
CN106291645B (zh) 适于高维gnss/ins深耦合的容积卡尔曼滤波方法
US4899117A (en) High accuracy frequency standard and clock system
CN103033186B (zh) 一种用于水下滑翔器的高精度组合导航定位方法
Weiss et al. AT2, a new time scale algorithm: AT1 plus frequency variance
CN101562451B (zh) 二级频标的精密驯服保持方法
CN110554597B (zh) 基于Vondark-Cepek滤波的氢铯时间尺度融合方法
CN109508510B (zh) 一种基于改进的卡尔曼滤波的铷原子钟参数估计算法
WO1995034850B1 (en) Assured-integrity monitored-extrapolation navigation apparatus
CN109581856B (zh) 一种基于高性能晶振频率校准的对时守时方法
US10033390B2 (en) Systems and methods for clock synchronization in a data acquisition system
CN111694040B (zh) 一种卫星/惯性组合导航***定位方法与装置
WO1985001807A1 (en) Selective parametric self-calibrating control system
US9456431B2 (en) Precise temperature and timebase ppm error estimation using multiple timebases
CN107607977B (zh) 一种基于最小偏度单形采样的自适应ukf组合导航方法
Tang et al. Complexity reduction of the Kalman filter-based tracking loops in GNSS receivers
CN104678343A (zh) 一种波形发生器频响特性校准方法、装置及***
Narasimhappa et al. An innovation based random weighting estimation mechanism for denoising fiber optic gyro drift signal
CN114567288A (zh) 基于变分贝叶斯的分布协同非线性***状态估计方法
CN111711446B (zh) 一种利用gps信号驯服晶振频率的方法、***及介质
CN112505386B (zh) 一种用于检定直流充电桩电流值的方法及***
CN109655057B (zh) 一种六推无人机加速器测量值的滤波优化方法及其***
Sarunic Development of GPS receiver Kalman Filter algorithms for stationary, low-dynamics, and high-dynamics applications
Griffin et al. Drift correction for active hydrogen MASERs
D'Appolito The evaluation of Kalman filter designs for multisensor integrated navigation systems

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