CN105577175A - Digital phase-locked loop atomic clock control method of Kalman filter and delayer - Google Patents

Abstract

The invention provides a digital phase-locked loop atomic clock control method of a Kalman filter and a delayer. A DPLL equivalent to the Kalman filter and the delayer is used to control an atomic clock, and specifically the method comprises the steps of: firstly deriving out a closed loop system transfer function and a closed loop error transfer function of the DPLL, and giving a realization structure thereof; and then obtaining the adjusting amount of the controlled atomic clock at each time, and giving a parameter selection method enabling optical frequency stability of output signals of the DPLL. On this basis, two DPLLs are cascaded to carry out two-level control on the atomic clock. Compared with a conventional atomic clock control algorithm, parameter selection is easier, the output signals are enabled to be synchronous with first-level reference input, and the optical frequency stability is ensured.

Description

Method controlled by the digital phase-locked loop atomic clock that a kind of Kalman filter adds delayer

Technical field

The present invention relates to temporal frequency, signal transacting field, devise a kind of atomic clock specifically and control method.

Background technology

Atomic clock is controlled technology and is played an important role in punctual laboratory and satellite navigation system.The main purpose controlled atomic clock has two: one to be that to make to be controlled atomic clock synchronous with referencing atom clock time, reduces the time deviation between them as far as possible; Two is promote to be controlled the long-term stability of atomic clock.

Typical atomic clock is controlled method and is comprised two large classes: method is controlled in open loop and closed loop controls method.The core that method is controlled in open loop is prediction algorithm reasonable in design.Closed loop is controlled method and is mainly comprised: Linear-Quadratic Problem (LGQ) control method, switch (Bang-Bang) control method and digital phase-locked loop (DPLL) method.But the parameter that these closed loops control method is not easy to choose, generally all for each embody rule, after quantity of parameters being compared by a large amount of emulation experiment, one group of optimized parameter to be selected.In addition, for DPLL method, if parameter choose is improper, not only cannot ensure to control performance, control system also can be caused unstable.

Summary of the invention

For the defect that prior art exists, the object of this invention is to provide the digital phase-locked loop atomic clock that a kind of Kalman filter adds delayer and control method.

Principle of the present invention is: the DPLL being equivalent to Kalman filter being added delayer by, for controlling atomic clock.The present invention has intactly derived the closed-loop system transfer function of DPLL and closed-loop error transfer function, give its implementation structure, with each for being controlled the adjustment amount of atomic clock, and give the parameter selection method of the frequency stability optimum that DPLL is outputed signal.On this basis, use two such DPLL to cascade up to carry out secondary to atomic clock and control.Theory analysis all shows with emulation experiment: this algorithm is compared conventional atom clock and controlled algorithm, and parameter choose is easier, and output signal and first order reference input can be made to keep synchronous, and ensures that frequency stability is optimum.

Technical scheme of the present invention is:

A method controlled by the digital phase-locked loop atomic clock that Kalman filter adds delayer, it is characterized in that comprising the steps:

S.1 derivation is equivalent to the transfer function that Kalman filter adds the DPLL of delayer;

S.1.1 observational equation and the state equation of Kalman filter is first provided;

State equation is expressed as:

$\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)$

Wherein, x _kand y _kbe two state variables, T is the sampling interval, u _kfor process noise.

Observational equation is expressed as:

z _k＝x _k+w _k(2)

Wherein, z _kfor observed quantity, w _kfor observation noise.

The form of these two equation matrixes is expressed as:

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

Wherein, s _k=[x _ky _k] ^t; J _k=[0u _k] ^t; $\mathrm{\φ}=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right];$ H=[10], the variance of process noise and observation noise is respectively: R=E [w _k ²], $Q=E\[J_k\·{J_k}^T\]=\left[\begin{array}{cc}0& 0\\ 0& E\[{u_k}^2\]\end{array}\right]=\left[\begin{array}{cc}0& 0\\ 0& Q_{22}\end{array}\right].$ Wherein, Q ₂₂be u _kvariance.

The observational equation of the Kalman filter S.1.2 S.1.1 provided according to step and state equation, the relation of deriving in Z territory between Kalman filter input and output, derives on this basis and provides the transfer function being equivalent to Kalman filter and adding the DPLL of delayer;

Kalman filter can be described by 5 steps below:

${\hat{s}}_{k,k-1}=\mathrm{\φ}\·{\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)

In above-mentioned 5 equations, the implication of each symbol is that the art is known, no longer describes its implication at this.

Wherein, K _kkalman gain matrix, P _k,kevaluated error matrix, P _{k, k-1}it is predicated error matrix.

The system that formula (3) defines is completely observable, therefore P _k,k, P _{k, k-1}and K _kall restrain; P _k,k, P _{k, k-1}and K _ksteady-state value be designated as respectively: Ps, Ps ^-and Ks.

By formula (4) and formula (7), when Kalman filter enters stable state, have:

$\{\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)$

Wherein, subscript ij represents K _ksteady-state value Ks matrix in i-th row jth row element.

