Sliding Time Master Digital Image Correlation Analyses of CubeSat Images for landslide Monitoring: The Rattlesnake Hills Landslide (USA)

: Landslide monitoring is a global challenge that can take strong advantage from opportunities offered by Earth Observation (EO). The increasing availability of constellations of small satellites (e.g., CubeSats) is allowing the collection of satellite images at an incredible revisit time (daily) and good spatial resolution. Furthermore, this trend is expected to grow rapidly in the next few years. In order to explore the potential of using a long stack of images for improving the measurement of ground displacement, we developed a new procedure called STMDA (Slide Time Master Digital image correlation Analyses) that we applied to one year long stack of PlanetScope images for back analyzing the displacement pattern of the Rattlesnake Hills landslide occurred between the 2017 and 2018 in the Washington State (USA). Displacement maps and time ‐ series of displacement of different portions of the landslide was derived, measuring velocity up to 0.5 m/week, i.e., very similar to velocities available in literature. Furthermore, STMDA showed also a good potential in denoising the time ‐ series of displacement at the whole scale with respect to the application of standard DIC methods, thus providing displacement precision up to 0.01 pixels.


Introduction
Space-borne Earth Observation (EO) is one of the most valuable solutions for Earth's surface spatial and temporal monitoring, with fundamental socio-economic benefits [1][2][3][4]. Since the early stage, EO has been considered a promising tool for better managing and preventing natural hazards due, especially, to the unique capability of looking simultaneously at large areas [5,6]. However, some evident limitations for monitoring purposes are connected with temporal and spatial resolution of the images [7][8][9].
Several successful applications of landslide mapping after extreme events such as rainfalls or earthquakes exists [10][11][12][13] and, especially in the last two decades, also monitoring small displacements caused by landslides using Synthetic Aperture Radar Interferometry has been possible [14][15][16][17]. In the last years, new image processing tools were developed and examples of landslide displacement monitoring, based on optical satellite images are available [7,12,[18][19][20][21]. However, only few authors attempted to derive time-series of displacement using long stack of satellite optical images, thus demonstrating the potential of such an approach to derive time of series of displacement along the plane normal to the line of sight with a sub-pixel accuracy [6,20,[22][23][24].
The most common optical satellite missions (e.g., QuickBird, SPOT, LANDSAT, Sentinel-2, WorldView, Pléiades, just to mention few of them) are characterized by spatial resolution ranging between 0.3 and 30 m and a revisit time of some days [25][26][27], therefore, they are not fully adequate to detect landslides at a spatial and temporal scale suitable for continuous monitoring. However, as stated in [4] and [28], the fusion of images collected by different optical satellite sensors has certainly increased the opportunities for cloud-free surface investigation, thus allowing also a general improvement of the temporal resolution.
In the last few years, EO has been affected by an astonishing sensor and platform development and improvement, thanks to small satellites. According to [27,29] and [30], among the main features of these satellites, properly named CubeSats, there are: i) the general small size and weight (a singleunit of CubeSat normally measures 10 × 10 × 11 cm and typically weights less than 1.5 kg), ii) high geometric resolution (~3-5 m/pixel) and iii) daily/near-daily revisit time. CubeSats represent a costeffective EO strategy, allowing unique chances for a wide variety of application fields [8,31].
Planet Labs Inc. has recently made available a large constellation (> 130 units) of optical 3-U CubeSats (10 × 10 × 30 cm), also called PlanetScope or, more commonly, "Doves" [27,30,[32][33][34]. With an average orbit height of about 475 km a.s.l., the main CCD array sensor is able to collect images of 7000 × 2000 pixels, with a footprint of approximately 24 × 8 km 2 . Each Doves collects one image/second along his orbit with a slight overlapping between consecutive scenes. In this way, taking advantage of the whole constellation, images at a daily sampling rate are achievable, thus allowing to obtain long stacks of multi-temporal images in a short period [8,34].
At present, only few authors used "Doves" images for the assessment of landslides [35] or earthquakes-induced displacement [8], while, no attempts of deriving time-series of displacements of landslides using stack of Doves images processed by Digital Image Correlation (DIC) are available in literature.
Recently, different authors presented interesting results derived from the application of DIC analysis on a stack of satellite images for landslide displacement monitoring [22,24] . Specifically, they developed a workflow adapted for satellite optical images with the aim of achieving the best signal/noise ratio by taking advantage of the data redundancy. In this paper, we developed and applied a similar processing procedure for the analyses of a long stack of Planet Scope Images aiming at providing time-series of displacement. The results derived from the application of our processing procedure, named STMDA (Sliding Time Master DIC Analysis), on a large landslide occurred in the Rattlesnake Hills (Washington State, USA) ( Figure 1), is herein presented. The Rattlesnake Hills landslide was selected as a nice example of a recent active landslide (i.e., allowing the availability of a good stack of Doves images) characterized by a relatively large size and amount of surface deformation. Furthermore, the availability of an extensive field monitoring network was a key factor for the validation of the results achieved in our study.

The Rattlesnake Hills Landslide
The Rattlesnake Hills landslide is located near Union Gap, about 10 km southern of the city of Yakima (WA, USA). According to [36][37][38] and the Washington State Department of Natural Resources, the landslide was activated in October 2017, when tension cracks were observed along the North side of the AK Anderson quarry managed by Columbia Asphalt and Gravel, Inc. The slide involved up to 4 million m 3 volume, affecting about 8*10 4 m 2 of the hillside, moving South toward the quarry ( Figure 2). The slope instability process has been controlled by sedimentary and tectonic features, which lead to a movement along a weak 10-15 degrees inclined sedimentary interbed within basalt flow sequence [36][37][38]. The Rattlesnake Hills landslide has been classified as a translational slide, significantly disjointed in the upper part as evidenced by huge tension cracks in continuous expansion [36,38]. Based on [37], the recent quarrying activities at the toe of the landslide may have weakened the whole slope system, thus activating a "wedge" of rock that in the first months of 2018 reached a velocity in the order of 0.5 m/week. Following the landslide triggering, different kinds of contact and remote sensing monitoring systems, such as GPS (Global Position System), seismometers, robotic total stations, multi-temporal UAV (Unmanned Aerial Vehicles) photogrammetric and TLS (Terrestrial Laser Scanner) surveys, TInSAR (Terrestrial Interferometric Synthetic Aperture Radar) and 24-hr webcam have been deployed [36][37][38][39]. According to [37] and [38], in situ displacement measurements confirmed that the Rattlesnake Hills landslide has moved toward S-SW direction, along the 10-15° dipping rupture surface discussed above. The extensive monitoring network allowed to measure the strong displacement trend as well as the landslide bounding features, direction and inclinations of the tension cracks. The landslide showed an accelerating trend between October 2017 and January 2018 with the highest displacement velocities showed in the central portion of the landslide mass in the period mid December 2017 to mid-January 2018 (with a peak of 56 cm/week horizontal velocity). The upper part of the landslide showed, in the monitoring period, velocities lower than the central portion (average horizontal velocity of about 12 cm/week) with a higher vertical component that caused large ground cracks, due to the stretching action. The lowest velocity was recorded in valley wedge mass. According to [36] and [38], this strain behavior can be explained by assuming the stabilizing role of the valley portion on the moving volume, and the subsidence trends of upper portion.

Material and methods
Thanks to the Planet Labs' image dissemination policy for research and education purposes, it has been possible to access a whole stack of 3 m/pixel PlanetScope analytic (ortho-rectified) images [34], over the Rattlesnake Hills landslide area. In total, 63 images of the 139 available have been selected for their low cloud coverage (maximum 20%) over a period of 12 months (i.e., from mid-September 2017 till mid-September 2018) ( Figure 3).

Sliding Time Master DIC Analyses (STMDA)
The first analyses of the available PlanetScope images was performed by the processing chain proposed by [18]. Specifically, the first two steps are included in the orthorectification and georeferencing elaboration made by Planet [8], while the third step (i.e., the correlation phase) is performed using an advanced stacking correlation scheme, similar to the one developed by [22,24] and extensively described below. Advanced stacking DIC analysis allowed to achieve, for each pixel, the value of displacement along the N-S and E-W directions and the correlation coefficient (i.e., an estimation of the "cross-correlation quality"), with an increase of signal to noise ratio [6,21,22,24]. Specifically, DIC analysis is here applied to all the possible image pairs of the selected dataset with a progressive correlation scheme showed in Figure 4. The progressive correlation scheme is applied iteratively using all the available images as the master, thus achieving a fully redundant stack of correlation maps characterized by different time-spans. . Progressive correlation scheme used to achieve the fully redundant stack of correlation maps. This process was iterated using all the 63 PlanetScope selected images as the master (images stack and correlation maps drawing modified from [6]).
The DIC analyses have been carried out by the open source plug-in of ENVI ® , named COSI-Corr (Co-registration of Optically Sensed Images and Correlation) that is based on a sub-pixel image correlation algorithm developed by [18] and [40]. COSI-Corr is a DIC software implemented in the IDL language for geoscience application, with a specific focus on geomorphological surface dynamics monitoring [41]. Specifically, all the correlation maps were obtained using the frequency correlation option available in COSI-Corr that consists in an iterative process made of a preliminary estimation of the pixel wise displacement between two patches followed by a final correlation to retrieve the subpixel displacement. Due to the huge number of image pairs, a batch correlation was performed by using a parameters selection scheme representing the best compromise between computational time and accuracy in displacement estimations for the Rattlesnake Hills case study. Correlation maps of about 6 m/pixel resolution resulting from original pixel size of 3 m, was investigated with a sliding window of 128 pixels and a correlation step size of 2 pixels. More than 3900 correlation maps have been achieved by the fully redundant correlation scheme of the 63 selected images. At this regard, it is worth to note that DIC displacement values are computed along the directions corresponding to the normal plane to the line of sight. For that reason, by using orthorectified images, we got displacements only along the horizontal directions (i.e., N-S and E-W directions).
The output of the above correlation process for each pixel and each component here defined as the "redundant correlation matrix" is represented by a number of time-series of displacement equal to the number of the master images (i.e., all the selected images) achieved by an inversion procedure similar to the one presented in [22,24]. The inversion allows to provide all the time-series (derived from the different master images) a "zero-displacement value" corresponding to time "zero" (i.e., the date of the oldest image of the stack).
For the analyses of the "redundant correlation matrix" an ad hoc processing chain was developed in MATLAB v. R2018a. The first step of the STMDA approach consists in the computation of the average displacement values for all the time-series available at each time instant. This procedure was iterated for each pixel, thus producing the denoised time-displacement stack that can be used for the time-series of displacement. Furthermore, the correlation coefficient has been used as a filtering tool by discarding values with Correlation Coefficient (CC) < 0.8. This procedure allows to denoising the displacement results, thus increasing the robustness of the measurement [22,24]. However, the above described procedure cannot solve the errors contribution due to localized changes in radiometric content that can occur for several causes at the scale of a single pixel, especially for spatially heterogeneous processes, like landslides. To this aim, an approach based on the analysis of Region-Of-Interest (ROI) made of adjacent pixels with a similar displacement behavior was developed. By analyzing the displacement values of each pixel at each time interval, from the redundant correlation matrix we derive, for each time instant and each time-series, one displacement value (for each component) that is statistically representative of the selected ROI (see at Figure 5). The measurements dispersion relative to the selected ROI in the time/displacement graph is reported in Figure 6A.
In order to remove the noisy time-series of displacement and find the best fitting between the derived data, we developed an algorithm for identifying the best I, II or III order polynomial functions that maximize the R 2 coefficient and fitting, at the same time, the maximum number of time-series.
Then, the displacement values from the redundant correlation matrix that do not fit the characteristics found after the above processing step are discarded ( Figure 6B). Following the above procedure, it is now possible to derive a new polynomial regression of the averaged time-series obtained from the filtered redundant correlation matrix, thus achieving the best results in computing the surface displacement in time. At this stage, the standard deviation of the instantaneous displacement resulting from the above described filtering procedure, can be assumed as the precision of the STMDA measurements at each time frame.
Residual error estimation of the STMDA procedure is then carried out analyzing ROIs not affected by displacement (i.e. area outside the landslide not affected by ground movement). In this area, we can assume that all the measured displacement values are due to systematic and random error contribution, hence, the average displacement values should be very close to zero and the dispersion of the measurements is representative of the errors.
With the aim of deriving the one-year long Rattlesnake Hills landslide kinematic history (September 2017 to September 2018) we applied the above described STMDA method to the selected stack of PlanetScope images, by identifying some ROIs that could describe the behavior of specific portions of the moving mass [39]. Specifically, three ROIs have been built on the landslide area with a medium size of about 100 pixels (corresponding to about 6400 m 2 ) and four ROIs outside the landslide area (i.e. stable ground) with a medium size of 120 pixels (corresponding to about 7600 m 2 ), for the error estimation. For each ROI on the landslide area the time-series of displacement are derived, while in the ROIs outside the landslide (stable areas), we performed a statistical analysis of displacement data for error estimations.

Results
By applying the STMDA processing method to the available stack of PlanetScope images, the multi-temporal horizontal displacements pattern of the Rattlesnake Hills landslide between September 2017 and September 2018 was achieved (Figure 7). The displacement vector field and the whole landslide boundary are consistent with the in situ observations from the WSDNR and others [36,37,38,39]. Time-series of displacement has been derived on the three ROIs located inside the landslide boundary showed in Figure 8. Specifically, the maximum computed surface deformation during the monitoring period is between 15 and 16 m in the ROI A and B, while is about 9 m in ROI C. In terms of temporal evolution, it is quite evident an acceleration phase in the first 200 days (corresponding to the period Sept. 2017 to March 2018) and, following, a decelerating phase (April to September 2018). These results are in line with the in situ data reported in [36][37][38][39]. The maximum velocity of about 0.5 m/week was reached in the period mid-February to March 2018 on ROI A and B, while a maximum velocity of about 0.2 m/week is reached in ROI C in March.
Looking at the error bars, we can observe two main evidences: 1) in the first 150 days the length of error bars is larger than in the following phase in all the ROIs; 2) there is a difference in the mean length of error bars between different areas (e.g., in ROI_A the error bars are smaller than in ROI_B and ROI_C). The average value of the error bars, even if slightly different for each ROI, is on the order of 0.01 m. It is worth to note that fracturing processes, tilting movements, vegetation changes, and other changes on the ground surface can generate errors on the correlation process. For example, ROI_B is characterized by lower precision as it corresponds to the upper part of the landslide where extensive fractures were observed, as demonstrated also by the decrease of correlation coefficient map showed in Figure 9.   Systematic error caused by the DIC processing tool (i.e., COSI-Corr) was also estimated by the pixelwise subtraction of displacements obtained from homologous image pairs (i.e., 1 vs. 2 and 2 vs. 1, 15 vs. 54 and 54 vs. 15, and so on), thus achieving differences in the order of 0.006 m.

