Unsupervised Learning for Thermophysical Analysis on the Lunar Surface

The Diviner instrument is a passive Vis-infrared (IR) radiometer on board NASA's LRO that has been collecting radiance measurements in nine different wavelength channels (0.3–200 μm) from its near-polar orbit since 2009 (Paige et al. 2010; Williams et al. 2017; Mazarico et al. 2018). Channels 6 through 9, located in the 12.5–200 μm range (middle to far-IR), are particularly sensitive to surface day and night temperatures. Each channel is defined by a linear array of 21 thermopile sensors (i.e., the instrument has a total of 189 sensors) and due to the limited field of view (FOV) of the sensors, Diviner's swath width is very narrow, resulting in quasi-point measurements along LRO's ground track (nadir push broom mapping). Thus, multiple orbits are required to completely cover the entire surface of the Moon and the number and effective FOV of the available measurements governs the spatial and temporal resolution of Diviner's final image or map products. The dominant influence on the size of the effective FOV of each point measurement is the spacecraft altitude, which has varied between ∼20 and 200 km over Diviner's operational period; at 30 km altitude, it is expected to be an ellipse approximately 160 × 120 m in size with its major axis oriented in the in-track direction, while at 200 km altitude, these dimensions are expected to double or triple (Sefton-Nash et al. 2017). Two other effects on the FOV include in-track smearing due to LRO's motion and footprint elongation caused by sensor memory (Paige et al. 2010; Williams et al. 2016). The in-track smearing can result in slight inaccuracies in the location of maximum probability of the source of back-scattered photons, but its significance weakens with increasing spacecraft altitude.

There are multiple Diviner data sets at different levels of processing publicly available from NASA's Planetary Data System. For this study, we collect the Diviner Level 1 Reduced Data Record (RDR) data (Paige et al. 2011), which consists of tables of calibrated radiance point measurements and their derived surface temperature values, as well as their associated ephemeris and geometry information, including the local lunar time at which each point measurement was recorded. We choose this data set because it is minimally processed from the raw Diviner telemetry data and contains the point measurements before any gridding is applied.

We download data from Diviner's entire operational period 2010 January–2019 March, which contains data recorded during LRO's nominal science mission phase as well as its three extended science mission phases. We use all available data for a number of reasons. First, because we are using a data-driven approach, we want to absorb as much data as possible and do not wish to bias our analysis to a narrow observational window. Second, as noted above, including more orbits increases the spatial-temporal resolution of the data. An important consideration is that LRO changed from a near-circular to an elliptical orbit in 2011 December, with its periapsis oriented over the lunar south pole. Thus, the effective FOV of point measurements is expected to vary with latitude after this date (Sefton-Nash et al. 2017).

