# Autonomous Orbit Determination for Lagrangian Navigation Satellite Based on Neural Network Based State Observer

Authors: Youtao Gao; Tanran Zhao; Bingyu Jin; Junkang Chen; Bo Xu
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

In order to improve the accuracy of the dynamical model used in the orbit determination of the Lagrangian navigation satellites, the nonlinear perturbations acting on Lagrangian navigation satellites are estimated by a neural network. A neural network based state observer is applied to autonomously determine the orbits of Lagrangian navigation satellites using only satellite-to-satellite range. This autonomous orbit determination method does not require linearizing the dynamical mode. There is no need to calculate the transition matrix. It is proved that three satellite-to-satellite ranges are needed using this method; therefore, the navigation constellation should include four Lagrangian navigation satellites at least. Four satellites orbiting on the collinear libration orbits are chosen to construct a constellation which is used to demonstrate the utility of this method. Simulation results illustrate that the stable error of autonomous orbit determination is about 10 m. The perturbation can be estimated by the neural network.


---

## Body

## 1. Introduction

The navigation of deep space probes is one of the main problems that restrict the deep space exploration. Generally, navigation support for deep space probes is primarily provided by the NASA’s Deep Space Network (DSN). In order to improve the navigation performance by efficiently determining angular position of interplanetary spacecraft, the delta differential one-way tanging technique is also employed by some missions. Besides that, some autonomous navigation strategies are also proposed to support future deep space exploration missions. Autonomous navigation is important for deep space probes to deal with communication delay as well as reducing the dependency on ground stations. Several autonomous navigation methods have been proposed for deep space probes. As early as 1968, the sextant had been used for autonomous navigation in Apollo program [1]. In 1999, Deep Space 1 achieved autonomous orbit determination by tracking small celestial bodies with an optical sensor [2]. The comet probe “Deep Impact” carried out its navigation automatically based on an optical navigation system [3]. The rotation period of X-ray pulsar is extremely stable; therefore, time and the location of spacecraft have been proposed to be determined by tracking several X-ray pulsars in [4, 5]. The Global Positioning System (GPS) can support the navigation of deep space probes when they are orbiting in low Earth orbits and medium Earth orbits. For deep space transfer orbits and deep space target orbits, due to limited visibility, extremely low signal-to-noise ratio, and poor relative geometry among sources and users, GPS is not adequate. In 2005, Hill suggested placing ​​navigation constellation on the periodic orbits in the vicinity of libration points of the Earth-Moon system to support ​​deep space navigation [6]. Similar to GPS, a high-precision satellite navigation constellation which consists of libration point satellites in the Earth-Moon system is introduced to provide navigation information for deep space probes, which can be called, accordingly, the Lagrangian point satellite navigation system. The satellites which construct the Lagrangian point satellite navigation system are called Lagrangian navigation satellites. Zhang and Xu analyzed the architecture and navigation performance of the Lagrangian point satellite navigation system [7–9]. The Lagrangian navigation constellation is introduced to navigate the deep space probes autonomously. Hence the navigation constellation itself should have the ability of autonomous orbit determination (AOD). In [10], the feasibility of AOD for satellites in quasiperiodic orbits about the Earth-Moon libration point was verified. Based on circular restricted three-body problem (CRTBP), Du et al. researched the autonomous orbit determination method of satellites in halo orbits, and only satellite-to-satellite range was used as observation [11]. For the Earth navigation satellite constellation, there is a rank deficiency problem when only satellite-to-satellite range is used to determine the orbit. However, the rank deficiency problem does not exist for the Lagrangian navigation satellites because of the special dynamics near the libration points [6]. Thus, the Lagrangian navigation satellites in the navigation constellation can autonomously determine their orbits using only satellite-to-satellite range. In [12], Gao et al. discussed the algorithm of autonomous orbit determination using only the satellite-to-satellite range measurement for Lagrangian navigation constellation. The current studies about AOD of Lagrangian navigation satellites are under the CRTBP model. However, the motion of the Moon around the Earth is eccentric. Therefore, the elliptic restricted three-body problem (ERTBP) is more accurate to describe the motion of the Lagrangian navigation satellite [13]. The ERTBP has been discussed in detail in [14–16]. Our purpose here is to extend the applicability of ERTBP to the study of the AOD of Lagrangian navigation satellites. In this study, we consider ERTBP with perturbation. The perturbation is estimated using a neural network. Meanwhile an observer is designed to determine the orbit of Lagrangian satellite. We reference the design of a reduced-order modified state observer which is introduced in [17]. However, in [17], the authors assume that the position can be measured. In our study, we improve the observer which can estimate the state of Lagrangian satellite with only satellite-to-satellite range.First the dynamical model of ERTBP with perturbation is introduced. Then we design a neural network based state observer to determine the orbit of Lagrangian satellites. Afterwards the stability of the observer is proved. Finally four Lagrangian satellites are chosen to validate the effectiveness of this AOD method.


