A method for estimating the interstellar arm length of three stars based on capacitive Kalman filtering

Abstract

This invention belongs to the aerospace field and relates to a method for estimating the inter-satellite arm length of three satellites based on capacitive Kalman filtering. The method involves: for the primary satellite and two secondary satellites in a three-satellite gravitational wave detection formation, employing formation-specific nonlinear dynamics modeling, estimating the relative states between the primary and secondary satellites based on capacitive Kalman filtering, and thus obtaining the arm length information between the satellites to achieve three-satellite inter-satellite arm length estimation. This invention effectively suppresses the numerical divergence problem caused by matrix ill-conditioning in traditional Kalman filtering by manipulating the square root of the covariance matrix instead of the original matrix. The capacitive rule achieves third-order accuracy approximation while requiring only a linearly increasing number of sampling points, significantly outperforming computationally expensive methods such as particle filtering. This method is particularly suitable for scenarios with strong nonlinearity and weak observability in satellite dynamics models.

Description

Technical Field

[0001] This invention belongs to the field of aerospace, specifically relating to a method for estimating the inter-star arm length of three stars based on capacitive Kalman filtering. Background Technology

[0002] Space-based gravitational wave detection is a landmark space science project. Currently, all space-based gravitational wave detection programs employ a three-satellite formation to retrieve the subtle effects of gravitational wave signals on spacetime through ultra-long-range laser interferometry ranging. Notable space-based gravitational wave detection programs include the European Space Agency's LISA program, the Chinese Academy of Sciences' Taiji program, and Sun Yat-sen University's Tianqin program. The three-satellite formation forms a giant equilateral triangle formation in space with sides ranging from hundreds of thousands to millions of kilometers. The three satellites are connected by six bidirectional laser links, creating a closed-loop measurement circuit to monitor distance fluctuations between satellites in real time. Each satellite carries an ultra-stable laser interferometer system, capable of achieving sub-picometer level displacement measurement accuracy.

[0003] After being launched in a single-launch, three-satellite configuration, the space gravitational wave detection formation undergoes a complex orbital transfer process to reach its predetermined orbit. During its deep-space cruise phase, it utilizes onboard ion thrusters for several months of flight and performs more than ten orbital corrections. Subsequently, it enters the formation deployment phase, gradually forming a precise triangular configuration through autonomous navigation and control. Finally, a stable six-channel laser link is established, enabling the system to enter scientific observation mode. The entire process places extremely high demands on orbital accuracy, requiring positional errors to be controlled within the hundreds of meters and velocity errors to be no more than centimeters per second. Simulation data from the LISA mission shows that relying solely on dynamic equations for state estimation results in an increase of approximately 1.2 meters in the semi-major axis error per day, leading to an overall positional deviation exceeding 12 meters after ten days. This level of error cannot meet the sub-meter orbital reference requirements of gravitational wave detection. Therefore, effectively utilizing onboard measurement data to achieve autonomous orbital accuracy has become a key research topic. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the existing technology and propose a three-star inter-star arm length estimation method based on capacitive Kalman filtering.

[0005] In view of this, the present invention proposes a method for estimating the interstellar arm length of three stars based on capacitive Kalman filtering, comprising:

[0006] For the primary satellite and two secondary satellites in the three-star gravitational wave detection array, a relative nonlinear dynamic model of the array is adopted. The relative state between the primary and secondary satellites is estimated based on the capacitive Kalman filter, and then the arm length information between the satellites is obtained, so as to realize the arm length estimation between the three stars.

[0007] As an improvement to the above method, the method specifically includes:

[0008] Step 1: Obtain the initial six-element numbers of the primary star and two secondary stars;

[0009] Step 2: Combine the effects of coordinate system rotation, centrifugal force, and differences in gravitational potential gradient to establish nonlinear relative motion equations;

[0010] Step 3: Based on inter-satellite interferometry parameters, solar pointing angle, and solar radial velocity, establish a noise and covariance model, and establish a multi-granularity observation model;

[0011] Step 4: State initialization. Within the set number of loops, repeat the time-series prediction stage and the measurement update stage. Update the process noise and measurement noise by calculating the square root of the error covariance until the number of loops is reached, and obtain the arm length information between satellites to achieve the estimation of the arm length between the three satellites.

