## Abstract

High spatio-temporal variability of atmospheric water vapor affects microwave signals of Global Navigation Satellite Systems (GNSS) and Interferometric Synthetic Aperture Radar (InSAR). A better knowledge of the distribution of water vapor improves both GNSS- and InSAR-derived data products. In this work, we present a collocation framework to combine and retrieve zenith and (relative) slant tropospheric delays. GNSS and InSAR meteorological products are combined aiming at a better retrieval of the atmospheric water vapor. We investigate the combination approach with synthetic and real data acquired in the Alpine region of Switzerland. Based on a closed-loop validation with simulated delays, a few mm accuracy is achieved for the GNSS-InSAR combination in terms of retrieved ZTDs. Furthermore, when real delays are collocated, the combination results are more congruent with InSAR computed products. This research is a contribution to improve the spatio-temporal mapping of tropospheric delays by combining GNSS-derived and InSAR-derived delays.

## 1 INTRODUCTION

The gases in the Earth’s atmosphere induce propagation delays and distortions to microwave satellite signals. Most of these delays (about 90%) are caused by the so-called dry gases and the remaining by water vapor. Although dry gases cause most of the effect, their small spatio-temporal variation allows one to predict the respective delays relatively accurately. On the other side, the water vapor content of the Earth’s atmosphere induces high spatiotemporal fluctuations in microwave signal delays that need to be quantified. These fluctuations in Global Navigation Satellite Systems (GNSS) and Interferometric Synthetic Aperture Radar (InSAR) techniques are quantified during the processing of the measurements. Thus, observations of GNSS and InSAR can provide valuable information about the quantity and the spatio-temporal distribution of the water vapor in the atmosphere.

A better knowledge of the amount of water vapor in the atmosphere is useful in geodetic applications of GNSS and InSAR where it is considered a source of uncertainty. However, this information can potentially enhance meteorological and climatological applications where water vapor is an important signal as the most abundant greenhouse gas and a vital component of the hydrological cycle.

The problem and the effect of “inclined troposphere” (i.e., horizontal gradient in the tropospheric refractivity field) has been discussed quite early (Geiger, 1987), and for several decades now, the use of GNSS for meteorological purposes has been investigated (Davis et al., 1996; Elgered et al., 1997; Kruse, 2004; Yuan et al., 1993) with its major application in meteorology, originally known as GPS meteorology (Bevis et al., 1992). The most basic GNSS meteorological product is the zenith total delay (ZTD), which represents the effect of the integrated amount of atmospheric gases in the zenith direction of the GNSS receiver. Other GNSS meteorological products include zenith delay gradients, measuring the atmospheric tilt in the east-west and north-south directions; slant total delay (STD), which represents the integrated atmospheric quantity in the direction from the satellite to the receiver; as well as 3D atmospheric water vapor content determined using GNSS tomography (Champollion et al., 2005; Heublein, 2019; Kruse, 2004; Möller, 2017; Perler, 2012; Troller, 2004).

One application of spaceborne SAR interferometry is the monitoring of geodynamic and geomorphological processes. Single-polarization SAR imagery contains information on the backscattering intensity and on the propagation time of the signal between the radar antenna and the Earth’s surface. This propagation time is hidden in the phase of the SAR signal, which also contains a phase term that depends on the backscatterer characteristics of the objects that lie in a specific resolution cell of a SAR image. Only by means of SAR interferometry, i.e., by exploiting the phase differences between co-registered SAR images, can the phase due to the propagation delay difference be revealed. The phase is composed of several components, which have to be separated to be able to retrieve the sought information (for instance ground displacement). The high variability of water vapor makes the tropospheric phase one of the most challenging phase delays to be estimated. The most typical approach to separate the tropospheric delay in an InSAR interferometric phase from the other components is Persistent Scatterer Interferometry (PSI) (Ferretti et al., 2001; Kampes, 2005), a technique where different spatio-temporal baselines and statistical assumptions about the spatio-temporal behavior of water vapor are assumed (Cavalié et al., 2007; Hooper et al., 2007). However, InSAR meteorology is a relatively new field. Indeed, (Hanssen, 2001; Hanssen et al., 1999; Zebker et al., 1997) were amongst the first published works that used InSAR for a meteorological purpose. Further contributions (Li et al., 2009; Puysségur et al., 2007) have investigated the calibration of atmospheric effects in InSAR based on external sources such as Medium Resolution Imaging Spectrometer (MERIS) or Moderate resolution Imaging Spectroradiometer (MODIS) or by integrating numerical weather prediction (NWP) models into InSAR applications such as Hobiger et al. (2010), Kinoshita et al. (2013), and Mateus et al. (2013). In more recent research (Bekaert, Walters, et al., 2015; Siddique et al., 2019; Wicks et al., 2002), the InSAR tropospheric delays are interpolated over the entire footprint by modeling their spatial dependence.

