Correcting InSAR Topographically Correlated Tropospheric Delays Using a Power Law Model Based on ERA-Interim Reanalysis

Tropospheric delay caused by spatiotemporal variations in pressure, temperature, and humidity in the lower troposphere remains one of the major challenges in Interferometric Synthetic Aperture Radar (InSAR) deformation monitoring applications. Acquiring an acceptable level of accuracy (millimeter-level) for small amplitude surface displacement is difficult without proper delay estimation. Tropospheric delay can be estimated from the InSAR phase itself using the spatiotemporal relationship between the phase and topography, but separating the deformation signal from the tropospheric delay is difficult when the deformation is topographically related. Approaches using external data such as ground GPS networks, space-borne spectrometers, and meteorological observations have been exploited with mixed success in the past. These methods are plagued, however, by low spatiotemporal resolution, unfavorable weather conditions or limited coverage. A phase-based power law method recently proposed by Bekaert et al. estimates the tropospheric delay by assuming the phase and topography following a power law relationship. This method can account for the spatial variation of the atmospheric properties and can be applied to interferograms containing topographically correlated deformation. However, the parameter estimates of this method are characterized by two limitations: one is that the power law coefficients are estimated using the sounding data, which are not always available in a study region; the other is that the scaled factor between band-filtered topography and phase is inverted by least squares regression, which is not outlier-resistant. To address these problems, we develop and test a power law model based on ERA-Interim (PLE). Our version estimates the power law coefficients by using ERA-Interim (ERA-I) reanalysis. A robust estimation technique was introduced in the PLE method to estimate the scaled factor, which is insensitive to outliers. We applied the PLE method to ENVISAT ASAR images collected over Southern California, US, and Tianshan, China. We compared tropospheric corrections estimated from using our PLE method with the corrections estimated using the linear method and ERA-I method. Accuracy was evaluated by using delay measurements from the Medium Resolution Imaging Spectrometer (MERIS) onboard the ENVISAT satellite. The PLE method consistently delivered greater standard deviation (STD) reduction after tropospheric corrections than both the linear method and ERA-I method and agreed well with the MERIS measurements.


Introduction
Interferometric synthetic aperture radar (InSAR) is a global all-weather remote sensing technique that can generate a high-precision surface deformation map with a high spatial resolution for large areas.Typically, a standard frame of ENVISAT ASAR in strip map mode covers approximately 100 km 2 with 30 m resolution.Atmospheric delay error, however, poses one of the major challenges to InSAR geodesy, especially for monitoring deformation signals with small magnitude and long wavelength.These deformation signals could be due to regional surface subsidence [1], landslides [2], volcanic eruptions [3], tidal loading [4], tectonic slow slip [5], and earthquakes [6].Atmospheric delays are typically introduced by the ionosphere and the troposphere when radio signals propagate through the atmosphere.Although the observed atmospheric delay in InSAR is due to these two terms, the major challenge mainly comes from tropospheric delay because of its high spatiotemporal variability.Spatiotemporal variations in the free electrons of the ionosphere (≥60 km) cause phase advances of the radio signals that are only significant at higher latitudes and for larger wavelength radar signals, such as the L and P bands [7].Tropospheric effects are due to variations in pressure, temperature, and relative humidity in the lower part of the troposphere (≤12 km) that often lead to phase delays up to a few decimeters in magnitude and overwhelm the deformation signal of interest [8,9].This study focuses on estimating topographically correlated tropospheric delay.We were not concerned with ionospheric effects because we used only C-band SAR images in our test.
Many methods have been developed to estimate the tropospheric delay on InSAR [10][11][12].In general, the approaches for mitigating tropospheric delay can be classified into two categories: those using ancillary delay measurements and empirical models based on the interferometric phase itself.
The methods based on ancillary datasets produce synthetic delay maps directly by external sources such as GPS wet delay measurements [13][14][15], multispectral imagery such as Moderate Resolution Imaging Spectroradiometer (MODIS) and Medium Resolution Imaging Spectrometer (MERIS) [16,17], GPS integrated with spectrometer data [18], and local or global meteorological reanalysis products [19][20][21].The advantages of these methods are their straightforwardness and the fact that they do not need a priori information about surface deformation and assumptions regarding the temporal correlation of the troposphere.However, methods using ancillary delay measurements are limited by the spatiotemporal resolution and coverage of the data source.Furthermore, weather conditions and sunlight dependency may reduce the quality of the multispectral imagery available for the study region.Following recent advances in the global atmospheric model (GAM), atmospheric models are now regarded as the most promising method of tropospheric correction for InSAR [22][23][24][25][26]. GAM provides atmospheric parameters (pressure, temperature, and humidity) that can be applied to calculate the refractivity equations and determine path delays.The distinctive features of GAMs include their availability under all weather conditions regardless of cloud presence and relatively high spatial resolution.The ERA-Interim (ERA-I) has a resolution of ~75 km [27] and the Weather Research and Forecasting (WRF) model has a resolution of ~7 km [28].Nevertheless, GAMs do not meet the high-resolution requirements for InSAR applications and are not acquired simultaneously with SAR data.Temporal, lateral, and vertical interpolations of atmospheric parameters are necessary for interferograms.Data interpolation also introduces uncertainties into the corrections.Time differences are not a problem for MERIS because they are acquired simultaneously with ASAR images.Thus, MERIS has the best accuracy and can be used to evaluate the performance of other tropospheric correction methods in ASAR interferograms.
The empirical methods calculate the tropospheric delay from the interferometric phase.The tropospheric delay ∆φ trop for an individual interferogram is calculated by assuming a linear model (∆φ trop = Kh + φ c ) [29,30] between the elevation h and the interferometric phase.K is the scaled factor that needs to be estimated and φ c represents a constant shift that is applied to the entire interferogram.Several techniques in linear models have been developed to separate the tropospheric signals from residual orbits, deformation signals, or other phase noises, such as removal of a preliminary estimate of the deformation signals [31,32], evaluation of a local phase-topography relationship by spatial band-filtering that is insensitive to deformation signals [33], and spatially variable consideration for a piecewise slope correction over multiple windows [34].The relationship between phase and topography using such empirical methods, however, depends on the spatial extent of the SAR scene, which sometimes leads to wrong estimates of spatial variation in the tropospheric stratification [26].A power law model (∆φ trop = K (h 0 − h) α ) developed by Bekaert et al. [35] can account for the spatially variable signal in tropospheric delay and can be applied to a region containing topographically correlated deformation.The two input coefficients h 0 and α are estimated from balloon sounding data, but these data are not always available for a region of interest.These coefficients are empirical and statistical results and not values estimated for each specific SAR acquisition.These limitations have been discussed by Bekaert et al. [35] who suggested using weather model data instead of sounding data to estimate the spatiotemporal variable coefficients, which could further improve the tropospheric delay estimation.Moreover, the success of the power law method is limited to cases where the relative delays can be evaluated using the average value of master and slave acquisitions.
In our study, we introduced ERA-I reanalysis products obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) instead of sounding data to estimate the coefficients in a power law model at each SAR acquisition.In our algorithm, which is called the PLE method, we assume a power law relationship between tropospheric delay and topography.In this instance, the coefficients α (interferometric tropospheric delay trend) and h 0 (constrained height) are sensitively simulated from the relative tropospheric delay on individual interferograms.Additionally, we introduced an expanded robust method to estimate the scaled factor between band-filtered topography and phase over multiple windows throughout each interferogram.As an outlier-resistant regression, the robust method is better than the least squares regression used in the original power law method.Finally, scaled factor values were multiple weighted and extrapolated to each pixel.We tested the proposed PLE method on ENVISAT ASAR images collected over Southern California, United States and Tianshan, China.Then, we compared these results with those from the conventional linear method and the independent ERA-I method.The accuracy was evaluated and validated with MERIS measurements.Our PLE method was derived from the power law method included in the toolbox for reducing atmospheric InSAR noise (TRAIN) [36].In addition, we used TRAIN to calculate the linear method, ERA-I and MERIS corrections.

