Assessing Earthquake-Induced Urban Rubble by Means of Multiplatform Remotely Sensed Data

Earthquake-induced rubble in urbanized areas must be mapped and characterized. Location, volume, weight and constituents are key information in order to support emergency activities and optimize rubble management. A procedure to work out the geometric characteristics of the rubble heaps has already been reported in a previous work, whereas here an original methodology for retrieving the rubble’s constituents by means of active and passive remote sensing techniques, based on airborne (LiDAR and RGB aero-photogrammetric) and satellite (WorldView-3) Very High Resolution (VHR) sensors, is presented. Due to the high spectral heterogeneity of seismic rubble, Spectral Mixture Analysis, through the Sequential Maximum Angle Convex Cone algorithm, was adopted to derive the linear mixed model distribution of remotely sensed spectral responses of pure materials (endmembers). These endmembers were then mapped on the hyperspectral signatures of various materials acquired on site, testing different machine learning classifiers in order to assess their relative abundances. The best results were provided by the C-Support Vector Machine, which allowed us to work out the characterization of the main rubble constituents with an accuracy up to 88.8% for less mixed pixels and the Random Forest, which was the only one able to detect the likely presence of asbestos.


Introduction
The enormous amount of rubble from partial or total collapse of buildings/structures caused by an earthquake hitting urbanized areas with vulnerable historical centres, like those of many Italian towns, must be rapidly mapped and characterized after catastrophic events. This activity is mandatory to support the emergency activities and ensure the accessibility, rescue and first assistance to population. During the post-emergency phase, reliable information about the distribution and characterization of seismic rubble in term of volume/weight and typology is fundamental for their proper management which involves handling, removal, pile up and transportation to recycling facilities or final disposal sites [1]. Despite this information being particularly valuable for optimizing the post emergency responses, there are not many methods to provide extensive and reliable estimates of the amount of seismic urban rubble in terms of volume/weight and typology distribution, including the detection of possibly dangerous materials. With this aim, remote sensing techniques, coupled with GIS, can provide valuable opportunities [2][3][4][5][6].
through the acquired VHR multispectral data was also tested [30,31]. To better discriminate the high heterogeneity of seismic rubble, Spectral Mixture Analysis (SMA) was preferred to commonly used hard classification schemes, which are more effective in cases of the prevalence of pixels with homogeneous covers, which can be univocally labelled according to their spectral signature. The SMA is based on the so-called soft classification approach, where each pixel is characterized by the presence of different pure materials (endmembers), each of which has its own spectral signature and contribute to the pixel reflectance spectrum proportionally to its relative abundance. SMA has been widely applied to overcome the mixed pixel problem, a typical issue associated with medium and coarse-resolution remote sensing imagery also in urban application [32][33][34][35][36]. The spectral endmember signatures can be obtained from the imagery or laboratory/field measurements. Generally, however, due to different spectral configuration, atmospheric noises and difference in acquisition conditions, these latter are not in agreement with the spectral response of airborne or satellite digital imagery. Therefore, it can be advantageous to extract endmembers directly from imagery. Usually, the derivation of spectral endmembers from multi/hyperspectral imagery is typically achieved through the implicit Pixel Purity Index (PPI) or explicit use of convex geometry in multidimensional spectral space where pixels reflectance spread out. In general, PPI is used with hyperspectral imagery after their Minimum Noise Fraction (MNF) transform to minimize noise in selected components, whereas convex geometry methods, like Sequential Maximum Angular Convex Cone (SMACC), focus on extreme pixels in spectral channel multidimensional distributions as good candidates for endmembers [32][33][34][35]. In addition to endmembers, the various SMA approaches also make it possible to obtain their abundances within the pixel, using its mixed reflectance responses through different algorithms, including SMACC [19,37].
Even if urban and suburban environments, characterized by complex landscapes and elements with different spectral signatures, may reduce the performance of traditional SMA approaches with a fixed set of endmembers, the improvement of radiometry, spatial and spectral resolution of the new generation of VHR sensors make it possible to obtain sub-meter performances [37,38]. In fact, although the SMA approach has been widely used for urban EO applications using the hyperspectral reflectance data, it also provides useful results by exploiting multispectral VHR satellite information coupled with advanced machine learning algorithms within classification and regression schema [24,33,38,39]. In this study, the SMA approach was applied to WorldView-3 (WV3) pansharpened data, after their geometric and atmospheric pre-processing, to provide the related spectral endmember sets and abundances through the SMACC algorithm [40].
The dramatically augmented availability of the VHR multi/hyperspectral and Synthetic Aperture Radar (SAR) data provided by increasing number of operative satellite remote sensing missions must be suitably coupled with advanced approaches based on data mining, machine learning and clustering scheme for properly monitoring and characterizing increasingly heterogeneous, complex and wide urban areas. In this context, the most recent machine learning algorithms based on Artificial Neural Network (ANN), Support Vector Machine (SVM), K-Nearest Neighbour (KNN) and Random Forest (RF) play a relevant role [41]. In our implemented methodology, once the in situ acquired spectral signatures had been resampled to the same sensor configuration, these algorithms were tested within a supervised learning scheme, for mapping (classifying) the endmembers provided by SMA into the rubble material of interest detected on site.

