InSAR Performance for Large-Scale Deformation Measurement

This article deals with the analysis of InSAR performance for large-scale deformation measurement. The study evaluates the use of models, especially numerical weather prediction reanalysis, to mitigate disturbances in SAR interferograms. The impact of such corrections is evaluated by analyzing short-time baseline phase variograms in order to derive a lower bound for the interferometric accuracy, especially at large distances. The variance is then propagated from single interferograms to deformation rates. Finally, using GNSS measurements, the predicted error bars are validated on a large Sentinel-1 data set.


I. INTRODUCTION
T HE Sentinel-1 mission [1] systematically provides SAR data suitable for interferometric applications with a swath width of 250 km. Future SAR satellites will further extend this to even larger swaths [2], [3]. The measurement of tectonic movements particularly benefits from the large-scale deformation measurements that SAR interferometry provides. Nevertheless, since the magnitude of the relative error increases with distance [4]- [6], the performance of relative deformation measurements between very distant points may not achieve the required accuracy for tectonic applications [7]. Typically, interferometry works well for applications such as infrastructure monitoring or urban subsidence since they involve relatively short scales (10-20 km). It is then common practice to remove the low-pass spatial frequencies of the estimated deformation signal by filtering or detrending. This is appropriate if the application is not the measurement of large-scale deformation. However, the rising scientific interest in using InSAR to measure this kind of deformation [8] requires preservation of all the spatial wavelengths of the deformation and proper treatment and knowledge of the error bars at such distances.
The atmosphere has always been recognized as one of the main error sources in interferometric measurements [4]. The main problem stems from its spatial characteristic that prevents complete separation from the deformation signal by spatial filtering. Temporal filtering, even if it improves the time series quality, is also not a general solution. This Manuscript received July 20, 2020; revised October 8, 2020; accepted November 14, 2020. Date of publication December 3, 2020; date of current version September 27, 2021. This work was supported in part by Helmholtz-Gemeinschaft under the Project Tectonics and Volcanoes in South America with InSAR (TecVolSA). (Corresponding author: Alessandro Parizzi.) The authors are with Remote Sensing Technology Institute, German Aerospace Center (DLR), 82234 Wessling, Germany (e-mail: alessandro.parizzi@dlr.de).
Digital Object Identifier 10.1109/TGRS.2020.3039006 is due to nonlinear deformation, uncompensated seasonal effects, and irregular sampling due to snow cover. The performance of the estimated deformation rates is still limited by the atmospheric delay. Therefore, external information is necessary to reconstruct and compensate for the unwanted delay in the interferograms [9]. Information from GNSS stations or other remote sensing instruments, such as MODIS or MERIS, are typically used despite the different resolutions [10]- [12]. The correction of tropospheric stratification was successful using both model-based and data-driven approaches [13]- [15]. In the past decade, the use of numerical weather prediction (NWP) became a valuable tool for systematic correction of both SAR group and phase delays, hence improving the performance and capabilities of the techniques [16]- [20]. Starting from this basis, this article first investigates interferometric performance at large distances when ECMWF ERA-5 NWP data were used for the corrections. Then, based on a simple theoretical model, a possible error description is provided and its validity is demonstrated on real data. Before going further, a brief discussion of other error sources affecting the measurements at long distances, and how they were considered, is warranted.
The ionosphere can generate huge errors at low frequency (L-band and P-band) SAR [21]. In C-band, its effect is less severe even though exceptions have been reported [22]. In such cases, compensation using the split spectrum technique can reduce the error to the level of the troposphere so that it is no longer the limiting factor [23]. Here, ionospheric effects have been corrected using CODE models in order to mitigate larger scale ionospheric effects [24]. The effect of solid earth tides (SETs) involves scales much larger than the SAR swath; however, their projection onto the SAR line of sight (LoS) along the swath can generate centimeter-level phase ramps. Therefore, SETs were compensated using the IERS 2010 convention [24], [25]. Finally, SAR missions are characterized by very stable oscillators [26] and very good orbit knowledge [27], [28]. This study does not consider the uncertainties related to the residuals of ionosphere and SET corrections nor the accuracy of the state vectors. Instead, it is implicitly assumed that the residual error of signals related to the ionosphere, SETs, and orbit state vectors is negligible compared with the tropospheric signal that becomes the limiting factor for interferometric performance at large distance, as verified in the Appendix. Violation of this assumption leads to an underestimation of the tropospheric component. Hence, this study is intended to outline the bound for the improvement related to tropospheric corrections. The use of NWP for mitigating the tropospheric signal in the interferometric phase is limited both in accuracy, due to the limitations of the model itself, and in resolution. The best sampling available is 9 km in the case of ECMWF products, several orders of magnitude larger than the SAR resolution, even after interferometric multilooking. The accuracy of NWP can be strongly site dependent since it depends on data availability and the fit of the model to actual tropospheric conditions. For the ECMWF products considered here, there is a strong dependence of accuracy on latitude [29]. However, it should be pointed out that since the interferograms are intrinsically relative, the impact of model bias in the simulated interferometric phase is limited. More important is the variance of the error between the real zenith delay z t and the modeled delay z t = z t − z t . This can be seen as a Gaussian process in time and space with mean μ , variance σ 2 , and spatial correlation function . The interferometric error between two points A and B for reference acquisition r and secondary acquisition s is where θ is the local incidence angle. The effect of the altered path delay due to topography is considered using the exponential mapping function as in [30] with h A and h B the height of the points with respect to the ellipsoid and H the so-called height of the atmosphere, typically fixed to 7 km. Assuming that the error is first-order stationary, the mean interferometric phase error is E[φ ] = 0 even if the residual error has nonzero mean, μ = 0 as visible from (1). Assuming that the atmosphere and the error are temporally uncorrelated, the interferometric phase error power is where is a correlation function that (A ≡ B) = σ 2 and (A, B) → 0 at large distance. The residual error in the interferometric phase between points A and B described in (2) is a parameterization of the sample variograms that can be calculated from data. Equation (2) approaches 0 as the distance d → 0 and is ∝ σ 2 as d → ∞. This means that at large distances (≥40-50 km), the interferometric measurements become limited by the accuracy of the models used for the correction, as verified in the Appendix.
In order to study the troposphere-induced error and its mitigation through the use of NWP, short-time interferograms are generated from every acquisition with the next. This approach should prevent bias due to a common reference as well as eventual seasonality present in the tropospheric delay 1 and the impact of deformation. The latter point requires clarification. The fast revisit time of the Sentinel-1 mission allows temporal baselines of 6, 12, or, in worse case, 24 days. Since the residual error, including effects such as tropospheric turbulence, is assessed to be on the order of centimeters, deformation rates of several tens of cm/y would be necessary in order to produce effects comparable to those from the troposphere. Such rates occur due to landslides or in mining areas that are typically restricted in time and space. It follows that variograms computed by averaging over the scene should not be strongly impacted by such rapid deformation types. The projection of tectonic plate motion onto the range varying LoS is a large-scale effect that could seriously bias variogram estimation. Such movements are mainly horizontal and can reach 6 and 7 cm/y. Projection onto the LoS (30 • to 45 • at near and far range, respectively, for Sentinel-1), combined with the previously mentioned revisit times, shows its impact to be negligible. The eventual presence of seismic events in the time series should be assessed and coseismic interferometric pairs avoided when estimating the variograms.
The tropospheric error is considered temporally uncorrelated. This hypothesis should hold at least for the troposphere corrected phase since all seasonal effects are included in the NWP and are hence compensated. The interferometric phase variogram between acquisitions n and m, n,m , is then the sum of the two variograms of the errors in acquisition n and m, n and m where φ is the phase at points A and B and d is the distance between the points. By averaging the set of short-time baseline variograms = E[], it is possible to derive the average behavior of the residual interferometric phase error. Performing this experiment before and after the NWP correction allows an analysis of the impact of the corrections at different scales. The interferometric data and NWP models have different resolutions, and hence, different behavior is expected at different scales. Fig. 1 shows an example of variogram behavior before and after correction using a long stripe of mosaicked Sentinel-1 interferograms (see also the Appendix and [24]). Although the spread of the gain is large at all scales, the temporal average of the variograms clearly shows the considerable gain due to the NWP correction that increases at large scales (≥40-50 km). This is particularly important in tectonic strain applications where high relative accuracy is required at large scales.
Assuming that the troposphere is the main source of error at large distances (≥40-50 km), one can now derive its impact on deformation rate measurements. The variogram of the linear deformation rate estimates is obtained by scaling by the linear regression formula [31] v (d) = 1 2 where the factor 1/2 accounts for the common reference acquisition, M is the number of acquisitions, and σ 2 t is the spread of the temporal sampling    2 shows an example of the impact of NWP corrections on the deformation rate estimates at different scales when varying the length of the time series. It is interesting to note how the typical 1-mm/y accuracy requirement is achieved after four years with the corrected data set, while almost eight years are required if no correction is applied. In this section, the rationale of the study has been briefly described accompanied by examples. In Section III, this approach is extended to a larger set of data worldwide to analyze and validate what has just been presented.

