Enhanced Troposphere Tomography: Integration of GNSS and Remote Sensing Data With Optimal Vertical Constraints

This article explores the enhancement of Global Navigation Satellite Systems (GNSS) tropospheric tomography by integrating remote sensing data and employing various vertical constraints. Wet refractivity modeling, critical for understanding atmospheric dynamics, has shown promising advancements. Leveraging tropospheric data from the Ocean and Land Color Instrument (OLCI), this research addresses the issue of empty voxels that impede GNSS-based tomography due to satellite and receiver geometries. Incorporating tropospheric data from remote sensing sensors mitigates empty voxels, enhancing retrieval accuracy for tropospheric water vapor. This study evaluates various vertical constraint functions in tropospheric tomography, presenting eight tomography schemes that utilize GNSS and OLCI data, highlighting their capacity to fill empty voxels without relying on empirical horizontal constraints. Results highlight the superiority of using OLCI observations in accuracy. Validation against radiosonde measurements and Weather Research and Forecasting model outputs affirms the reliability of this approach. Integrating OLCI observations with GNSS data reduces the average root mean square error by approximately 27%, with the Gaussian function exhibiting superior vertical constraint performance.


I. INTRODUCTION
W ATER vapor is regarded as a crucial component of the troposphere and is the most important parameter for wet refractivity, and its role in weather cycles and meteorological phenomena has been emphasized [1].However, the irregular spatial and temporal variations of water vapor and wet refractivity have posed challenges in its modeling and prediction [2], [3].To gain a comprehensive understanding and reconstruct the behavior of wet refractivity, obtaining a 3-D distribution becomes imperative, which can be accomplished through measurement methods like radiosonde and water vapor radiometer.
Nevertheless, these methods suffer from limitations in spatial and temporal resolution, as well as weather dependency [4].Despite these challenges, recent advancements have showcased the remarkable capabilities of Global Navigation Satellite Systems (GNSS) tropospheric products in presenting crucial atmospheric information.The ability of computed tropospheric delay from GNSS measurements has been proven in different fields of study such as machine learning-based subsidence prediction [5], drought monitoring [6], weather prediction [7], [8], and GNSS troposphere tomography [9].At present, GNSS tropospheric tomography stands out as a potent technique for obtaining 3-D atmospheric wet refractivity distribution based on the tropospheric delay of GNSS signals.The pioneering troposphere tomography work by Flores et al. [9] employed voxel-based tomography to estimate wet refractivity in empty voxels by applying horizontal and vertical constraints.Subsequent research has made significant progress in enhancing the accuracy of tomography by utilizing techniques, such as the damped least squares method [10], simulated data and atmospheric models [11], Gaussian weight function-based optimization of horizontal constraints [12], and Moore-Penrose pseudoinverse in local tomography over mountainous regions [13].In addition, researchers have explored the optimization of the GNSS tomography model, constructing correlations through robust Kalman filters [14], employing weighted equations for improved accuracy [15], developing novel algorithms considering GNSS signals from the side face of the tomographic area [16], [17], and optimizing the distribution of voxels horizontally and vertically to reduce the number of empty voxels [18].Further advancements include the reconstruction of precise ray paths through 3-D ray tracing methods [19], the adoption of various regularization techniques for inverse tropospheric tomography [20], and the introduction of function-based troposphere tomography, accompanied by a comparative analysis with voxel-based tomography, primarily focusing on positioning [21], [22], [23], [24].An alternative approach known as node-based tomography has also been explored, involving the estimation of tropospheric states at specific nodes [1], [25], [26].Moreover, the tomography results have found application in data analysis within the Weather Research and Forecasting (WRF) global prediction model under varying weather conditions [27].Dynamic tomography models have been employed [28], and WRF model data have been utilized to merge empty voxels, thereby reducing the number of unknown parameters [29].Furthermore, researchers have examined the impact of external observations on the tomography problem, employing Global Positioning System (GPS) observations and additional positioning systems like Galileo, GLONASS, and BeiDou, with outcomes varying [30], [31].They have also integrated data from radio occultation [32], [33], [34], interferometric synthetic aperture radar [35], [36], and atmospheric infrared sounder [2] into the tomography model and used GNSS reflectometry signals to address rank deficiency issues [37].In addition, some studies suggested the use of water vapor and wet refractivity products obtained from the Moderate-Resolution Imaging Spectroradiometer (MODIS) sensor, which has a higher spatial resolution [38], [39], [40].
Empty voxels present a significant challenge in voxel-based tomography, primarily due to the characteristics of GNSS satellite and receiver geometries.They lead to rank deficiency in the coefficient matrix of the tomography problem.This challenge becomes more pronounced in lower layers and near the Earth's surface, where numerous voxels lack ray passage, necessitating the imposition of additional empirical horizontal constraints on the problem.Empirical constraints may introduce errors in the final results.To overcome this challenge, we integrate wet refractivity data from the Ocean and Land Color Instrument (OLCI) as external observations in the tomography process, effectively covering empty voxels without relying on horizontal constraints.By utilizing integrated water vapor (IWV) as a dependable indicator, derived from a range of remote sensing sensors, it becomes possible to accurately represent the total water vapor content within the vertical atmospheric column.The conversion of IWV to slant wet delay (SWD) or slant water vapor (SWV) not only increases the degrees of freedom but also reduces the number of empty voxels, resulting in a favorable cone-shaped observation geometry that complements the GNSS setup.This combination of GNSS and remote sensing data is expected to enhance the retrieval of tropospheric water vapor.Another critical aspect of tomography involves the use of vertical constraints based on a function to create vertical separation in the tomography solution.In addition to evaluating the impact of using OLCI measurements to increase the number of signals and reduce empty voxels, this study investigates the effects of different equations as vertical constraints in troposphere tomography when using different types of input data.Finally, the obtained results are validated using radiosonde measurements and WRF model outputs.The following sections will delve into the theoretical background, the study area, the introduction of the dataset, and a discussion of the obtained results.

