# Evolution of the Selenopotential Model and Its Effects on the Propagation Accuracy of Orbits around the Moon

Authors: Young-Joo Song; Young-Rok Kim; Jonghee Bae; Bang-Yeop Kim
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

The current work analyzes the effect of applying different selenopotential models to the propagation of a lunar orbiting spacecraft. A brief evolution history of the selenopotential model is first presented; then, four representative selenopotential models are selected for force modeling. Expected propagation errors are presented with respect to three different circular polar orbits around the Moon. As a result, an expected but rather significant number of orbit propagation errors are discovered. Compared to the solutions obtained using the GRAIL1500E model, the overall 3D propagation errors for a 4-day period could reach up to several tens of kilometers (50 km altitude case with the GLGM2 model) and up to several hundreds of meters (50, 100, and 200 km altitude cases even with the GRAIL660B model). For each different orbiter’s altitude, the appropriate ranges of the degree and order of the gravitational harmonic coefficients are also suggested to yield the best propagation performances with respect to the performance obtained with the full harmonic coefficients using the GRAIL1500E model. The results of the current work are expected to serve as practical guidelines for the field of system budget analysis, mission design, mission operations, and the analysis of scientific results.


---

## Body

## 1. Introduction

An accurate understanding of the selenopotential model is very important not only for lunar scientific research but also for lunar spacecraft applications. Diverse areas of scientific study can be benefited, including the origin and evolution of the Moon, lunar interior physics, topographic mapping, and spacecraft flight dynamics in the vicinity of the Moon. The necessity of obtaining precise knowledge regarding the lunar gravity field was first raised during the Apollo era, with respect to bringing a man to the lunar surface and returning him safely to Earth, a purpose beyond purely theoretical scientific research [1]. Nevertheless, the associated outcomes have yielded numerous important scientific discoveries, such as the discovery of mascons (mass concentration) [2] and the discovery of lateral density variations inside the Moon, among others [1].Regarding more engineering-oriented aspects, knowledge of the lunar gravity field directly affects spacecraft orbit determination (OD) and orbit propagation (OP) knowledge. OD and OP knowledge are again closely related to the orbit prediction (OPD) performance, which considers OD uncertainties. Namely, the initial mean and covariance values of an orbit can be obtained based on the OD, and OPD and covariance prediction (CPD) can be conducted by considering the OP. For better OD knowledge, recent lunar exploration missions have devoted time and effort to developing a mission-oriented lunar gravity field model to achieve the best OD performance of each mission before the launch of the respective spacecraft, as a better understanding of the force model will eventually lead to better OD of the spacecraft. For example, a new global lunar gravity field model was generated during preparation of the Japan Aerospace Exploration Agency (JAXA) Selenological and Engineering Explorer (SELENE) mission [3], which was launched in 2007. Additionally, a lunar gravity model was developed for the National Aeronautics and Space Administration (NASA) Lunar Reconnaissance Orbiter (LRO) mission [4], launched in 2009. Recently, lunar mission complexity has increased; therefore, stricter OPD and OD requirements have been demanded, which again highlights the importance of adapting an improved lunar gravity field during the force modeling of a spacecraft. Indeed, the updated lunar gravity field model obtained from NASA’s Gravity Recovery and Interior Laboratory (GRAIL) mission, launched in 2011, improved the accuracy of the LRO’s precision OD results [5]. Information related to OPD directly affects mission planning as well as spacecraft pointing performances. From a mission planning point of view, an improved understanding of the lunar gravity field directly translates into accurate propellant gaging, which is used to control the orbit, and thus lower mission cost [1]. In addition, more precise contact time prediction ranging from an Earth-based ground station can be achieved. Regarding the spacecraft pointing knowledge, accurate OPD knowledge obtained from a gravity field is also an important factor; it may be used to establish the accuracy level of pointing knowledge needed by onboard instruments used for gravity field studies or for other mapping or scientific tasks to achieve the final mission goal. In addition to OPD and OD performance enhancement, a more critical application that requires better knowledge of force modeling is pinpoint unguided surface landing missions [6–9], as knowledge of the parking orbit is more critical than that of the gravity field perturbations induced during a relatively short descent phase [1].Since 1992, Korea has continuously operated more than ten Earth-orbiting satellites. Korea is now beginning to expand its interests to planetary missions. The Korean space program plans on launching a lunar orbiter and lander and on exploring Mars, asteroids, and deep space in the future. The first phase of the Korean Lunar Exploration Program (KLEP), announced in 2013 [10], has recently been initiated. For the first step, the Korea Pathfinder Lunar Orbiter (KPLO), an experimental lunar orbiter that will orbit the Moon approximately 100 km above the lunar surface for one year, is tentatively scheduled for launch, through an international collaboration, in or shortly after 2018. For the KPLO mission, a total of five instruments (four from Korea and one from NASA) will be installed on board to complete the primary mission objectives, including verification of the disruption tolerant network (DTN) methodologies. The KPLO instruments and developers are as follows: Lunar Terrain Imager (LUTI) by the Korea Aerospace Research Institute (KARI), Polarimetric Camera (PolCam) by the Korea Astronomy and Space Science Institute (KASI), KPLO Gamma Ray Spectrometer (KGRS) by the Korea Institute of Geoscience and Mineral Resources (KIGAM), KPLO Magnetometer (KMAG) by Kyung Hee University (KHU), and, finally, Shadow Camera (ShadowCam) by investigators at Arizona State University (ASU) and Malin Space Science Systems (MSSS) as a NASA contribution. Additionally, the Electronics and Telecommunications Research Institute (ETRI) of Korea and NASA are working together to develop and demonstrate DTN technologies.To ensure the successful launch of the KPLO, the KARI is now performing extensive prephase works, including joint efforts between KARI and NASA to validate the trajectory design and navigation performances. From the OD perspective, analyses of the following topics are ongoing to meet the given mission requirements: required minimum tracking arc, number of required ground sites and their locations, required media correction accuracy for ground OD system, minimal OD processing time with respect to overall ground operation flows, and so on. Additionally, expected OD performances, which depend on the lunar gravity field model selected with consideration of different levels of the degree and order of the gravitational harmonic coefficients, are being investigated. Regarding the OPD, not only ground OPD but also the KPLO onboard OPD requirements and expected errors are being extensively analyzed, investigated, and determined. During the analysis of OPD performance requirements, the issue of OP knowledge has been raised with respect to the successful operation of the onboard payloads. This issue is again related to the selection of an appropriate lunar gravity field model with an associated level of the degree and order of the gravitational harmonic coefficients. OP performances cannot be directly estimated by applying only the complete lunar gravity field model recently published because numerous tradeoff studies should be performed to satisfy mission requirements, which may slightly differ for each mission operation phase, with less burden on not only the spacecraft bus system but also the ground system.Several studies have already emphasized the importance of selecting an appropriate selenopotential model to improve OD results during the real flight operation of a lunar orbiter [5, 11–14]. Previous studies combined OD, OPD, and OP simultaneously during their analyses, an approach that is certainly realistic with respect to the mission operation phase. However, during the early system design and analysis phase, obtaining direct knowledge regarding the expected OP error is quite important, and direct insights can hardly be achieved via the approach used in previous studies. Therefore, the current work is more focused on obtaining direct knowledge related to OP performance. As already addressed, obtaining insights into the expected OP errors due to the selection of different lunar gravity field models is potentially important for a wide range of areas, particularly when establishing a spacecraft’s pointing budgets, which is closely related to the dynamic modeling errors. More efficient budget analysis can be conducted when the allowable uncertainties between the attitude control (AC), attitude determination (AD), and ground processed OD errors are more effectively shared, adjusted, and determined as a result of more knowledge regarding the expected OP errors. Those of efficient tradeoffs may finally lead not only to proper selection of the spacecraft hardware components but also to realization of the level of accuracy of the ground system.The main goal of the current work is to provide practical insights into the expected OP errors resulting from the selection of a lunar gravity field model, with the inclusion of a very recently released lunar gravity field model, and to provide guidelines for that selection. The results of this work can be used to estimate the dynamic model errors applied during OP. For a complete analysis of OPD knowledge, further analysis, including covariance analysis, should be performed using the initial OD errors, with more details regarding the dynamic model errors. The current work briefly presents the evolution history of the selenopotential model and analyzes its effects on a spacecraft’s OP accuracy, especially spacecraft orbiting the Moon at different altitudes with a circular polar orbit. Through this work, practical insights into the propagation errors of expected states can be obtained in cases where different selenopotential models are applied. Additionally, consideration of the minimum ranges of the degree and order of gravitational harmonic coefficients necessary to match the OP performance simulated using the full degrees and orders of gravitational harmonics is suggested with respect to different altitudes around the Moon. These insights will eventually aid in the establishment of an OPD and OD analysis strategy. In Section2, a brief history of the evolution of the lunar gravity field solution from the early 1960s to the present day is provided, and the associated characteristics are reviewed. The assumptions made and data used for the current simulations are described in Section 3, along with a brief summary of the method used to implement high-degree and high-order gravity fields. The analysis results of the expected OP errors when different selenopotential models are applied to a lunar orbiting spacecraft flying at different altitudes are presented in Section 4. Finally, the conclusions are presented in Section 5. Although the current work was conducted as a part of the KPLO system design phase activities, the results certainly offer various practical insights into the field of system budget analysis, mission design, mission operations, and the analysis of scientific results.


