Modeling Orbital Error in InSAR Interferogram Using Frequency and Spatial Domain Based Methods

Synthetic Aperture Radar Interferometry (SAR, InSAR) is increasingly being used for deformation monitoring. Uncertainty in satellite state vectors is considered to be one of the main sources of errors in applications such as this. In this paper, we present frequency and spatial domain based algorithms to model orbital errors in InSAR interferograms. The main advantage of this method, when applied to the spatial domain, is that the order of the polynomial coefficient is automatically determined according to the features of the orbital errors, using K-cross validation. In the frequency domain, a maximum likelihood fringe rate estimate is deployed to resolve linear orbital patterns in strong noise interferograms, where spatial-domain-based algorithms are unworkable. Both methods were tested and compared with synthetic data and applied to historical Environmental Satellite Advanced Synthetic Aperture Radar (ENVISAT ASAR) sensor and modern instruments such as Gaofen-3 (GF-3) and Sentinel-1. The validation from the simulation demonstrated that an accuracy of ~1mm can be obtained under optimal conditions. Using an independent GPS measurement that is discontinuous from the InSAR measurement over the Tohoku-Oki area, we found a 31.45% and 73.22% reduction in uncertainty after applying our method for ASAR tracks 347 and 74, respectively.


Introduction
The precise orbiting position information of space-borne Synthetic Aperture Radar (SAR) systems is of great significance in many Interferometric Synthetic Aperture Radar (InSAR) applications, especially in the case of ground motion monitoring [1][2][3].The InSAR technique uses two sensors, carried on satellites, with slightly different incidence angles, to measure a displacement along the radar line-of-sight (LOS) between two SAR acquisitions.The theoretical accuracy can be as high as a centimeter to as low as a millimeter.This technique has been widely used to identify many geophysical processes that usually cause long-wavelength crustal deformation, including strain accumulation along locked continental faults, coseismic deformation caused by the occurrence of faulting in the lithosphere, and postseismic deformation caused by afterslip and viscoelastic relaxation [4][5][6].
The radar, the ground, and the intervening medium create a compound observing system; within this, there are several potential sources of errors that have a detrimental impact on the accuracy of InSAR measurements.In particular, satellite orbital errors [7], temporal and spatial decorrelation effects [8,9], atmospheric screens [10], and high-deformation gradients are the main limitations [11].As far as orbital errors are concerned, they cause long-wavelength phase contributions to interferogram and are often referred to as 'the phase ramp'.The pattern of the phase ramp usually depends on the satellite-state-vector error, especially on the radial and cross-track components of the orbital error [7].The influence of orbital errors on the final accuracy of deformation products is largely contingent on the SAR instruments, i.e., the trajectory of the satellite orbit, the radar frequency, and the degree of overlap between the phase ramp and the deformation signals.
Several efforts have been made during the last two decades to mitigate the orbital error, and most of them start with spatial domain analyses of InSAR interferograms.A simple but very effective solution is de-ramping from the InSAR interferogram using a linear or quadratic surface fitting [3,12].The phase ramp is either subtracted from the original phase or used to refine the spatial baseline to infer a revised interferogram [13,14].Given a wrapped-phase pattern, the unwrapping operation is always required before fitting.The accuracy of the estimated coefficients in the polynomial model relies on the quality of the observations to be fitted.Error propagation in the unwrapping procedure over fast decorrelation areas with low coherence may distort the signal and, therefore, contribute to surface artifacts.To solve this problem, pixel offsets between bi-temporal-SAR data can be estimated using the cross-correlation algorithm; the baseline is then re-estimated to compensate for the residual phase in the interferogram [15].Despite many efforts, the accuracy of the estimated offsets, from meter to decimeter, restricts the correction of the baseline components.Kohlhase et al. (2003) suggest counting the fringes caused by orbital errors according to phase gradients in differential interferograms and to adjust the trajectories of satellite orbits [16].This does not work for coseismic scenes, where a stronger deformation causes dense fringes [7].
Advances related to orbital error correction refer to the time-series-SAR dataset.One such suggestion is to independently estimate the phase ramp on each interferogram using least squares scheme.To ensure consistency in the interferometric network, the polynomial coefficients are refined in a network sense by time-series inversions [17,18].Considering the inseparability of phase ramp and long-wavelength deformation signals in the spatial dimension, the alternative is to estimate both deformation and orbital errors by exploiting spatio-temporal characteristics of both quantities [17,19].However, resolving large, linear systems of algebraic equations, error propagation, and underdetermined problems have become the main challenges associated with these methods.
One promising alternative to separating long-wavelength displacement signals and orbit errors is to employ external data, for example, GPS measurement [20][21][22].GPS displacements located in the SAR scene are projected onto the line-of-sight direction according to the incidence angle and the heading of the SAR geometry.The phase ramp is then estimated by minimizing the residuals between the InSAR interferogram and the phase inverted by the GPS.The final accuracy of the InSAR deformation product depends on the accuracy of the GPS measurements, the number of GPS stations and their spatial distribution.
In summary, the current problems of orbital error removal are two-fold.The first problem originates from the fact that all long-wavelength signals are treated as orbital errors during the correction without ancillary data.Although a few studies try to use wavelet-multiresolution analysis to classify the different components in the unwrapped interferogram automatically [23], the method is still empirical, as the number of levels of wavelet decomposition must be specified manually according to the features of the interferogram.The second problem involves robust regression.As stated in Reference [24], the model selection is subjective and relies on the distribution and the density of the targets showing high coherence, while errors in unwrapping are likely to mislead surface fitting.In fact, the method of estimating the phase ramp in the frequency domain may avoid the unwrapping error.However, few studies discuss using this method for orbital error correction.In the following section, we focus our analyses on the individual interferogram, as it is essential for time-series inversion and many geophysical applications.
The remainder of the paper is organized as follows: In Section 2, we introduce the frequency and spatial domain based methods, respectively, from the image-processing viewpoint.In Section 3, the results and quantitative comparisons are presented using synthetic and real data.In Section 4, the feasibility of these methods, including their merits and faults, are investigated.In Section 5, conclusions are given.

Modeling Orbital Error in Frequency Domain
The fringe rate of long-wavelength signals in interferogram can be estimated by means of maximum likelihood (ML) frequency algorithms [25][26][27].The idea is that the orbital error has a dominant term in its Fourier expansion.Assuming original interferogram I can be expressed as: where φ de f , φ topo , φ orbit , φ other denote deformation phase, topographic phase, orbital ramp, and other phase components related to the atmosphere and decorrelation, respectively.D is the differential interferogram.After atmospheric correction and/or phase filtering, we rewrite D as a sinusoidal model: where x and y denote the pixel index in radar coordinate, f x and f y are true fringe frequencies along range and azimuth direction, respectively, and ρ is the residual phase.
The Fourier filter can be used to obtain the ML estimate of the average stripe rate in Equation (2).The peak location obtained by maximizing the 2D discrete Fourier transform (DFT) corresponds to the estimated frequencies fx and fy , respectively, and the phase ρ at the DFT peak is approximate to ρ.To better identify the exact frequencies, we increase the DFT frequency sampling by padding the signal in the window with enough zeros.The size of the padded zeros in both directions, n x and n y , can be determined by the Cramer-Rao bound of the variance of estimated frequency [26]: where SNR represents the signal-to-noise ratio, N x and N y are equivalent to each of the image sizes.The de-ramped interferogram can be obtained from the cross multiplication of the differential interferogram and the estimated phase ramp under the assumption fx ≈ f x , fy ≈ f y and ρ ≈ ρ.The significant advantage of this method is its simple computation without the need for phase unwrapping and higher reliability over low-coherence areas (see Section 3.1).However, this method can only estimate the linear terms of orbital error.For non-linear fringe patterns, which are usually present in long-strip images, polynomial fitting in the spatial domain is recommended.

Modeling Orbital Error in Spatial Domain
The most popular method used in orbital-error correction is surface fitting with linear or quadratic models, using unwrapped phase observations.Its general form can be defined as [23]: where z is the polynomial model, x and y indicate the pixel index in radar coordinates.To solve the equation for the unknown coefficients P 1 , P 2 and P 3 , the least squares scheme is used.The solution can be written as: where P T = [ P1 , P2 , P3 ] is the estimated vector in which the elements correspond to the coefficients in Equation (5).Z is the observations with size N, and A is the N × 3 matrix containing the coordinates x and y.However, there are several drawbacks to using this method.Firstly, Equation ( 6) is not robust to abnormal values, such as unwrapping errors.In addition, the observations contribute to the equivalent weight without the consideration of the interferogram qualities.Secondly, the local deformation will contribute to the polynomial coefficients with the increase of the dimension in P T .Finally, the determination of the order of the polynomial model is not clear yet and is usually empirical.
A more sophisticated method is therefore needed.

Preprocess: Multi-Looking and Manually Masking
An increase in looks can enhance the quality of long-wavelength signals and simultaneously mitigate the effect of the short-wavelength signals during the regression.On the other hand, the decreased dimension can reduce the computational burden for unwrapping and the iterative least squares scheme.For those regional deformation signals, manual masking should be undertaken.Although this step is simpler than wavelet decomposition [23], it is still effective without expert knowledge and computational complexity.In addition, masking is important for areas with a high phase gradient, where phase aliasing causes unexpected unwrapping errors [11].
The parameters n and m should be carefully determined to prevent overfitting in terms of the spatial distribution of observations.

Iteratively Reweighted Least Squares Fitting
The main disadvantage of ordinary least squares is that there is a constant deviation in the errors.For observations with different coherence magnitudes, the method of weighted least squares should be used.Moreover, iterative behavior can be employed to reduce the influence of outliers on regression.We designed an iterative version for the least-squares solution according to the phase standard deviation (STD) and robust Bi-square weight.The modification of the parameter estimate P in Equation (6) can be rewritten as: In this paper, we define weight W in two steps.The prior weight V = {v 1 , v 2 , . . . ,v N }, defined by interferometric phase STD, can be simply written as [27]: where L is the looks and γ i is the coherence observations.Compared with the coherence, phase STD σ i is a more appropriate weight, as it takes looks into account.The function model Equation ( 9) is a close approximation to the complicated stochastic model when looks are greater than four [27].After initial regression, the robust Bi-square weight B = {b 1 , b 2 , . . . ,b N }, which minimizes a weighted sum of squares, is used during the iteration [28]: where r i is least-squares residual for each observation from the previous iteration, c is a tuning constant (4.685), and s is the robust variance given by MAD/0.6745,where MAD is the absolute deviation of the residuals from their median; h i is leverage that adjust the residuals.The total weight function W, therefore, is: where l denotes the number of iterations.The iterations will stop if the fit converges (minimum change in coefficients <10 −5 ), or the maximum number of iterations allowed for fitting is reached (400).

Model Selection
The presence of the unwanted terms in Equation ( 7) and the overfitting lead to an unexpected bias or variance in interferograms.A trade-off should be made by taking advantage of K-cross validation, which has been widely used in machine learning [29].Specifically, we obtained some test data from the same distribution and picked appropriate m and n in Equation (7) to minimize the same sum-of-squares difference that we used for fitting to the training data.This can be achieved by making several different splits in data Z.Each subsample is used once for testing, and the rest is used for training.The arrangement of this method is described as follows.

1.
Split data Z into K subsamples with equivalent size N.

2.
For k = 1, 2, . . ., K, set validation data Z test to be the k th subsample, and training data Z train to be the other K − 1 subsamples.

3.
Fit each model to Z train and evaluate its performance on Z test through weighted root-mean-square error (WRMSE).

4.
Pick m and n that leads to minimum WRMSE by averaging K results.
Uniform sampling without replacement is used to separate data Z, and K = 10 is set throughout this paper.Considering the different qualities of the validation data, we define WRMSE as:

Synthetic Data
In this section, the performance of the orbital-error corrections are evaluated and compared using a simulation.For notional convenience, we refer to the method in Equations ( 3) and ( 4) as DFT, and the method developed in Section 2.2 as adaptive ADP-Poly.
Three phase components, including deformation, random noise, and orbital errors were simulated according to the ENVISAT ASAR geometry.The surface deformation caused by a finite rectangular source was first simulated by the Okada 85 model [30] and then inverted to an interferometric phase (Figure 1b).To add phase noise, the coherence map was simulated under three principle sources of error, i.e., the Doppler, spatial, and temporal decorrelations [25,31], in which the temporal decorrelation was created using an isotropic-2D-fractal surface with a power law spectrum [32].On the basis of the simulated coherence map under different looks (Figure 1a), the phase-STD map was obtained using the stochastic model (4.2.24) and (4.2.26) in [31].The phase-noise map was finally generated using pointwise multiplication of the phase-STD map and the standard normal distributed random numbers.The linear and non-linear phase ramps were added to the orbital errors (Figure 1c,d).The combination of all components in synthetic interferograms is given in Figure 2. To test the methods over fast decorrelation areas, we adjusted the coherence magnitude to an average value of γ = 0.2 under single look, and then added the phase-noise to the interferogram in Figure 2a.Likewise, we set γ = 0.4 under 2-looks to test the capability of both methods on moderate noise scenes, as shown in Figure 2f.For the frequency-domain-based method, the DFT is applied to noisy interferograms directly, while adaptive Goldstein filtering driven by phase STD [33] and minimum cost-flow unwrapping [34] were used before ADP-Poly surface fitting.
Visually, ADP-Poly could not detect the phase ramp in very noisy scenes and left phase residuals in Figure 2c.By contrast, DFT accurately captured the linear phase ramp.The residual between the truth and estimate confirms the potential of using DFT without phase unwrapping (Figure 2d).When the quality of interferogram became better (Figure 2f), ADP-Poly showed its superior performance at detecting nonlinear features (Figure 2h), while DFT could not completely remove the error and left nonlinear residuals in the interferogram (Figure 2g).
To make the results more statistically significant, we performed 500 simulations for each case and evaluated the phase residuals using root-mean-square error (RMSE).It can be seen from Figure 3 that DFT works well in linear orbital error correction over low coherence areas, and the averaged RMSE is up to 0.16 rad (4.99 rad for ADP-Poly).For nonlinear orbital features, ADP-Poly obtained RMSE of 0.10 rad (2.77 rad for DFT), which is equivalent to ~1 mm in deformation monitoring for C-band InSAR measurements.

Real Data
Four SAR datasets, covering the Datong and Tohoku-Oki areas, were used to validate the methods of orbital error correction, respectively.The information for the interferogram is summarized in Table 1.The pre-processing includes the co-registration to subpixel accuracy with cross correlation algorithm, multi-looking (3 × 15 looks for ASAR data set along range and azimuth, 4 × 4 looks for GF-3 with stripmap model, 9 × 3 looks for Sentinel-1 with TOPSAR mode), the topography-component removal using the 90-m Shuttle Radar Topography Mission v. 4.1 digital elevation model, and the adaptive Goldstein filtering.The atmospheric delay is not considered any further due to the lack of the external data.For this paper, we assume that this component, together with orbital error, comprises the linear/nonlinear phase ramps.

