CN104848862A - Precise and synchronous positioning and time-keeping method and system of Mars orbiting detector - Google Patents

Abstract

The invention discloses a precise and synchronous positioning and time-keeping method and system of a Mars orbiting detector. The precise and synchronous positioning and time-keeping method comprises the following steps: firstly inputting initial data including a simulated nominal orbit and relative parameters for filtering initialization; simulating observation data, and simulating the observed quantity of a pulsar according to a measurement equation of the pulsar, wherein measurement noises relative to the observed quantity comprise clock errors, and the clock errors are simulated by a clock error model; and performing adaptive Kalman filtering according to a filtering state equation and an observation equation of self-positioning and time keeping of the Mars orbiting detector to obtain orbits and the clock errors of the detector. According to the precise and synchronous positioning and time-keeping method and system disclosed by the invention, on the basis of establishing an observation model of an X-ray pulsar, the observation model, a detector kinetic model and the clock error model of a satellite-borne atomic clock are organically combined, and the position, speed and clock error parameters of the detector are also estimated so as to ensure that the clock error correction can be achieved, influences of the clock errors to the positioning accuracy can also be weakened, and the self-positioning accuracy can also be improved.

Description

The punctual method and system in a kind of ring fire detector precision synchronous location

Technical field

The present invention relates to mars exploration technical field, particularly relate to the punctual method and system in a kind of ring fire detector precision synchronous location.

Background technology

For performing the Mars orbit detector of task of exploring the Mars, in real time, the relevant informations such as the position of detector, speed and attitude are obtained accurately, being not only the basis of detector being carried out to precise navigation control, is also one of deciding factor guaranteeing that various predetermined scientific experiment and detection mission are implemented smoothly.

The complicated space environment that remote Distance geometry detector between Mars and the earth faces is a huge challenge for the conventional navigation method based on ground observing and controlling operating office, and the independent navigation of spacecraft will be one of gordian technique breaking through this bottleneck.X-ray pulsar navigation is a kind of novel autonomous navigation method, in Mars probes independent navigation, have application potential.

Kinetic model and the observations of pulsar model of ring fire detector is related to when utilizing X-ray pulsar observed quantity to position ring fire detector, can to be extrapolated by numerical integration by the kinetic model of ring fire detector but only can obtain the approximate orbit information of detector, in addition, due to technical conditions restriction, also there is error in the initial rail state of detector.

In X-ray pulsar navigational system, photon arrives time of detector and is obtained by satellite atomic clock measurement, and therefore the precision of satellite atomic clock and stability thereof are directly connected to the height of navigation accuracy.And due to the Mars probes flight time long, the atomic clock that detector carries, is subject to the restriction of various factors, will be difficult to keep degree of precision for a long time if irregularly correct, and so just time to cause damage to the precision of navigational system.

Current existing X-ray pulsar navigation is a kind of novel autonomous navigation method, the feature such as possess that reliability is high, independence is strong, accuracy is good and usable range is wide, and in deep space probe independent navigation field, application potential is huge.Many scholars and expert carry out studying and verifying at this novel autonomous navigation method in the world at present, a large amount of pertinent literatures is had to introduce it both at home and abroad, Suneel I is Use of Variable Celestial X-ray Sources for Spacecraft Navigation [D] S.2005.The: [Ph.D.] .Maryland:Department of Aerospace Engineering University of Maryland.Wei Erhu, Jin Shuanggen, Zhang Qi, et al.2013Autonomous navigation of Mars probe using X-ray pulsars:Modeling and results [J] .Advances in Space Research, 51 (2013): 849 – 857.Emadzadeh A A, RobertA, Speyer J L, et al.2011.Relative Navigation between Two Spacecraft Using X-ray Pulsars [J] .IEEE Transactions on Control Systems Technology, 19 (5): 1021-1035.Yang Tinggao.2008.Determination of X-ray pulsar pulse time of arrival at spacecraft [J] .Chin.J.Space Sci., 28 (4): 330-334.Chester T J, Butman S is Using X-ray Pulsars [R] .NASATechnical Reports N81-27129 A.1981.Navigation, 21-25.

Utilize X-ray pulsar observed quantity to carry out principle that detector independently locates as shown in Figure 1.In reference point (solar system barycenter SSB) coordinate system, X-ray detection equipment return pulse signal on spacecraft Probe, and then obtain time of arrival (toa) observed quantity, then by itself and datum time phase model calculated value poor, difference τ reflects the time delay that same pulse signal arrives spacecraft and datum, and the product of τ and light velocity c is c τ.This difference observed quantity can be expressed as the function of Space Vehicle position vector r, and under the prerequisite of known navigation pulsar direction vector n, the motion dynamics equations in conjunction with spacecraft filtering can estimate the position vector of spacecraft in reference point coordinate system.

Pulsar time phase model is based upon solar system barycenter SSB place usually under normal circumstances, in J2000.0 solar system geocentric coordinate system:

A certain epoch of observation, the radiation signal of pulsar i is to the direct range observed quantity ρ of spacecraft SC and Mars barycenter M _sCi, ρ _mifor:

ρ _SCi＝||r _SC-D _i||+d _Rel,SCi(1-1)

ρ _Mi＝||r _M-D _i||+d _Rel,Mi(1-2)

In formula, r _mand r _sCrepresent the position vector of Mars barycenter and detector respectively; D _irepresent the position vector of observation pulsar; d _{rel, Mi}and d _{rel, SCi}represent corresponding pulse delay signal observational error respectively.

Formula (1-2) deducts (1-1) and obtains single poor observed quantity δ ρ _i=ρ _mi-ρ _sCi" single difference measurements equation " is:

δρ _i＝(||r _M-D _i||-||r _SC-D _i||)+(d _Rel,Mi-d _Rel,SCi) (1-3)

Because the distance D of pulsar distance reference point _imagnitude much larger than the distance of spacecraft reference point, therefore formula (1-3) is launched, and ignores dimensionless above higher order term can obtain

δρ _i＝n _i·(r _SC-r _M)+(d _Rel,Mi-d _Rel,SCi) (1-4)

