An Accurate Method to Correct Atmospheric Phase Delay for InSAR with the ERA5 Global Atmospheric Model

Differential SAR Interferometry (DInSAR) has proven its unprecedented ability and merits of monitoring ground deformation on a large scale with centimeter to millimeter accuracy. However, atmospheric artifacts due to spatial and temporal variations of the atmospheric state often affect the reliability and accuracy of its results. The commonly-known Atmospheric Phase Screen (APS) appears in the interferograms as ghost fringes not related to either topography or deformation. Atmospheric artifact mitigation remains one of the biggest challenges to be addressed within the DInSAR community. State-of-the-art research works have revealed that atmospheric artifacts can be partially compensated with empirical models, point-wise GPS zenith path delay, and numerical weather prediction models. In this study, we implement an accurate and realistic computing strategy using atmospheric reanalysis ERA5 data to estimate atmospheric artifacts. With this approach, the Line-of-Sight (LOS) path along the satellite trajectory and the monitored points is considered, rather than estimating it from the zenith path delay. Compared with the zenith delay-based method, the key advantage is that it can avoid errors caused by any anisotropic atmospheric phenomena. The accurate method is validated with Sentinel-1 data in three different test sites: Tenerife island (Spain), Almería (Spain), and Crete island (Greece). The effectiveness and performance of the method to remove APS from interferograms is evaluated in the three test sites showing a great improvement with respect to the zenith-based approach.