## 2. Evolution of the Selenopotential Model

### 2.1. 1960s

The first study of the lunar gravity field was initiated when the Luna 10 Russian spacecraft launched in 1966. The initial distance of Luna 10 from the lunar surface was 350 km at perilune and 1015 km at apolune, and it transmitted its own signal until May 30, 1966. By using the Luna 10 spacecraft data, the first dynamical proof of the Moon’s oblateness was made, and a realistic oblateness coefficient of the Moon was obtained [15]. In August 1966, the first US lunar orbiter I (LO-I) was launched; four additional orbiters (LO-II, III, IV, and V) were placed in lunar orbit, with various inclinations and eccentricities and periapses of 50 to 100 km above the lunar surface, by August 1967. Based on the LO data, the spherical harmonics of the Moon (degree 4 in the tesserals and degree 8 in the zonals) were first published in 1968 [16]; later that year, the existence of large mascons under the lunar ringed maria was also discovered [2].In addition to the LO data, the Apollo program (initiated in 1963 and ended in 1972) also provided a vast amount of various types of data that could be used to model and to determine the structure of the Moon. Unlike LOs, the Apollo spacecraft was placed in near-circular orbits, although some tracking, from altitudes as low as 10 to 20 km, was acquired, with a mean altitude of 100 km having low inclinations. During the 1970s, diverse results regarding new lunar gravity solutions were announced by utilizing the Apollo data combined with various LO radio tracking data.


### 2.2. 1970s

In 1971, the 15th-degree and 8th-order global lunar gravity field solution was determined [17] using the tracking data of five LOs; also, a surface-layer representation of the lunar gravitational field was derived [18]. The 13th-degree and 13th-order spherical harmonic coefficients of the lunar gravitational field were determined in 1972, also by utilizing the tracking data from LOs [19].In addition to radio tracking data, Lunar Laser Ranging (LLR) has provided information regarding the lunar gravity field. For LLR experiments, retroreflectors were installed on the Moon during the Apollo program (11, 14, and 15) and the two Lunokhod missions of the Soviet Union [20]. LLR solutions related to lunar gravity field study became available in the early 1970s and were mainly used to improve the estimation of second-degree harmonics [1]. In 1973, the three principal lunar moments of inertia were determined based on LLR by utilizing the retroreflectors installed during the Apollo and Lunokhod-2 missions [21]. Lunar gravity solutions, especially the acceleration profile of mascons, were obtained from Apollo flight data until the mid-1970s [22–26]. Though the data from the Apollo spacecraft were short-arc and thus mainly used to produce the acceleration profile, they were also implemented in point-mass and surface-layer solutions [1]. Later, in 1977, a point-mass representation of the quasi-global gravity field of the Moon, consisting of 117 distributed point masses, was developed by processing the Apollo 15 and 16 subsatellite and LO-V radio tracking data [27]. A 16th-degree and 16th-order spherical harmonic lunar gravity field was also derived in 1977 based on the long-term Keplerian variations in the orbits of the Apollo subsatellites and LO-V [28]. This model was the first model that showed plausible correlations between far-side gravity anomalies and near-side surface features [28, 29].


### 2.3. 1980s to 1990s

In 1980, a 16th-degree spherical harmonic solution was obtained by combining most of the LO, Apollo, Apollo subsatellite, and LLR data available at the time [30]; this model was the best lunar gravity model for the next 13 years. A continuation of the LLR led to an improvement in accuracy of the lunar ephemeris of three orders of magnitude, an improvement in the measurement of the variations of the Moon’s rotation of several orders of magnitude, and the discovery of additional information regarding the Moon’s tidal acceleration and rotational dissipation [31]. New geodetic theory and computer technology development enabled the release of the 60th-degree and 60th-order lunar gravity field model, the Lun60D, based on the tracking of the LOs and the Apollo subsatellites [32]. This model has excellent performance when mapping the gravity anomalies on the lunar surface from an altitude of 100 km [33]. However, when the gravity anomalies are mapped onto the lunar surface using the Lun60D model, numerous artifacts become apparent because the high-degree terms have excessive power, thus making the model unsuitable for use in geophysical studies [33]. On January 24, 1994, the Clementine spacecraft was launched from Vandenberg Air Force Base and became the first US spacecraft to return to lunar orbit in almost 20 years [34]. Using the Clementine tracking data in addition to the same historic LO and Apollo data, a new lunar gravity field model of the 70th degree and order was developed and named the Goddard Lunar Gravity Model-1 (GLGM-1) [35]. However, GLGM-1 excluded the Apollo 16 subsatellite data, and the solutions were not calibrated [33]. In 1997, a final calibrated solution with the inclusion of the Apollo 16 subsatellite data was obtained, that is, GLGM-2 [33]. GLGM-2 provided an improvement in the low degree and sectoral terms (to the 20th degree) of the lunar gravity field. On January 6, 1998, NASA’s third lunar discovery mission, Lunar Prospector (LP), was launched. A series of maneuvers placed LP in a near-circular orbit at an altitude of 100 km with a 2-hour orbital period [36]. The first lunar gravity models obtained with the inclusion of the LP data were the 75th-degree and 75th-order models, that is, LP75D and LP75G [37], followed by the 100th-degree and 100th-order models, that is, LP100J and LP100K [38], which were obtained through inclusion of the extended mission data of the LP. All the LP gravity models also include all the data available from the previous LO I–V, Apollo 15 and 16, and Clementine missions. During that time, the LP100J and LP100K models were expected to provide the best OD accuracy versus computational time required to determine the orbits and for the initial operational use of future missions.


### 2.4. 2000s to Present Day

In 2001, the 165th-degree and 165th-order lunar gravity model, LP165P, which used data from all NASA orbiter missions existing at the time, including the LP data, was released [39]. The LP165P model described several large‐scale anomalies in the far-side low-latitude region [39]. Several years later, an improved version of LP165P, though limited to the 150th degree and 150th order, was released. This model, called the LP150Q spherical harmonic model, was considered the best available lunar gravity field model at that time [40]. Even though the LP150Q model provided high-resolution near-side gravity maps, it still suffered from poor resolution of the far-side solutions. The SELENE, launched in 2007, provided the first global data set of the Moon through the use of three satellites, a main orbiter and two relay subsatellites, and improved far-side gravity map solutions. The first lunar gravity field model that utilized the SELENE tracking data as well as the historic data was the 90th-degree and 90th-order SELENE Gravity Model 90d (SGM90d), presented in 2009 [41], followed the models of 100th and 150th degree and order (SGM100h, SGM100i, and SGM150) [42, 43]. In 2007, the first Chinese lunar exploration satellite, Chang’E-1, was launched and provided a lunar gravity model, that is, CEGM02, with a degree and order up to 100 [44]. Compared to the LP series model, the SGM series showed a significant improvement in the precision OD for the lunar far side and presented new gravity anomaly features of the lunar far side.To prepare for future human exploration of the Moon, that is, to identify safe landing sites, locate potential resources on the Moon, characterize the radiation environment, and demonstrate new technologies, NASA launched the LRO in June 2009. To show the best performances of OPD and OD for the LRO mission, NASA released a new gravity field model called GLGM-3, which included different strategies for data weighting and a regularization method that used the same LP150Q model. The performance of GLGM-3 was very similar to that of the LP150Q model, particularly with respect to orbital performance [4].In 2011, the GRAIL orbiter was launched. An unprecedented resolution of the lunar gravity field, achieved by placing two spacecrafts into the same low orbit around the Moon, was expected. Before the launch of the GRAIL mission, the SGM150 model was widely utilized for the OPD and OD analyses of the LRO mission. Shortly after the launch, two series of lunar gravity models designed for the LRO mission, that is, LRO Lunar Gravity Models 1 and 2 (LLGM-1 and LLGM-2), were released [11, 45]. The LLGM models were 150th-degree and 150th-order models derived based on the LO, Apollo subsatellite, LP, Clementine, and LRO’s Doppler and Lunar Orbiter Laser Altimeter (LOLA) tracking data. Almost concurrently with the release of the LLGM-2 model, the first lunar gravity model from the GRAIL mission was released. Meanwhile, the accuracy of the following lunar parameters was continuously improved by obtaining the LLR measurements: solid moment of inertia, fluid core moment of inertia, core oblateness, Love number determination, orbit and physical librations, inner core possibilities, and so on [46]. The first lunar gravity field model that utilized the GRAIL mission data had a degree and order of 150 and 270, respectively [12]. To determine the lunar gravity field based on the GRAIL mission data, independent but collaborative groups, one at the Jet Propulsion Laboratory (JPL) and the other at the Goddard Space Flight Center (GSFC), were involved. They differed in terms of the OD software, a priori models, data editing, and parameter estimation strategy used, among others. Later, global field models that utilized the GRAIL data, with a maximum degree and order of 420 [47], 660 [48, 49], 900 [50, 51], and 1200 [52], were released. Very recently, the highest resolution lunar gravity field model to date, with a maximum degree and order of 1500, was generated by analyzing the data from the primary and extended GRAIL missions, yielding a surface resolution of 3.6 km [53].


