Low-orbit satellite orbit determination method considering GPS elastic power

Abstract

A low-orbit satellite orbit determination method considering GPS elastic power comprises the steps that firstly, a carrier-to-noise ratio template is used for detecting the opening time period of the elastic power of a GPS satellite, then a non-geometric combination of satellite-borne observation values in the opening and closing periods of the elastic power is extracted, and a three-time 98% quantile is used for determining a non-geometric combination (GF) epoch difference threshold value in the opening and closing periods of the elastic power; for the signal-to-noise ratio jump, a signal-to-noise ratio jump detection amount DR based on the space between satellites is constructed according to a signal power propagation formula, then the distribution condition of the DR in the elastic power starting period is counted, a signal-to-noise ratio jump detection threshold value is determined by using a triple standard deviation criterion, and finally the ambiguity continuity is detected by combining the two methods. The method is designed for change characteristics of satellite-borne observation data of the low-orbit satellite under elastic power starting; the method is suitable for real-time and afterwards low-orbit satellite orbit determination, and can effectively improve the low-orbit satellite orbit determination precision in the GPS elastic power starting period.

Description

Technical Field

[0001] The present invention relates to the field of satellite positioning and navigation GNSS data and signal processing, and in particular to a low-orbit satellite orbit determination method that takes into account GPS elastic power for applications such as satellite orbit determination. Background Art

[0002] High-precision geometric orbit determination of low-orbit satellites plays an irreplaceable role in Earth science research, including inversion of the Earth's gravity field, detection of the Earth's magnetic field, and ocean altimetry. Geometric orbit determination of low-orbit satellites uses GPS carrier phase and pseudorange observations to determine the satellite's orbit through precise single-point positioning. The quality of geometric orbit determination depends primarily on the proper processing of phase observations, which must be accurate to the millimeter level. In addition to estimating common parameters of pseudorange and phase, such as receiver coordinates and receiver clock errors, the ambiguity parameters of the phase observations must be estimated with millimeter-level accuracy. Failure to accurately estimate these ambiguity parameters can significantly increase the error in orbit determination of low-orbit satellites. Therefore, accurate estimation of these ambiguity parameters is crucial for precise geometric orbit determination.

[0003] New GPS satellites, such as Block IIR-M, IIF, and III, can transmit signals at multiple frequencies and reallocate the transmit power of one or more frequencies over a specified area. This ability to adjust signal transmission power is called elastic power. To prevent signal interference, GPS uses elastic power to boost satellite signal transmission power, thereby improving its signal anti-interference capabilities. However, when GPS elastic power is enabled, the orbit quality of low-orbit satellites (such as GRACE-FO) is significantly degraded. Therefore, improving the orbit determination accuracy of low-orbit satellites during this period is crucial.

[0004] Currently, research has focused on methods for detecting elastic power and its impact on differential code bias and precise point positioning. However, less attention has been paid to its impact on precise orbit determination of low-orbit satellites. Furthermore, there are currently no techniques for mitigating the effects of elastic power on GPS observations when dealing with LEO satellite orbit determination under elastic power. Summary of the Invention

[0005] In response to the above technical problems, the present invention proposes a low-orbit satellite orbit determination method that takes into account the GPS elastic power. It solves the problems caused by elastic power from the GPS observation value level, improves the low-orbit satellite orbit determination accuracy during the period when the elastic power is turned on, and thus solves the low-orbit satellite orbit determination problem during the period when the GPS elastic power is turned on.

[0006] Technical Solution

[0007] A method for determining the orbit of a low-orbit satellite taking into account the elastic power of GPS comprises the following steps:

[0008] Step S1: obtaining GPS observation data from a low-orbit satellite onboard receiver, reading a GPS elastic power on period file, and classifying the GPS observation data into two categories: elastic power on and elastic power off, and performing data preprocessing on the classified GPS observation data;

[0009] Step S2: construct the geometry-free combination GF according to the spaceborne observations during the elastic power on and off periods;

[0010] Step S3: Calculate the epoch difference ΔGF between the two types of data GF when the elastic power is turned on and off, and set the detection threshold of ΔGF;

[0011] A month of low-orbit satellite observation data was collected, and the epoch difference ΔGF of the two types of data GF was calculated. The 98% quantile of the ΔGF set was obtained, and three times the 98% quantile was used as the ΔGF detection threshold. The detection threshold of ΔGF during the elastic power on and off periods was obtained.

[0012] Step S4: Collect observation data of low-orbit satellites for one month, and construct a signal-to-noise ratio jump test value DR based on the onboard observation values during the elastic power on period, the signal propagation formula, the distance between the GPS satellite and the low-orbit satellite, and the onboard signal-to-noise ratio.

