Visualizing convolutional neural network for classifying gravitational
  waves from core-collapse supernovae
In this study, we employ a convolutional neural network to classify
gravitational waves originating from core-collapse supernovae. Training is
conducted using spectrograms derived from three-dimensional numerical
simulations of waveforms, which are injected onto real noise data from the
third observing run of both Advanced LIGO and Advanced Virgo. To gain insights
into the decision-making process of the model, we apply class activation
mapping techniques to visualize the regions in the input image that are
significant for the model's prediction. The class activation maps reveal that
the model's predictions predominantly rely on specific features within the
input spectrograms, namely, the $g$-mode and low-frequency modes. The
visualization of convolutional neural network models provides interpretability
to enhance their reliability and offers guidance for improving detection
efficiency.
Seiya Sasaoka, Naoki Koyama, Diego Dominguez, Yusuke Sakai, Kentaro Somiya, Yuto Omae, Hirotaka Takahashi
Visualizing Convolutional Neural Network for Classifying Gravitational Waves from Core-Collapse Supernovae

###### Abstract

In this study, we employ a convolutional neural network to classify gravitational waves originating from core-collapse supernovae. Training was conducted using spectrograms derived from three-dimensional numerical simulations of waveforms, which were injected onto real noise data from the third observing run of both Advanced LIGO and Advanced Virgo. To gain insights into the model's decision-making process, we apply class activation mapping techniques to visualize the regions in the input image that are significant for the model's prediction. Our model distinguished between 9 different waveforms and noise with an accuracy of 98.4% at 1 kpc. Visualization through class activation mapping revealed that the model's predictions predominantly rely on specific features within the input spectrograms, namely the g-mode and low-frequency modes.

## I Introduction

The first detection of gravitational waves (GWs) from a binary black hole merger by the Advanced Laser Interferometer Gravitational-wave Observatory (Advanced LIGO) [1] in 2015 marked the beginning of GW astronomy [2]. Throughout three observing runs, O1, O2, and O3, Advanced LIGO and Advanced Virgo [3] reported 90 GW events [4; 5; 6; 7]. As of May 2023, the international GW network, now including KAGRA [8], has begun its fourth observing run (O4) with improved sensitivity.

All the GW events detected so far are exclusively from compact binary coalescences. However, short-duration GW bursts arising from core-collapse supernovae (CCSNe) are expected to be detected by the current and the next-generation GW detectors, such as the Einstein Telescope [9] and the Cosmic Explorer [10]. CCSNe, resulting from massive star explosions leading to neutron stars or stellar-mass black holes, stand as one of the most energetic astrophysical events in the universe, emitting electromagnetic waves, neutrinos, and GWs. While electromagnetic waves from CCSNe are frequently observed, neutrinos have only been detected from SN1987A [11; 12]. GWs are expected to carry information about the inner core's dynamics, providing vital insights into the explosion mechanism that remains elusive. The primary conundrum lies in discerning how a stalled shock wave is revived to cause a star to explode. Currently, there are two prevailing theories [13]: the neutrino-driven mechanism [14], in which shock waves are revived by beneath-stored neutrinos heating the surrounding matter, and the magnetorotational mechanism [15], in which rapid rotation of the progenitor causes explosions driven by strong magnetic fields. The typical GW detection range for neutrino-driven signals is expected to be around 10 kpc, while the detection range for magnetorotational signals is expected to be above 100 kpc [16].

Due to the stochastic nature of GW signals from CCSNe, the conventional matched filtering technique, which relies on specific waveform templates, is unsuitable. Alternative detection methods based on time-frequency representation have been devised in response. In particular, the coherent WaveBurst (cWB) pipeline [17; 18] detects and reconstructs burst GW signals by searching for excess power in a time-frequency map, with minimal reliance on a specific source model.

