Retrieval of Vertical Distribution of Tropospheric Refractivity through Ground-Based GPS Observation

Atmospheric refractivity is an important index for atmospheric stratification, which is mainly affected by pressure, temperature, and humidity. Variations in refractivity with respect to height can cause radio waves to bend or travel abnormally. These variations can also affect microwave communication and atmospheric detection. Therefore, atmospheric refractivity should be profiled to predict the propagation path of radio waves and the range of radar detection [Saastamoinen, 1972]. In addition to a refractive effect of the vertical structure of atmospheric refractivity on electromagnetic waves, variations in water vapor and temperature stratification are indirectly associated with this vertical structure that may indicate the evolution of a convective weather system and contribute to the improvement of short-term weather prediction (Roberts et al., 2008).

Since the proof-of-concept GPS Meteorology (GPS/ MET) experiment was launched in 1995 (Ware et al., 1996; Kursinski et al., 1997), the atmospheric limb sounding technique that uses radio signals transmitted by GPS satellites has become an emerging and promising approach in atmosphere remote sensing. In contrast to space-borne radio occultation techniques, ground-based GPS observation can focus on an individual region with a higher temporal resolution. Studies have been conducted on the capability of ground-based GPS observation to detect precipitable water vapor (Bevis et al., 1992; Rocken et al., 1993; Ware et al., 1997; Li et al., 1999; Wolfe and Gutman, 2000; Foelsche and Kirchengast, 2001; Wang and L¨u, 2005), and such systems have been assimilated in numerical weather forecast models (Rocken et al., 2003; Cucurull et al., 2004; Troller et al., 2006; Macpherson et al., 2008). However, ground-based GPS observation is rarely used in atmospheric profiling because receivers on platforms near the horizon fail to receive signals from negative elevations, producing ill-posed problems in Abel inversion [Lowry et al., 2002]. Therefore, other retrieval methods such as atmospheric tomography supported by a dense network of GPS receivers or low-elevation angle observations via an individual GPS receiver should be used in ground-based GPS observation to obtain the profile information of meteorological parameters [Rocken et al., 2003]. The tomography technique can also monitor temporal and spatial changes in the atmosphere in a local area (Flores et al., 2000; MacDonald and Xie, 2000). However, this technique requires a number of receivers to work and is relatively inflexible in terms of observation location. In contrast to the tomography technique, single-site ground-based profiling requires only one receiver and is more flexible in terms of the selection of observation platforms, e.g., building tops, vehicles, ships, and other sites where observations are made at low elevation angles.

Studies have been conducted to explore the capability of the ground-based GPS receiver to profile atmospheric refractivity. For instance, [Lowry et al., 2002] proposed a preliminary three-level vertical model to characterize the ducting conditions near the coast and stated that an excess phase path of GPS signals, i.e., the phase delayed by neutral atmosphere, in the case of ground-based observation contains more information than that of other observables, e.g., the Doppler shift. The authors further highlighted the ability of such a technique to detect the onset altitude of ducting from 1 to 2 km by setting a fixed vertical refractivity gradient of -160 km^-1 N-unit km^-1 (the critically refractive gradient). However, the performance of this method when used in other areas remains inadequately described. For instance, the use of a fixed refractive gradient inland is dangerous, considering that assumptions near the coast are seldom possible. Another study used a similar three-level model to perform a theoretical retrieval of the atmospheric refractivity near the boundary layer from a ground-based GPS bending angle and an atmospheric delay, but no practical observations were included. Single-site ground-based GPS observations, as a convenient and low-cost observation technique, should be further studied. The practical use of this technique under more common conditions and at a larger altitude range is also interesting.

In the work reported in the present paper, a three-level approach based on [Lowry et al., 2002] was considered. We present an improved retrieval approach that extends the retrieval altitude from the boundary layer to 10 km and removes the assumed refractive gradients. The overall aim of the study was to investigate the capability and limitations of single-site ground-based GPS measurements to perform profiling of tropospheric refractivity in a common background. To verify the effectiveness of the method, a theoretical retrieval was initially conducted. Practical experiments were then performed, in which neutral atmospheric delays obtained from the measurements of a dual-frequency ground-based GPS receiver were used to derive the vertical distribution of refractivity. The inversion results were compared with nearby radiosonde measurements, and a statistical analysis was conducted to estimate the practical performance of the retrieval technique.

