Abstract

The location and intensity of mesocyclone circulations can be tracked in real time by accumulating azimuthal shear values over time at every location of a uniform spatial grid. Azimuthal shear at low (0–3 km AGL) and midlevels (3–6 km AGL) of the atmosphere is computed in a noise-tolerant manner by fitting the Doppler velocity observations in the neighborhood of a pulse volume to a plane and finding the slope of that plane. Rotation tracks created in this manner are contaminated by nonmeteorological signatures caused by poor velocity dealiasing, ground clutter, radar test patterns, and spurious shear values. To improve the quality of these fields for real-time use and for an accumulated multiyear climatology, new dealiasing strategies, data thresholding, and multiple hypothesis tracking (MHT) techniques have been implemented. These techniques remove nearly all nonmeteorological contaminants, resulting in much clearer rotation tracks that appear to match mesocyclone paths and intensities closely.

1. Introduction

a. Motivation

To depict the paths and intensities of mesocyclone circulations as seen by radar, the National Severe Storms Laboratory (NSSL) creates products called rotation tracks. These rotation track fields are created by calculating the azimuthal shear fields (the azimuthal derivative of radial velocity) for each radar, merging the azimuthal shear data from multiple radars onto Cartesian grids and then accumulating the maximum values in those gridded fields over time onto one accumulated grid. An example is shown in Fig. 1.

Fig. 1.

Low-level (0–3 km AGL) rotation tracks from the 27 Apr 2011 tornado outbreak. Swaths of high maximum azimuthal shear approximate the movement and strength of mesocyclone circulations within the supercells that moved from southwest to northeast across the states of AL, MS, and TN between 1600 UTC 27 Apr and 0000 UTC 28 Apr.

Fig. 1.

Low-level (0–3 km AGL) rotation tracks from the 27 Apr 2011 tornado outbreak. Swaths of high maximum azimuthal shear approximate the movement and strength of mesocyclone circulations within the supercells that moved from southwest to northeast across the states of AL, MS, and TN between 1600 UTC 27 Apr and 0000 UTC 28 Apr.

These tracks can help alleviate some of the difficulty in interpreting and analyzing velocity fields. They can provide information about the spatial extent and strength of mesocyclone signatures over time and can be quite useful for conducting poststorm damage surveys. The American Red Cross of central Oklahoma uses NSSL rotation track products plotted on maps to determine where to deliver assistance after tornado events and the best routes to get there. After the 24 May 2011 tornado outbreak in central Oklahoma, the use of rotation tracks helped to significantly reduce the disaster assessment time (NOAA/NSSL 2011). These fields are also useful in real time as guidance for forecasters and have enormous data-mining potential. However, they are plagued by nonmeteorological signatures caused by poor velocity dealiasing, ground clutter, radar test patterns, and spurious shear values. As seen in Fig. 2, these artifacts can make the tracks difficult, if not impossible, to interpret meaningfully.

Fig. 2.

Low-level (0–3 km AGL) rotation tracks across VA, NC, and SC from the 16 Apr 2011 tornado outbreak generated using the default WSR-88D velocity dealiasing technique and without any quality control techniques. The spikes in the high azimuthal shear values are caused by poor velocity dealiasing along radials and can make data interpretation difficult, if not impossible, in some areas.

Fig. 2.

Low-level (0–3 km AGL) rotation tracks across VA, NC, and SC from the 16 Apr 2011 tornado outbreak generated using the default WSR-88D velocity dealiasing technique and without any quality control techniques. The spikes in the high azimuthal shear values are caused by poor velocity dealiasing along radials and can make data interpretation difficult, if not impossible, in some areas.

One of the goals of the Multi-Year Reanalysis of Remotely Sensed Storms (MYRORSS) project, a cooperative effort between National Oceanic and Atmospheric Administration's (NOAA) NSSL and the National Climatic Data Center (NCDC), is to create a conterminous United States (CONUS)-wide climatology of low- and midlevel rotation track fields for the lifetime of the Weather Surveillance Radar-1988 Doppler (WSR-88D) network.

Cintineo et al. (2011) developed an automated system to process level II radar data (Crum and Alberty 1993) at NSSL using a multiple-machine framework and the Warning Decision Support System—Integrated Information (WDSS-II; Lakshmanan et al. 2007b) suite of programs to process and quality control the data. A preliminary hail climatology using maximum estimated size of hail (MESH) grids was also created (Cintineo et al. 2012). Here, we extend the processing to velocity-based products, specifically to azimuthal shear accumulations.

The CONUS-wide rotation track climatology will provide an incredibly rich dataset with numerous potential climatological and severe weather applications. Track lengths, intensities, and other characteristics could be analyzed by geographical region and time of year. Potential relationships could also be discovered between rotation tracks and environmental parameters like convective available potential energy (CAPE) and storm-relative helicity (SRH). After correlating maximum azimuthal shear and maximum updraft helicity, rotation tracks could also be used as verification for high-resolution model-simulated maximum updraft helicity tracks like the ones discussed by Kain et al. (2010) and Clark et al. (2012, 2013).

The purpose of this paper is to discuss the special velocity dealiasing techniques, data thresholds, and multiple hypothesis tracking (MHT) techniques developed to isolate the rotation tracks in real time and for the MYRORSS climatology. Detailed explanations of each quality control effort will be given and the specific impacts of each step will be shown in example cases.

b. Background

Modern Doppler radars have the ability to provide high-resolution space and time measurements of storms that allow for the detection of mesocyclone-scale circulations. Couplets in radial velocity fields as well as hook echo signatures in reflectivity fields have been used to identify these circulations in the past, but methods of identifying mesocyclone or tornado circulations reliant solely on hook echo signatures have proven unreliable (Forbes 1981; Mitchell et al. 1998). Methods using radial velocity signatures such as the tornado detection algorithm (TDA; Mitchell et al. 1998) and the NSSL mesocyclone detection algorithm (MDA; Stumpf et al. 1998) have been more successful.