Mars orbit detector is relative to the position vector r of Mars barycenter _sC/M=r _sC-r _m, the direction vector of pulsar is n _i, then the single difference observation equation comprising error term is:

$\mathrm{\δ}{\widetilde{\mathrm{\ρ}}}_i=n_i\·\mathrm{\δ}{r}_{\mathrm{SC}/M}+{\mathrm{c\δt}}_{\mathrm{SC}}+{\mathrm{\η}}_i---(1-5)$

In formula, single poor reduction $\mathrm{\δ}{\widetilde{\mathrm{\ρ}}}_i=c(t_{\mathrm{Mi}}-{\tilde{t}}_{\mathrm{SCi}})-n_i\·{\tilde{r}}_{\mathrm{SC}/M}-d_{\mathrm{Rel},\mathrm{Mi}}+d_{\mathrm{Rel},\mathrm{SCi}},$ with δ r _sC/Mrepresent that detector is relative to the position vector of Mars barycenter and correction thereof respectively; t _mirepresent that Mars pulse signal arrives the time of Mars barycenter; Desired clock reading is t _sC; with δ t _sCbe respectively time and the clock error correction amount of satellite atomic clock reading; C is the light velocity; η _irepresent random measurement noise.

And in order to reduced data process in ring fire detector is independently located, usually replace solar system barycenter with Mars barycenter.

The clock correction parameter of atomic clock is contained in measurement model (1-5), and research and analyse and show that the atomic clock clock correction of 0.5us can bring the site error of about 300m and the velocity error of 0.03m/s, and the theoretical stability index that the theoretical stability index being applied to the atomic clock of satellite and spacecraft under current technical conditions is the atomic clock being applied to satellite and spacecraft is 5 × 10 ^-12/ s, for the deep space probe that duty cycle is longer, if irregularly carry out clock correction correction, the clock correction of atomic clock will add up to tens of microseconds even more, and this causes serious harm giving the navigation accuracy of detector.

The distance of pulsar and Mars is remote, and pulse signal is faint, and the accumulation of x-ray photon generally needs the time of 8 to 10 clocks, therefore in new observation data not collected period, needs to carry out Exact Forecast to detector track.The method of conventional numerical integration carries out the Orbit extrapolation of detector, every perturbative force that the method need be considered the center gravitation that detector is subject to and act on detector.

At J2000 Mars barycenter terrestrial equator coordinate system S _eCIin, the kinetic model of detector can be expressed as:

$\stackrel{\·\·}{r}=-\frac{{\mathrm{\μ}}_{M}r}{r^3}+a_r+W---(1-6)$

In formula, r=[x y z] ^tfor detector is at S _eCIin position vector, for its second derivative, r=||r||; μ _mit is Mars gravitational constant; W is the random noise meeting Gaussian distribution; a _rthen represent the perturbation acceleration summation that various perturbative forces that detector is subject to produce.Research shows, the J for middle high orbit Mars probes in the perturbation of Mars non-spherical shape ₂item is topmost impact.If only consider J ₂item perturbation, then have:

$a_r=a_{J_2}---(1-7)$

In formula

$a_{J_2}=\left[\begin{array}{c}-\frac{{\mathrm{\μ}}_{M}x}{r^3}{J}_2{\left(\frac{R_M}{r}\right)}^2\frac{3}{2}(1-5\frac{z^2}{r^2})\\ -\frac{{\mathrm{\μ}}_{M}y}{r^3}{J}_2{\left(\frac{R_M}{r}\right)}^2\frac{3}{2}(1-5\frac{z^2}{r^2})\\ -\frac{{\mathrm{\μ}}_{M}z}{r^3}{J}_2{\left(\frac{R_M}{r}\right)}^2\frac{3}{2}(3-5\frac{z^2}{r^2})\end{array}\right]---(1-8)$

In above formula, R _mfor Mars mean radius.

If obtain the detector Orbit State Parameters in a certain moment, just can utilize the dynamics of orbits model of foundation, by the track value in the method extrapolation next moment of numerical integration.

In sum, when utilizing ring fire detector dynamics of orbits model to carry out Orbit extrapolation, the precision of extrapolation track is directly by the accuracy of dynamics of orbits model and the impact of orbit integration method, if preliminary orbit parameter exists error in addition, extrapolation track will be dispersed very soon.Can see from formula (1-8), if irregularly carry out clock correction correction, the clock correction of atomic clock will cause damage to positioning precision.

Summary of the invention

For above problem, the present invention is setting up on X-ray pulsar observation model basis, the clock bias model of observation model, detector kinetic model and satellite atomic clock is organically combined, the position of detector, speed and clock correction parameter have been estimated simultaneously, realizing also weakening the impact of clock correction on positioning precision while clock correction corrects, improve autonomous positioning precision.

Technical scheme of the present invention provides a kind of ring fire detector precision synchronous to locate punctual method, comprises the following steps,

Step 1, input primary data, comprise simulation nominal track and for the initialized correlation parameter of filtering;

Step 2, carries out observation data simulation, and measure equation simulation observations of pulsar amount according to pulsar, the measurement noises that this observed quantity relates to comprises clock correction, and wherein clock correction is simulated by clock bias model,

If navigation pulsar combination comprises p pulsar, the position vector r=[x y z] of ring fire detector ^t, velocity v=[v _xv _yv _z] ^t, acceleration a=[a _xa _ya _z] ^t,

Described pulsar measure equation as shown in the formula,

$Z=\left[\begin{array}{c}\mathrm{\Δ}{t}_1\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{t}_i\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{t}_p\end{array}\right]=H\·X+\mathrm{\η}$

Wherein, Δ t _irepresent the observed quantity of i-th pulsar time delay, the value of i is 1,2 ..., p, p observations of pulsar amount composition observation vector Z, η is measurement noises, and H is observation equation matrix of coefficients, state parameter X=[x y z v _xv _yv _zx ₁x ₂x ₃] ^t;

Described clock bias model as shown in the formula,

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+\mathrm{\ϵ}\·x_2\left(t_k\right)+\frac{{\mathrm{\ϵ}}^2}{2}\·x_3\left(t_k\right)+W\left(t_k\right)$

