Toward automated hail disaster weather recognition based on spatio-temporal sequence of radar images

2.1 Situation and disaster in Hunan Province

According to incomplete records of hailstorm weather in Hunan Province, at around 10:10 p.m. on the evening of March 19, 2013, Jingzhou Miao and Dong Autonomous County in Hunan Province were suddenly hit by a very large hailstorm for nearly 30 min. A total of more than 14,800 households and 44,700 people were affected in the county, with direct economic losses of about 163 million yuan. Among them, 45,590 houses were damaged; 159,000 mu of crops were affected; 35 enterprises’ plants, equipment, raw materials, and finished products were severely damaged; and some electric power, water supply, transportation, water conservancy, commerce, storage, and other facilities were damaged. From March 4, to March 6, 2018, wind and hail damage occurred in some areas of Hunan Province. The wind and hail disaster caused a total of 173,000 people to be affected in Changsha, Zhuzhou, and other cities, with more than 710 people urgently relocating, 59 houses collapsing, more than 7,000 houses damaged to varying degrees, and direct economic losses of RMB 110 million. On May 11, 2021, some townships in Xintian County, Hunan Province, were hit by hailstorms, damaging tens of thousands of mu of baked tobacco in Xinxu Township, Xinlong Township, and Taoling Township. More than 6,700 mu was extinguished among them, with economic losses estimated at more than 35 million yuan.

To sum up, it is very important to minimize losses by forecasting in advance before the occurrence of extreme hailstorms. However, due to the sudden, localized, and short life cycle characteristics of the hailing process, its forecasting effect has been unsatisfactory for a long time. Even if it is only half an hour in advance of the proximity forecast, its progress is very slow, and the most serious difficulty lies in the identification of hail echoes. To address this problem, this article adopts the sequence data from radar stations to understand the weather trend and ensure the hail identification results’ correctness as much as possible.

2.2 New-generation Doppler weather radar

The data used in this study are the raw Doppler meteorological data acquired by the new-generation Doppler weather radar. The data include Julian time, body sweep mode, reflectivity factor, radial velocity, and spectral width. Among them, the reflectivity factor (in “dBZ”) can identify cumulus, stratus, and cumulonimbus cloud systems and weather features such as fronts and bob lines. Meteorological observers can determine the intensity of the echoes by the reflectivity factor and determine the weather of rain, snow, wind, lightning, hail, etc. Its characteristics can provide an important scientific guidance basis for weather forecasting [24]. In general, the higher the reflectivity intensity, the greater the possibility of strong convective weather. The reflectivity range is usually divided into 15 levels to facilitate observation. Each color corresponds to a range of reflectivity intensity, and black indicates invalid echoes.

In the specific representation of radar reflectivity images, the spatial distribution of radar reflectivity can be described according to the distribution of colors. Among them, green indicates the weaker reflectivity, yellow indicates the medium reflectivity, red indicates the stronger reflectivity, and purple indicates the strongest reflectivity. Therefore, we can judge the size of radar reflectivity according to the color to analyze the weather situation. For example, if the entire image is green or yellow, it indicates that the reflectivity in the region is relatively weak and the weather is relatively good. If there are red, blue, or purple areas in the image, it indicates that the reflectivity in the region is strong, and there may be weather phenomena such as heavy rain or hail. Additionally, the time series changes of the radar reflectivity image can be used to describe the evolution of weather. For example, suppose the color in the image appears to have a significant change. In that case, the weather conditions are changing. The reflectance intensity corresponds to the color, as shown in Figure 2.

Figure 2

Left: correspondence between basic colors of radar reflectivity map and reflectivity intensity; right: reflectance image.

The radar reflectivity image is derived from the reflectivity factor in the basic product, which takes the location of the radar as the origin, and the color taken at each coordinate point in the coordinate system is used to reflect the reflectivity intensity of the radar echo at the corresponding spatial location, as shown in Figure 2. During the processing of radar echo data, the reflectivity of the corresponding elevation angle is read from the radar base data. The point with the highest reflectivity is identified, and it is determined whether the reflectivity at that point is greater than 45 dBZ. If the reflectivity intensity in the area exceeds 45 dBZ, it can be initially concluded that hail may be present in the region. However, this is only a preliminary assessment and requires comprehensive analysis and judgment in conjunction with other factors.

3.1 Model architecture