Discussions
With the purpose of investigating the potential of using long stack of PlanetScope images for landslide monitoring purposes, we developed an original processing procedure called STMDA that was tested on the Rattlesnake Hills landslide. Figure 11 shows the quality improvement in the landslide displacement mapping using the STMDA stacking procedure with respect to standard DIC analysis (that compares a single couple of images), without applying any radiometric correction to the original images. This is evident in the northern sector of landslide that was characterized by the occurrence of ground cracks, thus causing a significant reduction of the correlation coefficient ( Figure  9). However, residual noise is present also in images A, B and C, mainly due by ground cracks caused by vertical movements (reported in [36][37][38][39]), that cannot be detected by DIC analyses ( Figure  9).
In Figure 12, three time-series of displacement corresponding to i) a single pixel, ii) average of pixels in ROI B and iii) results on ROI B achieved from the application of the complete STMDA procedure, are compared. The improvement in the displacement measurement performances using procedure ii) and, especially, procedure iii) (i.e., STMDA) is quite evident. The differences are probably due to the radiometric variations induced by the clouds and land cover coverage (especially in the images acquired in the winter season), that in a future perspective, could be partially solved by radiometric corrections applied to the original images. Time-series derived by the STMDA procedure on the three ROIs located in the landslide area ( Figure 10) revealed a whole activation/resting loop of the landslide with an initial acceleration phase (from October 2017 to April 2018) and the consequent deceleration phase (from May to September 2018) (Figure 8), thus showing the overall landslide surface kinematic pattern in the spatial and temporal domains. Specifically, the Rattlesnake Hills landslide showed different displacement magnitude with the lower values measured in the southern part (ROI_C), as confirmed by field measurements [36,38].
The increasing presence of fractures on some portions of the landslide surface represents the most important source of error for the DIC analyses in addition to the presence of snow in the winter season and lighting changes due to different sunlight incidence during the year. However, the STMDA procedure allowed us to mitigate these sources of error, whose residual contribution is described by the error bars showing, as one can expect, higher values during the winter season ( Figure 8).
Precision estimation has been performed by using the mean value of the error bars for each time in the diagrams presented in Figure 8. In such a way, it has been possible to achieve, for the ROI_A (located in correspondence of the GPS point characterized by the higher movement) a value of 0.01 m (i.e., less than 0.01 pixels), that is one order of magnitude lower than the results of similar studies available in the scientific literature [6,9,22].
The accuracy of displacement achieved by the STMDA procedure using the above described stack of PlanetScope images has been estimated by the comparison with traditional ground monitoring data available from literature (mainly achieved by total stations and GNSS) [38]. The comparison was performed on the ROI_A by sub-sampling the STMDA time-series in order to make the data comparable with the monthly displacement velocity data available from literature [36,38]. Results of the comparison (Figure 13) show the good correspondence between the displacement values with a few exceptions in the winter season and in September 2018. The general trend of displacement over time is very similar (1.48 cm/day average difference) while, as expected, the absolute values are under-estimated due to the neglection of the vertical component and the spatial average of the STMDA method. These results shed light to one of the limitations of the STMDA procedure, i.e., the capability of measuring displacement only along the plane normal to the line of sight that makes the procedure more effective for horizontal movements than for vertical movements if satellite images acquired with a Nadir view are used. However, we must point out that the STMDA procedure can be applied also to aerial or ground based images, that, having a different line of sight, could allow to measure also vertical movements. Figure 13. Comparison of the monthly highest velocity of the Rattlesnake Hills landslide measured by the STMDA procedure herein presented and the available data from literature [36][37][38][39].
DIC techniques are more popular for lab application but, in the last years, some authors presented interesting results also on natural environments, such as landslides, using a variety of optical sensors [6,9,12,22]. Staking DIC methods, like the STMDA procedure herein presented, have the potential to furtherly increase the capability of this technique taking advantage of a huge number of data with a high revisit time, in order to reduce the source of noise [6,[20][21][22]24]. The STMDA procedure based on satellite optical images is complementary to other satellite-based displacement monitoring techniques like Advanced Interferometric Synthetic Aperture Radar (A-DInSAR). A-DInSAR is able to measure movement along the Line of Sight (LOS) direction, providing best performance for vertical movements. However, it is worth to note that A-DInSAR monitoring of the Rattlesnake Hills landslide could have been affected by phase ambiguity problems, due to the high displacement rate. At this regard, it is worth mentioning that A-DInSAR and STMDA are complementary also in terms of displacement velocities. As a matter of fact, considering the classification of Cruden and Varnes [42] A-DInSAR is effective only for extremely slow and very slow landslides, while STMDA is effective for very slow to moderate landslides. Nevertheless, in the perspective of a widespread application of both A-DInSAR and STMDA techniques we cannot forget the main limitations related to the land cover, such as vegetation and snow.