Wherein, ε is sampling interval, x ₁(t _k), x ₂(t _k) and x ₃(t _k) represent moment t respectively _kthe clock correction of atomic clock, frequency drift and frequency drift rate of change, constitute clock correction parameter [x ₁x ₂x ₃]; W (t _k) be moment t _kcorresponding system noise, x ₁(t _k+1) represent moment t _k+1the clock correction of atomic clock;

Step 3, independently locates according to ring fire detector and the filter state equation of keeping time and observation equation, carries out adaptive Kalman filter,

Described observation equation adopts pulsar to measure equation,

If the state parameter X of moment t is designated as X (t), the first order derivative of X (t) is described filter state equation is as follows,

$\stackrel{\·}{X}\left(t\right)=F\left(t\right)\·X\left(t\right)+W\left(t\right)$

Wherein, W (t) is system state noise matrix,

State Equation Coefficients matrix F (t) is, $F\left(t\right)=\left[\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_x}{\∂y}& \frac{\∂a_x}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_y}{\∂y}& \frac{\∂a_y}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_z}{\∂x}& \frac{\∂a_z}{\∂y}& \frac{\∂a_z}{\∂z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$

Step 4, exports step 3 gained filter result, comprises track and the clock correction of filtering gained detector.

And the mode that step 3 carries out adaptive Kalman filter is as follows,

1) the correlation parameter initialization of filtering, comprises the state parameter x calculating initial time _kand variance matrix P _k, state error δ x _k;

2) one-step prediction value x is carried out _{k+1, k}with modified value δ x _{k+1, k}be calculated as follows,

x _k+1,k＝Φ _k+1,k·x _k

δx _k+1,k＝0

Wherein, Φ _{k+1, k}for state-transition matrix;

3) adaptive factor λ is carried out _k+1be calculated as follows,

$M_{k+1}=H_{k+1}{\mathrm{\Φ}}_{k+1,k}{P}_k{\mathrm{\Φ}}_{k+1,k}^T{H}_{k+1}^T+R_{k+1}$

δz _k+1＝z _k+1-H _k+1x _k+1,k

$N_{k+1}=\mathrm{\δ}{z}_{k+1}\mathrm{\δ}{z}_{k+1}^T$

${\mathrm{\λ}}_{k+1}=\mathrm{Max}\{1,\sqrt{\mathrm{sp}\left(N_{k+1}\right)/\mathrm{sp}\left(M_{k+1}\right)}\}$

In formula, H _k+1represent observation equation matrix of coefficients; Symbol M _k+1and N _k+1for calculating adaptive factor; z _k+1for observation vector; Prediction residual is δ z _k+1; R _k+1for observed reading varivance matrix; Sp (.) is for asking trace of a matrix symbol;

4) one-step prediction estimation error variance battle array P is carried out _{k+1, k}be calculated as follows,

$P_{k+1,k}={\mathrm{\λ}}_k{\mathrm{\Φ}}_{k+1,k}{P}_k{\mathrm{\Φ}}_{k+1,k}^T+Q_k$

In formula, Q _kfor system state noise variance matrix;

5) filter gain K is carried out _k+1be calculated as follows,

$K_{k+1}=P_{k+1,k}{H}_{k+1}^T{[H_{k+1}{P}_{k+1,k}{H}_{k+1}^T+R_{k+1}]}^{-1}$

6) measurement updaue is carried out as follows,

δx _k+1＝δx _k+1,k+K _k+1[z _k+1-H _k+1x _k+1,k]

$P_{k+1}=[I-K_{k+1}{H}_{k+1}]P_{k+1,k}{[I-K_{k+1}{H}_{k+1}]}^T+K_{k+1}{R}_{k+1}{K}_{k+1}^T$

In formula, P _k+1for filter state value x _k+1variance matrix, I represents the unit matrix of corresponding exponent number;

7) filter result x is carried out _k+1export as follows,

x _k+1＝x _k+1,k+δx _k+1

And, described one-step prediction value x _{k+1, k}numerical integration calculating is carried out by filter state equation by orbit integration device.

And, in advance according to observability and the combination of selection indicators value selection navigation pulsar of pulsar.

And the nominal track of simulation obtains according to STK emulation.

And, by the track of step 3 filtering gained detector and the nominal track of simulation poor, step 3 filtering gained clock correction and step 2 are simulated gained clock correction poor, obtain accuracy assessment result.

The present invention is also corresponding provides a kind of ring fire detector precision synchronous to locate Time keeping system, comprises with lower module,

Load module, for inputting primary data, comprise simulation nominal track and for the initialized correlation parameter of filtering;

Observation data analog module, for carrying out observation data simulation, measure equation simulation observations of pulsar amount according to pulsar, the measurement noises that this observed quantity relates to comprises clock correction, and wherein clock correction is simulated by clock bias model,

If navigation pulsar combination comprises p pulsar, the position vector r=[x y z] of ring fire detector ^t, velocity v=[v _xv _yv _z] ^t, acceleration a=[a _xa _ya _z] ^t,

Described pulsar measure equation as shown in the formula,

$Z=\left[\begin{array}{c}\mathrm{\Δ}{t}_1\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{t}_i\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{t}_p\end{array}\right]=H\·X+\mathrm{\η}$

Wherein, Δ t _irepresent the observed quantity of i-th pulsar time delay, the value of i is 1,2 ..., p, p observations of pulsar amount composition observation vector Z, η is measurement noises, and H is observation equation matrix of coefficients, state parameter X=[x y z v _xv _yv _zx ₁x ₂x ₃] ^t;

Described clock bias model as shown in the formula,

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+\mathrm{\ϵ}\·x_2\left(t_k\right)+\frac{{\mathrm{\ϵ}}^2}{2}\·x_3\left(t_k\right)+W\left(t_k\right)$

Wherein, ε is sampling interval, x ₁(t _k), x ₂(t _k) and x ₃(t _k) represent moment t respectively _kthe clock correction of atomic clock, frequency drift and frequency drift rate of change, constitute clock correction parameter [x ₁x ₂x ₃]; W (t _k) be moment t _kcorresponding system noise, x ₁(t _k+1) represent moment t _k+1the clock correction of atomic clock;

