Neural Network Fusion Processing and Inverse Mapping to Combine Multisensor Satellite Data and Analyze the Prominent Features

In the last decade, the increase of active and passive earth observation satellites has provided us with more remote sensing data. This fact has led to enhanced interest in the fusion of different satellite data since some of the satellites have properties complementary to others. Fusion techniques can improve the estimation in areas of interest by using the complementary information and inferring unknown parameters. They also have the potential to provide high-resolution detailed classification maps. Thus, we propose a neural network, which combines and analyzes the data obtained from synthetic aperture radar (SAR) and optical sensors to provide high-resolution classification maps. The neural network employs a novel activation function to construct a neural network explainability method termed as inverse mapping for prominent feature analysis. By applying inverse mapping to the data fusion neural network, we can understand which input features are the prominent contributors for which classification outputs. Inverse mapping realizes backward signal flow based on teacher-signal backpropagation dynamics, which is consistent with its forward processing. It performs the contribution analysis of the data pixel by pixel and class by class. In this article, we focus on earthquake damage detection by dealing with SAR and optical sensor data of the 2018 Sulawesi earthquake in Indonesia. The fusion-based results show increased classification accuracy compared to the results of independent sensors. Moreover, we observe that inverse mapping shows reasonable explanations in a consistent manner. It also indicates the contributions of features different from straightforward counterparts, namely, pre- and post-seismic features, in the detection of particular classes.

