- Email: [email protected]
A method to determine regional lunar gravity fields from earth-based satellite tracking data
ARTICLE IN PRESS
Planetary and Space Science 53 (2005) 1331–1340 www.elsevier.com/locate/pss
A method to determine regional lunar gravity fields from earth-based satellite tracking data S. Goossens,1, P.N.A.M. Visser, B.A.C. Ambrosius Department of Earth Observation and Space Systems, section Astrodynamics and Satellite Systems, Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The Netherlands Received 17 March 2005; received in revised form 23 June 2005; accepted 29 June 2005 Available online 18 August 2005
Abstract A method to determine regional gravity fields of the Moon from Earth-based Doppler and range satellite tracking data residuals of a low Moon-orbiting satellite has been developed and thoroughly tested in a controlled simulation environment. A short-arc approach, where one arc consists of the time it takes the satellite to cross the grid of interest on the lunar surface, is used in order to filter out most long-wavelength signal that can still be present in the residuals. Simulation results where the data are contaminated with either typical systematic or stochastic noise show that recovery of the local gravity field down to the level of several mGal is possible. The inclusion of extremely low-altitude data also means that regularisation in the sense of including a priori information in the form of a regularisation matrix is not necessary in order to obtain a good solution at high resolution. r 2005 Elsevier Ltd. All rights reserved. Keywords: Regional gravity field determination; Lunar gravity field; Gravity anomalies
1. Introduction The Moon and the Earth are in a spin-orbit resonance, which causes the Moon to always have the same face towards the Earth. This severely limits Earth-based tracking of a satellite orbiting the Moon, since no tracking data are available when the satellite flies over the far side of the Moon. This in its turn hampers the determination of the global gravity field of the Moon in terms of spherical harmonics, since less than one half of the Moon is lacking data coverage. By using the mathematical tool of regularisation, the ill-posed inverse problem can be transformed into a well-posed one, and a solution can be obtained. For the Moon, such solutions have been determined in the past, e.g. Lemoine et al. (1997) and Konopliv et al. (2001) using relatively Corresponding author.
E-mail addresses: [email protected], [email protected] (S. Goossens). 1 Currently at National Astronomical Observatory of Japan, Mizusawa. 0032-0633/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.pss.2005.06.009
recent data from the Clementine and Lunar Prospector spacecraft, respectively, and these solutions perform very well in terms of data fit and correlation with topographical features. In order to determine these solutions, a priori information in the form of a Kaula rule is used, which is based on the observation that the spectrum of the spherical harmonic coefficients of the gravity field behave according to b=l 2 , with l the degree and b a constant, see Kaula (1966). However, it has been shown by Floberghagen (2002) that the choice of the regularisation parameter or which Kaula constant to apply, is of great influence on the final result. When using global basis functions such as spherical harmonics, it cannot be avoided that the solution for the near side is influenced by that for the far side, and vice versa. Using global basis functions also limits the resolution that can be obtained. However, the current data set with excellent near side coverage with Lunar Prospector data suggests more information is present in the data than currently can be modelled using spherical harmonics alone (Konopliv et al., 2001). Due to the low
ARTICLE IN PRESS 1332
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
altitude of Lunar Prospector, at an average of 30 km in its Extended Mission, a resolution of 1 on the lunar surface should be achievable. Current spherical harmonics expansions are limited to degree 165 maximum, and this solution has been estimated in multiple steps (Konopliv et al., 2001). Thus, in order to exhaust the information present in the data and to overcome problems inherent to the use of global basis functions, a different representation needs to be used. Such a representation is expected to achieve a higher resolution for areas on the near side, which will aid the interpretation of selenophysical features. Gravity is represented in this paper by means of gravity anomalies on the lunar surface. A method to determine these anomalies from range and Doppler tracking data residuals has been developed and thoroughly tested. In previous efforts the data type used to determine gravity anomalies has often been line-of-sight (LOS) accelerations, which can be obtained by differentiating Doppler residuals with respect to time. Examples for the Moon include Muller and Sjogren (1968), who revealed the discovery of the mascons on the Moon, and Sjogren et al. (1972, 1974a–c) and Muller et al. (1974), who used LOS data due to coverage limitations and achieved a further improvement of knowledge of such features as Mare Serenitatis and Mare Crisium. Recent efforts for the Moon using LOS data of Lunar Prospector can be found in Sugano and Heki (2004a,b). In this paper however, the recovery will be done directly from range and Doppler data residuals. This is done in order to avoid problems with noise on the data which will be amplified in the differentiation process and to avoid problems coming from a possible lack of continuous data series. Similar problems were also noted by Barriot et al. (1997), who then also decided to avoid the delicate matters of numerical differentiation. This paper deals with the recovery method itself, and by means of extensive simulations in a controlled environment it is shown what level of recovery can be achieved. In a follow-on paper, real Lunar Prospector data will be dealt with. The simulations in this paper thus keep Lunar Prospector and its characteristics in mind. The paper is structured as follows: the recovery method itself will be explained in Section 2. In Section 3, simulation results with different error sources on the data will show the achievable recovery level. Section 4 contains the conclusions.
2. Recovery method The observation model that relates the tracking data residuals to gravity anomalies is formed by a linear variational approach. The approach used is an extension of the approach used by Vonbun et al. (1980). Let r_ be
the initial Doppler measurement. A residual Dr_ is obtained by comparing a computed measurement, created from reference models, with the actual observation. Subsequently, adjustments to the reference models can be made until a certain desired level of for instance data fit is achieved, and the final residuals are obtained. In this way, most effects on such a residual can be filtered, and the final residual itself can be thought to be caused by local gravity anomalies D~ g or a change in the initial state vector, given by D~ x0 . In a linear variational form, by using a first order Taylor expansion, the dependency of the data residuals on changes in the initial state vector and any kind of dynamical force model parameters ~ p becomes: qr_ q~ x qr_ q~ x Dr_ ¼ D~ p, (1) D~ x0 þ q~ x q~ x0 q~ x q~ p where ~ x is the state vector itself and ~ x0 is the initial state _ x are geometrical partials and vector. The partials qr=q~ can easily be formed from geometrical relationships by expressing the range-rate in terms of the coordinates of the satellite and ground station. The partials q~ x=q~ x0 and q~ x=q~ p are sensitivity coefficients. These coefficients can be obtained by integration of the variational equations, which are derived by differentiation of the equation of motion, see e.g. Montenbruck and Gill (2000). Let the equation of motion be expressed as d~ xðtÞ=dt ¼ f~ðt; ~ xðtÞ; ~ pÞ, with ~ xðtÞ the full 6 parameter state vector consisting of positions and velocities, and t the time. The function f~ describes the forces acting on the satellite, and it is dependent on the time, the state vector itself, and the force parameters. By differentiation of the equation of motion with respect to the parameters ~ p, the differential equation for q~ xðtÞ=q~ p reads: d q~ xðtÞ qf~ðt; ~ xðtÞ; ~ pÞ q~ xðtÞ qf~ðt; ~ xðtÞ; ~ pÞ ¼ þ . dt q~ p q~ xðtÞ q~ p q~ p
(2)
The difference between the formulation as used in this paper with the formulation used by Vonbun et al. (1980) is an extra term in the differential equation for the sensitivity coefficient q~ x=q~ p. This extra term is the gradient term qf~=q~ x of the forces acting on the satellite, and this term does not appear in the formulation used by Vonbun et al. (1980). There, it is mentioned that ancillary data types (basically, ground-based tracking of a satellite since they use satellite-to-satellite tracking to recover the anomalies) are used in order to overcome problems of indeterminacy of state vectors coupled with anomalies. The state vectors are thus estimated separately using these ancillary data types, and the anomalies are estimated from the residuals, much in the same way as will be done here. Yet this can be a reason for leaving out the gradient term, since it is only explicitly needed in case state vector adjustments are estimated from the residuals as well. Leaving out the term greatly simplifies the equations and speeds up the
ARTICLE IN PRESS S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
computation. Nevertheless, including this term greatly improves the quality of the solution, and it should be noted that for instance Kahn et al. (1982) estimated gravity anomalies from satellite-to-satellite tracking data based on the formulation of Rapp (1974), which already mentions the use of the gradient term. In the case presented here, the force parameters ~ p will be gravity field parameters. Again following Vonbun et al. (1980), the gravitational sensitivity coefficients are computed using Newton’s equation of planetary motion, i.e. ~ x€ ¼ qV =q~ x, with the potential V expressed as V ¼ U þ T, where U is the normal potential of the reference gravity field expressed in spherical harmonics, and T is the disturbing potential expressed by Stokes’ integral formula, see Heiskanen and Moritz (1984): ZZ R TðPÞ ¼ SðP; QÞ DgðQÞ dsðQÞ, (3) 4p sðQÞ
where P is the point of computation, Q denotes a point on the reference surface s, SðP; QÞ is the Stokes kernel function and R the reference radius (here taken as 1738 km). The reference gravity field as expressed by U is known, so partials of U with respect to gravity parameters are zero. This leaves only a sensitivity coefficient q~ xðtÞ=qD~ g to be evaluated, together with partials with respect to the gravity anomalies, when Newton’s equation is differentiated with respect to the gravity parameters. This leads to the equation of the sensitivity coefficient q~ xðtÞ=qD~ g being expressed as d q~ xðtÞ qf~ðt; ~ xðtÞ; D~ gÞ q~ xðtÞ qf~ðt; ~ xðtÞ; D~ gÞ ¼ þ . dt qD~ g q~ xðtÞ qD~ g qD~ g
(4)
In order to evaluate the derivative of the forcing function f~ðt; ~ xðtÞ; D~ gÞ for each separate gravity anomaly, first Eq. (3) needs to be discretised for numerical computations. To this extent, the sphere is divided into equi-angular grid cells of 1 1 ; this corresponds to the resolution that is aimed at using Lunar Prospector data. The integral over each grid cell is then replaced by the value of the Stokes function at the central point of the grid cell times the area of the grid cell times the mean value of the anomaly over the block. This discretisation then means that the midpoint rule is used, under the assumption that the integral can be replaced by an average value of the integrand over the grid cell. Simulations using both the forward problem (computation of accelerations at satellite altitude from grids of anomalies) and the inverse problem (the determination of anomalies from satellite tracking data) have shown this works very well with the chosen 1 1 grid. A subdivision of the grid cells, as mentioned in e.g. Kaula (1996) did not show much improvement, although it should be noted that this depends on the input signal as well. There are different ways of dealing with a changing Stokes function over the grid cell, such as using the
1333
aforementioned subdivision, by computing a solution for smaller blocks and then averaging them to larger blocks (e.g. Kahn et al., 1982), or by expanding the Stokes function into a Taylor series such that an analytical solution for the integral of the approximated function exists (Rudolph, 2000). However, since simulations showed that the 1 resolution is already small enough, no such subdivision is implemented here. Next, in order to compute the contribution of one single anomaly Dgi to qf~=qDg, the acceleration expressed in terms of anomalies is needed. This can be obtained by differentiating the discretised version of Eq. (3). For the contribution of one anomaly to the radial acceleration dri this leads to qdri ðPÞ R qSðP; Qi Þ dsi . ¼ qDgðQi Þ 4p qr
(5)
This is the core of the expression for qf~=qDg. Similar expressions as given in Eq. (5) can be obtained for the other two directions, thus completing the acceleration contributions in an orthogonal system. These can be transformed into the desired inertial coordinate system needed for integration of the variational equations. After the sensitivity coefficients are computed, a matrix system in the form of y ¼ Ax can be formed. The vectors x and y will contain the anomalies (and possibly state vector adjustments as well) and the measurements in terms of range and Doppler residuals, respectively. The entries in the design matrix A will be filled by subsequent multiplication of the geometrical partials with the sensitivity coefficients as shown in Eq. (1). The equation system can then be solved by standard tools of linear algebra.
3. Simulation results In order to evaluate this recovery method with respect to the accuracy that can be achieved, different simulation scenarios have been conducted. Together with the initial set-up, the results of these simulations will be discussed in this section. 3.1. Simulation set-up The tracking data residuals are generated by subtracting measurements of range and range-rate of two orbit generation runs with different force models. By subtracting the measurements from two orbit propagations using different force models, the simulated residuals will contain the information from the difference in force models. Both orbits will initially have the same starting position, assuring there will be no effects from a difference in state vector. The used force models depend on the effect under study. Both orbits can for instance use the same spherical harmonical model, up to
ARTICLE IN PRESS 1334
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
Table 1 A brief summary of the simulation set-up, including used force models, data, orbit and recovery grid Force models Reference orbit
LP75G, l ¼ 2 75, at 50 km altitude
New orbit Truth gravity
LP75G, l ¼ 2 75 þ extra force parameters that will be estimated in terms of anomalies Reference model LP75G þ anomalies LP150Q, l ¼ 76 150
Measurements Observations Sampling rate Stations Visibility Data weights Time period
2-way tracking, range and range-rate Every second, both range and range-rate data DSN stations 24, 42 and 66, coordinates as published by DSN Elevation check; cut-off angle for elevation of 10 One set using 1.0 m for range data and 1.0 m/s for range-rate data, called case A, and another set using 1.0 m for range data and 1.0 mm/s for range-rate data, called case B July 10, 1999, start at 00:00:00.0, one day duration
Orbit Initial elements Integrator Step size
a ¼ 1788 km, e ¼ 0:001, i ¼ 90 , O ¼ 101:696 , o ¼ 1:0 , M ¼ 0 10th order Adams–Bashfort–Moulton 1 second
Recovery grid Coverage Resolution Tracks
l ¼ 0 10 , f ¼ 0 10 (coordinates of central points) 1 1 , 121 anomalies 11 tracks over the grid; one over each column of anomalies
the same degree, or up to different degrees to study systematical errors on the data. Extra sets of anomalies can also be used for the second orbit. This will be explained case by case. Orbits are generated using a 10th order Adams–Bashfort–Moulton integrator (see Press et al., 1997) with a start-up routine to initialise the integrator (see Visser, 1992), using a 1-s step size. Data are generated each second. Only two-way tracking data are used. The stations used are from the DSN network, and they are stations 24, 42 and 66, using coordinates as published by DSN, see Thornton and Border (2000). The reference orbit will be initially at 50 km height above the lunar surface, simulating the low altitude of the Lunar Prospector mission. It will also be a polar orbit, with a very small eccentricity (e ¼ 0:001). It will be located directly above the grid under investigation on the lunar surface. For these simulation purposes, this grid will be in a 1 1 resolution, spanning fl; fg ¼ ½0 ; 10 , which are the coordinates of the central points. The reference force model that will be used is the lunar gravity field model LP75G (Konopliv et al., 1998), a spherical harmonical expansion up to degree and order 75. This force model will be used as starting model in the computation of the range and range-rate residuals, and also as the reference model for the computation of the variational equations. Input anomalies will be generated from the latest spherical harmonic product derived from Lunar Prospector data, LP150Q, complete up to degree and order 150. The truth model will then consist of the reference model LP75G plus anomalies from LP150Q using l ¼ 76 150. Since the reference model is the
same in the data generation and the recovery, the recovered anomalies must be those from LP150Q, l ¼ 76 150. This set-up information is also given in Table 1 for quick reference. Normal matrices will be formed for each track or arc, where a track or arc is defined as the period it takes the satellite to cross the geographical area under consideration. This means that only those residuals directly above the recovery grid will be used. Due to the localisation properties of anomalies, data points directly above the grid will also contain most information about these anomalies. This will also lead to a short-arc approach, since relatively small grids will be solved for. Such a short-arc approach is also explained in Rowlands et al. (2002) and used for the determination of local gravity in Rowlands et al. (2005), although it should be pointed out that they studied and used satellite-to-satellite tracking (SST) data (in the framework of the GRACE mission), which have a relatively larger signal at higher frequencies and thus are more correlated with local features than Lunar Prospector-type data are. The distinction should thus be made that in this work Earth-based tracking data are used, from which it is more difficult to extract the local gravity field. 3.2. Linearisation errors The observation equation as given in Eq. (1) is derived from a linear variational approach, and thus linearisation errors will arise. To assess these, simulations according to the set-up given in Table 1 have been conducted. No errors are included on the data. The
ARTICLE IN PRESS S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
extra force parameters in the new orbit are anomalies from the LP150Q model from l ¼ 76 150. No other parameters or errors are included. The input anomalies have an rms of 21.54 mGal. Fig. 1 shows the differences between recovered and input anomalies for the two different data weights as included in Table 1. This figure shows the recovery of the anomalies is remarkably good: the rms of differences is 0.030 mGal for data weight A and in case the rangerate data is given relatively more weight, the rms becomes as low as 1:5 10 4 mGal. These results show the linearisation errors are small, especially in the case range-rate data is emphasised in the recovery. It also shows that the linearisation errors for range data have a bigger impact on the recovery than the linearisation errors for range-rate , which is connected with the fact that range data are one more integration further from accelerations than range-rate data are. The biggest differences occur at the top of the grid; this is because
1335
the satellite flies from south to north, and the linearisation error builds up along the track. This can be seen more clearly by depicting the actual linearisation error lin along the track, which can be done by computing gtrue . lin ¼ y AD~
(6)
Fig. 2 shows this error, for both range and range-rate data. This figure shows that the linearisation errors indeed are small. The range data linearisation errors are almost at machine precision for the first half of the arcs. For range-rate data, the build-up along the track is shown very clearly. The influence of the gradient term as mentioned in Section 2 has also been assessed, by computing solutions with and without this term taken into account in the computation of the design matrix. This showed that this term greatly reduces the linearisation error in the solution. This is especially seen when data sets from different time periods are combined. When the term is
= 1.0 m, · = 1.0 m/s
= 1.0 m, · = 1.0 mm/s
x 10-6 Linearisation error for range data 3
x 10-8 Linearisation error for range-rate data 6
2
4
1
2
Error [m/s]
Error [m]
Fig. 1. Differences between recovered and input anomalies when two different weights for the data are used. The rms of differences for the results on the left is 0.030 mGal, for the results on the right it is 1:5 10 4 mGal. The input anomalies have an rms of 21.54 mGal.
0 -1
-2 -4
-2 -3
0
-6 0
2
4 6 Latitude [degrees]
8
10
-8
0
2
4 6 Latitude [degrees]
Fig. 2. Linearisation errors for range and range-rate data in the recovery method.
8
10
ARTICLE IN PRESS 1336
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
not taken into account, the linearisation error is bigger, and the combination of two such data sets amplifies in this case the recovery error. This shows the importance of the gradient term, since it keeps the linearisation error negligible. 3.3. Errors in the initial state vector An adjustment to the state vector can be estimated by adding the term q~ x=q~ x0 to the observation equation, as is done in Eq. (1). Here it is investigated how adding such a parameter interacts with the estimation of the anomalies. These state vector parameters will be arc dependent parameters, since for each arc a new set is estimated. Errors in the initial state vector are generated by adding a random number from a Gaussian distribution times a standard deviation to the initial state vector for the propagation of the new orbit which includes the anomalies as extra force parameters, see also Table 1. The initial state is then corrupted and this will lead to different observations. The level used here is 0.1 m for the position and 1.0 mm/s for the velocities, which was chosen to represent range and range-rate data errors. The rest of the set-up is the same as the one mentioned in Table 1; no other errors than the state vector errors are applied. With these errors, the rms of differences between input and recovered anomalies becomes 0.0138 mGal, which is still a very good recovery, but worse than the error-free solution which had an rms of differences of 1:50 10 4 mGal. In Fig. 3, the differences between the input and recovered state vector adjustments D~ x0 are shown, splitted into the norm of the position and velocity for each arc. This figure shows that indeed, as was already found by Vonbun et al. (1980), it is extremely difficult to estimate a state vector adjustment from such a short arc. The figure shows that for many tracks the recovery error is bigger than the initial norm
of D~ x0 . The recovery results might be improved by iterations, but this has not been implemented here. These results show the need for accurate orbits prior to the anomaly recovery. Such orbits can be determined by first processing the tracking data and doing precise orbit determination, and then using the residuals to determine the anomaly adjustments in a next step. The state vector part of the normal matrices for each track has also been investigated, and eigenvalue analysis showed that there are actually only two degrees of freedom: 3 or 4 eigenvalues are numerically close to zero. For the non-zero eigenvalues, analysis of the according eigenvectors revealed that the observable combinations are Dx þ Dy þ Dz and Dx_ þ Dy_ þ D_z. If these are estimated instead of a full six parameter state vector adjustment, however, the combinations are also not resolved well. It is interesting to note that in Rowlands et al. (2002) a similar approach is followed for low-low SST, yet more thoroughly exploited, where a transformation of 12 Cartesian coordinates (since two satellites are dealt with) into coordinates of the baseline vector of both satellites and the coordinate of the baseline vector midpoint are used in a spherical coordinate system, allowing a decoupling of the tasks of orbit determination and gravity field determination. Signal strength of the errors on the state vector is of great influence on the recovery error of the anomalies. To assess this influence, another test has been considered, where the state vector errors on position and velocities have been chosen according to the rule of thumb spos =svel / n, with n being the satellite’s mean motion in rad/s. At 50 km above the lunar surface, n ¼ 9:26 10 4 rad=s. A value of spos ¼ 10 m, chosen to represent positioning error, leads to svel ¼ 9:26 mm=s. With all 6 state vector parameters estimated, the rms of differences with the truth anomalies now is 0.284 mGal. The recovery error itself has a slight long-wavelength signal in it, as did the results using the smaller position and velocity errors of 0.1 m and 1.0 mm/s, but at a level
Fig. 3. Plot of the rms of D~ x0 for each arc, splitted for the position and velocity part. Shown are the norms of D~ x0 of the original, true state vector adjustments, the estimated adjustment, and the difference between these two.
ARTICLE IN PRESS S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
1337
below 1 mGal it is still very good. The results for the state vector adjustment estimates are similar as shown in Fig. 3. Comparing this result with the previous result shows the anomaly recovery indeed is directly linked to the size of the state vector error. One final way to assess the influence is to apply this error to the orbit used in the computation of the variational equations. This means a slightly different reference orbit is used in the computation of the design matrix. The data themselves remain uncorrupted, thus the state vector effect is not present there. The rms of differences with the truth anomalies then becomes 0.0663 mGal. Concluding, the state vector adjustments themselves cannot be recovered using this short-arc approach. Then again, this approach is not intended to be used for precise orbit determination. The results show that prior to the regional recovery, precise orbits are required. Even though state vector errors have their influence on the recovered anomalies, by taking care of them in the anomaly recovery process, the accuracy of the recovery is still very good at levels below 1 mGal.
359˚0˚ 1˚ 2˚ 3˚ 4˚ 5˚ 6˚ 7˚ 8˚ 9˚ 10˚11˚12˚13˚14˚15˚16˚17˚18˚19˚20˚21˚ 21˚ 21˚ 20˚ 20˚ 19˚ 19˚ 18˚ 18˚ 17˚ 17˚ 16˚ 16˚ 15˚ 15˚ 14˚ 14˚ 13˚ 13˚ 12˚ 12˚ 11˚ 11˚ 10˚ 10˚ 9˚ 9˚ 8˚ 8˚ 7˚ 7˚ 6˚ 6˚ 5˚ 5˚ 4˚ 4˚ 3˚ 3˚ 2˚ 2˚ 1˚ 1˚ 0˚ 0˚ -1˚ -1˚ 359˚0˚ 1˚ 2˚ 3˚ 4˚ 5˚ 6˚ 7˚ 8˚ 9˚ 10˚11˚12˚13˚14˚15˚16˚17˚18˚19˚20˚21˚
3.4. Near- and far-field contributions to the data
Fig. 4. Differences between the input and recovered anomalies in the case of far field contributions to the data. The rms of differences with the input anomalies from LP75G, l ¼ 61 75 is 12.88 mGal. The input anomalies themselves have an rms of 21.50 mGal.
Since the purpose is to recover regional gravity fields, the integration over the sphere as given in Eq. (3) is limited to a certain area of influence. Naturally, an omission error is introduced, since the data will also contain signal belonging to anomalies that are not accounted for in the recovery. This will be most clear on the boundaries of the grid. In order to show the boundary effect in the simulations more clearly, the recovery grid will be larger: it will span anomaly cells with their central points at fl; fg ¼ ½0 ; 20 . This means a set of 21 21 anomalies. The residuals will be computed using a reference orbit that is generated using the LP75G gravity field model up to degree and order 60, and a new orbit using l max ¼ 75. The recovered anomalies should thus represent the part missing from the reference orbit, i.e. the spherical harmonical spectrum from l ¼ 61 75. Anomalies from the LP75G model using l ¼ 61 75 have an rms of 21.50 mGal. Since spherical harmonics are global functions, there will be global information present in the data, that will be mapped into an area of limited extensions. In Fig. 4 the differences between the input and recovered anomalies are shown. This figure shows clearly the boundary effect. The standard way to deal with this is to leave out the boundary part, and thus retain only a central part. In Table 2 different sizes of such central areas are listed, together with the rms of differences. This shows that the initial rms of differences for the whole area of 12.88 mGal can be brought back to 1.01 mGal in case only the central part is taken along. So far, all solutions have been created without regularisation. However, the problem at hand is an illposed problem, mostly in the sense of stability, since
mgal -50 -40 -30 -20 -10
0
10 20 30 40 50
Table 2 Differences with the true solution for different central areas of the solution. Values are in mGal Area Input D~ g l ¼ 0–20 , l ¼ 3–17 , l ¼ 4–16 , l ¼ 5–15 ,
min/max/rms f ¼ 0–20 f ¼ 3–17 f ¼ 4–16 f ¼ 5–15
60.31/58.34/21.50 63.70/49.41/12.88 11.53/12.73/2.96 6.77/5.81/1.80 3.48/4.43/1.01
downward continuation of errors at satellite altitude will lead to amplified errors in the recovered anomalies on the reference surface. Thus, regularisation might be required to turn the problem into a well-posed one. Tikhonov regularisation is used, where the solution can be computed by means of: x^ a ¼ ðAT WA þ aIÞ 1 AT W~ y
(7)
with x^ a the estimated solution, W the weight matrix for the data, I the unity matrix and a a parameter that needs to be chosen. Other regularisation methods are also possible; especially worth mentioning is the use of a full variance–covariance matrix, as is applied in Barriot and Balmino (1992). Using a diagonal regularisation matrix assumes there are no correlations between the anomalies, while in reality, there are. If this is to be taken into account however, one will also have to prescribe a
ARTICLE IN PRESS 1338
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
certain structure of the signal, apart from the signal strength that is controlled when using a diagonal regularisation matrix. In the case of the Moon, prescribing this structure can be difficult, since the global models themselves rely heavily on the Kaula rule mentioned in Section 1. This is why here it is chosen to use only a diagonal matrix. Tests with a non-diagonal matrix have also been performed, where the anomalies are constrained by taking into account values of their direct neighbours, but these tests showed no improvement was to be gained by this. This was also found for instance by Ilk (1993) for the case of SST data for an Earth-orbiting satellite. The value for a needs to be determined, and in these simulations this can be done by trial-and-error, since the input truth anomalies are known. In reality, the truth is not known, but simulations like these can be used to estimate reasonable values for the regularisation parameter a. This results in a ¼ 10þ7 giving the smallest rms of differences with the input anomalies, namely 6.08 mGal, compared to the 12.88 mGal when no regularisation is used. These numbers are valid for the whole solution. In case the central part is examined, it turns out the regularised solution has an rms of differences of 1.58 mGal, compared to the 1.01 mGal of the unregularised solution. It should be noted that the regularisation matrix is also applied to the whole solution. In case the trial-and-error method is again applied, but now comparing the central part of the regularised solution with the central part of the input anomalies, it turns out that a ¼ 10þ6 leads to an rms of differences of 0.508 mGal, showing that indeed regularisation can help to improve the quality of the solution. Nevertheless, deriving a regularisation matrix which will only be applied to the central part can also be worthwhile to improve the solution even further. The results show that recovery down to mGal level is still possible, yet it should be noted this highly depends on the input signal as well. In this case, the new orbit was generated using l max ¼ 75, which means a resolution of roughly 180=75 ¼ 2:4 , which is larger than the 1 aimed for. Therefore, the numbers will be different in case there will be signal beyond the aimed-for resolution. Then one might also consider subdivision of the grid cells when forming the normal matrices as mentioned in Section 2. This has not been further investigated here. 3.5. Stochastic errors on the data Apart from possible systematic errors, such as farfield contributions, observation errors such as noise will obviously also be of influence on the recovery results. To simulate this, the data are contaminated by adding a random number from a Gaussian distribution times a certain standard deviation. In dealing with noise on the
-
-
Fig. 5. The rms of differences using all possible combinations of the 4 data sets versus the regularisation parameter when Tikhonov regularisation is applied. The average altitude of the satellite orbit is 50 km. The error signals used on the data are sr ¼ 0:5 m and sr_ ¼ 1:0 mm=s.
data coming from a random process, the more data are available, the better. For the simulations presented here, 4 data spans will be generated: one starting at June 10 1999 (990610), one at June 20 1999 (990620), one at July 10 1999 (990710) and one at July 20 1999 (990720). The noise on these data will have a realistic level: Konopliv et al. (2001) reported a data noise level of 0.5 m for range data and 1.0 mm/s for Doppler data, which are the levels used here. The recovery grid will again be the same as mentioned in Table 1, consisting of 11 11 anomalies. With such error levels, regularisation will be indispensable. Again, Tikhonov regularisation will be applied, with a free to choose. Fig. 5 shows the rms of differences with the input anomalies using all possible combinations of the data, for various values of the regularisation parameter a. This figure shows there actually is an optimal regularisation parameter, which is an important result when this method will be used on real Lunar Prospector data. The optimal parameter is a ¼ 10þ7 and produces an rms of differences of 13.52 mGal when all data sets are used. Compared to the rms of the input anomalies of 21.54 mGal, this result shows most of the signal can be recovered. The correlation of 0.78 between recovered anomalies and input anomalies confirms this. 3.6. Influence of altitude The results in the previous sub-section were all generated using an average satellite altitude of 50 km above the lunar surface. Even though this is already very low, Lunar Prospector flew at even lower altitudes
ARTICLE IN PRESS S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
1339
these data sets, the rms of differences without regularisation becomes 3.79 mGal, slightly better than the value of 3.86 mGal when using only the 20 km data. The most important result however is that the effects of downward continuation remain limited.
4. Conclusions
-
-
Fig. 6. Rms of differences between Tikhonov regularised solutions and the input anomalies versus the regularisation parameter a, for all possible combinations of the 4 data sets. The average altitude of the satellite is now 20 km, as opposed to 50 km in Fig. 5. The level of noise is the same.
during its Extended Mission. Therefore, the same simulations as described above have been conducted for an average altitude of 20 km. The data were again contaminated with noise with standard deviations sr ¼ 0:5 m and sr_ ¼ 1:0 mm=s, and regularised solutions were computed using different values for the regularisation parameter. The results are shown in Fig. 6. This figure again shows there is a minimum in the rms of differences between recovered and input anomalies, but what is more important is that the figure shows that when extremely low-altitude data are included, accurate solutions close to the optimum can be obtained without regularisation. This is because downward continuation becomes less and less of an issue with decreasing altitude. The best rms of differences is again obtained for a ¼ 10þ7 ; this is also an important result since the optimal regularisation parameter (in the sense of the smallest rms of differences) remains the same with varying altitude (it will however depend on the level of actual noise in the data). The best rms of differences is 3.60 mGal, which is already much better than the value of 13.52 mGal for the 50 km data. Without regularisation, the rms of differences is 3.86 mGal. This also shows that regularisation in the sense of including a priori information in the form of a regularisation matrix is not necessarily required in order to obtain a good solution at a resolution of 1 when low-altitude data are available. The improvement is not that much in terms of the rms of differences when using regularisation, and the cost would be to find the correct regularisation parameter. The combination of both the 50 km data and the 20 km data has also been investigated. When using both
A recovery method based on a linear variational approach, has been developed and thoroughly tested. This method allows gravity anomalies to be determined from Earth-based tracking data residuals. The method has been extensively tested by means of simulations. A short-arc approach has been used, where one arc or track consists of the period it takes the satellite to fly over the area of interest. It has been tried to estimate state-vector adjustments from such a short-arc, but this did not lead to satisfactory results. This means that prior to the determination of the anomalies, a precise orbit has to be determined, using elaborate reference models. In the presence of various error sources, recovery down to the level of several mGal is still possible. In the presence of noise on data taken at an average altitude of 50 km, regularisation will be indispensable to obtain a high-quality solution. An optimal regularisation parameter can be found, as is shown by comparing the recovered solutions with the input anomalies. The best rms of differences between input and recovered anomalies is 13.52 mGal, when using a regularisation parameter a ¼ 10þ7 . Compared to an input rms of 21.54 mGal, this shows most of the input signal can be recovered. When low altitude data are included, at an average of 20 km (which are available for the Lunar Prospector mission), regularisation is no longer required to obtain a good solution at a resolution of 1 , since the effects of downward continuation can be limited. The rms of differences without regularisation is 3.86 mGal, and with regularisation it is 3.60 mGal. This shows that by including low-altitude data, the recovery method becomes stable and hardly requires regularisation. These results show that the recovery method is capable of extracting the high-frequency gravity information from tracking data residuals. This method will therefore be used for processing real Lunar Prospector data, which will be the topic of a follow-on paper.
References Barriot, J., Balmino, G., 1992. Estimation of local planetary gravity fields using line of sight gravity data and an integral operator. Icarus 99 (1), 202–224. Barriot, J.-P., Balmino, G., Vale`s, N., 1997. Building reliable local models of the Venus gravity field from the Cycles 5 and 6 of the Magellan LOS gravity data. Geophysical Research Letters 24 (4), 477–480.
ARTICLE IN PRESS 1340
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
Floberghagen, R., 2002. Lunar Gravimetry. vol. 273, Astrophysics and Space Science Library. Kluwer Academic Publishers, Dordrecht, The Netherlands, 304 pp. Heiskanen, W.A., Moritz, H., 1984. Physical Geodesy. Institute of Physical Geodesy, Graz, Austria (original edition by W.H. Freeman and Company, San Francisco, 1967). Ilk, K., 1993. Regularization for High Resolution Gravity Field Recovery by Future Satellite Techniques. In: Anger, G., Gorenflo, R., Jochmann, H., Moritz, H., Webers, W. (Eds.), Inverse Problems: Principles and Applications in Geophysics, Technology and Medicine. vol. 74 of Mathematical Research. Akademie, Berlin, pp. 189–214. Kahn, W., Klosko, S., Wells, W., 1982. Mean gravity anomalies from a combination of Apollo/ATS 6 and GEOS 3/ATS 6 SST tracking campaigns. Journal of Geophysical Research 87 (B4), 2904–2918. Kaula, W., 1996. Regional gravity fields on Venus from tracking of Magellan cycles 5 and 6. Journal of Geophysical Research 101 (E2), 4683–4690. Kaula, W.M., 1966. Theory of Satellite Geodesy, Applications of Satellites to Geodesy. Blaisdell Publishing Company. Konopliv, A., Binder, A., Hood, L., Kucinkas, A., Sjogren, W., Williams, J., 1998. Improved gravity field of the moon from lunar prospector. Science 281 (5382), 1476–1480. Konopliv, A., Asmar, S., Carranza, E., Sjogren, W.L., Yuan, D., 2001. Recent gravity models as a result of the lunar prospector mission. Icarus 150, 1–18. Lemoine, F., Smith, D., Zuber, M., Neumann, G., Rowlands, D., 1997. A 70th degree lunar gravity model (GLGM-2) from Clementine and other tracking data. Journal of Geophysical Research 102 (E7), 16,339–16,359. Montenbruck, O., Gill, E., 2000. Satellite Orbits. Springer, Heidelberg. Muller, P., Sjogren, W., 1968. Mascons: Lunar Mass Concentrations. Science 161, 680–684. Muller, P., Sjogren, W., Wollenhaupt, W., 1974. Lunar Gravity: Apollo 15 doppler radio tracking. The Moon 10, 195–205. Press, W., Teukolsky, S., Vetterling, W., Flannery, B., 1997. Numerical Recipes in Fortran 77: The Art of Scientific Computing. Cambridge University Press. 2nd ed. vol. 1. Fortran Numerical Recipes. Rapp, R., 1974. Procedures and Results related to the direct determination of gravity anomalies from satellite and terrestrial
gravity data. Reports of the Department of Geodetic Science, vol. 211, Ohio State University. Rowlands, D., Ray, R., Chinn, D., Lemoine, F., 2002. Short-arc analysis of intersatellite tracking data in a gravity mapping mission. Journal of Geodesy 76, 307–316. Rowlands, D., Luthcke, S., Klosko, S., Lemoine, F., Chinn, D., McCarthy, J., Cox, C., Anderson, O., 2005. Resolving mass flux at high spatial and temporal resolution using GRACE intersatellite measurements. Geophysical Research Letters 32 (L04310) doi:10.1029/2004GL021908. Rudolph, S., 2000. Regionale und globale Gravitationsfeldanalyse hochauflo¨sender Satellitendaten mittels mehrgitterverfahren. Ph.D. thesis, Rheinische Friedrich-Wilhelms-Universita¨t Bonn, Mu¨nchen, deutsche Geoda¨tische Kommission,Dissertationen, Reihe C, Heft Nr. 528. Sjogren, W., Muller, P., Wollenhaupt, W., 1972. Apollo 15 Gravity Analysis from the S-band transponder experiment. The Moon 4, 411–418. Sjogren, W., Wimberly, R., Wollenhaupt, W., 1974a. Lunar Gravity: Apollo 16. The Moon 11, 35–40. Sjogren, W., Wimberly, R., Wollenhaupt, W., 1974b. Lunar Gravity: Apollo 17. The Moon 11, 41–52. Sjogren, W., Wimberly, R., Wollenhaupt, W., 1974c. Lunar Gravity via the Apollo 15 and 16 Subsatellites. The Moon 9, 115–128. Sugano, T., Heki, K., 2004a. High resolution gravity anomaly map from the lunar prospector line-of-sight acceleration data. Earth, Planets and Space 56, 81–86. Sugano, T., Heki, K., 2004b. Isostasy of the Moon from highresolution gravity and topography data: Implication for its thermal history. Geophysical Research Letters 31 (24) doi: 10.1024/ 2004GL022059. Thornton, C., Border, J., 2000. Radiometric Tracking Techniques for Deep-Space Navigation. Monograph 1, Deep-Space Communications And Navigation Series, Jet Propulsion Laboratory, JPL Publication 00–11. Visser, P., 1992. The Use of Satellites in Gravity Field Determination and Model Adjustment. Ph.D. Thesis, Delft University of Technology. Vonbun, F., Kahn, W., Wells, W., Conrad, T., 1980. Determination of 5 5 gravity anomalies using satellite-to-satellite tracking between ATS-6 and Apollo. Geophys. J. R. Astr. Soc. 61, 645–657.
Planetary and Space Science 53 (2005) 1331–1340 www.elsevier.com/locate/pss
A method to determine regional lunar gravity fields from earth-based satellite tracking data S. Goossens,1, P.N.A.M. Visser, B.A.C. Ambrosius Department of Earth Observation and Space Systems, section Astrodynamics and Satellite Systems, Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The Netherlands Received 17 March 2005; received in revised form 23 June 2005; accepted 29 June 2005 Available online 18 August 2005
Abstract A method to determine regional gravity fields of the Moon from Earth-based Doppler and range satellite tracking data residuals of a low Moon-orbiting satellite has been developed and thoroughly tested in a controlled simulation environment. A short-arc approach, where one arc consists of the time it takes the satellite to cross the grid of interest on the lunar surface, is used in order to filter out most long-wavelength signal that can still be present in the residuals. Simulation results where the data are contaminated with either typical systematic or stochastic noise show that recovery of the local gravity field down to the level of several mGal is possible. The inclusion of extremely low-altitude data also means that regularisation in the sense of including a priori information in the form of a regularisation matrix is not necessary in order to obtain a good solution at high resolution. r 2005 Elsevier Ltd. All rights reserved. Keywords: Regional gravity field determination; Lunar gravity field; Gravity anomalies
1. Introduction The Moon and the Earth are in a spin-orbit resonance, which causes the Moon to always have the same face towards the Earth. This severely limits Earth-based tracking of a satellite orbiting the Moon, since no tracking data are available when the satellite flies over the far side of the Moon. This in its turn hampers the determination of the global gravity field of the Moon in terms of spherical harmonics, since less than one half of the Moon is lacking data coverage. By using the mathematical tool of regularisation, the ill-posed inverse problem can be transformed into a well-posed one, and a solution can be obtained. For the Moon, such solutions have been determined in the past, e.g. Lemoine et al. (1997) and Konopliv et al. (2001) using relatively Corresponding author.
E-mail addresses: [email protected], [email protected] (S. Goossens). 1 Currently at National Astronomical Observatory of Japan, Mizusawa. 0032-0633/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.pss.2005.06.009
recent data from the Clementine and Lunar Prospector spacecraft, respectively, and these solutions perform very well in terms of data fit and correlation with topographical features. In order to determine these solutions, a priori information in the form of a Kaula rule is used, which is based on the observation that the spectrum of the spherical harmonic coefficients of the gravity field behave according to b=l 2 , with l the degree and b a constant, see Kaula (1966). However, it has been shown by Floberghagen (2002) that the choice of the regularisation parameter or which Kaula constant to apply, is of great influence on the final result. When using global basis functions such as spherical harmonics, it cannot be avoided that the solution for the near side is influenced by that for the far side, and vice versa. Using global basis functions also limits the resolution that can be obtained. However, the current data set with excellent near side coverage with Lunar Prospector data suggests more information is present in the data than currently can be modelled using spherical harmonics alone (Konopliv et al., 2001). Due to the low
ARTICLE IN PRESS 1332
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
altitude of Lunar Prospector, at an average of 30 km in its Extended Mission, a resolution of 1 on the lunar surface should be achievable. Current spherical harmonics expansions are limited to degree 165 maximum, and this solution has been estimated in multiple steps (Konopliv et al., 2001). Thus, in order to exhaust the information present in the data and to overcome problems inherent to the use of global basis functions, a different representation needs to be used. Such a representation is expected to achieve a higher resolution for areas on the near side, which will aid the interpretation of selenophysical features. Gravity is represented in this paper by means of gravity anomalies on the lunar surface. A method to determine these anomalies from range and Doppler tracking data residuals has been developed and thoroughly tested. In previous efforts the data type used to determine gravity anomalies has often been line-of-sight (LOS) accelerations, which can be obtained by differentiating Doppler residuals with respect to time. Examples for the Moon include Muller and Sjogren (1968), who revealed the discovery of the mascons on the Moon, and Sjogren et al. (1972, 1974a–c) and Muller et al. (1974), who used LOS data due to coverage limitations and achieved a further improvement of knowledge of such features as Mare Serenitatis and Mare Crisium. Recent efforts for the Moon using LOS data of Lunar Prospector can be found in Sugano and Heki (2004a,b). In this paper however, the recovery will be done directly from range and Doppler data residuals. This is done in order to avoid problems with noise on the data which will be amplified in the differentiation process and to avoid problems coming from a possible lack of continuous data series. Similar problems were also noted by Barriot et al. (1997), who then also decided to avoid the delicate matters of numerical differentiation. This paper deals with the recovery method itself, and by means of extensive simulations in a controlled environment it is shown what level of recovery can be achieved. In a follow-on paper, real Lunar Prospector data will be dealt with. The simulations in this paper thus keep Lunar Prospector and its characteristics in mind. The paper is structured as follows: the recovery method itself will be explained in Section 2. In Section 3, simulation results with different error sources on the data will show the achievable recovery level. Section 4 contains the conclusions.
2. Recovery method The observation model that relates the tracking data residuals to gravity anomalies is formed by a linear variational approach. The approach used is an extension of the approach used by Vonbun et al. (1980). Let r_ be
the initial Doppler measurement. A residual Dr_ is obtained by comparing a computed measurement, created from reference models, with the actual observation. Subsequently, adjustments to the reference models can be made until a certain desired level of for instance data fit is achieved, and the final residuals are obtained. In this way, most effects on such a residual can be filtered, and the final residual itself can be thought to be caused by local gravity anomalies D~ g or a change in the initial state vector, given by D~ x0 . In a linear variational form, by using a first order Taylor expansion, the dependency of the data residuals on changes in the initial state vector and any kind of dynamical force model parameters ~ p becomes: qr_ q~ x qr_ q~ x Dr_ ¼ D~ p, (1) D~ x0 þ q~ x q~ x0 q~ x q~ p where ~ x is the state vector itself and ~ x0 is the initial state _ x are geometrical partials and vector. The partials qr=q~ can easily be formed from geometrical relationships by expressing the range-rate in terms of the coordinates of the satellite and ground station. The partials q~ x=q~ x0 and q~ x=q~ p are sensitivity coefficients. These coefficients can be obtained by integration of the variational equations, which are derived by differentiation of the equation of motion, see e.g. Montenbruck and Gill (2000). Let the equation of motion be expressed as d~ xðtÞ=dt ¼ f~ðt; ~ xðtÞ; ~ pÞ, with ~ xðtÞ the full 6 parameter state vector consisting of positions and velocities, and t the time. The function f~ describes the forces acting on the satellite, and it is dependent on the time, the state vector itself, and the force parameters. By differentiation of the equation of motion with respect to the parameters ~ p, the differential equation for q~ xðtÞ=q~ p reads: d q~ xðtÞ qf~ðt; ~ xðtÞ; ~ pÞ q~ xðtÞ qf~ðt; ~ xðtÞ; ~ pÞ ¼ þ . dt q~ p q~ xðtÞ q~ p q~ p
(2)
The difference between the formulation as used in this paper with the formulation used by Vonbun et al. (1980) is an extra term in the differential equation for the sensitivity coefficient q~ x=q~ p. This extra term is the gradient term qf~=q~ x of the forces acting on the satellite, and this term does not appear in the formulation used by Vonbun et al. (1980). There, it is mentioned that ancillary data types (basically, ground-based tracking of a satellite since they use satellite-to-satellite tracking to recover the anomalies) are used in order to overcome problems of indeterminacy of state vectors coupled with anomalies. The state vectors are thus estimated separately using these ancillary data types, and the anomalies are estimated from the residuals, much in the same way as will be done here. Yet this can be a reason for leaving out the gradient term, since it is only explicitly needed in case state vector adjustments are estimated from the residuals as well. Leaving out the term greatly simplifies the equations and speeds up the
ARTICLE IN PRESS S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
computation. Nevertheless, including this term greatly improves the quality of the solution, and it should be noted that for instance Kahn et al. (1982) estimated gravity anomalies from satellite-to-satellite tracking data based on the formulation of Rapp (1974), which already mentions the use of the gradient term. In the case presented here, the force parameters ~ p will be gravity field parameters. Again following Vonbun et al. (1980), the gravitational sensitivity coefficients are computed using Newton’s equation of planetary motion, i.e. ~ x€ ¼ qV =q~ x, with the potential V expressed as V ¼ U þ T, where U is the normal potential of the reference gravity field expressed in spherical harmonics, and T is the disturbing potential expressed by Stokes’ integral formula, see Heiskanen and Moritz (1984): ZZ R TðPÞ ¼ SðP; QÞ DgðQÞ dsðQÞ, (3) 4p sðQÞ
where P is the point of computation, Q denotes a point on the reference surface s, SðP; QÞ is the Stokes kernel function and R the reference radius (here taken as 1738 km). The reference gravity field as expressed by U is known, so partials of U with respect to gravity parameters are zero. This leaves only a sensitivity coefficient q~ xðtÞ=qD~ g to be evaluated, together with partials with respect to the gravity anomalies, when Newton’s equation is differentiated with respect to the gravity parameters. This leads to the equation of the sensitivity coefficient q~ xðtÞ=qD~ g being expressed as d q~ xðtÞ qf~ðt; ~ xðtÞ; D~ gÞ q~ xðtÞ qf~ðt; ~ xðtÞ; D~ gÞ ¼ þ . dt qD~ g q~ xðtÞ qD~ g qD~ g
(4)
In order to evaluate the derivative of the forcing function f~ðt; ~ xðtÞ; D~ gÞ for each separate gravity anomaly, first Eq. (3) needs to be discretised for numerical computations. To this extent, the sphere is divided into equi-angular grid cells of 1 1 ; this corresponds to the resolution that is aimed at using Lunar Prospector data. The integral over each grid cell is then replaced by the value of the Stokes function at the central point of the grid cell times the area of the grid cell times the mean value of the anomaly over the block. This discretisation then means that the midpoint rule is used, under the assumption that the integral can be replaced by an average value of the integrand over the grid cell. Simulations using both the forward problem (computation of accelerations at satellite altitude from grids of anomalies) and the inverse problem (the determination of anomalies from satellite tracking data) have shown this works very well with the chosen 1 1 grid. A subdivision of the grid cells, as mentioned in e.g. Kaula (1996) did not show much improvement, although it should be noted that this depends on the input signal as well. There are different ways of dealing with a changing Stokes function over the grid cell, such as using the
1333
aforementioned subdivision, by computing a solution for smaller blocks and then averaging them to larger blocks (e.g. Kahn et al., 1982), or by expanding the Stokes function into a Taylor series such that an analytical solution for the integral of the approximated function exists (Rudolph, 2000). However, since simulations showed that the 1 resolution is already small enough, no such subdivision is implemented here. Next, in order to compute the contribution of one single anomaly Dgi to qf~=qDg, the acceleration expressed in terms of anomalies is needed. This can be obtained by differentiating the discretised version of Eq. (3). For the contribution of one anomaly to the radial acceleration dri this leads to qdri ðPÞ R qSðP; Qi Þ dsi . ¼ qDgðQi Þ 4p qr
(5)
This is the core of the expression for qf~=qDg. Similar expressions as given in Eq. (5) can be obtained for the other two directions, thus completing the acceleration contributions in an orthogonal system. These can be transformed into the desired inertial coordinate system needed for integration of the variational equations. After the sensitivity coefficients are computed, a matrix system in the form of y ¼ Ax can be formed. The vectors x and y will contain the anomalies (and possibly state vector adjustments as well) and the measurements in terms of range and Doppler residuals, respectively. The entries in the design matrix A will be filled by subsequent multiplication of the geometrical partials with the sensitivity coefficients as shown in Eq. (1). The equation system can then be solved by standard tools of linear algebra.
3. Simulation results In order to evaluate this recovery method with respect to the accuracy that can be achieved, different simulation scenarios have been conducted. Together with the initial set-up, the results of these simulations will be discussed in this section. 3.1. Simulation set-up The tracking data residuals are generated by subtracting measurements of range and range-rate of two orbit generation runs with different force models. By subtracting the measurements from two orbit propagations using different force models, the simulated residuals will contain the information from the difference in force models. Both orbits will initially have the same starting position, assuring there will be no effects from a difference in state vector. The used force models depend on the effect under study. Both orbits can for instance use the same spherical harmonical model, up to
ARTICLE IN PRESS 1334
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
Table 1 A brief summary of the simulation set-up, including used force models, data, orbit and recovery grid Force models Reference orbit
LP75G, l ¼ 2 75, at 50 km altitude
New orbit Truth gravity
LP75G, l ¼ 2 75 þ extra force parameters that will be estimated in terms of anomalies Reference model LP75G þ anomalies LP150Q, l ¼ 76 150
Measurements Observations Sampling rate Stations Visibility Data weights Time period
2-way tracking, range and range-rate Every second, both range and range-rate data DSN stations 24, 42 and 66, coordinates as published by DSN Elevation check; cut-off angle for elevation of 10 One set using 1.0 m for range data and 1.0 m/s for range-rate data, called case A, and another set using 1.0 m for range data and 1.0 mm/s for range-rate data, called case B July 10, 1999, start at 00:00:00.0, one day duration
Orbit Initial elements Integrator Step size
a ¼ 1788 km, e ¼ 0:001, i ¼ 90 , O ¼ 101:696 , o ¼ 1:0 , M ¼ 0 10th order Adams–Bashfort–Moulton 1 second
Recovery grid Coverage Resolution Tracks
l ¼ 0 10 , f ¼ 0 10 (coordinates of central points) 1 1 , 121 anomalies 11 tracks over the grid; one over each column of anomalies
the same degree, or up to different degrees to study systematical errors on the data. Extra sets of anomalies can also be used for the second orbit. This will be explained case by case. Orbits are generated using a 10th order Adams–Bashfort–Moulton integrator (see Press et al., 1997) with a start-up routine to initialise the integrator (see Visser, 1992), using a 1-s step size. Data are generated each second. Only two-way tracking data are used. The stations used are from the DSN network, and they are stations 24, 42 and 66, using coordinates as published by DSN, see Thornton and Border (2000). The reference orbit will be initially at 50 km height above the lunar surface, simulating the low altitude of the Lunar Prospector mission. It will also be a polar orbit, with a very small eccentricity (e ¼ 0:001). It will be located directly above the grid under investigation on the lunar surface. For these simulation purposes, this grid will be in a 1 1 resolution, spanning fl; fg ¼ ½0 ; 10 , which are the coordinates of the central points. The reference force model that will be used is the lunar gravity field model LP75G (Konopliv et al., 1998), a spherical harmonical expansion up to degree and order 75. This force model will be used as starting model in the computation of the range and range-rate residuals, and also as the reference model for the computation of the variational equations. Input anomalies will be generated from the latest spherical harmonic product derived from Lunar Prospector data, LP150Q, complete up to degree and order 150. The truth model will then consist of the reference model LP75G plus anomalies from LP150Q using l ¼ 76 150. Since the reference model is the
same in the data generation and the recovery, the recovered anomalies must be those from LP150Q, l ¼ 76 150. This set-up information is also given in Table 1 for quick reference. Normal matrices will be formed for each track or arc, where a track or arc is defined as the period it takes the satellite to cross the geographical area under consideration. This means that only those residuals directly above the recovery grid will be used. Due to the localisation properties of anomalies, data points directly above the grid will also contain most information about these anomalies. This will also lead to a short-arc approach, since relatively small grids will be solved for. Such a short-arc approach is also explained in Rowlands et al. (2002) and used for the determination of local gravity in Rowlands et al. (2005), although it should be pointed out that they studied and used satellite-to-satellite tracking (SST) data (in the framework of the GRACE mission), which have a relatively larger signal at higher frequencies and thus are more correlated with local features than Lunar Prospector-type data are. The distinction should thus be made that in this work Earth-based tracking data are used, from which it is more difficult to extract the local gravity field. 3.2. Linearisation errors The observation equation as given in Eq. (1) is derived from a linear variational approach, and thus linearisation errors will arise. To assess these, simulations according to the set-up given in Table 1 have been conducted. No errors are included on the data. The
ARTICLE IN PRESS S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
extra force parameters in the new orbit are anomalies from the LP150Q model from l ¼ 76 150. No other parameters or errors are included. The input anomalies have an rms of 21.54 mGal. Fig. 1 shows the differences between recovered and input anomalies for the two different data weights as included in Table 1. This figure shows the recovery of the anomalies is remarkably good: the rms of differences is 0.030 mGal for data weight A and in case the rangerate data is given relatively more weight, the rms becomes as low as 1:5 10 4 mGal. These results show the linearisation errors are small, especially in the case range-rate data is emphasised in the recovery. It also shows that the linearisation errors for range data have a bigger impact on the recovery than the linearisation errors for range-rate , which is connected with the fact that range data are one more integration further from accelerations than range-rate data are. The biggest differences occur at the top of the grid; this is because
1335
the satellite flies from south to north, and the linearisation error builds up along the track. This can be seen more clearly by depicting the actual linearisation error lin along the track, which can be done by computing gtrue . lin ¼ y AD~
(6)
Fig. 2 shows this error, for both range and range-rate data. This figure shows that the linearisation errors indeed are small. The range data linearisation errors are almost at machine precision for the first half of the arcs. For range-rate data, the build-up along the track is shown very clearly. The influence of the gradient term as mentioned in Section 2 has also been assessed, by computing solutions with and without this term taken into account in the computation of the design matrix. This showed that this term greatly reduces the linearisation error in the solution. This is especially seen when data sets from different time periods are combined. When the term is
= 1.0 m, · = 1.0 m/s
= 1.0 m, · = 1.0 mm/s
x 10-6 Linearisation error for range data 3
x 10-8 Linearisation error for range-rate data 6
2
4
1
2
Error [m/s]
Error [m]
Fig. 1. Differences between recovered and input anomalies when two different weights for the data are used. The rms of differences for the results on the left is 0.030 mGal, for the results on the right it is 1:5 10 4 mGal. The input anomalies have an rms of 21.54 mGal.
0 -1
-2 -4
-2 -3
0
-6 0
2
4 6 Latitude [degrees]
8
10
-8
0
2
4 6 Latitude [degrees]
Fig. 2. Linearisation errors for range and range-rate data in the recovery method.
8
10
ARTICLE IN PRESS 1336
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
not taken into account, the linearisation error is bigger, and the combination of two such data sets amplifies in this case the recovery error. This shows the importance of the gradient term, since it keeps the linearisation error negligible. 3.3. Errors in the initial state vector An adjustment to the state vector can be estimated by adding the term q~ x=q~ x0 to the observation equation, as is done in Eq. (1). Here it is investigated how adding such a parameter interacts with the estimation of the anomalies. These state vector parameters will be arc dependent parameters, since for each arc a new set is estimated. Errors in the initial state vector are generated by adding a random number from a Gaussian distribution times a standard deviation to the initial state vector for the propagation of the new orbit which includes the anomalies as extra force parameters, see also Table 1. The initial state is then corrupted and this will lead to different observations. The level used here is 0.1 m for the position and 1.0 mm/s for the velocities, which was chosen to represent range and range-rate data errors. The rest of the set-up is the same as the one mentioned in Table 1; no other errors than the state vector errors are applied. With these errors, the rms of differences between input and recovered anomalies becomes 0.0138 mGal, which is still a very good recovery, but worse than the error-free solution which had an rms of differences of 1:50 10 4 mGal. In Fig. 3, the differences between the input and recovered state vector adjustments D~ x0 are shown, splitted into the norm of the position and velocity for each arc. This figure shows that indeed, as was already found by Vonbun et al. (1980), it is extremely difficult to estimate a state vector adjustment from such a short arc. The figure shows that for many tracks the recovery error is bigger than the initial norm
of D~ x0 . The recovery results might be improved by iterations, but this has not been implemented here. These results show the need for accurate orbits prior to the anomaly recovery. Such orbits can be determined by first processing the tracking data and doing precise orbit determination, and then using the residuals to determine the anomaly adjustments in a next step. The state vector part of the normal matrices for each track has also been investigated, and eigenvalue analysis showed that there are actually only two degrees of freedom: 3 or 4 eigenvalues are numerically close to zero. For the non-zero eigenvalues, analysis of the according eigenvectors revealed that the observable combinations are Dx þ Dy þ Dz and Dx_ þ Dy_ þ D_z. If these are estimated instead of a full six parameter state vector adjustment, however, the combinations are also not resolved well. It is interesting to note that in Rowlands et al. (2002) a similar approach is followed for low-low SST, yet more thoroughly exploited, where a transformation of 12 Cartesian coordinates (since two satellites are dealt with) into coordinates of the baseline vector of both satellites and the coordinate of the baseline vector midpoint are used in a spherical coordinate system, allowing a decoupling of the tasks of orbit determination and gravity field determination. Signal strength of the errors on the state vector is of great influence on the recovery error of the anomalies. To assess this influence, another test has been considered, where the state vector errors on position and velocities have been chosen according to the rule of thumb spos =svel / n, with n being the satellite’s mean motion in rad/s. At 50 km above the lunar surface, n ¼ 9:26 10 4 rad=s. A value of spos ¼ 10 m, chosen to represent positioning error, leads to svel ¼ 9:26 mm=s. With all 6 state vector parameters estimated, the rms of differences with the truth anomalies now is 0.284 mGal. The recovery error itself has a slight long-wavelength signal in it, as did the results using the smaller position and velocity errors of 0.1 m and 1.0 mm/s, but at a level
Fig. 3. Plot of the rms of D~ x0 for each arc, splitted for the position and velocity part. Shown are the norms of D~ x0 of the original, true state vector adjustments, the estimated adjustment, and the difference between these two.
ARTICLE IN PRESS S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
1337
below 1 mGal it is still very good. The results for the state vector adjustment estimates are similar as shown in Fig. 3. Comparing this result with the previous result shows the anomaly recovery indeed is directly linked to the size of the state vector error. One final way to assess the influence is to apply this error to the orbit used in the computation of the variational equations. This means a slightly different reference orbit is used in the computation of the design matrix. The data themselves remain uncorrupted, thus the state vector effect is not present there. The rms of differences with the truth anomalies then becomes 0.0663 mGal. Concluding, the state vector adjustments themselves cannot be recovered using this short-arc approach. Then again, this approach is not intended to be used for precise orbit determination. The results show that prior to the regional recovery, precise orbits are required. Even though state vector errors have their influence on the recovered anomalies, by taking care of them in the anomaly recovery process, the accuracy of the recovery is still very good at levels below 1 mGal.
359˚0˚ 1˚ 2˚ 3˚ 4˚ 5˚ 6˚ 7˚ 8˚ 9˚ 10˚11˚12˚13˚14˚15˚16˚17˚18˚19˚20˚21˚ 21˚ 21˚ 20˚ 20˚ 19˚ 19˚ 18˚ 18˚ 17˚ 17˚ 16˚ 16˚ 15˚ 15˚ 14˚ 14˚ 13˚ 13˚ 12˚ 12˚ 11˚ 11˚ 10˚ 10˚ 9˚ 9˚ 8˚ 8˚ 7˚ 7˚ 6˚ 6˚ 5˚ 5˚ 4˚ 4˚ 3˚ 3˚ 2˚ 2˚ 1˚ 1˚ 0˚ 0˚ -1˚ -1˚ 359˚0˚ 1˚ 2˚ 3˚ 4˚ 5˚ 6˚ 7˚ 8˚ 9˚ 10˚11˚12˚13˚14˚15˚16˚17˚18˚19˚20˚21˚
3.4. Near- and far-field contributions to the data
Fig. 4. Differences between the input and recovered anomalies in the case of far field contributions to the data. The rms of differences with the input anomalies from LP75G, l ¼ 61 75 is 12.88 mGal. The input anomalies themselves have an rms of 21.50 mGal.
Since the purpose is to recover regional gravity fields, the integration over the sphere as given in Eq. (3) is limited to a certain area of influence. Naturally, an omission error is introduced, since the data will also contain signal belonging to anomalies that are not accounted for in the recovery. This will be most clear on the boundaries of the grid. In order to show the boundary effect in the simulations more clearly, the recovery grid will be larger: it will span anomaly cells with their central points at fl; fg ¼ ½0 ; 20 . This means a set of 21 21 anomalies. The residuals will be computed using a reference orbit that is generated using the LP75G gravity field model up to degree and order 60, and a new orbit using l max ¼ 75. The recovered anomalies should thus represent the part missing from the reference orbit, i.e. the spherical harmonical spectrum from l ¼ 61 75. Anomalies from the LP75G model using l ¼ 61 75 have an rms of 21.50 mGal. Since spherical harmonics are global functions, there will be global information present in the data, that will be mapped into an area of limited extensions. In Fig. 4 the differences between the input and recovered anomalies are shown. This figure shows clearly the boundary effect. The standard way to deal with this is to leave out the boundary part, and thus retain only a central part. In Table 2 different sizes of such central areas are listed, together with the rms of differences. This shows that the initial rms of differences for the whole area of 12.88 mGal can be brought back to 1.01 mGal in case only the central part is taken along. So far, all solutions have been created without regularisation. However, the problem at hand is an illposed problem, mostly in the sense of stability, since
mgal -50 -40 -30 -20 -10
0
10 20 30 40 50
Table 2 Differences with the true solution for different central areas of the solution. Values are in mGal Area Input D~ g l ¼ 0–20 , l ¼ 3–17 , l ¼ 4–16 , l ¼ 5–15 ,
min/max/rms f ¼ 0–20 f ¼ 3–17 f ¼ 4–16 f ¼ 5–15
60.31/58.34/21.50 63.70/49.41/12.88 11.53/12.73/2.96 6.77/5.81/1.80 3.48/4.43/1.01
downward continuation of errors at satellite altitude will lead to amplified errors in the recovered anomalies on the reference surface. Thus, regularisation might be required to turn the problem into a well-posed one. Tikhonov regularisation is used, where the solution can be computed by means of: x^ a ¼ ðAT WA þ aIÞ 1 AT W~ y
(7)
with x^ a the estimated solution, W the weight matrix for the data, I the unity matrix and a a parameter that needs to be chosen. Other regularisation methods are also possible; especially worth mentioning is the use of a full variance–covariance matrix, as is applied in Barriot and Balmino (1992). Using a diagonal regularisation matrix assumes there are no correlations between the anomalies, while in reality, there are. If this is to be taken into account however, one will also have to prescribe a
ARTICLE IN PRESS 1338
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
certain structure of the signal, apart from the signal strength that is controlled when using a diagonal regularisation matrix. In the case of the Moon, prescribing this structure can be difficult, since the global models themselves rely heavily on the Kaula rule mentioned in Section 1. This is why here it is chosen to use only a diagonal matrix. Tests with a non-diagonal matrix have also been performed, where the anomalies are constrained by taking into account values of their direct neighbours, but these tests showed no improvement was to be gained by this. This was also found for instance by Ilk (1993) for the case of SST data for an Earth-orbiting satellite. The value for a needs to be determined, and in these simulations this can be done by trial-and-error, since the input truth anomalies are known. In reality, the truth is not known, but simulations like these can be used to estimate reasonable values for the regularisation parameter a. This results in a ¼ 10þ7 giving the smallest rms of differences with the input anomalies, namely 6.08 mGal, compared to the 12.88 mGal when no regularisation is used. These numbers are valid for the whole solution. In case the central part is examined, it turns out the regularised solution has an rms of differences of 1.58 mGal, compared to the 1.01 mGal of the unregularised solution. It should be noted that the regularisation matrix is also applied to the whole solution. In case the trial-and-error method is again applied, but now comparing the central part of the regularised solution with the central part of the input anomalies, it turns out that a ¼ 10þ6 leads to an rms of differences of 0.508 mGal, showing that indeed regularisation can help to improve the quality of the solution. Nevertheless, deriving a regularisation matrix which will only be applied to the central part can also be worthwhile to improve the solution even further. The results show that recovery down to mGal level is still possible, yet it should be noted this highly depends on the input signal as well. In this case, the new orbit was generated using l max ¼ 75, which means a resolution of roughly 180=75 ¼ 2:4 , which is larger than the 1 aimed for. Therefore, the numbers will be different in case there will be signal beyond the aimed-for resolution. Then one might also consider subdivision of the grid cells when forming the normal matrices as mentioned in Section 2. This has not been further investigated here. 3.5. Stochastic errors on the data Apart from possible systematic errors, such as farfield contributions, observation errors such as noise will obviously also be of influence on the recovery results. To simulate this, the data are contaminated by adding a random number from a Gaussian distribution times a certain standard deviation. In dealing with noise on the
-
-
Fig. 5. The rms of differences using all possible combinations of the 4 data sets versus the regularisation parameter when Tikhonov regularisation is applied. The average altitude of the satellite orbit is 50 km. The error signals used on the data are sr ¼ 0:5 m and sr_ ¼ 1:0 mm=s.
data coming from a random process, the more data are available, the better. For the simulations presented here, 4 data spans will be generated: one starting at June 10 1999 (990610), one at June 20 1999 (990620), one at July 10 1999 (990710) and one at July 20 1999 (990720). The noise on these data will have a realistic level: Konopliv et al. (2001) reported a data noise level of 0.5 m for range data and 1.0 mm/s for Doppler data, which are the levels used here. The recovery grid will again be the same as mentioned in Table 1, consisting of 11 11 anomalies. With such error levels, regularisation will be indispensable. Again, Tikhonov regularisation will be applied, with a free to choose. Fig. 5 shows the rms of differences with the input anomalies using all possible combinations of the data, for various values of the regularisation parameter a. This figure shows there actually is an optimal regularisation parameter, which is an important result when this method will be used on real Lunar Prospector data. The optimal parameter is a ¼ 10þ7 and produces an rms of differences of 13.52 mGal when all data sets are used. Compared to the rms of the input anomalies of 21.54 mGal, this result shows most of the signal can be recovered. The correlation of 0.78 between recovered anomalies and input anomalies confirms this. 3.6. Influence of altitude The results in the previous sub-section were all generated using an average satellite altitude of 50 km above the lunar surface. Even though this is already very low, Lunar Prospector flew at even lower altitudes
ARTICLE IN PRESS S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
1339
these data sets, the rms of differences without regularisation becomes 3.79 mGal, slightly better than the value of 3.86 mGal when using only the 20 km data. The most important result however is that the effects of downward continuation remain limited.
4. Conclusions
-
-
Fig. 6. Rms of differences between Tikhonov regularised solutions and the input anomalies versus the regularisation parameter a, for all possible combinations of the 4 data sets. The average altitude of the satellite is now 20 km, as opposed to 50 km in Fig. 5. The level of noise is the same.
during its Extended Mission. Therefore, the same simulations as described above have been conducted for an average altitude of 20 km. The data were again contaminated with noise with standard deviations sr ¼ 0:5 m and sr_ ¼ 1:0 mm=s, and regularised solutions were computed using different values for the regularisation parameter. The results are shown in Fig. 6. This figure again shows there is a minimum in the rms of differences between recovered and input anomalies, but what is more important is that the figure shows that when extremely low-altitude data are included, accurate solutions close to the optimum can be obtained without regularisation. This is because downward continuation becomes less and less of an issue with decreasing altitude. The best rms of differences is again obtained for a ¼ 10þ7 ; this is also an important result since the optimal regularisation parameter (in the sense of the smallest rms of differences) remains the same with varying altitude (it will however depend on the level of actual noise in the data). The best rms of differences is 3.60 mGal, which is already much better than the value of 13.52 mGal for the 50 km data. Without regularisation, the rms of differences is 3.86 mGal. This also shows that regularisation in the sense of including a priori information in the form of a regularisation matrix is not necessarily required in order to obtain a good solution at a resolution of 1 when low-altitude data are available. The improvement is not that much in terms of the rms of differences when using regularisation, and the cost would be to find the correct regularisation parameter. The combination of both the 50 km data and the 20 km data has also been investigated. When using both
A recovery method based on a linear variational approach, has been developed and thoroughly tested. This method allows gravity anomalies to be determined from Earth-based tracking data residuals. The method has been extensively tested by means of simulations. A short-arc approach has been used, where one arc or track consists of the period it takes the satellite to fly over the area of interest. It has been tried to estimate state-vector adjustments from such a short-arc, but this did not lead to satisfactory results. This means that prior to the determination of the anomalies, a precise orbit has to be determined, using elaborate reference models. In the presence of various error sources, recovery down to the level of several mGal is still possible. In the presence of noise on data taken at an average altitude of 50 km, regularisation will be indispensable to obtain a high-quality solution. An optimal regularisation parameter can be found, as is shown by comparing the recovered solutions with the input anomalies. The best rms of differences between input and recovered anomalies is 13.52 mGal, when using a regularisation parameter a ¼ 10þ7 . Compared to an input rms of 21.54 mGal, this shows most of the input signal can be recovered. When low altitude data are included, at an average of 20 km (which are available for the Lunar Prospector mission), regularisation is no longer required to obtain a good solution at a resolution of 1 , since the effects of downward continuation can be limited. The rms of differences without regularisation is 3.86 mGal, and with regularisation it is 3.60 mGal. This shows that by including low-altitude data, the recovery method becomes stable and hardly requires regularisation. These results show that the recovery method is capable of extracting the high-frequency gravity information from tracking data residuals. This method will therefore be used for processing real Lunar Prospector data, which will be the topic of a follow-on paper.
References Barriot, J., Balmino, G., 1992. Estimation of local planetary gravity fields using line of sight gravity data and an integral operator. Icarus 99 (1), 202–224. Barriot, J.-P., Balmino, G., Vale`s, N., 1997. Building reliable local models of the Venus gravity field from the Cycles 5 and 6 of the Magellan LOS gravity data. Geophysical Research Letters 24 (4), 477–480.
ARTICLE IN PRESS 1340
S. Goossens et al. / Planetary and Space Science 53 (2005) 1331–1340
Floberghagen, R., 2002. Lunar Gravimetry. vol. 273, Astrophysics and Space Science Library. Kluwer Academic Publishers, Dordrecht, The Netherlands, 304 pp. Heiskanen, W.A., Moritz, H., 1984. Physical Geodesy. Institute of Physical Geodesy, Graz, Austria (original edition by W.H. Freeman and Company, San Francisco, 1967). Ilk, K., 1993. Regularization for High Resolution Gravity Field Recovery by Future Satellite Techniques. In: Anger, G., Gorenflo, R., Jochmann, H., Moritz, H., Webers, W. (Eds.), Inverse Problems: Principles and Applications in Geophysics, Technology and Medicine. vol. 74 of Mathematical Research. Akademie, Berlin, pp. 189–214. Kahn, W., Klosko, S., Wells, W., 1982. Mean gravity anomalies from a combination of Apollo/ATS 6 and GEOS 3/ATS 6 SST tracking campaigns. Journal of Geophysical Research 87 (B4), 2904–2918. Kaula, W., 1996. Regional gravity fields on Venus from tracking of Magellan cycles 5 and 6. Journal of Geophysical Research 101 (E2), 4683–4690. Kaula, W.M., 1966. Theory of Satellite Geodesy, Applications of Satellites to Geodesy. Blaisdell Publishing Company. Konopliv, A., Binder, A., Hood, L., Kucinkas, A., Sjogren, W., Williams, J., 1998. Improved gravity field of the moon from lunar prospector. Science 281 (5382), 1476–1480. Konopliv, A., Asmar, S., Carranza, E., Sjogren, W.L., Yuan, D., 2001. Recent gravity models as a result of the lunar prospector mission. Icarus 150, 1–18. Lemoine, F., Smith, D., Zuber, M., Neumann, G., Rowlands, D., 1997. A 70th degree lunar gravity model (GLGM-2) from Clementine and other tracking data. Journal of Geophysical Research 102 (E7), 16,339–16,359. Montenbruck, O., Gill, E., 2000. Satellite Orbits. Springer, Heidelberg. Muller, P., Sjogren, W., 1968. Mascons: Lunar Mass Concentrations. Science 161, 680–684. Muller, P., Sjogren, W., Wollenhaupt, W., 1974. Lunar Gravity: Apollo 15 doppler radio tracking. The Moon 10, 195–205. Press, W., Teukolsky, S., Vetterling, W., Flannery, B., 1997. Numerical Recipes in Fortran 77: The Art of Scientific Computing. Cambridge University Press. 2nd ed. vol. 1. Fortran Numerical Recipes. Rapp, R., 1974. Procedures and Results related to the direct determination of gravity anomalies from satellite and terrestrial
gravity data. Reports of the Department of Geodetic Science, vol. 211, Ohio State University. Rowlands, D., Ray, R., Chinn, D., Lemoine, F., 2002. Short-arc analysis of intersatellite tracking data in a gravity mapping mission. Journal of Geodesy 76, 307–316. Rowlands, D., Luthcke, S., Klosko, S., Lemoine, F., Chinn, D., McCarthy, J., Cox, C., Anderson, O., 2005. Resolving mass flux at high spatial and temporal resolution using GRACE intersatellite measurements. Geophysical Research Letters 32 (L04310) doi:10.1029/2004GL021908. Rudolph, S., 2000. Regionale und globale Gravitationsfeldanalyse hochauflo¨sender Satellitendaten mittels mehrgitterverfahren. Ph.D. thesis, Rheinische Friedrich-Wilhelms-Universita¨t Bonn, Mu¨nchen, deutsche Geoda¨tische Kommission,Dissertationen, Reihe C, Heft Nr. 528. Sjogren, W., Muller, P., Wollenhaupt, W., 1972. Apollo 15 Gravity Analysis from the S-band transponder experiment. The Moon 4, 411–418. Sjogren, W., Wimberly, R., Wollenhaupt, W., 1974a. Lunar Gravity: Apollo 16. The Moon 11, 35–40. Sjogren, W., Wimberly, R., Wollenhaupt, W., 1974b. Lunar Gravity: Apollo 17. The Moon 11, 41–52. Sjogren, W., Wimberly, R., Wollenhaupt, W., 1974c. Lunar Gravity via the Apollo 15 and 16 Subsatellites. The Moon 9, 115–128. Sugano, T., Heki, K., 2004a. High resolution gravity anomaly map from the lunar prospector line-of-sight acceleration data. Earth, Planets and Space 56, 81–86. Sugano, T., Heki, K., 2004b. Isostasy of the Moon from highresolution gravity and topography data: Implication for its thermal history. Geophysical Research Letters 31 (24) doi: 10.1024/ 2004GL022059. Thornton, C., Border, J., 2000. Radiometric Tracking Techniques for Deep-Space Navigation. Monograph 1, Deep-Space Communications And Navigation Series, Jet Propulsion Laboratory, JPL Publication 00–11. Visser, P., 1992. The Use of Satellites in Gravity Field Determination and Model Adjustment. Ph.D. Thesis, Delft University of Technology. Vonbun, F., Kahn, W., Wells, W., Conrad, T., 1980. Determination of 5 5 gravity anomalies using satellite-to-satellite tracking between ATS-6 and Apollo. Geophys. J. R. Astr. Soc. 61, 645–657.