Satellite on-board orbit determination method, device and computer equipment

Abstract

The application relates to the field of satellite TT&C (Telemetry, Tracking and Command) technology and discloses a satellite on-orbit orbiting method, device and computer equipment. The method comprises the following steps: acquiring current observation data of a target satellite in a current orbiting arc segment; performing state recursion calculation based on the current observation data to obtain original ephemeris data; wherein the original ephemeris data comprises original prediction data and original compensation data of all observation epochs in the current orbiting arc segment; performing residual analysis based on the original ephemeris data and the current observation data to obtain design matrix data and residual vector data; jointly solving the design matrix data and the residual vector data by using a scaling factor to determine state correction data; and updating the original prediction data by using the state correction data to obtain target orbiting data. Finally, the satellite can continuously and stably realize real-time high-precision orbiting without relying on ground post-processing by using limited on-board computing resources.

Description

Technical Field

[0001] This application relates to the field of satellite telemetry, tracking and command technology, and in particular to a satellite on-board orbit determination method, apparatus and computer equipment. Background Technology

[0002] In the field of spacecraft intelligence (such as low-Earth orbit communication constellations, high-resolution remote sensing satellites, and scientific exploration satellites), precise orbit determination is one of the core technologies for autonomous satellite operation. Its accuracy directly affects the application effects of remote sensing image geometric positioning, satellite orbit control, scientific data inversion, and inter-satellite laser communication.

[0003] However, most low-Earth orbit (LEO) satellite precision orbit determination methods in related technologies rely on post-processing using ground systems, which is insufficient to meet the application requirements of real-time telemetry, tracking, and command (TT&C), autonomous navigation, and rapid response. Therefore, on-board calculation methods for real-time orbit determination have emerged, but these are still limited by computational resources. Thus, a satellite-based orbit determination method capable of achieving real-time, high-precision orbit determination within the limited computing power of onboard computers is needed. Summary of the Invention

[0004] This application aims to at least partially solve one of the technical problems in related technologies. To this end, this application proposes a satellite on-board orbit determination method, apparatus, and computer equipment. The main technical solutions adopted in this application include: In a first aspect, this application provides an on-board orbit determination method for a satellite, comprising: acquiring current observation data of the target satellite in the current orbit determination arc; performing state recursive calculation based on the current observation data to obtain original ephemeris data; wherein the original ephemeris data includes original prediction data and original compensation data for all observation epochs in the current orbit determination arc; performing residual analysis based on the original ephemeris data and the current observation data to obtain design matrix data and residual vector data; jointly solving the design matrix data and residual vector data using a scaling factor to determine state correction data; wherein the scaling factor is a diagonal scaling matrix used to eliminate the difference in numerical magnitude between parameters of different dimensions; and updating the original prediction data using the state correction data to obtain target orbit determination data.

[0005] The current observation data of the target satellite within the current orbit determination arc is acquired to provide raw input for subsequent calculations. State recursion is used to generate raw ephemeris data containing raw prediction and compensation data, establishing a precise mathematical relationship between the observed values ​​and the parameters to be estimated. By introducing a scaling factor, the numerical ill-conditioning caused by parameters of different dimensions is effectively resolved. Stable and efficient least-squares solutions are achieved under limited onboard computing power and precision, ultimately realizing continuous real-time orbit determination output. This approach enables the satellite to continuously and stably achieve real-time high-precision orbit determination using limited onboard computing resources, without relying on ground-based post-processing.

[0006] Optionally, the design matrix data and residual vector data are jointly solved using a scaling factor to determine the state correction data, including: adjusting the magnitude of the design matrix data using a scaling factor to obtain a scaled design matrix; and obtaining the state correction data by solving the least squares method based on the scaled design matrix and residual vector data.

[0007] By introducing a scaling factor based on physical priors to adjust the magnitude of the design matrix, the numerical ill-conditioning caused by differences in dimensions is fundamentally solved without increasing the computational load. This significantly improves the stability and efficiency of least squares solution in a spaceborne environment with limited floating-point precision, while ensuring the physical validity of the solution results. It also adapts to the limited computing power constraints of spaceborne computers and provides algorithmic support for achieving real-time high-precision orbit determination on spacecraft.

[0008] Optionally, based on the scaled design matrix and residual vector data, the state correction data is obtained by solving the least squares method, including: constructing a system of linear equations based on the scaled design matrix and residual vector data and solving the least squares method to obtain the scaled correction data; and restoring the scaled correction data using the inverse transformation of the scaling factor to obtain the state correction data.

[0009] By using inverse transformation, the normalized intermediate correction values ​​are restored to state correction data with real physical meaning. This not only retains the numerical stability advantage brought by scaling, but also ensures the physical validity of the correction values, which can be directly used to update the satellite's orbital state.

[0010] Optionally, residual analysis is performed based on the original ephemeris data and the current observation data to obtain design matrix data and residual vector data. This includes: using the current observation data and the original prediction data for quality screening to determine the effective solution information of the effective observation epochs in the current orbit determination arc; wherein the effective solution information includes effective error data and effective residual data; and performing data assembly processing based on the effective solution information to obtain design matrix data and residual vector data.

[0011] Abnormal observation data was removed through quality screening, ensuring the reliability of the input data. By assembling the data to construct a standardized design matrix and residual vector, a precise linear relationship was established between the observation residuals and the estimated state correction, laying a solid data foundation for subsequent scaled least squares solutions and effectively improving the accuracy and stability of the orbit determination results.

[0012] Optionally, quality screening is performed using current observation data and raw prediction data to determine the effective solution information for valid observation epochs in the current orbit determination arc. This includes: for any observation epoch in the current orbit determination arc, calculating the difference between current observation data and raw prediction data to obtain position residual data; using the position residual data and a preset threshold to determine validity, obtaining a valid determination result; if the valid determination result indicates that the current observation epoch is a valid observation epoch, the raw compensation data of the valid observation epoch is used as the valid error data, and the position residual data of the valid observation epoch is used as the valid residual data.

[0013] By calculating the position residual and using a dynamic threshold to remove outliers, abnormal observation points caused by multipath effects, receiver noise, or signal obstruction can be effectively removed, avoiding the contamination of orbit determination accuracy by outliers. At the same time, high-quality observation data is retained for subsequent calculations, significantly improving the reliability of orbit determination results.

[0014] Optionally, data assembly processing is performed based on the effective solution information to obtain design matrix data and residual vector data, including: constructing an observation matrix based on the effective error data for each effective observation epoch; stacking the observation matrices of each effective observation epoch to obtain design matrix data; and stacking the effective residual data of each effective observation epoch to obtain residual vector data.

[0015] By constructing the observation matrix for each valid observation epoch and stacking it, the observation information scattered across various observation epochs can be integrated into a unified system of linear equations, providing a standard input data structure for subsequent least squares solutions and ensuring the standardization and efficiency of the solution process.

[0016] Optionally, state recursion calculations are performed based on the current observation data to obtain the original ephemeris data, including: determining the initial state vector in the current observation data; performing numerical integration based on the initial state vector and the preset dynamic model to generate the original ephemeris data; wherein, the preset dynamic model includes constraints on Earth's center of mass gravitational perturbation, Earth's non-spherical gravitational perturbation, atmospheric drag perturbation, solar radiation pressure perturbation, third-body gravitational perturbation, Earth's tidal deformation perturbation, and general relativistic effects.

[0017] By determining a reasonable initial state vector from the current observation data, an accurate starting point for numerical integration is provided. By constructing a complete dynamic model that includes constraints such as Earth's center of mass gravitational perturbation, Earth's non-spherical gravitational perturbation, atmospheric drag perturbation, solar radiation pressure perturbation, third-body gravitational perturbation, Earth's tidal deformation perturbation, and general relativistic effects, all major perturbations affecting the orbital motion of low-Earth orbit satellites are comprehensively covered, significantly improving the accuracy and reliability of orbit prediction and laying a solid physical foundation for achieving real-time high-precision orbit determination on satellites.