Preparation for InSAR Data and ERA-I Reanalysis
We investigated the correction capabilities of the PLE method on two test regions, Southern California and Tianshan.In California, we selected eight descending SAR images from the ENVISAT archives of the European Space Agency centered at 34 • N, 117.8 • W, covering an area of ~70 km by 100 km.These images were taken between 12 July 2008 and 25 September 2010 and were acquired on track 170 at 18:01 UTC (Figure 1a).The images extend from the San Gabriel Mountains (range of ~3 km elevation) in the north to the Los Angeles Basin area in the south.Significant and vertical surface displacements were revealed in the basin area (southern part of the interferogram) [37].For the Tianshan region, we selected seven descending ENVISAT ASAR images on track 191, spanning the period from September 2007 to June 2010 and covering an area of ~80 km by 200 km (Figure 1b).The topography ranges from approximately 1 km to 6 km, extending from the South Tianshan Orogenic Belt to the Kashi area.The Kashi Depression is one of the most complex active tectonic areas in the southern flank of Tianshan, China.In the interseismic deformation rate map showing the Kashi Depression area between 2003 and 2007 obtained by InSAR time series analysis [38], no significant earthquakes occurred during our study period.
By selecting two different study areas, we can test and assess our PLE method under different geographical regions and atmospheric conditions.The test region in Southern California lies to the east and northeast of the Pacific Ocean and therefore has a humid atmosphere rich with atmospheric water vapor, which leads to a dominant stratified delay in the interferogram.Meanwhile, the test region in Tianshan in Northwestern China is an inland plateau with a semi-arid environment and a complicated atmospheric process.The atmospheric condition in this region is more turbulent than that in Southern California.We will discuss the performance of our PLE method under these two atmospheric conditions in Section 3. In addition, the deformation signals in the two study regions have different features: the deformation in California is mainly local, small-scale (<5 km) land subsidence in Los Angeles whereas, in Tianshan, the deformation is mainly large-scale (>10 s km) tectonic deformation.The different deformation characteristics also affect the band filtering, leading to selection of different bands, as discussed in Section 2.3.
For interferometric processing, we used the Repeat Orbit Interferometry Package (ROI_PAC) to focus the raw data and used DORIS to generate interferograms [39,40].During processing, we used the precise ODR orbits from the Delft Institute for Earth-Oriented Space Research to minimize the orbital errors and the Shuttle Radar Topography Mission (SRTM) digital elevation model (30 m) to remove topographic phase [41].The differential interferogram was corrected with a plane trend in the range and azimuth directions for residual orbital errors.All interferograms were multi-looked by 40 looks in azimuth and 8 looks in range to reduce phase noise and enhance phase unwrapping accuracy.Before phase unwrapping, the adaptive power spectrum filter was applied to improve the signal-to-noise ratio of the interferograms [42].Finally, we obtained five interferograms in California and four interferograms in Tianshan for our case study and show the spatiotemporal baseline and coherence parameters in Table 1.We then used a coherence threshold of 0.2 to select coherent pixels, i.e., pixels with coherence values lower than 0.2 were discarded.Our test cases had small land deformation signals because of the short temporal baseline.Therefore, the dominant signals in the interferograms are contributed by tropospheric delay.water vapor, which leads to a dominant stratified delay in the interferogram.Meanwhile, the test region in Tianshan in Northwestern China is an inland plateau with a semi-arid environment and a complicated atmospheric process.The atmospheric condition in this region is more turbulent than that in Southern California.We will discuss the performance of our PLE method under these two atmospheric conditions in Section 3. In addition, the deformation signals in the two study regions have different features: the deformation in California is mainly local, small-scale (<5 km) land subsidence in Los Angeles whereas, in Tianshan, the deformation is mainly large-scale (>10 s km) tectonic deformation.The different deformation characteristics also affect the band filtering, leading to selection of different bands, as discussed in Section 2.3.
For interferometric processing, we used the Repeat Orbit Interferometry Package (ROI_PAC) to focus the raw data and used DORIS to generate interferograms [39,40].During processing, we used the precise ODR orbits from the Delft Institute for Earth-Oriented Space Research to minimize the orbital errors and the Shuttle Radar Topography Mission (SRTM) digital elevation model (30 m) to remove topographic phase [41].The differential interferogram was corrected with a plane trend in the range and azimuth directions for residual orbital errors.All interferograms were multi-looked by 40 looks in azimuth and 8 looks in range to reduce phase noise and enhance phase unwrapping accuracy.Before phase unwrapping, the adaptive power spectrum filter was applied to improve the signal-to-noise ratio of the interferograms [42].Finally, we obtained five interferograms in California and four interferograms in Tianshan for our case study and show the spatiotemporal baseline and coherence parameters in Table 1.We then used a coherence threshold of 0.2 to select coherent pixels, i.e., pixels with coherence values lower than 0.2 were discarded.Our test cases had small land deformation signals because of the short temporal baseline.Therefore, the dominant signals in the interferograms are contributed by tropospheric delay.
(a) California   Global and regional reanalysis products provide estimates of atmospheric parameters several times a day at different pressure levels.We used the free archived ERA-I model, the latest global atmospheric reanalysis product from ECMWF, which is produced based on a 4D variable assimilation of global surface and satellite meteorological data.ERA-I provides data (pressure, temperature, and relative humidity) at a spatial resolution of 0.7° grid (~75 km), at a 6-h interval (0:00, 6:00, 12:00 and 18:00 UTC) daily, and on 37 pressure levels from 1989 to the present [27].Figure 2 shows the distribution of ERA-I grid points in the two study areas.Global and regional reanalysis products provide estimates of atmospheric parameters several times a day at different pressure levels.We used the free archived ERA-I model, the latest global atmospheric reanalysis product from ECMWF, which is produced based on a 4D variable assimilation of global surface and satellite meteorological data.ERA-I provides data (pressure, temperature, and relative humidity) at a spatial resolution of 0.7 • grid (~75 km), at a 6-h interval (0:00, 6:00, 12:00 and 18:00 UTC) daily, and on 37 pressure levels from 1989 to the present [27].Figure 2 shows the distribution of ERA-I grid points in the two study areas.

Tropospheric Delay in InSAR Data
A tropospheric delay can be split into hydrostatic, wet, and liquid components depending on variations in atmospheric conditions.The wet component is the major limiting factor due to water vapor instability in space and time (approximately 20 mm to 300 mm) [43].Unlike the wet delay, temperature and pressure are smooth in space, leading to a better-resolved, longer wavelength hydrostatic delay, accounting for about 90% of the total tropospheric delay (approximately 2.3 m in the zenith direction), but it stabilizes in time.The liquid component only becomes significant in a saturated atmosphere and is expected to be considerably small, measuring less than 0.1 mm/km to 0.4 mm/km under usual weather conditions for the C-band SAR system [9].We focused only on the hydrostatic and wet delays, disregarding the effect of the water content in clouds [7].The total zenith one-path tropospheric delay can be derived by the integration of air refractivity between the ground z and the top of the troposphere top z (typically 15 km), above which the delay is assumed to be nearly negligible.The zenith hydrostatic and wet delays were projected to the line of sight (LOS) direction, which is expressed as follows [36]:  the incidence angle of the ray, P the total atmospheric pressure, e the partial pressure of water vapor in hPa, and T the temperature in Kelvin.
-1 -1 d 287.05J kg K R  and indicate the dry air and water vapor specific gas constants [44].The coefficients , and are empirical constants [45].Geopotential height, temperature, and relative humidity were extracted from ERA-I reanalysis at each SAR acquisition.
We then converted the geopotential height to a regular vertical metric grid by mapping 0 g (

Tropospheric Delay in InSAR Data
A tropospheric delay can be split into hydrostatic, wet, and liquid components depending on variations in atmospheric conditions.The wet component is the major limiting factor due to water vapor instability in space and time (approximately 20 mm to 300 mm) [43].Unlike the wet delay, temperature and pressure are smooth in space, leading to a better-resolved, longer wavelength hydrostatic delay, accounting for about 90% of the total tropospheric delay (approximately 2.3 m in the zenith direction), but it stabilizes in time.The liquid component only becomes significant in a saturated atmosphere and is expected to be considerably small, measuring less than 0.1 mm/km to 0.4 mm/km under usual weather conditions for the C-band SAR system [9].We focused only on the hydrostatic and wet delays, disregarding the effect of the water content in clouds [7].The total zenith one-path tropospheric delay can be derived by the integration of air refractivity between the ground z and the top of the troposphere z top (typically 15 km), above which the delay is assumed to be nearly negligible.The zenith hydrostatic and wet delays were projected to the line of sight (LOS) direction, which is expressed as follows [36]: ]dz (1) where δL total LOS (z) is the summation of the hydrostatic delay δL hyd LOS (z) and the wet delay δL wet LOS (z), θ the incidence angle of the ray, P the total atmospheric pressure, e the partial pressure of water vapor in hPa, and T the temperature in Kelvin.R d = 287.05J kg −1 K −1 and R v = 461.495J kg −1 K −1 indicate the dry air and water vapor specific gas constants [44].The coefficients k 1 = 0.776 K Pa −1 ,k 2 = 0.716 K Pa −1 , and k 3 = 3.75 × 10 3 K 2 Pa −1 are empirical constants [45].Geopotential height, temperature, and relative humidity were extracted from ERA-I reanalysis at each SAR acquisition.We then converted the geopotential height to a regular vertical metric grid by mapping g 0 (g 0 ≈ 9.8 ms −2 ) with latitude.The partial pressure e was estimated from the relative humidity by a mixed Clausius-Clapeyron law as shown in the following equations [21].
where e * (z) represents the partial pressure of saturation water vapor, R e is the relative humidity, e * w (z) Providing meteorological parameters such as temperature, pressure, and water vapor partial pressure, we can use Equations ( 1)-( 3) to calculate the relative tropospheric delay ∆L t m ,t s LOS (z) for interferograms generated from SAR acquisitions acquired at time t m and t s (Equation ( 6)).Thus, the relative tropospheric phase delay ∆φ t m ,t s LOS for the interferograms is easily derived by differentiating delays at acquisition times t m and t s , as shown in Equations ( 6) and (7).
In California, we used ERA-I outputs acquired at 18:00 UTC, which has a negligible difference with the SAR acquisition (18:01 UTC).In Tianshan, we selected ERA-I outputs at 00:00 UTC and 06:00 UTC, which are closest to the SAR acquisition (05:12 UTC).The linear interpolation in time was applied to obtain the corresponding delay at the SAR acquisition.By using Equations ( 1)-( 5), we compute the total tropospheric delay on each ERA-I grid point within the blue rectangle in Figure 2. The resulting vertical profiles of delay are vertically interpolated using a spline interpolant and horizontally interpolated to each pixel in the SAR scene using a bilinear interpolant.Finally, we obtain the relative tropospheric delay for the interferogram by differentiating delay maps at each SAR acquisition (Equations ( 6) and ( 7)).Time and space interpolations, which can introduce uncertainties into delay corrections, are needed because of the low spatiotemporal resolution of the ERA-I model.However, our proposed PLE method is a phase-based method that estimates atmospheric contributions from the interferometric phase itself.Thus, an interpolation such as that in the ERA-I method is not required, avoiding the uncertainties introduced by poor interpolation.The details of the derivational processes and our corresponding strategies for data processing are discussed in the following subsections.