[0013] Step S5: Counting the standard deviation of the signal-to-noise ratio jump test quantity DR epoch difference ΔDR, and taking three times the standard deviation as the signal-to-noise ratio jump threshold;

[0014] Step S6: Determine whether the elastic power is enabled in the current epoch according to the elastic power enable period file; if so, proceed to step S7;

[0015] If not, proceed to step S8;

[0016] Step S7: Using three times the 98% quantile of the elastic power on period ΔGF as the detection threshold, detect the continuity of the ambiguity and determine whether the ΔGF of the current epoch exceeds the limit;

[0017] If yes, go to step S10;

[0018] If not, proceed to step S9;

[0019] Step S8: Using three times the 98% quantile of the elastic power off period ΔGF as the detection threshold, detect the continuity of the ambiguity and determine whether the ΔGF of the current epoch exceeds the limit;

[0020] If yes, go to step S10;

[0021] If not, go to step S11;

[0022] Step S9: Using three times the standard deviation of ΔDR as the signal-to-noise ratio jump threshold, detect the continuity of the ambiguity and determine whether the ΔDR of the current epoch exceeds the limit;

[0023] If yes, go to step S10;

[0024] If not, go to step S11;

[0025] Step S10: Go to step S12 to initialize the ambiguity parameters of the satellite corresponding to the exceeded limit;

[0026] Step S11: Enter step S12, estimate the ambiguity parameters of the corresponding satellites that are not exceeded in a constant form; Step S12: Determine the low-orbit satellite orbit based on the ambiguity parameter settings.

[0027] Furthermore, in step S1, data preprocessing includes: GPS satellite Block type screening, GPS satellite cutoff elevation angle setting, and satellite signal transmission time iterative calculation;

[0028] Specifically, the iterative calculation formula for the satellite signal transmission time is:

[0029] P=c(T r -T s ) (1)

[0030] t s =T s -dt s =T r -P / c-dt s (2)

[0031] Wherein, subscripts r and s represent receiver and GPS satellite respectively, T r and T s represents the clock time of the receiver and GPS satellite respectively, c represents the speed of light, P represents the distance between the low-orbit satellite and the GPS satellite, t s is the time when the satellite signal is transmitted, dt s Indicates the satellite clock error.

[0032] Furthermore, step S3 includes:

[0033] Step S31: Collect one month of low-orbit satellite observation data and calculate the GF epoch difference ΔGF as follows:

[0034] ΔGF k =GF k -GF k-1 (3)

[0035] Wherein, the subscripts k and k-1 represent the kth and k-1th epochs respectively, GF k and GFk-1 denote the geometry-free combination GF, ΔGF of k and k-1 epochs respectively k represents the geometry-free epoch difference of the kth epoch;

[0036] Step S32: Take the absolute value of the GF epoch difference ΔGF as follows:

[0037] |ΔGF|=abs(ΔGF) (4)

[0038] Where |ΔGF| represents the set of absolute values of the GF epoch difference ΔGF, and abs(*) represents the absolute value function;

[0039] Step S33: Take the 98% quantile of the absolute value set of the GF epoch difference ΔGF, as follows:

[0040] 98% = 100% - (100% / n) * (x 98th -L) (5)

[0041] Among them, n represents the number of data sets, L represents the minimum value in the data set, and x 98th Indicates the quantile that the data needs to determine, i.e. the 98% quantile.

[0042] Calculate the 98% quantile x of the GF epoch difference during the elastic power on and off periods respectively 98th,on and x 98th,off .

[0043] Step S34: Using three times the 98% quantile as the detection threshold for the elastic power on and off periods of ΔGF, as follows:

[0044]

[0045] Among them, x 98th,on and x 98th,off are the 98% quantiles of the GF epoch difference ΔGF during the elastic power on and off periods, T ΔGF，on and T ΔGF，off are the ΔGF detection thresholds during the elastic power on and off periods, respectively.

[0046] Furthermore, the step S4 includes:

[0047] Step S41: Based on the precise ephemeris, the GPS satellite orbit coordinates are obtained by interpolation, as follows:

[0048]

[0049] Among them, n is the number of coordinates used for interpolation, that is, the number of nodes, the subscripts i and k are the node sequence numbers, and x i and x kis the time corresponding to the i-th and k-th nodes, y k is the function value corresponding to the kth node, and x is the time corresponding to any point in the interpolation interval.

[0050] Step S42: Calculate the geometric distance between the GPS satellite and the low-orbit satellite as follows:

[0051]

