CN100533065C - Interstellar cruising self-navigation method based on multi-star road sign - Google Patents

Abstract

The invention relates to the space cruising autonomous navigation method based on plural globe route sign. It includes the following steps: programming the navigation using minor planet table, processing navigation minor planet image, and determine the autonomous path. The center of figure of plural minor planet determines the path of craft. The error of craft path is very low. The location error could be closed to 100Km and the speed error could be bellow 0.3m/s. It could satisfy the path accuracy request of the detector.

Description

The interspace autonomous navigation method that cruises based on many celestial bodies road sign

Technical field

Autonomous Orbit in the present invention relates to that a kind of survey of deep space is interspace and cruising is determined method.

Background technology

The survey of deep space autonomous navigation technology is an important technology of Space Science and Technology development, and its basic goal is the independence that realizes that the survey of deep space track is determined to reduce the cost of Operating Complexity and reduction task.Particularly in interspace cruising phase, autonomous navigation technology is particularly important.Independent navigation can be divided into three major types by the method for its acquired information: 1) inertial navigation: continuously acceleration or the speed of measuring is obtained positional information to time integral; 2) celestial navigation/GPS: directly obtain positional information by metrical information; 3) optical guidance/radar graphic coupling navigation: handle metrical information and obtain navigation information.In existing optical navigation method, having of being associated with the interspace navigation of cruising is following several: 1) utilize star sensor to carry out independent navigation and starlight is simulated half things emulation; 2) utilize refraction of star sensor starlight or horizon instrument to measure the starlight angular distance and carry out orbit determination; 3) promptly utilizing directly gentle sensitively starlight to reflect responsive Horizon based on information fusion combines, satellite is carried out independent navigation; 4) based on the independent navigation of the image information of little celestial body or the autonomous navigation method of day ground month information.Existing optical navigation method has been obtained certain progress along with the propelling of time, but also has navigation procedure complexity, low, the not easy-operating shortcoming of precision.

Summary of the invention

The objective of the invention is also to have navigation procedure complexity, low, the not easy-operating problem of precision, a kind of interspace autonomous navigation method that cruises based on many celestial bodies road sign is provided for solving existing optical navigation method.The present invention has simple, the easy realization of navigation procedure, advantage of high precision.Technical scheme of the present invention is realized by following steps: the first, the planning navigation is tabulated with asteroid:

Navigation camera parameter and navigation asteroid choice criteria

With aircraft distance be 0～1e6km, with the relative velocity of aircraft be 0～7km/s,

Absolute magnitude is 0～12, and solar angle is ± 40 degree,

The tabulation of planning asteroid is based on the whole asteroidal tabulation of date that obtains from U.S. JPL laboratory, includes asteroidal ephemeris, the physical characteristics of magnitude.According to the asteroid ephemeris, adopt two track body computing method to obtain the position and the speed of asteroid current time.The distance of asteroid and aircraft and relative velocity are obtained by the nominal track and the asteroidal position and speed of aircraft.Thereby choose the celestial body that satisfies constraint requirements.The navigation initial time, aircraft is to the day orientation, and the body coordinate system overlaps with orbital coordinate system.Thereby adopt two vectors to decide the appearance method, the expectation attitude is obtained by position of aircraft and asteroid position.According to initial attitude and expectation attitude, and velocity of rotation calculates attitude switching time, and seclected time, many shortest asteroids were used as nautical star, formed the tabulation of navigation asteroid.The second, handle navigation asteroid image: its objective is in order to obtain the centre of form of asteroid image.The centre of form determines that method is in two steps:

At first determine in " the bright heart ".

$p_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}},$ $l_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}}$

In the formula, p _IjBe the brightness value (i is a pixel, and j is the picture line) of respective pixel, p _CbAnd l _CbBe the bright heart of navigation asteroid image that calculates, m * n=1024 * 1024th, the resolution of camera.

Determine " centre of form " according to the bright heart then.

p ₀＝p _cb-γR _ccosφ，l ₀＝l _cb-γR _csinφ

In the formula, $\mathrm{\φ}=a\tan \left(\frac{A_{\mathrm{cy}}}{A_{\mathrm{cx}}}\right),$ A _cBe that camera coordinates is the direction of the sun, by the inertia sight line vector A that points to the sun _IObtain through coordinate conversion; A _c=T _CIA _I, T wherein _CIBe tied to the transition matrix of camera coordinates system for inertial coordinate.R _cBe the pixel of asteroid radius correspondence, $R_c=\frac{{\mathrm{RK}}_{x}f}{\mathrm{\ρ}},$ R is asteroidal radius, and ρ is detector and asteroidal distance, and f is the focal length of camera, K _xIt is the conversion from the long measure to the pixel; γ is a deviation factors, value between 0 and 1, the centre of form deviation that expression radius of target error delta R causes, it embodies formula and is

