Advanced Grid Model of Weighted Mean Temperature Based on Feedforward Neural Network Over China

Weighted mean temperature (Tm) is a key parameter in Global Navigation Satellite System meteorology. In this study, European Centre for Medium‐Range Weather Forecasts Re‐Analysis product with a spatial resolution of 0.5° × 0.5° from 1999 to 2018 was used to study the spatiotemporal behaviors of Tm in China. Decomposed by Fast Fourier Transformation, Tm and lapse rate (β) variations are highly latitude‐dependent and exhibit periodicities on annual, semi‐annual, and diurnal scales. Meanwhile, Tm keeps increasing at a rate of 0.25 K per decade across China. Based on these discoveries, this study build a new grid Tm model based on feedforward neural network (FNN) with a spatial resolution of 0.5° × 0.5°, known as Grid‐FNN model. FNN is applied to each grid point to compensate the residual error of the corresponding periodic functions. And the fitting accuracy at each grid is improved by the FNN algorithm. ERA‐Interim product with a spatial resolution of 0.4° × 0.4° and Radiosonde data in 2018 are used to validate the new model, and the accuracy of Grid‐FNN model is proved 8.6% and 10.9% better than GPT2w‐Tm model, respectively. The Grid‐FNN model also shows better performance than IGPT2w, GTm‐III, and GTrop model in autumn and winter and in high‐altitude regions.