spectral, and temporal resolutions [1] and possess different peculiarities designed for their specific tasks. Earthquake damage detection is a field where remote sensing satellites are extensively helpful. They can assess an inaccessible region after a disaster and can cover a wider scene, because of their synoptic imaging capability, in particular when the event is located in a remote area or a region where communication infrastructures have been damaged. In such cases, a damage map is required for the establishment of the first step in the evacuation and mitigation plan [2]. Remote sensing technologies have shown to be effective damage evaluation tools in such cases by providing high amount of critical information spatially and temporally. Conventionally, optical satellites have been used for the direct detection of earthquake damages [3], [4], [5]. Optical satellites with their ever-increasing spatial resolutions can provide us with easy-to-understand general surface information.
Synthetic aperture radar (SAR) sensors, on the other hand, have the potential to detect subtle targets, such as ships and oil spills [6], [7]. SAR sensors are also prominent data sources for rapid disaster assessment applications, as they have the ability to view the earth's surface irrespective of the time of the day. In addition, they have the ability to detect surface deformations with high levels of accuracy [8]. Conventional methods for damage assessment using SAR include generating interferometric coherence changes from multitemporal SAR images for wide region mapping [9], [10], building level mapping [8], [11], and rapid damage mapping [12], [13], [14], [15]. Recent studies have focused on the fusion of SAR and optical remote sensing for forest mapping [16], [17], [18], [19], land-use land-cover applications [20], [21], [22], and earthquake damage detection [23], [24]. In such cases of fusion, SAR can provide information that is not visible in the optical data and vice versa. The two data sources are, thus, regarded as complementary, not competing. Multisensor fusion studies on earthquake damage detection [23], [24] focus on the detection of damage-stricken urban areas as they tend to have a higher fatality rate.
A recent paper [25] reported multisource data fusion using ensemble learning classifiers along with feature importance analysis. The study gave a comprehensive comparison of various combinations of SAR, optical, and digital elevation model (DEM) datasets for damage classification using three types of ensemble classifiers. It was concluded that DEM-and SAR-derived features contribute the most to the overall damage classification. The classification focuses primarily on building damage mapping and classifies the buildings as destroyed, damaged, possibly damaged, or no damage. The classification requires multiple stages of processing. Thus, the importance analysis becomes complicated and indirect, making it harder to relate the observations with the analysis decision. In addition, in order to provide more information about the damaged areas, it is necessary to focus on various other land-type classes as well, which was not accomplished in this article.
This article proposes the use of SAR and optical sensor data to generate high-resolution land classification maps for the accurate damage estimation of various landforms, not just urban areas. We construct a fully connected neural network with a novel activation function for the fusion of different features of various sensors to classify each pixel as one of the predefined land types. In addition, an approach for neural network explainability is also proposed.
Neural networks are very useful and effective in solving problems difficult for traditional feature extractors and classifiers [26]. To build trust in such intelligent systems and integrate them with our society, there should be clear transparency in the neural models that can explain why they predict what they predict. In the area of visual interpretability of neural networks, the attribution value method and feature visualization have gained significant attention [27]. Important pixels related to a particular output are highlighted with the use of the attribution method. The feature visualization method, on the other hand, provides an understanding of a single kernel or convolutional layer [28], [29], [30], [31]. Moreover, convolutional neural networks (CNNs), which are widely used nowadays, pay attention to spatial features, such as textures and shapes. However, in the CNN, the resolution of the data degrades because of pooling.
The layerwise relevance propagation is also a wellestablished and effective method in the domain of neural network explainability [32]. It is a method that computes the relevance scores for the image pixels and image regions to perform a backpropagation of the score layer by layer. These scores typically denote the impact of a particular image region on the prediction of the classifier [33]. In this article, we propose a method termed as "inverse mapping" that analyzes the prominent features by the direct use of the inversely propagating signals. While the layerwise relevance propagation follows the principle of relevance conservation and proportional decomposition [34] for describing the network's decision making, the inverse mapping focuses on the conductivity of the neural connections for the signals by quantitatively assessing the weights of the forward propagating signals. In the experiments below, we find that, for the current application, inverse mapping obtains more reasonable results.
In this article, an investigation is performed to understand the pixel-by-pixel and class-by-class contributions of each feature in the classification results. Thus, we propose inverse mapping dynamics. Inverse mapping dynamics determines the prominent contributors among the input features for each pixel in the test image by paying key attention to the signal flow in the neural network [35], [36]. The analysis of the inverse mapping is also class by class, which can correlate the output classes to the relevant input features, wherein the importance of the individual pixels is aggregated to estimate the feature importance on a class level. We find that the inverse mapping shows reasonable and consistent results and often clearly indicates the contribution of noncounterpart features in the neural network's class-by-class decisions.
In the subsequent sections, we first focus on the proposed neural network and inverse mapping in Section II. Section III presents the experiments and results of the application of the proposed neural network for the detection of earthquake damage in Sulawesi and the determination of prominent features using inverse mapping. Section IV presents the relevant discussion. Finally, Section V concludes this article.

A. Neural Network Structure and Forward Processing
Fig . 1 shows the construction of the proposed fully connected neural network, which consists of an input terminal layer, a hidden neuron layer, and an output neuron layer. Consider and y = [y 1 y 2 . . . y j . . . y J ] T to be the input and hidden signals, respectively, to be the biases for the hidden and output layers, respectively, and 1 W and 2 W to be the input-hidden weights and the hidden-output weights obtained after a learning phase, respectively. The values of z, where z = [z 1 z 2 . . . z k . . . z K ] T , at the output layer of the neural network can be calculated as where f represents the proposed modified logarithmic activation function working componentwise and defined as Fig. 2 shows the curve of the proposed activation function. We employ error backpropagation learning for this activation function, with one-hot encoded teacher signals.

