CN104006813A - Pulsar/starlight angle combination navigation method of high orbit satellite - Google Patents

Abstract

The invention discloses a pulsar/starlight angle combination navigation method of a high orbit satellite. The method comprises the following steps: establishing an orbital dynamic equation of a near earth satellite; establishing a pulsar navigation observation equation; establishing a starlight angle observation equation; using a dynamic filter to process a dynamic model and starlight angle distance information; and fusing the obtained filter result with the pulsar observation information by using a static filter. The combination navigation method uses the dynamic and static nonlinear filter to realize the optimal fusion of the original observation information of the pulsar with the starlight angle distance information, so compared with a pulsar navigation system, a combination navigation system has the advantages of small size, less power consumption short measurement period; and compared with CNS, the combination navigation system can obtain high precision navigation information.

Description

A kind of pulsar/starlight angular distance Combinated navigation method of high rail satellite

Technical field

The present invention is mainly concerned with air navigation aid field, refers in particular to a kind of high rail satellite pulsar/starlight angular distance Combinated navigation method.

Background technology

High rail satellite (orbit altitude is higher than 20000km) has been played the part of important role at aspects such as weather detection, disaster early warning, data relays.Accurate positional information is to guarantee the effectively key factor of work of high rail satellite.At present, the track of high rail satellite determines that work mostly realizes based on ground control station.Yet along with increasing of high rail number of satellite, the burden of ground control station increases gradually.Meanwhile, the day by day complicated of space environment also proposed active demand to the independence of high rail satellite.Having an autonomous navigation system highly reliable, high stable is the prerequisite that realizes high rail satellite autonomy.

Yet because orbit altitude is higher than Navsat, high rail satellite is difficult to receive four above navigation satellite signals simultaneously.Therefore, cannot directly utilize existing GLONASS (Global Navigation Satellite System) to complete the independent navigation work of high rail satellite.In addition, the celestial navigation system based on starlight angular distance (Celestial Navigation System, CNS) is a kind of autonomous navigation system of maturation.The method realizes navigator fix by measuring the sight line angle of the earth and astre fictif, and is applicable to whole terrestrial space.The method also has advantages of that navigation sensor volume is little, the metrical information sampling period is short.Yet the performance of the method reduces along with the increase of orbit altitude.For high rail satellite, CNS only can obtain the positioning precision that is better than 1km.

X ray pulsar is a kind of neutron star of High Rotation Speed, the electromagnetic radiation that it also can produce X-band away from the earth.Rotation period of X ray pulsar is steady in a long-term, the cycle stability degree of some millisecond pulsars current atomic clock that can match in excellence or beauty.The states such as America and Europe take to the method for utilizing X ray pulsar to realize earth satellite independent navigation from the seventies in last century.Because the signal whole day district of X ray pulsar is visible, the performance of X ray pulsar navigation system can not change along with the variation of satellite orbital altitude.Therefore, utilize X ray pulsar navigation system that the positional information of reliable, stable high rail satellite can be independently provided.Yet pulsar signal is faint, and discontinuous.Signal demand continues the accumulation of some minutes, just can obtain believable metrical information.Navigation error within the signal integration time cannot be ignored, and this error also can significantly reduce the performance of X ray pulsar navigation system.Current, mostly utilize area for 1m ²x-ray detector improve the signal to noise ratio (S/N ratio) of pulse signal.Yet, the volume that large area detector brings and power consumption burden, general satellite is difficult to bear.

In view of X ray pulsar navigation system and CNS can be complementary, the Chinese patent document that application number is 200710191527.9 has been announced a kind of Combinated navigation method merging based on much information (following brief note is method 1).Method 1 has merged four kinds of metrical informations, has wherein just comprised starlight angular distance and pulsar metrical information.But method 1 is not used the detector of small size, and do not use Nonlinear Dynamic static filter.Application number is that 200910063267.4 Chinese patent document has been announced a kind of Combinated navigation method (following brief note is method 2) based on pulsar/CNS.The observation data of the method based on adjacent epoch done difference method, adopts Unscented kalman filtering (Unscented Kalman filter, UKF)/H _∞wave filter weakens the impact of pulsar star catalogue site error and clock correction.Method 2 utilizes the starlight elevation angle to carry out integrated navigation, and utilizes the navigation results of federated filter fusion x-ray pulsar navigation system and CNS.

Summary of the invention

The technical problem to be solved in the present invention is: for current high rail satellite, be difficult to realize the problem of autonomous location, the pulsar/starlight angular distance Combinated navigation method of employing based on sound state nonlinear filter completes the Camera calibration work of high rail satellite.

For solving the problems of the technologies described above, the present invention by the following technical solutions:

High rail satellite pulsar/starlight angular distance Combinated navigation method, step is:

A, set up the dynamics of orbits equation of near-earth satellite;

B, set up pulsar navigation observation equation;

C, set up starlight angular distance observation equation;

