CN116125503A

CN116125503A - High-precision satellite orbit determination and prediction algorithm

Info

Publication number: CN116125503A
Application number: CN202211231390.4A
Authority: CN
Inventors: 张伟; 刘胶胶; 靳江旭; 梁博; 胡慧; 李亚男; 姜俊涛
Original assignee: Shaanxi Aerospace Technology Application Research Institute Co ltd
Current assignee: Shaanxi Aerospace Technology Application Research Institute Co ltd
Priority date: 2022-10-09
Filing date: 2022-10-09
Publication date: 2023-05-16

Abstract

The invention discloses a high-precision satellite orbit determination and prediction algorithm, which comprises the following steps: orbit determination, processing external measurement and GNSS internal measurement data, calculating the precise orbit of a satellite, and providing orbit numbers for injection parameters, orbit prediction and orbit control; and carrying out the track forecast, the ephemeris forecast is calculated by extrapolation by using a precise track, and on the basis, carrying out relevant forecast calculations such as the tracking condition forecast of measurement and control equipment, the equipment guiding forecast, the satellite point forecast, the ground moon shadow forecast and the like. The method can more quickly and accurately determine the instantaneous orbit parameters of a large number of satellites, and can more accurately forecast the orbits for a certain period of time in the future, thereby realizing accurate and effective control of the satellites and avoiding collision risks, and simultaneously providing reliable, effective and quick technical support for engineering tasks such as telemetry data reception, satellite platform state monitoring, load data transmission and the like of the satellites.

Description

High-precision satellite orbit determination and prediction algorithm

Technical Field

The invention belongs to the field of satellite application, and particularly relates to a high-precision satellite orbit determination and prediction algorithm.

Background

With the continuous development of the national commercial aerospace industry, more number and more functions of satellite constellations are sequentially launched and put into use, the number of in-orbit satellites is exploded, the orbit parameters of the satellites are different, especially in medium-low orbits, the orbit resources of the satellites are extremely preferential, and the satellite arrangement is crowded, so that the risks of mutual collision and interference among the satellites are greatly increased. Therefore, it is urgently needed to explore a high-precision satellite orbit determination and prediction algorithm, which can more quickly and accurately determine instantaneous orbit parameters of a large number of satellites and perform more accurate orbit prediction for a certain period of time in the future, so as to realize accurate and effective control of the satellites and avoid collision risks, and simultaneously provide reliable, effective and quick technical support for engineering tasks such as telemetry data reception, satellite platform state monitoring, load data transmission and the like of the satellites.

Disclosure of Invention

In order to solve the technical problems, the invention discloses a high-precision satellite orbit determination and prediction algorithm, which comprises the following steps:

s1, preprocessing data, namely preprocessing satellite orbit measurement data to generate a data file in a format required by orbit calculation;

s2, initial orbit calculation, namely calculating and generating initial orbit information of the satellite by using orbit measurement data for orbit improvement;

s3, performing track improvement, namely performing iterative improvement on the initial track by using measurement data based on a track dynamics method to generate a precise track;

s4, ephemeris forecast, according to the number of satellite orbits, extrapolating the state of the satellite in a period of time in the future;

s5, station tracking forecast, namely calculating station tracking forecast based on ephemeris forecast results and combining with geographic coordinates of the station measurement and control equipment;

s6, equipment guiding forecasting, namely calculating the guiding angle of the antenna tracking satellite of the measurement and control station equipment based on ephemeris forecasting results and combined with the geographic coordinates of the measurement and control station equipment;

s7, forecasting the satellite lower points, and calculating satellite lower point positions based on ephemeris forecasting results;

s8, earth and moon shadow forecasting, and calculating the time period of the satellite entering the earth or moon shadow area based on the ephemeris forecasting result.

Preferably, in the step S1, the specific step of data preprocessing includes:

the method comprises the steps of correcting troposphere delay;

ionosphere delay correction;

correcting the data time scale;

delay correction of equipment;

removing the wild value;

and (5) converting the data format.

Preferably, in S2, for different data sources, the initial track is calculated according to the classification of the type of the track-setting method.

Preferably, the method for calculating the primary track according to the classification treatment of the type of the track-setting method comprises the following steps:

the method comprises the steps of determining a track of external measurement data, and performing primary track calculation based on a Laplace method;

performing track determination on the GNSS internal measurement data, and performing primary track calculation by using a polynomial fitting method;

thirdly, forecasting ephemeris orbit determination, extracting specified epoch orbit data and performing initial orbit calculation.

