Applications of Image-Based Computer Vision for Remote Surveillance of Slope Instability

Landslides and slope failures represent critical hazards for both the safety of local communities and the potential damage to economically relevant infrastructure such as roads, hydroelectric plants, pipelines, etc. Numerous surveillance methods, including ground-based radar, InSAR, Lidar, seismometers, and more recently computer vision, are available to monitor landslides and slope instability. However, the high cost, complexity, and intrinsic technical limitations of these methods frequently require the design of alternative and complementary techniques. Here, we provide an improved methodology for the application of image-based computer vision in landslide and rockfall monitoring. The newly developed open access Python-based software, Akh-Defo, uses optical flow velocity, image differencing and similarity index map techniques to calculate land deformation including landslides and rockfall. Akh-Defo is applied to two different datasets, notably ground- and satellite-based optical imagery for the Plinth Peak slope in British Columbia, Canada, and satellite optical imagery for the Mud Creek landslide in California, USA. Ground-based optical images were processed to evaluate the capability of Akh-Defo to identify rockfalls and measure land displacement in steep-slope terrains to complement LOS limitations of radar satellite images. Similarly, satellite optical images were processed to evaluate the capability of Akh-Defo to identify ground displacement in active landslide regions a few weeks to months prior to initiation of landslides. The Akh-Defo results were validated from two independent datasets including radar-imagery, processed using state of the art SqueeSAR algorithm for the Plinth Peak case study and very high-resolution temporal Lidar and photogrammetry digital surface elevation datasets for the Mud Creek case study. Our study shows that the Akh-Defo software complements InSAR by mitigating LOS limitations via processing ground-based optical imagery. Additionally, if applied to satellite optical imagery, it can be used as a first stage preliminary warning system (particularly when run on the cloud allowing near real-time processing) prior to processing more expensive but more accurate InSAR products such as SqueeSAR.


INTRODUCTION
Numerous monitoring approaches, including ground-based radar (e.g., Rosenblad et al., 2015), InSAR (e.g., Ferretti et al., 2007;Ferretti, 2014), Lidar (e.g., Lin et al., 2013;Tortini et al., 2015;Kromer et al., 2017;Williams et al., 2018;Williams et al., 2019;Holst et al., 2021), seismometers (e.g., Suriñach et al., 2005), and more recently computer vision (O'Donovan, 2005;Zach et al., 2007;Wedel et al., 2009;Javier et al., 2013;Sánchez Pérez et al., 2013;Kim et al., 2020;Hermle et al., 2022), have been used to study landslides and slope instability. However, the elevated cost, complexity, and intrinsic technical limitations of these methods frequently necessitates the development of alternative and complementary techniques. At the regional scale, due to its ability to measure submillimetre deformation over large areas, InSAR has become one of the most powerful remote sensing monitoring techniques. However, a technical limitation is the restricted detection of measurement points over steep slopes as InSAR measures deformation along the line of sight (LOS) relative to the satellite orbit (i.e., ascending or descending) (Ferretti, 2014). In the case of more detailed and real-time monitoring of targeted mountain slopes, ground-based Radar, and Lidar (Kromer et al., 2017) can provide high resolution slope deformation monitoring (e.g., Tarchi, 2003;Vallee, 2019), but the relatively high operational costs of such techniques often limit their implementation, particularly in isolated areas with rugged topography. Hence, introducing a lower cost technique comparable to ground-based Lidar and Radar is of great interest.
Geophysical tools such as seismometers and infrasound sensors can also record signals generated by rockfalls and landslides. However, the applicability of such techniques strongly depends on the quality of the recording network, magnitude of the events, and the presence of technical expertise to properly interpretate the signals. An important limitation of these sensors is their inability to forecast landslides and rockfalls before their occurrence (e.g., Zimmer and Sitar, 2015;Ulivieri et al., 2020). Recently, a number of studies have been undertaken using computer vision codes and surveillance video cameras to forecast landslides (Amaki et al., 2019) and monitor and assess rockfalls (Kim et al., 2020); they have mainly focused on processing video frames with very short time lapse sequences. In remote and rugged mountainous areas, the challenges of providing both sufficient power and an effective means of transmitting continuous video, limits the use of realtime video datasets. These power and telemetry challenges can be addressed by collecting and processing daily static images with longer time-intervals (from hours to daily) while still enabling continuous monitoring.
In this paper, we provide an improved methodology for the application of computer vision in landslide and rockfall monitoring. The newly developed open access Python-based software, Akh-Defo (available via GitHub repository), uses optical flow velocity, image differencing and similarity index map techniques (Kim et al., 2020) to calculate land deformation (land displacement and rock fall) from a triplet of images (ground-based and satellite). Stacking of the processed triple image sets allows us to calculate the average displacement velocity over an extended period of time. Here, Akh-Defo is applied to two different datasets in British Columbia (BC), Canada and in California, United States of America (USA). The datasets include: 1) Ground-based camera images, with hourly to daily acquisitions, collected during the summer 2021 for the Plinth Peak slope of the Mount Meager Volcanic Complex (MMVC) in BC ( Figure 1); and 2) Daily orthorectified satellite images, provided by Planet Labs, for the Plinth Peak slope (July 28-13 August 2021) and the 20 May 2017, Mud Creek Landslide in California (April 5-31 May 2017) (Figure 1).

