Satellite interlink autonomous orbit determination method and system based on square root filter

Abstract

The application discloses a satellite inter-satellite link autonomous orbit determination method and system based on square root filtering, and relates to the technical field of satellite navigation. The method comprises the following steps: combining each of N satellites in a navigation system with each of the other satellites in the navigation system except the current satellite to obtain N*(N-1) satellite combinations; initializing the state parameters of a preset orbit of each satellite in the navigation system; constructing an inter-satellite ranging observation equation, a first satellite dynamic equation and a second satellite dynamic equation of each satellite combination; constructing a satellite state transition equation of the current satellite based on the state parameters of the preset orbit of the current satellite and the first satellite dynamic equation; and performing measurement updating and time updating on each satellite in the navigation system as the current satellite to complete the orbit determination of each satellite. The application improves the autonomous orbit determination accuracy and stability of the satellite inter-satellite link.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation technology, and in particular to a method and system for autonomous orbit determination of inter-satellite links based on square root filtering. Background Technology

[0002] Global Navigation Satellite Systems (GNSS) serve as crucial infrastructure providing real-time, precise location and time information to Earth's surface users, finding widespread application in critical military and civilian sectors such as defense, transportation, communications, power, and finance. Traditional GNSS systems consist of a ground control system, a satellite system, and a user system. They primarily utilize ground control system monitoring stations to acquire observational information, which is then processed by the ground control system's master control station to generate necessary navigation ephemeris parameters, such as orbits and clock biases, for the satellites. These parameters are then uploaded to the satellites, which relay them to users to provide navigation and positioning services. When the ground control system is destroyed due to natural or human-caused damage, the GNSS system loses its service capability. To address this, GNSS systems like GPS and BeiDou have added inter-satellite link payloads and enhanced onboard data processing capabilities. This allows GNSS systems to perform autonomous orbit determination and time synchronization data processing in orbit, and to autonomously update navigation ephemeris parameters, thus providing them with an autonomous navigation capability independent of the ground control system.

[0003] When using satellites for autonomous orbit determination and time synchronization, the computational power, storage capacity, and on-orbit reliability of onboard processors are significantly lower than those of ground-based processors. The centralized processing model typically used in ground-based operation and control systems cannot be directly applied to satellites, necessitating the design of new algorithms to simplify the computational load on a single satellite. To address this, the BeiDou Navigation Satellite System employs an approximate distributed filtering processing mode, referred to as the equivalent observation distributed processing method, in its autonomous navigation implementation. In this method, after acquiring inter-satellite ranging observations from linkable satellites, each satellite refines its position, velocity, and other state parameters and covariance information through Kalman filtering measurement updates and time updates. During Kalman filtering measurement updates, each satellite uses the predicted orbits of its linked satellites as known parameters, the predicted orbit positions of the linked satellites as reference values, and the orbital position covariance information as equivalent measurement errors in the autonomous orbit determination data processing. The aforementioned distributed filtering algorithm is simple to calculate and easy to implement. When each satellite can acquire a large number of inter-satellite link observations, it can obtain autonomous orbit determination results with a certain level of accuracy. However, the shortcomings of the above equivalent observation distributed processing method are also obvious. Since the covariance information between the measurable satellite state parameters is completely ignored during the filtering process, the filtering results have high requirements on the satellite measurement geometry. Especially when each satellite acquires few inter-satellite link observations, the filtering accuracy is poor and the results are difficult to converge stably.

