Lunar gravity field modeling critical analysis and challenges

Available online at www.sciencedirect.com

Advances in Space Research 45 (2010) 322–349 www.elsevier.com/locate/asr

Lunar gravity field modeling critical analysis and challenges Manoranjan Sinha a,*, N.S. Gopinath b, N.K. Malik c a

Department of Aerospace Engineering, IIT Kharagpur, India b Flight Dynamics Division, ISAC, ISRO, Bangalore, India c Mission Planning and Controls, ISAC, ISRO, Bangalore, India Received 19 January 2008; received in revised form 7 August 2009; accepted 8 October 2009

Abstract This paper summarizes and provides a critical analysis of the historical developments of lunar gravitational models from the earliest use of ground based tracking systems of the Lunar Orbiter to the Lunar Prospector mission. This encompasses a comprehensive and critical analysis of the various methods used in the estimation of the gravity coefficients and the processing of large batches of diverse measurements and data types. It has been shown that weakness exists in the current models of the lunar gravity field, which is primarily due to the lack of far side lunar tracking data information, which makes the lunar potential modeling difficult but expected to be overcome as data from SELENE satellite-to-satellite tracking becomes available. Comparisons of various lunar models reveal an agreement in the low order coefficients of the spherical harmonics. However, substantial differences in the models exist in the higher-order harmonics. A numerical comparison has been presented showing the performance of all the contemporary lunar gravitational models used within the astrodynamics community and available in public domain. Improvements to the current models are part of a continuing process and the recent model improvements and future possibilities in lunar gravity modeling are discussed. A brief review of the recent missions has been presented. It is hoped that this critical review will benefit the researchers by presenting the historical as well as state of the art in this field. Ó 2009 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: Lunar gravity; Spherical harmonic coefficients

1. Introduction 1.1. A brief description of the problem The accurate knowledge of the lunar gravity model is necessary for the long life of the lunar satellite as maneuver can be optimized for the altitude maintenance which can reduce the fuel spent on-board. Secondly, the mapping of the surface using a camera requires a constant altitude in successive passes which can be achieved efficiently if the gravity model is known accurately. Moreover, study of the underlying density distribution and its selenophysical *

Corresponding author. Tel.: +91 3222283024. E-mail addresses: [email protected], [email protected], [email protected], [email protected] (M. Sinha), nsgopi@ isac.gov.in (N.S. Gopinath).

implication requires a local gravity model which can be obtained on a global basis if the global gravity model is accurate enough. But there are few impediments in achieving the above objectives. These are (a) phase lock between Moon’s spin and orbital motion resulting in lack of observational data of satellite orbit on the far side of the Moon, (b) presence of mascons beneath the center of all ringed maria, (c) accuracy of observational data limits the features to be extracted, and (d) precision of the gravity model extracted depends on the Legendre’s coefficients in the gravity model equation which increases as the degree and order of the model is increased resulting in large amount of computation required which in turn limits the degree and order to which the gravity model can be considered. The remedy to the first problem is Satellite-to-Satellite Tracking (SST) (Floberghagen et al., 1996) which will be discussed later. The second problem has been tackled by

0273-1177/$36.00 Ó 2009 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2009.10.006

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

modeling the mascons as either a point mass or a disk mass placed at some predecided locations on the lunar surface (Wong et al., 1971; Ananda, 1975b, 1977), however, the higher degree and order of spherical harmonic models also capture this anomaly. The third problem has been eliminated by accurate Doppler tracking of the satellite at sub-millimeter level. The fourth problem has been tackled using different estimation techniques in the history of lunar research (Konopliv et al., 2001). A satellite-based gravity field modeling can be seen as an inverse problem formulation from the orbit perturbations due to all forces acting on the satellite. Four different ways of lunar gravity representation can be listed (a) spherical harmonic expansion that describes global deviations from a sphere (b) surface discrete mass distribution that describes surface/near surface mass anomalies like mascons, (c) Line of Sight (LOS) acceleration, and (d) elevation referred to a uniform density spheroid that is used to construct gravity maps. The spherical harmonic potential model (Kaula, 1966) which is mostly used for the gravity field modeling is given by " 1  ðnþ1Þ X n l X r ðC nm cosðmkÞ þ S nm sin ðmkÞÞ V ¼ R n¼0 R m¼0  ð1Þ  P nm ðsinð/ÞÞ ; ;

323

the mass of the spherical body and rest other symbols are same as in Eq. (1). Other discrete mass distribution models are also available in the literature, which are presented later. The LOS acceleration gravity modeling is used to plot the gravity map and is confined to near side of the moon where the direct observations are available. The fourth model (Sjogren, 1971) is used to construct gravity maps and referred to a uniform density spheroid is given as h¼R

N X n X

P nm ðsin /Þ

n¼1 m¼1

2n þ 1 ðC nm cosðmkÞ 3

þ S nm sinðmkÞÞ;

ð3Þ

where all the symbols have their usual meaning. Eq. (3) gives equivalent height variations from a spherical body of assumed uniform density. However, representation in terms of altitude requires the estimation of harmonic coefficients, and therefore, this is equivalent to spherical harmonic representation. Hence, only the first three gravity characterization schemes have been elaborated in context of various historical developments. In the next subsection an overview of gravity modeling, ranging from its foundation stone to some associated problems, has been described. 1.2. Lunar gravity modeling: an overview

where V is the potential function, R is the reference radius of the Moon, r is the distance between the point under consideration and the center of the heavenly body, l is the gravitational constant of the celestial body in consideration, k is the longitude of the point, / is the latitude of the point, n, m are the coefficients of summation, C nm ; S nm are the gravitational harmonic coefficients, and P nm is the associated Legendre polynomial. The main crux of this problem is the determination of coefficients C nm and S nm in Eq. (1) which describe the deviation of the Moon from a perfect sphere. The problem of determining the harmonic coefficients is fraught with difficulties. The number of undetermined harmonic coefficients grows dramatically with the degree and order of the potential model, and therefore, poses a problem of handling very large matrix inversion. Moreover, as the degree and order of model increases the effect of spherical harmonic coefficients show high correlation, i.e. their effects become indistinguishable, that create numerical problems in estimation by making the normal equation ill conditioned. The gravitational potential due to discrete surface mass model (discrete point masses near the surface) can be written as (Ananda, 1977)  N  l X i l ; ð2Þ V ¼ þ R i¼1 jr  r0 j where r and r0 are the position vectors of the satellite from the center of the celestial body and the point masses, respectively, and i is the ratio of the ith point mass to

The lunar gravity modeling has its root in satellite geodesy. Some early works in geodesy used the least squares technique and concluded that the Earth is not an ellipsoidal which led to the introduction of a shape based on the physics of the Earth: the equipotential surface of the gravity field. This was originally named the mathematical shape of the Earth by Gauss, was later renamed geoid by Floberghagen (2002). This proved to be a landmark in the history of geodesy and geophysics. The era of satellite geodesy heralded with the launch of Sputnik-1 on October 4, 1957. Various works done in 1960s laid a comprehensive foundation of orbit determination and estimation of other parameters, related to the motion of satellite in a suitable reference frame. Using special perturbation technique, extraction of gravitational parameters was realized. However, the early traces of satellite geodesy can be found in the work of Laplace, who in 1802 determined the dynamical flattening of the Earth from the motion of the lunar node, which was nothing but assuming the Moon to be an artificial satellite. A lot of work was done in the area of satellite geodesy prior to the launch of Sputnik-1. This provided a firm footing to the satellite geodesy and significant results were obtained soon after the launch of the first satellite. The success of this methodology can be judged from the fact that it provided an accurate determination of the Earth’s flattening, and thereby, the second spherical harmonic coefficient of the gravity field, from observations of Sputnik-2 and Explorer-1 (Floberghagen, 2002). With the launch of the Russian satellite Luna-10, the techniques

324

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

of geodesy got extended to selenodesy. Indeed, the results obtained from the Luna-10 mission in 1966, corroborated the fact that the oblateness of the Moon’s gravitational potential was larger than that estimated using hydrostatic equilibrium (Floberghagen, 2002; Akim, 1966). Moreover, it was concluded that lunar gravitational potential was pear shaped and also the center of figure and the center of mass did not match. These were some of the initial efforts towards the better understanding of the lunar gravitational characteristics. But the real quest to understand these characteristics commenced with the programs to land on the Moon. The first phase of efforts to obtain the Moon’s gravity model started with the Lunar Orbiter-1 and culminated with Lunar Orbiter-5 program (Lorell and Sjogren, 1968) with its launch on August 1, 1967 and later on using data available from subsatellites of Apollo 15 and 16. On the other hand, initial Russian efforts to model the gravity was concentrated on observation data available for Luna-10 (Akim, 1966). Although, these missions were mainly designed for the photographic purpose, the results from the analysis of the Doppler tracking data provided useful information about the lunar gravity field. It was hoped that these data would (1) resolve the long-standing question of whether or not the density increases toward the center of the Moon, and (2) establish the degree of north–south asymmetry (Lorell and Sjogren, 1968). These two questions were expected to be answered as soon as the numerical values of the 2nd and 3rd degree zonal harmonics were obtained. Some preliminary results were reported that included the values of these harmonics and concluded that the interior density of the Moon is considerably higher than the crust density (Lorell and Sjogren, 1968). Thereafter, continuous effort in this field yielded a number of new models especially during the late nineties and the beginning of this decade. In spite of these efforts, the existing models though very accurate but are still away from a truly global model. Some problems encountered in getting a truly global model has been briefly outlined in the following paragraphs. The main problem in the lunar gravity modeling is the lack of independent data such as in situ measurements to verify the gravity models obtained, therefore, the models developed till today cannot be termed as truly global models which can serve simultaneously both the selenodesy and selenophysical needs. Lack of far side data leads to unrealistic solution without a priori constraint. On the other hand, use of constraint makes the gravity results biased towards the constraint measure used. This problem is persisting for the lunar gravity modeling unlike for other planets. To circumvent this problem a number of regularization measures have been suggested but results will be as good as the regularization measure is. Further, the a priori constraint used as a regularization measure in the modeling of the lunar gravity is not tested through the in situ measurements to be true. Therefore, models recovered are very much dependent on the a priori constraint and the weights used to give relative weightage to the data. A detailed dis-

cussion on regularization measures has been elaborated in Floberghagen (2002). Efforts have been made to overcome such limitation by utilizing various ways of data weighting such as that formulated by Lerch (1991). It has been observed that not only the data weighting but also the method used for parameter estimation influences the quality of estimation, and in many cases divergence is also observed without a priori constraint. Therefore, the postestimation residual analysis is a key factor in assessment of the model performance for the orbit propagation purpose but must be cross checked with covariance and overlap analysis. And in fact, post-fit residual analysis is the best measure for the quality of estimation as far as orbit propagation is concerned. But such model may not be suitable for selenophysical purposes, as selenophysical modeling requires the local measures which are not reflected in the global gravity modeling for the good orbit propagation purpose due to high energy in the RMS magnitude spectrum. Still, models obtained for the orbit propagation purpose give reasonable information about the local features such as mascons and also reveal some internal structures such as crustal thickness etc. A brief review on lunar gravity modeling till 80s is available in Tuckness and Jost (1995), but it misses many of the important developments. In 2002, Floberghagen (2002) published a book on lunar gravimetry which provides an overview of lunar gravity modeling. Although, this does not describe many of the historical developments in the lunar gravity field modeling methodologies and the estimation techniques, but is an excellent treatment of the various issues related to the spherical harmonic model. This paper is basically a review work aimed at providing the modeling strategies, since 1960s till today, available in the lunar gravity field modeling. All the estimation techniques and gravity models including the recent ones have been included in this review. Therefore, the present work is divided into the various sections as follows. In Section 2, the parameterization schemes for gravity modeling and the formulation of normal equation for the extraction of gravity coefficients have been detailed. In Section 3, recent developments in gravity modeling have been discussed which essentially fall under spherical harmonic parameterization, and therefore, is an extension of Section 2 which deals with spherical harmonic parameterization. This has been done keeping in view interesting new results following rapid developments in computing technology. Besides, in Section 3 some new results from the analysis of existing models have been elaborated. Moreover, solution methodologies of normal equation have been discussed and the latest lunar missions such as SELENE and Chandrayaan-1 have been briefly presented. Finally, in Section 4, conclusion has been drawn. 2. Gravity modeling and data processing As discussed in Section 1, the gravity modeling process can be divided into three broad categories which depend

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

on the way the gravity field is parameterized, namely, the LOS acceleration modeling, the spherical harmonic modeling, and the discrete mass modeling. Integration of the equation of motion for propagating the state vector and the orbit determination are involved in all the three modeling strategies. However, only the latter two gravity modeling processes require the formulation of the normal equation using observation residuals in one or the other form, which is then solved using various strategies to extract the spherical harmonic coefficients or the discrete masses assumed in the perceived model. Thus the complete gravity modeling process consists of two stages: (a) parameterization of gravity model and (b) extraction of the parameters in the perceived model. The process of forming the normal equation or extracting the LOS acceleration involves integration of the equation of motion in its various forms which can be broadly divided into two categories: (1) direct method (short arc) and (2) indirect or averaged method (long arc) (Ferrari, 1977). The direct method involves direct integration of the equation of motion over short period of time while the indirect method involves integration of averaged equation of motion or use of averaged coordinates over a long period of time. Both these methods may use any of the gravity models such as the discrete mass model or the spherical harmonic model. The direct method involves weighted least square differential correction process with numerical integration of the spacecraft equations of motion, and the associated variational equations over a short period of time. The calculated values of the observable are derived from the numerical integration of the Cowell-type equation of motion of the spacecraft using a higher order predictor–corrector technique, starting with initial estimates of the spacecraft state at an initial epoch of a series of orbits. The partial derivatives, of the observables with respect to the set of parameters to be corrected, are obtained through numerical integration of the variational equations relating the change in the values of the observable to the change in the initial estimates of the parameters to be corrected. The differential correction procedure is applied to yield the estimates of the parameters through a number of iterations until the convergence is at an acceptable level. The indirect method is based on the results published by Kaula (1966) in mid-60s on the mathematical relations between the orbit of a satellite and its perturbations with respect to the orbital elements. He transformed the equation of motion in Cartesian variables to the equation of motion in orbital elements which is known as Lagrange’s bracket method (Roy, 1965). These equations were used later on to formulate the indirect method of extracting the spherical harmonic coefficients. The different ways of data processing that fall under this method were motivated by minimization of the amount of computation involved, accuracy achieved, and capturing the long term orbit evolution. Therefore, some were based on the integration of the Lagrange’s equations while others used numerical dif-

325

ferentiation to solve the problem. A brief outline of the indirect method has been presented in the following paragraphs. The indirect method uses either integration of the averaged Lagrangian equation to propagate the averaged orbital elements of the satellite over a long period of time or obviates the need of integration by fitting splines to the averaged orbital elements and then taking their derivatives in order to extract the harmonic coefficients. The latter one is numerically more efficient and uses Cartesian coordinates to propagate the satellite elements over one orbit for orbit determination. Mean orbital elements are calculated by averaging the osculating elements over many orbits and thereafter smoothed with patched cubic polynomials to determine their time derivatives to obtain the rate of orbital elements. A least squares fit is done to these rates after forming normal equation in order to estimate the spherical harmonic coefficients. The indirect method best captures the long period perturbation and therefore, estimates the lower degree and order coefficients with good precision. However, problems may be encountered while fitting to find the Kepler elements for orbit of low inclination or equatorial orbits where X and x may become indeterminate even after use of other combinations of orbital elements. Therefore, results obtained using low inclination orbits get affected. While doing such estimation it is necessary to estimate more number of coefficients than what are used to propagate the orbit. Any effort to estimate less number of coefficients results in aliasing. This model involves fitting to normal points obtained which obviously are not many in numbers, besides lacking in global coverage, and therefore, any estimate obtained will contain correlation for higher degree and order. Hence, the orbit propagation using such models produces secular drift for long time orbit propagation. With the advancement in computing speed, the Lagrange’s method of orbit propagation in terms of the orbital elements is sparingly used. Rather, using the higher order numerical integration techniques orbit is propagated very precisely in time in rectangular coordinates and the method elaborated by Ellis (1980) or some other suitable method (Lerch, 1991; Minter and Cicci, 1997; Tscherning, 2001) is used to determine the harmonic coefficients. The spacecraft states and other parameters are estimated using a square root information weighted least squares filter (Bierman, 1973; Lawson and Hanson, 1995) in selenocentric reference frame which is defined with respect to the coordinate system defined by the Earth’s mean equator at the epoch of J2000 (now the Geocentric Celestial Reference System (GCRS) the geocentric realization of International Celestial Reference System (ICRS) has replaced the reference frame of J2000). The global and the local parameters are separated out in the square root information matrix processing and observations are processed sequentially to save the computational effort and memory. This technique was illustrated by Kaula (1966) using partitioned matrices and has been further elaborated by Ellis (1980).

326

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

