Access the full text.
Sign up today, get DeepDyve free for 14 days.
Astrodynamics Vol. 5, No. 3, 217–236, 2021 https://doi.org/10.1007/s42064-020-0098-1 SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 1 1 2 2 William Gordon Ledbetter (B), Rohan Sood , James Keane , and Jeffrey Stuart 1. Astrodynamics and Space Research Laboratory, The University of Alabama, Tuscaloosa, AL 35487, USA 2. Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA ABSTRACT KEYWORDS A growing interest in small body exploration has motivated research into the rapid char- gravimetry acterization of near-Earth objects to meet economic or scientific objectives. Specifically, asteroid exploration knowledge of the internal density structure can aid with target selection and enables astrodynamics an understanding of prehistoric planetary formation to be developed. To this end, stochastic optimization multi-layer extensions to the polyhedral gravity model are suggested, and an inversion Research Article technique is implemented to present their effectiveness. On-orbit gravity gradiometry is Received: 26 May 2020 simulated and employed in stochastic and deterministic algorithms, with results that Accepted: 26 October 2020 imply robustness in both cases. © The Author(s) 2020 1 Introduction for model-predictive control [1], or, as is the topic of this paper, obtaining information about the internal structure Exploratory missions to small, unknown celestial bodies of the body. are limited in their ability to adequately characterize and Previous studies considered the use of gravimetry with investigate the resource potential of a body of interest. passive craft [2], and existing methods on polyhedral Multiple near-Earth objects (NEOs) have been identi- density estimation require an existing gravity model of a fied as potential targets for human exploration with the body prior to any analysis [3,4]. Taking inspiration from goal of harnessing available resources. However, limited the GOCE mission [5] and work by Park et al. [6], in knowledge of the physical characteristics of NEOs, such as this investigation, sensor-equipped probes derived from their shape, density, gravity field, and composition, poses the work of Carroll and Faber [7] enable direct gravity a challenge to any manned exploration. Additionally, measurements. Additionally, the probes are tracked by both asteroids and comets are of significant importance the chief as they move in the vicinity of a central body. to answer fundamental questions regarding the origin Carroll and Faber proposed using a bias-free accelerom- and evolution of our solar system. Driven by the need eter at a fixed distance from the center of mass of the for scientific exploration and considering our economic spacecraft to measure the tidal acceleration on the space- interest, the proposed strategy utilizes a swarm of Small- craft. Taking into account that the tidal acceleration is Sats to address specific strategic knowledge gaps (SKGs), a function of the gravity gradient (GG) acting on the enabling resource investigation and acting as a precursor spacecraft, GG was chosen as the measurement type for to human exploration. The algorithm proposes a multi- the recovery. vehicle approach to significantly enhance the quality of data and assist in close-proximity guidance of space- Our procedure requires only three sets of data to obtain craft. Specifically, a swarm of small probes, deployed a density map: a shape model, probe positions, and the corresponding gravity gradient measurements at those from the primary spacecraft, performs flybys of the ce- positions. This study focuses on the development and lestial object, enabling high-fidelity, in-situ gravimetry initial testing of a noise-free recovery algorithm, and also measurements by the chief. These measurements are then used to construct a gravity model, which can be used undertakes a stochastic analysis of positional uncertainty. B wgledbetter@crimson.ua.edu 218 W. G. Ledbetter, R. Sood, J. Keane, et al. To simulate this problem, a dynamical model was con- polyhedral model is given in Eq. (1): structed that consists of a swarm of probes, which is U = Gρ r¯ · E · r¯ · L e e e e ejected from the chief, and which flies by the central eϵ edges body while collecting gravity gradient data. We then exe- − Gρ r¯ · F · r¯ · ω (1) f f f f cute the gravity inversion algorithm, as detailed below, fϵ facets to obtain densities for each “node” of the shape model where G is the gravitational constant, ρ is the average of the central body. density of the body, r¯ is a vector from the point in space to a point on edge e, E is a matrix describing 2 Dynamical model the geometry of edge e, and L is the potential of a 1D “wire” expressed in terms of the distances to the edge Modern research in asteroid gravimetry utilizes the endpoints. In the second term, r¯ is a vector from the polyhedral gravity model developed by Werner and query point in space to a point on facet f, F is a matrix Scheeres [8]. The model is attractive for two primary describing the geometry and orientation of facet f, and reasons: ease of implementation and benefits over other ω is the signed solid angle subtended by facet f, from models. First, the model takes advantage of a three- the perspective of the point in space. Equations (2)–(5) dimensional (3D) shape reconstruction, data that can be detail the calculation of parameters used in Eq. (1): approximated from Earth via a light curve [9] or radar [10] observations. The existence of reliable methods for shape A B E = nˆ nˆ + nˆ nˆ (2) e A B 12 21 model construction enables the implementation of the F = nˆ nˆ (3) f f f polyhedral model in multiple scenarios [11]. Shape mo- r + r + e e,1 e,2 12 dels can also be used in conjunction with imaging for L = log (4) r + r − e e,1 e,2 12 navigation algorithms [12, 13], further motivating the r¯ · (r¯ × r¯ ) 1 2 3 shape model construction as a primary mission objective. ω = 2 arctan r r r + r r¯ · r¯ + r r¯ · r¯ + r r¯ · r¯ 1 2 3 1 2 3 2 3 1 3 1 2 Second, polyhedral models guarantee convergence of the (5) gravity field in close proximity to irregular bodies. The infinite sum of the spherical harmonics diverges within In Eqs. (2) and (3), nˆ is equivalent to nˆ and nˆ , f A B a sphere circumscribing the body, creating a significant all of which describe a unit outward normal vector for invalid zone for irregularly shaped objects. Although face f, A, or B, respectively. In Eq. (2), nˆ is a unit methods exist for mitigating this issue [14], the polyhe- vector that is in the plane of face A, is orthogonal to dral formulation presents itself as an appropriate vehicle the edge between faces A and B, and points out of face for investigating physical characteristics such as the shape A. Likewise, the unit vector nˆ is in the plane of face and density of an asteroid. B, orthogonal to the edge between B and A, and points out of face B. With the safe assumption that the shape 2.1 Heterogeneous polyhedral gravity model is static in the chosen reference frame, these values model can be pre-computed and need not be recalculated for The goal of this study is to use gravimetry techniques every evaluation of Eq. (1). to reconstruct the internal density maps. These recon- A simple modification of Eq. (1) via the distributive structions are often known as inversions because they property enables the quantification of heterogeneous den- progress in the direction opposite to that in which the sity, ρ , where each facet and edge has a unique density. equation is typically used. This practice necessitates the U = G r¯ · E · r¯ · L ρ initial expression of the gravity of a polyhedron with e e e e e eϵ edges some arbitrary internal structure as though the density − G r¯ · F · r¯ · ω ρ (6) were already known. The formulation begins with the e f f f f fϵ facets commonly used polyhedral gravity model, as detailed below. Equation (6) can then be converted into a vector opera- Werner and Scheeres’ formulation for a homogeneous tion: SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 219 U = u¯ · ρ¯ = [u¯ , u¯ ] · [ρ¯ , ρ¯ ] (7) However, they assumed complete heterogeneity, thus sim- e f e f plifying the layering assumptions without sacrificing the The vector u¯ encapsulates both the asteroid shape and scope of the expressible internal structures. The construc- the satellite position via Eqs. (8) and (9), which are tion of radial variability begins by assuming that each extracted from Eq. (6): internal layer is a concentric, scaled copy of the outermost u¯ = G[L (r¯ · E · r¯ ), L (r¯ · E · r¯ ),··· ] (8) e e1 e1 e1 e1 e2 e2 e2 e2 polyhedron. This assumption preserves the E and F e f matrices from the reference shape. Therefore, only the u¯ = −G[ω (r¯ · F · r¯ ), ω (r¯ · F · r¯ ),··· ] f f1 f1 f1 f1 f2 f2 f2 f2 (9) relative position vectors r¯ and r¯ , and the scalars L and e f e ω vary with each subsequent layer. The mathematical Each ρ and ρ are evaluated at the associated node, e f procedure necessary for layering is explained in Fig. 1 which is the midpoint of the respective facet or edge. and Eq. (10). In Eqs. (10)–(12), the size of each color Our notation convention is as follows: u¯ is a vector of all wheel is analogous to u¯ and the orientation of the color u and u , whereas u¯ is a vector of all u and u¯ is a e f e e f pattern is analogous to ρ¯ . In this study, the outermost vector of all u . The density vectors ρ¯ , ρ¯ , and ρ¯ are f e f shape and density pattern are assigned i = 1, with an paired with their corresponding u¯ . This convention is () increasing value of i for the inner layers. consistent for forthcoming derivations of the acceleration and gravity gradient, with the dimension of shape-based U = (u¯ − u¯ )ρ¯ + (u¯ − u¯ )ρ¯ + u¯ ρ¯ = u¯ · ρ¯ (10) 1 2 1 2 3 2 3 3 quantities increasing by one for each derivative (u¯ , (1×n) u¯ = [(u¯ − u¯ ), (u¯ − u¯ ), u¯ ] (11) 1 2 2 3 3 a , g¯ g ). (3×n) (3×3 ×n) ρ¯ = [ρ¯ , ρ¯ , ρ¯ ] (12) Thus far, the dynamical model is capable of calculating 1 2 3 the gravitational potential of a polyhedron with non- Acceleration is analytically obtained by taking the uniform density by assuming that each facet- and edge- gradient of Eq. (6), shown in Eq. (13), and can be refor- defined tetrahedron of the shape has a unique density. mulated as represented by Eq. (14): Each tetrahedron consists of one facet of the shape model X X and the origin of the coordinate system. This formulation A = −G r¯ · E · L ρ + G r¯ · F · ω ρ e e e e f f f f eϵ edges fϵ facets assumes variability only in the latitude and longitude (13) dimensions. The density of a facet represents the density A = a · ρ¯ = [a ,a ] · [ρ¯ , ρ¯ ] (14) of the volume of the tetrahedron, from facet to origin. e f e f In the derivation of Werner and Scheeres’ equation, each Here, u¯ is a n × 1 row vector, whereas a is a 3 × n edge must be counted twice, once for each facet it borders. matrix. As a result, A is a 3-element vector. The first In their formulation of the constant-density gravity field, step of the acceleration calculation entails building the a common terms cancel and combine such that each edge is matrix according to Eqs. (15) and (16): considered only once. For a heterogeneous body, the edge is still considered once, but with a density that is the a ¯ = −GL (r¯ · E ) (15) e e e e average of that of the two adjacent facets. As a result, a ¯ = Gω (r¯ · F ) (16) f f f f each tetrahedron is homogeneous from the center of the To maximize the computational efficiency, each a ¯ and body to the assumed surface. Radial variability is achieved by layering multiple poly- a ¯ vector is calculated in parallel and stored in the a hedra. In previous work, non-concentric homogeneous matrix, as is the case for u and u above. Considering e f layering was utilized to represent various morphologies [4]. a multi-layered body, the a matrices for each layer need Fig. 1 Visualization of polyhedral layering. The wheel size is analogous to the asteroid shape u¯ , and the color pattern is analogous to the density ρ¯ . By assuming concentric layering, the angular coordinates of vertices are consistent across layers. i 220 W. G. Ledbetter, R. Sood, J. Keane, et al. to be subtracted, as shown in Fig. 1 and Eq. (17): where X is a gravity measurement, ρ¯ is the density vector of the body, and x is the matrix function that relates A = (a − a )ρ¯ + (a − a )ρ¯ + ··· 1 2 1 2 3 2 the two vectors. In this case, the x matrix is a function + (a − a )ρ + a ρ¯ (17) n−1 n n−1 n n of shape and position, as explained in Eqs. (7)–(9). The a = [(a − a ), (a − a ),··· , (a − a ),a ] (18) recovery procedure was chosen to process the gravity gra- 1 2 2 3 n−1 n n dient as measurements to mimic the GOCE mission [5]. The acceleration (A) is the gradient of potential, whereas the gravity gradient (GG) is the gradient of 3.2 Observability and correlation acceleration, and can be represented by Eq. (19) [8]: X X In Ref. [6], an analysis is presented regarding finite-cube GG = −G E · L ρ + G F · ω ρ (19) e e e f f f and finite-sphere gravity models wherein the observability eϵ edges fϵ facets of internal density was estimated from on-orbit data. Again, the density is extracted via vectorization, resulting The method utilized a batch-least-squares filter that was in the g¯ g term, which is a vector of 3 × 3 matrices. applied to the range and range-rate measurements; no direct gravity measurements were considered. In recent GG = g¯ g · ρ¯ = [g¯ g ,g¯ g ] · [ρ¯ , ρ¯ ] (20) e f e f years, innovations in sensor technology have opened the gg = −GL E (21) e e e possibility of measuring the GG directly from small-scale gg = Gω F (22) f f f spacecraft [7]. The same mathematical technique used by Park et al. [6] can now be applied to predict the Each 3 × 3 GG matrix is then organized into a 6-element effectiveness of this new gravimetry technique. vector such that gg is a 6 × n matrix. Below-diagonal Assume that a gravity gradient measurement is taken components are omitted because the GG matrix is sym- metric: at a point in space, and the measurement can be defined by a polyhedral gravity model. The equation that would GG GG GG xx xy xz predict the measurement is GG GG GG yx yy yz GG GG GG zx zy zz ′ ′ ′ GG = gg · ρ¯ = [gg ,gg ] · [ρ¯ , ρ¯ ] (27) e f e f → − [GG GG GG GG GG GG ] (23) xx xy yy xz yz zz ′ Recall that GG is the vectorized gravity gradient ma- gg→ − g¯ g (24) trix, gg is a 6× n matrix of unwrapped gg matrices (e/f) The procedure for layering is identical to the potential for each edge and facet, and ρ¯ is the density vector descri- and acceleration: bing the body. Each matrix in the edge and facet vector ′ ′ ′ ′ ′ ′ is given in Eqs. (21) and (22), and is restructured accord- gg = [(gg − gg ), (gg − gg ),··· ,gg ] (25) 1 2 2 3 n ing to Eq. (23). The recovery problem assumes that the Because the gravitational potential must satisfy Laplace’s position and gravity gradient can be measured. The po- equation, only five of the six terms are independent. ′ sition is used to calculate gg , the gravity measurements Nonetheless, in consideration of future stochastic analy- are substituted for GG, and the goal of the problem is ses, all six components are simulated. to find the density vector, ρ¯ that most closely satisfies Eq. (27). 3 Procedure For a given position, the relationship between the den- sity and gravity gradient is linear. Subsequently, partial 3.1 Theory: mathematical recovery derivatives with respect to the density are found simply This section presents the development of a method that to be takes a known shape model and gravity measurements as ∂GG = gg (28) inputs, and outputs a density vector that, when paired ∂ρ¯ with the shape model, reproduces the gravity measure- Then, following Park et al., define ments. This relationship is expressed via matrix multi- ′T ′ plication: Λ = gg gg (29) X = xρ¯ (26) to obtain the covariance matrix: SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 221 −1 P = Λ (30) Algorithm 1 Adam algorithm (reproduced for reference) The covariance matrix given in Eq. (30) provides the cor- Input: α : Step size relations and standard deviations of elements in the den- Input: β , β ∈ [0, 1): Exponential decay rates 1 2 sity vector for a given set of measurement data. Analysis Input: i : Maximum iterations max of these values informs predictions about the accuracy Input: ∇f(x¯ ): Gradient function of a recovery attempt. For example, a high standard Input: x¯ : Initial guess 1: i ← 0 deviation indicates high uncertainty about the value of 2: m ¯ ← 0 a density node, and a high correlation with other nodes 3: v¯ ← 0 implies that the gravitational effects caused by that node 4: while i < i do max are difficult to distinguish from the effects of other nodes. 5: i ← i + 1 6: g¯ ← ∇f (x¯ ) i i i−1 3.3 Adam optimizer 7: m ¯ ← β m ¯ + (1 − β )g¯ i 1 i−1 1 i [2] 8: v¯ ← β v¯ + (1 − β )g¯ (square bracket i 2 i−1 2 The problem of high-resolution polyhedral density re- i exponent indicates element-wise power) construction bears a certain rudimentary resemblance 9: α ← α · 1 − β /(1 − β ) (correct for to the problem of training a neural network: a large initialization bias) number of parameters must be tuned to recreate a set of 10: x¯ ← x¯ − α m ¯ /( v¯ + ε) (element-wise i i−1 i i i noisy data. Methods that have delivered the best perfor- square root and division) mance in high-dimensional neural network training can 11: end while 12: return x¯ be classified as stochastic gradient methods. One such i algorithm, Adam, has quickly become one of the most popular methods since its introduction in 2014 [15]. By written in the same C++ class as the main optimizer, keeping the running average of the second moment of the and is set as the default mode. The objective of the gradient, Adam is able to scale the step size of each itera- least-squares optimization is tive cycle for empirically improved convergence. Certain min(y¯ − Ax¯ ) (31) situations result in slower convergence [16], and prior work has sought to modify the base algorithm to accom- where x¯ is the design vector, y¯ is the data vector, and A modate such cases [17]. In this investigation, the original relates the two, such that, for the optimal case, y¯ = Ax¯ . formulation based on adaptive moment estimation is im- Adam operates using only the gradient, which is plemented. The core loop, Algorithm 1, is shown for T T g¯ (x¯ ) = 2A Ax¯ − 2A y¯ (32) reference. Note that the loop control parameter, i , max limits the procedure to a fixed number of iterative cycles. Batching is often utilized when training deep neural This allows the user to terminate the computation and networks. In its most common form, batching consists investigate the initial behavior of the algorithm before of choosing a subset of data points to calculate the gra- investing significant time into a convergence-evaluated dient function, performing a few iterative cycles and test case. then randomly re-selecting those data points. This tech- The stochastic gradient function, ∇f, is minimized nique can prevent the optimizer from being trapped in by using the running averages of the mean, m ¯ , and the the non-convex regions of the solution space. For the uncentered variance, v¯ , to calculate an appropriate step least-squares formulation, the data points are re-sampled size. Kingma and Ba draw an analogy between the to select rows of the A matrix and the corresponding ratio m ¯ / v¯ and the signal-to-noise ratio (SNR) [15]. elements in the y¯ vector, and substituting these trun- i i A smaller SNR (larger v¯ ) indicates greater uncertainty cated values in Eq. (32). Two parameters that are often in the direction of the true gradient, and the step size is used to control batching are the batch size and batch appropriately scaled in line 10 of Algorithm 1. frequency. Each specification is fairly intuitive: the batch The optimization problem in this study is formulated size refers to the number of points extracted from matrix as a least-squares minimization. Because of its broad A and vector y¯ in each sample, and the batch frequency applicability, the framework for solving such problems is determines how often the data are resampled. 222 W. G. Ledbetter, R. Sood, J. Keane, et al. Table 1 Vertex file stores the coordinates of each vertex 3.4 Implementation X Y 3.4.1 Gravity model 1 1 1 A strategic implementation process was structured 2 2 3 3 6 2 around a sequence of events similar to that of an ac- 4 7 4 tual mission. Based on an object-oriented programming paradigm, a C++ class that contains the necessary proce- The corresponding facet file is provided in Table 2. dures for gravity calculations and recovery was developed to assist with investigations. First, the target is defined Table 2 Facet file stores the indices of the vertices ( V , V , 1 2 by constructing and importing a polyhedral model. A and V ) composing each facet polyhedron is fully defined by two sets of data: vertices V V V 1 2 3 and facets. Moreover, two derivative datasets are used 1 1 3 2 to provide information that is relevant for Eq. (1): edges 2 2 3 4 and facets-of-edge. Vertices are simply the indexed X– Y –Z coordinates of points on the surface of the body, The edge file (Table 3) defines all connected vertex with dimension numVert × 3 in three-dimensional space. pairs. Each facet is a set of three vertex indices that form a Table 3 Edge file stores the indices of the vertices ( V and counterclockwise triangle when viewed from outside the V ) comprising each edge body, with dimension numFacet × 3. The edge file is V V 1 2 a list of all pairs of points that are connected in the 1 1 2 facet definition, with dimensions of numEdge × 2. Each 2 1 3 edge borders two facets, and the facets-of-edge file defines 3 2 3 both the index of the adjacent facets and the associated 4 2 4 5 3 4 clockwise (CW) or counter-clockwise (CCW) direction of the edge relative to each facet. This file also has the dimension numEdge × 2. Finally, the facets-of-edge file (Table 4) defines which Considering the two-dimensional (2D) shape in Fig. 2, two facets border each edge, along with their respective the coordinates of points 1–4 are defined in the vertex orientations (with x representing a facet not shown in file in Table 1. Fig. 2). Fig. 2 2D polyhedron example. Planar points P1–P4 are used to define five edges and two facets. Tables 1–4 are used to define this shape in the C++ implementation. SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 223 Table 4 Facets-of-edge file stores the two facets adjacent to Recall the assumption that each interior layer is a scaled every edge, and their CCW or CW relationship copy of the outermost; thus, the λ and φ angles are CCW CW consistent for all layers. The discretization of λ and φ is 1 x 1 a function of the shape model and may not be uniformly 2 1 x distributed. 3 2 1 4 x 2 3.4.2 Matrix inversion recovery 5 2 x The GG and position data must now be recovered into a density vector. This procedure requires access to the Each text file is imported into a matrix and used to shape definition of the central body. The core of the initialize an object of the C++ class. The number of recovery is the solution to the linear algebra equation: layers is simply passed to the object as an integer. After importing the shape model, a density map is Ax¯ = b defined for each layer. This task is accomplished by where x¯ is density, and b is the GG measurement vector. using spherical harmonics to create deviations from the m m The number of unknowns is the product of the number average density. A set of Stokes coefficients ( C , S ) is l l of layers with the sum of edges and facets. passed to the object for each layer of density, enabling different patterns of variation for each layer. Within First, the A matrix is constructed. Recall Eq. (27), each layer, the density is only a function of the latitude where density extraction results in a vector of matrices ¯ ¯ and longitude. The radial variation must be defined by gg. By unwrapping each matrix in gg (Eqs. (21) and selecting different patterns for each layer. The use of (22)), as shown by Eq. (23), an n-dimensional vector of the spherical-harmonic-defined density enables the same symmetric 3 × 3 matrices becomes one 6 × n matrix. density map to be tested on different shape models. In The GG equation now resembles the acceleration formula addition, this method enforces reasonable continuity and with a height of 6 instead of 3: maintenance of the average. An option to explicitly dene fi GG = gg · ρ¯ (34) the density of each node via the vector input also exists. Subsequently, the object undergoes an initialization Note that the gg matrix is only a function of the position. process that calculates and stores certain parameters Attempting to invert Eq. (34) with a single sample point that are required for later calculations. These parameters is an underconstrained problem, given that the num- include a point in each facet (PIF), a point on each edge ber of density nodes, even for a polyhedron of medium (POE) as well as F , E , µ , and, if the spherical harmonic f e resolution, is much greater than 6. However, by accumu- density is used, the explicit density of each node. POE lating multiple gravity samples into a larger matrix, the is the midpoint of each edge, and PIF is the centroid of problem can be sufficiently constrained. For s samples, each facet. The matrices E and F are evaluated by e f we have GG and gg . Equation (35) outlines (6s×1) (6s×n) employing Eqs. (2) and (3), respectively. The explicit the fundamental requirements for the necessary number density of each node is calculated using the spherical of sample points and the number of probes required to coordinates of that node via PIF and POE. The latitude obtain those points: (φ) and longitude (λ ) of each PIF and POE are input s = p · k ⩾ n/6 (35) into Eq. (33) along with the appropriate set of density coefficients. where s is the total number of samples, p is the number of ∞ l probes, k is the number of samples taken on each probe, X X ρ (φ, λ) = ρ P (sin φ) i i,avg and n is the number of density nodes of the body. As l=0 m=0 is apparent from Eq. (35), as the number of probes in- m m · [C cos(mλ) + S sin(mλ)] (33) l,i l,i creases, a surplus of sample points is generated. Although m m ′ the option to use all the points does exist, making the gg The parameters ρ , C , and S are defined by the i,avg l,i l,i matrix very tall, full inclusion may significantly increase user to construct the body they wish to recover, and may be different for each layer. In this context, spherical the computational time necessary to obtain the solution. harmonics are used as a scalar function on a sphere. Therefore, a parameter termed the overconstraint factor 224 W. G. Ledbetter, R. Sood, J. Keane, et al. (OCF), which defines the height-to-width ratio of the rithm. Each body represents a different level of resolution matrix, is introduced. If p is the number of density coef- and, therefore, a different level of computational com- ficients, and c is the overconstraint factor, the number of plexity. Additionally, recall that the addition of layers sample points, N, in the inversion is defined by causes a multiplicative increase in the dimensionality of the problem. In the following analysis, the bodies are N = cp (36) used in the order of increasing fidelity. This choice allows The sample points are sorted by increasing the orbital the performance of the algorithm to be observed on rela- radius, and are sequentially pulled into the inversion tively fast tests before committing significant resources matrix until the overconstraint factor is met. Therefore, for more intensive analyses. Both bodies are assumed to the most distant points may be unused. All data points be non-rotating. in the matrix are assigned uniform weights. The linear A range of potential scenarios were simulated by vary- algebra problem is solved by executing a QR factorization ing certain parameters listed in Table 5. The number algorithm. of probes is varied to investigate the effects of ground Further analysis is then performed on the calculated track coverage and can directly affect the total number of density solution. One immediate test, performed as a sample points available in the inversion. The layers of all rough heuristic, was used to evaluate the validity of polyhedra are evenly spaced in the radial distance, and the density distribution. The C++ object is duplicated by considering more layers, the recovery procedure may to preserve all settings, and the output of the density return a more accurate representation of the central body. algorithm is set as its density map. Then, the potential Additionally, the number of layers affects the dimension- of the original and duplicate objects are evaluated at the ality of the inversion problem. The overconstraint factor same point, and each is displayed on the screen for user affects the number of samples that are accumulated in interpretation. If the inversion is successful, the outputs Table 5 Range of parameter variation will be identical. Although the test is not comprehensive, Property Range it provides a reasonable starting point for investigating unexpected outputs. Probe 1–25 Layer 1–3 Overconstraint factor 1–8 4 Results ◦ ◦ Ejection cone half-angle 22 –39 Batch size 8–8000 The two shape models, an octahedron and a disturbed Noise magnitude 0–10 m spheroid, shown in Fig. 3, were used to validate the algo- (a) Octahedral body: 20 nodes per layer, 10 km semimajor (b) Disturbed spheroid body: 560 nodes per layer, 15 km −4 3 2 −3 3 2 axis, µ ≈ 2.37 × 10 km /(kg·s ) semimajor axis, µ ≈ 2.35 × 10 km /(kg·s ) Fig. 3 Polyhedral bodies used for validation. SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 225 the information matrix. to the estimated density of asteroid Eros [18]. Recall The inertial velocity of the probes is then calculated that spherical harmonic coefficients are used to define such that the vectors are evenly spaced around a cone density variation. To obtain a valid density map, the of which the centerline connects the chief to the central first 6 × 6 coefficients from the 50 × 50 degrees and body. The half-angle of the ejection cone is conservatively order GEM-T3 Earth gravity model are utilized. In their chosen to prevent the probes from impacting the central normal context, these coefficients are used to define small body, and is proportional to the maximum radius of the gravitational perturbations; therefore, this set of coeffi- body. The resultant minimum altitudes range from 2 to cients also results in bounded variation from the mean 8 km, which is reasonable, because the probes are not density value. For the following tests, the sample rate recovered in this analysis. of the probes was arbitrarily selected to be 0.1 Hz. This Positional noise is assumed to follow a zero-mean Gaus- parameter directly affects the number of sample points sian distribution, the magnitude of which is controlled available for a given simulation duration. After manual by a user-defined standard deviation. All stochastic sce- testing, the ejection velocity was selected to approximate narios in the Results section have a positional standard the escape velocity. This choice, in conjunction with the deviation of 10 m. In reality, the positional noise depends cone half-angle, prevents the probes from impacting the on the navigation solution of the mission and can follow body for the range of explored parameters. non-standard distribution curves. For the purposes of Table 6 Constant parameters: the mean density [18], and this study, the Gaussian distribution is a reasonable test the density Stokes coefficients [19, 20] case because the use thereof is consistently able to cause Property Value Unit the matrix inversion technique to fail. 12 3 Mean density 2.67 ×10 kg/km The gravity gradient measurements are assumed to Density Stokes’ coeffs. GEM-T3 Earth model N/A originate from a system such as the tidal acceleration Sample rate 0.1 Hz gravity gradiometer proposed by Carroll and Faber [7]. Chief position (50, 0, 0) km Ejection velocity 1.4 µ/r km/s chief Unlike existing inversion techniques, this approach relies on on-board sensors as well as remote position determina- A typical simulation for a multiple-probe flyby of a tion. Additionally, even though we consider the position sample asteroid is illustrated in Fig. 4. In this case, five uncertainty directly, we currently assume perfect GG measurements. Recalling the matrix formulation of the probes are ejected from a relative position of (50, 0, 0) km from the asteroid, and are propagated forward for inversion, GG = gg · ρ¯ , the uncertainty in the position 6 h. The integration scheme is an Adams–Bashforth 10th perturbs the matrix gg, whereas the uncertainty in the order method, and the asteroid is assumed to be static. gravity measurement affects the vector GG. The values in Table 6 are held constant for all simu- In reality, an asteroid undergoes non-negligible rotation, lations. The average density was selected to correspond but the static assumption helps to verify the algorithm Fig. 4 Sample simulation using an Eros shape model. Five probes are ejected in a conical pattern and fly by the heterogeneous model of Eros. 226 W. G. Ledbetter, R. Sood, J. Keane, et al. and to isolate recovery behavior. The extension to the densities were accurate to machine precision for all 20 rotating case simply consists of a rotation transformation nodes in the body. Although this preliminary investiga- on the input data. The initial velocities for the probes tion instills confidence in the concept and validity of the are calculated to be evenly spaced around the body, implementation, higher fidelity modeling is required to resulting in more complete ground track coverage. The test the robustness of the strategy. Additionally, visual green and red trajectories are “behind” the body from representations of the density are difficult to interpret the given perspective. Owing to the inherent dynamical for such a low-resolution shape. nonlinearities, the probes do not maintain perfect spacing, The theoretical correlations for this case are given as evidenced by their non-singular reconvergence on the in Fig. 5. For each axis, node ID refers to the index right side of Fig. 4. However, they maintain sufficient of that node in the density vector, ρ¯ . Note that all spacing to achieve comprehensive ground track coverage, diagonal components were set to zero to ensure that the as was intended. color scale is more intuitive. The most evident difference Before noise is applied to the data, the theoretical ob- is the emergence of a visual distinction between edges servability is analyzed using Park’s method, as detailed and facets. Figure 5(b) shows consistent low-magnitude in Section 3.2. The standard deviations and the corre- negative correlations of edges with facets, and relatively lation matrix are saved in full precision and visualized small positive correlations between facets. The outliers later using applicable tools (MATLAB). Noise is sampled within the facets, such as (13, 18), are likely caused by from a Gaussian distribution with zero mean, where the the symmetry of the body. level of noise is defined by a standard deviation. A set The average standard deviation for the 10-probe case 5 3 is sampled from the defined distribution and is added to is 5.57 × 10 g/cm , an improvement of one order of nominal trajectories. magnitude over the average of the one-probe case of For Adam to operate, the problem must be cast as a 6 3 4.64 × 10 g/cm . However, because the average density least-squares matrix formulation. Thus, the objective is 3 of the body is 2.67 g/cm , both deviations suggest that ′ 2 a recovery would be imprecise. min(GG − gg · ρ)¯ (37) 4.1.2 Adding positional noise where GG is the measurement vector, ρ¯ is the den- Preliminary investigation suggested that optimization sity vector, and gg is the system matrix defined from with noise in the measurements may be nontrivial. How- Eq. (27). The noise associated with the position is con- ever, adding up to 10 m in positional uncertainty resulted sidered when calculating gg , and the noise in the gravity in convergence profiles that are nearly indistinguishable measurements is included in GG. Adam also requires an from the noise-free case. Zooming in to the level of Fig. 6 initial guess to start the procedure. In most cases, the is required to discern the difference. Note the scaling average body density is given as the initial guess. This of the y-axis. It is observed that in the 10 cm and 1 m assumption is reasonable because the average density of cases, the optimization outperformed that of the clean small bodies can often be estimated from the Earth. case. However, this result may be misleading: if the The optimization is initialized using the default values optimization is continued, the noisy cases will overshoot for exponential decay β = 0.9 and β = 0.999 [15]. 1 2 the target and settle with a higher final error. The step size, α = 20, is empirically evaluated, and is problem-specific (Kingma and Ba suggested an α of 4.1.3 Parameter tuning 0.001 for their application). Performance is quantified The concept of batching, as discussed in the section on the by calculating the L2-norm of the difference between the Adam Optimizer, is implemented in all simulations, and optimized density and the true density at each iterative two variables are introduced that control the process: the cycle. batch size and batch frequency. In the previous analyses, a large number of samples were available for use in the 4.1 Simulation 1: octahedron, 1 layer optimization, but because the batch size was constant 4.1.1 Noise-free case throughout, each iteration only utilized 800 samples. An In case one, a single probe flyby of a one-layer octahe- increase in the batch size enables more samples to be dron was performed. In this scenario, the reconstructed used in each iteration. The effects of varying the batch SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 227 (a) Case 1: 1 probe correlation matrix. No discernible (b) Case 1: 10 probe correlation matrix. Clear distinction patterns between edges and facets Fig. 5 Noise-free correlation comparison for case 1. The edges are top-left, and the facets are bottom-right, marked by bold lines. Symmetry is an intrinsic property of correlation matrices. Values on the diagonal are always equal to 1 by definition in all correlation matrices, but for intuitive coloration, these are plotted as zero. Fig. 6 Case 1: nal fi iterations of noise comparison. The effect of adding up to 10 m of error on the accuracy of the Adam algorithm after 50,000 iterative cycles is nearly negligible. size are illustrated in Fig. 7. The drastic improvement surprising. Despite the standard deviation being on the 5 3 in the convergence speed from a larger batch size can order of 10 g/cm , the densities are recovered to an −4 3 be explained by Kingma and Ba’s SNR analogy: a more accuracy of approximately 6 × 10 g/cm , a disparity comprehensive sample of the data results in a consistent of approximately 9 orders of magnitude. gradient vector. 4.2 Simulation 2: disturbed spheroid, 1 If the gradient is more stable, the algorithm takes larger layer steps. However, certain trade-offs exist with respect to 4.2.1 Noise-free case the runtime. Increasing the batch size is the same as increasing the height of the A matrix in Eq. (32); as such, To further continue the analysis and increase the com- the linear algebra operations may be computationally plexity, a “disturbed spheroid” replaces the octahedron expensive. as the central body. The spheroid has 560 density nodes Considering the standard deviation predictions from per layer, 28 times the resolution of the octahedron. The above, the result for the 8000 batch size is particularly simulations start with the same setup as the octahedron 228 W. G. Ledbetter, R. Sood, J. Keane, et al. Fig. 7 Case 1: batch size comparison. The inclusion of more data points per iteration results in more rapid convergence at the cost of robustness against non-convexity. —one layer and one probe. The increased resolution and because Fig. 8(c) is dominated by an incongruous spike geometric irregularities drastically impact the reconstruc- around (− 40, 10). The precise location of the spike does tion accuracy. A 3D representation of the orbit of the not seem to correspond with any particular feature of probe around the central body is represented in Fig. 8(a). the polyhedron, but it is notable that it occurs on the Figure 8(b) shows a lat/lon map of the true density hemisphere with no ground track coverage. Because of for one layer, with the ground track of the probe super- the apparent outlier in the reconstruction, the percent imposed in red. In the 3D graph, probes fly according error plot, Fig. 8(d), is shown on a common logarithmic to the arrows, and in lat/lon plots, probes start at the scale. A top–down visual inspection, Fig. 8(d), reveals an center and fly outwards. Recall that the ejection velocity interesting phenomenon with respect to the ground track is approximately equal to the escape velocity. Therefore, of the probe. The lowest errors occur directly below the for the short duration of these simulations, orbital cap- closest pass of the probe. Inspecting the error magnitude ture is not a concern. The minimum orbital radius for all via the color bar in Fig. 8(d) reveals that the error is disturbed spheroid simulations is approximately 23.7 km. too large for these data to be reliable. Three parameters Figures 8(c) and 8(d) have the same X (latitude) and Y may affect the error: the number of probes, data sample (longitude) axes as Fig. 8(b), but with different Z-axis rate, and overconstraint factor. First, we investigated the data (density or error). Figure 8(c) is the recovered den- effects of improved ground-track coverage by increasing sity, and Fig. 8(d) shows the error between the recovered the number of probes. and actual density. This collection of graphs forms the The correlation analysis can also provide insight into standard layout for subsequent analyses as well. The the poor accuracy of this recovery attempt and predict color conventions for the density plots are as follows: red any improvements that may arise by increasing the num- points are midpoints of edges, black points are midpoints ber of probes. Analysis of theoretical observability using of facets, blue points are vertices, yellow surface is for the same technique as case 1 produced a less intuitive higher z-value, and blue surface is for lower z-value. Some output. The standard deviation and correlation coeffi- points are unfortunately obscured by the 3D topography cient for many nodes was returned as NaN, implying that in the MATLAB plot; however, the colored interpolation these areas cannot be observed at all. These NaN nodes surface is an important analysis tool. Specific points are are illustrated as dark lines in Fig. 9(a). This unobserv- given to demonstrate the resolution of the body. ability is not surprising, considering the results in Fig. 8, With these details in mind, the oscillating pattern of but the complete absence of any apparent patterns was the true density generated by the spherical harmonic unexpected. definition is apparent in Fig. 8(b). These variations are However, for the 10-probe case in Fig. 9(b), the body is not immediately evident in the reconstruction, although, fully observable with correlations that generally approx- SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 229 (a) Orbit (b) Actual density (c) Recovered density (d) Percent density error Fig. 8 Simulation: disturbed spheroid, 1 probe, 1 layer. Color conventions: red points = edges, black points = facets, blue points = vertices, yellow surface = higher z-value, and blue surface = lower z-value. imate zero, predicting superior recovery performance. This graphical artifact may be misleading because an The average standard deviation for the 10-probe case is error still exists even though it was not captured owing 10 3 6.46 ×10 g/cm , which is significantly larger than the to truncation. mean of 2.67 g/cm ; however, previous results suggest 4.2.2 Adding positional noise that recovery may still be possible. Given the promising performance of the matrix inversion Keeping all other variables the same, the number of technique in the noise-free case, the Adam optimizer probes used in the recovery algorithm is increased to ten. is now employed in the same scenario with positional With the increase in ground coverage, the simulation uncertainty. A batch size of 8000 was selected (based on generates extremely promising results. The oscillating the results of Fig. 7), and 10 m is used as the standard pattern is clearly visible in the reconstruction, Fig. 10(a), deviation of the positional uncertainty. In Fig. 11, the and a side-view visual comparison between the truth and error is consistently reduced from the initial guess and result reveals the accurate reproduction of small graphical is unaffected by the number of probes in the simulation. nuances. The error plot, Fig. 10(c), confirms the visual This optimization proceeded significantly slower than the intuition, with large central zones vanishing completely. octahedral case, with 50,000 iterative cycles resulting in 230 W. G. Ledbetter, R. Sood, J. Keane, et al. (a) Case 2: 1 probe correlation matrix. Dark regions in- (b) Case 2: 10 probe correlation matrix. The near-zero correla- dicate unobservable regions. The absence of any patterns tions when using 10 probes indicate that each node is distin- indicates that accurate recovery is unlikely guishable from the others Fig. 9 Noise-free correlation comparison for case 2. Taken with Fig. 8, the improvement of the 10-probe case constitutes strong evidence of a relationship between the ground track coverage and recovery accuracy. −5 3 an approximately 3 × 10 g/cm improvement in the sion of the disturbed spheroid was set as the central density estimate. body. Density maps of the internal layers follow the Because these analyses do not converge, a different same pattern as the outermost layer, but with a reduced simulation is necessary to understand the convergence magnitude. The amplitude reduction is defined by behavior. Instead of using the average density as the l−1 [C , S ] = [C , S ]/5 (38) l l 1 1 initial guess, the correct answer (true density) is passed to the optimizer. The accuracy and convergence behavior where C and S are the Stokes coefficients, and l is the are interpreted from Fig. 12 based on the extent to which layer, with layer l = 1 being the outermost layer. Recall the algorithm is able to successfully maintain the solution, that the first 6 × 6 GEM-T3 Earth gravity model co- given the noisy data. efficients are used to define the outer-layer densities for It is evident that the number of probes has a stronger these simulations. Inspecting the error as a function of effect on the tail end of the optimization than at the layers, the outermost layer was found to have reduced outset. With the data from 20 probes, the Adam opti- in accuracy by two orders of magnitude relative to the −6 mizer is able to maintain the solution to within 1 × 10 single-layer body. The middle and innermost layers show g/cm , whereas the 5-probe simulation diverges much a further increase in error beyond the acceptable range. more rapidly. Noteworthy is that for any batch size less The addition of five more probes incrementally improved than the full dataset, the algorithm is not guaranteed the reconstruction for the outer and middle layers; how- to converge in the deterministic sense. Each iteration is ever, the central layer is largely unaffected even in a based on a different random subset of the available data, scenario with 20 probes. which will cause the stochastic directional changes seen Recall that the overconstraint factor (OCF) controls in Fig. 12. Because of these considerations, a limit in the the aspect ratio of the inversion matrix, and therefore number of iterative cycles is suggested as the stopping affects the number of sample points considered in the criterion instead of the step size tolerance. analysis. By increasing the OCF from 4 to 8, the number of utilized sample points doubles. This modification is 4.3 Simulation 3: disturbed spheroid, 3 lay- ers expected to have the same effect as increasing Adam’s 4.3.1 Noise-free case batch size parameter. For both the fifteen- and twenty- probe cases, increasing the overconstraint factor to 8 Next, multi-layer reconstruction of the central body was considered. Beginning with ten probes, a 3-layer ver- results in an error reduction of approximately an order SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 231 (a) Recovered density. The correct oscillating pattern was repro- (b) Recovered density: side view. The profile of the peaks and duced valleys matches that of the true density map (c) Percent density error. The central swath was recovered to (d) Percent density error: side view. Log errors below zero within operating precision indicate single-digit percent accuracy Fig. 10 Simulation: disturbed spheroid, 10 probes, 1 layer. First successful recovery of complex, irregular bodies. Table 7 Average magnitude of percent error for 3-layer of magnitude in layer 3. A full comparison of the average disturbed spheroid. The combination of additional probes percent error is given in Table 7, which demonstrates the with a greater number of considered data points, via the layer-3 average error reduction from 444% to 5.28%. OCF, results in a two-order-of-magnitude improvement in the 4.3.2 Adding positional noise central layer recovery accuracy In the stochastic analysis, the correlation matrix again Probes Layer OCF: 4 OCF: 8 predicts large areas that are not observable. For the 1 0.3009% 0.00094% areas that returned a number, the average standard devi- 15 2 16.6801% 0.2164% 14 3 3 443.9877% 85.1575% ation was 2.18 × 10 g/cm ; however, the results of case 1 suggest that a reasonably accurate answer can still be 1 0.1195% 0.00057% 20 2 4.2104% 0.0525% found. The NaN nodes are illustrated by dark lines in 3 296.6921% 24.0605% Fig. 13. Each layer of the polyhedron is represented by 1 0.0310% 0.00022% a third of the x- and y-axes, with the outermost layer 25 2 1.9908% 0.0133% shown on the top left. The concentration of dark lines 3 242.6755% 5.2795% in the bottom right implies that the inner layers are less 232 W. G. Ledbetter, R. Sood, J. Keane, et al. Fig. 11 1-Layer disturbed spheroid: effects of probe number on Adam optimization performance. No strong improvement is observed by adding more measurement data. Fig. 12 1-Layer disturbed spheroid: convergence behavior. The availability of more comprehensive measurements allows the algorithm to maintain the solution more steadily. discernible than the outer layer. to diverge faster. In Fig. 15, unlike previous cases, returns beyond those of the 15-probe simulations appear to be The results for this 3-layer case are consistent with the diminishing. behavior of the 1-layer case. Error reduction is consistent, if slow, and the algorithm produces an improvement of −5 3 5 Conclusions approximately 2.5 × 10 g/cm over 50,000 iterative cycles. It is also apparent that the number of probes is This paper proposes a new technique for small body not a major factor in the initial stages of recovery. The gravimetry and analysis. Based on recent developments asymptotic shape of Fig. 14 has become a common theme in sensor technology, a method capable of recovering for these initial plots. the internal density distribution of a polyhedral body The behavior close to the solution is also in good using only the position and gravity gradient data was correspondence with previous results. The addition of developed. A multi-layer gravity model was developed more probes enables the solution to be more effectively with reference to previous work, and an expression for maintained, whereas simulations with fewer probes tend the gravity gradient was derived as a linear function SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 233 Fig. 13 Case 3: correlations of 3-layer disturbed spheroid with 10 probes (dark blue = unobservable). The concentration of dark, unobservable regions in the bottom right coincides with the aforementioned results and Park et al.’s conclusion [6]. Fig. 14 3-Layer disturbed spheroid: effects of probe number on Adam optimization performance. Even for a more complex body, the number of probes has a minimal effect on the initial phases of optimization. of the density. This formulation was then used as a ability of a body and the empirical performance of the basis for recovering the gravity model from simulated Adam algorithm. Even with significant positional noise measurements. Both a GPU-based QR factorization applied to the dataset, Adam was able to consistently inversion algorithm and a stochastic method from the move from a decent initial guess to a more accurate solu- field of deep learning were employed to solve the matrix tion. In future, we intend investigating the disconnection problem. between the predicted performance and test results, and The most significant conclusion from the research is the aim to attempt to develop techniques and heuristics for apparent disagreement between the theoretical observ- more effective accuracy prediction. 234 W. G. Ledbetter, R. Sood, J. Keane, et al. Fig. 15 3-Layer disturbed spheroid: convergence behavior. The availability of more data points again helps to maintain the solution, but diminishing returns are observed beyond 15 probes. Acknowledgements [6] Park, R. S., Werner, R. A., Bhaskaran, S. Estimat- ing small-body gravity field from shape model and na- The researcher would like to thank the University of vigation data. Journal of Guidance, Control, and Dy- Alabama Graduate Council for financially supporting namics, 2010, 33(1): 212–221. this work. Additionally, Dr. Kieran A. Carroll, CTO [7] Carroll, K., Faber, D. Tidal acceleration gravity gra- of Gedex, is thanked for providing guidance at the 49th diometry for measuring asteroid gravity field from orbit. LPSC when this concept was still in development. A por- In: Proceedings of the 69th International Astronautical tion of this research was carried out at the Jet Propulsion Congress, 2018: IAC-18-A3.4B.3. Laboratory, California Institute of Technology, under a [8] Werner, R. A., Scheeres, D. J. Exterior gravitation of a polyhedron derived and compared with harmonic contract with the National Aeronautics and Space Ad- and mascon gravitation representations of asteroid 4769 ministration (No. 80NM0018D0004). Castalia. Celestial Mechanics and Dynamical Astronomy, 1996, 65(3): 313–344. References [9] Kaasalainen, M. Optimization methods for asteroid lightcurve inversion: I. shape determination. Icarus, [1] Liao-Mcpherson, D., Dunham, W. D., Kolmanovsky, 2001, 153(1): 24–36. I. Model predictive control strategies for constrained soft landing on an asteroid. In: Proceedings of the [10] Shepard, M. K., Richardson, J., Taylor, P. A., Rodriguez- AIAA/AAS Astrodynamics Specialist Conference, 2016: Ford, L. A., Conrad, A., de Pater, I., Adamkovics, M., 2016-5507. de Kleer, K., Males, J. R., Morzinski, K. M. et al. Radar observations and shape model of asteroid 16 Psyche. [2] Atchison, J., Mitch, R., Rivkin, A. Swarm flyby gravime- Icarus, 2017, 281: 388–403. try. Technical report. The Johns Hopkins University Applied Physics Laboratory, 2015. [11] Bercovici, B., McMahon, J. W. Improved shape determi- [3] Scheeres, D. J., Khushalani, B., Werner, R. A. Estimat- nation for autonomous state estimation. In: Proceedings of the 31st International Symposium on Space Technol- ing asteroid density distributions from shape and gravity ogy and Science, 2017. information. Planetary and Space Science, 2000, 48(10): 965–971. [12] Cocaud, C., Kubota, T. Autonomous navigation near [4] Takahashi, Y., Scheeres, D. J. Morphology driven density asteroid based on visual SLAM. In: Proceedings of the 23rd International Symposium on Space Flight Dynam- distribution estimation for small bodies. Icarus, 2014, 233: 179–193. ics, 2012. [5] Ditmar, P., Klees, R., Kostenko, F. Fast and accu- [13] Bercovici, B., McMahon, J. W. Point-cloud processing rate computation of spherical harmonic coefficients from using modified Rodrigues parameters for relative na- satellite gravity gradiometry data. Journal of Geodesy, vigation. Journal of Guidance, Control, and Dynamics, 2003, 76(11–12): 690–705. 2017, 40(12): 3167–3179. SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets 235 [14] Takahashi, Y., Scheeres, D. J. Small body surface gra- ratories. E-mail: wgledbetter@crimson.ua.edu. vity fields via spherical harmonic expansions. Celestial Mechanics and Dynamical Astronomy, 2014, 119(2): Rohan Sood is an assistant professor of 169–206. aerospace engineering at the University [15] Kingma, D. P., Ba, J. Adam: A method for stochastic of Alabama. He leads the Astrodynamics optimization. arXiv preprint, 2014, arXiv:1412.6980. and Space Research Laboratory where his [16] Reddi, S. J., Kale, S., Kumar, S. On the convergence of team investigates cost-effective, innova- Adam and beyond. In: Proceedings of the International tive, resilient trajectories and spacecraft Conference on Learning Representations, 2018. behavior. Some of the research topics are [17] Dozat, T. Incorporating Nesterov momentum into Adam. asteroid exploration, solar sailing, multi- body dynamics, spacecraft attitude control, and applications [18] Yeomans, D. K. Radio science results during the NEAR- of artificial intelligence in astrodynamics. Dr. Sood received shoemaker spacecraft rendezvous with Eros. Science, his bachelor of science degree in aerospace engineering from 2000, 289(5487): 2085–2088. the University at Buffalo. He received his Ph.D. degree in [19] Cazenave, A. Geoid, Topography and Distribution of aeronautics and astronautics with concentration in astrody- Landforms. Washington, D. C.: American Geophysical namics and space applications from Purdue University under Union, 1995: 32–39. the mentorship of Dr. Kathleen C. Howell and Dr. Henry [20] Lerch, F. J., Nerem, R. S., Putney, B. H., Felsentreger, (Jay) Melosh. Dr. Sood was a member of the GRAIL Science T. L., Sanchez, B. V., Klosko, S. M., Patel, G. B., Team. He has since collaborated with NASA MSFC, ARC, Williamson, R. G., Chinn, D. S., Chan, J. C. et al. and JPL on upcoming missions and to explore innovative Geopotential models of the Earth from satellite tracking, mission concepts. Dr. Sood is a pianist, an artist (oil on altimeter and surface gravity observations: GEM-T3 canvas, charcoal), and enjoys participating in competitive and GEM-T3S. NASA Technical Memorandum 104555. sports. E-mail: rsood@eng.ua.edu. NASA, Goddard Space Flight Center Greenbelt, MD, 1992. James Tuttle Keane hails from Cedar Rapids, Iowa. He received his bachelor William Gordon Ledbetter is in his degree in astronomy and geology from fourth year of Ph.D. study after receiving the University of Maryland, College his bachelor of science degree in aerospace Park, and his Ph.D. degree in planetary engineering from The University of Al- science from the University of Arizona, abama, summa cum laude. He has been under the mentorship of Dr. Isamu Mat- interested in space for as long as he can suyama. After a postdoctoral position remember, and learned how to code by at the California Institute of Technology, he started at JPL reading the TI-83 Plus manual. As the in 2020. Dr. Keane is a planetary scientist, studying the capstone of his senior year he participated in design projects interactions between orbital dynamics, rotational dynamics, performing mission architecture analysis for proposed missions and geologic processes on rocky and icy worlds across the to Mars. William is currently exploring autonomous mission solar system. He uses a combination of theoretical methods, design for swarm exploration of asteroids, including landing coupled with the analysis of spacecraft-derived datasets to and rendezvous. He is also investigating embedded systems for investigate the dynamics, structure, origin, and evolution of rapid, in-situ multi-agent optimal control. William received solar system bodies. He has extensive experience with NASA the 2017–2018 Graduate Council Fellowship and was named missions, including GRAIL, New Horizons, and the proposed the 2019 Aerospace Engineering graduate student of the year Io Volcano Observer. In addition to science, Dr. Keane at UA by the Engineering Council of Birmingham. William is an accomplished science illustrator and communicator, was selected as the 2019–2020 and 2020–2021 Alabama Space specializing in pen and pencil illustrations of planetary Grant Consortium Fellow. Beginning in the summer of 2020, science concepts, new results, and exciting missions. E-mail: William has an ongoing internship at Sandia National Labo- james.t.keane@jpl.nasa.gov. 236 W. G. Ledbetter, R. Sood, J. Keane, et al. Jeffrey Stuart received his B.S., M.S., mons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in and Ph.D. degrees in aeronautics and any medium or format, as long as you give appropriate credit astronautics from Purdue University in 2008, 2011, and 2014, respectively. He is to the original author(s) and the source, provide a link to currently technical staff in Mission Design the Creative Commons licence, and indicate if changes were made. & Navigation at the Jet Propulsion Labo- The images or other third party material in this article are ratory, California Institute of Technology. He is the MDNav Lead for the SunRISE included in the article’s Creative Commons licence, unless formation flying mission and works on a variety of projects indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and from early formulation to flight operations. His interests in- your intended use is not permitted by statutory regulation or clude automated trajectory design, advanced navigation tech- niques, formation flying, combinatorial optimization, and vi- exceeds the permitted use, you will need to obtain permission sual analytic methods. E-mail: Jeffrey.r.stuart@jpl.nasa.gov. directly from the copyright holder. To view a copy of this licence, visit http://creativecoorg/ licenses/by/4.0/. Open Access This article is licensed under a Creative Com-
Astrodynamics – Springer Journals
Published: Sep 1, 2021
Keywords: gravimetry; asteroid exploration; astrodynamics; stochastic optimization
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.