Method for achieving space-based autonomous navigation of global navigation satellite system (GNSS) satellites

Abstract

Disclosed is a method for achieving space-based autonomous navigation of global navigation satellite system (GNSS) satellites, and relates to the field of satellite navigation technologies. The method includes the following steps: optimizing a DRO, and establishing a dynamic model of an earth-moon space satellite orbit; establishing measurement links, by a low earth orbit (LEO) data relay satellite, with an earth-moon space DRO satellite and a GNSS respectively, and measuring an inter-satellite distance for modeling and linearization; adopting an extended Kalman filter (EKF) method to process inter-satellite measurement data, and autonomously determining a position and velocity of the global navigation satellite system without depending on the ground measurement and control support.

Claims

1. A method for achieving space-based autonomous navigation of global navigation satellite system (GNSS) satellites, comprising the following steps: S1: optimizing an earth-moon space DRO, and establishing a dynamic model of an earth-moon space satellite orbit; S2: establishing measurement links, by a low earth orbit (LEO) data relay satellite, with an earth-moon space DRO satellite and a GNSS satellite respectively, and measuring an inter-satellite distance for modeling and linearization; and S3: adopting a data processing method to process inter-satellite measurement data, and achieving autonomous navigation of the GNSS satellite without depending on the ground measurement and control support.

2. The method for achieving autonomous navigation of GNSS satellites according to claim 1, wherein in step S1, the optimizing an earth-moon space DRO specifically comprises: A1: deducing an earth-moon DRO family with periodic characteristics from a circular restricted three-body problem (CRTBP) model; and A2: calculating gravity asymmetry of each orbit, and selecting an orbit with high gravity asymmetry as a preferred earth-moon DRO.

3. The method for achieving autonomous navigation of GNSS satellites according to claim 1, wherein in step S1, the establishing a dynamic model of an earth-moon space satellite orbit specifically comprises: B 1: establishing satellite state vectors: at any time t, y.sub.D.sup.j(t), j.sub.M.sup.k(t) and y.sub.T.sup.l(t) denote position and velocity state vectors of a j-th DRO satellite, a k-th navigation satellite and an l-th LEO data relay satellite, respectively, namely y.sub.D.sup.j(t)=[r.sub.D.sup.j(t).sup.T,v.sub.D.sup.j(t).sup.T].sup.T, y.sub.M.sup.j(t)=[r.sub.M.sup.j(t).sup.T,v.sub.M.sup.j(t).sup.T].sup.T and y.sub.T.sup.j(t)=[r.sub.T.sup.j(t).sup.T,v.sub.T.sup.j(t).sup.T].sup.T; wherein the numbers of satellites participating in orbit determination are J, K and L respectively, j=0, . . . ,J−1, k=0, . . . ,K−1 and l=0, . . . ,L−1 , and the state vector of each type of satellites is denoted as:
y.sub.D(t)=[y.sub.D.sup.0(t).sup.T, . . . ,y.sub.D.sup.J−1(t).sup.T)].sup.T   (1) and
y.sub.M(t)=[y.sub.M.sup.0(t).sup.T, . . . ,y.sub.M.sup.K−1(t).sup.T)].sup.T   (2) and
y.sub.T(t)=[y.sub.T.sup.0(t).sup.T, . . . ,y.sub.T.sup.L−1(t).sup.T)].sup.T   (3) which are further combined into a state vector:
X(t)=[y.sub.D(t).sup.T,y.sub.M(t).sup.T,y.sub.T(t).sup.T].sup.T   (4); B2: a differential equation of an orbital motion state of any satellite is: $\begin{array}{cc}\left[\begin{array}{c}\stackrel{.}{r}\\ \stackrel{.}{v}\end{array}\right]=\left[\begin{array}{c}v\\ a_E+a_M+a_S+a_{\mathrm{SRP}}+a_D\end{array}\right]& \left(5\right)\end{array}$ wherein a.sub.E and a.sub.M are spherical harmonic gravity acceleration models of the earth and the moon respectively, a.sub.S is a gravity acceleration of the solar particle, a.sub.SRP is an acceleration of solar radiation light pressure, and a.sub.D is mainly an atmospheric drag acceleration suffered by an LEO satellite; and B3: calculating a state transition matrix Φ.sub.i□Φ(t.sub.i,t.sub.i−1)=∂X(t.sub.i)/∂X(t.sub.i−1) from any time t.sub.i−1 to any time t.sub.i.

4. The method for achieving autonomous navigation of GNSS satellites according to claim 3, wherein in step B2, a first-order dynamic differential equation used as a basis is:
{dot over (X)}(t)=f(X(t),t)   (6) and in step B3, the state transition matrix for state calculation follows a variational equation: $\begin{array}{cc}{\stackrel{.}{\Phi }}_i=\frac{\partial f\left(t_i,X\left(t_i\right)\right)}{\partial X\left(t_i\right)}{\Phi }_i.& \left(7\right)\end{array}$

