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

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 ²Method 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 χ ²Method of inspection detects system and the isolated fault subsystem that breaks down.

In the described step (1), foundation has comprised J ₂Item 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

$\stackrel{\·}{X}=f(X,t)+W\left(t\right)---\left(1\right)$

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)=\left[\begin{array}{c}{v}_x;v_y;v_z\\ -\frac{{\mathrm{\μ}}_{e}x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\μ}}_s(\frac{x_s-x}{|r_s-r|}-\frac{x_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{x_m-x}{|r_m-r|}-\frac{x_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\μ}}_{e}y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\μ}}_s(\frac{y_s-y}{|r_s-r|}-\frac{y_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{y_m-y}{|r_m-r|}-\frac{y_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\μ}}_{e}z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mathrm{\μ}}_s(\frac{z_s-z}{|r_s-r|}-\frac{z_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{z_m-z}{|r_m-r|}-\frac{z_m}{\left|r_m\right|})\end{array}\right]---\left(2\right)$

Wherein, μ _e, μ _s, μ _mBe respectively the gravitational constant of the earth, the sun, the moon, R _eBe earth radius, J ₂Be 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\mathrm{SNR}\sqrt{T}=\frac{1}{{5T}_{10}+{7T}_{50}}\sqrt{\frac{{5T}_{10}-T_{50}}{T_{50}(T_{10}-T_{50})}}---\left(3\right)$

In the formula, SNR is the signal to noise ratio (S/N ratio) of pulsar, and T is the pulsar swing circle, T ₅₀And T ₁₀Be 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

$[\mathrm{\Δ\φ}+\mathrm{\ΔN}]\·\mathrm{\λ}-\hat{n}\·r_E=\hat{n}\·r+c{\mathrm{\δt}}_{\mathrm{sat}}+{\mathrm{\Δ}}_{\mathrm{rel}}+v_k---\left(4\right)$

Wherein, Δ φ=φ _SSB-φ _Sat, Δ 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):

$\mathrm{\ΔN}=\mathrm{round}(\hat{n}\·(r_E+\tilde{r})/\mathrm{\λ}-\mathrm{\Δ\φ})---\left(5\right)$

Wherein, The satellite position predicted value that provides for the satellite dynamics equation;

GDOP value calculating method in the described step (2):

$\mathrm{GDOP}=\sqrt{\mathrm{trace}\left(C\right)}---\left(6\right)$

Wherein, trace represents to ask matrix trace, and C is the site error covariance matrix, calculate by following formula

$C=\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}\left[\begin{array}{ccc}{\mathrm{\λ}}_1^2{\mathrm{\σ}}_{\mathrm{\φ}1}^2& 0& 0\\ 0& {\mathrm{\λ}}_2^2{\mathrm{\σ}}_{\mathrm{\φ}2}^2& 0\\ 0& 0& {\mathrm{\λ}}_3^2{\mathrm{\σ}}_{\mathrm{\φ}3}^2\end{array}\right]{\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}}^T---\left(7\right)$

Wherein, H is the matrix of coefficients of measurement equation, λ _i(i=1,2,3) are the impulse radiation wavelength of i pulsar, σ _{φ i}It 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

