Modular pipeline for small bodies gravity field modeling: an efficient representation of variable density spherical harmonics coefficients
Antonio Rizza Carmine Buonagura Paolo Panicucci Francesco Topputo


75
Milan, Italy
14-18 October
2024
A3




*Antonio RizzaPhD Candidate,  Department of Aerospace Science and Technology (DAER), Politencico di Milano, antonio.rizza@polimi.it
Carmine BuonaguraPhD Candidate, Department of Aerospace Science and Technology (DAER), Politencico di Milano, carmine.buonagura@polimi.it
Paolo PanicucciAssistant Professor, Department of Aerospace Science and Technology (DAER), Politencico di Milano,paolo.panicucci@polimi.it
Francesco TopputoFull Professor, Department of Aerospace Science and Technology (DAER), Politencico di Milano, francesco.topputo@polimi.it



Proximity operations to small bodies, such as asteroids and comets, demand high levels of autonomy to achieve cost-effective, safe, and reliable Guidance, Navigation and Control (GNC) solutions. Enabling autonomous GNC capabilities in the vicinity of these targets is thus vital for future space applications. However, the highly non-linear and uncertain environment characterizing their vicinity poses unique challenges that need to be assessed to grant robustness against unknown shapes and gravity fields. In this paper, a pipeline designed to generate variable density gravity field models is proposed, allowing the generation of a coherent set of scenarios that can be used for design, validation, and testing of GNC algorithms. The proposed approach consists in processing a polyhedral shape model of the body with a given density distribution to compute the coefficients of the spherical harmonics expansion associated with the gravity field. To validate the approach, several comparison are conducted against analytical solutions, literature results, and higher fidelity models, across a diverse set of targets with varying morphological and physical properties. Simulation results demonstrate the effectiveness of the methodology, showing good performances in terms of modeling accuracy and computational efficiency. This research presents a faster and more robust framework for generating environmental models to be used in simulation and hardware-in-the-loop testing of onboard GNC algorithms.

Modular pipeline for small bodies gravity field modeling: an efficient representation of variable density spherical harmonics coefficients
    Ahmad T. Sheikh, Ali Shoker, Suhaib A. Fahmy and Paulo Esteves-Verissimo
CEMSE Division, King Abdullah University of Science and Technology (KAUST)

Thuwal 23955-6900, Kingdom of Saudi Arabia

{ahmad.sheikh, ali.shoker, suhaib.fahmy, paulo.verissimo{@kaust.edu.sa}}
    September 9, 2024
=============================================================================================================================================================================================================================================================================






§ INTRODUCTION

In the last twenty years, the space sector has experienced an unprecedented growth, marked by significant advancements and achievements both concerning Earth's orbiting and deep-space missions. The recent growing interest in small solar system bodies such as asteroids and comets for scientific inspection, exploitation of resources, and planetary defense reasons is pushing the development of innovative engineering solutions to better investigate these celestial bodies. Ground-based observations allow preliminary characterizations of small bodies in terms of bulk properties, such as mass and shape, orbit and rotational state. The major limitation of this methodology is the signal to noise ratio <cit.> which is acceptable only when the target is relatively close to the Earth  <cit.>. A drastic improvement in the body characterization can be obtained with in-situ observations with the use of specialized and instrumented probes. Several missions successfully performed proximity operations to these bodies such as the Near Earth Asteroid Rendezvous (NEAR) Shoemaker <cit.>, Dawn <cit.>, the Origins, Spectral Interpretation, Resource Identification, Security, Regolith Explorer (OSIRIS-REx) <cit.>, Hayabusa <cit.>, , , and the Double Asteroid Redirection Test (DART)<cit.>. Nowadays, space exploration is witnessing a transition towards the use of CubeSats, miniaturized platforms standardized in size and form factor, for the systematic exploration of the Solar System <cit.>. The idea is to use these platforms performing riskier tasks and operating in a cooperative multi-agent framework while cooping with limited resources <cit.>. 