Predicting GW signals from CCSNe remains a formidable challenge. However, recent advancements in theoretical research and multi-dimensional numerical simulations have revealed certain signal properties. For neutrino-driven CCSNe, the dominant emissions arise from the g-mode oscillation of the proto-neutron star (PNS) surface. These frequencies progressively increase over time, ranging from a few hundred Hz to a few kHz. Additionally, at low frequencies (\(\lesssim\) 200 Hz), GW emissions associated with hydrodynamics instabilities including neutrino-driven convection and standing accretion shock instability (SASI) [19] are observed in some simulations. Insights obtained from these simulations are pivotal in enhancing methods for CCSNe detection and
analysis.

In recent years, machine learning techniques, especially deep learning, has gained traction in a variety of scientific fields due to its capacity for recognizing intricate patterns and extracting meaningful features from large datasets. This ability has been especially noted in areas such as computer vision and natural language processing. Its application in GW research has followed, with numerous implementations and explorations, as highlighted in a comprehensive review by Ref. [20] and the foundational efforts by George and Huerta [21; 22]. In the field of CCSNe analysis, Astone _et al._[23] leveraged convolutional neural networks (CNNs) to detect CCSNe within Gaussian noise, using g-mode phenomenological waveforms for training, outperforming the cWB pipeline. Consecutive studies by Iess _et al._[24; 25] involved the training of both one- and two-dimensional CNNs and long short-term memory (LSTM) networks [26] to identify seven distinct CCSN waveforms embedded in real noise and glitches, with their models achieving 98% classification accuracy at 1 kpc with a three-detector network. Additionally, Chan _et al._[27] employed one-dimensional CNNs to investigate both magnetorotational and neutrino-driven signals in Gaussian noise, recording a true alarm probability of 80% for magnetorotational signals from sources at 60 kpc and 55% for neutrino-driven signals from sources at 10 kpc with a fixed false alarm probability of 10%. In another study, Edwards [28] used two-dimensional CNNs to classify 18 different equations of state (EOS) from pure magnetorotational CCSN signals, attaining an accuracy of 72%. Lopez _et al._[29] refined phenomenological waveforms originally used by Astone _et al._[23], achieving a 60% true alarm probability for signals located at 15 kpc with a 5% false alarm rate.

Although deep learning exhibits strong performance on a wide range of tasks, its intricate models, characterized by a large number of parameters, pose challenges in elucidating their decision-making processes. To address this, the field of explainable artificial intelligence (XAI) [30] has surged, aiming to make model decisions transparent and interpretable. Within the context of CNNs, efforts have been made to develop techniques that attempt to understand the decision-making process by reverse-mapping the output of the network into the input space to identify the specific input components that were discriminative in producing the output. Class activation mapping (CAM) [31] is one such method, which computes a weighted sum of the outputs of the last convolutional layer using the outputs of the global average pooling layer after the last convolutional layer as weights. It helps identify the regions in the input image that were important for a prediction, but the model needs to be modified to include a global average pooling layer, which may result in lower accuracy. Gradient-weighted class activation mapping (Grad-CAM) [32] was introduced as a solution to this limitation of CAM, offering the advantage of not requiring any modifications to the network architecture by using gradient information from the prediction for weighted parameters. Subsequently, Grad-CAM++ [33], a generalization of Grad-CAM, and Score-CAM [34], a gradient-free CAM method, were developed to generate more accurate saliency maps than Grad-CAM. These techniques to analyze deep learning models are commonly used in fields such as electrocardiogram signal analysis [35] and X-ray diagnosis [36], however, for GW analysis, they have only been used in Ref. [37] to the best of our knowledge.

In this study, we first take an approach similar to Ref. [25] and train a two-dimensional CNN model to classify CCSNe signals using short-time Fourier transformed spectrograms as input for simplicity. We use nine types of waveforms from recent three-dimensional numerical simulations and O3 real noise to train and validate our model. In the test, signals from sources between 1 and 10 kpc are considered, and the performance of the model for sources at each distance is discussed. To interpret the model, we use three CAM methods to generate saliency maps and evaluate them using two metrics: average drop and average increase. The best CAM method is then applied to correctly classified and also misclassified samples to visualize the regions in the input spectrogram that influence the predictions of our model.

The remainder of this paper is organized as follows. Section II describes our datasets, the CNN model, and the CAM techniques. In Sec. III, we discuss the classification performance of our model, and apply multiple visualization techniques to interpret the model. We summarize and conclude the paper in Sec. IV.