## 2.1. 1960s

The first study of the lunar gravity field was initiated when the Luna 10 Russian spacecraft launched in 1966. The initial distance of Luna 10 from the lunar surface was 350 km at perilune and 1015 km at apolune, and it transmitted its own signal until May 30, 1966. By using the Luna 10 spacecraft data, the first dynamical proof of the Moon’s oblateness was made, and a realistic oblateness coefficient of the Moon was obtained [15]. In August 1966, the first US lunar orbiter I (LO-I) was launched; four additional orbiters (LO-II, III, IV, and V) were placed in lunar orbit, with various inclinations and eccentricities and periapses of 50 to 100 km above the lunar surface, by August 1967. Based on the LO data, the spherical harmonics of the Moon (degree 4 in the tesserals and degree 8 in the zonals) were first published in 1968 [16]; later that year, the existence of large mascons under the lunar ringed maria was also discovered [2].In addition to the LO data, the Apollo program (initiated in 1963 and ended in 1972) also provided a vast amount of various types of data that could be used to model and to determine the structure of the Moon. Unlike LOs, the Apollo spacecraft was placed in near-circular orbits, although some tracking, from altitudes as low as 10 to 20 km, was acquired, with a mean altitude of 100 km having low inclinations. During the 1970s, diverse results regarding new lunar gravity solutions were announced by utilizing the Apollo data combined with various LO radio tracking data.


## 2.2. 1970s

In 1971, the 15th-degree and 8th-order global lunar gravity field solution was determined [17] using the tracking data of five LOs; also, a surface-layer representation of the lunar gravitational field was derived [18]. The 13th-degree and 13th-order spherical harmonic coefficients of the lunar gravitational field were determined in 1972, also by utilizing the tracking data from LOs [19].In addition to radio tracking data, Lunar Laser Ranging (LLR) has provided information regarding the lunar gravity field. For LLR experiments, retroreflectors were installed on the Moon during the Apollo program (11, 14, and 15) and the two Lunokhod missions of the Soviet Union [20]. LLR solutions related to lunar gravity field study became available in the early 1970s and were mainly used to improve the estimation of second-degree harmonics [1]. In 1973, the three principal lunar moments of inertia were determined based on LLR by utilizing the retroreflectors installed during the Apollo and Lunokhod-2 missions [21]. Lunar gravity solutions, especially the acceleration profile of mascons, were obtained from Apollo flight data until the mid-1970s [22–26]. Though the data from the Apollo spacecraft were short-arc and thus mainly used to produce the acceleration profile, they were also implemented in point-mass and surface-layer solutions [1]. Later, in 1977, a point-mass representation of the quasi-global gravity field of the Moon, consisting of 117 distributed point masses, was developed by processing the Apollo 15 and 16 subsatellite and LO-V radio tracking data [27]. A 16th-degree and 16th-order spherical harmonic lunar gravity field was also derived in 1977 based on the long-term Keplerian variations in the orbits of the Apollo subsatellites and LO-V [28]. This model was the first model that showed plausible correlations between far-side gravity anomalies and near-side surface features [28, 29].


## 2.3. 1980s to 1990s

In 1980, a 16th-degree spherical harmonic solution was obtained by combining most of the LO, Apollo, Apollo subsatellite, and LLR data available at the time [30]; this model was the best lunar gravity model for the next 13 years. A continuation of the LLR led to an improvement in accuracy of the lunar ephemeris of three orders of magnitude, an improvement in the measurement of the variations of the Moon’s rotation of several orders of magnitude, and the discovery of additional information regarding the Moon’s tidal acceleration and rotational dissipation [31]. New geodetic theory and computer technology development enabled the release of the 60th-degree and 60th-order lunar gravity field model, the Lun60D, based on the tracking of the LOs and the Apollo subsatellites [32]. This model has excellent performance when mapping the gravity anomalies on the lunar surface from an altitude of 100 km [33]. However, when the gravity anomalies are mapped onto the lunar surface using the Lun60D model, numerous artifacts become apparent because the high-degree terms have excessive power, thus making the model unsuitable for use in geophysical studies [33]. On January 24, 1994, the Clementine spacecraft was launched from Vandenberg Air Force Base and became the first US spacecraft to return to lunar orbit in almost 20 years [34]. Using the Clementine tracking data in addition to the same historic LO and Apollo data, a new lunar gravity field model of the 70th degree and order was developed and named the Goddard Lunar Gravity Model-1 (GLGM-1) [35]. However, GLGM-1 excluded the Apollo 16 subsatellite data, and the solutions were not calibrated [33]. In 1997, a final calibrated solution with the inclusion of the Apollo 16 subsatellite data was obtained, that is, GLGM-2 [33]. GLGM-2 provided an improvement in the low degree and sectoral terms (to the 20th degree) of the lunar gravity field. On January 6, 1998, NASA’s third lunar discovery mission, Lunar Prospector (LP), was launched. A series of maneuvers placed LP in a near-circular orbit at an altitude of 100 km with a 2-hour orbital period [36]. The first lunar gravity models obtained with the inclusion of the LP data were the 75th-degree and 75th-order models, that is, LP75D and LP75G [37], followed by the 100th-degree and 100th-order models, that is, LP100J and LP100K [38], which were obtained through inclusion of the extended mission data of the LP. All the LP gravity models also include all the data available from the previous LO I–V, Apollo 15 and 16, and Clementine missions. During that time, the LP100J and LP100K models were expected to provide the best OD accuracy versus computational time required to determine the orbits and for the initial operational use of future missions.


## 2.4. 2000s to Present Day

In 2001, the 165th-degree and 165th-order lunar gravity model, LP165P, which used data from all NASA orbiter missions existing at the time, including the LP data, was released [39]. The LP165P model described several large‐scale anomalies in the far-side low-latitude region [39]. Several years later, an improved version of LP165P, though limited to the 150th degree and 150th order, was released. This model, called the LP150Q spherical harmonic model, was considered the best available lunar gravity field model at that time [40]. Even though the LP150Q model provided high-resolution near-side gravity maps, it still suffered from poor resolution of the far-side solutions. The SELENE, launched in 2007, provided the first global data set of the Moon through the use of three satellites, a main orbiter and two relay subsatellites, and improved far-side gravity map solutions. The first lunar gravity field model that utilized the SELENE tracking data as well as the historic data was the 90th-degree and 90th-order SELENE Gravity Model 90d (SGM90d), presented in 2009 [41], followed the models of 100th and 150th degree and order (SGM100h, SGM100i, and SGM150) [42, 43]. In 2007, the first Chinese lunar exploration satellite, Chang’E-1, was launched and provided a lunar gravity model, that is, CEGM02, with a degree and order up to 100 [44]. Compared to the LP series model, the SGM series showed a significant improvement in the precision OD for the lunar far side and presented new gravity anomaly features of the lunar far side.To prepare for future human exploration of the Moon, that is, to identify safe landing sites, locate potential resources on the Moon, characterize the radiation environment, and demonstrate new technologies, NASA launched the LRO in June 2009. To show the best performances of OPD and OD for the LRO mission, NASA released a new gravity field model called GLGM-3, which included different strategies for data weighting and a regularization method that used the same LP150Q model. The performance of GLGM-3 was very similar to that of the LP150Q model, particularly with respect to orbital performance [4].In 2011, the GRAIL orbiter was launched. An unprecedented resolution of the lunar gravity field, achieved by placing two spacecrafts into the same low orbit around the Moon, was expected. Before the launch of the GRAIL mission, the SGM150 model was widely utilized for the OPD and OD analyses of the LRO mission. Shortly after the launch, two series of lunar gravity models designed for the LRO mission, that is, LRO Lunar Gravity Models 1 and 2 (LLGM-1 and LLGM-2), were released [11, 45]. The LLGM models were 150th-degree and 150th-order models derived based on the LO, Apollo subsatellite, LP, Clementine, and LRO’s Doppler and Lunar Orbiter Laser Altimeter (LOLA) tracking data. Almost concurrently with the release of the LLGM-2 model, the first lunar gravity model from the GRAIL mission was released. Meanwhile, the accuracy of the following lunar parameters was continuously improved by obtaining the LLR measurements: solid moment of inertia, fluid core moment of inertia, core oblateness, Love number determination, orbit and physical librations, inner core possibilities, and so on [46]. The first lunar gravity field model that utilized the GRAIL mission data had a degree and order of 150 and 270, respectively [12]. To determine the lunar gravity field based on the GRAIL mission data, independent but collaborative groups, one at the Jet Propulsion Laboratory (JPL) and the other at the Goddard Space Flight Center (GSFC), were involved. They differed in terms of the OD software, a priori models, data editing, and parameter estimation strategy used, among others. Later, global field models that utilized the GRAIL data, with a maximum degree and order of 420 [47], 660 [48, 49], 900 [50, 51], and 1200 [52], were released. Very recently, the highest resolution lunar gravity field model to date, with a maximum degree and order of 1500, was generated by analyzing the data from the primary and extended GRAIL missions, yielding a surface resolution of 3.6 km [53].


