Joint Exploitation of SAR and GNSS for Atmospheric Phase Screens Retrieval Aimed at Numerical Weather Prediction Model Ingestion

This paper proposes a simple and fast method to estimate Atmospheric Phase Screens (APSs) by jointly exploit a stack of Synthetic Aperture Radar (SAR) images and a dataset of GNSS-derived atmospheric product. The output of this processing is conceived to be ingested by Numerical Weather Prediction Models (NWPMs) to improve weather forecasts. In order to provide wide and dense area coverage and to respect requirements in terms of spatial resolution of ingestion products in NWPMs, both Permanent Scatterers (PSs) and Distributed Scatterers (DSs) are jointly exploited. While the formers are by definition stable targets, but unevenly distributed, the latter are ubiquitous but stable only within a certain temporal baseline that can vary depending on the operational frequency of the radar. The proposed method is thus particularly suited for C, L, and P band missions with low temporal baseline between two consecutive acquisitions of the same scene: these conditions, that are both necessary to provide the dense space-time coverage required by meteorologists, allow for a reliable and robust estimation of APSs thanks to the intrinsic limitation of temporal decorrelation. The proposed technique integrates Zenith Total Delay (ZTD) products computed on a very sparse grid from a network of GNSS stations to correct for SAR orbital errors and to provide the missing phase constant from the derived APS map. In this paper, the complete workflow is explained, and a comparison of the derived APSs is performed with phase screens derived from state-of-the-art SAR processing workflow (SqueeSAR®).


Introduction
A radar signal is affected by propagation delay when it passes through the atmosphere [1][2][3]. For displacement and elevation measurements (the most common products of InSAR processing), this delay is superimposed to the signal of interest and it needs to be removed or mitigated [4][5][6][7][8].
In recent years, several works have been done to integrate SAR measurements and Numerical Weather Prediction Models (NWPMs): the latter can provide an estimate of the atmospheric delay that can be used to correct SAR-derived products [9][10][11], while the first can deliver a measure useful as ingestion product into NWPMs where no other sources like ground-based radars or radio probes are available [12][13][14]. The measure that a SAR can provide is the delay induced by a variation of refractive index in the geometric path traveled by the electromagnetic wave, which, in turn, is related to the water vapor content in the atmosphere.
In this contribution, we will focus on the second topic, and in particular on the algorithms and processing needed for the extraction of the so-called Atmospheric Phase Screens (APSs).
As will be exhaustively discussed in Section 2, the quality of the estimated atmospheric maps is tightly conditioned to the target phase stability between acquisitions. The extracted APSs are, in fact, the difference between atmospheric conditions in two different time instants: this is an intrinsic condition in interferometric processing (InSAR).
The so-called Permanent Scatterers (PSs) are targets showing high stability over long periods and, thus, high accuracy in the estimate of the APS can be reached, but PSs are very often sparse or even not present in the scene (for example in forested or agricultural areas). This sparseness prevents the estimation of dense atmospheric maps and thus the fulfillment of spatial requirements for NWPM ingestion. This makes the use of DSs a mandatory choice. However, we then have to cope with the temporal stability of a DS scene, that may last from hours to months, depending on the scene vegetation and the wavelength [15][16][17].
The proposed method not only uses a combination of PSs and DSs, respecting in this way the constraints on spatial resolution, but it uses also GNSS-derived Zenith Total Delay (ZTD) measurements to properly correct for orbit errors: this further step is a necessary condition when working on large scales since an error in the determination of the orbit can generate low spatial frequency errors that are able to heavily corrupt the estimated APSs.
A boost in the so-called InSAR meteorology has been provided by the launch of the constellation Sentinel-1 A/B: large area coverage (in the order of 170 × 250 km in IW mode) and systematic acquisitions over land with an interferometric revisit of 6/12 days, that is well below the temporal decorrelation of a DS in C-Band [18].
Several works have been done in literature in the field of mapping Precipitable Water Vapor (PWV) using SAR acquisitions aimed at the ingestion into NWMPs [14,19,20].
These works generally exploit a set of interferograms formed between interferometric couples with short temporal baselines in order to reduce temporal decorrelation and thus noise in the estimate. While this technique is effective, it involves the integration over time of APS maps that in turn increases noise in the estimate. Moreover, orbital errors are usually not considered since they can be compensated locally by a stable ground control point. Such approach is not viable for us since it would remove APSs as well, which is the signal of interest.
The current state-of-the-art for the joint PS/DS processing is represented by advanced DInSAR (Differential InSAR) methods like SqueeSAR® [21]: while the processing is very similar to the proposed method, it is quite demanding on a computational level since it requires the exploitation of a lot of images in order to reduce effects that are not due to atmosphere (e.g., subsidence).
The proposed method relies on the spatio-temporal characteristics of the atmospheric phase and the processes that regulate the observable (the SAR phase) leading in this way to a robust and simple method for APS extraction.
This paper is organized into six sections: In Section 2, the requirements in term of spatial and temporal resolution for the ingestion into NWPMs are described and a review of the main properties of the atmospheric signal as seen by a SAR are reviewed, in Section 3, the characterization of the target is carried out by differentiating PS and DS and their decorrelation properties. In Section 4, the whole processing chain is explained accurately, assumptions are clarified and justified. In Section 5, the results derived from the processing of a real dataset are reported. The last section is about overall discussions and conclusions.

Requirements for NWPM Ingestion
A meteo-hydrological product must respect some requirements in terms of horizontal and temporal resolution in order to be useful in the process of ingestion into an NWPM.
The Observing Systems Capability Analysis and Review (OSCAR) tool, developed by the World Meteorological Organization (WMO) set the minimum horizontal resolution for Integrated Water Vapor (IWV) maps in high-resolution NWPMs at 20 km with a goal, over which further improvements are not necessary, of 0.5 km [22].
Concerning time resolution, the minimum requirement is 6 hours with a goal at every 15 minutes. The requirements are summarized in the Table 1. Table 1. Requirements for ingestion products into Numerical Weather Prediction Models (NWPMs).

Requirement for High Res NWPMs Threshold Breakthrough Goal
Temporal Resolution 6 h 60 min 15 min Spatial Resolution 20 km 5 km 0.5 km While GPS-derived Zenith Total Delay (ZTD), radio probes measurements and ground-based radar measurements are almost pointwise but continuous in time, a space-borne SAR system is only able to fulfill the horizontal resolution requirements at the goal level, while the temporal resolution is measured in days since it coincides with the acquisition frequency of the satellite over the same scene.
A possible innovation in this field could be the Geosynchronous SAR that provides revisit in the order of hours and not days allowing the fulfillment of both temporal and spatial resolution requirements [7,8].

