A real-time high-precision PNT service method based on space-ground joint observation resources

Abstract

The present invention discloses a real-time high-precision PNT service method based on joint observation resources of the earth and the sky, and belongs to the field of navigation and positioning. The method estimates the precise orbit and clock error products of GNSS satellites based on a small number of regionally distributed GNSS ground monitoring stations and a small number of low-orbit satellite-borne GNSS observation data; carries out GNSS satellite orbit and clock error prediction processing, fits into GNSS satellite ephemeris products, annotates to visible low-orbit satellites, and distributes to the entire low-orbit constellation; low-orbit satellites solve precise orbits and clock errors in real time, calculate predicted orbits, calculate predicted clock errors, fit predicted orbits and clock error products into low-orbit navigation messages and broadcast them to users; users carry out GNSS / low-orbit satellite joint precise single-point positioning processing, and solve user terminal positions and clock error information in real time. The present invention can effectively reduce the dependence on the global deployment of GNSS ground monitoring stations, and can realize real-time high-precision PNT services on a global scale based on a small number of regionally distributed ground monitoring networks.

Description

Technical Field

[0001] The present invention relates to the field of navigation and positioning, and specifically refers to a real-time high-precision PNT service method based on space-ground joint observation resources, which can meet the real-time high-precision PNT service needs of users under the condition of ground monitoring stations distributed in a small number of areas. Background Art

[0002] The Global Navigation Satellite System (GNSS) is based on ground monitoring network observation resources and can provide positioning, navigation and timing (PNT) services to users around the world. Limited by the accuracy of GNSS satellite broadcast ephemeris and pseudo-range measurement, the accuracy of GNSS standard single-point positioning service is at the meter level, which cannot meet the user's demand for real-time high-precision PNT services.

[0003] With the development of GNSS augmentation technology, the Ground-Based Augmentation System (GBAS) and the Satellite-Based Augmentation System (SBAS) estimate GNSS satellite orbits, clock errors, ionosphere and other correction information based on observation data from a certain number of ground monitoring stations, and broadcast them to the user end through the network or satellite, providing users in the service area with real-time high-precision PNT services at the decimeter to centimeter level. The GBAS system requires the deployment of a large number of ground monitoring stations within the service coverage area, making it difficult to achieve coverage of remote areas, deserts and ocean areas. The SBAS system, based on GEO satellites, can only broadcast wide-area augmentation information to specific areas, and still cannot meet the real-time high-precision PNT service needs of users in global regions.

[0004] With the construction and development of large-scale low-orbit Internet satellite constellations, researchers from commercial aerospace companies and scientific research institutes have proposed a satellite-based navigation enhancement technology based on low-orbit satellites, which uses a technical process of satellite-to-ground data transmission, ground information processing, and prediction of orbital clock error products to provide navigation enhancement services. Simulation studies have shown that low-orbit constellations that provide navigation enhancements must meet multiple coverage requirements to provide effective navigation enhancement services. If only the system operation mode of ground information processing is adopted, strict requirements will be placed on the performance of the satellite-to-ground link, the number of gateways, and information processing delay. The traditional large-capacity, long-chain processing mode will greatly increase the construction cost of the system and is not conducive to the long-term stable operation of the system. Summary of the invention

[0005] Aiming at the limitations of current GNSS augmentation services such as GBAS and SBAS and low-orbit navigation augmentation systems based on ground information processing, the present invention proposes a high-precision PNT service method based on space-ground joint observation resources. The method adopts a space-ground joint information processing mode and can provide real-time high-precision PNT services to users around the world.

[0006] The present invention is achieved through the following technical solutions:

[0007] A real-time high-precision PNT service method based on space-ground joint observation resources comprises the following steps:

[0008] (1) Based on the data of a small number of regionally distributed GNSS ground monitoring stations within 10 and a small number of low-orbit satellite-borne GNSS observation data within 10, the precise orbit and clock products of GNSS satellites are estimated by using the joint orbit determination and time synchronization processing method with loose constraints on the coordinates of the ground monitoring stations;

[0009] (2) Using the precise orbit and clock products of GNSS satellites, carry out orbit and clock error prediction processing, fit GNSS satellite ephemeris products based on the predicted orbit and clock error sequences of GNSS satellites, annotate GNSS satellite ephemeris products to visible low-orbit satellites through the satellite-to-ground link of the signal gateway station, and distribute them to the entire low-orbit constellation based on the inter-satellite link of low-orbit satellites;

[0010] (3) Based on GNSS satellite ephemeris products, onboard GNSS observations and LEO satellite inter-satellite link observation data, LEO satellites calculate precise orbits and clock errors in real time, use the dynamic orbit prediction method constrained by orbital altitude parameters to calculate the predicted orbit, and use the clock error modeling method that takes into account periodic terms to calculate the predicted clock error. The predicted orbit and clock error products of the LEO satellite are fitted into LEO navigation messages and then broadcast to users;

[0011] (4) Based on the downlink navigation signals and ephemeris products from GNSS and low-orbit satellites, users conduct joint precise single-point positioning processing of GNSS and low-orbit satellites to calculate the user terminal position and clock error information in real time.

[0012] Furthermore, the specific method of step (1) is:

[0013] (101) The observation equation is constructed based on the satellite and ground GNSS observation data:

[0014]

[0015]