Preferably, in the step S3, the specific step of track improvement includes:

modeling observation data;

selecting data and setting weights;

the dynamic model is selected and set;

and (5) selecting and setting the resolving parameters.

Preferably, in the step S3, the stability of the tracking calculation may be effectively improved by using a batch processing or sequential processing estimation method.

Preferably, the specific steps of the two estimation methods are as follows:

batch processing of

(1) Acquiring all observation data;

(2) constructing an observation equation;

(3) calculating covariance and residual error;

(4) solving a normal equation;

(5) an epoch optimal estimate;

sequentially processing of

(1) A current epoch state equation;

(2) state estimation;

(3) correcting the reference track;

(4) the next epoch state equation;

(5) and (5) optimally estimating.

Preferably, in the step S4, the dynamic model parameters are set in the ephemeris forecast calculation, and perturbation effects above the meter level need to be considered.

Preferably, in the step S5, the specific step of station tracking prediction includes:

forecasting the time of entering a station;

calculating the overhead time;

and 3, forecasting the outbound time.

Preferably, in the step S8, the specific step of ground shadow prediction includes:

calculating sun positions;

and (5) calculating shadow forecast.

The method has the beneficial effects that a high-precision satellite orbit determination and prediction method is realized, a mechanical model is given, the independent or combined use of external measurement and GNSS internal measurement data is supported, orbit improvement is carried out, the technical capability of orbit determination and prediction precision of a low-orbit spacecraft is broken through, team core competitiveness and independent intellectual property are formed, and service support and matched result conversion are provided for on-orbit management of a large number of satellites.

Drawings

FIG. 1 is a flow chart of a high-precision satellite orbit determination and prediction algorithm;

FIG. 2 is a schematic diagram of a track determination function interface of the present invention;

FIG. 3 is a block diagram of a track determination function of the present invention;

FIG. 4 is a flow chart of data preprocessing according to the present invention;

FIG. 5 is a flow chart of the initial track calculation process of the present invention;

FIG. 6 is a flow chart of the track improvement process of the present invention;

FIG. 7 is a flow chart of two estimation algorithms for track determination according to the present invention;

FIG. 8 is a schematic diagram of a track forecast function interface in accordance with the present invention;

FIG. 9 is a block diagram of a track forecast function in accordance with the present invention;

FIG. 10 is a flowchart of ephemeris forecast;

FIG. 11 is a flow chart of station tracking forecast in accordance with the present invention;

FIG. 12 is a flow chart of the device directed forecast of the present invention;

FIG. 13 is a schematic view of the undersea point of the present invention;

FIG. 14 is a flow chart of the present invention for the prediction of undersea points;

FIG. 15 is a flowchart of the present invention for moon image forecast;

fig. 16 is a plan view of a spacecraft in and out ground shadow of the present invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings.

The invention discloses a high-precision satellite orbit determination and prediction algorithm, which is shown in fig. 1 and comprises the following steps:

s1, S2 and S3 are the step contents of the orbit determination algorithm, process external measurement and GNSS internal measurement data, calculate the precise orbit of the satellite, and provide orbit numbers for injection parameters, orbit prediction and orbit control.

Precise orbit determination is a process of fitting a satellite orbit using precise dynamics through various kinds of observation data including observation errors. Obviously, the accuracy of orbit determination calculation is largely limited by the accuracy of measurement data, the accuracy of an observation model and the accuracy of a dynamics model, except for the observation geometry. The precise orbit determination relates to the processing of various observation data, including precise observation modeling and correction of various measurement errors, and the dynamic influence suffered by satellites of different orbit types is different, so that the dynamic models considered in orbit determination calculation are different. In addition, there is a certain error in both the dynamic model and the observation model, wherein the dynamic model is more remarkable, for example, the atmospheric damping force suffered by the low-orbit satellite is inevitably modeled by the atmospheric model. In summary, in the track calculation process, an appropriate mathematical model needs to be selected, and the measurement data is optimized, so that the accuracy of the calculation result is improved. The precision of satellite orbit determination calculation depends on a dynamic model and observation data, and the dynamic model to be considered for current precision orbit calculation comprises earth particle gravitation, earth non-spherical gravitation, solar optics, atmospheric damping force, post-Newton effect stress and the like.

Fig. 2 is a schematic diagram of an interface of the track determining function, and various information contents transmitted through the interface are shown in table 1:

table 1 track determination function interface information definition

Sequence number	Interface squareTo the direction of	Information type	Direction of transmission and reception	Delivery mode
					1.	User' s	Track determination command	Reception of	Operation interface
