CN103148856B - Swing-by probe autonomous celestial navigation method based on adaptive scale change - Google Patents

Abstract

The invention relates to a swing-by probe autonomous celestial navigation method based on adaptive scale change. The method comprises the following steps: establishing a deep space probe state model and an autonomous celestial navigation system measurement model; acquiring measurement quantity of the autonomous celestial navigation system, and employing an adaptive scale change Unscented Kalman filtering method, so that the time scale suitable for the current model is obtained, and the position and speed of the detector in an inertial coordinate system which takes a target object as the center are obtained; applying the obtained time scale, position and speed to navigation at the next moment. The method belongs to the technical field of aerospace navigation, the time scale can be adaptively changed along with the time change according to the gravity acceleration stressed by the probe, the navigation system model error caused by the fixed time scale time update is reduced, and the method is suitable for a swing-by deep space probe.

Description

A kind of swing-by flight detector celestial self-navigation method based on adaptive scale change

Technical field

The present invention relates to the autonomous navigation method based on the change of model trajectory adaptive scale when deep space probe swing-by flight, is a kind of autonomous navigation method being highly suitable for the deep space probe swing-by flight stage.

Background technology

Survey of deep space technology, as the key character of a national overall national strength and scientific technological advance level and mark, has caused the very big concern of countries in the world.Prelude has been pulled open in the fight of new round survey of deep space.At the beginning of 21 century, sight is focused to the deep space universe beyond the distance earth 380,000 kilometers by each spacefaring nation one after another.Main space flight group is proposed following survey of deep space plan in the interior world for the U.S., European Space Agency, Russia, Japan and India, will carry out manned or based on the unmanned probing of robot to each major planet and satellite thereof.Along with China's carrier rocket and the development of other survey of deep space technology and the raising of economic strength, China has possessed the ability of the even farther planets of the solar system of detection Mars.

Along with the expansion of investigative range and the increase of detection range, the emitted energy needed for deep space probe also significantly increases, and this is the maximum predicament that survey of deep space task faces.In the process of seeking head it off, swing-by flight technology owing to effectively can reduce the speed increment needed for emitted energy and detection mission, and becomes one of technology be most widely used in survey of deep space task.Swing-by flight refers to and utilizes the second gravitation body to change the energy of detector relative centre gravitation body, thus change size or the direction of detector speed, to save emitted energy, with less initial velocity, when without any the acceleration realized when power consumption detector, complete detection mission, even realize the detection mission only leaning on current carrier rocket to complete.

Although swing-by trajectory requires lower to carrier rocket emitted energy, but owing to needing to borrow power to celestial body in detector flight course, and there is strict index request to detector present position and travelling speed during the power of borrowing, therefore require that detector navigational system must can provide accurate navigation information in real time, if navigation information not in time or not accurate enough, the best opportunity of power borrowed by detector by missing navigation celestial body, correct speed correction cannot be carried out, even stray, finally causes mission failure.Therefore the navigation accuracy in swing-by flight stage is the important guarantee that swing-by flight Mission Success is implemented.

The impact of dynamics of orbits model accuracy is subject in the navigation accuracy in deep space probe swing-by flight stage.In independent navigation process, need that integration is carried out to gravitational acceleration and obtain speed and position.Therefore one of them key factor affecting dynamics of orbits model accuracy is exactly the integration step of dynamics of orbits.

Existing autonomous navigation system all adopts fixing integration step to carry out dynamics of orbits model integration, and this method has the following disadvantages: celestial body gravitation acceleration is thought of as constant value by the fixing integration step method that existing autonomous navigation system adopts in an integrating range.But in the swing-by flight survey of deep space task of actual motion, owing to moving closer to along with detector and Mars, detector is subject to Mars gravitation and increases gradually, and the gravitational acceleration of Mars is not constant value, changes greatly in time.The Mars gravitational acceleration changed greatly in time is thought of as constant value by existing fixed step size integration method in an integrating range, and therefore will introduce dynamics of orbits model error, this error increases with the increase of integration step.When the integration step that autonomous navigation system setting is less, at swing-by flight in earlier stage because detector motion changes slowly in time, this large compared with little integration step calculated amount, waste computational resource on detector star; When the integration step that autonomous navigation system setting is larger, when swing-by flight, navigation accuracy is low, cannot follow the tracks of the gravitation change of central body.

Summary of the invention

The technical problem to be solved in the present invention is: overcome existing fixed step size integration method large in swing-by flight calculated amount in early stage, in the low problem of swing-by flight stage navigation accuracy, for swing-by flight deep space probe provides a kind of autonomous navigation method of high-efficiency high-accuracy.

The technical solution adopted for the present invention to solve the technical problems is: to set up centered by target celestial body the deep space probe state model based on the sun and the eight major planets of the solar system gravitation in inertial coordinates system, target celestial body is obtained and angle information measurement amount between satellite and fixed star by optical guidance sensor, set up the measurement model of autonomous navigation system, and under existing time scale, use adaptive scale change Unscented kalman filter method, obtain position and the speed parameter of deep space probe relative target celestial body on the one hand; On the other hand according to the result of Orbit simulation, return suitable time scale, in time renewal process, adopt the time scale of gained to upgrade.Realize changing hours yardstick in time at model trajectory large, and time scale is little when model trajectory changes fast in time, finally for swing-by flight deep space probe provides the navigational parameter such as high precision position, speed efficiently.

Specifically comprise the following steps:

1. set up deep space probe based on the sun and the dynamic (dynamical) state model of the eight major planets of the solar system Attractive Orbit

Deep space probe is set up based on the sun and the dynamic (dynamical) state model of the eight major planets of the solar system Attractive Orbit, the i.e. state model of autonomous astronomical navigation system in inertial coordinates system centered by target celestial body barycenter;

X′(t)＝f ₁(X′(t),t)+w′(t) (5)

In formula, X ' (t)=[x ', y ', z ', v ' _x, v ' _y, v ' _z] ^tfor the state vector of t, x ', y ', z ', v ' _x, v ' _y, v ' _zbe respectively position and the speed of detector three axles in target celestial body barycenter inertial coordinates system, f ₁(X ' (t), t) is mission nonlinear continuous state transfer function, w ' (t)=[0,0,0, w ' _x, w ' _y, w ' _z] ^tfor the state model error of t, w ' _x, w ' _y, w ' _zit is the model error of three axle velocity differentials.

2. set up the autonomous astronomical navigation system measurements model based on starlight angular moment

The angle information θ of target celestial body and two satellites and three background fixed stars _1i, θ _2iand θ _3i(i=1,2,3) expression formula is:

$\left\{\begin{array}{c}{\mathrm{\θ}}_{1i}=\arccos (-{\stackrel{\→}{l}}_{p1}^I{\stackrel{\→}{s}}_{1i}^I)\\ {\mathrm{\θ}}_{2i}=\arccos (-{\stackrel{\→}{l}}_{p2}^I\·{\stackrel{\→}{s}}_{2i}^I)\\ {\mathrm{\θ}}_{3i}=\arccos (-{\stackrel{\→}{l}}_{p3}^I\·{\stackrel{\→}{s}}_{3i}^I)\end{array}\right.---\left(6\right)$

In formula, for in inertial coordinates system target celestial body to the unit vector of detector, for the unit vector of i-th fixed star in target celestial body image in inertial coordinates system; for in inertial coordinates system first satellite of target celestial body to the unit vector of detector, for the unit vector of i fixed star in target celestial body first satellite image in inertial coordinates system; for in inertial coordinates system target celestial body second satellite to the unit vector of detector, for the unit vector of i fixed star in target celestial body second satellite image in inertial coordinates system.

The expression formula can setting up autonomous astronomical navigation system measurements model according to formula (6) is:

Z′(t)＝h ₁[X(t),t]+v ₁(t) (7)

In formula, h ₁[X ' (t), t] measures function, Z ' (t)=[θ continuously for autonomous astronomical navigation mission nonlinear _11,θ ₁₂, θ ₁₃, θ ₂₁, θ ₂₂, θ ₂₃, θ ₃₁, θ ₃₂, θ ₃₃] ^tfor the autonomous astronomical navigation system quantities of t is measured, $v_1={[v_{{\mathrm{\θ}}_{11}},v_{{\mathrm{\θ}}_{12}},v_{{\mathrm{\θ}}_{13}},v_{{\mathrm{\θ}}_{21}},v_{{\mathrm{\θ}}_{22}},v_{{\mathrm{\θ}}_{23}},v_{{\mathrm{\θ}}_{31}},v_{{\mathrm{\θ}}_{31}},v_{{\mathrm{\θ}}_{32}},v_{{\mathrm{\θ}}_{33}}]}^T$ For autonomous astronomical navigation system measurements noise, be respectively and measure θ ₁₁, θ ₁₂, θ ₁₃, θ ₂₁, θ ₂₂, θ ₂₃, θ ₃₁, θ ₃₂, θ ₃₃observational error.

3. the state model in pair step 1 and step 2 and measurement model carry out discretize

X′(k)＝F ₁(X′(k-1),k-1)+W′(k) (8)

Z′(k)＝H ₁(X′(k),k)+V ₁(k) (9)

In formula, k=1,2 ..., F ₁(X ' (k-1), k-1) is f ₁(X ' (t), from kth-1 moment to the nonlinear state transfer function in kth moment after t) discrete, H ₁(X ' (k) k) is h ₁(X ' (t), t) the non-linear measurement function in discrete rear kth moment, W ' (k), V ₁k () is respectively w ' (t), v ₁the equivalent noise in (t) discrete rear kth moment, and W ' (k) and V ₁k () is uncorrelated mutually.

4. obtain measurement amount by autonomous astronomical navigation sensor

Obtained the image of target celestial body by autonomous astronomical navigation sensor, utilize image processing techniques, determine that the pixel of celestial body barycenter is as line coordinates; Through being tied to two-dimensional imaging plane coordinate system from pixel as line coordinates, being tied to three times of sensor surving coordinate system conversions from two-dimensional imaging planimetric coordinates, determine celestial body and the unit vector of background fixed star in sensor coordinate system thereof; Finally calculate the starlight angular distance between celestial body unit vector and background fixed star unit vector.

5. pair autonomous astronomical navigation system carries out adaptive scale change Unscented Kalman filtering

According to state model and the measurement model of the autonomous astronomical navigation system discrete form obtained by step 1 ~ step 3 and obtained the measurement amount obtained by celestial navigation sensor by step 4, carry out autonomous astronomical navigation system self-adaption dimensional variation Unscented Kalman filtering.Concrete steps are as follows:

A. the state vector X ' (n) of initialization (n+1)th filtering autonomous astronomical navigation system and the time scale h of (n+1)th filtering autonomous astronomical navigation system _n, the moment corresponding to n-th filtering n=0,1,2 ..., wherein X ' (n)= ;

As n=0,

${\hat{x}}_0^{\′}=E\left[x_0^{\′}\right],P_0^{\′}=E\left[(x_0^{\′}-{\hat{x}}_0^{\′}){(x_0^{\′}-{\hat{x}}_0^{\′})}^T\right],h_n=h_0---\left(10\right)$ In formula, for three shaft positions and the velocity original value of initial time detector in target celestial body barycenter inertial coordinates system, x ' ₀for initial time three shaft positions of detector and estimated value of speed in target celestial body barycenter inertial coordinates system, P ' ₀for the initial mean squared error matrix of state vector, h _nfor the time scale that (n+1)th filtering uses, h ₀for the time scale of initial time filtering.When n ≠ 0,

${\hat{x}}_n^{\′}=E\left[x_n^{\′}\right],P_n^{\′}=E\left[(x_n^{\′}-{\hat{x}}_n^{\′}){(x_n^{\′}-{\hat{x}}_n^{\′})}^T\right],h_n=h_{n-1}---\left(11\right)$

In formula, be (n+1)th filtering three shaft positions of detector and initial value of speed in target celestial body barycenter inertial coordinates system, x ' _nbe three shaft positions and the velocity estimation value of n-th filtering detector in target celestial body barycenter inertial coordinates system, P ' _nthe initial mean squared error matrix of state vector when being (n+1)th filtering, h _nfor the time scale that (n+1)th filtering uses, h _n-1for the time scale that n-th filtering uses, the moment corresponding to n-th filtering

B. the state vector after steps A initialization and time scale is utilized to carry out time scale self-adaptative adjustment

A) 12 right function K of interlude point in dimension self-adaption adjustment process computing time _stepi

K _step1＝F ₁(X′(n),n) (12)

$K_{\mathrm{stepi}}=F_1(X^{\′}\left({n+a}_{\mathrm{stepi}}h\right)+h\underset{\mathrm{stepj}=0}{\overset{\mathrm{stepi}-1}{\mathrm{\Σ}}}{b}_{\mathrm{stepi},\mathrm{stepj}}{K}_{\mathrm{stepj}},{n+a}_{\mathrm{stpi}}h,),---\left(13\right)$

In formula, the initial state vector in X ' (n) moment corresponding to (n+1)th filtering, F ₁(X ' (n), n) be from n-th filtering to (n+1)th filtering nonlinear state transfer discrete function, stepi=2,3..., 12, wherein b _{stepi, stepj}nonzero value be: b _1,0=2/27, b _2,0=1/36, b _2,1=1/12, b _3,0=1/24, b _3,2=1/8, b _4,0=5/12, b _4,2=-25/16, b _4,3=25/16, b _5,0=1/20, b _5,3=1/4, b _5,4=1/5, b _6,0=-25/108, b _6,3=125/108, b _6,4=-65/27, b _6,5=125/54, b _7,0=31/300, b _7,4=61/225, b _7,5=-2/9, b _7,6=13/900, b _8,0=2b _8,3=-53/6, b _8,4=704/45, b _8,5=-107/9, b _8,6=67/90, b _9,0=-91/108, b _9,3=23/108, b _9,4=-976/135, b _9,5=311/54, b _9,6=-19/60, b _9,7=17/6, b _9,8=-1/12, b _10,0=2383/4100, b _10,3=-341/164, b _10,4=4496/1024, b _10,5=-301/82, b _10,6=2133/4100, b _10,7=45/82, b _10,8=45/162, b _10,9=18/41, b _11,0=3/205, b _11,5=-6/41, b _11,6=-3/205, b _11,7=-3/41, b _11,8=3/41, b _11,9=6/41, b _12,0=-1777/4100, b _12,5=-289/82, b _12,6=2193/4100, b _12,7=51/82, b _12,8=33/164, b _12,9=12/41, b _12,11=1 , a _step0=0, a _step1=2/27, a _step2=1/9, a _step3=1/6, a _step4=5/12, a _step5=1/2, a _step6=5/6, a _step7=1/6, a _step8=2/3, a _step9=1/3, a _step10=1, a _step11=0, a _step12=1, X ' (n)= , h=h _n;

B) 7 rank in dimension self-adaption adjustment process computing time and 8 rank predicted state vectors

$X_{n+1}^{\′}=X_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{10}{\mathrm{\Σ}}}{c}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}---\left(14\right)$

${\tilde{X}}_{n+1}^{\′}=X_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{12}{\mathrm{\Σ}}}{\tilde{c}}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}---\left(15\right)$