Study Area
Central Italy is a rural and mountain landscape with several small ancient villages characterized by poorly engineered buildings (older residential units) and modern structures (newer buildings and life-lines made with anti-seismic criteria). The Apennines are a NW-SE oriented mountain belt affected by a multiphased contractioned and extensional tectonics [1]. It is one of the most seismically active areas in Italy and Europe, as African European continental margin pile-up during the Neogene formed a fault and thrust system belt. Earthquakes and aftershocks periodically hit Central Italy (where Lazio, Umbria, Marche and Abruzzo Regions are located). ISPRS Int. J. Geo-Inf. 2020, 9, 262 4 of 26 The study area was struck by a Mw 6.0 mainshock on 24 August 2016. About 300 people died and 20,000 homeless were left on a wide area. The seismic sequence was characterized by a subsequent Mw 5.9 mainshock on 26 October, followed by a Mw 6.5 near the town of Norcia on the 30 October. The last mainshock, with a maximum Mw of 5.5, occurred in the southern part on 18 January 2017, but it triggered an avalanche that wiped out the Rigopiano Hotel, causing another 29 victims.
The central Apennine is periodically affected by normal-faulting sequences. The area was already severely damaged by the great earthquakes of Naples and Aquila of January-February 1703 and previously the source of a strong earthquake in October 1639, parameterized in CPTI15 (Parametric Catalogue of Italian Earthquakes) with a Mw 6.2, and with a less extensive severe damage distribution and lesser intensity achieved compared to this of 2016. Struck at the time were some fractions of Amatrice located to the west of it (Figure 1b), while the earthquake was felt in L'Aquila, Ascoli Piceno and Rieti, but not in Rome. Overall, the earthquake of 2016 could be a twin event to that of 1639, generated-with the due differences in rupture length and directivity-from the same seismogenic source related to extensive tectonics. Therefore, such events are frequently followed by a long-lasting aftershock with many low-energy tremors. Between 2016 and 2017, more than 70,000 earthquakes were recorded in the study area [1]. The central Apennine is periodically affected by normal-faulting sequences. The area was already severely damaged by the great earthquakes of Naples and Aquila of January-February 1703 and previously the source of a strong earthquake in October 1639, parameterized in CPTI15 (Parametric Catalogue of Italian Earthquakes) with a Mw 6.2, and with a less extensive severe damage distribution and lesser intensity achieved compared to this of 2016. Struck at the time were some fractions of Amatrice located to the west of it (Figure 1b), while the earthquake was felt in L'Aquila, Ascoli Piceno and Rieti, but not in Rome. Overall, the earthquake of 2016 could be a twin event to that of 1639, generated-with the due differences in rupture length and directivity-from the same seismogenic source related to extensive tectonics. Therefore, such events are frequently followed by a long-lasting aftershock with many low-energy tremors. Between 2016 and 2017, more than 70,000 earthquakes were recorded in the study area [1]. The consequences of such a long series of seismic events is that building collapse can occur several days after being seriously damaged by the most intense earthquakes.  The consequences of such a long series of seismic events is that building collapse can occur several days after being seriously damaged by the most intense earthquakes.

Dataset
The following territorial data and information were used:

Data Processing
The whole methodology was developed according to the general scheme depicted in Figure 2 and described in the following paragraphs. The left side of Figure 2 (dashed box), regarding the procedures for retrieving the rubble heaps boundaries and the relative volumes, was actually the objective of a previous work [1]; hence, its description is reported only in the Appendix A, whose objective is to add crucial information regarding the main (surface) constituents of rubble heaps to this geometric information.

Determination of Heap Volume
The rubbles heaps boundaries were basically provided by the Copernicus EMS mapping, furtherly improved (according to our purpose) by visual interpretation of a 15 cm RGB orthophotos photointerpretation. For the volume estimation, Digital Elevation Model of terrain (DTM) and Digital Surface Model (DSM), produced from LiDAR surveys, were exploited: in particular pre-event DTM (DTMpre) and post-event DTM and DSM (DTMpost and DSMpost) were considered.
To each polygon (i.e., heap boundary) the building's base elevation value was manually assigned, deriving it from the DTMpost in correspondence with a visible point (e.g., portion of the surrounding streets free of rubbles), identified close to the considered rubble. The volume of a single rubble heap was estimated as the solid (3-D space) included between the bottom surface represented by the building's foundation developing at a constant altitude and the upper uneven surface represented by the DSMpost values over the polygonal area itself. In other words, the volume is assessed through the difference between the DSMpost and the DTMpost values (constant altitude) over the area identified by the heap's perimeter. An alternative approach has made it possible to perform the assignment of the elevation values of rubble's base automatically. The base values are automatically extracted by the pre-event DTMpre, which can be assumed as the bottom surface. Therefore, the volume calculation is directly obtained by the difference in elevation between post-event DSMpost and pre-event (DTMpre over the rubbles' perimeter). For further information about volume heap calculation of Cappucci et al. [1], see Appendix A. (2) WorldView-3 multispectral data; (3) hyperspectral signature library; (4) polygonal cover and related attributes of rubble heaps, from previously developed procedure [1], in input to SMACC and GIS processing. The rough outline of rubble volume estimates procedure (dashed box) has been included here only for completeness.

Determination of Heap Volume
The rubbles heaps boundaries were basically provided by the Copernicus EMS mapping, furtherly improved (according to our purpose) by visual interpretation of a 15 cm RGB orthophotos photointerpretation. For the volume estimation, Digital Elevation Model of terrain (DTM) and Digital Surface Model (DSM), produced from LiDAR surveys, were exploited: in particular pre-event DTM (DTMpre) and post-event DTM and DSM (DTMpost and DSMpost) were considered.
To each polygon (i.e., heap boundary) the building's base elevation value was manually assigned, deriving it from the DTMpost in correspondence with a visible point (e.g., portion of the (3) hyperspectral signature library; (4) polygonal cover and related attributes of rubble heaps, from previously developed procedure [1], in input to SMACC and GIS processing. The rough outline of rubble volume estimates procedure (dashed box) has been included here only for completeness.

In Situ Data Collection and Pre-Processing
The hyperspectral signatures of rubble materials were collected during a field survey carried out on the 6th and 7th of December 2016. Up to 150 hyperspectral signatures (each made up of about 2300 spectral narrow bands) were acquired. Most of them referred to the five principal building materials (rubble debris, concrete, tiles, bricks, natural stones) and the others to secondary/accessory materials (i.e., zinc/copper gutter, bitumen sheath). It is worth noting that the spectral signatures of the five principal building materials differed significantly among the various piles/urban places. Consequently, after the usual pre-processing (interpolation in the water vapor absorption bands, filtering and normalization in case of calibration faults), for each of the main five building materials, four different averaged signatures were worked out and stored in the spectral library. Therefore, the spectral library resulting from the in situ measurements included eight macro-classes (the abovementioned five main building materials plus the three accessory materials), characterized by twenty-three averaged signatures, with twenty (5 × 4) related to the main five building constituents and three to the accessory materials).
Additionally, two hyperspectral signatures of differently aged and colonized (by vegetal species) asbestos building materials were retrieved from the literature and included in the library, in order to test the possibility of also detecting this dangerous material through satellite data [42]. Therefore, in the end 25 hyperspectral signatures distributed into 9 macro-classes (considering one asbestos macro-class subdivided into the two different aging stages) were finally available.
To make these hyperspectral signatures compatible with the spectral configuration of the WV3 sensor, they were resampled using the related band spectral filters functions. The satellite sensors use electronic band-pass filters to achieve the different spectral band responses from the broad spectrum captured signal. Their analytic expression filter functions, indicated as f ij , (i, j indicating, respectively, wavelength and band) and ranging between 0 and 1, can be used to derive the multispectral signatures from the hyperspectral ones by convolving them as follows: where b j stands for multispectral response in band j, f ij is the filter function of j band of satellite sensor (available from the data provider) and R i the hyperspectral reflectance data. In this context, this resampling step is mandatory to make the in situ hyperspectral signatures compatible with the multispectral ones derivable from the WV3 data.

