CN102162731B - High-precision satellite independent navigation method based on pulse data of sun, earth and moon integrated sensor - Google Patents

Abstract

The invention discloses to a high-precision satellite independent navigation method based on pulse data of a sun, earth and moon integrated sensor, and relates to the field of satellite navigation. The high-precision satellite independent navigation method is provided in order to solve the problem that the photocenter and the centroid of the moon cannot be overlapped because the pulse data comprising noises, the oblateness of the earth, the first quarter and the last quarter of the moon is used as the original data of a navigation system. The method comprises the following steps of: 1, providing the pulse data for a navigation computer; 2, determining directions of the three celestial bodies according to the pulse data; 3, calculating an inner product of a vector of the heliocentric direction and a vector of the geocentric direction, and an inner product of a vector of the selenocenter direction and a vector of the geocentric direction; 4, performing double-vector rough posture fixing; 5, initializing navigation; 6, performing rough navigation calculation so as to obtain a navigation result; 7, finishing refining of the direction of the moon; 8, performing double-vector fine posture fixing; 9, correcting the vectors of a geocentric distance and the geocentric direction; and 10, performing fine navigation calculation so as to obtain a final navigation result. By the high-precision satellite independent navigation method, the influence because the photocenter and the centroid of the moon cannot be overlapped according to the noises, the oblateness of the earth, the first quarter and the last quarter of the moon which are comprised in the pulse data is eliminated.

Description

Satellite high-precision autonomous navigation method based on day ground month integrated sensor pulse data

Technical field

The present invention relates to the satellite navigation field, be specifically related to utilize day that day ground month integrated sensor measures,, month optics pulse data, position, speed and the attitude of satellite are carried out independently definite method.

Background technology

A day ground month autonomous navigation system is made up of navigation sensor and navigational computer.Navigational computer is handled the measured value of navigation sensor, through certain navigation algorithm, makes the position and the speed of satellite in real time, thereby realizes the independent navigation of satellite.Day ground month integrated sensor is as depicted in figs. 1 and 2; Scanning life sensor by the infrared double cone scanning type sensor of the earth and two fan-shaped slit visual fields is formed; The infrared double cone scanning type sensor of the earth has single optical scanning head, utilizes mirror structure to obtain two infrared visual fields, and the track of infrared visual field, scanning back is two coaxial circular cones; Optical head scanning one circle; Pyroelectric detector can detect four Horizons at most and pass through signal, can confirm the orientation of the earth's core direction vector with respect to satellite by the moment that signal occurs, and can try to achieve the distance of satellite to the earth's core.Two visible light sensors on the basis of the infrared double cone scanning type sensor of the earth, have been increased; In the scanning process of optical head; Fan-shaped slit visual field is inswept, and spherical zone is regional; Utilize the si-photodiode detector can be responsive to the sun and the moon, according to the moment that the sun, the moon occur in fan-shaped slit visual field can be in the hope of the orientation of its direction vector with respect to satellite.Detecting device has a plurality of light intensity threshold values, can distinguish the sun and moon signal, and can reject earth signal.Day ground month navigational system is exactly to utilize such sensor to confirm day orientation of the ground moon, thereby carries out the system of autonomous navigation of satellite.

At present to day ground month navigational system the research approach mainly be to be embodied in the angle information of utilizing; Study navigation algorithm; Do not consider measuring principle, these key factors that therefore yet can't embody compression of the earth, moon photocentre and barycenter do not overlap when going up lower edge and engineering is actual differs far away.

Pulse data is as the raw data of navigational system; The inside has comprised all information, and (noise, compression of the earth, last moon at the last quarter ball photocentre and barycenter do not overlap; Or the like), the objective of the invention is to study day measuring principle of ground month integrated sensor, and research is based on the high-precision ground month navigation algorithm of subsisting of pulse data; This algorithm can compensate compression of the earth, also can the moon pulse during the last lower edge be compensated.

Summary of the invention

The present invention has wherein comprised noise, compression of the earth, has gone up the problem that moon at the last quarter ball photocentre and barycenter do not overlap in the inside as the raw data of navigational system in order to solve pulse data in the existing method, and has proposed the satellite high-precision autonomous navigation method based on day ground month integrated sensor pulse data.A measuring principle that the objective of the invention is to day ground month integrated sensor, and based on the high-precision ground month navigation algorithm of subsisting of pulse data, this algorithm can compensate compression of the earth, also can the moon pulse during the last lower edge be compensated.

The step of the satellite high-precision autonomous navigation method of month integrated sensor pulse data is following with the present invention is based on day:

Step 1: the navigation sensor navigation computing machine of being made up of the scanning life sensor of the infrared double cone scanning type sensor of the earth and two fan-shaped slit visual fields provides pulse data;

Step 2: the orientation of carrying out three celestial bodies according to pulse data is confirmed: confirm the coordinate E of the earth's core direction vector under measurement coordinate system _SE, day heart direction vector coordinate S under measurement coordinate system _S|SEWith the coordinate M of moon heart direction vector under measurement coordinate system _S|SE, and the earth's core is apart from r;

Step 3: by the coordinate E of the earth's core direction vector under measurement coordinate system _SE, day heart direction vector coordinate S under measurement coordinate system _S|SEWith the coordinate M of moon heart direction vector under measurement coordinate system _S|SECalculate the inner product of day heart direction vector and the earth's core direction vector, and the inner product of month heart direction vector and the earth's core direction vector:

Step 4: carry out two vectors and slightly decide appearance, the vector of three pairwise orthogonals of definition is following:

v ₁：r _m×r _s/|r _m×r _s|

v ₂：r _m/|r _m| (1)

v ₃：v ₁×v ₂

Wherein, r _mBe star moon vector, r _sBe star day vector; Described star moon vector r _mWith star day vector r _sExpression under inertial coordinates system by the ground moon vector, day vector is approximate obtains, and ground month vector can be obtained by ephemeris with a ground day vector; If v ₁, v ₂, v ₃Under inertial coordinates system, be expressed as E ₁, E ₂, E ₃, under measurement coordinate system, be expressed as e ₁, e ₂, e ₃, then inertial coordinate is tied to the rotation matrix of measurement coordinate system:

R _SEI＝[e ₁，e ₂，e ₃][E ₁，E ₂，E ₃] ^-1 (2)

Step 5: the navigation initialization, computes is used in the initial position guestimate:

r ₀＝-rR _ISEE _SE (3)

is by the formula mistake wherein! Do not find Reference source.Obtain; The initial velocity guestimate then has the position of two adjacent moment to do approximate difference to obtain;

Step 6: adopt the least square navigation algorithm to carry out thick navigation operations and obtain navigation results;

Step 7: utilize the navigation results of the thick navigation in the step 6, revise star moon vector r _mValue, thereby obtain revised star moon vector r _{M repaiies}Accomplished becoming more meticulous of moon orientation;

Step 8: carry out two vectors essences and decide appearance, star moon vector r _mExpression under inertial coordinates system by on the revised star moon vector r that produces of a step _{M repaiies}Substitute star day vector r _sConstant, the vector of three pairwise orthogonals of definition is following:

v ₁': r _{M repaiies}* r _s/ | r _{M repaiies}* r _s|

v ₂': r _{M repaiies}/ | r _{M repaiies}|

v ₃′：v ₁×v ₂

If v ₁', v ₂', v ₃' under inertial system, be expressed as E ₁', E ₂', E ₃', under body coordinate system, be expressed as e ₁', e ₂', e ₃' also can ask; The inertial coordinate that two vectors of then revising are decided appearance is tied to the rotation matrix of measurement coordinate system:

R _SEI′＝[e ₁′，e ₂′，e ₃′][E ₁′，E ₂′，E ₃′] ^-1 (4)

Step 9: consider that inertial coordinate that compression of the earth is decided appearance based on two vectors of revising in the step 8 is tied to the rotation matrix R of measurement coordinate system _SEI' revise the earth's core apart from r and the earth's core orientation vector E _SE, obtain revised the earth's core apart from r _RepairWith revised the earth's core orientation vector E _{SE repaiies}

Step 10: adopt the least square navigation algorithm to carry out smart navigation operations and obtain final navigation results.

The present invention be according to day,, the optical pulse information of month integrated sensor confirms the method for position, speed and the attitude of aircraft.At first be day,, month the orientation confirm, two kinds of methods are arranged for orientation, the earth's core: confirm (being used for thick navigation) and confirm (smart navigation, ellipticity correction) based on the orientation, the earth's core of the spherical hypothesis of the earth based on the orientation, the earth's core of compression of the earth; Moon heart orientation for during the last lower edge is confirmed, has considered that photocentre and moon barycenter do not overlap the influence that brings, and have provided compensation scheme.Provide the least square navigation implementation algorithm of a cover then, carry out the simulation accuracy assessment at last based on compensation of pulse data ellipticity and last lower edge correction.

Sensor scanning rotating speed is ω _Rot=240r/min.Two semi-cone angle of circular cone Horizon sensor infrared horizon visual field circular cone are respectively γ ₁=38 ° and γ ₂=73 °.Scanning head is with respect to slit sensor plane of symmetry M ₁M ₂The hysteresis angle be B _R1=B _R2=0 °, the sensor measurement coordinate system is γ with respect to the setting angle of system _I=π/6, β _I=π/6, α _I=π/6.Fan-shaped life sensor is β with respect to the pitch angle with the scanning rotating shaft _s=16 ° of slit visual fields 1 are ahead of angular distance υ=4 ° of slit visual field 2 in the sensor equatorial plane.

The initial time of track is taken as 0: 0: 0 on the 6th May in 2012.Initial six key elements of true track are: major semi-axis is 6878km, and orbital eccentricity is 0, and orbit inclination is 92 °, and right ascension of ascending node is π/6, and argument of perigee is 0 °, and very near angle is 0 °.When simulated data produced, the attitude of supposing satellite was strict absolute orientation.The time of using the orbit measurement data is in 5000 seconds since 0: 0: 0 on the 12nd May in 2012.

Orbital data with STK produces is come analog pulse; Consider compression of the earth in the pulse generation mechanism; 3 σ of earth impulsive noise=6.5e-5s, 3 σ=0.1 of corresponding angle noise °, the earth sensor systematic error is that (the respective pulses noise is 2.3e-5s to 0.03 degree; 3 σ=the 3.5e-5s of the scanning sun (moon) impulsive noise, 3 σ=0.05 of corresponding angle noise °.Navigation model is still used the dynamics of orbits model among the MATLAB, carries out Monte Carlo simulation, and simulation times is 120 times.Track initial value estimated value and the true track filtering the initial value of the error location mean

location MSE

speed mean

speed standard deviation position error and velocity error 3σ ellipsoid in Figures 6 and 7.

Description of drawings

Fig. 1 is day mounting structure of ground month integrated sensor; A is a spacecraft, and B is the rotating shaft of the infrared double cone scanning type sensor of first earth, and C is the rotating shaft of the infrared double cone scanning type sensor of second earth, and D is the infrared visual fields of 38 degree, and E is the infrared visual fields of 73 degree, and F is the fan-shaped visual field of visible light; Fig. 2 is the visual field of the scanning life sensor of two fan-shaped slit visual fields; L is fan-shaped visual field, and M is the scanning rotating shaft, and N is an infrared probe; Fig. 3 is that earth sensor is measured geometric graph; Fig. 4 is a day ground month orientation synoptic diagram; Fig. 5 is a day ground month orientation synoptic diagram; Fig. 6 is 3 σ ellipsoid synoptic diagram of site error; Fig. 7 is 3 σ ellipsoid synoptic diagram of velocity error.

Embodiment

Embodiment one: combine Fig. 1 to Fig. 5 that this embodiment is described, this embodiment concrete steps are following:

Step 1: the navigation sensor navigation computing machine of being made up of the scanning life sensor of the infrared double cone scanning type sensor of the earth and two fan-shaped slit visual fields provides pulse data;

Step 2: the orientation of carrying out three celestial bodies according to pulse data is confirmed: confirm the coordinate E of the earth's core direction vector under measurement coordinate system _SE, day heart direction vector coordinate S under measurement coordinate system _S|SEWith the coordinate M of moon heart direction vector under measurement coordinate system _S|SE, and the earth's core is apart from r;

Step 3: by the coordinate E of the earth's core direction vector under measurement coordinate system _SE, day heart direction vector coordinate S under measurement coordinate system _S|SEWith the coordinate M of moon heart direction vector under measurement coordinate system _S|SECalculate the inner product of day heart direction vector and the earth's core direction vector, and the inner product of month heart direction vector and the earth's core direction vector:

Step 4: carry out two vectors and slightly decide appearance, the vector of three pairwise orthogonals of definition is following:

v ₁：r _m×r _s/|r _m×r _s|

v ₂：r _m/|r _m| (5)

v ₃：v ₁×v ₂

Wherein, r _mBe star moon vector, r _sBe star day vector; Described star moon vector r _mWith star day vector r _sExpression under inertial coordinates system by the ground moon vector, day vector is approximate obtains, and ground month vector can be obtained by ephemeris with a ground day vector; If v ₁, v ₂, v ₃Under inertial coordinates system, be expressed as E ₁, E ₂, E ₃, under measurement coordinate system, be expressed as e ₁, e ₂, e ₃, then inertial coordinate is tied to the rotation matrix of measurement coordinate system:

R _SEI＝[e ₁，e ₂，e ₃][E ₁，E ₂，E ₃] ^-1 (6)

Step 5: the navigation initialization, computes is used in the initial position guestimate:

r ₀＝-rR _ISEE _SE (7)

is by the formula mistake wherein! Do not find Reference source.Obtain; The initial velocity guestimate then has the position of two adjacent moment to do approximate difference to obtain;

Step 6: adopt the least square navigation algorithm to carry out thick navigation operations and obtain navigation results;

Step 7: utilize the navigation results of the thick navigation in the step 6, revise star moon vector r _mValue, thereby obtain revised star moon vector r _{M repaiies}Accomplished becoming more meticulous of moon orientation;

Step 8: carry out two vectors essences and decide appearance, star moon vector r _mExpression under inertial coordinates system by on the revised star moon vector r that produces of a step _{M repaiies}Substitute star day vector r _sConstant, the vector of three pairwise orthogonals of definition is following:

v ₁': r _{M repaiies}* r _s/ | r _{M repaiies}* r _s|

v ₂': r _{M repaiies}/ | r _{M repaiies}|

v ₃′：v ₁×v ₂

If v ₁', v ₂', v ₃' under inertial system, be expressed as E ₁', E ₂', E ₃', under body coordinate system, be expressed as e ₁', e ₂', e ₃' also can ask; The inertial coordinate that two vectors of then revising are decided appearance is tied to the rotation matrix of measurement coordinate system:

R _SEI′＝[e ₁′，e ₂′，e ₃′][E ₁′，E ₂′，E ₃′] ^-1 (8)

Step 9: consider that inertial coordinate that compression of the earth is decided appearance based on two vectors of revising in the step 8 is tied to the rotation matrix R of measurement coordinate system _SEI' revise the earth's core apart from r and the earth's core orientation vector E _SE, obtain revised the earth's core apart from r _RepairWith revised the earth's core orientation vector E _{SE repaiies}

Step 10: adopt the least square navigation algorithm to carry out smart navigation operations and obtain final navigation results.

Embodiment two: this embodiment is the coordinate E of the earth's core direction vector under measurement coordinate system in the step 2 with embodiment one difference _SEDefinite process, and the earth's core is following apart from definite process of r:

Two infrared probes on the infrared double cone scanning type sensor of the earth sweep terrestrial time:

α _1-in＝-((ω _rott _ref+B _R1)-ω _rott _1-in)

(9)

α _2-in＝-((ω _rott _ref+B _R2)-ω _rott _2-in)

ω wherein _RotBe the rotating speed of scanning sensor, B _R1, B _R2Be the drag angle of two infrared probes on the initial time double cone scanning type sensor with respect to benchmark x axle, t _RefFor the slit sensor plane of symmetry passes benchmark x axle to deserved pulse constantly, t _1-in, t _2-inThe pulse that forms when sweeping Horizon for two infrared probes on the double cone scanning type sensor;

Two infrared probes on the infrared double cone scanning type sensor of the earth scan out terrestrial time:

α _1-out＝ω _rott _1-out-(ω _rott _ref+B _R1)

(10)

α _2-out＝ω _rott _2-out-(ω _rott _ref+B _R2)

T wherein _1-out, t _2-outThe pulse that forms when scanning out Horizon for two infrared probes on the double cone scanning type sensor;

So the earth's core direction vector with respect to the position angle of measurement coordinate system is:

${\mathrm{\φ}}_{e1}=-(\frac{1}{2}(-{\mathrm{\α}}_{1-\mathrm{in}}+{\mathrm{\α}}_{1-\mathrm{out}})-{\mathrm{\α}}_{1-\mathrm{out}})$

(11)

${\mathrm{\φ}}_{e2}=-(\frac{1}{2}(-{\mathrm{\α}}_{2-\mathrm{in}}+{\mathrm{\α}}_{2-\mathrm{out}})-{\mathrm{\α}}_{2-\mathrm{out}})$

If do not consider compression of the earth, φ arranged _E1=φ _E2, and use φ _eExpression;

Thereby the string of the earth that the infrared double cone scanning type sensor of the earth obtains is wide be:

Ω _E1＝-α _1-in+α _1-out

(12)

Ω _E2＝-α _2-in+α _2-out

By the spherical trigonometry formula, have

$\cos \mathrm{\&rho;}=\cos {\mathrm{\&gamma;}}_1\cos \mathrm{\&eta;}+\sin {\mathrm{\&gamma;}}_1\sin \mathrm{\&eta;}\cos \frac{{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="1"\; label="E"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">E</figure-callout>}1}}{2}$

$\cos \mathrm{\&rho;}=\cos {\mathrm{\&gamma;}}_2\cos \mathrm{\&eta;}+\sin {\mathrm{\&gamma;}}_2\sin \mathrm{\&eta;}\cos \frac{{\mathrm{\&Omega;}}_{E2}}{2}$

Wherein ρ is the half angle with respect to the infrared radiation disk of the earth of satellite, γ ₁, γ ₂Be the semi-cone angle of earth sensor double cone, η is the angle between the earth's core direction and the sensor turning axle;

So

$\mathrm{\&eta;}={\cot }^{-1}\frac{{\sin \mathrm{\&gamma;}}_2\cos 0.5{\mathrm{\&Omega;}}_{E2}-{\sin \mathrm{\&gamma;}}_1\cos 0.5{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="1"\; label="E"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">E</figure-callout>}1}}{\cos {\mathrm{\&gamma;}}_1-{\cos \mathrm{\&gamma;}}_2}---\left(13\right)$

If do not consider compression of the earth, the coordinate of the earth's core direction vector under measurement coordinate system try to achieve into:

$E_{\mathrm{SE}}=\left[\begin{array}{c}\sin \mathrm{\&eta;}{\cos \mathrm{\&phi;}}_e\\ \sin \mathrm{\&eta;}\sin {\mathrm{\&phi;}}_e\\ \cos \mathrm{\&eta;}\end{array}\right]---\left(14\right)$

The earth's core is apart from by computes:

r＝R _E/sinρ (15)

Wherein, R _EBe the infrared radiation radius of a ball of the earth, ρ can be tried to achieve by the cosine law of spherical triangle:

$\cos \mathrm{\&rho;}=\cos \mathrm{\&gamma;}\cos \mathrm{\&eta;}+\sin \mathrm{\&gamma;}\sin \mathrm{\&eta;}\cos \frac{{\mathrm{\&Omega;}}_E}{2}---\left(16\right)$

Other step is identical with embodiment one.

Embodiment three: this embodiment is the coordinate S of day heart direction vector under measurement coordinate system in the step 2 with embodiment one difference _S|SEDefinite process following:

Position angle and the elevation angle of day heart orientation under the measurement body coordinate system is:

${\mathrm{\&phi;}}_s={\mathrm{\&omega;}}_{\mathrm{rot}}(\frac{{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">t</figure-callout>}}_{1\mathrm{sum}}+t_{2\mathrm{sum}}}{2}-t_{\mathrm{ref}})---\left(17\right)$

δ _s＝cot ^-1(sinσ _scosβ _s)

Wherein: t _1sun, t _2sunBe the responsive pulse that forms to sunshine time the in the first fan-shaped slit visual field, the second fan-shaped slit visual field, σ _s=(ω _Rott _2sun-ω _Rott _1sun-υ)/2, υ is the scanning life sensor of the first fan-shaped slit visual field is ahead of the scanning life sensor of the second fan-shaped slit visual field on the sensor equatorial plane a focal length, β _sBe the angle of slit visual field with respect to the scanning axes of rotation skew;

The coordinate of day heart direction vector under measurement coordinate system is:

$S_{S|\mathrm{SE}}=\left[\begin{array}{c}{\cos \mathrm{\&delta;}}_s{\cos \mathrm{\&phi;}}_s\\ \cos {\mathrm{\&delta;}}_s\sin {\mathrm{\&phi;}}_s\\ \sin {\mathrm{\&delta;}}_s\end{array}\right]---\left(18\right)$

Other step is identical with embodiment one.

Embodiment four: this embodiment is the coordinate M of moon heart direction vector under measurement coordinate system in the step 2 with embodiment one difference _S|SEDefinite process following:

Similar in the time of formula during the full moon the responsive moon time and responsive sunshine, position angle and the elevation angle of moon heart orientation under measurement coordinate system during the full moon is:

${\mathrm{\&phi;}}_m={\mathrm{\&omega;}}_{\mathrm{rot}}(\frac{{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">t</figure-callout>}}_{1\mathrm{moon}}+t_{2\mathrm{moon}}}{2}-t_{\mathrm{ref}})---\left(19\right)$

δ _m＝tan ^-1(sinσ _Mcosβ _s)

Wherein: t _1moon, t _2moonBe the responsive pulse that forms to the moon time in the first fan-shaped slit visual field, the second fan-shaped slit visual field, σ _M=(ω _Rott _2moon-ω _Rott _1moon-υ)/2;

Thereby the coordinate of moon heart direction vector under measurement coordinate system during the full moon is:

$M_{S|\mathrm{SE}}=\left[\begin{array}{c}{\cos \mathrm{\&delta;}}_m{\cos \mathrm{\&phi;}}_m\\ \cos {\mathrm{\&delta;}}_m\sin {\mathrm{\&phi;}}_m\\ \sin {\mathrm{\&delta;}}_m\end{array}\right]---\left(20\right)$

Moon heart orientation during the last lower edge confirm with full moon during a month heart orientation confirm that formula is the same, but individual impulse compensation module is arranged; Obtain the attitude information of satellite by thick navigation, can obtain sensor scan axis and x _mThe angle theta of (on the normal society face vertical axle) with ground moon direction _x, add a Δ t on the basis of the pulse that the pulse after the compensation should be original

$\mathrm{\&Delta;t}=\left\{\begin{array}{c}\frac{2\mathrm{\&delta;}{t}_{\max }}{\mathrm{\&pi;}}{\mathrm{\&theta;}}_x({\mathrm{\&theta;}}_x<\mathrm{\&pi;}/2)\\ \frac{2\mathrm{\&delta;}{t}_{\max }}{\mathrm{\&pi;}}(\mathrm{\&pi;}-{\mathrm{\&theta;}}_x){\mathrm{\&theta;}}_x\&GreaterEqual;\mathrm{\&pi;}/2\end{array}\right.$

Wherein, Δ t＞0 in the time of the moon at the first quarter, Δ t＜0 in the time of the moon at the last quarter; Again by the formula mistake! Do not find Reference source.Calculate φ _m, δ _m, by the formula mistake! Do not find Reference source.Moon heart orientation during the last lower edge after confirming to compensate.Other step is identical with embodiment one.

Embodiment five: this embodiment is step 7 correction star moon vector r with embodiment one difference _m, because the distance of the sun and the earth is far, so star day vector r _sWith the direction vector R that points to the sun from the earth's core _sCan regard as parallel.And the moon is not far with respect to the distance of the earth, star moon vector r _mWith ground day vector R _mCan not regard as parallel.So star moon vector r _mUnder inertial coordinates system expression month vector R practicably not _mReplace, need carry out the limited distance correction.

The preliminary definite position r of satellite under inertial coordinates system from thick navigation _b, and moon inertial position is tried to achieve by ephemeris, thus confirm star moon vector r _mExpression under inertial coordinates system, correction formula is following:

r _{M repaiies}=R _m-r _b

Embodiment six: combine Fig. 3 that this embodiment is described, this embodiment and embodiment one difference are in the step 9 to confirm to utilize based on the orientation, the earth's core of ellipticity correction the attitude information of step 8, to the earth's core apart from r and the earth's core orientation vector E _SERecomputate specific as follows:

By measuring geometric analysis, can know that the earth's core arrives the vector R of earth surface _EFor:

$R_E={R_{\mathrm{ISE}}}^{\&prime;}(r\hat{r}+l\widehat{\mathrm{\&rho;}})---\left(21\right)$

Wherein:

Be the direction vector (expression under measurement coordinate system) of sensing satellite from the earth's core,

With the unit vector (being projected under the measurement coordinate system) that scans sight line, make pulse t constantly respectively _1in, t _1out, t _2in, t _2outThe unit vector of corresponding scanning sight line does

If do not consider to measure noise, this tittle is accurately known, is respectively:

${\widehat{\mathrm{\&rho;}}}^1=\left[\begin{array}{c}\sin \left({\mathrm{\&gamma;}}_1\right)\cos ({\mathrm{\&omega;}}_{\mathrm{rot}}({\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">t</figure-callout>}}_{1\mathrm{in}}-t_{\mathrm{ref}})-B_{R1})\\ \sin \left({\mathrm{\&gamma;}}_1\right)\sin ({\mathrm{\&omega;}}_{\mathrm{rot}}({\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">t</figure-callout>}}_{1\mathrm{in}}-t_{\mathrm{ref}})-B_{R1})\\ \cos \left({\mathrm{\&gamma;}}_1\right)\end{array}\right]$

${\widehat{\mathrm{\&rho;}}}^2=\left[\begin{array}{c}\sin \left({\mathrm{\&gamma;}}_1\right)\cos ({\mathrm{\&omega;}}_{\mathrm{rot}}({\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">t</figure-callout>}}_{1\mathrm{out}}-t_{\mathrm{ref}})-B_{R1})\\ \sin \left({\mathrm{\&gamma;}}_1\right)\sin ({\mathrm{\&omega;}}_{\mathrm{rot}}({\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">t</figure-callout>}}_{1\mathrm{out}}-t_{\mathrm{ref}})-B_{R1})\\ \cos \left({\mathrm{\&gamma;}}_1\right)\end{array}\right]$

(22)

${\widehat{\mathrm{\&rho;}}}^3=\left[\begin{array}{c}\sin \left({\mathrm{\&gamma;}}_2\right)\cos ({\mathrm{\&omega;}}_{\mathrm{rot}}(t_{2\mathrm{in}}-t_{\mathrm{ref}})-B_{R2})\\ \sin \left({\mathrm{\&gamma;}}_2\right)\sin ({\mathrm{\&omega;}}_{\mathrm{rot}}(t_{2\mathrm{in}}-t_{\mathrm{ref}})-B_{R2})\\ \cos \left({\mathrm{\&gamma;}}_2\right)\end{array}\right]$

${\widehat{\mathrm{\&rho;}}}^4=\left[\begin{array}{c}\sin \left({\mathrm{\&gamma;}}_2\right)\cos ({\mathrm{\&omega;}}_{\mathrm{rot}}(t_{2\mathrm{out}}-t_{\mathrm{ref}})-B_{R2})\\ \sin \left({\mathrm{\&gamma;}}_2\right)\sin ({\mathrm{\&omega;}}_{\mathrm{rot}}(t_{2\mathrm{out}}-t_{\mathrm{ref}})-B_{R2})\\ \cos \left({\mathrm{\&gamma;}}_2\right)\end{array}\right]$

The earth's core is projected under the inertial coordinates system to the vector

on earth infrared radiation ball surface.

The ellipsoid equation is:

P=-1+1/ (1-ee) wherein ².

By getting in the formula (24):

(25)

${z_e}^2={\left(L_3\widehat{\mathrm{\&rho;}}\right)}^2{l}^2+2r{L}_3\hat{r}{L}_3\widehat{\mathrm{\&rho;}}l+r^2{\left(L_3\hat{r}\right)}^2$

Wherein: L ₃Be R _ISEThe third line.

With mistake! Do not find Reference source.A substitution mistake! Do not find Reference source.To obtaining quadratic equation:

(a+Δa)l ²+(b+Δb)l+(c+Δc)＝0 (26)

Wherein:

$\mathrm{\&Delta;a}=p{\left(L_3\widehat{\mathrm{\&rho;}}\right)}^2,\mathrm{\&Delta;b}=2\mathrm{pr}{L}_3\hat{r}{L}_3\widehat{\mathrm{\&rho;}},\mathrm{\&Delta;c}={\mathrm{pr}}^2{\left(L_3\hat{r}\right)}^2$

The corresponding scanning direction of visual lines vector

constantly of note pulse then has:

(b+Δb) ²-4(a+Δa)(c+Δc)＝0 (27)

Also promptly:

(28)

Each constantly in the navigation slightly can be solved from last equation as initial value revised

it should be noted that

also should satisfy constraint condition

4 unknown numbers, 5 equations approach with least-squares estimation Order

x ₀Be the iterative initial value of x, with F at x ₀Point single order Taylor expansion has:

$0=F\left(x\right)=F\left(x_0\right)+\frac{\&PartialD;F}{\&PartialD;x}{|}_{x=x_0}(x-x_0)$

Wherein: $\frac{\&PartialD;F}{\&PartialD;x}=\left[\begin{array}{cc}\frac{{\&PartialD;\mathrm{<figure-callout\; id="1"\; label="f"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">f</figure-callout>}}_1}{\&PartialD;r}& \frac{\&PartialD;{\mathrm{<figure-callout\; id="1"\; label="f"\; filenames="HDA0000043528470000031.png"\; state="\{\{state\}\}">f</figure-callout>}}_1}{\&PartialD;\hat{r}}\\ \frac{\&PartialD;f_2}{\&PartialD;r}& \frac{\&PartialD;f_2}{\&PartialD;\hat{r}}\\ \frac{\&PartialD;f_3}{\&PartialD;r}& \frac{\&PartialD;f_3}{\&PartialD;\hat{r}}\\ \frac{\&PartialD;f_4}{\&PartialD;r}& \frac{{\&PartialD;f}_4}{\&PartialD;\hat{r}}\\ \frac{\&PartialD;f_5}{\&PartialD;r}& \frac{\&PartialD;f_5}{\&PartialD;\hat{r}}\end{array}\right],$

And:

$\frac{{\&PartialD;f}_i}{\&PartialD;r}=2\frac{R_E^2}{r^3}(1+p{\left(L_3{\widehat{\mathrm{\&rho;}}}^i\right)}^2)i=\mathrm{1,2},\mathrm{3,4}$

Because

gets dimension is 5 * 4, makes

then to have:

The form of being write as iteration has:

x _kBeing exactly the optimal estimation of x, also is revised the earth's core distance and the earth's core direction vector

Embodiment seven: this embodiment is that with embodiment one difference the least square navigation algorithm in step 6 and the step 10 is: suppose N observation t constantly ₁, t ₂..., t _NObserved quantity do Described observed quantity is angle and the angle between day ground between month ground, is obtained by the arc cosine of step 3 inner product, and its true value does

Make residual error be:

$G({\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0^{*},t_i)=\hat{H}({\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0^{*},t_i)-H({\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0^{*},t_i)$

Not it should be noted that is unknown if measure noise , because true value

is unknown.

Set up the least square index:

Wherein:

$G({\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0,t_i)=\hat{H}({\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0^{*},t_i)-H({\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0,t_i)$

Obviously the time as

; J has minimal value, then

$\frac{\&PartialD;J}{{\&PartialD;\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0}{|}_{X_0={\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0^{*}}=0$

If

For

Estimated value, then with G (X ₀, t _i)

Near linearization can get:

$G({\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0,t_i)=G({\tilde{X}}_0,t_i)+\frac{\&PartialD;G}{{\&PartialD;\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0}{|}_{X_0={\tilde{X}}_0}(X_0-{\tilde{X}}_0)$

To go up in the equation substitution functional index and can get:

Following formula can be remembered and does:

$A({\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0^{*}-{\tilde{X}}_0)+B=0$

Wherein:

Thereby can in the hope of:

${\mathrm{<figure-callout\; id="0"\; label="X"\; filenames="HDA0000043528470000023.png,HDA0000043528470000031.png"\; state="\{\{state\}\}">X</figure-callout>}}_0^{*}={\tilde{X}}_0-A^{-1}B$

Corresponding iterative formula can be write:

X _k+1＝X _k-A ^-1B

Content of the present invention is not limited only to the content of above-mentioned each embodiment, and the combination of one of them or several embodiments equally also can realize the purpose of inventing.

Claims (4)

1. based on the satellite high-precision autonomous navigation method of day ground month integrated sensor pulse data, it is characterized in that its step is following:

Step 1: the navigation sensor navigation computing machine of being made up of the scanning life sensor of the infrared double cone scanning type sensor of the earth and two fan-shaped slit visual fields provides pulse data;

Step 2: the orientation of carrying out three celestial bodies according to pulse data is confirmed: confirm the coordinate E of the earth's core direction vector under measurement coordinate system _SE, day heart direction vector coordinate S under measurement coordinate system _S|SEWith the coordinate M of moon heart direction vector under measurement coordinate system _S|SE, and the earth's core is apart from r;

Step 3: by the coordinate E of the earth's core direction vector under measurement coordinate system _SE, day heart direction vector coordinate S under measurement coordinate system _S|SEWith the coordinate M of moon heart direction vector under measurement coordinate system _S|SECalculate the inner product of day heart direction vector and the earth's core direction vector, and the inner product of month heart direction vector and the earth's core direction vector:

Step 4: carry out two vectors and slightly decide appearance, the vector of three pairwise orthogonals of definition is following:

v ₁：r _m×r _s/|r _m×r _s|

v ₂：r _m/|r _m| (1)

v ₃：v ₁×v ₂

Wherein, r _mBe star moon vector, r _sBe star day vector; Described star moon vector r _mWith star day vector r _sExpression under inertial coordinates system by the ground moon vector, day vector is approximate obtains, and ground month vector can be obtained by ephemeris with a ground day vector; If v ₁, v ₂, v ₃Under inertial coordinates system, be expressed as E ₁, E ₂, E ₃, under measurement coordinate system, be expressed as e ₁, e ₂, e ₃, then inertial coordinate is tied to the rotation matrix of measurement coordinate system:

R _SEI＝[e ₁，e ₂，e ₃][E ₁，E ₂，E ₃] ^-1 (2)

Step 5: the navigation initialization, computes is used in the initial position guestimate:

r ₀＝-rR _ISEE _SE (3)

Wherein

obtained by formula (2); The initial velocity guestimate then has the position of two adjacent moment to do approximate difference to obtain;

Step 6: adopt the least square navigation algorithm to carry out thick navigation operations and obtain navigation results;

Step 7: utilize the navigation results of the thick navigation in the step 6, revise star moon vector r _mValue, thereby obtain revised star moon vector r _{M repaiies}Accomplished becoming more meticulous of moon orientation;

Step 8: carry out two vectors essences and decide appearance, star moon vector r _mExpression under inertial coordinates system by on the revised star moon vector r that produces of a step _{M repaiies}Substitute star day vector r _sConstant, the vector of three pairwise orthogonals of definition is following:

v ₁': r _{M repaiies}* r _s/ | r _{M repaiies}* r _s|

v ₂': r _{M repaiies}/ | r _{M repaiies}|

v ₃′：v ₁×v ₂

If v ₁', v ₂', v ₃' under inertial system, be expressed as E ₁', E ₂', E ₃', under body coordinate system, be expressed as e ₁', e ₂', e ₃' also can ask; The inertial coordinate that two vectors of then revising are decided appearance is tied to the rotation matrix of measurement coordinate system:

R _SEL′＝[e ₁′，e ₂′，e ₃′][E ₁′，E ₂′，E ₃′] ^-1 (4)

Step 9: consider that inertial coordinate that compression of the earth is decided appearance based on two vectors of revising in the step 8 is tied to the rotation matrix R of measurement coordinate system _SEI' revise the earth's core apart from r and the earth's core orientation vector E _SE, obtain revised the earth's core apart from r _RepairWith revised the earth's core orientation vector E _{SE repaiies}

Step 10: adopt the least square navigation algorithm to carry out smart navigation operations and obtain final navigation results.

2. the satellite high-precision autonomous navigation method based on day ground month integrated sensor pulse data according to claim 1 is characterized in that the coordinate E of the earth's core direction vector under measurement coordinate system in the step 2 _SEDefinite process, and the earth's core is following apart from definite process of r:

Two infrared probes on the infrared double cone scanning type sensor of the earth sweep terrestrial time:

α _1-in＝-((ω _rott _ref+B _R1)-ω _rott _1-in)

(5)

α _2-in＝-((ω _rott _ref+B _R2)-ω _rott _2-in)

ω wherein _RotBe the rotating speed of navigation sensor, B _R1, B _R2Be the drag angle of two infrared probes on the infrared double cone scanning type sensor of the initial time earth with respect to benchmark x axle, t _RefThe scanning life sensor plane of symmetry that is two fan-shaped slit visual fields passes benchmark x axle to deserved pulse constantly, t _1-in, t _2-inThe pulse that forms when sweeping Horizon for two infrared probes on the infrared double cone scanning type sensor of the earth;

Two infrared probes on the infrared double cone scanning type sensor of the earth scan out terrestrial time:

α _1-out＝ω _rott _1-out-(ω _rott _ref+B _R1)

(6)

α _2-out＝ω _rott _2-out-(ω _rott _ref+B _R2)

T wherein _1-out, t _2-outThe pulse that forms when scanning out Horizon for two infrared probes on the infrared double cone scanning type sensor of the earth;

So the earth's core direction vector with respect to the position angle of measurement coordinate system is:

(7)

If do not consider compression of the earth, φ arranged _E1=φ _E2, and use φ _eExpression;

Thereby the string of the earth that the infrared double cone scanning type sensor of the earth obtains is wide be:

Ω _E1＝-α _1-in+α _1-out

(8)

Ω _E2＝-α _2-in+α _2-out

By the spherical trigonometry formula, have

Wherein ρ is the half angle with respect to the infrared radiation disk of the earth of satellite, γ ₁, γ ₂Be the semi-cone angle of the infrared double cone scanning type sensor of earth double cone, η is the angle between the earth's core direction and the infrared double cone scanning type sensor of the earth turning axle;

So

If do not consider compression of the earth, the coordinate of the earth's core direction vector under measurement coordinate system try to achieve into:

The earth's core is apart from by computes:

r＝R _E/sinρ (11)

Wherein, R _EBe the infrared radiation radius of a ball of the earth, ρ can be tried to achieve by the cosine law of spherical triangle:

Wherein, γ is the semi-cone angle of the infrared double cone scanning type sensor of earth double cone.

3. the satellite high-precision autonomous navigation method based on day ground month integrated sensor pulse data according to claim 1 is characterized in that step 7 correction star moon vector r _m, the preliminary definite position r of satellite under inertial coordinates system from thick navigation _b, and moon inertial position is tried to achieve by ephemeris, thus confirm star moon vector r _mExpression under inertial coordinates system, correction formula is following:

r _mRepair=R _m-r _b

4. the satellite high-precision autonomous navigation method based on day ground month integrated sensor pulse data according to claim 1 is characterized in that the least square navigation algorithm in step 6 and the step 10 is: suppose N observation t constantly ₁, t ₂..., t _NObserved quantity do

Described observed quantity is angle and the angle between day ground between month ground, is obtained by the arc cosine of step 3 inner product, and its true value does

Make residual error be:

If do not measure noise

Set up the least square index:

Wherein:

Obviously the time as

; J has minimal value, then

If For

Estimated value, then with G (X ₀, t _i) Near linearization can get:

To go up in the equation substitution functional index and can get:

Following formula can be remembered and does:

Wherein:

Thereby can be in the hope of B

Corresponding iterative formula can be write:

X _k+1＝X _k-A ^-1B。

Images