In formula, X ' _n+1for 7 rank predicted state vectors in time scale self-adaptative adjustment process, for 8 rank predicted state vectors in time scale self-adaptative adjustment process, c _step0=41/840, c _step1=0, c _step2=0, c _step3=0, c _step4=0, c _step5=34/105, c _step6=9/35, c _step7=9/35, c _step8=9/280, c _step9=9/280, c _step10=41/840, ${\tilde{c}}_{\mathrm{step}0}=0,{\tilde{c}}_{\mathrm{step}1}=0,{\tilde{c}}_{\mathrm{step}2}=0,{\tilde{c}}_{\mathrm{step}3}=0,{\tilde{c}}_{\mathrm{step}4}=0,{\tilde{c}}_{\mathrm{step}5}=034/105,{\tilde{c}}_{\mathrm{step}6}=9/35,$ ${\tilde{c}}_{\mathrm{step}7}=9/35,{\tilde{c}}_{\mathrm{step}8}=9/280,{\tilde{c}}_{\mathrm{step}9}=9/280,{\tilde{c}}_{\mathrm{step}10}=0,{\tilde{c}}_{\mathrm{step}11}=41/840,{\tilde{c}}_{\mathrm{step}12}=41/840;$

C) the local truncation error Δ of computing mode vector

${\mathrm{\Δ}~X}_{n+1}-{\tilde{X}}_{n+1}~\frac{41}{810}(K_0+K_{10}-K_{11}{K}_{12})h---\left(16\right)$

D) judge that the step-length of adaptive scale change Unscented Kalman filtering is the need of adjustment

If Δ >1e-16, then adjust step-length h=h/2, return step a); Otherwise carry out step C;

C. the sampled point of adaptive scale change Unscented Kalman filtering is calculated

Moment corresponding to (n+1)th filtering of autonomous astronomical navigation system state vector near choose a series of sample point, average and the mean squared error matrix of these sample points are respectively with .If state vector is 6 × 1 dimensions, so 13 autonomous astronomical navigation systematic sample point χ ' _{0, n}..., χ ' _{12, n}and weights W ₀' ... W ' ₁₂as follows respectively:

${\mathrm{\χ}}_{0,n}^{\′}={\hat{x}}_n^{\′},W_0^{\′}=-1$

${\mathrm{\χ}}_{j,n}^{\′}={\hat{x}}_n^{\′}+\sqrt{3}{\left(\sqrt{P_n^{\′}}\right)}_j,W_j^{\′}=1/6---\left(17\right)$

${\mathrm{\χ}}_{j+6,n}^{\′}={\hat{x}}_n^{\′}-\sqrt{3}{\left(\sqrt{P_n^{\′}}\right)}_j,W_{j+6}^{\′}=1/6$

In formula, as P ' _n=A ' ^ta ' time, get the jth row of A ', as P ' _n=A ' A ' ^ttime, get the jth row of A ', sampled point χ ' when obtaining (n+1)th filtering _nuniform expression be:

D. the time upgrades

According to the step-length h of step B gained, to all sampled point χ ' of step C gained _nthe time of carrying out renewal, obtains the one-step prediction of autonomous astronomical navigation system state vector estimate mean squared error matrix one-step prediction measure estimate vector

K _step1＝F ₁(χ′ _n,n) (19)

$K_{\mathrm{stepi}}=F_1({\mathrm{\χ}}_{n+a_{\mathrm{stepi}}h}^{\′}+h\underset{\mathrm{stepj}=0}{\overset{\mathrm{stepi}-1}{\mathrm{\Σ}}}{b}_{\mathrm{stepi},\mathrm{stepj}}{K}_{\mathrm{stepj}},n+a_{\mathrm{stepi}}h,)---\left(20\right)$

${\mathrm{\χ}}_{n+1|n}^{\′}={\mathrm{\χ}}_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{10}{\mathrm{\Σ}}}{c}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}---\left(21\right)$

χ′ _n+1|n＝F ₁(χ′ _n,n) (22)

Result after the one-step prediction weighting of all sampled point state vectors of autonomous astronomical navigation system for:

${\hat{x}}_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}{{\mathrm{\χ}}^{\′}}_{j,n+1|n}---\left(23\right)$

In formula, W _j' be the weights of a jth sampled point;

The estimation mean squared error matrix one-step prediction of autonomous astronomical navigation system state vector for:

$P_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\hat{x}}_{n+1}^{\′-}]{[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\hat{x}}_{n+1}^{\′-}]}^T+Q_{n+1}^{\′}---\left(24\right)$

In formula, Q ' _n+1the state model error covariance matrix of autonomous astronomical navigation system when being (n+1)th filtering;

The measurement estimate vector Z ' that autonomous astronomical navigation systematic sampling point is corresponding _n+1|n

Z′ _n+1|n＝H ₁(χ′ _n+1|n,n+1) (25)

The all sampled points of autonomous astronomical navigation system measure estimates weighing vector

${\hat{z}}_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}{Z}_{j,n+1|n}^{\′}---\left(26\right)$

E. renewal is measured

Measure according to step D gained and estimate weighing vector result after the one-step prediction weighting of state vector to measurement mean squared error matrix, state vector measurement amount mean squared error matrix, filter gain, estimated state vector, estimate that mean squared error matrix upgrades:

Autonomous astronomical navigation system measurements mean squared error matrix for:

$P_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]{[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]}^T+R_{n+1}^{\′}---\left(27\right)$

In formula, R ' _n+1the measurement noise covariance matrix of autonomous astronomical navigation system when being (n+1)th filtering;

Autonomous astronomical navigation system state vector quantity measures mean squared error matrix :

$P_{{\hat{x}}_{n+1}{\hat{z}}_{n+1}}^{\′}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\hat{x}}_{n+1}^{\′-}]{[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]}^T---\left(28\right)$

Autonomous astronomical navigation system filter gain K ' _n+1for:

$K_{n+1}^{\′}=P_{{\hat{x}}_{n+1}{\hat{z}}_{n+1}}^{\′}{P}_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′-1}---\left(29\right)$

The estimated state vector of the different boat system of autonomous astronomy with estimation mean squared error matrix P ' _n+1:

${\hat{x}}_{n+1}^{\′}={\hat{x}}_{n+1}^{\′-}+K_{n+1}^{\′}(Z_{n+1}^{\′}-{\hat{z}}_{n+1}^{\′-})---\left(30\right)$

$P_{n+1}^{\′}=P_{n+1}^{\′-}-K_{n+1}^{\′}{P}_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′}{K}_{n+1}^{\′T}---\left(31\right)$

The most at last formula (30) and formula (31) obtain (n+1)th filtering time target celestial body barycenter inertial coordinates system estimated state vector with estimation mean squared error matrix P ' _n+1export, estimated state vector represent the position of detector in target celestial body barycenter inertial coordinates system, velocity information, the estimation mean squared error matrix P ' of output _n+1illustrate the performance that filtering is estimated, and these navigation informations are returned in autonomous astronomical navigation system respectively, for position during the n-th+2 times filtering, speed navigation information.

Principle of the present invention is: first according to the dynamics of orbits of the sun and the eight major planets of the solar system gravitation, set up System State Model in the inertial coordinates system centered by target celestial body; Then, set up the measurement model of autonomous astronomical navigation system, afterwards according to the principle of work of celestial navigation sensor, obtain and the measurement amount measured by celestial navigation sensor processed; After this, because the state model of detector and measurement model all present nonlinear characteristic, and system noise is non-Gaussian noise, therefore adopts Unscented kalman filter method, estimates the detector navigation information of two subsystems; Secondly, because System State Model is fast in detector swing-by flight phase change.Large according to larger set time yardstick solving state error in equation, be difficult to realize high precision swing-by flight, solve according to less set time yardstick, swing-by flight assesses the cost too large early stage, wastes the resource on star.Therefore adopt the Unscented kalman filter method of auto-adaptive time dimensional variation, obtain high precision position and the speed parameter of deep space probe relative target celestial body on the one hand; On the other hand according to the result of Orbit simulation, return suitable time scale, in time renewal process, adopt the time scale of gained to upgrade.Realize changing hours yardstick in time at model trajectory large, and time scale is little when model trajectory changes fast in time, finally for swing-by flight deep space probe provides the navigational parameter such as high precision position, speed efficiently.

