CN100575877C - Spacecraft shading device combined navigation methods based on many information fusion - Google Patents

Spacecraft shading device combined navigation methods based on many information fusion Download PDF

Info

Publication number
CN100575877C
CN100575877C CN200710191527A CN200710191527A CN100575877C CN 100575877 C CN100575877 C CN 100575877C CN 200710191527 A CN200710191527 A CN 200710191527A CN 200710191527 A CN200710191527 A CN 200710191527A CN 100575877 C CN100575877 C CN 100575877C
Authority
CN
China
Prior art keywords
satellite
earth
navigation
subsystem
landform
Prior art date
Application number
CN200710191527A
Other languages
Chinese (zh)
Other versions
CN101178312A (en
Inventor
乔黎
刘建业
熊智
赵伟
郑广楼
郁丰
李丹
张丽敏
王融
Original Assignee
南京航空航天大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 南京航空航天大学 filed Critical 南京航空航天大学
Priority to CN200710191527A priority Critical patent/CN100575877C/en
Publication of CN101178312A publication Critical patent/CN101178312A/en
Application granted granted Critical
Publication of CN100575877C publication Critical patent/CN100575877C/en

Links

Abstract

A kind of spacecraft shading device combined navigation methods based on many information fusion belongs to the spacecraft autonomous navigation method.This air navigation aid may further comprise the steps: set up the satellite motion state equation based on dynamics of orbits, foundation is based on the autonomous navigation of satellite subsystem of X-ray detector, foundation is based on the autonomous navigation of satellite subsystem of infrared horizon and star sensor, foundation is based on the autonomous navigation of satellite subsystem of radar altimeter, foundation is based on the autonomous navigation of satellite subsystem of three sensors of ultraviolet, information fusion is carried out in the output of above-mentioned four subsystems handle, the output navigation information adopts x 2Method of inspection detects system and the isolated fault subsystem that breaks down.This Combinated navigation method can be realized the high-precision independent navigation of deep space probe, the reliability height, and fault-tolerant ability is strong.

Description

Spacecraft shading device combined navigation methods based on many information fusion
Technical field
The invention belongs to the spacecraft shading device combined navigation technical field, be applicable to the high-precision independent navigation of the spacecraft that orbits the earth, also can be used for deep space probe, and the navigation of the high-precision independent of interplanetary flight spacecraft.
Background technology
The spacecraft autonomous navigation system can rely on even finally not rely under the situation that ground system supports few, determines in real time the position and the speed of spacecraft to realize autonomous operation at rail.For satellite system, independent navigation helps reducing the degree of dependence of satellite to ground, improves the system survival ability, for example wartime, when ground system suffers enemy's destruction and disturbs, still can finish determining and maintenance of track, this has very important significance to military satellite.In addition, independent navigation can also alleviate the burden of ground control station effectively, reduces ground and supports cost, thereby reduce the development cost of whole space program.
The autonomous navigation method of spacecraft mainly comprises GPS, inter-satellite link, magnetometer, radar and astronomical navigation method at present.Though wherein GPS and inter-satellite link can carry out real-time navigation and reach higher precision, spacecraft must rely on other spacecrafts measurement information is provided, can not complete at last autonomous air navigation aid.Use radar altimeter to carry out the independent navigation of satellite, be subjected to the influence that Terrain Elevation changes, can obtain navigation performance preferably in the time of still above the ocean.
Celestial navigation is a kind of autonomous positioning navigation method fully that celestial body (as planet, fixed star, the X ray pulse magnitude) information of utilizing heavenly body sensor to record is carried out the spacecraft location compute, and attitude, position, speed and temporal information can be provided simultaneously; Be applicable to low orbit satellite, high rail satellite and deep space probe, thereby enjoy favor.In the celestial navigation of satellite, according to the difference of object being observed attribute, the sensor of use mainly contains star sensor, infrared horizon, three sensors of ultraviolet and X-ray detector etc. at present.Wherein the object of observation of X-ray detector is the impulse radiation of x-ray source, obtains pulse arrival time information; Use X-ray detector to provide precision higher navigation information as satellite.For the combination of multiple uranometry sensor, because the difference of its object being observed attribute, thereby there is complementary possibility in its navigation attribute.
Along with the development of electronic technology, airmanship and control theory, have numerous navigation sensor modules on the satellite, possessed the possible condition that constitutes the multi-sensor combined navigation infosystem.Along with the requirements at the higher level that the development of Aero-Space cause, particularly modernized war propose the precision of navigator fix and reliability, single navigational system has been difficult to satisfy.Two or more navigational system are combined, application optimal estimation theory, form the optimum combination navigational system, help fully using the information of each navigational system to carry out message complementary sense and Cooperative For Information, the navigational system after the combination is all being increased aspect precision and the reliability.Automatic detection and trouble isolation serviceability when this also breaks down to satellite combined guidance system are simultaneously had higher requirement.
Summary of the invention
The objective of the invention is to: utilize a plurality of sensors to make up, the measurement data of each sensor is carried out information fusion to be handled, a kind of method that adopts the information fusion means that the high-precision independent navigational parameter is provided for satellite is provided, and automatic detectability and the trouble isolation serviceability of raising satellite when sensor breaks down, the method for reconstruct navigational system.
The technical solution adopted in the present invention is: based on the satellite combined guidance system of many information fusion, utilize the dynamics of orbits model of X-ray detector, star sensor, infrared horizon, radar altimeter and five kinds of sensors of three sensors of ultraviolet and satellite to constitute four autonomous navigation of satellite subsystems, satellite position, speed parameter to four subsystem outputs carry out Data Fusion, obtain the position and the speed parameter optimal estimation value of satellite, satellite is navigated according to described satellite position and the speed parameter that obtains by information fusion.Under the situation to subsystem fault, design error failure detects and partition method, isolated fault subsystem, the influence of fixing a breakdown.
Wherein four autonomous navigation of satellite subsystems are based on the autonomous navigation of satellite subsystem of X-ray detector respectively, autonomous navigation of satellite subsystem based on infrared horizon and star sensor, autonomous navigation of satellite subsystem based on radar altimeter, autonomous navigation of satellite subsystem based on three sensors of ultraviolet, be called subsystem 1, subsystem 2, subsystem 3, subsystem 4.
Specifically may further comprise the steps:
(1) foundation is based on the satellite motion state equation of dynamics of orbits;
(2) foundation is based on the autonomous navigation of satellite subsystem of X-ray detector;
(3) foundation is based on the autonomous navigation of satellite subsystem of infrared horizon and star sensor;
(4) foundation is based on the autonomous navigation of satellite subsystem of radar altimeter;
(5) foundation is based on the autonomous navigation of satellite subsystem of three sensors of ultraviolet;
(6) information fusion is carried out in the output of above-mentioned four subsystems and handle, the output navigation information;
(7) adopt χ 2Method of inspection detects system and the isolated fault subsystem that breaks down.
In the described step (1), foundation has comprised J 2Item non-spherical earth perturbation and the satellite orbit kinetics equation of solar-lunar perturbating, under J2000.0 Earth central inertial system, the satellite orbit kinetics equation is
X · = f ( X , t ) + W ( t ) - - - ( 1 )
In the formula, X=[x, y, z, v x, v y, v z] TBe state variable, x, y, z, v x, v y, v zBe respectively satellite at X, Y, the position and the velocity amplitude of three directions of Z; F (X, form t) is according to the difference of the set perturbation of satellite, only considers under the trisome gravitation situation of the non-spherical earth perturbation of J2 item and day, month,
f ( X , t ) = v x ; v y ; v z - μ e x r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 1.5 ) ] + μ s ( x s - x | r s - r | - x s | r s | ) + μ m ( x m - x | r m - r | - x m | r m | ) - μ e y r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 1.5 ) ] + μ s ( y s - y | r s - r | - y s | r s | ) + μ m ( y m - y | r m - r | - y m | r m | ) - μ e z r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 4.5 ) ] + μ s ( z s - z | r s - r | - z s | r s | ) + μ m ( z m - z | r m - r | - z m | r m | ) - - - ( 2 )
Wherein, μ e, μ s, μ mBe respectively the gravitational constant of the earth, the sun, the moon, R eBe earth radius, J 2Be terrestrial gravitation coefficient, r s=[x s, y s, z s] TAnd r m=[x m, y m, z m] TRepresent the position of the Sun and the Moon under Earth central inertial system respectively, x s, y s, z sBe the coordinate figure of the sun under Earth central inertial system, x m, y m, z mBe the coordinate figure of the moon under Earth central inertial system, x, y, z, v x, v y, v zBe respectively satellite at X, Y, the position and the velocity amplitude of three directions of Z, r is the distance that satellite is arrived in the earth's core, W (t) is the model error vector, has represented the influence of the perturbative force such as higher order term, solar radiation pressure perturbation atmospherical drag perturbation of perturbation of earths gravitational field, is assumed to be the zero-mean white Gaussian noise;
Described step (2) comprises the quality factor computing method of X ray pulsar, the computing method that concern between satellite position and the X ray pulsar rotation phase place, integer ambiguity computing method and GDOP value calculating method;
The computing method of the quality factor of the X ray pulsar in the described step (2) are
Q = 754.8 SNR T = 1 5 T 10 + 7 T 50 5 T 10 - T 50 T 50 ( T 10 - T 50 ) - - - ( 3 )
In the formula, SNR is the signal to noise ratio (S/N ratio) of pulsar, and T is the pulsar swing circle, T 50And T 10Be respectively the pulse width of pulsating wave peak intensity 50% and at 10% o'clock;
The computing method that concern between satellite position in the described step (2) and the X ray pulsar rotation phase place are
[ Δφ + ΔN ] · λ - n ^ · r E = n ^ · r + c δt sat + Δ rel + v k - - - ( 4 )
Wherein, Δ φ=φ SSBSat, Δ N is the integer ambiguity value between satellite and the solar system barycenter SSB, λ is the wavelength of the impulse radiation of pulsar, Be the direction of visual lines of pulsar, r EFor the earth relatively and the position vector of SSB, r is the position vector of satellite with respect to the earth's core, c is the light velocity, δ t SatBe the clock correction on the satellite, v kFor impulse phase is measured noise, Δ RelBe the relativistic effect correction member, comprise that mainly Roemer postpones to correct, Shapiro postpones to correct, the total delay of solar system planet is corrected;
The computing method of the integer ambiguity in the described step (2):
ΔN = round ( n ^ · ( r E + r ~ ) / λ - Δφ ) - - - ( 5 )
Wherein, The satellite position predicted value that provides for the satellite dynamics equation;
GDOP value calculating method in the described step (2):
GDOP = trace ( C ) - - - ( 6 )
Wherein, trace represents to ask matrix trace, and C is the site error covariance matrix, calculate by following formula
C = { ( H T H ) - 1 H T } λ 1 2 σ φ 1 2 0 0 0 λ 2 2 σ φ 2 2 0 0 0 λ 3 2 σ φ 3 2 { ( H T H ) - 1 H T } T - - - ( 7 )
Wherein, H is the matrix of coefficients of measurement equation, λ i(i=1,2,3) are the impulse radiation wavelength of i pulsar, σ φ iIt is the phase measurement accuracy of i pulsar;
In the described step (3), apart from the measurement equation that is observed quantity be with starlight angular distance and the earth's core
Z 1 = arccos ( - r · s r ) + v 1 - - - ( 8 )
Z 2 = x 2 + y 2 + z 2 + v 2 - - - ( 9 )
Wherein, Z 1Be the observed quantity of starlight angular distance, Z 2For the earth's core apart from observed quantity, r is satellite body system the earth's core vector down, s is the unit vector of satellite body system time starlight direction, v 1For the starlight angular distance is measured noise, v 2For the earth's core apart from measuring noise, (x, y z) be the satellite position coordinate figure of representing under Earth central inertial system;
In the described step (4), radar altimeter vertically be installed on satellite just below, the instrumented satellite platform is to the distance on the face of land, and geoid surface adopts the rich ellipsoidal model of gram Lay;
In the described step (4), the actual sea level elevation h of sub-satellite point (x, y) landform at random that is generated by accompanying software provides, and landform generation method is as follows at random:
Landform is decomposed into landform reference plane height h on short transverse 0With the topographic relief Z that on this plane, superposes (x, y), promptly (x, actual sea level elevation h y) (x y) is expressed as:
h(x,y)=h 0+Z(x,y)(10)
Artificial landform is produced by two-dimentional single order discrete autoregressive process:
Z(x i,y i)=a 1Z(x i-Δx,y i)+a 2Z(x i,y i-Δy)+a 3Z(x i-Δx,y i-Δy)+W(x i,y i)(11)
Wherein, Z (x i, y i) be (x i, y i) point the topographic relief height, Δ x, Δ y are respectively along x, the sampling interval of y direction, W (x i, y i) be the zero-mean white noise sequence, and W (x i, y i)~N (0, σ w 2);
Isotropy and smooth performance according to landform make coefficient a 1=a 2=a, a 3=b, and as a=exp (1/T Auc), b=-a 2, σ w 2 = ( d RMS ) 2 exp ( - ( i + j ) / T auc ) The time, bidimensional random function Z (x in the formula (11) i, y i) related function be
R ij = ( d MRS ) exp ( - i + j T auc ) - - - ( 12 )
The bidimensional stochastic process that obtains thus, its mean square deviation are d RMS, by exponential damping, the decorrelation coefficient is T Auc, artificially generated terrain is produced respectively by two one-dimensional random processes in the border landform of x direction and y direction, on the landform border of x direction, and y i=y 0, then have
Z(x i,y 0)=aZ(x i-Δx,y 0)+W(x i,y 0)(13)
On the landform border of y direction, x i=x 0, then have
Z(x 0,y i)=aZ(x 0,y i-Δy)+W(x 0,y i)(14)
By setting h 0, d RMS, T AucThree parameters utilize software programming to produce different landform at random, and the Terrain Elevation at random that produces is replaced actual sea level elevation h, and (x y) is added on the geoid surface;
In the described step (5), apart from the discrete form measurement equation of observed quantity be with the earth's core direction and the earth's core
z = u k r k + v k - - - ( 15 )
Wherein, u kBe the unit vector of the earth's core direction, r kBe the earth's core distance, v kBe measurement noise;
In the described step (2) (3) (4) (5), the subfilter of four subsystems adopts EKF method or Unscented kalman filter method, and the time of independently carrying out respectively upgrades and measures and upgrade;
In the described step (6), information fusion is carried out in the output of (5) four subsystems of above-mentioned steps (2) (3) (4) handle, the method for employing is as follows:
The measurement equation of the individual subsystem of i (i=1,2,3,4) is:
Z i(k)=H i(k)X(k)+V i(k),i=1,2,3,4(16)
The initial time global state is estimated as , its covariance matrix is P G0,, this information is passed through the information distribution factor-beta according to the information conservation principle iBe assigned to four subfilters and senior filter, distribution principle is as follows:
P i 0 - 1 ( k ) = β i P g 0 - 1 ( k ) X ^ i 0 ( k ) = X ^ g 0 ( k ) Q i 0 - 1 ( k ) = β i Q g 0 - 1 ( k ) i = 1,2,3,4 , m - - - ( 17 )
Wherein, i=1, four subfilters of 2,3,4 expressions, m represents senior filter, the information distribution factor-beta iSatisfy the information conservation principle: Σ i = 1 4 , m β i = 1 , Get β herein m=0, β 1234=1/4, each subfilter time of independently carrying out respectively upgrades, and utilizes the measurement information of its respective sensor to measure renewal, gets the partial estimation value of four subfilters With evaluated error covariance matrix P i(k) i=1,2,3,4; In senior filter, the local message that subfilter is exported carries out information fusion by following formula, obtains the global state estimated information and is
X ^ g = P g Σ i = 1 4 P i - 1 X ^ i P g - 1 = P 1 - 1 + P 2 - 1 + P 3 - 1 + P 4 - 1 - - - ( 18 )
Promptly carry out information fusion according to following formula, the output global state is estimated With its covariance matrix P g, the global state of output is estimated As the initial value of satellite orbit kinetics equation predicted satellite states, and feed back to the dynamics of orbits model of satellite, as next initial value of satellitosis constantly of satellite orbit dynamic forecasting;
In the described step (7), χ 2The method that method of inspection carries out system failure detection and isolation is as follows:
The χ that in each subfilter, all adds fault detect and isolation 2Method, constitute fault-tolerant federal wave filter, utilize the fault detection module of four subsystems, respectively to based on the autonomous navigation of satellite subsystem of X-ray detector, based on the autonomous navigation of satellite subsystem of infrared horizon and star sensor, based on the autonomous navigation of satellite subsystem of radar altimeter, whether detect based on the autonomous navigation of satellite subsystem fault of three sensors of ultraviolet;
Each subfilter is carried out fault detect and isolation after measuring and upgrading, and its method is at first to design the fault detect function D of four subsystems i(k) i=1,2,3,4,
d i ( k ) = z i ( k ) - H i ( k ) X ^ i ( k / k - 1 ) - - - ( 19 )
S i ( k ) = H i ( k ) P i ( k / k - 1 ) H i T ( k ) + R i ( k ) - - - ( 20 )
D i ( k ) = d i T ( k ) S i - 1 ( k ) d i ( k ) - - - ( 21 )
Carry out fault judgement: D (k)>T then DThe time, fault is arranged, D (k)<T DThe time, non-fault; T DBe predefined thresholding, can determine by the alert rate of mistake, as the alert rate P of mistake FaDuring=α, by P Fa=P[λ k>T D| H 0]=α solves thresholding T D
The present invention's advantage compared with prior art is:
(1) the present invention obtains the higher navigation information of precision with the measurement information fusion of multiple sensors, has the advantage that the high-precision independent navigation information is provided for satellite; The present invention simultaneously adopts the multisensor redundant configuration, and adopts system-level fault detect and partition method, has improved the reliability and the fault-tolerant ability of system.
(2) the autonomous navigation of satellite subsystem among the present invention based on X-ray detector, having provides the autonomous orbit determination of ten meters magnitudes precision ability, and is applicable to the high-precision independent navigation demand of deep space probe; Based on the autonomous navigation of satellite subsystem of radar altimeter, adopt the vertically arranged mode of radar altimeter, use rational sea level model and artificially generated terrain, improved the navigation accuracy of satellite.
(3) each independent navigation subsystem among the present invention, all can independently constitute the autonomous navigation of satellite system, thereby can under different situations, select for use different subsystems to make up, make up new integrated navigation system, also can only use a kind of subsystem wherein, as the autonomous navigation system of satellite.Be adapted under the different navigation accuracy requirement, to the control of navigational system cost and weight, volume.
(4) the present invention is applicable to the fly high-precision independent navigation in each stage of lunar orbiter and Mars probes.In each stage of lunar orbiter and Mars probes flight, star sensor, three sensors of ultraviolet, X-ray detector can both provide enough navigation observed quantities for navigational system, thereby provide the high precision navigation information for spacecraft.In around-the-moon flight with around the Mars mission phase, the radar altimeter ranging information provides the distance between spacecraft and the moon or the martian surface, thereby determine around the moon or around the orbit altitude of fire flight spacecraft, for spacecraft provides the high precision navigation information, realize the navigation of the complete autonomous type of survey of deep space aircraft.
Description of drawings
Fig. 1 is the fault-tolerant federal filter graph architecture of integrated navigation system of the present invention;
Fig. 2 among the present invention based on the pulsar navigation synoptic diagram of the autonomous navigation of satellite subsystem of X-ray detector.
Embodiment
Autonomous navigation of satellite system based on many information fusion is made of four subsystems and a senior filter, and as shown in Figure 1, specific implementation method of the present invention is as follows:
One, foundation is based on the satellite motion state equation of dynamics of orbits
(1) foundation contains the satellite orbit kinetics equation of solar-lunar perturbating
Under J2000.0 Earth central inertial system, satellite orbit kinetics equation (system state equation) is
X · = f ( X , t ) + W ( t ) - - - ( 1 )
In the formula, X=[x, y, z, v x, v y, v z] TBe state variable, x, y, z, v x, v y, v zBe respectively satellite at X, Y, the position and the velocity amplitude of three directions of Z.F (X, form t) is according to the difference of the set perturbation of satellite, under the trisome gravitation situation of the non-spherical earth perturbation of only considering the J2 item and day, month,
f ( X , t ) = v x ; v y ; v z - μ e x r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 1.5 ) ] + μ s ( x s - x | r s - r | - x s | r s | ) + μ m ( x m - x | r m - r | - x m | r m | ) - μ e y r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 1.5 ) ] + μ s ( y s - y | r s - r | - y s | r s | ) + μ m ( y m - y | r m - r | - y m | r m | ) - μ e z r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 4.5 ) ] + μ s ( z s - z | r s - r | - z s | r s | ) + μ m ( z m - z | r m - r | - z m | r m | ) - - - ( 2 )
Wherein, μ e, μ s, μ mBe respectively the gravitational constant of the earth, the sun, the moon, R eBe earth radius, J 2Be terrestrial gravitation coefficient, r s=[x s, y s, z s] TAnd r m=[x m, y m, z m] TRepresent the position of the Sun and the Moon under Earth central inertial system respectively, W (t) is the model error vector, has represented the influence of the perturbative forces such as higher order term, solar radiation pressure perturbation atmospherical drag perturbation of perturbation of earths gravitational field, is assumed to be the zero-mean white Gaussian noise.
Two, foundation is based on the autonomous navigation of satellite subsystem of X-ray detector
(2) the navigation system of selection of X ray pulsar
It is as follows with the selection principle of X ray pulsar to navigate: millisecond pulsar is main candidate's pulsar; The pulsar rotation period that is used to navigate is stable, and the pulsar with rotation transition should foreclose; And the background radiation of the power of consideration pulsar radiation signal, pulse profile shape, pulsar place direction, and select basic parameter to measure the high pulsar of precision; All pulsars evenly distribute in the space as far as possible.The pulsar quality factor are calculated formula:
Q = 754.8 SNR T = 1 5 T 10 + 7 T 50 5 T 10 - T 50 T 50 ( T 10 - T 50 ) - - - ( 3 )
In the formula, SNR is the signal to noise ratio (S/N ratio) of pulsar, and T is the pulsar swing circle, T 50And T 10Be respectively the pulse width of pulsating wave peak intensity 50% and at 10% o'clock.The figure of merit value of calculated candidate pulsar, the Q value is big more, and the navigation attribute of pulsar is good more, more can be as the navigation pulsar.
(3) setting up with pulsar rotation phase place (TOA) is the measurement equation of observed quantity
Fig. 2 is the pulsar navigation synoptic diagram around the ground satellite.Wherein, Sat, E, SSB represent satellite, the earth and solar system barycenter respectively, Be the direction of visual lines of certain pulsar, T is a recurrence interval; r EAnd r SSB_satRepresent the earth, the satellite position vector with respect to SSB respectively, r represents the position vector of satellite with respect to the earth's core.The pulsion phase place value that epoch of observation, the satellite place recorded during t is φ Sat, the phase value of solar system barycenter is φ SSB
The pulsar clock model that is defined in solar system barycenter is
Φ ( t ) = φ ( t 0 ) + f ( t - t 0 ) + Σ m = 1 M f ( m ) ( t - t 0 ) m + 1 ( m + 1 ) ! - - - ( 4 )
In the formula, φ (t 0) be with reference to t epoch 0The time the pulsar phase place, f is that pulsar is at t 0The time rotation frequency, f (m)M order derivative (generally getting m=1,2,3) for f.The pulsion phase place value Φ at solar system barycenter place in the time of can obtaining t by the forecast of (4) formula SSB, get its fraction part φ SSB(0≤φ SSB<1).
Integer ambiguity between satellite and the SSB is Δ N, and the wavelength of the impulse radiation of pulsar is λ (λ=cT, c are the light velocity), then measures equation and is
[ Δφ + ΔN ] · λ - n ^ · r E = n ^ · r + c δt sat + Δ rel + v k - - - ( 5 )
Wherein, Δ φ=φ SSBSat, v kFor impulse phase is measured noise.Δ RelBe the relativistic effect correction member, comprise that mainly Roemer postpones to correct, Shapiro postpones to correct, the total delay of solar system planet is corrected.Earth vector r EProvide by the JPL ephemeris.Use the observed quantity of three pulsars, then measurement equation is
( Δφ 1 + ΔN 1 ) · λ 1 - n ^ 1 · r E ( Δφ 2 + ΔN 2 ) · λ 2 - n ^ 2 · r E ( Δφ 3 + ΔN 3 ) · λ 3 - n ^ 3 · r E = H · r cδt sat + Δ rel 1 Δ rel 2 Δ rel 3 + V k - - - ( 6 )
Wherein H = n ^ 1 1 n ^ 2 1 n ^ 3 1 .
(4) calculate the integer ambiguity value
The satellite position predicted value that uses the satellite orbit precursor to provide , integer ambiguity for Δ N is
ΔN = round ( n ^ · ( r E + r ~ ) / λ - Δφ ) - - - ( 7 )
(5) calculate the GDOP value
When many X ray pulsars are visible, select three minimum pulsars of geometric dilution of precision (GDOP), as the navigation calculation pulsar.For the situation of observing three pulsars simultaneously, its GDOP value calculating method is as follows:
The site error covariance matrix is
cov ( P ) = { ( H T H ) - 1 H T } λ 1 2 σ φ 1 2 0 0 0 λ 2 2 σ φ 2 2 0 0 0 λ 3 2 σ φ 3 2 { ( H T H ) - 1 H T } T - - - ( 8 )
σ wherein φ i(i=1,2,3) are the phase measurement accuracy of i pulsar.Then the GDOP value is
GDOP = trace ( C ) = σ x 2 + σ y 2 + σ z 2 + σ t 2 - - - ( 9 )
Wherein, trace represents to ask matrix trace, σ x 2, δ y 2, σ z 2, σ t 2Be the diagonal entry of Matrix C, represent x respectively, y, the measuring accuracy of the distance of z axle and time variable t.
(6) filtering method of subfilter
Subfilter adopts the non-linear Kalman filtering method, and EKF method (EKF) or the Unscented Kalman filtering time of independently carrying out respectively of can considering upgrades and measures and upgrade.
The process that EKF is handled is as follows:
State equation and measurement equation are carried out discretize, and the discrete type linear disturbance equation and the discrete type measurement equation that obtain state equation are
δX k+1=φ k+1,kδX k+W k(10)
Z k+1=H k+1X k+1+V k+1(11)
φ in the formula K+1, kBe state-transition matrix, Z K+1For on three pulsar direction of visual lines between satellite and the solar system barycenter apart from observed quantity.The covariance matrix of state model noise is E [ W k W k T ] = Q ( k ) , The covariance matrix of measurement noise is E [ V k + 1 V k + 1 T ] = R ( k ) , W kWith V K+1Uncorrelated mutually.
Utilize state equation (10) and measurement equation (11) composition subfilter 1 as shown in Figure 1, adopt the EKF method to carry out information fusion and handle, satellite position, velocity information that the output subfilter obtains.The EKF fundamental equation is as follows:
X ^ k + 1 , k = X ^ k , k + f [ X ^ k , k , t k ] · T X ^ k + 1 / k + 1 = X ^ k + 1 / k + K k + 1 ( Z k + 1 - H k + 1 X k + 1 / k ) P k + 1 / k = φ k + 1 / k P k / k φ k + 1 / k T + Q k K k + 1 = P k + 1 / k H k + 1 T [ H k + 1 P k + 1 / k H k + 1 T + R k + 1 ] - 1 P k + 1 / k + 1 = ( I - K k + 1 H k + 1 ) P k + 1 / k ( I - K k + 1 H k + 1 ) T + K k + 1 R k + 1 K k + 1 T - - - ( 12 )
In the formula, f represents the satellite orbit kinetics equation, and T is the Kalman filtering cycle, φ K+1, kBe state-transition matrix, K K+1Be kalman gain coefficient, H K+1Be measurement equation matrix of coefficients, P K+1, kBe optimum prediction valuation error covariance matrix, P K+1, k+1Be optimal filtering value error covariance matrix, Be the one-step prediction value, Be the Kalman filtering value.
The process that the Unscented Kalman filtering is handled is as follows:
For system equation and the following discrete system of measurement equation
X k + 1 = F ( X k , k ) + W k Z k + 1 = H ( X k ) + V k - - - ( 13 )
The dimension of the system state variable of setting up departments is n * 1 dimension, and 2n+1 sampled point is so
χ 0 = X ^
χ i = X ^ + ( ( n + λ ) P k / k ) i , i = 1,2 , . . . , n - - - ( 14 )
χ i = X ^ - ( ( n + λ ) P k / k ) i , i = n + 1 , . . . , 2 n
The weight of each sampled point correspondence is
W 0 ( m ) = λ / ( n + λ )
W 0 ( c ) = λ / ( n + λ ) + ( 1 - α 2 + β ) - - - ( 15 )
W i ( m ) = W i ( c ) = λ / { 2 ( n + λ ) } , i = 1,2 , . . . , 2 n
Wherein, λ=α 2(n+k)-and n is a scalar parameter, constant alpha determines that sampled point centers on Distribution characteristics, be set to little positive number (1 〉=α 〉=10 usually -4), constant k is a scalar parameter, is set to 2 or 3-n usually.W i (m)Be the used weights of average weighting, W i (c)Be the used weights of covariance-weighted.The basic facilities process of Unscented Kalman filtering is as follows:
The calculating sampling point
χ k = X ^ k X ^ k + n + λ P k / k X ^ k - n + λ P k / k - - - ( 16 )
Time upgrades
χ k + 1 / k * = F ( χ k , k ) , X ^ k + 1 / k = Σ 0 2 n W i ( m ) χ k + 1 / k * P k + 1 / k = Σ i = 0 2 n W i ( c ) [ χ i , k + 2 / k * - X ^ k + 1 / k ] [ χ i , k + 1 / k * - X ^ k + 1 / k ] T + Q k + 1 χ k + 1 / k = X ^ k + 1 / k X ^ k + 1 / k + γ P k + 1 / k X ^ k + 1 / k - γ P k + 1 / k z k + 1 / k = H ( χ k + 1 / k ) , Z ^ k + 1 / k = Σ 0 2 n W i ( m ) z k + 1 / k - - - ( 17 )
Measure and upgrade
P Z ^ k + 1 , Z ^ k + 1 = Σ i = 0 2 n W i ( c ) [ z i , k + 1 / k - Z ^ k + 1 / k ] [ z i , k + 1 / k - Z ^ k + 1 / k ] T + R k + 1 P X ^ k + 1 , Z ^ k + 1 = Σ i = 0 2 n W i ( c ) [ χ i , k + 1 / k - X ^ k + 1 / k ] [ z i , k + 1 / k - Z ^ k + 1 / k ] T + R k + 1 K k + 1 = P X ^ k + 1 , Z ^ k + 1 P Z ^ k + 1 , Z ^ k + 1 - 1 , X ^ k + 1 / k + 1 = X ^ k + 1 / k + K k + 1 ( Z k + 1 - Z ^ k + 1 / k ) P k + 1 / k + 1 = P k + 1 / k - K k + 1 P Z ^ k + 1 , Z k + 1 K k + 1 T - - - ( 18 )
Wherein, the covariance matrix of state model noise is E = [ W k W k T ] = Q k , The covariance matrix of measurement noise is E [ V k + 1 V k + 1 T ] = R k , W kWith V K+1Uncorrelated mutually.
Utilize state equation (1) and measurement equation (6) to form subfilter 1, equation (1) and (6) are carried out discretize respectively, adopt EKF method or Unscented kalman filter method to carry out information fusion and handle, satellite position, velocity information that the output subfilter obtains.
Three, foundation is based on the autonomous navigation of satellite subsystem of infrared horizon and star sensor
(7) setting up with starlight angular distance and the earth's core distance is the measurement equation of observed quantity
Star sensor observation navigation fixed star can be determined the direction of starlight in the satellite body coordinate system; Utilize infrared horizon to record the direction of the earth's core vector in the satellite body coordinate system; Obtain the angle between starlight vector and the earth's core vector like this, i.e. the starlight angular distance.Infrared horizon observation earth infrared image obtains earth half angle ρ, earth radius R eFor known, thereby can calculate the earth's core distance r = R e sin ρ , Thereby apart from the measurement equation that is observed quantity be with starlight angular distance and the earth's core
Z 1 = arccos ( - r · s r ) + v 1 - - - ( 19 )
Z 2 = x 2 + y 2 + z 2 + v 2 - - - ( 20 )
Wherein, Z 1Be the observed quantity of starlight angular distance, Z 2For the earth's core apart from observed quantity, r is satellite body system the earth's core vector down, s is the unit vector of satellite body system time starlight direction, v 1For the starlight angular distance is measured noise, v 2For the earth's core apart from measuring noise, (x, y z) be the satellite position coordinate figure of representing under Earth central inertial system.
(8) the measurement equation group of measurement equation (19), (20) formation, utilize state equation (2) and measurement equation group structure subfilter 2 as shown in Figure 1, repeat the EKF method or the Unscented Kalman filtering disposal route of above-mentioned steps (6) and carry out data processing, position, the velocity estimation value of output satellite.
Four, foundation is based on the autonomous navigation of satellite subsystem of radar altimeter
(9) setting up with starlight angular distance and orbit altitude is the measurement equation of observed quantity
Radar altimeter records the distance of satellite platform to the face of land vertically according under satellite, i.e. the orbit altitude measurement information of satellite.Geoid surface adopts the rich ellipsoidal model of gram Lay, and the fluctuating on sea level is taken in as correction term.The rich level surface of Ke Lai can be expressed as following formula
Rich level surface can be expressed as following formula
Wherein, Be geocentric latitude, P 2For secondary is reined in to such an extent that allow polynomial expression, J 2Be humorous coefficient of earth zone, R is an earth mean radius value, parameter alpha eFor
α e = [ 2 3 J 2 ( α e R + m 2 ) ] [ 1 - 2 3 m ] m ≈ 0.003449963 - - - ( 22 )
The measurement equation that then with the orbit altitude is observed quantity is
Z 2 = r - R [ 1 - 2 α e 3 ( 0.5 ( z r ) 2 - 1 2 ) ] - h ( x , y ) + v 2 - - - ( 23 )
In the formula, Z 2Be the satellite orbital altitude that radar altimeter records, r is the distance of satellite to the earth's core, v 2For measuring noise, (x y) is the Terrain Elevation of substar to h, by landform software generation at random.With the starlight angular distance is the same equation of measurement equation (19) of observed quantity.
(10) landform generation method at random
Set up artificially generated terrain, the variation of simulation actual landform.Landform is decomposed into landform reference plane height h on short transverse 0With the topographic relief Z that on this plane, superposes (x, y), promptly (x, actual sea level elevation h y) (x y) can be expressed as:
h(x,y)=h 0+Z(x,y)(24)
Artificial landform is produced by two-dimentional single order discrete autoregressive process:
Z(x i,y i)=a 1Z(x i-Δx,y i)+a 2Z(x i,y i-Δy)+a 3Z(x i-Δx,y i-Δy)+W(x i,y i)(25)
Wherein, Z (x i, y i) be (x i, y i) point the topographic relief height, Δ x, Δ y are respectively along x, the sampling interval of y direction, W (x i, y i) be the zero-mean white noise sequence, W (x i, y i)~N (0, σ w 2).
Isotropy and smooth performance according to landform can make coefficient a 1=a 2=a, a 3=b, and as a=exp (1/T Auc), b=-a 2, σ w 2 = ( d RMS ) 2 exp ( - ( i + j ) / T auc ) The time, bidimensional random function Z (x in the formula (25) i, y i) related function be
R ij = ( d RMS ) exp ( - i + j T auc ) - - - ( 26 )
The bidimensional stochastic process that obtains thus, its mean square deviation are d RMS, by exponential damping, the decorrelation coefficient is T AucArtificially generated terrain is produced respectively by two one-dimensional random processes in the border landform of x direction and y direction.If on the landform border of x direction, y i=y 0, then have
Z(x i,y 0)=aZ(x i-Δx,y 0)+W(x i,y 0)(27)
If on the landform border of y direction, x i=x 0, then have
Z(x 0,y i)=aZ(x 0,y i-Δy)+W(x 0,y i)(28)
By setting h 0, d RMS, T AucThree parameters utilize software programming can produce different landform at random.(x y) is added on the geoid surface with the h of landform at random that produces.
(11) the measurement equation group of measurement equation (19), (23) formation, utilize state equation (2) and measurement equation group structure subfilter 3 as shown in Figure 1, repeat the EKF method or the Unscented Kalman filtering disposal route of above-mentioned steps (6) and carry out data fusion, position, the velocity estimation value of output satellite.
Five, foundation is based on the autonomous navigation of satellite subsystem of three sensors of ultraviolet
(12) foundation is based on the autonomous navigation of satellite subsystem of three sensors of ultraviolet
The responsive earth ultraviolet image of three sensors of ultraviolet extracts the earth's core direction in the celestial body system and the earth's core apart from information by Flame Image Process; Utilize the attitude matrix of satellite to be converted to the position vector of satellite in Earth central inertial system, thereby obtain the earth's core direction and the earth's core apart from observed quantity, the measurement equation of its discrete form is
z = u k r k + v k - - - ( 29 )
Wherein, u kBe the unit vector of the earth's core direction, r kBe the earth's core distance, v kBe measurement noise.
Utilize state equation (2) and measurement equation (29) structure subfilter 4 as shown in Figure 1, repeat the EKF method or the Unscented Kalman filtering disposal route of above-mentioned steps (6) and carry out data fusion, position, the velocity estimation value of output satellite.
Six, the information fusion of senior filter, fault detect and partition method
(13) information fusion of senior filter is handled
The measurement equation of the individual subsystem of i (i=1,2,3,4) is:
Z i(k)=H i(k)X(k)+V i(k),(i=1,2,3,4)(30)
The initial time global state is estimated as , its covariance matrix is P G0According to the information conservation principle, this information is passed through the information distribution factor-beta iBe assigned to four subfilters and senior filter, distribution principle is as follows:
P i 0 - 1 ( k ) = β i P g 0 - 1 ( k ) X ^ i 0 ( k ) = X ^ g 0 ( k ) Q i 0 - 1 ( k ) = β i Q g 0 - 1 ( k ) ( i = 1,2,3,4 , m ) - - - ( 31 )
Wherein, i=1, four subfilters of 2,3,4 expressions, m represents senior filter, the information distribution factor-beta iSatisfy the information conservation principle: Σ i = 1 4 m β i = 1 , Get β herein m=0, β 1234=1/4.Each subfilter time of independently carrying out respectively upgrades, and utilizes the measurement information of its respective sensor to measure renewal, gets the partial estimation value of four subfilters With evaluated error covariance matrix P i(k) (i=1,2,3,4), concrete measure such as equation (12) or (17) and (18).In senior filter, the local message that subfilter is exported carries out information fusion by following formula, obtains the global state estimated information and is
X ^ g = P g Σ i = 1 4 P i - 1 X ^ i P g - 1 = P 1 - 1 + P 2 - 1 + P 3 - 1 + P 4 - 1 - - - ( 32 )
Promptly carry out information fusion according to (32) formula, the output global state is estimated With its covariance matrix P gThe global state of output is estimated As the initial value of satellite orbit kinetics equation predicted satellite states, and feed back to the dynamics of orbits model of satellite, as next initial value of satellitosis constantly of satellite orbit dynamic forecasting.
In the federal filtering, senior filter and four subfilter associated working, information distribution is only carried out at initial time, and subfilter works alone, and does not receive the feedback information of senior filter.
(14) fault detect of integrated navigation and partition method
Adopt system-level fault detect and partition method, get rid of the sensor and the corresponding subsystem that break down.The χ that in each subfilter, all adds fault detect and isolation 2Method constitutes fault-tolerant federal wave filter, and its structure as shown in Figure 1.FDI 1, FDI 2, FDI 3, FDI 4 are respectively the fault detection module of four subsystems among Fig. 1, respectively whether subsystem 1,2,3,4 faults are detected, and whether decision is isolated, so that senior filter is recombinated to normal subsystem, after making part navigation sensor fault cause subsystem failure, integrated navigation system still can work on, the tool fault-tolerance.
Each subfilter time of independently carrying out respectively upgrades and measures and upgrades.Each subfilter is carried out fault detect and isolation after measuring and upgrading.Its method is at first to design the fault detect function D of four subsystems i(k) (i=1,2,3,4):
d i ( k ) = z i ( k ) - H i ( k ) X ^ i ( k / k - 1 ) - - - ( 33 )
S i ( k ) = H i ( k ) P i ( k / k - 1 ) H i T ( k ) + R i ( k ) - - - ( 34 )
D i ( k ) = d i T ( k ) S i - 1 ( k ) d i ( k ) - - - ( 35 )
Carry out fault judgement: D (k)>T then DThe time, fault is arranged; D (k)<T DThe time, non-fault.T DBe predefined thresholding, can determine by the alert rate of mistake.As the alert rate P of mistake FaDuring=α, P Fa=P[λ k>T D| H 0]=α solves thresholding T DIf judge the subsystem non-fault, then its filtering result is delivered to senior filter; If subsystem fault, then filtering result does not deliver to senior filter, and to its isolation, recovers normal until this subsystem.
Calculate separately partial estimation in each subfilter, and after carrying out fault detect, the common condition of the normal subfilter of output is carried out information fusion, obtain overall estimated information by senior filter.
(15) outgoing position, velocity information.
Carry out Computer Simulation according to above-mentioned steps (1)-(15), set up dynamics of orbits equation and measurement equation, utilize EKF method or Unscented kalman filter method that Filtering Processing is carried out in the observed quantity of four subsystems, the global state that adopts the information fusion means to obtain quantity of state is again estimated X ^ g = x ^ y ^ z ^ v ^ x v ^ y v ^ z T With its covariance matrix P g = p x p y p z p v , p v , p v T , Wherein Be respectively to position, the speed amount x of satellite in X, Y, three directions of Z, y, z, v under Earth central inertial system x, v y, v zOptimal estimation.
The content that is not described in detail in the instructions of the present invention belongs to this area professional and technical personnel's known prior art.

Claims (1)

1, a kind of spacecraft shading device combined navigation methods based on many information fusion is characterized in that: specifically may further comprise the steps:
(1) foundation is based on the satellite motion state equation of dynamics of orbits;
(2) foundation is based on the autonomous navigation of satellite subsystem of X-ray detector;
(3) foundation is based on the autonomous navigation of satellite subsystem of infrared horizon and star sensor;
(4) foundation is based on the autonomous navigation of satellite subsystem of radar altimeter;
(5) foundation is based on the autonomous navigation of satellite subsystem of three sensors of ultraviolet;
(6) information fusion is carried out in the output of (5) four subsystems of above-mentioned steps (2) (3) (4) and handle, the output navigation information;
(7) adopt χ 2Method of inspection detects system and the isolated fault subsystem that breaks down;
In the described step (1), foundation comprises J 2Item non-spherical earth perturbation and the satellite orbit kinetics equation of solar-lunar perturbating, under J2000.0 Earth central inertial system, the satellite orbit kinetics equation is
X · = f ( X , t ) + W ( t ) - - - ( 1 )
In the formula, X=[x, y, z, v x, v y, v z] TBe state variable, x, y, z, v x, v y, v zBe respectively satellite at X, Y, the position and the velocity amplitude of three directions of Z; F (X, form t) is according to the difference of the set perturbation of satellite, only considers under the trisome gravitation situation of the non-spherical earth perturbation of J2 item and day, month,
f ( X , t ) = v x ; v y ; v z - μ e x r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 1.5 ) ] + μ s ( x s - x | r s - r | - x s | r s | ) + μ m ( x m - x | r m - r | - x m | r m | ) - μ e y r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 1.5 ) ] + μ s ( y s - y | r s - r | - y s | r s | ) + μ m ( y m - y | r m - r | - y m | r m | ) - μ e z r 3 [ 1 - J 2 ( R e r ) 2 ( 7.5 z 2 r 2 - 4.5 ) ] + μ s ( z s - z | r s - r | - z s | r s | ) + μ m ( z m - z | r m - r | - z m | r m | ) - - - ( 2 )
Wherein, μ e, μ s, μ mBe respectively the gravitational constant of the earth, the sun, the moon, R eBe earth radius, J 2Be terrestrial gravitation coefficient, r s=[x s, y s, z s] T and r m=[x m, y m, z m] TRepresent the position of the Sun and the Moon under Earth central inertial system respectively, x s, y s, z sBe the coordinate figure of the sun under Earth central inertial system, x m, y m, z mBe the coordinate figure of the moon under Earth central inertial system, x, y, z, v x, v y, v zBe respectively satellite at X, Y, the position and the velocity amplitude of three directions of Z, r is the distance that satellite is arrived in the earth's core, W (t) is the model error vector, represented the influence of perturbative force of higher order term, solar radiation pressure perturbation, the atmospherical drag perturbation of the perturbation of earth aspherical, is made as the zero-mean white Gaussian noise;
Described step (2) comprises the computing method of the quality factor of X ray pulsar, the computing method that concern between satellite position and the X ray pulsar rotation phase place, the computing method of integer ambiguity and GDOP value calculating method;
The computing method of the quality factor of the X ray pulsar in the described step (2) are
Q = 754.8 SNR T 1 5 T 10 + 7 T 50 5 T 10 - T 50 T 50 ( T 10 - T 50 ) - - - ( 3 )
In the formula, SNR is the signal to noise ratio (S/N ratio) of pulsar, and T is the pulsar swing circle, T 50And T 10Be respectively the pulse width of pulsating wave peak intensity 50% and at 10% o'clock;
The computing method that concern between satellite position in the described step (2) and the X ray pulsar rotation phase place are
[ Δφ + ΔN ] · λ - n ^ · r E = n ^ · r + c δt sat + Δ rel + v k - - - ( 4 )
Wherein, Δ φ=φ SSBSat, Δ N is the integer ambiguity value between satellite and the solar system barycenter SSB, λ is the wavelength of the impulse radiation of pulsar, Be the direction of visual lines of pulsar, r EFor the earth relatively and the position vector of SSB, r is the position vector of satellite with respect to the earth's core, c is the light velocity, δ t SatBe the clock correction on the satellite, v kFor impulse phase is measured noise, Δ RelBe the relativistic effect correction member, comprise that mainly Roemer postpones to correct, Shapiro postpones to correct, the total delay of solar system planet is corrected;
The computing method of the integer ambiguity in the described step (2):
ΔN = round ( n ^ · ( r E + r ~ ) / λ - Δφ ) - - - ( 5 )
Wherein, The satellite position predicted value that provides for the satellite dynamics equation;
GDOP value calculating method in the described step (2):
GDOP = trace ( C ) - - - ( 6 )
Wherein, trace represents to ask matrix trace, and C is the site error covariance matrix, calculate by following formula
C = { ( H T H ) - 1 H T } λ 1 2 σ φ 1 2 0 0 0 λ 2 2 σ φ 2 2 0 0 0 λ 3 2 σ φ 3 2 { ( H T H ) - 1 H T } T - - - ( 7 )
Wherein, H is the matrix of coefficients of measurement equation, λ i, i=1,2,3 is the impulse radiation wavelength of i pulsar, σ φ iIt is the phase measurement accuracy of i pulsar;
In the described step (3), apart from the measurement equation that is observed quantity be with starlight angular distance and the earth's core
Z 1 = arccos ( - r · s r ) + v 1 - - - ( 8 )
Z 2 = x 2 + y 2 + z 2 + v 2 - - - ( 9 )
Wherein, Z 1Be the observed quantity of starlight angular distance, Z 2For the earth's core apart from observed quantity, r is satellite body system the earth's core vector down, s is the unit vector of satellite body system time starlight direction, v 1For the starlight angular distance is measured noise, v 2For the earth's core apart from measuring noise, (x, y z) be the satellite position coordinate figure of representing under Earth central inertial system;
In the described step (4), radar altimeter vertically be installed on satellite just below, the instrumented satellite platform is to the distance on the face of land, and geoid surface adopts the rich ellipsoidal model of gram Lay;
In the described step (4), the actual sea level elevation h of sub-satellite point (x, y) landform at random that is generated by accompanying software provides, and landform generation method is as follows at random:
Landform can be decomposed into landform reference plane height h on short transverse 0With the topographic relief Z that on this plane, superposes (x, y), promptly (x, actual sea level elevation h y) (x y) can be expressed as:
h(x,y)=h 0+Z(x,y) (10)
Artificial landform is produced by two-dimentional single order discrete autoregressive process:
Z(x i,y i)=a 1Z(x i-Δx,y i)+a 2Z(x i,y i-Δy)+a 3Z(x i-Δx,y i-Δy)+W(x i,y i) (11)
Wherein, Z (x i, y i) be (x i, y i) point the topographic relief height, Δ x, Δ y are respectively along x, the sampling interval of y direction, W (x i, y i) be the zero-mean white noise sequence, and W (w i, y i)~N (0, σ w 2);
Isotropy and smooth performance according to landform make coefficient a 1=a 2=a, a 3=b, and as a=exp (1/T Auc), b=-a 2, σ w 2 = ( d RMS ) 2 exp ( - ( i + j ) / T auc ) The time, bidimensional random function Z (x in the formula (11) i, y i) related function be
R ij = ( d RMS ) exp ( - i + j T auc ) - - - ( 12 )
The bidimensional stochastic process that obtains thus, its mean square deviation are d RMS, by exponential damping, the decorrelation coefficient is T Auc, artificially generated terrain is produced respectively by two one-dimensional random processes in the border landform of x direction and y direction, on the landform border of x direction, and y i=y 0, then have
Z(x i,y 0)=aZ(x i-Δx,y 0)+W(x i,y 0) (13)
On the landform border of y direction, x i=x 0, then have
Z(x 0,y i)=aZ(x 0,y i-Δy)+W(x 0,y i) (14)
By setting h 0, d RMS, T AucThree parameters utilize software programming to produce different landform at random, and the landform at random that produces is added on the geoid surface;
In the described step (5), apart from the discrete form measurement equation of observed quantity be with the earth's core direction and the earth's core
z = u k r k + v k - - - ( 15 )
Wherein, u kBe the unit vector of the earth's core direction, r kBe the earth's core distance, v kBe measurement noise;
In the described step (2) (3) (4) (5), the subfilter of four subsystems adopts EKF method or Unscented kalman filter method, and the time of independently carrying out respectively upgrades and measures and upgrade;
In the described step (6), information fusion is carried out in the output of (5) four subsystems of above-mentioned steps (2) (3) (4) handle, the method for employing is as follows:
I, i=1, the measurement equation of 2,3,4 subsystems is:
Z i(k)=H i(k)X(k)+V i(k),i=1,2,3,4 (16)
The initial time global state is estimated as , its covariance matrix is P G0,, this information is passed through the information distribution factor-beta according to the information conservation principle iBe assigned to four subfilters and senior filter, distribution principle is as follows:
P i 0 - 1 ( k ) = β i P g 0 - 1 ( k ) X ^ i 0 ( k ) = X ^ g 0 ( k ) Q 10 - 1 ( k ) = β i Q g 0 - 1 ( k ) , i = 1,2,3,4 , m - - - ( 17 )
Wherein, i=1, four subfilters of 2,3,4 expressions, m represents senior filter, the information distribution factor-beta iSatisfy the information conservation principle: Σ i = 1 4 , m β i = 1 , Get β herein m=0, β 1234=1/4, each subfilter time of independently carrying out respectively upgrades, and utilizes the measurement information of its respective sensor to measure renewal, gets the partial estimation value of four subfilters With evaluated error covariance matrix P i(k) i=1,2,3,4; In senior filter, the local message that subfilter is exported carries out information fusion by following formula, obtains the global state estimated information and is
X ^ g = P g Σ i = 1 4 P i - 1 X ^ i P g - 1 = P 1 - 1 + P 2 - 1 + P 3 - 1 + P 4 - 1 - - - ( 18 )
Promptly carry out information fusion according to following formula, the output global state is estimated With its covariance matrix P g, the global state of output is estimated As the initial value of satellite orbit kinetics equation predicted satellite states, and feed back to the dynamics of orbits model of satellite, as next initial value of satellitosis constantly of satellite orbit dynamic forecasting;
In the described step (7), χ 2The method that method of inspection carries out system failure detection and isolation is as follows:
The χ that in each subfilter, all adds fault detect and isolation 2Method, constitute fault-tolerant federal wave filter, utilize the fault detection module of four subsystems, respectively to based on the autonomous navigation of satellite subsystem of X-ray detector, based on the autonomous navigation of satellite subsystem of infrared horizon and star sensor, based on the autonomous navigation of satellite subsystem of radar altimeter, whether detect based on the autonomous navigation of satellite subsystem fault of three sensors of ultraviolet;
Each subfilter is carried out fault detect and isolation after measuring and upgrading, and its method is at first to design the fault detect function D of four subsystems i(k) i=1,2,3,4,
d i ( k ) = z i ( k ) - H i ( k ) X ^ i ( k / k - 1 ) - - - ( 19 )
S i ( k ) = H i ( k ) P i ( k / k - 1 ) H i T ( k ) + R i ( k ) - - - ( 20 )
D i ( k ) = d i T ( k ) S i - 1 ( k ) d i ( k ) - - - ( 21 )
Carry out fault judgement: D (k)>T then DThe time, fault is arranged, D (k)<T DThe time, non-fault; T DBe predefined thresholding, determine, as the alert rate P of mistake by the alert rate of mistake FaDuring=α, by P Fa=P[λ k>T D| H 0]=α solves thresholding T D
CN200710191527A 2007-12-12 2007-12-12 Spacecraft shading device combined navigation methods based on many information fusion CN100575877C (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN200710191527A CN100575877C (en) 2007-12-12 2007-12-12 Spacecraft shading device combined navigation methods based on many information fusion

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN200710191527A CN100575877C (en) 2007-12-12 2007-12-12 Spacecraft shading device combined navigation methods based on many information fusion

Publications (2)

Publication Number Publication Date
CN101178312A CN101178312A (en) 2008-05-14
CN100575877C true CN100575877C (en) 2009-12-30

Family

ID=39404655

Family Applications (1)

Application Number Title Priority Date Filing Date
CN200710191527A CN100575877C (en) 2007-12-12 2007-12-12 Spacecraft shading device combined navigation methods based on many information fusion

Country Status (1)

Country Link
CN (1) CN100575877C (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101915581A (en) * 2010-07-21 2010-12-15 中国人民解放军信息工程大学 Comet optical surface signal simulation method for deep space exploration
CN102353378A (en) * 2011-09-09 2012-02-15 南京航空航天大学 Adaptive federal filtering method of vector-form information distribution coefficients
CN102520461A (en) * 2011-12-08 2012-06-27 中国空间技术研究院 Method for determining interference from NGSO satellite earth detection system to deep space detection system
TWI459016B (en) * 2011-09-30 2014-11-01 Maishi Electronic Shanghai Ltd Device, method and receiver for determining mobile information

Families Citing this family (50)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101281248B (en) * 2008-05-20 2010-05-12 北京航空航天大学 Multi-fault recognizing method applied to combined satellite navigation system
CN101598557B (en) * 2009-07-15 2012-05-23 北京航空航天大学 Integrated navigation system applied to pilotless aircraft
CN101696885B (en) * 2009-11-05 2012-05-23 中国人民解放军国防科学技术大学 Method for improving data processing precision of star sensors
CN101706281B (en) * 2009-11-13 2011-08-31 南京航空航天大学 Inertia/astronomy/satellite high-precision integrated navigation system and navigation method thereof
CN101788295A (en) * 2010-02-26 2010-07-28 南京信息工程大学 Combined navigation system of small-scale underwater vehicle and method thereof
CN101793526B (en) * 2010-04-12 2011-10-26 哈尔滨工业大学 Autonomous relative navigation method for multi-information fusion formation spacecrafts
CN101915577B (en) * 2010-07-21 2011-11-23 中国人民解放军信息工程大学 Comet optical point signal simulation method for deep space exploration
CN101958010B (en) * 2010-09-03 2011-12-28 清华大学 Correlation function test method for arranged effect of aircraft movement measuring sensor
CN102175260B (en) * 2010-12-31 2012-11-14 北京控制工程研究所 Error correction method of autonomous navigation system
CN102221363B (en) * 2011-04-12 2012-12-19 东南大学 Fault-tolerant combined method of strapdown inertial integrated navigation system for underwater vehicles
CN102243311B (en) * 2011-04-15 2013-03-20 中国人民解放军国防科学技术大学 Pulsar selection method used for X-ray pulsar navigation
US9562778B2 (en) * 2011-06-03 2017-02-07 Robert Bosch Gmbh Combined radar and GPS localization system
CN102506876B (en) * 2011-12-08 2014-07-02 北京控制工程研究所 Self-contained navigation method for measurement of earth ultraviolet sensor
CN102506877B (en) * 2011-12-08 2014-01-15 北京控制工程研究所 Deep-space exploration navigation system filtering method with immunity to initial error
CN102564454B (en) * 2011-12-23 2015-02-11 北京控制工程研究所 Measuring data simulation method for autonomous navigation system based on sun-earth-moon azimuth information
CN102654407B (en) * 2012-04-17 2014-07-02 南京航空航天大学 Multiple-fault detecting device and detecting method for tightly-integrated inertial satellite navigation system
CN102679985B (en) * 2012-05-11 2016-11-02 北京航空航天大学 A kind of apply between star follow the tracks of the decentralized autonomous navigation method of spacecraft constellation
CN103017736A (en) * 2012-11-30 2013-04-03 北京控制工程研究所 Method for determining exposure moment of star sensor data of satellite
CN102997922B (en) * 2012-11-30 2015-10-21 北京控制工程研究所 A kind of pulse arrival time difference defining method utilizing optical navigation information
CN102944416B (en) * 2012-12-06 2015-04-01 南京匹瑞电气科技有限公司 Multi-sensor signal fusion technology-based fault diagnosis method for wind turbine blades
CN103033188B (en) * 2012-12-24 2016-01-20 中国科学院国家授时中心 The autonomous method for synchronizing time of Navsat based on synthetic aperture observation
CN103218482B (en) * 2013-03-29 2017-07-07 南京航空航天大学 The method of estimation of uncertain parameter in a kind of dynamic system
CN103234538B (en) * 2013-04-07 2015-08-05 北京理工大学 The final Approach phase autonomous navigation method of a kind of planet
CN103323031B (en) * 2013-06-28 2015-11-18 上海新跃仪表厂 A kind of method of the horizon instrument systematic error online compensation quick based on star
CN103674022B (en) * 2013-12-19 2016-08-17 中国空间技术研究院 A kind of fast-pulse star navigation Carrier Phase Ambiguity Resolution method
CN103697915B (en) * 2013-12-24 2016-05-04 北京控制工程研究所 A kind of satellite sensor failure diagnosticability evaluation method of considering interference effect
CN103697894B (en) * 2013-12-27 2016-05-25 南京航空航天大学 Multi-source information unequal interval federated filter method based on the correction of wave filter variance battle array
CN103868514B (en) * 2014-03-20 2016-08-17 北京航天自动控制研究所 A kind of at orbit aerocraft autonomous navigation system
CN103884340B (en) * 2014-03-31 2016-08-17 北京控制工程研究所 A kind of information fusion air navigation aid of survey of deep space fixed point soft landing process
CN104006813A (en) * 2014-04-03 2014-08-27 中国人民解放军国防科学技术大学 Pulsar/starlight angle combination navigation method of high orbit satellite
CN103940424B (en) * 2014-04-14 2016-09-14 中国人民解放军国防科学技术大学 A kind of X-ray pulsar navigation signal integer ambiguity detection and method of estimation
CN103900577B (en) * 2014-04-14 2016-08-17 武汉科技大学 A kind of Relative Navigation towards formation flight tests the speed and Combinated navigation method
CN103913173B (en) * 2014-04-18 2017-03-01 中国人民解放军国防科学技术大学 Single X-ray pulsar navigation sees star sequence selection method
CN104063537B (en) * 2014-05-30 2017-04-19 北京控制工程研究所 Multi-body dynamics parameter determination system based on distributive time trigger and method thereof
US9656593B2 (en) * 2014-06-26 2017-05-23 The Boeing Company Flight vehicle autopilot
CN104296754B (en) * 2014-10-10 2017-09-19 北京大学 Autonomous navigation system and its autonomous navigation method based on laser space communication terminal
CN104864875B (en) * 2015-04-03 2018-01-05 北京控制工程研究所 A kind of spacecraft autonomic positioning method based on non-linear H ∞ filtering
EP3158411B1 (en) * 2015-05-23 2020-07-29 SZ DJI Technology Co., Ltd. Sensor fusion using inertial and image sensors
CN104897155B (en) * 2015-06-05 2018-10-26 北京信息科技大学 A kind of individual's portable multi-source location information auxiliary revision method
CN105259787B (en) * 2015-11-03 2018-06-08 中国电子科技集团公司第五十四研究所 A kind of Integrated Navigation Semi-physical Simulation tests synchronisation control means
CN106197408A (en) * 2016-06-23 2016-12-07 南京航空航天大学 A kind of multi-source navigation data fusion method based on factor graph
CN106441589B (en) * 2016-09-07 2018-12-25 北京航空航天大学 A kind of planet infra-red radiation emulation mode based on sliding-model control
CN106679693A (en) * 2016-12-14 2017-05-17 南京航空航天大学 Fault detection-based vector information distribution adaptive federated filtering method
CN107421533B (en) * 2017-06-22 2019-07-30 北京航空航天大学 A kind of deep space probe X-ray pulsar TOA/DTOA Combinated navigation method
CN108313329A (en) * 2018-04-03 2018-07-24 上海微小卫星工程中心 A kind of satellite platform data dynamic fusion system and method
CN108548542A (en) * 2018-07-13 2018-09-18 北京航空航天大学 A kind of LEO based on atmospheric drag acceleration analysis determines method
CN109373999B (en) * 2018-10-23 2020-11-03 哈尔滨工程大学 Integrated navigation method based on fault-tolerant Kalman filtering
CN110648524B (en) * 2019-08-27 2020-08-07 上海航天控制技术研究所 Multi-probe star sensor data transmission fault monitoring and autonomous recovery method
CN110780583A (en) * 2019-10-29 2020-02-11 中国科学院国家天文台 Moon-based pulsar time reference generation system
CN110933597B (en) * 2019-12-06 2021-02-19 北京壹氢科技有限公司 Bluetooth-based multi-unmanned vehicle collaborative fault-tolerant navigation positioning method and system

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5841398A (en) * 1996-11-20 1998-11-24 Space Systems/Loral, Inc. Integrated navigation and communication satellite system
CN1851406A (en) * 2006-05-26 2006-10-25 南京航空航天大学 Gasture estimation and interfusion method based on strapdown inertial nevigation system
CN101033973A (en) * 2007-04-10 2007-09-12 南京航空航天大学 Attitude determination method of mini-aircraft inertial integrated navigation system
CN101059349A (en) * 2007-05-18 2007-10-24 南京航空航天大学 Minitype combined navigation system and self-adaptive filtering method
CN101082494A (en) * 2007-06-19 2007-12-05 北京航空航天大学 Self boundary marking method based on forecast filtering and UPF spacecraft shading device

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5841398A (en) * 1996-11-20 1998-11-24 Space Systems/Loral, Inc. Integrated navigation and communication satellite system
CN1851406A (en) * 2006-05-26 2006-10-25 南京航空航天大学 Gasture estimation and interfusion method based on strapdown inertial nevigation system
CN101033973A (en) * 2007-04-10 2007-09-12 南京航空航天大学 Attitude determination method of mini-aircraft inertial integrated navigation system
CN101059349A (en) * 2007-05-18 2007-10-24 南京航空航天大学 Minitype combined navigation system and self-adaptive filtering method
CN101082494A (en) * 2007-06-19 2007-12-05 北京航空航天大学 Self boundary marking method based on forecast filtering and UPF spacecraft shading device

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
基于北斗双星定位辅助的SAR/INS组合导航系统研究. 熊智,冷雪飞,刘建业.宇航学报,第28卷第1期. 2007 *
基于联邦滤波的惯性导航姿态组合算法. 赖际舟,刘建业,刘瑞华,华冰.天津大学学报,第39卷第3期. 2006 *
多信息融合组合导航半物理仿真系统设计与实现. 邢广华,刘建业,孙永荣,熊智.航天控制,第23卷第2期. 2005 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101915581A (en) * 2010-07-21 2010-12-15 中国人民解放军信息工程大学 Comet optical surface signal simulation method for deep space exploration
CN101915581B (en) * 2010-07-21 2012-01-25 中国人民解放军信息工程大学 Comet optical surface signal simulation method for deep space exploration
CN102353378A (en) * 2011-09-09 2012-02-15 南京航空航天大学 Adaptive federal filtering method of vector-form information distribution coefficients
CN102353378B (en) * 2011-09-09 2013-08-21 南京航空航天大学 Adaptive federal filtering method of vector-form information distribution coefficients
TWI459016B (en) * 2011-09-30 2014-11-01 Maishi Electronic Shanghai Ltd Device, method and receiver for determining mobile information
CN102520461A (en) * 2011-12-08 2012-06-27 中国空间技术研究院 Method for determining interference from NGSO satellite earth detection system to deep space detection system
CN102520461B (en) * 2011-12-08 2014-07-02 中国空间技术研究院 Method for determining interference from NGSO satellite earth detection system to deep space detection system

Also Published As

Publication number Publication date
CN101178312A (en) 2008-05-14

Similar Documents

Publication Publication Date Title
Grewal et al. Global navigation satellite systems, inertial navigation, and integration
Tapley et al. Precision orbit determination for TOPEX/POSEIDON
Woffinden et al. Relative angles-only navigation and pose estimation for autonomous orbital rendezvous
Koop Global gravity field modelling using satellite gravity gradiometry
Giorgini et al. Predicting the Earth encounters of (99942) Apophis
Hoots et al. History of analytical orbit modeling in the US space surveillance system
Milani et al. Gravity field and rotation state of Mercury from the BepiColombo Radio Science Experiments
Hasan et al. A review of navigation systems (integration and algorithms)
Wolff Direct measurements of the Earth's gravitational potential using a satellite pair
Ashby Relativity in the global positioning system
CN103017774B (en) Pulsar navigation method with single detector
US7197381B2 (en) Navigational system and method utilizing sources of pulsed celestial radiation
Doornbos et al. Neutral density and crosswind determination from arbitrarily oriented multiaxis accelerometers on satellites
Reigber et al. CHAMP Phase B: Executive Summary
Vetter Fifty years of orbit determination
Blanchet et al. Relativistic theory for time and frequency transfer to order
How et al. Orion-A low-cost demonstration of formation flying in space using GPS
Campbell et al. Gravity field of the Jovian system from Pioneer and Voyager tracking data
Christensen et al. Calibration of TOPEX/POSEIDON at platform harvest
Shapiro et al. Measurement of the de Sitter precession of the Moon: A relativistic three-body effect
Müller et al. High performance clocks and gravity field determination
Psiaki et al. Ground tests of magnetometer-based autonomous navigation (MAGNAV) for low-earth-orbiting spacecraft
Kim Simulation study of a low-low satellite-to-satellite tracking mission
Sheikh et al. Recursive estimation of spacecraft position and velocity using X-ray pulsar time of arrival measurements
Fong et al. Constellation deployment for the FORMOSAT-3/COSMIC mission

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
C17 Cessation of patent right
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20091230

Termination date: 20131212