Research Objectives
This work has two main objectives: 1) Test whether the Akh-Defo software can complement InSAR LOS limitations via processing of ground-based optical imagery to identify rockfalls and measure land displacement in steep-slope terrain. To achieve this, we chose to apply Akh-Defo on the Plinth Peak slope of MMVC due to the presence of steep slopes and the rugged alpine mountain setting. Independent multi-temporal radar imagery processed by TRE-Altamira scientists using SqueeSAR software (Ferretti et al., 2011) provides high precision measurements of ground displacement for validation.
2) Test whether the Akh-Defo software applied to satellite optical images can identify ground displacement in active landslide regions a few weeks to months prior to initiation of landslides. To achieve this, we investigate the 20 May 2017, Mud Creek catastrophic landslide which has been very well studied and includes very high-resolution pre-, syn-, and post-landslide Lidar and photogrammetry digital surface elevation datasets (Warrick et al., 2019 and references therein).

Plinth Peak Slope, Mount Meager
The slope of interest is located on the NNW flank of the MMVC in the Garibaldi Volcanic Belt of British Columbia, Canada, an area characterized by temperate rainforests, alpine glaciers, and recently observed volcanic fumaroles. The topography consists of steep, highly altered, and unstable volcanic rock slopes which are at least partially covered with snow; visibility is often limited due to fog or low-lying cloud ( Figure 1). The MMVC consists of overlapping volcanoes with a protracted history of basaltic to rhyodacitic volcanism generating lava flows and domes, pyroclastic deposits, and associated rock avalanche deposits from >1.9 Ma to present; these sit on Triassic to Tertiary Coastal Plutonic Complex intrusions and metasedimentary units (e.g., Read, 1990;Russell et al., 2021). The MMVC has erupted explosively at least twice in the past 25,000 years with the most recent event in 2360 B.P. (Russell et al., 2021); pyroclastic density currents from these eruptions were sourced at high elevation near the present-day Plinth Peak and since 2016, low level fumarolic degassing has been observed through glaciovolcanic caves on the west flank of Plinth Peak (e.g., Warwick et al., 2022).
The most recent and largest historic landslide in Canada occurred during the summer 2010 from the southern flank of Mount Meager; less than 3 km South of the current study area (Mokievsky-Zubok, 1977;Guthrie et al., 2012;Hetherington, 2014;Roberti et al., 2017;Roberti et al., 2018). Previous remote sensing and geomorphologic studies of the MMVC indicate that it hosts numerous historical and recent landslide deposits and numerous unstable areas have been identified, many of which are believed to be susceptible to future landslide events (e.g., Mokievsky-Zubok, 1977;Friele and Clague, 2004;Guthrie et al., 2012;Roberti et al., 2017;Roberti et al., 2018;Roberti et al., 2021). Hetherington (2014 and references therein) identify unstable areas on the western side of Plinth Peak and note that the eastern flank appears to have a very low factor of safety that can be presumed to be at the point of failure. The failure surface is described as being deep (Hetherington, 2014;Roberti, 2018;Roberti et al., 2021) such that should a landslide occur in this location, a major bulk of Plinth Peak could fail. Consequently, as with the 2010 landslide, it could dam the Lillooet River and cause widespread flooding downstream.

