Fusion of SAR Interferometry and Polarimetry Methods for Landslide Reactivation Study, the Bureya River (Russia) Event Case Study

In this paper, we demonstrate the estimation capabilities of landslide reactivation based on various SAR (Synthetic Aperture Radar) methods: Cloude-Pottier decomposition of Sentinel-1 dual polarimetry data, MT-InSAR (Multi-temporal Interferometric Synthetic Aperture Radar) techniques, and cloud computing of backscattering time series. The object of the study is the landslide in the east of Russia that took place on 11 December 2018 on the Bureya River. H-α-A polarimetric decomposition of C-band radar images not detected significant transformations of scattering mechanisms for the surface of the rupture, whereas L-band radar data show changes in scattering mechanisms before and after the main landslide. The assessment of ground displacements along the surface of the rupture in the 2019–2021 snowless periods was carried out using MT-InSAR methods. These displacements were 40 mm/year along the line of sight. The SBAS-InSAR results have allowed us to reveal displacements of great area in 2020 and 2021 snowless periods that were 30–40 mm/year along the line-of-sight. In general, the results obtained by MT-InSAR methods showed, on the one hand, the continuation of displacements along the surface of the rupture and on the other hand, some stabilization of the rate of landslide processes.


Introduction
Earth surface remote sensing methods and technologies, as well as ways to access these data, are being intensively developed now [1,2]. Study and monitoring of catastrophic natural processes are among the most relevant applications of these methods. These methods are successfully applied, for example, for the monitoring of seismically hazardous areas to identify earthquake precursors by analyzing anomalous variations in ionospheric parameters [3], recorded, including by signals from satellite navigation systems [4,5], by analyzing changes in lineament systems revealed during satellite images processing [6] and variations in thermal fields [7], etc. The efficiency of dangerous phenomena forecasting increases greatly when integrating satellite data analysis results with the approaches based, for example, on the application of geomechanical models [8,9] or seismic entropy [10] methods, etc. Remote sensing data are also applied to monitor wildfires [11,12] and their consequences [13], typhoons [14], and other natural disasters.
To register anomalous natural processes, it is promising to use radar remote sensing methods. The main advantage of such methods is their all-day, all-weather capabilities of imaging. Radar interferometry and polarimetry methods using synthetic aperture radar (SAR) data are now important tools for displacement estimation and change detection of the earth's surface structure due to various natural reasons. Publicly available data of

