CN107421550B - Earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging - Google Patents

Earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging Download PDF

Info

Publication number
CN107421550B
CN107421550B CN201710611461.6A CN201710611461A CN107421550B CN 107421550 B CN107421550 B CN 107421550B CN 201710611461 A CN201710611461 A CN 201710611461A CN 107421550 B CN107421550 B CN 107421550B
Authority
CN
China
Prior art keywords
earth
satellite
noise
lagrange
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710611461.6A
Other languages
Chinese (zh)
Other versions
CN107421550A (en
Inventor
杨静
卢帅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN201710611461.6A priority Critical patent/CN107421550B/en
Publication of CN107421550A publication Critical patent/CN107421550A/en
Application granted granted Critical
Publication of CN107421550B publication Critical patent/CN107421550B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/02—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by astronomical means

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Astronomy & Astrophysics (AREA)
  • Automation & Control Theory (AREA)
  • General Physics & Mathematics (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

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

Earth-Lagrange combined 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, and particularly relates to a method for realizing autonomous orbit determination of earth satellite constellations based on adaptive filtering by only utilizing inter-satellite ranging information in an earth-Lagrange (Lagrange) combined constellation.
Background
In recent decades, satellite navigation plays an increasingly important role in the field of national economy and military combat, and the improvement of the autonomous operation capability of the satellite has important significance in the aspects of reducing the ground measurement and control burden, reducing the satellite operation cost, improving the satellite survival capability, expanding the application potential of the spacecraft and the like. The autonomous orbit determination technology of the satellite is a precondition for realizing autonomous control of the satellite, plays a vital role in ensuring autonomous operation of the satellite, and becomes a development trend of the current satellite navigation and control technology.
The current autonomous orbit determination method of the constellation is mainly divided into two types of independent autonomous orbit determination methods of each satellite and integral autonomous orbit determination of the constellation. The former is to obtain the information such as the distance and direction of the satellite relative to other natural celestial bodies or navigation stars through various inertial devices, star sensors, navigation receivers and the like carried on the satellite to carry out on-line estimation on the self position, the system is simple in structure and small in operand, but the orbit determination precision is limited by the factors such as the irregularity degree of the surface of the natural celestial body and the measurement precision of the sensors, and the precision of the navigation level application is difficult to achieve by using the system alone; the latter makes full use of the relative measurement information among satellites in the constellation to realize the whole network orbit determination of the constellation, the short arc segment orbit determination has high precision, but the problem of 'loss rank' exists,namely, the overall direction rotation of the constellation in the inertial space cannot be determined only by using the inter-satellite relative measurement information, so that the orbit determination precision of the long arc segment of the constellation satellite gradually decreases along with the increase of time. In order to solve the problem of 'deficit rank', foreign scholars indicate that absolute information is introduced by 'uniqueness' of an orbit when a constellation or a satellite formation is highly asymmetric in a gravitational field, so that the problem of unobservable constellation overall rotation is solved. It has been found that the relative intensity of the third body's attractive force is greatest in many asymmetric perturbation gravitational fields, particularly at L1And L2Two Lagrange point neighbors, so the 'deficit rank' problem can be solved by joint autonomous orbit determination by combining the earth satellite constellation and the Lagrange point constellation. By using the autonomous orbit determination method, the complexity of an autonomous navigation system of the spacecraft can be greatly reduced, navigation measuring equipment is reduced, and the autonomous orbit determination precision is improved. Meanwhile, the Lagrange navigation constellation can also be used as a deep space exploration navigation base station to provide navigation information for other deep space detectors, so that the problems of the ground deep space exploration network in the aspects of deep space navigation accuracy, real-time performance, safety and the like are effectively solved, and therefore, the study on the autonomous orbit determination of the earth-Lagrange combined constellation is of great significance.
At present, the autonomous orbit determination method for earth-Lagrange joint constellation at home and abroad is mainly based on conventional algorithms such as batch processing algorithm, Extended Kalman Filter (EKF) algorithm and Unscented Kalman Filter (UKF) algorithm. However, batch processing algorithms are often used for post-processing and are not suitable for real-time tracking; although the conventional EKF and UKF algorithms adopt a recursion form which is convenient for real-time calculation, the prior information of system noise and measurement noise is required to be known, otherwise, the performance of the system is reduced, and filtering divergence is caused. Due to the fact that the actual stress condition of the in-orbit satellite is complex, particularly Lagrange satellites, accurate noise statistical characteristics are difficult to obtain, and irregular modeling errors caused by maneuvering and the like maintained by the satellite orbit can also cause changes of the statistical characteristics, so that the priori information is meaningless, and therefore the study of an autonomous navigation algorithm capable of being adaptively adjusted according to a measurement information sequence has important significance for improving the robustness of an autonomous orbit determination system and ensuring the autonomous orbit determination accuracy.
Disclosure of Invention
In order to solve the problems, the invention provides an earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging, which only utilizes satellite-to-satellite distance measurement information obtained by inter-satellite links, utilizes an Adaptive Nonlinear Filter (ANF) algorithm and adopts a centralized structure to fuse effective measurement information, thereby realizing the resolving of the position and the speed of each satellite in the earth-Lagrange combined constellation. The constellation orbit determination method comprises the steps of firstly, establishing a system model of earth-Lagrange combined constellation autonomous orbit determination by adopting a centralized structure, wherein a state equation of an earth navigation satellite is established mainly based on a subject two-body problem dynamic model under an earth center inertia rectangular coordinate system, the state equation of the Lagrange navigation satellite is established mainly based on a circular restrictive three-body problem dynamic model under a mass center convergence coordinate system, and a measurement equation adopts an inter-satellite distance measurement model under the earth center inertia rectangular coordinate system; then, effective measurement information is fused by adopting a self-adaptive nonlinear filtering algorithm to realize the autonomous orbit determination of the earth-Lagrange combined constellation; and finally, verifying the effectiveness of the method by taking a combined constellation consisting of 3 Medium Earth Orbit (MEO) satellites and 2 Lagrange satellites as an example.
An earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging specifically comprises the following steps:
the method comprises the following steps: under a centralized structure, establishing a state equation of an earth-Lagrange combined constellation autonomous orbit determination system;
step two: under a centralized structure, establishing a measurement equation of an earth-Lagrange combined constellation autonomous orbit determination system;
step three: under a centralized structure, determining a filtering method for realizing orbit parameter estimation according to the established earth-Lagrange combined constellation autonomous orbit determination system model;
step four: under a centralized structure, the earth-Lagrange combined constellation autonomous orbit determination method based on the selected filtering algorithm is specifically realized.
The invention has the advantages that:
(1) absolute information is provided for the system by introducing the Lagrange satellite, the problem of 'loss rank' is effectively solved, and the long-term autonomous orbit determination precision of the earth-Lagrange combined constellation is improved;
(2) only the inter-satellite distance measurement information is utilized, so that the complexity of system equipment is reduced;
(3) the coverage rate of the Lagrange navigation constellation to the earth navigation constellation is high, and the autonomous orbit determination precision is improved;
(4) the statistical characteristic of the system noise is estimated on line in real time by using the self-adaptive nonlinear filtering algorithm, and the precision and the stability of the filtering algorithm for autonomous orbit determination are improved, so that the autonomous orbit determination precision is improved.
Drawings
FIG. 1 is a schematic diagram of an earth-Lagrange joint constellation autonomous orbit determination system based on inter-satellite ranging;
FIG. 2 is a diagram of a method for autonomous orbit determination of a combined earth-Lagrange constellation based on inter-satellite ranging;
FIG. 3 is a flow chart of an autonomous navigation algorithm based on adaptive nonlinear filtering;
FIG. 4 shows the system accuracy as bQ、bC、bqA trend graph;
FIG. 5 is a diagram illustrating a position error curve of an exemplary GNSS satellite.
Detailed Description
The present invention will be described in further detail with reference to the accompanying drawings and examples.
The invention relates to an earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging, the system is shown in figure 1, the joint constellation orbit determination system is shown in the figure to be composed of an earth navigation constellation and inter-satellite link ranging thereof, a Lagrange navigation constellation and inter-satellite link ranging thereof and inter-satellite link ranging between two constellation satellites, wherein the earth navigation constellation is composed of m medium/high orbit satellites, and the serial numbers are 1,2, …, i and … m in sequence; lagrange's navigation constellation is operated at L by n pieces1And L2And the satellite composition on the Halo periodic orbit near the point is numbered as 1,2, …, k and … n in sequence. Earth-Lagrange combined constellation autonomous system based on inter-satellite rangingThe basic implementation process of the orbit determination method is shown in figure 2, which shows that the basic implementation process of the orbit determination method mainly comprises measurement information acquisition and constellation orbit parameter estimation, wherein the acquired measurement information comprises three types of distance measurement information including inter-satellite distance measurement information inside a terrestrial navigation constellation, inter-satellite distance measurement information inside a Lagrange navigation constellation and distance measurement information between a terrestrial navigation satellite and a Lagrange navigation satellite, ① establishes a state equation of a combined constellation orbit determination system based on a simplified orbit dynamics model of the terrestrial satellite and the Lagrange satellite when constellation orbit parameter estimation is carried out, ② establishes a measurement equation of three types of distance measurement information, ③ adopts a system model formed by the state equation and the measurement equation and adopts a proper self-adaptive nonlinear filtering algorithm to estimate orbit parameters, and finally outputs the position and the speed of the satellite in the estimated constellation, and adopts a rectangular coordinate description method (the position and speed components are used for representing the current orbit state of the satellite) under a geocentric inertial coordinate system and a mass center complex and a coordinate system, and the self-adaptive nonlinear filtering algorithm is combined with the earth orbit measurement information estimation method, and the measurement method is specifically comprises the steps of establishing a combined orbit determination system according to the Langrange orbit measurement equation:
the method comprises the following steps: under a centralized structure, establishing a state equation of an earth-Lagrange combined constellation autonomous orbit determination system;
the accurate system model is one of the main factors for guaranteeing the estimation precision of the constellation operation state parameters. The invention establishes a system state equation based on a satellite orbit dynamics model.
(1) Dynamic model of earth navigation satellite orbit
The earth navigation satellite is mainly a medium and high orbit satellite, is far more greatly influenced by the gravity of the center of mass of the earth than other acting forces, and generally adopts a dynamic model of a shot two-body problem. In order to meet the requirements of high efficiency and high precision of calculation, a main perturbation item is selected for modeling according to the perturbation force influence. For medium and high orbit satellites, the perturbation influence caused by the earth non-spherical symmetry can reach 10 in a plurality of perturbation items4Meter class, wherein the perturbation influence of gravitational field above 3 orders and J2The phase ratio is very smallMany, can be ignored; the influence of sun-moon gravity perturbation can reach 103Meter level; the perturbation effect of sunlight pressure can reach 102Meter level; while the effects of other perturbations such as tidal, relativistic and albedo perturbations are at 10-1The meter level is lower than the meter level and can be ignored; meanwhile, the atmospheric density at the high rail is very low, and the atmospheric resistance perturbation influence can be ignored. Therefore, the invention mainly considers the gravity of the center of mass of the earth and J2And an orbit dynamics model formed by three main perturbation forces of terms, daily gravitation and monthly gravitation is used for establishing a system state equation.
Under the earth-centered rectangular inertial coordinate system, establishing an orbit determination state equation for the earth satellite i to be detected in the combined constellation by combining an orbit dynamics model as follows:
Figure BDA0001359618900000041
wherein, the state vector corresponding to the earth satellite i to be measured is
Figure BDA0001359618900000042
Here [ x ]Ei,yEi,zEi]TAnd
Figure BDA0001359618900000043
respectively, a position vector and a velocity vector of the earth satellite. FEAs a function of the state of the system, WEi(t) is the system noise, satisfying the mean qEi(t) variance QEi(t) white Gaussian noise. Formula (1) can be written in further detail as:
Figure BDA0001359618900000051
wherein the content of the first and second substances,
Figure BDA0001359618900000052
and
Figure BDA0001359618900000053
is the noise component of x, y and z three-axis velocity,
Figure BDA0001359618900000054
and
Figure BDA0001359618900000055
is the acceleration equivalent noise component of the x, y and z three axes, including unmodeled perturbation term and noise term;
Figure BDA0001359618900000056
and
Figure BDA0001359618900000057
a modeled velocity vector and an acceleration vector corresponding to the earth satellite i, respectively, wherein the acceleration mainly takes into account the earth central gravitational acceleration a0,Ei、J2Acceleration of perturbation term
Figure BDA0001359618900000058
Acceleration a of sun-moon attractionMS,EiThe specific expression is as follows:
1) earth central gravitational acceleration a of satellite i considered by the invention0,EiSatisfies the following conditions:
Figure BDA0001359618900000059
in the formula (I), the compound is shown in the specification,
Figure BDA00013596189000000510
is the earth center distance of the satellite i to be measured, mu is the earth mass METhe product of the gravity constant G, the earth's gravity constant.
2)J2Acceleration of perturbation term
Figure BDA00013596189000000511
Satisfies the following conditions:
Figure BDA00013596189000000512
in the formula, ReIs the equatorial radius of the earth, J2Are second order band harmonic coefficients.
3) Perturbation acceleration a caused by sun-moon gravitationMS,EiSatisfies the following conditions:
Figure BDA0001359618900000061
in the formula (x)S,yS,zS) And (x)M,yM,zM) Respectively are three-dimensional position coordinates of the sun and the moon under a geocentric rectangular inertial coordinate system,
Figure BDA0001359618900000062
and
Figure BDA0001359618900000063
the distance from satellite i to the sun and moon respectively,
Figure BDA0001359618900000064
and
Figure BDA0001359618900000065
distance from the earth's center to the sun and moon, mu, respectivelySAnd muMThe solar attraction constant and the lunar attraction constant, respectively.
(2) Lagrange navigation satellite orbit dynamics model
Different from the earth navigation satellite, the Lagrange satellite runs near the Lagrange point of the earth-moon system and is subjected to the action of the earth central gravity and the action of the moon central gravity in a similar magnitude, so that the invention adopts a circular restrictive three-body model when describing the movement of the Lagrange satellite.
Under a centroid convergence coordinate system, a state equation for orbit determination of a lunar satellite k to be measured in a Lagrange constellation can be established by combining an orbit dynamics model as follows:
Figure BDA0001359618900000066
in the formula, the corresponding state vector of Lagrange lunar satellite k to be measured is
Figure BDA0001359618900000067
Here [ x ]Lk,yLk,zLk]TAnd
Figure BDA0001359618900000068
the dimensionless position vector and the velocity vector of the star respectively have characteristic quality, characteristic length and characteristic time shown in formula (7), wherein M isEAnd MMThe earth mass and the moon mass respectively, and L is the earth-moon distance. FLAs a function of the state of the system, WLk(t) is the system noise, satisfying the mean qLk(t) variance QLk(t) white Gaussian noise.
Figure BDA0001359618900000069
Equation (6) can be written in further detail as:
Figure BDA0001359618900000071
wherein
Figure BDA0001359618900000072
For the modeled velocity vector, μ, corresponding to Lagrange satellite k0＝MM/(MM+ME)，ΔE,Lk、ΔM,LkThe geocentric distance and the lunar centroidal distance of the Lagrange satellite k,
Figure BDA0001359618900000073
and
Figure BDA0001359618900000074
is the x, y and z three-axis velocity noise component of Lagrange satellite k,
Figure BDA0001359618900000075
and
Figure BDA0001359618900000076
is the x, y, z triaxial acceleration equivalent noise component of Lagrange satellite kIncluding unmodeled perturbation terms and noise terms.
(3) State equation of earth-Lagrange combined constellation autonomous orbit determination system
Assume that the entire autonomous orbiting system contains m earth satellites and n Lagrange satellites. According to the state equations (1) and (6) for single satellite orbit determination, the constellation orbit determination state equation with a centralized structure is established as follows:
Figure BDA0001359618900000077
namely:
Figure BDA0001359618900000078
in the formula (I), the compound is shown in the specification,
Figure BDA0001359618900000079
f is the state vector of the whole constellation orbit system, w (t) represents the system noise with mean value q (t) and variance q (t).
Step two: under a centralized structure, establishing a measurement equation of an earth-Lagrange combined constellation autonomous orbit determination system;
the observation information used by the earth-Lagrange combined constellation autonomous orbit determination system only has inter-satellite distance measurement, but because the dynamic models of two types of satellites are established in different coordinate systems, the inter-satellite distance measurement model of the whole constellation is divided into three types:
(1) inter-satellite distance measurement model between earth navigation satellites
For the earth navigation satellite, the observed distance between the satellite i and the satellite j is recorded as rhoEi,EjThen, there are:
ρEi,Ej＝hE(XEi,XEj)+vEi,Ej，i≠j (11)
wherein
Figure BDA0001359618900000081
vEi,EjIs the distance equivalence between satellite i and satellite jMeasurement noise, including modeling error and metrology noise.
(2) Inter-satellite distance measurement model between Lagrange navigation satellites
Let the observation of the distance between Lagrange satellites i and j be ρLi,LjThe method comprises the following steps:
ρLi,Lj＝hL(XLi,XLj)+vLi,Lj，i≠j (12)
wherein
Figure BDA0001359618900000082
vLi,LjNoise, including modeling error and metrology noise, is measured equivalently for the distance between Lagrange satellites i and j.
(3) Measurement model between earth navigation satellite and Lagrange satellite
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, and there are:
Figure BDA0001359618900000083
wherein the content of the first and second substances,
Figure BDA0001359618900000084
is the position vector of Lagrange satellite in the inertial coordinate system of earth center, which is the state vector X of Lagrange satellite in the coordinate system of center of mass convergenceLkA dimensionless position vector [ x ] ofLk,yLk,zLk]TThe conversion relationship is:
Figure BDA0001359618900000085
wherein, R represents a rotation matrix from a centroid convergence coordinate system to a geocentric inertia coordinate system, and the expression is as follows:
Figure BDA0001359618900000091
wherein u isM＝ωM+θM，iM、ΩM、ωMAnd thetaMRespectively showing the orbit inclination angle, the ascent crossing declination, the argument of the perigee and the true perigee angle of the orbit of the moon around the earth.
(4) Measurement equation of earth-Lagrange combined constellation autonomous orbit determination system
According to the above measurement equations (11), (12) and (13), the constellation orbit determination measurement equation for establishing the centralized structure is as follows:
Z(t)＝h[X(t),t]+V(t) (16)
namely:
Figure BDA0001359618900000092
where V (t) is zero mean and variance is R (t). The formula (17) gives an observation equation under a general condition, and in practical application, some distance measurement quantities are influenced by factors such as the shielding of the earth or the moon, so that the situation is not measurable, and corresponding distance measurement information needs to be removed from the formula (17).
Step three: determining a filtering method for realizing orbit parameter estimation according to the established earth-Lagrange combined constellation autonomous orbit determination system model based on inter-satellite distance measurement;
the state of the nonlinear system needs to be estimated using a nonlinear filtering method. Conventional nonlinear filtering methods, such as EKF, UKF, volume filtering, particle filtering, etc., can obtain a more ideal result under the condition that the prior statistical information of the noise is accurately known. In the present system, however, the equivalent system noise WkIn addition to containing unmodeled errors, discretization errors (and linearization errors) in the filtering algorithm are also contained, which makes the noise have time-varying non-gaussian statistical properties and difficult to acquire accurately. For this reason, when applying these nonlinear filtering algorithms, the noise is generally approximated as a stationary random white noise sequence, i.e. Q is consideredkQ and empirically determined. In order to ensure that the state estimation error can be converged to a stable value in the whole filtering process, the value of Q is usually large,namely, Q ≧ max { Q ≧ is satisfiedkWhich makes system performance difficult to optimize.
In order to solve the problem, the invention combines the idea of Sage-Husa adaptive algorithm to perform online estimation and correction on the statistical characteristic of the system noise on the basis of the conventional nonlinear filtering algorithm so as to improve the filtering precision and stability. In consideration of engineering practicability, the satellite orbit parameter estimation method adopts a self-adaptive EKF method or a self-adaptive Sigma Point Kalman Filter (SPKF) method to realize satellite orbit parameter estimation.
Firstly, when the system has higher requirement on the operation speed or the system storage space is smaller, the state estimation is carried out on the discretized nonlinear system by adopting a self-adaptive EKF method. The basic idea of the self-adaptive EKF is to utilize a Taylor series expansion method to carry out discretization and linearization processing on a continuous nonlinear state equation and an observation equation, then estimate and correct the statistical characteristics of system noise and observation noise in real time through a time-varying noise statistical estimator on the basis of classical Kalman filtering, and then carry out state estimation. The method comprises the following specific steps:
(1) time updating
Calculating a one-step prediction state:
Figure BDA0001359618900000101
in the formula (I), the compound is shown in the specification,
Figure BDA0001359618900000102
is an estimate of the mean value of the system noise at time k,
Figure BDA0001359618900000103
for one-step recursion of states using kinetic equations, one method is by calculation using equation (19) after discretization and linearization, where
Figure BDA0001359618900000104
For the state estimation at the time instant k,
Figure BDA0001359618900000105
estimating the one-step prediction state at the moment (k +1), wherein delta t is an integral step length; another method is directly obtained by numerical integration calculation of the continuous state equation (9) by a 4-order Runge-Kutta (RK4) method.
Figure BDA0001359618900000106
Calculating a state transition matrix:
Figure BDA0001359618900000107
in the formula phik+1/kIs the state transition matrix from time k to time (k + 1).
Calculating a one-step prediction error covariance matrix:
Figure BDA0001359618900000108
in the formula, Pk+1/kError covariance matrix, P, for one-step prediction of stateskThe error covariance matrix of the state is estimated for time k,
Figure BDA0001359618900000109
is the estimated value of the system noise variance matrix at the k moment.
Considering that the system state estimation error is not converged in a period of initial filtering and the samples used for estimating the system noise are few, the estimated noise statistical characteristic has large fluctuation and is not accurate enough, so that the filtering convergence is slow, and even the filtering divergence is caused in serious cases. Therefore, in a period of time (0-T) after the initial filtering, the algorithm does not apply the estimated noise statistical characteristics to the time updating process, but adopts the prior statistical information of the noise, and the time updating algorithm at this time is as follows:
Figure BDA0001359618900000111
Figure BDA0001359618900000112
(2) and (3) measurement updating:
calculating a gain matrix:
Figure BDA0001359618900000113
in the formula, Kk+1Is the gain matrix at time (k +1),
Figure BDA0001359618900000114
to predict state estimation based on one step
Figure BDA0001359618900000115
Jacobian matrix, R, of calculated observation vectors hk+1The noise covariance matrix is measured for time (k + 1).
Calculating an innovation vector:
Figure BDA0001359618900000116
in the formulak+1Is the innovation vector at time (k +1), Zk+1Is the measurement vector at time (k + 1).
Calculating a state estimation value:
Figure BDA0001359618900000117
calculating a state estimation error covariance matrix:
Figure BDA0001359618900000118
calculating an estimated value of the statistical characteristic of the system noise:
Figure BDA0001359618900000119
Figure BDA00013596189000001110
Figure BDA00013596189000001111
wherein the content of the first and second substances,
Figure BDA00013596189000001112
Figure BDA0001359618900000121
Figure BDA0001359618900000122
in the formula 0 < bq＜1，0＜bQ＜1，0＜bC< 1 are forgetting factors for calculating Q, Q, C, respectively. As the value of k is increased, it is,
Figure BDA0001359618900000123
Figure BDA0001359618900000124
and
Figure BDA0001359618900000125
respectively tend to (1-b)q)、(1-bQ) And (1-b)C) Thus b isq、bQAnd bCThe larger the proportion of the current information.
Figure BDA0001359618900000126
Focusing on the change of the noise mean value q of the tracking system, the overlarge value can cause
Figure BDA0001359618900000127
The value is small, and the filtering divergence can be caused if the value is too small;
Figure BDA0001359618900000128
emphasis on smoothing
Figure BDA0001359618900000129
If the value is too small, the smoothing effect is poor, and if the value is too large, the effect is poor
Figure BDA00013596189000001210
The filter divergence risk is increased due to small size;
Figure BDA00013596189000001211
the method is used for tracking the change of the variance C, the fluctuation of the estimated value is increased when the value is too small, and the tracking is delayed when the value is too large.
In practical applications, although the calculations of equations (29) and (30) have some smoothing effect, they still occur
Figure BDA00013596189000001212
In the case of negative determination, in order to ensure the stability of the filtering, the method of formula (29) is performed in the case of negative determination
Figure BDA00013596189000001213
The diagonal elements of the term take the absolute value and the off-diagonal elements take 0 to approximate.
Secondly, when the system pays more attention to the accuracy index and has no special limit on the operation speed and the occupied storage space, the state estimation is carried out on the continuous nonlinear system by adopting a self-adaptive Sigma Point Kalman Filter (SPKF) method. The basic idea of the self-adaptive SPKF is to directly transfer Sigma sampling points obtained according to a deterministic sampling strategy through a nonlinear function, then estimate and correct the statistical characteristics of system noise and observation noise in real time through a time-varying noise statistical estimator on the basis of a classic Kalman filtering algorithm framework, and then carry out state estimation. The method comprises the following specific steps:
1) selection of Sigma Point sampling strategy
Currently, a commonly used Sigma point sampling strategy has a plurality of sampling modes such as symmetric sampling, simplex sampling, central differential sampling, volume sampling and the like. Selecting a sampling strategy and calculating the weight coefficient W of the SPKFi m、Wi cWhere i is 0,1, …, L is the number of sampling points determined by the specific sampling strategy, Wi mFor mean value estimation, Wi cFor variance estimation.
2) And (3) time updating:
according to the sampling strategy selected in 1) and the state estimation value at the k moment
Figure BDA00013596189000001214
Sum estimation error covariance PkSigma point set chi for calculating k timei,k(i＝0,1,…,L)。
Calculate one-step predicted values for Sigma points:
χi,k+1/k＝F[χi,k,k]，i＝0,1,…,L (34)
wherein, χi,kIs the value of the ith Sigma point at time k, χi,k+1/kThe predicted value of the ith Sigma point at one step at the moment k +1 is obtained.
Calculating a one-step prediction state:
Figure BDA0001359618900000131
in the formula (I), the compound is shown in the specification,
Figure BDA0001359618900000132
is an estimate of the mean value of the system noise at time k,
Figure BDA0001359618900000133
the state estimate is predicted for one step at time (k + 1).
Calculating a state one-step prediction error covariance matrix:
Figure BDA0001359618900000134
in the formula, Pk+1/kFor the error covariance matrix of the one-step predicted state,
Figure BDA0001359618900000135
for the system noise variance at time kAn estimate of the array;
in the initial time 0-T of filtering, adopting the prior statistical information of noise, and updating the time as follows:
Figure BDA0001359618900000136
Figure BDA0001359618900000137
3) and (3) measurement updating:
according to the sampling strategy selected in 1) and
Figure BDA0001359618900000138
and Pk+1/kUpdating a Sigma point set;
estimating the measured value:
γi，k+1/k＝h(χi，k+1/k) (39)
Figure BDA0001359618900000139
wherein gamma isi，k+1/kIs an estimate of the measurement vector at time (k +1) calculated from the ith Sigma point,
Figure BDA00013596189000001310
is a weighted estimate of the measurement vector at time (k + 1).
Calculating a gain matrix:
Figure BDA00013596189000001311
Figure BDA00013596189000001312
Figure BDA00013596189000001313
wherein R isk+1For time measurement of (k +1)Noise covariance matrix, Py,k+1/kMeasuring the covariance of the vector estimates for the (k +1) th moment, Pxy,k+1/kIs the covariance of the (K +1) th measured vector estimate and the state one-step predictor, Kk+1Is the gain matrix at time (k + 1).
Calculating an innovation vector:
Figure BDA0001359618900000141
in the formulak+1Is the innovation vector at time (k +1), Zk+1Is the measurement vector at time (k + 1);
calculating a state estimate at time (k +1)
Figure BDA0001359618900000142
Figure BDA0001359618900000143
Computing a state estimation error covariance matrix P at time (k +1)k+1：
Figure BDA0001359618900000144
4) Calculating an estimated value of the statistical characteristic of the system noise:
Figure BDA0001359618900000145
Figure BDA0001359618900000146
Figure BDA0001359618900000147
wherein the content of the first and second substances,
Figure BDA0001359618900000148
Figure BDA0001359618900000149
Figure BDA00013596189000001410
in the formula 0 < bq＜1，0＜bQ＜1，0＜bC< 1 are forgetting factors for calculating Q, Q, C, respectively.
When it occurs
Figure BDA00013596189000001411
In the case of negative determination, in the case of equation (48)
Figure BDA00013596189000001412
The diagonal elements of the term take the absolute value and the off-diagonal elements take 0.
Step four: under a centralized structure, the earth-Lagrange combined constellation autonomous orbit determination algorithm based on a selected filtering method is specifically realized;
referring to fig. 4, a flow chart of the self-adaptive nonlinear filtering algorithm for implementing the earth-Lagrange joint constellation autonomous orbit determination algorithm provided by the invention is as follows:
the specific process is as follows:
(1) initializing data;
the initialized parameters include: k is 0, filtering step h, noise switching parameter T, initial state
Figure BDA0001359618900000151
Initial error covariance matrix P0Initial system noise mean
Figure BDA0001359618900000152
Initial system noise variance matrix
Figure BDA0001359618900000153
Measuring noise variance matrix R and forgetting factor bq、bQ、bC；
The noise switching parameter T is a time node for switching noise statistical characteristics in the time updating algorithm: when k is less than T, using the prior statistical property of the noise to update time; when k is larger than or equal to T, updating time by using the estimated noise statistical characteristics; initial state
Figure BDA0001359618900000154
Initial error covariance matrix P0Initial system noise variance matrix
Figure BDA0001359618900000155
Determining according to the requirement; initial system noise mean
Figure BDA0001359618900000156
Taking zero; the measurement noise variance matrix R is determined according to the performance of the measurement equipment; forgetting factor bq、bQThe value is between 0.99 and 1.0, bCThe value is between 0.9 and 0.99;
(2) if the self-adaptive EKF algorithm is adopted, directly turning to (4); if the self-adaptive SPKF algorithm is adopted, turning to (3);
(3) calculating weight coefficient W according to selected Sigma point sampling strategyi m、Wi c；
(4) If k is less than T, updating time by adopting the prior statistical property of the noise, and turning to (5); otherwise, updating time by adopting the estimated noise statistical characteristics, and turning to the step (6);
(5) time updating using a priori statistical properties of the noise; turning to (7) after completion;
state estimation using time k
Figure BDA0001359618900000157
And error covariance matrix PkSequentially calculating the state one-step predicted value as the state initial value and the error covariance matrix at the k moment
Figure BDA0001359618900000158
And the prediction error covariance matrix Pk+1/kImplementing adaptive Kalman filtering estimationAnd (4) updating the time of the method. When the self-adaptive EKF algorithm is adopted, the calculation formulas are respectively shown as a formula (22) and a formula (23); when the self-adaptive SPKF algorithm is adopted, the calculation formulas are respectively shown as a formula (37) and a formula (38);
(6) time updating using the estimated noise statistics;
state estimation using time k
Figure BDA0001359618900000159
And error covariance matrix PkSequentially calculating the state one-step predicted value as the state initial value and the error covariance matrix at the k moment
Figure BDA00013596189000001510
And the prediction error covariance matrix Pk+1/kAnd time updating of the self-adaptive Kalman filtering estimation method is realized. When the self-adaptive EKF algorithm is adopted, the calculation formulas are respectively an expression (18) and an expression (21); when the self-adaptive SPKF algorithm is adopted, the calculation formulas are respectively shown as the formula (35) and the formula (36);
(7) updating the measurement;
sequentially calculating a gain matrix K at (K +1) timek+1Innovation matrix vk+1State estimation value
Figure BDA0001359618900000161
And an error covariance matrix Pk+1And realizing the measurement updating of the self-adaptive Kalman filtering. When the self-adaptive EKF algorithm is adopted, the calculation formulas are respectively shown as formulas (24) to (27); when the self-adaptive SPKF algorithm is adopted, the calculation formulas are respectively expressed as formulas (43) to (46);
(8) estimating the statistical characteristics of system noise;
calculating mean value estimated value of system noise in turn
Figure BDA0001359618900000162
Sum covariance estimate
Figure BDA0001359618900000163
When the adaptive EKF algorithm is adopted, the calculation formulas are respectively shown as a formula (28) and a formula (29); when the adaptive SPKF algorithm is employed,the calculation formulas are respectively formula (47) and formula (48);
(9) judging whether the filtering process is finished or not; filtering is continued as necessary, k is k +1, and (4) is returned.
Firstly, researching the problem of 'deficit rank' existing in inter-satellite distance measurement, and introducing Lagrange satellites to provide an absolute reference; and then, under a centralized structure, realizing the autonomous orbit determination of the earth-Lagrange combined constellation by adopting a self-adaptive nonlinear filtering estimation method.
The invention provides an earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging, which comprises the steps of firstly establishing a state equation of an earth-Lagrange combined constellation in a centralized structure; then establishing a measurement equation of the earth-Lagrange combined constellation; and finally, realizing the earth-Lagrange combined autonomous orbit determination by adopting a self-adaptive nonlinear filtering method.
Example (b):
taking the geosphere-Lagrange constellation as an example, the constellation consists of 3 MEO satellites and 2 Lagrange satellites. MEO satellite number E1～E3(ii) a Lagrange satellite number L1～L2. In simulation, it is found that the orbit determination precision of each satellite in the MEO constellation is approximate, and the orbit determination precision of each satellite in the Lagrange constellation is also approximate, so that the number E is only used1And MEO satellite number L1The Lagrange satellite of (1) is taken as an example for illustration.
Simulation conditions are as follows:
and simulating the actual operation condition of the constellation by adopting the STK to generate nominal orbit data in the operation process of the constellation.
(1) Nominal data generation parameter settings
The simulation time is 2016 year 1 month 12:00: 00-2016 year 6 month 28 month 12:00:00, duration is 180d, step length T is 1s, an orbit dynamics model of an earth satellite adopts a High-precision orbit prediction model (HPOP), wherein The earth model adopts a WGS84-EGM96 model, 21-order aspheric perturbation is considered, and other perturbation items comprise solar attraction, lunar attraction and solar pressure (light pressure coefficient C)r1.0) and atmospheric resistance (atmospheric resistance coefficient C)d2.2); the ratio S/m of the cross-sectional area of the satellite to the mass of the satellite is 0.02m2In terms of/kg. The orbital dynamics model of the Lagrange satellite employs a circular restrictive trisomy problem model.
(2) Basic filter parameter setting
The state equation of the earth satellite also adopts a subject disomic problem model, but only the J2 term and the gravity of the sun and the moon are considered in the perturbation term; the Lagrange satellite state equation adopts a circular restrictive three-body problem model, and the filtering period T is 100 s.
Under a centralized structure, the initial value of the system noise variance matrix of the earth-Lagrange combined constellation autonomous orbit determination is Q0＝diag([QE1，0QE2，0QE3,0QL1,0QL2,0]) Wherein
Figure BDA0001359618900000171
Figure BDA0001359618900000178
Figure BDA0001359618900000172
σLx0＝σLy0＝σLz0＝10-10[L]，
Figure BDA0001359618900000179
Measure the variance matrix of noise as
Figure BDA0001359618900000173
σEi,Ej＝σEi,Lk＝σLk,Ll5 m; initial state estimation error covariance matrix P0＝diag([PE1,0PE2,0PE3,0PL1,0PL2,0])，
Wherein
Figure BDA0001359618900000174
σEx,0＝σEy,0＝σEz,0＝1m，
Figure BDA00013596189000001710
σLx,0＝σLy,0＝σLz,0＝1m，
Figure BDA00013596189000001711
(3) Evaluation method
Because the system state and the estimation error in the actual system are random quantities, after the estimation error of each simulation moment is obtained, the absolute orbit determination and relative positioning effects of the satellites in the constellation are evaluated by adopting a statistical method (calculating the root mean square error of error data). And counting errors by using the estimation state after the convergence of the filtering estimation as a sample. The root mean square error calculation formula is as follows:
Figure BDA0001359618900000176
wherein the content of the first and second substances,
Figure BDA0001359618900000177
is an estimate of time i, XiIs the true value at the ith moment, and n is the total data.
For the convenience of quantitative analysis, the distance error and the velocity error of satellite orbit are respectively taken as the sum of squares of three-axis position errors and the sum of squares of three-axis velocity errors, and the calculation formula of the root mean square error is as follows:
satellite distance estimation error:
Figure BDA0001359618900000181
satellite rate estimation error:
Figure BDA0001359618900000182
the implementation mode is as follows: the state vector of the centralized earth-Lagrange combined constellation autonomous orbit determination system model consists of the position and the speed of five stars, an inter-satellite distance measurement mode is adopted, and the improvement of the orbit determination precision by the self-adaptive EKF is mainly researched.
The method comprises the following steps: establishing state equation of earth-Lagrange combined constellation autonomous orbit determination system
Under the earth-centered rectangular inertial coordinate system, the state vector of the earth navigation satellite i is
Figure BDA0001359618900000183
Wherein [ x ]Ei,yEi,zEi]TAnd
Figure BDA0001359618900000184
and respectively establishing an orbit determination state equation of the satellite Ei by combining the position vector and the velocity vector of the satellite with an orbit dynamics model formula (1):
Figure BDA0001359618900000185
in the formula, the gravity of the earth center, J, of the satellite Ei2The acceleration caused by the perturbation force and the gravity of the sun and the moon is respectively a0,Ei、aJ2,EiAnd aMS，EiThe following are specifically defined by the following formulae (3) to (5). WEi(t) represents the system noise of the satellite Ei, including unmodeled perturbations, satisfying E [ W ]Ei(t)]＝0，E[WEi(t)WEi(τ)]＝QEi(t) (t- τ) and (t- τ) is the Dirac function, i.e.
Figure BDA0001359618900000186
Lagrange satellite k corresponds to a state vector of
Figure BDA0001359618900000187
Wherein [ x ]Lk,yLk,zLk]TAnd
Figure BDA0001359618900000188
the dimensionless position vector and the velocity vector of the satellite are respectively combined with the orbit dynamics model formula (6) to establish the satellite LkOrbit determination state equation:
Figure BDA0001359618900000191
in the formula, WLk(t) represents the system noise of the satellite Ei, including unmodeled perturbations, satisfying E [ W ]Lk(t)]＝0，E[WLk(t)WLk(τ)]＝QLk(t)(t-τ)。
According to the state equation of the single satellite, the five-star constellation state equation with a centralized structure can be established as follows:
Figure BDA0001359618900000192
in the formula (I), the compound is shown in the specification,
Figure BDA0001359618900000193
a state vector containing all states of five satellites in the constellation, wherein W (t) represents system noise and satisfies E [ W (t)]＝q(t)，E[W(t)W(τ)]＝Q(t)(t-τ)。
Step two: measurement equation for establishing earth-Lagrange combined constellation autonomous orbit determination system
Under a centralized structure, a five-star constellation adopts the following measurement equation of inter-satellite link measurement:
Figure BDA0001359618900000194
in the formula, hEIs a measurement equation between three earth satellites, hELFor the measurement equation between earth satellites and Lagrange satellites, hLSee equations (11) to (13) for the measurement equations between two Lagrange satellites.
V(t)＝[ΔEi,Ej…ΔEi,Lk…ΔLk,Ll…]TThe pseudo range observation noise vector is a pseudo range observation noise vector which is residual after error compensation of ionosphere delay, troposphere delay and the like, and is white noise with zero mean and standard deviation of R (t).
Step three: and under a centralized structure, determining a filtering method for realizing orbit parameter estimation according to the established earth-Lagrange combined constellation autonomous orbit determination system model.
The present example requires fast algorithm computation and low computational resource requirements, and therefore employs an adaptive EKF algorithm.
Step four: under a centralized structure, the combined autonomous orbit determination of the earth-Lagrange five-star constellation is realized by adopting a self-adaptive EKF estimation method, and the performance of the combined autonomous orbit determination of the earth-Lagrange five-star constellation is simulated and analyzed.
And (3) realizing the earth-Lagrange five-star constellation combined autonomous orbit determination by adopting a self-adaptive EKF estimation method by combining the earth-Lagrange five-star constellation centralized combined autonomous orbit determination flow chart (figure 4) and the step four of the invention content.
(1) The initialized parameters comprise k is 0, the filtering step length h is 100s, the threshold value T is (10 days × 86400 s/day)/100 s, and the initial state
Figure BDA0001359618900000201
Initial error covariance matrix P0Initial system noise mean
Figure BDA0001359618900000202
Initial system noise variance matrix
Figure BDA0001359618900000203
Measuring noise variance matrix R5 m and forgetting factor bq＝0.9999、bQ＝0.9999、bC＝0.90。
Note: because b isqThe filtering divergence is caused by over-small value, so that a larger value is taken as an initial value, and b is adjusted firstlyQAnd bCTo an optimum value, finally b is adjustedq。
(2) If k is less than T, updating time by adopting the prior statistical property of the noise, and turning to the step (3); otherwise, time updating is carried out by adopting the estimated noise statistical characteristics, and the step (4) is carried out.
(3) Temporal updating of the a priori statistical properties of the noise is employed. State estimation using time k
Figure BDA0001359618900000204
And error covariance matrix PkAs the initial state value and error covariance matrix at the time k, the sum of the equations (22) andequation (23) separately calculates the state one-step prediction values
Figure BDA0001359618900000205
And the prediction error covariance matrix Pk+1/kTime updating of the adaptive EKF estimation method is achieved. Go to (5).
(4) A temporal update of the estimated noise statistics is employed. State estimation using time k
Figure BDA0001359618900000206
And error covariance matrix PkAs the initial state value and error covariance matrix at the time k, the one-step predicted state values are calculated from the equations (18) and (21), respectively
Figure BDA0001359618900000207
And the prediction error covariance matrix Pk+1/kTime updating of the adaptive EKF estimation method is achieved.
(5) Adaptive EKF measurement updates. Sequentially calculating gain matrix K at time (K +1) according to equations (24) to (27)k+1Innovation matrix vk+1State estimation value
Figure BDA0001359618900000208
And an error covariance matrix Pk+1And realizing the measurement update of the EKF.
(6) The system noise statistics are estimated. The mean value estimation value of the system noise is calculated according to equations (28) and (29), respectively
Figure BDA0001359618900000209
Sum covariance estimate
Figure BDA0001359618900000211
(7) And judging whether the filtering process is finished or not. Filtering is continued as necessary, k is k +1, and (2) is returned.
Step five: adjusting a forgetting factor bQ、bC、bq。
Combining the effects and steps of the invention on three forgetting factors in step threeB is adjusted in sequence according to the approximate range of the forgetting factor given in the fourth stepQ、bC、bqAnd repeatedly executing the filtering process in the third step to enable the system to reach relative optimization, wherein the specific adjusting steps are as follows:
(1) fastening of bQ＝0.9999、bq0.9999, gradually increase bCDetermining b from the filtering resultCThe best value of (1)
Figure BDA0001359618900000212
(2) Fixing
Figure BDA0001359618900000213
bq0.9999, decreasing b graduallyQDetermining b from the filtering resultQThe best value of (1)
Figure BDA0001359618900000214
(3) Fixing
Figure BDA0001359618900000215
Gradually decrease bqDetermining b from the filtering resultqThe best value of (1)
Figure BDA0001359618900000216
System accuracy as bQ、bC、bqThe variation trend graph is shown in fig. 4, from which it can be determined that the optimal forgetting factor combination is: bq＝0.9993、bQ＝0.9997、bC＝0.96。
On the basis of analyzing the 'deficit rank' problem existing when only the inter-satellite distance measurement is utilized, Lagrange satellites are introduced to provide an absolute reference; and then, under a centralized structure, an adaptive EKF estimation method is adopted to realize the autonomous orbit determination of the earth-Lagrange joint constellation.
Under the condition that the inter-satellite range error is zero mean and the standard deviation is 5m, the following simulation results are obtained:
due to three in the constellationThe orbit determination accuracy of MEO satellites is similar, so that taking the satellite E1 as an example, the errors of the three-axis orbit determination distances of x, y and z are respectively 7.78m, 7.60m and 10.06m, the error curve is shown in FIG. 5, the distance estimation error is 14.81m, and the errors of the three-axis orbit determination speeds of x, y and z are respectively 1.43 × 10-3m/s，3.92×10-3m/s and 4.74 × 10-3m/s, speed estimation error 6.31 × 10-3m/s. similarly, the orbit determination accuracy of two Lagranges is similar, taking the satellite L1 as an example, the orbit determination distance errors of the three axes of x, y and z are respectively 2.66m, 3.60m and 3.39m, the distance estimation error is 5.62m, and the orbit determination speed errors of the three axes of x, y and z are respectively 2.11 × 10-5m/s，1.64×10-5m/s and 1.48 × 10-5m/s, velocity estimation error 3.06 × 10-5m/s。
The effectiveness of the method of the invention is verified.