[0018] Optionally, the original compensation data includes state transition data and parameter sensitivity data; numerical integration is performed in the following manner: based on the initial state vector and the preset dynamic model, the integration is performed from the starting point to the ending point of the current orbital segment, and the state equation and variational equation are solved in parallel to obtain the original prediction data, state transition data and parameter sensitivity data for each observation epoch.

[0019] By solving the state equations and variational equations in parallel, raw prediction data and raw compensation data for all observation epochs can be generated simultaneously during a single integration process. This avoids the increased computational load caused by repeated integrations, significantly improving computational efficiency. Simultaneously, it ensures the accuracy of the state transition data and parameter sensitivity data, providing accurate error propagation relationships for subsequent residual analysis and least squares solutions. Integrating the raw prediction data, state transition data, and parameter sensitivity data for each observation epoch constitutes complete raw ephemeris data, providing comprehensive data support for subsequent orbital state corrections.

[0020] Secondly, this application provides a satellite on-board orbit determination device, the device comprising: The data acquisition module is used to acquire the current observation data of the target satellite in the current orbit determination arc. The state recursion module is used to perform state recursion calculations based on the current observation data to obtain the original ephemeris data; the original ephemeris data includes the original prediction data and original compensation data of all observation epochs in the current orbit determination arc. The residual analysis module is used to perform residual analysis based on the original ephemeris data and the current observation data to obtain the design matrix data and residual vector data; The joint settlement module is used to jointly solve the design matrix data and residual vector data using a scaling factor to determine the state correction data; where the scaling factor is a diagonal scaling matrix used to eliminate the difference in numerical magnitude between parameters of different dimensions. The orbit determination module is used to update the original prediction data using state correction data to obtain the target orbit determination data.

[0021] Thirdly, this application also provides a computer device including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps of any of the above methods. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the specific embodiments of this application or the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0023] Figure 1 This is a flowchart of a satellite on-board orbit determination method according to an embodiment of this application; Figure 2 This is a flowchart of a method for determining ephemeris data according to an embodiment of this application; Figure 3 This is a flowchart of a method for determining residual data according to an embodiment of this application; Figure 4 This is a flowchart of a method for determining correction data according to an embodiment of this application; Figure 5 This is a structural block diagram of a satellite on-board orbit determination device according to an embodiment of this application; Figure 6 This is an internal structural diagram of a computer device provided according to an embodiment of this application. Detailed Implementation

[0024] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0025] It should be noted that, taking the dynamic method for precise satellite orbit determination as an example, this type of method processes observation data acquired by the onboard Global Navigation Satellite System (GNSS) receiver and combines it with a high-precision dynamic model to estimate the satellite's state. However, this implementation method typically relies on ground systems for post-processing. The processes of downlinking observation data, ground processing, and orbital injection can lead to delays of several hours or even days, making it impossible to meet real-time telemetry and control requirements. Furthermore, this implementation method is highly dependent on ground intervention and lacks autonomous anomaly detection, data preprocessing, and state recovery capabilities.

[0026] Furthermore, although onboard real-time orbit determination calculation methods have emerged in related technologies, they still have many limitations due to the strict constraints of onboard computing and storage resources.

[0027] First, spaceborne computers typically use aerospace-grade processors, whose computing power is far lower than that of ground-based workstations. Traditional batch processing algorithms are highly complex and require numerous iterations, making them difficult to run in real-time on satellite, thus posing a limitation on computing resources. Second, spaceborne memory capacity is limited, making it impossible to store long-term GNSS observation data. Orbit determination calculations must be completed within a limited data window, again highlighting the limitation on storage resources. Furthermore, taking low-Earth orbit satellites as an example, they are affected by various perturbations such as Earth's non-spherical gravity, atmospheric drag, the gravitational pull of the Sun, Moon, and other third bodies, solar radiation pressure, and relativistic effects. The computational load for accurate modeling is extremely large, further increasing the difficulty of real-time implementation on satellite.

[0028] Furthermore, the parameters to be estimated in dynamic orbit determination include the satellite's position, velocity, and atmospheric drag coefficient. The magnitudes of these parameters vary greatly, making it easy for the direct solution method to develop numerical ill-conditioning. Especially in an environment where the precision of onboard floating-point operations is limited, the algorithm struggles to converge stably.

[0029] Based on this, this application provides an embodiment of a satellite on-board orbit determination method, which can achieve rapid convergence, numerical stability, and high precision real-time orbit determination under conditions of limited onboard resources.

[0030] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0031] This embodiment provides a satellite on-board orbit determination method, such as... Figure 1 As shown, the method includes the following steps: S110: Obtain the current observation data of the target satellite in the current orbit determination arc.

[0032] The target satellite can refer to a spacecraft that requires real-time onboard orbit determination, such as a low-Earth orbit observation satellite, a communication satellite, or a scientific exploration satellite. The current orbit determination arc can refer to a continuous time interval selected for a batch orbit determination. All observation data generated by the target satellite within this time interval will be integrated into a single data unit for use in this orbit calculation.

[0033] It should be noted that the current orbit determination arc can consist of multiple consecutive observation epochs. That is, the current orbit determination arc is a continuous time range, while the observation epochs are different discrete sampling time points within that time range. Each observation epoch corresponds to an independent set of GNSS observation data. For example, if the current orbit determination arc is a continuous period of 60 minutes, assuming each observation epoch is spaced 30 seconds apart, then the arc contains a total of 120 observation epochs, and each observation epoch corresponds to a specific timestamp and a set of corresponding current observation data.

[0034] The current observation data can refer to raw measurement data that describes the target satellite's time information, spatial position, and motion state at each observation epoch. For example, the current observation data may include the target satellite's time data, position components, and velocity components at the current observation epoch. The time data can be a timestamp, such as a weekday second or Julian day in the GPS time system; the position components can be x, y, and z coordinate data in the Earth-fixed coordinate system; and the velocity components can be Vx, Vy, and Vz components in the Earth-fixed coordinate system.

[0035] Specifically, the raw observation data of the target satellite is first obtained, and then the raw observation data is processed to remove abnormal data and handle missing data in order to obtain the current observation data.

[0036] For example, raw observation data output by a spaceborne GNSS receiver can be obtained by reading a pre-recorded GNSS file (such as a common plain text table CSV file) along the target path. This GNSS file contains fields such as the timestamp, position component, and velocity component of the target satellite at each observation epoch, output by the spaceborne GNSS receiver.

[0037] During processing, the timestamps in the raw observation data are first parsed line by line, and then converted into a unified Coordinated Universal Time (UTC) date and time object according to a pre-determined time format (such as the time format specified in the configuration file during parsing). If the parsing of a timestamp fails, the line is directly discarded and an exception is logged.

[0038] Next, the integrity of each observation epoch is checked. This can be done line by line in each timestamp row, verifying the existence and completeness of the position and velocity components in the raw observation data at each observation epoch (time point). If the position and velocity components of the current observation epoch are incomplete (e.g., any field is missing), the observation epoch is marked as a missing observation epoch, and the corresponding missing data processing is performed.

[0039] Specifically, the missing epoch handling strategy can be as follows: If the number of missing observation epochs in the current orbit determination arc is within the first range (e.g., less than 5% of the total observation epochs in the current orbit determination arc), linear interpolation or cubic spline interpolation can be used to fill in the position and velocity components marked as missing observation epochs. If the number of missing observation epochs in the current orbit determination arc is within the second range (e.g., more than 5% but less than 50% of the total observation epochs in the current orbit determination arc), the processing of these missing observation epochs can be skipped. If the number of missing observation epochs in the current orbit determination arc is within the third range (e.g., more than 50% of the total observation epochs in the current orbit determination arc), the orbit determination processing for the current orbit determination arc is abandoned, and the current orbit determination process is terminated.

[0040] Next, the data obtained after missing data processing undergoes column deduplication. Specifically, if multiple observation data exist for the same timestamp, only the first valid record in the log is retained, resolving the issue of duplicate data generated during receiver data buffering or transmission. The deduplicated data is then sorted in ascending order by UTC timestamp to ensure that the entire data sequence monotonically increases in time, providing a correct time reference for subsequent numerical integration. Subsequently, the date and time objects under the UTC time scale are uniformly converted to a high-precision time representation under the GPS time system (such as GPS week seconds or Julian days), and the position and velocity components are saved as numerical arrays, ultimately yielding the standardized current observation data for the current orbital segment.

