CN107421550B - An Earth-Lagrange Joint Constellation Autonomous Orbit Determination Method Based on Inter-satellite Ranging - Google Patents

Abstract

The invention discloses an earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging, which comprises the following steps of: the method comprises the following steps: establishing a state equation of an earth-Lagrange combined constellation autonomous orbit determination system; step two: establishing a measurement equation of an earth-Lagrange combined constellation autonomous orbit determination system; step three: determining a filtering method for realizing orbit parameter estimation; step four: the earth-Lagrange combined constellation autonomous orbit determination method based on the selected filtering algorithm is specifically realized. According to the invention, by introducing Lagrange satellites, the problem of 'deficit rank' existing when autonomous orbit determination is carried out only by using inter-satellite ranging information is effectively solved, and the complexity of system equipment is reduced; the statistical characteristic of the system noise is estimated on line in real time through the self-adaptive nonlinear filtering algorithm, the requirement on noise prior information is low, the stability of the autonomous orbit determination filtering algorithm is improved, and the autonomous orbit determination precision is improved.

Description

An Earth-Lagrange Joint Constellation Autonomous Orbit Determination Method Based on Inter-satellite Ranging

technical field

The invention belongs to the field of autonomous orbit determination of earth satellite constellations, in particular, it is a method for realizing autonomous determination of earth satellite constellations based on adaptive filtering only by using inter-satellite ranging information in an earth-Lagrange joint constellation. rail method.

Background technique

In recent decades, satellite navigation has played an increasingly important role in the national economy and military struggle. Improving the autonomous operation of satellites is important for reducing the burden of ground measurement and control, reducing satellite operating costs, improving satellite survivability and expanding the application of spacecraft. potential, etc. are of great significance. The satellite's autonomous orbit determination technology is the premise of realizing the autonomous control of the satellite, which plays a vital role in ensuring its autonomous operation, and has also become the development trend of today's satellite navigation and control technology.

At present, the autonomous orbit determination methods of the constellation are mainly divided into two types: the independent orbit determination of each satellite and the overall autonomous orbit determination of the constellation. The former uses various inertial devices, star sensors and navigation receivers on the satellite to obtain information such as the distance and direction of the satellite relative to other natural celestial bodies or navigation stars to estimate its own position online. The system has a simple structure and requires less computation. It is small, but its orbit determination accuracy is limited by factors such as the irregularity of the surface of natural celestial bodies and the measurement accuracy of sensors. It is difficult to achieve the accuracy of navigation-level applications when used alone; the latter makes full use of the relative measurement between satellites in the constellation. The orbit determination of the constellation is based on information to realize the whole network orbit determination of the constellation. The orbit determination of the short arc segment has high accuracy, but there is a problem of "deficient rank", that is, only the relative measurement information between the satellites cannot determine the overall direction of the constellation in the inertial space. The orbit determination accuracy of the long arc segment decreases gradually with the increase of time. In order to solve the problem of "deficient rank", some foreign scholars have pointed out that when the gravitational field in which the constellation or satellite formation is located is highly asymmetric, absolute information will be introduced from the "uniqueness" of the orbit, thereby solving the unobservable problem of the overall rotation of the constellation. The study found that among many asymmetric perturbation gravitational fields, the relative strength of the third body gravitational force is the largest, especially in the area near the two Lagrange points L ₁ and L ₂ , so joint autonomy is carried out through the joint Earth satellite constellation and the Lagrange point constellation. Orbit determination can solve the "deficient rank" problem. Using this method of autonomous orbit determination can greatly reduce the complexity of the spacecraft's autonomous navigation system, reduce the number of navigation and measurement equipment, and improve the accuracy of autonomous orbit determination. At the same time, the Lagrange navigation constellation can also be used as a deep space exploration navigation base station to provide navigation information for other deep space probes, thereby effectively solving the problems of the ground deep space exploration network in deep space navigation accuracy, real-time, security, etc. Therefore, It is of great significance to study the autonomous orbit determination of the combined Earth-Lagrange constellation.

At present, the autonomous orbit determination methods for the Earth-Lagrange joint constellation at home and abroad mainly use conventional algorithms such as batch algorithm, Extended Kalman Filter (EKF) algorithm and Unscented Kalman Filter (UKF) algorithm. as the basis. However, batch processing algorithms are often used for post-processing and are not suitable for real-time orbit determination. Although conventional EKF and UKF algorithms use a recursive form that is convenient for real-time calculation, they often require prior information of system noise and measurement noise to be known. Otherwise, the system performance will be degraded and even filter divergence will be caused. Due to the complexity of the actual force of the satellites in orbit, especially the Lagrange satellites, it is difficult to obtain more accurate noise statistical characteristics, and irregular modeling errors caused by the maneuvers of satellite orbit maintenance will also cause changes in statistical characteristics, making the priori The information is meaningless, so the research on the autonomous navigation algorithm that can be adaptively adjusted according to the measurement information sequence is of great significance to improve the robustness of the autonomous orbit determination system and ensure the accuracy of the autonomous orbit determination.

SUMMARY OF THE INVENTION

In order to solve the above-mentioned problems, the present invention proposes a method for autonomous orbit determination of the Earth-Lagrange joint constellation based on inter-satellite ranging. The Adaptive Nonlinear Filter, ANF) algorithm adopts a centralized structure to fuse the effective measurement information to realize the calculation of the position and velocity of each satellite in the Earth-Lagrange joint constellation. The constellation orbit determination method first adopts a centralized structure to establish the system model of the Earth-Lagrange joint constellation autonomous orbit determination, in which the state equation of the earth navigation satellite is mainly established based on the received two-body problem dynamics model in the geocentric inertial rectangular coordinate system. , the state equation of the Lagrange navigation satellite is mainly established based on the dynamic model of the circular restricted three-body problem in the center of mass rendezvous coordinate system, and the measurement equation adopts the inter-satellite distance measurement model in the geocentric inertial rectangular coordinate system; The nonlinear filtering algorithm fuses the effective measurement information to realize the autonomous orbit determination of the Earth-Lagrange joint constellation. Finally, a joint constellation composed of three medium earth orbit (MEO) satellites and two Lagrange satellites is used as an example to verify the effectiveness of the method.

An Earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging, which specifically includes the following steps:

Step 1: Under the centralized structure, establish the state equation of the Earth-Lagrange joint constellation autonomous orbit determination system;

Step 2: Under the centralized structure, establish the measurement equation of the Earth-Lagrange joint constellation autonomous orbit determination system;

Step 3: Under the centralized structure, according to the above established model of the Earth-Lagrange joint constellation autonomous orbit determination system, determine a filtering method for realizing orbit parameter estimation;

Step 4: Concrete realization of the Earth-Lagrange joint constellation autonomous orbit determination method based on the selected filtering algorithm under the centralized structure.

The advantages of the present invention are:

(1) By introducing Lagrange satellites to provide absolute information for the system, the problem of "deficient rank" is effectively solved, and the long-term autonomous orbit determination accuracy of the Earth-Lagrange joint constellation is improved;

(2) Only the inter-satellite ranging information is used to reduce the complexity of the system equipment;

(3) The Lagrange navigation constellation has a high coverage rate of the earth navigation constellation, which improves the accuracy of autonomous orbit determination;

(4) Using the adaptive nonlinear filtering algorithm to estimate the statistical characteristics of the system noise online in real time, the accuracy and stability of the filtering algorithm for autonomous orbit determination are improved, thereby improving the accuracy of autonomous orbit determination.

Description of drawings

Figure 1 is a schematic diagram of the composition of the Earth-Lagrange joint constellation autonomous orbit determination system based on inter-satellite ranging;

Fig. 2 is the structure diagram of the method for autonomous orbit determination of the Earth-Lagrange joint constellation based on inter-satellite ranging;

3 is a flowchart of an autonomous navigation algorithm based on adaptive nonlinear filtering;