## 3. Simulation Setup

This section describes the method used throughout the current work while analyzing the effect of applying different selenopotential models on the propagation of a lunar orbiting spacecraft. First, the force modeling method used to implement high-degree and high-order gravity fields is briefly summarized. Then, several assumptions made for the current simulation are explained, including the shortcomings, numerical implications, and data used.


### 3.1. Force Modeling

The conventional perturbing acceleration acting on a spacecraft due to the nonspherical shape of the Moon can be described by a spherical-potential function,U, as shown in (1)Ur,φ,λ=GMr∑n=0∞∑m=0nRMrnC-nmcos⁡mλ+S-nmsin⁡mλP-nmsin⁡φ,where r is the distance between a spacecraft and the center of the Moon, φ and λ are the latitude and the longitude of a spacecraft, respectively, with respect to the principal axis (PA) reference system of the Moon, G is the universal gravitational constant of the Moon, M is the Moon mass, RM is the Moon equatorial radius, and P¯nm denotes the fully normalized Legendre polynomials with degree n and order m, computed based on the recursion relationship. In addition, C¯nm, S¯nm are the normalized degree and order harmonic coefficients, respectively. The spherical-potential function U shown in (1) can be rewritten with defined “lumped coefficients,” as shown in [54–56](2)Ur,φ,λ=GMr∑m=0NAm0cos⁡mλ+Bm0sin⁡mλ.The defined “lumped coefficients”Am(0) and Bm(0) in (2) can be expressed as follows [56]:(3a)Am0=∑n=mNRMrnC-nmP-nmφ,(3b)Bm0=∑n=mNRMrnS-nmP-nmφ,where N is the maximum finite degree of the spherical harmonic expansion [56]. More details regarding the force modeling method, for example, those related to implementation of the recursion algorithm and determination of the inertial acceleration due to the nonsphericity of the central body, expressed in body-fixed Cartesian coordinates, can be found in [54–57].


### 3.2. Data and Assumptions

To investigate the effect of applying different lunar gravity field models to determine the performance of the OP, four different representative selenopotential models, that is, GLGM2, LP150Q, GRAIL660B, and GRAIL1500E, are selected. GLGM2 was released in the late 1990s, LP150Q was released in the early 2000s, GRAIL660B was released in 2013, and GRAIL1500E, the most recent gravity model, was released in 2015. In the following discussions, the OP results simulated using the GRAIL1500E model, with the full degree and order of the harmonic coefficients, are always regarded as references, that is, the truth values, when comparing the OP performances. All these models are all available from NASA’s Planetary Data System (PDS) Geoscience Node [58]. For the initial orbital conditions, a fictitious lunar orbiter is assumed to orbit around the Moon with a circular (0 eccentricity) polar orbit (90 deg inclination) and at three different altitudes: 50, 100, and 200 km. Each orbiter is propagated for four days to allow time for uploading operation schedules during the weekends. Because the effect of applying different selenopotential models on the OP performance may cause the data to differ from that of the spacecraft’s actual lunar surface passage, five different Right Ascension of Ascending Node (RAAN) initial orbital element values are used: 0, 30, 60, 90, and 120 deg. With each of the five different RAANs, the spacecraft’s lunar ground track covers almost the entire lunar surface during each of the four days of propagation. The other remaining orbital elements, such as the Argument of Periapsis (AoP) and Argument of Latitude (AoL), are set at 0 deg.For a numerical integration, the Runge-Kutta-Fehlberg 7th-8th-order variable step-size integrator, with a truncation error tolerance ofε=1×10-13, is used. To accurately derive the ephemerides of the planets as well as all planetary constants, four different JPL ephemeris values are adapted for each different selenopotential model. The applied ephemerides and planetary constants are shown in Table 1. Additionally, the lunar orientation of the body-fixed coordinate frame specified by the JPL ephemeris series is used for each different selenopotential model for high-precision work. The initial orbital epoch is assumed to be January 01, 2020, 00:00:00 (UTC), corresponding to Korea’s first experimental lunar orbiter mission, the launch and subsequent nominal mission of which are set to occur in approximately 2019~2020. To numerically integrate the orbiter’s states, the proven performance of KARI’s lunar orbit propagator is used [57]. Note that external perturbing forces other than the force due to the lunar gravity field are not additionally considered to focus on the main objectives of the current study. As already discussed, the results presented in the current work cover only OP errors resulting from the selection of the selenopotential model. Therefore, the presented OP errors can be understood as errors due to incomplete dynamic modeling during OP and may differ from those of OPD performances, as OPD performance is closely related to OD knowledge.Table 1
Applied planetary ephemerides and constants with respect to different selenopotential models.
Selenopotential model
JPL Ephemeris used
Moon’s mean radius (km)
Moon’s gravitational constant (km3/s2)
GLGM2
DE403 [59]
1738.0 km
4902.8029500000
LP150Q
DE405 [60]
1738.0 km
4902.8010760000
GRAIL660B
DE421 [61]
1738.0 km
4902.8003055554
GRAIL1500E
DE430 [62]
1738.0 km
4902.8002221408


## 3.1. Force Modeling

The conventional perturbing acceleration acting on a spacecraft due to the nonspherical shape of the Moon can be described by a spherical-potential function,U, as shown in (1)Ur,φ,λ=GMr∑n=0∞∑m=0nRMrnC-nmcos⁡mλ+S-nmsin⁡mλP-nmsin⁡φ,where r is the distance between a spacecraft and the center of the Moon, φ and λ are the latitude and the longitude of a spacecraft, respectively, with respect to the principal axis (PA) reference system of the Moon, G is the universal gravitational constant of the Moon, M is the Moon mass, RM is the Moon equatorial radius, and P¯nm denotes the fully normalized Legendre polynomials with degree n and order m, computed based on the recursion relationship. In addition, C¯nm, S¯nm are the normalized degree and order harmonic coefficients, respectively. The spherical-potential function U shown in (1) can be rewritten with defined “lumped coefficients,” as shown in [54–56](2)Ur,φ,λ=GMr∑m=0NAm0cos⁡mλ+Bm0sin⁡mλ.The defined “lumped coefficients”Am(0) and Bm(0) in (2) can be expressed as follows [56]:(3a)Am0=∑n=mNRMrnC-nmP-nmφ,(3b)Bm0=∑n=mNRMrnS-nmP-nmφ,where N is the maximum finite degree of the spherical harmonic expansion [56]. More details regarding the force modeling method, for example, those related to implementation of the recursion algorithm and determination of the inertial acceleration due to the nonsphericity of the central body, expressed in body-fixed Cartesian coordinates, can be found in [54–57].


## 3.2. Data and Assumptions

To investigate the effect of applying different lunar gravity field models to determine the performance of the OP, four different representative selenopotential models, that is, GLGM2, LP150Q, GRAIL660B, and GRAIL1500E, are selected. GLGM2 was released in the late 1990s, LP150Q was released in the early 2000s, GRAIL660B was released in 2013, and GRAIL1500E, the most recent gravity model, was released in 2015. In the following discussions, the OP results simulated using the GRAIL1500E model, with the full degree and order of the harmonic coefficients, are always regarded as references, that is, the truth values, when comparing the OP performances. All these models are all available from NASA’s Planetary Data System (PDS) Geoscience Node [58]. For the initial orbital conditions, a fictitious lunar orbiter is assumed to orbit around the Moon with a circular (0 eccentricity) polar orbit (90 deg inclination) and at three different altitudes: 50, 100, and 200 km. Each orbiter is propagated for four days to allow time for uploading operation schedules during the weekends. Because the effect of applying different selenopotential models on the OP performance may cause the data to differ from that of the spacecraft’s actual lunar surface passage, five different Right Ascension of Ascending Node (RAAN) initial orbital element values are used: 0, 30, 60, 90, and 120 deg. With each of the five different RAANs, the spacecraft’s lunar ground track covers almost the entire lunar surface during each of the four days of propagation. The other remaining orbital elements, such as the Argument of Periapsis (AoP) and Argument of Latitude (AoL), are set at 0 deg.For a numerical integration, the Runge-Kutta-Fehlberg 7th-8th-order variable step-size integrator, with a truncation error tolerance ofε=1×10-13, is used. To accurately derive the ephemerides of the planets as well as all planetary constants, four different JPL ephemeris values are adapted for each different selenopotential model. The applied ephemerides and planetary constants are shown in Table 1. Additionally, the lunar orientation of the body-fixed coordinate frame specified by the JPL ephemeris series is used for each different selenopotential model for high-precision work. The initial orbital epoch is assumed to be January 01, 2020, 00:00:00 (UTC), corresponding to Korea’s first experimental lunar orbiter mission, the launch and subsequent nominal mission of which are set to occur in approximately 2019~2020. To numerically integrate the orbiter’s states, the proven performance of KARI’s lunar orbit propagator is used [57]. Note that external perturbing forces other than the force due to the lunar gravity field are not additionally considered to focus on the main objectives of the current study. As already discussed, the results presented in the current work cover only OP errors resulting from the selection of the selenopotential model. Therefore, the presented OP errors can be understood as errors due to incomplete dynamic modeling during OP and may differ from those of OPD performances, as OPD performance is closely related to OD knowledge.Table 1
Applied planetary ephemerides and constants with respect to different selenopotential models.
Selenopotential model
JPL Ephemeris used
Moon’s mean radius (km)
Moon’s gravitational constant (km3/s2)
GLGM2
DE403 [59]
1738.0 km
4902.8029500000
LP150Q
DE405 [60]
1738.0 km
4902.8010760000
GRAIL660B
DE421 [61]
1738.0 km
4902.8003055554
GRAIL1500E
DE430 [62]
1738.0 km
4902.8002221408


