Constraining the Thermal Properties of Planetary Surfaces using Machine Learning: Application to Airless Bodies

TL;DR

This paper introduces a machine learning-based method to determine surface thermal properties of airless planetary bodies from infrared flux measurements, validated on asteroid data and applicable to upcoming space missions.

Contribution

The study develops a neural network surrogate model combined with Bayesian inversion to accurately retrieve surface properties from infrared data, advancing remote sensing techniques for planetary surfaces.

Findings

Successfully retrieved surface properties of asteroids Bennu and Itokawa.

Separated contributions of regolith and rock components in thermal measurements.

Estimated rock abundance and thermal inertia consistent with previous in situ data.

Abstract

We present a new method for the determination of the surface properties of airless bodies from measurements of the emitted infrared flux. Our approach uses machine learning techniques to train, validate, and test a neural network representation of the thermophysical behavior of the atmosphereless body given shape model, illumination and observational geometry of the remote sensors. The networks are trained on a dataset of thermal simulations of the emitted infrared flux for different values of surface rock abundance, roughness, and values of the thermal inertia of the regolith and of the rock components. These surrogate models are then employed to retrieve the surface thermal properties by Markov Chain Monte Carlo Bayesian inversion of observed infrared fluxes. We apply the method to the inversion of simulated infrared fluxes of asteroid (101195) Bennu -- according to a geometry of observations similar to those planned for NASA's OSIRIS-REx mission -- and infrared observations of asteroid (25143) Itokawa. In both cases, the surface properties of the asteroid -- such as surface roughness, thermal inertia of the regolith and rock component, and relative rock abundance -- are retrieved; the contribution from the regolith and rock components are well separated. For the case of Itokawa, we retrieve a rock abundance of about 85% for pebbles larger than the diurnal skin depth, which is about 2 cm. The thermal inertia of the rock is found to be lower than the expected value for LL chondrites, indicating that the rocks on Itokawa could be fractured. The average thermal inertia of the surface is around 750 $J s^{-1/2} K^{-1} m^{-2}$ and the measurement of thermal inertia of the regolith corresponds to an average regolith particle diameter of about 10 mm, consistently with in situ measurements as well as results from previous studies.

Tables

Table 1. Table 1 : Description of the dataset of forward thermal simulations of asteroid (101955) Bennu. The values of thermal inertia and rock abundance are equally spaced in a range; as an example, the diciture “10 : 10 : 500 J ​ s − 1 / 2 ​ K − 1 ​ m − 2 𝐽 superscript 𝑠 1 2 superscript 𝐾 1 superscript 𝑚 2 Js^{-1/2}K^{-1}m^{-2}  ” means “values from 10 to 500 J ​ s − 1 / 2 ​ K − 1 ​ m − 2 𝐽 superscript 𝑠 1 2 superscript 𝐾 1 superscript 𝑚 2 Js^{-1/2}K^{-1}m^{-2}  with a step of 10 J ​ s − 1 / 2 ​ K − 1 ​ m − 2 𝐽 superscript 𝑠 1 2 superscript 𝐾 1 superscript 𝑚 2 Js^{-1/2}K^{-1}m^{-2} ”. Each thermal simulation uses a spin state: λ 𝜆 \lambda = 0 ∘ superscript 0 0^{\circ} , β 𝛽 \beta = 90 ∘ , P 𝑃 P = 4.288 hours. The asteroid is located at 1 au. The asteroid is thermalized, i.e., it attains thermal equilibrium with the environment, in over several tens of rotational periods before recording the infrared flux. Namely, we begin each simulation where each facet has a surface temperature equal to its mean temperature calculated over the full diurnal curve and assuming zero thermal inertia. The temperature in the sub-surface of the facet is initially the same of the surface temperature. Next the appropriate value of the thermal inertia for the simulation is set and the heat diffusion into depth is calculated using Equation 12 of Spencer et al. ( 1989 ) with the boundary conditions given as in Equation 13 of Spencer et al. ( 1989 ) . The simulations are started at an epoch before the one when infrared fluxes are needed to be modeled and the temperatures of the facet (or some facets) are monitored such that the diurnal temperature curve at the asteroid rotation N 𝑁 N is within a fraction (we use typically 0.25 K) of a K from the diurnal temperature at the rotation N + 1 𝑁 1 N+1 . The number of ”thermalization” rotations is determined by trial and error.

Dataset of forward simulations (OTES-like acquisitions)
Predictors	Values
Hapke angle	4, 12, 14, 21, 27, 30, 34, 40, 46, 50, 56, 62 degrees
Thermal inertia, reg.	10 : 10 : 500 $J{s}^{-1/2}{K}^{-1}{m}^{-2}$
Thermal inertia, rock	500 : 10 : 1000 $J{s}^{-1/2}{K}^{-1}{m}^{-2}$
Rock abundance	0 : 5 : 100 %
Responses	Values
7 infrared fluxes	[$Jy$], 5 to 50 $\mu m$ with a step of 1 $\mu m$

Table 2. Table 2 : Architecture and performance of the surrogate models (neural networks) for OTES observations. The performance of the neural networks are assessed in terms of Mean Square Error (MSE, the average squared difference between outputs and targets; lower values are better) and Regression R value (the correlation between outputs and targets; an R value of 1 means a close relationship). The reported values are relative to the testing phase. We used the Matlab Levenberg-Marquardt algorithm for training, with 10-fold crossvalidation. The dataset is split as: 70% training, 15% validation and 15% testing.

Surrogate Models: Triangular Equatorial Facet
Station n.	Hidden layers	MSE	R
3:20 AM	1 layer, 20 neurons	0.53	0.999
6:00 AM	1 layer, 20 neurons	0.42	0.999
10:00 AM	1 layer, 10 neurons	0.90	0.999
12:30 PM	1 layer, 10 neurons	2.00	0.999
3:00 PM	1 layer, 10 neurons	0.90	0.999
6:00 PM	1 layer, 10 neurons	0.30	0.999
8:40 PM	1 layer, 10 neurons	0.90	0.999
Surrogate Models: Disk-integrated
Station n.	Hidden layers	MSE	R
03:20 AM	1 layer, 10 neurons	0.085	0.999
06:00 AM	1 layer, 10 neurons	0.037	0.999
10:00 AM	1 layer, 10 neurons	0.053	0.999
12:30 PM	1 layer, 10 neurons	0.124	0.999
03:00 PM	1 layer, 10 neurons	0.182	0.999
06:00 PM	1 layer, 10 neurons	0.066	0.999
8:40 PM	1 layer, 10 neurons	0.082	0.999

Table 3. Table 3 : Mid-infrared observations of (25143) Itokawa from Müller et al. ( 2005 , 2014 ) . The epoch of observations are relative to the emission from the surface of the asteroid (i.e., light-time corrected).

Thermal observations of (25143) Itokawa
N.	Epoch [JD]	$\lambda$ [$\mu m$]	Flux [$Jy$]	$\sigma$ [$Jy$]
1	2451982.743	11.66	0.264	0.044
2	2452007.894	11.66	0.164	0.021
3	2452007.904	10.38	0.144	0.018
4	2452007.917	12.35	0.17	0.022
5	2452007.929	8.73	0.086	0.022
6	2452007.94	11.66	0.149	0.019
7	2452008.894	12.35	0.258	0.032
8	2452008.906	9.68	0.108	0.016
9	2452008.919	10.38	0.169	0.027
10	2452008.929	11.66	0.242	0.03
11	2452008.939	11.66	0.193	0.028
12	2453187.752	8.73	1.92	0.15
13	2453187.763	8.73	1.97	0.16
14	2453187.775	8.73	1.75	0.14
15	2453187.788	8.73	1.67	0.13
16	2453187.803	10.68	1.94	0.14
17	2453187.817	10.68	1.89	0.13
18	2453187.828	12.35	2.17	0.13
19	2453187.84	12.35	1.8	0.11
20	2453187.859	17.72	2.49	0.5
21	2453196.989	11.7	0.762	0.100
22	2453196.991	11.7	0.721	0.091
23	2453197.064	11.7	0.913	0.114
24	2453197.07	9.8	0.791	0.125
25	2453197.077	9.8	0.570	0.122
26	2454307.977	4.1	0.00032	0.00025
27	2454307.977	7.0	0.00469	0.00028
28	2454307.976	11.0	0.01422	0.00053
29	2454308.046	15.0	0.02137	0.00079
30	2454308.048	24.0	0.01947	0.00120

Table 4. Table 4 : Description of the datasets of forward simulations for the thermal model of asteroid (25143) Itokawa. The prime meridian is defined according to the IAU specifications (Gaskell et al. 2008 ) . Each thermal simulation uses a spin state: λ 𝜆 \lambda = 268.77 ∘ , β 𝛽 \beta = -89.53 ∘ , P 𝑃 P = 12.1324 hours. The asteroid is thermalized in over thirty rotational periods before recording the infrared flux.

2-D model (Itokawa)
Predictors	Values
Hapke angle	4, 12, 14, 21, 27, 30, 34, 40, 46, 50, 56, 62 degrees
Thermal inertia, av.	25: 25 : 1500 $J{s}^{-1/2}{K}^{-1}{m}^{-2}$
Responses	Values
Infrared flux	30x1 fluxes $\left[Jy\right]$

