Using WSR-88D Polarimetric Data to Identify Bird-Contaminated Doppler Velocities

Abstract

As an important part of Doppler velocity data quality control for radar data assimilation and other quantitative applications, an automated technique is developed to identify and remove contaminated velocities by birds, especially migrating birds. This technique builds upon the existing hydrometeor classification algorithm (HCA) for dual-polarimetric WSR-88D radars developed at the National Severe Storms Laboratory, and it performs two steps. In the first step, the fuzzy-logic method in the HCA is simplified and used to identify biological echoes (mainly from birds and insects). In the second step, another simple fuzzy logic method is developed to detect bird echoes among the biological echoes identified in the first step and thus remove bird-contaminated velocities. The membership functions used by the fuzzy logic method in the second step are extracted from normalized histograms of differential reflectivity and differential phase for birds and insects, respectively, while the normalized histograms are constructed by polarimetric data collected during the 2012 fall migrating season and sorted for bird and insects, respectively. The performance and effectiveness of the technique are demonstrated by real-data examples.

1. Introduction

Radar echoes from migrating birds can severely contaminate Doppler velocity measurements (Jungbluth et al. [1]; Gauthreaux et al. [2]). For meteorological applications, especially quantitative applications in radar data assimilation, it is necessary to remove bird-contaminated velocities using an automated identification technique. Such a technique was developed previously before the operational WSR-88D radars in the US were upgraded with dual-polarization capability (Zhang et al. [3]; Liu et al. [4]), but the technique could only crudely detect and remove bird-contaminated velocities volume-wise or tiltwise (for each sweep) because the usable input data from operational WSR-88D radars were limited to reflectivity and velocity measurements. A major drawback of this previous technique is that it rejects the entire volume (or tilt) of velocity observations even if only a fraction of the volume (or tilt) is contaminated by migrating birds. Now, almost all the operational WSR-88D radars in the US are upgraded with dual-polarization capability, which is essential for discriminating between meteorological and nonmeteorological scatterers and for distinguishing different hydrometeor types. In the hydrometeor classification algorithm (HCA) developed for polarimetric radars at the National Severe Storms Laboratory (NSSL) (Straka and Zrnić [5]; Zrnić and Ryzhkov [6]; Zrnić et al. [7]; Schuur et al. [8]; Ryzhkov et al. [9]; Park et al. [10]), birds and insects are not distinguished and labeled as biological scatterers entirely. This HCA is designed to remove biological scatterer contamination from meteorological scatterers for the purpose of quantitative precipitation estimation (QPE). It is thus timely and desirable to develop an improved technique that takes advantage of polarimetric measurements to provide a pixel-wise identification and removal of bird-contaminated velocities. As a continuation of the previous efforts of Zhang et al. [11, 12], this paper strives to develop such a desired technique by leveraging the NSSL existing HCA.

The most recent version of the HCA (Park et al. [10]) utilizes all available polarimetric variables to discern ten different classes of radar echoes with eight classes for various hydrometeors, one class for ground clutter (GC, including that due to anomalous propagation) and one class for biological scatterers (BS). This classification was designed mainly for using radar observations to improve the quantitative precipitation estimation. The GC class is useful for both reflectivity and velocity data quality controls, but the BS class is directly useful only for reflectivity data quality control. For radar velocity data quality control, it is necessary to differentiate insects from birds, because insects can be treated as passive tracers of air motions and thus are useful for wind measurements in most cases, while birds and especially migrating birds (mostly flying at night during the migrating seasons) can cause significant biases (to the order of 10 m s⁻¹) in radar measured Doppler velocities (Gauthreaux et al. [2]; Collins [13]; Bi et al. [14]). It has been known that the polarimetric signatures of insects and birds are different due to the shape and size differences between insects and birds (Riley [15]; Vaughn [16]; Zrnić and Ryzhkov [17]). In general, birds have higher differential phase and lower differential reflectivity than insects (Zhang et al. [11, 12]). These properties will be used in this paper to develop an automated technique for pixel-wise identification and removal of bird-contaminated velocities in two steps. In the first step, the fuzzy logic in the HCA is simplified to identify 3 classes: BS, GC, and MS, where MS represents meteorological scatterers and is essentially a combined category containing all meteorological hydrometeor categories from the original algorithm (see Section 2.1 of Schuur et al. [8]). In the second step, a new simple fuzzy logic method is developed and used to divide the BS class into two subclasses: (i) BS due to birds (including bats) especially migrating birds (BSb) and (ii) BS due to insects (BSi). The next section describes the simplified HCA developed for the first step. Section 3 presents the fuzzy logic method developed for the second step. The performance and effectiveness of the two-step method are demonstrated by real-data examples in Section 4. Conclusions follow in Section 5.

2. Simplified HCA for Identifying BS in the First Step