Definition:

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

Formula (10) is substituted into formula (9), obtains:

$\{\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)$

Formula (11) is expressed as in Z territory:

$\{\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)$

Wherein X, Y, V are respectively z-transformation.

Can be obtained by formula (12):

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

By formula (12) and formula (13), formula (10) is expressed as in Z territory:

$\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)$

Wherein, Z represents z _kz-transformation.

Obtained by formula (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)$

Definition:

$G^{\′}\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)$

By formula (13), formula (15) and formula (16), obtain

$\begin{array}{c}X=\frac{G^{\′}\left(z\right)}{1-{\mathrm{Ks}}_{11}+G^{\′}\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)$

Obviously, z _kas observed quantity, it is the input of Kalman filter; And as the estimated value of state variable, be the output of Kalman filter.So formula (17) gives the relation between the constrained input of steady-state Kalman filter in Z territory.

In order to make DPLL normally work, add a delayer z in the loop ^-1.Like this, the open cycle system transfer function of DPLL is expressed as:

$G\left(z\right)=\frac{z^{-1}}{1-{\mathrm{Ks}}_{11}}\·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}---\left(18\right)$

Closed-loop system transfer function is expressed as:

$\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)$

Closed-loop error transfer function is expressed as:

$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).$

Wherein, Ks ₁₁and Ks ₂₁be respectively the steady-state gain of Kalman filter, T is the sampling interval.

Formula (18), formula (19) and formula (20) intactly give this respectively and are equivalent to Kalman filter and add the open cycle system transfer function of the DPLL of delayer, closed-loop system transfer function and closed-loop error transfer function.

As can be seen from formula (19) and formula (20), when the sampling interval, T determined, steady state Kalman gain Ks ₁₁and Ks ₂₁determine the performance of DPLL completely.Further, steady state Kalman gain Ks ₁₁and Ks ₂₁completely by process-noise variance Q ₂₂determine with observation noise variance R.In practical application, can fixation procedure noise variance Q ₂₂=1s ²constant, such observation noise variance R just determines the performance of DPLL completely.

The advantage that this DPLL compares common second order 2 class DPLL maximum is: the parameter of common DPLL has multiple, chooses rational parameter and is not easy, and also needs the stability problem of consideration system in addition; DPLL of the present invention is equivalent to Kalman filter and adds delayer, and parameter only has 1, i.e. observation noise variance R, and therefore parameter choose is relatively easy; In addition, the system defined due to formula (3) is completely observable, so Kalman filter is stable.Because this DPLL is that Kalman filter adds retarder arrangement, after numerical simulation shows to add delayer, DPLL is also stability.

The transfer function of the DPLL S.2 S.1 provided according to step, provides the implementation structure of DPLL and each adjustment amount;

Integrating step S.1 in the open cycle system transfer function of DPLL that obtains, the implementation structure of DPLL as shown in Figure 1 can be obtained; Wherein, Cs represents caesium clock, and Hm represents hydrogen clock, and SteeredHm represents the hydrogen clock after controlling.

Obtained by Fig. 1:

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

Wherein, Cs represents the output signal of caesium clock, and symbol Hm represents the output signal of hydrogen clock, and the output signal of hydrogen clock is controlled in symbol Hm_steered representative.

Obtained by formula (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)$

