Satellite orbit position solution and error evaluation method and system

Abstract

This invention discloses a method and system for satellite orbital position calculation and error assessment. The method includes: acquiring ephemeris parameters; obtaining the satellite's three-dimensional coordinates in the BeiDou coordinate system based on the ephemeris parameters; converting the three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Geodetic Reference Frame (IGFR) coordinate system; acquiring the true three-dimensional coordinates; calculating the three-dimensional position deviation vector between the satellite's three-dimensional coordinates in the IGFR coordinate system and the true three-dimensional coordinates in the IGFR coordinate system; obtaining the satellite's radial orbit deviation, tangential orbit deviation, and normal orbit deviation based on the three-dimensional position deviation vector; acquiring the true clock error value from a ground tracking observation station; acquiring the clock error from ephemeris data; calculating the clock error deviation between the clock error and the true clock error value; and calculating the space signal ranging error based on the clock error deviation, radial orbit deviation, tangential orbit deviation, and normal orbit deviation.

Description

Technical Field

[0001] This invention relates to the field of satellite orbital position calculation and error assessment. More specifically, it relates to a method and system for satellite orbital position calculation and error assessment. Background Technology

[0002] Currently, with the improvement and iteration of global satellite navigation systems, countries in various fields such as national security, power, network communication, and navigation are gradually becoming reliant on the various services provided by satellite navigation systems. The improvement in service performance has a positive feedback effect on scenarios requiring high-precision timekeeping. In scenarios requiring remote time comparison and tracing, satellite common-view comparison technology, as a mature high-precision time and frequency transmission technology, works by synchronously observing the same navigation satellite from two remote observation stations. By eliminating common errors such as satellite clock bias and orbital errors, high-precision comparison of the clocks at the two locations is achieved. In traditional satellite common-view comparison processing, the accuracy of satellite orbital position calculation has a significant impact on the comparison of clocks at two locations. Typically, after the receiver obtains the broadcast ephemeris information broadcast by the satellite, it uses a 16-parameter broadcast ephemeris model based on Kepler elements for calculation. The Galileo system and the BDS system both use similar models. This model has a user ranging error of approximately 5 cm for fitting the orbit of a medium Earth orbit satellite. Therefore, it is urgent to further improve the accuracy of satellite orbit position calculation, taking into account the long-period term of orbital changes, thereby reducing calculation errors and improving the accuracy of clock comparison between two locations, thus achieving higher precision. Summary of the Invention

[0003] The purpose of this invention is to provide a method and system for satellite orbital position calculation and error assessment, so as to solve at least one of the problems existing in the prior art.

[0004] To achieve the above objectives, the present invention adopts the following technical solution: The first aspect of this invention provides a method for satellite orbital position calculation and error assessment, the method comprising: Obtain satellite ephemeris data, extract ephemeris parameters from the ephemeris data, and obtain the satellite's three-dimensional coordinates in the BeiDou coordinate system based on the ephemeris parameters; Transform the satellite's three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Earth Reference Frame coordinate system; The true three-dimensional coordinates of the satellite in the coordinate system of the International Earth Reference Frame are obtained from the ground tracking observation station. The three-dimensional position deviation vector between the three-dimensional coordinates of the satellite in the coordinate system of the International Earth Reference Frame and the true three-dimensional coordinates of the satellite in the coordinate system of the International Earth Reference Frame is calculated. Based on the three-dimensional position deviation vector, the radial orbit deviation, tangential orbit deviation, and normal orbit deviation of the satellite are obtained. The clock error is obtained from the ground tracking observation station, the clock error is obtained from the ephemeris data, the clock error deviation between the clock error and the clock error true value is calculated, and the space signal ranging error is calculated based on the clock error deviation, radial orbit deviation, tangential orbit deviation and normal orbit deviation.

[0005] Optionally, the ephemeris parameters include: ephemeris reference time, orbital semi-major axis, orbital eccentricity, orbital inclination, right ascension of the ascending node, perigee angular distance, mean perigee angle at the ephemeris reference time, mean angular velocity perturbation, rate of change of orbital inclination, deviation of the rate of change of right ascension of the ascending node, sine harmonic term of latitude argument, cosine harmonic term of latitude argument, sine harmonic term of orbital radius, cosine harmonic term of orbital radius, sine harmonic term of orbital inclination, cosine harmonic term of orbital inclination, rate of change of orbital semi-major axis, and second-order rate of change of right ascension of the ascending node.

[0006] Optionally, obtaining the satellite's three-dimensional coordinates in the BeiDou coordinate system based on ephemeris parameters includes: Calculate the time difference between the signal transmission time and the ephemeris reference time; The orbital semi-major axis at the current moment is calculated using the orbital semi-major axis, time difference, and rate of change of orbital semi-major axis. The average angular velocity at the current moment is calculated using the orbital semi-major axis at the current moment and the average angular velocity perturbation. The mean anomaly is calculated using the mean anomaly at the ephemeris reference time, the time difference, and the average angular velocity at the current time. The deviated anomaly is calculated using the mean anomaly, and the true anomaly is calculated using the deviated anomaly. Calculate the angular distance to the ascending node using the angular distance to the perigee and the true anterior angle; The correction amount for the ascending node distance is calculated using the sine and cosine harmonic terms of the latitude argument; the correction amount for the satellite radius length is calculated using the sine and cosine harmonic terms of the orbital radius; the correction amount for the orbital inclination is calculated using the sine and cosine harmonic terms of the orbital inclination; the corrected ascending node distance is calculated using the ascending node distance and its correction amount; the reference value for the satellite radius is calculated using the current orbital semi-major axis, orbital eccentricity, and perihelion angle; the corrected satellite radius length is calculated using the reference value for the satellite radius and its correction amount; the reference value for the orbital inclination is calculated using the orbital inclination, the rate of change of orbital inclination, and the time difference; and the corrected orbital inclination is calculated using the reference value for the orbital inclination and its correction amount. The polar coordinates of the satellite are obtained using the corrected ascending node angular distance and the corrected satellite radius vector length, and then the polar coordinates of the satellite are converted into two-dimensional coordinates of the satellite in the orbital plane. The corrected right ascension of the ascending node is calculated using the right ascension of the ascending node, the deviation of the rate of change of the right ascension of the ascending node, the second-order rate of change of the right ascension of the ascending node, and the time difference. The orbital inclination is corrected by rotating the two-dimensional coordinates around the x-axis and then the corrected right ascension of the ascending node is rotated around the z-axis to obtain the three-dimensional coordinates of the satellite in the BeiDou coordinate system.

[0007] Optionally, the calculation of the orbital semi-major axis at the current moment using the orbital semi-major axis, time difference, and rate of change of the orbital semi-major axis includes: The first calculation result is obtained by multiplying the time difference by the rate of change of the orbit's semi-major axis. The semi-major axis of the orbit is obtained by summing the result of the first calculation with the semi-major axis of the orbit at the current moment.

[0008] Optionally, the calculation of the corrected right ascension of the ascending node using the right ascension of the ascending node, the deviation of the rate of change of the right ascension of the ascending node, the second-order rate of change of the right ascension of the ascending node, and the time difference includes: The second calculation result is obtained by multiplying the second rate of change of right ascension of half the ascending node by the square of the time difference; The third calculation result is obtained by multiplying the deviation of the rate of change of right ascension of the ascending node with the time difference; The corrected right ascension of the ascending node is obtained by summing the second and third calculation results.