D, utilize kinetic filter to process kinetic model and starlight angular distance information.Concrete steps are:

1. init state amount with variance battle array P ₀

$\begin{array}{cc}{\hat{x}}_0=E\left(x_0\right)& P_0=E\left[(x_0-{\hat{x}}_0){(x_0-{\hat{x}}_0)}^T\right]\end{array}---\left(8\right)$

2. build state sampled point and weight

near set up a series of sampled point, average and the covariance of these sampled points are respectively and P _k-1.State variable is 6 * 1 dimensions, so 13 sampled point χ _{0, k-1}, χ _{1, k-1}..., χ _{13, k-1}and average weight with variance weight as follows:

$\left\{\begin{array}{c}\\ {\mathrm{\χ}}_{0,k-1}={\hat{x}}_{k-1}\\ {\mathrm{\χ}}_{i,k-1}={\hat{x}}_{k-1}+\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{P_{k-1}}\right)}_{i}i=\mathrm{1,2},L,n\\ {\mathrm{\χ}}_{i+n,k-1}={\hat{x}}_{k-1}-\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{P_{k-1}}\right)}_{i}i=n+1,n+2,L,2n\end{array}\right.---\left(9\right)$

$\left\{\begin{array}{c}{\mathrm{\ω}}_0^m=\frac{\mathrm{\ξ}}{n+\mathrm{\ξ}}\\ {\mathrm{\ω}}_0^c=\frac{\mathrm{\ξ}}{n+\mathrm{\ξ}}+(1-{\mathrm{\α}}^2+\mathrm{\β})\\ {\mathrm{\ω}}_i^m={\mathrm{\ω}}_i^c=\frac{1}{2(n+\mathrm{\ξ})}\end{array}\right.---\left(10\right)$

In formula (9) (10), the dimension that n is quantity of state, ξ=α ²(n+ κ)-n, wherein α is used for controlling sampling point distributions, its value between 0,1, κ=3-n.β is the parameter relevant with state prior distribution, and to Gaussian distribution, β is taken as 2.Work as P _k-1=A ^tduring A, the i that gets A is capable; Work as P _k-1=AA ^ttime, get the i row of A.

3. the kinetic filter time upgrades

The one-step prediction value χ of sampled point _k|k-1for

χ _k|k-1＝f(χ _k-1) (11)

Result after the one-step prediction value weighting of all sampled points for

$x_k^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{\mathrm{\χ}}_{i,k|k-1}---\left(12\right)$

The step estimation variance battle array that quantity of state is estimated for

$P_k^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^c[{\mathrm{\χ}}_{i,k|k-1}-x_k^{-}]{[{\mathrm{\χ}}_{i,k|k-1}-x_k^{-}]}^T+Q_k---\left(13\right)$

4. kinetic filter is measured and is upgraded

The starlight angular distance measured value Z of sampled point one-step prediction _k|k-1 ^stfor

Z _k|k-1 ^st＝h _st(χ _k|k-1) (14)

Result after the one-step prediction starlight angular distance measured value weighting of all sampled points for

$Z_{\mathrm{st},k}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{Z_{i,k|k-1}}^{\mathrm{st}}---\left(15\right)$

Starlight angular distance is measured variance battle array for

The covariance matrix of starlight angular distance measured value and quantity of state for

Filter gain K _{st, k}for

The state estimation value x of kinetic filter _{st, k} ⁺with variance battle array P _{st, k} ⁺be respectively

${x_{\mathrm{st},k}}^{+}=x_k^{-}+K_{\mathrm{st},k}(Z_{\mathrm{st},k}-{\hat{Z}}_{\mathrm{st},k}^{-})---\left(19\right)$

E, utilize static filter to merge kinetic filter result and observations of pulsar information.Concrete steps are

1. static filter initialization

The original state discreet value x of static filter _p,k ^-with variance battle array P _p,k ^-for

x _p,k ^-＝x _st,k ⁺P _p,k ^-＝P _st,k ⁺ (21)

2. build static filter sampled point

The sampled point ε of static filter _{0, k-1}, ε _{1, k-1}..., ε _{13, k-1}for

$\left\{\begin{array}{c}\\ {\mathrm{\ϵ}}_{0,k-1}={x_{p,k}}^{-}\\ {\mathrm{\ϵ}}_{i,k-1}={x_{p,k}}^{-}+\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{{P_{p,k}}^{-}}\right)}_{i}i=\mathrm{1,2},L,n\\ {\mathrm{\ϵ}}_{i+n,k-1}={x_{p,k}}^{-}-\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{{P_{p,k}}^{-}}\right)}_{i}i=n+1,n+2,L,2n\end{array}\right.---\left(22\right)$

Wherein, work as P _p,k ^-=A ^tduring A, the i that gets A is capable; Work as P _p,k ^-=AA ^ttime, get the i row of A.The definition of ξ is identical with formula (10).

3. static filter is measured and is upgraded