As explained in the introduction, the ten classes of radar echoes in the original HCA are consolidated into three classes using a simple fuzzy logic method that focuses on identifying BS, GC, and MS. Among the six polarimetric variables used by the original HCA, the simplified fuzzy logic method uses the following five variables: (i) the radar reflectivity 𝑍 (at horizontal polarization), (ii) the differential reflectivity 𝑍_DR, (iii) the cross-correlation coefficient 𝜌_hv between horizontally and vertically polarized radar returns, (iv) the texture parameter SD(𝑍) of 𝑍, and (v) the texture parameter SD(Φ_DP) of the differential phase Φ_DP. The specific differential phase 𝐾_DP is used in the original HCA for identifying different types of hydrometeors, so this variable is not used here.

As described in Park et al. [10], the parameter SD(𝑍) [or SD(Φ_DP)] characterizes the magnitude of the small-scale fluctuations of 𝑍 (or Φ_DP) along the radial. It is estimated by smoothing 𝑍 (or Φ_DP) data along the radial using a 1 km (or 2 km) running-average window, subtracting the smoothed 𝑍 (or Φ_DP) from the original values and calculating the root-mean-square (RMS) value of the residuals. Before applying the classification procedure, 𝑍, 𝑍_DR, and 𝜌_hv are smoothed along each radial using a 1 km averaging window for 𝑍 and a 2 km window for 𝑍_DR and 𝜌_hv. In addition, the Doppler velocity 𝑣_𝑟 is used to flag possible GC prior to the fuzzy logic classification procedure. In particular, if 𝑣_𝑟 < 1 m s⁻¹ (or ≥1 m s⁻¹) for a pixel, then this pixel will be (or not be) checked for GC in the fuzzy logic classification procedure. To account for attenuation in precipitation, the biases of 𝑍 and 𝑍_DR are estimated using heavily filtered Φ_DP for WSR-88D radars (Bringi et al. [18]): Δ𝑍^att(dBZ)=0.04Φ_DP^∘,Δ𝑍_DR^att(dB)=0.004Φ_DP^∘. (1)

From the above parameter values, the simplified fuzzy logic method computes the aggregation value 𝐴_𝑖 for the 𝑖th class of radar echoes at each pixel by 𝐴_𝑖=Σ_𝑗𝑊_𝑖𝑗𝑃^(𝑖)𝑉_𝑗Σ_𝑗𝑊_𝑖𝑗, (2) where 𝑃^(𝑖)(𝑉_𝑗) is the membership function that characterizes the distribution of the 𝑗th variable for the 𝑖th class and 𝑊_𝑖𝑗 is a weight between 0 and 1 assigned to the 𝑖th class and the 𝑗th variable. The classification decision is based on the maximal aggregation value. The formulation in (2) is similar to that in (3) of Park et al. [10] except that the confidence vector assigned to each variable is now simply set to 1.

Each membership function 𝑃^(𝑖)(𝑉_𝑗) in (2) has a trapezoidal shape described by four parameters: 𝑥₁, 𝑥₂, 𝑥₃, and 𝑥₄. This trapezoidal function increases linearly from 0 to 1 as 𝑉_𝑗 increases from 𝑥₁ to 𝑥₂, remains to be 1 until 𝑉_𝑗 increases 𝑥₃, and then decreases linearly from 1 to 0 as 𝑉_𝑗 further increases from 𝑥₃ to 𝑥₄ (see Figure 1 of Schuur et al. [8] or Park et al. [10]). The parameter values are listed in Table 1 for each of the 5×3 membership functions, where 𝑓₁ and 𝑓₂ are functions of 𝑍 (in dBZ) given by (4) of Park et al. [10] or originally by (7) and (8) of Schuur et al. [8] in the following forms: 𝑓₁=−0.50+2.50×10⁻³𝑍+7.50×10⁻⁴𝑍²,𝑓₂=0.08+3.64×10⁻²𝑍+3.57×10⁻⁴𝑍². (3)(4) Here, the five membership functions for GC are largely the same as those in Table  1 of Park et al. [10] or originally Table  1 of Schuur et al. [8] except for the following modifications according to our observations in the presence of GC: (i) the lower-limit value 𝑥₁ for 𝑍 is lowered from 15 to 5 dBZ, (ii) the lower-limit value 𝑥₁ for 𝑍_DR is increased from −4 to −3 dB , and (iii) the lower-shoulder value 𝑥₂ for 𝜌_hv is increased from 0.6 to 0.8 to reduce the possibility of misclassifying BS with 𝜌_hv between 0.5 to 0.8 as GC. The five membership functions for BS are mostly the same as those in Table 1 of Park et al. [10] or Table 1 of Schuur et al. [8], except that the upper-limit value 𝑥₄ for 𝜌_hv is increased from 0.83 to 1.01, because it has been found that 𝜌_hv occasionally reaches to 1.01 for BS in the operational polarimetric radar observations. The five membership functions for MS are mostly the same as those in Table 1 of Schuur et al. [8] except that the upper-limit value 𝑥₄ for 𝜌_hv is increased from 1.01 to 1.05 because incorrect values of 𝜌_hv>1 and even almost up to 1.05 have been seen from the polarimetric data so far collected in the presence of MS, and these incorrect values of 𝜌_hv(>1) are caused by excessive noise corrections (see (5) of Schuur et al. [8]).