A. Principal of Troposphere Tomography
The troposphere introduces delays and deflections to GNSS signals as they traverse it.The slant total delay encompasses both hydrostatic and wet components [41], [42].The slant hydrostatic delay signifies the hydrostatic portion's delay, whereas the SWD quantifies the delay attributed to the wet part of the tropospheric layer on GNSS signals [41].SWD serves as input for tomography, facilitating the reconstruction of wet refractivity.The fundamental tomography equation in this context is expressed as follows [9]: Rec.

N w ds
(1) where s denotes the ray's length and N w represents the wet refractivity value.Although (1) is applicable to continuous spatial problems, voxel-based tomography involves a discrete space divided into 3-D voxels.The corresponding equations for tomographic observations are also discretized as follows [9]: where i serves as the GNSS ray counter, j represents the number of voxels, Δs i j denotes the length traveled by ray i within voxel j, and N wj stands for the wet refractivity within voxel j.The matrix representation of this equation is provided as follows: where A GNSS denotes the coefficient matrix based on the GNSS data and N w represents the vector of unknowns.Given that tomography (3) is ill-posed and further compounded by the deficiency of rank in the coefficient matrix, to achieve vertical separation in the tomography solution, the application of horizontal and vertical constraints becomes necessary to estimate the unknowns in the tomography problem [9], [29], [30].After applying constraint equations, (3) will be transformed into the following forms: where A H and A V are horizontal and vertical constraints, respectively.Generally, tropospheric tomography is a large and ill-conditioned inverse problem due to the high number of observations and a wide area of modeling.Therefore, the use of regularization methods is necessary.In this article, the simultaneous algebraic reconstruction technique (SART) is used [43].SART belongs to the broader category of reconstruction techniques known as the simultaneous iterative reconstruction technique (SIRT).Both SART and SIRT methods are utilized in medical imaging and other fields to reconstruct images from projection data, such as those obtained in tomography.SIRT constitutes a general class of algorithms encompassing various iterative reconstruction techniques.SART, on the other hand, is a specific iterative algorithm within this class [43].While SART primarily focuses on simultaneous updates of the image and the data residual, other SIRT methods may employ different update strategies and regularization techniques.In summary, SART represents a particular instance of an iterative reconstruction method, falling under the broader umbrella of SIRT techniques.