Auto adapted filtering module, for the filter state equation of independently locating according to ring fire detector and keep time and observation equation, carries out adaptive Kalman filter,

Described observation equation adopts pulsar to measure equation,

If the state parameter X of moment t is designated as X (t), the first order derivative of X (t) is described filter state equation is as follows,

$\stackrel{\·}{X}\left(t\right)=F\left(t\right)\·X\left(t\right)+W\left(t\right)$

Wherein, W (t) is system state noise matrix,

State Equation Coefficients matrix F (t) is, $F\left(t\right)=\left[\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_x}{\∂y}& \frac{\∂a_x}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_y}{\∂y}& \frac{\∂a_y}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_z}{\∂x}& \frac{\∂a_z}{\∂y}& \frac{\∂a_z}{\∂z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$

Output module, for output adaptive filtration module gained filter result, comprises track and the clock correction of filtering gained detector.

The present invention considers to utilize X-ray pulsar observed quantity, in conjunction with the dynamics of orbits model of detector and the clock bias model of satellite atomic clock, by auto adapted filtering fast, estimate the position of ring fire detector, speed and clock correction parameter accurately, thus correction satellite atomic clock, weaken or eliminate clock correction to the impact of positioning precision, realize the precision synchronous location of ring fire detector and keep time.Invention also improves auto adapted filtering implementation, raise the efficiency further.The present invention completes under state natural sciences fund support, has important actual promotional value and application prospect, has very important effect to the development of national economy and the raising of living standards of the people.

Accompanying drawing explanation

Fig. 1 is the schematic diagram that prior art utilizes X-ray pulsar observed quantity to carry out detector independently to locate;

Fig. 2 is the process flow diagram of the embodiment of the present invention.

Embodiment

Below in conjunction with drawings and Examples, technical solution of the present invention is specifically described.

The present invention considers, the precision of X-ray pulsar observed quantity is relatively high, if kinetic model and measurement model are combined by filtering algorithm, the X-ray pulsar observed quantity utilizing precision higher is revised in real time to track extrapolated state value, just can estimate the Orbit State Parameters of ring fire detector accurately.

In ring fire detector is independently located, the kinetic model of Mars orbit detector can't accurately be set up, when numerical integration carries out Orbit extrapolation, the impact of some perturbation factors is left in the basket usually, directly cause extrapolation trajectory accuracy poor, and due to technical conditions restriction, preliminary orbit difficult parameters is accurately to obtain, and this also can cause dispersing of extrapolation track.How to utilize the X-ray pulsar observed quantity that precision is higher, quick, the high precision estimation that are realized detector's status parameter by filtering algorithm are one of matters of utmost importance realizing precision positioning.In addition, because the Mars probes flight time is long, if the atomic clock that it carries does not carry out regular calibration keep higher precision by being difficult to, from the measurement model of X-ray pulsar, the atomic clock clock correction of detector is very large to location precision, therefore how effectively to realize the estimation to clock correction, be keep atomic clock precision and improve that the autonomous positioning precision of detector faces have a key issue.Solution of the present invention is exactly the unified model establishing a kind of effective integration X-ray pulsar observation model, detector dynamics of orbits model and clock bias model, devise a kind of adaptive extended kalman filtering algorithm, synchronous, quick, the high precision that realize detector position, speed and clock correction parameter are resolved.Concrete scheme is as follows: (1) sets up ring fire detector accurate independently location and punctual model based on X-ray pulsar, comprises the clock bias model of measurement model, detector dynamics of orbits model and satellite atomic clock.

In conjunction with the concrete condition of ring fire detector, set up corresponding measurement model; Analyze ring fire detector stressing conditions, the basis of synthesis precision and counting yield is accepted or rejected some perturbation factors, then determines the dynamics of orbits model of detector; Understand principle of work and the performance estimating method of satellite atomic clock, on this basis modeling is effectively carried out to clock correction.

(2) adaptive extended kalman filtering algorithm is designed, the clock bias model of measurement model, dynamics of orbits model and atomic clock is merged utilization effectively, clock bias model parameter is incorporated in filtering algorithm, synchronously resolve with the position of detector and speed parameter, realize fast, hi-Fix and punctual.

(3) specific implementation

Comprise observation data simulation and adaptive extended kalman filtering resolve two large links, specific design is as follows:

During concrete enforcement, STK software emulation ring fire detector track can be utilized as " nominal track ", and suitable nautical star is selected from navigation pulsar database, according to the observed quantity of orbit measurement model (1-5) formula analog pulse star, the impacts such as clocking error, observational error, relativistic effect delay in simulation process, to be taken into account.

In conjunction with the concrete condition of ring fire detector, set up the X-ray pulsar observation equation of detector.

At J2000 Mars barycenter terrestrial equator coordinate system S _eCIin, if the same period has observed p pulsar, then ring fire detector has independently been located and with the measurement equation of Time keeping system has been:

$Z=\left[\begin{array}{c}\mathrm{\Δ}{t}_1\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{t}_i\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{t}_p\end{array}\right]=H\·X+\mathrm{\η}---(2-1)$

In formula, Δ t _irepresent the observed quantity of i-th pulsar time delay, the value of i is 1,2 ..., p, p observations of pulsar amount composition observation vector Z; η is measurement noises; State parameter X=[x y z v _xv _yv _zx ₁x ₂x ₃] ^t, observation equation matrix of coefficients H is:

$H=\left[\begin{array}{ccccccc}\frac{n_1}{c}& 0& 0& 0& 1& 0& 0\\ & & & \·& & & \\ & & & \·& & & \\ & & & \·& & & \\ \frac{n_i}{c}& 0& 0& 0& 1& 0& 0\\ & & & \·& & & \\ & & & \·& & & \\ & & & \·& & & \\ \frac{n_p}{c}& 0& 0& 0& 1& 0& 0\end{array}\right]---(2-2)$

In formula, n _iit is the direction vector of i-th pulsar; C is light velocity value.

Analyze the stressing conditions of ring fire detector, determine the perturbation factors will considered in orbit determination, set up the motion dynamics equations of ring fire detector, and utilize Taylor series to carry out linearization expansion.

At J2000 Mars barycenter terrestrial equator coordinate system S _eCIin, if the position vector r=of ring fire detector [x y z] ^t, velocity v=[v _xv _yv _z] ^t, acceleration a=[a _xa _ya _z] ^t; The state parameter vector of detector is x (t)=[x y z v _xv _yv _z] ^t, the first order derivative of x (t) is ω (t) is t state-noise, then its kinetics equation of t can be expressed as:

$\stackrel{\·}{x}\left(t\right)=f(x\left(t\right),t)+\mathrm{\ω}\left(t\right)---(2-3)$

Only considering Mars center gravitation and J ₂when item perturbation, this equation can do following expression:

$\stackrel{\·}{x}=f(x\left(t\right),t)=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\\ -{\mathrm{\μ}}_M\frac{x}{r_3}[1-\frac{3}{2}{J}_2{\left(\frac{R_M}{2}\right)}^2[\frac{5z^2}{r^2}-1\left]\right]\\ -{\mathrm{\μ}}_M\frac{y}{r^3}[1-\frac{3}{2}{J}_2{\left(\frac{R_M}{r}\right)}^2[\frac{5z^2}{r^2}-1\left]\right]\\ -{\mathrm{\μ}}_M\frac{z}{r^3}[1-\frac{3}{2}{J}_2{\left(\frac{R_M}{r}\right)}^2[\frac{5z^2}{r^2}-3\left]\right]\end{array}\right]---(2-4)$

