Detection of Urban Damage Using Remote Sensing and Machine Learning Algorithms: Revisiting the 2010 Haiti Earthquake

...................................................................................................ii ACKNOWLEDGMENTS...................................................................................iii TABLE OF CONTENTS....................................................................................iv LIST OF TABLES AND FIGURES........................................................................v PREFACE......................................................................................................vi MANUSCRIPT TITLE.......................................................................................


Introduction
Earthquakes accounted for over 60% of all natural disaster-related deaths from 2001 to 2011-a danger that will likely increase due to rapid global urbanization [1].Immediately after an earthquake occurs, satellite imagery is a critical component of damage mapping.Hussain et al. noted that "information derived from remote sensing data greatly helps the authorities in rescue and relief efforts, damage assessment, and the planning of remedial measures to safeguard such events effectively" [2].For immediate rescue operations, rapid damage maps derived from satellite imagery must be developed quickly.A study of the 1995 Kobe earthquake in Japan showed a drastic reduction of the total rescued and the proportion of survivors after the third day of recovery efforts [3,4].However, because rapid mapping is required to balance immediacy with in-depth analysis, early mapping efforts often yield coarse damage assessments [5].
Remote sensing has been used widely to map the effects of major disasters such as earthquakes.Numerous studies have utilized electro-optical (EO), synthetic aperture radar (SAR), light detection and ranging (LiDAR), ancillary data, or a combination thereof for post-earthquake damage detection [1,5,6].One technique for damage detection involves fusion of SAR and EO data in pixel-based damage detection.Stramondo et al. used a maximum likelihood (ML) classifier on SAR features derived from the European Remote Sensing mission in combination with EO data provided by the Indian Remote Sensing satellite in order to identify damaged structures following the 1999 Izmit, Turkey earthquake [7].A similar approach combined SAR from COSMO/SkyMed mission and very high resolution (VHR) EO data from the Quickbird satellite to improve damage detection at block level after combining the two datasets in a pixel-based classification following the 2009 L'Aquila earthquake [8].
As early as 1998, object-based image analysis (OBIA) has been used to detect earthquake damage from remote sensing [9].More recently, OBIA has been a continual focus in earthquake detection damage with many studies focusing on the use of unmanned aerial systems (UAS), LiDAR, and the popular image segmentation and classification software eCognition.Hussain et al. [2] fused GeoEye-1 VHR EO data and airborne LiDAR elevation models derived from the RIT-ImageCAT UAS for image segmentation using the Definiens (now eCognition) software suite.The data were classified using nearest neighbor and fuzzy membership sets to detect damaged buildings and rubble following the 2010 Haiti earthquake.Similarly, Pham et al. [6] used aerial VHR RGB composite and LiDAR data (also from the RIT-ImageCAT UAS) along with eCognition for object segmentation and damage detection.
The application of machine learning algorithms (MLAs) to earthquake damage detection is a relatively new area of study.MLAs actively adapt and learn the problem at hand, often mimicking natural or biological systems, instead of relying on statistical assumptions about data distribution [10].In addition to overall improved accuracy [11,12], MLAs have several advantages compared to traditional classification and change detection methods.MLAs work with nonlinear datasets [11,13], learn from limited training data [12,14], and successfully solve difficult-to-distinguish classification problems [15].
Ito et al. [16] used learning vector quantization (LVQ), a type of artificial neural network (ANN) to classify SAR features signifying damage after the 1995 Kobe earthquake.Li et al. [17] used a two-class support vector machine (SVM) on pre-and post-earthquake Quickbird imagery along with spatial relations derived from the local indicator of spatial association (LISA) index to detect structures damaged by the Wenchuan earthquake of 2008.Haiyang et al. utilized a SVM approach in combination with eCognition image segmentation on the RIT-ImageCAT RGB and LiDAR data, as well as the textural features of contrast, dissimilarity, and variance derived from the gray level co-occurrence matrix (GLCM) to detect urban damage in Port-Au-Prince [18].Kaya et al. [19] used OBIA in combination with support vector selection and adaptation (a type of SVM) on pansharpened Quickbird imagery to conduct damage detection for specific buildings within Port-au-Prince after the 2010 earthquake.While OBIA using SVMs have been researched extensively in the past, ANNs, particularly radial basis function neural networks (RBFNNs), and Random Forests (RF) have shown promise in pattern recognition and image classification [15,20] and have yet to be examined in the application of earthquake damage detection.All three algorithms require parameter-tuning process to achieve optimal performance and a cross-validation approach can be applied to automate the parameter-tuning.SVM has an advantage in dealing with small sample size problems due to its sparse characteristics.However, for applications where a large number of training samples are available, SVM often yields a large number of support vectors, resulting in unnecessary complexity and a long training time [21].
Evaluating structural dimensions such as the Laplacian of Gaussian (LoG) and object-based metrics in addition to spectral and textural information could greatly increase damage detection rates.LoG, a blob detection technique, has been used for medical applications in nuclei mapping [22] and for the detection of buildings in bitemporal images [23].As discussed earlier, OBIA has shown strong results in urban scenes and earthquake damage detection.Huang and Zhang had success applying the popular mean-shift segmentation algorithm for urban classification in hyperspectral scenes while statistical region merging (SRM) is another segmentation approach which is robust to noise and occlusions [24,25].Additionally, various metrics such as rectangular fit, morphological shadow index, and morphological building index can describe the structure of objects in the scene before or after segmentation [26,27].Applying structural descriptors such as a LoG filter and segmentation derived metrics to high resolution satellite imagery as an additional input to an MLA could evince damage in difficult to detect scenarios such as a pancake collapse [22,28].The robustness and generalizability of RF and ANN along with the additional dimensions of texture and structure may provide higher accuracies in the face of imperfect input data.
In past disasters, by the time an automated change detection scheme is ready for implementation, a crowdsourced team of visual interpreters is already mapping damaged buildings [29].Dong and Shan mention that while manual digitization of damaged structures requires trained image analysts and is unsuitable for large areas, "visual interpretation remains to be the most reliable and independent evaluation for automated methods" [1].Additionally, many previous studies suggest detection schema which require ancillary data such as UAS products, LiDAR, or GIS databases.Many of these products are unavailable in developing regions where the death toll is highest [30].Using MLAs (RF and ANN), a rapid damage map derived from readily available multispectral imagery could allow for a minimal compromise between time and accuracy and allow first responders to more rapidly allocate their resources in a crisis.RF and ANNs along with derived textural and structural features may provide improved balance between rapid and accurate damage detection.The main purposes of this study include:

