CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority under 35 U.S.C. § 119(b) to Chinese Application No. 201910263156.1, filed Apr. 2, 2019, the disclosure of which is incorporated by reference herein in its entirety.

TECHNICAL FIELD

The present invention relates to the field of satellite navigation technologies, and in particular, to a method for achieving space-based autonomous navigation of global navigation satellite system (GNSS) satellites.

BACKGROUD OF THE INVENTION

At present, a GNSS includes a GPS navigation constellation, a Beidou navigation constellation, a Glonass navigation constellation, and the like, and provided positioning and timing services have been widely used in various fields such as agriculture, industry and commerce, national defense, and people's livelihood. The satellite navigation system generally includes three parts: a space segment, a ground segment and a user segment. A main body of the space segment is a certain number of satellites distributed in a space orbit, usually located in a medium earth orbit (MEO), a geostationary orbit (GEO) or an inclined geosynchronous orbit (IGSO). A main function of the navigation satellites is to continuously transmit navigation signals to the near-earth space and the ground, so that at least four signals can be received at any position at any time. The ground segment is responsible for the overall management of the whole system. Its main functions include: tracking the whole constellation satellite, calculating a position state of the navigation satellite through measurement signals of the satellite and a ground station, and calculating a clock difference of a satellite clock to ensure the time synchronization; updating and uploading navigation message data, and the like. The user segment is to compute the position, velocity and time (PVT) of a user after receiving navigation signals by a GNSS receiver. The concept of autonomous navigation of navigation satellites mentioned in this patent mainly refers to the position determination (positioning) realized by navigation constellations independent of the ground station's measurement and control. And the timing of the navigations is not considered in this patent.

GNSS satellites has the requirement of autonomous navigation, moreover, there exist several problems:

1) Delay effects of a troposphere and an ionosphere on navigation signals

Navigation signals are subjected to delay effects in the atmospheric troposphere and ionosphere during the propagation to the ground. Although the delay amount is also calculated as a data product in the ground analysis process, the decrease in the measurement signal accuracy between satellites and the ground affects the positioning and timing accuracy of the navigation satellites to a certain extent.

2) The measurement geometry between satellites and the ground changes slowly.

An MEO satellite among the navigation satellites has an orbital period of 12 hours, and a measurement geometry (namely a relative geometric position relationship) between the MEO satellite and a ground station changes relatively slowly. An IGSO satellite among Beidou satellites has an orbital period of 24 hours, and even the Beidou GEO satellite is stationary with respect to the ground station. The slow change of the conventional satellite-ground measurement geometry is not conducive to the calculation of the position and clock difference of the navigation satellite.

BRIEF SUMMARY OF THE INVENTION

The present invention aims to provide a method for achieving autonomous navigation of GNSS satellites by using a DRO (Distant retrograde Orbit) and inter-satellite measurement.

In order to achieve the foregoing objective, the technical solution adopted by the present invention is as follows:

a method for achieving autonomous navigation of GNSS satellites, including the following steps:

S1: optimizing a cislunar DRO, and establishing a orbital dynamic model the cislunar space;

S2: establishing measurement links, by a low earth orbit (LEO) data relay satellite, with a DRO satellite and the GNSS satellites 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 autonomously determining the position and velocity of the GNSS satellites without using earth-based measurement.

For a near-earth two-body problem, only the size, shape , and the relative orientation of the orbits can be determined through inter-satellite links, and an absolute orientation of the orbit planes cannot be determined. If orbits of two satellites make the same rotation around the center of the earth, although an absolute orbit state changes, a relative series does not change. Therefore, in general, it is not observable for determining the absolute orbit of the satellite by using inter-satellite relative measurement in a symmetric gravitational field. However, there is a region in an earth-moon three-body orbital space, whose magnitudes influenced by the gravitational effect of the earth and the moon can be comparable, and the two are no longer different in order of magnitude. The distribution of the gravity field in this region is asymmetric. If the orbit of a satellite moves in this region, the asymmetry of the gravity field can provide absolute pointing information of the satellite orbit plane, and then the absolute orbit plane of the satellite can also be determined. This technology is called linked autonomous inter-planetarium satellite navigation (LiAISON), which makes it possible to determine the absolute orbit of a satellite by using inter-satellite relative measurement. Previous studies by scholars have mainly focused on autonomous navigation using aa halo orbit in cislunar space. However, the main problem of the halo orbit is unstable, and orbit maneuver is often required to maintain the orbit, so the halo orbit does not have advantageous for providing long-term navigation information for GNSS satellites. However, there is still a kind of DROs in the cislunar space, which is stable in the long-term even considering all the orbit perturbation errors. The first step of the present invention is to select one preferred DRO.

Preferably, in step S 1, the step of optimizing an earth-moon space DRO specifically includes:

A1: deducing an earth-moon DRO family with periodic characteristics from a circular restricted three-body problem (CRTBP) model.

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 asymmetry to evaluate the degree of asymmetry of the gravity field in different spatial positions.