In hailstorm weather image data processing, the newly developed Doppler weather radar is used for individual scanning. After pre-processing and preprocessing the collected data, the data were augmented through deep convolution generative adversarial networks (DCGANs), and the spatio-temporal data were generated through relatively formatted images. Then, LSTM-C3D was used to carry out feature vector extraction. In the information coding part, gradient highway unit (GHU) high-speed channel unit and causal LSTM computing unit are used to select and fuse long-time and short-time feature information, encode the complete sequence information spatio-temporally and map the long sequence to the low-dimensional feature space. The decoded information was further extracted using an improved 3D ConvNet. Finally, it was mapped to specific categories through global pooling and full connection. LSTM-C3D was used to analyze the characteristics of the data, and then, the candidate regions were divided. Based on this, a network structure with the focus mechanism as the core was constructed to obtain the key information of the hail weather characteristics to improve the accuracy of the model. Finally, to further improve the accuracy rate of hail occurrence, VGG, ResNet, and AlexNet were selected in this study to construct the triple-branch structure. The voting network was used to improve the universality of the model. The RPN regional recommendation algorithm was used to process the feature data in the regional convolutional network. After further obtaining the feature map, softmax was used to classify whether there was hail [25]. The structure diagram of this model is shown in Figure 3.

Figure 3

Structure of the LSTM-C3D model.

In hail identification problems, it is usually necessary to provide labeled learning sample number sets to train the model. When the classifier is equipped with classification capability after continuous training, it can accurately identify the presence or absence of hail using the latest input Doppler radar data. Figure 4 shows the meteorological image information processing flow of the whole hail using LSTM-C3D.

Figure 4

Schematic of the overall system.

3.2 Causal LSTM

To improve the accuracy of hail identification, this task uses sequence recognition to memorize the trend of radar data changes. Compared with the traditional ST-LSTM [26], the causal LSTM [27] can better solve the problem of recursive depth by adding a nonlinear layer, thus increasing the network depth and memory capacity, and its cell structure is shown in Figure 5. A causal LSTM cell contains dual memory, temporal memory C t k , and spatial memory M t k , where the subscript t indicates the time step and the superscript k indicates the kth hidden layer in the stacked causal LSTM network. The current temporal memory depends directly on its previous state C t − 1 k controlled by the forgetting gate f t , the input gate i t , and the input modulation gate g t . The current spatial memory M t k depends on the deep migration path M t k − 1 . For the bottom level ( k = 1 ), the uppermost spatial memory at ( t − 1 ) is allocated to M t k − 1 . The causal LSTM differs significantly from the original spatio-temporal LSTM in that it uses a cascade mechanism in which the spatial memory is a function of the temporal memory of another set of gate structures. This newly designed cascade memory outperforms the simple cascade structure of the spatio-temporal LSTM because of the significant increase in recursion depth along the spatio-temporal transition path.

Figure 5

causal LSTM cell structure.

The causal LSTM equation is expressed as follows, where * stands for convolution; σ is a sigmoid function; square brackets indicate a cascade of tensors; round brackets indicate a system of equations; x denotes the input vector; H denotes the output vector; C denotes the state vector; and g , i , and f are the valve vectors, representing the weights of the input and output information. C is added to the formulae to jointly determine the opening and closing of the gate, where C in the previous step of the first formula and C in the latter stage of the third formula form a cascade through a recursive relationship. In addition, tanh is adopted in the last step to add non-linearity:

In addition to short-term video dynamics, causal LSTM tends to suffer from difficulties in gradient backpropagation in the long term. Due to longer transitions, the temporal memory may forget obsolete features. Theoretical evidence shows that highway layers can efficiently pass gradients in very deep feedforward networks and prevent long-term gradients from disappearing rapidly in recurrent networks relying on GHU. Its cell structure is shown in Figure 6. The gradient highway unit (GHU) structure is introduced in the encoder’s first and second stacked layers. The input of the GHU is the output of the current lower layer and the input of the GHU of the previous moment, i.e., the GHU connects the current moment and the input of the previous moment, which can effectively solve the problem of gradient disappearance.

Figure 6

GHU cell structure.

The equation for GHU is expressed as follows, where W represents the convolution filter. Z t is the switch gate because it is capable of adaptive learning between converted input P t and hidden state Z t :

Considering the low frequency of hail in the actual meteorological process, this directly leads to the low amount of data available for training. Additionally, we made a tiny area crop of the original radar map during the data processing. Therefore, in the encoder part, the number of layers of vertical stacking is reduced from the stacking approach, which, on the one hand, relieves the pressure of the network too deep gradient difficult to backpropagate and, on the other hand, reduces the number of parameters of the network, thus speeding up network training and online extrapolation. Since small areas of radargrams are used as input, the shallower stacking is sufficient to extract spatially valid information. This approach prevents the network from over-fitting to changes in the values of single-pixel points and effectively focuses on the overall trend changes. Meanwhile, for the causal calculation unit of the encoder, the model is designed with two scales of calculation: large convolution kernels for capturing the overall radar echo change information and small convolution kernels for capturing the local echo generation and extinction. In particular, for the high-echo information that often accompanies hailstorms, the rate of change of its generation and extinction is correlated with the presence of hailstorms so that the smaller convolution kernels can be used to capture more of this detailed change sensitively.

