A Covariance-Based Approach to Merging InSAR and GNSS Displacement Rate Measurements

Abstract

This paper deals with the integration of deformation rates derived from Synthetic Aperture Radar Interferometry (InSAR) and Global Navigation Satellite System (GNSS) data. The proposed approach relies on knowledge of the variance/covariance of both InSAR and GNSS measurements so that they may be combined accounting for the spectral properties of their errors, hence preserving all spatial frequencies of the deformation detected by the two techniques. The variance/covariance description of the output product is also provided. A performance analysis is carried out on realistic simulated scenarios in order to show the boundaries of the technique. The proposed approach is finally applied to real data. Five Sentinel-1A/B stacks acquired over two different areas of interest are processed and discussed. The first example is a merged deformation map of the northern part of the Netherlands for both ascending and descending geometries. The second example shows the deformation at the junction between the North and East Anatolian Fault using three consecutive descending stacks.

1. Introduction and Motivation

Synthetic Aperture Radar Interferometry (InSAR) deformation rate measurements support a wide range of applications in the fields of geology, geophysics, and geohazards [1]. Since the nature of the interferometric measurements is inherently relative, the additive delays—like those of the troposphere and ionosphere—make the accuracy of such measurements strongly dependent on the distance [2,3,4]. Hence, InSAR performance varies according to the particular application. Typically, interferometry works well for applications such as infrastructure monitoring or urban subsidence, since they involve relatively short scales (10–20 km). However, there is also scientific interest in using InSAR to measure the strain accumulation of tectonic faults. This type of study requires very accurate measurements (1 mm/year at distances larger than 100 km) [5].

The new generation of Synthetic Aperture Radar (SAR) sensors—like Sentinel-1—provides systematically acquired data with a swath width of 250 km [6]; future missions plan to further extend this to 350 km [7]. Due to the atmospheric errors at such scales, the requirement of 1 mm/y could be a challenging goal for InSAR, especially if the available time series have a reduced observation time [8]. Therefore, in order to fully exploit the coverage capabilities of these missions, we are motivated to develop techniques that merge interferometric deformation rates with other geodetic measurements, not only to improve the performance of applications particularly affected by atmospheric effects (e.g., inter-seismic deformations), but also to remove the reference point, hence facilitating an easier integration with the other geometries or techniques.

The best candidate for this integration with SAR interferometry are the deformation rates estimated using a Global Navigation Satellite System (GNSS). This is due firstly to their well-known complementarity where InSAR coverage is combined with GNSS accuracy. Indeed, there is a long tradition of research in this area as evidenced by studies performed by geo-scientists; for example, in supporting interferometric phase unwrapping [9,10]. Combining InSAR and GNSS deformation rates has also been studied in order to calibrate the InSAR-estimated velocities [11,12,13], and to help the limited geometric sensitivity of the SAR system in measuring deformation [14,15,16]. The integration of InSAR with other instrumentation such as active transponders has been investigated in order to design optimal geodetic networks that enhance the spatial distribution of the deformation measurements [17,18,19].

However, it is important to highlight that the main argument of this work is not to claim the complementary InSAR/GNSS—which is well known, as previously mentioned—but rather to discuss the use of their error statistics in the combination, as well as its propagation through to the final product. In [8], the effect of the correction of the interferometric phase for tropospheric delays and solid earth tides using external models was shown. After such corrections, the mean variograms of the interferometric measurements error show a stationary behavior that can be well approximated by a covariance function. An optimal combination of the two measurement techniques is hence possible, since the spectral properties of their own errors can be taken into account. Moreover, no assumptions about the displacement pattern characteristics [19] and GNSS measurement density are necessary as in [11] and the data need not be filtered. Only the InSAR/GNSS differences, which can be statistically characterized, are estimated and used to compensate the original InSAR measurements. The error can then be propagated from the data to the results, allowing full characterization of the output uncertainty.