The present invention's advantage is compared with prior art: the Unscented kalman filter method adopting adaptive scale change, in time renewal process, according to the result of Orbit simulation, returns suitable time scale.Realize changing hours yardstick in time at model trajectory large, and time scale is little when model trajectory changes fast in time, finally for swing-by flight deep space probe provides the navigational parameter such as high precision position, speed efficiently

Accompanying drawing explanation

Fig. 1 is the autonomous astronomical navigation process flow diagram that the present invention is based on adaptive scale change Unscented Kalman filtering;

Fig. 2 is the imaging schematic diagram of autonomous astronomical navigation system of the present invention celestial navigation sensor used;

Fig. 3 is the process flow diagram of time scale self-adaptative adjustment in the present invention.

Embodiment

As shown in Figure 1, involved by preceding solution can be the interplanetary planets such as Mars, Venus, Jupiter, Saturn by means of power celestial body, below using Mars as embodiment, specific embodiment of the invention process is described:

1. set up deep space probe based on the sun and the dynamic (dynamical) state model of the eight major planets of the solar system Attractive Orbit

First initialized location, speed, if X '=[x ' y ' z ' v _x' v _y' v _z'] ^tfor the state vector in fiery heart inertial coordinates system, x ', y ', z ', v ' _x, v ' _y, v ' _zbe respectively position and the speed of detector three axles in fiery heart inertial coordinates system, above-mentioned variable is all the function relevant with t, and according to the Track desigh of detector, position and the speed initial value of choosing detector are X ' (0).

Deep space probe is set up based on the sun and the dynamic (dynamical) state model of the eight major planets of the solar system Attractive Orbit in fiery heart inertial coordinates system, consider that the sun and the sun such as Mars, the earth and the eight major planets of the solar system are to the graviational interaction of detector, choose fiery heart inertial coordinates system, can obtain the state model of deep space probe in fiery heart inertial coordinates system is:

$\left\{\begin{array}{c}{\stackrel{\·}{x}}^{\′}=v_x^{\′}\\ {\stackrel{\·}{y}}^{\′}=v_y^{\′}\\ {\stackrel{\·}{z}}^{\′}=v_z^{\′}\\ {\stackrel{\·}{v}}_x^{\′}=-{\mathrm{\μ}}_m\frac{x^{\′}}{r_{\mathrm{pm}}^{\′3}}-{\mathrm{\μ}}_s[\frac{x^{\′}-x_s^{\′}}{r_{\mathrm{ps}}^{\′3}}+\frac{x_s^{\′}}{r_{\mathrm{ms}}^{\′3}}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{x^{\′}-x_{i_c}^{\′}}{r_{{\mathrm{pi}}_c}^{\′3}}+\frac{x_{i_c}^{\′}}{r_{{\mathrm{mi}}_c}^3}]+w_x^{\′}\\ {\stackrel{\·}{v}}_y^{\′}=-{\mathrm{\μ}}_m\frac{y^{\′}}{r_{\mathrm{pm}}^{\′3}}-{\mathrm{\μ}}_s[\frac{y^{\′}-y_s^{\′}}{r_{\mathrm{ps}}^{\′3}}+\frac{y_s^{\′}}{r_{\mathrm{ms}}^{\′3}}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{y^{\′}-y_{i_c}^{\′}}{r_{{\mathrm{pi}}_c}^{\′3}}+\frac{y_{i_c}^{\′}}{r_{{\mathrm{mi}}_c}^3}]+w_y^{\′}\\ {\stackrel{\·}{v}}_z^{\′}=-{\mathrm{\μ}}_m\frac{z^{\′}}{r_{\mathrm{pm}}^{\′3}}-{\mathrm{\μ}}_s[\frac{z^{\′}-z_s^{\′}}{r_{\mathrm{ps}}^{\′3}}+\frac{z_s^{\′}}{r_{\mathrm{ms}}^{\′3}}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{z^{\′}-z_{i_c}^{\′}}{r_{{\mathrm{pi}}_c}^{\′3}}+\frac{z_{i_c}^{\′}}{r_{{\mathrm{mi}}_c}^3}]+w_z^{\′}\end{array}\right.---\left(1\right)$

In formula, x ', y ', z ' is detector three shaft positions in fiery heart inertial coordinates system, v ' _x, v ' _y, v ' _zfor detector three axle speed in fiery heart inertial coordinates system, for the differential of detector three shaft positions in fiery heart inertial coordinates system, for the differential of detector three axle speed in fiery heart inertial coordinates system, μ _s, μ _mwith be respectively the sun, Mars and i-th _cthe gravitational constant of planet; R ' _psfor day the heart to the distance of detector; R ' _pmfor Mars is to the distance of detector; R ' _msfor day the heart to the distance of the fiery heart; be i-th _cplanet is to the distance of detector; be i-th _cplanet barycenter is to the distance of the fiery heart; (x ' _s, y ' _s, z ' _s), be respectively the sun, i-th _cthe three shaft position coordinates of planet in fiery heart inertial coordinates system, can be obtained by planet ephemerides according to the time, w _x', w _y', w _z' be respectively the state model error of detector three axle velocity differentials in state model; i _cto represent in the sun and the eight major planets of the solar system i-th from the inside to the outside _cplanet, i _c=1,2,3..., N (i _c≠ 4), N=8, due to i _c=4 represent the 4th planet (Mars), are not therefore included in summation expression formula.

Each variable in formula (1) is all the variable relevant with time t, can be abbreviated as

${\stackrel{\·}{X}}^{\′}\left(t\right)=f_1(X^{\′}\left(t\right),t)+w^{\′}\left(t\right)---\left(2\right)$

In formula, X ' (t)=[x ', y ', z ', v ' _x, v ' _y, v ' _z] ^tfor the state vector of the autonomous astronomical navigation System State Model of t, t differential that () is X ' (t), f ₁(X ' (t), t) is the mission nonlinear continuous state transfer function of state model, w ' (t)=[0,0,0, w ' _x, w ' _y, w ' _z] ^tfor the state model error of t, w ' _x, w ' _y, w ' _zit is the model error of three axle velocity differentials.

2. set up the autonomous astronomical navigation system measurements model based on starlight angular distance;

Angle information θ between Mars and i-th background fixed star _1isize consistent with the expression formula in fiery heart inertial coordinates system in Mars sensor surving coordinate system, can be expressed as:

${\mathrm{\θ}}_{1i}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pm}}^S\·{\stackrel{\→}{s}}_{1i}^S)=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pm}}^I\·{\stackrel{\→}{s}}_{1i}^I)---\left(3\right)$

In formula, for in Mars sensor surving coordinate system from Mars to the unit vector of detector, for in inertial coordinates system, Mars, to the unit vector of detector, can be expressed as:

${\stackrel{\→}{l}}_{\mathrm{pm}}^I=\frac{1}{\sqrt{x^{\′2}+y^{\′2}+z^{\′2}}}\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]---\left(4\right)$

In formula, (x ', y ', z ') be detector three shaft position coordinates in fiery heart inertial coordinates system, for the unit vector of i-th background fixed star in Mars sensor surving coordinate system in Mars image, for the unit vector of i-th background fixed star starlight in inertial coordinates system, i=1,2,3, directly can be obtained by fixed star almanac data storehouse,

In like manner can obtain phobos and angle information θ between Deimos and its i-th background fixed star _2iand θ _3iexpression formula be:

${\mathrm{\θ}}_{2i}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pp}}^I\·{\stackrel{\→}{s}}_{2i}^I)---\left(5\right)$

${\mathrm{\θ}}_{3i}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pd}}^I\·{\stackrel{\→}{s}}_{3i}^I)---\left(6\right)$

The expression formula then can setting up autonomous astronomical navigation measurement model according to formula (4) ~ formula (6) is

Z′(t)＝h ₁[X′(t),t]+v ₁(t) (7)

