Unscented Batch Orbit Determination Method and System for Continuously Small-Thrust Maneuvering Satellites

Abstract

The present invention belongs to the field of aerospace technology, and relates to a method and system for orbit determination of a continuously small-thrust maneuvering satellite by unscented batch processing. The present invention obtains a radar observation vector according to radar measured data or satellite ephemeris data in combination with the radar station location coordinates; obtains the initial orbit state at the moment to be estimated according to the radar observation vector and by using the Lambert algorithm; generates an initial augmented state quantity according to the initial orbit state at the moment to be estimated; constructs a simplified orbit dynamics equation including the non-spherical perturbation of the earth, the atmospheric drag perturbation and the constant small-thrust acceleration in the RTN coordinate system in the ECI coordinate system; based on the simplified orbit dynamics equation, obtains the precise value of the augmented state quantity at the moment to be estimated through the initial augmented state quantity, the radar observation vector and the unscented batch processing filtering method; and substitutes the precise value of the augmented state quantity at the moment to be estimated into the simplified orbit dynamics equation for orbit prediction at other moments. The present invention significantly improves the orbit determination accuracy and reduces the influence of single-point observation errors.

Description

Technical Field

[0001] The present invention belongs to the technical field of aerospace, and relates to a method and system for orbit determination of a continuous small-thrust maneuvering satellite by unscented batch processing. Background Technique

[0002] In recent years, with the rapid development of low-Earth orbit (LEO) mega-constellations represented by "Starlink", the number of LEO space targets has increased rapidly, making the LEO region more crowded and posing a serious threat to the safety of space assets. Due to the large number of such satellites and their frequent maneuvers, higher requirements are also placed on space situational awareness. After the LEO mega-constellation satellites are launched into orbit, they usually use on-board electric thrusters to gradually raise the operating orbit to the target orbit in stages, and the satellite target group gradually spreads out under the action of the electric thrusters. Since the electric thrust is small and the maneuvering laws of non-cooperative targets are unknown, the maneuvering time in the orbit-raising and orbit-lowering segments may be as long as several months. Traditional non-thrust orbit determination and orbit prediction methods will bring large errors and cannot support effective space target cataloging and collision warning, etc.

[0003] In satellite orbit determination, the most commonly used estimation algorithms include batch processing methods (such as the least squares method) and sequential processing methods (such as Kalman filtering, etc.).

[0004] The main feature of the batch least squares method is the linearization approximation of the non-linear equations. The non-linear system and measurement equations are linearized and approximated through Taylor series expansion, and the higher-order terms are excluded. Then, the partial derivatives of the linearized equations including the system and measurement equations with respect to the state also need to be calculated. This method uses all (or most) of the observation data collected within a given observation period for one-time processing and solving to obtain the orbit parameters. The calculation process is relatively stable, and it has a good effect on orbit determination of uncontrolled (or with few maneuvers) space targets, and relatively stable and high-precision orbit determination results can be obtained. However, when the complexity of the non-linear equations and the number of state variables increase, the calculation of the partial derivatives becomes very complicated. When the system has strong non-linearity, the amount of observation data is small or sparse, the batch least squares method is prone to non-convergence, resulting in orbit determination failure.