[0009] Optionally, the transformation of the satellite's three-dimensional coordinates in the BeiDou coordinate system to its three-dimensional coordinates in the International Geodetic Reference Frame (IGFR) includes: Select at least three locations with their three-dimensional coordinates in the BeiDou coordinate system and their three-dimensional coordinates in the International Earth Reference Frame coordinate system, respectively; Based on the three-dimensional coordinates of at least three locations in the BeiDou coordinate system and the coordinate system of the International Earth Reference Frame, the optimal set of parameters in the Bursa seven-parameter coordinate transformation model is calculated using the least squares method, and the final Bursa seven-parameter coordinate transformation model is obtained. The final Bursa seven-parameter coordinate transformation model is used to transform the satellite's three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Earth Reference Frame coordinate system.

[0010] Optionally, obtaining the satellite's radial orbital deviation, tangential orbital deviation, and normal orbital deviation based on the three-dimensional position deviation vector includes: Establish the satellite's orbital coordinate system, which includes radial, tangential, and normal coordinates; The radial orbital deviation, tangential orbital deviation, and normal orbital deviation are obtained by projecting the three-dimensional position deviation vector onto the radial, tangential, and normal coordinates of the orbital coordinate system.

[0011] Optionally, the calculation of the spatial signal ranging error based on clock bias, radial orbital bias, tangential orbital bias, and normal orbital bias includes: The distance error is obtained by calculating the product of clock error and the speed of light. The spatial signal ranging error is calculated based on the distance error, radial orbit deviation, tangential orbit deviation, and normal orbit deviation.

[0012] Optionally, the calculation of the spatial signal ranging error based on distance error, radial orbit deviation, tangential orbit deviation, and normal orbit deviation includes: The fourth calculation result is obtained by multiplying the first weight by the radial trajectory deviation. The fifth calculation result is obtained by calculating the difference between the fourth calculation result and the distance error. The sixth calculation result is obtained by summing the squares of the tangential trajectory deviation and the normal trajectory deviation. The sixth calculation result is divided by the second weight to obtain the seventh calculation result; The sum of the fifth and seventh calculation results is used to obtain the eighth calculation result. The spatial signal ranging error is obtained by calculating the square root of the eighth calculation result.

[0013] A second aspect of the present invention provides a satellite orbital position calculation and error assessment system, the system comprising: The data acquisition module is used to acquire satellite ephemeris data, extract ephemeris parameters from the ephemeris data, and obtain the satellite's three-dimensional coordinates in the BeiDou coordinate system based on the ephemeris parameters. The coordinate transformation module is used to transform the satellite's three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Earth Reference Frame coordinate system. The orbital deviation calculation module is used to obtain the satellite's true three-dimensional coordinates in the coordinate system of the International Earth Reference Frame from the ground tracking observation station, calculate the three-dimensional position deviation vector between the satellite's three-dimensional coordinates in the coordinate system of the International Earth Reference Frame and the satellite's true three-dimensional coordinates in the coordinate system of the International Earth Reference Frame, and obtain the satellite's radial orbital deviation, tangential orbital deviation, and normal orbital deviation based on the three-dimensional position deviation vector. The space signal ranging error calculation module is used to obtain the true value of clock error from the ground tracking observation station, obtain the clock error from the ephemeris data, calculate the clock error deviation between the clock error and the true value of clock error, and calculate the space signal ranging error based on the clock error deviation, radial orbit deviation, tangential orbit deviation, and normal orbit deviation.

[0014] The beneficial effects of this invention are as follows: The technical solution of this invention provides a method for satellite orbit position calculation and error assessment based on the broadcast ephemeris 18 parameter model. This method can further improve the accuracy of satellite orbit position calculation by taking into account the long-period term of orbit change, thereby reducing calculation error and improving the accuracy of clock comparison between two locations, thus achieving higher accuracy, solving the calculation error caused by orbit change, and simultaneously assessing the error of the calculation results. Attached Figure Description

[0015] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings.

[0016] Figure 1 A flowchart illustrating the satellite orbit position calculation and error assessment method according to an embodiment of the present invention is shown.

[0017] Figure 2 This diagram illustrates a satellite orbit position calculation and error assessment system according to an embodiment of the present invention.

[0018] Figure 3 A flowchart illustrating another embodiment of the present invention, a method for satellite orbit position calculation and error assessment, is shown. Detailed Implementation

[0019] To more clearly illustrate the present invention, the following description, in conjunction with embodiments and accompanying drawings, further explains the invention. Similar components in the drawings are indicated by the same reference numerals. Those skilled in the art should understand that the specific description below is illustrative rather than restrictive and should not be construed as limiting the scope of protection of the present invention.

[0020] Currently, with the improvement and iteration of global satellite navigation systems, countries in various fields such as national security, power, network communication, and navigation are gradually becoming reliant on the various services provided by satellite navigation systems. The improvement in service performance has a positive feedback effect on scenarios requiring high-precision timekeeping. In scenarios requiring remote time comparison and tracing, satellite common-view comparison technology, as a mature high-precision time and frequency transmission technology, works by synchronously observing the same navigation satellite from two remote observation stations. By eliminating common errors such as satellite clock bias and orbital errors, high-precision comparison of the clocks at the two locations is achieved. In traditional satellite common-view comparison processing, the accuracy of satellite orbital position calculation has a significant impact on the comparison of clocks at two locations. Typically, after the receiver obtains the broadcast ephemeris information broadcast by the satellite, it uses a 16-parameter broadcast ephemeris model based on Kepler elements for calculation. The Galileo system and the BDS system both use similar models. This model has a user ranging error of approximately 5 cm for fitting the orbit of a medium Earth orbit satellite. Therefore, it is urgent to further improve the accuracy of satellite orbit position calculation, taking into account the long-period term of orbital changes, thereby reducing calculation errors and improving the accuracy of clock comparison between two locations, thus achieving higher precision.

[0021] In view of this, such as Figure 1 As shown, one embodiment of the present invention provides a method for satellite orbital position calculation and error assessment. The method includes: acquiring satellite ephemeris data, extracting ephemeris parameters from the ephemeris data, and obtaining the satellite's three-dimensional coordinates in the BeiDou coordinate system based on the ephemeris parameters; converting the satellite's three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Earth Reference Frame (IEFrame) coordinate system; acquiring the satellite's true three-dimensional coordinates in the IEFrame coordinate system from a ground tracking observation station, calculating the three-dimensional position deviation vector between the satellite's three-dimensional coordinates in the IEFrame coordinate system and the satellite's true three-dimensional coordinates in the IEFrame coordinate system, and obtaining the satellite's radial orbit deviation, tangential orbit deviation, and normal orbit deviation based on the three-dimensional position deviation vector; acquiring the clock error true value from the ground tracking observation station, acquiring the clock error from the ephemeris data, calculating the clock error deviation between the clock error and the clock error true value, and calculating the space signal ranging error based on the clock error deviation, radial orbit deviation, tangential orbit deviation, and normal orbit deviation.