B. Constraints
As previously mentioned, certain voxels within the tomography model are not intersected by any rays due to the geometric constraints imposed by the positioning of ground receivers and Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
satellites.It produces a rank-deficient coefficient matrix in tomography problems.One common approach to address this issue is the use of horizontal constraints.This technique only resolves the lack of structure in the coefficient matrix.Various constraint methods have been employed in tomography problems, and some of them are provided here.Equations incorporating constraints have been utilized [10], alongside techniques such as the Kalman filter and algebraic reconstruction to eliminate the null space within the coefficient matrix [9], [44].In addition, Gaussian weight functions have been employed as horizontal constraints [12], and wind current analysis has been introduced as an extra constraint [45].Another strategy involves combining adjacent voxels with similar water vapor values using WRF data to reduce the number of unknowns and unoccupied voxels [29].The horizontal constraint is utilized to estimate the water vapor content in vacant voxels, such that the content is the weighted average of neighboring voxels within the same height layer.This relationship can be expressed as follows [46]: where m denotes the number of voxels within each height layer, x m is the wet refractivity in the mth voxel, and w j represents the weighted horizontal coefficients of the jth voxel.These coefficients can be computed using the following formula: where d j,i signifies the distance between the center of the unoccupied voxel and the centers of the remaining voxels within the layer.This particular horizontal constraint is also employed in this study.
On the other hand, the vertical constraint serves to rationalize the problem and establish the overall water vapor pattern across different height layers.The inclusion of radiosonde and radio occultation data has been instrumental in achieving a singular solution, as proposed in previous studies [33].However, in most prior research, the well-known negative exponential function has been adopted as the vertical constraint, formulated as follows [47]: where h j , h k , N w j , and N w k refer to the height and wet refractivity in layers j and k, respectively, and H is the water vapor scale height, which is empirically 1-2 km.The vertical constraint is utilized to ascertain the general water vapor trend within the study region, considering the spatial and temporal irregularities present in tropospheric water vapor.The suitable function as vertical constraints depends on the geographical and weather conditions of the study area.Therefore, in this study, various functions have been explored for reconstructing the water vapor trend within the study area, ultimately identifying the most suitable function for use as a vertical constraint.The parameter n denotes the number of terms.In this article, we determine the value of n based on the best performance with varying numbers of terms for each function, as assessed by their fitness using meteorological data.The selection of these functions is grounded in their performance in local modeling.

C. GNSS-Based SWD Estimation
Within GNSS measurements, the computation of tropospheric delay occurs in the zenith direction, referred to as the zenith total delay (ZTD).This delay comprises two key components: the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD).Determination of ZHD is accomplished through the utilization of the Saastamoinen model [52] as follows: where P s represents the surface pressure in hPa, ϕis the latitude, and H s denotes the height of the station in meters.Subsequently, ZWD can be derived as the difference between ZTD and ZHD.The estimation of ZTD is achievable by processing GNSS signals via specialized software for GNSS signal processing.Alongside ZTD, the outcomes of GNSS processing encompass horizontal gradients.These gradients indicate the nonuniform distribution of atmospheric gases.Consequently, the specific SWD can be reconstructed using the following Equation [9]: where G W NS and G W EW represent the wet horizontal gradients in the north-south and east-west directions, respectively.Moreover, Az and α correspond to the azimuth and elevation angle of the satellite, respectively, R denotes the cleaning postfit residuals, mf wet characterizes the wet mapping function, and mf g represents the horizontal gradient function.