## 2. Elliptic Restricted Three-Body Problem with Perturbation

As shown in Figure1, P1 and P2 are the primaries in the three-body system. P1 and P2 are in elliptical orbits. P is the third body which is vanishingly small compared to the two primaries. Similar to the CRTBP, we study the motion of P in a rotating coordinates.Figure 1
The barycenter synodic coordinatesC-XYZ and the libration points.LetC-XYZ represent the barycentric synodic coordinate depicted in Figure 1; the x-axis of this frame is along the radius vector, which connects the primaries, positive in the direction pointing from P1 to P2. The y-axis of this frame is perpendicular to the x-axis, positive in the direction of the motion of P2. The XY frame rotates with an angular velocity equal to the instant motion of P2 with respect to P1. The z-axis is perpendicular to the orbital plane of the primaries.The equations of motion in ERTBP are dimensionless. The dimensionless units are as follows:(1)M=m1+m2L=a1-e21+ecos⁡fT=L3Gm1+m21/2=1+ecos⁡ff˙,where m1 and m2 are the masses of the two primaries. a, e, respectively, refer to the semimajor axis and eccentricity of the two primaries’ elliptical relative revolving orbit. f is the true anomaly of the secondary on the elliptic orbit. Then the motion of the Lagrangian satellite in the barycentric synodic coordinate system is governed by the following [18]:(2)X′′-2Y′=11+ecos⁡f∂Ω∂XY′′+2X′=11+ecos⁡f∂Ω∂YZ′′=11+ecos⁡f∂Ω∂Z.f is taken as the time-like independent variable. The first and second derivatives of the coordinate with respect to f are calculated as(3)X′=dXdf,X′′=d2Xdf2.Ω is the pseudo-potential function of the three-body problem; it is described as follows:(4)Ω=12X2+Y2+Z2+μ1-μ+1-μr1+μr2,where μ=m2/m1+m2, r1=X+μ2+Y2+Z2, r2=X-1+μ2+Y2+Z2.In this study, only the motions around the collinear libration pointL1 or L2 are investigated. Therefore the origin of the coordinate can be moved from the barycenter of system to the interested libration point for convenient analysis. As shown in Figure 2, the instantaneous distance between the libration point and its closest primary which is denoted as γ is chosen as the new length unit. The new reference frame is defined as the L1- or L2-centered synodic reference frame. The directions of the axes of this new coordinate system are coincident with the barycentric synodic reference frame.Figure 2
TheL1-centered synodic reference frame.Letx,y,z,x′,y′,z′ represent the state variables in L1- or L2-centered synodic reference frame. The transformation of coordinates between the barycentric synodic frame and the L1- or L2-centered synodic system is represented as follows [19]:(5)X=γx∓1+1-μY=γyZ=γz,where the upper sign refers to the L1 case and the lower one refers to the L2 case. The linearized equations of motion in the L1- or L2-centered synodic system can be formulated as follows [20]:(6)x′′-2y′-1+2c2x=∑i≥1-eicosif1+2c2xy′′+2x′-1-c2y=∑i≥1-eicosif1-c2yz′′+c2z=∑i≥1-eicosif1-c2z,where c2μ=1/γj3μ+1-μγj2/1∓γj3 and γj (j=1,2) is the instantaneous distance between Lj and its closest primary.We define a new variableX-=[x,y,z,x′,y′,z′]T; then (6) can be written as (7)X¯′=AX-,where(8)A=0001000000100000011+2c21+∑i≥1-eicosif0002001-c2∑i≥1-eicosif+10-20000000-c2+∑i≥1-eicosif1-c2In addition to the gravitational force from the two primaries, several perturbations can affect the motion of Lagrangian satellite. When these perturbations are considered, the motion equation of the Lagrangian satellite can be described as(9)X-′=AX-+BgX-,where(10)gX-=gxgygzis the nonlinear perturbing force on per unit mass. And (11)B=000000000100010001.


