High-resolution monitoring of landslides with UAS photogrammetry and digital image correlation

ABSTRACT

Periodically monitoring landslides is a key factor for supporting the realisation of hazard warning systems and risk reduction in the corresponding neighbourhood areas. Although satellite remote sensing solutions can be considered for low spatial resolution monitoring, this approach is still inappropriate for high spatial resolution investigations. Ground-based Radar Interferometry is also a widely used technique that allows for working at a proper spatial resolution, but it can often be an overbudget solution for most applications. Instead, photogrammetric surveys based on Unmanned Aerial System (UAS) imagery appear as a very interesting approach in terms of both spatial resolution and flexibility in temporally repeating the survey. Motivated by this observation, this work investigates the use of multi-temporal UAS surveys for landslide monitoring. To be more precise, Digital Image Correlation (DIC) has been applied to orthomosaics generated from different UAS photogrammetry surveys to compute the area’s deformation map. Compared with a reference GNSS survey, the results obtained using NHAZCA IRIS software and an in-house DIC approach show a deformation estimation accuracy of approximately 0.1 m, a reasonable accuracy for landslides moving at moderate velocity.

KEYWORDS:

• digital image correlation (DIC)
• unmanned aircraft system (UAS)
• GNSS
• deformation map
• monitoring

Introduction

Different techniques and sensors for landslide investigation, monitoring and early warning provide different kinds of information with different reliability (Baroň & Supper, 2013). The reliability of measurements is of paramount interest in any monitoring system (DiMatteo et al., 2017). Effective landslide monitoring is often obtained using RADAR SAR Interferometry (Lazecký et al., 2015; Leva et al., 2003; Luzi et al., 2004; Pieraccini & Miccinesi, 2019; Turrisi, 2017). Its flexibility, together with its long-range, accuracy and wide field of view, posed RADAR technology, both airborne (Notti et al., 2010; Wasowski & Bovenga, 2022) and ground-based (Antonello et al., 2008; Bellotti et al., 2014; Crosta et al., 2017), as a reference in this field. Such a technique makes building an accurate map of deformations representative of a wide area possible. However, this kind of technology, such as the Ground-Based RARAD, is limited by its high costs. Consequently, several alternative methods, often based on integrating different techniques, have been investigated to enable high-resolution landslide monitoring (Mucchi et al., 2018; Strozzi et al., 2010; Xiong et al., 2020). In this context, UAS photogrammetry is an emerging methodology that ensures high spatial resolution and flexible surveying periodicity, which are key factors in effectively monitoring an active landslide, considering its velocity (Angeli et al., 2000; Mucchi et al., 2018). Furthermore, UAS photogrammetry can conveniently map hard-to-reach places and impervious areas (Cefalo et al., 2011; Lindner et al., 2016; Sestras et al., 2021; Turner et al., 2015). Differently from interferometric RADAR techniques that straightforwardly produce a deformation map of the viewed area, in the photogrammetric case, successive 3D reconstructions of the area of interest should be appropriately compared. Once properly georeferenced, the derived 3D topographic models of the area can be used to obtain DEMs (Digital Elevation Model) that can be employed to investigate volume changes and profile variations (Blasone et al., 2014; Milan et al., 2011; Taddia et al., 2019) in the area, through the computation of DEM of Difference (DoD). Instead, a 2D deformation map of the studied area can be obtained by applying Digital Image Correlation (DIC) to data related to multi-temporal acquisitions (McCormick & Lord, 2010; Sutton et al., 2017). In particular, DIC has already been applied to orthomosaics obtained from photogrammetric surveys (Lucieer et al., 2014; Puniach et al., 2021; Shi & Liu, 2015) to obtain deformation maps of landslide areas. Nevertheless, the accuracy assessment of the obtained deformation results has only partially been investigated in such works, and an uncertainty indication on the obtained estimates is not provided.

This work aims to test the photogrammetric survey-DIC approach in a rather complex scenario, where large deformations have also been caused by human operator interventions, and develop a reasonable reliability/uncertainty index of the obtained deformation results. To this aim, first, orthomosaics of the case study area have been produced by well-established photogrammetric processing procedures of multi-temporal UAS imagery. Then, deformation estimates obtained through commercial software, namely NHAZCA IRIS and an in-house DIC approach (DIC-Flo) developed by the University of Florence, are compared. The proposed in-house solution represents a relatively standard realisation of multi-resolution image phase correlation with subpixel deformation estimation; however, it also provides estimates of the uncertainties of the determined deformations, which can also be exploited to improve the overall results in the multi-resolution case.

Materials and methods

This section describes the study area, the photogrammetric surveys, and the implemented DIC methods.