The TDA currently used with the WSR-88D system relies on high “gate-to-gate velocity difference” values to identify potentially tornadic circulations (Mitchell et al. 1998). Although termed a tornado detector, the algorithm identifies tornadic vortex signatures (which may or may not be associated with tornadoes) that are typically larger than a tornado owing to radar sampling resolution (e.g., Brown et al. 1978). The gate-to-gate difference represents the difference between velocity values at constant range from the radar between adjacent azimuths. These values can be affected adversely by the azimuthal offset of the radar beam center from the vortex, noisy data, and velocity aliasing (Wood and Brown 1997). Additionally, because the radar beam is much broader at far ranges than when near the radar-observed velocity peaks within vortices decrease in magnitude, allowing some vortices to be overlooked by the algorithm.

Liu et al. (2007) proposed a wavelet analysis technique to help mitigate these issues. The method examines region-to-region radial wind shears at a number of different spatial scales to more accurately determine the amount of shear present, reducing the number of false tornado detections.

Rather than using the gate-to-gate velocity difference, “peak to peak” methods of calculating rotational shear from Doppler radial velocity data or the wavelet analysis technique to detect a vortex, a two-dimensional (2D), local, linear least squares derivatives (LLSD) method can be used to reduce the impact of noise. Elmore et al. (1994) proposed this method of estimating the derivatives of radial velocity values by fitting a plane to the velocity field and finding its slope. The vertical vorticity field is estimated by the azimuthal derivative of the radial velocity field and is given by

 
formula

where is the radial velocity, is the coordinate in the azimuthal direction, is the arc length from the center point of the calculation to the point , is the radial velocity at point , and is the beamwidth at a given range. The coordinate is in the radial direction and is in the azimuthal direction. In addition, is a positive weight function that we set to 1 after determining that Cressman weight functions, among others, generated very small differences. The summation is performed over range gates in the neighborhood of the starting point of the calculation. This calculation of azimuthal shear is, here on, referred to as the LLSD method.

Smith and Elmore (2004) applied the LLSD calculation to simulated and observed circulations by first passing the velocity data through a 3 × 3 median filter (Lakshmanan 2012) to reduce speckle noise and then applying Eq. (1) to the filtered velocity data to estimate the azimuthal shear. The physical size of the neighborhood used in the calculation is held constant such that fewer radials are used in the calculation at far ranges from the radar. Typical sizes of radial and azimuthal neighborhoods are 750 and 2500 m, respectively.

The LLSD calculation helps to remove some of the dependence on radar location involved in rotation detection and also allows circulation signatures to be viewed in 3D space or as input to multisensor applications. It was shown by Smith and Elmore (2004) that LLSD shear values were reasonable estimations of actual shear values in simulated Rankine vortices when sampled by a theoretical WSR-88D radar (Brown et al. 2002) out to a range of approximately 140 km. The variance of these values was also much smaller compared to the peak-to-peak shear.

The use of the median filter in the LLSD shear technique can be a disadvantage, however. Whereas areas with large velocity gradients are preserved, this filter can also smooth out peaks in the velocity field. Although the median filter is beneficial when these peaks are associated with noisy data, the filter decreases the magnitude of peak velocities in mesocyclone signatures and completely eliminates tornado signatures in nearly all cases. For this reason, LLSD shear values may underestimate the actual azimuthal shear of circulations, especially for small circulations (Mitchell and Elmore 1998).

Defining a neighborhood size in the LLSD technique also forces a trade-off between spatial resolution and noise resistance. Smaller neighborhoods are more strongly affected by noise, whereas larger neighborhoods tend to underestimate actual shear values of circulations. The spatial scale of the LLSD calculation makes it most useful for detecting large mesocyclone-scale circulations; however, the increased noise resistance makes identification of small circulations more difficult.

Rotation tracks help to visualize the movement of mesocyclone circulations (or occasionally circulations associated with nearby large tornadoes) over time as seen by radar. To produce these tracks, radial velocity data are first dealiased using the default WSR-88D dealiasing algorithm (Eilts and Smith 1990). Next, a quality control neural network (Lakshmanan et al. 2007a) is used to remove echoes in the reflectivity field produced by biological targets, anomalous propagation, ground clutter, and test or interference patterns. The algorithm successfully removes nearly all nonmeteorological signatures from the reflectivity fields examined.

This quality-controlled reflectivity field, hereafter reflectivityQC, and the radial velocity field are then employed by the shear estimation algorithm of Smith and Elmore (2004) to compute the azimuthal shear (see example in Fig. 3). The 2D maximum azimuthal shear fields within low (0–3 km) and midlevel (3–6 km) layers above ground level (AGL) are also calculated using digital elevation model (DEM) data to determine the height of each point above the ground.

Fig. 3.

Radial velocity fields and the corresponding azimuthal shear fields from Jackson, MS (KDGX), at 2151 UTC 27 Apr 2011 created using different dealiasing strategies. Radial velocity dealiased using (a) the LED algorithm, (c) the LED algorithm with an environmental wind field from an input 20-km RUC point sounding, and (e) the Jing and Wiener (1993) technique with an environmental wind field from an input 20-km RUC point sounding. (b),(d),(f) Azimuthal shear fields associated with (a),(c), and (e), respectively.

Fig. 3.

Radial velocity fields and the corresponding azimuthal shear fields from Jackson, MS (KDGX), at 2151 UTC 27 Apr 2011 created using different dealiasing strategies. Radial velocity dealiased using (a) the LED algorithm, (c) the LED algorithm with an environmental wind field from an input 20-km RUC point sounding, and (e) the Jing and Wiener (1993) technique with an environmental wind field from an input 20-km RUC point sounding. (b),(d),(f) Azimuthal shear fields associated with (a),(c), and (e), respectively.

These single-radar 2D maximum azimuthal shear fields are then merged into a Cartesian multiradar grid using the intelligent agent formulation of Lakshmanan et al. (2006), accounting for varying radar beam geometry with range, vertical gaps between radar scans, and other issues. The maximum value of each pixel in the merged multiradar grid over a time interval (typically 60–120 min) is then used to produce the swaths of merged maximum azimuthal shear known as rotation tracks. Missing data in the rotation track fields (denoted by MD in the color bar) correspond to any data below the signal to noise threshold for the radar. Ideally, this should be shown as zero shear, not missing data. Figure 4 provides a flowchart showing the algorithms and fields used to create rotation tracks.

Fig. 4.

Flowchart showing how rotation track products are created. Gray boxes with dashed lines represent algorithms and white boxes with solid lines represent data fields.

Fig. 4.