The dynamics characterizing small bodies proximity operations is highly chaotic with solar radiation pressure and non-spherical gravity effects compromising the very existence of stable closed orbits. An accurate modeling of this environment becomes crucial for a safe trajectory design and to guarantee the satisfaction of mission objectives. Moreover, on-board operations requires fast computation to coop with limited available resources. A good compromise between accuracy and efficiency in modeling the gravity field, is typically found in an high-order expansion of the gravitational potential in the form of spherical harmonics. This model is typically valid only outside the Brillouin sphere of the target <cit.>, because of the convergence properties of the Legendre polynomials, used as functional base for the expansion <cit.>. An alternative formulation resolving the convergence issue inside the Brillouin sphere is the polyhedral gravity model proposed in Werner and Scheeres <cit.>. This model computes the gravity field from a constant density polyhedron in an analytical form. Given that the polyhedron is usually computed from images, errors are present in the polyhedral model. The sensitivity of the gravity field to perturbation in the polyhedral  shape are investigated by Bercovici et al. <cit.>, to assess the coupling between shape and gravity errors for Werner and Scheeres's <cit.> model. One of the main drawback of the polyhedral gravity model is that it is associated with a constant density distribution inside the body which sometimes results in contrast with gravity measurements <cit.>. Even a small density variation may induce large trajectory deviation with respect to a nominal path, leading the spacecraft to miss its target goals and potentially enter on impact or escape trajectories. A quantitative example of this deviation is discussed later in this work using the variable density gravity model generated with the proposed approach.   


Motivated by the need of comparing spherical harmonics coefficients from orbit determination and the ones gathered from constant-density shapes, Werner <cit.> shows how the spherical-harmonics coefficients can be analytically retrieved from a general shape polyhedral with uniform density model of the asteroid combining recursive formula with trinomial algebra.  The sensitivity of the spherical harmonics coefficients to shape variations is investigated by Panicucci et al. <cit.> to understand whether the shape can be a major source of error when comparing orbit determination coefficients to shape-deduced coefficients (e.g., see <cit.>)  To improve the fidelity of the forward gravity modeling of non-constant density polyhedron, Chen et al. <cit.> extends Werner <cit.>  to variable density under the assumption of trinomial density distribution. However, this hypothesis force the density field to be continuous inside the body forbidding the existence of density jumps. This may not be the case when the asteroid is formed by two different parent bodies or when denser nuclei and internal geological structures are considered.


Another widely-used gravity field is mascon model, introduced to model gravity anomalies on the Moon <cit.>. Mascon models <cit.> are typically used to model non-uniform mass distribution. An interesting comparison among different gravity models, spherical harmonics, mascon and polyhedral, is provided by Werner and Scheeres <cit.> for asteroid (4769) Castalia. To the authors knowledge there is currently no formulation of the spherical harmonics expansion under arbitrary variable density distribution and shape of the asteroid. This paper proposes a pipeline for retrieving such gravity model by performing analytical integration over a radially discretized polyhedron. This approach falls thus in between the methodologies proposed by Werner and Cheng <cit.> and the classical Mascon approach. 


The paper is structured as follows: Section <ref> describes the general formulation of the gravity field by means of a spherical-harmonics expansion, Section <ref> presents the methodology to compute the coefficients of the expansion with non-uniform density, Section <ref> applies the methodology to benchmark cases to test its performances and, finally, Section <ref> summarizes the major findings. 



§ GRAVITY FIELD MODELING

The acceleration due to the central gravity field, expressed in the asteroid fixed frame ℬ, can be computed as the gradient of the gravitational potential U(x⃗)


    a⃗_ℬ = ∇ U


Being due to a conservative field, the potential has to satisfy the Laplace equation outside the body, i.e.


    ∇^2 U (x⃗) = 0


A known solution to the Laplace equation, valid outside the Brouillon  of the asteroid, is expressed in spherical coordinates through an infinite expansion of the potential U(r,λ,ϕ) into a series of spherical harmonics projected onto a function space spanned by the associated Legendre . The coordinates r, λ and ϕ are respectively the radial distance from the asteroid, the longitude and the latitude of the spacecraft. To avoid numerical errors and improve accuracy for higher order terms this expansion is typically used in the normalized form


    U  =μr∑_n=0^∞∑_m=0^n(R_0r)^n nmu(nmcos(m λ)+.
     +.nmsin(m λ) )