The possibility of using GNSS for correcting InSAR delays has been investigated as well (Bekaert, Hooper, et al., 2015; Li et al., 2006; Onn, 2006). This requires a very dense network of GNSS stations, and moreover, interpolating GNSS delays to SAR points is very hard due to the high water vapor variability. Therefore, the combination of GNSS absolute measurements and InSAR relative measurements would be a more reasonable product with the opportunity to exploit the complementary spatio-temporal characteristics. Indeed, whilst a high temporal and sparse spatial resolution characterize GNSS measurements, the InSAR produces difference measurements with a high spatial density and low time resolution.

In this paper, we present a collocation framework that allows one to combine zenith and (relative) slant tropospheric delays, used to retrieve delays in the area of the investigation. The combination of the delays obtained by means of GNSS and InSAR techniques, with the goal of having a better retrieval of the atmospheric water vapor and of the respective delays than each technique individually, is investigated. The increasing interest in InSAR tropospheric delays and its complementary characteristics to GNSS are the main motivation behind this work. We integrate the GNSS and InSAR measurements using the least squares collocation software COMEDIE (Collocation of Meteorological Data for Interpretation and Estimation of Tropospheric Path Delays) (Eckert et al., 1992a, 1992b). COMEDIE has already been successfully deployed for decades to interpolate and extrapolate pointwise or integral atmospheric parameters produced by GNSS, radio sounding, numerical weather prediction models, or other meteorological information sources in four dimensions. For this work, it has been further upgraded to combine GNSS zenith total delays and InSAR difference slant total delays (and absolute slant total delays). Alshawaf (2013) and Heublein (2019) use “absolute” InSAR delays retrieved by combining GNSS and InSAR information, assuming spatio-temporal averaging/filtering by inverting difference maps using least squares. Here, we aim to combine GNSS ZTDs and InSAR double difference slant total delays (ddSTDs) from PSI to exploit their complementary characteristics in terms of resolution; we develop our method to combine GNSS and InSAR data to retrieve the common refractivity field and compute measurement delays experienced by microwave signals in the area of investigation.

In previous works, COMEDIE has been used to cross validate GNSS tropospheric error models from different stations (Hurter, 2014; Wilgan et al., 2017) and to detect stations with wrongly estimated tropospheric delays, i.e., improving the tropospheric models for those stations. Furthermore, as we show in this work, the parameters estimated in the collocation process can be used to produce zenith and/or (relative) slant tropospheric delays at any point of the investigated region. Therefore, they can be practically used for correcting tropospheric delays for any GNSS receiver in the area of investigation. The main limitation, for the time being, is that this correction can be applied in post-processing mode only. Indeed, producing atmospheric products from InSAR takes a long time due to the processing technique, which makes it impossible to produce real-time or near real-time delays.

Following a summary of tropospheric delays in GNSS and InSAR (Section 2), Section 3 describes the study area and the used datasets. Section 4 shows the principles of COMEDIE and the method used for the combination of GNSS and InSAR techniques’ tropospheric products, whilst in Section 5 we investigate important aspects of GNSS-InSAR combination with a set of simulated measurements. Section 6 shows the combination of real data in Valais, Switzerland, and finally, Section 7 sums up this paper.

## 2 TROPOSPHERIC DELAY IN MICROWAVE SIGNALS

### 2.1 Refractivity and path delay

The microwave signals experience a delay when travelling in the troposphere, which is the delay of the signal compared to a signal propagating with the speed of light in a vacuum. This delay is a function of the refractivity, *N*, along the path between the satellite and the station:

1

The refractivity is related to the refractive index, *n*:

2

The delay and the refractivity are caused by the dry gases and the water vapor. Therefore, it can be written as a function of dry and wet refractivity:

3

The refractivity itself is determined by meteorological quantities such as the air pressure, the partial water vapor pressure, and the temperature (Essen & Froome, 1951) as follows:

4

where the first term represents the dry refractivity and the last two terms represent the refractivity value induced by the water vapor. *P*_{tot}, *P*_{wet}, and *T* are the total pressure, the water vapor partial pressure, and the temperature, respectively. The coefficients *k*_{1}, *k*_{2}, and *k*_{3} have been reported in many publications. In this work, we have used the coefficients from Rueger (2002), where *k*_{1}, *k*_{2}, and *k*_{3} are 77.6890 KhPa^{−1}, 71.2952KhPa^{−1}, and 375463 K^{2}hPa^{−1}, respectively.

In case of simulated measurements, to compute the slant delay from meteorological parameters along the path, we have used the discrete formula where the integral is approximated by a sum as follows (Hurter, 2014):

5

with Δ*s* the geometrical distance between two successive points along the slant, *N*_{i} and *N*_{i + 1} the computed refractivities at two successive points, and m the total number of segments between successive points.

### 2.2 Tropospheric delay in GNSS