Atmospheric Contribution in SAR Images
A radar image is acquired by an active radar placed on-board airborne or space-borne platforms. The raw data acquired at time t is processed to form the final focused Single Look Complex (SLC) image where a pixel at position P is represented in polar form as: x(P, t) = a(P, t)e jφ(P,t) This image contains amplitude (a(P, t)) and phase (φ(P, t)) information: while the former is proportional to the intensity of the reflection, the latter instead contains the information about the delay time of the electromagnetic wave from radar to ground and back [23,24].
This delay is the equivalent one that considers both the geometrical distance between the satellite and ground and the excess distance induced by the earth's troposphere and ionosphere refractive index. While the former shows spatial and temporal variability due to turbulent mixing, the latter is smoother and less important for NWPM ingestion since it does not bring information about the water vapor content of the atmosphere [25].
The phase of a single pixel is not only containing the information about the distance between the sensing platform to the ground but contains also a term due to the target itself: since SAR is a coherent system, this phase term will be formed by the coherent superposition of all the elementary scattered signals inside the resolution cell.
The acquisition geometry is depicted in Figure 1. The phase of a point at time is then defined as: Where ( ) is the portion proportional to the optical path length between the satellite and the target, while ( ) is the target own phase. The first term can be written as: where is the speed of light in vacuum, ( ) is the speed of light that depends of the refractive index along the satellite to target path, τ is the one way satellite-target time of delay, 0 is the radar central frequency, ( , ) is the speed of light through a medium, ( , ) is the refraction index along the satellite-to-target path, 0 ( ) is the difference between the geometrical distances from the target of master and slave (it can be very well approximated by ∥ in Figure 1), and : is the equivalent extra path due to the presence of the atmosphere. Different conditions in pressure, temperature, and humidity in the path traveled by the electromagnetic wave change the refractive index of the medium and thus the phase of the signal. Since the phase of the target ( , ) is unknown, the phase of the SLC image cannot be used to estimate directly the component due to the atmosphere. An interferometric processing is therefore necessary, where differentiation is done between two images acquired in two separate time instants. The interferometric phase is now proportional to the difference in optical path length in the two measures: The phase of a point P at time t is then defined as: where φ R (t) is the portion proportional to the optical path length between the satellite and the target, while φ T (t) is the target own phase. The first term can be written as: where c is the speed of light in vacuum, c(n) is the speed of light that depends of the refractive index along the satellite to target path, τ is the one way satellite-target time of delay, f 0 is the radar central frequency, v(l, t) is the speed of light through a medium, n(l, t) is the refraction index along the satellite-to-target path, R 0 (P) is the difference between the geometrical distances from the target of master and slave (it can be very well approximated by B in Figure 1), and: is the equivalent extra path due to the presence of the atmosphere. Different conditions in pressure, temperature, and humidity in the path traveled by the electromagnetic wave change the refractive index of the medium and thus the phase of the signal.
Since the phase of the target φ T (P, t) is unknown, the phase of the SLC image cannot be used to estimate directly the component due to the atmosphere.
An interferometric processing is therefore necessary, where differentiation is done between two images acquired in two separate time instants. The interferometric phase is now proportional to the difference in optical path length in the two measures: The estimation of the APS is conditional on the following hypothesis: (1). ∆R 0 (P) can be removed easily from the interferometric phase by knowing the acquisition geometry and a Digital Elevation Model (DEM) of the scene. (2). The ionospheric contribution to the interferometric phase is dispersive, and then can be estimated by multi-frequency approaches [26,27]. (3). The phase of the target remains stable between the two acquisitions (∆φ T (P) ≈ 0).
The last hypothesis is the most stringent since the temporal behavior of a target is unknown. An indicator of the stability of the target is the so-called coherence. A high coherence (near to 1) means that the target remains stable over time, while low coherence (near to 0) means that the target has changed between two acquisitions and the interferometric phase is noisy [23,24,28]. An accurate target characterization is carried out in Section 3.

Orbit Accuracy Requirements
In order to properly remove ∆R 0 (P) in Equation (5) the baseline needs to be derived from the satellite's orbit state vectors (OSVs). An error in the determination of the baseline caused by and error in the determination of the position of the satellite generates a trend in the interferometric phase.
Bähr and Hanssen [29] demonstrated that, in the case of Sentinel-1 IW mode, an error in the determination of the normal baseline of just 11 cm generates a fringe in the range direction while an error in the derivative of the parallel baseline of just 1.1 mm/s generates a fringe in azimuth (see Table 2). A fringe is an error equivalent to half a wavelength that, at the nominal frequency of Sentinel-1, is about 2.8 cm. The Precise Orbital Products (POD) available after 20 days from the acquisition allows for accuracy in the 3D orbit determination of 5 cm RMS (Root Mean Square). This figure is referred to the single acquisition, so the accuracy in an interferometric framework is even worse. In order to work with large scale interferograms, it is important to estimate and remove the error induced by orbit indetermination and the procedure will be explained in the methodology section of this contribution.

Target Characterization
As stated in the previous section, a radar is an active system that records the echo coming from a scattering process happening on the ground. Two scattering models can be distinguished: point scattering and distributed scattering. In the point scattering case, a single object in the resolution cell is dominating the radar measurement, while in the distributed case a number of small reflections are summed up coherently to form the final measure.
Both scattering mechanisms can be, over time, coherent or incoherent. A measure of the correspondence between two complex measurements is the so-called coherence. When the scattering object on the surface does not change significantly between one acquisition and the other, the signal is said to be coherent in time. The opposite is total incoherence.
In past years, several techniques have been developed [4] to exploit a set of natural or artificial targets called Permanent Scatterers (PSs) that can remain stable for long periods. On these targets the term ∆φ T (P) in Equation (5) is negligible, driving in this way an accurate estimate of the tropospheric delay.
Even if the spatial density of PSs is very high in urban areas, this is not always enough in rural, vegetated or mountainous areas where instead the vast majority of scatterers are decorrelating and distributed (DS). To provide a wide and dense area coverage the processing of DSs is thus mandatory.
A model describing the decorrelation of a general distributed scatterer has been described in [18,30].
In the first place, they suppose that the decorrelation mechanism is primarily due to the motion of the scatterers as in the Brownian motion: for this reason, the autocorrelation (the coherence) is modeled as an exponential function in time in the form of: where ξ(t) is the complex reflectivity function and T the temporal sampling interval. For a Brownian motion in the look direction with a standard deviation of 1 mm 2 /day the decay constant τ turns out to be of about 40 days in C-Band. While this assumption is true for a general distributed target like short vegetation, it is no more valid for other kinds of DS such as snowbank or water bodies that show decorrelation time much lower than 40 days (hours for snow, milliseconds for water) [16]. This means that, with just temporal decorrelation due to the movement of the scatterers inside the resolution cell, the target loses most of its correlation after just 40 days. With a sampling interval T = 6 days the number of exploitable images is between 7 and 10.
The model can then be extended to consider also thermal noise: In the following section, we will exploit these considerations to provide a minimum number of looks to be used in order to retrieve the atmospheric phases with the proper accuracy.