$\mathrm{\γ}=\frac{3\mathrm{\π\ΔR}}{16}\left[\frac{\sin \mathrm{\α}(1+\cos \mathrm{\α})}{(\mathrm{\π}-\mathrm{\α})\cos \mathrm{\α}+\sin \mathrm{\α}}\right]$

α is sun phasing degree.

Carry out the calculating of centroid pixel at last,

Inertia sight line vector

, utilize transition matrix T _CIRotate to camera coordinate system,

${\stackrel{\→}{V}}_C=\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\\ V_{C3}\end{array}\right]={\mathrm{TC}}_I{\stackrel{\→}{V}}_I$

The sight line vector

Determine, be transformed in the camera focal plane of two dimension,

$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{V_{C3}}\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\end{array}\right]$

Wherein,

F camera focal length mm,

V _C1, V _C2, V _C3The sight line vector is component in camera system,

X, y sight line vector focal plane inner projection mm,

Then, by x and the y deviation delta x that optical distortion causes, Δ y,

$\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]=Q{\left[\begin{array}{cccccc}{v}_1& v_2& v_3& v_4& v_5& v_6\end{array}\right]}^T$

Wherein,

$Q=\left[\begin{array}{cccccc}-\mathrm{yr}& {\mathrm{xr}}^2& {-\mathrm{yr}}^3& {\mathrm{xr}}^4& \mathrm{xy}& x^2\\ \mathrm{xr}& {\mathrm{yr}}^2& {\mathrm{xr}}^3& {\mathrm{yr}}^4& y^2& \mathrm{xy}\end{array}\right]$

R=x ²+ y ², v ' s is the light distortion coefficient.The picture position x ' that revises, y ' is expressed as

$\left[\begin{array}{c}x\′\\ y\′\end{array}\right]=\left[\begin{array}{c}x+\mathrm{\Δx}\\ y+\mathrm{\Δy}\end{array}\right]$

At last, rectangular coordinate system is transformed into pixel and picture line,

$\left[\begin{array}{c}p\\ l\end{array}\right]=\left[\begin{array}{ccc}{K}_x& K_{\mathrm{xy}}& K_{\mathrm{xxy}}\\ K_{\mathrm{yx}}& K_y& K_{\mathrm{yxy}}\end{array}\right]\left[\begin{array}{c}x\′\\ y\′\\ x\′y\′\end{array}\right]+\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right]$

The K battle array be from the millimeter to the pixel/as the transition matrix of line, p ₀, l ₀Be center pixel and the picture line of CCD.

Three, carrying out autonomous Orbit based on many asteroids determines: utilize the asteroidal centre of form of above-mentioned many of obtaining to determine spacecraft orbit.

1) definition reference orbit parameter:

$X^{*}\left(t\right)={\left[\begin{array}{cccccc}x& y& z& \stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}\end{array}\right]}^T,$ Wherein, r=[x y z] ^TWith $\stackrel{\‾}{v}={\left[\begin{array}{ccc}\stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}\end{array}\right]}^T$ Be respectively position and the velocity of detector in heliocentric ecliptic coordinate system;

The orbit parameter of upgrading:

X ' (t)=X ^*(t)+and Δ X (t), wherein, Δ X (t) is for estimating the track correction.Because reference orbit parameter and true track are more or less the same, and are linear in a period of time inner orbit correction so, promptly at the orbit parameter correction amount X of a time point (t ₀) utilize state transition matrix to be mapped to any putting At All Other Times on the t linearly, i.e. Δ X (t)=Φ (t) Δ X (t ₀)

Wherein, preset time t state-transition matrix Φ satisfy

$\stackrel{\·}{\mathrm{\Φ}}=\frac{\∂\stackrel{\·}{X}\left(t\right)}{\∂X\left(t\right)}\frac{\∂X\left(t\right)}{\∂X\left(t_0\right)}=\mathrm{A\Φ}\left(t\right)$

In the formula, $A=\frac{\∂\stackrel{\·}{X}\left(t\right)}{\∂X\left(t\right)},$ Φ(t ₀)＝I _6×6。Obtain the function of time of Φ by the integration of these equations.

The state equation of system is written as:

$\stackrel{\&CenterDot;}{X}=\left[\begin{array}{c}\stackrel{\&OverBar;}{v}\\ -\frac{{\mathrm{\&mu;}}_s}{r^3}\stackrel{\&OverBar;}{r}+\underset{i=1}{\overset{n_p}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_i[\frac{{\stackrel{\&OverBar;}{r}}_{\mathrm{ri}}}{r_{\mathrm{ri}}^3}-\frac{{\stackrel{\&OverBar;}{r}}_{p1}}{{\mathrm{<figure-callout\; id="1"\; label="r\; p"\; filenames="C200610010105E00153.png"\; state="\{\{state\}\}">r</figure-callout>}}_p^1}]-\frac{\mathrm{AG}}{{\mathrm{mr}}^3}\stackrel{\&OverBar;}{r}+\frac{k}{m}\stackrel{\&OverBar;}{T}+\stackrel{\&OverBar;}{a}\end{array}\right]$

2) determine observing matrix:

For a certain epoch of observation, observing matrix is

$H=\left[\begin{array}{cccccc}\&PartialD;p/{\&PartialD;\mathrm{<figure-callout\; id="1"\; label="x"\; filenames="C200610010105E00153.png"\; state="\{\{state\}\}">x</figure-callout>}}_1& \&PartialD;p/{\&PartialD;x}_2& \&PartialD;p/{\&PartialD;x}_3& 0& 0& 0\\ \&PartialD;l/{\&PartialD;x}_1& \&PartialD;l/{\&PartialD;x}_2& \&PartialD;l/{\&PartialD;x}_3& 0& 0& 0\end{array}\right]$

In the formula, the partial differential relevant with speed component all is zero, and this is that l is only relevant with the asteroidal relatively position of the instantaneous detector of taking pictures owing to p, and irrelevant with the speed of detector.

3) unify epoch of observation:

Because many asteroid pixels, pixels obtaining are constantly different, so need utilize state-transition matrix that the observing matrix of each epoch of observation is transformed to same epoch of observation, have

$\tilde{H}=\mathrm{H\&Phi;}$

Wherein It is exactly the observation partial differential matrix of carving at a time.

4) realize that track is definite:

For the numerical stability that minimizes round-off error and guarantee algorithm, utilize the recursion weighted least square algorithm that decomposes based on the UD covariance to determine the track of detector.Measurement residual error Y wherein is that the pixel corresponding with calculating asteroid center that prediction obtains that obtain of Flame Image Process is poor, so just utilize many asteroidal observation datas (pixel at asteroid center), revise the orbit parameter of certain epoch of observation by the multistep recursion.At first provide a covariance matrix P ₀, weighting obsdervations matrix W, vector Y are calculated the difference of barycenter for observation barycenter and nominal trajectory forecast.Initial point state filtering equations is in batches found the solution vector

With formal covariance matrix P:

$\hat{q}={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}{\tilde{H}}^T\mathrm{WY}$

$P={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}$

Wherein,

$W=\left[\begin{array}{cc}1/{\mathrm{\&sigma;}}_0^2& 0\\ 0& 1/{\mathrm{\&sigma;}}_0^2\end{array}\right]$

${\mathrm{\&sigma;}}_0=\frac{{\tan }^{-1}(R/\mathrm{\&rho;})}{13\&times;{10}^{-6}}$

R asteroid supposition radius ρ asteroid scope.

Beneficial effect of the present invention:

For verifying effect of the present invention, be example to survey the asteroidal one section section of the cruising track of Ivar1627, carry out mathematical simulation, simulation parameter:

1) the nominal track of detector utilizes numerical integration to obtain.At J2000.0 day heart ecliptic inertial coordinates system,

The initial position [1.13984 * 10 of detector ¹¹-1.14516 * 10 ¹¹-6.73821 * 10 ⁶] m,

The initial velocity of detector is [2.71473 * 10 ⁴1.883397 * 10 ⁴-2.59150] m/s;

2) measuring accuracy: the attitude error variance is 10 ^-12Rad ², pixel error is: pixel 0.1, as line 0.1;

3) survey frequency: measurement data is exported with 450 seconds sampling interval, utilizes the asteroid image of selecting and planning to carry out track and determines;

4) initial error: detector position is 1 * 10 in the variance of all directions error ¹⁴m ², all directions velocity error variance is 10 ⁴m ²/ s ²

5) ephemeris error: asteroidal ephemeris is 1 * 10 in the variance of all directions error ¹⁰m ²

Reference locus state parameter initial value is

X _Reference=X _Nominal+ [1 * 10 ⁷1 * 10 ⁷1 * 10 ⁷-10 10 10] '

Covariance matrix

$P_0=\left[\begin{array}{cccccc}{10}^{14}& 0& 0& 0& 0& 0\\ 0& {10}^{14}& 0& 0& 0& 0\\ 0& 0& {10}^{14}& 0& 0& 0\\ 0& 0& 0& 100& 0& 0\\ 0& 0& 0& 0& 100& 0\\ 0& 0& 0& 0& 0& 100\end{array}\right]$

Uncertainty and weighting obsdervations matrix

$\begin{array}{cc}{\mathrm{\&sigma;}}_0=0.001& W=\left[\begin{array}{cc}1/{0.001}^2& 0\\ 0& 1/{0.001}^2\end{array}\right]\end{array}$

