Unknown target short-arc high-precision initial orbit determination method based on filter enhancement

Abstract

The application discloses a kind of unknown target short arc high-precision initial orbit determination method based on filter enhancement, first according to observation data type adaptation initial orbit determination method, after obtaining multiple groups of initial orbit candidate solution, determine initial orbit optimal solution according to the minimum principle of measurement residual, and convert into position, velocity vector.Combining target orbit type and physical characteristics, configure dynamic model and parameter;Based on the model, the state equation of the target is established, and the process noise compensation model error is introduced;Combined with the observation platform position information, the observation equation is constructed.To ensure the effective start of the filter, the orbit state and state error covariance, process noise and measurement noise parameters need to be calculated and initialized in advance.Finally, combined with the state equation and observation equation, the target orbit state is predicted and updated through EKF technology until all observation epochs are traversed, and high-precision target initial orbit information is obtained.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of space orbit dynamics and space target monitoring, and more particularly, relates to a short-arc high-precision initial orbit determination method for unknown targets based on filter enhancement. BACKGROUND

[0002] Initial orbit determination (IOD) refers to the determination of the orbit parameters of an unknown space target based on the obtained observation arc data under the condition of lacking prior information of the space target orbit. With the frequent development of space activities, the number of on-orbit space targets has increased dramatically, which has put forward higher requirements for space environment management and task safety guarantee. For space targets, especially non-cooperative spacecraft with maneuvering capability, accurate initial orbit determination is an important prerequisite for space target cataloging and collision warning. At the same time, high-precision initial orbit determination can also be used to support non-cooperative target anomaly detection and identification, maneuvering force inversion, disintegration tracing, and other tasks, and provide important technical support for space target situation awareness.

