A spatially variable power law tropospheric correction technique for InSAR data

Microwave signals traveling through the troposphere are subject to delays. These delays are mainly described by spatial and temporal variations in pressure, temperature, and relative humidity in the lower part of the troposphere, resulting in a spatially varying tropospheric signal in interferometric synthetic aperture radar (InSAR). Tropospheric correction techniques rely either on external data, often limited by spatial and temporal accuracy or can be estimated from the high‐resolution interferometric phase itself. However, current phase‐estimated correction techniques do not account for the spatial variability of the tropospheric properties and fail to capture tropospheric signals over larger regions. Here we propose and test a novel power law correction method that accounts for spatial variability in atmospheric properties and can be applied to interferograms containing topographically correlated deformation. The power law model has its reference fixed at the relative top of the troposphere and describes, through a power law relationship, how the phase delay varies with altitude. We find the power law model reduces tropospheric signals both locally (on average by ∼0.45 cm for each kilometer of elevation in Mexico) and the long‐wavelength components, leading to an improved fit to independent Global Navigation Satellite Systems data. The power law model can be applied in presence of deformation, over a range of different time periods and in different atmospheric conditions, and thus permits the detection of smaller‐magnitude crustal deformation signals with InSAR.


Introduction
Radar signals propagating through the atmosphere are affected by the medium in which they travel, which can result in a phase delay or advance. This is a major limiting factor when using Interferometric Synthetic Aperture Radar (InSAR ) [e.g., Wright et al., 2001;Hooper et al., 2012;Cavalie et al., 2013]. A phase delay through the atmosphere can be characterized by the refractivity, N, [Smith and Weintraub, 1953;Thayer, 1974;Davis et al., 1985], which can be split into hydrostatic, wet, liquid, and ionospheric components. For the troposphere, variations of temperature, pressure, and relative humidity lead to a spatially varying tropospheric phase delay, which partly correlates with topography. We focus only on the hydrostatic and wet components of the refractivity, N, as [Hanssen, 2001] the liquid component is expected to be small and will only become significant in case of a saturated atmosphere, and the ionospheric delays are assumed to be small for C band, leading to where P indicates total atmospheric pressure, T the temperature, and e the partial pressure of water vapor. The coefficients k 1 , k ′ 2 , and k 3 are constants estimated in literature as k 1 = 77.6 K hPa −1 , k ′ 2 = 23.3 K hPa −1 , and k 3 = 3.75 ⋅ 10 5 K 2 hPa −1 [Smith and Weintraub, 1953;Thayer, 1974]. The tropospheric phase delay (two-way), tropo , at a specific height, h = h 1 , corresponds to the integration of the refractivity between h 1 and the top of the troposphere h top along the radar line of sight [Hanssen, 2001].
where d tropo is the tropospheric delay (one-way), the incidence angle, the radar wavelength, and −4 ∕ a conversion factor to convert from pseudorange increase to phase delay.
BEKAERT ET AL. ©2015. The Authors.  Empirically it can be seen that at a specific altitude h 0 the relative delays converge to zero. We set h 0 to the altitude for which the standard deviation of the relative delays is smaller than 0.5 cm. Sounding measurements made during ascents of large inflated hydrogen or helium balloons provide a detailed vertical profile of atmospheric properties (e.g., pressure, temperature, relative humidity, and wind speed). How often these balloons are released, the measurement accuracy and sampling, and the maximum height vary strongly depending on location and operator [Parker et al., 1995]. We use balloon-sounding data provided by the University of Wyoming over Acapulco, to calculate refractivity, which when integrated according to equation (3) gives the tropospheric delay ( Figure 1a). It can be observed that delays increase at lower altitudes as the signal has traveled farther through the troposphere. With an interferogram being the phase difference between two SAR acquisitions, the atmospheric delay depends on the change in refractivity rather than the total refractivity. Figure 1b shows the relative tropospheric delays, i.e., the difference between all pairs of tropospheric delays shown in Figure 1a, with the darker color indicating more overlapping delay curves. Variations in interferograms in flat regions are mainly due to the wet component of the refractivity, i.e., the water vapor; the hydrostatic component often appears as a smooth signal due to the large spatial scale of high-and low-pressure fields [Hanssen, 2001]. In regions with significant topography, the hydrostatic component will also cause a correlation between the phase and topography [Elliott et al., 2008].
Different approaches have been documented to estimate the tropospheric signal in InSAR data, e.g. by using weather models [Pinel et al., 2011;Wadge et al., 2002;Walters et al., 2014], GNSS data [Onn and Zebker, 2006;Williams et al., 1998], spectrometer measurements [Li et al., 2006], or by combining weather models and spectrometer data [Walters et al., 2013], or GNSS and spectrometer measurements [Li et al., 2009;Puysségur et al., 2007]. While all of these techniques can, under certain conditions, reduce the tropospheric signal, they are often limited by the lower spatial resolution of the auxiliary data. Alternatively, tropospheric signals can be estimated from the correlation between the interferometric phase and the topography, in a non-deforming area [Cavalie et al., 2007;Wicks et al., 2002], or by filtering the InSAR data in space and time [e.g. Hooper et al., 2012]. The former method assumes a linear relation between the interferometric tropospheric delay and the topography (Δ tropo = K Δ h + Δ 0 ), estimated from data in the non-deforming region. K Δ is a constant relating the interferometric tropospheric phase to topography, and is used to compute the tropospheric signal throughout the full interferogram. Δ 0 can be neglected as it merely represents a constant shift applied to the whole interferogram.
There are two limitations to this linear height correction technique. First, a nondeforming region is required. To overcome this, Lin et al. [2010] developed a multiscale approach, which also assumes a linear relation, but which can be applied in a deforming region by using a spatial frequency band insensitive to deformation. This method makes use of the fact that the linear relationship holds for all spatial frequencies, but the deformation may only be significant in a certain frequency range. Alternatively, Elliott et al. [2008] remove a preliminary estimate of the deformation displacements prior to estimating K Δ . The second limitation follows from the assumption of a single relationship between phase and topography over the whole interferogram, as this does not account for the spatial variation of the tropospheric properties. This can be significant over large areas (greater than tens of kilometers), especially when the area includes different climatic zones. When applying a linear estimation over multiple small windows, Δ 0 will vary and cannot be neglected. However, Δ 0 is biased by other phase contributions and cannot be estimated from a spatial band insensitive to deformation. An exception may exist for the longest wavelength band, provided it is not contaminated by deformation, but this is not usually the case in areas of tectonic deformation.
In this paper we present a new tropospheric correction technique that can be applied to deforming regions and which allows for a spatially varying relationship between phase and topography. We apply our power law methodology to correct the interferograms in a tectonic case study over Mexico, where slow slip deformation is present at the longer spatial wavelengths (∼150 km) and where the tropospheric signals mask the tectonic signal. The tectonic implications of the results are presented in a companion paper [Bekaert et al., 2015]. Figure 2 shows a simplified schematic of the tropospheric delays computed from weather balloon data shown in Figures 1a and 1b. While the tropospheric delay continues to decrease with increasing height, the relative tropospheric variations of the delay curve are important only up to a certain altitude h 0 , typically around 7-13 km, after which the relative delays do not differ significantly ( <0.5 cm) and converge to zero. Consequently, the interferometric tropospheric phase at this reference height is approximately zero. In other words it is only necessary for the interferometric delays to describe the delay curve up to h 0 , as any delay accumulated at this height c will cancel out in the interferometric computation.