[0052] Among them, x s 、y s and z s is the XYZ coordinate of the GPS satellite, x r 、y r and z r The XYZ coordinates of the low-orbit satellite onboard receiver; GPS satellite coordinates are obtained through precise ephemeris, and the coordinates of the low-orbit satellite onboard receiver are obtained through single point positioning (SPP);

[0053] Step S43: Count the signal-to-noise ratio jump test value DR of the low-orbit satellite onboard observation data for one month. The signal-to-noise ratio jump test value DR is defined as follows:

[0054] DR=lg(d·SNR)(9)

[0055] Where d is the geometric distance between the GPS satellite and the low-orbit satellite, and SNR is the signal-to-noise ratio of the low-orbit satellite onboard receiver.

[0056] Furthermore, the step S5 includes:

[0057] Step S51: Calculate the epoch difference ΔDR of the signal-to-noise ratio jump test value DR as follows:

[0058] ΔDR k =DR k -DR k-1 (10)

[0059] Wherein, the subscripts k and k-1 represent the kth and k-1th epochs respectively, DR k and DR k-1 DR and ΔDR represent the signal-to-noise ratio jump test quantities of k and k-1 epochs respectively k The DR epoch difference represents the signal-to-noise ratio jump test quantity of the kth epoch;

[0060] Step S52: Calculate the average value of the signal-to-noise ratio jump test quantity DR epoch difference ΔDR set as follows:

[0061]

[0062] Wherein, subscript k is the number of the element in the ΔDR set, ΔDR k is the kth ΔDR, n is the number of elements in the ΔDR set, is the average value of the ΔDR set;

[0063] Step S53: Calculate the standard deviation of the ΔDR set as follows:

[0064]

[0065] Among them, σ ΔDR is the standard deviation of the ΔDR set, subscript k is the number of the element in the ΔDR set, ΔDR k is the kth ΔDR, n is the number of elements in the ΔDR set, is the average value of the ΔDR set;

[0066] Step S54: Using three times the standard deviation as the signal-to-noise ratio transition threshold, as follows:

[0067] T ΔDR =3·σ ΔDR (13)

[0068] Among them, σ ΔDR is the standard deviation of the ΔDR set, T ΔDR is the ΔDR detection threshold.

[0069] Furthermore, the step S12 includes:

[0070] Step S121: setting the variance-covariance matrix of the fuzzy parameters;

[0071]

[0072] Among them, Q amb is the variance of the ambiguity parameter. If the ambiguity needs to be initialized, the variance of the ambiguity parameter is set to (60m) 2 ; If the ambiguity is estimated as a constant, the variance of the ambiguity parameters will not be changed.

[0073] Step S122: construct the ionospheric-free combination of pseudorange and phase as follows:

[0074]

[0075] Among them, the superscript s represents the satellite, and the subscript r represents the satellite-borne receiver. and are the pseudorange and phase ionospheric-free combinations, is the distance from the GPS satellite mass center to the low-orbit satellite mass center; δt r and δt s are the receiver and satellite clock errors, c is the speed of light, Br，IF and are the ionospheric-free pseudorange hardware delays of the receiver and satellite, respectively, and D r，IF and are the ionospheric-free phase hardware delays of the receiver and satellite, respectively, IF and are the ionospheric-free combined wavelength and ionospheric-free combined ambiguity, respectively. and are the ionospheric-free combined pseudorange and phase observation noises, respectively.

[0076] Step S123: Kalman filter solves the parameters as follows:

[0077]

[0078] Among them, H k is the observation matrix, R k is the observed value variance-covariance matrix, K k is the Kalman filter gain matrix, and are the prior and posterior variance-covariance matrices of the parameters, and are parameter prediction and estimation vectors, respectively, z k is the observation vector, and I is the unit matrix. The observation values are obtained through the onboard observation data, and the observation matrix H k By linearizing the observation equation, R k By constructing a random model of altitude angle, The variance of the initial parameters is set to (60m) 2 , the parameters x to be estimated include the coordinates of the low-orbit satellite onboard receiver, the receiver clock error and the ambiguity parameters.

[0079] In summary, the present invention classifies the satellite observation values according to the GPS elastic power on period file, constructs the geometric combination epoch difference and determines its detection threshold, and establishes the signal-to-noise ratio jump detection model based on the signal power propagation formula, thereby realizing high-precision orbit determination of low-orbit satellites during the elastic power on period.

[0080] Beneficial effects

