Accuracy of a Model-Free Algorithm for Temporal InSAR Tropospheric Correction

: Atmospheric propagational phase variations are the dominant source of error for InSAR timeseries analysis, generally exceeding uncertainties from poor SNR or signal correlation. The spatial properties of these errors have been well studied, but their temporal dependence and correction have received much less attention to date. We present here an evaluation of the magnitude of tropospheric artifacts in derived time series after compensation using an algorithm that requires only the InSAR data themselves. The level of artifact reduction equals or exceeds that from many weather model based methods, while avoiding the need to access fine-scale atmosphere parameters globally at all times. Our method consists of identifying all points in an InSAR stack with consistently high correlation, and computing, then removing, a fit of the phase at each of these points with respect to elevation. Comparison with GPS truth yields a reduction of 3, from an rms misfit of 5-6 cm to ~2 cm over time. This algorithm can be readily incorporated into InSAR processing flows without need for outside information.


Introduction
Interferometric radar analyses follow from interpretation of the precise time delays and phase shifts in radar echoes as distance measurements, and then relating those geometric distances to topography or deformation of the surface. We generally assume that the signals propagate at known constant velocity to convert the delays to distance. However, if the signals propagate through the Earth's spatially and temporally inhomogeneous atmosphere, which has a slightly higher index of refraction than free space, the velocity is lowered slightly, leading to variable delays in space and time which contaminate the observations.
In modern radar systems, atmospheric propagation variations are usually the largest artifacts present in interferograms, dominating SNR and decorrelation noises. Typical errors in a single interferogram from the variable water vapor part of the troposphere are ~2 cm rms, and can be several times greater across a time sequence. Compare this to the uncertainty in distance from phase estimation: at C-band, an SNR of 20 dB yields a phase or distance error of about 3 mm, even if a minimal one look is used. Hence correcting for variations in the atmosphere is necessary in order to fully exploit the potential of spaceborne InSAR.
Many studies have addressed this problem specifically as regards InSAR analysis. Early work identified these phase artifacts [1][2][3][4], with [5][6][7][8] presenting the mapping from atmospheric parameters to InSAR phase. [5] also noted that the artifacts could most easily be compensated by stacking multiple observations, because models of the atmosphere at the time were not at sufficiently fine resolution to correct the measurements. Quite a few papers have examined the use of both water vapor models (e.g., [9][10][11][12][13]) or empirically derived corrections from GPS and other instrumentation [13][14][15][16][17][18] to compensate InSAR images for propagation path delays--a comprehensive review is found in [13]. At the risk of oversimplifying, the conclusions from these works are similar: the algorithms are effective except when they are not, and they tend to reduce the total artifact energy by about a factor of two.
Much of the effort in works addressing compensation for propagation artifacts has focused on spatial variation and how that may be minimized in individual interferograms. Temporal variations have been far less studied as regards InSAR analyses. Both [19] and [15] present sequences of GPS zenith wet delay that show variability with time, furthermore [15] and [20] compute temporal structure functions and power spectra of GPS zenith wet delay consistent with Kolmogorov-like -5/3 power law slopes. The spectra can grow quite large for long temporal intervals, consistent with decorrelation of the atmospheric delay over days-plus interferograms.
Some of the variations are well known to be readily compensated. The above and other [21][22][23] works generally show simple (approximately linear) relationships of atmospheric delay with elevation, hence proposing that subtraction of that dependence "flattens" deformation interferograms with respect to topography. [24] further proposed using that correction routinely in the production of interferograms, in their case using the atmospheric term resulting from persistent scatterer analysis. Here we validate and extend these efforts by assessing a set of Sentinel-1 InSAR data to verify and quantify the accuracy of temporal deformation observations using an algorithm that averages the elevation-residual phases in every interferogram and compensates for the long-time component of atmospheric variation. We show that a factor of three temporal artifact reduction can be obtained by an algorithm relating phase to elevation using only the InSAR data themselves, so that the problems of assembling global high-resolution water vapor models may be avoided. More is yet to be gained in future work, likely by combining our proposed approach with some of the model-based solutions, and by ensuring multiple data acquisitions to average away tropospheric variations.