Table 1 

Parameter values of the membership functions used by the simplified HCA for the three classes. Reflectivity-dependent parameters 𝑓₁ and 𝑓₂ are defined in (3) and (4), respectively.

The weights of the 5×3 membership functions are listed in Table 2. The weights of the five membership functions for GC are largely the same as those in Table 2 Park et al. [10] except for the following modifications: (i) the weight assigned to 𝑍 is increased from 0.2 to 0.4 to enhance and ensure the GC detection, (ii) the weight assigned to 𝜌_hv is decreased from 1.0 to 0.4, and (iii) the weight assigned to SD(𝑍) is fine tuned slightly from 0.6 to 0.5. The weights of the five membership functions for BS are the same as those in Table 2 of Park et al. [10]. The weights of the five membership functions for MS are the same as the maximum weights among the eight classes of hydrometeors in Table 2 of Park et al. [10].

Table 2 

Specified weights for the 5×3 membership functions in the simplified HCA.

3. Fuzzy-Logic Method for Identifying BSb in the Second Step

3.1. “Ground Truth” Data for Birds versus Insects

In the second step, the BS class in the HCA is divided into two subclasses: (i) BS due to birds especially migrating birds (BSb) and (ii) BS due mainly to insects (BSi). The goal is to develop a fuzzy logic method for further differentiating BSb from BSi for each BS pixel identified by the simplified HCA (described in Section 2). As reviewed in the introduction, insects and birds have different polarimetric signatures due to their differences in shape and size. In particular, birds have higher Φ_DP and lower 𝑍_DR than insects. These signatures are further quantified in Table 3 in comparison with those for MS. As an example, Figure 1 shows the color images of 𝑍_DR, Φ_DP, and 𝜌_hv from the KVNX 0.5° scan at 0102 UTC (local time 8 pm) on 27 October 2011. In this case, the following two types of pixels are readily identified by the simplified HCA: (i) BS pixels within the 120 km radial range around the radar and (ii) MS pixels far beyond 120 km radial range to the northwest of the radar. The identified BS pixels have relatively low 𝑍_DR, high Φ_DP, and low 𝜌_hv, so they should be further classified as BSb pixels. The identified MS pixels reveal a precipitation system to the northwest of the radar.

Table 3 

Typical ranges of 𝑍_DR, Φ_DP, and 𝜌_hv and height for BSb, BSi, and MS.

(a)

(b)

(c)

Figure 1 

Imagers of (a) 𝑍DR, (b) ΦDP, and (c) 𝜌hv from the KVNX 0.5° scan at 0102 UTC (local time 8 pm) on 27 October 2011. The spatial resolutions of the original data are 250 m in the radial direction and 0.5° in the azimuthal direction.

To develop a fuzzy-logic method to differentiate BSb from BSi, it is necessary to find proper membership functions that characterize the distributions of Φ_DP and 𝑍_DR for BSb versus BSi within the BS class based on the probability density function estimated from normalized histograms of Φ_DP and 𝑍_DR for the two subclasses. To construct these normalized histograms, it is necessary to collect large numbers of “ground truth” data for BSb versus BSi in the presence of BS, and this is done by the following selection procedure. First, it has long been recognized and well observed that BS echoes are caused dominantly by migrating birds during the nighttime but by insects during the daytime in the spring and fall migrating seasons, and the transition from birds to insects (or from insects back to birds) occurs around the sunrise time (or sunset time) almost every nonstormy day during a migrating season, especially over the southern great plain in the central United States (Zhang et al. [11]). Based on this well-observed fact (Dingle [19]; Drake and Gatehouse [20]), most BS pixels identified from nighttime (or daytime) radar scans during migrating seasons by the simplified HCA can be attributed to BSb (or BSi). Thus, we can sample midnight (or midday) BS data tilt-by-tilt and treat each sampled tilt of data as proxy “ground-truth” data for BSb (or BSi). Then, by visually examining color images of the sampled polarimetric data on each tilt, we can manually pick specific areas for BSb and BSi where the polarimetric signatures are distinct in Φ_DP and 𝑍_DR (with 𝜌_hv < 0.97) for BSb and BSi, respectively, according to Table 3. This completes the selection procedure.

As an example, Figures 2(a) and 2(b) show the 𝑍_DR--Φ_DP scattergrams for BS pixels identified by the simplified HCA from the KVNX 0.5° scans at 0706 UTC (local midnight) and 1901 UTC (local midday), respectively, on 17 October 2011. As shown, the midnight 𝑍_DR--Φ_DP scattergram in Figure 2(a) has relatively low 𝑍_DR and relatively high Φ_DP, and these are the expected characters for BSb. On the other hand, the midday 𝑍_DR--Φ_DP scattergram in Figure 2(b) has relatively high 𝑍_DR and relatively low Φ_DP, and these are the expected characters for BSi caused by insects. Thus, by plotting and examining their 𝑍_DR--Φ_DP scattergrams, we can further check and ensure the quality of the “ground truth” data.

