An improved GNSS remote sensing technique for 3D distribution of tropospheric water vapor

Water vapor plays an extremely important role in the monitoring and prediction of weather, and GNSS tomography can obtain 3D spatiotemporal change information and reliable water vapor profiles. In this paper, an improved global navigation satellite system (GNSS) tropospheric tomography technique using an ERA5 (the fifth generation ECMWF reanalysis) product is developed. First, the ERA5 product was adopted to analyze the spatiotemporal distribution of water vapor, and a water vapor density threshold defining the top of the tomography was determined; then, the height of the grid top (GT) of different seasons was obtained through linear fitting; finally, the water vapor value between GT and tropopause is constrained using the data of the ERA5 product as the initial value. The new method for using the ERA5 product to determine the height of the GT of the tomographic grid reduces the height of the top layer of the grid and increases the number of effective GNSS rays. Data from nine CORS stations in Hong Kong in 2019 were selected for experiments. The results showed that the new algorithm increased the number of effective satellite signals by 14%. In addition, the ERA5 data, the radiosonde data, and the COSMIC‐2 data were used as reference values. The accuracy of the water vapor density obtained by the algorithm was improved by 25%, 17% and 9%, respectively.


| INTRODUCTION
Water vapor is a crucial component of the atmosphere that undergoes significant changes and exerts a critical impact on the climate. Accurately detecting the 3D distribution and variation of water vapor is essential for weather forecasting and climate research. Traditional water vapor detection methods, such as radiosonde, microwave radiometer, infrared radiometer, and satellite remote sensing, have limitations such as low spatiotemporal resolution, high cost, and susceptibility to temperature changes. GNSS tomography developed in recent years not only makes up for the shortcomings of the traditional water vapor detection technology but also has the advantages of real-time, high-precision, all-weather detection, and other benefits (Chen & Liu, 2014;Dong & Jin, 2018;Flores et al., 2000).
Since the proposal of GNSS meteorology (Bevis et al., 1992), the accuracy of precipitable water vapor (PWV) obtained by GNSS tomography was first evaluated (Emardson et al., 1998;Rocken et al., 1993;Rocken, Anthes, et al., 1997), which verified the feasibility of GPS detecting water vapor. However, PWV can only reflect the 2D distribution of water vapor over the station and cannot reveal the vertical distribution characteristics of water vapor. With the processing of tomography technology, GNSS tomography can be used to invert the 3D change information for water vapor (Champollion et al., 2005;Emardson & Webb, 2002;Flores et al., 2000;Troller et al., 2002). GNSS water vapor tomography has developed in the last 20 years. Scholars worldwide have conducted extensive studies on several key technical problems in tomography technology, solved the rank deficiency problem of tomography observation equation (Flores et al., 2000;Song et al., 2006;Yu et al., 2010), and proposed different tomographic grid division methods (Chen & Liu, 2014;Song et al., 2006) as well as solution methods for tomographic equations (Bender et al., 2011;Perler et al., 2011).
At present, scholars worldwide focus on how to improve GNSS water vapor tomography and further improve the accuracy of water vapor tomography. The accuracy of wet delay information, the solution for the tomography equation, the division of tomographic grid, and the number of effective rays are four key parts in tomography accuracy. Wet delay information is generally obtained by subtracting the hydrostatic delay information from the total delay, and high-precision hydrostatic delay information can improve the accuracy of wet delay information Ye et al., 2016); additionally combining other observation data or reducing the error in signal propagation can also further improve the accuracy of wet delay information (Heublein et al., 2019;Möller & Landskron, 2019). In terms of solving tomographic equations, scholars have analyzed and improved the algebraic reconstruction method to improve the speed and accuracy of solution (He et al., 2015;Xia et al., 2013), and some scholars have also used the compressed sensing method to solve the formula (Heublein et al., 2019). A large number of experiments have proved as well that high-precision prior information can help improve the vertical and horizontal constraints and improve the accuracy of tomography (Benevides et al., 2018;Chen & Liu, 2016;Xia et al., 2018). Tomographic grid division is the current research focus, and its core concept is to reduce the number of voxels that no signal rays cross, thereby increasing the stability of tomographic equations. Nonuniform height stratification and tomographic boundaries combined with historical information are more in line with actual water vapor distribution Ye et al., 2016;Zhang et al., 2020). Some scholars have also proposed a dynamic tomographic area division method, which can reduce the number of grids at the bottom of the tomographic grid, making it possible to perform tomography at several stations (Ding, Zhang, Wu, Wang, Kealy, & Zhang, 2018;. Additionally, tomographic accuracy can be improved by increasing the number of effective rays. By introducing virtual signals, new height factors, and other methods, scholars have enabled side rays participation in the composition of the tomography equation, thereby increasing the accuracy of the tomography solution (Zhao et al., 2018. Furthermore, multisystem GNSS tomography greatly increases the number of effective rays compared with single-system tomography. The proposed method in this paper aims to improve the accuracy of GNSS water vapor tomography by reducing the number of tomographic unknowns and increasing the number of effective GNSS rays. To achieve this, the height of the top layer of the tomographic grid is reduced based on the temporal and spatial distribution information of water vapor obtained from the ERA5 product. In the actual tomography, the water vapor value from the top of the grid to the tropopause was tightly constrained by the ERA5 product. Continuous operational reference system (CORS) data from Hong Kong and the ERA5 product in 2019 were selected for experiments to verify the effectiveness of the new method.
The structure of this paper is divided into five sections. Section 1 is the introduction. Section 2 introduces the basic principle of tomography, including the tomography model and the tomographic grid division. Section 3 describes the data used in the experiments, as well as the determination of the grid top (GT) height, and evaluates the feasibility of using GT. Section 4 presents the comparative evaluation and analysis of the experimental results with the ERA5 product, radiosonde product, and COSMIC-2 product. Finally, Section 5 contains the conclusion.

| GNSS water vapor tomography model
In GNSS tomography, the first step is to obtain the zenith total delay (ZTD) through GNSS. ZTD includes two parts: zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD). The ZHD can be obtained by meteorological parameters and empirical models, and the ZWD related to water vapor can be obtained after ZHD is obtained. Then, one can use the mapping function to project ZWD to the ray direction of the satellite signal to determine the slant wet delay (SWD) and finally convert the SWD to the slant water vapor (SWV) through the wet conversion factor Π.
By constructing a tomographic grid, calculating the distance of the satellite signal crossing the grid, and combining the observation value matrix on the signal path, the observation equation for solving the water vapor density can be formulated. However, due to the geometric distribution characteristics of stations and satellites in tomographic grids, satellite signal rays cannot pass through all grids, which leads to the fact that the original GNSS tomographic equation is rank-deficient; therefore, it is necessary to add additional tomographic equation constraints. In the horizontal direction, considering that the horizontal distribution of water vapor is stable and continuous in a certain small range, the water vapor density of the solved grid is constrained with the surrounding grids using Gaussian distance weighting (Song et al., 2006): where ω is the constraint coefficient and x is the water vapor density. In the vertical direction, because the water vapor content gradually decreases with increasing height, a negative exponential function is used to constrain different heights: where h is the height of the tomographic grid and H is the water vapor scale height. From this, we can obtain the tomographic equation in the classical form: where A is the distance matrix of the satellite signal rays passing through the tomographic grid, W and V are the coefficient matrices of horizontal and vertical constraints, respectively, X is the water vapor density matrix to be solved, B is the SWV matrix, and Δ is the observation noise vector.

| Division of tomographic grid
The upper and lower boundaries of water vapor tomography height are tropopause and station height, respectively. In the range from a certain height below the tropopause to the tropopause, water vapor content is relatively small. In the process of using the tomographic model, the error source is often flooded with valid information, resulting in a negative value for the solution result in this region. Therefore, the height of the boundary layer on the grid must be re-calculated during the actual solution process. In this article, we defined this height as the GT. First, the initial grid top (IGT) based on the high-precision atmospheric product must be obtained, and the water vapor above the IGT was negligible because it was extremely small. Then, combined with the characteristic that the water vapor is mainly concentrated below 6 km, and the water vapor content above 3 km fluctuates more than the water vapor content below 3 km , a nonuniform division method was adopted for the vertical grid: taking 3 and 6 km as two nodes, the grid spacing of 3 km was 500-600 m, the grid spacing between 3 km and 6 km was 700 m-900 m, and the grid spacing above 6 km was greater than 1 km. From this, we could extrapolate the vertical grid under IGT. Figure 1 shows the division of vertical grid. Then, the actual grid top (AGT) was obtained through spatial-temporal analysis of water vapor on the reliable atmospheric product, and the AGT was combined with the vertical grid under the IGT, and thus the grid divided in this study was obtained. The water vapor under AGT was calculated by tomography, and the water vapor on AGT was calculated by reliable atmospheric product. Therefore, the water vapor density parameter that must be solved in traditional water vapor tomography can be divided into two parts: where X 0 is the water vapor density between AGT and IGT, calculated using atmospheric product, and X u is the unknown water vapor density parameter that needs to be solved. Therefore, the tomographic equation can be converted into the following formula: In this way, AGT reduces the height of tomographic inversion, increases the number of rays that completely cross the tomographic grid, and reduces the number of grids that have not been crossed by signals, which can improve the quality of tomographic inversion.

| DATA AND EXPERIMENTAL SCHEME
The determination of AGT requires reliable water vapor distribution information, so this section first introduces the atmospheric product used in this study, as well as the experimental data and experimental areas, then introduces the method and process of determining AGT, obtains the final experimental scheme, and finally analyzes the feasibility of AGT.

| Data description
The European Centre for Medium-Range Weather Forecasts (ECMWF) provides numerical weather reanalysis data on a global scale, which can be used to analyze water vapor content. ERA5 is the latest generation of reanalysis product after ERA-Interim, and its accuracy is further improved compared with ERA-Interim (Ssenyunzi et al., 2020;Zhang et al., 2019). Radiosonde can measure the wind direction, wind speed, temperature, humidity, and pressure at different heights. Its maximum height can typically reach 30 km, and it is the main tool for high-altitude meteorological observation. The Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) is a program jointly implemented by the United States and Taiwan, which can provide near real-time information about the Earth's atmosphere, including temperature, pressure, and water vapor (Rocken, Van Hove, & Ware, 1997). COSMIC-2 is a newly launched satellite constellation in June 2019, building on the success of COSMIC-1. Numerous studies have shown that the ERA5 product, radiosonde product, and COSMIC-2 product can obtain high-precision atmospheric information (Chen & Liu, 2016;Gui et al., 2017;Li et al., 2003;Ssenyunzi et al., 2020;Wang & Zhang, 2008;Zhang et al., 2019;Zhao et al., 2019) that can be used to assist in and analyze ground-based GNSS tomography (Benevides et al., 2018;Xia et al., 2018;Zhao et al., 2018).
Ground-based GNSS observations were taken from the Hong Kong Satellite Positioning Reference Station Network (SatRef ), and nine stations containing meteorological parameters were selected. Figure 2 shows the station distribution and horizontal grid division selected in this study. The horizontal grid division range is E113:85 -E114:40 and N22:20 -N22:60 . The longitude and latitude resolutions are 0:055 and 0:04 respectively.
F I G U R E 2 Horizontal grid division and station distribution. The black triangle is the SatRef station, and the grid line is the horizontal grid division boundary.
F I G U R E 3 Temporal and spatial distribution of water vapor density in Hong Kong in 2019. The horizontal coordinate is the water vapor density (unit: g/m 3 ), and the vertical coordinate is the vertical height (unit: km). The black curve is the annual average, and the red curve is the fluctuation range.