[0012] As an improvement to the above method, step 1 includes: the semi-major axis a, eccentricity e, orbital inclination i, and right ascension of the ascending node of the primary star and the two secondary stars. Argument of the pericentric point true anomaly of the orbit f .

[0013] As an improvement to the above method, the nonlinear relative motion equation established in step 2 is:

[0014]

[0015] In the formula, x, y, and z represent the three-dimensional position of the satellite. This represents the rate of change of acceleration of the satellite in the x, y, and z directions. Represents the true anterior angular acceleration. The rate of change of true approximate angular acceleration, This indicates the distance between the primary star SC1 and the central object. This represents the product of the Sun's mass and its gravitational constant. These represent the perturbation forces of the eight celestial bodies in the x, y, and z directions, respectively. These are Gaussian white noise perturbation terms in the x, y, and z directions, respectively, when let for w ( t When ), the following equation is satisfied:

[0016] ;

[0017] The superscript T indicates transpose. This indicates a time delay, where t represents time. Expressing expectations, This represents a 3×3 identity matrix. The variance of white noise is represented. This represents the standard deviation of white noise.

[0018] As an improvement to the above method, the noise and covariance model in step 3 includes:

[0019] All measurements are superimposed with Gaussian white noise, and the covariance matrix is... Let be a diagonal matrix, satisfying the following equation:

[0020]

[0021] in, , , , , These represent the standard deviations of inter-satellite distance, inter-satellite velocity, inter-satellite angle, solar pointing angle, and solar radial velocity, respectively.

[0022] As an improvement to the above method, the multi-granularity observation model in step 3 includes: granularity observation vectors for measuring inter-satellite distances and rates of change. Granular observation vector for inter-satellite angle measurement Granular observation vector of interstellar distances and rates of change superimposed with interstellar angle measurements. , and the granular observation vector of interstellar distances and rates of change and interstellar angles superimposed with solar azimuth measurements. They respectively satisfy the following formulas:

[0023]

[0024]

[0025]

[0026]

[0027] in, These represent the Euclidean distances between any two of the three stars. These represent the relative velocities between any two of the three stars; These represent the azimuth angles between each pair of the three stars. These represent the elevation angles between each pair of the three stars. These represent the equivalent line-of-sight vectors between each pair of the three stars; These represent the angles of the sun's direction between each pair of the three stars. These represent the solar radial velocities of the three stars, with the superscript T indicating transpose.

[0028] As an improvement to the above method, the state initialization in step 4 includes:

[0029] Initial state estimation given satellite position Initial covariance matrix for location uncertainty estimation That is, the covariance matrix Calculated using Cholesky decomposition:

[0030]

[0031] in, It is a lower triangular matrix.

[0032] As an improvement to the above method, the time-series prediction stage of step 4 includes:

[0033] Step S1: Volume point generation:

[0034]

[0035] in, For the set of volume points, for k- The square root factor of covariance at time 1 The weights of the volume point vectors, For prior state estimation;

[0036] Step S2: Nonlinear propagation:

[0037]

[0038] in, The set of volume points after propagation. Represents the nonlinear propagation function;

[0039] Step S3: State Prediction:

[0040]

[0041] in, Indicates by k The prediction at time -1 k Current state m Indicates the total number of volume points. i Indicates the first i One volume point;

[0042] Step S4: Covariance Update:

[0043] in, Indicates by k The prediction at time -1 k The set of volume points representing the state at any given time;

[0044] Step S5: Calculate the square root of the process noise. Combining the results, the square root factor of the predicted covariance at time k is obtained through QR decomposition. :

[0045] .

[0046] As an improvement to the above method, the measurement update phase of step 4 includes:

[0047] Step T1: Observation point generation:

[0048]

[0049] in, For the set of observation volume points, To predict the root covariance factor at time k, The weights of the volume point vectors, It is a posterior state estimate;

[0050] Step T2: Observe the propagation:

[0051]

[0052] in, For the observation volume point after propagation, Represents the observation model;

[0053] Step T3: Observation and Prediction:

[0054]

[0055] in, Indicates by k The prediction at time -1 k Observe constantly. m Indicates the total number of volume points. i Indicates the first i One volume point;

[0056] Step T4: Create a new covariance:

[0057]

[0058] in, Based on k The prediction at time -1 k Predict the average of observed values ​​at any given time; Indicates based on k The prediction at time -1 k Time of the first m Predicted values ​​for each sampling point;

[0059] Step T5: Covariance Correlation:

[0060]

[0061] in, The cross covariance matrix is... This is the state deviation vector. The observation prediction error is represented by T, which stands for transpose.

[0062] Step T6: Calculate the Kalman gain :

[0063]

[0064] in, This indicates that back substitution is used to solve the problem, avoiding direct inversion and making the numerical values ​​more stable.

[0065] Step T7: State Correction:

[0066]

[0067] in, It is the actual sensor measurement value at time k. Let k be the posterior state estimate at time k. Indicates by k The prediction at time -1 k Current state;

[0068] With the square root of the measured noise After combination, QR decomposition is performed to obtain the square root factor used for back-substitution solution.

[0069]

[0070] Step T8: Covariance Correction

[0071]

[0072] in, To predict the observed values, The posterior covariance matrix is... This is a QR decomposition used to calculate the square root of the updated covariance.

[0073] As an improvement to the above method, the inter-satellite arm length estimation includes:

[0074] After completing the state prediction for each satellite position, the inter-satellite arm length is calculated according to the following formula. :

[0075]

[0076] in,( x 1, y 1,z 1) and ( x 2, y 2, z 2) These represent the current positions of the two satellites respectively;

[0077] The arm length estimation error is obtained from the following formula. :

[0078]

[0079] in, This represents the estimated arm length value. This represents the known actual arm length value in the simulation experiment.

[0080] Compared with the prior art, the advantages of the present invention are:

[0081] The method of this invention effectively suppresses the numerical divergence problem caused by matrix ill-conditioning in traditional Kalman filtering by operating on the square root of the covariance matrix instead of the original matrix. The volume rule achieves third-order accuracy approximation while requiring only... The method achieves a linearly increasing number of sampling points, which is significantly better than computationally expensive methods such as particle filtering. This method is particularly suitable for scenarios with strong nonlinearity and weak observability in satellite dynamics models. Attached Figure Description

[0082] Figure 1 This is a flowchart of the square root volume Kalman filter;

[0083] Figure 2 It is the position state estimation error from the x component of satellite 1;

[0084] Figure 3 It is the velocity state estimation error from the x-component of satellite 1;

[0085] Figure 4 This is a diagram showing the estimation error of the arm length between the primary satellite and satellite 1.

[0086] Figure 5 This is a diagram showing the estimation error of the arm length between the primary satellite and satellite 2.

[0087] Figure 6 This is a diagram showing the error in estimating the arm length of satellite 1 and satellite 2. Detailed Implementation

[0088] I. Formation Relative Nonlinear Dynamics Modeling:

[0089] Within a heliocentric inertial frame, an LVLH (Local Vertical Local Horizontal) orbital coordinate system is constructed with the main satellite's center of mass as the origin. This coordinate system is defined as follows: the x-axis is along the instantaneous radial direction from the heliocenter to the main satellite; the z-axis is perpendicular to the orbital plane and in the same direction as the main satellite's angular momentum vector; and the y-axis forms an orthogonal coordinate system according to the right-hand rule. The coordinate system rotation angular velocity vector... ,in The rate of change of the true anomaly angle of the corresponding primary satellite, under the Keplerian orbit assumption, satisfies:

[0090]

[0091] In the formula The gravitational constant of the Sun, , These are the semi-major axis and eccentricity rate of the main satellite orbit, respectively. This is the time-varying radial distance. The main satellite's motion is determined by its orbital root numbers (semi-major axis a, eccentricity e, orbital inclination angle). i Argument of the pericentric point ω Longitude of the ascending node Ω and true anomaly angle f The expression for its position velocity is:

[0092]

[0093]

[0094] The radial velocity component Determined by the characteristics of an elliptical orbit:

[0095]

[0096] Define the relative position vector of the subordinate satellite:

[0097]

[0098] Relative velocity vector:

[0099]

[0100] Construct a 12-dimensional extended state vector:

[0101]

[0102] Considering the effects of coordinate system rotation, centrifugal force, and differences in gravitational potential gradient, a nonlinear equation of relative motion is established:

[0103]

[0104] In the formula The absolute position vector of the subordinate satellite in the LVLH system. This represents the disturbance force of the eight celestial bodies. The Gaussian white noise perturbation term satisfies:

[0105]

[0106] Formation observation system modeling:

[0107] The observation model for the Taiji formation is based on data fusion from multiple types of sensors, converting satellite states (position, velocity) into actual measurements, which is the core input for state estimation. The observation model is divided into two categories, specifically including inter-satellite relative distance measurements, inter-satellite relative velocity measurements, azimuth angle measurements, elevation angle measurements, equivalent line-of-sight vector measurements, and sun pointing measurements. The specific modeling is as follows:

[0108] 1. Inter-satellite measurements (relative status between satellites)

[0109] Relative distance: The Euclidean distance between satellites i and j directly reflects the formation geometry.

[0110]

[0111] Relative velocity: the instantaneous rate of change of distance between satellites, used to capture the dynamic characteristics of formation.

[0112]

[0113] Azimuth:

[0114] Pitch angle:

[0115] Equivalent line-of-sight vector:

[0116] 2. Measurement of Sun Pointing (Satellite's Relative State to the Sun)

[0117] Sun pointing angle: The unit vector pointing from the satellite's center of mass to the sun, used for absolute attitude determination.

[0118]

[0119] Solar radial velocity: The velocity of a satellite along the heliocentric radial direction, measured by Doppler frequency shift.

[0120]

[0121] 3. Noise and Covariance Modeling

[0122] All measurements are superimposed with Gaussian white noise, and the covariance matrix is... For diagonal matrices:

[0123]

[0124] 4. Mathematical Expression of Multi-Granularity Observation Model

[0125] As described above, satellite parameters include inter-satellite interferometry parameters (inter-satellite distance, rate of change of inter-satellite distance, inter-satellite angle), solar pointing angle, and solar radial velocity. Observation vector h The definition is as follows:

[0126]

[0127] In some cases, partial sensor failure or the inability to obtain some measurement parameter information may occur. Therefore, a multi-granularity observation vector is constructed. h These measurements can be categorized into four types: measurements of inter-satellite distances and rates of change only; measurements of inter-satellite angles only (azimuth, elevation, and equivalent line of sight); measurements of inter-satellite distances and rates of change combined with inter-satellite angles; and measurements of inter-satellite distances and rates of change combined with inter-satellite angles combined with solar azimuth. These four granularities are respectively denoted as... h 1. h 2. h 3 and h 4. As shown below

[0128]

[0129]

[0130]

[0131]

[0132] 5. Satellite State Estimation Method Based on Square Root Cumulative Kalman Filter Algorithm

[0133] This study employs the square root occlusive Kalman filter algorithm for satellite state estimation. As a nonlinear filtering method based on numerical stability design, this algorithm is suitable for strongly nonlinear system scenarios such as satellite state estimation. Its core idea is to maintain the square root form of the covariance matrix through Cholesky decomposition and combine it with occlusive rules to complete state prediction and updating. The specific algorithm flow is as follows:

[0134] State initialization: Given an initial state estimate The initial covariance matrix of uncertainty estimation That is, the covariance matrix The square root factor of covariance is calculated using Cholesky decomposition:

[0135]

[0136] in It is a lower triangular matrix.

[0137] Time series forecasting stage:

[0138] 1. Volume point generation: based on the current state dimension generate A symmetrical volume point. Passing through the square root factor. Point set Mapped to state space vector

[0139]

[0140] in, For the set of volume points, for k- The square root factor of covariance at time 1 The weights of the volume point vectors, For prior state estimation;

[0141]

[0142] Represents the identity matrix The i-th column.

[0143] 2. Nonlinear propagation: Inputting each volume point into the system's dynamic model. The predicted point set is obtained.

[0144]

[0145] in, Indicates by k The prediction at time -1 k Current state m Indicates the total number of volume points. i Indicates the first i One volume point;

[0146] 3. State Prediction: Calculate the mean of the predicted states:

[0147]

[0148] 4. Covariance Update: Constructing a Centralized Matrix

[0149]

[0150] and the square root of process noise Combining these elements, we obtain a new square root factor through QR decomposition:

[0151]

[0152] This operation directly incorporates process noise into the covariance update, avoiding explicit calculations. matrix. The matrix represents the process noise covariance matrix.

[0153] Measurement update phase:

[0154] 1) Observation point generation: Reusing the square root factor from the prediction stage Generate a new set of volume points:

[0155]

[0156] 2) Observation propagation: Input each point into the observation model The predicted observation points are obtained as follows:

[0157]

