一种Kalman滤波器加延迟器的数字锁相环原子钟驾驭方法

技术领域

本发明涉及时间频率、信号处理领域，具体的说是设计了一种原子钟驾驭方法。

背景技术

原子钟驾驭技术在守时实验室和卫星导航系统中发挥着重要作用。对原子钟进行 驾驭的主要目的有两个：一是使被驾驭原子钟与参考原子钟时间同步，尽可能缩小它们之 间的时间偏差；二是提升被驾驭原子钟的长期稳定度。

典型的原子钟驾驭方法包括两大类：开环驾驭方法和闭环驾驭方法。开环驾驭方 法的核心是设计合理的预测算法。闭环驾驭方法主要包括：线性二次型(LGQ)控制方法、开 关(Bang-Bang)控制方法和数字锁相环(DPLL)方法。但是，这些闭环驾驭方法的参数是不容 易选取的，一般都要针对每一个具体应用，通过大量仿真实验对大量参数进行比较后，选出 一组最优参数。另外，对于DPLL方法，如果参数选取不当，不但无法保证驾驭性能，还会引起 控制系统不稳定。

发明内容

针对现有技术存在的缺陷，本发明的目的是提供一种Kalman滤波器加延迟器的数 字锁相环原子钟驾驭方法。

本发明的原理是：通过一个等价于Kalman滤波器加延迟器的DPLL，用于对原子钟 进行驾驭。本发明完整地推导了DPLL的闭环系统传递函数和闭环误差传递函数，给出了其 实现结构，和每次的对于被驾驭原子钟的调整量，并给出了使DPLL输出信号的频率稳定度 最优的参数选取方法。在此基础上，使用两个这样的DPLL级联起来对原子钟进行二级驾驭。 理论分析和仿真实验都表明：该算法相比传统原子钟驾驭算法，参数选取更容易，可以使输 出信号与第一级参考输入保持同步，并保证频率稳定度最优。

本发明的技术方案是：

一种Kalman滤波器加延迟器的数字锁相环原子钟驾驭方法，其特征在于包括下述 步骤：

S.1推导等价于Kalman滤波器加延迟器的DPLL的传递函数；

S.1.1首先给出Kalman滤波器的观测方程和状态方程；

状态方程表示为：

$\left(\begin{array}{c}{x}_{k+1}=x_k+y_k·T\\ y_{k+1}=y_k+u_k\end{array}\right)---\left(1\right)$

其中，x_k和y_k为两个状态变量，T为采样间隔，u_k为过程噪声。

观测方程表示为：

z_k＝x_k+w_k(2)

其中，z_k为观测量，w_k为观测噪声。

这两个方程用矩阵的形式表示为：

