On the Use of Tri-Stereo Pleiades Images for the Morphometric Measurement of Dolines in the Basaltic Plateau of Azrou (Middle Atlas, Morocco)

: Hundreds of large and deep collapse dolines dot the surface of the Quaternary basaltic plateau of Azrou, in the Middle Atlas of Morocco. In the absence of detailed topographic maps, the morphometric study of such a large number of features requires the use of remote sensing techniques. We present the processing, extraction, and validation of depth measurements of 89 dolines using tri-stereo Pleiades images acquired in 2018–2019 (the European Space Agency (ESA) © CNES 2018, distributed by Airbus DS). Satellite image-derived DEMs were ﬁeld-veriﬁed using traditional mapping techniques, which showed a very good agreement between ﬁeld and remote sensing measures. The high resolution of these tri-stereo images allowed to automatically generate accurate morphometric datasets not only regarding the planimetric parameters of the dolines (diameters, contours, orientation of long axes), but also for what concerns their depth and altimetric proﬁles. Our study demonstrates the potential of using these types of images on rugged morphologies and for the measurement of steep depressions, where traditional remote sensing techniques may be hindered by shadow zones and blind portions. Tri-stereo images might also be suitable for the measurement of deep and steep depressions (skylights and collapses) on Martian and Lunar lava ﬂows, suitable targets for future planetary cave exploration. Pleiades-derived DEM, we carried out a ground-truthing on a subset of 89 different dolines (out of a total of 357 identiﬁed using low-resolution satellite data including LANDSAT 7 ETM+, ASTER, and ESA Sentinel 2). We found an overall good correlation ( p -value < 0.05 and correlation coefﬁcient >0.9) between the DEM-derived depths and those measured with telemeter in the ﬁeld. The largest differences were obtained on deep funnel-shaped dolines and on shallow bowl-shaped dolines, where disturbance is created by tall trees, human-built structures (dry walls), boulders, or a combination of these features. Our study shows that very-high-resolution satellite images such as Pleiades (0.5 m ground resolution), are extremely useful in obtaining precise planimetric measurements on large landforms such as dolines. DEMs extracted from Pleiades stereo images deliver good results on the vertical scale if doline ﬂoors are not masked by vegetation or disturbed by large debris (boulders or dry walls), with mean depth errors of 2 ± 2 m. These results will allow for extending the information extracted in the field to the dolines present in inaccessible areas of the plateau, allowing to obtain more information to complete their morphological profiling and propose a model on their unknown origin. The results of our work demonstrate that tri-stereo satellite images are suitable for the morphometric assessment of shadow-prone negative and positive landforms, and might also be used for the analysis of skylights and collapses on lava ﬂows on the Moon and on Mars, where vegetation or human-built structures are completely lacking, as preparatory surveys for future robotic and human explorations of extraterrestrial caves