Fig. 4 is the variation trend diagram of the system accuracy of the embodiment with b _Q , b _C , and b _q ;

FIG. 5 is a position error graph of an embodiment GNSS satellite.

Detailed ways

The present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

The present invention is an autonomous orbit determination method based on inter-satellite ranging based on the Earth-Lagrange joint constellation, and its system structure is shown in Figure 1. The figure shows that the joint constellation orbit determination system consists of the earth navigation constellation and its inter-satellite link ranging, The Lagrange navigation constellation and its inter-satellite link ranging, as well as the inter-satellite link ranging between the two constellations of satellites, consist of the earth navigation constellation consisting of m medium/high-orbit satellites, numbered 1,2,...,i, ...m; the Lagrange Navigation Constellation consists of n satellites operating in Halo periodic orbits near L ₁ and L ₂ , numbered 1, 2, ..., k, ... n in order. Figure 2 shows the basic implementation process of an Earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging, which shows that the basic implementation process of the orbit determination method mainly includes measurement information acquisition and constellation orbit parameter estimation. The acquired measurement information includes three types of distance measurement information: the inter-satellite distance measurement information within the Earth Navigation Constellation, the inter-satellite distance measurement information within the Lagrange Navigation Constellation, and the distance measurement information between the Earth Navigation Satellite and the Lagrange Navigation Satellite; In parameter estimation, ① the state equation of the joint constellation orbit determination system is established based on the simplified orbital dynamics model of the Earth satellite and the Lagrange satellite, ② the measurement equation of the three types of distance measurement information is established, and ③ the system model composed of the state equation and the measurement equation is used. , using a suitable adaptive nonlinear filtering algorithm to estimate the orbit parameters; finally, output the estimated position and velocity of the satellites in the constellation. Under the geocentric inertial coordinate system and the center of mass meeting coordinate system, the Cartesian coordinate description method is used (that is, the current operating state of the satellite is represented by the position and velocity components), and the adaptive nonlinear filtering estimation method is used. The system model established with the measurement information realizes the autonomous orbit determination of the Earth-Lagrange joint constellation. The specific steps are:

Step 1: Under the centralized structure, establish the state equation of the Earth-Lagrange joint constellation autonomous orbit determination system;

An accurate system model is one of the main factors to ensure the accuracy of constellation operating state parameter estimation. The present invention establishes the system state equation based on the satellite orbit dynamics model.

(1) Earth Navigation Satellite Orbit Dynamics Model

Earth navigation satellites are mainly medium and high-orbit satellites, and the gravitational force of the earth's center of mass is far greater than other forces. Generally, the dynamic model of the received two-body problem is adopted. In order to meet the requirements of computational efficiency and high precision, the main perturbation term is selected according to the influence of the perturbation force. For medium- and high-orbit satellites, among many perturbation terms, the perturbation effect caused by the non-spherical symmetry of the earth can reach the order of 10 ⁴ meters, of which the perturbation effect of the gravitational field of order 3 and above is much smaller than that of the J ₂ term. , can be neglected; the influence of sun-moon gravitational perturbation can reach 10 ³ meters; the influence of solar light pressure perturbation can reach 10 ² meters ^; Below ¹ meter level, it can be ignored; at the same time, the atmospheric density at high orbit is very small, and the influence of atmospheric resistance perturbation can be ignored. Therefore, the present invention mainly considers the orbital dynamics model composed of the earth's centroid gravitational force and the J ₂ term and three main perturbing forces of sun and moon gravitational force to establish the system state equation.

In the earth-centered rectangular inertial coordinate system, combined with the orbital dynamics model, the state equation for orbit determination of the earth satellite i to be measured in the joint constellation is established as follows:

这里[x_Ei,y_Ei,z_Ei]^T和

分别为该地球卫星的位置矢量和速度矢量。F_E为系统状态函数，W_Ei(t)为系统噪声，满足均值为q_Ei(t)，方差为Q_Ei(t)的高斯白噪声。式(1)可以进一步详细写为：In the formula, the state vector corresponding to the earth satellite i to be measured is

where [x _Ei , y _Ei , z _Ei ] ^T and

are the position vector and velocity vector of the earth satellite, respectively. F _E is the system state function, W _Ei (t) is the system noise, satisfying the Gaussian white noise with mean value q _Ei (t) and variance Q _Ei (t). Equation (1) can be further written as:

和

为x、y、z三轴速度噪声分量，

和

是x、y、z三轴加速度等效噪声分量，包括未建模摄动项和噪声项；

和

分别为地球卫星i对应的已建模速度矢量和加速度矢量，其中加速度主要考虑地球中心引力加速度a_0,Ei、J₂摄动项加速度

和日月引力摄动加速度a_MS,Ei，具体表达式如下：in,

and

is the three-axis velocity noise component of x, y, and z,

and

is the equivalent noise component of the x, y, and z three-axis acceleration, including the unmodeled perturbation term and noise term;

and

are the modeled velocity vector and acceleration vector corresponding to the earth satellite i respectively, in which the acceleration mainly considers the acceleration of the earth's center gravitational acceleration a _{0, Ei} , J ₂ perturbation term

and the sun-moon gravitational perturbation acceleration a _MS,Ei , the specific expression is as follows:

1) the earth center gravitational acceleration a _{0 of the satellite i considered by the present invention, Ei} satisfies:

为待测卫星i的地心距，μ为地球质量M_E与引力常数G的乘积，即地球引力常数。In the formula,

is the geocentric distance of the satellite _i to be measured, and μ is the product of the earth's mass ME and the gravitational constant G, that is, the earth's gravitational constant.

满足：2) J ₂ perturbation term acceleration

Satisfy:

In the formula, _Re is the equatorial radius of the earth, and J ₂ is the second-order harmonic coefficient.

3) The perturbation acceleration a _MS,Ei caused by the gravity of the sun and the moon satisfies:

和

分别为卫星i到太阳和月球的距离，

和

分别为地心到太阳和月球的距离，μ_S和μ_M分别为太阳引力常数和月球引力常数。where (x _S , y _S , z _S ) and (x _M , y _M , z _M ) are the three-dimensional position coordinates of the sun and the moon in the earth-centered rectangular inertial coordinate system, respectively,

and

are the distances from satellite i to the sun and the moon, respectively,

and

are the distances from the center of the earth to the sun and the moon, respectively, and μ _S and μ _M are the solar gravitational constant and the lunar gravitational constant, respectively.

(2) Lagrange Navigation Satellite Orbit Dynamics Model

Different from the earth navigation satellite, the Lagrange satellite operates near the Lagrange point of the Earth-Moon system, and the magnitude of the central gravitational action of the Earth and the central gravitational action of the moon is similar, so when describing its motion, the present invention adopts a circular restricted three-body model.

In the center-of-mass rendezvous coordinate system, combined with the orbital dynamics model, the equation of state for orbit determination of the lunar satellite k to be measured in the Lagrange constellation can be established as follows:

这里[x_Lk,y_Lk,z_Lk]^T和

分别为该星的无量纲位置矢量和速度矢量，其特征质量、特征长度和特征时间如式(7)所示，其中M_E和M_M分别为地球质量和月球质量，L为地月距离。F_L为系统状态函数，W_Lk(t)为系统噪声，满足均值为q_Lk(t)，方差为Q_Lk(t)的高斯白噪声。In the formula, the state vector corresponding to the Lagrange lunar satellite k to be tested is

where [x _Lk , y _Lk , z _Lk ] ^T and

are the dimensionless position vector and velocity vector of the star, respectively, and its characteristic mass, characteristic length and characteristic time are shown in Equation (7), where M _E and M _M are the mass of the earth and the moon, respectively, and L is the distance between the earth and the moon. _FL is the system state function, and W _Lk (t) is the system noise, satisfying the Gaussian white noise with mean value q _Lk (t) and variance Q _Lk (t).