(x (t), t) represents detector kinetics equation to f, i.e. detector model trajectory, is hereafter abbreviated as f in formula.

Formula (2-4) is nonlinear system, needs to carry out linearization at " equilibrium point " place by Taylor series expansion and obtains:

$\stackrel{\·}{x}\left(t\right)=G\left(t\right)\·x\left(t\right)+\mathrm{\ω}\left(t\right)---(2-5)$

In formula (2-5), t Jacobian matrix G (t) is abbreviated as G, its computing method:

$G=\frac{\∂f}{\∂x}=\frac{\∂}{\∂x}\left[\begin{array}{c}v\\ a\end{array}\right]=\left[\begin{array}{cc}\frac{\∂v}{\∂r}& \frac{\∂v}{\∂v}\\ \frac{\∂a}{\∂r}& \frac{\∂a}{\∂v}\end{array}\right]---(2-6)$

$\frac{\∂v}{\∂r}=\frac{\∂a}{\∂v}=0_{3\×3}---(2-7)$

$\frac{\∂v}{\∂v}=I_{3\×3}---(2-8)$

$\frac{\∂a}{\∂}=\left[\begin{array}{ccc}\frac{\∂a_x}{\∂x}& \frac{\∂a_x}{\∂y}& \frac{\∂a_x}{\∂z}\\ \frac{\∂a_y}{\∂x}& \frac{\∂a_y}{\∂y}& \frac{\∂a_y}{\∂z}\\ \frac{\∂a_z}{\∂x}& \frac{\∂a_z}{\∂y}& \frac{\∂a_z}{\∂z}\end{array}\right]---(2-9)$

In formula, 0 _{3 × 3}represent null matrix; I _{3 × 3}representation unit matrix.

Set up the clock bias model of atomic clock, and clock correction state equation is combined with kinetics equation forms filter state equation, using the position of clock bias model parameter and detector, speed parameter as filtering parameter.

Usually, the time deviation of atomic clock can represent by a second order polynomial:

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+\mathrm{\ϵ}\·x_2\left(t_k\right)+\frac{{\mathrm{\ϵ}}^2}{2}\·x_3\left(t_k\right)+W\left(t_k\right)---(2-10)$

In formula, ε is sampling interval, x ₁(t _k), x ₂(t _k) and x ₃(t _k) represent t respectively _kthe clock correction of moment atomic clock, frequency drift and frequency drift rate of change, they constitute clock correction parameter [x ₁x ₂x ₃]; W (t _k) be t _kmoment corresponding system noise, belongs to coloured noise; x ₁(t _k+1) represent t _k+1the clock correction of moment atomic clock.

Can be obtained by formula (2-10), with clock clock correction x ₁(t _k), clock frequency drift x ₂(t _k) and frequency drift rate of change x ₃(t _k) composition quantity of state ${\left[\begin{array}{ccc}{x}_1& x_2& x_3\end{array}\right]}_k^T,$ Then clock bias model is:

${\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]}_{k+1}={\mathrm{\Φ}}_{k+1,k}{\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]}_k+{\mathrm{\ω}}_k---(2-11)$

In formula, Φ _{k+1, k}for state-transition matrix, this matrix is relevant with sampling interval τ, and expression formula is:

${\mathrm{\Φ}}_{k+1,k}=\left[\begin{array}{ccc}1& \mathrm{\ϵ}& {\mathrm{\ϵ}}^2/2\\ 0.& 1& \mathrm{\ϵ}\\ 0& 0& 1\end{array}\right]---(2-12)$

In formula, ω _kfor t _kmoment clock correction noise vector.

If the state parameter X=of moment t is [x y z v _xv _yv _zx ₁x ₂x ₃] ^tbe designated as X (t), the first order derivative of X (t) is ring fire detector pulsar can be obtained in conjunction with ring fire detector motion dynamics equations (2-4) and clock bias model (2-11) independently to locate and with the state equation of Time keeping system be:

$\stackrel{\·}{X}\left(t\right)=F\left(t\right)\·X\left(t\right)+W\left(t\right)---(2-13)$

In formula, W (t) for system state noise matrix, State Equation Coefficients matrix F (t) is:

$F\left(t\right)=\left[\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_x}{\∂y}& \frac{\∂a_x}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_y}{\∂y}& \frac{\∂a_y}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_z}{\∂x}& \frac{\∂a_z}{\∂y}& \frac{\∂a_z}{\∂z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]---(2-14)$

Design adaptive extended kalman filtering device, utilizes X-ray pulsar observed quantity correction filter state parameter prediction value, resolves the position of detector, speed parameter and clock bias model parameter fast, accurately.

