Satellite cluster precise orbit determination method based on multi-source data fusion

Abstract

The application provides a satellite cluster precise orbit determination method based on multi-source data fusion, and belongs to the technical field of satellites. Basic parameters are set, a dynamic model and a state transition equation of a satellite are established, absolute positions and velocities of the satellite cluster are taken as observation quantities, an observation equation of the absolute state of the satellite is established, inter-satellite relative positions and velocities of the satellite cluster are taken as observation quantities, an observation model of the relative state of the satellite is established, a state quantity optimal estimation value is obtained by solving a cluster positioning state updating equation by using a Newton method, and the spatial position and velocity of the satellite cluster are calibrated. The method fuses the absolute state observation and the inter-satellite relative state observation of the satellite cluster, overcomes the limitation of a single measurement method in measurement accuracy, has expansibility, can fuse various observation methods, and further improves the orbit determination accuracy of the satellite cluster.

Description

Technical Field

[0001] This invention relates to the field of satellite technology, and in particular to a method for precise orbit determination of satellite constellations based on multi-source data fusion. Background Technology

[0002] A satellite constellation refers to a network of satellites forming a specific configuration, orbiting the Earth, maintaining a certain distance and orientation between them, communicating and coordinating with each other to form a unified whole and collaboratively accomplish a specific mission. Precise autonomous space positioning is essential to maintain this specific configuration and spatial position for effective mission completion.

[0003] The basic principle of autonomous orbit determination is to determine its own orbit without relying on a ground-based observation network, using onboard GNSS receivers, optical camera observations, radar ranging, laser ranging, and inter-satellite information exchange. Radar ranging and laser ranging require satellites to be equipped with specific instruments and communicate with the ground, which increases satellite mass, wastes satellite installation space, and increases energy consumption. Using a single GNSS positioning or inter-satellite optical observation positioning method is affected by satellite spatial position, relative distance, and Earth's obstruction, resulting in insufficient positioning accuracy for some satellites in a constellation at certain spatial locations. Summary of the Invention

[0004] The technical problem solved by this invention is that the increased burden on satellites and the lack of positioning accuracy caused by a single positioning method in the prior art.

[0005] To address the aforementioned problems, this invention provides a multi-source fusion positioning method for satellite constellations that utilizes GNSS positioning and inter-satellite optical observation, employing a Bayesian network model. In this scheme, a GNSS receiver acquires measurement data from GNSS satellites to measure the absolute position of individual satellites; inter-satellite optical observation equipment is used to measure the relative positions between satellites in the constellation. This scheme utilizes onboard monitoring cameras and other optical cameras, eliminating the need for additional dedicated cameras or ranging equipment, thus saving payload space and reducing satellite mass and energy consumption. By establishing a Bayesian network-based fusion positioning model, the optimal estimate of the satellite state is calculated.

[0006] The equipment required for the method described in this invention includes at least inter-satellite relative position observation equipment and satellite absolute position positioning equipment. The equipment used only needs to meet the functional requirements of relative and absolute positioning; it is not limited to a specific type of equipment and needs to be determined based on the satellite's payload space and functions in the actual application. In this invention, a typical optical camera is used as the equipment for inter-satellite relative position observation, and a GNSS receiver is used as the equipment for measuring the satellite's absolute position. For satellites with more positioning observation equipment, this method can be extended in the same form to make full use of more observation methods.

[0007] Satellites in a constellation observe neighboring satellites using optical cameras. By measuring the azimuth of the observed satellite within its field of view and combining this with the attitude information of the observing satellite, two azimuth angles relative to the measuring satellite can be obtained. Some optical observation methods, such as using binocular cameras or placing dimensional markers on the satellite, can obtain relative distances between satellites. For optical observation equipment lacking this capability, distance measuring devices such as laser rangefinders can be installed. For satellites equipped with GNSS receivers, the absolute azimuth of the satellite in inertial space can be determined by receiving positioning signals from GNSS satellites. Both observation methods inevitably introduce observation noise, affecting positioning accuracy. Therefore, an algorithm for obtaining the optimal estimation of the positioning state is needed to fuse the two or more observation methods and reduce the impact of random noise on positioning accuracy.

[0008] It should be noted that the algorithm in this method can be implemented in the onboard computer, or it can be implemented in the ground computer by sending the observation information to the ground station, and then sent to the satellite constellation. Where the computing power and energy consumption of the onboard computer allow, the former is preferred to reduce time delays in information transmission and eliminate unnecessary satellite payloads.

[0009] To address the above problems, the present invention provides the following solution.

[0010] A method for precise orbit determination of satellite constellations based on multi-source data fusion includes:

[0011] Step 1: Obtain basic parameters, acquire satellite cluster orbital state and attitude information, acquire state quantity estimates and observations, and calculate the covariance matrix of state prediction noise, the covariance matrix of GNSS measurement noise, and the covariance matrix of inter-satellite optical observation noise.

[0012] Step 2: Obtain the historical state data of the satellites. Based on the configuration and orbital parameters of the satellite constellation, establish the state transition equation of the satellite constellation. Based on the covariance matrix of the state prediction noise obtained in Step 1, establish the joint state transition equation.

[0013] Step 3: Using the absolute position and absolute velocity of the satellite cluster as observations, establish the observation equation for the absolute state of the satellites. Based on the covariance matrix of the GNSS measurement noise obtained in Step 1, establish the joint absolute observation equation and calculate the state update quantity of the absolute position observation of a single satellite using the GNSS method.

[0014] Step 4: Using the spatial relative position and relative velocity of the satellite cluster as observations, establish the observation equation for the relative state of the satellites. Based on the covariance matrix of the inter-satellite optical observation noise obtained in Step 1, establish the joint relative observation equation and calculate the state update amount of the inter-satellite relative position observed by the optical camera.

[0015] Step 5: Based on the state transition of the satellite cluster and the Bayesian network model of multi-source fusion observation, add the joint state transition equation obtained in Step 2, the joint absolute observation equation obtained in Step 3, and the joint relative observation equation obtained in Step 4 to obtain the cluster positioning state update equation. Use Newton's method to solve the cluster positioning state update equation to obtain the optimal estimate of the satellite state variables at the current time. Use the optimal estimate as the state update variable to update the satellite cluster state, realize the spatial position and velocity calibration of the satellite cluster, and thus determine the precise orbit of the satellite cluster.

[0016] As one aspect of the present invention, in step one, the basic parameters include the gravitational constant of the central celestial body and the configuration parameters of the satellite cluster. The configuration parameters of the satellite cluster are the orbital elements of each satellite in the cluster, and the central celestial body is the Earth.

[0017] The orbital status and attitude information of the satellite constellation includes: the orbital status of the reference satellite in the constellation, and the attitude information of each satellite in the constellation;