Satellite Data
Acquisition and pre-processing (geometric, radiometric and atmospheric corrections) of VHR satellite remotely sensed and auxiliary data.
The geometric pre-processing and the orthocorrection of WV3 frames were accomplished using suitable DEM and the rational polynomial coefficients (RPC) provided with the WV3 files with the support of detailed cartography in the WGS 84 / UTM zone 33N projection [43].
In general, satellite imagery is affected by light-wave scattering and attenuation from air molecules, haze, water vapor and particulates, depending on atmospheric turbidity. These noise factors may be present in any given scene and they are typically not uniformly distributed. Given the thickness of the atmosphere crossed by the satellite signals coming from earth surface, the degrading effects introduced may be relevant in many cases, even in some cases making remotely sensed imagery unusable. In our case, aiming at exploiting the field hyperspectral reflectance of specific building materials, the removal of atmospheric effects from satellite imagery is necessary in order to retrieve the corresponding reflectance at ground level. To this end we used the FLAASH (Fast Line-of-sight Atmospheric Analysis of Hypercubes) code, integrated in L3HARRIS Geospatial ENVI ® image analysis software, which allows the automatic retrieval of aerosol optical thickness on processed frames using well established algorithms [44,45]. It is worth noting that among the two options we tested-atmospherically correcting the WV3 images first and then performing pansharpening versus first performing the pansharpening and then the atmospheric correction-the latter provided the better results when checked visually.
Moreover, in order to improve the geometric accuracy of the final results, the panchromatic data provided by WV3 sensor at 0.33 m a.g.r. were exploited using the Gram-Schmidt pansharpening method. This choice represented a compromise to preserve both spectral and spatial information, so as to improve the geometric accuracy of final results [46]. Although its performances for spectral features are not at top, it was selected on the basis of preliminary visual interpretation of the results obtained using atmospherically corrected data [47]. In addition, considering that our interest was in improving spatial resolution and that the quality assessment of pansharpened images is still an open issue, the selection was determined also by easy to use and robust implementation of the algorithm within the most diffuse packages for satellite remotely sensed data processing [46,47].
In general, the pansharpening and atmospheric pre-processing of satellite data are performed to improve the results of subsequent classification (based on categorical labelling approach) or regression (based on continuous variables modelling) processes. Few authors have dealt with the evaluation of the suitability of the sequence in performing these two pre-processing steps (commutability) for these two different needs [48]. However, in agreement with the approach adopted here, it was reported that performing the atmospheric correction after the pansharpening, provides better results in continuous modelling approaches, [49].

SMA-SMACC
Data provided by the abovementioned multispectral synthetic image were used as input to the SMA procedure. Indeed, the composite nature of remotely sensed spectral information often hampers the detailed identification and mapping of targeted constituents of the earth's surface. SMA is a well-established and effective technique to address this spectral mixture problem. SMA models the pixel mixed spectrum as a linear or nonlinear combination of its constituent spectral components (or spectral endmembers) weighted by their subpixel fractional cover. By model inversion SMA provides subpixel endmember fractions [50]. SMA models can be classified into linear and nonlinear ones, depending on the way they model the complexity of light scattering in the spectral mixture analysis. In this paper, we consider linear spectral mixtures, which is considered more effective in the case of an absence of interaction between the component responses, and which may be expressed as: where: Y(i,k) = ith pixel spectral response in kth spectral band; X j (i,k) = ith pixel spectral response of jth endmember in kth spectral band; f j (i) = pixel fractional abundance of jth endmember; v(i,k) = pixel residual noise in the kth band; with n j f j (i) = 1 as basic constraint to fractional pixel reflectance contributions. In the Equation (2), X j (i,k) indicates the spectral signature of the j th endmember, Y(i,k) is pixel reflectance response, found through SMACC algorithm by applying consecutively the maximum convex cone constraint. The algorithm also assesses the related fractional abundance f j (i) by minimizing the residual term v(i,k). In the case of an a priori fixed number of endmembers, usually the n j=1 f (i) j < 1 constraint is assumed to improve the effectiveness of the algorithms by introducing also a shadow endmember. This latter includes the remaining of reflectance contribution to unity, arising for example from possible different illumination condition of the endmember portions within the pixel. Assuming the shadow endmember includes proportionally all the less illuminated parts of all others, estimated through SMA, their preventive normalization at pixel level is required, using their sum as normalization factor.
To retrieve the endmembers, the SMACC algorithm selects the extreme points at borders of the multidimensional distribution of pixels' spectral signatures establishing a convex cone delimited by the related extreme pixel spectral vectors (endmember spectral signatures), which are located and used as spectral endmembers. The algorithm starts with a single endmember with maximum albedo and incrementally increases their number using the oblique projection that maximize angle between the previous set of pixel extreme vectors to find the subsequent endmember. The cone grows by encompassing the new endmembers, and the iterative procedure repeats until an endmember is found that is already contained within the convex cone, or until the user-defined maximum number of endmembers has been reached. The convex cone formed by extreme vectors or endmembers contains the remaining data vectors that can be modelled by a linear combination of endmembers with the positive coefficients being identified as abundances of the endmembers. It is worth noting that the number of endmembers can ideally exceed that of used spectral channels of the sensor. This is a very useful feature especially when we deal with multispectral data having a limited number of reflectance bands. However, their number should be maintained at a similar level, even when a limited endmembers basis may represent/guarantee just a coarse approximation of the many materials and their spectral variations found in seismic rubble, in order to ensure that a well-conditioned and linearly independent model is obtained. In fact, endmember extraction from hyper/multispectral data remains a challenge task, considering the within-class spectral reflectance variance and noise contributions from various factors (sensor, atmosphere, illumination) to the endmember signatures. Sometimes an endmember may characterize a mixture of materials having very similar spectral signatures rather than represents a single material. More endmembers than pure covers are often retrieved in a remotely sensed imagery and frequently some of them, less important for scene description, can be coalesced in the above introduced shadow endmember.
The SMA through SMACC was performed on the data from the 8-band multilayer at 0.3 m. of a.g.r., obtained from WV3 pansharpening and atmospheric correction processes, considering only the pixels belonging to the identified rubble heaps. The area investigated through SMA corresponds with pixels belonging to the identified rubble heaps furtherly reduced by excluding those pixels whose vegetation cover was too relevant (e.g., rubbles spread on tree-covered scarps). This was achieved by selecting only pixels whose Normalised Difference Vegetation Index (NDVI) values were below a suitable threshold.
Taking into account the number of available bands and the number of the materials to be discriminated, the number of endmembers to be worked out by SMA was fixed at 10 (plus shadow), being aware that a small excess may be tolerated by the method without compromising its robustness too much. In addition, the SMA algorithms provided the percentage of these endmembers for each pixel within the heap surface.