Mud Creek Landslide
The Mud Creek Landslide is located between Monterey and Morro Bay on the Big Sur coast of California, USA ( Figure 1). Geomorphologically, Mud Creek coincides with the rugged Big Sur coast at the western edge of the Santa Lucia Mountains. This area includes numerous peaks higher than 1,000 m elevation, many of which are located less than 10 km from the coast. Geologically, the Mud Creek and surrounding mountain bedrock is composed of Mesozoic granitic and pre-Cretaceous metamorphic rocks, Miocene marine sedimentary rocks, and heterogeneous Mesozoic rocks which form the Franciscan Assemblage mélange overlying unconsolidated colluvial deposits (Warrick et al., 2019). Structurally, these rugged Creek 1 day before the failure derived from a photogrammetry model (after Warrick, et al., 2019). (H) Oblique aerial view of Mud Creek 7 days after the failure derived from a photogrammetry model (after Warrick, et al., 2019 mountains are related to the Big Sur Bend of the San Gregorio-Hosgri fault system, a region of transpression along the region's transform plate boundary (Warrick et al., 2019). This coastal mountain landslide is in an active landslide region, with coastal retreat rates during the 20th century averaging~0.3 m/year (Warrick et al., 2019).  Figure 1B for the location. Note the impact of Line of Sight (LOS) limitation to capture slope deformation for both ascending and descending satellite orbits. The graphs show the temporal deformation rate for 1) Eastwall Job Creek. 2) Affliction Creek 3) Mosaic Creek slope in the ascending orbit, and 4) Plinth Peak slope, and 5) East-wall Affliction Creek in the descending orbit.

Ground-Based Optical Images
At Mount Meager, with support from Weir-Jones Engineering, we installed a Stardot SD500BN camera, solar power, satellite telemetry, and weather station (all donated by Nupoint Systems Inc.), on a ridge approximately 1.5 km from the slope of interest ( Figure 1). The camera, facing the west-side of Plinth Peak, has a maximum resolution of 5 MP, 2,592 × 1,944 at 30 frames per second (fps). This was remotely programmed to take one image per hour during the daytime for a total of 12 images per day and data (both weather data and images) was transmitted via satellite to the Nupoint Systems and SFU data portals.
From November 2020 to September 2021, a total 1,315 images were collected (see online GitHub Repository for all the images). An automatic image-classification system (Bouti et al., 2020) was used to remove images with significant cloud, fog, and shadow (Supplementary Material). As the study area receives >3-4 m of snow during the winter, we were limited to processing only summer images that include the least amount of snow coverage. Therefore, during this study, we processed approximately 3 weeks of images from July 27 to 18 August 2021. For consistency across the dataset, processed images were taken at four different times (11 a.m., 12 p.m., 1 p.m., and 2 p.m.).

Satellite Optical Images
For both the Plinth Peak and Mud Creek study, high resolution (3 m spatial resolution) Planet Labs optical satellite images were acquired through a research and education license (Planet Team, 2017). The level 3B surface reflectance product of PlanetScope was chosen as it comes orthorectified and corrected for geometric, radiometric, and atmospheric noise. A total of nine cloud-free PlanetScope orthorectified images from July 28 to 13 August 2021, were used to calculate land-displacement rates for the Plinth Peak slope. Ten cloud-free scenes of PlanetScope

Satellite Radar Images
As part of ongoing collaborations between researchers at the Centre for Natural Hazards Research at Simon Fraser University (SFU) and TRE-Altamira, regional landslide and slope stability monitoring has been conducted since 2019 in the Garibaldi Volcanic Belt with a particular focus on the MMVC ( Figure 1). TRE-Altamira scientists processed temporal displacement data for the available SAR satellite (ERS, Sentinel-1, and Radarsat 1 and 2) images from 1992 to 2021 using the SqueeSAR algorithm (Ferretti et al., 2011). As the data were collected from different satellite platforms with different spatial and temporal resolution, we only present data with high temporal and spatial resolution from May 2019 to September 2021 ( Figure 2; Supplementary Table S1 in supplementary material). A total of 36 ascending and 40 descending Sentinel-1 images were processed to obtain 1D line-of-sight (LOS) and 2D true vertical and East-West displacements over the MMVC.

METHODOLOGY
The methodology of this paper is twofold ( Figure 3). The first part deals with optical image processing and the second part discusses processing Radar satellite images using TRE-Altamira's SqueeSAR algorithm. The Akh-Defo software developed in this study ( Figure 3) is applied to ground and satellite optical imagery. The Akh-Defo version for ground based optical imagery consists of two Jupyter Notebook files which can run on any Integrated Development and Learning Environment (IDE) with the ability to read Jupyter Notebook files (e.g., Visual Studio Code, Jupyter Lab, Jupyter Notebook). The first Jupyter Notebook file includes image-preprocessing (discussed in detail in the Supplementary Material) such as training the Deep-Learning Convolutional Neural Network (DLCNN) model, application of the trained model to classify and differentiate noisy from noise-free images, sorting of images based on hours of each day, and finally image-alignment analysis for quality assessment before image processing for change detection analysis ( Figure 3). The second Jupyter Notebook file includes a simple graphical user interface which performs the following steps: 1) image enhancement and definition of the area of interest; 2) static change detection; and 3) dynamic change detection (static and dynamic change detection described below). The Akh-Defo version for satellite optical imagery consists of one Jupyter Notebook file and only processes orthorectified imagery ( Figure 3). The software also includes a simple graphical user interface that performs the following steps: 1) image enhancement and definition of the area of interest; 2) static change detection, and 3) dynamic change detection. Unlike the ground-based version, the satellite version enables geocoding and orthorectification of the final products to real-world coordinates for use with any GIS platforms.

Static Change Detection
Static Change detection includes image differencing such as performing subtraction between two images to identify change of the visual scene or a similarity map to differentiate areas with no change from those with change. We call this static change detection because we can only accurately identify the changes if material (objects) is removed or added within the visual scene between two separate images (Figures 4, 5).

Image Differencing
This technique consists of a pixel-to-pixel value subtraction of an image (Im 1 ) at specific time, t, for a given date to the same pixel in an image for a different date, Im 2 , but at the same time (Figures 4A,5A): This process is performed sequentially for the temporal images from Plinth Peak. By default, image pixel subtraction produces a significant amount of noise, and the quality of the results depends on the time interval between subsequent images as the pixel intensity changes between images. In this study, the minimum time interval between images is about 24 h, but we compensate for the relatively long-time interval by comparing images at a fixed time in subsequent days.

Similarity Indices
The Mean Square error (MSE) and the Structural Similarity Index Mean (SSIM) are two Similarity indices that are frequently used to measure the quality of images (e.g., Palubinskas, 2014;Kim et al., 2020). MSE is considered as a measure of signal fidelity which compares two signals (images) and provides a quantitative score that measures the degree of similarity between them. We assume x {x i | i 1, 2, /, N} and y {y i | i 1, 2, /, N} in which N denotes the number Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 909078 of signal samples (i.e., the number of pixels in an image) such that x i and y i are values of pixels at locations x and y. Therefore, MSE can be calculated as follows (Wang and Bovik, 2009, and references therein): The SSIM technique is based on the human visual system which is particularly powerful in extracting structural information within visual scenes. Hence, to measure similarity, the preservation of signal (image) structure is essential, and it helps to measure the structural distortion (here the distortion related to deformation of the slope). Therefore, SSIM helps to distinguish between structural and nonstructural distortions. Nonstructural distortions include those from the ambient natural environment and instrumental conditions during image acquisition such as 1) Change of luminance or brightness, 2) Change of contrast, 3) Gamma distortion, and 4) Spatial shift of pixels (Palubinskas, 2014;Kim et al., 2020). The structural distortions caused by camera lens and The SSIM measures similarities of images based on the combination of three independent elements, i.e., the brightness values, the contrast values, and the structure of the visual image scene. Assuming that x and y are the locations of two areas within two images of the same size, the similarity of brightness value becomes l(x,y), the similarity of contrast is c(x,y), and the similarity of the area structures is s(x,y) (Palubinskas, 2014;Kim et al., 2020): where μ x and μ y are local sample means of x and y areas, respectively. σ x and σ y are local standard deviations of x and y areas; σ xy allocates the sample cross-correlation of x and y after removing the means of x and y areas. To avoid numerical instability in calculating sample means, variances and correlations, small positive constants C 1 , C 2 , and C 3 were introduced. The SSIM map is produced by computing a sliding window moving pixel-by-pixel across the entire image (Palubinskas, 2014;Kim et al., 2020). A single SSIM score for the whole image is then calculated by averaging the SSIM values across the entire image ( Figures 4B,5B).