## II Method

Our CNN model is trained to classify strains at three detectors LIGO Hanford (H1), LIGO Livingston (L1), and Virgo (V1) into 10 classes: noise and 9 different CCSN waveforms. In this section, we first provide an overview of the data used in this study, including the brief summary of the CCSN simulation data and the pre-processing strategy to generate our training, validation, and test sets. Subsequently, our CNN architecture and the theory of visualization technique of the model are explained.

### Dataset

#### ii.1.1 CCSN Waveforms

Modelling the stellar core collapse, bounce and the subsequent post-bounce evolution is very complicated and computationally expensive. However, remarkable advancements in three-dimensional numerical simulations of neutrino-driven explosions have been achieved by several groups in recent years. Specific details of the waveform depend on various properties of the progenitor, such
as mass, angular velocity, and EOS of the dense matter. Both the general relativity approximation and the handling of neutrino transport critically influence simulations. From the available simulation data under a variety of conditions, we selected nine types of waveforms from four recent three-dimensional numerical simulations [38; 39; 40; 41]. All of them allow us to compute the GW amplitude in any observer direction from the quadrupole moment.

Powell and Muller 2019 [38] performed simulations using the general relativistic neutrino hydrodynamics code CoCoNuT-FMT [42]. We used two waveforms from the models he3.5 and s18. The progenitor of he3.5 is an ultra-stripped star evolved from a helium star with an initial mass of 3.5 \(M_{\odot}\). The simulation is stopped at 0.7 s after core bounce. The GW is dominated by excitation of g-modes in the PNS with a peak frequency around 900 Hz. Model s18 is a single star with a zero-age main-sequence (ZAMS) mass of 18 \(M_{\odot}\). The simulation was stopped 0.89 s after core bounce. The GW emission is similar to model he3.5, with g-modes oscillations of the PNS with a peak frequency around 900 Hz.

Radice _et al._ 2019 [39] studied eight models using the Eulerian radiation-hydrodynamics code FORNAX [43]. We used waveforms from the models s13 and s25 corresponding to progenitors of 13 and 25 \(M_{\odot}\) ZAMS, respectively. The simulation is ended at 0.77 s in s13 and 0.62 s in s25 after bounce. Both waveforms are characterized by f- and g- modes with a peak frequency around 1400 Hz in s13 and 1100 Hz in s25. In addition, s25 waveform has a clear SASI mode around 100 Hz.

From the simulation by Powell and Muller 2020 [40], we used three models m39, y20, and s18np. Model m39 is a rapidly rotating 39 \(M_{\odot}\) Wolf-Rayet star with an initial surface rotation velocity of 600 \(\mathrm{kms}^{-1}\). It produces a neutrino-driven explosion without magnetic fields. The other two are non-rotating models of 20 \(M_{\odot}\) Wolf-Rayet star and an 18 \(M_{\odot}\) ZAMS star. The simulation is ended at 0.98 s, 1.2 s, 0.56 s after core bounce in models m39, y20, and s18np, respectively. All of the three models show GW emission associated with prompt convection shortly after bounce as well as f-mode oscillations of the PNS. In model s18np, the absence of strong perturbations from convective oxygen burning, in contrast to s18, prevents the shock from being revived and leads to the development of strong SASI activity, with a frequency reaching \(\sim\)400 Hz by the end of the simulation.

Powell _et al._ 2021 [41] performed simulations using three different EOS, LS220 [44], SFHx and SFHo [45]. The progenitor models are 85 and 100 \(M_{\odot}\) Population III ZAMS stars. We used z85_sfhx and z100_sfho models. The simulation is ended at 0.59 s in z85_sfhx and 0.62 s in z100_sfho model. Both waveforms show typical g-mode emission with a peak frequency of \(\sim\)700 Hz. In z85_sfhx model, the frequency of the SASI emission increases up to the point of shock revival reaching \(\sim\)200 Hz and decreases afterwards. In z100_sfho model, which does not explode, the frequency of the SASI emission continues to increase by the end of the simulation reaching \(\sim\)400 Hz.