Machine Learning
Then we adopted a supervised classification schema based on various machine learning algorithms, using the 25 resampled spectral signatures of the materials included in the (spectral) library as a training set in order to link these spectral signatures to the endmembers retrieved by SMACC. It is worth underlining that the association of the retrieved endmembers with the signatures of rubble materials acquired on field was rather complex, due to both the above-cited noise factors affecting remotely sensed techniques and the intrinsic spectral variability of field data, characterized indeed by multiple signatures for the same building material. Six typical machine learning algorithms were preliminarily selected for supervised classification of endmembers, taking into account their different capabilities and performance in various situation of noise, limited class samples and outlier presence in input data. The selected algorithms were: C4.5 (decision tree), C-SVC (Core Support Vector Machine-discrete class, continuous input), BVM (Ball Vector Machine), KNN (K-nearest neighbours), RF (Random Forest), ANN (Artificial Neural Network). The latter, though, preliminarily provided very poor results, so it was not considered anymore in the processing chain. C4.5 makes it possible to implement easy-to-interpret classification models by hierarchically splitting the training data set. This algorithm selects the best subset of attributes based on an entropy measure and organizes the classes in a decision tree rule-based structure. Each node of the tree relates to a split in the feature space, which is always orthogonal to its axes [50][51][52]. C-SVC, more often used in the per-pixel classification context, is a non-parametric supervised statistical learning technique with robustness against outliers and limited train. It is able to estimate a hyperplane in the feature space that minimizes misclassifications [41,52]. BVM is an improvement of standard SVM algorithm which combines techniques from computational geometry based on hyper-spheres to efficiently achieve a close-to-optimal solution [53]. The K-nearest neighbour algorithm (KNN) is a method for classifying objects based on the closest training examples in an n-dimensional feature space. Given an unknown feature pattern, the classifier searches the pattern space for the k training tuples that are closest to the unknown one [53]. The ANN nonparametric algorithms are based on the neural network concepts and work without assumptions about input data distribution and independency. They learn from the training dataset and build relationships (networks) between input (features) and output nodes (classes) through hidden neurons layer connection weights modulation. In our context, a critical issue for ANN effectiveness may be the amount of training occurrences [41]. The Random Forest classifier consists of a group of decision trees induced with different sub-sets of the training data. Each tree of the forest casts a vote for the class to which a given analysis unit (in this case, a given segment) should be associated. The class with most votes is the one associated with the segment [41,54].
These machine learning algorithms were tested as classifiers of endmembers with the objective of evaluating their capability to discriminate as much as possible among the 25 materials present in the library including asbestos. The selection was performed on the basis of error rate minimization in recognizing the classes of training set coupled with the assignment of the maximum number of main materials to endmembers.

Accuracy Assessment
The accuracy assessment of the thematic maps obtained from the usual hard classification scheme is a mandatory step to be accomplished using randomly selected samples, where the categorical class assignment provided by the classifier algorithm is verified through independent observations. However, this approach was not directly applicable to the results of SMACC SMA analysis, which provide continuous distributions of endmembers abundances within each pixel. Therefore, an original approach that also allowed us to assess accuracy for SMA-derived endmember abundance distributions was implemented. To this end, only the pixels where one endmember predominates the others and its abundance is superior to a predefined threshold were selected to be considered quasi pure in a spectral sense. To this end the pixel Predominant Endmember class Distribution (PED) was obtained by the maximum abundance distribution (PMA), in turn derived from the endmember distributions using different thresholds to exclude the pixels with lower PMA values. This allowed us to transform the multilayer of endmember abundances continuous distributions into a categorical distribution of predominant endmember classes (PED) to which a classical accuracy assessment scheme was applied.
Thus, the accuracy of the PED obtained from the results of the SMACC SMA, using the 0.3 m pan-sharpened WV3 data, was tested using a photointerpretation approach of random point samples reported on the RGB aero-photos (0.25 m a.g.r.), supported by preliminary unsupervised clustering analysis. To facilitate the photointerpretation approach, the five macro-classes of main rubble materials (Table 1) were grouped into three more visually discriminable categories: • 1st-debris-concrete (characterized by light grey shades); • 2nd-tiles-bricks (characterized by reddish shades); • 3rd-natural stone (characterized by dark grey shades). Asbestos asb Accordingly, the PED was recoded on the three previous groups (categories) through the classification output obtained by the machine learning algorithm to produce a three classes thematic map of the predominant endmembers. Then, the accuracy of this former was estimated using a randomly selected point sample whose classes were independently derived by means of photointerpretation on the corresponding RGB aero-photos. However, it is worth noting that excluding the few pixels exploited for deriving the SMACC endmembers (statistically insufficient for an effective accuracy assessment), most of rubble pixels are a mix of the three categories previously introduced. At the pixel level, the predominant (maximum abundance) endmembers are characterized by different abundances less than 1. In particular, the first category (debris-concrete) was prevalent in the majority of pixels and showed the highest sensible abundance values with respect to the other two classes.
Generally, the endmember maximum abundance values of the others two classes were lower, and these classes were prevalent in a minor number of pixels. Of course, greater is the abundance of the specific category in the pixel and more effective is the photointerpretation and consequent labelling. Therefore, increasingly rising abundance threshold values were adopted in order to provide a statistically sufficient number of samples for all three of the classes introduced above.
Accordingly, three related confusion matrices, one for each selected abundance threshold, were calculated, showing different overall accuracies and k-statistics parameters. Then, using the classification of the endmembers into the 25 classes of the library provided by machine learning classifiers, the percentage distribution of these classes was worked out.