Flowchart showing how rotation track products are created. Gray boxes with dashed lines represent algorithms and white boxes with solid lines represent data fields.

2. Methods

a. 2D velocity dealiasing

Because of the relationship between radar wavelength and the pulse repetition frequency a radar correctly measures the radial velocity given that it is in the range of where

 
formula

Here, is the Nyquist velocity and the true velocity is

 
formula

where is the measured velocity, is an unknown integer including zero, and must satisfy . Velocity dealiasing is the process of determining the correct value of n to successfully recover . In the cases where cannot be successfully recovered, the velocity is still aliased and can usually be identified in radial velocity fields by abrupt changes in values between neighboring measurements. Most first-generation dealiasing algorithms (e.g., Ray and Ziegler 1977; Bargen and Brown 1980) were one-dimensional and detected abrupt changes between single radials. For this reason, they were quite sensitive to noisy, incorrect data. Strong shear zones in these radial velocity fields sometimes cannot be dealiased without data from multiple dimensions. Merritt (1984), Boren et al. (1986), and Bergen and Albers (1988) took more sophisticated approaches and used 2D velocity data to dealias. These methods were costly in terms of computation time, however.

The local environment dealiasing (LED) algorithm (Eilts and Smith 1990) is the method currently used for WSR-88D data in real time. The scheme applies radial continuity constraints to remove local aliasing errors and azimuthal continuity checks to mitigate error. It also incorporates radial averages to determine n [see Eq. (3)] when continuity thresholds are not met. Each radial is processed individually and compared against the previously dealiased radials, allowing the algorithm to use less memory and process faster than other 2D algorithms. It can also ingest a vertical wind profile from an environmental sounding to produce initial values for each elevation scan and for isolated echoes. It is an efficient algorithm that begins by using simple checks and only moves on to more sophisticated techniques if needed. This approach performed very well on the cases of severe aliasing presented in Eilts and Smith (1990), but performed poorly when qualitatively examined in many of the tornadic cases examined for this study. Ingesting the vertical wind field from a 20-km Rapid Update Cycle model (RUC; Benjamin et al. 2004) point sounding at each radar site to use as an environmental estimate into the LED algorithm improved the dealiasing to some extent.

For this study, a sophisticated 2D dealiasing technique described by Jing and Wiener (1993) was implemented. The algorithm solves a linear system of equations that minimizes gate-to-gate shear in each isolated 2D region. Through using aliasing-induced discontinuity information, the correction values for all gates are found by solving a 2D least-mean-squares problem. Instead of making dealiasing decisions for each gate based on its neighbors, which can be subject to scattered incorrect data, this approach avoids local expansion of errors by attempting to find all dealiased values for a given dataset. Vertical profiles of horizontal wind data from the 20-km RUC point soundings were used as environmental wind estimates at the grid point nearest to each radar site. The calculated average is minimized by incrementing equally over the entire echo. The average local wind observed by radar is assumed to be less than .

A smooth environmental wind field with weak shear is assumed. This can be a poor assumption in isolated areas of strong wind shear associated with mesocyclones or microbursts. In these cases, relatively short falsely aliased border segments are detected and can typically be used to dealiase the field correctly. In more elongated regions of shear associated with strong gust fronts, for example, incorrect dealiasing is more likely. Example dealiased radial velocity fields and their corresponding azimuthal shear fields using both the LED and Jing and Wiener (1993) methods are shown in Fig. 3. Preliminary results from a study now under way to determine which velocity dealiasing method performs best on a set of case studies indicate that this 2D technique is more accurate than the LED technique.

b. Reflectivity quality control

As mentioned in the introduction, a quality control neural network (Lakshmanan et al. 2007b; Lakshmanan et al. 2010) is used to remove nonmeteorological echoes from the reflectivity field. The algorithm combines various measures from both the past literature (e.g., Steiner and Smith 2002; Kessinger et al. 2003; Fulton et al. 1998) and new measures to discriminate between precipitating and nonprecipitating echoes in the reflectivity data. A region-by-region classification is performed rather than a gate-by-gate basis. In addition, clear-air echoes from biological contamination are identified and removed using a two-stage intelligent machine algorithm while retaining echoes that correspond to precipitation (Lakshmanan et al. 2010).

c. Additions to LLSD shear algorithm

In addition to calculating the azimuthal shear fields as described earlier in this paper, extra operations have been added to the LLSD algorithm. These additions are discussed in the following sections.

1) Azimuthal shear range correction

A new azimuthal shear range-correction (Newman et al. 2013) algorithm is applied to the field in an effort to correct for range degradation resulting from radar beam widening. A multiple linear regression technique was used to create the equations based on comparisons between observed shear values and those computed using simulated Rankine vortices. First, the algorithm identifies significant circulations using reflectivity and LLSD shear criteria. To avoid applying the range correction to regions of noise in the shear fields collocated with low reflectivity values, only circulations in regions with reflectivity values greater than 20 dBZ and LLSD shear values exceeding 0.005 s−1 are identified. The peak-to-peak velocity differences and shear diameters of circulations satisfying the reflectivity and shear criteria are calculated next. A median filter is applied to the shear diameter field to provide potentially more accurate estimates of circulation size when circulations are larger than tornadic vortex signatures (TVSs; Brown et al. 1978). Then, new azimuthal shear values for each pixel in the circulations are calculated by inserting the associated shear diameter, maximum velocity measured, and range values into the appropriate regression equations. Newman et al. (2013) found that the algorithm increased tornadic shear values appropriately and aided in the differentiation between tornadic and nontornadic scans.

2) Data thresholds and removal

Significant vertical shear near the surface can cause false high azimuthal shear values very close to the radar. To prevent these high values from corrupting the multiyear climatology, azimuthal shear values within a 5-km radius of each radar site are set to “missing.” While this will remove some “good” data, it will also remove a great deal of anomalously high azimuthal shear values that could corrupt the climatology. An example is shown in Fig. 5. This near-radar data removal is not applied to the rotation tracks in this paper, nor will it be used in the generation of rotation tracks in real time. It is used only in the climatology to avoid accumulations at areas where the azimuthal shear values are known to be poor.

Fig. 5.