[0005] The unscented Kalman filter uses the unscented transform and utilizes a set of sampling points ( SigmaThe Gaussian distribution of a random variable is described by points, and then through the transmission of a non - linear function, the weighted statistical linear regression technique is used to approximate the posterior mean and variance of the non - linear function. This process does not require linearizing the equations and has good adaptability to systems with strong non - linearity. However, its filtering effect depends to a large extent on initial estimation assumptions such as the prior variance and process noise. In addition, this method updates the system state variables in real - time through the state equation and the observation equation, and combines the newly obtained measurement information to correct the system state variables, thereby recursively obtaining better estimation results. However, for space targets with continuous small - thrust maneuvers, when the time interval between adjacent visible arc segments in the observation data is relatively long (possibly up to more than a dozen hours), due to the unknown small - thrust acceleration, the state covariance updated using traditional formulas often spreads, easily leading to slower filtering convergence or even divergence, resulting in orbit determination failure.

[0006] In summary, traditional non - thrust orbit determination and orbit prediction methods will bring large errors and cannot support effective space target cataloging and collision warning, etc. Summary of the Invention

[0007] The purpose of the present invention is to provide an unscented batch orbit determination method, system equipment, and medium for a satellite with continuous small - thrust maneuvers to solve the technical problem of large errors in traditional non - thrust orbit determination and orbit prediction methods.

[0008] To achieve the above - mentioned purpose, the present invention adopts the following technical solutions:

[0009] In the first aspect, the present invention provides an unscented batch orbit determination method for a satellite with continuous small - thrust maneuvers, including the following steps:

[0010] When there are radar measured data and radar station coordinates, the radar observation vector is obtained according to the radar measured data combined with the radar station coordinates. When there are no measured data, the pseudo - radar observation vector is obtained based on the satellite ephemeris data;

[0011] The initial orbit state at the moment to be estimated is obtained according to the radar observation vector and by using the Lambert algorithm;

[0012] The initial extended state quantity is generated according to the initial orbit state at the moment to be estimated;

[0013] A simplified orbit dynamics equation including the non - spherical perturbation of the earth, the atmospheric drag perturbation, and the constant small - thrust acceleration in the RTN coordinate system is constructed in the ECI coordinate system;

[0014] Based on the simplified orbit dynamics equation, the precise value of the extended state quantity at the moment to be estimated is obtained through the initial extended state quantity, the radar observation vector, and the unscented batch filtering method;

[0015] Substitute the precise value of the extended state quantity at the moment to be estimated into the simplified orbit dynamics equation for orbit prediction at other moments.

[0016] In a second aspect, the present invention provides an unscented batch orbit determination system for a continuously small-thrust maneuvering satellite, including:

[0017] A radar observation vector acquisition module: used to obtain a radar observation vector according to the radar measured data in combination with the radar station coordinates when there are radar measured data and radar station coordinates, and obtain a pseudo-radar observation vector based on the satellite ephemeris data when there is no measured data;

[0018] An initial orbit state acquisition module: used to obtain the initial orbit state at the moment to be estimated according to the radar observation vector and using the Lambert algorithm;

[0019] An extended state quantity generation module: used to generate an initial extended state quantity according to the initial orbit state at the moment to be estimated;

[0020] An orbit dynamics equation construction module: used to construct a simplified orbit dynamics equation including the non-spherical perturbation of the earth, the atmospheric drag perturbation, and the constant small-thrust acceleration in the RTN coordinate system in the ECI coordinate system;

[0021] An extended state quantity precise value acquisition module: used to obtain the precise value of the extended state quantity at the moment to be estimated based on the simplified orbit dynamics equation, through the initial extended state quantity, the radar observation vector, and the unscented batch filtering method;

[0022] An orbit prediction module: used to substitute the precise value of the extended state quantity at the moment to be estimated into the simplified orbit dynamics equation for orbit prediction at other moments.

[0023] In a third aspect, the present invention provides an electronic device, including: a processor; a memory for storing computer program instructions; and for implementing the unscented batch orbit determination method for a continuously small-thrust maneuvering satellite when executing the computer program.

[0024] In a fourth aspect, the present invention provides a storage medium, the storage medium stores computer program instructions, and when the computer program instructions are loaded and run by a processor, the processor executes the unscented batch orbit determination method for a continuously small-thrust maneuvering satellite.

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

[0026] The method of the present invention obtains a radar observation vector according to the radar measured data combined with the radar site coordinates when there are radar measured data and radar site coordinates. When there are no measured data, a pseudo radar observation vector is obtained according to the satellite ephemeris data to provide a basis for subsequent orbit determination; the initial orbit state at the time to be estimated is obtained according to the radar observation vector and the Lambert algorithm is used to quickly provide an initial orbit estimate. An extended state quantity is generated according to the orbit state at the time to be estimated, so as to more comprehensively describe the motion state of the satellite. A simplified orbital dynamics equation containing the non-spherical perturbation of the earth, the atmospheric drag perturbation and the constant small thrust acceleration in the RTN coordinate system is constructed in the ECI coordinate system; based on the simplified orbital dynamics equation, the precise value of the extended state quantity at the time to be estimated is obtained by the initial extended state quantity, the radar observation vector and the traceless batch filtering method. By integrating the traceless transformation and the batch processing method, it is not necessary to calculate the partial derivatives of the state quantity of the system and the measurement equation, which reduces the complexity of the formula derivation and improves the adaptability to the nonlinear system. The precise value of the extended state quantity at the time to be estimated is brought into the simplified orbital dynamics equation to perform orbit prediction at other times. The present invention utilizes all the observation data collected in an observation period, processes them at one time, and obtains improved orbital parameters through iterative solution, thereby reducing the dependence on the initial state covariance and process noise estimation, and improving the stability of the calculation process. The orbit determination accuracy is significantly improved by combining untraceable transformation with batch filtering, and the influence of single-point observation errors is reduced by batch processing of observation data. The present invention can accurately estimate the equivalent continuous constant small thrust acceleration in the local coordinate system when the thrust acceleration value is unknown and the observation data is sparse, and effectively perform high-precision orbit determination and prediction for low-orbit small-thrust continuous maneuvering satellites.

[0027] The calculation process of the present invention does not involve the updating of the state covariance, thus avoiding the influence caused by the possible diffusion of the state covariance.

[0028] The system of the present invention includes a radar observation vector acquisition module, an initial orbit state acquisition module, an extended state variable generation module, an orbit dynamics equation construction module, an extended state variable precise value acquisition module, and an orbit prediction module; the radar observation vector acquisition module is used to obtain the radar observation vector according to the radar measured data in combination with the radar site coordinates when there are radar measured data and radar site coordinates, and obtain the pseudo-radar observation vector based on the satellite ephemeris data when there are no measured data; the initial orbit state acquisition module is used to obtain the initial orbit state at the moment to be estimated according to the radar observation vector and by using the Lambert algorithm; the extended state variable generation module is used to generate the initial extended state variables according to the initial orbit state at the moment to be estimated; the orbit dynamics equation construction module is used to construct a simplified orbit dynamics equation including the non-spherical perturbation of the earth, the atmospheric drag perturbation, and the constant small thrust acceleration in the RTN coordinate system in the ECI coordinate system; the extended state variable precise value acquisition module is used to obtain the precise value of the extended state variables at the moment to be estimated based on the simplified orbit dynamics equation, through the initial extended state variables, the radar observation vector, and the unscented batch filtering method; the orbit prediction module is used to substitute the precise value of the extended state variables at the moment to be estimated into the simplified orbit dynamics equation for orbit prediction at other moments. The system of the present invention can significantly improve the orbit determination accuracy and reduce the influence of single-point observation errors.

[0029] The device and medium of the present invention can also significantly improve the orbit determination accuracy and reduce the influence of single-point observation errors. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 It is the overall flowchart of continuous small thrust maneuver satellite orbit determination based on unscented batch filtering in an embodiment of the present invention;

[0031] Figure 2 It is the flowchart of the unscented batch filtering algorithm in an embodiment of the present invention;

[0032] Figure 3 It is the satellite position error diagram in the RTN coordinate system in an embodiment of the present invention;

[0033] Figure 4 It is the method flowchart in an embodiment of the present invention;

[0034] Figure 5 It is the system module diagram in an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0035] To enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the scope of protection of the present invention.

[0036] It should be noted that the terms "first", "second", etc. in the specification of the present invention and the above-mentioned drawings are used to distinguish similar objects, and do not necessarily need to be used to describe a specific order or sequence. It should be understood that such data can be interchanged under appropriate circumstances so that the embodiments of the present invention described here can be implemented in an order other than those illustrated or described here. In addition, the terms "comprising" and "having" and any variations thereof are intended to cover non-exclusive inclusion. For example, a process, method, system, product or device comprising a series of steps or units does not necessarily have to be limited to those steps or units clearly listed, but may include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0037] The present invention will be further described in detail below in conjunction with the accompanying drawings:

[0038] Embodiment 1:

[0039] See Figure 4 , this embodiment discloses an unscented batch orbit determination method for a continuously small-thrust maneuvering satellite, including the following steps:

[0040] S1. When there are radar measured data and radar station coordinates, obtain the radar observation vector according to the radar measured data in combination with the radar station coordinates. When there are no measured data, obtain the pseudo-radar observation vector based on the satellite ephemeris data to provide a basis for subsequent orbit determination, specifically as follows:

[0041] If there are radar measured data and radar station coordinates, extract the radar observation vector from the radar measured data according to the radar measured data and radar station coordinates. The radar observation vector includes the slant range, azimuth angle, and elevation angle of the satellite;

[0042] When there are no measured data, obtain the pseudo-radar observation vector based on the satellite ephemeris data, specifically as follows:

[0043] Given the radar station coordinates, extrapolate the position and velocity of the satellite based on the given radar station coordinates using the two-line orbital element data of the satellite, or obtain the position and velocity of the satellite through the precise orbit data file;

[0044] Calculate the positions and velocities of the satellite and the observation station in the geocentric coordinate system based on the position and velocity of the satellite, and calculate the position of the satellite relative to the observation station.

[0045] Transform the satellite position and velocity in the geocentric coordinate system and the position of the satellite relative to the observation station into the local north-east-down coordinate system of the observation station, and calculate the theoretical value of the observation vector.

[0046] Add random observation errors to the three parameters of slant range, azimuth angle, and elevation angle in the theoretical value of the observation vector, and impose detection visibility constraints on the slant range and elevation angle to simulate and generate a pseudo-radar observation vector.

[0047] S2. Obtain the initial orbit state at the moment to be estimated according to the radar observation vector and using the Lambert algorithm, as follows:

[0048] There are at least three radar observation vectors in the observation arc segment where the moment to be estimated is located.

[0049] When the radar observation vector at the moment to be estimated is the first or the last radar observation vector in the observation arc segment, obtain the position vectors of the satellite at the first and the last radar observation vectors in the observation arc segment respectively, and use the Lambert algorithm to obtain the orbit state of the satellite at the moment to be estimated.

[0050] When the radar observation vector at the moment to be estimated is between the first and the last radar observation vectors in the observation arc segment, use the first or the last radar observation vector in the observation arc segment, combine it with the radar observation vector at the moment to be estimated, and use the Lambert algorithm to obtain the orbit state of the satellite at the moment to be estimated.

[0051] The orbit state at the moment to be estimated includes the initial position and velocity vector of the satellite at the moment to be estimated.

[0052] It should be noted that the Lambert algorithm is a classic method for solving orbit transfer in the two-body problem, and its core is to calculate the transfer orbit of the spacecraft between two points by given initial and target positions and the transfer time.

[0053] S3. Generate the initial extended state quantity according to the initial orbit state at the moment to be estimated, as follows:

[0054] Add a constant small thrust acceleration vector in the RTN coordinate system to the orbit state at the moment to be estimated to obtain the extended state quantity.

[0055] The initial value of the small thrust acceleration vector is set to zero.

[0056] S4. Construct a simplified orbital dynamics equation in the ECI coordinate system that includes the Earth's non-spherical perturbation, atmospheric drag perturbation, and the constant small thrust acceleration in the RTN coordinate system, as follows:

[0057] The extended state variables include the position and velocity vectors of the satellite in the ECI coordinate system, and the constant small thrust acceleration vector of the satellite in the orbital RTN coordinate system. Among them, the thrust acceleration vector of the satellite in the orbital RTN coordinate system includes the radial component, along-track component, and normal component of the thrust acceleration in the RTN coordinate system;

[0058] It should be noted that the ECI coordinate system (Earth-Centered Inertial Frame) is a commonly used reference coordinate system for describing the motion and position of spacecraft, satellites, and other objects.

[0059] The specific form of the simplified orbital dynamics equation is as follows:

[0060]

[0061] In the formula, is the position vector of the satellite in the ECI coordinate system, are the velocity vectors of the satellite in the ECI coordinate system respectively, is the first derivative of the position vector of the satellite with respect to time, is the first derivative of the velocity vector of the satellite with respect to time, is the gravitational acceleration at the Earth's center, is the gravitational constant of the Earth; is the acceleration caused by the Earth's non-spherical perturbation; is the atmospheric drag perturbation acceleration, is the radial component of the thrust acceleration in the RTN coordinate system, is the along-track component of the thrust acceleration in the RTN coordinate system, is the normal component of the thrust acceleration in the RTN coordinate system. The thrust acceleration vector of the satellite in the orbital RTN coordinate system includes , and , is the first derivative with respect to time, is the first derivative with respect to time, is the first derivative with respect to time, is the coordinate transformation matrix from the RTN coordinate system to the ECI coordinate system.

[0062] S5. Based on the simplified orbital dynamics equation, the precise value of the extended state quantity at the moment to be estimated is obtained through the initial extended state quantity, the radar observation vector, and the unscented batch filtering method, as follows:

[0063] Construct a system nonlinear state function according to the orbital dynamics equation; construct a measurement function according to the conversion relationship between the radar measured data or satellite ephemeris and the observation vector;

[0064] Estimate the initial extended state quantity by combining the dimension of the extended state vector, the initial value of the extended state quantity and covariance at the moment to be estimated, and obtain the Sigma point set through unscented transformation;

[0065] It should be noted that the Sigma point set is a key concept in the Unscented Kalman Filter (UKF) algorithm. The Sigma point set approximates the probability distribution through nonlinear transformation, making the standard Kalman filter framework applicable to nonlinear systems.

[0066] Obtain the weight factor of the state quantity and the weight factor of the covariance according to the dimension of the extended state vector;

[0067] Based on the system nonlinear state function, each Sigma point is propagated to all measurement moments to obtain the value of the extended state quantity at the corresponding moment after propagation, and the theoretical value of the observation vector at the corresponding moment is calculated according to the measurement function;

[0068] At each visible measurement moment, the estimated value of the observation data at this moment is obtained by weighting the theoretical value of the observation vector calculated from the propagation values of all Sigma points combined with the weight factor of the state quantity. The estimated value vectors of the observation data at all measurement moments constitute the estimated value vector of the measurement value;

[0069] It should be noted that the visible measurement moment refers to the moment when effective observation data can be obtained during the filtering process.

[0070] According to the weight factor of the covariance, the set of theoretical values of the observation vector calculated from the propagation values of each Sigma point at all visible measurement moments, the estimated value vector of the measurement value, and the measurement noise covariance matrix, obtain the measurement value covariance matrix. The specific formula is as follows:

[0071]

[0072] According to the weight factor of the covariance, the extended state quantity Sigma point set at the moment to be estimated, the estimation of the extended state quantity at the moment to be estimated, and each SigmaThe theoretical value set of the observation vector calculated from the propagation values of the points at all visible measurement times and the estimated value vector of the measured values are used to obtain the cross-covariance matrix between the state quantity and the measured values. The specific formula is as follows:

[0073]

[0074] Among them, represents the transpose matrix, is the covariance matrix of the measured values, is the cross-covariance matrix between the state quantity and the observed values, is the weight factor of the covariance, is from the i th Sigma set of theoretical values of the observation vector calculated from the propagation values of all visible observation times of the point, is the vector composed of the estimated values of the observation data at all visible measurement times, is the matrix composed of the measurement noise covariance at all times, is Sigma the value of the point after propagation at the time to be estimated, is the value of the extended state quantity after propagation at the time to be estimated. At the selected epoch, the Sigma point and the propagation of the extended state quantity are equal to the previously estimated values. Therefore, , , is the Sigma point set at the time to be estimated, is the estimation of the extended state quantity at the time to be estimated;

[0075] The filtering gain is obtained according to the covariance matrix of the measured values and the cross-covariance matrix between the state quantity and the measured values;

[0076] The observation residual is obtained according to the radar observation vector and the estimated vector of the measured values;

[0077] The extended state quantity of the satellite at the time to be estimated is iterated according to the filtering gain and the observation residual until the iteration end condition is satisfied, and the precise value of the extended state quantity at the time to be estimated is output.

[0078] Preferably, the iteration of the extended state quantity of the satellite at the time to be estimated according to the filtering gain and the observation residual is carried out according to the following iteration formula:

[0079]

[0080] Among them, is the improved value of the extended state quantity, is the estimated value of the extended state quantity in the previous iteration, is the filtering gain, is the observation residual;

[0081] The iteration end conditions are as follows:

[0082]

[0083]

[0084] in, is the root mean square of the new observation residual, is the root mean square of the residual error of the previous observation, is a constant, is the matrix composed of the measurement noise covariance, is the total number of visible measurement moments, is the observed residual, is the transpose of the observed residuals.

[0085] S6, bring the precise value of the extended state quantity at the time to be estimated into the orbital dynamics equation to predict the orbit at other times.

[0086] The present invention utilizes all observation data collected within an observation period to perform one-time processing and iterative solution to determine the orbit, thereby reducing the dependence on initial state covariance and process noise estimation and improving the stability of the calculation process.

[0087] Based on the above method, this embodiment also discloses a traceless batch orbit determination system for continuous low-thrust maneuvering satellites, see Figure 5 ,include:

[0088] Radar observation vector acquisition module: used to obtain radar observation vector according to radar measured data or satellite ephemeris data combined with radar station coordinates;

[0089] Orbit initial state acquisition module: used to obtain the initial orbit state at the time to be estimated based on the radar observation vector and the Lambert algorithm;

[0090] Extended state quantity generation module: used to generate initial extended state quantity according to the initial orbit state at the time to be estimated;

[0091] Orbital dynamics equation construction module: used to construct simplified orbital dynamics equations in the ECI coordinate system, including the non-spherical perturbation of the earth, the atmospheric drag perturbation and the constant small thrust acceleration in the RTN coordinate system;

[0092] Extended state quantity precision value acquisition module: used to obtain the extended state quantity precision value at the time to be estimated through the initial extended state quantity, radar observation vector and untraceable batch filtering method based on the simplified orbital dynamics equation;

[0093] Orbit prediction module: It is used to substitute the precise value of the extended state quantity at the moment to be estimated into the simplified orbit dynamics equation for orbit prediction at other moments.

[0094] The present invention utilizes all the observation data collected within an observation period, conducts one-time processing and iterative solution to obtain the initial orbit state, reduces the dependence on the estimation of the initial state covariance and process noise, and improves the stability of the calculation process. By integrating the unscented transform and the batch processing method, it is not necessary to calculate the partial derivatives of the system and measurement equations with respect to the state quantity, reducing the complexity of formula derivation and improving the adaptability to nonlinear systems. The present invention can accurately estimate the equivalent continuous constant small thrust acceleration in the local coordinate system under the condition that the thrust acceleration value is unknown and the observation data is sparse, and effectively conduct high-precision orbit determination and prediction for low-earth orbit small thrust continuous maneuvering satellites. Through the combination of the unscented transform and batch processing filtering, the present invention significantly improves the orbit determination accuracy, and reduces the influence of single-point observation errors by batch processing the observation data.

[0095] Embodiment 2:

[0096] See Figure 1 , this embodiment discloses an unscented batch processing orbit determination method for a continuous small thrust maneuvering satellite, including: obtaining the radar observation vector within the visible arc segment of the selected time period, where the radar observation vector includes slant range, azimuth angle, and elevation angle; using the observation vector of the visible arc segment where the moment to be estimated is located for initial orbit determination; based on the orbit dynamics model of the continuous small thrust maneuvering satellite including the extended state quantity (position, velocity, and thrust acceleration vector), conducting orbit improvement through the radar observation vector and the unscented batch processing filtering method; and finally using the obtained precise value of the extended state quantity to conduct orbit prediction for the continuous small thrust maneuvering satellite. Figure 1 Shown is the overall flowchart for orbit determination of a continuous small thrust maneuvering satellite based on unscented batch processing filtering. Figure 2 Shown is the algorithm flow of the unscented batch processing filtering. The specific technical solution of the present invention is as follows:

[0097] S1. Obtain the radar observation vector , including slant range ( ), azimuth angle ( ), and elevation angle ( ).

[0098] If there are radar measured data and the radar station site coordinates (longitude, latitude, and altitude), extract the radar observation vector of the target from the measured data.

[0099] If there is no measured data, given the coordinates of the radar station site, the positions and velocities of the satellite in the Earth-Centered, Earth-Fixed (ECEF) coordinate system are calculated using the Two-Line Element (TLE) data of the satellite for extrapolation or the satellite positions and velocities provided by the precise orbit data file. Then, the position of the satellite relative to the observation station is calculated and transformed into the local Northeast-Up coordinate system of the observation station, and the theoretical value of the observation vector is further calculated. Combining with the current radar's observation ability for low-Earth orbit targets, random observation errors are added to the three parameters of slant range, azimuth angle, and elevation angle respectively, and detection visibility constraints are imposed on the slant range and elevation angle to simulate and generate pseudo-radar observation vectors.

[0100] S2, determination of the initial orbital state of the satellite: Determine the position and velocity vectors of the satellite in the Earth-Centered Inertial (ECI) coordinate system at the moment to be estimated

[0101] (included in the radar observation visibility time). The initial orbital state is obtained as follows: Assume that there are in total n groups of radar observation vectors in the observation arc segment containing the moment to be estimated , where

[0102] is a positive integer not less than 3. If = 1 or , then according to the first observation point and the last observation point of this visible observation arc segment, the position vectors of the satellite at the moments and are obtained respectively. Further, without considering the continuous small thrust acceleration, the initial orbital state of the satellite at the moment is obtained using the Lambert algorithm;

[0103] If ≠ 1 and ≠ , then using or and the observation data at the moments , without considering the continuous small thrust acceleration, the initial orbital state of the satellite at the moment is obtained using the Lambert algorithm.

[0104] Among them, is the radar observation vector of the first observation point within the observation arc segment, is the corresponding slant range, is the corresponding azimuth angle, is the corresponding pitch angle; the -th radar observation vector of the observation points within the observation arc segment, is the corresponding slant range, is the corresponding azimuth angle, is the corresponding pitch angle; the -th radar observation vector of the observation point corresponding to the time within the observation arc segment, is the corresponding slant range, is the corresponding azimuth angle, is the corresponding pitch angle; is the orbital state at the -th moment, is the initial position estimate at the -th moment, is the initial velocity estimate at the -th moment.

[0105] S3. Establish the orbital dynamics equation of the continuously small-thrust maneuvering satellite including the extended state variables. The extended state variables include the position and velocity vectors in the Earth-centered inertial (ECI) coordinate system, and the small-thrust acceleration vector in the orbital RTN coordinate system. The extended state variables , and are the position and velocity vectors of the satellite in the ECI coordinate system respectively, is the thrust acceleration vector of the satellite in the orbital RTN coordinate system, including , and .

[0106] For a continuously maneuvering low-Earth orbit satellite, the Earth's non-spherical perturbation and atmospheric drag perturbation are the main perturbation factors. Therefore, we adopt a simplified motion model that only considers these two perturbation factors and perform state augmentation, including the parameter to be estimated - the small-thrust acceleration vector in the orbital coordinate system:

[0107]

