CN106100609B - 单状态变量和两级Kalman滤波器时间尺度算法 - Google Patents

单状态变量和两级Kalman滤波器时间尺度算法 Download PDF

Info

Publication number
CN106100609B
CN106100609B CN201610386951.6A CN201610386951A CN106100609B CN 106100609 B CN106100609 B CN 106100609B CN 201610386951 A CN201610386951 A CN 201610386951A CN 106100609 B CN106100609 B CN 106100609B
Authority
CN
China
Prior art keywords
clock
time
time scale
group
kalman filter
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
CN201610386951.6A
Other languages
English (en)
Other versions
CN106100609A (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.)
Hunan Zhongdian Xinghe Electronics Co ltd
Original Assignee
National University of Defense Technology
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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN201610386951.6A priority Critical patent/CN106100609B/zh
Publication of CN106100609A publication Critical patent/CN106100609A/zh
Application granted granted Critical
Publication of CN106100609B publication Critical patent/CN106100609B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03HIMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00Adaptive networks
    • H03H21/0012Digital adaptive filters
    • H03H21/0025Particular filtering methods
    • H03H21/0029Particular filtering methods based on statistics
    • H03H21/003KALMAN filters
    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03HIMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00Adaptive networks
    • H03H21/0012Digital adaptive filters
    • H03H21/0043Adaptive algorithms

Landscapes

  • Physics & Mathematics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Compression, Expansion, Code Conversion, And Decoders (AREA)
  • Stabilization Of Oscillater, Synchronisation, Frequency Synthesizers (AREA)

Abstract

本发明提供了一种单状态变量和两级Kalman滤波器时间尺度算法,其中一种单状态变量Kalman滤波器时间尺度算法,包括以下步骤,首先建立原子钟模型;计算每台钟的过冲状态噪声方差;获取N‑1组观测钟差;用Kalman滤波器通过N‑1组观测钟差估计出N台钟的状态;计算时间尺度的权重和预测值;最后引入虚拟Kalman采样间隔,优化时间尺度TA。在该方法的基础上,本发明还提供了一种两级Kalman滤波器时间尺度算法,通过第一级Kalman滤波器和重构时差,获得N‑1组重构时差,然后,运行第二级的Kalman滤波器生成时间尺度TA。本发明可以通过选取虚拟Kalman采样间隔,使时间尺度的频率稳定度在任意某一个指定的平滑时间达到最优。

Description