(a) Low-level (0–3 km AGL) rotation track field from the 2 Mar 2012 outbreak from the Owensville, Indiana (KVWX) radar site before data removal near the radar site. Note the high values of azimuthal shear surrounding the radar site at the center of the image. (b) The same rotation track field after the removal of data within a 5-km radius of the radar site. This removal will only be performed for the climatology, not in real time.

Fig. 5.

(a) Low-level (0–3 km AGL) rotation track field from the 2 Mar 2012 outbreak from the Owensville, Indiana (KVWX) radar site before data removal near the radar site. Note the high values of azimuthal shear surrounding the radar site at the center of the image. (b) The same rotation track field after the removal of data within a 5-km radius of the radar site. This removal will only be performed for the climatology, not in real time.

When processing the 2D maximum azimuthal shear fields, all shear data not collocated with a given reflectivityQC value are removed so that only shear data associated with storm regions are retained. To retain meteorologically significant shear data in low-reflectivity hook-echo regions, a 5 × 5 dilation filter (Lakshmanan 2012) is applied to the reflectivityQC field. This operation assigns the maximum reflectivity value in a 5 × 5 moving window to each pixel, effectively expanding the areas of high-reflectivity values. The threshold operation is then performed on the dilated reflectivityQC field to help remove azimuthal shear associated with interference patterns, anomalous propagation, and other radar-related issues not successfully removed by the radar reflectivity quality control neural network. Setting this threshold at 40 dBZ retained the meteorological rotation signatures in the azimuthal shear fields while removing a great deal of shear collocated with nonmeteorological signatures (see Fig. 6).

Fig. 6.

(a) Reflectivity, (b) reflectivityQC, (c) radial velocity, and (d) azimuthal shear fields associated with a supercell over central AL at 2107 UTC 27 Apr 2011.

Fig. 6.

(a) Reflectivity, (b) reflectivityQC, (c) radial velocity, and (d) azimuthal shear fields associated with a supercell over central AL at 2107 UTC 27 Apr 2011.

d. Creation of rotation tracks

To better isolate the rotation tracks in the accumulated grid, new quality control strategies have been implemented on the input 2D maximum azimuthal shear fields for both low and midlevels of the atmosphere. Clusters of high azimuthal shear values at each time step are tested and removed if their sizes and/or data value distributions are not indicative of meteorological phenomena. If these remaining clusters in each time step are associated with lasting circulations, they make it into the final rotation track products.

1) Hysteresis segmentation

Before the circulation signatures in the 2D maximum azimuthal shear fields can be associated between time steps, they are isolated into clusters of high shear values using hysteresis segmentation. The term hysteresis (Jain 1989) refers to the lag observed between the application of an electromagnetic field and its subsequent effect on a substance. In image processing, the term refers to the lagging effect caused by using two thresholds: one to begin the thresholding process and another to end it. In this application, two data thresholds are maintained and a cluster is composed of contiguous pixels with values greater than the lower data threshold that contains at least one pixel with a data value greater than the higher threshold. The higher threshold identifies areas of high azimuthal shear associated with strong circulation and the lower hysteresis threshold grows the region around the high value to include all pixels associated with the circulation (see Fig. 7). Through experimentation on numerous tornadic and nontornadic case studies, it was determined that low and high data thresholds of 0.002 and 0.005 s−1, respectively, and a minimum size of 25 pixels performed well for isolating clusters of high azimuthal shear.

Fig. 7.

A schematic showing how hysteresis segmentation works. The first and last peaks of azimuthal shear values will be associated with clusters, while the middle peak will not be associated with a cluster because it contains no values above the higher data threshold.

Fig. 7.

A schematic showing how hysteresis segmentation works. The first and last peaks of azimuthal shear values will be associated with clusters, while the middle peak will not be associated with a cluster because it contains no values above the higher data threshold.

All pixels in the maximum 2D azimuthal shear layer fields not associated with identified clusters are then removed so that only the azimuthal shear clusters are accumulated over time to produce rotation tracks. This eliminates the low background azimuthal shear values not associated with circulations, making the tracks themselves more isolated, as seen in Fig. 8.

Fig. 8.

Low-level (0–3 km AGL) rotation tracks associated with the tornadic supercells that moved across MS and AL between 1600 UTC 27 Apr 2011 and 0000 UTC 28 Apr 2011. (a) Rotation tracks before hysteresis segmentation. (b) Rotation tracks after hysteresis segmentation is used to threshold the field.

Fig. 8.

Low-level (0–3 km AGL) rotation tracks associated with the tornadic supercells that moved across MS and AL between 1600 UTC 27 Apr 2011 and 0000 UTC 28 Apr 2011. (a) Rotation tracks before hysteresis segmentation. (b) Rotation tracks after hysteresis segmentation is used to threshold the field.

2) Multiple hypothesis tracking

Typically, clusters of high azimuthal shear values associated with mesocyclones persist through many time steps, whereas clusters associated with the remaining nonmeteorological signatures typically only appear sporadically. To isolate the meteorological rotation tracks, these nonmeteorological shear clusters need to be removed from the accumulated fields. MHT techniques (Cox and Hingorani 1996) are used to isolate the continuous tracks of azimuthal shear clusters associated with storms and remove the short-lived nonmeteorological azimuthal shear clusters.

MHT techniques have been used in the fields of video processing and military target tracking for years and recently have been adopted by the meteorology community (e.g., Root et al. 2011). The technique attempts to associate objects, in this case azimuthal shear clusters, throughout time. This approach is innovative as it considers time associations globally and makes association decisions that can be deferred until additional information is available. If the algorithm is not certain whether an existing track should be associated with cluster A or cluster B in the current time step, for example, it can create two hypotheses (see Fig. 9). Both possibilities are then propagated forward in time to when enough information should be available to determine which hypothesis is most likely. The clusters that do not meet the minimum longevity threshold (two time steps or roughly 10 min) are then retroactively pruned so that they are not admitted into the rotation track fields.

Fig. 9.

A schematic showing how the MHT algorithm tracks clusters and generates hypotheses. Solid shapes represent the position of azimuthal shear clusters at the given time step. Solid arrows show the movement of clusters between time steps. Dotted shapes show the projected locations of the azimuthal shear clusters in the next time steps. Dashed arrows show the projected movement of the clusters between time steps. Shape A shows the projected location of the original cluster at time t. Shape B represents the actual location of the original cluster at time t. The algorithm generates two hypotheses for time t + 1 (C and D). If the location of B is an error, then at time t + 1 the cluster should move to position D. If the location of B is not an error but a change in motion, then B should move to location C at time t + 1. At time t + 1, the discovery of the target at either C or D will confirm one hypothesis and disprove the other. The disproved cluster is deleted and the confirmed one continues on to the next time step. In this case, given the strange shape and size of B, it is likely that this cluster is associated with a nonmeteorological shear signature and will not persist in time t + 1.