[0041] Furthermore, since the raw observation data output by the spaceborne GNSS receiver only contains position and velocity information, but the dynamic model also contains non-conservative force factors that significantly affect the long-term orbital evolution of the satellite, the current observation data can also include an influence coefficient to quantify the impact of this factor. This influence coefficient can refer to an empirical parameter characterizing the strength of the non-conservative perturbation force, and can also participate in subsequent state deduction and correction as a state variable to be estimated. For example, taking a low-Earth orbit satellite as the target satellite, the influence coefficient can be the atmospheric drag coefficient (Cd), which can be determined by the prior Cd value provided in the configuration file, and together they serve as a component of the current observation data. That is, the current observation data includes a position component, a velocity component, and an atmospheric drag coefficient component.

[0042] S120. Perform state recursion calculations based on the current observation data to obtain the original ephemeris data.

[0043] The raw ephemeris data can refer to a complete dataset that describes the expected state of the target satellite at each observation epoch. Specifically, it can include the raw prediction data and raw compensation data for all observation epochs in the current orbit determination arc. The raw prediction data can refer to the theoretical predicted position and theoretical predicted velocity of the target satellite at each observation epoch, obtained based on the integration of the dynamic model.

[0044] It is important to understand that due to inherent errors in dynamic models and the difficulty in accurately modeling non-conservative perturbations such as atmospheric drag, theoretically calculated predictions often deviate from subsequent actual observations of the target satellite. Therefore, to meet the requirements of high-precision orbit determination, the raw ephemeris data can also include raw compensation data. This raw compensation data can serve as auxiliary data for propagating state deviations at any observation epoch to the corresponding time series, and can quantify the impact of positional perturbations and changes in dynamic parameters on the subsequent state of the target satellite.

[0045] For example, the original compensation data may include state transition data and parameter sensitivity data. State transition data can refer to matrix data describing the impact of initial position disturbances on the satellite state at any subsequent time. This can be understood as a satellite orbit deviation conversion table. If the target satellite has an initial error at the beginning of the current orbit determination arc, and this initial error changes continuously as the target satellite orbits the Earth, the state transition data can accurately calculate how much error this initial error will become at any observation epoch in the current observation arc. This unifies the state deviations of all dispersed observation epochs within the current observation arc with the estimated state at the end of the current observation arc, allowing the satellite state at any time in the entire arc to be deduced simply by solving for the state correction at the end of the current observation arc.

[0046] Parameter sensitivity data can refer to matrix data showing the impact of dynamic parameters (such as atmospheric drag coefficient) on the satellite's state at any given time. Taking a low-Earth orbit satellite as an example, parameter sensitivity data can specifically refer to matrix data showing the impact of changes in atmospheric drag coefficient on the satellite's state, which can quantify the changes in satellite position and velocity caused by each unit change in atmospheric drag coefficient.

[0047] Specifically, the initial state vector is first determined based on the current observation data, and numerical integration is performed based on the initial state vector and the preset dynamic model. The integration is performed from the starting point to the ending point of the current orbital arc, and the state equation and variational equation are solved in parallel to generate the original prediction data, state transition data and parameter sensitivity data for each observation epoch. The state transition data and parameter sensitivity data are used as the original compensation data, which together with the original prediction data constitute the original ephemeris data.

[0048] S130. Perform residual analysis based on the original ephemeris data and current observation data to obtain the design matrix data and residual vector data.

[0049] The design matrix data can refer to a matrix formed by stacking the observation matrices of each valid observation epoch, i.e., a coefficient matrix used to describe the linearized relationship between the state correction and the observation residuals. The residual vector data can refer to a vector formed by stacking the position residuals of each valid observation epoch, i.e., the difference between the observed position in the current observation data and the estimated position in the original prediction data.

[0050] It should be noted that, due to the influence of transmission status, the original prediction data may contain meaningless or erroneous parameters. Directly using these parameters for orbit determination correction can lead to error accumulation, slow convergence, or even divergence. To address this issue, residual analysis can be performed between the original prediction data and the current observation data to eliminate outlier observation epochs containing abnormal data. This allows for the construction of a design matrix and residual vector that reflect the information of valid observation epochs, providing numerically stable input for subsequent scaled least squares solutions, thereby significantly improving the convergence speed and accuracy of orbit determination.

[0051] Specifically, firstly, quality screening can be performed using current observation data and raw prediction data to determine the effective solution information for valid observation epochs in the current orbit determination arc. This effective solution information can include effective error data and effective residual data. Then, data assembly processing is performed based on the effective solution information; that is, for each valid observation epoch, an observation matrix is ​​constructed based on its effective error data. Finally, all observation matrices are stacked to obtain the design matrix data, and all effective residual data are stacked to obtain the residual vector data.

[0052] S140. Use scaling factors to jointly solve the design matrix data and residual vector data to determine the state correction data.

[0053] As mentioned earlier, spaceborne computers often employ aerospace-grade processors, which suffer from limited computing and storage resources. Therefore, to significantly reduce the computational burden, a scaling factor can be introduced to scale parameters of different dimensions in the design matrix to the same order of magnitude. This unifies the numerical magnitude without increasing computational load, fundamentally reducing the number of normal equations required for subsequent integration and improving the stability and efficiency of numerical solutions.

[0054] Understandably, the scaling factor could be a diagonal scaling matrix used to eliminate differences in numerical magnitude between parameters of different dimensions.

[0055] Specifically, the orbit determination process may involve the target satellite's three-dimensional position components (on the order of approximately 1e6 meters), three-dimensional velocity components (on the order of approximately 1e3 meters per second), and atmospheric drag coefficient Cd (on the order of approximately 1). The magnitudes of these components can differ by a factor of 1e6, resulting in significant differences in the magnitudes of the columns in the original design matrix data. Directly constructing and solving the normal equations would lead to severe numerical ill-conditioning. To address this issue, a diagonal scaling matrix is ​​introduced as a scaling factor, resembling an identity matrix, but its diagonal elements are taken as the reciprocals of the magnitudes of each component (or directly taking typical magnitude values ​​and then processing them through inverse operations). That is, each column of the design matrix data is multiplied by the corresponding scaling factor, normalizing the column vectors corresponding to position, velocity, and atmospheric drag coefficient to similar numerical ranges. This avoids numerical ill-conditioning during the least squares solution process, thereby enabling stable solution of the state correction data under the limited precision environment of spaceborne computing.

[0056] The state correction data can refer to the correction data used to update the original forecast data, which represents the specific adjustment to the current state vector to be estimated (including position, velocity and atmospheric drag coefficient).

[0057] Specifically, the design matrix data is scaled using a scaling factor to obtain a scaled design matrix. Based on the scaled design matrix and residual vector data, a system of linear equations is constructed and solved using least squares to obtain the scaled corrected data. Then, the scaled corrected data is restored using the inverse transformation of the scaling factor to obtain state corrected data with the original physical dimensions.

[0058] Finally, the validity of the obtained state correction data is checked to ensure that all elements are finite values ​​and the correction norm is within a reasonable range before it can be used as valid input for subsequent update steps.

[0059] S150. Update the original prediction data using the state correction data to obtain the target orbit determination data.

[0060] Among them, the target orbit determination data can refer to the precise orbit state obtained after iterative correction, that is, the precise position, precise velocity and estimated atmospheric drag coefficient of each observation epoch of the current orbit determination arc segment.

[0061] Specifically, the position correction, velocity correction, and atmospheric drag coefficient correction in the state correction data can be superimposed on the corresponding components of the original prediction data, and the atmospheric drag coefficient can be processed to ensure that it is within a reasonable range, so as to finally obtain the target orbit determination data for the current orbit determination arc.

[0062] Furthermore, after obtaining the target orbit determination data for the current orbit determination arc, in order to achieve continuous real-time orbit determination on the satellite, a sliding window batch processing method can be adopted to repeatedly execute all the above orbit determination processes on a continuously sliding data window.

