- Email: [email protected]
The Characteristics and Related Problems of the Orbits Around the Earth-Moon Libration Points
CHINESE ASTRONOMY AND ASTROPHYSICS Chinese Astronomy and Astrophysics 43 (2019) 278–291
The Characteristics and Related Problems of the Orbits Around the Earth-Moon Libration Points† LIU Lin1,2 1 2
TANG Jing-shi1,2
HOU Xi-yun1,2
School of Astronomy and Space Science, Nanjing University, Nanjing 210046
Institute of Space Environment and Astrodynamics, Nanjing University, Nanjing 210093
Abstract Due to the specific dynamics, the probes located at the halo orbits or Lissajous orbits around the Earth-Moon collinear libration point L1 or L2 are always studied in the synodic system to understand their trajectories. In fact, they are also orbiting the Earth in a distant Keplerian ellipse. Because of their intrinsic orbital instability, in the orbit prediction the initial errors propagate more prominently than those of the normal orbiting satellites, this requires special attention in the orbit design, maneuver, and control. Despite of all this, they are similar to the normal orbiting satellites in orbit determination and hardly require other special attentions. In this paper, the quantitative results of error propagation under the unstable dynamics, together with the theoretical analysis are presented. The results of precise orbit determination and short-arc orbit predictions are also shown, and compared with the results from the Beijing Aerospace Control Center. Key words celestial mechanics—Earth—Moon—orbit prediction 1.
INTRODUCTION
The orbit propagation (or orbit prediction) of the probe locating around the Earth-Moon collinear equilibrium points is not a new problem. The spacecraft of this type is actually a distant Earth satellite with specific orbital features. Although in theory its orbital period around the Earth is the same as the Moon’s, the lunar gravitation disturbance suffered by it is not an ignorable perturbation, but an external force that is comparable to the Earth Received 2017–12–25; revised version 2018–02–06
A translation of Acta Astron. Sin. Vol. 59, No. 3, pp. 29.1–29.12, 2018 [email protected]
[email protected]
c Elsevier 0275-1062/01/$-see front matter 2019 Elsevier B. V. All rights reserved. 0275-1062/01/$-see front matter © 2019 B. V. AllScience rights reserved. doi:10.1016/j.chinastron.2019.04.008 PII:
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
279
gravitation. Therefore, such a problem cannot be simply treated as a simple perturbed two-body problem. If being dealt in the Earth-Moon system as to construct the special trajectory (such as a halo trajectory) around the libration point, this problem cannot be handled as a simple disturbed two-body problem that is usually good enough for the study of a satellite’s motion. The difficulties in this problem arise from the following two aspects: (1) Even if we adopt the circular restricted three-body problem as the basic model[1−3] , the solutions are still complicated. If we use the elliptical restricted three-body problem by taking into account of the lunar orbit’s eccentricity, the solutions will be much more complicated, and more importantly, it does not help. Because the effect of solar gravitation perturbation is almost equivalent to the effect of lunar eccentricity, no matter which of the above-mentioned models is adopted, it is impossible to construct a perturbation solution as in the usual perturbed two-body problem. Thus, for the study of the motion of such probe, the model of perturbed two-body problem is useless in practice. (2) Even if we treat this problem as a perturbed two-body problem only formally, the difficulty still arises due to the comparable distances of the probe from the Moon and from the Earth. This is always a substantial difficulty in the perturbation theory of the third body. In view of the above analysis, for the purpose of ground TT & C, it is advisable to predict the orbits of the probe around the libration point of the Earth-Moon system in the J2000 geocentric celestial coordinate system using a numerical method. The mathematical model formally is a perturbed two-body problem, and the equation of motion is as follows: N μ r¨ = − 3 r + Fj (mj ) + Fε (ε) , r j=1
(1)
where r¨ is the acceleration of the probe, r is the distance from the probe to the Earth, μ = G (E + m), and G is the gravitation constant. The mass of probe m = 0 is ignorable, E is the Earth mass, and mj (j = 1, 2, · · ·) are the masses of perturbing objects among which the Moon (m1 ) and the Sun (m2 ) are the major ones. Fj stand for the perturbing accelerations from the third bodies, and N is the number of all the third bodies. In addition, Fε (ε) represents the perturbing acceleration caused by other forces (including the nonspherical gravitations of different celestial bodies, solar radiation pressure, etc.), and ε is a small perturbation parameter. 2.
FORCE ANALYSIS
The orbit prediction of spacecraft is a quantitative problem. The largely different dimensions and sizes of the involved physical quantities bring great inconvenience to the description and analysis of the problem, thus we propose the dimensionless unit (normalized unit) system. For the Earth satellite’s motion, including the distant objects such as the lunar probe around the L1 or L2 point of the Earth-Moon system, usually the normalized unit system is adopted
280
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
as follows. The units of length [L], mass [M ] and time [T ] are defined as: [L] = ae
[M ] = E
[T ] = [a3e /(GE)]1/2 ,
(2)
where GE is the geocentric gravitation constant, ae is the equatorial radius of the Earth reference ellipsoid. The unit of time [T ] is induced so that in this normalized system the gravitation constants G = 1 and μ = GE = 1. Taken the WGS-84 in the Earth-fixed coordinate system as the Earth gravitation model, the dynamical oblateness of the Earth is J2 = 1.082636022 × 10−3 . Adopting the above normalized units, the equation of motion of Eq. (1) becomes N 1 r¨ = − 3 r + Fj mj + Fε (ε) , r j=1
(3)
where mj are the dimensionless masses of the third bodies, such as the Moon and the Sun, i.e. m1 = GM /(GE), m2 = GS/(GE), · · · , (4) where GM and GS are the selenocentric and heliocentric gravitation constants, respectively. And, the acceleration due to the perturbation of the third body in Eq. (3) Fj mj reads j rj Δ Fj mj = −mj (5) + 3 , Δ3j rj j is the position vector of the probe with respect to the j-th object, Δ j = r − where Δ rj (j = 1, 2, · · ·), Δj is the distance from the probe to the j-th object, and rj represents the position vectors of the third bodies (the Moon, the Sun, etc.) in the geocentric celestial coordinate system. 2.1
Order of Magnitude Estimation for the Gravitational Perturbations from
Major Celestial Objects The gravitation constant of the Earth is GE=398600.4418 km3 /s2 , with respect to the Earth, the relative gravitation constants of the Moon, the Sun, Mercury, Venus, Mars, Jupiter and Saturn (given in GE) are 0.0123000383, 332946.050895, 0.055273598, 0.814998108, 0.107446732, 317.8942053 and 95.1574041, respectively. The orders of magnitude of the perturbations are estimated as εj = mj
r rj
3 ,
(6)
where mj expresses the mass ratios of the different major celestial objects with respect to the Earth, their values are the relative gravitation constants given above. The r and rj are the geocentric distances of the probe and perturbing object. In estimations, we can adopt the
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
281
average value of rj . But for Mercury, taking into account of its large eccentricity (e 0.2), the corresponding r will be set for the two situations corresponding to its perihelion and aphelion. The distance r from the Earth to the probe around the collinear libration point L1 or L2 in the Earth-Moon system is respectively ⎧ ⎨ r = 0.849065782EL = 51.2a 1 e , (7) ⎩ r2 = 1.167832664EL = 70.4ae where EL is the mean Earth-Moon distance. The order of magnitude for the perturbations from the above-mentioned major celestial objects on the probe around the L1 and L2 points can be estimated using Eq. (6) as: ⎧ ⎪ 0.8 × 10−2 , 2.0 × 10−2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 3.5 × 10−3 , 0.9 × 10−2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (3.8 ∼ 1.7) × 10−9 , (1.0 ∼ 0.45) × 10−8 ⎨
r 3 ⎪ ε = m = 4.0 × 10−7 , . (8) 1.0 × 10−6 ⎪ r ⎪ ⎪ ⎪ 0.8 × 10−8 , 2.0 × 10−8 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 4.5 × 10−8 , 1.2 × 10−7 ⎪ ⎪ ⎪ ⎪ ⎩ 1.1 × 10−9 , 2.9 × 10−9 2.2 Estimation of Earth Non-spherical Gravitation Perturbation The order of magnitude for the perturbation of the principal term (oblateness J2 ) is estimated as ⎧
⎨ 1.2 × 10−6 , L1 3 ε = ε (J2 ) = J2 = . (9) ⎩ 0.6 × 10−6 , L2 r2 2.3 Perturbation of Solar Radiation Pressure For the probe located around the L1 and L2 points in the Earth-Moon system with an ordinary scale (including mass and the equivalent cross section of solar irradiation), the order of magnitude of the solar radiation pressure can be estimated as: ⎧ ⎨ 1.2 × 10−5 , L κS 1 ε = ε (ρ ) = , (10) ρ r 2 = ⎩ 2.3 × 10−5 , L2 m where κ = 1.44, the surface area-to-mass ratio S/m = 109 , and ρ = 0.3169 × 10−17 is the solar radiation pressure at 1 au. According to the above order of magnitude estimation on the external force perturbations, if we consider only the perturbations larger than 10−6 for the motion of probe around the L1 and L2 points in the Earth-Moon system, the perturbations that should be taken
282
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
into account include the gravitational perturbations from the Moon, the Sun and Venus, the Earth’s non-spherical gravitation perturbation of oblateness J2 , and the solar radiation pressure. Among them, the most important perturbations are the mass-point gravitations from the Moon and the Sun. 3.
OVERVIEW OF PROBE ORBITS AROUND THE EARTH-MOON L1 AND L2 POINTS
Here we present two numerical simulations of probe orbits in the circular restricted threebody problem with the Earth and the Moon as primaries. The epoch is set as 201609-30UTC0:00:0.0, corresponding to TDT (terrestrial dynamical time) 57661.0007891667 (MJD). The probe is located on one orbit around each of the L1 and L2 points in the Earth-Moon system. The initial position (x, y, z), velocity (x, ˙ y, ˙ z), ˙ and the corresponding orbital elements in the J2000.0 geocentric celestial coordinate system for the two orbits can be calculated after a simple coordinate transformation, as listed in Table 1 and Table 2, respectively. Table 1 The initial positions and velocities of the orbits around L1 and L2 x/km
point
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
L1
–337774.810825 –337774.810825 16503.725924 –0.08734133 –0.793809190 –0.262191839
L2
–464586.522898
27663.672934
22699.764448 –0.120132126 –1.091831124 –0.360627231
Table 2 The orbital elements of the orbits around L1 and L2 point
a/km
e
i/◦
Ω/◦
ω/◦
M/◦
L1
242063.297
0.40063038
18.507748
4.970405
353.948868
353.948868
L2
1064951.700
0.56309178
18.507748
4.970405
166.026309
1.189763
As shown in Fig. 1 and Fig. 2, the two orbits are initially elliptical orbits with a large eccentricity around the Earth, and the probe locates at the apogee and perigee, respectively. The unit ae adopted in these figures is the equatorial radius of the Earth reference ellipsoid. For the system of Earth-Moon plus probe, these orbits are all the initial instantaneous orbits, and under the influence of lunar gravitation, the probe and the Moon perform the same circular orbital motion “synchronously”.
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
283
Fig. 1 The figure (in the equatorial plane) of initial orbit around L1 point in the J2000 reference system.
Fig. 2 The same as Fig. 1 but for the orbit around L2 point.
4. 4.1
ANALYSIS ON THE PREDICTION OF LIBRATION POINT ORBITS Mathematical Model for the Libration Point Orbit Analysis
As mentioned in Section 1, this problem should be considered in the J2000 geocentric celestial coordinate system and the orbit prediction should be performed via the numerical method. To show quantitatively the main features of the error propagation in the orbit prediction, we select the mass-point gravitation system consisting of Earth-Moon-Sun plus the probe, with the corresponding equation of motion being ¨ = F0 + F1 (m 1 ) + F2 (m 2 ) , r
(11)
284
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
where F0 = −r/r3 , is the acceleration due to the Earth’s gravitation, and m1 and m2 are the dimensionless masses of the Moon and the Sun, respectively (see Eq. (4)). The perturbation accelerations due to the Moon and the Sun F1 (m1 ) , F2 (m2 ) are expressed as 1 2 Δ Δ r r 1 2 (12) F1 (m1 ) = −m1 + 3 , F2 (m2 ) = −m2 + 3 . Δ1 3 Δ2 3 r1 r2 4.2 Types of Libration Point Orbits We consider three types of libration point orbits as follows: (1) Orbits initially located right at the L1 or L2 point in the Earth-Moon system. For short, they will be called L1 orbit or L2 orbit. (2) Orbits initially located on the halo orbit around the L1 or L2 point in the EarthMoon system. For short, hereafter they will be called L1 halo orbit or L2 halo orbit. (3) Orbits initially located on the Lissajous orbit around the L1 or L2 point in the Earth-Moon system. They will be called L1 Lissajous orbit or L2 Lissajous orbit. A preliminary design in the adopted mass-point gravitation system gives 6 orbits, whose initial conditions in the J2000 geocentric celestial coordinate system are described as follows: the initial epoch is 2016-09-30UTC0:00:0.0 (corresponding to TDT of MJD57661.0007891667), and the initial positions, velocities and orbital elements are listed in Table 3 and Table 4. In the tables, the orbit types 1, 2, · · · , 6 represent in sequence the L1 orbit, L1 halo orbit, L1 Lissajous orbit, L2 orbit, L2 halo orbit, and L2 Lissajous orbit, respectively. The same notation is adopted in Tables 4–8. In the following we will perform the orbit predictions of 7 d and 27 d for each of these 6 orbits, in order to obtain the quantitative profiles of the error propagations, then on this basis to make further qualitative analyse. Table 3 The initial positions and velocities of 6 orbits type
x/km
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
1
–337774.812386
20112.705545
16503.726001
–0.060238672
–0.795423014
–0.263516080
2
–345532.034293
23412.210001
18744.897406
18744.897406
–0.681345631
–0.223235504
3
–339800.358749
18612.308488
17838.793255
–0.084919606
–0.758511545
–0.244765666
4
–464586.522898
27663.672934
22699.764448
–0.120132126
–1.091831124
–0.360627231
5
–474559.546159
30334.673645
24471.171450
–0.108738242
–0.966961696
–0.318679777
6
–468100.537842
27030.333913
27058.570893
–0.116389398
–1.050980229
–0.346236851
4.3 Orbit Prediction of 7 d for 6 Libration Point Orbits For the simplicity (without loss of generality), we focus on the error propagation of semimajor axis (currently, the orbit determination has a precision of semi-major axis in the order of magnitude of 10 m). The results of orbit prediction of 7 d are listed in Table 5 and Table 6.
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
285
Table 4 The initial orbital elements of 6 orbits type
a/km
e
i/◦
Ω/◦
ω/◦
M/◦
1
2419 54.728
0.40015996
18.507748
4.970405
351.172034
180.000000
2
224124.112
0.54793918
18.406015
5.484455
351.389485
176.433825
3
234898.998
0.45201084
18.179965
6.049772
352.984029
173.757384
4
1064951.700
0.56309178
18.507748
4.970405
166.026309
1.189763
5
636901.958
0.25401430
18.487234
5.196394
162.260938
4.846848
6
867630.708
0.45956498
18.568585
6.587056
163.755904
1.916194
Table 5 The states of the libration point orbits propagated for 7 days type 1
2
3
4
5
6
4.4
model
a/km
e
i/◦
Ω/◦
λ/◦
standard
220999.515
0.415332
18.520032
5.002991
313.907599
Δa0 = 10 m
220999.532
0.415332
18.520032
5.002991
313.907501
standard
266442.683
0.267042
18.718604
4.923182
256.745900
Δa0 = 10 m
266443.184
0.267041
18.718603
4.923173
256.745426
standard
249439.450
0.364272
18.916482
5.301399
257.378240
Δa0 = 10 m
249439.933
0.364271
18.916479
5.301379
257.377737
standard
980939.283
0.536729
18.534082
5.083138
265.994812
Δa0 = 10 m
980939.634
0.536729
18.534082
5.083138
265.994765
standard
2623047.427
0.826111
18.614129
5. 248427
257.032873
Δa0 = 10 m
2623040.349
0.826110
18.614129
5. 248428
257.032827
standard
1223545.852
0.619375
19.078832
6.073193
257.502177
Δa0 = 10 m
1223545.763
0.619375
19.078833
6.073195
257.502134
Orbit Prediction of 27 d for 6 Libration Point Orbits
The location of the probe is only approximately at the libration point. In practice, continuous maneuver has to be made to maintain this kind of orbits. Therefore, the orbit prediction as long as 27 d is just for better understanding the dynamical features of this kind of orbits and the corresponding law of error propagation. With all initial errors being set on the semi-major axis, the results of orbit prediction of 27 d are listed in Table 7 and Table 8, in which the 2nd and 3rd type orbits (i.e. the L1 halo orbit and the L1 Lissajous orbit) have been propagated only for 22 days (see the reason in Section 4.5). As for the probe trajectory, the “long-term” or “short-term” in the position prediction indicates not only the long or short time interval, but also the length of the orbit arc. Therefore, in order to compare quantitatively the error propagations in the above-mentioned types of orbits, we need to know their orbital periods. The initial periods (TS ) of the six orbits
286
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
are TS =13.708768 d, 12.221658 d, 13.113512 d for the orbits around L1 and TS =126.58784 d, 58.847196 d, 93.089267 d for the orbits around L2 . Now we know the corresponding arc lengths of the six orbits for the 7 d and 27 d predictions. Table 6 The states of spatial positions of the libration point orbits propagated for 7 days type 1
2
3
4
5
6
model
x/km
y/km
z/km
r/km
standard
17349.352
–235429.496
–79071.509
248958.533
Δa0 = 10 m
17349.130
–235429.732
–79071.581
248958.763
standard
–52603.690
–316406.575
–105287.004
337587.947
Δa0 = 10 m
–52604.014
–316407.100
–105287.168
337588.541
standard
–61169.892
–317191.268
–106298.963
340075.765
Δa0 = 10 m
–61170.234
–317191.785
–106299.122
340076.358
standard
–75890.039
–430172.798
–141396.957
459130.737
Δa0 = 10 m
–75890.082
–430172.853
–141396.974
459130.801
standard
–96023.763
–424222.928
–139325.760
456724.559
Δa0 = 10 m
–96023.848
–424222.971
–139325.772
456724.621
standard
–89282.870
–433942.578
–145977.148
466462.131
Δa0 = 10 m
–89282.927
–433942.628
–145977.163
466462.193
Table 7 The states of the libration point orbits propagated for 27 days type 1
2
3
4
5
6
model
a/km
e
i/◦
Ω/◦
λ/◦
standard
212685.965
0.442498
18.517044
5.009819
211.007387
Δa0 = 10 m
212685.969
0.442498
18.517044
5.009819
211.007270
standard
222801.427
0.515663
18.570547
6.333758
88.715061
Δa0 = 10 m
247840.427
0.595452
18.616593
6.037248
56.512170
standard
235375.045
0.425546
18.494325
7.675744
85.777531
Δa0 = 10 m
378708.222
0.532445
18.639545
6.508178
37.783963
standard
1187150.224
0.614869
18.614143
5.245764
164.740828
Δa0 = 10 m
1187149.042
0.614869
18.614143
5.245764
164.740913
standard
690495.272
0.333894
18.560056
5. 215658
162.442776
Δa0 = 10 m
704825.483
0.346515
18.560165
5. 213139
161.458414
standard
860189.818
0.462166
17.811203
6.065852
165.371342
Δa0 = 10 m
875504.051
0.469931
17.813259
6.041064
164.128126
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
287
Knowing the dynamical features of the six orbits and the arc lengths of orbit predictions, we can summarize some dynamical laws reflected by the position errors listed in Table 5 as follows: (1) Orbits propagated for 7 days are all short arcs, with position errors smaller than 1 km. (2) For the orbits around the L1 and L2 points, the accumulation of position error in 7 d or 27 d is not serious, and the error propagation is just the reflection of the Kepler motion. Since the initial error (Δa0 = 10 m) is small, it can be regarded as a small perturbation. It will neither provoke the intrinsic instability nor change the effect of error accumulation in a short arc. On the contrary, the periodic effects of the error are more significant than the effect of long-term accumulation (see for example the error of Type 1 orbit in Table 5). (3) For the halo and Lissajous orbits around L1 and L2 , they are in fact “far away from” the unstable libration points L1 and L2 , while the strict designs of the halo orbit and Lissajous orbit are unable to realize, the locating of the probe at the libration points in space is only an approximation. For the given initial error of Δa0 = 10 m, it is rather a large disturbance to the unstable libration points than a small perturbation to the halo or Lissajous orbit. Even in a relatively not long predicting arc, the instability of the libration points may demonstrate themselves completely. See the orbits of Type 2, 3, 5 and 6 in Table 5, especially Types 2 and 3, relatively the arc length of 27 d seems to be even longer, their orbital eccentricities will go up quickly to e ≈ 1.0 after 22 days. Table 8 The states of spatial positions of the libration point orbits propagated for 27 days type 1
2
3
4
5
6
model
x/km
y/km
z/km
r/km
standard
–231950.200
–191657.303
–57161.753
306269.298
Δa0 = 10 m
–231950.225
–191657.248
–57161.734
306269.280
standard
–115564.838
294821.738
102728.456
332908.733
Δa0 = 10 m
–110340.192
314798.221
109363.932
351046.076
standard
–106522.282
294300.982
102316.285
329285.115
Δa0 = 10 m
–106322.537
324499.131
112814.962
359626.728
standard
–452814.170
55918.600
32699.152
457424.089
Δa0 = 10 m
–452813.996
55918.681
32699.174
457423.928
standard
–464361.937
21794.997
21461.005
465368.247
Δa0 = 10 m
–465954.469
22340.406
21685.299
466993.482
standard
–460599.098
39263.652
28181.440
463127.798
Δa0 = 10 m
–462359.879
40264.443
28501.149
464984.084
288
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
4.5 Features of Error Propagation of Libration Point Orbits The basic features of the position error propagation of the Lagrange point orbit can be found from the above results of orbit prediction of 6 orbits. We summarize the position errors arising from the initial error of Δa0 = 10 m in Table 9. Table 9 The quantitative state of position errors of the Lagrange point orbit prediction (unit: km) type
propagating arc/d 1
2
3
4
5
6
7
0.322
0.638
0.640
0.072
0.096
0.077
22
–
186181.8
–47023.6
–
–
–
27
0.063
–
–
0.193
1698.216
2050.401
In summary, the simple calculation and analysis on the error propagation of the libration point orbits show that it is difficult to achieve a maneuverless long-arc trajectory around the L1 and L2 points in the Earth-Moon system. This difficulty is intrinsic and cannot be avoided by orbit design. But for the short-arc orbit prediction, it is still easy to achieve a relatively high accuracy. Specific examples will be shown later in this section, and the verification by simulations will also be presented. 4.6 The Status of L1 Halo Orbit In order to further reveal the dynamical mechanism of error propagation in the libration point orbits, we present here a large retrograde elliptical orbit (corresponding to a remote Earth satellite) that is similar to the above-mentioned L1 halo orbit. The initial position, velocity and the corresponding orbital elements at the initial epoch 2016-09-30UTC0:00:0.0 are listed in Table 10 and Table 11, respectively. The initial period is 12.208831 d. To make the 27 d prediction on this orbit, the calculated results are listed in Table 12 and Table 13. Table 10 The position and velocity of the L1 halo orbit x/km
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
–345307.466263
21904.224767
19425.079239
–0.077571217
0.681599717
–0.223818048
Table 11 The orbital elements of the L1 halo orbit a/km
e
i/◦
Ω/◦
ω/◦
M/◦
223967.268
0.56141578
160.912742
160.912742
356.972610
356.972610
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
289
Table 12 The state of the L1 halo orbit propagated for 27 days model
a/km
e
i/◦
Ω/◦
λ/◦
standard
213214.203
0.653838
162.168040
350.586870
267.258205
Δa0 = 10 m
213214.209
0.653838
162.168042
347.586870
267.258158
Table 13 The position state of the L1 halo orbit propagated for 27 days model
x/km
y/km
z/km
r/km
standard
–229245.266
186181.824
–47023.643
299045.626
Δa0 = 10 m
–229245.393
186181.802
–47023.625
299045.707
Such a large elliptical retrograde orbit, in difference from the L1 halo orbit, is essentially an ordinary Kepler orbit. The inherent instability is not manifested in an arc of 27 d (only 2 periods). We present this trivial (in a sense) result to show from an other angle the inherent instability of the rapid error propagation of the libration point orbits in the Earth-Moon system. The motion and dynamics of the retrograde orbit however are beyond the scope of this paper. Readers can refer to Ref.[4] for details. 4.7
Determination of Libration Point Orbits and Verification of Short-arc Pre-
diction Accuracy Since the design of orbits strictly on the Lagrangian points is impossible in reality, the positioning of the probe on these points is only an approximation. The orbit maintenance must be continuously performed during operation, therefore we should consider only the short-arc problem for both the ground TT&C and on-orbit maneuver. Compared to the LEO and HEO satellites, the problems and difficulties met in the short-arc orbit determination and prediction of libration point orbits are the same. Using the orbit determination software developed by the Institute of Space Environment and Astrodynamics of Nanjing University and the USB (Unified S-band) observation data in China, without any other extra information, we have determined the orbit of the probe of the Chang’e 3 mission. The results are verified by comparing them with the post hoc data released by the Beijing Aerospace Flight Control Center (BAFCC). Based on this orbit determination, a very simple numerical extrapolation is adopted to predict the orbit. In the calculation, only the mass-point gravitations from the Earth, Moon and Sun, and a simple model of radiation pressure are taken into account, and the orbit determination given by the relevant mission is directly adopted as the initial orbital elements. A rather high accuracy is achieved without any difficulties. Omitting the unnecessary details, we list the results in Table 14 and Table 15. In Table 14 and Table 15, the orbit determinations A and B represent the results from the BAFCC and our calculations, respectively. Although the results in these tables
290
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
are basically able to explain the problem, we would like to give some further explanations as follows for a clearer understanding on the accuracies of the orbit determination and prediction for this kind of probes. (1) For the radiation pressure, in the condition without the specific details of the probe, we empirically use an equivalent planar model with an area-mass ratio obtained from the related estimates in our independent orbit determinations. (2) Although the details of both orbit determinations (including the observational data) are not specified, the calculation results in Tables 14–15 have demonstrated the accuracies of the orbit determination and prediction for such kind of orbits in this paper. In this two tables, the results of the 3 d and 7 d orbit predictions are compared with the results of post hoc precise orbit determinations, and the differences between the results of both orbit determinations are basically within 500 m. Such a comparison proves more strongly the reliability of the results of both orbit determinations and the rationality of the force model adopted by this paper. Table 14 The comparisons of two precise orbit determinations with the result of 3 d orbit prediction for the L2 halo orbit method
x/km
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
orbit determination A
–451897.422
138113.307
36463.951
–0.237782
–0.978285
–0.204902
orbit determination B
–451897.388
138112.927
36463.914
–0.237782
–0.978287
–0.204900
short arc propagation A
–451897.458
138113.145
36463.888
–0.237783
–0.978286
–0.204902
short arc propagation B
–451897.346
138112.816
36463.866
–0.237782
–0.978288
–0.204900
Table 15 The comparisons of two precise orbit determinations with the result of 7 d orbit prediction for the L2 halo orbit method
x/km
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
orbit determination A
–371590.823
–223419.646
–55369.560
0.740132
–0.975794
–0.327870
orbit determination B
–371590.387
–223420.334
–55369.450
0.740133
–0.975794
–0.327871
short arc propagation A
–371590.759
–223420.292
–55369.879
0.740132
–0.975795
–0.327871
short arc propagation B
–371589.646
–223420.816
–55369.813
0.740138
–0.975793
–0.327873
The above calculations and the supplementary explanations indicate thoroughly that no special difficulties exist in the short-arc orbit determination and high-precision orbit prediction for this kind of libration point orbits, although the propagation of initial errors in these special orbits is much more significant than that of usual encircling orbits.
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
291
4.8 Transformation of Libration Point Orbit between Two Coordinate Systems The J2000 geocentric celestial coordinate system is required for orbit determination and prediction, but for the purpose of orbital maneuver the Earth-Moon rotating coordinate system is needed, bring up the necessity of coordinate conversion between these two systems. This conversion is easy to implement. For a given space mission, a specific rotating coordinate system that is most suitable for the mission may be defined according to actual needs, we will not make further discussions. References 1
Liu L., Tang J. S., Theory and Application of Satellite Orbits, Beijing: Publishing House of Electronics Industry, 2015
2
Gomez G., Llibre J., Martinez R., et al., Dynamics and Mission Design Near Libration Points. Sin-
3
Liu L., Hou X. Y., Theory and Application of Deep Space Exploring Orbits, Beijing: Publishing
4
Zhao C. Y., Liu L., AcASn, 1994, 35, 434
gapore, New Jersey, London, Hong Kong: World Scientific, 2001 House of Electronics Industry, 2015
The Characteristics and Related Problems of the Orbits Around the Earth-Moon Libration Points† LIU Lin1,2 1 2
TANG Jing-shi1,2
HOU Xi-yun1,2
School of Astronomy and Space Science, Nanjing University, Nanjing 210046
Institute of Space Environment and Astrodynamics, Nanjing University, Nanjing 210093
Abstract Due to the specific dynamics, the probes located at the halo orbits or Lissajous orbits around the Earth-Moon collinear libration point L1 or L2 are always studied in the synodic system to understand their trajectories. In fact, they are also orbiting the Earth in a distant Keplerian ellipse. Because of their intrinsic orbital instability, in the orbit prediction the initial errors propagate more prominently than those of the normal orbiting satellites, this requires special attention in the orbit design, maneuver, and control. Despite of all this, they are similar to the normal orbiting satellites in orbit determination and hardly require other special attentions. In this paper, the quantitative results of error propagation under the unstable dynamics, together with the theoretical analysis are presented. The results of precise orbit determination and short-arc orbit predictions are also shown, and compared with the results from the Beijing Aerospace Control Center. Key words celestial mechanics—Earth—Moon—orbit prediction 1.
INTRODUCTION
The orbit propagation (or orbit prediction) of the probe locating around the Earth-Moon collinear equilibrium points is not a new problem. The spacecraft of this type is actually a distant Earth satellite with specific orbital features. Although in theory its orbital period around the Earth is the same as the Moon’s, the lunar gravitation disturbance suffered by it is not an ignorable perturbation, but an external force that is comparable to the Earth Received 2017–12–25; revised version 2018–02–06
A translation of Acta Astron. Sin. Vol. 59, No. 3, pp. 29.1–29.12, 2018 [email protected]
[email protected]
c Elsevier 0275-1062/01/$-see front matter 2019 Elsevier B. V. All rights reserved. 0275-1062/01/$-see front matter © 2019 B. V. AllScience rights reserved. doi:10.1016/j.chinastron.2019.04.008 PII:
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
279
gravitation. Therefore, such a problem cannot be simply treated as a simple perturbed two-body problem. If being dealt in the Earth-Moon system as to construct the special trajectory (such as a halo trajectory) around the libration point, this problem cannot be handled as a simple disturbed two-body problem that is usually good enough for the study of a satellite’s motion. The difficulties in this problem arise from the following two aspects: (1) Even if we adopt the circular restricted three-body problem as the basic model[1−3] , the solutions are still complicated. If we use the elliptical restricted three-body problem by taking into account of the lunar orbit’s eccentricity, the solutions will be much more complicated, and more importantly, it does not help. Because the effect of solar gravitation perturbation is almost equivalent to the effect of lunar eccentricity, no matter which of the above-mentioned models is adopted, it is impossible to construct a perturbation solution as in the usual perturbed two-body problem. Thus, for the study of the motion of such probe, the model of perturbed two-body problem is useless in practice. (2) Even if we treat this problem as a perturbed two-body problem only formally, the difficulty still arises due to the comparable distances of the probe from the Moon and from the Earth. This is always a substantial difficulty in the perturbation theory of the third body. In view of the above analysis, for the purpose of ground TT & C, it is advisable to predict the orbits of the probe around the libration point of the Earth-Moon system in the J2000 geocentric celestial coordinate system using a numerical method. The mathematical model formally is a perturbed two-body problem, and the equation of motion is as follows: N μ r¨ = − 3 r + Fj (mj ) + Fε (ε) , r j=1
(1)
where r¨ is the acceleration of the probe, r is the distance from the probe to the Earth, μ = G (E + m), and G is the gravitation constant. The mass of probe m = 0 is ignorable, E is the Earth mass, and mj (j = 1, 2, · · ·) are the masses of perturbing objects among which the Moon (m1 ) and the Sun (m2 ) are the major ones. Fj stand for the perturbing accelerations from the third bodies, and N is the number of all the third bodies. In addition, Fε (ε) represents the perturbing acceleration caused by other forces (including the nonspherical gravitations of different celestial bodies, solar radiation pressure, etc.), and ε is a small perturbation parameter. 2.
FORCE ANALYSIS
The orbit prediction of spacecraft is a quantitative problem. The largely different dimensions and sizes of the involved physical quantities bring great inconvenience to the description and analysis of the problem, thus we propose the dimensionless unit (normalized unit) system. For the Earth satellite’s motion, including the distant objects such as the lunar probe around the L1 or L2 point of the Earth-Moon system, usually the normalized unit system is adopted
280
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
as follows. The units of length [L], mass [M ] and time [T ] are defined as: [L] = ae
[M ] = E
[T ] = [a3e /(GE)]1/2 ,
(2)
where GE is the geocentric gravitation constant, ae is the equatorial radius of the Earth reference ellipsoid. The unit of time [T ] is induced so that in this normalized system the gravitation constants G = 1 and μ = GE = 1. Taken the WGS-84 in the Earth-fixed coordinate system as the Earth gravitation model, the dynamical oblateness of the Earth is J2 = 1.082636022 × 10−3 . Adopting the above normalized units, the equation of motion of Eq. (1) becomes N 1 r¨ = − 3 r + Fj mj + Fε (ε) , r j=1
(3)
where mj are the dimensionless masses of the third bodies, such as the Moon and the Sun, i.e. m1 = GM /(GE), m2 = GS/(GE), · · · , (4) where GM and GS are the selenocentric and heliocentric gravitation constants, respectively. And, the acceleration due to the perturbation of the third body in Eq. (3) Fj mj reads j rj Δ Fj mj = −mj (5) + 3 , Δ3j rj j is the position vector of the probe with respect to the j-th object, Δ j = r − where Δ rj (j = 1, 2, · · ·), Δj is the distance from the probe to the j-th object, and rj represents the position vectors of the third bodies (the Moon, the Sun, etc.) in the geocentric celestial coordinate system. 2.1
Order of Magnitude Estimation for the Gravitational Perturbations from
Major Celestial Objects The gravitation constant of the Earth is GE=398600.4418 km3 /s2 , with respect to the Earth, the relative gravitation constants of the Moon, the Sun, Mercury, Venus, Mars, Jupiter and Saturn (given in GE) are 0.0123000383, 332946.050895, 0.055273598, 0.814998108, 0.107446732, 317.8942053 and 95.1574041, respectively. The orders of magnitude of the perturbations are estimated as εj = mj
r rj
3 ,
(6)
where mj expresses the mass ratios of the different major celestial objects with respect to the Earth, their values are the relative gravitation constants given above. The r and rj are the geocentric distances of the probe and perturbing object. In estimations, we can adopt the
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
281
average value of rj . But for Mercury, taking into account of its large eccentricity (e 0.2), the corresponding r will be set for the two situations corresponding to its perihelion and aphelion. The distance r from the Earth to the probe around the collinear libration point L1 or L2 in the Earth-Moon system is respectively ⎧ ⎨ r = 0.849065782EL = 51.2a 1 e , (7) ⎩ r2 = 1.167832664EL = 70.4ae where EL is the mean Earth-Moon distance. The order of magnitude for the perturbations from the above-mentioned major celestial objects on the probe around the L1 and L2 points can be estimated using Eq. (6) as: ⎧ ⎪ 0.8 × 10−2 , 2.0 × 10−2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 3.5 × 10−3 , 0.9 × 10−2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (3.8 ∼ 1.7) × 10−9 , (1.0 ∼ 0.45) × 10−8 ⎨
r 3 ⎪ ε = m = 4.0 × 10−7 , . (8) 1.0 × 10−6 ⎪ r ⎪ ⎪ ⎪ 0.8 × 10−8 , 2.0 × 10−8 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 4.5 × 10−8 , 1.2 × 10−7 ⎪ ⎪ ⎪ ⎪ ⎩ 1.1 × 10−9 , 2.9 × 10−9 2.2 Estimation of Earth Non-spherical Gravitation Perturbation The order of magnitude for the perturbation of the principal term (oblateness J2 ) is estimated as ⎧
⎨ 1.2 × 10−6 , L1 3 ε = ε (J2 ) = J2 = . (9) ⎩ 0.6 × 10−6 , L2 r2 2.3 Perturbation of Solar Radiation Pressure For the probe located around the L1 and L2 points in the Earth-Moon system with an ordinary scale (including mass and the equivalent cross section of solar irradiation), the order of magnitude of the solar radiation pressure can be estimated as: ⎧ ⎨ 1.2 × 10−5 , L κS 1 ε = ε (ρ ) = , (10) ρ r 2 = ⎩ 2.3 × 10−5 , L2 m where κ = 1.44, the surface area-to-mass ratio S/m = 109 , and ρ = 0.3169 × 10−17 is the solar radiation pressure at 1 au. According to the above order of magnitude estimation on the external force perturbations, if we consider only the perturbations larger than 10−6 for the motion of probe around the L1 and L2 points in the Earth-Moon system, the perturbations that should be taken
282
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
into account include the gravitational perturbations from the Moon, the Sun and Venus, the Earth’s non-spherical gravitation perturbation of oblateness J2 , and the solar radiation pressure. Among them, the most important perturbations are the mass-point gravitations from the Moon and the Sun. 3.
OVERVIEW OF PROBE ORBITS AROUND THE EARTH-MOON L1 AND L2 POINTS
Here we present two numerical simulations of probe orbits in the circular restricted threebody problem with the Earth and the Moon as primaries. The epoch is set as 201609-30UTC0:00:0.0, corresponding to TDT (terrestrial dynamical time) 57661.0007891667 (MJD). The probe is located on one orbit around each of the L1 and L2 points in the Earth-Moon system. The initial position (x, y, z), velocity (x, ˙ y, ˙ z), ˙ and the corresponding orbital elements in the J2000.0 geocentric celestial coordinate system for the two orbits can be calculated after a simple coordinate transformation, as listed in Table 1 and Table 2, respectively. Table 1 The initial positions and velocities of the orbits around L1 and L2 x/km
point
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
L1
–337774.810825 –337774.810825 16503.725924 –0.08734133 –0.793809190 –0.262191839
L2
–464586.522898
27663.672934
22699.764448 –0.120132126 –1.091831124 –0.360627231
Table 2 The orbital elements of the orbits around L1 and L2 point
a/km
e
i/◦
Ω/◦
ω/◦
M/◦
L1
242063.297
0.40063038
18.507748
4.970405
353.948868
353.948868
L2
1064951.700
0.56309178
18.507748
4.970405
166.026309
1.189763
As shown in Fig. 1 and Fig. 2, the two orbits are initially elliptical orbits with a large eccentricity around the Earth, and the probe locates at the apogee and perigee, respectively. The unit ae adopted in these figures is the equatorial radius of the Earth reference ellipsoid. For the system of Earth-Moon plus probe, these orbits are all the initial instantaneous orbits, and under the influence of lunar gravitation, the probe and the Moon perform the same circular orbital motion “synchronously”.
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
283
Fig. 1 The figure (in the equatorial plane) of initial orbit around L1 point in the J2000 reference system.
Fig. 2 The same as Fig. 1 but for the orbit around L2 point.
4. 4.1
ANALYSIS ON THE PREDICTION OF LIBRATION POINT ORBITS Mathematical Model for the Libration Point Orbit Analysis
As mentioned in Section 1, this problem should be considered in the J2000 geocentric celestial coordinate system and the orbit prediction should be performed via the numerical method. To show quantitatively the main features of the error propagation in the orbit prediction, we select the mass-point gravitation system consisting of Earth-Moon-Sun plus the probe, with the corresponding equation of motion being ¨ = F0 + F1 (m 1 ) + F2 (m 2 ) , r
(11)
284
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
where F0 = −r/r3 , is the acceleration due to the Earth’s gravitation, and m1 and m2 are the dimensionless masses of the Moon and the Sun, respectively (see Eq. (4)). The perturbation accelerations due to the Moon and the Sun F1 (m1 ) , F2 (m2 ) are expressed as 1 2 Δ Δ r r 1 2 (12) F1 (m1 ) = −m1 + 3 , F2 (m2 ) = −m2 + 3 . Δ1 3 Δ2 3 r1 r2 4.2 Types of Libration Point Orbits We consider three types of libration point orbits as follows: (1) Orbits initially located right at the L1 or L2 point in the Earth-Moon system. For short, they will be called L1 orbit or L2 orbit. (2) Orbits initially located on the halo orbit around the L1 or L2 point in the EarthMoon system. For short, hereafter they will be called L1 halo orbit or L2 halo orbit. (3) Orbits initially located on the Lissajous orbit around the L1 or L2 point in the Earth-Moon system. They will be called L1 Lissajous orbit or L2 Lissajous orbit. A preliminary design in the adopted mass-point gravitation system gives 6 orbits, whose initial conditions in the J2000 geocentric celestial coordinate system are described as follows: the initial epoch is 2016-09-30UTC0:00:0.0 (corresponding to TDT of MJD57661.0007891667), and the initial positions, velocities and orbital elements are listed in Table 3 and Table 4. In the tables, the orbit types 1, 2, · · · , 6 represent in sequence the L1 orbit, L1 halo orbit, L1 Lissajous orbit, L2 orbit, L2 halo orbit, and L2 Lissajous orbit, respectively. The same notation is adopted in Tables 4–8. In the following we will perform the orbit predictions of 7 d and 27 d for each of these 6 orbits, in order to obtain the quantitative profiles of the error propagations, then on this basis to make further qualitative analyse. Table 3 The initial positions and velocities of 6 orbits type
x/km
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
1
–337774.812386
20112.705545
16503.726001
–0.060238672
–0.795423014
–0.263516080
2
–345532.034293
23412.210001
18744.897406
18744.897406
–0.681345631
–0.223235504
3
–339800.358749
18612.308488
17838.793255
–0.084919606
–0.758511545
–0.244765666
4
–464586.522898
27663.672934
22699.764448
–0.120132126
–1.091831124
–0.360627231
5
–474559.546159
30334.673645
24471.171450
–0.108738242
–0.966961696
–0.318679777
6
–468100.537842
27030.333913
27058.570893
–0.116389398
–1.050980229
–0.346236851
4.3 Orbit Prediction of 7 d for 6 Libration Point Orbits For the simplicity (without loss of generality), we focus on the error propagation of semimajor axis (currently, the orbit determination has a precision of semi-major axis in the order of magnitude of 10 m). The results of orbit prediction of 7 d are listed in Table 5 and Table 6.
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
285
Table 4 The initial orbital elements of 6 orbits type
a/km
e
i/◦
Ω/◦
ω/◦
M/◦
1
2419 54.728
0.40015996
18.507748
4.970405
351.172034
180.000000
2
224124.112
0.54793918
18.406015
5.484455
351.389485
176.433825
3
234898.998
0.45201084
18.179965
6.049772
352.984029
173.757384
4
1064951.700
0.56309178
18.507748
4.970405
166.026309
1.189763
5
636901.958
0.25401430
18.487234
5.196394
162.260938
4.846848
6
867630.708
0.45956498
18.568585
6.587056
163.755904
1.916194
Table 5 The states of the libration point orbits propagated for 7 days type 1
2
3
4
5
6
4.4
model
a/km
e
i/◦
Ω/◦
λ/◦
standard
220999.515
0.415332
18.520032
5.002991
313.907599
Δa0 = 10 m
220999.532
0.415332
18.520032
5.002991
313.907501
standard
266442.683
0.267042
18.718604
4.923182
256.745900
Δa0 = 10 m
266443.184
0.267041
18.718603
4.923173
256.745426
standard
249439.450
0.364272
18.916482
5.301399
257.378240
Δa0 = 10 m
249439.933
0.364271
18.916479
5.301379
257.377737
standard
980939.283
0.536729
18.534082
5.083138
265.994812
Δa0 = 10 m
980939.634
0.536729
18.534082
5.083138
265.994765
standard
2623047.427
0.826111
18.614129
5. 248427
257.032873
Δa0 = 10 m
2623040.349
0.826110
18.614129
5. 248428
257.032827
standard
1223545.852
0.619375
19.078832
6.073193
257.502177
Δa0 = 10 m
1223545.763
0.619375
19.078833
6.073195
257.502134
Orbit Prediction of 27 d for 6 Libration Point Orbits
The location of the probe is only approximately at the libration point. In practice, continuous maneuver has to be made to maintain this kind of orbits. Therefore, the orbit prediction as long as 27 d is just for better understanding the dynamical features of this kind of orbits and the corresponding law of error propagation. With all initial errors being set on the semi-major axis, the results of orbit prediction of 27 d are listed in Table 7 and Table 8, in which the 2nd and 3rd type orbits (i.e. the L1 halo orbit and the L1 Lissajous orbit) have been propagated only for 22 days (see the reason in Section 4.5). As for the probe trajectory, the “long-term” or “short-term” in the position prediction indicates not only the long or short time interval, but also the length of the orbit arc. Therefore, in order to compare quantitatively the error propagations in the above-mentioned types of orbits, we need to know their orbital periods. The initial periods (TS ) of the six orbits
286
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
are TS =13.708768 d, 12.221658 d, 13.113512 d for the orbits around L1 and TS =126.58784 d, 58.847196 d, 93.089267 d for the orbits around L2 . Now we know the corresponding arc lengths of the six orbits for the 7 d and 27 d predictions. Table 6 The states of spatial positions of the libration point orbits propagated for 7 days type 1
2
3
4
5
6
model
x/km
y/km
z/km
r/km
standard
17349.352
–235429.496
–79071.509
248958.533
Δa0 = 10 m
17349.130
–235429.732
–79071.581
248958.763
standard
–52603.690
–316406.575
–105287.004
337587.947
Δa0 = 10 m
–52604.014
–316407.100
–105287.168
337588.541
standard
–61169.892
–317191.268
–106298.963
340075.765
Δa0 = 10 m
–61170.234
–317191.785
–106299.122
340076.358
standard
–75890.039
–430172.798
–141396.957
459130.737
Δa0 = 10 m
–75890.082
–430172.853
–141396.974
459130.801
standard
–96023.763
–424222.928
–139325.760
456724.559
Δa0 = 10 m
–96023.848
–424222.971
–139325.772
456724.621
standard
–89282.870
–433942.578
–145977.148
466462.131
Δa0 = 10 m
–89282.927
–433942.628
–145977.163
466462.193
Table 7 The states of the libration point orbits propagated for 27 days type 1
2
3
4
5
6
model
a/km
e
i/◦
Ω/◦
λ/◦
standard
212685.965
0.442498
18.517044
5.009819
211.007387
Δa0 = 10 m
212685.969
0.442498
18.517044
5.009819
211.007270
standard
222801.427
0.515663
18.570547
6.333758
88.715061
Δa0 = 10 m
247840.427
0.595452
18.616593
6.037248
56.512170
standard
235375.045
0.425546
18.494325
7.675744
85.777531
Δa0 = 10 m
378708.222
0.532445
18.639545
6.508178
37.783963
standard
1187150.224
0.614869
18.614143
5.245764
164.740828
Δa0 = 10 m
1187149.042
0.614869
18.614143
5.245764
164.740913
standard
690495.272
0.333894
18.560056
5. 215658
162.442776
Δa0 = 10 m
704825.483
0.346515
18.560165
5. 213139
161.458414
standard
860189.818
0.462166
17.811203
6.065852
165.371342
Δa0 = 10 m
875504.051
0.469931
17.813259
6.041064
164.128126
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
287
Knowing the dynamical features of the six orbits and the arc lengths of orbit predictions, we can summarize some dynamical laws reflected by the position errors listed in Table 5 as follows: (1) Orbits propagated for 7 days are all short arcs, with position errors smaller than 1 km. (2) For the orbits around the L1 and L2 points, the accumulation of position error in 7 d or 27 d is not serious, and the error propagation is just the reflection of the Kepler motion. Since the initial error (Δa0 = 10 m) is small, it can be regarded as a small perturbation. It will neither provoke the intrinsic instability nor change the effect of error accumulation in a short arc. On the contrary, the periodic effects of the error are more significant than the effect of long-term accumulation (see for example the error of Type 1 orbit in Table 5). (3) For the halo and Lissajous orbits around L1 and L2 , they are in fact “far away from” the unstable libration points L1 and L2 , while the strict designs of the halo orbit and Lissajous orbit are unable to realize, the locating of the probe at the libration points in space is only an approximation. For the given initial error of Δa0 = 10 m, it is rather a large disturbance to the unstable libration points than a small perturbation to the halo or Lissajous orbit. Even in a relatively not long predicting arc, the instability of the libration points may demonstrate themselves completely. See the orbits of Type 2, 3, 5 and 6 in Table 5, especially Types 2 and 3, relatively the arc length of 27 d seems to be even longer, their orbital eccentricities will go up quickly to e ≈ 1.0 after 22 days. Table 8 The states of spatial positions of the libration point orbits propagated for 27 days type 1
2
3
4
5
6
model
x/km
y/km
z/km
r/km
standard
–231950.200
–191657.303
–57161.753
306269.298
Δa0 = 10 m
–231950.225
–191657.248
–57161.734
306269.280
standard
–115564.838
294821.738
102728.456
332908.733
Δa0 = 10 m
–110340.192
314798.221
109363.932
351046.076
standard
–106522.282
294300.982
102316.285
329285.115
Δa0 = 10 m
–106322.537
324499.131
112814.962
359626.728
standard
–452814.170
55918.600
32699.152
457424.089
Δa0 = 10 m
–452813.996
55918.681
32699.174
457423.928
standard
–464361.937
21794.997
21461.005
465368.247
Δa0 = 10 m
–465954.469
22340.406
21685.299
466993.482
standard
–460599.098
39263.652
28181.440
463127.798
Δa0 = 10 m
–462359.879
40264.443
28501.149
464984.084
288
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
4.5 Features of Error Propagation of Libration Point Orbits The basic features of the position error propagation of the Lagrange point orbit can be found from the above results of orbit prediction of 6 orbits. We summarize the position errors arising from the initial error of Δa0 = 10 m in Table 9. Table 9 The quantitative state of position errors of the Lagrange point orbit prediction (unit: km) type
propagating arc/d 1
2
3
4
5
6
7
0.322
0.638
0.640
0.072
0.096
0.077
22
–
186181.8
–47023.6
–
–
–
27
0.063
–
–
0.193
1698.216
2050.401
In summary, the simple calculation and analysis on the error propagation of the libration point orbits show that it is difficult to achieve a maneuverless long-arc trajectory around the L1 and L2 points in the Earth-Moon system. This difficulty is intrinsic and cannot be avoided by orbit design. But for the short-arc orbit prediction, it is still easy to achieve a relatively high accuracy. Specific examples will be shown later in this section, and the verification by simulations will also be presented. 4.6 The Status of L1 Halo Orbit In order to further reveal the dynamical mechanism of error propagation in the libration point orbits, we present here a large retrograde elliptical orbit (corresponding to a remote Earth satellite) that is similar to the above-mentioned L1 halo orbit. The initial position, velocity and the corresponding orbital elements at the initial epoch 2016-09-30UTC0:00:0.0 are listed in Table 10 and Table 11, respectively. The initial period is 12.208831 d. To make the 27 d prediction on this orbit, the calculated results are listed in Table 12 and Table 13. Table 10 The position and velocity of the L1 halo orbit x/km
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
–345307.466263
21904.224767
19425.079239
–0.077571217
0.681599717
–0.223818048
Table 11 The orbital elements of the L1 halo orbit a/km
e
i/◦
Ω/◦
ω/◦
M/◦
223967.268
0.56141578
160.912742
160.912742
356.972610
356.972610
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
289
Table 12 The state of the L1 halo orbit propagated for 27 days model
a/km
e
i/◦
Ω/◦
λ/◦
standard
213214.203
0.653838
162.168040
350.586870
267.258205
Δa0 = 10 m
213214.209
0.653838
162.168042
347.586870
267.258158
Table 13 The position state of the L1 halo orbit propagated for 27 days model
x/km
y/km
z/km
r/km
standard
–229245.266
186181.824
–47023.643
299045.626
Δa0 = 10 m
–229245.393
186181.802
–47023.625
299045.707
Such a large elliptical retrograde orbit, in difference from the L1 halo orbit, is essentially an ordinary Kepler orbit. The inherent instability is not manifested in an arc of 27 d (only 2 periods). We present this trivial (in a sense) result to show from an other angle the inherent instability of the rapid error propagation of the libration point orbits in the Earth-Moon system. The motion and dynamics of the retrograde orbit however are beyond the scope of this paper. Readers can refer to Ref.[4] for details. 4.7
Determination of Libration Point Orbits and Verification of Short-arc Pre-
diction Accuracy Since the design of orbits strictly on the Lagrangian points is impossible in reality, the positioning of the probe on these points is only an approximation. The orbit maintenance must be continuously performed during operation, therefore we should consider only the short-arc problem for both the ground TT&C and on-orbit maneuver. Compared to the LEO and HEO satellites, the problems and difficulties met in the short-arc orbit determination and prediction of libration point orbits are the same. Using the orbit determination software developed by the Institute of Space Environment and Astrodynamics of Nanjing University and the USB (Unified S-band) observation data in China, without any other extra information, we have determined the orbit of the probe of the Chang’e 3 mission. The results are verified by comparing them with the post hoc data released by the Beijing Aerospace Flight Control Center (BAFCC). Based on this orbit determination, a very simple numerical extrapolation is adopted to predict the orbit. In the calculation, only the mass-point gravitations from the Earth, Moon and Sun, and a simple model of radiation pressure are taken into account, and the orbit determination given by the relevant mission is directly adopted as the initial orbital elements. A rather high accuracy is achieved without any difficulties. Omitting the unnecessary details, we list the results in Table 14 and Table 15. In Table 14 and Table 15, the orbit determinations A and B represent the results from the BAFCC and our calculations, respectively. Although the results in these tables
290
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
are basically able to explain the problem, we would like to give some further explanations as follows for a clearer understanding on the accuracies of the orbit determination and prediction for this kind of probes. (1) For the radiation pressure, in the condition without the specific details of the probe, we empirically use an equivalent planar model with an area-mass ratio obtained from the related estimates in our independent orbit determinations. (2) Although the details of both orbit determinations (including the observational data) are not specified, the calculation results in Tables 14–15 have demonstrated the accuracies of the orbit determination and prediction for such kind of orbits in this paper. In this two tables, the results of the 3 d and 7 d orbit predictions are compared with the results of post hoc precise orbit determinations, and the differences between the results of both orbit determinations are basically within 500 m. Such a comparison proves more strongly the reliability of the results of both orbit determinations and the rationality of the force model adopted by this paper. Table 14 The comparisons of two precise orbit determinations with the result of 3 d orbit prediction for the L2 halo orbit method
x/km
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
orbit determination A
–451897.422
138113.307
36463.951
–0.237782
–0.978285
–0.204902
orbit determination B
–451897.388
138112.927
36463.914
–0.237782
–0.978287
–0.204900
short arc propagation A
–451897.458
138113.145
36463.888
–0.237783
–0.978286
–0.204902
short arc propagation B
–451897.346
138112.816
36463.866
–0.237782
–0.978288
–0.204900
Table 15 The comparisons of two precise orbit determinations with the result of 7 d orbit prediction for the L2 halo orbit method
x/km
y/km
z/km
x/(km/s) ˙
y/(km/s) ˙
z/(km/s) ˙
orbit determination A
–371590.823
–223419.646
–55369.560
0.740132
–0.975794
–0.327870
orbit determination B
–371590.387
–223420.334
–55369.450
0.740133
–0.975794
–0.327871
short arc propagation A
–371590.759
–223420.292
–55369.879
0.740132
–0.975795
–0.327871
short arc propagation B
–371589.646
–223420.816
–55369.813
0.740138
–0.975793
–0.327873
The above calculations and the supplementary explanations indicate thoroughly that no special difficulties exist in the short-arc orbit determination and high-precision orbit prediction for this kind of libration point orbits, although the propagation of initial errors in these special orbits is much more significant than that of usual encircling orbits.
LIU Lin et al. / Chinese Astronomy and Astrophysics 43 (2019) 278–291
291
4.8 Transformation of Libration Point Orbit between Two Coordinate Systems The J2000 geocentric celestial coordinate system is required for orbit determination and prediction, but for the purpose of orbital maneuver the Earth-Moon rotating coordinate system is needed, bring up the necessity of coordinate conversion between these two systems. This conversion is easy to implement. For a given space mission, a specific rotating coordinate system that is most suitable for the mission may be defined according to actual needs, we will not make further discussions. References 1
Liu L., Tang J. S., Theory and Application of Satellite Orbits, Beijing: Publishing House of Electronics Industry, 2015
2
Gomez G., Llibre J., Martinez R., et al., Dynamics and Mission Design Near Libration Points. Sin-
3
Liu L., Hou X. Y., Theory and Application of Deep Space Exploring Orbits, Beijing: Publishing
4
Zhao C. Y., Liu L., AcASn, 1994, 35, 434
gapore, New Jersey, London, Hong Kong: World Scientific, 2001 House of Electronics Industry, 2015