III. PERFORMANCE OF THE NWP TROPOSPHERIC CORRECTION
Section II demonstrated how compensation of the tropospheric delay reduced the spatially correlated noise. Variograms are a convenient way to characterize this noise component. By comparing the interferometric phase variograms before and after the correction, it should be possible to observe the gain at different scales. Since the data used for compensation are not provided at the same resolution of the interferograms, a significant reduction of the tropospheric contribution can only be expected at scales larger than the NWP model resolution. Here, ECMWF ERA-5 data are used with a spatial gridding of 30 km. It should be noted that the accuracy of the correction has been demonstrated to vary considerably worldwide [15], [29].
In order to comprehensively evaluate correction performance, 146 Sentinel-1 stacks 2 have been processed over various regions of the world. Each of these stacks spans at least four years and contains up to 200 acquisitions. As described in  Section II, only short-time baseline interferograms were used to ensure that the troposphere signal dominates. The average variogram was also calculated at different scales. A measure of the gain provided by the correction is provided by the metric the ratio, in decibels, of the average variograms with and without the correction. Fig. 3 shows G dB (d) for stacks from the German Ground Motion Service separately from all others.
An evaluation of the performance at large distances was also performed. The average variograms before and after corrections were fitted using typical variogram models (exponential, Gaussian, and spherical), hence reducing the variance of the variogram estimate. Using this information together with temporal sampling information and (4), the performance in measuring the deformation rates at 150 km was estimated with and without the tropospheric corrections. Fig. 4 shows the histograms of the accuracies in mm/y with the German Ground Motion Service shown separately from all others.
Mitigation of tropospheric effects and the correct characterization of their spatial characteristics are critical in developing error models for missions that allow [1] or will allow [2], [3] systematic interferometric processing on a global scale.
Figs. 3-5 show that a significant gain is attained when the tropospheric delay is corrected. As already pointed out in [29], the gain is region dependent. A very high gain is achieved where the model accuracy is also very good, such as in Europe or USA. Only one of the 146 stacks, located in Indonesia, was the corrections found to slightly worsen the performance at medium scale and is clearly visible in Fig. 4.
In Fig. 3, a smaller but nonnegligible gain is also observable at short scales (≤10 km). This value can vary from 0 up to  Fig. 6. The tropospheric models can correct stratification effects at scales smaller than the models' resolution since the topography, available at the reference DEM resolution, can be used to project the models to phase errors [15]. Therefore, such a gain can only be attributed to the compensation of tropospheric stratification since the resolution of the correction model is much coarser than the topographic variation.