The measurement of strain accumulation along hundreds of kilometers is a potential application for this framework. The integration of InSAR and GNSS does not remove any deformation component, but performs a weighted merging of the different spatial frequencies of the deformation detected by the two techniques according to their reliability. Moreover, the mathematical modeling, often applied on final measurements, requires a proper weighting of the different input data [20]. Other possible applications include the InSAR-based National Ground Motion Services [21]. Such projects are often required to provide a product that merges the InSAR-derived results and the results derived by the GNSS networks already deployed on the territory. Since these products are part of a Service, a consistent description and traceability of the uncertainties is strictly required.

In Section 2 the proposed methodology is described, separating the capability of retrieving the absolute motion from the calibration of the residual atmospheric errors. An analytical description of the error of the merged product is also provided. In Section 3, simulations are performed in order to provide reference numbers in terms of coverage and quality versus performance. Finally, the results of using real data in significant test cases are presented and discussed. The first example covers the northern part of the Netherlands and demonstrates merging on a large scale in order to provide a consistent product like those in national ground motion services. The second example is aimed at addressing geophysical applications, and shows the calibration of interferometric measurements over North Anatolian Fault to the Eurasian plate.

2. Methodology

Let us consider N locations in the processed area of interest, where two deformation rate measurements $v_D$ and ${\underline{v}}_G$ are performed using InSAR and GNSS, respectively. With $\underline{s}$ as the radar line of sight (LoS) it is possible to state the problem modeling the difference between the two velocities at the ith position ${\mathsf{\Delta }}_i$ as follows:

$${\mathsf{\Delta }}_i=v_{D,i}-{\underline{v}}_{G,i}^T\underline{s}=\delta (r_i,a_i)+v_{ref}+n(r_i,a_i),$$

(1)

where $v_{ref}$ is the velocity of the reference point used in InSAR processing, $\delta (r_i,a_i)$ is a space-variant error screen in the interferometric data, and $n(r_i,a_i)$ is the random error. The implicit hypothesis behind this definition is that the spatially correlated noise, basically related to the residual atmospheric delay, is present in the interferometric measurement only. This characteristic allows extraction, through the subtraction of the two velocities, of the parameters $\delta (r,a)$ and $v_{ref}$, which have to be estimated in order to properly merge the data. It is worth noticing that this operation can also be seen as a calibration of the InSAR data: It exploits GNSS reference rates so that the interferometric measurement is no longer relative to the reference point while also removing residual systematic effects. Figure 1 displays an overview of the steps of the proposed approach. In the first step, the error statistics are derived and the covariance matrix $\mathbf{R}$ is computed. The reference point velocity $v_{ref}$ is then estimated and subtracted from the offset vector $\underline{\mathsf{\Delta }}$. Finally, the residual error screen is interpolated at each InSAR measurement point, exploiting knowledge of the error statistics.

2.1. Error Description of the Input Data

A proper handling of the relation defined in Equation (1) requires a statistical description of the vector $\underline{\mathsf{\Delta }}$. Therefore, the covariance matrix $\mathbf{R}$ of the difference measurements has to be derived. It is possible to distinguish between three contributions: The random noise of the GNSS measurements ${\sigma }_{G,i}^2$ obtained by projecting the variance of the different components of the estimated rates onto the radar LoS, the random noise of the InSAR measurements ${\sigma }_{D,i}^2$, and spatially correlated noise due to residual atmospheric effects. The full covariance matrix can be defined as:

$$\begin{array}{cc}\hfill \mathbf{R}=\mathbb{E}\left[\underline{\mathsf{\Delta }}{\underline{\mathsf{\Delta }}}^T\right]& =diag\left({\sigma }_{G,0}^2,\dots {\sigma }_{G,N-1}^2\right)+\hfill \\ \hfill & +diag\left({\sigma }_{D,0}^2,\dots {\sigma }_{D,N-1}^2\right)+{\mathbf{C}}_{\alpha }.\hfill \end{array}$$

(2)

The first two contributions are diagonal matrices, since they represent the spatially uncorrelated errors of the independent GNSS and InSAR measurement processes, respectively. On the other hand, the derivation of the covariance matrix ${\mathbf{C}}_{\alpha }$ that represents the residual atmospheric error is particularly interesting [22]. A covariance function $\mathsf{\Gamma }\left(d\right)$ representing the residual error has to be robustly estimated and evaluated at the N positions where the ${\mathsf{\Delta }}_i$ are located; with $d_{i,j}$ as the distance between the ith and the jth measurement, this is:

$${\mathbf{C}}_{\alpha }(i,j)=\mathsf{\Gamma }\left(d_{i,j}\right).$$

(3)

The first two components of Equation (2) can be estimated using the supplied accuracy of the GNSS data and the interferometric temporal coherences, respectively. However, estimation of ${\mathbf{C}}_{\alpha }$ requires a more detailed discussion.

2.2. Estimation of the Residual Atmospheric Error Covariance Matrix

The covariance of the residual atmospheric error can be estimated from interferograms where the residual atmospheric delay is the dominant measured delay. The fast revisit time of the Sentinel-1 mission allows the computation of short time interferograms with temporal baselines $\delta t$ down to 6 days. Following the approach in [8], the short temporal baseline interferograms are used to compute variograms merely representing the residual atmospheric errors. Since such errors are expected to be on the order of magnitude of a centimeter, deformation rates of several tens of cm/y are necessary in order to bias the estimation of atmospheric errors. Such rates can be reached in landslides or mining areas that are typically restricted in coverage. Since the variograms are computed by averaging many different measurements over the whole scene, the effect of such areas should not strongly impact the estimation. It should be noted that since variogram estimation requires the unwrapped phase, this step must be performed at the end of the processing. At this stage, areas of very high deformation are apparent and can eventually be masked out. The effect of solid earth tides has been corrected in the interferometric phase using models [8]. However, the projection of tectonic plate motion onto the line of sight is a large-scale effect whose spatial characteristics cause a bias in variogram estimation. Such movements are mainly horizontal and can reach 6–7 cm/y. How large the horizontal motion must be in order to be comparable with the troposphere in a short temporal baseline interferogram can be easily estimated. For a 1 cm gradient along the radar swath, a horizontal motion of $1/(sin\left({\theta }_{near}\right)-sin\left({\theta }_{far}\right))$ is necessary where ${\theta }_{near}$ is the incidence angle at near range and ${\theta }_{far}$ the incidence angle at far range, leading to a horizontal motion of $\approx 5$ cm or 152 cm/y for a revisit time of 12 days. Therefore, keeping the maximum $\delta t<30$ days should result in negligible impact. It should be mentioned that the eventual presence of seismic events in the time series should also be assessed, and co-seismic interferometric pairs not be used to generate the variograms.

Given that such interferograms are almost deformation-free, it is possible to assume that the average of the variograms $\mathbb{E}\left[V\right]$ is a good estimator $\widehat{V}$ (the $\widehat{}$ symbol indicates an estimated parameter) of the covariance characteristics of the residual atmospheric delays. Since velocity estimation implies a linear regression on the phase, a scaling factor accounting for acquisitions’ time span and number must be applied to convert the single-phase measurement accuracy into deformation rate accuracy:

$${\widehat{V}}_{rate}\left(d\right)=\frac{{\lambda }^2}{16{\pi }^2}\frac{\widehat{V}\left(d\right)}{2}\frac{M}{M{\sum }_k{t}_k^2-{\left({\sum }_k{t}_k\right)}^2},$$

(4)

where d is the distance, $\lambda$ the radar wavelength, $t_k$ the acquisition times, and M the number of interferograms used for the linear regression. The analysis in [8] showed that if atmospheric phase corrections based on European Centre for Medium-Range Weather Forecasts (ECMWF) models are performed [23,24], the residual atmospheric effects after processing are well-modeled as stationary.

${\widehat{V}}_{rate}\left(d\right)$ can then be fitted using a covariance model in order to compute the model parameters, converted to $\widehat{\mathsf{\Gamma }}\left(d\right)$, and to estimate the covariance matrix ${\mathbf{C}}_{\alpha }$ in Equation (3). In this study, the exponential covariance model was used, as seen in Figure 2.

2.3. Estimation of the Reference Point Motion