Study Area
The object of the study is a large landslide that occurred on 11 December 2018 in the Far East of Russia (Khabarovsk Krai) on the left bank of the Bureya River from the slope opposite the mouth of its right tributary, the Sredny Sandar River. The location of the Bureya landslide (50 • 34 N, 131 • 29 E) is given in Figure 1. The peculiarities of the Bureya landslide are its occurrence in winter at a temperature of about -30 • C and the formation of a water-ice tsunami. The wave of that tsunami spread along the waterways of the Bureya Remote Sens. 2021, 13, 5136 3 of 20 and Sredny Sandar rivers for several kilometers, destroying the coastal forest on its way. Most of the sliding masses filled the Bureya River. reya and Sredny Sandar rivers for several kilometers, destroying the coastal forest o way. Most of the sliding masses filled the Bureya River.
Although the area where the Bureya landslide occurred is hard-to-reach and u habited (the distance to the nearest settlement is about 70 km), the event caused a mas public response. The landslide massif completely blocked the Bureya River, which in has been the Bureyskoe Reservoir since 2009. The river is the main source for energy duction by large Bureyskaya Hydropower Station, located 80 km downstream. The l slide also created a risk of flooding upstream settlements. Thus, blasting operation water channels were carried out, in January-February 2019. Then the channel wa panded naturally by spring floods and the gradual erosion by the river current. The geological and geographical analysis of the landslide, assessment of its cha teristics, possible mechanisms of its formation and slope failure were carried out in 58] based on in-situ studies. The total volume of displaced rocks at all formation st was estimated at 25 million m 3 . The Bureya landslide created a 50-70 m high and u 550 m wide dam on the river. According to the classification [59], it can be referred the type "fall".
Landslides quite often take place where they have already occurred, as a resu reactivation [60]. In June 2019, during an expedition to the Bureya landslide, it was covered that the main scarp and the landslide body underwent a transformation du the activity of various exogenous processes [56]. Activation of those processes is prob due to the thawing of permafrost. During the expedition in July 2020, organized by Russian Geographical Society, a geological and geomorphological survey of the land area was carried out [61]. The survey showed that various exogenous processes car on within the landslide circus and its vicinities. During the period of fieldwork in western section of the landslide scarp, repeated collapses of large blocks of rock in form of rockfalls were noted, especially after intense rainfalls. The intensive developm of erosion processes within the landslide circus with the formation of ravines up to meters deep was also revealed. Although the area where the Bureya landslide occurred is hard-to-reach and uninhabited (the distance to the nearest settlement is about 70 km), the event caused a massive public response. The landslide massif completely blocked the Bureya River, which in fact has been the Bureyskoe Reservoir since 2009. The river is the main source for energy production by large Bureyskaya Hydropower Station, located 80 km downstream. The landslide also created a risk of flooding upstream settlements. Thus, blasting operations for water channels were carried out, in January-February 2019. Then the channel was expanded naturally by spring floods and the gradual erosion by the river current.
The geological and geographical analysis of the landslide, assessment of its characteristics, possible mechanisms of its formation and slope failure were carried out in [56][57][58] based on in-situ studies. The total volume of displaced rocks at all formation stages was estimated at 25 million m 3 . The Bureya landslide created a 50-70 m high and up to 550 m wide dam on the river. According to the classification [59], it can be referred to as the type "fall".
Landslides quite often take place where they have already occurred, as a result of reactivation [60]. In June 2019, during an expedition to the Bureya landslide, it was discovered that the main scarp and the landslide body underwent a transformation due to the activity of various exogenous processes [56]. Activation of those processes is probably due to the thawing of permafrost. During the expedition in July 2020, organized by the Russian Geographical Society, a geological and geomorphological survey of the landslide area was carried out [61]. The survey showed that various exogenous processes carried on within the landslide circus and its vicinities. During the period of fieldwork in the western section of the landslide scarp, repeated collapses of large blocks of rock in the form of rockfalls were noted, especially after intense rainfalls. The intensive development of erosion processes within the landslide circus with the formation of ravines up to five meters deep was also revealed.  The Bureya River crosses the relatively low mountains of the Bureya area [57]. The absolute elevations of the adjoining dome-shaped watershed heights reach 600-640 m. The cross-section of the Bureya River valley, which had, before the creation of the reservoir, an elevation (at the site of the landslide formation) of 187 m above sea level, within the territory of interest is characterized by a trapezoidal shape with a steep left southern slope and a gentler right northern side. The full Bureyskoe Reservoir level is 256 m above sea level. The reservoir level on 11 December 2018 was 253 m. Because the landslide mass fell into the reservoir, this led to up to a 60 m high overwash, which completely destroyed the taiga on the opposite bank of the river valley in a strip of up to 550 m. The overwash extended to 3.6 km up the Sredny Sandar River valley, located opposite the slope failure zone, destroying the taiga in a strip up to 300 m wide. Upstream of the Bureya River, the overwash is traced at a distance of 7 km, and downstream-at a distance of 4 km, which is consistent with the north-north-east direction of the displacement of sliding masses, oriented more upstream [57]. Figure 3 presents 1 m-resolution panchromatic images (0.62-0.79 µ m spectral band) taken by Geoton-L1 optical sensor aboard the Russian Resurs-P No. 2 satellite on 2 May 2019 and 25 September 2019. These images demonstrate the erosion of the channel created through blasting during the time between these two images, as well as the clearing of the ice cover and the erosion of the channel in the Sredny Sandar River. The difference in width in the narrow part of the cleared channel of the Bureya River exceeds 100 m.
As we can see from Figure 3, the size of a hollow after the landslide is in the plane of the horizon along the sliding line 600-700 m (length), and 400-500 m across.
Hence, it seems relevant to study the zone of the Bureya landslide from the point of view of its possible relapse. The Bureya River crosses the relatively low mountains of the Bureya area [57]. The absolute elevations of the adjoining dome-shaped watershed heights reach 600-640 m. The cross-section of the Bureya River valley, which had, before the creation of the reservoir, an elevation (at the site of the landslide formation) of 187 m above sea level, within the territory of interest is characterized by a trapezoidal shape with a steep left southern slope and a gentler right northern side. The full Bureyskoe Reservoir level is 256 m above sea level. The reservoir level on 11 December 2018 was 253 m. Because the landslide mass fell into the reservoir, this led to up to a 60 m high overwash, which completely destroyed the taiga on the opposite bank of the river valley in a strip of up to 550 m. The overwash extended to 3.6 km up the Sredny Sandar River valley, located opposite the slope failure zone, destroying the taiga in a strip up to 300 m wide. Upstream of the Bureya River, the overwash is traced at a distance of 7 km, and downstream-at a distance of 4 km, which is consistent with the north-north-east direction of the displacement of sliding masses, oriented more upstream [57]. Figure 3 presents 1 m-resolution panchromatic images (0.62-0.79 µm spectral band) taken by Geoton-L1 optical sensor aboard the Russian Resurs-P No. 2 satellite on 2 May 2019 and 25 September 2019. These images demonstrate the erosion of the channel created through blasting during the time between these two images, as well as the clearing of the ice cover and the erosion of the channel in the Sredny Sandar River. The difference in width in the narrow part of the cleared channel of the Bureya River exceeds 100 m.
As we can see from Figure 3, the size of a hollow after the landslide is in the plane of the horizon along the sliding line 600-700 m (length), and 400-500 m across.
Hence, it seems relevant to study the zone of the Bureya landslide from the point of view of its possible relapse.