In formula, Z ' (t)=[θ ₁₁, θ ₁₂, θ ₁₃, θ ₂₁, θ ₂₂, θ ₂₃, θ ₃₁, θ ₃₂, θ ₃₃] ^tfor the autonomous astronomical navigation system quantities of t is measured, θ ₁₁, θ ₁₂, θ ₁₃, θ ₂₁, θ ₂₂, θ ₂₃, θ ₃₁, θ ₃₂, θ ₃₃be respectively the angle information of target celestial body and two satellites and three background fixed stars, h ₁[X ' (t), t] is the non-linear continuous measurement function of autonomous astronomical navigation, for autonomous astronomical navigation system measurements noise, be respectively and measure θ ₁₁, θ ₁₂, θ ₁₃, θ ₂₁, θ ₂₂, θ ₂₃, θ ₃₁, θ ₃₂, θ ₃₃observational error.

3. the state model in pair step 1 and step 2 and measurement model carry out discretize

X′(k)＝F ₁(X′(k-1),k-1)+W′(k) (8)

Z′(k)＝H ₁(X′(k),k)+V ₁(k) (9)

In formula, k=1,2 ..., F ₁(X ' (k-1), k-1) is f ₁(X ' (t), from kth-1 moment to the nonlinear state transfer function in kth moment after t) discrete, H ₁(X ' (k) k) is h ₁(X ' (t), t) the non-linear measurement function in discrete rear kth moment, W ' (k), V ₁k () is w ' (t), v ₁the equivalent noise in (t) discrete rear kth moment, and W ' (k) and V ₁k () is uncorrelated mutually.

4. obtain measurement amount by autonomous astronomical navigation sensor

Fig. 2 gives the imaging schematic diagram of autonomous astronomical navigation system Mars sensor used, and phobos sensor, Deimos sensor imaging process are similar with it.Mars sensor forms, at sensor surving coordinate system O ' X primarily of optical lens and two-dimensional imaging face battle array _cy _cz _cthe middle direction vector along Mars to detector mars sunlight reflection directive Mars sensor, now, the coordinate of Mars in Mars sensor surving coordinate system is (x _c, y _c, z _c); The optical lens of Mars sensor will be imaged in the battle array of two-dimensional imaging face with focal distance f after the light refraction of Mars, the image brightness signal impinged upon on each image-generating unit stores by two-dimensional imaging face battle array.According to the image-forming principle of sensor, the processing procedure that autonomous astronomical navigation system quantities is measured is as follows:

A. the process of celestial image

Because the image of Mars in the battle array of two-dimensional imaging face is not a point, but a circle, by image processing techniques determination Mars image centroids such as barycenter identifications at pixel as line coordinates system O _plx _ply _plin pixel as line (p, l).

B. center-of-mass coordinate is converted to the conversion of two-dimensional imaging plane coordinate system as line coordinates system from pixel

According to pixel as the relation between line coordinates system and two-dimensional imaging plane coordinate system, Mars center-of-mass coordinate is converted to two-dimensional imaging plane coordinate system from pixel as line coordinates system:

$\left[\begin{array}{c}{x}_{2d}\\ y_{2d}\end{array}\right]=K^{-1}(\left[\begin{array}{c}p\\ l\end{array}\right]-\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right])---\left(10\right)$

In formula, (x _2d, y _2d) for Mars barycenter is at two-dimensional imaging plane OX _2dy _2din coordinate, p and l be respectively Mars the pixel of Mars sensor two-dimensional imaging plane and picture line, $K=\left[\begin{array}{cc}{K}_x& K_{\mathrm{yx}}\\ K_{\mathrm{xy}}& K_y\end{array}\right]$ For being transferred to the sensor transition matrix of pixel by millimeter, p ₀and l ₀be respectively Mars sensor center at pixel as line coordinates system OX _ply _plin pixel and picture line.

C. center-of-mass coordinate is converted to the conversion of sensor surving coordinate system from two-dimensional imaging plane coordinate system

According to lens imaging principle, Mars center-of-mass coordinate is converted to sensor surving coordinate system from two-dimensional imaging plane coordinate system:

${\stackrel{\→}{l}}_{\mathrm{pm}}^S=\begin{array}{c}\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]\end{array}=\frac{1}{\sqrt{x_{2d}^2+y_{2d}^2+f^2}}\left[\begin{array}{c}{x}_{2d}\\ y_{2d}\\ -f\end{array}\right]---(11)$

In formula, (x _c, y _c, z _c) for Mars barycenter is at Mars sensor surving coordinate system O ' X _cy _cz _cin coordinate, f is the focal length of Mars sensor, for in Mars sensor surving coordinate system from Mars to the unit vector of detector.

In like manner can obtain the unit vector of i-th background fixed star in Mars sensor surving coordinate system in Mars image i=1,2,3.

D. starlight angular distance is calculated

According to Mars in Mars sensor surving coordinate system to the unit vector of detector with the unit vector of i-th background fixed star in Mars image calculate the starlight angular distance θ between two vectors _1i

${\mathrm{\θ}}_{1i}=\arccos [(-{\stackrel{\→}{l}}_{\mathrm{pm}}^S)\·{\stackrel{\→}{s}}_{1i}^S]---\left(12\right)$

In like manner can obtain phobos and its background fixed star, starlight angular distance θ between Deimos and its background fixed star _2i, θ _3i.

5. pair autonomous astronomical navigation system carries out adaptive scale change Unscented Kalman filtering

According to state model formula (8), autonomous astronomical navigation system measurements modular form (9), the measurement amount formula (12) that obtained by celestial navigation sensor, carry out autonomous astronomical navigation system self-adaption dimensional variation Unscented Kalman filtering.Concrete steps are as follows:

A. the state vector X ' (n) of initialization (n+1)th filtering autonomous astronomical navigation system and the time scale h of (n+1)th filtering autonomous astronomical navigation system _n, the moment corresponding to n-th filtering wherein

As n=0,

${\hat{x}}_0^{\′}=E\left[x_0^{\′}\right],P_0^{\′}=E\left[(x_0^{\′}-{\hat{x}}_0^{\′}){(x_0^{\′}-{\hat{x}}_0^{\′})}^T\right],h_n=h_0---\left(13\right)$

In formula, for three shaft positions and the velocity original value of initial time detector in fiery heart inertial coordinates system, x ' ₀for three shaft positions and the speed actual value of initial time detector in fiery heart inertial coordinates system, P ₀' be the initial mean squared error matrix of state vector, h _nfor the time scale that filtering uses, h ₀for the time scale of initial time filtering.

When n ≠ 0,

${\hat{x}}_n^{\′}=E\left[x_n^{\′}\right],P_n^{\′}=E\left[(x_n^{\′}-{\hat{x}}_n^{\′}){(x_n^{\′}-{\hat{x}}_n^{\′})}^T\right],h_n=h_{n-1}---\left(14\right)$

In formula, be (n+1)th filtering three shaft positions of detector and initial value of speed in fiery heart inertial coordinates system, x ' _nbe three shaft positions and the velocity estimation value of n-th filtering detector in fiery heart inertial coordinates system, P ' _nthe initial mean squared error matrix of state vector when being (n+1)th filtering, h _nfor the time scale that (n+1)th filtering uses, h _n-1for the time scale that n-th filtering uses, the moment corresponding to n-th filtering

B. utilize the state vector after steps A initialization and time scale to carry out time scale self-adaptative adjustment, process flow diagram as shown in Figure 3;

A) the right function K of beans-and bullets shooter in the middle of 12 in dimension self-adaption adjustment process computing time _stepi

K _step1＝F ₁(X′(n),n) (15)

$K_{\mathrm{stepi}}=F_1(X^{\′}(n+a_{\mathrm{stepi}}h)+h\underset{\mathrm{stepj}=1}{\overset{\mathrm{stepi}-1}{\mathrm{\Σ}}}{b}_{\mathrm{stepi},\mathrm{stepj}}{K}_{\mathrm{stepj}},n+a_{\mathrm{stepi}}h,)---\left(16\right)$