Dynamic Change Detection
Dynamic Change detection involves calculating optical flow velocity between triplets of images. The temporal period between each image can be anywhere from intervals of 1 day to several days but the images must be taken at the same time of day in order to minimize the effects of ambient change (e.g., shadows). Subsequently, we can calculate temporal and dynamic change detection from a stack of calculated triplets. We call this dynamic change detection because we can identify the change before (pre-cursor of rockfall and landslide), during (e.g., tracing rockfall and landslide flow path) and after (post-failure) material being removed or added within the scene between a number of subsequent images ( Figures 6C-F, 7, 8).

Optical Flow Velocity
Optical flow is represented as the apparent motion of objects and surfaces between two different timeframes caused by relative motion between the observer (camera) and the scene (e.g., objects, surfaces, and edges) (O'Donovan, 2005;Zach et al., 2007;Wedel et al., 2009;Hermle et al., 2022). The most important assumption in optical flow velocity calculations is constancy in image brightness. This assumption requires that for a short time interval, t 1 to t 2 , an object may change position, while the reflectivity and illumination will remain constant.
Here we process images using an interval of 24 h to a few days, however, it is more difficult to fulfill the assumption of brightness constancy without proper image preprocessing and appropriate image selection for processing. Longer time intervals were compensated by processing images at a fixed time as it helps to remove impact of the angular change in Sun illumination (Figures 7, 8). For processing optical satellite images, while repeated image acquisition is not possible at the same exact time for each day, the images acquired are near vertical (i.e., similar Sun illumination) and we used orthorectified images that have been pre-processed and corrected for atmospheric and geometric noise.  2005;Zach et al., 2007;Wedel et al., 2009;Hermle et al., 2022): By solving the above equation, we obtain the optical flow constraint equation: where ∇I is the spatial intensity gradient (I x , I y ), v = (u, v) is the flow vector for Δx and Δy and I t is gradient through time. Additionally, to better compensate for the image time interval, we calculate the optical flow for triplets of images (frames) instead of pairs of images. The sequence of image triplets is chosen in order to provide overlap among subsequent processed triplet datasets (Figure 7). For instance, (image n , image n+1 , and image n , image n+2 ) forms the first triplet set; the second triplet set is equal to (image n+1 , image n+2 , and image n+1 , image n+3 ) and so on for the rest of the image sequences. The calculated optical flow movements are then stacked for the horizontal Δx Δt , and vertical Δy Δt , obtained from each set of processed triplets. The stacking process allows us to calculate more accurate optical flow motion because it uses larger number of images and averages the results of the calculated triplet optical velocities. Additionally, it provides measurements over a longer time interval of months to even multiple years depending on the number of available processed triplets and the time-window of the triplet datasets ( Figures  6C-F, 7, 8).