IV. ACCURACY OF INSAR VELOCITIES AND VALIDATION
Knowledge of displacement rate precision is of central importance in geophysical studies of large-scale tectonic movements where the precision requirements are often very high [7]. In Section II, it has been shown how the residual tropospheric noise can be characterized from short-time baseline interferograms and how this can be translated into displacement rate accuracy. This approach involves assumptions about the residual error, i.e., the short-time baseline interferograms characterize the error, that could be questioned. For the sake of completeness, the behavior of the displacement rate as described in (4) must be verified in practice.
This was achieved by cross validation with GNSS measurements using the rationale developed in [32] for integrating GNSS and InSAR measurements. Consider N locations where two displacement rates v D (in LoS) and v G (in x, y, z) are obtained from InSAR and GNSS, respectively. With s as the LoS vector, the difference between the two velocities at position i , i , is where v ref is the velocity of the reference point used in InSAR processing, δ i represents spatially varying residual atmospheric zero-mean error from InSAR processing, and n i is the random error. The quantity i is obtained in practice by selecting a set of InSAR measurements located within a radius of a few hundred meters around the GNSS stations. In order to reduce the influence of clutter, the interferometric measurements are averaged. In this way, the random noise component n i can be considered almost completely due to GNSS. To use ( (6)), a statistical description of the vector is required, i.e., the covariance matrix R of the difference measurements. Three terms can be distinguished: 1) random noise from the GNSS measurements σ 2 G,i obtained by projecting the variance of the (x, y, z) components of the velocities onto the LoS; 2) random noise from the InSAR measurements σ 2 D,i ; and 3) spatially correlated noise due to residual atmospheric effects. The full covariance matrix is The first two terms are diagonal matrices since they represent the spatially uncorrelated errors of the independent GNSS and InSAR measurements, respectively. The last term, C δ , representing residual atmospheric error, still has to be found [33]. It can be determined from the average variogram corr (d) according to (4).
Given a full statistical characterization of the GNSS/InSAR differences and assuming the GNSS accuracies to be correct, we can now verify whether v actually represents the variogram of the velocity error. First, from the vector of N differences between GNSS and InSAR velocities, we compute all (N 2 − N)/2 unique pairwise differences between the elements. Their mean is expected to be zero according to (6) and their variance is The scales at which v can be verified by this method clearly depend on the spatial distribution of GNSS stations, with the range of scales being set by the smallest and largest distances between the stations.
A statistic defined as the standardized difference should follow a standard normal distribution if the error description is correct. To test the validity of the variance model in (8), a χ 2 -test can be performed on T . The confidence region at significance level α for the true σ 2 T is where the number of degrees of freedom is N − 1 andσ 2 T is the sample variance of T . The validation approach will focus on analyzing the distribution of T and the confidence interval for σ T .

