Deformation Retrievals for North America and Eurasia from Sentinel-1 DInSAR: Big Data Approach, Processing Methodology and Challenges

Abstract A fully automated processing system for measuring long-term ground deformation time series and deformation rates frame-by-frame using DInSAR processing technique was developed at the Canada Center for Remote Sensing. Ground deformation rates from 2017 to 2023 were computed over a large territory of North America and Eurasia from more than 220,000 readily available Sentinel-1 images, and the performance and shortcomings of the developed processing system were analyzed. Here, we present the processing methodology and several examples of deformation rate maps and time series produced with this automated system. Examples include the deformation of slow- moving deep-seated landslides in two regions of Canada, subsidence at the Komsomolskoe oil field in the Russian Arctic, the Tengiz oil field in Kazakhstan, multiple large subsiding regions and landslides in northwestern Iran, and two large subsiding regions in the Yellow River Delta and Xinjiang, China. Many deformation processes observed in these deformation rate maps, including large landslides, have previously been unknown to the research community. Systematic radar penetration depth changes were observed in multiple regions and were investigate in detail for 1 Eurasian region. Computed deformation rates for North America and Eurasia are available to the research community and can be downloaded from the data repository.


Introduction
The recent improvement in the availability of Synthetic Aperture Radar (SAR) data provides an opportunity for measuring ground deformation with high temporal and spatial resolution on a global scale.Processing techniques and systems were previously developed assuming data scarcity; these techniques are suboptimal for processing large SAR datasets containing hundreds of images (Zinno et al. 2016;Ronczyk et al. 2022;Samsonov 2022).We developed a prototype of a fully automated processing system suitable for processing large SAR datasets on High-Performance Computer (HPC) Linux clusters (Dudley andSamsonov 2020, 2021).Adopting a Big Data approach that aims to detect as many active ground deformation processes in complex environments as possible sometimes results in the tradeoff of reduced spatial and temporal resolutions and precision.To better understand the performance as well as shortcomings of the system, we computed ground deformation over a large area of North America and Eurasia; here we present results for a select subset of the study area.We observed a systematic noise source that limits measurement precision in multiple regions.
Differential Interferometric Synthetic Aperture Radar (DInSAR) is a processing technique for measuring ground deformation from 2 SAR images acquired by the same or similar satellites (Rodriguez and Martin 1992;Massonnet and Feigl 1998;Rosen et al. 2000).Given that the DInSAR interferogram captures a projection of the 3-dimensional (3D) ground deformation vector on the satellite Line-of-Sight (LOS) vector, having SAR data acquired on ascending and descending orbits allows us to observe ground deformation from different perspectives.Horizontal and vertical components of the 3D deformation vector can be computed from ascending and descending interferograms, speckle offsets (Manzo et al. 2006;Pepe and Calo 2017), e.g., Digital Elevation Model (DEM) assuming the Surface-Parallel Flow (SPF) constraint on deformation (Joughin et al. 1998;Mohr et al. 1998;Kumar et al. 2011;Samsonov 2019;Samsonov et al. 2020).
The time series of ground deformation is computed from multiple interferograms created from repeatedly acquired SAR data.The linear deformation rate is calculated by applying linear regression to the time series.
One needs to show cumulative deformation maps for each acquisition epoch to visualize the entire deformation history, which can be impractical for long time series.Instead, the linear deformation rate maps are used for visualizing an average ground deformation rate over the entire area, and time series are shown only for individual pixels of interest.Various techniques exist for computing deformation time series (Osmano glu et al. 2016), with the Small Baseline Subset (SBAS) technique being the most popular due to its conceptual simplicity (Lanari et al. 2004;Casu et al. 2006;Lanari et al. 2007;Sansosti et al. 2010;Cigna et al. 2011).Horizontal and vertical deformation time series can be computed from multiple ascending and descending interferograms, using, for example, the Multidimensional Small Baseline Subset (MSBAS) technique (Samsonov and d'Oreye 2012;Samsonov et al. 2013;Samsonov and d'Oreye 2017).
The C-band Sentinel-1 satellite constellation consists of 2 satellites launched in 2014 (Sentinel-1A) and 2016 (Sentinel-1B), with a third satellite (Sentinel-1C) expected to launch in 2023.The Sentinel-1B satellite has not been operational since December 2021.For many parts of the world, regularly acquired SAR data with an average temporal resolution of 12 days, sometimes on ascending and descending orbits, has been available since 2016/2017 and over Europe since 2014.Sentinel-1's Interferometric Wide Swath mode (IWS) with 250 km swath and 4 Â 14 m spatial resolution is implemented with a new ScanSAR mode called Terrain Observation with Progressive Scan (TOPS).Initial challenges in the interferometric processing of TOPS data were rapidly overcome, allowing Sentinel-1 to become the most utilized source of SAR data globally (Yague-Martinez et al. 2016).
Our first attempts to measure the long-term ground deformation in British Columbia and Alberta, Canada, revealed numerous challenges both in computation and interpretation.We successfully detected ground deformation produced by natural (e.g., landslides, glacier flow) (Choe et al. 2021a) and anthropogenic (e.g., caused by fracking, extraction of oil and groundwater, and mining) processes (Eyre et al. 2022).At the same time, we observed many signals later classified as processing artifacts that reduced measurement precision, making the interpretation challenging.To better understand the causes of these artifacts, we tested our processing technique outside Canada, in Eurasia.Southern Eurasia has many areas that experience anthropogenic ground deformation, so it provides a chance to examine rapidly changing deformation processes occurring over a short spatial scale.Northern Eurasia's land cover and climate are similar to Canada's, so it offers an opportunity to study the impact of precipitation and temperature seasonality.We then extended spatial coverage to test various computational aspects of our software, specifically the effectiveness of our parallelization routines and the hardware capacity of our HPC system, such as storage bandwidth.
Here, we describe our processing approach and discuss various challenges that we encountered.Deformation maps over these regions are provided in the data repository (Samsonov and Feng 2022); a select subset of interesting deformation signals is presented and discussed here.The results are presented frame by frame; methodologies for seamless mosaicking are not addressed.Our objectives are to stimulate discussion on developing an optimal fully automated processing system, improve the measurement quality of DInSAR at regional and global scales, and create a database of active deformation processes for regular operational monitoring.