(a)

(b)

Figure 2 

𝑍DR--ΦDP scattergrams for BS pixels identified by the simplified HCA from the KVNX 0.5° scans at (a) 0706 UTC (local midnight) and (b) 1901 UCT (local midday), respectively, on 17 October 2011.

3.2. Scattergrams and Normalized Histograms for Migrating Birds and Insects

To quantify and differentiate the polarimetric signatures between insects and birds, especially migrating birds, we started to monitor and sample real-time operational polarimetric observations from the operational KVNX radar since September 2011 (shortly after this radar was upgraded with dual-polarization capability, the first among the 159 operational WSR-88D radars in the US) and from the operational KICT radar since September 2012 (shortly after this radar was upgraded with dual-polarization capability). After polarimetric data are sampled from these radars during the fall and spring migrating seasons and processed by using the simplified HCA to identify BS pixels, “ground-truth” data are gathered for BSb (or BSi) through the two selection steps descried in Section 3.1. The “ground-truth” data gathered from KICT and KVNX radars for the 2012 fall migrating season have better qualities than those gathered early from KVNX radar for the 2011 fall migrating season (before KVNX radar was fully calibrated), so the scattergrams constructed from the “ground-truth” data gathered from KICT radar for the 2012 fall migrating season will be used to extract the membership functions for BSb and BSi in this subsection.

The red (or blue) dots in Figure 3 plot the 𝑍_DR--Φ_DP scattergram for BSb (or BSi) constructed by 1692359 (or 1385090) pixels of “ground-truth” data gathered from the KICT observations during the nighttime (or daytime) on 2 October 2012. Here, again as expected, the scattergram for BSb (or BSi) has relatively low (or high) 𝑍_DR and relatively high (or low) Φ_DP, and the main differences between the two scattergrams for BSb and BSi are similar to those showed in Figure 3 of Zrnić and Ryzhkov [17], although the core area of their scattergram for BSb is relatively narrow (roughly within −2 dB  ≤ 𝑍_DR ≤ 3 dB) and the core area of their scattergram for BSi is not bounded by 𝑍_DR = 8 dB. Note that the input Φ_DP is obtained by subtracting the radar system Φ_DP from the radar measured Φ_DP, but the radar system Φ_DP can be estimated and recorded only through a rainy event. The true system Φ_DP often drifts slowly away from the recorded value, but the latter will be still used as the true system Φ_DP until it is estimated and recorded through the next rainy event. The slow drift of the true system Φ_DP can cause uncertainties in the input Φ_DP (obtained by subtracting the recorded system Φ_DP from the radar measured Φ_DP). For the scattergrams in Figure 3, the subtracted system Φ_DP is 14.81°. Scattergrams constructed from “ground-truth” data gathered from the KICT observations for other days during the 2012 fall migrating season are similar to those in Figure 3 for BSb and BSi, respectively, but the core areas often shift to slightly different Φ_DP ranges due to the above explained uncertainties in the input Φ_DP.

Figure 3 

𝑍DR--ΦDP scattergrams for BSb (plotted by red dots) and BSi (plotted by blue dots) constructed by 1692359 and 1385090 pixels of “ground-truth” data gathered from the KICT observations during the nighttime and daytime, respectively, on 2 October 2012.

As we can see from Figure 3, the two scattergrams for BSb and BSi are not completely separated. Nevertheless, their overlapped area is small which may imply a small degree of true inseparability between BSb and BSi. The normalized histograms of 𝑍_DR (or Φ_DP) that estimate the probability density functions of 𝑍_DR (or Φ_DP) are extracted from the scattergrams in Figure 3 and plotted for BSb by the solid curve and for BSi by the dashed curve in Figures 4(a) (or 4(b)). As shown in Figure 4(a), BSb is dominant and the probability for BSi is almost zero when −4 dB ≤ 𝑍_DR < 4 dB, whereas BSi is dominant and the probability for BSb is almost zero when 𝑍_DR > 4 dB. As shown in Figure 4(b), BSb is dominant and the probability for BSi is nearly zero when Φ_DP > 50°, whereas BSi is dominant and the probability for BSb is much smaller than that for BSi when Φ_DP < 35°, so the situation becomes ambiguous (with no clear dominance in either BSb or BSi) only when 30° < Φ_DP < 50°. Normalized histograms of 𝜌_hv are also plotted and examined in pairs for BSb and BSi (not shown), and the dominant range of 𝜌_hv for BSb is found to be slightly lower than that for BSi. Since 𝜌_hv has a significantly overlapped range between BSb and BSi, it will not be used as an input variable for the fuzzy logic method in the next subsection. Normalized histograms of other variables, such as 𝑍, 𝑣_𝑟, SD(𝑍), SD(Φ_DP), and SD (𝑣_𝑟), exhibit largely overlapped ranges between BSb and BSi (not shown), so they are not very useful here. Thus, only 𝑍_DR and Φ_DP will be used as input variables for the fuzzy logic method in the next subsection.