Datong Area
The deformation feature in the GF-3 interferogram is locally distributed because of the dynamics of land subsidence caused by underground mining activities over the Datong area, China.Errors in the GF-3 satellite state vectors can be observed in the differential interferogram (Figure 4c) and cause a linear phase ramp along the azimuth direction.Because of the relatively low coherence over the area (Figure 4b), we employed DFT to remove the long-wavelength signals (Figure 4d) without the phase unwrapping.The results in the amplified sub-regions (1-3) show that mining-related surface deformation could be monitored after removal of orbital errors.The multi-looking GF-3 SAR image is shown in Figure 4a.In contrast to GF-3, the better orbital qualities of Sentinel-1 instruments mitigate the phase ramp in the interferogram (Figure 5c).After burst and sub-swath merging, a nonlinear phase ramp with a relatively large spatial coverage can be observed in the upper right part of image.To remove the phase ramp, ADP-Poly is motivated by the high coherence magnitude (Figure 5b).A ten-fold cross validation was performed to select the optimal order of polynomial model.It can be seen from Figure 6a that the overfitting causes a larger bias for regression with the increase in order in polynomial, while the polynomial coefficient with n = 3 and m = 3 reaches the minimum WRMSE (2.15 rad, corresponding to ~9 mm).The final de-ramping interferogram is shown in Figure 5d, where the spatially distributed subsidence bowls can be clearly seen.The multi-looking sentinel-1 SAR image is shown in Figure 5a.

Tohoku-Oki Area
For long-wavelength displacement signals covering the whole image, external GPS measurements should be used to correct the displacement before phase-ramp estimation.To this end, GPS observations, including total 120 stations, were collected from the Advanced Rapid Imaging and Analysis (ARIA) team at Jet Propulsion Laboratory (JPL) and Caltech [20].During the process, the GPS observations were first projected into the LOS direction using the central unit vector for east, north and vertical directions in descending ASAR (Table 2).Then, all GPS locations were transformed into radar coordinates.The nearest SAR neighbors, with respect to each GPS measurement, were averaged using 3 × 3 boxcar windows, and the weights were determined by means of averaged phase STD with a factor of 1/3.When the phase ramp fitted using ADP-Poly, 90% of the GPS observations located in each SAR track were used to correct the deformation in unwrapped InSAR interferogram, and the remaining 10% were used to evaluate the discrepancies between the GPS observations and the orbit-corrected InSAR coseismic measurements.We present the results of before and after ADP-Poly orbital error correction in Figures 6 and 7 and Table 2.It can be seen that ADP-Poly suggests a quadratic model to fit the ramps for both tracks (Figure 6b,c).The RMSE values, before and after the orbital error correction, are 35.55 cm and 9.52 cm for track 347, and 12.24 cm and 8.39 cm for track 74, respectively, demonstrating the usefulness of the ADP-poly method developed by the authors.After removal of phase ramp, the displacement caused by earthquakes is rewrapped to 11.8 cm for each fringe, as shown in Figure 7.As stated in previous studies [17,20,35], the discontinuity between GPS measurements and InSAR displacement can be attributed to unwrapping errors over very low coherence area and the aftershock and postseismic deformation that were not covered by the GPS observation period.

Discussion
Two methods of modeling orbital error in InSAR interferograms that deploy DFT and the polynomial model with adaptive order determination were presented.DFT is a robust estimator of linear components of the long-wavelength signals, even in strong noise environment.It transfers spatial phase ramps into the frequency spectrum and finds the fringe rate at the peak location of the Fourier transformation.The interpolation of the DFT of the interferogram can better identify the exact fringe rate.A significant advantage of this method is that phase wrapping is not required, and, therefore, the errors caused by the poor quality of the interferograms are avoided.For nonlinear fringe patterns, the polynomial model (ADP-Poly) can compensate for nonlinear components that cannot be estimated with DFT.The development of adaptive order determination for the polynomial model using K-cross validation ensures the model can follow the features of orbital errors more accurately.This step is missed in previous studies, which usually determine the order empirically.
Two problems in orbital error correction should be further pointed out.Firstly, the approach has no capacity to distinguish between orbital errors and long-wavelength atmospheric delays.To our knowledge, this is a common problem in most related studies [12][13][14]20,23].We therefore made the assumption that the long-wavelength signal only represents the orbital error, implying that part of the tropospheric delay will also be removed from the interferogram.A good alternative is to remove the atmospheric delays prior to orbital error correction, for example, by using the Generic Atmospheric Correction Online Service, applying the Iterative Tropospheric Decomposition model [34].
Secondly, manual masking is used to remove local deformation signals, leading to gaps in the local region, and so, no observation will yield contributions during de-ramping.To solve this problem, a deformation free interferogram can be inverted using a source geophysical model, e.g., Mogi's formulation for an infinitesimal spherical source in an elastic half-space, by finding the global minimum in the misfit space using a simulated annealing algorithm [36].To do this, General Inversion of Phase Technique, which estimates parameters in a quantitative model directly from the wrapped-phase data, can be employed [36].