[0004] To address the issue that equivalent observation distributed processing methods do not adequately consider the covariance information of satellite state parameters, P. A. Ferguson proposed a Schmidt-Kalman distributed filtering method. In this method, all satellite state parameters are divided into two categories: the satellite state parameters that actually need to be estimated (each satellite's own state parameters) and the associated parameters that do not need to be estimated (the state parameters of other satellites). During the Schmidt-Kalman distributed filtering measurement update data processing, each satellite only estimates and updates its own state parameters, while only considering the covariance information of associated satellite parameters. The Schmidt-Kalman distributed filtering method comprehensively considers the influence of covariance information among satellite state parameters. Compared to equivalent observation distributed processing, which only considers the satellite's own covariance and ignores other information, it uses more comprehensive covariance information, theoretically has better parameter estimation accuracy, and produces more stable estimation results. However, when the Schmidt-Kalman filtering method is directly applied to the processing of autonomous orbit determination data from inter-satellite measurements, the near-rank deficiency of the observation equations for autonomous orbit determination causes the filter covariance matrix to easily lose its positive definiteness, leading to rapid divergence of the Schmidt-Kalman distributed filtering method. Although adding auxiliary ground station observations can improve the rank deficiency of the observation equations to some extent and thus slow down the filter divergence rate, it cannot fundamentally solve the problem.

[0005] Therefore, existing autonomous orbit determination methods for inter-satellite links not only have low orbit determination accuracy but also poor stability. Summary of the Invention

[0006] The purpose of this invention is to provide a method and system for autonomous orbit determination of inter-satellite links based on square root filtering, which improves the accuracy and stability of autonomous orbit determination of inter-satellite links.

[0007] To achieve the above objectives, the present invention provides the following solution:

[0008] A method for autonomous orbit determination of satellites via inter-satellite links based on square root filtering, the method comprising:

[0009] Each of the N satellites in the navigation system is paired with a satellite in the navigation system other than the current satellite to obtain N×(N-1) satellite combinations; each satellite combination includes the current satellite and the link-establishing satellite, and each of the N satellites in the navigation system is considered as a current satellite, with one current satellite corresponding to N-1 satellite combinations; N>1;

[0010] Initialize the state parameters of the preset orbit of each satellite in the navigation system; the state parameters include position, velocity, and solar radiation pressure parameters;

[0011] For each satellite combination, construct the inter-satellite ranging observation equation, the first satellite dynamics equation, and the second satellite dynamics equation; the first satellite dynamics equation is the satellite dynamics equation of this satellite, and the second satellite dynamics equation is the satellite dynamics equation of the linked satellite; the inter-satellite ranging observation equation is a function of the actual state parameters of this satellite, the state parameters of the preset orbit corresponding to this satellite, the actual state parameters of the linked satellite, and the state parameters of the preset orbit corresponding to the linked satellite; the satellite dynamics equation is a function of the actual state parameters of the satellite, Earth's gravity, Sun's gravity, Moon's gravity, solar radiation pressure, and residual perturbation force;

[0012] Based on the state parameters of the satellite's preset orbit and the first satellite dynamic equation, the satellite state transition equation is constructed.

[0013] Each satellite in the navigation system is treated as the local satellite for measurement and time updates to complete the orbit determination of each satellite;

[0014] For any given satellite, measurement and time updates are performed to complete the satellite's orbit determination, specifically including:

[0015] Determine if the current time is the preset final time;

[0016] If so, the final state parameters of the satellite are determined by measuring and updating the correction amount based on the N-1 state parameters of the satellite at the current moment and the state parameters of the satellite's preset orbit, thus completing the orbit determination of the satellite.

[0017] If not, at the current time, the Kalman filter method is used to update the measurements of each of the N-1 satellite combinations corresponding to the current satellite, to obtain the current state parameter measurement update correction amount and the current total covariance information square root matrix update value of the current satellite. Then, the Kalman filter method is used to update the time of each of the N-1 satellite combinations corresponding to the current satellite based on the satellite state transition equation, until the current time is the preset final time.

[0018] Specifically, at any given time, measurements are updated for any N-1 satellite combinations corresponding to the current satellite, including:

[0019] Based on the state parameters of the preset orbits of the satellite and the preset orbits of the linked satellites in the d-th satellite combination, a total covariance information square root matrix is ​​constructed for the d-th satellite combination. The total covariance information square root matrix is ​​a matrix relating to a first autocovariance information square root matrix, a second autocovariance information square root matrix, and a cross-covariance information square root matrix. The first autocovariance information square root matrix is ​​the autocovariance information square root matrix of the state parameters of the satellite; the second autocovariance information square root matrix is ​​the autocovariance information square root matrix of the state parameters of the linked satellite; and the cross-covariance information square root matrix is ​​the cross-covariance information square root matrix of the state parameters of the satellite and the linked satellite. d = 1, 2, ..., N-1.

[0020] Based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite and the state parameters of the corresponding preset orbit of this satellite, determine the measurement update correction amount of the d-th state parameter of this satellite;

[0021] Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, determine the updated value of the square root matrix of the total covariance information of the d-th satellite combination;

[0022] Determine if d is less than N-1;

[0023] If so, the measurement update correction amount of the d-th state parameter is determined as the state parameter of the preset orbit of the satellite in the (d+1)-th satellite combination, and the d-th satellite combination is replaced with the (d+1)-th satellite combination. The result is "Construct the square root matrix of the total covariance information of the d-th satellite combination based on the state parameters of the preset orbit of the satellite in the d-th satellite combination and the state parameters of the preset orbit of the linked satellite".

[0024] If not, then the updated correction values ​​of the N-1 state parameters of this satellite and the updated value of the square root matrix of the total covariance information of the N-1 satellite combinations corresponding to this satellite are obtained.

[0025] Specifically, for time e, time updates are performed on the N-1 satellite combinations corresponding to each current satellite, where e > 1. This includes:

[0026] Using the satellite state transition equation, the state parameters at time e are determined based on the measurement update corrections of the N-1 state parameters of the current satellite at time e-1 and the state parameters of the current satellite's preset orbit; wherein, the measurement update corrections of the N-1 state parameters of the current satellite at the initial time are the preset corrections.

[0027] Optionally, based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit of this satellite, the measurement update correction amount of the d-th state parameter of this satellite is determined, specifically including:

[0028] Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, the first square root filtering intermediate vector, the second square root filtering intermediate vector, and the measurement weight inverse of the d-th satellite combination are determined; the first square root filtering intermediate vector is the square root filtering intermediate vector of this satellite, and the second square root filtering intermediate vector is the square root filtering intermediate vector of the linked satellite.

[0029] The measurement update correction amount of the d-th state parameter of this satellite is determined based on the square root matrix of the total covariance information of the d-th satellite combination, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the inverse of the measurement weight, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit.

[0030] Optionally, based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, the updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined, specifically including:

[0031] Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, the first square root filter intermediate vector, the second square root filter intermediate vector, the measurement weight inverse, and the square root filter covariance attenuation coefficient of the d-th satellite combination are determined; the first square root filter intermediate vector is the square root filter intermediate vector of this satellite, and the second square root filter intermediate vector is the square root filter intermediate vector of the link-established satellite.

[0032] The updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined based on the square root matrix of the total covariance information, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the inverse of the measurement weight, and the square root filter covariance attenuation coefficient.

[0033] Optionally, the measurement update correction amount for the d-th state parameter of the satellite is determined based on the square root matrix of the total covariance information of the d-th satellite combination, the first square root filter intermediate vector, the second square root filter intermediate vector, the measurement weight inverse, the actual state parameters of the satellite, and the state parameters of the corresponding preset orbit. Specifically, this includes:

[0034] Based on the square root matrix of the total covariance information of the d-th satellite combination, the first square root filter intermediate vector, the second square root filter intermediate vector, and the measurement weight inverse, the state parameter update gain matrix of the d-th satellite combination is determined; the state parameter update gain matrix includes a first state parameter update gain matrix and a second state parameter update gain matrix, the first state parameter update gain matrix is ​​the state parameter update gain matrix of this satellite, and the second state parameter update gain matrix is ​​the state parameter update gain matrix of the linked satellite;

[0035] The measurement update correction amount of the d-th state parameter of this satellite is determined based on the actual state parameters of this satellite, the state parameters of the corresponding preset orbit, and the state parameter update gain matrix of the d-th satellite combination.

[0036] Optionally, the updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined based on the square root matrix of the total covariance information, the first square root filter intermediate vector, the second square root filter intermediate vector, the inverse of the measurement weights, and the square root filter covariance attenuation coefficient. Specifically, this includes:

[0037] Based on the square root matrix of the total covariance information of the d-th satellite combination, the first square root filter intermediate vector, the second square root filter intermediate vector, and the measurement weight inverse, the state parameter update gain matrix of the d-th satellite combination is determined; the state parameter update gain matrix includes a first state parameter update gain matrix and a second state parameter update gain matrix, the first state parameter update gain matrix is ​​the state parameter update gain matrix of this satellite, and the second state parameter update gain matrix is ​​the state parameter update gain matrix of the linked satellite.

[0038] The updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined based on the square root matrix of the total covariance information, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the square root filter covariance attenuation coefficient, and the state parameter update gain matrix.

[0039] An autonomous orbit determination system for inter-satellite links based on square root filtering, the system comprising:

[0040] The satellite combination determination module is used to combine each of the N satellites in the navigation system with other satellites in the navigation system except the current satellite, to obtain N×(N-1) satellite combinations; each satellite combination includes the current satellite and the link-establishing satellite, and the N satellites in the navigation system are all considered as the current satellite, and one current satellite corresponds to N-1 satellite combinations; N>1;

[0041] The preset orbit initialization module is used to initialize the state parameters of the preset orbit of each satellite in the navigation system; the state parameters include position, velocity, and solar radiation pressure parameters;

[0042] The first equation construction module is used to construct the inter-satellite ranging observation equation, the first satellite dynamics equation, and the second satellite dynamics equation for each satellite combination. The first satellite dynamics equation is the satellite dynamics equation for this satellite, and the second satellite dynamics equation is the satellite dynamics equation for the linked satellite. The inter-satellite ranging observation equation is a function of the actual state parameters of this satellite, the state parameters of the preset orbit corresponding to this satellite, the actual state parameters of the linked satellite, and the state parameters of the preset orbit corresponding to the linked satellite. The satellite dynamics equation is a function of the actual state parameters of the satellite, Earth's gravity, Sun's gravity, Moon's gravity, solar radiation pressure, and residual perturbation force.

[0043] The second equation construction module is used to construct the satellite state transition equations based on the state parameters of the satellite's preset orbit and the first satellite dynamic equations.

[0044] The orbit determination module is used to perform measurement updates and time updates on each satellite in the navigation system as if it were the local satellite, thereby completing the orbit determination of each satellite.

[0045] Specifically, in terms of performing measurement and time updates for any given satellite to complete its orbit determination, the orbit determination module is used for:

[0046] Determine if the current time is the preset final time;

[0047] If so, the final state parameters of the satellite are determined by measuring and updating the correction amount based on the N-1 state parameters of the satellite at the current moment and the state parameters of the satellite's preset orbit, thus completing the orbit determination of the satellite.

[0048] If not, at the current time, the Kalman filter method is used to update the measurements of each of the N-1 satellite combinations corresponding to the current satellite, to obtain the current state parameter measurement update correction amount and the current total covariance information square root matrix update value of the current satellite. Then, the Kalman filter method is used to update the time of each of the N-1 satellite combinations corresponding to the current satellite based on the satellite state transition equation, until the current time is the preset final time.

[0049] Specifically, at any given time, measurements are updated for any N-1 satellite combinations corresponding to the current satellite, including:

[0050] Based on the state parameters of the preset orbits of the satellite and the preset orbits of the linked satellites in the d-th satellite combination, a total covariance information square root matrix is ​​constructed for the d-th satellite combination. The total covariance information square root matrix is ​​a matrix relating to a first autocovariance information square root matrix, a second autocovariance information square root matrix, and a cross-covariance information square root matrix. The first autocovariance information square root matrix is ​​the autocovariance information square root matrix of the state parameters of the satellite; the second autocovariance information square root matrix is ​​the autocovariance information square root matrix of the state parameters of the linked satellite; and the cross-covariance information square root matrix is ​​the cross-covariance information square root matrix of the state parameters of the satellite and the linked satellite. d = 1, 2, ..., N-1.

[0051] Based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite and the state parameters of the corresponding preset orbit of this satellite, determine the measurement update correction amount of the d-th state parameter of this satellite;

[0052] Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, determine the updated value of the square root matrix of the total covariance information of the d-th satellite combination;

[0053] Determine if d is less than N-1;

[0054] If so, the measurement update correction amount of the d-th state parameter is determined as the state parameter of the preset orbit of the satellite in the (d+1)-th satellite combination, and the d-th satellite combination is replaced with the (d+1)-th satellite combination. The result is "Construct the square root matrix of the total covariance information of the d-th satellite combination based on the state parameters of the preset orbit of the satellite in the d-th satellite combination and the state parameters of the preset orbit of the linked satellite".

[0055] If not, then the updated correction values ​​of the N-1 state parameters of this satellite and the updated value of the square root matrix of the total covariance information of the N-1 satellite combinations corresponding to this satellite are obtained.

[0056] Specifically, for time e, time updates are performed on the N-1 satellite combinations corresponding to each current satellite, where e > 1. This includes:

[0057] Using the satellite state transition equation, the state parameters at time e are determined based on the measurement update corrections of the N-1 state parameters of the current satellite at time e-1 and the state parameters of the current satellite's preset orbit; wherein, the measurement update corrections of the N-1 state parameters of the current satellite at the initial time are the preset corrections.

[0058] According to specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0059] This invention discloses a satellite inter-satellite link autonomous orbit determination method and system based on square root filtering. It utilizes the square root matrix of the self-covariance information and the square root matrix of the cross-covariance information in the satellite combination for measurement and time updates, realizing distributed autonomous orbit determination based on inter-satellite measurements. It solves the problem that the covariance matrix is ​​prone to lose positive definiteness in the traditional SCHMIDT distributed KALMAN filtering calculation process due to the near-rank deficiency of the inter-satellite ranging observation equation, which causes the filtering results to easily diverge. This invention significantly improves the accuracy and stability of data processing for distributed autonomous orbit determination of inter-satellite links. Attached Figure Description

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

