Unsupervised Spectral Unmixing For Telluric
Correction Using A Neural Network Autoencoder

Rune D. Kjærsgaard
DTU Computerdokj@dtu.dk

Aaron Bello-ArufeDTU Space
aarb@space.dtu.dk

Lars A. Buchhave
DTU Computebuchhave@space.dtu.dk

Alexander D. Rathcke
DTU Spacerathcke@space.dtu.dk
Line K. H. Clemmensen
DTU Computelkhc@dtu.dk

Abstract
The absorption of light by molecules in the atmosphere of Earth is a complication
for ground-based observations of astrophysical objects. Comprehensive information on various molecular species is required to correct for this so called telluric
absorption. We present a neural network autoencoder approach for extracting atelluric transmission spectrum from a large set of high-precision observed solar
spectra from the HARPS-N radial velocity spectrograph. We accomplish this by
reducing the data into a compressed representation, which allows us to unveil
the underlying solar spectrum and simultaneously uncover the different modes ofvariation in the observed spectra relating to the absorption of H2 O and O2 in theatmosphere of Earth. We demonstrate how the extracted components can be used
to remove H2 O and O2 tellurics in a validation observation with similar accuracy
and at less computational expense than a synthetic approach with molecfit.

1

Introduction

Absorption of light in the atmosphere of Earth, called telluric absorption, can hinder astrophysicalobservations by partially obscuring the object of interest. Various methods have been introduced
to remove the effects of this absorption from observed spectra. One such acknowledged method,
called molecfit [15, 12], relies on computing a synthetic transmission spectrum of the atmosphereof Earth. Synthetic methods are inherently reliant on external factors to an observation, such asatmospheric measurements and molecular line lists. Another realm of methods take a data-driven
approach attempting to exploit the modes of variation in a number of observed spectra to uncover theunderlying components. By analysing such variation, the telluric absorption can be modelled without
relying on external factors to an observation. One such approach, based on principal component
analysis (PCA), has been explored in the literature [2]. PCA methods are however ineffective on verylarge data sets, where the entire data can not be stored in memory. Another approach called wobble[3] uses a linear model with a convex objective to model the telluric component of observed spectra.
We present a new data-driven approach using a neural network autoencoder. Autoencoders have seenuse in the literature for decades [4, 10] and have long been known to discover effective compresseddata representations through dimensionality reduction [9]. To demonstrate the approach we analyse
1257 observed solar spectra [7, 5] from the high-precision spectrograph HARPS-N [6] with the aimof disentangling the observed spectra into an underlying solar component and high accuracy telluric
components from H2 O and O2 . The extracted components could aid in the detection of radial velocity
signals of planetary systems by quickly and accurately removing tellurics from observed spectra,
leading to an increase in observation quality and hereby a reduction in observing time and cost.
Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

2

Proposed approach

Our approach has roots in spectral unmixing, which seeks to unmix distinct endmember spectra
and their weights from an observed spectral image by constructing a mixing model of the problem.
Endmember unmixing from spectral data is a rich discipline with many existing approaches [16, 8].
We consider the solar spectrum and telluric spectrum as endmember spectra with associated abundance
weights and use these components to construct a linear mixing model in log-space [17] describing
the observed spectra:

xn =

R
X

wr,n mr = M wn + ϵn ,

(1)

r=1

where xn is the nth observed spectrum from a finite set of N observed spectra, mr is the rth
endmember spectrum of R endmembers with individual endmembers r = 1, ..., R. Furthermore,
wr,n is the abundance of endmember r for observation n, M is the endmember matrix having
endmembers as columns, wn is the abundance vector of the nth observation and ϵn is an error
term. The HARPS-N spectra cover the optical wavelength range and for this reason we consider the
combined telluric spectrum to be comprised of the two strongest absorbing molecules in this region,
namely H2 O and O2 . From this we get R = 3 with r = 1 representing the solar endmember, r = 2
representing the H2 O endmember and r = 3 representing the O2 endmember.
The goal is to extract the endmember matrix M and abundance vector wn by training the neuralnetwork on a set of preprocessed observed solar spectra, with the purpose of applying the extracted
telluric spectrum to non-solar observations. The endmember matrix M is extracted for training
regions with P pixels. Figure 1 shows a graphical representation of the approach.
In [14] they present an autoencoder for blind unmixing of hyperspectral images (HSI). We buildon this idea by adapting the network architecture to the domain of astrophysical spectral data. Thisrequires various structural changes and the introduction of several new constraints on the network.
The training data consists of 1257 observations of the solar spectrum all split into 69 apertures with
P = 4096 pixels each. This gives a total of 69 × 4096 pixels per observation. We use solar data fortraining since this data has a very high signal to noise ratio and does not take away observing timefrom night time observations. The observations are from 2020 spanning a month of observations from
October 22 through November 19. We filter out low flux and high airmass observations. Subsequently
we interpolate all observations to a common wavelength grid, apply the natural logarithm andcontinuum normalise the spectra. The described procedures leave N = 838 spectra of which 50%
are reserved for training and 50% are used for validation.
Endmembers