[0081] This invention addresses the issues caused by elastic power at the GPS observation level. This method rationally processes carrier phase observations, accurately checks the continuity of GPS ambiguity under elastic power, and improves the orbit determination accuracy of low-orbit satellites during the elastic power period. This improved orbit determination accuracy during the elastic power period has significant application value for LEO satellite orbit determination under elastic power. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1This is a flow chart of a method for determining the orbit of a low-orbit satellite taking into account GPS elastic power according to an embodiment of the present invention;

[0083] Figure 2 This is a schematic diagram of a specific process of step S3 in an embodiment of the present invention;

[0084] Figure 3 This is a schematic diagram of a specific process of step S4 in an embodiment of the present invention;

[0085] Figure 4 This is a schematic diagram of a specific process of step S5 in an embodiment of the present invention;

[0086] Figure 5 3 is a schematic diagram of a specific flow chart of step S12 in an embodiment of the present invention. DETAILED DESCRIPTION

[0087] The following describes a specific embodiment of the present invention in more detail with reference to schematic diagrams. The advantages and features of the present invention will become more apparent from the following description and claims. It should be noted that the drawings are greatly simplified and not to exact scale, and are intended solely to facilitate and clearly illustrate the embodiments of the present invention.

[0088] To ensure high-precision orbit determination of low-orbit satellites during the period when elastic power is on, the present invention proposes a method for orbit determination of low-orbit satellites that takes into account GPS elastic power. The basic idea is: first, read the GPS elastic power on period file, determine whether the onboard observation value is enhanced by the elastic power signal based on the file, and then divide the onboard observation value into two categories: elastic power on period and elastic power off period. Use the onboard observation data to construct the geometry-free combined epoch difference ΔGF, and calculate the 98% quantile of the geometry-free combined epoch difference of the two types of data. Three times the 98% quantile is used as the detection threshold T for the two types of data. ΔGF，on and T ΔGF，off The GPS satellite orbital coordinates are calculated using precise ephemeris, and then the geometric distance between the GPS satellite and the low-orbit satellite is calculated. The signal-to-noise ratio jump test quantity DR is constructed by the inter-satellite geometric distance and the signal-to-noise ratio of the onboard observation data, and the epoch difference ΔDR of the signal-to-noise ratio jump test quantity is calculated. The standard deviation of the epoch difference of the signal-to-noise ratio jump test quantity is calculated, and three times the standard deviation is used as the detection threshold T of the signal-to-noise ratio jump test quantity. ΔDR Based on the above two types of indicators, the low-orbit satellite onboard observation data is tested. When the elastic power is turned off, the T ΔGF，off and the ΔGF of the current epoch, if ΔGF is greater than T ΔGF,off , then initialize the fuzzy parameters; if ΔGF is less than T ΔGF,off , then the fuzzy parameters are estimated in constant form. When the elastic power is turned on, the T obtained during the elastic power on period is compared. ΔGF,on and the ΔGF of the current epoch, if ΔGF is greater than TΔGF，on , then initialize the fuzzy parameters; if ΔGF is less than T ΔGF，on , compared with T ΔDR and the ΔDR of the current epoch, if ΔDR is greater than T ΔDR , then initialize the ambiguity parameters; if ΔDR is less than T ΔDR , the ambiguity parameters are estimated as constants.

[0089] refer to Figure 1 A method for determining the orbit of a low-orbit satellite taking into account the elastic power of GPS comprises the following steps:

[0090] Step S1:

[0091] Obtain GPS observation data from low-orbit satellite onboard receivers, read the GPS elastic power on period file, and classify the GPS observation data into two categories: elastic power on and elastic power off. Perform data preprocessing on the classified GPS observation data.

[0092] Data preprocessing includes but is not limited to: satellite block type screening, satellite cutoff elevation angle setting, and iterative calculation of satellite signal transmission time;

[0093] Preferably, the satellite signal transmission time iteration formula used is:

[0094] P=c(T r -T s ) (1)

[0095] t s =T s -dt s =T r -P / c-dt s (2)

[0096] Wherein, subscripts r and s represent receiver and satellite respectively, T r and T s represents the clock time of the receiver and the satellite respectively, c represents the speed of light, P represents the distance between the satellite and the earth, t s is the time when the satellite signal is transmitted, dt s Indicates the satellite clock error.

[0097] Step S2: construct the geometry-free combination GF according to the satellite observation values during the elastic power on and off periods;

[0098] Step S3: Calculate the epoch difference ΔGF between the two types of data GF, obtain the 98% quantile of the ΔGF set, and use three times the 98% quantile as the detection threshold of ΔGF;

[0099] Specifically, refer to Figure 2 , step S3 includes:

[0100] Step S31: Collect one month of low-orbit satellite observation data and calculate the GF epoch difference ΔGF as follows:

[0101] ΔGF k =GF k -GF k-1 (3)

[0102] Wherein, the subscripts k and k-1 represent the kth and k-1th epochs respectively, GF k and GF k-1 denote the geometry-free combination GF, ΔGF of k and k-1 epochs respectively k represents the geometry-free epoch difference of the kth epoch;

[0103] Step S32: Take the absolute value of the GF epoch difference ΔGF as follows:

[0104] |ΔGF|=abs(ΔGF) (4)

[0105] Where |ΔGF| represents the set of absolute values of the GF epoch difference ΔGF, and abs(*) represents the absolute value function;

[0106] Step S33: Take the 98% quantile of the absolute value set of the GF epoch difference ΔGF, as follows:

[0107] 98% = 100% - (100% / n) * (x 98th -L) (5)

[0108] Among them, n represents the number of data sets, L represents the minimum value in the data set, and x 98th Indicates the quantile that the data needs to determine, i.e. the 98% quantile.

[0109] Calculate the 98% quantile x of the GF epoch difference during the elastic power on and off periods respectively 98th,on and x 98th,off .

[0110] Step S34: Using three times the 98% quantile as the detection threshold for the elastic power on and off periods of ΔGF, as follows:

[0111]

[0112] Among them, x 98th,on and x 98th，off are the 98% quantiles of the GF epoch difference ΔGF during the elastic power on and off periods, T ΔGF，on and T ΔGF，off are the ΔGF detection thresholds during the elastic power on and off periods, respectively.

[0113] Step S4: Based on the onboard observation values during the elastic power on period, based on the signal propagation formula, using the distance between the GPS satellite and the low-orbit satellite and the onboard signal noise ratio, construct the signal-to-noise ratio jump test value DR;

[0114] Specifically, refer to Figure 3 , step S4 includes:

[0115] Step S41: Based on the precise ephemeris, the GPS satellite orbit coordinates are obtained by interpolation, as follows:

[0116]

[0117] Among them, n is the number of coordinates used for interpolation, that is, the number of nodes, the subscripts i and k are the node sequence numbers, and x i and x k is the time corresponding to the i-th and k-th nodes, y k is the function value corresponding to the kth node, and x is the time corresponding to any point in the interpolation interval.

[0118] Step S42: Calculate the geometric distance between the GPS satellite and the low-orbit satellite as follows:

[0119]

[0120] Among them, x s 、y s and z s is the XYZ coordinate of the GPS satellite, x r 、y r and z r The XYZ coordinates of the low-orbit satellite onboard receiver; GPS satellite coordinates are obtained through precise ephemeris, and the coordinates of the low-orbit satellite onboard receiver are obtained through single point positioning (SPP);

[0121] Step S43: Count the signal-to-noise ratio jump test value DR of the low-orbit satellite onboard observation data for one month. The signal-to-noise ratio jump test value DR is defined as follows:

[0122] DR=lg(d·SNR) (9)

[0123] Where d is the geometric distance between the GPS satellite and the low-orbit satellite, and SNR is the signal-to-noise ratio of the low-orbit satellite onboard receiver.

[0124] Step S5: Counting the standard deviation of the signal-to-noise ratio jump test quantity DR epoch difference ΔDR, and taking three times the standard deviation as the signal-to-noise ratio jump threshold;

[0125] Specifically, refer to Figure 4 , step S5 includes:

[0126] Step S51: Calculate the epoch difference ΔDR of the signal-to-noise ratio jump test value DR as follows:

[0127] ΔDR k =DR k -DR k-1 (10)

[0128] Wherein, the subscripts k and k-1 represent the kth and k-1th epochs respectively, DR k and DR k-1 DR and ΔDR represent the signal-to-noise ratio jump test quantities of k and k-1 epochs respectively k The DR epoch difference represents the signal-to-noise ratio jump test quantity of the kth epoch;

[0129] Step S52: Calculate the average value of the signal-to-noise ratio jump test quantity DR epoch difference ΔDR set as follows:

[0130]

[0131] Wherein, subscript k is the number of the element in the ΔDR set, ΔDR k is the kth ΔDR, n is the number of elements in the ΔDR set, is the average value of the ΔDR set;

[0132] Step S53: Calculate the standard deviation of the ΔDR set as follows:

[0133]

[0134] Among them, σ ΔDR is the standard deviation of the ΔDR set, subscript k is the number of the element in the ΔDR set, ΔDR k is the kth ΔDR, n is the number of elements in the ΔDR set, is the average value of the ΔDR set;

[0135] Step S54: Using three times the standard deviation as the signal-to-noise ratio transition threshold, as follows:

[0136] T ΔDR =3·σ ΔDR (13)