In formula, the initial state vector in X ' (n) moment corresponding to (n+1)th filtering, F ₁(X ' (n), n) be from n-th filtering to (n+1)th filtering nonlinear state transfer discrete function, stepi=2,3..., 12, wherein b _{stepi, stepj}nonzero value be: b _1,0=2/27, b _2,0=1/36, b _2,1=1/12, b _3,0=1/24, b _3,2=1/8, b _4,0=5/12, b _4,2=-25/16, b _4,3=25/16, b _5,0=1/20, b _5,3=1/4, b _5,4=1/5, b _6,0=-25/108, b _6,3=125/108, b _6,4=-65/27, b _6,5=125/54, b _7,0=31/300, b _7,4=61/225, b _7,5=-2/9, b _7,6=13/900, b _8,0=2b _8,3=-53/6, b _8,4=704/45, b _8,5=-107/9, b _8,6=67/90, b _9,0=-91/108, b _9,3=23/108, b _9,4=-976/135, b _9,5=311/54, b _9,6=-19/60, b _9,7=17/6, b _9,8=-1/12, b _10,0=2383/4100, b _10,3=-341/164, b _10,4=4496/1024, b _10,5=-301/82, b _10,6=2133/4100, b _10,7=45/82, b _10,8=45/162, b _10,9=18/41, b _11,0=3/205, b _11,5=-6/41, b _11,6=-3/205, b _11,7=-3/41, b _11,8=3/41, b _11,9=6/41, b _12,0=-1777/4100, b _12,5=-289/82, b _12,6=2193/4100, b _12,7=51/82, b _12,8=33/164, b _12,9=12/41, b _12,11=1 a _step0=0, a _step1=2/27, a _step2=1/9, a _step3=1/6, a _step4=5/12, a _step5=1/2, a _step6=5/6, a _step7=1/6, a _step8=2/3, a _step9=1/3, a _step10=1, a _step11=0, a _step12=1, h=h _n;

B) 7 rank in dimension self-adaption adjustment process computing time and 8 rank predicted state vectors

$X_{n+1}^{\′}=X_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{10}{\mathrm{\Σ}}}{c}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}---\left(17\right)$

${\tilde{X}}_{n+1}^{\′}=X_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{12}{\mathrm{\Σ}}}{\tilde{c}}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}---\left(18\right)$

In formula, X ' _n+1for 7 rank predicted state vectors in time scale self-adaptative adjustment process, for 8 rank predicted state vectors in time scale self-adaptative adjustment process, c _step0=41/840, c _step1=0, c _step2=0, c _step3=0, c _step4=0, c _step5=34/105, c _step6=9/35, c _step7=9/35, c _step8=9/280, c _step9=9/280, c _step10=41/840, ${\tilde{c}}_{\mathrm{step}0}=0,{\tilde{c}}_{\mathrm{step}1}=0,{\tilde{c}}_{\mathrm{step}2}=0,{\tilde{c}}_{\mathrm{step}3}=0,{\tilde{c}}_{\mathrm{step}4}=0,{\tilde{c}}_{\mathrm{step}5}=34/105,{\tilde{c}}_{\mathrm{step}6}=9/35,$ ${\tilde{c}}_{\mathrm{step}7}=9/35,{\tilde{c}}_{\mathrm{step}8}=9/280,{\tilde{c}}_{\mathrm{step}9}=9/280,{\tilde{c}}_{\mathrm{step}10}=0,{\tilde{c}}_{\mathrm{step}11}=41/840,{\tilde{c}}_{\mathrm{step}12}=41/840;$

C) the local truncation error Δ of computing mode vector

$\mathrm{\Δ}~X_{n+1}-{\tilde{X}}_{n+1}~\frac{41}{810}(K_0+K_{10}-K_{11}-K_{12})h---\left(19\right)$

D) judge that the step-length of adaptive scale change Unscented Kalman filtering is the need of adjustment

If Δ >1e-16, then adjust step-length h=h/2, return step a); Otherwise carry out step C

C. the sampled point of adaptive scale change Unscented Kalman filtering is calculated

Moment corresponding to (n+1)th filtering of autonomous astronomical navigation system state vector near choose a series of sample point, average and the mean squared error matrix of these sample points are respectively with P ' _n.If state vector is 6 × 1 dimensions, so 13 autonomous astronomical navigation systematic sample point χ ' _{0, n}..., χ ' _{12, n}and weights W ' ₀w ' ₁₂as follows respectively:

${\mathrm{\χ}}_{0,n}^{\′}={\hat{x}}_n^{\′},W_0^{\′}=-1$

${\mathrm{\χ}}_{j,n}^{\′}={\hat{x}}_n^{\′}+\sqrt{3}{\left(\sqrt{P_n^{\′}}\right)}_j,W_j^{\′}=1/6---\left(20\right)$

${\mathrm{\χ}}_{j+6,n}^{\′}={\hat{x}}_n^{\′}-\sqrt{3}{\left(\sqrt{P_n^{\′}}\right)}_j,W_{j+6}^{\′}=1/6$

In formula, as P ' _n=A ' ^ta ' time, get the jth row of A ', as P ' _n=A ' A ' ^ttime, get the jth row of A ', sampled point χ ' when obtaining (n+1)th filtering _nuniform expression be:

D. the time upgrades

According to the step-length h of step B gained, to all sampled point χ ' of step C gained _nthe time of carrying out renewal, obtains the one-step prediction of autonomous astronomical navigation system state vector estimate mean squared error matrix one-step prediction measure estimate vector

$K_{\mathrm{step}1}=F_1({\mathrm{\χ}}_n^{\′},n)---\left(22\right)$

$K_{\mathrm{stepi}}=F_1({\mathrm{\χ}}_{n+a_{\mathrm{stepi}}h}^{\′}+h\underset{\mathrm{stepj}=0}{\overset{\mathrm{stepi}-1}{\mathrm{\Σ}}}{b}_{\mathrm{stepi},\mathrm{stepj}}{K}_{\mathrm{stepj}},n+a_{\mathrm{stepi}}h,),---\left(23\right)$

${\mathrm{\χ}}_{n+1|n}^{\′}={\mathrm{\χ}}_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{10}{\mathrm{\Σ}}}{c}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}---\left(24\right)$

χ′ _n+1|n＝F ₁(χ′ _n,n) (25)

Result after the one-step prediction weighting of all sampled point state vectors of autonomous astronomical navigation system for:

${\hat{x}}_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}{\mathrm{\χ}}_{j,n+1|n}^{\′}---\left(26\right)$

In formula, W _j' be the weights of a jth sampled point;

The estimation mean squared error matrix one-step prediction of autonomous astronomical navigation system state vector for:

$P_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\widehat{\mathrm{\χ}}}_{n+1}^{\′-}]{[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\hat{x}}_{n+1}^{\′-}]}^T+Q_{n+1}^{\′}---\left(27\right)$

In formula, Q ' _n+1the state model error covariance matrix of autonomous astronomical navigation system when being (n+1)th filtering;

The measurement estimate vector Z ' that autonomous astronomical navigation systematic sampling point is corresponding _n+1|n

Z′ _n+1|n＝H ₁(χ′ _n+1|n,n+1) (28)

The all sampled points of autonomous astronomical navigation system measure estimates weighing vector

${\hat{z}}_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}{Z}_{j,n+1|n}^{\′}---\left(29\right)$

E. renewal is measured

Measure according to step D gained and estimate weighing vector , state vector one-step prediction weighting after result to measurement mean squared error matrix, state vector measurement amount mean squared error matrix, filter gain, estimated state vector, estimate that mean squared error matrix upgrades:

Autonomous astronomical navigation system measurements mean squared error matrix for:

$P_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]{[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]}^T+R_{n+1}^{\′}---\left(30\right)$