Study area

The investigated area is a portion of the well-studied Ca’Lita landslide (Cervi et al., 2012). Several studies have been conducted, and different monitoring systems have been installed during the last 20 years (Mulas et al., 2020a). Previous studies on the Ca’Lita landslide show that its velocity range varies from 10 m to 10 cm per month; from moderate to slow, according to Hungr (Hungr et al., 2014). Hence, the photogrammetric surveys in this research have been tailored to the landslide velocity characteristics. The Ca’Lita landslide (sized 0.9 Km²) is inside the administrative boundary of Baiso, in the Reggio Emilia Province (Italy) (Figure 1). The landslide is in the Secchia River Valley and is located in the North-East mountainside of the northern Apennines (Cervi et al., 2012). Ca’Lita has a longitudinal length of approximately 3 km and a maximum width of about 1.4 km (Mulas et al., 2020a). It is a large (Lindner et al., 2016) compound landslide composed of a roto-translational rockslide in the head zone (Cruden & Varnes, 1993; Hungr et al., 2014), with flysch rock masses and in the downslope, it has a translational earth slide-earthflow acting clayey complexes and debris material from the degradation of flysch rock masses (Borgatti et al., 2006).

Figure 1. The study area (red region) (Google Earth, 2021).

The landslide of Ca’Lita has a long succession of reactivations that have threatened the safety of infrastructures and people in the crowning and valley areas. In April 2004, an initial paroxysmal reactivation caused significant retrogression and advancement of the landslide (Mulas et al., 2020b).

Photogrammetric surveys

In the study area, characterised by a diffuse, relatively short vegetation coverage, eight targets have been permanently installed on the ground to guarantee a long-lasting reference for tracking landslide movement over the years. Figure 2 shows the locations of the permanent targets. The main landslide movement direction, i.e. the main flow direction, in the considered area, is approximately from North-West to South-East. The landslide velocity is usually faster on its central part and vice versa on its sides. Consequently, the targets have been distributed approximately along a line orthogonal to the main flow direction to check the velocity variations along its width properly.

Figure 2. Locations of the permanent targets.

Two different photogrammetric surveys were performed at two different epochs, using two UASs and cameras (camera characteristics are reported in Table 1):

• October 2020, DJI Phantom Pro, using the built-in camera

• April 2021, DJI Matrice 300 RTK, mounting DJI P1 camera.

Table 1. Camera characteristics.

Characteristic	DJI Phantom 4	Zenmuse P1
Sensor	CMOS 1”	Full Frame
Pixel	20 MP	45 MP
Lens	FOV 84.00° 8.80 mm (equivalent to 24 mm)	FOV 63.50° 35 mm
Shutter speed	1/2000 sec	1/2000 sec
ISO range	100–3200 (automatic) 100–12800 (manual)	100–25600
Range aperture	f/2.8 - f/11	f/2.8 – f/16

More than 20 ground targets were temporarily installed in the considered area during both campaigns: 14 Ground Control Points (GCPs) and 7 Check Points (CkPs) in 2020 (see Figure 3(a)), 12 GCPs and 11 CkPs in 2021 (see Figure 3(b)).

Figure 3. Orthophotos of the case study area in (a) 2020, (b) 2021. Coordinates are reported in a local reference system.

GCPs and CkPs were quite homogeneously distributed over the case study area, even if not all the landslide locations were accessible. The reference coordinates of the targets were determined with Network Real-Time Kinematic (NRTK) GNSS measurements, receiving real-time corrections from the HxGN SmartNet (HxGN, 2022) network. In both cases, the imagery has been processed using Agisoft Metashape, with all the settings at a level of accuracy “High”, performing camera self-calibration (obtained reprojection error: 0.51 pix and 0.30 pix, respectively). In order to limit the self-calibration error, photogrammetric acquisitions were performed with both nadiral and oblique camera views, ensuring orthogonal roll angles (acquisitions over a double grid) and a reasonably high redundant network (most of the tie points visible by at least nine images). Table 2 summarises the characteristics of the photogrammetric surveys.

Table 2. Summary of the characteristics of the photogrammetric surveys.

Characteristic	October 2020	April 2021
UAS	DJI Phantom 4 Pro	DJI Matrice 300 RTK
Flying altitude (above ground)	40.2 m	64.3 m
Number of images	1413	2951
GSD	1.1 cm	0.8 cm
Forward overlap	80%	80%
Side overlap	60%	60%
Reprojection error	0.51 pix	0.30 pix

Orthomosaics and DEMs have been produced by means of Agisoft Metashape for both the surveys and at two different spatial resolutions (5 cm and 1 cm, named, respectively, low- and high-resolution hereafter).