## 3. The Observation

In this paper, satellite-to-satellite range is the observation. The relationship between the observation and the state is(12)Y-=ρX-,t=x-x22+y-y22+z-z22,where xyzT and x2y2z2T are the position variable of Lagrangian satellite 1 and Lagrangian satellite 2, respectively, in L1- or L2-centered synodic reference frame.Expanding (12) in a Taylor series with respect to estimated state yields (13)Y-=x-x22+y-y22+z-z22=Y-^+∂Y-∂xx=x^x-x^+∂Y-∂yy=y^y-y^+∂Y-∂zz=z^z-z^+⋯.We get the linearized relation between observation and states by neglecting the higher order terms. (14)Y-=Y-^+∂Y-∂xx=x^x-x^+∂Y-∂y1y=y^y-y^+∂Y-∂zz=z^z-z^=Y-^+x^-x2ρx-x^+y^-y2ρy-y^+z^-z2ρz-z^.Now define the state estimate error as(15)X~=X--X-^.Further, the residual error of observation is defined as (16)Y~=Y--Y-^.Then we get (17)Y~=CX~,where (18)C=x^-x2ρy^-y2ρz^-z2ρ000.Usually we only can get the estimated value of the position of Lagrangian satellite 2; thus the matrix C can be calculated by using the following:(19)C=x^-x^2ρy^-y^2ρz^-z^2ρ000.


## 4. Design of the Adaptive Observer Based on Neural Network