The historic developments of lunar gravity modeling can be listed as given in Tables 1–5. In the following subsections various gravity parameterization methods and their outcomes are described briefly. 2.1. LOS gravity modeling The LOS gravity modeling involves direct method of data processing. In the 60s, estimation of higher degree and order spherical harmonic coefficients posed formidable problem considering the computational power and amount of data then available. Therefore, Muller and Sjogren (1968) instead of doing spherical harmonic modeling directly extracted the local acceleration, for plotting the gravitational map of the near side, of the satellite by processing the Doppler residuals of the observations of the Lunar Orbiter-5 spacecraft taken from a number of orbits. A triaxial lunar model was used and gravitational perturbation of planets was included. Doppler residuals were generated after orbit determination was done to which a patched cubic polynomial was fitted and derivatives were then extracted to estimate the LOS acceleration. The LOS acceleration was then normalized to represent the acceleration at an altitude of 100 km under the assumption that the typical mascons were at a depth of 50 km below the surface, that is, 2 anorm ¼ acomp  ðh þ 50Þ =1502 ;

ð4Þ

where anorm is the normalized acceleration, acomp is the computed acceleration, and h is the altitude in kilometers. Various relations for finding the velocity residuals from the Doppler residuals are given below. _ ð2qÞ ; k Doppler residual  c ; Dq_ ¼  2  frequency

fd ¼ 

ð5Þ ð6Þ

where fd is the Doppler frequency shift, k is the wavelength of the tracking signal, c is the speed of light, and q_ is the range-rate. Using the extracted acceleration Muller and Sjogren concluded that large rate of change in acceleration over the five ringed maria (Imbrium, Serenitatis, Crisium, Nectaris, and Humorum) indicated presence of small mas-

cons (50–200 km) (large mascons will cover large area and therefore acceleration change will be slow). It was further concluded by them that this was a confirmation of the Urey–Gilbert theory of lunar history which predicted the presence of such large-scale mass concentration below these maria, which they termed as mascons (Muller and Sjogren, 1968). This method was a breakthrough at that time and provided big impetus to research on lunar gravity modeling. However, LOS gravity not being the vertical acceleration exactly, and correction for the least squares filtering effects (Gottlieb, 1970) not being included which reduces the larger amplitudes (both positive and negative) by 20–30%, the produced gravity model did not represent the gravity field exactly. This is evident from the fact that the LOS gravity model derived from residuals based on spherical harmonic model of LP165P truncated at 110° (Tables 1–3 in Konopliv et al. (2001)) shows differences in magnitudes for different mascons though their sign matches, thereby, confirming the shortcomings in this model. This inconsistency is expected because while carrying out the acceleration normalization at an altitude of 100 km it was assumed that the mascons were located at a depth of 50 km below the surface which was an ad hoc assumption. Thus equating the LOS acceleration with the vertical gravity and ad hoc assumption of the location of the mascons below the surface caused impaired values for the gravity maximum and minimum obtained. Thus even though the LOS gravity revealed many surface features but magnitude of these features do not conform to the values predicted in Konopliv et al. (2001). Moreover, this could not be used for long term orbit propagation being in the form of gravity map. Therefore, need for a more consistent model was felt which could be utilized for orbit propagation and would simultaneously provide the local gravity fields and surface features besides revealing the interiors of the Moon. This led to the development of the spherical harmonic and the discrete mass models which are described later. In continuation of efforts along this direction, with the availability of data from Apollo 14–17 new LOS gravity models were produced. The technique used in Muller and Sjogren (1968) was used for the data processing and getting LOS acceleration. When the Apollo 14 was in low periapsis

Table 1 Lunar gravity modeling using LOS parameterization. Reference

Satellite data; gravity model; method of state propagation

Results

Muller and Sjogren (1968)

Lunar Orbiter-5; triaxial model; direct method

Sjogren et al. (1972a), Sjogren and Muller (1972b) Muller et al. (1974) Sjogren et al. (1974a)

Apollo 14 and 15, respectively; triaxial model; direct method

Line-of-sight acceleration derived from the Doppler measurements; discovery of lunar mascons; local gravity plot Line-of-sight acceleration

Sjogren et al. (1974b,c)

Apollo 15; triaxial model; direct method Apollo 15 and 16 subsatellites; low degree spherical harmonic model; direct method Apollo 16 and Apollo 17, respectively; direct method; overlaps with discrete mass model

Line-of-sight acceleration Line-of-sight acceleration Line-of-sight acceleration

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

327

Table 2 Lunar gravity modeling using spherical harmonics. Reference

Satellite data; gravity model; method of state propagation

Results

Lorell and Sjogren (1968) Lorell (1970)

Lunar Orbiters 1–4; spherical harmonic model; indirect method Lunar Orbiter 1–5; spherical harmonic model; indirect method Lunar Orbiter 1–5; spherical harmonic model; indirect method

4  4 spherical harmonic solution and coefficients C50, C60, C70, and C80, a total of 25 coefficients; homogeneous moon concluded from C20 and C22 coefficients 8  4 spherical harmonic solution; tesseral harmonics have negligible effect on the models of lower order and degree (n < 6), finite difference method used 15  8 spherical harmonic solution (full 8  8 and zonal harmonics through 9–15); predicted orbital parameters better than any model till that date; gravity map of the moon; some far side features resolved 4  4 spherical harmonic solution; gravity map; identified by Hertzsprung and Korolev 13  13 spherical harmonic solution; using low altitude data derived coefficients up to 13th degree and order; lunar maria Imbrium, Serenitatis, Crisium confirmed

Liu and Laing (1971) Sjogren (1971) Michael and Blackshear (1972) Ferrari (1972) Bryant and Williamson (1974) Ferrari (1975)

Blackshear and Gapcynski (1977)

Lunar Orbiter-4; spherical harmonic model; direct method Lunar Orbiters 1–5; spherical harmonic model; direct method Lunar Orbiter-1,2,3,5; spherical harmonic model; indirect method Explorer 49; spherical harmonic model; indirect method

4  4 spherical harmonic solution; gravity map; equipotential surface revealed trial axial mass distribution with large areas of mass concentration 3  3 spherical harmonic solution

Apollo-15, 16 and Lunar Orbiter-5; spherical harmonic model; indirect method

16  16 spherical harmonic solution; solution arrived at using normal point elements by fitting a patched cubic polynomial. Long term applicability is questionable as selenographic system induces small perturbations in the mean element of the satellite which are not consistent with the conventional form of the Lagrangian equation J 2  J 6 zonal coefficients only; lunar moment of inertia determination

Explorers 35 and 49; spherical harmonic model; indirect method

Table 3 Lunar gravity modeling using spherical harmonics. Reference

Satellite data; gravity model; method of state propagation

Results

Ferrari and Ananda (1977) Ferrari (1977)

Lunar Orbiter-5 and Apollo 15, 16 subsatellites; spherical harmonic model; indirect method

16  16 spherical harmonic solution; tried to model far side gravity

Lunar Orbiters-5 and Apollo 15, 16 subsatellites; spherical harmonic model; indirect method

Ferrari et al. (1980)

Doppler tracking data from LO-4 and lunar laser ranging data; spherical harmonic model; direct method

Akim and Vlasova (1977, 1983) Konopliv (1993)

Luna 10, 12, 14, 15 and 22 for 1977 model; spherical harmonic model; direct method

16  16 spherical harmonic solution; showed correlation between far side surface features (Moscoviense, Mendeelev, Korolev) and surface gravity anomalies 5  5 spherical harmonic solution; lunar principal moment; lunar gravitational constant; earth/moon mass ratio; obliquity of lunar pole; selenocentric coordinates of lunar retro reflector etc 7  7 spherical harmonic solution; could not produce any localized anomaly

Lunar Orbiters 1–5 and Apollo 15, 16 subsatellites; spherical harmonic model; direct method

Lemoine et al. (1994) Lemoine et al. (1997)

Lunar Orbiters 2–5, and Apollo 15 subsatellite, Clementine; spherical harmonic model; direct method Lunar Orbiters 1–5, Apollo subsatellites 15 and 16, Clementine; spherical harmonic model; direct method

Konopliv et al. (1998)

Lunar Orbiters 1–5, Apollo subsatellites, Clementine, Lunar Prospector; spherical harmonic model; direct method Lunar Orbiters 1–5, Apollo subsatellites, Clementine, Lunar Prospector; spherical harmonic model; direct method Lunar Orbiters 1–5, Apollo subsatellites, Clementine, Lunar Prospector; spherical harmonic model; direct method

Konopliv and Yuan (1999) Konopliv et al. (2001)

60  60 spherical harmonic solution (Lun60d); surface features obtained in agreement with the previous line-of-sight reductions, far side features revealed, good orbit prediction as compared to all earlier models, high degree terms have excessive power, gravitational anomalies mapped on to the lunar surface showed artifacts 70  70 spherical harmonic solution (GLGM-1); internal structure of the moon revealed 70  70 spherical harmonic solution (GLGM-2); confirms major features of gravity fields from earlier researches but also reveals significantly more details 75  75 spherical harmonic solution (LP75D and LP75G); reveals some mascons on far side earlier identified to be negative features 100  100 spherical harmonic solution (LP100J and LP100K); reveals some mascons on far side earlier identified to be negative features 165  165 (LP165P multi-step), 150  150 (LP150Q one-step) spherical harmonic solution

328

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

Table 4 Lunar gravity modeling using discrete mass. Reference

Satellite data; gravity model; method of state propagation

Results

Wong et al. (1971) Ananda (1977)

Data for point mass model – Lunar Orbiter 4 and 5, for disk mass model Lunar Orbiters 1–5, Apollo 8, and Apollo 12; harmonic model and point or disk mass model; direct method Apollo 15 and 16 subsatellite, Lunar Orbiter 5; used 117 point masses; indirect method

Discrete (point or disk model) + spherical harmonic modeling; discrete masses mapped to 20  20 spherical harmonic model; confirmed the results obtained by Sjogren (1971) Gravity map resolved all the near side mascons, in agreement with LOS model; 117 discrete masses were mapped to 20  20 spherical harmonic coefficients

Table 5 Lunar gravity models using hybrid method. Reference

Data/model used

Results

Bills and Ferrari (1980) Sagitov et al. (1986)

Lunar laser ranging data, Lunar Orbiters 1–5, Apollo 15, 16 subsatellites, Apollo 8–12; hybrid model

16  16 spherical harmonic solution; provided high accuracy gravity map of the near side and also shows reasonable correlation with the most significant mascons on the near side 16  16 spherical harmonic solution; results similar to Bills–Ferrari model

Based on five different data sources including Luna and Apollo 14, 15, 16, 17; hybrid model

orbit (16 km) the gathered data had high frequency information about the lunar gravity. Based on the available data from Apollo 14, Sjogren et al. (1972a) produced a LOS gravity model, and another model from Apollo 15 tracking data (Sjogren and Muller, 1972b). They concluded that Theophilus, Hipparchus, and Ptolemaeus are negative features whereas Mare Nectaris is a large positive region indicative of a broad disk near the surface rather than a deeply buried spherical body. These features are also confirmed in the latest studies (Konopliv et al., 2001) proving the validity of the results obtained by Sjogren et al. (1972a). It was further concluded by them based on a surface disk mass simulation that the acceleration profiles (simulated) at the Mare Serenitatis and Mare Crisium locations matched with that of the observed LOS gravity profile while a simulation based on an embedded disk showed a shift in acceleration profile at the Mare Serenitatis and Mare Crisium locations confirming that lunar mascons are essentially surface features with a mass distribution of approximately 800 kg=cm2 unlike the assumption made by Muller and Sjogren (1968) in normalizing the LOS gravity that mascons are located at a depth of 50 km. Further some new analysis of data from Apollo 15 (Muller et al., 1974) and 16 (Sjogren et al., 1974a,b) confirmed the results obtained in the previous analyses of Apollo 14 observations. The modeling technique for this model differed from the previous ones in that surface mass distribution was used to model the limb areas of the Moon. Therefore, there is a overlap of the LOS gravity modeling and the discrete mass modeling to be discussed later. They concluded that volcanic type features contributed significantly in forming the local gravity highs. The sharpest gravity lows were associated with the large unfilled craters. Results due to Konopliv et al. (Table 3 given in Konopliv et al. (2001)) also show these craters as the centers of gravity lows and were formed perhaps due to the ejection of mass following an impact of the asteroid with the lunar surface. Data analysis from Apollo 17, Sjogren et al.

(1974c) confirmed some anomalies like Grimaldi, Copernicus, Sinus Aestuum, Apennines, Mare Serenitatis, Littrow and Mare Crisium. However, shift was seen in the optical centers of these anomalies and the acceleration derived contours. Sjogren et al. (1974c) concluded that this was perhaps due to equating the line of sight acceleration with the vertical acceleration. However, a major lacking of such representation is the artificial fixing of the position of the discrete masses. Even though the results of this approach conforms to some of the surface features on the limb area but may not be taken to be the final reality. Moreover, the LOS gravity being the local feature and not having a representation in terms of equation limits the use of LOS method to selenophysical purpose only. However, these models led to a better understanding of the lunar gravity. From the above observations it can be concluded that although LOS gravity models revealed the near side features but lacked severely in the magnitude determination and moreover, their locations differed slightly from the present locations reported in Konopliv et al. (2001) which is based on large amount of data and precise modeling. 2.2. Gravity modeling using spherical harmonic coefficients Gravity modeling using spherical harmonic coefficients has been outlined in the introduction. In fact, this modeling strategy is the mostly used method in the recent times. The modeling strategy and the results obtained using this method, since nineties, are elaborated in Section 3 including some new analysis of the existing gravity data. A brief review of the gravity modeling in terms of spherical harmonic coefficients which were carried out till the end of 80s are presented in this subsection. Although gravity modeling in terms of spherical harmonic coefficients require solution of the normal equation based on observation residuals but is preceded by integration of the equations of motion using direct or indirect method. Therefore, this

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

modeling strategy has been elaborated under the headings of direct and indirect methods as follows. 2.2.1. Direct method The spherical harmonic solution for lunar gravity modeling was first attempted in 1968 (Lorell and Sjogren, 1968), however, indirect method of data processing was used for this work. The first effort to use the direct method for lunar gravity modeling was reported in 1971 by Sjogren (1971) who integrated the equation of motion and then estimated spherical harmonic coefficients up to 4th degree and order using 727 Doppler observations from Lunar Orbiter-4 high altitude near polar orbits (2700–6000 km). The coefficients thus estimated did not have any signature of high frequency because of high altitude, and this allowed the model to fit to the noise level of 1 mm/s as compared to the models derived from low spacecraft altitudes (100 km) data where systematic residuals of tens of millimeters per second occurred. Sjogren (1971) reported that the results obtained were consistent with the solution obtained by Langley Research Center (LRC), 13th order estimate (Michael et al., 1969), and the JPL 15th order zonal and 8th order tesseral model (Liu and Laing, 1971). He further concluded that solutions for C 20 and C 22 were in good agreement with other independent estimates, and therefore, implied a homogeneous Moon. Besides, summation of the coefficients C 20 ; C 22 ; C 31 and C 33 showed 572 m bulge at the equator on the far side and a 408 meter bulge on the near side. He concluded that this was due to the far side positive anomaly spread out along the equator. These results were derived based on the uniform density moon which does not hold in reality. Values of bulge on the near side and the far side remained uncertain because it requires conformation with the topographic study especially in the view of large uncertainty in the estimation of harmonic coefficients. The post-fit residuals were at the level of 1 mm/s but its capacity for mission planning cannot be assessed because of lack of any data on long term orbit propagation. If a comparative study with the LP165P model (Konopliv et al., 2001) coefficients is done it shows some offset in the values of the above-mentioned coefficients. This may be due to the improved quality of tracking, modeling, and use of large number of spherical harmonic coefficients in the LP165P model. Further improvement in the coefficients modeling through 13th degree and order was reported, in 1972, by Michael and Blackshear (1972) based on 20,148 Doppler tracking data spanning over 80 days of tracking coverage from Lunar Orbiters 1–5 and using direct integration of equation of motion and associated variational equation, as an extension of their work in 1970 (Michael et al., 1970). Michael and Blackshear (1972) reported that the gravity contour plotted for an altitude of 100 km above the mean lunar surface confirmed the lunar maria Imbrium, Serenitatis, and Crisium, however, the contour in the polar region and for the far side showed correlation with the regions of less dense data coverage indicating that con-

329