$\{\begin{array}{c}{s}_{k+1}=\phi ·s_k+J_k\\ z_k=H·s_k+w_k\end{array}---\left(3\right)$

其中，s_k＝[x_ky_k]^T；J_k＝[0u_k]^T；$\phi =\left(\begin{array}{cc}1& T\\ 0& 1\end{array}\right);$H＝[10]，过程噪声和观测噪声的方差 分别为：R＝E[w_k²]，$Q=E[J_k·{J_k}^T]=\left(\begin{array}{cc}0& 0\\ 0& E\left[{u_k}^2\right]\end{array}\right)=\left(\begin{array}{cc}0& 0\\ 0& Q_{22}\end{array}\right).$其中，Q₂₂即为u_k的方差。

S.1.2根据步骤S.1.1给出的Kalman滤波器的观测方程和状态方程，推导在Z域中 Kalman滤波器输入和输出之间的关系，在此基础上推导给出等价于Kalman滤波器加延迟器 的DPLL的传递函数；

Kalman滤波器可以用下面5个步骤进行描述：

${\hat{s}}_{k,k-1}=\phi ·{\hat{s}}_{k-1,k-1}---\left(4\right)$

P_k,k-1＝φP_k-1,k-1φ^T+Q(5)

K_k＝P_k,k-1·H^T(H·P_k,k-1·H^T+R)^-1(6)

${\hat{s}}_{k,k}={\hat{s}}_{k,k-1}+K_k·(z_k-H·{\hat{s}}_{k,k-1})---\left(7\right)$

P_k,k＝(I-K_k·H)·P_k,k-1(8)

上述5个方程中各符号的含义是本技术领域所公知的，在此不再描述其含义。

其中，K_k是Kalman增益矩阵，P_k,k是估计误差矩阵，P_k,k-1是预测误差矩阵。

式(3)定义的系统是完全可观测的，因此P_k,k，P_k,k-1和K_k都收敛；把P_k,k、P_k,k-1和K_k的 稳态值分别记为：Ps、Ps^-和Ks。

由式(4)和式(7)，当Kalman滤波器进入稳态时，有：

$\{\begin{array}{c}{\hat{x}}_k={\hat{x}}_{k-1}+{\hat{y}}_{k-1}·T+{\mathrm{Ks}}_{11}·(z_k-{\hat{x}}_{k-1}-{\hat{y}}_{k-1}·T)\\ {\hat{y}}_k={\hat{y}}_{k-1}+{\mathrm{Ks}}_{21}·(z_k-{\hat{x}}_{k-1}-{\hat{y}}_{k-1}·T)\end{array}---\left(9\right)$

其中，下标ij表示K_k的稳态值Ks矩阵中的第i行第j列的元素。

定义：

$v_k=(z_k-{\hat{x}}_{k-1}-{\hat{y}}_{k-1}·T)---\left(10\right)$

将式(10)代入式(9)，得到：

$\{\begin{array}{c}{\hat{x}}_k={\hat{x}}_{k-1}+{\hat{y}}_{k-1}·T+{\mathrm{Ks}}_{11}·v_k\\ {\hat{y}}_k={\hat{y}}_{k-1}+{\mathrm{Ks}}_{21}·v_k\end{array}---\left(11\right)$

式(11)在Z域中表示为：

$\{\begin{array}{c}Z=z^{-1}·X+z^{-1}·T·Y+{\mathrm{Ks}}_{11}·V\\ Y=z^{-1}·Y+{\mathrm{Ks}}_{21}·V\end{array}---\left(12\right)$

其中X,Y,V分别为的Z变换.

由式(12)可以得到：

$X=(\frac{{\mathrm{Ks}}_{11}}{1-z^{-1}}+\frac{{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2})·V---\left(13\right)$

由式(12)和式(13)，式(10)在Z域中表示为：

$\begin{array}{c}V=Z-z^{-1}·X-T·z^{-1}·Y=Z-X+{\mathrm{Ks}}_{11}·V\\ =Z-(\frac{{\mathrm{ks}}_{11}}{1-z^{-1}}+\frac{{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2})·V+{\mathrm{Ks}}_{11}·V\end{array}---\left(14\right)$

其中，Z代表z_k的Z变换。

由式(14)得到：

$(1-{\mathrm{Ks}}_{11}+\frac{{\mathrm{Ks}}_{11}}{1-z^{-1}}+\frac{{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2})V=Z---\left(15\right)$

定义：

$G^{\prime }\left(z\right)=\frac{{\mathrm{Ks}}_{11}}{1-z^{-1}}+\frac{{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2}=\frac{{\mathrm{Ks}}_{11}·(1-z^{-1})+{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2}.---\left(16\right)$

由式(13)、式(15)和式(16)，得到

$\begin{array}{c}X=\frac{G^{\prime }\left(z\right)}{1-{\mathrm{Ks}}_{11}+G^{\prime }\left(z\right)}·Z\\ =\frac{{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-1}}{{(1-z^{-1})}^2+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-1}}·Z\end{array}.---\left(17\right)$

显然，z_k作为观测量，是Kalman滤波器的输入；而作为状态变量的估计值，是 Kalman滤波器的输出。于是，式(17)给出了在Z域中稳态Kalman滤波器的输入与输出之间的 关系。

为了使DPLL正常工作，在环路中加入一个延迟器z^-1。这样，DPLL的开环系统传递函 数表示为：

$G\left(z\right)=\frac{z^{-1}}{1-{\mathrm{Ks}}_{11}}·G^{\prime }\left(z\right)=\frac{{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·z^{-1}(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-2}}{{(1-z^{-1})}^2}---\left(18\right)$

闭环系统传递函数表示为：

$\begin{array}{c}H\left(z\right)=\frac{G\left(z\right)}{1+G\left(z\right)}\\ =\frac{{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·z^{-1}·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-2}}{{(1-z^{-1})}^2+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·z^{-1}·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-2}}\end{array},---\left(19\right)$

闭环误差传递函数表示为：

$He\left(z\right)=\frac{1}{1+G\left(z\right)}=\frac{{(1-z^{-1})}^2}{{(1-z^{-1})}^2+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·z^{-1}·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-2}}---\left(20\right).$

其中，Ks₁₁和Ks₂₁分别为Kalman滤波器的稳态增益，T为采样间隔。

式(18)，式(19)和式(20)分别完整地给出了该等价于Kalman滤波器加延迟器的 DPLL的开环系统传递函数、闭环系统传递函数和闭环误差传递函数。

由式(19)和式(20)可以看出，在采样间隔T确定的情况下，稳态Kalman增益Ks₁₁和 Ks₂₁完全决定了DPLL的性能。进一步地，稳态Kalman增益Ks₁₁和Ks₂₁完全由过程噪声方差Q₂₂和观测噪声方差R决定。实际应用中，可以固定过程噪声方差Q₂₂＝1s²不变，这样观测噪声方 差R就完全决定了DPLL的性能。

该DPLL相比普通的二阶2类DPLL最大的优势在于：普通DPLL的参数有多个，选取合 理的参数并不容易，另外还需要考虑系统的稳定性问题；本发明的DPLL等价于Kalman滤波 器加延迟器，参数只有1个，即观测噪声方差R，因此参数选取相对容易；另外，由于式(3)定 义的系统是完全可观测的，所以Kalman滤波器是稳定的。由于该DPLL是Kalman滤波器加延 迟器结构，数值仿真表明加入延迟器后DPLL也是稳定度。

S.2根据步骤S.1给出的DPLL的传递函数，给出DPLL的实现结构和每次的调整量；

结合步骤S.1中得到的DPLL的开环系统传递函数，可以得到如图1所示的DPLL的实 现结构；其中，Cs代表铯钟，Hm代表氢钟，SteeredHm代表驾驭后的氢钟。

由图1得到：

Hm_steered(z)＝G(z)·(Cs(z)-Hm_steered(z))+Hm(z)(21)

其中，Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号Hm_steered代表 被驾驭氢钟的输出信号。

由式(21)得到：

$\begin{array}{c}{\mathrm{Hm}}_{steered}\left(z\right)=\frac{G\left(z\right)}{1+G\left(z\right)}·Cs\left(z\right)+\frac{1}{1+G\left(z\right)}·Hm\left(z\right)\\ =H\left(z\right)·Cs\left(z\right)+He\left(z\right)·Hm\left(z\right)\end{array}---\left(22\right)$

将式(18)代入式(21)，得到在Z域中每次对于氢钟的调整量为 $({\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·\frac{z^{-1}}{1-z^{-1}})·\left(Cs\right(z)-{\mathrm{Hm}}_{steered}(z\left)\right).$

把驾驭误差，即铯钟和驾驭后氢钟之间的偏差Cs-Hm_steered记为Err。于是，在时域 中，驾驭后氢钟的时差与自由振荡氢钟的时差的关系表示为：

$Hm\_steered(i+1)=Hm\left(i\right)+\underset{j=1}{\overset{i}{\Sigma }}({\mathrm{ks}}_{11}·Err(j)+\underset{k=1}{\overset{j-1}{\Sigma }}{\mathrm{Ks}}_{21}·T·Err(k\left)\right)/(1-{\mathrm{Ks}}_{11})---\left(23\right)$

其中，Err(j)为第j时刻铯钟与被驾驭氢钟之间的时差Cs(j)-Hm_steered(j)。

等号右边第二项即为在时域中每次对于氢钟的调整量；对于氢钟的调整通常是利 用相位微跃计来实现的。氢钟经过相位微跃计后输出的信号即为驾驭后氢钟的信号。

S.3根据步骤S.1得到的DPLL的传递函数，确定参数优化DPLL输出的频率稳定度， 具体方法是固定过程噪声方差Q₂₂＝1s²不变，调整观测噪声方差R。

S.3.1，铯钟和氢钟的相位噪声分别表示为：

$L_{Cs}\left(f\right)=10log(0.5·\frac{{f_0}^2}{f^2}·\underset{-2}{\overset{2}{\Sigma }}{h}_{\left(Cs\right)i}·{f_{\left(Cs\right)}}^i),---\left(24\right)$

和

$L_{Hm}\left(f\right)=10log(0.5·\frac{{f_0}^2}{f^2}·\underset{-2}{\overset{2}{\Sigma }}{h}_{\left(Hm\right)i}·{f_{\left(Hm\right)}}^i),---\left(25\right)$

其中f₀为载波频率，h_i为噪声系数，f为边带频率，i用于指明噪声类型。

通过求解方程L_Cs(f)-L_Hm(f)＝0，可以求得铯钟相位噪声曲线和氢钟相位噪声曲 线交点处的频率，记为f’。

S.3.2，使用近似变化z＝e^j2πf·T，代入式(19)和式(20)，得到：

$\begin{array}{c}H\left(f\right)\\ =\frac{{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}·(1-e^{-j2\pi f·T})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}}{(1-e^{-j2\pi f·T})+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}·(1-e^{-j2\pi f·T})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}}\end{array}---\left(26\right)$

和

$\begin{array}{c}He\left(f\right)\\ =\frac{{(1-e^{-j2\pi f·T})}^2}{{(1-e^{-j2\pi f·T})}^2+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}·(1-e^{-j2\pi f·T})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}}\end{array}---\left(27\right)$

式(26)和式(27)为近似的APLL的闭环系统传递函数和闭环误差传递函数。显然， 式(26)所示的闭环系统传递函数相当于一个低通滤波器，而式(27)所示的闭环误差传递函 数相当于一个高通滤波器，所以它们幅频响应曲线相交于一点。把交点的频率记为f”。

S.3.3，固定过程噪声方差Q₂₂＝1s²不变，调整观测噪声方差R的值，运行Kalman滤 波器，得到稳态Kalman增益Ks₁₁和Ks₂₁的值。观察式(26)和式(27)所示APLL闭环系统传递函 数和闭环误差传递函数的交点频率f”。当f”＝f’时，对应的R记为R’，为最优参数。

S.3.4，取R’和对应的稳态Kalman增益Ks₁₁和Ks₂₁，由式(28)，计算得到了驾驭后氢 钟的单边带相位噪声；由式(23)，得到在时域中每次对于氢钟的调整量；对氢钟进行调整， 得到了驾驭后的氢钟；然后计算得到驾驭后氢钟的Allan偏差。

由式(22)，驾驭后氢钟的单边带相位噪声为：

$\begin{array}{c}{L}_{Hm\_steered}\left(f\right)=|\frac{G\left(e^{j2\pi f·T}\right)}{1+G\left(e^{j2\pi f·T}\right)}{|}^2·L_{Cs}\left(f\right)+|\frac{1}{1+G\left(e^{j2\pi f·T}\right)}{|}^2·L_{Hm}\left(f\right)\\ =|H\left(e^{j2\pi f·T}\right){|}^2·L_{Cs}\left(f\right)+|H\left(e^{j2\pi f·T}\right){|}^2·L_{Hm}\left(f\right)\end{array}---\left(28\right)$

本发明根据步骤S.1得到的DPLL的传递函数以及步骤S.2和S.3得到的参数选取方 法，得到了一种两级DPLL驾驭算法，具体方法是采用两个这样的DPLL级联起来，第一级DPLL 的输出作为第二级DPLL的输入，对原子钟进行驾驭。每个DPLL的参数选取方法和步骤S.2相 同。两级驾驭算法的原理图如图2所示，其中G₁(z)和G₂(z)分别为两个DPLL的开环系统传递 函数。

本发明可以应用于设计锁相振荡器和建立卫星导航系统(GNSS)的系统时间。本发 明在上述应用中的原理和步骤都是相同的。

本发明的有益效果：

1.可以有效地对原子钟进行驾驭，驾驭性能(驾驭误差、频率稳定度)非常好；

2.本发明的参数选取方法相对比较简单；

3.两级驾驭可以综合利用两级输入的频率稳定度性能，使输出信号的短期、中期、 长期稳定度都得到优化。

附图说明

图1图解说明本发明的实现结构图；

图2图解说明两级驾驭算法的原理；

图3图解说明氢钟和铯钟的时差以及Allan偏差，其中图3a图解说明氢钟和铯钟的 时差，图3b图解说明氢钟和铯钟的Allan偏差；

图4图解说明氢钟和铯钟的单边带相位噪声；

图5图解说明DPLL闭环系统传递函数和闭环误差传递函数的幅频响应；

图6图解说明氢钟、铯钟、驾驭后氢钟的单边带相位噪声；

图7图解说明氢钟、铯钟、驾驭后氢钟的Allan偏差；

图8图解说明氢钟、铯钟、驾驭后氢钟的时差；

图9图解说明NCO、氢钟、铯钟的时差和Allan偏差，其中图9a图解说明NCO、氢钟、铯 钟的时差，图9b图解说明NCO、氢钟、铯钟的Allan偏差；

图10图解说明NCO、氢钟、铯钟的单边带相位噪声；

图11图解说明铯钟、氢钟、NCO，和驾驭后NCO的时差和Allan偏差，其中图11a图解 说明铯钟、氢钟、NCO，和驾驭后NCO的时差，图11b图解说铯钟、氢钟、NCO，和驾驭后NCO的 Allan偏差；

具体实施方式

下面将结合附图和实施例对本发明做进一步的说明。

本发明包括以下步骤。

S.1推导等价于Kalman滤波器加延迟器的DPLL的传递函数；

S.1.1首先给出Kalman滤波器的观测方程和状态方程；

对于一个二状态变量的系统，其状态方程表示为：

$\left(\begin{array}{c}{x}_{k+1}=x_k+y_k·T\\ y_{k+1}=y_k+u_k\end{array}\right)---\left(1\right)$

其中，x_k和y_k为两个状态变量，T为采样间隔，u_k为过程噪声。

观测方程表示为：

z_k＝x_k+w_k(2)

其中，z_k为观测量，w_k为观测噪声。

这两个方程用矩阵的形式表示为：

$\{\begin{array}{c}{s}_{k+1}=\phi ·s_k+J_k\\ z_k=H·s_k+w_k\end{array}---\left(3\right)$

其中，s_k＝[x_ky_k]^T；J_k＝[0u_k]^T；$\phi =\left(\begin{array}{cc}1& T\\ 0& 1\end{array}\right);$H＝[10]，过程噪声和观测噪声的方差 分别为：R＝E[w_k²]，$Q=E[J_k·{J_k}^T]=\left(\begin{array}{cc}0& 0\\ 0& E\left[{u_k}^2\right]\end{array}\right)=\left(\begin{array}{cc}0& 0\\ 0& Q_{22}\end{array}\right).$其中，Q₂₂即为u_k的方差。

S.1.2结合Kalman滤波器的观测方程和状态方程，推导在Z域中Kalman滤波器输入 和输出之间的关系，在此基础上给出等价于Kalman滤波器加延迟器的DPLL的传递函数。

Kalman滤波器可以用下面5个步骤进行描述：

${\hat{s}}_{k,k-1}=\phi ·{\hat{s}}_{k-1,k-1}---\left(4\right)$

P_k,k-1＝φP_k-1,k-1φ^T+Q(5)

K_k＝P_k,k-1·H^T(H·P_k,k-1·H^T+R)^-1(6)

${\hat{s}}_{k,k}={\hat{s}}_{k,k-1}+K_k·(z_k-H·{\hat{s}}_{k,k-1})---\left(7\right)$

P_k,k＝(I-K_k·H)·P_k,k-1(8)

上述5个方程中各符号的含义是本技术领域所公知的，在此不再描述其含义。其 中，K_k是Kalman增益矩阵，P_k,k是估计误差矩阵，P_k,k-1是预测误差矩阵。

可以证明式(3)定义的系统是完全可观测的，因此P_k,k，P_k,k-1和K_k都收敛.把P_k,k、 P_k,k-1和K_k的稳态值分别记为：Ps、Ps_-和Ks。

由式(4)和式(7)，当Kalman滤波器进入稳态时，有：

$\{\begin{array}{c}{\hat{x}}_k={\hat{x}}_{k-1}+{\hat{y}}_{k-1}·T+{\mathrm{Ks}}_{11}·(z_k-{\hat{x}}_{k-1}-{\hat{y}}_{k-1}·T)\\ {\hat{y}}_k={\hat{y}}_{k-1}+{\mathrm{Ks}}_{21}·(z_k-{\hat{x}}_{k-1}-{\hat{y}}_{k-1}·T)\end{array}---\left(9\right)$

其中，下标ij表示K_k的稳态值Ks矩阵中的第i行第j列的元素。

定义：

$v_k=(z_k-{\hat{x}}_{k-1}-{\hat{y}}_{k-1}·T)---\left(10\right)$

将式(10)代入式(9)，得到：

$\{\begin{array}{c}{\hat{x}}_k={\hat{x}}_{k-1}+{\hat{y}}_{k-1}·T+{\mathrm{Ks}}_{11}·v_k\\ {\hat{y}}_k={\hat{y}}_{k-1}+{\mathrm{Ks}}_{21}·v_k\end{array}---\left(11\right)$

式(11)在Z域中表示为：

$\{\begin{array}{c}Z=z^{-1}·X+z^{-1}·T·Y+{\mathrm{Ks}}_{11}·V\\ Y=z^{-1}·Y+{\mathrm{Ks}}_{21}·V\end{array}---\left(12\right)$

其中X,Y,V分别为的Z变换.

由式(12)可以得到：

$X=(\frac{{\mathrm{Ks}}_{11}}{1-z^{-1}}+\frac{{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2})·V---\left(13\right)$

由式(12)和式(13)，式(10)在Z域中表示为：

$\begin{array}{c}V=Z-z^{-1}·X-T·z^{-1}·Y=Z-X+{\mathrm{Ks}}_{11}·V\\ =Z-(\frac{{\mathrm{ks}}_{11}}{1-z^{-1}}+\frac{{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2})·V+{\mathrm{Ks}}_{11}·V\end{array},---\left(14\right)$

其中，Z代表z_k的Z变换。

由式(14)得到：

$(1-{\mathrm{Ks}}_{11}+\frac{{\mathrm{Ks}}_{11}}{1-z^{-1}}+\frac{{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2})V=Z---\left(15\right)$

定义：

$G^{\prime }\left(z\right)=\frac{{\mathrm{Ks}}_{11}}{1-z^{-1}}+\frac{{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2}=\frac{{\mathrm{Ks}}_{11}·(1-z^{-1})+{\mathrm{Ks}}_{21}·T·z^{-1}}{{(1-z^{-1})}^2}.---\left(16\right)$

由式(13)、式(15)和式(16)，得到

$\begin{array}{c}X=\frac{G^{\prime }\left(z\right)}{1-{\mathrm{Ks}}_{11}+G^{\prime }\left(z\right)}·Z\\ =\frac{{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-1}}{{(1-z^{-1})}^2+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-1}}·Z\end{array}.---\left(17\right)$

显然，z_k作为观测量，是Kalman滤波器的输入；而作为状态变量的估计值，是 Kalman滤波器的输出。于是，式(17)给出了在Z域中稳态Kalman滤波器的输入与输出之间的 关系。

为了使DPLL正常工作，在环路中加入一个延迟器z^-1。由式(16)和式(17)，该DPLL的 开环系统传递函数表示为：

$G\left(z\right)=\frac{z^{-1}}{1-{\mathrm{Ks}}_{11}}·G^{\prime }\left(z\right)=\frac{{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·z^{-1}(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-2}}{{(1-z^{-1})}^2}---\left(18\right)$

闭环系统传递函数表示为：

$\begin{array}{c}H\left(z\right)=\frac{G\left(z\right)}{1+G\left(z\right)}\\ =\frac{{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·z^{-1}·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-2}}{{(1-z^{-1})}^2+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·z^{-1}·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-2}}\end{array},---\left(19\right)$

闭环误差传递函数表示为：

$He\left(z\right)=\frac{1}{1+G\left(z\right)}=\frac{{(1-z^{-1})}^2}{{(1-z^{-1})}^2+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·z^{-1}·(1-z^{-1})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·z^{-2}}---\left(20\right).$

式(18)，式(19)和式(20)分别完整地给出了该等价于Kalman滤波器加延迟器的 DPLL的开环系统传递函数、闭环系统传递函数和闭环误差传递函数。

由式(19)和式(20)可以看出，在采样间隔T确定的情况下，稳态Kalman增益Ks₁₁和 Ks₂₁完全决定了DPLL的性能。进一步地，稳态Kalman增益Ks₁₁和Ks₂₁完全由过程噪声方差Q₂₂和观测噪声方差R决定。实际应用中，可以固定过程噪声方差Q₂₂＝1s²不变，这样观测噪声方 差R就完全决定了DPLL的性能。

该DPLL相比普通的二阶2类DPLL最大的优势在于：普通DPLL的参数有多个，选取合 理的参数并不容易，另外还需要考虑系统的稳定性问题；本文的DPLL等价于Kalman滤波器 加延迟器，参数只有1个，即观测噪声方差R，因此参数选取相对容易；另外，由于式(3)定义 的系统是完全可观测的，所以Kalman滤波器是稳定的。由于该DPLL是Kalman滤波器加延迟 器结构，数值仿真表明加入延迟器后DPLL也是稳定度。

由式(19)和式(20)可以看出，在采样间隔T确定的情况下，稳态Kalman增益Ks₁₁和 Ks₂₁完全决定了DPLL的性能。进一步地，稳态Kalman增益Ks₁₁和Ks₂₁完全由过程噪声方差Q₂₂和观测噪声方差R决定。实际应用中，可以固定过程噪声方差Q₂₂＝1s²不变，这样观测噪声方 差R就完全决定了DPLL的性能。

S.2根据步骤S.1给出的DPLL的传递函数，给出DPLL的实现结构和每次的调整量；

由DPLL的开环系统传递函数，可以得到DPLL的实现结构。图1以用铯钟驾驭氢钟为 例描述了DPLL的实现结构图。图1中，Cs代表铯钟，Hm代表氢钟，SteeredHm代表驾驭后的氢 钟。

下面结合图1描述DPLL的实现结构。

由图1得到：

Hm_steered(z)＝G(z)·(Cs(z)-Hm_steered(z))+Hm(z)(21)

由式(21)得到：

$\begin{array}{c}{\mathrm{Hm}}_{steered}\left(z\right)=\frac{G\left(z\right)}{1+G\left(z\right)}·Cs\left(z\right)+\frac{1}{1+G\left(z\right)}·Hm\left(z\right)\\ =H\left(z\right)·Cs\left(z\right)+He\left(z\right)·Hm\left(z\right)\end{array}---\left(22\right)$

将式(18)代入式(21)，得到在Z域中每次对于氢钟的调整量为$({\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·\frac{z^{-1}}{1-z^{-1}})·\left(Cs\right(z)-{\mathrm{Hm}}_{steered}(z\left)\right).$

把驾驭误差，即铯钟和驾驭后氢钟之间的偏差Cs-Hm_steered记为Err。于是，在时域 中，驾驭后氢钟的时差与自由振荡氢钟的时差的关系表示为：

$Hm\_steered(i+1)=Hm\left(i\right)+\underset{j=1}{\overset{i}{\Sigma }}({\mathrm{ks}}_{11}·Err(j)+\underset{k=1}{\overset{j-1}{\Sigma }}{\mathrm{Ks}}_{21}·T·Err(k\left)\right)/(1-{\mathrm{Ks}}_{11})---\left(23\right)$

其中，Cs代表铯钟的输出信号，符号Hm代表氢钟的输出信号，符号Hm_steered代表 被驾驭氢钟的输出信号，Err(j)为第j时刻铯钟与被驾驭氢钟之间的时差Cs(j)-Hm_ steered(j)。

等号右边第二项即为在时域中每次对于氢钟的调整量。对于氢钟的调整通常是利 用相位微跃计来实现的。氢钟经过相位微跃计后输出的信号即为驾驭后氢钟的信号。

S.3根据步骤S.1得到的DPLL的传递函数，确定参数优化DPLL输出的频率稳定度， 具体方法是固定过程噪声方差Q₂₂＝1s²不变，调整观测噪声方差R。

首先，铯钟和氢钟的相位噪声分别表示为：

$L_{Cs}\left(f\right)=10log(0.5·\frac{{f_0}^2}{f^2}·\underset{-2}{\overset{2}{\Sigma }}{h}_{\left(Cs\right)i}·{f_{\left(Cs\right)}}^i),---\left(24\right)$

和

$L_{Hm}\left(f\right)=10log(0.5·\frac{{f_0}^2}{f^2}·\underset{-2}{\overset{2}{\Sigma }}{h}_{\left(Hm\right)i}·{f_{\left(Hm\right)}}^i),---\left(25\right)$

其中f₀为载波频率，h_i为噪声系数，f为边带频率，i用于指明噪声类型。

通过求解方程L_Cs(f)-L_Hm(f)＝0，可以求得铯钟相位噪声曲线和氢钟相位噪声曲 线交点处的频率，记为f’。

第二，使用近似变化z＝e^j2πf·T，代入式(19)和式(20)，得到：

$\begin{array}{c}H\left(f\right)\\ =\frac{{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}·(1-e^{-j2\pi f·T})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}}{(1-e^{-j2\pi f·T})+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}·(1-e^{-j2\pi f·T})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}}\end{array}---\left(26\right)$

和

$\begin{array}{c}He\left(f\right)\\ =\frac{{(1-e^{-j2\pi f·T})}^2}{{(1-e^{-j2\pi f·T})}^2+{\mathrm{Ks}}_{11}/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}·(1-e^{-j2\pi f·T})+{\mathrm{Ks}}_{21}·T/(1-{\mathrm{Ks}}_{11})·e^{-j2\pi f·T}}\end{array}---\left(27\right)$

式(26)和式(27)为近似的APLL的闭环系统传递函数和闭环误差传递函数。显然， 式(26)所示的闭环系统传递函数相当于一个低通滤波器，而式(27)所示的闭环误差传递函 数相当于一个高通滤波器，所以它们幅频响应曲线相交于一点。把交点的频率记为f”。

为了使DPLL输出信号的频率稳定度最优，本文的方法是调整R，使f”＝f’。在仿真 实验部分，将展示该方法使近似的APLL的输出单边带相位噪声最优。

最终，由式(22)，驾驭后氢钟的单边带相位噪声为：

$\begin{array}{c}{L}_{Hm\_steered}\left(f\right)=|\frac{G\left(e^{j2\pi f·T}\right)}{1+G\left(e^{j2\pi f·T}\right)}{|}^2·L_{Cs}\left(f\right)+|\frac{1}{1+G\left(e^{j2\pi f·T}\right)}{|}^2·L_{Hm}\left(f\right)\\ =|H\left(e^{j2\pi f·T}\right){|}^2·L_{Cs}\left(f\right)+|H\left(e^{j2\pi f·T}\right){|}^2·L_{Hm}\left(f\right)\end{array}---\left(28\right)$

下面以用铯钟驾驭氢钟作为实施例进行说明该DPLL驾驭算法的性能。

1)按照文献(KasdinN.J.,DiscreteSimulationofColoredNoiseand StochasticProcessesand1/fPowerLawNoiseGeneration[J],1995,Proceedings oftheIEEE,83,5,pp:802-827)的方法生成一台氢钟和一台铯钟。氢钟的参数为：h_(Hm)0＝1 ×10^-24，h_(Hm)-2＝8×10^-31。铯钟的参数为：h_(Cs)0＝5×10^-23，h_(Cs)-2＝6×10^-32。采样间隔T＝ 1s。每台钟都含有200000个数据点。图3描述了仿真氢钟和仿真铯钟的时差，以及Allan偏 差。设仿真氢钟和铯钟的载波频率为f₀＝10MHz。

2)由式(24)和式(25)，作出氢钟和铯钟的单边带相位噪声曲线，如图4所示。从图4 中可以看出，这两条曲线交点的频率f’大约为10^-3.9Hz。

3)固定过程噪声方差Q₂₂＝1s²不变，调整观测噪声方差R的值，运行Kalman滤波器， 得到稳态Kalman增益Ks₁₁和Ks₂₁的值。观察式(26)和式(27)所示APLL闭环系统传递函数和 闭环误差传递函数的交点频率f”。当f”＝f’时，对应的R记为R’，为最优参数。经过实验，发 现当R’＝2×10^-14时，APLL闭环系统传递函数和闭环误差传递函数的幅频响应的交点处的 频率f”近似为10^-3.9Hz，近似等于f’，如图5所示。于是，取R’＝2×10^-14。

4)取R’＝2×10^-14和对应的稳态Kalman增益Ks₁₁和Ks₂₁，由式(28)，计算得到了驾 驭后氢钟的单边带相位噪声；由式(23)，得到在时域中每次对于氢钟的调整量；对氢钟进行 调整，得到了驾驭后的氢钟；然后计算得到驾驭后氢钟的Allan偏差。图6和图7分别描述了 氢钟、铯钟、和驾驭后氢钟的单边带相位噪声和Allan偏差。图8描述了氢钟、铯钟、和驾驭后 氢钟的时差。

图6和图7表明：驾驭后氢钟的稳定度综合了氢钟的中期频率稳定度和铯钟的长期 频率稳定度；从而验证了该方法可以使DPLL输出信号的频率稳定度最优。图8表明：驾驭后 氢钟与铯钟保持同步，验证了驾驭算法的有效性。

根据步骤S.1得到的DPLL的传递函数以及步骤S.2和S.3得到的参数选取方法，得 到了一种两级DPLL驾驭算法，具体采用以下步骤：

采用两个这样的DPLL级联起来，第一级DPLL的输出作为第二级DPLL的输入，对原 子钟进行驾驭。G₁(z)和G₂(z)分别为两个DPLL的开环系统传递函数。G₁(z)和G₂(z)的参数是 不一样的，它们的参数选取方法都需要采用步骤S.3描述的方法。

两级驾驭算法的原理图如图2所示；具体以一个有两路参考输入(其中一路为氢 钟，一路为铯钟)的锁相振荡器为实施例进行说明两级驾驭算法的原理：锁相振荡器的内部 比相器可以获取氢钟、铯钟、和数控振荡器(NCO)两两之间的偏差；该锁相振荡器使用两个 DPLL对数控振荡器(NCO)进行驾驭；第一个DPLL用铯钟驾驭氢钟，形成一个纸面时间，即第 一个DPLL的输出，这作为第一级驾驭；然后第二个DPLL用该纸面时间驾驭NCO，作为第二级 驾驭；从理论上分析，该锁相振荡器的输出信号将综合NCO的短期频率稳定度，氢钟的中期 频率稳定度，铯钟的长期频率稳定度，并且与铯钟保持时间同步。相比传统锁相振荡器只利 用了一路参考输入(铯钟或是氢钟)，该方法可以同时综合氢钟和铯钟的频率稳定度，具有 明显优势。

下面以通过仿真实验来验证两级DPLL驾驭算法的性能。实验步骤如下：

1)按照文献(KasdinN.J.,DiscreteSimulationofColoredNoiseand StochasticProcessesand1/fPowerLawNoiseGeneration[J],1995,Proceedings oftheIEEE,83,5,pp:802-827)的方法生成一个NCO。NCO的参数为：h_(NCO)0＝2×10^-25， h_(NCO)-2＝5×10^-30。采样间隔T＝1s。该NCO共含有200000个数据点。氢钟和铯钟采用本文第2 节生成的相同的氢钟和铯钟。图9描述了仿真NCO、氢钟、铯钟的时差和Allan偏差。图10描述 了仿真NCO、氢钟、铯钟的单边带相位噪声。

2)两级驾驭算法中，每个DPLL的参数选取采用步骤S.3描述的方法。由第2节所示， 对于DPLL1，当R’＝2×10^-14时，近似有f”＝f’。于是对于DPLL1，取R’＝2×10^-14。由图10可以 看到，氢钟与NCO的单边带相位噪声曲线的交点频率约为10^-2.64Hz。通过实验，发现当R’＝ 1.6×10^-9时，DPLL2的闭环系统传递函数和闭环误差传递函数的交点频率近似为f”＝10 ^-2.64Hz，即近似有f”＝f’。于是，对于DPLL2，取R’＝1.6×10^-9。于是可以由式(22)计算得到 每次DPLL1对于氢钟和DPLL2对于NCO的调整量。

3)运行如图2所示的这两个级联的DPLL，通过计算得到的调整量，分别对氢钟和 NCO进行调整，最终得到驾驭后的NCO的时差和Allan偏差，如图11所示。