${\mathrm{<figure-callout\; id="1"\; label="Z"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">Z</figure-callout>}}_1=\arccos (-\frac{r\&CenterDot;s}{r})+v_1---\left(8\right)$

$Z_2=\sqrt{x^2+{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">y</figure-callout>}}^2+z^2}+v_2---\left(9\right)$

Wherein, Z ₁Be the observed quantity of starlight angular distance, Z ₂For 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 ₁For the starlight angular distance is measured noise, v ₂For 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 ₀With 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 ₀+Z(x，y)(10)

Artificial landform is produced by two-dimentional single order discrete autoregressive process:

Z(x _i，y _i)＝a ₁Z(x _i-Δx，y _i)+a ₂Z(x _i，y _i-Δy)+a ₃Z(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 ²);

Isotropy and smooth performance according to landform make coefficient a ₁=a ₂=a, a ₃=b, and as a=exp (1/T _Auc), b=-a ², ${\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="w"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">w</figure-callout>}}^2={\left(d_{\mathrm{RMS}}\right)}^2\exp (-(i+j)/T_{\mathrm{auc}})$ The time, bidimensional random function Z (x in the formula (11) _i, y _i) related function be

$R_{\mathrm{ij}}=\left(d_{\mathrm{MRS}}\right)\exp (-\frac{i+j}{T_{\mathrm{auc}}})---\left(12\right)$

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 ₀, then have

Z(x _i，y ₀)＝aZ(x _i-Δx，y ₀)+W(x _i，y ₀)(13)

On the landform border of y direction, x _i=x ₀, then have

Z(x ₀，y _i)＝aZ(x ₀，y _i-Δy)+W(x ₀，y _i)(14)

By setting h ₀, 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=\left[\begin{array}{c}{u}_k\\ r_k\end{array}\right]+v_k---\left(15\right)$

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:

$\left\{\begin{array}{c}{P}_{i0}^{-1}\left(k\right)={\mathrm{\&beta;}}_i{P}_{g0}^{-1}\left(k\right)\\ {\hat{X}}_{i0}\left(k\right)={\hat{X}}_{g0}\left(k\right)\\ Q_{i0}^{-1}\left(k\right)={\mathrm{\&beta;}}_i{Q}_{g0}^{-1}\left(k\right)\end{array}\right.i=\mathrm{1,2,3,4},m---\left(17\right)$

Wherein, i=1, four subfilters of 2,3,4 expressions, m represents senior filter, the information distribution factor-beta _iSatisfy the information conservation principle: $\underset{i=1}{\overset{4,m}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i=1,$ Get β herein _m=0, β ₁=β ₂=β ₃=β ₄=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

$\left\{\begin{array}{c}{\hat{X}}_g=P_g\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{P}_i^{-1}{\hat{X}}_i\\ P_g^{-1}=P_1^{-1}+P_2^{-1}+P_3^{-1}+P_4^{-1}\end{array}\right.---\left(18\right)$

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), χ ²The 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 ²Method, 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\left(k\right)=z_i\left(k\right)-H_i\left(k\right){\hat{X}}_i(k/k-1)---\left(19\right)$

$S_i\left(k\right)=H_i\left(k\right)P_i(k/k-1)H_i^T\left(k\right)+R_i\left(k\right)---\left(20\right)$

$D_i\left(k\right)=d_i^T\left(k\right)S_i^{-1}\left(k\right)d_i\left(k\right)---\left(21\right)$

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 ₀]=α 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

$\stackrel{\&CenterDot;}{X}=f(X,t)+W\left(t\right)---\left(1\right)$

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)=\left[\begin{array}{c}{v}_x;v_y;v_z\\ -\frac{{\mathrm{\&mu;}}_{e}x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\&mu;}}_s(\frac{x_s-x}{|r_s-r|}-\frac{x_s}{\left|r_s\right|})+{\mathrm{\&mu;}}_m(\frac{x_m-x}{|r_m-r|}-\frac{x_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\&mu;}}_{e}y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\&mu;}}_s(\frac{y_s-y}{|r_s-r|}-\frac{y_s}{\left|r_s\right|})+{\mathrm{\&mu;}}_m(\frac{y_m-y}{|r_m-r|}-\frac{y_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\&mu;}}_{e}z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mathrm{\&mu;}}_s(\frac{z_s-z}{|r_s-r|}-\frac{z_s}{\left|r_s\right|})+{\mathrm{\&mu;}}_m(\frac{z_m-z}{|r_m-r|}-\frac{z_m}{\left|r_m\right|})\end{array}\right]---\left(2\right)$

Wherein, μ _e, μ _s, μ _mBe respectively the gravitational constant of the earth, the sun, the moon, R _eBe earth radius, J ₂Be 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\mathrm{SNR}\sqrt{T}=\frac{1}{{5T}_{10}+{7T}_{50}}\sqrt{\frac{{5T}_{10}-T_{50}}{T_{50}(T_{10}-T_{50})}}---\left(3\right)$

In the formula, SNR is the signal to noise ratio (S/N ratio) of pulsar, and T is the pulsar swing circle, T ₅₀And T ₁₀Be 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_sat}Represent 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

$\mathrm{\&Phi;}\left(t\right)=\mathrm{\&phi;}\left(t_0\right)+f(t-t_0)+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{f^{\left(m\right)}{(t-t_0)}^{m+1}}{(m+1)!}---\left(4\right)$