A. Cross Validation of the German Ground Deformation Service Data Set Using GNSS Stations
The cross-validation scheme was performed on the German deformation map. The data set includes 41 stacks acquired in both ascending and descending geometries. This region was chosen due to easy and open access to a very dense GNSS network from Nevada Geodetic Laboratories [34], [35] with all stacks containing sufficient GNSS stations.
The confidence intervals for σ T at 5% significance level are shown in Fig. 7(a). Of course, the dependence on the number of GNSS stations, shown in Fig. 7(b), is strong. However, for almost all stacks, the confidence intervals include σ T = 1 indicating acceptance of the null hypothesis H 0 : σ T = 1 versus the alternative H 1 : σ T = 1. Fig. 8 shows the distribution of T pooled over all stacks, along with the nominal standard normal distribution under the null. The pooled value ofσ T is 1.03, showing very good agreement between the model and GNSS data. Fig. 9(a) shows the error variogram from InSAR data compared to all possible pairwise differences ( i − j ) 2 in order to show the fit of the variogram model at various scales. Since the accuracy of GNSS measurements also plays a role, the GNSS error contribution is included as per (8). The power of the GNSS contribution can be close to that of the InSAR one, especially at smaller scales where GNSS may distort the validation. For this purpose, a special validation was designed.

B. Validation at Local Scales for the German Deformation Map
Given that there is only limited large-scale deformation over German territory, the experiment was repeated using a Fig. 9. Comparison between variograms (blue line) and all pairwise differences ( i − j ) 2 (red dots) for an ascending stack of the German Ground Motion Service. In order to highlight the plausible area for the measurements, the blue dashed lines correspond to (3σ ) 2 variograms and the bold one to σ 2 . The green dots correspond to the pairwise differences after averaging over 2-km bins, as in (8). The variograms are shown (a), including the GNSS error component as in Section IV-A and (b) using the zero-velocity synthetic network as in Section IV-B. network of synthetic GNSS having no motion generated on a 0.1 • × 0.1 • grid. The well-known large-scale deformation due to mining, gas extraction, and so forth were excluded. Since continental drift can generate ramps over range, the motion of the Eurasian plate was also removed from the InSAR velocities [36]. It is noted that since the motion model is constant over the whole plate, residual effects may be present. However, since the analysis is on relatively small scales (≤45 km), it is assumed that such effects are also small.
Confidence intervals for σ T are shown in Fig. 10 for each stack. The intervals are much smaller than in Fig. 7 since the amount of available data is much larger. In general, the confidence intervals lie close to the null value σ T = 1, but the agreement is not as great as previously. This may be related to the hypothesis of zero motion failing or other residual effects. Fig. 11 shows the histogram of T pooled over all stacks, the pooled value ofσ T is 1.15. The variogram of the pairwise differences is also shown for a single stack in Fig. 9, in Fig. 9(a) using real and in Fig. 9(b) synthetic GNSS data.