2.	Track database	Track measurement data	Reception of	Database access
					3.	Track database	Rail-fixing calculation result	Transmitting	Database access

The track determination function is shown in fig. 3.

the data preprocessing flow is shown in fig. 4, and when the system is started, a database interface is initialized; acquiring corresponding track measurement data from a database according to a track determination command sent by a user; executing a series of data processing processes such as data correction, outlier rejection, data format conversion and the like; and after finishing the data preprocessing, outputting a data preprocessing result.

The specific steps of data preprocessing include:

the method comprises the steps of correcting troposphere delay;

and based on temperature, humidity and pressure data measured at the measuring station, finishing troposphere refraction correction on satellite ranging and speed measurement through an atmospheric refraction correction model.

Ionosphere delay correction;

and finishing ionosphere delay correction on satellite ranging and speed measurement based on the ionosphere grid data.

Correcting the data time scale;

correcting the time marks of various measurement data.

Delay correction of equipment;

the transponder delay of the ground based ranging data is corrected.

Removing the wild value;

and performing outlier rejection on the external measurement data and the GNSS internal measurement data according to specific standards.

Converting the data format;

and performing format conversion after the track measurement data correction processing to generate data content in a format required by track determination calculation.

the satellite initial orbit information is generated by using the external measurement data, the internal measurement data or the forecast ephemeris file for orbit improvement.

The initial track calculation processing flow is shown in fig. 5, and classification processing is carried out according to the designated track determination method type when the system is started; respectively carrying out primary track calculation processes of different types, such as external measurement data track determination, GNSS internal measurement data track determination, forecast ephemeris track determination and the like; and generating and outputting an initial track calculation result.

For different data sources, the initial track is calculated according to the classification processing of the type of the track-setting method, and the method comprises the following steps:

based on the Laplace method, the satellite initial orbit calculation is completed by using ranging, azimuth and altitude. The dynamic model adopted in the primary orbit calculation comprises earth particle attraction, J2 perturbation and solar-lunar perturbation, and is suitable for low orbit and medium-high orbit satellites of the earth. The foundation ranging precision of the current measurement and control system is about 2m, the speed measurement precision is 0.01m/s, the precision of azimuth angle and altitude angle is 0.01 degrees, and the initial orbit precision which is better than 500m can be realized for medium-low orbit spacecrafts, so that the initial value requirement of follow-up orbit improvement is met.

and calculating the initial orbit information of the satellite by using satellite-borne GNSS self-positioning data and a polynomial fitting method.

For earth orbit satellites (orbit height is below 1000 km), the precision of the current satellite-borne GNSS self-positioning data is better than 12m, the speed precision is better than 0.05m/s, and the initial orbit precision depends on the ephemeris forecast precision.

Thirdly, forecasting ephemeris orbit determination, extracting specified epoch orbit data and performing initial orbit calculation;

and extracting satellite position and speed information of the appointed epoch from the forecast ephemeris file, or interpolating to calculate the position and speed information of the appointed epoch to generate satellite initial orbit information.

track improvement is based on a track dynamics model, and the initial track is iteratively improved by using measured data to generate a precise track. The differences between the different orbit types are mainly reflected in the different choice of the perturbation model. Corresponding power model usage and data usage settings need to be set in the configuration file for a particular type of track.

The orbit improvement is based on prior position and speed information provided by initial orbit calculation, an accurate dynamic model is generally adopted, and considered mechanical factors comprise earth particles, earth non-sphericity, sun-moon and large planet gravitation, atmospheric damping force, solar pressure and the like. For GNSS internal measurement data, since the self-positioning precision is better than 15m, the orbit improvement can realize orbit precision better than 12m through dynamic fitting; for ground measurement, the accuracy of measurement data and tracking arc segments are limited, and for a low-orbit satellite, the position accuracy of external measurement data adopting 4 circles of arc segments a day is better than 200 meters, and the speed accuracy is better than 0.05m/s; for medium-high orbit satellites, three stations are adopted for 24-hour time-division tracking, the position precision is better than 200 meters, and the speed precision is better than 0.05m/s.

The track improvement processing flow is shown in fig. 6, and a series of process initialization processes are executed when the system is started; performing iterative improvement based on the initial track calculation result and the measurement data; and after the iteration improvement convergence is completed, generating and outputting a precise orbit determination result.

The specific steps of track improvement comprise:

modeling observation data;

and carrying out accurate observation modeling on observation data such as external measurement, internal measurement, inter-satellite measurement and the like.