## 4. Results and Analysis

### 4.1. Orbit Propagation Performance Comparisons

In this section, the OP performances are compared, and the expected OP errors due to selecting different selenopotential models are investigated. To compare the OP performances, the full degree and order of the harmonic coefficients are considered during the simulations of each selected lunar gravity model. A 70th-degree and 70th-order are considered for GLGM2, 150th for LP150Q, 660th for GRAIL660B, and 1500th for the GRAIL1500E model. By regarding the simulated OP results of the GRAIL1500E model as references, Figures1–3 compare the 3D position differences when the same OP is used for different lunar gravity models. Figure 1 shows the comparison between GRAIL1500E and GLGM2, Figure 2 shows that between GRAIL1500E and LP150Q, and Figure 3 shows that between GRAIL1500E and GRAIL660B. Additionally, in each of the figures, subfigure (a) presents the comparison results for the 50 km spacecraft orbital altitude case, (b) presents the 100 km case, and (c) presents the 200 km case.Figure 1
OP performance comparison between the GRAIL1500E and GLGM2 models. Subfigure (a) presents the 50 km spacecraft altitude case, (b) presents the 100 km case, and (c) presents the 200 km case.
(a)
(b)
(c)Figure 2
OP performance comparison between the GRAIL1500E and LP150Q models. Subfigure (a) presents the 50 km spacecraft altitude case, (b) presents the 100 km case, and (c) presents the 200 km case.
(a)
(b)
(c)Figure 3
OP performance comparison between the GRAIL1500E and GRAIL660B models. Subfigure (a) presents the 50 km spacecraft altitude case, (b) presents the 100 km case, and (c) presents the 200 km case.
(a)
(b)
(c)Figures1–3 show that a significant number of OP errors are induced simply by considering different selenopotential models. Additionally, as expected, when the spacecraft altitude decreases, the number of associated OP errors tends to increase. For example, the comparison of the 3D positions of OP errors between the GRAIL1500E and GLGM2 models shows a maximum of approximately 44.233 km errors when the spacecraft’s altitude is 50 km, approximately 28.812 km errors when the altitude is 100 km, and approximately 11.724 km errors when the altitude is 200 km. In addition, the number of expected OP errors has different trends based on the given initial orbital element of RAAN. For the GLGM2 comparison case, the maximum number of OP errors was observed for a RAAN of 0 deg, as was already discussed.However, the minimum number of OP errors was observed for an initial RAAN of 90 deg, and 3D position errors of approximately 6.759 km were found. This is due to the uncertainties of the lunar gravity model present during the early phase of the OP time span. Indeed, if a RAAN of 0 deg is given, the initial location of the spacecraft is determined to be at the far side of the Moon; by contrast, a RAAN of 90 deg indicates an initial location at the near side of the Moon. These results indirectly show the modeling accuracy of lunar far-side gravity fields of the GLGM2 model. Nevertheless, the OP errors caused solely by applying different lunar gravity field models are found to be rather significant, even though GLGM2 is a very old lunar gravity field model (having been released in the late 1990s). For the comparison case with the LP150Q model, which is shown in Figure2, OP errors are still found within orders of several kilometers. The expected maximum 3D-position OP error for the four days of propagation is found to be approximately 6.516 km at an altitude of 50 km, 4.476 km at an altitude of 100 km, and 1.550 km at an altitude of 200 km. Unlike the comparison case with the GLGM2 model, the maximum 3D-position OP errors are found for a RAAN of 30 deg for the LP150Q model, with the spacecraft still initially located on the far side of the Moon. However, when the RAAN is 90 and 120 deg, indicating an initial location on the near side of the Moon, the 3D OP errors of the LP150Q model at an altitude of 50 km are approximately 113 m and 95 m, respectively. The OP comparison results for the case with the GRAIL660B model are presented in Figure 3. Unlike the other two cases (GLGM2 and LP150Q), the OP errors simulated using the GRAIL660B model remained within several hundreds of meters with very similar trends regardless of the altitude. Very minor differences due to the altitude are actually observed, but they are not as significant as those observed in the former two cases. Additionally, the initially given RAANs did not seriously affect the OP performances. This result may be due to the model resolution improvement (surface resolution ranging from 9.3 to 14.6 km for the GRAIL660B model [47] and of 3.6 km for GRAIL1500E [53]) rather than the far-side modeling accuracy improvement, as occurred between LP150Q/GLGM-2 and the GRAIL660B model.Because Figures1–3 only plot the total 3D position differences, details of the results, including the position components of the Radial-Transverse-Normal (RTN) direction, are provided separately in Tables 2–4. Table 2 presents the results for the 50 km spacecraft altitude case, Table 3 is for the 100 km case, and Table 4 is for the 200 km case. The RTN direction shown in Tables 2–4 is defined as follows. The unit direction of the radial component is always from the center of the Moon to the center of the spacecraft, and the normal component always lies in the direction of the spacecraft’s orbital angular momentum vector. The transverse component, which always lies in the velocity direction because a circular orbit is assumed in the current study, is defined to complete the orthogonal set.Table 2
Details of the expected OP errors when different lunar gravity models are applied for a spacecraft orbiting the Moon at an altitude of 50 km.
Applied gravity model
RAAN(deg)
R component error (km)
T component error (km)
N component error (km)
Total positionerror (km)
Max. in anti-R direction
Max. inR direction
Max. in anti-T direction
Max. inT direction
Max. in anti-N direction
Max. inN direction
3D Max.
3D mean
GLGM2
0
24.616
25.063
22.132
26.978
43.977
26.978
44.233
12.675
30
5.211
16.529
16.906
22.578
31.077
22.578
35.052
10.952
60
19.782
18.041
18.649
10.578
9.793
10.578
24.476
8.178
90
1.896
4.988
4.896
0.899
2.007
0.899
6.759
2.960
120
9.915
10.956
20.364
7.303
15.142
7.303
20.926
6.556
LP150Q
0
0.530
0.557
1.118
0.993
1.053
0.059
1.135
0.506
30
4.080
6.124
5.607
6.330
6.036
0.398
6.516
3.239
60
0.816
0.741
0.811
0.419
0.206
0.630
0.870
0.422
90
0.100
0.100
0.113
0.000
0.053
0.052
0.113
0.052
120
0.084
0.085
0.085
0.049
0.075
0.028
0.095
0.042
GRAIL-660B
0
0.137
0.167
0.297
0.298
0.002
0.299
0.307
0.168
30
0.201
0.125
0.196
0.174
0.023
0.189
0.212
0.090
60
0.083
0.094
0.056
0.086
0.096
0.010
0.105
0.070
90
0.245
0.244
0.001
0.251
0.136
0.108
0.252
0.140
120
0.144
0.134
0.088
0.152
0.027
0.132
0.157
0.062Table 3
Details of the expected OP errors when different lunar gravity models are applied for a spacecraft orbiting the Moon at an altitude of 100 km.
Applied gravity model
RAAN(deg)
R component error (km)
T component error (km)
N component error (km)
Total positionerror (km)
Max. in anti-R direction
Max. inR direction
Max. in anti-T direction
Max. inT direction
Max. in anti-N direction
Max. inN direction
3D Max.
3D mean
GLGM2
0
16.080
16.672
14.425
16.957
28.174
3.596
28.812
8.199
30
2.315
11.018
11.445
15.702
20.724
4.114
22.831
6.555
60
10.045
9.997
9.277
6.317
5.486
9.696
12.822
4.235
90
1.162
3.072
3.308
0.056
1.129
1.295
4.184
1.980
120
6.346
7.072
12.343
4.706
9.357
3.837
12.805
4.199
LP150Q
0
0.621
0.619
1.249
1.211
1.255
0.030
1.295
0.616
30
2.926
4.118
4.031
4.299
4.097
0.283
4.476
2.230
60
0.212
0.188
0.205
0.113
0.049
0.157
0.219
0.109
90
0.116
0.114
0.116
0.001
0.059
0.059
0.120
0.060
120
0.052
0.056
0.064
0.033
0.048
0.013
0.066
0.028
GRAIL-660B
0
0.126
0.157
0.279
0.275
0.001
0.279
0.287
0.157
30
0.183
0.121
0.185
0.169
0.022
0.169
0.197
0.084
60
0.077
0.087
0.053
0.083
0.090
0.009
0.097
0.066
90
0.227
0.229
0.020
0.232
0.127
0.100
0.236
0.131
120
0.135
0.121
0.086
0.136
0.023
0.117
0.146
0.058Table 4
Details of the expected OP errors when different lunar gravity models are applied for a spacecraft orbiting the Moon at an altitude of 200 km.
Applied gravity model
RAAN(deg)
R component error (km)
T component error (km)
N component error (km)
Total positionerror (km)
Max. in anti-R direction
Max. inR direction
Max. in anti-T direction
Max. inT direction
Max. in anti-N direction
Max. inN direction
3D Max.
3D mean
GLGM2
0
6.448
6.779
5.267
6.204
11.404
2.117
11.724
3.340
30
1.643
4.440
4.243
5.848
8.024
2.631
8.903
2.471
60
3.268
3.615
3.380
2.371
2.415
3.387
4.252
1.302
90
0.936
1.365
1.796
0.007
0.515
0.493
1.838
0.972
120
3.100
3.046
5.061
2.355
3.912
1.707
5.284
2.031
LP150Q
0
0.053
0.052
0.115
0.096
0.113
0.010
0.118
0.050
30
0.962
1.442
1.350
1.453
1.421
0.099
1.550
0.768
60
0.080
0.076
0.095
0.052
0.025
0.081
0.102
0.047
90
0.088
0.082
0.099
0.005
0.051
0.054
0.101
0.045
120
0.060
0.058
0.074
0.034
0.055
0.017
0.075
0.032
GRAIL-660B
0
0.110
0.136
0.243
0.247
0.002
0.248
0.251
0.137
30
0.165
0.101
0.160
0.147
0.020
0.151
0.171
0.073
60
0.068
0.077
0.048
0.074
0.079
0.008
0.085
0.058
90
0.202
0.203
0.013
0.204
0.115
0.089
0.207
0.115
120
0.120
0.112
0.073
0.126
0.023
0.105
0.128
0.051For the 50 km altitude case, the OP value obtained when using the GLGM2 model could be, at most, approximately 25.063 km in the radial direction, 26.978 km in the transverse direction, and approximately 43.977 km in the antinormal direction. If the LP150Q model is applied, OP errors are reduced to the order of several kilometers in each direction, and with the GRAIL660B model, errors are further reduced to within several hundreds of meters. Similar to the total 3D position error case, the magnitude of expected errors for each component is greatly reduced if the spacecraft’s altitude is increased to 100 and 200 km, as shown in Tables3 and 4, respectively. Note that the transverse direction error is closely related to the latitude of the spacecraft’s observation point, the normal direction error is related to the longitude, and the radial direction error is related to the altitude of the spacecraft, which again leads to the definition of the requirements for the relevant payload, especially the optical system, of the spacecraft during the design phase. Therefore, OP errors not only eventually lead to errors in expected observation target sites on the lunar surface but will also require more effort to calibrate the resolution of the images. Because of the improved knowledge regarding the selenopotential model obtained over the last 50 years, the OP performance has greatly improved, and an appropriate gravitational field model can be properly selected based on the requirements of the given mission. For example, suppose that a mission analyst intends to use the full degree and order of the gravitational harmonic coefficients while analyzing the spacecraft bus pointing budgets for an orbit at an altitude of 100 km. If the given tolerance of the OP errors for the mission itself is within several hundred of meters with respect to the truth model (regarded as the OP results of the GRAIL1500E model in the current work), then applying the GRAIL660B model instead of the GRAIL1500E model will be acceptable for the analysis. In addition, the results presented in this section can also enable a more in-depth analysis of scientific results.