Processing Scheme
In Figure 2, the whole processing scheme for the extraction of differential APS maps is depicted. Each step will be explained in this section of the paper. A model describing the decorrelation of a general distributed scatterer has been described in [18,30].
In the first place, they suppose that the decorrelation mechanism is primarily due to the motion of the scatterers as in the Brownian motion: for this reason, the autocorrelation (the coherence) is modeled as an exponential function in time in the form of: where ( ) is the complex reflectivity function and the temporal sampling interval. For a Brownian motion in the look direction with a standard deviation of / the decay constant turns out to be of about 40 days in C-Band. While this assumption is true for a general distributed target like short vegetation, it is no more valid for other kinds of DS such as snowbank or water bodies that show decorrelation time much lower than 40 days (hours for snow, milliseconds for water) [16].
This means that, with just temporal decorrelation due to the movement of the scatterers inside the resolution cell, the target loses most of its correlation after just 40 days. With a sampling interval = 6 days the number of exploitable images is between 7 and 10.
The model can then be extended to consider also thermal noise: In the following section, we will exploit these considerations to provide a minimum number of looks to be used in order to retrieve the atmospheric phases with the proper accuracy.

Processing Scheme
In Figure 2, the whole processing scheme for the extraction of differential APS maps is depicted. Each step will be explained in this section of the paper.

Coregistration and Ionospheric Phase Compensation
As a first step, a preprocessing procedure is applied to the data: this workflow consists of the simple coregistration of all the images with respect to a single master image and subsequent debursting of the images.
Since the only interesting contribution to the atmospheric phase is the tropospheric one, ionospheric phase compensation needs to be carried out. The well-known Split Spectrum Method [27] is employed in this scenario to properly estimate and remove the ionospheric contribution to the interferometric phase.

Coregistration and Ionospheric Phase Compensation
As a first step, a preprocessing procedure is applied to the data: this workflow consists of the simple coregistration of all the images with respect to a single master image and subsequent debursting of the images.
Since the only interesting contribution to the atmospheric phase is the tropospheric one, ionospheric phase compensation needs to be carried out. The well-known Split Spectrum Method [27] is employed in this scenario to properly estimate and remove the ionospheric contribution to the interferometric phase.

Topographic Phase Compensation
It is possible to demonstrate that, if the satellite looks at the same spot on ground with slightly different angles as in Figure 1, a phase term related to topography arises in the interferometric phase: where B ⊥ is the normal baseline as in Figure 1, R is the slant range distance, θ is the local incidence angle, and q is the height of the target over the reference plane. This phase term is included in the phase term related to ∆R 0 (P) in Equation (5). This term can be used to generate Digital Elevation Models (DEMs) or since the final output product of the processing is not a DEM, it must be compensated using a known DEM for the area of interest. In our case, a high-resolution SRTM DEM (Shuttle Radar Topography Mission Digital Elevation Model) is used.
It is useful to recognize that, even if the DEM used is not very accurate both in horizontal and vertical resolution (in our case it is below 10m RMS), the error induced on the phase is negligible since the baseline are usually in the order of a few tens of meters, while the slant range distance R for a space-borne platform is the hundreds of thousands of meters. In this scenario, the parameter q in (9) is the height with respect to the reference DEM and since k z is very small, the resulting error ψ topo will be small.

Phase Estimation via Phase Linking
The main algorithm behind the proposed method is the so-called Phase Linking [31]. The basic approach is simple: let us suppose that there is a stack of N focused and coregistered SAR images.
With such a number of images, it is possible to form N(N-1)/2 independent interferograms. Each interferometric phase is modeled exactly as in Equation (5). This set of phases form and overdetermined systems where the unknown are the best N-1 phases that fit with minimum error the observations (the original N(N-1)/2 interferometric phases).
The best set of N-1 phases is found by maximizing with respect to ψ [30]: A closed-form solution is not available, and a stable solution is thus found in an iterative way. It is useful to notice that, with N images, we can form a complex N × N covariance matrix C(ψ). If we normalize the powers on the diagonal, we obtain the so-called coherence matrix.
The cell of the coherence matrix at row i and column j now contains a complex coherence: the phase of this number is the (averaged) interferometric phase between image i and j while the absolute value is a number between 0 and 1 expressing the quality of the interferometric phase itself as already explained in Section 1.
Without loss of consistency from now on, we will refer to C as the coherence matrix and not the covariance matrix.
It is important to access the expected performances of this method in comparison with other techniques used in literature.
We performed a Monte Carlo simulation using three different estimators: (1). The phase linking estimator just explained.
(2). The AR(1) estimator that consists in integrating the phases of interferograms formed using consecutive acquisitions:ψ where C(i, i + 1) is the complex element inside the sample coherence matrix at row i a column i + 1. This method takes advantage of the high coherence of interferograms formed with consecutive acquisitions. It is also possible to show that, if in equation (7) we force γ 0 = 1, the AR(1) estimator is the exact estimator (in Maximum Likelihood sense).
The classic DInSAR estimator consists of evaluating interferograms at increasing temporal baselines.ψ whereĈ(1, 1 . . . N) denotes the first row and all the columns of the sample coherence matrix.
In Figures 3 and 4 two different scenarios, represented by the sample coherence matrix,Ĉ are simulated. For all the 10.000 Monte Carlo runs we used 100 independent looks to estimate the coherence matrix. The temporal baseline between two acquisitions is set to 6 days in order to mimic the Sentinel-1 repeated geometry interval that will be used in the result section of this paper.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 26 acquisitions. It is also possible to show that, if in equation (7) we force = , the AR(1) estimator is the exact estimator (in Maximum Likelihood sense).
The classic DInSAR estimator consists of evaluating interferograms at increasing temporal baselines.
Where ̂( , … ) denotes the first row and all the columns of the sample coherence matrix. In Figure 3 and Figure 4 two different scenarios, represented by the sample coherence matrix, ̂ are simulated. For all the 10.000 Monte Carlo runs we used 100 independent looks to estimate the coherence matrix. The temporal baseline between two acquisitions is set to 6 days in order to mimic the Sentinel-1 repeated geometry interval that will be used in the result section of this paper.  In the case of Figure 3a, a PS-like target is generated where the target remains stable for all the 10 images with a coherence γ = . .
The AR(1) method is badly performing due to the accumulation of the noise in the integration acquisitions. It is also possible to show that, if in equation (7) we force = , the AR(1) estimator is the exact estimator (in Maximum Likelihood sense).
The classic DInSAR estimator consists of evaluating interferograms at increasing temporal baselines.
Where ̂( , … ) denotes the first row and all the columns of the sample coherence matrix.
In Figure 3 and Figure 4 two different scenarios, represented by the sample coherence matrix, ̂ are simulated. For all the 10.000 Monte Carlo runs we used 100 independent looks to estimate the coherence matrix. The temporal baseline between two acquisitions is set to 6 days in order to mimic the Sentinel-1 repeated geometry interval that will be used in the result section of this paper.  In the case of Figure 3a, a PS-like target is generated where the target remains stable for all the 10 images with a coherence γ = . .
The AR(1) method is badly performing due to the accumulation of the noise in the integration step. The phase variance increases linearly with the temporal baseline. The DInSAR and Phase In the case of Figure 3a, a PS-like target is generated where the target remains stable for all the 10 images with a coherence γ = 0.6.
The AR(1) method is badly performing due to the accumulation of the noise in the integration step. The phase variance increases linearly with the temporal baseline. The DInSAR and Phase Linking estimator performs similarly and in accordance with the decorrelation model.
In the case of Figure 4a, a target showing an exponential decorrelation is generated with a decay constant of 35 days (i.e., after 35 days the coherence is about 0.37) [18].
In this scenario, the DInSAR estimator is poorly performing due to the loss of correlation between interferograms at increasing temporal baselines. On the other end, the AR(1) and Phase Linking estimators are performing similarly and with good performances.
It can be said that the phase linking can automatically adapt to the scenario (i.e., the coherence structure of the target) and thus returning the best minimum variance estimate of the interferometric phases.

