High-Resolution and Accurate Topography Reconstruction of Mount Etna from Pleiades Satellite Data

The areas characterized by dynamic and rapid morphological changes need accurate topography information with frequent updates, especially if these are populated and involve infrastructures. This is particularly true in active volcanic areas such as Mount (Mt.) Etna, located in the northeastern portion of Sicily, Italy. The Mt. Etna volcano is periodically characterized by explosive and effusive eruptions and represents a potential hazard for several thousands of local people and hundreds of tourists present on the volcano itself. In this work, a high-resolution, high vertical accuracy digital surface model (DSM) of Mt. Etna was derived from Pleiades satellite data using the National Aeronautics and Space Administration (NASA) Ames Stereo Pipeline (ASP) tool set. We believe that this is the first time that the ASP using Pleiades imagery has been applied to Mt. Etna with sub-meter vertical root mean square error (RMSE) results. The model covers an area of about 400 km2 with a spatial resolution of 2 m and centers on the summit portion of the volcano. The model was validated by using a set of reference ground control points (GCP) obtaining a vertical RMSE of 0.78 m. The described procedure provides an avenue to obtain DSMs at high spatial resolution and elevation accuracy in a relatively short amount of processing time, making the procedure itself suitable to reproduce topographies often indispensable during the emergency management case of volcanic eruptions.


Introduction
Topography is the result of many Earth internal and external processes that interact and influence each other [1,2]. An accurate three-dimensional (3D) reconstruction of the Earth's surface is often essential to better understand some of these processes and to predict and model potential natural hazards [3][4][5][6][7][8][9]. Among the areas most affected by natural hazards, active volcanic areas are some of the most fascinating. In fact, they are considered as natural laboratories to study the phenomena connected In this work, the Pleiades stereo imagery from July 28, 2015 was processed to generate a DSM over Mt. Etna with submeter vertical RMSE using the National Aeronautics and Space Administration (NASA) Ames Stereo Pipeline (ASP). The ASP is a suite of open source tools intended for processing satellite, aerial, and robotic rover imagery [29]. After a short geo-vulcanological introduction of the study area, we present the datasets and the methodology used to derive the DSM, followed by the validation of the model, and, finally, the discussion and conclusions highlighting and summarizing the main results.

Study Area
Mt. Etna is the highest active continental volcano in Europe (3329.6 m above sea level as measured by an airborne laser scanning survey in the summer of 2007) [20], and it emerges along the northeastern coast of Sicily (Southern Italy), covering an area of about 1200 km 2 considering the basal boundary ( Figure 1). The area is densely populated with hundreds of thousands of local people and tourists living near the volcano and visiting it. Mt. Etna is an active quaternary composite volcano formed by the accumulation of lava and pyroclastic products erupted during the last 500 thousand years [30] and characterized mainly by Na-alkaline magmatism. During the last 15 years, the volcano has had a frequent eruptive activity both effusive [31] and explosive [13,32]. Since 2005, the most remarkable eruptive events were recorded between 2006-2007 and 2008-2009 (this is the longest flank eruption of Etna after the 1991-1993 eruption), 2011,2014,2015,2017, and December 2018. Since 2011, the summit activity of Etna has been characterized by more than 50 paroxysmal episodes, accompanied by lava flow output from one of its major craters [33], and minor effusive eruptions [24] from other craters. On the whole, the lava flows and effusive eruptions affected mostly the South Eastern Crater, in the side of which a new crater developed, named the New Southeast Crater. Additionally, this dynamic activity modified other few adjacent craters such as Voragine, Bocca Nuova, and the Northeast Crater [33,34]. This caused significant and rapid morphological changes of the summit area, and, in some cases, produced pyroclastic density currents [35]. The frequency of eruptive phenomena has particularly increased in the last few years, confirming the importance of continual updates of the topography of Mt. Etna, a densely urbanized volcano [36]. The area of Mt. Etna investigated in this work covers approximately 400 km 2 (yellow box in Figure 1) and covers the summit and the Valle del Bove to the east, the main area affected by eruptions. measured by an airborne laser scanning survey in the summer of 2007) [20], and it emerges along the northeastern coast of Sicily (Southern Italy), covering an area of about 1200 km 2 considering the basal boundary ( Figure 1). The area is densely populated with hundreds of thousands of local people and tourists living near the volcano and visiting it. Mt. Etna is an active quaternary composite volcano formed by the accumulation of lava and pyroclastic products erupted during the last 500 thousand years [30] and characterized mainly by Na-alkaline magmatism. During the last 15 years, the volcano has had a frequent eruptive activity both effusive [31] and explosive [13,32]. Since 2005, the most  remarkable eruptive events were recorded between 2006-2007 and 2008-2009 (this is the longest flank  eruption of Etna after the 1991-1993 eruption), 2011,2014,2015,2017, and December 2018. Since 2011, the summit activity of Etna has been characterized by more than 50 paroxysmal episodes, accompanied by lava flow output from one of its major craters [33], and minor effusive eruptions [24] from other craters. On the whole, the lava flows and effusive eruptions affected mostly the South Eastern Crater, in the side of which a new crater developed, named the New Southeast Crater. Additionally, this dynamic activity modified other few adjacent craters such as Voragine, Bocca Nuova, and the Northeast Crater [33,34]. This caused significant and rapid morphological changes of the summit area, and, in some cases, produced pyroclastic density currents [35]. The frequency of eruptive phenomena has particularly increased in the last few years, confirming the importance of continual updates of the topography of Mt. Etna, a densely urbanized volcano [36]. The area of Mt. Etna investigated in this work covers approximately 400 km 2 (yellow box in Figure 1) and covers the summit and the Valle del Bove to the east, the main area affected by eruptions.