To easily integrate the performed survey with other technical maps, the output has been computed according to the European Terrestrial Reference Frame ETRF2000 (ETRS89).

Deformation map by means of digital image correlation

DIC is an image processing technique to quantify pixel displacements between two digital images. The rationale is computing displacements by maximising a properly defined functional, usually related to the correlation between two windows in the considered images. When dealing with georeferenced deformation map generation, displacements should be calculated by comparing accurately georeferenced (co-registered) images collected at different epochs (Pan et al., 2008).

In this work, the performance of IRIS software, commercialised by NHAZCA S.r.l., is compared with that of an in-house DIC solution (DIC-Flo) implemented at the University of Florence. In both cases, the inputs of the DIC analysis are two orthomosaics (at the same spatial resolution) associated with the October 2020 and the April 2021 photogrammetric surveys. Figure 4 shows the two surveys’ corresponding portions of the orthomosaics (spatial resolution = 5 cm). A visual comparison between (a) and (b) allows to find both some common graphical elements and some different features, e.g. the human interventions visible on the top-right of the image in Figure 4(b).

Figure 4. Corresponding portions of the orthomosaics used for the DIC analysis: (a) October 2020, (b) April 2021.

In-house DIC implementation (DIC-Flo)

The proposed in-house solution (DIC-Flo) is a quite standard realisation of DIC, implemented in the frequency domain, i.e. image phase correlation. It supports multi-resolution investigations and subpixel deformation estimation. In addition, it also provides estimates of the estimation uncertainties. This subsection first summarises the main characteristics of DIC and image phase correlation and then provides the mathematical details of the more specific aspects of the proposed approach. Readers familiar with DIC and phase correlation may consider skipping the following description until equation (5)(5) $\begin{array}{c}\widehat{\mathrm{\Delta }x}=i,\mathrm{i}\mathrm{f}i\le 2^{k-1},\\ \widehat{\mathrm{\Delta }x}=i-2^k,\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{a}\mathrm{d}\end{array}$(5) .

Let I₁ and I₂ be co-registered (accurately georeferenced) orthomosaics produced, at the same spatial resolution, from imagery collected in two different epochs. Let W₁ be a relatively small portion of I₁ and W₂ be a moving window of the same size as W₁ in I₂. Then, the rationale of DIC is that the movement of the central pixel in W₁ can be reasonably estimated by searching in I_2, which is the position of W₂ that maximises the zero-normalised cross-correlation between W₁ and W₂:

(1) $r_{\mathrm{\Delta }x,\mathrm{\Delta }y}=\frac{{\sum }_{\left(x,y\right)\in W_1}[I_1\left(x,y\right)-{\bar{I}}_1\left][I_2\left(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y\right)-{\bar{I}}_2\right]}{{\sigma }_{I_1}{\sigma }_{I_2}}$(1)

Where ${\sigma }_{I_1}$, ${\sigma }_{I_2}$, and ${\bar{I}}_1$, ${\bar{I}}_2$ are standard deviations and averages of I₁ and I₂, respectively. It is worth noticing that different choices for the functional in (1) can also be considered, but an exhaustive examination of the possible alternatives is out of the scope of this work. The reader is referred to (Hoyt et al., 2006); Hsieh et al. (2008) for more detailed descriptions and comparisons of the possible alternatives.

The use of (1) for deformation mapping should be as follows:

(2) $\left(\widehat{\mathrm{\Delta }x},\widehat{\mathrm{\Delta }y}\right)=arg\underset{\mathrm{\Delta }{x}^{\prime },\mathrm{\Delta }{y}^{\prime }}{max}{r}_{\mathrm{\Delta }{x}^{\prime },\mathrm{\Delta }{y}^{\prime }}$(2)

where$\left(\widehat{\mathrm{\Delta }x},\widehat{\mathrm{\Delta }y}\right)$ is the estimated movement of the central pixel in W₁. Let (x₀,y₀) be the coordinates of such pixel in I₁; the estimated coordinates of its corresponding point in I₂ are (x₀+$\widehat{\mathrm{\Delta }x}$, y₀+$\widehat{\mathrm{\Delta }y}$). (2) should be iteratively used for (x₀,y₀) all over the image in order to compute the entire deformation map. Transforming the deformation map in pixels to a metric scale can be done by straightforward multiplication by the pixel size.

The result of such a procedure depends on the size of the considered windows, where a larger size typically corresponds to a higher stability of the results, despite the detection of some small deformations may be partially lost, while a smaller window size could induce lower result stability.