5. The method for achieving autonomous navigation of GNSS satellites according to claim 1, wherein step S2 specifically comprises: S21: at any time t.sub.i, an inter-satellite measurement model of an l-th LEO satellite and a j-th DRO satellite and a k-th navigation satellite being as follows: $\begin{array}{cc}\{\begin{array}{c}{z}_{\mathrm{DT},i}^{jl}=h_{DT}^{jl}\left(y_D^j\left(t_i\right),y_T^l\left(t_i\right)\right)+{\epsilon }_{DL}=\left|r_D^j\left(t_i\right)-r_T^l\left(t_i\right)\right|+{\epsilon }_{DL}\\ z_{\mathrm{MT},i}^{kl}=h_{MT}^{kl}\left(y_M^k\left(t_i\right),y_T^l\left(t_i\right)\right)+{\epsilon }_{ML}=\left|r_M^k\left(t_i\right)-r_T^l\left(t_i\right)\right|+{\epsilon }_{ML}\end{array}& \left(8\right)\end{array}$ wherein z is a measured value, h is a model value, and ε is a measurement noise; S22: X.sub.i.sup.−[(y.sub.D,i.sup.−).sup.T,(y.sub.M,i.sup.−).sup.T,(y.sub.T,i.sup.−).sup.T].sup.T being a known initial state of a satellite, and linearizing the measurement model (8) near the state: $\begin{array}{cc}{z}_{\mathrm{DT},i}^{jl}=h_{DT}^{jl}\left(y_{D,i}^{j,-},y_{T,i}^{l,-}\right)+H_{\mathrm{DD},i}^{\mathrm{jl}}.Math.\Delta y_{D,i}^j+H_{\mathrm{DT},i}^{\mathrm{jl}}.Math.\Delta y_{T,i}^l+{\epsilon }_{\mathrm{DT}}& \left(9\right)\end{array}$ $\begin{array}{cc}{z}_{\mathrm{MT},i}^{\mathrm{kl}}=h_{\mathrm{MT}}^{\mathrm{kl}}\left(y_{M,i}^{k,-},y_{T,i}^{j,-}\right)+H_{\mathrm{MM},i}^{\mathrm{kl}}.Math.\Delta y_{M,i}^k+H_{\mathrm{MT},i}^{\mathrm{kl}}.Math.\Delta y_{T,i}^l+{\epsilon }_{\mathrm{MT}}& \end{array}$ $\begin{array}{cc}\mathrm{wherein}& \end{array}$ $\begin{array}{cc}{H}_{\mathrm{DD},i}^{\mathrm{jl}}▯\frac{\partial h_{\mathrm{DT},i}^{\mathrm{jl}}}{\partial y_{D,i}^{j,-}},H_{\mathrm{DT},i}^{\mathrm{jl}}▯\frac{\partial h_{\mathrm{DT},i}^{\mathrm{jl}}}{\partial y_{T,i}^{l,-}},H_{\mathrm{MM},i}^{\mathrm{kl}}▯\frac{\partial h_{\mathrm{MT},i}^{\mathrm{kl}}}{\partial y_{M,i}^{k,-}},H_{\mathrm{MT},i}^{\mathrm{kl}}▯\frac{\partial h_{\mathrm{MT},i}^{\mathrm{kl}}}{\partial y_{T,i}^{l,-}},& \end{array}$  and there are relations H.sub.DT,i.sup.jl=−H.sub.DD,i.sup.jl and H.sub.MT,i.sup.kl=−H.sub.MM,i.sup.kl; ΔX.sub.i=[(Δy.sub.D,i.sup.−).sup.T,(Δy.sub.M,i.sup.−).sup.T,(Δy.sub.T,i.sup.−).sup.T].sup.T is set to be an increment of a to-be-estimated state with respect to an initial value state, and then ΔX.sub.i=X.sub.i−X.sub.i.sup.−; S23: combining all measured values into z.sub.i=[z.sub.DT,i.sup.T,z.sub.MT,i.sup.T].sup.T, with a component denoted as: $\begin{array}{cc}{z}_{\mathrm{DT},i}={\left[z_{\mathrm{DT},i}^{00},.Math.,z_{\mathrm{DT},i}^{J-1,0},.Math.,z_{\mathrm{DT},i}^{0,L-1},.Math.,z_{\mathrm{DT},1}^{J-1,L-1}\right]}^T{z}_{\mathrm{MT},i}={\left[z_{\mathrm{MT},i}^{00},.Math.,z_{\mathrm{MT},i}^{K-1,0},.Math.,z_{\mathrm{MT},i}^{0,L-1},.Math.,z_{\mathrm{MT},i}^{K-1,L-1}\right]}^T& \left(10\right)\end{array}$ wherein a corresponding model value h.sub.i(X.sub.i)=[h.sub.DT,i.sup.T,h.sub.MT,i.sup.T].sup.T, and a measurement error is ε=[(ε.sub.DL).sup.T,(ε.sub.ML).sup.T]; standard deviations of microwave measurement and laser measurement errors are σ.sub.DL and σ.sub.ML respectively; and S24: at the time t.sub.i, a linearization measurement model comprising all measured values is: $\begin{array}{cc}\left[\begin{array}{c}{z}_{\mathrm{DT},i}\\ z_{\mathrm{MT},i}\end{array}\right]=\left[\begin{array}{c}{h}_{\mathrm{DT},i}\left(X_i^{-}\right)\\ h_{\mathrm{MT},i}\left(X_i^{-}\right)\end{array}\right]+\left[\begin{array}{ccc}{H}_{\mathrm{DD},i}& 0& H_{\mathrm{DT},i}\\ 0& H_{\mathrm{MM},i}& H_{\mathrm{MT},i}\end{array}\right].Math.\left[\begin{array}{c}\Delta y_{D,i}\\ \Delta y_{M,i}\\ \Delta y_{T,i}\end{array}\right]+\left[\begin{array}{c}{\epsilon }_D\\ {\epsilon }_M\end{array}\right]& \left(11\right)\end{array}$ which is further simplified as:
z.sub.i=h.sub.i(X.sub.i.sup.−)+H.sub.i.Math.ΔX.sub.i+ε  (12)

6. The method for achieving autonomous navigation of GNSS satellites according to claim 1, wherein the data processing method in step S3 is an extended Kalman filter (EKF) method, which specifically comprises: S31: adding a process noise u(t) to a dynamic equation:
{dot over (X)}(t)=f(X(t),u(t),t)  (13) wherein .sup.u(.sup.t) is a white Gaussian noise with zero mean; correspondingly, a process noise conversion matrix is: $\begin{array}{cc}\Gamma \left(t_i,t_{i-1}\right)=\frac{\partial X\left(t_i\right)}{\partial u\left(t_{i-1}\right)};& \left(14\right)\end{array}$  and S32: given an initial state and time of the satellite, processing, by an extended Kalman filter (EKF), measurement data based on continuous measurement epochs, and updating a to-be-estimated state.

7. The method for achieving autonomous navigation of GNSS satellites according to claim 6, wherein step S32 specifically comprises: A. time update wherein a satellite orbit state and a state estimation covariance are advanced from the previous time t.sub.i−1 to the current observation time t.sub.i:
X.sub.i.sup.−=X(t.sub.i;X(t.sub.i−1)=X.sub.i−1.sup.+)+Γ(t.sub.i,t.sub.i−1)u(t.sub.i−1)
P.sub.i.sup.−=Φ(t.sub.i,t.sub.i−)P.sub.i−1.sup.+Φ.sup.T(t.sub.i,t.sub.i−1)+Γ(t.sub.i,t.sub.i−1)Q.sub.cΓ.sup.T(t.sub.i,t.sub.i−1)   (15) wherein P is an orbit state covariance matrix, Q.sub.c is a process noise covariance matrix and is a (3J+3K+3L)-dimensional diagonal matrix, a diagonal value is related to dynamic model accuracy of each satellite, and a variance of an acceleration error is set as σ.sub.c.sup.2; a format of the process noise conversion matrix Γ is: $\begin{array}{cc}\Gamma =\left[\begin{array}{cccc}{\Gamma }_0& & & 0\\ & {\Gamma }_1& & \\ & & \ddots & \\ 0& & & {\Gamma }_{J+K+L‐1}\end{array}\right]& \left(16\right)\end{array}$ wherein a diagonal submatrix is denoted as: $\begin{array}{cc}{\Gamma }_k={\left[\begin{array}{c}\Delta t^2/2\times I_{3\times 3}\\ \Delta t\times I_{3\times 3}\end{array}\right]}_{6\times 3}& \left(17\right)\end{array}$ where Δt=t.sub.i−t.sub.i−1, I is an identity matrix, 0 is a zero matrix, and k=0,1, . . . ,J+K+L−1. B. measurement update where the solution of a satellite orbit state and update of an orbit state covariance are completed by using measurement data:
K.sub.i=P.sub.i.sup.−H.sub.i.sup.T(R.sub.i+H.sub.iP.sub.i.sup.−H.sub.i.sup.T).sup.−1
X.sub.i.sup.+=X.sub.i.sup.−+K.sub.i(z.sub.i−h(X.sub.i.sup.−))
P.sub.i.sup.+=(I−K.sub.iH.sub.i)P.sub.i.sup.−(I−K.sub.iH.sub.i).sup.T+K.sub.iR.sub.iK.sub.i.sup.T   (18) wherein K.sub.i is a Kalman gain, H.sub.i is a design matrix and is calculated by equation (12), I is an identity matrix, R.sub.i is a measuring noise covariance diagonal matrix, a diagonal element value corresponding to laser measurement data is σ.sub.DL.sup.2, and a diagonal element value corresponding to microwave measurement data is σ.sub.ML.sup.2.

8. The method for achieving autonomous navigation of GNSS satellites according to claim 1, wherein the method is applicable to all navigation systems, comprising at least one of a GPS system, a Beidou system, a Glonass system, a Galileo system and a quasi-zenith navigation system.

9. The method for achieving autonomous navigation of GNSS satellites according to claim 2, wherein the method is applicable to all navigation systems, comprising at least one of a GPS system, a Beidou system, a Glonass system, a Galileo system and a quasi-zenith navigation system.

10. The method for achieving autonomous navigation of GNSS satellites according to claim 3, wherein the method is applicable to all navigation systems, comprising at least one of a GPS system, a Beidou system, a Glonass system, a Galileo system and a quasi-zenith navigation system.

11. The method for achieving autonomous navigation of GNSS satellites according to claim 4, wherein the method is applicable to all navigation systems, comprising at least one of a GPS system, a Beidou system, a Glonass system, a Galileo system and a quasi-zenith navigation system.

12. The method for achieving autonomous navigation of GNSS satellites according to claim 5, wherein the method is applicable to all navigation systems, comprising at least one of a GPS system, a Beidou system, a Glonass system, a Galileo system and a quasi-zenith navigation system.

13. The method for achieving autonomous navigation of GNSS satellites according to claim 6, wherein the method is applicable to all navigation systems, comprising at least one of a GPS system, a Beidou system, a Glonass system, a Galileo system and a quasi-zenith navigation system.

14. The method for achieving autonomous navigation of GNSS satellites according to claim 7, wherein the method is applicable to all navigation systems, comprising at least one of a GPS system, a Beidou system, a Glonass system, a Galileo system and a quasi-zenith navigation system.

Description

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

(1) The foregoing summary, as well as the following detailed description of the invention, will be better understood when read in conjunction with the appended drawings. For the purpose of illustrating the invention, there are shown in the drawings embodiments which are presently preferred. It should be understood, however, that the invention is not limited to the precise arrangements and instrumentalities shown.

(2) In the drawings:

(3) FIG. 1 is a schematic flowchart of implementation of autonomous navigation of a navigation constellation;

(4) FIG. 2 is a schematic diagram of space-based measurement and control autonomous navigation of a GNSS constellation based on a DRO and an LEO data relay orbit;

(5) FIG. 3 is a DRO family in a (CRTBP) model;

(6) FIG. 4 is a preferred earth-moon long-distance DRO in Embodiment 1;

(7) FIG. 5 shows spatial distribution of gravity asymmetry, where a dashed line is a periodic orbit under the CRTBP model (the amplitude is 49000 km), and a solid line is an orbit under a real ephemeris of the solar system;

(8) FIG. 6 shows stability analysis of a DRO under different ephemeris models of the solar system;

(9) FIG. 7 shows position filtering uncertainty of a GPS satellite (PRN 1, a), a DRO satellite (b) and an LEO satellite (c), indicating that X/Y/Z is the X axis, Y axis and Z axis of coordinates respectively in a J2000 coordinate system;

(10) FIG. 8 shows position residuals of the GPS satellite (PRN 1, a), the DRO satellite (b) and the LEO satellite (c), indicating changes with time in the J2000 coordinate system, where X/Y/Z is the X axis, Y axis and Z axis of coordinates respectively;

(11) FIG. 9 shows position filtering accuracy of GPS constellation satellites when inter-satellite measurement accuracy is 0.01 m; and

(12) FIG. 10 shows position filtering accuracy of GPS constellation satellites when inter-satellite measurement accuracy is 0.1 m.

DETAILED DESCRIPTION OF THE INVENTION

(13) To make the objective, technical solutions and advantages of the present invention clearer, the following further describes the present invention in detail with reference to the accompanying drawings. It should be understood that the specific implementations described herein are merely illustrative of the present invention and are not intended to limit the present invention.

(14) Explanation of professional terms:

(15) DRO family: Earth-moon DRO family

(16) GNSS : Global Navigation Satellite System

(17) Gravity asymmetry: It is defined as the ratio of the amplitude of the gravity acceleration of a third body received at any spatial position to the sum of all accelerations received at the position. When near the earth space, the moon is the third body; and in contrast, when near the moon space, the earth is classified as the third body.

Embodiment 1

(18) This embodiment provides a method for achieving GNSS-based autonomous navigation of a space segment satellite, as shown in FIG. 1, including the following steps.

(19) S1: Optimize an earth-moon space DRO, and establish a dynamic model of an earth-moon space satellite orbit.

(20) S2: Establish measurement links, by an LEO data relay satellite, with an earth-moon space DRO satellite and a GNSS satellite respectively, and measure an inter-satellite distance for modeling and linearization.

(21) S3: Adopt an extended Kalman filter (EKF) method to process inter-satellite measurement data, and autonomously determine a position and velocity of the GNSS without depending on the ground measurement and control support.

(22) Specifically, earth-moon DROs with periodic characteristics are deduced from a CRTBP model. DROs are usually expressed in an earth-moon rotating coordinate system, the origin is the earth-moon center of mass, the X axis points from the center of the earth to the center of the moon, the Y axis is perpendicular to the orbit plane of the moon around the earth, and the Z axis completes the right-handed coordinate system, as shown in FIGS. 3 and 4. FIG. 3 shows a family of DROs, and the amplitude in the orbit plane ranges from 49000 km to 149000 km. The amplitude in the orbit plane is defined as the distance from the center of mass of the moon to the DRO (pointing in the direction of the center of the earth). The amplitude of the DRO family outside the orbit plane is generally relatively small. Arrows indicate directions of motion of DROs, which are opposite to the direction of motion of the moon. DU is the average distance between the earth and the moon, which is 384000 km. FIG. 4 shows a DRO with an amplitude of 49000 km in an earth-centered inertia frame, which is a resonance orbit close to 3:1. That is, when the moon orbits the earth once, the DROs the moon three times with an orbital period of about 8 days.

(23) In fact, some orbits in the DRO family are also located in regions with strong asymmetry of the earth-moon three-body gravity field. Hill gave a concept of gravity field asymmetry to evaluate the degree of asymmetry of the gravity field in different spatial positions. The gravity field asymmetry is defined as the ratio of the amplitude of the acceleration of a third body received at any spatial position to the sum of all accelerations received at the position. When near the earth space, the moon's gravity is classified as the third body; and in contrast, when near the moon space, the earth's gravity is classified as the third body. FIG. 5 shows a diagram of gravity field asymmetry distribution of the moon's orbit plane (namely the orbit plane of the moon around the earth). The darker the color is, the higher the gravity field asymmetry is, and the stronger the asymmetry of the field asymmetry is, the better the autonomous determination ability of the orbit in the period is. The solid line is the DRO with an amplitude of 49000 km in the orbit plane, which can be seen to be in the region with relatively strong asymmetry of the gravitational field, and is the preferred orbit for this method. The ideal DRO is deduced from a CRTBP model. However, under the real solar system ephemeris, the distance between the earth and the moon is time-varying. Besides the influence of the earth-moon gravity, the distance between the earth and the moon is also influenced by the perturbation of solar gravity and the solar radiation light pressure. Under the real ephemeris, the DROs do not coincide strictly any more in each period, as shown by dashed lines in FIG. 5.

(24) If a satellite is deployed in the DRO, it is connected to an earth GNSS satellite through an inter-satellite link, which is expected to realize autonomous navigation of the GNSS satellite without depending on ground measurement and control. Compared with a halo orbit, the DRO has the greatest advantage of long-term stability. As shown in FIG. 6, Bezrouk gave the stability analysis of DROs with different periods on the 500-year time scale, and the orbit in a region covered by solid lines is stable. Even if all perturbative forces in the solar system are taken into account (the seventh case), the DRO used herein (with a period of about 8 days) still has a stability of 500 years and is qualified to provide long-term navigation information for GNSS satellites.

(25) Specifically, the establishing a dynamic model of an orbit in step S1 includes the following steps.

(26) B 1: At any time t, y.sub.D.sup.j(t), y.sub.M.sup.k(t) and y.sub.T.sup.l(t) denote position and velocity state vectors of a j-th DRO satellite, a k-th navigation satellite and an l-th LEO data relay satellite, respectively, namely and y.sub.D.sup.j(t)=[r.sub.D.sup.j(t).sup.T,v.sub.D.sup.j(t).sup.T].sup.T, y.sub.M.sup.j(t)=[r.sub.M.sup.j(t).sup.T,v.sub.M.sup.j(t).sup.T].sup.T) and y.sub.T.sup.j(t)=[r.sub.T.sup.j(t).sup.T,v.sub.T.sup.j(t).sup.T].sup.T.