Equation (6) can be further written as:

为Lagrange卫星k对应的已建模速度矢量，μ₀＝M_M/(M_M+M_E)，Δ_E,Lk、Δ_M,Lk分别为Lagrange卫星k的地心距和月心距，

和

为Lagrange卫星k的x、y、z三轴速度噪声分量，

和

是Lagrange卫星k的x、y、z三轴加速度等效噪声分量，包括未建模摄动项和噪声项。in

is the modeled velocity vector corresponding to the Lagrange satellite k, μ ₀ =M _M /(M _M +M _E ), Δ _E,Lk , Δ _M,Lk are the geocentric and lunar distances of the Lagrange satellite k, respectively,

and

is the x, y, and z three-axis velocity noise components of the Lagrange satellite k,

and

is the equivalent noise component of the x, y, and z three-axis acceleration of the Lagrange satellite k, including the unmodeled perturbation term and the noise term.

(3) State equation of the Earth-Lagrange joint constellation autonomous orbit determination system

It is assumed that the entire autonomous orbit determination system includes m Earth satellites and n Lagrange satellites. According to the state equations (1) and (6) of the above single satellite orbit determination, the state equation of the constellation orbit determination with a centralized structure is established as follows:

which is:

为状态向量，F为整个星座定轨系统的状态函数矢量，W(t)表示均值为q(t)，方差为Q(t)的系统噪声。In the formula,

is the state vector, F is the state function vector of the entire constellation orbit determination system, W(t) represents the system noise with mean q(t) and variance Q(t).

Step 2: Under the centralized structure, establish the measurement equation of the Earth-Lagrange joint constellation autonomous orbit determination system;

The observation information used by the Earth-Lagrange joint constellation autonomous orbit determination system is only the inter-satellite distance measurement, but since the dynamic models of the two types of satellites are established in different coordinate systems, the inter-satellite distance measurement models of the entire constellation are divided into Three categories:

(1) Inter-satellite distance measurement model between earth navigation satellites

For earth navigation satellites, denote the observed distance between satellite i and satellite j as ρ _Ei,Ej , then we have:

ρ _Ei,Ej =h _E (X _Ei ,X _Ej )+v _Ei,Ej , i≠j (11)

v_Ei,Ej为卫星i和卫星j之间的距离等效测量噪声，包括建模误差和量测噪声。in

v _Ei,Ej is the distance equivalent measurement noise between satellite i and satellite j, including modeling error and measurement noise.

(2) Inter-satellite distance measurement model between Lagrange navigation satellites

Denote the observed distance between Lagrange satellites i and j as ρ _Li,Lj , we have:

ρ _Li,Lj =h _L (X _Li ,X _Lj )+v _Li,Lj , i≠j (12)

v_Li,Lj为Lagrange卫星i和j之间的距离等效测量噪声，包括建模误差和量测噪声。in

v _{Li, Lj} is the distance equivalent measurement noise between Lagrange satellites i and j, including modeling error and measurement noise.

(3) Measurement model between earth navigation satellites and Lagrange satellites

In order to represent the distance relationship between the GNSS satellite i and the Lagrange satellite k, the position parameters of the two satellites need to be represented in the same coordinate system, as follows:

是Lagrange卫星在地心惯性坐标系中的位置向量，它是关于Lagrange卫星在质心会合坐标系中的状态向量X_Lk中的无量纲位置向量[x_Lk,y_Lk,z_Lk]^T的函数，其转换关系为：in,

is the position vector of the Lagrange satellite in the geocentric inertial coordinate system, which is a function of the dimensionless position vector [x _Lk , y _Lk , z _Lk ] ^T in the state vector X _Lk of the Lagrange satellite in the center of mass rendezvous coordinate system, Its conversion relationship is:

Among them, R represents the rotation matrix from the convergent coordinate system of the center of mass to the inertial coordinate system of the earth's center, and the expression is as follows:

where u _M = ω _M + θ _M , i _M , Ω _M , ω _M and θ _M represent the orbital inclination of the moon's orbit around the Earth, the right ascension of the ascending node, the argument of perigee and the true perigee angle, respectively.

(4) Measurement equation of the Earth-Lagrange joint constellation autonomous orbit determination system

According to the above measurement equations (11), (12) and (13), the constellation orbit determination measurement equation with a centralized structure is established as follows:

Z(t)=h[X(t),t]+V(t) (16)

which is:

where V(t) is zero mean and the variance is the equivalent measurement noise of R(t). Equation (17) gives the observation equation in general. In practical application, some distance measurements will be affected by factors such as the occlusion of the earth or the moon, resulting in unmeasurable situations. ), the corresponding ranging information will be eliminated.

Step 3: According to the above-established Earth-Lagrange joint constellation autonomous orbit determination system model based on inter-satellite distance measurement, determine a filtering method for realizing orbit parameter estimation;

The state of the nonlinear system needs to be estimated using nonlinear filtering methods. Conventional nonlinear filtering methods, such as EKF, UKF, volumetric filtering and particle filtering, can achieve ideal results when the prior statistical information of noise is accurately known. However, in this system, the equivalent system noise _Wk not only contains the unmodeled error, but also contains the discretization error (and linearization error) in the filtering algorithm, which makes the noise have time-varying non-Gaussian statistical characteristics, and it is difficult to Accurately obtain. For this reason, when applying these nonlinear filtering algorithms, the noise is usually approximated as a stationary random white noise sequence, that is, Q _k ≡ Q, and Q is determined empirically. In order to ensure that the state estimation error can converge to a stable value during the entire filtering process, the value of Q is often too large, that is, to satisfy _Q≥max {Qk}, which makes it difficult to achieve optimal system performance.

In order to solve this problem, the present invention, based on the conventional nonlinear filtering algorithm, combines the idea of Sage-Husa adaptive algorithm to estimate and correct the statistical characteristics of the system noise online, so as to improve the filtering accuracy and stability. Considering the practicality of the project, the present invention adopts the adaptive EKF method or the adaptive Sigma Points Kalman Filter (SPKF) method to realize the estimation of satellite orbit parameters.

First, when the system has high requirements on the operation speed or the system storage space is small, the adaptive EKF method is used to estimate the state of the discretized nonlinear system. The basic idea of the adaptive EKF is to use the Taylor series expansion method to discretize and linearize the continuous nonlinear state equation and observation equation, and then use the time-varying noise statistical estimator on the basis of the classical Kalman filter. , estimate and correct the statistical characteristics of system noise and observation noise in real time, and then perform state estimation. Specific steps are as follows:

(1) Time update

Compute the one-step predicted state:

为k时刻的系统噪声均值的估计值，

为利用动力学方程进行的状态一步递推，一种方法是通过离散化和线性化后的方程(19)计算得到，其中

为k时刻的状态估计，

为(k+1)时刻的一步预测状态估计，Δt为积分步长；另一种方法是通过4阶Runge-Kutta(RK4)方法直接对连续状态方程(9)数值积分计算得到。In the formula,

is the estimated value of the mean value of the system noise at time k,

For the one-step recurrence of states using the kinetic equation, one method is computed by discretizing and linearizing equation (19), where

is the state estimate at time k,

is the one-step predicted state estimation at (k+1) time, and Δt is the integration step size; another method is to directly numerically integrate the continuous state equation (9) by the fourth-order Runge-Kutta (RK4) method.

Compute the state transition matrix:

In the formula, Φ _k+1/k is the state transition matrix from time k to time (k+1).

Compute the one-step forecast error covariance matrix:

为k时刻系统噪声方差阵的估计值。In the formula, P _k+1/k is the error covariance matrix of the one-step prediction state, P _k is the error covariance matrix of the estimated state at time k,

is the estimated value of the system noise variance matrix at time k.

Considering that the system state estimation error has not converged in the initial period of filtering and there are few samples used to estimate the system noise, the estimated noise statistical characteristics fluctuate greatly and are not accurate enough at this time, which will make the filtering converge slowly, and even lead to severe cases. Filter divergence. Therefore, during the initial period of filtering (0~T), this algorithm does not apply the estimated noise statistical characteristics to the time update process, but uses the prior statistical information of the noise. The time update algorithm at this time is as follows:

(2) Measurement update:

Calculate the gain matrix:

为基于一步预测状态估计

计算的观测向量h的Jacobian矩阵，R_k+1为(k+1)时刻测量噪声协方差阵。In the formula, K _k+1 is the gain matrix at time (k+1),

is a state estimation based on one-step prediction

The Jacobian matrix of the calculated observation vector h, R _k+1 is the measurement noise covariance matrix at the time (k+1).

Compute the innovation vector:

In the formula, υ _k+1 is the innovation vector at time (k+1), and Z _k+1 is the measurement vector at time (k+1).

Compute the state estimate:

Compute the state estimation error covariance matrix:

Compute an estimate of the noise statistics of the system:

in,

和

的值分别趋向于(1-b_q)、(1-b_Q)和(1-b_C)，因此b_q、b_Q和b_C越大，当前信息的比重越小。

侧重于跟踪系统噪声均值q的变化，取值过大会使

偏小，取值过小则会引起滤波发散；

侧重于平滑

的取值，取值过小平滑效果较差，取值过大则使

偏小，增大了滤波发散的风险；

用来跟踪方差C的变化，取值过小会使估值的波动变大，取值过大则会使跟踪产生延迟。In the formula, 0<b _q <1, 0<b _Q <1, and 0<b _C <1 are the forgetting factors for calculating q, Q, and C, respectively. As k increases,

and

The values of are tending to (1-b _q ), (1-b _Q ) and (1-b _C ) respectively, so the larger b _q , b _Q and b _C are, the smaller the weight of current information is.

Focus on tracking the change of the system noise mean q, if the value is too large, it will make the

If the value is too small, it will cause filter divergence;

focus on smoothness

The value of , if the value is too small, the smoothing effect is poor, and if the value is too large, the

It is too small, which increases the risk of filter divergence;

It is used to track the change of variance C. If the value is too small, the fluctuation of the valuation will increase, and if the value is too large, the tracking will be delayed.

负定的情况，为了保证滤波的稳定，在负定的情况下对式(29)中的

项的对角线元素取绝对值，非对角线元素取0来近似处理。In practical applications, although the calculations of equations (29) and (30) have a certain smoothing effect, there will still be

In the case of negative definite, in order to ensure the stability of filtering, in the case of negative definite

The diagonal elements of the item take the absolute value, and the off-diagonal elements take 0 for approximation.

Second, when the system pays more attention to the accuracy index, and there is no special restriction on the operation speed and occupied storage space, the adaptive Sigma Points Kalman Filter (SPKF) method is used to estimate the state of the continuous nonlinear system. The basic idea of adaptive SPKF is to directly transfer the Sigma sampling points obtained according to the deterministic sampling strategy through a nonlinear function, and then, on the basis of the classical Kalman filter algorithm framework, through the time-varying noise statistical estimator, real-time estimation and correction Statistical characteristics of system noise and observation noise, and then state estimation. Specific steps are as follows:

1) Select the Sigma point sampling strategy

At present, the commonly used Sigma point sampling strategies include symmetric sampling, simplex sampling, central difference sampling and volume sampling. Select a sampling strategy and calculate the weight coefficients W _i ^m and W _i ^c of SPKF, where _i =0,1,...,L, L is the number of sampling points determined by the specific sampling strategy, and Wi ^m is used for the mean value estimate, W _ic is ^used for variance estimation.

2) Time update:

和估计误差协方差P_k计算k时刻的Sigma点集χ_i,k(i＝0,1,…,L)。According to the sampling strategy selected in 1) and the state estimate value at time k

and the estimated error covariance P _k to calculate the Sigma point set χ _i,k (i=0,1,...,L) at time k.

Calculate the one-step predicted value of the Sigma point:

χ _i,k+1/k =F[χ _i,k ,k], i=0,1,...,L (34)

Among them, χ _i,k is the value of the i-th Sigma point at time k, and χ _i,k+1/k is the one-step prediction value of the i-th Sigma point at time k+1.

Compute the one-step predicted state:

为k时刻的系统噪声均值的估计值，

为(k+1)时刻的一步预测状态估计。In the formula,

is the estimated value of the mean value of the system noise at time k,

is the one-step predicted state estimate at time (k+1).

Compute the state one-step prediction error covariance matrix:

为k时刻系统噪声方差阵的估计值；In the formula, P _k+1/k is the error covariance matrix of the one-step prediction state,

is the estimated value of the system noise variance matrix at time k;

Within the initial time 0 to T of the filter, the prior statistical information of the noise is used, and the time update is as follows:

3) Measurement update:

和P_k+1/k更新Sigma点集；According to the sampling strategy selected in 1) and

Update the Sigma point set with P _k+1/k ;

Estimated measurements:

γ _{i, k+1/k} = h(χ _{i, k+1/k} ) (39)

为第(k+1)时刻测量向量的加权估计值。where γ _i,k+1/k is the estimated value of the measurement vector at the (k+1)th time calculated by the i-th Sigma point,

is the weighted estimate of the measurement vector at the (k+1)th time.

Calculate the gain matrix:

Among them, R _k+1 is the measurement noise covariance matrix at the time (k+1), P _y,k+1/k is the covariance of the estimated value of the measurement vector at the (k+1)th time, P _{xy,k+1/ k} is the covariance between the estimated value of the measurement vector and the state one-step predicted value at the (k+1)th time, and K _k+1 is the gain matrix at the (k+1) time.

Compute the innovation vector:

where υ _k+1 is the innovation vector at (k+1) time, and Z _k+1 is the measurement vector at (k+1) time;

Calculate the state estimate at time (k+1)

Calculate the state estimation error covariance matrix P _k+1 at time (k+1):

4) Calculate the estimated value of the statistical characteristics of the system noise:

in,

In the formula, 0<b _q <1, 0<b _Q <1, and 0<b _C <1 are the forgetting factors for calculating q, Q, and C, respectively.

负定的情况时，对式(48)中的

项的对角线元素取绝对值，非对角线元素取0。when it appears

In the case of negative definiteness, for Eq. (48)

The diagonal elements of the item take the absolute value, and the off-diagonal elements take 0.

Step 4: Concrete realization of the Earth-Lagrange joint constellation autonomous orbit determination algorithm based on the selected filtering method under the centralized structure;

The flowchart of realizing the autonomous orbit determination algorithm of the Earth-Lagrange joint constellation according to the adaptive nonlinear filtering algorithm proposed by the present invention is shown in Figure 4, and the specific process is as follows:

The specific process is as follows:

(1) Data initialization;

初始误差协方差阵P₀、初始系统噪声均值

初始系统噪声方差阵

测量噪声方差阵R以及遗忘因子b_q、b_Q、b_C；The initialized parameters include: k=0, filter step size h, noise switching parameter T, initial state

Initial error covariance matrix P ₀ , initial system noise mean

Initial system noise variance matrix

Measure the noise variance matrix R and forgetting factors b _q , b _Q , b _C ;

初始误差协方差阵P₀、初始系统噪声方差阵

根据需要确定；初始系统噪声均值

取零；测量噪声方差阵R根据测量设备的性能确定；遗忘因子b_q、b_Q取值介于0.99～1.0之间，b_C取值介于0.9～0.99之间；Among them, the noise switching parameter T is the time node for switching the statistical characteristics of noise in the time update algorithm: when k<T, use the prior statistical characteristics of noise for time update; when k≥T, use the estimated statistical characteristics of noise to perform time update update; initial state

Initial error covariance matrix P ₀ , initial system noise variance matrix

Determined as needed; initial system noise mean

Take zero; the measurement noise variance matrix R is determined according to the performance of the measurement equipment; the forgetting factors b _q and b _Q are between 0.99 and 1.0, and the value of b _C is between 0.9 and 0.99;

(2) If the adaptive EKF algorithm is adopted, go directly to (4); if the adaptive SPKF algorithm is adopted, go to (3);

(3) Calculate the weight coefficients W _i ^m and W _i ^c according to the selected Sigma point sampling strategy;

(4) If k<T, use the priori statistical characteristics of noise to perform time update, and go to (5); otherwise, use the estimated noise statistical characteristics to perform time update, and go to (6);

(5) Time update using a priori statistical characteristics of noise; go to (7) after completion;

及误差协方差阵P_k作为k时刻的状态初值及误差协方差阵，依次计算状态一步预测值

及预测误差协方差矩阵P_k+1/k,实现自适应Kalman滤波估计方法的时间更新。当采用自适应EKF算法时，计算公式分别为式(22)和式(23)；当采用自适应SPKF算法时，计算公式分别为式(37)和(38)；Using state estimation at time k

and the error covariance matrix P _k as the initial state value and the error covariance matrix at time k, and calculate the state one-step prediction value in turn

and the prediction error covariance matrix P _k+1/k , to realize the time update of the adaptive Kalman filter estimation method. When the adaptive EKF algorithm is adopted, the calculation formulas are formula (22) and formula (23) respectively; when the adaptive SPKF algorithm is adopted, the calculation formulas are formulas (37) and (38) respectively;

(6) Using the time update of the estimated noise statistics;

及误差协方差阵P_k作为k时刻的状态初值及误差协方差阵，依次计算状态一步预测值

及预测误差协方差矩阵P_k+1/k,实现自适应Kalman滤波估计方法的时间更新。当采用自适应EKF算法时，计算公式分别为式(18)和式(21)；当采用自适应SPKF算法时，计算公式分别为式(35)和(36)；Using state estimation at time k

and the error covariance matrix P _k as the initial state value and the error covariance matrix at time k, and calculate the state one-step prediction value in turn

and the prediction error covariance matrix P _k+1/k , to realize the time update of the adaptive Kalman filter estimation method. When the adaptive EKF algorithm is used, the calculation formulas are formula (18) and formula (21) respectively; when the adaptive SPKF algorithm is used, the calculation formulas are formulas (35) and (36) respectively;

(7) Measurement update;

及误差协方差矩阵P_k+1，实现自适应Kalman滤波的测量更新。当采用自适应EKF算法时，计算公式分别为式(24)～(27)；当采用自适应SPKF算法时，计算公式分别为式(43)～(46)；Calculate the gain matrix K _k+1 , the innovation matrix υ _k+1 , and the state estimation value at (k+1) time in turn

and the error covariance matrix P _k+1 to realize the measurement update of the adaptive Kalman filter. When the adaptive EKF algorithm is adopted, the calculation formulas are respectively equations (24) to (27); when the adaptive SPKF algorithm is adopted, the calculation formulas are respectively equations (43) to (46);

(8) Estimate the statistical characteristics of system noise;

和协方差估计值

当采用自适应EKF算法时，计算公式分别为式(28)和(29)；当采用自适应SPKF算法时，计算公式分别为式(47)和(48)；Calculate the mean estimate of the system noise in turn

and covariance estimates

When the adaptive EKF algorithm is used, the calculation formulas are formulas (28) and (29) respectively; when the adaptive SPKF algorithm is used, the calculation formulas are formulas (47) and (48) respectively;

(9) Determine whether the filtering process is over; if it is necessary to continue filtering, k=k+1, and return to (4).

The present invention firstly studies the problem of "deficient rank" when only using inter-satellite distance measurement, and introduces Lagrange satellites to provide an absolute reference; Autonomous track.

A method for autonomous orbit determination of the Earth-Lagrange joint constellation based on inter-satellite ranging proposed by the present invention firstly establishes the state equation of the Earth-Lagrange joint constellation under a centralized structure; then establishes the measurement equation of the Earth-Lagrange joint constellation ; Finally, an adaptive nonlinear filtering method is used to realize the autonomous orbit determination of the Earth-Lagrange joint.

Example:

Take an Earth-Lagrange constellation as an example, the constellation consists of 3 MEO satellites and 2 Lagrange satellites. The MEO satellites are numbered E ₁ to E ₃ ; the Lagrange satellites are numbered L ₁ to L ₂ . In the simulation, it is found that the orbit determination accuracy of each satellite in the MEO constellation is similar, and the orbit determination accuracy of each satellite in the Lagrange constellation is also similar. Therefore, only the MEO satellite numbered E ₁ and the Lagrange satellite numbered L ₁ are used as examples in this paper. illustrate.

Simulation conditions:

STK is used to simulate the actual operation of the constellation, and to generate the nominal orbit data during the operation of the constellation.

(1) Nominal data generation parameter setting

The simulation time is from 12:00:00 on January 1, 2016 to 12:00:00 on June 28, 2016, which lasted 180d, and the step size was T=1s. The orbital dynamics model of the earth satellite adopts the high-precision orbit propagator (HPOP), in which the earth model adopts the WGS84-EGM96 model, which considers 21×21 order aspheric perturbations; other perturbation terms include solar gravity, Lunar gravity, solar light pressure (light pressure coefficient C _r =1.0) and atmospheric drag (atmospheric drag coefficient C _d =2.2); the ratio of satellite cross-sectional area to satellite mass S/m=0.02m ² /kg. The orbital dynamics model of the Lagrange satellite adopts a circular restricted three-body problem model.

(2) Basic filter parameter settings

The state equation of the earth satellite also adopts the received two-body problem model, but the perturbation term only considers the J2 term and the gravitational force of the sun and the moon; the state equation of the Lagrange satellite adopts the circular restricted three-body problem model, and the filter period T=100s.

σ_Lx0＝σ_Ly0＝σ_Lz0＝10^-10[L]，

测量噪声方差阵为

σ_Ei,Ej＝σ_Ei,Lk＝σ_Lk,Ll＝5m；初始状态估计误差协方差阵P₀＝diag([P_E1,0 P_E2,0 P_E3,0 P_L1,0 P_L2,0])，Under the centralized structure, the initial value of the system noise variance matrix for autonomous orbit determination of the Earth-Lagrange joint constellation is Q ₀ =diag([Q _E1,0 Q _E2,0 Q _E3,0 Q _L1,0 Q _L2,0 ]) ,in

σ _Lx0 =σ _Ly0 =σ _Lz0 =10 ⁻¹⁰ [L],

The measurement noise variance matrix is

σ _Ei,Ej =σ _Ei,Lk =σ _Lk,Ll =5m; initial state estimation error covariance matrix P ₀ =diag([P _E1,0 P _E2,0 P _E3,0 P _L1,0 P _L2,0 ]),

σ_Ex,0＝σ_Ey,0＝σ_Ez,0＝1m，

σ_Lx,0＝σ_Ly,0＝σ_Lz,0＝1m，

in

σ _Ex,0 =σ _Ey,0 =σ _Ez,0 =1m,

σ _Lx,0 =σ _Ly,0 =σ _Lz,0 =1m,

(3) Evaluation method

Because the system state and estimation error in the actual system are both random quantities, after obtaining the estimation error at each simulation time, a statistical method (calculating the root mean square error of the error data) is used to evaluate the absolute orbit determination and relative relative orbit determination of the satellites in the constellation. positioning effect. The estimated state after the convergence of the filter estimation is used as a sample to calculate the error. The formula for calculating the root mean square error is as follows:

为第i时刻的估计值，X_i为第i时刻的真值，n为数据总数。in,

is the estimated value at the i-th moment, X _i is the true value at the i-th moment, and n is the total number of data.