Given the above observation, a computationally efficient implementation of the procedure is typically needed to enable the use of reasonably large window sizes. To this aim, it is worth noticing that the cross-correlation operation between can be easily expressed by means of convolution, and hence, thanks to the convolution theorem (Oppenheim & Schafer, 1975), it can be equivalently computed in the frequency domain, i.e. through the Fourier transform. When dealing with signals with discrete domains, as in this case, the above considerations should be slightly modified by substituting the cross-correlation with its circular version. From a practical point of view, the implementation in the frequency domain is usually much faster, thanks to the use of the Fast Fourier Transform (FFT) (Oppenheim & Schafer, 1975).

The approach in the frequency domain, applying an element-wise magnitude normalisation to the cross-power spectrum coefficients, leads to the phase correlation approach formulation: let F₁ and F₂ be the results of the Fourier transform applied to the two corresponding windows in the spatial domain W₁ and W₂, where the window size is set to 2^k, being k an integer number, and k ≥ 1, in order to use the FFT efficiently. Then, the element (i,j) of the normalised cross-power spectrum $R^{\prime }$ between F₁ and F₂ is computed as

(3) $R{ }_{ij}^{\prime }=\frac{F_{1,ij}{F}_{2,ij}^{\ast }}{\left|F_{1,ij}{F}_{2,ij}^{\ast }\right|}$(3)

where $F_{m,ij}$ is the (i,j) element in F_m, for m={1, 2}, and $\ast$ stands for the complex conjugate operator. Similarly, the element (i,j) of the cross-power spectrum $R$ between F₁ and F₂ can be computed as

(4) $R_{ij}=F_{1,ij}{F}_{2,ij}^{\ast }$(4)

The cross-correlation r can be computed by applying the inverse Fourier transform to $R$, and, as in (2), the estimated movement $\left(\widehat{\mathrm{\Delta }x},\widehat{\mathrm{\Delta }y}\right)$ for the central pixel in W₁ is set according to the index of the maximum value in $R$. In a similar manner, the movement can be assessed by considering the maximum point of $r^{\prime }$,which is the inverse Fourier transform of $R{ }^{\mathrm{\prime }}$, as well. The phase correlation approach, i.e. using $R{ }^{\mathrm{\prime }}$ and $r^{\mathrm{\prime }}$, is usually computationally more efficient than spatial-domain cross-correlation computations thanks to its frequency domain implementation. Furthermore, it is usually less sensitive to certain factors, such as noise and occlusions. Consequently, the movement is assumed to be determined hereafter using the phase correlation method. It is worth to notice that several alternatives to phase correlation have also been investigated in the literature, however a detailed analysis of such methods is out of the scope of this paper. The reader is referred, for instance, to Bickel et al. (2018) and Hsieh et al. (2008) and the references therein for such analysis.

Let the element index start from 0, and let (i,j) be the position of the maximum value in $r^{\prime }$, then the detected movement is:

(5) $\begin{array}{c}\widehat{\mathrm{\Delta }x}=i,\mathrm{i}\mathrm{f}i\le 2^{k-1},\\ \widehat{\mathrm{\Delta }x}=i-2^k,\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{a}\mathrm{d}\end{array}$(5)

and

(6) $\begin{array}{c}\widehat{\mathrm{\Delta }y}=j,\mathrm{i}\mathrm{f}j\le 2^{k-1},\\ \widehat{\mathrm{\Delta }y}=j-2^k,\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{a}\mathrm{d}.\end{array}$(6)

Similarly to the previously presented case, the computation of $\left(\widehat{\mathrm{\Delta }x},\widehat{\mathrm{\Delta }y}\right)$ shall be iteratively repeated all over the I₁ image to compute the entire deformation map.

From equations (5)(5) $\begin{array}{c}\widehat{\mathrm{\Delta }x}=i,\mathrm{i}\mathrm{f}i\le 2^{k-1},\\ \widehat{\mathrm{\Delta }x}=i-2^k,\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{a}\mathrm{d}\end{array}$(5) and (6(6) $\begin{array}{c}\widehat{\mathrm{\Delta }y}=j,\mathrm{i}\mathrm{f}j\le 2^{k-1},\\ \widehat{\mathrm{\Delta }y}=j-2^k,\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{a}\mathrm{d}.\end{array}$(6) ), it is quite clear that the maximum detectable movement is limited by (half of) the window size. Furthermore, to make the computation quite reliable, choosing a window size larger than the (expected) maximum deformation value is suggested.

When dealing with high spatial resolution orthomosaics and significant landslide movements between two considered epochs, the above consideration may take to using large window sizes, which leads to a high computational burden. A multi-resolution approach has been developed to reduce the computational burden of the above conditions. The proposed method is composed of the following two steps:

• First, a low-resolution analysis (orthomosaic at 5 cm spatial resolution) is executed, determining initial estimates of the movements.

• Then, a high-resolution analysis (orthomosaic at 1 cm spatial resolution) is executed to determine only small variations (≤15 cm) on the previously estimated movements. In our current implementation, the threshold on the maximum detectable variations can be set as a function of the uncertainty of the estimate at the previous step, hence spatially adapting its value depending on the reliability of the already available information.

The steps above are summarised in Figure 5.

Figure 5. DIC-Flo two-step workflow.

Multi-resolution analysis can easily be generalised to a larger number of steps, if needed. Since the landslide can deform at meter level between the two epochs, the window size should be set to an order of magnitude larger than the expected deformation at the first step. Differently, in the second step, the window at the second epoch is already centred at the previously estimated shift; hence, the window size should consider only a possible relatively small variation in the estimated value due to the use of higher spatial resolution data.

In both steps, subpixel maximum estimation is obtained by locally interpolating the values of $r^{\prime }$ around its maximum point. Local interpolation is done by means of a quadratic polynomial (second-order degree polynomial in two variables). Then, the subpixel shift (≤1 in absolute value on both axes) from the previously found maximum point is computed by determining the position of the maximum value of the interpolated quadratic polynomial. Local interpolation in our current implementation is done on a 3 × 3 pixel neighbourhood to reduce the computational burden.

In order to determine the uncertainty on the computed shift, the value of the full-width half maximum (FWHM) of $r^{\prime }$ close to its maximum is computed, along both the axes, obtaining the values FWHM_x and FWHM_y. Then, the uncertainties are set as follows: σ_x = FWHM_x/2, and σ_y = FWHM_y/2.

Such uncertainty values only consider the very local shape of $r^{\prime }$ close to its maximum value. Nevertheless, several local maximums may be detectable in $r^{\prime }$: noise, repetitive patterns (e.g. grass) and quite significant variations in the considered area can also lead to detecting wrong maximum locations. Generally, when the two highest peaks (separated, not in adjacent pixels) have quite similar values, identifying the correct maximum is not always reliable. Let p₁ and p₂ be the values of the two highest peaks, with p₁ ≥ p₂; then the detected shift is set to unreliable if the p₂/p₁ ratio is smaller than the threshold μ. In our current implementation μ = 0.9. Furthermore, a relatively large, detected shift concerning the window size typically indicates an unreliable estimate. Any detected shift larger than ¼ of the window size in our implementation is marked as unreliable.

The approach presented above clearly just allows to assess the 2D landslide movements. Nevertheless, the acquired imagery can also be used to produce DEMs. If DEMs are also available, the Up coordinate variation can be simply obtained by subtracting the DEM height (on the second epoch) evaluated on the new point positions from the initial point heights (determined from the first DEM). Then, the error on the determined height variation is given by the combination of the errors on the DEMs and on the propagation of the error on the localisation of the point on the second epoch orthomosaic to its height. For simplicity, hereafter, the first factor (the error on the computed DEMs) is neglected, even if in certain cases it may lead to a significant contribution, whereas let us just focus on the second term, in particular on the propagated error on the height of the point at the second epoch. Error propagation is implemented by taking into consideration of the 2D point position uncertainties (i.e. σ_x, σ_y and an estimate of the correlation σ_xy), which are assumed to have already been computed, as previously explained. Then, the derivatives of the terrain height on the considered point are numerically computed (e.g. from terrain height variations), and the uncertainty σ_z is obtained by propagating σ_x, σ_y to the Up direction according to the terrain height derivatives.

NHAZCA IRIS software

Pre-processing operations are usually performed by IRIS to standardise the images (named Master and Slave) and to facilitate movement detection, e.g. mean normalisation and smoothing were used to minimise the effects of different lighting and vegetation that affected the Master and Slave orthomosaics. After the pre-processing step, displacement computation with the “phase correlation” algorithm was performed. A median filter (Mikolajczak & Peksinski, 2016) is applied by IRIS at the end of the process to decrease the noise and improve visualisation. A multi-resolution approach has been considered in this case as well.

Results

The validation results of the photogrammetric reconstruction are reported in Table 3 for the two considered surveys (October 2020, April 2021) of the Ca’Lita landslide, where the Root Mean Square Error (RMSE) on control and check points is shown. Error is shown in projected coordinates (E = Easting, N = Northing, U = Up).

Table 3. Photogrammetric reconstruction: RMSE on GCPs and CkPs.

	E [cm]	N [cm]	U [cm]	Total [cm]
	October 2020