Selecting data and setting weights;

and setting a certain data type of parameter orbit determination calculation or a combination of several types of data types, an arc segment range of the observed data and weights used by various observed data through configuration files.

The dynamic model is selected and set;

according to the satellite orbit type, a specified dynamics model is set to be used, and optional dynamics models comprise: the order of non-spherical gravitation, solid tide, sea tide perturbation, atmospheric damping force, solar pressure, solar-lunar attraction and other large body perturbation. For satellites of different orbit types, different dynamics models can be used according to profile settings, the differences in dynamics models for low orbit and high orbit satellites are mainly manifested in atmospheric damping forces, orders of earth asphericity, and third body perturbation. Under the premise of not considering efficiency, the high-order non-spherical shape, the atmospheric damping force and the third body perturbation can be set to be used, and the track calculation requirements of the low track, the medium track and the high track can be met simultaneously.

Selecting and setting resolving parameters;

setting a resolving parameter according to the requirement, wherein the optional resolving parameter comprises: (1) position, speed; (2) atmospheric damping force; (3) solar light pressure coefficient; (4) constant empirical force in RTN direction; (5) RTN direction is periodically empirical.

The orbit determination is based on the principle that the best estimate of the satellite state is obtained by using measurement data containing measurement errors in combination with the dynamic constraints imposed on the satellite. For precise orbit determination, the traditional method is called orbit improvement, but because some dynamics related to orbit or physical parameters on measurement can be determined while the satellite orbits, the simple orbit improvement in the traditional sense is expanded, which is now called precise orbit determination, and the stability of orbit determination calculation can be effectively improved by adopting a proper estimation method.

The stability of the orbit determination calculation is effectively improved by adopting two estimation methods of batch processing or sequential processing, and the specific steps of the two estimation methods are as follows:

batch processing of

The batch processing utilizes all the observation data to obtain the track parameters at a certain moment, and is suitable for post-processing. The method is characterized in that more data and high solving precision are obtained, but the storage amount is large, all observation data are required to be stored, and after the data sampling is finished, the calculation can be carried out, and the method comprises the following steps:

(1) acquiring all observation data;

(2) constructing an observation equation;

(3) calculating covariance and residual error;

(4) solving a normal equation;

(5) an epoch optimal estimate;

sequentially processing of

The sequential processing is to calculate the satellite orbit at the moment as soon as the observed data exist, and the problem of large storage capacity is avoided, so that the method is widely applied to quick and real-time orbit determination. The main problem is that recursion diverges easily, and as time goes by, the calculated trajectory gradually deviates from the true trajectory, a stress model needs to be introduced to compensate algorithm errors. The method comprises the following steps:

(1) a current epoch state equation;

(2) state estimation;

(3) correcting the reference track;

(4) the next epoch state equation;

(5) an optimal estimation;

the specific process flow of the two estimation methods is shown in fig. 7.

In general, the state differential equation for satellite motion can be written as follows:

wherein X represents a state quantity of the satellite

And other parameters to be estimated->

To give X in the above formula in a given reference state

The place is unfolded and recorded->

The higher order term may be omitted and represented as a linear equation,

wherein, the liquid crystal display device comprises a liquid crystal display device,

its solution can be expressed as

x(t)＝Φ(t ₀ ,t)x(t ₀ )

The above formula of phi (t) ₀ ,t ₀ ) Referred to as the state transition matrix, I represents the identity matrix.

The observed quantity has certain relation with the satellite state quantity and can be expressed as

Y＝G(t,X)+ε

Where Y represents the observed quantity and ε represents the measurement noise.

Also, expanding on the above at the reference state, omitting higher order terms has

In the method, in the process of the invention,

represents the observed partial derivative of the observed quantity with respect to the state quantity at the observation epoch, and H represents the observed partial derivative with respect to the state quantity of the improvement epoch. For time series observations, a formula can be expressed as

y _i ＝H _i x ₀ +ε _i

Estimation of the state equation is typically accomplished in two ways: batch processing and sequential processing.

The batch processing is to calculate all the observed data to be processed at one time, and generally takes the weight matrix W of the observation to satisfy:

omitting the deriving process and estimating an estimate of the available batch based on the linear unbiased minimum variance

Is that

Wherein R is _i To observe the square of the noise level of the data,

is the a priori covariance of the state quantity to be estimated.

Thereby calculating an optimal estimate of the epoch to be estimated

The method comprises the following steps:

the corresponding covariance matrix is:

sequential processing is a recursive algorithm. Let t be known _k Estimation at time (k=0, 1,2, …, n)

And covariance matrix->

Then t _k+1 The state quantity at the moment is:

the formulas constitute a recursive formula for sequential processing. t is t _k+1 The time-of-day corrected state estimate may be expressed as,

s4 to S8 are the step contents of the track prediction algorithm, the track prediction is performed by using a precise track to perform extrapolation calculation ephemeris prediction, and then related prediction calculations such as measurement and control equipment tracking condition prediction (three-point prediction: inbound, overhead and outbound), equipment guiding prediction, satellite point prediction, ground moon shadow prediction and the like are performed on the basis of the precise track.

The prediction precision of the low orbit spacecraft is limited by the calculation of atmospheric damping force, the existing atmospheric model has 10-15% model error, the atmospheric density error during the magnetic storm can be several times, and for the low orbit satellite (orbit height is less than 1000 km), the prediction precision of 12 hours is better than 100m, and the orbit prediction precision of the medium-high orbit satellite for 7 days is better than 200 m.

As previously described, orbit predictions include ephemeris predictions, station tracking predictions, equipment guiding predictions, satellite spot predictions, ground moon predictions, etc. The invention has the advantages that the codes of the track forecasting part are all independently researched and developed, the C++ language is adopted for coding, the reasonable data structure is designed, the object-oriented programming thought is adopted, the algorithm flow is re-carded and optimized, the storage and access of variables are especially optimized, the matrix and vector operation increases the basic functions covering the linear algebra, a good operation basis is provided for processing differential equations and numerical calculation, the double-precision floating point operation is adopted for data storage, and the efficiency is higher and better. Track prediction is performed for 24 hours, and the calculation time is better than 30 seconds; the single-star multi-station 8-arc segment data orbit determination calculation is better than 1 minute.

Fig. 8 is a schematic diagram of an interface of the track prediction function, and various information contents transmitted through the interface are shown in table 2:

table 2 track forecast function interface information definition

Sequence number	Interface orientation	Information type	Direction of hair	Delivery mode
					1.	User' s	Track forecast command	Reception of	Operation interface
2.	Track database	Satellite orbit number	Reception of	Database access
					3.	Track database	Forecast calculation results	Transmitting	Database access

The track forecast function is shown in fig. 9.

the orbit root number supports the Kepler six root number and two rows of orbit root numbers, wherein the calculation model of the Kepler six root number is an HPOP high-precision orbit prediction model. The forecasting model of the number of the two lines of tracks adopts the SGP4/SDP4 model of NOARD, which are commonly adopted in the industry at present, and the accuracy and the efficiency of the model are verified by a large number of actual engineering tests, so that the model can well support the bottom algorithm service of the testing operation control service.

The ephemeris forecast process flow is shown in figure 10, and when the system is started, the database connection is initialized; acquiring the corresponding satellite orbit number from a database according to an orbit forecast command sent by a user; setting dynamic model parameters; calculating satellite ephemeris forecast; after the forecast calculation is completed, storing an ephemeris forecast result;

for the motion law of satellites, taking the position vector r and the velocity vector into consideration under the geocentric celestial coordinate system

To describe.

The motion equation is the initial value problem P

Right function a can be divided into a0 and a _ε Two parts, a0 is the acceleration of the spacecraft under the two-body problem, a _ε For the sum of other various dynamic accelerations, it can be expressed as follows

Considering the perturbation effects above the meter scale, satellites are affected including: non-spherical gravitational force of the earth; attraction between sun and moon; solar pressure perturbation; atmospheric resistance perturbation; solid tide perturbation, etc.

Non-spherical gravitational perturbation of the earth

In the geodetic coordinate system, the bit function of the global non-spherical gravitational perturbation is:

wherein:

and->

Normalized Legendre and associated Legendre functions, respectively

Harmonic term coefficients for normalized earth gravitational field

And->

Tian Xiexiang coefficients for normalized earth gravitational field

GM _E Is the constant of the gravitational force

R _E For the average equatorial radius of the earth

Is longitude and latitude

r is the distance between the centers of earth

The normalization process is as follows:

the recurrence relation of (2) is as follows:

the recurrence relation of (2) is as follows:

/>

the global non-spherical perturbation acceleration can thus be obtained as follows:

wherein:

wherein:

the partial derivative of the legendre polynomial is calculated using the following derivation formula:

/>

the associated Legendre polynomial partial derivative recurrence formula in the fan harmonic term and the field harmonic term is as follows:

the acceleration of the spacecraft in the earth fixed coordinate system is converted into a corresponding epoch geocentric celestial coordinate system.

N-shaped gravitational perturbation

The sun, moon and large planets of the solar system all have attractive effects on the aircraft, the effect of which can be described approximately as point mass perturbation P. In the epoch geocentric celestial coordinate system, the following formula can be used to calculate N-body perturbation:

wherein: GMBjB is the gravitational constant of the jth perturbation; r is R _j The position vector of the jth perturbation body in the geocentric inertial system; delta _j For spaceflightThe director of the jth perturbation. Calendar generated by the us jet propulsion laboratory (Jet Propulsion Laboratory, JPL) provides accurate values for sun, moon and large planet positions.

Solar radiation pressure perturbation

Solar radiation pressure perturbation is P due to solar radiation impinging on the spacecraft creating pressure on the spacecraft. Also referred to as light pressure perturbation. Thus, there is only a perturbation in light pressure when the spacecraft position is outside the shadows of the earth and moon.

The solar radiation perturbation acceleration of the body portion is calculated as follows:

wherein:

f: shadow factors, which can be calculated from the position of the sun, month, earth and spacecraft and cylindrical or conical ground shadow models, divide the shadows of the earth and moon into penumbra, penumbra and pseudo penumbra.

ρ _SR : solar pressure constant in the vicinity of the spacecraft.

η: reflection coefficient of the illuminated surface of the spacecraft. The value range of eta is as follows: 0 < η < 1, typically η=0.5 is preferred. η may be used as a correction factor for solar radiation pressure perturbation, estimated in the orbit estimation process.

m: spacecraft mass.

Δ _S : spacecraft to sun vector.

A: perpendicular to

Is a cross-sectional area of the spacecraft.

The different angles of installation of their solar panels and the manner in which the orientation of the sun is controlled may be different. The direction of solar radiation against the solar panel pressure is also different. Thus, solar radiation perturbation accelerations are calculated according to the specifics of each spacecraft solar panel. Calculating the perturbation acceleration of the solar sailboard:

wherein:

A _P : area of solar sailboard.

Beta: the reflection coefficient of solar sailboard to solar radiation can be used as estimated value.

Atmospheric resistance perturbation

The low orbit spacecraft flying around the earth is in an environment other than vacuum, and the atmosphere generates resistance P to the spacecraft. Atmospheric drag perturbation may be described as:

wherein: CBDB is the atmospheric drag coefficient; a is the projection of the cross-sectional area of the aircraft on a plane perpendicular to the orbit; m is the mass of the aircraft; ρbab is the atmospheric density at the aircraft; v (V) _r Is the velocity of the aircraft relative to the atmosphere.

For a spacecraft with a solar sailboard, the atmospheric resistance generated by the solar sailboard needs to be considered, and the formula is similar to the formula above:

wherein: ABpB is the area of the solar sailboard; the CBDPB is an atmospheric resistance coefficient suitable for a solar sailboard; gamma is the angle between the normal direction of the solar sailboard and the speed of the aircraft; i A _P cos γ is the effective area of the solar panel in a plane perpendicular to the track.

The above-mentioned atmospheric density ρ _a The calculation may be performed using atmospheric density models, including: exponential models, harris-Priester, jacchia 77, dtm (Drag Temperature Model) and MSIS series, etc., see the relevant data.

The disturbance of the solid tide

Since the earth is not a rigid body, the land portion of the earth is elastically deformed P by the lunar and solar gravitation sites. This deformation is known as solid tide. The solid tide makes the fluctuation amplitude of the earth shell reach 20-30 cm. The solid tide causes the mass distribution within the earth to change over time, and thus the gravitational field of the earth to also change over time, an effect known as dynamic tides. This change in the gravitational attraction caused by the solid tide is called the solid tide additional bit. Typically, the level of the induced moisture is represented by a harmonic function. The LOVE numbers and the load deformation coefficient values corresponding to the tide level of different orders are different. Only the second order moisture guiding level is considered for solid tide at present. The change in the earth's attraction caused by the solid tide can be expressed as a change in the spherical harmonic coefficient of the earth's attraction:

wherein:

the change of the gravitational potential coefficient caused by the gravitational attraction of the sun and the moon; HBkB, Θ _k Amplitude and phase of k divided tidal wave; />

The k-order Love number; delta _0m Is a Kronecker operator; a, a _e Is the average equatorial radius of the earth.

The calculation of the perturbation acceleration caused by the solid tide additional bit is similar to the calculation of the global non-spherical gravitational perturbation acceleration, and will not be described here.