[0003] Currently, ground-based radar and space / ground-based optical observation are important means to obtain the orbit data of near-earth space targets. For non-cooperative space targets, due to the limitations of observation resource allocation, target motion characteristics and other factors, only short-time observation data can be obtained, and the arc time span is usually less than 1% of the observed target orbit period, which is considered as ultra-short-arc data. Using such data for initial orbit determination faces the challenges of limited data information and weak geometric configuration. The classical initial orbit determination algorithms (such as Laplace method, Gauss method, double-R iteration method, etc.) can achieve good results when dealing with longer arc data, but in the case of short-arc / ultra-short-arc, the success rate and precision of initial orbit determination cannot be guaranteed. In modern IOD technology, the Gooding method, as an internationally recognized initial orbit determination algorithm, has good convergence and stability; the distance search method optimizes the initial orbit determination capability of optical observation data by adding distance constraints based on target type (see Zhang Pin's paper "Low-orbit space debris initial orbit determination method based on distance search" published in 2017). These two methods can maintain a high orbit determination success rate in short-arc / ultra-short-arc observation scenarios, but the initial orbit determination precision cannot meet the needs of high-precision space applications, and further orbit refinement is needed.

[0004] The above IOD technology can provide the initial value of the orbit of an unknown target, providing necessary prior information for further improvement of the orbit. However, most existing orbit refinement methods are based on batch least squares (BLS) algorithm, which is prone to ill-conditioned normal equation when facing short-arc / ultra-short-arc data, and an effective orbit refinement path needs to be explored. SUMMARY

[0005] To overcome the above deficiencies of the prior art, the application provides a short-arc high-precision initial orbit determination method for unknown targets based on filter enhancement. The method adopts a two-step method of "initial orbit determination + filter enhancement": first, an initial orbit is determined by an algorithm with high orbit determination accuracy and success rate such as the Gooding method and the distance search method, to provide reliable initial values to accelerate filter convergence; and then, based on a high-precision orbit dynamics model, the initial orbit is refined through a filter algorithm for all available observation values, thereby effectively improving the accuracy and practicability of initial orbit determination based on short-arc / radar / optical observation data.

[0006] According to an aspect of the present application, a short-arc high-precision initial orbit determination method for unknown targets based on filter enhancement is provided, comprising:

[0007] According to the type of observation data, a method suitable for the type is used to determine the initial orbit, and an optimal initial orbit solution is determined from a plurality of initial orbit candidate solutions obtained;

[0008] The dynamics model is configured in combination with the type and physical characteristics of the target orbit;

[0009] The target state equation is established based on the dynamics model;

[0010] The observation equation is constructed in combination with the observation data and the position information of the observation platform;

[0011] In combination with the target state equation and the observation equation, the target orbit state is predicted and updated through extended Kalman filtering starting from the optimal initial orbit solution and the initial orbit state covariance, until all observation epochs are traversed, to obtain high-precision initial orbit information of the target.

[0012] The above technical solution should first be pre-processed according to the type of observation data and then determined by a method suitable for the type when performing subsequent steps, and the orbit is refined using the dynamics model and extended Kalman filtering technology to improve the accuracy of the initial orbit.

[0013] Further, various types of observation data can be used in initial orbit determination and orbit refinement, including but not limited to azimuth, elevation, two-way range, angle rate and range rate in the geocentric coordinate system, right ascension, declination and three-dimensional position of space targets in the inertial coordinate system, etc.

[0014] When determining the optimal initial orbit solution, various solution pool optimization methods can be used, and the principle of minimum measurement residual is used to determine the optimal initial orbit solution, i.e. . Wherein, is the semi-major axis, is the eccentricity, is the orbit inclination, is the longitude of the ascending node, is the argument of perigee, is the true anomaly, is the total number of observation data points, represents the group of observation values, is the theoretical value of the group of observation values. It should be pointed out that the solution pool optimization method refers to determining the optimal initial orbit parameters by using a certain method from a plurality of solutions obtained by using the initial orbit determination method.

[0015] The above technical solution sets the dynamic model parameters for the initial orbit refinement model, and the parameters should be configured according to the physical characteristics and the orbit type of the space target.

[0016] As a further technical solution, when the observation data type is space-based optical observation data, the process of determining the optimal initial orbit solution includes:

[0017] Observation data points are selected at the beginning and end of the arc segment, respectively, and are combined in pairs to form observation data point pairs;

[0018] The direction of the range vector is determined according to the right ascension and declination in the inertial coordinate system, and the initial and final ranges are determined by searching the large and small step distances according to the orbit type of the space target and the position information of the observation platform. Then, the initial and final range vectors are calculated by combining the direction of the range vector and the initial and final ranges, and a Lambert problem is constructed.

[0019] The Lambert problem is solved by using an initial orbit determination algorithm to obtain a plurality of groups of initial orbit candidate solutions;

[0020] The optimal initial orbit solution is determined based on the principle of minimum measurement residual, and is converted into an initial orbit state vector.

[0021] As a further technical solution, when the observation data type is ground-based optical observation data, the azimuth and elevation angles in the topocentric coordinate system obtained by ground-based observation are converted into the right ascension and declination in the inertial coordinate system, and the same method as for space-based optical observation data is used to determine the optimal initial orbit solution.

[0022] As a further technical solution, when the observation data type is ground-based radar observation data, the process of determining the optimal initial orbit solution includes:

[0023] Observation data points are selected at the beginning and end of the arc segment, respectively, and are combined in pairs to form observation data point pairs;

[0024] The satellite initial and final range vectors are calculated according to the station position, ranging and angular information, and a Lambert problem is constructed;

[0025] The Lambert problem is solved by using an initial orbit determination algorithm to obtain a plurality of groups of initial orbit candidate solutions;

[0026] The optimal solution for the initial orbit is determined based on the principle of minimizing the measurement residual, and then converted into the initial orbit state vector.

[0027] As a further technical solution, a target state equation is established based on the aforementioned dynamic model, including:

[0028] Based on the aforementioned dynamic model, the motion equations of the space target are constructed;

[0029] The motion equations of the space target are discretized, and Taylor expansion is performed at the initial orbit optimal solution to obtain the discretized state equations.

[0030] When constructing the target state equation, this technical solution generally selects a second-order differential equation. The specific expression is determined by the selected dynamic model, and white noise should be introduced into the equation. .in, The covariance is mainly determined by the unconsidered maximum perturbation force level, and is used to compensate for random errors such as nonlinearity errors, simplification of dynamic models and measurement models.

[0031] As a further technical solution, observation equations are constructed by combining the location information of the observation platform, including:

[0032] set up Observation vector at time It consists of angle measurement data, or angle measurement and distance measurement data at the same moment, satisfying

[0033] ,

[0034] in, Indicates measurement noise. and These are the observation vectors with respect to the instantaneous state vector. and initial orbital state Theoretical value;

[0035] residual vector Approximate value of the optimal solution for the initial trajectory The approximation is

[0036] ,

[0037] in, , , For model observation pairs The Jacobian matrix of the orbital state vector at time t=0 satisfies

[0038] , ,

[0039] in, For observation epochs The approximate value of the orbital state at time step is calculated from the target state equation. To update the state prediction results, the measurement residual vector is... Represented as time The function is calculated by the following formula. The New Message of Time:

[0040] .

[0041] When constructing the observation equations, the above technical solution should consider the prior orbital state at the observation epoch. Linearization is performed at the point, and the linearized residual vector satisfies .

[0042] As a further technical solution, combining the target state equation and the observation equation, starting from the optimal solution of the initial orbit and the covariance of the initial orbit state, the target orbit state is predicted and updated through extended Kalman filtering, including:

[0043] State prediction: Based on discretized state equations, utilizing... Time-based filtering improved state vector and state covariance matrix forecast Prior orbital state at time and prior state covariance matrix In this context, the superscript "-" indicates the prior orbital value calculated from the state equation, and the superscript "+" indicates the orbital state and covariance updated based on the observation data.

[0044] State update: Utilizing prior state Prior state covariance matrix and measurement noise covariance matrix Calculate Kalman gain ,use Obtain the updated target orbit state and the target state covariance matrix ;

[0045] Filtering from and Start-up, given process noise Observation vector and measurement noise covariance matrix , Repeat the state prediction and state update steps until all observation epochs are traversed to obtain a refined high-precision initial orbit.

[0046] The above technical solution requires the filter initialization results when refining the trajectory based on extended Kalman filtering. and The process begins by using state equations to predict the target orbital state, combined with observational data. The orbit is updated to obtain and .

[0047] Furthermore, when performing orbit refinement based on the extended Kalman filter technique, the termination condition for orbit refinement is to traverse all observation times and finally output the refined orbit state and corresponding covariance matrix for all observation times. It is generally believed that the orbit state at the last observation time has the highest accuracy and is used as the result of orbit refinement. The state can be orbital elements, or position and velocity vectors.

[0048] As a further technical solution, before predicting and updating the target orbit state through extended Kalman filtering, the following filter initialization is also included:

[0049] The covariance matrix of the white noise process is determined based on the accuracy of the dynamic model.

[0050] The optimal solution of the initial orbit As the initial orbital state of the filter, the initial orbital state covariance matrix is ​​set based on the initial orbital determination accuracy. ;

[0051] Initialize the measurement noise covariance matrix based on the accuracy of the observations. .

[0052] The above technical solution should use the initial trajectory optimal solution during filter initialization. As the initial orbital state of the filter, and based on experience or statistical data, the initial state covariance matrix is... Assign values; process noise covariance matrix Primarily determined by the simplification error of the force model; measurement noise covariance matrix The accuracy is determined by the measurement precision; if there is no specific measurement precision, it is given based on experience.

[0053] According to one aspect of the present invention, a device for determining the high-precision initial trajectory of an unknown target using a short arc based on enhanced filtering is provided, comprising a memory and a processor. The memory stores program instructions that are executed by the processor, and the processor invokes the program instructions to execute the method for determining the high-precision initial trajectory of an unknown target using a short arc based on enhanced filtering.

[0054] According to one aspect of the present invention, a non-transitory computer-readable storage medium is provided, the non-transitory computer-readable storage medium storing computer instructions that cause the computer to execute the described method for determining the high-precision initial trajectory of an unknown target using short arc based on filter enhancement.

[0055] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0056] Compared to traditional IOD algorithms and existing BLS refinement methods, the technical solution proposed in this invention can flexibly integrate various types of observation data, fully mine the information in the observation data, and incorporate more complex dynamic system models. Crucially, this technology can specifically overcome the ill-conditioned problem of the normal equations in the BLS method, making it suitable for high-precision initial orbit determination of short-arc / ultra-short-arc orbits for various orbital types of space targets. This invention employs a two-step method: using high-accuracy and high-success-rate initial orbit determination algorithms such as the Gooding method and distance search method for initial orbit determination, followed by orbit refinement using a filtering enhancement algorithm based on a high-precision orbital dynamic model. This significantly improves the accuracy of initial orbit parameter calculation. Attached Figure Description