GCPs	1.1	0.7	0.9	1.6
CkPs	1.6	1.7	1.4	2.7
	April 2021
GCPs	1.2	0.6	1.9	2.3
CkPs	1.1	1.3	1.6	2.4

Then, Table 4 shows the reference displacements (measured with GNSS) for the set of eight points, shown in Figure 2, used to validate the DIC-based deformation estimation.

Table 4. Reference point displacements (GNSS-based).

Point	E [m]	N [m]	U [m]	2D [m]
01	0.018	0.008	−0.085	0.020
02	0.015	0.012	−0.069	0.019
03	1.774	−0.789	−0.231	1.942
04	1.586	−0.873	−0.055	1.810
05	1.716	−0.795	−0.227	1.891
06	0.244	0.085	−0.007	0.258
07	0.158	0.020	0.030	0.159
08	0.006	0.011	−0.063	0.013

Figure 6 shows the deformation map of all considered areas computed with the NHAZCA IRIS software, using multi-resolution analysis, applied to the low-resolution orthomosaic. Similar results have also been obtained by means of the DIC-Flo approach, despite the results on some areas have been marked as unreliable. In the latter case, DIC-Flo has been applied at first as a single resolution analysis.

Figure 6. Deformation map from DIC, computed with NHAZCA IRIS software, and position of the reference points in Table 4.

A comparison between the numerical results obtained by using NHAZCA IRIS and DIC-Flo is reported in Table 5.

Table 5. Error on the estimated displacements: comparison between estimation errors of NHAZCA IRIS and DIC-Flo (orthomosaic resolution = 5 cm).

	IRIS			DIC-Flo
Point	eE [cm]	eN [cm]	e2D [cm]	eE | σE [cm]	eN | σN [cm]	eU | σU [cm]	e2D [cm]
01	0.3	4.2	4.2	6.2 | 5.4	4.6 | 4.1	1.8 | 1.2	7.7
02	−31.5	−11.1	33.4	5.4 | 3.6	4.1 | 4.0	−2.9 | 0.1	6.8
03	0.2	8.9	8.9	−7.9 | 9.7	6.2 | 3.7	−7.6 | 36.9	10.1
04	−4.8	7.3	8.7	−10.9 | 15.2	6.9 | 5.4	−11.4 | 52.5	12.9
05	7.2	−0.6	7.2	5.5 | 5.1	12.7 | 5.1	0.6 | 5.1	13.8
06	−4.4	−3.4	5.6	-	-	-	-
07	12.7	7.3	14.6	-	-	-	-
08	0.5	8.9	8.9	6.3 | 5.9	6.7 | 5.1	5.3 | 1.8	9.3

To be more precise, Table 5 shows the best results obtained by applying the two software to the photogrammetric products at 5 cm spatial resolution: best results have been obtained by using a 64 pixel × 64 pixel window for NHAZCA IRIS and a 256 pixel × 256 pixel window for DIC-Flo. Errors (obtained by subtracting the reference values, GNSS-based, to the estimated ones) are reported along the three axes (projected coordinates), namely Easting, Northing, Up.

With a slight abuse of notation, the uncertainties previously indicated with σ_x, σ_y, σ_z are now named σ_E, σ_N, σ_U.

DIC-Flo marked the results on two of the considered points as unreliable, hence they have not been reported in the corresponding cells in Table 5. On the remaining points, DIC-Flo obtained an average 2D error of 10.1 cm, whereas IRIS average error on the same points was 11.9 cm. Maximum 2D error is 33.4 cm for IRIS and 13.8 for DIC-Flo. IRIS also provided quite reasonable estimates for the shifts of the two points marked by DIC-Flo as unreliable.

The second step of DIC-Flo has also been executed, using a 256 pixel × 256 pixel window, in order to exploit the photogrammetric results at their full spatial resolution: the obtained results are reported in Table 6.

Table 6. Error on the displacements: comparison between estimation errors of NHAZCA IRIS and DIC-Flo (orthomosaic resolution = 1 cm).

	DIC-Flo
Point	eE cm]	eN [cm]	eU [cm]	e2D [cm]
01	−2.9	−6.5	2.0	7.1
02	−0.9	1.1	−3.0	1.4
03	−13.9	−6.9	−52.3	15.6
04	−16.7	3.8	−40.4	17.1
05	0.5	5.7	2.0	5.7
06	-	-	-	-
07	-	-	-	-
08	−1.0	2.1	6.5	2.4

