Ground-Based GNSS Precipitable Water Vapor Retrieval With the BeiDou B2B Service

1. Introduction

Water vapor is an important variable constituent of the global atmosphere. Accurate measurement of water vapor is essential for meteorology, hydrology and climatology on weather forecasting, disaster warning, and climate prediction [1,2,3]. Many observation techniques have been employed for the water vapor measurement. The traditional radiosonde is expensive, normally twice a day in low temporal resolution. Microwave radiometer and satellite remote sensing are often in poor accuracy in cloudy or rainy weather. However, the emergence of GPS-based PWV retrieval has gained increasing attention from scientists with all-weather capability, high accuracy, and cost-effectiveness [4,5] Many meteorological agencies have used the technique in the ground GNSS meteorological network solutions in the US, Japan, Germany and China, to name a few [6,7,8]. High quality water vapor observation from GNSS is validated with the accuracy of about 2mm, which can match with, even better than the Radiosonde [9,10].

Double-difference (DD) network solution and precise point positioning (PPP) are two widely used approaches to process ground-based GNSS measurements. The DD solution is a basic method used to calculate the precise positions or satellite orbits and other parameters from the global or regional GNSS network [11]. On the contrary, PPP fixes satellite orbits and clocks and assumes that the mathematical models are consistent with those applied in the processing of the global reference network used to derive the satellite products [12]. In the DD algorithm, receiver clock biases are eliminated and the integer ambiguity resolution is easy, while in the PPP and data screening for outliers and cycle slip detection and editing is more challenging. Thus the DD performs better in the systems stability and is wider used than the PPP in global and regional GNSS PWV analysis in recent years.

The DD network solution requires that all stations transfer their GNSS raw observations to the center for centralized processing. It is a challenging work to obtain high tempo-resolution of PWV with short time latency. With the high reliability and availability of satellite orbits and clocks products these years, the PPP technique acquire much progress in its accuracy and performance. It demonstrates its advantages of computational efficiency and timeliness over the DD network solution in data processing. With timely and high tempo-resolution PWV, severe very short weather in summer could be detected and targeted.

The BeiDou Navigation Satellite System (BDS-3) provides five public service signals of B1I, B1C, B2a, B2b, and B3I in three frequency bands of B1(central frequency 1575.42MHz), B2(1176.45MHz), and B3(1268.52MHz). Since 2020, China's BeiDou has begun to provide an initial real-time precise point positioning (PPP) service via B2b signals for the Asia–Pacific region. The service uses 3 GEO satellites (C59–C61) to broadcast corrections of orbit and clock offset for the BeiDou and GPS navigation satellites [13].

After its service, precise orbit determination [14,15,16], satellite clock offset estimation [17,18,19], tropospheric delay estimation [20], precise positioning [21,22,23,24] and time transfer [25,26,27,28], are studied. For precise applications, two representative technologies are Real-time Kinematic (RTK) and PPP [29]. But few studies are involved in the real-time water vapor retrieval with the BeiDou B2b Service. In the paper, an experiment was conducted and real-time PWV was retrieved. Validated by the radiosonde and ERA5, performance of BDS real-time PWV was assessed.

2. Materials and Methods

2.1. Experiment description

In order to test the performance of real-time GNSS PWV retrieval based on the BeiDou B2b service, a temporary BeiDou GNSS meteorological station , shown in Figure 1, was setup at the New Technology Base of Chinese Academy of Sciences, Dengzhuang South Road, Haidian District, Beijing on January 10, 2023. The station was located at 116.16 E in longitude , 40.04 N latitude and 94.78 m above the sea level, in very good observational environment with little line-of-sight obstruction. Configured with a BeiDou GNSS receiver AirPro, a choke ring antenna HX-CGX611A, a set of automatic weather station WTS400 and a set of BeiDou data transmission terminal, the observation lasted for 17 days.

2.2. The B2b real-time precise ephemeris of the Beidou and GPS satellites

Messages for the correction of satellite orbit , clock offset, differential code bias (DCB) and user range accuracy index (URAI) are broadcasted in the PPP-B2b signal. The Orbit corrections and clock offset are used to correct the broadcast ephemeris.

The orbit correction of a satellite is written as