[0022] In a specific example, a method for satellite orbital position calculation and error assessment based on broadcast ephemeris is provided. The satellite common-view comparison system receives broadcast ephemeris information, calculates the satellite orbital position using an 18-parameter broadcast ephemeris model, performs the conversion from the BeiDou coordinate system to the ITRF2020 coordinate system using the Bursa seven-parameter model, and finally uses the true values ​​obtained from the precise ephemeris file to perform comparison calculations. The spatial signal ranging error is obtained through RTN projection. This method can obtain the broadcast ephemeris or clock difference anomaly status of a single satellite, which is of great significance for satellite selection and accuracy improvement in remote common-view comparison.

[0023] In a specific example, a method is provided to calculate satellite orbital positions using a broadcast ephemeris 18-parameter model, and to transform the BeiDou coordinate system to the ITRF2020 coordinate system using a Bursa 7-parameter model to obtain the satellite position errors in the broadcast ephemeris and precise ephemeris. Projecting the coordinates to the RTN coordinate system and combining the clock bias obtained from the satellite common-view receiver, the ranging errors of various satellite spatial signals can be calculated separately. This error is an important component in evaluating the positioning accuracy of the Global Navigation Satellite System, reflecting the status of individual satellites in a timely manner and providing a reliable data source for remote satellite common-view comparisons.

[0024] The technical solution of this invention provides a method for satellite orbit position calculation and error assessment based on the broadcast ephemeris 18 parameter model. This method can further improve the accuracy of satellite orbit position calculation by taking into account the long-period term of orbit change, thereby reducing calculation error and improving the accuracy of clock comparison between two locations, thus achieving higher accuracy, solving the calculation error caused by orbit change, and simultaneously assessing the error of the calculation results.

[0025] In one possible implementation, the ephemeris parameters include: ephemeris reference time, orbital semi-major axis, orbital eccentricity, orbital inclination, right ascension of the ascending node, perigee angular distance, mean perigee angle at the ephemeris reference time, mean angular velocity perturbation, rate of change of orbital inclination, deviation of rate of change of right ascension of the ascending node, sine harmonic term of latitude argument, cosine harmonic term of latitude argument, sine harmonic term of orbital radius, cosine harmonic term of orbital radius, sine harmonic term of orbital inclination, cosine harmonic term of orbital inclination, rate of change of orbital semi-major axis, and second-order rate of change of right ascension of the ascending node.

[0026] In a specific example, the steps of a satellite orbital position calculation and error assessment method based on broadcast ephemeris include, firstly, establishing a satellite common-view comparison system. A satellite common-view comparison system for remote time synchronization includes a rubidium clock, a satellite common-view receiver, a common-view calculation module, and a computational control core module. The rubidium clock serves as the frequency source, providing a stable frequency signal for the system; the satellite common-view receiver receives broadcast ephemeris information broadcast in real time by the satellite; the computational control core module calculates the acquired satellite data and uses the data to calculate the satellite orbital position and the time difference between the satellite and BeiDou; the common-view calculation module performs remote satellite common-view comparison and conducts time tracing between the two locations.

[0027] In a specific example, such as Figure 2 As shown, the system includes external interfaces, such as an antenna—an external satellite receiving antenna responsible for capturing radio frequency signals from navigation satellites such as BeiDou and GPS, and inputting the signals to the common-view receiver. Other interfaces include a data interface, the device's external communication port, used to output common-view comparison results and time synchronization data, and can also receive external configuration commands.

[0028] Furthermore, the internal core modules include: Rubidium Clock 1, a rubidium atomic clock, which serves as the device's local high-stability frequency reference. It outputs a 10MHz standard frequency signal to provide a reference clock for the common-view receiver; simultaneously, it acts as a local time source, receiving disciplined control from the computational control core module to continuously align with the standard time. Common-view receiver 2 is the satellite navigation receiving unit. Using the 10MHz signal from the rubidium clock as a frequency reference, it acquires, tracks, and demodulates the satellite signals received by the antenna, calculating raw satellite data such as satellite ephemeris, pseudorange observations, and satellite time, and outputting it to the data processing module. Data processing module 4 is the common-view data calculation unit. It receives the raw satellite data output from the receiver, performs data parsing, error correction, and common-view clock bias calculation to obtain the deviation between the local time and the satellite system time; it is also responsible for outputting the final comparison results externally through the data interface. The computational control core module 3 is the device's control and closed-loop correction center. It receives the clock difference results output by the data processing module, calculates the frequency and phase deviations of the local rubidium clock, generates a discipline control signal to feed back to the rubidium clock, and adjusts the output parameters of the rubidium clock in real time, forming a closed-loop control of observation, calculation, and correction. Power supply 5 is the equipment's power supply unit, providing a stable operating power supply to all internal modules.

[0029] Furthermore, a highly stable 10MHz standard frequency signal generated by a rubidium clock is supplied as the reference frequency, providing a unified time and frequency reference for the common-view receiver and ensuring the accuracy of satellite signal processing. Satellite signal reception and processing: Under the 10MHz reference clock, the common-view receiver receives navigation satellite signals through its antenna, demodulates raw observation data such as satellite ephemeris and pseudorange, and uses this as the raw input for common-view comparison. Common-view clock difference resolution: The data processing module processes the satellite observation data, combines parameters such as satellite position and atmospheric correction, and calculates the difference between the local rubidium clock time and the satellite system time (local clock difference). Closed-loop discipline correction: Based on the clock difference result, the arithmetic control core module calculates the frequency and phase adjustments required by the rubidium clock, outputs a discipline control signal to correct the rubidium clock output, and continuously converges the local time to the satellite standard time, enabling the rubidium clock to simultaneously possess short-term high stability and long-term time traceability capabilities. External output: The data processing module outputs the final time synchronization results and common-view comparison data through a data interface for use by external systems.

[0030] In a specific example, such as Figure 3As shown, this method includes inputting a data source to obtain satellite broadcast ephemeris information, i.e., ephemeris parameters (Kepler orbit, perturbation correction parameters) publicly broadcast by navigation satellites to users. This is the data source used by user terminals to actually calculate satellite positions and is also the object of this evaluation. Next, it obtains precise satellite ephemeris information, high-precision satellite orbit products (such as IGS precise ephemeris) obtained through joint post-calculation by global ground monitoring stations. The orbital accuracy can reach the centimeter level and is used as the true value benchmark in this process to measure the magnitude of the broadcast ephemeris error. This data plays two roles in the process: first, it participates in the Bursa model solution to align with the coordinate benchmark; second, it is used in the RTN projection stage to calculate the difference between the true value and the broadcast ephemeris.

[0031] Furthermore, the method includes calculating the satellite's orbital position. Based on the Kepler parameters and perturbation parameters of the broadcast ephemeris, and following the standard ephemeris calculation process (mean anomaly, off-anomaly, true anomaly, orbital plane coordinates, and geocentric-geostatic coordinates), the three-dimensional spatial position of the satellite at each signal transmission time is calculated to obtain the orbital result of the broadcast ephemeris.