[0063] Specifically, the length of the sliding window and the sliding step size parameters are configured according to the satellite orbital altitude, the computing power of the onboard computer and the orbit determination accuracy requirements. An initial value transfer mechanism between the old and new orbit determination arcs is established. The orbit determination process is executed cyclically on the continuously received GNSS observation data stream to achieve real-time continuous output of the orbit determination results.

[0064] For example, taking a low-Earth orbit satellite with an altitude of approximately 500 km as an example, the orbit determination window length is set to 60 minutes, and the sliding step size is 10 minutes. The initial value of the state to be estimated for the new arc segment adopts the endpoint state of the orbit determination result of the previous arc segment, that is, the precise position, velocity, and estimated atmospheric drag coefficient after the iteration convergence of the previous window. This initial value transfer mechanism ensures that the initial state of the new arc segment is close to the true value, and convergence is usually achieved in only 2 to 3 iterations.

[0065] The continuous orbit determination process is as follows: In the initial stage, after the system starts up, it waits to accumulate the first complete 60-minute window of GNSS observation data. After the data accumulation is completed, the complete orbit determination process is executed to complete the first orbit determination.

[0066] During the stabilization phase, every 10 minutes the window slides, the system appends the newly received 10 minutes of observation data to the end of the window, while discarding old data at the beginning of the window that exceeds 60 minutes in length, forming a new orbit determination arc. The orbit determination result of the previous arc is used as the initial value to perform orbit determination updates.

[0067] In the output phase, after each orbit determination is completed, the precise orbital state of the arc segment endpoint and the estimated atmospheric drag coefficient are output. To ensure real-time performance, the calculation time for a single orbit determination can be controlled within 5 minutes. This is achieved through measures such as reducing the number of iterations by initial value propagation, improving solution efficiency by scaling least squares, and accelerating state recursion through parallel integration, all of which work together to ensure time constraints are met.

[0068] In the above implementation, current observation data of the target satellite within the current orbit determination arc is acquired to provide raw input for subsequent calculations. State recursion calculations are used to generate raw ephemeris data containing raw prediction data and raw compensation data, establishing a precise mathematical relationship between the observed values ​​and the parameters to be estimated. By introducing a scaling factor, the numerical ill-conditioning caused by parameters of different dimensions is effectively resolved. Stable and efficient least-squares solutions are achieved under limited onboard computing power and precision, ultimately realizing continuous real-time orbit determination output. This scheme enables the satellite to continuously and stably achieve real-time high-precision orbit determination using limited onboard computing resources, without relying on ground-based post-processing.

[0069] In some implementation methods, please refer to the appendix. Figure 2 Based on the current observation data, a state recursive calculation is performed to obtain the original ephemeris data, including: S210. Determine the initial state vector from the current observation data.

[0070] It should be noted that dynamic integration is a recursive calculation process, which requires starting from a known reference time and using the integral dynamic equation to predict the satellite state at each time point (observation epoch) within the entire prediction interval (current orbit determination arc). Therefore, it is necessary to first determine the state vector at that reference time as the starting point for integration.

[0071] Here, the initial state vector can refer to the set of the target satellite's position, velocity, and estimated dynamic parameters at the start of the current orbit determination arc. For example, it can be a seven-dimensional vector of the form x = [rx, ry, rz, vx, vy, vz, cd]ᵀ, containing the three-dimensional position components (rx, ry, rz), three-dimensional velocity components (vx, vy, vz), and atmospheric drag coefficient (cd) at the start of the orbit determination arc.

[0072] Specifically, since this scheme estimates the state at the end of the current orbit determination arc, the position and velocity components of the last epoch in the preprocessed current observation data are taken and combined with the prior atmospheric drag coefficient provided in the configuration file to form a 7-dimensional initial state vector. Although the direction of numerical integration is from the start to the end of the arc, the state transition data can be used to uniformly correlate the state deviations of each observation epoch with the end of the arc. Therefore, using the last epoch as the initial value simplifies the subsequent construction of the normal equations and improves computational efficiency.

[0073] S220. Numerical integration is performed based on the initial state vector and the preset dynamic model to generate the original ephemeris data.

[0074] The preset dynamic model can refer to a mathematical model that describes the relationship between the sum of accelerations of various forces acting on the target satellite in the inertial frame and the satellite's state, i.e., a set of equations used to calculate the satellite's acceleration at any given time.

[0075] Specifically, the pre-defined dynamic model includes constraints on gravitational perturbation by the Earth's center of mass, gravitational perturbation by the Earth's non-spherical shape, atmospheric drag perturbation, solar radiation pressure perturbation, gravitational perturbation by the third body, tidal deformation perturbation, and general relativistic effects.

[0076] For example, the preset dynamic model can be represented by the following equation: In the formula, The acceleration vector representing the target satellite; The position vector of the target satellite; The distance from the target satellite to the Earth's center; The gravitational constant of Earth; For non-spherical perturbation acceleration of the Earth; For the perturbation force of the three bodies; Atmospheric drag; Solar radiation pressure; For other perturbation terms; m is the mass of the target satellite.

[0077] Regarding the constraint of gravitational perturbation by the Earth's center of mass, which serves as the fundamental perturbation force for precise satellite orbit determination, accurate modeling directly determines the baseline accuracy of orbit determination. For example, it can be defined as: In the formula, The position vector of the target satellite; The distance from the target satellite to the Earth's center; is the Earth's gravitational constant.

[0078] To address the constraint of Earth's non-spherical gravitational perturbation, Newton's law of universal gravitation simplifies Earth as a homogeneous sphere. However, Earth is actually an oblate spheroid, flattened in the north and south and wider in the middle. This causes significant variations in the gravitational force exerted on a satellite at different latitudes and longitudes, leading to substantial errors if Newton's law of universal gravitation is applied. Therefore, Earth's gravitational potential U can be defined as follows: In the formula, This represents the Earth's reference radius, which is approximately 6378.137 km. and These are the geocentric latitude and longitude of the target satellite, respectively. and These are the correlation harmonic coefficients, which are related to the Earth's mass distribution. It is an nth-order m-degree associated Legendre polynomial, defined as follows: in Let be an nth-order Legendre polynomial, defined as follows: Furthermore, due to the unknown distribution of mass within the Earth, and It cannot be calculated, but it can be obtained through satellite measurement. , , , Once matched, the satellite's non-spherical perturbation acceleration vector can be calculated using the following formula: in, This represents the gradient of Earth's gravitational potential.

[0079] For atmospheric drag perturbation constraints, atmospheric drag is a dissipative perturbation force generated by the interaction between the atmosphere and the satellite. Its magnitude is mainly related to atmospheric density, the relative velocity between the target satellite and the atmosphere, and the shape of the satellite, and its impact is more significant on low-Earth orbit spacecraft. Its expression is: In the formula, The aerodynamic coefficient, Where S is the atmospheric density and S is the windward area. The relative velocity between the target satellite and the atmosphere. for The model, using the static atmosphere assumption, is expressed as follows: in, The ground velocity vector of the target satellite in the J2000 coordinate system. This represents the Earth's rotational angular velocity vector in the J2000 coordinate system. This is the position vector of the target satellite in the J2000 coordinate system.

[0080] Regarding the constraint of solar radiation pressure perturbation, also known as solar radiation perturbation, it is mainly generated by the radiation of photons and other particles emitted by the sun onto the satellite surface, and its expression is: in, As the shadow factor, The solar radiation pressure constant near Earth. The surface reflectance of the target satellite is typically taken as 0.5. For the rate of change over time, The time interval from the start time, For the target satellite cross-section, This is the vector from the sun to the target satellite. The magnitude of solar radiation pressure perturbation is typically in the range of 10. -8 about.

[0081] Regarding the constraint of gravitational perturbation by a third body, since the satellite is subject not only to the Earth's gravity but also to the gravity of other celestial bodies such as the Sun and Moon, the specific expression for the three-body perturbation force is as follows: in, For the mass of the third body, For the mass of the target satellite, This is the position vector of the third body in the current coordinate system. This represents the position vector of the target satellite. This module considers low-Earth orbit spacecraft, primarily taking into account the three-body perturbation of the Sun and Moon.