tour plot for the far side and polar regions are not much reliable. This problem is merely an indication of the inaccuracies in the model extracted and may be due to a small number of harmonic coefficients used to model the gravity which consists of an infinite series and the lack of wide coverage especially the far side. Finally, in 1980, a 5  5 spherical harmonic model which was really meant for the geophysical parameters estimation of the Earth–Moon system was published by Ferrari et al. (1980) using direct method through simultaneous reduction of lunar laser ranging and Lunar Orbiter-4 tracking data unlike earlier separate data processing which were combined later to get the estimates of inertia. The value of the principal polar moment from this model is C=MR2 ¼ 0:3905  0:0023 as against that from model LP75G of 0:3932  0:0002. The high uncertainty in the model due to Ferrari et al. may be due to the small amount of observations used as compared to the LP75G model. Basically, this work confirmed the results of Gapcynski et al. (1975) who derived the moment of inertia of the Moon using indirect method based on data from Explorers 35 and 49. Besides this, research done by Soviet scientists Akim and Vlasova (1977, 1983) led to a 7th degree and order gravity model based on all-Soviet-data. Even though this model was independent of other models, but the low inclination of the Luna spacecraft as well as the low degree and order of the model made it unable to reproduce even any localized anomalies. In fact, the direct method of data processing which processes short arc of observations is dominated by the short period perturbations caused by the nearside surface anomalies, and therefore, lack in the information about the far side surface features especially if the polar region coverage is poor. Hence the models extracted using this method before 80s suffered from unrealistic revelation about the far side surface features. Besides, unmodeled propulsive forces used for the attitude maneuver also affected the estimated spherical harmonic coefficients. Therefore, the models extracted using this method though successfully revealed many near side surface features but failed to predict the far side features. Besides, they could not be used successfully for the long term orbit propagation because of lack of precision in the estimation of low degree and order coefficients which dominate the orbit evolution over a long period of time unlike the models extracted using the indirect method to be discussed in the following subsubsection. 2.2.2. Indirect method The indirect method was devised for accurately extracting the low degree and order coefficients in the gravity model for mission planning purpose. The indirect method, as elaborated earlier, tries to average out the long term evolution of the satellite orbit which are affected by the low degree and order spherical harmonic coefficients. The first attempt to model the lunar gravity using this method was reported, in 1968, by Lorell and Sjogren (1968) who pro-

330

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

cessed the observation data from Lunar Orbiters 1–4 applying the method of averages (Lorrel et al., 1964) to extract the fourth degree and order coefficients including C 50 ; C 60 ; C 70 , and C 80 summing together a total of 25 coefficients. It was shown by them that C 20 was 3% less than that for a homogeneous Moon indicating interior density moderately higher than crust density. However, the amount of data used in processing being not sufficient especially due to incomplete coverage the coefficients thus estimated were not very accurate. A comparison of the normalized values of these coefficients with the LP165P shows offsets. For example reported value of normalized C 20 by Lorell and Sjogren (1968) is 0:9062  104 while that from LP165P is 0:9089  104 . Similarly the normalized value of the C 22 is 0:3394  104 while that from LP165P is 0:3463  104 . Other coefficients also differ a lot. These differences are enough to cause large error in long term orbit propagation. However, considering the amount of data used in this model it seems that the results obtained were commendable as a first effort which laid the background for further improvement. In 1970, Lorell (1970) published results which were improvement over his earlier work done in 1967 and 1968 and differed from them in that the method of numerical integration of the equation of motion using Lagrange’s form and the harmonic expansion derived by Kaula were applied for the formulation. Averaging method was applied over the two orbits to get the average value of the orbital element as described in Bogoliubov and Mitropolosky (1961). He estimated coefficients up to fourth-degree tesseral and eighth-degree zonal harmonics and concluded that tesseral harmonics have negligible effect on the models of degree lower than 6. In 1971, Liu and Laing (1971) extended the degree and order of the estimated model to 8  8 including zonal from degree 9 to 15 for lunar mission navigation purpose using data from all the five Lunar Orbiters. This was an extension and refinement of the method used by Lorell (1970), but unlike Lorell, Liu and Laing used variational equations instead of finite differences in order to complete the least square fit (Liu and Laing, 1971). To overcome the problem of large scale data handling, the radar data were compressed into normal points consisting of five mean orbital parameters a, e, i, X, and x averaged over an anomalistic period. Normal points were associated with a weighting matrix which contained information about statistics and correlations. They generated a global model for the Moon which was consistent with the gross features derived from physical libration data (Campbell et al., 1969) and which did not contradict the fine features (mascons) derived by direct mapping of the Doppler data by Muller and Sjogren (1968). To conclude, an independent check of the model’s validity was its ability to predict orbital parameters better than any model till that date. In fact, Liu and Laing reached the limit of improvement that could have been achieved using this method and the available Lunar Orbiters data and was the best model till that date. However, due to its inability to predict real orbits over a longer per-

iod and a very limited correlation with surface features, the model was declared unrealistic later on. Thus even though the indirect method is designed to capture coefficients affecting the long term evolution of the satellite, but due to limited coverage especially on the far side and approximations introduced while applying the method of averages, and small number of normal points thus obtained made the estimated coefficients not very realistic and therefore, these models could not predict long term orbit and also failed to satisfy the resolution required for selenophysical study. Further research during seventies and early 80s were motivated by the requirement of long term orbit propagation which resulted in various modeling strategies to reduce the amount of computation involved. Along this line, in 1972, Ferrari (1972) worked out a 4  4 spherical harmonic model where in he assumed a solution to the Lagrange’s equation valid for a day to be six dimensional time series given as ~kðtÞ ¼ K e 1t þ K e 2 t2 þ d~k  þ d~k  þ d~k sr ; e0 þ K

ð7Þ

e 1; K e 2 are e 0; K where ~k is the vector of orbital elements, K Keplerian constants determined in the first step of least squares fitting of Doppler tracking data, and d~k  ; d~k  ; and d~k sr are the orbital element variation due to the Earth, sun gravity and solar radiation pressure perturbations. These values were calculated using numerical integration. Once the Keplerian parameters were available through least squares estimation, Eq. (7) was differentiated to obtain the rates. The lunar gravity field determination was performed using 199 normal points in a second weighted least squares processor which used as input Kepler element rates, Kepler elements, and a weighting matrix (Ferrari, 1972). To conclude, the Keplerian parameters estimated in the fist step showed high correlation which implies sub-quality of the estimates. This was perhaps due to the use of big chunk of data coming from single ground station tracking which makes the separability of orbital elements in estimation process difficult. It is expected that the coefficients estimated in the second step would have got affected due to estimates of normal points in the first step. Thus this model could not prove to be fit for long time orbit propagation. The equipotential surfaces from this gravity field was indicative of a triaxial moon and indicated three major areas of mass concentration the bigger one lying on the front side while the smaller two on the far side. In the perspective of the lunar gravity map by Muller and Sjogren (1968) the results of this model seem to be unrealistic and this is expected as this model was conceived for the long term orbit propagation. This model was followed by a 3  3 model presented by Bryant and Williamson (1974) using 234 days of Explorer-49 data from 1065 km altitude orbit and was used to improve the lunar moment of inertia obtained by Gapcynski et al. (1975) processing 600 days of the time variation of ascending node of the spacecraft (Blackshear and Gapcynski, 1977). Use of this model in orbit propagatin is not known.

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

From 1975 to 1977, Ferrari (1975, 1977) and Ferrari and Ananda (1977) published results based on data for Apollo sub-satellites obtained from Johnson Space Center (JSC). Using his earlier methods (Ferrari, 1972) he improved the gravitational model up to 16th degree and order and the corresponding coefficients can be found in Ferrari (1975, 1977). Ferrari chose to use the orbital elements in describing the perturbations but obviated the need to integrate Lagrange’s equation (Ferrari, 1972, 1973) in order to avoid the approximation introduced in calculating the orbit on the far side. He differentiated the normal points (compressed data) elements by fitting a patched cubic polynomial (Lawson and Hanson, 1974) and obtained the rates (this ensured effects of unmodeled spacecraft maneuver is minimized) which modeled the effects of harmonic coefficients (Sun, Earth and solar radiation pressure effects were removed), where the mean orbital elements were obtained using numerical integration of the equation of motion in Cartesian form. The results were encouraging reflecting a good fit for all the elements, thereby, represented a good lunar gravity model. But unlike his previous work (Ferrari, 1972) he used a solution of type   1  1   C  ð8Þ p þ F T W 1 k_ ; p ¼ F T W 1 F þ C1 where W is a weighting matrix for each set of elements rates (assumed diagonal), C is the a priori covariance associated with the initial estimates  p of the harmonics, and k_ is a vector of mean element rates. This equation stems from minimizing the performance index defined by Eq. (9) with a priori constraint.  T   p W 1 k_  F  p : ð9Þ 2 ¼ k_  F  In order to enhance the numerical stability of the least squares solution, the square root form (Bierman, 1973) of these equations was used. The a priori statistics for generating C were those given by Kaula (1969) and described by Eq. (10), and the a priori estimates of p were assumed to be zero. 3:5  104 N ðn; mÞ ; rðC nm ; S nm Þ ¼ n2

ð10Þ

where N ðn; mÞ ¼

ðn  mÞ! ð2n þ 1Þð2  d0m Þ ðn þ mÞ!

1=2 ;

ð11Þ

and d0m is the Kronecker delta. This essentially caused the gravity potential due to the harmonics to have a 1=n2 spectrum as did those of the Earth. They showed through simulation that it is possible to extract the far side features by fitting to mean elements derived from near side tracking (Ferrari and Ananda, 1977). This was carried out separately using harmonic modeling and another one using five point masses first to generate the orbit and thereafter, recovering the spherical harmonic coefficients and the masses, respectively (Ferrari and Ananda, 1977). However,

331

validity of discrete mass simulation depends on the separability of the effect of discrete masses at the altitude of the orbit. In reality, the wide separation used for locating the discrete masses does not hold true in reality and therefore, though simulation showed positive outcome but the results on real data are not be very promising. The resolution of far side features from the near side observation data were similarly shown to be possible through simulation (Ferrari and Ananda, 1977). However, this possibility is controlled by a number of factors such as interval of data sampling, number of data points in the sample. Sampling rate determines the possible resolution in terms of degree and order of the spherical harmonic coefficients. Hence the model extracted by this process will be restricted by the sampling rate which in turn will limit the resolution up to which features can be extracted. Generally, the feature extraction of areas out of direct coverage is an ill conditioned problem and therefore, even though some of the features on the far side could be resolved does not imply the results obtained to be reliable. This is corroborated from the fact that the features such as Korolev, Mendeelev, Moscoviense which were found to be negative features in the model extracted by Ferrari (1977) turned out to be positive features in the LP165p model (Konopliv et al., 2001). Thus the theoretical simulation by Ferrari and Ananda (1977) though showed positive outcome about extraction of far side model from the nearside coverage and the results due to Ferrari (1977) supported this fact but has been negated by the results of LP165P model (Konopliv et al., 2001) (even spherical harmonic models of late nineties showed these features to be negative ones). Besides, the selenographic system is known to induce small perturbations in the mean elements of the satellite (Ferrari and Heffron, 1973) which are not consistent with the conventional form of the Lagrangian equations making the wide and long term applicability of indirect method questionable. 2.3. Discrete mass modeling The results obtained by Muller and Sjogren (1968) and others showed that harmonic coefficients were not well defined as these were correlated, while the gravity field appeared to be well defined, Therefore, need of local gravity model was felt. Following the efforts of Gottlieb (1970) and Sjogren (1971) to use mascon model, Wong et al. (1971) developed a discrete mass model based on the comprehensive foundation of discrete mass modeling already put forward by Koch and Morrison (1970) in 1970 where in the discrete masses magnitudes were differentially corrected to minimize the quadratic sum of the Doppler residuals. The equation of motion was integrated using 8th order predictor corrector Gauss–Jackson differencing method (aerospace orbit program Trace-66) and is given as ! N 1 X €r ¼ rV ¼ Gmr mi F i ; þ ð12Þ r i¼1

332

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

where Gm is the product of gravitational constant and lunar mass, mi stands for the magnitude of ith mass in units of lunar mass, r is the position vector of the satellite in selenocentric reference frame, and F i denotes the ith potential function due to either the point mass or the disk as discussed later. Here, the masses mi were unknown and to be determined through differential correction using the Doppler data. The differential correction coefficients were obtained by integrating a set of variational equations given as follows " # N X 1 r ¼ Gm dðr Þ þ r dmi F i : ð13Þ d€ r i¼1

where DC is the observation residual vector, R is the data variance matrix (assumed diagonal), DS denotes the initial state vector deviations and R0 is the a priori state variance diagonal matrix. This method can be categorized under direct method. In fact, Wong et al. (1971) discussed two models: the point mass model (1969 solution) and the disk mass model (1970) solutions. Later in 1977, Ananda (1977) published results based on 117 point masses using indirect method and mapped the model thus obtained to 20  20 spherical harmonic coefficients which are described briefly in the following sub-subsections.

only polar orbits data seem to be the prime reason for the under performance of this model. Later efforts were made by Ananda (1977) to incorporate the new data from Apollo mission. He obtained the mean orbital element rates as in Ferrari and Ananda (1977) but the second step of Ferrari of deriving element rates was achieved using equations of motion including the perturbation functions with 117 point masses (Wong et al., 1971) (this he has elaborated in his paper Ananda, 1975b). An alternate form of Lagrange’s equation was used by applying the averaging method developed by Bogoliubov and Mitropolosky (1961) to average the partial derivatives of the perturbing function potential and substituting in Lagrange’s equations (Kaula, 1966). The averaging was done over the period of the orbit and this was achieved by averaging the partials with respect to the mean anomaly and thus eliminating all short periodic terms. He concluded that mapping of 117 discrete masses to 20  20 gravity field showed agreement for low order coefficients with the earlier models. However, the gravity maps obtained for the two sides though described the known main features but they lacked precision and the contours were scattered. This was supported by the lack of similarity of the coefficients with the past research except some lower order coefficients. The artificial framing of the model using point masses placed at predecided locations could not make this model promising in spite of the fact that the extra data for Apollo 15 and 16 subsatellites were in low inclination orbit. Thus point mass modeling seems to be a naive approach for the complicated gravity field modeling especially for the Moon where the lack of far side data draws the limit.

2.3.1. Point mass model The point mass dynamical model consisted of point masses located at chosen grid points, on the mean lunar surface. Both the initial states of each arc and the discrete masses were estimated to avoid effects of any state error which might have got aliased with the discrete masses in spite of the fact that initial state components and point mass parameters have very different Doppler perturbation signature (a mascon produces a high-frequency, localized oscillation in contrast to the smooth, long wavelength effects of initial-position or velocity changes). Because of uncertainties in the estimation of position and the velocity components due to Earth-based Doppler observation over short arc, a priori constraints were imposed to ensure fast convergence of estimation process. The radial acceleration from the point mass model at an altitude of 100 km confirms some of the major mascons on the near side, however, there is a shift of position and magnitude if these are compared to the LP165P model. This model fitted best to the polar orbits and this is expected since data from Lunar Orbiters in polar orbits were used to extract this model. However, it failed for the equatorial orbits as a consequence of the limited data used as mentioned above. Thus artificial positioning of the discrete masses and the use of

2.3.2. Disk mass model The field on the surface of moon due to point mass becomes singular indicating one of the shortcomings of the point mass model. To avoid the singularity problem due to the point mass model the simplest solution was moving the masses below the earth surface and this did not demand change in the computer program, however, Wong et al. (1971) preferred disk mass model which also conformed to the finite size of the mascons as postulated by Kaula (1969). For large distance the disk mass served as the point mass. The Root Mean Square (RMS) data fit to the low altitude equatorial orbits (60 km periapse altitude) gave much less Doppler residuals indicating better modeling of gravity as compared to the point mass model. They observed that the residuals of polar orbits from the disk solution were degraded relative to the point mass model, while the equatorial orbit residuals were much improved. Besides, they found that, in general, the point mass model always fitted better than the triaxial model on polar orbits and worse on equatorial ones; whereas, the disk model did better than the triaxial model for both polar and equatorial orbits. The improved results from the disk mass model is expected because the data used for this model involved Lunar Orbiter data which ranged

From the computed variations dr, the differentials of the observables were obtained using the differentiation chain rule and was technically consistent because over short arcs of about 1 h, numerical computation errors were essentially quite small as compared to the data noise. The quantity minimized was the quadratic sum Q ¼ DC T R1 DC þ DS T R1 0 DS;

ð14Þ

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

from 12° to 85° orbit inclination besides data from Apollo 8 and 12 in 12° and 15° inclination orbits, respectively. Thus not only the orbits of low inclinations were added but also the amount of data used was around double as compared to the point mass model. Of course the disk model is a better representation of the mascon, however, the above-mentioned factors also played role in the improved result using disk model. In the following paragraph composite model of lunar gravity is briefly outlined to make the review of discrete mass modeling complete. Wong et al. (1971) combined spherical harmonic expansion for the back side of the Moon by Michael et al. (1969) and the discrete mass model for the front side and proposed composite model. The resultant field was then analyzed to obtain a set of spherical harmonic coefficients using the observational data. They observed that the coefficients derived did not perform well in orbit prediction because for the long term prediction the lower order coefficients are required most accurately which was not the case for this model as these were based on short arc limited data. For the low degree expansion, they adopted the coefficients derived by Liu and Laing (1971). The model was capable of predicting ephemerides of low altitude equatorial orbits for a day with Doppler residual of less than 20–25 Hz. To conclude it appears that discrete mass modeling with arbitrary framing of masses and far side modeling using spherical harmonic coefficients for which the modeling uncertainty is lot due to lack of far side data does not add much advantage for extracting the spherical harmonic coefficients accurately. 2.4. Hybrid model From the discussions in the previous subsections it is obvious that different gravity parameterization schemes with direct or indirect methods of data processing did not perform satisfactorily for the both the selenophysical and the orbit propagation purposes. Therefore, to overcome this shortcoming Bills and Ferrari (1980) combined data from different parameterization schemes with different data processing methods to take advantages of each data type and extracted a 16th (16th degree and order spherical harmonic model). The first data set was obtained by relating the variations in Keplerian rates to the gravity parameters as proposed in Ferrari and Ananda (1977) which involved integration of the equation of motion and averaging of osculating orbital elements over one orbital period and removal of third body and solar radiation perturbations. The resulting orbital elements variations were related to lunar gravity features as explained in Ferrari and Ananda (1977) and is given by the following equation: ð15Þ d 1 ¼ A1 p  r1 ;

 _ _ X; x_ is the four-vector Kepler element where d 1 ¼ e_ ; I; rates defined at each instant, the parameter vector p is the 289-vector of harmonic coefficients in the gravitational