The pulsar measured value Z of sampled point one-step prediction _k|k-1 ^pfor

Z _k|k-1 ^p＝h _p(ε _k|k-1) (23)

Result after the one-step prediction pulsar measured value weighting of all sampled points for

$Z_{p,k}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{Z_{i,k|k-1}}^p---\left(24\right)$

Pulsar is measured variance battle array for

The covariance matrix of starlight angular distance measured value and quantity of state for

Filter gain K _p,kfor

The state estimator of static filter with variance battle array be respectively

${\hat{x}}_k=x_k^{-}+K_{p,k}(Z_{p,k}-Z_{p,k}^{-})---\left(28\right)$

By the quantity of state obtaining with variance battle array return to kinetic filter, for k+1 constantly, k=1,2 ...Finally obtain each navigation state estimation value constantly with variance battle array each navigation state estimation value is constantly its positional information in this moment that high rail satellite Autonomous obtains, thereby has realized the independent navigation function of high rail satellite.

Preferably, steps A of the present invention comprises:

At Earth central inertial, be in J2000.0, the dynamics of orbits equation of setting up near-earth satellite is

Wherein, x=[r ^t, v ^t] ^tthe state vector of spacecraft, w=[w _r ^tw _v ^t] ^tfor kinetic model noise, can be modeled as zero-mean white Gaussian noise.The variance battle array of this noise is Q.A=a _tB+ a _nS+ a _t+ a _h.O.Tbe the acceleration that spacecraft is subject to, comprise following several.

1) a _tB=-μ _er/|r| ³the earth disome gravitational acceleration that spacecraft is subject to, wherein μ _eit is the gravitational constant of the earth.

2) it is the non-spherical Gravitational perturbation of the earth.U _nSEcan be expressed as

$U_{\mathrm{NSE}}=-\frac{{\mathrm{\μ}}_E}{\left|r\right|}\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}{\left(\frac{R_e}{\left|r\right|}\right)}^n{J}_n{P}_n\sin \mathrm{\φ}-\frac{{\mathrm{\μ}}_E}{\left|r\right|}\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{R_e}{\left|r\right|}\right)}^n{J}_{n,m}{P}_{\mathrm{nm}}\sin \mathrm{\φ}\cos m(\mathrm{\λ}-{\mathrm{\λ}}_{n,m})---\left(2\right)$

Wherein, R _ebe the radius of the earth, φ and λ are respectively longitude and latitude, P _nand P _nmlegendre polynomial, J _nto be with humorous coefficient, λ _n,mhumorous of field, J _n,mfor humorous of field coefficient.

3) it is trisome perturbation acceleration.μ _ithe gravitational constant of i celestial body, r _ithat i celestial body is with respect to the position vector of the earth.

4) a _h.O.Tit is the higher order term that affects spacecraft acceleration.With respect to the perturbation acceleration of modeling, the impact of these higher order terms can be ignored.

Formula (1) can be noted by abridging and is

x＝f(x)+w (3)

Preferably, step B of the present invention comprises:

The observation equation of pulsar navigation is

Z _p＝h _P(x)+V _p (4)

Wherein, V _pfor observation noise, h _p(x)=[h _p ¹(x) L h _p ^j(x) L h _p ⁿ(x)] ^t, h _p ^j(x) be the observation equation of j pulsar of observation, its expression formula is

Wherein, be the measured value of j pulsar direction vector, for the position of the earth by ephemeris predict, D ₀ ^jbe j the pulsar distance of sun barycenter apart, b is that sun barycenter is with respect to the position of solar system barycenter, μ _sfor the gravitational constant of the sun, c is the light velocity.

In formula (4), V _pcan regard zero-mean white Gaussian noise as, its standard variance is determined by following formula.

${\mathrm{\σ}}_p=\frac{W\sqrt{[B_X+F_X(1-p_f)]d+F_X{p}_f}}{{2F}_X{p}_f\sqrt{{\mathrm{At}}_m}}---\left(6\right)$

Wherein, W is pulse width, B _xfor cosmic background radiation flow, F _xfor pulsar flow, p _ffor the ratio of pulse length to the total cycle length of pulsar, d is pulse width and the ratio of recurrence interval, the area that A is X-ray detector, t _mfor pulsar measuring period.

Preferably, step C of the present invention comprises:

Starlight angular distance observation equation is

$Z_{\mathrm{st}}=h_{\mathrm{st}}\left(x\right)+v_{\mathrm{\α}}=\arccos (-\frac{\mathrm{rgs}}{\left|\right|r\left|\right|})+v_{\mathrm{\α}}---\left(7\right)$

Wherein, the direction vector that s is astre fictif.V _αcan regard zero-mean white Gaussian noise as, its standard variance is determined by the precision of star sensor and optical camera.

Compared with prior art, the invention has the advantages that:

(1) adopt pulsar/starlight angular distance Combinated navigation method, can, in performance pulsar navigation system and CNS advantage, make up the deficiency of two kinds of single operations of navigational system.From formula (6), the measuring accuracy of pulsar navigation and spacecraft orbit are irrelevant, only depend on physical characteristics, Satellite-borne Detector area and the observations of pulsar cycle of pulsar.Therefore, pulsar navigation still can obtain desirable navigation effect to high rail satellite, but corresponding, and the observation cycle of pulsar needs to extend, and needs large-area X-ray detector.From formula (7), the precision of CNS only depends on the precision of navigation sensor, irrelevant with observation time.Therefore, the observation cycle of CNS is much smaller than pulsar navigation system.In addition, according to current technological level, the volume power consumption of star sensor and optical camera all much smaller than X-ray detector (Wang Peng. the autonomous navigation of satellite based on star sensor and attitude are determined method research [D]. Harbin Institute of Technology, 2008.).But because the angle information that CNS measures is relevant with the orbit altitude of spacecraft, CNS is on the low side to the navigation accuracy of high rail satellite.The integrated navigation system that this patent proposes merges the metrical information of having utilized pulsar navigation system, CNS.Than pulsar navigation system, this navigational system can adopt CNS metrical information to navigate at an observations of pulsar in the cycle, thus the measuring period of having reduced whole navigational system.Moreover because fusion has utilized the metrical information of CNS, under the prerequisite allowing in navigation accuracy, the area of X-ray detector also can suitably reduce, thereby can reduce volume and the power consumption of whole navigational system.On the other hand, the relatively high-precision pulsar metrical information of this integrated navigation system utilization merges the low precision measure information of CNS, can obtain the navigation accuracy that is better than CNS.

(2) adopt sound state nonlinear filter can realize the optimum fusion of the original observation information of pulsar and starlight angular distance information.First the conventional information blending algorithms such as existing Federated Filters utilize kinetic model information respectively with the original observation information of pulsar, the calculating of navigating of starlight angular distance information.Then utilize two groups of navigation results to carry out information fusion.Can find out, existing this kind of way can cause the information of kinetic model to be reused, thereby caused having correlativity, information fusion result suboptimum between two groups of navigation results of final fusion.Yet, from formula (9)-(21), can find out, in sound state nonlinear filter, the information of kinetic model is only used in kinetic filter, and static filter is that the result based on kinetic filter is carried out work.So, in sound state nonlinear filter, there is not the situation of information recycling, can realize the optimum fusion of the original observation information of pulsar and starlight angular distance information.

(3) in sound state nonlinear filter, embed the impact that reduces linearization error without mark conversion (unscented transformation).Traditional wave filter for linear system derive (Yang Yuanxi. Multiple Source Sensor is moving, Static Filtering merges navigation [J]. Wuhan University Journal (information science version), 2003, 28 (4): 386-388.), although author utilizes EKF will move static filter and is applied to nonlinear situation, but the linearization error that EKF exists can affect the performance (Chang of moving static filter, L.B., Hu, B.Q., Chang, G.B, et al., Multiple Outliers Suppression Derivative-Free Filter Based on Unscented Transformation[J], Journal of Guidance, Control, and Dynamics, 2012, 35 (6), 1902-1906.).The moving static filter that this patent proposes has adopted the thought without mark conversion in formula (9) and formula (22), does not have linearization procedure, thereby has reduced the impact of linearization error.

Because high rail satellite is being played the part of very important role in daily life, the Camera calibration function of the high rail satellite that the present invention realizes, at least has the advantage of following aspect: 1, can significantly reduce human input and the fund input of ground observing and controlling system; If 2 ground observing and controlling systems lost efficacy, high rail satellite still can obtain the positional information of self by autonomous location, thereby has guaranteed that high rail satellite function realizes.

Accompanying drawing explanation

Fig. 1 is the schematic flow sheet of the concrete application example of the present invention.

Fig. 2 is the effect contrast figure of the concrete application example of the present invention.

Embodiment

As shown in Figure 1, a kind of based on high rail Autonomous Navigation Algorithm in pulsar and starlight angular distance, concrete implementation step is as follows:

(1) set up the dynamics of orbits equation of near-earth satellite;

At Earth central inertial, be in J2000.0, the dynamics of orbits equation of setting up near-earth satellite is

Wherein, x=[r ^t, v ^t] ^tthe state vector of spacecraft, w=[w _r ^tw _v ^t] ^tfor kinetic model noise, can be modeled as zero-mean white Gaussian noise.The variance battle array of this noise is Q.A=a _tB+ a _nS+ a _t+ a _h.O.Tbe the acceleration that spacecraft is subject to, comprise following several.

1) a _tB=-μ _er/|r| ³the earth disome gravitational acceleration that spacecraft is subject to, wherein μ _eit is the gravitational constant of the earth.

2) it is the non-spherical Gravitational perturbation of the earth.UNSE can be expressed as