[0082] Regarding the constraint of Earth's tidal deformation perturbation, since Earth is a non-rigid, elastic sphere, its shape and mass distribution are affected by the gravitational pull of other celestial bodies such as the Sun and Moon, resulting in changes. This small-scale, periodic deformation of the solid Earth is called solid tides. In addition, the gravitational pull of the Sun and Moon causes ocean tides, thus altering Earth's gravitational potential. Both solid tides and ocean tides cause changes in Earth's gravitational potential, thereby affecting the motion of target satellites.

[0083] Therefore, in precise satellite orbit determination, it is necessary to consider the effect of tides on the acceleration of the satellite. The change in Earth's gravitational potential caused by tides can be expressed as the change in the coefficients of the spherical harmonic function of Earth's gravitational field, which can be described by correcting the change in the coefficients of the spherical harmonic function of Earth's gravitational potential.

[0084] The formula for calculating the correction of the spherical harmonic coefficient caused by tidal perturbations of the solid Earth is shown below: in, The Love coefficients are of order n. , These are the latitude and longitude of the Earth-fixed system of the perturbed celestial body, respectively. , Here, represents the gravitational constant of the perturbed celestial body and its distance from the Earth's center, respectively. Tidal perturbations of the solid Earth cause a decrease in satellite acceleration by at least a quarter order of magnitude. The amplitude of the influence of ocean tides is an order of magnitude smaller than that of solid Earth tides, and the formula for calculating its spherical harmonic coefficient correction is shown below: in, , These are the seawater density and the load deformation coefficient, respectively. , , and These are the tidal coefficients of the tidal component S, respectively. The sum is a weighted sum of six Doodson variables, which are the fundamental arguments of the Sun and Moon's orbits and are closely related to the arguments of the nutation sequence.

[0085] To address the constraints imposed by general relativistic effects, satellite orbit determination should be solved using general relativistic theory, as defined below: in, and These represent the satellite's position and velocity unit vectors, respectively. c It is the speed of light in a vacuum.

[0086] After constructing the pre-defined dynamic model, it can be integrated. Specifically, numerical integration can be performed in the following ways: Based on the initial state vector and the preset dynamic model, the state equation and variational equation are solved in parallel from the starting point to the ending point of the current orbit segment to obtain the original prediction data, state transition data and parameter sensitivity data for each observation epoch.

[0087] It should be noted that since this method estimates the state at the end of the arc segment, and the integration direction is from the starting point to the end point, it is necessary to calculate the state transition matrix and the parameter sensitivity matrix simultaneously during the integration process. The state transition matrix is ​​used to linearly propagate the state deviations from each observation epoch to the end of the arc segment, and the parameter sensitivity matrix is ​​used to describe the influence of changes in the atmospheric drag coefficient on the state. Therefore, it is necessary to construct the state equation and the variational equation.

[0088] The state equation can be defined as: In the formula, Represents the target satellite at time The state vector (including position, velocity, and atmospheric drag coefficient), i.e., the target satellite's state vector at time [time value missing]. The state vector to be estimated.

[0089] The state transition data for each observation epoch can then be determined by the following variational equation: In the formula, A is the system Jacobian matrix of the state equation with respect to the state vector, and we have f is the state derivative; Used to propagate the initial state error to the observation time. This refers to the state transition data, used to transfer the initial time... The state error propagates to the observation time ; I It is an identity matrix.

[0090] The parameter sensitivity data for each observation epoch can also be determined by the following equation: In the formula, is parameter sensitivity data, describing the partial derivative of the state with respect to the atmospheric drag coefficient; B represents the orbital state correction matrix, which is the partial derivative matrix of the state derivative with respect to the atmospheric drag coefficient, used to describe the degree of influence of changes in the atmospheric drag coefficient on the satellite state.

[0091] Specifically, during the numerical integration process, starting from the beginning of the current orbit determination arc, a high-precision numerical integrator (such as the variable-step Runge-Kutta method) is used to integrate the state equation and variational equation in parallel. Within each integration step, the satellite state, state transition matrix, and parameter sensitivity matrix are updated simultaneously until the integration reaches the end of the arc. After integration is completed, the raw prediction data (position, velocity), state transition data, and parameter sensitivity data for each observation epoch (i.e., each discrete output time) are obtained.