Formula (18) is substituted into formula (21), and the adjustment amount obtained at every turn for hydrogen clock in Z territory is $({\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).$

Controlling error, i.e. caesium clock and the deviation Cs-Hm that controls between rear hydrogen clock _steeredbe designated as Err.So in the time domain, the relation controlling the time difference of rear hydrogen clock and the time difference of free oscillation hydrogen clock is expressed as:

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

Wherein, Err (j) is for jth moment caesium clock and by time difference Cs (j)-Hm_steered (j) of controlling between hydrogen clock.

On the right of equal sign, Section 2 is each adjustment amount for hydrogen clock in the time domain; Adjustment for hydrogen clock normally utilizes phase microstepper to realize.The signal of hydrogen clock after the signal that hydrogen clock exports after phase microstepper is and controls.

The transfer function of the DPLL S.3 S.1 obtained according to step, determines that the frequency stability that parameter optimization DPLL exports, concrete grammar are fixation procedure noise variance Q ₂₂=1s ²constant, adjustment observation noise variance R.

S.3.1, the phase noise of caesium clock and hydrogen clock is expressed as:

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

With

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

Wherein f ₀for carrier frequency, h _ifor noise factor, f is sideband frequency, and i is used for indicating noise type.

By solving equation L _cs(f)-L _hmf ()=0, in the hope of the frequency at caesium clock phase noise curve and hydrogen clock phase noise intersections of complex curve place, can be designated as f '.

S.3.2, approximate change z=e is used ^{j2 π fT}, substitute into formula (19) and formula (20), obtain:

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

With

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

Formula (26) and formula (27) are closed-loop system transfer function and the closed-loop error transfer function of approximate APLL.Obviously, the closed-loop system transfer function shown in formula (26) is equivalent to a low pass filter, and the closed-loop error transfer function shown in formula (27) is equivalent to a high pass filter, so their amplitude-frequency responses intersect at a point.The frequency of intersection point is designated as f ".

S.3.3, fixation procedure noise variance Q ₂₂=1s ²constant, the value of adjustment observation noise variance R, runs Kalman filter, obtains steady state Kalman gain Ks ₁₁and Ks ₂₁value.The intersection frequency f of APLL closed-loop system transfer function and closed-loop error transfer function shown in observation type (26) and formula (27) ".As f "=f ' time, corresponding R is designated as R ', is optimized parameter.

S.3.4, R ' and corresponding steady state Kalman gain Ks is got ₁₁and Ks ₂₁, by formula (28), calculated the output power single sideband phase noise controlling rear hydrogen clock; By formula (23), obtain each adjustment amount for hydrogen clock in the time domain; Hydrogen clock is adjusted, obtains the hydrogen clock after controlling; Then the Allan deviation controlling rear hydrogen clock is calculated.

By formula (22), the output power single sideband phase noise controlling rear hydrogen clock is:

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

The parameter selection method that S.2 and S.3 the transfer function of the DPLL that S.1 the present invention obtains according to step and step obtain, obtain a kind of two-stage DPLL and control algorithm, concrete grammar is that the such DPLL of employing two cascades up, the output of first order DPLL, as the input of second level DPLL, is controlled atomic clock.The parameter selection method of each DPLL is S.2 identical with step.Two-stage controls the schematic diagram of algorithm as shown in Figure 2, wherein G ₁(z) and G ₂z () is respectively the open cycle system transfer function of two DPLL.

The present invention can be applied to design phase-locked oscillator and set up the system time of satellite navigation system (GNSS).The principle of the present invention in above-mentioned application and step are all identical.

Beneficial effect of the present invention:

1. effectively can control atomic clock, control performance (controlling error, frequency stability) very good;

2. parameter selection method of the present invention is relatively simple;

3. two-stage controls the frequency stability performance that can fully utilize two-stage input, and the short-term of output signal, mid-term, long-term stability are all optimized.

Accompanying drawing explanation

Fig. 1 graphic extension implementation structure figure of the present invention;

Fig. 2 graphic extension two-stage controls calculation ratio juris;

The time difference of Fig. 3 graphic extension hydrogen clock and caesium clock and Allan deviation, the wherein time difference of Fig. 3 a graphic extension hydrogen clock and caesium clock, the Allan deviation of Fig. 3 b graphic extension hydrogen clock and caesium clock;

The output power single sideband phase noise of Fig. 4 graphic extension hydrogen clock and caesium clock;

The amplitude-frequency response of Fig. 5 graphic extension DPLL closed-loop system transfer function and closed-loop error transfer function;

Fig. 6 graphic extension hydrogen clock, caesium clock, control the output power single sideband phase noise of rear hydrogen clock;

Fig. 7 graphic extension hydrogen clock, caesium clock, control the Allan deviation of rear hydrogen clock;

Fig. 8 graphic extension hydrogen clock, caesium clock, control time difference of rear hydrogen clock;

The time difference of Fig. 9 graphic extension NCO, hydrogen clock, caesium clock and Allan deviation, the wherein time difference of Fig. 9 a graphic extension NCO, hydrogen clock, caesium clock, the Allan deviation of Fig. 9 b graphic extension NCO, hydrogen clock, caesium clock;

The output power single sideband phase noise of Figure 10 graphic extension NCO, hydrogen clock, caesium clock;

Figure 11 graphic extension caesium clock, hydrogen clock, NCO, and the time difference and the Allan deviation of controlling rear NCO, wherein Figure 11 a graphic extension caesium clock, hydrogen clock, NCO, and control the time difference of rear NCO, Figure 11 b schemes to explain orally caesium clock, hydrogen clock, NCO, and controls the Allan deviation of rear NCO;

Embodiment

Below in conjunction with drawings and Examples, the present invention is described further.

The present invention includes following steps.

S.1 derivation is equivalent to the transfer function that Kalman filter adds the DPLL of delayer;

S.1.1 observational equation and the state equation of Kalman filter is first provided;

For the system of a two-state variable, its state equation is expressed as:

$\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)$

Wherein, x _kand y _kbe two state variables, T is the sampling interval, u _kfor process noise.

Observational equation is expressed as:

z _k＝x _k+w _k(2)

Wherein, z _kfor observed quantity, w _kfor observation noise.

The form of these two equation matrixes is expressed as:

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

Wherein, s _k=[x _ky _k] ^t; J _k=[0u _k] ^t; $\mathrm{\φ}=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right];$ H=[10], the variance of process noise and observation noise is respectively: R=E [w _k ²], $Q=E\[J_k\·{J_k}^T\]=\left[\begin{array}{cc}0& 0\\ 0& E\[{u_k}^2\]\end{array}\right]=\left[\begin{array}{cc}0& 0\\ 0& Q_{22}\end{array}\right].$ Wherein, Q ₂₂be u _kvariance.