Fig. 9.

A schematic showing how the MHT algorithm tracks clusters and generates hypotheses. Solid shapes represent the position of azimuthal shear clusters at the given time step. Solid arrows show the movement of clusters between time steps. Dotted shapes show the projected locations of the azimuthal shear clusters in the next time steps. Dashed arrows show the projected movement of the clusters between time steps. Shape A shows the projected location of the original cluster at time t. Shape B represents the actual location of the original cluster at time t. The algorithm generates two hypotheses for time t + 1 (C and D). If the location of B is an error, then at time t + 1 the cluster should move to position D. If the location of B is not an error but a change in motion, then B should move to location C at time t + 1. At time t + 1, the discovery of the target at either C or D will confirm one hypothesis and disprove the other. The disproved cluster is deleted and the confirmed one continues on to the next time step. In this case, given the strange shape and size of B, it is likely that this cluster is associated with a nonmeteorological shear signature and will not persist in time t + 1.

An association cost matrix is constructed so an entry indicates the cost of matching cluster at one time step with cluster at the next time step. Each association has a computed cost based on cluster sizes, ages, proximities to clusters from previous time steps, and other characteristics. The associations with the lowest cost are made. Enumeration of all the hypothesis matrices to find the lowest costs can increase exponentially with each time step, so a technique based on Murty (1968) is used to prune the set to retain only the k best associations at each time step. The algorithm is illustrated in Fig. 10. For more details and to see quantitative improvements in simulated fields through the use of MHT, the reader is referred to Lakshmanan et al. (2013).

Fig. 10.

MHT flowchart adapted from Root et al. (2011).

Fig. 10.

MHT flowchart adapted from Root et al. (2011).

In an effort to remove any lingering nonmeteorological clusters, a data value distribution threshold was set. It was observed that the majority of clusters associated with meteorological clusters exhibited unimodal distributions of azimuthal shear data values with central tendencies while the nonmeteorological clusters typically exhibited uniform distributions of very high azimuthal shear values (see examples in Fig. 11). To account for this, clusters are pruned if more than 80% of their pixel values are greater than or equal to 0.02 s−1. This has only an incremental impact on the field itself since most nonmeteorological clusters are removed before this point.

Fig. 11.

(a) An example of a cluster of high azimuthal shear values associated with a mesocyclone circulation and (b) the histogram of its data values after hysteresis segmentation. (c) An example of a cluster of high azimuthal shear values associated with a nonmeteorological artifact and (d) the histogram of its data values after hysteresis segmentation. Note the nearly uniform distribution of very high values.

Fig. 11.

(a) An example of a cluster of high azimuthal shear values associated with a mesocyclone circulation and (b) the histogram of its data values after hysteresis segmentation. (c) An example of a cluster of high azimuthal shear values associated with a nonmeteorological artifact and (d) the histogram of its data values after hysteresis segmentation. Note the nearly uniform distribution of very high values.

Bands of high azimuthal shear values associated with linear meteorological phenomena like outflow boundaries and bow echoes also appear in the rotation track fields. In an effort to isolate the mesocyclone signatures from these bands of shear, all clusters were fit to ellipses and their aspect ratios were calculated. After testing many different thresholds, it was determined that size, data value distribution, and aspect ratio information could not be successfully used to discriminate between mesocyclone clusters and shear band clusters. Because of this, the band signatures remain in the rotation track fields for now.

3. Results

The quality control techniques discussed in the methods section were developed and tuned through testing on a variety of tornadic and nontornadic cases. The specific impacts of each technique will now be discussed and demonstrated in this section using cases that were not part of this training dataset.

a. New velocity dealiasing techniques

Prior to this study, the LED dealiasing technique (Eilts and Smith 1990), the default method used for real-time processing of WSR-88D data, was used in the creation of rotation tracks. Recently, it was found that using the vertical profile of horizontal wind from the 20-km RUC point sounding at radar sites as estimates of the environmental wind in the algorithm helped to alleviate some of the dealiasing issues, though many still persisted.

The 2D dealiasing technique described by Jing and Wiener (1993) was tested and appears to perform much better at properly dealiasing the velocity fields because of the large reduction in radial spikes and nonmeteorological velocity signatures. Using the RUC wind profiles as environmental estimates made some additional improvements as well. As seen by the representative example in Fig. 12, the 2D technique dealiases correctly many of the areas that the LED technique did not. Almost all of the noisy, high azimuthal shear values associated with dealiasing issues are removed by using the Jing and Wiener technique, making the rotation tracks much easier to interpret. A quantitative study to determine the best velocity dealiasing techniques is on going.

Fig. 12.

Low-level (0–3 km AGL) rotation tracks associated with the 27 Apr 2011 tornado outbreak across MS and AL produced using different velocity dealiasing techniques, no thresholds, and no MHT. (a) Tracks made using velocity dealiased with the LED algorithm. (b) Tracks made using velocity dealiased with the LED algorithm with 20-km RUC input sounding. (c) Tracks made using velocity dealiased with the Jing and Wiener (1993) technique. (d) Tracks made using velocity dealiased with the Jing and Wiener (1993) technique with 20-km RUC input sounding.

Fig. 12.

Low-level (0–3 km AGL) rotation tracks associated with the 27 Apr 2011 tornado outbreak across MS and AL produced using different velocity dealiasing techniques, no thresholds, and no MHT. (a) Tracks made using velocity dealiased with the LED algorithm. (b) Tracks made using velocity dealiased with the LED algorithm with 20-km RUC input sounding. (c) Tracks made using velocity dealiased with the Jing and Wiener (1993) technique. (d) Tracks made using velocity dealiased with the Jing and Wiener (1993) technique with 20-km RUC input sounding.

b. Range correction and reflectivity thresholds