| Determination of GT
In order to determine GT, this study first selected the ERA5 product from Hong Kong in 2019 to analyze the spatial and temporal distribution of water vapor in this region. Figure 3 shows the temporal and spatial distribution of water vapor density in Hong Kong in 2019 calculated by the ERA5 product. It can be seen from the figure that water vapor content above 10 km is very small. And relevant scholars have used long-term radiosonde and COSMIC historical data to obtain the tropopause of Hong Kong as 10 km Zhao et al., 2018), so the IGT was set as 10 km. At the same time, some scholars have selected 8 and 8.5 km as the tropopause in Hong Kong . However, we found that the fluctuation range of water vapor density above 8 km was much larger than the average value. And according to the original ERA5 data, the height containing data between 8 km and 10 km is about 8.5 km. At the same time, the water vapor density fluctuation range obtained from the ERA5 data at 8.5 km is at most 1 g=m 3 . Therefore, in this study, the water vapor density was 1 g=m 3 as the threshold, and the height corresponding to the whole year was fitted. The fitting function adopted a polynomial function, and the fitting criterion is the least square criterion. Figure 4 shows the distribution of the original height sequences and the fitting curves obtained under the quadratic function to the 11th power function.
It can be seen in Figure 4 that the overall variation trend is in line with the seasonal variation in water vapor content. The corresponding height is 4-5 km in winter, about 6 km in spring and autumn, and about 7-8 km in summer. At the same time, as the power of the fitting function gradually increases, the curve of the fitting function gradually tends toward consistency. After the calculations, it was found that when a function higher than the 11th power was used to fit, there would be an overfitting effect. Figure 5 shows the statistical results of the sum of squares for the residuals of each fitting function. The abscissa is the sum of squares for the residuals of different fitting functions, and the ordinate is the sum of squares for the residuals.
It can be seen from Figure 5 that the sum of the squares for the residuals from the quartic function to the eleventh power function is significantly smaller than other functions and is relatively close, so the quartic function image was selected as the fitting result in this study. Then, combining the seasonal variation law obtained in Figure 4 and the grid division method under IGT in Figure 1, the AGT in different seasons was obtained. The ERA5 data from 2016 to 2019 were used to fit AGT. And taking into account the division method of the vertical grid, the final AGT was determined. Table 1 shows the AGT in different years and seasons, finally with a minimum of 4.4 km in winter, 6.1 km in spring and autumn, and a maximum of 7.3 km in summer.