Data and Methods
The Google Earth Engine platform [63] was used for the analysis of Sentinel-1B Cband backscattering variations. The radar acquired images at a descending orbit with a heading angle of 192° and an incidence angle of 36.4°.
As mentioned above, PSI and SBAS-InSAR methods, based on time series processing of radar images, are widely used to study landslides using SAR [64][65][66][67][68]. Methods based on the analysis of interferometric coherence [22] and polarimetry [47][48][49] are also used. However, unfortunately, the use of full polarimetric analysis is currently inaccessible due to the lack of data obtained from fully polarimetric radars [22]. Thus, because of the limited data set, an analysis of the consequences of the landslide was undertaken using various methods: • Analyzing the results of cloud computing of backscattering time series allowing for sensing geometry and relief features; • Cloude-Pottier decomposition of Sentinel-1 dual polarimetry data; • PSI and SBAS-InSAR with Sentinel-1 data.
MT-InSAR (PSI and SBAS) methods and Cloude-Pottier polarimetric decomposition are used to analyze Sentinel-1B dual (VV and VH) polarimetry data. Open access images [69] were obtained in the IW interferometric mode, with a swath of 250 km at a spatial resolution of 5 m × 14 m. In order to eliminate the effect of snow cover on interferometric measurements, we used images obtained during a snowless period from May to mid-October. Polarimetric and interferometric data were processed using the SARScape software package.
Let us present a brief overview of polarimetric [47][48][49] and interferometric [35][36][37][38][39][40][41] methods. Polarimetric measurements are available from polarimetric radars which transmit both vertically (V) and horizontally (H) polarized waves and receive co-polarized (VV and HH) and cross-polarized (VH and HV) backscattered fields. The result can be summarized in a matrix: which is referred to as the scattering matrix or Sinclair matrix. The elements of this matrix are the complex scattering amplitudes, which describe the polarization states of a radar

Data and Methods
The Google Earth Engine platform [63] was used for the analysis of Sentinel-1B C-band backscattering variations. The radar acquired images at a descending orbit with a heading angle of 192 • and an incidence angle of 36.4 • .
As mentioned above, PSI and SBAS-InSAR methods, based on time series processing of radar images, are widely used to study landslides using SAR [64][65][66][67][68]. Methods based on the analysis of interferometric coherence [22] and polarimetry [47][48][49] are also used. However, unfortunately, the use of full polarimetric analysis is currently inaccessible due to the lack of data obtained from fully polarimetric radars [22]. Thus, because of the limited data set, an analysis of the consequences of the landslide was undertaken using various methods:

•
Analyzing the results of cloud computing of backscattering time series allowing for sensing geometry and relief features; • Cloude-Pottier decomposition of Sentinel-1 dual polarimetry data; • PSI and SBAS-InSAR with Sentinel-1 data.
MT-InSAR (PSI and SBAS) methods and Cloude-Pottier polarimetric decomposition are used to analyze Sentinel-1B dual (VV and VH) polarimetry data. Open access images [69] were obtained in the IW interferometric mode, with a swath of 250 km at a spatial resolution of 5 m × 14 m. In order to eliminate the effect of snow cover on interferometric measurements, we used images obtained during a snowless period from May to mid-October. Polarimetric and interferometric data were processed using the SARScape software package.
Let us present a brief overview of polarimetric [47][48][49] and interferometric [35][36][37][38][39][40][41] methods. Polarimetric measurements are available from polarimetric radars which transmit both vertically (V) and horizontally (H) polarized waves and receive co-polarized (VV and HH) and cross-polarized (VH and HV) backscattered fields. The result can be summarized in a matrix: which is referred to as the scattering matrix or Sinclair matrix. The elements of this matrix are the complex scattering amplitudes, which describe the polarization states of a radar backscattering. Note, that the first subscript of the elements refers to the polarization of the scattered wave while the second subscript refers to the polarization of the incident wave. The basis of polarimetry is the dependence of the polarization state of a radar signal on the physical mechanism of backscattering. Various scattering matrix transformations allow us to reveal these mechanisms and, thereby, to identify different earth's cover types that appear in radar image data. SAR can also operate in dual polarization (dual-pol) mode, for example, Sentinel-1A/B. It employs a single polarization transmission (V) and dual (V and H) coherent reception. The dual-pol mode provides two complex scattering amplitudes, S VV and S VH : The Freeman-Durden and Cloude-Pottier decompositions are most widely used for polarimetric analysis. The polarimetric decomposition methods divide the backscattered field into a sum of several basic components with different scattering mechanisms. In the Freeman-Durden decomposition [47], a covariance matrix is represented as a contribution of three scattering mechanisms, i.e., surface (single bounce) scattering, double bounce dihedral scattering, and volume scattering. This decomposition refers to incoherent decomposition, which uses coherency or covariance matrices. The Cloude-Pottier (entropy/alpha/anisotropy) decomposition [48] is an eigenvalue/eigenvector decomposition, where entropy (H), anisotropy (A), and alpha angle (α) parameters are introduced. Entropy is a measure of the randomness of scattering. Anisotropy can be defined as a normalized difference between the appearance probabilities of scattering components. The parameter α is an indicator of a type of scattering. Using the Cloude-Pottier decomposition in the H-α plane, one can determine nine scattering mechanisms. A dual-polarized version of the entropy/alpha decomposition method is developed by Cloude in [49].
Satellite radar interferometry (InSAR) is a powerful remote sensing technique able to measure surface deformations and displacements with millimeter accuracy.
Classical differential radar interferometry (DInSAR) is based on measurements of phase differences (interferograms) of signals that are proportional to the difference in the distances from the satellite to the earth's surface from close points of two orbits of successive satellite radar flights [15,[35][36][37]. The main DInSAR limitation is the lack of coherence between images obtained when sensing from quite separated orbits (spatial decorrelation), or in the case of significant changes in the time between surveys (temporal decorrelation).
The development of the DInSAR method is the persistent scatterer method [38] and the small baseline subset method [39], based on the analysis of multi-temporal series of images and combined with their modifications under the general name of Multi-temporal InSAR (MT-InSAR). The Persistent Scatterer Interferometry (PSI) method allows us to calculate the time behavior of point persistent scatterers with a quite intensive and stable reflected signal. A multitude of interferometric phase relationships is calculated for these discrete scatterers with respect to one reference image. These phase relationships allow us to evaluate deformation magnitude and velocity at certain time intervals. The Small BAseline Subset (SBAS) approach involves the use of a series of data with small spatial (less than a threshold value) baselines, at the time of surveys, and short time intervals between surveys. Deformations for distributed objects (set of pixels) can be assessed at relatively small coherence values than for PSI.
An application of radar data time-series to generate multiple interferograms allows us to more accurately (compared to DInSAR) evaluate the deformations with help from discrete persistent scatterers (PS) or distributed scatterers (SBAS). The processing and analysis of dozens of interferograms can significantly reduce the influence of the atmosphere, inaccuracies in the reference digital elevation model (DEM), and orbital errors.