[0137] Among them, σ ΔDR is the standard deviation of the ΔDR set, T ΔDR is the ΔDR detection threshold.

[0138] Step S6: Determine whether the elastic power is enabled in the current epoch according to the elastic power enable period file; if so, proceed to step S7;

[0139] If not, proceed to step S8;

[0140] Step S7: Using three times the 98% quantile of the elastic power on period ΔGF as the detection threshold, detect the continuity of the ambiguity and determine whether the ΔGF of the current epoch exceeds the limit;

[0141] If yes, go to step S10;

[0142] If not, proceed to step S9;

[0143] Step S8: Using three times the 98% quantile of the elastic power off period ΔGF as the detection threshold, detect the continuity of the ambiguity and determine whether the ΔGF of the current epoch exceeds the limit;

[0144] If yes, go to step S10;

[0145] If not, go to step S11;

[0146] Step S9: Using three times the standard deviation of ΔDR as the signal-to-noise ratio jump threshold, detect the continuity of the ambiguity and determine whether the ΔDR of the current epoch exceeds the limit;

[0147] If yes, go to step S10;

[0148] If not, go to step S11;

[0149] Step S10: Go to step S12 to initialize the ambiguity parameters of the satellite corresponding to the exceeded limit;

[0150] Step S11: Enter step S12, estimate the ambiguity parameters of the corresponding satellites that are not exceeded in a constant form; Step S12: Determine the low-orbit satellite orbit based on the ambiguity parameter settings.

[0151] Specifically, refer to Figure 5 , step S12 includes:

[0152] Step S121: Setting the variance-covariance matrix of the fuzziness parameters, which can be divided into two cases:

[0153]

[0154] Among them, Q amb is the variance of the ambiguity parameter. If the ambiguity needs to be initialized, the variance of the ambiguity parameter is set to (60m) 2 ; If the ambiguity is estimated as a constant, the variance of the ambiguity parameters will not be changed.

[0155] Step S122: construct the ionospheric-free combination of pseudorange and phase as follows:

[0156]

[0157] Among them, the superscript s represents the satellite, and the subscript r represents the satellite-borne receiver. and are the pseudorange and phase ionospheric-free combinations, is the distance from the GNSS satellite center of mass to the LEO satellite center of mass; δt r and δt s are the receiver and satellite clock errors, c is the speed of light, B r，IF and are the ionospheric-free pseudorange hardware delays of the receiver and satellite, respectively, and D r，IF and are the ionospheric-free phase hardware delays of the receiver and satellite, respectively, IF and are the ionospheric-free combined wavelength and ionospheric-free combined ambiguity, respectively. and are the ionospheric-free combined pseudorange and phase observation noises, respectively.

[0158] Step S123: Kalman filter solves the parameters as follows:

[0159]

[0160] Among them, H k is the observation matrix, R k is the observation variance-covariance matrix, K k is the Kalman filter gain matrix, and are the prior and posterior variance-covariance matrices of the parameters, and are parameter prediction and estimation vectors, respectively, z k is the observation vector, and I is the unit matrix. The observation values are obtained through the onboard observation data, and the observation matrix H k By linearizing the observation equation, R k By constructing a random model of altitude angle, The variance of the initial parameters is set to (60m) 2 , the parameters x to be estimated include the coordinates of the low-orbit satellite onboard receiver, the receiver clock error and the ambiguity parameters.

[0161] The processing scenarios and processing methods corresponding to steps S6-S12 are shown in the following table.

[0162]

[0163]

[0164] In summary, the present invention classifies the satellite observation values according to the GPS elastic power on period file, constructs the geometric combination epoch difference and determines its detection threshold, and establishes the signal-to-noise ratio jump detection model based on the signal power propagation formula, thereby realizing high-precision orbit determination of low-orbit satellites during the elastic power on period.

[0165] Using one month of onboard data, the orbit determination errors of low-orbit satellites using the present invention are compared with those without, as shown in the table below. With the present invention, the 3D orbit determination errors of both GRACE-C and GRACE-D LEO satellites are less than 5 cm, representing improvements of 36% and 21%, respectively, compared to results without the present invention.

[0166]

[0167] The above description is only a description of the preferred embodiments of the present application and does not limit the scope of the present application. Any changes or modifications made by any person skilled in the art based on the above disclosed technical content should be regarded as equivalent valid embodiments and fall within the scope of protection of the technical solution of the present application.

Claims