Similarly, NHAZCA IRIS was applied to the full resolution orthomosaic, however the results obtained in our tests, with error still at decimetre level, were generally worse than those presented in Table 5, hence they are not reported in Table 6. The obtained average 2D error of DIC-Flo is of 8.2 cm for the six points whose estimated displacements were previously marked as reliable. It is worth to notice that the use of a quite small spatial window on a high resolution orthomosaic typically leads to the presence of several local maxima in the cross-correlation, hence the DIC-Flo uncertainty estimates in this case provide too optimistic values and are not reported in the table. On the other hand, using a large spatial window remarkably increases the computational burden, in particular for high-resolution orthomosaics, making the computation slow and quite difficult to be completed if the computer resources are not sufficient. A similar reasoning, along with the significant appearance variations in certain of the considered areas, can also be used to explain the deterioration on the software performance at certain points, in particular for NHAZCA IRIS.

Discussion

This work aimed at investigating the effectiveness of photogrammetry and DIC as tools for monitoring land deformations in complicated, vast and risk-prone areas. The whole Ca’Lita landslide area, of about 0.9 km², is partially covered by vegetation, mostly short grass (70%), some bushed area (2%) from 0.5 to 1.5 meters tall, and very sporadic areas with tall vegetation (<0.3%). In such a context, deploying a reliable and robust monitoring system could be a costly initiative, e.g. the costs related to the maintenance of a complex monitoring system, such as a RADAR SAR, could be too demanding for small administrations such as local municipalities. In this framework, the use of UAS photogrammetric surveying and DIC analysis can represent a less expensive and more flexible solution, alternative to more complex system.

Both the photogrammetric surveys considered in this work have been planned to generate products at approximately 1 cm spatial resolution (GSD = 1.1 cm and 0.8 cm, respectively). Image overlapping was 80% (forward)-60% (side), 14 and 12 control points have been introduced in the photogrammetric workflow for the two reconstructions, and both nadir and oblique views have been used during the imagery acquisition, in order to ensure a more robust camera network. Processing has been performed by means of Agisoft Metashape, with settings at “High” accuracy level. Despite the “Highest” accuracy option is also available, the computation in such case crashed in the used workstation, probably due to insufficient computational/memory resources given the quite high number of images to be processed. The validation of the photogrammetric results, reported in Table 3, shows that the obtained errors are quite reasonable (error on GCPs and CkPs of 1.6 cm, 2.7 cm in the 2020 survey, and 2.3 cm, 2.4 cm in the 2021 one, with error on the Up coordinate around 2 cm in both the cases on CkPs). Despite in this work GCPs have been fully employed in the photogrammetric processing procedure, it is worth to notice that nowadays high-level UASs are often provided with RTK GNSS receivers, reducing the need for GCPs, and hence easing the use of UAS photogrammetry also in areas hard to reach (hence particularly useful for instance in case of natural hazards).

For what concerns the use of DIC for deformation map generation, the comparison between NHAZCA IRIS and the in-house DIC-Flo (Table 5) showed that in both the cases movements can be tracked with a 2D accuracy around 10 cm (using orthomosaics at 5 cm spatial resolution), with slightly better numerical results for DIC-Flo on 6 of the considered points (average and maximum 2D errors are 11.9 cm and 33.4 cm for IRIS, 10.1 cm and 13.8 for DIC-Flo). Nevertheless, NHAZCA IRIS allowed to obtain reasonable shift measurements also on points with estimates marked as unreliable from DIC-Flo. Overall, NHAZCA IRIS allowed to track deformations on a larger number of locations, probably also thanks to the post-processing operations applied by IRIS to regularize the obtained results.

Running the software on the orthomosaics at their full resolution allowed to slightly improve the 2D deformation estimates in the DIC-Flo case, but the improvement in this case is limited by the local changes in the appearance of the different areas, due both to deformations and to the vegetation growth (visible on certain reference point neighbourhoods, as shown in Figure 7), and by the reconstruction error.

Figure 7. A and b respectively represent point 03 and point 04 surrounding areas in April 2020; a’and b’ respectively represent point 03 and point 04 in October 2021.

The deformations visible in certain areas can cause a degradation of the performance of the software, in particular when applied to high spatial resolution orthomosaic, i.e. when the window size can hardly be significantly increased.

Interestingly, a comparison of the results reported on the “e_E | σ_E”|“e_E | σ_E” and “e_N | σ_N”|“e_N | σ_N” columns of Table 5 shows that the uncertainties σ_E and σ_N are usually quite realistic given the values of the errors (e_E and e_N) with respect to the reference GNSS measurements.

Despite having gaps, such as in the DIC-Flo case, in the produced deformation map is clearly something that should be avoided, if possible, this straightforwardly reflects the low reliability of the estimates provided in such locations. This is, for example, the case of point 07 in Table 5, which is actually located very close to an area involved in human interventions in the second epoch (see Figure 8, top-right of Figures 5(b) and 6).