[0057] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the accompanying drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0058] Figure 1 This is a flowchart illustrating the high-precision initial trajectory determination of an unknown target using short-arc / ultra-short-arc observation data and based on filtering enhancement technology, as described in an embodiment of the present invention.

[0059] Figure 2 This is a map showing the orbital altitude distribution of the experimental satellite calculated from precise ephemeris during the visible period of the station on the first day of the simulated observation period (i.e., November 28, 2024) in an embodiment of the present invention.

[0060] Figure 3 In this embodiment of the invention, a typical filtered three-dimensional position error convergence effect diagram calculated based on 180s short-arc ground-based radar observation data is shown, displaying the three-dimensional position errors of ultra-short arc / short arc at 15s, 30s, 60s, 120s, and 180s, respectively, and converted to the Along-track, Cross-track, and Radial RSW coordinate systems. Detailed Implementation

[0061] The terms “comprising” and “having”, and any variations thereof, in the specification, claims, and accompanying drawings of this invention are intended to cover a non-exclusive inclusion, such as a process, method, system, product, or apparatus that includes a series of steps or units, not necessarily limited to those explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0062] This invention provides a high-precision initial trajectory determination method for unknown targets using short-arc radar based on enhanced filtering. First, the initial trajectory is determined using an initial trajectory determination algorithm (preferably the Gooding method / range search method for radar / optical data). Then, the results are refined using extended Kalman filtering technology, which can effectively improve the accuracy and practicality of the initial trajectory determination results based on short-arc / ultra-short-arc radar / optical observation data.

[0063] The specific steps are as follows:

[0064] Step 1: Calculate the initial orbital elements using the appropriate IOD method based on the data type of the observation data (such as radar observation data, optical observation data, etc.).

[0065] 1.1 Space-based observation equipment typically records angle information as the platform's line-of-sight direction to a space target in the TOD (True of Date) inertial coordinate system, i.e., right ascension (RA) and declination (DEC). Ground-based observation equipment, on the other hand, records angle information as azimuth (Az) and elevation (El) in the station-centered coordinate system. During orbit calculation, the latter requires combining the station's precise position information in the TOD coordinate system to convert the angle information into right ascension and declination in the TOD coordinate system. Optical observation data only contains the aforementioned angle information, while ground-based radar observation data adds range information to the angle information, expressed as a two-way range (R).