Equation (1) shows how $\underline{\mathsf{\Delta }}$ can be seen as an observation of $v_{ref}$ with superimposed random and colored noise. Since the noise statistics of $\underline{\mathsf{\Delta }}$ exhibit stationary behavior, the derivation of the reference point velocity can be carried out by averaging the observed $\underline{\mathsf{\Delta }}$. Therefore, according to the considerations in the previous Section 2.1, ${\widehat{v}}_{ref}$ can be found by taking into account the computed covariance matrix $\mathbf{R}$:

$${\widehat{v}}_{ref}={\left(\underline{u}{\mathbf{R}}^{-1}{\underline{u}}^T\right)}^{-1}{\underline{u}}^T{\mathbf{R}}^{-1}\underline{\mathsf{\Delta }},$$

(5)

where $\underline{u}$ is a unitary vector that maps, one to one, the measurements $\underline{\mathsf{\Delta }}$ with the unknown scalar $v_{ref}$. The estimated ${\widehat{v}}_{ref}$ represents the motion of the reference point in its local LoS and has to be added to the InSAR velocities in order to make them absolute.

2.4. Estimation of the 0-Mean Calibration Screen

In the previous Section 2.3, the overall offset between InSAR and GNSS deformation rates was estimated while accounting for the covariance matrix. Now, the space-variant error screen between the GNSS and InSAR velocities can be estimated by performing a covariance-based interpolation (Kriging) of the residual offsets ${\underline{\mathsf{\Delta }}}_0=\underline{\mathsf{\Delta }}-{\widehat{v}}_{ref}$. A set of coefficients $\underline{c}$ that, combined with the vector ${\underline{\mathsf{\Delta }}}_0$, allows the reconstruction of $\delta (r,a)$ everywhere must be estimated.

$$\widehat{\delta }(r,a)={\underline{c}}^T{\underline{\mathsf{\Delta }}}_0$$

(6)

According to theory, this can be obtained by imposing the condition that the interpolation/prediction error be uncorrelated:

$${ϵ}_n=\widehat{\delta }(r_n,a_n)-\delta (r_n,a_n)$$

(7)

with the data $\underline{{\mathsf{\Delta }}_0}$ [25]. Substituting Equation (6) into Equation (7) and imposing the uncorrelatedness condition with $\underline{{\mathsf{\Delta }}_0}$ gives:

$$\underline{c}={\mathbf{R}}^{-1}\underline{\rho },$$

(8)

where $\underline{\rho }=\mathbb{E}\left[\delta (r,a){\underline{\mathsf{\Delta }}}_0\right]$ is the vector representing the correlation between the data vector ${\underline{\mathsf{\Delta }}}_0$ and the error screen $\delta$ at the current position.

2.5. Variance and Covariance of the Results

To enhance usability, the error of the final product should be characterized. The final result is obtained by compensating the InSAR velocities for the estimated ${\widehat{v}}_{ref}$ and $\widehat{\delta }(r,a)$. The reference point deformation rate is a bias added to all points. Its error is therefore a constant value in the final product. Its contribution can be derived by computing the variance of the linear system inversion:

$${\sigma }_{{\widehat{v}}_{ref}}^2={\left(\underline{u}{\mathbf{R}}^{-1}{\underline{u}}^T\right)}^{-1}.$$

(9)

The error of the estimated $\widehat{\delta }$ will be space variant. The variance and covariance of the estimated error screen $\widehat{\delta }(r,a)$ can be derived as in [25].

$$\mathbb{E}\left[{ϵ}_n^2\right]=\mathbb{E}\left[{(\delta (r_n,a_n)-{\underline{c}}^T{\underline{\mathsf{\Delta }}}_0)}^2\right]=\mathsf{\Gamma }\left(0\right)-{\underline{{\rho }_n}}^T{\mathbf{R}}^{-1}\underline{{\rho }_n}$$

(10)

where $\mathsf{\Gamma }\left(0\right)$ represents the original error due to the spatially correlated signal, and the second part of the equation represents a “mitigation factor” that reduces such variance according to the distance from the GNSS data. Analogously, the covariance is:

$$\mathbb{E}\left[{ϵ}_n{ϵ}_m\right]=\mathsf{\Gamma }\left(d_{m,n}\right)-{\underline{{\rho }_m}}^T{\mathbf{R}}^{-1}\underline{{\rho }_n}.$$