Figure 1 shows the amplitude spectral density of the plus mode of each waveform at 1 kpc from the polar direction. The amplitude of the m39 waveform is the largest of these waveforms. Waveforms that have peaks around 100 Hz such as the s25 and the s13 indicate that they have SASI induced GW mode.

#### ii.1.2 Data Processing

From the simulation data presented in the previous section, we calculate the amplitude of the GW for generating the datasets. The directions of radiation \((\theta,\phi)\) are uniformly sampled and the plus and cross polarizations of each GW are calculated using the formulae

\[h_{+} =\frac{1}{D}\frac{2G}{c^{4}}(\ddot{Q}_{\theta\theta}-\ddot{Q}_{ \phi\phi}), \tag{1}\] \[h_{\times} =\frac{1}{D}\frac{G}{c^{4}}\ddot{Q}_{\theta\phi}, \tag{2}\]

where \(Q\) is the traceless quadrupole moment, and \(D\) is the distance between a source and Earth. As the sampling of simulation is usually not uniform in time, we resampled data uniformly with a sampling rate of 4096 Hz. A high pass filter with a cutoff frequency of 11 Hz and a Tukey window with \(\alpha=0.1\) are applied to the resampled signals. Each signal is then truncated or padded with zeros to make the length one second. In order to make the model robust, we randomly time-shifted the signals so that the time of core bounce is between 0 and 0.15 s. For the training and validation sets, the signals are scaled using optimal matched filter signal-to-noise ratio (SNR), as the amplitudes of the simulated signals are

Figure 1: Amplitude spectral density of the plus mode of each waveform at 1 kpc. The observer is in the polar direction.
quite different. SNR is defined as

\[\rho=\sqrt{4\int_{f_{\rm{min}}}^{f_{\rm{max}}}\frac{|\tilde{h}(f)|^{2}}{S_{n}(f)} \mathrm{d}f}, \tag{3}\]

where \(\tilde{h}(f)\) is the Fourier transform of the signal and \(S_{n}(f)\) is the one-side power spectral density of the noise. The network SNR of the detectors H1, L1, and V1, given by

\[\rho_{\rm{net}}=\sqrt{\rho_{\rm{H1}}^{2}+\rho_{\rm{L1}}^{2}+\rho_{\rm{V1}}^{2}}, \tag{4}\]

is used to scale the signals. We generated samples with network SNRs from 20 to 50 for training and validation sets. For the test set, the signals are scaled to have the distances between 1 and 10 kpc. Sky location is also randomly selected and GW amplitude \(h(t)\) is computed, taking into account the antenna pattern functions \(F_{+}\) and \(F_{\times}\), and the delay in arrival time of each detector with the following equation:

\[h(t)=F_{+}(\alpha,\delta,\psi,t)h_{+}(t+\Delta t)\] \[\qquad\qquad+F_{\times}(\alpha,\delta,\psi,t)h_{\times}(t+\Delta t), \tag{5}\]

where \(\alpha\) is the right ascension, \(\delta\) is the declination and \(\psi\) is the polarization angle. \(\Delta t\) is the delay in arrival time between the detector and the center of Earth. We used PyCBC software library [46] to carry out these computations.

Noise used in this study are O3 real data of Advanced LIGO and Advanced Virgo, obtained from Gravitational Wave Open Science Center (GWOSC) [47]. Data from GPS time 1238265720 to 1238252308 was used for the training set, 1238265720 to 1238354855 was used for the validation set, and 1238404064 to 1238457121 was used for the test set. Data around the event time reported in the second Gravitational Wave Transient Catalog (GWTC-2) [6] are excluded. After a signal is injected in noise, each sample is whitened with the power spectral density computed using Welch's method [48] and then short-time Fourier transformed with a window size of 0.0625 seconds to produce a spectrogram. The spectrogram is normalized to [0, 1] before input to the network.

We generated 60,000 samples as each of training and validation sets, and 100,000 samples as test set. The test set has 1,000 samples for each class and each distance. Sample spectrograms in the training set are shown in Fig. 2.