$$\mathsf{\delta }\mathsf{\rho }={\left[\begin{array}{ccc}{\mathsf{\delta }\mathsf{\rho }}_{\mathrm{r}}& {\mathsf{\delta }\mathsf{\rho }}_{\mathrm{t}}& {\mathsf{\delta }\mathsf{\rho }}_{\mathrm{a}}\end{array}\right]}^{\mathrm{T}}$$

(1)

where $\mathsf{\delta }\mathsf{\rho }$ is the position correction, the subscripts r, t, and a denote the radial, the tangential, and the axial direction respectively. The superscript T is the transpose of the vector. The precise satellite position $X_C^S$can be corrected from the broadcast ephemeris $X_{brdc}^S$ with the following equation:

$$X_c^s=X_{brdc}^s-\left[e_r{e}_t{e}_a\right]\mathsf{\delta }\mathsf{\rho }$$

(2)

where $e_r=\overrightarrow{r_s}$/$\left|\overrightarrow{r_s}\right|$，$e_t=\overrightarrow{r_s}$·$\overrightarrow{v_s}$/$\left|\overrightarrow{r_s}\times \overrightarrow{v_s}\right|$，$e_a=e_t\times e_r$，and $\overrightarrow{r_s}$ is the vector of the satellite’s position, $\overrightarrow{v_s}$ is the vector of the satellite’s speed.

The satellite clock is corrected by this equation:

$$d{t}_{\mathrm{c}}^s=d{t}_{\mathrm{b}\mathrm{r}\mathrm{d}\mathrm{c}}^s-\frac{\delta {\tau }^s}{c}$$

(3)

where $d{t}_{\mathrm{c}}^s$ is the corrected clock offset, $d{t}_{\mathrm{b}\mathrm{r}\mathrm{d}\mathrm{c}}^s$ is the clock offset of broadcast ephemeris, $\delta {\tau }^s$ is the clock offset correction, and c is the speed of light in vacuum.

B2b correction message can be matched by two steps. In the first step the IODN (Issue of Data for Navigation) is used to select the corresponding broadcast ephemeris. For BDS-3 and GPS, the broadcast ephemeris here are navigation messages CNAV1 and LNAV1 respectively. In the second step the orbit and the clock corrections are matched by IOD CORR (Issue of Data for Correction).

Since the B2b correction messages can be delayed by propagation, the validity of the correction must be ensured. The nominal validity delays for the orbit correction, the clock offset correction, the DCB, and the URAI are 96s, 12s, 86400s and 96s respectively.

2.3. Real Time PPP for ZTD with BeiDou B2b Ephemeris

The PPP model is the dual frequency ionosphere-free linear combination (IFLC) model, in which the pseudorange and carrier observations can be expressed as

$$\begin{array}{c}{P}_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s={\rho }_r^s+\left(d{t}_r+d_{r,{\mathrm{I}\mathrm{F}}_{jk}}\right)-\left(d{t}^s+d_{,{\mathrm{I}\mathrm{F}}_{jk}}^s\right)+T_r^s+{\epsilon }_P \\ L_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s={\rho }_r^s+\left(d{t}_r+b_{r,{\mathrm{I}\mathrm{F}}_{jk}}\right)-\left(d{t}^s+b_{,{\mathrm{I}\mathrm{F}}_{jk}}^s\right)+T_r^s+{\lambda }_{{\mathrm{I}\mathrm{F}}_{jk}}{N}_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s+{\epsilon }_L\end{array}$$

(4)

where $P_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s$ $L_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s$are observations of pseudorange and carrier phase respectively. The superscript s represents the satellite, the subscript r represents the receiver, and j and k represent the signal frequency. $d{t}_r$ and $d{t}^s$ denote the clock offset of the receiver and the satellite. $T_r^s$ is the tropospheric delay, $d_{r,{\mathrm{I}\mathrm{F}}_{jk}}$and $b_{r,{\mathrm{I}\mathrm{F}}_{jk}}$ are the pseudorange and carrier hardware delay of the receiver. $d_{,{\mathrm{I}\mathrm{F}}_{jk}}^s$and $b_{,{\mathrm{I}\mathrm{F}}_{jk}}^s$ are the pseudorange and carrier hardware delay of the satellite. $N_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s$ and ${\lambda }_{{\mathrm{I}\mathrm{F}}_{jk}}$ are the ambiguity and the wavelength. ${\epsilon }_P$ and ${\epsilon }_L$ are the noises in pseudorange and carrier observations.