Introduction
Geological landscapes and structures can be easily studied starting from high spatial resolution Digital Elevation Models (DEMs) and a set of useful tools able to extract specific landforms and quantify geomorphological parameters [1][2][3]. Moreover, topographical maps of sufficient scale and precision are often not available for many regions in the world. In these areas, the combination of remote sensing data from different sources into Geographic Information Systems (GIS) allows to obtain new information and integrate them with field-based data.
The detection of karst morphologies from DEMs processing is a well-known methodology, largely applied for geomorphological studies, risk and geohazard assessment, and land use management [4][5][6][7][8][9][10]. Studies based on the comparison between DEM and fieldbased data have demonstrated that the use of digital images with an appropriate spatial resolution generally provides reliable results in the morphometric characterization of sinkholes [11,12].
The elevation model extracted from optical satellite images, however, must guarantee a quality standard. Accurate elevation and morphological details are critical points for geomorphological mapping purposes if satellite data are to be used [2]. Terrain roughness, ground resolution, base-to-height ratio, and image contrast affect the accuracy of the geomorphological mapping purposes if satellite data are to be used [2]. Terrain roughness, ground resolution, base-to-height ratio, and image contrast affect the accuracy of the produced DEMs from optical space images [12,13]. In the last years, several very-high resolution satellites were launched (e.g., QuickBird, OrbView, GeoEye-1, and the series of WorldView platforms, as well as Pleiades), with the capacity of providing images with a ground resolution of up to 0.3 m.
Conventional techniques for DEM creation, based on terrestrial surveying and aerial images and improved using GPS, have been integrated by the production of elevation models created with satellite stereoscopic images [14][15][16]. Optical and Synthetic Aperture Radar (SAR) are consolidated methods to obtain DEMs with different resolution and accuracy at the local, regional, and global scales. Among them, the stereoscopic triplet can be considered an improvement of the simple stereo-couple, allowing for a better retrieval of altitudes over impervious terrains, steep slopes, and shadows [17].
The Pleiades mission was the first satellite system that introduced the stereoscopic triplet, adopting a new acquisition mode, based on the nearly simultaneous acquisition of three images, one backward looking (B), one forward looking (F), plus a third near-nadir image (N) [18]. The Pleiades constellation consists of two satellites: 1A spacecraft was launched on 16 December 2011, and Pleiades 1B in late 2012. Panchromatic and multispectral data are acquired simultaneously at a nominal resolution of 0.5 m and 2 m. The products are mainly used to produce 50-cm ortho-imagery, and the acquisition of synchronous images of the same area with different looking angle, allows for obtaining very high-resolution DEMs. In particular, the use of a near-nadir (N) image improves the signal of the classical forward (F) and backward (B) looking stereo pairs, allowing a better retrieval of heights over terrains with high roughness and steep slopes characterized by large areas with shadows ( Figure 1). DEMs extracted from these data proved to be well-suited for the analyses of elevation in different geomorphological environments such as mountain glaciers, volcanic ranges, and in urban areas [13,19,20]. Moreover, the use of these data is focused on the multitemporal analysis of DTMs based on multiple acquisitions with the possibility to extract surface and volume changes [13,17]. Stereoscopic cover capabilities over the studied surfaces with deep morphological depressions: on the left, the red region represents the area that cannot be covered by the standard stereo images; on the right, the nadir image allows to complete the visualization of the steep walls and the floor of the depression (in green).
The aim of this study is to verify the use of DEMs extracted from Pleiades tri-stereo optical data for the evaluation of morphometric parameters in the Middle Atlas in Morocco, where high resolution topographic maps are not available. In particular, this study addresses the problem of evaluating the depth of sinkhole-like depressions dotting the basaltic plateau of Azrou. Several of these "dolines" are deep and have sub-vertical walls, or are vegetated by tall Atlas cedars, and are therefore affected by shadow problems. Stereoscopic cover capabilities over the studied surfaces with deep morphological depressions: on the left, the red region represents the area that cannot be covered by the standard stereo images; on the right, the nadir image allows to complete the visualization of the steep walls and the floor of the depression (in green).
The aim of this study is to verify the use of DEMs extracted from Pleiades tristereo optical data for the evaluation of morphometric parameters in the Middle Atlas in Morocco, where high resolution topographic maps are not available. In particular, this study addresses the problem of evaluating the depth of sinkhole-like depressions dotting the basaltic plateau of Azrou. Several of these "dolines" are deep and have sub-vertical walls, or are vegetated by tall Atlas cedars, and are therefore affected by shadow problems.
In this study, Pleiades tri-stereo images were used for the first time to generate DEMs and to extract morphometric parameters of these steep and deep depressions with different degrees of disturbances (shadow zones, tall trees, human-built structures, etc.). A DEM with a resolution of 5 m/pixel has been extracted from Pleiades tri-stereo imagery to obtain morphometric parameters of the doline fields in the Azrou Plateau (Table 1). Whereas planimetric measurements can be obtained at very-high resolution, less is known about the accuracy of vertical measurements, and especially the depth of such steep and rough depressions. The depth measures of the sinkholes calculated from the Pleiadesderived DEM have been compared with the same parameter collected in the field for a set of these landforms (n = 89). The morphometric measurements extracted in the GIS environment were thus used to verify the vertical accuracy of the remote-sensing derived data, and especially to assess if vertical measurements (i.e., depth of dolines) are reliable. Therefore, this study proposes a possible remote sensing workflow for the use of Pleiades-derived DEMs and the morphometric study of steep and rugged shadow-prone terrains.