A2: calculating gravity asymmetry of each orbit, and selecting an orbit with high gravity asymmetry as a preferred earth-moon DRO.

Preferably, in step S 1, the step of establishing a dynamic model of an earth-moon space satellite orbit specifically includes:

B 1: establishing satellite state vectors: at any time t, yDj(t), yMk(t) and yTl(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 yDj(t)32[rDj(t)T,vDj(t)T]T, yMj(t)=[rMj(t)T,vMj(t)T]Tand yTj(t)=[rTj(t)T,vTj(t)T]T;

where 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:
yD(t)=[yD0(t)T, . . . ,yDJ−1(t)T)]T

and
yM(t)=[yM0(t)T, . . . ,yMK−1(t)T)]T

and
yT(t)=[yT0(t)T, . . . ,yTL−1(t)T)]T

which are further combined into a state vector:
X(t)=[yD(t)T,yM(t)T,yT(t)T]T

B2: extrapolating an orbit state of any satellite according to a first-order dynamic differential equation:

[r.v.]=[vaE+aM+aS+aSRP+aD]

where aEand aMm are spherical harmonic gravity acceleration models of the earth and the moon respectively, aSis a gravity acceleration of the solar particle, aSRPis an acceleration of solar radiation light pressure, and aDis mainly an atmospheric drag acceleration suffered by an LEO satellite; and

B3: calculating a state transition matrix Φi□Φ(ti,ti−1)=∂X(ti)/∂X(ti−1) from any time ti−1to the next time ti.

Preferably, in step B2, the first-order dynamic differential equation can be described as:
{dot over (X)}(t)=f(X(t),t)

in step B3, the state transition matrix n follows a variational equation:

Φ.i=∂f⁡(ti,X⁡(ti))∂X⁡(ti)⁢Φi.

Preferably, step S2 specifically includes:

S21: at any time ti, 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:

{zDT,ijl=hDTjl(yDj(ti),yTl(ti))+εDL=❘"\[LeftBracketingBar]"rDj(ti)-rTl(ti)❘"\[RightBracketingBar]"+εDLzMT,ikl=hMTkl(yMk(ti),yTl(ti))+εML=❘"\[LeftBracketingBar]"rMk(ti)-rTl(ti)❘"\[RightBracketingBar]"+εML

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

S22: Xi−[(yD,i−)T,(yM,i−)T,(yT,i−)T]Tbeing a known initial state of a satellite, and linearizing the measurement model (8) near the state:

zDT,ijl=hDTjl(yD,ij,-,yT,il,-)+HDD,ijl·Δ⁢yD,ij+HDT,ijl·Δ⁢yT,il+εDT⁢zMT,ikl=hMTkl(yM,ik,-,yT,ij,-)+HMM,ikl·Δ⁢yM,ik+HMT,ikl·Δ⁢yT,il+εMT⁢where⁢HDD,ijl⁢▯⁢∂hDT,ijl∂yD,ij,-,HDT,ijl⁢▯⁢∂hDT,ijl∂yT,il,-,HMM,ikl⁢▯⁢∂hMT,ikl∂yM,ik,-,HMT,ikl⁢▯⁢∂hMT,ikl∂yT,il,-,
and there are relations HDT,ijl=−HDD,ijland HMT,ikl=−HMM,ikl;

ΔXi=[(ΔyD,i−)T,(ΔyM,i−)T,(ΔyT,i−)T]Tis set to be an increment of a to-be-estimated state with respect to an initial value state, and then ΔXi=Xi−Xi−;

S23: combining all measured values into

zi=[zDT,iT,zMT,iT]T,
with a component denoted as:

zDT,i=[zDT,i00,…,zDT,iJ-1,0,…,zDT,i0,L-1,…,zDT,iJ-1,L-1]T⁢zMT,i=[zMT,i00,…,zMT,iK-1,0,…,zMT,i0,L-1,…,zMT,iK-1,L-1]T

where a corresponding model value is hi(Xi)=[hDT,iT,hMT,iT]T, and a measurement error is ε=[(εDL)T,(εML)T]; standard deviations of microwave measurement and laser measurement errors are σDLand σMLin respectively; and

S24: at the time ti, a linearization measurement model including all measured values is:

[zDT,izMT,i]=[hDT,i(Xi-)hMT,i(Xi-)]+[HDD,i0HDT,i0HMM,iHMT,i]·[Δ⁢yD,iΔ⁢yM,iΔ⁢yT,i]+[εDεM]

which is further simplified as:
zi=hi(Xi−)+Hi·ΔXi+ε

Further, there are two types of measurement links used in establishing the measurement model. As shown inFIG. 2, one type is laser measurement links of DRO satellites and LEO satellites, which need to be provided with laser transceivers respectively. The other type is microwave measurement links of navigation constellations and LEO satellites, which are LEO satellites. The data measured between the two types of satellites is estimated in orbit by a Kalman filter method, and orbit states of all satellites in the link are at the same time. ephemeris of the navigation constellation is injected into the navigation satellite to maintain the ability of the navigation satellite to provide high-accuracy navigation services.

Preferably, step S3 specifically includes:

S31: adding a process noise u(t) to a dynamic equation:
{dot over (X)}(t)=f(X(t),u(t),t)

where u(t) is a white Gaussian noise with zero mean value; correspondingly, a process noise transform matrix is:

Γ⁡(ti,ti-1)=∂X⁡(ti)∂u⁡(ti-1)

S32: given an initial state and time of the satellite, processing, by an EKF, measurement data based on continuous measurement epochs, and updating a to-be-estimated state.

Preferably, step S32 specifically includes:

A. time update

where a satellite orbit state and a state estimation covariance are advanced from the previous time ti−1(add “+” to a right superscript) to the current observation time ti, (add “−” to a right superscript):
Xi−=X(ti;X(ti−1)=Xi−1+)+Γ(ti,ti−1)u(ti−1)
Pi−=Φ(ti,ti−)Pi−1+ΦT(ti,ti−1)+Γ(ti,ti−1)QcΓT(ti,ti−1)

where P is an orbit state covariance matrix, Qcis 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 σc2; a format of the process noise conversion matrix Γ is:

Γ=[Γ00Γ1⋱0ΓJ+K+L-1]

where a diagonal submatrix is denoted as:

Γk=[Δ⁢t2/2×I3×3Δ⁢t×I3×3]6×3

where Δt=ti−ti−1, I is an identity matrix,0is 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:
Ki=Pi−HiT(Ri+HiPi−HiT)−1
Xi+=Xi−+Ki(zi−h(Xi−))
Pi+=(I−KiHi)Pi−(I−KiHi)T+KiRiKiT

where Kiis a Kalman gain, Hiis a design matrix and is calculated by equation (12), I is an identity matrix, Riis a measuring noise covariance diagonal matrix, a diagonal element value corresponding to laser measurement data is σDL2, and a diagonal element value corresponding to microwave measurement data is σML2.

It is worth noting that there is currently no effective solution for autonomous navigation of navigation constellation satellites. Some scholars have proposed using X-ray pulsar signals to realize autonomous navigation of navigation satellites, but the solution has the following shortcomings. First, a timing measurement error of the pulsar signals is relatively poor, which is far from meeting the demand of navigation satellites for high-accuracy positioning. At present, the SEXTANT project in the United States is carrying out pulsar navigation experiments on the International Space Station, with an in-orbit positioning accuracy of only 5 km. Second, the solution requires each navigation satellite to carry a heavy and expensive X-ray detector, which does not conform to the current operation mode of the navigation satellite. The links between navigation satellites currently being studied cannot realize autonomous navigation of satellites either. Due to the symmetry of the gravity field in the near-earth space, the rotation of the constellation cannot be sensed, and navigation satellite positioning still relies mainly on ground measurement and control stations. The main objective of the inter-satellite link is to keep the time synchronization between satellites. Therefore, there is no effective solution at present to meet both high-accuracy and long-term autonomous positioning demands.

There are two key points for realizing this technical solution. One key point is that some orbits in an earth-moon space retrograde orbit family are located in a space region with strong asymmetry of an earth-moon three-body gravity field; and if the orbits establish inter-satellite measurement links with earth satellites, it is expected to determine absolute positions of all satellites at both ends of the link in an inertial space through inter-satellite relative measurement. The other key point is that in view of the objective reality that a DRO cannot directly establish a high-accuracy measurement link with GNSS satellite orbit, a concept of a LEO data relay satellite is proposed innovatively, and indirect high-accuracy measurement of the DRO and the GNSS satellite orbit is achieved.

The present invention has the following beneficial effects:

This patent proposes an innovative method to implement an indirect link between a DRO and navigation constellation through a LEO data relay satellite. Suppose the LEO data relay satellite carries a GNSS receiver and a laser measuring device. A laser link is established between a DRO satellite and the LEO data relay satellite, and a microwave link is established between the navigation constellation and the LEO data relay satellite, which has the following advantages:

(1) The existing operation mode of the navigation constellation is not affected, and the LEO data relay satellite can receive high-accuracy microwave measurement data of the navigation satellite.

(2) The selected DRO is 330,000 to 440,000 km away from the earth, with a large observation field of view. At the same time, the DRO has long-term stability and is expected to greatly reduce orbital maintenance costs.

(3) The period of LEO satellite is short, the geometric range between the LEO satellite the GNSS navigation constellation and the DRO changes significantly, and the observation geometry is better, which is conducive to navigation calculation.

(4) Space laser and microwave links are not affected by the earth's atmosphere and ionosphere delay, and the measurement data processing accuracy is higher.