D. Unveiling SWD Computation Through Remote Sensing Data
The OLCI instrument functions as a multispectral imaging spectrometer, installed on the Copernicus Sentinel-3 satellites, namely, Sentinel-3A and 3B.This advanced sensor has been developed using the opto-mechanical and imaging framework originally utilized in the Medium Resolution Imaging Spectrometer sensor, which was part of the ENVISAT satellite [53].This particular instrument provides three distinct levels of data accessibility to the public audience.For the purposes of this research, we exclusively utilize data from the Sentinel-3 OLCI Level 2 Land Full Resolution (OL_2_LFR) products, which include invaluable IWV observations.Within the capabilities of this sensor, the production of the OLCI-IWV product is achieved through the water vapor absorption channel with a spatial resolution of 300 m [54].The retrieval of water vapor from OLCI measurements relies on the differential absorption technique, which involves comparing the measured radiance at a nonabsorbing spectral band with the reference water vapor absorption band.While the measured radiance in each spectral band depends on factors such as solar irradiance, atmospheric transmittance, and surface albedo, the differential absorption technique allows us to derive IWV as [55] where L 18 and L 19 denote the measured spectral radiance, whereas k 19 represents the mass extinction coefficient of water vapor.
Due to the coexistence of multiple observations within a voxel column, it becomes imperative to amalgamate these diverse observations into a singular entity.Consequently, a solitary observation equation emerges, signifying that the count of observation equations resulting from OLCI-IWV equals the voxel count within a given layer.It is important to underline the distinction between GNSS observations, which align with the satellite's line of sight, and OLCI-IWV values that exclusively pertain to the zenith direction.The latter solely provides vertical distances within the voxel column where the OLCI-IWV pixel is positioned.Thus, for the utilization of OLCI-IWV values, a preliminary transformation is essential to convert them into an oblique form, directed between the OLCI sensor and the corresponding pixel.This process can be likened to the methodology employed in estimating SWD from GNSS observations.Fig. 2 depicts the general framework of advanced troposphere tomography in the presence of remote sensing data, along with a 3-D schematic scheme for estimating SWD through its utilization.The reconstruction of oblique rays between the OLCI sensor and the pixels hinges on the determination of elevation angle and azimuth, ultimately enabling the estimation of water vapor content along these rays.It is noteworthy that the conversion of OLCI-IWV values occurs in two sequential steps: first to ZWD and subsequently from ZWD to SWD.The subsequent equations elucidate the mechanisms of these transformations as [39], [56] with where Γ w (α, Az) refers to the wet delay caused by the horizontal gradient, R v is the gas constant for water vapor, k 2 and k 3 are the experimental constants, and T m is the average atmospheric temperature, which is calculated using the local surface temperature T s and empirical formula T m = 70.2+ 0.72 T s [41], [57].
It is worth mentioning that when using remote sensing observations, no voxels remain empty of signals, and the coefficient matrix does not have a rank deficiency.Therefore, empirical horizontal constraints can be eliminated from the problem.After applying SWD from remote sensing observations in the tomography problem, the matrix form of the problem will be as follows: where A RemoteSensing is the coefficient matrix obtained from remote sensing observations.

III. STUDY AREA AND DATASET
To explore the impact of incorporating remote sensing information, a distinct region within North America was selected for analysis.To comprehensively assess the effectiveness of the proposed approach, a collection of observations from 21 GNSS stations was compiled over a 7-day period, spanning Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.multiple months in the year 2021 and encompassing a wide range of weather conditions.This particular time frame was chosen deliberately due to the dynamic variations in humidity indices observed at the core of the study area.Elaborate representations detailing the spatial distribution of GNSS stations and the topographical attributes of the study area can be found in Fig. 3.
The Sentinel-3A mission was launched on February 16, 2016, followed by the launch of Sentinel-3B on April 25, 2018.Initially, Sentinel-3B was operated in tandem with Sentinel-3A until it reached its nominal orbit in December 2018.The repeat cycle of the Sentinel-3A/-3B satellites is 27 days.However, global coverage in 2/3 days is enabled by OLCI due to its large swath of 1270 km.In this article, the better satellite for obtaining IWV data from the OLCI mission has been selected based on the synchronization of tomographic epochs and the visiting time of Sentinel-3A and Sentinel-3B.A set of OLCI data from track numbers 10, 11, and 12 for the same epochs of tomography processing has been utilized to extract IWV information.
Previous studies show the ability of the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data including ERA-Interim and ERA5 data in different fields of geodesy and remote sensing such as tropospheric delay computation using different kinds of ray tracing techniques and reducing the tropospheric effect in radar observations [58], [59], [60].In this study, the ERA5 data from ECMWF have been used to extract meteorological data needed for computing ZHD and to estimate the coefficient of vertical constraints [61].
The accuracy and credibility of the obtained outcomes were fortified by exclusively incorporating data retrieved from a radiosonde station situated within the study locale.Leveraging radiosonde technology facilitated the vertical profiling of tropospheric parameters.To thoroughly assess the tomographic results across distinct sections of the tomography model, the WRF model outputs were used.Fueling the WRF model and facilitating the prediction of meteorological conditions involved utilizing inputs drawn from the Global Forecast System (GFS) analysis data.Administered by the National Centers for Environmental Prediction, the GFS operates as a numerical weather forecast model that provides a comprehensive array of atmospheric and land-soil variables.Distinguished for its advanced data assimilation system and adept parallel processing capabilities, the WRF model stands as a widely embraced tool for conducting mesoscale numerical weather prediction.Its utility extends to both atmospheric research and the realm of operational weather forecasting.The model's versatility enables simulations spanning a spectrum of meteorological scenarios, effectively accommodating resolutions ranging from the sub-tens-of-meter scale to expansive scales reaching thousands of kilometers [62].