[0158] 3) Observational Prediction: Calculate the predicted observation mean:

[0159]

[0160] 4) Create new covariance: Construct the observation-centered matrix

[0161]

[0162] With the square root of the measured noise QR decomposition after combination

[0163]

[0164] 5) Covariance Correlation: Calculate the cross-covariance between state and observation:

[0165]

[0166] in It is a state-centered matrix, and its calculation method is similar to that of the observation-centered matrix.

[0167] 6) Calculate the Kalman gain: Solve for the Kalman gain through two back-substitution operations:

[0168]

[0169] This method avoids direct inversion and improves numerical stability.

[0170] 7) State correction: using actual observations Updated state estimate:

[0171]

[0172] 8) Covariance Correction: Update the square root factor through QR decomposition:

[0173]

[0174] This step directly incorporates the effect of Kalman gain into the covariance update, preserving the orthogonal triangulation property of the matrix. The algorithm flow is as follows: Figure 1 As shown.

[0175] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and embodiments.

[0176] Example 1

[0177] The following explanation, with reference to the accompanying diagram, uses a three-satellite gravitational wave detection formation with an arm length of 3 million kilometers as an example. For ease of description, the three satellites are named SC1, SC2, and SC3. SC1 is assumed to be the master satellite, and the others are slave satellites. Assuming the state of SC1 is known, we employ relative nonlinear dynamics modeling of the formation's relative motion to estimate the relative states between the master and slave satellites, thereby obtaining the arm length information between them.

[0178] Step 1: Determine the initial six-element number of the satellite orbit

[0179] In the simulation experiment, it is assumed that the initial orbital root numbers of the three satellites are as shown in Table 1:

[0180] Table 1. Initial orbital numbers of the three satellites

[0181]

[0182] Where a represents the semi-major axis, e represents the eccentricity, and i represents the orbital inclination angle. Indicates the right ascension of the ascending node of the orbit. Indicates the argument of the pericentric point. f This represents the true anomaly angle of the orbit, with all angles expressed in radians.

[0183] Step 2: Formation Relative Nonlinear Dynamics Modeling

[0184] Considering the effects of coordinate system rotation, centrifugal force, and differences in gravitational potential gradient, a nonlinear equation of relative motion is established:

[0185]

[0186] In the formula Let be the absolute position vectors of satellites SC2 and SC3 in the LVLH system. This represents the disturbance force of the eight celestial bodies. The Gaussian white noise perturbation term satisfies:

[0187]

[0188] Step 3: Modeling the Formation Observation System

[0189] A formation observation system model was constructed, incorporating inter-satellite interferometry parameters (inter-satellite distance, rate of change of inter-satellite distance, inter-satellite angle), solar pointing angle, and solar radial velocity. Noise and covariance models were established. Based on these, a multi-granularity observation model was constructed.

[0190] Relative distance: The Euclidean distance between satellites i and j directly reflects the formation geometry.

[0191]

[0192] Relative velocity: the instantaneous rate of change of distance between satellites, used to capture the dynamic characteristics of formation.

[0193]

[0194] Azimuth:

[0195] Pitch angle:

[0196] Equivalent line-of-sight vector:

[0197] Sun pointing angle:

[0198]

[0199] Solar radial velocity:

[0200]

[0201] Noise and covariance model:

[0202] All measurements are superimposed with Gaussian white noise, and the covariance matrix is... For diagonal matrices:

[0203]

[0204] Multi-granularity observation model:

[0205] In some cases, partial sensor failure or the inability to obtain some measurement parameter information may occur. Therefore, a multi-granularity observation vector is constructed. h These measurements can be categorized into four types: measurements of inter-satellite distances and rates of change only; measurements of inter-satellite angles only (azimuth, elevation, and equivalent line of sight); measurements of inter-satellite distances and rates of change combined with inter-satellite angles; and measurements of inter-satellite distances and rates of change combined with inter-satellite angles combined with solar azimuth. These four granularities are respectively denoted as... h 1. h 2. h 3 andh 4. As shown below

[0206]

[0207]

[0208]

[0209]

[0210] Assuming the simulation experiment parameters are set as shown in Table 2 in this use case:

[0211] Table 2 Simulation Experiment Parameters

[0212]

[0213] Step 4: Satellite State Estimation Method Based on Square Root Cumulative Kalman Filter Algorithm