Endmember Abundances

Reconstruction

Input

Encoder

LatentDecoder

Figure 1: The figure shows the architecture of the autoencoder. Observed spectra xn are given as
input and passed through the encoder into a lower dimensional latent space, which is subsequently
decoded into the reconstruction x̂n . After training by minimising the reconstruction error througha gradient descent algorithm, the endmember matrix M is extracted as the weights of the decoder
and the abundance vector wn is extracted as the latent representation hn . P is the number of pixels
for each aperture in the observed spectrum. The network is illustrated for R = 3 endmembers
representing the solar (orange, top), H2 O (blue, middle) and O2 (green, bottom) endmembers.

2

2.1

Neural network autoencoder

The network uses fully connected layers and consists of an encoder function hn = f (xn ), whichmaps the input data xn ∈ RP to an internal latent representation hn ∈ RR . This representation isthen passed through the decoder function g(hn ) = x̂n , which seeks to reconstruct the input data xn
with the reconstruction x̂n ∈ RP . The decoder is constructed without bias terms and performs thefollowing affine transformation:x̂n = W hn ,(2)where W ∈ RP ×R are the weights of the decoder, which are extracted as the endmember spectra M ,and hn ∈ RR is the latent representation, which can be extracted as the endmember abundances wn .
The network features a utility layer responsible for normalisation of endmember abundance variation,
which ensures a fixed abundance of the solar component, and a utility layer clamping the solar decoder
weights to interval [0, 1] and the telluric decoder weights to interval [−1, 0]. Moreover, the network
contains a utility layer Doppler shifting the solar decoder weights to account for the spectral shift
caused by the rotation and elliptical orbit of the Earth. We perform this shift using the barycentric
Earth radial velocity of each observation. The network also features a batch normalisation layer[13, 11]. We use the validation set to determine hyperparameters based on a Tree-structured Parzen
Estimator approach carried out with optuna [1]. Training is performed with stochastic gradientdescent to minimise the mean squared error (MSE). The network is trained on non-stitched spectra
for the 69 apertures separately to retain the high fidelity of the observed spectra and to avoid thecomplications involved in stitching spectra. Training on all apertures takes approximately 3 hours
and 15 minutes on an Intel 6 core i7, UHD 630 CPU laptop.

3

Results

The extracted endmembers for aperture 60 are shown in Figure 2.

Aperture 60 Extracted Endmembers

Normalized Log Flux

0.0
0.1
0.0
0.1
1.0
0.9
1.0
0.9

H2O
6270

6280

6290

6300

6310

6320

6330

6270

6280

6290

6300

6310

6320

6330

6270

6280

6290

6300

6310

6320

6330

6340

O2
6340

Solar
6340

Observed

6270

6280

6290

6300
6310
Wavelength [Å]

6320

6330

6340

Figure 2: The figure illustrates extracted endmembers in addition to an observed solar spectrum foraperture 60. The extracted endmembers represent from top to bottom the H2 O (blue), O2 (green) andsolar (orange) components of the observed spectrum (black, bottom).
To validate extracted endmembers we perform telluric correction using our autoencoder tellurics andcompare with a molecfit synthetic telluric spectrum. Telluric correction aims to remove tellurics bydividing the observed spectrum with a telluric transmission spectrum. We perform the comparisonby correcting a solar observation with strong tellurics from HARPS-N, which the autoencoder has
not been trained on. We compute the molecfit telluric spectrum on a HPC cluster using a 10 coreIntel Xeon E5-2660v3, Huawei XH620 V3 node and utilise atmospheric measurements from the timeof the observation in addition to a fit to the stitched version of the observation from the HARPS-N
pipeline. We compute the autoencoder correction using the extracted telluric components, which havebeen converted back from log-space to represent standard transmission spectra. Autoencoder telluric
abundance weights are found using a least squares fit to known telluric lines in the spectrum. We
interpolate the observed spectrum to the telluric wavelength axes of molecfit and the autoencoder
before the corrections. The comparison can be seen in Figure 3.
3