S.1.2 in conjunction with observational equation and the state equation of Kalman filter, the relation of deriving in Z territory between Kalman filter input and output, provides the transfer function being equivalent to Kalman filter and adding the DPLL of delayer on this basis.

Kalman filter can be described by 5 steps below:

${\hat{s}}_{k,k-1}=\mathrm{\φ}\·{\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)

In above-mentioned 5 equations, the implication of each symbol is that the art is known, no longer describes its implication at this.Wherein, K _kkalman gain matrix, P _k,kevaluated error matrix, P _{k, k-1}it is predicated error matrix.

Can prove that the system that formula (3) defines is completely observable, therefore P _k,k, P _{k, k-1}and K _kall restrain. P _k,k, P _{k, k-1}and K _ksteady-state value be designated as respectively: Ps, Ps _-and Ks.

By formula (4) and formula (7), when Kalman filter enters stable state, have:

$\{\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)$

Wherein, subscript ij represents K _ksteady-state value Ks matrix in i-th row jth row element.

Definition:

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

Formula (10) is substituted into formula (9), obtains:

$\{\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)$

Formula (11) is expressed as in Z territory:

$\{\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)$

Wherein X, Y, V are respectively z-transformation.

Can be obtained by formula (12):

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

By formula (12) and formula (13), formula (10) is expressed as in Z territory:

$\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)$

Wherein, Z represents z _kz-transformation.

Obtained by formula (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)$

Definition:

$G^{\′}\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)$

By formula (13), formula (15) and formula (16), obtain

$\begin{array}{c}X=\frac{G^{\′}\left(z\right)}{1-{\mathrm{Ks}}_{11}+G^{\′}\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)$

Obviously, z _kas observed quantity, it is the input of Kalman filter; And as the estimated value of state variable, be the output of Kalman filter.So formula (17) gives the relation between the constrained input of steady-state Kalman filter in Z territory.

In order to make DPLL normally work, add a delayer z in the loop ^-1.By formula (16) and formula (17), the open cycle system transfer function of this DPLL is expressed as:

$G\left(z\right)=\frac{z^{-1}}{1-{\mathrm{Ks}}_{11}}\·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}---\left(18\right)$

Closed-loop system transfer function is expressed as:

$\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)$

Closed-loop error transfer function is expressed as:

$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).$

Formula (18), formula (19) and formula (20) intactly give this respectively and are equivalent to Kalman filter and add the open cycle system transfer function of the DPLL of delayer, closed-loop system transfer function and closed-loop error transfer function.

As can be seen from formula (19) and formula (20), when the sampling interval, T determined, steady state Kalman gain Ks ₁₁and Ks ₂₁determine the performance of DPLL completely.Further, steady state Kalman gain Ks ₁₁and Ks ₂₁completely by process-noise variance Q ₂₂determine with observation noise variance R.In practical application, can fixation procedure noise variance Q ₂₂=1s ²constant, such observation noise variance R just determines the performance of DPLL completely.

The advantage that this DPLL compares common second order 2 class DPLL maximum is: the parameter of common DPLL has multiple, chooses rational parameter and is not easy, and also needs the stability problem of consideration system in addition; DPLL is herein equivalent to Kalman filter and adds delayer, and parameter only has 1, i.e. observation noise variance R, and therefore parameter choose is relatively easy; In addition, the system defined due to formula (3) is completely observable, so Kalman filter is stable.Because this DPLL is that Kalman filter adds retarder arrangement, after numerical simulation shows to add delayer, DPLL is also stability.

As can be seen from formula (19) and formula (20), when the sampling interval, T determined, steady state Kalman gain Ks ₁₁and Ks ₂₁determine the performance of DPLL completely.Further, steady state Kalman gain Ks ₁₁and Ks ₂₁completely by process-noise variance Q ₂₂determine with observation noise variance R.In practical application, can fixation procedure noise variance Q ₂₂=1s ²constant, such observation noise variance R just determines the performance of DPLL completely.

The transfer function of the DPLL S.2 S.1 provided according to step, provides the implementation structure of DPLL and each adjustment amount;

By the open cycle system transfer function of DPLL, the implementation structure of DPLL can be obtained.Fig. 1 is to control the implementation structure figure that hydrogen clock describes DPLL with caesium clock.In Fig. 1, Cs represents caesium clock, and Hm represents hydrogen clock, and SteeredHm represents the hydrogen clock after controlling.

The implementation structure of DPLL is described below in conjunction with Fig. 1.

Obtained by Fig. 1:

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

Obtained by formula (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)$