The path delay experienced by the GNSS receiver is in the slant direction; however, the basic tropospheric parameter used by the GNSS community is the ZTD. The zenith delay is estimated in the GNSS data analysis as an unknown parameter. Therefore, the individual slant delays from each visible satellite are mapped into the zenith direction (Boehm et al., 2006) and separated into a dry and a wet part, which have values of approximately 2.3 m and 0–0.4 m at sea level, respectively (Walpersdorf et al., 2001). In many publications, the zenith delay is separated into a hydrostatic and a non-hydrostatic part where the hydrostatic delay also contains contributions from the water vapor, although many works consider the hydrostatic delay as the dry delay. From the zenith delay and meteorological parameters, the column-integrated water vapor (IWV) in the atmosphere can be derived (Askne & Nordius, 1987).

Several models of the tropospheric delay in GNSS have been developed. Initial contributions from Hopfield (1971) and Saastamoinen (1973) are still the most common models used by the GNSS community. The Saastamoinen model (Saastamoinen, 1973) is based only on surface meteorological parameters and is defined as follows:

6

where z is the off-zenith angle, *ϕ* is the geographical latitude, *h*_{ellip} is the ellipsoidal height, and the other parameters are as defined in the previous section.

More recent works account for improved meteorological parameters and improved mapping functions to model the tropospheric delay such as Collins (1999), Leandro et al. (2009), and Möller et al. (2014). Leandro et al. (2009) is more appropriate for North America, whilst Möller et al. (2014) is based on numerical weather data.

### 2.3 Tropospheric delay in InSAR

Satellite radar signals (e.g., InSAR), typically used to generate maps of Earth’s topography and measure the Earth’s surface deformation, are affected by the tropospheric delay similarly to GNSS. The InSAR interferometric phase component contains several sources that contribute to the propagation delay. Typically, the interferometric phase signal is modeled as 1) a geometric phase component equivalent to propagation in the vacuum space (further split into a reference surface and a topographic component), 2) a line-of-sight displacement component, 3) an atmospheric path delay (further split into ionosphere and troposphere depending on the wavelength of the microwave signal), plus a phase noise component. Therefore, the interferometric phase is

7

with the right side components of Equation (7) being respectively: ∅_{topo} the topographic phase component, ∅_{disp} the phase displacement component of the surface of the Earth between two InSAR acquisitions, ∅_{tro} differential phase delay due to the troposphere, ∅_{iono} differential phase delay due to the ionosphere, ∅_{orbit} phase components due to satellite orbit inaccuracies, ∅_{flat} uncertainties in the Earth’s ellipsoidal reference surface, and ∅_{noise} system thermal noise and loss of coherence between two acquisitions (Ferretti et al., 2001).

The most common technique used to separate InSAR atmospheric phases is PSI. In PSI, persistent scatterers, which are pixels that exhibit long-term temporal coherence, are identified. The PSI-derived APS (Atmospheric Phase Screen) is a model-based retrieval where the LOS (line-of-sight) displacement is typically modeled as a linear function of 1) velocity (e.g., deformation velocity is assumed constant), 2) a height correction, plus 3) a phase constant. The parameters are iteratively estimated by linear regression based on the interferometric phase differences. The residual phase of the regression contains atmospheric phase, non-linear deformation phase, and baseline errors as well as phase noise. The phase standard deviation of the regression residuals is used as a criterion whether to accept or not accept the regression solution (i.e., to keep this PS (persistent scatterer) point or not). Spatially low-frequency filtered and unwrapped residual phases are considered as atmospheric contributions and are taken into account in the phase modeling of the next iteration of the PSI regression analysis. Thus, atmospheric effects can be isolated at PSI locations. Here, the estimated tropospheric delay is not an absolute delay as the interferograms in InSAR processing are built with respect to a single reference image, called the reference image, and moreover, the interferometric phases are calculated relative to a reference point (one of the identified PSI points). Therefore, the InSAR tropospheric delay is a double difference measurement. The differential tropospheric phase can be converted into differential tropospheric delay by using the signal wavelength *λ*_{wavelength}:

8

In other terms, this InSAR delay can be rewritten as (Fornaro et al., 2015)

9

where the dSTD represents the relative slant total delay between any InSAR point and the reference point; *x*, *x*_{ref} are respectively the coordinate vectors of the pixel and reference pixel; and *t*, *t*_{ref} are the time of acquisition of a PS point and the reference image acquisition time, respectively.

## 3 DATASETS

The area of investigation in this work is the Alpine region in Valais, Switzerland, which has an interesting altitude variation of about 3 km. Due to this large variation of altitudes, microwave signals (with relatively close observation sites on Earth) can easily be affected differently by the atmospheric gases, due to the very different atmospheric portions that they travel. For this work, a set of simulated data and a set of real data are utilized. The investigated area is illustrated in Figure 1, where the available GNSS stations are marked in red and purple and the InSAR footprint is represented by the yellow rectangle.

### 3.1 Simulated data