Establish ring fire detector by formula (2-1) to formula (2-14) independently to locate and the filter state equation of keeping time and observation equation, then need to adopt adaptive extended kalman filtering device, the filtering carrying out being correlated with calculates.The present invention proposes further, on existing adaptive extended kalman filtering device (AEKF wave filter) basis, the account form of adaptive factor and error covariance matrix can also be improved, to improve filtering speed of convergence, suppress filtering divergence, see following formula (2-18) to formula (2-22).

Systematic state transfer matrix Φ after linearization _{k+1, k}, observation equation matrix H _k+1, state equation and observation equation are respectively:

$\left\{\begin{array}{c}{x}_{k+1}={\mathrm{\Φ}}_{k+1,k}\·x_k\\ z_{k+1}=H_{k+1}\·x_{k+1}\end{array}\right.---(2-15)$

T wherein before the corresponding discretize of k, k+1 _kand t _k+1moment.

The design of embodiment adaptive extended kalman filtering device computation process is as follows:

● filtering parameter initialization

According to existing technical conditions, provide the state parameter x calculating initial time _kand variance matrix P _k, state error is δ x _k.

● one-step prediction value x _{k+1, k}with modified value δ x _{k+1, k}calculate

x _k+1,k＝Φ _k+1,k·x _k(2-16)

δx _k+1,k＝0 (2-17)

● adaptive factor λ _k+1calculating

$M_{k+1}=H_{k+1}{\mathrm{\Φ}}_{k+1,k}{P}_k{\mathrm{\Φ}}_{k+1,k}^T{H}_{k+1}^T+R_{k+1}---(2-18)$

δz _k+1＝z _k+1-H _k+1x _k+1,k(2-19)

$N_{k+1}=\mathrm{\δ}{z}_{k+1}\mathrm{\δ}{z}_{k+1}^T---(2-20)$

${\mathrm{\λ}}_{k+1}=\mathrm{Max}\{1,\sqrt{\mathrm{sp}\left(N_{k+1}\right)/\mathrm{sp}\left(M_{k+1}\right)}\}---(2-21)$

In formula, H _k+1represent observation equation matrix of coefficients; Symbol M _k+1and N _k+1for calculating adaptive factor; z _k+1for observation vector; Prediction residual is δ z _k+1; R _k+1for observed reading varivance matrix; Sp (.) is for asking trace of a matrix symbol.

● one-step prediction estimation error variance battle array P _{k+1, k}calculating

$P_{k+1,k}={\mathrm{\λ}}_k{\mathrm{\Φ}}_{k+1,k}{P}_k{\mathrm{\Φ}}_{k+1,k}^T+Q_k---(2-22)$

In formula, Q _kfor system state noise variance matrix.

● filter gain K _k+1calculating

$K_{k+1}=P_{k+1,k}{H}_{k+1}^T{[H_{k+1}{P}_{k+1,k}{H}_{k+1}^T+R_{k+1}]}^{-1}---(2-23)$

● measurement updaue

δx _k+1＝δx _k+1,k+K _k+1[z _k+1-H _k+1x _k+1,k] (2-24)

$P_{k+1}=[I-K_{k+1}{H}_{k+1}]P_{k+1,k}{[I-K_{k+1}{H}_{k+1}]}^T+K_{k+1}{R}_{k+1}{K}_{k+1}^T---(2-25)$

In formula, P _k+1for filter state value x _k+1variance matrix, I represents the unit matrix of corresponding exponent number.

● filter result x _k+1export

x _k+1＝x _k+1,k+δx _k+1(2-26)

When given detector original state and initial clock correction parameter, above-mentioned adaptive extended kalman filtering device just can be utilized to position and the calculating of parameter of keeping time.

Method provided by the present invention can adopt software mode to realize automatically running.See Fig. 2, the realization flow of embodiment is as follows: step 1, input primary data, comprise simulation nominal track and for the initialized correlation parameter of filtering.

The nominal track of simulation can obtain according to STK emulation, and when specifically implementing, the ring fire detector orbital data file that the STK that can be prepared in advance emulates, to obtain nominal track as preliminary orbit parameter for STK simulation.Candidate impulse star Parameter File can also be prepared in advance, according to navigation pulsar number, utilize existing pulsar preferable procedure to select optimal navigation pulsar to combine according to the observability of pulsar and selection indicators value.

Step 2, carries out observation data simulation, and measure equation simulation observations of pulsar amount according to pulsar, the measurement noises that this observed quantity relates to comprises clock correction, and wherein clock correction is simulated by clock bias model.

The observed quantity of equation (2-1) formula analog pulse star is measured according to pulsar, this observed quantity generally will add random observational error, clock correction, Roemer postpones and Shapiro postpones, wherein clock correction is simulated by clock bias model (2-10) formula, based on the moment t preset _kclock correction initial parameter, moment t can be simulated _k+1the clock correction of atomic clock.Other errors can be preset according to prior art.

Step 3, independently locates according to ring fire detector and the filter state equation of keeping time and observation equation, carries out adaptive Kalman filter.

Then adaptive Kalman filter and accuracy assessment module is entered, filter state equation and observation equation are respectively formula (2-13) and (2-1), concrete filtering flow process is undertaken by formula (2-16) ~ (2-26), wherein calculating (2-16) formula of one-step prediction value can carry out numerical integration calculating by orbit integration device by formula (2-13), obtains and surveys data prediction.Orbit integration device can adopt prior art, and it will not go into details in the present invention.

After filtering completes can by the track of the detector calculated with 1. in the STK nominal track of simulate poor, the clock correction that filtering calculates, with 3. the middle clock correction simulated is poor, completes accuracy assessment.

Step 4, output filtering result, comprises track and the clock correction of filtering gained detector.

During concrete enforcement, also by poor for the nominal track of the track of step 3 filtering gained detector and simulation, step 3 filtering gained clock correction and step 2 can be simulated gained clock correction poor, obtain accuracy assessment result.

During concrete enforcement, modular mode also can be adopted to provide a kind of ring fire detector precision synchronous to locate Time keeping system, and system that embodiment provides comprises with lower module:

Load module, for inputting primary data, comprise simulation nominal track and for the initialized correlation parameter of filtering;

Observation data analog module, for carrying out observation data simulation, measure equation simulation observations of pulsar amount according to pulsar, the measurement noises that this observed quantity relates to comprises clock correction, and wherein clock correction is simulated by clock bias model,

If navigation pulsar combination comprises p pulsar, the position vector r=[x y z] of ring fire detector ^t, velocity v=[v _xv _yv _z] ^t, acceleration a=[a _xa _ya _z] ^t,

Described pulsar measure equation as shown in the formula,

$Z=\left[\begin{array}{c}\mathrm{\Δ}{t}_1\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{t}_i\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{t}_p\end{array}\right]=H\·X+\mathrm{\η}$

Wherein, Δ t _irepresent the observed quantity of i-th pulsar time delay, the value of i is 1,2 ..., p, p observations of pulsar amount composition observation vector Z, η is measurement noises, and H is observation equation matrix of coefficients, state parameter X=[x y z v _xv _yv _zx ₁x ₂x ₃] ^t;

Described clock bias model as shown in the formula,

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+\mathrm{\ϵ}\·x_2\left(t_k\right)+\frac{{\mathrm{\ϵ}}^2}{2}\·x_3\left(t_k\right)+W\left(t_k\right)$

Wherein, ε is sampling interval, x ₁(t _k), x ₂(t _k) and x ₃(t _k) represent moment t respectively _kthe clock correction of atomic clock, frequency drift and frequency drift rate of change, constitute clock correction parameter [x ₁x ₂x ₃]; W (t _k) be moment t _kcorresponding system noise, x ₁(t _k+1) represent moment t _k+1the clock correction of atomic clock;

Auto adapted filtering module, for the filter state equation of independently locating according to ring fire detector and keep time and observation equation, carries out adaptive Kalman filter,

Described observation equation adopts pulsar to measure equation,

If the state parameter X of moment t is designated as X (t), the first order derivative of X (t) is described filter state equation is as follows,

$\stackrel{\·}{X}\left(t\right)=F\left(t\right)\·X\left(t\right)+W\left(t\right)$

Wherein, W (t) is system state noise matrix,

State Equation Coefficients matrix F (t) is, $F\left(t\right)=\left[\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_x}{\∂y}& \frac{\∂a_x}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_y}{\∂y}& \frac{\∂a_y}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_z}{\∂x}& \frac{\∂a_z}{\∂y}& \frac{\∂a_z}{\∂z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$

Output module, for output adaptive filtration module gained filter result, comprises track and the clock correction of filtering gained detector.

Each module specific implementation is corresponding to method step, and it will not go into details in the present invention.

Above embodiment is used for illustrative purposes only, but not limitation of the present invention, person skilled in the relevant technique, without departing from the spirit and scope of the present invention, various conversion or modification can also be made, therefore all equivalent technical schemes also should belong within category of the present invention, should be limited by each claim.

Claims (7)

1. the punctual method in ring fire detector precision synchronous location, is characterized in that: comprise the following steps,

Step 1, input primary data, comprise simulation nominal track and for the initialized correlation parameter of filtering;

Step 2, carries out observation data simulation, and measure equation simulation observations of pulsar amount according to pulsar, the measurement noises that this observed quantity relates to comprises clock correction, and wherein clock correction is simulated by clock bias model,

If navigation pulsar combination comprises p pulsar, the position vector r=[x y z] of ring fire detector ^t, velocity v=[v _xv _yv _z] ^t, acceleration a=[a _xa _ya _z] ^t,

Described pulsar measure equation as shown in the formula,

$Z=\left[\begin{array}{c}\mathrm{\Δ}{t}_1\\ .\\ .\\ .\\ \mathrm{\Δ}{t}_i\\ .\\ .\\ .\\ \mathrm{\Δ}{t}_p\end{array}\right]=H\·X+\mathrm{\η}$

Wherein, Δ t _irepresent the observed quantity of i-th pulsar time delay, the value of i is 1,2 ..., p, p observations of pulsar amount composition observation vector Z, η is measurement noises, and H is observation equation matrix of coefficients, state parameter X=[x y z v _xv _yv _zx ₁x ₂x ₃] ^t;

Described clock bias model as shown in the formula,

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+\mathrm{\ϵ}\·x_2\left(t_k\right)+\frac{{\mathrm{\ϵ}}^2}{2}\·x_3\left(t_k\right)+W\left(t_k\right)$

Wherein, ε is sampling interval, x ₁(t _k), x ₂(t _k) and x ₃(t _k) represent moment t respectively _kthe clock correction of atomic clock, frequency drift and frequency drift rate of change, constitute clock correction parameter [x ₁x ₂x ₃]; W (t _k) be moment t _kcorresponding system noise, x ₁(t _k+1) represent moment t _k+1the clock correction of atomic clock;

Step 3, independently locates according to ring fire detector and the filter state equation of keeping time and observation equation, carries out adaptive Kalman filter,

Described observation equation adopts pulsar to measure equation,

If the state parameter X of moment t is designated as X (t), the first order derivative of X (t) is described filter state equation is as follows,

$\stackrel{.}{X}\left(t\right)=F\left(t\right)\·X\left(t\right)+W\left(t\right)$

Wherein, W (t) is system state noise matrix,

State Equation Coefficients matrix F (t) is, $F\left(t\right)=\left[\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_x}{\∂y}& \frac{\∂a_x}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_y}{\∂y}& \frac{\∂a_y}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_z}{\∂x}& \frac{\∂a_z}{\∂y}& \frac{\∂a_z}{\∂z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$

Step 4, exports step 3 gained filter result, comprises track and the clock correction of filtering gained detector.

2. ring fire detector precision synchronous locates punctual method according to claim 1, it is characterized in that: the mode that step 3 carries out adaptive Kalman filter is as follows,

1) the correlation parameter initialization of filtering, comprises the state parameter x calculating initial time _kand variance matrix P _k, state error δ x _k;

2) one-step prediction value x is carried out _{k+1, k}with modified value δ x _{k+1, k}be calculated as follows,

x _k+1,k＝Φ _k+1,k·x _k

δx _k+1,k＝0

Wherein, Φ _{k+1, k}for state-transition matrix;

3) adaptive factor λ is carried out _k+1be calculated as follows,

$M_{k+1}=H_{k+1}{\mathrm{\Φ}}_{k+1,k}{P}_k{\mathrm{\Φ}}_{k+1,k}^T{H}_{k+1}^T+R_{k+1}$

δz _k+1＝z _k+1-H _k+1x _k+1,k

$N_{k+1}=\mathrm{\δ}{z}_{k+1}\mathrm{\δ}{z}_{k+1}^T$

${\mathrm{\λ}}_{k+1}=\mathrm{Max}\{1,\sqrt{\mathrm{sp}\left(N_{k+1}\right)/\mathrm{sp}\left(M_{k+1}\right)}\}$

In formula, H _k+1represent observation equation matrix of coefficients; Symbol M _k+1and N _k+1for calculating adaptive factor; z _k+1for observation vector; Prediction residual is δ z _k+1; R _k+1for observed reading varivance matrix; Sp (.) is for asking trace of a matrix symbol;

4) one-step prediction estimation error variance battle array P is carried out _{k+1, k}be calculated as follows,