1. A method for determining the orbit of a low-orbit satellite taking into account the elastic power of GPS, characterized in that: The following steps are involved: Step S1: obtaining GPS observation data from a low-orbit satellite onboard receiver, reading a GPS elastic power on period file, and classifying the GPS observation data into two categories: elastic power on and elastic power off, and performing data preprocessing on the classified GPS observation data; Step S2: construct the geometry-free combination GF according to the spaceborne observations during the elastic power on and off periods; Step S3: Calculate the epoch difference ΔGF between the two types of data GF when the elastic power is turned on and off, and set the detection threshold of ΔGF; A month of low-orbit satellite observation data was collected, and the epoch difference ΔGF of the two types of data GF was calculated. The 98% quantile of the ΔGF set was obtained, and three times the 98% quantile was used as the ΔGF detection threshold. The detection threshold of ΔGF during the elastic power on and off periods was obtained. Step S4: Collect observation data of low-orbit satellites for one month, and construct a signal-to-noise ratio jump test value DR based on the onboard observation values during the elastic power on period, the signal propagation formula, the distance between the GPS satellite and the low-orbit satellite, and the onboard signal-to-noise ratio. Step S5: Counting the standard deviation of the signal-to-noise ratio jump test quantity DR epoch difference ΔDR, and taking three times the standard deviation as the signal-to-noise ratio jump threshold; Step S6: judging whether the elastic power is enabled in the current epoch according to the elastic power enable period file; If yes, go to step S7; If not, proceed to step S8; Step S7: Using three times the 98% quantile of the elastic power on period ΔGF as the detection threshold, detect the continuity of the ambiguity and determine whether the ΔGF of the current epoch exceeds the limit; If yes, go to step S10; If not, proceed to step S9; Step S8: Using three times the 98% quantile of the elastic power off period ΔGF as the detection threshold, detect the continuity of the ambiguity and determine whether the ΔGF of the current epoch exceeds the limit; If yes, go to step S10; If not, go to step S11; Step S9: Using three times the standard deviation of ΔDR as the signal-to-noise ratio jump threshold, detect the continuity of the ambiguity and determine whether the ΔDR of the current epoch exceeds the limit; If yes, go to step S10; If not, go to step S11; Step S10: Go to step S12 to initialize the ambiguity parameters of the satellite corresponding to the exceeded limit; Step S11: Go to step S12, estimate the ambiguity parameters of the corresponding satellites that are not exceeded in a constant form; Step S12: Determine the low-orbit satellite orbit according to the ambiguity parameter settings.

2. A method for determining the orbit of a low-orbit satellite taking into account the elastic power of GPS according to claim 1, characterized in that: In step S1, the GPS observation data preprocessing includes: satellite Block type screening, satellite cutoff elevation angle setting and satellite signal transmission time iterative calculation; The iterative calculation formula for satellite signal transmission time is: P=c(T r -T s ) (1) t s =T s -dt s =T r -P / c-dt s (2) Wherein, subscripts r and s represent receiver and GPS satellite respectively, T r and T s represents the clock time of the receiver and GPS satellite respectively, c represents the speed of light, P represents the distance between the low-orbit satellite and the GPS satellite, t s is the time when the satellite signal is transmitted, dt s Indicates the satellite clock error.

3. A method for determining the orbit of a low-orbit satellite taking into account the elastic power of GPS according to claim 1, characterized in that: The step S3 comprises: Step S31: Collect one month of low-orbit satellite observation data and calculate the GF epoch difference ΔGF as follows: ΔGF k =GF k -GF k-1 (3) Wherein, the subscripts k and k-1 represent the kth and k-1th epochs respectively, GF k and GF k-1 denote the geometry-free combination GF, ΔGF of k and k-1 epochs respectively k represents the geometry-free epoch difference of the kth epoch; Step S32: Take the absolute value of the GF epoch difference ΔGF as follows: |ΔGF|=abs(ΔGF) (4) Where |ΔGF| represents the set of absolute values of the GF epoch difference ΔGF, and abs(*) represents the absolute value function; Step S33: Take the 98% quantile of the absolute value set of the GF epoch difference ΔGF, as follows: 98%=100%-(100% / n)*(x 98th -L) (5) Among them, n represents the number of data sets, L represents the minimum value in the data set, and x 98th Indicates the quantile that the data needs to determine, that is, the 98% quantile; Calculate the 98% quantile x of the GF epoch difference during the elastic power on and off periods respectively 98th,on and x 98th,off ; Step S34: Using three times the 98% quantile as the detection threshold for the elastic power on and off periods of ΔGF, as follows: Among them, x 98th,on and x 98th,off are the 98% quantiles of the GF epoch difference ΔGF during the elastic power on and off periods, T ΔGF，on and T ΔGF，off are the ΔGF detection thresholds during the elastic power on and off periods, respectively.