A. Analysis of GNSS Data
The processing of GNSS data, which encompassed observations from both GPS and GLONASS observations, was carried Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.out using Bernese software version 5.2, known for its recognized reliability and adaptability [63].This software enables data processing in both static and kinematic modes, supporting zero-difference and double-difference observables through relative and precise point positioning (PPP) approaches.In this study, the PPP technique was employed.In our data processing approach, consideration was given to the ionospheric-free linear combination, and ZTD was calculated at 15-min intervals, incorporating gradients in both north-south and east-west directions, all with a 1-h resolution.To convert the zenith direction into a slant direction, the global mapping function was applied [64].Processing parameters to calculate the tropospheric delay as input to the tomography methods have been considered based on Table II.
Examples of the ZTD values obtained through this process are illustrated in Fig. 4. Variations in delay observed at these stations are attributed to diverse topographical features and changing  weather conditions.ZWD was computed by subtracting the ZHD calculated using the Saastamoinen model from the ZTD.The Saastamoinen model was applied using ERA5 data.In the subsequent phase, the ZWD was projected along the line of sight of the satellites.Subsequently, the SWV derived from the SWD was considered as input for the tomography technique.It should be noted that in the context of SWD estimation, postfit residuals were not incorporated.In some of the previous studies, researchers tried to improve the accuracy of SWD values by adding postfit residuals of the GNSS analysis to include the unmodelled parts of the delays [68], [69], [70].However, recent studies show that the postfit residuals do not have a significant impact on increasing the accuracy of SWD [71].Some studies Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.even mention that errors of the GNSS measurements, e.g., multipath, could remain in the residuals [26], [72].Therefore, in this study, the postfit residuals have not been considered.

B. Analysis of Remote Sensing Data
The IWV data have found within the Sentinel-3 OL_2_LFR product.The term "full resolution" is indicative of a pixel size of 300 m in the OLCI IWV imagery.To acquire the OLCI imagery and determine the IWV values, the Copernicus Open Access Application Programming Interface (API) hub was utilized.Initially, the Sentinelsat module was employed to conduct a spatial query through the Copernicus API for each relevant satellite pass.The goal was to identify Sentinel products that met the specified criteria, including mission, product type, and the study area's boundary.The identified Sentinel-3 products were subsequently automatically downloaded and organized into separate folders.Within each Sentinel-3 product folder, an examination was carried out on the IWV value, its associated error, geographical coordinates, date, and time stamps.Fig. 5 displays the extracted IWV maps from OLCI data in tomographic epochs.
After extracting IWV information, the ZWD and SWD data have been reconstructed using the formula mentioned in Section II.Fig. 6 shows the average root mean square error (RMSE) of ZWD obtained from OLCI data and GNSS measurements at the location of the GNSS stations.As can be seen, the range of RMSE is between 4 and 20 mm, proving that the obtained ZWD from OLCI is confident and reasonable.
It should be noted that the necessary gradient for computing SWD from IWV data has been extracted from the GRAD dataset [73] based on the standard gradient model by Chen and Herring [65].