The remainder of the paper is organized as follows. Section 2 presents a brief introduction to the ray-tracing model, which was used to simulate the atmospheric delays from refractivity profiles. Section 3 describes the retrieval technique and the theoretical retrieval experiment to verify this technique. Section 4 presents the practical retrieval. Section 5 discusses the inversion results, as well as the potential and limitations of the practical retrieval. Section 6 summarizes the key findings of the study.

This section provides a very brief description of the ray-tracing model that we used in the retrieval process.

GPS satellite signals are refracted as these signals travel through the atmosphere because of the refractive gradients along their paths. Ray-tracing is a simplified geometrical model, in which the atmospheric delays are defined as the difference between the actual path of the GPS signal and the straight-line path that the signal could possibly take in the absence of the atmosphere. The ray-tracing model can be used to simulate the delays induced by the neutral atmosphere when the atmosphere is a refractive medium assumed as spherical and radially symmetric (Eresmaa et al., 2008; Nievinski and Felipe, 2009).

Considering local spherical symmetry, the signal propagation path in the neutral atmosphere of any azimuth is expressed as an integral expression along ray path [Davis et al., 1995]:

where and n is the refractive index. In Fig. 1, R₁ is the distance between the local center of the curvature and the ground-based GPS receiver; R₂ is the distance between the local center of the curvature and the upper bound of the neutral atmosphere; l is the real path; r is the distance from the center of the curvature (in practice, the center of the spheroid is defined as the local center of the curvature); dr=dl×cosα;α is the zenith angle; β is the angle of elevation; and θ is the angle between R₁ and R₂.

The signal path in a vacuum is expressed as follows:

Thus, the simulated atmospheric delay ΔS for a given refractivity profile n(r) is expressed as follows:

ΔS=S-S’ (3)

Figure 1.  Geometry of GPS signal propagation.

At an altitude of ≥40 km, the neutral atmosphere barely affects the GPS signal path. Therefore, we assumed that 40 km is the upper bound of the neutral atmosphere in the following retrieval.

A more detailed derivation process of ΔS as well as calculation of β and α can be found in Lowry’s paper [Lowry et al., 2002]. It is important to note that the angle of elevation βand the complementary angle of the zenith angle (π/2-α) are different, and the gap between β and (π/2-α) increases as βdecreases. At a certain β, a (and the only) corresponding α is defined if the refractivity profile is specified.

Figure 2.  Comparison of (a) COSMIC and CIRA refractivity profiles, for which the COSMIC refractivity profile was a measurement on 21 March 2011 in the area of (40.0°N, 122.0°E) to (41.4°N, 121.9°E); and (b) neutral atmosphere delays calculated with COSMIC and CIRA refractivity profiles.

Figure 3.  (a) Search space for the “theoretical retrieval”; (b) inversion result compared with initial (CIRA) profile and “true” (COSMIC) profile; (c) relative error of the initial profile and the inversion result compared with the “true” profile.

A previous study showed that refractivity is retrieved on a three-level model (from the surface to approximately 2 km) with fixed vertical gradients between two adjacent levels. The temporal and spatial variations in temperature, pressure, and water vapor in the lower and middle troposphere (below approximately 10 km in tropical and mid-latitude regions) lead to a variable vertical structure of refractivity, which cannot be precisely characterized by climatology data. Therefore, we extended the vertical range to 10 km. We also defined 11 constant levels linearly from the surface to 10 km instead of defining the fixed vertical refractive gradient, which is not always necessarily true. We determined the refractivity values on each of the levels that eventually constitute a “best-fit” refractivity profile.

The “exhaustive search” method is implemented on each level to determine the best solution, which is a method used to ensure that a solution to a discrete problem is provided. This procedure also ensures that the search process covers all of the possible values at any level and generates an overall optimal solution. However, the exhaustive search may require substantial computing resources given that the number of potential solutions is exponential to the number of levels [Nievergelt, 2000]. To expedite the exhaustive search, the search space was reduced. To reduce the search space, we used five-year data of the COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate) mission [Rocken et al., 2000] from 2007 to 2011 to estimate the reasonable range of the initial profile in advance.