For our tests, the NWP model deployed by MeteoSwiss, COSMO-1 (Consortium for Small-scale Modeling), was used as reference to simulate synthetic measurements for GNSS and InSAR. COSMO-1 is used to produce the weather forecast in the challenging Alpine region (MeteoSwiss 2019). It is based on physical laws (energy conservation, momentum, mass, and simulating phenomena such as phase transitions of water and radiation processes), and it needs correct initial and boundary fields to appropriately predict future events (MeteoSwiss 2019). Atmospheric parameters such as pressure, temperature, and water vapor partial pressure are calculated on a three-dimensional grid (plus time), where vertical distances between the grid points increase at higher altitudes. More specifically, for the simulation of GNSS and InSAR observations, the COSMO-1 model is used with a horizontal grid distance of 1.1 km and vertical altitude of about 11 km with 65 vertical layers. Atmospheric parameters of the COSMO-1 model (pressure, temperature, and water vapor partial pressure) were interpolated with a resolution of 300 m along the path of GNSS and InSAR measurements and were used to compute refractivity as shown in Equation (4). The refractivity values were therefore integrated to produce COSMO-1 based simulated GNSS and InSAR measurements using Equation (5). For the tropospheric medium above the COSMO-1 model, the Saastamoinen model (Saastamoinen, 1973) was used to calculate the remaining path delay as shown in Equation (6).

Therefore, GNSS measurements for the GNSS AGNES/COGEAR networks, as well as InSAR ddSTDs for Valais, were simulated. In this case, we simulated InSAR images with observation angles varying between 24.4° and 28.5°. In total, we produced 60 maps of double difference delays (i.e., one per day) for August and September 2016. In reality, the time resolution of real SAR observations is usually lower than one day; it varies from a few days to weeks (or even months), which is enough for the monitoring of geodynamic processes. However, atmospheric delays are not correlated after one day, and therefore, the simulation of one tropospheric map per day is still realistic for our application.

### 3.2 Real data

For the real data, we have used 29 acquisitions from the COSMO-SkyMed satellite as interferometric measurements for a period of 5 years from 2008 to 2013 (acquired during the June to October time span when it does not snow). The COSMO- SkyMed satellite operates at X-band (*λ* = 3.12 cm), frequencies at which the ionospheric effects are negligible. Meanwhile, during the same time span, GNSS measurements from the Federal Office of Topography (swisstopo) AGNES network (Swisstopo 2019) (Figure 1, red dots) and COGEAR network of ETH Zürich (Figure 1, purple dots) were available hourly.

The InSAR data were processed using the interferometric point target analysis (IPTA) module of the Gamma software (Wegmuller et al., 2004; Werner et al., 2003; Wegmuller et al., 2010), whilst GNSS measurements are ZTDs obtained from processing GNSS observables using the Bernese software (Dach et al., 2015) in a double difference post-processing mode.

## 4 GNSS-InSAR COMBINATION

The main objective of this work is to exploit the complementary spatio-temporal characteristics of GNSS and InSAR (as shown in Figure 2) in order to improve the accuracy of path delays retrieved by these techniques (and estimated water vapor). We propose the combination of GNSS and InSAR measurements in a collocation framework. We use the least squares collocation software COMEDIE (Eckert et al., 1992a, 1992b) to directly combine GNSS ZTDs and InSAR atmospheric ddSTDs produced from PSI.

### 4.1 Least squares collocation software COMEDIE

COMEDIE has been developed and deployed for more than 20 years by the Institute of Geodesy and Photogrammetry at ETH Zurich. It can interpolate and extrapolate pointwise or integral atmospheric parameters produced by GNSS, radio sounding, numerical weather prediction models, or other meteorological information sources in four dimensions (Eckert et al., 1992a, 1992b; Hurter, 2014; Troller, 2004; Wilgan et al., 2017). Equipped with a least squares collocation algorithm, the model is based on a functional and a stochastic component where systematic parts are determined and separated from measurement noise. The main equation of the collocation is (for instance from Troller, 2004):

10

where *l* is the measurement. Its model contains a functional part *f*(*u*, *x*, *t*), representing models of meteorological parameters based on physical realities with *u*, *x*, and *t*, respectively; the state vector to be estimated; the coordinates where the measurements are acquired; and time. The functional part is dependent on height as well as on horizontal and temporal gradients.

The residuals of the measurements are separated into a signal part *s*(*C*_{ss},*x*,*t*) and a noise part *ε*. The first part is the systematic components of the measurements (depending on a modeled covariance *C*_{ss}), whilst the noise is stochastically uncorrelated Gaussian.

#### 4.1.1 Collocation of GNSS ZTD measurements

In case GNSS ZTDs are processed via the collocation model, the deterministic part of the measurements can be described as (from Troller, 2004)

11

The state vector is u = (*ZTD*_{0},*a*_{ZTD}, *b*_{ZTD}, *c*_{ZTD}, *H*_{ZTD}), where *ZTD*_{0} is the ZTD at reference position (*x*_{0}, *y*_{0}, *h*_{0}, *t*_{0}) (an arbitrary chosen point inside the GNSS network), *a*_{ZTD} gradient parameters in the x coordinate, *b*_{ZTD} gradient parameters in the y coordinate, *c*_{ZTD} gradient parameters in time, and *H*_{ZTD} the scale height. *x*, *y*, *h*, *t* are the coordinates and time of a measured point.

