Pre-analysis of GNSS tomography solution using the concept of spread of model resolution matrix

GNSS tropospheric tomography is one of the applications of the Global Navigation Satellite Systems (GNSS) signals which attracts more and more interest in the field of meteorology. This method can reconstruct the water vapour of the atmosphere, which has a considerable effect on weather forecasting and early warning systems of severe weather. In GNSS tomography, traditionally, a regular spaced 3D grid stretches from the GNSS network to the effective height of the troposphere in the area of interest. Therefore, this method can offer a permanent monitoring service for water vapour and wet refractivity fields at a low cost and a reasonable spatial resolution compared to conventional observations, like radiosonde and radio occultation profiles. Nevertheless, the quality of the reconstructed field is still one of the challenges in the GNSS tomography. In this research, we propose the concept of spread as a mathematical tool to provide a quality measure without using the reference field and calculating statistical measures like RMSE and Bias. Thereby, two synthetic and one real datasets (part of Germany and Czechia) covering overlapping periods between 29 May and 14 Jun of the year 2013 (DoY 149–165; DoY 160–165; DoY 160–165, 2013) have been tested to investigate the proposed method. According to the obtained results, the proposed tool shows a strong correlation (up to 0.81 for synthetic and 0.72 for real observations) with the standard deviation of the reconstructed wet refractivity with respect to the radiosonde profile reference. The correlation between spread and the Bias of the retrieved wet refractivity field is also significant. However, there is no clear picture depending on the applied spread computation models. Therefore, the spread of the resolution matrix can be used as a proxy for the accuracy of the tomography reconstruction field based on the quality of the observations, the initial field, and the design matrix.

atmospheric layers. One of the most effective layers during the transmission of GNSS signals to a ground receiver is the troposphere, which causes a bending impact on the satellite's path because of the propagation properties of electromagnetic waves. Moreover, due to the physical characteristics of the troposphere, the refraction index in this layer highly depends on temperature, pressure, and water vapour content. Therefore, modelling the refraction index is complicated, especially concerning the wet component of the troposphere, because of dynamic changes of tropospheric conditions (Seeber 1993). Consequently, the tropospheric error in the GNSS signals can be considered as a valuable data source for meteorological studies by an estimation of parameters such as integrated water vapour (IWV) or precipitable water vapour (PWV) (Bevis et al. 1992;Boccolari et al. 2002;Chiarabba and Amato 1997;Emardson et al. 1998;Heublein 2019;Paziak 2012;Rocken et al. 1997;Toomey and Foulger 1989). Nowadays, GNSS tomography is one of the most promising methods to reconstruct the spatio-temporal structure of the troposphere. In this method, a three-dimensional (3D) grid in horizontal and vertical directions stretches from the ground based GNSS network of interest to an effective tropospheric height. However, as the amount of GNSS rays is insufficient at a single epoch, a tomography time window should be determined to accumulate observations within a specific period. Then, using the inversion techniques, an average of the refractivity field caused by the water vapour structure is reconstructed according to the 3D voxel-based method over the defined tomography window. In recent years, numerous studies have proved the potential of GNSS tomography to reconstruct the tropospheric structure Benevides et al. 2018;Flores et al. 2000a, b;Guo et al. 2016;Heublein et al. 2019;Nilsson et al. 2007b;Rohm and Bosy 2009;Trzcina and Rohm 2019). Besides, some tomography models have been set up in different countries to estimate the 3D wet refractivity or water vapour fields based on partly dense GNSS networks (Adavi and Mashhadi 2014;Brenot et al. 2018;Champollion et al. 2005;Rohm et al. 2014;Sá 2018). To evaluate the reliability of tomography solutions, independent sources like radiosonde, water vapour radiometer (WVR), and numerical weather models have been applied (Adavi et al. 2020;Bastin et al. 2007;Brenot et al. 2018;Elgered et al. 1991;Gradinarsky and Jarlemark 2004;Hanna et al. 2019;Nilsson et al. 2007a;Notarpietro et al. 2011;Troller 2004;Van Baelen et al. 2011). Nevertheless, the quality of the reconstructed wet part of the troposphere suffers from hurdles such as the unbalanced and insufficient distribution of observations deeply linked to the geometrical distribution of GNSS satellites and receivers, the applied regularisation method as well as the choice of the a priori field. Up to now, ground-based GNSS observations (e.g. Möller 2017), radiosonde profiles (e.g. Hanna et al. 2019), and numerical weather models (e.g. Brenot et al. 2020) have been used for assessing the accuracy of the reconstructed field.
One of the well-known tools to deal with inversion problems, e.g. in seismology, is the spread of the model resolution matrix to evaluate the quality of the estimated solution. This quantity measures the quality of the model parameters by considering the goodness of data, model resolution matrix, and covariance matrix. The main concept behind this value is the model space resolution matrix, which contains valuable information about the design matrix and observation quality. To the authors' best knowledge, this research is the first published study to analyse the quality of GNSS tomography model parameters using the spread of the resolution matrix. It should be highlighted that the covariance matrix of the model parameters also contains valuable information about the quality of the retrieved solution and dependency between different model parameters. This matrix is one effective parameter in the spread calculation, but there are also other effective parameters like the initial field and observation covariance matrix. Furthermore, the cofactor matrix is another important matrix that can be used for a prior analysis of the model parameters. However, both matrices are not unique indicators in predicting the overall model quality. In consequence, we hypothesise that a careful validation of the spread could provide enough information on the GNSS tomography model quality before the computation of the wet refractivity field in the inversion process.
The paper is structured as follows. Section 2 describes the theoretical basis for the GNSS tomography and then the spread of the resolution matrix in Sect. 2.1 and Sect. 2.2, respectively. Moreover, the general formula for computing the quality of the reconstructed wet refractivity is presented in Sect. 2.3. In Sect. 3, the test case, which is located in western parts of the Czech Republic and South Germany with 72 GNSS stations, is described. The numerical results of the study are presented in Sect. 4. Section 4.1 analyses the accuracy analysis of the reconstructed wet refractivity profile over all experimental periods. An assessment of spread of resolution matrix as a prior quality indicator for the tomography solution is discussed in Sect. 4.2. Finally, Sect. 5 summarizes the obtained results and associated conclusions are illustrated.