[0108] In the formula, is the position vector of the satellite in the ECI coordinate system, are the velocity vectors of the satellite in the ECI coordinate system respectively, is the first-order derivative of the position vector of the satellite with respect to time, is the first-order derivative of the velocity vector of the satellite with respect to time, is the gravitational acceleration at the Earth's center, is the gravitational constant of the Earth; is the acceleration caused by the Earth's non-spherical perturbation; is the acceleration perturbation due to atmospheric drag, is the radial component of the thrust acceleration in the RTN coordinate system, is the along-track component of the thrust acceleration in the RTN coordinate system, is the normal component of the thrust acceleration in the RTN coordinate system. The thrust acceleration vector of the satellite in the orbital RTN coordinate system includes , and . is the first derivative of with respect to time, is the first derivative of with respect to time, is the coordinate transformation matrix from the RTN coordinate system to the ECI coordinate system, is expressed as:

[0109]

[0110] Denote and as the system non-linear state function and the measurement function respectively. The system non-linear state function is the established orbital dynamics model of the continuous low-thrust maneuvering satellite, and the measurement function is to calculate the theoretical observations based on the extended state variables , and . Since this method uses a batch processing technique to provide state estimates at the selected moments using a set of measurement data, and the batch processing process converges through iteration, there is no need for process noise compensation. Therefore, the non-linear system can be expressed as:

[0111]

[0112] where is the state vector at time with covariance , and is the measurement vector. is the additional measurement noise vector, which is a zero-mean Gaussian distribution with covariance , is the system non-linear state function, is the measurement function.

[0113] S4, according to the moment to be estimated The initial position and velocity vectors of the satellite are set, and the initial acceleration vector is set to 0 to obtain the estimated extended state variables of the satellite at the moment , and a suitable covariance estimate of the initial extended state variables is selected . According to the initial filtering values, the initial extended state variables are calculated through unscented transformation in the following manner Sigma points, Sigma The points are sample points, and the initial filtering values include the estimated extended state variables and the covariance estimate;

[0114]

[0115]

[0116]

[0117] where is the estimated extended state variables of the satellite at the moment, is the moment to be estimated According to generated Sigma point set, L is the dimension of the extended state vector, λ The definition of is: , α is the scaling factor, , controlling α The value of can control Sigma the distribution degree of the point set near , usually set to a very small positive number, = 3 – L .

[0118] The non-recursive unscented batch filter uses all the measurement data sets to estimate the extended state variables at the selected epoch moments. The propagation of the Sigma points at the selected epoch is equal to the previously estimated value, and the propagated state and covariance at the selected epoch can be directly set as follows:

[0119] ,

[0120]

[0121]

[0122] In the formula, the "-" symbol on each parameter represents the value of the corresponding parameter after propagation, is Sigma the value of the point after propagation at the moment to be estimated, For the value of the propagated state quantity at the moment to be estimated, For the value of the propagated covariance estimate of the state quantity to be estimated at the moment to be estimated.

[0123] Set as the weight factor of the state quantity, as the weight factor of the covariance, and it is calculated according to the following formula:

[0124]

[0125]

[0126] Wherein, β is the scaling factor, which reflects the high-order characteristics of the state history information. For the Gaussian distribution β = 2 is optimal.

[0127] S5, use the selected state quantity to be estimated at the moment point set Sigma of the ( i = 0, 1, 2, …, 2 L ) and the established orbital dynamics model of the continuous small-thrust maneuvering satellite, and propagate each Sigma point to all measurement moments ( , j = 1, 2, …, N , j ≠ k , N is the total number of visible measurement moments), and obtain the value of the propagated state quantity at the corresponding moment . = When, record . Calculate the theoretical value of the observation vector at the corresponding moment according to the propagated value of the state quantity and calculate the estimated value of the observation data by weighting , record as the set of theoretical values of the observation vector calculated from the propagated values of all visible observation moments of the i th Sigma point, is the vector composed of the estimated values of the observation data at all visible measurement moments, that is, the estimated value of the measurement value. is the measurement noise covariance matrix, and further calculate the measurement value covariance matrix and the cross-covariance matrix of the state quantity and the observed value :

[0128]

[0129]

[0130]

[0131]

[0132] wherein, is the covariance matrix of measurement values, is the cross-covariance matrix between the state quantity and the observation value.

[0133] S6. Calculate the filtering gain according to the covariance matrix of measurement values and the cross-covariance matrix between the state quantity and the observation value , specifically as follows:

[0134]

[0135] S7. Calculate the observation residual according to all the radar observation data at visible times obtained in step 1 and the estimated value of the measurement value obtained in step 5 , specifically as follows:

[0136]

[0137] wherein, is the calculated observation residual, is the radar observation vector, is the estimated value of the measurement value.

[0138] S8. Iteratively improve the extended state quantity of the satellite to be estimated at the time according to the filtering gain and the observation residual

[0139] until the iteration end condition is satisfied, and output the precise value of the extended state quantity at the epoch to be estimated.

[0140] ;

[0141] wherein, is the improvement value of the extended state quantity, is the value of the extended state quantity estimate after propagation, is the filtering gain, is the observation residual

[0142] When the absolute value of the difference between the root mean square (RMS) value of the measurement residual and the previous update result is less than a certain determined convergence criterion ε , the iteration ends:

[0143]

[0144]

[0145] Among them, is the root mean square of the new observation residuals, is the root mean square of the previous observation residuals, is a constant, is the matrix composed of the measurement noise covariance, is the total number of visible measurement times, is the observation residual, is the transpose of the observation residual.

[0146] S9. Substitute the precise value of the extended state quantity at the epoch to be estimated into the established orbit dynamics model to perform orbit prediction on the continuously small-thrust maneuvering satellite.

[0147] Compared with the prior art, the present invention has the following technical effects:

[0148] This method combines the unscented transform and the batch processing method. On the one hand, it reduces the complexity of formula derivation (without the need to calculate the partial derivatives of the system and measurement equations with respect to the state quantity) and improves the adaptability to the nonlinear system; on the other hand, it uses all the observation data collected during an observation period for one-time processing and iterative solution to obtain the orbit parameters, reducing the dependence on the initial state covariance and process noise estimation and improving the stability of the calculation process; at the same time, the calculation process does not involve the update of the state covariance, avoiding the influence brought by the possible diffusion of the state covariance. In short, this method can accurately estimate the equivalent continuous constant small-thrust acceleration in the local coordinate system under the conditions of unknown thrust acceleration value and sparse observation data, and effectively perform high-precision orbit determination and prediction on the low-orbit small-thrust continuously maneuvering satellite.

[0149] Based on the orbit dynamics equation of the continuously small-thrust maneuvering satellite including the extended state quantity (position, velocity, and thrust acceleration vector), the present invention improves the adaptability to the system nonlinearity, reduces the dependence on the initial state covariance and process noise estimation, and avoids the influence brought by the possible diffusion of the state covariance through the radar observation vector and the unscented batch filtering method, and can achieve high-precision orbit determination and prediction for the low-orbit continuously small-thrust maneuvering satellite.