Phase Linking for APS Estimation
Phase Linking can be applied on any interferometric stack gathered at any arbitrary carrier frequency and with any repetition interval [31][32][33][34], but some practical precautions are needed in order to apply this algorithm for atmospheric phase screen retrieval: (1). The total temporal extent of the stack must be carefully tuned. This precaution has two main purposes: a.
It reduces the effects of subsidences on the interferometric phase. The model of the interferometric phase in Equation (5) can present another term due to linear displacement in the line of sight direction that is equal to: where ∆t = t 2 − t 1 . In the presence of a normal subsidence rate in the order of 10 mm/year (i.e., not the one induced by an earthquake or a seismic event), if the stack is kept short (let us say eight images with a temporal extent of about 50 days), the error will be of less than 2 mm that is tolerable for our purposes. b.
Following the decorrelation model explained in Section 3, we can say that the stack temporal extent needs to take into consideration the average "life" of a distributed scatterer. With phase linking, we form all the possible interferograms with N images and from them, we estimate N-1 phases, if the coherence of the interferograms with very long temporal baseline is very low, they will bring noise into the final estimate. A solution is then to reduce the maximum temporal baseline by considering the decorrelation model. It is useful to remember that the decorrelation time depends on the wavelength used for the measure: In [18], Rocca made the example of 40 days in C-Band, but the reasoning can be easily extended in L or P Band where the average decorrelation time is much higher and thus a larger dataset can be used.
(2). The estimation window must be carefully tuned. The quality of the interferometric phase is tightly related to the number of looks used to estimate it. In [23], the expression of the variance on the estimate of the interferometric phase is given for a single interferometric couple: It is interesting to notice how the quality depends on the coherence and on the number of looks L used for estimating it. A distributed scatterer with low coherence can indeed generate a good phase estimate provided that a sufficient number of looks is used.
In the case of the phase linking procedure, where a joint stack of N images is used, the formula for the interferometric phase variance is different and it is directly derived by the expression of the Cramer-Rao Bound. In particular, the minimum attainable variance is given by the diagonal of the inverse of the Fisher Information Matrix (FIM): where χ is the N × N FIM expressed as [30]: where C is the coherence matrix, I N is the identity matrix of size N × N and L is again the number of looks. Note that now the variance is a vector of size N.
In order to evaluate the minimum number of looks required, we model that the matrix C to show an exponential decay as in Figure 4a with a decorrelation constant of 40 days and a sampling interval of 6 days as for the C-Band Sentinel-1 constellation. The required standard deviation is set to be 0.5 mm: this requirement is set for the last phase in the stack, namely the one that is noisier among all the estimated phases (as shown in Figure 4b).
With these constraints, we obtain a number of looks L = 250 ( Figure 5) where is the coherence matrix, is the identity matrix of size × and L is again the number of looks. Note that now the variance is a vector of size .
In order to evaluate the minimum number of looks required, we model that the matrix to show an exponential decay as in Figure 4a with a decorrelation constant of 40 days and a sampling interval of 6 days as for the C-Band Sentinel-1 constellation. The required standard deviation is set to be 0.5 mm: this requirement is set for the last phase in the stack, namely the one that is noisier among all the estimated phases (as shown in Figure 4b).
With these constraints, we obtain a number of looks = ( Figure 5). Sentinel-1 shows a nominal resolution of × in IW acquisition mode in range and azimuth respectively. The images are properly oversampled to about . × . In order to obtain a number of independent looks of 250 we need roughly 1250 pixels: a × window can be used spanning an area of × . Figure 5. Standard deviation of the estimated phase versus the number of independent looks used for the estimate. For this processing we used a covariance matrix modelled as in (7), with 0 = 0.7, = 0.975, = 6 as suggested in [8].
In our studies carried over several datasets with different location and scatterer characteristics and with temporal baseline of 6 days, we found that a good estimation window is in the order of × : big enough for allowing reliable evaluation of the APS, but small enough to comply with OSCAR requirements for spatial resolution. The discrepancy with respect to the theoretical window size could be explained by the fact that, in some places, the data is not respecting the decorrelation model and indeed the decorrelation time is much lower than the nominal 40 days predicted. This condition calls for a bigger estimation window to compensate for the poor coherence.
In this contribution, we decided to use the twin Satellites Sentinel-1 A/B: their C-Band (5.40 GHz) payload, together with the 6 days repeated geometry and closely spaced and accurate orbits make the constellation the perfect instrument for the generation of differential SAR APSs.
A very interesting perspective for the future is the possibility to have a geostationary or geosynchronous platform that shows a repetition interval in the order of hours and not days: this will heavily limit temporal decorrelation giving also the possibility to map troposphere with an unprecedent spatial and temporal resolution. Figure 5. Standard deviation of the estimated phase versus the number of independent looks used for the estimate. For this processing we used a covariance matrix modelled as in (7), with γ 0 = 0.7, ρ = 0.975, T = 6 as suggested in [8].