Power Law Model Based on ERA-I Reanalysis (PLE)
The relationship between the absolute tropospheric delay and elevation can be illustrated by a power law model, as shown in Equation ( 8) [35].Similar to Bekaert et al. [35,36], we estimate h c and α from the tropospheric delay curves.α is the slope of the mean delay curve in the log-log domain and h c is the constrained height where the standard deviation and average value of the relative delays are both smaller than 1.0 rad (~0.45 cm). Figure 3 shows the vertical profile of the tropospheric delay in the California test region estimated using ERA-I reanalysis on six grid points located in the green where K is a scaled factor that relates tropospheric delay to topography, α is the coefficient describing the power law model estimated by ERA-I reanalysis, and φ c is the phase delay at h c , indicating the relative phase-delay offset on an individual interferogram.
where K is a scaled factor that relates tropospheric delay to topography,  is the coefficient  Figure 3c shows that at a specific altitude (~5 km) the relative delays converged but the delays at master and slave acquisitions did not.The difference in  was found to be considerably small and the average value (~1.2) was often used when using the power law method.However, in many cases, the average  could not indicate the relationship between the phase and topography because the interferometric phase is the phase difference between the master and slave acquisitions.Thus, tropospheric delay depends on the change in refractivity rather than the absolute refractivity in each acquisition.The differences in  among different interferograms was found to be as high as 0.2 in our case, leading to a relative tropospheric delay change in 10 s of rad.To address this problem, we directly estimated  on each interferogram from the ERA-I reanalysis.Similar to the absolute tropospheric delay at SAR acquisition, the relative tropospheric delay on the interferogram also decreases with elevation and does not change significantly and converged when a certain height c h is reached.Thus, the tropospheric delay above this constrained height is negligible.
After estimating the power law relationship, we set c h at a fixed value by controlling the offset c  . was also estimated by plotting the slope of the log-log relationship between phase and topography, which were the constants for a given interferogram.However, the scaled factor is not a constant for an entire interferogram and is not easy to estimate due to the spatially variable tropospheric delay.Interferometric phases are a superposition of multiple signals, including tropospheric delay, deformation, and residual orbit errors, in which wavelengths contain a multi-scale dependent spectrum and may have different sensitivities to one another.For example, tropospheric signals contain short wavelength-scale (few kilometers) components and are introduced by turbulent components in the troposphere.Long wavelength-scale (tens of kilometers) components are Figure 3c shows that at a specific altitude (~5 km) the relative delays converged but the delays at master and slave acquisitions did not.The difference in α was found to be considerably small and the average value (~1.2) was often used when using the power law method.However, in many cases, the average α could not indicate the relationship between the phase and topography because the interferometric phase is the phase difference between the master and slave acquisitions.Thus, tropospheric delay depends on the change in refractivity rather than the absolute refractivity in each acquisition.The differences in α among different interferograms was found to be as high as 0.2 in our case, leading to a relative tropospheric delay change in 10 s of rad.To address this problem, we directly estimated α on each interferogram from the ERA-I reanalysis.Similar to the absolute tropospheric delay at SAR acquisition, the relative tropospheric delay on the interferogram also decreases with elevation and does not change significantly and converged when a certain height h c is reached.Thus, the tropospheric delay above this constrained height is negligible.
After estimating the power law relationship, we set h c at a fixed value by controlling the offset φ c .
α was also estimated by plotting the slope of the log-log relationship between phase and topography, which were the constants for a given interferogram.However, the scaled factor is not a constant for an entire interferogram and is not easy to estimate due to the spatially variable tropospheric delay.Interferometric phases are a superposition of multiple signals, including tropospheric delay, deformation, and residual orbit errors, in which wavelengths contain a multi-scale dependent spectrum and may have different sensitivities to one another.For example, tropospheric signals contain short wavelength-scale (few kilometers) components and are introduced by turbulent components in the troposphere.Long wavelength-scale (tens of kilometers) components are introduced by hydrostatic and wet delays due to weather condition variations and are partly correlated with topography.In addition, longer wavelength-scale (hundreds of kilometers) components introduced by ground deformation such as tidal loading, tectonic slow slip, or earthquakes may be more sensitive to other sources of confounding noises.Therefore, the selection of the band filter is not always trivial, as it requires empirical information about multi-scale dependent spectrums.However, we took advantage of the fact that the tropospheric signal is present at all wavelength scales to estimate the scaled factor from a spatial frequency band that was relatively insensitive to the other signals [33].Similar to the original power law method [35,36], band filtering was applied to the phase and topography to select an appropriate band such that other phase contributions were negligible.Thus, we can modify Equation ( 8) into the following [35]: where ∆φ trop, f ilt is the tropospheric delay fitted by the band-filtered phase ∆φ f ilt and reference height (h c − h) α .h in Equation ( 9) is the elevation of pixels extracted from SRTM DEM and in Equation ( 8) is the elevation of the ERA-I outputs.h c is the constrained height.K l is the local scaled factor estimated in a moving window for an interferogram.K m is the spatially variable factor calculated by multi-weighted K l .The band filtering would not change the scaled factor estimates only if there are contaminated signals.∆φ trop is the final estimated tropospheric delay on interferograms.α was estimated based on ERA-I reanalysis.
To introduce band-filtering artifacts conveniently, we removed the initial shift from interferometric phases.This shift was a long wavelength component that was filtered out later.We estimated the shift from interferometric phases over multiple windows, which varied in space and was influenced by other phase contributions.Therefore, in many cases, the shift cannot exactly represent the tropospheric delay offset and cannot be estimated from a spatial band that is insensitive to other phase signals.
Finally, we used Equation ( 9) to estimate the local factor K l between the band-filtered phase and topography, where the band was selected such that the contribution from other signals was negligible.Accounting for spatially variable tropospheric delays, we applied a moving window to the interferogram.The window was moving from the bottom to the top in the azimuth direction and from the left to the right in the range direction (Figure 1a).Over each window, the local band-filtered phase-topography relationship was analyzed and K l was estimated.The local robust estimation of K l plays a critical role in increasing the accuracy of tropospheric delays.An expanded robust technique was used to estimate K l in each window, which can reduce the sensitivity to outliers (e.g., due to residual deformation signal or other measurement error sources) [46].The error equation can be constructed based on Equation (9) as follows: where v i , ∆φ trop, f ilt,i and φ l represent the random error, band-filtered phase, and a constant shift for the ith pixel on a local window, respectively.N represents the number of pixels on the local window.
The parameter estimator can be derived based on the adjustment principle of least squares, as follows: X = (A T PA) −1 A T PL (12) where P represents the equivalent weight matrix of the vector L. The ith diagonal element of P is calculated by p i = p i w ii .w ii is the corresponding weight factor calculated by the function of IGGIII [47].
The parameters X can be robustly and reliably estimated by an expanded robust estimation [46] given that the initial value of the equivalent weight matrix was set.When the parameters are estimated, the covariance matrix of the estimates can be obtained through variance propagation in Equation ( 12) as follows: where σ2 0 is the corresponding variance factor of the unit weight.When the K l of each window was estimated, a multiple weighted and extrapolated method was applied to all pixels in the interferograms.The multi-weighted function was constructed with the STD of the K l estimate and the distance from the window centers to each pixel (Equations ( 14)-( 16)).
where S represents the reformed weights of the STD from the robust estimation (Equations ( 13)).G represents the Gaussian distribution based on the distance.s i is the ith STD of K l in the ith window.
The multiple weights W are a N × n matrix.K m is the final derived scaled factor for all pixels in Equation (10).N is the number of pixels, and n is the number of windows for the entire interferogram.
Figure 4 shows the flow chart of the PLE method on the interferogram.The PLE method can account for lateral variations in space by estimating K m in multiple windows and can be applied in regions without removing a preliminary deformation [33].The robust estimation technique is always available in each study area for a given day at a closet SAR acquisition.However, the selection of the band filter is not always trivial because it requires empirical information on multi-scale dependent spectrums.A potential option is to perform a statistical analysis of multiple spatial bands to select the appropriate bands.We compared the mitigation results of multiple bands with the corresponding unwrapped interferograms, expecting the largest reduction in STD if this band exhibits the best performance.The STD for each interferogram before and after correction using the PLE method by different spatial band filters is shown in Figure 5.The STD significantly varies between different bands in California (mean STD of 8.6 mm) than in Tianshan (mean STD of 7.1 mm).The arrows in Figure 5, which were selected for each interferogram, represent the bands with the largest reduction in STD.The spatial band-filtering approach relies on the manifestation of tropospheric delay at all spatial scales; often, the deformation signal may be at specific spatial wavelengths.However, determining the spatial bands, which are the representatives of the tropospheric delay, can be challenging.Apart from the deformation signal, turbulent and coherent tropospheric features, DEM errors, and orbital errors, which contaminate the spatial bands, can also be observed.Therefore, different deformation features, atmospheric conditions, possible DEM errors and orbit errors for different interferograms, as well as different topographic variations for different study regions lead to the different spatial bands, as selected in Figure 5.
The arrows in Figure 5, which were selected for each interferogram, represent the bands with the largest reduction in STD.The spatial band-filtering approach relies on the manifestation of tropospheric delay at all spatial scales; often, the deformation signal may be at specific spatial wavelengths.However, determining the spatial bands, which are the representatives of the tropospheric delay, can be challenging.Apart from the deformation signal, turbulent and coherent tropospheric features, DEM errors, and orbital errors, which contaminate the spatial bands, can also be observed.Therefore, different deformation features, atmospheric conditions, possible DEM errors and orbit errors for different interferograms, as well as different topographic variations for different study regions lead to the different spatial bands, as selected in Figure 5.  14)-( 16)); and (7) calculating the final tropospheric delay using Equation (10).(5) expanded robust estimation for local K l (9)) and D X (Equation ( 13)) in each window; (6) local K l values are extrapolated to each pixel using the multi-weighted method (Equations ( 14)-( 16)); and (7) calculating the final tropospheric delay using Equation (10).
The arrows in Figure 5, which were selected for each interferogram, represent the bands with the largest reduction in STD.The spatial band-filtering approach relies on the manifestation of tropospheric delay at all spatial scales; often, the deformation signal may be at specific spatial wavelengths.However, determining the spatial bands, which are the representatives of the tropospheric delay, can be challenging.Apart from the deformation signal, turbulent and coherent tropospheric features, DEM errors, and orbital errors, which contaminate the spatial bands, can also be observed.Therefore, different deformation features, atmospheric conditions, possible DEM errors and orbit errors for different interferograms, as well as different topographic variations for different study regions lead to the different spatial bands, as selected in Figure 5.