[0150] Embodiment 3:

[0151] In one embodiment, as Figure 1 shown, an unscented batch orbit determination method for a continuously small-thrust maneuvering satellite is provided, including the following steps:

[0152] S1. Obtain the radar observation vector , including the slant range (ρ), azimuth angle (Az), and elevation angle (El).

[0153] Select the precise orbit data of the "Starlink" satellite 58826 during the in-orbit climbing phase for one day starting from the epoch time 2024-02-07 21:16:42 UTC as the "true value" of the orbit. The site coordinates (longitude, latitude, and altitude) of the measurement station are [75.99 deg, 39.47 deg, 0 km]; combined with the current radar's observation ability for low-Earth orbit targets, random observation errors with standard deviations of 100 m, 0.028°, and 0.028° (100 arcseconds) are added to the three parameters of slant range, azimuth angle, and elevation angle respectively. It is assumed that the radar can only detect when the slant range does not exceed 2000 km and the elevation angle is greater than 5°. Pseudo-radar observation vectors are simulated at a data rate of 1 / 30 (i.e., one set of observation data is generated every 30 seconds). A total of 7 arcs (see Table 1) and 70 sets of observation data are obtained. The time interval between the first 6 adjacent visible arcs is approximately 1.5 hours, and the time interval between the 7th and 6th visible arcs is approximately 15 hours.

[0154] Table 1. Simulation results of the observation data of satellite 58826 for one day starting from 2024-02-07 21:16:42 UTC:

[0155]

[0156] S2. Use the radar observation vector to determine the initial orbit state of the target at the selected moment to be estimated.

[0157] Select the first visible moment as the moment to be estimated, that is, = 2024-02-07 22:06:42. According to the observation vectors at the beginning and end points of the first visible arc, use the Lambert algorithm to obtain the initial orbit state of the target at The moment is:

[0158]

[0159] = [-5767427.02857054, 739786.320264708, 3439130.13587933] m,

[0160] = [2111.39988515649, -5661.00221865180, 4746.77719102793] m / s.

[0161] S3. Establish a simplified perturbed orbit dynamics equation, system state equation, and measurement equation that include continuous small thrust vectors.

[0162] Model the orbit dynamics equation considering only the Earth's non-spherical perturbation, atmospheric drag perturbation, and constant continuous small thrust as:

[0163]

[0164] In the formula, is the position vector of the satellite in the ECI coordinate system, are the velocity vectors of the satellite in the ECI coordinate system respectively, is the gravitational acceleration at the Earth's center; is the acceleration caused by the Earth's non-spherical perturbation; is the atmospheric drag perturbation acceleration, 、 、 are the radial, along-track, and normal components of the thrust acceleration in the RTN coordinate system respectively, is the coordinate transformation matrix from the RTN coordinate system to the ECI coordinate system, expressed as:

[0165]

[0166] Denote and as the system nonlinear state function (i.e., the orbital dynamics equation of the continuously small-thrust maneuvering satellite established by us) and the measurement function (calculating the theoretical observation vectors ρ, Az, and El according to the state quantity x) respectively. Then the nonlinear system can be expressed as:

[0167]

[0168] In the formula, is the state vector at time with covariance and is the measurement vector. is the additional measurement noise vector, which is a zero-mean Gaussian distribution with covariance . Taking the standard deviation of the measurement noise as 100 m for the slant range, 0.028° for the azimuth angle, and 0.028° (100 arcseconds) for the elevation angle, then

[0169]

[0170] S4, determine the initial value of the extended state quantity estimation and the covariance estimation of the satellite at the moment to be estimated. Calculate the Sigma points of the initial extended state quantity through the unscented transformation, and assign the weight factors of the state quantities and covariances of each point.

[0171] According to the fact that the orbital plane of the Starlink satellite remains almost unchanged during the climbing stage, it can be known that there is no thrust in the normal direction of the orbital plane. Therefore, in this embodiment, only the small-thrust accelerations in the radial and along-track directions 、 It is included as an extended item in the state quantity for estimation. Since the thrust acceleration of the non-maneuvering target is unknown, the initial value of the acceleration is taken as , combined with the estimated time calculated in step 2 the initial position and velocity vector of the satellite, to obtain the estimation of the extended state quantity of the satellite at time -4 m / s 2 Select the standard deviations of the errors in the directions of the position, velocity, and acceleration vectors to be 1000 m, 10 m / s, and 10 .

[0172] According to the initial filtering values ( and ), calculate the Sigma points of the initial extended state quantity through unscented transformation in the following manner:

[0173]

[0174]

[0175]

[0176] where is the dimension of the extended state vector. In this embodiment = 8, the scaling factor α is taken as 0.0001, = 3 – = -5, .

[0177] At the time to be estimated, the propagation of the Sigma point is equal to the previously estimated value, and the propagated state and covariance at the time can be directly set as follows:

[0178] ,

[0179]

[0180]

[0181] In the formula, the "-" symbol on each parameter represents the value of the corresponding parameter after propagation.

[0182] Set as the weight factor of the state quantity, as the weight factor of the covariance, and calculate according to the following formula:

[0183]

[0184]

[0185] Among them, β is the scaling factor, and in this embodiment, β = 2 is taken.

[0186] S5. According to the Sigma point set of the initial extended state quantity at the moment to be estimated and the orbit dynamics model, obtain the Sigma propagation values of the extended state vectors and the corresponding theoretical values of the observation vectors at all radar visible measurement moments corresponding to each point, and calculate the estimated values of the observation vectors at all visible measurement moments through weighted calculation.

[0187] Using the extended state quantity Sigma point set ( i = 0, 1, 2, …, 2 L ) and the orbit dynamics equation, propagate each Sigma point to all measurement moments ( t j , j = 1, 2, …, N , j ≠ k , N where is the total number of visible measurement moments), to obtain the values of the extended state quantity at the corresponding moments after propagation . Calculate the theoretical values of the observation vectors at the corresponding moments according to the propagation values of the extended state quantity and calculate the estimated values of the observation data through weighted calculation. Denote Sigma as the set of theoretical values of the observation vectors calculated from the propagation values at all visible observation moments of the first point, and

[0188]

[0189]

[0190] as the vector composed of the estimated values of the observation data at all visible measurement moments:

[0191] Calculate the covariance matrix of the measurement values and the cross-covariance matrix of the state quantity and the observation values, and then calculate the filtering gain. , , , obtained in step 5, and the measurement noise covariance matrix , and the measurement noise covariance matrix , calculate the covariance matrix of the measurement values and the cross-covariance matrix of the state quantity and the observation values

[0192]

[0193]

[0194] S6, Calculate the filtering gain :

[0195] 。

[0196] S7, Calculate the observation residual, and use the filtering gain to iteratively improve the satellite extended state vector at the