333

potential model, the sensitivity matrix A1 comprises of partial derivatives ðod 1 =opÞ as given by Kaula (1966), and r1 is a vector consisting of four Kepler rate errors. The second model used was that due to Wong et al. (1971) which computed the vertical components of acceleration at the nodes of a 5°  5° grid on a spherical surface at 100-km altitude using observation data of Apollo 8 and 12 at low altitude, low inclination, and in nearly circular orbit; and Lunar Orbiter 1, 2, 3, 4, and 5 observations which were in highly eccentric orbit. The effect of acceleration and the parameters involved in this model are given by the following normal equation: d 2 ¼ A2 p  r2 ;

ð16Þ

where r2 is the vector of data errors, d 2 is the 1008 – vector of acceleration, and A2 is the 1008  289 matrix of partial derivatives and in explicit form it can be written as  r ðnþ2Þ od 2 ðr; h; /Þ GM ¼ 2 ðn þ 1Þ Knm ðh; /Þ: ð17Þ A2 ¼ oGnm R R Here Knm ðh; /Þ is a spherical harmonic function of degree n and order m, normalized so that Z 2p Z þ1 K2nm ðl; /Þ dl d/ ¼ 4p; l ¼ sinðhÞ; ð18Þ 0

1

and Gnm is a two-vector of Fourier–Legendre coefficients. If the total lunar mass is correctly estimated and the origin of the coordinate system coincides with the lunar center of mass, then the harmonic coefficients of degrees zero and one vanishes identically. However, they included these terms in the processing of data using 600 discrete masses (Wong et al., 1971) which had a nonzero total mass and predominantly located on the lunar near side in order to remove the biases. Otherwise, the biases appear in the estimates of other harmonic coefficients due to the fact that when integrated over the sphere with variable weighting, the spherical harmonics no longer remain orthogonal. Thus the terms of degrees zero and one were purely a product of modeling procedure used by Wong et al. and the estimates obtained for these terms have no real physical significance (Bills and Ferrari, 1980). The third model that is due to Ferrari et al. (refer to Bills and Ferrari (1980)) which is sensitive to the low degree gravitational harmonics unlike the first two due to the inclusion of Doppler tracking data from high altitude phase of the Lunar Orbiter 4, and the laser ranging data to the lunar surface. The ranging data include the effect of lunar gravitational field through its interaction with lunar physical librations. The normal equations for this data set was given by Ferrari et al. (refer to Bills and Ferrari (1980)) d 3 ¼ A3 p  r3 ;

ð19Þ

where the matrix A3 contains the effect of the Doppler data of all gravity harmonics from degrees two through five plus the sectoral harmonics ðn ¼ mÞ of degrees six and seven. Thus the combined data from these different models provided high resolution coverage of the complete near

334

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

side as well as the near equatorial far side and weak resolution coverage of the remaining far side of the Moon. These data were constrained by augmenting them with a priori covariance model (Bills and Ferrari, 1980). Using two different a priori variance spectra, the first spectrum having an asymptotic form V ðG; nÞ ¼ Oðn3 Þ based on terrestrial gravity (Kaula, 1963) and the second one V ðG; nÞ ¼ Oðn4 Þ, constraining the same data set, they extracted two different solutions. The second spectrum functional form corresponds to the Bouguer correction term for topography with a Vening–Meinesz type spectrum (Bills and Ferrari, 1980). According to Bills and Ferrari even after combination of data from different models, few basic problems persisted because of localization of data over few regions especially the far side which still remained unsurveyed. This creates problems, in the event of higher degree and order coefficients were estimated, by rendering the eigenfunctions of the problem to be a linear combination of many different harmonic coefficients and because of use of least squares method wherein the squared residuals are weighted by the reciprocal of the data error variances. Therefore, the model produces large excursion on the far side regions with large data errors which is manifested as large power in the variance spectrum corresponding to the poorly determined coefficients. They reported that the field was fairly smooth and also showed reasonable correlation with the most significant mascons on the near-side. As a result of inclusion of various gravity parameterization schemes, which included observation data over a wide range of inclinations, in conjunction with indirect or direct data processing the consistency of the model obtained was much superior to the other models which preceded it. Besides, the use of a priori covariance, which reduces correlation among the estimated coefficients, led this model to produce a high accuracy near side gravity map without creating large amplitude variations on the far side. It was concluded by Bills and Ferrari that the produced model was more consistent with the Vening–Meinesz type spectral constraint. However, in the next section it is shown that the spectral power of this model matches with the Kaula’s type spectral constraint given as 1:5  104 =n2 . Thus extraction of this model using above-mentioned Kaula’s type constraint might have produced better results. This is a fact that up to degree and order of 8 the spectral power of the estimated coefficients with and without a priori constraint matches asserting the fact that coefficients are better estimated in this range. The use of Vening–Meinesz type spectral constraint resulted in excess power over the spectral constraint after the degree and order of 5 (Bills and Ferrari, 1980) which were corroborated to the effect of the constraint itself. However, the spectral power of this model as stated above coincides with the unconstrained model up to degree and order of 8–9, therefore, effect of constraint from degree and order 6 and above alone does not seem to be reasonable. Rather constraint is more effective after degree and order of 8–9. Thus it appears that the excess spectral power over the constraint shown by this

model from degree and order of 6 is perhaps combined effect of uncertainty in the estimation of the coefficients as a result of less amount of data used as compared to what is used today besides some effect from the constraint. In spite of all these differences this model was the best available one till 80s. Sagitov et al. (1986), in 1986, published results of a 16th degree and order spherical harmonic model. The model is based on five different data sources: Akim and Vlasova (1977, 1983), Ferrari (4th degree and order) (Ferrari, 1972) and Bills–Ferrari (16th degree and order) (Bills and Ferrari, 1980) models as well as a mascon model developed from observed LOS accelerations. In addition, near-side observations of Apollo 14, 15, 16 and 17 were used. Sagitov 16  16 model behaves very similarly to Bills and Ferrari (1980) (16th degree and order). 3. Recent efforts in spherical harmonic modeling A decade of lull followed the research in the lunar gravity modeling after the efforts made in the early 80s. In the late 80s, the President Bush’s proposal of a Space Exploration Initiative (in 1989) sparked a new interest in lunar mission planning. Renewed interest in lunar gravity modeling bloomed along with the great advancement in the computing technology that led to the development of the 60th degree and order lunar gravity model “Lun60d” by Konopliv (1993) by processing the observations for Lunar Orbiters (1–5) and Apollo 15 and 16 subsatellite missions. This was a first step towards higher degree spherical harmonic modeling which ultimately led to a step wise 165  165 spherical harmonic solution LP165P and a 150  150 one step spherical harmonic solution LP150Q. In the following subsection, a comparative picture has been presented about the developments starting from Lun60D spherical harmonic solution. 3.1. Lun60D solution based on Lunar Orbiters 1–5 and Apollo subsatellites data Lun60D (Konopliv et al., 2001) was the first high degree and order spherical harmonic solution for the lunar gravity after Bills and Ferrari (1980), and was made possible, firstly, because of availability of high computational power which facilitated large size matrix inversion, and secondly, the use of Kaula’s rule ð1:5  104 =n2 Þ which forces the higher degree and order spherical harmonics to be grouped together. The ability of this lunar gravitational model to predict the known lunar surface features which were in agreement with the line of sight results of Muller and Sjogren (1968), and also ability to predict the general far side features for which direct observations were not available, besides accurate orbit propagation generated much confidence as compared to the other models developed in seventies and 80s. The observations from Lunar Orbiters 1–5 and Apollo 15 and 16 subsatellites were processed using Double Precision Orbit Determination Program (DPODP)

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

developed by JPL which uses Square Root Information (SRIF) weighted least squares estimator (Ellis, 1980) in the coordinate system defined by the Earth’s Mean Equator at the epoch of J2000.0 (EME2000) and Developmental ephemeris (DE200) (Konopliv, 1993). The orientation of the Moon was obtained from Davies et al. (1992) in the J2000.0 inertial coordinate system. The Apollo 16 subsatellite crashed into the Moon’s surface 35 days after it was released from the Central Service Module (CSM). Konopliv (1993) compared the periapsis altitude of Apollo 16 subsatellite using the “Lun60d” gravity model truncated at degree and order 16 (the 60th degree and order (Lun60d)), Bills and Ferrari (1980) and Liu and Laing (1971) models with the actual observed altitude. They found that for Liu–Laing model, which did not use Apollo 16 data in gravity modeling, impact was not observed even after 60 days. The epoch value of the orbit was generated using first day of Doppler data from Apollo 16 using Lun60d gravity model. It should be noted that Lun60d uses Apollo 16 data as one day arcs. Bills–Ferrari model included Apollo data as mean element variations over one arc spanning 35 days and the simulation showed impact in 35 days. The Lun60d predicted impact in 34 days approximately which was found quite close to the actual value, and the truncated Lun60d (16  16) predicted the impact in 37 days. This is expected from fact that the variance spectrum of Lun60d contains higher power which is suitable for orbit propagation purpose. As reported in Konopliv (1993), the map of the lunar gravity anomalies appeared smooth when mapped on a surface at 100 km altitude, however, when the gravitational anomalies were mapped onto the lunar surface disparity with the reality became apparent rendering this model unsatisfactory for selenophysical studies. Konopliv et al. showed that the LOS acceleration obtained using this model provided the same information as the previous independent local estimates (Muller and Sjogren, 1968; Sjogren et al., 1974a,b,c; Sjogren and Smith, 1976) with coefficients of degree and order 60, besides the gravity contours followed the surface features much closer than previous harmonic solutions. In spite of this fact, the work done by Konopliv (1993) proved to be a landmark in the history of lunar gravity modeling, was best till that date, which provided impetus for further improvement. The efforts to improve the gravity model got boost with the availability of the data from Clementine and Lunar Prospector as described later, and resulted in a number of gravity models one after another of increasing degree and order. In 1994, the problem of lunar parameters estimation (moment of inertia) using Lunar Laser Ranging (LLR) data were again published (Dickey et al., 1994). LLR improved the lunar ephemeris by three-orders-of-magnitude. This also provided a number of useful informations for astronomical purposes. Besides added information on Earth’s precession, the Moon’s tidal acceleration, and lunar rotational dissipation.

335

3.2. Gravity models from Clementine and previous missions – GLGM-1 and 2 Lemoine et al. (1994, 1997) released their models Goddard Lunar Gravity Model-1 (GLGM-1) (70  70) and GLGM-2 (70  70), with a priori law of the form 1:5  104 =n2 , in 1994 and 1997, respectively. However, this field differs from Lun60d in that it used data from Clementine high altitude polar orbit of 430  2950 km, besides, data from previous lunar missions. The high altitude orbit of Clementine did not help in resolving the high frequency information, however, it provided constraint on the lower degree and order coefficients and improved it which were already estimated satisfactorily using data of previous missions only. This is evident from the fact the lower degree coefficients matched within some uncertainty from different models and their estimation without a priori constraint did not pose any problem. Moreover, the power spectrum of these unconstrained values was coinciding almost up to degree and order of eight to the solution with a priori constraints. Lemoine et al. reported that this field (actually evaluated at the lunar surface) showed 10% more anomalies when evaluated at an altitude of 100 km as compared to Lun60d and resolved some high latitude far side surface structures, which were not present in Lun60d, and correlated with the major mascons on the near side. The presence of low altitude data from previous missions made this model capture the short wavelength information and revealed, as reported by Lemoine et al., many mascons ringed by negative anomalies outside the mare basins suggesting the flexure of the lunar lithosphere. The details of tracking data used for GLGM-2 model can be referred to Lemoine et al. (1997). GLGM-2 is complete to degree and order of 70 which corresponds to a half wavelength resolution of 80 km. Lemoine et al. (1997) reported that this model has free-air gravity anomalies in a range 294 to +358 mGal. This was relatively better than that due to Lun60d which has a range from 454 to 529 mGal, this was due to the fact that power in the high degree terms of GLGM-2 were constrained due to a stiff Kaula’s constraint and down weighting of data in the calibration process. This implies that GLGM-2 under estimates the power in the lunar gravity field at short wavelengths. Besides, they also reported that this model was able to resolve major near side mascons and the far side basins such as Hertzsprung, Mendeleev, Tsiolkovsky, Moscoviense, and Freundlich–Sharonov which are gravity low consistent with the earlier gravity results (Ananda, 1977; Ferrari, 1977) (however, this result has been contradicted by the LP100J, LP100K and LP165P results (Konopliv et al., 2001) to be elaborated later). Although, data from Lunar Orbiters and Apollo subsatellites helped in capturing high frequency information but Clementine data alone captured the major mascon basins as reported by them (Lemoine et al., 1997). This is reflected in the fact that the many of the lows surrounding the positive highs

336

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

which were present in previous models were also present in the Clementine only free-air anomaly. In Figs. 1 and 2, the RMS magnitude spectra for GLGM-2, and Bills–Ferrari models are shown. Fig. 1 is plotted for full GLGM-2 model (up to degree 70), while Fig. 2 is plotted for GLGM-2 truncated to degree 16. The a priori law used for Bills–Ferrari model is given in the variance form as 1:4  106 =ð2n þ 1Þ2 ðnÞðn þ 1Þ. It is evident from Fig. 2 that GLGM-2 and Bills–Ferrari RMS magnitude spectra almost coincide except at degree

RMS Magnitude

1e−04

glgm2 glgm2−sigma a priori power law 1.5e−4/n^2 bills−ferrari−16 bills−ferrari−sigma

1e−05

1e−06

1e−07

1e−08

1e−09

0

10

20

30

40

50

60

70

Degree Fig. 1. Comparison of coefficient degree and coefficient sigma degree variances for GLGM-2 and Bills–Ferrari 16  16 spherical harmonic solutions.

1e−04

glgm2−16 glgm2−sigma a priori power law 1.5e−4/n^2 bills−ferrari−16 bills−ferrari−sigma

RMS Magnitude

1e−05

1e−06

1e−07

1e−08

0

2

4

6

8

10

12

14

16

Degree Fig. 2. Comparison of coefficient degree and sigma degree variances for GLGM-2 (truncated to degree 16) and Bills–Ferrari 16  16 solution.

5, 7, and 15 where Bills–Ferrari RMS model is having higher power. GLGM-2 is more reliable than Bills–Ferrari model because of large amount of data used. The lower power in GLGM-2 is a result of down weighting of data in the calibration procedure as mentioned earlier. Fig. 2 reflects that Bills–Ferrari model is close to the truncated GLGM-2 but a close look of the unnormalized coefficients in Table 6 shows offsets in the coefficients values. The coefficient C 31 has been underestimated as compared to other models while S 31 has been over estimated, the combined result is that RMS spectrum at degree 3 does not differ much from that of GLGM-2. Nonetheless, Bills–Ferrari model predicted the lunar features on the near side which matched with the features extracted using Lun60d model, and also showed smooth far side gravity field. But some of the negative anomaly region such as Freundlich–Sharanov which is located around 210° east and 15–18° north in the Bills and Ferrari (1980) gravity anomaly map has turned out to be a mascon surrounded by negative region as has been shown in Konopliv et al. (2001). Not only Bills–Ferrari model but other models up to GLGM-2 (Lemoine et al., 1997) also exhibit this region to be a negative one. This may be due to the higher altitude orbits, lesser data and lesser spherical harmonic coefficients, as compared to LP165P (Konopliv et al., 2001), estimated which render the anomaly negative in the stated region. Still GLGM-2 is considered to be an improvement over its predecessors as this model includes a lot of data from Clementine high altitude orbits and shows improvement in the lower order coefficients. The sigma variances of GLGM-2 and Bills–Ferrari models differ due to different data weighting schemes and larger uncertainty in the harmonic coefficients of Bills–Ferrari. Table 6 shows the values of low degree and order spherical harmonic coefficients for Bills–Ferrari, Lun60d, LLR and GLGM-2. The Lun60d and GLGM-2 coefficients are comparable while those due to LLR and Bills–Ferrari differ a lot. The major differences between the GLGM-1 and GLGM-2 are the weighting methods applied in the data processing and the use of Apollo subsatellite 16 data which was not used for GLGM-1. As stated in the previous paragraph, the high altitude of the Clementine-1 spacecraft does not provide any advantages for high degree gravity coefficients estimation, however, provided a lot of improvement

Table 6 Comparison of low degree and order spherical harmonics (unnormalized) from satellite tracking and lunar laser ranging (value  106). Coefficients

