Retrieving Accurate Precipitable Water Vapor Based on GNSS Multi‐Antenna PPP With an Ocean‐Based Dynamic Experiment

As an attractive technique for measuring water vapor, the Global Navigation Satellite System (GNSS) faces additional challenges in dynamic applications such as in the open sea. We present a new method of retrieving precipitable water vapor (PWV) based on GNSS multi‐antenna precise point positioning (PPP), which uses GNSS data from multiple antennas and incorporates the constraints of known baseline vector and common tropospheric delay. The 4‐day shipborne dynamic experiment along the China coast demonstrates that the baseline vector constraint shortens the convergence time of positioning and atmospheric parameters, and also slightly improves their accuracies. The common tropospheric delay constraint helps to provide compromised, more robust, and sometimes more accurate PWV estimates. An evaluation with radiosonde‐derived PWVs shows that the combination of the two constraints achieves the best accuracy reaching 4.2 mm. This method helps to expand GNSS meteorology to the vast ocean and benefits satellite altimetry and weather forecasting.

10.1029/2023GL102982 2 of 9 instruments, the accuracy of GNSS-derived precipitable water vapor (PWV) has been extensively investigated, reaching few-millimeter level for land-based static GNSS stations (M. Wang et al., 2007;Yuan et al., 2014;Yu et al., 2017). The applicability of GNSS meteorology for weather forecasting/nowcasting has been proved by many studies (Gendt et al., 2004;Karabatić et al., 2011;Yan et al., 2009;Zhao et al., 2022). As a result, many land-based GNSS networks for meteorological applications have been established and rapidly developing, with notable examples of the American NOAA/FSL, the European E-GVAP, and the MAGIC for Mediterranean (Karabatić et al., 2011).
There are two main techniques associated with GNSS meteorology, namely relative positioning and precise point positioning (PPP). Relative positioning is widely adopted in GNSS meteorology studies, usually using well-developed software packages such as GAMIT and Bernese to produce network solutions (Benevides et al., 2015;Karabatić et al., 2011). However, relative positioning has the drawback that at least one static reference station is required to form one or more baselines. This makes relative positioning unsuitable for remote areas such as the open ocean. The PPP technique uses undifferenced observations from a single station for high-accuracy parameter estimation with precise satellite orbit/clock products and elaborate error modeling/ correction (Zumberge et al., 1997). Since no reference station is needed, it has the advantages of high efficiency and flexibility over relative positioning, and is quite suitable for remote areas. Many studies used PPP to obtain PWV and achieved comparable accuracy with relative positioning (Gendt et al., 2004;Byun & Bar-Sever, 2009;. PPP was also applied to obtain oceanic PWV on floating platforms such as buoys and ships, which is helpful to expand GNSS meteorology to the vast ocean and benefits satellite altimetry and weather forecasting (Bosser et al., 2021(Bosser et al., , 2022Fujita et al., 2008;Männel et al., 2021;Rocken et al., 2005;Wang et al., 2019). However, PPP generally requires a long convergence or reconvergence time, say, 20-30 min, before reaching cm-level positioning accuracy (Kouba & Héroux, 2001). Moreover, in dynamic environments such as in the open sea, it is faced with the additional challenge that coordinates and the tropospheric delay should be simultaneously estimated, which may potentially further degrade the performance, especially under rough sea conditions with severe wave and tidal actions.
To overcome the major drawbacks of PPP, some researchers introduced the constraint from multiple antennas mounted on the same platform to realize improved GNSS parameter estimation. Teunissen (2012) presented the concept of array-aided PPP which uses the known antenna-array geometry to improve the estimation of position, attitude, time and atmospheric delays. Array-aided PPP was illustrated to improve positioning accuracy, reduce convergence time, and achieve higher success rates of ambiguity resolution. Guan et al. (2021) proposed a baseline constrained PPP method, which applies the fixed baseline length constraint of two antennas to the Kalman filter, to improve the positioning accuracy of PPP.
Inspired by previous studies, we present a new method for retrieving PWV based on GNSS multi-antenna PPP. This method exploits the increased redundance provided by multiple closely spaced antennas, which is expected to shorten the convergence time and improve the accuracy of GNSS-derived PWV. The effectiveness of this method is demonstrated in ocean-based dynamic environments with a 4-day shipborne experiment.

Multi-Antenna PPP Method for PWV Retrieval
The multi-antenna PPP uses GNSS data from multiple antennas closely spaced and mounted on the same platform. Compared to traditional PPP, multi-antenna PPP can incorporate additional constraints of known baseline vector and common tropospheric delay for closely spaced antennas. By incorporating these constraints, multi-antenna PPP is expected to shorten the convergence time and improve the accuracy of estimated parameters including positioning coordinates and zenith wet delay (ZWD) (Chen et al., 2020;Guan et al., 2021;Teunissen, 2012).
Assuming a two-antenna configuration, which can be readily generalized to three-or more-antenna cases, and bearing in mind the ionosphere-free (IF) observation model of traditional PPP, the observation model for multi-antenna PPP can be described as follows: 10.1029/2023GL102982 where IF is the IF combination of GNSS carrier phase observations after proper error corrections, and IF is the corresponding wavelength; the subscript ( = 1, 2) stands for the th antenna or station; is the geometric distance from the receiver to the satellite, with ( , , ) and ( , , ) representing the receiver and satellite coordinates, respectively; is the receiver clock error and is the light speed; IF is the ambiguity of the IF combination; and IF represents the sum of residual errors.
The baseline vector constraint can be expressed as: where [Δ 12 Δ 12 Δ 12] moving−baseline is the baseline vector coordinates from station 1 to 2, which can be precisely determined using GNSS relative positioning with the moving-baseline mode (Takasu, 2013) in a real-time manner. Previous studies use nonlinear constraints of known baseline length or direction, where linearization errors may lead to estimation performance degradation or even divergence, especially for very short baselines (Chen et al., 2020;Guan et al., 2021). By contrast, the baseline vector constraint of Equation 2 presented in this study is more advantageous because it contains full baseline information and is free from any linearization errors. This baseline vector constraint can be used as pseudo-observations in the parameter estimation with a given weight , similar to the practice of Ge et al. (2006).
The tropospheric delay in the observations of Equation 1 can be corrected by first calculating the zenith tropospheric delay using an empirical model (see, e.g., Zhao et al., 2023), and then mapping it to the satellite elevation angle. The zenith tropospheric delay can be separated into the zenith hydrostatic delay and the ZWD, the latter of which is related to the water vapor in the troposphere. Considering that ZWD is highly variable and difficult to model, an additional parameter, namely the residual ZWD (RZWD), is usually required in PPP. In traditional PPP, RZWD is independently parameterized for each station, for example, RZWD 1 for station 1 and RZWD 2 for station 2, as demonstrated in Equation 1 with being the mapping function. In our multi-antenna PPP, however, we take into account the correlation of RZWDs, which can be regarded as equivalent for closely spaced stations. Therefore, we set a common parameter RZWD for different stations expressed as By adding the estimate of RZWD to the a priori ZWD, we can obtain an accurate estimate of ZWD, which can be converted to PWV using the linear regression of Bevis et al. (1994). We have implemented the multi-antenna PPP method for PWV retrieval on the basis of the PANDA software developed by Wuhan University (Liu & Ge, 2003). The evaluation of this method with a 4-day shipborne experiment will be described in the following sections.

Shipborne Experiment and Data Processing
We conducted a 4-day shipborne experiment on 24-27 July 2019 along the China coast, with the ship trajectory mainly shown in Figure 1. Two GNSS antennas with the model of TRM159900.00 were mounted on the 2000ton ship "XiangYangHong-18." The two antennas, denoted by antennas 1 and 2, were located on the left and right sides of the ship, respectively. Each antenna was connected to a Trimble Alloy GNSS receiver. During the experiment, multi-constellation and multi-frequency GNSS data were collected at a sampling rate of 100 Hz. The satellite numbers and Position Dilution of Precision (PDOP) values considering GPS only are shown in Figure S1 in Supporting Information S1. The mean satellite numbers are 7, 8, 8 and 8, and the mean PDOP values are 2.6, 2.5, 2.2 and 2.4 on 24-27 July, respectively.
We processed the GNSS data with the modified PANDA software in kinematic PPP mode, using the IF combination of dual-frequency GPS observations. The processing time interval was 30 s and the elevation cut-off angle was set as 10°. The a priori tropospheric delay was calculated using the Saastamoinen model (Saastamoinen, 1973) with the global mapping function (Boehm et al., 2006). The RZWD and tropospheric gradients were parameterized as piecewise constants with interval length of 1 and 6 hr, respectively. The station coordinates and receiver clock errors were estimated on an epoch-by-epoch basis. A single float ambiguity was estimated for each continuously tracked orbital arc of any satellite. Precise satellite orbit and clock products were obtained from the International GNSS Service (IGS) with time intervals of 15 min and 30 s, respectively.
To investigate the effects of the two additional constraints in our multi-antenna PPP method on the derived PWV, we have processed the GNSS data with four strategies. The four strategies include the Conventional one using conventional PPP without constraints, the Bl one with baseline vector constraint of Equation 2, the Com one with common ZTD constraint of Equation 3, and the Bl + Com one with both baseline vector and common ZTD constraints of Equations 2 and 3. The reference values for the GNSS-derived PWVs come from the ERA5 hourly reanalysis data of ECMWF with a spatial resolution of 0.25°. They are interpolated at the location of the ship and serve the purpose of our evaluation. We should note that as a reanalysis product, the spatial resolution of ERA5 has a different implication from that of in-situ observations. The topography in the reanalysis model is rougher than the spatial resolution of output values, generally introducing additional errors.

Results and Analysis
We will first investigate the effect of the baseline vector constraint by comparing the Conventional and Bl solutions. During the dynamic experiment, although the coordinates and relative orientation of the two receivers changed over time, the baseline length kept fixed all the time. We compute the instantaneous baseline length from coordinate estimates of the two receivers, and show its time series in the left column of Figure 2 for the two solutions, together with the reference value of 12.314 m.
At the beginning of each day, the baseline length of the conventional solution exhibits a relatively long convergence process. After the first convergence, there exist different noticeable behaviors with sharp jumps and quickly damped variations, except for 25 July. These noticeable behaviors indicate several re-convergence processes probably caused by temporary loss of lock in satellite tracking or detected cycle slips. The convergence periods are shown in gray bars for further analysis. As shown in Figure 2, by adding the baseline vector constraint, the Bl solution can greatly shorten the convergence time, from generally 1 hr or longer to almost instantaneous convergence. After convergence, baseline lengths of both solutions are very close to the reference value. As can be seen roughly from the fluctuations and widths of the lines in Figure 2, The Bl solution agrees with the reference value better with a lower noise level. This comparison proves that the baseline vector constraint works well and enhances the positioning accuracy in dynamic conditions.
The ZWD estimates within the first 2.5 hr at stations 1 and 2 are given in the middle and right columns of Figure 2, respectively. The forward filtered results, represented in blue or red and denoted as rt, are regarded as simulated real-time estimates. The backward smoothed results, represented in black and denoted as smt, are regarded as post-processed estimates and serve as reference values. By comparing the convergence speed of the real-time estimates, we can roughly say that the Bl solution achieves similar or better performance. We compute the convergence time and list them in Table S1 in Supporting Information S1. According to Table S1 in Supporting Information S1, the convergence time of the Bl solution is shortened under challenging environments with original convergence time exceeding 1 hr, but could be prolonged under friendly environments with original convergence time within 1 hr. An average improvement of 11% in ZWD convergence speed is still achieved with the baseline vector constraint.
We will then evaluate the performance of GNSS PWV derived with conventional and multi-antenna PPP. The time series of PWV from the Conventional, Bl, Com and Bl + Com solutions, together with the ERA5 PWV as the reference, are shown in Figure 3 for 24 July and Figures S2-S4 in Supporting Information S1 for 25-27 July.
We also compute the differences between the GNSS and ERA5 PWVs and show them in Figure S5 in Supporting Information S1. The curves of ERA5 PWV exhibit a smooth shape with different trends for the 4 days, with ranges of 51. 2-57.4, 54.4-66.8, 62.5-66.4, and 66.3-71.8 mm on 24-27 July, respectively. At the beginning and a few follow-up periods, the GNSS PWV values exhibit large jumps and obviously deviate from the reference, especially for the Conventional solution. These periods correspond to the convergence and re-convergence processes for ZWD parameters, parallel to the convergence periods of baseline length labeled by gray bars.
As shown by panels a and b in Figure 3 and Figures S2-S4 in Supporting Information S1, the Bl solution generally reduces these large jumps compared to the Conventional solution. This proves that the baseline vector constraint can accelerate the convergence process for ZWD parameters, which is consistent with above results regarding the ZWD time series in Figure 2 and Table S1 in Supporting Information S1. Panel c in Figure 3 and Figures S2-S4 in Supporting Information S1 shows an interesting pattern for the Com solution, whose values are almost always between those at stations 1 and 2. Rather than a crude average, this is a compromise in a statistical sense, for example, at 0 and 10 hr in Figure 3c. This implies that the common ZTD constraint helps to obtain a more robust estimates than a single-antenna solution. By combining the two constraints together, the Bl + Com solution shows a very similar pattern to the Com Solution. However, during convergence periods, the Bl + Com solution has smaller jumps compared to the Com solution, which is attributed to the baseline vector constraint. This phenomenon is clearer by comparing these two solutions in Figure S5 in Supporting Information S1. After convergence (namely the periods without gray bars), the four GNSS PWV solutions are close to each other and well consistent with ERA5 PWVs. The differences between the GNSS and ERA5 PWVs fluctuate around zero with magnitudes of only a few mm.
To quantitatively analyze the performance of GNSS-derived PWV, we compute the root mean square errors (RMSEs) of the differences between GNSS and ERA5 PWVs and list them in Table 1. To exclude the interference of convergence periods, we use the backward smoothed GNSS PWVs when computing RMSEs. For the Conventional solution, the RMSEs during the 4 days range within 2.15-5.66 and 2.46-5.14 mm at stations 1 and 2, respectively. By inspecting the Bl solution, we find that the baseline vector constraint can generally slightly reduce the RMSEs for both stations. Exceptions occur on 26 and 27 July at station 1, where the Bl solution has negligibly larger RMSEs. The average improvements of the Bl solution are 2% and 6% for stations 1 and 2, respectively. Unlike the Bl solution, the Com solution provides RMSEs either in between or smaller than those of the Conventional solution at the two stations, numerically demonstrating that the common ZTD constraint helps to provide compromised and more robust PWV estimates. The RMSEs of the Bl + Com solution are almost the same as those of the Com solution, indicating that the common ZTD constraint predominates the results when convergence is already achieved.
To evaluate the GNSS PWVs more objectively, we have collected the radiosonde-derived PWVs from the Integrated Global Radiosonde Archive (IGRA). The IGRA data set involves more than 1,500 global distributed stations, typically observing at 00UTC and 12UTC (Durre et al., 2006). The  GNSS PWVs are compared with the radiosonde-derived PWVs of nearest IGRA stations at common epochs, with an average separation of 230.72 km between GNSS and IGRA stations. The distribution of involved IGRA stations is shown in Figure 1. The GNSS PWVs after backward smoothing during the 4-day shipborne experiment are shown in Figure 4, with the radiosonde-derived PWVs shown as a reference. The RMSEs of the differences between the GNSS-and radiosonde-derived PWVs are listed in Table S2 in Supporting Information S1. According to Figure 4 and Table S2 in Supporting Information S1, compared to the Conventional solution, the constrained solutions achieve slightly higher accuracies. The Bl + Com solution which combines the two constraints achieves the best accuracy reaching 4.2 mm. Such a few-millimeter-level accuracy is comparable with the reported performances for ground-based static GNSS, which well meets the requirements of meteorological applications over coastal or oceanic areas (Bosser et al., 2021;Rocken et al., 2005).
From the above analysis, we become aware that the two constraints in multi-antenna PPP contribute to PWV retrieval in somewhat different ways. The baseline vector constraint mainly helps to shorten the convergence time of positioning and atmospheric parameters, and also slightly improves their accuracies. The improved convergence performance makes multi-antenna PPP very valuable for near real-time meteorological applications, which increases the availability of PWV results. This is particularly true in ocean-based dynamic environments, where reconvergences frequently occur due to the disturbance of wave and tidal actions and severe multipath effects. Unlike the baseline vector constraint, the common ZTD constraint helps to provide compromised, more robust, and sometimes more accurate PWV estimates. This is also valuable because when one station encounters bad observation or failure, this station is automatically degraded and the other station(s) are fully used to still provide accurate PWV results. The combination of the two constraints suppresses PWV jumps during convergence periods, and achieves a negligibly higher accuracy compared to the common ZTD constraint after convergence. This shows the advantage of combining the two constraints over any single-constraint solution.

Conclusions
This study presents a new method of retrieving PWV based on GNSS multi-antenna PPP. Through a 4-day shipborne dynamic experiment, the baseline vector constraint is shown to shorten the convergence time of positioning and atmospheric parameters, and also slightly improves their accuracies. The common tropospheric delay constraint helps to provide compromised, more robust, and sometimes more accurate PWV estimates. An evaluation with radiosonde-derived PWVs shows that the combination of the two constraints achieves the best accuracy reaching 4.2 mm. This method helps to expand GNSS meteorology to challenging environments such as in the open sea, which will benefit satellite altimetry, and weather forecasting and nowcasting.