Table 5. Table 5 : Metrics of performance of the surrogate models (neural networks) for the thermal observations of (25143) Itokawa, deployed on the 2-D dataset and on the 4-D dataset. For the 2-D case, all the networks have a single hidden layer with 30 neurons. For the 4-D case, all the networks have a single hidden layer with 10 neurons, except networks n. 17 and 20, that have 20 neurons. The reported values are relative to the testing phase. The average Mean Squared Error of the surrogate models is within the noise level of the observations. We used the Matlab Levenberg-Marquardt algorithm for training, with 10-fold crossvalidation. The dataset is used as: 70% training, 15% validation and 15% testing.

Surrogate models for (25143) Itokawa
N.	MSE (2D)	R (2D)	MSE (4D)	R (4D)
1	1.0E-07	0.999	1.5E-05	0.998
2	3.0E-07	0.999	3.0E-05	0.998
3	6.0E-08	0.999	1.5E-05	0.998
4	5.0E-08	0.999	1.5E-05	0.998
5	3.8E-08	0.999	5.4E-06	0.991
6	5.0E-08	0.999	9.0E-06	0.991
7	1.4E-07	0.999	3.9E-05	0.991
8	5.0E-07	0.999	1.2E-05	0.991
9	1.3E-07	0.999	1.5E-05	0.991
10	1.3E-07	0.999	1.5E-05	0.991
11	2.6E-07	0.999	1.0E-05	0.999
12	7.8E-08	0.999	4.0E-04	0.999
13	4.0E-07	0.999	6.0E-04	0.999
14	6.0E-07	0.999	3.0E-04	0.999
15	3.0E-07	0.999	1.0E-04	0.999
16	2.9E-07	0.999	1.0E-04	0.999
17	3.0E-07	0.999	4.0E-04	0.999
18	2.0E-07	0.999	5.0E-04	0.999
19	9.7E-08	0.999	1.0E-04	0.999
20	2.0E-07	0.999	5.0E-05	0.999
21	4.0E-07	0.999	9.8E-06	0.999
22	8.0E-08	0.999	8.5E-06	0.999
23	3.8E-08	0.999	5.0E-05	0.999
24	1.2E-07	0.999	4.0E-05	0.999
25	4.0E-08	0.999	1.0E-05	0.999
26	3.9E-09	0.999	7.0E-08	0.999
27	8.2E-09	0.999	3.0E-08	0.999
28	3.0E-09	0.999	1.0E-08	0.991
29	1.4E-08	0.999	3.0E-07	0.910
30	1.3E-08	0.999	2.0E-07	0.900

Table 6. Table 6 : A measurement of the thermal inertia of the regolith can be converted into mean radius of the grain. Such a conversion requires the knowledge of the listed parameters. Some of these parameters – noticeably the specific heat, c p subscript 𝑐 𝑝 c_{p} , and the particle thermal conductivity, k 𝑘 k – are temperature dependent and have been therefore computed at a mean surface temperature of 300 K, relative to the portion of the asteroid in the field of view (the heliocentric distance is about 1 au for each observation).

Property	Value	Source
Particle Density	3220 $kg{m}^{-3}$	Consolmagno et al. (2008), average value for LL chondrites
Specific Heat	682 $Jk{g}^{-1}{K}^{-1}$	Wach et al. (2013), LL chondrite NWA 4560
Particle Thermal Conductivity	1.5 $W{m}^{-1}{K}^{-1}$	Opeil et al. (2012); Flynn et al. (2017); assuming microporosity of 0.082.
Gravity	8.6e-5 $m{s}^{-2}$	Scheeres et al. (2006)
Young’s Modulus	$7.80{10}^{10}$ $Pa$	Schultz (1995); Gundlach and Blum (2013)
Poisson’s Ratio	0.25	Schultz (1995); Gundlach and Blum (2013)
Surface Energy	0.02 $J{m}^{-2}$	Gundlach and Blum (2013)
Emissivity	1	Gundlach and Blum (2013); Sakatani et al. (2018) (assumed)

Equations

$$\delta =\sqrt{\frac{k}{\rho c}}\sqrt{\frac{P}{\pi }}=\frac{\Gamma }{\rho c}\sqrt{\frac{P}{\pi }}$$

$$F(\lambda ,t)=f_{rock}(\lambda ,t)RA+f_{regolith}(\lambda ,t)(1-RA)$$

$$\chi^{2}_{r}=\frac{1}{\nu}\sum^{N}_{i=1}\frac{\big{(}O-M(\textbf{x})\big{)}^{2}}{\sigma^{2}}$$

$$\chi^{2}_{r}~{}\in~{}\big{[}~{}\chi^{2}_{r,min},~{}\chi^{2}_{r,min}+\frac{\sqrt{2\nu}}{\nu}~{}\big{]}$$

$$f_L(\mathbf{x})=\sum _{i=1}^L{\mathbf{\beta }}_i{g}_i(\mathbf{x})=\sum _{i=1}^L{\mathbf{\beta }}_{i}G({\mathbf{a}}_i,b_i,\mathbf{x})$$

$$G({\mathbf{a}}_i,b_i,\mathbf{x})=g({\mathbf{a}}_i^T\mathbf{x}+b_i)$$

$$h(y)=\frac{2}{1+\exp (-2y)}-1$$

$$ϵ({\mathbf{a}}_i,b_i)=\frac{1}{N}\sum _{i=1}^N(f_L({\mathbf{x}}_i)-\mathbf{F}({\mathbf{x}}_i){)}^2$$

$$\pi (\mathbf{x}|\mathbf{m})=\frac{{\pi }_{pr}(\mathbf{x})L(\mathbf{m}|\mathbf{x})}{\pi (\mathbf{m})}$$

$$\pi (\mathbf{m})={\int }_x\pi (\mathbf{x},\mathbf{m})dx={\int }_x{\pi }_{pr}(x)L(\mathbf{m}|\mathbf{x})dx$$

$${\Gamma }_{average}={\Gamma }_{rock}RA+{\Gamma }_{regolith}(1-RA)$$

$$\rho =\sqrt{f\frac{ln(1-2S)-2S(S-1)}{4S(S-1)}}$$

Full text

††journal: Icarus

1 Introduction

Surfaces of asteroids carry information about the composition, geology, formation and evolution of these bodies (Murdoch et al. 2015, and references therein). All the asteroids observed so far have been found to be covered in regolith, which is a layer of fragmented and unconsolidated material (e.g., McKay et al. 1991) on an airless world. Small asteroids typically display substantial regolith with a variety of features. NASA’s NEAR-Shoemaker mission to asteroid (433) Eros imaged debris aprons, fine-grained “ponded” deposits, talus cones, and bright and dark streamers on steep slopes indicative of efficient downslope movement of the regolith (Robinson et al. 2002). JAXA’s Hayabusa mission found a considerable amount of regolith distributed nonuniformly on the surface of the small ($\sim$300-m diameter) near-Earth asteroid (25143) Itokawa (Miyamoto et al. 2007), which led to the question of how the regolith is segregated. In the case of larger airless bodies like the Moon, size-sorted regolith is locally concentrated, typically covering a bedrock. In the case of small bodies there is evidence of unconsolidated surface materials extending several hundred meters deep (Robinson et al. 2002; Vincent et al. 2012) and, in many cases, it extends throughout the asteroid, leading to the concept that they are rubble piles held together by gravity and cohesive forces (Richardson et al. 2002). These cohesive forces, including van der Waals interactions between the constituent grains, could play important roles in shaping the geomorphology of small rubble-pile asteroids (Scheeres et al. 2010; Rozitis and Green 2014) and their responses to collisions (Asphaug et al. 2015). The presence of unconsolidated surface material, however, is not diagnostic of a rubble pile configuration; a thick layer of regolith has been recorded also on differentiated or partially differentiated bodies, such as (4) Vesta (Keil 2002, and references therein) and (21) Lutetia (Vincent et al. 2012; Weiss et al. 2012).

The properties of asteroid regoliths and in general the geomorphology of the surfaces of these bodies are diverse, and vary with their sizes and compositions. Only a few asteroids have been examined in close physical detail, by spacecraft, and a systematic exploration is being conducted by Earth-based telescopes and radars (Murdoch et al. 2015, and references therein). A few general conclusions, however, can already be drawn. It appears that large asteroids with sizes in the range of some hundreds $km$ are found to have very fine regolith grains of the order of tens to some hundreds of $\mu m$, while $km$-sized and smaller asteroids have coarser regoliths, whose typical grains can reach $cm$ or $dm$ sizes (Vernazza et al. 2012; Gundlach and Blum 2013). Itokawa regolith is gravel-sized in texture, while Eros regolith is dusty. This is not unexpected (Hartzell and Scheeres 2011), since larger asteroids have more gravity and can hold onto more of their near-surface fragmentation products.