### CNN Model

Our CNN model consists of two convolutional layers of kernel size 3, each followed by a max-pooling layer of size 2 and a rectified linear unit (ReLU) layer. The outputs of these layers are fed into two fully connected layers, and finally the softmax layer outputs a size 10 vector whose elements represent a probability of each class. The model has 427,378 trainable parameters in total. This model is shallower than the one used in Ref. [25]. However, its classification performance is comparable to the previous study, prompting us to adopt this model. Reducing the number of layers also helps us generate higher-resolution CAM maps.

The model is trained using categorical cross entropy as the loss function and Adam optimizer [49] with a learning rate of \(5\times 10^{-4}\) to update the weights. In the training, we adopt curriculum learning [50] as a strategy to enhance the model and accelerate the training by starting from inputting high SNR samples and gradually adding lower

Figure 2: Sample whitened spectrograms of each class at the H1 detector in the training set. Each signal sample is observed in the polar direction and scaled to have an SNR of 40. The bounce time is fixed at 0.1 s.
SNR samples. We trained the model on a single GPU (NVIDIA GeForce RTX3090) for 120 epochs with a mini-batch size of 128.

### Visualization

After training the model, we use CAM techniques to generate saliency maps. These maps show the regions in the input that influenced the model's prediction. In this study, we selected three CAM methods: Grad-CAM, Grad-CAM++, and Score-CAM, which are widely used today to interpret CNN models. All of these CAM techniques are applied to the convolutional layer prior to the final max-pooling layer in our model.

#### ii.3.1 Grad-CAM

Grad-CAM is a gradient-based visualization technique that highlights the important regions of an input image that the model is looking at while making a prediction. Suppose that for a given input, the prediction score for class \(c\) before the softmax layer of the trained model is \(y^{c}\), and the \(k\)-th output matrix of the last convolutional layer is \(A^{k}\). To obtain the Grad-CAM map of class \(c\), we first compute the gradients of the score \(y^{c}\) with respect to the \((i,j)\) component of the \(k\)-th feature map \(A^{k}\). We then take the global average of these gradients:

\[\alpha_{k}^{c}=\frac{1}{Z}\sum_{i,j}\frac{\partial y^{c}}{\partial A_{ij}^{k}}, \tag{6}\]

where \(Z\) is the number of pixels in \(A^{k}\). This weight \(\alpha_{k}^{c}\) represents the importance of the feature map \(k\) for the class \(c\).

The Grad-CAM map of the class \(c\) is computed as linear sum of \(A^{k}\) with \(\alpha_{k}^{c}\) as weights. The ReLU function is applied to extract only features that have a positive contribution to the prediction score. The resulting map of class \(c\) is expressed as

\[L_{\text{Grad-CAM}}^{c}=\text{ReLU}\Bigg{(}\sum_{k}\alpha_{k}^{c}A^{k}\Bigg{)}. \tag{7}\]

Since convolutional layers and pooling layers make the size of the feature map smaller than the input, Grad-CAM map is finally interpolated to make it the same size as the input.

#### ii.3.2 Grad-CAM++

While Grad-CAM takes a global average of the gradient matrix when calculating the weight \(\alpha_{k}^{c}\) in Eq. (6), Chattopadhay _et al._[33] proposed a method to fully include the importance of each pixel in the gradient matrix by taking its weighted average for the weight:

\[\alpha_{k}^{c}=\sum_{i,j}\alpha_{ij}^{kc}\ \text{ReLU}\Bigg{(}\frac{\partial y^{c}} {\partial A_{ij}^{k}}\Bigg{)}. \tag{8}\]

The ReLU function is used to account for features that increase the activation of the output neuron rather than suppress the activation of the output neuron. The weights \(\alpha_{ij}^{kc}\) can be theoretically derived using higher-order derivatives:

\[\alpha_{ij}^{kc}=\frac{\frac{\partial^{2}y^{c}}{(\partial A_{ij}^{k})^{2}}}{2 \frac{\partial^{2}y^{c}}{(\partial A_{ij}^{k})^{2}}+\sum_{a,b}A_{ab}^{k} \frac{\partial^{3}y^{c}}{(\partial A_{ij}^{k})^{2}}}. \tag{9}\]