Materials and Methods
Denote the complex amplitude of a unit intensity plane wave at position x in a medium with variable refractive index n(x) and wavelength λ by where the wavenumber k = 2π/λ. The accumulated phase along the propagation path is nearly For a radio signal propagating through free space, n(x) = 1, and the phase shift is directly proportional to path length: The phase shift depends only on the free-space wavelength λ and the distance propagated x. If the signal propagates instead through an atmosphere, n(x) is not constant and an additional phase shift follows. For the Earth's neutral atmosphere, n(x) is always real and slightly greater than one, so we can expand n(x) as 1 + 10 -6 N(x), where N(x), called the refractivity, is the additional refractive index due to the atmosphere. The change in n(x) from the free space value of 1 is very small for Earth's atmosphere, hence the factor of 10 -6 in the definition. In this case eq. (3) becomes where ∆x = 10 -6 N(x) • x represents additional phase shift as a change in effective path length ∆x. This value is often broken into two parts (Goldhirsch and Rowland, 1982): ∆x = (∆x)dry + (∆x)wet (6) where (∆x)dry and (∆x)wet represent the contributions to path length from the "dry" atmosphere and from water vapor. Empirical measurement of these effects (Smith and Weintraub, 1953) shows that (6) may be approximated by where X is the total path length through the atmosphere, P is the atmospheric pressure in millibars, T is the temperature in Kelvins, and e is the partial pressure of water vapor in millibars -all of these quantities may vary with location along the propagation path. The constants preceding each integral are generally valid to within about 0.5% for frequencies up to 30 GHz and normal variations in pressure, temperature, and humidity. Previous work [5] noted that P and T change fairly slowly with location, while the inhomogeneity of water vapor (e) remains significant at small scales [8]. Furthermore, the path length X through the atmosphere depends on topography as propagation to sea level traverses more atmosphere than that to the top of a mountain.
Interferograms for deformation studies are acquired at different points in time, and it is also their temporal phase (delay) difference that presents as tropospheric artifacts. While pressure or temperature may vary slowly with location, because InSAR component images are acquired days apart these quantities can be quite different, and this difference will depend on elevation. The wet troposphere delay will vary both spatially and temporally, with the longer spatial scale components behaving like P and T and the short spatial scale artifacts persisting as highly variable and obvious noises. None of the present water vapor models are very accurate at m-scale resolution, hence only averaging can reduce these.
Fortunately, most of the signal power in the phase distortions is at long scales due to their power law expression. The atmosphere exhibits a rotationally averaged spatial power spectrum of phase delay that decreases approximately as a power law [6,8] with slope -5/3 at length scales from an outer scale the size of the radar scene (~100 km) to the wet troposphere scale height of 1-2 km [6], and of slope -8/3 down to an inner scale of the SAR resolution (~5 m). Integrating the power spectrum of the artifacts using these values shows that 97% of the energy is at the long scales with only 3% at the short scales, hence the value of the correction procedures. As discussed next, the goal of most tropospheric correction algorithms is to remove, as much as is feasible, the longer scale artifacts, and with the smaller scale distortions reduced via averaging [5].

Results
3.1 Correction algorithm. Before presenting our correction algorithm, it is useful to note the steps needed for interrelating a stack of radar interferograms. The processing flow we use here is the following: 1. Geocode and phase compensate single-look complex (SLC) images 2. Form interferogram pairs as desired and unwrap 3. Choose reference point(s) to remove artifacts due to unwrapping zero phase point 4. Remove troposphere using height regression 5. Compute time series using SBAS Steps 3 and 4 define our atmospheric correction algorithm. Step 3 is needed for time series analysis because each interferogram may contain an arbitrary phase constantphase unwrapping consists of computing phase gradients and then integrating the result, thus all derived InSAR phases are relative to the choice of the reference point. Also, if the reference point itself is moving all inferred time series will reflect this motion.
Step 4 involves correcting the set of interferograms for elevation-dependent temporal variations in pressure, temperature, and the long spatial scale component of water vapor.
Step 5, the temporal dependence estimator, can incorporate the final atmosphere-reducing procedure, namely averaging multiple scenes to lessen the variable water vapor artifacts at all scales.
3.2 Reference point selection. If we are to analyze successive temporal measurements of InSAR phase, we must separate motions of the surface from the propagational phases due to the variable atmosphere. Again, unwrapping the interferogram phases in each image loses the constant term because each interferogram consists only of relative values. In most cases a location assumed to non-moving is chosen because the entire derived interferogram would appear to move in unison if the reference point itself is not stationary.
Even if the selected reference point does not move, or has known motion, the reference will undergo apparent motion because it too is imaged through the variable atmosphere. Thus there are two challenges to selecting the reference: we need outside knowledge of the scene to know where the nonmoving points may be, and the reference exhibits apparent motion from phase artifacts. The latter effect may be minimized if we select multiple reference points and average the results. This will compensate for the short length-scale variations, and result in an average phase offset for the image. By subtracting the average offset from its interferogram in the next step, the long wavelength portion will be removed as well. Our approach is to automatically select reference points by finding all of the points in a scene that are reliably estimated, specifically those with high correlation values. We identify these by examining the correlation at each point in all interferograms forming a stack, and choose only those points whose correlation exceeds a threshold in every interferogram. It is important that the references in every interferogram have reliably estimated phase, or else the scene will add arbitrary temporal noise to the stack. Because terrains differ, the appropriate threshold will vary in different stacks. For better correlated scenes, a value of 50-60% is adequate, while for more highly vegetated scenes a value as low as 20-30% may be needed. We typically examine thresholds that identify several hundred to several thousand reference points, ensuring a wide enough spatial distribution to average out spatial artifacts.
We illustrate in fig. 1 reference point selection in a set of Sentinel-1 data acquired over the island of Hawaii. The background image is the average correlation, which ranges from near zero in the water and in the densely forested areas, to nearly unity on the bare lava surfaces. The green squares are the chosen reference points, which are located primarily in the high correlation regions and are well distributed spatially.