单状态变量和两级Kalman滤波器时间尺度算法
技术领域
本发明涉及时间频率、信号处理领域,具体的说是设计了一种单状态变量Kalman滤波器时间尺度算法和一种两级Kalman滤波器时间尺度算法。
背景技术
时间尺度算法(time scale algorithms)在守时实验室和全球导航卫星系统(GNSS)中发挥着重要作用。时间尺度算法的目的是为了综合钟组中所有的原子钟,生成一个可靠性、频率稳定度更高的自由的纸面时间。单台钟和纸面时间都可以视为是一个时间尺度;没有被驾驭的时间尺度称为自由的时间尺度。所以说,使用时间尺度算法是为了生成一个时间尺度,该时间尺度一般通过它与某一台钟之间的时间差来表示;算法通常不考虑该时间尺度的频率准确度,将提升其频率准确度的任务留给驾驭算法来完成。
对于一个守时实验室,该自由的纸面时间通常记为TA(k)(其中k为实验室代号);可以用于监测物理钟的性能;也可以在使用钟差预测算法使TA(k)与UTC保持同步后,对物理钟进行驾驭,得到UTC的本地物理实现,记为UTC(k)。
对于GNSS,该自由的纸面时间即为自由的“系统时间”,记为自由的“GNSST”;例如对于北斗系统,记为自由的“BDT”。之所以打引号,是因为这时的“BDT”是自由的时间尺度,还没有与UTC(BSNC)(其中BSNC为北京卫星导航中心的代号)同步,不能称之为真正的系统时间。通过UTC(BSNC)对其进行驾驭,得到BDT,此时BDT与UTC(BSNC)在时间和频率上保持了同步。
时间尺度算法通常可以分为两大类:1)一类是经典的加权平均算法,例如:国际计量局(BIPM)采用的加权平均算法(称为:ALGOS算法)、美国国家标准与技术研究院(NIST)采用的加权平均算法(称为:AT1算法);对于任何加权平均算法,生成的时间尺度都可以用时间尺度基本方程来表示;2)另一类是采用Kalman滤波器实现的算法。由于AT1算法和Kalman滤波器算法所生成的时间尺度是实时的,比较适用于GNSS。
上世纪80年代NIST首次采用Kalman滤波器建立时间尺度。在一段时期内,该算法和AT1算法同时在NIST运行。该算法的不足在于:1)由于系统不完全可观测,造成估计误差矩阵发散,表明状态估计的误差不断变大;2)实验表明:对于铯钟和氢钟组成的钟组,时间尺度TA主要由长期频率稳定度最好的钟决定,从而没有利用到钟组中短期频率稳定度较好的钟的短期频率稳定度性能。
发明内容
为解决现有技术存在的问题,本发明提出一种单状态变量Kalman滤波器时间尺度算法;在单状态变量Kalman滤波器时间尺度算法的基础上,本发明还提供一种两级Kalman滤波器时间尺度算法。
本发明的技术方案是:
一种单状态变量Kalman滤波器时间尺度算法,包括以下步骤:
S1、建立原子钟模型
其中,X1(t)、X2(t)分别为两个钟差的状态变量;其中X1(t)代表瞬时时差,X2(t)代表不含频率白噪声(WFM)的瞬时频差;x0和y0分别代表瞬时时差和瞬时频差的初值;d代表频漂;W1(t)和W2(t)分别代表两个独立的维纳过程,并且有W(t)~N(0,t),即每个维纳过程服从均值为0,方差为时间t的正态分布;σ1和σ2分别是这两个维纳过程的扩散系数(diffusioncoefficients),用于表明噪声的强度。
扩散系数与Allan方差的关系表示为:
其中,代表平滑时间为τ时的Allan方差(是个约定俗成的表示Allan方差的方式,下标y的含义就是表明了Allan方差。),τ为平滑时间。
S2、计算每台钟的过冲状态噪声方差
钟组的状态方程为:
X(k+1)=φ·X(k)+J(k) (3)
其中,X(k)=[x1(k),x2(k),…,xN(k)]T是钟组在k时刻的时差;xi(k)(i=1,2,…,N)是钟组中第i台钟的时差;是状态转移矩阵;J(k)是钟组的过程噪声,表示为J(k)=[J1(k),J2(k),…,JN(k)]T,其中Ji(k)(i=1,2,…,N)是钟组中第i台钟的过程噪声。
由式(1),Ji(k)(i=1,2,…,N)可以表示为:
其中,W1(t),W2(t),σ1和σ2的含义和式(1)中相应的符号相同,下标i代表第i台钟。
第i台钟的过程状态噪声方差即为式(4)的方差,表示为:
其中,T为采样间隔;
由式(5)和式(2),得到:
其中,表示第i台钟的Allan方差。
于是,钟组的过程噪声方差表示为:
S3、获取输入量即N-1组观测钟差
钟组的观测方程为
Z(k)=H·X(k)+V(k) (8)
其中,H表示为:
V是测量噪声,表示为V(k)=[v1(k),v2(k),…,vN-1(k)]T,其中vi(k)(i=1,2,…,N-1)是V每一组钟差的观测噪声,它们之间相互独立;vi(k)的噪声协方差为Ri,所以钟组的测量噪声V的协方差表示为:
由式(8)和式(9),观测钟差可以表示为:
Z(k)=[x1(k)-x2(k)+v1(k) … x1(k)-xN(k)+vN-1(k)]T (11)
其中的每一个量代表在k时刻第1台钟与第i台钟(i≠1)之间的钟差。
S4、使用Kalman滤波器通过N-1组观测钟差估计出N台钟的状态
Pk,k-1=φ·Pk-1,k-1·φT+Q (13)
Kk=Pk,k-1·HT(H·Pk,k-1·HT+R)-1 (14)
Pk,k=(I-Kk·H)·Pk,k-1 (16)
其中:代表了钟组在k-1时刻的时差的估计值、代表了钟组在k-1时刻对k时刻的时差的预测值、Pk-1,k-1是钟组时差的估计误差矩阵、Pk,k-1是钟组时差的预测误差矩阵、Kk是Kalman滤波器增益。
通过式(12)至至式(16)五个方程得到是第i台钟相对于时间尺度TA的时差的估计值。
S5、计算时间尺度的权重和预测值
将式(12)代入式(15),得到:
式(17)可以写成如下形式:
其中:K11(k)表示Kk矩阵的第1行第1列的元素;K12(k)表示Kk矩阵的第1行第2列的元素;依次类推,KN(N-1)(k)表示Kk矩阵的第N行第N-1列的元素;
将每台钟相对于时间尺度TA的时间差减去其时间差的估计值得到即为校准钟,由式(18)得到第1台校准钟表示为:
当观测噪声为零时,式(19)满足时间尺度基本方程,即:
其中,时间尺度的权重表示为:
时间尺度的预测值表示为:
当测量噪声为零时,由式(18),任何一台“校准钟”都满足时间尺度基本方程,即:
由式(18),当i≠1时,时间尺度的权重表示为:
其中:Kij(k)表示Kk矩阵的第i行第j列的元素;
而时间尺度的预测值的表达式与式(22)相同。
S6、引入虚拟Kalman采样间隔,优化时间尺度TA,方法为:
引入虚拟Kalman采样间隔Tv用于计算得到每台钟的虚拟状态噪声协方差,表示为:
按照步骤S2至S5的方法,TA将表示为:
本发明还提供一种两级Kalman滤波器时间尺度算法,包括以下步骤:
S1、通过N-1个Kalman滤波器分别对N-1组钟差进行状态估计,得到每组钟差的两个状态的估计量,记为
N-1个Kalman滤波器组成第一级Kalman滤波器,第一级Kalman滤波器其输入量分别为N-1组观测钟差,把其中每组钟差分别记为Z(1,n)(k),其中上标(1,n)表示第1台钟与第n台钟之间的钟差,第i台钟的时差记为xi(k),通过这N-1个Kalman滤波器,估计得到N-1组钟差的两个状态的估计量;
S2、对时差进行重构
对于每一组钟差估计值,令
序列的初始值取为零,上式可以写为:
于是,式(28)就是通过估计量进行重构,得到的第1台钟和第n台钟之间的重构时差
把第i台钟的重构钟的时差记为于是有其中为重构时差;
S3、第二级Kalman滤波器采用如权利要求1所述的单状态变量Kalman滤波器时间尺度算法生成时间尺度TA;
输入量为N-1组重构时差输出为N台重构钟的时差估计值,记为代表了第i台重构钟的时差与时间尺度TA之间的偏差,于是得到:
本发明提出一种单状态变量Kalman滤波器时间尺度算法。其中,每台钟只有1个状态变量,即每台钟的时差。在获得了权重的解析表达式的基础上,引入了虚拟Kalman采样间隔,用于计算过程噪声协方差;并证明在这种情况下,权重将反比于平滑时间为虚拟Kalman采样间隔的Allan方差;这样,可以通过选取虚拟Kalman采样间隔,使时间尺度的频率稳定度在任意某一个指定的平滑时间达到最优。
在此基础上,本发明还提出了一种两级Kalman滤波器时间尺度算法。第一级的Kalman滤波器用于估计钟差的两个状态。其中,第一个状态中包含了频率白噪声(WFM)和频率随机游走噪声(RWFM)两种噪声;而第二个状态中只包含了RWFM。第二级的Kalman滤波器实际上是一个单状态变量Kalman滤波器时间尺度算法,用于生成时间尺度;其输入为第一级Kalman滤波器的第二个状态估计值的“重构时差”;这样,算法相当于对只含有RWFM的时间序列进行加权,生成的时间尺度中也只含有RWFM这一种噪声,中短期频率稳定度更高。
附图说明
图1图解说明采用第一种算法时,在第一种情况下,每台钟、TA的时差(a)和Allan偏差(b);
图2图解说明采用第一种算法时,在第二种情况下,每台钟、TA的时差(a)和Allan偏差(b);
图3图解说明两级Kalman滤波器算法原理;
图4图解说明采用第二种算法时,每台“重构钟”、TA的时差(a)和Allan偏差(b);
具体实施方式
下面结合附图对本发明的具体实施方式作出说明。
首先,时间尺度基本方程表示为:
其中,TA(t)是生成的纸面上的时间尺度;N为钟组中原子钟的数量;ωi为第i台原子钟的权重;xi(t)代表第i台原子钟的时差;xi’(t)代表第i台原子钟的时差的预测值。
时间尺度基本方程可以这样理解:对于加权平均时间尺度算法,权重用于优化时间尺度的频率稳定度,而预测值则消除了每台钟确定性趋势项的影响;所以TA的物理意义是对去除确定性趋势项后只含有噪声分量的原子钟时差的加权平均值。
事实上,不论是ALGOS算法还是AT1算法,其生成的时间尺度都可以进行化简,最终使用时间尺度基本方程来表示。
本发明提供了一种单状态变量Kalman滤波器时间尺度算法,包括以下步骤:
S1、建立原子钟模型
其中,X1(t)、X2(t)分别为两个钟差的状态变量;其中X1(t)代表瞬时时差,X2(t)代表不含频率白噪声(WFM)的瞬时频差;x0和y0分别代表瞬时时差和瞬时频差的初值;d代表频漂;W1(t)和W2(t)分别代表两个独立的维纳过程,并且有W(t)~N(0,t),即每个维纳过程服从均值为0,方差为时间t的正态分布;σ1和σ2分别是这两个维纳过程的扩散系数(diffusioncoefficients),用于表明噪声的强度。
式(1)表明:W1(t)和W2(t)对时间t的积分在状态变量X1上分别表现为WFM和RWFM;即式(1)的第一个方程中,最后第二项代表WFM,最后一项代表RWFM。
扩散系数与Allan方差的关系表示为:
其中,代表平滑时间为τ时的Allan方差,τ为平滑时间.式(2)中等号后面第一项的斜率为-1,第二项斜率为1,这两项分别对应WFM和RWFM,这说明在对数Allan方差图中,WFM的斜率为-1,RWFM的斜率为1。
S2、计算每台钟的过冲状态噪声方差
钟组的状态方程为:
X(k+1)=φ·X(k)+J(k) (3)
其中,X(k)=[x1(k),x2(k),…,xN(k)]T是钟组在k时刻的时差;xi(k)(i=1,2,…,N)是钟组中第i台钟的时差;是状态转移矩阵;J(k)是钟组的过程噪声,表示为J(k)=[J1(k),J2(k),…,JN(k)]T,其中Ji(k)(i=1,2,…,N)是钟组中第i台钟的过程噪声;由式(1),Ji(k)(i=1,2,…,N)可以表示为:
第i台钟的过程状态噪声方差即为式(4)的方差,表示为:
由式(5)和式(2),得到:
其中,表示第i台钟的Allan方差;
于是,钟组的过程噪声方差表示为:
S3、获取输入量即N-1组观测钟差
钟组的观测方程为
Z(k)=H·X(k)+V(k) (8)
其中,H表示为:
V是测量噪声,表示为V(k)=[v1(k),v2(k),…,vN-1(k)]T,其中vi(k)(i=1,2,…,N-1)是V中每一组钟差的观测噪声,它们之间相互独立;vi(k)的噪声协方差为Ri,所以钟组的测量噪声V的协方差表示为:
由式(8)和式(9),观测钟差可以表示为:
Z(k)=[x1(k)-x2(k)+v1(k) … x1(k)-xN(k)+vN-1(k)]T (11)
其中的每一个量代表在k时刻第1台钟与第i台钟(i≠1)之间的钟差;
S4、使用Kalman滤波器通过N-1组观测钟差估计出N台钟的状态
Kalman滤波器应用于时间尺度建立时,Kalman滤波器的作用是通过N-1组观测钟差估计出N台钟的状态,记为和所有Kalman滤波器一样,此时Kalman滤波器的步骤可以由下面5个方程来描述。
Pk,k-1=φ·Pk-1,k-1·φT+Q (13)
Kk=Pk,k-1·HT(H·Pk,k-1·HT+R)-1 (14)
Pk,k=(I-Kk·H)·Pk,k-1 (16)
其中:代表了钟组在k-1时刻的时差的估计值、代表了钟组在k-1时刻对k时刻的时差的预测值、Pk-1,k-1是钟组时差的估计误差矩阵、Pk,k-1是钟组时差的预测误差矩阵、Kk是Kalman滤波器增益。
通过式(12)至至式(16)五个方程得到是第i台钟相对于时间尺度TA的时差的估计值。也就是说:Kalman估计得到的量是每台钟相对于时间尺度TA的时间差。
S5、计算时间尺度的权重和预测值
将式(12)代入式(15),得到:
式(17)可以写成如下形式:
将每台钟相对于时间尺度TA的时间差减去其时间差的估计值得到即为校准钟,由式(18)得到第1台校准钟表示为:
当观测噪声为零时,式(19)满足时间尺度基本方程,即:
其中,时间尺度的权重表示为:
时间尺度的预测值表示为:
当测量噪声为零时,由式(18),任何一台“校准钟”都满足时间尺度基本方程,即:
由式(18),当i≠1时,时间尺度的权重表示为:
而时间尺度的预测值的表达式与式(22)相同。
S6、引入虚拟Kalman采样间隔,优化时间尺度TA,方法为:
引入虚拟Kalman采样间隔Tv用于计算得到每台钟的虚拟状态噪声协方差,表示为:
按照步骤S2至S5的方法,TA将表示为:
式(26)表明:在采用虚拟Kalman采样间隔Tv时,算法可以使时间尺度TA在平滑时间为Tv的频率稳定度达到了最优。因此,在实际运用中,可以直接选取一个Tv的值,这样生成的TA就在指定的平滑时间Tv的频率稳定度达到了最优。
采用上述方法进行仿真验证:
首先仿真生成2台氢钟和1台铯钟组成钟组(每台钟共5000个数据点,相邻数据点之间的采样间隔为3600s)。每台氢钟如式(1)所示的平方扩散系数设为σ1 2=5.0×10-25s和σ2 2=1.3×10-35s-1;铯钟的平方扩散系数设为σ1 2=4.8×10-23s和σ2 2=1.9×10-36s-1。这些扩散系数符合典型氢钟和铯钟的指标。设置观测噪声为零。为了简化分析,设每台钟如式(1)所示的确定性分量为x0=y0=d=0。
然后,按照上述方法生成时间尺度TA。分为两种情况。第一种情况,设置虚拟Kalman采样间隔Tv=432000s(5天);第二种情况,设置Tv=2592000s(30天)。当观测噪声为零时,所有“校准钟”都是相等的,所以不妨把第1台“校准钟”作为TA。
表1列出了在这两种情况下,由式(26)计算得到的每台钟权重的理论值。图1描绘了在第一种情况下,这三台钟、TA的时差(a)和Allan偏差(Allan方差的平方根)(b)。图2描绘了在第二种情况下,它们的时差(a)和Allan偏差(b)。
表1两种情况下每台钟的权重
由表1、图1和图2,可以看出:在Tv较小时,氢钟的权重较大,于是TA主要由两台氢钟决定;而在Tv较大时,铯钟的权重较大,于是TA主要由铯钟决定。
本发明提供的一种单状态变量Kalman滤波器时间尺度算法实质上是一种加权平均算法;当测量噪声为零时,权重反比于每台钟在平滑时间为Tv时的Allan偏差。所以,可以根据实际需求,灵活地设置Tv的值,使TA在该平滑间隔的频率稳定度达到最优。
本发明还提供了一种两级Kalman滤波器时间尺度算法,其原理如图3所示,包括以下步骤:
S1、通过N-1个Kalman滤波器分别对N-1组钟差进行状态估计,得到每组钟差的两个状态的估计量,记为
N-1个Kalman滤波器组成第一级Kalman滤波器,第一级Kalman滤波器其输入量分别为N-1组钟差观测值(如式11所示),把其中每组钟差分别记为Z(1,n)(k),其中上标(1,n)表示第1台钟与第n台钟之间的钟差,第i台钟的时差记为xi(k),;如果没有观测噪声,显然有Z(1,n)(k)=x1(k)-xn(k)。
通过这N-1个Kalman滤波器分别对这N-1组钟差进行状态估计,得到每组钟差的如式(1)所示的两个状态的估计量,记为由式(1)可知,时间序列中的噪声为WFM和RWFM;而时间序列中将只含有RWFM。
S2、对时差进行重构
对于每一组钟差估计值,令
序列的初始值取为零,上式可以写为:
于是,式(28)就是通过估计量进行重构,得到的第1台钟和第n台钟之间的重构时差显然,“重构时差”中也只含有RWFM。
把第i台钟的重构钟的时差记为于是有其中为重构时差。
通过第一级Kalman滤波器和重构时差这两个环节,相当于滤除了输入时间序列(即观测钟差)中的WFM,得到了只含有RWFM的“重构时差”。能达到这样的效果,是利用了以下特性:Kalman滤波器用于状态估计时,按照原子钟模型即式(1),估计得到的第1个状态变量中含有WFM和RWFM,而第2个状态变量中只含有RWFM。
S3、第二级Kalman滤波器采用如权利要求1所述的单状态变量Kalman滤波器时间尺度算法生成TA;
输入量为N-1组重构时差输出为N台重构钟的时差估计值,记为代表了第i台重构钟的时差与TA之间的偏差,于是得到:
采用上述方法进行仿真验证:
实际应用中的步骤应该为:1)首先使用第一级的Kalman滤波器,并对时差进行重构,获得N-1组“重构时差”,即:H1-H2和H1-Cs的“重构时差”;2)然后,运行如第二级的Kalman滤波器,得到的估计量为每台“重构钟”与TA之间的时间差。
在本仿真实验中,假设单台钟的时差是已知的。
在本次仿真中,首先使用第一级的3个Kalman滤波器,分别对3台钟进行状态估计,对时差进行重构,获得3台“重构钟”的时差,显然也就得到了2组“重构时差”。然后,运行第二级的Kalman滤波器,其输入为两组“重构时差”,输出为估计得到的每台“重构钟”与TA之间的时间差,于是也就生成了TA。其中,设置虚拟Kalman采样间隔Tv=432000s(即5天)。
表2每台“重构钟”的权重
Table 2the weights of each“reconstructed clock”
表2列出了每台“重构钟”的权重。从表2中可以看出:“重构铯钟”获得的权重更大,权重反比于σ2 2的值。图4描述了每台“重构钟”、TA的时差(a)和Allan偏差(b)。从图4中可以看出:TA中也只含有RWFM;它的中短期频率稳定度相比第3节生成的TA得到了明显的提高。
综上所述,虽然本发明已以较佳实施例揭露如上,然其并非用以限定本发明,任何本领域普通技术人员,在不脱离本发明的精神和范围内,当可作各种更动与润饰,因此本发明的保护范围当视权利要求书界定的范围为准。