Synthetic Aperture Radar
Synthetic Aperture Radar (SAR) satellites acquire images of the Earth's surface by emitting electromagnetic waves and analyzing the reflected signal. As SAR satellites are continuously orbiting the globe, a stack of SAR images can be collected for the same area over time to extract information about changes of the Earth's surface. Each SAR acquisition contains two fundamental properties: amplitude and phase. The amplitude is related to the energy of the backscattered signal and is used in Speckle/Pixel Tracking applications and change detection mapping. Phase is related to the sensor-totarget distance, and it is this which is used in interferometric applications.

InSAR
Interferometric Synthetic Aperture Radar (InSAR), also referred to as Interferometric SAR, is the measurement of signal phase change or interference over time. When a point on the ground moves, the distance between the sensor and the ground target changes affecting the phase value recorded by the SAR sensor. Figure 2 shows the relationship between ground movement and the corresponding shift in signal phase between two SAR signals acquired over the same area. This relationship is expressed as: where Δφ is the change in signal phase, λ is the wavelength and ΔR is the displacement. Consequently, any displacement of a radar target along the satellite Line of Sight (LOS) creates a phase shift in the radar signal that can be detected by comparing the phase values of two SAR images acquired at different times. Due to the acquisition geometry, any phase variation is detected along the satellite LOS. In general, due to the known parameters of SAR satellite orbits and acquisition geometry (approximate North-South orbital direction and side and tilted view), InSAR can detect motion components predominantly in E-W and vertical directions. In contrast, the sensitivity to N-S motion is minimal, even in the presence of strong movement. The comparison of the phase of two SAR images can only be used when the normal baseline (i.e., distance between the satellite tracks during the first and second acquisitions) is no more than the so-called "critical baseline", a parameter depending on the SAR sensor in use.