[0016] In the formula, ρ LEOi and ρ GROUNDi Respectively represent the satellite and ground GNSS at time t i Observations, including pseudorange and phase observations; is the geometric distance of the low-orbit vehicle relative to the GNSS satellite, Where (x L ,y L ,z L ) and (x s ,y s ,z s) are the three-dimensional position coordinates of the antenna phase centers of the low-orbit vehicle and the GNSS satellite respectively; is the geometric distance between the ground station and the GNSS satellite, Where (x g ,y g ,z g ) and (x s ,y s ,z s ) are the three-dimensional position coordinates of the antenna phase center of the ground station and the GNSS satellite respectively; c is the speed of light in vacuum; δt r and δt s are receiver clock correction and GNSS satellite clock correction respectively; δ ion is the ionospheric delay correction; δ rel is the relativistic delay correction; δ trop is the tropospheric correction; δ tide is the tidal correction; δ wind-up Correction for phase wrapping; Correction of the phase center of ground station or low-orbit aircraft receiver; is the GNSS satellite phase center correction; LEOi and GNSSi is the observation noise of the corresponding observation;

[0017] (102) The observation equation is linearized and expressed as:

[0018] ρ LEOi =H LG (X GNSSi ,X LEOi ,X LGi ,t i )+ξ LEOi

[0019] ρ GROUNDi =H GG (X GNSSi ,X GGi ,t i )+ξ GNSSi

[0020] Where, X LEOi and X GNSSi They represent the low-orbit satellite and GNSS satellite at time t i The state vector, X LGi and X GGi In addition to X LEOi and X GNSSi Other parameters to be estimated besides ;

[0021] Linearize the nonlinear observation equation:

[0022]

[0023]

[0024] In the formula, Respectively represent the correction number of the corresponding parameters, l LEOi and l GROUNDi It represents the difference between the observed value and the theoretical value;

[0025] The GNSS and LEO satellite correction numbers solved by the above linearized equations are used to obtain the precise orbit and clock products of the GNSS and LEO satellites through orbital integration.

[0026] Furthermore, the specific method of step (2) is:

[0027] (201) Using the GNSS satellite precise orbit sequence estimated in step (1) as a virtual observation value, the initial position vector is estimated using the dynamic parameters. The estimation process is as follows:

[0028] Assume the initial position and velocity of the satellite are Use orbital integration to get the position and velocity of the satellite at time t

[0029]

[0030] In the formula, is the state vector of the perturbation parameters of the satellite, including parameters such as solar pressure, from which the position vector of the satellite is written:

[0031]

[0032] Assume that the position of the satellite in the precise ephemeris corresponding to time t is (X o Y o Z o ), then the observation equation for dynamic parameter estimation is expressed as:

[0033]

[0034] The observation equation is estimated by least squares to obtain the initial position, velocity and perturbation parameter state vector correction of the satellite, thereby realizing the estimation of the dynamic parameters of the satellite orbit;

[0035] The predicted orbit product of the GNSS satellite is obtained by using the dynamic parameters and the initial state of the satellite to perform orbital integral calculations;

[0036] (202) Using the precise clock error product of the GNSS satellite estimated in step (1), the satellite clock error prediction is carried out based on the quadratic polynomial model. The process is as follows:

[0037] The satellite clock error model is expressed as

[0038] clk i =a 0 +a 1 (tt c )+a 2 (tt c ) 2

[0039] Where clk i is the satellite clock error at time t, t c is the reference time of the clock error model, a 0 t c The satellite clock error at time a 1 t c The clock speed parameter of the clock at this moment, a 2 t c The frequency drift parameter of the clock at this moment;

[0040] Based on the preprocessed GNSS satellite precise clock error sequence, the quadratic polynomial model parameter a of the precise clock error sequence is estimated. 0 , a 1 , a 2 ; Then the GNSS satellite clock error during the forecast period is calculated using the quadratic polynomial model parameters;

[0041] (203) The least squares method is used to fit the orbit of the GNSS satellite. The process is as follows:

[0042] Assume that the satellite ephemeris model parameter at the reference epoch is σ(t 0 ), then the observation equation can be expressed as

[0043]

[0044] In the formula, is the position component of the satellite; linearizing the above equation, we get

[0045] y=HΔσ

[0046] In the formula,

[0047] According to the least squares estimation principle, the optimal estimate of Δσ is:

[0048] Δσ=(H T H) -1 H T y

[0049] Using Δσ to correct σ(t 0) to obtain the satellite ephemeris model parameters of the GNSS satellite reference epoch; use the quadratic polynomial model to fit the GNSS predicted clock error into the clock error parameters of the reference epoch; compile the GNSS satellite ephemeris model parameters and the GNSS clock error parameters into a broadcast ephemeris, which is uploaded to the visible low-orbit satellites through the satellite-to-ground link of the signal gateway station, and distributed to the entire low-orbit constellation based on the low-orbit satellite inter-satellite link.

[0050] Furthermore, the specific method of step (3) is:

[0051] (301) Based on GNSS satellite ephemeris products, onboard GNSS observations and low-orbit inter-satellite link observation data, low-orbit satellites can calculate precise orbits and clock errors in real time. The observation equation is written as:

[0052] Y i =G(X i ,t i )+ε i

[0053] Among them, Y i is the i-th observation, G(X i ,t i ) is t i The function of the state quantity at the moment, that is, the calculation formula of the theoretical value of the observation quantity, ε i is the observation noise;

[0054] Linearize the above equation and write it as:

[0055] y i =H i x 0 +ε i

[0056] Among them, y i t i The difference between the observed value and the theoretical value of the observed quantity at a moment, H i is the coefficient matrix of the error observation equation, expressed as Where Φ(t i ,t 0 ) is the state transfer matrix;

[0057] According to the least squares principle, the estimated parameters to be solved are:

[0058]

[0059] The parameters to be estimated are is m 0 Global parameters, including satellite orbit state parameters, station coordinates, ambiguity parameters, atmosphere, and solar radiation pressure parameters; is the ith observation epoch m iclock error parameters to be estimated, i = 1, 2, ... N;

[0060] (302) The predicted orbit is calculated using a dynamic orbit prediction method constrained by orbital altitude parameters; the specific method is:

[0061] The gravitational potential of the Earth's non-spherical gravitational perturbation force is expressed as

[0062]

[0063] Among them, GM E is the earth's gravitational constant, r is the distance from the satellite to the earth's center, a e is the semi-major axis of the Earth, N is the order of the non-spherical gravitational perturbation of the Earth, and λ are the geocentric latitude and longitude, are all normalized quantities;

[0064] The order of the non-spherical gravitational perturbation of the Earth constrained by orbital altitude parameters is used to achieve efficient dynamical orbit prediction. The specific parameter constraints are set as follows:

[0065] H≤500km,N=90

[0066] 500km≤H≤1000km,N=60

[0067] H≥1000km,N=30

[0068] Where H is the orbital height of the low-orbit satellite, and N is the order of the non-spherical gravitational perturbation of the Earth;

[0069] (303) The forecast clock error is calculated by adopting the clock error modeling method taking into account the periodic term; the specific method is as follows:

[0070] After the spectrum analysis of the satellite clock error series is carried out and its obvious periodic terms are determined, the satellite clock error modeling taking into account the periodic terms is used to carry out clock error prediction. The function model is:

[0071]

[0072] Where t is the epoch time; A 0 , A i , B i is the parameter to be estimated, i=1,2,...,n; n is the order of the trigonometric series; ω i is the frequency component of the low-orbit satellite clock, obtained by spectrum analysis; ε is the random error term;

[0073] The clock error prediction sequence is fitted based on the corresponding clock error model. After the clock error parameters are obtained, they are compiled into a low-orbit satellite message together with the precise ephemeris product of the low-orbit satellite. They are then distributed to the entire low-orbit constellation based on the low-orbit satellite inter-satellite link and broadcast to users.

[0074] Furthermore, the specific method of step (4) is:

[0075] (401) The user-side GNSS and low-orbit satellite pseudorange and carrier phase observation equations are established as follows:

[0076]

[0077] In the formula, the subscript r represents the user terminal number, i and j represent the frequency identifiers of the GNSS and low-orbit satellite downlink navigation signals, respectively, the superscript GNSS and m represent the GNSS satellite system and its satellite number, respectively, and the superscript LEO and n represent the low-orbit satellite system and its satellite number, respectively; and They represent the pseudorange and carrier phase observation values ​​of the frequency of satellite i in the GNSS system m observed by the user end, and They represent the pseudorange and carrier phase observation values ​​of the frequency j of the n-th satellite in the low-orbit system observed by the user end respectively; and Represent the geometric distances between the user terminal and the GNSS satellite m and the low-orbit satellite n, respectively, and are expressed as and Where X GNSS,m , X LEO,n and X r are the position vectors of GNSS satellite, low-orbit satellite and user terminal respectively; dt r , dT GNSS,m and dT LEO,n Represent the clock errors of the user end, GNSS satellite and low-orbit satellite respectively, represents the tropospheric delay of GNSS satellite m, represents the ionospheric delay of GNSS satellite m at frequency i, represents the ionospheric delay of low-orbit satellite n at frequency j, represents the wavelength of the navigation signal of GNSS satellite m frequency i, represents the wavelength of the navigation signal of frequency j of low-orbit satellite n, and Represents the carrier phase ambiguity of the downlink navigation signal from GNSS and low-orbit satellites, Represent the measurement noise of pseudorange and carrier phase of GNSS and LEO satellite respectively;

[0078] (402) Performing the user-side position and clock error solution based on the joint precise point positioning observation equation of GNSS and low-orbit satellite:

[0079] Calculate the position X of GNSS and LEO satellites based on the broadcast ephemeris of GNSS and LEO satellites GNSS,m and XLEO,n , and the clock difference dT between GNSS and LEO satellites GNSS,m and dT LEO,n ;

[0080] Read the pseudo-range and carrier phase observation data of GNSS and low-orbit satellites collected by the user end, and perform gross error elimination and cycle slip detection;

[0081] Calculate the floating-point solution of carrier phase ambiguity of GNSS and LEO satellites based on least squares or Kalman filter estimation method;

[0082] The LAMBDA method is used to fix the ambiguity parameters of GNSS and low-orbit satellites, thereby calculating the position parameter X of the user end. r and clock difference dt r .

[0083] Compared with the prior art, the present invention has the following beneficial effects:

[0084] (1) Based on the observation resources of space-based low-orbit satellites and ground-based GNSS monitoring networks, the present invention provides users with broadcast ephemeris products under a unified time and space reference framework through a joint processing method. Compared with the traditional ground-based observation resource mode, the PNT service performance is greatly improved;

[0085] (2) Based on the predicted orbit products of the ground processing center and the onboard GNSS observation data, low-orbit satellites estimate high-precision real-time clock difference products based on the joint observation data of the ground and the sky, which can effectively ensure the contribution of low-orbit satellites to the improvement of PNT service performance;

[0086] (3) The space-ground joint information processing mode is adopted. The precise orbit and clock error products of GNSS satellites are generated by processing data from a small number of ground stations and a small number of low-orbit satellites on the ground. The precise products of low-orbit satellites are solved on the satellite. This method effectively alleviates the requirements on the performance of the satellite-ground link, the number of gateway stations and the information processing delay, greatly reduces the construction cost of large-scale constellation low-orbit navigation augmentation systems, and is conducive to the long-term stable operation of the system.

[0087] (4) Space-based observation resources can effectively reduce the dependence on the global deployment of GNSS ground monitoring stations, and can provide real-time, high-precision PNT services on a global scale based on a small number of regionally distributed ground monitoring networks. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 The present invention provides a flow chart of a real-time high-precision PNT service method based on space-ground joint observation resources in an embodiment of the present invention. DETAILED DESCRIPTION

[0089] In order to better illustrate the purpose and advantages of the present invention, the technical solution of the present invention is further described below in conjunction with the accompanying drawings.

[0090] A real-time high-precision PNT service method based on space-ground joint observation resources. The method is based on a small number of regionally distributed GNSS ground monitoring stations and a small number of low-orbit satellite-borne GNSS observation data. The method uses a joint orbit determination and time synchronization processing method with loose constraints on the coordinates of the ground monitoring station to estimate the precise orbit and clock error products of the GNSS satellite. The GNSS satellite orbit and clock error prediction processing are carried out using a dynamic model and a quadratic polynomial model respectively, and are fitted into GNSS satellite ephemeris products. The products are then annotated to visible low-orbit satellites using the satellite-to-ground link of a gateway station. Based on the inter-satellite link analysis of the low-orbit satellite, the GNSS satellite ephemeris products are Sent to the entire low-orbit constellation; the low-orbit satellite calculates the precise orbit and clock error in real time based on the GNSS satellite ephemeris products, onboard GNSS and low-orbit inter-satellite link observation data, calculates the predicted orbit using the dynamic orbit prediction method constrained by orbital altitude parameters, and calculates the predicted clock error using the clock error modeling method that takes into account periodic terms. The predicted orbit and clock error products are fitted into a low-orbit navigation message and then broadcast to users; users carry out GNSS / low-orbit satellite joint precise single-point positioning processing based on the downlink navigation signals and ephemeris products of the GNSS / low-orbit satellite, and calculate the user terminal position and clock error information in real time.

[0091] like Figure 1 As shown, the specific steps include:

[0092] (1) Based on a small number of regionally distributed GNSS ground monitoring station data and a small number of low-orbit satellite-borne GNSS observation data, a joint orbit determination and time synchronization processing method with loose constraints on the ground monitoring station coordinates is used to estimate the precise orbit and clock products of the GNSS satellite;

[0093] (2) Utilize the precise ephemeris and clock products of GNSS satellites to carry out orbit and clock error prediction processing. Based on the predicted orbit and clock error sequences of GNSS satellites, fit the GNSS satellite ephemeris products. Use the satellite-to-ground link of the gateway station to annotate the GNSS satellite ephemeris products to visible low-orbit satellites, and distribute them to the entire low-orbit constellation based on the inter-satellite links of low-orbit satellites.

[0094] (3) Based on GNSS satellite ephemeris products, onboard GNSS observations and low-orbit inter-satellite link observation data, low-orbit satellites calculate precise orbits and clock errors in real time, use the dynamic orbit prediction method constrained by orbital altitude parameters to calculate the predicted orbit, and use the clock error modeling method that takes into account periodic terms to calculate the predicted clock error. The predicted orbit and clock error products are fitted into low-orbit navigation messages and then broadcast to users;

[0095] (4) Based on the downlink navigation signals and ephemeris products of GNSS / LEO satellites, users carry out GNSS / LEO satellite joint precise point positioning (PPP) processing to calculate the user terminal position and clock error information in real time.

[0096] The specific method of step (1) is as follows:

[0097] (101) The observation equation is constructed based on the satellite and ground GNSS observation data:

[0098]

[0099]

[0100] In the formula, ρ LEOi and ρ GROUNDi Respectively represent the satellite and ground GNSS at time t i Observations, including pseudorange and phase observations; is the geometric distance of the low-orbit vehicle relative to the GNSS satellite, Where (x L ,y L ,z L ) and (x s ,y s ,z s ) are the three-dimensional position coordinates of the antenna phase centers of the low-orbit vehicle and the GNSS satellite respectively; is the geometric distance between the ground station and the GNSS satellite, Where (x g ,y g ,z g ) and (x s ,y s ,z s ) are the three-dimensional position coordinates of the antenna phase center of the ground station and the GNSS satellite respectively; c is the speed of light in vacuum; δt r and δt s are receiver clock correction and GNSS satellite clock correction respectively; δ ion is the ionospheric delay correction; δ rel is the relativistic delay correction; δ trop is the tropospheric correction; δ tide is the tidal correction; δ wind-up Correction for phase wrapping; Correction of the phase center of ground station or low-orbit aircraft receiver; is the GNSS satellite phase center correction; LEOi and GNSSi is the observation noise of the corresponding observation.

[0101] (102) The observation equation is linearized and expressed as:

[0102] ρ LEOi =H LG (X GNSSi ,X LEOi ,XLGi ,t i )+ξ LEOi