Pleiades Satellite Data
The Pleiades high-resolution satellite data are acquired by panchromatic and multispectral sensors onboard the polar orbiting Pleiades 1A and Pleiades 1B satellites in operation since 2012. The sensors work in the visible (VIS) and the near-infrared (NIR) spectral range. The panchromatic sensor acquires images at 50 cm ground resolution while the multispectral sensor acquires at 2 m. Pleiades was initially developed as part of the French-Italian Optical and Radar Federated Earth Observation program (ORFEO). Other European partners contributed to the sensor construction such as Austria, Belgium, Spain, and Sweden. The main characteristics of the Pleiades constellation are detailed in the Pleiades Imagery User Guide [37].
The two pairs of stereo Pleiades imagery used for this study cover approximately 240 and 400 km 2 , respectively and were acquired on July 28, 2015 as a multiband VIS-NIR processed at the sensor level and pansharpened at 1 m resolution. The first image in the stereo pair acquisition (1&2) covers only 20 by 12 km (blue box, Figure 2), while the second stereo pair (2&3) images covers~20 by 20 km (yellow box, Figure 2).

Pleiades Satellite Data
The Pleiades high-resolution satellite data are acquired by panchromatic and multispectral sensors onboard the polar orbiting Pleiades 1A and Pleiades 1B satellites in operation since 2012. The sensors work in the visible (VIS) and the near-infrared (NIR) spectral range. The panchromatic sensor acquires images at 50 cm ground resolution while the multispectral sensor acquires at 2 m. Pleiades was initially developed as part of the French-Italian Optical and Radar Federated Earth Observation program (ORFEO). Other European partners contributed to the sensor construction such as Austria, Belgium, Spain, and Sweden. The main characteristics of the Pleiades constellation are detailed in the Pleiades Imagery User Guide [37].
The two pairs of stereo Pleiades imagery used for this study cover approximately 240 and 400 km 2 , respectively and were acquired on July 28, 2015 as a multiband VIS-NIR processed at the sensor level and pansharpened at 1 m resolution. The first image in the stereo pair acquisition (1&2) covers only 20 by 12 km (blue box, Figure 2), while the second stereo pair (2&3) images covers ~20 by 20 km (yellow box, Figure 2).  landscape considering all natural and anthropogenic elements present in the investigated area, while the second reproduces only the elevation surface of bare ground. Both models are characterized by high vertical and horizontal accuracy [18] and were processed in the World Geodetic System 1984 (WGS84) geoid projected into Universal Transverse Mercator (UTM) 33 North. Each lidar DEM covers an area of 620 km 2 (red box in Figure 1) with a spatial resolution of 2 m and a vertical RMSE of 0.24 m. In this work, the lidar DTM published in [18] (Figure 2) was used to assign orthometric heights to the Pleiades generated 3D point clouds through data alignment and to decrease vertical errors. The lidar DSM was used as the base to select the GCPs from the geodetic monitoring network of Mt. Etna for Pleiades DSM validation in areas with minimum deformation.