Conclusions
This paper demonstrates the potential of the high temporal frequency satellite optical constellations (in this specific case PlanetScope nanosatellite) for landslide monitoring using the herein proposed Sliding Time Master Digital image correlation Analyses (STMDA) processing method. It has been possible to identify and measure the Rattlesnake Hills landslide kinematic behavior between September 2017 and September 2018, achieving a very good fitting with the displacement data collected by local authorities during multi-temporal field surveys. Precision of measurements are up to 1/100 of the pixel size (3 m in this case), that is one order of magnitude higher than those achieved by other authors using ground based, satellite and aerial optical imagery with traditional DIC approaches. The accuracy estimation remarks small differences (in the order of 0.01 m/day) between the velocity achieved using the herein proposed method and the results from in situ measurements, thus opening incredible perspective for systematic satellite-based landslide monitoring at global level. In addition to landslides, STMDA procedure applied to CubeSats images can be useful for other surface displacement processes like glaciers and earthquakes and it is not affected by problems related to high displacement rate as it occurs with A-DInSAR. However, as for any other remote sensing satellite-based displacement monitoring technique, some drawbacks that can affect the STMDA monitoring results must be considered like the vegetation and land cover and the slope inclination. Nevertheless, STMDA has the potential to be applied also to images acquired by others source such as, aerial and ground based sensors, thus allowing to solve most of the above described issues.