| Verification and experimental scheme
After obtaining the AGT and before performing tomographic analysis, we needed to calculate the SWV value on the AGT. Therefore, this study first used the ERA5 product to calculate the SWV value above the AGT for feasibility analysis, and the calculation results were compared and analyzed using radiosonde data as the F I G U R E 4 Function fitting diagram of different degrees. The horizontal coordinate is the annual cumulative day sequence in 2019 (day), and the vertical coordinate is the vertical height (km). The blue dots are the original height sequences, and the different colored lines are different fitting functions. true value. Table 2 gives the statistical results of the calculations, including bias, mean absolute error (MAE), and root mean square error (RMSE).
It can be seen from the table that the maximum deviation value is about 4.5 mm, the MAE is less than 0.6 mm, and the RMSE is less than 0.8 mm, indicating that the SWV value calculated using the ERA5 product above the AGT and the SWV value calculated from the radiosonde data have a strong correlation. Consistency also demonstrates the feasibility of using ERA5 data to constrain water vapor density above the AGT.

| RESULTS AND DISCUSSION
For the AGT obtained in different seasons, this study carried out experimental verification with the optimized method. At the same time, in order to test the improvement effect after reducing tomographic height, the traditional method of IGT was used as the control experiment. GNSS data (GPS and BDS) from nine SatRef stations containing meteorological sensors were used from Hong Kong (Figure 2). The sampling rate is 30 s and the processing software is GAMP. The Saastamoinen model and meteorological parameters were used to obtain ZHD, and the NMF mapping function was used to calculate SWD, and finally SWV was calculated by wet conversion factor. In this section, the number of effective ray signals under the AGT and IGT was counted and analyzed, and finally the tomography results were compared and discussed using the ERA5 product, radiosonde product, and COSMIC-2 product.