Ground Control Points (GCPs) Dataset
The dataset used to validate the Pleiades DSM consisted of 41 GCPs (Figure 2) selected from the geodetic monitoring network of Mt. Etna [38], according to the following criteria: (i) spatial distribution of the points as homogeneous as possible over the entire overlapping area between the Pleiades and lidar DSM; (ii) GCP sites not affected by morphological changes or land use variations between 2005 and 2015; and (iii) clear recognition of the same GCP sites in Pleiades and lidar DSM. The GCPs have 5 cm coordinate accuracy and are established on six main ground classes [39]: lava flow ground, lava flow ground with sparse and low vegetation, scoriae ground, scoriae ground with sparse and low vegetation, anthropogenic features, and sparse vegetated ground. Points in Valle del Bove were not taken into consideration for this dataset ( Figure 2) due to several lava flows deposited in this area between 2005 to 2015. This dataset has elevation recorded in both ellipsoid and orthometric heights, and thus can be used to validate DSMs in both vertical datums.

Workflow and Methodology
This section illustrates the workflow and methodology used to reconstruct, for the first time, an orthometric high resolution DSM of Mt. Etna by using Pleiades satellite imagery and the Ames Stereo Pipeline (ASP) software, a suite of open source automated geodesy and stereo-photogrammetry tools implemented by NASA. The workflow is described in Figure 3 and differs from the standard processing workflow described by [40] in the sequence of some steps, optional functions used, and the use of different feature matching algorithms and procedures to improve the vertical accuracy. Additionally, it is the first time that this software was methodically and exhaustively tested using Pleiades satellite imagery. The ASP supports imagery that uses the rational polynomial coefficient (RPC) camera models. Previously, ASP was tested and calibrated for extraterrestrial data and later for Earth imagery, specifically DigitalGlobe data. The workflow proposed here included three main phases ( Figure 3) that are detailed in the following subsections. The Mt. Etna 2005 airborne lidar survey acquired 250 * 10 6 3D points on September 29-30 in order to reconstruct two DEMs: the DSM and DTM. The first reproduces the elevation surface of the landscape considering all natural and anthropogenic elements present in the investigated area, while the second reproduces only the elevation surface of bare ground. Both models are characterized by high vertical and horizontal accuracy [18] and were processed in the World Geodetic System 1984 (WGS84) geoid projected into Universal Transverse Mercator (UTM) 33 North. Each lidar DEM covers an area of 620 km 2 (red box in Figure 1) with a spatial resolution of 2 m and a vertical RMSE of 0.24 m. In this work, the lidar DTM published in [18] (Figure 2) was used to assign orthometric heights to the Pleiades generated 3D point clouds through data alignment and to decrease vertical errors. The lidar DSM was used as the base to select the GCPs from the geodetic monitoring network of Mt. Etna for Pleiades DSM validation in areas with minimum deformation.

Ground Control Points (GCPs) Dataset
The dataset used to validate the Pleiades DSM consisted of 41 GCPs ( Figure 2) selected from the geodetic monitoring network of Mt. Etna [38], according to the following criteria: (i) spatial distribution of the points as homogeneous as possible over the entire overlapping area between the Pleiades and lidar DSM; (ii) GCP sites not affected by morphological changes or land use variations between 2005 and 2015; and (iii) clear recognition of the same GCP sites in Pleiades and lidar DSM. The GCPs have 5 cm coordinate accuracy and are established on six main ground classes [39]: lava flow ground, lava flow ground with sparse and low vegetation, scoriae ground, scoriae ground with sparse and low vegetation, anthropogenic features, and sparse vegetated ground. Points in Valle del Bove were not taken into consideration for this dataset ( Figure 2) due to several lava flows deposited in this area between 2005 to 2015. This dataset has elevation recorded in both ellipsoid and orthometric heights, and thus can be used to validate DSMs in both vertical datums.

Workflow and Methodology
This section illustrates the workflow and methodology used to reconstruct, for the first time, an orthometric high resolution DSM of Mt. Etna by using Pleiades satellite imagery and the Ames Stereo Pipeline (ASP) software, a suite of open source automated geodesy and stereo-photogrammetry tools implemented by NASA. The workflow is described in Figure 3 and differs from the standard processing workflow described by [40] in the sequence of some steps, optional functions used, and the use of different feature matching algorithms and procedures to improve the vertical accuracy. Additionally, it is the first time that this software was methodically and exhaustively tested using Pleiades satellite imagery. The ASP supports imagery that uses the rational polynomial coefficient (RPC) camera models. Previously, ASP was tested and calibrated for extraterrestrial data and later for Earth imagery, specifically DigitalGlobe data. The workflow proposed here included three main phases ( Figure 3) that are detailed in the following subsections.