The covariance of the signal in this case is (from Troller, 2004)

12

with (*x*_{i}, *y*_{i}, *h*_{i}, *t*_{i}) and (*x*_{j}, *y*_{j}, *h*_{j}, *t*_{j}) space and time coordinates of two measured points; Δ*x*_{0}, Δ*y*_{0}, Δ*h*_{0}, Δ*t*_{0} the correlation lengths in space and time; σ_{0} *a priori* covariance of signal; and *h*_{0} the scale height modifying the correlation lengths as a function of heights. In Hurter (2014), it is explained how the parameters of the covariance matrix can be determined from the data.

#### 4.1.2 Collocation of slant measurements

For this study, COMEDIE has been further developed with the goal to include in the collocation software slant delays. This would not only be beneficial for combining GNSS ZTDs with InSAR slant measurements, but in the future work, we aim to use other slant measurements provided from GNSS, spectrometers, or any other technique.

A slant delay is expressed as the zenith delay mapped into the direction of the satellite with an appropriate mapping function, which is a geometric factor describing the dependence of the delay from the elevation angle. Therefore, the functional and stochastic parts of slant measurements at a 4-D location (*x*, *y*, *h*, *t*) are modeled as follows:

13

14

with the state vector equal to u = (*ZTD*_{0}, *a*_{ZTD}, *b*_{ZTD}, *c*_{ZTD}, *H*_{ZTD}), same as for zenith measurements. The covariance between two slant delays is modeled by the covariance of the zenith delays at the slant positions *C*_{ss}(*ZTD*_{1}, *ZTD*_{2}) and the respective mapping functions *MF*_{1}, *MF*_{2}.

#### 4.1.3 Collocation of InSAR measurements

For InSAR, the inputs to the collocation process are not absolute slant measurements, but are difference delays as shown in Equation (9). The PSI estimated atmospheric phase describes a double difference slant delay, with respect to a reference PSI point at location *p*_{ref} = (*x*_{ref}, *y*_{ref}, *h*_{ref}) and a reference acquisition image at time *t*_{ref}. Therefore, the functional part can be rewritten as

15

In this case, two sets of unknown parameters have to be estimated in the state vector:

16

where the subscripts (ref and 1, respectively) denote the batches of the reference acquisition image at time *t*_{ref} and another image acquired by the InSAR technique at time *t*1.

The covariance function between two ddSTDs is computed by performing variance (covariance) propagation of the variance (covariance) of slant delays forming the double difference slants. It results in the sum of the covariances of absolute STDs used to describe the double differences for the first and second measurements with its final form as follows:

17

#### 4.1.4 Collocation of GNSS and InSAR measurements

The correlations between the GNSS ZTDs and PSI ddSTDs, are described by the covariance function as follows:

18

or for a GNSS ZTD in the second timespan of the InSAR ddSTD (same time batch as the reference acquisition, i.e., *t*_{ref} in Equation (9)):

19

with the state vector the same as the one in Equation (16) for InSAR observations.

## 5 SIMULATED MEASUREMENTS: EVALUATIONS AND RESULTS

For this research, in order to validate our method to combine zenith delays and/or slant delays, we performed some extensive tests based on simulated measurements. In this context, the scope of this section is to investigate mainly the algorithms/models used in COMEDIE. Real GNSS and InSAR measurements contain errors related to algorithms, specificity of the scenario (such as Alpine region), and equipment/technology issues (for example, antenna phase center variations or orbital/clock errors not perfectly modeled). Therefore, considering especially the new field of InSAR meteorology (path-retrieval), we use simulated delays for which we have more control. It must be pointed out that the simulated delays also contain errors due to the simulation chain, such as quantization errors, coordinate approximation errors, and errors in the non-perfect numerical weather prediction models. However, these errors are more obvious (and smaller) than the ones in real data.

### 5.1 Comparison of different mapping functions

Here, we propose the consideration of slant delays into the collocation process. Our model considers an important element, which is a geometric factor mapping zenith delays into slant delays. This mapping function has to model accurately the electrical path length at geometric elevation angle *ε* to the electrical path length in the zenith direction.

The most basic mapping function is the sine mapping function. It assumes a planar atmosphere, which has been proven to be an accurate approximation for high elevation angles. More accurate mapping functions have been developed, which can be classified into navigation mapping functions and geodetic ones. Mapping functions used for geodetic purposes require knowledge of meteorological parameters, while the ones used for real-time navigation are simpler and less accurate (Guo & Langley, 2003). In this work, we have investigated the impact of several mapping functions. We have considered the sine function, the geometrical mapping function (Hirter, 1998), and the B&E (Black and Eisner) mapping function (Black & Eisner, 1984), which are very practical for navigation purposes. Moreover, the Niell mapping function (Niell, 1996), UNB3m mapping function (Leandro et al., 2006), and the global mapping function (GMF) presented in Boehm et al. (2006) are as well investigated. They are given as a three-term continued fraction and require additional parameters such as station location, day of the year, or the hemisphere in which the station is located. All the investigated mapping functions are summed up in Table 1. More details about the coefficients of the Niell, UNB3m, and GMF mapping functions can be found in the respective papers, Leandro et al. (2006), Niell (1996), and Boehm et al. (2006).

