COMPARISON AND EVALUATION OF DIFFERENT APPROACHES FOR EFFICIENT PROCESSING OF LONG GROUND-BASED SAR TIMES SERIES

Ground-based Synthetic Aperture Radar (GB-SAR) is a monitoring tool which, once installed, acquires a large amount of data autonomously. For the IBIS-FM system, approximately 760 SAR images per day are acquired, which corresponds to more than 23 000 scenes per month. Therefore, this paper analysis different strategies for the interferometric processing of such large data stacks to find a compromise between accuracy, computational effort and the ability to (re-)process specific time intervals independently. This study compares the single master approach with the sequential approach and in addition two block-wise approaches. Moreover, a new baseline configuration called Daily Baseline Subset (DBAS) is compared which uses interferograms having a multiple of one day as temporal baseline. We evaluate them on a data stack of 30 000 images, acquired at Enguri Dam in Georgia. We check the unwrapping errors and the quality of the displacement estimation to compare the different configurations. We found that block-wise approaches show the best results considering unwrapping errors and Root Mean Square Error, while in our study the DBAS approach shows to the most plausible displacement map which is also dependent on the individual reduction of atmospheric noise.


INTRODUCTION
Interferometric Synthetic Aperture Radar (InSAR) methods have been developed originally for observing motions of or on the terrain fron satellite SAR sensors. Bamler and Hartl (1998) explains in detail the measurement concept, the processing chain and limitations of spaceborne InSAR. This technique has then been incorporated into ground-based SAR (GB-SAR) systems in Tarchi et al. (2000) to overcome some of these limitations and allow for flexible acquisitions with higher repetition rate. Since then, it has been applied in several case studies and application scenarios, such as DEM formation,  landslide monitoring Tarchi et al. (2003), Herrera et al. (2009), mine monitoring Cao et al. (2021), glaciers monitoring Luzi et al. (2007). It can also be used for dam monitoring as demonstrated for example in Alba et al. (2008), DiPasquale et al. (2018) or Wang et al. (2020). An overview of the applications of GB-SAR for deformation measurements is available in Monserrat et al. (2014).
(GB-SAR) systems have several advantages over classical geodetic methods. Once installed, it takes SAR acquisitions at a regular sampling rate autonomously. Typically, the sampling interval ranges between 10 seconds and a few minutes. In our study we use the IBIS-FM from IDS Georadar with a sampling period of approximately 2 minutes. The acquired SAR image corresponds to a spatial sampling of the observed grid with a coherent imaging system. Each pixel provides information on amplitude and phase of the backscattered signal. Interferometric processing of SAR images results in phase differences -the so-called interferometric phase -which is proportional to scatterer motions (deformations) in the radar beam's line-of-sight up the phase's 2π ambiguity. Multiple differences can be set-up, when expanding interferometry along the temporal and/or spa- * Corresponding author: matthieu.rebmeister@kit.edu tial dimension. According to Monserrat et al. (2014), the precision of the GB-SAR depends on the target characteristics and its distance to the radar head and ranges from sub-millimetres to a few millimetres.
Adapted from interferometric processing of spaceborne SAR imagery, the processing of GB-SAR time series faces the same main processing steps: (i) Interferogram formation of two acquisitions, to obtain phase differences with respect to a given temporal and spatial reference. Arbitrary acquisition combinations are possible, however, every image needs to be included in at least one interferogram. (ii) Spatial or 3D phase unwrapping to obtain real-valued phase values from complex-valued wrapped interferometric observations. (iii) Temporal integration to transform double differential observations to single differential observations with respect to a given spatial reference. There are various processing chains for GB-SAR available in the literature such as Monserrat (2011), Rodelsperger (2011), Hu et al. (2019 or Wang et al. (2019). In the work of Monserrat (2011), N − 1 sequential interferograms are calculated from N GB-SAR images, from which phases are integrated after spatial phase unwrapping. After that, the results of the phase unwrapping are checked with a temporal phase check based on a temporal redundant network. However, it can be time consuming for large time series. In the work of Rodelsperger (2011), a real time processing framework is proposed, using a Kalman filter for temporal unwrapping. The number of interferograms processed is also N − 1. The other two approaches are based on block or unit processing, where interferogram formation and unwrapping are done in each block independently. In Hu et al. (2019) the single master approach is used while in Wang (2018), the Small Baseline Subsets (SBAS) from Berardino et al. (2002) are preferred.
The processing of long-term GB-SAR time-series face many challenges, such as strong variations of the atmospheric con-ditions or the decrease of coherence over time. Moreover, the amount of data needs to find a compromise because it would be unfeasible and useless to compute all possible interferograms, though redundancy is helpful to increase estimation accuracy and detect errors during phase unwrapping.
In this study we test and compare several processing strategies with different baselines configurations. We evaluate the approaches by applying them to a data stack acquired at the Enguri Dam in the southern Caucasus Mountains in Georgia. In this case GB-SAR is used to study the long-term deformation of the dam with geodetic accuracy, unlike in most applications where the GB-SAR is used as an early-warning system to monitoring areas for strong and spontaneous surface displacements.

METHODOLOGY
In the remainder of this section n denotes an index in the set [|1; N |] and tn the time of the n th acquisition. For each point in interferogram m, the signal is considered as the sum of the displacement (ϕ disp ), atmospheric (ϕatmo) and noise (ϕ noise ) components as follows: To reduce the amount of data to process, a first preprocessing step containing a coarse pixel selection is made on the first 30 scenes of the dataset. Pixels which have an average coherence above 0.85 and an amplitude dispersion below 0.35 are considered for the rest of the processing. They are considered as Persistent Scatterers Candidates (PSC). The coherence is computed on sequential interferograms with a 5x5 window. The amplitude dispersion DA is defined as in Ferretti et al. (2001), as the point-wise quotient of the mean amplitude and the standard deviation of the amplitude.
The preprocessing also contains the geocoding which is required to determine the elevation dependent component of the atmospheric phase as described in Iglesias et al. (2014).
The minimal baseline configuration are the Sequential (Seq.) and Single Master (SM) stacks. For N acquired scenes, the number of interferograms in both cases is N − 1. For the SM approach, one scene at time tM is considered as the master. The middle scene is the adequate choice because it provides the shortest temporal baselines. The main drawback of this technique is that the coherence decreases with time as it can be seen in Figure 3. Moreover, the cumulative displacement and atmospheric effects introduce a strong signal which is difficult to unwrap for large temporal baselines. This approach is commonly used in satellite InSAR such as in Hooper (2006).
The sequential approach enables a real-time monitoring and does not suffer from the decorrelation problem. However, our experience is that with a strict sequential integration small errors due to processing propagate through the whole time series and sum up to large total errors. While this may not be a drawback for real time event detection in early warning systems, it is a severe shortcoming for applications in long-term geodetic monitoring. A block-wise processing synergises the advantages of short baselines but having a link to the master acquisition. The Block-Wise Single Master (BW-SM) approach is descried in detail in Hu et al. (2019). We compute single master interferograms for each day with an overlap of at least 10% the size of the block. To make the process a little more robust, we test another approach, which uses the same block-wise processing, but with the Small Baseline Subsets (SBAS) concept. For each block, all interferograms with a temporal baseline smaller than a given time period are used. This approach is the Block-Wise SBAS approach (BW-SB) and is fully described in Wang et al. (2019). The main disadvantage of this method is the long processing time.
The classical Small Baseline approach (SBAS) can not be applied for large GB-SAR stacks because of the amount of interferograms it will generate. Indeed, with an acquisition every 2 minutes or even less, it is possible to create a huge amount of small baseline interferograms, depending on what is considered as small temporal baseline. We adapt the SBAS approach such that we make use of the redundancy by limiting the total number of interferograms to a reasonable amount at the same time.
Taking into account the large contribution by daily variations of the atmospheric phase delay, a promising approach is a reduced SBAS interferogram stack, where all interferograms have one day or integer multiples of one day as temporal baseline. This approach is expected to reduce the diurnal atmospheric signal in each interferogram significantly. Moreover, those interferograms also show higher values for (ϕ disp ) with lower SNR, than interferograms with short temporal baelines. We refer to this baseline configuration as Daily Baseline Subset (DBAS) approach.
For the remainder of this paper, we will denote with SM the set of scenes used as master for the interferogram formation and SV the set of scenes used as slaves for the interferogram formation. Both sets have the same length. Those sets depend on the baseline configuration chosen. For example, in the SM cases, SM = {1, 1, ...1} and SV = {2, 3, ..., N }.

Figure 1 recaps the different baseline approaches. For the BW-
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France SM case, the represented block-size is 4 scenes and the overlap in this case is 1 scene. Overlap and block-size can be set by the user. We alternated full and dashed lines in this illustration to indicate consecutive blocks. In the case of the BW-SB, we only represented one block to reduce the amount of interferograms.
In Figure 1, a block size of 5 scenes and a maximal baseline of 3 scenes is represented as example. Full lines are for interferograms made with consecutive scenes. Dashed and dotted lines are used to represent interferograms with a baseline of 2 and 3 scenes respectively. For the last configuration, DBAS, we also did not plot all the interferograms to enhance the visibility. In this case, full lines represents interferograms with a baseline of one day and dashed lines represent interferograms with a baseline of two days. To compare only the baselines configuration, we used similar processing chains with identical parameters, where possible.
The first step is the selection of the stable pixels in the whole stack with the amplitude dispersion threshold as fast estimator of the phase noise. The points satisfying this criterion are called Persistent Scatterers (PS). For the Seq., SM, and DBAS approach, this criterion is applied over the whole stack. For the BW-SM and BW-SBAS, this criterion is applied on each block, meaning a different number of points per block. The advantage is that we can process points that are not always stable.
Then the interferograms are constructed depending on the baseline configuration and the phase of the selected reference area (considered as stable) is subtracted to all interferograms.
All interferograms are filtered with a 2D Goldstein filter Goldstein and Werner (1998) to make the unwrapping with the software SNAPHU Chen and Zebker (2001) easier. In configurations with small baselines (i.e. Seq. and BW-SBAS), some processing time could be saved by only unwrapping interferograms with potential phase jumps. If the time between two acquisitions is small, the atmospheric conditions are comparable and almost no deformation occurs, the whole interferogram will be close to zero and would not require phase unwrapping. To identify the interferograms requiring a 2D phase unwrapping, we use statistical quantities described in Mardia et al. (2000): the mean resultant lengthR from Equation 2 and the circular variance V from Equation 3.
where N PS is the number of PS in the interferogram θj denotes the wrapped phase of point j If the circular variance of the interferogram is below 0.1, the phase of the interferogram is accepted without further unwrapping. The threshold was determined based on simulations with synthetic data. It was chosen conservatively, to be sure to unwrap all necessary interferograms.
The atmospheric phase delay is corrected using a spatial model on the unwrapped interferograms. The simplest model is a linear range dependent phase delay. Noferini et al. (2005) enhanced this model by fitting a polynomial signal in range. This model is also used in Zhang et al. (2011). Iglesias et al. (2014) derived that the atmospheric signal can be described, in first approximation of the Taylor expansion of the atmospheric phase expression, by a range-height dependant model as follows: where a0, a1, a2 are the coefficients to estimate r(P ) = range of point P h(P ) = height of point P This model is estimated in a first step with all available points. In a second step, the residuals of the estimation are analysed and the adjustment is made one more time, only for points having residuals smaller than a given threshold. More details are given in Iglesias et al. (2014).
Once we corrected the atmospheric phase component ϕatmo, it is possible to estimate the displacement rate of each point. The classical choice is to consider the temporal baseline of each interferogram. But other choices could be more practical in some cases of displacement, in particular when it is nonlinear in time. In our case, we monitor the elastic response of the dam to changes in the water level and therefore we use the water level as baseline in order to checking the deformation estimates for plausibility and expressing themotion rates in line of sight as displacement per meter water level. Let Φ being the vector of unwrapped interferometric phases. Then Φ = (ψ(ti) − ψ(tj))i∈S M ,j∈S V , where ψ denotes the integrated phase and ti the acquisition time of image i.
With these notations, we create the design matrix A as a single column containing the coefficients − 4π λ (B(ti) − B(tj)), where λ is the wavelength of the radar system and B the baseline of choice. A standard least squares approach is used to estimate the displacement rate d. The term rate here depends on the type of baseline. In our cases, we use the water table height of the reservoir as baseline and estimate the displacement rate as displacement per meter water level change in mm/m, the denominator corresponding to the water level change. This inversion enables to compute for each point its Root Mean Square Error (RMSE) defined as: The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France

vector of residuals M = number of interferograms computed
Another parameter evaluated is the mean absolute deviation, denoted with MAD and given by: This parameter is also evaluated as the RMSE can strongly increase in case of unwrapping failures, leading to residuals in the neighborhood of 2π.
For the BW-SM and the BW-SBAS approach, this processing chain is applied for each block and we then have a displacement rate for each block. To combine all those results, we find the points which are processed in all blocks and we compute a weighted meand for each point. The weights used are the inverse mean RMSE of each block.
whered(P ) = weighted mean of the displ. rate of point P RMSEi = mean RMSE of block i di(P ) = displacement rate in block i for point P N blocks = Number of blocks

EXPERIMENTAL DATASET
Our experimental dataset consist of 30.000 acquisitions at the Enguri Dam. The Enguri dam is 270 m high and is characterized by its strong periodic water level variations throughout a year which can be up to 100 m. The observation period is from April 20 to May 30 2021. The time between two GB-SAR acquisitions is 113 s. The preprocessing includes the image focusing and the preselection of 4938 stable persistent scatterer candidates (PSC). We used all available acquisitions including those with considerable noise, to also compare the robustness of the methods.
The observation standpoint is shown in Figure 2. A shelter was constructed later to protect the system from wind, rain and snow.
For the BW-SM approach, we used a block size of one day, i.e. 760 acquisitions per block and an overlap of 10% of the block length, which corresponds to 76 scenes. The same block size was used for the BW-SBAS approach. A constraint of 4 interferograms per scene was added, summing up to 3030 interferograms per block for the BW-SBAS approach.
For the DBAS approach, one would have to process more than 560.000 interferograms when considering all possible interferogram combinations with a baseline of one day and multiples of one day. To keep processing time at a reasonable level, we restricted the temporal baselines of interferograms to a maximum of 4 days.

RESULTS
Considering the first scene as master, we computed firstly the mean coherence of all the PSC and present the decay over time on Figure 3. For this time interval we can estimate an almost linear loss of the mean coherence of around 0.05 per day. This is a first indicator showing the limitation of the single master approach for long time processing.
Using very large baselines in the single master approach lead to severe unwrapping errors. Figure 4 is an interferogram between two acquisitions separated by 15 days. The sides of the dam are easy to unwrap but it is not the case of its central part. This will lead to a series of some interferograms appropriately unwrapped and some interferograms with unwrapping failures.
The results of the main evaluation metrics are presented in  2. We used eight parameters, five regarding the processing and three to evaluate the quality of the approaches. The first line (RT/BU/PP) indicates the real-time capability of the respective approach, i.e. if the baseline configuration enables a real-time monitoring (RT) with instant and independent update after each acquisition, or block update (BU) for block-wise approaches which need to collect a certain number of acquisitions in a block, or post processing (PP) for the approaches SM and DBAS, which include temporal baselines over the whole observation time. Even if it would be possible to have a realtime monitoring with the single master approach, we categorized 'PP' because of the decreasing coherence using a fixed master. In long-term applications the master acquisition would be chosen to be somewhere in the middle of the time series, and a block-wise processing scheme with long blocks would be applied. The same category 'PP' applies for the new DBAS approach, because it is currently programmed for a fixed data stack.
The second line in table 2 lists N ifg , the total number of interferograms used with each approach. Tproc. corresponds to the processing time in minutes. The absolute values should not be overinterpreted, as the processing approaches have not yet been tuned with respect to computation time. However the comparison between the methods gives a good first idea of which approaches are more CPU-intensive. In the whole processing chain, the unwrapping is the step requiring the most time, with nearly 90% of the total processing time.
The value Nps refers to the number of processed points satisfying the amplitude dispersion threshold. For block wise approaches, we computed the number of points present in all blocks (intersection). Not surprisingly, this number is lower than for other approaches. That is why we also show the mean number of points in a block with N ps,block , which in turn is higher than the number for the other approaches.
The first of the quality measures in table 2 is the percentage of unwrapping failures in the whole stack (Uwf. sp.). In Wang et al. (2020), this quantity is calculated by computing the phase misclosures in closed temporal loops. This quantity can only be computed in the case of a redundant network. In our case, it is not possible in the Sequential, SM and BW-SM case and we therefore did not compute the percentage of unwrapping failures in this way. Instead, we computed an alpha shape of our scene with a radius of 15 m. The alpha shape has the advantage over the Delaunay triangulation that it avoids long edges Edelsbrunner et al. (1983). We then calculated the phase differences of each edge in the triangulation for each unwrapped interferogram. If an absolute difference is larger than a given threshold (set to π in our experiment, it is considered as an unwrapping error. To compare all the methods, we compute the percentage of wrong edges among the number of all edges. The second and third parameters used are the Root Mean Square Error and the Mean Absolute Deviation respectively defined in Equation 5 and 6. In this comparison, the sequential approach has by far the poorest performance, as it has the largest errors among all the methods. This is related to the large number of unwrapping errors and can be related to the fact that the processing chain implies a 2D filtering and unwrapping and then the integration. This procedure does not reduce noise over time but is prone to create points with large phase jumps which are not really expected. We can also see this in Figure 6 (A) and (B), where (A) represents the spatially (2D) filtered, wrapped interferograms before the phase unwrapping. The mean value is close to 0 but there are periods affected by strong noise and/or atmospheric signals. Panel (B) shows the phase after phase unwrapping, phase integration and atmospheric correction. We can clearly see that around noisy areas, phase jumps occur and lead to incorrect trend estimation. According to the RMSE, the SM and the DBAS give similar results. However, regarding the respective deformation in Figure 5, it seems that the SM approach is still affected by an atmospheric pattern on the right side of the dam which is not the case of the DBAS approach. The final maps of estimated displacement rates are illustrated in 5. As can be seen, the displacement rates resulting from by the sequential approach are not satisfying. They behave noisy and the values are completely different from that of other techniques. The results of all other approaches show more or less the same trend with a negative deformation on the left of the dam and a positive increasing deformation on the right. However, the amplitudes of the deformation is not exactly the same which means that the temporal baseline configuration has an influence on the results. The panel on the bottom left shows the location of the points chosen for the time series display in figure 6. Figure 6 displays the time series, for the points P1, ..., P4. Panel (A) represents the 2D filtered times series. We separated each acquisition with a constant of 2π for better visualisation. Panel Figure 5. Deformation map in Cartesian radar coordinates from the five different approaches and selected points for the time series plot (see Fig. 6).
(B) shows the times series after unwrapping and atmospheric correction. The wrapped time series of the SM approach are shown in (C) and finally, (D) shows the unwrapped and atmospheric corrected time series of the SM approach. We can clearly see the unwrapping errors, especially at begin and end of the times series. The black lines represent the estimated trend with the least squares approach. For points P1 and P3, we clearly see that the curve fits better to the points than what would be achieved by using temporal baselines and estimating a straight line representing a linear deformation in mm/year.

CONCLUSION AND OUTLOOK
In this paper, we choose several approaches to process long time series of Ground-Based SAR data. We choose a classical processing chain with tools known for their efficiency and reliability. After applying each baseline configuration to a data stack of 30000 GB-SAR images, acquired at the Enguri dam, we showed that the choice of the interferogram configuration has a strong impact on the obtained results. With simple tools quantifying the inner quality of the estimation, we can conclude that block wise approaches seems to give better results. However, to quantify the accuracy of the methods, it would be required to compare the obtained deformation maps with an existing known method, such as a classical geodetic deformation monitoring. Figure 6. Times series wrapped and unwrapped form the Seq. and SM. approaches The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France