Formula (18) is substituted into formula (21), and the adjustment amount obtained at every turn for hydrogen clock in Z territory is $({\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).$

Controlling error, i.e. caesium clock and the deviation Cs-Hm that controls between rear hydrogen clock _steeredbe designated as Err.So in the time domain, the relation controlling the time difference of rear hydrogen clock and the time difference of free oscillation hydrogen clock is expressed as:

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

Wherein, Cs represents the output signal of caesium clock, symbol Hm represents the output signal of hydrogen clock, and the output signal of hydrogen clock is controlled in symbol Hm_steered representative, and Err (j) is for jth moment caesium clock and by time difference Cs (j)-Hm_steered (j) of controlling between hydrogen clock.

On the right of equal sign, Section 2 is each adjustment amount for hydrogen clock in the time domain.Adjustment for hydrogen clock normally utilizes phase microstepper to realize.The signal of hydrogen clock after the signal that hydrogen clock exports after phase microstepper is and controls.

The transfer function of the DPLL S.3 S.1 obtained according to step, determines that the frequency stability that parameter optimization DPLL exports, concrete grammar are fixation procedure noise variance Q ₂₂=1s ²constant, adjustment observation noise variance R.

First, the phase noise of caesium clock and hydrogen clock is expressed as:

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

With

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

Wherein f ₀for carrier frequency, h _ifor noise factor, f is sideband frequency, and i is used for indicating noise type.

By solving equation L _cs(f)-L _hmf ()=0, in the hope of the frequency at caesium clock phase noise curve and hydrogen clock phase noise intersections of complex curve place, can be designated as f '.

The second, use approximate change z=e ^{j2 π fT}, substitute into formula (19) and formula (20), obtain:

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

With

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

Formula (26) and formula (27) are closed-loop system transfer function and the closed-loop error transfer function of approximate APLL.Obviously, the closed-loop system transfer function shown in formula (26) is equivalent to a low pass filter, and the closed-loop error transfer function shown in formula (27) is equivalent to a high pass filter, so their amplitude-frequency responses intersect at a point.The frequency of intersection point is designated as f ".

The frequency stability outputed signal to make DPLL is optimum, and method is herein adjustment R, makes f "=f '.In emulation experiment part, will show that the method makes the output output power single sideband phase noise of the APLL be similar to optimum.

Finally, by formula (22), the output power single sideband phase noise controlling rear hydrogen clock is:

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

Be described as embodiment the performance that this DPLL controls algorithm below to control hydrogen clock with caesium clock.

1) according to document (KasdinN.J., DiscreteSimulationofColoredNoiseandStochasticProcessesan d1/fPowerLawNoiseGeneration [J], 1995, ProceedingsoftheIEEE, 83,5, pp:802-827) method generate a hydrogen clock and a caesium clock.The parameter of hydrogen clock is: h _{(Hm) 0}=1 × 10 ^-24, h _(Hm)-2=8 × 10 ^-31.The parameter of caesium clock is: h _{(Cs) 0}=5 × 10 ^-23, h _(Cs)-2=6 × 10 ^-32.Sampling interval T=1s.Every table clock all contains 200000 data points.Fig. 3 describes the time difference of emulation hydrogen clock and emulation caesium clock, and Allan deviation.If the carrier frequency of emulation hydrogen clock and caesium clock is f ₀=10MHz.

2) by formula (24) and formula (25), the output power single sideband phase noise curve of hydrogen clock and caesium clock is made, as shown in Figure 4.As can be seen from Figure 4, the frequency f of these two intersections of complex curve ' be approximately 10 ^-3.9hz.

3) fixation procedure noise variance Q ₂₂=1s ²constant, the value of adjustment observation noise variance R, runs Kalman filter, obtains steady state Kalman gain Ks ₁₁and Ks ₂₁value.The intersection frequency f of APLL closed-loop system transfer function and closed-loop error transfer function shown in observation type (26) and formula (27) ".As f "=f ' time, corresponding R is designated as R ', is optimized parameter.Through experiment, find to work as R '=2 × 10 ^-14time, the frequency f of the point of intersection of the amplitude-frequency response of APLL closed-loop system transfer function and closed-loop error transfer function " be approximately 10 ^-3.9hz, is approximately equal to f ', as shown in Figure 5.So, get R '=2 × 10 ^-14.

4) R '=2 × 10 are got ^-14with the steady state Kalman gain Ks of correspondence ₁₁and Ks ₂₁, by formula (28), calculated the output power single sideband phase noise controlling rear hydrogen clock; By formula (23), obtain each adjustment amount for hydrogen clock in the time domain; Hydrogen clock is adjusted, obtains the hydrogen clock after controlling; Then the Allan deviation controlling rear hydrogen clock is calculated.Fig. 6 and Fig. 7 respectively describes hydrogen clock, caesium clock and controls output power single sideband phase noise and the Allan deviation of rear hydrogen clock.Fig. 8 describes hydrogen clock, caesium clock and controls the time difference of rear hydrogen clock.

Fig. 6 and Fig. 7 shows: the stability controlling rear hydrogen clock combines the frequency stability and the long-term frequency stability of caesium clock in mid-term of hydrogen clock; Thus it is optimum to demonstrate the frequency stability that the method can make DPLL output signal.Fig. 8 shows: control rear hydrogen clock and keep synchronous with caesium clock, demonstrate the validity controlling algorithm.

The parameter selection method that S.2 and S.3 the transfer function of the DPLL S.1 obtained according to step and step obtain, obtains a kind of two-stage DPLL and controls algorithm, specifically adopt following steps:

Adopt two such DPLL to cascade up, the output of first order DPLL, as the input of second level DPLL, is controlled atomic clock.G ₁(z) and G ₂z () is respectively the open cycle system transfer function of two DPLL.G ₁(z) and G ₂z the parameter of () is different, their parameter selection method all needs the method adopting step S.3 to describe.

Two-stage controls the schematic diagram of algorithm as shown in Figure 2; Specifically the phase-locked oscillator of two-way reference input (wherein a road is for hydrogen clock, and a road is caesium clock) is had to be that embodiment is described two-stage and controls calculation ratio juris with one: the inside phase comparator of phase-locked oscillator can obtain hydrogen clock, caesium clock and digital controlled oscillator (NCO) deviation between any two; This phase-locked oscillator uses two DPLL to control digital controlled oscillator (NCO); First DPLL caesium clock controls hydrogen clock, and form a paper time, the i.e. output of first DPLL, this controls as the first order; Then second DPLL controls NCO with this paper time, controls as the second level; Analyze theoretically, the output signal of this phase-locked oscillator by the short-term frequency stability of comprehensive NCO, the frequency stability in mid-term of hydrogen clock, the long-term frequency stability of caesium clock, and synchronous with the caesium clock retention time.Compare traditional locks oscillator phase and only make use of a road reference input (caesium clock or hydrogen clock), the method can the frequency stability of simultaneously comprehensive hydrogen clock and caesium clock, has clear superiority.

Below to be verified that by emulation experiment two-stage DPLL controls the performance of algorithm.Experimental procedure is as follows:

1) according to document (KasdinN.J., DiscreteSimulationofColoredNoiseandStochasticProcessesan d1/fPowerLawNoiseGeneration [J], 1995, ProceedingsoftheIEEE, 83,5, pp:802-827) method generate a NCO.The parameter of NCO is: h _{(NCO) 0}=2 × 10 ^-25, h _(NCO)-2=5 × 10 ^-30.Sampling interval T=1s.This NCO is altogether containing 200000 data points.The identical hydrogen clock that hydrogen clock and caesium clock adopt Section 2 herein to generate and caesium clock.Fig. 9 describes emulation NCO, hydrogen clock, time difference of caesium clock and Allan deviation.Figure 10 describes the output power single sideband phase noise of emulation NCO, hydrogen clock, caesium clock.

2) two-stage is controlled in algorithm, the method that the parameter choose of each DPLL adopts step S.3 to describe.Shown in Section 2, for DPLL1, work as R '=2 × 10 ^-14time, be similar to and have f "=f '.So for DPLL1, get R '=2 × 10 ^-14.Can be seen by Figure 10, the intersection frequency of the output power single sideband phase noise curve of hydrogen clock and NCO is about 10 ^-2.64hz.By experiment, find to work as R '=1.6 × 10 ^-9time, the closed-loop system transfer function of DPLL2 and the intersection frequency of closed-loop error transfer function are approximately f "=10 ^-2.64hz, is namely similar to and has f "=f '.So, for DPLL2, get R '=1.6 × 10 ^-9.So each DPLL1 can be calculated for hydrogen clock and the DPLL2 adjustment amount for NCO by formula (22).

3) run the DPLL of these two cascades as shown in Figure 2, by the adjustment amount calculated, respectively hydrogen clock and NCO are adjusted, finally obtain the time difference and the Allan deviation of the NCO after controlling, as shown in figure 11.

Found out by Figure 11, this embodiment demonstrates the conclusion of theory analysis, that is: the output signal of phase-locked oscillator combines the short-term frequency stability of NCO, frequency stability a middle or short term of hydrogen clock, the secular frequency stability of caesium clock, and input (caesium clock) retention time with the first order synchronous.

Control algorithm according to above-described two-stage, describe this algorithm and setting up the application in the GNSS system time.

GNSS is actually a clock synchronization system, all master stations, outer station, and the atomic clock on satellite all will keep synchronous with system time.So it is necessary for setting up a GNSS system time.In addition, the GNSS system time needs again to keep synchronous with international coordination universal time UTC.

For Beidou satellite navigation system, the system time BDT of the current Big Dipper is produced by master station master clock, and the local physics realization (being designated as: UTC (CMTC)) of the international coordination universal time (being designated as UTC) maintained with military time center frequency (being designated as CMTC) is controlled BDT, makes BDT and UTC (CMTC) retention time synchronous.

The mode of this generation BDT also has the space of improving.In future, calculated BDT will be the same with GPST, be a paper time; BDT is by the atomic clock on comprehensive master station, outer station, satellite, and service time, scaling algorithm, calculated BDT on paper and the every table clock deviation relative to paper time BDT; The frequency stability of such BDT compares the BDT produced by separate unit master clock at present with reliability, will be greatly improved.In addition, along with the foundation of satellite two-way pumping station link, the local physics realization (being designated as: UTC (NTSC)) of the international coordination universal time (being designated as UTC) that can maintain with international time service center (NTSC) further and UTC (CMTC), mutually as backup, control BDT.