In this work, we compare the mapping functions to understand their impact on the accuracy of our slant delay models (Equations (15), (17)–(19)) at different elevation angles in our area of investigation. To carry out this evaluation, we have used the numerical weather prediction model COSMO-1. From COSMO-1 data, we have simulated slant delays (at different elevation and azimuth angles for 13 days during August and September 2016) and mapped them to the zenith direction. Thus, we compare the mapped ZTDs with the reference one calculated from COSMO-1 meteorological parameters.

Figure 3 shows the impact of the sine mapping function on the total, dry, and wet delays, respectively. From this figure, we can notice the uniformity of the dry delay for different azimuth angles and the high spatial variability of the wet delay. For the total delay, its similarity to the dry delay (it is about 90% of the total delay) is visible at a first glance, but high wet delay variations are reflected if we look more carefully (wet delay is about 10% of the total delay).

Figure 4 shows the impact of the six investigated mapping functions on the wet delay. The high variability of the error (difference between mapped ZTDs and COSMO calculated ones) for wet delays at different azimuth angles is confirmed for all functions, caused by the high spatial variability (and therefore unpredictability) of water vapor. We can notice a similar performance of the three geodetic mapping functions, whilst the simpler functions show less agreement among each other or with the geodetic ones. The differences are more visible at low elevation angles, whilst at high elevation angles (**≥**45°) all mapping functions have an absolute error below 20 mm.

The errors of the ZTDs for all mapping functions are plotted in Figure 5. We can clearly notice an improvement with more complex mapping functions at low elevation angles for the total and dry delays. However, at high elevation angles (i.e., **≥**45°), all six mapping functions perform similarly with error statistics of about 1 mm. Therefore, at InSAR elevation angles (which are larger than 60° in our study), it is not necessary to increase the complexity of the mapping factor more than the sine function. It must be noted that for possible slant measurements of lower elevation angles (such as GNSS STDs), more sophisticated models should be used.

### 5.2 Results for simulated measurements

The results presented here consider a dataset of synthetic GNSS and InSAR measurements as explained in Section 3. Simulated GNSS ZTDs and InSAR ddSTDs were combined in COMEDIE, outputting collocation parameters that can be used to compute tropospheric (slant or zenith) delays at any point in the area of investigation. In this work, ZTDs were interpolated at COSMO-1 grid coordinates. We have also calculated delays in the COSMO-1 grid, using the NWP model meteorological parameters. Therefore, the ZTDs produced by GNSS-InSAR combination are directly compared with these delays, considered reference (or “true”) delays. The flowchart of the work for simulated measurements is provided in Figure 6, where we are proposing a closed-loop validation of our results.

We estimated (interpolated) ZTDs with COMEDIE at COSMO-1 grid points from the combined InSAR and GNSS measurements. Then, we computed the differences of these estimated ZTDs with COSMO-1 calculated ZTDs. Figure 7 presents the standard deviation over time of these differences, which is plotted following the topography of the terrain. The figure shows that the differences are generally smaller at higher altitudes and larger at lower ones. At higher altitudes, the total delay is generally smaller (smaller portion of atmosphere traveled), and therefore, the error of estimated delays is smaller; however, the error to total delay ratio is similar at all altitudes. Figure 8 shows the time series of these differences.

As shown in Figure 7, using InSAR-GNSS combination, the standard deviations of the estimated ZTD errors (compared to COSMO-1-calculated ZTDs) have values between 1.3 mm and 3.8 mm. The differences of estimated ZTDs (over all grid points and epochs) vary in the interval [−19 mm; 12 mm], as illustrated in Figure 8. The mean standard deviation is 2.1 mm, whilst the mean bias over all points is insignificant, with a spatial standard deviation of 1.8 mm. These statistics have been computed over 531 COSMO-1 grid points, which are inside the InSAR image area.

Additional to the comparison of the COMEDIE produced delays (after the combination) with COSMO-1 calculated delays, we have also investigated the residuals of our used ZTD measurements. They report how well the functional model fits with the measurements. In Figure 9, the respective residuals (top of figure) as well as the corresponding standard deviation (bottom of figure) are displayed. The residuals vary in the interval [−13.1 mm; 20.7 mm] and have standard deviations between 1.8 mm and 7.7 mm for the investigated timespan.

Figure 10 displays the experienced ZTDs and ZWDs at the AGNES/COGEAR stations (shown in Figure 1). The corresponding variability of these delays (in terms of standard deviation) is depicted in Figure 11. The zenith wet delays vary between 10 mm and 220 mm with corresponding standard deviations between 20 mm and 50 mm. As for the total delays, they have a high variance range of about 30 mm similar to the one of the wet component (with values in the interval [235 mm; 262 mm]). The main part of the standard deviation for absolute ZTDs comes from the dry gases in the atmosphere, which cause different dry delays mainly due to the different topographies of the GNSS stations. The residuals of the ZTDs have a mean standard deviation of about 4 mm leading to a 3*σ* of 12 mm while the ZWDs vary in a 200-mm interval. This shows a good fitting of the model used in Equation (11), which is able to capture a good part of the variability of ZWDs.