Data Pre-Processing
This phase selects and prepares the data for the feature matching procedure. The two pairs of Pleiades images called 1&2 and 2&3 were selected (for details see Section 2.2.1) and processed eschewing the bundle adjustment step because the Pleiades sensor ephemeris/attitude information is precise [37]. Otherwise, the sensor ephemeris/attitude information needs to be calculated by applying the bundle adjustment algorithm described in [41]. After the stereo pair selection, band 1 was extracted from each image and aligned using the affine-epipolar algorithm [42]. This pre-aligns the images by detecting tie-points between them. The tie-points are used to rotate the sensor position in order to have pairs of conjugate epipolar lines collinear and parallel to one of the image axes. The ASP software does not work with multi-band imagery, so we selected band 1 because it is characterized by a wide range of radiance data with good visual contrast useful in discriminating image features from their shadow. In the active summit zone, this band allows for the identification of areas affected by volcanic plume, fumaroles, and clouds, as they were recognized by Landsat 8, band 9, acquired on the same date as the Pleiades data [43].

Feature Matching
The feature matching phase prepares the data for 3D point cloud generation (Figure 3). Two main algorithms, one using a local search window (LSW), and the other using a more global matching (MGM) procedure, a modification of the semi-global matching (SGM) algorithm, were used [44,45]. The LSW algorithm estimates the apparent motion of the scene points, or disparity by comparing the local windows around each pixel in the stereo pair. This method does not address ambiguities in the image texture such as smoothness or repetitive patterns. Instead, the MGM algorithm compensates for ambiguities by using the smoothness of the disparity map and reducing the two-dimensional smoothness constraint to the average of the one-dimensional line optimization problem. Each stereo pair was processed with both algorithms resulting in relative disparity maps.

Topography Reconstruction
The stereo triangulation method ( Figure 3) is used to combine spacecraft ephemeris/attitude information, sensor model, and the disparity map to obtain 3D locations (x, y, z) of the closest intersection between lines that connect the sensor orbital position to all matched pixels in both the left and right images. The result is a 4-band raster format where band 1 stores triangulated x-coordinates (m), band 2 y-coordinates (m), band 3 the ellipsoid z-coordinates (m), and band 4 a triangulation error metric (i.e., the distance in meters between the two rays at the closest intersection). Band 4 permits the evaluation of the quality of the disparity matches, the sensor model, and ephemeris/attitude data. Subsequently, the result was transformed into a 3D point cloud geo-coded into the WGS84 UTM 33N coordinate system with ellipsoid heights.
In addition to the ASP standard procedure [40], a supplementary step was included to assign orthometric heights to the output 3D point cloud in post-processing (3D reconstruction in Figure 3) by testing two different procedures. The first is based on the alignment of the 3D point cloud to a more accurately positioned (if potentially sparser) dataset; in this case, the lidar DTM at a 2 m spatial resolution. The alignment procedure uses the iterative closest point (ICP) algorithm [46] that eliminates extreme outliers when their displacement is greater than a defined threshold. This threshold assumes the mean value calculated on 75% of the smallest displacements. For Pleiades derived point cloud data, the displacement value was 50 m after investigating the error values reported by the alignment algorithm using a much larger initial displacement value. The final displacement value (50 m) was calculated to be almost twice the reported 75% of the smallest displacement errors, which was more aggressive in eliminating outliers while none of the legitimate points were rejected. Since the external dataset used for alignment was in the same coordinate system as the Pleiades derived 3D point cloud, but in orthometric heights, as a by-product, the aligned result elevation values were comparable with orthometric heights with the added benefit of eliminating some of the vertical elevation bias.
Subsequently, all the 3D derived point clouds (four unaligned ellipsoid models), the aligned 3D point cloud, and the orthometric height converted 3D point cloud using the Institute Geographic Military (IGM) transformation equations [47] were gridded according to the procedure in [48] to obtain the final DSMs with a spatial resolution of 2 m.