Since BDS broadcast their ephemeris based on B3I frequency, the corrected real-time satellite clock offset should include the hardware delay of B3I signal. $\mathrm{L}\mathrm{e}\mathrm{t} \mathsf{\delta }{\widehat{\mathrm{t}}}^{\mathrm{s}}=\mathrm{d}{\mathrm{t}}^{\mathrm{s}}+{\mathrm{d}}_{,\mathrm{B}3\mathrm{I}}^{\mathrm{s}}$，and ${\mathsf{\delta }\mathrm{t}}_{\mathrm{r}}=\mathrm{d}{\mathrm{t}}_{\mathrm{r}}+{\mathrm{d}}_{{\mathrm{r},\mathrm{I}\mathrm{F}}_{\mathrm{j}\mathrm{k}}}$, then

$$\begin{array}{c}{P}_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s={\rho }_r^s+\delta {\widehat{t}}_r-\delta {\widehat{t}}^s+T_r^s-\left(d_{,{\mathrm{I}\mathrm{F}}_{ij}}^s-d_{,B3I}^s\right)+{\epsilon }_P\\ L_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s={\rho }_r^s+\delta {\widehat{t}}_r-\delta {\widehat{t}}^s+T_r^s+{\lambda }_{{\mathrm{I}\mathrm{F}}_{jk}}{B}_{r,{\mathrm{I}\mathrm{F}}_{jk}}+{\epsilon }_L\end{array}$$

(5)

Where ${\lambda }_{{\mathrm{I}\mathrm{F}}_{jk}}{B}_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s={\lambda }_{{\mathrm{I}\mathrm{F}}_{jk}}{N}_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s+\left(b_{r,{\mathrm{I}\mathrm{F}}_{jk}}-d_{r,{\mathrm{I}\mathrm{F}}_{jk}}\right)-\left(b_{,{\mathrm{I}\mathrm{F}}_{jk}}^s-d_{,B3I}^s\right)$, and

$$\begin{array}{c}{d}_{,{\mathrm{I}\mathrm{F}}_{jk}}^s-d_{,B3I}^s=\left(\frac{f_j^2}{f_j^2-f_k^2}{d}_{,j}^s-\frac{f_k^2}{f_j^2-f_k^2}{d}_{,k}^s\right)-d_{,B3I}^s\\ =\frac{f_j^2}{f_j^2-f_k^2}{\mathrm{D}\mathrm{C}\mathrm{B}}_{,j}^s-\frac{f_k^2}{f_j^2-f_k^2}{\mathrm{D}\mathrm{C}\mathrm{B}}_{,k}^s\end{array}$$

(6)

${\mathrm{D}\mathrm{C}\mathrm{B}}_{,j}^s$ and ${\mathrm{D}\mathrm{C}\mathrm{B}}_{,k}^s$ can be directly used to correct the IFLC pseudorange observation because the corrections of DCB broadcast by PPP-B2b contain the pseudorange hardware delay of j/k frequency relative to B3I frequency. Unlike BDS, the broadcast ephemeris clock offset of GPS is estimated by the IFLC observations of L1 and L2, so the users do not need to correct DCB when using L1 and L2 dual frequency for positioning.

The real-time precise orbit and clock offset are brought into Equation (6), and then linearized to obtain the following observation equation:

$$\left\{\begin{array}{c}{\tilde{P}}_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s={\gamma }_r^{s}dx+\delta {\widehat{t}}_r+g_r^s{\tau }_r+{\epsilon }_P \\ {\tilde{L}}_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s={\gamma }_r^{s}dx+\delta {\widehat{t}}_r+g_r^s{\tau }_r+{\lambda }_{{\mathrm{I}\mathrm{F}}_{jk}}{B}_{r,{\mathrm{I}\mathrm{F}}_{jk}}^s+{\epsilon }_L\end{array}\right.$$

(7)

Then Kalman filter method is used for parameters estimation and Zenith Total Delay (ZTD) can be acquired.

2.5. Convert ZTD to PWV

ZTD is the sum of zenith hydrostatic delay (ZHD) and Zenith Wet Delay (ZWD). ZHD can be calculated by Saastamoinen empirical formula [30].

$$ZHD=10^{-3}{k}_1{R}_d\frac{P_s}{9.784[1000-2.66\mathit{cos}(2\lambda )-0.28H]}$$

(8)

Then, ZWD can be obtained by subtracting ZHD from the ZTD. It has been proved that the ZWD is related to PWV as follows [4],

$$\mathrm{P}\mathrm{W}\mathrm{V}=\prod ·\mathrm{Z}\mathrm{W}\mathrm{D}$$

(9)

Π is a mapping function that can be calculated by the following equation [4],

$$\prod =\frac{{10}^6}{{[\mathrm{R}}_{\mathrm{v}}(\frac{{\mathrm{k}}_3}{{\mathrm{T}}_{\mathrm{m}}}+{\mathrm{k}}_2^{\mathrm{\text{'}}})]}$$

(10)

where ${\mathrm{T}}_{\mathrm{m}}$ is the weighted mean temperature, ${\mathrm{T}}_{\mathrm{m}}=\frac{\int \left(\raisebox{1ex}{${\mathrm{P}}_{\mathrm{w}}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{T}$}\right.\right)\mathrm{d}\mathrm{z}}{\int \left(\raisebox{1ex}{${\mathrm{P}}_{\mathrm{w}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{T}}^2$}\right.\right)\mathrm{d}\mathrm{z}}$,${\mathrm{P}}_{\mathrm{w}}$ is the pressure of the water vapor, T is the Kelvin temperature, ${\mathrm{T}}_{\mathrm{m}}$ is often calculated empirically with the surface temperature Ts. ${\mathrm{T}}_{\mathrm{m}}=\mathrm{a}+\mathrm{b}·{\mathrm{T}}_{\mathrm{s}}$ [4].

2.6. Radiosonde and ERA 5 reanalysis data

There is an L-band radio-sounde station (#54511) in Beijing, which observes twice a day to obtain meteorological profiles of temperature, pressure, humidity, wind with a resolution of every minute. The observation errors are ±0.3℃ for temperature, ±0.5hPa for pressure, and ±5% for relative humidity.

ERA-5 is the latest reanalysis dataset from ECMWF's fifth generation [31], with finer spatial grids (79 km to 31 km), higher temporal resolution (3-6 hours to 1 hour), higher vertical resolution (60-137), new numerical weather prediction (NWP) models, and increased assimilation data volume compared to its predecessor ERA-Interim [32]. The dataset covers the period from 1950 to the present. As a three-dimensional atmospheric parameter reanalysis dataset providing global hourly resolution, ERA-5 serves as a good baseline for comparison and validation with other remote sensing observations. The precipitable water vapor can be calculated with the same integration method as for radiosonde.

2.7. The PWV accuracy evaluation

When it comes to assess the accuracy of PWV, three metrics, the average deviation (bias),the root mean square error (RMSE), and correlation coefficient (r) are often used for the validation, their formula are written as follows:

$$bias={\sum }_{i=1}^n\frac{\left(Y_i-Y_i^{\mathrm{\text{'}}}\right)}{n}$$

(11)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}={\left[{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Y}}_{\mathrm{i}}-{\mathrm{Y}}_{\mathrm{i}}^{\mathrm{\text{'}}}\right)}^2/\mathrm{n}\right]}^{\frac{1}{2}}$$

(12)

$$r={\sum }_{i=1}^n\left[\left(Y_i-{\overline{Y}}_i\right)\times \left(Y_i-{\overline{Y}}_i^{\mathrm{\text{'}}}\right)\right]/{\left[{\sum }_{i=1}^n\left[{\left(Y_i-{\overline{Y}}_i\right)}^2\times {\sum }_{i=1}^n{\left(Y_i-{\overline{Y}}_i^{\mathrm{\text{'}}}\right)}^2\right]\right]}^{\frac{1}{2}}$$

(13)

Where $Y_i,{\overline{Y}}_i$ are the values and the mean values of PWV, respectively, and n represents the sample number.

3. Results

3.1. Validation data for the GNSS PWV

Radiosonde are often used for comparison and validation of GNSS PWV [33]. However, radiosonde data as a reference also has the following disadvantages that the profile is not vertical at the specific site, as the sounding balloon flow due to the wind. What is more, it takes about 3 hours to finish a profiling, which makes it in poor temporal resolution and timeliness. Although it is often used to validate the GNSS PWV traditionly, the ERA5 is more appropriate in time and space matching.

Meteorological reanalysis uses advanced numerical prediction models and data assimilation systems to integrate model forecasts and historical observation data to obtain long-series historical weather data with rich variables, complete spatial coverage, and stable time uniformity. It enables meteorological agencies to possess more detailed and timely data in weather analysis and study. The ERA5 reanalysis is currently one of the best data for weather and climate research, though errors may exists in some region and some time [34,35].

Figure 2 shows the PWVs calculated from Radiosonde and ERA5 during the experiment. Between the radiosonde and ERA5, the bias of PWV is -0.27 mm, and the root mean square error is 0.08 mm. In terms of relative values, the relative bias is -6%, and the relative RMSE is 2%. This complies with the relative humidity error of the radiosonde, which indicates that ERA5 has a very good analysis quality and ERA5 should be a reliable data for PWV validation.

3.2. High tempo-resolutional PWV with the PPP solution based on B2b service

In China Meteorological Administration, the tradional network processing method is used by the GAMIT software. The zenith tropospheric delay and water vapor content are estimated hourly after a latency of 50 minutes with the sliding window (12 h) method [36].

Figure 3 shows the contrast of PWVs from the traditional network solution and the PPP solution with BDS on January 15, 2023. The temperal resolution of the traditional network solution is 1hour. On the contrary, with the PPP solution based on BeiDou's B2B service, PWV can be obtained every 2 minutes near real time. Very good agreement is shown for PWVs from the two solutions, but more details of the PWV variation are revealed from the PPP resolution.

3.3. Assessment of PWVs processing with single BeiDou and single GPS

Figure 4 shows that the PWVs processing with single BeiDou and single GPS contrast with the PWV calculation of ERA5. Since the ERA5 data are in 1 hour resolution, it changes smoother, where many peaks and great values are lost. Between PWVs from single BeiDou and ERA, there are a bias of 0.7 mm and a RMSEof 1.8 mm. While for the GPS case, the bias is 0.73 mm and RMSE is 1.97 mm. Both the BeiDou and the GPS processing performs very well in the estimation of PWV.

The ZTDs shows very similar situation as shown in Figure 5. Referenced with the ERA5 calculation, for the Beidou processing, the bias is 9.27 mm and the RMSE is 10.82 mm. For the GPS processing, the bias is 11.63 mm and the RMSE is 12.68 mm.

3.4. PWV from the BeiDou- GPS combined solution

In the PPP processing, we treat BeiDou MEO Satellites just as GPS and obtained the BeiDou -GPS combined solution. As shown in Figure 6, it is very similar to figure 4. Referenced with the ERA5 calculation, the bias is 0.74 mm and the RMSE is 1.97 mm for the BeiDou-GPS combined processing. Little improvement is found for the assistance of GPS to BeiDou and vice versa. This may be due to the poor handling of two different kinds of errors for the GPS and the BeiDou systems in our PPP software.

We optimize the PWV estimation after PWVs estimated from both GPS and BeiDou by the Bayesian model averaging [37]. Their weights is the inversions of their fitting errors. Table 1 is the biases and RMSEs for single GPS and BeiDou solutions, the combined solution and the fusion of PWVs.

4. Conclusions and Discussion

Three solutions are designed with PPP using the BeiDou system only, the GPS system only and the BeiDou-GPS combined system. Real-time PWVs for the 3 solutions are validated with the radiosonde and the ERA5 reanalysis data. The biases are less than 1 mm and RMSEs are less than 2mm. The BeiDou System performs very well, matching with the GPS system. But little improvement is found for the assistance of GPS to BeiDou and vice versa, which might be due to the poor handling of two different kinds of errors between the GPS and the BeiDou in our PPP software. Further improvement will be needed in the future. A better PWV estimation using the two systems is to estimate PWVs with single system only at first, and then to obtain the weighted mean estimation.

The PPP-based algorithm is much faster than the double difference network solution with very high temporal resolution for atmospheric water vapor estimation. Referenced with ERA5, the deviation and the RMSE are relatively small. Therefore, it has an advantage over the traditional network solution for severe weather analysis and warnings.

This testing experiment was conducted in Beijing in January. The water vapor content is relative small in this season. So it was a short experiment for the PWV validation. Long experiment will be needed in different seasons and different places for a comprehensive assessment in the future.