Fitting a triaxial ellipsoid to a set of quasi-selenoidal points

2.1 Description of the lunar gravity field models

The lunar gravity field spherical harmonic models used vary in degree, with GLGM-1 [15] being the earliest model available and of the lowest degree, GL0660B [16] and GRGM660PRIM [17] being the models of the highest degree, and AIUB-GRL350A and AIUB-GRL350B [18] being the latest models available for use.

Both models GLGM-1 and GLGM-2 were published by Lemoine et al. in 1994 and 1995, respectively, and are both of 70° and order. GLGM-1 uses S-band Doppler tracking data from the Clementine mission and historical data from the Lunar Orbiter Satellites 1 to 5 and the Apollo 15 subsatellite, while GLGM-2 also includes historical data from the Apollo 16 subsatellite. The Lunar Orbiter Satellites were launched from 1966 through 1967 and were placed in elliptical orbits about the Moon. The Apollo 15 subsatellite was launched in July 1971 and the Apollo 16 in April 1972, both of which were placed in near-circular orbits and were the first to achieve the acquirement of an extensive amount of data from a lunar satellite in low-altitude near-circular orbit. The spatial coverage of the tracking data acquired from both aforementioned missions was incomplete and neither provided direct tracking data over large portions of the lunar farside due to the Moon’s synchronous rotation. The Clementine spacecraft was launched on January 25, 1994, and was placed into an elliptical orbit about the Moon. The tracking data provided by the Clementine mission allowed the improvement of the low-degree coefficients of the gravity field. Moreover, the GLGM-2 model improved its predecessor, GLGM-1, not only by the inclusion of the Apollo 16 subsatellite’s data, but by also applying a different weight to the data and calibrating the solution, and managed to resolve the gravitational signatures of both farside and nearside basins and mascons, even though direct tracking was still not available on a large part of the lunar farside [19].

The JGL models, namely JGL075D1, JGL075G1, JGL100J1, JGL100K1, JGL150Q1 and JGL165P1, were published by [20] and their degree varies from 75 to 165, as listed in Table 1. These gravity models use two-way coherent Doppler tracking data [20] from the Lunar Prospector mission launched on January 7, 1998, by the NASA Discovery Program, along with historical Doppler data from previous missions, such as Clementine. The Lunar Prospector was placed in the polar orbit of the Moon and at a low altitude. As a result, it provided high-resolution measurements of the crustal gravity field across the lunar near side, long-wavelength gravity anomalies over the entire Moon, and further constraints on the lunar moment of inertia [20], thus improving significantly the measurements provided by previous missions, such as the Apollo program or Clementine, both on the lunar near side and the lunar farside, despite still not providing direct tracking data on the latter.

The Gravity Recovery and Interior Laboratory (GRAIL) mission – the first dedicated gravity mission in planetary science [21] – whose primary phase was launched on September 10, 2011, was the first to acquire direct tracking data of high accuracy of the lunar farside by using two identical spacecraft, Ebb or GRAIL-A and Flow or GRAIL-B, placed in the near-polar circular orbits [16, 17]. By using DSN two-way S-band Doppler, DSN one-way X-band Doppler, and instantaneous interspacecraft Ka-band range-rate (KBRR) data [16, 17], the Jet Propulsion Laboratory (JPL) group developed a spherical harmonic lunar gravity field model of degree 660, GL0660B [16], which significantly improved the lunar gravity field and provided new information concerning the Moon’s interior. The Goddard Space Flight Center (GSFC) derived the GRGM660PRIM [17] model using the same type of data provided by the GRAIL primary mission. The GRGM660PRIM is also of degree 660, and greatly improved the global root-mean-square (RMS) error in the gravity anomalies compared to pre-GRAIL models. In addition, the free-air gravity anomalies are presented with a dramatic increase, especially over the lunar farside [17]. These two spherical harmonic models are counterparts of each other and are important both because of their high resolution, and because they are the first to contain highly accurate information about the lunar farside.

The AIUB-GRL200A and AIUB-GRL200B models, both of maximum degree and order 200, were derived by Arnold et al. and are also based on the data of GRAIL’s primary mission phase [22]. For these models, the Celestial Mechanics Approach (CMA) [23] was used via a development version of the Bernese GNSS Software [22]. In more detail, the gravity field recovery was treated as an extended orbit determination problem [22]; Ka-band range data series from GRAIL’s primary mission was used as observations, while GNI1B satellite positions of the GRAIL probes – as provided by NASA JPL – were used as pseudo-observations. For the AIUB-GRL200A, the latest lunar gravity field model GRGM900C [24] was used as the a priori field up to degree and order 200, while for the AIUB-GRL200B, the same model was used up to degree and order 660. However, independence from the a priori field is established by presenting a solution that converges on the AIUB-GRL200A with a pre-GRAIL lunar gravity field model, and specifically JGL165P1 [25], as the a priori field. The CMA, along with the availability of the GNI1B positions, create the possibility of recovering the lunar gravity field with sufficient accuracy without prior processing of DSN Doppler data, although implementing them into the Bernese GNSS software would allow the development of a fully stand-alone solution [22].