where u = sin(ϕ). The function nmu refers to the normalized associated Legendre polynomials, while, the scaling factor R_0 is a reference distance used to compute the normalized coefficients nm, and nm. The polynomials nmu are expressed as a function of the Legendre polynomials P_n as <cit.> 


    nmu = N_nm1-u^2^m2d^m P_ndu^m

with 

    N_nm = √((n-m) !  (2 n +1) (2-δ(m) )( n+m) !)


and δ(·) indicating the Dirac delta operator. 


    δ(m)  = {[         1     m = 0;         0 otherwise ].




To speed up the computation of the fully-normalized Legendre polynomials in Equation <ref> the recursive formula introduced by Rapp <cit.> are used. The sequence is anchored on the first three elements as


    00u = 1        
    
       10u = √(3) u       
    
       22u = √(3)√(1-u^2)


and then the terms are recursively computed as in Equation <ref>.


    For  m  =    n,   n > 2,
    
    nnu =    √(2n+12n)√(1-u^2)n-1n-1u
    
    for  m =    n-1,
    
    nn-1u =   √(2n+3) u n-1n-1u
    
    for  m <    n-1,
    
    nmu =   Γ_n,mn-1mu +
       - Γ_n,mΓ_n-1,mn-2mu



with


    Γ_n,m = √((2n+1) (2n-1)(n-m) (n+m))


Expressing the spacecraft position in Cartesian coordinates, i.e. r⃗_ℬ = x i⃗+y j⃗+z k⃗, the gradient of the pontetial can be expressed as <cit.>


    a⃗_ℬ = ∇ U =[∂ U∂r⃗]^T = 
    
    
            =    ∂ U∂ r[∂ r∂r⃗]^T +  ∂ U∂λ[∂λ∂r⃗]^T +  ∂ U∂ϕ[∂ϕ∂r⃗]^T = 
    
    
       =    [ (1r∂ U∂ r - zr^2 η∂ U ∂ϕ) x - (1η^2∂ U∂λ) y] i⃗ + 
    
     
    		+   [ (1r∂ U∂ r - zr^2 η∂ U ∂ϕ) y + (1η^2∂ U∂λ) x] j⃗ + 
    
    
    		+   [ 1r∂ U∂ r z + ηr^2∂ U∂ϕ] k⃗


with η = √(x^2+y^2). The partial derivatives ∂ U∂ r and ∂ U∂λ are given by



    ∂ U∂ r = -μr^2∑_n=0^∞∑_m=0^n(R_0r)^n (n+1 ) nmu
    (nmcos(m λ) + nmsin(m λ) )




    ∂ U∂λ = μr∑_n=0^∞∑_m=0^n(R_0r)^n (n+1 ) nmu m 
    (-nmsin(m λ)+nmcos(m λ) )


The derivative ∂ U∂ϕ is more challenging since it involves the derivatives of the associated Legendre polynomial


    ∂ U∂ϕ = μr∑_n=0^∞∑_m=0^n(R_0r)^n ∂nmsin(ϕ)∂ϕ
    (nmcos(m λ)+ nmsin(m λ) )


There are several recursive formulas to compute the derivative of the fully normalized Legendre . For this implementation this is achieved by combining the definition of the normalized Legendre polynomials, in Equation <ref>, with the set of recursive formula in Equation <ref> leading to


    ∂nmsin(ϕ)∂ϕ = [∂nmu∂ u]  cosϕ = 
       
        =  -m tanϕnmu + Kn,mnm+1u


with 


    K_n,m  = {[ √(n-mn+m+1)       m > 0;                        ;   √(n n+12)       m = 0 ].









§ NORMALIZED COEFFICIENTS COMPUTATION

The coefficients nm and nm are a function of the asteroid shape and density distribution. In particular they are obtained integrating a set of trinomial shape functions c_nm and s_nm over the body <cit.>


    [ nm; nm ] =    ∭_B [ c_nm(x,y,z); s_sm(x,y,z) ] dm = 
     =    ∭_B ρ[ c_nm(x,y,z); s_sm(x,y,z) ] dx dy dz


 If the integration is performed over a polyhedron, split in a collection of tethraedra, the full integral is equivalent to the summation of the integrals over each 


    [ nm; nm ] =   ∑_s = 1^n_s( ∭_s ρ[ c_nm(x,y,z); s_sm(x,y,z) ] dx dy dz)


Lien and Kajiya <cit.> shows that, under the hypothesis of constant density, each of these integrals, can be solved analytically if a proper change of coordinates is performed.


A similar approach is followed here to derive the formulation with variable density distribution. Consider the single polyhedron shown in Figure <ref>a


The general position inside the tetrahedron can be identified as a linear combination of the vertex coordinates C⃗_1, C⃗_2 and C⃗_3 


    [ x; y; z ] =  [ C⃗_1 C⃗_2 C⃗_3 ][ X; Y; Z ]


or better expressed in compact form as


    x⃗ =J⃗_s X⃗

where the matrix J⃗_s is the Jacobian of the transformation. This change of coordinates allows to perform the integration in Equation <ref> over a standard simplex, as shown in Figure <ref>b. From Equation <ref> thus follows


    [ nm; nm ] =    ∑_s = 1^n_s∭_ss( ρ[ c_nm(X⃗); s_sm(X⃗) ] det(J⃗_s))  dX⃗


Since c̅_nm(x⃗) and s̅_nm(x⃗) are trinomials of order n, and the transformation in Equation <ref> is linear, c̅_nm(X⃗) and s̅_nm(X⃗) will also be trinomials of order n. The shape functions can then be expressed as 


    [ c_nm(X⃗); s_sm(X⃗) ] =    ∑_i+j+k = n([ α_ijk; β_ijk ] X^i Y^j Z^k) = 
    
      =    ∑_i+j+k = n(p⃗_ijk X^i Y^j Z^k)


The components α_ijk and β_ijk of vector p⃗_ijk can be retrieved through a series of recursive relations anchored on certain initial conditions. In this paper the details of this procedure are omitted and the coefficients can be considered to be known. The interested reader can find a comprehensive discussion of this process in <cit.>. 


Combining Equations <ref> and <ref> the steps in  hold. Note that this does not contain any assumption on the internal density distribution of the body.  



Consider now the uniform discretization of the polyhedron, shown in Figure <ref>, in n_q segments each of them at a constant density ρ_q = ρ(x⃗_q) with x⃗_q being the geometric center of segment q. 



The integration over the standard simplex can then be split further more in a summation of integrals over each segment, leading to


    [ nm; nm ] =   ∑_s = 1^n_s( det(J⃗_s) ∑_i+j+k = n([ α_ijk; β_ijk ]∑_q=1^n_qρ_q . .
       . . (∫_0^q^+   g_ijk(Z)  dZ-∫_0^q^- g_ijk(Z) dZ)))


with 


    g_ijk(Z) = ∫_0^1-Z∫_0^1-Z-Y  X^i Y^j Z^k dX dY dZ


Lien and Kajiya <cit.> shows that this integral can be computed analytically as 


    g_ijk(Z) = β(j+1,i+2)i+1(1-Z)^(i+j+2)Z^k


where β is the Beta function, also known as Euler's integral, with tabulated values. By substituting the definition of g_ijk(Z) in Equation <ref> follows that 


    [ nm; nm ] =    ∑_s = 1^n_s( det(J_s) ∑_i+j+k = n([ α_ijk; β_ijk ]. . 
       . . ∑_q=1^n_q(ρ_q (h_q^+,i,j,k. .
            - h_q^-,i,j,k) ) ) )


with 


    h_q^+,i,j,k = β(j+1,i+2)i+1β(q,k+1,i+j+3)


and β̃ is the incomplete beta function, also known and tabulated.



§ VALIDATION

 The expression provided in Equation <ref> provides an exact way of computing the normalized coefficients nm and nm starting from a generic shape model of the asteroid and a generic density distribution ρ(x⃗). Apart from defining the shape of the gravity field, the normalized coefficients provide also important information on the inertia properties of the body. For example, the Center of Mass (CoM) of the body can be derived as: 


    r⃗_CoM,ℬ = [ 11; 11; 10 ]√(3) R_0


This can be used to validate the proposed methodology. Five test cases are investigated: the first three are toy problems designed to validate the approach comparing the CoM position computed from the coefficients with known analytical solutions while the other two presents a more realistic application case to asteroid (433) Eros in which the acceleration computed with spherical harmonics model is compared against the one obtained with a Mascon model.


The first case is the one of a uniform density sphere with density  and radius R_0 = 500 m. Being a sphere with constant density the CoM is clearly in its center, as illustrated in Figure <ref>(a). The second case considers the same body but introduces a density gradient between two hemispheres, see Figure <ref>(b). In particular the density is assumed to be


    ρ( x⃗) = ρ_1 = 3204  kg/m^3      x ≥ 0 
    ρ_2 = 1335  kg/m^3      x < 0

In this case the CoM will only have an x component that can be analytically computed as 


    x_CoM = 38 R_0 ( ρ_1 - ρ_2ρ_1 + ρ_2)


The third case considers instead the asteroid (486958) Arrokoth, selected for its peculiar shape, see . A uniform density distribution is assumed here, with a density of 235  kg/m^3 and a center of mass position , both estimated by Keane et al.<cit.>. 




Normalized spherical harmonics coefficients are computed for these three cases using n_q = 10 in the methodology presented in this section, and the CoM is computed as in Equation <ref>. Table <ref> shows the accuracy achieved in the three cases. The table shows that an accuracy always below 1 meter is achieved even when the CoM location is not analytically computed but retrieved from literature as in the Arrokoth case. 



The last two examples refers to asteroid Eros. Thanks to extensive mapping performed by NEAR <cit.>, Eros is the small body for which higher fidelity models are available. A good description of its gravity field is provided by Garmier et al. <cit.>. While the study reconstruct the asteroid density to be almost uniform inside the body, gravitational anomalies are detected that could be due to the existence of denser areas in the central .


To test the accuracy of the gravity model the acceleration obtained from Equation <ref> is compared with the one computed with a variable density Mascon model <cit.> with the same density distribution assumption. To make the computation affordable, the model is derived starting from a 50k polyhedron model[<https://3d-asteroids.space/asteroids/433-Eros>] and performing a mesh simplification to obtain a polyhedron with 12290 vertices and 24576 faces used to compute the gravitational coefficients. 


First, a uniform density model is considered with a density 2670 kg/m^3. Figure <ref> shows the error in the acceleration among the two models over an ellipsoid with semi-major axis 30 km× 20 km× 20 km seen from the asteroid north pole. The figure shows an error that is always below 1  mgal = 1e-5 m/s^2.  This results to be less than one order of magnitude lower than the perturbation due non spherical gravity at that distance from Eros as shown in Rizza et al. <cit.>. Also note that the picks in the error are clearly linked with the closest regions to the Brillouin sphere.





Then a variable density distribution with an inner core with a density 10% larger than the nominal values is considered. In particular


    ρ( r⃗) = ρ_1 = 2937  kg/m^3      r ≤ 5  km
    ρ_2 = 2670  kg/m^3      r > 5  km


This density profile, shown in Figure <ref>, leads to an acceleration error that is slightly larger than the one of the previous case, see Figure <ref> but still always below 1  mgal.





To show the relevance of density variation on the spacecraft trajectory a comparison is shown in the following. Considering the initial condition shown on , the trajectory is propagated for 24 h using the gravity field generated from the two previously described density models. 






Inertial and asteroid fixed trajectories are shown in Figure <ref> while the error on position and velocity is reported in Figure <ref>.





The observed effect of the variable density gravity field is to induce a cumulative error of more than 10 km and around 0.7 m/s at the end of the considered time horizon. 




§ CONCLUSIONS

In this paper a pipeline for semi-analytical computation of normalized spherical-harmonics coefficients is proposed. The algorithm takes as input only a shape model of the asteroid and a generic density distribution function ρ(x⃗). Different levels of accuracy can be achieved simply by increasing the fidelity of the shape model or the number of radial discretization intervals n_q. The methodology has been tested against analytical results, literature findings and numerical simulations, showing good accuracy. The set of generated coefficients could be used on-board to efficiently compute the effects of non-uniform density distribution without relying on computationally heavier gravity models which are difficult to run on-board.  


§ REFERENCES

[heading=none]