Claims (4)

1. An earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging comprises the following steps:
the method comprises the following steps: under a centralized structure, establishing a state equation of an earth-Lagrange combined constellation autonomous orbit determination system;
the method specifically comprises the following steps:
(1) establishing earth navigation satellite orbit dynamics model
Under the earth-centered rectangular inertial coordinate system, the state equation for orbit determination by combining the earth satellite i to be detected in the constellation is as follows:
Figure FDA0002569392000000011
wherein, the state vector corresponding to the earth satellite i to be measured is
Figure FDA0002569392000000012
[xEi,yEi,zEi]TAnd
Figure FDA0002569392000000013
respectively, a position vector and a velocity vector of the earth satellite;FEas a function of the state of the system, WEi(t) is the system noise, satisfying the mean qEi(t) variance QEi(t) white gaussian noise; the formula can be written in further detail as:
Figure FDA0002569392000000014
wherein the content of the first and second substances,
Figure FDA0002569392000000015
and
Figure FDA0002569392000000016
is the noise component of x, y and z three-axis velocity,
Figure FDA0002569392000000017
and
Figure FDA0002569392000000018
is the acceleration equivalent noise component of the x, y and z three axes, including unmodeled perturbation term and noise term;
Figure FDA0002569392000000019
and
Figure FDA00025693920000000110
a modeled velocity vector and an acceleration vector corresponding to the earth satellite i, respectively, wherein the acceleration mainly takes into account the earth central gravitational acceleration a0,Ei、J2Acceleration of perturbation term
Figure FDA00025693920000000111
Acceleration a of sun-moon attractionMS,EiThe specific expression is as follows:
1) acceleration a of earth's center gravity of satellite i0,EiSatisfies the following conditions:
Figure FDA00025693920000000112
in the formula (I), the compound is shown in the specification,
Figure FDA00025693920000000113
mu is the product of the earth mass and the gravity constant G, namely the gravity constant of the earth;
2)J2acceleration of perturbation term
Figure FDA0002569392000000021
Satisfies the following conditions:
Figure FDA0002569392000000022
in the formula, ReIs the equatorial radius of the earth, J2Second order band harmonic coefficients;
perturbation acceleration a caused by sun-moon gravitationMS,EiSatisfies the following conditions:
Figure FDA0002569392000000023
in the formula (x)S,yS,zS) And (x)M,yM,zM) Respectively are three-dimensional position coordinates of the sun and the moon under a geocentric rectangular inertial coordinate system,
Figure FDA0002569392000000024
and
Figure FDA0002569392000000025
the distance from satellite i to the sun and moon respectively,
Figure FDA0002569392000000026
and
Figure FDA0002569392000000027
distance from the earth's center to the sun and moon, mu, respectivelySAnd muMAre respectively the sun leadingA force constant and a lunar gravity constant;
(2) lagrange navigation satellite orbit dynamics model
Establishing a state equation for orbit determination of the lunar satellite k to be measured in the Lagrange constellation under a centroid convergence coordinate system:
Figure FDA0002569392000000028
in the formula, the corresponding state vector of Lagrange lunar satellite k to be measured is
Figure FDA0002569392000000029
[xLk,yLk,zLk]TAnd
Figure FDA00025693920000000210
the dimensionless position vector and the velocity vector of the star respectively have characteristic quality, characteristic length and characteristic time shown in formula (7), wherein MEAnd MMThe earth mass and the moon mass are respectively, and L is the earth-moon distance; fLAs a function of the state of the system, WLk(t) is the system noise, satisfying the mean qLk(t) variance QLk(t) white gaussian noise;
Figure FDA0002569392000000031
equation (6) is written as:
Figure FDA0002569392000000032
wherein the content of the first and second substances,
Figure FDA0002569392000000033
for the modeled velocity vector, μ, corresponding to Lagrange satellite k0＝MM/(MM+ME)，ΔE,Lk、ΔM,LkThe earth center distance and the moon center distance of the satellite respectively,
Figure FDA0002569392000000034
and
Figure FDA0002569392000000035
is the noise component of x, y and z three-axis velocity,
Figure FDA0002569392000000036
and
Figure FDA0002569392000000037
is the acceleration equivalent noise component of the x, y and z three axes, including unmodeled perturbation term and noise term;
(3) state equation of earth-Lagrange combined constellation autonomous orbit determination system
Assuming that the whole autonomous navigation system comprises m earth satellites and n Lagrange satellites, establishing a constellation orbit determination state equation of a centralized structure according to a state equation for determining orbit of an earth satellite i to be detected and a state equation for determining orbit of a lunar satellite k to be detected as follows:
Figure FDA0002569392000000038
namely:
Figure FDA0002569392000000039
in the formula (I), the compound is shown in the specification,
Figure FDA00025693920000000310
f is a state vector, F is a state function vector of the whole constellation orbit determination system, w (t) represents system noise with a mean value q (t) and a variance q (t);
step two: establishing a measurement equation of an earth-Lagrange combined constellation autonomous orbit determination system;
step three: under a centralized structure, determining a filtering method for realizing orbit parameter estimation according to the established earth-Lagrange combined constellation autonomous orbit determination system model based on inter-satellite distance measurement;
step four: under a centralized structure, earth-Lagrange combined constellation autonomous orbit determination based on a selected filtering algorithm.
2. The earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging as claimed in claim 1, wherein the second step specifically comprises:
(1) inter-satellite distance measurement model between earth navigation satellites
For the earth navigation satellite, the observed distance between the satellite i and the satellite j is recorded as rhoEi,EjThen, there are:
ρEi,Ej＝hE(XEi,XEj)+vEi,Ej，i≠j (11)
wherein
Figure FDA0002569392000000041
vEi,EjEquivalent measurement noise between the satellite i and the satellite j, including modeling error and measurement noise;
(2) inter-satellite distance measurement model between Lagrange navigation satellites
Let the observation of the distance between Lagrange satellites i and j be ρLi,LjThe method comprises the following steps:
ρLi,Lj＝hL(XLi,XLj)+vLi,Lj，i≠j (12)
wherein
Figure FDA0002569392000000042
vLi,LjEquivalent measurement noise between Lagrange satellites i and j, including modeling error and measurement noise;
(3) measurement model between earth navigation satellite and Lagrange satellite
When a GNSS satellite i and a Lagrange satellite k are represented in the same coordinate system, there are:
Figure FDA0002569392000000043
wherein the content of the first and second substances,
Figure FDA0002569392000000044
is the coordinate of Lagrange satellite in the inertial coordinate system of earth center, which is the dimensionless coordinate X of Lagrange satellite in the coordinate system of center-of-mass convergenceLkThe conversion relationship is:
Figure FDA0002569392000000051
wherein, R represents a rotation matrix from a centroid convergence coordinate system to a geocentric inertia coordinate system, and the expression is as follows:
Figure FDA0002569392000000052
wherein u isM＝ωM+θM，iM、ΩM、ωMAnd thetaMRespectively representing the orbit inclination angle, the ascent intersection declination, the argument of the perigee and the true perigee angle of the orbit of the moon around the ground;
(4) measurement equation of earth-Lagrange combined constellation autonomous orbit determination system
According to the above measurement equations (11), (12) and (13), the constellation orbit determination measurement equation for establishing the centralized structure is as follows:
Z(t)＝h[X(t),t]+V(t) (16)
namely:
Figure FDA0002569392000000053
where V (t) is the equivalent measurement noise with a mean of zero and a variance of R (t).
3. The earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging according to claim 1, wherein the third step employs an EKF method or an SPKF method, specifically:
(1) when the state estimation is carried out on the discretized nonlinear system by adopting the self-adaptive EKF method, the specific steps are as follows:
1) time updating
Calculating a one-step prediction state:
Figure FDA0002569392000000061
in the formula (I), the compound is shown in the specification,
Figure FDA0002569392000000062
is an estimate of the mean value of the system noise at time k,
Figure FDA0002569392000000063
for one-step recursion of states using kinetic equations, one method is by calculation using equation (19) after discretization and linearization, where
Figure FDA0002569392000000064
For the state estimation at the time instant k,
Figure FDA0002569392000000065
estimating the one-step prediction state at the moment (k +1), wherein delta t is an integral step length; the other method is obtained by directly integrating and calculating the numerical value of the continuous state equation (9) by a 4-order Runge-Kutta method;
Figure FDA0002569392000000066
calculating a state transition matrix:
Figure FDA0002569392000000067
in the formula phik+1/kIs the state transition matrix from time k to time (k + 1);
calculating a one-step prediction error covariance matrix:
Figure FDA0002569392000000068
in the formula, Pk+1/kError covariance matrix, P, for one-step prediction of stateskThe error covariance matrix of the state is estimated for time k,
Figure FDA0002569392000000069
the estimated value of the system noise variance matrix at the moment k is obtained;
in the initial time 0-T of filtering, adopting the prior statistical information of noise, and updating the time as follows:
Figure FDA00025693920000000610
Figure FDA00025693920000000611
2) and (3) measurement updating:
calculating a gain matrix:
Figure FDA00025693920000000612
in the formula, Kk+1Is the gain matrix at time (k +1),
Figure FDA00025693920000000613
to predict state estimation based on one step
Figure FDA00025693920000000614
Jacobian matrix, R, of calculated observation vectors hk+1Measuring a noise covariance matrix for the (k +1) time;
calculating an innovation vector:
Figure FDA00025693920000000615
in the formulak+1Is (k +1)An innovation vector of time, Zk+1Is the measurement vector at time (k + 1);
calculating a state estimation value:
Figure FDA0002569392000000071
calculating a state estimation error covariance matrix:
Figure FDA0002569392000000072
3) calculating an estimated value of the statistical characteristic of the system noise:
Figure FDA0002569392000000073
Figure FDA0002569392000000074
Figure FDA0002569392000000075
wherein the content of the first and second substances,
Figure FDA0002569392000000076
Figure FDA0002569392000000077
Figure FDA0002569392000000078
in the formula 0 < bq＜1，0＜bQ＜1，0＜bC< 1 are forgetting factors for calculating Q, Q, C, respectively, and as k increases,
Figure FDA0002569392000000079
Figure FDA00025693920000000710
and
Figure FDA00025693920000000711
respectively tend to (1-b)q)、(1-bQ) And (1-b)C)；
When it occurs
Figure FDA00025693920000000712
In the case of negative determination, in the case of equation (29)
Figure FDA00025693920000000713
The diagonal elements of the terms take absolute values, and the non-diagonal elements take 0;
(2) when the state estimation is carried out on the continuous nonlinear system by adopting the self-adaptive SPKF method, the method comprises the following specific steps:
1) selection of Sigma Point sampling strategy
Calculating a weight coefficient W of a Sigma point Kalman filteri m、Wi cWhere i is 0,1, …, L is the number of sampling points, Wi mFor mean value estimation, Wi cFor variance estimation;
2) and (3) time updating:
according to the sampling strategy selected in 1) and the state estimation value at the k moment
Figure FDA00025693920000000714
Sum estimation error covariance PkSigma point set chi for calculating k timei,k；
Calculate one-step predicted values for Sigma points:
χi,k+1/k＝F[χi,k,k]，i＝0,1,…,L (34)
wherein, χi,kIs the value of the ith Sigma point at time k, χi,k+1/kA one-step predicted value of the ith Sigma point at the moment of k +1 is obtained;
calculating a one-step prediction state:
Figure FDA0002569392000000081
in the formula (I), the compound is shown in the specification,
Figure FDA0002569392000000082
is an estimate of the mean value of the system noise at time k,
Figure FDA0002569392000000083
one-step prediction state estimation for time (k + 1);
calculating a state one-step prediction error covariance matrix:
Figure FDA0002569392000000084
in the formula, Pk+1/kFor the error covariance matrix of the one-step predicted state,
Figure FDA0002569392000000085
the estimated value of the system noise variance matrix at the moment k is obtained;
in the initial time 0-T of filtering, adopting the prior statistical information of noise, and updating the time as follows:
Figure FDA0002569392000000086
Figure FDA0002569392000000087
3) and (3) measurement updating:
according to the sampling strategy selected in 1) and
Figure FDA0002569392000000088
and Pk+1/kUpdating a Sigma point set;
estimating the measured value:
γi,k+1/k＝h(χi,k+1/k) (39)
Figure FDA0002569392000000089
wherein: gamma rayi,k+1/kIs an estimate of the measurement vector at time (k +1) calculated from the ith Sigma point,
Figure FDA00025693920000000810
a weighted estimate of the measurement vector at time (k + 1);
calculating a gain matrix:
Figure FDA00025693920000000811
Figure FDA00025693920000000812
Figure FDA0002569392000000091
wherein R isk+1Measuring the noise covariance matrix for the (k +1) timey,k+1/kMeasuring the covariance of the vector estimates for the (k +1) th moment, Pxy,k+1/kIs the covariance of the (K +1) th measured vector estimate and the state one-step predictor, Kk+1A gain matrix for time (k + 1);
calculating an innovation vector:
Figure FDA0002569392000000092
in the formulak+1Is the innovation vector at time (k +1), Zk+1Is the measurement vector at time (k + 1);
calculating a state estimate at time (k +1)
Figure FDA0002569392000000093
Figure FDA0002569392000000094
Computing a state estimation error covariance matrix P at time (k +1)k+1：
Figure FDA0002569392000000095
4) Calculating an estimated value of the statistical characteristic of the system noise:
Figure FDA0002569392000000096
Figure FDA0002569392000000097
Figure FDA0002569392000000098
wherein the content of the first and second substances,
Figure FDA0002569392000000099
Figure FDA00025693920000000910
Figure FDA00025693920000000911
in the formula 0 < bq＜1，0＜bQ＜1，0＜bCLess than 1 is used for calculating forgetting factors of Q, Q and C respectively;
when it occurs
Figure FDA00025693920000000912
In the case of negative determination, in the case of equation (48)
Figure FDA00025693920000000913
The diagonal elements of the term take the absolute value and the off-diagonal elements take 0.
4. The earth-Lagrange joint constellation autonomous orbit determination method based on inter-satellite ranging as claimed in claim 1, wherein the fourth step specifically is:
(1) initializing data;
the initialized parameters include: k is 0, filtering step h, noise switching parameter T, initial state
Figure FDA0002569392000000101
Initial error covariance matrix P0Initial system noise mean
Figure FDA0002569392000000102
Initial system noise variance matrix
Figure FDA0002569392000000103
Measuring noise variance matrix R and forgetting factor bq、bQ、bC；
The noise switching parameter T is a time node for switching noise statistical characteristics in the time updating algorithm: when k is less than T, using the prior statistical property of the noise to update time; when k is larger than or equal to T, updating time by using the estimated noise statistical characteristics; initial state
Figure FDA0002569392000000104
Initial error covariance matrix P0Initial system noise variance matrix
Figure FDA0002569392000000105
Determining according to the requirement; initial system noise mean
Figure FDA0002569392000000106
Taking zero; measuring noise variance matrix R according to measuring apparatusPerformance determination; forgetting factor bq、bQThe value is between 0.99 and 1.0, bCThe value is between 0.9 and 0.99;
(2) if the self-adaptive EKF algorithm is adopted, directly turning to (4); if the self-adaptive SPKF algorithm is adopted, turning to (3);
(3) calculating weight coefficient W according to selected Sigma point sampling strategyi m、Wi c；
(4) If k is less than T, updating time by adopting the prior statistical property of the noise, and turning to (5); otherwise, updating time by adopting the estimated noise statistical characteristics, and turning to the step (6);
(5) time updating using a priori statistical properties of the noise; turning to (7) after completion;
state estimation using time k
Figure FDA0002569392000000107
And error covariance matrix PkSequentially calculating the state one-step predicted value as the state initial value and the error covariance matrix at the k moment
Figure FDA0002569392000000108
And the prediction error covariance matrix Pk+1/kTime updating of the self-adaptive Kalman filtering estimation method is realized; when the self-adaptive EKF algorithm is adopted, the calculation formulas are respectively shown as a formula (22) and a formula (23); when the self-adaptive Sigma point Kalman filtering algorithm is adopted, the calculation formulas are respectively shown as a formula (37) and a formula (38);
(6) time updating using the estimated noise statistics;
state estimation using time k
Figure FDA0002569392000000109
And error covariance matrix PkSequentially calculating the state one-step predicted value as the state initial value and the error covariance matrix at the k moment
Figure FDA00025693920000001010
And the prediction error covariance matrix Pk+1/k"ShiUpdating the time of the current self-adaptive Kalman filtering estimation method; when the self-adaptive EKF algorithm is adopted, the calculation formulas are respectively an expression (18) and an expression (21); when the self-adaptive SPKF algorithm is adopted, the calculation formulas are respectively shown as the formula (35) and the formula (36);
(7) updating the measurement;
sequentially calculating a gain matrix K at (K +1) timek+1Innovation matrix vk+1State estimation value
Figure FDA00025693920000001011
And an error covariance matrix Pk+1Realizing the measurement update of the self-adaptive Kalman filtering; when the self-adaptive EKF algorithm is adopted, the calculation formulas are respectively shown as formulas (24) to (27); when the self-adaptive SPKF algorithm is adopted, the calculation formulas are respectively expressed as formulas (43) to (46);
(8) estimating the statistical characteristics of system noise;
calculating mean value estimated value of system noise in turn
Figure FDA0002569392000000111
Sum covariance estimate
Figure FDA0002569392000000112
When the adaptive EKF algorithm is adopted, the calculation formulas are respectively shown as a formula (28) and a formula (29); when the self-adaptive SPKF algorithm is adopted, the calculation formulas are respectively shown as a formula (47) and a formula (48);
(9) judging whether the filtering process is finished or not; filtering is continued as necessary, k is k +1, and (4) is returned.
CN201710611461.6A 2017-07-25 2017-07-25 Earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging Active CN107421550B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710611461.6A CN107421550B (en) 2017-07-25 2017-07-25 Earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710611461.6A CN107421550B (en) 2017-07-25 2017-07-25 Earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging

Publications (2)

Publication Number Publication Date
CN107421550A CN107421550A (en) 2017-12-01
CN107421550B true CN107421550B (en) 2020-08-28

Family

ID=60431049

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710611461.6A Active CN107421550B (en) 2017-07-25 2017-07-25 Earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging

Country Status (1)

Country Link
CN (1) CN107421550B (en)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108287334A (en) * 2018-02-06 2018-07-17 西安四方星途测控技术有限公司 A kind of Spin Satellite Attitude method of estimation and system based on rcs measurement data
CN109031349B (en) * 2018-04-20 2022-04-08 南京航空航天大学 Intelligent autonomous operation system of GEO satellite
CN109752006B (en) * 2018-11-23 2022-09-09 中国西安卫星测控中心 Method for using incomplete external measurement data in real-time filtering
CN109682383B (en) * 2018-11-23 2022-11-04 中国西安卫星测控中心 Real-time filtering positioning method for measuring distance and data by using deep space three-way
CN109933847B (en) * 2019-01-30 2022-09-16 中国人民解放军战略支援部队信息工程大学 Improved active segment trajectory estimation algorithm
CN109917431B (en) * 2019-04-02 2021-03-23 中国科学院空间应用工程与技术中心 Method for realizing GNSS satellite autonomous navigation by utilizing DRO orbit and inter-satellite measurement
CN109975839B (en) * 2019-04-10 2023-04-07 华砺智行(武汉)科技有限公司 Joint filtering optimization method for vehicle satellite positioning data
CN110793528B (en) * 2019-09-27 2021-04-13 西安空间无线电技术研究所 Low-orbit satellite-based anchoring-based Beidou navigation constellation autonomous orbit determination method
CN115390109A (en) * 2020-04-30 2022-11-25 中国科学院微小卫星创新研究院 Beidou satellite centralized constellation autonomous navigation method
CN111947667B (en) * 2020-06-24 2022-08-12 火眼位置数智科技服务有限公司 Low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination
CN112504281B (en) * 2020-11-16 2023-06-27 中国人民解放军63921部队 Spacecraft orbit determination method based on Beidou inter-satellite unidirectional link
CN113008243B (en) * 2021-02-23 2024-03-19 重庆两江卫星移动通信有限公司 Autonomous orbit determination method and system for low orbit satellite constellation
CN113204917B (en) * 2021-04-25 2021-12-07 中国科学院国家空间科学中心 Space-based optical angle measurement arc section initial orbit determination method for GEO target and correlation method

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9091552B2 (en) * 2011-10-25 2015-07-28 The Boeing Company Combined location and attitude determination system and methods
WO2014172675A1 (en) * 2013-04-18 2014-10-23 California Institute Of Technology Real-time and post-processed orbit determination and positioning
US9939260B2 (en) * 2014-08-28 2018-04-10 The Boeing Company Satellite transfer orbit search methods
CN105716612B (en) * 2016-02-29 2017-05-10 武汉大学 Method for designing strapdown inertial navigation system simulator
CN106338753B (en) * 2016-09-22 2019-03-12 北京航空航天大学 One kind being based on earth station/inter-satellite link/GNSS combined measurement geostationary orbit constellation orbit determination method
CN106885577B (en) * 2017-01-24 2020-01-21 南京航空航天大学 Autonomous orbit determination method for Lagrange navigation satellite

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"导航卫星自主导航关键技术研究";肖寅;《中国博士学位论文全文数据库（电子期刊）》;20160115;正文第24-29页 *