Conclusions
This paper has explored orbital error correction in InSAR interferogram using frequency (DFT) and spatial domain (ADP-Poly) based methods.The results reveal that the modeling of orbital error should follow the features of fringe patterns and the quality of the interferogram.To detect linear features of orbital error, the method based on maximum likelihood fringe rate estimation should be deployed.The main advantage of this technique is that it is workable over fast decorrelation areas without the need of unwrapping.However, only linear components can be estimated using this method.For nonlinear features, ADP-Poly, which determines the order of polynomial automatically using K-cross validation, should be used.We evaluated both methods using historical and modern sensors, including ENVISAT ASAR, GF-3 and Sentinel-1.The experiments show that both historical and modern systems can benefit from our method to varying degrees, leading to clearer ground deformation signals, both locally and globally.

Figure 1 .
Figure 1.Simulated InSAR parameters (a) coherence map; (b) noise-free differential interferogram.The black frame denotes the masked area for surface fitting in spatial domain; (c) linear orbital error and (d) nonlinear orbital error with non-uniform phase gradient.

Figure 2 .
Figure 2. Simulated linear and nonlinear orbit errors in InSAR interferograms.(a) noise added interferogram with averaged coherence γ = 0.2 and single look; (b) orbit-corrected interferogram using DFT; (c) orbit-corrected interferogram using ADP-Poly; the difference between truth and phase ramp estimated from DFT (d) and ADP-Poly (e); (f) noise added interferogram with averaged coherence γ = 0.4 and 2-looks; (g) orbit-corrected interferogram using DFT; (h) orbit-corrected interferogram using ADP-Poly; the difference between truth and phase ramp estimated from DFT (i) and ADP-Poly (j).

Figure 3 .
Figure 3. Accuracy assessment for DFT and ADP-Poly under different noise degree and phase ramp patterns.

Figure 4 .
Figure 4. Orbital error correction for GF-3 sensor over Datong area: (a) SAR intensity acquired on 1 April 2017; (b) corresponding coherence map; (c) original interferogram with linear phase ramp; (d) de-ramping interferogram estimated from DFT.The sub-image 1-3 corresponds to frames shown in Figure 4d.

Figure 5 .
Figure 5. Orbital error correction for Sentinel-1 sensor over Datong area: (a) SAR intensity acquired on 15 October 2015; (b) corresponding coherence map; (c) original interferogram with nonlinear phase ramp; (d) de-ramping interferogram estimated from ADP-Poly.The sub-image 1-2 corresponds to frames shown in Figure 5d.

Figure 7 .
Figure 7. Coseismic deformation interferograms generated using data from ASAR descending tracks 347 and 74, before (a) and after (b) orbital error correction.The locations of the cities (khaki square) and GPS stations (green triangle) that are used to remove the orbital errors.The arrows indicate the biases in the InSAR results as calculated from the GPS observations.Each color cycle represents 11.8 cm of LOS displacement toward or away from the satellite.The epicenters of the magnitude 9.0 main shock and a magnitude 7.4 aftershock are marked as red stars.

Table 2 .
Deformation validation and RMSE before and after orbital error correction over Tohoku-Oki area.