Geometric Parameters
Results of radar sensing depend on survey geometry and object geometry. Radar interferometry allows us to measure only one-dimensional v LOS velocities of an object along the line of sight (LOS), i.e., satellite side-looking direction that is orthogonal to the satellite's flight. In fact, the landslide velocity is a vector and consists of the eastern v E , northern v N , and vertical v Z components. The projections of these vectors in the LOS direction depend on the satellite geometry shown in Figure 4. In this figure θ is the incidence angle (the angle between the vertical direction and the radar LOS);r is the LOS unit vector pointing from the ground to the radar; α is the angle between the satellite heading and north direction; ϕ is the slope angle; γ is the aspect angle;û is the (opposite) steepest slope unit vector, along the downward/upward steepest slope direction, assumed positive upward.
The velocity vector projection onto the direction of ̂ vector to the radar operating in the right side-looking mode is as follows [37]: The argument in brackets ( −   3   2 ) is an azimuth look direction.
This expression is valid for both ascending and descending orbits; however, in the latter case, the more descriptive equation is: The Equations (1) and (2) are reduced to the following: = − sin cos + sin cos + cos (3) Figure 4. Descending orbit acquisition geometry for a point located on a slope (modified from [70]).
From Equations (1)-(3) it follows that in a general case, to determine the total velocity vector, three independent orbits with their own sensors are required. In the case of landslides with a known direction of movement, only one equation is needed to determine the movement velocity [70].
= cos (4) where β is the angle between the steepest slope direction and the LOS. The value of cos β is determined by the product of ̂ and ̂ unit vectors. The velocity vector projection onto the direction ofr vector to the radar operating in the right side-looking mode is as follows [37]: The argument in brackets α − 3π 2 is an azimuth look direction. This expression is valid for both ascending and descending orbits; however, in the latter case, the more descriptive equation is: The Equations (1) and (2) are reduced to the following: From Equations (1)-(3) it follows that in a general case, to determine the total velocity vector, three independent orbits with their own sensors are required. In the case of landslides with a known direction of movement, only one equation is needed to determine the movement velocity [70].
where β is the angle between the steepest slope direction and the LOS. The value of cos β is determined by the product ofr andû unit vectors.
cos β = sin θ cos α sin γ cos ϕ − sin θ sin α cos γ cos ϕ + cos θ sin ϕ= sin θ cos ϕ sin(γ − α) + cos θ sin ϕ, From (5) it follows that if the heading angle is equal to the aspect of the slope or their difference equal to π (the satellite movement and the projection of the landslide movement on the horizontal plane are parallel or antiparallel), the value of cos β is determined only by θ and ϕ angles.

Backscattering
A digital elevation model based on TerraSAR-X/TanDEM-X data obtained after the landslide (28 September 2020), was created for the analysis of possible ground movements in the landslide zone. A 3D visualization of this model (5 m resolution) is given in Figure 5a. For the comparison, Figure 5b presents a 30 m-resolution 3D DEM [71] based on Shuttle Radar Topography Mission data (SRTM, 11-22 February 2000). The red lines show the boundaries of the failed slope. In Figure 5b the river level during the SRTM mission was 60 m lower than the current level that was formed when the Bureyskoye Reservoir was filled in 2009 due to the construction of the Bureyskaya Hydropower Station dam. cos = sin cos sin cos − sin sin cos cos + cos sin = = sin cos sin( − ) + cos sin , From (5) it follows that if the heading angle is equal to the aspect of the slope or their difference equal to π (the satellite movement and the projection of the landslide movement on the horizontal plane are parallel or antiparallel), the value of cos β is determined only by θ and φ angles.