V. CONCLUSION
This study discussed the potential and the limitations of InSAR for large-scale displacement measurements. The main concluding points are as follows.
1) The measurement of large-scale displacement signals precludes spatial filtering since all displacement components must be preserved in the final product. In order to provide sufficient accuracy at large scales, tropospheric corrections are recommended. 2) Performance may be site dependent, but the target of 1 mm/y at 100 km is achievable in five years given good weather models (in Europe or USA) or in general where the errors are small due to reduced atmospheric delay, i.e., over the Tibetan Plateau. 3) Short-time baseline variograms allow a reasonable evaluation of displacement rate performance providing valuable error bars to users. Such measurements must be included in the product since residual tropospheric noise is the main limiting factor at large distances. 4) In the Appendix, a lower bound on the interferometric phase accuracy at large distance was outlined and related to the accuracy of the models used for the corrections. Finally, as NWP models are continuously being improved, follow-on effects for InSAR can be expected.
Future work shall address the correction of tropospheric stratification in more detail. In this study, the problem has been only partially discussed since the focus was more on large scales. The designed framework does not allow for a quantitative validation as the one carried out for Germany since a dense network of GNSS installed both on ridges and valleys would be needed. Therefore, a different approach must be implemented.

APPENDIX VERIFYING THE LOWER BOUND ON CORRECTED PHASE ACCURACY
Equation (1) models the residual tropospheric error in NWP corrected interferograms. The variance of the error is obtained by expanding The temporal correlation of the reference and secondary acquisitions errors can be considered 0 everywhere since they are separated by at least six days E[ r s ] ≡ 0. On the other hand, a spatial correlation E[ A B ] = (A, B) exists for both reference and secondary acquisitions, leading to (2). According to this model, the spatial correlation function at large distances should tend to 0 (∞) = 0, and therefore, the error variance in such a regime should only depend on the accuracy of the tropospheric corrections used. Moreover, since this verification experiment is performed on a large area characterized by flat or slightly varying topography, the exponential term can be considered constant and hence absorbed by the measurement error In order to verify this, a study was carried out using Sentinel-1 interferograms with very long azimuth extension over Germany. The unwrapped phases from multiple slices were mosaicked allowing variograms at scales of up to 300 km. The variogram samples were chosen along (almost) isorange lines in order that the incidence angle not varies and the measurements projected onto the zenith direction. In this way, the effects of incidence angle could be factored out.
Analogously to [29], the zenith path delay predicted by ERA-5 was compared to GNSS-based estimates. Since the NWP model accuracy varies with geographic location, a subset of nine GNSS stations in Germany were used. The difference between the GNSS and ERA-5 ZPDs should approximate the NWP error, given that the GNSS values represent the true delay. This comparison was performed for ascending and descending geometries at the local solar times of 18:00 and 06:00, respectively. Statistics were calculated over one year of data for each geometry.
In an ERA-5 corrected interferogram, the error between two points separated by a distance d should approach the saturation value σ 2 Z = 4 σ 2 in (12) since (d) → 0 for large d and θ A = θ B = 0 at zenith. Hence, it is sufficient to confirm that the variograms saturate on average at the value 4 σ 2 . The results in Fig. 12 show an estimate of 4 σ 2 along with the variograms for ascending and descending geometries. It is interesting to note that the performance of these two geometries is comparable at 300 km. One would expect better performance for the descending pass, made in the morning when the atmosphere contains less energy. This is partially true considering midscales from 50 to 80 km where stronger turbulence due to more atmospheric power for the ascending pass leads to a larger average power and spread around the average. However, once the distance reaches the scale of the NWP model, one is limited by the error of the NWP model independently of the time of day.
An important implication of this is that a global analysis of ERA-5/GNSS deviations could provide a lower bound for the achievable accuracy in measuring large-scale displacement using InSAR.