A low-orbit satellite Doppler frequency offset accurate prediction method, storage medium and system

Abstract

This invention discloses a method, storage medium, and system for accurate prediction of Doppler frequency offset of low-Earth orbit satellites, which runs on a general-purpose computer platform. The method employs the following complementary and synergistic technical elements: (A) orbital propagation based on Kepler orbital elements and analytical Newton iteration; (B) a complete implementation of the IAU1980 precession-nutation model retaining no fewer than five main nutation terms; and (C) dynamic updating of the Greenwich Mean Time (GMT) based on the Julian day at the signal transmission time. 2000 →TOD→ECEF coordinate transformation; (D) an optical time convergence process with the receiving time as the initial value and at least 2 iterations, and the coordinate transformation described in feature (C) is called independently in each optical time iteration; (E) a numerical implementation strategy using 64-bit IEEE 754 double-precision floating-point arithmetic throughout; This invention is applied to the technical field of satellite telemetry and control and ground receiving signal processing.

Description

Technical Field

[0001] This invention relates to satellite telemetry, tracking, and command (TT&C) and ground-based signal processing, and particularly to a method, storage medium, and system for accurate prediction of Doppler frequency offset of low-Earth orbit satellites running on a general-purpose computer platform. Background Technology

[0002] Low-Earth orbit satellites (orbital altitude 200-2000 km) can move at speeds relative to the Earth's center of gravity up to approximately 7.8 km / s. When passing over ground stations, they generate significant azimuth / elevation angular velocities (peak values ​​can reach several degrees per second) and Doppler frequency offsets (the Doppler peak at a Ku-band carrier frequency of approximately 17.8 GHz is on the order of ±400 kHz, far exceeding the typical ±40 kHz acquisition bandwidth of engineering receivers; therefore, pre-acquisition frequency pre-compensation based on predicted Doppler values ​​is necessary). Ground receiving systems need to acquire high-precision azimuth, elevation, and Doppler frequency offset time series before the satellite passes over the ground to drive antenna servo tracking and carrier frequency compensation before signal acquisition.

[0003] Constructing low Earth orbit (LEO) satellite transit prediction algorithms based on Kepler orbital elements is a common technical approach in this field. However, implementing this approach into a stable, reliable, and accurate prediction system that aligns point-by-point with professional simulation software (such as STK) requires careful handling of several technical factors. Improper handling of any one of these factors will significantly reduce the final prediction accuracy: Element 1: Completeness of Precession and Nutation Modeling. The Earth's rotation axis exhibits long-term precession and short-period nutation in inertial space. Nutation includes multiple frequency components such as lunar period, semi-lunar period, and annual period. If only a polynomial approximation is made for the precession while ignoring the nutation, or if only a simplified model of the first few terms of the nutation sequence is taken, the coordinate transformation from the inertial frame to the Earth-fixed frame will accumulate angular errors on the order of sub-arcsecond to arcsecond.

[0004] Element Two: Earth's Rotation Correction During Signal Propagation Time. The propagation time (optical travel time τ) of a signal from the satellite to the ground station is approximately 1.3 ms to 3.3 ms in the LEO scenario, while the Earth's rotational angular velocity ω... E ≈15.041 arcseconds / second, during which time the Earth rotates about 0.020 to 0.050 arcseconds around its polar axis. If the Greenwich Mean Time (GMT) at the time of reception is used consistently in the coordinate transformation... Completing the TOD→ECEF rotation of the satellite position at launch time is equivalent to projecting the "signal transmitted by the satellite at time t" onto the "Earth orientation corresponding to the ground station at time t+τ". The semantic confusion of time introduces a systematic position offset related to orbital geometry. This offset corresponds to an order of approximately 0.67 to 1.67 meters (sub-meter to meter level) at the LEO orbital altitude (geocenter distance of approximately 6900 km). Although the order of magnitude is limited, as a systematic deviation related to orbital geometry, it accumulates with other error sources in the engineering delivery scenario through the chain transmission of radial velocity ρ and Doppler frequency offset Δf, and still constitutes a non-negligible link in accuracy loss.

[0005] Element 3: Iterative convergence of the optical travel time itself. The optical travel time τ depends on the satellite position at the time of transmission, and the time of transmission in turn depends on τ. If τ is directly approximated by "position at the time of reception / c" without iteration, then τ has a zero-order residual that is proportional to the radial velocity of the satellite, which further couples and amplifies the position error and the Doppler error.

[0006] Element Four: High-precision floating-point representation. Julian day values ​​in J... 2000 The approximate location is 2451545, with a time precision in the second range corresponding to approximately 1.16 × 10⁻⁶ Julian days. -5 The absolute order of magnitude, with microsecond precision, corresponds to approximately 1.16 × 10⁻⁶. -11 Single-precision IEEE 754 floating-point numbers have approximately 7 decimal bits, which cannot maintain microsecond precision at this level. Mixing single and double precision can also easily introduce truncation errors when adding or subtracting large numbers.

[0007] Element 5: Correct distinction of time scales. The polynomial parameter of precession nutation is defined in Earth time TT, while the polynomial of sidereal time GMST is defined in UT1 (usually approximated by UTC in engineering). The current cumulative deviation between the two is 69.184 seconds. If UTC is used as the polynomial time parameter for all terms, observable phase deviations will accumulate in the long term.

[0008] In existing known technologies, there are implementation schemes that simplify the processing of the above elements to varying degrees: for example, using only J2 perturbation to approximate the complete nutation, using only the GAST at the receiving time to complete the coordinate rotation, not performing optical travel time iteration, and using single-precision floating point, etc. These simplifications are acceptable in low-precision demand scenarios, but they are difficult to meet the accuracy requirements in engineering delivery scenarios that require point-by-point alignment with professional simulation software and in high-frequency (Ku / Ka) Doppler pre-compensation scenarios.

[0009] The technical problem to be solved by this invention is how to properly handle the above five elements on a general-purpose computer platform and make them work together to form a complete engineering-level accurate forecasting scheme. Summary of the Invention

[0010] The technical problem to be solved by the present invention is to overcome the shortcomings of the prior art and provide a method, storage medium and system for accurate prediction of Doppler frequency offset of low-Earth orbit satellites running on a general computer platform, which mainly solves the engineering-level accuracy problem of prediction of azimuth, elevation angle and Doppler frequency offset during the transit of low-Earth orbit satellites.

[0011] The technical solution adopted in this invention is as follows: This invention includes a method for accurate prediction of Doppler frequency offset of low-Earth orbit satellites, a storage medium, and a system. The method for accurate prediction of Doppler frequency offset of low-Earth orbit satellites includes the following steps: S1. Obtain satellite orbit parameters, orbit epoch times, and ground station coordinate parameters; S2. Determine the Julian day (JD) corresponding to the predicted time based on the orbital epoch time. receive ; S3. Based on the satellite orbit parameters and the Julian day (JD) of the receiving time. receive By solving the Kepler equations, the satellite's position in J... 2000 The position vector r of the inertial coordinate system J2000 and velocity vector v J2000 ; S4a. The received time is Julian day JD. receive Let the initial signal be the initial value of the light travel time. τ 0 =0, initial value of Julian day at launch time =JD receive Perform the optical travel time iteration step; S4b. In each of the aforementioned optical travel time iteration steps, based on the Julian date of the signal transmission time corresponding to the previous iteration step k. Using the IAU1980 precession-nutation model as the time reference, the nutation longitude ΔΨ and nutation angle Δε are calculated by calling no fewer than five main nutation parameters, and based on the Julian date of the signal transmission time. Independent calculation of Greenwich Mean Sidereal Time Greenwich True Sidereal Time According to the nutation parameters and Greenwich Mean Time Position vector r at time J2000 and velocity vector v J2000 After TOD coordinate transformation and then ECEF coordinate transformation, the satellite's ground-fixed position vector r at the launch time is obtained. T and the satellite's ground-to-solid velocity vector v at launch time T ; S4c. Based on the satellite's ground-based position vector r at the launch time. T When updating the optical path with the ground station coordinate parameters, and according to the Julian day JD of the reception time. receive The Julian date for the signal transmission time is iteratively updated according to the optical travel time. Return to the optical travel iteration step until the preset iteration termination condition is met, and the number of iterations is not less than 2. S5. Based on the satellite ground-to-solid position vector r obtained at the launch time when the preset iteration termination condition is met. T Satellite ground-to-solid velocity vector v at launch time T And the radial velocity ρ is determined using the ground station coordinate parameters; S6. Determine the Doppler frequency offset Δf based on the radial velocity ρ; All numerical calculations in steps S1-S6 above are performed using a 64-bit IEEE 754 double-precision floating-point data type to avoid the loss of significant bits in single-precision floating-point numbers near Julian large numbers.

[0012] Furthermore, in step S1, the satellite orbital parameters are the six orbital elements, which include the semi-major axis a, eccentricity e, orbital inclination i, right ascension of the ascending node Ω, argument of perigee ω, and mean perigee M0. The ground station coordinate parameters include latitude φ, longitude λ, and elevation h. In step S2, the Julian day of the receiving time for each forecast time is calculated recursively according to a preset step size, based on the epoch time. receive ; Step S3 involves calculating the mean angle from the current forecast time. Where M is the mean anomaly angle at the current time t; M0 is the mean anomaly angle at epoch t0; and n is the mean angular velocity; additionally... GM is the Earth's gravitational constant. The Kepler equations are solved iteratively. We obtain the near-point angle E, and then obtain the satellite's position at J. 2000 The position vector r of the inertial coordinate system J2000 and velocity vector v J2000 .

[0013] Furthermore, step S4b specifically includes: using the signal transmission time corresponding to the current iteration step k. Using time as the reference, the following transformations are performed sequentially: (i) from

[0014] Calculating Earth Time Julian Day ,in (ii) Using the Julian Day based on Coordinated Universal Time (UTC), TAI is International Atomic Time, UTC is Coordinated Universal Time, TAI-UTC, to obtain the cumulative leap seconds; (iii) Using the IAU1980 precession-nutation model that retains no less than five main nutation parameters, calculate the nutation longitude ΔΨ and nutation angle Δε; (iii) Using the IAU1976 precession matrix P and nutation matrix N, r J2000v J2000 Rotate to the true equator and true vernal equinox coordinate system to obtain r TOD v TOD (iv) by Independent calculation of Greenwich Mean Sidereal Time Greenwich True Sidereal Time (v) with Greenwich true stars Corresponding Z-axis rotation matrix r TOD v TOD Rotate to the ECEF coordinate system and superimpose the Earth's rotational entrainment velocity ω. E ×r T The satellite's ground-to-solid position vector r at the launch time is obtained. T and the satellite's ground-to-solid velocity vector v at launch time T ; Step S4c uses the satellite's ground-based position vector r at the launch time. T and the satellite's ground-to-solid velocity vector v at launch time T calculate ;r R Here is the location of the ground station, c is the speed of light in a vacuum, and the updated launch time is in Julian days. , ; Among them, JD epoch For reference to the Julian era; t receive_offset To receive the time offset relative to the epoch, and in Return to step S4b and re-execute the complete coordinate transformation, with at least two iterations.

[0015] Furthermore, the radial velocity ρ in step S5 is calculated as follows: the line-of-sight vector is calculated in the ECEF coordinate system. and radial velocity , The dot product of the line-of-sight vector and the velocity vector. Line-of-sight distance; The calculation of the Doppler frequency offset Δf in step S6 is as follows: f0 is the carrier frequency, i.e., the nominal center frequency. The negative sign indicates a positive frequency offset when approaching the center frequency. Let c be the radial velocity and c be the speed of light in a vacuum. Following step S6 is step S7: Azimuth and elevation calculation: ... Project the coordinates onto the Northeast-Sky (ENU) coordinate system with the ground station as the origin, and calculate the azimuth angle. Where the azimuth angle Az is measured clockwise from true north; E0 is the eastward component of the line-of-sight vector; N is the northward component, and the elevation angle is... The pitch angle El is measured upwards from the horizontal plane, and U is the celestial component of the line-of-sight vector. This represents the horizontal distance component.

[0016] Furthermore, in step S4b , , r T and v T The calculation formula is as follows: , The launch time for the k-th iteration is the Julian day; JD epoch For reference to the Julian era; t receive_offset The offset of the received time relative to the epoch, in seconds; τ k The current optical fiber valuation; GMST expansion formula is: ; for Distance J 2000 Julian days of the epoch; T is the Julian century number, used for higher-order precession correction. in ; For distance J 2000 The number of Julian days in the calendar era; 2451545.0 is J. 2000 The corresponding Julian days, T is the distance from J 2000 The Julian century number; 36525.0 is the number of days per Julian century; the time variable used for higher-order correction terms in the GMST expansion; , ε0 is Greenwich Mean Sidereal Time; ε0 is the obliquity of the ecliptic, i.e., the average obliquity of the ecliptic without nutation correction; ΔΨ×cosε0 is the right ascension nutation term, which corrects the mean sidereal time to true sidereal time, in degrees, and needs to be converted to radians before participating in coordinate rotation; the aforementioned and The calculation is performed in each optical travel time iteration using the data updated in the current iteration. Executed independently, and the entire process must not be based on the Julian date of receipt (JD). receive Corresponding fixed sidereal time value substitution; , The satellite's Earth-fixed position vector is obtained by rotating from the TOD coordinate system to the ECEF coordinate system; R z r is the rotation matrix about the Z-axis; TOD The satellite position in the TOD coordinate system; v T v is the satellite's ground-to-solid velocity vector at the time of launch. TOD The velocity in the TOD coordinate system; ω E×r T The correction term for the entrainment velocity caused by the Earth's rotation. Let ω be the rotation matrix about the Z-axis. E This is the Earth's rotational angular velocity. r T The satellite's ground-based position at the launch time obtained in step S4b; the ω E ×r T The term is the velocity correction term caused by the Earth's rotation, and its physical meaning is to make v T The correct expression of the satellite's velocity relative to the Earth's fixed reference frame ensures that the subsequent calculation of the radial velocity ρ reflects the satellite's true relative motion with respect to the ground station.

[0017] Furthermore, in step S4b, the nutation matrix N is constructed using the IAU1980 truncation model, retaining no fewer than five main nutation parameters. These nutation parameters are phase parameters, which are linear combinations of five Delaunay basic variables {l,l',F,D,Ω}, where l,l',F,D,Ω are the mean anomaly angle of the moon, the mean anomaly angle of the sun, the ascending node distance of the moon, the mean distance between the sun and the moon, and the ecliptic longitude of the ascending node of the moon, respectively.

[0018] Furthermore, in step S2, the time system adopts the following correction chain: the deviation between UTC and TAI is taken as the cumulative number of seconds jumps published by the International Bureau of Time; the deviation between TAI and TT is always 32.184 seconds, where TT is the relativistic time scale; the method uses JD when performing calculations of precession, nutation, obliquity of the ecliptic ε0, and polynomial time variable T. TT JD is used when performing GMST and GAST calculations. UTC This allows for the correct distinction between the applicable scope of TT and the Earth's rotation timescale UT, avoiding systematic errors introduced by mixing timescales.

[0019] Furthermore, a computer-readable storage medium stores computer program instructions that, when executed by a processor, implement all the steps of the method.

[0020] Furthermore, a low-Earth orbit satellite Doppler frequency offset accurate prediction system includes the following five functional modules that cooperate and work together: The orbit propagation module is used to obtain the satellite's orbital position in the J orbit by iteratively solving the Kepler equations based on the Kepler orbital parameters and orbital epoch times. 2000 Initial position vector r in the inertial coordinate system J2000 and velocity vector v J2000 ; The precise coordinate transformation module includes an IAU1980 complete precession-nutation submodule and a sidereal time calculation submodule. The IAU1980 complete precession-nutation submodule retains no fewer than five principal nutation terms to calculate the nutation longitude ΔΨ and nutation angle Δε. The sidereal time calculation submodule uses the Julian date of the signal transmission time. As input, Greenwich Mean Time is calculated independently. Greenwich True Sidereal Time and with Complete the TOD coordinate transformation and then proceed to the ECEF coordinate transformation; The optical timing iteration module iteratively solves for the signal transmission time τ using the reception time as the initial value. Each iteration uses the current τ as the starting point. k Updated The precise coordinate transformation module is invoked, and the iteration count is no less than 2 times; The Doppler calculation module, in the ECEF coordinate system, uses the satellite velocity v at the time of launch. T Calculate the radial velocity ρ, where , and according to Output Doppler frequency offset Δf, where f0 is the carrier frequency; c is the speed of light in vacuum; the negative sign indicates a positive frequency offset when approaching; The data output module outputs azimuth, elevation, and Doppler frequency offset prediction sequences to standard output or a text file. The functional modules use 64-bit IEEE754 double-precision floating-point data type to express internal state quantities. The time state is passed between modules through the time offset parameter relative to the epoch, avoiding the loss of floating-point effective bits caused by directly using the Julian day as the parameter. The system does not maintain a global variable time state. The calculation of each forecast time is independent of each other and can support parallel expansion.

[0021] Furthermore, in each iteration step k of the optical time iteration module, the sidereal time calculation submodule calculates independently. For input, The launch time (Julian day) is calculated independently for the k-th iteration, and each iteration takes the updated τ. k Recalculate; do not substitute fixed reception time with sidereal time; calculate Greenwich Mean Time independently. Greenwich Mean Time The value changes dynamically with different iteration steps k, and must not be based on a fixed GAST corresponding to the signal reception time throughout the process. receive Value substitution; the Greenwich true sidereal time The dynamic changes, together with the complete IAU1980 nutation implementation and full double-precision data type of the other modules in the system, work together to ensure the final accuracy of the Doppler frequency offset prediction.

[0022] In summary, this invention achieves the following five technical elements working together in a synergistic effect: Technical Element (A) – Orbital Propagation: Based on the six roots of the Keplerian orbit (a,e,i,Ω,ω,M0) and the epoch time, using the average angular velocity

[0023] The mean anomaly angle M is calculated at each predicted time, and the off-anomaly angle E is obtained by iteratively solving the Kepler equation Ee×sinE=M. Then, the satellite J is obtained from the orbital plane position / velocity through three rotations. 2000 The position vector r of the inertial frame J2000 and velocity vector v J2000 ; Technical Element (B) – IAU1980 Complete Nutation Model: Retains no fewer than five main nutation parameters, using a linear combination of five Delaunay basic variables {l,l',F,D,Ω} as phase parameters to calculate nutation longitude ΔΨ and nutation angle Δε with an accuracy better than 0.001 arcseconds; Technical Element (C) – Coordinate Transformation for Dynamic Update of GAST at Launch Time: J is executed at each predicted time. 2000 →TOD→ECEF coordinate transformation is performed using the Julian date of signal transmission. (rather than the receiving time JD) receive As Greenwich Mean Time Greenwich True Sidereal Time The calculation benchmark, and rotated matrix around the Z-axis. Complete the Z-axis rotation from TOD to ECEF; Technical Element (D) – Iterative Convergence of Optical Travel Time: with τ 0 =0 is the initial value, according to

[0024] Iteration, each iteration using the updated Re-execute the complete coordinate transformation described in feature (C); after at least two iterations, the position residual can converge to the sub-nanosecond level; Feature (E) – Full Double Precision Numerical Implementation: All numerical calculations (including Julian day, nutation parameters, sidereal time angle, matrix elements, matrix multiplication, radial velocity, and Doppler frequency offset) are performed using 64-bit IEEE 754 double-precision floating-point data type, avoiding the loss of significant bits near large Julian day numbers in single-precision floating-point calculations.

[0025] The synergistic relationship of the above five technical elements can be summarized as follows: Element (A) provides the initial state of the satellite in the inertial frame; Element (B) provides accurate nutation quantities ΔΨ and Δε for subsequent coordinate transformations; Element (C) depends on the high-precision nutation results provided by Element (B) and on the successively updated nutation quantities provided by Element (D). The iterative convergence process of feature (D) fully invokes the coordinate transformations defined by features (B) and (C) at every step; the full double precision strategy of feature (E) is applied throughout all numerical steps (A), (B), (C), and (D). Any missing or simplified feature will result in a degradation of overall accuracy.

[0026] The beneficial effects of this invention are: 1. Significantly improved accuracy: After point-by-point comparison using simulation software (2146 data points, 7.2 minutes of transit arc, including near-overhead segment), the mean absolute error of the azimuth angle is 0.0014°, the peak absolute error is 0.0308°, the mean absolute error of the pitch angle is 0.0017°, the peak absolute error is 0.0031°, and the mean error of the Doppler frequency offset is 3.20Hz / peak 4.32Hz (approximately 0.24ppb relative to the 17.8GHz carrier frequency, approximately 10.8ppm for the Doppler peak at 401kHz, and approximately 0.005% of the typical ±40kHz capture bandwidth); 2. High versatility: All algorithms are implemented on general-purpose computers using standard C language double-precision floating-point operations, with no specific hardware dependencies, facilitating cross-platform deployment and subsequent porting to embedded platforms; 3. Verifiable in engineering: The continuous output mode provides point-by-point comparison data compatible with simulation software formats, which can be directly used in the accuracy acceptance stage of engineering delivery; 4. Strong numerical stability: The entire process uses double-precision floating-point operations to eliminate single-precision truncation errors, which is especially important for Doppler forecasting in high-frequency bands such as Ku / Ka. 5. Clear architecture: Each functional module is decoupled with a relative epoch time offset as the interface parameter, with no globally variable time state, and can support parallel expansion. Attached Figure Description

[0027] Figure 1 This is a block diagram of the overall architecture of the forecasting system of this invention; Figure 2 A comparative diagram illustrating the dynamic update principle of GAST at launch time; in, Figure 2 (a) For the existing simplified implementation: fixed GAST receive plan, Figure 2 (b) is the dynamics of this invention. Differences in position projection of the proposed schemes under LEO satellite overpass scenarios; Figure 3 This is a schematic diagram of the optical travel time iterative convergence process; Figure 4 This is a flowchart of the forecasting method of the present invention; Figure 5 The figure shows the point-by-point comparison results between the embodiment of the present invention and the STK simulation, including three types of error curves: azimuth angle, elevation angle, and Doppler frequency deviation. The test conditions are: carrier frequency 17799.96MHz, Shijiazhuang ground station, and typical transit arc segment on 2026-01-09.

[0028] in, Figure 5 (a) is the azimuth error curve. Figure 5 (b) is a pitch angle error curve. Figure 5 (c) is the Doppler frequency offset error curve. The test conditions were a carrier frequency of 17799.96MHz, Shijiazhuang ground station, and a typical transit arc on January 9, 2026. Detailed Implementation

[0029] like Figures 1-5 As shown, in this embodiment, the present invention includes a method for accurate prediction of Doppler frequency offset of low-Earth orbit (LEO) satellites, a storage medium, and a system. The accurate prediction method for LEO satellites is executed by a computer and is applicable to LEO satellite transit prediction scenarios with orbital altitudes of 200-2000 km. The method employs the following complementary and synergistic technical features: (A) orbital propagation based on Kepler satellite orbital parameters and analytical Newton iteration; (B) a complete implementation of the IAU1980 precession-nutation model retaining no fewer than five main nutation parameters; and (C) using the Julian date of signal transmission. Dynamically update Greenwich Mean Time J 2000 →TOD→ECEF coordinate transformation; (D) an optical time convergence process with the receiving time as the initial value and at least two iterations, and the coordinate transformation described in feature (C) is independently called in each optical time iteration; (E) a numerical implementation strategy using 64-bit IEEE 754 double-precision floating-point operations throughout; the above features (A), (B), (C), (D), and (E) work together to ensure that the relative error of the Doppler frequency offset prediction relative to the carrier frequency under Ku-band carrier frequency conditions is no greater than 10. -9 The magnitude (i.e., the ratio of the absolute Doppler frequency offset error to the carrier frequency, corresponding to an absolute error of no more than approximately 18 Hz at approximately 17.8 GHz in the Ku band); the method includes the following steps: S1. Parameter loading: Obtain satellite orbit parameters, orbit epoch time, and ground station coordinate parameters. The satellite orbit parameters are the six orbital elements, which include the semi-major axis a, eccentricity e, orbital inclination i, right ascension of the ascending node Ω, argument of perigee ω, and mean perigee M0. The ground station coordinate parameters include latitude φ, longitude λ, and elevation h. Before conducting Doppler frequency offset prediction, necessary basic data needs to be prepared. Satellite orbital parameters are key information describing the satellite's trajectory, typically including six basic orbital elements: semi-major axis (a), eccentricity (e), orbital inclination (i), right ascension of the ascending node (Ω), argument of perigee (ω), and mean perigee (M0). These parameters collectively determine the shape of the satellite's orbit, its spatial orientation, and its specific position in orbit. The semi-major axis (a) defines the size of the elliptical orbit, and its value directly affects the satellite's orbital period and velocity. The eccentricity (e) describes the flattening of the ellipse; when the eccentricity is zero, the orbit is circular, and the larger the eccentricity, the flatter the ellipse. The orbital inclination (i) represents the angle between the orbital plane and the Earth's equatorial plane; this parameter determines the satellite's operating range and ground coverage area. The right ascension of the ascending node (Ω) indicates the position of the ascending node on the equatorial plane; it is the point where the satellite's orbit crosses the equatorial plane from the Southern Hemisphere into the Northern Hemisphere. The perigee argument ω describes the angular position of the perigee relative to the ascending node, and this angle is measured in the orbital plane. The mean perigee angle M0 represents the satellite's motion in its orbit, and this parameter can be used to calculate the satellite's orbital position at any given time.

[0030] The orbital epoch time is the reference time corresponding to the above orbital parameters. Because satellite orbits are affected by various perturbations such as Earth's non-spherical gravity, atmospheric drag, and solar radiation pressure, orbital parameters accurately describe the satellite's state only at specific moments. Orbital epoch times are usually expressed using the Julian day, a continuous timekeeping system that begins counting from noon on January 1, 4713 BC, avoiding the complexity of month and year changes in a conventional calendar. The Julian day facilitates time difference calculations and time system conversions. In practical applications, satellite operators regularly update orbital parameters and epoch times to ensure the accuracy of orbital predictions. For low Earth orbit (LEO) satellites, due to the significant impact of atmospheric drag, orbital parameters may be updated daily; for geostationary orbit satellites, the orbit is relatively stable, and the update cycle can be extended to several days or weeks.

[0031] Ground station coordinates describe the location of the receiving equipment on the Earth's surface. Common representations include geodetic coordinates and geocentric-fixed coordinates. The geodetic coordinate system uses three parameters to define the ground station's position: latitude φ, longitude λ, and elevation h. Longitude λ represents the east-west position of the ground station relative to the prime meridian, latitude φ represents its north-south position relative to the equatorial plane, and elevation h represents its vertical height relative to a reference ellipsoid. The geocentric-fixed coordinate system uses the Earth's center of mass as its origin and represents the ground station's position using three-dimensional rectangular coordinates. This coordinate system rotates with the Earth's rotation. In Doppler frequency offset calculations, the ground station coordinates need to be transformed to a unified geocentric-fixed coordinate system for vector operations with the satellite position. Accurate ground station positioning directly affects forecast accuracy; positioning errors can lead to deviations in the calculated relative velocity between the satellite and the ground station, thus affecting the Doppler frequency offset forecast results.

[0032] S2. Time Calculation: Determine the Julian Day (JD) corresponding to the predicted time based on the orbital epoch time. receive That is, using the epoch time as a reference, the Julian day for the reception of each forecast time is recursively calculated according to a preset step size. receive ; Doppler frequency offset forecasting requires targeting a specific future time, known as the forecast time. The forecast time is typically specified by the user based on specific application needs and can be a time several minutes, hours, or even days in the future. For ease of subsequent calculations, the forecast time needs to be converted to Julian day format. This conversion process first requires determining the calendar date and time of the forecast time, including year, month, day, hour, minute, and second. Then, calculations are performed according to the definition rules of the Julian day to convert the calendar time into a consecutive number of days starting from the Julian day's starting point.

[0033] In satellite communication, the signal received by the ground station was transmitted by the satellite at a certain point in the past. It takes a certain amount of time for the signal to travel from the satellite to the ground station; this time is called optical travel time. Therefore, the predicted time actually corresponds to the signal reception time, while the actual transmission time is earlier than the reception time. The Julian day of the reception time is the Julian day value corresponding to the predicted time, and this time serves as the reference for subsequent iterative calculations. When determining the Julian day of the reception time, it is necessary to ensure the consistency of the time system. Satellite orbital parameters are usually given based on Coordinated Universal Time (UTC) or International Atomic Time (IAT), but in practical applications, local time or other time systems may be used. Therefore, appropriate time system conversions are required to ensure that all time quantities use a consistent time reference.

[0034] S3. Orbital Propagation: Based on the satellite orbital parameters and the Julian day (JD) of the receiving time... receive By solving the Kepler equations, the satellite's position in J... 2000 The position vector r of the inertial coordinate system J2000 and velocity vector vJ2000 Specifically, the mean anomaly angle M = M0 + n × (t - t0) is calculated from the current forecast time. The average angular velocity

[0035] Iteratively solving the Kepler equation Ee×sinE=M yields the anomalous angle E, which in turn gives the satellite's position at point J. 2000 The position vector r of the inertial coordinate system J2000 and velocity vector v J2000 ; The motion of a satellite in space follows Kepler's laws. By solving Kepler's equations, the satellite's orbital parameters can be used to calculate its position in space. 2000 The position vector r of the inertial coordinate system J2000 and velocity vector v J2000 The position vector r of the inertial coordinate system J2000 and velocity vector v J2000 An inertial coordinate system is a coordinate system that does not rotate with the Earth's rotation. A commonly chosen coordinate system is the Earth's center-of-mass inertial coordinate system, where the origin is located at the Earth's center of mass, the fundamental plane is the Earth's equatorial plane, and the fundamental direction points towards the vernal equinox. In this coordinate system, the satellite's trajectory remains stable, facilitating orbital mechanics calculations.

[0036] S31, based on the mean apogee angle M and average motion parameters in the satellite orbit parameters, and the Julian day JD at the time of reception. receive The time difference with the orbital epoch is used to determine the target mean apogee angle at the current forecast time.

[0037] The mean anomaly angle M describes the satellite's position progression in its orbit, increasing linearly with time. Mean motion parameters define the angular velocity of a satellite orbiting the Earth, and their value is closely related to the semi-major axis of its orbit. According to Kepler's third law, the square of a satellite's orbital period is proportional to the cube of its semi-major axis, from which the formula for calculating the mean motion parameters can be derived. Mean motion parameters are typically expressed in radians per second or orbits per day, reflecting the time required for a satellite to complete one orbit.

[0038] When calculating the target's mean anterior angle M, it is first necessary to determine the Julian day JD at the time of reception. receive The time difference between the reference time and the predicted time. This time difference represents the length of time elapsed from the reference time of the orbital parameters to the predicted time, usually in days. Multiplying the time difference by the average motion parameters yields the increment of the mean anomaly angle of the satellite during this period. The JD at the epoch time... receiveAdding this increment yields the target mean anomaly angle M0 corresponding to the predicted time. Since the mean anomaly angle M is a periodic parameter, its value ranges from zero to 360 degrees or from zero to twice pi. When the calculated result exceeds this range, a modulo operation is required to normalize the angle to the standard range. This calculation process assumes that the satellite moves along an ideal Keplerian orbit and ignores the influence of various perturbations. For short-term forecasts and medium-to-high orbit satellites, this assumption has sufficient accuracy.

[0039] S32, based on the target's mean apogee angle M and the eccentricity e in the satellite orbital parameters, the apogee angle E is obtained by iteratively solving the Kepler equation.

[0040] Kepler's equations describe the relationship between the mean anomaly M, the off-anomaly E, and the eccentricity e. These equations are transcendental equations and cannot be directly solved analytically; they require numerical iteration methods. The off-anomaly E is an auxiliary angular parameter in orbital mechanics, used to accurately calculate the satellite's actual position in its orbit.

[0041] The iterative solution process typically employs Newton's iteration method or a simple iteration method. Newton's iteration method converges quickly and is suitable for orbits with various eccentricities. At the start of the iteration, the target mean anomaly angle M can be used as the initial estimate of the eccentricity angle E. Subsequently, the function value and derivative of the Kepler equation are calculated based on the current estimated eccentricity angle E, and the estimated eccentricity angle E is updated using Newton's iteration formula. This process is repeated until the difference between two consecutive iterations is less than a preset convergence threshold; at this point, the iteration is considered converged, and an accurate value of the eccentricity angle E is obtained. For near-circular orbits, the eccentricity e is close to zero, and the iteration process usually converges within three to five iterations. For high eccentricity orbits, more iterations may be required to achieve the desired accuracy. In practical applications, to prevent non-convergence and computational loops, a maximum iteration limit is usually set. When the number of iterations exceeds this limit, the computation terminates and an anomaly is reported.

[0042] S33. Based on the anomalous point angle E, determine the satellite's position and velocity within the orbital plane, and then transform the satellite's position and velocity within the orbital plane to the inertial coordinate system using a coordinate rotation matrix to obtain the position vector r. J2000 and velocity vector v J2000 .

[0043] S4a. Optical travel time iterative initialization: using the Julian day JD of the received time... receive Let τ be the initial signal and the initial value of the light travel time. 0 =0, initial value of Julian day at launch time =JD receive Perform the optical travel time iteration step; In satellite communication systems, there is a propagation delay in the signal from satellite transmission to ground station reception; this delay is called optical travel time. The magnitude of optical travel time depends on the distance between the satellite and the ground station, as well as the speed of electromagnetic wave propagation in a vacuum. For geostationary orbit satellites, optical travel time is approximately 0.1-2 seconds, while for low Earth orbit satellites, it is typically between a few milliseconds and tens of milliseconds. Accurate calculation of optical travel time is crucial for Doppler frequency offset prediction because Doppler frequency offset reflects the relative motion between the satellite and the ground station at the moment of signal transmission, rather than the state at the moment of signal reception.

[0044] Since the optical travel time itself depends on the satellite's position at the time of transmission, and the satellite's position in turn depends on the transmission time, an interdependent relationship is formed, requiring an iterative solution. At the start of the iteration, since the exact signal transmission time is not yet determined, the Julian date for the reception time is used. receive The initial estimate of the Julian day of the signal transmission time is used as the initial value. This initial assumption is equivalent to assuming the optical travel time is zero, i.e., ignoring signal propagation delay. In subsequent iterations, the optical travel time and the Julian day of the signal transmission time are continuously updated to gradually approach the true value. The optical travel time iteration step is the core component of the entire Doppler frequency offset prediction method; the convergence and accuracy of the iteration directly affect the accuracy of the final prediction result.

[0045] Specifically, the optical travel time iteration process includes several sub-steps such as coordinate system transformation, distance calculation, optical travel time update, and transmission time update. Each iteration requires recalculating the satellite's position based on the currently estimated signal transmission time (Julian day), then calculating a new optical travel time estimate based on that position, and subsequently updating the signal transmission time (Julian day). This iterative process typically needs to be repeated two to three times to converge to the required accuracy because the optical travel time is very small relative to the satellite's orbital period, and the changes in position and time decrease rapidly after each iteration.

[0046] S4b. Complete coordinate transformation of transmission time: In each of the aforementioned optical travel time iteration steps, the signal transmission time corresponding to the Julian date of the previous iteration step k is transformed... Using the IAU1980 precession-nutation model as the time reference, the nutation longitude ΔΨ and nutation angle Δε are calculated by calling no fewer than five main nutation parameters, and based on the Julian date of the signal transmission time. Independent calculation of Greenwich Mean Sidereal Time Greenwich True Sidereal Time According to the nutation parameters and Greenwich Mean Time Position vector r at time J2000 and velocity vector v J2000 After TOD coordinate transformation and then ECEF coordinate transformation, the satellite's ground-fixed position vector r at the launch time is obtained. T and the satellite's ground-to-solid velocity vector v at launch timeT ; To calculate the relative position and velocity between the satellite and the ground station, the satellite's position needs to be transformed from an inertial coordinate system to a geostationary coordinate system. The geostationary coordinate system is a coordinate system fixed to the Earth, rotating with the Earth's rotation, in which the ground station's coordinates remain constant. The transformation from inertial to geostationary coordinates involves corrections for several physical effects, including precession, nutation, rotation, and polar motion. Nutation, in particular, is a small, periodic oscillation of the Earth's axis of rotation relative to its average position caused by the gravitational pull of the Moon and Sun on the Earth's equatorial bulge. Although the nutation effect is small in magnitude, it has a significant impact on high-precision satellite positioning and frequency prediction.

[0047] S4c. Optical Travel Time Update and Iteration: Based on the satellite's ground-to-solid position vector r at the aforementioned launch time... T When updating the optical path with the ground station coordinate parameters, and according to the Julian day JD of the reception time. receive The Julian date for the signal transmission time is iteratively updated according to the optical travel time. Return to the optical travel iteration step until the preset iteration termination condition is met, and the number of iterations is not less than 2. S5. Radial velocity calculation: Based on the satellite ground-fixed position vector r obtained at the launch time when the preset iteration termination condition is met. T Satellite ground-to-solid velocity vector v at launch time T And the radial velocity ρ is determined using the ground station coordinate parameters; Radial velocity ρ is the projected component of the satellite's velocity along the line-of-sight direction, reflecting the satellite's approaching or receding velocity relative to the ground station. A positive radial velocity ρ indicates the satellite is moving away from the ground station, with a negative Doppler frequency offset and a receiving frequency lower than the transmitting frequency. A negative radial velocity ρ indicates the satellite is approaching the ground station, with a positive Doppler frequency offset and a receiving frequency higher than the transmitting frequency.

[0048] Radial velocity is calculated based on the dot product operation. The dot product of the line-of-sight vector and the velocity vector is equal to the product of the magnitudes of the two vectors multiplied by the cosine of the angle between them. Physically, this means the projection of the velocity along the line-of-sight direction multiplied by the magnitude of the line-of-sight vector. Therefore, dividing the dot product by the magnitude of the line-of-sight vector yields the projection component of the velocity along the line-of-sight direction, which is the radial velocity.

[0049] S6. Doppler frequency offset output: Determine the Doppler frequency offset Δf based on the radial velocity ρ; The physical essence of Doppler frequency offset Δf is the shift in the received frequency relative to the transmitted frequency caused by the relative motion between the transmitter and receiver. In satellite communication, Doppler frequency offset Δf is equal to the negative radial velocity ρ multiplied by the carrier frequency and then divided by the speed of light in a vacuum. The carrier frequency is the operating frequency of the satellite signal. Common satellite communication frequency bands include those around 1.5 GHz, 4 GHz to 8 GHz, 12 GHz to 18 GHz, and 26 GHz to 40 GHz. The Doppler frequency offset Δf varies significantly across different frequency bands, with higher frequency bands exhibiting even larger magnitudes of offset. All numerical calculations in steps S1-S6 above are performed using a 64-bit IEEE 754 double-precision floating-point data type to avoid the loss of significant bits in single-precision floating-point numbers near Julian large numbers.

[0050] In this embodiment, step S4b specifically includes: using the signal transmission time corresponding to the current iteration step k. Using time as the reference, the following transformations are performed sequentially: (i) from

[0051] Calculating the Julian Day of Earth Time, among which For Julian days based on Coordinated Universal Time (UTC), TAI is International Atomic Time, UTC is Coordinated Universal Time, TAI-UTC yields the cumulative leap seconds, 32.184 is the fixed difference between Earth Time and International Atomic Time in seconds, and the denominator 86400 is the number of seconds in a day; (ii) call the IAU1980 precession-nutation model, which retains no less than five main nutation parameters, to calculate the nutation longitude ΔΨ and nutation angle Δε; (iii) call the IAU1976 precession matrix P and nutation matrix N, and r J2000 v J2000 Rotate to the true equator and true vernal equinox coordinate system to obtain r TOD v TOD (iv) by Independent calculation of Greenwich Mean Sidereal Time Greenwich True Sidereal Time (v) with Greenwich true stars Corresponding Z-axis rotation matrix r TOD v TOD Rotate to the ECEF coordinate system and superimpose the Earth's rotational entrainment velocity ω. E ×r T The satellite's ground-to-solid position vector r at the launch time is obtained. T and the satellite's ground-to-solid velocity vector v at launch time T ; Step S4c uses the satellite's ground-based position vector r at the launch time. Tand the satellite's ground-to-solid velocity vector v at launch time T calculate ;r R Here is the location of the ground station, c is the speed of light in a vacuum, and the updated launch time is in Julian days. , ; Among them, JD epoch For reference to the Julian era; t receive_offset To receive the time offset relative to the epoch, and in Return to step S4b and re-execute the complete coordinate transformation, with at least two iterations.

[0052] In this embodiment, the radial velocity ρ in step S5 is calculated as follows: the line-of-sight vector is calculated in the ECEF coordinate system. and radial velocity , The dot product of the line-of-sight vector and the velocity vector. Line-of-sight distance; The calculation of the Doppler frequency offset Δf in step S6 is as follows: f0 is the carrier frequency, i.e., the nominal center frequency. Let c be the radial velocity and c be the speed of light in a vacuum. Following step S6 is step S7: Azimuth and elevation calculation: ... Project the coordinates onto the Northeast-Sky (ENU) coordinate system with the ground station as the origin, and calculate the azimuth angle. Where the azimuth angle Az is measured clockwise from true north; E0 is the eastward component of the line-of-sight vector; N is the northward component, and the elevation angle is... The pitch angle El is measured upwards from the horizontal plane, and U is the celestial component of the line-of-sight vector. This represents the horizontal distance component.

[0053] Azimuth (Az) and elevation (El) are two angular parameters describing the satellite's position relative to a ground station. Azimuth (Az) is typically defined as the angle measured clockwise from true north to the satellite's projected direction, ranging from zero to 360 degrees. Elevation (El) is defined as the angle between the line-of-sight vector and the horizontal plane, ranging from -90 to +90 degrees; a positive value indicates the satellite is above the horizontal plane.

[0054] This azimuth and elevation information can be used to drive the ground tracking antenna to align with the satellite, and in conjunction with Doppler frequency offset prediction, achieve efficient reception and tracking of satellite signals.

[0055] Among them, technical feature (A) provides satellite orbital parameters based on Kepler orbital elements and analytical Newton iterations in J... 2000 The initial position vector r in the inertial coordinate system J2000 and velocity vector vJ2000 As the input basis for nutation modeling, coordinate transformation, and optical travel time iteration defined in subsequent technical features (B), (C), and (D); the complete IAU1980 nutation model of technical feature (B) provides the basis for the accurate calculation of ΔΨ and ε in the coordinate transformation at the launch time; the dynamic update of the launch time GAST of technical feature (C) depends on the high-precision nutation quantity provided by feature (B); the optical travel time iteration of technical feature (D) is necessary for the dynamic update of technical feature (C). Value; The full double precision strategy of technical element (E) runs through all numerical steps of technical elements (A), (B), (C) and (D), especially ensuring that no significant bits are lost in the numerical expression of precision-sensitive quantities such as JD large numbers.

[0056] The five features work synergistically to eliminate multiple error sources in existing simplified implementations, such as missing nutation modeling, ambiguity in coordinate transformation timing, lack of iteration in optical travel time, and precision truncation. Figure 5 The results, verified point-by-point by professional simulation software (carrier frequency 17799.96 MHz Ku band, Shijiazhuang ground station, typical LEO transit arc on January 9, 2026, arc length 7.2 minutes, 2146 data points, elevation angle range 10°-85°): azimuth mean absolute error 0.0014°, peak absolute error 0.0308°, elevation mean absolute error 0.0017°, peak absolute error 0.0031°, Doppler frequency offset mean error 3.20Hz, peak error 4.32 Hz (approximately 0.24 ppb relative to 17.8 GHz carrier frequency, approximately 10.8 ppm of the Doppler peak), and position residual less than 1 nanosecond after two iterations of optical travel time.

[0057] In this embodiment, in step S4b , , r T and v T The calculation formula is as follows: , The launch time for the k-th iteration is the Julian day; JD epoch For reference to the Julian era; t receive_offset The offset of the received time relative to the epoch, in seconds; τ k This is an estimate based on the current optical travel time. (GMST expansion; unit is degrees;) for Distance J 2000 Julian days in an epoch; T is the Julian century (used for higher-order precession correction), in degrees; in ( For distance J2000 The number of Julian days in the calendar era; 2451545.0 is J. 2000 The corresponding Julian days), (T is the distance from J) 2000 The Julian century number; 36525.0 is the number of days per Julian century; the time variable used for higher-order correction terms in the GMST expansion). (ε0 is the obliquity of the ecliptic, i.e., the average obliquity of the ecliptic without nutation correction; ΔΨ×cosε0 is the nutation term of right ascension, correcting mean sidereal time to true sidereal time; the unit is degrees, which need to be converted to radians before participating in coordinate rotation); and The calculation is performed in each optical travel time iteration using the data updated in the current iteration. Executed independently, and the entire process must not be based on the Julian date of receipt (JD). receive Corresponding fixed sidereal time value substitution; ( The satellite's Earth-fixed position vector is obtained by rotating from the TOD coordinate system to the ECEF coordinate system; R z r is the rotation matrix about the Z-axis; TOD (Satellite position in TOD coordinate system) (v) T v is the satellite's ground-to-solid velocity vector at the time of launch. TOD The velocity in the TOD coordinate system; ω E ×r T (Correction term for the speed of the Earth's rotation) in Let ω be the rotation matrix about the Z-axis. E The angular velocity of Earth's rotation. (ω) E (r is a scalar value of the Earth's rotational angular velocity, in rad / s) T The satellite's ground-based position at the launch time obtained in step S4b; the ω E ×r T The term is the velocity correction term caused by the Earth's rotation, and its physical meaning is to make v T The correct expression of the satellite's velocity relative to the Earth's fixed reference frame ensures that the subsequent calculation of the radial velocity ρ reflects the satellite's true relative motion with respect to the ground station.

[0058] In this embodiment, the nutation matrix N in step S4b is constructed using the IAU1980 truncation model, retaining no less than five main nutation parameters. The nutation parameters are phase parameters, which are linear combinations of five Delaunay basic variables {l,l',F,D,Ω}, where l,l',F,D,Ω are the mean anomaly of the moon, the mean anomaly of the sun, the ascending node distance of the moon, the mean distance between the sun and the moon, and the ecliptic longitude of the ascending node of the moon, respectively.

[0059] In this embodiment, the time system in step S2 adopts the following correction chain: the deviation between UTC and TAI is taken as the cumulative number of seconds jumps published by the International Bureau of Time; the deviation between TAI and TT is always 32.184 seconds, where TT is the relativistic time scale; the method uses JD when performing calculations of precession, nutation, obliquity of the ecliptic ε0, and polynomial time variable T. TT JD is used when performing GMST and GAST calculations. UTC This allows for the correct distinction between the applicable scope of TT and the Earth's rotation timescale UT, avoiding systematic errors introduced by mixing timescales.

[0060] In this embodiment, a computer-readable storage medium stores computer program instructions that, when executed by a processor, implement all the steps of the method described.

[0061] In this embodiment, a low-Earth orbit satellite Doppler frequency offset accurate prediction system includes the following five functional modules that cooperate and work together: The orbit propagation module is used to obtain the satellite's orbital position in the J orbit by iteratively solving the Kepler equations based on the Kepler orbital parameters and orbital epoch times. 2000 The initial position vector r in the inertial coordinate system J2000 and velocity vector v J2000 ; The precise coordinate transformation module includes an IAU1980 complete precession-nutation submodule and a sidereal time calculation submodule. The IAU1980 complete precession-nutation submodule retains no fewer than five principal nutation terms to calculate the nutation longitude ΔΨ and nutation angle Δε. The sidereal time calculation submodule uses the Julian date of the signal transmission time. As input, Greenwich Mean Time is calculated independently. Greenwich True Sidereal Time and with Complete the TOD coordinate transformation and then proceed to the ECEF coordinate transformation; The optical timing iteration module iteratively solves for the signal transmission time τ using the reception time as the initial value. Each iteration uses the current τ as the starting point. k Updated The precise coordinate transformation module is invoked, and the iteration count is no less than 2 times; The Doppler calculation module, in the ECEF coordinate system, uses the satellite velocity v at the time of launch. T Calculate the radial velocity component , and according to Output Doppler frequency offset Δf; The data output module outputs azimuth, elevation, and Doppler frequency offset prediction sequences to standard output or a text file. The functional modules use 64-bit IEEE754 double-precision floating-point data type to express internal state quantities. The time state is transmitted between modules through the time offset parameter relative to the epoch (in seconds), avoiding the loss of floating-point effective bits caused by directly using the Julian day as the parameter. The system does not maintain a global variable time state. The calculation of each forecast time is independent of each other and can support parallel expansion.

[0062] In this embodiment, in each iteration step k of the optical time iteration module, the sidereal time calculation submodule calculates independently. Use as input to independently calculate Greenwich Mean Time (GMT). Greenwich Mean Time The value changes dynamically with different iteration steps k, and must not be based on a fixed GAST corresponding to the signal reception time throughout the process. receive Value substitution; the Greenwich true sidereal time The dynamic changes, together with the complete IAU1980 nutation implementation and full double-precision data type of the other modules in the system, work together to ensure the final accuracy of the Doppler frequency offset prediction.

[0063] Overall Implementation Process This invention is implemented in standard C language on a general-purpose x86-64 computer platform, conforming to the ISO C99 standard, and can be compiled by GCC, Clang, or MSVC. The overall implementation process corresponds to steps S1-S7 of the method, and follows... Figure 4 The complete flowchart shown is executed. The following section uses two pseudocode snippets to illustrate the synergistic effect of features (C) and (D).

[0064] [Typical Simplified Implementation] (This is only used to illustrate the necessity of features (C) and (D), and is not part of the present invention): gast_fixed = calc_GAST(JD_receive); / / Only once, out-of-iteration computation tau=0; for(k=0;k <N;k++){ r_j2000=propagate(t_receive-tau); r_ecef=R_z(-gast_fixed)*NP*r_j2000; / / Time-based semantic confusion tau=|r_ecef-r_R| / c; } The position residual achieved in the above simplified implementation is proportional to the satellite's radial velocity, and there is still an unavoidable GAST angle difference when the number of iterations approaches infinity.

[0065] [Implementation of this invention] (corresponding to the call relationship between the source code obs_geometry_lighttime and j2k_state_to_ecef): tau=0; for(k=0;k<=N;k++){ t_tx_off = t_receive - tau; jd_tx = JD_epoch + t_tx_off / 86400.0; / / Updates with iteration state=orbital_to_j2000_state(sat,t_tx_off); j2k_state_to_ecef(state,jd_tx,&r_T,&v_T); / / Internally computes GAST_tx independently. tau=|r_T-r_R| / c; } 5.2 Specific Implementation of the IAU1980 Nutation Model (Feature B) The nutation sequence of feature (B) adopts the first few main terms of the IAU1980 standard (referencing the IERS Conventions (1996) nutation model specification and the nutation implementation conventions of the IAUSOFA mathematical standard library), wherein the embodiment retains ten terms (NUTATION_TERMS=10), each term using a linear combination of integer coefficients of the five Delaunay basic variables {l, l', F, D, Ω} as phase parameters, corresponding to the nutation longitude amplitude Ap and the nutation angle amplitude Be. Each Delaunay basic variable is calculated according to the following formula (in arcseconds, T is J). 2000 (the starting point of the Julian century)

[0066] The nutation longitude ΔΨ and nutation angle Δε are obtained by summing the amplitudes of each term according to the following formula (the first ten terms are retained in the example, and the amplitude unit is 0.0001 arcseconds): ΔΨ=Σ k (Ap k +Adp k ×T)×sin(Arg k ),Δε=Σ k (Be k +Bdek ×T)×cos(Arg k ) Arg k =nl k ×l+nl' k ×l'+nF k ×F+nD k ×D+nΩ k ×Ω Note: Retaining ten items (example values) is sufficient to control the nutation modeling residuals below 0.001 arcseconds, corresponding to an azimuth / elevation angle error contribution of less than 0.001°; "no less than five items" is the lower limit of the protection range, and the NUTATION_TERMS parameter can be flexibly configured according to the accuracy requirements during engineering implementation.

[0067] 5.3 Detailed Explanation of the Full-Double Precision Numerical Strategy (Feature E) This embodiment uses 64-bit IEEE 754 double-precision floating-point (C language double type, relative precision approximately 2.22 × 10⁻⁶) throughout. -16 The accuracy requirements of the main computational workloads and the implementation strategy of this invention are compared in Table 1: Table 1. Accuracy Requirements and Implementation Strategies for Each Computational Quantity in This Invention

[0068] The accuracy of the Julian day is the most critical: a 1-microsecond time error corresponds to approximately 1.2 × 10⁻⁶ Julian days. -11 Heaven, single-precision IEEE 754 floating-point (float, approximately 7 decimal significant digits, relative precision approximately 1.19 × 10⁻⁶) -7 The relative truncation precision near the Julian day order of 2451545 is 10. -7 The magnitude is insufficient to express the 10^10 microseconds of time precision required. -11 The relative magnitude is far from meeting the requirements. This invention uses 64-bit IEEE 754 double-precision floating-point throughout (double, with a relative precision of approximately 2.22 × 10⁻⁶). -16 This guarantees that JD's relative accuracy is better than 10. -15 Absolute time accuracy can reach 10. -9 Up to 10 -11 Massive engineering needs, covering microsecond-level time precision.

[0069] 5.4 Experimental Verification Results of Examples and Simulation Software Test conditions: Shijiazhuang ground station (38.0559°N, 114.3554°E, altitude 101.5m), carrier frequency 17799.96MHz (Ku band), satellite orbital elements (double-precision reference implementation): semi-major axis 6826302.5884m (altitude approximately 448km), eccentricity 6.4781×10 -4 The simulation parameters are: orbital inclination 55.01617°, perigee argument 16.85885°, ascending node right ascension 13.97521°, mean anomaly 96.57491°, epoch: January 9, 2026, 19:42:47.000 UTC, and the transit arc from 19:43:43.000 to 19:50:52.000 UTC on January 9, 2026 (arc length approximately 7.2 minutes, elevation range 10.02°-84.87°, including the near-overpass segment). The simulation step size is 0.2 seconds, with a total of 2146 valid data points. Each point is compared with the simulation software output.

[0070] Table 2. Point-by-point comparison results between the embodiments of the present invention and the simulation software.

[0071] Note 1: The above error values ​​are obtained by comparing the measured output of the C source code of this invention with the simulation software point by point, covering the complete transit arc of 2146 data points. The Doppler frequency offset error is positive (1.15-4.32Hz) and exhibits a one-sided distribution (standard deviation of about 1.0Hz, mean of about 3.2Hz), which is a systematic deviation rather than random noise. The relative error of this deviation to the 17.79996GHz carrier frequency is about 0.24ppb, which is about 10.8ppm of the Doppler peak (about 401kHz), much smaller than the typical signal acquisition bandwidth (±40kHz) and frequency search step size (typically 1kHz), which fully meets the requirements of engineering applications. For details of its physical source, please refer to Note 2 in Section 5.5.

[0072] Note 2: The mean absolute value of the azimuth error is 0.0014°, the peak absolute error is 0.0308°, and the RMS is 0.0044°; the mean absolute value of the elevation error is 0.0017°, the peak absolute error is 0.0031°, and the RMS is 0.0018°. The Az / El peaks all occur near the top of the overpass arc (elevation angle above 80°). This is because, after enabling optical travel time iteration in this invention, the azimuth / elevation angles are calculated based on the "satellite geometric position at the time of transmission" (i.e., the antenna is pointing towards the "position of the satellite at the time of signal transmission," which is correct in engineering semantics). This differs from the simulation software's default output of azimuth / elevation angles based on the "geometric position at the time of reception." satThe angle definition difference of ×τ≈35m / range≈0.03° is not a loss of algorithm accuracy; the mean azimuth / elevation error of the remaining arc segments is only about 0.0002°, exhibiting stable weak noise characteristics, proving that the synergistic effect of the IAU1980 complete nutation modeling, GAST dynamic update at the launch time, optical travel time iterative convergence, and full double-precision floating-point arithmetic defined by features (B), (C), (D), and (E) of this invention effectively guarantees the angle accuracy. The accuracy of the entire arc segment meets the corresponding engineering objectives (azimuth target <0.05°, elevation target <0.01°). The optical travel time iterative convergence performance is outstanding, converging to the sub-nanosecond position residual level in 2 iterations, and further iterations do not bring additional accuracy improvement, proving that the "2 iterations" adopted in this invention is the optimal engineering value.

[0073] 5.5 Quantitative Comparison of the Necessity of Synergistic Effect of Five Features To illustrate the irreplaceable nature of the "combined synergy" of the five features of this invention, this section presents the accuracy degradation data compared with simulation software after removing each feature in turn. The simulation conditions are exactly the same as in Section 5.4: the same number of orbital elements (semi-major axis 6826302.5884m, eccentricity...). 6.4781×10 -4 The data includes the following parameters: tilt angle 55.01617°, perigee argument 16.85885°, right ascension of ascending node 13.97521°, mean anomaly 96.57491°, epoch 2026-01-09 19:42:47 UTC; the same transit arc (19:43:43-19:50:52 UTC, 2146 data points); and the same carrier frequency and base station. Table 3 presents a point-by-point comparison between the Python reference implementation output and the simulation software output under the corresponding configuration. The data in Table 3 are consistent with those in Table 2 and can be directly compared.

[0074] Table 3 Comparison of the synergistic necessity of feature combinations

[0075] Note 1: The data listed in Table 3 are the absolute error statistics obtained by comparing the point-by-point output of each configuration in the transit arc segment (a total of 2146 data points) from 19:43:43 to 19:50:52 UTC on 2026-01-09 with the output of the simulation software, which is consistent with the scope of Table 2; the complete configuration (i.e., baseline row) of the Python reference implementation used in this table corresponds to each function in the C source code of this invention, and the Az / El / Dop / slant range values ​​output by the two under the same input are completely equivalent at the algorithm level. The measured point-by-point difference is on the order of magnitude: the peak azimuth angle is approximately 3.5 × 10⁻⁶. -12 The peak pitch angle is approximately 4.4 × 10⁻⁶ degrees. -13 The slope and angle are completely identical (difference is 0), and the Doppler peak value is approximately 1.0 × 10⁻⁶. -5Hz, peak optical travel time approximately 1.6 × 10⁻⁶ -17 The accuracy of the C implementation and the Python reference implementation is at the 64-bit IEEE 754 Floating-Point Last Bit (ULP) level, demonstrating that they can serve as precision benchmarks for each other.

[0076] Note 2: The non-zero error in the baseline (the complete invention) is not a loss of algorithm accuracy, but rather a result of the two-body Kepler orbital propagation model used in this invention and the high-fidelity perturbation model (including J2, J3) used in the simulation software. 3、 This is due to differences in the mechanical models between (such as atmospheric drag); this difference manifests as a stable systematic deviation of approximately 3.2 Hz, which is unrelated to the accuracy of the algorithm.

[0077] Note 3: The absence of feature D (without iterative optical travel time) causes the absolute error of the Doppler mean to jump from 3.20Hz at the baseline to 7.00Hz (+119%), and the peak value to jump from 4.32Hz to 11.37Hz (+163%), a significant change in magnitude. This is because the lack of optical travel time iteration simultaneously destroys two key benchmarks: the position at the time of transmission and the GAST at the time of transmission. Note: When D is missing, the peak values ​​of azimuth / elevation angles decrease from 0.0308° / 0.0031° to 0.0150° / 0.0018°—this is because after enabling optical travel time iteration in this invention, the azimuth / elevation angles are calculated based on the satellite position at the time of transmission, while the simulation software defaults to outputting the azimuth / elevation angles based on the geometric position at the time of reception. There is a difference in the azimuth / elevation angles during the overhead segment. sat The angle definition difference is on the order of ×τ≈35m / range≈0.03°; this 0.03° difference represents the difference in physical meaning between the two angle definitions, not a loss of algorithmic accuracy (the transmission time position definition adopted in this invention corresponds to "where the antenna should be pointing to receive the current signal" in engineering, and its physical semantics are correct); when D is missing, the algorithm degenerates into the reception time position, which happens to match the definition in the simulation software, so the Az / El peak value actually decreases. Table 3, with its Doppler direction data, reflects the true difference in algorithmic accuracy and is the core evidence for the necessity of D.

[0078] Note 4: Missing feature C (using the GAST of the receiving time) receive Alternate launch time This resulted in an increase in the peak azimuth angle from 0.0308° to 0.0318° (+3.2%), the mean Doppler frequency from 3.20Hz to 3.25Hz (+1.6%), and the peak Doppler frequency from 4.32Hz to 4.34Hz; all five indicators deteriorated across the board without any directional reversal. This, based on measured data, demonstrates that the dynamic update of GAST at the transmission time (characteristic C) has a quantifiable accuracy advantage over the simplified scheme of GAST at a fixed reception time. The physical mechanism is as follows: when C is missing, the algorithm will update the frequency at the transmission time... The GAST evaluation was incorrectly performed at the receiving time t. receive Evaluation, introducing approximately ω EA constant phase shift of ×τ≈0.05 arcseconds (corresponding to an ECEF tangential offset of approximately 1 meter) is amplified by geometric projection in the over-the-top segment into an azimuth degradation of approximately 0.001°, and transmitted by radial velocity projection into a mean Doppler frequency shift degradation of approximately 0.05 Hz; simultaneously, the peak azimuth angle deteriorates from 0.0308° to 0.0318° in an increment of 0.001°, corresponding to ω. E The angular magnitude of ×τ matches, representing a true reflection of the phase shift in geometric projection. In long arc segments, low-orbit high-speed scenarios, or under different orbital geometries, this shift will be further amplified, with degradation reaching several times that of the test scenario. Feature C's assertion of "launch time value" is a constraint on the consistency of the invention's engineering implementation—without this constraint, independent encoding by different implementers may yield values ​​containing ω. E Inconsistent output with a deviation on the order of ×τ; Feature C and Feature D work closely together (D provides the necessary dynamic updates). For C to use for recalculation Both are indispensable and directly reflect the "inseparable dependency between features," which is also one of the core arguments for the "combined synergy" of the five features proposed in this invention.

[0079] Note 5: The absence of feature B (nutation sequence truncated to 2 terms) increases the absolute error of the Doppler mean from 3.20Hz to 3.42Hz (on the order of +7% and +0.22Hz) and the Doppler peak value from 4.32Hz to 4.36Hz (+0.04Hz). The final perturbation of the azimuth / elevation peak value—B—provides a refinement guarantee outside of the nutation main term (accounting for about 5% of the weight of the complete IAU1980 sequence), contributing about ±0.2Hz to the Doppler accuracy. The absence of feature E (GMST key quantity downgraded to IEEE 754 single precision) caused the peak azimuth angle to increase from 0.0308° to 0.0312° (+1.3%), and the final position of the peak elevation angle was disturbed. The Doppler direction remained at the baseline level because the chain accumulation of Julian heliotropy + GMST single precision truncation under the 2146 data points of this test had not yet fully manifested. The engineering necessity of E lies in: single precision floating point (approximately 1.19 × 10⁻⁶). -7 The relative precision (to the order of Julian day ≈ 2451545) is 10. -7 The order of magnitude is insufficient to express the 10^10 microseconds of time precision required. -11 Relative magnitude, in long arc segments (>10) 4In scenarios such as point-to-point testing, embedded platforms, and fine frequency offset acceptance, uncontrollable ULP-level cumulative deviations will inevitably occur. The test conditions in this paper are just below the E influence manifestation threshold and do not mean that E can be omitted in engineering use. The results from Tables 2 and 3 show that the complete invention (A+B+C+D+E) has achieved an engineering-grade accuracy level that matches the simulation software to the final level of azimuth / elevation angles (Az peak 0.0308°, E1 peak 0.0031°, including the contribution of the angle definition difference in the overpass segment) and sub-10 Hz (Dop peak -4 Hz). Turning off any one of the technical features (B), (C), (D), and (E) will cause an accuracy degradation ranging from sub-millidegrees to tens of hertz, and the direction of degradation is highly consistent (except for feature (D) which causes Az / E1 to reverse due to the angle definition difference, as explained in Note 3, the other three technical features (B), (C), and (E) all show a comprehensive degradation of all indicators or a final disturbance). There are no counterexamples of "turning off a certain feature actually improves accuracy".

[0080] 5.6 Typical Application Scenarios 1. Ground station antenna guidance: Generate azimuth / elevation time series before satellite transit to drive antenna servo controller for high-precision tracking; 2. Pre-acquisition Doppler compensation: The Doppler forecast is sent to the receiver to compress the carrier frequency search range from ±40kHz to ±100Hz, thereby improving the probability of first acquisition and acquisition delay performance. 3. Accuracy Verification Benchmark: The accuracy of the algorithm is compared point by point with professional simulation software such as STK in continuous output mode to quantify the accuracy of the algorithm and provide a quantitative basis for the acceptance of engineering projects. 4. Embedded Porting Reference: The double-precision reference implementation of this invention on a general-purpose computer can serve as a numerical correctness benchmark for the embedded version (FPGA / DSP) to verify the precision loss of the ported code.

[0081] Although the embodiments of the present invention are described with reference to actual solutions, they do not constitute a limitation on the meaning of the present invention. Modifications to the embodiments and combinations with other solutions based on this specification will be obvious to those skilled in the art.

Claims

1. A method for accurate prediction of Doppler frequency offset of low-Earth orbit satellites, characterized in that, It includes the following steps: S1. Obtain satellite orbit parameters, orbit epoch times, and ground station coordinate parameters; S2. Determine the Julian day (JD) corresponding to the predicted time based on the orbital epoch time. receive ; S3. Based on the satellite orbit parameters and the Julian day (JD) of the receiving time. receive By solving the Kepler equations, the satellite's position in J... 2000 The position vector r of the inertial coordinate system J2000 and velocity vector v J2000 ; S4. Perform the optical timing iteration step, calling the IAU1980 precession-nutation model with no fewer than five main nutation parameters, based on the Julian date of the signal transmission time. Dynamically update Greenwich Mean Time In J 2000 The satellite's Earth-fixed position vector r at launch time is obtained by performing a TOD coordinate transformation followed by an ECEF coordinate transformation on the inertial coordinate system. T and the satellite's ground-to-solid velocity vector v at launch time T And according to the Julian date of receipt JD receive The Julian date for the signal transmission time is iteratively updated during the light-time process. Return to the optical travel iteration steps until the preset iteration termination condition is met; S5. Based on the satellite ground-fixed position vector r obtained at the launch time when the preset iteration termination condition is met. T Satellite ground-to-solid velocity vector v at launch time T And the radial velocity ρ is determined using the ground station coordinate parameters; S6. Determine the Doppler frequency deviation Δf based on the radial velocity ρ; All numerical calculations in steps S1-S6 above are performed using a 64-bit IEEE 754 double-precision floating-point data type to avoid the loss of significant bits in single-precision floating-point numbers near Julian large numbers.

2. The method for accurate prediction of Doppler frequency offset of low-orbit satellites according to claim 1, characterized in that: In step S1, the satellite orbit parameters are the six orbital elements, which include the semi-major axis a, eccentricity e, orbital inclination i, right ascension of the ascending node Ω, argument of perigee ω, and mean perigee M0. The ground station coordinate parameters include latitude φ, longitude λ, and elevation h. In step S2, the Julian day of the receiving time for each forecast time is calculated recursively according to a preset step size, based on the epoch time. receive ; Step S3 involves calculating the mean angle from the current forecast time. Where M is the mean anomaly angle at the current time t; M0 is the mean anomaly angle at epoch t0; and n is the average angular velocity; additionally... ， GM is the Earth's gravitational constant; the Kepler equations are solved iteratively. We obtain the near-point angle E, and then obtain the satellite's position at J. 2000 The position vector r of the inertial coordinate system J2000 and velocity vector v J2000； Step S4 includes the following steps: S4a. The received time is Julian day JD. receive Let τ be the initial signal and the initial value of the light travel time. 0 =0, initial value of Julian day at launch time =JD receive Perform the optical travel time iteration step; S4b. In each of the aforementioned optical travel time iteration steps, based on the Julian date of the signal transmission time corresponding to the previous iteration step k. Using the IAU1980 precession-nutation model as the time reference, the nutation longitude ΔΨ and nutation angle Δε are calculated by calling no fewer than five main nutation parameters, and based on the Julian date of the signal transmission time. Independent calculation of Greenwich Mean Sidereal Time Greenwich True Sidereal Time According to the nutation parameters and Greenwich Mean Time Position vector r at time J2000 and velocity vector v J2000 After TOD coordinate transformation and then ECEF coordinate transformation, the satellite's ground-fixed position vector r at the launch time is obtained. T and the satellite's ground-to-solid velocity vector v at launch time T ; S4c. Based on the satellite's ground-based position vector r at the launch time. T When updating the optical path with the ground station coordinate parameters, and according to the Julian day JD of the reception time. receive The Julian date for the signal transmission time is iteratively updated according to the optical travel time. Then, return to execute the optical time iteration steps until the preset iteration termination condition is met, and the number of iterations is not less than 2.

3. The method for accurate prediction of Doppler frequency offset of low-orbit satellites according to claim 2, characterized in that: Step S4b specifically includes: using the signal transmission time corresponding to the current iteration step k. Using time as the reference, the following transformations are performed sequentially: (i) from ， Calculating Earth Time Julian Day ,in (ii) Using the Julian Day based on Coordinated Universal Time (UTC), TAI is International Atomic Time, UTC is Coordinated Universal Time, TAI-UTC, to obtain the cumulative leap seconds; (iii) Using the IAU1980 precession-nutation model that retains no less than five main nutation parameters, calculate the nutation longitude ΔΨ and nutation angle Δε; (iii) Using the IAU1976 precession matrix P and nutation matrix N, r J2000 v J2000 Rotate to the true equator and true vernal equinox coordinate system to obtain r TOD v TOD (iv) by Independent calculation of Greenwich Mean Sidereal Time Greenwich True Sidereal Time (v) with Greenwich true stars Corresponding Z-axis rotation matrix r TOD v TOD Rotate to the ECEF coordinate system and superimpose the Earth's rotational entrainment velocity ω. E ×r T The satellite's ground-to-solid position vector r at the launch time is obtained. T and the satellite's ground-to-solid velocity vector v at launch time T ; Step S4c uses the satellite's ground-based position vector r at the launch time. T and the satellite's ground-to-solid velocity vector v at launch time T calculate ;r R Here is the location of the ground station, c is the speed of light in a vacuum, and the updated launch time is in Julian days. , ； Among them, JD epoch For reference to the Julian era; t receive_offset To receive the time offset relative to the epoch, and in Return to step S4b and re-execute the complete coordinate transformation, with at least two iterations.

4. The method for accurate prediction of Doppler frequency offset of low-Earth orbit satellites according to claim 3, characterized in that: The radial velocity ρ in step S5 is calculated as follows: the line-of-sight vector is calculated in the ECEF coordinate system. and radial velocity , The dot product of the line-of-sight vector and the velocity vector. Line-of-sight distance; The calculation of the Doppler frequency offset Δf in step S6 is as follows: f0 is the carrier frequency, i.e., the nominal center frequency. Where is the radial velocity, and c is the speed of light in vacuum; the negative sign indicates a positive frequency offset when approaching. Following step S6 is step S7: Azimuth and elevation calculation: ... Project the coordinates onto the Northeast-Sky (ENU) coordinate system with the ground station as the origin, and calculate the azimuth angle. Where the azimuth angle Az is measured clockwise from true north; E0 is the eastward component of the line-of-sight vector; N is the northward component, and the elevation angle is... The pitch angle El is measured upwards from the horizontal plane, and U is the celestial component of the line-of-sight vector. This represents the horizontal distance component.

5. The method for accurate prediction of Doppler frequency offset of low-orbit satellites according to claim 3, characterized in that, In step S4b , , r T and v T The calculation formula is as follows: ， The launch time for the k-th iteration is the Julian day; JD epoch For reference to the Julian era; t receive_offset τ is the offset of the receiving time relative to the epoch. k The current optical fiber valuation; GMST expansion formula is: ； for Distance J 2000 Julian days in an epoch; T is the Julian century number, used for higher-order precession correction. in For distance J 2000 The number of Julian days in the calendar era; 2451545.0 is J. 2000 The corresponding Julian days, T is the distance from J 2000 The Julian century number; 36525.0 is the number of days per Julian century; the time variable used for higher-order correction terms in the GMST expansion; ; ε0 is Greenwich Mean Sidereal Time; ε0 is the obliquity of the ecliptic, i.e., the average obliquity of the ecliptic without nutation correction; ΔΨ×cosε0 is the right ascension nutation term, which, when correcting mean sidereal time to true sidereal time, needs to be converted to radians before participating in coordinate rotation; the aforementioned and The calculation is performed in each optical travel time iteration using the data updated in the current iteration. Executed independently, and the entire process must not be based on the Julian date of receipt (JD). receive Corresponding fixed sidereal time value substitution; , The satellite's Earth-fixed position vector is obtained by rotating from the TOD coordinate system to the ECEF coordinate system; R z r is the rotation matrix about the Z-axis; TOD The satellite position in the TOD coordinate system; ;where v TOD Let ω be the velocity in the TOD coordinate system. E ×r T The correction term for the entrainment velocity caused by the Earth's rotation. Let ω be the rotation matrix about the Z-axis. E This is the Earth's rotational angular velocity. r T The satellite's ground-based position at the launch time obtained in step S4b; the ω E ×r T The term is the velocity correction term caused by the Earth's rotation, and its physical meaning is to make v T The correct expression of the satellite's velocity relative to the Earth's fixed reference frame ensures that the subsequent calculation of the radial velocity ρ reflects the satellite's true relative motion with respect to the ground station.

6. The method for accurate prediction of Doppler frequency offset of low-orbit satellites according to claim 1, characterized in that: In step S4b, the nutation matrix N is constructed using the IAU1980 truncation model, retaining no fewer than five main nutation parameters. The nutation parameters are phase parameters, which are linear combinations of five Delaunay basic variables {l,l',F,D,Ω}, where l,l',F,D,Ω are the mean anomaly angle of the moon, the mean anomaly angle of the sun, the ascending node distance of the moon, the mean distance between the sun and the moon, and the ecliptic longitude of the ascending node of the moon, respectively.

7. The method for accurate prediction of Doppler frequency offset of low-orbit satellites according to claim 1, characterized in that: In step S2, the time system adopts the following correction chain: the deviation between UTC and TAI is taken as the cumulative number of seconds jumps published by the International Bureau of Time; the deviation between TAI and TT is always 32.184 seconds, where TT is the relativistic time scale; the method uses JD when performing calculations of precession, nutation, obliquity of the ecliptic ε0, and polynomial time variable T. TT JD is used when performing GMST and GAST calculations. UTC This allows for the correct distinction between the applicable scope of TT and the Earth's rotation timescale UT, avoiding systematic errors introduced by mixing timescales.

8. A computer-readable storage medium storing computer program instructions that, when executed by a processor, implement all the steps of the method as claimed in any one of claims 1 to 7.

9. A low-Earth orbit satellite Doppler frequency offset accurate prediction system applied to the method described in any one of claims 1-7, characterized in that: It includes the following five functional modules that work together in a coordinated manner: The orbit propagation module is used to obtain the satellite's orbital position in the J orbit by iteratively solving the Kepler equations based on the Kepler orbital parameters and orbital epoch times. 2000 Initial position vector r in the inertial coordinate system J2000 and velocity vector v J2000 ; The precise coordinate transformation module includes the IAU1980 complete precession-nutation submodule and the sidereal time calculation submodule; The optical timing iteration module iteratively solves for the signal transmission time τ using the reception time as the initial value. Each iteration uses the current τ as the starting point. k Updated The precise coordinate transformation module is invoked, and the iteration count is no less than 2 times; The Doppler calculation module, in the ECEF coordinate system, uses the satellite velocity v at the time of launch. T Calculate the radial velocity ρ; The data output module outputs azimuth, elevation, and Doppler frequency offset prediction sequences to standard output or a text file. The functional modules use 64-bit IEEE754 double-precision floating-point data type to express internal state quantities throughout, and time states are passed between modules through time offset parameters relative to the epoch.

10. A precise prediction system for low-Earth orbit satellite Doppler frequency offset according to claim 9, characterized in that: In each iteration step k of the optical time iteration module, the sidereal time calculation submodule calculates independently. Use as input to independently calculate Greenwich Mean Time (GMT). Greenwich Mean Time The value changes dynamically with different iteration steps k, and must not be based on a fixed GAST corresponding to the signal reception time throughout the process. receive Value substitution; The IAU1980 complete precession-nutation submodule retains no fewer than five principal nutation terms to calculate the nutation longitude ΔΨ and nutation angle Δε; the sidereal time calculation submodule uses the Julian day of the signal transmission time. As input, Greenwich Mean Time is calculated independently. Greenwich True Sidereal Time and with Complete the TOD coordinate transformation and then proceed to the ECEF coordinate transformation.

Images