(a)

(b)

Figure 4 

Normalized histograms of (a) 𝑍DR and (b) ΦDP exacted from the scattergrams in Figure 3 and plotted by the solid curve for BSb and by the dashed curve for BSi.

3.3. Fuzzy-Logic Method for Identifying Migrating Birds

The fuzzy logic formulation for identifying migrating birds is similar to that in (2) but is further simplified by using only 2 membership functions to characterize the distributions of the two variables 𝑍_DR and Φ_DP for the subclass BSb. Based on the normalized histogram of 𝑍_DR for BSb in Figure 4(a), the four parameters 𝑥₁, 𝑥₂, 𝑥₃, and 𝑥₄ are set to −5, −3, 2, and 4 dB, respectively, for the membership functions of 𝑍_DR for BSb. Based on the normalized histogram of Φ_DP for BSb in Figure 4(b), the four parameters 𝑥₁, 𝑥₂, 𝑥₃, and 𝑥₄ are set to 0°, 40°, 120^°, and 150°, respectively, for the membership function of Φ_DP for BSb. The weights are set to 1 and 0.8 for the membership functions of 𝑍_DR and Φ_DP, respectively. Because we only need to identify BSb (versus BSi) among BS, it is convenient and sufficient to compute a single aggregation value 𝐴_𝑖 in (2) for 𝑖 = BSb only and then check the computed 𝐴_BSb with a properly tuned threshold value (which is 0.3 in this case). Thus, if 𝐴_BSb is larger than this threshold value for a BS pixel, then this pixel is identified as a BSb pixel.

Since the method is developed as a part of the radar velocity data quality control for operational data assimilation applications (Liu et al. 2009 [21]), the goal is to detect and remove all pixels contaminated by migrating birds. To achieve this goal, the threshold is purposely tightened to a relatively low value (0.3) to enhance the detection of BSb. In this case, some BSi pixels may be also removed as they are incorrectly identified as “BSb” pixels due to the tightened threshold, and this is a price paid for the enhanced detection and removal of BSb to ensure that the processed and accepted velocity data are free of bird contaminations. Thus, the method is designed and tasked to detect and remove all true BSb pixels rather than to identify and retain all true BSi pixels, and this task is implemented effectively by considering only the subclass BSb in the fuzzy logic formulation.

4. Effectiveness of the Method

The two-step method developed in this paper has been tested with polarimetric data collected from the operational KVNX and KICT radars during the 2011 fall and 2012 spring and fall migrating seasons. The results indicate that the method is as effective as, but more efficient than, the original HCA for identifying GC and BS in the first step. The method is also demonstrated to be effective and efficient for identifying BSb in the second step. Examples are presented below to show the effectiveness of this two-step method.

4.1. Effectiveness of the Simplified HCA

As mentioned in Section 2, the simplified HCA is designed to identify and discriminate BS from GC and MS. An example is given in this subsection to illustrate that this simplified HCA is effective for this designed purpose in comparison with the original HCA. Figures 5(a) and 5(b) show the images of reflectivity 𝑍 and radial velocity 𝑣_𝑟, respectively, at the lowest tilt (0.5°) scanned by the operational KVNX radar at 013037 UTC on 28 October 2011. The local time was 8:30 pm. As shown, the radar echoes were weak over most of the area within 100 km radial range. 𝑍_DR and Φ_DP images (not show) indicate that these weak echoes are BS from migrating birds. They were also GC identified empirically (with |𝑣_𝑟| < 1 m s⁻¹ as shown in Figure 5(b)) around 𝑟 = 50 km   to the south-southeast and west-northwest of the radar. Around 𝑟 = 120 km to the north of radar, there was an isolated area of weak echo that lasted only for a few volume scans (not shown) around sunrise and sunset, and these weak echoes were BS of local birds, as identified by human expertise. The scattered areas with relatively high reflectivity (around 20 dBZ) outside 𝑟 = 150 km to the southeast of radar are clearly MS, again as identified with human expertise. Figure 5(c) shows the classification results from the simplified HCA. Here, the BS, GC, and MS are colored in cyan, gray, and green, respectively. The results are consistent with those identified above based on human expertise. The classification results from the original HCA are plotted in Figure 5(d), where the BS and GC are again in cyan and gray, respectively, and the eight different classes of MS are in other different colors. By comparing Figure 5(c) with Figure 5(d), we can see that the BS and GC areas identified by the simplified HCA cover the BS and GC areas identified by the original HCA, respectively, and the MS areas identified by the simplified HCA cover the areas of hydrometeors identified by the original HCA. Thus, the BS pixels identified by the simplified HCA in the first step can be used to further identify BSb pixels in the second step.