SqueeSAR
SqueeSAR ™ is the proprietary multi-interferogram technique developed by TRE Altamira which provides high precision measurements of ground displacement by processing multitemporal SAR images acquired over the same area.
SqueeSAR ™ typically needs a dataset of at least 15-20 SAR images to be applied, acquired over the same area with the same acquisition mode and geometry. By statistically exploiting the image stack, SqueeSAR ™ singles out measurement points (MP) on the ground that display stable amplitude and coherent phase throughout every image of the dataset (Ferretti et al., 2011). The MPs belong to two different families: 1) Permanent Scatterers (PS) point-wise radar targets characterized by high stable radar signal return (e.g., buildings, rocky outcrops, linear structures, etc.) and 2) Distributed Scatterers (DS) such as areas of ground exhibiting a lower but homogenous radar signal return (e.g., uncultivated land, debris, deserted areas, etc.). The density and distribution of the MP is related to the resolution of the SAR images used and the local surface characteristics and topography. In general, the MP density increases with the satellite resolution and the presence of man-made structures and decreases with the presence of vegetation. The highest density is reached over urban and bare areas, while it is lower over areas heavily vegetated areas, affected by strong reflectivity changes with time ( Figure 2 and Supplementary  The reference point is assumed to be motionless and selected for its radar properties and motion behavior. This corresponds to a radar target with low phase noise in all the scenes of the imagery, not affected by displacement rate variations (non-linear movement or cyclical deformations) in the covered period. As such, the reference point can still be affected by a linear movement related to regional deformation phenomena. Any regional component can be identified via an independent check with other surface monitoring techniques, such as a GPS network. The selection of the reference point is imagery dependent: if the imagery changes (e.g., number of scenes in the stack or time span), the MP selected as the reference point can consequently change.

Plinth Peak, Mount Meager, Canada
Results from the high temporal and spatial resolution InSAR data (summer 2019-2021), processed with SqueeSAR, identified at least five highly unstable regions on the North flank of the Mount Meager massif. These regions were observed with the ascending SAR satellite geometry and show annual average deformation rates of 5.2, 2.2 and 3.8 cm for Mosaic Glacier slope, Affliction Creek, and East-Wall Job Creek, respectively ( Figure 2). The unstable regions observed via descending SAR satellite geometry show annual average deformation rates of 7.8 and 7.3 cm for East-Wall Affliction Creek and Plinth Peak slopes, respectively ( Figure 2). We processed both ground-based optical images and orthorectified satellite images for Plinth Peak within the same time-window. Ground-based imagery processed with the Akh-Defo dynamic change detection technique measured monthly deformation rates on Plinth Peak from 27 July to 18 August 2021, of between 7 and 27 mm (Table 1; Figure 9). Low velocity rockfall (less than 600 mm/day) ( Figure 8) but frequent ( Figure 10) were observed during the four different daily time periods between 27 July and 18 August 2021. The total detected rockfalls at each time window are as follows: 15 events at 11 a.m., 15 events at 12 p.m., seven events at 1 p.m., and 10 events at 2 p.m. (Figure 8). A small landslide was identified on 11 August 2021, with an area of more than 3,500 m 2 (Figure 9). By contrast, PlanetLab orthorectified imagery processed with the Akh-Defo dynamic change detection was unable to quantify the rockfall activity on Plinth Peak. This is particularly related to the geometry of the camera viewing angle.
The land displacement results are validated from the calculated displacement rates from SqueeSAR for the equivalent period (Table 1) and at a location where we have active displacement data available from both techniques. The validation includes normalization of SqueeSAR annual displacement rates to match displacement rates for the same period of Akh-Defo dataset of 22 days (Table 1). At the location shown in Figure 9E, for the 22 days between 27 July and 18 August 2021, Akh-Defo stacked processing obtained a total displacement of 20 mm (equivalent 27 mm per month) (downward relative to the camera view; see Figure 1 for the camera orientation and location). Comparatively, at the same site and time window (Figure 9G), the SqueeSAR deformation results obtained a ground displacement of 16.9 mm (22.5 mm per month), along SAR satellite LOS (Descending Orbit). Additionally, the PlanetLab.
Our ground-based imagery results show that the Plinth Peak slope is subject to active and frequent rockfalls albeit with slow slope displacement and that the rockfalls are mainly concentrated along very-steep channels within the slope (Figure 10 rockfall frequency and Figure 11 for sources of rockfalls). Rockfall trajectories were validated from the animated image time-lapse (see the online GitHub in the gif_dir folder dataset).
The Akh-Defo static change detection (using Image Differencing and Structural Similarity Index Map) was effective in identifying one small landslide (rockslide) on the Plinth Peak slope, approximately 30-m-wide and 120-m-long (at least 3,600 m 2 ) between 10 and 11 August 2021, although with different sensitivities (Figure 12). For the same landslide, the SSIM visually produced clearer results compared with those using Image Differencing (Figure 12). In fact, we also captured a precursor of this event using dynamic change detection 7 days before failure where the total displacement ranged from 1 to 10 cm from 3 to 10 August 2021 ( Figure 12G). Additionally, we reconstructed the rockslide trajectory and velocity of the rockslide at >60 cm/day between 10 and 12 August 2021 ( Figure 12H).
The mechanism of the rockfalls and the rockslide can be further quantified based on the morphology of the Plinth Peak slope. The slope is composed of a large talus cone (Rapp and Fairbridge, 1997) formed due to the accumulation of rock-debris at its base ( Figure 11). The bedrock of the slope is composed of rhyodacite lavas which may have facilitated relatively uniform weathering across the West part of the Plinth Peak resulting in a single talus slope rather than a series of cones. According to Rapp and Fairbridge (1997), rates of creep movement in active talus cones can be~10 cm/year on the upper part of the cone decreasing to zero towards the base. These rates are similar to the SqueeSAR displacement rates for Plinth Peak of 7.3 cm/year ( Figure 2).
The rockfalls detected by Akh-Defo on the Plinth Peak slope can be categorized into two types of talus deposits, rockfall talus and talus creep (Rapp and Fairbridge, 1997). Rockfall talus related to individual block failing, such as the 11 August rockslide, are characterized by shattering, rolling, and bouncing caused by freeze-thaw, and heavy rain ( Figures 10C,D, 11, 12). One week before the 11 August 2021, rockslide, the camera weather station recorded a sharp drop in temperature as well as heavy rainfall ( Figure 10).
Although, at this time there is no published information regarding the structural and kinematic behavior of the Plinth Peak slope, we performed a preliminary assessment in order to map linear features from the Lidar-based hillshade map. The Plinth Peak slope encompasses at least three sets of lineaments, striking NW-SE, ESE-WNW, and NNE-SSW. Our results show that most of the rockfalls and slope deformation occur within and along the trends of mapped lineaments ( Figure 13) with rockfalls particularly concentrated along the NNE-SSW striking lineaments.

Mud Creek Landslide, California
Akh-Defo processing of orthorectified PlanetScope satellite imagery identified displacement between 50 and 300 cm (Figures 6, 7), 41 days prior to the 20 May 2017, Mud Creek landslide. The timestamp-sequential displacement ( Figure 7) and average displacement ( Figure 6) show two main unstable zones. The upper zone (referred to as the depletion zone by Warrick et al., 2019) showed a maximum of 100 cm displacement on 10 April 2017, 40 days before the catastrophic landslide failure. In comparison, the lower zone (referred to as the accumulation zone by Warrick et al., 2019) showed displacement rates of as little as 60-80 cm for the same time period 40 days before the failure ( Figure 6).
This failure event was followed by accelerating displacement in the lower zone from 80 cm to more than 300 cm at 35 and 6 days, respectively, before the catastrophic failure (Figure 7). Four days before the event, both the upper and lower zones became extremely unstable and showed displacement of more than 300 cm. The above pattern of accelerating displacement in the lower zone at days 35 to 6, preceding the upper zone, might have facilitated the occurrence of the catastrophic failure. This would kinematically validate the loss of slope piedmont support for material from the upper slope zone. Warrick et al. (2019) show that although the lateral flank of the failure moved significantly during the failure, it did not fail completely ( Figures 6G,H). We  Table 2) which increased to a few meters between May 1 and 16, 2017, just 4 days before the catastrophic failure (Figure 7). Although there is limited data and information available regarding the cause of the Mud Creek 2017 landslide, according to Handwerger et al. (2019), the region underwent a 5-year drought before 2017 and the failure occurred following a period of persistent rainfall lasting several days. Simple 1D hydrological modelling suggests that significant increases in pore-fluid pressure led to rapid destabilization of the slope resulting in the catastrophic Mud Creek landslide.