[0103] ρ GROUNDi =H GG (X GNSSi ,X GGi ,t i )+ξ GNSSi

[0104] Where, X LEOi and X GNSSi They represent the low-orbit satellite and GNSS satellite at time t i The state vector, X LGi and X GGi In addition to X LEOi and X GNSSi Other parameters to be estimated.

[0105] Linearize the nonlinear observation equation and get:

[0106]

[0107]

[0108] In the formula, Respectively represent the correction number of the corresponding parameters, l LEOi and l GROUNDi Represents the difference between the observed value and the theoretical value. Using the GNSS and LEO satellite corrections solved by the above linearized equations, the precise orbit and clock products of GNSS and LEO satellites are calculated through orbit integration.

[0109] The specific method of step (2) is as follows:

[0110] (201) Using the GNSS satellite precise orbit sequence estimated in step (1) as a virtual observation value, the initial position vector is estimated with dynamic parameters. The estimation process is as follows:

[0111] Assume the initial position and velocity of the satellite are Use orbital integration to get the position and velocity of the satellite at time t

[0112]

[0113] In the formula is the state vector of the perturbation parameters of the satellite, mainly including the solar pressure parameter, from which the position vector of the satellite can be written:

[0114]

[0115] Assume that the position of the satellite in the precise ephemeris corresponding to time t is (X o Y o Z o ), then the observation equation for dynamic parameter estimation can be expressed as:

[0116]

[0117] The observation equation is estimated by least squares to obtain the initial position, velocity and perturbation parameter state vector correction of the satellite, and the dynamic parameters of the satellite orbit can be estimated. The predicted orbit product of the GNSS satellite is obtained by orbit integration calculation using the dynamic parameters and the initial state of the satellite.

[0118] (202) Using the precise clock error product of the GNSS satellite estimated in step (1), the satellite clock error prediction is carried out based on the quadratic polynomial model:

[0119] The satellite clock error model can be expressed as

[0120] clk i =a 0 +a 1 (tt c )+a 2 (tt c ) 2

[0121] Where clk i is the satellite clock error at time t, t c is the reference time of the clock error model, a 0 t c The satellite clock error at time a 1 t c The clock speed (frequency deviation) parameter of the clock at this moment, a 2 t c The frequency drift parameter of the clock at this moment.

[0122] Based on the preprocessed GNSS satellite precise clock error sequence, the quadratic polynomial model parameter a of the precise clock error sequence is estimated. 0 , a 1 , a 2 ; Then the GNSS satellite clock error during the forecast period is calculated using the quadratic polynomial model parameters.

[0123] (203) The orbit fitting of GNSS satellites adopts the least square method, and its process is as follows:

[0124] Assume that the satellite ephemeris model parameter at the reference epoch is σ(t 0 ), then the observation equation can be expressed as

[0125]

[0126] In the formula, is the satellite position component. After linearizing the above equation, we get

[0127] y=HΔσ

[0128] In the formula, Δσ=σ-σ'. According to the least squares estimation principle, the optimal estimate of Δσ is:

[0129] Δσ=(H T H) -1 H T y

[0130] Using Δσ to correct σ(t 0 ) can obtain the satellite ephemeris model parameters of the GNSS satellite reference epoch. The predicted clock error of GNSS is fitted into the clock error parameters of the reference epoch using a quadratic polynomial model. The ephemeris model parameters of the GNSS satellite and the clock error parameters of GNSS are compiled into a broadcast ephemeris in a certain format, which is then uploaded to the visible low-orbit satellites through the satellite-to-ground link of the gateway station and distributed to the entire low-orbit constellation based on the inter-satellite link of the low-orbit satellite.

[0131] The specific method of step (3) is as follows:

[0132] (301) Based on GNSS satellite ephemeris products, onboard GNSS observations and low-orbit inter-satellite link observation data, low-orbit satellites can calculate precise orbits and clock errors in real time. The observation equation can be written as:

[0133] Y i =G(X i ,t i )+ε i

[0134] where Y i is the i-th observation, G(X i ,t i ) is t i The function of the state quantity at the moment, that is, the calculation formula of the theoretical value of the observation quantity, ε i is the observation noise.

[0135] After linearization, the above equation can be written as:

[0136] y i =H i x 0 +ε i

[0137] in:

[0138] According to the least squares principle, the estimated parameters to be solved are:

[0139]

[0140] Among them: The parameters to be estimated are are m0 global parameters, including satellite orbit state parameters, station coordinates, ambiguity parameters, atmosphere, solar radiation pressure parameters, etc.; is the i-th (i=1,2,…N) observation epoch m i The clock error parameters to be estimated.

[0141] (302) The predicted orbit is calculated using the dynamic orbit prediction method constrained by orbital altitude parameters. The non-spherical gravitational perturbation force of the Earth is the largest error term affecting the accuracy of orbit prediction. The gravitational potential of this perturbation force can be expressed as

[0142]

[0143] Among them, GM E is the earth's gravitational constant, r is the distance from the satellite to the earth's center, a e is the semi-major axis of the Earth, N is the order of the non-spherical gravitational perturbation of the Earth, and λ are the geocentric latitude and longitude, All are normalized quantities. This method uses the Earth's non-spherical gravitational perturbation order constrained by orbital height parameters to achieve efficient dynamical orbit prediction. The specific parameter constraints are set as follows:

[0144] H≤500km,N=90

[0145] 500km≤H≤1000km,N=60

[0146] H≥1000km,N=30

