[发明专利]一种Kalman滤波器加延迟器的数字锁相环原子钟驾驭方法

权利要求书

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

S.1 等价于Kalman滤波器加延迟器的DPLL在Z域中的开环系统传递函数，闭环系统传递函数和闭环误差传递函数；

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

状态方程表示为：

x k + 1 = x k + y k · T y k + 1 = y k + u k - - - ( 1 ) ]]>

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

观测方程表示为：

z_k＝x_k+w_k (2)

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

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

s k + 1 = φ · s k + J k z k = H · s k + w k - - - ( 3 ) ]]>

其中，s_k＝[x_k y_k]^T；J_k＝[0 u_k]^T；H＝[1 0]，过程噪声和观测噪声的方差分别为：R＝E[w_k²]，其中，Q₂₂即为u_k的方差；

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

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

s ^ k , k - 1 = φ · s ^ k - 1 , k - 1 - - - ( 4 ) ]]>

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)

s ^ k , k = s ^ k , k - 1 + K k · ( z k - H · s ^ k , k - 1 ) - - - ( 7 ) ]]>

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

其中，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滤波器进入稳态时，有：

x ^ k = x ^ k - 1 + y ^ k - 1 · T + Ks 11 · ( z k - x ^ k - 1 - y ^ k - 1 · T ) y ^ k = y ^ k - 1 + Ks 21 · ( z k - x ^ k - 1 - y ^ k - 1 · T ) - - - ( 9 ) ]]>

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

定义：

v k = ( z k - x ^ k - 1 - y ^ k - 1 · T ) - - - ( 10 ) ]]>

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

x ^ k = x ^ k - 1 + y ^ k - 1 · T + Ks 11 · v k y ^ k = y ^ k - 1 + Ks 21 · v k - - - ( 11 ) ]]>

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

X = z - 1 · X + z - 1 · T · Y + Ks 11 · V Y = z - 1 · Y + Ks 21 · V - - - ( 12 ) ]]>

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

由式(12)可以得到：

X = ( Ks 11 1 - z - 1 + Ks 21 · T · z - 1 ( 1 - z - 1 ) 2 ) · V - - - ( 13 ) ]]>

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

V = Z - z - 1 · X - T · z - 1 · Y = Z - X + Ks 11 · V = Z - ( Ks 11 1 - z - 1 + Ks 21 · T · z - 1 ( 1 - z - 1 ) 2 ) · V + Ks 11 · V - - - ( 14 ) ]]>

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

由式(14)得到：

( 1 - Ks 11 + Ks 11 1 - z - 1 + Ks 21 · T · z - 1 ( 1 - z - 1 ) 2 ) V = Z - - - ( 15 ) ]]>

定义：

G ′ ( z ) = Ks 11 1 - z - 1 + Ks 21 · T · z - 1 ( 1 - z - 1 ) 2 = Ks 11 · ( 1 - z - 1 ) + Ks 21 · T · z - 1 ( 1 - z - 1 ) 2 . - - - ( 16 ) ]]>

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

X = G ′ ( z ) 1 - Ks 11 + G ′ ( z ) · Z = Ks 11 / ( 1 - Ks 11 ) · ( 1 - z - 1 ) + Ks 21 · T / ( 1 - Ks 11 ) · z - 1 ( 1 - z - 1 ) 2 + Ks 11 / ( 1 - Ks 11 ) · ( 1 - z - 1 ) + Ks 21 · T / ( 1 - Ks 11 ) · z - 1 · Z . - - - ( 17 ) ]]>

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

在DPLL的开环系统中加入一个延迟器z^-1；该DPLL的开环系统传递函数表示为：

G ( z ) = z - 1 1 - Ks 11 · G ′ ( z ) = Ks 11 / ( 1 - Ks 11 ) · z - 1 ( 1 - z - 1 ) + Ks 21 · T / ( 1 - Ks 11 ) · z - 2 ( 1 - z - 1 ) 2 - - - ( 18 ) ]]>

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