The models GrazLGM300c, and GrazLGM420a, GrazLGM420b, GrazLGM420b+ are up to degree and order 300 and 420, respectively, with degree 420 being equal to the average global resolution of GRAIL primary mission data [16], and were derived within the project GRAZIL. The latter aims for a highly accurate lunar gravity field recovery by using inter-satellite Ka-band ranging (KBR) measurements of GRAIL’s primary mission [26]. The first step to doing so is the precise orbit determination of the twin GRAIL satellites. The first two Graz models relied on the GNI1B positions provided by NASA JPL, while the two latter utilized DSN Doppler S-band radiometric tracking data to determine the absolute positions of the two satellites, thus making the GrazLGM420b the first completely independent lunar gravity field model [21]. For the lunar gravity field recovery, an integral equation approach using short orbital arcs is applied [26]. For both parts, the Graz team uses internally developed software packages (ORCA, GROOPS) [21].

Lastly, the AIUB-GRL350A and AIUB-GRL350B models were developed by [18] and are both of degree and order 350. For these models, data from the GRAIL primary mission were used, and orbit determination was conducted in the Bernese GNSS software by using DSN doppler tracking data, and a combination of Doppler and Ka-band range-rate (KBRR) data [18]. It is important to note that for both orbit determination and lunar gravity field recovery, different parametrizations were used during the derivation of the above models, thus providing different results. Also, GRAIL and pre-GRAIL (Selene) lunar gravity field models were used as a priori fields (GRGM900C, [24], and SGM150J, [27]).

2.2 Data processing and statistics

The data sets of points in space, corresponding to points of the quasi-selenoid, at particular values of geodetic coordinates φ and λ, using values of the height anomaly ζ, a.k.a. the height of the quasi-selenoid provided by the ICGEM service for the lunar gravity field models [9], in the following reference system (MOON):

1. Reference radius (Radius): R = 1,738,140 m

2. Flattening (Flat): 3240

3. Product of gravitational constant and mass of the Earth: Gm = 4.90280107 ⋅ 10¹² m³ s⁻²

4. Angular rotational speed (Omega): ω = 2.661621 ⋅ 10⁻⁶ rad s⁻¹

5. Grid: φ = − 90 ° , 90 ° , λ = − 180 ° , 180 ° , while the Grid Step was determined by the maximum degree (n) of the gravity field model as 180°/n. For example, the model GLGM-1 has n = 70 and hence resolution 180°/n ≈ 2.5°.

6. Height over ellipsoid: h = 0 m

Also, Gaussian filtering was not applied and the zero-degree term was taken into consideration, maintaining the model’s defined tidal system. Subsequently, statistical analysis for each data set was executed and the results are presented in Table 1. Specifically, we compute the statistical indexes mean ζ , min ζ , max ζ and the Root Mean Square (RMS) by:

where k is the number of points. From the results of Table 1, we observe that the different quasi-selenoid models have a very small impact on the statistics of the data sets. This is due to the common reference surface which is a biaxial ellipsoid (oblate spheroid). Finally, the geodetic coordinates φ i , λ i , ζ i of a point on the quasi-selenoid are related to the corresponding Cartesian coordinates x i , y i , z i by well-known expressions, where major-semiaxis a = R, flattening f = 1/3240 and ellipsoidal height h _i ≡ ζ _i .

In all numerical computations that are presented in this work, we used a personal computer running a 64 bit Windows operating system. The main characteristics of the hardware were: Intel Core i3-1011OU CPU (clocked at 2.1 GHz) and 8 GB of RAM. All algorithms were coded in MATLAB, which provides a precision of 16 digits.

4.1 General case

Using the 18 data sets described in Section 2, we perform the adjustment computations and obtain the results for the nine ellipsoid parameters and their standard deviations (σ – square roots of diagonal elements of the v–c matrix of the derived ellipsoidal parameters), in the general case (G), as presented in Figures 1 –3.

Figure 1:

The translation parameters t _x , t _y and t _z for all models along with a vertical error bar at each result.

Figure 2:

The rotation angles θ _x , θ _y and θ _z for all models along with a vertical error bar at each result.

Figure 3:

The ellipsoid semiaxes a _x , a _y and a _z for all models along with a vertical error bar at each result.

More specifically, Figure 1 shows the values of the translation parameters t _x , t _y and t _z , and their standard deviations in the form of an error bar for each model. From the results, it is apparent the improvement of precision of the estimated parameters in the cases of the models which characterized by higher degree. However, the values of translation parameters fluctuate around zero meters and are in agreement considering their statistical error bounds.