[0032] Furthermore, this method involves substituting the broadcast ephemeris position and the precise ephemeris position at the same time as common points into the Bursa seven-parameter model to solve for coordinate transformation parameters. Its core function is to eliminate potential coordinate reference differences between the two ephemeris sets (such as systematic deviations between different Earth reference frames), unify the coordinate reference, and prevent deviations in the coordinate system itself from being miscalculated as orbital errors.

[0033] Furthermore, the method includes coordinate system transformation. Using the obtained Bursa seven parameters, the satellite coordinates of the broadcast ephemeris are completely transformed to the coordinate reference of the precise ephemeris, so that the two sets of orbital data are in the same coordinate system, providing a basis for error comparison.

[0034] Furthermore, this method includes RTN coordinate system projection, which projects the broadcast ephemeris position after standardization and the precise ephemeris true position onto the satellite orbital local coordinate system (RTN coordinate system) to obtain position error components in three directions: radial (R), tangential (T), and normal (N). The RTN coordinate system is defined as follows: radial along the satellite-Geocentric line, tangential along the orbital motion direction, and normal perpendicular to the orbital plane. SISRE is an error index for user-oriented ranging; the weight of position errors in different directions on ranging varies greatly, and it is necessary to decompose it into local orbital directions for accurate calculation.

[0035] Furthermore, the method includes SISRE error calculation, which calculates the spatial signal ranging error value at that moment by weighting the orbital error components in the three directions of RTN and the satellite clock error according to the SISRE standard formula, and thus completing the quantitative assessment of the ranging accuracy of the broadcast ephemeris.

[0036] This embodiment clarifies the complete set of ephemeris parameters required for orbit calculation, including basic Kepler orbit parameters, as well as long-period correction parameters such as mean angular velocity perturbation and orbit inclination rate of change, six short-period harmonic perturbation parameters, and higher-order extrapolation parameters such as semi-major axis rate of change and second-order rate of change of right ascension of ascending node. The parameter dimensions are complete and fully compatible with the parameter system of BeiDou broadcast ephemeris, providing a complete and standardized input basis for high-precision satellite orbit calculation, while supporting the maintenance of orbit accuracy in long-term extrapolation scenarios.

[0037] Furthermore, the method includes a second step: calculating the satellite's orbital position using broadcast ephemeris data. From the broadcast ephemeris data transmitted by the satellite, the satellite's operational status information, including a reference time, can be extracted. Six Kepler orbital parameters, including the semi-major axis deviation. Orbital eccentricity Track inclination Right ascension of ascending node Angular distance from perigee The mean anterior angle at the moment of signal transmission Three linear correction parameters, including the average angular velocity perturbation. Rate of change of orbital inclination Deviation value of the rate of change of right ascension at the ascending node There are 6 sine and cosine harmonic correction parameters, including the cosine harmonic term for latitude argument. Sine harmonic term of latitude argument The cosine harmonic term of the orbital radius The sinusoidal harmonic term of the orbital radius Cosine harmonic term of orbital inclination Sine harmonic term of orbital inclination angle ; and the newly added orbital semi-major axis variation rate based on the 16-parameter model. and mean angular velocity perturbation rate of change .