Secondary in this paper is controlled algorithm and can be directly applied to and set up BDT.Secondary is controlled in algorithm, and first order DPLL controls BDT for using UTC (CMTC) or UTC (NTSC), make BDT and UTC (CMTC) or UTC (NTSC) retention time synchronous; Second level DPLL controls master station master clock for using the BDT after controlling, and obtains the physics realization of BDT, is designated as BDT (MC).

Be the local physics realization of UTC because UTC (CMTC) or UTC (NTSC) can find out, can think that " locking " is in UTC, so the same with UTC, there are higher long-term frequency stability; BDT combines the atomic clock on the ground such as multiple stage hydrogen clock, caesium clock, rubidium clock and star, therefore has higher secular frequency stability; Master station master clock is an active hydrogen clock, has higher frequency stability a middle or short term.To sum up, control algorithm at application secondary and set up in the process of Big Dipper time reference, UTC (CMTC) or UTC (NTSC) is equivalent to the caesium clock in Fig. 2, and BDT is equivalent to the hydrogen clock in Fig. 2, and BDT (MC) is equivalent to the NCO in Fig. 2.Final utilization secondary controls BDT and BDT (MC) that algorithm obtains all will be synchronous with UTC (CMTC) or UTC (NTSC) retention time; BDT is by the secular frequency stability of comprehensive BDT self, and the long-term frequency stability of UTC (CMTC) or UTC (NTSC); And BDT (MC) is by frequency stability a middle or short term of comprehensive master station master clock, the secular frequency stability of BDT self, and the long-term frequency stability of UTC (CMTC) or UTC (NTSC).

In sum; although the present invention (controls hydrogen clock by caesium clock with preferred embodiment; NCO is controlled by caesium clock, hydrogen clock two-stage) disclose as above; so itself and be not used to limit the present invention; any those of ordinary skill in the art; without departing from the spirit and scope of the present invention, when doing various change and retouching, the scope that therefore protection scope of the present invention ought define depending on claims is as the criterion.

Claims (2)

1. a method controlled by the digital phase-locked loop atomic clock that Kalman filter adds delayer, it is characterized in that: comprise the steps:

S.1 the open cycle system transfer function of DPLL in Z territory that Kalman filter adds delayer is equivalent to, closed-loop system transfer function and closed-loop error transfer function;

Open cycle system transfer function is expressed as:

$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}---\left(1\right)$

Closed-loop system transfer function is expressed as:

$\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(2\right)$

Closed-loop error transfer function is expressed as:

$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(3\right)$

Wherein, Ks ₁₁and Ks ₂₁be respectively the steady-state gain of Kalman filter, T is the sampling interval;

The transfer function of the DPLL S.2 S.1 provided according to step, provides the implementation structure of DPLL and each adjustment amount;

Adopt caesium clock to control hydrogen clock, integrating step S.1 in the open cycle system transfer function of DPLL, obtain the implementation structure of DPLL;

Hm _steered(z)＝G(z)·(Cs(z)-Hm _steered(z))+Hm(z)(4)

Wherein, Cs represents the output signal of caesium clock, and symbol Hm represents the output signal of hydrogen clock, and the output signal of hydrogen clock is controlled in symbol Hm_steered representative;

Obtained by formula (4):

$\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(5\right)$