[0092] For example, state transition data Satisfying the above variational equation, at the initial time The state transition matrix is ​​the identity matrix; parameter sensitivity data Satisfying the above differential equation, at the initial time The parameter sensitivity matrix is ​​a zero matrix. By performing parallel numerical integration on these two differential equations, the state transition data and parameter sensitivity data corresponding to each observation epoch can be obtained. That is, assuming there are m observation epochs in the current orbit determination arc, with times as follows: (in As the starting point of the arc segment, (For the endpoint of the arc segment). Through parallel integration, the following data can be generated for each epoch: raw prediction data Each epoch State transition data: (i.e., from) To the finish line (state transition matrix) and each epoch Parameter sensitivity data: (That is, the partial derivative of the atmospheric drag coefficient with respect to the state of that epoch).

[0093] By solving the state equations and variational equations in parallel, raw prediction data and raw compensation data for all observation epochs can be generated simultaneously during a single integration process. This avoids the increased computational load caused by repeated integrations, significantly improving computational efficiency. Simultaneously, it ensures the accuracy of the state transition data and parameter sensitivity data, providing accurate error propagation relationships for subsequent residual analysis and least squares solutions. Integrating the raw prediction data, state transition data, and parameter sensitivity data for each observation epoch constitutes complete raw ephemeris data, providing comprehensive data support for subsequent orbital state corrections.

[0094] Thus, after numerical integration based on the initial state vector and the preset dynamic model, the raw ephemeris data containing the original prediction data, state transition data, and parameter sensitivity data for each observation epoch can be obtained.

[0095] In the above implementation, by determining a reasonable initial state vector from the current observation data, an accurate starting point for numerical integration is provided. By constructing a complete dynamic model that includes constraints such as Earth's center of mass gravitational perturbation, Earth's non-spherical gravitational perturbation, atmospheric drag perturbation, solar radiation pressure perturbation, third-body gravitational perturbation, Earth's tidal deformation perturbation, and general relativistic effects, all major perturbation forces affecting the orbital motion of low-Earth orbit satellites are comprehensively covered, significantly improving the accuracy and reliability of orbit prediction and laying a solid physical foundation for achieving real-time high-precision orbit determination on satellites.

[0096] In some implementation methods, please refer to the appendix. Figure 3 Residual analysis is performed based on the original ephemeris data and current observation data to obtain the design matrix data and residual vector data, including: S310. Use the current observation data and the original prediction data to perform quality screening in order to determine the effective solution information of the effective observation epochs in the current orbit determination arc.

[0097] It should be noted that the raw observation data may contain outliers (outliers) caused by receiver noise, multipath effects, or communication anomalies. These outliers can severely distort the construction of the design matrix and residual vector, leading to divergence in the least squares solution or a decrease in orbit determination accuracy. Therefore, it is necessary to remove these outliers through quality screening to retain high-quality, valid data at effective observation epochs. These effective observation epochs refer to observation epochs where the data quality is acceptable, i.e., reliable observation sampling points that can be used for subsequent least squares solutions.

[0098] Specifically, an effective observation epoch has effective solution information, which can refer to the set of relevant observation data possessed by the effective observation epoch, including effective error data and effective residual data.

[0099] Among them, effective error data can refer to the original compensation data corresponding to the effective observation epoch, including state transition data and parameter sensitivity data, while effective residual data can refer to the position residual data of that epoch, that is, the difference between the observed position in the current observation data and the estimated position in the original prediction data.

[0100] Specifically, quality screening is performed using current observation data and raw prediction data to determine the effective solution information for valid observation epochs in the current orbit determination arc. This includes: firstly, for any observation epoch in the current orbit determination arc, the difference between the current observation data and the raw prediction data is calculated to obtain position residual data; then, the position residual data is compared with a preset threshold to determine validity, resulting in a valid determination result; subsequently, if the valid determination result indicates that the current observation epoch is a valid observation epoch, the raw compensation data of the valid observation epoch is used as the valid error data, and the position residual data of the valid observation epoch is used as the valid residual data.

[0101] Among them, the location residual data can refer to the vector difference between the observed position in the current observation data and the estimated position in the original prediction data at the same observation epoch. It can be obtained by subtracting the estimated position from the observed position. Its mathematical form is a three-dimensional column vector, which is used to quantify the deviation between the observed value and the model prediction value.

[0102] Specifically, for any observation epoch in the current orbit determination arc segment First, the estimated position of the epoch is obtained from the original ephemeris data in the original prediction data. The estimated position is defined as the theoretical position of the target satellite obtained through state recursion calculation, and the observed position at that epoch is obtained from the current observation data. That is, the actual measured position output by the GNSS receiver.

[0103] Then calculate the estimated location. and observation location The difference is used to obtain the location residual data. This difference directly reflects the degree of deviation between the current dynamic model and the initial state and the actual motion.

[0104] After obtaining the location residual data, the location residual data and the preset threshold can be used to make an effective judgment and obtain an effective judgment result.

[0105] The preset threshold can refer to the numerical limit used to determine whether the positional residual is within an acceptable range. It can be dynamically determined based on the statistical characteristics of the observed data. For example, it can be 3 times the mean error (3σ criterion), or it can be calculated using robust statistical methods such as the MAD method or the IQR method.

[0106] The valid determination result can be a binary result used to identify whether the observation epoch can participate in the subsequent orbit determination calculation. It can include two cases: valid and invalid, which respectively indicate whether the observation epoch can be adopted for subsequent calculation.

[0107] Specifically, the Euclidean norm of the residual data at each observation epoch is first calculated, and then the norm is compared with a preset threshold. If the norm exceeds the preset threshold, the observation epoch is determined to be an invalid observation epoch, which is then removed and an anomaly is recorded in the log. If the norm does not exceed the preset threshold, the observation epoch is determined to be a valid observation epoch.

[0108] For example, suppose the current arc segment has m observation epochs. The dynamic threshold is determined based on the statistical characteristics of prior observation noise or the number of epochs. For example, one or three times the standard error of the residual norm at all locations can be taken. For each epoch... ,calculate ,like > If the epoch is abnormal, it is recorded and skipped; otherwise... ≤ If the epoch is 1, it is determined to be a valid observation epoch.

[0109] Understandably, this embodiment employs the 3σ criterion: first, the standard deviation σ of the residuals at all epochs is calculated, and then a threshold is set. =3σ. Furthermore, for observational data with numerous extreme outliers, the MAD method can be used, with the threshold calculation formula being median(|residual - median(residual)|)×k, where k is typically taken as 1.4826; alternatively, the standard IQR method in statistics can be used, with the threshold calculation formula being Q3 + 1.5×IQR, where IQR is the difference between the upper quartile Q3 and the lower quartile Q1, to enhance resistance to outliers.

[0110] Furthermore, if the valid determination result indicates that the current observation epoch is a valid observation epoch, the original compensation data of the valid observation epoch is used as the valid error data, and the position residual data of the valid observation epoch is used as the valid residual data.

[0111] That is, for each epoch marked as valid Extract the original compensation data corresponding to this epoch from the original ephemeris data, including state transition data. and parameter sensitivity data This will be used as valid error data. Simultaneously, the position residuals calculated for this epoch will be... This serves as valid residual data. This valid data will be used to subsequently construct the observation matrix and residual vector.

[0112] By calculating the position residual and using a dynamic threshold to remove outliers, abnormal observation points caused by multipath effects, receiver noise, signal obstruction, and other factors can be effectively removed, avoiding the contamination of orbit determination accuracy by outliers. At the same time, high-quality observation data is retained for subsequent calculations, significantly improving the reliability of orbit determination results.

[0113] S320. Based on the effective solution information, perform data assembly processing to obtain design matrix data and residual vector data.

[0114] Specifically, data assembly processing is performed based on effective solution information to obtain design matrix data and residual vector data, including: first, constructing an observation matrix based on effective error data for each effective observation epoch; then stacking the observation matrices of each effective observation epoch to obtain design matrix data; and subsequently stacking the effective residual data of each effective observation epoch to obtain residual vector data.

[0115] For example, assuming that after quality screening, the current orbit determination arc has m valid observation epochs, and for the i-th valid epoch, its valid error data includes state transition data. (Size 6×6) and parameter sensitivity data (Size is 6×1).

[0116] Since the observation data only involves the position component of the target satellite (not its velocity), it is necessary to extract the position-related part from the above data.

[0117] Specifically, the first three rows (corresponding positional components) of the state transition data are taken to form a submatrix. (Size is 3×6). Take the first 3 rows of the parameter sensitivity data to form a subvector. (Size is 3×1). Then, these two parts are horizontally concatenated to obtain the observation matrix for the i-th valid epoch. Its size is 3×7, corresponding to the 7 components of the state vector to be estimated (position increment, velocity increment, atmospheric drag coefficient correction).

[0118] Subsequently, the observation matrix of all m valid observation epochs is... Stacked chronologically from top to bottom, forming a total design matrix H with dimensions (3m) × 7. Simultaneously, the effective residual data from all valid observation epochs are... (Each is a three-dimensional column vector) Stacked from top to bottom in the same order to form a total residual vector y with a dimension of (3m)×1. At this point, the assembly of the design matrix data and the residual vector data is complete.

[0119] By constructing the observation matrix for each valid observation epoch and stacking it, the observation information scattered across various observation epochs can be integrated into a unified system of linear equations, providing a standard input data structure for subsequent least squares solutions and ensuring the standardization and efficiency of the solution process.

[0120] In the above implementation, abnormal observation data was removed through quality screening, ensuring the reliability of the input data. By assembling the data to construct a standardized design matrix and residual vector, a precise linear relationship was established between the observation residual and the estimated state correction, laying a solid data foundation for the subsequent scaling least squares solution and effectively improving the accuracy and stability of the orbit determination results.

[0121] In some implementation methods, please refer to the appendix. Figure 4 The design matrix data and residual vector data are jointly solved using a scaling factor to determine the state correction data, including: S410. Adjust the magnitude of the design matrix data using a scaling factor to obtain a scaled design matrix.

[0122] The scaled design matrix can refer to the matrix obtained by normalizing the dimensions of the design matrix data, that is, by scaling each column of the design matrix data according to the typical magnitude of the corresponding parameter.

[0123] It should be noted that a diagonal scaling factor, pre-set based on the typical physical magnitude of the parameter to be estimated, can be used to multiply each column of the design matrix data by the reciprocal of the corresponding scaling factor, so that the column vectors corresponding to position, velocity and atmospheric drag coefficients are normalized to a similar numerical range, thereby eliminating the difference in numerical magnitude between parameters of different dimensions.

[0124] Specifically, the scaling factor S can be defined by the following formula: In the formula, The dimensions of different components in the design matrix data are represented by i, which can be 7. It includes 7 parameters, namely the three-dimensional position components (rx, ry, rz), the three-dimensional velocity components (vx, vy, vz), and the atmospheric drag coefficient (cd).

[0125] For example, the design matrix data is H (3m×7 dimension, where m is the number of effective observation epochs), and the scaling factor is a 7×7 diagonal matrix S, whose diagonal elements are set according to the typical magnitude of the parameters to be estimated: S=diag([1e6,1e6,1e6,1e3,1e3,1e3,1]) The first three elements, 1e6, correspond to the position component (meter-level dimension), the middle three elements, 1e3, correspond to the velocity component (meter / second dimension), and the last element, 1, corresponds to the atmospheric drag coefficient (dimensionless).

[0126] Then, the design matrix after scaling Calculated using the following formula: After the above scaling process, the numerical magnitudes of each column of the design matrix are unified to around 1, and the normal equations... The condition number is reduced from the original 1e12 to the 1e2~1e3 order of magnitude, resulting in a fundamental improvement in numerical stability.

[0127] S420. Based on the scaled design matrix and residual vector data, the state correction data is obtained by solving the least squares method.

[0128] Specifically, based on the scaled design matrix and residual vector data, the state correction data is obtained by solving the least squares method, including: first, based on the scaled design matrix and residual vector data, a system of linear equations is constructed and solved by least squares to obtain the scaled correction data; then, the scaled correction data is restored by using the inverse transformation of the scaling factor to obtain the state correction data.

[0129] Scaling correction data can refer to intermediate state correction quantities after magnitude normalization, that is, dimensionless correction vectors where the magnitude of each component is uniformly around 1.

[0130] Specifically, firstly, a system of linear equations is constructed based on the scaled design matrix and the original residual vector data. This system is then solved using the least squares method to obtain normalized scaled correction data. Next, the inverse matrix of the scaling factor is used to restore the dimensions of the scaled correction data, converting it into state correction data with the original physical dimensions. Subsequently, it is checked whether all elements of the state correction data are finite values, and the norms of the position correction and velocity correction are calculated. Finally, the atmospheric drag coefficient correction is physically limited to ensure it remains within a reasonable range.

[0131] For example, the residual vector data y can first be represented as: In the formula, H represents the design matrix data; This represents the initial correction data.

[0132] Introducing a scaling factor S allows for the correction of the state data to be solved. Defined as: That is, the initial correction data have: Substituting this into the original least squares problem (i.e., y = H· From the approximate relationship, we can obtain the scaled system of linear equations: Based on this, because Since the magnitudes of all columns are uniform, the system of equations can be solved directly using the standard least squares method. That is, the least squares problem can be defined as: Further transforming this into solving a scaled problem, we obtain the following equation: Finally, by iteratively solving the above equation, the final state correction data can be obtained.

[0133] After obtaining the scaled corrected data, the original state corrected data can be obtained by inverse transformation of the scaling factor.

[0134] For example, the restoration can be performed using the following formula: Based on the above formula, the intermediate corrected data after scaling and restoration can be obtained.

[0135] Furthermore, it is possible to further check whether all elements in the intermediate correction data are finite values ​​that are not NaN or Inf, and to calculate the position correction norm dr = ∥Δx[0:3]∥ and the velocity correction norm dv = ∥Δx[3:6]∥, while limiting the atmospheric drag coefficient to the range of [0.5, 10.0]. The final corrected state data is then obtained.

[0136] By using inverse transformation, the normalized intermediate correction values ​​are restored to state correction data with real physical meaning. This not only retains the numerical stability advantage brought by scaling, but also ensures the physical validity of the correction values, which can be directly used to update the satellite's orbital state.

[0137] Understandably, after obtaining the state correction data, it is superimposed on the current state vector to be estimated, and then the state recursion calculation, residual analysis, and scaled least squares solution are re-executed, forming an iterative optimization process. Iteration continues until the convergence condition is met or the maximum number of iterations is reached: the convergence condition is that the position correction norm is less than 1e-3 meters and the velocity correction norm is less than 1e-6 meters / second; the maximum number of iterations is set to 10. If the convergence condition is met, the final state correction data is output; if convergence is not achieved after reaching the maximum number of iterations, the current optimal state correction data is output and a convergence warning is issued.

[0138] In the above implementation, by introducing a scaling factor based on physical priors to adjust the magnitude of the design matrix, the numerical ill-conditioning caused by the difference in dimensions is fundamentally solved without increasing the amount of computation. This significantly improves the stability and efficiency of least squares solution under the limited floating-point precision environment on the satellite, while ensuring the physical validity of the solution results. It also adapts to the limited computing power constraints of the onboard computer and provides algorithmic support for realizing real-time high-precision orbit determination on the satellite.

[0139] It should be understood that although the steps in the flowchart above are shown sequentially as indicated by the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowchart above may include multiple steps or stages, which are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages in other steps.

[0140] This specification also provides an on-board satellite orbit determination device 500, such as... Figure 5 As shown, it includes: a data acquisition module 510, a state recursion module 520, a residual analysis module 530, a joint settlement module 540, and a track determination module 550, wherein: The data acquisition module 510 is used to acquire the current observation data of the target satellite in the current orbit determination arc.

[0141] The state recursion module 520 is used to perform state recursion calculations based on the current observation data to obtain the original ephemeris data; wherein, the original ephemeris data includes the original prediction data and original compensation data of all observation epochs in the current orbit determination arc.

[0142] The residual analysis module 530 is used to perform residual analysis based on the original ephemeris data and the current observation data to obtain the design matrix data and residual vector data.

[0143] The joint settlement module 540 is used to perform joint calculations on the design matrix data and residual vector data using a scaling factor to determine the state correction data; wherein, the scaling factor is a diagonal scaling matrix used to eliminate the difference in numerical magnitude between parameters of different dimensions.

[0144] The orbit determination module 550 is used to update the original prediction data using state correction data to obtain the target orbit determination data.

[0145] In some implementations, the joint settlement module 540 is also used to adjust the magnitude of the design matrix data using a scaling factor to obtain a scaled design matrix; based on the scaled design matrix and the residual vector data, the state correction data is obtained by solving the least squares method.

[0146] In some implementations, the joint settlement module 540 is also used to construct a system of linear equations and perform least squares solution based on the scaled design matrix and residual vector data to obtain scaled corrected data; and to restore the scaled corrected data by using the inverse transformation of the scaling factor to obtain state corrected data.

[0147] In some implementations, the residual analysis module 530 is also used to perform quality screening using the current observation data and the original prediction data to determine the effective solution information of the effective observation epoch in the current orbit determination arc; wherein, the effective solution information includes effective error data and effective residual data; and data assembly processing is performed based on the effective solution information to obtain design matrix data and residual vector data.

[0148] In some implementations, the residual analysis module 530 is further configured to calculate the difference between the current observation data and the original prediction data for any observation epoch in the current orbit determination arc to obtain position residual data; to make a valid determination by comparing the position residual data with a preset threshold to obtain a valid determination result; and if the valid determination result indicates that the current observation epoch is a valid observation epoch, to use the original compensation data of the valid observation epoch as the valid error data and the position residual data of the valid observation epoch as the valid residual data.

[0149] In some implementations, the residual analysis module 530 is further configured to construct an observation matrix based on the effective error data for each effective observation epoch; stack the observation matrices of each effective observation epoch to obtain design matrix data; and stack the effective residual data of each effective observation epoch to obtain residual vector data.

[0150] In some implementations, the state recursion module 520 is also used to determine an initial state vector in the current observation data; and to perform numerical integration based on the initial state vector and a preset dynamic model to generate raw ephemeris data; wherein the preset dynamic model includes constraints on Earth's center of mass gravitational perturbation, Earth's non-spherical gravitational perturbation, atmospheric drag perturbation, solar radiation pressure perturbation, third-body gravitational perturbation, Earth's tidal deformation perturbation, and general relativistic effects.

[0151] In some implementations, the original compensation data includes state transition data and parameter sensitivity data; the state recursion module 520 is also used to integrate from the starting point to the ending point of the current orbital segment based on the initial state vector and the preset dynamic model, and solve the state equation and variational equation in parallel to obtain the original prediction data, state transition data and parameter sensitivity data for each observation epoch.

[0152] For specific limitations regarding an on-board satellite orbit determination device, please refer to the limitations regarding an on-board satellite orbit determination method described above, which will not be repeated here. Each module in the aforementioned on-board satellite orbit determination device can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device in hardware form, or stored in the memory of a computer device in software form, so that the processor can call and execute the operations corresponding to each module.

[0153] In this embodiment, a satellite on-board orbit determination device is presented in the form of a functional unit. Here, a unit refers to an ASIC (Application Specific Integrated Circuit) circuit, a processor and memory that execute one or more software or fixed programs, and / or other devices that can provide the above functions.

[0154] Please see Figure 6 , Figure 6 This is a schematic diagram of the structure of a computer device provided in an embodiment of this application, such as... Figure 6 As shown, the computer device includes one or more processors 10, memory 20, and interfaces for connecting the components, including high-speed interfaces and low-speed interfaces. The components communicate with each other via different buses and can be mounted on a common motherboard or otherwise installed as needed. The processor can process instructions executed within the computer device, including instructions stored in or on memory to display graphical information of a GUI on an external input / output device (such as a display device coupled to the interface). In some alternative implementations, multiple processors and / or multiple buses can be used with multiple memories and multiple memory modules, if desired. Similarly, multiple computer devices can be connected, each providing some of the necessary operations. Figure 6 Take a processor 10 as an example.

[0155] Processor 10 may be a central processing unit, a network processor, or a combination thereof. Processor 10 may further include a hardware chip. The hardware chip may be an application-specific integrated circuit (ASIC), a programmable logic device (PLD), or a combination thereof. The programmable logic device may be a complex programmable logic device (CAMP), a field-programmable gate array (FPGA), a general-purpose array logic (GDA), or any combination thereof.

[0156] The memory 20 stores instructions executable by at least one processor 10 to cause the at least one processor 10 to perform the method shown in the above embodiments.

[0157] The memory 20 may include a program storage area and a data storage area. The program storage area may store the operating system and applications required for at least one function; the data storage area may store data created based on the use of the computer device. Furthermore, the memory 20 may include high-speed random access memory and may also include non-transitory memory, such as at least one disk storage device, flash memory device, or other non-transitory solid-state storage device. In some alternative embodiments, the memory 20 may optionally include memory remotely located relative to the processor 10, and these remote memories may be connected to the computer device via a network. Examples of such networks include, but are not limited to, the Internet, intranets, local area networks, mobile communication networks, and combinations thereof.

[0158] The memory 20 may include volatile memory, such as random access memory; the memory may also include non-volatile memory, such as flash memory, hard disk or solid-state drive; the memory 20 may also include a combination of the above types of memory.

[0159] The computer device also includes an input device 30 and an output device 40. The processor 10, memory 20, input device 30, and output device 40 can be connected via a bus or other means. Figure 6 Taking the example of a connection between China and Israel via a bus.

[0160] Input device 30 can receive input numerical or character information, and generate key signal inputs related to user settings and function control of the computer device, such as a touchscreen, keypad, mouse, trackpad, touchpad, joystick, one or more mouse buttons, trackball, joystick, etc. Output device 40 may include display devices, auxiliary lighting devices (e.g., LEDs), and haptic feedback devices (e.g., vibration motors). The aforementioned display devices include, but are not limited to, liquid crystal displays, light-emitting diodes, displays, and plasma displays. In some alternative embodiments, the display device may be a touchscreen.

[0161] Although embodiments of this application have been described in conjunction with the accompanying drawings, those skilled in the art can make various modifications and variations without departing from the spirit and scope of this application, and all such modifications and variations fall within the scope defined by the appended claims.

[0162] The apparatus, module, or unit described in the above embodiments can be implemented by a computer chip or entity, or by a product having a certain function. A typical implementation device is a computer. Specifically, the computer can be, for example, a personal computer, laptop computer, cellular phone, camera phone, smartphone, personal digital assistant, media player, navigation device, email device, game console, tablet computer, wearable device, or any combination of these devices.

[0163] For ease of description, the above devices are described separately by function as various units. Of course, in implementing this application, the functions of each unit can be implemented in one or more software and / or hardware.

[0164] Those skilled in the art will understand that embodiments of this application can be provided as methods, apparatus, or computer devices. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application 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.

[0165] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus, and computer devices according to embodiments of this application. 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... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0166] 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.

[0167] 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.

[0168] It should also be noted that the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0169] The various embodiments in this specification are described in a progressive manner. Similar or identical parts between embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the apparatus embodiments are basically similar to the method embodiments, so the description is relatively simple; relevant parts can be referred to the descriptions of the method embodiments.

