Abstract

Quality control (QC) is among the most important steps in any data processing. These steps are elaborated for high-vertical-resolution radiosonde datasets that were gathered and analyzed to study atmospheric winds. The database is composed of different radiosonde wind-finding systems (WFSs), including radio theodolite, Loran C, and GPS. Inspection of this database, particularly for wind, wind shear, and ascent height increments (dz), showed a nonnegligible amount of outliers in radio theodolite data as compared to the two other WFSs, thus denoting quality differences between the various systems. An effective statistical QC (SQC) is then developed to isolate and eliminate outliers from the more realistic observations. Improving the accuracy of the radio theodolite WFS is critical to the derivation of the vertical motion and the vertical gradients of the horizontal wind—that is, wind shear—mainly because of the direct dependence of these quantities on dz. Based on the climatological distribution of the quality-controlled dz, a new approach is suggested to estimate these wind quantities for radio theodolite data. The approach is validated with the high-quality modern WFSs (Loran C and GPS). Although initially of reduced quality, applying SQC and using the climatological mean dz of 12-s smoothed radio theodolite profiles shows very good improvement in the climatological wind analyses of radio theodolite WFSs. Notably, the climatologies of ascent rate, vertical motion, horizontal wind, and vertical shear now look comparable for the various WFSs. Thus, the SQC processing steps prove essential and may be extended to other variables and measurement systems.

1. Introduction

Radiosonde wind vectors (speed and direction) are derived from the successive horizontal positions of the radiosonde balloon. Thus, the accuracy of the observed wind is intimately related to the accurate measurement of the radiosonde position. Position calculations are performed using different wind-finding systems (WFSs). Here we focus mainly on three WFS types—radio theodolite, long-range navigation (loran), and global positioning system (GPS)—that are used to collect the high-resolution radiosonde database used in the context of the preparation of the European Space Agency’s (ESA) Atmospheric Dynamic Mission Aeolus (ADM-Aeolus; Stoffelen et al. 2005). This database, illustrated in Fig. 1, is utilized following the quality control presented in Houchi et al. (2010) and Marseille et al. (2011). An overview of the measurement principle of the three WFSs is given in section 2b.

Fig. 1.

The geographical locations of analyzed high-resolution radiosonde datasets: SPARC (circles), BADC (hexagrams), AMMA (diamonds), and De Bilt (square) over different climate regions: tropics, subtropics, midlatitudes, and polar. Note the orography in the map that may explain the appearance of lee waves for some stations, in particular over the Rocky Mountains. The right legend bar indicates the altitude of the earth surface above sea level.

Fig. 1.

The geographical locations of analyzed high-resolution radiosonde datasets: SPARC (circles), BADC (hexagrams), AMMA (diamonds), and De Bilt (square) over different climate regions: tropics, subtropics, midlatitudes, and polar. Note the orography in the map that may explain the appearance of lee waves for some stations, in particular over the Rocky Mountains. The right legend bar indicates the altitude of the earth surface above sea level.

Although Stratosphere–Troposphere Processes and their Role in Climate (SPARC) radio theodolite data are subject to a strict quality control (NCDC 1998) at the University Corporation for Atmospheric Research and the Joint Office for Science Support (UCAR-JOSS), a nonnegligible amount of unrealistic observations (outliers) are observed in these data, as compared to Loran C and GPS radiosonde WFSs. This indicates mainly the difference in the accuracy between these three radiosonde WFSs. SPARC data are widely used in many studies (Wang and Geller 2003; Gong et al. 2008; among others) but with limited references to quality control (QC) aspects. The goal in this paper is to present a statistical QC (SQC) method for winds, wind shear, and ascent height increments (dz) of high-resolution radiosondes, and to demonstrate its effectiveness in isolating and removing outliers from a realistic probability density distribution (PDF) of the observations. Because of the moderate accuracy of the radio theodolite data as compared to Loran C and GPS WFS data, the derivation of the meteorological quantities depending on dz is very challenging. Thus, to estimate such quantities for radio-theodolite-based SPARC data, two new approaches are proposed in this manuscript. The first approach is to estimate the vertical motion based on the ascent rate statistics, and the second is to establish wind shear climatology. Both approaches, which use the statistics (mean and percentiles) of the QC-ed dz, are verified and validated against more accurate measurements of Loran C and GPS WFSs from Africa Monsoon Multidisciplinary Analysis (AMMA) projects, the British Atmospheric Data Centre (BADC), and De Bilt station (Fig. 1).

The different high-resolution datasets used and their related WFSs are presented in section 2, this includes an overview of the initial QC of SPARC data applied at UCAR-JOSS. In section 3, subsequent analysis of the different radiosonde WFS types shows the pronounced differences in the distribution of vertical differences (du, ) of the horizontal wind components (u, υ) and dz (i.e., the sampling interval). Then, sets of typical observed raw profiles of horizontal wind and their corresponding shear profiles for various WFSs are shown. In section 4, the developed SQC method is described in detail, while it is graphically illustrated in section 5 for the different analyzed quantities: u and υ components, their vertical shears, and dz. This includes the comparison of z and dz statistics for the different WFSs. In section 6, approaches are proposed to statistically estimate (i) the ascent rate and vertical motion and (ii) wind shear. Overall results and benefits of the SQC processing steps are summarized in section 7.

2. High-resolution radiosonde database

In this section, we present the high-resolution radiosonde database that we utilized and the different WFSs used to track and collect data from balloon radiosondes. Special attention is paid to the different WFSs; thus, derivation methods of vector winds for each WFS are described. An overview of the quality of the ensemble datasets is presented.

a. High-resolution radiosonde data and QC status

The Royal Netherlands Meteorological Institute (KNMI) high-resolution database is being developed in the framework of the ESA Aeolus mission with radiosondes, initially from the SPARC project dataset (Chanin 1995) and subsequently extended to more datasets from the BADC, the AMMA, and De Bilt station. But it is still being extended with other data sources, for example, from Météo-France and Finland. In this study, the first four datasets have been exploited, that is, SPARC, BADC, AMMA, and De Bilt station. The SPARC dataset is based on radio theodolite, while there is mixing of Loran C and GPS WFSs for the other datasets. Radiosonde data are available at high temporal resolution: 6 s (~30 m) for SPARC data, 2 s (~10 m) for BADC and AMMA, and 10 s (~50 m) for De Bilt in 2006. Values in parentheses denote the corresponding spatial vertical resolution, by considering the approximate radiosonde balloon ascent rate of 5 m s–1. The balloon is designed to rise at a more or less constant rate, because as the air density and thus friction per area decreases, its cross section gets larger (Brock and Richardson 2001, 223–230; WMO 2008b).