Influence of relativistic theory

Due to the generalized relativistic effect, the equations of motion of an aircraft in a local inertial coordinate system with the earth centroid as the origin will differ from equation of motion P when only newton gravitational fields are considered. This difference can be seen as an additional perturbation to the aircraft. Generalized relativistic perturbation of a near-earth aircraft may be described as

Wherein: c is the speed of light; r is (r) _E ,

Position and velocity vectors for the aircraft in the geocentric inertial system; GM (GM) _e Is the gravitational constant of the earth; beta, gamma are relativistic parameters, and the value of the relativistic parameters is 1 for Einstein.

based on ephemeris forecast results and combined with geographical coordinates of measurement and control station equipment, three-point forecast, namely station tracking forecast, is calculated, wherein forecast contents comprise station measurement time, azimuth angle, altitude angle, distance and distance change rate of three time points of station entering, station exiting, station passing and the like.

The flow of station tracking and forecasting is shown in figure 11, and when the system is started, the database connection is initialized; acquiring a corresponding satellite ephemeris forecast from a database according to an orbit forecast command sent by a user; sequentially calculating the overhead time of each angle; calculating other parameters such as the distance corresponding to the overhead time; and after the forecast calculation is completed, the station tracking forecast result is stored.

The specific steps of station tracking forecast include:

forecasting the time of entering a station;

and calculating satellite arrival time and azimuth according to the tracking height cut-off angle, and calculating distance and distance variability.

Calculating the overhead time;

and calculating the satellite overhead time, calculating the altitude angle and the azimuth angle, and calculating the distance and the distance variability.

Forecasting the outbound time;

based on ephemeris forecast results and combined with the geographic coordinates of the measurement and control station equipment, the guiding angle of the antenna tracking satellite of the measurement and control station equipment is calculated. A coordinate system (ENU, northeast coordinate system) is usually established on the local ground level where the antenna is located, and is used for describing the motion state of the satellite in the antenna coordinate system, and the description of the state can be represented by three basic physical quantities, namely an azimuth angle (Az, a value range of 0-360 °), a pitch angle (El, a value range of 0-90 °), and a distance (R).

The processing flow of the device guiding forecast is shown in fig. 12, and when the system is started, the database connection is initialized; acquiring a corresponding satellite ephemeris forecast from a database according to an orbit forecast command sent by a user; in the equipment tracking time period, calculating the azimuth angle and the pitch angle of the corresponding equipment antenna at fixed time intervals; after the forecast calculation is completed, the storage device guides the forecast result.

based on the ephemeris forecast, the satellite lower point position is calculated. The position of the satellite point can reflect the running track of the satellite relative to the earth surface, and provides decision basis for the satellite to develop application service, such as coverage condition and revisit period of a certain area. The satellite lower point can be characterized by three basic physical quantities, namely longitude (Lon, value range-180 degrees, east longitude is positive), latitude (Lat, value range-90, north latitude is positive) and altitude (Alt). The satellite under-satellite points are connected into a whole in a certain period of time, so that a track, namely an under-satellite point track, is formed on the ground. Assuming that the earth is a rotational ellipsoid, the normal line of the earth is taken as the center of mass of the satellite, the intersection point of the normal line and the ellipsoid is defined as the undersea point, i.e. the vertical projection of the satellite on the ellipsoid is defined as the undersea point, as shown in fig. 13.

The processing flow of the satellite lower point forecast is shown in figure 14, and when the system is started, the database connection is initialized; acquiring a corresponding satellite ephemeris forecast from a database according to an orbit forecast command sent by a user; calculating the corresponding position of the satellite point at fixed time intervals in a specified forecast time period; and after the forecast calculation is completed, storing the forecast result of the point under the satellite.

S8, forecasting the earth and moon shadow, and calculating the time period of the satellite entering the earth or moon shadow area based on the ephemeris forecasting result;

based on the ephemeris forecast, a time period for the satellite to enter the earth or moon shadow area is calculated. Certain satellite loads normally work under specific illumination conditions, so that the relative position relation of the satellite relative to the sun and moon is required to be considered, the forecasting results of the shadow surface and the sun surface are given, and decision basis is provided for the arrangement of a mission plan.

The processing flow of the ground moon shadow forecast is shown in figure 15, and when the system is started, the database connection is initialized; acquiring a corresponding satellite ephemeris forecast from a database according to an orbit forecast command sent by a user; calculating the position of the sun in a specified forecast time period; respectively calculating the position of the earth and the moon in a specified forecast time period; calculating the time when the satellite enters and exits the earth or moon shadow area in a specified forecast time period; and after the forecast calculation is completed, storing the ground moon shadow forecast result.