H ( z ) = G ( z ) 1 + G ( z ) = Ks 11 / ( 1 - Ks 11 ) · z - 1 · ( 1 - z - 1 ) + Ks 21 · T / ( 1 - Ks 11 ) · z - 2 ( 1 - z - 1 ) 2 + Ks 11 / ( 1 - Ks 11 ) · z - 1 · ( 1 - z - 1 ) + Ks 21 · T / ( 1 - Ks 11 ) · z - 2 , - - - ( 19 ) ]]>

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

H e ( z ) = 1 1 + G ( z ) = ( 1 - z - 1 ) 2 ( 1 - z - 1 ) 2 + Ks 11 / ( 1 - Ks 11 ) · z - 1 · ( 1 - z - 1 ) + Ks 21 · T / ( 1 - Ks 11 ) · z - 2 - - - ( 20 ) . ]]>

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

式(18)，式(19)和式(20)分别为等价于Kalman滤波器加延迟器的DPLL的开环系统传递函数、闭环系统传递函数和闭环误差传递函数；

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

采用铯钟驾驭氢钟，结合步骤S.1中的DPLL的开环系统传递函数，得到DPLL的实现结构；

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

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

由式(21)得到：

Hm s t e e r e d ( z ) = G ( z ) 1 + G ( z ) · C s ( z ) + 1 1 + G ( z ) · H m ( z ) = H ( z ) · C s ( z ) + H e ( z ) · H m ( z ) - - - ( 22 ) ]]>

将式(18)代入式(21)，得到在Z域中每次对于氢钟的调整量为

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

H m _ s t e e r e d ( i + 1 ) = H m ( i ) + Σ j = 1 i ( Ks 11 · E r r ( j ) + Σ k = 1 j - 1 Ks 21 · T · E r r ( k ) ) / ( 1 - Ks 11 ) - - - ( 23 ) ]]>

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

等号右边第二项即为在时域中每次对于氢钟的调整量；

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

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

L C s ( f ) = 10 l o g ( 0.5 · f 0 2 f 2 · Σ - 2 2 h ( C s ) i · f ( C s ) i ) , - - - ( 24 ) ]]>

和

L H m ( f ) = 10 l o g ( 0.5 · f 0 2 f 2 · Σ - 2 2 h ( H m ) i · f ( H m ) i ) , - - - ( 25 ) ]]>

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

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

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

H ( f ) = Ks 11 / ( 1 - Ks 11 ) · e - j 2 π f · T · ( 1 - e - j 2 π f · T ) + Ks 21 · T / ( 1 - Ks 11 ) · e - j 2 π f · 2 T ( 1 - e - j 2 π f · T ) 2 + Ks 11 / ( 1 - Ks 11 ) · e - j 2 π f · T . ( 1 - e - j 2 π f · T ) + Ks 21 · T / ( 1 - Ks 11 ) · e - j 2 π f · 2 T - - - ( 26 ) ]]>

和

H e ( f ) = ( 1 - e - j 2 π f · T ) 2 ( 1 - e - j 2 π f · T ) 2 + Ks 11 / ( 1 - Ks 11 ) · e - j 2 π f · T · ( 1 - e - j 2 π f · T ) + Ks 21 · T / ( 1 - Ks 11 ) · e - j 2 π f · 2 T - - - ( 27 ) ]]>

式(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)，驾驭后氢钟的单边带相位噪声为：

L H m _ s t e e r e d ( f ) = | G ( e j 2 π f · T ) 1 + G ( e j 2 π f · T ) | 2 · L C s ( f ) + | 1 1 + G ( e j 2 π f · T ) | 2 · L H m ( f ) = | H ( e j 2 π f · T ) | 2 · L C s ( f ) + | H e ( e j 2 π f · T ) | 2 · L H m ( f ) - - - ( 28 ) . ]]>