Moreover, SPARC data available for the period 1998–2006 are quasi regularly given at 0000 and 1200 UTC in ASCII files. The data records contain 22 fields including time from release, pressure, temperature, dewpoint, relative humidity, u and υ wind components, wind speed and direction, ascent rate, balloon position data, and altitude. BADC and AMMA radio soundings of 2-s raw resolution constitute a collection of about 21 stations each. The former are mainly dedicated to operational weather prediction, while the latter are dedicated to a campaign study of the African monsoon. Contrary to the analyzed SPARC period (1998–2006), these two datasets appear very irregular and, for example, unavailable for several years. Thus, after assessment of the distribution in time of all available BADC and AMMA radiosoundings, only few stations appear to be useful for climatology studies—that is, those sounding stations that are regular and continuous in time. Over 2006, for instance, only nine BADC and six AMMA stations (available at both 1200 and 0000 UTC) were selected as valuable. At 0600 and 1800 UTC, the number of useful stations is even much lower in 2006—two to six stations for BADC and only one station for AMMA—depending on the range of the targeted period. Also, the data from station De Bilt are analyzed, which are regular. This station uses both Loran C (during the night) and GPS (during the day) WFSs.

For SPARC radio theodolite datasets of both radiosonde brands (Sippican and Vaisala), used to derive winds at 6-s sampling, it is observed that elevation and azimuth angles are biased by elevation angle oscillations, which occur occasionally. Consequently, this leads to large oscillations in wind velocity, at low elevation angles in particular. Therefore, these measured parameters, together with the pressure, temperature, and relative humidity (PTU) quantities, are QC-ed at the UCAR-JOSS to remove outliers before deriving the wind components. A full description of the QC method is available in NCDC (1998), but we succinctly summarize the most important information on the joint UCAR-JOSS QC. The data underwent a two-stage QC process as follows:

  • QC-1. An internal consistency check that includes two inspections: “reasonable” limit checks on all parameters as defined by WMO (DiMego et al. 1985) and rate-of-change checks on temperature, pressure, and ascent rate (Schwartz and Govett 1992; DiMego et al. 1985).

  • QC-2. Each sounding endured a visual examination verifying those parameters that are too variable for automatic checks: wind speed, wind direction, and relative humidity. In addition, the verification of the former automatic QC flag (QC-1) is allowed at this stage.

Regarding the QC of the other WFS datasets, based on Loran C and GPS, no quality control has been applied on these raw data. Nevertheless, at first view their quality appears to be much better than the radio theodolite SPARC data, as can be seen in the section 3. So, the limited accuracy in most of data based on radio theodolite called for quality enhancement measures and the development of a new SQC to get rid of outliers. Although initially developed for SPARC data, this statistical method will also be applied to the sparser non-SPARC datasets.

b. Radiosonde wind-finding systems