B. Inverse Mapping
Inverse mapping is a dynamics that traces which input is prominent for which output by paying attention to the signal flow in the network. It is based on teacher-signal backpropagation dynamics where the teacher signal propagates backward from the output layer toward the input layer [37], [38]. Though the dynamics shows learning results similar to those of usual error backpropagation in an ordinary (real-valued) neural network, they provide us with a direct suggestion on the inverse mapping as follows.  For understanding the inverse mapping, it is helpful to consider how the neural network determines a winner neuron, by focusing on the role of the weights of the forward processing network. In Fig. 1, the input signals propagate forward through 1 W and 2 W. The values of respective weights represent the "conductivity" of the neural connections for the signals, resulting in a set of outputs. Then, we can consider a reciprocal signal flow by taking transposed matrices 1 W T and 2 W T as the inverse signal propagation, and no additional learning is required for weight optimization for inverse mapping.   3 illustrates the inverse mapping process. The value obtained at the winning nodek in the forward processing is fed to the same node at the input of the inverse mapping, while all the other nodes are fed with zeros. This is done to suppress any influence that the nonprominent classes could introduce. The inverse mapping is represented by where z = 0 0 . . . zk . . . 0 T is the modified output values of the forward processing network fed as an input to the inverse mapping network, x is the output of the inverse mapping network, and the inverse activation function f −1 is defined as Thus, we employ the inverse of the function f −1 in the inverse mapping since the system deals with signal magnitude quantitatively [26]. The transpose weight matrices 1 W T and 2 W T realize signal propagation backward. For the validation of the modified logarithmic activation function in the proposed forward processing neural network, the modified logarithmic activation function was compared with the traditional activation functions. It was observed that, in comparison with the traditional activation functions, the proposed activation function realizes similar results for data classification. 1 However, in the case of inverse mapping, the expected range of inversely input values is (−∞, +∞). When an existing activation function, such as sigmoid, tanh, or rectified linear unit (ReLU), is used in forward processing, in inverse mapping, its inverse function (i.e., logit, inverse tangential hyperbolic, or inverse ReLU, respectively) is used. Logit, inverse tangential hyperbolic, or inverse ReLU have a limited inverse range, i.e., [0, 1], [−1, 1], or [0, +∞), respectively. Thus, the pixels having inverse input values beyond the defined range are not successfully mapped in the inverse mapping. On the other hand, the range of the proposed modified logarithmic function defined in (2) is (−∞, +∞), which makes it possible to deal with any large backward propagating value, which is ideal for inverse mapping.
The inversely obtained value x of the inverse mapping in (3) suggests the prominent factors to determine the class of the input pixel data. A large absolute value x i indicates the ith observation data contributing to the winner class determination for that pixel. For each pixel, all the elements of x are normalized by the maximum absolute value among the elements.

A. Study Area
Sulawesi Island, in eastern Indonesia, is at the triple junction of the Sunda, Philippine, and Australian plates [39]. On September 28, 2018, an earthquake of magnitude 7.5 struck Sulawesi. This earthquake was responsible for generating multiple landslides, soil liquefaction, and a tsunami.
The earthquake's epicenter was Donggala Regency (0.20 S, 119.89 E), 80 km north of the provincial capital city Palu. The most significant contributors to the casualties were the landslides within the alluvial valley to the south of Palu. Over 4340 people went missing from the villages, which got buried completely due to the landslides [2], [40].
In this article, we focus on three major landslide areas that were affected by the earthquake, namely, Sidera, Petobo, and Balaroa. The locations of these landslides can be seen in Fig. 4. These areas of interest (AOIs) have high geospatial diversity as Sidera AOI primarily consists of nonurban areas, while the Balaroa AOI consists of primarily urban areas and the Petobo AOI consists of an equal mix of both urban and nonurban areas. It was also observed that, despite being damaged by the same earthquake, the cause of damage was different in the case of Sidera/Petobo and Balaroa. Sidera and Petobo landslides originated as lateral spreads on the gentle slopes directly beneath the Gumbasa aqueduct. In contrast, the Balaroa landslide has been attributed to falling along the secondary fault line [40]. In addition, the analysis in this article is pixel by pixel, which also ensures that the data have high diversity.

B. Dataset
We map the damaged areas by using the L-band SAR sensor Advanced Land Observation Satellite-2 (ALOS-2) of the Japan Aerospace Exploration Agency and the Sentinel-2 optical sensor of the European Space Agency (ESA). The SAR datasets and optical datasets that were required for the damage assessment in this study have been described in Table I.