Methodology
A processing system is designed to achieve specific goals.The objective of our fully-automated processing system is to detect as many deformation processes over as wide (i.e., 250 km swath of Sentinel-1) an area as possible, one frame at a time.We accept that some deformation processes occurring in challenging DInSAR environments may be measured with a reduced precision due to low coherence or other factors.After an active deformation process is detected, the processing strategy can be modified to improve that specific process's measurement precision, but this is beyond the scope of this study.
The processing system is designed to run on HPC Linux clusters containing 1 or more computer nodes, and the GAMMA software (Wegmuller and Werner 1997) is used for computing interferograms (Dudley andSamsonov 2020, 2021).A node is similar to an individual computer with shared memory space and several CPU cores.The nodes used by our system contain 128 GB of RAM and two 20-cores CPUs.A typical job uses 10 nodes (i.e., 400 cores) to produce deformation time series from over a 100 Sentinel-1 SLC images of a specific frame and path (Figure 1).
A job is subdivided into tasks that use 8 cores each to generate a deformation map from 2 SAR images, i.e., an unwrapped and geocoded interferogram with LOS displacement measured in meters.Below we use the 2 terms, deformation map and interferogram, interchangeably.At most, 50 tasks (i.e., 400 available cores/8 cores per task) are run in parallel by a single job.Our tests show that using more than 8 cores per task to compute a single deformation map does not noticeably reduce processing time.Since 5 tasks (i.e., 40 cores per node/8 cores per task) can run in parallel on each node, the maximal amount of RAM available for each task is about 25 GB (i.e., 128 GB RAM per node/5 tasks per node).After all deformation maps are computed, the MSBAS version 10 software computes the deformation time series.The MSBAS software, coded in Cþþ, uses Open Multi-Processing (OpenMP) and Message Passing Interface (MPI) parallelization methods to distribute the processing between the allocated nodes.OpenMP allows the application to use multiple cores of the same node, and MPI allows using several nodes.Our HPC cluster contains over 100 nodes, so at least 10 jobs can be run simultaneously.
Slow-moving deep-seated landslides are one of the most common deformation processes observed with DInSAR (Choe et al. 2021a;Samsonov et al. 2020;Samsonov and Blais-Stevens 2023).The spatial extent (e.g., length) of such landslides can range from tens of meters to kilometers, so the spatial resolution of deformation products must be sufficiently high to capture small landslides.The spatial resolution of deformation products is controlled by the spatial resolution of SAR imagery and multilooking, which is applied during the interferometric processing.Multilooking improves measurement precision by averaging phase values at nearby pixels and reducing the interferogram's size, making phase unwrapping less computationally intensive.When processing Sentinel-1 TOPS data with about 4 Â 14 m spatial resolution, we apply 12 Â 3 multilooking resulting in about 50 Â 50 m ground resolution.Such multilooking produces interferograms with approximately 6,000 Â 6,000 pixels that are too large for efficient phase unwrapping.Unwrapping large interferograms is slow and consumes a large amount of RAM.In our system, each 8-core task has only 25 GB of RAM, which is insufficient for unwrapping large interferograms.Therefore, before phase unwrapping, we subsample the interferograms to approximately 2,000 Â 2,000 pixels and then, after coarse unwrapping, oversample them to the original resolution by copying values from the wrapped interferograms after the 2p ambiguity was resolved.Such an approach significantly reduces processing time and required RAM while sometimes missing minor deformation processes.Alternatively, it is possible to divide an interferogram into several parts that can be unwrapped separately, but this often introduces errors when separately unwrapped interferogram parts are recombined.
The topographic phase is removed using the 30 m resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Digital Elevation Model (DEM) (Abrams et al. 2020).Differential interferograms are filtered using adaptive filtering with a filtering function based on the local fringe spectrum (Goldstein and Werner 1998) and unwrapped using the minimum cost flow algorithm (Costantini 1998).We use precise orbital information and do not apply an orbital correction.For the time series analysis with MSBAS, we resample all the deformation maps to a common lat/long grid with the Geospatial Data Abstraction Library (GDAL, https://gdal.org)gdalwarp program.
The MSBAS version 10 software is an improved version of the MSBAS versions 3 and 6 software, previously used in various studies (Samsonov and d'Oreye 2012;Samsonov et al. 2013;Samsonov and d'Oreye 2017) now modified to run on HPC Linux clusters.It uses the singular value decomposition (SVD) algorithm from the Linear Algebra PACKage (LAPACK, https://www.netlib.org/lapack)library and Tikhonov regularization (Tikhonov and Arsenin 1977) for inverting the time matrix for each pixel.The explicit definition of the time matrix can be found in (Samsonov and d'Oreye 2017;Samsonov et al. 2021a).This study uses the first order Tikhonov regularization that minimizes differences between consecutive deformation rates, and overall, we are satisfied with using the first order rather than zeroth or second order.
We compute 3 deformation maps using each SAR image as primary and 3 immediately acquired images as secondary.Building an optimal network from these deformation maps for inversion is not a simple task (Smittarello et al. 2022) because each deformation map has unique spatial coverage defined by its coherence, i.e., pixels with coherence below a particular threshold are masked out.We use a low coherence threshold of 0.4 to avoid missing potential deformation signals and accept that some deformation maps may contain incorrectly unwrapped areas, reducing the overall precision of the final deformation products.Spatial coverage can also be affected by the different number of bursts in SAR images of the same frame.
To select deformation maps for inversion, we assign to each ith deformation map a weight based on its coherence and temporal baseline equal to c Dt i i , where Dt i is the temporal baseline, and c i is the average coherence for the entire deformation map.Then, we select deformation maps connecting the first and the last SAR images that maximize Q c Dt i i : A more sophisticated weighting scheme of c Dt a i i with a < 1 could be applied for penalizing interferograms with small temporal baselines, but it is not used here.
Practically, the network is built using the igraph package (https://igraph.org)implemented in R; Python and C implementations of this package are also available.For this, we assume that each deformation map is an edge of a graph connecting 2 vertices representing SAR images.The distance (or weight in graph theory) between 2 vertices is equal to ÀDt i ln c i : The minus sign is introduced for reformulating the problem in terms of finding the shortest path (i.e., minimization problem) of the graph, a well-researched problem in graph theory.The overall distance between the first and the last vertices of the graph (i.e., first and last acquired SAR images) of a specific path is then P ÀDt i ln c i , and many paths connecting the first and last vertices can be constructed.Selecting deformation maps for inversion is then reduced to selecting the shortest path of the constructed graph, implemented in the igraph package.
Additionally, we compute the average weight of all selected deformation maps and add to the network those previously unselected deformation maps whose weight is lower than the average, which allows us to introduce some redundancy.
Theoretically, an optimal network should be built for each pixel separately using coherence values at that pixel, but this process would need to be repeated 36 M times for an area of 6,000 Â 6,000 pixels, which is prohibitively slow.Instead, we do it only once using an average for the coherence values of the deformation map and then allow MSBAS to handle the inversion of partially coherent pixels.MSBAS records the time matrix rank at each pixel before adding the regularization terms.If the network at that pixel is fully connected, the time matrix rank is equal to the number of SAR images minus 1 (Samsonov and d'Oreye 2017).When 2 consecutive SAR images are unconnected in the network, the time matrix rank decreases by 1.Therefore, by analyzing the rank file created by MSBAS, we can determine the overall network quality at each pixel.Potentially, we can use this methodology to allow the MSBAS software to build its network at each pixel, skipping the graph's shortest path approach entirely.Practically, we still select displacement maps using the graph's shortest path approach because we want to limit the quality and the number of displacement maps provided to MSBAS for analysis to reduce processing time and then allow MSBAS to refine the network pixel by pixel.To improve spatial coverage reduced by decorrelation, we call the GDAL GDALFillNodata function from MSBAS to fill gaps in the deformation maps by interpolating within a small radius.Upon MSBAS completion, we retain only pixels with a rank above the selected threshold value, chosen as half of the maximum value of the entire rank file.Significantly, this process is performed after completing MSBAS processing, which is time-and resource-consuming and does not need to be repeated for testing different threshold values (Samsonov et al. 2021a).
Figure 2 helps understand how MSBAS uses the time matrix rank for assessing network quality.For example, assume that there are 3 SAR images acquired at times t 1 , t 2 , t 3 and 3 interferograms I 1 , I 2 , I 3 are computed from these images, spanning Dt 12 , Dt 23 , Dt 12 þ Dt 23 : Note that the distribution of coherent pixels in each interferogram can vary.All 3 interferograms are supplied to MSBAS with the objective of estimating the instantaneous velocities V 12 , V 23 so that the cumulative deformation time series 0, V 12 Dt 12 , V 12 Dt 12 þ V 23 Dt 23 can be calculated at each pixel.If at an arbitrary pixel, all 3 interferograms contain coherent values (Figure 2b), then the time matrix rank is 2 (i.e., 3 SAR images minus 1), and the system of the equations is overdetermined (i.e., 3 equations and 2 unknowns).If at another pixel (any of 3 possibilities in Figures 2b and  2c), only 2 interferograms have coherent values, then the time matrix rank is also 2, and the system of the equations is well-posed (i.e., 2 equations and 2 unknowns).If at another pixel (any of 3 possibilities in Figures 2f, 2g, and 2h) only 1 interferogram has a coherent pixel, then the time matrix rank is 1, and the system of the equations is underdetermined (or illposed).In this case, Tikhonov regularization comes into play by setting V 12 equal to V 23 or vice versa.Therefore, computed with MSBAS deformation time series and deformation rate contain valid values at all these pixels.In post-processing, we can select pixels with the desired network quality by analyzing the produced rank file at each pixel.In the previous versions of MSBAS, having just 1 incoherent pixel would have resulted in excluding this pixel from the time series, significantly limiting overall spatial coverage.
When the network has 1 or more discontinuities and the time matrix rank is reduced, the first order Tikhonov regularization interpolates the gaps in the time series in the temporal domain.We find that this is a valuable behavior in many cases, and it rarely produces the undesirable effects described below.
In addition to standard DInSAR time series analysis, our fully automated processing system can produce azimuth and range deformation time series computed with speckle offset tracking (SPO) (Samsonov 2022).With minimal manual intervention, combining ascending and descending DInSAR or range/azimuth SPO data, the system can compute 2D-4D deformation time series using advanced processing methods implemented in MSBAS (Samsonov et al. 2021a(Samsonov et al. , 2021b) ) and which are not discussed in this study.

Data
Sentinel-1 data in Single Look Complex (SLC) format is downloaded from the Alaska Satellite Facility (ASF, https://asf.alaska.edu) 1 frame at a time, and only frames with more than 100 SAR images are selected.The system can potentially concatenate images from the adjacent frames, but in the case of Sentinel-1, it is suboptimal because the number of available images per frame along the same path can vary significantly.We use the 30 m resolution ASTER DEM for removing the topographical phase.This DEM is chosen because it provides uniform world coverage, but it is possible to use different DEMs, including SRTM (Farr and Kobrick 2000), ArcticDEM (Showstack 2017) or custom LiDAR-derived.The system can fill gaps in DEM over water bodies with the geoid values.This option is beneficial when computing deformation maps of coastal sea/river ice (Choe et al. 2021b) or water-terminating glaciers.Precise Sentinel-1 orbit files are also downloaded from the ASF.

Results
We processed more than 220,000 ascending and descending SAR images from about 2,200 frames (Table 1, Figure 3).Here, we present processing results that include deformation rate maps and deformation time series for 8 select regions.A single Sentinel-1 frame covers each region.Active ground deformation is observed in 7 regions, while significant processing artifacts are observed in 1 region.This region is described in detail as similar artifacts are observed in about half of all the computed deformation rate maps in the repository.
Figure 4a shows a deformation rate map of a region in western Canada computed from 149 descending Sentinel-1 images from path 42 frame 391.In this and other figures, unless specified otherwise, the color scale is clipped in a range [À0.1;0.1]m/year for consistency.The insert in the top-left corner shows the location of this region in the world, and the insert in the top-right corner shows deformation rate contour lines plotted over an optical image from Google Earth.The solid contour lines range from 0.06 (red) to 0.12 (purple) m/year, and the dashed contour lines range from À0.05 (purple) to À0.02 (red) m/year.The observed deformation signal is due to 2 previously unknown slow-moving deep-seated landslides, one of which is likely the largest in Canada.The larger landslide moves toward the east (i.e., positive LOS displacements), and the smaller landslide moves toward the west (i.e., negative LOS displacements).The network of deformation maps used for the inversion is shown in Figure 4b, and the deformation time series for the 5 Â 5 pixel area marked with a red star in the top-right insert is shown in Figure 4c.Note that the multiple regions of apparent subsidence in Figure 4a that are processing artifacts are discussed below.
Figure 5a shows a deformation rate map of a region in northern Canada computed from 128 ascending Sentinel-1 images from path 108 frame 221.The insert in the top-right corner shows the location of this region in the world, and another insert shows deformation rate contour lines plotted over the Sentinel-2 image.The solid contour lines range from À0.03 (red) to À0.21 (purple) m/year.The observed deformation signal is another large, previously unknown slow-moving deep-seated landslide in Canada.The sliding motion is toward the east (i.e., negative LOS displacements).The network of deformation maps used for the inversion is shown in Figure 5b, and the deformation time series for the 5 Â 5 pixel area marked with a red star in the top-right insert is shown in Figure 5c.
Figure 6a shows a deformation rate map of a region in Russia computed from 117 descending Sentinel-1 images from path 151 frame 377.As in the previous figures, the insert in the top-left corner shows the location of this region in the world, and the insert in the top-right corner shows deformation rate contour lines plotted over an optical image from Google Earth.The solid contour lines range from À0.24 (purple) to À0.07 (red) m/year.The observed deformation signal is believed to be subsidence at the Komsomolskoe oil field.It is over 10 km in diameter and affects the mining and transportation infrastructure.The network of deformation maps used for the inversion is shown in Figure 6b, and the deformation time series for the 5 Â 5 pixel area marked with a red star in the top-right insert is shown in Figure 6c.
Figure 7a shows a deformation rate map of a region in China computed from 206 Sentinel-1 images from path 69 frame 119, and Figure 7b shows a deformation rate map of another region in China computed from 147 ascending Sentinel-1 images from path 41 frame 140.The insert shows locations of the region in the world, the network of deformation maps used for the inversion, and the deformation time series for the 5 Â 5 pixel area pointed to with black line(s).The Yellow River Delta (Figure 7a) is one of the largest river deltas in China, where a population of over 5 million people live in a region with an area of 18,000 sq.km.Due to increasing anthropogenic activities, including petroleum extraction, significant subsidence has long been observed (Zhang et al. 2015(Zhang et al. , 2019;;Liu et al. 2021).Specifically, Zhang et al. ( 2019) used ENVISAT and Sentinel-1 SAR data to retrieve ground deformation over 2 periods; the latest deformation rate map computed from 2015 to 2018 Sentinel-1 data revealed a similar deformation pattern to that produced by our automated system.The ground deformation observed in Xinjiang (Figure 7b) is likely related to extensive farming in this region.
Figure 8a shows a deformation rate map of a region in northwestern Iran computed from 181 Sentinel-1 images from path 79 frame 464.Several large subsiding areas with a maximum deformation rate of about 0.15 m/year are observed.The subsidence is likely due to the extraction of groundwater.Multiple landslides, labeled with the letter L, are also observed.Similar landslide signals are observed in many areas around the world.
Finally, Figure 8b shows a deformation rate map of a region in Kazakhstan computed from 135 Sentinel-1 images from path 35 frame 441.A large, over 20 km in diameter, subsidence is observed at the Tengiz oil field.Note that different color scale ranges are used in these figures.
Additional deformation signals over the studied areas can be observed in the deformation maps provided in the data repository (Samsonov and Feng 2022).

Processing artifacts
Figure 9a shows a deformation rate map computed from 134 descending Sentinel-1 images in Siberia from path 149 frame 398.The inserts show locations of the region in the world, the network of deformation maps used for the inversion, and the deformation time series for the 5 Â 5 pixel area pointed to with a black line.The deformation time series in the insert show apparent subsidence; similar seasonal behavior is observed each year, with moderate apparent uplift from summer to early fall followed by apparent rapid subsidence from late fall to late spring.The apparent seasonal subsidence amplitude increases with elevation, resembling a topography-correlated pattern.
For testing purposes, we recompute time series using the previously used deformation maps and a few newly computed interannual deformation maps with approximately 1-year span that are partially coherent.We accept only pixels with a coherence of 0.8 and greater, and do not apply Tikhonov regularization.This processing produces deformation time series with significantly reduces spatial and temporal coverages.The deformation time series are then extracted at 2 regions close to point P from Figure 9a.The time series of the 5 Â 5 pixel area computed from the disconnected network (rank reduced by 2) is shown in Figure 9b.In this region, the interannual deformation maps and those spanning from late May to early June (red lines in the figure) do not have coherent pixels.The time series computed from the fully connected network is shown in Figure 9c.In this time series, the signal from late May to early June is reconstructed from the coherent at those pixels' interannual deformation maps.The resolved apparent uplift-like deformation signal is approximately equal to the total subsidence of the preceding months.Unfortunately, interannual deformation maps have low overall coherence and cannot be reliably used to reconstruct the seasonal cycle over the entire area.However, we expect that the signal observed in the time series in Figure 9c reflects the actual pattern we would have observed in the entire region if interannual deformation maps were fully coherent.
Because such a signal is observed in many deformation rate maps worldwide and individual deformation maps with temporal baselines of about 1 year (from summer-to-summer next year, not used in inversion) do not show any deformation, the presence of active deformation processes is ruled out.Likely, a similar process is responsible for the apparent subsidence observed in Figure 4 in the northern part of the map and for broader subsidence in the center of the map in Figure 6.We believe the seasonal deformation rate variability observed in Figures 7b and 8a is the actual subsidence, but it needs further validation.

Discussion and conclusions
A fully automated system was developed for computing from SAR data long-term ground deformation time series and deformation rates using the DInSAR processing methodology.This system was tested on over 220,000 Sentinel-1 images acquired over North America and Eurasia, and the results were analyzed.The average precision of the deformation rates was 0.018 m/year.Many novel deformation signals were detected, which will be addressed in detail in future studies.Several of these signals are shown and discussed.We were able to resolve many ground deformation signals in densely vegetated areas and areas that experience significant seasonal changes (e.g., Figures 4  and 6) by utilizing displacement maps acquired over short temporal baselines, often only 6-12 days.The most commonly observed deformation signals were landslides (e.g., Figures 4, 5, and 8a) and subsidence caused by natural and anthropogenic processes, such as mining and oil/gas/groundwater extraction (e.g., Figures 6, 7, and 8b).Fast glacial flow almost always was not resolved, so glaciers appear incoherent in our results.Coherent deformation signals were sometimes observed over frozen small water bodies.Validation of the observed deformation signals is beyond the scope of this study due to the vast amount of data in the repository; however, in-depth analysis of ground deformation detected with this system over a remote region in Canada where many landslides have been observed near pipelines can be found in (Samsonov and Blais-Stevens 2023).Deformation signals observed in the Yellow River Delta (Figure 7a) have been previously described in (Zhang et al. 2015(Zhang et al. , 2019;;Liu et al. 2021).The single most significant source of the noise observed in almost half of all the frames is an annual periodic signal similar to the one shown in Figure 9.It can be due to the seasonal variations in temperature, pressure and water vapor pressure in the atmosphere, the seasonal change in radar penetration depth, or both.However, because this noise is observed only at elevations below 2,500 m, it is believed to be mainly due to the seasonal change in radar penetration depth.The atmospheric noise can be reduced by applying various available atmospheric corrections (Jolivet et al. 2011;Bekaert et al. 2015;Wang et al. 2021), but it remains  unknown if the available atmospheric corrections perform well in remote areas, such as the Arctic.Comparing the backscatter intensity time series with the deformation time series, we observe some similarities, particularly during the apparent rapid subsidence period from late fall to late spring.Interferograms computed during this period retain high coherence suggesting that dry snow, if present, remains transparent to SAR.We attribute the apparent deformation signal to a change in radar penetration depth; however, it requires further investigation and field validation.During summer to early fall, penetration depth fluctuates with precipitation.In late fall, surface temperature decreases below zero, causing water to freeze and the radar penetration depth to increase.This process continues all winter and early spring until a rapid surface melt occurs.If the deformation map that spans the surface melting remained coherent and was accurately unwrapped, it would have captured a rapid uplift, completing the annual cycle.In areas with significant topographic relief, seasonal temperature decrease occurs earlier at higher elevations, which increases radar penetration depth relative to low elevation, thus producing the apparent subsiding signal at higher elevations.
The seasonal transition from ice/snow to water in spring is often incoherent.Therefore, rapid penetration depth change during that period is not resolved, resulting in the missing half of the annual penetration-depth-change cycle.The first order Tikhonov regularization, implemented in MSBAS, then interpolates the temporal gap, producing an erroneous interannual trend (Figure 9a) instead of the periodic annual signal (Figure 9c) and absence of the long-term deformation trend.Even though this process is well understood (Wegmuller 1990;Nolan et al. 2003;Nolan and Fatland 2003;Baghdadi et al. 2018), to the best of our knowledge, there are no readily available solutions to mitigating it.Two approaches seem promising.The first approach is using interannual deformation maps.For example, using summerto-next-summer deformation maps can help to resolve this signal.However, such deformation maps usually have low coherence.The second approach is to use SAR backscattering intensity.The previous studies showed that SAR backscattering intensity is sensitive to surface temperature, which is also apparent in Figure 9.For example, Mironov and Muzalevsky (2013) showed that the average temperature, in a range from À30 to 25 degrees Celsius in the active topsoil can be derived from the radar backscattering coefficient, which then potentially can be used for computing penetration depth.
Computing deformation time series and deformation rates using the DInSAR þ MSBAS methodology requires significant processing resources.The processing software must be parallelized for work on a High-Performance Computer, and the parallelization algorithm has been proposed.The system performed well during processing 220,000 Sentinel-1 images over North America and Eurasia.The typical processing time required to complete 1 job was about 24 hours.About one-third of this time was used for computing individual interferograms with GAMMA, and the remaining time was used for computing deformation time series and deformation rates with MSBAS.Occasionally the processing was interrupted by hardware failure.After restarting processing, the system could reuse already computed deformation maps, but the MSBAS processing had to restart from scratch.Several frames with a large fraction of coverage over the water failed to process during the coregistration step.Further development and optimization of the system for reducing processing time is desirable.For example, building a DInSAR network for MSBAS processing requires knowledge of the average interferogram coherence, thus, all the interferograms need to be computed first.It is suboptimal to compute many interferograms, measure their average coherences and discard most of them when their coherence is too low.The better approach would be to estimate expected coherence based on a priori known parameters, such as spatial and temporal baselines, wavelength, polarization, geographic location, and land cover.Work remains to be done for optimal stitching of adjacent deformation maps.Having a continuous deformation map would allow for resolving long-wavelength deformation signals due to tectonic motion and postglacial rebound.Several processing steps (outlined in pink in Figure 1b) have not been parallelized yet.During the execution of those processing steps, only a fraction of the allocated resources are utilized, leaving most of the resources idle and unavailable to others.Further work needs to be done to develop an efficient post-processing strategy, for example, selecting an optimal time matrix rank threshold value and type and strength of spatial filter.Several methodologies for extracting deformation signals from noisy data have already been tested, such as the hotspot analysis based on the Getis-Ord Gi Ã statistics (Samsonov and Blais-Stevens 2023) that performed very well when used for automated landslide detection.Computed deformation rates for North America and Eurasia are available to the research community and can be downloaded from the data repository (Samsonov and Feng 2022).These deformation rates were high pass filtered to remove noise due to the seasonal penetration depth change, so the detection of active deformation processes is simplified.The observed artifacts, especially their spatial distribution and seasonal pattern, also have a significant value as they provide new insight into radar interaction with the earth's surface.Many deformation processes observed in these deformation rate maps, including large landslides, have been previously unknown to the research community.

Funding
The work was partially supported by the Canada Centre for Remote Sensing, NRCan under the Status and Trend Mapping Program.The work was also partly supported by Guangdong Province Introduced Innovative R&D Team of Geological Processes and Natural Disasters around the South China Sea (2016ZT06N331), the Guangdong Provincial Natural Science Fund (2020A1515011078) and the National Key R&D Program of China (2021YFC3001000).

Figure 1 .
Figure1.(a) Schematics of processing system in terms of resource allocation.Job (red dashed outline) produces deformation time series from over 100 Sentinel-1 SLCs of specific path/frame.Job is subdivided into tasks running on individual task units (25 GB and 8 cores each, outlined in black) that compute deformation maps from 2 Sentinel-1 images and which cannot be distributed between different nodes (blue dotted outline).Five tasks run on single node simultaneously.MSBAS software uses all 10 nodes available for job.Cluster made of 100 nodes can run up to 10 jobs simultaneously.(b) Schematics of processing system in terms of processing steps.Steps 1 and 3 require only 1 task unit, thus leaving most of resource idle.Step 5 uses all 10 nodes as 1 task unit.

Figure 2 .
Figure 2. Schematics for estimating matrix rank assuming case with 3 SAR images (shown as black points) and 3 interferograms (shown as black lines).Various scenarious comprise overdetermined (b), well-posed (c-e) and underdetermined (f-h) systems of equations, respectively.MSBAS version 10 produces valid pixels in computed deformation time series and rates in all these scenarious (b-h).

Figure 3 .
Figure 3. Spatial coverage of Sentinel-1 ascending (a) and descending (b) SAR frames with over hundred SLC images per frame used in this study.

Figure 4 .
Figure 4. (a) Deformation rate map computed from 149 Sentinel-1 images from path 42 frame 391.Color scale is clipped in range [À0.1; 0.1] m/year.Insert in top-left corner shows location of this region in world, and insert in top-right corner shows deformation rate contour lines plotted over optical image from Google Earth in m/year: solid lines from 0.06 (red) to 0.12 (purple) and dashed lines from À0.05 (purple) to À0.02 (red).White arrows show direction of motion.(b) Network of interferograms used for inversion.(c) Deformation time series for 5 Â 5 pixel area marked with red star in top-right insert.

Figure 5 .
Figure 5. (a) Deformation rate map computed from 128 Sentinel-1 images from path 108 frame 221.Color scale is clipped in range [À0.1;0.1]m/year.Insert in top-right corner shows location of this region in world, and another insert shows deformation rate contour lines plotted over Sentinel-2 image in m/year: in range from À0.03 (red) to À0.21 (purple).White arrow shows direction of motion.(b) Network of interferograms used for inversion.(c) Deformation time series for 5 Â 5 pixel area marked with red star in insert.

Figure 7 .
Figure 7. (a) Deformation rate map computed from 206 Sentinel-1 images from path 69 frame 119.(b) Deformation rate map computed from 147 Sentinel-1 images from path 41 frame 140.Inserts show locations of these regions in world, network of interferograms used for inversion and deformation time series for 5 Â 5 pixel area pointed to with black line(s).

Figure 8 .
Figure 8.(a) Deformation rate map computed from 181 Sentinel-1 images from path 79 frame 464.(b) Deformation rate map computed from 135 Sentinel-1 images from path 35 frame 441.Inserts show locations of these regions in world, network of interferograms used for inversion and deformation time series for 5 Â 5 pixel area pointed to with black line(s).Several landslides in (a) are marked with L.

Figure 9 .
Figure 9. (a) Deformation rate map computed from 134 descending Sentinel-1 images from path 149 frame 398.Inserts show locations of this region in world, network of interferograms used for inversion and deformation and backscatter intensity time series for 5 Â 5 pixel area pointed to with black line.(b and c) Deformation and backscatter intensity time series are manually reprocessed for closely located to P areas without Tikhonov regularization.

Table 1 .
Sentinel-1 ascending (A) and descending (D) SAR frames with over 100 SLC images per frame used in this study.