Backscattering
A digital elevation model based on TerraSAR-X/TanDEM-X data obtained after the landslide (28 September 2020), was created for the analysis of possible ground movements in the landslide zone. A 3D visualization of this model (5 m resolution) is given in Figure  5a. For the comparison, Figure 5b Figure 6b shows the elevation profiles from DEM SRTM (red line) and TerraSAR-X/TanDEM-X (black line). Measurements in the graph in Figure 6b are shown from west to east along the yellow line in Figure 6a. Figure 6c presents the Sentinel-1B backscattering coefficient measured along this profile in 2018-2021. Similar results for the red line from Figure 6a are given in Figure 6d,e measured from north to south. The backscattering coefficient data were obtained using Google Earth Engine [63]. The 10 m-resolution S1_GRD collection was used for computations. In each time interval (from May to October of each year), the data were averaged and stored as a single image, which was used to construct the profile. In total, there were at least 25 Sentinel-1 images in each time interval.   Figure 6b shows the elevation profiles from DEM SRTM (red line) and TerraSAR-X/TanDEM-X (black line). Measurements in the graph in Figure 6b are shown from west to east along the yellow line in Figure 6a. Figure 6c presents the Sentinel-1B backscattering coefficient measured along this profile in 2018-2021. Similar results for the red line from Figure 6a are given in Figure 6d,e measured from north to south. The backscattering coefficient data were obtained using Google Earth Engine [63]. The 10 m-resolution S1_GRD collection was used for computations. In each time interval (from May to October of each year), the data were averaged and stored as a single image, which was used to construct the profile. In total, there were at least 25 Sentinel-1 images in each time interval.
As we can see from Figure 6b As we can see from Figure 6b, d, the depth of the displaced mass reached 140-150 m. After the failure, the slope steepness of the main scarp was about 70° and approximately 20°-25° in the middle part of the surface of rupture. The slope steepness of the eastern scarp was about 22°, while the western scarp in the upper part is close to sheer and its steepness exceeds 70°.
Sentinel-1 backscattering coefficient measurement results, given in Figure 6c,e, allow us to assess the values averaged for the 2018 summer months (before the landslide event) and for 2019-2021 (after the event). The Sentinel-1B SAR acquiring images at a descending orbit with a heading angle of 192° from the north in the right-side looking observation mode, directed mainly to the west, with an incidence angle of 36.6°. For the western slope, after a landslide, both foreshortening and layover are observed, when the radar beam first reaches the upper part of the slope, since the slope angle, in this case, is greater than the incidence angle (Figure 6b). Reflection from a steep slope results in a maximum in the scattered signal, and this maximum is stable over time (Figure 6c). In contrast to the western flank, in the eastern flank, before failure, there was a slope facing the incident radiation (Figure 6b), which led to the maximum signal level in 2018 (red curve). The destruction of most parts of the slope did not cause this maximum decrease. Note, that backscattering from the surface not affected by the landslide remained stable. Variations of the scattered signal in the surface of the rupture after the slope failure in different years behave in a rather similar way (within 2-3 dB) and do not allow drawing a conclusion about interannual changes.
For clarity, Figure 7 also shows RGB composite image obtained using Google Earth Engine computing platform [63] from Sentinel-1B cross-polarized C-band radar data averaged over 2018, 2019, and 2020 (90 images in total). In the RGB image, red (R) represents the averaged values of radar backscatter in 2018, green (G) those for 2019, and blue (B) for 2020 (see inset). The qualitative interpretation is as follows. Red predominates in areas where backscattering decreased in 2019 and 2020 compared to 2018, i.e., in this case, there was a forest. Cyan, which is the result of mixing green and blue, shows an increase in backscattering in 2019-2020, i.e., landslide masses that filled the Bureya River. We can see red, cyan, and white colors on the surface of the rupture which is evidence of the complicated structure of that surface. Moreover, this qualitative interpretation indicates different mechanisms of backscattering transformation and small fragments remaining unchanged. Polarimetry and interferometry methods will be used further for the qualitative assessment of the surface of the rupture state and dynamics. Sentinel-1 backscattering coefficient measurement results, given in Figure 6c,e, allow us to assess the values averaged for the 2018 summer months (before the landslide event) and for 2019-2021 (after the event). The Sentinel-1B SAR acquiring images at a descending orbit with a heading angle of 192 • from the north in the right-side looking observation mode, directed mainly to the west, with an incidence angle of 36.6 • . For the western slope, after a landslide, both foreshortening and layover are observed, when the radar beam first reaches the upper part of the slope, since the slope angle, in this case, is greater than the incidence angle (Figure 6b). Reflection from a steep slope results in a maximum in the scattered signal, and this maximum is stable over time (Figure 6c). In contrast to the western flank, in the eastern flank, before failure, there was a slope facing the incident radiation (Figure 6b), which led to the maximum signal level in 2018 (red curve). The destruction of most parts of the slope did not cause this maximum decrease. Note, that backscattering from the surface not affected by the landslide remained stable. Variations of the scattered signal in the surface of the rupture after the slope failure in different years behave in a rather similar way (within 2-3 dB) and do not allow drawing a conclusion about interannual changes.
For clarity, Figure 7 also shows RGB composite image obtained using Google Earth Engine computing platform [63] from Sentinel-1B cross-polarized C-band radar data averaged over 2018, 2019, and 2020 (90 images in total). In the RGB image, red (R) represents the averaged values of radar backscatter in 2018, green (G) those for 2019, and blue (B) for 2020 (see inset). The qualitative interpretation is as follows. Red predominates in areas where backscattering decreased in 2019 and 2020 compared to 2018, i.e., in this case, there was a forest. Cyan, which is the result of mixing green and blue, shows an increase in backscattering in 2019-2020, i.e., landslide masses that filled the Bureya River. We can see red, cyan, and white colors on the surface of the rupture which is evidence of the complicated structure of that surface. Moreover, this qualitative interpretation indicates different mechanisms of backscattering transformation and small fragments remaining unchanged. Polarimetry and interferometry methods will be used further for the qualitative assessment of the surface of the rupture state and dynamics.