C. Preprocessing and Feature Extraction
Preprocessing operations are meant to correct the images for sensor-specific and platform-specific radiometric and geometric distortions and convert them to intelligible form. Preprocessing steps involve geometric and radiometric calibration of the raw data. In the case of SAR data, speckle filtering is additionally required. To avoid mismatch in the spatial resolution, resampling operation is performed on the Sentinel-2 data, following which all the bands of Sentinel-2 are sampled to 10-m resolution. The preprocessing of the satellite data is majorly performed by using ESA's Sentinel Application Platform (SNAP) [41]. From the preprocessed SAR and the optical datasets, we choose features that have been studied to show sensitivity to the damage caused by disasters. As shown in Fig. 5 from SAR data, we obtained the intensity (in dB) and the coherence values, where both these features are obtained in horizontal-horizontal (HH) and horizontal-vertical (HV) polarizations. From the optical data, we extract the Normalized Difference Vegetation Index (NDVI), the Normalized Difference Built-up Index (NDBI), and the Build-Up Index (BUI). These features are defined as follows.
1) Intensity refers to the backscattering received by the SAR sensor. Intensity values can provide useful information about the damage [42] as there is a change in the backscattering coefficient before and after the disaster. It is obtained as the square of the amplitude of the SAR raw image as Intensity = 20 log 10 (Amplitude). 2) Coherence is a feature obtained from SAR data that depicts the stability of an area. It is based on both the amplitude and phase of the backscattered signals. It shows where changes have occurred and detects even subtle changes with high levels of accuracy ranging from several millimeters to centimeters [24]. Coherence can be calculated as where P and S represent the complex amplitude obtained in primary and secondary observations, respectively, (·) * above them represents the complex conjugate transpose, and the bracket denotes the spatial average for 3 × 3 window. Coherence values range from 0 to 1. 3) The NDVI is an index obtained from the optical sensors.
The NDVI is often used to classify an area into vegetation or nonvegetation class [23]. This feature can be determined from the contribution of the satellite bands corresponding to visible red and near-infrared (NIR) wavelengths [43] corresponding to Band-4 (B4) and Band-8 (B8) of Sentinel-2, respectively, and is given by their power values as 4) The NDBI is an index that highlights the urban areas and is given by the normalization of the short-wave infrared (SWIR) and NIR wavelengths [44] corresponding to Band-11 (B11) and Band-4 (B4) of Sentinel-2 and can be calculated by their power values as 5) The BUI is the index used for the analysis of urban patterns using NDVI and NDBI. A higher value of a pixel in this index indicates a greater possibility that it represents a built-up area [45], which allows built-up areas to be mapped effectively. The BUI is given their power values as Tables II-IV illustrate the normalized features extracted from the multiple sensors. A total of 14 features were extracted from the given datasets. Seven features were extracted from the preseismic scenes, and seven features were extracted from the postseismic scenes.
To ensure geometric accuracy and good spatial alignment between the 14 extracted feature images, we perform collocation by using the SNAP toolbox. In the collocation operation, for each pixel in the primary product, the corresponding pixel is looked up (by geolocation) in the secondary product, and the primary and secondary images are mapped to the same geographical location. After that, the datasets are clipped very precisely to the coordinates of the three AOIs, namely, Sidera, Petobo, and Balaroa. All the feature values except for coherence HH and HV are then normalized in the range of −1 and 1 linearly.
When we observe intensity in Tables II-IV, we can see that in the areas where the landslide occurred, there was an increase in  II  FEATURES EXTRACTED FOR THE SIDERA AOI  TABLE III  FEATURES EXTRACTED FOR THE PETOBO AOI  TABLE IV  FEATURES EXTRACTED FOR THE  the brightness of the intensity in both HH and HV polarizations. However, in the coherence, it can be observed that after the landslide, there is a loss of coherence at the places where the landslide occurred. A similar change in the characteristic can be seen for the three optical features NDVI, NDBI, and BUI. In the case of NDVI, there is a loss in the NDVI due to the fact that the area had a loss of vegetation. However, for the NDBI and the BUI, the forest and urban areas were converted to barren lands, and thus, there is an increase in these indices for the affected regions.