Results and Discussion
We applied our proposed PLE method to the Southern California and Tianshan study areas and compared tropospheric corrections estimated from the linear and ERA-I Accuracy was evaluated by using the unwrapped SAR interferograms and validated by comparing the delay measurements from the passive multispectral imager MERIS onboard the ENVISAT satellite.

Estimation of PLE Parameters
The study by Bekaert et al. [35] assumed that  and c h are constants for the study region and for each interferogram and estimated them by using sounding data.In our study, we assumed that the coefficients change with time (i.e., they are different for each interferogram), and instead of sounding data, we used ERA-I reanalysis to estimate these two coefficients, as shown in Figure 6.The relative delay curves can be fit by a power law function, and the coefficients vary with different interferograms.The success of the PLE method is limited to the cases in which the relative delays can be reasonably approximated by a power law function.

Results and Discussion
We applied our proposed PLE method to the Southern California and Tianshan study areas and compared tropospheric corrections estimated from the linear and ERA-I methods.Accuracy was evaluated by using the unwrapped SAR interferograms and validated by comparing the delay measurements from the passive multispectral imager MERIS onboard the ENVISAT satellite.

Estimation of PLE Parameters
The study by Bekaert et al. [35] assumed that α and h c are constants for the study region and for each interferogram and estimated them by using sounding data.In our study, we assumed that the coefficients change with time (i.e., they are different for each interferogram), and instead of sounding data, we used ERA-I reanalysis to estimate these two coefficients, as shown in Figure 6.The relative delay curves can be fit by a power law function, and the coefficients vary with different interferograms.The success of the PLE method is limited to the cases in which the relative delays can be reasonably approximated by a power law function.The sounding data have extremely low spatiotemporal resolution (only provide data twice per day).ERA-I reanalysis provides data at a relatively high spatiotemporal resolution (at six-hour intervals daily).Therefore, we can estimate the two coefficients from the statistical results of the tropospheric delay difference based on spatial variation, the values of which are estimated for each specific SAR acquisition.We used ERA-I reanalysis to calculate the relative delay curves at six grid points (the red points inside the green rectangle in Figure 2a), and determine c h by constraining the STD and average value.Then, we used the mean relative curves to estimate  by fitting a power law function.As shown in Figure 7, we set c h from the statistical results of the combinations of relative delay curves (difference of delay at two different days).Relative tropospheric delays (hydrostatic delay + wet delay, shown in Figure 7a) and wet delays (Figure 7b) are computed from ERA-I reanalysis in California for 120 days (18:00 UTC) between June and October 2009.The data could empirically represent the delays for interferogram 20090627/20090905 (example in Figure 3).The relative delays did not change significantly at a certain height (red arrows) and converged (red dashed lines) at a specific height c h of 12.8 km (Figure 7a); wet delays converged at the height of 6.8 km (Figure 7b).The constrained height for total tropospheric delays was higher than that for wet delays.This case shows that hydrostatic delays affected the estimation of the constrained height and cannot be simply ignored in the interferograms.However, the relative delays have an obvious tendency to converge at a certain height (indicated by red arrows in Figure 7a); this height value is extremely close to the convergent height for the wet delays, as shown in Figure 7b.Therefore, we set the constrained height to ~5 km, which is lower than the temporal statistics of the tropospheric delays, but near the wet delays.In our experiment, the constrained heights, which are specifically estimated from ERA-I reanalysis by spatial statistics, are generally lower than the global results from temporal statistics.Both coefficients were estimated by using ERA-I.We set h c to the altitude for which the STD and average value of the relative delay offsets are both less than 1.0 rad (~0.45 cm).
The sounding data have extremely low spatiotemporal resolution (only provide data twice per day).ERA-I reanalysis provides data at a relatively high spatiotemporal resolution (at six-hour intervals daily).Therefore, we can estimate the two coefficients from the statistical results of the tropospheric delay difference based on spatial variation, the values of which are estimated for each specific SAR acquisition.We used ERA-I reanalysis to calculate the relative delay curves at six grid points (the red points inside the green rectangle in Figure 2a), and determine h c by constraining the STD and average value.Then, we used the mean relative curves to estimate α by fitting a power law function.As shown in Figure 7, we set h c from the statistical results of the combinations of relative delay curves (difference of delay at two different days).Relative tropospheric delays (hydrostatic delay + wet delay, shown in Figure 7a) and wet delays (Figure 7b) are computed from ERA-I reanalysis in California for 120 days (18:00 UTC) between June and October 2009.The data could empirically represent the delays for interferogram 20090627/20090905 (example in Figure 3).The relative delays did not change significantly at a certain height (red arrows) and converged (red dashed lines) at a specific height h c of 12.8 km (Figure 7a); wet delays converged at the height of 6.8 km (Figure 7b).The constrained height for total tropospheric delays was higher than that for wet delays.This case shows that hydrostatic delays affected the estimation of the constrained height and cannot be simply ignored in the interferograms.However, the relative delays have an obvious tendency to converge at a certain height (indicated by red arrows in Figure 7a); this height value is extremely close to the convergent height for the wet delays, as shown in Figure 7b.Therefore, we set the constrained height to ~5 km, which is lower than the temporal statistics of the tropospheric delays, but near the wet delays.In our experiment, the constrained heights, which are specifically estimated from ERA-I reanalysis by spatial statistics, are generally lower than the global results from temporal statistics.To verify the effectiveness of our robust estimation method, we tested the approach on InSAR data from interferogram 20081025/20081129 in California.The dataset is shown in Table 1.The interferometric phase and the reference height are both band filtered in the 2-32 km spatial The local robust estimation of the scaled factor K l in each window plays a critical role in increasing the accuracy of tropospheric delays.As described in Section 2.3, we used a robust technique to estimate K l to avoid the impact from outliers.We can account for the spatial variation of tropospheric delay by applying the PLE method.The spatial variation maps of K m are shown in Figure 8.The PLE method splits the study region into 16 blocks with an overlap area of 50% in our study.The partially overlapping windows can reduce edge effects and ensure adjacent consistency.To verify the effectiveness of our robust estimation method, we tested the approach on InSAR data from interferogram 20081025/20081129 in California.The dataset is shown in Table 1.The interferometric phase and the reference height are both band filtered in the 2-32 km spatial To verify the effectiveness of our robust estimation method, we tested the approach on InSAR data from interferogram 20081025/20081129 in California.The dataset is shown in Table 1.The interferometric phase and the reference height are both band filtered in the 2-32 km spatial wavelength.We calculated the STD of the scaled factor estimates by assuming a linear model using the least squares and robust estimation methods.Although the dominant signals in the original band-filtered phases are contributed by topographically-correlated tropospheric delays, we can still detect some outlier rejection phases (72, indicated by red points in Figure 9a).The improvement in the STD values is ~10%.As shown in Figure 9b, we applied these two methods to the original band-filtered data with extra simulated outliers (200), which are randomly introduced to the locations.We detected 248 outlier rejection phases (recognition rate is ~91%, indicated by red points in Figure 9b), which do not contribute in estimating the scaled factor.The improvement of STD is ~15%.Figure 9c shows the relationship between improvement of STD and pollution rate of the data.Before a specific outlier percentage (~22%, indicated by the arrow in Figure 9c), the improvement of STD consistently increased with increasing outlier percentage.Then, when the outlier percentage is high, the improvement of STD decreases.When the outlier percentage is close to ~50%, the least squares regression and robust estimation would be ineffective.These results suggest that the robust estimation method yields the scaled factor robustly and reliably.
Remote Sens. 2017, 9, 765 15 of 23 wavelength.We calculated the STD of the scaled factor estimates by assuming a linear model using the least squares and robust estimation methods.Although the dominant signals in the original bandfiltered phases are contributed by topographically-correlated tropospheric delays, we can still detect some outlier rejection phases (72, indicated by red points in Figure 9a).The improvement in the STD values is ~10%.As shown in Figure 9b, we applied these two methods to the original band-filtered data with extra simulated outliers (200), which are randomly introduced to the locations.We detected 248 outlier rejection phases (recognition rate is ~91%, indicated by red points in Figure 9b), which do not contribute in estimating the scaled factor.The improvement of STD is ~15%.Figure 9c shows the relationship between improvement of STD and pollution rate of the data.Before a specific outlier percentage (~22%, indicated by the arrow in Figure 9c), the improvement of STD consistently increased with increasing outlier percentage.Then, when the outlier percentage is high, the improvement of STD decreases.When the outlier percentage is close to ~50%, the least squares regression and robust estimation would be ineffective.These results suggest that the robust estimation method yields the scaled factor robustly and reliably.