$U_{\mathrm{NSE}}=-\frac{{\mathrm{\μ}}_E}{\left|r\right|}\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}{\left(\frac{R_e}{\left|r\right|}\right)}^n{J}_n{P}_n\sin \mathrm{\φ}-\frac{{\mathrm{\μ}}_E}{\left|r\right|}\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{R_e}{\left|r\right|}\right)}^n{J}_{n,m}{P}_{\mathrm{nm}}\sin \mathrm{\φ}\cos m(\mathrm{\λ}-{\mathrm{\λ}}_{n,m})---\left(2\right)$

Wherein, R _ebe the radius of the earth, φ and λ are respectively longitude and latitude, P _nand P _nmlegendre polynomial, J _nto be with humorous coefficient, λ _n,mhumorous of field, J _n,mfor humorous of field coefficient.

3) it is trisome perturbation acceleration.μ _ithe gravitational constant of i celestial body, r _ithat i celestial body is with respect to the position vector of the earth.

4) a _h.O.Tit is the higher order term that affects spacecraft acceleration.With respect to the perturbation acceleration of modeling, the impact of these higher order terms can be ignored.

Formula (1) can be noted by abridging and is

x＝f(x)+w (3)

(2) set up pulsar navigation observation equation;

The observation equation of pulsar navigation is

Z _p＝h _P(x)+V _p (4)

Wherein, V _pfor observation noise, h _p(x)=[h _p ¹(x) L h _p ^j(x) L h _p ⁿ(x)] ^t, h _p ^j(x) be the observation equation of j pulsar of observation, its expression formula is

Wherein, be the measured value of j pulsar direction vector, for the position of the earth by ephemeris predict, D ₀ ^jbe j the pulsar distance of sun barycenter apart, b is that sun barycenter is with respect to the position of solar system barycenter, μ _sfor the gravitational constant of the sun, c is the light velocity.

In formula (4), Vp can regard zero-mean white Gaussian noise as, and its standard variance is determined by following formula.

${\mathrm{\σ}}_p=\frac{W\sqrt{[B_X+F_X(1-p_f)]d+F_X{p}_f}}{{2F}_X{p}_f\sqrt{{\mathrm{At}}_m}}---\left(6\right)$

Wherein, W is pulse width, B _xfor cosmic background radiation flow, F _xfor pulsar flow, p _ffor the ratio of pulse length to the total cycle length of pulsar, d is pulse width and the ratio of recurrence interval, the area that A is X-ray detector, t _mfor pulsar measuring period.

(3) set up starlight angular distance observation equation;

Starlight angular distance observation equation is

$Z_{\mathrm{st}}=h_{\mathrm{st}}\left(x\right)+v_{\mathrm{\α}}=\arccos (-\frac{\mathrm{rgs}}{\left|\right|r\left|\right|})+v_{\mathrm{\α}}---\left(7\right)$

Wherein, the direction vector that s is astre fictif.V _αcan regard zero-mean white Gaussian noise as, its standard variance is determined by the precision of star sensor and optical camera.

(4) utilize kinetic filter to process kinetic model and starlight angular distance information.

Concrete steps are

1. init state amount with variance battle array P ₀

$\begin{array}{cc}{\hat{x}}_0=E\left(x_0\right)& P_0=E\left[(x_0-{\hat{x}}_0){(x_0-{\hat{x}}_0)}^T\right]\end{array}---\left(8\right)$

2. build state sampled point and weight

near set up a series of sampled point, average and the covariance of these sampled points are respectively and P _k-1.State variable is 6 * 1 dimensions, so 13 sampled point χ _{0, k-1}, χ _{1, k-1}..., χ _{13, k-1}and average weight with variance weight as follows:

$\left\{\begin{array}{c}\\ {\mathrm{\χ}}_{0,k-1}={\hat{x}}_{k-1}\\ {\mathrm{\χ}}_{i,k-1}={\hat{x}}_{k-1}+\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{P_{k-1}}\right)}_{i}i=\mathrm{1,2},L,n\\ {\mathrm{\χ}}_{i+n,k-1}={\hat{x}}_{k-1}-\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{P_{k-1}}\right)}_{i}i=n+1,n+2,L,2n\end{array}\right.---\left(9\right)$

