High‐resolution digital elevation model from tri‐stereo Pleiades‐1 satellite imagery for lava flow volume estimates at Fogo Volcano

Resolving changes in topography through time using accurate high‐resolution digital elevation models (DEMs) is key to understanding active volcanic processes. For the first time in a volcanic environment, we utilize very high‐resolution tri‐stereo optical imagery acquired by the Pleiades‐1 satellite constellation and generate a 1 m resolution DEM of Fogo Volcano, Cape Verde—the most active volcano in the Eastern Atlantic region. Point cloud density is increased by a factor of 6.5 compared to conventional stereo imagery, and the number of 1 m2 pixels with no height measurements is reduced by 43%. We use the DEM to quantify topographic changes associated with the 2014–2015 eruption at Fogo. Height differences between the posteruptive Pleiades‐1 DEM and the preeruptive topography from TanDEM‐X give a lava flow volume of 45.83 ± 0.02 × 106 m3, emplaced over an area of 4.8 km2 at a mean rate of 6.8 m3 s−1.


Introduction
Volcanic activity is among the fastest processes causing changes to the Earth's surface, whether new volumes are added during effusive eruptions or catastrophically removed by explosive events and triggered landslides. Precise, quantitative analyses of topographic changes associated with volcanic eruptions provide the means to infer key parameters for the assessment of hazards associated with the volcanic activity (e.g., magma discharge rate in effusive events) [Pinel et al., 2014]. The topographic approach, which constrains the changes in topography by differentiating preeruptive, coeruptive, and posteruptive digital elevation models (DEMs) [Stevens et al., 1999], can nowadays be considered the most suitable method to accurately quantify the volume of new volcanic deposits, especially when data are acquired by spaceborne Earth Observation (EO) platforms [e.g., Lu et al., 2003;Rowland et al., 2003;Bignami et al., 2013;Poland, 2014;Albino et al., 2015;Kubanek et al., 2015;Martino et al., 2015]. These data sets provide densely spaced measurements of heights at relatively high temporal frequency (hours to days, if combined) and meter level vertical accuracy, without the need for direct field measurements. The cost of spaceborne EO data is also far lower than that of dedicated airborne or ground-based campaigns (e.g., Light Detection And Ranging, LiDAR).
In this study, we evaluate the use of very high resolution (VHR) tri-stereo optical imagery from the Pleiades-1 satellite constellation for volcanological applications. The data set is 100% cloud free and was acquired after the 2014-2015 eruption at Fogo Volcano, Cape Verde. The eruption, the tenth since 1761 [Amelung and Day, 2002], lasted 78 days (23 November 2014 to 8 February 2015) and emplaced highly destructive lava flows in the Chã das Caldeiras (Figure 1) [González et al., 2015;Cappello et al., 2016;Richter et al., 2016]. We apply photogrammetric techniques to the Pleiades-1 tri-stereo data set and generate a 1 m resolution posteruptive DEM of the summit and eastern flank of Fogo.
Introduced with the Pleiades mission, the stereoscopic triplet is a novel and enhanced acquisition mode compared to that of previous optical sensors (e.g., Quickbird, IKONOS). The triplet is composed of three nearly simultaneously acquired images, one backward looking, one forward looking, plus a third near-nadir image [Gleyzes et al., 2012]. This configuration, in particular the use of the near-nadir image, allows a better retrieval of heights over terrains where the performance of classic photogrammetry with forward backward looking stereo pairs may be limited (e.g., high roughness, steep slopes, and shadows). The Pleiades-1 mission is programmed to continue for several years [Centre National d'Etudes Spatiales, 2016], and the tri-stereo (or multistereo) acquisition mode is set to become a common feature in future satellite missions. Pleiades-1 tri-stereo imagery has recently proven to be valuable in estimating submeter to meter scale elevation changes associated with natural phenomena, such as earthquakes [Zhou et al., 2015 [Berthieratelliet al., 2014]. Fogo Volcano represents an optimal location for the assessment of the use of Pleiades-1 tri-stereo imagery to study volcanic landforms and processes. The volcanic edifice is characterized by varying slopes (0°to >50°) and variable land coverage, factors that could influence the accuracy or prevent the retrieval of heights in a natural environment.
From the Pleiades-1 posteruption topography we subtract heights from a preeruptive DEM, obtained using spaceborne synthetic aperture radar (SAR) data from the TanDEM-X mission. Such differencing provides the means to estimate the volume of the 2014-2015 lava flow with an unprecedented accuracy for Fogo Volcano [e.g., Richter et al., 2016]. The volume estimate allows us to constrain the mean output rate of the 2014-2015 eruption and infer a minimum magma supply rate since the last eruption. We also precisely map the lava flow areal extent utilizing Pleiades-1 orthorectified VHR multispectral images. Finally, using SAR data acquired by the Sentinel-1a satellite, we apply SAR interferometry (e.g., InSAR) to measure the lava flow subsidence due to cooling and contraction in the months after its emplacement and compare this to the measured lava flow thickness.

Data Processing and Analysis
We processed the preeruptive TanDEM-X DEM and the Sentinel-1 InSAR data using standard methods and algorithms described in the supporting information (Texts S1 and S2).

Pleiades-1 Tri-Stereo Optical Imagery and Processing
The Pleaides-1 system consists of a constellation of two satellites for VHR panchromatic (PA) and multispectral (  other, which allows a minimum revisit time of 24 h. The PA and XS images are acquired simultaneously at a nominal resolution of 0.5 m and 2 m, respectively [de Lussy et al., 2012]. The Pleiades system is the first of its kind capable of acquiring three or more nearly synchronous images of the same area with a stereo angle varying between~6°and~28°. After the end of the 2014-2015 eruption at Fogo Volcano, we tasked the acquisition of Pleiades-1 tri-stereo optical imagery to study the new topography. PHR 1A acquired a set of images on 20 June 2015, with a total areal coverage of~100 km 2 (red outline in Figure 1). Along-track incidence angles of the three images are À11.3°, 1.0°, and 11.2°for the forward (F), near-nadir (N), and backward (B) viewing geometries, respectively, while the across-track angle varies between À2°and 3°.
We processed the Pleiades-1 tri-stereo data using the ERDAS IMAGINE Photogrammetry toolbox. We tested different combinations of parameters for the photogrammetric processing and opted for those that maximized the number of matched points. We first identified~120 tie points over the three PA images to calculate shift parameters between images at each tie point, obtaining an overall root-mean-square error of less than 0.1 pixels (<0.05 m). IMAGINE Photogrammetry toolbox implements a hierarchical cross-correlation matching algorithm for the automatic DEM extraction. An image pyramid was created with nine levels of detail. After a rigid shift of the images, we computed the normalized cross-correlation index with a window size of 9 × 9 pixels, a minimum correlation coefficient of 0.2 for the highest pyramid level and a maximum coefficient of 0.7 for the last pyramid. The photogrammetric algorithm calculates point clouds (X, Y, Z) from each image pair (F-B, F-N, and N-B), for the tri-stereo triplet (F-N-B), and from the merging of the four data sets. The merged point cloud was then filtered using a median filter with a cell size of 1 m and gridded to a spatial resolution of 1 × 1 m 2 pixel size using a continuous curvature spline in tension (tension factor 0.75 [Smith and Wessel, 1990]).

Pleiades-1 Point Clouds
The tri-stereo triplet (F-N-B) provides the largest point cloud (~110.5 million points). Stereo pairs using the near-nadir image (N) produce point clouds with~96.2 (N-B) and~91.8 (F-N) million points, corresponding to 13% and 18% fewer points, respectively, than the tri-stereo combination. Finally, the classic stereo-pair (F-B) produces a point cloud with~54.3 million points, 51% less than the tri-stereo triplet. When averaged over the entire image, we obtained mean point densities of 1.11, 0.96, 0.91, and 0.54 points/m 2 for the F-N-B, N-B, F-N, and F-B combinations, respectively. By merging the four point clouds, we obtain a data set with~352.2 million points and a mean point density~6.5 times larger (3.52 points/m 2 ) than that from the classic F-B stereo pair.
Assuming that the points in a point cloud are independent, the higher their density, the more robust the height estimate for a given pixel in the derived DEM. Density is variable; therefore, the mean density is not representative of the accuracy of the estimated heights. In Figures 2a-2d we show a comparison of the density of points for each combination of images, measured across a 10 × 10 m 2 regular grid. It can be clearly seen that point cloud density is not homogenous across the imaged area. To better understand this distribution, we compare the point cloud density maps with the surface slope ( Figure 2e) and the geological map of the same area ( Figure 2f). The surface slope does not directly correlate with point cloud density (correlation coefficient: <0.1), and variable distribution can be observed in the flat areas as well as on the steep slopes. On the other hand, there is a correlation between point cloud density and ground cover. In decreasing order of density are areas covered by lava flows, the summit of Pico do Fogo where lava and scoria are present, cinder cones, and, lastly, ash-covered regions (e.g., the mid/lower portion of Pico do Fogo). These last two land types are characterized by loose, fine, granular material at repose angle, which appear largely homogeneous and featureless to an optical sensor, making cross correlation less successful. In this case the surface slope may have an indirect effect by controlling the stability and remobilization of the loose material, resulting in featureless surfaces. The use of the triplet (F-N-B) or the merged data set, however, reduces the number of 1 × 1 m 2 pixels with no measurements by 27% and 43%, respectively, when compared to classic stereo (F-B).

DEM Horizontal and Vertical Accuracy
At the time of writing, no other high-accuracy DEMs of Fogo are publicly available to assess the accuracy (closeness to the true value) of the heights and their horizontal position in the Pleiades-1 DEM. We therefore used differential GPS solutions for the X, Y, Z position of 19 ground control points (GCPs, Figure 1)  and their clearly visible positions in a pansharpened image, formed using the near-nadir PA/XS images and orthorectified using the DEM. The mean offsets are À7.6 m and À1.3 m in the east and north directions, respectively (standard deviations, 0.4 m and 0.3 m), lower than the nominal absolute location accuracy of PHR 1A imagery (8.5 m [Oh and Lee, 2014]). The estimated horizontal shift was then applied to the Pleiades-1 DEM, and the vertical accuracy was estimated at the 19 GCPs. The same assessment was also performed, for completeness, on the 5 m resolution preeruptive TanDEM-X DEM and other two lower resolution DEMs of Fogo (the 1 arc second resolution Shuttle Radar Topographic Mission, SRTM, and the 30 m resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Map, GDEM) (see Figure S1 in supporting information). The standard deviation of heights in the Pleiades-1 DEM when compared to GPS is 0.51 m, similar to that reported by Zhou et al. [2015] from the comparison of a Pleiades-1 DEM with an airborne LiDAR DEM. The other three data sets, TanDEM-X, SRTM, and ASTER-GDEM have standard deviations of 1.12 m, 3.64 m, and 5.74 m, respectively (see Table S3 in the supporting information). Part of the contribution to all standard deviations comes from errors in the GPS data, but the analysis demonstrates that the higher resolution DEMs are more accurate, and that Pleiades 1 has a standard deviation lower by a factor of at least 2 than TanDEM-X and at least 7 than SRTM and ASTER.
To accurately constrain the 2014-2015 lava flow volume, the relative difference in heights between the preeruptive TanDEM-X and the posteruptive Pleiades-1 DEMs must be adjusted [e.g., Poland, 2014]. To accomplish this, we selected three 600 × 600 m 2 control areas located near the lava flow that should not have had any change in topography between the acquisitions of the two data sets (Figure 3a). We oversampled the TanDEM-X 5 m resolution data set to the resolution of the Pleiades-1 DEM (1 m). The mean values in the control areas are À3.03, À3.05, and À3.01 m, while standard deviations are 0.76, 0.73, and 0.48 m. We therefore added the weighted mean value of 3.03 m to the Pleiades-1 DEM before differencing the two DEMs. The mean values in the control areas after the adjustment are 0.00, À0.02, and 0.01 m, and the standard deviation, σ ΔZ , of all three areas is 0.67 m (Figure 3b). Height differences in the control areas have a Gaussian distribution, but spatial correlation exists between neighboring pixels. We therefore estimated the covariance as function of distance and approximated it as a 1-D exponential function: with a = 0.2703, b = À0.2145, and r is the distance between pixels in meters. With this function, the covariance reduces to approximately zero at a distance of~15 m.

Area and Volume Estimates
The extent of the 2014-2015 lava flow was calculated by digitizing the flow perimeter on orthorectified VHR XS images. Using geographic information system software (QGIS, http://www.qgis.org) we calculated that, at the end of the eruption, the flow covered a surface of 4.8 km 2 .
A mask derived from the lava flow perimeter was applied to the difference between the preeruptive and posteruptive DEMs to constrain the lava flow volume (Figure 3a). The sum of the values of height changes within the lava flow perimeter, multiplied by the area of a single pixel (1 m 2 ) provides the estimate of the total volume [e.g., Favalli et al., 2010]. We estimated the error variance of the volume calculation through propagation of errors where A is the area of a pixel and C ij is the covariance between the height error at pixel i and pixel j. In practice we assumed that the covariance between all pixels >15 m apart was zero; using equation (2), we calculated the covariance, C r , between a single pixel and the n pixels up to 15 m from it. Assuming all pixels are >15 m from the lava edge gives where is N the total number of pixels within the flow perimeter. The contribution to the total variance from pixels <15 m from the edge is thus slightly exaggerated, but this effect is small.
Based on the method and the assumptions described above, we estimate that the 2014-2015 lava flow at Fogo Volcano had a total volume of 45.83 ± 0.02 × 10 6 m 3 as of 20 June 2015. Note that this ignores any subsidence of the substrate due to elastic or inelastic processes. A maximum thickness of >55 m was emplaced on the NW side of the eruptive fissure, where a scoria/cinder cone formed during the eruption. The lava flow reached its maximum thickness~650 m SW of the emission point. A thickness of~26 m is also measured more than 2500 m away from the fissure, in an area of preeruptive topographic low, W of the town of Portela (Figure 3a, profile A-A′ in Figure 3c) [e.g., Richter et al., 2016].

Comparison Between Lava Flow Thickness and Subsidence
A comparison between the lava flow thickness estimate and the cumulative lava flow subsidence from cooling and contraction, measured using InSAR data spanning~9 months after the end of the eruption (May 2015 to March 2016), shows a clear spatial correlation (Figures 4a and 4b). In fact, the maximum subsidence is recorded in those areas where lava flow thickness is also maximum. Other smaller regions, however, show large cumulative subsidence and deviate from the overall linear correlation between thickness and subsidence (Figure 4c). The larger subsidence in these areas is most likely related to the highly compactable nature of the substrate onto which the lava flow was emplaced (e.g., Pico do Fogo ash under the southern branch and talus at the terminal part of the western branch of the flow or west of the road in the northern branch; Figures 4b and 4c). The correlation between lava flow thickness and subsidence likely explains the relationship between thickness and phase decorrelation often observed in InSAR imagery [e.g., Dietterich et al., 2012;Albino et al., 2015].

Discussion
The use of VHR Pleiades-1 imagery and radar data from TanDEM-X has allowed an accurate estimate of the areal extent and volume of the 2014-2015 lava flow at Fogo Volcano. The level of accuracy of our estimate is comparable, or higher, to that achieved by the most recent studies at other volcanoes that used data from spaceborne sensors [e.g., Poland, 2014;Albino et al., 2015;Kubanek et al., 2015]. The areal extent (4.8 km 2 ) is 0.6 km 2 smaller than that calculated by Cappello et al. [2016] using Landsat-8 OLI imagery but matches very well that calculated by Richter et al. [2016] using coherence maps of TerraSAR-X InSAR pairs (4.85 km 2 ). Similarly, our volume estimate (45.83 ± 0.02 × 10 6 m 3 ) agrees very well with estimates by Richter et al.
The erupted volume during the 2014-2015 eruption is very similar to that emplaced in the previous eruption at Fogo in 1995 (~46 × 10 6 m 3 [Amelung and Day, 2002]) and approximately half of that erupted in 1951 (~110 × 10 6 m 3 [Cappello et al., 2016]). The eruption lasted 78 days, and by dividing the total erupted volume by the duration of the eruption we can infer the bulk mean output rate (MOR BULK ) [e.g., Harris et al., 2007] to be equal to~588,000 m 3 d À1 or~6.8 m 3 s À1 . The MOR BULK , however, includes vesicles and gaps between lava blocks. The dense rock equivalent MOR (MOR DRE ) can be determined if the vesicularity is known, but this parameter is not constrained in the case of the most recent lava flow at Fogo. We consider therefore an average vesicularity for basaltic aa lava flows of 25% [e.g., Wolfe et al., 1987;Poland, 2014] and adjust the total volume to be~34.37 ± 0.02 × 10 6 m 3 . The MOR DRE is therefore~441,000 m 3 d À1 or~5.1 m 3 s À1 . This value is the mean discharge rate throughout the entire eruption. HOTSAT thermal sensors [Cappello et al., 2016] and SO 2 emissions [Barrancos et al., 2015] indicate that the discharge rate was much higher (>25 m 3 s À1 ) during the initial phase of the eruption.
By adding the erupted volume to that of the intrusion that fed the eruption we can also estimate the minimum magma supply rate since the last eruption, on the assumption that all the magma supplied since the previous eruption in 1995 was either intruded or erupted in 2014-2015. González et al. [2015] estimated that the total volume of the intrusion was~3 ± 2 × 10 6 m 3 . Adding this to the DRE lava flow volume gives 37.4 ± 2 × 10 6 m 3 of magma. Although magma accumulation at depths greater than 15 km may go undetected by geodetic measurements [Amelung and Day, 2002], no significant deformation related to magma accumulation/withdrawal from shallow reservoirs has been recorded during, or between, the two most recent eruptions [González et al., 2015]. By dividing the supplied volume by 19.5 years (time interval between the two eruptions) we obtain a minimum magma supply rate of 1.9 ± 0.1 × 10 6 m 3 yr À1 or 0.06 m 3 s À1 . Although speculative, this number is similar to that inferred by Amelung and Day [2002] over 337 years (>1.7 × 10 6 m 3 yr À1 ). This rate of magma supply is 2 orders of magnitude lower than that of Kilauea, Hawaii (3-4 m 3 s À1 [Poland, 2014]) and~10% of that of all Galápagos volcanoes (~0.6 m 3 s À1 [Bagnardi, 2014]), both of which are also related to hot spot ocean island volcanism.

Conclusions
We demonstrate that Pleiades-1 tri-stereo imagery can be successfully used to derive high-resolution DEMs with submeter vertical accuracy to study volcanic processes. Its performance is best over lava flows and lower in featureless regions covered with loose granular deposits (e.g., ash and cinder). When compared with the classic stereo approach, the use of tri-stereo imagery highly enhances the ability of photogrammetric techniques to estimate heights through increasing the point cloud density-by a factor of 6.5 in the example of Fogo Volcano-and by reducing the number of 1 m 2 pixels with no measurements (by 43% at Fogo).

10.1002/2016GL069457
Tri-stereo (or multistereo imagery) is likely to become a common feature in future optical satellite missions, and its application to volcanic studies will provide the means to study topographic changes and volcanic features at high resolution.
By subtracting from the Pleiades-1 posteruptive DEM, the preeruptive topography measured using TanDEM-X radar data, we estimated the volume of the 2014-2015 lava flow at Fogo Volcano to be 45.83 ± 0.02 × 10 6 m 3 (34.37 ± 0.04 m 3 , when converted to DRE), emplaced over an area of 4.8 km 2 at a mean rate of 6.8 m 3 s À1 (MOR DRE : 5.1 m s À1 ). Assuming no long-term crustal magma storage, this implies a minimum magma supply rate to Fogo of 0.06 m 3 s À1 between 1995 and 2014, which, although speculative, is much lower than that measured at other ocean island volcanoes (e.g., Hawaii and Galápagos).