C. Tomography Model
The skeletal structure employed for tomographic processing in this research is visible in Fig. 7.In alignment with previous research findings, the optimal horizontal resolution for the voxel-based tomography model has been ascertained through an examination of the resolution matrix.This matrix, a pivotal component of the coefficient matrix, embodies both the resolution and geometry of the tomography model.It should be emphasized that the insignificance of any diagonal elements within this matrix has an adverse impact on the accuracy of parameter estimation.The quest for an optimal design of the tomography model has yielded a resolution matrix that closely approximates the identity matrix, a phenomenon corroborated by some of the previous studies [19], [22], [30].
In this particular context, the horizontal resolution has been fixed at 0.2°for the voxel-based tomography model.Regarding vertical resolution, in contrast to prior studies where various methods, such as the exponential model, were utilized to partition vertical layers [1], our experimentation in this study area has unveiled that maintaining a constant resolution in the lower segment of the tomography model yields superior outcomes.As a result, a high yet consistent resolution has been adopted up to a specified height above ground level.Specifically, the vertical resolution is established at 500 m for the initial six layers and subsequently transitions to 1000 m for the ensuing seven layers.Considering these predefined horizontal and vertical resolutions, the model encompasses a total of 325 voxels and unknowns.These resolutions, in conjunction with the topography data, constitute the foundational framework of our tomography study.

D. Results
Using the aforementioned methods, tropospheric tomography results have been obtained in eight different schemes based on GNSS, OLCI, horizontal, and vertical constraints.The defined schemes are visible in Table III.It is worth mentioning that when remote sensing observations are used alongside GNSS observations, all of the empty voxels are filled with signals.However, in cases where remote sensing observations are not utilized, nearly 20% of the voxels remain empty, necessitating the use of horizontal constraints to compensate for the rank deficiency of the coefficients matrix and to solve the tomography problem.Regarding the vertical constraints, it should be noted that ERA5 data were used to estimate the coefficients of different functions.It was not possible to use radiosonde measurements and WRF model outputs for this purpose because these data are reserved for validating the final results and must not be used in the tomography processing.This ensures that the validation data and obtained results remain independent.The processes have been carried out in seven epochs during different months of the year 2021.An example of reconstructed wet refractivity using the GO_gau scheme can be seen in Fig. 8.
Based on Fig. 8, it is evident that variations in wet refractivity during 7 months of the year, under varying weather conditions, are observable.In April and May, the lowest amount of wet refractivity is calculated.In contrast, during June, July, and August, the wet refractivity levels are higher, especially in August.This phenomenon can be attributed to the generally higher average humidity during these months.Furthermore, the trend of wet refractivity variation along the vertical direction exhibits a decreasing trend, which is normal since wet refractivity decreases with increasing altitude.It should be noted that all the schemes used successfully reconstructed 3-D wet refractivity with a decreasing trend along the vertical direction and higher values during more humid months.