| Analysis of the number of satellite signals
Among most current tropospheric tomography studies, only satellite signals that pierced from the top of the tomographic grid into the range of the tomographic grid can be used for tomographic calculation. The number of available satellite signals determines the number of equations in the tomographic solution. The more the signals, the more the equations that can be solved, and the more accurate the tomographic results. Therefore, firstly, the number of effective rays with the optimized method proposed in this paper and the traditional method were F I G U R E 5 Sum of squares for residuals of different fitting functions. The vertical coordinate is the sum of squares for the residuals of different fitting functions, and the horizontal coordinate is the sum of squares for the residuals.
T A B L E 1 AGT in different years and seasons (unit: km). counted. Figure 6 shows the sequence diagram of the number of effective satellite signal rays between the optimized method and the traditional method in different seasons. It can be seen from the figure that the red line is significantly higher than the blue line; that is, the number of effective rays using the optimized method is larger than that using the traditional method. Table 3 also gives the specific statistics of the number of effective rays in different seasons, including the average number and the improvement rate. It can be seen from the table that in the four seasons, the optimized method has an improvement of more than 10% compared with the traditional method. It increases by about 18% in winter when the AGT is the lowest, while the improvement effect decreases in spring, autumn, and summer with the increase of AGT.

| Comparison of ERA5 data
The tomographic results below the AGT can also use the ERA5 product as reference data. Table 4 presents the MAE and RMSE statistics between the two methods with the ERA5 product. It can be seen from the table that the optimized method is a significant improvement in the overall tomographic accuracy compared with the traditional method. Winter has the best improvement effect, with an RMSE reduction of about 38%, followed by spring and autumn, around 27% and 25%, and summer has a slightly worse improvement effect of about 24%. Figure 7 shows the RMSE statistical values at different tomographic heights when compared adopting the ERA5 product. It can be seen from the figure that compared with the traditional method, the optimized method in this paper effectively reduces the RMSE value at each tomography height.