[0066] 1.2 For given short arc / ultra-short arc observation data, select observation data points at the beginning and end of the arc segment respectively. and ( Pairs are combined to form several observation data point pairs. This is used to construct Lambert problems. For example, if 5 observation points are selected for the initial and final arc segments, a total of 25 Lambert problems can be constructed.

[0067] 1.3 Constructing the Lambert Problem:

[0068] 1.3.1 For ground-based radar observation data, the satellite's initial and final slant range vectors are calculated based on the station's location and ranging and angular information. and Construct the Lambert problem;

[0069] 1.3.2 For optical observation data, the direction of the slant range vector is first determined by right ascension and declination. and Next, the orbital type of the space target (such as LEO orbit, MEO orbit, etc.) is obtained, and based on the orbital type and the location information of the observation platform, the initial and final slant ranges of the target are determined sequentially through large and small step distance searches. and Finally, the initial and final slant range vectors of the space target are calculated using the station's location, slant range information, and angular measurement data. and This leads to the construction of the Lambert problem.

[0070] 1.4. Use Gooding's method to solve the Lambert problem for all combinations constructed in step 1.3, and obtain multiple sets of candidate solutions for the initial orbit. , ( That is, the initial orbital solution set. , , , , and These represent the semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and angle of apogee, respectively.

[0071] 1.5 Determining the optimal solution for the initial trajectory based on the principle of minimizing measurement residuals And convert it into the corresponding orbital state vector. This is used as the initial value of the orbital state for the orbital refinement step. , These are the position and velocity vectors of the space target, respectively.

[0072] Step 2: Set the dynamic model parameters for the initial orbit refinement model.

[0073] 2.1 The forces acting on a space target are divided into conservative and non-conservative forces. Conservative forces include the Earth's central gravity, non-spherical gravity, the gravitational forces of the Sun, Moon, and planets, Earth's tides, and general relativistic effects. Non-conservative forces include atmospheric drag, solar radiation pressure, Earth's radiation pressure, and mechanical forces. Calculating the perturbations acting on a space target involves obtaining the ephemeris data of the Sun, Moon, and planets, as well as selecting models such as the Earth's gravity field model, atmospheric mass density model, solar radiation pressure model, and ocean tides and solid tides. Furthermore, non-conservative forces are closely related to satellite velocity, attitude, and surface characteristics, requiring the determination of the atmospheric drag coefficient. (Typical values ​​range from 1.5 to 3.0), solar radiation pressure coefficient (Typical values ​​range from 1.2 to 1.9) and key parameters such as the surface-to-mass ratio of space targets. The surface-to-mass ratio... Cross-sectional area of ​​the space target With quality For the ratio, it is recommended to take 0.01 for unknown targets.

[0074] It should be noted that this invention is based on a unified dynamic equation and force model framework constructed using a geocentric inertial coordinate system. On this basis, the force model can be modified in a targeted manner for targets of different orbit types (such as weakening atmospheric drag in high orbits and adjusting solar radiation pressure in low orbits). It can also incorporate all perturbation force models as needed, thereby achieving accurate modeling of various orbit targets (covering different orbit types such as LEO, MEO, and GEO) while ensuring the consistency of dynamic physics. This provides a flexible and high-precision dynamic foundation for the subsequent initial orbit refinement process.

[0075] Step 3: Construct the state equations.

[0076] 3.1 Orbital state of a space target in an inertial coordinate system Change over time The following equations of motion can be used to represent it:

[0077] (1)

[0078] in, , , These are the position and velocity vectors of the space target, respectively. This is a white noise process used to compensate for random errors such as nonlinearity errors, simplification of dynamic and measurement models. Its expected value is... Covariance satisfy

[0079] (2)

[0080] (3)

[0081] in It is a positive definite matrix. It is Kronek function.

[0082] 3.2 Since the observations are conducted at discrete time points, it is necessary to discretize the motion equation (i.e., the state equation) in equation (1). Let the time span of the single-arc observation data be... seconds, each observation time is , , ..., In approximate values , If we perform a Taylor expansion and ignore higher-order terms, then we have:

[0083] (4)

[0084] , (5)

[0085] in , The equations of motion of equation (1) are respectively... The Jacobian matrix of the orbital state vector and process noise at each time step.

[0086] make , representing the update amount of the target orbit prediction state, then we have

[0087] (6)

[0088] The solution to the above equation is

[0089] (7)

[0090] In the formula The state transition matrix defines the epoch time. arrive The relationship between state vectors can be obtained by solving the following variational equations through numerical integration.

[0091] (8)

[0092] (9)

[0093] in, The identity matrix; the state transition matrix satisfy

[0094] (10).

[0095] To simplify the expression, let , Then equation (7) becomes

[0096] (11)

[0097] In the above formula The system noise in the state equation has a variance. Defined as

[0098] (12).

[0099] When the time interval between adjacent observation epochs When →0, it is available Moment Alternative Therefore, equation (12) above can be further simplified to:

[0100] (13).

[0101] Furthermore, target status at any time and state covariance It can be represented as:

[0102] (14)

[0103] (15).

[0104] It should be noted that a numerical integrator is used to solve the variational equations (8)-(9) to obtain the state transition matrix. In such cases, a high-order fixed-step or adaptive-step integrator is typically chosen. When the integrator step size is inconsistent with the observation epoch interval, the integrator is advanced to the observation epoch time. Afterwards, according to The orbital prediction state update is obtained by interpolating the integral state before and after time step. Its specific expression is shown in equation (11). Finally, based on equation (6), the prediction is obtained. orbital state at any given time and state covariance matrix (i.e., the prior state of filtering) and covariance matrix ).

[0105] Step 4: Construct the observation equation.

[0106] 4.1, Assuming single-arc observation data is in Observation vector at time If the observations consist solely of simultaneous angle measurements, or angle and distance measurements, then the initial orbital state... Satisfy the following functions:

[0107] (16)

[0108] In the formula The variance represents the measurement noise. It can be determined by the accuracy of the observed values, or it can be given directly based on experience; and All are the first The theoretical values ​​of the observed vectors are expressed as follows: Instantaneous orbital state at a given moment and initial orbital state The function.

[0109] 4.2 Approximate value of the residual vector in the initial orbital state The location can be approximated as

[0110]

[0111] (17)

[0112] in, , For model observation pairs The partial derivative Jacobian matrix of the orbital state vector at time step is given by:

[0113] , (18).

[0114] In the EKF algorithm, in order to update the state prediction results in step 3, it is necessary to change the measurement residual vector. Represented as time The function is used to calculate the result. The New Message of Time:

[0115] (19).

[0116] Step 5, filter initialization.

[0117] 5.1 Determine the white noise process based on the accuracy of the dynamic model. covariance matrix .

[0118] 5.2 Using the initial trajectory optimal solution The initial orbital state is used as the initial orbital state, and the initial orbital state covariance matrix is ​​set based on the initial orbital determination accuracy. It should be noted that the accuracy of the initial trajectory determination can be given based on experience, and more accurate initial trajectory state covariance information will accelerate the convergence speed of the filter.

[0119] 5.3 Initialize the measurement noise covariance matrix based on the accuracy of the observation data. .

[0120] Step 6: Use EKF technology to refine the initial trajectory and determine the results.

[0121] 6.1 Since the filtering process involves the prediction and updating of the orbital state, for ease of description, the superscript "-" indicates the prior orbital state and covariance information obtained from the state equation prediction, and the superscript "+" indicates the orbital information after updating the state and covariance based on the observation data.

[0122] 6.2 State Prediction. Based on formulas (4)-(6) in step 3.2, using... Time-based filtering improved state vector and state covariance matrix Forecast received Prior state vector at time t and state covariance matrix .

[0123] 6.3 State Update. First, the prior information calculated in step 6.2 is used. , and measurement noise covariance matrix Calculate Kalman gain As shown in equation (20). Next, use Update target orbit status Covariance Matrix For example, equations (21) and (22):

[0124] (20)

[0125] (twenty one)

[0126] (twenty two).

[0127] 6.5, the filter from and Start-up, given process noise Observation vector and measurement noise covariance ( Repeat steps 6.2-6.3 until all observation times have been covered, and the above techniques can be used to complete the refinement process of the initial trajectory of the unknown target. After trajectory refinement, the accuracy of the single-arc IOD result is significantly improved compared with the initial trajectory parameters before filtering.

[0128] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined to form new technical solutions. Such combinations are not bound by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0129] This invention uses single-station ground-based radar short-arc / ultra-short-arc observation data (station parameters are shown in Table 1) as an example to verify and explain the initial orbit refinement method. The observation data is generated based on SpaceX's publicly available Starlink satellite precise ephemeris simulation. The example uses a batch of Starlink satellites (NORAD IDs: 62032 to 62055, a total of 24 satellites) launched on November 21, 2024. The observation period is from November 28, 2024 to December 4, 2024 (7 days, observing 1023 observation arcs), with an orbital altitude spanning 300 km to 340 km. It should be noted that the satellite's orbital period is approximately 94.27 minutes. Therefore, this embodiment defines arcs within 60 seconds (less than 1% of the orbital period) as ultra-short-arc data, and assumes that the experimental object is a new target with completely unknown physical characteristics and orbital parameters. Within a single observation arc, the orbital altitude change of the experimental satellite is within 100 meters. Figure 2 The orbital altitude distribution of the experimental satellite at the simulated observation time is shown. Figure 3 This paper presents a typical filtered three-dimensional position error convergence effect diagram calculated based on 180s short-arc ground-based radar observation data, including orbit determination accuracy for ultra-short arcs (15s, 30s, 60s, 120s, and 180s), and transforms it to the RSW coordinate system. It should be noted that the simulation data generation process in this embodiment is only a means of verifying the technical effect and does not constitute the core technical solution claimed in this invention; this embodiment is only used to exemplarily explain the technical principles and implementation methods of this invention and does not limit the scope of protection of this invention.

[0130] Table 1 Simulation parameter settings for observed values

[0131] .

[0132] Figure 1This is a schematic diagram illustrating the high-precision initial trajectory determination process for unknown targets based on short-arc / ultra-short-arc observation data and filtering enhancement technology provided by this invention. The following will be combined with... Figure 1 The embodiments shown further illustrate the implementation steps and core innovations of the technical solution of the present invention.

[0133] Step 1: Determine the initial orbit based on short arc / ultra-short arc observation data.

[0134] 1.1 In this embodiment, single-station ground-based radar observation data is used, and the observation vectors are azimuth, elevation, and two-way range. The station (114.3°E, 30.5°N, altitude 39.3m) excluded observations with elevation angles less than 5° (see Table 1 "Cutting-out Elevation Angle"), and the time period was truncated. The experiment was conducted using arc segments of seconds. Each arc segment contained... One observation point, Taking 15, 30, 60, 120, and 180, the corresponding time intervals for the arc segments are: .

[0135] 1.2, in Ten observation points are taken at the beginning and end of the arc segment (first segment: End section: ), forming 100 observation pairs Calculate the initial and final slant distance vectors for each observation pair. and The Lambert equations (see Valladolid's 2007 monograph "Fundamentals of Astrodynamics and Applications") were used to obtain 100 sets of orbital elements. , ( .

[0136] 1.3 Calculate the observation residuals corresponding to each group of orbital elements, and select the group with the smallest sum of squared residuals as the optimal initial orbit solution: The results showed that the initial orbit determination was successful for all arc segments, achieving a 100% success rate for Initial Orbit Determination (IOD), and yielding a total of 1023 optimal initial orbit solutions. Furthermore, the optimal initial orbit solutions were converted into initial state vectors. ( .

[0137] 1.4 High-precision state vector extracted from SpaceX's precise ephemeris Assuming a true orbit, evaluate the initial orbit to determine the state vector error. ( And calculate the initial orbital state error RMS:

[0138] (twenty three)

[0139] (twenty four)

[0140] in, This represents the total number of arcs that were successfully used in IOD. ; and These represent the RMS values ​​for position and velocity errors, respectively.

[0141] Based on the IOD statistics of the 60s measurement arc segment, the position error is shown. Speed ​​error Therefore, this embodiment selects a diagonal matrix. The covariance matrix serves as the initial orbital state for filtering. It should be noted that the initial orbital accuracy of the IOD method varies across different time segments, and in practical applications, it can be flexibly adjusted according to the time span of the observed data. Furthermore, accurate initial orbital covariance information can accelerate filter convergence; this embodiment employs a unified [covariance matrix]. This is solely for verifying the effectiveness of the algorithm.

[0142] Step 2, setting the parameters of the filter dynamics model. The maximum magnitude of the non-conservative force perturbation experienced by the LEO satellite is... (See Montenbruck and Gill's 2000 monograph "Satellite Orbits: Models, Methods, and Applications"). Experiments show that during short-arc observation periods (e.g., 100s arc length), the impact of non-conservative force perturbations, including atmospheric drag, on satellite orbits is not significant, with position errors less than meters. Furthermore, non-conservative forces are affected by the physical parameters of the space target; for targets with completely unknown physical information, the accuracy of non-conservative forces is difficult to guarantee. Therefore, for Starlink satellites in LEO orbits, the filtering process mainly considers conservative forces, specifically including Earth's central gravity, Earth's non-spherical gravity, and three-body gravity. The force model used is a 10×10 order JGM-3 Earth gravity field model, and the influence of three-body gravitational perturbations is calculated using the DE406 planetary ephemeris.

[0143] Step 3, constructing the state equations.

[0144] 3.1 Once the dynamic model is determined, the equations of motion for the Starlink satellites in the inertial coordinate system can be obtained:

[0145] (1).

[0146] To compensate for nonlinear errors, errors caused by simplification of the dynamic model, etc., take ,in for Identity matrix (see Montenbruck and Gill's 2000 monograph "Satellite Orbits: Models, Methods, and Applications").

[0147] 3.2 Discretize the motion equation (1) and obtain the initial trajectory solution. , Taylor unfolds, and we get

[0148] (6)

[0149] in, , , .

[0150] 3.3 Furthermore, an 11th-order Cowell numerical integrator with a step size of 30s is used to solve the variational equations, and the state transition matrix is ​​derived. Thus obtain Update amount of orbital prediction state at time:

[0151] (11)

[0152] In the above formula, The system noise has a variance. The sampling interval of the observations is calculated using the following formula (13). s.

[0153] (13).

[0154] It should be noted that, due to the integrator step size and the observation sampling interval... Inconsistent, therefore the following definition ( Numerical integration time, Refers to the time of each observation epoch. To obtain the observation epoch... Update amount of orbit prediction status The integrator needs to be advanced to the observation epoch first. Then, according to Interpolate the integral states before and after. Based on Further calculations will yield the following results. orbital state at any given time and the corresponding state covariance matrix .in, , .

[0155] Step 4: Construction of observation equations.

[0156] 4.1, given Time-based radar observation vector ,but Time-observation and state approximation satisfy Among them, the measurement noise covariance matrix Determined by the measurement accuracy of simulated observations, i.e. ,in It is used to convert arcseconds into radians.

[0157] 4.2, Further Calculations Time-residual vector:

[0158] (19)

[0159] in, , indicating the observation pair The partial derivative of the orbital state vector at time step Jacobian matrix.

[0160] Step 5, from and To begin, the EKF algorithm is used to refine the initial trajectory and determine the result.

[0161] 5.1 State Prediction. Utilizing The updated state vector at each time step Covariance Matrix The forecast received Prior state vector at time t and prior covariance Among them, the state vector , Prior covariance of state vectors System noise covariance It is calculated from equation (13).

[0162] 5.2, State Update. Kalman Gain ,use Updated orbital state vector Covariance Matrix for:

[0163] (twenty one)

[0164] (twenty two)

[0165] In equation (22), It is an identity matrix.

[0166] 5.3 Repeat steps 5.1 and 5.2 until the filtering is completed at the last observation time, and obtain the final high-precision initial orbit.

[0167] Step 6: Using the precise ephemeris released by SpaceX as a reference, evaluate the filtered improved initial orbit (EKFIOD) state vector. ( and orbital elements The accuracy of the initial orbit is shown in Tables 2 and 3. Tables 2 and 3 present the statistical results of the orbital elements and external coincidence error of the orbital state before and after initial orbit refinement for short-arc / ultra-short-arc radar observation arcs.

[0168] The following conclusions can be drawn from the table:

[0169] (1) In the ultra-short arc scenario, the present invention overcomes the ill-conditioned problem of the normal equation in the BLS method and significantly improves the accuracy of the semi-major axis: in the 15s arc segment, the RMS of the semi-major axis error is improved by 40.2% compared with the original; in the short arc segment scenario, the accuracy of the semi-major axis after filtering is improved by about 10% compared with the accuracy of the IOD semi-major axis.

[0170] (2) Under various arc length observation conditions, the present invention significantly improves the tilt angle and eccentricity. Specifically, under the 120s arc segment condition, the tilt angle error RMS decreased from 47.76″ to 6.94″, an improvement of approximately 85.5%; the eccentricity error RMS decreased from... Down to This improved the accuracy by approximately 62.7%, effectively optimized the orbital geometry, and enabled accurate inversion of initial orbital parameters under short-arc observation conditions.

[0171] (3) At the orbital state level, the present invention has significant optimization effects on position and velocity. Under the 15s ultra-short arc observation condition, the position error RMS decreased from 0.64km to 0.24km, an improvement of 62.5%, and the velocity error RMS decreased from 67.06m / s to 28.21m / s, an improvement of about 57.9%. Under the 120s arc segment, the position and velocity errors RMS further decreased to 0.10km and 1.39m / s, respectively, an improvement of about 84.8% and 79.9% compared with the original, with a significant improvement in accuracy.

[0172] The above conclusions show that the present invention breaks through the accuracy limitation of determining the initial orbit of unknown targets in short arc / ultra-short arc scenarios, and achieves significant performance improvements in both orbital elements and state vectors, providing key technical support for high-precision space applications such as spatial target arc segment association, anomaly detection and identification, high-precision orbit cataloging, and disassembly tracing.

[0173] Based on the same inventive concept as the foregoing embodiments, this embodiment of the invention also provides a device for determining the high-precision initial trajectory of an unknown target with a short arc based on filter enhancement, including a memory and a processor. The memory stores program instructions that are executed by the processor, and the processor calls the program instructions to execute the method for determining the high-precision initial trajectory of an unknown target with a short arc based on filter enhancement.

[0174] Based on the same inventive concept as the foregoing embodiments, this embodiment of the invention also provides a non-transitory computer-readable storage medium storing computer instructions that cause the computer to execute the aforementioned method for determining the high-precision initial trajectory of an unknown target using short arcs based on enhanced filtering.

[0175] Table 2. Statistical results of the external coincidence error of the orbital elements for the initial orbit determination of 1023 short-arc / ultra-short-arc radar data.

[0176] .

[0177] Table 3. Statistical results of orbit state out-of-range coincidence error for initial orbit determination of 1023 short-arc / ultra-short-arc radar data.

[0178] .

[0179] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0180] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0181] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0182] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0183] In summary, this invention discloses a high-precision initial orbit determination method for unknown targets with short arcs based on enhanced filtering, belonging to the field of aerospace orbital dynamics and space target monitoring. This method can determine the initial orbit of near-Earth space targets using ground-based radar or space / ground-based optical monitoring equipment, employing high-success-rate initial orbit determination (IOD) methods such as the Gooding method and range search method. Then, extended Kalman filtering (EKF) technology is used to refine the initial orbit, thereby improving the accuracy of the initial orbit parameters.

[0184] This invention first adapts the initial orbit determination method to the data type of the observations. After obtaining multiple candidate initial orbit solutions, the optimal initial orbit solution is determined according to the principle of minimizing measurement residuals and converted into position and velocity vectors. Next, a dynamic model and parameters are configured based on the target orbit type and physical characteristics. A target state equation is established based on this model, while process noise is introduced to compensate for model errors. An observation equation is constructed using the position information of the observation platform. To ensure effective filter activation, parameters such as orbit state, state error covariance, process noise, and measurement noise need to be calculated and initialized in advance. Finally, by combining the state equation and the observation equation, the target orbit state is predicted and updated using EKF technology until all observation epochs are traversed, obtaining high-precision target initial orbit information.

[0185] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A method for determining the high-precision initial trajectory of an unknown target in a short arc based on enhanced filtering, characterized in that: include: Based on the type of observation data, a suitable method is used to determine the initial orbit, and the optimal initial orbit solution is determined from the multiple sets of initial orbit candidate solutions obtained. Configure a dynamic model based on the target orbit type and physical characteristics; Establish the target state equation based on the aforementioned dynamic model; The observation equations are constructed by combining observation data and observation platform location information, including: set up Observation vector at time It consists of angle measurement data, or angle measurement and distance measurement data at the same moment, satisfying ， in Indicates measurement noise. and These are the observation vectors with respect to the instantaneous state vector. and initial orbital state Theoretical value; residual vector Approximate value of the optimal solution for the initial trajectory The approximation is ， in, , , For model observation pairs The Jacobian matrix of the orbital state vector at time t=0 satisfies ， in, For observation epochs The approximate value of the orbital state at time t is used to update the state prediction results, and the measured residual vector is... Represented as time The function is calculated by the following formula. The New Message of Time: ； Combining the target state equation and the observation equation, starting from the optimal solution of the initial orbit and the initial orbit state covariance, the target orbit state is predicted and updated through extended Kalman filtering until all observation epochs are traversed, thus obtaining high-precision target initial orbit information.

2. The method for determining the high-precision initial trajectory of an unknown target using a short arc based on filter enhancement according to claim 1, characterized in that, When the observed data type is space-based optical observation data, the process of determining the optimal initial orbit solution includes: Observation data points are selected at the beginning and end of the arc segment, and then combined in pairs to form observation data point pairs; The direction of the slant range vector is determined based on the right ascension and declination in the inertial coordinate system. Based on the orbit type of the space target and the location information of the observation platform, the initial and final slant ranges are determined through large and small step distance searches. Then, the initial and final slant range vectors are calculated by combining the direction of the slant range vector and the initial and final slant ranges, thus constructing the Lambert problem. The Lambert problem was solved using the initial orbit determination algorithm, and multiple sets of candidate solutions for the initial orbit were obtained. The optimal solution for the initial orbit is determined based on the principle of minimizing the measurement residual, and then converted into the initial orbit state vector.

3. The method for determining the high-precision initial trajectory of an unknown target using a short arc based on filter enhancement according to claim 2, characterized in that, When the data type of the observation is ground-based optical observation data, the azimuth and elevation angles obtained from the ground-based observation in the station center coordinate system are converted into right ascension and declination in the inertial coordinate system, and then the optimal solution for the initial orbit is determined in the same way as for the space-based optical observation data.

4. The method for determining the high-precision initial trajectory of an unknown target using a short arc based on filter enhancement according to claim 1, characterized in that, When the data type of the observation is ground-based radar observation data, the process of determining the optimal initial orbit solution includes: Observation data points are selected at the beginning and end of the arc segment, and then combined in pairs to form observation data point pairs; Based on the station location, distance measurement, and angle measurement information, calculate the satellite's initial and final slant range vectors and construct the Lambert problem; The Lambert problem was solved using the initial orbit determination algorithm, and multiple sets of candidate solutions for the initial orbit were obtained. The optimal solution for the initial orbit is determined based on the principle of minimizing the measurement residual, and then converted into the initial orbit state vector.

5. The method for determining the high-precision initial trajectory of an unknown target using a short arc based on filter enhancement according to claim 1, characterized in that, The target state equation is established based on the aforementioned dynamic model, including: Based on the aforementioned dynamic model, the motion equations of the space target are constructed; The motion equations of the space target are discretized, and Taylor expansion is performed at the initial orbit optimal solution to obtain the discretized state equations.

6. The method for determining the high-precision initial trajectory of an unknown target using a short arc based on filter enhancement according to claim 1, characterized in that, Combining the target state equation and the observation equation, starting from the initial orbit optimal solution and the initial orbit state covariance, the target orbit state is predicted and updated using extended Kalman filtering, including: State prediction: Based on discretized state equations, utilizing... Time-based filtering improved state vector and state covariance matrix forecast Prior orbital state at time and prior state covariance matrix In this context, the superscript "-" indicates the prior orbital value calculated from the state equation, and the superscript "+" indicates the orbital state and covariance updated based on the observation data. State update: Utilizing prior state Prior state covariance matrix and measurement noise covariance matrix Calculate Kalman gain ,use Get the updated target state and the target state covariance matrix ; Filtering from and Start-up, given process noise Observation vector and measurement noise covariance matrix , Repeat the state prediction and state update steps until all observation epochs are traversed to obtain refined high-precision initial orbit information.

7. The method for determining the high-precision initial trajectory of an unknown target using a short arc based on filter enhancement according to claim 6, characterized in that, Before predicting and updating the target orbit state using the extended Kalman filter, the following filter initialization is also included: The covariance matrix of the white noise process is determined based on the accuracy of the dynamic model. The optimal solution of the initial orbit As the initial orbital state of the filter, the initial orbital state covariance matrix is ​​set based on the initial orbital determination accuracy. ; Initialize the measurement noise covariance matrix based on the accuracy of the observations. .

8. A device for determining the short-arc high-precision initial trajectory of an unknown target based on filter enhancement, characterized in that, The method includes a memory and a processor, wherein the memory stores program instructions that are executed by the processor, and the processor invokes the program instructions to execute the method for determining the high-precision initial trajectory of an unknown target based on filtering enhancement as described in any one of claims 1 to 7.

9. A non-transitory computer-readable storage medium, characterized in that, The non-transitory computer-readable storage medium stores computer instructions that cause the computer to execute the method for determining the high-precision initial trajectory of an unknown target using short arc based on filtering enhancement, as described in any one of claims 1 to 7.

Images