In formula, R ' _n+1the measurement noise covariance matrix of autonomous astronomical navigation system when being the (n+1)th subwave;

Autonomous astronomical navigation system state vector quantity measures mean squared error matrix

$P_{{\hat{x}}_{n+1}{\hat{z}}_{n+1}}^{\′}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\hat{x}}_{n+1}^{\′-}]{[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]}^T---\left(31\right)$

Autonomous astronomical navigation system filter gain K ' _n+1for:

$K_{n+1}^{\′}=P_{{\hat{x}}_{n+1}{\hat{z}}_{n+1}}^{\′}{P}_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′-1}---\left(32\right)$

The estimated state vector of autonomous astronomical navigation system with estimation mean squared error matrix P ' _n+1for

${\hat{x}}_{n+1}^{\′}={\hat{x}}_{n+1}^{\′-}+K_{n+1}^{\′}(Z_{n+1}^{\′}-{\hat{z}}_{n+1}^{\′-})---\left(33\right)$

$P_{n+1}^{\′}=P_{n+1}^{\′-}-K_{n+1}^{\′}{P}_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′}{K}_{n+1}^{\′T}---\left(34\right)$

The most at last formula (33) and formula (34) obtain (n+1)th filtering time in fiery heart inertial coordinates system estimated state vector with estimation mean squared error matrix P ' _n+1export, estimated state vector represent the position of detector in fiery heart inertial coordinates system, velocity information, the estimation mean squared error matrix P ' of output _n+1illustrate the performance that filtering is estimated, and these navigation informations are returned in autonomous astronomical navigation system respectively, for position during the n-th+2 times filtering, speed navigation information.

The content be not described in detail in instructions of the present invention belongs to the known prior art of professional and technical personnel in the field.

Claims (1)

1. the swing-by flight detector celestial self-navigation method based on adaptive scale change, it is characterized in that: state model and the measurement model of first setting up deep space probe, autonomous astronomical navigation system is utilized to obtain the measurement amount of relative target celestial body, estimate to obtain the position of detector centered by target celestial body in inertial coordinates system and speed by adaptive scale change Unscented filtering, specifically comprise the following steps:

1. deep space probe is set up based on the sun and the eight major planets of the solar system Attractive Orbit dynamics state model;

The state model of deep space probe in target celestial body barycenter inertial coordinates system is:

${\stackrel{.}{X}}^{\′}\left(t\right)=f_1(X^{\′}\left(t\right),t)+w^{\′}\left(t\right)$

In formula, X ' (t)=[x ', y ', z ', v ' _x, v ' _y, v ' _z] ^tfor the autonomous astronomical navigation system state vector of t, for the differential of X ' (t), x ', y ', z ', v ' _x, v ' _y, v ' _zbe respectively detector in the position of three axles in target celestial body barycenter inertial coordinates system and speed, f ₁(X ' (t), t) is mission nonlinear continuous state transfer function, w ' (t)=[0,0,0, w ' _x, w ' _y, w ' _z] ^tfor the state model error of t, w ' _x, w ' _y, w ' _zit is the model error of three axle velocity differentials;

2. the autonomous astronomical navigation system measurements model based on starlight angular distance is set up;

Z′(t)＝h ₁[X′(t),t]+v ₁(t) (1)

In formula, Z ' (t)=[θ ₁₁, θ ₁₂, θ ₁₃, θ ₂₁, θ ₂₂, θ ₂₃, θ ₃₁, θ ₃₂, θ ₃₃] ^tfor the autonomous astronomical navigation system quantities of t is measured, θ ₁₁, θ ₁₂, θ ₁₃, θ ₂₁, θ ₂₂, θ ₂₃, θ ₃₁, θ ₃₂, θ ₃₃be respectively the angle information of target celestial body and two satellites and three background fixed stars, h ₁[X ' (t), t] is the non-linear continuous measurement function of autonomous astronomical navigation, for autonomous astronomical navigation system measurements noise, be respectively and measure θ ₁₁, θ ₁₂, θ ₁₃, θ ₂₁, θ ₂₂, θ ₂₃, θ ₃₁, θ ₃₂, θ ₃₃observational error;

3. to step 1. with step 2. in state model and measurement model carry out discretize, obtain state model and the measurement model of autonomous astronomical navigation system discrete form;

X′(k)＝F ₁(X′(k-1),k-1)+W′(k)

Z′(k)＝H ₁(X′(k),k)+V ₁(k)

In formula, k=1,2 ..., F ₁(X ' (k-1), k-1) is f ₁(X ' (t), from kth-1 moment to the nonlinear state transfer function in kth moment after t) discrete, H ₁(X ' (k) k) is h ₁(X ' (t), t) the non-linear measurement function in discrete rear kth moment, W ' (k), V ₁k () is respectively w ' (t), v ₁the equivalent noise in (t) discrete rear kth moment, and W ' (k) and V ₁k () is uncorrelated mutually;

4. measurement amount is obtained by autonomous astronomical navigation sensor;

5. according to state model and the measurement model of autonomous astronomical navigation system discrete form that are 3. obtained by step and 4. obtained the measurement amount obtained by celestial navigation sensor by step, carry out autonomous astronomical navigation system self-adaption dimensional variation Unscented Kalman filtering, concrete steps are as follows:

A. the state vector of initialization (n+1)th filtering autonomous astronomical navigation system with the time scale h of (n+1)th filtering _n, the moment corresponding to n-th filtering n=0,1,2 ..., wherein

B. the state vector after steps A initialization and time scale is utilized to carry out time scale self-adaptative adjustment;

A) 12 right function K of interlude point in dimension self-adaption adjustment process computing time _stepi

K _step1＝F ₁(X′(n),n)

$K_{\mathrm{stepi}}=F_i(X^{\′}(n+a_{\mathrm{stepi}}h)+h\underset{\mathrm{stepj}=0}{\overset{\mathrm{stepi}-1}{\mathrm{\Σ}}}{b}_{\mathrm{stepi},\mathrm{stepj}}{K}_{\mathrm{stepj}},n+a_{\mathrm{stepi}}h,)$

In formula, the initial state vector in X ' (n) moment corresponding to (n+1)th filtering, F ₁(X ' (n), n) be from n-th filtering to (n+1)th filtering nonlinear state transfer discrete function, stepi=2,3..., 12, wherein b _{stepi, stepj}nonzero value be:

b _1,0＝2/27,b _2,0＝1/36,b _2,1＝1/12,b _3,0＝1/24,b _3,2＝1/8,b _4,0＝5/12,

b _4,2＝-25/16,b _4,3＝25/16,b _5,0＝1/20,b _5,3＝1/4,b _5,4＝1/5,b _6,0＝-25/108,b _6,3＝125/108,

b _6,4＝-65/27,b _6,5＝125/54,b _7,0＝31/300,b _7,4＝61/225,b _7,5＝-2/9,b _7,6＝13/900,b _8,0＝2

b _8,3＝-53/6,b _8,4＝704/45,b _8,5＝-107/9,b _8,6＝67/90,b _9,0＝-91/108,b _9,3＝23/108,b _9,4＝-976/135,

b _9,5＝311/54,b _9,6＝-19/60,b _9,7＝17/6,b _9,8＝-1/12,b _10,0＝2383/4100,b _10,3＝-341/164,

b _10,4＝4496/1024,b _10,5＝-301/82,b _10,6＝2133/4100,b _10,7＝45/82,b _10,8＝45/162,

b _10,9＝18/41,b _11,0＝3/205,b _11,5＝-6/41,b _11,6＝-3/205,b _11,7＝-3/41,b _11,8＝3/41,b _11,9＝6/41,

b _12,0＝-1777/4100,b _12,5＝-289/82,b _12,6＝2193/4100,b _12,7＝51/82,b _12,8＝33/164,b _12,9＝12/41,