| Comparison of radiosonde data
With the use of radiosonde data, water vapor density profile with a more precise and higher vertical resolution can be obtained, which can be used as the control value to evaluate the accuracy of the tomographic water vapor density results.   Figure 8 shows the water vapor density profiles for the two methods in the four seasons in Hong Kong in 2019 and also gives the water vapor density profiles derived from the radiosonde data at the corresponding time. It can be seen from the figure above that both methods can invert the distribution and change of water vapor density effectively, but the red curve is closer to the black curve than the blue curve, especially at the bottom of the tomographic grid. This shows that the optimized method proposed in this paper can effectively improve the accuracy of bottom layer tomography. Table 5 also contains statistics for the MAE and RMSE between the tomography results with the radiosonde data. It can be seen from the table that the optimized method proposed in this paper leads to a certain reduction in MAE and RMSE compared with the traditional method, and the overall decrease is around 17%, which is in line with the conclusion of 4.2 on the improvement effect of different seasons. Figure 9 shows the RMSE values at different heights when compared using the radiosonde data. As can be seen from Figure 9, when compared using the radiosonde data, the optimized method proposed in this paper also significantly improves the tomographic accuracy of each layer. Above the AGT using the ERA5 product, its RMSE value is much smaller than that with the traditional method.

| Comparison of COSMIC-2 data
By acquiring and calculating occultation data in Hong Kong, high-precision water vapor density values can be obtained. We selected the occultation events that occurred in or near Hong Kong and carried out the corresponding time experiments according to the time of the occultation events. Figure 10 presents the water vapor profile and RMSE statistics at different heights compared adopting COSMIC-2 data.
It can be seen from the figure that compared with those of the traditional method, the results of the optimized method are also closer to the water vapor density profile of COSMIC-2, and the RMSE values at different heights have also improved to a certain extent. At the same time, we also found that the water vapor density value in the bottom layer has a large gap with the COSMIC-2 data, and the RMSE value of the water vapor density in the bottom layer is much larger than that in other layers. As a result, we analyzed the vertical extent of the COSMIC-2 data. Figure 11 shows the statistical graph for the lowest height in the vertical range of COSMIC-2 data.
It can be seen from the figure that the vertical lower limit of COSMIC-2 is higher, much larger than that in F I G U R E 7 RMSE statistics at different heights compared using ERA5 data. The vertical coordinate is the RMSE value and the horizontal coordinate is the height of different layers. Blue shows the traditional scheme of using IGT, and red shows the optimized method of using AGT. the bottom layer of the grid divided in this study. Ho and Chen analyzed the COSMIC-2 data and found that COSMIC-2 was significantly improved compared with COSMIC-1 (Chen et al., 2021;Ho et al., 2020). They found the penetration rate of COSMIC-2 below 1 km has reached 80%, but there are still 20% that cannot reach below 1 km, which is also in line with the statistical results given in Figure 10. They also pointed out that compared with the radiosonde data, COSMIC-2 data will have a large deviation below 2 km. Therefore, Table 6 shows the statistical results of the overall water vapor density for scheme 1 with the lowest layer of the tomographic grid and scheme 2 without the lowest layer. As can be seen from the table, compared with those under scheme 1, the MAE value and RMSE value under scheme 2 are reduced by 22% and around 26%, respectively. It can also be seen that under the two statistical schemes, the optimized method still has an 11% improvement in the MAE value and a 9% improvement in the RMSE value compared with the traditional method.