(a)

(b)

(c)

(d)

Figure 5 

Images of (a) 𝑍, (b) 𝑣𝑟, (c) classification results produced by the simplified HCA, and (d) classification results produced by the original HCA from the KVNX 0.5° scan at 013150 UTC on 28 October 2011. The 10 classes in the original HCA are: (1) ground clutter including that due to anomalous propagation (GC/AP); (2) biological scatterers (BS); (3) dry aggregated snow (DS); (4) wet snow (WS); (5) crystals of various orientations (CR); (6) graupel (GR); (7) “big drops” (BD); (8) light and moderate rain (RA); (9) heavy rain (HR); and (10) a mixture of rain and hail (RH). The spatial resolutions of the original data are the same as in Figure 1.

4.2. Effectiveness of the Fuzz-Logic Method in the Second Step

4.2.1. Effectiveness for Identifying BSb

After BS pixels are identified by the simplified HCA in the first step, the fuzz-logic method described in Section 3.3 is applied to discriminate BSb from BSi among the BS pixels. An example is given in this subsection to illustrate the effectiveness of the method for identifying and removing BSb contaminations. Figures 6(a)–6(e) show the images of reflectivity 𝑍, radial velocity 𝑣_𝑟, cross-correlation coefficient 𝜌_hv, differential reflectivity 𝑍_DR, and differential phase Φ_DP, respectively, at the lowest tilt (0.5°) scanned by the operational KICT radar at 090142 UTC on 30 September 2012. The local time was around 4 am during the fall migrating season, and the sky over the 150 km  radial range from the KICT radar was covered by elongated rapidly-moving patches of thin clouds aloft with no precipitation or, at least, no significant precipitation based on GOES-14 infrared observations as shown in Figure 7. The presence of migrating birds is apparent and prevailing as judged by human expertise from the radar reflectivity image in Figure 6(a) and the GOES-14 infrared image in Figure 7. In particular, as shown in Figure 6(a), the gray colored reflectivity (<18 dBZ) over the broad area within 150 km radial range, which is identified as BS by the simplified HCA in the first step, was caused by migrating birds, whereas the green colored reflectivity (>18 dBZ) beyond 150 km  to the south of radar, which is identified as MS in the first step, was caused by precipitation. Figure 6(c) shows that the 𝜌_hv values (mostly around 1) in the MS area are larger than the 𝜌_hv values (mostly below 0.7) in the broad BSb area. The 𝑍_DR field in Figure 6(d) and Φ_DP field in Figure 6(e) are smoother and more uniform in the MS area than in the BSb area. Moreover, the 𝑍_DR values are mostly below 3 dB and the Φ_DP values are larger than 60° in the broad BSb area. These polarimetric features further confirm that the BS pixels identified in the first step are indeed BSb pixels rather than BSi pixels. After MS and BS pixels are identified by the simplified HCA in the first step, BSb pixels are further identified from the BS pixels by the fuzz-logic method in the second step. The final results are shown in Figure 6(f), where the MS pixels are in green and the BSb pixels are in yellow. These results are consistent with those identified above by human expertise. This exemplifies the effectiveness of the method for identifying BSb and removing BSb-contaminated velocities.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Images of (a) 𝑍, (b) 𝑣𝑟, (c) 𝜌hv, (d) 𝑍DR, (e) ΦDP, and (f) final classification results produced in the second step from the KICT 0.5° scan at 090142 UTC (local time around 4 am) on 30 September 2012. The spatial resolutions of the original data are the same as in Figure 1.

Figure 7 

Infrared images from GOES-14 at 084500 UTC (upper panel) and 091500 UTC (lower panel) on 30 September 2012 over Kansas State and its neighboring states. The state boundaries are plotted in purple. The three cyan circles plot the 50, 100, and 200 km  range rings cocentered at the KICT radar site (shown by the cyan dot). The dark areas show relatively warm brightness temperature remotely sensed from the ground, indicating clear sky over each dark area. The green areas highlight cloud tops colder than 240°K (≈−30°C around 7 km  height).

From the classification results in Figure 6(f), we can also see that the method is not perfect and its performance is affected by noises in polarimetric measurements. For example, ideally, all the pixels should be identified as MS in the green-dominated areas beyond the 175 km  radial range to the south and southeast of the radar in Figure 6(f), but the method identifies 37562 MS pixels among the total of 39452 pixels in these green-dominated area, so the hit rate is 95.21% slightly below the perfect 100%. The remaining pixels are 4.79% of the total, and they are identified as GC, BSi, and BSb in 0.68%, 0.25%, and 3.86%, respectively. Thus, the correctly retained (MS and BSi) pixels are 95.46% of the total, and the incorrectly rejected (GC and BSb) pixels are merely 4.54% of the total. Similarly, all the pixels should be identified as BSb in the yellow-dominated area within the 150 km  radial range in Figure 6(f), but the method identifies 372021 BSb pixels among the total of 392649 pixels in this yellow-dominated area, so the hit rate is 94.75%. The remaining pixels are 5.25% of the total, and they are identified as GC, BSi, and MS in 0.56%, 1.95%, and 2.74%, respectively. Thus, the correctly rejected (GC and BSb) pixels are 95.31% of the total, and the incorrectly retained (MS and BSi) pixels are merely 4.69% of the total. These retained (MS and BSi) pixels are scattered and mostly isolated from each other, so they can be easily removed by a simple continuity check after the classification.