E. Validation and Discussion
To check the effect of using remote sensing data and different functions as vertical constraints, it is necessary to validate the results of different schemes.First, the results have been validated using observations from only radiosonde stations in the study area.Fig. 9 shows the vertical profile obtained from tomography schemes and radiosonde measurements at the location of the radiosonde station.Table IV displays the average statistical parameters of this comparison, which include the RMSE, Nash-Sutcliffe efficiency (NSE), R 2 , and mean absolute error (MAE) [74].Among the various tomography schemes evaluated in this study, the one with the lowest RMSE is the GO_gau scheme, which exhibited an RMSE value of 1.0585 ppm.A lower RMSE value indicates a better agreement between the model-generated profiles and the actual measurements.The MAE values for the different schemes, ranging from 0.9289 ppm (for GO_gau) to 1.2859 ppm (for GH_pow), consistently indicated that all schemes achieved a high level of accuracy in their estimates when compared to the radiosonde measurements.This suggests that the tomographic schemes, including GO_gau, effectively captured the vertical distribution of tropospheric properties, with minimal deviations from the observed data.NSE, a measure of how well the model replicates the variability in the observed data, demonstrated high values across the board.The NSE values for the various tomography schemes ranged from 0.9797 (for GH_pow) to 0.9892 (for GO_gau).These high NSE values indicated that the schemes skillfully reproduced the temporal and spatial variability in the tropospheric profiles, and GO_gau remained at the forefront in terms of performance.The R 2 further reinforced the findings, with R 2 values consistently exceeding 0.84 for all schemes.Specifically, GO_gau achieved an R 2 value of 0.8489, underlining its strong capability to explain the variability in the radiosonde measurements.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.There is only one radiosonde station in the study area, making it impossible to validate the results in other parts of the tomography model.Therefore, the obtained results have been evaluated using WRF model outputs in addition to radiosonde measurements.Table V lists the physics schemes used in the WRF model.Fig. 10 shows an example of the difference between the reconstructed 3-D wet refractivity using different schemes and the WRF model outputs.In addition, Fig. 11 illustrates the comparison between WRF data and the wet refractivity obtained from the GH_gau and GO_gau schemes.To facilitate a comprehensive comparison of different tomography schemes and WRF data, average statistical parameters are presented in Table VI.
The GO_gau method demonstrates exceptional performance in terms of the lowest RMSE (2.2710 ppm), signifying its superior accuracy in tropospheric profiling.This method effectively combines GNSS and OLCI observations, leveraging the complementary information from both sources.In addition,     accurate atmospheric profiling.GH methods display slightly lower NSE values but still maintain acceptable accuracy levels.R 2 measures the proportion of variance in observed data that is explained by the model.Higher R 2 values suggest a better model fit.In this aspect, GO methods consistently perform exceptionally well.GO_gau stands out with an impressive R 2 of 0.9645, indicating its excellent ability to explain the variance in atmospheric profiles.GO_exp also exhibits a high R 2 value, emphasizing its suitability for accurate tropospheric profiling.While GH methods exhibit slightly lower R 2 values, they still provide meaningful explanations of variance in tropospheric profiles.In this analysis, OLCI consistently outperforms horizontal constraints in both RMSE and MAE across all vertical constraint functions.This suggests that OLCI provides more accurate results in predicting atmospheric parameters.Moving on to NSE and R 2 , which measure goodness of fit and model efficiency, OLCI again demonstrates superior performance compared to horizontal constraints.This implies that OLCI not only offers accuracy but also excels in representing observed data, making it a more reliable choice for tomography modeling.Furthermore, the choice of vertical constraint functions significantly impacts the accuracy and reliability of both OLCI and horizontal constraint methods.Examining the performance of these functions within the OLCI method.The Gaussian vertical constraint function produces commendable results, with relatively low RMSE and MAE values, indicating high accuracy.It also yields competitive NSE and R 2 values, showcasing its suitability for modeling vertical constraints within the OLCI method.The exponential function, while reasonable, generally falls slightly behind the Gaussian function in terms of RMSE, MAE, NSE, and R 2 .However, it can still be a viable choice depending on specific modeling requirements.The polynomial function performs well in the OLCI method, particularly in terms of NSE and R 2 .It demonstrates a strong capacity to fit observed data.The Power function, while adequate, tends to have relatively higher RMSE and MAE values compared to Gaussian and polynomial functions within the OLCI method.However, it still provides reasonably good NSE and R 2 values, suggesting its effectiveness in certain scenarios.The statistical analysis of tomography results and WRF model outputs indicates that the OLCI method, when combined with appropriate vertical constraint functions, outperforms the horizontal constraints approach across all key metrics: RMSE, MAE, NSE, and R 2 .The choice of vertical constraint function within the OLCI method should align with specific modeling objectives, but Gaussian and polynomial functions emerge as strong contenders for their overall better performance.Ultimately, the selection between OLCI and horizontal constraints, as well as the choice of vertical constraint function, should be made based on the specific goals and constraints of the tomography modeling project, while considering the tradeoffs between accuracy and computational complexity.

Ⅴ. CONCLUSION
In this article, remote sensing data and various vertical constraints were employed to enhance GNSS tropospheric tomography.Wet refractivity modeling, critical for understanding atmospheric dynamics, was an area where promising advancements were made.Leveraging tropospheric data from the OLCI, this research addressed the issue of empty voxels that were impeding GNSS-based tomography due to satellite and receiver geometries.Incorporating tropospheric data from remote sensing sensors was a key element in mitigating empty voxels, enhancing retrieval accuracy for tropospheric water vapor.The study evaluated various vertical constraint functions in tropospheric tomography, presenting eight tomography schemes that utilized GNSS and OLCI data, highlighting their capacity to fill empty voxels without relying on empirical horizontal constraints.The results highlighted the superiority of using OLCI observations in accuracy.Validation against radiosonde measurements and WRF model outputs affirmed the reliability of this approach.Integrating OLCI observations with GNSS data reduced the average RMSE by approximately 27%, with the Gaussian function exhibiting superior vertical constraint performance.