Regolith evolves in time in response to processes such as impacts, thermal cycling and space weathering. Despite their low gravity, it is believed that most asteroids form from the accumulation of catastrophic disruption fragments (Benz and Asphaug 1999; Michel and Richardson 2013) resulting in the creation of surfaces composed of boulders and smaller remnants. There is also indication that boulders are subsequently broken up into the smaller particles constituting the regolith by micrometeoroid impacts (Hoerz et al. 1975; Hörz and Cintala 1997; Basilevsky et al. 2013), a fraction of them held by the asteroid and a fraction of them escaping the asteroid. Boulders are broken more gently by thermal fatigue cracking, caused by the numerous temperature variations between day and night (Delbo et al. 2014; Molaro and Byrne 2011; El Mir et al. 2016).

A realization of planetary sciences of the last $\sim$20 years is that the properties of the surfaces (and therefore of regolith) affect the orbital evolution of asteroids: the finite thermal inertia of the regolith causes the asteroid to emit more thermal photons in the afternoon compared to the morning emission. This imbalanced emission produces an acceleration that has a components along the orbital velocity of the asteroid, causing body’s orbital semi-major axis to change with time, i.e., the so-called Yarkovsky effect (Vokrouhlický 1998; Bottke et al. 2006; Vokrouhlický et al. 2015). A related mechanism, known with the acronym of YORP effect, is related to the force torque due to reflection and emission of visibile and thermal photons from the asteroid surface. YORP produces a variation of the rotation period and direction of the spin axis. Interestingly, the tangential component of the YORP effect is strongly affected by the amount of rocks with respect to that of fines (Golubov and Krugly 2012).

Knowledge of surface temperatures, determinations of thermal inertia, and the presence and abundance of visible rocks and their size distributions, are crucial in defining the thermal and mechanical environments on an asteroid’s surface. These properties are used to estimate “sampleability” of a surface, to plan near-surface, landing, and sampling operations of active space missions, such as JAXA’s Hayabusa 2 (Abe et al. 2012) and NASA’s OSIRIS-REx (Lauretta et al. 2014), and, in the future, JAXA’s MMX mission to Phobos, and human interactions with asteroids and the small moons of Mars. As an example, the sampling device of OSIRIS-REx, TAGSAM, which is designed to collect and return at least 60g of the regolith of the asteroid (101955) Bennu, cannot collect grains larger than $\sim$2 cm, which illustrates how accurate estimates of grain size in asteroid regolith influence sampleability.

The present work focuses on the development of a new method to determine regolith and rock physical properties on airless bodies, coupling machine learning techniques and Bayesian inversion of observed infrared fluxes. In the following Section 2 we introduce the state-of-the-art techniques and the proposed new methodology. In Section 3.1 we invert simulated thermal infrared observations by NASA’s OSIRIS-REx. In Section 3.2 we interpret thermal infrared observations of the near-Earth asteroid Itokawa to retrieve surface properties (namely: surface roughness, rock abundance, thermal inertia of the rock and regolith components). We conclude the paper with a discussion of the potentialities of the method and future applications.

2 Methods

2.1 Infrared remote sensing of planetary surfaces

A powerful technique to estimate the properties of regoliths on airless bodies consists in the analysis of the heat emitted in response to the changing diurnal insolation. This heat is typically measured by means of remote-sensing observations in the thermal infrared, which allows the determination of physical quantities such as the thermal inertia ($\Gamma=\sqrt{k\rho c}$, where $k$ is the thermal conductivity of the surface material, $\rho$ is the bulk density of the granular material, and $c$ is the mean specific heat), the degree of roughness, and the effective grain size of the regolith (Gundlach and Blum 2013; Delbo et al. 2015; Davidsson et al. 2015). In particular, surfaces with low thermal inertia respond to change of the illumination with instantaneous changes of their temperature, causing large temperature differences between day and night (Harris and Lagerros 2002; Delbo et al. 2015, and references therein), whereas in the case of higher values of the thermal inertia, the surface responds more slowly to changes of the insolation, resulting in smoother diurnal temperature curves, with smaller differences in temperature between the day and the night. The granularity, porosity, degree of compaction, and lateral and vertical variations in these parameters within the field of view of the infrared instrument also affect the effective thermal inertia of the regolith. Regolith temperatures are influenced by the presence and amount of large rocks, in addition to factors such as the thermophysical properties of the regolith fines, latitude and local slopes, and radiative heating from adjacent crater walls.

The infrared radiance spectrum of a multi-material surface cannot be matched to a single blackbody at a single temperature but rather is a combination of Planck curves in radiance space. The infrared flux of the surface is contributed by the emission from both material (rock and regolith), weighted according to their abundance. It is common experience that rocks have higher thermal inertia of fines, such as sand. Rocks remain warmer than regolith during the night, due to their higher thermal inertia, so that mutual heating by exchange of radiation between warm and cold surfaces may increase regolith temperatures in rocky areas. Our definition of a rock is every boulder that is larger than the thermal skin depth of the surface, defined as a function of the average thermal conductivity $k$, bulk density $\rho$, specific heat $c$ of the surface layer, and the spin rotation period $P$:

[TABLE]

where $\Gamma$ is the thermal inertia of the surface material. Under this assumption, the infrared flux is the linear superposition of two spectra: the emission of fine regolith and the emission of solid rocks. The coefficient of the sum that multiplies each spectrum is the areal rock abundance (hereafter called rock abundance, or RA), for the rock and the regolith components, $RA$ and $(1-RA)$ respectively, i.e.,

[TABLE]

where $\lambda$ is the wavelength and $t$ is the epoch of the observation (function of observation geometry and illumination conditions). Assuming the above scheme, maps of the rock abundance of lunar terrain were calculated using the discrepancy between nighttime brightness temperatures from Diviner’s thermal channels on board of the Lunar Reconnaissance Orbiter (LRO), for fields of view that contain mixtures of warm rocks and cooler regolith. Diviner’s short-wavelength thermal infrared detectors measure higher temperatures than the longer-wavelength detectors do for a given scene containing both rocks and fine-grained materials. Through the analysis of three separate wavelengths, Bandfield et al. (2011) solved simultaneously for: the areal fraction of each scene occupied by exposed rocks $\sim$1 m and larger; and the temperature of the fines, using modeled rock temperatures as inputs. Theoretically, also the rock temperature could have been estimated given the availability of three channels; practically, the solution of inverse problem is not unique. Because of that, lunar rock temperatures were modeled a priori assuming the properties for vesicular basalt described in Horai and Simmons (1972). The thermal inertia was set to be 1570 $Js^{-1/2}K^{-1}m^{-2}$  at 200 K; the rock had albedo 0.15 and emissivity 0.95 and was modeled as an infinitely thick, level, and laterally continuous layer. The analysis of lunar surface resulted in a very low rock abundance (0.4% global average within $\pm 60^{\circ}$ latitude) which confirms the dominance of fine regolith due to micrometeorite bombardment and space weathering, expected also from thermal inertia measurements.

On Mars, infrared remote sensing data have been used for decades to make interpretation of rock abundance (chapter 18 in Bell 2008, and references therein, for a review), starting in earnest with analyses of the martian surface by infrared radiometers aboard the Mariner probes (Neugebauer et al. 1971), the IRTM (Infrared Thermal Mapper) instruments on the Viking orbiters (Kieffer et al. 1976). Although the physics of heat flow through regolith on Mars differ from airless bodies in that conduction by gas in the pore spaces is quite importance, the same general techniques to interpret infrared emission data may still be applied. As an example, the rock abundance and the physical properties of the regolith fines on Mars have been analyzed on a global scale using IRTM data (Christensen 1986; Jakosky and Christensen 1986) and Thermal Emission Spectrometer (TES) (Mellon et al. 2000; Nowicki and Christensen 2007) by comparing short and long wavelength channels (7 and 20 $\mu m$. and 9 and 30 $\mu m$ for IRTM and TES respectively) for brigthness temperature derived using the KRC 1D model (Kieffer 2013, for a review). Global thermophysical maps at 100 m/pixel resolution have also been derived from THEMIS (Thermal Emission Imaging System) data (Fergason et al. 2006b), although rock abundance has not been estimated due to the limited number of spectral bands (i.e., lack of a high-resolution long-wavelength channel). In-situ, high-resolution thermal analyses have also been conducted on Mars by mini-TES (Fergason et al. 2006a) and REMS (Hamilton et al. 2014) on board the MER rovers and the MSL rover, respectively. In general, these analyses are in good agreement with their orbital counterpart.

On asteroids, the thermal inertia of regolith and rock components, and relative rock abundance, have never been measured from infrared observations. Infrared data of asteroids are usually acquired from ground-based telescopes. This does not allow to observe the night-side of the asteroid, whose infrared flux is mostly dominated by the emission from rocks. The observed fluxes are matched to surface properties by means of a thermophysical model of the asteroid, which maps input parameters (surface properties, observation geometry, illumination conditions, instrument performances, shape model of the asteroid) into an infrared flux (Delbo et al. 2015). The canonical approach is to generate a large number of cases, i.e., entries of the type {surface properties, infrared flux}, and then look for the solution which minimizes the residual error between the predicted and the observed fluxes. Common practice is to look for the minimum of the reduced $\chi^{2}$, defined as:

[TABLE]

where ”O” is the observed flux, ”$M(\textbf{x})$” is the modeled flux corresponding to the surface properties (array x), $\sigma$ is the measurement error on the observed fluxes, N is the number of observations and $\nu$ = N - DOF, where DOF is the number of degrees of freedom (the number of parameters to be estimated). According to the $\chi^{2}_{r}$ statistics, all the solutions with $\chi^{2}_{r}$ in the range:

[TABLE]

are acceptable. If the dimensionality of the problem is enhanced, however, the canonical $\chi^{2}_{r}$ approach is unsuccessful. It is indeed computationally expensive to map the parameter space at a resolution sufficient to rule out local minima. To appreciate this, we can think of a model with four inputs (e.g., roughness of the surface, rock abundance, thermal inertia regolith and rock) used to fit infrared fluxes acquired by a nearby spacecraft observing both the day- and night-side of the asteroid. A primary source of computational complexity arises from modelling the surface roughness. This is done by carving craters on the surface for different values of semiaperture angle and crater surface density; the global effect is an effective surface tilt, also called the macroscopic roughness angle or Hapke angle (Hapke 1984). The calculator needs to compute a full heat exchange solution (i.e., considering all the view factors within the concavities of both the shape model and craters), as the usual Lagerros (1998) approximation cannot be used to fit the night-side observations. With this simulation setting, the asteroid attains simulated thermal equilibrium with the environment in about 3.5 minutes (on a single processor 2.8 GHz Intel Core i7, for a shape model of the asteroids with 2292 facets). Let us sample the parameter space using 25 points for roughness, 25 for the thermal inertia of regolith and 25 for that of the rock components; the solutions can be linearly combined to study the effect of various rock abundances (Equation 2). The exercise requires to run 1250 full heat exchange thermal simulations (about 3 days). This sampling of the parameter space, however, is too coarse to identify local minima or saddle points of the $\chi^{2}_{r}$ function. Therefore, building a look-up table can be misplaced effort, unless the functional relationship between parameters and infrared thermal response is learned and generalized in some way. As we show in the following, the dataset can be instead used to train a neural network capable of predicting the infrared fluxes; as opposed to the “parent” model, this tool can enable a fine – and fast – mapping of the parameter space into the corresponding infrared flux within a known level of accuracy .

2.2 Surrogate models for thermophysical applications

The new approach consists in approximating the thermophysical function $\boldsymbol{y}=\boldsymbol{F}(\boldsymbol{x})$ ($\boldsymbol{x}$: surface properties; $\boldsymbol{y}$: infrared flux) using neural networks trained on a dataset of thermal simulations of the type: {surface properties; infrared flux}. After training, the networks can predict an infrared flux at at highly reduced computational time with respect to the “parent” model (running in less than a second versus minutes on a single processor 2.8 GHz Intel Core i7), thus enabling Markov Chain Monte Carlo (MCMC) Bayesian inference of surface properties from observed fluxes. MCMC requires thousands of runs of the forward model to sample the unknown posterior distributions – thus the need for a fast and accurate predictor, i.e., a trained neural network.

A dataset of thermophysical simulations is generated using the TPM code (Delbo et al. 2015), for different combinations of surface properties, given the shape model of the asteroid, the observational geometry and the illumination condition; 2. 2.

The dataset is used to train, validate and test a neural network (Section 2.2.1). The trained network is a surrogate model which predicts the asteroid’s infrared flux at a highly reduced computational speed with respect to the ”parent” thermophysical model (i.e., less than a second on a single processor 2.8 GHz Intel Core i7). 3. 3.

A series of blind tests is performed to approve the methodology. Each blind test consists in the generation of a case TPM flux by one of the co-author; the case flux is sent to another co-author, who does not know the values of the surface properties used to generate the case. The second co-author performs a Markov Chain Monte Carlo Bayesian inversion of the case TPM flux using the surrogate model (Section 2.2.2). The posterior distributions of the surface properties are sent back to the first co-author – the only person who knows the values of the surface properties of the case flux – who verifies the success of the inversion. 4. 4.

Once the methodology is validated, we use it to invert observed infrared flux in order to characterize the surface properties of the body, e.g., asteroid (25143) Itokawa in Section 3.2.

2.2.1 Neural Networks

Neural Networks (NN) are mathematical models that take inspiration from the fundamental structure of the brain (Schmidhuber 2015). NNs are able to learn (i.e., improve the performance of a specific tasks) from data, commonly subdivided in training, validation and testing subsets. NNs are very successful in modeling the functional relationship between inputs and outputs. Such relationship is generally represented by a set of examples (training set) and learned during the training process. NNs can be employed as physically-based fast surrogate models for parameter estimation, replacing computationally expensive numerical models.

NNs consist of many mathematical units called neurons that are connected via synapses. Neurons are organized by layers and communicate in a parallel fashion through weights that represent the strength of the corresponding synapses. A typical shallow network employed to include surrogate models, comprises one input layer, one hidden layer and an output layer. The hidden layer is assumed to have a specified number of neurons $L$ which are the basic processing units for the network. The overall process begins with a summation of each input with the correspondent weights (synapses) and then further processing by an activation function. In regression problems, the overall NN output function is typically represented as follows:

[TABLE]

where $\boldsymbol{x}\in\ \mathbb{R}^{d}$ and $\boldsymbol{\beta}_{i}\in\mathbb{R}^{m}$. Neurons can usually have different activation functions $g(\cdot)$. Additive nodes with activation functions have the following structure:

[TABLE]

An activation function which works well in shallow neural networks is the tanh-sigmoid function

[TABLE]

where $\boldsymbol{a}_{i}\in\mathbb{R}^{m}$ and $b_{i}\in\mathbb{R}$. The weights $\boldsymbol{a}_{i}$ and biases $b_{i}$ are usually determined during the training process which implies minimization of a loss function. For regression problems which are usually solved to generate surrogate models, the typical loss function is the Mean Square Error (MSE), i.e.:

[TABLE]

where $\{(\boldsymbol{x}_{i},\boldsymbol{F}(\boldsymbol{x}_{i})\}$ is the associated training set. In our application, $\boldsymbol{x}_{i}$ is the vector of surface parameters and $\boldsymbol{F}(\boldsymbol{x}_{i})$ is the corresponding infrared fluxes. The accuracy of the network is also quantified by the correlation coefficient between predictions and labels. MSE and correlation coefficient are indicator of the reliability with which the network ”mimics” (approximate) the ”parent” model (e.g., the TPM code).

The dataset is split into a training set (70%), validation set (15%) and testing set (15%). The training set is employed for proper network training, which involves the fitting of the network parameters (weight $\boldsymbol{a}_{i}$ and biases $b_{i}$) by minimizing a cost (loss) function (Equation 8). Loss minimization is implemented by standard backpropagation and stochastic gradient descent which are the most common modern approach to training (Schmidhuber 2015). The overall process is not a simple interpolation of the training set, but rather it involves the search for classes of parametric functions (i.e., the neural network) that globally fit the data (generalization). It is indeed undesirable for any algorithm to match exactly the training set (e.g., through spline interpolation), as the method would then perform poorly on unseen data. The validation set is employed to determine the optimal network hyperparameters (e.g., the learning rate). Finally, the testing set is employed during the training process to evaluate the behavior of the MSE on an unseen ensemble of data, for an independent assessment of the generalization capabilities of the network. Importantly, both validation and testing sets are employed to protect against overfitting of the training set (Bishop et al. 1995). Indeed, properly trained networks ensure that the data in the validation and the testing sets follow the same probability distribution of the data in the training dataset. Training occurs as a succession of training epochs. At each epoch, the network processes the full set of training data (i.e., it computes the response for all of the data in the training set). Consequently, backpropagation of the residuals between targets and prediction is executed to update the network parameters (weights and biases) that reduce the MSE. Additionally, at every epoch, the MSE for validation and testing is computed; the training is completed (convergence) if the validation MSE does not decrease further after 6 consequent epochs.

2.2.2 Markov Chain Monte Carlo Bayesian inversion

The general goal of inverse problems is to use a physical or surrogate model to retrieve properties that characterize a physical system. In the Bayesian approach, the parameters to be retrieved are treated as random variables. The Bayesian approach to inverse problem generally yields a probability distribution (a-posteriori or posterior distribution) for each of the retrieved variables. Such distribution is obtained by virtue of the Bayes’ theorem which combines a-priori distribution and likelihood given the observed data to obtain an estimation of the parameters distribution (Stuart 2010; Aster et al. 2011). Indeed, once observations are collected, one can compute the conditional probability distribution $L(\boldsymbol{m}|\boldsymbol{x})$, i.e. given a certain set of parameters to be retrieved $\boldsymbol{x}$, the set of observations $\boldsymbol{m}$ are observed. As mentioned above, the Bayes’ theorem is employed to combine such $likelihood$ conditional probability distribution, $L(\boldsymbol{m}|\boldsymbol{x})$, with the prior distribution $\pi_{pr}(\boldsymbol{x})$ to compute the conditional posterior distribution for the model parameters $\pi(\boldsymbol{x}|\boldsymbol{m})$. Solving inverse problems within the Bayesian framework relies on computing the conditional posterior probability distribution, i.e.:

[TABLE]

where $\pi(\boldsymbol{m})$ is the probability distribution of the observed data:

[TABLE]

Generally, the inverse problem is regularized by using the prior distribution $\pi_{pr}(\boldsymbol{x})$ which encodes prior knowledge or belief about the parameters. As stated in Stuart (2010), since all variables are assumed to be random, they reflects the uncertainty of the true values and, importantly, the degree of uncertainty is encoded in the probability distribution of the parameters to be retrieved. Typical estimates using $\pi(\boldsymbol{x}|\boldsymbol{m})$ are the 1) Maximum A-Posteriori (MAP) distribution and 2) Mean of the posterior distribution. If the posterior distribution is normal, the MAP and posterior mean will coincide (Stuart 2010; Aster et al. 2011). Generally, analytical evaluations of such integrals are quite impossible. When $n$ is large, traditional numerical quadrature schemes are not practical and one has to resort to Monte Carlo (MC) integration methods. MC algorithms are applied by sampling $N$ points on $\pi(\boldsymbol{x}|\boldsymbol{m})$, ${\boldsymbol{x}_{k},k=1,2,...,N}$. In this work, a Monte Carlo Markov Chain (MCMC) approach is employed for sampling. There are many algorithms for MCMC; one of the most widely used is the Metropolis–Hastings (MH) algorithm that performs a Markov chain with a specified limiting distribution. Other used algorithms for MCMC are the Adaptive Metropolis-Hastings (AM), Delayed rejection (DR), or their combination called DRAM (Haario et al. 2006).