$\left\{\begin{array}{c}{\mathrm{\ω}}_0^m=\frac{\mathrm{\ξ}}{n+\mathrm{\ξ}}\\ {\mathrm{\ω}}_0^c=\frac{\mathrm{\ξ}}{n+\mathrm{\ξ}}+(1-{\mathrm{\α}}^2+\mathrm{\β})\\ {\mathrm{\ω}}_i^m={\mathrm{\ω}}_i^c=\frac{1}{2(n+\mathrm{\ξ})}\end{array}\right.---\left(10\right)$

In formula (9) (10), the dimension that n is quantity of state, ξ=α ²(n+ κ)-n, wherein α is used for controlling sampling point distributions, its value between 0,1, κ=3-n.β is the parameter relevant with state prior distribution, and to Gaussian distribution, β is taken as 2.Work as P _k-1=A ^tduring A, the i that gets A is capable; Work as P _k-1=AA ^ttime,

3. the kinetic filter time upgrades

The one-step prediction value χ of sampled point _k|k-1for

χ _k|k-1＝f(χ _k-1) (11)

Result after the one-step prediction value weighting of all sampled points for

$x_k^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{\mathrm{\χ}}_{i,k|k-1}---\left(12\right)$

The step estimation variance battle array that quantity of state is estimated for

$P_k^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^c[{\mathrm{\χ}}_{i,k|k-1}-x_k^{-}]{[{\mathrm{\χ}}_{i,k|k-1}-x_k^{-}]}^T+Q_k---\left(13\right)$

4. kinetic filter is measured and is upgraded

The starlight angular distance measured value Z of sampled point one-step prediction _k|k-1 ^stfor

Z _k|k-1 ^st=h _st(χ _k|k-1) result after the one-step prediction starlight angular distance measured value weighting of (14) all sampled points for

$Z_{\mathrm{st},k}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{Z_{i,k|k-1}}^{\mathrm{st}}---\left(15\right)$

Starlight angular distance is measured variance battle array for

The covariance matrix of starlight angular distance measured value and quantity of state for

Filter gain K _{st, k}for

The state estimation value x of kinetic filter _{st, k} ⁺with variance battle array P _{st, k} ⁺be respectively

${x_{\mathrm{st},k}}^{+}=x_k^{-}+K_{\mathrm{st},k}(Z_{\mathrm{st},k}-{\hat{Z}}_{\mathrm{st},k}^{-})---\left(19\right)$

(5) utilize static filter to merge kinetic filter result and observations of pulsar information.Concrete steps are

1. static filter initialization

The original state discreet value x of static filter _p,k ^-with variance battle array P _p,k ^-for

x _p,k ^-＝x _st,k ⁺P _p,k ^-＝P _st,k ⁺ (21)

2. build static filter sampled point

The sampled point ε of static filter _{0, k-1}, ε _{1, k-1}..., ε _{13, k-1}for

$\left\{\begin{array}{c}\\ {\mathrm{\ϵ}}_{0,k-1}={x_{p,k}}^{-}\\ {\mathrm{\ϵ}}_{i,k-1}={x_{p,k}}^{-}+\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{{P_{p,k}}^{-}}\right)}_{i}i=\mathrm{1,2},L,n\\ {\mathrm{\ϵ}}_{i+n,k-1}={x_{p,k}}^{-}-\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{{P_{p,k}}^{-}}\right)}_{i}i=n+1,n+2,L,2n\end{array}\right.---\left(22\right)$

Wherein, work as P _p,k ^-=A ^tduring A, the i that gets A is capable; Work as P _p,k ^-=AA ^ttime, get the i row of A.The definition of ξ is identical with formula (10).

3. static filter is measured and is upgraded

The pulsar measured value Z of sampled point one-step prediction _k|k-1 ^pfor

Z _k|k-1 ^p＝h _p(ε _k|k-1) (23)

Result after the one-step prediction pulsar measured value weighting of all sampled points for

$Z_{p,k}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{Z_{i,k|k-1}}^p---\left(24\right)$

Pulsar is measured variance battle array for

The covariance matrix of starlight angular distance measured value and quantity of state for

Filter gain K _p,kfor

The state estimator of static filter with variance battle array be respectively

${\hat{x}}_k=x_k^{-}+K_{p,k}(Z_{p,k}-Z_{p,k}^{-})---\left(28\right)$

By the quantity of state obtaining with variance battle array return to kinetic filter, for k+1 constantly, k=1,2 ...Finally obtain each navigation state estimation value constantly with variance battle array

The high rail satellite of a certain earth of take is example, and preliminary orbit radical is: semi-major axis 26570.08km, excentricity 0.0047,63.15 ° of orbit inclinations, 136.31 ° of right ascension of ascending node, 211.64 ° of inbreeding point angular distances, 304.82 ° of mean anomalys.The navigation time is 7 days.Fig. 2 has provided the performance comparison situation of the integrated navigation system that this patent proposes, starlight angular distance navigational system, pulsar navigation system.As can be seen from Figure 2, along with the increase of time, the position estimation error curve of integrated navigation system, pulsar navigation system, starlight angular distance navigational system (CNS) all converges near 0 gradually.But can find out, the position estimation error curve convergence of integrated navigation system is fastest, and obtainable final position evaluated error is minimum, and navigation accuracy is the highest.Relative, the graph of errors speed of convergence of pulsar navigation is better than starlight angular distance navigational system, and the final navigation accuracy of pulsar navigation system is also better than starlight angular distance navigational system.Therefore, as shown in Figure 2, the performance of integrated navigation system is better than simple starlight angular distance navigational system and pulsar navigation system.