(11)

Some conclusions concerning the final product can now be drawn. It provides absolute measurements that can be characterized by a variance directly derived from the applied methodology. It still contains a residual spatial correlation, though mitigated by the removal of $\widehat{\delta }(r,a)$. However, this residual error is no longer stationary. Its covariance depends on the two considered points, as highlighted by Equation (11). The variogram or covariance representation of the error—as for the input InSAR velocities—is hence no longer possible, but must must to be computed locally, accounting for the position $(r,a)$ with respect to the GNSS stations.

3. Simulations

Simulations were performed in order to assess the validity of the method and evaluate its performance. The scope of these simulations is to understand and quantify the effects of the spatial distribution (density) and the quality of the GNSS measurements in different processing scenarios, and to provide some numbers that summarize the achievable accuracies in the retrieval of ${\widehat{v}}_{ref}$ and $\widehat{\delta }$.

A set of randomly distributed points with zero deformation was simulated over a surface of $175\times 250$ km${}^2$ (Sentinel-1 slice size) in order to represent the InSAR measurements. Random noise (clutter) and spatially correlated noise (atmospheric residuals) were added to the points. For the generation of spatially correlated noise, an exponential covariance $\mathsf{\Gamma }\left(d\right)={\sigma }_{\alpha }^{2}exp(-d/L_c)$ was used. A reduced set of GNSS zero-deformation data were also simulated at random positions within the scene. According to the model only random noise was added to the simulated GNSS velocities. For the sake of simplicity, positions were uniformly distributed within the scene. The reconstruction depends on how well the available samples are able to represent the error spectrum. In practice, the more high-pass the error, the more samples will be needed. Of course, a regular sampling is desirable. If the data are concentrated in an area, the re-construction will be good in this area. When moving away from data points, the estimator will extrapolate.

The simulations were performed while varying the two main parameters, the number of reference GNSS measurements, and their accuracy. Two different scenarios were tested with low and high atmospheric residual power ${\sigma }_{\alpha }^2$. In the first scenario, ${\sigma }_{\alpha }^2=2$${\mathrm{mm}}^2/{\mathrm{y}}^2$; for the Sentinel-1 mission, this is comparable to having a long time series of $T_{obs}\approx 3$ y in the case of applied atmospheric corrections [8]. In the second scenario, ${\sigma }_{\alpha }^2=9$${\mathrm{mm}}^2/{\mathrm{y}}^2$, representing the case of a short time series of $T_{obs}\approx 1$ y [8]. A realistic correlation length, according to real data, was used for the simulations—$L_{corr}=60$ km in both cases. This comes from direct experience with the data. We observed that the average variograms after tropospheric corrections exhibited values of 40–100 km when fitted with an exponential model. An example of the simulation framework is shown in Figure 3.

Observing the simulation results in Figure 4 and Figure 5 it can be concluded that, due to the spatial correlation of the error, higher densities of GNSS measurements improves the estimation of $\delta$ since it provides a better sampling of the error field, but does not help much in retrieving the absolute velocity $v_{ref}$, due to the “data redundancy” introduced by the spatial covariance.

3.1. Performance in Retrieving ${\widehat{v}}_{ref}$

The accuracy of ${\widehat{v}}_{ref}$ is that of the mean computed from a set of correlated samples. The level of dependence of the given dataset is controlled by the spatial density of the latter and the correlation length $L_{corr}$ of the superimposed noise. We observed that it is logical to expect that, given a $L_{corr}$, increasing the number of data improves estimation of ${\widehat{v}}_{ref}$ up to a certain level, since as the density of measurements increases, so to does the amount of correlation between measurement, hence limiting its impact on final performance. In a typical scenario ($T_{obs}\approx$ 3 y, ${\sigma }_{\alpha }^2=$2 mm${}^2$/y${}^2$), the accuracy can easily be brought below 1 mm/y, even with a limited set of GNSS stations (Figure 4).

3.2. Performance in Retrieving $\widehat{\delta }$