Flux
Flux

12500010000075000
5945.0150000Autoencoder vs. Molecfit Telluric Correction

5945.55946.05946.5

1000006290.06290.56291.06291.56292.0150000Flux

5947.06293.0

Molecfit Corrected Spectrum
Autoencoder Corrected Spectrum
Observed Spectrum

1000006542.06292.56542.56543.06543.56544.0

Wavelength [Å]6544.56545.0Figure 3: The figure displays a comparison of corrected spectra computed either using autoencoder
extracted tellurics or molecfit tellurics for H2 O and O2 lines in three spectral regions of highinterest to the study of exoplanetary atmospheres. For pure telluric lines, an ideal correction would
result in a flat line. In the middle plot, the autoencoder and molecfit agree on the presence of solarfeatures around 6290.5 Å, 6291 Å and 6292.75 Å, which remain in the corrected spectra.

4

Discussion

High-resolution spectroscopy is limited to ground-based observations, and consequently highresolution ground truth solar spectra are not available. This makes evaluating the performancedifference between the two approaches complicated. A natural area of comparison however lies
around corrected tellurics, where an accurate correction will leave no trace of the telluric absorption.
As illustrated in Figure 3, the autoencoder correction removes tellurics in the observed spectrum tocontinuum level, while molecfit leaves slight traces of the correction. The sinusoidal shape of themolecfit correction around 6544 Å in Figure 3 could indicate imprecision on the exact wavelength
location of the telluric line centre used in molecfit. The molecfit correction takes approximately
30 minutes to compute, while the autoencoder correction takes about 0.2 seconds. This differencein compute time is significant and makes the autoencoder much more feasible for correction ofmultiple spectra. Additionally, the autoencoder approach can be used as a complementary data-driven
validation tool to inspect the accuracy of synthetic approaches like molecfit.
While the results show that the autoencoder correction of the validation observation is performed
with similar accuracy to molecfit, the question of how well the extracted endmembers generalise toobservations dissimilar to the training data still remains. Future work includes exploring this potential
limitation by performing corrections on numerous non-solar observations and inspecting the impact
on radial velocity extraction, as well as the consistency of retrieved exoplanetary atmosphere signals.
This paper demonstrates the approach applied to the HARPS-N spectrograph, but the telluric autoencoder is designed as a general tool, which can be trained on solar data from any spectrograph andwavelength range and subsequently perform corrections on new observations (including non-solar
observations) from the given spectrograph. This is an advantage as solar data is already gathered andready to train on for many spectrographs.

5

Conclusion