$P_{k+1,k}={\mathrm{\λ}}_k{\mathrm{\Φ}}_{k+1,k}{P}_k{\mathrm{\Φ}}_{k+1,k}^T+Q_k$

In formula, Q _kfor system state noise variance matrix;

5) filter gain K is carried out _k+1be calculated as follows,

$K_{k+1}=P_{k+1,k}{H}_{k+1}^T{[H_{k+1}{P}_{k+1,k}{H}_{k+1}^T+R_{k+1}]}^{-1}$

6) measurement updaue is carried out as follows,

δx _k+1＝δx _k+1,k+K _k+1[z _k+1-H _k+1x _k+1,k]

$P_{k+1}=[I-K_{k+1}{H}_{k+1}]P_{k+1,k}{[I-K_{k+1}{H}_{k+1}]}^T+K_{k+1}{R}_{k+1}{K}_{k+1}^T$

In formula, P _k+1for filter state value x _k+1variance matrix, I represents the unit matrix of corresponding exponent number;

7) filter result x is carried out _k+1export as follows,

x _k+1＝x _k+1,k+δx _k+1。

3. ring fire detector precision synchronous locates punctual method according to claim 2, it is characterized in that: described one-step prediction value x _{k+1, k}numerical integration calculating is carried out by filter state equation by orbit integration device.

4. the punctual method in ring fire detector precision synchronous location according to claim 1 or 2 or 3, is characterized in that: in advance according to observability and the combination of selection indicators value selection navigation pulsar of pulsar.

5. the punctual method in ring fire detector precision synchronous location according to claim 1 or 2 or 3, is characterized in that: the nominal track of simulation obtains according to STK emulation.

6. the punctual method in ring fire detector precision synchronous location according to claim 1 or 2 or 3, it is characterized in that: by the track of step 3 filtering gained detector and the nominal track of simulation poor, step 3 filtering gained clock correction and step 2 are simulated gained clock correction poor, obtain accuracy assessment result.

7. a ring fire detector precision synchronous location Time keeping system, is characterized in that: comprise with lower module,

Load module, for inputting primary data, comprise simulation nominal track and for the initialized correlation parameter of filtering;

Observation data analog module, for carrying out observation data simulation, measure equation simulation observations of pulsar amount according to pulsar, the measurement noises that this observed quantity relates to comprises clock correction, and wherein clock correction is simulated by clock bias model,

If navigation pulsar combination comprises p pulsar, the position vector r=[x y z] of ring fire detector ^t, velocity v=[v _xv _yv _z] ^t, acceleration a=[a _xa _ya _z] ^t,

Described pulsar measure equation as shown in the formula,

$Z=\left[\begin{array}{c}\mathrm{\Δ}{t}_1\\ .\\ .\\ .\\ \mathrm{\Δ}{t}_i\\ .\\ .\\ .\\ \mathrm{\Δ}{t}_p\end{array}\right]=H\·X+\mathrm{\η}$

Wherein, Δ t _irepresent the observed quantity of i-th pulsar time delay, the value of i is 1,2 ..., p, p observations of pulsar amount composition observation vector Z, η is measurement noises, and H is observation equation matrix of coefficients, state parameter X=[x y z v _xv _yv _zx ₁x ₂x ₃] ^t;

Described clock bias model as shown in the formula,

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+\mathrm{\ϵ}\·x_2\left(t_k\right)+\frac{{\mathrm{\ϵ}}^2}{2}\·x_3\left(t_k\right)+W\left(t_k\right)$

Wherein, ε is sampling interval, x ₁(t _k), x ₂(t _k) and x ₃(t _k) represent moment t respectively _kthe clock correction of atomic clock, frequency drift and frequency drift rate of change, constitute clock correction parameter [x ₁x ₂x ₃]; W (t _k) be moment t _kcorresponding system noise, x ₁(t _k+1) represent moment t _k+1the clock correction of atomic clock;

Auto adapted filtering module, for the filter state equation of independently locating according to ring fire detector and keep time and observation equation, carries out adaptive Kalman filter,

Described observation equation adopts pulsar to measure equation,

If the state parameter X of moment t is designated as X (t), the first order derivative of X (t) is described filter state equation is as follows,

$\stackrel{.}{X}\left(t\right)=F\left(t\right)\·X\left(t\right)+W\left(t\right)$

Wherein, W (t) is system state noise matrix,

State Equation Coefficients matrix F (t) is, $F\left(t\right)=\left[\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_x}{\∂y}& \frac{\∂a_x}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_x}{\∂x}& \frac{\∂a_y}{\∂y}& \frac{\∂a_y}{\∂z}& 0& 0& 0& 0& 0& 0\\ \frac{\∂a_z}{\∂x}& \frac{\∂a_z}{\∂y}& \frac{\∂a_z}{\∂z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$

Output module, for output adaptive filtration module gained filter result, comprises track and the clock correction of filtering gained detector.