Study Area
The Quaternary Azrou volcanic plain in the Middle Atlas of Morocco covers more than 400 km 2 [21], and well-preserved cones and calderas are still clearly visible ( Figure 2) [22]. Martin surveyed over 100 large collapses in the 70-80s on the basalt plain using 1:100,000 maps and some aerial photographs [23]. Since the basalts cover the Jurassic limestones, it is not well known if these voids are karstic in origin, or if they are parts of lava tube systems. Williams defined them as "caprock collapse sinkholes" in his chapter on dolines in the Encyclopedia of Caves and Karst [24], thus suggesting the presence of karstic voids underneath. However, some small lava tube segments are known in this area [23], and some of the collapses are clearly aligned, or located at the center of the lava flows. This plateau thus offers the unique opportunity of studying collapses formed by the supposed presence of lava tubes and of karst voids (ancient doline fields or collapsed caves). Recent studies have shown a good agreement between remote-sensing derived lineaments and karst features [25,26]. These studies used LANDSAT 7 ETM+ (spatial resolution of 30 m, 15 m in panchromatic), Sentinel-2 (spatial resolution greater than 10 m), and ASTER GDEM data (spatial resolution of 30 m). The possible presence of well-developed karst systems, covered by basalts but reactivated after the emplacement of these volcanic rocks, should be taken into account when considering the vulnerability of this significant regional aquifer, with important karst springs used for drinking water purposes (e.g., Oum Er-Rbia) [26,27]. However, currently, none of the above-mentioned studies in the area used high resolution images (e.g., Pleiades).
TER GDEM data (spatial resolution of 30 m). The possible presence of well-developed karst systems, covered by basalts but reactivated after the emplacement of these volcanic rocks, should be taken into account when considering the vulnerability of this significant regional aquifer, with important karst springs used for drinking water purposes (e.g., Oum Er-Rbia) [26,27]. However, currently, none of the above-mentioned studies in the area used high resolution images (e.g., Pleiades).

Satellite Imagery
During the first steps of this study, the photointerpretation of low-resolution satellite data (LANDSAT 7 ETM+, ASTER, and ESA Sentinel 2) allowed us to map 357 possible sinkholes in the Azrou plateau ( Figure 3). The low resolution of these satellite images, however, did not allow us to extract detailed morphometric information on these collapses. In this study, we tested the Pleiades tri-stereo product (ground resolution of panchromatic data 0.5 m), provided by Airbus Defence and Space, to produce a very highresolution DEM for a detailed morphological analysis.

Satellite Imagery
During the first steps of this study, the photointerpretation of low-resolution satellite data (LANDSAT 7 ETM+, ASTER, and ESA Sentinel 2) allowed us to map 357 possible sinkholes in the Azrou plateau ( Figure 3). The low resolution of these satellite images, however, did not allow us to extract detailed morphometric information on these collapses.
In this study, we tested the Pleiades tri-stereo product (ground resolution of panchromatic data 0.5 m), provided by Airbus Defence and Space, to produce a very high-resolution DEM for a detailed morphological analysis.
As Pleiades tri-stereo data of the area of interest were not available, we planned an acquisition on-demand of a dedicated tri-stereo product. In order to use the high-resolution imagery for the photointerpretation of geomorphological features, ortho-images completed the list of data used in this study (Table 2)  As Pleiades tri-stereo data of the area of interest were not available, we planned an acquisition on-demand of a dedicated tri-stereo product. In order to use the high-resolution imagery for the photointerpretation of geomorphological features, ortho-images completed the list of data used in this study (Table 2). Data were acquired in 2018 and 2019 in the framework of an ESA-funded project (European Space Agency © CNES 2018, distributed by Airbus DS).  Pleiades products were delivered in DIMAP V2 format. This file format provides, in addition to the Image file (JPEG 2000 or GeoTIFF), the RPCs (Rational Polynomial Coefficients), allowing to extract the DEM, and to apply the orthorectification and geometric processing. Moreover, a KMZ file for the localization in Google Earth environment, and the masks on data quality and cloud cover complete the data.