m＝500Kg；k＝1；A＝50(m ²)；G＝4.65e-6(N/m ²)

Under the situation that does not have measuring error and asteroid ephemeris error, the relation that observes between asteroid number and the track evaluated error is as shown in table 1.As can be seen, along with the increase that observes the asteroid number, the track estimated accuracy progressively improves; For can observe 7 asteroids the time, after 4 Post Orbits were determined, positional precision can reach 4.3m, and velocity accuracy 0.0013m/s is very near the track true value; For can observe 8 asteroids the time, 2 Post Orbits determine that the back is just very near the track true value; But for can observe 6 following asteroids the time, track determines that speed of convergence is slow, and can not precision determine track.These show, need to observe 7 asteroids at least and just can carry out track and determine.

The relation of table 1 observation asteroid number and track evaluated error

Detector track under the situation that does not have measuring error and asteroid ephemeris error is determined error as shown in Figure 2, as can be seen, can fully accurately determine the track of detector based on the autonomous optical navigation algorithm of many asteroid images.Autonomous Orbit under the emulation assumed conditions is determined the result as shown in Figure 3, and as can be seen, site error is near 100km, and velocity error can satisfy the requirement of the detector section of cruising to trajectory accuracy in the 0.3m/s scope.

Description of drawings

Fig. 1 is the geometrical constraint coordinate diagram that camera is taken pictures, and Fig. 2 is that the track under no measuring error and the ephemeris error situation is determined error curve diagram, and Fig. 3 utilizes the track of many asteroid images to determine error curve diagram.

Embodiment

Embodiment one: (referring to Fig. 1～Fig. 3) step of present embodiment is as follows:

The first, the planning navigation is tabulated with asteroid:

1) nautical star choice criteria,

Navigation camera parameter and navigation asteroid choice criteria

Choice criteria	Scope
		With the aircraft distance	0～1000000km
Relative velocity with aircraft	0～7km/s
		Absolute magnitude	0～12
Solar angle	± ± 40 degree

2) the planning asteroid tabulation first step: can obtain the asteroid tabulation from U.S. JPL laboratory, include asteroidal ephemeris, physical characteristicss such as magnitude.Adopt two track body computing method to obtain current position of asteroid and speed.Can obtain by nominal track and asteroidal ephemeris with the distance and the relative velocity of aircraft.

3) second step of planning asteroid tabulation: aircraft is to the day orientation, and the body coordinate system overlaps (seeing accompanying drawing 1) with orbital coordinate system.Adopt two vectors to decide the appearance method, the expectation attitude is obtained by position of aircraft and asteroid position.Initial attitude and expectation Attitude Calculation go out switching time, seclected time the shortest 12;

The second, handle navigation asteroid image:

Its objective is in order to obtain the asteroid image centre of form.The centre of form determines that method is in two steps:

The first step is determined in " the bright heart ".

$p_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{ip}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}},$ $l_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{jp}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}}$

In the formula, p _IjBe the brightness value (i is a pixel, and j is the picture line) of respective pixel, p _CbAnd l _CbBe the bright heart of navigation asteroid image that calculates, m * n=1024 * 1024th, the resolution of camera.

In second step, determine " centre of form ".

p ₀＝p _cb-γR _ccosφ，l ₀＝l _cb-γR _csinφ

In the formula, $\mathrm{\&phi;}=a\tan \left(\frac{A_{\mathrm{cy}}}{A_{\mathrm{cx}}}\right),$ A _cBe that camera coordinates is the direction of the sun, can be by the inertia sight line vector A that points to the sun _IObtain through coordinate conversion; A _c=T _CIA _I, T wherein _CIBe tied to the transition matrix of camera coordinates system for inertial coordinate.R _cBe the pixel of asteroid radius correspondence, $R_c=\frac{{\mathrm{RK}}_{x}f}{\mathrm{\&rho;}},$ R is asteroidal radius, and ρ is detector and asteroidal distance, and f is the focal length of camera, K _xIt is the conversion from the long measure to the pixel; γ is a deviation factors, value between 0 and 1, the centre of form deviation that expression radius of target error delta R causes, it embodies formula and is

$\mathrm{\&gamma;}=\frac{3\mathrm{\&pi;\&Delta;R}}{16}\left[\frac{\sin \mathrm{\&alpha;}(1+\cos \mathrm{\&alpha;})}{(\mathrm{\&pi;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}+\sin \mathrm{\&alpha;}}\right]$

α is sun phasing degree.

The 3rd step, the calculating of pixel,

If obtained an inertia sight line vector

Utilize transition matrix T _CIRotate to camera coordinate system, ${\stackrel{\&RightArrow;}{V}}_C$

${\stackrel{\&RightArrow;}{V}}_C=\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\\ V_{C3}\end{array}\right]={\mathrm{TC}}_I{\stackrel{\&RightArrow;}{V}}_I$