### 5.3 InSAR parameters evaluation

Although the number of identified PSI points may be very large, the PSI points are not distributed over all the footprints of the image with the same sample spacing, and furthermore, some InSAR image footprints are larger than others depending on the mission. Tables 2 and 3 summarize the results of the InSAR ddSTD simulations with different PSI points sample spacing and InSAR images’ sizes, respectively. We have simulated InSAR ddSTD maps in the SAR foot print of Figure 1, where we assume a PSI point is available every 400 m, 800 m, 1200 m, 1600 m, 2400 m, and 3600 m. Moreover, we decreased the footprint size by 1.5, 2.4, 4.4, 11, and 31 times.

Tables 2 and 3 show that considering a smaller PSI points sample spacing or ddSTD map size does not considerably affect the results. Considering that PSI points are not distributed everywhere with a high sample spacing and that in many cases the PSI points are located in one part of the footprint this is an interesting result, which needs further confirmation in future simulations or tests.

## 6 COMBINATION OF REAL DATA

The results presented in this section are derived from the processing of InSAR observables from the COSMO-SkyMed satellite and GNSS measurements from AGNES+COGEAR stations (see Figure 1) for the dates of InSAR acquisitions between 2008 and 2013, as explained in Section 3. In this case, we do not have reference measurements (as in the case of simulated data), but we compare qualitatively the InSAR-GNSS combination outputs with ddSTDs estimated from GNSS and InSAR separately. The ddSTDs estimated by GNSS are obtained by processing AGNES GNSS ZTDs in COMEDIE and interpolating them to the coordinates of PS locations. In addition, they are mapped into the slant direction. It must be pointed out that for InSAR we have more than 300,000 PSI points, but only few points (about 1,000) are used for the combination. During the iterative PSI processing, the APS is obtained by repeated spatial filtering and unwrapping of the residual phase obtained from linear regression at each PS point. At each iteration, the updated atmospheric screen is then taken into account in the phase modeling of the regression. Therefore, the atmospheric phase signal is spatially smooth, and it is safe to select only a subset of PS points. Moreover, this is also supported by the results in the previous section where we showed that we do not benefit a lot from using all the InSAR data. Using only a subset of points helps to accelerate the measurement processing in COMEDIE. Therefore, the interpolation and comparison of ddSTDs is done at InSAR PSI points not taken into account in the combination. The flowchart for this set of data is provided in Figure 12.

From GNSS and InSAR ddSTD comparisons (Wilgan et al., 2019), we noticed that there are days when the agreement between InSAR estimated ddSTDs and GNSS estimated ddSTDs is relatively good and other days when this agreement is poor. This is displayed in Figure 13 where we show 2D maps of ddSTDs computed by GNSS, InSAR, and their combination for two selected days. In September 2011, we can notice a good agreement between GNSS and InSAR separately computed ddSTDs. However, in September 2008, the variability of estimated 2D ddSTD maps is very low and much higher for the GNSS and InSAR techniques, respectively. It must be pointed out that the high ddSTD variations shown in Figure 13 are not only due to spatial water vapor variability. Indeed, for this mountainous area, large variations of ddSTDs are expected as the terrain varies quite notably (Figure 1). Thus, the height component directly constitutes an important part in the ddSTD variability. This is important as in our models in COME-DIE we model the measurements (and their correlations) dependent on a height component. It is important to point out that after removing the height component from the InSAR and GNSS-derived ddSTDs shown in Figure 13 the agreement (depending on the dates) between the residual delays does not change compared to the ones shown here. This means that the ddSTDs’ height component (important in Alpine regions) is not the main reason of disagreement.

Figure 14 shows the histograms of estimated ddSTDs for the two cases. In the case where the ddSTDs, separately estimated from GNSS and InSAR, agree relatively well, the histogram of their combination lies in the middle of their individual histograms. However, it is interesting to notice that when their agreement is poor, the histogram shape of the combination follows mostly the histogram of the PSI estimated double delays. Indeed, this is the case due to a much higher density of InSAR points (several hundred up to few thousands of PSI points, compared to few tens of GNSS observations), which directly leads to a higher weight in the estimation process. As expected, high spatial variations that cannot be captured and interpolated accurately by the sparse GNSS network can be modeled more accurately by high density InSAR measurements.

Figure 15 shows the time series of all estimated ddSTDs from GNSS, InSAR, and their combination. We can confirm that for all the images, the estimated delays from the combination follow more closely the InSAR estimated delays, but there is clearly a smoothing by the contribution of GNSS data in the combination.

## 7 DISCUSSION AND CONCLUSIONS