Estimation of Tropospheric Delays
We evaluated the performance of our PLE method through a comparison with the conventional linear and ERA-I methods.In addition, the accuracy was validated with the MERIS in California [16].We implemented the linear, ERA-I, and MERIS corrections, which are included in TRAIN.No MERIS data are available for comparison in Tianshan.MERIS-derived delay prediction is only used when

Estimation of Tropospheric Delays
We evaluated the performance of our PLE method through a comparison with the conventional linear and ERA-I methods.In addition, the accuracy was validated with the MERIS in California [16].We implemented the linear, ERA-I, and MERIS corrections, which are included in TRAIN.No MERIS data are available for comparison in Tianshan.MERIS-derived delay prediction is only used when the cloud coverage is less than 20% of SAR scene.The map of tropospheric wet delay was calculated from the MERIS precipitable water vapor (PWV) as follows: where Π is conversion factor given by Bevis et al. [48] as where ρ w is the density of liquid water, w is the ratio of molecular masses of water vapor and dry air (w ≈ 0.668), and T m is the weighted average of the temperature between the ground and reference height, as follows: We evaluated the T m at each pixel of interferogram using the atmospheric outputs from ERA-I to produce a map of Π.We then calculated the tropospheric wet delay by multiplying the MERIS PWV by Π, as shown in Equation (17).Values of Π that typically range from 5 to 7 with an average value of 6.2 are usually adopted [16].To validate our approach, we also added the hydrostatic delay derived from ERA-I to the MERIS wet delay to determine the total tropospheric delay.
Figure 10 shows examples of the tropospheric delays produced from the four methods for each test region.The interferograms are mainly caused by tropospheric delays, which are topographically correlated.The tropospheric delay is the sum of two components, namely, stratified delay and turbulent delay.The stratified delay is the result of different vertical refractivity profiles from the two SAR images; it mainly affects mountainous terrain and shows a strong correlation with topography.The turbulent delay results from turbulent processes in the atmosphere and affects flat terrain as well as mountainous terrain.The spatiotemporal filtering in InSAR time series methods cannot estimate the tropospheric delay because the seasonality of delay is underestimated [21].Therefore, we used the PLE method to estimate the stratified delay for each individual interferogram.The PLE method, however, similar to all phase-based methods, is not suitable in cases where significant turbulence mixing occurs.Therefore, we only focus on estimating the topographically correlated tropospheric delay that often dominates the tropospheric signal on an interferogram.The linear method (∆φ trop = Kh + φ c ) assumes a simple phase-topography relationship for the entire interferogram; therefore, it does not allow for spatially variable tropospheric delay.
Remote Sens. 2017, 9, 765 16 of 23 the cloud coverage is less than 20% of SAR scene.The map of tropospheric wet delay was calculated from the MERIS precipitable water vapor (PWV) as follows: where  is conversion factor given by Bevis et al. [48] as where  w is the density of liquid water, w is the ratio of molecular masses of water vapor and dry air ( 0.668  w ), and m T is the weighted average of the temperature between the ground and reference height, as follows: We evaluated the m T at each pixel of interferogram using the atmospheric outputs from ERA-I to produce a map of  .We then calculated the tropospheric wet delay by multiplying the MERIS PWV by  , as shown in Equation (17).Values of  that typically range from 5 to 7 with an average value of 6.2 are usually adopted [16].To validate our approach, we also added the hydrostatic delay derived from ERA-I to the MERIS wet delay to determine the total tropospheric delay.
Figure 10 shows examples of the tropospheric delays produced from the four methods for each test region.The interferograms are mainly caused by tropospheric delays, which are topographically correlated.The tropospheric delay is the sum of two components, namely, stratified delay and turbulent delay.The stratified delay is the result of different vertical refractivity profiles from the two SAR images; it mainly affects mountainous terrain and shows a strong correlation with topography.The turbulent delay results from turbulent processes in the atmosphere and affects flat terrain as well as mountainous terrain.The spatiotemporal filtering in InSAR time series methods cannot estimate the tropospheric delay because the seasonality of delay is underestimated [21].Therefore, we used the PLE method to estimate the stratified delay for each individual interferogram.The PLE method, however, similar to all phase-based methods, is not suitable in cases where significant turbulence mixing occurs.Therefore, we only focus on estimating the topographically correlated tropospheric delay that often dominates the tropospheric signal on an interferogram.The linear method (  Figure 10 shows that all these methods indicate delays with similar spatial distribution and magnitude.Discrepancies can also be observed between the four estimated tropospheric delays in California and Tianshan.MERIS-estimated tropospheric delays show a strong resemblance to the unwrapped interferograms because they were acquired simultaneously with ASAR images with high precision.The ERA-I method shows delays with an accurate location, but not always with a similar order of magnitude (e.g., interferogram 20100405/20100614 in Tianshan).This condition is partly due to an inaccurate interpolation of the ERA-I reanalysis outputs.The phase-based methods (the linear and PLE) assume a phase-topography relationship; thus, they are incapable of capturing turbulent signals.Although the ERA-I method can capture the spatial trends of tropospheric delay accurately in a region, this method cannot capture the higher gradients in the turbulent delay (e.g., interferogram 20080712/20080816 in California).MERIS delays, shown in Figure 10a, also contain the hydrostatic component estimated from ERA-I, allowing for comparison with the phase-based methods.The hydrostatic delay component is related to the ratio of atmospheric pressure P and the surface temperature T at the scale of an interferogram (typically 100 km).The spatial variations of pressure are usually small, typically at an order of magnitude of 1 hPa, whereas larger variations in water vapor partial pressure are common.Consequently, the differential wet delay usually overwhelms the differential hydrostatic delay on the interferogram.Therefore, most correction methods focused on estimating the wet delay.Doin et al. [49] validated the seasonal change in hydrostatic delay with temperature in seasonal fluctuations and pointed out that hydrostatic delay must be considered when the change in the local temperature is more than 10° within one year.Figure 11 shows the relation between height and hydrostatic delay estimated using ERA-I method and the magnitude of hydrostatic delays can reach 1 rad.This condition demonstrates that the hydrostatic delay cannot be simply ignored on the interferograms.Figure 10 shows that all these methods indicate delays with similar spatial distribution and magnitude.Discrepancies can also be observed between the four estimated tropospheric delays in California and Tianshan.MERIS-estimated tropospheric delays show a strong resemblance to the unwrapped interferograms because they were acquired simultaneously with ASAR images with high precision.The ERA-I method shows delays with an accurate location, but not always with a similar order of magnitude (e.g., interferogram 20100405/20100614 in Tianshan).This condition is partly due to an inaccurate interpolation of the ERA-I reanalysis outputs.The phase-based methods (the linear and PLE) assume a phase-topography relationship; thus, they are incapable of capturing turbulent signals.Although the ERA-I method can capture the spatial trends of tropospheric delay accurately in a region, this method cannot capture the higher gradients in the turbulent delay (e.g., interferogram 20080712/20080816 in California).MERIS delays, shown in Figure 10a, also contain the hydrostatic component estimated from ERA-I, allowing for comparison with the phase-based methods.The hydrostatic delay component is related to the ratio of atmospheric pressure P and the surface temperature T at the scale of an interferogram (typically 100 km).The spatial variations of pressure are usually small, typically at an order of magnitude of 1 hPa, whereas larger variations in water vapor partial pressure are common.Consequently, the differential wet delay usually overwhelms the differential hydrostatic delay on the interferogram.Therefore, most correction methods focused on estimating the wet delay.Doin et al. [49] validated the seasonal change in hydrostatic delay with temperature in seasonal fluctuations and pointed out that hydrostatic delay must be considered when the change in the local temperature is more than 10 • within one year.Figure 11 shows the relation between height and hydrostatic delay estimated using ERA-I method and the magnitude of hydrostatic delays can reach 1 rad.This condition demonstrates that the hydrostatic delay cannot be simply ignored on the interferograms.Figure 10 shows that all these methods indicate delays with similar spatial distribution and magnitude.Discrepancies can also be observed between the four estimated tropospheric delays in California and Tianshan.MERIS-estimated tropospheric delays show a strong resemblance to the unwrapped interferograms because they were acquired simultaneously with ASAR images with high precision.The ERA-I method shows delays with an accurate location, but not always with a similar order of magnitude (e.g., interferogram 20100405/20100614 in Tianshan).This condition is partly due to an inaccurate interpolation of the ERA-I reanalysis outputs.The phase-based methods (the linear and PLE) assume a phase-topography relationship; thus, they are incapable of capturing turbulent signals.Although the ERA-I method can capture the spatial trends of tropospheric delay accurately in a region, this method cannot capture the higher gradients in the turbulent delay (e.g., interferogram 20080712/20080816 in California).MERIS delays, shown in Figure 10a, also contain the hydrostatic component estimated from ERA-I, allowing for comparison with the phase-based methods.The hydrostatic delay component is related to the ratio of atmospheric pressure P and the surface temperature T at the scale of an interferogram (typically 100 km).The spatial variations of pressure are usually small, typically at an order of magnitude of 1 hPa, whereas larger variations in water vapor partial pressure are common.Consequently, the differential wet delay usually overwhelms the differential hydrostatic delay on the interferogram.Therefore, most correction methods focused on estimating the wet delay.Doin et al. [49] validated the seasonal change in hydrostatic delay with temperature in seasonal fluctuations and pointed out that hydrostatic delay must be considered when the change in the local temperature is more than 10° within one year.Figure 11 shows the relation between height and hydrostatic delay estimated using ERA-I method and the magnitude of hydrostatic delays can reach 1 rad.This condition demonstrates that the hydrostatic delay cannot be simply ignored on the interferograms.

Statistical Comparison of Tropospheric Correction Methods
We performed a statistical analysis to compare tropospheric delays estimated from the four correction methods (MERIS, the linear, ERA-I, and PLE methods) to the unwrapped interferograms.In Figure 12, we performed the comparison by plotting the interferometric phase before and after correction against the elevation for interferograms presented in Figure 10. Figure 12a shows that all these methods were effective in mitigating topographically correlated tropospheric delays in California.The slope between the corrected phases and elevation reduced in all correction methods.However, limited improvements were found in Tianshan, as shown in Figure 12b because the tropospheric delays do not show a strong relationship with topography, and might have mainly resulted from the turbulent delay.After using ERA-I for corrections, the results even deteriorated due to the low spatial resolution of the reanalysis, indicating that ERA-I could introduce atmospheric estimation error into the final displacement results because of the uncertainties produced by the interpolation operation.However, the PLE method always yielded effective delay mitigation, as shown by STD reduction after tropospheric correction using the three different tropospheric techniques.We expect a high STD reduction when a method performs efficiently.As shown in Figure 13, the left panel shows the STD of the unwrapped interferograms, whereas the right panel shows the residual STD after correction using the different tropospheric estimation methods.

Statistical Comparison of Tropospheric Correction Methods
We performed a statistical analysis to compare tropospheric delays estimated from the four correction methods (MERIS, the linear, ERA-I, and PLE methods) to the unwrapped interferograms.In Figure 12, we performed the comparison by plotting the interferometric phase before and after correction against the elevation for interferograms presented in Figure 10. Figure 12a shows that all these methods were effective in mitigating topographically correlated tropospheric delays in California.The slope between the corrected phases and elevation reduced in all correction methods.However, limited improvements were found in Tianshan, as shown in Figure 12b because the tropospheric delays do not show a strong relationship with topography, and might have mainly resulted from the turbulent delay.After using ERA-I for corrections, the results even deteriorated due to the low spatial resolution of the reanalysis, indicating that ERA-I could introduce atmospheric estimation error into the final displacement results because of the uncertainties produced by the interpolation operation.However, the PLE method always yielded effective delay mitigation, as shown by STD reduction after tropospheric correction using the three different tropospheric techniques.We expect a high STD reduction when a method performs efficiently.As shown in Figure 13, the left panel shows the STD of the unwrapped interferograms, whereas the right panel shows the residual STD after correction using the different tropospheric estimation methods.In California, all methods are effective in mitigating tropospheric delays, leading to STD reduction.MERIS, on average, provides better STD reduction (~41%) than the linear method (~33%) and ERA-I (~23%).The MERIS results in STD reduction are approximately equal to the MERIS ~1.8 rad accuracy level computed using the measured ~1.1 mm PWV accuracy (brown symbols in Figure 13a) [16].Moreover, MERIS delays containing hydrostatic component from ERA-I showed better reduction (average value was ~2%) than those without hydrostatic component.The average STD reduction of the residuals after correction using MERIS without and with hydrostatic delay is ~39% In California, all methods are effective in mitigating tropospheric delays, leading to STD reduction.MERIS, on average, provides better STD reduction (~41%) than the linear method (~33%) and ERA-I (~23%).The MERIS results in STD reduction are approximately equal to the MERIS ~1.8 rad accuracy level computed using the measured ~1.1 mm PWV accuracy (brown symbols in Figure 13a) [16].Moreover, MERIS delays containing hydrostatic component from ERA-I showed better reduction (average value was ~2%) than those without hydrostatic component.The average STD reduction of the residuals after correction using MERIS without and with hydrostatic delay is ~39% and ~41%, respectively.Therefore, we cannot simply ignore the hydrostatic delay effect.Unlike ERA-I, no uncertainty is introduced through temporal interpolation for MERIS because MERIS measurements were acquired simultaneously with ASAR images.Figure 13a shows that ERA-I yielded an average STD reduction of ~23% (blue symbols), which was the worst among the four methods.This result was probably due to the low spatial resolution.The phase-based methods, mainly designed to resolve the topographically correlated delay, can retrieve the tropospheric delay; the PLE and the linear methods provide an STD reduction of ~42% and ~33%, respectively.Thus, the PLE method performed better than the linear method because the PLE method can account for the spatial variation of tropospheric delay, whereas the conventional linear model cannot.The tropospheric delays in the region of California are mainly the stratified delays that are strongly correlated with topography (also shown in Figure 12a).Thus, the phase-based methods performed well in this area.The PLE method given that it is applied locally is more subjected to contamination from various tropospheric delay components present at different spatial scales than the global linear method.For example, for the interferogram 20081025/20081129, the STD reduction is ~37% after the linear correction, whereas the reduction is ~29% after the PLE method.Considering that both phase-based methods have the potential to introduce inaccurate signals in the presence of serious turbulence and coherent short-scale tropospheric signals, neither of them should be applied in those cases.
In Tianshan, MERIS-derived tropospheric delays were not available for comparison.We did not observe an improvement when using ERA-I except for interferogram 20090209/20090420 (blue symbols in Figure 13b).On average, STD increased (increase ~57%) rather than decreased after correction using ERA-I.The two possible reasons for this deterioration are as follows.First, the requirement for temporal, lateral, and vertical interpolations of the atmospheric parameters introduces uncertainties into the delay estimations.Second, the turbulent process in the atmosphere in Tianshan often appears due to the high altitude and complicated terrain.Parts of the tropospheric delays on the interferograms are mainly the small-scale turbulent delay, which cannot be captured by ERA-I.For phase-based methods, we found slight improvements, with an average reduction in STD of ~10% and ~13% for the linear and PLE methods, respectively.
Remote Sens. 2017, 9, 765 19 of 23 and ~41%, respectively.Therefore, we cannot simply ignore the hydrostatic delay effect.Unlike ERA-I, no uncertainty is introduced through temporal interpolation for MERIS because MERIS measurements were acquired simultaneously with ASAR images.Figure 13a shows that ERA-I yielded an average STD reduction of ~23% (blue symbols), which was the worst among the four methods.This result was probably due to the low spatial resolution.The phase-based methods, mainly designed to resolve the topographically correlated delay, can retrieve the tropospheric delay; the PLE and the linear methods provide an STD reduction of ~42% and ~33%, respectively.Thus, the PLE method performed better than the linear method because the PLE method can account for the spatial variation of tropospheric delay, whereas the conventional linear model cannot.The tropospheric delays in the region of California are mainly the stratified delays that are strongly correlated with topography (also shown in Figure 12a).Thus, the phase-based methods performed well in this area.The PLE method given that it is applied locally is more subjected to contamination from various tropospheric delay components present at different spatial scales than the global linear method.For example, for the interferogram 20081025/20081129, the STD reduction is ~37% after the linear correction, whereas the reduction is ~29% after the PLE method.Considering that both phasebased methods have the potential to introduce inaccurate signals in the presence of serious turbulence and coherent short-scale tropospheric signals, neither of them should be applied in those cases.
In Tianshan, MERIS-derived tropospheric delays were not available for comparison.We did not observe an improvement when using ERA-I except for interferogram 20090209/20090420 (blue symbols in Figure 13b).On average, STD increased (increase ~57%) rather than decreased after correction using ERA-I.The two possible reasons for this deterioration are as follows.First, the requirement for temporal, lateral, and vertical interpolations of the atmospheric parameters introduces uncertainties into the delay estimations.Second, the turbulent process in the atmosphere in Tianshan often appears due to the high altitude and complicated terrain.Parts of the tropospheric delays on the interferograms are mainly the small-scale turbulent delay, which cannot be captured by ERA-I.For phase-based methods, we found slight improvements, with an average reduction in STD of ~10% and ~13% for the linear and PLE methods, respectively.
(a) California  In conclusion, the PLE method outperforms the linear and ERA-I methods in the regions where tropospheric delays are correlated with topography, and agrees well with MERIS measurements.Moreover, the success of our PLE method is limited to the cases in which the relative delays can be reasonably approximated by a power law function.The main limitation of the robust estimation technique in the PLE method is that the algorithm should be implemented iteratively; thus, its computational efficiency is lower than that of least squares regression.From our analysis, each method has its limitations, which vary with the region of application and with time.In future work, we will identify suitable metrics and combine different tropospheric estimation methods, which could further improve the accuracy of InSAR processing.

Conclusions
We developed and tested a method to estimate the tropospheric delay in radar interferogram using a power law model based on ERA-I reanalysis.The PLE method is a modification of the power law method proposed by Bekaert et al., which mainly has two novelties: one is that the power law coefficients are estimated using ERA-I reanalysis instead of the sounding data, which is usually available in a study region; the other is that the scaled factor is estimated using the robust estimation technique instead of least squares regression, which is outlier-resistant.To separate deformation and tropospheric signals, a frequency band insensitive to deformation is tested and selected in the study.The PLE method splits a region into multiple blocks (approximately 100 km 2 in our study) and estimates the scaled factor robustly to capture the spatial variations of tropospheric delay in lager regions.We tested the PLE method on ENVISAT ASAR images in Southern California and Tianshan, which have different geographical regions and atmospheric conditions.In our study regions, the iterative calculation approach applied in the PLE method is relatively straightforward and efficient when estimating tropospheric delay.
We performed a comparison of the tropospheric corrections estimated from the conventional linear, independent ERA-I, and PLE methods.The accuracy was evaluated by corrected SAR interferograms and was validated by using MERIS measurements.After correction, our PLE method had greater reduction in the topographically correlated tropospheric signals than the conventional linear and ERA-I methods.Moreover, our method improved the agreement with delay measured from MERIS.The success of the PLE method is limited to cases in which the relative delays can be reasonably approximated by a power law relationship between topography and phase.The main limitation of the robust estimation technique in the PLE method is that the algorithm should be implemented iteratively; thus, its computational efficiency is lower than that of least squares regression.Each method has its limitations in different regions and at different times.Therefore, in In conclusion, the PLE method outperforms the linear and ERA-I methods in the regions where tropospheric delays are correlated with topography, and agrees well with MERIS measurements.Moreover, the success of our PLE method is limited to the cases in which the relative delays can be reasonably approximated by a power law function.The main limitation of the robust estimation technique in the PLE method is that the algorithm should be implemented iteratively; thus, its computational efficiency is lower than that of least squares regression.From our analysis, each method has its limitations, which vary with the region of application and with time.In future work, we will identify suitable metrics and combine different tropospheric estimation methods, which could further improve the accuracy of InSAR processing.

Conclusions
We developed and tested a method to estimate the tropospheric delay in radar interferogram using a power law model based on ERA-I reanalysis.The PLE method is a modification of the power law method proposed by Bekaert et al., which mainly has two novelties: one is that the power law coefficients are estimated using ERA-I reanalysis instead of the sounding data, which is usually available in a study region; the other is that the scaled factor is estimated using the robust estimation technique instead of least squares regression, which is outlier-resistant.To separate deformation and tropospheric signals, a frequency band insensitive to deformation is tested and selected in the study.The PLE method splits a region into multiple blocks (approximately 100 km 2 in our study) and estimates the scaled factor robustly to capture the spatial variations of tropospheric delay in lager regions.We tested the PLE method on ENVISAT ASAR images in Southern California and Tianshan, which have different geographical regions and atmospheric conditions.In our study regions, the iterative calculation approach applied in the PLE method is relatively straightforward and efficient when estimating tropospheric delay.
We performed a comparison of the tropospheric corrections estimated from the conventional linear, independent ERA-I, and PLE methods.The accuracy was evaluated by corrected SAR interferograms and was validated by using MERIS measurements.After correction, our PLE method had greater reduction in the topographically correlated tropospheric signals than the conventional linear and ERA-I methods.Moreover, our method improved the agreement with delay measured from MERIS.The success of the PLE method is limited to cases in which the relative delays can be reasonably approximated by a power law relationship between topography and phase.The main limitation of the robust estimation technique in the PLE method is that the algorithm should be implemented iteratively; thus, its computational efficiency is lower than that of least squares regression.Each method has its limitations in different regions and at different times.Therefore, in some cases, different tropospheric estimation methods should be combined due to the complicated geological structure of the earth and the variability in the atmospheric circulation.

Figure 1 .
Figure 1.InSAR data coverage for: (a) California; and (b) Tianshan.The background is the SRTM DEM.Black polygons indicate the illuminated ground area of the entire SAR image and red polygons represent our study regions.The black circular markers indicate the center of multiple windows used for the PLE method.The blue dashed and solid lines indicate the moving windows.Windows are arranged in an ascending manner and start in the lower left corner.These windows within a tropospheric region constrain the delay estimation by a 50% overlap.

Figure 1 .
Figure 1.InSAR data coverage for: (a) California; and (b) Tianshan.The background is the SRTM DEM.Black polygons indicate the illuminated ground area of the entire SAR image and red polygons represent our study regions.The black circular markers indicate the center of multiple windows used for the PLE method.The blue dashed and solid lines indicate the moving windows.Windows are arranged in an ascending manner and start in the lower left corner.These windows within a tropospheric region constrain the delay estimation by a 50% overlap.

Figure 2 .
Figure 2. The distribution of ERA-I grid points for: (a) California; and (b) Tianshan.Black circular markers indicate the location of the ERA-I grid points at ~75 km spatial resolution.The red points within the circular markers are the points used for the ERA-I method.Red solid lines outline our test regions.The background is the SRTM topography.The points within the green rectangles were used for our PLE method.

Figure 2 .
Figure 2. The distribution of ERA-I grid points for: (a) California; and (b) Tianshan.Black circular markers indicate the location of the ERA-I grid points at ~75 km spatial resolution.The red points within the circular markers are the points used for the ERA-I method.Red solid lines outline our test regions.The background is the SRTM topography.The points within the green rectangles were used for our PLE method.
rectangle in Figure 2a.The absolute tropospheric delay at SAR acquisitions 20090627 and 20090905 are shown in Figure 3a,b.The relative tropospheric delay (∆φ t m ,t s LOS ) on the interferogram formed by these two images is shown in Figure 3c.By plotting the mean relative tropospheric delay (black solid line) against the reference elevation [(h c − h) α ], we can determine the power law relationship (red dashed line) indicated by the linear pattern.
Figure 2a.The absolute tropospheric delay at SAR acquisitions 20090627 and 20090905 are shown in Figure 3a,b.The relative tropospheric delay ( two images is shown in Figure 3c.By plotting the mean relative tropospheric delay (black solid line) against the reference elevation [ ( ) c h h   ], we can determine the power law relationship (red dashed line) indicated by the linear pattern.
describing the power law model estimated by ERA-I reanalysis, and c  is the phase delay at c h , indicating the relative phase-delay offset on an individual interferogram.

Figure 3 .
Figure 3. Mean delay and fitted power law with respect to the elevation.The black solid line represents the corresponding phase delay calculated from the ERA-I data at six grid points (18:00 UTC) on interferogram 20090627/20090905 in California.The points are shown within the green rectangle in Figure 2a.The red dashed line shows the best fit line using the power law model.The black dashed line indicates the maximum height of the InSAR region.(a,b) The tropospheric delays for the master and slave acquisitions show a near power law function with  being 1.18 and 1.19.The delay offset at the height of 4.7 km was added back after a log-log linear estimation.(c) The relative tropospheric delay shows a power law relationship (   1.39).

Figure 3 .
Figure 3. Mean delay and fitted power law with respect to the elevation.The black solid line represents the corresponding phase delay calculated from the ERA-I data at six grid points (18:00 UTC) on interferogram 20090627/20090905 in California.The points are shown within the green rectangle in Figure 2a.The red dashed line shows the best fit line using the power law model.The black dashed line indicates the maximum height of the InSAR region.(a,b) The tropospheric delays for the master and slave acquisitions show a near power law function with α being 1.18 and 1.19.The delay offset at the height of 4.7 km was added back after a log-log linear estimation.(c) The relative tropospheric delay shows a power law relationship relationship (α = 1.39).

Figure 4 .
Figure 4. Flow chart of tropospheric delay estimation using the PLE method for an interferogram: (1) log-log fitting Equation (8) to obtain  and c h based on ERA-I reanalysis; (2) extracting elevation h for all pixels from SRTM DEM; (3) and (4) band filtering in the corresponding spatial wavelength;(5) expanded robust estimation for local

Figure 4 .
Figure 4. Flow chart of tropospheric delay estimation using the PLE method for an interferogram: (1) log-log fitting Equation (8) to obtain α and h c based on ERA-I reanalysis; (2) extracting elevation h for all pixels from SRTM DEM; (3) and (4) band filtering in the corresponding spatial wavelength;(5) expanded robust estimation for local K l (9)) and D X (Equation (13)) in each window;(6) local K l values are extrapolated to each pixel using the multi-weighted method (Equations (14)-(16)); and (7) calculating the final tropospheric delay using Equation(10).

Figure 4 .
Figure 4. Flow chart of tropospheric delay estimation using the PLE method for an interferogram: (1) log-log fitting Equation (8) to obtain  and c h based on ERA-I reanalysis; (2) extracting elevation h for all pixels from SRTM DEM; (3) and (4) band filtering in the corresponding spatial wavelength;(5) expanded robust estimation for local

Figure 5 .
Figure 5.Comparison of mitigation results of different spatial bands for each interferogram in: (a) California; and (b) Tianshan.The different color bars represent different interferograms.Upper panel represents STD of the unwrapped interferograms (hollow bars) compared to the residual STD after correcting the tropospheric delay derived from the PLE method (filled bars).Lower panel represents reduction of STD after delay correction of different spatial band filters.Arrows were selected for each interferogram to indicate the bands with the largest reduction in STD.

Figure 5 .
Figure 5.Comparison of mitigation results of different spatial bands for each interferogram in: (a) California; and (b) Tianshan.The different color bars represent different interferograms.Upper panel represents STD of the unwrapped interferograms (hollow bars) compared to the residual STD after correcting the tropospheric delay derived from the PLE method (filled bars).Lower panel represents reduction of STD after delay correction of different spatial band filters.Arrows were selected for each interferogram to indicate the bands with the largest reduction in STD.

Figure 6 .
Figure 6.Examples of constrained height c h and power law coefficient  for our study area in California.Both coefficients were estimated by using ERA-I.We set c h to the altitude for which the STD and average value of the relative delay offsets are both less than 1.0 rad (~0.45 cm).

Figure 6 .
Figure 6.Examples of constrained height h c and power law coefficient α for our study area in California.Both coefficients were estimated by using ERA-I.We set h c to the altitude for which the STD and average value of the relative delay offsets are both less than 1.0 rad (~0.45 cm).

Figure 7 .Figure 8 .
Figure 7. (a) Relative tropospheric delays computed from ERA-I reanalysis in California for 120 days (18:00 UTC) between June and October 2009.Relative delay curves are 800 combinations of the tropospheric delay difference at different days, covering the acquisition for interferogram 20090627/20090905 (example in Figure 3).The relative delays did not change significantly at certain height (red arrows) and converged (red lines) at a specific height c h of 12.8 km.The STD of the relative curves is less than 1 rad.(b) Wet delays converged at the height of 6.8 km.The delays are calculated by integrating the refractivity from surface upward (15 km) and projecting to LOS direction.The local robust estimation of the scaled factor l K in each window plays a critical role in increasing the accuracy of tropospheric delays.As described in Section 2.3, we used a robust technique to estimate l K to avoid the impact from outliers.We can account for the spatial variation of tropospheric delay by applying the PLE method.The spatial variation maps of m K are shown in Figure 8.The PLE method splits the study region into 16 blocks with an overlap area of 50% in our study.The partially overlapping windows can reduce edge effects and ensure adjacent consistency.

Figure 7 .
Figure 7. Relative tropospheric delays computed from ERA-I reanalysis in California for 120 days (18:00 UTC) between June and October 2009.Relative delay curves are 800 combinations of the tropospheric delay difference at different days, covering the acquisition for interferogram 20090627/20090905 (example in Figure 3).The relative delays did not change significantly at certain height (red arrows) and converged (red lines) at a specific height h c of 12.8 km.The STD of the relative curves is less than 1 rad.(b) Wet delays converged at the height of 6.8 km.The delays are calculated by integrating the refractivity from surface upward (15 km) and projecting to LOS direction.

Figure 7 .Figure 8 .
Figure 7. (a) Relative tropospheric delays computed from ERA-I reanalysis in California for 120 days (18:00 UTC) between June and October 2009.Relative delay curves are 800 combinations of the tropospheric delay difference at different days, covering the acquisition for interferogram 20090627/20090905 (example in Figure 3).The relative delays did not change significantly at certain height (red arrows) and converged (red lines) at a specific height c h of 12.8 km.The STD of the relative curves is less than 1 rad.(b) Wet delays converged at the height of 6.8 km.The delays are calculated by integrating the refractivity from surface upward (15 km) and projecting to LOS direction.The local robust estimation of the scaled factor l K in each window plays a critical role in increasing the accuracy of tropospheric delays.As described in Section 2.3, we used a robust technique to estimate l K to avoid the impact from outliers.We can account for the spatial variation of tropospheric delay by applying the PLE method.The spatial variation maps of m K are shown in Figure 8.The PLE method splits the study region into 16 blocks with an overlap area of 50% in our study.The partially overlapping windows can reduce edge effects and ensure adjacent consistency.

Figure 8 .
Figure 8. Examples of spatial maps of

Figure 9 .
Figure 9. Examples to test the effectiveness of robust estimation.(a) Outlier rejection results of robust estimation for the original band-filtered data obtained from interferogram 20081025/20081129 in California.(b) Results of robust estimation for the original band-filtered data with extra simulated outliers, which are randomly introduced to the locations.(c) Improvement of STD by robust estimation compared to least squares regression.The average STD improvement (red curve) is calculated from the datasets of 10 simulations.The arrow in (c) indicates the outlier percentage (~22%) with the largest improvement of STD.Gray points in (a) and (b) represent the band-filtered phases, which are considered as "clean observations" without containing outliers.Blue points are suspected to contain outliers (corresponding weights are reduced).Red points are considered to contain outliers (corresponding weights are set to zero), which do not contribute in estimating the scaled factor.

Figure 9 .
Figure 9. Examples to test the effectiveness of robust estimation.(a) Outlier rejection results of robust estimation for the original band-filtered data obtained from interferogram 20081025/20081129 in California.(b) Results of robust estimation for the original band-filtered data with extra simulated outliers, which are randomly introduced to the locations.(c) Improvement of STD by robust estimation compared to least squares regression.The average STD improvement (red curve) is calculated from the datasets of 10 simulations.The arrow in (c) indicates the outlier percentage (~22%) with the largest improvement of STD.Gray points in (a) and (b) represent the band-filtered phases, which are considered as "clean observations" without containing outliers.Blue points are suspected to contain outliers (corresponding weights are reduced).Red points are considered to contain outliers (corresponding weights are set to zero), which do not contribute in estimating the scaled factor.

Figure 10 .
Figure 10.Examples of tropospheric delay estimated by different methods in: (a) California; and (b) Tianshan.Columns represent the unwrapped interferograms, and the estimated tropospheric delays using MERIS, the linear, ERA-I, and PLE methods (from left to right).All observations are converted to displacements in LOS.

Figure 11 .
Figure 11.Relationship between height and hydrostatic delay estimated from ERA-I method in California.

Figure 10 .
Figure 10.Examples of tropospheric delay estimated by different methods in: (a) California; and (b) Tianshan.Columns represent the unwrapped interferograms, and the estimated tropospheric delays using MERIS, the linear, ERA-I, and PLE methods (from left to right).All observations are converted to displacements in LOS.

Figure 10 .
Figure 10.Examples of tropospheric delay estimated by different methods in: (a) California; and (b) Tianshan.Columns represent the unwrapped interferograms, and the estimated tropospheric delays using MERIS, the linear, ERA-I, and PLE methods (from left to right).All observations are converted to displacements in LOS.

Figure 11 .
Figure 11.Relationship between height and hydrostatic delay estimated from ERA-I method in California.

Figure 11 .
Figure 11.Relationship between height and hydrostatic delay estimated from ERA-I method in California.

Figure 12 .
Figure 12.Examples of final results after different tropospheric corrections in: (a) California; and (b) Tianshan.MERIS delays contain hydrostatic delay component from ERA-I.

Figure 12 .
Figure 12.Examples of final results after different tropospheric corrections in: (a) California; and (b) Tianshan.MERIS delays contain hydrostatic delay component from ERA-I.

Figure 13 .
Figure 13.Comparisons between different tropospheric correction methods in: (a) California; and (b) Tianshan.Left panel: STD of the unwrapped interferograms.Right panel: STD after correction of the tropospheric delay with different methods (colored symbols).Other illustrations show the STD reduction.Point, circle, cross, and square symbols correspond to different interferograms, and triangular symbols represent the average value.Wet delay components are indicated by (w) and hydrostatic delay components by (h).

Figure 13 .
Figure 13.Comparisons between different tropospheric correction methods in: (a) California; and (b) Tianshan.Left panel: STD of the unwrapped interferograms.Right panel: STD after correction of the tropospheric delay with different methods (colored symbols).Other illustrations show the STD reduction.Point, circle, cross, and square symbols correspond to different interferograms, and triangular symbols represent the average value.Wet delay components are indicated by (w) and hydrostatic delay components by (h).

California Date (yyyymmdd) Perpendicular Baseline (m) Temporal Baseline (Days) Doppler Central Baseline (Hz) Average Coherence
provides the Clausius-Clapeyron law "over water", and e * i (z) gives the Clausius-Clapeyron law "over ice".Over T 0 = 273.16K, e * w (z) is used, and under T i = 250.16K, e * i (z) is used and a weighted average of both is used in between.The parameter values are a 1 = 611.21hPa, a 3,w = 17.502, a 4,w = 32.19K, a 3,i = 22.587 K, and a 4,i = −0.7 K.