We demonstrate the potentialities of the proposed method by applying it to two case studies. The first application (Section 3.1) discusses the inversion of simulated infrared fluxes of asteroid (101955) Bennu as acquired by the OSIRIS-REx mission. The second application (Section 3.2) is about the inversion of observed infrared fluxes of asteroid (25143) Itokawa. In both cases, the method is able to accomplish the inversion and to provide insights about the surface properties of the bodies.

3 Results and Discussion

3.1 Blind tests on simulated infrared flux of (101955) Bennu

Once at the target asteroid (101955) Bennu, the OSIRIS-REx mission will perform a detailed survey of the surface in order to constrain its properties and select the sampling site (Lauretta et al. 2017). OSIRIS-REx requirement on sampling is to collect regolith not hotter than 350K, and TAGSAM should be used to collect fine regolith, e.g., grain size less than 2 cm. Furthermore, OSIRIS-REx should avoid rocks (which are hazardous) during the sampling phase. These mission requirements will constrain the choice of the latitude of the sample selection area, the local time, and the arrival date at the asteroid (Delbo et al. 2015). The OSIRIS-REx Thermal Emission Spectrometer (OTES, Christensen et al. 2018) will collect the infrared fluxes emitted by the surface of the asteroid covering the spectral range $5.71-100~{}\mu m$ $(1750-100~{}cm^{-1})$ with a spectral sample interval of $8.66~{}cm^{-1}$. The thermal properties of the surface will inform the team about the sampleability of the different regions of the asteroid. The abundance of rock and surface roughness will be used by the mission engineers to define the hazard associated with sampling a certain region; the thermal inertia of the regolith can be converted in average grain size (Gundlach and Blum 2013; Sakatani et al. 2017).

3.1.1 Surrogate models for OTES acquisition

In the case of Bennu, we apply the methodology described in Section 2 as follows. A 4-dimensional dataset of infrared fluxes (4 surface properties: surface roughness, thermal inertia of the regolith, thermal inertia of the rock, rock abundance) is populated with simulated OTES-like acquisitions of the infrared flux from (101955) Bennu, using the TPM (Delbo et al. 2015). We work with two shape models: a single triangular equatorial facet, and the full shape model of the asteroid by Nolan et al. (2013). The shape model is represented by a mesh of 2292 triangular facets. Facets can have any number of vertices, in principle; however, the temperature and flux of the facet depends only on its area and normal direction. On the surface, we use the full heat exchange solution (i.e., considering all the view factors within the concavities of both the shape model and craters), as the Lagerros (1998) approximation is not suitable to fit the night-side observations. In the simulations, the observer is always looking nadir, i.e. the unity vector from the center of the asteroid to the observer and the normal to the surface have scalar product equal to one, and it is located at infinite distance on the equatorial plane. The heliocentric distance is assumed to be 1 au. The model output are 7 fluxes in the wavelength range between 5 and 50 $\mu m$, each corresponding to one of the different O-REx’s equatorial stations of the detailed survey (03:20 AM, 06:00 AM, 10:00 AM, 12:30 PM, 03:00 PM, 06:00 PM, 08:40 PM). We model the surface with surface roughness values between 4 and 62 degrees (in terms of Hapke angle, Hapke 1984). We assume constant thermal inertia on the surface, for values varying between 0 and 1000 $Js^{-1/2}K^{-1}m^{-2}$. The fluxes corresponding to values of thermal inertia below 500 $Js^{-1/2}K^{-1}m^{-2}$ are assumed to be emitted by a surface covered in regolith ($f_{regolith}$). The fluxes corresponding to values of thermal inertia above 500 $Js^{-1/2}K^{-1}m^{-2}$ are instead assumed to be emitted by a surface covered in boulders ($f_{rock}$). The final dataset is composed by all the possible combinations of $f_{regolith}$ and $f_{rock}$, weighted in terms of rock abundance (see Equation 2, Table 1). Building the dataset requires to run 1200 TPM simulations (12 different values of roughness for 100 different values of thermal inertia). Each simulation runs in 0.602 seconds for a single facet and 3.5 minutes for the disk-integrated case, corresponding to about 12 minutes and 3 days for the overall building of the dataset, respectively. The fluxes for regolith and rock are then linearly combined in terms of rock abundance (Equation 2).

The dataset of infrared fluxes is used to train seven (7) surrogate models (one for each of the equatorial station) which mimic the acquisition of infrared fluxes at Bennu by OTES. Each of the surrogate model is derived from the deployment of shallow neural networks on the 4-dimensional dataset of infrared simulated fluxes (Table 1). We refer to Section 2.2.1 for the training, validation and testing procedure, and to Table 2 for the architecture and performance of the networks. The whole training procedure requires a computational effort of about 10 minutes. The trained networks run in less than a second. All the run time are computed on a single processor 2.8 GHz Intel Core i7.

Each neural network maps the flux into a 4-Dimensional parameter space (surface roughness, thermal inertia of the regolith, thermal inertia of the rock, rock abundance). As an example, Figure 1 is a contour plot of the peak value of the infrared flux, as predicted by the neural network and observed on the night-side (equatorial station: 3:20 am). The free parameters are surface roughness and the average thermal inertia in lieu of rock abundance (RA), e.g.,

[TABLE]

where $\Gamma_{rock}$ and $\Gamma_{regolith}$ are the thermal inertias of the rock and the regolith. In the example of Figure 1, $\Gamma_{rock}$ and $\Gamma_{regolith}$ are kept constant – and equal to 1000 and 150 $Js^{-1/2}K^{-1}m^{-2}$, respectively – to allow a 2-D visualization of the parameter space. The peak flux increases with the average thermal inertia because rocks remain warmer on the night-side and, the higher the rock abundance, the higher the average thermal inertia, the higher the peak flux. An increase in the flux is also observed with roughness. Since roughness is modeled by carving craters on the surface, we observe an increase in flux because of self-heating by re-absorption of thermal radiation from other parts of the crater (Figure 6.8 in Delbo 2004).

3.1.2 MCMC Bayesian inversion

The reduced computational time of the neural networks and their capability of generalizing a response for predictors in between the training nodes make these models statistically invertible (Section 2.2.2). The methodology is here blind-tested against a set of case fluxes (step 3 in Section 2.2). The case flux is generated using the TPM code and mimic a ”real” observed flux that could be acquired by OTES. The fluxes’ Signal-to-Noise Ratio (SNR) is modeled such that it resembles the real performance of the spectrometer (Christensen et al. 2018, and priv. comm.). We use the trained networks as forward models in a Markow Chain Monte Carlo Bayesian inversion of the simulated infrared fluxes. For this exercise, the inversion is not informed about any a-priori information (i.e., we use uninformative a-priori distributions). In real application, it is reasonable to have some a-priori information on the parameters, e.g., the OSIRIS-REx Visible and Infrared Spectrometer (OVIRS, Reuter et al. 2018) would also provide information about the material on the surface, thus providing information about the thermal inertia. Using proper a-priori distributions improves the accuracy of the inversion and helps achieving convergence. As we show below, however, for this application the method is able to correctly perform the inversion also without any a-priori information.

The blind test for the triangular equatorial facet uses a case TPM flux corresponding to surface properties: $\bar{\theta}=46^{\circ}$, $\Gamma_{regolith}=60$ $Js^{-1/2}K^{-1}m^{-2}$, $\Gamma_{rock}=820$ $Js^{-1/2}K^{-1}m^{-2}$, RA = $25.0\%$ (top panel in Figure 2). The MCMC inversion of the blind test provides the posterior distributions in Figure 2, bottom panel. The same blind-test procedure, but for the disk-integrated flux, used a case flux corresponding to surface properties: $\bar{\theta}=34^{\circ}$, $\Gamma_{regolith}=50$ $Js^{-1/2}K^{-1}m^{-2}$, $\Gamma_{rock}=850$ $Js^{-1/2}K^{-1}m^{-2}$, RA = $30.0\%$ (top panel in Figure 3); the posterior distributions from the inversion are in Figure 3, bottom panel. In both cases, the blind tests are successful as the values of the surface properties for the case fluxes belong to the posterior distributions; the difference between the reference and the statistical mean is less than 1–$\sigma$. An exception is the thermal inertia of the rock abundance for the single triangular facet, for which the reference value is within 2–$\sigma$ of the estimated mean value.

