CN111650617B - 基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质 - Google Patents
基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法、***及介质 Download PDFInfo
- 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
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表示输入向量,建立频差的状态方程和量测方程:
上式中,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表示量测噪声向量的转置;确定频差状态向量的初始状态:
确定频差状态向量的扩维处理后的状态向量初始状态:
上式中,表示扩维后的状态向量初值,表示扩维后的状态向量均值,表示表示扩维后的状态向量,表示频差状态向量初值,表示扩维后的状态向量协方差初值,表示求取扩维后状态向量初值的协方差,P0表示频差向量的协方差向量初值,Q表示k时刻的过程噪声向量ω(k)的协方差矩阵,R表示k时刻的量测噪声向量γ(k)的协方差矩阵;
2.4)针对频差状态向量采用对称sigma采样,得到k时刻状态的sigma点集i=1,…,L,其中表示第i个维数为a的sigma点,L为sigma点采样个数,L=2a+1,此时状态维数a=n+q+m,其中n表示频差状态向量的维数,q表示过程噪声向量的维数,m表示量测噪声向量的维数;且sigma点集的函数表达式如下式所示:
2.5)确定k时刻状态的sigma点集中各sigma点对应的权重:
上式中,Wi m表示均值权值,a为状态维数,λ=α2(a+κ)-a,κ为比例系数;
可选地,步骤3)中针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程时,以表示粒子的前n维组成的列向量;表示粒子的n+1维到n+q维组成的列向量;表示粒子的n+q+1维到n+q+m维组成的列向量,建立任意k时刻的预测过程的状态方程如式(12)~(14)所示,建立任意k时刻的更新方程如式(15)~(21)所示;
式(12)~(14)中,表示k+1时刻的sigma点预测值,表示k时刻状态向量sigma点的估计值,u(k)表示输入向量,表示k时刻过程噪声向量的sigma点估计值;表示k+1时刻状态向量预测值,Wi m表示均值权值,L为sigma点采样个数;P(k+1,k)表示k+1时刻的协方差预测值;
式(15)~(21)中,zi(k+1,k)表示k+1时刻的量测向量估计值,表示k时刻状态向量sigma点的估计值,u(k)为输入向量,表示k+1时刻量测噪声向量的sigma点估计值;表示k+1时刻量测向量的预测值,Wi m表示均值权值,L为sigma点采样个数;Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差,αk+1为自适应调整因子,Wi c表示协方差权值;Pxz(k+1,k)表示k+1时刻状态向量和量测向量估计偏差的互协方差,表示k+1时刻的sigma点预测值,表示k+1时刻的状态向量预测值;K(k+1)表示k+1时刻的滤波增益;表示k+1时刻状态向量的最优估计值,z(k+1)表示k+1时刻的量测值;P(k+1,k+1)表示k+1时刻协方差最优估计值。
可选地,自适应调整因子αk+1的计算函数表达式如下式所示:
上式中,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为量测时间更新残差值,
可选地,步骤4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中的函数表达式如式(23)~(27)所示:
Dk+1=P(k+1,k)+Rk+1 (27)
式(23)~(27)中,表示k+1时刻状态向量的最优估计值,表示k+1时刻状态向量的预测值,K(k+1)表示k+1时刻的滤波增益值,φk+1表示修正函数,τk+1表示衡量k+1时刻量测误差大小的变量,ek+1表示k+1时刻的量测误差;z(k+1)表示k+1时刻的量测值;表示k+1时刻量测向量预测值;τk+1服从自由度为m的χ2分布,当出现野值时τk+1增大、φk+1(·)减小;Ck+1为选取的门限值或常数序列,λk+1为矩阵的最大特征值,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表示输入向量,建立频差的状态方程和量测方程:
上式中,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表示量测噪声向量的转置;确定频差状态向量的初始状态:
确定频差状态向量的扩维处理后的状态向量初始状态:
上式中,表示扩维后的状态向量初值,表示扩维后的状态向量均值,表示表示扩维后的状态向量,表示频差状态向量初值,表示扩维后的状态向量协方差初值,表示求取扩维后状态向量初值的协方差,P0表示频差向量的协方差向量初值,Q表示k时刻的过程噪声向量ω(k)的协方差矩阵,R表示k时刻的量测噪声向量γ(k)的协方差矩阵;
2.4)针对频差状态向量采用对称sigma采样,得到k时刻状态的sigma点集i=1,…,L,其中表示第i个维数为a的sigma点,L为sigma点采样个数,L=2a+1,此时状态维数a=n+q+m,其中n表示频差状态向量的维数,q表示过程噪声向量的维数,m表示量测噪声向量的维数;且sigma点集的函数表达式如下式所示:
上式中,表示状态向量的均值,Pxx表示状态向量的协方差,γ为系数,a为状态维数,λ=α2(a+κ)-a,κ为比例系数(比例系数κ通常可取值为0或者3-a,可用于调整sigma点和x均值之间的距离),γPxx的平方根矩阵的第i行或列记为
2.5)确定k时刻状态的sigma点集中各sigma点对应的权重:
上式中,Wi m表示均值权值,a为状态维数,λ=α2(a+κ)-a,κ为比例系数;
本实施例步骤3)中针对k时刻的sigma点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程时,以表示粒子的前n维组成的列向量;表示粒子的n+1维到n+q维组成的列向量;表示粒子的n+q+1维到n+q+m维组成的列向量,建立任意k时刻的预测过程的状态方程如式(12)~(14)所示,建立任意k时刻的更新方程如式(15)~(21)所示;
式(12)~(14)中,表示k+1时刻的sigma点预测值,表示k时刻状态向量sigma点的估计值,u(k)表示输入向量,表示k时刻过程噪声向量的sigma点估计值;表示k+1时刻状态向量预测值,Wi m表示均值权值,L为sigma点采样个数;P(k+1,k)表示k+1时刻的协方差预测值;
式(15)~(21)中,zi(k+1,k)表示k+1时刻的量测向量估计值,表示k时刻状态向量sigma点的估计值,u(k)为输入向量,表示k+1时刻量测噪声向量的sigma点估计值;表示k+1时刻量测向量的预测值,Wi m表示均值权值,L为sigma点采样个数;Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差,αk+1为自适应调整因子,Wi c表示协方差权值;Pxz(k+1,k)表示k+1时刻状态向量和量测向量估计偏差的互协方差,表示k+1时刻的sigma点预测值,表示k+1时刻的状态向量预测值;K(k+1)表示k+1时刻的滤波增益;表示k+1时刻状态向量的最优估计值,z(k+1)表示k+1时刻的量测值;P(k+1,k+1)表示k+1时刻协方差最优估计值。
本实施例中,自适应调整因子αk+1的计算函数表达式如下式所示:
上式中,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为量测时间更新残差值,
在实际应用中,由于GPS信号受噪声干扰还会产生较大的跳变野值,带有野值的频差信号将会使滤波结果发生偏移甚至发散。为了解决上述不敏卡尔曼算法的抗野值问题,有学者将新息序列加权的方法引入卡尔曼滤波来消除野值的影响。当观测数据不包含野值时,滤波算法能够正常运行,对观测噪声进行滤除;当观测数据包含野值时,能够克服或减少野值对于滤波结果的影响,使滤波结果更加准确。本实施例中,步骤4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中的函数表达式如式(23)~(27)所示:
Dk+1=P(k+1,k)+Rk+1 (27)
式(23)~(27)中,表示k+1时刻状态向量的最优估计值,表示k+1时刻状态向量的预测值,K(k+1)表示k+1时刻的滤波增益值,φk+1表示修正函数,τk+1表示衡量k+1时刻量测误差大小的变量,ek+1表示k+1时刻的量测误差;z(k+1)表示k+1时刻的量测值;表示k+1时刻量测向量预测值;τk+1服从自由度为m的χ2分布,当出现野值时τk+1增大、φk+1(·)减小;Ck+1为选取的门限值或常数序列,λk+1为矩阵的最大特征值,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点集确定自适应不敏卡尔曼滤波算法的状态方程和更新方程时,以表示粒子的前n维组成的列向量;表示粒子的n+1维到n+q维组成的列向量;表示粒子的n+q+1维到n+q+m维组成的列向量,建立任意k时刻的预测过程的状态方程如式(12)~(14)所示,建立任意k时刻的更新方程如式(15)~(21)所示;
式(12)~(14)中,表示k+1时刻的sigma点预测值,表示k时刻状态向量sigma点的估计值,u(k)表示输入向量,表示k时刻过程噪声向量的sigma点估计值;表示k+1时刻状态向量预测值,Wi m表示均值权值,L为sigma点采样个数;P(k+1,k)表示k+1时刻的协方差预测值;
式(15)~(21)中,zi(k+1,k)表示k+1时刻的量测向量估计值,表示k时刻状态向量sigma点的估计值,u(k)为输入向量,表示k+1时刻量测噪声向量的sigma点估计值;表示k+1时刻量测向量的预测值,Wi m表示均值权值,L为sigma点采样个数;Pzz(k+1,k)表示k+1时刻量测向量估计偏差的协方差,αk+1为自适应调整因子,Wi c表示协方差权值;Pxz(k+1,k)表示k+1时刻状态向量和量测向量估计偏差的互协方差,表示k+1时刻的sigma点预测值,表示k+1时刻的状态向量预测值;K(k+1)表示k+1时刻的滤波增益;表示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表示输入向量,建立频差的状态方程和量测方程:
上式中,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表示量测噪声向量的转置;确定频差状态向量的初始状态:
确定频差状态向量的扩维处理后的状态向量初始状态:
上式中,表示扩维后的状态向量初值,表示扩维后的状态向量均值,表示表示扩维后的状态向量,表示频差状态向量初值,表示扩维后的状态向量协方差初值,表示求取扩维后状态向量初值的协方差,P0表示频差向量的协方差向量初值,Q表示k时刻的过程噪声向量ω(k)的协方差矩阵,R表示k时刻的量测噪声向量γ(k)的协方差矩阵;
2.4)针对频差状态向量采用对称sigma采样,得到k时刻状态的sigma点集i=1,…,L,其中表示第i个维数为a的sigma点,L为sigma点采样个数,L=2a+1,此时状态维数a=n+q+m,其中n表示频差状态向量的维数,q表示过程噪声向量的维数,m表示量测噪声向量的维数;且sigma点集的函数表达式如下式所示:
2.5)确定k时刻状态的sigma点集中各sigma点对应的权重:
上式中,Wi m表示均值权值,a为状态维数,λ=α2(a+κ)-a,κ为比例系数;
5.根据权利要求3所述的基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法,其特征在于,步骤4)整合新息序列加权到自适应不敏卡尔曼滤波算法的状态方程和更新方程中的函数表达式如式(23)~(27)所示:
Dk+1=P(k+1,k)+Rk+1 (27)
式(23)~(27)中,表示k+1时刻状态向量的最优估计值,表示k+1时刻状态向量的预测值,K(k+1)表示k+1时刻的滤波增益值,φk+1表示修正函数,τk+1表示衡量k+1时刻量测误差大小的变量,ek+1表示k+1时刻的量测误差;z(k+1)表示k+1时刻的量测值;表示k+1时刻量测向量预测值;τk+1服从自由度为m的χ2分布,当出现野值时τk+1增大、φk+1(·)减小;Ck+1为选取的门限值或常数序列,λk+1为矩阵的最大特征值,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中任意一项所述基于新息加权自适应不敏卡尔曼滤波的晶振频率驯服方法的计算机程序。
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)
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)
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 | 国网湖南省电力有限公司 | 基于并网变换器的电力***振荡自适应抑制方法、***及介质 |
-
2020
- 2020-06-10 CN CN202010524159.9A patent/CN111650617B/zh active Active
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 |