[0038] Furthermore, using the aforementioned 18 parameters, the planning time is calculated sequentially. The average angular velocity of the satellite The mean anterior angle at the moment of signal transmission Near the point angle True close angle angular distance of ascending node Perturbation correction item ( The perturbation correction amount for the ascending node angular distance at time k. The satellite radius vector length perturbation correction amount at time k, (This is the orbital inclination perturbation correction at time k), and the ascending node angular distance after perturbation correction. Satellite radius length and orbital inclination The position of the satellite in the orbital plane , The first coordinate of the satellite's position in the orbital plane. This is the second coordinate of the satellite's position in the orbital plane. Calculate the corrected right ascension of the ascending node. Then, the three-dimensional coordinates of MEO and IGSO satellites in the BeiDou coordinate system can be calculated. , This is the first coordinate of the three-dimensional coordinates of MEO and IGSO satellites in the BeiDou coordinate system. This is the second coordinate of the three-dimensional coordinates of MEO and IGSO satellites in the BeiDou coordinate system. This is the third coordinate of the three-dimensional coordinates of MEO and IGSO satellites in the BeiDou coordinate system.

[0039] In one possible implementation, obtaining the satellite's three-dimensional coordinates in the BeiDou coordinate system based on ephemeris parameters includes: calculating the time difference between the signal transmission time and the ephemeris reference time; calculating the orbital semi-major axis at the current time using the orbital semi-major axis, time difference, and rate of change of the orbital semi-major axis; calculating the average angular velocity at the current time using the orbital semi-major axis and the average angular velocity perturbation; calculating the mean anomaly angle using the mean anomaly angle at the ephemeris reference time, time difference, and the average angular velocity at the current time; calculating the off-anomaly angle using the mean anomaly angle; calculating the true anomaly angle using the off-anomaly angle; calculating the ascending node distance using the perigee distance and the true anomaly angle; calculating the correction amount for the ascending node distance using the sine and cosine harmonic terms of the latitude argument; calculating the correction amount for the satellite's radius vector length using the sine and cosine harmonic terms of the orbital radius; and calculating the orbital inclination using the sine and cosine harmonic terms of the orbital inclination. The inclination correction is calculated using the ascending node angular distance and its correction value. The satellite's radius reference value is calculated using the current orbital semi-major axis, orbital eccentricity, and perihelion angle. The corrected satellite radius length is calculated using the satellite radius reference value and its correction value. The orbital inclination reference value is calculated using the orbital inclination angle, the rate of change of orbital inclination, and the time difference. The corrected orbital inclination is calculated using the orbital inclination reference value and its correction value. The satellite's polar coordinates are obtained using the corrected ascending node angular distance and the corrected satellite radius length. These polar coordinates are then converted to two-dimensional coordinates in the orbital plane. The corrected ascending node right ascension is calculated using the ascending node right ascension deviation, the second-order rate of change of ascending node right ascension, and the time difference. The corrected ascending node right ascension is obtained by rotating the two-dimensional coordinates around the x-axis to correct the orbital inclination, and then rotating the corrected ascending node right ascension around the z-axis. This yields the satellite's three-dimensional coordinates in the BeiDou coordinate system.

[0040] This embodiment employs a step-by-step calculation logic of time extrapolation, angle calculation, perturbation correction, and coordinate transformation. The process is rigorous and standardized, conforming to the calculation rules of satellite orbital mechanics. By performing independent short-period perturbation corrections on the ascending node angular distance, satellite radius vector, and orbital inclination, the periodic deviations in orbit caused by perturbations such as the Earth's non-spherical gravity and the gravitational forces of the Sun and Moon can be effectively offset. At the same time, it completes a precise rotational transformation from two-dimensional coordinates in the orbital plane to three-dimensional coordinates at the Earth's core and solid surface, ensuring the accuracy of the satellite's three-dimensional position calculation in the BeiDou coordinate system.

[0041] In one possible implementation, the calculation of the orbital semi-major axis at the current moment using the orbital semi-major axis, time difference, and orbital semi-major axis change rate includes: calculating the product between the time difference and the orbital semi-major axis change rate to obtain a first calculation result; and calculating the sum between the orbital semi-major axis and the first calculation result to obtain the orbital semi-major axis at the current moment.

[0042] This embodiment uses a linear extrapolation method to calculate the real-time orbital semi-major axis, which has low computational complexity and high operational efficiency. At the same time, it corrects the long-term drift characteristics of the orbital dimensions by using the semi-major axis change rate, which makes up for the error accumulation problem when extrapolating with a fixed semi-major axis, and can effectively improve the accuracy of orbital dimension calculation for periods other than the ephemeris reference time.

[0043] In one possible implementation, the calculation of the corrected right ascension of the ascending node using the right ascension of the ascending node, the deviation of the rate of change of the right ascension of the ascending node, the second rate of change of the right ascension of the ascending node, and the time difference includes: calculating the product between half of the second rate of change of the right ascension of the ascending node and the square of the time difference to obtain a second calculation result; calculating the product between the deviation of the rate of change of the right ascension of the ascending node and the time difference to obtain a third calculation result; and calculating the sum of the right ascension of the ascending node, the second calculation result, and the third calculation result to obtain the corrected right ascension of the ascending node.

[0044] This embodiment introduces a quadratic term correction of the second-order rate of change based on the first-order linear extrapolation of the right ascension of the ascending node. This can more accurately fit the nonlinear drift law of the right ascension of the ascending node, reduce the calculation error of the azimuth of the ascending node in long-term extrapolation scenarios, improve the accuracy of the solution of the spatial orientation of the orbital plane, and thus ensure the overall accuracy of the satellite's three-dimensional coordinates.

[0045] In one possible implementation, the transformation of the satellite's three-dimensional coordinates in the BeiDou coordinate system to its three-dimensional coordinates in the International Earth Reference Frame (IEFrame) coordinate system includes: selecting at least three locations with their respective three-dimensional coordinates in the BeiDou coordinate system and the IEFrame coordinate system; calculating the optimal set of parameters in the Bursa seven-parameter coordinate transformation model using the least squares method based on the three-dimensional coordinates of the at least three locations in the BeiDou coordinate system and the IEFrame coordinate system, thereby obtaining the final Bursa seven-parameter coordinate transformation model; and using the final Bursa seven-parameter coordinate transformation model to transform the satellite's three-dimensional coordinates in the BeiDou coordinate system to its three-dimensional coordinates in the IEFrame coordinate system.

[0046] Furthermore, the method includes a third step: coordinate transformation using the Bursa seven-parameter model. The definition of the BeiDou coordinate system (BDCS) conforms to the International Earth Rotation Service (IERS) standard and is consistent with the definition of the China Geodetic Coordinate System 2000 (CGCS2000), having identical reference ellipsoid parameters. The geometric center of the reference ellipsoid coincides with the Earth's center of mass, and the rotation axis of the reference ellipsoid coincides with the Z-axis. To facilitate comparison with precise ephemeris files, the acquired satellite orbital positions need to be transformed from the BeiDou coordinate system to the ITRF2020 coordinate system. This invention uses the Bursa seven-parameter model for this transformation. Based on error theory and the fundamentals of measurement adjustment, the coordinate transformation model between the two spatial rectangular coordinate systems is expressed as follows:

[0047] In the formula, These are the coordinates in the target space Cartesian coordinate system (ITRF), and These are the coordinates in the source space Cartesian coordinate system (BDCS), where ΔX, ΔY, and ΔZ are the first, second, and third translation parameters for transforming from the source space Cartesian coordinate system to the target space Cartesian coordinate system. , , These are the first rotation parameters, second rotation parameters, and third rotation parameters for transforming the source space Cartesian coordinate system to the target space Cartesian coordinate system. m The scale parameter is used to transform the source space rectangular coordinate system to the target space rectangular coordinate system; To rotate the coordinate system about the Z-axis The resulting rotation matrix; To rotate the coordinate system about the Y-axis The resulting rotation matrix; To rotate the coordinate system about the X-axis The resulting rotation matrix.

[0048] Furthermore, the calculation formulas for the three rotation matrices are as follows:

[0049]

[0050]

[0051] Furthermore, after processing, the Bursa seven-parameter coordinate transformation model can be obtained:

[0052] In the formula, Let be the first coordinate in the target space rectangular coordinate system (ITRF). The second coordinate in the target space rectangular coordinate system (ITRF) is... The third coordinate in the target space rectangular coordinate system (ITRF) is... Let be the first coordinate in the source space rectangular coordinate system (BDCS). The second coordinate in the source space rectangular coordinate system (BDCS) is... As shown above, the Bursa transformation model has a total of 7 transformation parameters: 3 translation parameters (ΔX, ΔY, ΔZ) and 3 rotation parameters (ΔX, ΔY, ΔZ). This represents the third coordinate in the source space rectangular coordinate system (BDCS). , , ) and 1 scale parameter m.

[0053] Furthermore, when performing coordinate transformation, an indirect adjustment model needs to be constructed based on the above formula:

[0054] In the formula, The first component of the coordinate transformation residual of the first common point. The second component of the coordinate transformation residual of the first common point. The third component of the coordinate transformation residual of the first common point. The first component of the coordinate transformation residual of the Nth common point. The second component of the coordinate transformation residual of the Nth common point. The third component of the coordinate transformation residual of the Nth common point. Let the first coordinate of the first common point be... The second coordinate of the first common point. The third coordinate of the first common point Let the first coordinate of the Nth common point be... Let the coordinates of the Nth common point be two coordinates. Let be the third coordinate of the Nth common point. Solving this model requires 3 or more common points. If we assume there are N common points, then the total number of observations n is 3N, the necessary number of observations is 7, and the degrees of freedom r is 3N-7. In this case, the least squares method can be used directly to solve the indirect adjustment model.

[0055] Furthermore, the above formula can be simplified to...

[0056] In the formula, For coordinate transformation residual vector, This is the coefficient matrix of the Bursa model. The optimal estimates for the seven parameters are... This is the coordinate difference observation vector.

[0057] Furthermore, after simplification, the parameter estimates are obtained as follows: Substituting the seven solved parameters into the Bursa model yields the error of each coordinate, and the accuracy is evaluated according to the following formula.

[0058]

[0059] In the formula, The error is the unit weight standard error after verification. The observation weight matrix, The total number of observations. This is the number of necessary observations.

[0060] This embodiment uses the least squares method to solve for the optimal Bursa seven parameters by using the common point coordinates. This can globally and optimally eliminate systematic deviations such as translation, rotation, and scaling between the BeiDou coordinate system and the International Earth Reference Frame, achieving high-precision alignment of the two coordinate references. It can also avoid the systematic errors of the coordinate system itself being mistakenly included in the orbital position error, ensuring the objectivity and reliability of subsequent orbital deviation and ranging error assessment results.

[0061] In one possible implementation, obtaining the satellite's radial orbital deviation, tangential orbital deviation, and normal orbital deviation based on the three-dimensional position deviation vector includes: establishing the satellite's orbital coordinate system, which includes radial, tangential, and normal axes; and projecting the three-dimensional position deviation vector onto the radial, tangential, and normal axes of the orbital coordinate system to obtain the radial orbital deviation, tangential orbital deviation, and normal orbital deviation.