Enhanced
Troposphere Tomography: Integration of GNSS and Remote Sensing Data With Optimal Vertical Constraints Saeed Izanlou , Saeid Haji-Aghajany , and Yazdan Amerian Abstract-This article explores the enhancement of Global Navigation Satellite Systems (GNSS) tropospheric tomography by integrating remote sensing data and employing various vertical constraints.Wet refractivity modeling, critical for understanding atmospheric dynamics, has shown promising advancements.Leveraging tropospheric data from the Ocean and Land Color Instrument (OLCI), this research addresses the issue of empty voxels that impede GNSS-based tomography due to satellite and receiver geometries.Incorporating tropospheric data from remote sensing sensors mitigates empty voxels, enhancing retrieval accuracy for tropospheric water vapor.This study evaluates various vertical constraint functions in tropospheric tomography, presenting eight tomography schemes that utilize GNSS and OLCI data, highlighting their capacity to fill empty voxels without relying on empirical horizontal constraints.Results highlight the superiority of using OLCI observations in accuracy.Validation against radiosonde measurements and Weather Research and Forecasting model outputs affirms the reliability of this approach.Integrating OLCI observations with GNSS data reduces the average root mean square error by approximately 27%, with the Gaussian function exhibiting superior vertical constraint performance.Index Terms-Global navigation satellite systems (GNSS), Ocean and Land Color Instrument (OLCI), troposphere tomography, vertical constraints.

Fig. 1
presents illustrative examples showcasing the diverse characteristics exhibited by these functions based on the meteorological data from one epoch.According to Fig. 1, power and Gaussian functions exhibit more variations from the surface to a height of 10 km.

Fig. 2 .
Fig. 2. Geometry of GNSS and remote sensing observations (left) and schematic of SWD estimation using remote sensing observations (right).Here, 1 and n refer to signal numbers.

Fig. 3 .
Fig. 3. Topography of the study area and 3-D distribution of GNSS receiver (circles) and radiosonde station (square).

Fig. 4 .
Fig. 4. Computed ZTD for 21 GNSS stations on one of the processing days.

Fig. 9 .
Fig. 9. Comparison between two schemes of tomography results and radiosonde measurements.
capture the same level of detail as the GO methods, which make use of both GNSS and OLCI data.Moving on to the MAE, it complements RMSE by measuring the average magnitude of errors.Once again, GO methods demonstrate better performance compared to GH methods, with lower MAE values.Notably, GO_gau excels with the lowest MAE (1.8143 ppm), reinforcing its precision in tropospheric profiling.GO_pol and GO_exp maintain competitive MAE values, indicating their ability to provide accurate atmospheric profiles.Although GH methods exhibit slightly higher MAE values, they still offer respectable atmospheric profiling accuracy.NSE evaluates the efficiency of the tomographic methods in replicating observed data.Values close to 1 indicate a close match between the model and observations.Here, GO methods outperform GH methods, with GO_gau leading with an NSE of 0.9205, showcasing its excellent capability to replicate observed data.GO_exp, GO_pol, and GO_pow also achieve high NSE values, underlining their suitability for Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 10 .
Fig. 10.Difference between results obtained from GO_gau scheme and WRF model outputs.

TABLE I DIFFERENT
FUNCTIONS INVESTIGATED IN THIS STUDY FOR CONSTRUCTING THE VERTICAL CONSTRAINT Fig. 1.Illustrative examples of different functional behaviors.

TABLE III DIFFERENT
SCHEMES USED IN THIS STUDY

TABLE V PHYSICS
SCHEMES USED IN THE WRF MODEL CONFIGURATIONGO_pol and GO_exp also exhibit competitive RMSE values, with 2.4256 ppm and 2.47012 ppm, respectively.These results underscore the effectiveness of integrating GNSS and OLCI data in enhancing tropospheric profile accuracy.Conversely,

TABLE VI AVERAGE
STATISTICAL PARAMETERS OF COMPARISON BETWEEN TOMOGRAPHY RESULTS AND WRF MODEL OUTPUTS GH methods show slightly higher RMSE values, with GH_gau having the lowest RMSE (2.7301 ppm) among them.While GH methods incorporate horizontal constraints, they may not