[0147] (303) The clock error modeling method taking into account the periodic term is used to calculate the predicted clock error. Since the performance of the onboard atomic clock of the low-orbit satellite is not as good as that of the GNSS satellite clock, this method first performs a spectrum analysis on the satellite clock error sequence. After determining its obvious periodic term, the satellite clock error modeling taking into account the periodic term is used to carry out clock error prediction. The function model is:

[0148]

[0149] Where t is the epoch time; A 0 , A i , B i is the parameter to be estimated, where i = 1, 2, ..., n; n is the order of the trigonometric series, ω i is the frequency component of the low-orbit satellite clock, obtained by spectrum analysis; ε is the random error term.

[0150] The clock error prediction sequence is fitted based on the corresponding clock error model. After the clock error parameters are obtained, they are compiled into a low-orbit satellite message together with the precise ephemeris product of the low-orbit satellite. Then, they are distributed to the entire low-orbit constellation based on the low-orbit satellite inter-satellite link and broadcast to users.

[0151] The specific method of step (4) is as follows:

[0152] (401) The user-side GNSS and low-orbit satellite pseudorange and carrier phase observation equations are as follows:

[0153]

[0154] Wherein, the subscript r represents the user terminal number, i and j represent the frequency identifiers of the GNSS and low-orbit satellite downlink navigation signals, respectively, the superscript GNSS and m represent the GNSS satellite system and its satellite number, respectively, and the superscript LEO and n represent the low-orbit satellite system and its satellite number, respectively. and They represent the pseudorange and carrier phase observation values ​​of the frequency of satellite i in the GNSS system m observed by the user end, and They respectively represent the pseudorange and carrier phase observation values ​​of the j-frequency of the n-satellite in the low-orbit system observed by the user end. and Represent the geometric distances between the user terminal and the GNSS satellite m and the low-orbit satellite n, respectively, which can be expressed as and Where X GNSS,m , X LEO,n and X r They are the position vectors of GNSS satellite, low-orbit satellite and user terminal respectively. r , dT GNSS,m and dT LEO,n Represent the clock errors of the user end, GNSS satellite and low-orbit satellite respectively, represents the tropospheric delay of GNSS satellite m, represents the ionospheric delay of GNSS satellite m at frequency i, represents the ionospheric delay of low-orbit satellite n at frequency j, represents the wavelength of the navigation signal of GNSS satellite m frequency i, represents the wavelength of the navigation signal of frequency j of low-orbit satellite n, and Represents the carrier phase ambiguity of the downlink navigation signal from GNSS and low-orbit satellites, represent the measurement noise of pseudorange and carrier phase of GNSS and LEO satellite respectively.

[0155] (402) The method for solving the user-side position and clock error based on the joint precise point positioning observation equation of GNSS and low-orbit satellite is as follows:

[0156] Based on the broadcast ephemeris of GNSS and LEO satellites, the known quantities that can be calculated include the positions X of GNSS and LEO satellites. GNSS,m and X LEO,n , the clock difference between GNSS and LEO satellite dT GNSS,m and dT LEO,n .

[0157] Therefore, the parameters to be estimated in the above equation include the user terminal position parameter X r 、User-side clock error dt r , GNSS satellite tropospheric delay parameter T GNSS,m And the carrier phase integer ambiguity of GNSS and low-orbit satellites and Since the GNSS and low-orbit satellite dual-frequency joint precise single-point positioning algorithm is adopted, it is also necessary to additionally estimate the inter-system bias ISB and inter-frequency bias IFB parameters between GNSS and low-orbit satellites.

[0158] First, the pseudorange and carrier phase observation data of GNSS and low-orbit satellites collected by the user end are read to perform gross error elimination and cycle slip detection. Then, based on the least squares or Kalman filter estimation method, the floating-point solution of the carrier phase ambiguity of GNSS and low-orbit satellites is calculated. Finally, the ambiguity parameters of GNSS and low-orbit satellites are fixed using methods such as LAMBDA, so as to calculate the position parameter X of the user end. r and clock difference dt r .

[0159] In summary, the present invention adopts a space-ground joint information processing mode. First, based on the onboard GNSS observation data of a small number of regionally distributed ground monitoring stations and a small number of low-orbit satellites, the precise orbit and clock error products of the GNSS satellite are generated, and then injected to the low-orbit constellation through the satellite-to-ground link, and distributed to the entire low-orbit constellation by the low-orbit inter-satellite link. Based on the precise orbit and clock error products of the GNSS satellite, the onboard GNSS and the low-orbit inter-satellite link observation data, each low-orbit satellite solves the precise orbit and clock error in real time, and fits them into a low-orbit navigation message, which is broadcast to the ground together with the downlink navigation signal of the low-orbit satellite, providing real-time and high-precision PNT services to users around the world.

Claims