[0197] time until the iteration end condition is satisfied. and the measurement value estimate obtained in step 5 Calculate the observation residual as follows:

[0198] 。

[0199] S8, Update the satellite extended state vector according to the following formula

[0200] 。

[0201] When the absolute value of the difference between the root mean square (RMS) value of the observation residual and the previous update result is less than a certain determined convergence criterion ε, the iteration ends:

[0202]

[0203]

[0204] In this embodiment, ε = 10 -4 。

[0205] Obtain the precise value of the extended state quantity at the

[0206] = [-5768806.58821858, 738521.548084285, 3438720.37494993] m;

[0207] = [2116.53368902988, -5656.95526541309, 4749.96272962505] m / s.

[0208] Check the true ephemeris value at this time as:

[0209] = [-5768677.37003690, 738604.560328600, 3438684.61668840] m;

[0210] = [2116.10198940000, -5657.01619310000, 4750.24608580000] m / s.

[0211] The error between the precise value of the extended state quantity and the true ephemeris value is:

[0212] = [-129.218181683682, -83.0122443146538, 35.7582615278661] m;

[0213] = [0.431699629875311, 0.0609276869117821, -0.283356174951223] m / s.

[0214] The estimated radial and along-track small thrust accelerations are:

[0215] = -0.000147363463861147 m / s2

[0216] = 7.08971836251527e-05 m / s2

[0217] Substitute the precise value of the extended state quantity at the moment into the established orbital dynamics equation to predict the orbit of the continuously small-thrust maneuvering satellite. The error between the prediction result for two days and the target true ephemeris value is as Figure 3 shown. Figure 3 The satellite position error in the RTN coordinate system. It can be seen that the maximum position error within two days does not exceed 2000 m.

[0218] An electronic device, comprising: a processor; a memory for storing computer program instructions; and for implementing the unscented batch orbit determination method for continuously small-thrust maneuvering satellites when executing the computer program.

[0219] A storage medium storing computer program instructions, when the computer program instructions are loaded and run by a processor, the processor executes the unscented batch orbit determination method for continuously small-thrust maneuvering satellites.

[0220] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk memory, CD-ROM, optical memory, etc.) that contain computer-usable program code.

[0221] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to embodiments of the present invention. It should be understood that each flow and / or block in the flowchart and / or block diagram, and the combination of flows and / or blocks in the flowchart and / or block diagram, can be realized by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, such that the instructions executed by the processor of the computer or other programmable data processing devices generate means for realizing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0222] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing devices to work in a specific manner, such that the instructions stored in the computer-readable memory generate a manufactured article including instruction means that realizes the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0223] These computer program instructions can also be loaded onto a computer or other programmable data processing devices, such that a series of operation steps are executed on the computer or other programmable devices to generate a computer-implemented process, and thus the instructions executed on the computer or other programmable devices provide steps for realizing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0224] The above content is only to illustrate the technical idea of the present invention and cannot be used to limit the protection scope of the present invention. Any modification made on the basis of the technical solution according to the technical idea proposed by the present invention falls within the protection scope of the present invention.

Claims

1. An unscented batch orbit determination method for continuously small-thrust maneuvering satellites, characterized in that, It includes the following steps: When there are radar measured data and radar site coordinates, obtain the radar observation vector based on the radar measured data combined with the radar site coordinates; when there are no measured data, obtain the pseudo-radar observation vector based on the satellite ephemeris data; Obtain the initial orbital state at the moment to be estimated according to the radar observation vector and using the Lambert algorithm; Generate the initial extended state variables according to the initial orbital state at the moment to be estimated; Construct a simplified orbital dynamics equation in the ECI coordinate system, including the non-spherical perturbation of the earth, the atmospheric drag perturbation, and the constant small thrust acceleration in the RTN coordinate system, as follows: The extended state variables include the position and velocity vectors of the satellite in the ECI coordinate system, and the constant small thrust acceleration vector of the satellite in the orbital RTN coordinate system. Among them, the thrust acceleration vector of the satellite in the orbital RTN coordinate system includes the radial component, the along-track component, and the normal component of the thrust acceleration in the RTN coordinate system; The simplified orbital dynamics equation is as follows: In the formula, is the position vector of the satellite in the ECI coordinate system, is the velocity vector of the satellite in the ECI coordinate system, is the first-order derivative of the position vector of the satellite with respect to time, is the first-order derivative of the velocity vector of the satellite with respect to time, is the gravitational acceleration at the Earth's center, is the gravitational constant of the Earth; is the acceleration caused by the Earth's non-spherical perturbation; is the perturbation acceleration of atmospheric drag, is the radial component of the thrust acceleration in the RTN coordinate system, is the along-track component of the thrust acceleration in the RTN coordinate system, is the normal component of the thrust acceleration in the RTN coordinate system, is the first-order derivative with respect to time, is the first-order derivative with respect to time, is the first-order derivative with respect to time, is the coordinate transformation matrix for converting from the RTN coordinate system to the ECI coordinate system; Based on the simplified orbital dynamics equation, obtain the precise value of the extended state variables at the moment to be estimated through the initial extended state variables, the radar observation vector, and the unscented batch filtering method; Substitute the precise value of the extended state variables at the moment to be estimated into the simplified orbital dynamics equation for orbit prediction at other moments.

2. The unscented batch orbit determination method for a continuously small-thrust maneuvering satellite according to claim 1, wherein When there are radar measured data and radar site coordinates, obtain the radar observation vector based on the radar measured data combined with the radar site coordinates, specifically as follows: If there are radar measured data and radar site coordinates, extract the radar observation vector from the radar measured data according to the radar measured data and the radar site coordinates. The radar observation vector includes the slant range, azimuth angle, and elevation angle of the satellite; When there are no measured data, obtain the pseudo-radar observation vector based on the satellite ephemeris data, specifically as follows: Given the radar site coordinates, extrapolate the position and velocity of the satellite using the two-line orbital element data of the satellite based on the given radar site coordinates, or obtain the position and velocity of the satellite through the precise orbit data file; Calculate the position and velocity of the satellite and the observation station in the geodetic coordinate system according to the position and velocity of the satellite, and calculate the position of the satellite relative to the observation station; Convert the position and velocity of the satellite in the geodetic coordinate system and the position of the satellite relative to the observation station to the local northeast celestial coordinate system of the observation station, and calculate the theoretical value of the observation vector; Add random observation errors to the three parameters of the slant range, azimuth angle, and elevation angle in the theoretical value of the observation vector respectively, and impose the detection visibility limit conditions on the slant range and elevation angle to simulate and generate the pseudo-radar observation vector.