[0062] This embodiment projects the three-dimensional position deviation vector onto the radial, tangential, and normal directions of the local coordinate system of the track, thereby achieving directional decoupling of track errors. It can distinguish the differentiated impact of position errors in different directions on user ranging, providing a decomposition basis for the accurate weighted calculation of subsequent spatial signal ranging errors and avoiding the distortion of results caused by directly using the three-dimensional total deviation assessment.

[0063] In one possible implementation, the calculation of spatial signal ranging error based on clock bias, radial orbital bias, tangential orbital bias, and normal orbital bias includes: calculating the distance error by multiplying the clock bias and the speed of light; and calculating the spatial signal ranging error based on the distance error, radial orbital bias, tangential orbital bias, and normal orbital bias.

[0064] In a specific example, the method includes a fourth step: calculating the space signal ranging error. User equivalent range error (UERE) is a component of the positioning accuracy of a Global Navigation Satellite System (GNSS), while space signal ranging error (SISRE) is a major component of UERE, reflecting the impact of deviations in the navigation messages broadcast by the satellites on user ranging.

[0065] Furthermore, based on the aforementioned Bursa model, the satellite orbital position is transformed from the BeiDou coordinate system to the ITRF2020 coordinate system, and then the position deviation vector is calculated. As shown in the formula below:

[0066] In the formula, Calculate the position vector for broadcast ephemeris. This is the satellite position vector in the precise ephemeris.

[0067] Furthermore, the basis vectors of the RTN coordinate system, including the radial vector, are calculated. Normal vector and tangential vector The formula is as follows:

[0068]

[0069]

[0070] In the formula, This is a precise ephemeris satellite three-dimensional position vector; This is the satellite velocity vector, which can be obtained from a precise ephemeris file.

[0071] Furthermore, the RTN deviation is obtained by projecting the basis vectors of the RTN coordinate system.

[0072]

[0073] In the formula, For radial orbital deviation, For tangential trajectory deviation, This represents the deviation in the normal orbit.

[0074] This embodiment converts clock bias into distance dimension by using the speed of light, and after unifying the dimensions with orbital position bias, they are used together in the ranging error calculation, with clear and consistent physical meaning; at the same time, it fully covers the two core sources of space signal ranging error: satellite orbital position error and satellite clock bias error, ensuring the completeness of error assessment.

[0075] In one possible implementation, the calculation of the spatial signal ranging error based on the distance error, radial orbit deviation, tangential orbit deviation, and normal orbit deviation includes: calculating the product between the first weight and the radial orbit deviation to obtain a fourth calculation result; calculating the difference between the fourth calculation result and the distance error to obtain a fifth calculation result; calculating the sum between the square of the tangential orbit deviation and the square of the normal orbit deviation to obtain a sixth calculation result; dividing the sixth calculation result by the second weight to obtain a seventh calculation result; calculating the sum between the fifth and seventh calculation results to obtain an eighth calculation result; and calculating the square root of the eighth calculation result to obtain the spatial signal ranging error.

[0076] In a specific example SISRE For all points within the satellite's Earth surface coverage area, instantaneous SISRE The RMS statistical value can be calculated using the satellite's radial, tangential, and normal orbital errors and clock biases according to the following formula:

[0077] In the formula, SISRE is the spatial signal ranging error; c is the speed of light. This is clock bias. For radial orbital deviation, For tangential trajectory deviation, This represents the normal orbital deviation. The values ​​of the first parameter α and the second parameter β are shown in Table 1 below.

[0078] Table 1. Values ​​of the first parameter α and the second parameter β

[0079] In a specific example, the algorithm flow of a satellite orbital position calculation and error assessment method based on broadcast ephemeris includes the common-view receiver of the satellite common-view comparison system sending the broadcast ephemeris information broadcast by the satellite into the data processing module to calculate the satellite's orbital position. Within the module, the Bursa seven-parameter model solution, coordinate system transformation, RTN coordinate system projection, and SISRE error calculation are performed sequentially. The precise ephemeris file information used in the calculation can be post-processed through the data interface. After obtaining the broadcast ephemeris information, 18 state parameters of the satellite's operation are extracted.

[0080] Furthermore, the method includes calculating the planning time. .

[0081]

[0082] In the formula, The time of signal transmission. For ephemeris reference time, for a valid ephemeris, The value should be in Between the two hours before and after, that is The absolute value must be less than 7200s.

[0083] Furthermore, the method includes calculating the major semi-axis. .

[0084]

[0085]

[0086] In the formula, It is the semi-major axis of the reference time. This is the reference value for the semi-major axis. The reference value for MEO satellites is 27,906,100 m, and the reference value for IGSO satellites is 42,162,200 m. This is the deviation value of the semi-major axis of the track; This represents the rate of change of the semi-major axis of the orbit.

[0087] Furthermore, the method includes calculating the satellite's average angular velocity. According to Kepler's third law, the average angular velocity of a satellite can be derived from the following formula:

[0088] in, GM The gravitational constant of Earth; constant =3.986005×10¹⁴ ; The corrected average angular velocity of a hypothetical satellite orbiting in a circular orbit. for:

[0089]

[0090] In the formula, For the first k The average angular velocity perturbation correction at time t. The average angular velocity perturbation at the reference time. This represents the rate of change of the average angular velocity perturbation.

[0091] Furthermore, the method includes calculating the mean anterior angle at the moment of signal transmission. .

[0092]

[0093] In the formula, The value should be between 0 and between, The value is 3.1415926535898. The reference angle is the mean aperitoneal angle.

[0094] Furthermore, the proximity angle at the time of signal transmission is calculated. The relationship between the anomalous and plan anomalies can be given by Kepler's equations:

[0095] In the formula, This is the estimated near-point angle for the j-th iteration. The angle of approach at the current moment, For satellite orbital eccentricity, It represents the near-point angle value in the (j-1)th iteration.

[0096] Furthermore, an iterative method is used to adjust the near-point angle. Solve from the Kepler equations, let initial value of iteration Set as ,Right now , , The result of two iterations can be regarded as the near-point angle value, where, Let be the angle of approach at time k. The near-point angle obtained in the first iteration. The near-point angle is obtained from the second iteration.

[0097] Furthermore, the method includes calculating the true anterior angle at time k. .

[0098]

[0099]

[0100] Furthermore, the true anterior angle value can be obtained from the sine and cosine values:

[0101] In the formula, the true anterior angle value ranges from .

[0102] Furthermore, the angular distance of the rising node at the moment of signal transmission is calculated. .

[0103]

[0104] In the formula, Let be the true perimeter angle at time k. This is the argument of perigee.

[0105] Furthermore, the method includes calculating the perturbation correction term at the signal transmission time. , , .

[0106]

[0107]

[0108]

[0109] In the formula, The perturbation correction amount for the angular distance of the ascending intersection at time k. The coefficient of the sine term of the perturbation of the angular distance of the ascending node. The coefficient of the cosine term of the perturbation at the angular distance of the ascending intersection is... Let $k$ be the satellite radius perturbation correction value at time $k$. The coefficients of the sinusoidal term of the radius perturbation are... The coefficients of the cosine term of the radius perturbation are... This is the correction amount for the orbital inclination perturbation at time k. The coefficients of the sinusoidal term for the orbital inclination perturbation are... The coefficient of the cosine term of the orbital inclination perturbation.