1. A real-time high-precision PNT service method based on joint space-ground observation resources, It is characterized in that The following steps are involved: (1) Based on the data of a small number of regionally distributed GNSS ground monitoring stations within 10 and a small number of low-orbit satellite-borne GNSS observation data within 10, the precise orbit and clock products of GNSS satellites are estimated by using the joint orbit determination and time synchronization processing method with loose constraints on the coordinates of the ground monitoring stations; (2) Using the precise orbit and clock products of GNSS satellites, carry out orbit and clock error prediction processing, fit GNSS satellite ephemeris products based on the predicted orbit and clock error sequences of GNSS satellites, annotate GNSS satellite ephemeris products to visible low-orbit satellites through the satellite-to-ground link of the signal gateway station, and distribute them to the entire low-orbit constellation based on the inter-satellite link of low-orbit satellites; (3) Based on GNSS satellite ephemeris products, onboard GNSS observations and inter-satellite link observation data of low-orbit satellites, the low-orbit satellites solve the precise orbit and clock error in real time, use the dynamic orbit prediction method constrained by orbital altitude parameters to calculate the predicted orbit, and use the clock error modeling method taking into account the periodic term to calculate the predicted clock error. The predicted orbit and clock error products of the low-orbit satellite are fitted into low-orbit navigation messages and then broadcast to users. The specific methods are as follows: (301) Based on GNSS satellite ephemeris products, onboard GNSS observations and low-orbit inter-satellite link observation data, low-orbit satellites can calculate precise orbits and clock errors in real time. The observation equation is written as: Y i =G(X i ,t i )+ε i Among them, Y i is the i-th observation, G(X i ,t i ) is t i The function of the state quantity at the moment, that is, the calculation formula of the theoretical value of the observation quantity, ε i is the observation noise; Linearize the above equation and write it as: y i =H i x 0 +ε i Among them, y i t i The difference between the observed value and the theoretical value of the observed quantity at a moment, H i is the coefficient matrix of the error observation equation, expressed as Where Φ(t i ,t 0 ) is the state transfer matrix; According to the least squares principle, the estimated parameters to be solved are: The parameters to be estimated are is m 0 Global parameters, including satellite orbit state parameters, station coordinates, ambiguity parameters, atmosphere, and solar radiation pressure parameters; is the ith observation epoch m i clock error parameters to be estimated, i = 1, 2, ... N; (302) The predicted orbit is calculated using a dynamic orbit prediction method constrained by orbital altitude parameters; the specific method is: The gravitational potential of the Earth's non-spherical gravitational perturbation force is expressed as Among them, GM E is the earth's gravitational constant, r is the distance from the satellite to the earth's center, a e is the semi-major axis of the Earth, N is the order of the non-spherical gravitational perturbation of the Earth, and λ are the geocentric latitude and longitude, are all normalized quantities; The order of the non-spherical gravitational perturbation of the Earth constrained by orbital altitude parameters is used to achieve efficient dynamical orbit prediction. The specific parameter constraints are set as follows: H≤500km,N=90 500km≤H≤1000km,N=60 H≥1000km,N=30 Where H is the orbital height of the low-orbit satellite, and N is the order of the non-spherical gravitational perturbation of the Earth; (303) The forecast clock error is calculated by adopting the clock error modeling method taking into account the periodic term; the specific method is as follows: After the spectrum analysis of the satellite clock error series is carried out and its obvious periodic terms are determined, the satellite clock error modeling taking into account the periodic terms is used to carry out clock error prediction. The function model is: Where t is the epoch time; A 0 , A i , B i is the parameter to be estimated, i=1,2,...,n; n is the order of the trigonometric series; ω i is the frequency component of the low-orbit satellite clock, obtained by spectrum analysis; ε is the random error term; Based on the corresponding clock error model, the clock error prediction sequence is fitted. After the clock error parameters are obtained, they are compiled into a low-orbit satellite message together with the precise ephemeris products of the low-orbit satellite. Then, they are distributed to the entire low-orbit constellation based on the low-orbit satellite inter-satellite link and broadcast to users. (4) Based on the downlink navigation signals and ephemeris products from GNSS and low-orbit satellites, users conduct joint precise single-point positioning processing of GNSS and low-orbit satellites to calculate the user terminal position and clock error information in real time.

2. A real-time high-precision PNT service method based on space-ground joint observation resources according to claim 1, It is characterized in that The specific method of step (1) is: (101) The observation equation is constructed based on the satellite and ground GNSS observation data: In the formula, ρ LEOi and ρ GROUNDi Respectively represent the satellite and ground GNSS at time t i Observations, including pseudorange and phase observations; is the geometric distance of the low-orbit vehicle relative to the GNSS satellite, Where (x L ,y L ,z L ) and (x s ,y s ,z s ) are the three-dimensional position coordinates of the antenna phase centers of the low-orbit vehicle and the GNSS satellite, respectively; is the geometric distance between the ground station and the GNSS satellite, Where (x g ,y g ,z g ) and (x s ,y s ,z s ) are the three-dimensional position coordinates of the antenna phase center of the ground station and the GNSS satellite respectively; c is the speed of light in vacuum; δt r and δt s are receiver clock correction and GNSS satellite clock correction respectively; δ ion is the ionospheric delay correction; δ rel is the relativistic delay correction; δ trop is the tropospheric correction; δ tide is the tidal correction; δ wind-up Correction for phase wrapping; Correction of the phase center of ground station or low-orbit aircraft receiver; is the GNSS satellite phase center correction; LEOi and GNSSi is the observation noise of the corresponding observation; (102) The observation equation is linearized and expressed as: r LEOi =H LG (X GNSSi ,X LEOi ,X LGi ,t i )+ξ LEOi r GROUNDi =H GG (X GNSSi ,X GGi ,t i )+ξ GNSSi In the formula, X LEOi and X GNSSi They represent the low-orbit satellite and GNSS satellite at time t i The state vector, X LGi and X GGi In addition to X LEOi and X GNSSi Other parameters to be estimated besides ; Linearize the nonlinear observation equation: In the formula, Respectively represent the correction number of the corresponding parameters, l LEOi and l GROUNDi It represents the difference between the observed value and the theoretical value; The GNSS and LEO satellite correction numbers solved by the above linearized equations are used to obtain the precise orbit and clock products of the GNSS and LEO satellites through orbital integration.