Similarly, subsequently, in the decoder part, the spatio-temporal information of both scales is stored in the encoder. Thus, the corresponding feature data are effectively extracted in the decoder and, after 3DCNN block, finally mapped to categories for classification using global pooling and full concatenation.

3.3 C3D

In the design of the decoder, inspired by the design of the C3D network [28], the feature extraction structure design of the network is similar to that of the C3D network but has been improved on the basis of the original. There are eight convolution layers (Conv1, Conv2, Conv3a, Conv3b, Conv4a, Conv4b, Conv5a, and Conv5b), five pooled layers (Pool1, Pool2, Pool3, Pool4, and Pool5), three fully connected layers (Fc6, Fc7, and Fc8), and one softmax classification layer, in which each convolution layer and fully connected layer is connected with ReLu nonlinear layer. The fully connected layers (Fc6, Fc7, and Fc8) with the ReLu nonlinear layer use dropout technology to prevent overfitting in network training. The number of convolution cores in eight convolution layers (Conv1, Conv2, Conv3a, Conv3b, Conv4a, Conv4b, Conv5a, and Conv5b) is 64, 64, 128, 256, 256, 512, 512, 512, and 512 in turn, and the size of convolution cores is (3, 3, 3). The convolution step is (1, 1, 1). Among the five pooling layers (Pool1, Pool2, Pool3, Pool4, and Pool5), except the first pooling layer (Pool1), the sampling window size is (2, 2, 1), and the sampling step is (2, 2, 1). The sampling window size of the other pooling layers (Pool2, Pool3, Pool4, and Pool5) is (2, 2, 2), and the sampling step is (2, 2, 2), in order to ensure that the time-domain information can be retained as early as possible. After several convolution pooling operations, the features are sampled down to obtain more globally representative features. Each fully connected layer (Fc6, Fc7, Fc8) contains (8,192, 4,096, and 4,096) output units. The last layer of the network is the softmax classification layer, and the output is the hail zone and non-hail zone. For the adaptation task, the network input channel of C3D is changed to 3, and the category output dimension is set to 2. Similarly, subsequently, in the decoder part, the spatio-temporal information of both scales is stored in the encoder. Thus, the corresponding feature data are effectively extracted in the decoder and, after 3DCNN block, finally mapped to categories for classification by means of global pooling and full concatenation.

3.4 Attention mechanism

Most hail images are based on intersection-over-union (IoU) matching method to determine positive and negative samples, which only obtains the features that make up the larger part of the data while neglecting the identification and localization of other features that are obvious and have an impact factor, resulting in the inability to understand the determination of hail generation fully. In addition, some data features disappear after deep learning-based multi-layer convolution, resulting in a shallow network structure in the recognition model that cannot extract semantic information. To address these issues, an attention-focused mechanism is used to overcome these problems and obtain more accurate hail data features. In this article, we design the network structure of the spatial attention mechanism as shown in Figure 7, using the residual layer for pooling operation to reduce the image feature vector dimension and then extracting the image spatial feature information through a convolution kernel of 3 × 3 size and step size (padding) of 1. This information is processed by the sigmoid function and weighted with the pooled image feature vector information in the scale layer.

Figure 7

Structure of attention mechanism.

3.5 Voting networks

Voting networks are typically used to process the output results of multiple classifiers in order to obtain more accurate classification results. In the voting network, the data commonly used are the prediction results of multiple classifiers on the same input data. For example, the same image can be fed into multiple different CNNs, each network outputting a probability distribution for a class, and these probability distributions are then passed to the voting network for integration.