Results
Using two pairs of stereo images and applying the LSW and MGM feature matching algorithms to each, four different DSMs in ellipsoid heights were derived from the Pleiades data with a spatial resolution of 2 m in the WGS84 UTM 33N projection (Figure 4). as the Pleiades derived 3D point cloud, but in orthometric heights, as a by-product, the aligned result elevation values were comparable with orthometric heights with the added benefit of eliminating some of the vertical elevation bias. Subsequently, all the 3D derived point clouds (four unaligned ellipsoid models), the aligned 3D point cloud, and the orthometric height converted 3D point cloud using the Institute Geographic Military (IGM) transformation equations [47] were gridded according to the procedure in [48] to obtain the final DSMs with a spatial resolution of 2 m.

Results
Using two pairs of stereo images and applying the LSW and MGM feature matching algorithms to each, four different DSMs in ellipsoid heights were derived from the Pleiades data with a spatial resolution of 2 m in the WGS84 UTM 33N projection (Figure 4).  Both Pleiades DSMs replicate the natural and anthropogenic features present in the scene (e.g., boundaries and channels of lava flows, vegetation on scoria cones and roads), but with different degrees of detail. All of the elements reproduced in the MGM DSM appeared more detailed than the same identified elements in the LSW DSM, being more similar to those observed in the lidar derived as the Pleiades derived 3D point cloud, but in orthometric heights, as a by-product, the aligned result elevation values were comparable with orthometric heights with the added benefit of eliminating some of the vertical elevation bias. Subsequently, all the 3D derived point clouds (four unaligned ellipsoid models), the aligned 3D point cloud, and the orthometric height converted 3D point cloud using the Institute Geographic Military (IGM) transformation equations [47] were gridded according to the procedure in [48] to obtain the final DSMs with a spatial resolution of 2 m.

Results
Using two pairs of stereo images and applying the LSW and MGM feature matching algorithms to each, four different DSMs in ellipsoid heights were derived from the Pleiades data with a spatial resolution of 2 m in the WGS84 UTM 33N projection (Figure 4).  Both Pleiades DSMs replicate the natural and anthropogenic features present in the scene (e.g., boundaries and channels of lava flows, vegetation on scoria cones and roads), but with different degrees of detail. All of the elements reproduced in the MGM DSM appeared more detailed than the same identified elements in the LSW DSM, being more similar to those observed in the lidar derived Both Pleiades DSMs replicate the natural and anthropogenic features present in the scene (e.g., boundaries and channels of lava flows, vegetation on scoria cones and roads), but with different degrees of detail. All of the elements reproduced in the MGM DSM appeared more detailed than the same identified elements in the LSW DSM, being more similar to those observed in the lidar derived Remote Sens. 2019, 11, 2983 9 of 17 DSM. Additionally, the 3D point cloud generated by the LSW algorithm has a larger dispersion and more data gaps over which the gridding procedure has to interpolate than the MGM 3D point cloud.
In order to assess the most accurate model, an elevation validation was computed for both ellipsoid and orthometric height DSMs using the GCPs belonging to the geodetic monitoring network of Mt. Etna [37]. From these, 23 GCPs were used to validate the 1&2 DSMs and 41 GCPs were used to validate the 2&3 DSMs, respectively. Table 2 summarizes, for each DSM, the main statistic parameters calculated for validation, height differences (errors) between the estimated data (Pleiades DSM) and the measured data (GCPs from the geodetic monitoring network of Mt. Etna). The distribution of errors calculated for each DSM is described by the split violin plots in Figure 6. A violin plot is, to a large extent, a box plot that incorporates information about the underlying distribution of the data, adding a rotated kernel density plot on each side [49,50]. In comparison, a box plot shows only the summary statistics such as mean/median and interquartile ranges, whereas the violin plot shows the full distribution of the data. A split violin plot shows two different groupings for a certain category, in our case, the different matching algorithms for the same pair of Pleiades images (LSW and MGM, Figure 6a) or different Pleiades pair images for the same matching algorithm (1&2 and 2&3, Figure 6b).
The position errors between the estimated values and the measured values were described using the mean error (ME), mean absolute error (MAE), RMSE, and standard deviation (SD). The ME indicates if there is any bias in the data, the MAE represents the average magnitude of the errors without considering their direction, and the RMSE penalizes large errors more than the others, since RMSE increases with the variance of the frequency distribution of error values. The SD indicates what could be the minimum RMSE if the bias in the data is eliminated. The vertical RMSE ranged from 6.68 to 4.80 m ( Table 2) and the SD from 1.52 to 0.88 m, indicating that only the 2&3 MGM model had the capacity to achieve a submeter RMSE, if the bias is eliminated. Only the 2&3 DSM obtained using the LSW matching algorithm was converted from ellipsoid to orthometric heights using the Institute Geographic Military (IGM) transformation equations [47]. Usually in Italy, the transformation from ellipsoid to orthometric heights is not free of charge. Table 2 shows no significant differences between the validation in ellipsoid and orthometric heights, so the other three DSM models were not converted into orthometric heights. Specifically, Table 2 highlights that the 2&3 DSM results were significantly better than the 1&2 DSM results. These validated Pleiades DSMs were obtained using only information from the satellite imagery and its metadata without any external GCPs or previous external DEM. In this case, their vertical RMSE are comparable with results obtained from RADAR (TanDEM -X 2012) and raster cartography (ATLAS) (see Table 1).   Since only the 2&3 MGM 3D point cloud had the smallest SD (0.88 m, see Table 2), and hence the potential of submeter RMSE, the least amount of data gaps, and produced the best visually pleasing DSM ( Figure 5), we decided that only this model was to be aligned with the pre-existing high accuracy lidar DEM to save computation time. In order to validate the final aligned 2&3 Pleiades MGM DSM, 41 GCP points were selected. The vertical validation statistics and error distributions are presented in Table 3 and Figure 7. The vertical RMSE decreased from 4.80 m for the not-aligned model to 0.78 m for the aligned model.  In addition, 60 points were randomly selected along roads and parking lots in the vicinity of Rifugio Sapienza, located to the south of the main summit crater (Figure 8)   In addition, 60 points were randomly selected along roads and parking lots in the vicinity of Rifugio Sapienza, located to the south of the main summit crater (Figure 8) Table 4. From these road points, over 76% have vertical alignment differences between −0.5 to +0.5 m, and only three points out of 60 had differences outside the +/− 1 m range. The low alignment difference RMSE (48 cm) and small negative bias (15 cm) indicate that the alignment procedure is adequate.