Phase Unwrapping
In our studies carried over several datasets with different location and scatterer characteristics and with temporal baseline of 6 days, we found that a good estimation window is in the order of 350 × 350 m: big enough for allowing reliable evaluation of the APS, but small enough to comply with OSCAR requirements for spatial resolution. The discrepancy with respect to the theoretical window size could be explained by the fact that, in some places, the data is not respecting the decorrelation model and indeed the decorrelation time is much lower than the nominal 40 days predicted. This condition calls for a bigger estimation window to compensate for the poor coherence.
In this contribution, we decided to use the twin Satellites Sentinel-1 A/B: their C-Band (5.40 GHz) payload, together with the 6 days repeated geometry and closely spaced and accurate orbits make the constellation the perfect instrument for the generation of differential SAR APSs.
A very interesting perspective for the future is the possibility to have a geostationary or geosynchronous platform that shows a repetition interval in the order of hours and not days: this will heavily limit temporal decorrelation giving also the possibility to map troposphere with an unprecedent spatial and temporal resolution.

Phase Unwrapping
Once the phases are estimated properly, they are still only of known modulo 2π (i.e., they are wrapped). In order to unwrap the phases, we rely on the fact that the spatial variation of such phase field is smooth from resolution cell to resolution cell with respect to the ambiguity, which is defined as the length that makes the phase wrap. This measure, in a SAR framework, is half the operational wavelength of the system: in the case of Sentinel-1, it is 2.8 cm. The advantage of using a DS processing for APS retrieval is that the phase field is spatially smooth and can be easily unwrapped without doing differences (arcs) in the space domain that will lead to noise accumulation after re-integration. The operation of phase unwrapping is scientifically mature and it is out of the scope of this contribution.

Orbit Correction: GNSS Processing and Integration
In order to properly compensate an interferogram for flat earth contribution and for topography, orbit state vectors must be known precisely, otherwise, a residual phase is present in the interferometric phase.
This error is usually neglected if the orbits are known with very high precision or if the area of interest is small. Requirements for orbit accuracy are described in Section 2.3 and in [35].
Since the 3D precision of the positioning of Sentinel-1 is in the order of 5 cm RMS and since we want to produce wide-area APS maps, we must consider and correct those errors to prevent the generation of low-frequency errors in the estimated atmospheric product.
An approach used in literature is to model as a plane the phase induced on the interferogram by an error in the orbit: this solution is not feasible in our circumstances since the signal of interest (the APS) may be modeled in the exact same way and the result would be to invalidate the atmospheric map generated. We rely instead on a network of GNSS stations to properly compensate for orbit error.
The meteorological applications of GNSS are well known [36,37]. A GNSS receiver is able to determine its coordinates in a global reference frame by using simultaneous observation of its distance from a number of satellites of known position. The distances are derived from the time it takes to the GNSS signals to cover the satellite-receiver paths and by assuming that they travel with the constant velocity of light in vacuum.
This assumption results in an observed distance different from the geometrical one by an amount called atmospheric delay. The delay contains an ionospheric component, which can be removed by a proper combination of the dual-frequency GNSS signals and a tropospheric component.
The tropospheric delay has to be accounted for in the GNSS data processing to get accurate positioning.
To this aim, the tropospheric delays affecting the signals received by a GNSS station from all the satellites simultaneously in view are expressed as a common delay in the zenith direction above the receiver. That is, each slant delay is projected onto the zenith direction with a known mapping function. The common zenith projection is then estimated in the general adjustment of GNSS observations, for instance by modeling it as a random walk process, resulting in a high temporal resolution time series of zenith delays for each considered station. In this work, the free and open-source software goGPS [38] was used to estimate ZTD time series from the GNSS raw observations.
The procedure for obit correction that exploits both SAR and GNSS observations is explained in the following.
The absolute interferometric phase can be expressed as in (5) as: where B (P, t 1 , t 2 ) is the parallel baseline as seen from the point P on ground and between two acquisitions at time t 1 and t 2 . From now on, the pixel P will be associated with a looking angle θ(off-nadir) and a zero-Doppler acquisition time t.
In order to access the sensibility of the interferometric phase to the orbit error, we differentiate Equation (16) leading to: where now B ε, (P, t 1 , t 2 ) is the error of parallel baseline and ψ ε (P, t 1 , t 2 ) is the phase error induced by the orbital error. Expanding Equation (16) in Taylor series with respect to θ and t and stopping at the first term we obtain: This is the equation of a plane in the time-angle domain where the parameters to be estimated are dB ε, (θ 0 ,t 0 ) dt and B ε,⊥ (θ 0 , t 0 ), namely the variation of the parallel baseline with time and the normal baseline. In other words, we first estimate the APS with SAR and with GNSS stations in the location of the stations themselves. We then remove the APSs generated with GNSS from the APSs generated with SAR. The only remaining signal (apart from thermal noise and clutter) is the orbital error. Now we robustly estimate (in L1 norm) the parameters of Equation (17), we project the estimated parameters over the whole scene by evaluating Equation (17) over the whole grid. What we find is a plane that can be corrected from the unwrapped interferogram.
As a final note, we can also state that the SAR measurements can be used to qualify the GNSS data by identifying outliers. It is sufficient, in fact, to access the residuals after the inversion of (17): where ψ ε,d (θ, t) is the SAR data with the APSs removed using GNSS data and ψ ε (θ, t) is the reconstructed data with the estimated parameters. If the residual of a station at a particular time instant is high, it means that the data of that GNSS station was not useful for the inversion (and thus discarded by robust inversion). Numeric results are given in the next section of this paper.

Dataset
In this section, we report the results obtained with data from the European Space Agency (ESA) satellite constellation Sentinel-1 A/B.
As already discussed, the unprecedent repetition frequency over the same scene and the C-Band payload make the Sentinel-1 constellation the perfect instrument for extracting atmospheric products.
While several experiments have been carried out over different regions with different orographic conditions, in this section, we report the experiment carried out over northern Italy.
Other locations where experiments have been performed are central Italy, southern Italy, and Southern Africa.
The full dataset is represented in the Table 3 and it is composed of 7 IW Swath images gathered over northern Italy. The total time span is 42 days in order to prevent high temporal decorrelation and subsidences bias in the estimate as explained in the methodology section of this paper.