This method is known as Grad-CAM++, since it can be considered as a generalization of Grad-CAM. The saliency map for Grad-CAM++ is expressed in the same way as for Grad-CAM, using weights in Eq. (8) and feature maps, as

\[L_{\text{Grad-CAM++}}^{c}=\text{ReLU}\Bigg{(}\sum_{k}\alpha_{k}^{c}A^{k} \Bigg{)}. \tag{10}\]

#### ii.3.3 Score-CAM

Wang _et al._[34] proposed a gradient-free CAM method called Score-CAM. It solves the problem of gradient-based CAM methods, namely the gradient is unstable, easily disturbed by noise, and can vanish or explole in deep networks. To generate a Score-CAM map, feature maps are used to mask an input image. Let \(H^{k}\) be the \(k\)-th feature map, up-sampled to the same size as the input and normalized to [0, 1]. Given an input image \(X\), the weight for the \(k\)-th feature map is computed as the difference between the score of the masked image \(X\circ H^{k}\) and the score of the baseline image \(X_{\text{b}}\):

\[\alpha_{k}=f(X\circ H^{k})-f(X_{\text{b}}), \tag{11}\]

where \(f(\cdot)\) denotes the output of the CNN and \(\circ\) denotes the Hadamard product. A black image is used as a baseline image. The Score-CAM map of class \(c\) is then computed as linear sum of the \(c\)-th value of \(\alpha_{k}\) and the feature map \(A^{k}\) as

\[L_{\text{Score-CAM}}^{c}=\text{ReLU}\Bigg{(}\sum_{k}\alpha_{k}^{c}A^{k}\Bigg{)}. \tag{12}\]

## III Results and Discussion

### Classification Performance

Classification accuracy is defined as the proportion of correctly classified samples out of the total number of
samples. After training, the model achieved a classification accuracy of 97.8% on a validation set consisting of uniformly sampled signals with SNRs between 20 and 50. On the test set, our model shows an accuracy of 98.4% for signals with sources from 1 kpc, which is comparable to the results of the previous study [25], despite some differences in the condition that we used O3 noise instead of O2 noise and performed 10-class classification instead of 8-class classification. In Fig. 3, we plot a true positive rate (TPR) for each waveform in the test set against distance. A TPR, also known as the sensitivity of a class \(c\), is defined as the ratio of the number of samples correctly classified into class \(c\) to the number of samples of class \(c\) in the test set. For signals from sources at 1 kpc, each waveform has a TPR greater than 90%, and this decreases monotonically with distance of the source, having an average TPR of 26.1% at 10 kpc. For the m39 waveform, because the amplitude of the strain is much larger than the others due to its rapid rotation and high explosion energy, the TPR for sources at 10 kpc is 99.2%.

The performance of a multi-class classifier is also expressed by a confusion matrix, which shows the number of samples classified into each class. Figure 4 plots the confusion matrices normalized for each class and the distribution of the network SNR for signals from sources at 1, 5, and 10 kpc. We can see from the confusion matrices that as the distance increases, the amplitude of the signal becomes smaller and the number of samples misclassified as noise increases. The accuracy for signals at 10 kpc is 33.2%, and our model cannot identify most of these signals, except for the m39 waveforms, whose SNR is much higher than others with a median value of 47.9.

### Dimensionality Reduction

Before implementing CAM techniques, we used the t-distributed Stochastic Neighbor Embedding (t-SNE) [51] algorithm to see if the convolutional layers in the model could extract the features in the input to classify samples. The t-SNE algorithm is a dimensionality reduction technique that minimizes the Kullback-Leibler divergence between two probability distributions: one representing pairwise similarities between data points in the original high-dimensional space and another representing pairwise similarities in a lower-dimensional space. In our CNN model, each sample is compressed into a vector with a length of 2112 before the dense layers. The t-SNE algorithm is used to map this vector into two-dimensional space to make it interpretable for humans. We visualized the dimensionally reduced feature maps of the test set, whose signals are coming from sources at 1 kpc, for which our model showed a good classification accuracy. The visualized data are shown in Fig. 5. We can clearly see that there are 10 clusters in the dataset and our model could extract meaningful features to classify these samples into 10 classes. The fact that some signal samples are also found in the noise cluster and that s13 samples are found in other clusters, especially in the noise cluster, is consistent with the results of the confusion matrix in Fig. 4(a).