Statistics, Meters Alignment Errors on Roads-60 Random Points
Mean Error −0.15 The final Pleiades DSM in orthometric height derived from the 2&3 stereo pair matched with the MGM algorithm and aligned to lidar DTM is shown in Figure 9 as a 3D perspective image and covers an extent of 400 km 2 with elevation ranges from 382 m to 3327 m. After the satellite images are available, it can take less than 48 hours to produce the final product.  Table 4. From these road points, over 76% have vertical alignment differences between −0.5 to +0.5 m, and only three points out of 60 had differences outside the +/− 1 m range. The low alignment difference RMSE (48 cm) and small negative bias (15 cm) indicate that the alignment procedure is adequate. The final Pleiades DSM in orthometric height derived from the 2&3 stereo pair matched with the MGM algorithm and aligned to lidar DTM is shown in Figure 9 as a 3D perspective image and covers an extent of 400 km 2 with elevation ranges from 382 m to 3327 m. After the satellite images are available, it can take less than 48 hours to produce the final product. Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 19 Figure 9. Pleiades digital surface model (DSM) 2&3 pair in the 3D perspective view.

Discussion
In this section, we discuss three key points: (i) the suitability of NASA ASP software to accurately process Pleiades satellite imagery for high relief areas; (ii) the selection of the feature matching algorithm to create a DSM as accurate as possible in ellipsoid heights; and (iii) a more suitable procedure to obtain a DSM in orthometric heights with high vertical accuracy and eliminate biases in the data.
The first point was realized by using different combinations of Pleiades satellite imagery over the summit of Mt. Etna with an elevation range in excess of 2900 m, with different correlation algorithms. The results were validated using high accuracy GCPs selected from the geodetic monitoring network of Mt. Etna.
The second point was achieved by testing two different algorithms (LSW and MGM) in the same area, but using different pairs of satellite imagery. The Pleiades DSM in ellipsoid elevations (not aligned to a lidar DTM) were positively biased to 4.5 to 6.5 m, with an error RMSE between 4.6 to 6.8 m ( Figure 6 and Table 2). The violin plots displayed a noticeable bi-modality for pair 1&2 for both algorithms, and for pair 2&3 with the LSW processing algorithm, or perhaps a weak tri-modality. The last lower elevation mode only has one or two GCP points in that population, so we can consider the lower elevation mode as under-sampled. The higher elevation mode contained approximately 45 to 55% of the total GCP points used in the validation.
The error distribution for the Pleiades pair 2&3 with the MGM processing algorithm displayed a moderate tri-modality, being the single model with a SD below 1 m, and consequently, the single model that theoretically can attain a submeter vertical RMSE if the bias is removed ( Figure 6 and Table 2). There is no statistically significant difference (Welch two-sample t-test, p-value >> 0.05) between the results of the two algorithms used for each Pleiades pair (Figure 6a), although the MGM algorithm generally creates a smoother and more detailed surface, having less data gaps to interpolate over. However, there was a statistically significant difference between the results produced by the two pairs of Pleiades satellite imagery (1&2, and 2&3), respectively, for the same processing algorithm (Welch two-sample t-test, p-value << 0.05) (Figure 6b). This dissimilarity could be attributed in part to the different incidence angle of pair 1&2 (12 and 16 degrees, respectively) versus pair 2&3 (16 degrees for both images).
The 2&3 LSW DSM had the smallest vertical RMSE (4.66 m), but its SD at 1.18 m indicates that even if the bias is removed successfully, the vertical RMSE will not be submeter. The direct transformation from ellipsoid to orthometric heights using the IGM equations does not eliminate the bias in the data, although the SD decreased from 1.18 m to 1.03 m, respectively ( Table 2). The model obtained from pair 2&3 using the MGM algorithm not only had the smallest SD, suggesting the