[0061] Figure 1 This is a schematic diagram of the autonomous orbit determination method for inter-satellite links based on square root filtering provided in Embodiment 1 of the present invention;

[0062] Figure 2 A schematic diagram of the SCHMIDT distributed autonomous orbit determination filtering algorithm improved with covariance square root information;

[0063] Figure 3 The data flow diagram of the SCHMIDT distributed autonomous orbit determination filtering algorithm improved with covariance square root information. Detailed Implementation

[0064] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0065] The purpose of this invention is to provide a method and system for autonomous orbit determination of inter-satellite links based on square root filtering, aiming to improve the accuracy and stability of autonomous orbit determination of inter-satellite links.

[0066] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0067] Example 1

[0068] Figure 1This is a schematic diagram of the autonomous orbit determination method for inter-satellite links based on square root filtering provided in Embodiment 1 of the present invention. Figure 1 As shown, the autonomous orbit determination method for inter-satellite links based on square root filtering in this embodiment includes:

[0069] Step 101: Combine each of the N satellites in the navigation system with the satellites in the navigation system other than the current satellite, to obtain N×(N-1) satellite combinations; each satellite combination includes the current satellite and the link-establishing satellite. The N satellites in the navigation system are all considered as the current satellite, and one current satellite corresponds to N-1 satellite combinations; N>1.

[0070] Step 102: Initialize the state parameters of the preset orbit of each satellite in the navigation system; the state parameters include position, velocity and solar radiation pressure parameters.

[0071] Step 103: Construct the inter-satellite ranging observation equation, the first satellite dynamics equation, and the second satellite dynamics equation for each satellite combination; the first satellite dynamics equation is the satellite dynamics equation for this satellite, and the second satellite dynamics equation is the satellite dynamics equation for the linked satellite; the inter-satellite ranging observation equation is a function of the actual state parameters of this satellite, the state parameters of the preset orbit corresponding to this satellite, the actual state parameters of the linked satellite, and the state parameters of the preset orbit corresponding to the linked satellite; the satellite dynamics equation is a function of the actual state parameters of the satellite, Earth's gravity, Sun's gravity, Moon's gravity, solar radiation pressure, and residual perturbation force.

[0072] Specifically, for any satellite combination, the inter-satellite ranging observation equation is:

[0073] ρ ij =[(x i -x j ) 2 +(y i -y j ) 2 +(z i -z j ) 2 ] 1 / 2 .

[0074] Where, ρ ij The measured inter-satellite distance between this satellite and the linked satellite is [x] i ,y i ,z i [x] represents the actual position of this satellite (i.e., satellite i), [x] j ,y j ,z j [ ] represents the actual position of the link-establishing satellite (i.e., satellite j).

[0075] For any satellite, the satellite dynamics equations are:

[0076]

[0077] in, The actual position of the satellite (e.g., when the satellite is this satellite (i.e., satellite i)). ), where t is time, The gravitational pull of Earth on the satellite, The gravitational pull of the sun on the satellite, The gravitational pull of the moon on the satellite, The solar radiation pressure experienced by the satellite The perturbations acting on the satellite (such as forces that are not yet known to humankind, such as relativity and tidal forces).

[0078] Step 104: Based on the state parameters of the satellite's preset orbit and the first satellite dynamic equation, construct the satellite state transition equation.

[0079] Specifically, each satellite utilizes pre-given a priori dynamic state parameters such as initial satellite position, velocity, and solar radiation pressure parameters, as well as information such as satellite mass and planetary ephemeris, to solve the satellite dynamic equations using numerical integration methods. This allows for the calculation of the state transition matrix between the preset orbit and the satellite's position, velocity, and other state parameters relative to the preset orbital state parameters.

[0080] For this satellite, the solution to the satellite dynamics equations based on prior information can be expressed as:

[0081]

[0082] in, Let be the prior position parameters of satellite i (i.e., the position of satellite i's preset orbit). Let p be the prior velocity parameter of satellite i (i.e., the velocity of satellite i's preset orbit). i Let be the prior solar radiation pressure parameters of satellite i (i.e., the solar radiation pressure parameters of the preset orbit of satellite i).

[0083] The linearized form of the inter-satellite ranging observation equation obtained by numerical integration using the pre-defined orbit of satellite i is as follows:

[0084]

[0085] Where, Δρ ij The difference between the measured intersatellite spacing and the intersatellite spacing calculated based on the satellite's preset orbit using the integration of prior dynamic state parameters is given. To calculate the interstellar distance using prior dynamic state parameters, This is the correction amount for the actual state parameters of satellite i relative to the prior state parameters corresponding to the preset orbit. H represents the correction amount of the actual state parameters of satellite j relative to the prior state parameters corresponding to the preset orbit, ε is the measurement noise, and H is the value of H. i Let H be the coefficient matrix of the observation equation for satellite i. j Let be the coefficient matrix of the observation equation for satellite j, specifically...