The search begins from the initial profile provided by the CIRA86aQ_UoG model [Kirchengast et al., 1999]. Variable search steps are employed instead of using a fixed search step, e.g., 5 N-unit on each level. Given that refractivity decreases exponentially with height, a fixed search step suitable for refractivity at a lower altitude probably results in a large relative error at a high altitude. The atmospheric delay is an integral effect of the refractive atmosphere along the ray path. Refractivity and refractive gradients are larger in the lower troposphere than those in the upper troposphere. The variations in these two variables indicate that refractive effects at a lower altitude contribute more to overall ΔS. Thus, retrieval accuracy should be guaranteed by using more precise search steps at lower levels. The search steps on levels 1-7 (0-6 km) are set at 1% of the initial refractivity on that level, and the search steps on levels 8-11 (7-10 km) are set at 2%.

The cost function is introduced to determine the “best-fit” solution. The cost function is defined as the sum of squares of the residuals:

Where ΔS is the observed delay and ΔS’ is the delay simulated by the ray-tracing model from every refractivity profile during the exhaustive search. A corresponding α is determined when the observation experiment is performed for each of the observations. num is the total number of α and the number of observations used for retrieval. K=1,2,3,…, num. The refractivity profile that makes the minimum F is considered as the “best-fit” solution.

In this section, we report a theoretical retrieval that was conducted using the methodology described in section 3.1. The refractivity profile provided by the COSMIC measurement on 20 March 2011 at (40.0°N, 122.0°E) and (41.4°N, 121.9°E) (COSMIC refractivity profile) was considered as the “true” profile. The delays simulated by the ray-tracing from the “true” profile were used to constrain the retrieval [as ΔS in Eq. (4)]. The initial refractivity profile was calculated from the zonal average pressure and the temperature of the CIRA86aQ_UoG data in March by using the following expression:

where P is the air pressure in hPa, T is the temperature in K, and e is the water vapor pressure in hPa.

Figure 2a shows the “true” profile and the initial profile. Figure 2b shows the delays simulated from the two respective profiles based on the ray-tracing model as a function of elevation. As can be seen, the atmospheric delays increased sharply as the elevation angles decreased, and the most significant differences were observed in the signals received at low elevation angles. Thus, only those measurements were used in the retrieval we now describe. We used five years of COSMIC data from March 2007 to 2011 at (35°-45°N, 113°-127°E) to constrain the exhaustive search space (Fig. 3a). α was given from 84.50° to 89.50° for every 0.05°, and the corresponding β was from 0.04° to 5.36°.

The results of the theoretical retrieval are shown in Fig. 3b. The relative error is shown in Fig. 3c. Figure 3 also shows that the retrieved refractivity profile was similar to the COSMIC refractivity profile.

We performed 100 cases to evaluate the retrieval error statistically and verify the accuracy of the retrieval. Figure 4 shows the statistical relative error of these 100 cases. The relative errors ranged from -2% and 2%, and the range of errors did not change much with altitude. In these cases, the “true” refractivity profile was determined in advance; thus, no observation error or error induced by mismatch was found between observed and simulated delays. Therefore, this retrieval error was mainly attributed to the retrieval method, suggesting that this retrieval method is adequate in practical retrieval.

Figure 4.  Statistical relative error of the “theoretical retrieval”. Triangles indicate the average of relative error on each level and the black line and dots indicate the error range.

An observation experiment was performed in May 2011 to collect GPS measurements by using a dual-frequency GPS receiver (Leica RS500) and a choke ring antenna (Leica T504). These instruments were placed on top of the Institute of Atmospheric Physics, Chinese Academy of Sciences (39.98°N, 116.38°E; altitude = 83.91 m). GPS signals can be received at that location from almost all directions without evident hindrances. The GPS data were sampled every 30 s (1/30 Hz). The observation data (O-file) collected by the GPS receiver were in RINEX format. The O-file data were the pseudo-ranges of GPS signal frequencies (L1 = 1575.42 MHz, L2 = 1227.60 MHz) and carrier phases of the C/A code (coarse/acquisition code) on L1 and the P code (precise code) on L2.

Precise satellite orbits and clock data were downloaded from the Scripps Orbit and Permanent Array Center (N-file). These data were then interpolated at a ratio of 1/30 Hz. The surface temperature and humidity were collected using a meteorology sensor near the antenna. The pressure data were provided by the National Meteorological Information Center, China Meteorological Administration. The refractivity at the receiver (the first level in the practical retrieval) was calculated with these surface meteorology data by using Eq. (4). In the following retrieval process, the refractivity at the first level of the 11-level model was fixed on the refractivity at the receiver. These meteorology data were interpolated at a ratio of 1/30 Hz and written in a RINEX format (M-file).

The atmospheric delays were estimated using GAMIT/ GLOBK version 10.4. We used this procedure to derive the total slant of the atmospheric delays from the GPS observations by using GAMIT/GLOBK. We initially obtained the pressure at the site to determine the hydrostatic zenith delays by using the Saastamoinen model, and the zenith wet delays were derived according to the same step used in partial derivatives. The mapping function was then used to map the zenith delays to slant directions, and the total slant delays used later were the slant hydrostatic delays plus the slant wet delays. The observed atmospheric delays varied from approximately 2 m to 100 m, and such delays increased as the angle of elevation β decreased. The signals received at low β would have travelled through significant gradients of tropospheric refractivity, carrying most of the profile information on atmospheric refractivity, particularly <5°[Gaikovich and Sumin, 1986].

The VMF1 mapping function [Boehm et al., 2006] was employed in the mapping process. The VMF1 mapping function recommended by the International GNSS Service (IGS) [Steigenberger et al., 2009] is derived from ECMWF (European Centre for Medium-Range Weather Forecasts) operational analysis data. This function is one of the most accurate and reliable mapping functions to date. However, the use of this function at very low β(<3°) is dangerous [Boehm et al., 2007]. The ray-tracing model is a simplified mathematical model that does not involve the horizontal gradients of refractivity, the phenomenon of multipath, or atmospheric diffraction. This representative error is particularly significant; that is, the gap between simulated atmospheric and real delays is large at low β. Figure 5 shows an example of the inconsistency between the observed atmospheric delays and the delays simulated by the ray-tracing model from nearby radiosonde measurements. In this case, the radiosonde measurements were collected at the same time via GPS observations. The observed atmospheric delays and the delays simulated from radiosonde measurements should be the same because such delays represent the same atmospheric conditions. However, the observed delays increased more rapidly than the simulated delays at very low βbecause only observations at βfrom 3° to 5° were used in the practical retrieval.

Figure 5.  An example of inconsistency between the observation of the GPS receiver and delays simulated by the ray-tracing model.

The retrieval process for the practical retrieval is summarized in Fig. 6. In the practical retrieval, the initial refractivity profile was calculated from the CIRA86aQ_UoG data in May, and the search space was restrained using five years of COSMIC measurements from May 2007 to 2011 at (35°-45°N, 110°-120°E) (Fig. 7).

The results from four days (7, 8, 20 and 22 May 2011) of observations are presented as examples. On each day, we conducted two retrievals at 1130 UTC and 2330 UTC when the radiosonde was released. To obtain the hourly mean of the atmospheric refractivity, we allotted 1 hour of GPS measurements at 1130 UTC and 2330 UTC.

Figure 6.  Flowchart of the retrieval process.

Figure 7.  The initial value and the search space for the practical retrieval.

Figure 8.  Examples of practical retrieval results compared with the initial profile and nearby radiosonde measurements.

Figure 8 shows the inversion results of the eight cases and a comparison with the nearby radiosonde measurements (39.93°N, 116.28°E). The initial profiles differed from the real refractivity distributions, particularly in the near-surface layers. The initial profiles were then improved at different degrees after the retrieval. The inversion results generally represented the real status of the vertical refractivity structure. In these cases, a minor error was found above 7 km, and a larger error possibly appeared at lower altitudes where abrupt changes in vertical refractive gradients occurred. Figure 9 shows examples of simulated atmospheric delays of the CIRA profile, the inversion result, and the corresponding observations. The simulated delays of the CIRA profile were 3-5 m smaller than the observations. By contrast, the simulated delays of the inversion result had slight differences from the observations.

We studied another 68 cases using observations from May and June 2011 to estimate retrieval accuracy statistically. The relative error of these 76 cases (Fig. 10) ranged from -3% to 5%, which was reasonably larger than that of the error in theoretical cases. The retrieved refractivity was likely larger than the radiosonde measurements (as indicated by the positive average relative error in Fig. 10). The lower troposphere, particularly below 4 km, indicated a larger error than the upper levels.

Diurnal variations in the refractivity profile per day with a temporal resolution of 1 hour can be shown based on error statistics. To demonstrate these diurnal variations, we initially performed retrievals for 38 days in May and June 2011 and obtained the hourly atmospheric refractivity profile. We then calculated the mean refractivity from all of those profiles. For each profile, the deviation from the mean was calculated. The deviations were grouped based on local time (0000 LST to 2300 LST) to derive the mean deviation at each hour. The average deviation of local refractivity is shown in Fig. 11. The largest variations appeared at altitudes of 0 km to 4 km, and the largest gap between the maximum deviation and the minimum deviation was approximately 20 N-units at the surface. These trends were reasonable because temperature and water vapor can change dramatically at altitudes on and near the surface. The results are also consistent with those of a previous study [Jie et al., 2009], in which the diurnal variations in refractivity at the tropical ocean from 0 to 10 km were analyzed based on COSMIC data. However, some differences were observed. The variations in our study showed greater magnitude, possibly because of higher temperatures and greater water vapor amplitudes on the continents than in the oceans. Another reason can be attributed to the difference between the sample number and database used in the two studies. This difference suggests that the retrieval accuracy is adequate to demonstrate the hourly variation in local refractivity profiles per day.

In the practical retrieval, the error shown in Fig. 8 was mainly caused by two aspects other than the inconsistencies of time and the locations of the GPS receiver and the radiosonde station. The first aspect indicates that the error can be brought in by the retrieval method. A higher vertical resolution and a smaller search step yielded more accurate results. However, this procedure requires additional computation hours, which may reduce the temporal resolution. The second aspect considers the characteristics of atmospheric delays. Atmospheric delays are an integral effect of refractivity gradients along a signal path. For this reason, complex vertical gradients of atmospheric refractivity along signal paths probably result in an atmospheric delay similar to that caused by a refractivity profile with relatively “medium” gradients. Figure 12 illustrates the two refractivity profiles and their respective atmospheric delays. The profile indicated by the dashed line greatly differs from that indicated by the solid line at altitudes of 2.5 to 3 km. However, the delays these profiles caused were almost the same, indicating that the atmospheric delays may not be sufficient to retrieve a complicated vertical distribution of refractivity with a positive gradient. This phenomenon can also be observed in Fig. 8. For example, Fig. 8a (2) and 8b (2) show that a positive gradient that appeared at approximately 3 km resulted in a large error because the gradient offset the excess signal path induced by the negative gradient. In Fig. 8a (1) (at heights from 0 to 4 km) and 8d (2) (at heights from 0 to 2 km), a “medium” refractivity profile was derived from the retrieval method instead of a profile with a positive gradient.

Figure 9.  Examples of atmospheric delays simulated from the initial profile and the inversion results on 20 March 2011, compared with measurements of GPS receivers.

Figure 10.  Statistical relative error of the practical retrieval. Triangles indicate the average of relative error on each level and the black line and dots indicate the error range.

Figure 11.  Diurnal variation statistics: (a) deviation from mean refractivity profile from 0000 to 1100 LST; (b) deviation from mean refractivity profile from 1200 to 2300 LST.

Figure 12.  An example of the deficiency of atmospheric delays in retrieving a positive refractivity gradient: (a) shows the radiosonde observation at 2330 UTC 7 March 2011 and an artificially constructed refractivity profile; (b) shows the atmospheric delays simulated from the two profiles in (a).

Errors can also be introduced by the observations (the data in O-files and M-files sent to GAMIT) and the calculation processes (Saastamoinen model, mapping function etc.) in the processing software. [Lowry et al., 2002] quantified the error contributions from the GPS processing software (e.g., estimation of hydrostatic and wet tropospheric delay parameters) and the GPS receiver (e.g., receiver clock and multipath), which are independent of the retrieval method and could thus be used in our study.

The present reported study focused on estimating the possibility of deriving a local refractivity profile with a high temporal resolution from ground-based GPS measurements. The study also focused on the limitations of the method in practical use. We established an 11-level stratification model at altitudes of 0-10 km with a vertical resolution of 1 km. We used an exhaustive search method to retrieve the vertical distribution of tropospheric refractivity. After verification, the methodology was used to retrieve the refractivity profiles from the observations of a ground-based GPS receiver. The inversion results were compared with radiosonde measurements, and the statistical error was evaluated. The error sources were discussed.

With ground-based GPS observations at low elevation angles, the retrieval method successfully demonstrated the general vertical distribution of local refractivity and its daily variation. This retrieval method is also a low-cost, stable, and high-temporal-resolution approach to detecting atmospheric refractivity despite the ambiguities in reflecting positive refractive gradients in a vertical direction. Further studies are needed to resolve the complicated vertical refractive gradients (particularly positive gradients) and reduce the inversion error in the lower troposphere.