GLOBAL POSITIONING SYSTEM PRECIPITABLE WATER VAPOR INTERPOLATION USING INVERSE MULTIQUADRIC, ARTIFICIAL NEURAL NETWORK AND INVERSE DISTANCE WEIGHTED

: Precipitable water vapor (PWV) is one of the most critical data in many meteorological departments. This component has great spatial and temporal changes, so the global positioning system (GPS) always seeks to increase the accuracy of estimating the water vapor component in the troposphere. The waves sent by the satellites of this system are delayed due to passing through atmospheric layers such as the troposphere. In this paper, interpolation methods are used to estimate precipitable water vapor. Inverse multiquadric (IMQ) interpolation which is based on radial basis functions, artificial neural network (ANN) method, and inverse distance weighted (IDW) which are the most common method of interpolation in meteorology. A region in North America with 23 GPS stations was randomly selected. Then, the interpolation of precipitable water vapor on a summer day is done using GPS data. The root mean square error value (RMSE) for the IMQ method was the lowest compared to other methods and was equal to 2.11 mm. Finally, using the IMQ interpolation method, a dense map of Precipitable water vapor changes in the troposphere layer is developed for the study area.


INTRODUCTION
With the advent of GNSS and advancements in geoscience studies and weather forecasting, further analysis of the atmosphere and its parameters has become particularly important.Furthermore, it has caused the emergence of different methods to achieve higher accuracy for estimating atmospheric parameters (Yang et al., 2019).The troposphere layer is essential in atmospheric studies, for the reason that GNSS signals reach the ground in a bent and delayed manner due to refraction caused by the troposphere layer (Nilsson et al., 2013).The tropospheric delay of the GNSS signal consists of dry and wet parts.The prediction of the zenith wet delay and the interpolation of its parameters is difficult due to the rapid temporal and spatial change of water vapor compared to the dry part (Smith and Weintraub, 1953).besides, the parameter of PWV is a function of ZWD and has many uses in weather forecasting programs (Chen et al., 2018).using interpolation, the desired parameters in the troposphere can be obtained.In general, interpolation is a process by which the value of a quantity at points with known coordinates can be obtained using the value of the same quantity at other points (Davis, 1975).Therefore, to get the amount of PWV that can be precipitated in the troposphere layer, we can use the interpolation methods such as inverse multiquadric, artificial neural network, and inverse distance weighted.To do that it is necessary to have the PWV value in the desired reference stations.Therefore, first, the tropospheric delay in the zenith direction should be calculated using the precise point positioning method (PPP) from the GPS observation.To calculate ZHD, pressure observations are needed, which can be obtained from equation (1) by using the dry experimental model of saastamoinen, which is one of the most accurate global prediction models (Mendes, 1999, Saastamoinen, 1972).0.002277 ZHD 1 0.00266 cos(2 ) 0.00000028 where P is surface pressure ( hPa ),  is latitude and, h is altitude (m).The value of the ZWD is calculated using equation (2) (Bevis et al., 1992): Then, PWV is obtained from the following equation : is the atmospheric weighted average temperature.Using the Bevis formula do that as follows (Bevis et al., 1992): where T s is the surface temperature is in Kelvin at the location.

Radial basis functions
The radial basis function (RBF) is a universal interpolation method (Leyla, 2011).This method is defined as the weighted sum of transformations of the radially symmetric basis function.
In this method, a set of distinct point data and corresponding observation values are considered, and the primary radial basis function used for interpolation is expressed as follows (Sharifi, 2016): where n is the number of reference stations and U (r) i is a function of the kernel interpolator.This function depends on the Euclidean distance.r is the Euclidean distance and w i are the unknown expansion coefficients and can be calculated with the input values of the kernel interpolator through the linear system.The generalized version of the radial basis function of Hardy, which was developed for the preparation of topographic maps, is in the form of the following equation (Leyla, 2011)  are unknown polynomial coefficients The second term of the equation (6) represents the polynomial function whose coefficients are as follows: As a result of the linear interpolation system, the number of unknowns is more than the number of observations, and conditions are added for its solution.Therefore, m u more equations are needed to form a certain linear system.Consequently, equation (8) must hold (Leyla, 2011, Sharifi, 2016).
In other words, in this article, according to the data in the 3 dimensions of longitude, latitude and altitude, the relationships The combination of interpolation conditions and additional equations forms a linear system of equations, expressed as follows (Leyla, 2011).
where A is the design matrix and contains the kernel of the radial basis function, is a square matrix according to the number of reference stations and is defined as follows: And the matrix P is defined as follows: There are different kernels for radial basis function interpolation, in this paper we select the inverse multiquadric (IMQ) interpolator kernel.