D. Neural Network Training
To understand the effect of the fusion of SAR and optical features on the neural decision accuracy, we perform three different types of neural network processing, where the inputs are the following: 1) only SAR features; 2) only optical features; 3) both SAR and optical features. The input layer is fed with pixel-by-pixel values of the feature images shown in Tables II-IV obtained after the preprocessing steps illustrated in Section III-C. The weights of the neural network are randomly initialized and are updated by the backpropagation process.
At the output layer of the neural network, the network classifies each pixel as one of the following seven landform classes: damaged forest, damaged irrigation field, damaged urban, undamaged forest, undamaged irrigation field, undamaged urban, and undamaged dryland regions. Since the dryland region was not affected by the disaster, "damaged dryland class" was not included in our assessment. Table V illustrates the use of the AOIs in this study. We adopt a pixel-by-pixel image analysis method. Fig. 6 shows the ground truth of the Sidera AOI along with the training and validation pixels. We train the neural networks on randomly selected 3176 pixels (enclosed in black boxes) and validate them on 700 pixels (enclosed in blue boxes) of the ground-truth data of Sidera AOI. 2 Table VI shows the neural parameters that were determined to ensure that the number of free variables was similar in the three neural networks. The neural networks trained with the training data of the Sidera region are also applied to other unseen areas in Sidera as well as Petobo and Balaroa AOIs.

E. Neural Network Training Results
After the training on the selected pixels, the neural networks are applied to the unseen pixels of the various AOIs. They are expected to estimate which of the seven classes each pixel belongs to. The neural networks generate classification maps for the Sidera and the Petobo AOIs, which are compared with the ground-truth maps shown in Fig. 8(d) and (h), respectively. These ground-truth maps were derived from the land classification published in [40] and modified for the damaged region based on the damage maps provided by the Copernicus emergency management service following the disaster [46]. Moreover, since there is no ground truth for the Balaroa region, the results are compared visually with the RGB image given in Fig. 8(l). Fig. 8(a)-(c) are learning curves showing the training and validation accuracy for the neural networks trained on SARonly features, optical-only features, and both SAR and optical features, respectively. All of them illustrate smooth learning. In Fig. 7, we can observe an increase in the accuracy for the neural network trained on both SAR and optical features. Fig. 8(e), (i), and (m) shows the classification results for Sidera, Petobo, and Balaroa regions, respectively, for the neural network trained on SAR-only parameters. The classification is poor, and not much information can be obtained from this damage map. However, it can be seen that some of the undamaged dryland pixels and the urban pixels have been detected. Fig. 8(f), (j), and (n) illustrates the results for the neural network trained on only optical features. The classification is much better than the SAR-only classification, particularly because the shape of the landslide is detected well. However, with a close visual inspection, it is observed that some classes, specifically the urban regions, are not detected at all. Fig. 8(g), (k), and (o) shows the classification for the neural network trained on both SAR and optical features. The classification appears similar to the optical analysis results. However, unlike the optical-only classification, the urban class has been detected, and the identification of the other classes is better as well. Thus, it can be stated that the results of this neural network are much closer to the ground-truth data as compared to the classifications by the other neural networks.