[0086] Using the pre-defined orbit and satellite dynamics equations, the linearized satellite state transition equations can be obtained as follows:

[0087]

[0088] in, Let $\mathbf{k}$ be the correction amount for the satellite state parameters of satellite $i$ at time $k$. This represents the correction amount for the satellite state parameters at time k-1. Let ω be the state transition matrix of satellite i. i Let be the dynamic noise matrix of satellite i.

[0089] Step 105: Perform measurement and time updates on each satellite in the navigation system as if it were the local satellite, and complete the orbit determination for each satellite.

[0090] For any given satellite, measurement and time updates are performed to complete the satellite's orbit determination, specifically including:

[0091] Determine whether the current time is the preset final time.

[0092] If so, the final state parameters of the satellite are determined by measuring and updating the correction amount of the N-1 state parameters of the satellite at the current moment and the state parameters of the satellite's preset orbit, thus completing the orbit determination of the satellite.

[0093] If not, at the current time, the Kalman filter method is used to update the measurements of each of the N-1 satellite combinations corresponding to the current satellite, to obtain the current state parameter measurement update correction and the current total covariance information square root matrix update value of the current satellite. Then, the Kalman filter method is used to update the time of each of the N-1 satellite combinations corresponding to the current satellite based on the satellite state transition equation, until the current time is the preset final time.

[0094] Specifically, at any given time, measurements are updated for any N-1 satellite combinations corresponding to the current satellite, including:

[0095] Based on the state parameters of the preset orbits of the satellite and the preset orbits of the linked satellites in the d-th satellite combination, a total covariance information square root matrix is ​​constructed for the d-th satellite combination. The total covariance information square root matrix is ​​a matrix relating to the first self-covariance information square root matrix, the second self-covariance information square root matrix, and the cross-covariance information square root matrix. The first self-covariance information square root matrix is ​​the self-covariance information square root matrix of the state parameters of the satellite, the second self-covariance information square root matrix is ​​the self-covariance information square root matrix of the state parameters of the linked satellite, and the cross-covariance information square root matrix is ​​the cross-covariance information square root matrix of the state parameters of the satellite and the linked satellite. d = 1, 2, ..., N-1.

[0096] Specifically, utilizing the inter-satellite communication function of the inter-satellite link payload, after two satellites in a satellite constellation complete an inter-satellite measurement communication, each satellite can obtain the latest corrections to its state parameters, such as position, velocity, and dynamic model parameters, as well as the corresponding square root covariance information matrix of the satellite with which it established the inter-satellite measurement link. The square root covariance information matrix of the state parameters of the two satellites corresponding to a satellite constellation that achieves inter-satellite measurement adopts a lower triangular matrix form:

[0097]

[0098] in, The square root matrix of the total covariance information; This is the correction amount for the satellite state parameters of satellite i. Let be the correction amount for the prior position parameters of satellite i. Δp is the correction amount for the prior velocity parameters of satellite i. i This is the correction amount for the prior solar radiation pressure parameters of satellite i; Let J be the correction amount for the satellite state parameters of satellite j. Let be the correction amount for the prior position parameters of satellite j. Δp is the correction amount for the prior velocity parameters of satellite j. j Δ is the correction factor for the prior solar radiation pressure parameters of satellite j. i Let Δ be the square root matrix of the autocovariance information of satellite i (i.e., the square root matrix of the first autocovariance information); j Let Δ be the square root matrix of the autocovariance information of satellite j (i.e., the square root matrix of the second autocovariance information); i，j Let be the square root matrix of the cross-covariance information of satellite i and satellite j. In single-satellite distributed filtering, each satellite only uses the state parameter correction amount of its own satellite. Δp i As a parameter to be estimated, the autocovariance Δ of the satellite's state parameters is also updated. iThe square root matrix Δ of the cross-covariance information between the state parameters of this satellite and the linked satellites ij The state parameter corrections and autocovariance information matrix Δ of the linked satellite j The filtering process uses the forecast value updated after the last filtering time, and does not perform filtering on this satellite.

[0099] Based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit of this satellite, determine the measurement update correction amount of the d-th state parameter of this satellite.

[0100] Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, determine the updated value of the square root matrix of the total covariance information of the d-th satellite combination.

[0101] Determine if d is less than N-1.

[0102] If so, the measurement update correction amount of the d-th state parameter is determined as the state parameter of the preset orbit of the satellite in the (d+1)-th satellite combination, and the d-th satellite combination is replaced with the (d+1)-th satellite combination. The result is "Construct the square root matrix of the total covariance information of the d-th satellite combination based on the state parameters of the preset orbit of the satellite in the d-th satellite combination and the state parameters of the preset orbit of the linked satellite".

[0103] If not, then the updated values ​​of the N-1 state parameters of this satellite and the updated values ​​of the square root matrix of the total covariance information of the N-1 satellite combinations corresponding to this satellite are obtained.

[0104] Specifically, for time e, time updates are performed on the N-1 satellite combinations corresponding to each current satellite, where e > 1. This includes:

[0105] Using the satellite state transition equation, the state parameters at time e are determined based on the measurement update corrections of the N-1 state parameters of the current satellite at time e-1 and the state parameters of the current satellite's preset orbit; wherein, the measurement update corrections of the N-1 state parameters of the current satellite at the initial time are the preset corrections.

[0106] Specifically, through processing, each satellite can acquire the filtered measurement update results of satellite state variables and their autocovariance and cross-covariance matrices for a single inter-satellite measurement link. Since each satellite can establish inter-satellite measurement links with multiple satellites in each processing batch, each satellite needs to sequentially perform the above steps on each measurement link in its batch. After traversing all inter-satellite ranging data for the current satellite and batch, the measurement update processing of all inter-satellite link observations for a single satellite and single processing batch can be achieved, obtaining the corrected satellite state variables and covariance information matrix for all observation data in the current satellite and batch. Each satellite in the entire constellation processes data simultaneously in a distributed mode, enabling the measurement update of each satellite's state parameters and covariance.

[0107] After implementing distributed measurement and update processing for the entire satellite constellation, KALMAN filtering time update processing is required for the state variables and corresponding covariance information of all satellites in the constellation. This allows for extrapolation to obtain prior values ​​of the filter state variables and covariance information for the next processing batch, thus achieving sequential processing. The calculation formula for the time update of all satellite filter state variables is as follows:

[0108]

[0109] in, Let i be the state transition matrix of satellite i. The updated state variable correction is measured for the satellite at time k-1. This is the state quantity correction amount after the time update at time k.

[0110] The formula for calculating the time update of filter state quantity covariance information is as follows:

[0111]

[0112] in:

[0113]

[0114]

[0115]

[0116]

[0117] Where Δ is the square root matrix of the covariance of the entire constellation state vector after time updates. Let ω be the square root matrix of the covariance of the entire constellation's state vector after measurement updates, and let ω be the entire constellation's dynamic noise matrix. For each satellite, ω is the state transition matrix. i Let i be the dynamic noise matrix of satellite i. To filter and measure the updated covariance matrix and Δ of the state parameters of satellite i and satellite ji,j Let T be the covariance matrix of the state parameters of satellite i and satellite j after the filtering time update, and T be the transpose.