Whereas the azimuthal shear range correction does not make many visible improvements to the rotation track fields, the azimuthal shear values in storm circulations are more accurate estimates of the actual shear values. The reflectivityQC threshold below which all collocated shear data are removed is increased from 20 to 40 dBZ. This “stamping out” of azimuthal shear by the dilated reflectivityQC field helps to further isolate the azimuthal shear signatures associated with storms, as seen in Fig. 13.

Fig. 13.

Low-level (0–3 km AGL) rotation tracks from the 27 Apr 2011 tornado outbreak. (a) Tracks created using the Jing and Wiener (1993) velocity dealiasing method with the input 20-km RUC sounding, without the azimuthal shear range correction, and without the increased reflectivityQC threshold. (b) As in (a), but with the azimuthal shear range correction and the increased reflectivityQC threshold. Note that the circled area of high azimuthal shear values associated with dealiasing errors over central TN is much less prominent in (b).

Fig. 13.

Low-level (0–3 km AGL) rotation tracks from the 27 Apr 2011 tornado outbreak. (a) Tracks created using the Jing and Wiener (1993) velocity dealiasing method with the input 20-km RUC sounding, without the azimuthal shear range correction, and without the increased reflectivityQC threshold. (b) As in (a), but with the azimuthal shear range correction and the increased reflectivityQC threshold. Note that the circled area of high azimuthal shear values associated with dealiasing errors over central TN is much less prominent in (b).

c. Hysteresis segmentation

The rotation tracks are even further isolated from any background azimuthal shear values through the use of hysteresis segmentation. The tracks are more isolated and easier to interpret (Fig. 14) since only the azimuthal shear clusters are accumulated over time rather than the entire field, including the low background values.

Fig. 14.

Low-level (0–3 km AGL) rotation tracks from the 27 Apr 2011 tornado outbreak. (a) Tracks created using the Jing and Wiener (1993) velocity dealiasing method with the input 20-km RUC sounding, with range correction, increased reflectivityQC thresholds, and without hysteresis segmentation. (b) As in (a), but with hysteresis segmentation implemented.

Fig. 14.

Low-level (0–3 km AGL) rotation tracks from the 27 Apr 2011 tornado outbreak. (a) Tracks created using the Jing and Wiener (1993) velocity dealiasing method with the input 20-km RUC sounding, with range correction, increased reflectivityQC thresholds, and without hysteresis segmentation. (b) As in (a), but with hysteresis segmentation implemented.

d. Multiple hypothesis tracking

After using hysteresis segmentation to form azimuthal shear clusters, the MHT algorithm is used to isolate the persistent clusters associated with storm-scale circulations from the nonmeteorological clusters associated with any remaining poor velocity dealiasing signatures. Figure 15 illustrates the removal of some small, leftover circulations from the 27 April 2011 case by the MHT algorithm.

Fig. 15.

Low-level (0–3 km AGL) rotation tracks across MS and AL from the 27 Apr 2011 tornado outbreak. (a) Tracks created using the Jing and Wiener (1993) velocity dealiasing method with the input 20-km RUC sounding, with range correction, increased reflectivityQC thresholds, hysteresis segmentation, and without MHT. (b) As in (a), but with MHT. Note the several small clusters in west-central MS are removed.

Fig. 15.

Low-level (0–3 km AGL) rotation tracks across MS and AL from the 27 Apr 2011 tornado outbreak. (a) Tracks created using the Jing and Wiener (1993) velocity dealiasing method with the input 20-km RUC sounding, with range correction, increased reflectivityQC thresholds, hysteresis segmentation, and without MHT. (b) As in (a), but with MHT. Note the several small clusters in west-central MS are removed.

Figure 16 illustrates the differences between the original rotation track fields and the rotation track fields after the quality control efforts were implemented on four recent tornadic cases. Because of space constraints, only the overall improvements are shown. Radial spikes, which were especially problematic in the 16 April 2011 case over Virginia and North Carolina, were successfully removed. The low background azimuthal shear values and nearly all azimuthal shear values not associated with storms are removed in all four cases. Broad areas of nonmesocyclone shear still exist near some radar sites, but overall a significant improvement is seen in the quality of the data and the ease of interpretation.

Fig. 16.

The impacts (left) before and (right) after quality control efforts on low-level (0–3 km AGL) rotation track products associated with four recent tornado events: (a),(b) the 27 Apr 2011 event, (c),(d) the 16 Apr 2011 event, (e),(f) the 24 May 2011 event, and (g),(h) the 2 Mar 2012 event.

Fig. 16.

The impacts (left) before and (right) after quality control efforts on low-level (0–3 km AGL) rotation track products associated with four recent tornado events: (a),(b) the 27 Apr 2011 event, (c),(d) the 16 Apr 2011 event, (e),(f) the 24 May 2011 event, and (g),(h) the 2 Mar 2012 event.

To get a more quantitative idea of how MHT impacts rotation track fields, a cluster-tracking algorithm described in Lakshmanan et al. (2003) was used to identify and track the number of clusters in the 24 May 2011 event before and after implementing MHT. Before MHT was implemented, 62 different clusters were identified and tracked, whereas after MHT, only 41 clusters were identified and tracked in the 8-h case. The low-level rotation tracks associated with these post-MHT clusters compare well to reported tornado tracks in the several cases visually examined. Figure 17 shows how the low-level rotation track products associated with tornado damage paths from the 24 May 2011 event across central Oklahoma compare to the enhanced Fujita scale (EF) ratings.

Fig. 17.

(a) Low-level (0–3 km AGL) rotation tracks over a 145-min period associated with the 24 May 2011 tornado outbreak across central OK. (b) Plotted tornado tracks and EF-scale intensity ratings associated with tornadoes that occurred during that same time period as in (a). (Survey data provided courtesy of K. Ortega, B. Smith, and G. Garfield.)

Fig. 17.

(a) Low-level (0–3 km AGL) rotation tracks over a 145-min period associated with the 24 May 2011 tornado outbreak across central OK. (b) Plotted tornado tracks and EF-scale intensity ratings associated with tornadoes that occurred during that same time period as in (a). (Survey data provided courtesy of K. Ortega, B. Smith, and G. Garfield.)

4. Conclusions