Polarimetric Decomposition
Earlier in [53], the authors estimated changes in the physical mechanisms of radar backscattering using fully polarimetric data of ALOS-2 PALSAR-2 SAR. The key attribute of the landslide zone was the mechanism of transformation of the dominant backscattering before and after the event. The polarimetric images of the landslide area were processed using the Freeman-Durden and Cloude-Pottier decompositions. The results obtained using the Freeman-Durden method [47] are shown in Figure 8a,b as color images representing different contributions of a double-bounce (Pd, red color), volume (Pv, green color) and a surface (Ps, blue color) scattering components. A red dashed line highlights the surface of the rupture, a pink solid line delineates landslide masses in the river. The black dash-dot line shows the part of the opposite shore where the forest was washed away by a tsunami wave. The site before the event was represented by forest cover with a relevant dominance of the volumetric scattering component (84.1%). After the landslide, despite the absence of trees, the volume scattering decreased by 13% but remained dominant. Such an effect was interpreted by the presence of chaotically oriented, as well as randomly located in space, scatterers (small-scale elevation inhomogeneities, piles of trees), acting like volumetric inhomogeneities. An increase of more than 13% in surface scattering means the appearance of rough surfaces with irregularities of the order of the wavelength after the event. Figure 9 shows the result of the H-α-A polarimetric decomposition technique applied to the same ALOS-2 PALSAR-2 polarimetric data set for the surface of the rupture, as described above.

Polarimetric Decomposition
Earlier in [53], the authors estimated changes in the physical mechanisms of radar backscattering using fully polarimetric data of ALOS-2 PALSAR-2 SAR. The key attribute of the landslide zone was the mechanism of transformation of the dominant backscattering before and after the event. The polarimetric images of the landslide area were processed using the Freeman-Durden and Cloude-Pottier decompositions. The results obtained using the Freeman-Durden method [47] are shown in Figure 8a,b as color images representing different contributions of a double-bounce (Pd, red color), volume (Pv, green color) and a surface (Ps, blue color) scattering components. A red dashed line highlights the surface of the rupture, a pink solid line delineates landslide masses in the river. The black dash-dot line shows the part of the opposite shore where the forest was washed away by a tsunami wave. The site before the event was represented by forest cover with a relevant dominance of the volumetric scattering component (84.1%). After the landslide, despite the absence of trees, the volume scattering decreased by 13% but remained dominant. Such an effect was interpreted by the presence of chaotically oriented, as well as randomly located in space, scatterers (small-scale elevation inhomogeneities, piles of trees), acting like volumetric inhomogeneities. An increase of more than 13% in surface scattering means the appearance of rough surfaces with irregularities of the order of the wavelength after the event. Figure 9 shows the result of the H-α-A polarimetric decomposition technique applied to the same ALOS-2 PALSAR-2 polarimetric data set for the surface of the rupture, as described above.
It is interesting that the upper (right) part of the point cloud, falling into zones 2 and 5, and characterizing the volumetric scattering, practically did not change. This can be explained by the fact that after the landslide many fallen trees remained, and there was a large amount of soil and fragments of rock, which determined the nature of multiple scattering, similar to the crowns of trees in a forest. Further, from Figure 9 it follows that after the landslide, as expected, the fraction (zone 6) of surface scattering increased. Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 21  It is interesting that the upper (right) part of the point cloud, falling into zones 2 and 5, and characterizing the volumetric scattering, practically did not change. This can be explained by the fact that after the landslide many fallen trees remained, and there was a large amount of soil and fragments of rock, which determined the nature of multiple scattering, similar to the crowns of trees in a forest. Further, from Figure 9 it follows that after the landslide, as expected, the fraction (zone 6) of surface scattering increased.
The comparison of Freeman-Durden and Cloude-Pottier decompositions have shown their general agreement, however, some regions showed significant differences due to the specifics of the basic model assumptions. Overall, it was shown that before the landslide (28 November 2018), there were three main mechanisms of radar signal scattering in the landslide area: single surface scattering, volumetric scattering, and double scattering. After the landslide (12 December 2018), the single scattering component increased in this area, which is typical for a reflecting surface without vegetation, so the landslide zone can be confidently recognized.
Like with the results obtained using completely polarimetric ALOS-2 PALSAR-2 data, we use the H-α-A polarimetric decomposition technique to analyze backscatter mechanism variations after the landslide in 2019-2020 based on dual-polarization Senti-  It is interesting that the upper (right) part of the point cloud, falling into zones 2 and 5, and characterizing the volumetric scattering, practically did not change. This can be explained by the fact that after the landslide many fallen trees remained, and there was a large amount of soil and fragments of rock, which determined the nature of multiple scattering, similar to the crowns of trees in a forest. Further, from Figure 9 it follows that after the landslide, as expected, the fraction (zone 6) of surface scattering increased.
The comparison of Freeman-Durden and Cloude-Pottier decompositions have shown their general agreement, however, some regions showed significant differences due to the specifics of the basic model assumptions. Overall, it was shown that before the landslide (28 November 2018), there were three main mechanisms of radar signal scattering in the landslide area: single surface scattering, volumetric scattering, and double scattering. After the landslide (12 December 2018), the single scattering component increased in this area, which is typical for a reflecting surface without vegetation, so the landslide zone can be confidently recognized.
Like with the results obtained using completely polarimetric ALOS-2 PALSAR-2 data, we use the H-α-A polarimetric decomposition technique to analyze backscatter mechanism variations after the landslide in 2019-2020 based on dual-polarization Senti- The comparison of Freeman-Durden and Cloude-Pottier decompositions have shown their general agreement, however, some regions showed significant differences due to the specifics of the basic model assumptions. Overall, it was shown that before the landslide (28 November 2018), there were three main mechanisms of radar signal scattering in the landslide area: single surface scattering, volumetric scattering, and double scattering. After the landslide (12 December 2018), the single scattering component increased in this area, which is typical for a reflecting surface without vegetation, so the landslide zone can be confidently recognized.
Like with the results obtained using completely polarimetric ALOS-2 PALSAR-2 data, we use the H-α-A polarimetric decomposition technique to analyze backscatter mechanism variations after the landslide in 2019-2020 based on dual-polarization Sentinel-1 data. In Figure 10a,b for comparison, the results of H-α-A polarimetric decomposition data are presented in pairs for the entire surface of the rupture for 15 May and 31 August 2019, as well as for 15 May 2019 and 15 August 2020. It can be seen that the scattering mechanisms remained practically unchanged, which is consistent with the similarity of the amplitude profiles along the surface of the rupture for the snowless periods of 2019-2021. This generally reflects the fact that the change in inhomogeneities in the surface of the rupture retains the overall roughness for the C-band due to the roughness of the rock and, apparently, is associated with the insufficient ability of the combination of polarizations VV and VH (HH and HV) to distinguish between scattering mechanisms [72].
nel-1 data. In Figure 10a, b for comparison, the results of H-α-A polarimetric decom tion data are presented in pairs for the entire surface of the rupture for 15 May an August 2019, as well as for 15 May 2019 and 15 August 2020. It can be seen that the tering mechanisms remained practically unchanged, which is consistent with the sim ity of the amplitude profiles along the surface of the rupture for the snowless perio 2019-2021. This generally reflects the fact that the change in inhomogeneities in the face of the rupture retains the overall roughness for the C-band due to the roughne the rock and, apparently, is associated with the insufficient ability of the combinati polarizations VV and VH (HH and HV) to distinguish between scattering mechan [72].