[0118] Considering that the square root of the covariance matrix is ​​a triangular matrix, in actual calculations, it is not necessary to calculate all covariances according to the above formula. Only the column HOUSEHOLDER transformation of the following matrix is ​​needed to obtain the square root matrix Δ of the covariance information after the filter time update.

[0119]

[0120] Where Q is the HOUSEHOLDER orthogonal transformation matrix.

[0121] Since the satellite state transition matrix is ​​a diagonal matrix, the above covariance matrix calculation process can be decomposed for each satellite. Each satellite only needs to perform state transition calculations on its own autocovariance and cross-covariance matrices and implement a Householder orthogonal transformation to obtain the square root matrix of the autocovariance and cross-covariance information of its own related satellite state parameters after the corresponding time update.

[0122] As an optional implementation, based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit of this satellite, the measurement update correction amount of the d-th state parameter of this satellite is determined, specifically including:

[0123] Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, the first square root filtering intermediate vector, the second square root filtering intermediate vector, and the measurement weight inverse of the d-th satellite combination are determined; the first square root filtering intermediate vector is the square root filtering intermediate vector of this satellite, and the second square root filtering intermediate vector is the square root filtering intermediate vector of the linked satellite.

[0124] The measurement update correction amount of the d-th state parameter of this satellite is determined based on the square root matrix of the total covariance information of the d-th satellite combination, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the inverse of the measurement weight, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit.

[0125] As an optional implementation, the updated value of the total covariance information square root matrix of the d-th satellite combination is determined based on the total covariance information square root matrix of the d-th satellite combination and the inter-satellite ranging observation equation, specifically including:

[0126] Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, the first square root filter intermediate vector, the second square root filter intermediate vector, the measurement weight inverse, and the square root filter covariance attenuation coefficient of the d-th satellite combination are determined; the first square root filter intermediate vector is the square root filter intermediate vector of this satellite, and the second square root filter intermediate vector is the square root filter intermediate vector of the link-established satellite.

[0127] The updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined based on the square root matrix of the total covariance information, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the inverse of the measurement weight, and the square root filter covariance attenuation coefficient.

[0128] Specifically, after each batch of inter-satellite measurement communication is completed, the distributed filtering measurement update processing procedure is initiated. The filtering measurement update processing completes the processing of a single link at a time. All inter-satellite ranging observations for this satellite in a single batch are sequentially processed link by link until the ranging data of all inter-satellite links in this batch are processed. For each measurement link, the square root filtering intermediate vector, the inverse of the inter-satellite measurement weights, and the filtering covariance attenuation coefficient for each link are calculated first. The square root filtering intermediate vector corresponds to a pair of inter-satellite ranging observations. The calculation formula is:

[0129]

[0130] H=(H i H j ), Let be the filtered intermediate vector of satellite i (i.e., the first square root filtered intermediate vector). This is the intermediate filtering vector of satellite j (i.e., the second square root intermediate filtering vector).

[0131] Single-link inter-satellite measurement weight inverse q ij The calculation formula is:

[0132]

[0133] Where b is the covariance of the measurement noise ε.

[0134] Single-link square root filter covariance attenuation coefficient γ ij The calculation formula is:

[0135]

[0136] As an optional implementation, the measurement update correction amount for the d-th state parameter of the satellite is determined based on the square root matrix of the total covariance information of the d-th satellite combination, the first square root filtered intermediate vector, the second square root filtered intermediate vector, the measurement weight inverse, the actual state parameters of the satellite, and the state parameters of the corresponding preset orbit. Specifically, this includes:

[0137] Based on the square root matrix of the total covariance information of the d-th satellite combination, the first square root filter intermediate vector, the second square root filter intermediate vector, and the measurement weight inverse, the state parameter update gain matrix of the d-th satellite combination is determined. The state parameter update gain matrix includes the first state parameter update gain matrix and the second state parameter update gain matrix. The first state parameter update gain matrix is ​​the state parameter update gain matrix of this satellite, and the second state parameter update gain matrix is ​​the state parameter update gain matrix of the linked satellite.

[0138] Specifically, the obtained single-link square root filter intermediate vector is used. By combining the square root lower triangular matrix of the covariance information of the linked satellites obtained from inter-satellite communication, the gain matrix K for updating satellite state parameters of single-link measurement KALMAN filtering can be calculated. The calculation formula is as follows:

[0139]

[0140] Among them, K i Let K be the filter gain matrix (i.e., the first state parameter update gain matrix) corresponding to the state variables to be estimated for satellite i. j These are the filter gain matrices (i.e., the second state parameter update gain matrices) corresponding to the state variables to be estimated for satellite j.

[0141] The measurement update correction amount of the d-th state parameter of this satellite is determined based on the actual state parameters of this satellite, the state parameters of the corresponding preset orbit, and the state parameter update gain matrix of the d-th satellite combination.

[0142] Specifically, satellite i updates the gain matrix K using the acquired KALMAN filtered satellite state parameters. i (i.e., the first state parameter update gain matrix), combined with the measured inter-satellite distance obtained from the inter-satellite link observation equation, the difference Δρ between the measured inter-satellite distance and the prior parameters is calculated. ij Based on the prior values ​​of the satellite state parameter corrections, the single-link satellite state quantity measurement update correction can be calculated. The calculation formula is as follows:

[0143]

[0144] in, This refers to the measurement update correction amount for the state parameters of satellite i in this single-link measurement correction. It should be noted that in the distributed filtering measurement update processing, the state parameters of the linked satellites corresponding to this satellite are only used as reference values ​​in this satellite's processing, and their state parameters are not measured and updated in this satellite's processing.

[0145] As an optional implementation, the updated value of the total covariance information square root matrix of the d-th satellite combination is determined based on the total covariance information square root matrix, the first square root filter intermediate vector, the second square root filter intermediate vector, the measurement weight inverse, and the square root filter covariance attenuation coefficient. Specifically, this includes:

[0146] Based on the square root matrix of the total covariance information of the d-th satellite combination, the first square root filter intermediate vector, the second square root filter intermediate vector, and the measurement weight inverse, the state parameter update gain matrix of the d-th satellite combination is determined. The state parameter update gain matrix includes the first state parameter update gain matrix and the second state parameter update gain matrix. The first state parameter update gain matrix is ​​the state parameter update gain matrix of this satellite, and the second state parameter update gain matrix is ​​the state parameter update gain matrix of the linked satellite.

[0147] The updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined based on the square root matrix of the total covariance information, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the square root filter covariance attenuation coefficient, and the state parameter update gain matrix.

[0148] Specifically, the intermediate vector of the single-link square root filter is used. Square root filter covariance attenuation coefficient γ ij The measurement update value of the square root matrix of the autocovariance and mutual covariance information of the state variables of the local satellite and the linked satellite, along with the KALMAN filter gain matrix of the single-link satellite state parameters, can be calculated using the formula:

[0149]

[0150] in, This is the updated covariance information matrix after inter-satellite link measurement update (i.e., the updated value of the square root matrix of total covariance information).

[0151] Specifically, such as based on Figure 2 The principle shown above, the process described above can be summarized as follows: Figure 3 The process shown:

[0152] (1) Each satellite in the constellation carries an inter-satellite link payload. Real-time distance measurement information between satellites can be obtained by using inter-satellite measurements, and satellite orbit determination observation equations based on inter-satellite ranging can be constructed. By modeling the main perturbation forces affecting satellite orbital motion, satellite orbit dynamics equations can be constructed.