Methodology
This section first introduces the concept of GNSS tropospheric tomography of the wet refractivity. Then, the concept of the spread using the model space resolution matrix is detailed. Finally, statistical measures like RMSE and Bias, which are commonly applied to evaluate the accuracy of the tomography model, are presented.

GNSS tropospheric tomography
GPS tomography applies slant wet delays (SWDs [unit: mm]) inside a defined 3D model to retrieve the spatio-temporal Nw field [unit: ppm] for each model element (voxel) using the following functional relation (Flores, et al. 2000a, b;Heublein 2019): where d i j is the length of the ith ray inside the jth voxel [unit: km] in the east (1 : L 1 ), north (1 : L 2 ) and height (1 : L 3 ) directions. We can reformulate Eq. 1 in matrix notation as follows (Flores, et al. 2000a, b): Here, SW D is the observation vector of length m, where m is calculated from the number of GNSS stations and visible satellites in the defined time window of the tomography model. N w represents the refractivity field vector of size n where n (L 1 × L 2 × L 3 ) is the number of voxels in the tomography model. Consequently, A is a m × n matrix with the responsibility of mapping the unknown space onto the measurement space. A is called a design matrix with the following definition (Heublein 2019;Rohm and Bosy 2009): According to the geometry between ground stations and satellite constellation within the investigated time span, the design matrix A contains many zero elements. Moreover, as some voxels are not passed by any rays and some others passed by just only a few rays or a vast number of rays, Eq. 2 is mixed-determined within the defined temporal resolution (Menke 2012). Therefore, the exact solution of the tomography problem cannot be estimated directly through Eq. 2 since the inversion of the design matrix is incalculable. In other words, the null space of the design matrix is non-trivial N ( A) {0} as the singular values of the design matrix decay gradually towards zero without any noticeable gap between non-zero and zero singular values (Aster et al. 2013;Hansen 1998;Menke 2012). In order to find a unique reconstructed field of wet refractivity using Eq. 2, different constraints and also regularisation techniques should be applied. One of these methods is to consider the measurement covariance matrix (C obs ) and also the a priori covariance matrix of the unknown parameters (C m ) to reconstruct the wet refractivity field using damped least-squares inversion methods (Brenot et al. 2018;Tarantola 2005) as detailed below: with whereby N w0 is an a priori wet refractivity field and could be extracted from the numerical weather model and G −1 is the inverse of G. In this research, C m is a n × n diagonal matrix and can be expressed as (Brenot et al. 2018;Champollion et al. 2005): where δ m is a damping coefficient and defined within the range of (0 1). In order to estimate the SW D, we used the Bernese GNSS software to estimate zenith total delay (ZTD), north gradient delay (G N ) and east gradient delay (G E ) (Dach et al. 2015). Then, the zenith wet delay (ZWD) is derived by subtracting the hydrostatic part of the delay (ZHD) using the hydrostatic formulation (Saastamoinen 1973): Now, SWD can be obtained from Eq. 5 (Kačmařík et al. 2017;: Here, ε and α are the elevation and azimuth angles of the visible satellite, and M w is the non-hydrostatic mapping function which is computed by means of the Vienna Mapping Function (VMF) (Böhm et al. 2006;Landskron and Böhm 2018). Moreover, a gradient mapping function M G is calculated as follows (Chen and Herring 1997):

The concept of spread of resolution matrix
One of the fundamental challenges in the tomography technique is to appraise the quality of the reconstructed model parameters. According to the solution quality, we can identify the regions, which are fairly-well resolved or unresolved. One of the most common methods in this regard is the assessment of the model resolution matrix R m . Using this valuable mathematical tool, the quality of a generalised inverse solution is characterised (Aster et al. 2013;Menke 2012). To identify the model resolution matrix, we imagine the true model solution, here N true w , and the estimated solution, here N est w . The true model parameters are determined by Eq. 2 (SW D obs AN true w ). In practice, the data sources contain some error sources, and therefore, the solution of Eq. 2 can be calculated using Eq. 4 as below: By simplifying Eq. 10, we obtain (Brenot et al. 2018;Menke, 2012;Tarantola, 2005): Equation 11 defines how closely the estimated model parameters fit a true model. Therefore, according to Eq. 11, the model resolution matrix R m is given by Eq. 12 as follows (Aster et al. 2013;Menke 2012): According to Eq. 12, if the null space of the design matrix is trivial N ( A) {0} and, in other words, if the resolution matrix is identity (R m I), then, all voxels will be recovered exactly with high accuracy. R m can capture the experiment physics and geometry of the inversion problem because the resolution matrix is gained from the design matrix (Aster et al. 2013). Moreover, it reflects valuable information on the reconstruction quality of each inversion parameter (Menke 2012). Based on that, large diagonal elements with small off-diagonal elements indicate the desired parameter could be adequately resolved by the current geometry and quality of observation data. Large off-diagonal elements lead to low quality of the determined parameters.
Therefore, Menke (2012) applied the spread of the diagonal and off-diagonal elements of the resolution matrix to measure the goodness of the resolution in inverse problems using Frobenius norm as follows: where δ i j are the elements of the identity matrix I. Equation 13 is sometimes called the Dirichlet spread function. However, the Dirichlet spread function is not a proper measure of the goodness of resolution because the off-diagonal elements of this matrix are all weighted equally, despite whether they are close or far from the main diagonal (Menke 2012). Therefore, a weighting factor W i j can be added to solve this issue. Consequently, Eq. 13 is rewritten as follows (Menke 2012;Miller and Routh 2007;Toomey and Foulger 1989): Here, W i j is the physical distance between elements in [km]. The new spread function is frequently called the Backus-Gilbert spread (BG) function (Kaltenbacher et al. 2008;Piretzidis and Sideris 2016). In this research, W i j is calculated based on the Gaussian inverse distance (GID) between voxels in the same layer, namely horizontally, and therefore W i j is zero for voxels in different layers. From now on, we use the BGH abbreviation for the BG spread due to considering GID for every horizontal layer. According to Eq. 14, a well-resolved parameter corresponds to a smaller spread of the corresponding row of the model resolution matrix, whereas a poorly resolved parameter corresponds to a larger one.
Another way to define the spread is as follow (Maercklin 2004;Michelini and McEvilly 1991): where R m i is the L 2 norm for the ith row of the resolution matrix R m .This definition for the spread is known as Michelini spread function (Mich). In Eq. 15, D ji is the spatial distance between the ith and jth model parameters. Here, in contrast to BGH, the spread is negative, and a larger negative value corresponds to a well-resolved parameter.

Quality of tomography model
The accuracy of the reconstructed tomography model is typically evaluated using RMSE, Mean Bias, and standard deviation (Std) (Rohm and Bosy 2009;Shangguan et al. 2013;Xiaoying et al. 2014;Zhang et al. 2017). These statistics tools can be calculated by the following equations (Guerova 2003;Xiaoying et al. 2014): Here, n is the number of height levels and N wtomo and N wre f are the computed wet refractivity field from the tomography model and from the reference observation data like radiosonde profiles. In this study, we consider only voxels intersected by radiosonde profiles in different height levels as reference.

Case study
We performed a two steps validation approach to our spread framework: (1) firstly, we used a so-called 'synthetic dataset' [part of the collected dataset by COST Action ES1206 (https://www.cost.eu/actions/ES1206/)] composed of raytraced slant delays passing through the reference NWP fields (please refer to Sect. 3 of (Kačmařík et al. 2017) for more details), (2) and secondly real GNSS observations for the exact same area and a slightly different time span. The use of synthetic data could provide a reasonable judgment compared to the real dataset due to eliminating the main GNSS error sources in that. Thereby, two different datasets have been considered here to assess the values of spread with and without data noises in the COST benchmark dataset. The COST GNSS network is mostly placed in western parts of the Czech Republic and South Germany with 72 GNSS stations (see Fig. 1). Moreover, the tomography model ranges from 10.15°to 13.30°in longitude, 49°-51.70°in latitude and 0-15 km in WGS84 ellipsoidal height. More information about this dataset is detailed in (Douša et al. 2016). The full-time period of the synthetic dataset was chosen to cover European floods between 29 May and 14 Jun  in the year 2013. Especially the period DoYs 150-155 was characterised by heavy rainfalls. Moreover, a subperiod between DoYs 160-165 was analysed both with real and synthetic data. Therefore, three case studies are presented in the remaining part of this manuscript: (1) synthetic DoY 149-165, (2) synthetic DoY 160-165, and (3) real DoY 160-165. Furthermore, the tomographic model was designed with 50 km horizontal resolution and extended by 5°i n all directions to guarantee all rays pass completely through the model. According to previous researches, an exponential model was applied in the vertical direction (Manning 2013;Perler 2011). Furthermore, the ALADIN numerical weather model with 1 h time resolution was used as the initial field in this case study (Douša et al. 2016;Farda et al. 2010).

Numerical results
To investigate the efficiency of the spread as an indicator for the model accuracy, the correlation between Std, Bias, and • Loose constraints (LC): damping coefficient 0.1 (δ m 0.1) refer to Eqs. 6 and 12 • Tight constraints (TC): DAMPING coefficient 0.9 (δ m 0.9) refer to Eqs. 6 and 12 Therefore, the correlation between the values of spread and statistical parameters could be calculated according to the following diagram (see Fig. 2).
In the following, first, the accuracy of the reconstructed tomography profiles is investigated in both synthetic and GNSS datasets over the experimental periods. Next, the correlation between spreads (Mich and BGH) and statistical parameters (Std and Bias) are obtained to assess the spread as a quality indicator for GNSS tropospheric tomography.

Validation of tomography solution using NWM and RS data
The reference wet refractivity field for the GNSS dataset was calculated from the Kummersbruk (RS10771), and Meiningen (RS10548) radiosonde station gathered at hours 00:00  Figure 3 shows the agreement between the reconstructed wet refractivity profile, the RS10771 profile, and the NWM profile for two random days at midnight. According to this figure, the discrepancy between the synthetic tomography solution and the NWM profile is quite small. For the GNSS dataset, the inconsistency between the recovered tomography wet refractivity field and RS profile is visible, but the behaviour of both wet refractivity profiles is similar.
In Table 1, we give the average values of RMSE, Std, and Bias over all experimental periods. According to the obtained results, the retrieved tomography profiles using the synthetic dataset are generally underestimating the wet refractivity derived by the NWM model (ALADIN). In contrast, the reconstructed wet refractivity using the GNSS dataset is overestimated with respect to the wet refractivity obtained by the RS measurements. Moreover, due to the absence of observation errors, the quality of the reconstructed wet refractivity solution using the synthetic observations is better than for the GNSS observations.
The same procedure has been done for the location of the tomography model, intersected by RS10548. Figure 4 shows the retrieved wet refractivity profile performance compared to the reference profiles in two selected days on DoYs 162 and 164 at midnight. According to this figure, the agreement between the tomography solution using synthetic data and NWM data on DoY 162 is better than for DoY 164. However, in general, the estimated wet refractivity field using the synthetic dataset behaves like the NWM profile. For the GNSS dataset, the discrepancy between the tomography solution and the RS profile is smaller in the upper layers compared to the lower layers.
In addition, Table 2 gives the statistical results of comparing the tomography solutions and the reference profiles over the experimental period for both synthetic and GNSS datasets. According to the obtained results, the reconstructed wet refractivity profiles in all datasets and periods are underestimated with respect to the corresponding reference profiles. Moreover, similar to RS10771, the performance of the estimated tomography profiles in the synthetic dataset is better than for the GNSS dataset.

Correlation analysis of spread of resolution matrix
In order to analyse the correlation between spread and two statistical parameters, namely Std and Bias, first, the summation of the spread of voxels crossed by the radiosonde profiles was computed. Then, to gain a better interpretation of the dependency between statistical parameters and the spread, all values were normalised by the following formula, and then the differences between those values were considered.
whereby X and X n are the original value and the normalised value, respectively. Finally, Pearson correlation has been utilised to compute the similarity coefficient between statistical parameters (Std and Bias) and spread. Figure 5 shows the difference between spreads (BGH and Mich) and Std for three different time spans of the synthetic and GNSS dataset at RS10771. It can be seen from Fig. 5, that the difference between BGH spread and Std is smaller than for the Mich spread. Therefore, the time series of BGH spread follows more closely the Std variations compared to the Mich spread.
To better interpret Fig. 5, Table 3 gives the correlation between spread and Std time series for the different schemes and datasets. We must consider that small BGH spread and large Mich spread point to a well-resolved wet refractivity field due to their specification. Therefore, the computed correlation for BGH and Mich spreads are positive and negative, respectively. For the synthetic dataset, the negative correlation of the Mich spread in the long period is higher than for the short period. However, inspecting the results gained with the GNSS dataset, the correlation of the Mich spread is almost comparable to the synthetic dataset. According to the obtained results, the BGH spread shows the considerable positive correlation in all datasets. Altogether, both spread types show a considerable and promising correlation with the Std of the recovered wet refractivity for almost all investigated periods. However, applying LC results in a better match between spread and Std in comparison to TC. Consequently, the selection of C m significantly affects the obtained results.
In Table 4, the correlation between the Bias of the recovered refractivity field and the spread is given. According to these results, the correlation of BGH spread with Bias (Eq. 16) is highest in all datasets. However, as the Bias can be positive or negative, the correlation just reflects the tendency for a considerable absolute spread to a large bias. Moreover, the correlation of Mich spread shows reasonable large numbers between − 0.6 to − 0.4 for the GNSS dataset and synthetic dataset (149-165). In addition, compared to  TC, LC provides higher coherency with the spread based on the obtained correlation. The same analysis has been performed for radiosonde location RS10548. Figure 6 shows the differences between spread and the corresponding Std. According to that, the overall match between BGH spread and Std in the synthetic schemes is slightly better in comparison to Mich spread, which is consistent with the results obtained for radiosonde location RS10771 results. Table 5 presents the correlation between Std and the spread. Based on the obtained correlations, the general performance of the BGH spread is almost similar to Mich spread. Both spreads have a correlation of up to 0.8 with respect to Std for the GNSS dataset (Ref: RS wet refractivity) and synthetic dataset (Ref: NMW wet refractivity). TC schemes show similar performance in comparison to LC schemes. Unfortunately, the smallest correlation is obtained with the real GNSS observation dataset. The improved coherence of Mich spread with Std at this RS location is caused by the fact that the investigated model column pertains to boundary voxels of the model area.
As shown in Table 6, the correlation between refractivity biases and BGH spread is not considerable for the RS10548 location. Nevertheless, Mich spread shows a reasonable correlation with the Bias. Besides, the LC solution again performs better in comparison to TC.  6 Differences between spread and Std for RS10548. The left shows the solution by applying tight constraints on the a priori field, whereas the right column shows loose constraints In order to achieve a better interpretation of the obtained results, the absolute differences between Spread and Std of the studied time series could be compared. The smaller difference shows high similarity, whereas the larger difference indicates low similarity between time series. Therefore, the average of accumulated differences between spread and Std were computed. Figure 7 presents the comparison of the average accumulated differences during the periods of interest for the locations of RS10771 (hatched bars) and RS10548 (non-hatched bars). According to this figure, by applying loose constraints, spread provides a better characterisation of the accuracy of the tomography model in most cases. In addition, BGH spread has a higher consistency with the variations of Std in comparison to Mich spread. Consequently, BGH spread calculated with loose constraints on the a priori refractivity field are the recommended quantities to predict the accuracy of the retrieved tomography field.

Conclusions
Assessing the quality of the model reconstruction is one of the most important challenges in the GNSS tomography due to its application in storm now-casting, severe weather, and forecasting research. Previous tomography studies used statistical measures like RMSE and Std to assess the quality of the tomography profile in comparison to the external datasets like radiosonde or NWM. As a prerequisite, the whole tomography model has to be processed even though this procedure is considerably time-consuming. Consequently, defining a new approach to validate the tomography solution as well as predict a model accuracy would be essential, especially for increasing the applicability of GNSS tomography in meteorological studies. As a result, we investigated the performance of the spread of the resolution matrix to access if this parameter can be an indicator to predict the accuracy of the tomography model in this study. For this purpose, two different datasets, synthetic and GNSS, were considered to provide a reasonable judgment in the presence of various error sources in the measurements and without that. The case study was located in western parts of the Czech Republic and East Germany with 72 GNSS stations. Moreover, the simulation dataset covered the European floods between DoYs 149-165 (between 29 May and 14 Jun) in the year 2013 and DoYs 160-165 in the same year for the GNSS dataset.
In this study, two spread definitions, denoted Michelini (Mich) and Backus-Gilbert (BGH), were used in order to analyse the correlation of the model accuracy. Due to the impact of the observation covariance matrix as well as the quality of the initial field, the damped least-squares method was applied to calculate a reasonable resolution matrix. This is due to the key responsibility of the resolution matrix in spread computations. For the first implementation, we have only considered C obs as a diagonal matrix because calculating the full-populated measurements covariance matrix is quite challenging. However, this could be improved in the future by applying the turbulence theory to estimate off-diagonal elements of the covariance observation matrix. Besides, the a priori covariance matrix of the unknown parameters was defined by considering low and high damping coefficients, called LC and TC, respectively. Nevertheless, introducing the 'GNSS accuracy' of the wet refractivity model extracted from the NWM would lead to more realistic results. This can be achieved if the standard deviation of temperature, pressure, and water vapour pressure fields of the NWM is accessible.
To investigate the success percentage of spread as a proxy for the GNSS tomography, the correlation coefficient between these values and statistical measures, Std and Bias, were calculated over the experimental period for the synthetic and GNSS datasets. According to the obtained results, the correlation between spread and Std is considerable. However, the Bias shows different behaviour with respect to spread due to the unpredictable nature of this parameter. In addition, LC has a higher correlation in comparison to the TC solutions with respect to Std. Moreover, the absolute differences of BGH spread with respect to Std are generally smaller than Mich spread. Hence, we could conclude that applying BGH spread with LC weighting is the promising method to investigate the accuracy of a tomography model. Nevertheless, it should be noted that calculating a realistic prior covariance matrix of unknowns is necessary to achieve acceptable results.
In consequence, this work confirms the high correlation of the spread (up to 0.81 for synthetic and 0.72 for GNSS) with the Std of the retrieved refractivity field. Therefore, this parameter could be used in future as an a priori quality index for the tomography solution to analyse the performance of the tomography model before the reconstruction process, especially for the Near Real-Time (NRT) applications. In addition, this factor can be employed to recheck the reconstructed tomography solution to assure the quality of all parts of the model, which is essential for now-casting and forecasting applications. Nevertheless, further studies in different case studies and time periods are encouraged to assess the performance of spread.