DEM Extraction
The IMAGINE Photogrammetry toolbox available in ERDAS IMAGINE 2014 (Version 14.00, © 1990-2013 Intergraph Corporation)) was used to process the Pleiades-1 tri-stereo data. As discussed in [13,17,28], in order to obtain the final DEM, point clouds were extracted from each image pair (F-N, F-B, and N-B), and from the merging of the resulting datasets. The application of a semi-global matching (SGM), and the use of the RPC files delivered together with the original scene, allow the orientation and georectification for the imagery. Tie points (TP) for the exterior orientation were automatically identified and improved by adding 12 ground control points (GCP). These points represent 6 sinkholes: for each one, the point of acquisition along the border and the point of maximum depth were identified in each Pleiades image. These points were added to the TPs to input the real field measures of the depth to the dataset. An overall root mean square error of less than 0.1 pixels shift between images at each TP was obtained.
The resulting merged point cloud was then gridded to a spatial resolution of 5 × 5 m/pixel size. The elevation value for each cell was calculated as the average elevation of all the points within one cell as proposed in [13].

Dolines Identification and Morphometric Measurements
In order to highlight the depressions, and thus facilitate their semi-automatic recognition in the DEM, a morphometric processing was applied. Following the methodology applied for geomorphological mapping in [29], a DEM-based morphometric map was produced ( Figure 4). The map represents the combination of profile curvatures, with classes from negative to positive values, and the slope values. The resulting values vary from negatives (concave pixels) to positives (convex pixels).
The contours of the dolines were extracted automatically using GIS and checked with the photointerpretation of the Pleiades orthoimages. In general, the border of most collapse dolines is rather abrupt, with sinkholes opening in a surrounding flat terrain (e.g., ancient lava flow). The border thus corresponds to an abrupt change in angle, from near horizontal to steep or near vertical, and can easily be traced in the DEMs. In more gentle terrains, the knickpoints (changes from convex to concave profile) were extracted from the DEM, and elevation profiles were traced across the dolines. The position of the change of the convex curvature was compared with the existing x,y position on the shape of the same feature. Moreover, contours were extracted from the DEM and compared with the morphologies of the depressions in the orthoimage ( Figure 5).

Figure 4.
The morphometric map shows the overlay of the sinkholes as mapped in the field on the depressions extracted from the combination of the profile curvature and slope maps. The green color in the legend represents pixels with a negative curvature value ranging from low to high value of slope (from −5 to −1); red pixels have a positive value of curvature value, ranging from low to high value of slope (from 1 to 5). The zero value, in bright yellow, represents flat areas. Black bold numbers identify the codes given to name the mapped dolines.
The contours of the dolines were extracted automatically using GIS and checked with the photointerpretation of the Pleiades orthoimages. In general, the border of most collapse dolines is rather abrupt, with sinkholes opening in a surrounding flat terrain (e.g., ancient lava flow). The border thus corresponds to an abrupt change in angle, from near horizontal to steep or near vertical, and can easily be traced in the DEMs. In more gentle terrains, the knickpoints (changes from convex to concave profile) were extracted from the DEM, and elevation profiles were traced across the dolines. The position of the change of the convex curvature was compared with the existing x,y position on the shape of the same feature. Moreover, contours were extracted from the DEM and compared with the morphologies of the depressions in the orthoimage ( Figure 5). . The morphometric map shows the overlay of the sinkholes as mapped in the field on the depressions extracted from the combination of the profile curvature and slope maps. The green color in the legend represents pixels with a negative curvature value ranging from low to high value of slope (from −5 to −1); red pixels have a positive value of curvature value, ranging from low to high value of slope (from 1 to 5). The zero value, in bright yellow, represents flat areas. Black bold numbers identify the codes given to name the mapped dolines.