[0153] (2) Using the satellite dynamic state parameters such as the prior position, velocity and dynamic model parameters of each satellite, the reference orbit and orbit state transition matrix of each satellite can be calculated by numerical integration method. By using the satellite reference orbit to linearize the satellite orbit determination observation equation and dynamic equation, the linearized autonomous orbit determination observation equation and dynamic state transition equation with the satellite position, velocity and dynamic model parameters as the state parameters to be estimated can be obtained.

[0154] (3) Each satellite can obtain prior values ​​of state parameters such as position, velocity and dynamic model parameters of the satellite to be linked and measured by inter-satellite link communication, as well as the corresponding error covariance information. Combined with the covariance information of the satellite's own state parameters, the square root matrix of the auto-covariance and mutual covariance information of the satellite and the linked satellite can be obtained.

[0155] (4) Using the square root covariance matrix of a pair of linked satellite state parameters and the coefficient matrix of the single-link inter-satellite measurement observation equation, the intermediate vector of the single-link inter-satellite measurement square root filter can be calculated. Using inter-satellite measurement square root filter intermediate vector The measurement noise can be used to calculate the inverse of the inter-satellite measurement weights q for a single link. ij Using the inverse weight q ij The square root filter covariance attenuation coefficient γ can be calculated from the measured noise. ij .

[0156] (5) Using the inverse of the single-link inter-satellite measurement weight q ij The square root matrix of the covariance information of the two corresponding satellites, and the intermediate vector of the square root filter for the corresponding inter-satellite measurement. It can calculate the gain matrix K of the satellite state parameters update for single-link measurement using KALMAN filtering.

[0157] (6) By using the single-link measurement filter to measure and update the gain matrix K, the prior values ​​of satellite state parameters, the measurement residuals, and the coefficient matrix of the observation equation, and using the predicted state of the linked satellite as a reference, the improved satellite state parameters of this satellite can be calculated by the single-link measurement.

[0158] (7) Utilize the single-link measurement filter to measure and update the gain matrix K and the single-link inter-satellite measurement square root filter intermediate vector. Sum of square root filter covariance attenuation coefficient γ ijBy combining the square root matrices of the autocovariance and cross-covariance of the satellite state parameters of the two corresponding satellites, the square root matrices of the autocovariance and cross-covariance information of the satellite and the satellite established in this link can be calculated for the single-link measurement improvement.

[0159] (8) For each satellite, all inter-satellite link observations acquired in a single epoch are processed in steps (4) to (7) sequentially for each measurement link to complete the KALMAN filter measurement update of all observations in a single satellite in a single epoch, and obtain the improved value of the state parameter to be estimated for this satellite, as well as the square root matrix of the self-covariance and mutual covariance information of this satellite and the linked satellite after the filter measurement update.

[0160] (9) Using the dynamic state transition matrix of each satellite and the satellite state parameters updated by measurement, the predicted values ​​of the satellite state parameters after filtering time can be calculated. Using the satellite state transition matrix, satellite dynamic noise information, and the square root matrix of the satellite autocovariance and cross-covariance information updated by filtering measurement, the square root matrix of the covariance information after KALMAN filtering time can be obtained by using the HOUSHOLDER orthogonal transformation method.

[0161] (10) The orbital state parameters of this satellite after time-updated KALMAN filtering, the square root matrix of the self-coherence and mutual coherence variance information between this satellite and the linked satellite, etc., can be sent to other linked satellites through inter-satellite links, and used as input for step (3) to enter the next filtering processing loop. Repeat steps (3)-(10) to complete all filtering processing and realize autonomous orbit determination processing.

[0162] Example 2

[0163] The autonomous orbit determination system for inter-satellite links based on square root filtering in this embodiment includes:

[0164] The satellite combination determination module is used to combine each of the N satellites in the navigation system with other satellites in the navigation system except the current satellite, to obtain N×(N-1) satellite combinations. Each satellite combination includes the current satellite and the link-establishing satellite. The N satellites in the navigation system are all considered as the current satellite, and one current satellite corresponds to N-1 satellite combinations. N>1.

[0165] The preset orbit initialization module is used to initialize the state parameters of the preset orbit of each satellite in the navigation system; the state parameters include position, velocity and solar radiation pressure parameters.

[0166] The first equation construction module is used to construct the inter-satellite ranging observation equation, the first satellite dynamics equation, and the second satellite dynamics equation for each satellite combination. The first satellite dynamics equation is the satellite dynamics equation for this satellite, and the second satellite dynamics equation is the satellite dynamics equation for the linked satellite. The inter-satellite ranging observation equation is a function of the actual state parameters of this satellite, the state parameters of the preset orbit corresponding to this satellite, the actual state parameters of the linked satellite, and the state parameters of the preset orbit corresponding to the linked satellite. The satellite dynamics equation is a function of the actual state parameters of the satellite, Earth's gravity, Sun's gravity, Moon's gravity, solar radiation pressure, and residual perturbation force.

[0167] The second equation construction module is used to construct the satellite state transition equations based on the state parameters of the satellite's preset orbit and the first satellite dynamic equation.

[0168] The orbit determination module is used to perform measurement and time updates on each satellite in the navigation system as if it were the local satellite, thus completing the orbit determination for each satellite.

[0169] Specifically, in determining the orbit of any given satellite by performing measurement and time updates, the orbit determination module is used for:

[0170] Determine whether the current time is the preset final time.

[0171] If so, the final state parameters of the satellite are determined by measuring and updating the correction amount of the N-1 state parameters of the satellite at the current moment and the state parameters of the satellite's preset orbit, thus completing the orbit determination of the satellite.

[0172] If not, at the current time, the Kalman filter method is used to update the measurements of each of the N-1 satellite combinations corresponding to the current satellite, to obtain the current state parameter measurement update correction and the current total covariance information square root matrix update value of the current satellite. Then, the Kalman filter method is used to update the time of each of the N-1 satellite combinations corresponding to the current satellite based on the satellite state transition equation, until the current time is the preset final time.

[0173] Specifically, at any given time, measurements are updated for any N-1 satellite combinations corresponding to the current satellite, including:

[0174] Based on the state parameters of the preset orbits of the satellite and the preset orbits of the linked satellites in the d-th satellite combination, a total covariance information square root matrix is ​​constructed for the d-th satellite combination. The total covariance information square root matrix is ​​a matrix relating to the first self-covariance information square root matrix, the second self-covariance information square root matrix, and the cross-covariance information square root matrix. The first self-covariance information square root matrix is ​​the self-covariance information square root matrix of the state parameters of the satellite, the second self-covariance information square root matrix is ​​the self-covariance information square root matrix of the state parameters of the linked satellite, and the cross-covariance information square root matrix is ​​the cross-covariance information square root matrix of the state parameters of the satellite and the linked satellite. d = 1, 2, ..., N-1.

[0175] Based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit of this satellite, determine the measurement update correction amount of the d-th state parameter of this satellite.

[0176] Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, determine the updated value of the square root matrix of the total covariance information of the d-th satellite combination.

[0177] Determine if d is less than N-1.

[0178] If so, the measurement update correction amount of the d-th state parameter is determined as the state parameter of the preset orbit of the satellite in the (d+1)-th satellite combination, and the d-th satellite combination is replaced with the (d+1)-th satellite combination. The result is "Construct the square root matrix of the total covariance information of the d-th satellite combination based on the state parameters of the preset orbit of the satellite in the d-th satellite combination and the state parameters of the preset orbit of the linked satellite".

[0179] If not, then the updated correction values ​​of the N-1 state parameters of this satellite and the updated value of the square root matrix of the total covariance information of the N-1 satellite combinations corresponding to this satellite are obtained.