[0170] The above description is merely an embodiment of this application and is not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

[0171] Although embodiments of this application have been described in conjunction with the accompanying drawings, those skilled in the art can make various modifications and variations without departing from the spirit and scope of this application, and such modifications and variations all fall within the scope defined by the appended claims.

Claims

1. A satellite on-board orbit determination method, characterized in that, The method includes: Acquire current observation data of the target satellite in the current orbit determination arc; Based on the current observation data, a state recursive calculation is performed to obtain the original ephemeris data; wherein, the original ephemeris data includes the original prediction data and original compensation data of all observation epochs in the current orbit determination arc; Residual analysis is performed based on the original ephemeris data and the current observation data to obtain design matrix data and residual vector data; The design matrix data is adjusted by a scaling factor to obtain a scaled design matrix; wherein the scaling factor is a diagonal scaling matrix used to eliminate the difference in numerical magnitude between parameters of different dimensions. Based on the scaled design matrix and the residual vector data, a system of linear equations is constructed and solved by least squares to obtain the scaled corrected data. The scaled and corrected data is restored by using the inverse transformation of the scaling factor to obtain the state corrected data; The original predicted data is updated using the state correction data to obtain the target orbit determination data.

2. The method according to claim 1, characterized in that, The residual analysis based on the original ephemeris data and the current observation data to obtain design matrix data and residual vector data includes: The current observation data and the original prediction data are used for quality screening to determine the effective solution information of the effective observation epochs in the current orbit determination arc; wherein, the effective solution information includes effective error data and effective residual data; Data assembly processing is performed based on the effective solution information to obtain the design matrix data and the residual vector data.