GIS Processing
Finally, the characterization of rubbles within the test area was obtained through zonal GIS processing of the heaps' polygonal cover. In particular, for each heap, the mean values of the fractional abundancies of the three main constituents were assessed. Of course, by multiplying each average value by the number of pixels of the heaps, the surface percentage of the corresponding (surface) materials abundancy, at heap level, is calculated.

Results
Eight principal constituents (macro-classes) of urban rubble were identified in the field during the two-day survey, during which 150 hyperspectral signatures (about 2300 narrow bands for each spectral signature) were acquired. In Table 1, the nine related macro-classes, including the asbestos macro-class, whose spectral signatures were derived from the literature, are reported.
The majority of these field data refer to the more abundant rubble materials, which correspond to the first five macro-classes, whose spectral signatures differed significantly among the various piles and urban places. Therefore, after a preliminary analysis and in order to preserve the spectral variability of these five materials, four different averaged signatures were recorded for each of them.
In Figure 3, the four graphs show: the in situ hyperspectral signatures acquired on concrete (top left graph) and brick (top right graph) after the first pre-processing phase but before smoothing by filtering; the corresponding filtered spectral signatures of the concrete (lower left graph) and brick lower right graph) categories. The bottom image of Figure 3 represents some samples of the eight types of rubble materials. The graph in the upper part of Figure 4 shows the multispectral signatures resampled using the WV3 sensor filter function (the asbestos signatures were omitted for clearness). The multispectral signatures of the 10 endmembers, extracted via SMACC from WV3 data of rubble heaps, are reported in the lower graph of Figure 4. The graph in the upper part of Figure 4 shows the multispectral signatures resampled using the WV3 sensor filter function (the asbestos signatures were omitted for clearness). The multispectral signatures of the 10 endmembers, extracted via SMACC from WV3 data of rubble heaps, are reported in the lower graph of Figure 4. From the WV3 pansharpened pre-processed imagery a small number (10) of pure material signatures, i.e., endmembers, were derived through the SMACC algorithm to be used as reference for the main material classes of urban rubble, by matching them with the resampled multispectral signatures of the in situ materials using machine learning classifiers. The SMACC algorithm thus provided a multilayer distribution of ten endmembers, plus shadow, whose basic statistics are shown in Table 2.
The first five endmembers except end3, plus shadow, are the most abundant, on average, as reported in Table 2. In particular, the shadow shows high abundance values in a relevant number of pixels, as shown in Figure 5, where the distributions of the abundance of the shadow and first five endmembers are depicted.  From the WV3 pansharpened pre-processed imagery a small number (10) of pure material signatures, i.e., endmembers, were derived through the SMACC algorithm to be used as reference for the main material classes of urban rubble, by matching them with the resampled multispectral signatures of the in situ materials using machine learning classifiers. The SMACC algorithm thus provided a multilayer distribution of ten endmembers, plus shadow, whose basic statistics are shown in Table 2.
The first five endmembers except end3, plus shadow, are the most abundant, on average, as reported in Table 2. In particular, the shadow shows high abundance values in a relevant number of pixels, as shown in Figure 5, where the distributions of the abundance of the shadow and first five endmembers are depicted.   After the spectral resampling, each of the 25 resampled signatures, corresponding to one of the nine main macro-classes of urban rubble materials identified on field, was matched (classified) with those of EMs derived from the WV3 pansharpened data using several machine learning algorithms. The related results are reported in Table 3Table 3. Please note that the x suffix, x = 1 to 4, in the first five class labels and x=1 to 2 in the asbestos labels, refers to the four averaged signatures stored in the library for each of the five main rubble materials and to the two signatures available for asbestos, respectively. In the first 25 rows of Table 3Table 3, the train classes corresponding to the 25 spectral signatures stored in the library ("Training" column) and the related class labels assigned by the exploited classifiers (3rd to 7th column) are reported, respectively.  After the spectral resampling, each of the 25 resampled signatures, corresponding to one of the nine main macro-classes of urban rubble materials identified on field, was matched (classified) with those of EMs derived from the WV3 pansharpened data using several machine learning algorithms. The related results are reported in Table 3. Please note that the x suffix, x = 1 to 4, in the first five class labels and x = 1 to 2 in the asbestos labels, refers to the four averaged signatures stored in the library for each of the five main rubble materials and to the two signatures available for asbestos, respectively. In the first 25 rows of Table 3, the train classes corresponding to the 25 spectral signatures stored in the library ("Training" column) and the related class labels assigned by the exploited classifiers (3rd to 7th column) are reported, respectively.