[0214] State initialization: Given an initial state estimate The initial covariance matrix of uncertainty estimation The square root factor of covariance is calculated using Cholesky decomposition:

[0215]

[0216] in It is a lower triangular matrix.

[0217] Time series forecasting stage:

[0218] 1. Volume point generation:

[0219]

[0220] 2. Nonlinear propagation:

[0221]

[0222] 3. State prediction:

[0223]

[0224] 4. Covariance Update:

[0225]

[0226] Measurement update phase:

[0227] 1) Observation point generation:

[0228]

[0229] 2) Observation of propagation:

[0230]

[0231] 3) Observation and prediction:

[0232]

[0233] 4) Create a new covariance:

[0234]

[0235] 5) Covariance correlation:

[0236]

[0237] 6) Calculate the Kalman gain:

[0238]

[0239] 7) Status Correction:

[0240]

[0241] 8) Covariance Correction:

[0242]

[0243] Step 5: Assessment of Arm Length Estimation Error

[0244] The arm length estimation error is calculated using the following formula:

[0245]

[0246] in, L error This indicates the error in arm length estimation. This represents the estimated arm length value. L true This represents the actual arm length value.

[0247] After completing the state prediction for each satellite position, calculate according to the following formula. :

[0248]

[0249] in,( x 1, y 1, z 1) and ( x 2, y 2, z 2) These represent the current positions of the two satellites respectively.

[0250] In this use case, a total of [number] simulations were performed. The time interval is seconds, and the result of the relative state estimation between the satellite and the main satellite is as follows: Figure 2 , Figure 3 As shown.

[0251] The histogram of the arm length results between the primary satellite (SC1) and secondary satellite 1 (SC2) is shown below. Figure 4 As shown, the standard deviation of the measurement error is 0.9998m, and the standard deviation of the Kalman filter estimate is 0.0777m.

[0252] The histogram of the arm length results between the primary satellite (SC1) and secondary satellite 2 (SC3) is shown below. Figure 5 As shown, the standard deviation of the measurement error is 1.0002m, and the standard deviation of the Kalman filter estimate is 0.0187m.

[0253] The histogram of the arm length results between satellite 1 (SC2) and satellite 2 (SC3) is shown below. Figure 6 As shown, the standard deviation of the measurement error is 0.9999m, and the standard deviation of the Kalman filter estimate is 0.1116m.

[0254] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to the embodiments, those skilled in the art should understand that modifications or equivalent substitutions to the technical solutions of the present invention do not depart from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A three-satellite inter-satellite arm length estimation method based on cubature Kalman filter, comprising: Step 1: obtaining the initial orbital elements of the primary satellite and two secondary satellites in a three-satellite gravitational wave detection formation; Step 2: combining the effects of coordinate system rotation, centrifugal force terms and differences in gravitational potential gradients, establishing a nonlinear relative motion equation; Step 3: based on the inter-satellite interferometric measurement parameters, the sun pointing angle and the sun radial velocity, establishing a noise and covariance model, and establishing a multi-granularity observation model; Step 4: state initialization, repeating the time series prediction stage and the measurement update stage based on the cubature Kalman filter within a set number of cycles, updating the process noise and measurement noise through the calculated error covariance square root, until the number of cycles is reached, obtaining the arm length information between the satellites, and realizing the three-satellite inter-satellite arm length estimation.

2. The cubature Kalman filter based three-satellite inter-satellite limb length estimation method according to claim 1, characterized in that, The step 1 includes: the semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and true anomaly of the primary satellite and two secondary satellites.

3. The cubature Kalman filter based three-satellite inter-satellite limb length estimation method according to claim 1, characterized in that, The nonlinear relative motion equation established in step 2 is: ； where x, y, z represent three-dimensional positions of the satellite, representing a rate of change of acceleration of the satellite in the x, y, z directions, representing a true anomaly acceleration, a rate of change of true anomaly acceleration, representing a distance of the primary star SC1 from the central celestial body, representing a product of the mass of the sun and the gravitational constant, representing perturbation forces of the eight major celestial bodies in the x, y, z directions, respectively, representing Gaussian white noise perturbation terms in the x, y, z directions, respectively, when = 0, w ( t ) are satisfied. ； where the upper index T denotes the transpose, denotes a time delay, t denotes time, denotes expectation, denotes a 3 x 3 identity matrix, denotes the variance of the white noise, denotes the standard deviation of the white noise.