### 4.2. Minimum Degree and Order of Harmonics Requirement

In the previous section, the expected OP errors were investigated while applying different selenopotential models with consideration of the full degree and order of their gravitational harmonics. However, in real flight operations, not only the processing time of the ground system but also the accuracy of forces modeled during OP as well as OD should be regarded. These processing times and accuracies of modeled forces are greatly affected by the degree and order of the gravitational harmonics considered during the simulation. Therefore, the current section analyzes the minimum required degree and order of harmonics that will match the OP performances simulated with consideration of the full degree and order of the gravitational harmonics. Unlike the results presented in Section4.1, which could aid a mission analyst in the selection of an appropriate selenopotential model to satisfy mission-oriented requirements, the results of the current section are expected to serve as practical guidelines for the selection of the proper degree and order of the harmonic coefficients with respect to different spacecraft orbital altitudes.To investigate the minimum degree and order of the harmonics requirements, the GRAIL1500E selenopotential model is applied, and four days of OP results of the 1500th degree and 1500th order are considered the reference OP solutions. Then, a total of five different degree and order cases, from the 100th to the 500th degree and order in steps of 100, are compared. The OP performances at different altitudes, 50, 100, and 200 km, are also compared. In Figures4–6, the expected 3D total position errors, with consideration of different degrees and orders of the harmonic coefficients, are depicted. Figure 4 represents the case in which the spacecraft’s altitude is 50 km, Figure 5 represents 100 km, and Figure 6 represents 200 km.Figure 4
Expected 3D total position OP errors when different degrees and orders of harmonic coefficients are applied for a spacecraft with an altitude of 50 km.Figure 5
Expected 3D total position OP errors when different degrees and orders of harmonic coefficients are applied for a spacecraft with an altitude of 100 km.Figure 6
Expected 3D total position OP errors when different degrees and orders of harmonic coefficients are applied for a spacecraft with an altitude of 200 km.According to Figure4, if the altitude of the spacecraft orbiting the Moon is planned to be approximately 50 km, then a minimum of a 200th degree and 200th order of the gravitational harmonic coefficients is required to keep the OP errors less than several tens (specifically, less than 20 m) of meters over the 4 days. If the OP is performed using 100th-degree and 100th-order harmonic coefficients, the expected maximum 3D position error is approximately 301.24 m. If a degree and order of 300 are considered, a 3D OP error of less than 0.5 m is expected with respect to the results derived by considering the full degree and order of the harmonic coefficients. If the spacecraft’s altitude is increased to 100 km, a degree and order of 100 are determined to be sufficient to keep the 3D positions of OP errors less than several meters, that is, 6.62 m for this case, while they must be less than several centimeters if the considered degree and order are increased to 200, as shown in Figure 5. If the spacecraft’s altitude is increased to 200 km, as expected, the consideration of a degree and order of less than 100 is sufficient to keep the OP accuracy within several centimeters. Regarding the results discussed above, considering the full degree and order of the harmonic coefficients may not always be the best solution; rather, appropriate ranges of the degree and order of the harmonic coefficients may be determined during the early mission design phases with consideration of the mission success criteria, expected orbit altitude, instrument performances, expected knowledge regarding the spacecraft attitude, mission-oriented OP, OPD and OD requirements, permitted ground processing time, and so on.Unlike the OP errors for different lunar gravity models shown in Section4.1, the errors shown in the current section are much smaller. This difference occurs because the current simulation is performed using the same selenopotential model but different degrees and orders of the gravitational harmonic coefficients, as setting the degree and order of the harmonic coefficients to lower values may not degrade the modeling error of the selected gravitational field itself. Indeed, the number of expected OP errors for a particular selenopotential model (as shown in Section 4.1) will likely be maintained regardless of the degree and order of the gravitational harmonic coefficients selected under the condition that the OP errors computed based on the selected degree and order and the full potential model are negligible. Actually, for the case of the GRAIL660B model, the minimum degree and order of harmonics required to match the OP performances simulated using the truth model (which considers the full degree and order) displayed a behavior very similar to that of the results obtained using the GRAIL1500E model. With the GRAIL660B model, a minimum of a 200th degree and 200th order of the gravitational harmonics for the 50 km case and a degree and order less than the 100th degree and 100th order for the 100 km and 200 km cases are required to restrict the OP errors to within several tens of meters and less than several meters, respectively. Then, the expected OP errors obtained using the GRAIL1500E and GRAIL660B models are again compared. For this simulation, the gravitational harmonics considered for both models are as follows: 200th degree and 200th order for 50 km altitude case and 100th degree and 100th order for 100 and 200 km altitude cases. As expected, the OP errors observed from this simulation almost match the results shown in Figure 3 in the previous section. Once again, note that the current simulation results correspond only to the OP errors, which may offer useful insights when establishing spacecraft budgets or scientific results analysis, both of which are closely related to dynamic modeling errors. Regarding the OPD performances, obtained using a combination of different OD knowledge bases, they are influenced more by the selenopotential model chosen than by the particular degrees and orders of harmonic coefficients applied as long as the lowest degrees and orders of the harmonic coefficients are sufficiently considered for the relevant analysis. Given that the current lunar missions usually require stricter system requirements, the tendency of expected OP errors discovered in the current study should be considered during the early system design phase.


## 4.1. Orbit Propagation Performance Comparisons