Discussion
In this section, we discuss three key points: (i) the suitability of NASA ASP software to accurately process Pleiades satellite imagery for high relief areas; (ii) the selection of the feature matching algorithm to create a DSM as accurate as possible in ellipsoid heights; and (iii) a more suitable procedure to obtain a DSM in orthometric heights with high vertical accuracy and eliminate biases in the data.
The first point was realized by using different combinations of Pleiades satellite imagery over the summit of Mt. Etna with an elevation range in excess of 2900 m, with different correlation algorithms. The results were validated using high accuracy GCPs selected from the geodetic monitoring network of Mt. Etna.
The second point was achieved by testing two different algorithms (LSW and MGM) in the same area, but using different pairs of satellite imagery. The Pleiades DSM in ellipsoid elevations (not aligned to a lidar DTM) were positively biased to 4.5 to 6.5 m, with an error RMSE between 4.6 to 6.8 m ( Figure 6 and Table 2). The violin plots displayed a noticeable bi-modality for pair 1&2 for both algorithms, and for pair 2&3 with the LSW processing algorithm, or perhaps a weak tri-modality. The last lower elevation mode only has one or two GCP points in that population, so we can consider the lower elevation mode as under-sampled. The higher elevation mode contained approximately 45 to 55% of the total GCP points used in the validation.
The error distribution for the Pleiades pair 2&3 with the MGM processing algorithm displayed a moderate tri-modality, being the single model with a SD below 1 m, and consequently, the single model that theoretically can attain a submeter vertical RMSE if the bias is removed ( Figure 6 and Table 2). There is no statistically significant difference (Welch two-sample t-test, p-value >> 0.05) between the results of the two algorithms used for each Pleiades pair (Figure 6a), although the MGM algorithm generally creates a smoother and more detailed surface, having less data gaps to interpolate over. However, there was a statistically significant difference between the results produced by the two pairs of Pleiades satellite imagery (1&2, and 2&3), respectively, for the same processing algorithm (Welch two-sample t-test, p-value << 0.05) (Figure 6b). This dissimilarity could be attributed in part to the different incidence angle of pair 1&2 (12 and 16 degrees, respectively) versus pair 2&3 (16 degrees for both images).
The 2&3 LSW DSM had the smallest vertical RMSE (4.66 m), but its SD at 1.18 m indicates that even if the bias is removed successfully, the vertical RMSE will not be submeter. The direct transformation from ellipsoid to orthometric heights using the IGM equations does not eliminate the bias in the data, although the SD decreased from 1.18 m to 1.03 m, respectively ( Table 2). The model obtained from pair 2&3 using the MGM algorithm not only had the smallest SD, suggesting the potential for the smallest vertical RMSE if the bias is eliminated, but also had one of the lowest vertical RMSE and the best visual representation of morphology (Figures 4 and 5). This implies that the MGM algorithm is more suitable to match images where the topography is very irregular and characterized by frequent changes in heights. For this reason, we chose to use the MGM algorithm in the third step of the processing workflow for alignment ( Figure 3) as being more suitable for high relief topographies characterized by strong roughness such as Mt. Etna area.
The third point is illustrated by comparing the validation statistics between the orthometric model obtained without the use of any GCPs or a previous high-accuracy DTM and the model vertically aligned to the lidar DTM. In this case, the validation vertical RMSE decreased from 4.8 m to 0.78 m for the aligned MGM 2&3 DSM model and the bias from +4.72 m to −0.52 m, respectively ( Table 3). The advantage of the alignment to a pre-existing DEM in orthometric heights is that the aligned results will have comparable orthometric heights and the bias is minimized, although not entirely eliminated. This procedure is suitable to update preexistent topographies and therefore it can be used for areas affected by frequent morphological changes such as active volcanic zones. Figure 7 shows similar data distributions for vertical validation errors and vertical alignment differences on the 41 GCP locations. There was no statistical difference (Welch two-sample t-test, p-value >> 0.05) between the aligned Pleiades DSM validation errors (Table 3, column 3) and the alignment differences on the same GCP points (Table 3, column 5). While the Pleiades DSM was negatively biased to the GCP values by approximately half a meter, the lidar DSM was positively biased at approximately 17 cm for the same GCP points (Table 3). This translates into a negative bias between the Pleiades DSM and the lidar DSM of 64 cm (Table 3) on the same 41 GCP locations.
Alternatively, orthometric heights can be achieved using the transformation equations implemented by the IGM. In Italy, that change according to the latitude and longitude of the investigated zone [47]. Converting ellipsoid heights into orthometric heights using the IGM transformation equations [47] does not achieve the same vertical accuracy obtained by the alignment procedure since it does not eliminate the vertical bias (Table 2). If the vertical bias can be removed after the orthometric conversion, the final RMSE in this case could be as low as 1.03 m for the LSW model (Table 2). This is important since we want the data to be part of a reference system that can be used to compare and conduct change analysis with, against other sets of data. The workflow we propose using NASA ASP needs to be repeatable and robust, so that results without data alignment (in both ellipsoid and orthometric heights) will have their limitations understood.