[0180] Specifically, for time e, time updates are performed on the N-1 satellite combinations corresponding to each current satellite, where e > 1. This includes:

[0181] Using the satellite state transition equation, the state parameters at time e are determined based on the measurement update corrections of the N-1 state parameters of the current satellite at time e-1 and the state parameters of the current satellite's preset orbit; wherein, the measurement update corrections of the N-1 state parameters of the current satellite at the initial time are the preset corrections.

[0182] Example 3

[0183] An electronic device, comprising:

[0184] One or more processors.

[0185] A storage device on which one or more programs are stored.

[0186] When one or more programs are executed by one or more processors, the one or more processors implement the autonomous orbit determination method for inter-satellite links based on square root filtering as described in Example 1.

[0187] Example 4

[0188] A storage medium storing a computer program, wherein when the computer program is executed by a processor, it implements the autonomous orbit determination method for inter-satellite links based on square root filtering as described in Example 1.

[0189] Technical effects of the present invention:

[0190] (1) When implementing distributed KALMAN filtering data processing, this invention not only utilizes the autocovariance information of the link-established satellite state variables, but also the mutual covariance information. Therefore, the accuracy of the filter results and the stability of the filtering process are improved compared with the equivalent observation distributed filtering method used in traditional inter-satellite link autonomous orbit determination.

[0191] (2) In the KALMAN filter measurement update and time update data processing, the present invention uses the square root of covariance matrix to ensure the positive definiteness of the covariance matrix in the filtering process, and avoids the problem of easy filter divergence caused by the near-rank deficiency of the inter-satellite link observation equation when the SCHMIDT distributed filtering algorithm is directly used.

[0192] (3) This invention uses the square root matrix of covariance information to implement KALMAN filtering for autonomous orbit determination. Under the same processing hardware conditions, this invention has more significant digits in the covariance information. This invention can reduce the impact of computer rounding errors on the accuracy of filtering and improve data processing accuracy. In addition, this invention adopts a link-by-link filtering measurement update method, which helps to make full use of the sparsity of the coefficient matrix of the inter-satellite measurement observation equation to improve the processing algorithm, which facilitates further reduction of data processing computation and improves the computational efficiency of the autonomous orbit determination algorithm.

[0193] (4) The SCHMIDT distributed filtering algorithm based on square root information proposed in this invention can not only be applied to autonomous orbit determination and time synchronization data processing of inter-satellite links, but also be widely applied to various applications of orbit determination, positioning and time synchronization of spacecraft using inter-satellite measurement and communication systems, such as autonomous orbit determination of formation satellites, precise positioning of formation spacecraft, and precise orbit determination of lunar orbit satellites.

[0194] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the systems disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the descriptions are relatively simple; relevant parts can be referred to the method section.

[0195] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A satellite inter-satellite link autonomous orbit determination method based on square root filtering, characterized in that, The method includes: Each of the N satellites in the navigation system is paired with a satellite in the navigation system other than the current satellite to obtain N×(N-1) satellite combinations; each satellite combination includes the current satellite and the link-establishing satellite, and each of the N satellites in the navigation system is considered as a current satellite, with one current satellite corresponding to N-1 satellite combinations; N>1; Initialize the state parameters of the preset orbit of each satellite in the navigation system; the state parameters include position, velocity, and solar radiation pressure parameters; For each satellite combination, construct the inter-satellite ranging observation equation, the first satellite dynamics equation, and the second satellite dynamics equation; the first satellite dynamics equation is the satellite dynamics equation of this satellite, and the second satellite dynamics equation is the satellite dynamics equation of the linked satellite; the inter-satellite ranging observation equation is a function of the actual state parameters of this satellite, the state parameters of the preset orbit corresponding to this satellite, the actual state parameters of the linked satellite, and the state parameters of the preset orbit corresponding to the linked satellite; the satellite dynamics equation is a function of the actual state parameters of the satellite, Earth's gravity, Sun's gravity, Moon's gravity, solar radiation pressure, and residual perturbation force; Based on the state parameters of the satellite's preset orbit and the first satellite dynamic equation, the satellite state transition equation is constructed. Each satellite in the navigation system is treated as the local satellite for measurement and time updates to complete the orbit determination of each satellite; For any given satellite, measurement and time updates are performed to complete the satellite's orbit determination, specifically including: Determine if the current time is the preset final time; If so, the final state parameters of the satellite are determined by measuring and updating the correction amount based on the N-1 state parameters of the satellite at the current moment and the state parameters of the satellite's preset orbit, thus completing the orbit determination of the satellite. If not, at the current time, the Kalman filter method is used to update the measurements of each of the N-1 satellite combinations corresponding to the current satellite, to obtain the current state parameter measurement update correction amount and the current total covariance information square root matrix update value of the current satellite. Then, the Kalman filter method is used to update the time of each of the N-1 satellite combinations corresponding to the current satellite based on the satellite state transition equation, until the current time is the preset final time. Specifically, at any given time, measurements are updated for any N-1 satellite combinations corresponding to the current satellite, including: Based on the state parameters of the preset orbits of the satellite and the preset orbits of the linked satellites in the d-th satellite combination, a total covariance information square root matrix is ​​constructed for the d-th satellite combination. The total covariance information square root matrix is ​​a matrix relating to a first autocovariance information square root matrix, a second autocovariance information square root matrix, and a cross-covariance information square root matrix. The first autocovariance information square root matrix is ​​the autocovariance information square root matrix of the state parameters of the satellite; the second autocovariance information square root matrix is ​​the autocovariance information square root matrix of the state parameters of the linked satellite; and the cross-covariance information square root matrix is ​​the cross-covariance information square root matrix of the state parameters of the satellite and the linked satellite. d = 1, 2, ..., N-1. Based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite and the state parameters of the corresponding preset orbit of this satellite, determine the measurement update correction amount of the d-th state parameter of this satellite; Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, determine the updated value of the square root matrix of the total covariance information of the d-th satellite combination; Determine if d is less than N-1; If so, the measurement update correction amount of the d-th state parameter is determined as the state parameter of the preset orbit of the satellite in the (d+1)-th satellite combination, and the d-th satellite combination is replaced with the (d+1)-th satellite combination. The result is "Construct the square root matrix of the total covariance information of the d-th satellite combination based on the state parameters of the preset orbit of the satellite in the d-th satellite combination and the state parameters of the preset orbit of the linked satellite". If not, then the updated correction values ​​of the N-1 state parameters of this satellite and the updated value of the square root matrix of the total covariance information of the N-1 satellite combinations corresponding to this satellite are obtained. Specifically, for time e, time updates are performed on the N-1 satellite combinations corresponding to each current satellite, where e > 1. This includes: Using the satellite state transition equation, the state parameters at time e are determined based on the measurement update corrections of the N-1 state parameters of the current satellite at time e-1 and the state parameters of the current satellite's preset orbit; wherein, the measurement update corrections of the N-1 state parameters of the current satellite at the initial time are the preset corrections.

2. The autonomous orbit determination method for inter-satellite links based on square root filtering according to claim 1, characterized in that, Based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit of this satellite, the measurement update correction of the d-th state parameter of this satellite is determined, specifically including: Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, the first square root filtering intermediate vector, the second square root filtering intermediate vector, and the measurement weight inverse of the d-th satellite combination are determined; the first square root filtering intermediate vector is the square root filtering intermediate vector of this satellite, and the second square root filtering intermediate vector is the square root filtering intermediate vector of the linked satellite. The measurement update correction amount of the d-th state parameter of this satellite is determined based on the square root matrix of the total covariance information of the d-th satellite combination, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the inverse of the measurement weight, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit.