Figure 8. Comparison of the neighbourhood area of point 07 in 2020 (a) and 2021 (b).

Finally, for what concerns the assessment of the movement on the Up direction, this is in general affected by both the error on the photogrammetric product generation and on the propagation of the DIC deformation assessment error. Consequently, the error with respect to the reference measurements is usually larger in this case with respect to the 2D one. Nevertheless, it is worth to notice that the larger error and uncertainty on the Up direction on certain points (e.g. 03 and 04) is mostly due to a quite fast terrain height variation in their neighbourhoods, e.g. due to landslide movement and terrain moved by bulldozers for the rearrangement of hydraulic drainage works (Figure 8).

The investigation conducted in this work showed a very good potential in determining landslide variations up to decimetre level, as shown by the results obtained by both the considered software. It is worth to notice that the use of very high-resolution orthomosaics in ideal working conditions may allow to obtain even a higher accuracy of the estimated deformation map. However, the presence of remarkable changes in the considered area may lead to unreliable estimates, in particular when dealing with very high-resolution investigations, when the spatial window size is often limited by its associated computational requirements. Despite the NHAZCA IRIS software provided very reasonable results in the low-resolution case, being able to compute deformations even on quite critical areas, the results obtained in the high-resolution case did not improve those at low resolution. In accordance with the obtained results, providing a reliability estimate of the obtained results, as in the DIC-Flo case, showed to be useful in particular when the estimate uncertainty is large: on the one hand, this allows to discarding unreliable estimates and, on the other hand, this information can be used in order to determine to size of the neighborhood to which the searching area can be confined for the higher resolution steps.

Conclusions

The presented paper investigated the performance of the integration between UAS photogrammetry and DIC in deformation monitoring on a landslide body in a real case study. The used UAS was the user grade DJI Phantom 4 Pro and DJI Matrice 300 RTK. Ca’Lita landslide is an active landslide, with a displacement rate up to 50 meters per year. Installing ordinary landslide monitoring systems such as extensimeters end inclinometers could be time demanding, and, for a moderate velocity landslide, as this one, the life of an inclinometer can be very short. Other more performing instruments such as interferometric RADAR could be expensive and need to be permanently installed and maintained during all the monitoring period.

This study shows that the integration of UAS photogrammetry and DIC can be an effective deformation monitoring solution, ensuring an average accuracy on 2D movement assessment around 10 cm in our tests, and allowing the generation of a detailed deformation map, which is a proper tool for interpreting landslide dynamics.

Investigating the use of GNSS-controlled triangulation, implemented using RTK GNSS measurements from the receiver onboard of the UAS, without the need (or with a reduced number) of GCPs will be considered in our future work: this kind of solution can be of particular interest in case of natural hazards and in all the cases where the access to the area to be monitored is difficult or risky.

The considered commercial software, NHAZCA IRIS, provided very good performance in the considered low-resolution case. Nevertheless, the results obtained in our tests in the high resolution one, generally worse than in the low-resolution case, prove the usefulness of providing an index of reliability of the obtained results, which could be used, as by DIC-Flo, also to confine the optimization searching area, hence reducing the chances of wrong maximum selection. The proposed DIC-Flo approach provided quite interesting results, in particular for what concerns estimating the uncertainty on the provided deformation measurements; however, differently from the NHAZCA IRIS case, gaps are present in its generated deformation map (when the estimates are considered unreliable). Its performance shall be tested and optimized in a much wider number of scenarios in our future investigations.

Acknowledgments

The authors would like to thank Dr. Giovanni Bertolini, Regione Emilia-Romagna - Agenzia per la Sicurezza Territoriale e Protezione Civile Servizio Sicurezza Territoriale e Protezione Civile di Reggio Emilia, Via Emilia Santo Stefano 25, 42100 Reggio Emilia; Prof. Alessandro Corsini, Dr. Marco Mulas, Dr. Giuseppe Ciccarese; Department of Chemical and Geological Sciences, Via Università 4, 41121 Modena; Prof. Grazia Tucci, Department of Civil and Environmental Engineering, University of Florence, via di S. Marta, 3 - 50139 Firenze, Italy.

The authors want to thank NHAZCA in a special way, for making the IRIS software available to perform the presented tests and for giving the necessary support (https://www.nhazca.it/).

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

Data are available on request from the authors. The data that support the findings of this study are available from the corresponding author [Mugnai F.] upon reasonable request (https://drive.google.com/drive/folders/1EfSI1VtoWHugYD5XWCK8MvTtbkMZzvJ2?usp=sharing).

Additional information

Funding

This research received no external funding