In this paper, we presented a collocation framework to combine and produce (absolute and relative) zenith and slant tropospheric delays. More specifically, we investigated the case of combining GNSS-derived and InSAR-derived estimates of tropospheric path delays. The increasing interest in InSAR tropospheric products and the complementary spatio-temporal characteristics of GNSS and InSAR are the main motivations behind this case study. We integrated the (relative) path delays individually retrieved from GNSS and InSAR/PSI in the least squares software COMEDIE, further upgraded to process GNSS zenith total delays and InSAR difference slant total delays simultaneously. The combination of GNSS and InSAR delays in a collocation approach is a new aspect of this research.

The area of investigation is located in the Swiss Alps, over the Matter Valley, which has a high variability of altitude; therefore, also a high variability of atmospheric delays. This is quite a challenging scenario. Moreover, the time of investigation corresponds to a period between June and October. Therefore, most of the delays are estimated for hot and humid summer days, where the situation is turbulent and complicated, leading to high variability of the derived tropospheric delays.

In this work, we used a set of simulated measurements in order to validate our method to access only the algorithms/models and assumptions used in COME-DIE. Reprocessed data from the NWP model COSMO-1 were used for the simulation of GNSS and InSAR delays as well as for comparison purposes of retrieved delays. Based on simulated data, a few mm accuracy (in terms of mean and standard deviation) is achieved in the case of GNSS-InSAR combined techniques. This accuracy is obtained in terms of retrieved ZTDs, compared to calculated ZTDs from the NWP model COSMO-1. Furthermore, from tests performed in simulation mode, we investigated important aspects of GNSS/InSAR delay integration in COMEDIE. Initially, we report that lower SAR footprint size or PSI points’ spatial sampling does not significantly change the results. Moreover, the investigated mapping functions (sine, geometrical, B&E, Niell, UNB3m, and GMF) have similar error accuracy at InSAR elevation angles (above 45°); however, for slant measurements with lower elevation angles, using a more sophisticated model is necessary.

In addition, we displayed the combination of real delays from GNSS and InSAR. When the tropospheric estimations of the two techniques are integrated and the individual outputs of GNSS- and InSAR-estimated ddSTDs are different, the combination products are more congruent with the InSAR output. Indeed, due to its very high pixel resolution, InSAR measurements can capture more accurately the spatial variations of the troposphere than sparse GNSS networks. Due to the high spatial density, they are automatically weighted more in the estimation process. It must be pointed out that in the work presented here, we assume that InSAR delays have a relatively good quality. Behind this assumption lies the fact that there is no ground truth to better access the quality of these products. However, in case there is more information about the quality of the InSAR delays, we can discard and/or directly change the measurement weight in COMEDIE.

### 7.1 Further considerations

InSAR meteorology is a relatively new field of investigation. Therefore, InSAR-derived tropospheric products are not fully standardized. The advancing of InSAR processing algorithms to derive tropospheric products is a good motivation for researchers, who aim to exploit this technique and its capabilities to sense the atmosphere.

However, the framework presented here can be adopted for different available measurements; the GNSS and InSAR combination is one case. For instance, using only GNSS ZTDs and/or STDs, (relative) zenith and slant delays can be retrieved everywhere in the area of investigation and can be used to correct/check tropospheric delays for any GNSS receiver in the area of investigation.

Moreover, the products of this framework can also be used for meteorological purposes, i.e., assimilated in numerical weather models, to ameliorate forecast prediction shortcomings. Slant delays are very interesting as they result from travel through a bigger portion of the troposphere and experience horizontal variations of atmospheric water vapor (this is not the case for zenith delays).

## HOW TO CITE THIS ARTICLE

Shehaj E, Wilgan K, Frey O, Geiger A. A collocation framework to retrieve tropospheric delays from a combination of GNSS and InSAR. *NAVIGATION*. 2020;67:823–842. https://doi.org/10.1002/navi.398

## ACKNOWLEDGMENTS

The authors would like to thank the Swiss National Science Foundation (SNSF) for financing this work (project number 200021E-168952), swisstopo for providing the GNSS AGNES/COGEAR measurements, and MeteoSwiss for the COSMO-1 data. The data stack used in this work is a COSMO-SkyMed Product ©ASI – Agenzia Spaziale Italiana - (2013). The PSI processing of the interferometric data was performed by Tazio Strozzi, Gamma Remote Sensing. The work presented in this paper is also part of the AtmoWater project, a collaboration between the Institute of Geodesy and Photogrammetry at ETH Zürich (IGP-ETHZ) and several institutes (GIK, IPF and IKM-IFU) at Karlsruhe Institute of Technology (KIT).

## Footnotes

↵∗ Karina Wilgan is now at Technische Universitat Berlin (TUB), Berlin, Germany and GFZ German Research Centre for Geosciences, Potsdam, Germany.

**Funding information**Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung, Grant/Award Number: 200021E-168952

- Received December 20, 2019.
- Revision received August 9, 2020.
- Accepted August 15, 2020.

- Copyright © 2020 The Authors. NAVIGATION published by Wiley Periodicals LLC on behalf of Institute of Navigation.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.