Processing
While all the processing has been carried out full swath, thus on an area of about 42.000 km 2 , the comparison with other techniques and datasets in the next subsections is carried out only a small subset of about 10,000 km 2 .
In Figure 6 an example of a single full swath atmospheric phase screen is proposed. The missing portions are on water bodies where there is not enough coherence to retrieve a reliable interferometric phase.

Processing
While all the processing has been carried out full swath, thus on an area of about 42.000 km 2 , the comparison with other techniques and datasets in the next subsections is carried out only a small subset of about 10,000 km 2 .
In Figure 6 an example of a single full swath atmospheric phase screen is proposed. The missing portions are on water bodies where there is not enough coherence to retrieve a reliable interferometric phase. In Figure 7a all the interferometric phases are shown: the interferogram between image and image is represented in the row and column. In Figure 7b, instead, the Phase-Linked interferometric phases are shown. This matrix is formed as the outer product of the vector found by the linking procedure : ̂= (20) It can be shown that this is the best rank-1 approximation of the coherence matrix . In Figure 7a all the interferometric phases are shown: the interferogram between image k and image l is represented in the k th row and l th column. In Figure 7b, instead, the Phase-Linked interferometric phases are shown. This matrix is formed as the outer product of the vector found by the linking procedure: It can be shown that this is the best rank-1 approximation of the coherence matrix C.  In Figure 7c, the residual between the two images is presented: it is clear by looking at the cells far away from the diagonal, that, increasing the temporal baseline, the standard interferograms becomes noisy, while the phase-linked reconstruction preserves spatial smoothness. This is even clearer in Figure 8 where the standard interferograms with increasing temporal baseline are shown in the top row, while the phase-linked interferograms are in the middle row and the residual between the two is in the bottom row. The areas with no residuals (light blue in the bottom row) correspond to cities in the scene where the temporal decorrelation is not significant since the vast majority of targets can be represented as Permanent Scatterers (PS) and thus the phase linking method and the PS-Like one perform similarly. In Figure 7c, the residual between the two images is presented: it is clear by looking at the cells far away from the diagonal, that, increasing the temporal baseline, the standard interferograms becomes noisy, while the phase-linked reconstruction preserves spatial smoothness. This is even clearer in Figure 8 where the standard interferograms with increasing temporal baseline are shown in the top row, while the phase-linked interferograms are in the middle row and the residual between the two is in the bottom row. The areas with no residuals (light blue in the bottom row) correspond to cities in the scene where the temporal decorrelation is not significant since the vast majority of targets can be represented as Permanent Scatterers (PS) and thus the phase linking method and the PS-Like one perform similarly. On the other end, on vegetated areas, the phase-linked interferograms show less noise as predicted in the methodology part of this paper (Figure 3 and Figure 4).
In the bottom row of Figure 8 spatial variograms at short distances are shown: while the formal definition of variogram will be carried out in Section 5.3, it is sufficient to notice that the standard DInSAR processing shows higher semi-variance with respect to Phase Linking proving that the interferograms are indeed noisier (or equivalently that the phase linked phases are smoother).

Orbit Correction
As described in Section 3 and Section 4, an orbit error compensation is mandatory if the objective of the processing is to derive wide-area APS maps. A linear phase trend in range and azimuth is present if the orbit shows an error in the normal baseline and/or in the derivative w.r.t. time of the parallel baseline (Figure 1). This processing step aims at retrieving those errors from the interferometric phase itself provided that the only contribution the phase is the orbit error itself (and obviously the noise). On the other end, on vegetated areas, the phase-linked interferograms show less noise as predicted in the methodology part of this paper (Figures 3 and 4).
In the bottom row of Figure 8 spatial variograms at short distances are shown: while the formal definition of variogram will be carried out in Section 5.3, it is sufficient to notice that the standard DInSAR processing shows higher semi-variance with respect to Phase Linking proving that the interferograms are indeed noisier (or equivalently that the phase linked phases are smoother).

Orbit Correction
As described in Section 2.3, an orbit error compensation is mandatory if the objective of the processing is to derive wide-area APS maps. A linear phase trend in range and azimuth is present if the orbit shows an error in the normal baseline and/or in the derivative w.r.t. time of the parallel baseline ( Figure 1). This processing step aims at retrieving those errors from the interferometric phase itself provided that the only contribution the phase is the orbit error itself (and obviously the noise).
In order to achieve this quality of the interferometric phase, we used a network of GNSS stations in order to derive absolute Zenith Total Delay and then we differentiated w.r.t the ZTD of the GNSS at the acquisition time of the SAR master image. In Figure 9 the network of GNSS stations used for the case study is shown.
In order to achieve this quality of the interferometric phase, we used a network of GNSS stations in order to derive absolute Zenith Total Delay and then we differentiated w.r.t the ZTD of the GNSS at the acquisition time of the SAR master image. In Figure 9 the network of GNSS stations used for the case study is shown. The estimation is carried out in a robust L1 norm in order to obtain a reliable estimation from noisy (and possibly with outliers) data.
For the dataset of interest (in Table 3) the estimated parameters are shown in Figure 10 it is important to notice that, since we work in an interferometric framework, the derived error is the sum of the individual errors of the master and slave images. We can then decide to attribute all the errors to the slave, to the master, or half to the master and half to the slave. In our case, we decided to adopt the first approach. The values represented in Figure 10 are consistent with the nominal values for orbit accuracy: Precise Orbital Products (POD) allows for an accuracy in the 3D orbit determination of 5 cm RMS (1 sigma). The estimation is carried out in a robust L1 norm in order to obtain a reliable estimation from noisy (and possibly with outliers) data.
For the dataset of interest (in Table 3) the estimated parameters are shown in Figure 10 it is important to notice that, since we work in an interferometric framework, the derived error is the sum of the individual errors of the master and slave images. We can then decide to attribute all the errors to the slave, to the master, or half to the master and half to the slave. In our case, we decided to adopt the first approach. The values represented in Figure 10 are consistent with the nominal values for orbit accuracy: Precise Orbital Products (POD) allows for an accuracy in the 3D orbit determination of 5 cm RMS (1 sigma).
In Figure 11, instead, the residuals for each GNSS station is presented: it is clear that the average residual is quite low as expected, so the model used fits well the data. There are just two exceptions represented by the stations with latitude and longitude equals to (45.6892, 9.61201) and (45.4447, 11.0024). This anomaly can be either an error on the processing of GNSS-derived data, an area with poor coherence in the SAR data or malfunctioning of a GNSS station.
After the estimation of the parameters using only the points with a GNSS station, the direct problem is solved for all the points in the grid and an orbital plane is derived. A profile along the direction of maximum slope of the orbital plane error is shown in Figure 12a where it is clear that, after the removal of the orbit error, the profile has no more a linear trend with negative slope. In Figure 11, instead, the residuals for each GNSS station is presented: it is clear that the average residual is quite low as expected, so the model used fits well the data. There are just two exceptions represented by the stations with latitude and longitude equals to (45.6892, 9.61201) and (45.4447, 11.0024). This anomaly can be either an error on the processing of GNSS-derived data, an area with poor coherence in the SAR data or malfunctioning of a GNSS station. After the estimation of the parameters using only the points with a GNSS station, the direct problem is solved for all the points in the grid and an orbital plane is derived. A profile along the direction of maximum slope of the orbital plane error is shown in Figure 12a where it is clear that, after the removal of the orbit error, the profile has no more a linear trend with negative slope.  In Figure 11, instead, the residuals for each GNSS station is presented: it is clear that the average residual is quite low as expected, so the model used fits well the data. There are just two exceptions represented by the stations with latitude and longitude equals to (45.6892, 9.61201) and (45.4447, 11.0024). This anomaly can be either an error on the processing of GNSS-derived data, an area with poor coherence in the SAR data or malfunctioning of a GNSS station. After the estimation of the parameters using only the points with a GNSS station, the direct problem is solved for all the points in the grid and an orbital plane is derived. A profile along the direction of maximum slope of the orbital plane error is shown in Figure 12a where it is clear that, after the removal of the orbit error, the profile has no more a linear trend with negative slope.

Variograms and Radially Average Spectra
In order to statistically characterize the derived atmospheric maps, spatial variograms and spectra have been computed. In a first approximation, the total slant range delay is the sum of

Variograms and Radially Average Spectra
In order to statistically characterize the derived atmospheric maps, spatial variograms and spectra have been computed. In a first approximation, the total slant range delay is the sum of hydrostatic and wet delay. The first one is usually much smoother than the latter and it usually manifests as a trend with weak slope [28,39].
The wet delay, on the other end, has a magnitude that is much smaller than the hydrostatic delay (typically less than 10% of the total slant delay), but its fluctuations are larger both in space and time. The wet delay is responsible for the temporal and spatial decorrelation between APS maps.
The wet delay is usually statistically modeled using a variogram [40]: where P is a location in the APS map, ∆P is a displacement w.r.t. P, K is related to the strength of the turbulence and ρ is related to the smoothness. Turbolence-induced delay has a scale-variant power law behavior with a transition region at approximately 1.5 km. Below this scale a proper value for ρ is 5/3 while above it is approximately 2/3 [41,42]. In Figure 13 the variograms computed from the data are compared with the theoretical model. For this computation, we assumed that the APS is isotropic. A good fitting is shown in the average of all the variograms derived for every single data. The same can be said about the radially averaged spectra where we assumed again that the underneath process is isotropic. In this case, we computed the 2D Fourier transform of each APS map and averaged radially with respect to the center of the transformed domain.
In Figure 14 the results from the derived APS maps are shown: the derived spectra fit really well the theoretical value provided for the turbulent part of the APS by Hanssen [28] and Tatarski [40] as: where = / for radially averaged power spectra. Even in this case, good accordance with the theoretical model is shown by the data. A good fitting is shown in the average of all the variograms derived for every single data. The same can be said about the radially averaged spectra where we assumed again that the underneath process is isotropic. In this case, we computed the 2D Fourier transform of each APS map and averaged radially with respect to the center of the transformed domain.
In Figure 14 the results from the derived APS maps are shown: the derived spectra fit really well the theoretical value provided for the turbulent part of the APS by Hanssen [28] and Tatarski [40] as: where α = 5/3 for radially averaged power spectra. Even in this case, good accordance with the theoretical model is shown by the data.  Table 3. (blue line) mean of all the previous spectra. (red dashed line) theoretical spectrum computed with α = 5/3.

Comparison with Reference APS Maps from SqueeSAR®
In the methodology section of this paper, we explained that in order to prevent the effect of subsidences in the estimated phases we kept the total temporal length of the stack quite small. Subsidences, in fact, usually are in the order of a few mm/year and if you kept the stack very large (for example 1 year of measurements) we had to estimate and remove the subsidences, while by keeping the total stack extend to a maximum of 40 days we can completely discard this effect (in normal situations, not in the case of earthquake or big terrain velocities).
In order to access this assumption, we confronted the APS maps derived with our methodology with a stack of 46 APSs computed with SqueeSAR® [21].
SqueeSAR® is an advanced InSAR technique mainly applied in deformation estimation and displacement monitoring. The main algorithm involves, as the already presented Phase Linking, the estimation of N-1 optimal phases from the N(N-1)/2 observations. The workflow is however composed by several steps: the processing starts with an identification of statistically homogeneous pixels (SHPs) by means of the Kolmorogov-Smirnov test. This test is applicable with high reliability to data-stacks with as little as 8 images.
When the SHPs have been identified it performs a space adaptive filtering to generate sample coherence matrices over which the optimal phases are estimated.
Once the optimal solution has been obtained, the quality of the estimated phase is assessed: if the "goodness of fit" is over a certain threshold, the original phases are substituted with their optimized values.
SqueeSAR® includes a standard Permanent Scatterer (PS) processor that allows for the removal of the portion of the interferometric phase related to subsidences and residual DEM: in order to  Table 3. (blue line) mean of all the previous spectra. (red dashed line) theoretical spectrum computed with α = 5/3.