PSI and SBAS-InSAR
To assess possible landslide reactivation using multitemporal C-band Sentin data we used PSI and SBAS-InSAR methods. To exclude a possible impact of sno differential interferometry results [73], SAR data for 3 May-6 October 2019 (14 scen May-12 October 2020 (13 scenes) and 4 May-25 September 2021 (12 scenes) were use co-polarized (VV) images in total. Figure 11a presents spatial segmentation of PS distribution superimposed on the tinel-2 (1 June 2021) image: right Bureya River bank, landslide body in the river and face of the rupture. The given data show that the greatest PS number (blue color in F 11a) was on the right Bureya River's bank and along the shores of the Sredny Sandar where the tsunami wave had passed in the area of the destroyed forest.

PSI and SBAS-InSAR
To assess possible landslide reactivation using multitemporal C-band Sentinel-1B data we used PSI and SBAS-InSAR methods. To exclude a possible impact of snow on differential interferometry results [73], SAR data for 3 May-6 October 2019 (14 scenes), 9 May-12 October 2020 (13 scenes) and 4 May-25 September 2021 (12 scenes) were used; 39 co-polarized (VV) images in total. Figure 11a presents spatial segmentation of PS distribution superimposed on the Sentinel-2 (1 June 2021) image: right Bureya River bank, landslide body in the river and surface of the rupture. The given data show that the greatest PS number (blue color in Figure 11a) was on the right Bureya River's bank and along the shores of the Sredny Sandar River where the tsunami wave had passed in the area of the destroyed forest. A small number of PSs is observed in the surface of the rupture (red color in Figure  11a). This may indicate a small number of natural point scatterers and/or permanent ground movements in the surface of the rupture, the displacement velocities of which exceed the velocities that can be discovered by the PSI method. The latter assumption is confirmed by the low coherence of interferometric pairs with a time baseline greater than 12 days. In particular, the application of the MAI InSAR method [74] showed no interferometric coherence at long time intervals (more than 12-24 days), and at 12-24 day time A small number of PSs is observed in the surface of the rupture (red color in Figure 11a). This may indicate a small number of natural point scatterers and/or permanent ground movements in the surface of the rupture, the displacement velocities of which exceed the velocities that can be discovered by the PSI method. The latter assumption is confirmed by the low coherence of interferometric pairs with a time baseline greater than 12 days. In particular, the application of the MAI InSAR method [74] showed no interferometric coherence at long time intervals (more than 12-24 days), and at 12-24 day time baselines, the shifts were indistinguishable. The PSs within the surface of the rupture are located at a distance from its flanks and, apparently, make a translatory movement in the longitudinal direction of the landslide. Figure 11b presents the displacement values for three PS segments highlighted with blue, yellow (landslide body in the river), and red colors, respectively. A digital elevation model (DEM) obtained from TerraSAR-X/TanDEM-X radar images on 28 September 2020 was used as a reference. The spatial baseline was 130.3 m, 2π ambiguity in height-45.7 m.
With a general negative displacement of all PSs in a given area of the terrain, a high rate of displacement from the radar along LOS was observed for PSs in the surface of the rupture, reaching 40 mm/year.
We assume that PSs were moving along the direction of the fastest landsliding with the values ϕ ≈ 22º and γ ≈ 20º obtained from DEM. Thus, from (4) and (5) equations we obtain cos β = 0.23. Note, that due to the proximity of the direction of the fastest sliding to the "south-north" line, the value of cos β turned out to be less than the threshold value of 0.3, stated in [75] to eliminate anomalous situations. Using this threshold value, we obtain the following: Thus, the PSI method shows that after the landslide, there was a noticeable displacement velocity of discrete PS along the surface of the rupture.
PSI assessment for a snowless period of 2019-2021 has shown that there were few persistent scatterers on the surface of the rupture (16 PSs were grouped in the center), so the reliability of the results with their use may raise some doubts. For a more detailed analysis of the situation, the deformations were calculated using the SBAS-InSAR method ( Figure 12) according to Sentinel-1B data separately for snowless periods of time in 2019, 2020, and 2021.
Distributed scatterers' location and velocity (mm/year) of displacement over the surface with a coherence of more than 0.2 calculated using the SBAS-InSAR method are given in Figure 12. The right part of this figure presents the temporal and spatial baseline distributions of the radar interferograms from the Sentinel-1B data sets; the green squares represent the valid acquisitions, the red one marks the selected super master image of the small baseline subset and the solid lines connects interferometric pairs for SBAS-InSAR. The results show a general increase of the area of regions with a coherence of more than 0.2 over time for snowless periods in 2019, 2020, and 2021. This indicates a decrease in surface changes and significant deformations along the surface the rupture in the period under consideration. The spatial distribution of the areas of the largest (red-and-yellow) displacements along the LOS has changed. In 2019, the displacements were closer to the center; in 2020 and 2021, deformations are noted at the edge of the lateral eastern edge of slope failure. Max displacement velocity reduced from 40 mm/year in 2019 (red color in Figure 12a) down to 30 mm/year in 2020-2021 (yellowish-red in Figure 12b Distributed scatterers' location and velocity (mm/year) of displacement over the surface with a coherence of more than 0.2 calculated using the SBAS-InSAR method are given in Figure 12. The right part of this figure presents the temporal and spatial baseline distributions of the radar interferograms from the Sentinel-1B data sets; the green squares rep- In general, the obtained SBAS-InSAR results show, on the one hand, continuing displacements along the surface of the rupture, on the other hand, some stabilization (decrease) in the rate of landsliding processes. Distributed scatterers, as noted, are located on the eastern slope of the surface of the rupture. From the DEM analysis, it follows that for this slope, the following estimates for the geometric parameters can be taken approximately: ϕ ≈ 22º and γ ≈ 250º. Therefore, it follows from the Equations (4) and (5): Comparison with PSI velocities shows that the displacement movement velocity along the main slope was 2.5-fold higher than the displacement movement velocity along the lateral eastern slope, determined by the SBAS-InSAR method. Since the difference in displacement velocities was observed on different slopes, we can suppose that this difference is because of the structures of each of the slopes. This assumption is supported by the fact that persistent (point) scatterers are observed only on the main slope and distributed (areal) scatterers are observed on the eastern slope.

Conclusions
The capabilities of PSI and SBAS-InSAR methods, polarimetric decomposition techniques, and SAR backscattering multitemporal series estimation have been demonstrated for the case study of the landslide reactivation on the Bureya River in Russia. Freely available Sentinel-1B SAR interferometric data for 2019-2021 snowless periods were used [67]. Multitemporal series of backscattering data have been analyzed with help of cloud computing [61] taking into account radar sensing geometry and relief elevation.
It is shown that the spatial variations of the backscattered signal in the surface of the rupture in different years behave in a rather similar way (within 2-3 dB) and do not allow make a conclusion about interannual changes. The results of polarimetric analysis over the entire surface of the rupture based on dual polarization C-band data showed that the scattering mechanisms remained practically unchanged. This is consistent with the similarity of the amplitude profiles along the surface of the ruptures for the snowless periods of 2019-2021. Overall, both results show that the change in the inhomogeneities in the surface of the rupture retains the overall roughness comparable with the C-band wavelength. Apparently, this is also due to the insufficient ability of the combination of VV and VH polarizations to distinguish scattering mechanisms [72]. At the same time, a qualitative interpretation of false-color composites formed from year-averaged SAR backscattering images indicates different mechanisms of backscattering transformation and small fragments without changes.
A quantitative assessment of the total changes for the snowless periods of 2019, 2020 and 2021, carried out using the PSI method, showed an overall negative displacement of all PSs for the studied area. A high rate of displacement along the radar LOS is observed for persistent scatterers in the surface of the rupture, reaching 40 mm/year. At the same time, the difference in the magnitude of the displacements of the main group of persistent scatterers on the opposite bank and PS in the surface of the rupture in 2021 reached 20 mm.
The results of the Small baseline subsets (SBAS-InSAR) method made it possible to identify the overall increase in the area of regions with a coherence of more than 0.2 in the period from 2019 to 2021. This indicates a decrease in the processes of surface changes and the absence of significant deformations in this area of the surface of the rupture during the period under consideration. The spatial distribution of the regions of greatest displacement along the LOS changed. It was found that in 2019 the displacements were closer to the center, and in 2020 and 2021, deformations at the edge of the lateral eastern edge of slope failure were revealed. Max displacement velocity reduced from 40 mm/year in 2019 down to 30 mm/year in 2020-2021. PSI method results have shown that the displacement velocity along the main slope exceeded the displacement velocity along the eastern slope 2.5-fold determined by the SBAS-InSAR method.
Thus, the studies carried out have demonstrated the effectiveness of the combined application of MT-InSAR methods for identifying the processes of landslide reactivation. At the same time, the PSI method allows one to estimate the total displacements of individual persistent scatterers over several years in snowless periods of time, and the SBAS-InSAR method allows one to obtain estimates of displacements over large areas, but at short time intervals.