Claims (3)

1.一种单状态变量Kalman滤波器时间尺度算法,其特征在于,包括以下步骤:
S1建立原子钟模型
其中,X1(t)、X2(t)分别为观测钟差的两个状态变量;其中X1(t)代表瞬时时差,X2(t)代表不含频率白噪声的瞬时频差;x0和y0分别代表瞬时时差和瞬时频差的初值;d代表频漂;W1(t)和W2(t)分别代表两个独立的维纳过程,并且有W(t)~N(0,t),即每个维纳过程服从均值为0,方差为时间t的正态分布;σ1和σ2分别是这两个维纳过程的扩散系数,用于表明噪声的强度;
扩散系数与Allan方差的关系表示为:
其中,代表平滑时间为τ时的Allan方差,τ为平滑时间;
S2计算每台钟的过冲状态噪声方差
钟组的状态方程为:
X(k+1)=φ·X(k)+J(k) (3)
其中,X(k)=[x1(k),x2(k),…,xN(k)]T是钟组在k时刻的时差;xi(k)(i=1,2,…,N)是钟组中第i台钟的时差;是状态转移矩阵;J(k)是钟组的过程噪声,表示为J(k)=[J1(k),J2(k),…,JN(k)]T,其中Ji(k)(i=1,2,…,N)是钟组中第i台钟的过程噪声;由式(1),Ji(k)(i=1,2,…,N)可以表示为:
第i台钟的过程状态噪声方差即为式(4)的方差,表示为:
其中,T为采样间隔;
由式(5)和式(2),得到:
其中,表示第i台钟的Allan方差;
于是,钟组的过程噪声方差表示为:
S3获取输入量即N-1组观测钟差
钟组的观测方程为
Z(k)=H·X(k)+V(k) (8)
其中,H表示为:
V是测量噪声,表示为V(k)=[v1(k),v2(k),…,vN-1(k)]T,其中vi(k)(i=1,2,…,N-1)是V中每一组钟差的观测噪声,它们之间相互独立;vi(k)的噪声协方差为Ri,所以钟组的测量噪声V的协方差表示为:
由式(8)和式(9),观测钟差可以表示为:
Z(k)=[x1(k)-x2(k)+v1(k)…x1(k)-xN(k)+vN-1(k)]T (11)
其中的每一个量代表在k时刻第1台钟与第i台钟(i≠1)之间的钟差;
S4使用Kalman滤波器通过N-1组观测钟差估计出N台钟的状态
Kk=Pk,k-1·HT(H·Pk,k-1·HT+R)-1 (14)
Pk,k=(I-Kk·H)·Pk,k-1 (16)
其中:代表了钟组在k-1时刻的时差的估计值、代表了钟组在k-1时刻对k时刻的时差的预测值、Pk-1,k-1是钟组时差的估计误差矩阵、Pk,k-1是钟组时差的预测误差矩阵、Kk是Kalman滤波器增益;
通过式(12)至至式(16)五个方程得到 是第i台钟相对于时间尺度TA的时差的估计值;
S5计算时间尺度的权重和预测值
将式(12)代入式(15),得到:
式(17)可以写成如下形式:
其中:K11(k)表示Kk矩阵的第1行第1列的元素;K12(k)表示Kk矩阵的第1行第2列的元素;依次类推,KN(N-1)(k)表示Kk矩阵的第N行第N-1列的元素;
将每台钟相对于时间尺度TA的时间差减去其时差的估计值得到即为校准钟,由式(18)得到第1台校准钟表示为:
当观测噪声为零时,式(19)满足时间尺度基本方程,即:
其中,时间尺度的权重表示为:
时间尺度的预测值表示为:
当测量噪声为零时,由式(18),任何一台“校准钟”都满足时间尺度基本方程,即:
由式(18),当i≠1时,时间尺度的权重表示为:
其中:Kij(k)表示Kk矩阵的第i行第j列的元素;
而时间尺度的预测值的表达式与式(22)相同。
2.根据权利要求1所述的单状态变量Kalman滤波器时间尺度算法,其特征在于,还包括步骤S6,引入虚拟Kalman采样间隔,优化时间尺度TA,方法为:
引入虚拟Kalman采样间隔Tv用于计算得到每台钟的虚拟状态噪声协方差,表示为:
按照步骤S2至S5的方法,时间尺度TA将表示为:
3.一种两级Kalman滤波器时间尺度算法,其特征在于,包括以下步骤:
S1、通过N-1个Kalman滤波器分别对N-1组钟差进行状态估计,得到每组钟差的两个状态的估计量,记为
N-1个Kalman滤波器组成第一级Kalman滤波器,第一级Kalman滤波器其输入量分别为N-1组观测钟差,把其中每组钟差分别记为Z(1,n)(k),其中上标(1,n)表示第1台钟与第n台钟之间的钟差,第i台钟的时差记为xi(k),通过这N-1个Kalman滤波器,估计得到N-1组钟差的两个状态的估计量;
S2、对时差进行重构
对于每一组钟差估计值,令
序列的初始值取为零,上式可以写为:
于是,式(28)就是通过估计量进行重构,得到的第1台钟和第n台钟之间的重构时差
把第i台钟的重构钟的时差记为于是有其中为重构时差;
S3、第二级Kalman滤波器采用如权利要求1所述的单状态变量Kalman滤波器时间尺度算法生成时间尺度TA;
输入量为N-1组重构时差输出为N台重构钟的时差估计值,记为代表了第i台重构钟的时差与TA之间的偏差,于是得到:
CN201610386951.6A 2016-06-02 2016-06-02 单状态变量和两级Kalman滤波器时间尺度算法 Active CN106100609B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610386951.6A CN106100609B (zh) 2016-06-02 2016-06-02 单状态变量和两级Kalman滤波器时间尺度算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610386951.6A CN106100609B (zh) 2016-06-02 2016-06-02 单状态变量和两级Kalman滤波器时间尺度算法