由图11看出，该实施例验证了理论分析的结论，即：锁相振荡器的输出信号综合了 NCO的短期频率稳定度，氢钟的中短期频率稳定度，铯钟的中长期频率稳定度，并且与第一 级输入(铯钟)保持时间同步。

根据上面描述的两级驾驭算法，描述了该算法在建立GNSS系统时间中的应用。

GNSS实际上是一个时间同步系统，所有主控站，外场站，卫星上的原子钟都要和系 统时间保持同步。所以，建立一个GNSS系统时间是必须的。此外，GNSS系统时间又需要和国 际协调世界时UTC保持同步。

以北斗卫星导航系统为例，目前北斗的系统时间BDT由主控站主钟产生，并用军用 时间频率中心(记为CMTC)维持的国际协调世界时(记为UTC)的本地物理实现(记为：UTC (CMTC))对BDT进行驾驭，使BDT与UTC(CMTC)保持时间同步。

这种产生BDT的方式还有改进的空间。未来，计划中的BDT将和GPST一样，是一个纸 面时间；BDT将综合主控站、外场站、卫星上的原子钟，使用时间尺度算法，计算得到纸面上 的BDT和每台钟相对于纸面时间BDT的偏差；这样BDT的频率稳定度和可靠性相比目前由单 台主钟产生的BDT，将得到很大的提升。另外，随着卫星双向比对链路的建立，可以进一步用 国际授时中心(NTSC)维持的国际协调世界时(记为UTC)的本地物理实现(记为：UTC(NTSC)) 和UTC(CMTC)互相作为备份，对BDT进行驾驭。