DISCUSSION
Aside from data preprocessing, the main task of Akh-Defo is to perform static and dynamic change detection. Dynamic change detection is a powerful method for slope deformation monitoring as it provides pre-failure warning signals-for example, 7 days before the 11 August 2021, rockslide on Plinth Peak, the software measured >2 cm of slope movement at the rockslide source area ( Figures 12G,H); additionally, 41 days before the 20 May 2017, Mud Creek landslide, the Akh-Defo software measured more than 50 cm of slope movement (Figures 6, 7).
Static change detection consists of image difference maps through subtraction of subsequent images and similarity maps through measuring image similarities based on the combination of three independent elements: 1) brightness values, 2) contrast values, and 3) structure of the visual image scene. Unlike similarity maps, image differencing solely relies on pixel intensity change between subsequent images; hence, it is sensitive to light conditions and produces relatively more noise (Figure 14). For instance, changes in solar elevation angle produce different Sun illumination patterns and different shade patterns on ground surfaces at different hours of the day. We minimized light condition variability related to daily change of Sun elevation angle and Sun azimuth at different hours of the day by processing images acquired at the same time of the day ( Figure 14).
The Dynamic Change-Detection approach can detect slow (subcm) and fast (more than a few metres) slope movement. The user needs to apply other complementary techniques to differentiate the source of calculated displacement such as rockfall, rockslide or slow slope deformation. In this study we used a number of methods to distinguish and interpret the source of slow and fast displacement rates including SqueeSAR, Static Change-Detection and visual observation of the time-lapse images. SqueeSAR (InSAR) is a valuable independent method to validate slow slope displacement (Table 1; Figure 9). Static change detection and visual observation of the image time-lapse both are important to quantify fast slope movement such as rockfalls and rockslides (Figures 10, 12).
Through installation of a low-cost fixed camera and telemetry system, we were able to collect more frequent images (daily to hourly) of the Plinth Peak slope compared to exclusively relying on InSAR ( Figures 9G,E). Statistical comparison between SqueeSAR results and Akh-Defo shows less than 20% difference in displacement rate from Akh-Defo relative to InSAR (Table 1). SqueeSAR processes radar images and relies on the change of distance between the SAR satellite and the target on the ground. SqueeSAR can processes data collected at all times and under various weather conditions such as cloud, fog, and rain. However, SqueeSAR cannot calculate sufficient measurement points in very steep slopes or if the target invisible to the LOS.
Alternatively, Akh-Defo processes optical images and relies on the change of the intensity pattern within the optical image but is restricted to only data collected during the daytime and under clear weather conditions. In this study, the Akh-Defo software was applied to the Planet Labs pre-processed orthorectified optical images to measure the Plinth Peak slope as well as the pre-Mud Creek landslide movements. The calculated land movement from orthorectified satellite optical images are more accurate and easier to project to the real-world trajectories relative to the ground-based optical images, mainly due to the geometric characteristics of the satellite images relative to the ground images. For instance, in the case of orthorectified satellite images, the pre-processed optical images have already been corrected (using digital elevation models) for external image distortion caused by the difference in depth of objects within the image scene (i.e., any distortions caused by topography have been corrected). Akh-Defo uses optical flow algorithms to calculate slope movement. The measurement accuracy depends on the datacollection system such as resolution of the raw image, weather conditions during capturing the image, stability of the camera system (in case of ground-based images) and even time of the day the subsequent images were taken. In our methodology, we were careful to select images for processing to obtain acceptable results for validation with other independent systems such as SqueeSAR (Plinth Peak case study) and Lidar DEM change detection (Mud Creek case study). For instance, we only processed images taken during clear weather conditions (no fog, snow, or rain), we performed image alignment to compensate for the stability of the ground camera system (see Supplementary Material for image alignment process) and we only processed images taken at fixed hour of subsequent days. In addition to careful selection of images and image pre-processing steps such as image alignment and image enhancement, we also integrated similarity thresholds (comparable to InSAR coherence thresholds). The advantage of using this parameter is that it enables calculation of displacement velocity only for pixels with known stable texture and light condition (Figures 15B,D).
Current limitations of the Akh-Defo software include the inability to accurately identify the displacement vector field for the calculated displacement from the ground-based images processed for the Plinth Peak slope. This limitation is inherent to all computer vision optical flow codes and its known as the Aperture problem (Xue et al., 2015). The significance of the aperture problem depends on the size of the object of interest; for instance, if the object of interest is larger than what can be seen through the camera scene view (such as the Plinth Peak slope), we can only estimate an approximate vector field for objects within the slope. In other words, while we cannot identify the movement FIGURE 15 | Deformation rates in millimeter per month for Plinth Peak slope. (A) Deformation rate calculated from ground-based optical images from July 27 to 18 August 2021; note, the colour-bar is shown as a logarithmic scale, deformation rates larger than 20 mm represents traces of slow but continuous rockfalls within the processed time window period. (B) Deformation rate calculated from orthorectified optical satellite images from July 28 to 13 August 2021, deformation map processed with similarity threshold equal to 0.75. (C) Deformation rate calculated from radar images for summer 2019, 2020, and 2021 using TRE-ALTAMIRA SqueeSAR Software. SqueeSAR by default calculate annual deformation rate, we rescaled to a monthly deformation rate. (D) Deformation rate calculated from orthorectified optical satellite images from July 28 to 13 August 2021, deformation map processed with similarity threshold equal to 0.2. of the entire slope relative to Mount Meager itself, we can calculate the displacement rate and vector of smaller landslides and rockfalls withing the field of view of our installed camera.
Another limitation is that the optical flow algorithm can only measure horizontal displacement in a 2-dimensional image space such as (left, right, top, and bottom). This becomes more challenging in the case of ground-based optical image processing to translate the calculated displacement magnitude vector into a real-world vector motion. In contrast, for optical satellite images the solution becomes relatively more practical by processing orthorectified optical images ( Figure 3); we can then reproject the calculated velocity magnitude and vector back to the real-world projection using available digital elevation models (Figures 5C,D).
Our results show that ground-based optical images are more suitable to monitor rockfalls compared to orthorectified satellite optical images. This is mainly related to the different angles of view as the orthorectified satellite images are acquired at nearly nadir while ground-based optical images have a nearly horizontal view ( Figure 15).