Figure 8 shows the vertical profiles of absolute value and azimuthal angle for the following two types of horizontal velocities: (i) the horizontal velocity of air plus birds (plotted by the solid profiles) produced by the velocity azimuth display (VAD) analysis (which is a by-product of the VAD-based dealiasing of Xu et al. [22]) from the same volume of velocity data as that displayed (on 0.5° tilt only) in Figure 6(b) but within the 150 km radial range to cover BSb pixels only and (ii) the horizontal velocity of air (plotted by the dashed profiles) above the KICT radar site at 090142 UTC interpolated from the hourly analyzed wind fields produced by the operational Rapid Refresh (RAP) model (see http://rapidrefresh.noaa.gov/) on 30 September 2012. Subtracting the second velocity from the first gives the differential velocity caused by migrating birds in the BSb area. The absolute value and direction of the differential velocity are plotted by the dotted vertical profiles in Figure 8. As shown, when the height increases from 250 to 800 m, the absolute value of the differential velocity increases from 4.5 to 6.3 m s⁻¹ and the azimuthal angle changes slightly from 190° to 176° (that is, within ±10° of the southward direction). This result is consistent with the well-recognized fact that migrating birds fly southward at speeds of 5~10 m s⁻¹ nighttime during the fall migrating season over the southern great plain in central United States. The differential velocities diagnosed in Figure 8 are also in the range of the previously documented VAD wind biases caused by migrating birds (Jungbluth et al. [1]; Gauthreaux et al. [2]; Collins [13]). Clearly, the differential velocities caused by migrating birds are not small. It is thus necessary to identify BSb pixels and remove them from velocity data, and this task can be performed effectively by the method in the second step.

(a)

(b)

Figure 8 

Vertical profiles of (a) absolute value (m s−1) and (b) azimuthal angle (clockwise from the north) for two types of horizontal velocities: (i) the horizontal velocity of air plus birds (plotted by the solid profiles) produced by the VAD analysis from the KICT-scanned volumetric velocity data at 090142 UTC but within the 150 km radial range to cover BSb pixels only and (ii) the horizontal velocity of air (plotted by the dashed profiles) interpolated from RAP hourly wind analyses to the KICT radar location at 090142 UTC on 30 September 2012. The absolute value and azimuthal angle of the vector difference of the above two velocities (the first minus the second velocity) are plotted by the dotted vertical profiles in panels (a) and (b), respectively.

4.2.2. Effectiveness for Identifying BSi

The method is also effective for identifying and retaining BSi pixels in the second step, and this is illustrated by the example in Figure 9 where the images were from the daytime KICT 0.5° scan at 215544 UTC (4 pm local time) on 30 September 2012. As shown in Figure 9(a), most pixels within the 150 km radial range are gray colored with 𝑍 < 15 dBZ. As these gray-colored low-reflectivity pixels are already classified as BS by the simplified HCA in the first step, they can be further identified as BSi pixels judged by human expertise based on the time of the day (4 pm local time). Furthermore, as we can see from Figures 9(c)–9(e), on these low-reflectivity BS pixels, the 𝜌_hv values are mostly below 0.97 (but occasionally become unrealistically larger than 1 and even reach 1.05 for the reason explained in Section 2), the 𝑍_DR values are mostly above 6 dB, and the Φ_DP values are mostly below 60°. These polorimetric features further confirm that the low-reflectivity BS pixels can be identified as BSi pixels according to Table 3. From Figure 9(a), we can also see scattered areas of green-colored pixels (with 𝑍 > 15 dBZ) near the 100 km range to the north of radar, around the 70 km range to the northeast of radar, and inside the 50 km range to the northwest of the radar. These scattered areas (𝑍 > 15 dBZ) are also already classified as BS by the simplified HCA in the first step. In these scattered areas (with 𝑍 > 15 dBZ), again as shown in Figures 9(b)–9(e), the observed radial velocities are highly nonuniform with large-amplitude irregular fluctuations, the 𝜌_hv values are mostly below 0.7, the 𝑍_DR values are mostly below 3 dB, and the Φ_DP values are mostly above 60°. These features indicate that the BS pixels in the scattered areas of 𝑍 > 15 dBZ can be further identified as BSb pixels according to Table 3, but these BSb pixels indicate the presence of local nesting birds rather than migrating birds as identified by human expertise based on the local time (around 4 pm) and the highly fluctuated radial velocities in these scattered areas (𝑍 > 15 dBZ). The BSb and BSi areas identified above by human expertise are also well identified automatically by the method in the second step as shown in Figure 9(f) by the yellow and cyan arrears, respectively.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9 

As in Figure 6 but from the KICT 0.5° scan at 215542 UTC (local time around 4 pm) on 30 September 2012.

Figure 10 shows the vertical profiles of absolute value and azimuthal angle for (i) the VAD velocity (plotted by the solid profiles) computed from the same volume of velocity data as that displayed (on 0.5° tilt only) in Figure 9(b) but with all the BSb pixels removed, (ii) the interpolated velocity from RAP hourly wind analyses (plotted by the dashed profiles) at 215542 UTC on 30 September 2012, and (iii) the differential velocity (plotted by the dotted vertical profiles) obtained by subtracting the second velocity from the first. As shown, the differential velocity is small, and its absolution value decreases from 3.4 to 1.6 m s⁻¹ and its azimuthal angle is nearly constant around 150~155° when the height increases from 250 to 800 m . These diagnosed differential velocities could be due to the RAP hourly wind analysis error plus insects’ flying speeds (Drake and Gatehouse [20]). This further confirms that BSi pixels are correctly identified by the method in the second step. As the differential velocities caused by insects are small, the observed radial velocities on BSi pixels are retained by the method for data assimilation applications.

(a)

(b)

Figure 10 

As in Figure 8 except that the VAD velocity (plotted by the solid profiles) is produced from the KICT-scanned volumetric velocity data at 215542 UTC with all the BSb pixels removed and the interpolated velocity from RAP hourly wind analyses (plotted by the dashed profiles) is also at 215542 UTC on 28 October 2011.

5. Conclusions

In this paper, a two-step method is developed to identify and remove contaminated velocities by birds, especially migrating birds, in addition to those contaminated by ground clutter (GC, including that due to anomalous propagation). In the first step, the existing hydrometeor classification algorithm (HCA) developed for polarimetric radars at the NSSL (Zrnić et al. [7]; Schuur et al. [8]; Ryzhkov et al. [9]; Park et al. [10]) is simplified to identify three classes of radar echoes: (i) GC, (ii) biological scatterers (BS), and (iii) meteorological scatterers (MS). In the second step, a fuzzy-logic method is developed and used to further identify scatterers of birds (including bats) especially migrating birds (BSb) versus scatterers of insects (BSi) among BS. The simplified HCA in the first step uses five input polarimetric variables to identify BS, GC, and MS, and the associated membership functions (listed in Table 1) are condensed from those used in the original HCA. Each BS pixel identified in the first step is used as a starting point to further identify BSb versus BSi on that BS pixel in the second step. The fuzzy-logic method in the second step uses two input polarimetric variables, that is, 𝑍_DR and Φ_DP, and the associated membership functions are extracted from normalized histograms of 𝑍_DR and Φ_DP that estimate the probability density functions of 𝑍_DR and Φ_DP, respectively, while the normalized histograms are constructed from “ground truth” data selected empirically by human expertise (as described and exemplified in Section 3.1).

The method has been tested with polarimetric data collected from the operational KVNX and KICT radars during the 2011 fall and 2012 spring and fall migrating seasons. The simplified HCA in the first step is found to be as effective as the original HCA, and the effectiveness is shown by the example in Section 4.1. This simplified HCA has been incorporated into the operational radar reflectivity data quality control package (Liu et al. [21]) for radar data assimilation application at the national centers for environmental prediction (NCEP). The fuzzy-logic method in the second step is also found to be effective in further identifying BSb versus BSi for each BS pixel identified in the first step, and the effectiveness is shown by the examples in Section 4.2. The two-step method has been continuously applied to real-time polarimetric data from the operational KVNX and KICT radars with the classification results verified indirectly by comparing the interpolated RAP wind analyses with the two types of VAD velocities produced from volumetric velocity data (i) on BSb pixels only and (ii) on non-BSb pixels only (as illustrated by the examples in Figures 8 and 10). Through this real-time run, more and better “ground truth” data will be accumulated and used to further improve the method in the second step. The improved method will be incorporated into the operational radar velocity data quality control package (Liu et al. [21]) for radar data assimilation application at NCEP.

Acknowledgments

The authors are thankful to Bob Rabin at NSSL for kindly providing the GOES-14 infrared images in Figure 7, Terry Schuur at NSSL, and anonymous reviewers for their comments and suggestions that improved the representation of the results. The research was supported by the ONR Grant N000141010778 to the University of Oklahoma (OU). Funding was also provided by NOAA/OAR under NOAA-OU Cooperative Agreement NA17RJ1227 U.S. Department of Commerce, by the China’s National Natural Science Foundation Grant 41175038, by the National Key Program for Developing Basic Sciences Grant 2012CB417202, and by the Chinese Academy of Meteorological Sciences Basic Scientific and Operational Project “Observation and retrieval methods of micro-physics and dynamic parameters of cloud and precipitation with multi-wavelength Remote Sensing.”