Taxonomic analysis of asteroids with artificial neural networks

We constructed an ANN of an input layer, one hidden layer, and an output layer (see Fig.2). Neurons in the input layer take the elements of an input data or features of spectra, represented by a vector $X=[x_{i}],\;i=1,\ldots,R$. The number of neurons in the input layer $R$ is the number of features in the input spectral data. The $K$ neurons in the output layer present a discrete probability distribution $Y=[y_{j}],\;j=1,\ldots,K$ for $K$ classes ($j$ is the class index, and neuron $y_{j}$ in the output layer presents the probability of the input spectral data belonging to class $j$).

Here, we temporarily set the number of neurons in the input layer as $R=N_{in}=151$ considering the wavelength range of available training data (i.e., SMASS II, 0.435 $\sim$ 0.925 µm ), and a lower limitation of the resolution of 200 of the CSST slitless spectroscopic units. The concrete element $\lambda_{i}$ of input data of our ANNs is computed by Eq.(1), in which $\lambda_{0}$ is the shortest wavelength in the training data or test data, $(1+\frac{1}{R})^{i}$ ($i$, the index of element of input data) is the interval distance of two adjacent elements of $i$ and $i-1$. For a wavelength range of spectral data from 0.435 to 0.925 µm, the maximum of $i$ is 151:

The number of neurons in the output layer is set as $K=10$, namely, only ten most popular classes of asteroids are included in the presented ANNs. Detailed information about the chosen asteroid classes is described in Sec. 3.

The number of neurons in the hidden layer is usually hard to determine accurately due to its flexible property. Often, it is determined by iterating over different choices. Generally, more neurons in the hidden layer allow more complicated problems to be solved by the ANN. One should take care when choosing the number of neurons in the hidden layer, because too few neurons in the hidden layer would result in under-fitting the input data, while too large number of neurons can lead to over-fitting. In other words, for the under-fitting case, the ANN will give a low accuracy of prediction to the input data, and for the over-fitting case, the ANN could give a high accuracy to training data, but a low accuracy to the test data. In our case, we found $N_{h}=45$ neurons in the hidden layer to work well.

As a feed-forward neural network, the input information is transferred to neurons in the hidden layer, and the output of the hidden layer is continued to be transferred to the neurons of the output layer, and the output of the output layer offers us the results. The transfer, or the connections between the neurons in two consecutive layers, is done using weights, biases, and activation functions.

In detail, the input layer receives the information/data by storing the features/elements of input in the neurons. Each neuron $h$ in the hidden layer receives information/signals from all neurons $[x_{i}],i=1,N_{in}$ in the input layer, multiplies that with the weights $v_{i,h}$ and adds bias $b^{1}_{h}$, finally outputting $a_{h}$ after a nonlinear activation function (see Eq. (2)). In here we use the ReLU function which can improve the convergence rate of the estimation of weights and biases:

The output of the hidden layer $a_{h},h=1,\ldots,N_{h}$ is taken as the input of the output layer. With a similar procedure but using a different activation function, the information continue to be passed to the neurons in the output layer with the weights $w_{h,j}$ and the biases $b^{out}_{j}$. With a classification problem as in our work, the so-called softmax function is taken as the activation function, denoted as $f_{2}$, between the hidden layer and the output layer. After the transfer function $f_{2}$, the output of the output layer $Y=[y_{j}\epsilon[0,1]],\;j=1,\ldots,K$ (see Eq.(3)) is the predicted probability distribution for the classes:

For the goal of the ANN to classify asteroids using spectral data, the output $y_{j}$ of this ANN should be close to that of the known labeled class. If the weights and biases are properly fitted, the predicted values of the ANN should be consistent with the labeled values of the input data. Otherwise, these weights and biases involved in the hidden layer and the output layer need to be adjusted. In our case, there are $(N_{in}+1)N_{h}+(N_{h}+1)K=7300$ parameters (weights and biases) that need to be optimized to achieve accurate prediction of the input data classes. The procedure to find the best values for these involved parameters is also called a ’learning’ or ’training’ procedure.

The structure of our ANN provided, the back-propagation algorithm is applied to optimize the weights and biases. The cross-entropy measure (see Eq. (4)) is used as the cost function of our ANN (Boer et al., 2005), which measures the difference between the predicted values and the labeled values in the training data. Given the labels of input data, the smaller value of loss function means a better performance of this network. If there are $N_{sample}$ training samples included in a training set, the loss function can be written as:

where vectors $P=([p_{(i,n)}],i=1,\ldots,K;n=1,\ldots,N_{sample}$) and $Q=[q_{(i,n)}]$ are the labeled and predicted values, respectively. The subscript $(i,n)$ refer to class $i$ and sample $n$.

The learning procedure is actually optimization of the ANN where the weights and biases are adjusted in order to minimize the loss function. The optimization method we applied here is the mini-batch gradient descent (MBGD) algorithm (Goodfellow et al., 2016). It is one variant of the stochastic gradient descents (SGD), and has both the merits of the SGD and the batch gradient descent (BGD). The MBGD involved a small part of samples picked randomly from the whole training data set, called a batch of samples. With a batch of samples, the weights and biases are iteratively corrected until the cross-entropy reaches a minimum. The most-used approach in the optimization of the ANN is to move along the gradient direction of the loss function. Let a vector $z=[w,b]^{T}$ represent the parameters with $w$ as the weights and $b$ as the biases. The parameter vector is iteratively updated along to the negative gradient of the loss function with Eq. (5).

From Eq. (5), new step of update to the parameters $z_{k+1}$ is related to the previous values of parameters $z_{k}$ and the directions of gradient decrease of the loss function $g1,g2$. The quantity $\alpha$ ($>0$), called the learning rate, is usually set considering an efficient convergence of the training procedure.

Then, a new batch of samples is picked up from the rest of the training data set, and the parameters are continuously updated until each sample in the training data set is used. In practice, the entire training data set is randomly split into some batches according to the number of samples in a batch. The round when all batches are used in training the ANN is called a training epoch. Continuously, the training data set is split randomly and the next epoch of training of the ANN is done, until some conditions are satisfied, for examples, the variation of the loss function is smaller than a threshold or the number of epoch reaches the maximum limit. Based on our experience in training the ANN, the MBGD has a stable convergence speed, and is more efficient than the BGD.

The SMASS data are the first spectral data of asteroids obtained with the spectroscopic instrument using a CCD detection, so more features than the broad-band spectrophotometric colors are detected. The SMASS II database contains spectral data of 1447 asteroids covering a wavelength range of 0.435–0.925 µm. Each spectrum is composed of 182 data points with a sampling interval of 2.5 nm in wavelength.

To inherit from the Tholen taxonomy, Bus & Binzel (2002a) constructed a feature-based asteroid taxonomy with a PCA method based on the SMASS II data and information. It is called the Bus-Binzel or SMASS II taxonomic system, in which 26 classes (A, B, C, Cb, Cg, Cgh, Ch, D, K, L, Ld, O, Q, R, S, Sa, Sk, Sl, Sq, Sr, T, V, X, Xc, Xk and Xe) are identified from the spectra of 1447 asteroids.

The spectral data of SMASS II we used here is downloaded from a website¹¹1https://sbnapps.psi.edu/ferret/. Among the 1447 taxonomic asteroids in the SMASS II system, 1139 asteroids belong to C, S, and X-complexes, each of which are further divided into subclasses according to features in spectral data. The rest are identified as the special classes: A, B, D, K, L, O, Q, R, T, V. The occurrence of subclasses in the Bus-Binzel system reflect the connection or transition to its closely neighboring classes (Bus & Binzel, 2002a). We consider that it is also important to understand the spectral variations due to heterogeneous surfaces and due to different observational geometry and epoch.

As the first stage of the analysis to the spectral data using CSST and ground-based telescopes, 10 classes (A, B, C, D, K, L, S, T, V, and X-type) are included in our ANN to have enough samples in each class. The classification to subclasses of asteroids will be considered in future work. In all, 834 samples are picked from the SMASS II database. The detailed numbers of samples for each class are listed in Table 1. For convenience, we attributed integer numbers $0,\ldots,9$ to the classes we use here (A, B, C, D, K, L, T, V, X, and S-type).