Bills–Ferrari

LLR

Lun60d

GLGM-2

J2 J3 C 22 C 31 S 31 C 32 S 32 C 33 S 33

202.431 ± 1.140 8.889 ± 1.508 22.263 ± 0.129 23.719 ± 1.123 7.16 ± 1.285 4.829 ± 0.048 1.626 ± 0.037 2.213 ± 0.143 0.341 ± 0.143

204.0 ± 1.0 8.66 ± 0.16 22.50 ± 0.11 32.4 ± 2.4 4.67 ± 0.73 4.869 ± 0.025 1.696 ± 0.009 1.73 ± 0.050 0.28 ± 0.020

203.805 ± 0.057 8.25 ± 0.06 22.372 ± 0.011 28.618 ± 0.018 5.87 ± 0.02 4.891 ± 0.009 1.646 ± 0.008 1.719 ± 0.003 0.211 ± 0.003

203.986 ± 0.131 9.99 ± 0.24 22.227 ± 0.023 28.290 ± 0.098 5.54 ± 0.07 4.886 ± 0.030 1.656 ± 0.029 1.756 ± 0.007 0.270 ± 0.006

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

337

Table 7 Comparison of low degree and order spherical harmonics (unnormalized) from satellite tracking (value  106). Coefficient

LP75G

LP100J

LP100K

LP150Q

LP 165P

J2 J3 C 22 C 31 S 31 C 32 S 32 C 33 S 33

203.428 ± 0.018 8.423 ± 0.215 22.348 ± 0.003 28.521 ± 0.008 5.869 ± 0.008 4.869 ± 0.003 1.665 ± 0.003 1.717 ± 0.001 0.246 ± 0.001

203.216 ± 0.012 8.407 ± 0.016 22.351 ± 0.002 28.461 ± 0.007 5.903 ± 0.005 4.850 ± 0.002 1.669 ± 0.002 1.713 ± 0.001 0.246 ± 0.001

203.199 ± 0.012 8.411 ± 0.016 22.350 ± 0.002 28.460 ± 0.007 5.905 ± 0.005 4.849 ± 0.002 1.668 ± 0.002 1.713 ± 0.001 0.246 ± 0.001

203.261 ± 0.012 8.474 ± 0.016 22.358 ± 0.002 28.452 ± 0.008 5.901 ± 0.005 4.845 ± 0.002 1.669 ± 0.202 1.713 ± 0.001 0.248 ± 0.001

203.236 ± 0.012 8.476 ± 0.017 22.357 ± 0.002 28.436 ± 0.008 5.902 ± 0.005 4.846 ± 0.002 1.671 ± 0.002 1.713 ± 0.001 0.249 ± 0.001

in the lower degree coefficients. This model (GLGM-2) failed to exhaust the signal in the tracking data though it performed better than its predecessors particularly with respect to the selenophysical structures (Lemoine et al., 1997). This is expected because of less power in this model at higher degree and order which makes it more suitable for selenophysical purpose than for orbit propagation. Lemoine et al. (1997) suggested alternate gravity analysis techniques such as the mapping of LOS accelerations for extracting further information. As stated earlier that this model intended to serve the purpose of selenophysical studies, but the difference in the far side features especially the negative anomaly which turned out to be positive one (Table 2 in Konopliv et al. (2001)) makes the far side gravity unreliable. In spite of this anomaly, the near side features were predicted correctly posing confidence in the model.

resolution became available which contained a lot of informations about the higher degree and order spherical harmonic coefficients. This resulted in two spherical harmonic solutions LP75D and LP75G. Moreover, the extended missions of Lunar Prospector which were in very low orbits 30km provided a lot of inputs about the high frequency signature of the lunar gravity and as a consequence of this, higher degree models such as LP100J, LP100K, one step solution LP150Q, and multi-step solution LP165P (Konopliv et al., 2001) became available. In this section gravity solutions based on tracking data from Lunar Orbiters, Apollo subsatellites, Clementine and Lunar Prospector are described. The LP75D and LP75G (Konopliv et al., 1998) models were published in 1998. The lower degree coefficients for the above models are summarized in Table 7. A summary of data used in various models is given in Table 8. A plot of RMS magnitude for GLGM-2 and LP75G is presented in Fig. 3 while that for LP75D and LP75G is presented in Fig. 4. In Fig. 3, a priori law ð1:5  104 =n2 Þ used for the GLGM-2 model and the exhibited spectral power of

3.3. Models based on Lunar Prospector data With the launch of Lunar Prospector satellite, for the first time tracking data for low lunar polar orbit with high