Introduction
Tropospheric delay is a major error source in Global Navigation Satellite System (GNSS) positioning and navigation. Zenith total delay consists of zenith hydrostatic delay and zenith wet delay (ZWD). Bevis et al. (1992) proposed the concept of GNSS meteorology by deriving a formula converting ZWD to precipitable water vapor (PWV). Subsequent research showed good consistency between PWV and sounding observations. As a result, the GNSS technology has gradually become an effective method to monitor atmospheric water vapor ever since. In general, PWV can be calculated by multiplying ZWD with a parameter obtained from weighted mean temperature (T m ) (Chen, 1998). To accurately calculate T m , vertical profiles of temperature and water vapor pressure are required but difficult to obtain. Such demands further motivate the research interest of T m modeling. At present, mainstream T m models can be classified into two categories depending on whether meteorological parameters are used, such as surface temperature (T s ). Bevis et al. (1992) analyzed 8,718 records of sounding data from 13 radiosonde stations in North American. T s was found correlated with T m , so he suggested a linear regression formula for mid-latitude regions, known as the Bevis model. This model has been widely used in ground-based GNSS water vapor monitoring. However, after investigating 23-year data from 53 radiosonde stations around the world, Ross and Rosenfeld (1997) concluded that the relationship between T s and T m could vary temporally and spatially. Yao et al. (2015) also confirmed the existence of nonlinearity by studying radiosonde data. Wang et al. (2016) studied the applicability of Bevis model in various latitude and found biases. In tropical and subtropical regions, the Bias could fluctuate from −1 to −6 K. In middle and high latitude area, the magnitude could raise to 2-5 K. As a general conclusion, spatiotemporal impacts should be taken into consideration for Bevis model. Researchers have been studying T m and T s in specific regions to establish localized Bevis models using radiosonde data since 2000, which achieved fairly good results (Baltink et al., 2002;Bokoye et al., 2003;Jade & Vijayan, 2008;Liou et al., 2001;Y. Wang et al., 2007;Yi et al., 2017;Zhu & Hu, 2018). This kind of model can achieve fair accuracy but has limited applicability due to high dependence on surface meteorological parameters.
The second type of T m model is independent of surface meteorological parameters. Yao et al. (2012) established an empirical model named GWMT model in a global scale using radiosonde data of 135 stations in [2005][2006][2007][2008][2009]. Using multi-source data, Yao et al. (2013Yao et al. ( , 2014 further improved its accuracy and global applicability. The Global Pressure and Temperature 2 Wet (GPT2w) model was established based on the data of European Centre for Medium-Range Weather Forecasts Re-Analysis (ERA-Interim) fields with a horizontal resolution of 1° × 1°, which has been most commonly used. Users of GPT2w only need to input date and 3D coordination to obtain a T m estimation. The model searches the nearest four grid points then calculates T m by bilinear interpolation (Böhm et al., 2015). However, the elevation difference between the model's reference level and the station was somehow ignored. The Root Mean Square (RMS) error for China was 4.43 K in 2015 . Zhang et al. (2017) studied the influence of elevation on T m in four quarters of China to enhance two T m models.  proposed a new grid model known as GGT m , which took latitude and altitude changes into consideration. Its RMS in a global scale was 3.54 K, which was 8% higher than the GPT2w model.  developed the GPT2w model for China by taking the T m lapse rate (β) into consideration, which served as the vertical correction of T m and the adjusted model showed superior performance. Using the ERA-Interim reanalysis data from 1979 to 2017, Sun et al. (2019) proposed a new model (GTrop-T m ) considering annual and semi-annual T m fluctuations as well as T m spatial variations. The new model could provide higher accuracy of T m at high-altitudes. These models have accelerated the development of GNSS meteorology. However, it is necessary to optimize localized T m model for China due to the complexity of terrain and climate, in order to improve the PWV accuracy.
T m reflects the characteristics of tropospheric atmospheric profile. It can exhibit significant variations with altitude and latitude. Different annual and semi-annual periodic features have been identified at different latitudes, which can hardly be described simply by linear or periodic functions (Ding, 2020). The neural network (NN) model is a massive parallel distributed processor composed of simple processing units, which grants NN strong generalization ability (Hu, 2006). When dealing with complex nonlinear problems, it is likely to obtain higher accuracy than conventional approaches. Taking T s into consideration, some researchers improved the first category T m models by NN. Ding (2018) established a nonlinear mapping model of T m based on NN. The model only need to input value of T s , empirical T m , and the latitude to obtain an RMS of 3.3 K on a global scale. A NN model was developed in Ding (2020) to estimate T m from surface to the top of the troposphere, which showed 29.1% accuracy improvement compared to the Gtrop-T m model over different heights of the troposphere on a global scale as the T s was considered. These models proved the effectiveness of the NN algorithm. Therefore, this study aims to establish a second-type T m grid model based on NN algorithm.
In this study, periodic characteristics and long-term linear trend of T m and β are investigated using ERA-Interim products. Meanwhile, models of T m and β are established by combining the grid model and NN model, then validated by ERA-Interim products and radiosonde data. Finally, the main findings are summarized in order to demonstrate the effectiveness of the new modeling method.

Definition of T m and Its Lapse Rate (β)
T m can be approximated calculated as (Equation 1) (Davis et al., 1985) 1 where P v is the partial pressure (in hPa) of water vapor, T i is the atmospheric temperature (in Kelvin) at the ith level. P v can be calculated by Equation 2 (Ross & Elliott, 1996)   where RH is relative humidity. These parameters can be provided by the ERA-Interim product or radiosonde data.
Although there have been many models for T m calculation, vertical adjustment method is still desired because height differences exist between GNSS station and the selected reference level, especially in vast regions like China with complex terrain. To explore the lapse rate of T m , its profiles of each grid over 20 years (January 1, 1999∼December 31, 2018 derived from the ERA-Interim data are calculated over China. Based on previous research, a linear Equation 3 can be used to describe the relationship between T m and height (Zhang et al., 2017).
where β is the T m lapse rate (K/km), h is the geopotential height (km). T m0 is the value of T m at the reference level. In this study, reference level is the average sea level instead of the real ground surface, so the geopotential height at reference level is 0 m. T mh is the value of T m at stations of certain geopotential heights.

Data Sets
European Centre for Medium-Range Weather Forecasts (ECMWF) has accumulated reanalysis data of atmosphere and oceans for decades. The ERA-Interim product has better quality with higher spatial resolutions. It provides data of temperature, potential and relative humidity of 37 layers from 1,000 to 1 hPa, with different spatial resolution at 0000, 0600, 1200, and 1800 UT each day. ERA-Interim products have made remarkable achievements in multiple fields so far.
Radiosonde is a major data source of atmospheric detection. It can measure data such as pressure, temperature, and relative humidity, etc., at a frequency of twice per day with a radio telemeter and sensitive components on air balloons. In this study, ERA-Interim products with a spatial resolution of 0.5° × 0.5° during the period of January 1, 1999-December 31, 2017 in China are used for modeling. ERA-Interim data with a spatial resolution of 0.4° × 0.4° and radiosonde data of 78 stations in 2018 are used for validating.
T m can be directly determined from ERA-Interim by linearly interpolating data from four nearest grids surrounding the radiosonde station (Wang et al., 2016). In this study, T mh of the four grids are calculated using T m0 and β by Equation 3, and then T m over the radiosonde station can be obtained by bilinear interpolation. The correlation coefficient of two T m series is 0.99. The corresponding Bias and RMS are 0.70 and 1.54 K, respectively. The statistical parameters suggest that two T m series show no significant difference. Radiosonde data is believed to be more convincing because of direct measurement. However, the radiosonde observations are usually incomplete, and the distribution of stations is sparse and inhomogeneous, so it may be difficult to acquire an accurate T m trend. In contrast, ERA-Interim have a higher spatial-temporal resolution with rare missing data. Therefore, ERA-Interim product can be used to analysis the temporal and spatial varieties of T m and further modeling. Cooley and Tukey (1965) proposed an algorithm for calculating Fourier coefficients, which required much less computation compared with traditional discrete Fourier transformation. It is well known as the fast Fourier transform (FFT). FFT provides a method to study signals in frequency domain. It has long been used to analyze linear systems and identify frequency components. To analyze the temporal characteristics of T m0 and β, the study area is divided into six subregions according to the topography, as shown in Figure 1. 27°N and 40°N are considered as two critical latitudes. Between them lies the Tibet Plateau, which has an entirely different terrain. 105°E is treated as the boundary dividing east and west China. The topography referred to ETOPO1.

Analysis of T m0
The relationship between T m0 and T mh can be linearly expressed by Equation 3. As a result, T m0 and T mh share identical amplitude and phase. The frequency and amplitude figures of the six subregions obtained by the FFT algorithm are shown in Figure 2 with numerical details in Table 1.
It can be inferred that in all the six subregions, the time series of T m0 exhibit annual, semi-annual, and daily variations, as well as certain spatial patterns in amplitude.
To study the long-term behavior of T m0 , linear fitting equation (in the form of y = ax + b) is used in all subregions with the result shown in Table 2, where slope "a" represents the long-term change, and y represents T m0 . Comparing the p-value with the significance level  (0.05 used in this study), the significance of "a" can be determined.
The p-value of Region I(A) is greater than 0.05, indicating that the T m0 long-term trend is not significant. But the other five sub-regions are close to 0, so a long-term increasing trend exist in most Chinese region. For the spatial variation of "a", it increases with latitude in eastern China, while in the west it shows an opposite tendency. To conclude, analysis of T m0 series suggests obvious annual, semi-annual and diurnal periodic variations with certain geographical distribution. Considering the longterm behavior, the average T m0 in China has an increasing rate of 0.25 K ZHU ET AL.
10.1029/2020EA001458 4 of 13 where i a (i = 1,2…4), i b (i = 1,2, 3), and c are the coefficients of the periodic function. JD refers to the modified Julian date.

Analysis of the Lapse Rate (β)
FFT based temporal analysis of β is preformed and shown in Figure 3, with details tabulated in Table 3   To study the long-term behavior of β, linear fitting equation (in the form of y = ax + b) is used on all six subregions with their details listed in Table 4.
For lapse rate β, the p-values of all the western subregions are greater than 0.05, indicating that the linear relationship between β and time in these regions is not significant. p-values in all the eastern subregions were close to 0, so the linearity was significant. Slope "a" for  represents its long-term changes. In eastern China, it shows a value less than 0.04 K/km/decade, while in west China it's almost zero, which suggest little changes in long-term period. To conclude, β shows certain annual, semi-annual and daily cyclic changes. Thus, the periodic variations of β can be approximatively expressed by Equation 5. sin .
where a i  , b i  (i = 1, 2, 3) and c are the periodic function coefficients. JD refers to the modified Julian date.

Feedforward Neural Network (FNN)
In high latitudes regions, the relationship between T m and T s tends to be nonlinear (Mendes et al., 2000). T m can show a complicated periodic variety depending on the terrain and location. The type of NN model used in this study is feedforward neural network (FNN). It is widely used in function approximation, pattern recognition, classification and data compression. More details about FNN and its applications can be found in (Simon, 1999). In this study, FNN is applied to each grid point to improve the accuracy of T m0 and β estimation. A FNN model has a three-layer structure, which are input layer, hidden layer, and output layer. The back-propagation algorithm, which is the most popular learning technique, is used in the model training. Hyperbolic function, as shown in Equation 6, is used as the transfer functions of the neuron in the hidden layer.
FNN algorithm is implemented using the toolbox of Matlab R2018b.

The Grid-FNN Model
From Section 3.1, it is obvious that T m0 and β have annual, semi-annual and daily periodical characteristics with spatial variation. Data from ERA-Interim from January 1, 1999 to December 31, 2017 is used to establish a Grid-FNN model with a resolution of 0.5° × 0.5°, as shown in Figure 4. Equations 4 and 5 are applied to each grid points to compute the approximate T m0 and β, while FNN algorithms are used to calculate the corresponding residual errors ( s  ) of the two equation.
Subfigure in the right red rectangle in Figure 4 is an enlargement of the FNN structure. The day of year (doy) and the target parameter, that is, the approximate T m0 or β calculated by Equation 4 or Equation 5 are treated as two neurons in the input layer. The hidden layer includes four neurons. The output layer is the corresponding  of the approximate T m0 or β. The modeling process mainly consists of two steps:

Determining the coefficients of Equations 4 and 5 by fitting the T m0
and β series at each grid point using least squared method; ZHU ET AL.   All the input parameters for Grid-FNN model are longitude (in degree), latitude (in degree), geopotential height (in meter) of the user and modified Julian Day, while the output is estimated T m . The Grid-FNN model operates with the following steps: 1. Based on the station's location, the model can locate its four nearest grids; 2. The approximate T m0 and β of the four grid points (the lower black circle in Figure 4) can be calculated with Equations 4 and 5 by input the Julian Day; 3. Input the converted date (Doy) from the modified Julian Day together with the calculated approximate T m0 or β in step two into each FNN model to get the corresponding s, then T m0 or β can be adjusted by add its ; 4. Calculate the T mh at the station elevation (the black dot in Figure 4) with Equation 3, then the station's T m can be acquired by linear interpolation.

Results
The ERA-Interim data with a spatial resolution of 0.4° × 0.4° and the radiosonde data of China in 2018 are used to validate the accuracy of Grid-FNN model. RMS and Bias are used as the accuracy criteria, which can be calculated by Equations 7 and 8.
where  i MP is the reference value from ERA-Interim data or radiosonde data, i MP is the value estimated by the model and N is the number of observations. Bias and RMS comparison is made between Grid-FNN model and the traditional GPT2w-T m model (Böhm et al., 2015), IGPT2w-T m model , GT m -III model (Yao et al., 2014), GTrop-T m Model (Sun et al., 2019)

Validation of the Model by ERA-Interim Data
The Bias and RMS validated by ERA-Interim product are shown in Table 5.
Judging from Table 5, the mean Bias in all five models is close to 0 K, which suggests that there are no systematic errors. The maximum Bias of GPT2w model is 14.86 K, while the minimum value is −12.93 K. The range is 27.62 K. In contrast, the range for IGPT2w, GT m -III, 7.41,4.66,and 3.92 K, respectively, which further prove the accuracy and stability of Grid-FNN model. For the RMS, GPT2w model has the greatest RMS in China. The other models with elevation adjustment show better performance than GPT2w. Considering the max and mean RMS, Grid-FNN model shows the minimum values. And the mean RMS of Grid-FNN is 10.9% improvement over GPT2w. Figure 5 shows the Bias distribution of the five models calibrated with ERA-Interim data.
For GPT2w model, the Bias is comparatively greater due to the lack of elevation adjustment. As a result, the positive and negative Bias show a complicated distribution in high altitude region. For IGPT2w and GT m -III models, positive Bias exist in the northwest and Tibet Plateau, as well as GTrop model. Figure 5e suggests that Grid-FNN model is the best among all five models. Figure 6 shows that the accuracy of Grid-FNN model increases from north to south. For GPT2w model, its accuracy shows similar pattern with Grid-FNN model but lower accuracy in western China due to its negligence of the elevation differences. In regions with high altitude, Grid-FNN model also performs better than other models.

Validation of the Model by Radiosonde Data
Bias and RMS of the five models from Radiosonde data in 2018 are shown in Table 6.
It can be seen from Table 6 that radiosonde data-based validation shows that the mean Bias of GT m -III model is 1.5 K, while other models are all smaller than 1K. For Grid-FNN, it is only 0.1 K, which denied the ZHU ET AL.  possibility of systematic error. Grid-FNN model also exhibit the minimum range of Bias among all models, which is 9.38 K, while GTrop and GPT2w are 10.5 and 14.64 K, respectively, which further prove the accuracy and stability of Grid-FNN model. Validating with radiosonde data, the mean RMS of Grid-FNN shows 8.6% improvement over GPT2w. When comparing with other models, Grid-FNN model has the highest accuracy, and it can be inferred that FNN algorithm can effectively improve the accuracy of the model. Figure 7 is the Bias distribution of five model validated by radiosonde data. It can be seen that all models exhibit greater positive Bias at 113.30°E, 23.48°N. Figure 8 plot the RMS distribution of five models validated by radiosonde data. Similar trend is found in all the models, generally decreasing from north to south.
For GPT2w model, larger RMS is found in stations with elevations of 1-2 km and over 2 km. All models have greater RMS at 113.30°E, 23.48°N and 110.00°E, 23.56°N. For Grid-FNN model, its RMS is smaller than 5 K in 67% of the stations. Meanwhile, the percentages for GPT2w, IGPT2w, GT m -III and GTrop are 50%, 62%, 55%, and 63%, respectively. In general, Grid-FNN model shows the best elevation compatibility with an RMS of 4.47 K.  Table 7. It can be observed from Table 8 that all the five models exhibit low accuracy in winter and high accuracy in summer when validated by ERA-Interim data. In January, February, November, and December, Grid-FNN model is significant better than other models, while the other models only exhibit slightly improvement compared with GPT2w-T m model. In March, April, June, and September, Grid-FNN model shows little lower accuracy than GTrop, but still much better than GPT2w-T m . GT m -III model and IGPT2w model show the best accuracy in August. In general, Grid-FNN model has the best accuracy and stability.

Temporal Analysis of Bias and RMS of Grid-FNN Model
ZHU ET AL.

Accuracy Variation of the Model in Vertical Direction
In order to verify the accuracy variation with altitude, the height from 0 to 7 km is divided into seven layers (the interval is 1 km). The Bias and RMS of each layer validated by 2018 ERA-Interim data is shown in Figure 9: It can be seen from Figure 9 that except Grid-FNN model, the other four models show larger positive Bias in high area. Grid-FNN model exhibited negative Bias from 2 to 5 km. No significant RMS difference is found among all models from 0 to 1 km. But as altitude increased, GPT2w-T m model shows greater RMS due to the absence of elevation adjustment. Although IGPT2w, GT m -III and GTrop-T m consider such adjustment, the RMSs still show tendency of increasing where the elevation is greater than 5 km. Meanwhile, the accuracy of Grid-FNN remain stable. In all elevation except 3-4 km, Grid-FNN model exhibit the minimum RMS among all models.

Verification of the Grid-FNN Model With Different Hidden Node
The difficulty of the application of NN is that there are limited examples to follow to determine the neuron number in hidden layer. In this study, comparisons are made from 3 to 12 hidden layer neurons in order to discuss the accuracy of the model. Calibration is made with ERA-Interim, as shown in Table 9.
It can be seen from Table 9 that as more neurons are added into the FNN algorithm, the Bias and RMS change slightly, varying around −0.17 and 3.76 K, respectively. The high accuracy suggests good stability and effectiveness of FNN algorithm. Grid-FNN model weaken the influence of spatial factors on FNN approach, overcoming its stability and generalization, as well as the impact of number of hidden layer neurons. Meanwhile, FNN algorithm can compensate the residue of traditional periodic function in the grid points, improving its accuracy.

Conclusions
As a key parameter to transform ZWD into PWV, weighted mean temperature is crucial in GNSS meteorology. In this study, the periodic characteristics of T m0 and β are analyzed by FFT. 1. A linear trend of T m0 of 0.25 K per decade is discovered and is probably a consequence of global warming. However, no significant linear trend behavior is found for β in China. The variation is less than 0.04 K/ km per decade. 2. Based on the analysis of T m0 and β, a 0.5° × 0.5° Grid-FNN model without meteorological parameters is established. Compared with GPT2w-T m , Grid-FNN model increases the accuracy by 10.9% and 8.6% when validated by ERA-Interim with a resolution of 0.4° × 0.4° and radiosonde data in 2018, respectively. When compared with other models, Grid-FNN model also shows superior performance. 3. The accuracy of Grid-FNN model shows an increasing trend from north to south. Unlike other models, its accuracy do not decrease as the elevation increased. In January, February, November, and December, Grid-FNN model shows significant better performance than other models. These are the advantages of the Grid-FNN model among other models. 4. In this study, we first combine the grid model and the FNN model to build a novel T m model. A simple FNN is used to compensate the residual error of the periodic function at each grid point, while the grid model guarantee the stability of the network. The influence of the number of hidden layer neurons on the FNN network is also weakened. By using the FNN algorithm, Grid-FNN model improves the fitting accuracy of the traditional periodic function, and has obvious advantages in situations where the accuracy of traditional modeling methods is not satisfactory (in regions with high altitude or in autumn and winter).
In this study, we proposed a new grid model combining grid model with artificial intelligence algorithm to improve the fitting accuracy of periodic function. The Grid-FNN model shows superior performance comparing with other models, especially at the station with high altitude or in autumn and winter. The conclusions provide new modeling ideas which are beneficial to improve the accuracy of GNSS-based real-time PWV.
ZHU ET AL.