In this article, we used AlexNet, ResNet50, and VGG deep learning neural network benchmark models to build a three-branch structural fusion network model, the structure of which is shown in Figure 8. Since the size of the convolutional kernel is an important hyperparameter, it determines the size of the receptive field of the convolutional layer when extracting image features. Different sizes of convolutional kernels can capture features of different scales, thus improving the feature extraction ability of the network. Specifically, using different sizes of convolutional kernels can help the network extract feature information of different scales, thereby enhancing the feature expression ability of the network. Therefore, the first branch uses a 1 × 1 convolution kernel to extract the candidate regions of the feature map. The smaller convolution kernel can capture more local details and improve the diversity and stability of the features. The second branch uses a 3 × 3 convolutional kernel to extract candidate regions from the feature map, which can capture relatively larger image features. The third branch uses two 3 × 3 convolutional kernels in series to extract candidate regions from the feature map, and its receptive field is equivalent to a 5 × 5 convolutional kernel. Specifically, the first 3 × 3 convolutional kernel can capture local features of the input image, while the second 3 × 3 convolutional kernel can combine these local features to form a more global feature representation. This not only reduces the number of model parameters and improves the computational efficiency of the model but also enhances the nonlinear expression ability of the network, improving the accuracy and robustness of feature extraction. The final model can adapt to multi-scale image feature extraction, increase the generality of the model, and improve the accuracy of the model determination.

Figure 8

Structure of the voting network.

For the voting network, we adopt the idea of a soft voting classifier, which takes the average probability of all models predicting samples to be a certain class as the standard, and the type corresponding to the highest probability is the final prediction result.

4.1 Experimental data and implementation

4.1.2 Data pre-processing

For a certain city/county/district with hail, it is necessary to find the radar station echo map with its correlation and slice it. The specific steps are as follows:

First, the actual distance between the radar station and the real station is calculated according to the latitude and longitude of the two stations, and the corresponding relationship between the real station and the radar station is determined.

Second, find the corresponding radar station data according to the hail data. For the case that the live data are one moment, the processing method is as follows: from the data of the radar station whose recording time is nearest to the live time, a total of five data are taken as the input sequence, i.e., a total of 30 min. If the live data are for a period of time, the handling method is as follows: first, determine the recording time of the radar station nearest to the start time and the end time, respectively. Start with the data of the radar station nearest to the end time and save five consecutive data. Then, start with the data of the radar station at the previous time and protect five consecutive data successively until the data of the radar station nearest to the start time are obtained. All data will be stored in a group of five in the way of radar station serial number – live station serial number coding to facilitate the subsequent data processing.

Third, the last data in the sequence are selected to read the reflectance data in the radar base data according to the 2.4 ∘ elevation angle and save it in the format of NumPy matrix. The specific value of the azimuth angle is held in the 0th column of the matrix, which is convenient for slicing the reflectance matrix in the next step. Figure 9 shows the evolution trend of radar reflectance in five frames of a radar detection area in Hunan.

Figure 9

Visualization of radar reflectivity evolution in a region over 30 min: (a) start, (b) 6th minute, (c) 12th minute, (d) 18th minute, and (e) 24th minute.

Finally, the program calculates the azimuth angle between the radar station and the live station and cuts the azimuth angle ± 3 0 ∘ with the live station as the center, 48 km at the proximal end and 48 km at the centrifugal end, i.e., a 60 × 96 area, and finds the point with the highest reflectivity among them. It is judged whether the point is greater than 45 dBZ; if so, the azimuth and radial bank numbers of the point are saved for slicing the other data and elevation angles in the sequence. For all six elevation angles, the azimuth angle ± 1 6 ∘ , 16 km at the proximal end, and 48 km at the centrifugal end, i.e., 32 × 64 , are selected as the center of the point. The sliced data of the six elevation angles are stacked together in the first dimension in a way to obtain data of size ( 6 × 32 × 64 ), finally (the length of the sequence is 5), and the data are given the label “Hail exists” label. The same process is applied to the radar data based on the live data (historical heavy precipitation data) without hail, and a label of “no hail” is given to this area data.

After forming the sample set with labeled hail data, the data are expanded using DCGANs based on Gaussian-distributed probability density random numbers as DCGAN input, and the sample set is partitioned according to 7:2:1 to form the train set, validation set, and test set.

4.2 Performance indicators

In the field of weather recognition and prediction, accuracy (A), precision (P), recall (R), and F 1 scores are commonly used as model performance metrics measurements, and this article also uses the aforementioned evaluation metrics to evaluate the performance of hail hazard weather recognition models. Applying the aforementioned evaluation metrics to the performance test of the hail hazard weather recognition model, two types of result images can be obtained: predicted results with hail and predicted results without hail. Among them, true positive (TP), false positive (FP), true negative (TN), and false negative (FN) are the key metrics used to describe the accuracy. Specifically, TP hail refers to the case where there is hail in the real situation and hail in the detection result at the same time, FP hail indicates the case where there is no hail in the real situation but hail in the detection result, TN hail refers to the case where there is no hail in the real situation and hail in the detection result at the same time, indicating that the model detects correctly, and FN hail indicates the case where there is hail in the real situation and hail in the detection result at the same time. The calculation formulas for accuracy, precision, recall, and F 1 score [29] are as follows:

Accuracy (A) is the percentage of correctly predicted results to the total sample and is given by the following equation:

Precision (P), also known as check accuracy, is for the predicted outcome, and it means the probability of the actual positive sample among all the predicted positive samples and its formula is as follows:

Recall (R), also known as check-all rate, is for the original sample, and it means the probability of a positive sample being predicted among all samples that are actually positive, and its formula is as follows:

F1 score is the summed average of precision rate and recall rate, which integrally reflects the model performance, and its formula is as follows:

4.3 Recognition results and analysis

4.3.1 Validation experiments

To analyze the prediction effect more intuitively, three sets of individual cases were selected from the test data set for analysis. The baseline model itself is used to generate future 10-frame predictions from five frames of historical radar maps. The experimental results are shown in Figure 10 (0–1 h).

Figure 10

Radar echoes at 6 min intervals (5 frames): (a) model input, (b) from the 6 to 30 min, (c) prediction from 6 to 30 min, (d) from the 36 to 60 min, and (e) prediction from 36 to 60 min.

As can be seen from the results in Figure 10, the proposed method can always capture the evolution of different intensity echoes and their locations more accurately, indicating that the model can “learn” and “understand” the development of strong convective weather and that the evolution of strong convective features learned from the training set can be effectively applied to new “future” scenarios with good generalization capability. Moreover, the model can effectively identify and remove the “noise” (radar clutter) from the input radar image (the noise in the upper-left part of the live image), thanks to the advanced filtering features inherent in the model’s convolution operation algorithm. It is worth noting that the model also consistently and accurately tracks the movement of the nascent echoes just entering the study area in the lower-right part of the figure. This is because the model encounters similar situations many times in the training set of massive data, so it can learn the movement trend of such radar echoes appearing on the boundary and give reasonable prediction results. However, the edge details of the different intensity echoes are gradually lost over time. The reason for this problem may be that the decoder network continues to generate the predicted output iteratively, and the difficulty of prediction is gradually increased so that the confidence in the expected value is gradually reduced and fed back to the input of the next moment, which eventually leads to the accumulation of the uncertainty of the prediction. After visualization, the performance of the image is reflected in the edge details, which slowly disappear and become blurred. Figure 11 shows the results of the prediction of the model versus the actual situation.

Figure 11

Model predictions vs actual conditions: a.1–c.1 for 0.5 h and a.2–c.2 for 1 h.

As seen in Figure 11, the model can accurately track the movement trend of convective storms in strong convective weather. The position and development situation of strong-echoes are closer to the corresponding actual radar observations, but there are some shortcomings as well. Compared with the real situation, the edge contours of different intensity echoes in the predicted radar echo maps are blurred and some details are lost. For example, if the real image shows several single thunderstorm echoes at similar locations, which are difficult to be distinguished by the model, they merge into a larger echo in the predicted image and eventually create a false alarm problem. One possible reason for this phenomenon is the inherent uncertainty in this weather prediction problem (nonlinear rapid evolution of strong convection), which propagates rapidly with time and makes it almost impossible to accurately predict the exact location of future echoes. Therefore, to reduce prediction errors, the models tend to produce fuzzy forecasts to improve the “hit rate” of the forecasts. These details are more “copied” from the original radar echoes, which are always changing dynamically, and this way of retaining details brings more false alarms (“false alarm rate” increases), and the position of the echoes it tracks is not accurate (the forecast position is somewhat ahead). In general, the 0.5 and 1 h proximity forecasts of the baseline model for different areas of strong convective echoes are in better agreement with the actual conditions.

4.3.2 Ablation experiments

Figure 12 shows the optimization results of the parameters. It can be seen that the accuracy of the model gradually increases with the increase of the training period and eventually tends to be stable due to the model having reached the fitting after the training period is greater than 800. By comparing different result models, it can be seen that after adding the attention mechanism and voting network structure, the accuracy of the model has been significantly improved under the premise of the same training cycle. At the same time, in the actual training, the final accuracy is not significantly improved by increasing the number of branch structures of the voting network, but the training time and prediction time are increased, so the final hail prediction system adopts the three-branch structure model. When the training period is 800, the three-branch structure voting network is used to compare the accuracy of the model, the baseline model, the additional attention focus structure, the additional attention focus structure, and the voting network. The accuracy is 91.3, 93.4, and 97.8%, respectively. The accuracy of the model can be improved by adding the attention focus mechanism. The voting network was added to enhance the performance of the model further.

Figure 12

Comparison of LSTM-C3D models with different modules added.