The sight line vector

In case determine, need be transformed in the camera focal plane of two dimension,

$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{V_{C3}}\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\end{array}\right]$

Wherein,

F camera focal length mm,

V _C1, V _C2, V _C3The sight line vector is component in camera system,

X, y sight line vector focal plane inner projection mm,

Then, by x and the y deviation delta x that optical distortion causes, Δ y,

$\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\end{array}\right]=Q{\left[\begin{array}{cccccc}{v}_1& v_2& v_3& v_4& v_5& v_6\end{array}\right]}^T$

Wherein,

$Q=\left[\begin{array}{cccccc}-\mathrm{yr}& {\mathrm{xr}}^2& {-\mathrm{yr}}^3& {\mathrm{xr}}^4& \mathrm{xy}& x^2\\ \mathrm{xr}& {\mathrm{yr}}^2& {\mathrm{xr}}^3& {\mathrm{yr}}^4& y^2& \mathrm{xy}\end{array}\right]$

R=x ²+ y ², v ' s is the light distortion coefficient.The picture position x ' that revises, y ' is expressed as:

$\left[\begin{array}{c}x\&prime;\\ y\&prime;\end{array}\right]=\left[\begin{array}{c}x+\mathrm{\&Delta;x}\\ y+\mathrm{\&Delta;y}\end{array}\right]$

At last, rectangular coordinate system is transformed into pixel and picture line,

$\left[\begin{array}{c}p\\ l\end{array}\right]=\left[\begin{array}{ccc}{K}_x& K_{\mathrm{xy}}& K_{\mathrm{xxy}}\\ K_{\mathrm{yx}}& K_y& K_{\mathrm{yxy}}\end{array}\right]\left[\begin{array}{c}x\&prime;\\ y\&prime;\\ x\&prime;y\&prime;\end{array}\right]+\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right]$

The K battle array be from the millimeter to the pixel/as the transition matrix of line, p ₀, l ₀Be center pixel and the picture line of CCD.

Three, carrying out autonomous Orbit based on many celestial bodies road sign (many asteroids) determines: the attitude of utilizing the asteroidal centre of form of above-mentioned many of obtaining and background fixed star to determine, determine track.

1) definition reference orbit parameter:

$X^{*}\left(t\right)={\left[\begin{array}{cccccc}x& y& z& \stackrel{\&CenterDot;}{x}& \stackrel{\&CenterDot;}{y}& \stackrel{\&CenterDot;}{z}\end{array}\right]}^T,$ Wherein, r=[x y z] ^TWith $\stackrel{\&OverBar;}{v}={\left[\begin{array}{ccc}\stackrel{\&CenterDot;}{x}& \stackrel{\&CenterDot;}{y}& \stackrel{\&CenterDot;}{z}\end{array}\right]}^T$ Be respectively position and the velocity of detector in heliocentric ecliptic coordinate system;

The orbit parameter of upgrading:

X ' (t)=X ^*(t)+and Δ X (t), wherein, Δ X (t) is for estimating the track correction.If reference orbit parameter and true track are more or less the same, and are linear in a period of time inner orbit correction so, promptly at the orbit parameter correction amount X of a time point (t ₀) can utilize state transition matrix to be mapped to any putting At All Other Times on the t linearly, i.e. Δ X (t)=Φ (t) Δ X (t ₀)

Wherein, preset time t state-transition matrix Φ satisfy

$\stackrel{\&CenterDot;}{\mathrm{\&Phi;}}=\frac{\&PartialD;\stackrel{\&CenterDot;}{X}\left(t\right)}{\&PartialD;X\left(t\right)}\frac{\&PartialD;X\left(t\right)}{\&PartialD;X\left(t_0\right)}=\mathrm{A\&Phi;}\left(t\right)$

In the formula, $A=\frac{\&PartialD;\stackrel{\&CenterDot;}{X}\left(t\right)}{\&PartialD;X\left(t\right)},$ Φ(t ₀)＝I _6×6。The state equation of system can be written as:

$\stackrel{\&CenterDot;}{X}=\left[\begin{array}{c}\stackrel{\&OverBar;}{v}\\ -\frac{{\mathrm{\&mu;}}_s}{r^3}\stackrel{\&OverBar;}{r}+\underset{i=1}{\overset{n_p}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_i[\frac{{\stackrel{\&OverBar;}{r}}_{\mathrm{ri}}}{r_{\mathrm{ri}}^3}-\frac{{\stackrel{\&OverBar;}{r}}_{p1}}{{\mathrm{<figure-callout\; id="1"\; label="r\; p"\; filenames="C200610010105E00153.png"\; state="\{\{state\}\}">r</figure-callout>}}_p^1}]-\frac{\mathrm{AG}}{{\mathrm{mr}}^3}\stackrel{\&OverBar;}{r}+\frac{k}{m}\stackrel{\&OverBar;}{T}+\stackrel{\&OverBar;}{a}\end{array}\right]$

In the formula:

R detector position vector in heliocentric coordinates; V detector velocity in heliocentric coordinates; r _PiI the perturbation position vector of planetary in heliocentric coordinates considered in the emulation, water, gold,, the Perturbation Effect of fire, wooden five stars and the moon, r _RiThe position vector of relative i the perturbation planetary of detector, wherein r _Ri=r _Pi-r; μ _δSolar gravitation constant GM; μ _iThe gravitational constant of i perturbation celestial body; n _pThe quantity of perturbation celestial body; The useful area of A detector; G day luminous flux constant; T propulsion system thrust vectoring; K thrust scale-up factor is approximately 1; M detector quality; The coriolis acceleration that a detector is total;

In following formula, first on the right is the acceleration that right centrosome solar gravitation causes; Second is the summation of trisome gravitational acceleration; The 3rd is sun optical pressure; The 4th is the propulsive acceleration of propulsion system; Last expression acts on other coriolis accelerations on the detector.The effect of preceding two gravitational accelerations is directly, but the non-gravitation factor that acts on the detector is to be worth discussing.For sun optical pressure, adopted simple spherical model for detector in the formula, but in practice, the feature area of detector is mainly by two solar array decisions, and detector body only accounts for a very little part.

2) determine observing matrix:

For a certain epoch of observation, observing matrix is

$H=\left[\begin{array}{cccccc}\&PartialD;p/{\&PartialD;x}_1& \&PartialD;p/{\&PartialD;x}_2& \&PartialD;p/{\&PartialD;x}_3& 0& 0& 0\\ \&PartialD;l/{\&PartialD;x}_1& \&PartialD;l/{\&PartialD;x}_2& \&PartialD;l/{\&PartialD;x}_3& 0& 0& 0\end{array}\right]$

In the formula, the partial differential relevant with speed component all is zero, and this is that l is only relevant with the asteroidal relatively position of the instantaneous detector of taking pictures owing to p, and irrelevant with the speed of detector.

3) unify epoch of observation:

Because many asteroid pixels, pixels obtaining are constantly different, so need utilize state-transition matrix that the observing matrix of each epoch of observation is transformed to same epoch of observation, have

$\tilde{H}=\mathrm{H\&Phi;}$

4) realize that track is definite:

For the numerical stability that minimizes round-off error and guarantee algorithm, utilize the recursion weighted least square algorithm that decomposes based on the UD covariance to determine the track of detector.Measurement residual error Y wherein is that the pixel corresponding with calculating asteroid center that prediction obtains that obtain of Flame Image Process is poor, so just can utilize many asteroidal observation datas (pixel at asteroid center), revise the orbit parameter of certain epoch of observation by the multistep recursion.At first provide a covariance matrix P ₀, weighting obsdervations matrix W, vector Y are calculated the difference of barycenter for observation barycenter and nominal trajectory forecast.Initial point state filtering equations is in batches found the solution vector

With formal covariance matrix P:

$\hat{q}={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}{\tilde{H}}^T\mathrm{WY}$

$P={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}$

Wherein,

$W=\left[\begin{array}{cc}1/{\mathrm{\&sigma;}}_0^2& 0\\ 0& 1/{\mathrm{\&sigma;}}_0^2\end{array}\right]$

${\mathrm{\&sigma;}}_0=\frac{{\tan }^{-1}(R/\mathrm{\&rho;})}{13\&times;{10}^{-6}}$

R asteroid supposition radius ρ asteroid scope.

Claims (1)

1, a kind of interspace autonomous navigation method that cruises based on many celestial bodies road sign, realized by following steps based on the technical scheme of the interspace autonomous navigation method that cruises of many celestial bodies road sign: the first, the planning navigation is tabulated with asteroid;

Navigation camera parameter and navigation asteroid choice criteria: with the aircraft distance be 0～1e6km, with the relative velocity of aircraft be 0～7km/s, absolute magnitude is 0～12, solar angle be ± 40 to spend;

The tabulation of planning asteroid is based on the whole asteroidal tabulation of date that obtains from U.S. JPL laboratory, include asteroidal ephemeris, the physical characteristics of magnitude, according to the asteroid ephemeris, adopt two track body computing method to obtain the position and the speed of asteroid current time, the distance of asteroid and aircraft and relative velocity are obtained by the nominal track and the asteroidal position and speed of aircraft, thereby choose the celestial body that satisfies constraint requirements, the navigation initial time, aircraft is to the day orientation, and the body coordinate system overlaps with orbital coordinate system, thereby adopts two vectors to decide the appearance method, and the expectation attitude is obtained by position of aircraft and asteroid position, according to initial attitude and expectation attitude, and velocity of rotation calculates attitude switching time, and seclected time, many shortest asteroids were used as nautical star, formed the tabulation of navigation asteroid; It is characterized in that second, handle navigation asteroid image: the centre of form determines that method is in two steps:

At first determine in " the bright heart ",

$p_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{ip}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}},$ $l_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{jp}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}}$

In the formula, p _IjBe the brightness value of respective pixel, i is a pixel, and j is the picture line, p _CbAnd l _CbBe the bright heart of navigation asteroid image that calculates, m * n=1024 * 1024th, the resolution of camera;

Determine " centre of form " according to the bright heart then,

p ₀＝p _cb-γR _ccosφ，l ₀＝l _cb-γR _csinφ

In the formula, $\mathrm{\&phi;}=a\tan \left(\frac{A_{\mathrm{cy}}}{A_{\mathrm{cx}}}\right),$ A _cBe that camera coordinates is the direction of the sun, by the inertia sight line vector A that points to the sun _IObtain through coordinate conversion; A _c=T _CIA _I, T wherein _CIFor inertial coordinate is tied to the transition matrix that camera coordinates is, R _cBe the pixel of asteroid radius correspondence, $R_c=\frac{{\mathrm{RK}}_{x}f}{\mathrm{\&rho;}},$ R is asteroidal radius, and ρ is detector and asteroidal distance, and f is the focal length of camera, K _xIt is the conversion from the long measure to the pixel; γ is a deviation factors, value between 0 and 1, the centre of form deviation that expression radius of target error delta R causes, it embodies formula and is

$\mathrm{\&gamma;}=\frac{3\mathrm{\&pi;\&Delta;R}}{16}\left[\frac{\sin \mathrm{\&alpha;}(1+\cos \mathrm{\&alpha;})}{(\mathrm{\&pi;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}+\sin \mathrm{\&alpha;}}\right]$

α is sun phasing degree,

Carry out the calculating of centroid pixel at last,

Inertia sight line vector

Utilize transition matrix T _CIRotate to camera coordinate system,

${\stackrel{\&RightArrow;}{V}}_C=\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\\ V_{C3}\end{array}\right]=T_{\mathrm{CI}}{\stackrel{\&RightArrow;}{V}}_I$

The sight line vector

Determine, be transformed in the camera focal plane of two dimension,

$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{V_{C3}}\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\end{array}\right]$

Wherein,

F camera focal length mm,

V _C1, V _C2, V _C3The sight line vector is component in camera system,

X, y sight line vector focal plane inner projection mm,

Then, by x and the y deviation delta x that optical distortion causes, Δ y,

$\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\end{array}\right]=Q{\left[\begin{array}{cccccc}{v}_1& v_2& v_3& v_4& v_5& v_6\end{array}\right]}^T$

Wherein,

$Q=\left[\begin{array}{cccccc}-\mathrm{yr}& {\mathrm{xr}}^2& -{\mathrm{yr}}^3& {\mathrm{xr}}^4& \mathrm{xy}& x^2\\ \mathrm{xr}& {\mathrm{yr}}^2& {\mathrm{xr}}^3& {\mathrm{yr}}^4& y^2& \mathrm{xy}\end{array}\right]$

R=x ²+ y ², v ₁v ₂v ₃v ₄v ₅v ₆It is the light distortion coefficient; The picture position x ' that revises, y ' is expressed as

$\left[\begin{array}{c}x\&prime;\\ y\&prime;\end{array}\right]=\left[\begin{array}{c}x+\mathrm{\&Delta;x}\\ y+\mathrm{\&Delta;y}\end{array}\right]$

At last, rectangular coordinate system is transformed into pixel and picture line,

$\left[\begin{array}{c}p\\ l\end{array}\right]=\left[\begin{array}{ccc}{K}_x& K_{\mathrm{xy}}& K_{\mathrm{xxy}}\\ K_{\mathrm{yx}}& K_y& K_{\mathrm{yxy}}\end{array}\right]\left[\begin{array}{c}x\&prime;\\ y\&prime;\\ x\&prime;y\&prime;\end{array}\right]+\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right]$

The K battle array be from the millimeter to the pixel/as the transition matrix of line, p ₀, l ₀Be center pixel and the picture line of CCD;

Three, carrying out autonomous Orbit based on many asteroids determines: utilize the asteroidal centre of form of above-mentioned many of obtaining to determine spacecraft orbit;

1) definition reference orbit parameter:

$X^{*}\left(t\right)={\left[\begin{array}{cccccc}x& y& z& \stackrel{.}{x}& \stackrel{.}{y}& \stackrel{.}{z}\end{array}\right]}^T,$ Wherein, r=[x y z] ^TWith $\stackrel{\&OverBar;}{v}={\left[\begin{array}{ccc}\stackrel{.}{x}& \stackrel{.}{y}& \stackrel{.}{z}\end{array}\right]}^T$ Be respectively position and the velocity of detector in heliocentric ecliptic coordinate system;