Table 8 Summary of tracking data in gravity models (Konopliv et al., 2001; http://pds-geosciences.wustl.edu/missions/lunarp/shadr.html). Mission

Gravity fields

LO-1 LO-2 LO-3 LO-4 LO-5 Apollo 15 subsatellite Apollo 16 subsatellite Clementine LP 100-km nominal mission January 11, 1998–February 17, 1998 LP 100-km nominal mission January 11, 1998–April 12, 1998 LP 100-km nominal mission January 11, 1998–December 19, 1998

Used in Used in Used in Used in Used in Used in Used in Used in LP75D

LP 40-km and 30-km extended mission December 19, 1998–February 8, 1999 LP 40-km and 30-km extended mission December 19, 1998–July 30, 1999

all all all all all all all all

LP75G LP100J LP100K LP165P LP150Q LP100J24 LP150Q LP165P LP100K LP165P LP150Q

Number of arcs

Arc length (days)

Number of observations

58 82 48 5 41 21 7 29 23

1 1 1 3 1 3 4 3 2

37,651 69,827 56,472 9,309 39,752 45,438 25,475 97,055 250,520

36

2

604,997

176

2

2,282,094

2

111

306,909

2

1,366,759

338

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349 1e−04

glgm2 glgm2−sigma a priori power law 1.5e−4/n^2 lp75g lp75g−sigma power law 2.5e−4/n^2

RMS Magnitude

1e−05

1e−06

1e−07

1e−08

1e−09

0

10

20

30

40

50

60

70

80

Degree Fig. 3. Comparison of coefficient degree and coefficient sigma degree variances for GLGM-2 and LP75G.

1e−04

lp75d lp75d−sigma lp75g lp75g−sigma power law 2.5e−4/n^2

RMS Magnitude

1e−05

1e−06

1e−07

1e−08

1e−09

0

10

20

30

40

50

60

70

80

Degree Fig. 4. Comparison of coefficient degree and coefficient sigma degree variances for LP75D and LP75G.

all the Lunar Prospector models, i.e. 2:5  104 =n2 is also plotted. The a priori law for the GLGM-2 and other models are given in Table 9. The RMS magnitude spectrum of GLGM-2 and LP75G both almost coincide up to degree 11. Thereafter, GLGM-2 shows lower power while LP75G follows the 2:5  104 =n2 power law. Aliasing is visible for the LP75G from 70° to 75°. LP75G differs from LP75D in that it contains some more tracking data of Lunar Prospector as shown in Table 8. LP75D shows little less power than the power law 2:5  104 =n2 after degree 40, whereas LP75G follows 2:5  104 =n2 (power content is 1:2  104 =n1:8 Konopliv et al., 1998). It has been reported by Konopliv et al. that the power spectrum for LP50pnoap (no a priori model LP gravity model up to degree and order 50) shows Table 9 A priori power law used in various gravity models (value  104) (in case of LP165P a priori law of the form ð2–5  104 =n2 Þ was used for individual elements around 122°, 145° and 165°).

n = 2–25 n = 26–30 n ¼ 31  nmax

GLGM-2

LP75G

LP100J, LP100K, LP150Q, LP165P

1.5 1.5 1.5

2.5 2.0 1.5

3.6 3.6 3.6

unexpected power beyond degree 16–20 (Konopliv et al., 2001; Konopliv et al., 1998) (data not available for LP50pnoap). It should be noted that earlier gravity models show unexpected power beyond degree 6–8 without a priori constraint. Therefore, LP50pnoap stability up to degree and order of 16 may be due to the large amount of data used both from the high and low attitude orbits besides better coverage as compared to any of the previous models. The signal to noise ratio for the LP75G model is higher than one implying under conditioning of the spherical harmonic solution, for higher degree, through a relaxed a priori spectral constraint and is true that the constraint for the LP75G was relaxed as compared to the GLGM-2 (Table 9). However, LP75G performs well in terms of orbit propagation, besides, revealing high frequency features in orbit perturbation. The LP75G model even though it shows higher spectral power still provides a lot of informations about the selenophysical structures. It is reported in Konopliv et al. (1998) that processing of the Lunar Prospector data revealed seven new large mascons with maximum field at the center of the basin surrounded by comparatively weak field in the surrounding of the basin. The three were located are located on the near side of the moon beneath the impact basins Mare Humboltianum, Mendel–Rydberg, and Schiller–Zucchius, where the latter basin has no visible mare fill. The far side four mascons are in the large far side basins of Hertzsprung, Coulomb–Sarton, Freundlich–Sharonov, and Mare–Moscoviense (these features appeared to be negative anomalies in GLGM-2 and previous models). It was concluded by Konopliv et al. (1998) due to lack of far side data the far side mascons showed large uncertainties. They also concluded that in spite of all these shortcomings better global coverage obtained through Lunar Prospector improved the coherence which is a measure of correlation of gravity with the topography for the mid-wavelength frequencies ðn ¼ 20–50Þ, while the admittance of topography and the gravity was also more apparent than the previous models (Konopliv et al., 1998). Another most important result from LP75G is the normalized polar moment of inertia of the Moon. The solution from Lunar Laser Ranging (LLR) data are C=MR2 ¼ 0:394  0:002, (Dickey et al., 1994; Konopliv et al., 1998) and while gravity solution from earlier lunar spacecrafts gives C=MR2 ¼ 0:393  0:001 (Konopliv et al., 1998; Williams et al., 1996). With the availability of the Lunar Prospector gravity models the new solution turns out to be C=MR2 ¼ 0:3932  0:0002 and the average moment I=MR2 ¼ 0:3931  0:0002 (Konopliv et al., 1998). The uncertainty in the polar moment of inertia and the average moment is five times the formal error, which enters through the use of J 2 and C 22 coefficients having uncertainty of five times the formal error. However, the resulting uncertainty in C=MR2 is reduced by about a factor of 5 over previous estimates. As a whole it appears that LP75G is better than GLGM-2 not only for the orbit propagation purpose but for selenophysical studies too. The follow on Lunar Prospector gravity models (LP100J, LP100K, and LP165P) were published in Konop-

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349 1e−04

1e−04 lp75g1 lp100k lp100k−sigma power law 2.5e−4/n^2 lp75g1−sigma

lp100k lp165p power law 2.5*e−4/n^2 lp165p−sigma lp100k−sigma

1e−05

RMS Magnitude

RMS Magnitude

1e−05

1e−06

1e−07

1e−08

1e−09

339

1e−06

1e−07

1e−08

0

10

20

30

40

50

60

70

80

90

1e−09

100

0

20

40

60

Degree

80

100

120

140

160

180

Degree

Fig. 5. Comparison of coefficient degree and coefficient sigma degree variances for LP75G and LP100K.

Fig. 7. Comparison of coefficient degree and coefficient sigma degree variances for LP100K and LP165P.

liv et al. (2001). The data source for these Lunar Prospector gravity models are shown in Table 8. In contrast with LP75G model which was constrained using the power constraint of the form 1:5  104 =n2 , a relaxed power constraint was used for LP100J 3:6  104 =n2 (Konopliv et al., 2001; http://pds-geosciences.wustl.edu/missions/lunarp/shadr.html; Konopliv and Yuan, 1999). Data tracked from Pomonkey tracking station was not used in this analysis because of tracking file conversion problem. The improved results obtained by LP100J model are reflected in the errors for the equatorial near side about 20 mGal and as large as 100+ mGal for far side short wavelength mascon features as reported by Konopliv et al. (2001). LP100J was followed by LP100K gravity model which contained the rest of the extended mission data of low altitude orbits (Konopliv et al., 2001) (Table 8). The LP100J and LP100K models are probably the best for the mission planning purpose as they are not compute intensive as other models which followed it. RMS magnitude plots for Lunar Prospector models LP75G, LP100J, LP100K, and LP165P are shown in Figs. 5–7. The LP100K model shows little more spectral power as compared to LP75G and aliasing is seen in both the models. LP100K uncertainty approaches the signal asymptotically, whereas for LP75G it remains

less than the signal. These are results of the data weighting and a priori law used. The a priori law 3:6  106 =n2 used for the LP100K which is relaxed as compared to that for the LP75G of 1:5–2:5  104 =n2 makes the LP100K model more realistic. The LP100J and LP100K models are almost similar because of the same processing method adopted for both of them. Both are suitable for the orbit propagation purpose but the latter one contains more data from Lunar Prospector extended mission, therefore, this becomes choice for the orbit prediction purpose. In Fig. 7, a comparison between LP100K and LP165P is shown. LP165P which contains data from extended LP mission provides a smooth solution without aliasing up to degree 110 (versus degree 90 for the 100th degree models) and resolves numerous craters (Konopliv et al., 2001). However, LP165P has been superseded by one step spherical harmonic solution LP150Q (http://pds-geosciences.wustl.edu/missions/lunarp/shadr.html). A multistep solution (Konopliv and Yuan, 1999; Konopliv et al., 1999) was utilized to obtain the LP165P gravity solution with spectral constraint (Konopliv et al., 2001) that resulted in strong aliasing at the end and beginning of the cutoff degrees 122 and 145 as these two cutoffs were the steps involved in estimation process (Fig. 7), besides high degree 2.5e−07

1e−04

lp100j lp100k lp100k−sigma power law 2.5e−4/n^2 lp100j−sigma

2e−07

RMS Magnitude

RMS Magnitude

1e−05

1e−06

1e−07

1.5e−07

1e−07

5e−08

1e−08

1e−09

rmsdiff−100j−165p sigma−variance−100j sigma−variance−165p

0

10

20

30

40

50

60

70

80

90

100

Degree Fig. 6. Comparison of coefficient degree and coefficient sigma degree variances for LP100J and LP100K.

0

0

10

20

30

40

50

60

70

80

90

100

Degree Fig. 8. Comparison of RMS difference between coefficients of LP100J and LP165P and sigma-variance of LP100J and LP165P.

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

Table 10 RMS difference between coefficients of LP100J and LP165P; and sigmavariance of LP100J (value  109). Degree

RMS difference of coefficients

Sigma variance LP100J

Sigma-variance LP165P

2 3 4 5

8.66543739 15.99033479 18.12240589 18.88296197

4.12964752 5.69389585 7.47661400 10.30862341

4.64975771 6.19025912 7.96803196 10.84707080

1e−04

lp150q lp100k lp150q−sigma lp100k−sigma power law 2.5e−4/n^2

1e−05

RMS Magnitude

coefficients were more correlated and had to be suppressed. As compared to LP100J or LP100K, the LP165P fits to the Doppler residuals much better. The radial accuracy for LP100J or LP100K for the extended mission is about 2 m; while perpendicular to it is 20 m which is a lot more than that for LP165P in the radial direction for the nominal mission (0.5 m) (Carranza et al., 1999). Thus LP165P is better than LP100J or LP100K in terms of orbit propagation. The method used for the gravity modeling of LP165 is due to Ellis (1980). A spectral constraint of the form 3:6  104 =n2 was used to restrict the power content of the gravity field at higher degree. The observed power of LP165P is the same as for the other Lunar Prospector models, i.e. 2:5  104 (Konopliv et al., 2001). The normalized polar moment of inertia obtained by this model is consistent with the previous ones (LP100J and LP100K, C=MR2 ¼ 0:3932  0:0002) and the lunar core constraints reported for the LP75G model (Konopliv et al., 1998) implying the accuracy or the consistency of these models. It was reported by Konopliv et al. (2001) that the value of this parameter is sensitive to the relative data weight between the Lunar Prospector data and the Clementine data which is perhaps due to the lack of the far side observations. However, it seems that normalized polar moment of inertia which depends on the J 2 and C 22 coefficients play important role. The estimates of these parameters get influenced by the high altitude orbits of Clementine in the weighting process. Another contribution from LP165P gravity map may be existence of mascons in six large far side basins (Konopliv et al., 2001) (four of these were also identified as mascons by LP75G model) which were earlier identified as negative anomalies (Ananda, 1977; Konopliv, 1993; Lemoine et al., 1997). The RMS of differences between the coefficients of the LP100J and LP165P models are larger than the sigma-variance (or formal error) of LP100J by about a factor of 2 for coefficients less than degree 20 while for coefficients larger than 20 are quite small comparatively as shown in Fig. 8 and Table 10. It is evident from Fig. 8 that uncertainty in the estimation of the coefficients for LP100J and LP165P is lower than the RMS differences of coefficients or in other words the RMS of differences of coefficients of the two models is about twice the formal error. Uncertainty levels for both the models are same (plots for sigma-variance overlap). It has been stated earlier for the Bills–Ferrari and GLGM-2 spherical harmonic solutions that even

1e−06

1e−07

1e−08

1e−09

0

20

40

60

80

100

120

140

160

Degree Fig. 9. Comparison of coefficient degree and coefficient sigma degree variances for LP100K and LP150Q.

though their spectral variance values matched at certain degrees but the coefficients themselves differed in the two models (either C nm over/under estimated and S nm under/ over estimated with respect to each other). For the case of LP100J and LP165P the same is true. Due to the use of more data in the case of LP165P from the extended mission as compared to LP100J the estimates of coefficients differ in these two models. Moreover, degree of the model being estimated also influences these coefficients though these are more evident as aliasing effect for the last few degrees of the model. Since the RMS of differences of the coefficients is about twice the formal error, therefore, both the models are consistent. However, more data in the LP165P model gives more confidence in it as compared to LP100J. It is described later that using LP100K model, which is similar to LP100J, Doppler residual of 5 mm/s for near side while 33 mm/s for the far side was obtained in orbit fit for SELENE SST data. Thus even though the current models seem to fit well to the near side data and predict some of the selenophysical features, still these are away from its reality. LP150Q is a one step solution up to degree and order of 150 and therefore, the best till date, but is compute intensive for normal orbit propagation purposes. The LP100K 1e−04

lp150q lp165p power law 2.5e−4/n^2 lp165p−sigma lp150q−sigma

1e−05

RMS Magnitude

340

1e−06

1e−07

1e−08

1e−09

0

20

40

60

80

100

120

140

160

180

Degree Fig. 10. Comparison of coefficient degree and coefficient sigma degree variances for LP150Q and LP165P.

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349 2.5e−07

3.5e−07

rmsdiff−150q−165p sigma−variance−150q sigma−variance−165p

rmsdiff−100k−150q sigma−variance−100k sigma−variance−150q

3e−07

RMS Magnitude

2e−07

RMS Magnitude

341

1.5e−07

1e−07

2.5e−07 2e−07 1.5e−07 1e−07

5e−08 5e−08 0

0

20

40

60

80

100

120

140

0

160

0

10

20

30

40

Degree Fig. 11. Comparison of RMS difference between coefficients of LP150Q and LP165P and sigma-variance of LP150Q and LP165P.

alternatively can be used for orbit propagation purpose thereby avoiding large computation while compromising little on the accuracy of orbit propagation. Figs. 9 and 10 show the power spectrum for LP100K, LP165P, and LP150Q. In Fig. 9, a strong aliasing can be seen for both LP100K and LP150Q gravity models. The LP150Q model shows less power in its RMS spectrum as compared to LP165P. Moreover, the RMS of difference of coefficients of LP150Q and LP165P is less than the sigma variance in any of the two models (Fig. 11). RMS difference exceeds the sigma variance only after degree 110 up to which the RMS spectrum of LP165P is smooth. In the light of above facts it seems that LP150Q may be more reliable than the LP165P model, however, further investigations are required using the orbit modeling or selenophysical interpretation in conjunction with topography study to resolve it. Fig. 12, shows the RMS difference of LP100J and LP150Q. The RMS difference of coefficients is about twice larger than the formal error indicating the consistency of the result. In Fig. 13, the RMS difference of LP150Q and LP100K is shown. Fig. 14 shows the RMS difference of LP100J and LP100K models. It is obvious from these fig-

60

70

80

90

100

Fig. 13. Comparison of RMS difference between coefficients of LP100K and LP150Q and sigma-variance of LP100K and LP150Q.

ures that LP100J and LP100K differ slightly due to some more data used in LP100K even though they are of the same order and the same data processing and weighting was used. Thus it is evident that amount of data changes the spectral power of the model even though their sigma variances remain same. RMS difference of LP100K and LP150Q is of the similar nature as that for LP100J and LP150Q. Therefore, it may conclude that amount of data, degree of solution sought, type of constraint, and the weight used decide the quality of the coefficients extracted. Alone looking as the spectral power of the coefficients may be misleading, any conclusion must be justified in conjunction with the RMS difference magnitude plot, residuals from the not used tracking data etc. In Table 11, GM values from different models have been listed. Clearly, the uncertainty in the LP150Q model is the least and is an indication of consistency of this model. 3.4. Some other methods in lunar gravity modeling and future possibilities The role and methods of satellite altimetry in gravity field modeling has been discussed in Andersen and Knudsen (1995), Andersen et al. (1996) and Sandwell and Smith

3.5e−07

2.2e−07

rmsdiff−100j−150q sigma−variance−100j sigma−variance−150q

3e−07

2e−07

rmsdiff−100j−100k sigma−variance−100j sigma−variance−100k

1.8e−07 1.6e−07

RMS Magnitude

2.5e−07

RMS Magnitude

50

Degree

2e−07 1.5e−07 1e−07

1.4e−07 1.2e−07 1e−07 8e−08 6e−08 4e−08

5e−08

2e−08 0

0

10

20

30

40

50

60

70

80

90

100

Degree Fig. 12. Comparison of RMS difference between coefficients of LP100J and LP150Q and sigma-variance of LP100J and LP150Q.

0

0

10

20

30

40

50

60

70

80

90

100

Degree Fig. 14. Comparison of RMS difference between coefficients of LP100J and LP100K and sigma-variance of LP100J and LP100K.

342

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

Table 11 Comparison of GM for Moon from different solutions. Analysis

Description

GM (km3/s2)

Ferrari et al. (1980)

LLR and LO-4 Lun60d (60  60) GLGM-1 (70  70) GLGM-2 (70  70) LP75G (75  75) LP100K (100  100) LP150Q (150  150) LP165P (165  165)

4902.7993 ± 0.0029

Konopliv (1993) Lemoine et al. (1994, 1997) Lemoine et al. (1997) Konopliv et al. (1998) http://pds-geosciences.wustl.edu/ missions/lunarp/shadr.html http://pds-geosciences.wustl.edu/ missions/lunarp/shadr.html Konopliv et al. (2001)

4902.79781 ± 0.00122 4902.80263 ± 0.0095 4902.80295 ± 0.00224 4902.80027 ± 0.00023 4902.800238 ± 0.000206 4902.80107 ± 0.00008 4902.801056 ± 0.000217

(1997). The same model has been used by Andersen and Knudsen for modeling gravity field in coastal areas (Andersen and Knudsen, 2000) but it determines the local features not the global ones. Further, use of satellite altimetry for gravity modeling of Mars to degree and order of 80 (Goddard Mars Model 2B (GMM-2B)) using X band tracking data of Mars Global Surveyor (MGS) and altimetry crossovers formed from the Mars Orbiter Laser Altimeter (MOLA) has been reported in Lemoine et al. (2001). On the local scale gravity modeling, a new method has been proposed in Goossens et al. (2005a,b) for the recovery of regional lunar gravity fields. A short-arc approach has been adopted which consists of the time a satellite takes to cross the grid of interest on the lunar surface, and is used to filter out most long-wavelength signal that can be present in the residuals. Simulation showed that contamination of data either with the systematic error or stochastic noise did not much affect the result and local gravity field was recovered down to the level of several mGal (Goossens et al., 2005a). This method consists of a linear variational approach where in the tracking data residuals are mapped to the gravity anomalies. Based on this method a high resolution gravity anomalies on the lunar surface was recovered using the Lunar Prospector tracking data residuals (Goossens et al., 2005b). The results for Mare Serenitatis with model were found to be consistent and comparable with the global recovery results with much less computation. The possibility for better results for spherical harmonic coefficients exists in Satellite-to-Satellite Tracking (SST) and Satellite Gravity Gradiometry (SGG). The SST technique involves flying a formation in the low orbit where the motion of satellite gets affected by the higher degree and order coefficients which tend to get smeared at higher altitudes. A stable formation can capture the gravity variation of the far side because of the visibility for tracking the master by the slave satellite over a long period of time. Or at least one satellite should be in the low orbit which is

supposed to capture the variation in gravity (Floberghagen et al., 1996). Low-low configuration has its own advantages over the high-low configuration, in that it contains much richer information about high frequency information of the gravity field (Floberghagen et al., 1996; Milani et al., 1996; Matsumoto et al., 1999; Namiki et al., 1999; Sasaki et al., 2003). Details of this technique can be referred to Floberghagen et al. (1996) and Floberghagen (2002). The Satellite Gravity Gradiometry (SGG) can provide drastic improvement in the gravity model. The SGG involves measuring the gravity tensor, i.e. the second order derivatives of the gravitational potential. Till today no SGG study has been done for mapping the lunar gravity field because of accuracy required in the measurements and the complexity of the experiment together with the cost involved. In fact, the two best candidates for the accurate gravity modeling, especially for the Moon where the far side direct observations are not available, are SST and SGG. Work done in 80s on SGG are reported in Wagner (1983), Kaula (1983), Colombo (1984) and Rummel and Colombo (1985). The satellite gradiometer senses the gravitational field in several directions to measure the complete curvature structure of the local gravitational field or as the case may be. Two variants of SGG has been reported for the estimation of the spherical harmonic coefficients, the time-wise approach, and the space-wise (Ditmar et al., 2003). The time-wise approach has been elaborated in Colombo (1986) and Koop (1993). In the time-wise approach each measurement is considered as an element of a time series along the satellite orbit. Coordinates of observations are expressed in terms of orbital elements and in this case inclination functions are used to represent the spherical harmonics as given in Kaula (1966). The simplicity of representation of the time-wise approach and the analytical expressions that directly relate the gravity field and the SGG tensor (provided that the z-axis of Local Orbit Reference Frame (LORF) points strictly in the radial direction) makes it amenable to implementation. Moreover, colored noise in data can also be treated easily. But its rigorous implementation is time consuming. Variants of this methods are also available such as that based on Fourier transform but its application to large data set is perhaps fraught with difficulties and is too restrictive in the sense of various assumptions made (Schrama, 1991). In the space-wise approach coordinates of an observation is expressed using associated Legendre functions using geographical coordinates. The observation equation formed in terms of gravity tensor and gravity coefficients are essentially linear in gravity coefficients and therefore, all the methods discussed for the gravity modeling are applicable. Space-wise method is particularly suitable for the satellite based gravimetry. Various methods used to solve for this case have been discussed in Ditmar et al. (2003) and a new method has been proposed based on conjugate gradients with preconditioning and fast algorithms. However, a lunar mission with gradiometer is awaited before its result can be analysed. Various literatures mainly

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

tackle the issue of handling large amount of data arising out of SGG.

3.5. A brief of current estimation techniques in gravity modeling Primarily the gravity model can be extracted either in the form of LOS acceleration, discrete masses, or spherical harmonic coefficients. All the three have been elaborated in previous sections. Essentially, the solution technique for spherical harmonic coefficients involves processing of observation data to get the information array. The data processing for estimating the spherical harmonic coefficients can be broadly categorized as (a) those aiming to provide better regularization to constrain the solution in the absence of far side data, and (b) those aiming to solve the normal equation more efficiently and accurately. The first one is concerned with better modeling of the problem itself by providing realistic constraints to obtain better solution such as working with spectral constraints (Eq. (10)) or spatial constraint modeling (Konopliv et al., 1999). Besides, weighting of data from various sources can also be seen as the way of regularization. While the second one comprises of different least squares methods or their variants which aim to reduce the fitting errors and tackle the ill conditioning of the normal matrix. In fact, if the normal matrix is not ill condition its solution by different methods should be same within some computational uncertainty. However, ill conditioning may make the matrix poorly invertible then the way of solving becomes important and may result in solutions which may differ considerably. However, there is a overlap between these two approaches as the regularization also helps in overcoming the illposedness of the problem, and therefore, partially addresses the issues involved in the second approach. Therefore, an integrated approach to gravity modeling merges the two methods, i.e. first formulating the objective function constrained with various regularization schemes and followed by solution of the resulting normal equation. The different ways of regularization and solution of normal equation has been elaborated in Golub and Loan (2007). Besides, regularization has been elaborately discussed in the context of gravity modeling in Floberghagen (2002). The spectral constraint as a way of regularization has been discussed in the previous section. Here, spatial constraint as given in Konopliv and Yuan (1999) is briefly presented for the sake of completeness. This constraint is based on the evaluation of the radial acceleration ðan Þ from all coefficients of degree n and its uncertainty on the reference sphere, an is given by GM an ¼ 2 ðn þ 1Þ ae n X   P nm ðsin /Þ C nm cos mk þ S nm sin mk ;  m¼0

ð20Þ

343

where ae is the radius of the reference sphere (spatial constraining of the signal did not work for the lunar gravity case). The expected acceleration profile is then given by GM pffiffiffiffiffiffiffiffi ð21Þ ðan ÞRMS ¼ 2 K 2=n; ae and is derived from Kaula’s rule "Pn  #1=2 2 2 m¼0 C nm þ S nm ¼ K=n2 ; 2n þ 1

ð22Þ

which gets reduced to Eq. (21) for n 1. On the reference sphere at each point, Eq. (21) defines the expected “signal” for the acceleration (Konopliv and Yuan, 1999). Alternatively, for different regions RMS of the acceleration can be computed and this can be defined. The uncertainty for the coefficients of degree n is written as rðan Þ and is given by Eq. (23) rðan Þ ¼ rða2;n Þ  rða2;n1 Þ;

ð23Þ

where rða2;n Þ is given by Eq. (24) rða2;n Þ ¼

oaT2;n oa2;n P noapð2;nÞ ; e e 2;n o G 2;n oG

ð24Þ

and oa2;n GM ¼ 2 ðn þ 1ÞP nm ðsin /Þ cos mk; oC n;m ae oa2;n GM ¼ 2 ðn þ 1ÞP nm ðsin /Þ sin mk; oS n;m ae

ð25Þ ð26Þ

e 2;n is a vector containing all the normalized gravity where G coefficients of degree 2–n, and P noapð2;nÞ is the corresponding covariance. This method basically constrains the noise of the gravity field to 0 with some uncertainty when the noise exceeds the signal (Konopliv and Yuan, 1999). The acceleration due to all the spherical harmonic coefficients of degree greater than or equal to the degree strength (degree strength of gravity field is the crossing point of Kaula’s signal with the acceleration uncertainty for a particular latitude and longitude) is constrained to zero, at the surface, with an uncertainty approximately equal to the expected signal at the degree strength. The acceleration (observation) of degree strength D is given by aD;N ¼

N GM X ðn þ 1Þ a2e n¼D



n X

  P nm ðsin /Þ C nm cos mk þ S nm sin mk

m¼0

and therefore, the linearized observation equation is given by zi ¼ Ai~x þ mi ; oaD;N Ai ¼ e oG

344

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

where zi is the difference between the nominal and the observed value (the accumulated acceleration at the surface for degrees D to N from the a priori gravity, observed value is zero in this case), Ai is the row vector of observation partials, ~x is a vector containing the parameters to be estimated, and mi is the observation error. The observations thus generated are merged with the unconstrained gravity square root information array using Householder transformation as described in Ellis (1980). In the normal form then the constrained gravity estimate ~x is given by Eq. (27) h i1 h i T T ~x ¼ P 1 P 1 z ; ð27Þ noap þ A WA noap xnoap þ A W ~ where P noap is the unconstrained covariance of the gravity coefficients, A is observation matrix containing observation partials, W is the diagonal weight matrix, ~xnoap is the unconstrained gravity estimate, and ~z is a vector of linearized observation. It showed very high power for the higher frequencies greater than 15, and therefore, not good for selenophysical implications. Therefore, the spatial constraining technique has not been used for the lunar gravity modeling problem. The process of solving the normal equation can be carried out using various methods as described in Golub and Loan (2007). All the methods do not find their use in lunar gravity modeling. Some of the methods used for solving the normal equation which are widely used are based on Householder transform, Cholesky decomposition, singular value decomposition, ridge type estimation, least squares collocation etc. However, their utility depends on the number of parameters to be estimated, quality of data being handled and the computational power available including the size of Random Accessible Memory (RAM). The recent gravity models are based on Householder Transform based processing as elaborated in Ellis (1980) using Kaula’s spectral constraint. This is not elaborated here as this method is quite widely used and well appreciated. The ridge type estimation (Minter and Cicci, 1997) is briefly presented here. The ridge type estimation (Minter and Cicci, 1997) was used for the Mars gravity field, to optimally weight a set of a priori constraints based upon Kaula’s rule, in order to obtain more accurate estimates. This was based on Householder transformation based ridge type method which he published in Cicci (1992) and is elaborated in the following paragraph. The minimum variance estimate of the normal equation y ¼ Hx þ ;

ð28Þ

is given by the  1 ^x ¼ H T R1 H H T R1 y;

ð29Þ

where y is an ðm  1Þ vector of observation residuals, H is an ðm  nÞ mapping matrix, x is ðn  1Þ vector of corrections to nominal state, and  is an ðm  1Þ vector of observation errors. R is the covariance matrix of the observation errors. The covariance of the estimate ^x is given by

1

P ¼ ðH T R1 H Þ ;

ð30Þ

when one or multiple eigen values of the normal matrix ðH T R1 H Þ has very small eigen values, then the model variables are linearly dependent making the problem ill conditioned. The solution obtained for this type of problem may give small residuals but the estimate may be far from the reality. Therefore, such solutions are constrained by using a priori covariance information in the normal matrix formation. The minimum variance estimate with a priori information is given by  1 ^x ¼ H T R1 H þ P 1 H T R1 y; ð31Þ where P is the a priori covariance matrix and serves as a measure of the uncertainty in the nominal state. The addition of the a priori covariance matrix removes the ill conditioning of the normal matrix while making the magnitude of the eigenvalues high. If there is uncertainty in the value of matrix P , estimates may by highly erroneous. In Eq. (31) the a priori covariance term is equally weighted against the observational data. Therefore, an erroneous a priori covariance leads to an inaccurate estimate. To eliminate this problem, ridge type method is used. The estimate can be written as ^x ¼ D^xN ;

ð32Þ

where ^xN is given by  1 ^xN ¼ H TN H N þ DP 1 D DH T R1 y;

ð33Þ

and H TN H N ¼ DH T R1 HD:

ð34Þ

In Eq. (34), D is a square matrix of size n whose elements are equal to the inverse of the square root of the diagonal elements of ðH T R1 H Þ, and P is assumed to be diagonal and positive definite. In Eq. (33) the covariance term can be a diagonal matrix of size n of biasing parameters, K. This gives the ridge type solution as given by Eq. (35)  1 ^xN ¼ H TN H N þ KDP 1 D DH T R1 y: ð35Þ The value of K is found by minimizing the Mean Square Error function MSE ¼

n X i¼1

1

þ 2 ð1 þ k i Þ

2 n X ðk i xN i Þ i¼1

ð1 þ k i Þ

2

;

ð36Þ

where k i is ith diagonal element of the K matrix and xN i is ith element of the vector xN . The presence of small amount of bias in ill-conditioned problem serves to reduce the total variance of the solution so that the mean square error is reduced. Since MSE is a function of the true solution, an iterative solution must be used to approximate the value of each k i which minimizes the MSE function. The ridge type Householder transformation is used to solve Eq. (35) instead of matrix inversion or the Cholesky decomposition which can be referred to Minter and Cicci (1997) and Cicci (1992). It was concluded by Minter and Cicci (1997) that

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

use of ridge type estimation method to with proper weighting of Kaula’s constraint improved the accuracy of the Mars gravity field model Mars50c derived from Viking and Mariner 9 Doppler tracking data. They also concluded that ridge-type estimation method fine-tunes the solution and was obvious when statistical comparison was made with other solutions. However, improvement in higher degree terms was less visible due to aliasing and other effects in least squares method. The gravity model obtained by this method appeared to be comparable in accuracy to Mars50c in that total variance and Sum of the Squares of Residuals for Mars50c and ridge type were similar. However, weighting was not determined before the solution was obtained because the biasing parameters enable the a priori covariance weighting to be independent of the solution process. Besides, least squares collocation technique (Tscherning, 2001) also exists but for lunar gravity modeling this has not been reported. This provides good estimate if the a priori information is very good, in the absence of which estimate derived is not reliable. 3.6. Recent lunar missions and gravity modeling activities Two lunar satellites were launched in 2007, SELENE and Chang’e-1. Status of functionality of Chang’e-1 is not known at present, however, gravity modeling activities using SST data from SELENE is going on. The Chang’e-1 orbiter, the first of Chinese missions to the Moon, was launched on October 24, 2007. The orbiter has a mass of 2350 kg, containing around 1180 kg of propellant and scientific payloads of 130 kg. It was deployed in a 127 min period, 200 km altitude, circular high-inclination orbit. This mission planned lifetime was 1 year during which sent data about the lunar environment and surface regolith. But more details and any information about the gravity modeling with the data of this spacecraft is not available. Chandrayaan-1 has been launched on October 22, 2008 and many other missions are waiting. In the following sections a brief account is given for SELENE and Chandrayaan-1. 3.6.1. SELENE mission On September 14th, 2007, SELENE was launched which deployed two small spin stabilized subsatellites, Rstar and Vstar in October 2007 in the lunar orbit. The three satellites were tracked using RSAT, a SST Doppler tracking sub-system, and VRAD (artificial radio sources for VLBI). Multifrequency differential VLBI tracking between Rstar and Vstar for the far side gravity and the 4-way Doppler tracking between the main satellite and Rstar and were utilized. The Doppler residual obtained using LP100K lunar gravity model over the near side was 5 mm/s while the far side residuals were as large as 30 mm/s reflecting the fact that currently modeled gravitational fields are not representing the far side gravity properly (Sasaki et al., 2008). Drastic improvement in gravity model is expected using SELENE SST data. A simulation run with Lunar Prospector track-

345

ing data for SELENE is available in Goossens and Matsumoto (2007). It appears from the previous review that all the models have some deficiency to small or large extent irrespective of the gravity model or the estimation technique used due to the under sampling of data for the far side. Gravity coefficients obtained by satellite tracking are not very good as when they are reduced to the height of the earth surface show numerous artifacts. In spite of all the efforts to compensate for the under sampling of data using regularization methods, the uncertainties in the gravity coefficients remain significant and show considerable anomalies when reduced to the surface of the Moon. Moreover, from the previous paragraph it is apparent that even the Doppler residuals are high for the far side with LP100K solution indicating the shortcoming present in these models. The launch of SELENE has brought a lot of hope in improving the gravity solutions using the SST method. 3.6.2. Chandrayaan-1 The Chandrayaan-1 was launched on October 22, 2008 through PSLV C11 (Polar Satellite Launch Vehicle) from the Satish Dhawan Space Center Sriharikota (SHAR) and had a planned life time of 2 years. This is an experimental satellite to enhance the technological capabilities in space and giving some important scientific information about the lunar surface. The Lunar Laser Ranging Instrument (LLRI) was supposed to map the surface topography. However, at present due to failure of the onboard star sensor LLRI has been shut down following the altitude raising from 100 km to 200 km. Tracking is going on but due to altitude increase the sensitivity of the orbit for the higher degree harmonic coefficients is lost. However, it is hoped that data gathered before the sensor failure would help in gravity and topography modeling to some extent till the Chandrayaan-2 is launched. In the following paragraph, the reference system and the methods of gravity modeling to be used is briefly discussed. The recent lunar gravity fields (Lunar Prospector series) were developed using the JPL planetary ephemerides DE403. This ephemeris program calculates the orientation of the Moon with respect to EME2000 in terms of three Euler angles (a) a rotation measured from the x-axis of EME2000, /, about the z-axis of EME2000 to the ascending node of the lunar equator, (b) a rotation about the new x-axis generated due to the above rotation to coincide with the lunar equator by an angle h, and (c) the rotation by an angle w along the lunar equator (actually principal axes being defined). These angles were generated using numerical integration suggested in Newhall and Williams (1997). Besides, this also provides the ephemeris of different planets and the sun which is used for the third body perturbation computation. The gravity modeling for Chandrayaan1 has been formulated using ephemerides from DE405 which refers to International Celestial Reference System (ICRS) as the inertial reference frame (the latest model DE418 and DE421 are also available and update to these

346

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

models is being considered). The Earth rotation, precession and nutation are computed using the IAU2006/IAU2000A precession and nutation models (IAU2006 P03 precession model and IAU2000A nutation model) using the Celestial Intermediate Origin (CIO) based transformation. The orientation of the Moon is derived from DE405 ephemeris program. Since the DE405 is oriented on to the ICRS, therefore, the ICRS has been chosen for the complete simulation of observation vector of the lunar satellite. This is for the first time that these models are being implemented for the lunar gravity modeling problem. A brief account is given in the next paragraph. Usually, the reference frame chosen for the gravity modeling is J2000.0 (EME2000) as referred to in various literatures (Konopliv et al., 2001). But the current JPL planetary ephemeris DE405 is oriented on to the ICRS. From the literature, it is not apparent whether J2000.0 offsets with ICRS has been accounted for or not in the modeling. Therefore, it was felt that the latest developments in the Earth orientation and measurement should be utilized and use the ephemerides generated by JPL DE405 which is a dynamical realization of the ICRS. The Earth orientation which is given by precession-nutation that locates the changing direction of the Earth’s pole and an adopted origin of right ascension on the celestial sphere. The coordinate system for this celestial sphere is the ICRS and the two directions are the Celestial Intermediate Pole (CIP) and Celestial Intermediate Origin (CIO). Thus, the Earth’s pole orientation with respect to the ICRS can be described in terms of rotations which take care of (a) offset (bias) of Mean Equator of J2000.0 pole with respect to the ICRS, (b) precession, and (c) nutation. There are two equivalent bias-precession-nutation transformations from ICRS to International Terrestrial Reference System (ITRS), (a) CIO and (b) equinox based. The equinox can be defined geometrically on the surface of a sphere, while the CIO is defined kinematically using the concept of Non Rotating Origin (NRO) of Guinot (1979) and Capitaine et al. (1986). With more precise observations available, the IAU 2000 precession-nutation model was adopted by the IAU (Resolution B1.6) to replace the IAU 1976 precession model (Lieske et al., 1977) and the IAU 1980 theory of nutation (Wahr, 1981; Seidelmann, 1982) in order to better describe the orientation of the CIP. The IAU 2000 precession-nutation model consists of IAU 2000A and IAU 2000B formulation, former being more accurate uses MHB 2000 model provided in Mathews et al. (2002). The differences of the origin of the x-axis in the ICRS and the equinox of J2000.0 have been properly modelled in the IAU2000 model (in the IERS 2003 convention) which was not the case for the previous precession-nutation models. Based on the concept of the NRO a combined representation of the precession-nutation in the form of X and Y, coordinates of Celestial Intermediate Pole (CIP), was suggested (Capitaine, 2000) which is consistent with the IAU 2000A precession-nutation model. With the availability of the accurate measurements of the sub-daily polar motion,

the nutation and the polar motion have been separated out in the modeling in the high frequency domain. Therefore, all the terms with period less than 2 days are considered as the polar motion. The IAU2006 (Resolution B1) has replaced the IAU 2000A precession model with P03 precession theory (Capitaine et al., 2003c) because it was not consistent with the dynamical theory (Hilton et al., 2006). Moreover, terms like lunisolar precession and planetary precession has been replaced by (Resolution B1) precession of the equator and precession of the ecliptic, respectively (Hilton et al., 2006). The procedure to implement the IAU 2006 based precession-nutation model has been published in Wallace and Capitaine (2006). Thus in this new framework, the observation vectors are simulated which is followed by gravity modeling. The estimation problems in the gravity modeling framework is being tackled using various methods such as Householder transformation (Ellis, 1980), Householder transformation based on ridge type method (Minter and Cicci, 1997), and the method proposed by Lerch (1991). The method proposed by Lerch has been used by Lemoine et al. (1997). This method uses a subset least-squares solutions of the satellite data contained within the complete solution and requires that the differences of the parameters of the subset solutions and the complete solution to be in agreement with their error estimates by adjusting the data weights. Lerch reported that estimated data weights are markedly smaller than the weights implied by the formal uncertainties of the measurements and the orbital test arcs as well as surface gravity comparisons showed significant improvements for solutions when more realistic data weighting is achieved. In fact, in the present scenario Doppler residual on the untrained data set from various orbits (high, low, polar, and equatorial) must be small if the model extracted really represents the global model. Doppler residuals on the trained orbit (used in gravity modeling) do not represent the universal representation. Therefore, it is required that not only data weight selection be improved (as it is one of the major factors influencing the model) but data clustering may also need further improvement as it is also an important factor in the extraction of a global gravity model. The latest models LP100J, LP100K, LP165P, and LP150Q by Konopliv et al. follow the method proposed by Ellis (1980). The main problem of the gravity modeling is the solution of the normal equation with less amount of computation, high accuracy, and consistency. Proper data weighting has to be taken into account for any such processing besides constraining the solution with Kaula’s type spectral constraint. However, constraining the solution makes the coefficient estimates biased and solution will depend on the constraint used. Therefore, a method based on data mining technique such as clustering can be first used to cluster the data of different types into different groups and then apply the fuzzy weighting to give weightage to different clusters thus made. Fuzzy sets in this case may be the inherent choice which can provide relation/

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

dependencies between data in qualitative sense or in terms of qualitative and quantitative measures. A good survey of fuzzy data mining can be referred to Mitra et al. (2002). A number of references are available on fuzzy modeling but reference (Ross, 2005) provides a good background in illustrative way. Currently, such formulation is under way. Moreover, neural network sensitivity analysis is also planned for lunar gravity modeling and has been successfully utilized for the aircraft parameter estimation. The key issue here will be extraction of spherical harmonic coefficients using the neural sensitivity analysis method as described in Garhwal et al. (2007). The beauty of this method (Garhwal et al., 2007) is that it does not require any a priori values of the aircraft parameters as other methods do and without which they may diverge and produce wrong results. The results obtained were very encouraging as it provided the same level of accuracy without a priori constraint as other methods did with a priori constraint (derived from wind tunnel testing) as described in Jategaonkar (2006). Further work on this method has shown unexpectedly good results. Neural sensitivity analysis works well for the linear problems, however, its real potential can be realized for the nonlinear problems hence the whole problem of coefficients determination has to be recasted in the form of a nonlinear equation rather than in the form of normal equation. For example, the LOS gravity map can be extracted in a much better way under the same assumption as followed in Muller and Sjogren (1968). Though current status of Chandrayaan-1 has given the setback but it is hoped that all the proposed methods of gravity modeling will get materialized in short time with the launch of Chandrayaan-2.

347

though continuous improvement in the gravity modeling has taken place but absence of the observation data of the far side of the Moon creates anomalies in the recovered model in the sense that some of them perform good at sampling altitude while fail at the lunar surface though the latest model LP150Q shows quite improved results. While others, which perform well at the surface level show skewed results at the higher altitudes. It has been shown using the RMS difference of the coefficients that amount of tracking data, weighting method, and degree of the solution sought directly affects the quality of the solution. Doppler residuals only for the near side are not enough to judge the performance of the model which became evident from the SELENE satellite to satellite tracking data residuals. Therefore, it is concluded that it is mandatory to have the far side observation data if we want to have a global gravity model. Besides, we need to cross check the validity of results through Doppler residual, surface features recognition, RMS difference of the coefficient determined on subset and full set of data. Thus satellite to satellite tracking in the low lunar orbit which will have the much richer signature of the lunar gravity field or gradiometry is indispensable for the modeling purpose in the case of the Moon. A brief review of Chandrayaan-1 has also been presented which uses the ICRS reference frame for all the inertial calculations. Acknowledgements This study is supported by ISRO, India, for Chandrayaan-1 mission. Authors are grateful to anonymous reviewers whose comments helped in enormously improving this manuscript.

4. Conclusion References In this paper, a review of development of gravity modeling during late 60s to early 80s and thereafter since nineties has been presented. These modeling techniques have been broadly categorized based on the way gravity has been parameterized, i.e. LOS modeling, discrete mass modeling and spherical harmonic modeling. The method of orbit propagation which is essentially the foundation work in gravity modeling has been divided as direct and indirect method and have been elaborated while discussing different parameterization schemes. Direct method integrates the equation of motion directly while the indirect method uses averaging technique to average the equation of motion which is then integrated to propagate the averaged orbital elements. Besides, hybrid formulation has also been mentioned which tries to capture the direct and indirect methods using data from these developments to form the normal matrix and thereby, solving for the spherical harmonic coefficients. The gravity modeling activity since nineties has also been described in details which have been categorized under spherical harmonic modeling using direct method because direct integration of equation of motion has been used. From the review it is obvious that

Akim, E.L. Determination of the gravity field of the Moon by the motion of Luna. Doklady Akademii Nauk USSR, 170, 1966. Akim, E.L., Vlasova, Z.P. Model of the lunar gravitational field, derived from Luna 10, 12, 14, 15 and 22 tracking data. DAN SSSR 235, 38–41, 1977. Akim, E.L., Vlasova, Z.P. Studies of the lunar gravitational field involving Doppler tracking of Soviet artificial satellites. Kosmicheskaya Issledovaniya 21, 499–511 (in Russian), 1983. Ananda, M.P. Mean rates of the orbital elements of a satellite perturbed by a lens shaped mass concentration. Celestial Mechanics 12, 495–511, 1975b. Ananda, M.P. Lunar gravity: a mass point model. Journal of Geophysical Research 82, 3049–3064, 1977. Andersen, O.B., Knudsen, P. Global gravity field from the ERS-1 geodetic mission. Earth Observation Quarterly 47, 1–5, ESRIN, 1995. Andersen, O.B., Knudsen, P. The role of satellite altimetry in gravity field modeling in coastal areas. Physics and Chemistry of the Earth (A) 25 (1), 17–24, 2000. Andersen, O.B., Knudsen, P., Tscherning, C.C. Investigation of methods for global gravity field recovery from dense ERS-1 geodetic mission altimetry. In: Rapp, C., Nerem, I.A.G. (Eds.), Global Gravity Field and Its Temporal Variations, Symposium No. 116, Springer, Berlin, pp. 218–226, 1996. Bierman, G.J. A comparison of discrete linear filtering algorithms. IEEE Transaction on Aerospace Electronic Systems 9, 28–37, 1973.

348

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349

Bills, B.G., Ferrari, A.J. A harmonic analysis of lunar gravity. Journal of Geophysical Research 85 (B2), 1013–1025, 1980. Blackshear, W.T., Gapcynski, J.P. An improved value of the lunar moment of inertia. Journal of Geophysical Research 82, 1699–1701, 1977. Bogoliubov, N.N., Mitropolosky, Y.A. Asymptotic Methods in the Theory of Nonlinear Oscillations (Translated from Russian). Hindustan Publishing Corp., New Delhi, 1961. Bryant, W.C., Williamson, R.G. Lunar gravity analysis results from Explorer 49, AIAA Paper 74-810. In: Presented at the AIAA Mechanics and Control of Flight Conference, American Institute of Aeronautics and Astronautics, Anaheim, California, August 1974. Campbell, M.J., O’Leary, B.T., Sagan, C. Moon: two new mascon basin. Science 164, 1273–1275, 1969. Capitaine, N. Definition of the celestial ephemeris pole and the celestial ephemeris origin. In: Johnston, K.J., McCarthy, D.D., Lazum, B.J., Kaplan, G.H. (Eds.), Towards Models and Constants for Submicroarcsecond Astrometry, US Naval Observatory, pp. 153–163, 2000. Capitaine, N., Guinot, B., Souchay, J. A non-rotating origin on the instantaneous equator: definition, properties and use, Celestial Mechanics 39, 283–307, 1986. Capitaine, N., Wallace, P.T., Chapront, J. Expressions for IAU2000 precession quantities. Astronomy and Astrophysics 412 (2), 567–586, 2003c. Carranza, E., Konopliv, A., Ryne, M. Lunar Prospector orbit determination uncertainties using high resolution lunar gravity models. In: AAS/AIAA Astrodynamics Specialist Conference, Girwood, AK, August, AAS paper 99-325, 1999. Cicci, D.A. Application of the householder transformation to ridge-type estimation methods. Applied Mathematics and Computation 51 (2 and 3), 159–165, 1992. Colombo, O.L. The global mapping of gravity with two satellites. Netherlands Geodetic Commission, New Series, 7, 3, Delft, 1984. Colombo, C.L. Notes on the mapping of gravity field using satellite data, in: Sun¨kel, H. (Ed.), Mathematical and Numerical Techniques in Physical Geodesy, Lecture Notes in Earth Sciences, vol. 7. Springer, Berlin, pp. 261–316, 1986. Davies, M.E., Abalakin, A., Brahic, A., Bursa, M., Chovitz, B.H., Lieske, J.H., Seidelmann, P.K., Sinclair, A.T., Tjuflin, Y.S. Report of the IAU/IAG/CSPAR working group on cartographic coordinates and rotational elements of the planets and satellites: 1991. Celestial Mechanics 53, 377–397, 1992. Dickey, J.O., Bender, P.L., Faller, J.E., Newhall, X.X., Ricklefs, R.L., Ries, J.G., Shelus, P.J., Veillet, C., Whipple, A.L., Wiant, J.R., Williams, J.G., Yoder, C.F. Lunar laser ranging: a continuing legacy of the Apollo program. Science 265, 482–490, 1994. Ditmar, P., Klees, R., Kostenko, F. Fast and accurate computation of spherical harmonic coefficients from satellite gravity gradiometry data. Journal of Geodesy 76, 690–705, 2003. Ellis, J. Large scale state estimation algorithms for DSN tracking station location determination. Journal of Astronautical Sciences 28, 15–30, 1980. Ferrari, A.J. An empirically derived lunar gravity field. Moon 5, 390–410, 1972. Ferrari, A.J. Lunar gravity derived from long period satellite motion – a proposed method. Celestial Mechanics 7, 46–47, 1973. Ferrari, A.J. Lunar gravity: the first far side map. Science 188, 1297–1300, 1975. Ferrari, A.J. Luhnar gravity: a harmonic analysis. Journal of Geophysical Research 82 (20), 3065–3084, 1977. Ferrari, A.J., Ananda, M.P. Lunar gravity: long-term Keplerian rate method. Journal of Geophysical Research 82, 3085–3097, 1977. Ferrari, A.J., Heffron, W.G. Effects of physical librations of Moon on the orbital elements of a lunar satellite. Celestial Mechanics 8, 111–120, 1973. Ferrari, A.J., Sinclair, W.S., Sjogren, W.L., Williams, J.G., Yoder, C.F. Geophysical parameters of the Earth–Moon system. Journal of Geophysical Research 85 (B7), 3939–3951, 1980.

Floberghagen, R. Lunar Gravimetry: Revealing the Far-Side. Kluwer Academic Publishers, The Netherlands, pp. 1–2, 2002. Floberghagen, R., Noomen, R., Visser, P.N.A.M., Racca, G.D. Global lunar gravity recovery from satellite-to-satellite tracking. Planetary Space Science 44 (10), 1081–1097, 1996. Gapcynski, J.P., Blackshear, W.T., Tolson, R.H., Compton, R. A determination of the lunar moment of inertia. Geophysical Research Letter 2 (8), 353–356, 1975. Garhwal, R., Halder, A., Sinha, M. Aircraft parameter estimation using neural sensitivity analysis. In: IEEE Conference, Kuala Lumpur, Malaysia, November 2007. Golub, G.H., Loan, C.F.V. Matrix Computations, third ed Hindustan Book Agency, 2007. Goossens, S., Matsumoto, K. Lunar satellite orbit determination analysis and quality assessment from Lunar Prospector tracking data and SELENE simulations. Advances in Space Research 40, 43–50, 2007. Goossens, S., Visser, P.N.A.M., Ambrosius, B.A.C. A method to determine regional lunar gravity fields from Earth-based satellite tracking data. Planetary and space Science 53, 1331–1340, 2005a. Goossens, S., Visser, P.N.A.M., Heki, K., Ambrosius, B.A.C. Local gravity from Lunar Prospector data: results for Mare Serenitatis. Earth Planets Space 57, 1127–1132, 2005b. Gottlieb, P. Estimation of local gravity features. Radio Science 5, 303– 312, 1970. Guinot, B. Basic problems in the kinematics of the rotation of the Earth. In: McCarthy, D.D., Pilkington, J.D. (Eds.), Proceedings of IAU symposium, Time and the Earth’s Rotation, D. Reidel, Dordrecht, pp. 7–18, 1979. Hilton, J.L., Capitaine, N., Chapront, J., Ferrandiz, J.M., Fienga, A., Fukushima, T., Getino, J., Mathews, P., Simon, J.-L., Soffel, M., Vondrak, J., Wallace, P., Williams, J. Progress report of the IAU Working Group on precession and the ecliptic. Celestial Mechanics and Dynamical Astronomy 94 (3), 351, 2006. Available from: . Jategaonkar, R.V. Flight Vehicle System Identification: A Time Domain Methodology. AIAA Inc., Reston, VA, 2006. Kaula, W.M. The investigation of the gravitational fields of the Moon and planets with artificial satellites. Advances in Space Science Technology 5, 210–230, 1963. Kaula, W.M. Theory of Satellite Geodesy. Blaisdell Publishing Company, USA, 1966. Kaula, W.M. The gravitational field of the Moon. Science 166, 1591–1598, 1969. Kaula, W.M. Inferences satellitesatellite of variations in the gravity field from satellite-to-satellite range rate. Journal of Geophysical Research 88 (B10), 8345–8349, 1983. Koch, K.R., Morrison, A. A simple layer model of the geopotential from a combination of satellite and gravity data. Journal of Geophysical Research 75 (8), 1483–1492, 1970. Konopliv, A.S., Sjogren, W.L., Wimberly, R.N., Vijayaraghavan, A. A high resolution lunar gravity field and predicted orbit behavior. In: AAS/AIAA Astrodynamics Specialist Conference, American Astronautical Society, Victoria, BC, Canada, 16–19 August 1993, AAS paper 93-622, 1993. Konopliv, A.S., Yuan, D.N. Lunar Prospector 100th degree gravity model development. In: 30th Lunar and Planetary Science Conference, Houston, TX, USA, 15 March 1999. Konopliv, A.S., Binder, A.B., Hood, L.L., Kucinskas, A.B., Sjogren, W.L., Williams, J.G. Improved gravity field of the Moon from Lunar Prospector. Science 281, 1476–1480, 1998. Konopliv, A.S., Banerdt, W.B., Sjogren, W.L. Venus gravity: 180th degree and order model. Icarus 139, 3–18, 1999. Konopliv, A.S., Asmar, S.W., Carranza, E., Sjogren, W.L., Yuan, D.N. Recent gravity models as a result of the lunar prospector mission. Icarus 150 (1), 1–18, 2001.

M. Sinha et al. / Advances in Space Research 45 (2010) 322–349 Koop, R. Global gravity field modeling using satellite gravity gradiometry, Publ Geodesy, New series, No. 38, Nedersland Geodetic Commission, Delft, 1993. Lawson, C.L., Hanson, R.J. Solving Least Square Problems. PrenticeHall, Englewood Cliffs, NJ, pp. 222–231, 1974. Lawson, C.L., Hanson, R.J. Solving least squares problems. In: SIAM Classics in Applied Mathematics, vol. 15, Society for industrial and Applied Mathematics, 1995. Lemoine, F.G., Smith, D.E., Zuber, M.T. Goddard lunar gravity model-1 (GLGM-1): a 70th degree and order gravity model for the Moon (abstract), Eos Trsns. AGU, Fall Meeting Supplementary, 75(44), 400, 1994. Lemoine, F.G., Smith, D.E., Zuber, M.T., Neumann, G.A., Rowland, D.D. A 70th degree lunar gravity model (GLGM-2) from Clementine and other tracking data. Journal of Geophysical Research 102, 16339– 16359, 1997. Lemoine, F.G., Smith, D.E., Rowlands, D.D., Zuber, M.T., Neumann, G.A., Chinn, D.S., Palvis, D.E. An improved solution of the gravity field of Mars (GMM-2B) from Mars Global Surveyor. Journal of Geophysical Research 106 (E10), 23359–23376, 2001. Lerch, F.J. Optimum data weighting and error calibration for the estimation of gravitational parameters. Bulletin Geodesique 65, 44– 52, 1991. Lieske, J.H., Lederle, T., Fricke, W., Morando, B. Expression for the precession quantities based upon the IAU (1976) system for astronomical constants. Astronomy and Astrophysics 58, 1–16, 1977. Liu, A., Laing, P. Lunar gravity analysis from long-term effects. Science 173, 1017–1020, 1971. Lorell, J. Lunar orbiter gravity analysis. Moon 1, 190–231, 1970. Lorell, J., Sjogren, W.L. Lunar gravity: preliminary estimates from lunar orbit. Science 159, 625–627, 1968. Lorrel, J., Anderson, J.D., Lass, H. Application of the method of averages to celestial mechanics. Technical Report 32-482, JPL, Pasadena, California, 1964. Mathews, P.M., Herring, T.A., Buffett, B.A. Modeling of nutationprecession: new nutation series for nonrigid Earth, and insights into the Earth’s interior. Journal of Geophysical Research 107, B4, doi:10.1029/2001JB00390, 2002. Matsumoto, K., Heki, K., Rowlands, D. Impact of far side satellite tracking on gravity estimation in the SELENE project. Advances in Space Research 23 (11), 1809–1812, 1999. Michael, W., Blackshear, W. Recent results on the mass, gravity field, and moments of inertia of the Moon. Moon 3, 388–402, 1972. Michael Jr., W.H., Blackshear, W.T., Gapcynski, J.P. Dynamics of satellite. In: Morando, Bruno (Ed.) Proceedings of the Prague 12th Plenary of COSPAR and 10th International Space Science Symposium, May 11–24, 1969, pp. 42–56, 1969. Michael Jr., W.H., Blackshear, W.T., Gapcynski, J.P., in: Morando, Bruno (Ed.), Results on the Mass and Gravitational Field of the Moon as Determined from Dynamics of Lunar Satellites (1969). Springer, Berlin, pp. 2–56, 1970. Milani, A., Luise, M., Scortecci, F. The lunar subsatellite experiment of the ESA MORO mission: goals and performances. Planet and Space Sciences 44 (10), 1065–1076, 1996. Minter, C.F., Cicci, D.A. Improvement in mars gravity field determination using ridge-type estimation methods. In: ASS 97-702, pp. 1549–1566, 1997. Mitra, S., Pal, S.K., Mitra, P. Data mining in soft computing framework: a survey. IEEE Transactions on Neural Networks 13 (1), 2002. Muller, P.M., Sjogren, W.L. Mascons: Lunar mass concentrations. Science 161, 680–684, 1968. Muller, P.M., Sjogren, W.L., Wollenhaupt, W.R. Lunar gravity: Apollo 15 radio tracking. Moon 10, 195–205, 1974. Namiki, N., Hanada, H., Tsubokawa, T., Kawano, N., Ooe, M., Heki, K., Iwata, T., Ogawa, M., Takano, T., and RSAT/VRAD/LALT mission

349

groups. Selenodetic experiments of SELENE: relay satellite, differential VLBI, and laser altimeter, Advances in Space Research, 23(11), 1817–1820, 1999. Newhall, X.X., Williams, J.G. Estimation of the lunar physical libration. Celestial Mechanics and Dynamical Astronomy 66, 21–30, 1997. Ross, T.J. Fuzzy logic: with engineering applications. Wiley (Asia), 2 Clementi Loop #02-01, Singapore, 2005. Roy, A. The Foundations of Astrodynamics. The Macmillan Company, Collier-Macmillan Limited, London, 1965. Rummel, R., Colombo, O.L. Gravity field determination from satellite gradiometry. Bulletin Geodesique 59, 233–246, 1985. Sagitov, M.U., Bodri, B., Nazarenko, V.S., Tadzhidinov, Kh.G. Lunar Gravimetry, 35th of International Geophysics Series. Academic Press, New York, 1986. Sandwell, D.T., Smith, W.H.F. Marine gravity anomaly from Geosat and ERS-1 satellite altimetry. Journal of Geophysical Research 102, 10039–10054, 1997. Sasaki, S., Iijiama, Y., Tanaka, K., Kato, M., Hashimoto, M., Mizutani, H., Takizawa, Y. The SELENE mission: goals and status. Advances in Space Research 31 (11), 2335–2340, 2003. Sasaki, S., Namiki, N., Handa, H., Iwata, T., Kawano, N., Matsumoto, K. New lunar gravity filed from SELENE (KAGUYA) gravity experiment using two subsatellites. Geophysical Research Abstracts, 10, EUG2008-A-05821, 2008. Schrama, E. Gravity field error analysis: applications of GPS receivers and gradiometry on low orbiting platforms. Journal of Geophysical Research B96 (B12), 20041–20051, 1991. Seidelmann, P.K. 1980 IAU nutation: the final report of the IAU working group on nutation. Celestial Mechanics 27, 79–106, 1982. Sjogren, W.L. Lunar gravity estimate: independent confirmation. Journal of Geophysical Research 76, 7021–7026, 1971. Sjogren, W.L., Muller, P.M. Wollenhaupt, Apollo 15 gravity analysis from S band transponder experiment. Moon 4, 411–418, 1972b. Sjogren, W.L., Smith, J.C. Quantitative mass distribution models for Mare Orientale. In: Proceedings of Lunar Science Conference, 7th, Houston, Texas, pp. 2639–2648, 1976. Sjogren, W.L., Gottlieb, P., Muller, P.M., Wollenhaupt, W. Lunar gravity via Apollo 14 Doppler radio tracking. Science 175, 165–168, 1972a. Sjogren, W.L., Wimberly, R.N., Wollenhaupt, W.R. Lunar gravity via the Apollo 15 and 16 subsatellites. Moon 9, 115–128, 1974a. Sjogren, W.L., Wimberly, R.N., Wollenhaupt, W.R. Lunar gravity: Apollo 16. Moon 11, 35–40, 1974b. Sjogren, W.L., Wimberly, R.N., Wollenhaupt, W.R. Lunar gravity: Apollo 17. Moon 11, 41–52, 1974c. Tscherning, C.C. Computation of spherical harmonic coefficients and their error estimates using least-squares collocation. Journal of Geodesy 75, 12–18, 2001. Tuckness, D.G., Jost, B. A critical analysis of the numerical and analytical methods used in the construction of the lunar gravity potential. Applied Mathematics and Computation 71, 69–90, 1995. Wagner, C.A. Direct determination of gravitational harmonics from lowlow GRAVSAT data. Journal of Geophysical Research 88 (B12), 10309–10321, 1983. Wahr, J.M. The forced nutation of an elliptical, rotating, elastic and oceanless Earth. Geophysics Journal of Royal Astronomical Society 64, 705–727, 1981. Wallace, P.T., Capitaine, N. Precession-nutation procedures consistent with IAU 2006 resolutions. Astronomy & Astrophysics 459, 981–985, 2006. Williams, J.G., Newhall, X.X., Dickey, J.O. Lunar moments, tides, orientation, and coordinate frames. Planetary and Space Science 44 (10), 1077–1080, 1996. Wong, L., Buechler, G., Doens, W., Sjogren, W., Muller, P., Gottlieb, P. A surface layer representation of the lunar gravity field. Journal of Geophysical Research 76, 6220–6236, 1971.