3. The unscented batch orbit determination method for a continuously small-thrust maneuvering satellite according to claim 1, characterized in that Obtain the initial orbital state at the moment to be estimated according to the radar observation vector and using the Lambert algorithm, specifically as follows: There are at least three radar observation vectors in the observation arc where the moment to be estimated is located; When the radar observation vector at the moment to be estimated is the first or the last radar observation vector in the observation arc, obtain the position vectors of the satellite at the first and the last radar observation vectors in the observation arc respectively according to the first and the last radar observation vectors in the observation arc, and use the Lambert algorithm to obtain the initial orbital state of the satellite at the moment to be estimated; When the radar observation vector at the to-be-estimated moment is between the first radar observation vector and the last radar observation vector in the observation arc segment, based on the first radar observation vector or the last radar observation vector within the observation arc segment, combined with the radar observation vector at the to-be-estimated moment and using the Lambert algorithm, the initial orbital state of the satellite at the to-be-estimated moment is obtained; The initial orbital state at the to-be-estimated moment includes the initial position and velocity vector of the satellite at the to-be-estimated moment.

4. The unscented batch orbit determination method for a continuously small-thrust maneuvering satellite according to claim 1, wherein The generation of the initial extended state quantity according to the initial orbital state at the to-be-estimated moment is specifically as follows: Adding a constant small thrust acceleration vector in the RTN coordinate system to the orbital state at the to-be-estimated moment to obtain the extended state quantity; The initial value of the small thrust acceleration vector is set to zero.

5. The unscented batch orbit determination method for a continuously small-thrust maneuvering satellite according to claim 1, characterized in that, Based on the simplified orbital dynamics equation, through the initial extended state quantity, the radar observation vector, and the unscented batch filtering method, the precise value of the extended state quantity at the to-be-estimated moment is obtained, specifically as follows: Constructing a system nonlinear state function according to the simplified orbital dynamics equation; constructing a measurement function according to the conversion relationship between the radar measured data or satellite ephemeris and the observation vector; According to the dimension of the augmented state vector, the augmented state quantity at the moment to be estimated, and the initial value of the covariance, an initial augmented state quantity is obtained by combining the unscented transform and estimating the Sigma point set; Obtaining the weight factor of the state quantity and the weight factor of the covariance according to the dimension of the extended state vector; Based on the system's non-linear state function, each Sigma point is propagated to all measurement times to obtain the values of the extended state variables at the corresponding times after propagation, and the theoretical values of the observation vectors at the corresponding times are calculated according to the measurement function. At each visible measurement moment, the theoretical value of the observation vector calculated from all Sigma point propagation values is weighted by the weight factor of the state quantity to obtain the estimated value of the observation data at this moment. The estimated values of the observation data at all measurement moments form the estimated value vector of the measurement value; The theoretical value set of the observation vector, the estimated value vector of the measured value, and the measurement noise covariance matrix calculated according to the weight factor of the covariance and the propagation values of each point at all visible measurement times are used to obtain the measurement value covariance matrix; Sigma ​ Weight factor according to covariance, extended state quantity at the moment to be estimated Sigma Point set, estimation of extended state quantity at the moment to be estimated, each Sigma The theoretical value set of the observation vector calculated from the propagation values of each point at all visible measurement moments and the estimated value vector of the measured value are used to obtain the cross-covariance matrix between the state quantity and the measured value; Obtaining the filtering gain according to the measurement value covariance matrix and the cross-covariance matrix between the state quantity and the measurement value; Obtaining the observation residual according to the radar observation vector and the estimated value vector of the measurement value; Iterating the extended state quantity of the satellite at the to-be-estimated moment according to the filtering gain and the observation residual until the iteration end condition is satisfied, and outputting the precise value of the extended state quantity at the to-be-estimated moment.

6. The unscented batch orbit determination method for a continuously small-thrust maneuvering satellite according to claim 5, characterized in that The iterative improvement of the extended state quantity of the satellite at the to-be-estimated moment according to the filtering gain and the observation residual, and the iterative formula is as follows: Among them, is the improved value of the extended state variable, is the estimated value of the extended state variable, is the filtering gain, is the observation residual; The iteration end condition is as follows: wherein, is the root mean square of the new observation residual, is the root mean square of the observation residual of the previous iteration, is a constant, is the measurement noise covariance matrix, is the total number of visible measurement times, is the observation residual, is the transpose of the observation residual.

7. A UKF batch orbit determination system for a continuously small-thrust maneuvering satellite, characterized in that, Including: Radar observation vector acquisition module: used to obtain the radar observation vector according to the radar measured data combined with the radar site coordinates when there are radar measured data and radar site coordinates, and obtain the pseudo radar observation vector according to the satellite ephemeris data when there is no measured data; Orbital initial state acquisition module: used to obtain the initial orbital state at the to-be-estimated moment according to the radar observation vector and using the Lambert algorithm; Extended state quantity generation module: used to generate the initial extended state quantity according to the initial orbital state at the to-be-estimated moment; Orbital dynamics equation construction module: used to construct a simplified orbital dynamics equation including the non-spherical perturbation of the earth, the atmospheric drag perturbation, and the constant small thrust acceleration in the RTN coordinate system in the ECI coordinate system, specifically as follows: The extended state quantity includes the position and velocity vector of the satellite in the ECI coordinate system, and the constant small thrust acceleration vector of the satellite in the orbital RTN coordinate system, where the thrust acceleration vector of the satellite in the orbital RTN coordinate system includes the radial component, the along-track component, and the normal component of the thrust acceleration in the RTN coordinate system; The simplified orbital dynamics equation is specifically as follows: In the formula, is the position vector of the satellite in the ECI coordinate system, is the velocity vector of the satellite in the ECI coordinate system, is the first derivative of the position vector of the satellite with respect to time, is the first derivative of the velocity vector of the satellite with respect to time, is the gravitational acceleration at the Earth's center, is the gravitational constant of the Earth; is the acceleration caused by the Earth's non-spherical perturbation; is the perturbation acceleration of atmospheric drag, is the radial component of the thrust acceleration in the RTN coordinate system, is the component of the thrust acceleration in the along-track direction in the RTN coordinate system, is the normal component of the thrust acceleration in the RTN coordinate system, is the first derivative with respect to time, is the first derivative with respect to time, is the first derivative with respect to time, is the coordinate transformation matrix for converting from the RTN coordinate system to the ECI coordinate system; Extended state quantity precise value acquisition module: It is used to obtain the precise value of the extended state quantity at the moment to be estimated based on the simplified orbital dynamics equation, through the initial extended state quantity, radar observation vector, and unscented batch filtering method; Orbit prediction module: It is used to substitute the precise value of the extended state quantity at the moment to be estimated into the simplified orbital dynamics equation for orbit prediction at other moments.

8. An electronic device, comprising: A processor; a memory, and the electronic device is used to store computer program instructions; it is characterized in that when the computer program is executed, it realizes the unscented batch orbit determination method for a continuously small-thrust maneuvering satellite as described in any one of claims 1-6.

9. A storage medium storing computer program instructions, characterized in that, When the computer program instructions are loaded and run by the processor, the processor executes the unscented batch orbit determination method for a continuously small-thrust maneuvering satellite as described in any one of claims 1-6.

Images