F. Confusion Matrix
The classification maps in Fig. 8(g), (k), and (o) generated by the fusion-based neural network are compared for Sidera and Petobo AOIs with the ground-truth data shown in Fig. 8(d) and (h). Fig. 9 presents confusion matrices that compare the classification maps pixel by pixel with the ground-truth data. It shows what percentage of the AOI pixels is correctly identified in the resulting classification. A two-class comparison has also been performed, in which we group all damaged pixels as a single damaged class and all undamaged pixels as a single undamaged class.
The confusion matrices for Sidera region are shown in Fig. 9(a) and (b) for the seven-class and two-class classifications, respectively. Similarly, the confusion matrices for the Petobo region are shown in Fig. 9(c) and (d). The classifications have also been quantitatively evaluated by using precision scores (PSs). Each individual PS has been calculated as
In Fig. 9, the PS of each cell is shown in red. The confusion matrices indicate a high PS for most classes. However, the damaged urban class in Sidera AOI and additionally damaged irrigation field class in Petobo AOI show low PSs. Fig. 10 illustrates row-normalized confusion matrices. Here, each cell (i, j) represents what percentage of prediction of class i pixels was predicted as class j. In Figs. 9 and 10, it can be seen that though most classes have good scores, some classes have poor scores. The cause of the low PSs for these classes is attributed to the mismatch in the spatial resolutions. That is, the ground-truth data have a higher spatial resolution compared to the datasets used in this study. Thus, these precision issues exist. Confusion matrices were generated for the results from the only-SAR and only-optically trained neural networks as well for evaluation using an overall accuracy score. Fig. 11 presents the comparison of the accuracy score for the confusion matrices of all three neural networks. The seven-class results illustrate that the classification accuracy increases with the fusion of the data. From the two-class results, we can see that the classification accuracy is similar for optical-only and SAR and optical classifications. Thus, it can be postulated that, for optical-only classification, though the pixels were correctly identified as damaged or undamaged, the land type could not be accurately identified as compared to the fusion-based classification. This means that optical data detect less information compared to fusion-based data.

G. Inverse Mapping Result Using Both SAR and Optical Features
Section II-B explained the inverse mapping process. Here, in Figs. 12-15, we illustrate the inverse mapping results. In this experiment, if occasionally the value of the forward winner output zk in (1) is negative, indicating that the network is not confident of its estimation, the pixel is discarded in inverse mapping. The inverse mapping output x is plotted in box-whisker plots for all the pixels. In individual boxplots, the density of the pixels closer to −1 or 1 indicates the prominence of the corresponding features. The results of the inverse mapping illustrate the impact of the 14 features on the classification results. Table VII shows the inference based on the box and whisker plots in Figs. 12-15. Figs. 12 and 14 show the inverse mapping results for the pixels evaluated in the damaged category for Sidera and Petobo AOIs, respectively. In both the figures, we can observe that the damaged forest class was determined by post-seismic NDVI, pre-seismic NDBI, and pre-seismic NDVI. The damaged irrigation field class was selected by post-seismic NDVI, post-seismic BUI, post-seismic NDBI, and also pre-seismic coherence HV. Finally, the damaged urban class was chosen by pre-seismic BUI, pre-seismic coherence HV, and post-seismic NDBI.
Figs. 13 and 15 illustrate the inverse mapping results for the undamaged pixels. Post-seismic NDBI is considered a prominent feature for undamaged forest decision. Pixels labeled as undamaged irrigation fields are attributed to the post-seismic NDBI and pre-seismic BUI. Post-seismic coherence HH and pre-seismic NDVI influenced the undamaged urban class decision. Post-seismic NDVI and post-seismic BUI are prominent features for the undamaged dryland class. We can observe consistency in the results for the different AOIs. The results are reasonable for each class when we consider the land type and physical landslide dynamics. Furthermore, it can also be noted that features that independently contribute to the decision of a class worked collaboratively for the class, as observed for the damaged urban class.
Since the analysis including nonlinear effect may depend on the signal magnitude, additional experiments were conducted by replacing W T in (3) with W T |W| or with W T |W| 2 . However, the inverse mapping of both the experiments resulted in observations similar to that of the original inverse mapping, which confirmed the robustness of our method.