[0018] The method for obtaining the state quantity estimates and observations is as follows: Obtain the state quantities within a certain period preceding the position calibration solution time, as the state estimates. The length of this preceding period is determined based on the specific scenario's accuracy requirements. The state estimate for satellite i at time k is... in These represent the first, second, and third position components of the satellite in the geocentric inertial coordinate system. The first, second, and third velocity components of the satellite in the geocentric inertial coordinate system are used; the observation records stored in the satellite's onboard computer are used as the observation values.

[0019] The methods for calculating the covariance matrix of state prediction noise, the covariance matrix of GNSS measurement noise, and the covariance matrix of inter-satellite optical observation noise are as follows:

[0020] Based on the satellite orbital state data stored in the satellite's onboard computer, the state prediction noise vector W is obtained. kThen calculate the covariance matrix of the state prediction noise:

[0021]

[0022] Where k is the k-th time, Q k Let W be the covariance matrix of the state prediction noise. k For state prediction noise, W k (i) is W k The i-th element in the vector, dim, represents the state variable X. k Number of elements in the middle For mathematical expectation, E((W) k (1)) 2 To calculate the variance of the first element in the state prediction noise, E((W) k (2)) 2 To calculate the variance of the second element in the state prediction noise, E((W) k (dim)) 2 To calculate the variance of the last element in the state prediction noise, E(W) k (1)·W k (2)) and E(W k (2)·W k (1) both calculate the covariance of the first and second elements in the state prediction noise, E(W k (1)·W k (dim)) and E(W k (dim)·W k (1) both refer to calculating the covariance of the first and last elements in the state prediction noise, E(W k (2)·W k (dim)) and E(W k (dim)·W k (2) Both are to find the covariance of the second and last elements in the state prediction noise. The covariance matrix of the state prediction noise is a symmetric matrix.

[0023] Based on GNSS observation records, the GNSS observation noise was obtained. Vector, calculation The variance of each element s in Where s is the navigation satellite number, the covariance matrix of GNSS measurement noise is obtained:

[0024]

[0025] Where i is the satellite number, k represents the k-th time, and N k This represents the number of navigation satellites that can receive signals at the current moment. Represents the covariance matrix of GNSS measurement noise. This represents the variance of the first element in the GNSS observation noise vector. This represents the variance of the second element in the GNSS observation noise vector. This represents the variance of the last element in the GNSS observation noise vector. The covariance matrix of the GNSS measurement noise is a diagonal matrix. Repeat the above process to calculate the covariance matrix of the state prediction noise for each satellite.

[0026] Based on inter-satellite optical observation records, the observation noise of the direct relative position between satellite i and satellite j was calculated. Calculate the covariance matrix of inter-satellite optical observation noise:

[0027]

[0028] Where i and j are the numbers of the two satellites that are being observed relative to each other, and k is the time at which the observation occurs. Let E(·) be the covariance matrix of inter-satellite optical observation noise, and E(·) be the expected value. To find the variance of the first element in the relative position observation noise, To find the variance of the second element in the relative position observation noise, To find the variance of the third element in the relative position observation noise, and Both are used to calculate the covariance of the first and second elements in the relative position observation noise. and Both are used to calculate the covariance of the first and third elements in the relative position observation noise. and The covariance of the second and third elements in the relative position observation noise is calculated. The above process is repeated to calculate the covariance matrix of the relative inter-satellite optical observation noise for each pair of satellites (i,j) with relative observation relationship.

[0029] As one aspect of the present invention, the historical records of satellite state variables are obtained, and a state transition equation for the satellite constellation is established based on the configuration and orbital parameters of the satellite constellation. A joint state transition equation is established based on the covariance matrix of the state prediction noise obtained in step one, including:

[0030] S021: Obtain the historical records of satellite state variables and establish the state transition equations for the satellite constellation;

[0031] The orbital dynamics equations of a satellite under the influence of perturbation are as follows:

[0032]

[0033] Where a is the semi-major axis of the orbit, e is the eccentricity, Ω is the right ascension of the ascending node, i is the orbital inclination, ω is the argument of perigee, M is the mean perigee, p is the semi-major radius of the orbit, μ is the gravitational constant of Earth, θ is the true perigee, t is time, r is the distance from the Earth's center, and f r f u and f h These are the x, y, and z components of the perturbation force in the satellite's first orbital coordinate system, respectively; the parameters of the orbital dynamics equations are the configuration and orbital parameters of the satellite constellation.

[0034] The state transition equations for the satellite constellation are:

[0035]

[0036] in, C represents the number of satellites in the cluster, X k Let Xk-1 be the state variable of all satellites in the cluster at time k, and XtoE be the state variable of all satellites in the cluster at time k-1. k-1 f is the transformation function from Cartesian coordinates to Kepler orbital features. k Let f be the state transition equation at time k. e,k The state transition equation in terms of Kepler orbital elements, W k The state prediction noise is Gaussian white noise with zero mean, E k E represents the Kepler orbital features of all satellites in the cluster at time k. k-1 This represents the Kepler orbital features of all satellites in the cluster at time k-1. The rate of change of the satellite's orbital elements over time;

[0037] Linearizing the state transition equations yields the state transition equations for the satellite constellation:

[0038] X k =F k (X k-1 )+W k ,

[0039] Among them, X k Let X be the state variable of all satellites in the cluster at time k. k-1 F represents the state variables of all satellites in the cluster at time k-1. k Let W be the state transition matrix of the satellite constellation at time k. k The noise vector is used to predict the state.

[0040] S022: Establish the joint state transition equations:

[0041] Let the state update vector be... Where k represents the k-th time, Ws Indicates a time window. The initial linearization point for iterative computation. This represents the updated state quantity estimate, calculated by combining the joint prediction matrix and the state expectation vector based on the satellite's state transition matrix and the covariance matrix of the state prediction noise at each time point within the time window.

[0042] Get the entire time window W s The joint state transition equation corresponding to the cluster fusion positioning is:

[0043]

[0044] Here, the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where t represents time, k represents the current time of the solution, and W represents the state update vector. s Indicates a time window. This is the joint prediction matrix for each time t within the time window. Let be the state expectation vector at each time t within the time window.

[0045] As one aspect of the present invention, using the absolute spatial position and absolute velocity of the satellite constellation as observations, an observation equation for the absolute state of the satellites is established. A joint absolute observation equation is established based on the covariance matrix of the GNSS measurement noise obtained in step one. The state update quantities for the absolute position observations of a single satellite using the GNSS method are calculated, including:

[0046] S031: Obtain the absolute spatial position of the satellite constellation via the onboard GNSS receiver, including pseudorange observations received by the satellites from the navigation satellites. Record the pseudorange observations received by the satellites as... At time k, satellite i receives the pseudorange observation value from navigation satellite number 1. At time k, satellite i receives pseudorange observations from navigation satellite number 2. At time k, satellite i receives a signal numbered N. k pseudorange observations of navigation satellites, N k Let be the number of navigation satellites received by satellite i at time k;

[0047] The observation equations for establishing the absolute state of satellites for GNSS positioning are as follows:

[0048]

[0049] Where i is the satellite number, k is the k-th time, and Xk For satellite state variables, h i Let V be the observation function. i k For observing noise, the observation function for satellite i is:

[0050]

[0051]

[0052] Where i is the satellite number, k is the k-th time, and s is the navigation satellite number. Let h be the state variable of satellite i at time k. i Let h be the observation function of satellite i. GNSS(i1) Let N be the observation function of satellite i relative to navigation satellite 1, and so on. k Let k be the number of navigation satellites that the satellite can receive at time k. Let i be the position of satellite i at time k. Let be the position of navigation satellite s at time k;

[0053] S032: Establish the joint absolute observation equations and calculate the state update quantities of single-satellite absolute position observations using the GNSS method:

[0054] Let the state update vector be... Where k represents the k-th time. The initial linearization point for iterative computation. To obtain the updated state quantity estimates, the joint absolute observation matrix and the absolute observation expectation vector are calculated based on the GNSS observation equations and the covariance matrix of the GNSS measurement noise at each time point within the time window.

[0055] Based on the above results, the entire time window W is obtained. s The following joint absolute observation equation corresponds to cluster fusion localization:

[0056]

[0057] Here, the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where i represents the satellite number, s represents the navigation satellite number, t represents the time, k represents the current time of the solution, and W... s Indicates a time window. This represents the joint absolute observation matrix of satellite i relative to navigation satellite s at each time t within the time window. Let t be the absolute observation expectation vector of satellite i relative to navigation satellite s at each time point t within the time window;

[0058] Based on the entire time window W s The joint absolute observation equation corresponding to cluster fusion positioning is used to calculate the state update of single-satellite absolute position observations using the GNSS method.

[0059] As one aspect of the present invention, using the spatial relative position and relative velocity of the satellite constellation as observations, an observation equation for the relative state of the satellites is established. Based on the covariance matrix of the inter-satellite optical observation noise obtained in step one, the relative observation equation is jointly calculated to determine the state update quantity of the inter-satellite relative position observed by the optical camera, including:

[0060] S041: Observe other satellites in the cluster using the onboard optical camera to obtain relative observation information, including relative distance L. ij and two angles measured only α ij β ij This refers to the spatial relative position and relative velocity of the satellite constellation. The spatial relative position and relative velocity of the satellite constellation are used as the observations, and these observations are... The observation equation for obtaining the relative state of the satellite is:

[0061]

[0062]

[0063] Where i and j are the satellite numbers of two satellites with a relative observation relationship, k is the time k, and X k Let k be the state variable of the satellite at time k. These are the coordinate components of satellite j in the first orbital coordinate system of satellite i at time k;

[0064] Linearize the observation equations and express them in matrix form:

[0065]

[0066] Where i and j are the satellite numbers of two satellites with a relative observation relationship, and k is the time at time k. V is the relative observation matrix. k ij For relative observation noise, X k For the state variables of the satellite constellation;

[0067] S042: Establish the joint relative observation equation and calculate the state update of the inter-satellite relative positions observed by the optical camera;

[0068] Let the state update vector be... Where k represents the k-th time. The initial linearization point for iterative computation. This represents the updated state quantity estimate, calculated based on the observation equations of the satellite's relative state at each moment within the time window and the covariance matrix of the inter-satellite optical observation noise, along with the joint relative observation matrix and the relative observation expectation vector.

[0069] Get the entire time window W s The following joint relative observation equation corresponds to cluster fusion positioning:

[0070]

[0071] Here, the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where i and j represent satellite numbers, t represents time, k represents the current solution time, and W represents the current solution time. s Indicates a time window. This represents the joint relative observation matrix of satellite i relative to satellite j at each time t within the time window. Let t be the relative observation expectation vector of satellite i relative to satellite j at each time t within the time window.

[0072] As one aspect of the present invention, based on the Bayesian network model of satellite cluster state transition and multi-source fusion observation, the joint state transition equation obtained in step two, the joint absolute observation equation obtained in step three, and the joint relative observation equation obtained in step four are added to obtain the cluster positioning state update equation. The Newton-Raphson method is used to solve the cluster positioning state update equation to obtain the optimal estimate of the satellite state variables at the current moment. The optimal estimate is used as the state update variable to update the satellite cluster state, thereby achieving spatial position and velocity calibration of the satellite cluster and determining the precise orbit of the satellite cluster, including:

[0073] S051: Add the joint state transition equation, the joint absolute observation equation, and the joint relative observation equation to obtain the cluster positioning state update equation for all times within the time window:

[0074]

[0075] Where, the symbol minarg is the parameter used to find the minimum value of the expression on the right. Let Σ be the norm 2, and let Σ be the summation. Update the state vector. Let W be the state update vector that minimizes the expression, where i and j are satellite numbers, s is the navigation satellite number, t is the time, k is the current time of the solution, and W is the time of the solution.s For time window, This is the joint prediction matrix for each time t within the time window. Let be the state expectation vector at each time t within the time window. This represents the joint absolute observation matrix of satellite i relative to navigation satellite s at each time t within the time window. Let be the absolute observation expectation vector of satellite i relative to navigation satellite s at each time t within the time window. This represents the joint relative observation matrix of satellite i relative to satellite j at each time t within the time window. Let t be the relative observation expectation vector of satellite i relative to satellite j at each time t within the time window;

[0076] S052: The method for iteratively solving the cluster positioning state update equation using Newton's method to obtain the optimal estimate of the satellite state variables at the current moment is as follows: Determine the initial state estimate based on the original observation values, use Newton's optimization algorithm to solve the cluster positioning state update equation, calculate the update vector for the current step, add the update vector to the initial state estimate, update the state estimate, update the cluster positioning state update equation, and repeat the above process; when the difference between the state variable estimates before and after iteration is less than the accuracy requirement required in the specific scenario, or when the number of iterations exceeds the set maximum number of iterations, stop the iteration, and use the final result as the optimal estimate of the satellite state variables at the current moment. Use the optimal estimate as the state update value to update the satellite cluster state, realize the spatial position and velocity calibration of the satellite cluster, and thus determine the precise orbit of the satellite cluster.

[0077] S053: Let k = k + 1, return to step S051, and solve for the optimal estimate of the satellite state variables at the next moment.

[0078] The beneficial effects of this invention are as follows:

[0079] This invention integrates the observation of the absolute position and velocity of a satellite constellation with the observation of the relative position and velocity between satellites, thus combining multiple sources of observation. This overcomes the limitations of a single measurement method in terms of measurement accuracy. At the same time, the method is scalable and can integrate more possible observation methods to further improve the accuracy of spatial position and velocity calibration of satellite constellations. Attached Figure Description

[0080] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0081] The present invention will now be described in further detail with reference to the accompanying drawings.

[0082] The equipment required for the method described in this invention includes at least inter-satellite relative position observation equipment and satellite absolute position positioning equipment. The equipment used only needs to meet the functional requirements of relative and absolute positioning; it is not limited to a specific type of equipment and needs to be determined based on the satellite's payload space and functions in the actual application. In this invention, a typical optical camera is used as the equipment for inter-satellite relative position observation, and a GNSS receiver is used as the equipment for measuring the satellite's absolute position. For satellites with more positioning observation equipment, this method can be extended in the same form to make full use of more observation methods.

[0083] Satellites in a constellation observe neighboring satellites using optical cameras. By measuring the azimuth of the observed satellite within its field of view and combining this with the attitude information of the observing satellite, two azimuth angles relative to the measuring satellite can be obtained. Some optical observation methods, such as using binocular cameras or placing dimensional markers on the satellite, can obtain relative distances between satellites. For optical observation equipment lacking this capability, distance measuring devices such as laser rangefinders can be installed. For satellites equipped with GNSS receivers, the absolute azimuth of the satellite in inertial space can be determined by receiving positioning signals from GNSS satellites. Both observation methods inevitably introduce observation noise, affecting positioning accuracy. Therefore, an algorithm for obtaining the optimal estimation of the positioning state is needed to fuse the two or more observation methods and reduce the impact of random noise on positioning accuracy.

[0084] It should be noted that the algorithm in this method can be implemented in the onboard computer, or it can be implemented in the ground computer by sending the observation information to the ground station, and then sent to the satellite constellation. Where the computing power and energy consumption of the onboard computer allow, the former is preferred to reduce time delays in information transmission and eliminate unnecessary satellite payloads.

[0085] This embodiment describes a method for precise orbit determination of satellite constellations based on multi-source data fusion, such as... Figure 1 As shown, it includes:

[0086] Step 1: Obtain basic parameters, acquire satellite cluster orbital state and attitude information, acquire state quantity estimates and observations, and calculate the covariance matrix of state prediction noise, the covariance matrix of GNSS measurement noise, and the covariance matrix of inter-satellite optical observation noise.

[0087] Understandably, in step one, the basic parameters include the gravitational constant of the central celestial body and the configuration parameters of the satellite cluster. The configuration parameters of the satellite cluster are the orbital elements of each satellite in the cluster, and the central celestial body is the Earth.

[0088] The orbital and attitude information of the satellite constellation includes: the orbital status of the reference satellite in the constellation, and the attitude information of each satellite in the constellation; a geocentric inertial coordinate system is constructed, with the origin at the Earth's center, the xoy plane coinciding with the equatorial plane, and the x-axis pointing to the vernal equinox, usually using the direction of the vernal equinox at 12:00 on January 1, 2000 as the standard, hence also called the J2000 coordinate system; the z-axis is perpendicular to the equatorial plane and points to the North Pole; the y-axis direction is determined by the x-axis and z-axis using the right-hand rule; this coordinate system does not rotate with the Earth's rotation and has a fixed direction in space, and is considered an inertial coordinate system; the absolute position measurement and motion state description of the satellites are usually performed in this coordinate system;

[0089] A first orbital coordinate system is constructed, with the origin at the satellite's center of mass. The x-axis points from the Earth's center of mass to the satellite's center of mass, the z-axis is in the same direction as the orbit's angular momentum and perpendicular to the orbital plane, and the y-axis is determined by the right-hand rule through the x-axis and z-axis. This coordinate system rotates with the satellite's motion and is a non-inertial coordinate system. The relative position observations obtained from inter-satellite optical observations will be described in this coordinate system.

[0090] The method for obtaining the state quantity estimates and observations is as follows: Obtain the state quantities within a certain period preceding the position calibration solution time, as the state estimates. The length of this preceding period is determined based on the specific scenario's accuracy requirements. The state estimate for satellite i at time k is... in These represent the first, second, and third position components of the satellite in the geocentric inertial coordinate system. The first, second, and third velocity components of the satellite in the geocentric inertial coordinate system are used; the observation records stored in the satellite's onboard computer are used as the observation values.

[0091] Understandably, the longer the length and the higher the precision of the previous period, the greater the amount of computation.

[0092] The methods for calculating the covariance matrix of state prediction noise, the covariance matrix of GNSS measurement noise, and the covariance matrix of inter-satellite optical observation noise are as follows:

[0093] Based on the satellite orbital state data stored in the satellite's onboard computer, the state prediction noise vector W is obtained. k Then calculate the covariance matrix of the state prediction noise:

[0094]

[0095] Where k is the k-th time, Q k Let W be the covariance matrix of the state prediction noise. k For state prediction noise, W k (i) is W kThe i-th element in the vector, dim, represents the state variable X. k Number of elements in the middle For mathematical expectation, E((W) k (1)) 2 To calculate the variance of the first element in the state prediction noise, E((W) k (2)) 2 To calculate the variance of the second element in the state prediction noise, E((W) k (dim)) 2 To calculate the variance of the last element in the state prediction noise, E(W) k (1)·W k (2)) and E(W k (2)·W k (1) both calculate the covariance of the first and second elements in the state prediction noise, E(W k (1)·W k (dim)) and E(w) k (dim)·W k (1) both refer to calculating the covariance of the first and last elements in the state prediction noise, E(W k (2)·W k (dim)) and E(W k (dim)·W k (2) Both are to find the covariance of the second and last elements in the state prediction noise. The covariance matrix of the state prediction noise is a symmetric matrix.

[0096] Based on GNSS observation records, the GNSS observation noise was obtained. Vector, calculation The variance of each element s in Where s is the navigation satellite number, the covariance matrix of GNSS measurement noise is obtained:

[0097]

[0098] Where i is the satellite number, k represents the k-th time, and N k This represents the number of navigation satellites that can receive signals at the current moment. Represents the covariance matrix of GNSS measurement noise. This represents the variance of the first element in the GNSS observation noise vector. This represents the variance of the second element in the GNSS observation noise vector. This represents the variance of the last element in the GNSS observation noise vector. The covariance matrix of the GNSS measurement noise is a diagonal matrix. Repeat the above process to calculate the covariance matrix of the state prediction noise for each satellite.

[0099] Based on inter-satellite optical observation records, the observation noise of the direct relative position between satellite i and satellite j was calculated. Calculate the covariance matrix of inter-satellite optical observation noise:

[0100]

[0101] Where i and j are the numbers of the two satellites that are being observed relative to each other, and k is the time at which the observation occurs. Let E(·) be the covariance matrix of inter-satellite optical observation noise, and E(·) be the expected value. To find the variance of the first element in the relative position observation noise, To find the variance of the second element in the relative position observation noise, To find the variance of the third element in the relative position observation noise, and Both are used to calculate the covariance of the first and second elements in the relative position observation noise.

[0102] and Both are used to calculate the covariance of the first and third elements in the relative position observation noise. and The covariance of the second and third elements in the relative position observation noise is calculated. The above process is repeated to calculate the covariance matrix of the relative inter-satellite optical observation noise for each pair of satellites (i,j) with relative observation relationship.

[0103] Step two: Obtain the historical state variables of the satellites. Based on the configuration and orbital parameters of the satellite constellation, establish the state transition equations for the satellite constellation. Based on the covariance matrix of the state prediction noise obtained in step one, establish the joint state transition equations, including:

[0104] S021: Obtain the historical records of satellite state variables and establish the state transition equations for the satellite constellation;

[0105] The orbital dynamics equations of a satellite under the influence of perturbation are as follows:

[0106]

[0107] Where a is the semi-major axis of the orbit, e is the eccentricity, Ω is the right ascension of the ascending node, i is the orbital inclination, ω is the argument of perigee, M is the mean perigee, p is the semi-major radius of the orbit, μ is the gravitational constant of Earth, θ is the true perigee, t is time, r is the distance from the Earth's center, and f r f u and f h These are the x, y, and z components of the perturbation force in the satellite's first orbital coordinate system, respectively; the parameters of the orbital dynamics equations are the configuration and orbital parameters of the satellite constellation.

[0108] The state transition equations for the satellite constellation are:

[0109]

[0110] in, C represents the number of satellites in the cluster, X k Let Xk-1 be the state variable of all satellites in the cluster at time k, and XtoE be the state variable of all satellites in the cluster at time k-1. k-1 f is the transformation function from Cartesian coordinates to Kepler orbital features. k Let f be the state transition equation at time k. e,k The state transition equation in terms of Kepler orbital elements, W k The state prediction noise is Gaussian white noise with zero mean, E k E represents the Kepler orbital features of all satellites in the cluster at time k. k-1 This represents the Kepler orbital features of all satellites in the cluster at time k-1. Let be the rate of change of the satellite's orbital elements over time; the formula for calculating the rate of change of the satellite's orbital elements over time is:

[0111]

[0112] In subsequent calculations, the calculation of the optimal estimate of the state variable at a certain moment needs to take into account the known state variable storage data for a period of time before the current moment;

[0113] Linearizing the state transition equations yields the state transition equations for the satellite constellation:

[0114] X k =F k (X k-1 )+W k ,

[0115] Among them, X k Let X be the state variable of all satellites in the cluster at time k. k-1 F represents the state variables of all satellites in the cluster at time k-1. k Let W be the state transition matrix of the satellite constellation at time k. k The state prediction noise vector is Gaussian white noise with zero mean;

[0116] Optionally, the state transition functions of all satellites in the cluster can be obtained by combining the state transition function of the reference satellite with the CW equations of relative motion of the slave satellites.

[0117] S022: Establish the joint state transition equations:

[0118] Let the state update vector be... Where k represents the k-th time, W s Indicates a time window. The initial linearization point for iterative computation. This represents the updated state quantity estimate, calculated by combining the joint prediction matrix and the state expectation vector based on the satellite's state transition matrix and the covariance matrix of the state prediction noise at each time point within the time window.

[0119] The squared Mahalanobis distance d is defined as follows:

[0120]

[0121] Where, Σ l Let l be the covariance matrix of column vector l, with the sign of l. It is a norm 2;

[0122] The factors corresponding to the cluster state equation are:

[0123]

[0124] The cluster state matrix within the time window is represented as follows: Selecting the initial linearization point for iterative computation Let the state update vector be Calculated variables:

[0125]

[0126]

[0127] Where, k0 = kW s +1, Q k Let be the covariance matrix of the state prediction noise;

[0128] Get the entire time window W s The joint state transition equation corresponding to the cluster fusion positioning is as follows:

[0129]

[0130] Here, the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where t represents time, k represents the current time of the solution, and W represents the state update vector. s Indicates a time window. This is the joint prediction matrix for each time t within the time window. Let be the state expectation vector at each time t within the time window;

[0131] Step 3: Using the absolute spatial position and absolute velocity of the satellite constellation as observations, establish the observation equations for the absolute state of the satellites. Based on the covariance matrix of the GNSS measurement noise obtained in Step 1, establish the joint absolute observation equations, and calculate the state update quantities for the absolute position observations of individual satellites using the GNSS method, including:

[0132] S031: Obtain the absolute spatial position of the satellite constellation via the onboard GNSS receiver, including pseudorange observations received by the satellites from the navigation satellites. Record the pseudorange observations received by the satellites as... At time k, satellite i receives the pseudorange observation value from navigation satellite number 1. At time k, satellite i receives pseudorange observations from navigation satellite number 2. At time k, satellite i receives a signal numbered N. k pseudorange observations of navigation satellites, N k Let be the number of navigation satellites received by satellite i at time k;

[0133] The observation equations for establishing the absolute state of satellites for GNSS positioning are as follows:

[0134]

[0135] Where i is the satellite number, k is the k-th time, and X k For satellite state variables, h i Let V be the observation function. i k For observing noise, the observation function for satellite i is:

[0136]

[0137]

[0138] Where i is the satellite number, k is the k-th time, and s is the navigation satellite number. Let h be the state variable of satellite i at time k. i Let h be the observation function of satellite i. GNSS(i1) Let N be the observation function of satellite i relative to navigation satellite 1, and so on. k Let k be the number of navigation satellites that the satellite can receive at time k. Let i be the position of satellite i at time k. Let be the position of navigation satellite s at time k;

[0139] Observation noise V k i The covariance matrix is ​​expressed as:

[0140]

[0141] S032: Establish the joint absolute observation equations and calculate the state update quantities of single-satellite absolute position observations using the GNSS method:

[0142] Let the state update vector be... Where k represents the k-th time. The initial linearization point for iterative computation. To obtain the updated state quantity estimates, the joint absolute observation matrix and the absolute observation expectation vector are calculated based on the GNSS observation equations and the covariance matrix of the GNSS measurement noise at each time point within the time window.

[0143] The error factor for absolute position observation is:

[0144]

[0145] Calculate the following variables:

[0146]

[0147]

[0148] in:

[0149]

[0150]

[0151]

[0152]

[0153] Based on the above results, the entire time window W is obtained. s The following joint absolute observation equation corresponds to cluster fusion localization:

[0154]

[0155] Where the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where i represents the satellite number, s represents the navigation satellite number, t represents the time, k represents the current time of the solution, and W... s Indicates a time window. This represents the joint absolute observation matrix of satellite i relative to navigation satellite s at each time t within the time window. Let t be the absolute observation expectation vector of satellite i relative to navigation satellite s at each time point t within the time window;

[0156] Based on the entire time window W s The joint absolute observation equation corresponding to cluster fusion positioning is used to calculate the state update of single-satellite absolute position observations using the GNSS method;

[0157] Step four: Using the spatial relative position and relative velocity of the satellite constellation as observations, establish observation equations for the relative state of the satellites. Based on the covariance matrix of inter-satellite optical observation noise obtained in step one, establish joint relative observation equations, and calculate the state update quantities of inter-satellite relative positions observed through optical cameras, including:

[0158] S041: Observe other satellites in the cluster using the onboard optical camera to obtain relative observation information, including relative distance L. ij and two angles measured only α ij β ij This refers to the spatial relative position and relative velocity of the satellite constellation. The spatial relative position and relative velocity of the satellite constellation are used as the observations, and these observations are... The observation equation for obtaining the relative state of the satellite is:

[0159]

[0160]

[0161] Where i and j are the satellite numbers of two satellites with a relative observation relationship, k is the time k, and X k Let k be the state variable of the satellite at time k. These are the coordinate components of satellite j in the first orbital coordinate system of satellite i at time k;

[0162] Linearize the observation equations and write them in matrix form:

[0163]

[0164] Where i and j are the satellite numbers of two satellites with a relative observation relationship, and k is the time at time k. V is the relative observation matrix. k ij For relative observation noise, X k For the state variables of the satellite constellation;

[0165] S042: Establish the joint relative observation equation and calculate the state update of the inter-satellite relative positions observed by the optical camera;

[0166] Let the state update vector be... Where k represents the k-th time. The initial linearization point for iterative computation. This represents the updated state quantity estimate, calculated based on the GNSS observation equations and the covariance matrix of GNSS measurement noise at each moment within the time window, along with the joint relative observation matrix and the relative observation expectation vector.

[0167] The error factor for relative position observation is:

[0168]

[0169] Calculate the following variables:

[0170]

[0171]

[0172] in This is the covariance matrix of the relative position observation error;

[0173] Get the entire time window W s The joint relative observation equation corresponding to the cluster fusion positioning is as follows:

[0174]

[0175] Here, the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where i and j represent satellite numbers, t represents time, k represents the current solution time, and W represents the current solution time. s Indicates a time window. This represents the joint relative observation matrix of satellite i relative to satellite j at each time t within the time window. Let t be the relative observation expectation vector of satellite i relative to satellite j at each time t within the time window;

[0176] Step 5: Based on the state transition of the satellite cluster and the Bayesian network model of multi-source fusion observations, the joint state transition equation obtained in Step 2, the joint absolute observation equation obtained in Step 3, and the joint relative observation equation obtained in Step 4 are added to obtain the cluster positioning state update equation. The Newton-Raphson method is used to solve the cluster positioning state update equation to obtain the optimal estimate of the satellite state variables at the current moment. The optimal estimate is used as the state update variable to update the satellite cluster state, thereby achieving spatial position and velocity calibration of the satellite cluster and determining the precise orbit of the satellite cluster, including:

[0177] The Bayesian network model based on satellite state and observations is obtained using maximum a posteriori probability (MAP) inference:

[0178]

[0179] Where p represents probability, X 0:k ={X0,X1,...,X k} represents the state of the satellite constellation, Z 1:k ={Z1,Z2,...,Z k} represents the observed values ​​of the satellite's state. in Let represent the nodes that are adjacent nodes with a relative observation relationship to satellite i; the MAP problem in this formula is equivalent to:

[0180] ,

[0181] The optimization problem corresponding to this formula is equivalent to a nonlinear least-squares (NLS) problem. The first term corresponds to the state transition equation, the second term corresponds to the absolute position measurement of a single satellite, and the third term corresponds to the inter-satellite relative measurement; f t It is the state transition equation, h t i and h t ij It is the observation equation, X t It is a state variable, Z t ij and Z t i It is an observation, W s For time windows;

[0182] S051: Add the joint state transition equation, the joint absolute observation equation, and the joint relative observation equation to obtain the cluster positioning state update equation for all times within the time window:

[0183]

[0184] Where, the symbol minarg is the parameter used to find the minimum value of the expression on the right. Let Σ be the norm 2, and let Σ be the summation. Update the state vector. Let W be the state update vector that minimizes the expression, where i and j are satellite numbers, s is the navigation satellite number, t is the time, k is the current time of the solution, and W is the time of the solution. s For time window, This is the joint prediction matrix for each time t within the time window. Let be the state expectation vector at each time t within the time window. This represents the joint absolute observation matrix of satellite i relative to navigation satellite s at each time t within the time window. Let be the absolute observation expectation vector of satellite i relative to navigation satellite s at each time t within the time window. This represents the joint relative observation matrix of satellite i relative to satellite j at each time t within the time window. Let t be the relative observation expectation vector of satellite i relative to satellite j at each time t within the time window;

[0185] S052: The method for iteratively solving the cluster positioning state update equation using Newton's method to obtain the optimal estimate of the satellite state variables at the current moment is as follows: Determine the initial state estimate based on the original observation values, use Newton's optimization algorithm to solve the cluster positioning state update equation, calculate the update vector for the current step, add the update vector to the initial state estimate, update the state estimate, update the cluster positioning state update equation, and repeat the above process; when the difference between the state variable estimates before and after iteration is less than the accuracy requirement required in the specific scenario, or when the number of iterations exceeds the set maximum number of iterations, stop the iteration, and use the final result as the optimal estimate of the satellite state variables at the current moment. Use the optimal estimate as the state update value to update the satellite cluster state, realize the spatial position and velocity calibration of the satellite cluster, and thus determine the precise orbit of the satellite cluster.

[0186] S053: Let k = k + 1, return to step S051, and solve for the optimal estimate of the satellite state variables at the next moment.

[0187] The parts of this invention not described in detail are well-known to those skilled in the art.

Claims

1. A method for precise orbit determination of satellite constellations based on multi-source data fusion, characterized in that, include: Step 1: Obtain basic parameters, acquire satellite cluster orbital state and attitude information, acquire state quantity estimates and observations, and calculate the covariance matrix of state prediction noise, the covariance matrix of GNSS measurement noise, and the covariance matrix of inter-satellite optical observation noise. Step 2: Obtain the historical state data of the satellites. Based on the configuration and orbital parameters of the satellite constellation, establish the state transition equation of the satellite constellation. Based on the covariance matrix of the state prediction noise obtained in Step 1, establish the joint state transition equation. Step 3: Using the absolute position and absolute velocity of the satellite cluster as observations, establish the observation equation for the absolute state of the satellites. Based on the covariance matrix of the GNSS measurement noise obtained in Step 1, establish the joint absolute observation equation and calculate the state update quantity of the absolute position observation of a single satellite using the GNSS method. Step 4: Using the spatial relative position and relative velocity of the satellite cluster as observations, establish the observation equation for the relative state of the satellites. Based on the covariance matrix of the inter-satellite optical observation noise obtained in Step 1, establish the joint relative observation equation and calculate the state update amount of the inter-satellite relative position observed by the optical camera. Step 5: Based on the state transition of the satellite cluster and the Bayesian network model of multi-source fusion observation, add the joint state transition equation obtained in Step 2, the joint absolute observation equation obtained in Step 3, and the joint relative observation equation obtained in Step 4 to obtain the cluster positioning state update equation. Use Newton's method to solve the cluster positioning state update equation to obtain the optimal estimate of the satellite state variables at the current time. Use the optimal estimate as the state update variable to update the satellite cluster state, realize the spatial position and velocity calibration of the satellite cluster, and thus determine the precise orbit of the satellite cluster.

2. The method for precise orbit determination of satellite constellations based on multi-source data fusion according to claim 1, characterized in that, In step one, the basic parameters include the gravitational constant of the central celestial body and the configuration parameters of the satellite cluster. The configuration parameters of the satellite cluster are the orbital elements of each satellite in the cluster, and the central celestial body is the Earth. The orbital status and attitude information of the satellite constellation includes: the orbital status of the reference satellite in the constellation, and the attitude information of each satellite in the constellation; The method for obtaining the state quantity estimates and observations is as follows: Obtain the state quantities within a certain period preceding the position calibration solution time, as the state estimates. The length of this preceding period is determined based on the specific scenario's accuracy requirements. The state estimate for satellite i at time k is... in These represent the first, second, and third position components of the satellite in the geocentric inertial coordinate system. The first, second, and third velocity components of the satellite in the geocentric inertial coordinate system are used; the observation records stored in the satellite's onboard computer are used as the observation values. The methods for calculating the covariance matrix of state prediction noise, the covariance matrix of GNSS measurement noise, and the covariance matrix of inter-satellite optical observation noise are as follows: Based on the satellite orbital state data stored in the satellite's onboard computer, the state prediction noise vector W is obtained. k Then calculate the covariance matrix of the state prediction noise: Where k is the k-th time, Q k Let W be the covariance matrix of the state prediction noise. k For state prediction noise, W k (i) is W k The i-th element in the vector, dim, represents the state variable X. k Number of elements in the middle For mathematical expectation, E((W) k (1)) 2 To calculate the variance of the first element in the state prediction noise, E((W) k (2)) 2 To calculate the variance of the second element in the state prediction noise, E((W) k (dim)) 2 To calculate the variance of the last element in the state prediction noise, E(W) k (1)·W k (2)) and E(W k (2)·W k (1) both calculate the covariance of the first and second elements in the state prediction noise, E(W k (1)·W k (dim)) and E(W k (dim)·W k (1) both refer to calculating the covariance of the first and last elements in the state prediction noise, E(W k (2)·W k (dim)) and E(W k (dim)·W k (2) Both are to find the covariance of the second and last elements in the state prediction noise. The covariance matrix of the state prediction noise is a symmetric matrix. Based on GNSS observation records, the GNSS observation noise was obtained. Vector, calculation The variance of each element s in Where s is the navigation satellite number, the covariance matrix of GNSS measurement noise is obtained: Where i is the satellite number, k represents the k-th time, and N k This represents the number of navigation satellites that can receive signals at the current moment. Represents the covariance matrix of GNSS measurement noise. This represents the variance of the first element in the GNSS observation noise vector. This represents the variance of the second element in the GNSS observation noise vector. This represents the variance of the last element in the GNSS observation noise vector. The covariance matrix of the GNSS measurement noise is a diagonal matrix. Repeat the above process to calculate the covariance matrix of the state prediction noise for each satellite. Based on inter-satellite optical observation records, the observation noise of the direct relative position between satellite i and satellite j was calculated. Calculate the covariance matrix of inter-satellite optical observation noise: Where i and j are the numbers of the two satellites that are being observed relative to each other, and k is the time at which the observation occurs. Let E(·) be the covariance matrix of inter-satellite optical observation noise, and E(·) be the expected value. To find the variance of the first element in the relative position observation noise, To find the variance of the second element in the relative position observation noise, To find the variance of the third element in the relative position observation noise, and Both are used to calculate the covariance of the first and second elements in the relative position observation noise. and Both are used to calculate the covariance of the first and third elements in the relative position observation noise. and The covariance of the second and third elements in the relative position observation noise is calculated. The above process is repeated to calculate the covariance matrix of the relative inter-satellite optical observation noise for each pair of satellites (i,j) with relative observation relationship.

3. The method for precise orbit determination of satellite constellations based on multi-source data fusion according to claim 1, characterized in that, Obtain the historical state variables of the satellites. Based on the configuration and orbital parameters of the satellite constellation, establish the state transition equations for the satellite constellation. Establish a joint state transition equation based on the covariance matrix of the state prediction noise obtained in step one, including: S021: Obtain the historical records of satellite state variables and establish the state transition equations for the satellite constellation; The orbital dynamics equations of a satellite under the influence of perturbation are as follows: Where a is the semi-major axis of the orbit, e is the eccentricity, Ω is the right ascension of the ascending node, i is the orbital inclination, ω is the argument of perigee, M is the mean perigee, p is the semi-major radius of the orbit, μ is the gravitational constant of Earth, θ is the true perigee, t is time, r is the distance from the Earth's center, and f r f u and f h These are the x, y, and z components of the perturbation force in the satellite's first orbital coordinate system, respectively; the parameters of the orbital dynamics equations are the configuration and orbital parameters of the satellite constellation. The state transition equations for the satellite constellation are: in, C represents the number of satellites in the cluster, X k Let X be the state variable of all satellites in the cluster at time k. k-1 Let XtoE(X) be the state variable of all satellites in the cluster at time k-1. k-1 f is the transformation function from Cartesian coordinates to Kepler orbital features. k Let f be the state transition equation at time k. e,k The state transition equation in terms of Kepler orbital elements, W k The state prediction noise is Gaussian white noise with zero mean, E k E represents the Kepler orbital features of all satellites in the cluster at time k. k-1 This represents the Kepler orbital features of all satellites in the cluster at time k-1. The rate of change of the satellite's orbital elements over time; Linearizing the state transition equations yields the state transition equations for the satellite constellation: X k =F k (X k-1 )+W k , Among them, X k Let X be the state variable of all satellites in the cluster at time k. k-1 F represents the state variables of all satellites in the cluster at time k-1. k Let W be the state transition matrix of the satellite constellation at time k. k The noise vector is used to predict the state. S022: Establish the joint state transition equations: Let the state update vector be... Where k represents the k-th time, W s Indicates a time window. The initial linearization point for iterative computation. This represents the updated state quantity estimate, calculated by combining the joint prediction matrix and the state expectation vector based on the satellite's state transition matrix and the covariance matrix of the state prediction noise at each time point within the time window. Get the entire time window W s The joint state transition equation corresponding to the cluster fusion positioning is: Here, the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where t represents time, k represents the current time of the solution, and W represents the state update vector. s Indicates a time window. This is the joint prediction matrix for each time t within the time window. Let be the state expectation vector at each time t within the time window.

4. The method for precise orbit determination of satellite constellations based on multi-source data fusion according to claim 1, characterized in that, Using the absolute spatial position and absolute velocity of the satellite constellation as observations, an observation equation for the absolute state of the satellites is established. Based on the covariance matrix of the GNSS measurement noise obtained in step one, a joint absolute observation equation is established. The state update quantities for the absolute position observations of individual satellites using the GNSS method are calculated, including: S031: Obtain the absolute spatial position of the satellite constellation via the onboard GNSS receiver, including pseudorange observations received by the satellites from the navigation satellites. Record the pseudorange observations received by the satellites as... At time k, satellite i receives the pseudorange observation value from navigation satellite number 1. At time k, satellite i receives pseudorange observations from navigation satellite number 2. At time k, satellite i receives a signal numbered N. k pseudorange observations of navigation satellites, N k Let be the number of navigation satellites received by satellite i at time k; The observation equations for establishing the absolute state of satellites for GNSS positioning are as follows: Where i is the satellite number, k is the k-th time, and X k For satellite state variables, h i For the observation function, For observing noise, the observation function for satellite i is: Where i is the satellite number, k is the k-th time, and s is the navigation satellite number. Let h be the state variable of satellite i at time k. i Let h be the observation function of satellite i. GNSS(i1) Let N be the observation function of satellite i relative to navigation satellite 1, and so on. k Let k be the number of navigation satellites that the satellite can receive at time k. Let i be the position of satellite i at time k. Let be the position of navigation satellite s at time k; S032: Establish the joint absolute observation equations and calculate the state update quantities of single-satellite absolute position observations using the GNSS method: Let the state update vector be... Where k represents the k-th time. The initial linearization point for iterative computation. To obtain the updated state quantity estimates, the joint absolute observation matrix and the absolute observation expectation vector are calculated based on the GNSS observation equations and the covariance matrix of the GNSS measurement noise at each time point within the time window. Based on the above results, the entire time window W is obtained. s The following joint absolute observation equation corresponds to cluster fusion localization: Here, the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where i represents the satellite number, s represents the navigation satellite number, t represents the time, k represents the current time of the solution, and W... s Indicates a time window. This represents the joint absolute observation matrix of satellite i relative to navigation satellite s at each time t within the time window. Let t be the absolute observation expectation vector of satellite i relative to navigation satellite s at each time point t within the time window; Based on the entire time window W s The joint absolute observation equation corresponding to cluster fusion positioning is used to calculate the state update of single-satellite absolute position observations using the GNSS method.

5. The method for precise orbit determination of satellite constellations based on multi-source data fusion according to claim 1, characterized in that, Using the spatial relative position and relative velocity of the satellite constellation as observations, an observation equation for the relative state of the satellites is established. Based on the covariance matrix of the inter-satellite optical observation noise obtained in step one, the relative observation equation is used to calculate the state update quantities of the inter-satellite relative positions observed through the optical camera, including: S041: Observe other satellites in the cluster using the onboard optical camera to obtain relative observation information, including relative distance L. ij and two angles measured only α ij β ij This refers to the spatial relative position and relative velocity of the satellite constellation. The spatial relative position and relative velocity of the satellite constellation are used as the observations, and these observations are... The observation equation for obtaining the relative state of the satellite is: Where i and j are the satellite numbers of two satellites with a relative observation relationship, k is the time k, and X k Let k be the state variable of the satellite at time k. These are the coordinate components of satellite j in the first orbital coordinate system of satellite i at time k; Linearize the observation equations and express them in matrix form: Where i and j are the satellite numbers of two satellites with a relative observation relationship, and k is the time at time k. For relative observation matrices, For relative observation noise, X k For the state variables of the satellite constellation; S042: Establish the joint relative observation equation and calculate the state update of the inter-satellite relative positions observed by the optical camera; Let the state update vector be... Where k represents the k-th time. The initial linearization point for iterative computation. This represents the updated state quantity estimate, calculated based on the observation equations of the satellite's relative state at each moment within the time window and the covariance matrix of the inter-satellite optical observation noise, along with the joint relative observation matrix and the relative observation expectation vector. Get the entire time window W s The following joint relative observation equation corresponds to cluster fusion positioning: Here, the symbol minarg represents the parameter that minimizes the expression on the right. It is a 2-norm, and the symbol Σ represents summation. Represents the state update vector. This represents the state update vector that minimizes the expression, where i and j represent satellite numbers, t represents time, k represents the current solution time, and W represents the current solution time. s Indicates a time window. This represents the joint relative observation matrix of satellite i relative to satellite j at each time t within the time window. Let t be the relative observation expectation vector of satellite i relative to satellite j at each time t within the time window.

6. The method for precise orbit determination of satellite constellations based on multi-source data fusion according to claim 1, characterized in that, Based on the state transition of the satellite constellation and the Bayesian network model of multi-source fusion observations, the joint state transition equation obtained in step two, the joint absolute observation equation obtained in step three, and the joint relative observation equation obtained in step four are added to obtain the constellation positioning state update equation. The Newton-Raphson method is used to solve the constellation positioning state update equation to obtain the optimal estimate of the satellite state variables at the current moment. This optimal estimate is used as the state update variable to update the satellite constellation state, thereby achieving spatial position and velocity calibration of the satellite constellation and determining its precise orbit, including: S051: Add the joint state transition equation, the joint absolute observation equation, and the joint relative observation equation to obtain the cluster positioning state update equation for all times within the time window: Where, the symbol minarg is the parameter used to find the minimum value of the expression on the right. Let Σ be the norm 2, and let Σ be the summation. Update the state vector. Let W be the state update vector that minimizes the expression, where i and j are satellite numbers, s is the navigation satellite number, t is the time, k is the current time of the solution, and W is the time of the solution. s For time window, This is the joint prediction matrix for each time t within the time window. Let be the state expectation vector at each time t within the time window. This represents the joint absolute observation matrix of satellite i relative to navigation satellite s at each time t within the time window. Let be the absolute observation expectation vector of satellite i relative to navigation satellite s at each time t within the time window. This represents the joint relative observation matrix of satellite i relative to satellite j at each time t within the time window. Let t be the relative observation expectation vector of satellite i relative to satellite j at each time t within the time window; S052: The method for iteratively solving the cluster positioning state update equation using Newton's method to obtain the optimal estimate of the satellite state variables at the current moment is as follows: Determine the initial state estimate based on the original observation values, use Newton's optimization algorithm to solve the cluster positioning state update equation, calculate the update vector for the current step, add the update vector to the initial state estimate, update the state estimate, update the cluster positioning state update equation, and repeat the above process; when the difference between the state variable estimates before and after iteration is less than the accuracy requirement required in the specific scenario, or when the number of iterations exceeds the set maximum number of iterations, stop the iteration, and use the final result as the optimal estimate of the satellite state variables at the current moment. Use the optimal estimate as the state update value to update the satellite cluster state, realize the spatial position and velocity calibration of the satellite cluster, and thus determine the precise orbit of the satellite cluster. S053: Let k = k + 1, return to step S051, and solve for the optimal estimate of the satellite state variables at the next moment.

Images