In this section, the OP performances are compared, and the expected OP errors due to selecting different selenopotential models are investigated. To compare the OP performances, the full degree and order of the harmonic coefficients are considered during the simulations of each selected lunar gravity model. A 70th-degree and 70th-order are considered for GLGM2, 150th for LP150Q, 660th for GRAIL660B, and 1500th for the GRAIL1500E model. By regarding the simulated OP results of the GRAIL1500E model as references, Figures1–3 compare the 3D position differences when the same OP is used for different lunar gravity models. Figure 1 shows the comparison between GRAIL1500E and GLGM2, Figure 2 shows that between GRAIL1500E and LP150Q, and Figure 3 shows that between GRAIL1500E and GRAIL660B. Additionally, in each of the figures, subfigure (a) presents the comparison results for the 50 km spacecraft orbital altitude case, (b) presents the 100 km case, and (c) presents the 200 km case.Figure 1
OP performance comparison between the GRAIL1500E and GLGM2 models. Subfigure (a) presents the 50 km spacecraft altitude case, (b) presents the 100 km case, and (c) presents the 200 km case.
(a)
(b)
(c)Figure 2
OP performance comparison between the GRAIL1500E and LP150Q models. Subfigure (a) presents the 50 km spacecraft altitude case, (b) presents the 100 km case, and (c) presents the 200 km case.
(a)
(b)
(c)Figure 3
OP performance comparison between the GRAIL1500E and GRAIL660B models. Subfigure (a) presents the 50 km spacecraft altitude case, (b) presents the 100 km case, and (c) presents the 200 km case.
(a)
(b)
(c)Figures1–3 show that a significant number of OP errors are induced simply by considering different selenopotential models. Additionally, as expected, when the spacecraft altitude decreases, the number of associated OP errors tends to increase. For example, the comparison of the 3D positions of OP errors between the GRAIL1500E and GLGM2 models shows a maximum of approximately 44.233 km errors when the spacecraft’s altitude is 50 km, approximately 28.812 km errors when the altitude is 100 km, and approximately 11.724 km errors when the altitude is 200 km. In addition, the number of expected OP errors has different trends based on the given initial orbital element of RAAN. For the GLGM2 comparison case, the maximum number of OP errors was observed for a RAAN of 0 deg, as was already discussed.However, the minimum number of OP errors was observed for an initial RAAN of 90 deg, and 3D position errors of approximately 6.759 km were found. This is due to the uncertainties of the lunar gravity model present during the early phase of the OP time span. Indeed, if a RAAN of 0 deg is given, the initial location of the spacecraft is determined to be at the far side of the Moon; by contrast, a RAAN of 90 deg indicates an initial location at the near side of the Moon. These results indirectly show the modeling accuracy of lunar far-side gravity fields of the GLGM2 model. Nevertheless, the OP errors caused solely by applying different lunar gravity field models are found to be rather significant, even though GLGM2 is a very old lunar gravity field model (having been released in the late 1990s). For the comparison case with the LP150Q model, which is shown in Figure2, OP errors are still found within orders of several kilometers. The expected maximum 3D-position OP error for the four days of propagation is found to be approximately 6.516 km at an altitude of 50 km, 4.476 km at an altitude of 100 km, and 1.550 km at an altitude of 200 km. Unlike the comparison case with the GLGM2 model, the maximum 3D-position OP errors are found for a RAAN of 30 deg for the LP150Q model, with the spacecraft still initially located on the far side of the Moon. However, when the RAAN is 90 and 120 deg, indicating an initial location on the near side of the Moon, the 3D OP errors of the LP150Q model at an altitude of 50 km are approximately 113 m and 95 m, respectively. The OP comparison results for the case with the GRAIL660B model are presented in Figure 3. Unlike the other two cases (GLGM2 and LP150Q), the OP errors simulated using the GRAIL660B model remained within several hundreds of meters with very similar trends regardless of the altitude. Very minor differences due to the altitude are actually observed, but they are not as significant as those observed in the former two cases. Additionally, the initially given RAANs did not seriously affect the OP performances. This result may be due to the model resolution improvement (surface resolution ranging from 9.3 to 14.6 km for the GRAIL660B model [47] and of 3.6 km for GRAIL1500E [53]) rather than the far-side modeling accuracy improvement, as occurred between LP150Q/GLGM-2 and the GRAIL660B model.Because Figures1–3 only plot the total 3D position differences, details of the results, including the position components of the Radial-Transverse-Normal (RTN) direction, are provided separately in Tables 2–4. Table 2 presents the results for the 50 km spacecraft altitude case, Table 3 is for the 100 km case, and Table 4 is for the 200 km case. The RTN direction shown in Tables 2–4 is defined as follows. The unit direction of the radial component is always from the center of the Moon to the center of the spacecraft, and the normal component always lies in the direction of the spacecraft’s orbital angular momentum vector. The transverse component, which always lies in the velocity direction because a circular orbit is assumed in the current study, is defined to complete the orthogonal set.Table 2
Details of the expected OP errors when different lunar gravity models are applied for a spacecraft orbiting the Moon at an altitude of 50 km.
Applied gravity model
RAAN(deg)
R component error (km)
T component error (km)
N component error (km)
Total positionerror (km)
Max. in anti-R direction
Max. inR direction
Max. in anti-T direction
Max. inT direction
Max. in anti-N direction
Max. inN direction
3D Max.
3D mean
GLGM2
0
24.616
25.063
22.132
26.978
43.977
26.978
44.233
12.675
30
5.211
16.529
16.906
22.578
31.077
22.578
35.052
10.952
60
19.782
18.041
18.649
10.578
9.793
10.578
24.476
8.178
90
1.896
4.988
4.896
0.899
2.007
0.899
6.759
2.960
120
9.915
10.956
20.364
7.303
15.142
7.303
20.926
6.556
LP150Q
0
0.530
0.557
1.118
0.993
1.053
0.059
1.135
0.506
30
4.080
6.124
5.607
6.330
6.036
0.398
6.516
3.239
60
0.816
0.741
0.811
0.419
0.206
0.630
0.870
0.422
90
0.100
0.100
0.113
0.000
0.053
0.052
0.113
0.052
120
0.084
0.085
0.085
0.049
0.075
0.028
0.095
0.042
GRAIL-660B
0
0.137
0.167
0.297
0.298
0.002
0.299
0.307
0.168
30
0.201
0.125
0.196
0.174
0.023
0.189
0.212
0.090
60
0.083
0.094
0.056
0.086
0.096
0.010
0.105
0.070
90
0.245
0.244
0.001
0.251
0.136
0.108
0.252
0.140
120
0.144
0.134
0.088
0.152
0.027
0.132
0.157
0.062Table 3
Details of the expected OP errors when different lunar gravity models are applied for a spacecraft orbiting the Moon at an altitude of 100 km.
Applied gravity model
RAAN(deg)
R component error (km)
T component error (km)
N component error (km)
Total positionerror (km)
Max. in anti-R direction
Max. inR direction
Max. in anti-T direction
Max. inT direction
Max. in anti-N direction
Max. inN direction
3D Max.
3D mean
GLGM2
0
16.080
16.672
14.425
16.957
28.174
3.596
28.812
8.199
30
2.315
11.018
11.445
15.702
20.724
4.114
22.831
6.555
60
10.045
9.997
9.277
6.317
5.486
9.696
12.822
4.235
90
1.162
3.072
3.308
0.056
1.129
1.295
4.184
1.980
120
6.346
7.072
12.343
4.706
9.357
3.837
12.805
4.199
LP150Q
0
0.621
0.619
1.249
1.211
1.255
0.030
1.295
0.616
30
2.926
4.118
4.031
4.299
4.097
0.283
4.476
2.230
60
0.212
0.188
0.205
0.113
0.049
0.157
0.219
0.109
90
0.116
0.114
0.116
0.001
0.059
0.059
0.120
0.060
120
0.052
0.056
0.064
0.033
0.048
0.013
0.066
0.028
GRAIL-660B
0
0.126
0.157
0.279
0.275
0.001
0.279
0.287
0.157
30
0.183
0.121
0.185
0.169
0.022
0.169
0.197
0.084
60
0.077
0.087
0.053
0.083
0.090
0.009
0.097
0.066
90
0.227
0.229
0.020
0.232
0.127
0.100
0.236
0.131
120
0.135
0.121
0.086
0.136
0.023
0.117
0.146
0.058Table 4
Details of the expected OP errors when different lunar gravity models are applied for a spacecraft orbiting the Moon at an altitude of 200 km.
Applied gravity model
RAAN(deg)
R component error (km)
T component error (km)
N component error (km)
Total positionerror (km)
Max. in anti-R direction
Max. inR direction
Max. in anti-T direction
Max. inT direction
Max. in anti-N direction
Max. inN direction
3D Max.
3D mean
GLGM2
0
6.448
6.779
5.267
6.204
11.404
2.117
11.724
3.340
30
1.643
4.440
4.243
5.848
8.024
2.631
8.903
2.471
60
3.268
3.615
3.380
2.371
2.415
3.387
4.252
1.302
90
0.936
1.365
1.796
0.007
0.515
0.493
1.838
0.972
120
3.100
3.046
5.061
2.355
3.912
1.707
5.284
2.031
LP150Q
0
0.053
0.052
0.115
0.096
0.113
0.010
0.118
0.050
30
0.962
1.442
1.350
1.453
1.421
0.099
1.550
0.768
60
0.080
0.076
0.095
0.052
0.025
0.081
0.102
0.047
90
0.088
0.082
0.099
0.005
0.051
0.054
0.101
0.045
120
0.060
0.058
0.074
0.034
0.055
0.017
0.075
0.032
GRAIL-660B
0
0.110
0.136
0.243
0.247
0.002
0.248
0.251
0.137
30
0.165
0.101
0.160
0.147
0.020
0.151
0.171
0.073
60
0.068
0.077
0.048
0.074
0.079
0.008
0.085
0.058
90
0.202
0.203
0.013
0.204
0.115
0.089
0.207
0.115
120
0.120
0.112
0.073
0.126
0.023
0.105
0.128
0.051For the 50 km altitude case, the OP value obtained when using the GLGM2 model could be, at most, approximately 25.063 km in the radial direction, 26.978 km in the transverse direction, and approximately 43.977 km in the antinormal direction. If the LP150Q model is applied, OP errors are reduced to the order of several kilometers in each direction, and with the GRAIL660B model, errors are further reduced to within several hundreds of meters. Similar to the total 3D position error case, the magnitude of expected errors for each component is greatly reduced if the spacecraft’s altitude is increased to 100 and 200 km, as shown in Tables3 and 4, respectively. Note that the transverse direction error is closely related to the latitude of the spacecraft’s observation point, the normal direction error is related to the longitude, and the radial direction error is related to the altitude of the spacecraft, which again leads to the definition of the requirements for the relevant payload, especially the optical system, of the spacecraft during the design phase. Therefore, OP errors not only eventually lead to errors in expected observation target sites on the lunar surface but will also require more effort to calibrate the resolution of the images. Because of the improved knowledge regarding the selenopotential model obtained over the last 50 years, the OP performance has greatly improved, and an appropriate gravitational field model can be properly selected based on the requirements of the given mission. For example, suppose that a mission analyst intends to use the full degree and order of the gravitational harmonic coefficients while analyzing the spacecraft bus pointing budgets for an orbit at an altitude of 100 km. If the given tolerance of the OP errors for the mission itself is within several hundred of meters with respect to the truth model (regarded as the OP results of the GRAIL1500E model in the current work), then applying the GRAIL660B model instead of the GRAIL1500E model will be acceptable for the analysis. In addition, the results presented in this section can also enable a more in-depth analysis of scientific results.