### Saliency Maps

To quantitatively evaluate different CAM methods, we use two metrics, average drop and average increase [33], which focus on the change in a model's score caused by the explanation map. An explanation map for a target class \(c\) is generated as element-wise multiplication of a saliency map \(L^{c}\) with an original image \(X\):

\[E^{c}=L^{c}\circ X. \tag{13}\]

Average drop measures the percentage decrease in a model's score for a target class \(c\) when inputting only the explanation map, instead of the original image. It is expressed as

\[\text{Average drop}=100\cdot\frac{1}{N}\sum_{i}\frac{\max(0,y_{i}^{c}-o_{i}^{ c})}{y_{i}^{c}}, \tag{14}\]

where \(y_{i}^{c}\) is the score for class \(c\) on \(i\)-th original image and \(o_{i}^{c}\) is the score on the explanation map. The lower this value, the more effective the visualization method, since the explanation map includes more of the relevant information for making a correct prediction.

Average increase measures the number of samples in the dataset, the model's confidence increased when providing only the explanation map as input. It is expressed as

\[\text{Average increase}=100\cdot\frac{1}{N}\sum_{i}\Theta(o_{i}^{c}-y_{i}^{c}), \tag{15}\]

Figure 3: True positive rate of each waveform in the test set against source distance.
Figure 4: Confusion matrices of the test set (left) and violin plots of network SNR of each waveform (right) from sources with distances of 1, 5, and 10 kpc.
where \(\Theta\) is the Heaviside step function. Unlike the previous metric, the higher this value is, the more effective the visualization method will be because there are more samples that score higher when given the explanation map than when given the original image.

For the three visualization methods Grad-CAM, Grad-CAM++, and Score-CAM, the two metrics described above were computed using signals from sources at 1 kpc in the test set. The results are summarized in Tab. 1. Score-CAM showed the best results in both metrics, meaning that it is the best visualization technique for our model among the three CAM methods considered in this study. We also qualitatively compared these methods by visualizing some samples. One example is shown in Fig. 6. The input image is represented by a color image, with the red, green, and blue channels corresponding to the H1, L1, and V1 spectrograms, respectively. All the three saliency maps take large values around the SASI mode around 100 Hz. At high frequencies, the Grad-CAM and Grad-CAM++ maps only take slightly larger values around 1 kHz, whereas Score-CAM has g-mode-like arch shapes around 1 kHz. It suggests that the visualization by Score-CAM captures more of the input features that were discriminative for the prediction.

As discussed above, we determined that the Score-CAM is the optimal method for generating saliency maps for our model. We produced saliency maps by Score-CAM for input of each class, which can be seen in Fig. 7. In the input images, as in the previous figure, the red, green, and blue channels corresponding to the H1, L1, and V1 data, respectively. This means, for example, that in the reddish image such as the m39 sample in Fig. 7, the SNR at the H1 detector is smaller than that at the L1 and V1 detectors. All of the plotted signal samples are scaled to have SNR of 40 and were correctly classified by our model. We plotted several CAM maps for noise samples in addition to the one in this figure, but the regions that the model sees to predict them to be noise were random. In the he3.5 and s13 samples, we can see that the model focuses on g-mode arch shape, especially their low and high frequencies area. In the s18 and y20 model, the model sees the entire g-modes. In the s18np model, the CAM map indicates that the model considers not only g-mode but also prompt convection and SASI. The s25 model has SASI activity, but its amplitude is not too large, and the CAM map shows that the model's prediction is based on the prompt convection and the high-frequency g-mode. In the m39, z85_sfhx, and z100_sfho model, the CAM maps take large values at high frequencies in g-mode. In addition, in z85_sfhx and z100_sfho, the low-frequency SASI mode, whose frequency increases with time, is also visible in the CAM maps. To summarize these outcomes, we found that the model looks at g-mode in all signal waveforms, and also looks at SASI and prompt convection in some signal waveforms when classified.