Formula (1) is substituted into formula (4), and the adjustment amount obtained at every turn for hydrogen clock in Z territory is $({\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);$

Controlling error, i.e. caesium clock and the deviation Cs-Hm that controls between rear hydrogen clock _steeredbe designated as Err; So in the time domain, the relation controlling the time difference of rear hydrogen clock and the time difference of free oscillation hydrogen clock is expressed as:

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

Wherein, Err (j) is for jth moment caesium clock and by time difference Cs (j)-Hm_steered (j) of controlling between hydrogen clock;

On the right of equal sign, Section 2 is each adjustment amount for hydrogen clock in the time domain;

The transfer function of the DPLL S.3 S.1 obtained according to step, determines that the frequency stability that parameter optimization DPLL exports, concrete grammar are fixation procedure noise variance Q ₂₂=1s ²constant, adjustment observation noise variance R;

S.3.1, the phase noise of caesium clock and hydrogen clock is expressed as:

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

With

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

Wherein f ₀for carrier frequency, h _ifor noise factor, f is sideband frequency, and i is used for indicating noise type;

By solving equation L _cs(f)-L _hmf ()=0, in the hope of the frequency at caesium clock phase noise curve and hydrogen clock phase noise intersections of complex curve place, can be designated as f ';

S.3.2, approximate change z=e is used ^{j2 π fT}, substitute into formula (2) and formula (3), obtain:

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

With

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

Closed-loop system transfer function shown in formula (9) is equivalent to a low pass filter, and the closed-loop error transfer function shown in formula (10) is equivalent to a high pass filter, so their amplitude-frequency responses intersect at a point; The frequency of intersection point is designated as f ";

S.3.3, fixation procedure noise variance Q ₂₂=1s ²constant, the value of adjustment observation noise variance R, runs Kalman filter, obtains steady state Kalman gain Ks ₁₁and Ks ₂₁value; The intersection frequency f of APLL closed-loop system transfer function and closed-loop error transfer function shown in observation type (9) and formula (10) "; As f "=f ' time, corresponding R is designated as R ', is optimized parameter;

S.3.4, R ' and corresponding steady state Kalman gain Ks is got ₁₁and Ks ₂₁, by formula (11), calculated the output power single sideband phase noise controlling rear hydrogen clock; By formula (6), obtain each adjustment amount for hydrogen clock in the time domain; Hydrogen clock is adjusted, obtains the hydrogen clock after controlling; Then the Allan deviation controlling rear hydrogen clock is calculated;

By formula (5), the output power single sideband phase noise controlling rear hydrogen clock is:

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

2. method controlled by the digital phase-locked loop atomic clock that Kalman filter according to claim 1 adds delayer, it is characterized in that: step S.1 in, be equivalent to the open cycle system transfer function of DPLL in Z territory that Kalman filter adds delayer, closed-loop system transfer function and closed-loop error transfer function, the derivation of these three transfer functions is:

S.1.1 observational equation and the state equation of Kalman filter is first provided;

State equation is expressed as:

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

Wherein, x _kand y _kbe two state variables, T is the sampling interval, u _kfor process noise;

Observational equation is expressed as:

z _k＝x _k+w _k(13)

Wherein, z _kfor observed quantity, w _kfor observation noise;

The form of these two equation matrixes is expressed as:

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

Wherein, s _k=[x _ky _k] ^t; J _k=[0u _k] ^t; $\mathrm{\φ}=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right];$ H=[10], the variance of process noise and observation noise is respectively: R=E [w _k ²], $Q=E\[J_k\·{J_k}^T\]=\left[\begin{array}{cc}0& 0\\ 0& E\[{u_k}^2\]\end{array}\right]=\left[\begin{array}{cc}0& 0\\ 0& Q_{22}\end{array}\right];$ Wherein, Q ₂₂be u _kvariance;

The observational equation of the Kalman filter S.1.2 S.1.1 provided according to step and state equation, the relation of deriving in Z territory between Kalman filter input and output, derives on this basis and provides the transfer function being equivalent to Kalman filter and adding the DPLL of delayer;

Kalman filter can be described by 5 steps below:

${\hat{s}}_{k,k-1}=\mathrm{\φ}\·{\hat{s}}_{k-1,k-1}---\left(15\right)$

$P_{k,k-1}={\mathrm{\φP}}_{k-1,k-1}{\mathrm{\φ}}^T+Q---\left(16\right)$

$K_k=P_{k,k-1}\·H^T{(H\·P_{k,k-1}\·H^T+R)}^{-1}---\left(17\right)$

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

P _k,k＝(I-K _k·H)·P _k,k-1(19)

Wherein, K _kkalman gain matrix, P _k,kevaluated error matrix, P _{k, k-1}it is predicated error matrix;

The system that formula (14) defines is completely observable, therefore P _k,k, P _{k, k-1}and K _kall restrain; P _k,k, P _{k, k-1}and K _ksteady-state value be designated as respectively: Ps, Ps ^-and Ks;

By formula (15) and formula (18), when Kalman filter enters stable state, have:

$\left\{\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}}_{11}\·(z_k-{\hat{x}}_{k-1}-{\hat{y}}_{k-1}\·T)\end{array}\right.---\left(20\right)$

Wherein, subscript ij represents K _ksteady-state value Ks matrix in i-th row jth row element;

Definition:

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

Formula (21) is substituted into formula (20), obtains:

$\left\{\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}}_{11}\·v_k\end{array}\right.---\left(22\right)$

Formula (22) is expressed as in Z territory:

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

Wherein X, Y, V are respectively v _kz-transformation.

Can be obtained by formula (23):

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

By formula (23) and formula (24), formula (21) is expressed as in Z territory:

$\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(25\right)$

Wherein, Z represents z _kz-transformation;

Obtained by formula (25):

$(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(26\right)$

Definition:

$G^{\′}\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(27\right)$

By formula (24), formula (26) and formula (27), obtain

$\begin{array}{c}X=\frac{G^{\′}\left(z\right)}{1-{\mathrm{Ks}}_{11}+G^{\′}\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}}\end{array}.---\left(\mathrm{28}\right)$

Obviously, z _kas observed quantity, it is the input of Kalman filter; And as the estimated value of state variable, be the output of Kalman filter; So formula (28) gives the relation between the constrained input of steady-state Kalman filter in Z territory;

A delayer z is added in the open cycle system of DPLL ^-1; The open cycle system transfer function of this DPLL is expressed as:

$G\left(z\right)=\frac{z^{-1}}{1-{\mathrm{Ks}}_{11}}\·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}$

Closed-loop system transfer function is expressed as:

$\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},$

Closed-loop error transfer function is expressed as:

$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}}.$