The future slitless spectral data of the CSST will cover a wavelength range of 0.255 to 1.0 µm with a lower limitation of the resolution of 200 (Zhan, 2021). Before the construction of the ANN for the CSST spectra, the training data, meaning the SMASS II spectra, need to be re-sampled to match the CSST observations. Considering the shortest wavelength of $\lambda_{0}=0.435$ µm and the longest wavelength of 0.925 µm, the re-sampled data of SMASS II are composed of 151 data points(Resampling is completed according to expression (1)).

In the actually observed spectral data, as the wavelength increases, the number of data points will become more and more sparse. Using expression (1) to resample will be closer to the actual observation data. The reflectance value corresponding to each re-sampled wavelength is obtained by cubic spline fitting to the spectral data. Considering the wavelength interval of original spectra of SMASS II, we set 0.02 µm as the bin interval of cubic spline fitting, so at least 10 data points are involved in fitting each piece-wise spline function. Most re-sampled spectra are very close to the original ones, only some very noisy data are smoothed in the spline fitting procedure (see Fig. 3). The re-sampling procedure gives us 834 spectra, each having 151 features/channels.

From Table 1, it is easy to notice that the sample counts of different classes are not balanced, for example, the D-type has only 9 samples whereas the S-type has the largest number of samples, 382. If such data is used to train the ANN, the neural network could optimize and learn the features of the classes with large numbers of samples only. To balance the samples in different classes, we added some synthetic spectra into those classes with few samples. The synthetic spectra of a class are generated by a clone method based on all real data in that class.

The spectra of asteroids of the same class show similar shapes with some variation. To obtain the synthetic spectra, first, the representative spectrum of each class is found by averaging all re-sampled spectra in the class (the red line in Fig. 4). In detail, we calculated the mean reflectance at wavelengths $\lambda_{i}$ (i = 1, …, 151) and the corresponding standard deviations $\sigma^{i}_{k}$ (with $k$ as class index, $k=0,\ldots,9$) over all re-sampled spectra in the class $k$. Then, a clone spectrum of class $k$ is simulated by adding Gaussian noise with standard deviation of $\sigma^{i}_{k}$ into the representative spectrum of the $k$ class. Repeatedly, cloned spectra are generated until the total number of spectra in each class reaches 400. The cloned and real spectra for each classes are shown in Fig. 4. The $10\times 400$ spectra form the dataset used in training, validating, and testing our ANN in the next step. During the construction of the ANN, the dataset is divided into training, validation, and test data sets.

We have an input dataset of 4000 spectra to train and test the ANN, and a ratio of 4:1 is set to divide randomly the 4000 spectra into the training and test datasets. The training procedure of the ANN starts with initializing the weights and biases by random values, and then those weights and biases are optimized iteratively by the training procedure with the following steps: (1) separate the input dataset into training and test data set with the ratio 4:1; (2) split randomly the training data set into batches according to the batch size of 64; (3) correct the weights and the biases by the gradient descent method to minimize the loss function, in which the back-propagation method is used to implement the optimization; (4) repeat to update the weights and the biases with another new batch of data until all batches are used; that is one epoch of iteration to construct the ANN; (5) repeat (1)–(4) to reach the required number of epochs; (6) check and test the ANN performance with the test dataset.

Additionally, we introduced a dropout layer in our ANN to avoid an over-fitting problem. The ratio of dropout layer is set as 0.2, which means that random 20 % of neurons in the hidden layer are not connected to the neurons in the output layer in each round of training.

The optimization with gradient decrease (GD) family methods often finds a local minimum directly related to the initial values of the optimized parameters. To overcome this, five ANNs with different initial values for the parameters are constructed, and the final classification result to an input data is determined by voting among the five ANNs.

Using 834 real spectra of the SMASS II dataset, the classification performance of the trained ANN is assessed using the classification accuracy, the count ratio of correctly predicted cases to the entire sample. When using the softmax function as the activation function of the output layer, the ANN outputs a probability distribution of the taxonomic classes. The final predicted label of a spectrum is the class with the maximum probability in the probability distribution. To make the prediction result more robust, the vote result of five predicted labels suggested by five ANNs, initialized with different random numbers, is taken as the final predicted label for the input sample.