Additionally, we plotted saliency maps of the misclassified samples. Figure 8 shows the spectrogram of the s25 signal sample at each detector and the Score-CAM map, which the model classified as s18np. The SNR of this signal is 85, which is quite large, and the g-mode and the prompt convection are visible at H1 and L1 detectors, but there is a glitch in the strain at L1 detector. The Score-CAM map shows that the model focuses on the prompt convection and the glitch, which are used to determine that the signal is s18np. Another example is plotted in Fig. 9. This sample contains a s13 signal with an SNR of 48, and there are no glitches, but the model classified it as y20. We can see SASI induced GW mode around 100 Hz from 0.2 to 0.4 s especially in the L1 spectrogram, but the Score-CAM map indicates that the model sees g-mode and does not see the low frequency mode.

From the misclassified samples and the Score-CAM maps, it was found that the performance of the model is sometimes affected by glitches, and does not fully take advantage of the characteristics of the signals. The former could be resolved by generating the training sets containing more glitches, and the latte

\begin{table}
\begin{tabular}{l c c c} Method & Grad-CAM & Grad-CAM++ & Score-CAM \\ \hline Ave. Drop (\%) & & & \\ (_Lower is better_) & 30.70 & 17.58 & **9.61** \\ \hline Ave. Increase (\%) & & & \\ (_Higher is better_) & 1.30 & 1.40 & **1.96** \\ \end{tabular}
\end{table}
Table 1: Results for evaluation of the explanations generated by Grad-CAM, Grad-CAM++, and Score-CAM on the test set.

Figure 5: Features of the test samples at 1 kpc extracted by CNN and mapped into two-dimensional space by the t-SNE algorithm.
Figure 6: Qualitative comparison of three CAM maps for the s25 sample at 1kpc.

Figure 7: Input spectrograms and Score-CAM maps of correctly classified samples. The SNR of each signal sample is 40.
using better time-frequency representation able to reflect the various features that the CCSN signals have.

## IV Conclusions

In this study, we trained a two-dimensional CNN model to classify CCSN GW signals immersed in real noise of O3 observation data. Our model achieved a high accuracy of 98.4% for signals from sources with distances of 1 kpc, but the model struggled to correctly identify most of the signals from sources at 10 kpc.

To interpret our model, we used t-SNE algorithm and mapped the extracted features by the convolutional layers into a two-dimensional space. The dimension-reduced features show that the convolutional filters could extract meaningful features that are significant for classifying the signals.

To gain insights into the decision-making process of the model, we applied the CAM technique to visualize the regions in the inputs that were influential to the predictions. Three methods, Grad-CAM, Grad-CAM++, and Score-CAM were considered and we concluded that the Score-CAM is the best for our model in terms of the average drop and average increase metrics. The Score-CAM maps of correctly classified signal samples revealed that the model's predictions were heavily affected by a part of the entire g-mode in the spectrogram of each signal. In some waveform models such as s18np, s25, z85_sfho, and z100_sfho, their CAM maps suggest that the prompt convection or SASI induced GW mode also affects the model's prediction.

In this analysis, a time-frequency map was created from the short-time Fourier transform, but its resolution is limited by the uncertainty relationship between time and frequency. In future studies, we would like to improve the accuracy of the CNN model by using methods such as the Hilbert-Huang transform [52], which can generate higher resolution time-frequency maps, and to confirm that the CNN can also utilize several more GW modes to classify CCSN signals.

###### Acknowledgements.

 The authors would like to thank J. Powell for providing us gravitational wave simulation data. This research was supported in part by JSPS Grant-in-Aid for Scientific Research [No. 22H01228 (K. Somiya), and Nos. 19H01901, 23H01176 and 23H04520 (H. Takahashi)]. This research was also supported by the Joint Research Program of the Institute for Cosmic Ray Research, University of Tokyo and Tokyo City University Prioritized Studies. This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation, as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan.