b _12,11＝1 a _step0＝0,a _step1＝2/27,a _step2＝1/9,a _step3＝1/6,a _step4＝5/12,a _step5＝1/2,a _step6＝5/6,a _step7＝1/6,a _step8＝2/3,a _step9＝1/3,a _step10＝1,a _step11＝0,a _step12＝1， h＝h _n；

B) 7 rank in dimension self-adaption adjustment process computing time and 8 rank predicted state vectors

$X_{n+1}^{\′}=X_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{10}{\mathrm{\Σ}}}{c}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}---\left(2\right)$

${\tilde{X}}_{n+1}^{\′}=X_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{12}{\mathrm{\Σ}}}{\tilde{c}}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}---\left(3\right)$

In formula, for 7 rank predicted state vectors in time scale self-adaptative adjustment process, for 8 rank predicted state vector c in time scale self-adaptative adjustment process _step0=41/840, c _step1=0, c _step2=0, c _step3=0, c _step4=0, c _step5=34/105, c _step6=9/35, c _step7=9/35, c _step8=9/280, c _step9=9/280, c _step10=41/840, ${\tilde{c}}_{\mathrm{step}0}=0,{\tilde{c}}_{\mathrm{step}1}=0,{\tilde{c}}_{\mathrm{step}2}=0,{\tilde{c}}_{\mathrm{step}3}=0,{\tilde{c}}_{\mathrm{step}4}=0,{\tilde{c}}_{\mathrm{step}5}=34/105,{\tilde{c}}_{\mathrm{step}6}=9/35,$ ${\tilde{c}}_{\mathrm{step}7}=9/35,{\tilde{c}}_{\mathrm{step}8}=9*280,{\tilde{c}}_{\mathrm{step}9}=9/280,{\tilde{c}}_{\mathrm{step}10}=0,{\tilde{c}}_{\mathrm{step}11}=41/840,{\tilde{c}}_{\mathrm{step}12}=41/840;$

C) the local truncation error Δ of computing mode vector

$\mathrm{\Δ}~X_{n+1}^{\′}-{\tilde{X}}_{n+1}^{\′}~\frac{41}{810}(K_0+K_{10}-K_{11}-K_{12})h---\left(4\right)$

D) judge that the step-length of adaptive scale change Unscented Kalman filtering is the need of adjustment

If Δ >1e-16, then adjust step-length h=h/2, return step a); Otherwise carry out step C;

C. the sampled point of adaptive scale change Unscented Kalman filtering is calculated;

Moment corresponding to (n+1)th filtering of autonomous astronomical navigation system state vector near choose a series of sample point, average and the mean squared error matrix of these sample points are respectively with P ' _nif state vector is 6 × 1 dimensions, so 13 autonomous astronomical navigation systematic sample point χ ' ₀, n ..., χ ' ₁₂, n and weights W thereof ' ₀w ' ₁₂as follows respectively:

${\mathrm{\χ}}_{0,n}^{\′}={\hat{x}}_n^{\′},$ W′ ₀＝-1

${\mathrm{\χ}}_{j,n}^{\′}={\hat{x}}_n^{\′}+\sqrt{3}{\left(\sqrt{P_n^{\′}}\right)}_j,$ W′ _j＝1/6

${\mathrm{\χ}}_{j+6,n}^{\′}={\hat{x}}_n^{\′}-\sqrt{3}{\left(\sqrt{P_n^{\′}}\right)}_j,$ W′ _j+6＝1/6

In formula, as P ' _n=A ^{' T}a ' time, get the jth row of A ', as P ' _n=A ' A ^{' T}time, get the jth row of A ', sampled point χ ' when obtaining (n+1)th filtering _nuniform expression be:

D. the time upgrades;

According to the step-length h of step B gained, to all sampled point χ ' of step C gained _nthe time of carrying out renewal, obtains the one-step prediction of autonomous astronomical navigation system state vector estimate mean squared error matrix one-step prediction measure estimate vector

K _step1＝F ₁(χ′ _n,n)

$K_{\mathrm{stepi}}=F_1({\mathrm{\χ}}_{n+a_{\mathrm{stepi}}h}^{\′}+h\underset{\mathrm{stepj}=0}{\overset{\mathrm{stepi}-1}{\mathrm{\Σ}}}{b}_{\mathrm{stepi},\mathrm{stepj}}{K}_{\mathrm{stepj}},n+a_{\mathrm{stepi}}h,)$

${\mathrm{\χ}}_{n+1|n}^{\′}={\mathrm{\χ}}_n^{\′}+h\underset{\mathrm{stepi}=0}{\overset{10}{\mathrm{\Σ}}}{c}_{\mathrm{stepi}}{K}_{\mathrm{stepi}}$

χ′ _n+1|n＝F ₁(χ′ _n,n)

Result after the one-step prediction weighting of all sampled point state vectors of autonomous astronomical navigation system for:

${\hat{x}}_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}{{\mathrm{\χ}}^{\′}}_{j,n+1|n}$

In formula, W ' _jfor the weights of a jth sampled point;

The estimation mean squared error matrix one-step prediction of autonomous astronomical navigation system state vector for:

$P_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\hat{x}}_{n+1}^{\′-}]{[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\hat{x}}_{n+1}^{\′-}]}^T+Q_{n+1}^{\′}$

In formula, Q ' _n+1the state model error covariance matrix of autonomous astronomical navigation system when being (n+1)th filtering;

The measurement estimate vector Z ' that autonomous astronomical navigation systematic sampling point is corresponding _n+1|n

Z′ _n+1|n＝H ₁(χ′ _n+1|n,n+1)

The all sampled points of autonomous astronomical navigation system measure estimates weighing vector

${\hat{z}}_{n+1}^{\′-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_J^{\′}{Z}_{j,n+1|n}^{\′}$

E. renewal is measured;

Measure according to step D gained and estimate weighing vector result after the one-step prediction weighting of state vector to measurement mean squared error matrix, state vector measurement amount mean squared error matrix, filter gain, estimated state vector, estimate that mean squared error matrix upgrades;

Autonomous astronomical navigation system measurements mean squared error matrix for:

$P_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]{[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]}^T+R_{n+1}^{\′}$

In formula, R ' _n+1the measurement noise covariance matrix of autonomous astronomical navigation system when being (n+1)th filtering;

Autonomous astronomical navigation system state vector quantity measures mean squared error matrix

$P_{{\hat{x}}_{n+1}{\hat{z}}_{n+1}}^{\′}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j^{\′}[{\mathrm{\χ}}_{j,n+1|n}^{\′}-{\hat{x}}_{n+1}^{\′-}]{[Z_{j,n+1|n}^{\′}-{\hat{z}}_{n+1}^{\′-}]}^T$

Autonomous astronomical navigation system filter gain K ' _n+1for:

$K_{n+1}^{\′}=P_{{\hat{x}}_{n+1}{\hat{z}}_{n+1}}^{\′}{P}_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′-1}$

The estimated state vector of autonomous astronomical navigation system with estimation mean squared error matrix P ' _n+1for

${\hat{x}}_{n+1}^{\′}={\hat{x}}_{n+1}^{\′-}+K_{n+1}^{\′}(Z_{n+1}^{\′}-{\hat{z}}_{n+1}^{\′-})$

$P_{n+1}^{\′}=P_{n+1}^{\′-}-K_{n+1}^{\′}{P}_{{\hat{z}}_{n+1}{\hat{z}}_{n+1}}^{\′}{K}_{n+1}^{\′T}$

Estimated state vector during (n+1)th filtering of final acquisition in fiery heart inertial coordinates system with estimation mean squared error matrix P ' _n+1export, estimated state vector represent the position of detector in fiery heart inertial coordinates system, velocity information, the estimation mean squared error matrix P ' of output _n+1illustrate the performance that filtering is estimated, and these navigation informations are returned in autonomous astronomical navigation system respectively, for position during the n-th+2 times filtering, speed navigation information.