The perturbationg(X-) in (9) can be estimated using a neural network over a compact set DX- as follows [17]:(20)gX-=WTϕX-+εX-ε≤ε∗,X-∈DX-,where W is a matrix of unknown neural network weights, ϕX- is a known vector of bounded basis functions, and ε∗ is a uniform bound on the approximation error.In order to estimate the state of Lagrangian satellite using the satellite-to-satellite range, an observer is designed as(21)X-^′=AX-^+Bg^X-^-υf+KY~,where υf is a robust term that will be determined later. f is the time-like independent variable. K is a user-specified gain matrix. g^(X-^) is the estimated uncertainty vector which is given by(22)g^X-^=W^TϕX-^.A weight estimate update W^˙ is designed as(23)W^˙=FϕX-^Y~T-σW^,where F is a positive-definite symmetric matrices that will be determined later. And(24)υf=-DY~Y~-εmaxY~.ϕX-^ is calculated as follows [17]:(25)ϕX-^=1x^R⊕y^R⊕z^R⊕x^y^R⊕2x^z^R⊕2y^z^R⊕2.Now we have completed the design of the observer. However when this method is used to determine the orbits of the Lagrangian satellite in a navigation constellation, the number of the Lagrangian navigation satellites must satisfy a constraint condition. In the following, we will derive the minimum number of the Lagrangian navigation satellites which is needed in this method.In order to prove that the designed observer in (21) can estimate the state of Lagrangian satellites accurately, we define the errors in the neural network weights and the basis vector ϕ as follows:(26)W~=W-W^ϕ~=ϕX--ϕX-^.And a candidate Lyapunov functions is chosen as [17](27)V=12X~TPX~+12trW~TF-1W~,where P and F are positive-definite symmetric matrices. By differentiating (27), the Lyapunov function’s derivative is illustrated as(28)V˙=X~TPX~˙+trW~TF-1W~˙=X~TPA-KCX~+BWTϕX--W^TϕX-^+υf+BεX-+trW~TF-1W~˙.We define (29)Acl=A-KC.Then it yields (30)V˙=X~TPAclX~+BWTϕX--WTϕX-^+WTϕX-^-W^TϕX-^+υf+BεX-+trW~TF-1W~˙=X~TPAclX~+BWTϕ~+W~TϕX-^+υf+εX-+trW~TF-1W~˙=12X~TPAcl+AclTPX~+X~TPBWTϕ~+W~TϕX-^+υf+εX-+trW~TF-1W~˙.If PB=CT, it will yield (31)V˙=12X~TPAcl+AclTPX~+CX~TWTϕ~+W~TϕX-^+υf+εX-+trW~TF-1W~˙=12X~TPAcl+AclTPX~+trWTϕ~+W~TϕX-^+υf+εX-Y~T+W~TF-1W~˙.Here the Meyer-Kalman-Yakubovich (MKY) lemma [17, 21] is used to derive PB=CT. Noting that P~6×6, B~6×3, if PB=CT, C should satisfy C~6×3. It means that three satellite-to-satellite ranges should be provided. In other words, the navigation constellation should include four Lagrangian satellites which can generate three satellite-to-satellite ranges.If four Lagrangian satellites construct the navigation constellation,Ccan be derived as(32)C=x^-x^2ρ2y^-y^2ρ2z^-z^2ρ2000x^-x^3ρ3y^-y^3ρ3z^-z^3ρ3000x^-x^4ρ4y^-y^4ρ4z^-z^4ρ4000,where x^3y^3z^3T and x^4y^4z^4T are the estimated position variable of Lagrangian satellite 3 and Lagrangian satellite 4, respectively.Define matrixQ as(33)-Q=PAcl+AclTP,Q>0.The Lyapunov function’s derivative becomes(34)V˙=-12X~TQX~+trWTϕ~+W~TϕX-^+υf+εX-Y~T+W~TF-1W~˙.Since W~˙=-W^˙ and W^˙=FϕX-^Y~T-σW^, (35)V˙=-12X~TQX~+trWTϕ~+υf+εX-Y~T+W~TϕX-^Y~T-F-1FϕX-^Y~T-σW^=-12X~TQX~+trWTϕ~Y~T-trDY~Y~TY~-trεmaxY~Y~T+trεX-Y~T+trσW~TW^=-12X~TQX~+trWTϕ~Y~T-trDY~-trεmaxY~Y~T+trεX-Y~T+trσW~TW^.From [17], it is noted that(36)trWTϕ~Y~T≤αY~.We also can get (37)-12X~TQX~≤-λminQX~2≤-λminQY~2trεmaxY~Y~T=εmaxY~2trσW~TW^≤WmaxW^F-W^F2.·F is calculated as the Frobenius norm.Substituting (36)-(37) into (35), we can get (38)V˙≤-λminQ+εmaxY~2+α-D+ε∗Y~+σWmaxW^F-σW^F2≤-λminQ+εmaxY~2+α-D+ε∗Y~-σW^F-Wmax22+σWmax24.Thus, from (38), V˙<0 when(39)Y~≥α-D+ε∗+α-D+ε∗2+λminQ+εmaxσWmax22λminQ+εmax,or(40)W^F≥Wmax+Wmax2+α-D+ε∗2/σλminQ+εmax2.Note that parameters D, εmax, σ are user selected so they can be appropriately chosen to minimize the estimation error.


## 5. Simulation and Analysis