In order to facilitate quantitative analysis, the sum of squares of three-axis position errors and the sum of squares of three-axis velocity errors are used as the distance error and rate error of satellite orbit determination. The calculation formula of the root mean square error is as follows:

Satellite distance estimation error:

Satellite rate estimation error:

Implementation: The state vector of the centralized Earth-Lagrange joint constellation autonomous orbit determination system model is composed of the positions and speeds of five stars, and the inter-satellite distance measurement method is used, focusing on the improvement of the orbit determination accuracy of the adaptive EKF.

Step 1: Establish the state equation of the Earth-Lagrange joint constellation autonomous orbit determination system

其中[x_Ei,y_Ei,z_Ei]^T和

分别为该星的位置矢量和速度矢量，结合轨道动力学模型式(1)可建立卫星Ei定轨的状态方程：In the earth-centered rectangular inertial coordinate system, the state vector of the earth navigation satellite i is

where [x _Ei , y _Ei , z _Ei ] ^T and

are the position vector and velocity vector of the satellite, respectively, combined with the orbital dynamics model formula (1), the state equation of the satellite Ei orbit determination can be established:

In the formula, the acceleration caused by the earth's center gravity, J ₂ perturbation force and sun-moon gravity of satellite Ei are a _0,Ei , a _J2,Ei and a _MS,Ei , respectively. For the specific definitions, please refer to equations (3) to (5) . W _Ei (t) represents the system noise of satellite Ei, including unmodeled perturbation, satisfying E[W _Ei (t)]=0, E[W _Ei (t)W _Ei (τ)]=Q _Ei (t) δ(t-τ), δ(t-τ) is the Dirac function, that is

其中[x_Lk,y_Lk,z_Lk]^T和

分别为该星的无量纲位置矢量和速度矢量，结合轨道动力学模型式(6)可建立卫星L_k定轨的状态方程：The state vector corresponding to the Lagrange satellite k is

where [x _Lk , y _Lk , z _Lk ] ^T and

are the dimensionless position vector and velocity vector of the star, respectively. Combined with the orbital dynamics model formula (6), the state equation for the orbit determination of the satellite L _k can be established:

In the formula, W _Lk (t) represents the system noise of satellite Ei, including unmodeled perturbation, satisfying E[W _Lk (t)]=0, E[W _Lk (t)W _Lk (τ)]=Q _Lk (t)δ(t-τ).

According to the state equation of the above single satellite, the state equation of the five-star constellation with a centralized structure can be established as follows:

为包含星座中五颗卫星所有状态的状态向量，W(t)表示系统噪声，满足E[W(t)]＝q(t)，E[W(t)W(τ)]＝Q(t)δ(t-τ)。In the formula,

is the state vector containing all states of the five satellites in the constellation, W(t) represents the system noise, which satisfies E[W(t)]=q(t), E[W(t)W(τ)]=Q(t )δ(t-τ).

Step 2: Establish the measurement equation of the Earth-Lagrange joint constellation autonomous orbit determination system

Under the centralized structure, the measurement equation of the five-star constellation using the inter-satellite link measurement is as follows:

In the formula, h _E is the measurement equation between the three Earth satellites, h _EL is the measurement equation between the Earth satellite and the Lagrange satellite, and h _L is the measurement equation between the two Lagrange satellites, see equations (11) to (13) .

V(t)=[ _ΔEi,Ej … _ΔEi,Lk … _ΔLk,Ll …] ^T is the residual pseudorange observation noise vector after error compensation such as ionospheric and tropospheric delays, with zero mean and standard deviation R (t) white noise.

Step 3: Under the centralized structure, according to the above established model of the Earth-Lagrange joint constellation autonomous orbit determination system, determine a filtering method for realizing orbit parameter estimation.

In this example, the algorithm requires fast calculation speed and less computing resource requirements, so the adaptive EKF algorithm is adopted.

Step 4: Under the centralized structure, the adaptive EKF estimation method is used to realize the joint autonomous orbit determination of the Earth-Lagrange five-star constellation, and the performance of the joint autonomous orbit determination of the Earth-Lagrange five-star constellation is simulated and analyzed.

Combined with the Earth-Lagrange five-star constellation centralized joint autonomous orbit determination flowchart (Fig. 4) and the fourth step of the invention, the adaptive EKF estimation method is used to realize the Earth-Lagrange five-star constellation joint autonomous orbit determination.

初始误差协方差阵P₀、初始系统噪声均值

初始系统噪声方差阵

测量噪声方差阵R＝5m以及遗忘因子b_q＝0.9999、b_Q＝0.9999、b_C＝0.90。(1) Data initialization. The initialized parameters include: k=0, filter step size h=100s, threshold T=(10 days×86400s/day)/100s, initial state

Initial error covariance matrix P ₀ , initial system noise mean

Initial system noise variance matrix

Measurement noise variance matrix R=5m and forgetting factors b _q =0.9999, b _Q =0.9999, b _C =0.90.

Note: Because the value of b _q is too small, it will cause filter divergence, so a larger value is taken as the initial value here, first adjust b _Q and b _C to the best value, and finally adjust b _q .

(2) If k<T, use the prior statistical characteristics of noise to perform time update, and go to (3); otherwise, use the estimated noise statistical characteristics to perform time update, and go to (4).

及误差协方差阵P_k作为k时刻的状态初值及误差协方差阵，根据式(22)和式(23)分别计算状态一步预测值

及预测误差协方差矩阵P_k+1/k,实现自适应EKF估计方法的时间更新。转到(5)。(3) Temporal update using a priori statistical properties of noise. Using state estimation at time k

and the error covariance matrix P _k as the initial state value and the error covariance matrix at time k, and calculate the state one-step predicted value according to Equation (22) and Equation (23) respectively

and the prediction error covariance matrix P _k+1/k to realize the time update of the adaptive EKF estimation method. Go to (5).

及误差协方差阵P_k作为k时刻的状态初值及误差协方差阵，根据式(18)和式(21)分别计算状态一步预测值

及预测误差协方差矩阵P_k+1/k,实现自适应EKF估计方法的时间更新。(4) Temporal update with estimated noise statistics. Using state estimation at time k

and the error covariance matrix P _k as the initial state value and the error covariance matrix at time k, and calculate the state one-step predicted value according to Equation (18) and Equation (21) respectively

and the prediction error covariance matrix P _k+1/k to realize the time update of the adaptive EKF estimation method.

及误差协方差矩阵P_k+1，实现EKF的测量更新。(5) Adaptive EKF measurement update. Calculate the gain matrix K _k+1 , the innovation matrix υ _k+1 , and the state estimation value at the time (k+1) sequentially according to equations (24) to (27).

and the error covariance matrix P _k+1 to realize the measurement update of the EKF.

和协方差估计值

(6) Estimate the statistical characteristics of the system noise. Calculate the mean estimated value of system noise according to equations (28) and (29) respectively

and covariance estimates

(7) Determine whether the filtering process is over. To continue filtering, k=k+1, and return to (2).

Step 5: Adjust forgetting factors b _Q , b _C , b _q .

Combined with the functions of the three forgetting factors in step 3 of the invention and the approximate range of the forgetting factors given in step 4, adjust b _Q , b _C , b _q in turn, and repeat the filtering process in step 3, so that the system reaches the relative maximum. Excellent, the specific adjustment steps are as follows:

(1) Fix b _Q = 0.9999, b _q = 0.9999, increase b _C gradually, and determine the optimal value of b _C according to the filtering result, denoted as

b_q＝0.9999，逐渐减小b_Q，根据滤波结果确定b_Q的最佳取值，记为

(2) Fixed

b _q =0.9999, gradually reduce b _Q , determine the optimal value of b _Q according to the filtering result, denoted as

逐渐减小b_q，根据滤波结果确定b_q的最佳取值，记为

(3) Fixed