Inverse Multiquadric:
IMQ method is based on radial basis functions.For n distinct data, The IMQ is expressed as follows: is the coordinates of the station where the interpolation is done, and (φ ,λ ,h ) i i i is the coordinates of the reference stations.
The IMQ interpolator kernel is defined as (Sturgill, 2009) where r is the Euclidean distance and c is the shape parameter determining the interpolation accuracy.The shape parameter is obtained by the cross-validation method (Sturgill, 2009).
Interpolation based on the IMQ method requires discrete data, therefore, the following equation must hold (Leyla, 2011): In this paper calculation is based on having the geodetic coordinates of the points from equation ( 16): ( ) where ss 2 0 A sin 2 besides R is half the greatest length of the earth and s are the data points.
The position of the control points, which in this paper is in the form of longitude, latitude, and altitude, and the information, is expressed as follows: The unknown coefficients  ,  of equation ( 13) are calculated via the following equation: Then by using equation ( 13), PWV values in other places can be obtained.

Artificial neural network
ANN is an information processing system that mimics the brain's biological neural network by connecting many artificial neurons.Artificial neural networks can train and generalize.The work of a neural network is to create a specific output pattern according to a learning process based on the input pattern provided to the system.There are different types of neural networks.Perceptron neural network is used in this paper.because it is mainly used for interpolation (Zhang et al., 1998).

Perceptron neural network:
Multilayer perceptron neural networks (MLPNNs) is the most common type of neural network due to their fast performance, simple implementation process, and small training set requirement (Kocyigit et al., 2008).According to Figure (1), MLPNN consists of three consecutive layers: input, hidden, and output layers.The hidden layer processes the input information and sends it to the output layer.There is no way to determine the number of neurons in the hidden layer.Therefore, the number of neurons is selected only by trial and error.If the network does not converge to the desired value, the number of hidden layer neurons are increased (Kocyigit et al., 2008, Subasi, 2005).Perceptron neural network training is based on the backpropagation algorithm.The weights in perceptron neural network training are adjusted in a way that the error between the desired output and the output is reduced.
The values of the main inputs are actual values, where our inputs are the three variables of latitude, longitude, and geodetic height.According to the importance the entries are weighted.Each neuron j in the hidden layer sums its input signals xi impinging onto it after multiplying them by their respective connection weights wji.The output of each neuron is described as follows (Orhan et al., 2011, Sharifi, 2016): where f is the activation function using the weighted sum of inputs, ReLU activation function is used in this paper, is as defined as by Eq (20).and is of the form shown by

Inverse distance weighted (IDW)
One of the most common spatial interpolation methods in earth sciences is the IDW method.In this method, the position of each point is considered separately, and the relative position of the points is not considered.In this method, the weight is a function of the inverse of the distance, and the closer points have a more significant effect on the estimation of the unknown point so that the points with the same distance from the point with unknown values are given the same weight.The closest points are given highest weight, and it follows as: (Shepard, 1968, Sarma, 2010).
where w i are the weights and are obtained from the following equation: where r is the Euclidean distance between the sample stations and the location to be interpolated.And n is the number of reference stations.In addition, four stations of this network are considered as control stations to check the accuracy of interpolation.The details of these stations are shown in table (1).Besides, in this paper, the reanalyzed data of ERA-5 is used for accurate estimation of surface pressure and surface temperature for PWV calculation.The temporal resolution of ERA-5 data is 1 hour, and its spatial resolution is about 30 kilometers and provides data in 37 pressure layers.