3. A real-time high-precision PNT service method based on space-ground joint observation resources according to claim 1, It is characterized in that The specific method of step (2) is: (201) Using the GNSS satellite precise orbit sequence estimated in step (1) as a virtual observation value, the initial position vector is estimated using the dynamic parameters. The estimation process is as follows: Assume the initial position and velocity of the satellite are Use orbital integration to get the position and velocity of the satellite at time t In the formula, is the state vector of the perturbation parameters of the satellite, including parameters such as solar pressure, from which the position vector of the satellite is written: Assume that the position of the satellite in the precise ephemeris corresponding to time t is (X o Y o Z o ), then the observation equation for dynamic parameter estimation is expressed as: The observation equation is estimated by least squares to obtain the initial position, velocity and perturbation parameter state vector correction of the satellite, thereby realizing the estimation of the dynamic parameters of the satellite orbit; The predicted orbit product of the GNSS satellite is obtained by using the dynamic parameters and the initial state of the satellite to perform orbital integral calculations; (202) Using the precise clock error product of the GNSS satellite estimated in step (1), the satellite clock error prediction is carried out based on the quadratic polynomial model. The process is as follows: The satellite clock error model is expressed as clk i =a 0 +a 1 (t-t c )+a 2 (t-t c ) 2 Where clk i is the satellite clock error at time t, t c is the reference time of the clock error model, a 0 t c The satellite clock error at time a 1 t c The clock speed parameter of the clock at this moment, a 2 t c The frequency drift parameter of the clock at this moment; Based on the preprocessed GNSS satellite precise clock error sequence, the quadratic polynomial model parameter a of the precise clock error sequence is estimated. 0 , a 1 , a 2 ; Then the GNSS satellite clock error during the forecast period is calculated using the quadratic polynomial model parameters; (203) The least squares method is used to fit the orbit of the GNSS satellite. The process is as follows: Assume that the satellite ephemeris model parameter at the reference epoch is σ(t 0 ), then the observation equation can be expressed as In the formula, is the position component of the satellite; linearizing the above equation, we get y=HΔσ In the formula, ' Δσ=σ-σ; According to the least squares estimation principle, the optimal estimate of Δσ is: Δσ=(H T H) -1 H T y Using Δσ to correct σ(t 0 ), obtain the satellite ephemeris model parameters of the GNSS satellite reference epoch time; use the quadratic polynomial model to fit the predicted GNSS clock error into the clock error parameters of the reference epoch time; The ephemeris model parameters of the GNSS satellite and the clock error parameters of the GNSS are compiled into a broadcast ephemeris, which is then uploaded to the visible low-orbit satellites using the satellite-to-ground link of the gateway station and distributed to the entire low-orbit constellation based on the inter-satellite link of the low-orbit satellites.

4. A real-time high-precision PNT service method based on space-ground joint observation resources according to claim 1, It is characterized in that The specific method of step (4) is: (401) The user-side GNSS and low-orbit satellite pseudorange and carrier phase observation equations are established as follows: In the formula, the subscript r represents the user terminal number, i and j represent the frequency identifiers of the GNSS and low-orbit satellite downlink navigation signals, respectively, the superscript GNSS and m represent the GNSS satellite system and its satellite number, respectively, and the superscript LEO and n represent the low-orbit satellite system and its satellite number, respectively; and They represent the pseudorange and carrier phase observation values ​​of the frequency of satellite i in the GNSS system m observed by the user end, and They represent the pseudorange and carrier phase observation values ​​of the frequency j of the n-th satellite in the low-orbit system observed by the user end respectively; and Represent the geometric distances between the user terminal and the GNSS satellite m and the low-orbit satellite n, respectively, and are expressed as and Where X GNSS,m , X LEO,n and X r are the position vectors of GNSS satellite, low-orbit satellite and user terminal respectively; dt r , dT GNSS,m and dT LEO,n Represent the clock errors of the user end, GNSS satellite and low-orbit satellite respectively, represents the tropospheric delay of GNSS satellite m, represents the ionospheric delay of GNSS satellite m at frequency i, represents the ionospheric delay of low-orbit satellite n at frequency j, represents the wavelength of the navigation signal of GNSS satellite m frequency i, represents the wavelength of the navigation signal of frequency j of low-orbit satellite n, and Represents the carrier phase ambiguity of the downlink navigation signal from GNSS and low-orbit satellites, Represent the measurement noise of pseudorange and carrier phase of GNSS and LEO satellite respectively; (402) Performing the user-side position and clock error solution based on the joint precise point positioning observation equation of GNSS and low-orbit satellite: Calculate the position X of GNSS and LEO satellites based on the broadcast ephemeris of GNSS and LEO satellites GNSS,m and X LEO,n , and the clock difference dT between GNSS and LEO satellites GNSS,m and dT LEO,n ; Read the pseudo-range and carrier phase observation data of GNSS and low-orbit satellites collected by the user end, and perform gross error elimination and cycle slip detection; Calculate the floating-point solution of carrier phase ambiguity of GNSS and LEO satellites based on least squares or Kalman filter estimation method; The LAMBDA method is used to fix the ambiguity parameters of GNSS and low-orbit satellites, thereby calculating the position parameter X of the user end. r and clock difference dt r .

Images