A Generic Approach to Parameterize the Turbulent Energy of Single-Epoch Atmospheric Delays From InSAR Time Series

The observed phase in time series of interferometric synthetic aperture radar (InSAR) products is a superposition of various components. Differential topography, line-of-sight displacements, and differential atmospheric delays are the main contributions and need to be disentangled to derive accurate digital elevation model (DEM), deformation, or atmospherical products from InSAR. However, isolating the atmospheric component has been proven difficult as it is spatiotemporally highly dynamic and a superposition of two atmospheric states. Here, we propose an approach to parameterize the stochastic properties of the single-epoch atmospheric delay field as a way to define the atmospheric signal. We found that the atmospheric signal of a time series of interferograms can be characterized by structure functions, which can be used to isolate the single-epoch structure functions. Due to the scaling properties of the atmospheric signal, it is then possible to construct a parametric function per SAR acquisition, using two isotropic and three anisotropic parameters. In particular, the isotropic parameters for the short-distance variation and long-distance variation in atmospheric delay can be used to characterize the atmospheric signal. For a test set of 151 Sentinel-1 acquisitions, this results in an atmospheric energy range of about 10 for short-distance scales and about 50 for long-distance scales. Our parameterization demonstrates that we can describe the spatiotemporal variability of InSAR atmospheric delays, which provides a measure for atmospheric noise for individual epochs in deformation time series based on distance and azimuth.

A Generic Approach to Parameterize the Turbulent Energy of Single-Epoch Atmospheric Delays From InSAR Time Series Gert Mulder , Freek J. van Leijen , Member, IEEE, and Ramon F. Hanssen , Senior Member, IEEE Abstract-The observed phase in time series of interferometric synthetic aperture radar (InSAR) products is a superposition of various components. Differential topography, line-of-sight displacements, and differential atmospheric delays are the main contributions and need to be disentangled to derive accurate digital elevation model (DEM), deformation, or atmospherical products from InSAR. However, isolating the atmospheric component has been proven difficult as it is spatiotemporally highly dynamic and a superposition of two atmospheric states. Here, we propose an approach to parameterize the stochastic properties of the single-epoch atmospheric delay field as a way to define the atmospheric signal. We found that the atmospheric signal of a time series of interferograms can be characterized by structure functions, which can be used to isolate the single-epoch structure functions. Due to the scaling properties of the atmospheric signal, it is then possible to construct a parametric function per SAR acquisition, using two isotropic and three anisotropic parameters. In particular, the isotropic parameters for the short-distance variation and long-distance variation in atmospheric delay can be used to characterize the atmospheric signal. For a test set of 151 Sentinel-1 acquisitions, this results in an atmospheric energy range of about 10 for short-distance scales and about 50 for long-distance scales. Our parameterization demonstrates that we can describe the spatiotemporal variability of InSAR atmospheric delays, which provides a measure for atmospheric noise for individual epochs in deformation time series based on distance and azimuth.

I. INTRODUCTION
A TMOSPHERIC disturbances are dominant contributors of noise in interferometric synthetic aperture radar (InSAR) deformation time-series analysis [1], [2]. Over the years, different algorithms have been developed to suppress this noise, including low-pass filtering in time and different types of spatial filters [3], [4], [5]. Beside the turbulent effect of the atmosphere, there is also an elevation-dependent effect due to stratification of the atmosphere [1]. To mitigate the stratification effect, the relation between elevation and atmospheric delays is used [6], [7].
Apart from estimating the disturbing signal from the InSAR data by means of filters and correlations with elevation, also numerous studies exist that correct the InSAR signal by means of weather models [8], [9], [10], [11], global navigation satellite system (GNSS) measurements [12], moderate resolution imaging spectroradiometer (MODIS), or medium resolution imaging spectrometer (MERIS) data [13]. Although these methods do show some improvement in the deformation estimates, they are limited in the further improvement of atmosphere correction. First, the resolution of the models is much lower than that of the radar data [14]. In some studies, this problem is partly solved by increasing the model resolution [15], but if these models are not backed up by accurate measurements, this does not yield reliable results. Second, the timing of the model realizations is often not correct, which leads to corrections of weather fronts and other large systems to be shifted [16]. Third, weather models themselves are not tuned on the exact location of convective atmospheric disturbances, but mainly whether they can predict that certain weather events will occur in a certain region. This makes the correction of the InSAR data with these models suboptimal.
In most contemporary persistent scatterer (PS)-InSAR techniques, statistical techniques are used to reduce the influence of the atmosphere. To quantify atmospheric noise in time, variance component estimation (VCE) is used [5], while the atmospheric delay per epoch is isolated using least-squares collocation [17] or estimation of the reference (previously denoted as "master") atmosphere first [4], [18]. However, these techniques do not offer a way to discern between spatially correlated tropospheric and deformation patterns, and removing the reference atmosphere only works for long time series and is hindered by slow steady-state deformation and decorrelation. Also, the quality of deformation estimates is often determined by a posteriori statistical analysis of residuals [5], which gives only quality information in retrospect. These methods are therefore difficult to use for early warning systems.
Here, we propose a statistical way to estimate and characterize the atmospheric signal in InSAR time series. By computing an interconnected subset of interferometric combinations, we find anisotropic turbulence information for single epochs, which characterizes the atmospheric disturbance dependent on distance and azimuth for every epoch in InSAR deformation time series. These estimated atmospheric disturbances can then be used to describe the expected noise and apply weights to the individual epochs in a deformation time series. Meanwhile, this offers a means to improve deformation time series that is based on physics and does not remove any nonatmospheric signals. We show how this approach can offer a physics-based characterization of the atmospheric state for individual images using only five parameters per single epoch. The approach involves: 1) the creation of an interconnected subset of interferometric combinations with short temporal baselines; 2) calculation of anisotropic second-order structure functions for all interferometric combinations; 3) estimation of single-epoch structure functions based on the results of the interferometric structure functions; and 4) parameterization of the isotropic and anisotropic structure functions for single-epoch SAR images.
The final parameterizations can then be used to describe InSAR tropospheric noise in, for example, deformation time series or change detection in early warning systems.