PROCESSES AND RESULTS
In Bernese software, the data of GPS stations have been processed by the PPP method, and ZTD has been calculated in each station with a time resolution of 1 hour.Besides ZHD has been calculated according to the latitude and atmospheric pressure obtained from ERA-5 data at the stations using the Saastamoinen dry model.Then, with the difference of ZHD from ZTD, the troposphere delay is estimated, and using equation ( 3), PWV has been obtained for each station in millimeters.For interpolation by IMQ, using the crossvalidation method, the optimal shape has been determined so that RMSE be having the smallest value.RMSE values according to different shape parameters are shown in Figure ( 4).
The shape parameter 5000 has the lowest RMSE value and is selected as the optimal shape parameter.Then, the Euclidean distance between the reference stations and the control stations is calculated.Using the design matrix formed by the IMQ interpolator kernel, the unknown coefficients in equation ( 11) have been obtained.Finally, the interpolation has been done by the IMQ method for the control stations.By changing the number of hidden layers and trial and error, the number of hidden layers with a good answer is obtained.The number of neurons is also selected only by trial and error.Having too few neurons can reduce the stability of the system.Using more hidden layers, instead of increasing the number of hidden neurons, increases the processing power of the system.So, in this paper, an input layer, two hidden layers, and an output layer are used for perceptron neural network interpolation.Furthermore, the activation function of both layers is considered as the ReLU function.Using this function, modeling time is reduced.The first hidden layer has twelve neurons and the second hidden layer has eight neurons.30% of  1538.67 122.217349 46.17813 P700 1197.74 120.724609 46.950925 SC00 183.25 119.544194 45.953335 P450 .74 120.724609 46.950925 SC00 183.25 119.544194 45.953335 P450 the data have been selected as test data, and 70% of the data have been selected as training data.The inputs, which include the latitude, longitude and altitude of the stations, are multiplied by their corresponding weight, and their sum is In the training step, the difference between the actual PWV values and the desired PWV is calculated.In this step, the weights are modified, and the network is trained.Finally, the latitude and longitude, and height of the control stations are entered into the system to obtain the PWV values.In the next step, the interpolation is done by the IDW using the data weighting based on the Euclidean distance between the reference stations and the control stations and using the mentioned equation ( 21).RMSE and relative root mean square value (RRMSE) are the evaluation criteria and comparison between the interpolation methods performed in this paper as follows: where R PWV i is the real value and PWV the interpolated value at each control station.And n is the number of stations.Table (2   In the last step, a dense PWV map is prepared using the IMQ method in the region.For this purpose, interpolation points have been selected on a regular grid in the area with proper distribution and a distance of about 6 km from each other.Then, interpolation using the IMQ method has been used to estimate PWV values at forecast points.Dense maps of PWV to investigate PWV changes with the best accuracy in 3 different hours of the day are displayed in Figure ( 7).The amount of PWV is changing during the day.According to the figure (7), at 00:02, the amount of PWV is from 10.1 to 19.1 mm, at 14:00 from 3.2 to 21.2 mm, and at 22:00 from 3.7 to 24 mm.

CONCLUSION
In this paper, to improve the accuracy of PWV estimation in the troposphere, The IMQ method was used and compared with ANN and IDW methods.Four stations among the points in the area were considered as control stations.The RMSE value obtained for the IMQ method is 2.11 mm less than other methods, which means that the IMQ method has better accuracy than other methods.The RRMSE value obtained for the IMQ method is 16/9% mm, which is 9% better than the ANN method and 17% better than the IDW.As a result, the IMQ interpolation method is more suitable for sparse data, and the dense map created by the IMQ method shows the PWV changes very well.

Figure 2 .
Figure 2. ReLU activation function (Dubey and Jain, 2019) 3. STUDY AREA AND DATAAn area of 15448 square kilometers was selected in the state of Washington, USA.In this paper, 23 GPS stations of the UNAVCO network, a research institute related to earth sciences, have been used.The location and distribution of the stations are shown in figure (3).

Figure 3 .
Figure 3. Map of the study area and the satellite positioning reference stations

Figure 4 .
Figure 4. Determination of the shape parameter ) shows the comparison of the results.The results obtained from RMSE of the interpolated models with a resolution of 1 hour are shown in Figure (5).Also, the comparison between RMSE between PWV obtained from GPS and PWV obtained from the IMQ interpolation method in each control station is shown in figure (6).

Figure 6 .Figure 5 .
Figure 6.Comparison between PWV obtained from GPS signals and PWV obtained from IMQ interpolation of the four control stations

Table 1 .
Names and coordinates of control stations

Table 2 .
Statistical comparison according to the results obtained