Figure 2 presents the values of the rotation angles θ _x , θ _y and θ _z and an error bar for each one. The standard deviations of the rotation angles do not follow the same pattern as those of the translation parameters, as they don’t necessarily decrease as the model’s degree increases. We attribute this change to the numerical precision attained in the implementation of the respective complicated formulae involved in the application of the law of propagation of variances, see [11]. However, the values of parameters fluctuate around zero degrees and are in agreement considering their statistical error bounds. Therefore, the translation parameters, as well as the rotation angles, have very small impact on the statistics of the fitting.

Figure 3 exhibits the values of the three semiaxes a _x , a _y and a _z , along with their standard deviations, for each model. From the results, we note that the value of each of the semiaxes stabilizes over time and the standard deviations decrease as the model’s degree increases, while the difference between the subsequent results tends to be within their error bound. Obviously, the improvement of the lunar gravity field models resulted in more uniform values and smaller uncertainties of the triaxiality parameters.

Lastly, Figure 4 presents the values of the a posteriori standard deviation of unit weight for each model. From the results, we observe the goodness of fit of each model to data due to the common reference surface which is a triaxial ellipsoid. As a comparison, for the model GL0660B, the respective values for a biaxial ellipsoid and sphere are 1.01 ⋅ 10⁻⁴ and 1.78 ⋅ 10⁻⁴, respectively.

Figure 4:

The a posteriori standard deviation of unit weight for all models.

In Table 2, we present the resulting numerical values of the nine ellipsoidal parameters and statistics of fitting for the models GLGM-1, GL0660B, GrazLGM420b+, and ULCN, in the general case (G).

We choose to present the earliest available model GLGM-1 to compare the results with those obtained through the models GL0660B, a model of the highest degree, and GrazLGM420b+, one of the recent models with high precision, as indicated from the resulting values of standard deviations of the ellipsoidal parameters, see Figures 1–3.

From the estimates in Table 2, in the case of ULCN, we note their consistency and the very small difference in the results of the corresponding values of semiaxes and translation parameters presented in [6] (see Table 4). Also, the difference in the values of the rotation angles is due to the different definitions of angles between the two works. On the other hand, we also note the uniform values of the parameters and the improvement of the deviations in the case of the quasi-selenoidal points of the lunar gravity field models.

To facilitate a more exhaustive analysis between the two reference surfaces (quasi-selenoid and lunar topography) represented by their triaxial ellipsoids of fitting, we use the 9 parameters estimated for the model ULCN and compute the RMS value for the height anomaly of the points of the model GL0660B. The resulting value of RMS = 1745 m indicates the different positions in space of the two ellipsoids and the worse fit of the one ellipsoid on the other surface.

4.2 Special and degenerate cases

From the results presented in Figures 1 and 2, and Table 2, it is clear that the values of the translation parameters and the rotation angles are close to zero meters and degrees, respectively, and smaller than their standard deviations in all cases. Thus, we can use the special and degenerate cases discussed in Section 3 for the fitting moving forward. Moreover, we use only the GL0660B and GrazLGM420b+ models for the fitting in the special case (T) and in the two degenerate cases (B and S), the results of which are presented in Table 3.

From the resulting values of the parameters and statistics in the special case (T), we justify that the effect of the translation parameters and rotation angles is insignificant. As an additional experiment to examine the role of rotational angles, for the points of the model GL0660B and their three semiaxes, we use the three rotational angles from the model ULCN, i.e. θ _x = 7.48°, θ _y = 27.69° and θ _z = 23.96° and compute the RMS value of height anomaly. The resulting value of RMS = 175 m indicates the high impact of the orientation parameters when they are different from the adjustment values. Furthermore, the resulting values of the parameters in the degenerate case (B), are in agreement with those of the reference system “MOON” used by the ICGEM platform (i.e. a = 1,738,140 m, b = 1,737,603.5 m).

4.3 Results from other works

Table 4 presents the estimates of the ellipsoidal parameters reported in other works. In [3] the parameters are related to the equipotential surface (selenoid) while the set of values reported by [5], [6], [7], [8] represent the best fitting ellipsoids to the surface of the lunar topography using different sources of data and models, as described in the Introduction.

The resulting values of the parameters refer to different surfaces (selenoid and lunar topography) from a quasi-selenoid, therefore are not comparable with the values of quasi-selenoids in Table 2. However, we observe that the parameters of semiaxes of the selenoid credited to Burša (1971) [3] have the same differences as that of quasi-selenoid from this work, in order of 200 and 400 m, respectively. The difference in absolute values is due to the scale factor for semiaxes being computed based on absolute gravity measurement made by the first lunar-landing mission Apollo 11 at the landing site under the assumption of sufficient accuracy [3]. Furthermore, as studied by [28] the difference between the quasi-selenoid and the selenoid is in the order of ±3 m. Therefore, the resulting parameters, in this work, may be assumed quite similar to that of a selenoid.