To test the robustness of the method, we additionally perform an inversion using the observations relative to the daytime stations only (10:00 AM, 12:30 AM, 3:00 AM; disk integrated case). On the night-side, the infrared flux is mainly defined by the rocky material because of its inherently higher thermal inertia and no shadow is present. Therefore the nighttime observations are important to constrain the presence of rocky material and to get a more complete characterization of the surface. The posterior distributions for this case are the dashed curves in Figure 3, bottom panel. The inversion of the disk-integrated case using only the daytime stations is still successful. The main difference with respect to the full inversion is observed on the posterior distribution of the thermal inertia of the rock; this result was expected as the rocks appear brighter in the infrared on the night-side (as it cools down slower than the regolith). We conclude that the inversion method is able to retrieve the correct surface properties both for a single facet and for the disk-integrated case without using any a-priori information. The (small) deviations between the reference and retrieved values are attributed to the effect of noise from the OTES instrument.

3.2 Inversion of observed infrared fluxes of (25143) Itokawa

The near-Earth asteroid (25143) Itokawa was characterized in great detail by both ground-based observations and the JAXA’s Hayabusa mission. Here we use the proposed methodology to invert the available mid-infrared photometric observations and we compare the results with pre and post-Hayabusa findings, e.g., surface roughness, thermal inertia and grain size of the regolith (Kitazato et al. 2008; Müller et al. 2014; Sakatani et al. 2017). The mid-infrared observations are summarized in Table 3. These consist in 25 ground-based observations and five measurements from the Japanese infrared astronomical satellite AKARI (Müller et al. 2005, 2014). All the epoch of observations are the time of emission by the asteroid. The position vectors of the asteroid and the Earth have been queried to the JPL Horizons Ephemeridis system (https://ssd.jpl.nasa.gov/horizons.cgi) and converted in the asteroid-centric reference frame.

3.2.1 Thermal surrogate model of (25143) Itokawa

The steps for the design of the surrogate model in the case of (25143) Itokawa are analogous to those for (111995) Bennu (Section 3.1.1), with some proper adaptations which are described in the followings. For (25143) Itokawa we use the TPM code to populate two datasets: a 2-D dataset and a 4-D dataset, Table 4. The infrared fluxes in the 2-D dataset are modeled with two surface properties: Hapke angle and average thermal inertia of the surface. The infrared fluxes in the 4-D dataset are modeled with four surface properties: Hapke angle, thermal inertia of the regolith, thermal inertia of the rock, rock abundance. For the 4-D dataset, the thermal infrared radiation is computed as superposition of the contribution from the regolith and the rock (Equation 2). Each of the infrared fluxes is generated using the shape model by (Gaskell et al. 2008) at a resolution of 4096 facets. The surface roughness is represented using craters and the fluxes have been computed using the Lagerros’ approximation (Lagerros 1998), as no night-side observations are available. Each run simulates the acquisition of the flux corresponding to the wavelength and epoch of the observation in Table 3; the asteroid is observed from Earth. The run time for a simulation is 15.2 seconds if the Lagerros (1998) approximation is employed.

Each surrogate model is tailored to the specific observation geometry and spectral range and derives from the deployment of shallow neural networks on the two datasets described in Table 4. For each observation in Table 3, a neural network is trained, validated and tested on the 2-D dataset, composed by 720 predictions of infrared fluxes corresponding to different combinations of surface properties. Therefore, the total computation cost to populate the 2-D dataset is about 4 days ($720\times 30\times 15.2$ seconds). The same procedure is followed for the 4-D dataset, which is composed by 1200 simulations per observation (total run time of 5 days). We refer to Section 2.2.1 for the training, validation and testing procedure. The whole training procedure requires a computational effort of about 1 hour. Table 5 indicates the metrics of performance of the networks, which run in less than a second. All the time are computed on a single processor 2.8 GHz Intel Core i7.

3.2.2 Bayesian inversion of the 2-D dataset

The 2-D surrogate models (inputs: surface roughness and average thermal inertia) are used as forward models in the Markov Chain Monte Carlo Bayesian inversion of the observed infrared fluxes in Table 3. The inversion of the 2-D dataset allows constraining the Hapke angle of the surface and the average thermal inertia, which is both contributed by the regolith and the rock components. The posterior distributions for the inversion of the 2-D problem are in Figure 4, top panel. We retrieve an average thermal inertia between 700 and 800 $Js^{-1/2}K^{-1}m^{-2}$, while the distribution of the Hapke angle is bi-modal.

The presence of two peaks for the roughness parameter is intriguing because of what is known about the surface of Itokawa. Analyses of Hayabusa data assessed the presence of both a smoothed terrain and a rough terrain, which have been found to correlate with low- and high- land regions respectively (e.g., Yano et al. 2006; Susorney et al. 2018, and references therein). While the average thermal inertia is well constrained ($\Gamma_{av.}=721\pm 45$ $Js^{-1/2}K^{-1}m^{-2}$), the posterior distribution for the Hapke angle is broad and does not provide any constraint on the roughness. Interestingly, however, our estimate of the Hapke angle is $26~{}\pm~{}18^{\circ}$ (statistical mean and standard deviation). This is consistent both with the results from (Kitazato et al. 2008), who found a photometric roughness of $26~{}\pm~{}1^{\circ}$ in terms of Hapke angle for the near-infrared range 0.76 - 2.25 $\mu m$, and with the average value derived for main-belt asteroids (Muller and Lagerros 1999). Müller et al. (2005) found an average thermal inertia of the surface of Itokawa between 500 and 1000 $Js^{-1/2}K^{-1}m^{-2}$, and therefore proposed a value of 750 $Js^{-1/2}K^{-1}m^{-2}$. The Hapke angle, by contrast, remained unconstrained, and the ”canonical” value for Main Belt Asteroids was assumed (Muller and Lagerros 1999), in terms of surface slope root mean square $\rho$ and a surface crater density $f$: $\rho=0.7$, $f=60\%$. The definition of surface slope root mean square, $\rho$, is that of Lagerros (1998):

[TABLE]

where $S=0.5~{}(1-\cos{\gamma_{c}})$, with $\gamma_{c}$ equal to the semi-aperture angle of the crater. Values of $\rho=0.7$, $f=60\%$ correspond to a Hapke angle of about 26∘.

To better investigate the bimodality, we perform a classical reduced $\chi^{2}$ analysis using the surrogate models in the forward mode, to explore the ($\bar{\theta}$, $\Gamma_{av.}$) phase space. The capability of the networks to generalize the responses for predictors different from the training set enables a fine mapping of the parameter space. After a first iteration, five observations show residuals bigger than 2.5 the standard deviation of the measurement ($\sigma$-s in column 3, Table 3). The observed fluxes could be affected by unknown or/and unreported errors, as the data are an ensemble of observations performed in different conditions and at different epoch. Applying the Chauvenet criterion for data rejection, we remove the five observations (n. 4, 6, 21, 22, 23 in Table 4) and repeat the survey. The phase space $(\bar{\theta},~{}\Gamma_{av.})$ of the selected data is in Figure 4, bottom panel. The $\chi^{2}_{r}$ function shows multiple minima, at different Hapke angles, but at a consistent value of thermal inertia around 750 $Js^{-1/2}K^{-1}m^{-2}$, in agreement with previous findings (Müller et al. 2005). The solution (statistics of the absolute minimum at $\chi^{2}_{r}$ = 1.5) is used to inform the Bayesian inversion. The resulting posterior distributions are in Figure 5; the estimated surface properties are $\bar{\theta}$ = $4~{}\pm~{}1^{\circ}$ and $\Gamma_{av.}$ = $721~{}\pm~{}45$ $Js^{-1/2}K^{-1}m^{-2}$.

3.2.3 Bayesian inversion of the 4-D dataset

The 4-D scheme models the surface in terms of 4 parameters: Hapke angle, thermal inertia of the regolith, thermal inertia of the rock, rock abundance. The last three surface properties contribute to the average thermal inertia of the surface according to their weights (i.e. rock abundance). The 4-D surrogate models are used as forward models in the Markov Chain Monte Carlo inversion of the observed infrared fluxes in Table 3. The results from the inversion of the 2-D dataset inform the inversion of the 4-D model in two ways. A priori, the Hapke angle is informed by the distribution found in Section 3.2.2. A posteriori, if the distributions of thermal inertia of the regolith ($\Gamma_{regolith}$), thermal inertia of the rock ($\Gamma_{rock}$) and rock abundance are linearly combined to give the disk-integrated thermal inertia, i.e., Equation 11, the post-processed distribution must be consistent with the analogous result of the invertion of the 2-D dataset.

We firstly perform a fully unconstrained Bayesian inversion of all the fluxes using only the a-priori information on the Hapke angle (derived from the inversion of the 2-D dataset). Unfortunately, the fully unconstrained approach fails in splitting the regolith’s and rock’s contributions to the thermal inertia. Therefore we need to find admissible ranges for the parameters to support the inversion process. We recall some information from the retrievals of the Hayabusa spacecraft, as discussed below. Using these constraints, we repeat the inversion and successfully retrieve the surface properties.

A-priori information: rock abundance. The size distribution of blocks on the surface of Itokawa, derived from Hayabusa/AMICA images (Mazrouei et al. 2014), allows setting a lower limit for the rock abundance. The integral of the size distribution is 16%, indicating the abundance of rock bigger of around 5 meters (the limiting resolution). This value is therefore a lower limit for the rock abundance, as the infrared flux is contributed also by smaller boulders.

A-priori information: thermal inertia. The thermal inertia of the rock ranges between the lowest thermal inertia recorded for stony meteorites – around 700 $Js^{-1/2}K^{-1}m^{-2}$ (Cold Bokkeveld, CM2 chondrite, Opeil et al. 2010), and 2500 $Js^{-1/2}K^{-1}m^{-2}$. A constraint on the thermal inertia of the regolith can be derived by exploring the relationship between regolith’s average particle size and thermal conductivity.

To interpret the fine-particulate component thermal inertia derived from this work in terms of particle size, we utilize the methodology of Gundlach and Blum (2013). This is involves comparing theoretical values for regolith thermal conductivity at different particle sizes to thermal conductivity values derived from the predicted thermal inertia, leaving regolith porosity as a free parameter. In lieu of the particulate thermal conductivity model developed and utilized by Gundlach and Blum (2013), we instead utilize the improved model by Sakatani et al. (2017) along with updated model tunable parameter values from the experimental work by Ryan (2018). The model by Sakatani calculates the solid and radiative components of the thermal conductivity of particulates as a function of relevant material and environmental parameters (Table 6). Parameters $\zeta$ and $\xi$ are used to tune the radiative and solid conduction components of the model, respectively, in order to fit the model to experimental measurements of particulate thermal conductivity. Sakatani et al. (2017, 2018) provide values for particulates up to approximately 1 $mm$ in diameter. Recently, Ryan (2018) obtained values for the parameters from conductivity experiments for particles up to 1 $cm$ in diameter. The values presented in that work can reasonably be related to particle diameter, $D$, with a power function, where $\zeta=0.149D^{-0.3052}$ and $\xi=0.899D^{0.2588}$. Values for all other model parameters and their sources are provided in Table 6. Nakamura et al. (2011) report that the grains retrieved from Itokawa by Hayabusa are identical to those of thermally metamorphosed LL chondrites. Thus, average material properties for LL chondrites are used where available.

In order to directly relate the thermal inertia to particle size, a range of possibility regolith macroporosity values must be considered. A very loose, random packing of spheres has a maximum porosity of about 0.45 under terrestrial conditions where gravitational forces exceed electrostatic attractive forces between particles (Dullien 2012). Some angular basaltic samples used by Sakatani et al. (2018), however, had porosity values up to 0.66. Kiuchi and Nakamura (2014) develop a relationship between particle size and porosity on small bodies with very small surface accelerations where inter-particle forces might instead dominate. Using their model and the nominal values for silicate particles provided in their paper, and a maximum surface acceleration of 0.086 $mms^{-2}$ for Itokawa (Scheeres et al. 2006), we estimate the maximum porosity for different particle sizes. The computed values for particles with diameters of 0.5 $mm$, 1 $mm$, 5 $mm$, 1 $cm$, and 5 $cm$ are 0.86, 0.83, 0.74, 0.7, and 0.58 (using a cleanliness ratio of unity, which is recommended for particles in space). Given the likely angular nature of the particles on Itokawa, we assume that porosities as high as 0.66 are possible for all particle size ranges. Dense random packings of spheres have a minimum porosity of about 0.36. However, we consider such a dense packing unlikely for angular particles and instead adapt a more modest minimum porosity of 0.4.

The assumption that the regolith can be treated as homogeneous is no longer valid once the particles size exceeds the length of the diurnal skin depth (chapter by Mellon et al. in Bell 2008). For each value of thermal conductivity predicted by the Sakatani et al. (2017) model, a corresponding diurnal skin depth is calculated (Equation 1). If the particle size exceeds the skin depth, these values are rejected. Although the procedure tends to lead to the rejection of larger particle sizes, this does not preclude the possibility of the existence of large particles; it simply means that the model cannot be used to confidently make interpretations within this particle size range. Figure 7 shows the relationship between particle size, porosity and thermal conductivity. The thermal conductivity curves are dotted if the particle size exceeds the calculated skin depth; where the curves are dotted, these are expected to deviate upwards from the apparent power law trend to eventually reach a constant value for thermal conductivity that corresponds to the conductivity of solid rock. The maximum value predicted with confidence by the Sakatani model is to 0.093 $W~{}m^{-1}~{}K^{-1}$ for a porosity $\phi=0.6$; this corresponds to a value of the thermal inertia of the regolith of about 285 $Js^{-1/2}K^{-1}m^{-2}$. We use this value to inform the inversion regarding the upper limit of thermal inertia of the regolith.

As previously discussed, existing models for regolith thermal conductivity can only be reliably used for particle sizes less than a diurnal skin depth. As such, there potentially exist a range of particle sizes larger than the diurnal skin depth, but still small enough that their thermal signature is below that of an infinite bedrock. We choose to avoid interpretations within this range by utilizing the hard boundary for regolith thermal conductivity, but we do acknowledge that this topic needs additional research. In doing so, we make the assumption that the two extremes in thermal cycling (i.e. low amplitude thermal cycling of rock and high amplitude thermal cycling of the finest regolith) will dominate the observed signature, and that the presence of an intermediate particle size regolith within this ”transition” zone would not significantly impact the end results of our interpretations. A smooth transition from rock into regolith should be expected, as the result of comminution by thermal fatigue fragmentation (Delbo et al. 2014; El Mir et al. 2016; Molaro and Byrne 2011) and micrometeoroid bombardment (Basilevsky et al. 2013; Hoerz et al. 1975; Hörz and Cintala 1997).

Using the above constrains, we perform a detailed $\chi^{2}_{r}$ survey on all the observation but those rejected in the inversion the 2-D simulation. The a-priori information available from literature allows focusing the search on admissible sub-regions of the $\chi^{2}_{r}$ hyperspace while looking for the absolute minimum. The solution (statistics of the absolute minimum at $\chi^{2}_{r}$ = 1.6) is used to inform the Bayesian inversion of the 4-D model. Figure 6 shows the posterior distributions for the parameters; the best estimates have statistics: $\bar{\theta}=4\pm 1^{\circ}$, $\Gamma_{regolith}=203\pm 36$ $Js^{-1/2}K^{-1}m^{-2}$, $\Gamma_{rock}=894\pm 122$ $Js^{-1/2}K^{-1}m^{-2}$, RA = $84~{}\pm 9\%$.

3.2.4 Discussion

We compare the results of the inversion of the 4-D dataset to those of the 2-D dataset as well as to in-situ measurement of regolith size. The average thermal inertia of the surface is re-constructed as the combination of the posterior distributions of the thermal inertias (regolith and rock) and the rock abundance. Figure 8 shows the average thermal inertia from the inversion of the 2-D dataset (solid line) versus the post-processed average thermal inertia from the inversion of the 4-D dataset (dashed line). Both curves suggest a value of thermal inertia around $750$ $Js^{-1/2}K^{-1}m^{-2}$ ($\Gamma_{av.,~{}2D}$ = $721~{}\pm~{}45$ $Js^{-1/2}K^{-1}m^{-2}$ and $\Gamma_{av.,~{}4D}=785~{}\pm~{}119$ $Js^{-1/2}K^{-1}m^{-2}$, for the 2-D case and 4-D case respectively). The two curves peak at approximately the same value of $\Gamma_{av.}$ and the estimated values overlap at less than 1–$\sigma$. Both results are consistent with previous findings (Müller et al. 2005). The 4-D solution to average thermal inertia, however, has a broader posterior distribution with respect to the 2-D solution; we attribute this degradation to error propagation.

The value of thermal inertia of the regolith derived in Section 3.2.3 is converted to thermal conductivity, using the mean density and specific heat of LL chrondrites, for different porosity values (horizontal lines in Figure 7). Where corresponding porosity lines intersect, a prediction for particle size is made, following the method by (Gundlach and Blum 2013). Solutions exist for porosity equal to 0.4, 0.5 and 0.6 (squared dots in Figure 7). After error propagation, the mean particle diameter is found to be: $9^{+6}_{-4}$ $mm$ ($\phi=0.4$); $10^{+7}_{-4}$ $mm$ ($\phi=0.5$); $11^{+7}_{-5}$ $mm$ ($\phi=0.6$). These results are in agreement with in-situ findings by Hayabusa (Yano et al. 2006; Kitazato et al. 2008).

The value of thermal inertia of the rock derived at the end of Section 3.2.3 is converted to thermal conductivity, again using the mean density and specific heat of LL chondrites, for a microporosity of 0.082. After error propagation, we retrieve a value of $0.4\pm 0.1~{}W~{}m^{-1}~{}K^{-1}$. This value is substantially lower than the thermal conductivity of coherent LL chondritic rocks, i.e., around 1.5 $~{}W~{}m^{-1}~{}K^{-1}$ (Figure 2 in Opeil et al. 2012); we therefore conclude that the rocks on Itokawa could be fractured.

Finally, the minimum $\chi^{2}_{r}$ for the 4-D model is 1.6. Within 1–$\sigma$ ($\sigma_{\chi^{2}_{r}}=0.3$), this value is the same of the 2-D model and, in this respect, the two models are equivalent. In addition to the 2-D model, however, the 4-D model is able to split the contribution from the regolith and the rock to the average thermal inertia, and therefore is able to inform us about the rock abundance on the surface of the asteroid without any assumption on the thermal inertia of the components (besides the discussed upper and lower limits).

4 Conclusion

Infrared fluxes emitted from asteroids carry information about the properties of the surface. The customary approach to the inverse problem (i.e., backing up surface properties from observed infrared fluxes) is to run thermophysical simulations, in which the asteroid thermal properties are mapped into an infrared flux. The flux is then compared to the observed data following $\chi^{2}$ minimization (see Section 2). Using this approach, past studies (e.g., Harris and Lagerros 2002; Müller et al. 2005; Emery et al. 2006, 2014; Delbo and Tanga 2009; Delbo et al. 2015; Rozitis and Green 2014) have been successful in retrieving size, albedo and/or average thermal inertia of asteroids, but others important properties, such as surface roughness and rock abundance, remained elusive and unconstrained. Many difficulties arise when multiple surface properties are included in the game. The number of simulations needed to densely populate the parameter space grows exponentially with the dimension N of the input array. The objective function to be minimized (e.g., $\chi^{2}$ function in Equation 3), may show local minima or saddle points also for relatively low values of N, as shown in Figure 4 for the interpretation of disk-resolved thermal infrared data of asteroid (25143) Itokawa using the 2-D dataset. In order to split the contribution of regolith and rock to thermal inertia (e.g., Figure 6 for Itokawa), we want to increase the dimensionality of the input array to N = 4. For this task, the use of the TPM code results in a coarse mapping of the parameter space unless we employ significant computational resources, or we generalize the sparse simulations, i.e., by training neural networks.

Our new methodology relies on the concept of (thermal) surrogate model, which is a neural network representation of a thermophysical scheme. The “parent” model is used in forward mode to populate a dataset of thermal simulations, which are then used to train, validate and test the networks (Section 2.2). The surrogate model assimilates the model and predicts answers at a known level of accuracy (with respect to the testing examples) and in an extremely low computational time. These properties enable detailed mapping of the parameter space, as it was not possible before. Because of the above, the surrogate models are ”statistically invertible”; they can be used to sample the (unknown) posterior distributions of the input parameters which is associated to an observed output. In our study, we adopted a Bayesian approach to inversion, using a Markov Chain Monte Carlo technique. A MCMC algorithm runs ten of thousands of simulations, inputting guesses of the solution associated with an observed set of data. The output is then compared to the observed data, and the residual is used to compute the probability that the guessed set of parameters is the optimal solution.

The developed methodology has been tested on two different cases. Firstly, we simulated a remote sensing problem for asteroid (101955) Bennu, the target of the NASA OSIRIS-REx mission, of which the infrared camera on-board the spacecraft (OTES, Christensen et al. 2018) will retrieve the emitted infrared fluxes. We used the TPM model (Delbo et al. 2015) to generate a grid of thermal simulation, and then we used it to train a set of neural networks. The peculiar geometry of the observations (which allows observing also the night side of the asteroid), the wide range of wavelengths that will be observed, as well as the expected high quality of the observations (characterized by the signal-to-noise-ratio) allow retrieving the surface properties of the asteroid - surface roughness, thermal inertia of the regolith and rock, and relative rock abundance - at a high accuracy, without suggesting any a-priori information to the Bayesian algorithm (Section 3.1). The accuracy of the method is expected to further increase if additional information about surface properties provided by other OSIRIS-REx instruments (e.g., OVIRS, OCAMS) are recalled.

The use of uninformative a-priori distributions, however, is not always possible. When data are partial (e.g., ground-based observations, which cannot image the night-side of the asteroid), and/or we lack information at long wavelengths, the use of well-educated guesses about surface properties of the asteroid must be employed. This necessity clearly emerged during our second exercise, in which we inverted observed infrared flux of asteroid (25143) Itokawa, the target of the Hayabusa mission (Section 3.2). Data from the spacecraft have been coupled with models to firstly inform the Bayesian inversion and therefore verify the consistency of the retrieved quantities. Following the proposed procedure (creation of the training dataset, training of the surrogate model, Bayesian inversion), we constrained the surface properties of the asteroid. In Section 3.2.2, we discuss an apparent bimodality in the posterior distribution of the surface roughness, which qualitatively matches with in-situ observations of the presence of both smoothed and rough terrains on the asteroid (e.g., Yano et al. 2006; Susorney et al. 2018, and reference therein.). In Section 3.2.3, we retrieve the thermal inertia of the regolith and the rock components, and relative rock abundance. The rock abundance on the surface (i.e., abundance of material whose size is larger than the diurnal skin depth) is found to be around 85%. The thermal inertia of the regolith suggests a mean particle size around 10 $mm$ in diameter, in agreement with previous findings (Yano et al. 2006; Kitazato et al. 2008). The value for the rock indicates a low thermal conductivity (if compared to LL chondrites), which could suggest the presence of fractured rocks. If the retrieved parameters are combined (Equation 11), we get an average value of thermal inertia equal to $785\pm 119$ $Js^{-1/2}K^{-1}m^{-2}$, consistent with previous findings by Müller et al. (2005), who suggested a value around 750 $Js^{-1/2}K^{-1}m^{-2}$ (Figure 8).

The present study is optimized towards the interpretation of infrared fluxes by spacecraft mission operating in proximity of airless bodies. In such scenario, the asteroid properties (e.g., albedo and size) and morphology are well characterized by in-situ measurements, and the infrared response of each facet of the shape model can be characterized by means of an associated surrogate model, leading to a global mapping of the surface properties. On the other hand, the interpretation of disk-integrated data could be affected by heterogeneity in rock abundance as well as presence of large boulders of the surface; the inversion of the surrogate models informs about the co-presence of regolith and boulders on the surface, while it does not resolve their relative distribution over the surface. Additionally, when shape models have non-negligible uncertainty in the size, the optimization of size or albedo as free parameters should be included. We will discuss how to overcome the application limit of our method to disk-resolved data in a future work.

We foresee the use of the proposed methodology to a variety of infrared remote sensing problem. Immediate applications are the analysis of (101955) Bennu by means of forthcoming data from the OSIRIS-REx spacecraft (this case has been already modeled in Section 3.1) and of (162173) Ryugu, the target of the JAXA’s Hayabusa 2 mission (Kuninaka and Hayabusa 2 Project 2013). Both spacecrafts will acquire infrared fluxes (by using an infrared spectrometer and an infrared camera, respectively) in close proximity of the asteroid, therefore providing infrared fluxes also of the night-side. Our method may be used to characterize the fine-scale composition, which is important for shedding lights on the collisional history, as well as the surface evolution, of the target asteroids. In case of application to spacecraft data collected during an extended period of time, the change in heliocentric distance can be taken into account while training the networks. During the detailed survey for (101955) Bennu, however, the heliocentric distance of the asteroid will vary between 1.0 and 1.25 au, resulting in a change of the surface temperature of at most 10 K. We consider this a very small effect probably obscured by other sources of errors. In any case, we expect no difference in the accuracy of the determination of the surface physical parameters as long as the thermo-physical parameters are temperature independent. We will study temperature depended thermal inertia in a future work. Finally, the finding that rocks on Itokawa could be fractured is also relevant to the interpretation of future infrared observations of (101195) Bennu and (162173) Ryugu. Given the small size of these targets and their high rock abundance, differences in thermal inertia could be representative of more or less fractured rocks, rather than indicating the presence of regolith material (i.e., pebbles with size smaller than the diurnal skin depth).

The concept of thermal surrogate model is not limited to asteroid studies, but it can be applied whenever a thermal model is available. Measurements of comet’s surface temperature (e.g., by ROSETTA/VIRTIS at comet 67P/C-G, Coradini et al. 1999) can also be inverted in their regolith and rock contributions, providing insights in the processes of evolution of cometary surfaces. Improved rock abundance maps on Mars and the Moon can be generated, thanks to the large availability of thermal data from previous and current missions. In future, infrared fluxes of Jupiter’s moon Europa and Ganymede will be provided by the NASA’s Europa Clipper and the ESA’s JUICE missions; a surrogate thermal model of the outer ice shell behaviour would allow the characterization of the surface properties, informing future landing mission about sampleability at the diurnal skin depth level.

Acknowledgements

MD and AR acknowledge support from the Academies of Excellence on Complex Systems and Space, Environment, Risk and Resilience of the Initiative d’EXcellence ”Joint, Excellent, and Dynamic Initiative” (IDEX JEDI) of the Université Côte d’Azur, as well as the Centre National d’Etudes Spatiales (CNES). The authors thank the anonymous referees for the precious comments and suggestions that improved this manuscript.

5 References