Endmembers
Moreover, in the first row of Table 3, beneath each classifier label, the related error rates (number of wrong classified samples/total number of samples) achieved are reported, indicating C-SVC and RF as the best performing algorithms (minimum error rate; note that RF's 100% suitability against the training set is a priori guaranteed by definition). In the remaining row, the classification of the ten endmembers (end1÷end10), by means of each classifier, are shown.
Taking into account the abovementioned results, the PED categorical distributions corresponding to different minimum PMA thresholds (0.25, 0.35 and 0.45), selected on the basis of the PMA's histogram (Figure 7), were calculated to be used for the accuracy assessment. In this way the three rightmost fractions of the histogram distribution ( Figure 7) were selected to generate the maps whose accuracy was assessed. Taking into account the abovementioned results, the PED categorical distributions corresponding to different minimum PMA thresholds (0.25, 0.35 and 0.45), selected on the basis of the PMA's histogram (Figure 7), were calculated to be used for the accuracy assessment. In this way the three rightmost fractions of the histogram distribution ( Figure 7) were selected to generate the maps whose accuracy was assessed. In Figure 8 the recoded PED distribution obtained using a 0.25 PMA threshold is displayed.  In Figure 8 the recoded PED distribution obtained using a 0.25 PMA threshold is displayed. The accuracy assessment was then performed by comparing the PED classes at randomly selected points with the classes independently derived (for the same points) by photointerpretation. The accuracy assessment results obtained using the three thresholds for PMA, are reported in the following confusion matrices with the related accuracy parameters (Tables 4-6). The accuracy assessment was then performed by comparing the PED classes at randomly selected points with the classes independently derived (for the same points) by photointerpretation. The accuracy assessment results obtained using the three thresholds for PMA, are reported in the following confusion matrices with the related accuracy parameters (Tables 4-6). concrete-debris brick-tile nat. stone Producer Accuracy User Accuracy    The random sample was about 1-2% of the entire population, being at its maximum when the 0.25 threshold was used and decreasing sensibly with the threshold rise ( Figure 6). However, it must be underlined that the samples number effectively exploited for 0.25 PMA threshold is lower than that related to the 0.35 one, due to the exclusion of samples found contaminated by vegetation.
According to the RF classification results, endmembers n. 5 (proxy of the asbestos' distribution) was worked out and reported in Figure 9. The different colours represent the asbestos relative abundance on rubble heaps. The random sample was about 1-2% of the entire population, being at its maximum when the 0.25 threshold was used and decreasing sensibly with the threshold rise ( Figure 6). However, it must be underlined that the samples number effectively exploited for 0.25 PMA threshold is lower than that related to the 0.35 one, due to the exclusion of samples found contaminated by vegetation.
According to the RF classification results, endmembers n. 5 (proxy of the asbestos' distribution) was worked out and reported in Figure 9. The different colours represent the asbestos relative abundance on rubble heaps. It must be highlighted that, probably due to misclassification of aged dark-red tiles erroneously assigned to worn asbestos, which is also caused by diffuse colonization of local roofs by micro flora, the asbestos distribution initially showed unlikely highest concentrations (abundance > 0.3) corresponding to intact tile roofs (rubble heaps n. [23][24][25][26][27][28][29]. Therefore, after a visual analysis, a maximum accepted value of 0.25 was applied to avoid this overestimate, providing a more realistic distribution of asbestos within rubble. Then, as shown in Figure 10, the n. 2, 4, 8, 9, 10, 42 heaps are those characterized by areas with the highest asbestos' concentrations. It must be highlighted that, probably due to misclassification of aged dark-red tiles erroneously assigned to worn asbestos, which is also caused by diffuse colonization of local roofs by micro flora, the asbestos distribution initially showed unlikely highest concentrations (abundance > 0.3) corresponding to intact tile roofs (rubble heaps n. [23][24][25][26][27][28][29]. Therefore, after a visual analysis, a maximum accepted value of 0.25 was applied to avoid this overestimate, providing a more realistic distribution of asbestos within rubble. Then, as shown in Figure 10, the n. 2, 4, 8, 9, 10, 42 heaps are those characterized by areas with the highest asbestos' concentrations. Finally, the zonal processing of the raster map derived from the classification implementing the per pixel distribution of endmembers percentages, provided us with the relative percentages of the three main categories used for the accuracy assessment within the rubble heaps. The estimated distributions of the main materials' typologies are shown in Figure 9. The most abundant rubble heap typology (macro class) is represented by concrete/debris (range: 40-50%), followed by brick/tiles (range: 25-30%) and natural stones (range: 7-17%), with other materials (metal, plastic, bitumen sheath, etc.) deriving from the remaining endmembers completing the list (range: 3-25%).
In Table 7 an example of the quantitative results provided by the whole procedure is reported for several rubble heaps. Indeed, for each heap the procedure makes it possible to determine the volume and the main constituents' percentages, including the information regarding the presence of asbestos. Finally, the zonal processing of the raster map derived from the classification implementing the per pixel distribution of endmembers percentages, provided us with the relative percentages of the three main categories used for the accuracy assessment within the rubble heaps. The estimated distributions of the main materials' typologies are shown in Figure 9. The most abundant rubble heap typology (macro class) is represented by concrete/debris (range: 40-50%), followed by brick/tiles (range: 25-30%) and natural stones (range: 7-17%), with other materials (metal, plastic, bitumen sheath, etc.) deriving from the remaining endmembers completing the list (range: 3-25%). In Table 7 an example of the quantitative results provided by the whole procedure is reported for several rubble heaps. Indeed, for each heap the procedure makes it possible to determine the volume and the main constituents' percentages, including the information regarding the presence of asbestos.

Discussion
The selection and exploitation of machine learning algorithms for matching (classifying) the SMACC endmembers' spectral signatures to those derived from the field hyperspectral signatures was an important and critical step of the work.
The global assessment of the algorithms' performance was accomplished taking into account both the right classification of training materials, summarized by the related error rate (Table 3), and the number/typology of main materials (including asbestos) detected through the endmembers.
The endmembers labelled by BVM classifier accounted for the four main materials (cem, matc, lat, matn) plus gua (see Table 3), but failed to discriminate the training samples. In addition, the results showed an unreliable assignment of the asbestos class. Conversely, the C-SVC performed the best with regard to training sample classification, and was also able to associate the four endmembers with the four main materials. Although the KNN performance in classifying the training set was slightly better than that of BVM, it wasn't able to avoid unrealistic asbestos labelling for the three training samples with respect to other rubble materials. However, this algorithm was able to account for four main rubble materials through the endmembers so as the RF. The latter was, in particular, the only one able to detect the likely presence of asbestos (end5). The related endmember (n. 5) showed an average abundance of few percent but with the highest standard deviation meaning a significant concentration of this dangerous material in specific heaps (Table 5). Finally, the C4.5 was the worst performing, both in classifying the training set and in the number of main materials associated with endmembers.
Following these criteria, the C-SVC and RF were the best performing. Moreover, these two algorithms agreed with the labelling of the five endmembers. They were able to discriminate four main materials with differences in terms of natural stones (matn), which was detected by C-SVC but not by RF, which instead was able to classify rubble debris (calc), which was missed by the other classifier. Therefore, the ten endmembers provided by the SMA from pansharpened WV3 data were satisfyingly mapped to the in situ spectral signatures of the rubble materials using C-SVC machine learning algorithm, which made it possible to obtain the most effective discrimination of the main building materials (first 5 endmembers of Table 1).
In agreement with the general schema of the mapping model, which includes validation and accuracy assessment based on independent measurements, an original method for accomplishing this necessary test on the obtained continuous distribution of rubble typologies was implemented. It is based on the pixel maximum abundance (PMA) and predominant endmember (PED) distributions, derived from the SMA outputs. The PMA histogram distribution ( Figure 6) confirms once again the extreme spectral heterogeneity of rubble, with a maximum of less than 0.3. Such an approach allowed us to assess various PED distributions (obtained using different PMA thresholds) using photointerpretation method, i.e., by comparing corresponding locations on the RGB high-resolution aerophotos. According to the typical thematic accuracy assessment it was possible to evaluate overall accuracy (O.A.) and K parameters through related confusion matrices. As expected, the overall accuracy (O.A.) assessed for the three PMA rising thresholds (Tables 4-6) increases with the threshold increase due to endmember predominance increasing at pixel level, which in turn makes correct photointerpretation easier. Indeed, both O.A. and K accuracy global parameters are maximum (i.e., 0.888 and 0.754, respectively) at the 0.45 threshold, which is in substantial agreement with the photointerpretation results. Although this relative abundance value is low, it seems sufficient to allow the pixel main material to be recognized by the independent expert interpreter.
Additionally, the PED maps obtained with the two lowest thresholds (0.35 and 0.25) correspond to a poorer/moderate agreement with the photointerpretation according to the decrease in endmember predominance [55].
The classification of materials present on the surface of rubble heaps is one of the most important results achieved in the present study. Among these, asbestos is of great interest due to its implications for human health (for both earthquake survivor population and rescue teams). The present work showed that concrete debris and brick tiles are more abundant (Table 7 and Figure 11a) within the study area and about 50% of the rubble heaps may contain asbestos (in different quantities; Table 7 and Figure 11b).
It does not mean that 50% of rubbles contains asbestos in concentrations that could be harmful to human health. However, the present research demonstrate that it is possible to identify the potential presence of this hazardous substance and support management of waste and rubbles to be removed. Sustainable land reclamation is one of the main goals of ENEA as it has a great socioeconomic importance with a high impact on economic development [56]. Results of the present study are of extreme practical and logistical relevance, and it could be replicated to improve the effectiveness of management in post-earthquake emergencies around the world.
to a poorer/moderate agreement with the photointerpretation according to the decrease in endmember predominance [55].
The classification of materials present on the surface of rubble heaps is one of the most important results achieved in the present study. Among these, asbestos is of great interest due to its implications for human health (for both earthquake survivor population and rescue teams). The present work showed that concrete debris and brick tiles are more abundant (Table 7 and Figure 11a) within the study area and about 50% of the rubble heaps may contain asbestos (in different quantities; Table 7 and Figure 11b). It does not mean that 50% of rubbles contains asbestos in concentrations that could be harmful to human health. However, the present research demonstrate that it is possible to identify the potential presence of this hazardous substance and support management of waste and rubbles to be removed. Sustainable land reclamation is one of the main goals of ENEA as it has a great socioeconomic importance with a high impact on economic development [56]. Results of the present study are of extreme practical and logistical relevance, and it could be replicated to improve the effectiveness of management in post-earthquake emergencies around the world.

Conclusions
The radiometric and geometric resolution parameters of the remote sensors play an important role in allowing suitable mapping of the spectral-spatial heterogeneity typical of urban and industrial/infrastructure areas. In particular, they are fundamental for the discrimination of the seismic rubble at the required scale. To achieve this goal, the panchromatic and multispectral information provided by WV3 sensors were fused using the Gram-Schmidt pansharpening method. Subsequently, they were atmospherically corrected in order to obtain a suitable input at 31 cm of a.g.r. for the SMA procedure. In general, even though widely exploited for their improved geometric resolution, the quantitative use of synthetic spectral reflectance of pansharpened products is not widespread, since the degradation of the quality during the fusing processes does not make it possible to fully preserve the original spectral properties of the acquisition channels. In this case, rather effective results were obtained by performing the pansharpening process using the spectrally conservative Gram-Schmidt method prior to atmospheric correction. The inverse procedure (atmospheric correction carried out before pansharpening process) did not produce the quality level in the results.
The implemented methodology provided the relative abundance distribution of the considered rubble materials. The results are in agreement with the information that has been made available until now with respect to transportation and final disposal. Indeed, concrete and rubble debris were the most abundant seismic rubble typologies. These findings are consistent with the typical Figure 11. Data of Table 7 show the main constituents' percentages of m 3 (a). Even if asbestos is detected within five of the fourteen considered rubbles heaps, it can be found within almost 50% of the entire volume ((b); Table 7).

Conclusions
The radiometric and geometric resolution parameters of the remote sensors play an important role in allowing suitable mapping of the spectral-spatial heterogeneity typical of urban and industrial/infrastructure areas. In particular, they are fundamental for the discrimination of the seismic rubble at the required scale. To achieve this goal, the panchromatic and multispectral information provided by WV3 sensors were fused using the Gram-Schmidt pansharpening method. Subsequently, they were atmospherically corrected in order to obtain a suitable input at 31 cm of a.g.r. for the SMA procedure. In general, even though widely exploited for their improved geometric resolution, the quantitative use of synthetic spectral reflectance of pansharpened products is not widespread, since the degradation of the quality during the fusing processes does not make it possible to fully preserve the original spectral properties of the acquisition channels. In this case, rather effective results were obtained by performing the pansharpening process using the spectrally conservative Gram-Schmidt method prior to atmospheric correction. The inverse procedure (atmospheric correction carried out before pansharpening process) did not produce the quality level in the results.
The implemented methodology provided the relative abundance distribution of the considered rubble materials. The results are in agreement with the information that has been made available until now with respect to transportation and final disposal. Indeed, concrete and rubble debris were the most abundant seismic rubble typologies. These findings are consistent with the typical proportion of building materials within the centre of Amatrice and the propensity of debris to spread over heaps surface in case of building collapses (Figure 7, Figure 9). The clean-up and removal of rubble has not yet been completed in the central Apennines, and relevant information obtained in the future will be useful for refining the implemented methodology. For instance, one of the crucial feedbacks will be the correspondence between the material distribution retrieved at the surface and the inner part of the rubble heap.
Although detailed information about the asbestos distribution over rubble heaps was not available, the global presence of this dangerous material, mainly as a constituent of roof covers, has been reported by various independent observations made in the field and by transportation operators in charge of clearing the rubble from the centre of Amatrice. In this context, it was not possible to test the accuracy of the asbestos distribution retrieved through RF machine learning algorithm.
In the context of post-event activities related to seismic events, the effective management of the huge amount of resulting rubble and debris in urban areas represents one of the major issues and can be usefully supported by EO-based monitoring and mapping applications, also exploiting the last generation of multispectral satellite VHR sensors. The surveys from the aerial platform can provide detailed up-to-date information, but are extremely expensive for the continuous monitoring of emergency areas required for controlling the evolution of the earthquake urban scenario, which instead can be repetitively mapped by available sensors on board orbiting satellite platforms, which can provide an increasing amount of territorial data with continuously improving spatial and spectral resolution. Results obtained at the present time through the innovative methodology developed in this work show the applicability of the WV3 satellite VHR data for providing effective surface estimates of the distribution of seismic rubble types in heaps distributed within urban areas. This information can be advantageously integrated with that from related volumes derived from aerial surveys carried out in the aftermath of the earthquake (Tables 2 and 7; Figures 8 and 10), providing suitably support to the management activities related to the urban rubble. The stereoscopic acquisition capabilities of VHR satellite sensors could also support the estimation of volumes in monitoring the evolution of the emergency scenario, possibly with a "one-off" calibration based on the LiDAR information deriving from an initial aerial survey. The preliminary analysis of rubble heaps with EO tools can also have a relevant impact on cost of post emergency action like rubble management and the definition of priority in order to guarantee human health protection.
Thus, this original EO-based methodology allowed us to produce realistic estimates of the types and volumes of seismic rubble in the historical centre of Amatrice. It can be generalized with a high degree of automation and can be easily applied to similar situations, even with possible improvements deriving from the integration of the SAR (Synthetic Aperture Radar) and hyperspectral data provided by ongoing EO missions (Copernicus Sentinel-1, Cosmo SkyMed, PRISMA, EnMap). Future activities will include a greater exploitation of high-resolution and hyperspectral satellite data, for assessing both volumes and typologies of earthquake rubble heaps, also investigating the possibility to detect hazardous rubble materials (i.e., asbestos). recognizable as a well-defined single structure, a unique rubble heap was digitized for each single building; a pooled polygon referring to a group of (not distinguishable) buildings was sketched otherwise.
In addition, a classification of building's typology was performed (assigning a specific class, e.g., "religious", "residential", "educational", "medical", etc.), to each polygon, supported by Copernicus EMS data, as well as other web sources like Google Maps.

4.
To each polygon, the building's base elevation value (B) was assigned: it was derived from the DTM post in correspondence of a visible point (e.g., portion of the surrounding streets free of rubbles), identified close to the considered building block/rubble. In case of a significant slope of the streets contiguous to the building block, the minimum altitude value among different possible points nearby the area of interest was selected, to not underestimate the actual elevation of rubble heaps ( Figure A1).

5.
The basis of the previously described steps, the volume of a single rubble heap can be estimated as the solid ( This calculation was carried out in GIS environment by means of the ERDAS Imagine "zonal attribute" function (Zonal Statistics to Polygon Attributes: this ERDAS Imagine (by Hexagon Geospatial) operator extracts the zonal statistics of the background image of a vector feature layer and save them as vector attributes), which allows to assign the calculated volume value as an attribute of the corresponding GIS polygon representing the rubble heap (in shapefile format). 3. In addition, a classification of building's typology was performed (assigning a specific class, e.g. "religious", "residential", "educational", "medical", etc.), to each polygon, supported by Copernicus EMS data, as well as other web sources like Google Maps.
4. To each polygon, the building's base elevation value (B) was assigned: it was derived from the DTMpost in correspondence of a visible point (e.g., portion of the surrounding streets free of rubbles), identified close to the considered building block/rubble. In case of a significant slope of the streets contiguous to the building block, the minimum altitude value among different possible points nearby the area of interest was selected, to not underestimate the actual elevation of rubble heaps (FigFigure A1).
5. The basis of the previously described steps, the volume of a single rubble heap can be estimated as the solid (3-D space) included between the flat horizontal polygon, i.e. the bottom surface (B) represented by the building's foundation developing at a constant altitude (selected DTMpost value), and the upper uneven surface represented by the DSMpost values over the polygonal area itself. So, the volume is assessed through the difference between the DSMpost and the selected DTMpost value (constant altitude B) over the area identified by the heap's perimeter, as expressed by the formula (1): where B is the elevation value of the rubble heap's base and h is the elevation value (in m) extracted from the DSMpost for each pixel, n = number of pixels. Volume values are expressed in m 3 .
This calculation was carried out in GIS environment by means of the ERDAS Imagine "zonal attribute" function 1 , which allows to assign the calculated volume value as an attribute of the corresponding GIS polygon representing the rubble heap (in shapefile format).

Method 2
Method 2 follows the same steps 1 and 2 of M1 but different step 3, since the assignment of the elevation values of rubble's base is not manually performed. Rather, the base values are automatically extracted by the pre-event DTMpre, which can be assumed as the bottom surface.
Therefore, in M2 processing chain, the volume calculation is directly obtained by the difference in elevation between post-event DSMpost and pre-event (DTMpre over the rubbles' perimeter (step 3), as described in equation (2): Zonal Statistics to Polygon Attributes: this ERDAS Imagine (by Hexagon Geospatial) operator extracts the zonal statistics of the background image of a vector feature layer and save them as vector attributes.

Method 2
Method 2 follows the same steps 1 and 2 of M1 but different step 3, since the assignment of the elevation values of rubble's base is not manually performed. Rather, the base values are automatically extracted by the pre-event DTM pre , which can be assumed as the bottom surface.
Therefore, in M2 processing chain, the volume calculation is directly obtained by the difference in elevation between post-event DSM post and pre-event (DTM pre over the rubbles' perimeter (step 3), as described in Equation (2) The above methodology must be cited as [1].