NSSL rotation track products are valuable tools in disaster response situations. They allow users to quickly assess both the spatial extent of mesocyclone circulations as seen by radar over time and assess their relative intensities. While these products were useful, a great deal of contamination initially was present because of poor velocity dealiasing, ground clutter, radar test patterns, and spurious shear values. These nonmeteorological signatures made the tracks nearly impossible to see in some extreme cases.

To mitigate these problems for both real-time use and for a multiyear rotation track climatology as part of the MYRORSS project, quality control strategies were developed and implemented. A 2D velocity dealiasing technique using 20-km RUC wind data as input made large visual improvements in the quality of the initial radial velocity field. An azimuthal shear range correction algorithm and some simple data thresholds were added to the LLSD shear algorithm. Hysteresis segmentation was used to isolate clusters of high azimuthal shear associated with mesocyclone circulations in each time step of the 2D maximum azimuthal shear fields and MHT techniques were used to associate them throughout time. Any clusters that did not persist for at least 10 min (2 time steps) or were composed of an unrealistic distribution of high azimuthal shear values were pruned. The remaining clusters were kept and used to create the rotation track fields.

While a few issues like the broad shear signatures around some radar sites remain, overall the rotation track fields show a great deal of qualitative improvement after the implementation of the quality control efforts. Whereas the tracks associated with the mesocyclones initially were diluted by background noise and nonmeteorological signatures, they are now isolated and easier to interpret. Incorporating these improvements, processing of the MYRORSS rotation track climatology should begin in the near future.

Acknowledgments

Funding for the authors was provided under NOAA–OU Cooperative Agreement NA17RJ1227. The authors thank Kiel Ortega, Kevin Manross, John Cintineo, and Jennifer Newman for all their help and advice. The authors also thank Gabe Garfield, Kiel Ortega, and Brandon Smith for allowing us to use their damage survey data for the 24 May 2011 tornadoes.

REFERENCES