本文提出的二级驾驭算法可以直接应用于建立BDT。二级驾驭算法中，第一级DPLL 用于使用UTC(CMTC)或UTC(NTSC)对BDT进行驾驭，使BDT与UTC(CMTC)或UTC(NTSC)保持时间 同步；第二级DPLL用于使用驾驭后的BDT对主控站主钟进行驾驭，得到BDT的物理实现，记为 BDT(MC)。

由于UTC(CMTC)或是UTC(NTSC)可以看出是UTC的本地物理实现，可以认为“锁定” 于UTC，所以和UTC一样，具有较高的长期频率稳定度；BDT综合了多台氢钟、铯钟、铷钟等地 面和星上的原子钟，因此具有较高的中长期频率稳定度；主控站主钟是一台主动型氢钟，具 有较高的中短期频率稳定度。综上，在应用二级驾驭算法建立北斗时间基准的过程中，UTC (CMTC)或UTC(NTSC)相当于图2中的铯钟，BDT相当于图2中的氢钟，BDT(MC)相当于图2中的 NCO。最终使用二级驾驭算法得到的BDT和BDT(MC)都将与UTC(CMTC)或UTC(NTSC)保持时间 同步；BDT将综合BDT自身的中长期频率稳定度，和UTC(CMTC)或UTC(NTSC)的长期频率稳定 度；而BDT(MC)将综合主控站主钟的中短期频率稳定度，BDT自身的中长期频率稳定度，和 UTC(CMTC)或UTC(NTSC)的长期频率稳定度。

综上所述，虽然本发明已以较佳实施例(通过铯钟驾驭氢钟，通过铯钟、氢钟两级 驾驭NCO)揭露如上，然其并非用以限定本发明，任何本领域普通技术人员，在不脱离本发明 的精神和范围内，当可作各种更动与润饰，因此本发明的保护范围当视权利要求书界定的 范围为准。