The orbit parameter of upgrading:

X ' (t)=X ^*(t)+and Δ X (t), wherein, Δ X (t) because reference orbit parameter and true track are more or less the same, is linear in a period of time inner orbit correction, promptly at the orbit parameter correction amount X of a time point (t for estimating the track correction so ₀) utilize state transition matrix to be mapped to any putting At All Other Times on the t linearly, i.e. Δ X (t)=Φ (t) Δ X (t ₀)

Wherein, preset time t state-transition matrix Φ satisfy

$\stackrel{.}{\mathrm{\&Phi;}}=\frac{\&PartialD;\stackrel{.}{X}\left(t\right)}{\&PartialD;X\left(t\right)}\frac{\&PartialD;X\left(t\right)}{\&PartialD;X\left(t_0\right)}=\mathrm{A\&Phi;}\left(t\right)$

In the formula, $A=\frac{\&PartialD;\stackrel{.}{X}\left(t\right)}{\&PartialD;X\left(t\right)},$ Φ (t ₀)=I _{6 * 6}, obtaining the function of time of Φ by the integration of these equations, the state equation of system is written as:

$\stackrel{.}{X}=\left[\begin{array}{c}\stackrel{\&OverBar;}{v}\\ -\frac{{\mathrm{\&mu;}}_s}{r^3}\stackrel{\&OverBar;}{r}+\underset{i=1}{\overset{n_p}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_i[\frac{{\stackrel{\&OverBar;}{r}}_{\mathrm{ri}}}{r_{\mathrm{ri}}^3}-\frac{{\stackrel{\&OverBar;}{r}}_{\mathrm{pi}}}{r_{\mathrm{pi}}^3}]-\frac{\mathrm{AG}}{{\mathrm{mr}}^3}\stackrel{\&OverBar;}{r}+\frac{k}{m}\stackrel{\&OverBar;}{T}+\stackrel{\&OverBar;}{a}\end{array}\right]$

In the formula: r detector position vector in heliocentric coordinates; V detector velocity in heliocentric coordinates; r _PiThe position vector of i perturbation planetary in heliocentric coordinates, r _RiThe position vector of relative i the perturbation planetary of detector, wherein r _Ri=r _Pi-r; μ _δSolar gravitation constant GM; μ _iThe gravitational constant of i perturbation celestial body; n _pThe quantity of perturbation celestial body; The useful area of A detector; G day luminous flux constant; T propulsion system thrust vectoring; K thrust scale-up factor is approximately 1; M detector quality; The coriolis acceleration that a detector is total;

2) determine observing matrix:

For a certain epoch of observation, observing matrix is

$H=\left[\begin{array}{cccccc}\&PartialD;p/\&PartialD;x_1& \&PartialD;p/\&PartialD;x_2& \&PartialD;p/\&PartialD;x_3& 0& 0& 0\\ \&PartialD;l/\&PartialD;x_1& \&PartialD;l/\&PartialD;x_2& \&PartialD;l/\&PartialD;x_3& 0& 0& 0\end{array}\right]$

In the formula, the partial differential relevant with speed component all is zero;

3) unify epoch of observation:

Utilize state-transition matrix that the observing matrix of each epoch of observation is transformed to same epoch of observation, have

$\tilde{H}=\mathrm{H\&Phi;}$

Wherein

It is exactly the observation partial differential matrix of carving at a time;

4) realize that track is definite:

The recursion weighted least square algorithm that utilization is decomposed based on the UD covariance is determined the track of detector, measurement residual error Y wherein is that the pixel corresponding with calculating asteroid center that prediction obtains that obtain of Flame Image Process is poor, so just utilize many asteroidal observation datas, it is the pixel at asteroid center, revise the orbit parameter of certain epoch of observation by the multistep recursion, at first provide a covariance matrix P ₀, the weighting obsdervations matrix W, vector Y is the difference of observation barycenter and nominal trajectory forecast calculating barycenter, initial point state filtering equations is in batches found the solution vector

With formal covariance matrix P:

$\hat{q}={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}{\tilde{H}}^T\mathrm{WY}$

$P={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}$

Wherein,

$W=\left[\begin{array}{cc}1/{\mathrm{\&sigma;}}_0^2& 0\\ 0& 1/{\mathrm{\&sigma;}}_0^2\end{array}\right]$

${\mathrm{\&sigma;}}_0=\frac{{\tan }^{-1}(R/\mathrm{\&rho;})}{13\&times;{10}^{-6}}$

R asteroid supposition radius, ρ asteroid scope.

Images