REFERENCES
Bargen
,
D. W.
, and
R. C.
Brown
,
1980
: Interactive radar velocity unfolding. Preprints, 19th Conf. on Radar Meteorology, Miami, FL, Amer. Meteor. Soc., 278–283.
Benjamin
,
S. G.
, and
Coauthors
,
2004
:
An hourly assimilation–forecast cycle: The RUC
.
Mon. Wea. Rev.
,
132
,
495
518
.
Bergen
,
W. R.
, and
S. C.
Albers
,
1988
:
Two- and three-dimensional dealiasing of Doppler radar velocities
.
J. Atmos. Oceanic Technol.
,
5
,
305
319
.
Boren
,
T. A.
,
J. R.
Cruz
, and
D. S.
Zrnić
,
1986
: An artificial intelligence approach to Doppler weather radar velocity dealiasing. Proc. 23rd Conf. on Radar Meteorology, Snowmass, CO, Amer. Meteor. Soc., 107–110.
Brown
,
R. A.
,
L. R.
Lemon
, and
D. W.
Burgess
,
1978
:
Tornado detection by pulsed Doppler radar
.
Mon. Wea. Rev.
,
106
,
29
38
.
Brown
,
R. A.
,
V. T.
Wood
, and
D.
Sirmans
,
2002
:
Improved tornado detection using simulated and actual WSR-88D data with enhanced resolution
.
J. Atmos. Oceanic Technol.
,
19
,
1759
1771
.
Cintineo
,
J.
,
T.
Smith
,
V.
Lakshmanan
, and
S.
Ansari
,
2011
: An automated system for processing the Multi-Year Reanalysis of Remotely Sensed Storms (MYRORSS). Preprints, 27th Conf. on Interactive Information Processing Systems (IIPS), Seattle, WA, Amer. Meteor. Soc., J9.3. [Available online at https://ams.confex.com/ams/91Annual/webprogram/Paper182332.html.]
Cintineo
,
J.
,
T.
Smith
,
V.
Lakshmanan
,
H. E.
Brooks
, and
K. L.
Ortega
,
2012
:
An objective high-resolution hail climatology of the contiguous United States
.
Wea. Forecasting
,
27
,
1235
1248
.
Clark
,
A. J.
,
J. S.
Kain
,
P. T.
Marsh
,
J.
Correia
Jr.
,
M.
Xue
, and
F.
Kong
,
2012
:
Forecasting tornado pathlengths using a three-dimensional object algorithm applied to convection-allowing forecasts
.
Wea. Forecasting
,
27
,
1090
1113
.
Clark
,
A. J.
,
J.
Gao
,
P.
Marsh
,
T.
Smith
,
J.
Kain
,
J.
Correia
,
M.
Xue
, and
F.
Kong
,
2013
:
Tornado pathlength forecasts from 2010 to 2011 using ensemble updraft helicity
.
Wea. Forecasting
,
28
,
387
407
.
Cox
,
I.
, and
S. L.
Hingorani
,
1996
:
An efficient implementation of Reid's multiple hypothesis tracking algorithm and its evaluation for the purpose of visual tracking
.
IEEE Trans. Pattern Anal. Mach. Intell.
,
18
,
138
150
.
Crum
,
T. D.
, and
R. K.
Alberty
,
1993
:
The WSR-88D and the WSR-88D Operational Support Facility
.
Bull. Amer. Meteor. Soc.
,
74
,
1669
1687
.
Eilts
,
M. D.
, and
S. D.
Smith
,
1990
:
Efficient dealiasing of Doppler velocities using local environment constraints
.
J. Atmos. Oceanic Technol.
,
7
,
118
128
.
Elmore
,
K. M.
,
E. D.
Albo
,
R. K.
Goodrich
, and
D. J.
Peters
,
1994
: NASA/NCAR airborne and ground-based wind shear studies. Final Rep., Contract NCC1-155, 343 pp.
Forbes
,
G. S.
,
1981
:
On the reliability of hook echoes as tornado indicators
.
Mon. Wea. Rev.
,
109
,
1457
1466
.
Fulton
,
R.
,
D.
Breidenback
,
D.
Miller
, and
T.
O'Bannon
,
1998
:
The WSR-88D rainfall algorithm
.
Wea. Forecasting
,
13
,
377
395
.
Jain
,
A.
,
1989
: Fundamentals of Digital Image Processing. Prentice Hall, 569 pp.
Jing
,
Z.
, and
G.
Wiener
,
1993
:
Two-dimensional dealiasing of Doppler velocities
.
J. Atmos. Oceanic Technol.
,
10
,
798
808
.
Kain
,
J. S.
,
S. R.
Dembek
,
S. J.
Weiss
,
J. L.
Case
,
J. J.
Levit
, and
R. A.
Sobash
,
2010
:
Extracting unique information from high-resolution forecast models: Monitoring selected fields and phenomena every time step
.
Wea. Forecasting
,
25
,
1536
1542
.
Kessinger
,
C.
,
S.
Ellis
, and
J.
Van Andel
,
2003
: The radar echo classifier: A fuzzy logic algorithm for the WSR-88D. Preprints, Third Conf. on Artificial Applications to the Environmental Sciences, Long Beach, CA, Amer. Meteor. Soc., P1.6. [Available online at https://ams.confex.com/ams/pdfpapers/54946.pdf.]
Lakshmanan
,
V.
,
2012
: Neighborhood and window operations. Automating the Analysis of Spatial Grids: A Practical Guide to Data Mining Geospatial Images for Human and Environmental Applications, V. Lakshmanan, Ed., Springer, 129–171.
Lakshmanan
,
V.
,
R.
Rabin
, and
V.
DeBrunner
,
2003
:
Multiscale storm identification and forecast
.
Atmos. Res.
,
67–68
,
367
380
,
doi:10.1016/S0169-8095(03)00068-1
.
Lakshmanan
,
V.
,
T.
Smith
,
K.
Hondl
,
G. J.
Stumpf
, and
A.
Witt
,
2006
:
A real-time, three-dimensional, rapidly updating, heterogeneous radar merger technique for reflectivity, velocity, and derived products
.
Wea. Forecasting
,
21
,
802
823
.
Lakshmanan
,
V.
,
A.
Fritz
,
T.
Smith
,
K.
Hondl
, and
G. J.
Stumpf
,
2007a
:
An automated technique to quality control radar reflectivity data
.
J. Appl. Meteor. Climatol.
,
46
,
288
305
.
Lakshmanan
,
V.
,
T.
Smith
,
G. J.
Stumpf
, and
K.
Hondl
,
2007b
:
The Warning Decision Support System—Integrated Information
.
Wea. Forecasting
,
22
,
596
612
.
Lakshmanan
,
V.
,
J.
Zhang
, and
K.
Howard
,
2010
:
A technique to censor biological echoes in radar reflectivity data
.
J. Appl. Meteor. Climatol.
,
49
,
453
462
.
Lakshmanan
,
V.
,
M.
Miller
, and
T.
Smith
,
2013
:
Quality control of accumulated fields by applying spatial and temporal constraints
.
J. Atmos. Oceanic Technol.
,
30
,
745
758
.
Liu
,
S.
,
M.
Xue
, and
Q.
Xu
,
2007
:
Using wavelet analysis to detect tornadoes from Doppler radar radial-velocity observations
.
J. Atmos. Oceanic Technol.
,
24
,
344
359
.
Merritt
,
M. W.
,
1984
: Automatic velocity dealiasing for real-time applications. Proc. 22nd Conf. on Radar Meteorology, Zurich, Switzerland, Amer. Meteor. Soc., 528–533.
Mitchell
,
E. D.
, and
K. L.
Elmore
,
1998
: A technique for identifying regions of high shear associated with mesocyclones and tornadic vortex signatures. Preprints, 14th Int. Conf. on Interactive Information and Processing Systems (IIPS) for Meteorology, Oceanography, and Hydrology, Phoenix, AZ, Amer. Meteor. Soc., 312–315.
Mitchell
,
E. D.
,
S. V.
Vasiloff
,
G. J.
Stumpf
,
A.
Witt
,
M. D.
Eilts
,
J. T.
Johnson
, and
K. W.
Thomas
,
1998
:
The National Severe Storms Laboratory tornado detection algorithm
.
Wea. Forecasting
,
13
,
352
366
.
Murty
,
K. G.
,
1968
:
An algorithm for ranking all the assignments in order of increasing cost
.
Oper. Res.
,
16
,
682
687
.
Newman
,
J. F.
,
V.
Lakshmanan
,
P. L.
Heinselman
,
M. B.
Richman
, and
T. M.
Smith
,
2013
:
Range-correcting azimuthal shear in Doppler radar data
.
Wea. Forecasting
,
28
,
194
211
.
NOAA/NSSL
, cited
2011
: NOAA technology helps American Red Cross respond faster. [Available online at https://secure.nssl.noaa.gov/briefings/2011/06/noaa-technology-helps-american-red-cross-respond-faster/.]
Ray
,
P. S.
, and
C.
Ziegler
,
1977
:
De-aliasing first-moment Doppler estimates
.
J. Appl. Meteor.
,
16
,
563
564
.
Root
,
B.
,
M.
Yeary
, and
T.-Y.
Yu
,
2011
: Novel storm cell tracking with multiple hypothesis tracking. Preprints, 27th Conf. on Interactive Information Processing Systems (IIPS), Seattle, WA, Amer. Meteor. Soc., 8B.3. [Available online at https://ams.confex.com/ams/91Annual/webprogram/Paper184250.html.]
Smith
,
T. M.
, and
K. L.
Elmore
,
2004
: The use of radial velocity derivatives to diagnose rotation and divergence. Preprints, 11th Conf. on Aviation, Range, and Aerospace, Hyannis, MA, Amer. Meteor. Soc., P5.6. [Available online at https://ams.confex.com/ams/pdfpapers/81827.pdf.]
Steiner
,
M.
, and
J. A.
Smith
,
2002
:
Use of three-dimensional reflectivity structure for automated detection and removal of nonprecipitating echoes in radar data
.
J. Atmos. Oceanic Technol.
,
19
,
673
686
.
Stumpf
,
G. J.
,
A.
Witt
,
E. D.
Mitchell
,
P. L.
Spencer
,
J. T.
Johnson
,
M. D.
Eilts
,
K. W.
Thomas
, and
D. W.
Burgess
,
1998
:
The National Severe Storms Laboratory mesocyclone detection algorithm for the WSR-88D
.
Wea. Forecasting
,
13
,
304
326
.
Wood
,
V. T.
, and
R. A.
Brown
,
1997
:
Effects of radar sampling on single-Doppler velocity signatures of mesocyclones and tornadoes
.
Wea. Forecasting
,
12
,
928
938
.