DEM Validation
In order to check the accuracy of the resulting DEM, and use it to extend the morphometric data collected in the field for a limited number of sinkholes with our elevation model, two separated methodologies were applied to extract the horizontal and vertical accuracies: (a) a qualitative evaluation of the horizontal accuracy was performed using the Pleiades orthoimages acquired in 2019, where the shapes of the sinkholes were overlaid ( Figure 5); (b) the depth of the dolines was measured in the field using a Lietz Sokkisha (now Sokkia, Tokyo) rangefinder with 50 cm base and an operating range of 10-250 m and an accuracy of around 5%. The data collected in the field were compared with the same parameters extracted from the DEM using statistical tests. The position of the dolines was acquired using a handheld rugged GPS (Trimble Juno equipped with the extension ESRI's ArcPad), allowing us to identify the previously numbered dolines (from Landsat and ASTER images), and add the missing ones (not visible in the low-resolution satellite images). A total of 89 dolines were field checked.

DEM Validation
In order to check the accuracy of the resulting DEM, and use it to extend the morphometric data collected in the field for a limited number of sinkholes with our elevation model, two separated methodologies were applied to extract the horizontal and vertical accuracies: (a) a qualitative evaluation of the horizontal accuracy was performed using the Pleiades orthoimages acquired in 2019, where the shapes of the sinkholes were overlaid ( Figure 5); (b) the depth of the dolines was measured in the field using a Lietz Sokkisha (Now Sokkia, Tokyo) rangefinder with 50 cm base and an operating range of 10-250 m and an accuracy of around 5%. The data collected in the field were compared with the same parameters extracted from the DEM using statistical tests. The position of the dolines was acquired using a handheld rugged GPS (Trimble Juno equipped with the extension ESRI's ArcPad), allowing us to identify the previously numbered dolines (from Landsat and ASTER images), and add the missing ones (not visible in the low-resolution satellite images). A total of 89 dolines were field checked. The depth validation was performed making use of a non-parametric signed-rank Wilcoxon test for non-normal distributed samples using R software [30]. This test compared the field-based and the DEM-derived depth measurements to assess whether their distributions are significantly different. Additionally, Pearson's correlation test and linear regression models were calculated to compare the two datasets. The correlation coefficient is a statistical measure of the strength of the relationship between the relative movements of two variables. The correlation coefficient ranges between −1.0 (negative correlation) and +1.0 (positive correlation), with absolute values of 1 representing perfect correlations. These statistical tests were performed on the original field-based dataset (n = 89) and on a subset of the population after stripping of the outliers identified through R software (n = 84).
Extracted depressions range in depth from a minimum of 1 to a maximum of 62 m. Doline perimeters vary between 49 and 868 m, whereas their plan surfaces range between 989 and 61,677 m 2 .
To evaluate the accuracy of the extracted depth values, statistical tests were performed as described in Section 2.5. The results (p-values) of the tests are reported in Table 4. The signed-rank Wilcoxon tests indicate no significant variations between the two paired distributions (p-value > 0.05) both for the entire population and the subset without outliers. The differences between the depth measured in the field and those extracted from the DEM (∆D) for each doline are shown in the boxplots of Figure 7, classified by type. The summary of the basic statistics calculated for the datasets is reported in Table 5. The difference ∆D represents a measure of the vertical accuracy of the DEM extracted from Pleiades tri-stereo images. The boxplots in Figure 8 testify a good accuracy of the population without outliers (mean ∆D value of 2 ± 2 m). The dataset with outliers shows a mean ∆D value of 3 ± 4 m. The main differences represented by the 5 dolines identified as outliers range from 9 to 28 m, whereas the maximum error in the dataset without outliers is 7 m. As shown in the right boxplots of Figure 8, the main errors (both mean and median >2 m) are referred to funnel-shaped and bowl-shaped dolines.
Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 18 The boxplots in Figure 8 testify a good accuracy of the population without outliers (mean ΔD value of 2 ± 2 m). The dataset with outliers shows a mean ΔD value of 3 ± 4 m. The main differences represented by the 5 dolines identified as outliers range from 9 to 28 m, whereas the maximum error in the dataset without outliers is 7 m. As shown in the right boxplots of Figure 8, the main errors (both mean and median >2 m) are referred to funnel-shaped and bowl-shaped dolines.  Table 4. The resulting p-value (<<0.05) and correlation coefficients (>0.9) indicate a good positive correlation between the values derived from field and DEM measurements, with an excellent correlation coefficient of 0.97 for the dataset without outliers.
Furthermore, the linear regression best-fit models calculated for the entire population ( Figure 9A) and the subset without outliers ( Figure 9B Table 4. The resulting p-value (<<0.05) and correlation coefficients (>0.9) indicate a good positive correlation between the values derived from field and DEM measurements, with an excellent correlation coefficient of 0.97 for the dataset without outliers.
Furthermore, the linear regression best-fit models calculated for the entire population ( Figure 9A) and the subset without outliers ( Figure 9B) show high R 2 values of 0.83 and 0.95, respectively. The slope of the trend lines (a value of 1 would mean a perfect match between field and DEM derived measurements) is 0.91 confirming a good accuracy of our DEM. The dolines with the highest errors are those presenting dense vegetation cover or human-built structures that obscure the view of the real bottom surface (Figure 10).

Conclusions
The use of high-resolution satellite imagery and digital elevation models (DEM) for the extraction of morphometric parameters of different positive and negative landforms is extremely useful especially over large and remote areas where detailed topographic maps are not available. The Pleiades satellite constellation delivers images with a ground resolution of 0.5 m, and their use in a GIS environment allows to automatically extract very precise planimetric measurements. We used Pleiades images to extract morphomet-

Conclusions
The use of high-resolution satellite imagery and digital elevation models (DEM) for the extraction of morphometric parameters of different positive and negative landforms is extremely useful especially over large and remote areas where detailed topographic maps are not available. The Pleiades satellite constellation delivers images with a ground resolution of 0.5 m, and their use in a GIS environment allows to automatically extract very precise planimetric measurements. We used Pleiades images to extract morphomet-

Conclusions
The use of high-resolution satellite imagery and digital elevation models (DEM) for the extraction of morphometric parameters of different positive and negative landforms is extremely useful especially over large and remote areas where detailed topographic maps are not available. The Pleiades satellite constellation delivers images with a ground resolution of 0.5 m, and their use in a GIS environment allows to automatically extract very precise planimetric measurements. We used Pleiades images to extract morphometric measurements of the doline fields opening in the basaltic plateau of Azrou, in the Middle Atlas (Morocco). To verify the accuracy of the vertical measurements calculated from the Pleiades-derived DEM, we carried out a ground-truthing on a subset of 89 different dolines (out of a total of 357 identified using low-resolution satellite data including LANDSAT 7 ETM+, ASTER, and ESA Sentinel 2). We found an overall good correlation (p-value < 0.05 and correlation coefficient >0.9) between the DEM-derived depths and those measured with telemeter in the field. The largest differences were obtained on deep funnel-shaped dolines and on shallow bowl-shaped dolines, where disturbance is created by tall trees, human-built structures (dry walls), boulders, or a combination of these features.
Our study shows that very-high-resolution satellite images such as Pleiades (0.5 m ground resolution), are extremely useful in obtaining precise planimetric measurements on large landforms such as dolines. DEMs extracted from Pleiades stereo images deliver good results on the vertical scale if doline floors are not masked by vegetation or disturbed by large debris (boulders or dry walls), with mean depth errors of 2 ± 2 m.
These results will allow for extending the information extracted in the field to the dolines present in inaccessible areas of the plateau, allowing to obtain more information to complete their morphological profiling and propose a model on their unknown origin.
The results of our work demonstrate that tri-stereo satellite images are suitable for the morphometric assessment of shadow-prone negative and positive landforms, and might also be used for the analysis of skylights and collapses on lava flows on the Moon and on Mars, where vegetation or human-built structures are completely lacking, as preparatory surveys for future robotic and human explorations of extraterrestrial caves [31][32][33].