(27) The numbers of satellites participating in orbit determination are J, K and L respectively, j=0, . . . ,J−1, k=0, . . . ,K−1 and l=0, . . . ,L−1, and the state vector of each type of satellites is denoted as:
y.sub.D(t)=[y.sub.D.sup.0(t).sup.T, . . . ,y.sub.D.sup.J−1(t).sup.T)].sup.T

(28) and
y.sub.M(t)=[y.sub.M.sup.0(t).sup.T, . . . ,y.sub.M.sup.K−1(t).sup.T)].sup.T

(29) and
y.sub.T(t)=[y.sub.T.sup.0(t).sup.T, . . . ,y.sub.T.sup.L−1(t).sup.T)].sup.T

(30) which are further combined into a state vector:
X(t)=[y.sub.D(t).sup.T,y.sub.M(t).sup.T,y.sub.T(t).sup.T].sup.T

(31) An orbit state of a satellite is extrapolated according to a first-order dynamic differential equation:
{dot over (X)}(t)=f(X(t),t)

(32) B2: Extrapolate a position and velocity state of any satellite of (5) according to a first-order dynamic differential equation, where the differential equation is denoted as:

(33) $\left[\begin{array}{c}\stackrel{.}{r}\\ \stackrel{.}{v}\end{array}\right]=\left[\begin{array}{c}v\\ a_E+a_M+a_S+a_{\mathrm{SRP}}+a_D\end{array}\right]$