A. Parameterization of Atmospheric Signal in InSAR
The relation between the SAR signal delay and atmospheric parameters is defined by the integration of the atmospheric refractivity N from the scatterer on the ground to the satellite along the path s where d is the one-way atmospheric delay. The largest part of the absolute atmospheric delay is therefore related to the length of the integration path in (1), i.e., dependent on the topographic elevation of the scatterers. The refractivity is [19] where T is the temperature in kelvin, e is the partial pressure of water vapor, P d is the partial pressure of dry air, L is the liquid water content, n e is the electron density per cubic meter, f is the radar frequency, and W is the liquid water content. The values of the constants are k 1 = 77.6, k 2 = 71.6, k 3 = 3.75 × 10 5 , k 4 = −4.028 × 10 7 , and k 5 = 1.4 [20]. The first term in this equation represents the hydrostatic part and the second and third parts represent the wet part [21]. The fourth term represents the ionospheric delay, while the fifth term is the liquid water delay. From these, the ionospheric part can be modeled accurately [22] and the liquid water delay is negligible small [1]. The double-difference nature of the InSAR observations yields a sensitivity to the spatial variability of N , rather than its absolute value. Therefore, we need to focus on the spatial variability of P d , T , and e in (2). These three parameters are not uncorrelated. Most importantly, the maximum value of the water vapor pressure e is exponentially dependent on temperature [23] e max (T C ) = m 1 exp where T C is the temperature in • C and m 1 = 6.1094, m 2 = 17.625, and m 3 = 243.04 are constants. This means that although most delay variations in InSAR images are driven by water vapor fluctuations, the maximum amount of water vapor, and therefore the magnitude of the signal, is defined by temperature.
The delay variations of the tropospheric delay for a single SAR image can be described by a second-order structure function [24] where D d (l) is the structure function for distance l for the spatial variation of InSAR delay and p is the location of any pixel in the SAR image.

B. Scaling of Structure Functions
As with many other atmospheric measurements, InSAR delay shows a clear scaling behavior [1], [25], related to energy cascades of atmospheric turbulence and leading to the often observed 2/3 and 5/3 scaling of atmospheric data [26], [27]. Based on the turbulence theory, also predictions have been made for the scaling of tropospheric delay of microwave radiation [24], yet it is not undisputed whether the assumed scaling that is observed from small-scale wind measurements can be translated into similar scaling of InSAR delays that are dependent on refractivity. This theoretical scaling has also been compared with InSAR measurements [1], albeit for distances up to ∼50 km. Moreover, due to significant orbit errors in satellite missions such as ERS-1 and 2, InSAR images were always detrended, which led to the removal of spectral energy at large distances and therefore unreliable scaling values at these distances.
With new wide-swath SAR missions such as Sentinel-1 and more precise orbits, scaling of tropospheric InSAR delay up to 100 km can reliably be observed and investigated. Fig. 1 gives an example of a structure function derived from a Sentinel-1 interferogram. An average structure function is calculated based on the structure functions of 1450 interferograms from a time series of 151 Sentinel-1 SAR acquisitions over The Netherlands (see Fig. 2).
The curve(s) can be approximated by two exponential functions, one for local distances and one for regional distances. The location of the transition between these two functions is given by the transition point distance l t i ,t j T , where t i and t j are used to denote the interferometric combination at times t i and t j , respectively. This gives with the distance of the transition point l t i ,t j T , energy factors C t i ,t j s and C t i ,t j w , and scaling exponents ζ s and ζ w , for the local and regional scales, respectively. Because the scaling exponents are constant for local and regional scales, this results in two connecting straight lines in log-to-log space, with the energy factors C t i ,t j s and C t i ,t j w defining the location of the lines at 1 km, and scaling exponents ζ s and ζ w defining the slopes of the lines. Hence, Fig. 2 shows a different scaling for local and local distances. Note that in other studies, the energy factors are denoted as C 2 [1], [24], but we use C for brevity. Because the same local and regional scaling is consistently observed in the InSAR images in this study, they will be the core principle on which we base the parameterization of structure functions. The observed scaling exponents in Fig. 2 are ζ s = 0.67 ± 0.01 for small (i.e., local) distances of up to a number of kilometers, and ζ w = 1.34 ± 0.01 at large (i.e., regional) distances, which are equivalent to a 2/3 and 4/3 scaling, suggesting a connection with Kolmogorov turbulence theory. Note that the observed scaling exponents contrast with findings in other studies [24], but are consistent for many individual interferograms over different seasons. We will parameterize single-epoch structure functions instead of theoretical predictions because the resulting parameterized functions will resemble the observed structure functions much more closely. The observed scaling at local and regional scales will therefore be the main base for the proposed parameterization in this study.

A. InSAR Data
We use a stack of 151 SAR acquisitions and create a set of interferograms spanning a maximum lag of B t,m = 60 days. All interferograms are coregistered and resampled using the Sentinel-1 precise orbits, corrected for the topographic phase based on the SRTM digital elevation model (DEM) [28], and georeferenced. After projection using an oblique Mercator projection, square grid cells of 500 m are defined, in which the complex values are averaged (multilooking), to limit the computational cost and improve phase signal quality. Then, the geocoded and multilooked interferograms are unwrapped [29]. To prevent unwrapping errors, slow steady-state deformation, and decorrelation of large areas in our study, we limit the temporal baseline to B t,m = 60 days. As our focus is on tropospheric disturbances, acquisitions with large ionospheric trends and interferograms with evident extreme deformation patterns are removed from the analysis. Alternatively, ionospheric signals [22] could be corrected for. Fig. 1 gives an example of InSAR atmospheric delay over the region of interest. To derive zenith delays from the unwrapped phase data, we use where d t i ,t j is the one-way delay for interferogram (t i , t j ), φ t i ,t j is the unwrapped phase, λ is the radar wavelength, and θ is the elevation angle.
To create structure functions from these interferograms, we mask out water areas and use a coherence threshold of γ = 0.1. By differencing the unwrapped delay values over different distances and azimuths, the structure function per interferogram can then be obtained.

B. Calculation of Structure Functions
Based on the theory from Section II-A, we can now develop methods to calculate structure functions based on a time series of interferograms, obtained using methods from Section III-A. Individual structure function values for different distance and azimuth bins can be estimated by averaging found InSAR phase values for shifts in the east and north directions where D d (a, b) is the structure function for a pixels in the east direction x and b pixels in the north direction y. N x and N y are the number of pixels in the xand y-directions, respectively, and is the interferogram resolution, which is 500 m in our case. The distance l and azimuth α of each structure function value are based on the integer steps a and b. The variance of the structure function value D d (a, b) is given by By calculating D d (a, b) and σ 2 D d (a, b) for all a and b, we can cover the full range of distances l and azimuths α. Because the combination (−a, −b) will yield the same results as (a, b), we use only a half-space with the range of a, . All combinations of a and b are then mapped to pairs of distance l and azimuth α, creating a vector of distance and azimuth values,

C. Estimation of Structure Functions for Single Epochs
To characterize the atmospheric signal for individual epochs in an InSAR time series, we need to disentangle the combined interferometric structure functions into single-epoch structure functions. The structure function of a summation (or difference) of two independent atmospheric states is formed by the addition of the contributing variances [30]. Under the assumption that the InSAR signal only contains an atmospheric signal, we can use D where t i and t j are the two epochs of the interferogram. Due to this relationship, we can find the single-epoch values for a specific distance l using a set of interferometric combinations at that specific distance l. Therefore, the conversion to single-epoch structure function is done for all epochs at once, but separately per individual distance, l. After this conversion, the parameters of the single-epoch functions are estimated for all distances, but separately per time step t (see Section III-D). To estimate the structure functions per epoch, we create an interconnected subset of interferometric combinations with a maximum temporal lag of B t,m = 60 days. A system of observation equations is then formulated per individual distance l as Here, the N × 1 observation vector y, where N is the total number of interferometric combinations with maximum temporal baseline B t,m , contains the structure functions values D Matrix A defines all these interferometric combinations, i.e., and vector s contains the single-epoch structure function values at a distance l The used covariance matrix Q y l contains the variances σ t i ,t j D (l) for the interferometric combinations at distance l on the diagonal, see (7), and zeros at the off-diagonal elements.
After least-squares inversion, we obtain both the estimated single-epoch vectorŝ l and its variance-covariance matrix Qˆs l , which will be used as input for the parameter estimation procedure in Section III-D. Fig. 3 shows the results for two single-epoch structure functions compared with the original structure function for an interferogram. To evaluate the validity of the used model [see (12)], we perform an overall model test using the residualsê l = y l −ŷ l [31] A value ofσ 2 > 1 can be interpreted as an underestimation of the a priori variances, whileσ 2 < 1 can be interpreted as an overestimation the a priori variances. Using the test statistiĉ σ 2 , we then adjust Q y l to [5] Q cal y l =σ 2 Q y l .
Using least-squares inversion, we can subsequently obtain the calibrated variance-covariance matrix of the estimated single-epoch values Q cal s l . Accommodating for the potential anisotropy of the signal, the single-epoch structure function may need to become directional as where α is the azimuth direction of the structure function.
We can now estimate the directional single-epoch structure functions in the same way as derived in (12)-(15), resulting in a new parameter for s l,α for one particular distance l and azimuth α Fig. 1(c) gives the result of this directional structure function. Its orientation varies from north-south to east-west and back to north-south. Often these values can be related to image-wide trends in tropospheric delays. Fig. 1 shows anisotropy for distances larger than 5 km in contrast with a more isotropic behavior for smaller scales, which corresponds to the local and regional turbulence regimes. After least-squares inversion, we obtain the estimates of the anisotropic single-epoch structure functionsŝ l,α and its variance-covariance matrix Qˆs l,α . Combining all estimates over all distances l and azimuths α into one matrix where the columns s l,α represent all n single-epoch solutions per (l, α)-combination and the rows s t each single-epoch solution for all (l, α)-combinations. The rows ofŜ will be used to estimate the descriptive parameters per single-epoch function in the next section. To find the variance-covariance matrix Q s t for every individual s t vector, we use the variance for time step t, of all variance-covariance matrices Q cal s l of the (l, α)-combinations. These variances are then used as the diagonal elements of matrix Q s t and all off-diagonal values are set to zero.

D. Parameterization of Single-Epoch Estimates
As described in Section II-A, we can divide the atmospheric refractivity signal into three components, one related to the local turbulence regime, one to the regional atmospheric patterns, and one to DEM differences. As the DEM component can be corrected for, we focus here on the two other regimes by choosing a study area over The Netherlands, with minimal topography. The remaining local and regional regimes can both be characterized by one variable, C t w and C t s . To find those values for every isotropic solution, we use the structure function values s t for the single-epoch solutions obtained by methods described in Section III-C. The used isotropic parameterization of the structure function is where C t s is the parameter representing the strength of the local component, C t w is the parameter representing the regional component, and ζ s = 0.67 and ζ w = 1.34, their respective scaling parameters. Both components are combined by squaring the components and taking the square root the combined values, to get a transition between local and regional component that is close to the observed structure functions. In logarithmic space, this results in a line with a 0.67 slope for the local scaling and a 1.34 slope for the regional scaling (see Fig. 4). This line can then be fit using the estimated structure functions for single epochs s t as observation values. However, when the single-epoch values and variances are used directly in the cost function, the final estimation becomes very sensitive to outliers at larger distances due to the logarithmic scaling of the structure function. Therefore, the function will be fit by minimizing the residuals in log space where s t are the estimated single-epoch structure functions values, f t D ( p t ) is the parameterized structure function from (21), and Q By minimizing the residuals in (22), we obtain the vector of function parametersp t = [C t s , C t w ] T and its variance-covariance matrix Q loĝ p t Although the data fit for the isotropic solution can provide an estimate of the strength of the local and regional atmospheric signal, it does not provide any information on a potential direction-dependent, i.e., anisotropic part of the signal. To model the anisotropic component, we need to find the direction at which the structure function is strongest and whether and how the anisotropy develops over distance. The local turbulence signal appears to be almost isotropic, which The dotted blue, green, and red represent the local turbulence factor C t s and the regional turbulence factor C t w . Adding these two lines together results in the black line as described in (21). The transition point T t from the local to regional regime is given by the location where the red and blue dotted lines cross, which is at a distance l t T of about 8 km in this case.
is also clear from Fig. 1. Therefore, an anisotropic component can be added to (21) as To limit the number of used variables, this is done by adding three variables: the direction α t max , related to the strongest signal, the anisotropy scaling factor r t l , and by splitting the regional scaling parameter C t w in two parameters, one parameter in the x-direction C t w,max along azimuth α t max and one parameter in the y-direction C t w,min orthogonal to x. By solving the cost function given in (22), we then obtain the parameter vectorpp and its variance-covariance matrix in log space Q loĝ p t . Fig. 5 gives an example of the anisotropic single-epoch fit and the relative error between the estimated and fitted anisotropic single-epoch structure function. To clearly show the contribution of the anisotropic part of the structure function fit, also a decomposition in a separate isotropic and anisotropic component is given in Fig. 6.

E. Transition Point
To give an indication of the point where the local and regional components have equal strength, we introduced the transition point T t earlier (see Fig. 4). As the transition point T t is defined as the point where the isotropic local and regional component are the same, the distance of the transition point, l t T can be derived by C t w l t ζw T = C t s l t ζs The log-scale standard deviation of the estimated l t T is then obtained using error propagation of the variances and the covariance of C t s and C t w from the variance-covariance matrix Q loĝ p t .

F. Residuals of Single Epoch and Parameter Estimation
To analyze the residuals of the single epoch and parameter estimation, we calculate and combine the residuals similar to matrixŜ [see (20)] to estimate the residuals over specific time steps t, distances l, and azimuths α. The residuals for the derivation of single-epoch functionsÊ s and parameter estimationÊ p are given bŷ E s = ŷ l 1 ,α 1ŷ l 1 ,α 2 . . .ŷ l z ,α k − y l 1 ,α 1 y l 1 ,α 2 . . . y l z ,α k Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. where y l,α is the vector of interferometric structure function values for a distance and azimuth pair and s t is the vector of single-epoch structure function values for a specific epoch. By replacing all vectors s l,α forŝ l,α to derive y l,α , we can find the combined residuals for both the single epoch and the parameter estimation,Ê ĉ where s l,α is the vector of single-epoch structure function values for a distance and azimuth pair. These residuals can then be binned and or averaged for different azimuths α, distances l, and epochs t to find a root-mean-square error (RMSE). However, this value will be heavily skewed to the larger values because the size of the residuals is often a fraction of the structure function value. We will therefore use the root-meansquare relative error (RMSRE). For example, the RMSRE of the parameter estimation as a function of distance and azimuth is given by whereŝ t andÊ t p are the row vectors fromŜ andÊ p for specific time t, distance l, and azimuth α values. Similar to (30), the RMSRE over distance and azimuth for the single-epoch estimation RMSRE t s (l, α) and the combined single epoch and parameter estimation RMSRE t c (l, α) can be calculated. The summation can be done over time t but also over distance l or azimuth α, providing us estimates for RMSRE l,α p (t) and RMSRE t,α p (l).

IV. RESULTS AND DISCUSSION
In this section, we first discuss the results of the parameterization of the isotropic and anisotropic single-epoch structure functions for a time series of 151 Sentinel-1 images over The Netherlands. We discuss the transition point between the local and the regional component and the effect of phase noise due to deformation and decorrelation in InSAR data on the estimated structure functions. Finally, we discuss the residuals of the single epoch and parameter estimation. The used Sentinel-1 data stack in this study contains 1445 interferograms between early 2017 and late 2019, with a repeat pass of six-day and satellite overpass around 7 A.M. local time.

A. Isotropic Estimates for Small and Large Distances
For the isotropic parameterization, the results can be summarized in a time series of the local and the regional component, C t s and C t w , respectively. Fig. 7(a) shows these parameters over time, with a moving average filter of two months. From (5), it follows that these parameter values correspond to the local and regional component curves at 1 km. Both parameters show a yearly cycle, with high values in summer and low values in winter. Seasonal temperature variation affects the maximum water vapor pressure and, consequently, the integrated refractivity. Although temperature in itself has a negative linear correlation with refractivity [see (2)], the maximum water vapor pressure increases exponentially with increasing temperatures [see (3)], which causes the observed higher parameter values in summer than in winter. Fig. 7(a) also shows that the peak of the regional parameter lags behind the peak of the local parameter. This suggests that the local and regional variations are driven by different mechanisms. The energy for the local component peaks in early summer, which coincides with the yearly peak of incoming solar radiation that drives small-scale turbulent processes. These small-scale turbulence develop due to heating of the surface by solar radiation and cause variations in both temperature and water vapor, which affects the InSAR delay variations. The regional component peaks at the end of summer when the mean air temperatures are highest [32], which mainly affects the total delay differences due to large-scale weather patterns.
The error bars in Fig. 7(a) indicate the precision of the local and regional parameter estimates in log space, σ 2 C t s Fig. 7. Variation of the fit local parameter C t s and regional parameter C t w for single-epoch solutions. These values are equal to the local and regional components C t s l 0.67 and C t w l 1.34 at 1 km l = 1, see The blue and red lines show the two-month moving average. This shows that both variables follow a yearly cycle, with low values in winter and high values in summer, although the local parameters show a larger peak. (b) Distribution of local and regional parameters. This shows a maximum relative difference of 10 for local parameters and 50 for regional parameters. and σ 2 C t w , which indicates more certainty for local parameters than for regional parameters. This is partly due to the lower precision of structure function values for distances greater than 50 km, which is because the uneven distribution and strong decrease of available pixel pairs for distances close to the satellite swath width. This leaves only a short reliable interval to fit the regional component between the transition from local to regional scaling around 10 km and the maximum reliable distance of 50 km. Also, due to the longer distances, few atmospheric features with the same distance size are covered in one image, which makes the estimation less robust as it is more sensitive to specific cases. Therefore, if the regional parameter for a single epoch could not be reliably estimated, it has been excluded from the analysis. Using a maximum uncertainty interval of exp σ 2 p t 2 > 1.5, this led to exclusion of 20% of the cases, which are also removed in Fig. 7(a). Fig. 7(b) shows the total distribution of the local and regional parameter, which shows a smaller spread of the local parameters than for the regional parameters. Under the assumption that both are log-normally distributed, the respective logarithmic mean and standard deviation are 10 −5.5 m 2 and 0.60 for the local parameters and 10 −6.3 m 2 and 0.84 for the regional parameters, which is equivalent to the maximum relative difference of 10 for local parameters and 50 for regional parameters.
This parameter estimation assumes that the signal is stationary over the whole scene. Although, in most cases, the weather state will be similar over the entire scene and the values in Fig. 7 are representative. In specific cases, however, there can be nonstationarity, for example, if there is a weather front present with different weather types at either side of the weather front. This spatial variation can cause variation as large as the expected seasonal variations, i.e., the vertical bandwidth in Fig. 7(a), leading to a local overestimation or underestimation. Yet, the estimated value for these cases is still valuable as an estimate of the average atmospheric signal for the whole scene.

B. Anisotropic Solution
The level of anisotropy a t (l) at distance l for every structure function is defined based on (24) where a t (l) is the anisotropic factor as a function of distance l. The magnitude of the anisotropic factor is defined by the strength of the scaling parameter r t l , which describes the difference in scaling between the azimuth direction with the highest variability, i.e., α t max , and its orthogonal complement, and the difference in strength between the two directions C t w,max and C t w,min . Cases with a dominant anisotropic signal are usually connected with a strong delay gradient over the original InSAR image. In these cases, the regional regime is relatively dominant compared to the local regime. Fig. 8(a) shows the relationship between the level of anisotropy and the location of transition point l t T , which scales with the ratio of the regional and local parameters C t s and C t w [see (26)]. The level of anisotropy is shown at 50 km because the anisotropy is more distinct at larger distances, and until 50 km, the anisotropic single-epoch structure function can still be reliably estimated. The level of anisotropy at 50 km ranges between 1, which means no anisotropy, to about 16, with most cases around a value of 2. Values with low accuracies are those where the regional component of the structure function is not well defined (see Fig. 7), which also leads to lower accuracies for the anisotropy values. Fig. 8(b) shows the distribution of the distances of the transition point l t T between the regional and local regimes. High values indicate a dominant local regime, while low values indicate a dominant regional regime. A dominant local regime is related to a very turbulent troposphere, i.e., strong convective processes. A dominant regional regime indicates wide-scale delay trends, i.e., related to passing weather fronts. T , which scales with the ratio of the local C t s and regional C t w parameters [see (26) and Fig. 4]. The level of anisotropy is given as the ratio between the direction with the highest and lowest structure function value at 50 km [see (31)]. The correlation between the level of anisotropy and location of the transition point indicates that when large-distance weather systems are dominant, the level of anisotropy for these systems also increases. (b) Distribution of the ratio of the local C t s and regional C t w parameter. The standard deviation per epoch σ T,t is represented by the shade of blue, with darker blue for lower standard deviation. This shows that although there is a large distribution of l t T , the high accuracy estimation is clustered around the 5-10-km range.

C. Transition Point
Because the precision of the single-epoch structure function deteriorates at distances over 50 km due to the InSAR image size, the estimated distances of the transition points l t T become less precise if they approach 50 km. This can be seen from Fig. 8(a), but also Fig. 8(a) where transition point l t T values with low precision are indicated in light blue. Fig. 8(a) shows that there is a clear relationship between the level of anisotropy and the location of the transition point, with a smaller distance of the transition point with increasing anisotropy. This trend indicates that when the regional component is more dominant, which results in a lower estimate of l t T , the level of anisotropy is also larger. This is possibly related to cases with a less turbulent character of the delay, where the delays are mainly caused by large-scale trends in pressure and temperature, which leads to a strong anisotropic component.
The position of the transition point is often related to the effective height of the wet troposphere [24], [33], which would lead to different scaling behavior between smaller atmospheric processes that are not limited by the height of the atmosphere and larger processes that are limited by the height of the atmosphere [34]. However, the estimated values of l t T of tenths of kilometers are too high to represent the effective height of the wet troposphere, and also, the range between the lowest and highest values of l t T indicates that a direct conversion between the transition point and the effective height of the wet troposphere cannot be made. Interestingly, there are also no gaps between local and regional curves, which indicates that both processes are present at all scales, but that one is dominant at local distances and the other at regional distances. This is also implicated in (21) and (24) and can be seen from the compounding effect around the transition point (see Fig. 4). Similar scaling values for local and regional distances and location of the transition points are found for wind measurements [35], [36], [37], which could represent the same processes as we see in InSAR data although the underlying measurement is very different. Also, another explanation for the observed scaling behavior has been hypothesized [38], [39], [40], but how these would result in similar scaling in InSAR data is unclear.

D. Bias in Estimated Structure Functions
In addition to the atmospheric components, a part of the InSAR delay variations consists of deformation and decorrelation noise, which is only partly mitigated by masking out highly decorrelated areas and removal of interferograms with evident extreme deformation patterns. Because we expected that the noise term would increase with temporal baseline B t , only interferograms with a maximum temporal baseline of 60 days were included in the analysis, to limit the noise in our final estimate (see Section III-C). To quantify the increase in noise with temporal baseline, we analyze the results for the single-epoch estimation for varying temporal baselines.
The noise increase with larger temporal baselines can be shown by the ratio between the interferometric structure functions and the adjusted interferometric structure functions y l /ŷ l for specific temporal baselines B t [see (13)]. The resulting ratios are given in Fig. 9, which shows the structure function values of the observations y l=5 and adjusted observationsŷ l=5 at 5 km distance as a function of the date of the primary SAR acquisition and the temporal baseline. Fig. 9(c) shows a clear increase in this ratio with larger temporal baselines, which indicates increasing noise values for larger temporal baselines. There is also a clear correlation between larger ratios, shown in red in Fig. 9(c) and low structure function values, which are blue in Fig. 9(a) and (b). This is likely because the relative contribution of deformation and decorrelation noise becomes larger with less atmospheric turbulence.
So far, Fig. 9 has shown the ratios for 5 km, but these are expected to decrease with distance as phase variations due to temporal decorrelation and deformation are generally only correlated over short distances, while the atmospheric delays scale with distance. Fig. 10 shows the dependence of the noise parameter on distance. Similar to Fig. 9, we show the ratios y l /ŷ l , but now for a specific temporal baseline t Because the absolute values vary a lot between seasons, we use a ratio instead of absolute values. This shows that as the temporal baseline increases, the structure function results become more contaminated with nonatmospheric noise due to decorrelation and deformation. The primary date of the interferogram is the most recent date of the two dates of an interferogram.
instead of a specific distance y l,B t = t /ŷ l,B t = t . This gives a clear view on the average structure function values per distance for the different temporal baselines B t = [6, 12, . . . , B t,60 ]. This indicates that structure function values at B t = 60 days are more than 50% larger than those at B t = 6 days for local distances. As expected, the noise contributions decrease with distance. Due to the used weights for the single-epoch estimation, the average noise in the adjusted interferograms is about the same as the measured 18-day interferograms. This means that in estimated single-epoch solutions, there is a larger noise contribution than in the structure functions for six-and 12-day interferograms.
To minimize the noise contributions in the single-epoch estimation, we could therefore increase the weight for structure functions with short temporal baselines, but this could cause problems in case one SAR acquisition has a bad quality, due to, for example, decorrelation caused by snow cover. Also, the improvement will be limited to local distances, while the new solutions for regional distances will degrade as they are much less affected by noise and provide better solutions with more equal weights for different temporal baselines. Another solution would be to use a network of PSs instead of the used spatial averaging technique. However, this would cause a very uneven sampling, with most data points in build up areas, resulting in a bias of the estimated structure functions. Different weighting techniques could therefore help to mitigate decorrelation and deformation noise but would introduce other errors in the single-epoch estimates.
The change in grid cell size and maximum temporal baseline will also reduce the bias, which becomes clear from Fig. 10. If the maximum temporal baseline is decreased, the more biased interferograms, shown as the top lines in Fig. 10, are removed from the analysis, resulting in a reduced bias in the single-epoch estimate. Similarly, a larger grid cell will cut off the x-axis in Fig. 10, removing the distances that are mostly affected by biases. Fig. 10. Comparison of the relative residuals y l,Bt = t /ŷ l,Bt = t for different temporal baselines t. This shows that relative residuals increase with temporal baseline, especially at smaller distances, which indicates that the relative contribution of nontropospheric noise is larger with larger temporal baselines and smaller distances. The increase with shorter distances is largely due to the smaller tropospheric signal, while the nontropospheric noise remains constant as its spatial correlation is small. Due to the used weighting techniques, the noise level of the final solution averages around the value for 18-day interferograms.

E. Residuals of Structure Function Parameterization
To quantify how well the estimated isotropic and anisotropic models can represent the statistical characteristics of the calculated structure functions and estimated single-epoch structure functions, the RMSRE is used (see Section III-F). The RMSRE gives the average relative deviation of the model from the observed value as a fraction, for different epochs, distances, or directions (see Fig. 11). Most found RMSRE values are between 0.1 and 0.3, which corresponds to a relative error of 10%-30%, but there are some cases with larger residuals. Fig. 11(a), (d), and (g) shows the average RMSRE over distance; Fig. 11 The residuals over distance for the single-epoch estimates in Fig. 11(a) show that for almost every distance bin, the anisotropic solution shows a better performance. This difference is smallest for distances lower than ∼5 km, as the structure functions are much more isotropic over these distances. The anisotropic curve in Fig. 11(a), (d), and (g) first decreases until about 5-10 km and increases again at distances larger than ∼10 km. The higher residuals at small distances are mainly due to biases caused by deformation and phase noise, as discussed in Section IV-D. This affects the single-epoch estimation because (11) only holds when the temporal baseline does not affect the structure function of an interferogram, which does hold for the tropospheric delays but not for deformation and phase noise.
The increase of RMSRE values for distances over ∼10 km is likely also because the isotropic condition in (11) is not met. This is because the structure function becomes more anisotropic at these distances, which is also reflected in the increase in RMSRE of the isotropic relative to the anisotropic solution in Fig. 11(a), (d), and (g). The further increase around 50 km for the parameter estimation in Fig. 11(d) is due to This shows that residuals for the anisotropic solution are generally better, but mainly at higher distances. Lowest residuals are found at the 5-10-km range for the anisotropic case, with about two times larger residuals at the smallest distances and about five larger higher residuals at the largest distances. The overall mean RMSRE is 0.18 for the anisotropic and 0.31 for the isotropic case. the uneven distribution and strong decrease of the number of available pixel pairs for distances close to the satellite swath width. Fig. 11(b), (e), and (h) shows that isotropic residuals are higher than the anisotropic ones, the latter having an RMSRE of lower than 0.2 for most individual epochs. The distribution for the isotropic case is also more skewed, because for cases with a high level of anisotropy, the precision of the isotropic parameter estimation decreases. Fig. 11(c) shows the average RMSRE in both distance and azimuth, which shows largely the same pattern as the anisotropic curve in Fig. 11(a), but shows a lower precision around the dominant west-south-west wind direction for larger distances, which is likely due to stronger anisotropic InSAR delay trends in this direction. Fig. 11(g), (h), and (i) shows the combined error of both the single epoch and parameter estimates. Interestingly, the RMSRE for the combined residuals in Fig. 11(g) is for some distances lower than the RMSRE distribution for the parameter estimation in Fig. 11(a) and (d). This is due to a "correction" of low-precision single-epoch estimations at large distances in the parameter estimation. The overall mean RMSRE value derived from the distribution in Fig. 11(h) is 0.31, with a standard deviation of 0.15 for the isotropic case and 0.18 with a standard deviation of 0.07 for the anisotropic case.

V. CONCLUSION
In this article, we developed a new method to estimate and parameterize single-epoch structure functions, using the scaling properties over local and regional distances. Using this method, we showed that we can accurately describe the strength and statistical characteristics of the tropospheric signal for single-epoch SAR acquisitions, using only two isotropic and three anisotropic parameters.
To derive these results, the structure functions of an interconnected subset of interferograms need to be calculated from unwrapped InSAR interferograms and used to derive single-epoch structure functions. These single-epoch structure function are then parameterized using an isotropic and anisotropic function fit in log space.
The results show a clear distinction between a local scaling up to 3-30 km and a regional scaling for larger distances with increasing anisotropy. Both in local and regional scaling, a clear seasonal cycle was observed due to temperature and water vapor variations, with peak values in early summer for local and late summer for regional distances. In our study case over The Netherlands, this resulted in a maximum variation with a factor of ∼10 for the local and ∼50 for the regional parameters during the study period.
The resulting single-epoch structure functions can be used as a physics-based method to quantify the atmospheric signal in the InSAR time series. This can then be used as a precision and weighting factor in model estimation for InSAR deformation time series, as it provides a measure for atmospheric noise for every individual SAR scene based on distance and azimuth. In addition, the gained insights can be used to assimilate InSAR data in weather models as a meteorological measurement.