CN106100609B - 单状态变量和两级Kalman滤波器时间尺度算法 - Google Patents
单状态变量和两级Kalman滤波器时间尺度算法 Download PDFInfo
- 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
- formula
- 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
Links
- 238000000034 method Methods 0.000 claims abstract description 24
- 238000005070 sampling Methods 0.000 claims abstract description 17
- 239000011159 matrix material Substances 0.000 claims description 18
- 238000009499 grossing Methods 0.000 claims description 11
- 238000005259 measurement Methods 0.000 claims description 10
- 230000005653 Brownian motion process Effects 0.000 claims description 9
- 238000009792 diffusion process Methods 0.000 claims description 9
- 230000007704 transition Effects 0.000 claims description 3
- 238000005457 optimization Methods 0.000 abstract 1
- 229910052792 caesium Inorganic materials 0.000 description 7
- TVFDJXOCXUVLDH-UHFFFAOYSA-N caesium atom Chemical compound [Cs] TVFDJXOCXUVLDH-UHFFFAOYSA-N 0.000 description 7
- UFHFLCQGNIYNRP-UHFFFAOYSA-N Hydrogen Chemical compound [H][H] UFHFLCQGNIYNRP-UHFFFAOYSA-N 0.000 description 6
- 229910052739 hydrogen Inorganic materials 0.000 description 6
- 239000001257 hydrogen Substances 0.000 description 6
- 238000004088 simulation Methods 0.000 description 4
- 230000001360 synchronised effect Effects 0.000 description 3
- 238000012795 verification Methods 0.000 description 2
- MRMBZHPJVKCOMA-YJFSRANCSA-N biapenem Chemical compound C1N2C=NC=[N+]2CC1SC([C@@H]1C)=C(C([O-])=O)N2[C@H]1[C@@H]([C@H](O)C)C2=O MRMBZHPJVKCOMA-YJFSRANCSA-N 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000005295 random walk Methods 0.000 description 1
Classifications
-
- 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
Landscapes
- Physics & Mathematics (AREA)
- Probability & Statistics with Applications (AREA)
- Stabilization Of Oscillater, Synchronisation, Frequency Synthesizers (AREA)
- Compression, Expansion, Code Conversion, And Decoders (AREA)
Abstract
本发明提供了一种单状态变量和两级Kalman滤波器时间尺度算法，其中一种单状态变量Kalman滤波器时间尺度算法，包括以下步骤，首先建立原子钟模型；计算每台钟的过冲状态噪声方差；获取N‑1组观测钟差；用Kalman滤波器通过N‑1组观测钟差估计出N台钟的状态；计算时间尺度的权重和预测值；最后引入虚拟Kalman采样间隔，优化时间尺度TA。在该方法的基础上，本发明还提供了一种两级Kalman滤波器时间尺度算法，通过第一级Kalman滤波器和重构时差，获得N‑1组重构时差，然后，运行第二级的Kalman滤波器生成时间尺度TA。本发明可以通过选取虚拟Kalman采样间隔，使时间尺度的频率稳定度在任意某一个指定的平滑时间达到最优。
Description
技术领域
本发明涉及时间频率、信号处理领域，具体的说是设计了一种单状态变量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之间的偏差，于是得到：
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)
| 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)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103148856A (zh) * | 2013-03-04 | 2013-06-12 | 北京航空航天大学 | 一种基于自适应尺度变化的借力飞行探测器自主天文导航方法 |
-
2016
- 2016-06-02 CN CN201610386951.6A patent/CN106100609B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103148856A (zh) * | 2013-03-04 | 2013-06-12 | 北京航空航天大学 | 一种基于自适应尺度变化的借力飞行探测器自主天文导航方法 |
Non-Patent Citations (3)
| 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 |
|---|---|---|
| Vaccaro et al. | Statistical modeling of rate gyros | |
| Trumbore et al. | Radiocarbon nomenclature, theory, models, and interpretation: Measuring age, determining cycling rates, and tracing source pools | |
| Friedrichs | A data assimilative marine ecosystem model of the central equatorial Pacific: Numerical twin experiments | |
| Karspeck et al. | A global coupled ensemble data assimilation system using the Community Earth System Model and the Data Assimilation Research Testbed | |
| Bustos et al. | Health Monitoring of Lithium-Ion Batteries Using Dual Filters | |
| Possolo et al. | Evaluating the uncertainty of input quantities in measurement models | |
| Rahimifard et al. | Adaptive smooth variable structure filter strategy for state estimation of electric vehicle batteries | |
| CN106100609B (zh) | 单状态变量和两级Kalman滤波器时间尺度算法 | |
| Abdelfatah et al. | FPGA-based real-time embedded system for RISS/GPS integrated navigation | |
| Mu et al. | Design and implementation of vector tracking loop for high-dynamic gnss receiver | |
| Ramillien et al. | Recovery of rapid water mass changes (RWMC) by Kalman filtering of GRACE observations | |
| Quaglino et al. | High-dimensional and higher-order multifidelity Monte Carlo estimators | |
| Mamoutos et al. | A high-resolution numerical model of the North Aegean Sea aimed at climatological studies | |
| Afthanorhan et al. | Modelling the multimediator on motivation among youth in higher education institution towards volunteerism program | |
| Wang | On the latent state estimation of nonlinear population dynamics using Bayesian and non-Bayesian state-space models | |
| Cortés et al. | Evaluation of low-complexity adaptive full direct-state Kalman filter for robust GNSS tracking | |
| Song et al. | Improved multiple-model adaptive estimation method for integrated navigation with time-varying noise | |
| Li et al. | Rainfall erosivity characteristics during 1961–2100 in the Loess Plateau, China | |
| CN111623764B (zh) | 微纳卫星姿态估计方法 | |
| Desouky et al. | A Recursive Approach for Magnetic Field Estimation in Spacecraft Magnetic Attitude Control | |
| Artstein et al. | Analysis and computation of a discrete KdV-Burgers type equation with fast dispersion and slow diffusion | |
| Neves et al. | The robustness of a simple dynamic model of island biodiversity to geological and eustatic change | |
| Alken et al. | Improved horizontal wind model HWM07 enables estimation of equatorial ionospheric electric fields from satellite magnetic measurements | |
| Nickles et al. | Leveraging river network topology and regionalization to expand SWOT-derived river discharge time series in the Mississippi River basin | |
| CN111950123B (zh) | 一种陀螺仪误差系数曲线拟合预测方法及系统 |
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 |
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 |
|
| CP01 | Change in the name or title of a patent holder | ||
| 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 |
|
| TR01 | Transfer of patent right |