We investigated the accuracy of the classification ANN by comparing the final predicted labels for the 834 chosen real spectra from SMASS II to their labels in the Bus-Binzel taxonomic system. Generally, 764 out of 834 asteroids are correctly classified, which gives an accuracy of 92 %. Considering the accuracy of each class, the accuracies of A, B, L, and S-class are higher than 94 %, whereas the C, D, T, X, and K-classes have slightly lower accuracies. The detailed values for the 10 classes are listed in Table 3. From the confusion matrix (see Fig. 5), we can see where the problematic cases are. The diagonal elements of the confusion matrix are the number of samples with predicted labels consistent with that in the Bus-Binzel taxonomy, and conflicting labels are the off-diagonal elements of the matrix. It is easy to note that the top conflicting labels are between the C and X-class, and between the S and K-class. For example, 25 C-class asteroids are predicted to be X-class, and 14 S-class asteroids were labeled as K-class with our ANN. Among the 25 C-class samples in the Bus-Binzel system labeled as X-class by our ANN, we found 9 to show a gentle convex behaviour between 0.6–0.8 µm, indicating that they might belong to the Xc-class or the transition between the C and X-classes. The 14 S-class samples in the Bus-Binzel system are predicted as K-class by our ANN probably because the classes show very similar spectral shapes.

As additional test of our ANN, we selected spectra of 415 asteroids from the Small Solar System Objects Spectroscopic Survey (S3OS2) (Lazzaro et al., 2007, for short, S3OS2), which are obtained with the 1.52-m telescope of the European Southern Observatory during the period of November 1996 and September 2001. The wavelength of the S3OS2 data range 0.49-0.92 µm. Lazzaro et al. (2007) classified the S3OS2 data refering to the Bus taxonomy system, the analysis results denoted here as S3OS2(B). Among the 415 samples picked up from the S3OS2 dataset, 80 are of S-type, 205 of B, C, or X-type, and 130 of A, D, K, L, T, or V-type.

With a similar procedure, those test data are firstly re-sampled according to the format of the input data of our ANN. Then, they are analyzed with our ANN. If comparing our classification results to those in S3OS2(B), 362 out of 415 samples have the same taxonomy labels. The ratio of the matched samples to the whole data is $\sim$ 87 %. In detail, the accuracies for S-type, the B/C/X-type group, and the A/D/K/L/T/V-type group are 0.89, 0.87, and 0.86, respectively.

The five nights of spectroscopic observations were carried out in three nights in 2006 (December 26th, 27th, and 28th) and in two nights in 2007 (November 17th and 19th). The spectral data of the selected asteroids were obtained by an OMR cardigan spectrometer and a PI $1340\!\times\!400$ CCD detector. A grating of 300 lines/mm of the OMR and a long slit of 2.5” in width gives a resolution of $R=200$. The orientation of the long slit of the OMR instrument is along the South-North direction. At each night of observations, at least two solar analog stars (selected from HR4030, HR3538, HR3951, HD28099, HD191854, and HR996) were observed one or two times. For each observation of asteroids and solar-like stars, one spectrum of a HeAr lamp installed in the OMR system was followed for the wavelength calibration. In order to eliminate the effect of atmospheric extinction, we arranged the spectroscopic observation of asteroids at the smallest possible zenith distance.

Because of the motion of an asteroid relative to the stars in the same field of view, it is possible for the asteroid to closely pass by a star, and that could result in stellar light entering the slit, which would lead to contamination of the spectral data of the asteroid. We checked the observation images of asteroids, and contaminated data are picked up and rejected. Also, part of the spectral data on Dec. 28th, 2006 are not used due to the abnormal operation of the OMR. Finally, 64 spectra of 42 asteroids are involved in the analysis. In Table 4, the detailed observational information, including the time of observation, exposure time, airmass, and the solar analog stars used, is listed for each spectrum.

The data reduction of the spectroscopic observations of asteroids were done according to the standard procedure (Bus et al., 2002) with the IRAF software. Firstly, the systematic errors of the spectroscopic image, i.e., bias, dark current, and flat-field effects were corrected, then cosmic rays in the image were identified by a threshold of 4 times the standard deviation of the sky background and removed. Secondly, one-dimensional spectra of asteroids, solar-like stars, and HeAr lamp were extracted by an optimal aperture and a fitted background. The wavelength calibration of a target’s spectrum was done with the aid of the spectrum of the HeAr lamp after the target. As for the atmospheric extinction effect, we used an average extinction curve of the Xinglong site (Code 327).