Conclusions
In this work, a submeter vertical RMSE for Mt. Etna high-resolution DSM derived from Pleiades satellite imagery using NASA ASP open source software was achieved for the first time. We tested two correlation algorithms, LSW and MGM, with two different datasets and several parameters. The results indicate that (i) for the same stereo pair, there was no statistical significance between the vertical errors obtained using the LSW or the MGM algorithms, although the MGM algorithm always produced less noisy results with less data gaps; (ii) there was a marked statistical difference between the vertical errors obtained from different image combination pairs. The difference in vertical errors is likely to be due to the incidence angle difference. The 2&3 pair (incidence angle of 16 degrees) performed better than the 1&2 pair (incidence angle of 12 and 16 degrees, respectively); and (iii) data alignment to a pre-existing high-accuracy DEM referenced to orthometric heights allowed the model to achieve a vertical RMSE below 1 m (78 cm). Most of the alignment differences (over 76% of data) measured on roads and parking lots near the Rifugio Sapienza area located south of the summit have a vertical RMSE of 48 cm.
The NASA ASP supports imagery that uses RPC camera models, and previously tested and calibrated for DigitalGlobe satellite imagery. This research investigates the NASA ASP performance in terms of absolute vertical accuracy when using Pleiades satellite imagery over terrain with large elevation ranges and roughness. Although lidar data are always preferable when available, due to high acquisition costs and longer post-processing times-sometimes up to several months-they are not always suitable for frequent and repeated surveys in areas with a dynamic morphology. In contrast, Pleiades satellite imagery has a repeat cycle of 26 days, can cover 400 km 2 in one scene, and can be processed in less than 48 hours to obtain a submeter vertical RMSE. The submeter vertical RMSE is achieved in ideal conditions such as imagery free of clouds, with equal angles of incidence, and a very accurate if potentially sparse external DEM for alignment. Considering the high spatial resolution (2 m), the vertical accuracy (78 cm RMSE) and the relative short processing times to generate the DSM, this procedure is appropriate to produce updated topographies for active volcanic areas. The derived orthometric aligned DSM model and associated metadata are available as an electronic supplement to this work.
Future work will investigate other approaches to increase the accuracy of the digital model including the use of other bands from the multispectral stereo imagery when the stereo panchromatic band is not available.