In order to illustrate the effectiveness of the method introduced in this paper, we design a constellation around the Earth-Moon L1 libration point which includes four Lagrangian satellites. The initial states of the four satellites are listed in Table1. The data are normalized. Quasiperiodic orbits are sensitive to initial states. Thus the initial states must be accurate enough. This is why the data in Table 1 retain many valid digits after the dot. We calculate quasiperiodic orbits with different initial states. The conditions in Table 1 may be different from the actual task situation. Therefore the initial conditions and measurement noise here are assumed conditions based on the characteristics of CRTBP.Table 1
Initial state of the Lagrangian satellites.
Satellites
S1
S2
S3
S4
x
L
-
0.009158653890369
-
0.008426021835493
-
0.011852455058632
-
0.012153563498421
y
L
0.057117380302627
0.063116578748419
0.018916778218648
-
0.000027358900568
z
L
-
0.005441690073741
-
0.009538268525430
0.000050795432543
-
0.001861700547480
x
˙
L
/
T
0.071242540915778
0.077724037910168
0.043377272215699
0.030872891241568
y
˙
L
/
T
0.025237253422237
0.044432703977482
0.002938265403217
-
0.000128391119819
z
˙
L
/
T
0.024768142966689
0.092125880718452
-
0.000042323140915
0.011290787920073Four quasiperiodic orbits are shown in Figure3. The ephemeris of the Moon and Earth is obtained by DE405. In order to illustrate the AOD error, the “true trajectory” of Lagrangian satellites should be provided. In this paper, the “true trajectory” is calculated by using the precise dynamical equations as follows:(41)R¨=-μ3R3R-∑i=1nμiΔiΔi3+RiRi3,where R is the position vector of the Lagrangian satellite in the J2000 coordinate. Δi is the position vector of the ith celestial with respect to the Moon. Ri is the position vector of the ith celestial with respect to the Earth. μi is the gravitational parameter of ith celestial. The quasiperiodic orbits are calculated by using multiple shooting method. The initial targets are obtained by using approximate analytic solutions. The initial targets are expressed in the barycentric synodic coordinate. Therefore, we should translate the initial targets into J2000 coordinate. Then the initial target is modified according to the task constraint, and the ideal quasiperiodic trajectory is obtained. Since the Lagrangian satellites in the Earth-Moon system is far away from the Earth, the perturbation caused by atmospheric resistance, temporal change in gravity, oblateness, and mass distribution is not considered. The main perturbations come from the gravity of the Sun and other planets. Therefore, the above equation is appropriate to describe the precise motion of the Lagrangian satellites.Figure 3
The quasiperiodic orbits around the L1 libration.The various tuning constants used in this scenario were chosen as follows:(42)K=1.5×10-30001.5×10-30001.5×10-36.5×10-40006.5×10-40006.5×10-4F=diag⁡1.5×10-5,3.2×10-5,1.5×10-5,3.2×10-5,2.25×10-5,1.5×10-5,1.5×10-5.The parameters are designed as εmax=1×10-12, D=1×10-9. The step of AOD is 1 hour. Three scenarios with different initial state errors and observation errors as shown in Table 2 are analyzed. The gravitational attraction from the Sun is the only perturbation considered in the simulation.Table 2
The simulation conditions.
Initial state error (m)
Observation error (m)
Scenario 1
1
1
Scenario 2
10
1
Scenario 3
100
1For scenario 1, we give the simulation results of two Lagrangian satellites S1 and S2. The AOD errors in 30 days are presented in Figures4(a) and 4(b). It can be seen that the AOD error of the both Lagrangian satellites is less than 10 m. The RMS of the AOD error of S1 in x-, y-, and z-axis are 2.3 m, 2.6 m, and 2.1 m, respectively. From Figures 5(a) and 5(b), we can see that this method can estimate the perturbation with small errors. The red curves are the theoretical calculating value of the Sun gravitational perturbation. The blue curves are the estimation of the Sun’s gravity. The RMS of the estimate error of perturbation of S1 in x-, y-, and z-axis are 6.1×10-6 m/s2, 5.9×10-6 m/s2, and 4.3×10-7 m/s2.Figure 4
(a)
AOD error of S1 in Scenario 1
(b)
AOD error of S2 in Scenario 1Figure 5
(a)
True perturbation and the estimated perturbation of S1 in Scenario 1
(b)
True perturbation and the estimated perturbation of S2 in Scenario 1In Figures6(a), 6(b), 7(a), and 7(b), we give the results for Scenario 2 and Scenario 3. As we can see, the maximum AOD errors increase with the initial state errors. All the stable AOD errors are less than 10 m. As shown in Figure 7(a), the AOD errors can converge to less than 10 m in 5 days, even the initial errors 100 m.Figure 6
(a)
AOD error of S1 in Scenario 2
(b)
AOD error of S2 in Scenario 2Figure 7
(a)
AOD error of S1 in Scenario 3
(b)
AOD error of S2 in Scenario 3


## 6. Conclusions

We analyzed the AOD for Lagrangian navigation constellation based on accurate dynamical mode. A neural network based method is applied to estimate the solar gravitational perturbation and the state of the Lagrangian satellites. The update rate of neural network weights is obtained by convergence of Lyapunov function. Therefore, the stability of the method introduced in this paper is proven. We derived that the constellation must include more than four satellites when this method is used. The simulation results show that this method can achieve a stable AOD error about 10 m only using satellite-to-satellite range as observation. Since the CRTBP is an approximate model for the Earth-Moon system, the AOD based on perturbed ERTBP is much useful for the future Lagrangian navigation constellation. We only test a Lagrangian navigation constellation with four satellites around L1 libration in this paper. Future work will be focused on a more complex constellation with satellites around different librations.



---

*Source: 9734164-2017-06-21.xml*