In the formula, φ (t ₀) be with reference to t epoch ₀The time the pulsar phase place, f is that pulsar is at t ₀The 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

$[\mathrm{\&Delta;\&phi;}+\mathrm{\&Delta;N}]\&CenterDot;\mathrm{\&lambda;}-\hat{n}\&CenterDot;r_E=\hat{n}\&CenterDot;r+c{\mathrm{\&delta;t}}_{\mathrm{sat}}+{\mathrm{\&Delta;}}_{\mathrm{rel}}+v_k---\left(5\right)$

Wherein, Δ φ=φ _SSB-φ _Sat, 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

$\left[\begin{array}{c}({\mathrm{\&Delta;\&phi;}}_1+{\mathrm{\&Delta;N}}_1)\&CenterDot;{\mathrm{\&lambda;}}_1-{\hat{n}}_1\&CenterDot;r_E\\ ({\mathrm{\&Delta;\&phi;}}_2+{\mathrm{\&Delta;N}}_2)\&CenterDot;{\mathrm{\&lambda;}}_2-{\hat{n}}_2\&CenterDot;r_E\\ ({\mathrm{\&Delta;\&phi;}}_3+{\mathrm{\&Delta;N}}_3)\&CenterDot;{\mathrm{\&lambda;}}_3-{\hat{n}}_3\&CenterDot;r_E\end{array}\right]=H\&CenterDot;\left[\begin{array}{c}r\\ {\mathrm{c\&delta;t}}_{\mathrm{sat}}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\&Delta;}}_{\mathrm{rel}1}\\ {\mathrm{\&Delta;}}_{\mathrm{rel}2}\\ {\mathrm{\&Delta;}}_{\mathrm{rel}3}\end{array}\right]+V_k---\left(6\right)$

Wherein $H=\left[\begin{array}{cc}{\hat{n}}_1& 1\\ {\hat{n}}_2& 1\\ {\hat{n}}_3& 1\end{array}\right].$

(4) calculate the integer ambiguity value

The satellite position predicted value that uses the satellite orbit precursor to provide

, integer ambiguity for Δ N is

$\mathrm{\&Delta;N}=\mathrm{round}(\hat{n}\&CenterDot;(r_E+\tilde{r})/\mathrm{\&lambda;}-\mathrm{\&Delta;\&phi;})---\left(7\right)$

(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

$\mathrm{cov}\left(P\right)=\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}\left[\begin{array}{ccc}{\mathrm{\&lambda;}}_1^2{\mathrm{\&sigma;}}_{\mathrm{\&phi;}1}^2& 0& 0\\ 0& {\mathrm{\&lambda;}}_2^2{\mathrm{\&sigma;}}_{\mathrm{\&phi;}2}^2& 0\\ 0& 0& {\mathrm{\&lambda;}}_3^2{\mathrm{\&sigma;}}_{\mathrm{\&phi;}3}^2\end{array}\right]{\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}}^T---\left(8\right)$

σ wherein _{φ i}(i=1,2,3) are the phase measurement accuracy of i pulsar.Then the GDOP value is

$\mathrm{GDOP}=\sqrt{\mathrm{trace}\left(C\right)}=\sqrt{{\mathrm{\&sigma;}}_x^2+{\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">y</figure-callout>}}^2+{\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">z</figure-callout>}}^2+{\mathrm{\&sigma;}}_t^2}---\left(9\right)$

Wherein, trace represents to ask matrix trace, σ _x ², δ _y ², σ _z ², σ _t ²Be 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, k}Be 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\left[W_k{W}_k^T\right]=Q\left(k\right),$ The covariance matrix of measurement noise is $E\left[V_{k+1}{V}_{k+1}^T\right]=R\left(k\right),$ 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:

$\left\{\begin{array}{c}{\hat{X}}_{k+1,k}={\hat{X}}_{k,k}+f[{\hat{X}}_{k,k},t_k]\&CenterDot;T\\ {\hat{X}}_{k+1/k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-H_{k+1}{X}_{k+1/k})\\ P_{k+1/k}={\mathrm{\&phi;}}_{k+1/k}{P}_{k/k}{\mathrm{\&phi;}}_{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\end{array}\right.---\left(12\right)$

In the formula, f represents the satellite orbit kinetics equation, and T is the Kalman filtering cycle, φ _{K+1, k}Be state-transition matrix, K _K+1Be kalman gain coefficient, H _K+1Be measurement equation matrix of coefficients, P _{K+1, k}Be optimum prediction valuation error covariance matrix, P _{K+1, k+1}Be 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

$\left\{\begin{array}{c}{X}_{k+1}=F(X_k,k)+W_k\\ Z_{k+1}=H\left(X_k\right)+V_k\end{array}\right.---\left(13\right)$

The dimension of the system state variable of setting up departments is n * 1 dimension, and 2n+1 sampled point is so

${\mathrm{\&chi;}}_0=\hat{X}$

${\mathrm{\&chi;}}_i=\hat{X}+{\left(\sqrt{(n+\mathrm{\&lambda;})P_{k/k}}\right)}_i,i=\mathrm{1,2},...,n---\left(14\right)$

${\mathrm{\&chi;}}_i=\hat{X}-{\left(\sqrt{(n+\mathrm{\&lambda;})P_{k/k}}\right)}_i,i=n+1,...,2n$

The weight of each sampled point correspondence is

$W_0^{\left(m\right)}=\mathrm{\&lambda;}/(n+\mathrm{\&lambda;})$

$W_0^{\left(c\right)}=\mathrm{\&lambda;}/(n+\mathrm{\&lambda;})+(1-{\mathrm{\&alpha;}}^2+\mathrm{\&beta;})---\left(15\right)$

$W_i^{\left(m\right)}=W_i^{\left(c\right)}=\mathrm{\&lambda;}/\left\{2(n+\mathrm{\&lambda;})\right\},i=\mathrm{1,2},...,2n$

Wherein, λ=α ²(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

${\mathrm{\&chi;}}_k=\left[\begin{array}{ccc}{\hat{X}}_k& {\hat{X}}_k+\sqrt{n+\mathrm{\&lambda;}}\sqrt{P_{k/k}}& {\hat{X}}_k-\sqrt{n+\mathrm{\&lambda;}}\sqrt{P_{k/k}}\end{array}\right]---\left(16\right)$

Time upgrades

$\left\{\begin{array}{c}{\mathrm{\&chi;}}_{k+1/k}^{*}=F({\mathrm{\&chi;}}_k,k),{\hat{X}}_{k+1/k}=\underset{0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_i^{\left(m\right)}{\mathrm{\&chi;}}_{k+1/k}^{*}\\ P_{k+1/k}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_i^{\left(c\right)}[{\mathrm{\&chi;}}_{i,k+2/k}^{*}-{\hat{X}}_{k+1/k}]{[{\mathrm{\&chi;}}_{i,k+1/k}^{*}-{\hat{X}}_{k+1/k}]}^T+{\mathrm{<figure-callout\; id="1"\; label="Q\; k\; +"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">Q</figure-callout>}}_{k+}\\ {\mathrm{\&chi;}}_{k+1/k}=\left[\begin{array}{ccc}{\hat{X}}_{k+1/k}& {\hat{X}}_{k+1/k}+\mathrm{\&gamma;}\sqrt{P_{k+1/k}}& {\hat{X}}_{k+1/k}-\mathrm{\&gamma;}\sqrt{P_{k+1/k}}\end{array}\right]\\ z_{k+1/k}=H\left({\mathrm{\&chi;}}_{k+1/k}\right),{\hat{Z}}_{k+1/k}=\underset{0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_i^{\left(m\right)}{z}_{k+1/k}\end{array}\right.0---\left(17\right)$

Measure and upgrade

$\left\{\begin{array}{c}{P}_{{\hat{Z}}_{k+1},{\hat{Z}}_{k+1}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_i^{\left(c\right)}[z_{i,k+1/k}-{\hat{Z}}_{k+1/k}]{[z_{i,k+1/k}-{\hat{Z}}_{k+1/k}]}^T+R_{k+1}}\\ P_{{\hat{X}}_{k+1},{\hat{Z}}_{k+1}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_i^{\left(c\right)}[{\mathrm{\&chi;}}_{i,k+1/k}-{\hat{X}}_{k+1/k}]{[z_{i,k+1/k}-{\hat{Z}}_{k+1/k}]}^T+R_{k+1}}\\ K_{k+1}=P_{{\hat{X}}_{k+1},{\hat{Z}}_{k+1}}{P}_{{\hat{Z}}_{k+1},{\hat{Z}}_{k+1}}^{-1},{\hat{X}}_{k+1/k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-{\hat{Z}}_{k+1/k})\\ P_{k+1/k+1}=P_{k+1/k}-K_{k+1}{P}_{{\hat{Z}}_{k+1},Z_{k+1}}{K}_{k+1}^T\end{array}\right.---\left(18\right)$

Wherein, the covariance matrix of state model noise is $E=\left[W_k{W}_k^T\right]=Q_k,$ The covariance matrix of measurement noise is $E\left[V_{k+1}{V}_{k+1}^T\right]=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=\frac{R_e}{\sin \mathrm{\&rho;}},$ Thereby apart from the measurement equation that is observed quantity be with starlight angular distance and the earth's core

${\mathrm{<figure-callout\; id="1"\; label="Z"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">Z</figure-callout>}}_1=\arccos (-\frac{r\&CenterDot;s}{r})+v_1---\left(19\right)$

${\mathrm{<figure-callout\; id="2"\; label="Z"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">Z</figure-callout>}}_2=\sqrt{x^2+{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">y</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">z</figure-callout>}}^2}+v_2---\left(20\right)$

Wherein, Z ₁Be the observed quantity of starlight angular distance, Z ₂For 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 ₁For the starlight angular distance is measured noise, v ₂For 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 ₂For secondary is reined in to such an extent that allow polynomial expression, J ₂Be humorous coefficient of earth zone, R is an earth mean radius value, parameter alpha _eFor

${\mathrm{\&alpha;}}_e=\frac{\left[\frac{2}{3}{J}_2(\frac{{\mathrm{\&alpha;}}_e}{R}+\frac{m}{2})\right]}{[1-\frac{2}{3}m]m}\&ap;0.003449963---\left(22\right)$

The measurement equation that then with the orbit altitude is observed quantity is

${\mathrm{<figure-callout\; id="2"\; label="Z"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">Z</figure-callout>}}_2=r-R[1-\frac{2{\mathrm{\&alpha;}}_e}{3}(0.5{\left(\frac{z}{r}\right)}^2-\frac{1}{2})]-h(x,y)+v_2---\left(23\right)$

In the formula, Z ₂Be the satellite orbital altitude that radar altimeter records, r is the distance of satellite to the earth's core, v ₂For 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 ₀With 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 ₀+Z(x，y)(24)

Artificial landform is produced by two-dimentional single order discrete autoregressive process:

Z(x _i，y _i)＝a ₁Z(x _i-Δx，y _i)+a ₂Z(x _i，y _i-Δy)+a ₃Z(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 ²).

Isotropy and smooth performance according to landform can make coefficient a ₁=a ₂=a, a ₃=b, and as a=exp (1/T _Auc), b=-a ², ${\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="w"\; filenames="C2007101915270022H1.png"\; state="\{\{state\}\}">w</figure-callout>}}^2={\left(d_{\mathrm{RMS}}\right)}^2\exp (-(i+j)/T_{\mathrm{auc}})$ The time, bidimensional random function Z (x in the formula (25) _i, y _i) related function be

$R_{\mathrm{ij}}=\left(d_{\mathrm{RMS}}\right)\exp (-\frac{i+j}{T_{\mathrm{auc}}})---\left(26\right)$

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 ₀, then have

Z(x _i，y ₀)＝aZ(x _i-Δx，y ₀)+W(x _i，y ₀)(27)

If on the landform border of y direction, x _i=x ₀, then have

Z(x ₀，y _i)＝aZ(x ₀，y _i-Δy)+W(x ₀，y _i)(28)

By setting h ₀, 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=\left[\begin{array}{c}{u}_k\\ r_k\end{array}\right]+v_k---\left(29\right)$

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:

$\left\{\begin{array}{c}{P}_{i0}^{-1}\left(k\right)={\mathrm{\&beta;}}_i{P}_{g0}^{-1}\left(k\right)\\ {\hat{X}}_{i0}\left(k\right)={\hat{X}}_{g0}\left(k\right)\\ Q_{i0}^{-1}\left(k\right)={\mathrm{\&beta;}}_i{Q}_{g0}^{-1}\left(k\right)\end{array}\right.(i=\mathrm{1,2,3,4},m)---\left(31\right)$

Wherein, i=1, four subfilters of 2,3,4 expressions, m represents senior filter, the information distribution factor-beta _iSatisfy the information conservation principle: $\underset{i=1}{\overset{4m}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i=1,$ Get β herein _m=0, β ₁=β ₂=β ₃=β ₄=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

$\left\{\begin{array}{c}{\hat{X}}_g=P_g\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{P}_i^{-1}{\hat{X}}_i\\ P_g^{-1}=P_1^{-1}+P_2^{-1}+P_3^{-1}+P_4^{-1}\end{array}\right.---\left(32\right)$

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 ²Method 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\left(k\right)=z_i\left(k\right)-H_i\left(k\right){\hat{X}}_i(k/k-1)---\left(33\right)$

$S_i\left(k\right)=H_i\left(k\right)P_i(k/k-1)H_i^T\left(k\right)+R_i\left(k\right)---\left(34\right)$

$D_i\left(k\right)=d_i^T\left(k\right)S_i^{-1}\left(k\right)d_i\left(k\right)---\left(35\right)$

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 ₀]=α 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 ${\hat{X}}_g={\left[\begin{array}{cccccc}\hat{x}& \hat{y}& \hat{z}& {\hat{v}}_x& {\hat{v}}_y& {\hat{v}}_z\end{array}\right]}^T$ With its covariance matrix $P_g={\left[\begin{array}{cccccc}{p}_x& p_y& p_z& p_v,& p_v,& p_v\end{array}\right]}^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 χ ²Method of inspection detects system and the isolated fault subsystem that breaks down;

In the described step (1), foundation comprises J ₂Item 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

$\stackrel{\&CenterDot;}{X}=f(X,t)+W\left(t\right)---\left(1\right)$

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)=\left[\begin{array}{c}{v}_x;v_y;v_z\\ -\frac{{\mathrm{\&mu;}}_{e}x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\&mu;}}_s(\frac{x_s-x}{|r_s-r|}-\frac{x_s}{\left|r_s\right|})+{\mathrm{\&mu;}}_m(\frac{x_m-x}{|r_m-r|}-\frac{x_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\&mu;}}_{e}y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\&mu;}}_s(\frac{y_s-y}{|r_s-r|}-\frac{y_s}{\left|r_s\right|})+{\mathrm{\&mu;}}_m(\frac{y_m-y}{|r_m-r|}-\frac{y_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\&mu;}}_{e}z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mathrm{\&mu;}}_s(\frac{z_s-z}{|r_s-r|}-\frac{z_s}{\left|r_s\right|})+{\mathrm{\&mu;}}_m(\frac{z_m-z}{|r_m-r|}-\frac{z_m}{\left|r_m\right|})\end{array}\right]---\left(2\right)$

Wherein, μ _e, μ _s, μ _mBe respectively the gravitational constant of the earth, the sun, the moon, R _eBe earth radius, J ₂Be 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\mathrm{SNR}\sqrt{T}\frac{1}{{5T}_{10}+7T_{50}}\sqrt{\frac{{5T}_{10}-T_{50}}{T_{50}(T_{10}-T_{50})}}---\left(3\right)$

In the formula, SNR is the signal to noise ratio (S/N ratio) of pulsar, and T is the pulsar swing circle, T ₅₀And T ₁₀Be 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

$[\mathrm{\&Delta;\&phi;}+\mathrm{\&Delta;N}]\&CenterDot;\mathrm{\&lambda;}-\hat{n}\&CenterDot;r_E=\hat{n}\&CenterDot;r+c{\mathrm{\&delta;t}}_{\mathrm{sat}}+{\mathrm{\&Delta;}}_{\mathrm{rel}}+v_k---\left(4\right)$

Wherein, Δ φ=φ _SSB-φ _Sat, Δ 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):

$\mathrm{\&Delta;N}=\mathrm{round}(\hat{n}\&CenterDot;(r_E+\tilde{r})/\mathrm{\&lambda;}-\mathrm{\&Delta;\&phi;})---\left(5\right)$

Wherein,

The satellite position predicted value that provides for the satellite dynamics equation;

GDOP value calculating method in the described step (2):

$\mathrm{GDOP}=\sqrt{\mathrm{trace}\left(C\right)}---\left(6\right)$

Wherein, trace represents to ask matrix trace, and C is the site error covariance matrix, calculate by following formula

$C=\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}\left[\begin{array}{ccc}{\mathrm{\&lambda;}}_1^2{\mathrm{\&sigma;}}_{\mathrm{\&phi;}1}^2& 0& 0\\ 0& {\mathrm{\&lambda;}}_2^2{\mathrm{\&sigma;}}_{\mathrm{\&phi;}2}^2& 0\\ 0& 0& {\mathrm{\&lambda;}}_3^2{\mathrm{\&sigma;}}_{\mathrm{\&phi;}3}^2\end{array}\right]{\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}}^T---\left(7\right)$

Wherein, H is the matrix of coefficients of measurement equation, λ _i, i=1,2,3 is the impulse radiation wavelength of i pulsar, σ _{φ i}It 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 (-\frac{r\&CenterDot;s}{r})+v_1---\left(8\right)$

$Z_2=\sqrt{x^2+y^2+z^2}+v_2---\left(9\right)$

Wherein, Z ₁Be the observed quantity of starlight angular distance, Z ₂For 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 ₁For the starlight angular distance is measured noise, v ₂For 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 ₀With 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 ₀+Z(x，y) (10)

Artificial landform is produced by two-dimentional single order discrete autoregressive process:

Z(x _i，y _i)＝a ₁Z(x _i-Δx，y _i)+a ₂Z(x _i，y _i-Δy)+a ₃Z(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 ²);

Isotropy and smooth performance according to landform make coefficient a ₁=a ₂=a, a ₃=b, and as a=exp (1/T _Auc), b=-a ², ${\mathrm{\&sigma;}}_w^2={\left(d_{\mathrm{RMS}}\right)}^2\exp (-(i+j)/T_{\mathrm{auc}})$ The time, bidimensional random function Z (x in the formula (11) _i, y _i) related function be

$R_{\mathrm{ij}}=\left(d_{\mathrm{RMS}}\right)\exp (-\frac{i+j}{T_{\mathrm{auc}}})---\left(12\right)$

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 ₀, then have

Z(x _i，y ₀)＝aZ(x _i-Δx，y ₀)+W(x _i，y ₀) (13)

On the landform border of y direction, x _i=x ₀, then have

Z(x ₀，y _i)＝aZ(x ₀，y _i-Δy)+W(x ₀，y _i) (14)

By setting h ₀, 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=\left[\begin{array}{c}{u}_k\\ r_k\end{array}\right]+v_k---\left(15\right)$

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:

$\left\{\begin{array}{c}{P}_{i0}^{-1}\left(k\right)={\mathrm{\&beta;}}_i{P}_{g0}^{-1}\left(k\right)\\ {\hat{X}}_{i0}\left(k\right)={\hat{X}}_{g0}\left(k\right)\\ Q_{10}^{-1}\left(k\right)={\mathrm{\&beta;}}_i{Q}_{g0}^{-1}\left(k\right)\end{array}\right.,i=\mathrm{1,2,3,4},m---\left(17\right)$

Wherein, i=1, four subfilters of 2,3,4 expressions, m represents senior filter, the information distribution factor-beta _iSatisfy the information conservation principle: $\underset{i=1}{\overset{4,m}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i=1,$ Get β herein _m=0, β ₁=β ₂=β ₃=β ₄=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

$\left\{\begin{array}{c}{\hat{X}}_g=P_g\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{P}_i^{-1}{\hat{X}}_i\\ P_g^{-1}=P_1^{-1}+P_2^{-1}+P_3^{-1}+P_4^{-1}\end{array}\right.---\left(18\right)$

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), χ ²The 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 ²Method, 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\left(k\right)=z_i\left(k\right)-H_i\left(k\right){\hat{X}}_i(k/k-1)---\left(19\right)$

$S_i\left(k\right)=H_i\left(k\right)P_i(k/k-1)H_i^T\left(k\right)+R_i\left(k\right)---\left(20\right)$

$D_i\left(k\right)=d_i^T\left(k\right)S_i^{-1}\left(k\right)d_i\left(k\right)---\left(21\right)$

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 ₀]=α solves thresholding T _D

Images