4. The cubature Kalman filter based three-satellite inter-satellite limb length estimation method according to claim 1, characterized in that, The noise and covariance model of step 3 includes: All measurements are superimposed with Gaussian white noise, covariance matrix is a diagonal matrix satisfying the following equation: ； where, , , , , respectively represent the standard deviation of interstellar distance, interstellar velocity, interstellar angle, solar pointing angle and solar radial velocity.

5. The cubature Kalman filter based three-satellite inter-satellite limb length estimation method according to claim 1, characterized in that, The multi-particle size observation model of step 3 includes: a particle size observation vector only with interstellar distance and rate of change measurement a particle size observation vector only with interstellar angle measurement a particle size observation vector only with interstellar distance and rate of change superimposed interstellar angle measurement a particle size observation vector only with interstellar distance and rate of change and interstellar angle superimposed solar azimuth measurement respectively satisfy the following formula: ； ； ； ； where denote the Euclidean distances between the three stars two by two, denote the relative velocities between the three stars two by two; denote the azimuth angles between the three stars two by two, denote the elevation angles between the three stars two by two, denote the equivalent line-of-sight vectors between the three stars two by two; denote the solar pointing angles between the three stars two by two, denote the solar radial velocities of the three stars, with the upper index T denoting the transpose.

6. The cubature Kalman filter based three-satellite inter-satellite limb length estimation method according to claim 4, characterized in that, The state initialization of step 4 includes: Initial state estimate given satellite position Initial covariance matrix for position uncertainty estimate i.e. covariance matrix computed by Cholesky decomposition: ； wherein is a lower triangular matrix.

7. The cubature Kalman filter based three-satellite inter-satellite limb length estimation method according to claim 1, characterized in that, The time series prediction stage of step 4 includes: Step S1: Cubature point generation: ； wherein, is a set of volume points, is k- is a covariance square root factor at time 1, is a volume point vector weight, is a prior state estimate; Step S2: Nonlinear propagation: ； wherein, is the set of volume points after propagation, denotes a non-linear propagation function; Step S3: State prediction: ； wherein, represents the volume point at the time k predicted at time -1 k state at time m represents the total number of volume points, i represents the volume point at the time i th volume point; Step S4: Covariance update: ； wherein, represents a volume point set of the state at time k predicted at time -1 k predicted at time -1 Step S5: Combining with process noise square root Combining, by QR decomposition, the k-time prediction covariance square root factor : 。 8. The cubature Kalman filter based three-satellite inter-satellite limb length estimation method according to claim 7, characterized in that, The measurement update stage of step 4 includes: Step T1: Observation point generation: ； wherein, is a set of observation volume points, is a k-time predicted covariance square root factor, is a volume point vector weight, is a posterior state estimate; Step T2: Observation propagation: ； wherein, is the propagated observation volume point, denotes the observation model; Step T3: Observation prediction: ； wherein, represents the volume point at the time of observation, k predicted at the time of -1, k observed at the time of 0, m represents the total number of volume points, i represents the volume point at the time of observation, i the volume point at the time of observation. Step T4: Create new covariance: ； wherein, is based on k the prediction obtained at time k the predicted observation mean at time denotes the prediction based on k the prediction obtained at time k the predicted value at the m th sampling point at time With the measurement noise square root Combining and QR decomposition to get square root factors for back substitution solution : ； Step T5: Covariance association: ； wherein is the cross-covariance matrix, is the state bias vector, is the observation prediction error, T denotes transpose; Step T6: Calculate Kalman gain : ； Step T7: State correction: ； wherein, is the actual sensor measurement at time k, is the posterior state estimate at time k, denotes the prediction from k -1 time prediction k time state; Step T8: Covariance correction: ； where, is the predicted observation, is the posterior covariance matrix, is the QR decomposition used to compute the square root of the updated covariance.

9. The cubature Kalman filter based three-satellite inter-satellite limb length estimation method according to claim 1, characterized in that, The inter-satellite arm length estimation includes: After the state prediction for each satellite position is completed, the satellite inter-arm length is calculated according to the following formula : ； wherein, x 1, y 1, z 1) and x 2, y 2, z 2) represent the current positions of the two satellites, respectively; The arm length estimation error is obtained according to the following formula : ； wherein, represents the known true arm length value in the simulation experiment.

Images