[0110] Furthermore, the method includes calculating the perturbation-corrected ascending node angular distance. Satellite radius length and orbital inclination .

[0111]

[0112]

[0113]

[0114] In the formula, This is the corrected semi-major axis of the orbit at the current moment. For reference time, the orbital inclination angle, This represents the rate of change of the orbital inclination angle.

[0115] Furthermore, the method includes calculating the satellite's position in the orbital plane at the moment of signal transmission. Polar coordinates can be converted using the following formula. Transform into a Cartesian coordinate system on the orbital plane coordinates in ,

[0116]

[0117] In the formula, Let k be the first coordinate of the orbital plane at time k. Let k be the satellite radius vector length after perturbation correction at time k. The angular distance of the ascending intersection point after perturbation correction at time k. Let be the second coordinate of the orbital plane at time k.

[0118] Furthermore, the method includes calculating the right ascension of the ascending node at the moment of signal transmission. .

[0119]

[0120] In the formula, The right ascension of the ascending node is the ephemeris reference time. The rate of change of right ascension at the ascending node. This is the Earth's rotational angular velocity. This is the ephemeris reference time. Note that after the right ascension of the ascending node is determined, it must be equivalent to 0 to 1. between.

[0121] Furthermore, the method includes calculating coordinates in the BeiDou coordinate system. .

[0122]

[0123]

[0124]

[0125] In the formula, Let k be the first coordinate in the BeiDou geocentric coordinate system at time k. The second coordinate in the BeiDou geocentric coordinate system at time k is... The third coordinate in the BeiDou geocentric coordinate system at time k is... The orbital inclination angle after perturbation correction at time k.

[0126] Furthermore, the satellite orbital position calculation and error assessment methods, utilizing the Bursa seven-parameter model and RTN coordinate system projection, are key steps in evaluating satellite orbital position calculations, providing a data foundation for assessing space signal ranging errors. The computational control core module obtains the post-processed precise ephemeris file through a data interface, extracting the required precise satellite velocity vector from it. Satellite position vector This is called the truth value.

[0127] Furthermore, the satellite common-view receiver receives and transmits the data to the data processing module. The satellite orbital position information obtained by the 18-parameter broadcast ephemeris model is compared with the true value. Three or more common points are selected from them, and an indirect adjustment model is constructed. The coordinate system transformation parameters can then be obtained by solving the least squares method.

[0128] Furthermore, the transformed satellite orbital position vector is expanded and vector operations are performed with the true value to obtain the RTN coordinate system basis vectors, including the radial vector. Normal vector and tangential vector The RTN deviation is obtained by projecting the basis vectors of the RTN coordinate system.

[0129] Furthermore, by substituting the satellite's radial, tangential, and normal orbital errors and clock bias into the SISRE error calculation formula, the instantaneous error assessment value of a single satellite can be obtained.

[0130] This embodiment addresses the varying degrees of impact of radial, tangential, and normal errors on user ranging by employing differentiated weighted calculation rules: radial and clock errors directly affect the ranging results and are given high weights; tangential and normal errors are factored in through geometric factors before being included in the total error. This calculation method better reflects the actual geometric characteristics of user ranging, accurately reflects the actual impact of broadcast ephemeris errors on terminal ranging accuracy, and improves the accuracy and engineering application value of spatial signal ranging error assessment results.

[0131] Another embodiment of the present invention provides a satellite orbital position calculation and error assessment system, which includes: a data acquisition module for acquiring satellite ephemeris data, extracting ephemeris parameters from the ephemeris data, and obtaining the satellite's three-dimensional coordinates in the BeiDou coordinate system based on the ephemeris parameters; a coordinate transformation module for transforming the satellite's three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Geodetic Frame (IGFR) coordinate system; an orbital deviation calculation module for acquiring the satellite's true three-dimensional coordinates in the IGFR coordinate system from a ground tracking observation station, calculating the three-dimensional position deviation vector between the satellite's three-dimensional coordinates in the IGFR coordinate system and the satellite's true three-dimensional coordinates in the IGFR coordinate system, and obtaining the satellite's radial orbital deviation, tangential orbital deviation, and normal orbital deviation based on the three-dimensional position deviation vector; and a space signal ranging error calculation module for acquiring the clock error true value from the ground tracking observation station, acquiring the clock error from the ephemeris data, calculating the clock error deviation between the clock error and the clock error true value, and calculating the space signal ranging error based on the clock error deviation, radial orbital deviation, tangential orbital deviation, and normal orbital deviation.

[0132] The technical solution of this invention provides a method for satellite orbit position calculation and error assessment based on the broadcast ephemeris 18 parameter model. This method can further improve the accuracy of satellite orbit position calculation by taking into account the long-period term of orbit change, thereby reducing calculation error and improving the accuracy of clock comparison between two locations, thus achieving higher accuracy, solving the calculation error caused by orbit change, and simultaneously assessing the error of the calculation results.

[0133] In the description of this invention, it should be noted that the terms "upper," "lower," etc., indicating the orientation or positional relationship are based on the orientation or positional relationship shown in the accompanying drawings, and are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Unless otherwise expressly specified and limited, the terms "installed," "connected," and "linked" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; they can refer to the internal communication between two elements. For those skilled in the art, the specific meaning of the above terms in this invention can be understood according to the specific circumstances.

[0134] It should also be noted that in the description of this invention, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0135] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the implementation of the present invention. For those skilled in the art, other variations or modifications can be made based on the above description. It is impossible to exhaustively list all the implementation methods here. All obvious variations or modifications derived from the technical solutions of the present invention are still within the protection scope of the present invention.

Claims

1. A method for satellite orbital position calculation and error assessment, characterized in that, The method includes: Obtain satellite ephemeris data, extract ephemeris parameters from the ephemeris data, and obtain the satellite's three-dimensional coordinates in the BeiDou coordinate system based on the ephemeris parameters; Transform the satellite's three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Earth Reference Frame coordinate system; The true three-dimensional coordinates of the satellite in the coordinate system of the International Earth Reference Frame are obtained from the ground tracking observation station. The three-dimensional position deviation vector between the three-dimensional coordinates of the satellite in the coordinate system of the International Earth Reference Frame and the true three-dimensional coordinates of the satellite in the coordinate system of the International Earth Reference Frame is calculated. Based on the three-dimensional position deviation vector, the radial orbit deviation, tangential orbit deviation, and normal orbit deviation of the satellite are obtained. The clock error is obtained from the ground tracking observation station, the clock error is obtained from the ephemeris data, the clock error deviation between the clock error and the clock error true value is calculated, and the space signal ranging error is calculated based on the clock error deviation, radial orbit deviation, tangential orbit deviation, and normal orbit deviation.

2. The satellite orbital position calculation and error assessment method according to claim 1, characterized in that, The ephemeris parameters include: ephemeris reference time, orbital semi-major axis, orbital eccentricity, orbital inclination, right ascension of the ascending node, perigee angular distance, mean perigee angle at the ephemeris reference time, mean angular velocity perturbation, rate of change of orbital inclination, deviation of rate of change of right ascension of the ascending node, sine harmonic term of latitude argument, cosine harmonic term of latitude argument, sine harmonic term of orbital radius, cosine harmonic term of orbital radius, sine harmonic term of orbital inclination, cosine harmonic term of orbital inclination, rate of change of orbital semi-major axis, and second-order rate of change of right ascension of the ascending node.