We have demonstrated an approach for computing a compressed representation of the solar data fromHARPS-N with a constrained neural network autoencoder. This representation can be used to extract
endmembers that directly relate to the solar spectrum as well as the transmission spectra of H2 O andO2 . After the autoencoder representation has been computed for a given spectrograph, the extracted
components can be used to perform very quick and accurate telluric correction on any observations
from the same spectrograph. The autoencoder approach and the detailed extracted telluric spectrum
could aid in the detection of faint radial velocity signals and atmospheric features of Earth analogue
exoplanets observed from ground-based telescopes.
4References[1] Akiba, T., Sano, S., Yanase, T., Ohta, T. and Koyama, M. [2019], ‘Optuna: A next-generation
hyperparameter optimization framework [arxiv]’, Arxiv p. 10 pp.
[2] Artigau, É., Astudillo-Defru, N., Delfosse, X., Bouchy, F., Bonfils, X., Lovis, C., Pepe, F.,
Moutou, C., Donati, J.-F., Doyon, R. et al. [2014], Telluric-line subtraction in high-accuracy
velocimetry: a pca-based approach, in ‘Observatory Operations: Strategies, Processes, andSystems V’, Vol. 9149, International Society for Optics and Photonics, p. 914905.
[3] Bedell, M., Hogg, D. W., Foreman-Mackey, D., Montet, B. T. and Luger, R. [2019], ‘wobble:
a data-driven analysis technique for time-series stellar spectra’, The Astronomical Journal158(4), 164.
[4] Bourlard, H. and Kamp, Y. [1988], ‘Auto-association by multilayer perceptrons and singularvalue decomposition’, Biological cybernetics 59(4), 291–294.
[5] Dumusque, X., Cretignier, M., Sosnowska, D., Buchschacher, N., Lovis, C., Phillips, D., Pepe,F., Alesina, F., Buchhave, L., Burnier, J. et al. [2021], ‘Three years of harps-n high-resolution
spectroscopy and precise radial velocity data for the sun’, Astronomy & Astrophysics 648, A103.
[6] HARPS-N Instrument Page [n.d.], http://www.tng.iac.es/instruments/harps/.
[7] HARPS-N Published Solar Observations [n.d.], https://dace.unige.ch/dashboard/.
[8] Heylen, R., Parente, M. and Gader, P. [2014], ‘A review of nonlinear hyperspectral unmixing
methods’, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing
7(6), 1844–1868.
[9] Hinton, G. E. and Salakhutdinov, R. R. [2006], ‘Reducing the dimensionality of data with neuralnetworks’, science 313(5786), 504–507.
[10] Hinton, G. E. and Zemel, R. S. [1994], ‘Autoencoders, minimum description length, andhelmholtz free energy’, Advances in neural information processing systems 6, 3–10.
[11] Ioffe, S. and Szegedy, C. [2015], Batch normalization: Accelerating deep network training by
reducing internal covariate shift, in ‘International conference on machine learning’, PMLR,
pp. 448–456.
[12] Kausch, W., Noll, S., Smette, A., Kimeswenger, S., Barden, M., Szyszka, C., Jones, A., Sana,
H., Horst, H. and Kerber, F. [2015], ‘Molecfit: A general tool for telluric absorption correction-ii.
quantitative evaluation on eso-vlt/x-shooterspectra’, Astronomy & Astrophysics 576, A78.
[13] LeCun, Y. A., Bottou, L., Orr, G. B. and Müller, K.-R. [2012], Efficient backprop, in ‘Neuralnetworks: Tricks of the trade’, Springer, pp. 9–48.
[14] Palsson, B., Sigurdsson, J., Sveinsson, J. R. and Ulfarsson, M. O. [2018], ‘Hyperspectral
unmixing using a neural network autoencoder’, IEEE Access 6, 25646–25656.
[15] Smette, A., Sana, H., Noll, S., Horst, H., Kausch, W., Kimeswenger, S., Barden, M., Szyszka,C., Jones, A., Gallenne, A. et al. [2015], ‘Molecfit: A general tool for telluric absorption
correction-i. method and application to eso instruments’, Astronomy & Astrophysics 576, A77.
[16] Somers, B., Asner, G. P., Tits, L. and Coppin, P. [2011], ‘Endmember variability in spectral
mixture analysis: A review’, Remote Sensing of Environment 115(7), 1603–1616.
[17] Zhao, H. and Zhao, X. [2019], ‘Nonlinear unmixing of minerals based on the log and continuum
removal model’, European Journal of Remote Sensing 52(1), 277–293.

5Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [Yes] See Section 4.
(c) Did you discuss any potential negative societal impacts of your work? [No] Our
approach incurs low computational costs when performing numerous corrections. This
results in a net societal benefit compared to the continued use of more computationallyexpensive methods like molecfit.
(d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A](b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [No] The project is awork in progress. The code and data will be published when the project is complete.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [Yes] See Section 2.
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [N/A](d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes] See Section 2.1, Section 3 andSection 4.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] See Section 1.
(b) Did you mention the license of the assets? [No] The HARPS-N solar data is public. We
reference both the paper detailing the data release [5] as well as the site to download
the data [7].
(c) Did you include any new assets either in the supplemental material or as a URL? [No]The project is a work in progress. The code and data will be published when the projectis complete.
(d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A] The HARPS-N solar data is public and available for everyone.
(e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A] The astrophysical data does not containpersonally identifiable information or offensive content.
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]

6