Also Published As

Publication number Publication date
CN107421550A (en) 2017-12-01

Similar Documents

Publication Publication Date Title
CN107421550B (en) Earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging
Xiangyu et al. An autonomous optical navigation and guidance for soft landing on asteroids
CN109255096B (en) Geosynchronous satellite orbit uncertain evolution method based on differential algebra
CN107797130A (en) Low orbit spacecraft multiple spot multi-parameter track upstream data computational methods
Takahashi et al. Spin state and moment of inertia characterization of 4179 Toutatis
CN111552003B (en) Asteroid gravitational field full-autonomous measurement system and method based on ball satellite formation
CN112161632B (en) Satellite formation initial positioning method based on relative position vector measurement
CN111522037A (en) Autonomous navigation method and navigation system for constellation co-orbital plane satellite
Hesar et al. Lunar far side surface navigation using linked autonomous interplanetary satellite orbit navigation (LiAISON)
CN104048664A (en) Autonomous orbit determination method of navigation satellite constellation
CN104848862A (en) Precise and synchronous positioning and time-keeping method and system of Mars orbiting detector
CN109752005A (en) A kind of Method of Spacecraft Initial Orbit Determination based on precise orbit model
Lee et al. Vision-based relative state estimation using the unscented Kalman filter
Ke et al. Pico-satellite autonomous navigation with magnetometer and sun sensor data
Qian et al. Sun–Earth–Moon autonomous orbit determination for quasi-periodic orbit about the translunar libration point and its observability analysis
Gui et al. A time delay/star angle integrated navigation method based on the event-triggered implicit unscented Kalman filter
CN113589832B (en) Constellation rapid design method for stable observation coverage of ground surface fixed area target
Babcock CubeSat attitude determination via Kalman filtering of magnetometer and solar cell data
Lou et al. A consider unscented particle filter with genetic algorithm for UAV multi-source integrated navigation
Baroni Attitude determination by unscented Kalman filter and solar panels as sun sensor
Liu et al. Guidance and control technology of spacecraft on elliptical orbit
Nordkvist et al. Attitude feedback tracking with optimal attitude state estimation
Liuqing et al. An autonomous navigation study of Walker constellation based on reference satellite and inter-satellite distance measurement
CN111854765B (en) Medium-orbit navigation satellite orbit long-term forecasting method
Olson Sequential estimation methods for small body optical navigation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant