Real-time data communication
A 10M + 2M dual-link fiber network is established to realize the communication between the GNSS/MET operational network and the national integrated data processing center of BGAN (as shown in Fig. 3). And OSPF (Open Shortest Path First) protocol is employed with defining different verification methods for areas, and exchanging information through interactive verification of routes, improving the security of the data communication process effectively. 175 reference stations were built for the BGAN, of which CMA is responsible for building 35 stations. The real-time data communication has enhanced and integrated the existing ground-based GNSS/MET network, which is conducive to building an intensive application system, establishing unified data application standards, and improving data application efficiency.
Data processing for meteorological service
The ground-based GNSS/MET uses navigation satellite data to calculate the ZTD (Total Tropospheric Delay), and then separates ZHD (Zenith Hydrostatics Delay). The difference between ZTD and ZHD is the ZWD (Zenith Wet Delay), and then PWV can be derived with the conversion function. As shown in Fig. 3, data processing includes runtime monitoring and meteorological product inversion, and the functions are deployed in the real-time monitoring area and core operation area respectively. The real-time monitoring area realizes functions such as data reception monitoring, product monitoring, system monitoring, and fault alarm. The core operation area achieved real-time processing of data from 175 BGAN stations nationwide and about 1000 GNSS/MET operational network stations of CMA, and then retrieves water vapor products with time resolution of 1-hour, meeting the needs of water vapor monitoring and numerical weather prediction.
The key technology of data processing is the solution to combining multi-GNSS of BDS and GPS. In the water vapor solution of multi-GNSS data, BDS and GPS share a common tropospheric zenith delay parameter. Compared to single GNSS system, combination processing can theoretically improve the accuracy and reliability of calculation by adding redundant observations. However, compared to single satellite system data processing, BDS and GPS data combination must consider the impact of different hardware delay biases between different satellite systems (Lu et al. 2016), and usually all computed biases of other system are obtained relative to the biases for the GPS observations (Li et al. 2015).
Solution to combining multi-GNSS of BDS and GPS
The commonly used GNSS/ZTD solutions include PPP (Precise Point Positioning) and DD (the network Double Difference) modes (Li et al. 2015; El-Mowafy et al. 2016). The PPP mode utilizes precise satellite orbits and clock products to directly obtain ZTD, which is superior to the DD mode in terms of observation requirements and computational efficiency. For the real-time PPP processing, precise satellite orbits and clocks have to be first determined using the observation data from ground-based GNSS network. With the development of GNSS, PPP mode has become an important development trend in the application of combination of multi-GNSS. A lot of GNSS meteorological applications have confirmed that the two processing modes provide ZTD results with similar quality (Xu et al, 2013).
Generally, in order to eliminate the first-order effects of the ionosphere, the ionosphere-free combination of dual-frequency carrier phase and pseudo-range are used in PPP processing when multiple GNSS systems are combined and applied (Lu et al. 2016). And precise ephemeris can also eliminate the satellite code delay biases in ionosphere-free combination model. (Li et al., 2015; Montenbruck et al. 2014; Jiang et al. 2015; Wanninger 2011).
The inter system bias (ISB) is introduced for signal code bias correction with the GPS receiver clock used as the reference. With satellite orbit, satellite clock, station position, and model correction such phase center offset, solid tide, phase wind-up, relativistic effects, etc. (Kouba 2009), the combined observation equation can be expressed as follows.
$$\left\{\begin{array}{c}{P}_{r}^{C}\left(i\right)={\rho }_{r}^{C}\left(i\right)+c\cdot d{t}_{r}\left(i\right)+c\cdot {d}_{rC}+{m}_{r}^{C}\left(i\right)\cdot zt{d}_{r}\left(i\right)+{\epsilon }_{r,P}^{C}\left(i\right) \\ {P}_{r}^{G}\left(i\right)={\rho }_{r}^{G}\left(i\right)+c\cdot d{t}_{r}\left(i\right)+c\cdot {d}_{rG}+{m}_{r}^{G}\left(i\right)\cdot zt{d}_{r}\left(i\right)+{\epsilon }_{r,P}^{G}\left(i\right)\\ {\varphi }_{r}^{C}\left(i\right)={\rho }_{r}^{C}\left(i\right)+c\cdot d{t}_{r}\left(i\right)+{\lambda }_{C}({b}_{rC}-{\text{b}}^{C})+{m}_{r}^{C}\left(i\right)\cdot zt{d}_{r}\left(i\right)-{\lambda }_{C}{N}_{r}^{C}+{\epsilon }_{r,\varphi }^{C}\left(i\right)\\ {\varphi }_{r}^{G}\left(i\right)={\rho }_{r}^{G}\left(i\right)+c\cdot d{t}_{r}\left(i\right)+{\lambda }_{G}({b}_{rG}-{\text{b}}^{G})+{m}_{r}^{G}\left(i\right)\cdot zt{d}_{r}\left(i\right)-{\lambda }_{G}{N}_{r}^{G}+{\epsilon }_{r,\varphi }^{G}\left(i\right)\end{array}\right. \left(1\right)$$
where \(i\) is the epoch number, C and G represent BDS and GPS satellites respectively, and r is the station number;\(P\) is the pseudo-range, \(\varphi\) is the phase observation༛\(\rho\) is the geometric distance between the measuring station and the satellite༛\(d{t}_{r}\) represents the clock deviation of station, \(c\) i s the speed of light in vacuum༛\(m\) indicates the global mapping function of the ZTD (Böhm et al. 2006); \(\lambda\) is the wavelength of the observation phase observed, \(N\) is the phase ambiguity; \({d}_{rC}\) and \({d}_{rG}\) are the code biases of the receivers receiving BDS and GPS signals respectively, while \({b}_{rG}\) and \({b}_{rG}\) are the phase delay biases of the receivers; \({b}^{C}\)and \({b}^{G}\) are the phase delay biases of the BDS and GPS satellites respectively༛\(\epsilon\) represents other errors. To simplify the formula, we merge them, and the parameters to be estimated in the equation are followed.
\(d{t}_{k}\) : the clock deviation of station, one parameter to be estimated per epoch per station.
\({N}_{k}^{j}\) : the phase ambiguity, number of the parameters to be estimated depends on the number of satellites observed and the status of ambiguity.
\(zt{d}_{k}\) : ZTD, one parameter to be estimated for a period of data processed (such as 15 minutes, half an hour, one hour, etc.).
The phase delay biases can be merged into the corresponding ambiguity parameter, and the estimated signal code bias between BDS and GPS is ISB. And the formula (1) can be expressed as following.
$$\left\{\begin{array}{c}{P}_{r}^{C}\left(i\right)={\rho }_{r}^{C}\left(i\right)+c\cdot d{t}_{r}\left(i\right)+c\cdot IS{B}_{r}+{m}_{r}^{C}\left(i\right)\cdot zt{d}_{r}\left(i\right)+{\epsilon }_{r,P}^{C}\left(i\right)\\ {P}_{r}^{G}\left(i\right)={\rho }_{r}^{G}\left(i\right)+c\cdot d{t}_{r}\left(i\right)+{m}_{r}^{G}\left(i\right)\cdot zt{d}_{r}\left(i\right)+{\epsilon }_{r,P}^{G}\left(i\right)\\ {\varphi }_{r}^{C}\left(i\right)={\rho }_{r}^{C}\left(i\right)+c\cdot d{t}_{r}\left(i\right)-{\lambda }_{C}{\stackrel{\sim}{N}}_{r}^{C}+c\cdot IS{B}_{r}+{m}_{r}^{C}\left(i\right)\cdot zt{d}_{r}\left(i\right)+{\epsilon }_{r,\varphi }^{C}\left(i\right)\\ {\varphi }_{r}^{G}\left(i\right)={\rho }_{r}^{G}\left(i\right)+c\cdot d{t}_{r}\left(i\right)-{\lambda }_{G}{\stackrel{\sim}{N}}_{r}^{G}+{m}_{r}^{G}\left(i\right)\cdot zt{d}_{r}\left(i\right)+{\epsilon }_{r,\varphi }^{G}\left(i\right)\end{array}\right. \left(2\right)$$
where \({\stackrel{\sim}{N}}_{r}^{C}={N}_{r}^{C}+{\text{b}}^{C}-{b}_{rC}\), \({\stackrel{\sim}{N}}_{r}^{G}={N}_{r}^{G}+{b}^{G}-{b}_{rG}\). In the formula, ISB is generally estimated over a period (such as a day or a period of data processed). It is found that ISB is extremely stable with time, and it is the normally distributed approximately.
PWV inversion
\(ZWD=ZTD-ZHD\) , and PWV can be converted to when ZWD is accurately estimated.
$$\text{P}\text{W}\text{V}=\prod \bullet ZWD \left(3\right)$$
$$\prod ={10}^{6}\times {\left[{{\rho }}_{\text{w}\text{a}\text{t}\text{e}\text{r}}{\text{R}}_{\text{w}}\left({\text{k}}_{2}^{{\prime }}+\frac{{k}_{3}}{{\text{T}}_{\text{m}}}\right)\right]}^{-1} \left(4\right)$$
where \({k}_{1},{ k}_{2},{ k}_{3}\) is the physical constant, \({\text{k}}_{2}^{{\prime }}={k}_{2}-{k}_{1}({R}_{d}/{R}_{w})\);\({R}_{d}\) and \({R}_{w}\) are gas constant of dry air and the water vapor respectively. \({\text{T}}_{\text{m}}\) can be estimated based on the ground temperature \({\text{T}}_{s}\), and usually \({\text{T}}_{\text{m}}=\text{a}+\text{b}{\text{T}}_{s}\). (Davis et al.1985; Bevis et al. 1992)
ZTD estimation is an important part of ground-based GNSS/PWV inversion. To ensure the integrity and continuity of PWV, prior models such as GPT or GPT2 models are used in ZTD estimation, and a prior ZHD is calculated with air pressure from the prior models. ZHD can be obtained more accurately (better than 1mm) through models such as Saastamoinen model with the precise air pressure observation (Saastamoinen 1972), and a pressure error of 1 hPa can cause ZHD error of nearly 3 mm. Altitude difference between meteorological equipment and GNSS antenna is considered in product inversing, which usually causes deviations in meteorological parameters.