Introduction
Differential Interferometric Synthetic Aperture Radar (DInSAR) and its derived Persistent Scatters Interferometry (PSI) techniques have already proven to be an extraordinary geodetic approach for large-scale and high accuracy ground deformation measurement.They can be used to monitor deformations related to earthquakes [1][2][3], tectonic movements [4,5], volcanic actions [6][7][8], landslides [9,10], and land subsidence [11,12].With the unprecedented development of SAR missions, such as Sentinel-1, TerraSAR-X, and TanDEM-X, ALOS-2, COSMO-SkyMed, RADARSAT-2, PAZ, Gaofen-3, and the planned NISAR and TanDEM-L, a large amount of SAR data are or will be available with short repeat cycles and wide swath modes.It is therefore believed that DInSAR techniques can provide an efficient and near-real-time way for monitoring dynamic processes on the Earth's surface on a global scale [13][14][15].However, atmospheric perturbations, hereafter called the Atmospheric Phase Screen (APS), remain one of the major limitations for the retrieval of geophysical signals from Interferometric Synthetic Aperture Radar (InSAR) stacks.
Initial studies [10,[16][17][18] related to atmospheric artifacts in InSAR demonstrated that atmospheric perturbations are caused by temporal and spatial changes of the atmospheric state, such as temperature, pressure, and humidity.Consequently, the refractivity of the atmosphere is not constant.When microwaves pass through the different atmospheric layers, the speed of propagation changes as the refractivity index does, and the two-way trip delay is affected as well.Meanwhile, the state of the atmosphere is temporally dependent, and it will be usually different at the times the two images of the interferogram were acquired.As a result, additional fringes would appear in the interferogram not related to either topography or deformation.Atmospheric delays thus cannot be ignored because they can often be comparable, or even much larger, in magnitude to the geophysical signals of interest.For example, the work in [17] showed that spatial and temporal changes of just 20% in relative humidity can result in 10-cm errors in the deformation products.The work in [18] reported that the atmospheric delays can be of the order of several centimeters.
To date, numerous methods have been explored to mitigate APS, all of which generally can be summarized into two categories: The first takes advantage of the APS statistical properties.On the one hand, based on the different behavior of APS and deformation signals in interferograms, i.e., turbulent atmospheric artifacts are usually highly correlated in space and uncorrelated in time, while the deformation phase component can be assumed to behave just the opposite.Therefore, a combination of spatial and temporal filters is used to separate the different phase components in PSI processing [19][20][21][22].Improved filtering approaches using auxiliary data to optimize the filtering parameters have also been proposed [23,24].Besides, the work in [25] studied the spatial autocorrelation of APS and its impact when retrieving geophysical signals.Although this approach is efficient with large datasets of SAR images under the assumption that atmospheric artifacts present a Gaussian distribution, it is still challenging in terms of filtering parameter optimizations in accordance with the atmospheric characteristics.On the other hand, in some mountainous scenarios, atmospheric artifacts can be strongly correlated with the topography [26], which is known as stratified atmosphere [18].Strategies based on empirical models have been used and further improved to remove this stratified APS, which include a linear model [27][28][29][30] and a power-law model [31].However, the other phase components often contaminate the model coefficients' estimation.The work in [32] proposed an improved linear model, which can minimize the influence from turbulent APS.Although empirical models have shown effectiveness in some case studies, they fail in situations where the APS does not follow the specific model [30][31][32].If multiple models were used, it would be almost impossible to determine automatically the boundaries of the areas where the different models should be applied.
The second category of mitigation methods utilizes various auxiliary information sources, which include Zenith Total Delay (ZTD) from GPS measurement (e.g., [33][34][35][36]), multi-spectral observations (e.g., [37,38]), and either local meteorological models (e.g., [39,40]) or Global Atmospheric Models (GAM) (e.g., [29,[41][42][43]).The main drawbacks reported from the above-mentioned references are the low spatial and/or temporal resolution and the precision of the external datasets.Although [44] proposed a framework to correct APS routinely using a GPS-based method, GPS stations are sparsely distributed or even absent in many regions.Additionally, some GPS datasets are still not freely available to the public.
Fortunately, with the development of the numerical weather prediction, GAM are able to provide accurate and higher resolution parameters for characterizing the atmospheric state, such as ERA5 data generated using Copernicus Climate Change Service Information [45].Moreover, the GAM datasets (e.g., ERA-Interim, ERA5, MERRA-2) are freely available to the public on a global scale.Consequently, GAM-based methods are considered as promising and practical for APS mitigation.
Previous studies have validated the potential ability of using GAM [29,43,46], in which the zenith path delay is estimated first.The phase delay along the Line-of-Sight (LOS) direction is then obtained geometrically by considering the incidence angle.However, in real situations, the calculation based on converting the Zenith to LOS direction (Z-LOS) can produce biases if the atmospheric delay is spatially anisotropic, especially in cases where the incidence angle is large.Some other studies [41,42,47] have introduced a Direct integration method along the LOS direction (D-LOS) based on GAM.However, comprehensive comparison analyses between the Z-LOS and D-LOS method have not been carried out in previous studies.As a result of this, the merit of the D-LOS approach has not been realized by researchers.Consequently, the Z-LOS method, which is an approximate calculation strategy compared to the D-LOS method, is still widely used in recent studies [43,48,49].
In this paper, we intend to answer whether an approximate calculation (Z-LOS), substituting the realistic D-LOS method, is acceptable or not.Firstly, detailed procedures of the D-LOS method implementation were presented, which were not provided in previous papers.After that, both Z-LOS and D-LOS methods based on the latest GAM ERA5 data were implemented for the Tenerife (Spain), Almería (Spain), and Crete (Greece) test sites with Sentinel-1 data.The performance results were compared in-depth to highlight the limitations of the Z-LOS approach.

Traditional Methods Based on GAM
Variations of temperature, pressure, and humidity in the troposphere make the refractivity index, N, not constant.It is well known that APS is caused by the velocity change of electromagnetic waves, and it is linearly dependent on N, when passing through the troposphere along the LOS direction.The tropospheric medium is mainly made up of a wet component (water vapor) and a dry (hydrostatic) component, so the refractivity index N includes both.N dry depends on the partial pressure of dry air, P, while N wet depends on the water vapor partial pressure, e.Both refractivity components depend on the temperature, T. The relation of N and the different atmospheric parameters can be expressed as [50]: where k 1 = 0.776 KPa −1 , k 2 = 0.716 KPa −1 , and k 3 = 3.75e3K 2 Pa −1 .
The atmospheric phase delay term, φ atm , of a two-way traveling path between two points, r 1 and r 2 , is proportional to the integration of the total refractivity N(T, P, e) along the line joining them: where λ is the radar wavelength.For an orbital SAR image, r 2 will be the location of the satellite when the image was acquired, r sat , and r 1 will be the ground location of the pixel for which the atmospheric artifacts are calculated, r ground .For an interferogram, the differential atmospheric delay ∆φ atm is the difference between both acquisitions.With Equation (2), APS correction methods based on meteorological models can calculate phase delays by using atmospheric parameters from local or global weather models.Global and regional atmospheric reanalysis data provide atmospheric variables of temperature, geopotential height, pressure, and humidity with relatively high temporal and spatial resolution.To calculate phase delays from these data, for each grid point, the geopotential variable is first converted to a regular vertical metric grid parameter, the geopotential height.Atmospheric parameters are then interpolated from pressure levels to altitude levels based on the geopotential height parameter.
The different methods found in the literature [40,43,46,48,51] determine the atmospheric phase delays with an approximation.Instead of integrating along the LOS direction with Equation (2), the phase delays are calculated by integrating along the zenith one.Therefore, Equation (2) becomes Equation (3), where θ is the local incidence angle (derived from the LOS direction), h ground the height of the ground location, and h top the maximum zenithal height.
Once atmospheric phase delays are integrated along a zenithal path on the grid points using Equation (3), a cubic spline interpolation in the vertical direction and a bilinear interpolation in the horizontal one are performed to obtain the phase delays for the entire SAR scene from the sparse grid points.Once having the ZTD along the vertical direction, the cosine of the incidence angle is accounted for by back-projecting the result to the LOS direction.The effectiveness of the compound method is validated in several geographic and tectonic settings [40,43,48,51].

Drawbacks of the ZTD Approach
The way of converting delay from the zenithal path to LOS is to divide it by a factor, cosθ, that depends on the local incidence angle, θ, under the assumption that the atmosphere is isotropic.In reality, this is not always true as the atmosphere is a dynamically complex system that is rarely isotropic in nature.Figure 1 provides a simple example to illustrate the heterogeneity of the atmosphere.The cloud located in vortex P clearly affects the path along LOS, while it is not present in the zenithal path.Therefore, the projection of the zenithal delay is not a good approximation of the LOS delay.In practice, water vapor, pressure, and temperature parameters can change spatially under naturally chaotic weather conditions, making isotropic assumption unacceptable.Besides, the new generation of SAR satellites can operate Terrain Observation by Progressive Scans (TOPS) acquisition mode, such as Sentinel-1 and TerraSAR-X.TOPS employs a rotation of the antenna in both the along track and cross track directions.Consequently, the incidence angle varies in a relatively large range when compared with using conventional stripmap mode.Figure 2 shows the Interferometric Wide Swath (IW) and Extra Wide Swath (EW) modes of Sentinel-1, which are implemented as three and five subswaths TOPS SAR modes.The incidence angle range for IW and EW is 29.1 • -46.0 • and 18.9 • -47.0 • , respectively [52].Hence, atmospheric heterogeneity along the zenithal and LOS direction occurs with high probability in such a big incidence angle.Consequently, the projection of ZTD to LOS delays is not a reliable approach.The incidence angle varies in a relatively big range.

Precise LOS Phase Delay Calculation
In order to avoid possible errors caused by the approximation in the Z-LOS method, the tropospheric delays actually can be estimated more accurately by integrating the refractivity along the LOS direction using Equation ( 2).
The steps of the D-LOS algorithm are the following: (1) Determination of the sampling locations along the LOS path: The locations of both ground points and zero-Doppler satellite locations are expressed in Cartesian coordinates using the World Geodetic System 1984 (WGS84) ellipsoid with ellipsoidal heights; see the geometry in Figure 3.For a general point, firstly, the three-dimensional Cartesian coordinates for the ground-range point r ground and its corresponding satellite zero-Doppler location r sat are determined.An external Digital Elevation Model (DEM) is back-projected to slant-range coordinates using the precise orbits to determine for each pixel of the image its corresponding geodetic coordinates.The LOS vector, − − → LOS, is calculated from the Cartesian coordinates of the ground point and the satellite, After calculating the LOS vector, the direction vector ûLOS , which describes the LOS direction, is calculated from the unit vector: Once the direction vector is determined, the coordinates (P ix , P iy , P iz ) for any arbitrary point P i along LOS can be easily obtained.The required sampling, ∆LOS, depends on the resolution of the available GAM data.In practice, for an efficient computation and with no loss of reliability, a LOS spatial sampling of 200 m has been set.
with i an integer.
With the above equations, locations along LOS for each ground-range point can be determined.After that, an interpolation strategy is designed to obtain the atmospheric parameters on these locations from the available GAM data.
(2) Interpolation of atmospheric parameters: What GAM datasets provide are a grid in geographic coordinates of meteorological parameters, which include pressure, temperature, and humidity.For example, with ERA5 data, sampling values of the grids are 37 pressure levels in altitude and 31 km in the horizontal.In order to obtain the required finer spatial sampling of the atmospheric state parameters on the desired locations, a cubic spline interpolation along the vertical direction and a bilinear interpolation in the horizontal plane are utilized.Additionally, as refractivity N above a reference elevation is relatively stable, locations with elevations higher than a specific threshold can be considered negligible.The determination of the reference elevation is discussed in Section 4.1.
When dealing with the vertical interpolation, the atmospheric parameters need to be interpolated with the same datum of the DEM, WGS84 in this paper.The different height types, geoidal or ellipsoidal, and all geodetic datums should be taken into consideration.One parameter included in GAM data is the geopotential, Φ(h), which is widely used in the meteorology community.It can be easily converted to geopotential height, H p , simply dividing by the gravity constant at mean sea level, g, Then, the geometric altitude H g referring to Earth's mean sea level is a function of the geopotential altitude, H p : in which E is the Earth's radius.Finally, as the DEM used considers ellipsoid heights, the ellipsoid altitude H e is obtained by adding the geoid height N obtained from the Earth Gravitational Model 1996 (EGM96) geopotential model.
Once all altitudes are ellipsoidal referred to the WGS84 datum, the atmospheric parameters can be interpolated along the vertical direction.It has to be noticed that the precise orbits of the satellites are also in Cartesian coordinates referred to WGS84, and so, all calculations are unified to the same framework.

Tenerife Island, Spain
Tenerife island is situated in the Atlantic Ocean opposite the northwestern coast of Africa.It is a rugged and volcanic island sculpted by successive eruptions throughout its history.Figure 4a shows the location of the test site.The center of the island is dominated by the volcano Teide.With an elevation of 3718 m above sea level, it is the third largest volcano in the world.A total of 15 eruptions have been historically recorded on the island [53].Among various geodetic monitoring tools spread on this island, InSAR techniques have demonstrated their capacity for volcano monitoring purposes [54][55][56].However, as the topography variation (from sea level to 3718 m) on Tenerife is significant jointly with the proximity of the sea, most differential interferograms are affected by both stratified and turbulent atmospheric artifacts.The work in [57] applied a method using the Weather Research and Forecast (WRF) model [58] to mitigate atmospheric effects, but without considering the realistic phase delay calculation along LOS.Their results showed WRF as a promising tool for APS mitigation while failing at anisotropic atmospheric situations.

Crete, Greece
The island of Crete is located in the Eastern Mediterranean (Figure 4), which lies within the outer fore-arc of the largest and most active subduction zone in Europe.PSI techniques have been explored to detect the active vertical surface deformation patterns across Crete [59].However, its offshore location and significant topography incur atmospheric artifacts, which need to be considered.A previous study [60] showed that Crete is an ideal study area, because it represents a complicated set of tropospheric conditions involving maritime-terrestrial heating and moisture contrasts combined with extensive topographical variations.

Almería, Spain
Almería basin lies in one of the most arid parts of Europe, receiving an average rainfall of only 250 mm per year [61].Given historical patterns of agricultural and urban development, the water demand has been largely supplied from ground water, leading to the situation of water overexploitation.Consequently, Almería basin suffers ground subsidence induced by aquifer exploitation.PSI has become the optimum technique for periodically monitoring subsidence.A study from [61] by using DInSAR revealed that between 2003 and 2009, this region suffered a subsidence rate ranging from 1.7-5 mm/year.Unfortunately, interferograms in this area are also vulnerable to APS contamination due to the strong topography and its coastal location (Figure 4c), which may lead to unreliable estimations of subsidence patterns and values.

European Center for Medium-Range Weather Forecasts (ECMWF) Datasets
In this paper, the advanced D-LOS algorithm was validated using the ERA5 dataset, which is a global atmospheric reanalysis product provided by the ECMWF.ERA5 is the latest generation of atmospheric reanalysis data after ERA-Interim generated using Copernicus Climate Change Service Information [45].This product covers the period from 1979 to the present as ERA-Interim covers.Compared with the 6-h temporal and 79-km spatial sampling of ERA-Interim data, ERA5 provides a much higher resolution in both time and space.Hourly reanalysis data are available at a horizontal resolution of 31 km with ERA5 data.The higher temporal resolution minimizes the errors due to the time difference between SAR acquisition and meteorological data.The nearest ERA5 data to the SAR images acquisition times were selected.Another advantage of ERA5 data is that uncertainty estimates for all parameters are available at three-hour intervals and at a horizontal resolution of 62 km.The uncertainties can be used as reliability weights when interpolating the atmospheric parameters, as well as potentially applied to evaluate the trustworthiness of the calculated APS maps.

SAR Datasets and Data Processing
The performance of the Z-LOS and D-LOS algorithms were evaluated based on Sentinel-1 SAR datasets from the European Space Agency (ESA).In the three test sites, the selected SAR images were acquired under complex weather conditions from different seasons in order to validate the general applicability of the corrections.Differential interferograms were generated using SUBSIDENCE-GUI software, developed by UPC [22], with USGS SRTM DEMs and precise orbits from the Copernicus POD (Precise Orbit Determination) Service to remove topographic phase components.As values in a wrapped interferogram are limited to −π and π, a phase-elevation trend cannot be clearly seen using wrapped phases.In this paper, unwrapped phases, i.e., after solving phase ambiguities (2π jumps), were used.All differential interferograms were unwrapped using the SNAPHU package [62].

Algorithm Validation in Tenerife Island
The D-LOS algorithm was validated demonstratively in Tenerife island using Sentinel-1 and ERA5 data, with explicit comparison with Z-LOS.As can be seen in the following sections, atmospheric artifacts presented in interferograms are not correlated with topography in a typical linear or power-law model, empirical model methods cannot estimate the APS correctly.Besides, it is the stratified APS, rather than the turbulent component, that dominates the atmospheric artifacts; therefore, classical filtering methods have a limited ability to compensate the APS.The SAR images in Tenerife were acquired at 07:02 in the morning, very close to the hourly ERA5 reanalysis data.The atmospheric conditions between SAR images and ERA5 data were assumed functionally simultaneous, as each was separated by only two minutes.

Implementation of the Improved Method
(1) Cubic spline interpolation of atmospheric parameters: Although ECMWF provides 37 pressure levels reanalysis data higher than 40 km above the sea level, in practice beyond a certain elevation, the parameters can be considered constant.A cubic spline interpolation was applied on grid points from sea level to a reference elevation for the temperature, the pressure, and the relative humidity along the vertical direction.Figure 5 visualizes the vertical interpolation results for different parameters at Grid Point A (see Figure 4a).It is clear that the curves for all parameters were smooth, with no abrupt changes, and so, the cubic spline interpolation was adequate.The interpolation results at Point A represent a general situation of the low-pass behavior of the parameters.The interpolated data along the vertical direction will further provide input data for the horizontal interpolation on the desired locations along the LOS direction.(2) Maximum altitude determination: As mentioned above, the refractivity of tropospheric changes in both temporal and spatial domains.However, since the interferometric phase is the pixel-wise difference between two SAR acquisitions, the atmospheric delays depend on the refractivity variations rather than their singly-imaged absolute values.Other meteorological measuring tools such as sounding measurements provide detailed vertical profile of the atmospheric state.The work in [18] analyzed radiosonde data in De Bilt, the Netherlands, and found that the variability of the refractivity above 5 km can be considered negligible.The work in [43] assumed negligible effects of atmospheric stratification above 30 km.The work in [31] studied weather balloon data and found that relative delays above a 7-13 km height did not differ significantly.In the Tenerife island case and for Grid Point A, 153 relative phase delays obtained from all possible combinations of 18 ERA5 data from different days are displayed as a scatter plot in Figure 6.It can be clearly seen that the relative phase delays converged to zero with the increasing of altitude.In order to optimize the integration of refractivity along the LOS direction, it is not necessary to consider all elevations, but only those that affect the relative phase delays.The maximum altitude to consider, h top , can be set to those for which the standard deviation of the relative delays is smaller than 1 rad.The maximum altitude value to be considered in the integration for the Tenerife case can be set to 28 km.

Case Study Discussion of Tenerife Island, Spain
In order to reduce as much as possible the phase contributions caused by DEM inaccuracies and deformation in the interferometric phase, the spatial perpendicular baselines (separation of the orbits) and temporal baselines (time difference of two acquisitions) were limited to be less than 50 m and 36 days, respectively.In addition, as can be seen in Table 1, the interferograms were all made from two acquisitions covering two consecutive months to minimize the temporal baseline and, at the same time, maximize the phase quality.The one-year time span also guaranteed a diversity of climate and weather conditions.Figure 7 shows all the selected differential interferograms covering the period from September 2015 to September 2016.At first glance, the phases were strongly correlated with topography.This phase pattern can also be caused either by errors on the DEM or the orbital information.In order to exclude this factor, the relationship between the number of fringes and the spatial baseline was analyzed.Figure 8 shows a scatter plot of the interferograms, in which it is evident that the number of fringes did not follow a linear trend with the spatial baselines.To underscore the complexity of atmospheric artifacts in Tenerife, a scatter plot relating the height and unwrapped original phase of 105,137 selected Persistent Scatters (PSs) for the interferogram 20160721-20160814 (Ifg11) is shown in Figure 9f.From the scatter plot, it is clear that there is no strict linear or apparent power-law correlation between atmospheric artifacts and topography; therefore, correction methods based on empirical models may fail for interferograms having similar characteristics as this one.Figure 9 shows the correction results of Ifg11 by using different methods, in which the first and second rows present the original interferogram, the predicted APSs, and the corresponding residuals.The unwrapped phase-elevation scatter plots are shown in the third and fourth rows.Atmospheric artifacts can be clearly seen in the unwrapped original interferogram (Figure 9a), with phase patterns increasing from the coast (low elevation) towards the center of the area (top of the mountain).The short temporal baseline also ensures that the phase disturbances due to localized terrain deformation are extremely small and therefore do not affect the analysis of the overall behavior of the atmospheric compensation.With the traditional zenith to LOS method (Z-LOS), i.e., zenith path delay is first calculated and then back-projected to LOS divided by a factor cosθ, APS can be mitigated to some extent.Figure 9b,c shows the estimated Z-LOS APS and its residual phase after APS correction, respectively.From the predicted APS, a similar delay trend was observed across the whole interferogram.After correcting with Z-LOS delays, long-wavelength tropospheric delays can be removed.As can be seen from the residual phase, the phase dispersion was reduced.However, the Z-LOS method failed in places where the topography changed abruptly.As a result, the residual phase (Figure 9c) exhibited a large variation around the top of the mountain, and a significant residual phase appeared along the SW to NE direction in the middle of the interferogram.Besides, the corresponding phase-elevation scatter plots of predicted APS (Figure 9g) and the residual phase (Figure 9h) also proved the limitations of Z-LOS method in estimating delays, especially at elevations between 2000 m and 3700 m, where the phases varied significantly.Looking at the residual phase-elevation plots (Figure 9h), phases fluctuated around zero, but with a wide variation.From the statistical point of view, the Z-LOS method reduced for this particular case the Standard Deviation (SD) from 6.79 to 3.83 rad.
With the D-LOS method, both the predicted APS and residual phase showed a better performance when compared with Z-LOS.Specifically, D-LOS predicted APS presented a very similar phase pattern with the original interferogram in terms of magnitude and location, especially at the center of the interferogram.The uniformity of the residual phase in Figure 9e further demonstrated its excellent performance.Compared with the Z-LOS results, the residual interferogram presented a smaller phase dispersion.The clear improvement can also be seen when comparing the estimated APS phase-elevation scatter plots (Figure 9i) with the original one (Figure 9f), as they both showed an almost identical trend.A further inspection of the phase residues in Figure 9j showed that the residual phases were distributed around zero with a narrow variation.After applying the D-LOS correction, the SD was reduced by a 61%, down to 2.62 rad.

Statistical Analyses Discussion
To evaluate whether the improvements shown in the above case can also be achieved in a more general situation, statistical analyses of 12 short temporal baseline interferograms in Tenerife covering an entire year were carried out.In Table 1, the Original column shows a large phase SD for each differential interferogram, due to the strong APS contamination.After applying the Z-LOS correction method, all interferograms except Ifg9 presented a phase SD reduction.For Ifg9, the SD even increased from 2.98 to 3.62 rad.In other words, the Z-LOS method cannot correct the APS properly in this specific case.In the best case (Ifg4), the Z-LOS approach achieved a 71% SD reduction.On the other side, the D-LOS method's performance was comprehensively better than the Z-LOS method.For all pairs, D-LOS can achieve a considerable SD reduction.The best case occurred in Ifg4, reducing the SD of 77% from 5.99 to 1.40 rad, while even in the worst case (Ifg9), it provided a 23% reduction.The statistical SD values of the original, Z-LOS residual, and D-LOS residual in Table 1 are visualized in Figure 10.Once again, it is clear that both Z-LOS and D-LOS can mitigate APS delay, but D-LOS consistently outperformed the Z-LOS approach.

Applications to Different Sites
In order to further validate the performance of the D-LOS method for APS mitigation, it was tested on another two sites: Crete (Greece) and Almería (Spain).

Results and Discussion in Crete, Greece
In Crete, strong variation of topography and an extensive coastline are prone to suffering tropospheric phase delays.From the one-year dataset of Sentinel-1 images, eight representative interferograms, each one covering two consecutive months (see Table 2) were generated.They represented almost one year of climate conditions except the period from December to March.The interferograms from December to March were discarded as the area was affected by snow on the Lefka Ori mountain, making the differential phases useless.It is also noted that in this area, the time difference between Sentinel-1 acquisition and ERA5 reanalysis was 24 min.This is acceptable compared with other global reanalysis data (e.g., ERA-Interim, MERRA) with a 6-h interval.In Figure 11, one interferogram (Ifg8), formed from two acquisitions during the wet season, was selected to illustrate the performance of the D-LOS algorithm.The short temporal and spatial baseline precluded deformation and DEM phase errors, so the primary phase component can be considered as atmospheric artifacts in the original unwrapped interferogram (Figure 11a).The main visible feature was the phase trend from low elevations to the top of the mountains.This feature was more visible looking at the phase-elevation scatter plots in Figure 11f.
With the Z-LOS method, this feature can be partly estimated (Figure 11b), and so, the dominant trend can be removed with the correction, as can be seen in Figure 11c.However, paying attention to certain areas, the Z-LOS method over-predicted the tropospheric delay locally.The scatter plots of Z-LOS predicted and corrected phases shown in Figure 11g and Figure 11h further proved this.From these plots, it is clear that the Z-LOS predicted phase delay showed a large variance, producing relatively large fluctuations in the residual phases.Statistically, the phase SD was reduced by 19% from 2.88 to 2.34 rad.With the D-LOS algorithm, the patterns of predicted APS in Figure 11d showed a good correspondence with the original interferogram, both in location, as well as amplitude.The D-LOS corrected interferogram (Figure 11e) showed relatively small phase values, fluctuating around zero.The observed D-LOS phase-elevation pattern shown in Figure 11i matched the features of the original interferogram.Additionally, it is important to highlight the smaller variance in the residual phase compared with the Z-LOS one.The good correspondence between the original interferogram and D-LOS prediction was also confirmed by a significant SD reduction in the phase residues, by 28% from 2.88 to 2.07 rad.Additionally, statistical analyses for all eight interferograms were carried out and listed in Table 2.In this test site, both Z-LOS and D-LOS methods can achieve variance reductions for all interferograms.Moreover, once again, the D-LOS algorithm outperformed the Z-LOS approach for each individual pair in terms of SD reduction.For Z-LOS, the best correction case was achieved for Ifg2, showing a 48% variance reduction; while Ifg1 was the best case for D-LOS, where the phase SD was reduced by 74%.For both methods, the worst correction was for Ifg3, with reductions of 5% and 15%, respectively.Figure 12 shows graphically the performance of both methods detailed in Table 2.

Results and Discussion in Almería, Spain
Almería is a region located in the southeast of the Iberian Peninsula.Its southeast coastline opens to the western Mediterranean Sea, while its interior northwest is mountainous.It is interesting to highlight that, in the center of the study area, there is a basin surrounded by mountains, so the atmospheric conditions were expected to be complicated due to the Mediterranean climate and the presence of mountains.To evaluate the performance of APS correction, twelve interferograms covering one year with short temporal spans were selected in this test site.Among them, the pair 20171031-20171124 (Ifg11) was explicitly analyzed as a typical example (Figure 13).
The time difference between ERA5 and Sentinel-1 acquisition was less than 10 min in this case.Again, it can be assumed that the atmospheric conditions would be very similar in such a short time difference.Figure 13a shows how the differential interferometric phase due to the atmospheric artifacts was strongly correlated with the topography.The Z-LOS predicted phase, shown in Figure 13b, was able to retrieve partially the trend of the original interferogram, especially in the flat areas, but it failed in some of the mountainous locales.The bad performance of the estimated APS was better highlighted with the residual phases shown in Figure 13c.The phase SD was reduced only by 8%, from 2.54 to 2.33 rad, with the Z-LOS correction.A significant improvement in the performance can be observed with the D-LOS approach.Specifically, the D-LOS prediction shown in Figure 13d exhibited a very similar phase pattern with the original one in terms of the shape and amplitude.As Figure 13e shows, the tropospheric delays all over the interferogram were reduced except in some localized areas located near the top of some mountains.With the D-LOS method, the phase SD was reduced by 39% to 1.79 rad.Both methods were applied to all the interferograms, and the performance results are listed in Table 3 and graphically illustrated in Figure 14.It is clear from Figure 14 that, for most interferograms, the Z-LOS mitigation reduced the tropospheric delays, except for Ifg3 and Ifg7.For these two cases, the SD values actually increased.After correcting with the D-LOS method, all interferograms showed an obvious SD reduction.Overall, compared with the average 14% mitigation achieved by the Z-LOS method, the D-LOS algorithm was able to achieve a 29% reduction on average.

Discussion
Overall, the realistic D-LOS method outperformed the Z-LOS method in all three test sites.On the one hand, from the statistical point of view for all the interferograms, D-LOS approach presented more significant variance reduction than using the Z-LOS method.In Tenerife island, D-LOS can achieve a 23-77% SD reduction for twelve interferograms with an average level of 50%, while Z-LOS achieved only 14-71% and an average decrease of 37%.In the Crete, Greece site, after applying D-LOS correction, the SD decreased around 15-74%, with an average reduction of 38%, compared with 5-48% and the corresponding average level of 22% for the Z-LOS method for eight interferograms.Over Almería, the Z-LOS method achieved a SD reduction of 14% on average, while the D-LOS strategy showed a better performance, reaching 29% on average for twelve interferograms.On the other hand, D-LOS approach succeeded to mitigate atmospheric artifacts in some specific interferograms, where the SD of residual phases even increased with Z-LOS approach.
In principle, the APS correction methods based on GAM data had the feasibility to compensate APS globally under all weather conditions.This was proven in the Tenerife, Almería test site with the interferograms covering one entire year and Crete with interferograms covering eight months.The reason for discarding the interferograms from December to March in Crete was that the interferometric phases themselves in these interferograms were random.The random characteristics were due to the scattering mechanism of the snow coverage in SAR images.Although the GAM-based methods were promising for APS compensation, it is still noted that the performance was dependent on the quality of the weather forecast data.In practice, it is difficult to assess whether the weather reanalysis data are reliable enough or not.For the APS compensation purpose, SD reduction can be used as a metric to evaluate its performance.
It is noted that both Z-LOS and D-LOS strategies can remove long-wavelength tropospheric delays partially.Consequently, SD reductions can be seen with both methods for most interferograms.Moreover, the D-LOS approach demonstrated its robustness in correcting atmospheric artifacts in some localized areas (e.g., areas with high-elevation and/or surface heterogeneity) where the Z-LOS method failed.This merit was proven by exhibiting more SD reductions after D-LOS correction than the Z-LOS method.

Conclusions
In this paper, a realistic calculation strategy based on atmospheric reanalysis data was implemented to mitigate atmospheric artifacts for InSAR.Compared to the conventional Z-LOS method, the advanced D-LOS approach calculated the tropospheric delay along the real LOS path rather than an approximation.Detailed comparison analyses between the two methods were carried out over Tenerife, Crete, and Almería with Sentinel-1 data.The advanced D-LOS method in all three test sites demonstrated its robustness and better performance than the Z-LOS method.
The APS correction technique based on GAM data can be applied to all weather conditions globally without any empirical model assumption.With the development of GAM, on account of its global coverage, both high temporal and spatial resolution and more accurate atmospheric parameters, the compensation method is very promising for near real-time InSAR applications.
As atmospheric artifacts from differential interferograms are in great agreement with APS from atmospheric reanalysis data, it demonstrates a potential ability to study atmospheric dynamics using InSAR.On the one hand, it is possible to retrieve atmospheric water vapor mapping from differential interferograms.Compared with GPS data, InSAR is a promising technique to retrieve the water vapor product with global coverage and higher resolution.On the other hand, the InSAR technique can also provide potential enrichment of datasets used for research using numerical weather models, especially in very localized areas for turbulent atmosphere research.Even now, the latest suite of numerical weather prediction models have difficulty in accurately predicting meso-and micro-scale atmospheric dynamics, due mainly to computational discretization and observational scarcity.It is believed that atmospheric delays from InSAR can provide a promising observation dataset for weather data assimilation.

Figure 1 .
Figure 1.An example of realistic atmospheric heterogeneities along the LOS path.The cloud clearly affects the LOS direction, but not in the zenithal one.ZTD, Zenith Total Delay.

Figure 2 .
Figure 2. Sentinel-1 Interferometric Wide Swath (IW) and Extra Wide Swath (EW) acquisition modes.The incidence angle varies in a relatively big range.

Figure 3 .
Figure 3. Geometry of the LOS path in Cartesian coordinate system.

Figure 4 .
Figure 4. Overview of geographic locations of the test sites: Tenerife island, Spain, Almería, Spain, and Crete, Greece.(a), (b), and (c) show the local topography of each study area.

Figure 5 .
Figure 5. Spline interpolation for temperature, pressure, and water vapor pressure along the vertical direction.

Figure 6 .
Figure 6.Relative tropospheric phase delay along the vertical direction.The darker the color is, the more the scatters overlap.

Figure 7 .Figure 8 .
Figure 7. Original wrapped differential interferograms covering one entire year over Tenerife island, Spain.Each interferogram connects two consecutive months.Detailed temporal and spatial baselines information for each interferogram can be found in Table1.

Figure 9 .
Figure 9. 20160721-20160814 (Ifg11) case study of APS compensation in Tenerife.(a-e) represent the original interferogram, Z-LOS predicted APS, residue phase after Z-LOS correction, D-LOS predicted APS, and the residue phase after D-LOS correction, respectively.(f-j) are the phase-elevation scatters for the original phase, Z-LOS predicted APS, residue phase after Z-LOS correction, D-LOS predicted APS, and the residue phase after D-LOS correction, respectively.

Figure 10 .
Figure 10.SD comparison of original interferograms (magenta bars), residual phases corrected from the Z-LOS method (blue bars), and residual phases after D-LOS correction (orange bars) on all twelve pairs in Tenerife.

Figure 11 .
Figure 11.20171126-20171208 (Ifg8) case study of APS compensation in Crete.(a-e) represent the original interferogram, Z-LOS predicted APS, residue phase after Z-LOS correction, D-LOS predicted APS, and residue phase after D-LOS correction, respectively.(f-j) are the phase-elevation scatters for the original phase, Z-LOS predicted APS, the residue phase after Z-LOS correction, D-LOS predicted APS, and the residue phase after D-LOS correction, respectively.

Figure 12 .
Figure 12.SD comparison of original interferograms (magenta bars), residual phases corrected from Z-LOS method (blue bars), and residual phases after D-LOS correction (orange bars) on all eight pairs in Crete.

Figure 13 .
Figure 13.20171031-20171124 (Ifg11) case study of APS compensation in Almería.(a) The original interferogram.(b,c) The Z-LOS predicted APS and corrected interferogram with Z-LOS method.(d,e) The D-LOS predicted APS and corrected interferogram with the D-LOS method.

Figure 14 .
Figure 14.SD comparison of original interferograms (magenta bars), residual phases corrected from the Z-LOS method (blue bars), and residual phases after D-LOS correction (orange bars) on all twelve pairs in Almería.

Table 1 .
Statistical comparison of different mitigation methods for 12 Sentinel-1 interferograms used in the Tenerife island.Values in parentheses are the correction percentage of (original-residual)/original.Ifg, Interferogram.

Table 2 .
Statistical comparison of different mitigation methods for 8 Sentinel-1 interferograms used in Crete, Greece.Values in parentheses are the correction percentage of (original-residual)/original.

Table 3 .
Statistical comparison of different mitigation methods for 12 Sentinel-1 interferograms used in Almería, Spain.Values in parentheses are the correction percentage of (original-residual)/original.