The above is only the preferred embodiment of the present invention, and protection scope of the present invention is also not only confined to above-described embodiment, and all technical schemes belonging under thinking of the present invention all belong to protection scope of the present invention.Should propose, for those skilled in the art, improvements and modifications without departing from the principles of the present invention, these improvements and modifications also should be considered as protection scope of the present invention.

Claims (4)

1. high rail satellite pulsar/starlight angular distance Combinated navigation method, is characterized in that, comprises step:

A, set up the dynamics of orbits equation of near-earth satellite;

B, set up pulsar navigation observation equation;

C, set up starlight angular distance observation equation;

D, utilize kinetic filter to process kinetic model and starlight angular distance information, wherein:

1. init state amount with variance battle array P ₀

$\begin{array}{cc}{\hat{x}}_0=E\left(x_0\right)& P_0=E\left[(x_0-{\hat{x}}_0){(x_0-{\hat{x}}_0)}^T\right]\end{array}$

2. build state sampled point and weight

near set up a series of sampled point, average and the covariance of described sampled point are respectively and P _k-1; State variable is 6 * 1 dimensions, 13 sampled point χ _{0, k-1}, χ _{1, k-1}..., χ _{13, k-1}and average weight with variance weight as follows:

$\left\{\begin{array}{c}\\ {\mathrm{\χ}}_{0,k-1}={\hat{x}}_{k-1}\\ {\mathrm{\χ}}_{i,k-1}={\hat{x}}_{k-1}+\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{P_{k-1}}\right)}_{i}i=\mathrm{1,2},L,n\\ {\mathrm{\χ}}_{i+n,k-1}={\hat{x}}_{k-1}-\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{P_{k-1}}\right)}_{i}i=n+1,n+2,L,2n\end{array}\right.$

$\left\{\begin{array}{c}{\mathrm{\ω}}_0^m=\frac{\mathrm{\ξ}}{n+\mathrm{\ξ}}\\ {\mathrm{\ω}}_0^c=\frac{\mathrm{\ξ}}{n+\mathrm{\ξ}}+(1-{\mathrm{\α}}^2+\mathrm{\β})\\ {\mathrm{\ω}}_i^m={\mathrm{\ω}}_i^c=\frac{1}{2(n+\mathrm{\ξ})}\end{array}\right.$

In above-mentioned two formulas, the dimension that n is quantity of state, ξ=α ²(n+ κ)-n, wherein α is used for controlling sampling point distributions, its value between 0,1, κ=3-n; β is the parameter relevant with state prior distribution, and to Gaussian distribution, β is taken as 2; Work as P _k-1=A ^tduring A, the i that gets A is capable; Work as P _k-1=AA ^ttime, get the i row of A;

3. the kinetic filter time upgrades

The one-step prediction value χ of sampled point _k|k-1for

χ _k|k-1＝f(χ _k-1)

Result after the one-step prediction value weighting of all sampled points for

$x_k^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{\mathrm{\χ}}_{i,k|k-1}$

The step estimation variance battle array that quantity of state is estimated for

$P_k^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^c[{\mathrm{\χ}}_{i,k|k-1}-x_k^{-}]{[{\mathrm{\χ}}_{i,k|k-1}-x_k^{-}]}^T+Q_k$

4. kinetic filter is measured and is upgraded

The starlight angular distance measured value Z of sampled point one-step prediction _k|k-1 ^stfor

Z _k|k-1 ^st＝h _st(χ _k|k-1)

Result after the one-step prediction starlight angular distance measured value weighting of all sampled points for

$Z_{\mathrm{st},k}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{Z_{i,k|k-1}}^{\mathrm{st}}$

Starlight angular distance is measured variance battle array for

The covariance matrix of starlight angular distance measured value and quantity of state for

Filter gain K _{st, k}for

The state estimation value x of kinetic filter _{st, k} ⁺with variance battle array P _{st, k} ⁺be respectively

${x_{\mathrm{st},k}}^{+}=x_k^{-}+K_{\mathrm{st},k}(Z_{\mathrm{st},k}-{\hat{Z}}_{\mathrm{st},k}^{-})$

E, utilize static filter to merge kinetic filter result and observations of pulsar information, wherein:

1. static filter initialization

The original state discreet value x of static filter _p,k ^-with variance battle array P _p,k ^-for

x _p,k ^-＝x _st,k ⁺P _p,k ^-＝P _st,k ⁺

2. build static filter sampled point

The sampled point ε of static filter _{0, k-1}, ε _{1, k-1}..., ε _{13, k-1}for

$\left\{\begin{array}{c}\\ {\mathrm{\ϵ}}_{0,k-1}={x_{p,k}}^{-}\\ {\mathrm{\ϵ}}_{i,k-1}={x_{p,k}}^{-}+\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{{P_{p,k}}^{-}}\right)}_{i}i=\mathrm{1,2},L,n\\ {\mathrm{\ϵ}}_{i+n,k-1}={x_{p,k}}^{-}-\sqrt{n+\mathrm{\ξ}}\·{\left(\sqrt{{P_{p,k}}^{-}}\right)}_{i}i=n+1,n+2,L,2n\end{array}\right.$

