An algorithm for automatic storm identification, tracking, and nowcasting using tree structure representation of radar reflectivity images is proposed. The algorithm aims to track and nowcast different kinds of storm objects (stratiform regions, convective storms, and storm cells) simultaneously and to preserve their spatial relationships in the tracking and nowcasting processes. The algorithm applies a region tree structure to represent intensity regions and their spatial relationships in radar reflectivity images. Storm objects are identified by clustering regions within the region tree structure. Storm tracking is accomplished using an iterative region tree matching algorithm. Storm nowcasting applies the tree structure to the nowcasting of the internal structures of storm objects. Using eight cases with different storm types, a comparative evaluation with the enhanced Thunderstorm Identification, Tracking, Analysis, and Nowcasting (ETITAN) method and the Storm Cell Identification and Tracking (SCIT) method has shown that the proposed tree-based storm-tracking algorithm achieves better performance in storm tracking and nowcasting. The critical success index (CSI) value of storm association is 78.16% for the tree-based method, as compared with 74.88% for SCIT and 74.71% for ETITAN. The CSI value of an 18-min nowcast is 29.02% for the tree-based method, as compared with 24.98% for SCIT and 24.44% for ETITAN. The evaluation also shows that the tree-based method is able to nowcast the internal structure of storms and therefore produces small mean absolute errors (MAE).
Storm tracking, using weather radar, is an important component of weather nowcasting systems. A storm-tracking algorithm is usually used to extract storm tracks and other time series of storm attributes, such as volume, height, and vertically integrated liquid (VIL) water. The output of a storm-tracking algorithm can be used to predict the position and starting time of severe weather, as well as study the mechanisms of storm evolution.
Radar-based nowcasting systems also include area-based, statistical, and probabilistic nowcasting approaches (Ruzanski et al. 2011). Area-based approaches (e.g., Rinehart and Garvey 1978; Rinehart 1981; Rood 1987; Tuttle and Foote 1990; Li et al. 1995; Wolfson et al. 1999; Lai 1999; Germann and Zawadzki 2002; Turner et al. 2004) estimate motion vectors within radar reflectivity echoes. For example, Rinehart and Garvey (1978) estimated motion vectors by determining the locations of maximum cross-correlation coefficients between subgrids of successive radar fields. The area-based approaches extrapolate the radar image and produce a forecast using the motion vectors. The advantage of area-based nowcasting approaches is that they can provide accurate motion vector within each grid. The statistical and probability approaches aim at providing a measure of uncertainty in nowcasting. The statistical approaches (e.g., Xu et al. 2005; Seed 2003; Fox and Wikle 2005) represent atmospheric evolution using statistical models at various scales. The probabilistic approaches (e.g., Germann and Zawadzki 2004; Megenhardt et al. 2004; Schmeits et al. 2008; Dance et al. 2010) predict the probability of encountering a precipitation rate above a certain threshold for each grid. When compared with the above-mentioned approaches, object-based storm-tracking methods are often used to forecast and analyze convective storm systems.
The object-based methods, summarized in Fig. 1a, are usually composed of three parts: identifying, tracking, and forecasting of storm objects. The term storm objects often refers to convective storms or individual convective cells in previous publications. A variety of storm object identification methods exist. For example, the Thunderstorm Identification, Tracking, Analysis, and Nowcasting (TITAN) method (Dixon and Wiener 1993) identified three-dimensional (3D) regions with intensity higher than a threshold as convective storms. Later, Han et al. (2009) adopted multiple thresholds in their enhanced TITAN (ETITAN) algorithm to identify convective storms and adopted dilation and erosion operations in mathematical morphology to mitigate the false merger problem during storm identification. In addition, Lakshmanan et al. (2003) introduced a k-means cluster method to identify storms at different scales. Rossi et al. (2015) identified individual storm systems by clustering convective cells. The Storm Cell Identification and Tracking (SCIT) algorithm (Johnson et al. 1998) used multiple thresholds to identify three-dimensional storm cells. Handwerker (2002) initially applied a low threshold to detect regions of intense precipitation (ROIP), then determined the maximal reflectivity within each ROIP, and finally defined a storm cell as a contiguous region with reflectivity not lower than the maximal reflectivity minus a threshold (e.g., 10 dBZ) within the ROIP.
After identification, the storm objects between volume scans are associated with each other and grouped to form storm tracks. For example, TITAN matched storm sets using a combinatorial optimization method that joined storms of similar characteristic and produced short tracks. SCIT associated storm cells using a centroid projection method that first predicted the locations of previous storm cells based on their motion vectors, and then assigned a previous storm cell to the current cell that was closest to its predicted location. Rosenfeld (1987) and Kyznarová and Novák (2009) matched storm cells using the area-overlapping method. In this method, a storm cell at the current time matched with a storm cell at a previous time with which it had maximum overlap. ETITAN first matched storms using an area-overlapping technology, and then used a combinatorial optimization method to match the remaining storms. Lakshmanan and Smith (2010) proposed a hybrid tracking method that combined the centroid projection method and the combinatorial optimization method. Moreover, the method restricted the matching process within an area-based search radius and considered long-lived cells first in the matching process. Rossi et al. (2015) clustered convective cells that were identified in the current and previous scans to track convective storms.
The object-based storm-tracking approaches can nowcast storm objects based on motion vectors derived from storm tracks or some other methods. For example, the SCIT and TITAN methods estimated storm velocities by linear fitting of storm centroid tracks; the ETITAN method first estimated grid motion vectors within the entire radar field using the Tracking Radar Echoes by Correlation (TREC) algorithm (Rinehart and Garvey 1978), and then calculated the velocity of a storm as the average grid motion vector within it.
The object-based storm-tracking approaches have two limitations. One limitation is that object-based methods often rely on the outputs of their identification methods; that is, a tracking method is suitable only for storm objects (convective storms or convective cells) that are identified by the method itself. It would be useful to develop a method that could track different kinds of storm objects at different scales simultaneously. Another limitation with most object-based storm nowcasting methods is that they regard a storm object as a rigid body. In other words, the distance between any two given points of a storm object remains constant when the storm object is nowcasted. However, storms always exhibit internal relative motion and changing shapes. For example, the motion vector of a convective storm could be different from that of its internal storm cell.
Considering the above-mentioned problems, we propose a tree-based storm-tracking method in this paper as diagramed in Fig. 1b. The tree-based storm-tracking method is built on the foundation that a radar reflectivity image is represented using a region tree structure. Each node in the tree corresponds to a reflectivity region, and each link denotes the overlapping relationship between various reflectivity regions. Identification of storm objects is based on combination of regions within the tree; thus, different kinds of storm objects could be identified by applying different combinatorial methods simultaneously. This storm-tracking approach does not track storm objects directly. Instead, it first associates the regions in the tree, and then deduces the matching results of storm objects. In this way, the storm-tracking method does not depend on the storm identification results and could track different kinds of storm objects simultaneously. In the nowcasting stage, regions within the tree are nowcasted individually and then combined into the nowcasted storm objects.
When compared with the traditional object-based approaches, the goals of tree-based storm nowcasting approach are 1) different kinds of storm objects, which have various scales, are identified, tracked, and nowcasted simultaneously; 2) the spatial relationships among storm objects are also identified; 3) the tracking and nowcasting methods are independent from the storm identification results; and 4) the internal relative motion and change within a storm object are considered in the storm object nowcasting method.
The tree-based storm-tracking approach not only nowcasts storm objects but also nowcasts the entire radar field. In this way, the tree-based nowcasting approach is similar with that of Wolfson et al. (1999). They nowcasted storm envelopes at two scales using an area-based approach. As compared with their approach, the tree-based approach in the present study could nowcast storm envelopes at different scales.
The region tree structure is a common tool that is employed during image analysis (e.g., Peura 2001). Peura (1999) proposed a property tree–based cell tracking method. He also represented a radar reflectivity image using a tree structure, but he did not apply the tree structure during storm identification and nowcasting procedures. Moreover, his tree matching method cannot handle storm splits or mergers. The tree-based storm-tracking approach in the present study combines a region tree structure representation method and an object-based storm-tracking method, and could be regarded as an extension of Peura’s work.
The remainder of the present paper is organized as follows: Section 2 illustrates the radar data and data preprocessing. Section 3 first introduces the region tree structure to represent radar reflectivity data. Then, the methods of identifying, tracking, and nowcasting of storm objects based on region tree structure representation are described. Section 4 describes experimental work, and finally section 5 provides a summary.
2. Data and preprocess
The archived Weather Surveillance Radar-1988 Doppler (WSR-88D) level II dataset used in this study was obtained from the National Climatic Data Center (NCDC). Eight cases, which included different storm types [e.g., isolated nonsevere, isolated severe, minisupercell, and mesoscale convective system (MCS)] and came from five WSR-88D radar sites, were selected to form a test dataset. Each case started at 1800 UTC and ended at 2100 UTC (Table 1). These cases were chosen because they had been employed by Johnson et al. (1998) and Lakshmanan and Smith (2010) in their storm-tracking test. Storm numbers in these cases were verified by human observers.
The radar data in the dataset have a range resolution of 1 km and an azimuth increment of 1°; the data cover a 230 km × 230 km domain centered on the radar site. The intensity of radar reflectivity ranges from −25 to 75 dBZ and is discretized with a precision of 1 dBZ. The temporal resolutions of radar data depend on the volume coverage pattern (VCP) mode (see Table 1): 6 min for VCP 21 and 5 min for VCP 11.
The radar reflectivity data from the lowest elevation angle (0.5°) are used in this work and are shown in a plan position indicator (PPI) display. The radar reflectivity data are transformed into Cartesian coordinates before the algorithm begins. The spatial resolution of the Cartesian grid is about 1 km × 1 km. Median filtering, which replaces each grid with the median of the neighboring grid, is performed in 3 pixel × 3 pixel neighborhoods to remove noise in the images. In addition to a median filter, radar reflectivity measurement below 20 dBZ is set to 0 dBZ because the radar reflectivity of each convective storm, storm cell, or stratiform region is larger than 20 dBZ.
a. Region tree structure representation of radar image
The first step of this work is to represent the radar image using a tree structure. A tree structure can be defined as a collection of nodes together with links between the nodes: , where is the set of nodes and is the set of edges linking . The terms listed in Table 2 are often used when referring to a tree structure.
In the image processing community, tree structure is often used to represent the hierarchical structure of regions at different threshold levels in an image. Each node in the tree corresponds to a region and each link represents the overlapping relationships between two regions with different intensities (Houlahan and Scalo 1992). Steps involved in generating a tree structure from an image P are given as follows.
Select a set of thresholds, , where . Let , , and .
Threshold the image P using ; that is, a pixel in the P is set to 0 if its intensity is less than . The output image is denoted by .
Identify all regions of connected nonzero pixels in using a region growth method (Robert and Shapiro 1992, 28–48) and compute the properties of these regions, such as area, center coordinate, and average intensity.
Create a node for each region r to record the properties of the region r and coordinates of the pixels within the region r. Add to the set .
If , then find all possible region pairs , where and are regions in and , respectively, and contains . Add node pair to the set .
Let . If , then the process is terminated. Otherwise, return to step 2.
Figure 2 illustrates how to represent an image using a region tree structure. The P is segmented using three thresholds (30, 35, and 40 dBZ), resulting in , , and , respectively. When the regions in , , and are identified and the overlapping relationships between these regions are determined, then the region tree structure is constructed.
The thresholds used to segment radar reflectivity images in the method are 0, 20, 25, 30, 35, 40, 45, 50, 55, and 60 dBZ. The threshold 0 dBZ is applied because we want to represent the entire image using one tree structure to simplify the later tree matching process. When more thresholds are used, the constructed tree structure will be more integrated, but the trade-off is the cost of calculation. In practice, the thresholds can be selected as needed. For example, if one is more interested in the individual storm cells, then the thresholds can be selected as 0, 40, 42, 44, …, 58, 60 dBZ. In this way, the intervals between the thresholds are smaller, and thus the constructed region tree structure could describe storm cells more precisely.
b. Identification of storm objects in rainfall systems
In the present study, a storm object is defined as a region of high-intensity precipitation in a radar reflectivity image, where each pixel is larger than a threshold . The term storm object refers to any of the following three concepts: a convective storm, a convective cell, or a stratiform region. Here, a convective cell is defined as an individual updraft core inside a storm, and a convective storm is a larger convective system that may encompass several individual convective cells. The goal of storm object identification is to identify different kinds of storm objects simultaneously and to recognize the spatial relationships among them. Stratiform regions are also identified because they are useful during storm classification (Parker and Johnson 2000).
1) Convective storms
A convective storm in a radar reflectivity image is a region of high intensity (Lakshmanan et al. 2009). It can be identified by thresholding the radar image based on a physically reasonable value such as 30 or 35 dBZ (Dixon and Wiener 1993). However, determination of the threshold is difficult because a lower threshold will result in false mergers of storms and a higher threshold will miss storms that are growing. It would be better to use a threshold that varies throughout the image rather than a single global threshold.
In the present study, the method used identifies convective storms using two thresholds, and , which are set to 30 and 35 dBZ, respectively. Figure 3 illustrates how to identify convective storms using these two thresholds. First, the method segments all low-intensity regions (LIRs) using . An LIR is a region of connected pixels where each pixel is larger than . Then, for each LIR, the method identities high-intensity regions (HIRs) within each LIR using . An HIR is a region of connected pixels where each pixel is larger than . If a specific LIR contains more than one HIR, it may correspond to a cluster of convective storms (see in Fig. 3). The method outputs all HIRs within the LIR as convective storms. If an LIR contains less than or equal to one HIR, it may correspond to a cellular storm (see in Fig. 3) or a growing storm (see in Fig. 3). The method outputs the LIR as an identified convective storm.
A region tree structure will simplify the above-mentioned storm identification method because reflectivity regions and their overlapping relationships have been recorded in the region tree structure. When a radar reflectivity image P is represented using a region tree structure T, each convective storm in P corresponds to a subtree structure that meets one of the following two conditions: 1) the intensity of is equal to 30 dBZ and or 2) the intensity of is equal to 35 dBZ and . Figure 4 provides examples of subtrees that each meet one of these two conditions, where subtrees 3 and 4 meet condition 2, and subtree 5 meets condition 1. To identify convective storms, the method in the present study first finds all subtrees that meet one of the above-mentioned two conditions, and then combines regions within subtrees to form convective storms.
2) Storm cells
In a radar reflectivity image, a storm cell is a high-intensity region whose intensity should be larger than 40 or 45 dBZ (Dixon and Wiener 1993). The procedures for identifiying storm cells are similar to those used for identifying convective storms. Two thresholds, 40 and 45 dBZ, are used to identify storm cells. Let T be a region tree structure of image P. The method first finds all that satisfy one of the following two conditions: 1) the intensity of is equal to 40 dBZ and or 2) the intensity of is equal to 45 dBZ and . Then, the storm cell corresponding to is obtained by combining all regions within .
3) Stratiform regions
A region in a radar reflectivity image is identified as a stratiform region if at least 70% of the region contains reflectivity between 20 and 40 dBZ (Gagne et al. 2009). When a radar reflectivity image P is represented using a region tree structure T, each stratiform region in P corresponds to a subtree that meets the following two conditions: 1) the intensity of is equal to 20 dBZ; and 2) , where and are the sums of areas of regions with intensities of 40 and 20 dBZ, respectively.
4) Spatial relationships among storm objects
After identification of storm objects, their spatial relationships are then identified. Let be a tree structure of P. Storm objects identified from P are recorded in three datasets: , , and , where each storm object corresponds to a subtree structure of T. Suppose and are two storm objects and their corresponding subtree structures are and , respectively. Two kinds of spatial relationships between and are defined as follows: 1) contains if ; and 2) and are neighbors if a storm object exists, and , so that and .
Based on the above-mentioned two definitions of spatial relationships, the spatial relationships among the storm objects can be easily determined. These spatial relationships include contained relationships between and , neighbor relationships between , contained relationships between and , and neighbor relationships between . Once the spatial relationships among the storm objects are obtained, the storm objects could also be organized using an object tree structure. For example, a region tree is shown in Fig. 4a and storm objects in the tree are labeled using rectangular boxes (Fig. 4a). An object tree structure in Fig. 4b is used to represent spatial relationships among storm objects. Each node in the object tree denotes a storm object and each link in the object tree presents contained relationships between two storm objects. Two storm objects are neighbors if they are contained together in the same storm object.
c. Tracking of storm objects
Suppose and are two region trees of radar reflectivity images and at successive scans, and and are sets of identified storm objects from images and , respectively. The tracking method first employs an iterative tree matching algorithm to obtain a mapping M from nodes in to nodes in ; next, it uses the mapping M to associate storm objects in and .
The tree matching algorithm is the key to this tracking method. Before the tree matching algorithm begins, regions corresponding to nodes in are first extrapolated based on their velocities. Section 3d provides the process of velocity estimation. Nodes in and are divided into subsets according to their depth, such that and , where denotes a subset of nodes in , and each node in has a depth of k. Let m be the minimum depth of and , and be the matching result of nodes in and . The iterative tree matching procedures are described as follows:
In the first iteration (), match nodes in and . Since contains only and contains only , .
- In the kth iteration, the obtained mapping result is denoted as and are the numbers of nodes in and , respectively, while denotes a pair of associated nodes. For each pair of , their child node sets, which are denoted as and , respectively, are found from the tree structures.
- For each node pair between and , their overlapping ratio is calculated as is the area of the region corresponding to the node υ and is the area of the overlapped region of nodes and . If is larger than a threshold of minimum overlap ratio , which is determined in the parameter sensitivity test, then nodes and are associated. The node pair is preserved in .
- Nodes matched in step 1 are removed from and . The remaining nodes in the sets are matched using a combinatorial optimization method (Dixon and Wiener 1993). The objective function we use is and are the indices of the nodes in and , respectively. Term is the cost function between nodes and , and is given by is the distance between the predicted centroid location of and the location of , and is the difference in area between nodes and . The cost function is similar with that adopted by Dixon and Wiener (1993). Selection of the cost function is based on the assumption that correct storm associations will join storms with similar characteristic (e.g., area) and will produce short storm tracks.
For matched regions in step 2, region velocities are calculated as the difference between centroid locations divided by the time interval between consecutive radar frames. The matched regions are preserved in if their corresponding region velocity is less than a threshold of the maximum region velocity , which is determined in a parameter sensitivity test.
If is empty or , then and the process is terminated. Otherwise, let and go to step 2.
Once the mapping M from nodes to nodes is obtained, matching of storm objects within and will be simple. Recall that each storm object corresponds to a subtree. Let and be two storm objects of the same kind (e.g., they are both convective storms). If any node pair exists, where and , then the storm objects and are associated.
Figure 5 illustrates the iterative process of the matching algorithm using an example. Figures 5a–c show two successive radar reflectivity images and their overlapping images. The tree structures corresponding to Figs. 5a,b are given in Fig. 5d. Nodes in trees and their corresponding regions are labeled with the same number. The matching results are shown in Fig. 5d using dashed lines. In this example, node pairs (,), (,), (,), and (,) are associated using the area-overlapping method; other pairs are matched using the combinatorial optimization method. Note that overlaps both and ; it is therefore considered that splits into and . When compared with TITAN and ETITAN, the tree-based algorithm also adopts a similar area-overlapping-based method to detect splits and mergers; however, the tree-based method detects splits and mergers between regions rather than between storm objects.
Figure 6 shows an example of storm object matching. Given two tree structures at times and , the dashed arrows indicate tree matching results and rectangular boxes represent identified storm objects. Storm object matching results are shown using solid arrows. Two storm objects are associated if any nodes exist that belong to them and are matched.
d. Nowcasting of storm objects
In this paper, we introduce a tree-based nowcasting method (see Fig. 7). When a storm object (see Fig. 7a) is represented using a region tree (see Fig. 7b), the tree-based nowcasting method initially extrapolates each region within the tree structure based on its velocity and the forecast time (see Fig. 7c), and then combines all the extrapolated regions in the tree to obtain the nowcasted storm object (see Fig. 7d). This tree-based nowcasting method is capable of nowcasting different kinds of storm objects simultaneously, because the method does not nowcast storm objects directly, but it nowcasts all regions that are components of all storm objects.
One problem with this tree-based nowcasting approach lies in estimating the velocity of each region within the tree structure. In this work, two methods are combined to estimate region velocities. One is an area-based method (Han et al. 2009) that initially estimates the motion vectors within the entire radar reflectivity field using the TREC method (Rinehart and Garvey 1978) and then defines the velocity of a region to be the average motion vector within it. Parameters used in the TREC method (Table 3) were adopted from the work of Rinehart and Garvey (1978).
The second one is a linear regression method (Johnson et al. 1998) that estimates the motion vector of a region based on its centroid track, defined as , where is the present location is the location one time step in the past, and so on, and is the length of the region track. Let be a measure of time. A linear regression is performed between and , yielding an equation of a straight line, . The slope of the line is the velocity of the region in the X direction (). A similar procedure is performed to obtain , the velocity of the region in the Y direction. Thus, the estimated motion vector of the region is . The linear regression method needs to handle region splits and mergers (Dixon and Wiener 1993). When region tracks merge, the merged track is a weighted combination of the parent centroid tracks, where the weights are the ratio of the region area for each parent to the sum of the areas of all parents. If a region splits into two or more regions, each new region will inherit the centroid track of the old region. Moreover, the inherited centroid track should add an offset, which is calculated as , where is the centroid position of the new region and is the average centroid position of all new regions.
For each region in the tree, if its area is larger than a threshold , then its velocity is estimated using the area-based method, otherwise its velocity is estimated using the linear regression method. The threshold is set to 250 km2 in a parameter sensitivity test (see section 4c).
In the region tree structure, there are regions that are smaller than and these do not have centroid tracks. Velocities of these regions cannot be estimated using the above-mentioned two methods. In this case, their velocities are estimated according to their ancestors in the tree. In practice, velocities of regions are first estimated in the tree using the area-based and linear regression methods. Then, the method traverses the region tree from root to leaves. If the velocity for a region in the tree is zero, but the velocity of its parent is not zero, then the velocity of its parent is assigned to it.
4. Test and results
a. Case example
This section shows individual case results of the proposed tree-based storm-tracking method. Figure 8a shows the case of an MCS at 2011 UTC (from case 3 in Table 1). Figure 8b shows the corresponding region tree structure of Fig. 8a. Figure 8c shows identified storm objects that are organized using a tree structure. Links between storms and storm cells, and links between storms and stratiform regions are represented using blue lines and blue dotted lines, respectively. Figure 8c indicates that the proposed method can identify different kinds of storm objects simultaneously and preserve the spatial relationships between them. Figure 8d presents identified storm objects using their contours; here one can see that storm objects are identified using adaptive thresholds. For example, the largest linear convective storm has an intensity larger than 35 dBZ, and the small convective storms on its left have intensities larger than 30 dBZ.
Figure 9 shows tracks of identified convective storms and cells from case 3 at two times, 2006 and 2011 UTC. Dots and black lines indicate the centroid positions and tracks of the storm objects, respectively. In this case, the largest linear convective storm moved from west to east in the radar field, split at 2011 UTC. Storm cells moved from southwest to northeast and exhibited many merges. As we can see, convective storms and storm cells within them are tracked simultaneously.
Figure 10 shows an example of results of the tree-based nowcasting method. Figure 10a shows an 18-min nowcast of the radar reflectivity image in Fig. 8d. Figure 10b shows the actual radar reflectivity image at the nowcast time (2029 UTC). Storm objects in both the nowcast and the actual radar reflectivity images are indicated using their contours. As we can see, the tree-based nowcasting method could nowcast the entire radar image, as well as nowcast different kinds of storm objects simultaneously. A storm object can be regarded as a group of regions with different intensities overlaid together. The tree-based method extrapolates each region within the storm object individually, and then combines all extrapolated regions to form the nowcasted storm object. Figure 10a also shows that there are situations where the internal part of storms could move out of the boundary of the storm. In the future, this problem needs be solved to obtain better nowcasting results.
b. Quantitative evaluation of storm-tracking performance
To evaluate the performance of the tree-based storm-tracking algorithm, it was compared with two other storm-tracking methods: SCIT and ETITAN. Parameters of the three storm-tracking methods are shown in Table 3. It should be mentioned that a direct comparison between different storm-tracking methods is impossible because these methods track different storm objects (Wilson et al. 1998). For example, SCIT and ETITAN are designed to track 3D storm cells and 3D convective storms, respectively. Therefore, the evaluation compares different storm association heuristics applied only in these specific storm-tracking methods. To perform a fair comparison, different storm-tracking techniques were run against convective storms identified by the tree-based storm identification approach.
In this evaluation, storm cases listed in Table 1 were categorized into two classes: isolated convective storm (cases 1, 2, 5, and 6) and MCS (cases 3, 4, 7, and 8). Convective storms in these cases were identified using the tree-based storm identification method (see section 3b). Then, identified convective storms were tracked using the three tracking methods (SCIT, ETITAN, and tree based). Note that the tree-based method is able to identify three kinds of storm objects, but only convective storm detections were applied during the evaluation. Figure 11 shows truncated area histograms of identified storm objects. Storm cells were not evaluated because they exhibit similar area distributions with those of convective storms but their number is much fewer than convective storms. Stratiform regions were not applied because they are not considered by conventional tracking algorithms, such as SCIT and ETITAN.
1) Evaluation of storm assignment
In the first test, the three storm-tracking methods were evaluated by comparing storm association results with those determined by humans. To overcome the subjectivity of human analysis, we gained the cooperation of three operational meteorologists in the Tinjin Meteorological Bureau—Jia Huizhen, Xie Yiyang, and Luo Kai—for the test. These three people produced three unique storm assignment datasets. Based on these datasets, a true assignment dataset was calculated; an assignment was called “true” only if it was contained in at least two assignment datasets.
The three storm-tracking methods produced three automated storm assignment datasets. By comparing each dataset of automated assignment with the true assignment dataset, we created three kinds of assignments: an assignment was called a “hit,” a “miss,” or a “false alarm” if we found both in the true assignment dataset and the algorithm output dataset, if we found in the true assignment dataset but did not find in the algorithm output dataset, or if we found in the algorithm output dataset but did not find in the true assignment dataset, respectively. Let the number of hits, misses, and false alarms be X, Y, and Z, respectively; three kinds of scores are defined: the probability of detection (POD) = , the false alarm rate (FAR) = , and the critical success index (CSI) = .
Figure 12 shows the scores of the three storm-tracking methods. For all categories, the tree-based tracking method achieves the highest CSI values. The average CSI value for the tree-based tracking method is 78.16% as compared to 74.88% and 74.71% for SCIT and ETITAN, respectively. However, the tree-based method has a relatively higher average FAR value of 11.98% than either SCIT or ETITAN with FAR values of 9.01% and 7.50%, respectively. This high FAR value can be attributed to incorrect storm associations when a storm splits or merges. To get the maximum CSI value, the parameter of the minimum overlap rate is set relatively low. Recall that the tree-based tracking method uses the region overlapping method to identify region splits and mergers (see section 3c). A lower would result in more false mergers and splits during storm tracking.
Figure 12 also shows that the tree-based tracking method achieves similar CSI scores (78.64% for isolated storms, 77.68% for MCS, and 78.16% for average) for all categories. In contrast, SCIT has a higher CSI value than ETITAN for isolated storms, while ETITAN achieves better performance than SCIT for MCSs. On average, SCIT and ETITAN achieve almost the same performance.
2) Evaluation of storm tracks
In the second test, an objective storm-tracking evaluation method that had been proposed by Lakshmanan and Smith (2010) was carried out. This method evaluates a storm-tracking algorithm in terms of three factors (the durations of the tracks, the linearity of the tracks, and the preservation of a storm attribute) of the storm tracks produced by that storm-tracking algorithm. In general, a storm-tracking algorithm that generates longer, more linear tracks where the storm attributes are relatively constant is better.
Once storm tracks produced by a storm-tracking algorithm are obtained, three statistics related to each track are initially calculated; they are the duration of the track (), the standard deviation of the area along the track (), and the root-mean-square error () of the centroid locations from their linear least squares fit. Then, the centroid tendencies of the above-mentioned three statistics are computed on all storm tracks:
Term is the median duration of tracks in the dataset.
The mismatch error is the mean on the tracks, which have a duration larger than .
The linearity error is the mean on all tracks with a duration larger than .
The evaluation method uses to measure the duration of storm tracks because this parameter is not affected by outliers (the tracks of duration are less than two frames). The statistic is a measure of the preservation of the storm area along the track, which is lower if there are fewer incorrect storm associations. The statistic is a measure of the linearity of a storm track, which is lower if the storm track is more linear.
Figure 13 shows score values of the three storm-tracking methods (SCIT, ETITAN, and tree-based) that ran on the eight cases listed in Table 1. For the isolated storm category, the differences in the score results are not very large when compared with the MCS category. On average, the tree-based algorithm produces the smallest mismatch error, the second longest median duration, and the largest linearity error. This indicates that the tree-based tracking method performs well on measures of mismatch error and median duration, but poorly on the measure of linearity error. The poor performance on linearity error may be caused by the use of the area-overlapping method during the tracking process. The area-overlapping method does not consider centroid distance between regions and tends to associate regions with similar areas. As a result, the tree-based tracking method produces tracks with small mismatch errors and relatively large linearity errors.
Figure 13 also shows that no tracking algorithm performs well on all three measures. For example, SCIT performs worst in terms of mismatch error: it achieves the largest mismatch error on average; ETITAN has the worst performance in terms of median duration.
3) Evaluation of nowcasts
The third evaluation used the three storm-tracking methods (SCIT, ETITAN, and tree based) to create short time nowcasts of convective storms and then compared the nowcasts with actual convective storms at forecast time. In total, 442 frames of radar reflectivity images from eight cases listed in Table 1 were evaluated during this test. The verification data were from a 0.5° radar reflectivity image at the forecast time. Convective storms in the verification field, which were identified using the tree-based storm identification method, were used as reference data. A contingency table based on a grid-to-grid comparison of the nowcasts and the reference data was used to calculate the POD, FAR, and CSI. In the calculations, the resolution of the grid was 1 km × 1 km in horizontal. At each grid point, a success/hit occurred when both the nowcast and the reference data contained storm events, a miss occurred when the reference data contained a storm event but the nowcast did not, and a false alarm occurred when the nowcast contained a storm event but the reference data did not.
Figure 14 shows score values of the three nowcasting methods for 18- and 30-min nowcasts. The values of POD, FAR, and CSI are the average values for all the cases belonging to each category. For 18-min nowcasts, the tree-based nowcasting method achieves an average CSI of 29.02%, which is higher than both 24.98% for SCIT and 24.44% for ETITAN. Figure 14 also shows that large differences occur for the category of MCS: as compared with SCIT and ETITAN, the tree-based method achieves improvements of CSI values by 6.05% and 5.81%, respectively. For 30-min nowcasts, the tree-based method achieves an average CSI of 18.46% as compared with SCIT (CSI = 15.96%) and ETITAN (CSI = 15.51%).
As described in the previous section related to the nowcasting of storm objects, the tree-based nowcasting method is a combination of the area-based nowcasting method (e.g., ETITIAN) and the linear regression–based nowcasting method (e.g., SCIT). This experiment demonstrates that the combination of these two methods performs better in predicting the future positions of storm objects.
The CSI is a measure used to predict a storm’s location. The CSI alone is insufficient in this evaluation, because in addition to location nowcasting, the tree-based storm nowcasting method has the ability of nowcasting the internal evolution of a storm, but the calculation of the CSI is not affected by storm’s internal evolution. Therefore, to evaluate this ability, other verification measurement is needed to supplement the CSI.
In this evaluation, we introduced another scalar measurement, mean absolute error (MAE). The MAE is defined as
where and are the ith of N pairs of nowcasts and reference data, respectively. In general, if a nowcasting algorithm could predict the internal evolution of a storm, it should produce a small MAE value.
Figure 14d shows MAE values of the three nowcasting methods. For 18-min nowcasts, the tree-based nowcasting method achieves an average MAE of 4.63 dBZ, which is lower than 5.22 dBZ for SCIT and 4.80 dBZ for ETITAN. For 30-min nowcasts, the tree-based method also produces a smaller MAE value (MAE = 4.86 dBZ) than both SCIT (MAE = 5.56 dBZ) and ETITAN (MAE = 5.44 dBZ). This result is excepted because a storm always exhibits internal evolution. Nowcasting algorithms such as SCIT or ETITAN that ignore these evolutions would produce relatively large MAE values. Figure 14d indicates that the tree-based nowcasting method has the ability to nowcast the internal evolution of a storm.
4) Running time evaluation
The running times of storm-tracking methods are very important because these methods need to run online. The three storm-tracking methods were implemented in MATLAB 2014b on an Intel Corel2 1.86-GHz processor with 2 GB of RAM. The average storm-tracking and nowcasting running times for the three storm-tracking methods and the average tree-based storm identification running time are shown in Table 4. The average storm identification running time (965.25 ms) is much larger than the average tracking and nowcasting running times for the three storm-tracking methods. However, when compared with the sampling rate of WSR-88D radar (5 or 6 min), the average storm identification running time for the tree-based method is acceptable.
c. Determination of parameters
The performances of the three storm-tracking methods are affected by their parameters (Table 3). The three storm-tracking methods ran on four cases (1–4) listed in Table 1 using different combinations of parameters to maximize CSI values and determine the optimal parameters.
Three parameters of the tree-based storm-tracking method were tested: two parameters—the minimum overlap rate and the maximum region velocity—were used to associate region trees; in addition, the region area threshold was used in motion vector estimation. The first two parameters were tested using the CSI value of the storm assignment. The third parameter was determined according to the CSI value of an 18-min nowcast. Table 5 shows the optimal minimum overlap rate lies between 0.05 and 0.2 and the optimal maximum region velocity lies between 100 and 200 km h−1. Figure 15 shows the optimal region area threshold lies between 200 and 300 km2.
Performance of SCIT is influenced by two parameters: 1) the maximum search distance used for associating previous storms and current storms, and 2) the number of previous locations used in motion vector estimation. The CSI value of the storm assignment was employed to test these two parameters. Table 6 shows the CSI values of the storm assignment with respect to different combinations of these two parameters. As we can see, the optimum maximum search distance lies between 10 and 20 km and the CSI value of the storm assignment is not sensitive to the number of previous locations.
The parameter of the minimum overlap rate was tested for the ETITAN algorithm. Figure 16 shows the variation of the CSI values related to the storm assignment with respect to the minimum overlap rate. The optimum value of the minimum overlap rate is about 0.1. Note that this value is different from that (0.5) suggested in Han et al. (2009).
This paper proposed a tree-based storm-tracking method that aims to identify, track, and nowcast different kinds of storm objects (stratiform regions, convective storms, and individual cells) based on radar weather data. When compared with previous object-based storm-tracking methods, the tree-based storm-tracking method is able to 1) identify, track, and nowcast different kinds of storm objects simultaneously; 2) recognize spatial relationships between these storm objects; and 3) nowcast the internal movements and changes within storm objects. In addition, the storm object identification process is independent of the storm object tracking and nowcasting processes, so storm object tracking and nowcasting are not affected by storm object identification.
The tree-based storm-tracking method describes a radar-based reflectivity image using a region tree structure where each node corresponds to a radar reflectivity region and each link denotes an overlapping relationship between two radar reflectivity regions. When a region tree structure of a radar reflectivity image is constructed, regions in the tree structure are clustered to form different kinds of storm objects. The output of storm identification is a tree structure of storm objects that represents the spatial relationships among the storm objects.
During the tracking stage, regions in the tree are first associated using an iterative tree matching method that combines many existing association techniques, such as the cost function method applied by TITAN, the area-overlapping method used by ETITAN, and the centroid project method applied by SCIT. After that, two storm objects are matched if the regions contained by them were associated in the tree matching method.
In the nowcasting stage, regions in the tree are first individually extrapolated according to their velocities, and then combined to form nowcasted storm objects. The velocities of the regions in the tree structure are estimated using two methods: If a region is larger than 250 km2, then its velocity is estimated using the Tracking Radar Echoes by Correlation (TREC) algorithm (Rinehart and Garvey 1978), otherwise its velocity is estimated using the linear regression method (Dixon and Wiener 1993; Johnson et al. 1998).
Using eight cases of different kinds of convective storm events, performances of the three storm-tracking methods (SCIT, ETITAN, and tree based) were quantitatively evaluated in three tests. The first test evaluated a storm-tracking method using CSI scores of storm association. The results indicate the tree-based method achieved better performance than the other two methods and had a CSI score of 78.16% as compared to 74.88% for SCIT and 74.71% for ETITAN. The second test evaluated a storm-tracking method using statistics of storm tracks produced by the storm-tracking method. The results show that the tree-based storm-tracking method performs well on two out of three measures; it produced long tracks with small mismatch errors. The third test evaluated a storm-tracking method through CSI and MAE scores of storm nowcasts. The results of the third test show that the tree-based method outperforms both SCIT and ETITAN in nowcasting positions of storms and also has the ability to nowcast the internal evolution of a storm. For 18-min nowcasting scores, the tree-based method achieved an average CSI of 29.02%, compared with 24.98% for SCIT and 24.44% for ETITAN. Also, the tree-based method produced an average MAE of 4.63 dBZ, compared with 5.22 dBZ for SCIT and 4.80 dBZ for ETITAN.
Although the tree-based storm-tracking method performs better than the other methods (SCIT and ETITAN), many refinements are still needed. For example, the region tree structure used in this work is able to represent only two-dimensional radar images; however, the three-dimensional structure of a storm object is very important in weather nowcasting. Therefore, a three-dimensional region tree structure representation method could be explored in the future.
Funding for this work was provided by the China Special Fund for Meteorological Research in the Public Interest (GYHY201406004) and the Natural Science Foundation of Tianjin, China (14JCYBJC21800).