Evaluating the performance of $\widehat{\delta }$ is not a simple task. The most natural way would be to compare the variograms of the results before and after application of the technique. This approach would deliver a deeper insight into the achieved gain, since it would be possible to show this as a function of the scale. However, as mentioned in Section 2.5, the error of the merged results is no longer stationary making a variogram representation questionable. In order to avoid this issue, a simpler but more robust approach was followed. The performances were evaluated in terms of mean square error (MSE) in $dB$, over the whole scene. In practice, the mean power of the residual deformation signal after removal of the estimated error screen $\widehat{\delta }$ from the measured rates $v_{meas}$ is:

$$MS{E}_{dB}=10Lo{g}_{10}\left(\mathbb{E}\left[{(v_{meas}-\widehat{\delta })}^2\right]\right).$$

(12)

The results displayed in Figure 5a,b show that, in both cases, the performance cannot improve beyond a certain level by improving the quality of the GNSS measurements only. An improvement in GNSS coverage is also necessary in order to better compensate the higher wave-numbers of the error screen.

4. Results

For the study, two Sentinel-1A/B datasets covering the northern part of the Netherlands (two stacks, ascending and descending) and the junction between the North Anatolian Fault and East Anatolian Fault (three stacks, descending) were used. The interferograms were computed, corrected for tropospheric delays using ECMWF ERA-5 (ECMWF Re-Analysis) data [8] and processed using the PSInSAR technique [26]. In order to preserve all wave-numbers of the deformation signal, no spatial high-pass filtering or polynomial detrending was performed on the final data. The technique described above was applied to the estimated deformation rates. As visible in Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11, the method is able to retrieve the absolute motion while also mitigating the undesired residual atmospheric error, hence displaying all available spectral components in the deformation signal. Figure 7a, Figure 9a, and Figure 11a show that the residual zero-mean error ${\mathsf{\Delta }}_0$ is quite small, varying between $\pm 2$ mm/y. This can already be considered as a kind of validation of the InSAR data, whose accuracy is close to 1 mm/y.

For the reference measurements, the GNSS data processed by Nevada Geodetic Laboratories [27,28] were used. The correspondence GNSS/PSs for the calculation of vector $\underline{\mathsf{\Delta }}$ was implemented by averaging all of the PS rates within a radius of 250 m from the GNSS station.

A further elucidation needs to be made. The reference system that describes SAR geometry (state vectors) does not account for continental drift making plate movement visible in the interferometric measurements, projected along the LoS. The GNSS data used for the calibration are also not continental drift compensated. Hence, the two data-sets can be considered “compatible”. After merging (also adding $v_{ref}$), the data contains continental drift that could dominant the visualialization. Therefore, plate movement was compensated using the model in [29] to show movement relative to a specific plate (Eurasian). The removed motion is basically an overall offset of several mm/y, plus a light ramp due to the LoS projection of the horizontal motion. This is generally not necessary since it is only a representation issue, but it facilitates the interpretation of the results. For the sake of completeness, the removed model is also shown in Figure 7b, Figure 9b, and Figure 11b.

4.1. Netherlands Datasets

The systematic generation of deformation maps on a national scale using SAR interferometry has become more feasible, as missions like Sentinel-1 provide global coverage on a regular basis. Often, these kinds of services must be combined with existing geodetic measurements in order to make them integrated and comparable [21]. In this work, a typical example is provided where InSAR velocities covering a large area—acquired in descending and ascending passes—are combined with a consistent number of GNSS-derived velocities. The considered datasets are two Sentinel-1A/B stacks covering the northern part of the Netherlands. The area includes many examples of man-induced subsidences, such as in Groeningen [30]. The details of the dataset are described by

Table 1. Netherlands Dataset.

Track	Orbit	Acquisitions	Time [years]	GNSS Stations
088	ASCE	140	3.4	47
037	DESCE	122	3.0	37

. Observing the raw result of the PSI-processing, some spatially correlated deformation signals are visible when observing between $\pm 4$ mm/y. Since no large-scale deformation phenomena are expected in this area, such patterns must therefore be related to the residual atmospheric signal. In the area of interest, a considerably high number of GNSS data are available (see Table 1). The previously described methodology was applied to the data using the available GNSS measurements [27]. The results in Figure 6 and Figure 8 show how the technique permits the further reduction of the residual atmospheric error and the relation of every single PS measurement to the Eurasian plate only. Further information about the measured ${\mathsf{\Delta }}_0$ and the removed continental drift model is also available in Figure 7 and Figure 9. Since the ascending/descending data are now absolute, it would be possible to combine them directly using the error description provided in Equations (9)–(11).