4. A method for determining the orbit of a low-orbit satellite taking into account the elastic power of GPS according to claim 1, characterized in that: The step S4 comprises: Step S41: Based on the precise ephemeris, the GPS satellite orbit coordinates are obtained by interpolation, as follows: Among them, n is the number of coordinates used for interpolation, that is, the number of nodes, the subscripts i and k are the node sequence numbers, and x i and x k is the time corresponding to the i-th and k-th nodes, y k is the function value corresponding to the kth node, and x is the time corresponding to any point in the interpolation interval; Step S42: Calculate the geometric distance between the GPS satellite and the low-orbit satellite as follows: Among them, x s 、y s and z s is the XYZ coordinate of the GPS satellite, x r 、y r and z r is the XYZ coordinate of the low-orbit satellite onboard receiver; GPS satellite coordinates are obtained through precise ephemeris, and the low-orbit satellite onboard receiver coordinates are obtained through single point positioning (SPP); Step S43: Counting the signal-to-noise ratio jump test value DR of the low-orbit satellite onboard observation data for one month; the signal-to-noise ratio jump test value DR is defined as follows: DR=lg(d·SNR) (9) Where d is the geometric distance between the GPS satellite and the low-orbit satellite, and SNR is the signal-to-noise ratio of the low-orbit satellite onboard receiver.

5. The method for determining the orbit of a low-orbit satellite taking into account the elastic power of GPS according to claim 1, wherein: The step S5 comprises: Step S51: Calculate the epoch difference ΔDR of the signal-to-noise ratio jump test value DR as follows: ΔDR k =DR k -DR k-1 (10) Wherein, the subscripts k and k-1 represent the kth and k-1th epochs respectively, DR k and DR k-1 DR and ΔDR represent the signal-to-noise ratio jump test quantities of k and k-1 epochs respectively k The DR epoch difference represents the signal-to-noise ratio jump test quantity of the kth epoch; Step S52: Calculate the average value of the signal-to-noise ratio jump test quantity DR epoch difference ΔDR set as follows: Wherein, subscript k is the number of the element in the ΔDR set, ΔDR k is the kth ΔDR, n is the number of elements in the ΔDR set, is the average value of the ΔDR set; Step S53: Calculate the standard deviation of the ΔDR set as follows: Among them, σ ΔDR is the standard deviation of the ΔDR set, subscript k is the number of the element in the ΔDR set, ΔDR k is the kth ΔDR, n is the number of elements in the ΔDR set, is the average value of the ΔDR set; Step S54: Using three times the standard deviation as the signal-to-noise ratio transition threshold, as follows: T ΔDR =3·s ΔDR (13) Among them, σ ΔDR is the standard deviation of the ΔDR set, T ΔDR is the ΔDR detection threshold.

6. A method for determining the orbit of a low-orbit satellite taking into account the elastic power of GPS according to claim 1, characterized in that: The step S12 includes: Step S121: setting the variance-covariance matrix of the fuzzy parameters; Among them, Q amb is the variance of the ambiguity parameter. If the ambiguity needs to be initialized, the variance of the ambiguity parameter is set to (60m) 2 ; If the ambiguity is estimated as a constant, the variance of the ambiguity parameter does not change; Step S122: construct the ionospheric-free combination of pseudorange and phase as follows: Among them, the superscript s represents the satellite, and the subscript r represents the satellite-borne receiver. and are the pseudorange and phase ionospheric-free combinations, is the distance from the GNSS satellite center of mass to the LEO satellite center of mass; δt r and δt s are the receiver and satellite clock errors, c is the speed of light, B r，IF and are the ionospheric-free pseudorange hardware delays of the receiver and satellite, respectively, and D r，IF and are the ionospheric-free phase hardware delays of the receiver and satellite, respectively, IF and are the ionospheric-free combined wavelength and ionospheric-free combined ambiguity, respectively. and are the ionospheric-free combined pseudorange and phase observation noises, respectively; Step S123: Kalman filter solves the parameters as follows: Among them, H k is the observation matrix, R k is the observed value variance-covariance matrix, K k is the Kalman filter gain matrix, and are the prior and posterior variance-covariance matrices of the parameters, and are parameter prediction and estimation vectors, respectively, z k is the observation value vector, I is the unit matrix; the observation value is obtained through the onboard observation data, and the observation matrix H k By linearizing the observation equation, R k By constructing a random model of altitude angle, The variance of the initial parameters is set to (60m) 2 , the parameters x to be estimated include the coordinates of the low-orbit satellite onboard receiver, the receiver clock error and the ambiguity parameters.

Images