At each observational night, multiple bias frames were obtained at the beginning and end of observation of that night, and were combined into a synthetic bias image by averaging, called the best bias image. All other spectroscopic images, i.e., those of the targets, flats, and HeAr lamp, were subtracted with the best bias image to correct the bias effect of the CCD detector. Similarly, a combined flat image was generated from multiple observed flat images which actually are spectroscopic images of an arc lamp in the telescope. The illumination and response effect on the combination flat was simulated and then removed from the combination flat image. The final flat image is a normalized combination flat image without the illumination and response effect. The flat-field effect of the CCD detector on all spectroscopic images are corrected by dividing with the final flat image.

When extracting one-dimensional spectra of celestial objects, a trace procedure was applied. That is to trace the peak flux along the dispersion direction in the two-dimensional spectroscopic image with an optimal width (or aperture), after which all flux within the aperture are summed along the vertical dispersion direction to form the one-dimensional spectrum. The optimal aperture of the extraction of the one-dimensional spectrum is determined by balancing the maximum target flux and the minimum of background flux in the aperture.

For the wavelength calibration of a target spectrum, we used a Legendre function with the order of three or five. The calibration is derived by fitting the pixel positions of the emission lines for the spectrum of the HeAr lamp following the target to their wavelength values. On average, the wavelength range of our spectra is between 0.40 and 0.83 µm.

To reduce the influence of extinction as much as possible, we scheduled the observation epochs for each asteroid for the time when they are as high in zenith as possible. In practice, we applied an average extinction curve of the Xinglong site of National Astronomical Observatory: the correction of extinction for each source spectrum is done according to the airmass of the time of the observation.

To obtain the reflectance of an asteroid, we conducted spectroscopic observations for the selected solar analog stars in the same night of the asteroid observations. Both the spectra of the asteroids and the solar analog stars are normalized to unity at the wavelength of 0.55 µm, then the reflectance of an asteroid is that of the normalized spectrum of the asteroid divided by the normalized spectrum of the selected solar analog star. Sometimes multiple spectra of single solar analog star are obtained during one night, the averaged spectrum of those is then used to derive the reflectance of the asteroids. If spectra of multiple solar analog stars are obtained, we prefer to use the spectrum of a relatively stable star, or the averaged spectrum of those stars with no obvious activity.

The sample of 55 reflectance spectra for 41 asteroids is shown in Fig. 6 and an example with 9 reflectance spectra of asteroid (469) Argentina is shown in Fig. 7.