4.2. North Anatolian Fault Dataset

Applications related to the measurement of tectonic movements are the most challenging for InSAR. The requirement of an accuracy of 1 mm/y at more that 100 km (<10 nstrain) pushes the technique to its limits [5]. Notwithstanding, it has been demonstrated that, since actual SAR missions are characterized by very stable oscillators [31] and very good orbit knowledge [32,33], such numbers can be considered achievable. The major limitation is atmospheric effects. If the interferometric time series is not long enough, integration with other data would be necessary in order to fulfill the requirements at very large distances. Moreover, the possibility of referencing the measurements to standard reference systems (e.g., Eurasian Plate, etc.) would extend the usability of the final results.

The technique must then also be demonstrated on an appropriate example for tectonics. The North Anatolian Fault was chosen, as it is a typical case study investigated by many geo-scientists using SAR interferometry [34,35]. The area of interest is covered by a Sentinel-1A/B stripe extending for more than 600 km in the along-track direction, see

Table 2. North Anatolian Fault Dataset.

Track	Orbit	Acquisitions	Time [years]	GNSS Stations
$123/2$	DESCE	134	3.3	5
$123/3$	DESCE	134	3.3	9
$123/4$	DESCE	134	3.3	4

. In order to make the processing feasible, the stripe was divided into three frames. The three frames were processed independently using the PSInSAR technique. The area is very tectonically active and it is very important to preserve all scales of the measured deformation signal. Figure 10 displays the results together with a mapping of the main path of the North and East Anatolian Fault (dashed white line).

5. Conclusions

The results show how the technique is able to both remove dependence on the reference point and mitigate residual atmospheric errors present in the final InSAR results. The described methodology works under some straightforward assumptions on the input data that have to be considered:

• The GNSS/InSAR time series must overlapped in space and time,
• the 3D GNSS velocities and their variances are known, and
• the motion of the GNSS station should be representative of the motion of the area.

The method was studied on both simulated and real data. Simulations show that even with a reduced set of GNSS stations it is possible to retrieve the absolute motion with good accuracy. On the other hand, the capability of the approach to mitigate residual atmospheric errors depends both on the GNSS station density and the spatial correlation of the errors. In order to provide some numbers, the error was modeled as an $AR\left(1\right)$ process with a correlation length of 60 km and variable overall variance. Such numbers were used because they match with the experiments. Two real test cases were then processed to demonstrate the approach by trying to address applications that should benefit from the technique. The results are promising and open the door to an accurate cross-validation between GNSS/InSAR. The extension of the approach to time series would also be interesting, merging InSAR phases with continuous GNSS measurements. The approach could be conceptually equivalent, since the deformation rates are just a scaling of the single-phase measurements, but the spatial covariance of each interferogram would be needed.

For the sake of precision, it is noted that estimation of ${\widehat{v}}_{ref}$ is also generally possible using only InSAR data. The key parameter to be considered here is indeed the length of the time series $T_{obs}$. Theoretically speaking, with a long time series, it is possible to robustly estimate the motion of the reference point by comparing the interferometric deformation rates with the rates estimated from group delays, using co-registration shifts or the PS positions [36,37]. Moreover, if $T_{obs}$ is large enough, the residual atmospheric effects are also sufficiently small to fulfill even the strict requirements [8] set by applications aimed at measuring tectonic movements [5]. Notwithstanding, if $T_{obs}$ is not large enough to fulfill the accuracy requirements at large distances, then merging with GNSS is necessary.

In general, it is possible to conclude that the developed method helps to increase, when necessary, the accuracy of InSAR over large distances [8]. The final product has the spatial coverage of InSAR, and is an optimal combination of InSAR and GNSS information based on the knowledge of the error statistics. This preserves all scales of the deformation signal and also allows characterization of the final product error.