CONCLUSION
In this study, we developed software to process ground-based optical images and orthorectified Planet Lab satellite optical images. The ground-based optical images were used to test the capability of the developed software to complement InSAR LOS limitations for the steep slopes of Plinth Peak. Our ground-based camera analysis, using stacks of optical images from 27 July to 18 August 2021, gives comparable and complementary results (less than 20% difference; see Table 1) for the same time-period with SqueeSAR. In addition to defining deformation rates, we were able to characterize velocity of rockfall trajectories spatially within the slope itself. Additionally, we used two techniques, Image Differencing and Structural Similarity Index (SSIM), to identify areas of significant slope change such as landslide areas. Results from ground-based camera monitoring detected movements of the slope and defined areas of highest rockfall on the western flank of Plinth Peak. Our results show that ground-based optical images are more suitable to monitor rockfalls compared to orthorectified satellite optical images.
We also tested the Akh-Defo program to calculate the pre-failure movement 41 days before the catastrophic 20 May 2017 Mud Creek Landslide from the orthorectified Planet Lab satellite imagery and validated our results with those from Warrick et al. (2019). We calculated timestamp (Figure 7) and average ( Figures 6C-G) displacement for the period 40 days before the catastrophic landslide occurred.
In this paper, the following conclusions are presented to enhance cost-effective and near real-time early warning landslide and rockfall monitoring: 1) Radar satellite datasets are a highly effective source to identify unstable areas as a baseline knowledge if followed by less-expensive but continuous ground-based and optical satellite monitoring. 2) Computer vision codes such as optical-flow velocity are highly effective for calculating land deformation semi-quantitively if used with a proper methodological workflow which includes choosing the appropriate type of image dataset and reasonable image pre-processing preparation. 3) Structural similarity maps appear to be more effective in identifying rockfall events and areas of change compared to classic image differencing. 4) The pre-failure Mud Creek displacement results calculated by Akh-Defo software indicate that Akh-Defo program (based on optical flow velocity code) can be used as a first stage preliminary warning system (particularly if run on the Cloud) prior to processing more expensive but more accurate InSAR products such as SqueeSAR.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found on the below Online GitHub Repository https://github.com/mahmudsfu/AkhDefo. This repository includes the following materials: A copy of the developed Akh-Defo Python-based software for Ground-based Image processing. A copy of the modified Akh-Defo program to process Planet Labs orthorectified images Ground-based Plinth Peak image dataset evaluated with the Akh-Defo software Mud-Creek Planet Lab orthorectified satellite images evaluated with the modified Akh-Defo software. Plinth Peak Planet Lab orthorectified satellite images evaluated with the modified Akh-Defo software. A step-by-step instruction manual on how to run Akh-Defo software on Jupyter notebook program.

AUTHOR CONTRIBUTIONS
All authors contributed to writing, processing data and reviewing this manuscript. MM developed the python code and methodology supervised by GW-J and DS. GF and RT processed Radar images using SqueeSAR.