IV. DISCUSSION
This article focuses on the fusion of features from two satellites by using a fully connected neural network, which employs a modified logarithmic function. This study indicates that the accuracy of land-cover estimation increases with the fusion of data from the SAR and optical sensors.
For the classification using only SAR parameters, it can be inferred that the selected parameters were insufficient, and thus, the SAR-only combinations were not that impactful. Though SAR features did not show significant results when used alone, they increased the accuracy and interpretability of the damage classification results when used synergetically with the optical data. Certain classes, such as urban areas that were not detected by the optical sensor, were detected by the SAR sensor. This suggests that the coherence features of SAR carry information that complements the optical data. In the present study, it was also found that the impact of intensity features was not as prominent as that of coherence features. This could be because for building-detection tasks, coherence is a stronger feature compared to intensity.
We conducted additional experiments using the conventional layerwise relevance propagation method. The results showed deviation from those of inverse mapping. 3 Fig. 19(a)  that in the case of the damaged forest class, the layerwise relevance propagation shows the contribution of three pre-seismic optical parameters for the damage mapping. However, this is not reasonable because post-seismic features must be present for the interpretation of damaged classes. By comparing the inverse mapping and layerwise relevance propagation results for the damaged classes, we find that the results of inverse mapping are more reasonable. In inverse mapping, the combination of modified logarithmic activation function and W T deals with the nonlinearity effectively, providing more quantitatively convincing results.
The pixel-by-pixel approach introduced in this article can be extended to other multisensor data fusion tasks, such as change detection [47] and classification [48]. In addition, the use of inverse mapping can assist in a more application-specific choice of features, bands, or sensors.
In this experiment, we aim to realize high-resolution classification by paying attention to the physical reflectance values of individual pixels; thus, approaches related to texture or shape identification have not been considered. The pixel-by-pixel analysis using a fully connected neural network introduced does not hamper the spatial resolution. CNNs are not helpful because the use of CNNs degrades the resolution of the data by pooling. Unlike CNNs, which identify the texture or the shape from images, the fully connected neural network introduced here deals with the physical values directly.
We also focus on determining the prominent features for each class using inverse mapping. Consistent results were obtained for the three AOIs, and the results were reasonable for the particular classes. It was also found that features that individually contributed to the detection of any class worked together to detect that class. For example, in the case of damaged urban class, NDBI, BUI, and coherence HV were determined to be prominent, and these features had been previously studied to be individually good indicators of detecting changes in urban regions [24], [44], [45]. Traditionally, for change detection, counterparts such as pre-seismic and post-seismic data of a single feature are used [8], [13]. In our study, it was found that such counterparts do not necessarily complement each other when used with other features. Investigation of this observation can be the future scope of this work.

V. CONCLUSION
In this article, we proposed an explainable neural network with a novel modified logarithmic activation function that combines the capabilities of various damage-indicating features obtained from SAR and optical sensors. The combination of two sensors was found effective for land-type and damage detection and gave higher accuracy decisions. In this study, we introduced inverse mapping, using which we were able to identify consistent and reasonable prominent features for damage assessment. It was observed that features that independently contribute to damage detection of a particular class also determine collaboratively the specific class when used together. The prominent features belong to both SAR and optical sensors, highlighting the necessity of sensor integration. This study illustrated that the synergetic use of multiple sensors can provide us with significant information to generate effective postdisaster maps in an explainable manner.

APPENDIX
Here, we present Fig. 16, which shows the forward processing classification results for modified logarithmic function, tanh, sigmoid, and ReLU activation functions. The classification performances of all the activation functions are similar to one another. In addition, in Figs. 19 and 20, we show the results of layerwise relevance propagation for damaged and undamaged classes, respectively. The results for the layerwise relevance propagation show that a set of only pre-seismic features is responsible for the damaged forest class, which is not reasonable    for change detection. In contrast, the inverse mapping gives more reasonable results for each class.