3. The satellite inter-satellite link autonomous orbit determination method based on square root filtering according to claim 1, characterized in that, Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, the updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined, specifically including: Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, the first square root filter intermediate vector, the second square root filter intermediate vector, the measurement weight inverse, and the square root filter covariance attenuation coefficient of the d-th satellite combination are determined; the first square root filter intermediate vector is the square root filter intermediate vector of this satellite, and the second square root filter intermediate vector is the square root filter intermediate vector of the link-established satellite. The updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined based on the square root matrix of the total covariance information, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the inverse of the measurement weight, and the square root filter covariance attenuation coefficient.

4. The autonomous orbit determination method for inter-satellite links based on square root filtering according to claim 2, characterized in that, Based on the square root matrix of the total covariance information of the d-th satellite combination, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the inverse of the measurement weight, the actual state parameters of this satellite, and the state parameters of the corresponding preset orbit, the measurement update correction amount of the d-th state parameter of this satellite is determined, specifically including: Based on the square root matrix of the total covariance information of the d-th satellite combination, the first square root filter intermediate vector, the second square root filter intermediate vector, and the measurement weight inverse, the state parameter update gain matrix of the d-th satellite combination is determined; the state parameter update gain matrix includes a first state parameter update gain matrix and a second state parameter update gain matrix, the first state parameter update gain matrix is ​​the state parameter update gain matrix of this satellite, and the second state parameter update gain matrix is ​​the state parameter update gain matrix of the linked satellite; The measurement update correction amount of the d-th state parameter of this satellite is determined based on the actual state parameters of this satellite, the state parameters of the corresponding preset orbit, and the state parameter update gain matrix of the d-th satellite combination.

5. The satellite inter-satellite link autonomous orbit determination method based on square root filtering according to claim 3, characterized in that, Based on the square root matrix of the total covariance information of the d-th satellite combination, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the inverse of the measurement weights, and the square root filter covariance attenuation coefficient, the updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined, specifically including: Based on the square root matrix of the total covariance information of the d-th satellite combination, the first square root filter intermediate vector, the second square root filter intermediate vector, and the measurement weight inverse, the state parameter update gain matrix of the d-th satellite combination is determined; the state parameter update gain matrix includes a first state parameter update gain matrix and a second state parameter update gain matrix, the first state parameter update gain matrix is ​​the state parameter update gain matrix of this satellite, and the second state parameter update gain matrix is ​​the state parameter update gain matrix of the linked satellite. The updated value of the square root matrix of the total covariance information of the d-th satellite combination is determined based on the square root matrix of the total covariance information, the intermediate vector of the first square root filter, the intermediate vector of the second square root filter, the square root filter covariance attenuation coefficient, and the state parameter update gain matrix.

6. A satellite inter-satellite link autonomous orbit determination system based on square root filtering, characterized in that, The system includes: The satellite combination determination module is used to combine each of the N satellites in the navigation system with other satellites in the navigation system except the current satellite, to obtain N×(N-1) satellite combinations; each satellite combination includes the current satellite and the link-establishing satellite, and each of the N satellites in the navigation system is considered as a current satellite, with one current satellite corresponding to N-1 satellite combinations; N>1; The preset orbit initialization module is used to initialize the state parameters of the preset orbit of each satellite in the navigation system; the state parameters include position, velocity, and solar radiation pressure parameters; The first equation construction module is used to construct the inter-satellite ranging observation equation, the first satellite dynamics equation, and the second satellite dynamics equation for each satellite combination. The first satellite dynamics equation is the satellite dynamics equation for this satellite, and the second satellite dynamics equation is the satellite dynamics equation for the linked satellite. The inter-satellite ranging observation equation is a function of the actual state parameters of this satellite, the state parameters of the preset orbit corresponding to this satellite, the actual state parameters of the linked satellite, and the state parameters of the preset orbit corresponding to the linked satellite. The satellite dynamics equation is a function of the actual state parameters of the satellite, Earth's gravity, Sun's gravity, Moon's gravity, solar radiation pressure, and residual perturbation force. The second equation construction module is used to construct the satellite state transition equations based on the state parameters of the satellite's preset orbit and the first satellite dynamic equations. The orbit determination module is used to perform measurement updates and time updates on each satellite in the navigation system as if it were the local satellite, thereby completing the orbit determination of each satellite. Specifically, in terms of performing measurement and time updates for any given satellite to complete its orbit determination, the orbit determination module is used for: Determine if the current time is the preset final time; If so, the final state parameters of the satellite are determined by measuring and updating the correction amount based on the N-1 state parameters of the satellite at the current moment and the state parameters of the satellite's preset orbit, thus completing the orbit determination of the satellite. If not, at the current time, the Kalman filter method is used to update the measurements of each of the N-1 satellite combinations corresponding to the current satellite, to obtain the current state parameter measurement update correction amount and the current total covariance information square root matrix update value of the current satellite. Then, the Kalman filter method is used to update the time of each of the N-1 satellite combinations corresponding to the current satellite based on the satellite state transition equation, until the current time is the preset final time. Specifically, at any given time, measurements are updated for any N-1 satellite combinations corresponding to the current satellite, including: Based on the state parameters of the preset orbits of the satellite and the preset orbits of the linked satellites in the d-th satellite combination, a total covariance information square root matrix is ​​constructed for the d-th satellite combination. The total covariance information square root matrix is ​​a matrix relating to a first autocovariance information square root matrix, a second autocovariance information square root matrix, and a cross-covariance information square root matrix. The first autocovariance information square root matrix is ​​the autocovariance information square root matrix of the state parameters of the satellite; the second autocovariance information square root matrix is ​​the autocovariance information square root matrix of the state parameters of the linked satellite; and the cross-covariance information square root matrix is ​​the cross-covariance information square root matrix of the state parameters of the satellite and the linked satellite. d = 1, 2, ..., N-1. Based on the square root matrix of the total covariance information of the d-th satellite combination, the inter-satellite ranging observation equation of the d-th satellite combination, the actual state parameters of this satellite and the state parameters of the corresponding preset orbit of this satellite, determine the measurement update correction amount of the d-th state parameter of this satellite; Based on the square root matrix of the total covariance information of the d-th satellite combination and the inter-satellite ranging observation equation, determine the updated value of the square root matrix of the total covariance information of the d-th satellite combination; Determine if d is less than N-1; If so, the measurement update correction amount of the d-th state parameter is determined as the state parameter of the preset orbit of the satellite in the (d+1)-th satellite combination, and the d-th satellite combination is replaced with the (d+1)-th satellite combination. The result is "Construct the square root matrix of the total covariance information of the d-th satellite combination based on the state parameters of the preset orbit of the satellite in the d-th satellite combination and the state parameters of the preset orbit of the linked satellite". If not, then the updated correction values ​​of the N-1 state parameters of this satellite and the updated value of the square root matrix of the total covariance information of the N-1 satellite combinations corresponding to this satellite are obtained. Specifically, for time e, time updates are performed on the N-1 satellite combinations corresponding to each current satellite, where e > 1. This includes: Using the satellite state transition equation, the state parameters at time e are determined based on the measurement update corrections of the N-1 state parameters of the current satellite at time e-1 and the state parameters of the current satellite's preset orbit; wherein, the measurement update corrections of the N-1 state parameters of the current satellite at the initial time are the preset corrections.

Images