As mentioned previously, radiosonde WFSs use different ways to determine the horizontal position of the balloon and, consequently, the derivation of the wind components (u, υ, and w) is also different. For the vertical position, apart from GPS WFS, the balloon vertical ascent height is generally determined by using the hypsometric relation [Eq. (1):

 
formula

where is the geopotential height at level i ( in m), R is the specific gas constant for air [287.06 (J kg–1 K–1) for dry air], is the mean virtual temperature (K) for the layer between two successive isobaric surfaces at pressure levels and , and g is the gravitational acceleration (m s–2), with a standard value of 9.806 65 according to WMO technical regulations. Note that Eq. (1) allows the determination of the ascent rate and vertical wind after removing the estimated net balloon lifting (see end of section 5).

The theodolite WFS is one of the first and oldest WFS. Initially the instrument was optical, which makes the height to which the balloon can be tracked dependent on many factors, such as the speed of winds aloft; obstruction in the light path between the balloon and the observer, mainly cloud; and the bursting point of the balloon. During a night observation, radiosonde tracking is accomplished by attaching a small self-illuminating light stick or a small battery-powered lighting unit to the balloon (OFCM 2006). Fortunately, since the introduction of the radar, balloon tracking has become much more simple and accurate—hence, the name radio theodolite. The radio theodolite system requires the determination of the azimuth angle and two of the following three quantities—height, elevation angle, and range—during the balloon ascent through the air until it bursts. With this information, the balloon position at different atmospheric levels can be determined and, consequently, the wind components. However, the measurement of the range is only possible when the radiosonde is equipped with a reflector, to be tracked by the radar.

Loran is a globally distributed terrestrial radio-navigation system using low-frequency radio transmitters, introduced in 1958. The detection system uses multiple transmitters, known as hyperbolic positioning or multilateration. In case of the radiosonde, three transmitters are generally used to determine the position and the speed of the balloon receiver (Fang 1986). Loran C operates in the low-frequency range of 90–110 kHz, and it is widely used in many countries to determine, for instance, the location of an aircraft or a ship. Loran, which started to lose popularity with the invention and growing use of GPS, is becoming more popular again, mainly to serve as a backup and land-based alternative to GPS or any other global navigation satellite system (GNSS). Loran C is the latest version before the recent introduction of the modernized version, named eLoran (enhanced loran). A prototype of eLoran was introduced in 2008 (Basker and Williams 2010). The interest for this type of navigation system grows with the introduction of the new modernized version, eLoran, which is independent, dissimilar, and complementary to GNSS (International Loran Association 2007).

GPS is currently the most widely used radiosonde WFS. The principle of GPS is quite similar to Loran C but with the multilateration method applied spatially in three dimensions. GPS receivers locate at least three satellites and estimate the range to each. This information is then used to deduce the receiver location in Earth’s atmosphere, with an accuracy of 0.5 m in the horizontal and 30 m in the vertical (Brock and Richardson 2001, 223–230). In the case of radiosonde balloons, the procedure is repeated at regular time intervals to determine the successive locations of a balloon and thus derive the wind vector. Many available networks, using radiosonde WFSs different from GPS, are switching progressively to this modern GPS WFS. In the United States, for instance, the entire radio theodolite radiosonde network was progressively replaced with modern GPS radiosondes in recent years, as planned during the Replacement Network System Project (Facundo 2005). More detailed information on all existing WFSs can be found in WMO (2008a, 241–307). The vertical wind shear vector s is defined as the gradient of the horizontal wind vector v = (u, υ) with respect to vertical height z [Eq. (2)]:

 
formula

with i indicating the level number of wind or wind shear (each shear level is in between two wind levels); 1 = 1, …, m − 1 is for shear levels; and m is the total number of vertical wind levels, which depends on the vertical sampling used; and u and υ are the zonal and meridional wind components, respectively. Wind shear levels are plotted as a function of shear heights, .

3. Comparison of wind-finding systems

a. PDF profile of wind and height differences

To get an idea about the measurement accuracy of each WFS used in our database, we derive the difference between successive values in the ascent profiles of zonal (du) and meridional () winds and dz, denoted here as “difference profiles.” Then, the PDFs of individual profiles of these differences are determined, using bin sizes of 1 m s–1 for du and , and 10 m for dz. For a fair comparison, both the BADC and AMMA datasets of 2 s are resampled at the vertical resolution of SPARC data of 6 s. However, De Bilt data are kept at their raw vertical resolution of 10 s. The right plots of Figs. 2a,b illustrate these PDFs, where the relative frequencies of the occurrence of u, υ, and z differences within the mentioned range bin sizes are depicted. From these figures, we can clearly see a large dispersion in SPARC radio theodolite data as compared to the other WFSs (e.g., shown for GPS in De Bilt in Fig. 2b) for the three analyzed quantities. We can see in particular that only in SPARC radio theodolite data do outliers exist, that is, many values that fall far away from the main mode of the PDF. Though some dispersion may be attributed to meteorological events—such as extreme weather conditions in some climate regions, gravity waves at high altitude, etc.—many values appear unrealistic, that is, as outliers. It is important to mention that these findings for each WFS are very typical for most individual profiles. This indicates that radio theodolite data suffer from a relatively large number of outliers, and therefore an additional QC is necessary. It appears that above about 22 km, the accuracy in dz measurements is severely decaying due to a digitization issue in the SPARC radio theodolite data, as show in Fig. 2a. It is verified that this digitization issue originates from the pressure profiles, which are used together with temperature profiles to derive geopotential heights (ascent heights), z [see Eq. (1)]. The radio theodolite WFS is further known to have a limited accuracy for elevation angles below 17° (Vaisala 2002). This is due to the interference between signals received directly from the radiosondes and those received by reflection from adjacent surfaces, generally termed multipath interference (WMO 2008b). The amount of multipath interference depends very critically on the positioning of the antenna relative to adjacent reflecting surfaces, for example, whether the radio theodolite is positioned on a roof or on the ground.

Fig. 2.

(a) Typical radiosonde difference profile (left column) for (top to bottom) zonal (du) and meridional () winds, and ascent height increment (dz) and (next column) their corresponding relative distributions over range bins of 1 m s−1 for du and and 10 m for dz. (b) As in (a), but for a GPS WFS from De Bilt data, with a vertical resolution of 10 s (~50 m). This result is typical for a single profile for modern WFSs (GPS and Loran C) sampled at 6 s (~30 m) or more.

Fig. 2.

(a) Typical radiosonde difference profile (left column) for (top to bottom) zonal (du) and meridional () winds, and ascent height increment (dz) and (next column) their corresponding relative distributions over range bins of 1 m s−1 for du and and 10 m for dz. (b) As in (a), but for a GPS WFS from De Bilt data, with a vertical resolution of 10 s (~50 m). This result is typical for a single profile for modern WFSs (GPS and Loran C) sampled at 6 s (~30 m) or more.

b. Distributions of wind and wind shear profiles

We show some examples of the common raw observation profiles of zonal wind and wind shear as function of ascent height, at z and , respectively, first for SPARC radio theodolite and subsequently for the other types, Loran C and GPS at station De Bilt, and for the BADC and AMMA campaigns. Figure 3a shows a profile set of wind (top) and wind shear (bottom) for the complete year of 2006 at the tropical stations. The corresponding collocated ECMWF model data are also shown alongside for comparison; the model data are the 12-h forecast model fields interpolated to the radiosonde launch location. The ECMWF forecast field is preferred to the analysis to avoid incestuous comparison, since the latter contains the current radiosonde data. Large spikes in SPARC data denote the presence of outliers, which should be removed before generating any statistics. In Fig. 3a, we can see a zoomed-in version. These plots demonstrate more clearly the global shape of the wind and wind shear profiles for both SPARC radio theodolite–radiosonde and the ECMWF model, which look pretty similar in some aspects for the wind profiles. However, the profile set of wind shear is substantially broader for the radiosondes. This is not due to the outliers but suggests that vertical wind variations are underestimated in the model (Houchi et al. 2010). The figure also provides a rough estimate of reasonable thresholds necessary for wind and wind shear quality controls. Figure 3b is similar to Fig. 3a and its zoomed-in version, but now for modern WFSs’ (Loran C and GPS) midlatitude BADC data in 2006. This dataset shows very few outliers in wind and shear values; only few values are probably unrealistic near the surface and in the upper troposphere. This suggests that the Loran C and GPS WFSs are of much better quality than the radio theodolite WFS. This is also seen in individual profiles based on Loran C and GPS WFSs, for example, for AMMA and De Bilt (not shown).

Fig. 3.

(a) Tropical (first row) zonal wind and (second row) wind shear profiles for (left column) the 11 SPARC radiosonde stations for 2006 at a vertical resolution of 6 s (~30 m) and (next column) the collocated ECMWF 12-h short-range forecast. The number of collocated radiosonde profiles used is ~80% [3217 profiles/(11 stations × 365 days)]. Notice that the ECMWF model resolution is interpolated to the radiosonde resolution. (rows three and four) As in (first and second row), but zoomed horizontally. (b) As in (a) (first and second row) for six BADC midlatitude stations over 2006 at a vertical resolution of 2 s (~10 m). The number of available profiles represents more than 83% [1834/(6 ×365)]. Notice the very low number of extreme outliers as compared to (a).

Fig. 3.

(a) Tropical (first row) zonal wind and (second row) wind shear profiles for (left column) the 11 SPARC radiosonde stations for 2006 at a vertical resolution of 6 s (~30 m) and (next column) the collocated ECMWF 12-h short-range forecast. The number of collocated radiosonde profiles used is ~80% [3217 profiles/(11 stations × 365 days)]. Notice that the ECMWF model resolution is interpolated to the radiosonde resolution. (rows three and four) As in (first and second row), but zoomed horizontally. (b) As in (a) (first and second row) for six BADC midlatitude stations over 2006 at a vertical resolution of 2 s (~10 m). The number of available profiles represents more than 83% [1834/(6 ×365)]. Notice the very low number of extreme outliers as compared to (a).

4. Statistical quality control procedure

a. Preamble

Although a strict QC is applied on SPARC radio theodolite data at UCAR-JOSS (see section 2), a visual check of the data shows several points that appear to be unrepresentative and what we suspect to be outliers. Some extreme values of wind and shear exceed the tolerance limit as defined by WMO (DiMego et al. 1985) and some other values have unrealistic changes within a 30-m vertical interval. The outliers cause generally large spikes in the profile distributions (seen in Fig. 3a) and its zoomed-in version (Fig. 4, 5a, 7a). To segregate and remove these outliers from the PDF, a SQC procedure, based on percentiles [Eq. (4)] and mean values of the PDF at each vertical level, has been developed. The percentile method sorts in ascending order the values of the observed variables collected at each 1-km vertical bin. If is a vector of n sorted values, then the values can be specified by percentile ranks , using Eq. (3). There are various definitions for percentiles in the literature; in this paper we adopt the definition as used in the MATLAB software (MATLAB 2010):

 
formula

with i = 1, …, n

Fig. 4.

Percentile profiles of (left) zonal wind and (right) wind shear for high-resolution radiosondes (top) before and (middle) after the automated QC, and (bottom) a zoomed-in version after QC. This is based on radio theodolite data, which contain nine SPARC tropical stations as shown in Fig. 1. The percentiles are monotonically distributed from 0 to 100 as indicated in the legend (starting from the left curve profile and going to the right profile curve). In each plot, only the added mean profile (very black) can overlap with the median profile 9 (gray) at certain altitudes, and these two latest compare well after the QC. Notice that at each 1-km vertical bin, only 0.1%–0.2% of data at maximum are removed by the SQC method.

Fig. 4.

Percentile profiles of (left) zonal wind and (right) wind shear for high-resolution radiosondes (top) before and (middle) after the automated QC, and (bottom) a zoomed-in version after QC. This is based on radio theodolite data, which contain nine SPARC tropical stations as shown in Fig. 1. The percentiles are monotonically distributed from 0 to 100 as indicated in the legend (starting from the left curve profile and going to the right profile curve). In each plot, only the added mean profile (very black) can overlap with the median profile 9 (gray) at certain altitudes, and these two latest compare well after the QC. Notice that at each 1-km vertical bin, only 0.1%–0.2% of data at maximum are removed by the SQC method.

Fig. 5.

(a) As in Fig. 4 (top and bottom), but for 37 SPARC (radio theodolite) subtropical stations (see Fig. 1). Notice the increased number of outliers compared to Fig. 4. (b) As in (a), but for mixed profiles of Loran C and GPS WFS from six BADC midlatitude stations.

Fig. 5.

(a) As in Fig. 4 (top and bottom), but for 37 SPARC (radio theodolite) subtropical stations (see Fig. 1). Notice the increased number of outliers compared to Fig. 4. (b) As in (a), but for mixed profiles of Loran C and GPS WFS from six BADC midlatitude stations.

For any given percentage p between 50/n% and 100–50/n%, the percentile value is defined using a linear interpolation between the closest ranks as follows:

 
formula

knowing that the minimum and maximum values in X correspond to the minimum and maximum percentiles, p = 50/n% and p = 100–50/n%, respectively.

The mean profiles are also calculated, since they can serve as a double check on the median profiles to represent the average state of the analyzed quantity. This may also give an indication on the persistence of the outliers, since the mean is more sensitive to outliers than the median. Globally, the SQC consists of determining thresholds for the observed variables to isolate and eliminate the outliers, that is, those observations falling outside the threshold values in the PDF. The threshold values are determined statistically at each 1-km vertical bin of the atmosphere using test criteria (see Fig. 6a). The different steps of the SQC method are detailed and illustrated with graphics in the next section.

Fig. 6.

(a) Flowchart of the SQC algorithm. (b) Illustration of the SQC on a set of zonal wind profiles for the tropical stations in 2006. Gray points at the edges of the different distributions with the altitude are the rejected wind values, while the black are the accepted points. Dots are used in the profiles in order to highlight the number of removed data points by the SQC.

Fig. 6.

(a) Flowchart of the SQC algorithm. (b) Illustration of the SQC on a set of zonal wind profiles for the tropical stations in 2006. Gray points at the edges of the different distributions with the altitude are the rejected wind values, while the black are the accepted points. Dots are used in the profiles in order to highlight the number of removed data points by the SQC.

b. Description of the statistical QC algorithm

Here we describe the most important steps of the SQC algorithm, which are summarized in the flowchart in Fig. 6a. Initially, some very extreme outliers are removed by using the tolerance limits as defined by WMO and cited previously. Following the flowchart, the first step in the SQC procedure is to collect a sufficient data amount (ideally one year or more) at each 1-km range bin and then compute and plot the percentiles [Eq. (4)] and means. This is done for both wind components u and υ, their corresponding shear values, and dz. The bin size may be modified–adapted to get the desired amount of data at each vertical bin for the QC statistics over a number of stations in a set time period, but here we keep it fixed to1 km. This is based on our observation that the number of data is almost constant throughout the troposphere and decreases only at the highest levels, because of the difference in the balloon-bursting height from one sounding to another (not shown).

Before performing the SQC on the wind profiles, initial wind shear profiles are approximated as a gradient of the horizontal wind over a constant vertical interval (dz), in order to avoid the digitization issue seen in the radio theodolite WFS. The dz is fixed here to 30 m, which corresponds to ~6-s time sampling and a mean ascent rate of 5 m s–1. So, this fixed height at each vertical level scales all shear values and percentile shear values in the same way and thus does not affect the SQC method based on percentile analysis. The exact wind shear is then recomputed afterward, using the obtained QC-ed products for wind components u and υ and mean dz. For example, QC-ed u is shown in Fig. 6b. In addition, when wind shear values are rejected, the corresponding underlying wind values are also rejected and vice versa. Once the final QC-ed wind profiles are obtained, they are used to establish the statistics, for example, shown in the right-bottom plots of Figs. 4, 5a for wind shear, successively, for 9 tropical and 37 midlatitude stations, respectively.

The density of the number of percentile profiles is enhanced at the extremes when approaching the 0 and 100 percentiles for each of the analyzed quantities, for example, u, υ, dz, etc. The subsequent percentiles are 17 discrete values in our case, as shown in Table 1, below:

Table 1.

Percentiles corresponding to profile numbers 1 to 17 (and plotted from left to right in the figures).

Percentiles corresponding to profile numbers 1 to 17 (and plotted from left to right in the figures).
Percentiles corresponding to profile numbers 1 to 17 (and plotted from left to right in the figures).

So, the percentiles represented by values from 1 to 9 (from 0 to 50) increase in value and cover the left side of distribution, while the percentiles for values of 17 to 9 (from 100 to 50) decrease in value and cover the right side of distribution; see the legends of Figs. 5, 7, for instance. The comparison of subsequent differences of percentile values will allow for determining thresholds N at each vertical level—for example, and for the left and right values of the distribution, respectively—defined as condition 1. Condition 2 and so forth will be the next subsequent differences, as shown in the flowchart in Fig. 6a. The N denotes the tolerated magnitude of the relative jump between two subsequent percentiles at the edges of the PDF. The N is statistically dependent on the data volume (time period and number of stations) present at each level, thus on the width of the percentiles. But, it is more critically dependent on the outlier likelihood, existing in each used dataset or subset.

Fig. 7.

(a) The dz profiles (left) before and (middle) after the SQC, and (right) the latter zoomed-in at each 1-km bin (as shown in dots for mean profiles) for radio theodolite SPARC profiles of 6 s (~30 m). (b),(c) As in (a), but for De Bilt (10 s, ~50 m) and AMMA (2 s, ~10 m), respectively, both based on Loran C and GPS WFSs. Legend as in Figs. 5, 6b for the percentile and mean profiles. Notice also the very small number of outliers in De Bilt and AMMA data indicating the excellent data quality.

Fig. 7.

(a) The dz profiles (left) before and (middle) after the SQC, and (right) the latter zoomed-in at each 1-km bin (as shown in dots for mean profiles) for radio theodolite SPARC profiles of 6 s (~30 m). (b),(c) As in (a), but for De Bilt (10 s, ~50 m) and AMMA (2 s, ~10 m), respectively, both based on Loran C and GPS WFSs. Legend as in Figs. 5, 6b for the percentile and mean profiles. Notice also the very small number of outliers in De Bilt and AMMA data indicating the excellent data quality.

The set value of N determines the trade-off between the rejection of representative extreme values (false alarm rate) and of unrepresentative outliers (outlier detection probability). These two characteristics are clearly illustrated in Figs. 4, 5a. An increase in the number of real outliers increases the magnitude of the jump in value between the extreme percentiles. But since the number of outliers is arbitrarily different from one dataset–subset to another, even for the same WFS, N has to be evaluated heuristically (involving trial and error and visual checks) at least once for each processed dataset. Once N is set for each dataset, the procedure of determining the percentile threshold and removing outliers at each atmospheric vertical level is automatic. We note, however, that the estimation of N should be reconsidered if the adopted bin size is different than the 1 km used here, since the underlying wind and outlier PDFs change with bin size. For one station, we have 365 profiles per year and about 12 000 points over a 1-km bin. For the tropical and polar stations, we get good results for N = 2, as can be seen in Fig. 6b. On the other hand, for other datasets with multiple stations, the proportion of data in each percentile profile is higher and outliers stand out more clearly at the extreme values, as seen in Fig. 5a as compared to Fig. 4. In these cases, condition 1 is sufficient to eliminate the outliers, and the factor N should be larger to restrict the passing from condition 1 to condition 2 in the QC procedure and avoid the elimination of realistic observations. So, for datasets where we have more than eight stations, the tests show that for N = 10, we remove only the points that are most probably outliers.

By setting pQC as the percentile that defines the threshold, it denotes the level of rejection at both sides of the PDF distribution. Since outliers appear on both sides of the distributions in similar amounts, we considered a symmetric QC for the determination of the percentile threshold values at each vertical level. As such, we relate pQC to a left and right PDF percentile counter, mm and pp, respectively, defined as follows: mm = 1 + pQC (for the left percentile) and pp = 17 − pQC (for the right percentile). So, for a given pQC, mm sets the percentile position number in ascending threshold value order (minimum accepted PDF value), while pp in descending threshold order (maximum accepted PDF value). For instance, when pQC is initiated at pQC = 0, the corresponding percentile position numbers are 1 (0%) and 17 (100%) and all values in the PDF are accepted. So, pQC is either equal to 0 (no outliers), 1 (0.1% outliers on each side of the PDF; i.e., 1/500 points removed), 2 (2/500 points removed), or 3 and so forth. Its value depends on the fulfillment or not of the conditions defined in the flowchart in Fig. 6a. In our tests, we observed that from pQC > 2, we start to remove many realistic observations at some range bins. So, the extreme three percentile levels—that is, accepting all, and 0.1–99.9 and 0.2–99.8 percentiles—are sufficient to reach our goal to remove the most harmful outliers from the statistics. So, in our datasets a maximum threshold level of pQC = 2 (i.e., 0.2% outliers on both sides) was sufficient to eliminate the outliers causing large and unreasonable jumps in the distributions. The distribution of the QC-ed products is plotted in the same figure as the raw data. The SQC method is applied for both SPARC and non-SPARC datasets, but it results in a very small amount of outliers in the latter. Overall, the QC results are illustrated graphically and discussed below in more detail for the winds, shears, and dz.

5. SQC results

In this section, the SQC applied to the horizontal wind components, their vertical shear, and dz is graphically illustrated and discussed accordingly. The SQC method can be run for different vertical aggregations. To reduce random noise and for a better comparison, most of the results are established with a smoothed 12-s (~60 m) vertical resolution rather than 6 s (~30 m) for each analyzed quantity. However, the QC-ed u, υ, and dz can be used to compute other important quantities, such as the drift, the vertical air motion, etc. The final QC-ed products are used in this paper to derive more reliable wind shear characteristics, including estimation of the vertical motion from the ascent rate statistics (section 6).

a. Zonal and meridional winds

Figure 4 illustrates the SQC method applied to the 2006 tropical wind and shear data of nine SPARC radiosounding stations. We can see from this figure that the outer percentiles of 0.1 and (99.9) isolate in most cases the values that present a very large jump (outliers). The 0.2 and (99.8) percentiles’ threshold is rarely used for rejection in this case, and generally it is only near the surface and lower stratosphere. Similarly, Figs. 5a,b illustrate the SQC of wind and shear data for 37 SPARC subtropical and 6 BADC midlatitude stations in 2006, respectively. Except in non-SPARC (modern WFSs) data shown in Fig. 5b, all SPARC radio theodolite often exhibited a large increase in values between subsequent percentile profiles at the edges of the percentile distributions between (0, and 0.1, and 99.9 and 100). This is not due to the increase in the number of profiles—thus, the increased amount of data collected at each 1-km vertical range bin—since this would better sample the percentile and presumable result in a smooth vertical line. Near the edges, the large jump between subsequent percentile profiles is rather due to the frequency of occurrence of extreme outliers, which is highly dependent on the data quality. This can be seen from Figs. 4, 5a for SPARC radio theodolite data as compared to Fig. 5b for modern WFSs. Indeed, according to these figures SPARC radio theodolite data show notable and typically a high occurrence of the frequency of outliers (big spikes). We found this to differ from one station to another and from one climate region to another, etc., even for the same data type.

Moreover, Figs. 4, 5a show not only that the size of the dataset used (number of profiles and thus the data amount collected at each vertical bin) affects the determination of the threshold values in the SQC procedure, but it also affects the magnitude of the percentile jumps at the extremes of the PDF and the associated SQC parameter, N (defined previously). This is due to the fact that the fraction of outliers does not increase systematically with the number of data at each vertical bin. Note that we presume that the PDF exists of a unimodal part with real values and a more wide (symmetric) PDF of outliers. As more data are used, both PDFs become better defined. Therefore, the percentile profiles generally become smoother and less erratic as more input data are used. However, note that big jumps between two outer percentiles are not only related to the amount of data but also due to the number of cumulative outliers at each level bin. So, for large datasets the N factor is mainly dependent on the data quality (measurement accuracy), which can also differ from one station to another and from one level to another for the same dataset. This justifies the necessity to determine this factor heuristically for each processed dataset, as already mentioned. Figure 5b shows the SQC for non-SPARC data (mixed Loran C and GPS) for six BADC midlatitude stations of the BADC dataset. In comparison with SPARC data, these results do not show a large jump between the very closely separated outer percentiles at most atmospheric vertical levels. Consequently, it is seen during the SQC of the non-SPARC data that the percentile thresholds pQC remain mostly equal to zero, thereby indicating the good data quality and negligible amount of outliers. Figure 6b shows an example of the final QC-ed product for zonal wind profiles (black), excluding the removed outliers (gray) from the raw data (black and gray). Globally, the results show that most outliers fall in the outer 0.1 (99.9) percentile profile for the east (west) wind and shear values at most vertical levels, with (N) set to 2 here.

b. Ascent height increments dz

We illustrate here the results of the SQC on dz and on the ascent height (z). These results are used in section 6 to estimate more reliably the wind shear, ascent rate, and vertical motion. Figure 7 shows the percentile statistics of dz before and after SQC at each 1-km vertical bin, for SPARC, De Bilt, and AMMA data. The 37 subtropical stations are used over 2006 for the SPARC dz and z plots, while only one station each for De Bilt and AMMA is used. The statistics of radio theodolite radiosonde data (SPARC) show many unrepresentative values (outliers) compared to Loran C and GPS radiosonde data (De Bilt and AMMA), in particular near the surface. Moreover, the digitization issue is observed at high altitude in SPARC data, and that originates from pressure profiles as stated in section 3a. Although the SQC method is effective in removing outliers at most atmospheric levels in dz profiles, above about 22 km the effect of the pressure digitization; which is amplified by the gravity waves occurrence, persists. To reduce this digitization problem and to lower the random noise and the gravity wave effects, the SPARC data are degraded from 6 s (~30 m) to 12 s (~60 m) vertical resolution, as shown in Fig. 8a for the statistics of dz and z.

Fig. 8.

Mean (dots) and percentile (lines) profiles for the QC-ed (left) dz and (right) z for (a) SPARC radio theodolite, (b) mixed Loran C, and (c) GPS datasets at De Bilt and AMMA GPS datasets; these statistics are established over each 1-km vertical bin. For a better comparison, profiles of each WFS dataset are resampled to about similar vertical resolution, i.e., De Bilt at 10 s (~50 m), SPARC at 12 s (~60 m), and AMMA at 12 s (~60 m). Notice the residual perturbations in SPARC data above ~22 km, even after smoothing, because of the digitization effect.

Fig. 8.

Mean (dots) and percentile (lines) profiles for the QC-ed (left) dz and (right) z for (a) SPARC radio theodolite, (b) mixed Loran C, and (c) GPS datasets at De Bilt and AMMA GPS datasets; these statistics are established over each 1-km vertical bin. For a better comparison, profiles of each WFS dataset are resampled to about similar vertical resolution, i.e., De Bilt at 10 s (~50 m), SPARC at 12 s (~60 m), and AMMA at 12 s (~60 m). Notice the residual perturbations in SPARC data above ~22 km, even after smoothing, because of the digitization effect.

In these last plots, we still see some residual perturbations in SPARC z (thus in dz) above ~22 km even after smoothing. This indicates once more that the systematic digitization effect is not easy to remove, and thus there remains large oscillations in the percentile statistics in dz plots. Also, some imperfections are observed close to the surface. These are likely due to the known radio theodolite WFS limitation for elevation angles below 17° (Vaisala 2002). However, similar dz and z analyses for the two other WFS datasets of De Bilt and AMMA (Figs. 8b,c) show a quasi-linear trend in all the percentiles and mean profiles, suggesting a better measurement quality. Apart from the abovementioned situations, we can see at most atmospheric levels that after the SQC, radio theodolite SPARC dz PDFs compare well with the dz PDFs of the modern WFS (Loran C and GPS) used in AMMA and De Bilt.

6. Vertical motion and wind shear estimation

In this section, we propose a scheme to estimate wind-derived quantities more reliably—the ascent rate, the vertical air motion, and wind shear.

a. Ascent rate and vertical air motion

Following our discussion in section 5b on the statistics of dz and z shown in Fig. 8, we propose here an approach to approximate vertical air motion based on the ascent rate statistics. The percentile statistics plots for dz may be seen as the balloon ascent rate statistics, since the latter is defined as the dz over a fixed time sampling Δt—for example, 12 s for SPARC data—dz depends on the vertical air motion w, the balloon ascent rate in still air AR, and the time sampling Δt: , such that

 
formula

This is similar to the model proposed by Wang et al. (2009). Their model enables the extraction of the air vertical velocity from radiosonde data, by considering the balloon ascent rate as a result of two contributions, one representing the balloon ascent in still air and the second the vertical air motion. So, if the balloon ascent rate in the absence of vertical winds may be assumed constant, then air vertical velocity is obtained by subtracting the ascent rate in still air from the actual ascent rate.

So, the spread in the dz PDF around the median is proportional to the spread in the vertical motion PDF. The AMMA and De Bilt radiosondes show small spread in the stratosphere between 15- and 20-km height. Since vertical air motion at these levels is limited, this spread may be essentially due to variations in balloon lift speed. In turn, such variations at one particular station are typically caused by differences in Helium (He) amounts, as the balloon type and thus launch mass is generally fixed. Variations in He amount cause similar variations in speed at all levels, while the mean ascent speed varies smoothly by a few tens of meters per second with height. The typical spread (1σ value) in the dz PDFs represented in Fig. 8 is about 3 m, corresponding to a vertical motion spread of 0.25 m s–1. Note, however, that the 90% percentile in the troposphere deviates relatively much from the other percentiles, particularly for the tropical AMMA data. This PDF skewness toward larger dz and thus larger upward motion is likely associated with the diabatic heating of ascending air due to condensation. The tropical tropospheric deviation of the 90% percentile of dz from the median is about 10 m, corresponding with an upward air motion of roughly 0.8 m s–1. Most other known lift effects either broaden the PDF on both sides, such as turbulence, or slow down the balloon due to precipitation loading, and do not explain a upward-only motion anomaly. Following this argument, tropospheric PDF broadening due to turbulence and moist convection may be visible in the Loran C and GPS dz PDFs. The radio theodolite PDFs show more variation in the stratosphere, indicating larger variation in balloon mass and/or He. Therefore, diabatic effects appear less asymmetric and less pronounced in the SPARC climatology.

b. Wind shear

The direct computation of wind shear using the standard definition [Eq. (2)] is not reliable for SPARC radio theodolite data as mentioned previously, even after the SQC of winds and dz. So, we propose here an approach for the vertical wind shear computation for the radio theodolite data by using the climatological mean dz, rather than an instantaneous dz profile. This method is verified with modern WFS data (loran and GPS), since for the latter the standard calculation of wind shear appears accurate. The comparison of the two methods is tested and shown for De Bilt data over 2008 (Loran C and GPS data) in Fig. 9. Absolute shear values are considered for the statistics in this figure and plotted in decadic logarithmic scales (along the horizontal axis) as a function of shear height . It is clear that this approximation ignores variability in dz due to vertical motion [Eq. (5)] of the balloon, which corresponds to about 5%. Moreover, diabatic effects may occasionally further enhance the deviation from a nominal balloon lift. However, the results show that wind shear computed with the two methods leads to very similar results, supporting therefore our approach for wind shear estimation for radio theodolite data.

Fig. 9.

Wind shear statistics [mean (dots) and percentile (lines)] profiles for De Bilt radiosonde 10-s (~50 m) resolution computed with two methods: (left) according to the standard definition [Eq. (2)] and (right) using the climatological mean of the QC-ed dz at each 1-km bin. Notice the logarithmic scales along the x axis.

Fig. 9.

Wind shear statistics [mean (dots) and percentile (lines)] profiles for De Bilt radiosonde 10-s (~50 m) resolution computed with two methods: (left) according to the standard definition [Eq. (2)] and (right) using the climatological mean of the QC-ed dz at each 1-km bin. Notice the logarithmic scales along the x axis.

7. Conclusions

In this paper we present an efficient statistical QC method to isolate and eliminate outliers from a large distribution of “realistic” profiles of wind, wind shear, and ascent height. The SQC method uses fine percentiles—that is, very small percentage intervals—at the extremities of the PDFs of each analyzed quantity at each atmospheric level (1-km vertical bin size is used in this manuscript). The SQC method is developed for radiosondes based on the radio theodolite WFS, but it has also been applied to loran and GPS data. For the two latter WFSs, very few outliers are detected, because of the good quality of these measurements. Globally, the SQC method shows a beneficial impact on the different analyses (on individual profiles and statistics) applied on the meteorological variables as studied in this paper. After QC, the climatology analyses for radio theodolite WFSs and more modern WFSs data appear to be comparable. However, above 25 km the altitude comparison is more complex, since the data amount can decrease substantially due to the different heights reached by the individual soundings and the increased presence of gravity waves.

A statistical approach is introduced in this paper to correct radio theodolite WFS data. In terms of weather, radiosondes based on radio theodolite may be difficult to interpret, as seen for instance here for the ascent rate, vertical air motion, and wind shear. This is due to the poor reporting accuracy of such WFS, particularly in the ascent height increments due to the pressure measurement digitization. We estimate the vertical air motion from the ascent rate climatology and in addition derive wind shear profiles and climatology. In all cases the climatological statistics of the QC-ed ascent height increments are used. The posteriori height increment estimation method using climatology has been validated with data from the more accurate WFSs, that is, GPS and Loran C. With this approach, the climatologies of the ascent rate, the horizontal wind components, and the vertical shear for radio theodolite WFSs data become comparable with the Loran C and GPS WFS results. However, a skewness of the tropospheric vertical motion PDF, probably due to latent heating, is clearly visible in GPS and Loran C but much less so in the radio theodolite WFS data. The SQC approach may also be beneficial in the case of individual profiles—for instance, to characterize physical phenomena, such as turbulence, gravity waves, and so forth. A similar approach may also be applied to the other meteorological variables provided by radiosonde instruments, that is, pressure, temperature, and humidity. But new threshold criteria have to be derived for these variables. In terms of applications, these variables are of particular interest for high-vertical-resolution cloud analysis of radiosonde measurements. This would, for example, be useful for preparing the ESA ADM-Aeolus Doppler wind lidar mission, in which wind profile quality will depend on the combined wind and cloud structure occurrence.

Acknowledgments

We thank first the people who helped us to gather the high-resolution radiosonde data and who provided support: Stefan Liess (SPARC data), Guillaume Brissebrat (AMMA data), and David Hooper and his team (BADC, Met Office); and Marc Allaart, Jitze van der Meulen, Henk Klein Baltink, and Ge Verver (De Bilt data, KNMI). This work is supported by a grant from the Dutch Ministry of Transport, Public Works and Water Management. SPARC high-resolution radiosondes data in this paper were provided by the Canadian Stratospheric Processes and their Role in Climate Network (C-SPARC), funded by the Canadian Foundation for Climate and Atmospheric Sciences and the Canadian Space Agency. BADC high-resolution radiosonde data are provided courtesy of the Met Office (United Kingdom) through the British Atmospheric Data Centre. Based on a French initiative, AMMA was built by an international scientific group and was funded by a large number of agencies, especially from France, the United Kingdom, the United States, and Africa. It has been the beneficiary of a major financial contribution from the European Community’s Sixth Framework Programme. Detailed information on scientific coordination and funding is available on the AMMA International website (http://www.amma-international.org).

REFERENCES

REFERENCES
Basker
,
S.
, and
P.
Williams
,
2010
: Navigating eLoran: Challenges and the way forward. Proc. XVIIth IALA Conf., Cape Town, South Africa, International Association of Marine Aids to Navigation and Lighthouse Authorities.
Brock
,
F. V.
, and
S. J.
Richardson
,
2001
:
Meteorological Measurement Systems.
Oxford University Press, 304 pp.
Chanin
,
M.-L.
,
1995
:
The SPARC programme
.
Phys. Chem. Earth
,
20
, 21–32, doi:.
DiMego
,
J. G.
,
P. A.
Phoebus
, and
J. E.
McDonnel
,
1985
: Data processing and quality control for optimum interpolation analysis at the National Meteorological Center. NMC Office Note 306, 38 pp.
Fang
,
B. T.
,
1986
:
Trilateration and extension to Global Positioning System navigation
.
J. Guid. Control Dyn.
,
9
,
715
717
, doi:.
Facundo
,
J.
,
2005
: Update on the implementation of the National Weather Service’s radiosonde replacement system. Ninth Symp. on Integrated Observing and Assimilation Systems for the Atmosphere, Oceans, and Land Surface, San Diego, CA, Amer. Meteor. Soc., 11.4. [Available online at https://ams.confex.com/ams/Annual2005/techprogram/paper_84775.htm.]
Gong
,
J.
,
M. A.
Geller
, and
L.
Wang
,
2008
:
Source spectra information derived from U.S. high-resolution radiosonde data
.
J. Geophys. Res.
,
113
, D10106, doi:.
Houchi
,
K.
,
A.
Stoffelen
,
G.
Marseille
, and
J.
De Kloe
,
2010
:
Comparison of wind and wind shear climatologies derived from high-resolution radiosondes and the ECMWF model
.
J. Geophys. Res.
,
115
, D22123, doi:.
International Loran Association
,
2007
: Enhanced Loran (eLoran) definition document. Version 0.1, Rep., 17 pp.
Marseille
,
G.
,
K.
Houchi
,
J. de
Kloe
, and
A.
Stoffelen
,
2011
:
The definition of an atmospheric database for Aeolus
.
Atmos. Meas. Tech.
,
4
,
67
88
, doi:.
MATLAB
,
2010
: Descriptive statistics. Statistics toolbox for use with Matlab. User’s guide, MathWorks, 3-1–3-20. [Available online at http://www.pi.ingv.it/~longo/CorsoMatlab/OriginalManuals/stats.pdf.]
NCDC
,
1998
: Data and documentation for rawinsonde 6-second dataset. Tech. Rep. TD6211, 12 pp. [Available online at ftp://atmos.sparc.sunysb.edu/pub/sparc/hres/6_second/TD6211.pdf.]
OFCM
,
2006
: Processing sounding data. Rawinsonde and pibal observations, FCM-H3-1997, Federal Meteorological Handbook 3, FCM-H3-1997, 5-1–5-10. [Available online at http://www.ofcm.gov/fmh3/pdf/00-entire-FMH3.pdf.]
Schwartz
,
B.
, and
M.
Govett
,
1992
: A hydrostatically consistent North American radiosonde data base at the Forecast Systems Laboratory. NOAA Tech. Rep. ERL FSL-4, 81 pp.
Stoffelen
,
A.
, and Coauthors
,
2005
:
The Atmospheric Dynamics Mission for global wind field measurement
.
Bull. Amer. Meteor. Soc.
,
86
,
73
87
, doi:.
Vaisala
,
2002
: Vaisala radiotheodolite RT20. Product brochure, 6 pp. [Available online at http://www.sankotsusho.co.jp/products/kisho/kousoukansoku/pdf/rt20.pdf.]
Wang
,
J.
,
J.
Bian
,
W. O.
Brown
,
H.
Cole
,
V.
Grubišic
, and
K.
Young
,
2009
:
Vertical air motion from T-REX radiosonde and dropsonde data
.
J. Atmos. Oceanic Technol.
,
26
,
928
942
, doi:.
Wang
,
L.
, and
M. A.
Geller
,
2003
:
Morphology of gravity-wave energy as observed from 4 years (1998–2001) of high vertical resolution U.S. radiosonde data
.
J. Geophys. Res.
,
108
, 4489, doi:.
WMO
,
2008a
: Measurement of meteorological variables. Guide to meteorological instruments and methods of observation, World Meteorological Organization CIMO Guide, 7th ed. WMO-8.
WMO
,
2008b
: Measurement of upper wind. Guide to meteorological instruments and methods of observation, World Meteorological Organization CIMO Guide, 7th ed. WMO-8, I.13-1–I-13.22.