Publications (2)

Publication Number Publication Date
CN106100609A CN106100609A (zh) 2016-11-09
CN106100609B true CN106100609B (zh) 2018-08-31

Family

ID=57446848

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610386951.6A Active CN106100609B (zh) 2016-06-02 2016-06-02 单状态变量和两级Kalman滤波器时间尺度算法

Country Status (1)

Country Link
CN (1) CN106100609B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108229747B (zh) * 2018-01-12 2020-12-15 中国计量科学研究院 校准控制方法、装置及时间信号产生系统
CN110865530A (zh) * 2019-11-27 2020-03-06 国网思极神往位置服务(北京)有限公司 一种原子时计算方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103148856A (zh) * 2013-03-04 2013-06-12 北京航空航天大学 一种基于自适应尺度变化的借力飞行探测器自主天文导航方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103148856A (zh) * 2013-03-04 2013-06-12 北京航空航天大学 一种基于自适应尺度变化的借力飞行探测器自主天文导航方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
基于小波降噪的 Kalman 滤波预报卫星钟差;王继刚;《测绘科学》;20120930;全文 *
基于消参数双向Kalman 滤波的GEO 卫星定轨;贺凯飞等;《测绘科学》;20110331;全文 *
导航卫星原子钟Kalman 滤波中噪声方差-协方差的确定;郭海荣;《测绘学报》;20100430;全文 *