3. The method according to claim 2, characterized in that, The process of using the current observation data and the original prediction data for quality screening to determine the effective solution information for the effective observation epochs in the current orbit determination arc includes: For any observation epoch in the current orbit determination arc, the difference between the current observation data and the original prediction data is calculated to obtain the position residual data; The positional residual data is used to make an effective determination with a preset threshold to obtain an effective determination result; If the valid determination result indicates that the current observation epoch is the valid observation epoch, the original compensation data of the valid observation epoch is used as the valid error data, and the position residual data of the valid observation epoch is used as the valid residual data.

4. The method according to claim 3, characterized in that, The data assembly process based on the effective solution information to obtain the design matrix data and the residual vector data includes: For each of the effective observation epochs, an observation matrix is ​​constructed based on the effective error data; The observation matrices of each effective observation epoch are stacked to obtain the design matrix data; The effective residual data of each effective observation epoch are stacked to obtain the residual vector data.

5. The method according to claim 1, characterized in that, The process of performing state recursion calculations based on the current observation data to obtain the original ephemeris data includes: Determine the initial state vector from the current observation data; Numerical integration is performed based on the initial state vector and the preset dynamic model to generate raw ephemeris data; wherein, the preset dynamic model includes constraints on gravitational perturbation of Earth's center of mass, gravitational perturbation of Earth's non-spherical shape, atmospheric drag perturbation, solar radiation pressure perturbation, third-body gravitational perturbation, Earth's tidal deformation perturbation, and general relativistic effects.

6. The method according to claim 5, characterized in that, The original compensation data includes state transition data and parameter sensitivity data; The numerical integration process is performed in the following manner: Based on the initial state vector and the preset dynamic model, the state equation and variational equation are solved in parallel from the starting point to the ending point of the current orbital segment to obtain the original prediction data, state transition data and parameter sensitivity data for each observation epoch.

7. A satellite on-board orbit determination device, characterized in that, The device includes: The data acquisition module is used to acquire the current observation data of the target satellite in the current orbit determination arc. The state recursion module is used to perform state recursion calculations based on the current observation data to obtain the original ephemeris data; wherein, the original ephemeris data includes the original prediction data and original compensation data of all observation epochs in the current orbit determination arc; The residual analysis module is used to perform residual analysis based on the original ephemeris data and the current observation data to obtain design matrix data and residual vector data; The joint settlement module is used to adjust the magnitude of the design matrix data using a scaling factor to obtain a scaled design matrix; based on the scaled design matrix and the residual vector data, a system of linear equations is constructed and solved by least squares to obtain scaled corrected data; the scaled corrected data is restored using the inverse transformation of the scaling factor to obtain state corrected data; wherein, the scaling factor is a diagonal scaling matrix used to eliminate the difference in numerical magnitude between parameters of different dimensions. The orbit determination module is used to update the original prediction data using the state correction data to obtain the target orbit determination data.

8. A computer device, characterized in that, include: A memory and a processor, the memory and the processor being communicatively connected to each other, the memory storing computer instructions, the processor executing the computer instructions to perform the method of any one of claims 1 to 6.

Images