(34) where a.sub.E and a.sub.M are spherical harmonic gravity acceleration models of the earth and the moon respectively, a.sub.S is a gravity acceleration of the solar particle, a.sub.SRP is an acceleration of solar radiation light pressure, and a.sub.D is mainly an atmospheric drag acceleration suffered by an LEO satellite.

(35) B3: The calculation of a state transition matrix Φ.sub.i□Φ(t.sub.i,t.sub.i−1)=∂X(t.sub.i)/∂X(t.sub.i−1) from any time t.sub.i−1 to any time t.sub.i follows a variational equation:

(36) ${\stackrel{.}{\Phi }}_i=\frac{\partial f\left(t_i,X\left(t_i\right)\right)}{\partial X\left(t_i\right)}{\Phi }_i$

(37) The deduction of the satellite orbit state and the state transition matrix is performed by numerical integration according to equations (6) and (7).

(38) Specifically, step S2 specifically includes:

(39) S21: At any time t.sub.i, an inter-satellite measurement model of an l-th LEO satellite and a j-th DRO satellite and a k-th navigation satellite is as follows:

(40) $\{\begin{array}{c}{z}_{\mathrm{DT},i}^{\mathrm{jl}}=h_{\mathrm{DT}}^{\mathrm{jl}}\left(y_D^j\left(t_i\right),y_T^l\left(t_i\right)\right)+{\epsilon }_{\mathrm{DL}}=.Math.r_D^j\left(t_i\right)-r_T^l\left(t_i\right).Math.+{\epsilon }_{\mathrm{DL}}\\ z_{\mathrm{MT},i}^{\mathrm{kl}}=h_{\mathrm{MT}}^{\mathrm{kl}}\left(y_M^k\left(t_i\right),y_T^l\left(t_i\right)\right)+{\epsilon }_{\mathrm{ML}}=.Math.r_M^k\left(t_i\right)-r_T^l\left(t_i\right).Math.+{\epsilon }_{\mathrm{ML}}\end{array}$

(41) where z is a measured value, h is a model value, and ε is a measurement noise;

(42) S22: X.sub.i.sup.−[(y.sub.D,i.sup.−).sup.T,(y.sub.M,i.sup.−).sup.T,(y.sub.T,i.sup.−).sup.T].sup.T being a known initial state of a satellite, and linearizing the measurement model (8) near the state:

(43) $z_{\mathrm{DT},i}^{\mathrm{jl}}=h_{\mathrm{DT}}^{\mathrm{jl}}\left(y_{D,i}^{j,-},y_{T,i}^{l,-}\right)+H_{\mathrm{DD},i}^{\mathrm{jl}}.Math.\Delta y_{D,i}^j+H_{\mathrm{DT},i}^{\mathrm{jl}}.Math.\Delta y_{T,i}^l+{\epsilon }_{\mathrm{DT}}{z}_{\mathrm{MT},i}^{\mathrm{kl}}=h_{\mathrm{MT}}^{\mathrm{kl}}\left(y_{M,i}^{k,-},y_{T,i}^{j,-}\right)+H_{\mathrm{MM},i}^{\mathrm{kl}}.Math.\Delta y_{M,i}^k+H_{\mathrm{MT},i}^{\mathrm{kl}}.Math.\Delta y_{T,i}^l+{\epsilon }_{\mathrm{MT}}\mathrm{where}{H}_{\mathrm{DD},i}^{\mathrm{jl}}▯\frac{\partial h_{\mathrm{DT},i}^{\mathrm{jl}}}{\partial y_{D,i}^{j,-}},H_{\mathrm{DT},i}^{\mathrm{jl}}▯\frac{\partial h_{\mathrm{DT},i}^{\mathrm{jl}}}{\partial y_{T,i}^{l,-}},H_{\mathrm{MM},i}^{\mathrm{kl}}▯\frac{\partial h_{\mathrm{MT},i}^{\mathrm{kl}}}{\partial y_{M,i}^{k,-}},H_{\mathrm{MT},i}^{\mathrm{kl}}▯\frac{\partial h_{\mathrm{MT},i}^{\mathrm{kl}}}{\partial y_{T,i}^{l,-}},$
and there are relations H.sub.DT,i.sup.jl=−H.sub.DD,i.sup.jl and H.sub.MT,i.sup.kl=−H.sub.MM,i.sup.kl;

(44) ΔX.sub.i=[(Δy.sub.D,i.sup.−).sup.T,(Δy.sub.M,i.sup.−).sup.T,(Δy.sub.T,i.sup.−).sup.T].sup.T is set to be an increment of a to-be-estimated state with respect to an initial value state, and then ΔX.sub.i=X.sub.i−X.sub.i.sup.−;

(45) S23: combining all measured values into z.sub.i=[z.sub.DT,i.sup.T,z.sub.MT,i.sup.T].sup.T, with a component denoted as:

(46) $z_{\mathrm{DT},i}={\left[z_{\mathrm{DT},i}^{00},.Math.,z_{\mathrm{DT},i}^{J-1,0},.Math.,z_{\mathrm{DT},i}^{0,L-1},.Math.,z_{\mathrm{DT},1}^{J-1,L-1}\right]}^T{z}_{\mathrm{MT},i}={\left[z_{\mathrm{MT},i}^{00},.Math.,z_{\mathrm{MT},i}^{K-1,0},.Math.,z_{\mathrm{MT},i}^{0,L-1},.Math.,z_{\mathrm{MT},i}^{K-1,L-1}\right]}^T$

(47) where a corresponding model value is h.sub.i(X.sub.i)=[h.sub.DT,i.sup.T, h.sub.MT,i.sup.T].sup.T, and a measurement error is ε=[(ε.sub.DL).sup.T,(ε.sub.ML).sup.T]; standard deviations of microwave measurement and laser measurement errors are σ.sub.DL and σ.sub.ML in respectively; and

(48) S24: at the time t.sub.i, a linearization measurement model including all measured values is:

(49) $\left[\begin{array}{c}{z}_{\mathrm{DT},i}\\ z_{\mathrm{MT},i}\end{array}\right]=\left[\begin{array}{c}{h}_{\mathrm{DT},i}\left(X_i^{-}\right)\\ h_{\mathrm{MT},i}\left(X_i^{-}\right)\end{array}\right]+\left[\begin{array}{ccc}{H}_{\mathrm{DD},i}& 0& H_{\mathrm{DT},i}\\ 0& H_{\mathrm{MM},i}& H_{\mathrm{MT},i}\end{array}\right].Math.\left[\begin{array}{c}\Delta y_{D,i}\\ \Delta y_{M,i}\\ \Delta y_{T,i}\end{array}\right]+\left[\begin{array}{c}{\epsilon }_D\\ {\epsilon }_M\end{array}\right]$

(50) which is further simplified as:
z.sub.i=h.sub.i(X.sub.i.sup.−)+H.sub.i.Math.ΔX.sub.i+ε

(51) The technical solution of this patent is that LEO data relay satellites receive measurement signals of DRO satellites and navigation satellites respectively. The measurement model can be simplified to an instantaneous position measurement value between any two satellites without considering signal propagation time, and the signal propagation time can be calculated through an iterative process.

(52) The data processing method in this embodiment is an EKF method. In the EKF method, in order to compensate for a dynamic modeling error, a process noise u(t) is usually added to a dynamic equation:
{dot over (X)}(t)=f(X(t),u(t),t)

(53) u(t) is a white Gaussian noise with zero mean. Correspondingly, a process noise conversion matrix is:

(54) $\Gamma \left(t_i,t_{i-1}\right)=\frac{\partial X\left(t_i\right)}{\partial u\left(t_{i-1}\right)}$

(55) The EKF is generally used in real-time application scenarios. Given an initial state and time of the satellite, an EKF processes measurement data based on continuous measurement epochs, and a to-be-estimated state is updated.

(56) The process includes two parts. The first step is “time update”. A satellite orbit state and a state estimation covariance are advanced from the previous time t.sub.i−1 (add “+” to a right superscript) to the current observation time t.sub.i (add “−” to a right superscript):
X.sub.i.sup.−=X(t.sub.i;X(t.sub.i−1)=X.sub.i−1.sup.+)+Γ(t.sub.i,t.sub.i−1)u(t.sub.i−1)
P.sub.i.sup.−=Φ(t.sub.i,t.sub.i−)P.sub.i−1.sup.+Φ.sup.T(t.sub.i,t.sub.i−1)+Γ(t.sub.i,t.sub.i−1)Q.sub.cΓ.sup.T(t.sub.i,t.sub.i−1)

(57) where P is an orbit state covariance matrix, Q.sub.c is a process noise covariance matrix and is a (3J+3K+3L)-dimensional diagonal matrix, a diagonal value is related to dynamic model accuracy of each satellite, and a variance of an acceleration error is set as σ.sub.c.sup.2. A format of the process noise conversion matrix Γ is:

(58) $\Gamma =\left[\begin{array}{cccc}{\Gamma }_0& & & 0\\ & {\Gamma }_1& & \\ & & \ddots & \\ 0& & & {\Gamma }_{J+K+L‐1}\end{array}\right]$

(59) where a diagonal submatrix is denoted as:

(60) ${\Gamma }_k={\left[\begin{array}{c}\Delta t^2/2\times I_{3\times 3}\\ \Delta t\times I_{3\times 3}\end{array}\right]}_{6\times 3}$

(61) where Δt=t.sub.i−t.sub.i−1, I is an identity matrix, 0 is a zero matrix, and k=0, 1, . . . J+K+L−1.

(62) The process noise covariance matrix can appropriately slow down the convergence rate of the orbit state covariance to maintain the continuous correction of the satellite state by measured values at the end of filtering. Especially in the case of considering the error of the orbital dynamic model, the process noise covariance matrix can adjust the weight of the orbital dynamic model and the measured values to ensure the stable convergence of the navigation process. The state covariance matrix obtained at any measurement time is an estimation of the uncertainty of the satellite state accuracy at this time, which reflects the convergence and accuracy of the orbit state.

(63) The second step is “measurement update”. The solution of a satellite orbit state and update of an orbit state covariance are completed by using measurement data:
K.sub.i=P.sub.i.sup.−H.sub.i.sup.T(R.sub.i+H.sub.iP.sub.i.sup.−H.sub.i.sup.T).sup.−1
X.sub.i.sup.+=X.sub.i.sup.−+K.sub.i(z.sub.i−h(X.sub.i.sup.−))
P.sub.i.sup.+=(I−K.sub.iH.sub.i)P.sub.i.sup.−(I−K.sub.iH.sub.i).sup.T+K.sub.iR.sub.iK.sub.i.sup.T

(64) where K.sub.i is a Kalman gain, H.sub.i is a design matrix and is calculated by equation (12), I is an identity matrix, R.sub.i is a measuring noise covariance diagonal matrix, a diagonal element value corresponding to laser measurement data is σ.sub.DL.sup.2, and a diagonal element value corresponding to microwave measurement data is σ.sub.ML.sup.2 .

Embodiment 2

(65) In this embodiment, the steps of the method in Embodiment 1 are used to perform simulation scene experiment determination:

(66) The simulation time is set at 0:00 on Jan. 1, 2018 (UTC time), and the simulation time is 14 days. A DRO preferred in Embodiment 1 is selected, with a plane amplitude of 49000 km in the earth-moon rotating coordinate system. Table 1 gives an initial state of the DRO in an earth-centered inertial J2000 coordinate system. The LEO orbit selected is a nearly circular orbit with an altitude of 400 km and an inclination angle of 43°. The initial state is shown in Table 2.

(67) TABLE-US-00001 TABLE 1 Initial state of DROs in a J2000 coordinate system X [km] Y [km] Z [km] Vx [km/s] Vy [km/s] Vz [km/s] −20018.812 298708.811 126160.320 1.263656 0.198710 0.182677

(68) TABLE-US-00002 TABLE 2 Initial value state of orbits of LEO data relay satellites Right Semi- ascension Argument major Inclination of ascending of True axis Eccen- angle node perigee anomaly [km] tricity [deg] [deg] [deg] [deg] LEO 6778.140 0.002 43.0 0.0 0.0 0.0

(69) Taking a GPS navigation constellation as an example, according to a GPS almanac, 31 GPS satellites are selected and distributed in 6 orbit planes. The initial state is shown in Table 3.

(70) TABLE-US-00003 TABLE 3 Initial value state of orbits of GPS constellation satellites Right ascension Semi- of Argument major Inclination ascending of True axis Eccen- angle node perigee anomaly PRN [km] tricity [deg] [deg] [deg] [deg] 1 26559.562 0.005 55.3 92.0 28.3 237.0 2 26559.447 0.015 54.0 88.9 232.7 272.9 3 26558.923 0.001 55.0 153.5 197.6 28.9 4 26559.346 0.012 54.0 89.6 63.9 232.5 5 26559.819 0.004 54.3 152.5 26.0 42.5 6 26558.616 0.001 55.2 91.6 159.7 45.7 7 26559.960 0.009 55.4 272.8 206.1 187.3 8 26558.963 0.001 55.1 32.1 261.1 105.6 9 26560.176 0.001 54.8 211.5 135.4 314.5 11 26559.120 0.016 51.3 32.2 84.3 219.2 12 26559.110 0.005 56.7 337.5 36.7 189.4 13 26559.960 0.005 55.7 218.7 117.7 237.6 14 26559.829 0.008 55.3 216.4 248.3 351.3 15 26560.443 0.008 53.4 208.0 24.2 290.2 16 26559.709 0.008 56.7 338.5 18.1 54.9 17 26558.903 0.011 55.9 35.5 245.2 24.9 18 26559.482 0.016 53.1 150.8 249.4 95.7 19 26575.076 0.011 55.7 37.9 37.6 119.9 20 26559.663 0.005 53.1 147.7 72.0 319.7 21 26560.287 0.023 53.6 89.5 253.0 123.6 22 26559.855 0.008 52.9 150.8 242.9 54.6 23 26560.921 0.010 54.2 211.5 209.7 248.5 24 26559.497 0.004 54.4 269.6 16.8 258.5 25 26559.316 0.005 56.0 334.2 41.5 146.5 26 26559.729 0.000 55.0 333.2 349.8 109.5 27 26560.136 0.003 55.6 32.6 17.2 341.1 28 26560.267 0.020 56.7 338.9 265.1 26.1 29 26559.814 0.001 55.9 35.9 331.1 137.0 30 26559.824 0.002 54.6 275.1 181.1 187.2 31 26559.230 0.008 55.7 273.6 330.7 211.3 32 26559.839 0.011 54.3 157.8 9.1 215.1

(71) Inter-satellite measurement is performed between an LEO data relay satellite (the number of the satellite is 1, L=1) and a DRO satellites (the number of the satellite is 1, J=1), and a standard deviation of a measurement error is set to 0.01 m (16). Inter-satellite measurement is performed between an LEO satellite and each navigation constellation satellite (the number of the satellites is 31, K=31), and a standard deviation of a measurement error is also set to 0.01 m (1σ). There is no measurement link between navigation satellites. The measurement accuracy setting of 0.01 m is basically in line with the actual space laser and microwave measurement accuracy. Planetary occlusion is taken into account in the measurement. At the same time, it is assumed that an antenna of an LEO satellite navigation receiver points to the zenith direction, and measurement signals can be received when the antenna elevation angle is 5° or more. For the sake of simplicity, all satellites adopt the same orbital dynamic model setting. The gravity fields of the earth and the moon both adopt a 70-order spherical harmonic model. Taking the solar particle gravity, the solar radiation light pressure and a near-earth atmospheric drag model into account, the area-to-mass ratio is set to a fixed value of 0.002 m.sup.2/kg, and the light pressure coefficient C.sub.R and drag coefficient C.sub.D are 1.3 and 2.5 respectively. See Table 4 for specific simulation settings.

(72) TABLE-US-00004 TABLE 4 Dynamic model and measurement parameter settings Ephemeris model JPL DE405 Planetary ephemeris JPL DE405 Lunar libration angle Reference coordinate J2000 frame Mechanical model Substantially Gravity field of the earth 70 × 70 GGM02C Gravity field of the moon 70 × 70 GL0660B Three-body gravity Solar particle gravity Solar radiation light Area-to-mass ratio: 0.002 m.sup.2/kg, pressure C.sub.R:1.3 Near-earth atmospheric Area-to-mass ratio: 0.002 m.sup.2/kg, drag C.sub.D:2.5 Numerical integrator Runge-Kutta-4, with a step size of 60 s Measurement model Instantaneous position White Gaussian noise measurement Measurement period Average: 0 m Standard deviation: 0.01 m (1σ) Δt = 60 s

(73) The dynamic model in the navigation calculation process also adopts the settings in Table 4. An initial position error of each satellite is 1000 m (1σ) and an initial velocity error is 0.1 m/s. Diagonal elements in a noise matrix during extended Kalman filter are all set to 0.01.sup.2 m.sup.2, which is consistent with a measurement error. A diagonal element σ.sub.c in the process noise matrix Q.sup.c of equation (15) corresponding to a DRO satellite state is set to 1×10.sup.−9 m/s.sup.2, and σ.sub.c corresponding to navigation satellite and LEO satellite states is set to 1×10.sup.−8 m/s.sup.2.

(74) The position uncertainty is calculated from diagonal elements of the state covariance matrix in the filter process, which can characterize the convergence and accuracy of the position filtering process. FIG. 7 shows changes of position filtering uncertainty of GPS satellites (taking PRN 1 as an example), DRO satellites and LEO satellites with time. Only simulation results of GPS satellites with PRN of 1 are given, and other GPS satellites have similar results. It can be seen from the figure that the LEO satellite converges fastest, the GPS satellite comes second, the DRO satellite is slowest, but in the end all satellite position states converge to higher accuracy. The GPS satellite converges to 0.1 m, the DRO satellite converges to 1 m, and the LEO satellite converges to cm level. FIG. 8 shows changes of the position residual (namely the difference between the position filtering solution and the true position value) with time. It can be seen that the convergence rate of position filtering for all satellites is very fast, and the final positioning accuracy is very high. When the filtering process converges and then the positioning accuracy is counted, the 3D RMS of DRO position residual is 1 m, the 3D RMS of the LEO is 0.01 m, which reaches the inter-satellite measurement accuracy, and the 3D RMS of the GPS satellite (PRN 1) is 0.04 m, which also achieves very high accuracy. FIG. 9 shows the positioning accuracy of all satellites of the GPS constellation when the inter-satellite measurement accuracy is 0.01 m, which shows that the positioning accuracy of all satellites is better than 0.06 m. This simulation result is approximately equivalent to the accuracy of the high-precision ephemeris of the GPS constellation.

(75) When other simulation conditions remain unchanged and the inter-satellite measurement accuracy is relaxed to 0.1 m, as shown in FIG. 10, the positioning accuracy of the DRO satellite is 1.6 m, and the positioning accuracy of the LEO satellite is 0.03 m. The positioning accuracy of most satellites of the GPS constellation can reach 0.15 m.

(76) Two sets of simulation results are given above. When the inter-satellite laser and microwave measurement accuracy is set to 0.01 m, the autonomous navigation position accuracy of GPS constellation satellites reaches 0.06 m. When the measurement accuracy is set to 0.1 m, the GPS autonomous navigation position accuracy reaches 0.15 m. The simulation results show that if the inter-satellite measurement accuracy is better than 0.1 m, the technical solution of autonomous navigation of the navigation constellation satellites proposed in this patent is feasible. Even if there is only one DRO satellite and one LEO satellite, the navigation constellation can realize high-precision autonomous navigation without depending on ground measurement and control.

(77) The foregoing descriptions are only preferred implementations of the present invention. It should be noted that for a person of ordinary skill in the art, several improvements and modifications may further be made without departing from the principle of the present invention. These improvements and modifications should also be deemed as falling within the protection scope of the present invention.