According to the format of the input data of our ANN tool (151 channels ranging from 0.43 to 0.92 µm with an increase of the span of $(1+\frac{1}{R})^{i}$ $(i=1,\ldots,151)$, the 64 spectra of 42 asteroids are re-sampled using the cubic spline algorithm. Because the wavelength range of our data is around between 0.40 and 0.83 µm, the shortage of data beyond 0.83 µm is complemented by linearly extrapolated values.

These 64 re-sampled spectra are analysed with our ANN tool. The classification results for the 64 spectra of the 42 observed asteroid are listed in Table 5. Statistically, among those 64 spectra, 17 are predicted as C-class, 13 as B-class, 11 as X-class, 14 as T-class, 7 as D-class, and S and K-classes have one sample each.

Among the 42 observed asteroids, 27 have only a single spectrum, 14 have two spectra, and (469) Argentina has nine spectra. The nine spectra of (469) Argentina obtained in three adjacent nights in 2006, obtained T-class, except one of a label of C-class (see Fig. 7). For the case of two spectra, spectral data of ten asteroids (34, 35, 199, 381, 567, 613, 680, 772, 780, and 1303) were obtained in the same night, while four asteroids (95, 326, 696, and 1004) were obtained in two adjacent nights. The analysis results with our ANN tool show that asteroids with two spectra obtained consistent labels except for asteroids (696) and (780). These two asteroids obtained C and X-class spectra (see Fig. 8).

We compared our results of observed targets to their existing classification (see Column 5 in Table 5). Among the 42 observed asteroids, 31 have taxonomic type in the Tholen system (Tholen, 1984) and/or the Bus-Binzel system (Bus & Binzel, 2002a). For these 31 labeled asteroids, 24 are of C-type or sub-class of the C-complex in Tholen and/or Bus-Binzel system, three are of E/M/P-type and another three are of P-type, and one is X-type.

Asteroid (199), labeled as X-type in the Bus-Binzel system, is classified as T-type with our ANN tool based on its two spectra. Three P-type asteroids (613, 1512, and 1911) in the Tholen system are identified here as D-class (613 and 1911) and T-class (1512). For the three E/M/P-type asteroids in the Tholen system, (469) is identified as T-class, (514) as K-class, and (663) as S-class. Considering the large albedo of E/M/P-type asteroids, they usually have a moderate spectral slope in the visible band. The new labels for the three asteroids are reasonable.

The 24 asteroids labeled as C-type or sub-class of C in the Tholen system are classified as C, B, and X-class (in detail, 8 asteroids as C, 11 as B, 3 as X and 2 as T) by our ANN tool. Checking the data of the Bus-Binzel taxonomic system, we found that 16 of the 42 observed asteroids are included, and fourteen asteroids were labeled as sub-class in the C-complex (C, Ch, Cb, or Cg-type), and two asteroids (199 and 907) as X-complex. For these fourteen C-complex asteroids, six are classified as B-type, five as C-type, two as X-type, and one as T-type by our ANNs tool. The asteroid (199) is classified as T-type and the asteroid (907) as B-type by our ANN tool.

Generally, asteroids of C-type in the Tholen system or C-complex in the Bus-Binzel system are classified as one of the C, B, X, and T-types with our ANN tool. We think this is due to the featureless spectra of the B, C, X, and T-types. Some features, i.e., the slope between 0.55 and 0.86 µm, shallow absorption in UV, absorption around 0.70 µm, or maybe absorption around 0.85 µm (see Fig. 4) are key points to distinguish the above four types. The reasons of confusion between the four types could be the low signal-to-noise ratio in the UV-band and the lack of data above 0.83 µm of our data. Anyway, we could not rule out the spectral variations of an asteroid rising from heterogeneous composition over its surface and/or change of the observational geometry, e.g., solar phase angle.

In our data, 14 asteroids have two spectra, and one asteroid has 9 spectra obtained in three adjacent nights. For the aim of classification of asteroids, we built five ANNs to determine the type for each spectra by voting. To check the final type of each spectrum for those asteroids with more than one spectrum, 13 asteroids received a consistent type, two asteroids have different labels (see Table 5). To understand the subtle spectral variation of an asteroid at different rotational phase, we choose one ANN from our five ANNs which gives the most similar result as compared to the voted result. Here, we compared the maximum probability output by the chosen ANN of each spectrum of asteroids with multiple spectra.

For the 14 asteroids with two spectra for each, 12 asteroids are classified as the same type for their two spectra, while two asteroids (696) and (780) received different labels (see Fig. 8). By the maximum probability that a spectrum belong to a certain type ( the column 3 of Table 6),the extent of variation of spectra of an asteroid, or different asteroids but belong to the same type, could be investigated. As for the two asteroids (696) and (780), we checked the extinction effect on their observations. We found that the airmasses for the two spectra of (696) in the two nights are 1.64 and 1.22, respectively, and that the airmasses of (780) in the same nights are large (2.61 and 2.63). So, we could not rule out the possibility that the diverse labels of the two asteroids rise from the atmospheric effects, although the extinction correction has been done with the averaged extinction curve of the Xinglong station.

For the asteroid (469) Argentina, we compared its nine spectra obtained in three adjacent nights. Considering the maximum probability of each spectra, eight spectra are labeled as T-type with high probabilities, and one spectrum as C-type (see Table 6). The photometric data of Argentina show a more complicated lightcurve shape than that of a sinusoid with two peaks, implying multiple potential periods (17.57 h, 12.84 h, and 8.79 h) of brightness variations (Colazo et al., 2021; Wang et al., 2005; Warner, 2007). Here, we selected a spin period of 8.79 h (this period gives a double-peaked lightcurve over a full rotation) to calculate the rotational phase (zero phase at JD2454096.14340). Figure 7 shows nine spectra of Argentina arranged from bottom to top with increasing rotational phase angle. From Fig. 7, the most significant variation occurs at the phase of 0.067 deg, and spectral variations at other rotational phases could be related to the values of maximum probability of each spectrum in the ANN’s output layer (see Table 6).