Decrease b _q gradually, and determine the optimal value of b _q according to the filtering result, denoted as

Figure 4 shows the variation trend of system accuracy with b _Q , b _C , and b _q . According to the figure, it can be determined that the optimal forgetting factor combination is: b _q =0.9993, b _Q =0.9997, and b _C =0.96.

On the basis of analyzing the "deficient rank" problem when only using inter-satellite distance measurement, the present invention introduces Lagrange satellites to provide an absolute reference; then, under the centralized structure, an adaptive EKF estimation method is adopted to realize the Earth-Lagrange joint The autonomous orbit of the constellation.

When the inter-satellite ranging error is zero mean and the standard deviation is 5m, the following simulation results are obtained:

Since the orbit determination accuracy of the three MEO satellites in the constellation is similar, taking satellite E1 as an example, the orbit determination distance errors of the x, y and z axes are 7.78m, 7.60m and 10.06m respectively. The error curves are shown in Figure 5. The distance estimation error is 14.81m; the orbit determination speed errors of the x, y, and z axes are 1.43×10 ^-3 m/s, 3.92×10 ^-3 m/s and 4.74×10 ^-3 m/s, respectively. The estimated error is 6.31×10 ^-3 m/s. Similarly, the orbit determination accuracy of the two Lagranges is similar. Taking satellite L1 as an example, the orbit determination distance errors of the x, y and z axes are 2.66m, 3.60m and 3.39m respectively, and the distance estimation error is 5.62m; The orbit determination velocity errors of the three axes and z are 2.11×10 ^-5 m/s, 1.64×10 ^-5 m/s and 1.48×10 ^-5 m/s respectively, and the velocity estimation error is 3.06×10 ^-5 m/s.

The effectiveness of the method of the present invention is verified.

Claims (4)

1. An Earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging, comprising the following steps: Step 1: Under the centralized structure, establish the state equation of the Earth-Lagrange joint constellation autonomous orbit determination system; Specifically: (1) Establish a dynamic model of the earth navigation satellite orbit In the earth-centered rectangular inertial coordinate system, the state equation for orbit determination of the earth satellite i to be measured in the joint constellation is as follows:

In the formula, the state vector corresponding to the earth satellite i to be measured is

[x _Ei ,y _Ei ,z _Ei ] ^T and

are the position vector and velocity vector of the earth satellite, respectively; F _E is the system state function, W _Ei (t) is the system noise, satisfying the Gaussian white noise with mean value q _Ei (t) and variance Q _Ei (t); It can be written in further detail as:

in,

and

is the three-axis velocity noise component of x, y, and z,

and

is the equivalent noise component of the x, y, and z three-axis acceleration, including the unmodeled perturbation term and noise term;

and

are the modeled velocity vector and acceleration vector corresponding to the earth satellite i respectively, in which the acceleration mainly considers the acceleration of the earth's center gravitational acceleration a _{0, Ei} , J ₂ perturbation term

and the sun-moon gravitational perturbation acceleration a _MS,Ei , the specific expression is as follows: 1) The earth center gravitational acceleration a _0,Ei of satellite i satisfies:

In the formula,

is the geocentric distance of the satellite i to be measured, and μ is the product of the mass of the earth and the gravitational constant G, that is, the gravitational constant of the earth; 2) J ₂ perturbation term acceleration

Satisfy:

where _Re is the equatorial radius of the earth, and J ₂ is the second-order harmonic coefficient; The perturbation acceleration a _MS,Ei caused by the gravity of the sun and the moon satisfies:

where (x _S , y _S , z _S ) and (x _M , y _M , z _M ) are the three-dimensional position coordinates of the sun and the moon in the earth-centered rectangular inertial coordinate system, respectively,

and

are the distances from satellite i to the sun and the moon, respectively,

and

are the distances from the center of the earth to the sun and the moon, respectively, μ _S and μ _M are the solar gravitational constant and the lunar gravitational constant, respectively; (2) Lagrange Navigation Satellite Orbit Dynamics Model The state equation for orbit determination of the lunar satellite k to be measured in the Lagrange constellation is established in the coordinate system of the center of mass convergence:

In the formula, the state vector corresponding to the Lagrange lunar satellite k to be tested is

[x _Lk , y _Lk , z _Lk ] ^T and

are the dimensionless position vector and velocity vector of the star, respectively, and its characteristic mass, characteristic length and characteristic time are shown in formula (7), where M _E and M _M are the mass of the earth and the moon, respectively, and L is the distance between the earth and the moon; _FL is the system state function, W _Lk (t) is the system noise, satisfying the Gaussian white noise with mean value q _Lk (t) and variance Q _Lk (t);

Formula (6) is written as:

in,

is the modeled velocity vector corresponding to the Lagrange satellite k, μ ₀ =M _M /(M _M +M _E ), Δ _E , _Lk , Δ _M,Lk are the satellite’s geocentric and lunar distances, respectively,

and

is the three-axis velocity noise component of x, y, and z,

and

is the equivalent noise component of the x, y, and z three-axis acceleration, including the unmodeled perturbation term and noise term; (3) State equation of the Earth-Lagrange joint constellation autonomous orbit determination system Assuming that the entire autonomous navigation system includes m earth satellites and n Lagrange satellites, according to the state equation of the orbit determination of the earth satellite i to be measured and the state equation of the orbit determination of the lunar satellite k to be measured, a centralized structure of the constellation orbit determination state is established The equation is as follows:

which is:

In the formula,

is the state vector, F is the state function vector of the entire constellation orbit determination system, W(t) represents the system noise with mean value q(t) and variance Q(t); Step 2: Establish the measurement equation of the Earth-Lagrange joint constellation autonomous orbit determination system; Step 3: Under the centralized structure, according to the above established model of the Earth-Lagrange joint constellation autonomous orbit determination system based on inter-satellite distance measurement, determine a filtering method for realizing orbit parameter estimation; Step 4: Under the centralized structure, the Earth-Lagrange joint constellation based on the selected filtering algorithm autonomously determines the orbit. 2. a kind of Earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging according to claim 1, described step 2, is specially: (1) Inter-satellite distance measurement model between earth navigation satellites For earth navigation satellites, denote the observed distance between satellite i and satellite j as ρ _Ei,Ej , then we have: ρ _Ei,Ej =h _E (X _Ei ,X _Ej )+v _Ei,Ej , i≠j (11) in

v _{Ei, Ej} is the equivalent measurement noise between satellite i and satellite j, including modeling error and measurement noise; (2) Inter-satellite distance measurement model between Lagrange navigation satellites Denote the observed distance between Lagrange satellites i and j as ρ _Li,Lj , we have: ρ _Li,Lj =h _L (X _Li ,X _Lj )+v _Li,Lj , i≠j (12) in

v _{Li, Lj} is the equivalent measurement noise between Lagrange satellites i and j, including modeling error and measurement noise; (3) Measurement model between earth navigation satellites and Lagrange satellites Representing GNSS satellite i and Lagrange satellite k in the same coordinate system, there are:

in,

is the coordinate of the Lagrange satellite in the geocentric inertial coordinate system, which is a function of the dimensionless coordinate X _Lk of the Lagrange satellite in the center of mass convergence coordinate system, and its conversion relationship is:

Among them, R represents the rotation matrix from the convergent coordinate system of the center of mass to the inertial coordinate system of the earth's center, and the expression is as follows:

where u _M = ω _M + θ _M , i _M , Ω _M , ω _M and θ _M represent the orbital inclination of the moon's orbit around the Earth, the right ascension of the ascending node, the argument of perigee and the true perigee angle, respectively; (4) Measurement equation of the Earth-Lagrange joint constellation autonomous orbit determination system According to the above measurement equations (11), (12) and (13), the constellation orbit determination measurement equation with a centralized structure is established as follows: Z(t)=h[X(t),t]+V(t) (16) which is:

where V(t) is the equivalent measurement noise with zero mean and variance R(t). 3. a kind of Earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging according to claim 1, described step 3, adopts EKF method or SPKF method, concrete: (1) When the adaptive EKF method is used to estimate the state of the discretized nonlinear system, the specific steps are as follows: 1) Time update Compute the one-step predicted state:

In the formula,

is the estimated value of the mean value of the system noise at time k,

For the one-step recurrence of states using the kinetic equation, one method is computed by discretizing and linearizing equation (19), where

is the state estimate at time k,

is the one-step predicted state estimation at time (k+1), Δt is the integration step size; another method is to directly integrate the continuous state equation (9) numerically by the fourth-order Runge-Kutta method;

Compute the state transition matrix:

In the formula, Φ _k+1/k is the state transition matrix from time k to time (k+1); Compute the one-step forecast error covariance matrix:

In the formula, P _k+1/k is the error covariance matrix of the one-step prediction state, P _k is the error covariance matrix of the estimated state at time k,

is the estimated value of the system noise variance matrix at time k; Within the initial time 0 to T of the filter, the prior statistical information of the noise is used, and the time update is as follows:

2) Measurement update: Calculate the gain matrix:

In the formula, K _k+1 is the gain matrix at time (k+1),

is a state estimation based on one-step prediction

The Jacobian matrix of the calculated observation vector h, R _k+1 is the measurement noise covariance matrix at time (k+1); Compute the innovation vector:

where υ _k+1 is the innovation vector at (k+1) time, and Z _k+1 is the measurement vector at (k+1) time; Compute the state estimate:

Compute the state estimation error covariance matrix:

3) Calculate the estimated value of the statistical characteristics of the system noise:

in,

In the formula, 0<b _q <1, 0<b _Q <1, 0<b _C <1 are the forgetting factors for calculating q, Q, and C, respectively. As k increases,

and

The values tend to be (1-b _q ), (1-b _Q ) and (1-b _C ), respectively; when it appears

In the case of negative definiteness, for Eq. (29)

The diagonal elements of the item take the absolute value, and the off-diagonal elements take 0; (2) When the adaptive SPKF method is used to estimate the state of the continuous nonlinear system, the specific steps are as follows: 1) Select the Sigma point sampling strategy Calculate the weight coefficients W _i ^m and W _i ^c of the Kalman filtering of Sigma points, where _i =0,1,...,L, L is the number of sampling points, Wi ^m is used for mean estimation _, and Wi ^c is used for variance estimation; 2) Time update: According to the sampling strategy selected in 1) and the state estimate value at time k

Calculate the Sigma point set χ _{i,k at time k with the estimated error covariance P k ;} Calculate the one-step predicted value of the Sigma point: χ _i,k+1/k =F[χ _i,k ,k], i=0,1,...,L (34) Among them, χ _i,k is the value of the i-th Sigma point at time k, and χ _i,k+1/k is the one-step prediction value of the i-th Sigma point at time k+1; Compute the one-step predicted state:

In the formula,

is the estimated value of the mean value of the system noise at time k,

is the one-step prediction state estimate at time (k+1); Compute the state one-step prediction error covariance matrix:

In the formula, P _k+1/k is the error covariance matrix of the one-step prediction state,

is the estimated value of the system noise variance matrix at time k; Within the initial time 0 to T of the filter, the prior statistical information of the noise is used, and the time update is as follows:

3) Measurement update: According to the sampling strategy selected in 1) and

Update the Sigma point set with P _k+1/k ; Estimated measurements: γ _i,k+1/k =h(χ _i,k+1/k ) (39)

Where: γ _i,k+1/k is the estimated value of the measurement vector at the (k+1)th time calculated by the i-th Sigma point,

is the weighted estimate of the measurement vector at the (k+1)th moment; Calculate the gain matrix:

Among them, R _k+1 is the measurement noise covariance matrix at the time (k+1), P _y,k+1/k is the covariance of the estimated value of the measurement vector at the (k+1)th time, P _{xy,k+1/ k} is the covariance of the estimated value of the measurement vector and the state one-step predicted value at the (k+1)th time, and K _k+1 is the gain matrix at the (k+1) time; Compute the innovation vector:

where υ _k+1 is the innovation vector at (k+1) time, and Z _k+1 is the measurement vector at (k+1) time; Calculate the state estimate at time (k+1)

Calculate the state estimation error covariance matrix P _k+1 at time (k+1):

4) Calculate the estimated value of the statistical characteristics of the system noise:

in,

In the formula, 0<b _q <1, 0<b _Q <1, 0<b _C <1 are the forgetting factors for calculating q, Q, and C, respectively; when it appears

In the case of negative definiteness, for Eq. (48)

The diagonal elements of the item take the absolute value, and the off-diagonal elements take 0. 4. a kind of Earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging according to claim 1, described step 4, is specially: (1) Data initialization; The initialized parameters include: k=0, filter step size h, noise switching parameter T, initial state

Initial error covariance matrix P ₀ , initial system noise mean

Initial system noise variance matrix

Measure the noise variance matrix R and forgetting factors b _q , b _Q , b _C ; Among them, the noise switching parameter T is the time node for switching the statistical characteristics of noise in the time update algorithm: when k<T, use the prior statistical characteristics of noise for time update; when k≥T, use the estimated statistical characteristics of noise to perform time update update; initial state

Initial error covariance matrix P ₀ , initial system noise variance matrix

Determined as needed; initial system noise mean

Take zero; the measurement noise variance matrix R is determined according to the performance of the measurement equipment; the forgetting factors b _q and b _Q are between 0.99 and 1.0, and the value of b _C is between 0.9 and 0.99; (2) If the adaptive EKF algorithm is adopted, go directly to (4); if the adaptive SPKF algorithm is adopted, go to (3); (3) Calculate the weight coefficients W _i ^m and W _i ^c according to the selected Sigma point sampling strategy; (4) If k<T, use the prior statistical characteristics of noise to perform time update, and go to (5); otherwise, use the estimated noise statistical characteristics to perform time update, and go to (6); (5) Time update using a priori statistical characteristics of noise; go to (7) after completion; Using state estimation at time k

and the error covariance matrix P _k as the initial state value and the error covariance matrix at time k, and calculate the state one-step prediction value in turn

and the prediction error covariance matrix P _k+1/k to realize the time update of the adaptive Kalman filter estimation method; when the adaptive EKF algorithm is used, the calculation formulas are equation (22) and equation (23) respectively; When the Sigma point Kalman filtering algorithm is used, the calculation formulas are equations (37) and (38) respectively; (6) Using the time update of the estimated noise statistics; Using state estimation at time k

and the error covariance matrix P _k as the initial state value and the error covariance matrix at time k, and calculate the state one-step prediction value in turn

and the prediction error covariance matrix P _k+1/k , to realize the time update of the adaptive Kalman filter estimation method; when the adaptive EKF algorithm is adopted, the calculation formulas are equation (18) and equation (21) respectively; In the SPKF algorithm, the calculation formulas are equations (35) and (36) respectively; (7) Measurement update; Calculate the gain matrix K _k+1 , the innovation matrix υ _k+1 , and the state estimation value at (k+1) time in turn

and the error covariance matrix P _k+1 to realize the measurement update of the adaptive Kalman filter; when the adaptive EKF algorithm is used, the calculation formulas are formulas (24) to (27) respectively; when the adaptive SPKF algorithm is used, the calculation formula are formulas (43) to (46) respectively; (8) Estimate the statistical characteristics of system noise; Calculate the mean estimate of the system noise in turn

and covariance estimates

When the adaptive EKF algorithm is used, the calculation formulas are formulas (28) and (29) respectively; when the adaptive SPKF algorithm is used, the calculation formulas are formulas (47) and (48) respectively; (9) Judging whether the filtering process is over; if it is necessary to continue filtering, k=k+1, and return to (4).

Images