A number of preprocessing steps are carried out before inserting this data into our unsupervised learning workflow. First, for each RDR data table, we only keep measurements that have instrument activity flags equal to 110 (indicating "on moon" orientation, standard nadir observation, and a "nominal" instrument mode), calibration flags equal to 0 (indicating an in-bounds interpolated measurement), geometry flags equal to 12 (indicating tracking data was used to generate geometry information), miscellaneous flags equal to 0 (indicating no miscellaneous observations), and emission angles less than 10° (the angle between the vector from the surface FOV center to Diviner and the "normal" vector to the Moon's surface). Second, for simplicity, in this study we only use data from channel 7, which records radiance in the 25–50 μm band, is the most sensitive to temperature values in the range 69–178 K, and has a high single-to-noise ratio over the majority of surface temperatures (Vasavada et al. 2012; Williams et al. 2017; Feng et al. 2020). Finally, we sort the data into 05 × 05 latitude and longitude bins, which allows temperature values at each point in space to be easily extracted in our subsequent workflows. This preprocessing step takes considerable computational resources; 400 CPU cores are used to process over 425 billion point measurements contained in 35 TB of uncompressed Diviner RDR data. Example preprocessed data are shown in Figure 1.

Zoom In Zoom Out Reset image size

Deep learning is a subfield of machine learning, which itself a subfield of AI, and offers a class of machine-learning algorithms that use neural networks with many layers to learn complex, nonlinear relationships in data. A quintessential example is the deep feed-forward network (Goodfellow et al. 2016). Given a set of training samples ${ \mathcal D }=\{({x}_{1},{y}_{1}),...,({x}_{n},{y}_{n})\}$ from a function $y={f}^{* }(x):X\to Y,$ where X and Y are vector spaces, the deep feed-forward network defines a mapping y = f(x; θ) and attempts to learn the values of the parameters θ that best approximate f^*. The mapping is typically composed of multiple functions, or layers, of the form $f(x)={f}^{(d)}\circ ...\,\circ {f}^{(3)}\circ {f}^{(2)}\circ {f}^{(1)}(x),$ where the output dimensionality (or "width") of each function and the number of functions (or "depth") can vary. One standard form is the fully connected network, where each layer function has the form ${f}^{(i)}(x;{W}^{(i)},{b}^{(i)})=g({W}^{(i)\top }x\,+\,{b}^{(i)}),$ where the weight matrix W⁽ⁱ⁾ and bias vector b⁽ⁱ⁾ contain learnable parameters and g is an activation function, often chosen to be the rectified linear unit, given by $g(x)=\max \{0,x\}$. Rather remarkably, the universal approximation theorem (Hornik et al. 1989) states that a fully connected feed-forward network with 2 or more layers can approximate any continuous function for inputs within a specific range. It is this flexibility that potentially allows deep neural networks to learn across many different problems.

CNNs are a specific type of deep neural network that have become hugely popular in recent years. In a CNN, each component in a layer's output vector (or "neuron") shares the same weight values as the other output components. Furthermore, each neuron has a limited FOV of the layer's input vector. Mathematically, the layer function is defined as ${f}^{(i)}(x;{w}^{(i)},{b}^{(i)})=g({w}^{(i)}\star x+{b}^{(i)}),$ where ⋆ denotes cross-correlation and w⁽ⁱ⁾ is a vector of a fixed size, or filter, which does not depend on the size of the layer's input, and b⁽ⁱ⁾ is a scalar. CNNs are nearly always extended to allow multiple filters (or "channels") per layer and 2D and 3D cross-correlation over its inputs. An important feature of CNNs is that they are equivariant to translations of their inputs, which is extremely valuable when learning about data that contains spatial correlations, such as images. CNNs have been successfully applied in computer vision, where they achieve state-of-the-art performance (Russakovsky et al. 2015), as well as many other areas.

VAEs use deep neural networks to carry out unsupervised learning (Kingma & Welling 2014; Rezende et al. 2014; Doersch 2016). Rather than estimating a function, they use deep neural networks to carry out approximate inference in a latent probabilistic model. The goal is to (implicitly) learn a distribution $p(x),x\in X$, where X is a vector space, which is as similar as possible to an unknown target distribution p*(x), given a set of samples from this distribution ${ \mathcal D }=\{{x}_{1},...,{x}_{n}\}$. VAEs make the assumption that x depends on some set of underlying latent variables $z\in Z$, which are distributed according to some known distribution p(z), and use deep neural networks to model the likelihood $p(x| z)$ and an approximation of the posterior distribution $p(z| x)$. The latent representation is learned by training the VAE to reconstruct the training data. One typically chooses the dimensionality of the latent representation to be much lower than the dimensionality of X, forcing the VAE to learn a salient representation of the input data. Furthermore, in a standard VAE, these latent variables are assumed to be independent and normally distributed, which has the effect of keeping input samples that are similar in X close together in the latent space Z and is very useful when one wants to learn "meaningful" representations that vary smoothly over the input space. In computer science, VAEs have successfully been applied to disentangle hair styles, poses, and facial expressions from images of people's faces (Larsen et al. 2015; Higgins et al. 2017). By using deep neural networks to infer $p(x| z)$ and approximate $p(z| x)$, VAEs are able to learn highly nonlinear mappings between the input and latent space. A detailed mathematical description of VAEs is included in Appendix A.

Our basic idea is to use a VAE to learn the underlying factors of variation in the Diviner surface temperature measurements. The input to the VAE is the observed surface temperature variations over the lunar day at each point on the surface, and we learn a latent representation that is able to reconstruct this input. In doing so, we learn a model of the thermophysical data, without the use of an a priori physics model. We choose a VAE framework (as specified in Equation (A3)–(A6)) because it attempts to learn a model that (1) accurately reconstructs the input profiles, (2) disentangles independent factors of variation in the surface temperature, and (3) has a latent representation, which varies smoothly over the input space. We stress that these are very general assertions and therefore this method could be applied to other similar space science and exploration tasks.

A major benefit of this approach is that it is fully data driven. In general, methods that rely on physical models to analyze physical data may be unreliable if one is not confident on the accuracy of the model (e.g., Hayne et al. 2017). Furthermore, we learn salient features directly from the data; other traditional approaches that uses handcrafted features (such as averaged or evening temperature maps; e.g., Bandfield et al. 2011) to interpret physical data risk discarding important information. We interpret the VAE by comparing its latent variables to parameters estimated using inversion with an existing thermophysical model. Once trained, the VAE is run over the entire lunar surface to map its latent variables and these maps are used to search for thermal anomalies. We name the VAE "Lunar ORbiter-Derived VAE," or LordVAEder.

Our VAE workflow is shown in more detail in Figure 2. The input to the VAE is the 1D profile of the surface temperature variation over the lunar day extracted at a single point on the lunar surface. The VAE has multiple "encoder" convolutional layers that reduce this profile into its latent representation and multiple "decoder" convolutional layers that reconstruct it, and it is trained using example profiles from many different locations on the lunar surface.

Zoom In Zoom Out Reset image size

To train the VAE, millions of example 1D profiles of the surface temperature variation over the lunar day are extracted from 48 areas of interest (AOI). We curate the AOIs so that they encapsulate a wide range of surface temperature variation. The AOIs include impact craters of different ages, large pyroclastic deposits, and swirls (magnetic anomalies), with areas ranging from 1 to 10,000 km². They are listed in Table B1 and described in more detail in Appendix B. The locations of all the AOIs are shown in Figure B1.

We extract the 1D temperature profiles by locally projecting and binning the observed Diviner surface temperature measurements at each AOI onto a 200 × 200 m grid and sorting the points in each bin by the local lunar time at which they were recorded, obtaining a profile at each bin location. An orthographic projection centered over each AOI with a radius of R_moon = 1737.4 × 10³ m is used to define each grid. To ensure sufficient sampling of each temperature profile, we reject profiles whose maximum spacing in local lunar time between points is less than 4 hr. We also reject some profiles from background areas to ensure they do not dominate the training set. This workflow results in 1,985,693 training profiles extracted over all AOIs and examples of extracted profiles are shown in Figure 1.

Before inputting the profiles into our VAE workflow, we interpolate their temperature measurements onto a regular grid by using Gaussian process (GP) regression (Rasmussen & Williams 2005). A GP with a Matern-1.5 kernel, a maximum length scale of 6 hr, and an assumption of 10 K standard deviation of uncertainty in each temperature measurement is fit to each profile. We sample every 0.2 hr, resulting in a 1D input profile with 120 samples. Example GP interpolations are also shown for the profiles in Figure 1.

A convolutional encoder-decoder design is used for our VAE architecture, shown in Figure 2, which consists of nine convolutional layers in the encoder and eight convolutional layers in the decoder. The encoder takes an input profile of dimensionality 120 × 1 and slowly reduces it to two output vectors each of dimensionality 1 × n, which represent the mean, g, and standard deviation, h, of the approximate latent posterior distribution in Equation (A5). The number of latent variables, n, is a hyperparameter of the VAE model we test. The decoder takes a sample from the latent posterior distribution and slowly expands it to output a reconstruction of the input profile.

We use 32 hidden channels in the encoder and decoder and batch normalization (Ioffe & Szegedy 2015) after all hidden layers to help training converge. For n = 4, the total number of trainable parameters in the model is 28,713. A complete specification of the network is given in Table B2 and the network architecture is described in more detail in Appendix B. We design the network architecture by hand to maximize its reconstruction performance.

During training, the mean and standard deviation vectors are used to generate a single sample (or vector of latent values) from the approximate posterior distribution, which is fed into the decoder so that the expectation term in the loss function (Equation (A6)) can be estimated. During inference, only the mean vector is passed to the decoder and used to define the latent values. We use a β-value of 0.2 and normalize the training profiles to have a collective mean of 0 and a standard deviation of 1 before inputting them to the VAE. The network is trained using the Adam stochastic gradient descent optimization algorithm (Kingma & Ba 2014), using an exponentially decaying learning rate starting at $1\times {10}^{-3}$ and a mini-batch size of 200 profiles, reserving 20% of the example data set for cross-validation. We use the PyTorch framework (Paszke et al. 2018) to define and train the VAE.

A downside of our approach is that the VAE model lacks a clear physical interpretation. To aid its interpretation, we compare the VAE's latent representation to a physics-based inversion. We modify the thermophysical model proposed by Hayne et al. (2017), which consists of numerically solving the heat equation over time t and depth z in a 1D solid medium:

$$\begin{eqnarray}&&\rho (z,t){C}_{p}(z,t)\displaystyle \frac{\partial T(z,t)}{\partial t}=\displaystyle \frac{\partial }{\partial z}\left(K(z,t)\ \displaystyle \frac{\partial T(z,t)}{\partial z}\right),\end{eqnarray} \tag{ 1 }$$

where ρ(z, t) is the mass density, C_p(z, t) is the specific heat capacity, T(z, t) is the temperature, and K(z, t) is the thermal conductivity. The boundary conditions and functions ρ(z, t), C_p(z, t), and K(z, t) must be known in order to solve this differential equation for the temperature profile. We make a few important alterations to the model used by Hayne et al. (2017). First, we assume that the density ρ is constant and equal to the surface density ρ_s given by Hayne et al. (2017). Second, we use the phonon conductivity K_c as a free fitting parameter, rescaled with respect to the surface conductivity K_s given by Hayne et al. (2017). For this purpose, we introduce the thermal conductivity factor as K_c/K_s. Third, we assume that the solar incidence angle approximately coincides with the hour angle and absorb the effect of the latitude and slope on solar incidence in an effective Bond albedo, which we use as another free fitting parameter. Hence, our model treats the lunar regolith as a solid medium with a constant density but with a varying effective albedo and thermal conductivity factor.

We run inversion using this model by tuning the thermal conductivity factor, albedo, and onset time of solar radiation to best match each input interpolated profile. Bayesian optimization with an expected improvement acquisition function and a L2 loss function is used as the inversion algorithm (Frazier 2018). We then compare the estimated parameters from the inversion to the latent variables generated by the VAE.

As a final step, we run the trained VAE with n = 4 over the entire lunar surface and generate maps of its latent variables to search for thermal anomalies. This allows us to compute global statistics of the learned latent variables and interpret these new features alongside existing physical maps.

The first step in our mapping workflow is to run VAE inference and to obtain latent values over the Moon's entire sphere. To do this, we tessellate the sphere into many latitude–longitude boxes of approximately equal areas. For each box, we bin the observed Diviner surface temperature measurements onto a 200 × 200 m grid, using an orthographic projection centered at the center of each box, obtaining a profile at each bin location. Each profile is interpolated using GP interpolation, rejecting profiles whose maximum spacing in local lunar time between points is less than 4 hr. The VAE is run on each interpolated profile, obtaining a set of latent values at each bin location, and the central coordinate of each bin is unprojected to obtain a set of latent value latitude–longitude points for each box.

The second step is to project and interpolate the latent value latitude–longitude points onto a global map. We use the standard plate-carrée map projection (equirectangular) and linearly interpolate the points onto a 100 × 100 m grid. Similar to the first step, we do this by tessellating the sphere into many latitude–longitude boxes of approximately equal area. For each box, we project the latent value latitude–longitude points onto the global map and use linear interpolation with Delaunay triangulation to interpolate the points onto the global grid. Interpolated points are set to null values where their Delaunay triangle has a maximum distance between its vertices on the surface of the Moon's sphere greater than 2000 m to ensure that the map does not include over-interpolated points.

It is important to note that geometrical error in our workflow is introduced in the first step when binning temperature measurements onto the 200 × 200 m grid, because the bins in the orthographic projection are not exactly equally sized when unprojected onto the surface of the sphere. However, we choose small enough box sizes when tessellating the sphere so that this effect is negligible; 379,448 boxes of approximately 10 × 10 km shape on the Moon's surface are used for the first step and 10,322 boxes of approximately 61 × 61 km shape are used for the second step. For both steps, we ensure the borders of the boxes are sufficiently overlapping during the processing of each box to ensure that no edge effects are present. Part of the tessellation used for the second step is plotted in Figure C1. Both steps require considerable computational resource, but tessellation allows the data to be processed in parallel using 100 CPU cores. Our final output products are four latent maps and one quality control map that shows the L1 reconstruction loss of the VAE. The maps are cropped between −70° and 70° latitude to avoid excessive deformation and areas of high VAE reconstruction loss observed at the poles (see Appendix C and Figure C2 for more details).

In total 425,745,195,120 Diviner point measurements are extracted, of which 326,190,541,311 (77%) remain after the first preprocessing filtering step, and 36,483,372,924 remain after filtering channel 7, requiring 2.9 TB of storage. We plot the number of points after preprocessing in each latitude–longitude bin in Figure D1. There is a higher point density closer to the poles, which is expected because LRO's orbit means the Diviner instrument samples the poles more frequently than the equator.

Example data after preprocessing over the Tycho crater are plotted in Figure 1. We observe a wide variation in temperature around this crater. Outside of the crater, the surface temperature increases to around 350 K at midday and decreases to 100 K at night. On the eastern slope of the crater the profiles are shifted in local time by approximately 2 lunar hours, which is consistent with the slope's shadow delaying the incidence time of the Sun's radiation. On the northern slope the peak temperature is lower at 270 K, which is consistent with the slope receiving less thermal radiation per unit area from the Sun's elevation at Tycho's latitude. In the center of the crater we observe a profile which has a higher evening temperature around 150 K. Figure 3 (right) plots our physics model simulation when varying thermal conductivity factor; one explanation is that the material at Tycho's center has a higher thermal conductivity than its surroundings, but it should be noted that its central peak has complex topography and surface roughness that are likely to influence the measured temperature profiles.

Zoom In Zoom Out Reset image size

While training the VAE, the training and validation reconstruction loss values converge to similar values, suggesting that the VAE generalizes well and does not over-fit to the training data set. Training takes approximately 1 hr using one Nvidia Tesla T4 GPU. Example reconstructed profiles from the validation data set after training are plotted in Figure 4 (top middle panel). We find that the VAE is able to accurately reconstruct the profiles from the validation data set. Figure 4 (top left panel) shows the average reconstruction loss over the validation data set after training the VAE when varying the number of latent variables. The reconstruction loss asymptotes at n = 4, suggesting that only four latent variables are needed to accurately represent the input profiles.

Zoom In Zoom Out Reset image size

Figure 5 shows example profiles generated by the VAE when varying each latent variable independently and fixing the others to zero for n = 4. We observe that latent variable 1 responds to the onset of the profile, latent 2 responds to the peak temperature, latent 3 responds to the evening temperature "sidelobe," and latent 4 responds to the width of the temperature profile. To aid our interpretation of the VAE further, we compare its latent variables with estimated parameters from the physics inversion. The VAE and physics inversion are run over the validation data set and the latent values are plotted against the estimated physics parameters in Figure 4 (bottom). We find that latent variables 1, 2, and 3 have a strong linear correlation with the onset of solar radiation, effective albedo, and the logarithm of the thermal conductivity factor, respectively. This correlation can also be seen in Figure 3; scanning the physics simulation across the physical parameters results in similar profile variations to those seen when scanning the latent variables in Figure 5.

Zoom In Zoom Out Reset image size

A straight line fit is shown for each correlation plot in Figure 4 (bottom panel). As thermal inertia is proportional to the square root of thermal conductivity, we expect that applying the transform $\hat{I}={e}^{{{mz}}_{3}/2}$ to the latent 3 values will obtain a quantity that is linearly correlated with the thermal inertia at constant density and temperature, under the assumptions of our physics model. Here m is the gradient of the latent 3 and log thermal conductivity factor correlation and is estimated to be 0.93.

Maps of the latent variables and the VAE L1 reconstruction loss over Tycho are plotted in Figure 6. The VAE takes approximately 0.001 s to compute latent values for each profile when running on a single CPU, while each physics inversion takes of the order of 60 s to run using the same CPU. This significant efficiency improvement (factor of 60,000 or 4 orders of magnitude) allows the VAE to be run over large areas.

Zoom In Zoom Out Reset image size

The global maps produced by our mapping workflow are shown in Figures 7 and 8. We find that our observations of the latent variables are consistent on a global scale; the latent 1 map consistently shows higher values on the eastern sides of craters and lower values on their western sides, indicating that latent 1 responds to the onset delay of solar radiation as a function of slope (east–west) aspect. The latent 2 map shows higher values at the equator and lower values closer to the poles, suggesting that it represents a combination of an effective incidence angle and surface albedo. The latent 3 map mostly shows a background value with some craters and recent impact basins being either anomalously bright (e.g., the Tycho crater) or showing a brightness level slightly above the background (e.g., Mare Orientale) and could plausibly indicate changes in thermal conductivity. The latent 4 map is generally darker within smaller and narrower craters, particularly toward the poles. Besides their dominant trends, both the latent 1 and 2 maps show minor, local variations in the equatorial regions (for example, around the main maria (centered at 315° lon./ 10° lat. or 45° lon./ 0° lat.)) that could represent the influence of terrain roughness on latent 1 and varying albedo on latent 2. We define an anomalously high value of the latent 3 variable to be one that is greater than two standard deviations from its global mean and find that approximately 4% of the surface area of the Moon (within −70° to 70° latitude) contains these values.

Zoom In Zoom Out Reset image size

The global map of the VAE L1 reconstruction loss is shown in Figure 9. The loss is low (typically less than 10 K) in the range of −70° to 70° latitude, while it increases dramatically closer to the poles. In general, the loss slightly increases along topographic gradients such as the rims of craters and in the more heavily cratered highland regions. Locally, the loss increases at east–west-facing slopes of some impact features such as the Chaucer or Rydberg craters. A regional exemption is the north-western and central portion of Mare Smythii (85° lon./ 0° lat.), which show higher than average loss values. Some observations in the loss map are discussed in more detail in Appendix C.

Approximately 80% of the surface area of the Moon is covered by our extraction workflow. The missing data are where the extracted profiles do not meet our minimum temporal sampling criteria and are mostly located around the equator where the density of point measurements is lowest due to LRO's near-polar orbit.