Troposphere height regression.
Once we have selected the set of reference points, we compute the phase correction function. A geographically distributed set of reference points averages spatial variations of delay, but the points lie at many different elevation levels, and the propagation path length through the atmosphere will vary. We noted that the correction for spatially diverse variations in P, T, and e will be dependent on elevation, so we compute a regression of the phase at each point against its topographic height rather than simply average all reference phases. Our experience is that a simple linear fit suffices. We compute the elevation fit for each interferogram in the stack, and subtract that function from the observed phases. In this way each interferogram is referenced to a common zero level. This approach will be valid as long as the image area undergoing significant deformation is small compared to the full scene. Figure 2. Sequence of 12 example interferograms (Jan-May 2018) from Hawaii stack, before (left) and after (right) regression correction. The corrected images display much less variation with time than the uncorrected images, hence the derived time series will be more accurate.
The results of the correction can be seen in fig. 2, where we compare several interferograms from the Hawaii stack described above before and after the elevation correction. The corrected interferograms (right) show much more consistent phases than the uncorrected values (left). The corrected interferograms then form the input data set for an SBAS temporal solution.  We compute radar line of sight GPS timeseries for the year 2018, and plot them along with the InSAR timeseries at each GPS location for several SBAS analysis parameters. Fig. 4 shows three example comparisons, one at each of the three regions studied. Complete results are tabulated in Table 1, which also shows the average rms error for each processing case. Solutions are shown both with and without the troposphere correction algorithm. The solutions denoted without default to the standard method of picking one reference point to compensate for the arbitrary offset of each interferogram, in this case near the center of the image. We do not observe significant deformation at this location.  Table 1.  fig. 4) both (with) using the atmosphere correction described here, and (without) using the standard single reference point algorithm. Results are shown for SBAS solutions with maximum temporal baselines of 12, 30, and 100 days. 100 day solutions reduce the atmosphere more through averaging, but lose some signal from undersampling the phase gradient (see text).
In nearly all cases the corrected solution errors (denoted with) are smaller than the single reference approach. Furthermore, the improvement is greatest for the sites near the Mauna Loa summit, demonstrating the need for elevation dependence in the correction. Looking at the mean rms values for each processing stream, it is clear that the gain is greater for shorter maximum baselines. For a longer baseline limit, more interferograms span a given time, and their averaging in the SBAS reduction helps mitigate the atmospheric artifacts. The overall corrected error is not as sensitive to baseline parameters because the longer baselines underestimate deformation in large-motion areas (e.g. stations OUTL and PUOC) due to temporal aliasing of the raw interferometer fringes at the 60 m spatial sampling used here (personal communication, K. Pepin, Stanford University). This added smoothing compensates a bit for the atmospheric noise, albeit at loss of signal.

Discussion
The dominant source of error for InSAR timeseries imaging, phase artifacts from the variable atmosphere, can be compensated with an algorithm that requires only the InSAR data themselves. Comparison with GPS truth suggests that a reduction of a factor of 3 in rms misfit is often possible, which is equal to or better than many weather model based methods, and avoids the need to access fine-scale atmosphere parameters globally at all times. Our method is automatic in the sense that it consists of identifying all points in an InSAR stack with consistently high correlation, and computing a fit of the phase at each of these points in each interferogram with respect to elevation prior to temporal analysis. Thus, this algorithm can be readily incorporated into InSAR processing flows without need for outside information. Future work promises to reduce these noises further if the advantages of local weather data can be reliably added as a final step.
Funding: This research was funded by NASA, grant numbers NNX17AE03G and NSSC19K1485.
Data Availability Statement: All Sentinel data used here are available from the Alaska Satellite Facility as part of the NASA radar DAAC. GPS data may be obtained from the Nevada Geodetic Laboratory at the University of Nevada, Reno.