•
Assess the performance of neural networks (to include radial basis function neural networks), and Random Forests on very high resolution satellite imagery in earthquake damage detection

•
Investigate the usefulness of structural feature identifiers to include the Laplacian of Gaussian and rectangular fit in identifying damaged regions

Study Area and Data
The earthquake on 12 January 2010 near Port-au-Prince, Haiti, was an exceptionally devastating event.The initial shock of 7.0 M w caused an astounding death toll: between 217,000 and 230,000 were reported dead by the Haitian government [2] and the official estimate has grown to 316,000 [31].Additionally, the earthquake and its subsequent aftershocks caused extensive damage in and around Port-au-Prince including numerous landmark buildings such as the National Palace [2].Because of the high density of collapsed and damaged structures available for training and validation of MLAs, the 2010 Haiti earthquake is an ideal case study to evaluate automated damage detection methods.
For this research, we obtained high resolution multispectral and panchromatic remote sensing data from the DigitalGlobe Foundation.A pre-disaster panchromatic image was acquired in December of 2009 by the WorldView-1 satellite (accessed via DigitalGlobe's EnhancedView system) and post-disaster multispectral and panchromatic images were acquired on 15 January 2010 by the Quickbird 2 satellite.All datasets (see Table 1) were resampled to 0.6 meters using the nearest neighbor technique.Images were atmospherically corrected to top of atmosphere (TOA) reflectance [32] and clipped to the study area of central Port-Au-Prince.Pre-and post-earthquake imagery were coregistered using 15 control points and a third order polynomial transformation with a root mean square error (RMSE) of 0.55 m.
Table 1.Data used for the analysis.Note the similar look angles between the two satellite images used and the relatively short gap between the pre-and post-earthquake images.

Sensor Native Resolution Acquisition Date Look Angle
WorldView-1 0. The study area (Figure 1) was divided into training and validation regions.Training samples for the algorithms were selected by manually digitizing polygons over damaged structures, ensuring that training data represented the input space of the entire study area.A total of 1,214,623 undamaged and 134,327 damaged pixels were used for training.For validation purposes, over 900 buildings were digitized and marked as damaged or undamaged according to the UNOSAT Haiti damage assessment [33].This dataset classifies damage as points classified according to the European Macroseismic Scale developed in 1998 [34].This classification schema defines five damage levels ranging from minor/no damage to destruction [1,5,6,35].In order to simplify the damage detection problem, the UNOSAT points belonging to the most severe damage class (EMS grade V) were extracted for the validation sets and used for the damaged building class.This action is required primarily because visual distinction between classes is difficult even with sub-meter pixel resolution [1].While these validation points have been previously used in several studies [6,36], it is important to note that the UNOSAT/UNITAR dataset was derived through manual interpretation of satellite and aerial imagery and very few points were ground verified.However, the assessment remains the standard for the validation of damage that occurred in the 2010 Haiti earthquake.Macroseismic Scale developed in 1998 [34].This classification schema defines five damage levels ranging from minor/no damage to destruction [1,5,6,35].In order to simplify the damage detection problem, the UNOSAT points belonging to the most severe damage class (EMS grade V) were extracted for the validation sets and used for the damaged building class.This action is required primarily because visual distinction between classes is difficult even with sub-meter pixel resolution [1].While these validation points have been previously used in several studies [6,36], it is important to note that the UNOSAT/UNITAR dataset was derived through manual interpretation of satellite and aerial imagery and very few points were ground verified.However, the assessment remains the standard for the validation of damage that occurred in the 2010 Haiti earthquake.

Texture and Structure
Because several studies have shown that classification and change detection performance can be increased with the addition of spatial information [11,20], textural and structural information was extracted from the pre-and post-earthquake panchromatic images.Figure 2 shows the impact of earthquake damage on two selected textural and structural features.Entropy, energy, dissimilarity, and homogeneity are all second order texture features derived from the GLCM that have been correlated with damage or used as proxies for damage in previous studies [35,[37][38][39].In order to reduce dimensionality and eliminate redundancy, the two consistently correlated GLCM features of Entropy (a measure of gray level randomness) and Dissimilarity (a measure of gray level difference (the square of contrast)) were chosen as texture inputs.
In order to reduce noise, a Gaussian filter (σ = 1) was applied to the panchromatic images before measuring image Entropy and Dissimilarity.A 7 × 7 sliding window was used to compute

Texture and Structure
Because several studies have shown that classification and change detection performance can be increased with the addition of spatial information [11,20], textural and structural information was extracted from the pre-and post-earthquake panchromatic images.Figure 2 shows the impact of earthquake damage on two selected textural and structural features.Entropy, energy, dissimilarity, and homogeneity are all second order texture features derived from the GLCM that have been correlated with damage or used as proxies for damage in previous studies [35,[37][38][39].In order to reduce dimensionality and eliminate redundancy, the two consistently correlated GLCM features of Entropy (a measure of gray level randomness) and Dissimilarity (a measure of gray level difference (the square of contrast)) were chosen as texture inputs.
In order to reduce noise, a Gaussian filter (σ = 1) was applied to the panchromatic images before measuring image Entropy and Dissimilarity.A 7 × 7 sliding window was used to compute the 0 • GLCM and the corresponding texture values were calculated for both the before and after panchromatic images.
To define structural features, a Laplacian of Gaussian (LoG) filter was applied as it is one of the most commonly implemented methods of blob detection.A 2-dimensional LoG filter of size (x, y) can be constructed using: Because the Laplacian computes the second derivative, abruptly changing regions of an image will be highlighted.When combined with a Gaussian smoothing filter, blobs can be detected at different scales defined by σ = (r − 1)/3 where r is the radius of a blob of interest [22].In order to detect buildings of different sizes, a multiscale approach is required.A total of 50 separate convolutions of the LoG filter were applied to the pre-and post-earthquake imagery with sigma adjusted at equal intervals between the values of 15 and 35.Additionally, the LoG filter size was increased from 10 × 10 to a 25 × 25 window in order to accommodate the larger sigma values.The minimum response in scale space was assigned to each pixel due to the fact that LoG filters produce negative responses for bright areas (the majority of buildings produced high reflectance in the panchromatic images).
Remote Sens. 2016, 8, 868 5 of 17 the 0° GLCM and the corresponding texture values were calculated for both the before and after panchromatic images.
To define structural features, a Laplacian of Gaussian (LoG) filter was applied as it is one of the most commonly implemented methods of blob detection.A 2-dimensional LoG filter of size (x, y) can be constructed using: Because the Laplacian computes the second derivative, abruptly changing regions of an image will be highlighted.When combined with a Gaussian smoothing filter, blobs can be detected at different scales defined by σ = (r − 1)/3 where r is the radius of a blob of interest [22].In order to detect buildings of different sizes, a multiscale approach is required.A total of 50 separate convolutions of the LoG filter were applied to the pre-and post-earthquake imagery with sigma adjusted at equal intervals between the values of 15 and 35.Additionally, the LoG filter size was increased from 10 × 10 to a 25 × 25 window in order to accommodate the larger sigma values.The minimum response in scale space was assigned to each pixel due to the fact that LoG filters produce negative responses for bright areas (the majority of buildings produced high reflectance in the panchromatic images).Other structural identifiers can be derived through object-based image analysis (OBIA).While a number of geometric indices are possible through OBIA, this study chose rectangular fit as an input feature due to its obvious possibility for correlation with building shape.The analysis of equal area rectangles, drawn according to object moment, has been used as a more robust version of minimum bounding rectangle comparison [26].Rectangular fit is defined as Other structural identifiers can be derived through object-based image analysis (OBIA).While a number of geometric indices are possible through OBIA, this study chose rectangular fit as an input feature due to its obvious possibility for correlation with building shape.The analysis of equal area rectangles, drawn according to object moment, has been used as a more robust version of minimum bounding rectangle comparison [26].Rectangular fit is defined as where A O is the area of the original object, A R is the area of the equal rectangle (A O = A R ) and A D is the overlaid difference between the equal rectangle and the original object [26].For this study, image segmentation was performed on both the before and after panchromatic images using statistical region merging (SRM) as posited by Nock and Nielsen because of its fast and simple implementation [40].
Setting the scale parameter Q = 512 enabled the many small buildings in the scene to be distinguished.Each object's rectangular fit was calculated and included as an input dimension.

Machine Learning Algorithms
Once all preprocessing steps were taken and textural and structural features were derived, the training and validation datasets were transformed into 14-dimensional arrays consisting of pre-panchromatic; pre-entropy; pre-dissimilarity pre-LoG; pre-rectangular fit; post-panchromatic; post-entropy; post-dissimilarity; post-LoG; post-rectangular fit; blue; green; red, and near infrared multispectral layers (see Table 2).Values in these layers were normalized to fall between the values of −1 and 1 in order to standardize the input and validation layers.While cross-validation was implemented for each MLA, exhaustive parameter analysis was not performed due to the large number of training samples available.All MLA design, algorithm implementation, training and validation were performed using MATLAB release 2015a on a 3.5 GHz quad core processor with 64 GB RAM.
Table 2. List of the input features and their corresponding data sources which were used as predictors for earthquake damage.

Neural Networks
Artificial neural networks (ANN) are of continued interest due to their ability to approximate any function given sufficient neurons in each layer, the flexibility of data distribution, and their reduced computational complexity compared with statistical classification methods [11,13,20,41].A feedforward ANN (also known as multilayer perceptron) assigns small random multiplicative weights and additive biases to input neurons and iteratively adjusts them with each additional training data input, minimizing the surface of the performance function representing the error between the training data and the network's output in order to make the best classification [11,13,41].Neural network design is difficult due to the number of free parameters and ambiguous requirements for network depth and complexity which depend upon the problem at hand [41].As a result, neurons were grown and pruned in various combinations within 1 to 3 hidden layers until the network performance was maximized.The resulting ANN consisted of two hidden layers of 20 neurons each with sigmoid transfer functions, and a binary softmax output layer (see Figure 3).Training was accomplished using backpropagation and the scaled conjugate gradient algorithm to minimize cross-entropy and two-fold cross validation via early stopping (when decrease of cross-entropy in a validation subset ceased) was used to prevent overfitting [41,42].Finally, variable importance was measured by retraining the network an additional 14 times with one of each of the input dimensions set to zero for all data points and recording the change in cross-entropy after a set number of iterations (400).An additional ANN which has been used for classification is the radial basis function neural network (RBFN).The first layer of RBFNs consists of basis functions, centered throughout the input space and the second layer consists of a transfer function which combines the results from the basis functions and classifies them into the categories of interest [41].Similar to the backpropagation network design, the number of basis functions was determined through trial and error and grown until performance was maximized.The radial basis layer consisted of 150 Gaussian functions centered using an unsupervised k-Means approach, which intelligently assigns cluster centers in areas of the input space where high activity exists [43,44].The radial basis outputs fed into the second layer, which consisted of a softmax transfer function, trained by minimizing cross-entropy using the scaled conjugate gradient with early stopping as before.

Random Forests
Random Forests, which was originally proposed by Breiman [45], is an ensemble classifier consisting of a large number of classification and regression trees (CART), where final classification is performed by a winner-takes-all voting system.The algorithm trains each tree on an independently drawn subset of the original data (bootstrap aggregating or bagging) and determines the number of features to be used at each node by the evaluation of a random vector [45,46].Because RF is an ensemble classifier, the Law of Strong Numbers dictates that RF will converge without overfitting the model, so that the computationally optimal number of trees can be found through testing the algorithm.Additional benefits include robustness to outliers and noise and built-in estimates of error and variable importance [45,47].Using MATLAB's TreeBagger function, training was accomplished using 400 (additional tree growth resulted in no further decrease in out-of-bag error) classification trees grown by selecting three variables (at random) for each node split.Additionally, the out-of-bag error was collected and variable importance was measured for comparison with the ANN approach.

Testing and Accuracy Assessement
After training, validation testing was performed on two different datasets.The first dataset consisted of an area including the National Palace just to the south of the training area where building footprints were previously digitized.For this testing area, output from the algorithms was converted to polygons and a building by building accuracy assessment was performed on the first dataset, resulting in a simple confusion matrix.This approach was taken in order to ensure that our accuracy measurements were based on actual objects in the scene rather than a pixel-by-pixel assessment.The An additional ANN which has been used for classification is the radial basis function neural network (RBFN).The first layer of RBFNs consists of basis functions, centered throughout the input space and the second layer consists of a transfer function which combines the results from the basis functions and classifies them into the categories of interest [41].Similar to the backpropagation network design, the number of basis functions was determined through trial and error and grown until performance was maximized.The radial basis layer consisted of 150 Gaussian functions centered using an unsupervised k-Means approach, which intelligently assigns cluster centers in areas of the input space where high activity exists [43,44].The radial basis outputs fed into the second layer, which consisted of a softmax transfer function, trained by minimizing cross-entropy using the scaled conjugate gradient with early stopping as before.

Random Forests
Random Forests, which was originally proposed by Breiman [45], is an ensemble classifier consisting of a large number of classification and regression trees (CART), where final classification is performed by a winner-takes-all voting system.The algorithm trains each tree on an independently drawn subset of the original data (bootstrap aggregating or bagging) and determines the number of features to be used at each node by the evaluation of a random vector [45,46].Because RF is an ensemble classifier, the Law of Strong Numbers dictates that RF will converge without overfitting the model, so that the computationally optimal number of trees can be found through testing the algorithm.Additional benefits include robustness to outliers and noise and built-in estimates of error and variable importance [45,47].Using MATLAB's TreeBagger function, training was accomplished using 400 (additional tree growth resulted in no further decrease in out-of-bag error) classification trees grown by selecting three variables (at random) for each node split.Additionally, the out-of-bag error was collected and variable importance was measured for comparison with the ANN approach.

Testing and Accuracy Assessement
After training, validation testing was performed on two different datasets.The first dataset consisted of an area including the National Palace just to the south of the training area where building footprints were previously digitized.For this testing area, output from the algorithms was converted to polygons and a building by building accuracy assessment was performed on the first dataset, resulting in a simple confusion matrix.This approach was taken in order to ensure that our accuracy measurements were based on actual objects in the scene rather than a pixel-by-pixel assessment.The second validation dataset included a larger area in order to test the algorithms' performance when based on kernel density map rank matching.This secondary assessment followed the procedures used by Tiede et al. [36] and Pham et al. [6] in which the centroid points of the damage polygons are computed and kernel density raster maps are generated for both the test damage points as well as the UNOSAT/UNITAR control data points.The damage density for each dataset was projected onto a 20 m raster grid and the rank value of each cell was computed according to the quartile of density that the cell belonged to.The two maps were overlaid and accuracy was assessed by subtracting the UNOSAT/UNITAR kernel density map from the test kernel density maps, with final output values ranging from −3 (omission error) to +3 (commission error).

Results
Table 3 (supplemented by the confusion matrices in Appendix A) compares overall performance of three MLAs using building-by-building and kernel density accuracy assessment.Our findings showed that the multilayer feedforward ANN outperformed both the RBFNN and random forests with an omission error rate of 37.7%.Figure 4 matches the spatially explicit locations of damage detected by the ANN algorithm with the digitized buildings marked as damaged or undamaged (please refer to Appendix B for the other algorithms' damage maps).Both RBFNN and RF had higher overall accuracies, yet drastically underestimated damage.RBFNN and RF created a high user's accuracy for the damaged building class, but a lower producer's accuracy.The feedforward ANN also had the shortest runtime, an advantage that was primarily gained through the extremely fast implementation of the ANN for testing.The kappa value for each of the algorithms indicated that the distribution of damaged and undamaged buildings could not be accounted for by random chance, however both of the ANNs clearly outperformed Random Forests in this measure as well.For the kernel density accuracy assessment, a good result was measured as a value in the comparison map that ranged between −1 and 1 in accordance with the Tiede et al. [36] approach.RBFNN reached a 90% kernel density map match (shown in Figure 5), outperforming the standard ANN and RF, which were not far behind.This indicates that the radial basis function ANN may be able to generalize to a larger area with greater success than either a feedforward network or Random Forests.Even so, each of the algorithms performed at higher accuracies when generalizing the distribution of damage over a wide area instead of detecting individual, building level damage.
overall accuracies, yet drastically underestimated damage.RBFNN and RF created a high user's accuracy for the damaged building class, but a lower producer's accuracy.The feedforward ANN also had the shortest runtime, an advantage that was primarily gained through the extremely fast implementation of the ANN for testing.The kappa value for each of the algorithms indicated that the distribution of damaged and undamaged buildings could not be accounted for by random chance, however both of the ANNs clearly outperformed Random Forests in this measure as well.For the kernel density accuracy assessment, a good result was measured as a value in the comparison map that ranged between −1 and 1 in accordance with the Tiede et al. [36] approach.RBFNN reached a 90% kernel density map match (shown in Figure 5), outperforming the standard ANN and RF, which were not far behind.This indicates that the radial basis function ANN may be able to generalize to a larger area with greater success than either a feedforward network or Random Forests.Even so, each of the algorithms performed at higher accuracies when generalizing the distribution of damage over a wide area instead of detecting individual, building level damage.As well as examining an algorithm's wide area generalizability, kernel density accuracy assessment also allowed for investigation into areas of common error of both omission and commission.One of the more interesting results was that these areas were common in all three algorithms.A common error of commission occurred in the north-center of the test area and coincided with the development of an internally displaced persons (IDP) camp (see Figure 6).According to the algorithms, this camp broke up the "structure" of the underlying field and increased the randomness of the texture, which led it to misclassifying the area as a damaged building.Figure 7 shows a common area of omission error contained within the map in the central region of the testing area.The underlying cause of error in this region is the scene complexity and high density of small As well as examining an algorithm's wide area generalizability, kernel density accuracy assessment also allowed for investigation into areas of common error of both omission and commission.One of the more interesting results was that these areas were common in all three algorithms.A common error of commission occurred in the north-center of the test area and coincided with the development of an internally displaced persons (IDP) camp (see Figure 6).According to the algorithms, this camp broke up the "structure" of the underlying field and increased the randomness of the texture, which led it to misclassifying the area as a damaged building.Figure 7 shows a common area of omission error contained within the map in the central region of the testing area.The underlying cause of error in this region is the scene complexity and high density of small structures before the earthquake occurred.This preexisting randomness and highly variable structure and texture were difficult for the algorithm to interpret, leading to an error of omission.Even so, with kernel density matching occurring in nearly 90% of cells for each algorithm, a wide area damage density classification was successful for both of the algorithms tested.Finally, variable importance shows similar trends for ANN and RF algorithms.There were a few variables that showed rather different utilization between the algorithms.Figure 8 shows the change in error between each variable and was developed by averaging the assigned variable importance between the pre-and post-earthquake datasets.Overall, the multispectral variables had lower changes in out-of-bag error and cross entropy in comparison to the textural and structural values, however the panchromatic images and the near infrared band was useful to each of the algorithms.Of the two texture measures, entropy was utilized more than dissimilarity in all three algorithms.Finally, variable importance shows similar trends for ANN and RF algorithms.There were a few variables that showed rather different utilization between the algorithms.Figure 8 shows the change in error between each variable and was developed by averaging the assigned variable importance between the pre-and post-earthquake datasets.Overall, the multispectral variables had lower changes in out-of-bag error and cross entropy in comparison to the textural and structural values, however the panchromatic images and the near infrared band was useful to each of the algorithms.Finally, variable importance shows similar trends for ANN and RF algorithms.There were a few variables that showed rather different utilization between the algorithms.Figure 8 shows the change in error between each variable and was developed by averaging the assigned variable importance between the pre-and post-earthquake datasets.Overall, the multispectral variables had lower changes in out-of-bag error and cross entropy in comparison to the textural and structural values, however the panchromatic images and the near infrared band was useful to each of the algorithms.Of the two texture measures, entropy was utilized more than dissimilarity in all three algorithms.Rectangular fit was marked as important for both the ANN and RF however the Laplacian of Gaussian filter was more impactful on the RBFNN assessment and also had a moderately high impact on RF performance.

Discussion
It is difficult to determine outright which algorithm is better.While the multilayer ANN outperformed RF in building by building assessments and required slightly less overall training time, it required a large number of training samples in order to perform well.This is not necessarily the case for the RF algorithm.Also, it is rather easy to overfit an ANN to the training data, which can be avoided using RF due to its nature as an ensemble classifier [41,45].Finally, ANNs can also become stuck at local minima in the performance surface without reaching the global solution, yielding an insufficient result [41].However, in our study, the multilayer ANN had the lowest rates of error in detecting damaged buildings, without sacrificing much performance in wide area generalization or overall accuracy.SVMs, while not examined here, have shown promise in earthquake detection in previous studies [17][18][19].While our study focused on ANNs and RF, as little research has been done on their applicability to earthquake detection problems, future studies may investigate the performance of these algorithms (to include SVM) with respect to training sample size.
Beyond damage detection performance, practical considerations require an investigation into time complexity, particularly when considering any kind of operationalization of an algorithm in automatic damage detection.RF took much longer than either of the ANNs to train and test the datasets.The time complexity of a single classification and regression tree is O(mn•logn) where m is the total number of variables and n is the number of samples [48].Because RF is an ensemble classifier, the overall time complexity of Random Forests can be summarized as O(M(mn•logn)) where M is the number of trees grown.For a large number of samples with moderate dimensionality, this can be quite time consuming.In contrast, neural network complexity is highly dependent on network architecture.Time complexity for the scaled conjugate descent algorithm is often polynomial, overall complexity is determined by the problem and the number of free parameters (weights, biases, or basis functions (in the case of RBFNN)) required to describe that problem [49].As such, testing showed that the ANNs trained faster and tested faster, which is important to consider given the requirement to process a potentially large amount of imagery in an operational context.
As previously discussed, a number of preprocessing steps are required to develop each of the

Discussion
It is difficult to determine outright which algorithm is better.While the multilayer ANN outperformed RF in building by building assessments and required slightly less overall training time, it required a large number of training samples in order to perform well.This is not necessarily the case for the RF algorithm.Also, it is rather easy to overfit an ANN to the training data, which can be avoided using RF due to its nature as an ensemble classifier [41,45].Finally, ANNs can also become stuck at local minima in the performance surface without reaching the global solution, yielding an insufficient result [41].However, in our study, the multilayer ANN had the lowest rates of error in detecting damaged buildings, without sacrificing much performance in wide area generalization or overall accuracy.SVMs, while not examined here, have shown promise in earthquake detection in previous studies [17][18][19].While our study focused on ANNs and RF, as little research has been done on their applicability to earthquake detection problems, future studies may investigate the performance of these algorithms (to include SVM) with respect to training sample size.
Beyond damage detection performance, practical considerations require an investigation into time complexity, particularly when considering any kind of operationalization of an algorithm in automatic damage detection.RF took much longer than either of the ANNs to train and test the datasets.The time complexity of a single classification and regression tree is O(mn•logn) where m is the total number of variables and n is the number of samples [48].Because RF is an ensemble classifier, the overall time complexity of Random Forests can be summarized as O(M(mn•logn)) where M is the number of trees grown.For a large number of samples with moderate dimensionality, this can be quite time consuming.In contrast, neural network complexity is highly dependent on network architecture.Time complexity for the scaled conjugate descent algorithm is often polynomial, overall complexity is determined by the problem and the number of free parameters (weights, biases, or basis functions (in the case of RBFNN)) required to describe that problem [49].As such, testing showed that the ANNs trained faster and tested faster, which is important to consider given the requirement to process a potentially large amount of imagery in an operational context.
As previously discussed, a number of preprocessing steps are required to develop each of the textural and structural dimensions.Also, a k-means unsupervised clustering was used as part of the RBFNN algorithm to intelligently center the basis functions before training.Each of these steps adds time and complexity to the final product.For future studies, the parallelization of many of these processes is one way to greatly reduce computational time.Our data were gathered using serial processing (primarily because a parallel implementation of k-means was not immediately available) in order to establish a baseline and fairly assess each algorithm, however parallel implementations (which include graphics processing units [GPUs]) of both ANN and RF training are readily available for use and will greatly speed up training and implementation of these machine learning algorithms.
The actual results from this study mirrored what was expected quite well; the areas of imagery where texture and structure were broken up were often identified as damage, as one would expect.As mentioned in the results section, one interesting finding was that each of the algorithms erroneously identified IDP camp areas as building damage.These IDP camps are ad hoc structures (tents, tarps, and shanties) built primarily on open spaces.As these IDP camps were erected, they broke up the coherent texture and structure of the underlying terrain, causing the algorithm to mark them as damage.While this is technically an error of commission, it is nevertheless a useful result in showing the power of MLAs in seeking out patterns as well as their ability to simultaneously detect damage and displaced persons.In an operational context, the MLA results in combination with a priori knowledge of building distribution via a GIS database would allow first responders and emergency planners to easily distinguish damaged structures from these IDP camps.
As the experiment on variable importance showed, the textural and structural features were some of the most important factors which allowed both ANNs and RF to detect damage and IDP camps.Stramondo et al. [7] also used important textural features for earthquake damage detection although in their study a maximum likelihood classifier was used.This line of thinking, paramount to computer vision applications, is expanded here by using more intelligent algorithms and readily available data.The importance of the panchromatic features along with the texture and structure products derived exclusively from that panchromatic imagery presupposes that future implementations of machine learning may be able to perform earthquake damage detection from panchromatic imagery alone.One reason that the multispectral imagery was not a driving variable is simply due to resolution; the native 2.4 m is too coarse to detect many of the features associated with earthquake damage.Interestingly, the only multispectral product which was found to be an important variable for each of the algorithms was the near infrared band, which may have resulted from a correlation between exposed rubble and a higher near infrared reflectance.These findings may guide future research in determining which variables to focus on in earthquake damage studies.
Our focus on simple panchromatic imagery is a departure from many previous earthquake studies.The state-of-the-art focuses on LiDAR, SAR, unmanned aerial vehicles, and the software suite eCognition [1].However, the access and availability of these additional data requirements may be limited in the aftermath of a destructive earthquake in a less developed region such as Haiti.A return to easily accessible data products such as multispectral or even panchromatic imagery alone could allow a MLA (potentially even one that is pre-trained) to detect imagery without the requirement of ancillary data.One potential disadvantage of the reliance on bitemporal VHR imagery is the requirement for precise coregistration.Different look angles can cause problems in classification and change detection.While image registration is still important to our study, a small look angle difference may not be critical due to our use of textural and structural features rather than the VHR imagery alone.Additionally, registration errors can be seen as a source of noise in the system; each of the MLAs used has been shown to be robust to noise [41,45].The difference in our look angles (~7 • ) did not appear to cause any damage detection errors in a visual survey of our results.Future research may investigate the limits of acceptable look angle differences or use a complex coregistration approach to eliminate the issue altogether [50].
The future of earthquake damage detection may lie in deep convolutional neural networks (DCNN) coupled with high performance computing and GPUs.DCNNs have already pushed the boundaries of artificial intelligence and image recognition; rather than being told which textural, structural, or spectral inputs should be used, these networks automatically learn and identify the defining features (convolutions) of the problem at hand in order to classify future samples [51,52].Initial results are promising.Using MatConvNet (a deep learning library for Matlab), we experimented with training a DCNN with the VGG-F architecture following the approach and using the hyperparameters described by Chatfield et al. [52,53].We segmented the post-earthquake pansharpened image using SRM and trained the DCNN on each labeled, extracted object.The DCNN was not only able to detect buildings at a comparable rate (>55% detection rate), it was able to distinguish between damage and IDP camps (>65% detection rate) and did so using an after-only pansharpened image, reducing data requirements and eliminating the need for coregistration.A pre-trained DCNN optimized specifically for earthquake detection may offer a robust and operationally implementable solution to the much studied topic of earthquake damage detection.

Conclusions
This study analyzed the use of machine learning algorithms to include feedforward neural networks, radial basis function neural networks, and Random Forests in detecting earthquake damage caused by the 12 January 2010 event near Port-au-Prince, Haiti.The algorithms' efficacy was improved by providing coregistered 0.6 m multitemporal imagery, texture features (dissimilarity and entropy), and structure features (Laplacian of Gaussian and rectangular fit) as inputs.Detection results were assessed on a structure specific basis by digitizing more than 900 buildings and comparing each MLA's response to the UNITAR/UNOSAT validation set.For a wide area generalization, a kernel density map comparison was performed between each of the algorithms' classification results and the UN damage validation points.
The feedforward ANN consisting of two hidden layers had the lowest error rate (<40%) without sacrificing much performance in overall accuracy or generalization to wider area damage estimates.Additionally, textural and structural features derived from panchromatic imagery were shown to be more important than spectral variables in the algorithms' classification process.Each algorithm had common errors of commission occurring around ad hoc IDP camps that were spontaneously formed in open spaces following the earthquake; this technically incorrect result can be easily integrated with a GIS layer containing building footprints.
The results of this study show that not only do MLAs have potential for use in earthquake damage detection, but that panchromatic or pansharpened imagery can be the exclusive data source for training and testing.Measures of variable importance found that the panchromatic derived texture and structure products are the main drivers behind the success of these "shallow" machine learning algorithms.Future research into an operationally implementable machine learning method is warranted.An attractive next step is to transition into deep learning where convolutional neural networks move beyond pixel-based or object-based paradigms and begin to detect remotely sensed features in ways akin to natural image recognition.
for supporting this publication.Finally, we would like to thank Earl and Marion Nutter; without their entrusted scholarship this study would not be possible.

Figure 1 .
Figure 1.Study area in Port-au-Prince showing training, building test and kernel density test sites (satellite image courtesy of the DigitalGlobe Foundation).

Figure 1 .
Figure 1.Study area in Port-au-Prince showing training, building test and kernel density test sites (satellite image courtesy of the DigitalGlobe Foundation).

Figure 2 .
Figure 2.This figure comparing building damage shows the effects of an entropy filter and a LoG filter on pre-and post-earthquake imagery considering two different collapse scenarios: (a) a normal structural collapse; (b) a pancake collapse.Pre-earthquake satellite images © copyright DigitalGlobe.Post-earthquake satellite images courtesy of the DigitalGlobe Foundation.

Figure 2 .
Figure 2.This figure comparing building damage shows the effects of an entropy filter and a LoG filter on pre-and post-earthquake imagery considering two different collapse scenarios: (a) a normal structural collapse; (b) a pancake collapse.Pre-earthquake satellite images © copyright DigitalGlobe.Post-earthquake satellite images courtesy of the DigitalGlobe Foundation.
Remote Sens. 2016, 8, 868 7 of 17set to zero for all data points and recording the change in cross-entropy after a set number of iterations (400).

Figure 3 .
Figure 3.A simplified design layout of the feedforward ANN used for training and testing.For input layer detail, see Table2.In reality each of the twenty neurons in Hidden Layer 1 is connected to each of the twenty neurons in Hidden Layer 2.

Figure 3 .
Figure 3.A simplified design layout of the feedforward ANN used for training and testing.For input layer detail, see Table2.In reality each of the twenty neurons in Hidden Layer 1 is connected to each of the twenty neurons in Hidden Layer 2.

Figure 4 .
Figure 4. Results of the feedforward ANN (highest performer) on the building test dataset.Direct output of the algorithm (pink) overlays the building categories are represented by shaded polygons (satellite image courtesy of the DigitalGlobe Foundation).

Figure 4 .
Figure 4. Results of the feedforward ANN (highest performer) on the building test dataset.Direct output of the algorithm (pink) overlays the building categories are represented by shaded polygons (satellite image courtesy of the DigitalGlobe Foundation).

Figure 5 .
Figure 5. Results of the RBFNN algorithm on the kernel density test dataset.Each 20 m cell is marked from −3 to 3, where an adequate density match falls between −1 and 1 (satellite image courtesy of the DigitalGlobe Foundation).

Figure 5 .
Figure 5. Results of the RBFNN algorithm on the kernel density test dataset.Each 20 m cell is marked from −3 to 3, where an adequate density match falls between −1 and 1 (satellite image courtesy of the DigitalGlobe Foundation).

Figure 6 .
Figure 6.Results of the RBFNN on the kernel density test dataset overlaid before (a); and after (b) panchromatic imagery.This ad hoc IDP caused an error of commission among each of the algorithms.Pre-earthquake satellite image copyright © DigitalGlobe.Post-earthquake satellite image courtesy of the DigitalGlobe Foundation.

Figure 7 .
Figure 7. Results of the RBFNN on the kernel density test dataset overlaid before (a); and after (b) panchromatic imagery.This area of complex small structures was a common error of omission among the algorithms.Pre-earthquake satellite image copyright © DigitalGlobe.Post-earthquake satellite image courtesy of the DigitalGlobe Foundation.

Figure 6 .Figure 6 .
Figure 6.Results of the RBFNN on the kernel density test dataset overlaid before (a); and after (b) panchromatic imagery.This ad hoc IDP caused an error of commission among each of the algorithms.Pre-earthquake satellite image copyright © DigitalGlobe.Post-earthquake satellite image courtesy of the DigitalGlobe Foundation.

Figure 7 .
Figure 7. Results of the RBFNN on the kernel density test dataset overlaid before (a); and after (b) panchromatic imagery.This area of complex small structures was a common error of omission among the algorithms.Pre-earthquake satellite image copyright © DigitalGlobe.Post-earthquake satellite image courtesy of the DigitalGlobe Foundation.

Figure 7 .
Figure 7. Results of the RBFNN on the kernel density test dataset overlaid before (a); and after (b) panchromatic imagery.This area of complex small structures was a common error of omission among the algorithms.Pre-earthquake satellite image copyright © DigitalGlobe.Post-earthquake satellite image courtesy of the DigitalGlobe Foundation.

Figure 8 .
Figure 8. Chart of variable importance for both ANNs and RF.Overall, structure and texture played a larger role in classification than spectral information.Random Forests use the ΔOOB Error measure while ANN and RBFNN use ΔCross Entropy.

Figure A2 .
Figure A2.Results of the RF algorithm on the building test dataset.Direct output of the algorithm (pink) overlays the building categories are represented by shaded polygons (satellite image courtesy of the DigitalGlobe Foundation).

Table 3 .
Results from the test datasets.

Table 3 .
Results from the test datasets.
Chart of variable importance for both ANNs and RF.Overall, structure and texture played a larger role in classification than spectral information.Random Forests use the ∆OOB Error measure while ANN and RBFNN use ∆Cross Entropy.