## 4.2. Minimum Degree and Order of Harmonics Requirement

In the previous section, the expected OP errors were investigated while applying different selenopotential models with consideration of the full degree and order of their gravitational harmonics. However, in real flight operations, not only the processing time of the ground system but also the accuracy of forces modeled during OP as well as OD should be regarded. These processing times and accuracies of modeled forces are greatly affected by the degree and order of the gravitational harmonics considered during the simulation. Therefore, the current section analyzes the minimum required degree and order of harmonics that will match the OP performances simulated with consideration of the full degree and order of the gravitational harmonics. Unlike the results presented in Section4.1, which could aid a mission analyst in the selection of an appropriate selenopotential model to satisfy mission-oriented requirements, the results of the current section are expected to serve as practical guidelines for the selection of the proper degree and order of the harmonic coefficients with respect to different spacecraft orbital altitudes.To investigate the minimum degree and order of the harmonics requirements, the GRAIL1500E selenopotential model is applied, and four days of OP results of the 1500th degree and 1500th order are considered the reference OP solutions. Then, a total of five different degree and order cases, from the 100th to the 500th degree and order in steps of 100, are compared. The OP performances at different altitudes, 50, 100, and 200 km, are also compared. In Figures4–6, the expected 3D total position errors, with consideration of different degrees and orders of the harmonic coefficients, are depicted. Figure 4 represents the case in which the spacecraft’s altitude is 50 km, Figure 5 represents 100 km, and Figure 6 represents 200 km.Figure 4
Expected 3D total position OP errors when different degrees and orders of harmonic coefficients are applied for a spacecraft with an altitude of 50 km.Figure 5
Expected 3D total position OP errors when different degrees and orders of harmonic coefficients are applied for a spacecraft with an altitude of 100 km.Figure 6
Expected 3D total position OP errors when different degrees and orders of harmonic coefficients are applied for a spacecraft with an altitude of 200 km.According to Figure4, if the altitude of the spacecraft orbiting the Moon is planned to be approximately 50 km, then a minimum of a 200th degree and 200th order of the gravitational harmonic coefficients is required to keep the OP errors less than several tens (specifically, less than 20 m) of meters over the 4 days. If the OP is performed using 100th-degree and 100th-order harmonic coefficients, the expected maximum 3D position error is approximately 301.24 m. If a degree and order of 300 are considered, a 3D OP error of less than 0.5 m is expected with respect to the results derived by considering the full degree and order of the harmonic coefficients. If the spacecraft’s altitude is increased to 100 km, a degree and order of 100 are determined to be sufficient to keep the 3D positions of OP errors less than several meters, that is, 6.62 m for this case, while they must be less than several centimeters if the considered degree and order are increased to 200, as shown in Figure 5. If the spacecraft’s altitude is increased to 200 km, as expected, the consideration of a degree and order of less than 100 is sufficient to keep the OP accuracy within several centimeters. Regarding the results discussed above, considering the full degree and order of the harmonic coefficients may not always be the best solution; rather, appropriate ranges of the degree and order of the harmonic coefficients may be determined during the early mission design phases with consideration of the mission success criteria, expected orbit altitude, instrument performances, expected knowledge regarding the spacecraft attitude, mission-oriented OP, OPD and OD requirements, permitted ground processing time, and so on.Unlike the OP errors for different lunar gravity models shown in Section4.1, the errors shown in the current section are much smaller. This difference occurs because the current simulation is performed using the same selenopotential model but different degrees and orders of the gravitational harmonic coefficients, as setting the degree and order of the harmonic coefficients to lower values may not degrade the modeling error of the selected gravitational field itself. Indeed, the number of expected OP errors for a particular selenopotential model (as shown in Section 4.1) will likely be maintained regardless of the degree and order of the gravitational harmonic coefficients selected under the condition that the OP errors computed based on the selected degree and order and the full potential model are negligible. Actually, for the case of the GRAIL660B model, the minimum degree and order of harmonics required to match the OP performances simulated using the truth model (which considers the full degree and order) displayed a behavior very similar to that of the results obtained using the GRAIL1500E model. With the GRAIL660B model, a minimum of a 200th degree and 200th order of the gravitational harmonics for the 50 km case and a degree and order less than the 100th degree and 100th order for the 100 km and 200 km cases are required to restrict the OP errors to within several tens of meters and less than several meters, respectively. Then, the expected OP errors obtained using the GRAIL1500E and GRAIL660B models are again compared. For this simulation, the gravitational harmonics considered for both models are as follows: 200th degree and 200th order for 50 km altitude case and 100th degree and 100th order for 100 and 200 km altitude cases. As expected, the OP errors observed from this simulation almost match the results shown in Figure 3 in the previous section. Once again, note that the current simulation results correspond only to the OP errors, which may offer useful insights when establishing spacecraft budgets or scientific results analysis, both of which are closely related to dynamic modeling errors. Regarding the OPD performances, obtained using a combination of different OD knowledge bases, they are influenced more by the selenopotential model chosen than by the particular degrees and orders of harmonic coefficients applied as long as the lowest degrees and orders of the harmonic coefficients are sufficiently considered for the relevant analysis. Given that the current lunar missions usually require stricter system requirements, the tendency of expected OP errors discovered in the current study should be considered during the early system design phase.


## 5. Conclusions

In this work, the effect of applying different lunar gravity field models on a propagating lunar orbiter’s state is analyzed. Obtaining insights into the expected OP errors resulting from the selection of a particular lunar gravity field model is typically important when establishing spacecraft budgets, which are closely related to dynamic modeling errors. The history of selenopotential model evolution is briefly presented, and four representative lunar gravity field models (GLGM2, LP150Q, GRAIL660B, and GRAIL1500E) are selected for the analysis. Three different fictitious lunar circular polar orbits with altitudes of 50 km, 100 km, and 200 km are generated, and the expected OP errors are investigated during four days of propagation. As a result, a significant number of OP errors are discovered to have occurred simply because of the use of different selenopotential models. Compared to the solutions obtained using the GRAIL1500E model, a very recent lunar gravity field model, the overall 3D propagation errors could reach up to several tens of kilometers if the GLGM2 model is applied for a 50 km altitude and up to several hundreds of meters for the 50, 100, and 200 km orbiting cases even when the GRAIL660B model is applied. These propagation errors are also found to be affected by the spacecraft’s initial conditions with respect to the Moon’s geometric location, that is, the far or near side of the Moon, as they reflect the modeling accuracies of the applied lunar gravity field. With the GRAIL1500E model as a basis, appropriate ranges of the degrees and orders of harmonics, which must be chosen to match the OP performances simulated using the full degrees and orders of gravitational harmonics, are also investigated. To offer the best OP performances, a degree and order of at least 200 should be considered to maintain errors of less than several tens of meters for the 50 km altitude case. If the spacecraft flying altitude is increased to 100 km, a degree and order of at least 100 should be considered to maintain similar propagation accuracy. For the 200 km case, a degree and order of less than 100 are sufficient to yield centimeter-level accuracy. The results of the current work are expected to serve as practical guidelines for the field of system budget analysis, mission design, mission operations, and the analysis of scientific results and will eventually aid in the further establishment of an OPD and OD analysis strategy.



---

*Source: 5493679-2017-08-10.xml*