| Discussion
With the comparison among the ERA5 data, the radiosonde data, and the COSMIC-2 data, it is obvious that the proposed tomography method can significantly improve the accuracy of water vapor tomography, especially the accuracy of water vapor density in the middle and lower parts of the tomographic grid. For different reference F I G U R E 8 Water vapor density profile diagrams on experimental days in different seasons. The vertical coordinate is height, and the horizontal coordinate is water vapor density. The black curve represents the radiosonde data, the blue curve is the result of the traditional method, and the red curve means the result of the optimized method.
T A B L E 5 Statistics compared using radiosonde data (unit: g=m 3 ). F I G U R E 9 RMSE statistics at different heights compared using radiosonde data. The vertical coordinate is the RMSE value and the horizontal coordinate is the height of different layers. Blue shows the traditional scheme of using IGT, and red shows the optimized method of using the AGT.

MAE RMSE
F I G U R E 1 0 Water vapor density profile and RMSE statistics at different heights compared using COSMIC-2 data. The vertical coordinate is the height, the left horizontal coordinate is the water vapor density, and the right horizontal coordinate is the RMSE value compared using COSMIC-2 data. data, there are two main reasons for the inconsistent accuracy. First, ERA5 data was selected as the initial iterative value, so accuracy would be higher when using the ERA5 data as the reference data. Second, in order to increase the available occultation data, part of the occultation data outside the tomographic area was adopted when screening the available occultation data, so that accuracy was the worst when using COSMIC-2 data as the reference data.
But no matter what kind of reference data, the optimized method proposed in this paper has obvious improvement compared with the traditional method.

| CONCLUSION
Based on the principle that increasing the number of effective rays can improve tomographic accuracy, this paper proposes an improved tomographic method. Through the adoption of reliable atmospheric product to reduce the height of the tomography layer, the number of effective rays is increased and tomographic accuracy is improved. The main contents and contributions of this paper are as follows: The first is to analyze the temporal and spatial distribution of water vapor in the tomographic range through the ERA5 product, determine IGT, and conclude that water vapor content in the top layer fluctuates greatly. Then, the water vapor density of the top layer is counted, and the water vapor density reference value for calculating the AGT is obtained. Finally, taking the obtained water vapor density reference value as the threshold, the height change of the threshold is obtained by fitting with least square criterion and polynomial functions. It is found that the change is seasonal and repeated, with the lowest in winter, similar in spring and autumn, and the highest in summer. Based on this, this study obtains the AGT in different seasons. The second is to use the measured satellite data and meteorological data of SatRef in Hong Kong to carry out experiments in different seasons, which proves the feasibility of this method. Firstly, by counting the changes in the number of available satellites, it is shown that this method can significantly improve the number of effective rays passing through the grid, so as to enhance the stability and reliability of the tomography equation. Secondly, based on the true value of water vapor density calculated from ERA5 data and radiosonde data, it is found that the method proposed in this paper can significantly improve the accuracy of tomography inversion compared with the traditional tomography method, and the RMSE value of retrieved water vapor density is reduced by about 20%; considering statistics for tomography results at different heights, it is found that the proposed method can significantly improve tomography accuracy at the bottom of the tomographic grid. Finally, the feasibility of the method proposed in this paper is further verified by comparison adopting COSMIC-2 data. It is also found that COSMIC-2 has the disadvantage of low accuracy at low altitude.
In addition, the proposed method for reducing the height of AGT requires reliable water vapor information. ERA5 product is selected for research in this paper because its spatial resolution is more in line with the requirements of medium-and small-scale tomography, and the results also prove its reliability.
AUTHOR CONTRIBUTIONS Ankang Long: Data curation (equal); formal analysis (equal); investigation (equal); methodology (equal); writingoriginal draft (equal); writingreview and editing (equal). Shirong Ye: Conceptualization (equal); methodology (equal); resources (equal); validation (equal); writingreview and editing (equal). Pengfei Xia: Conceptualization (equal); data curation (equal); methodology (equal); resources (equal); supervision (equal). F I G U R E 1 1 Lower boundary height sequences and mean value of COSMIC-2 data. The vertical coordinate is the lowest altitude of the COSMIC-2 data, and the horizontal coordinate is the epoch of the COSMIC-2 event. The blue dots and lines are the height sequences, and the red line is the mean value. UCAR for providing COSMIC-2 data. The authors also would like to thank the Survey and Mapping Office (SMO) of Lands Department, Hong Kong, for providing GNSS data.