Also Published As

Publication number Publication date
CN106100609A (zh) 2016-11-09

Similar Documents

Publication Publication Date Title
Blasques et al. In-sample confidence bands and out-of-sample forecast bands for time-varying parameters in observation-driven models
Diao et al. Analysis and compensation of MEMS gyroscope drift
Zhou et al. Critical issues on Kalman filter with colored and correlated system noises
FR3028031A1 (fr) Procede d'estimation d'un etat de navigation contraint en observabilite
Nakamori et al. Recursive estimators of signals from measurements with stochastic delays using covariance information
CN106100609B (zh) 单状态变量和两级Kalman滤波器时间尺度算法
Wu et al. A Kalman filter approach based on random drift data of fiber optic gyro
Guerrier et al. Wavelet-based improvements for inertial sensor error modeling
Rasoulzadeh et al. Implementation of A low-cost multi-IMU hardware by using a homogenous multi-sensor fusion
Padilla Recursive identification of continuous-time systems with time-varying parameters
DeMars et al. Modular framework for implementation and analysis of recursive filters with considered and neglected parameters
Lee et al. State and parameter estimation using measurements with unknown time delay
Cox et al. Estimation of LPV-SS models with static dependency using correlation analysis
Yakimov et al. Digital estimation of correlation function moments using analog-stochastic sign quantization of a random process
CN106096787B (zh) 一种自适应匹配的地球自转参数预报方法
CN111950123B (zh) 一种陀螺仪误差系数曲线拟合预测方法及系统
Mohiuddin et al. Continuous-time Kalman filtering with implicit discrete measurement times
Cordeiro et al. Boot. EXPOS in NNGC competition
Greenhall Reduced Kalman filters for clock ensembles
Sobolewski Predicting the Lithuanian Timescale UTC (LT) by means of GMDH neural network
CN112577478A (zh) 微机电系统陀螺随机噪声的处理方法及处理装置
Lindsey On the Kalman filter and its variations
CN112698368B (zh) 一种导航接收机的导航信号解析方法及计算机存可读介质
Al-Smadi Automatic identification of ARMA systems
CN117406253A (zh) 基于LSTM和二次多项式组合模型的Multi-GNSS卫星钟差建模预报方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CP01 Change in the name or title of a patent holder
CP01 Change in the name or title of a patent holder

Address after: 410073 Hunan province Changsha Kaifu District, Deya Road No. 109

Patentee after: National University of Defense Technology

Address before: 410073 Hunan province Changsha Kaifu District, Deya Road No. 109

Patentee before: NATIONAL University OF DEFENSE TECHNOLOGY

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220628

Address after: 410000 block a, building 1, Changsha National Security Industrial Park, No. 699 Qingshan Road, Yuelu District, Changsha City, Hunan Province

Patentee after: Hunan Institute of advanced technology

Address before: 410073 Hunan province Changsha Kaifu District, Deya Road No. 109

Patentee before: National University of Defense Technology

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20221213

Address after: Building 4, Hunan Military civilian Integration Science and Technology Innovation Industrial Park, No. 699, Qingshan Road, Changsha Hi tech Development Zone, 410000, Hunan

Patentee after: Hunan Zhongdian Xinghe Electronics Co.,Ltd.

Address before: 410000 block a, building 1, Changsha National Security Industrial Park, No. 699 Qingshan Road, Yuelu District, Changsha City, Hunan Province

Patentee before: Hunan Institute of advanced technology