等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法
等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法
【技术领域】
[0001] 本发明涉及时间频率、信号处理领域，具体的说是设计了一种等价于二状态变量 Kalman滤波器的数字锁相环原子钟驾驭方法。
【背景技术】
[0002] 原子钟驾驭技术在守时实验室和卫星导航系统中发挥着重要作用。对原子钟进行 驾驭的主要目的有两个：一是使被驾驭原子钟与参考原子钟时间同步，尽可能缩小它们之 间的时间偏差；二是提升被驾驭原子钟的长期稳定度。
[0003] 典型的原子钟驾驭方法包括两大类：开环驾驭方法和闭环驾驭方法。开环驾驭方 法的核心是设计合理的预测算法。闭环驾驭方法主要包括：线性二次型（LG?控制方法、开 关（Bang-Bang)控制方法和数字锁相环（DPLL)方法。但是，这些闭环驾驭方法的参数是不 容易选取的，一般都要针对每一个具体应用，通过大量仿真实验对大量参数进行比较后，选 出一组最优参数。另外，对于DPLL方法，如果参数选取不当，不但无法保证驾驭性能，还会 引起控制系统不稳定。
【发明内容】
[0004] 为解决现有技术中存在的问题，本发明的目的在于提供一种等价于二状态变量 Kalman滤波器的数字锁相环原子钟驾驭方法。
[0005] 本发明的原理是：通过一个等价于二状态变量Kalman滤波器的DPLL，用铯钟对氢 钟进行驾驭。本发明推导了该DPLL的闭环系统传递函数，证明该DPLL为二阶2类DPLL。 本发明推导了由过程噪声方差和观测噪声方差比值近似决定的Kalman增益。于是，由DPLL 的闭环传递函数可知，该比值，以及环路增益，采样间隔这3个参数决定了该DPLL的性能。 进一步，本发明给出了一个简单有效的选取参数的方法。通过等效变换，得到了等效于该 DPLL的模拟锁相环（APLL)。在频域中，通过调整参数，使该等价APLL的系统传递函数和误 差传递函数的交点频率等于氢钟和铯钟的SSB相位噪声曲线的交点频率，最终获取最优参 数。简而言之，使用参考原子钟（或时间尺度）来驾驭被驾驭原子钟（或时间尺度），使被 驾驭原子钟（或时间尺度）与参考原子钟（或时间尺度）保持同步。
[0006] 本发明采用的技术方案为：
[0007] -种等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法，采用以下 步骤：
[0008] S. 1推导等价于二状态变量Kalman滤波器的DPLL的传递函数；
[0009] S. 1. 1使用DPLL，用铯钟来驾驭氢钟。图1描述了该DPLL的原理结构，概括为： 获取铯钟和驾驭后氢钟的时间偏差；用环路滤波器对该时间偏差进行滤波，滤除高频噪声； 然后用滤波后的时间偏差对氢钟的时间偏差进行调整，得到驾驭后氢钟的时间偏差。图1 中：Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢 钟的输出信号，LF代表环路滤波器。
[0010]在Z域中，把DPLL的环路滤波器LF的传递函数记为G(z);于是，该DPLL的输出 为：
[0012] 显然，
分别为DPLL的系统传递函数和误 差传递函数；Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号Steered Hm代表 被驾驭氢钟的输出信号。
[0013] 通过z = e^f'T(其中，f代表频率，T代表采样周期）的近似变化，将该DPLL等 价为相应的APLL ;于是，式（2)用SSB相位噪声的形式来表示，写成：
[0015] 式中，符号L代表SSB相位噪声，符号G代表传递函数；下标Cs代表铯钟的输出信 号，符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢钟的输出信号；
[0016] S. 1. 2根据步骤S. 1. 1的系统传递函数和误差传递函数，推导等价于二状态变量 Kalman滤波器的DPLL的系统传递函数和误差传递函数的表达式；
[0017] 最终，DPLL的系统传递函数表示为：
[0019] 其中，心为环路增益。Ks:1和Ks22分别为Kalman滤波器的稳态增益；T为采样间 隔。
[0020] 该DPLL的误差传递函数表示为：
[0022] 至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传 递函数。由式（22)和式（23):参数（KsmKsmT，!^)决定了该DPLL的性能。
[0023] S. 2根据步骤S. 1得到的DPLL的系统传递函数和误差传递函数，推导系统传递函 数和误差传递函数中KSldP Ks 22的由过程噪声方差Q22、观测噪声方差R、采样间隔T来表达 的近似值。
[0024] 当Kalman滤波器进入稳态时，得到
[0028]由式（30)和式（35)，Kalman增益由Q22和R之间的比例关系决定。
[0029] 令Q22= Is 2,将式（30)和式（35)代入式（22)和式（23)，得到系统传递函数和误 差传递函数，分别表示为：
[0033] 由式（36)和式（37)可以看出：当〇22=182时，参数〇?，1'，1( (|)决定了0?1^的性 能。把式（36)和式（37)代入式（1)，得到DPLL输入和输出的关系为：
[0035] S. 3根据步骤S. 2得到的系统传递函数和误差传递函数，确定参数（R，T，，优化 DPLL输出的频率稳定度，具体采用以下步骤：
[0036] 首先，铯钟和氢钟的相位噪声分别表不为：
[0040]其中符号L代表SSB相位噪声，&为载波频率，hi为噪声系数，f为边带频率，i用 于指明噪声类型。
[0041] 通过求解方程LjfVLa^f) = 0,可以求得铯钟相位噪声曲线和氢钟相位噪声曲 线交点处的频率，记为f'。
[0042] 第二，将z = e^f'T代入式（36)和式（37)，得到了等价APLL的系统传递函数和 误差传递函数，分别表示为：
[0046]由式（41)和（42)，求得系统传递函数与误差传递函数的幅频响应201og|H(f) |和 201og | He (f) |，并求得这两条曲线交点处的频率，记为f"。
[0047] 第三，对于给定的T，首先给&设一个值（为了方便分析，可以先设为1，2等整数， 然后再逐渐设为其它数值），在该aig条件下，通过改变r的值，来改变f"的值，使f"= f'。此时，两交点的频率相等。此时，在给定T情况下，可以求得相应的R。于是，得到了一 组使两交点的频率相等的参数，记为（R'，IV )。当两交点的频率相等时，DPLL最大程度地 滤除了铯钟的短期噪声和氢钟的长期噪声，使DPLL的输出的相位噪声达到了最优，所以此 时DPLL输出的频率稳定度达到了最优。
[0048] 第四，对已给定的T，改变&的值，使用同样的方法，可以得到许多组最优参数 (R'lV)。这些参数（R'lV)都是最优参数，都保证了 DPLL输出的频率稳定度达到了最 优。
[0049] 具体地，本发明步骤S. 1. 2中，式（22)和式（23)的推导过程如下：
[0050] 对于一个二状态变量的系统，其状态方程表示为：
[0052] 其中，xk和y k为两个状态变量，T为采样间隔，u k为过程噪声；
[0053] 观测方程表示为：
[0054] zk= x k+wk (4)
[0055] 其中，zk为观测量，w k为观测噪声；
[0056] 状态方程和观测方程这两个方程用矩阵的形式表示为：
[0058]其中，sk= [x k yk]T;J k= [0 u k]T;
； H = [1 0]；
[0059]过程噪声和观测噪声的方差分别为：R= E[wk2]， 其中，符号E代表取数学期望；
[0060] Kalman滤波器用下面5个公式进行描述：
[0065] pk k= (I-Kk.H)卟…（10)
[0066] 其中，是状态变量的在k时刻的估计值，Pk， k是估计误差矩阵，4,*-i是在k_l时 刻对k时刻状态变量的单步预测值，Puh是单步预测误差矩阵，K k是Kalman增益矩阵； [0067]能够证明式（5)定义的系统是完全可观测的，因此PkfpPu和Kk都收敛，把Pk，k、 K k的稳态值分别记为：Ps、PjT和Ks ;
[0068]由式（6)和式（9)，当Kalman滤波器进入稳态时，有：
[0070] 其中，毛是状态变量x的在k时刻的估计值，九是状态变量y的在k时刻的估计 值，表示Ks矩阵中的第i行第j列的元素；
[0071]定义：
[0073]将式（12)代入式（11)，得到：
[0075]式（13)在Z域中表示为：
[0077] 其中X，Y，V分别为H v*的Z变换；
[0078] 由式（14)得到：
[0080]由式（14)和式（15)，式（12)在Z域中表示为：
[0082] 其中，Z代表zk的Z变换；
[0083] 由式（16)，得到：
[0085]定义：
[0087]由式（15)、式（17)和式（18)，得到
[0089] 显然，zk是Kalman滤波器的输入，而元:是Kalman滤波器的输出；于是，式（19)表 明了稳态Kalman滤波器输入与输出之间的关系。
[0090] 于是，系统传递函数表示 为：
[0092]由式（20)，Kalman滤波器的传递函数与二阶2类DPLL的闭环系统传递函数是相 同的。于是，Kalman滤波器等价为二阶2类DPLL。
[0093]由于Ksn<< 1，于是式（20)近似写为：
[0095] 为了使DPLL正常工作，引入一个延迟器z'于是DPLL的系统传递函数表示为：
[0097] 其中，K。为环路增益；
[0098] 该DPLL的误差传递函数表示为：
[0100] 至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传 递函数。
[0101] 本发明S. 2中，式（30)和式（35)的推导过程如下：
[0102] 当Kalman滤波器进入稳态时，式（7)写为：
[0104]式⑶写为：
[0106]式（10)写为：
[0108] 式（25)-(27)中，下标ij表示对应矩阵的第i行第j列的元素；
[0109] 由式（25)第3式，式（26)第2式，式（27)第3式，得到：
[0110] Ks2l.Ps-l2=Q22 C28)
[0111] 由式（28)和式（26)第2式，得到：
[0113]由于％ ? A，由式（29)得到
[0115]由式（25)-(27)，得到
[0117]由式（31)，得到
[0119]由于 Ks21 ? T << Ksn，由式（32)，得到：
[0121] 把式（33)代入式（26)，得到：
[0123]由式（34)，得到：
[0125] 本发明可以应用于：用一台原子钟（例如：铯钟）驾驭另一台原子钟（例如：氢 钟）；用一台原子钟驾驭一个数控振荡器（NC0);用一个纸面时间尺度（例如：TA(k))驾驭 一台原子钟；用一个纸面时间尺度（例如：UTC)驾驭另一个纸面时间尺度（例如：TA(k)); 用一台原子钟（例如：UTC(NTSC))驾驭一个纸面时间尺度（例如：BDT)。
[0126] 在上述每一个应用中，都有一个参考原子钟（或时间尺度），和一个被驾驭原子钟 (或时间尺度）。上述这些应用中，参考原子钟（或时间尺度）分别为：一台原子钟（如：铯 钟），一台原子钟，一个纸面时间尺度，一个纸面时间尺度，一台原子钟；被驾驭原子钟（或 时间尺度）分别为：另一台原子钟（如：氢钟），一个数控振荡器，一台原子钟，另一个纸面 时间尺度，一个纸面时间尺度。
[0127] 本发明的原理和步骤在所有这些应用中都是相同的。技术方案中以用铯钟驾驭氢 钟为例来描述本发明。实际上，只要把铯钟替换为其它参考原子钟（或时间尺度），把氢钟 替换为其它被驾驭原子钟（或时间尺度），就可以认为是本发明在上述其它应用背景下的 应用。
[0128] 本发明通过一个等价于二状态变量Kalman滤波器的DPLL，用铯钟对氢钟进行驾 驭。本发明推导了该DPLL的闭环系统传递函数，证明该DPLL为二阶2类DPLL。本发明推 导了由过程噪声方差和观测噪声方差比值近似决定的Kalman增益。于是，由DPLL的闭环传 递函数可知，该比值，以及环路增益，采样间隔这3个参数决定了该DPLL的性能。进一步， 本发明给出了一个简单有效的选取参数的方法。通过等效变换，得到了等效于该DPLL的模 拟锁相环（APLL)。在频域中，通过调整参数，使该等价APLL的系统传递函数和误差传递函 数的交点频率等于氢钟和铯钟的SSB相位噪声曲线的交点频率，最终获取最优参数。
[0129] 本发明的有益效果为：
[0130] 1.可以有效地用铯钟对氢钟进行驾驭，驾驭性能（驾驭误差、频率稳定度）非常 好；
[0131] 2.本发明的参数选取方法相对比较简单。
【附图说明】
[0132] 图1图解说明本发明的实现原理；
[0133] 图2图解说明本发明的算法原理；
[0134] 图3图解说明氢钟和铯钟的SSB相位噪声曲线
[0135] 图4图解说明DPLL的系统传递函数和误差传递函数的幅频响应；
[0136] 图5图解说明氢钟、铯钟和被驾驭氢钟的SSB相位噪声曲线；
[0137] 图6图解说明当氢钟没有频漂时，氢钟、铯钟和被驾驭氢钟的时差；
[0138] 图7图解说明当氢钟没有频漂时，氢钟、铯钟和被驾驭氢钟的Allan偏差；
[0139] 图8图解说明当氢钟有频漂时，氢钟、铯钟和被驾驭氢钟的Allan偏差；
[0140] 图9图解说明当氢钟没有频漂时的驾驭误差；
[0141]
【具体实施方式】
[0142] 下面结合附图对本发明的【具体实施方式】作出说明。
[0143] 本发明一种等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法，包 括以下步骤。
[0144] S. 1推导等价于二状态变量Kalman滤波器的DPLL的传递函数；
[0145] S. 1. 1本发明的实现原理图如图1所示。在图1中描述了使用DPLL，用铯钟来驾 驭氢钟的原理。图中，符号Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号 Steered Hm代表被驾驭氢钟的输出信号，LF代表环路滤波器。
[0146] 在Z域中分析该DPLL的性能，把该DPLL的环路滤波器LF的传递函数记为G (z)。 于是，该DPLL的输出为：
[0148] 显然，
_分别为该DPLL的系统传递函数和 误差传递函数。
[0149] 通过z = e^f'T(其中，f代表频率，T代表采样周期）的近似变化，可以将该DPLL 等价为相应的APLL。于是，式（2)用SSB相位噪声的形式来表示，写成：
[0151]式中，符号L代表SSB相位噪声，符号G代表传递函数；下标中符号Cs代表铯钟的 输出信号，符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢钟的输出信号；
[0152] S. 1. 2根据步骤S. 1. 1的系统传递函数和误差传递函数，推导等价于二状态变量 Kalman滤波器的DPLL的系统传递函数和误差传递函数的表达式。
[0153] 对于一个二状态变量的系统，其状态方程表示为：
[0155] 其中，xk和y k为两个状态变量，T为采样间隔，u k为过程噪声。
[0156] 观测方程表示为：
[0157] zk= x k+wk (4)
[0158] 其中，zk为观测量，wk为观测噪声。
[0159] 这两个方程用矩阵的形式表示为：
[0161]其中，sk= [x k yk]T;J k= [0 u k]T;
; H = [1 0]。
[0162]过程噪声和观测噪声的方差分别为：R= E[wk2]，
其中，符号E代表取数学期望。
[0163] Kalman滤波器可以用下面5个步骤进行描述：
[0168] Pkk= (I-Kk.H) (10)
[0169] 其中，l:是状态变量的在k时刻的估计值，Pk， k是估计误差矩阵，<是在k_l时 刻对k时刻状态变量的单步预测值，Puh是单步预测误差矩阵，K k是Kalman增益矩阵。
[0170] 可以证明式（5)定义的系统是完全可观测的，因此和Kk都收敛。把P k，k、 K k的稳态值分别记为：Ps、Ps _和Ks。
[0171]由式（6)和式（9)，当Kalman滤波器进入稳态时，有：
[0173] 其中，4是状态变量x的在k时刻的估计值，A是状态变量y的在k时刻的估计 值，表示Ks矩阵中的第i行第j列的元素。
[0174]定义：
[0176]将式（12)代入式（11)，得到：
[0178]式（13)在Z域中表示为：
[0180] 其中X，Y, V分别为H v*的Z变换.
[0181] 由式（14)可以得到：
[0183]由式（14)和式（15)，式（12)在Z域中表示为：
[0185] 其中，Z代表zk的Z变换。
[0186] 由式（16)，得到：
)
[0188]定义：
[0190]由式（15)、式（17)和式（18)，得到
[0192] 显然，zk是Kalman滤波器的输入，而元是Kalman滤波器的输出。于是，式（19)表 明了稳态Kalman滤波器输入与输出之间的关系。
[0193] 于是，系统传递函数表示为：
[0195] 由式（20)，Kalman滤波器的传递函数与二阶2类DPLL的闭环系统传递函数是相 同的。于是，Kalman滤波器可以等价为二阶2类DPLL。
[0196] 通过实验，发现Ksn<< 1，于是式（20)可以近似写为：
[0198]为了使DPLL正常工作，引入一个延迟器z'于是DPLL的系统传递函数表示为：
[0200] 其中，K。为环路增益。
[0201] 该DPLL的误差传递函数表示为：
[0203] 至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传 递函数。由式（22)和式（23)可以看出：参数（KsmKsmT，!^)决定了该DPLL的性能。
[0204] 图2描述了该DPLL的算法原理。从图2可以看出，在每一个采样 间隔内，对氢钟 的控制量在Z域中表示为
。于是，在时域中，被驾 驭氢钟的时差和氢钟的时差的关系表示为：
[0206] 式（25)中，Err (j)为第j时刻铯钟与被驾驭氢钟之间的时差Cs (j) -Hmst_d (j)。
[0207] S. 2根据步骤S. 1得到的DPLL的系统传递函数和误差传递函数，推导系统传递函 数和误差传递函数中KSlJP Ks 22的由Q 22、R、T来表达的近似值。
[0208] 实际上，运行Kalman滤波器，待Kalman滤波器进入稳态时，能够得到 真实值。本步骤的意义在于可以不运行Kalman滤波器，通过Q22、R、T的值来计算得到K Sll 和1(822的近似值。仿真实验表明：近似值与真实值的误差很小。于是，可以用近似值来代替 真实值，将KsjP Ks 22的近似值的表达式代入式（23)和式（24)，得到由（Q 22, R，T，心）进行 表达的系统传递函数和误差传递函数。于是，参数（Q22，R，T，ig决定了该DPLL的性能。
[0209] 当Kalman滤波器进入稳态时，式（7)可以写为：
[0211] 式⑶可以写为：
[0213]式（10)可以写为：
[0215] 式（25)_(27)中，下标ij表示对应矩阵的第i行第j列的元素。
[0216] 由式（25)第3式，式（26)第2式，式（27)第3式，得到：
[0218]由式（28)和式（26)第2式，得到：
[0220] 由于％ ?及，由式（29)得到
[0222]由式（25)-(27)，得到
[0224]由式（31)，得到
[0226]由于 Ks21 ? T < < Ksn，由式（32)，得到：
[0228]把式（33)代入式（26)，得到：
[0230]由式（34)，得到：
[0232] 由式（30)和式（35)，可以看出：Kalman增益近似由Q22和R之间的比例关系决定。
[0233] 仿真实验表明这种近似的误差很小。例如：当Q22= ls2, R = 6X102°s2时，由式 (30)和式（35)，计算得到近似值 Ksn~ 5. 42X10 _4和 Ks22~ 4. 08X10 ―11;而运行 Kalman 滤波器，可以得到真实值Ksn= 5. 32X 1(T4和Ks 22= 3. 95X KT11;显然，近似值与真实值之 间的误差小于3%。于是，可以用式（30)和式（35)的计算得到的近似值来近似代替真实 值。
[0234] 令Q22= Is 2,将式（30)和式（35)代入式（22)和式（23)，得到系统传递函数和误 差传递函数，分别表示为：
[0238] 由式（36)和式（37)可以看出：当〇22=182时，参数〇?，1'，1( (|)决定了0?1^的性 能。把式（36)和式（37)代入式（1)，得到DPLL输入和输出的关系为：
[0240] S. 3根据步骤S. 2得到的系统传递函数和误差传递函数，确定参数（R，T，，优化 DPLL输出的频率稳定度。
[0241] 频率稳定度可以在时域中通过Allan方差来表示，也可以通过SSB相位噪声在频 域中进行表示。在频域中比在时域中更容易确定参数。通过变换z = e^f'T将DPLL等效 为一个等价的APLL。于是可以分析该等价APLL在频域的性能。
[0242] 在频域中，选取参数（R，T，K。）的目的主要是优化该等价APLL输出信号的SSB相 位噪声性能，方法为：通过选取参数（R，T，ig使系统传递函数与误差传递函数交点处的频 率等于铯钟与氢钟SSB相位噪声曲线交点处的频率。
[0243] 首先，铯钟和氢钟的相位噪声分别表示为：
[0247] 其中&为载波频率，h i为噪声系数，f为边带频率，i用于指明噪声类型。
[0248] 通过求解方程LjfVLa^f) = 0,可以求得铯钟相位噪声曲线和氢钟相位噪声曲 线交点处的频率，记为f'。
[0249] 第二，将z = e^f'T代入式（36)和式（37)，得到了等价APLL的系统传递函数和 误差传递函数，分别表示为：
[0253] 由式（41)和（42)，可以求得系统传递函数与误差传递函数的幅频响应 201og|H(f) |和201og|He(f) |，并求得这两条曲线交点处的频率，记为f"。
[0254] 第三，对于给定的T，首先给&设一个值，在该（T，K J条件下，通过改变R的值，来 改变f"的值，使f"=f'。此时，两交点的频率相等。此时，在给定T情况下，可以求得相应 的R。于是，得到了一组使两交点的频率相等的参数，记为0?'，&'）。当两交点的频率相 等时，DPLL最大程度地滤除了铯钟的短期噪声和氢钟的长期噪声，使DPLL的输出的相位噪 声达到了最优，所以此时可以认为DPLL输出的频率稳定度达到了最优。
[0255] 第四，对已给定的T，改变&的值，使用同样的方法，可以得到许多组最优参数 (R'lV)。这些参数（R'lV)都是最优参数，都保证了 DPLL输出的频率稳定度达到了最 优，但是实际上的驾驭性能还是有细微的区别。通过仿真实验，可以分析不同组（R'，心'） 条件下的驾驭性能。
[0256] 通过仿真生成一台氢钟和一台铯钟。该氢钟和铯钟主要包含频率白噪声和频率随 机游走噪声。氢钟的噪声系数为h ani)(l= 6Xl(T2f^Ph ani)_2= 2X1CT36,铯钟的噪声系数为 h(Cs)0= 3X 10 _23和h (Cs)_2= 3. 33X 10 _38。氛钟和艳钟的序列各有10000个采样点，采样间 隔为 T = 3600s。
[0257] 图2画出了氢钟和铯钟当采样频率为fQ= 10MHz时的理论SSB相位噪声曲线。从 图2中可以看出交点处的频率约为10_6 6Hz。首先选取心=12,发现当R = 6X10 2(1时，幅 频响应201og|H(f) |和201og|He(f) |曲线的交点处频率近似等于KT6_6Hz。于是，获取 了一组参数（R'，心'）=(6X10 2°s2, 12)。通过同样的方法，获取了其它一些参数类似于 (3X 1022s2, 30)，（5X1017s2, 2)等。图 3 画出了当（R'，K。'）= (6X102°s2, 12)时的系统传 递函数和误差传递函数的幅频响应。
[0258] 由式（3)，式（41)和式（42)，可以计算得到在参数（R'，K。'）和T = 3600s的情 况下，等价APLL输出的理论SSB相位噪声steOTed(f)。回到Z域和式（36)所示的系统传 递函数。在第1种情况下，当（R，，V ) = (5X1017s2,2)时，极点为0.9968±lX10_8i ;在 第2种情况下，当（R'，V ) = (6 X 102°s2, 12)时，极点为0. 9997和0. 9938 ;在第3种情况 下，当（R，，K。'）= (3X 1022 s2, 30)时，极点为0? 9999和0? 9940。在这3种情况下，由于极 点都在单位圆内，所以系统是稳定的。第1种情况下，系统为欠阻尼控制系统。第2和第3 种情况下，系统为过阻尼控制系统，但是第2种情况下系统更灵敏。综上，选取第2种情况 的参数（R，，V ) = (6X 102〇s2, 12)。
[0259] 图4展示了当参数（R'，K。'）= (6X102°s2，12)时的氢钟、铯钟以及驾驭后氢钟的 理论SSB相位噪声。从图4中可以看出，驾驭后氢钟的SSB相位噪声性能非常好。
[0260] 由式（38)，参数（R'，K。'）= (6X102°s2，12)，以及生成的氢钟和铯钟序列，可以得 到DPLL的输出序列，即驾驭后氢钟的序列。
[0261] 图5展示了前4800个数据点的氢钟、铯钟、驾驭后氢钟的时差。图6展示了氢钟、 铯钟、驾驭后氢钟的A1 lan偏差。从图5中可以看出，该DPLL可以有效地使用铯钟驾驭氢 钟，驾驭误差较小。从图6可以看出：驾驭后氢钟的频率稳定度体现了氢钟的短期稳定度和 铯钟的长期稳定度，其短期稳定度和长期稳定度相比氢钟和铯钟分别只有很小的恶化。
[0262] 在原有的氢钟序列上加入幅度为3X 10_21^的频漂。图7展示了当氢钟存在频漂 时，氢钟、铯钟、驾驭后氢钟的Allan偏差。对比图6和图7,可以看出频漂并不会影响驾驭 后氣钟的频率稳定度。
[0263] 为了能更清楚地展示驾驭误差（即驾驭后氢钟相对铯钟的偏差），重新生成了氢 钟和铯钟序列，每个序列的参数和原氢钟和原铯钟的参数相同，共有200000个采样点。通 过使用相同的DPLL方法，获得了驾驭后氢钟的序列。图8展示了驾驭误差。从图8中可以 看出：最大驾驭误差为30ns，标准偏差为8. 99ns。实际上，驾驭误差是由于未跟踪的氢钟和 铯钟的相位噪声引起的，并且驾驭误差的标准偏差，即抖动，可以通过控制理论精确地计算 得到。因此可以推测：在通过UTC来驾驭纸面时间这类应用中，由于纸面时间和UTC都更 稳定，包含了更少的等效SSB相位噪声分量，因此驾驭误差也会相应地更小。
[0264] 在原有的氢钟序列上加入幅度为SXKT21。 1的频漂。图9展示了当氢钟含有频漂 时的驾驭误差。对比图8和图9,可以看出：当氢钟存在频漂时，会造成一个稳态跟踪误差， 幅度大约为20ns。
[0265] 综上，本发明可以有效地用参考原子钟（或时间尺度）来驾驭被驾驭原 子钟（或 时间尺度），驾驭误差和驾驭后频率稳定度的性能都较优。被驾驭原子钟（或时间尺度）的 频漂会引起一个稳态跟踪误差，这个稳态跟踪误差的值可以精确地计算得到，并通过其它 算法来进行补偿。
[0266] 综上所述，虽然本发明已以较佳实施例（通过铯钟驾驭氢钟）揭露如上，然其并非 用以限定本发明，任何本领域普通技术人员，在不脱离本发明的精神和范围内，当可作各种 更动与润饰，因此本发明的保护范围当视权利要求书界定的范围为准。
【主权项】
1. 一种等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法，其特征在于， 采用以下步骤： S. 1推导等价于二状态变量Kalman滤波器的DPLL的传递函数； S. I. 1使用DPLL，用铯钟来驾驭氢钟；DPLL的原理结构为：获取铯钟和驾驭后氢钟的时 间偏差；用环路滤波器对该时间偏差进行滤波，滤除高频噪声；然后用滤波后的时间偏差 对氢钟的时间偏差进行调整，得到驾驭后氢钟的时间偏差，其中：Cs代表铯钟的输出信号， 符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢钟的输出信号，LF代表环路 滤波器； 在Z域中，把DPLL的环路滤波器LF的传递函数记为G(Z);于是，该DPLL的输出为：(1) 显然：分别为DPLL的系统传递函数和误差传递 函数；其中：符号Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号Steered Hm 代表被驾驭氢钟的输出信号； 通过z = #π?'τ的近似变化，其中，f代表频率，T代表采样周期；将该DPLL等价为相 应的APLL ;于是，式（2)用SSB相位噪声的形式来表示，写成：式中，符号L代表SSB相位噪声，符号G代表传递函数；下标中符号Cs代表铯钟的输出 信号，符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢钟的输出信号； S. 1.2根据步骤S. 1.1的系统传递函数和误差传递函数，推导等价于二状态变量 Kalman滤波器的DPLL的系统传递函数和误差传递函数的表达式； 最终，DPLL的系统传递函数表示为：其中，Ktl为环路增益；Ks Ks 22分别为Kalman滤波器的稳态增益；T为采样间隔； 该DPLL的误差传递函数表示为：至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传递函 数；由式（22)和式（23):参数（Ksn，Ks22, T，Ktl)决定了该DPLL的性能； S. 2根据步骤S. 1得到的DPLL的系统传递函数和误差传递函数，推导系统传递函数和 误差传递函数中Ks1JP Ks 22的由过程噪声方差Q22、观测噪声方差R、采样间隔T来表达的近 似值； 当Kalman滤波器进入稳态时，得到由式（30)和式（35)，Kalman增益由Q22和R之间的比例关系决定； 令Q22= Is 2,将式（30)和式（35)代入式（22)和式（23)，得到系统传递函数和误差传 递函数，分别表示为：由式（36)和式（37)可以看出：当Q22= Is 2时，参数（RJJtl)决定了 DPLL的性能；把 式（36)和式（37)代入式（1)，得到DPLL输入和输出的关系为：S. 3根据步骤S. 2得到的系统传递函数和误差传递函数，确定参数（R，T，Ktl)，优化DPLL 输出的频率稳定度，具体采用以下步骤： 首先，铯钟和氢钟的相位噪声分别表不为：其中符号L代表SSB相位噪声，&为载波频率，h i为噪声系数，f为边带频率，i用于指 明噪声类型； 通过求解方程L&(f)-LHm(f) = 0,可以求得铯钟相位噪声曲线和氢钟相位噪声曲线交 点处的频率，记为f' ； 第二，将z = #π?'τ代入式（36)和式（37)，得到了等价APLL的系统传递函数和误差 传递函数，分别表示为：和由式（41)和（42)，求得系统传递函数与误差传递函数的幅频响应201og|H(f) I和 201og I He (f) I，并求得这两条曲线交点处的频率，记为f" ； 第三，对于给定的T，首先给Ktl设一个值（，在该（T，K ^条件下，通过改变R的值，来改 变f"的值，使f" = f' ；此时，两交点的频率相等；此时，在给定T情况下，可以求得相应的 R;于是，得到了一组使两交点的频率相等的参数，记为（R'，Kc/ );当两交点的频率相等时， DPLL最大程度地滤除了铯钟的短期噪声和氢钟的长期噪声，使DPLL的输出的相位噪声达 到了最优，所以此时DPLL输出的频率稳定度达到了最优； 第四，对已给定的T，改变Ktl的值，使用步骤三中相同的方法，能够得到多组最优参数 (R'，Kc/ )，这些参数（R'，Kc/ )都是最优参数，都保证了 DPLL输出的频率稳定度达到了最 优。2.根据权利要求1所述的等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭 方法，其特征在于，步骤S. 1. 2中等价于二状态变量Kalman滤波器的DPLL的系统传递函数 和误差传递函数的推导过程如下： 对于一个二状态变量的系统，其状态方程表示为：(3) 其中，xk和yk为两个状态变量，T为采样间隔，Uk为过程噪声； 观测方程表示为： zk= x k+wk ⑷ 其中，zk为观测量，wk为观测噪声； 状态方程和观测方程这两个方程用矩阵的形式表示为：(5) 其中：过程噪声和观测噪声的方差分别为：其中，符号E代表取数学期望； Kalman滤波器用下面5个公式进行描述：其中，I是状态变量的在k时刻的估计值，Ptk是估计误差矩阵，&?是在k-1时刻对 k时刻状态变量的单步预测值，？0_1是单步预测误差矩阵，K k是Kalman增益矩阵； 能够证明式（5)定义的系统是完全可观测的，因此Ptk^Ptk和Kk都收敛，把P k^Pk^ 和Kk的稳态值分别记为：Ps、Ps _和Ks ; 由式（6)和式（9)，当Kalman滤波器进入稳态时，有：(11) 其中，毛是状态变量X的在k时刻的估计值，λ是状态变量y的在k时刻的估计值，Ksij表示Ks矩阵中的第i行第j列的元素； 定义： vA =(-Α -·ν?->ν?*Γ) (12) 将式（12)代入式（11)，得到：(13) 式（13)在Z域中表示为：(14) 其中X，Y，V分别为H 的ζ变换； 由式（14)得到：(15) 由式（14)和式（15)，式（12)在Z域中表示为：(16) 其中，Z代表Zk的Z变换； 由式（16)，得到：由式（15)、式（17)和式（18)，得到显然，Zk是Kalman滤波器的输入，而毛是Kalman滤波器的输出；于是，式（19)表明了 稳态Kalman滤波器输入与输出之间的关系； 于是，系统传递函数表示为：由式（20)，Kalman滤波器的传递函数与二阶2类DPLL的闭环系统传递函数是相同的； 于是，Kalman滤波器等价为二阶2类DPLL ; 显然，zk是Kalman滤波器的输入，而L是Kalman滤波器的输出；于是，式（19)表明了 稳态Kalman滤波器输入与输出之间的关系； 于是，系统传递函数表示为：由式（20)，Kalman滤波器的传递函数与二阶2类DPLL的闭环系统传递函数是相同的； 于是，Kalman滤波器等价为二阶2类DPLL ; 由于Ksn<< 1，于是式（20)近似写为：为了使DPLL正常工作，引入一个延迟器厂\于是DPLL的系统传递函数表示为：(22) 其中，Ktl为环路增益； 该DPLL的误差传递函数表示为：(23) 至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传递函 数。3.根据权利要求2所述的等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭 方法，其特征在于，S. 2中式（30)和式（35)的推导过程如下： 当Kalman滤波器进入稳态时，式（7)写为：式（25)-(27)中，下标ij表示对应矩阵的第i行第j列的元素； 由式（25)第3式，式（26)第2式，式（27)第3式，得到： Ks2l*Psl2 =Q22 (28 ) 由式（28)和式（26)第2式，得到：(29) 由于〃vn<^，由式（29)得到(30) 由式（25)-(27)，得到(31) 由式（31)，得到(32) 由于Ks21 · T << Ksn，由式（32)，得到：(33) 把式（33)代入式（26)，得到：(34) 由式（34)，得到：(35)
【专利摘要】一种等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法，通过一个等价于二状态变量Kalman滤波器的DPLL，用铯钟对氢钟进行驾驭。首先推导了该DPLL的闭环系统传递函数，证明该DPLL为二阶2类DPLL。然后推导了由过程噪声方差和观测噪声方差比值近似决定的Kalman增益。于是，由DPLL的闭环传递函数可知，该比值，以及环路增益，采样间隔这3个参数决定了该DPLL的性能。进一步，本发明给出了一个简单有效的选取参数的方法。通过等效变换，得到了等效于该DPLL的模拟锁相环(APLL)。在频域中，通过调整参数，使该等价APLL的系统传递函数和误差传递函数的交点频率等于氢钟和铯钟的SSB相位噪声曲线的交点频率，最终获取最优参数。其有效地用铯钟对氢钟进行驾驭，驾驭性能非常好。
【IPC分类】G04F5/14
【公开号】CN104898404
【申请号】CN201510353792
【发明人】龚航, 伍贻威, 黄新明, 刘文祥, 欧钢, 朱祥维, 李垣陵, 黄仰博
【申请人】中国人民解放军国防科学技术大学
【公开日】2015年9月9日
【申请日】2015年6月25日
【技术领域】
[0001] 本发明涉及时间频率、信号处理领域，具体的说是设计了一种等价于二状态变量 Kalman滤波器的数字锁相环原子钟驾驭方法。
【背景技术】
[0002] 原子钟驾驭技术在守时实验室和卫星导航系统中发挥着重要作用。对原子钟进行 驾驭的主要目的有两个：一是使被驾驭原子钟与参考原子钟时间同步，尽可能缩小它们之 间的时间偏差；二是提升被驾驭原子钟的长期稳定度。
[0003] 典型的原子钟驾驭方法包括两大类：开环驾驭方法和闭环驾驭方法。开环驾驭方 法的核心是设计合理的预测算法。闭环驾驭方法主要包括：线性二次型（LG?控制方法、开 关（Bang-Bang)控制方法和数字锁相环（DPLL)方法。但是，这些闭环驾驭方法的参数是不 容易选取的，一般都要针对每一个具体应用，通过大量仿真实验对大量参数进行比较后，选 出一组最优参数。另外，对于DPLL方法，如果参数选取不当，不但无法保证驾驭性能，还会 引起控制系统不稳定。
【发明内容】
[0004] 为解决现有技术中存在的问题，本发明的目的在于提供一种等价于二状态变量 Kalman滤波器的数字锁相环原子钟驾驭方法。
[0005] 本发明的原理是：通过一个等价于二状态变量Kalman滤波器的DPLL，用铯钟对氢 钟进行驾驭。本发明推导了该DPLL的闭环系统传递函数，证明该DPLL为二阶2类DPLL。 本发明推导了由过程噪声方差和观测噪声方差比值近似决定的Kalman增益。于是，由DPLL 的闭环传递函数可知，该比值，以及环路增益，采样间隔这3个参数决定了该DPLL的性能。 进一步，本发明给出了一个简单有效的选取参数的方法。通过等效变换，得到了等效于该 DPLL的模拟锁相环（APLL)。在频域中，通过调整参数，使该等价APLL的系统传递函数和误 差传递函数的交点频率等于氢钟和铯钟的SSB相位噪声曲线的交点频率，最终获取最优参 数。简而言之，使用参考原子钟（或时间尺度）来驾驭被驾驭原子钟（或时间尺度），使被 驾驭原子钟（或时间尺度）与参考原子钟（或时间尺度）保持同步。
[0006] 本发明采用的技术方案为：
[0007] -种等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法，采用以下 步骤：
[0008] S. 1推导等价于二状态变量Kalman滤波器的DPLL的传递函数；
[0009] S. 1. 1使用DPLL，用铯钟来驾驭氢钟。图1描述了该DPLL的原理结构，概括为： 获取铯钟和驾驭后氢钟的时间偏差；用环路滤波器对该时间偏差进行滤波，滤除高频噪声； 然后用滤波后的时间偏差对氢钟的时间偏差进行调整，得到驾驭后氢钟的时间偏差。图1 中：Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢 钟的输出信号，LF代表环路滤波器。
[0010]在Z域中，把DPLL的环路滤波器LF的传递函数记为G(z);于是，该DPLL的输出 为：
[0012] 显然，
分别为DPLL的系统传递函数和误 差传递函数；Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号Steered Hm代表 被驾驭氢钟的输出信号。
[0013] 通过z = e^f'T(其中，f代表频率，T代表采样周期）的近似变化，将该DPLL等 价为相应的APLL ;于是，式（2)用SSB相位噪声的形式来表示，写成：
[0015] 式中，符号L代表SSB相位噪声，符号G代表传递函数；下标Cs代表铯钟的输出信 号，符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢钟的输出信号；
[0016] S. 1. 2根据步骤S. 1. 1的系统传递函数和误差传递函数，推导等价于二状态变量 Kalman滤波器的DPLL的系统传递函数和误差传递函数的表达式；
[0017] 最终，DPLL的系统传递函数表示为：
[0019] 其中，心为环路增益。Ks:1和Ks22分别为Kalman滤波器的稳态增益；T为采样间 隔。
[0020] 该DPLL的误差传递函数表示为：
[0022] 至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传 递函数。由式（22)和式（23):参数（KsmKsmT，!^)决定了该DPLL的性能。
[0023] S. 2根据步骤S. 1得到的DPLL的系统传递函数和误差传递函数，推导系统传递函 数和误差传递函数中KSldP Ks 22的由过程噪声方差Q22、观测噪声方差R、采样间隔T来表达 的近似值。
[0024] 当Kalman滤波器进入稳态时，得到
[0028]由式（30)和式（35)，Kalman增益由Q22和R之间的比例关系决定。
[0029] 令Q22= Is 2,将式（30)和式（35)代入式（22)和式（23)，得到系统传递函数和误 差传递函数，分别表示为：
[0033] 由式（36)和式（37)可以看出：当〇22=182时，参数〇?，1'，1( (|)决定了0?1^的性 能。把式（36)和式（37)代入式（1)，得到DPLL输入和输出的关系为：
[0035] S. 3根据步骤S. 2得到的系统传递函数和误差传递函数，确定参数（R，T，，优化 DPLL输出的频率稳定度，具体采用以下步骤：
[0036] 首先，铯钟和氢钟的相位噪声分别表不为：
[0040]其中符号L代表SSB相位噪声，&为载波频率，hi为噪声系数，f为边带频率，i用 于指明噪声类型。
[0041] 通过求解方程LjfVLa^f) = 0,可以求得铯钟相位噪声曲线和氢钟相位噪声曲 线交点处的频率，记为f'。
[0042] 第二，将z = e^f'T代入式（36)和式（37)，得到了等价APLL的系统传递函数和 误差传递函数，分别表示为：
[0046]由式（41)和（42)，求得系统传递函数与误差传递函数的幅频响应201og|H(f) |和 201og | He (f) |，并求得这两条曲线交点处的频率，记为f"。
[0047] 第三，对于给定的T，首先给&设一个值（为了方便分析，可以先设为1，2等整数， 然后再逐渐设为其它数值），在该aig条件下，通过改变r的值，来改变f"的值，使f"= f'。此时，两交点的频率相等。此时，在给定T情况下，可以求得相应的R。于是，得到了一 组使两交点的频率相等的参数，记为（R'，IV )。当两交点的频率相等时，DPLL最大程度地 滤除了铯钟的短期噪声和氢钟的长期噪声，使DPLL的输出的相位噪声达到了最优，所以此 时DPLL输出的频率稳定度达到了最优。
[0048] 第四，对已给定的T，改变&的值，使用同样的方法，可以得到许多组最优参数 (R'lV)。这些参数（R'lV)都是最优参数，都保证了 DPLL输出的频率稳定度达到了最 优。
[0049] 具体地，本发明步骤S. 1. 2中，式（22)和式（23)的推导过程如下：
[0050] 对于一个二状态变量的系统，其状态方程表示为：
[0052] 其中，xk和y k为两个状态变量，T为采样间隔，u k为过程噪声；
[0053] 观测方程表示为：
[0054] zk= x k+wk (4)
[0055] 其中，zk为观测量，w k为观测噪声；
[0056] 状态方程和观测方程这两个方程用矩阵的形式表示为：
[0058]其中，sk= [x k yk]T;J k= [0 u k]T;
； H = [1 0]；
[0059]过程噪声和观测噪声的方差分别为：R= E[wk2]， 其中，符号E代表取数学期望；
[0060] Kalman滤波器用下面5个公式进行描述：
[0065] pk k= (I-Kk.H)卟…（10)
[0066] 其中，是状态变量的在k时刻的估计值，Pk， k是估计误差矩阵，4,*-i是在k_l时 刻对k时刻状态变量的单步预测值，Puh是单步预测误差矩阵，K k是Kalman增益矩阵； [0067]能够证明式（5)定义的系统是完全可观测的，因此PkfpPu和Kk都收敛，把Pk，k、 K k的稳态值分别记为：Ps、PjT和Ks ;
[0068]由式（6)和式（9)，当Kalman滤波器进入稳态时，有：
[0070] 其中，毛是状态变量x的在k时刻的估计值，九是状态变量y的在k时刻的估计 值，表示Ks矩阵中的第i行第j列的元素；
[0071]定义：
[0073]将式（12)代入式（11)，得到：
[0075]式（13)在Z域中表示为：
[0077] 其中X，Y，V分别为H v*的Z变换；
[0078] 由式（14)得到：
[0080]由式（14)和式（15)，式（12)在Z域中表示为：
[0082] 其中，Z代表zk的Z变换；
[0083] 由式（16)，得到：
[0085]定义：
[0087]由式（15)、式（17)和式（18)，得到
[0089] 显然，zk是Kalman滤波器的输入，而元:是Kalman滤波器的输出；于是，式（19)表 明了稳态Kalman滤波器输入与输出之间的关系。
[0090] 于是，系统传递函数表示 为：
[0092]由式（20)，Kalman滤波器的传递函数与二阶2类DPLL的闭环系统传递函数是相 同的。于是，Kalman滤波器等价为二阶2类DPLL。
[0093]由于Ksn<< 1，于是式（20)近似写为：
[0095] 为了使DPLL正常工作，引入一个延迟器z'于是DPLL的系统传递函数表示为：
[0097] 其中，K。为环路增益；
[0098] 该DPLL的误差传递函数表示为：
[0100] 至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传 递函数。
[0101] 本发明S. 2中，式（30)和式（35)的推导过程如下：
[0102] 当Kalman滤波器进入稳态时，式（7)写为：
[0104]式⑶写为：
[0106]式（10)写为：
[0108] 式（25)-(27)中，下标ij表示对应矩阵的第i行第j列的元素；
[0109] 由式（25)第3式，式（26)第2式，式（27)第3式，得到：
[0110] Ks2l.Ps-l2=Q22 C28)
[0111] 由式（28)和式（26)第2式，得到：
[0113]由于％ ? A，由式（29)得到
[0115]由式（25)-(27)，得到
[0117]由式（31)，得到
[0119]由于 Ks21 ? T << Ksn，由式（32)，得到：
[0121] 把式（33)代入式（26)，得到：
[0123]由式（34)，得到：
[0125] 本发明可以应用于：用一台原子钟（例如：铯钟）驾驭另一台原子钟（例如：氢 钟）；用一台原子钟驾驭一个数控振荡器（NC0);用一个纸面时间尺度（例如：TA(k))驾驭 一台原子钟；用一个纸面时间尺度（例如：UTC)驾驭另一个纸面时间尺度（例如：TA(k)); 用一台原子钟（例如：UTC(NTSC))驾驭一个纸面时间尺度（例如：BDT)。
[0126] 在上述每一个应用中，都有一个参考原子钟（或时间尺度），和一个被驾驭原子钟 (或时间尺度）。上述这些应用中，参考原子钟（或时间尺度）分别为：一台原子钟（如：铯 钟），一台原子钟，一个纸面时间尺度，一个纸面时间尺度，一台原子钟；被驾驭原子钟（或 时间尺度）分别为：另一台原子钟（如：氢钟），一个数控振荡器，一台原子钟，另一个纸面 时间尺度，一个纸面时间尺度。
[0127] 本发明的原理和步骤在所有这些应用中都是相同的。技术方案中以用铯钟驾驭氢 钟为例来描述本发明。实际上，只要把铯钟替换为其它参考原子钟（或时间尺度），把氢钟 替换为其它被驾驭原子钟（或时间尺度），就可以认为是本发明在上述其它应用背景下的 应用。
[0128] 本发明通过一个等价于二状态变量Kalman滤波器的DPLL，用铯钟对氢钟进行驾 驭。本发明推导了该DPLL的闭环系统传递函数，证明该DPLL为二阶2类DPLL。本发明推 导了由过程噪声方差和观测噪声方差比值近似决定的Kalman增益。于是，由DPLL的闭环传 递函数可知，该比值，以及环路增益，采样间隔这3个参数决定了该DPLL的性能。进一步， 本发明给出了一个简单有效的选取参数的方法。通过等效变换，得到了等效于该DPLL的模 拟锁相环（APLL)。在频域中，通过调整参数，使该等价APLL的系统传递函数和误差传递函 数的交点频率等于氢钟和铯钟的SSB相位噪声曲线的交点频率，最终获取最优参数。
[0129] 本发明的有益效果为：
[0130] 1.可以有效地用铯钟对氢钟进行驾驭，驾驭性能（驾驭误差、频率稳定度）非常 好；
[0131] 2.本发明的参数选取方法相对比较简单。
【附图说明】
[0132] 图1图解说明本发明的实现原理；
[0133] 图2图解说明本发明的算法原理；
[0134] 图3图解说明氢钟和铯钟的SSB相位噪声曲线
[0135] 图4图解说明DPLL的系统传递函数和误差传递函数的幅频响应；
[0136] 图5图解说明氢钟、铯钟和被驾驭氢钟的SSB相位噪声曲线；
[0137] 图6图解说明当氢钟没有频漂时，氢钟、铯钟和被驾驭氢钟的时差；
[0138] 图7图解说明当氢钟没有频漂时，氢钟、铯钟和被驾驭氢钟的Allan偏差；
[0139] 图8图解说明当氢钟有频漂时，氢钟、铯钟和被驾驭氢钟的Allan偏差；
[0140] 图9图解说明当氢钟没有频漂时的驾驭误差；
[0141]
【具体实施方式】
[0142] 下面结合附图对本发明的【具体实施方式】作出说明。
[0143] 本发明一种等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法，包 括以下步骤。
[0144] S. 1推导等价于二状态变量Kalman滤波器的DPLL的传递函数；
[0145] S. 1. 1本发明的实现原理图如图1所示。在图1中描述了使用DPLL，用铯钟来驾 驭氢钟的原理。图中，符号Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号 Steered Hm代表被驾驭氢钟的输出信号，LF代表环路滤波器。
[0146] 在Z域中分析该DPLL的性能，把该DPLL的环路滤波器LF的传递函数记为G (z)。 于是，该DPLL的输出为：
[0148] 显然，
_分别为该DPLL的系统传递函数和 误差传递函数。
[0149] 通过z = e^f'T(其中，f代表频率，T代表采样周期）的近似变化，可以将该DPLL 等价为相应的APLL。于是，式（2)用SSB相位噪声的形式来表示，写成：
[0151]式中，符号L代表SSB相位噪声，符号G代表传递函数；下标中符号Cs代表铯钟的 输出信号，符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢钟的输出信号；
[0152] S. 1. 2根据步骤S. 1. 1的系统传递函数和误差传递函数，推导等价于二状态变量 Kalman滤波器的DPLL的系统传递函数和误差传递函数的表达式。
[0153] 对于一个二状态变量的系统，其状态方程表示为：
[0155] 其中，xk和y k为两个状态变量，T为采样间隔，u k为过程噪声。
[0156] 观测方程表示为：
[0157] zk= x k+wk (4)
[0158] 其中，zk为观测量，wk为观测噪声。
[0159] 这两个方程用矩阵的形式表示为：
[0161]其中，sk= [x k yk]T;J k= [0 u k]T;
; H = [1 0]。
[0162]过程噪声和观测噪声的方差分别为：R= E[wk2]，
其中，符号E代表取数学期望。
[0163] Kalman滤波器可以用下面5个步骤进行描述：
[0168] Pkk= (I-Kk.H) (10)
[0169] 其中，l:是状态变量的在k时刻的估计值，Pk， k是估计误差矩阵，<是在k_l时 刻对k时刻状态变量的单步预测值，Puh是单步预测误差矩阵，K k是Kalman增益矩阵。
[0170] 可以证明式（5)定义的系统是完全可观测的，因此和Kk都收敛。把P k，k、 K k的稳态值分别记为：Ps、Ps _和Ks。
[0171]由式（6)和式（9)，当Kalman滤波器进入稳态时，有：
[0173] 其中，4是状态变量x的在k时刻的估计值，A是状态变量y的在k时刻的估计 值，表示Ks矩阵中的第i行第j列的元素。
[0174]定义：
[0176]将式（12)代入式（11)，得到：
[0178]式（13)在Z域中表示为：
[0180] 其中X，Y, V分别为H v*的Z变换.
[0181] 由式（14)可以得到：
[0183]由式（14)和式（15)，式（12)在Z域中表示为：
[0185] 其中，Z代表zk的Z变换。
[0186] 由式（16)，得到：
)
[0188]定义：
[0190]由式（15)、式（17)和式（18)，得到
[0192] 显然，zk是Kalman滤波器的输入，而元是Kalman滤波器的输出。于是，式（19)表 明了稳态Kalman滤波器输入与输出之间的关系。
[0193] 于是，系统传递函数表示为：
[0195] 由式（20)，Kalman滤波器的传递函数与二阶2类DPLL的闭环系统传递函数是相 同的。于是，Kalman滤波器可以等价为二阶2类DPLL。
[0196] 通过实验，发现Ksn<< 1，于是式（20)可以近似写为：
[0198]为了使DPLL正常工作，引入一个延迟器z'于是DPLL的系统传递函数表示为：
[0200] 其中，K。为环路增益。
[0201] 该DPLL的误差传递函数表示为：
[0203] 至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传 递函数。由式（22)和式（23)可以看出：参数（KsmKsmT，!^)决定了该DPLL的性能。
[0204] 图2描述了该DPLL的算法原理。从图2可以看出，在每一个采样 间隔内，对氢钟 的控制量在Z域中表示为
。于是，在时域中，被驾 驭氢钟的时差和氢钟的时差的关系表示为：
[0206] 式（25)中，Err (j)为第j时刻铯钟与被驾驭氢钟之间的时差Cs (j) -Hmst_d (j)。
[0207] S. 2根据步骤S. 1得到的DPLL的系统传递函数和误差传递函数，推导系统传递函 数和误差传递函数中KSlJP Ks 22的由Q 22、R、T来表达的近似值。
[0208] 实际上，运行Kalman滤波器，待Kalman滤波器进入稳态时，能够得到 真实值。本步骤的意义在于可以不运行Kalman滤波器，通过Q22、R、T的值来计算得到K Sll 和1(822的近似值。仿真实验表明：近似值与真实值的误差很小。于是，可以用近似值来代替 真实值，将KsjP Ks 22的近似值的表达式代入式（23)和式（24)，得到由（Q 22, R，T，心）进行 表达的系统传递函数和误差传递函数。于是，参数（Q22，R，T，ig决定了该DPLL的性能。
[0209] 当Kalman滤波器进入稳态时，式（7)可以写为：
[0211] 式⑶可以写为：
[0213]式（10)可以写为：
[0215] 式（25)_(27)中，下标ij表示对应矩阵的第i行第j列的元素。
[0216] 由式（25)第3式，式（26)第2式，式（27)第3式，得到：
[0218]由式（28)和式（26)第2式，得到：
[0220] 由于％ ?及，由式（29)得到
[0222]由式（25)-(27)，得到
[0224]由式（31)，得到
[0226]由于 Ks21 ? T < < Ksn，由式（32)，得到：
[0228]把式（33)代入式（26)，得到：
[0230]由式（34)，得到：
[0232] 由式（30)和式（35)，可以看出：Kalman增益近似由Q22和R之间的比例关系决定。
[0233] 仿真实验表明这种近似的误差很小。例如：当Q22= ls2, R = 6X102°s2时，由式 (30)和式（35)，计算得到近似值 Ksn~ 5. 42X10 _4和 Ks22~ 4. 08X10 ―11;而运行 Kalman 滤波器，可以得到真实值Ksn= 5. 32X 1(T4和Ks 22= 3. 95X KT11;显然，近似值与真实值之 间的误差小于3%。于是，可以用式（30)和式（35)的计算得到的近似值来近似代替真实 值。
[0234] 令Q22= Is 2,将式（30)和式（35)代入式（22)和式（23)，得到系统传递函数和误 差传递函数，分别表示为：
[0238] 由式（36)和式（37)可以看出：当〇22=182时，参数〇?，1'，1( (|)决定了0?1^的性 能。把式（36)和式（37)代入式（1)，得到DPLL输入和输出的关系为：
[0240] S. 3根据步骤S. 2得到的系统传递函数和误差传递函数，确定参数（R，T，，优化 DPLL输出的频率稳定度。
[0241] 频率稳定度可以在时域中通过Allan方差来表示，也可以通过SSB相位噪声在频 域中进行表示。在频域中比在时域中更容易确定参数。通过变换z = e^f'T将DPLL等效 为一个等价的APLL。于是可以分析该等价APLL在频域的性能。
[0242] 在频域中，选取参数（R，T，K。）的目的主要是优化该等价APLL输出信号的SSB相 位噪声性能，方法为：通过选取参数（R，T，ig使系统传递函数与误差传递函数交点处的频 率等于铯钟与氢钟SSB相位噪声曲线交点处的频率。
[0243] 首先，铯钟和氢钟的相位噪声分别表示为：
[0247] 其中&为载波频率，h i为噪声系数，f为边带频率，i用于指明噪声类型。
[0248] 通过求解方程LjfVLa^f) = 0,可以求得铯钟相位噪声曲线和氢钟相位噪声曲 线交点处的频率，记为f'。
[0249] 第二，将z = e^f'T代入式（36)和式（37)，得到了等价APLL的系统传递函数和 误差传递函数，分别表示为：
[0253] 由式（41)和（42)，可以求得系统传递函数与误差传递函数的幅频响应 201og|H(f) |和201og|He(f) |，并求得这两条曲线交点处的频率，记为f"。
[0254] 第三，对于给定的T，首先给&设一个值，在该（T，K J条件下，通过改变R的值，来 改变f"的值，使f"=f'。此时，两交点的频率相等。此时，在给定T情况下，可以求得相应 的R。于是，得到了一组使两交点的频率相等的参数，记为0?'，&'）。当两交点的频率相 等时，DPLL最大程度地滤除了铯钟的短期噪声和氢钟的长期噪声，使DPLL的输出的相位噪 声达到了最优，所以此时可以认为DPLL输出的频率稳定度达到了最优。
[0255] 第四，对已给定的T，改变&的值，使用同样的方法，可以得到许多组最优参数 (R'lV)。这些参数（R'lV)都是最优参数，都保证了 DPLL输出的频率稳定度达到了最 优，但是实际上的驾驭性能还是有细微的区别。通过仿真实验，可以分析不同组（R'，心'） 条件下的驾驭性能。
[0256] 通过仿真生成一台氢钟和一台铯钟。该氢钟和铯钟主要包含频率白噪声和频率随 机游走噪声。氢钟的噪声系数为h ani)(l= 6Xl(T2f^Ph ani)_2= 2X1CT36,铯钟的噪声系数为 h(Cs)0= 3X 10 _23和h (Cs)_2= 3. 33X 10 _38。氛钟和艳钟的序列各有10000个采样点，采样间 隔为 T = 3600s。
[0257] 图2画出了氢钟和铯钟当采样频率为fQ= 10MHz时的理论SSB相位噪声曲线。从 图2中可以看出交点处的频率约为10_6 6Hz。首先选取心=12,发现当R = 6X10 2(1时，幅 频响应201og|H(f) |和201og|He(f) |曲线的交点处频率近似等于KT6_6Hz。于是，获取 了一组参数（R'，心'）=(6X10 2°s2, 12)。通过同样的方法，获取了其它一些参数类似于 (3X 1022s2, 30)，（5X1017s2, 2)等。图 3 画出了当（R'，K。'）= (6X102°s2, 12)时的系统传 递函数和误差传递函数的幅频响应。
[0258] 由式（3)，式（41)和式（42)，可以计算得到在参数（R'，K。'）和T = 3600s的情 况下，等价APLL输出的理论SSB相位噪声steOTed(f)。回到Z域和式（36)所示的系统传 递函数。在第1种情况下，当（R，，V ) = (5X1017s2,2)时，极点为0.9968±lX10_8i ;在 第2种情况下，当（R'，V ) = (6 X 102°s2, 12)时，极点为0. 9997和0. 9938 ;在第3种情况 下，当（R，，K。'）= (3X 1022 s2, 30)时，极点为0? 9999和0? 9940。在这3种情况下，由于极 点都在单位圆内，所以系统是稳定的。第1种情况下，系统为欠阻尼控制系统。第2和第3 种情况下，系统为过阻尼控制系统，但是第2种情况下系统更灵敏。综上，选取第2种情况 的参数（R，，V ) = (6X 102〇s2, 12)。
[0259] 图4展示了当参数（R'，K。'）= (6X102°s2，12)时的氢钟、铯钟以及驾驭后氢钟的 理论SSB相位噪声。从图4中可以看出，驾驭后氢钟的SSB相位噪声性能非常好。
[0260] 由式（38)，参数（R'，K。'）= (6X102°s2，12)，以及生成的氢钟和铯钟序列，可以得 到DPLL的输出序列，即驾驭后氢钟的序列。
[0261] 图5展示了前4800个数据点的氢钟、铯钟、驾驭后氢钟的时差。图6展示了氢钟、 铯钟、驾驭后氢钟的A1 lan偏差。从图5中可以看出，该DPLL可以有效地使用铯钟驾驭氢 钟，驾驭误差较小。从图6可以看出：驾驭后氢钟的频率稳定度体现了氢钟的短期稳定度和 铯钟的长期稳定度，其短期稳定度和长期稳定度相比氢钟和铯钟分别只有很小的恶化。
[0262] 在原有的氢钟序列上加入幅度为3X 10_21^的频漂。图7展示了当氢钟存在频漂 时，氢钟、铯钟、驾驭后氢钟的Allan偏差。对比图6和图7,可以看出频漂并不会影响驾驭 后氣钟的频率稳定度。
[0263] 为了能更清楚地展示驾驭误差（即驾驭后氢钟相对铯钟的偏差），重新生成了氢 钟和铯钟序列，每个序列的参数和原氢钟和原铯钟的参数相同，共有200000个采样点。通 过使用相同的DPLL方法，获得了驾驭后氢钟的序列。图8展示了驾驭误差。从图8中可以 看出：最大驾驭误差为30ns，标准偏差为8. 99ns。实际上，驾驭误差是由于未跟踪的氢钟和 铯钟的相位噪声引起的，并且驾驭误差的标准偏差，即抖动，可以通过控制理论精确地计算 得到。因此可以推测：在通过UTC来驾驭纸面时间这类应用中，由于纸面时间和UTC都更 稳定，包含了更少的等效SSB相位噪声分量，因此驾驭误差也会相应地更小。
[0264] 在原有的氢钟序列上加入幅度为SXKT21。 1的频漂。图9展示了当氢钟含有频漂 时的驾驭误差。对比图8和图9,可以看出：当氢钟存在频漂时，会造成一个稳态跟踪误差， 幅度大约为20ns。
[0265] 综上，本发明可以有效地用参考原子钟（或时间尺度）来驾驭被驾驭原 子钟（或 时间尺度），驾驭误差和驾驭后频率稳定度的性能都较优。被驾驭原子钟（或时间尺度）的 频漂会引起一个稳态跟踪误差，这个稳态跟踪误差的值可以精确地计算得到，并通过其它 算法来进行补偿。
[0266] 综上所述，虽然本发明已以较佳实施例（通过铯钟驾驭氢钟）揭露如上，然其并非 用以限定本发明，任何本领域普通技术人员，在不脱离本发明的精神和范围内，当可作各种 更动与润饰，因此本发明的保护范围当视权利要求书界定的范围为准。
【主权项】
1. 一种等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法，其特征在于， 采用以下步骤： S. 1推导等价于二状态变量Kalman滤波器的DPLL的传递函数； S. I. 1使用DPLL，用铯钟来驾驭氢钟；DPLL的原理结构为：获取铯钟和驾驭后氢钟的时 间偏差；用环路滤波器对该时间偏差进行滤波，滤除高频噪声；然后用滤波后的时间偏差 对氢钟的时间偏差进行调整，得到驾驭后氢钟的时间偏差，其中：Cs代表铯钟的输出信号， 符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢钟的输出信号，LF代表环路 滤波器； 在Z域中，把DPLL的环路滤波器LF的传递函数记为G(Z);于是，该DPLL的输出为：(1) 显然：分别为DPLL的系统传递函数和误差传递 函数；其中：符号Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号Steered Hm 代表被驾驭氢钟的输出信号； 通过z = #π?'τ的近似变化，其中，f代表频率，T代表采样周期；将该DPLL等价为相 应的APLL ;于是，式（2)用SSB相位噪声的形式来表示，写成：式中，符号L代表SSB相位噪声，符号G代表传递函数；下标中符号Cs代表铯钟的输出 信号，符号Hm代表氢钟的输出信号，符号Steered Hm代表被驾驭氢钟的输出信号； S. 1.2根据步骤S. 1.1的系统传递函数和误差传递函数，推导等价于二状态变量 Kalman滤波器的DPLL的系统传递函数和误差传递函数的表达式； 最终，DPLL的系统传递函数表示为：其中，Ktl为环路增益；Ks Ks 22分别为Kalman滤波器的稳态增益；T为采样间隔； 该DPLL的误差传递函数表示为：至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传递函 数；由式（22)和式（23):参数（Ksn，Ks22, T，Ktl)决定了该DPLL的性能； S. 2根据步骤S. 1得到的DPLL的系统传递函数和误差传递函数，推导系统传递函数和 误差传递函数中Ks1JP Ks 22的由过程噪声方差Q22、观测噪声方差R、采样间隔T来表达的近 似值； 当Kalman滤波器进入稳态时，得到由式（30)和式（35)，Kalman增益由Q22和R之间的比例关系决定； 令Q22= Is 2,将式（30)和式（35)代入式（22)和式（23)，得到系统传递函数和误差传 递函数，分别表示为：由式（36)和式（37)可以看出：当Q22= Is 2时，参数（RJJtl)决定了 DPLL的性能；把 式（36)和式（37)代入式（1)，得到DPLL输入和输出的关系为：S. 3根据步骤S. 2得到的系统传递函数和误差传递函数，确定参数（R，T，Ktl)，优化DPLL 输出的频率稳定度，具体采用以下步骤： 首先，铯钟和氢钟的相位噪声分别表不为：其中符号L代表SSB相位噪声，&为载波频率，h i为噪声系数，f为边带频率，i用于指 明噪声类型； 通过求解方程L&(f)-LHm(f) = 0,可以求得铯钟相位噪声曲线和氢钟相位噪声曲线交 点处的频率，记为f' ； 第二，将z = #π?'τ代入式（36)和式（37)，得到了等价APLL的系统传递函数和误差 传递函数，分别表示为：和由式（41)和（42)，求得系统传递函数与误差传递函数的幅频响应201og|H(f) I和 201og I He (f) I，并求得这两条曲线交点处的频率，记为f" ； 第三，对于给定的T，首先给Ktl设一个值（，在该（T，K ^条件下，通过改变R的值，来改 变f"的值，使f" = f' ；此时，两交点的频率相等；此时，在给定T情况下，可以求得相应的 R;于是，得到了一组使两交点的频率相等的参数，记为（R'，Kc/ );当两交点的频率相等时， DPLL最大程度地滤除了铯钟的短期噪声和氢钟的长期噪声，使DPLL的输出的相位噪声达 到了最优，所以此时DPLL输出的频率稳定度达到了最优； 第四，对已给定的T，改变Ktl的值，使用步骤三中相同的方法，能够得到多组最优参数 (R'，Kc/ )，这些参数（R'，Kc/ )都是最优参数，都保证了 DPLL输出的频率稳定度达到了最 优。2.根据权利要求1所述的等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭 方法，其特征在于，步骤S. 1. 2中等价于二状态变量Kalman滤波器的DPLL的系统传递函数 和误差传递函数的推导过程如下： 对于一个二状态变量的系统，其状态方程表示为：(3) 其中，xk和yk为两个状态变量，T为采样间隔，Uk为过程噪声； 观测方程表示为： zk= x k+wk ⑷ 其中，zk为观测量，wk为观测噪声； 状态方程和观测方程这两个方程用矩阵的形式表示为：(5) 其中：过程噪声和观测噪声的方差分别为：其中，符号E代表取数学期望； Kalman滤波器用下面5个公式进行描述：其中，I是状态变量的在k时刻的估计值，Ptk是估计误差矩阵，&?是在k-1时刻对 k时刻状态变量的单步预测值，？0_1是单步预测误差矩阵，K k是Kalman增益矩阵； 能够证明式（5)定义的系统是完全可观测的，因此Ptk^Ptk和Kk都收敛，把P k^Pk^ 和Kk的稳态值分别记为：Ps、Ps _和Ks ; 由式（6)和式（9)，当Kalman滤波器进入稳态时，有：(11) 其中，毛是状态变量X的在k时刻的估计值，λ是状态变量y的在k时刻的估计值，Ksij表示Ks矩阵中的第i行第j列的元素； 定义： vA =(-Α -·ν?->ν?*Γ) (12) 将式（12)代入式（11)，得到：(13) 式（13)在Z域中表示为：(14) 其中X，Y，V分别为H 的ζ变换； 由式（14)得到：(15) 由式（14)和式（15)，式（12)在Z域中表示为：(16) 其中，Z代表Zk的Z变换； 由式（16)，得到：由式（15)、式（17)和式（18)，得到显然，Zk是Kalman滤波器的输入，而毛是Kalman滤波器的输出；于是，式（19)表明了 稳态Kalman滤波器输入与输出之间的关系； 于是，系统传递函数表示为：由式（20)，Kalman滤波器的传递函数与二阶2类DPLL的闭环系统传递函数是相同的； 于是，Kalman滤波器等价为二阶2类DPLL ; 显然，zk是Kalman滤波器的输入，而L是Kalman滤波器的输出；于是，式（19)表明了 稳态Kalman滤波器输入与输出之间的关系； 于是，系统传递函数表示为：由式（20)，Kalman滤波器的传递函数与二阶2类DPLL的闭环系统传递函数是相同的； 于是，Kalman滤波器等价为二阶2类DPLL ; 由于Ksn<< 1，于是式（20)近似写为：为了使DPLL正常工作，引入一个延迟器厂\于是DPLL的系统传递函数表示为：(22) 其中，Ktl为环路增益； 该DPLL的误差传递函数表示为：(23) 至此，得到了等价于二状态变量Kalman滤波器的DPLL的系统传递函数和误差传递函 数。3.根据权利要求2所述的等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭 方法，其特征在于，S. 2中式（30)和式（35)的推导过程如下： 当Kalman滤波器进入稳态时，式（7)写为：式（25)-(27)中，下标ij表示对应矩阵的第i行第j列的元素； 由式（25)第3式，式（26)第2式，式（27)第3式，得到： Ks2l*Psl2 =Q22 (28 ) 由式（28)和式（26)第2式，得到：(29) 由于〃vn<^，由式（29)得到(30) 由式（25)-(27)，得到(31) 由式（31)，得到(32) 由于Ks21 · T << Ksn，由式（32)，得到：(33) 把式（33)代入式（26)，得到：(34) 由式（34)，得到：(35)
【专利摘要】一种等价于二状态变量Kalman滤波器的数字锁相环原子钟驾驭方法，通过一个等价于二状态变量Kalman滤波器的DPLL，用铯钟对氢钟进行驾驭。首先推导了该DPLL的闭环系统传递函数，证明该DPLL为二阶2类DPLL。然后推导了由过程噪声方差和观测噪声方差比值近似决定的Kalman增益。于是，由DPLL的闭环传递函数可知，该比值，以及环路增益，采样间隔这3个参数决定了该DPLL的性能。进一步，本发明给出了一个简单有效的选取参数的方法。通过等效变换，得到了等效于该DPLL的模拟锁相环(APLL)。在频域中，通过调整参数，使该等价APLL的系统传递函数和误差传递函数的交点频率等于氢钟和铯钟的SSB相位噪声曲线的交点频率，最终获取最优参数。其有效地用铯钟对氢钟进行驾驭，驾驭性能非常好。
【IPC分类】G04F5/14
【公开号】CN104898404
【申请号】CN201510353792
【发明人】龚航, 伍贻威, 黄新明, 刘文祥, 欧钢, 朱祥维, 李垣陵, 黄仰博
【申请人】中国人民解放军国防科学技术大学
【公开日】2015年9月9日
【申请日】2015年6月25日
最新回复(0)