Power Law Tropospheric Correction
By plotting the log-log relationship between the mean tropospheric delay and the topography constrained by the reference height (h 0 − h) (black solid line in Figure 1c), we find the fit to be approximately linear indicating a power law relationship. We can find the power law exponent, , from the slope of this linear fit. We can therefore approximate the tropospheric phase delay tropo as So we can rewrite equation (4) as, where K ′ Δ is a spatially varying and unknown coefficient, describing the relation between topography and tropospheric phase. Both h 0 and can be estimated empirically from balloon sounding data or weather model data and are assumed to be constant for a given area. In time, seasonal and long-term variations can be observed for both h 0 and , shown in Figure 3 for Acapulco in Mexico. In this study an average over the whole time period is assumed, alternatively each interferogram can be evaluated using the average h 0 and for the master and slave acquisition dates. The success of the power law method is limited to the cases in which the relative delays can be reasonably approximated by a power law function. From a regional analysis of weather model data for our test region in Mexico ( Figure S1 in the supporting information), we found that, on average, ∼80% of the relative delay curves can be fit by a power law function with a correlation coefficient larger than 0.9.
The interferometric phase is a superposition of multiple signals including deformation, tropospheric, and ionospheric delays, and errors due to, for example, incorrect orbits and unwrapping errors. Estimating K ′ Δ from an interferogram where all these signals are present is not trivial. However, following Lin et al. [2010], we can take advantage of the fact that the tropospheric signal is present at all wavelength scales to estimate K ′ Δ from a spatial frequency band relatively insensitive to these other signals. In other words, we can estimate K ′ Δ by applying equation (6) to band-filtered phase and topography, where the band should be selected such that the contribution from the other signals is negligible. We allow lateral variability by splitting the study area into multiple small windows and estimating K ′ Δ locally. Figure 4 shows a step-by-step example of the tropospheric delay estimation for an interferogram over Mexico spanning 31 December 2004 to 16 December 2005. We used a mean power law decay, , of 1.4 and a reference height, h 0 , of 7 km, estimated from sounding data at Acapulco between 2004 and 2010 (http://weather.uwyo.edu/upperair/sounding.html), to scale the topography (step 1). We found the mean of h 0 and estimated at Acapulco to be consistent with estimates of multiple sounding stations at distances up to 1000 km from Acapulco, with h 0 in the range 6.89 ± 0.02 km and in the range 1.35 ± 0.02. This supports our simplification of using a single mean value for h 0 and over our InSAR region. A 2-8 km bandfilter was chosen to avoid contamination of the tectonic slow slip deformation signal (step 2, Figure 4), which from an initial InSAR and GNSS analyses we found to be around 150 km wavelength [Bekaert et al., 2015]. We estimate K ′ Δ locally in approximately 50 km square windows that overlap by 50%. An example of the K ′ Δ estimation for an individual window is shown in the scatterplot (step 2, Figure 4), where K ′ Δ corresponds to the slope of the red dashed line. Once K ′ Δ is estimated for each window we interpolate to all points of the interferogram (step 3, Figure 4); the interpolated value for each point is calculated as the weighted average of the estimated K ′ Δ values, with the weighting being inversely proportional to the uncertainty of the K ′ Δ estimate and the distance of the window to the point. The tropospheric signal is then calculated by resubstituting the estimated values into equation (6), (step 3, Figure 4). Note that in this example the polarity of K ′ Δ switches between the center of the image compared to the north and the south. A single K Δ , estimated over the whole interferogram (i.e., the linear method), would not be capable of capturing this variation in space.
BEKAERT ET AL.  Step 1 after bandfiltering in the 2-8 km spatial wavelength with a window over which K ′ Δ is estimated locally, represented by the slope of the red dashed line. (step 3) The spatial map of K ′ Δ obtained after estimation in local windows and extrapolation to all points and the corresponding tropospheric delay. The last window shows the tropospheric delay estimate when assuming a linear relationship. The input interferogram and the estimated tropospheric delay are referenced with respect to the mean of the whole region.

Results
We tested the power law correction individually on 17 Envisat (C band) interferograms acquired over Mexico ( Figure 5), for which the perpendicular, temporal, and Doppler centroid baselines are summarized in Table 1, and the unwrapped time-series interferograms are shown in Figure 6a. All SAR data have been focused using the ROI_PAC software [Rosen et al., 2004]. Interferometric processing was done using the DORIS software [Kampes et al., 2003]. We applied time series InSAR processing using the StaMPS software [Hooper, 2008;Hooper et al., 2012] to reduce noise and estimate a digital elevation model (DEM) correction. Note that time series processing is not required, as the power law correction method is applied to individual interferograms. A summary of the applied processing parameters are contained in Table 2. Tropospheric delays introduced by water vapor variations in the lower part of the troposphere and topographic variations, up to ∼3.5 km elevation within our test region, mask the much smaller tectonic slow slip displacement signal [Hooper et al., 2012], which should be visible in images acquired after December 2006. Figure 7 shows the relationship between the unwrapped phase and the topography, with K Δ estimated in 50 km segments along the satellite flight direction. The relationship varies spatially and reverses in sign in some cases, indicating the need for a spatially varying tropospheric correction.
On average we find a strong local relationship |K Δ | of 2.8 rad/km (1.5 cm/km). The conventional linear tropospheric correction as shown in Figure 6b, was estimated assuming a simple phase-topography relationship for each interferogram and therefore does not allow for spatial variability. After applying the linear correction, we find that 5 out of 17 interferograms show a decrease in local correlation with topography Δ|K Δ | (blue solid line in Figure 7). A maximum reduction of 2.5 rad/km (∼1.13 cm/km) can be observed for 15 April 2005, for example. While the linear correction can work well in specific situations, we find that on average it increases the local correlation with topography by ∼0.4 rad/km (∼0.18 cm/km). While the local correlation increases, the correlation based on the full interferogram and the topography on average decreases from 3.3 rad/km before correction to −0.7 rad/km after the linear correction. The increase in local correlation is due to the spatial variation of the troposphere, which cannot be accounted for using the linear method. By applying the power law correction method locally, we are able to account for spatial variation, as can be seen from the variation in the spatial maps of K ′ Δ shown in Figure 6c.    [Kampes et al., 2003] and StaMPS (S) [Hooper et al., 2012].
The corresponding power law tropospheric signal we estimate for each interferogram is shown in Figure 6d. After tropospheric correction using the power law method (Figure 6e), we find an average reduction in local correlation with topography of 1 rad/km (∼0.45 cm/km) (red solid line in Figure 7), about one third of the signal, with the strongest reductions for our 11 March 2005 (2.7 rad/km or 1.2 cm/km) and 31 December 2004 (2.4 rad/km or 1.1 cm/km) interferograms. In total, we find a decrease for 15 interferograms (average decrease of 1.2 rad/km or 0.5 cm/km), while for two interferograms a small increase can be observed (average increase of 0.25 rad/km or 0.1 cm/km), outperforming the conventional linear correction method. Based on the whole interferogram, we find the power law to decrease the correlation with the topography to an average of −0.2 rad/km, 0.5 rad/km better than the linear method. We find that the power law method does not capture the full signal near the coast, which is possibly due to the smaller topographic height range available to estimate the linear relation between the interferometric phase and the power law scaled heights. Also, in the north of our study region, rapid subsidence of Mexico City contains deformation in the spatial band of 2-8 km that could contaminate the tropospheric estimates.
After tropospheric correction we find that the RMSE of the residuals between the displacements and those predicted by the time series model decrease (Figure 8), where the time series model is defined as a linear trend due to interseismic loading, and a jump d SSE at the time of the slow slip event. A detailed description is included in Bekaert et al. [2015]. Close to the reference area, A, we find minimal change, as expected, with a slight RMSE increase of ∼0.1 cm; farther away we find a strong decrease in RMSE of ∼0.6 cm and ∼1 cm for regions B and C, respectively. During the tropospheric correction, we make no a priori assumptions about the form of the time series model.
Figures 9 (left) and 9 (middle) compare the LOS GNSS and InSAR slow slip estimates from the time series model along the P1 profile before and after tropospheric correction. With respect to our estimated slow slip displacement signal, we find the gradient SSE between the slow slip displacements and the topography to decrease on average by 1.1 cm/km (Figure 9c), where d SSE = SSE h + constant, especially toward the center and north of our InSAR study region.

Discussion
Our time series estimate of the slow slip deformation signal after tropospheric correction compares well with another InSAR study by Cavalie et al. [2013]. While the same InSAR track was analyzed, Cavalie et al. [2013] applied a different processing strategy by stacking eight interferograms spanning the slow slip event, to reduce tropospheric effects. Stacking is limited by the number of independent interferograms and how uniform the acquisitions are sampled in time; it assumes that the random component of the atmosphere is reduced by averaging of interferograms, with a noise reduction approximately ∼ √ n for n independent interferograms [Zebker et al., 1997].  Comparing the local slope K Δ , where Δ = K Δ h + constant, between the interferograms and the topography, along profile P1 in 50 km segments, using all data. K Δ is estimated before (black solid line) and after correction for tropospheric signals using both a linear (blue solid line) and power law (red solid line) correction method ( Figure 6). For the linear method, 5 out of 17 interferograms show a decrease in local slope Δ|K Δ |, with a maximum reduction of 2.5 rad/km (∼1.13 cm/km) on 15 April 2005. For the power law method, 15 out of 17 interferograms show a decrease in slope, with a maximum reduction of 2.7 rad/km (∼1.2 cm/km) for the 11 March 2005 interferogram. On average the linear correction does worse than the power law with an increase in local slope of ∼0.4 rad/km (∼0.18 cm/km) compared to a decrease in local slope of 1 rad/km (∼0.45 cm/km). (a-d) The scatterplots for the local correlation between phase and topography before and after correction for the linear and power law methods.    Figure 9. Comparison of the estimated slow slip surface displacements, accumulated over 8-9 months, (left) before and (middle) after tropospheric correction with the estimates from GPS (red 2 error bars). (right) The local gradient SSE , where d SSE = SSE h + constant, between the slow slip deformation signal and topography before and after tropospheric correction, estimated in 50 km segments along profile P1, using all data. Before correction, the signal is strongly correlated in the middle and north of the InSAR study region. After correction, this local correlation is decreased, with on average a reduction of 1.1 cm/km. Light markers (Figure 9, middle) indicate the region around Mexico City, where subsidence contaminates the tropospheric correction estimates at the selected spatial frequency band of 2-8 km, which was chosen to avoid contamination of the long wavelength (∼150 km) slow slip signal.
We define the power law reference height h 0 to be the lowest height above which the relative delays do not differ significantly ( <0.5 cm). A tighter constraint of increases the altitude of h 0 . For example, when using <0.05 cm we find h 0 and to be respectively 9 km and 1.5. We investigated the impact of wrongly assumed power law coefficients, by comparing the estimated power law delay used in this study (h 0 = 7 km and = 1.4) with that estimated from coefficients with a tighter constraint (h 0 = 9 and = 1.5). We found the impact to be small, with an average RMS difference between the two corrections of 2 rad or ∼0.9 cm, similar to the ∼1 cm line-of-sight accuracy of MERIS [Walters et al., 2013].
In order to generalize the proposed power law relationship for performing tropospheric corrections in InSAR data (individual or time series) in other study regions, both the reference height, h 0 , and the power law coefficient, , should be reestimated as they might vary somewhat depending on the climate, season and time of day. While sounding data give an accurate representation of the atmosphere, they might not always be available near a given study region for a given day. Alternatively, weather model data like the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-I or the North American Regional Reanalysis (NARR) product can be used. These would also allow for spatial variation of the power law coefficients, which could further improve the tropospheric delay estimation.
To separate deformation and tropospheric signals, a frequency band insensitive to deformation is required. As the troposphere is present at all spatial wavelength scales, a potential option is to perform a statistical analysis of multiple spatial bands to identify those bands that are consistent and likely to represent the troposphere. Alternatively, one can compare the delay maps as estimated from power law for each spatial band with an independent estimate from MERIS or weather models to identify the best bands that represent the troposphere. Once identified, the bands can be applied throughout the whole data set. Failing to select a correct frequency band could cause deformation to leak into the estimated atmospheric correction.
The selection of the window size is a balance between spatial variability of the troposphere properties and the reliability of the estimate for K ′ Δ . Smaller windows are likely to have a smaller topographic range, making it more difficult to get a robust estimate of K ′ Δ , while a too large window will not allow K ′ Δ to vary  Figure 10a shows that the "optimal" window size is a balance between a window small enough that allows for a spatially varying K ′ Δ , while large enough to provide a reliable estimate of K ′ Δ . Figure 10b gives the delays assuming a fixed window size of ∼90 km and shows that RMSE decreases with increasing topographic range. As input, we use a simulated interferogram containing only tropospheric delay computed from the Weather Research and Forecasting model, at the acquisition times of the 16 March 2007 interferogram. The colored lines represent the different spatial band filters (in meters) used in the power-law estimation.
in space. The selection of the window size will therefore be specific depending on the region of application. We tested the sensitivity to window size over our study region, by simulating the relative delay for the 16 March 2007 interferogram using the Weather Research and Forecasting (WRF) model, and then by estimating the delay from the simulated interferogram using different window sizes. The RMSE between the estimated and simulated delay is plotted in Figure 10a, showing a minimum for a window size of ∼90 km. Figure 10b shows the RMSE when the window size is held constant at 90 km, and the topographic range is allowed to vary instead. An increased RMSE can be observed with decreasing topographic range. Note that we assumed the same spatial variability of atmospheric conditions for all topographic ranges, whereas in reality, topography influences the weather conditions. We did not find clear trends with respect to the band filter size. In general, the largest bandwidth is limited by the spatial extent of the study region, while the smallest bandwidth is limited by the resolution of the data. In Mexico, the delay differences are dominated by the wet component, which are on average 9.4 cm in the radar line of sight, compared to 0.8 cm for the hydrostatic component. In other regions, it may be possible to include a more complex function that accounts for different scale heights and exponents for the hydrostatic and wet components. This will be at the cost of adding additional degrees of freedom. In doing so, both components would need to be untangled as each of them partially correlate with the topography, and we suspect that this would be challenging in most circumstances.
We compared the power law method with the conventional linear method and found the power law method to reduce the local correlation with the topography by on average 1 rad/km (∼0.45 cm/km), compared to an increase of ∼0.4 rad/km (∼0.18 cm/km) for the linear method. While different tropospheric correction methods exist, each with its own limitations, the integration of multiple tropospheric observations or correction techniques can be used to complement each other. For example, by using the temporal information certain observations can be used to constrain the others. In case of Envisat, GNSS observations and measurements from the Medium Resolution Imaging Spectrometer onboard Envisat can be used to constrain weather models and phase correction methods, while in return the phase-based methods could be used to cope with cloud cover. Alternatively a weighted inversion of the different observations can be used to estimate the tropospheric properties (pressure, temperature, and humidity), or a base pressure and an integrated quantity like precipitable water vapor, which then can be used to compute the tropospheric delay.

Conclusions
Previously, tropospheric corrections estimated from interferometric phase have been estimated over a full interferogram, in local windows in a band insensitive to deformation, or using a nondeforming region. While a local linear correction allows for the estimation of the slope-relating interferometric phase and topography, the intercept will be biased by other signals. Failing to include this intercept will lead to an incorrect estimate of the tropospheric delay. We have developed a power law relationship that can be applied locally, allowing us to estimate a spatially varying signal from the interferometric phase. Allowing for spatial variation in the tropospheric properties becomes more important for larger spatial data sets (greater than tens of kilometers) as pressure, temperature, and relative humidity can change significantly over these length scales. We tested the power law technique over Mexico, where interseismic and slow slip deformation occurs on a long wavelength (∼150 km). After correction, we find a reduction in the topography-correlated signals and an improved correlation between the GNSS and InSAR estimated slow slip surface displacements.