Comparison with Reference APS Maps from SqueeSAR®
In the methodology section of this paper, we explained that in order to prevent the effect of subsidences in the estimated phases we kept the total temporal length of the stack quite small. Subsidences, in fact, usually are in the order of a few mm/year and if you kept the stack very large (for example 1 year of measurements) we had to estimate and remove the subsidences, while by keeping the total stack extend to a maximum of 40 days we can completely discard this effect (in normal situations, not in the case of earthquake or big terrain velocities).
In order to access this assumption, we confronted the APS maps derived with our methodology with a stack of 46 APSs computed with SqueeSAR® [21].
SqueeSAR®is an advanced InSAR technique mainly applied in deformation estimation and displacement monitoring. The main algorithm involves, as the already presented Phase Linking, the estimation of N-1 optimal phases from the N(N-1)/2 observations. The workflow is however composed by several steps: the processing starts with an identification of statistically homogeneous pixels (SHPs) by means of the Kolmorogov-Smirnov test. This test is applicable with high reliability to data-stacks with as little as 8 images.
When the SHPs have been identified it performs a space adaptive filtering to generate sample coherence matrices over which the optimal phases are estimated.
Once the optimal solution has been obtained, the quality of the estimated phase is assessed: if the "goodness of fit" is over a certain threshold, the original phases are substituted with their optimized values.
SqueeSAR®includes a standard Permanent Scatterer (PS) processor that allows for the removal of the portion of the interferometric phase related to subsidences and residual DEM: in order to estimate and compensate for this effect in an accurate and reliable way it is mandatory exploit a stack of tens of images resulting in high memory usage.
When all the unwanted contribution to the interferometric phase has been estimated and compensated, the only relevant component that is left in the interferometric phase is the APS. A statistical interpolation can now be performed to fit those areas showing poor quality for the optimal estimated phase.
This processing pipeline is surely effective, but it involves several steps that are generally very time consuming if performed over a large spatial scale such as the one of interest for APS estimation.
The presented method, on the other end, exploits statistical characteristics of the atmospheric signal, performing filtering over large spatial windows to retain the only portion of the interferometric phase that plays a role in such large scales: an SHP test is thus not required. It also does not involve the estimation and removal of subsidences since the largest temporal baseline in the stack is kept small reducing in this way unwanted effects due to linear deformation of the terrain.
While SqueeSAR®provides accurate results at the expense of a large computational complexity, the simple approach proposed does so at the expenses of using a small data stack leading in this way to a much shorter time series of APSs.
The proposed algorithm, as already explained, concerns also the integration with GNSS-derived atmospheric products for orbit correction. Future developments of the technique will also allow for the extraction of absolute water vapor maps.
Since the SqueeSAR®dataset was provided without ionospheric compensation, we skipped this step even in our processing to make them comparable.
The comparison algorithm runs as follows: the SqueeSAR®derived APS are taken in a subset corresponding to the APS derived with our methodology, then we reference all the phases to a common (unknown) master. As final steps, we resample both maps on a common grid and we generate residual maps.
In Figure 15 an example of comparison is shown. For visualization simplicity, we decided to omit all the maps with the comparison and instead synthesized the results in the Table 4.
From the image and the numerical results in the Table 4 is clear that the two APSs computed are qualitatively and quantitatively comparable. The mean standard deviation is well below 3 mm for all the processed interferograms.
This confirms the hypothesis that, by using a limited total temporal extent of the stack, not only we limit the effects of temporal decorrelation, but we limit also the effect of subsidences that have power is the same spatial scales of the atmosphere.  This confirms the hypothesis that, by using a limited total temporal extent of the stack, not only we limit the effects of temporal decorrelation, but we limit also the effect of subsidences that have power is the same spatial scales of the atmosphere.

A Note about GNSS and NWMP Comparison and NWMP Ingestion of SAR-Derived APS
A comparison between SAR-derived atmospheric phase screen and GNSS derived atmospheric phases can be done easily by simply differentiating the Zenith Total Delay of the GNSS in the date of the slave of the SAR acquisitions with respect to the master of the SAR stack.
In this paper, however, we decided to avoid this comparison since all the GNSS stations were used for the orbit error correction. In the presence of a larger number of stations, it is indeed possible to divide them into two subsets using one for orbit correction and one as a validation set.
In order to allow the ingestion of SAR-derived delay maps into NWPMs, an absolute phase extraction is necessary.
All the processing done up to now deals with differential maps, thus an unknown "master" ZTD maps needs to be derived. While this argument is not the topic of this paper it is useful to give a possible solution to the problem and hints for a future research topic. The extraction of the absolute phases from a set of differences (interferometric phases) is an undetermined problem with infinite solutions. The null space of the problem is a constant value that is missing and impossible to estimate from the given data. A possible solution is to have a-priori information for the extraction that can be

A Note about GNSS and NWMP Comparison and NWMP Ingestion of SAR-Derived APS
A comparison between SAR-derived atmospheric phase screen and GNSS derived atmospheric phases can be done easily by simply differentiating the Zenith Total Delay of the GNSS in the date of the slave of the SAR acquisitions with respect to the master of the SAR stack.
In this paper, however, we decided to avoid this comparison since all the GNSS stations were used for the orbit error correction. In the presence of a larger number of stations, it is indeed possible to divide them into two subsets using one for orbit correction and one as a validation set.
In order to allow the ingestion of SAR-derived delay maps into NWPMs, an absolute phase extraction is necessary.
All the processing done up to now deals with differential maps, thus an unknown "master" ZTD maps needs to be derived. While this argument is not the topic of this paper it is useful to give a possible solution to the problem and hints for a future research topic. The extraction of the absolute phases from a set of differences (interferometric phases) is an undetermined problem with infinite solutions. The null space of the problem is a constant value that is missing and impossible to estimate from the given data. A possible solution is to have a-priori information for the extraction that can be used to regularize the problem. This a-priori information can be the absolute ZTD phases from the NWPM itself or pointwise ZTD from a GNSS network.
A possible cost function to minimize can be: J(φ SAR ) = φ SAR − φ NWPM 2 2 +λ T (ψ SAR − Hφ SAR ) (23) where φ SAR is the unknown vector of N phases to be estimated, φ NWPM is the known a priori vector of absolute ZTD from GNSS measurements of NWPM run, λ is the Lagrange multiplier, ψ SAR is the vector containing the N-1 phase linked atmospheric differential phases (the data) and H is the design matrix that generates all the differences of absolute phases. The minimization of this cost function allows reconstructing the absolute phases that are most similar to the ones that are given as a-priori, without losing data fidelity thanks to a constrained optimization.

Conclusions
This paper proposes a technique aimed at estimating in a fast, robust, and reliable way atmospheric phase screens from a stack of focused and coregistered SAR images, for the ingestion into Numerical Weather Prediction Models (NWPMs).
The requirements, in terms of spatial resolution for the ingestion of such APS maps into NWPMs, imply the exploitation of both Permanent and Distributed Scatterers in order to have wide and dense area coverage.
A target characterization has been proposed and, together with the theoretical performances of the estimator, a minimum number of looks required for an accurate estimation have been suggested.
In this processing, GNSS-derived atmospheric products have been used as an a-priori information to estimate and compensate the orbital plane induced by an error in the orbit determination: this correction is mandatory when we deal with wide-area APS estimation since an orbital plane may create a phase trend that is indeed an error in the low spatial frequencies of the derived map.
A case study and comparison over northern Italy has been performed also with a stack of SqueeSAR®derived APS: these maps are computed from a dataset composed of 46 images where the effect of subsidences has been removed, while the Phase Linked maps exploit just seven images. The numerical comparison confirms the hypothesis that, even if the stack processed with the proposed method is smaller and the effect of subsidences has not been removed, the quality of the APS map is not compromised.