Wherein, work as P _p,k ^-=A ^tduring A, the i that gets A is capable; Work as P _p,k ^-=AA ^ttime, get the i row of A; The definition of ξ with in step D 2.;

3. static filter is measured and is upgraded

The pulsar measured value Z of sampled point one-step prediction _k|k-1 ^pfor

Z _k|k-1 ^p＝h _p(ε _k|k-1)

Result after the one-step prediction pulsar measured value weighting of all sampled points for

$Z_{p,k}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^m{Z_{i,k|k-1}}^p$

Pulsar is measured variance battle array for

The covariance matrix of starlight angular distance measured value and quantity of state for

Filter gain K _p,kfor

The state estimator of static filter with variance battle array be respectively

${\hat{x}}_k=x_k^{-}+K_{p,k}(Z_{p,k}-Z_{p,k}^{-})$

By the quantity of state obtaining with variance battle array return to kinetic filter, for k+1 constantly, k=1,2, Obtain each navigation state estimation value constantly with variance battle array

2. a kind of high rail satellite pulsar/starlight angular distance Combinated navigation method according to claim 1, is characterized in that, described steps A comprises step:

At Earth central inertial, be in J2000.0, the dynamics of orbits equation of setting up near-earth satellite is

Wherein, x=[r ^t, v ^t] ^tthe state vector of spacecraft, w=[w _r ^tw _v ^t] ^tfor kinetic model noise, be modeled as zero-mean white Gaussian noise; The variance battle array of this noise is Q.A=a _tB+ a _nS+ a _t+ a _h.O.Tbe the acceleration that spacecraft is subject to, comprise following several:

1) a _tB=-μ _er/|r| ³the earth disome gravitational acceleration that spacecraft is subject to, wherein μ _eit is the gravitational constant of the earth;

2) it is the non-spherical Gravitational perturbation of the earth; U _nSEcan be expressed as

$U_{\mathrm{NSE}}=-\frac{{\mathrm{\μ}}_E}{\left|r\right|}\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}{\left(\frac{R_e}{\left|r\right|}\right)}^n{J}_n{P}_n\sin \mathrm{\φ}-\frac{{\mathrm{\μ}}_E}{\left|r\right|}\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{R_e}{\left|r\right|}\right)}^n{J}_{n,m}{P}_{\mathrm{nm}}\sin \mathrm{\φ}\cos m(\mathrm{\λ}-{\mathrm{\λ}}_{n,m})$

Wherein, R _ebe the radius of the earth, φ and λ are respectively longitude and latitude, P _nand P _nmlegendre polynomial, J _nto be with humorous coefficient, λ _n,mhumorous of field, J _n,mfor humorous of field coefficient;

3) it is trisome perturbation acceleration; μ _ithe gravitational constant of i celestial body, r _ithat i celestial body is with respect to the position vector of the earth;

4) a _h.O.Tit is the higher order term that affects spacecraft acceleration.

3. a kind of high rail satellite pulsar/starlight angular distance Combinated navigation method according to claim 1, is characterized in that, described step B comprises step:

The observation equation of pulsar navigation is:

Z _p＝h _P(x)+V _p

Wherein, h _p(x)=[h _p ¹(x) L h _p ^j(x) L h _p ⁿ(x)] ^t,

H _p ^j(x) be the observation equation of j pulsar of observation, its expression formula is

Wherein, be the measured value of j pulsar direction vector, for the position of the earth by ephemeris predict, D ₀ ^jbe j the pulsar distance of sun barycenter apart, b is that sun barycenter is with respect to the position of solar system barycenter, μ _sfor the gravitational constant of the sun, c is the light velocity;

V _pfor zero-mean white Gaussian noise, its standard variance is determined by following formula:

${\mathrm{\σ}}_p=\frac{W\sqrt{[B_X+F_X(1-p_f)]d+F_X{p}_f}}{{2F}_X{p}_f\sqrt{{\mathrm{At}}_m}}$

Wherein, W is pulse width, B _xfor cosmic background radiation flow, F _xfor pulsar flow, p _ffor the ratio of pulse length to the total cycle length of pulsar, d is pulse width and the ratio of recurrence interval, the area that A is X-ray detector, t _mfor pulsar measuring period.

4. a kind of high rail satellite pulsar/starlight angular distance Combinated navigation method according to claim 1, is characterized in that, described step C comprises step:

Starlight angular distance observation equation is

$Z_{\mathrm{st}}=h_{\mathrm{st}}\left(x\right)+v_{\mathrm{\α}}=\arccos (-\frac{\mathrm{rgs}}{\left|\right|r\left|\right|})+v_{\mathrm{\α}}$

Wherein, the direction vector that s is astre fictif.