3. The satellite orbital position calculation and error assessment method according to claim 1, characterized in that, The process of obtaining the satellite's three-dimensional coordinates in the BeiDou coordinate system based on ephemeris parameters includes: Calculate the time difference between the signal transmission time and the ephemeris reference time; The orbital semi-major axis at the current moment is calculated using the orbital semi-major axis, time difference, and rate of change of orbital semi-major axis. The average angular velocity at the current moment is calculated using the orbital semi-major axis at the current moment and the average angular velocity perturbation. The mean anomaly is calculated using the mean anomaly at the ephemeris reference time, the time difference, and the average angular velocity at the current time. The deviated anomaly is calculated using the mean anomaly, and the true anomaly is calculated using the deviated anomaly. Calculate the angular distance to the ascending node using the angular distance to the perigee and the true anomaly. The correction amount for the ascending node distance is calculated using the sine and cosine harmonic terms of the latitude argument; the correction amount for the satellite radius length is calculated using the sine and cosine harmonic terms of the orbital radius; the correction amount for the orbital inclination is calculated using the sine and cosine harmonic terms of the orbital inclination; the corrected ascending node distance is calculated using the ascending node distance and its correction amount; the reference value for the satellite radius is calculated using the current orbital semi-major axis, orbital eccentricity, and perihelion angle; the corrected satellite radius length is calculated using the reference value for the satellite radius and its correction amount; the reference value for the orbital inclination is calculated using the orbital inclination, the rate of change of orbital inclination, and the time difference; and the corrected orbital inclination is calculated using the reference value for the orbital inclination and its correction amount. The polar coordinates of the satellite are obtained using the corrected ascending node angular distance and the corrected satellite radius vector length, and then the polar coordinates of the satellite are converted into two-dimensional coordinates of the satellite in the orbital plane. The corrected right ascension of the ascending node is calculated using the right ascension of the ascending node, the deviation of the rate of change of the right ascension of the ascending node, the second-order rate of change of the right ascension of the ascending node, and the time difference. The orbital inclination is corrected by rotating the two-dimensional coordinates around the x-axis and then the corrected right ascension of the ascending node is rotated around the z-axis to obtain the three-dimensional coordinates of the satellite in the BeiDou coordinate system.

4. The satellite orbital position calculation and error assessment method according to claim 3, characterized in that, The calculation of the orbital semi-major axis at the current moment using the orbital semi-major axis, time difference, and rate of change of orbital semi-major axis includes: The first calculation result is obtained by multiplying the time difference by the rate of change of the orbit's semi-major axis. The semi-major axis of the orbit is obtained by summing the result of the first calculation with the semi-major axis of the orbit at the current moment.

5. The satellite orbital position calculation and error assessment method according to claim 3, characterized in that, The corrected right ascension of the ascending node, calculated using the right ascension of the ascending node, the deviation of the rate of change of the right ascension of the ascending node, the second-order rate of change of the right ascension of the ascending node, and the time difference, includes: The second calculation result is obtained by multiplying the second rate of change of right ascension of half the ascending node by the square of the time difference; The third calculation result is obtained by multiplying the deviation of the rate of change of right ascension of the ascending node with the time difference; The corrected right ascension of the ascending node is obtained by summing the second and third calculation results.

6. The satellite orbital position calculation and error assessment method according to claim 1, characterized in that, The process of transforming the satellite's three-dimensional coordinates in the BeiDou coordinate system to its three-dimensional coordinates in the International Geodetic Reference Frame includes: Select at least three locations, respectively, in the three-dimensional coordinates of the BeiDou coordinate system and the three-dimensional coordinates of the International Earth Reference Frame coordinate system; Based on the three-dimensional coordinates of at least three locations in the BeiDou coordinate system and the coordinate system of the International Earth Reference Frame, the optimal set of parameters in the Bursa seven-parameter coordinate transformation model is calculated using the least squares method, and the final Bursa seven-parameter coordinate transformation model is obtained. The final Bursa seven-parameter coordinate transformation model is used to transform the satellite's three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Earth Reference Frame coordinate system.

7. The satellite orbital position calculation and error assessment method according to claim 1, characterized in that, The process of obtaining the satellite's radial orbital deviation, tangential orbital deviation, and normal orbital deviation based on the three-dimensional position deviation vector includes: Establish the satellite's orbital coordinate system, which includes radial, tangential, and normal coordinates; The radial orbital deviation, tangential orbital deviation, and normal orbital deviation are obtained by projecting the three-dimensional position deviation vector onto the radial, tangential, and normal coordinates of the orbital coordinate system.

8. The satellite orbital position calculation and error assessment method according to claim 1, characterized in that, The calculation of spatial signal ranging error based on clock bias, radial orbital bias, tangential orbital bias, and normal orbital bias includes: The distance error is obtained by calculating the product of clock error and the speed of light. The spatial signal ranging error is calculated based on the distance error, radial orbit deviation, tangential orbit deviation, and normal orbit deviation.

9. The satellite orbital position calculation and error assessment method according to claim 8, characterized in that, The calculation of spatial signal ranging error based on distance error, radial orbit deviation, tangential orbit deviation, and normal orbit deviation includes: The fourth calculation result is obtained by multiplying the first weight by the radial trajectory deviation. The fifth calculation result is obtained by calculating the difference between the fourth calculation result and the distance error. The sixth calculation result is obtained by summing the squares of the tangential trajectory deviation and the normal trajectory deviation. The sixth calculation result is divided by the second weight to obtain the seventh calculation result; The sum of the fifth and seventh calculation results is used to obtain the eighth calculation result. The spatial signal ranging error is obtained by calculating the square root of the eighth calculation result.

10. A satellite orbital position calculation and error assessment system, characterized in that, The system includes: The data acquisition module is used to acquire satellite ephemeris data, extract ephemeris parameters from the ephemeris data, and obtain the satellite's three-dimensional coordinates in the BeiDou coordinate system based on the ephemeris parameters. The coordinate transformation module is used to transform the satellite's three-dimensional coordinates in the BeiDou coordinate system to the satellite's three-dimensional coordinates in the International Earth Reference Frame coordinate system. The orbital deviation calculation module is used to obtain the satellite's true three-dimensional coordinates in the coordinate system of the International Earth Reference Frame from the ground tracking observation station, calculate the three-dimensional position deviation vector between the satellite's three-dimensional coordinates in the coordinate system of the International Earth Reference Frame and the satellite's true three-dimensional coordinates in the coordinate system of the International Earth Reference Frame, and obtain the satellite's radial orbital deviation, tangential orbital deviation, and normal orbital deviation based on the three-dimensional position deviation vector. The space signal ranging error calculation module is used to obtain the true value of clock error from the ground tracking observation station, obtain the clock error from the ephemeris data, calculate the clock error deviation between the clock error and the true value of clock error, and calculate the space signal ranging error based on the clock error deviation, radial orbit deviation, tangential orbit deviation, and normal orbit deviation.

Images