The specific steps of the ground moon shadow forecast include:

calculating sun positions;

the ephemeris is read to calculate the sun position.

Shadow predictors;

and calculating the time for the satellite to enter and exit the ground shadow according to the sun position.

Ground moon shadow calculation principle: since the sun is far away, the sun is at infinity (i.e., neglecting parallax), the solar rays on earth are parallel light, as can be seen from FIG. 16, the angle between the station position vector R and the sun direction L2 (unit vector) is greater than 90+|beta ₀ I, i.e. in the second quadrant, this condition is as follows:

on the ground, the ground consists of 90 DEG++ |beta ₀ The boundaries defined by the are called solar lines, which divide the earth's surface into two areas, namely "daytime" and "nighttime", optical observations can be made when the station is at night.

The optical observation also requires that the surface of the spacecraft is illuminated by the sun, and the device can observe the reflected light rays, which is determined by the relative position relationship between the spacecraft and the sun, namely that the spacecraft must be outside the ground shadow, and for this condition, only a simple cylindrical ground shadow needs to be considered, as shown in fig. 16.

As can be seen from fig. 16, the following conditions are satisfied on the ground shadow boundary:

the right end of cos psi is given a negative sign, because psi is more than 90 DEG, when the included angle between the spacecraft and the sun direction is less than psi, the spacecraft is not in the ground shadow and can be illuminated by the sun, the intersection point of the connecting line of the sun and the earth center and the ellipsoid is called a sunward point, the area with the sunward point as the center radius psi on the ground is the area where the spacecraft can be illuminated, the boundary of the area is called a ground shadow line, and the size of the area without the ground shadow can be estimated by replacing r with the radius of the track half length, namely

For optical observations of a spacecraft, all three conditions mentioned above must be met, the observable range being constant with respect to the measuring station for one spacecraft, while the line of sight and the line of earth are moved from east to west as a function of the position of the sun, for other observations, such as laser observations, radio observations, etc., these latter two conditions being theoretically not necessary.

Claims

1. The high-precision satellite orbit determination and prediction algorithm is characterized by comprising the following steps of:

2. The high-precision satellite orbit determination and prediction algorithm according to claim 1, wherein in S1, the specific step of preprocessing the data comprises:

the method comprises the steps of correcting troposphere delay;

ionosphere delay correction;

correcting the data time scale;

delay correction of equipment;

removing the wild value;

and (5) converting the data format.

3. The high-precision satellite orbit determination and prediction algorithm according to claim 1, wherein in S2, for different data sources, the initial orbit is calculated according to the orbit determination method type classification process.

4. A high accuracy satellite orbit determination and prediction algorithm according to claim 3, in which the method of calculating the primary orbit is based on an orbit determination method type classification process:

5. The high-precision satellite orbit determination and prediction algorithm according to claim 1, wherein in S3, the specific step of orbit improvement comprises:

modeling observation data;

selecting data and setting weights;

the dynamic model is selected and set;

and (5) selecting and setting the resolving parameters.

6. The high-precision satellite orbit determination and prediction algorithm according to claim 1, wherein in S3, the stability of orbit determination calculation can be effectively improved by using a batch processing or sequential processing estimation method.

7. The high-accuracy satellite orbit determination and prediction algorithm according to claim 6, wherein the specific steps of the two estimation methods are as follows:

batch processing of

(1) Acquiring all observation data;

(2) constructing an observation equation;

(3) calculating covariance and residual error;

(4) solving a normal equation;

(5) an epoch optimal estimate;

sequentially processing of

(1) A current epoch state equation;

(2) state estimation;

(3) correcting the reference track;

(4) the next epoch state equation;

(5) and (5) optimally estimating.

8. The high-precision satellite orbit determination and prediction algorithm according to claim 1, wherein in S4, dynamic model parameters are set in ephemeris prediction calculation, and perturbation effects above the meter level need to be considered.

9. The high-precision satellite orbit determination and prediction algorithm according to claim 1, wherein in S5, the specific step of station tracking prediction comprises:

forecasting the time of entering a station;

calculating the overhead time;

and 3, forecasting the outbound time.

10. The high-precision satellite orbit determination and prediction algorithm according to claim 1, wherein in S8, the specific step of ground shadow prediction comprises:

calculating sun positions;

and (5) calculating shadow forecast.