*Article* **Precipitation Nowcasting Based on Deep Learning over Guizhou, China**

**Dexuan Kong 1,2 , Xiefei Zhi 1,\*, Yan Ji 1, Chunyan Yang 2,\*, Yuhong Wang 3, Yuntao Tian 1, Gang Li <sup>4</sup> and Xiaotuan Zeng <sup>5</sup>**


**Abstract:** Accurate precipitation nowcasting (lead time: 0–2 h), which requires high spatiotemporal resolution data, is of great relevance in many weather-dependent social and operational activities. In this study, we are aiming to construct highly accurate deep learning (DL) models to directly obtain precipitation nowcasting at 6-min intervals for the lead time of 0–2 h. The Convolutional Long Short-Term Memory (ConvLSTM) and Predictive Recurrent Neural Network (PredRNN) models were used as comparative DL models, and the Lucas–Kanade (LK) Optical Flow method was selected as a traditional extrapolation baseline. The models were trained with high-quality datasets (resolution: 1 min) created from precipitation observations recorded by automatic weather stations in Guizhou Province (China). A comprehensive evaluation of the precipitation nowcasting was performed, which included consideration of the root mean square error, equitable threat score (ETS), and probability of detection (POD). The evaluation indicated that the reduction of the number of missing values and data normalization boosted training efficiency and improved the forecasting skill of the DL models. Increasing the time series length of the training set and the number of training samples both improved the POD and ETS of the DL models and enhanced nowcasting stability with time. Training with the Hea-P dataset further improved the forecasting skill of the DL models and sharply increased the ETS for thresholds of 2.5, 8, and 15 mm, especially for the 1-h lead time. The PredRNN model trained with the Hea-P dataset (time series length: 8 years) outperformed the traditional LK Optical Flow method for all thresholds (0.1, 1, 2.5, 8, and 15 mm) and obtained the best performance of all the models considered in this study in terms of ETS. Moreover, the Method for Object-Based Diagnostic Evaluation on a rainstorm case revealed that the PredRNN model, trained well with high-quality observation data, could both capture complex nonlinear characteristics of precipitation more accurately than achievable using the LK Optical Flow method and establish a reasonable mapping network during drastic changes in precipitation. Thus, its results more closely matched the observations, and its forecasting skill for thresholds exceeding 8 mm was improved substantially.

**Keywords:** precipitation forecast; nowcasting; deep learning; ConvLSTM; PredRNN

### **1. Introduction**

Precipitation nowcasting means forecasting precipitation with a lead time of 0–2 h, focusing more on mesoscale–microscale weather systems at a high spatiotemporal resolution [1–3]. Highly accurate precipitation nowcasting is vital in support of various operational activities, e.g., disaster relief, energy management, and flood early warning

**Citation:** Kong, D.; Zhi, X.; Ji, Y.; Yang, C.; Wang, Y.; Tian, Y.; Li, G.; Zeng, X. Precipitation Nowcasting Based on Deep Learning over Guizhou, China. *Atmosphere* **2023**, *14*, 807. https://doi.org/10.3390/ atmos14050807

Academic Editor: Stefano Federico

Received: 19 February 2023 Revised: 22 April 2023 Accepted: 26 April 2023 Published: 28 April 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

systems. Consequently, the accuracy of precipitation nowcasting has a critical socioeconomic impact.

Currently, the primary means of improving the accuracy of precipitation nowcasting involve the integration of extrapolation and ensemble numerical weather prediction (NWP) [4–8]. However, literature shows that global coarse-resolution NWP models have challenges in generating accurate precipitation with a lead time of 0–2 h [9,10], affected by the spin-up issue and the difficulties in non-Gaussian data assimilation. Although convection-permitting models with high resolution could greatly improve performance, these models are computationally expensive. Traditional integration of extrapolation methods, such as cross-correlation tracking and optical flow, has been used widely in operational weather forecasting [11–19]. Still, they exhibit a marked reduction in forecasting skills with increasing lead time owing to their linear operating limits and their deficiencies in forecasting storm growth and decay [20].

In recent years, the application of deep learning (DL) [21] has become popular in the field of meteorology [22–24] owing to its nonlinear operation, and it has achieved substantial advances in terms of quantitative precipitation nowcasting using radar echo extrapolation [25–32]. In comparison with the traditional integration of extrapolation and NWP, DL models can realize relatively accurate precipitation nowcasting. For example, Shi et al. [27] involved the convolution operation in input-to-state and state-to-state transitions in the transformation of a two-dimensional image into a three-dimensional tensor and proposed the Convolutional Long Short-Term Memory (ConvLSTM) model. This approach, which was shown capable of effectively capturing spatial correlations and further realizing the extrapolation of spatiotemporal sequences, was applied to achieve quantitative precipitation nowcasting. Subsequently, they involved the convolution operation in the recurrent gated unit (GRU) and used a subnetwork to output a location-variant connection structure before performing state transitions. Then, they proposed the trajectory GRU that could handle the spatial correlations better and perform more accurately than previous methods [28]. Wang et al. [29] designed the Spatiotemporal LSTM (ST-LSTM) model that can transfer memories vertically and horizontally, and then they proposed the Predictive Recurrent Neural Network (PredRNN) model. To strengthen the power for modeling short-term dynamics and to alleviate the vanishing gradient problem, they improved the PredRNN model to PredRNN++ [30], which incorporates a cascaded dual-memory structure and a gradient highway unit. Ji et al. [31] exploited the advantages of different DL model architectures in combination with the ConvLSTM unit in the U-Net generator and proposed the CLGAN model, which can better capture the precipitation object and its characteristics. Chen et al. [32] compared and analyzed the extrapolation prediction of radar echoes using the ConvGRU method, cross-correlation method, optical flow method, and particle filter method and found that the output of the ConvGRU method more closely matched the location, intensity, and shape of actual radar echoes. Many other studies have modified existing techniques to obtain relatively better nowcasting performance in comparison with that achieved using the integration of extrapolation [33–47]. However, precipitation calculated from radar echoes is based on the Z–R relationship. We believe that the Z–R relationship cannot describe the nonlinear relationship between a radar echo and precipitation. Moreover, it is difficult to avoid calculation errors derived from the Z–R relationship. Factors such as radar model, detection range, and clutter interference will lead to poor universality of a DL model. Consequently, the objective of this study was to construct DL models for precipitation nowcasting based on precipitation observation data, which could directly obtain precipitation nowcasting.

The reason we use DL models is their nonlinear operation, and the ConvLSTM model is one of the classical models used in the precipitation nowcasting problem. It can extract the spatial characteristics while capturing the time characteristics efficiently. PredRNN improved memory information by transferring memories vertically and horizontally, and it can perform better for forecast time compared with ConvLSTM. The Lucas–Kanade Optical Flow (LK Optical Flow) is a good method that improved optical flow vector calculation for traditional optical flow [48], and it was selected as a traditional baseline. We hope to compare the two typical DL models with LK Optical Flow and try to construct highly accurate models based on precipitation observation data.

In this study, we formed two high-quality datasets: a pre-processed precipitation dataset (Pre-P dataset) and a heavy precipitation dataset (Hea-P dataset). Then, we constructed precipitation DL models that could directly obtain precipitation nowcasting based on precipitation observation data. Furthermore, we improved the forecasting skill of the DL models by increasing the length of the time series of the training set and by training with the Hea-P dataset. The study area selected was Guizhou Province in China (29.2–24.7◦ N, 103.6–109.6◦ E), which features highly complex terrain and frequently experiences geological disasters and floods [49,50]. Moreover, the region is also lacking in terms of high-accuracy precipitation nowcasting.

The remainder of this manuscript is organized as follows. Section 2 describes the data and methods used in the study. Section 3 comprehensively evaluates both the DL models trained with different datasets and the LK Optical Flow method. Finally, the conclusions and a discussion are presented in Section 4.

#### **2. Data and Method**

### *2.1. Guizhou Automatic Weather Station (AWS) Observations*

High-quality real-time precipitation observation data (resolution: 1 min) recorded over a nine-year period from 20:00 (all times UTC+8) on 31 December 2012 to 20:00 on 31 December 2021 (observation day: from 20:00 on the previous day to 20:00 on the same day) by 1835 automatic weather stations (AWS) in Guizhou Province (as shown in Figure 1) were obtained from the Guizhou Meteorological Information Center.

**Figure 1.** Distribution of the 1835 automatic weather stations in Guizhou Province, China.

### *2.2. Pre-Processed Precipitation Dataset (Pre-P Dataset)*

The number of missing values in the observation data was reduced using the neighborhood method (where the time series before and after the missing values and the data of nearby stations were complete). Then, the bilinear interpolation method was used to interpolate station data within the region 24.5–29.5◦ N, 103.5–110◦ E to the gridded field (resolution: 0.125◦ × 0.125◦). Data normalization resulted in a high-quality dataset with a range of variation of [0, 1]. Ultimately, the pre-processed(Pre-P) dataset was formed, comprising 788,400 samples.

### *2.3. Heavy Precipitation Dataset (Hea-P Dataset)*

Data in the Pre-P dataset from May–September in each of the nine years from 2013–2021 were selected, and samples with little or no precipitation were removed. Consequently, the heavy precipitation (Hea-P) dataset was determined, comprising 332,640 samples.

### *2.4. Lucas–Kanade (LK) Optical Flow Method*

The optical flow method is defined as a method for calculating the intensity of image pixel points over time to infer the speed and direction of object movement. It finds the correspondence between the previous and the current frame based on changes in the pixels of the image in the time series and the correlation between adjacent frames, from which it calculates the motion information of objects between adjacent frames. This study employed the LK Optical Flow method [48], the process of which can be divided into two steps:

In the first step, the optical flow vector field is calculated. The instantaneous velocity two-dimensional vector field of the changing trend of grayscale at each point on the image and the optical flow vector field are solved using the following three basic assumptions.

(1) The brightness is constant, i.e., when the same target moves between different frames, the brightness will not change in any way. According to this assumption, the basic equation of the optical flow method can be obtained:

$$I(\mathbf{x}, y, t) = I(\mathbf{x} + d\mathbf{x}, y + dy, t + dt) \tag{1}$$

where *x* and *y* are the coordinate positions of a pixel in the image, *t* is the time series in which the image is located, and (*dx*, *dy*) is the distance moved to the next frame, using *dt* time, i.e., the light intensity of the pixel before and after motion has remained constant.

(2) The movement in continuous time is a 'small movement', which means that the change of time will not cause a drastic change in the target position. Expanding the Taylor section of Equation (1) yields:

$$I(\mathbf{x}, y, t) = I\left(\mathbf{x}, y, t + \frac{\partial I}{\partial \mathbf{x}} d\mathbf{x} + \frac{\partial I}{\partial y} dy + \frac{\partial I}{\partial t} dt + \varepsilon\right),\tag{2}$$

where *ε* represents a second-order infinitesimal term that is negligible. Assuming the velocity vectors *u* and *v* of the optical flow in the *x*-axis and *y*-axis directions, the formulas will be as follows:

$$Iu = \frac{d\mathbf{x}}{dt}, \ v = \frac{d\mathbf{y}}{dt}.\tag{3}$$

The partial derivatives of the gray scale of pixel points in the image on *x*, *y*, and *t* are as follows:

$$I\_{\chi} = \frac{\partial I}{\partial x'} \; I\_{\mathcal{Y}} = \frac{\partial I}{\partial y'} \; I\_{\mathcal{Y}} = \frac{\partial I}{\partial t'} \; \tag{4}$$

Then, substituting Equation (2) into Equation (1) and dividing by *dt* gives the following:

$$I\_{\mathbf{x}}\mathbf{u} + I\_{\mathbf{y}}\mathbf{v} + I\_{\mathbf{l}} = \mathbf{0},\tag{5}$$

where *Ix*, *Iy*, and *Iz* can be obtained from image data, and (*u*, *v*) is the optical flow vector.

(3) Spatial consistency (or spatial continuity) means that some points which around a specific point have the same optical flow field or velocity.

There are two unknowns in only one constrained equation. For analyzing the *n* × *n* region around a pixel, the pixel motion of the local region is assumed to be consistent. Then, the *n* × *n* equations can be established in the following matrix form:

$$
\begin{bmatrix} I\_{\ge 1} & I\_{y1} \\ I\_{\ge 2} & I\_{y2} \\ \vdots & \vdots \end{bmatrix} \quad \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} I\_{t1} \\ -I\_{t2} \\ \vdots \end{bmatrix} \tag{6}
$$

and written as vector:

$$A\stackrel{\rightarrow}{u} = \stackrel{\rightarrow}{b}\_{\prime} \tag{7}$$

where <sup>→</sup> *u* is the velocity vector:

$$
\stackrel{\rightarrow}{u} = (A^T A)^{-1} A^T \stackrel{\rightarrow}{b}\_{\prime} \tag{8}
$$

Then, the equations are solved using the least squares method.

⎡ ⎢ ⎣

The second step is to use the real-time precipitation field at the initial time and the optical flow vector field calculated in the first step to extrapolate the precipitation field. In this study, the semi-Lagrangian advection scheme was used to extrapolate the precipitation field [48]. The semi-Lagrangian advection formula can be expressed as follows:

$$F(t\_0 + \tau, \mathbf{x}) = F(t\_0, \mathbf{x} - \mathbf{a}),\tag{9}$$

where *F* is the extrapolated forecast precipitation field, *F* is the real-time precipitation observation at the initial time, *t*<sup>0</sup> is the lead time, *τ* is the forecast time, and *α* is the displacement vector within the lead time.

The semi-Lagrangian advection scheme divides the lead time *τ* into *N* steps for extrapolation, in this study, the lead time Δ*t* was 6 min. The displacement amount *α* can be obtained from the following iteration:

$$
\alpha^{(n+1)} \Delta t \stackrel{\rightarrow}{V} (t\_{0\prime} \ge -\frac{\alpha^{(n)}}{2}),
\tag{10}
$$

where <sup>→</sup> *V <sup>t</sup>*0, *<sup>x</sup>* <sup>−</sup> *<sup>α</sup>*(*n*) 2 is the velocity vector of the precipitation at point *x* − *<sup>α</sup>* <sup>2</sup> , and <sup>→</sup> *V* is the optical flow vector field, assuming the initial value of *α* is 0. The total displacement in the lead time is the sum of the displacements in *N* steps, and *n* is the number of iterations.

The specific extrapolation scheme uses the remapped particle-mesh Semi-Lagrangian scheme improved by Reich [51].

### *2.5. Convolutional Long Short-Term Memory (ConvLSTM) Model*

The ConvLSTM model can extract the two-dimensional spatial characteristics of precipitation while capturing the time characteristics. The convolution operation and pooling operation are integrated into the LSTM, which adds the "input gate," "output gate," and "forget gate" based on the RNN model and can record historical data over a longer time and improve forecast accuracy. Then, realizing the spatiotemporal sequence prediction of the precipitation, the main calculation equations in the model are as follows [27]:

$$\mathbf{i}\_{l} = \sigma(\mathcal{W}\_{\text{xi}} \ast \mathcal{X}\_{l} + \mathcal{W}\_{\text{hi}} \ast \mathcal{H}\_{t-1} + \mathcal{W}\_{\text{ci}} \odot \mathcal{C}\_{t-1} + \mathcal{b}\_{i}),\tag{11}$$

$$f\_t = \sigma \left( \mathcal{W}\_{xf} \* \mathcal{X}\_t + \mathcal{W}\_{hf} \* \mathcal{H}\_{t-1} + \mathcal{W}\_{cf} \odot \mathcal{C}\_{t-1} + b\_f \right),\tag{12}$$

$$\mathbb{C}\_{t} = f\_{t} \odot \mathbb{C}\_{t-1} + i\_{t} \odot \tanh(\mathcal{W}\_{\text{xc}} \ast X\_{t} + \mathcal{W}\_{\text{hc}} \ast H\_{t-1} + b\_{\text{c}}),\tag{13}$$

$$O\_l = \sigma(\mathcal{W}\_{\text{xo}} \* X\_l + \mathcal{W}\_{\text{ho}} \* H\_{t-1} + \mathcal{W}\_{\text{co}} \odot \mathcal{C}\_t + b\_o), \tag{14}$$

$$H\_t = O\_t \odot \tanh(\mathcal{C}\_t) \tag{15}$$

where *Xt* and *Ht* represent input data and output data, respectively; *it*, *ft*, and *Ot* represent the input gate, forget gate, and output gate, respectively; *Ct* is the cell state; "∗" denotes the convolution operator; and " " denotes the Hadamard product.

### *2.6. Predictive Recurrent Neural Network (PredRNN)*

For more accurate transmission of time information, the new ST-LSTM was proposed, which can transfer memories vertically and horizontally. Then, the PredRNN model was obtained based on the ConvLSTM model, and the complete formula of the ST-LSTM model can be expressed as follows [29]:

$$\mathbf{i}\_{l} = \sigma \left( \mathbf{W}\_{\text{xi}} \ast \mathbf{X}\_{l} + \mathbf{W}\_{\text{hi}} \ast \mathbf{H}\_{t-1}^{l} + \mathbf{b}\_{l} \right), \tag{16}$$

$$f\_t = \sigma \left( \mathcal{W}\_{\mathbf{x}f} \* \mathcal{X}\_t + \mathcal{W}\_{\mathbf{h}f} \* H\_{t-1}^l + b\_f \right), \tag{17}$$

$$\mathbf{C}\_{t}^{l} = f\_{t} \odot \mathbf{C}\_{t-1}^{l} + i\_{\mathrm{f}} \odot \tanh\left(\mathcal{W}\_{\mathrm{xg}} \ast \mathbf{X}\_{\mathrm{f}} + \mathcal{W}\_{\mathrm{hg}} \ast H\_{t-1}^{l} + b\_{\mathrm{g}}\right),\tag{18}$$

$$\dot{\mathbf{u}}'\_{t} = \sigma \left( \mathcal{W}'\_{xi} \* \mathcal{X}\_{t} + \mathcal{W}\_{mi} \* \mathcal{M}\_{t}^{l-1} + \mathcal{b}'\_{i} \right) \tag{19}$$

$$f\_t' = \sigma \left( \mathcal{W}\_{xf}' \* \mathcal{X}\_t + \mathcal{W}\_{mf} \* \mathcal{M}\_t^{l-1} + b\_f' \right), \tag{20}$$

$$M\_t^l = f\_i' \odot M\_t^{l-1} + \mathbf{i}\_t' \odot \tanh\left(\mathcal{W}\_{\text{xg}}' \ast \mathbf{X}\_t + \mathcal{W}\_{\text{mg}} \ast M\_t^{l-1} + b\_{\mathcal{S}}'\right),\tag{21}$$

$$O\_l = \sigma \left( \mathcal{W}\_{\text{xo}} \ast \mathcal{X}\_l + \mathcal{W}\_{\text{ho}} \ast \mathcal{H}\_{t-1}^l + \mathcal{W}\_{\text{co}} \ast \mathcal{C}\_t^l + \mathcal{W}\_{\text{mo}} \ast \mathcal{M}\_t^l + b\_o \right), \tag{22}$$

$$H\_t^l = O\_l \odot \tanh\left(W\_{1 \times 1} \* \left[\mathbf{C}\_{t\prime}^l M\_t^l\right]\right),\tag{23}$$

*Wxi*, *Whi*, *Wmi*, *Wx f* , *Wh f* , *Wm f* , *Wxg*, *Whg*, *Who*, *Wmo*, *Wxi*, *Wx f* , *Wxg* are all weight parameters, and *bi*, *bf* , *bg*, *bo*, *b i* , *bf* , *bg* are learnable offset parameters. With *C<sup>l</sup> <sup>t</sup>* representing time memory, *M<sup>l</sup> <sup>t</sup>* representing spatial memory, *H<sup>l</sup> <sup>t</sup>* representing the value of the hidden layer, subscript *t* representing the time step, and superscript representing the kth hidden layer existing in the ST-LSTM network.

For the parameters of DL models in this paper, we set the initial learning rate and the epochs as 0.001 and 20, respectively. The root mean square error (RMSE) was selected as the loss function.

### *2.7. Verification Metrics*

2.7.1. Root Mean Square Error (RMSE)

The expression for the RMSE is as follows:

$$RMSE = \frac{1}{n} \sum\_{i=1}^{n} [(f\_i - O\_i)^2]^{\frac{1}{2}} \,\,\,\,\tag{24}$$

where *n* represents the number of grid samples of the space field, *fi* represents the forecast value of the *i*–th sample, and *Oi* represents the observation of the *i*–th sample. The smaller the *RMSE* value, the smaller the difference between the forecast value and the observed value, i.e., the smaller the forecast error.

2.7.2. Probability of Detection (POD), False Alarm Ratio (FAR), Probability of False Detection (POFD), and Equitable Threat Score (ETS)

Expressions for the probability of detection (POD), false alarm ratio (FAR), probability of false detection (POFD), and equitable threat score (ETS) are as follows [52–54]:

$$\text{POD} = \frac{N\_A}{N\_A + N\_\mathbb{C}} \,\text{}\tag{25}$$

FAR <sup>=</sup> *NB NA* + *NB* , (26)

$$\text{POFD} = \frac{N\_{\text{C}}}{N\_{B} + N\_{E}} \,\text{\,\,\,}\tag{27}$$

$$E = \frac{(N\_A + N\_B)(N\_A + N\_C)}{N\_A + N\_B + N\_C - E} \,\text{}\,\tag{28}$$

$$\text{ETS} = \frac{N\_A - E}{N\_A + N\_B + N\_C - E} \,\text{<}\tag{29}$$

where *NA* represents the number of forecast events that correspond to observed events (forecast has, observation has); *NB* represents the number of events that are null (forecast has, observation has not); *NC* represents the number of events that are missed (forecast has not, observation has); *NE* represents the number of events that both forecast, and observation do not occur. POD refers to the proportion of the predicted actual precipitation area in the total actual precipitation area; the larger the value, the higher the forecast accuracy. FAR refers to the proportion of the area with no actual precipitation in the forecast precipitation area in the total forecast precipitation area; the smaller the value, the smaller the forecast null rate. POFD refers to the proportion of the area that is missed in the actual precipitation area in the area where all actual precipitation does not occur; the smaller the value, the lower the FAR of the forecast. Additionally, the higher the ETS score, the better the forecast performance.

In this study, we applied five thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm) to calculate the ETS of the hourly graded precipitation. The forecasting skill of 6-min precipitation was calculated based on the threshold of 0.01 mm.

### 2.7.3. Method for Object-Based Diagnostic Evaluation (MODE)

The Method for Object-Based Diagnostic Evaluation (MODE) is a spatial evaluation method based on object attribute characteristics, which mainly evaluates a prediction by comparing and analyzing the attributes and similarities of the main observed and forecast areas of precipitation [55,56].

Through the operation of spatial convolution and then under different precipitation thresholds (i.e., 0.01, 0.1 and 2.5 mm), the important areas of precipitation to be studied are calculated and identified as follows [57]:

$$\mathcal{C}(\mathbf{x}, y) = \sum\_{\boldsymbol{\mu}, \boldsymbol{\upsilon}} \varphi(\boldsymbol{\mu}, \boldsymbol{\upsilon}) f(\mathbf{x} - \boldsymbol{\mu}, y - \boldsymbol{\upsilon}), \tag{30}$$

$$\varphi(u,v) = \begin{cases} \frac{1}{\pi R^2} & u^2 + v^2 \le R^2\\ 0, & u^2 + v^2 > R^2 \end{cases} \tag{31}$$

where *f* represents the original data field, *C* represents the convolutional field, *ϕ* represents the filter function, and (*x*,*y*) and (*u*,*v*) represent grid coordinates. The mask field *M* is obtained by threshold control of the convolutional field, i.e., the precipitation area in the convolutional field where precipitation intensity is greater than or equal to threshold *T* is calculated [55,56]:

$$M(\mathbf{x}, y) = \begin{cases} 1, & \mathcal{C}(\mathbf{x}, y) \le T \\ 0, & \mathcal{C}(\mathbf{x}, y) > T \end{cases} \tag{32}$$

By assigning grid points in the continuous region of *M* = 1 to the value of the corresponding grid points in the original precipitation field, the reconstruction field *F* can be obtained, which not only retains most of the original precipitation information of each object (precipitation without convolution processing), but also identifies the main area of falling precipitation when the precipitation threshold is reached. The formula for calculating *F* is as follows [55,56]:

$$F(\mathbf{x}, y) = M(\mathbf{x}, y) f(\mathbf{x}, y), \tag{33}$$

Then, according to Equation (33), certain important properties of the observation field and the precipitation field in these aeras of falling precipitation are calculated. In this study, area, angle, aspect ratio, centroid position of longitude, and centroid position of latitude were selected as the attributes for assessment and analysis of the area of falling precipitation.

According to Davis et al. [55,56], the matching rule proposed by uses the calculated attributes to match objects, i.e., all matches in the process come from the important areas of falling precipitation in the forecast field and in the observation field identified as needing study, and the calculation formula used in this process is as follows [55–58]:

$$D < \frac{Area\_O^{\frac{1}{2}} + Area\_f^{\frac{1}{2}}}{2},\tag{34}$$

where *Areao* and *Areaf* are the areas of the main areas of falling precipitation identified in the observation and forecast fields, and *D* is the centroid distance between them.

Finally, according to the weight and the confidence factor of the attribute, the total similarity between the observation field and the important area of falling precipitation in the precipitation field is calculated using the fuzzy logic method [57]:

$$I = \frac{\sum\_{i=1}^{n} \omega\_{i} c\_{i} G\_{i}}{\sum\_{i=1}^{n} \omega\_{i} c\_{i}} , \tag{35}$$

where *ci* and *ω<sup>i</sup>* represent the confidence factor and the weight for property *i*, respectively, which depend only on the specific properties of the subject of tax reduction, and where the confidence factor varies with the area size and distance of the area of falling precipitation; and *Gi* is the similarity factor of the *i*-th attribute that is a monotonous recursive function with a value range of between 1 and 0.

### **3. Results**

### *3.1. Data Quality Control (DQC) Evaluation*

We selected May–September 2019 as the verification period because several extreme precipitation events occurred during this period. For example, intense rainstorms that occurred during 5–11 June 15–19 June, 20–25 June, and 27 June to 1 July caused notable geological and flood disasters [59].

Reducing the number of missing values in the original AWS data produced a complete time series and improved the spatial resolution of the dataset, thereby making the training samples more representative of the observed precipitation. Data normalization narrowed the range of extreme precipitation and improved the learning efficiency of the DL models with respect to the features of precipitation observations during the training process. Both processes helped DL models that exhibited improved performance in terms of ETS. The Hea-P dataset before and after data quality control (DQC) was divided into a training set and a validation set with a ratio of 8:1, and ETS evaluations were performed for the verification period.

The ETSs of the DL models for 6-min precipitation nowcasting before and after DQC are shown in Figure 2. The average ETS of the ConvLSTM model and of the PredRNN model increased by 0.1458 and 0.0660, respectively. The ConvLSTM model showed greater improvement in comparison with the PredRNN model; however, the PredRNN model also exhibited reasonable improvement beyond the 12-min lead time.

**Figure 2.** Comparison of the ETSs of the DL models before and after data quality control (DQC).

Currently, there is no unified evaluation metric for minute-level precipitation. Consequently, we calculated the ETS of hourly graded precipitation accumulated from 6-min precipitation for different thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm), as listed in Table 1. The forecasting skill of the ConvLSTM model with thresholds of 0.1 and 1 mm was well improved, and the ETS for the threshold of 0.1 mm for the 1-h lead time increased by 0.2286. Overall, the PredRNN model exhibited greater improvement in comparison with that of the ConvLSTM model, i.e., the average ETS for thresholds exceeding 2.5 mm for the 1-h lead time and for thresholds below 2.5 mm for the 2-h lead time increased by more than 0.09 and by 0.1305, respectively.


**Table 1.** ETSs of hourly graded precipitation for DL models before and after DQC for different thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm).

### *3.2. DL Models Trained Using Datasets with Different Time Series Lengths*

Pre-P datasets with DQC and different time series lengths (i.e., 1, 3, 5, and 8 years) were divided into training sets, testing sets, and validation sets with ratios of 2:1:1, 6:1:1, 10:2:1, and 16:3:1, respectively, and the verification period for all evaluations was May– September 2019. Then, DL models were constructed using the training sets with different time series lengths (i.e., 1, 3, 5, and 8 years). Verifications including RMSE, ETS, POD, FAR, and POFD were applied to assess the 6-min precipitation nowcasting for the threshold of 0.01 mm generated using the DL models, as shown in Figure 3. With increasing length of the time series of the training set, the DL models showed a marked increase in terms of POD and ETS, and the stability of the nowcasts with time was enhanced, especially for the ConvLSTM model trained with the training set with the 8-year time series length, i.e., the average POD increased by 0.1066. Additionally, the PredRNN model trained with the training set with the 8-year time series length performed better than the ConvLSTM model, i.e., the RMSE decreased sharply, the FAR was low, and the ETS increased by 0.0763.

**Figure 3.** Forecasting skill of DL models trained with different time series lengths: (**a**) RMSE, (**b**) ETS, (**c**) POD, (**d**) FAR, and (**e**) POFD.

Taking the same verification period as above, we calculated the ETS of hourly graded precipitation (Table 2). The ETSs of the DL models were improved with the increasing length of the time series of the training set. For the training set with the 8-year time series length, the ETS of the ConvLSTM model for the threshold of 0.1 mm increased sharply by 0.2484 for the 2-h lead time, while that for thresholds below 2.5 mm for the 1-h lead time increased by more than 0.05. However, the PredRNN model showed better improvement in terms of ETS. For example, the average ETS for the 0–2-h lead time for different thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm) increased by 0.0852, while that for the 1-h lead time increased by 0.1067, and that for the threshold of 2.5 mm increased by 0.1464.

**Table 2.** Hourly graded precipitation ETSs for the two DL models trained with Pre-P datasets with different time series lengths for different thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm).


To summarize, increasing the length of the training set time series and the number of the training samples both improved the forecasting skill of the DL models. Overall, the PredRNN model exhibited greater improvement in comparison with the ConvLSTM model, especially for the 1-h lead time, and the PredRNN model trained with the Pre-P dataset with the 8-year time series length obtained the best performance in terms of ETS.

### *3.3. Deep Learning Model Training Using the Heavy Precipitation (Hea-P) Dataset*

Although increasing the length of the time series of the training set could improve the forecasting skill of the DL models, the model with the best performance still exhibited no improvement in terms of the ETS for hourly graded precipitation when compared with the results obtained using the LK Optical Flow method (Table 3). Even though the ETS of the PredRNN model with the best performance was slightly higher than that of the traditional LK Optical Flow method for the thresholds of 0.1, 2.5, and 8.0 mm for the 1-h lead time, the ETS for the threshold of 15.0 mm for the 1-h lead time was lower than that derived using the LK Optical Flow method. Moreover, the ETSs of the PredRNN model with different thresholds (i.e., 0.1, 2.5, 8, and 15 mm) decreased rapidly with time, resulting in a lower overall ETS in comparison with that derived using the LK Optical Flow method.

**Table 3.** ETSs of hourly graded precipitation for the LK Optical Flow method and the DL models trained using the two types of datasets for different thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm).


To further improve the forecasting skill of the DL models, we removed samples with little or no precipitation and left the extreme precipitation to form the Hea-P dataset, which we hoped would represent a more effective training sample. Taking the same verification period as described in Section 3.2, Hea-P datasets with different time series lengths (i.e., 1, 3, 5, and 8 years) were divided into training sets, and validation sets with a ratio of 1:1, 3:1, 5:1, and 8:1, respectively. We constructed DL models using the training sets with time series of different lengths. Then, evaluation of the 6-min precipitation nowcasting was performed based on the RMSE, ETS, POD, FAR, and POFD, which allowed comparative assessment of the LK Optical Flow method and the DL models trained with the two types of datasets (i.e., Pre-P dataset and Hea-P dataset), as shown in Figure 4.

**Figure 4.** Forecasting skill of the LK Optical Flow model and the PredRNN model trained with two types of datasets: (**a**) RMSE, (**b**) ETS, (**c**) POD, (**d**) FAR, and (**e**) POFD.

With a sharp increase in POD and reduction in RMSE for the 8-year time series length, the DL models trained with the Hea-P dataset performed better than those trained with the Pre-P dataset in terms of 6-min precipitation nowcasting; specifically, the PredRNN model outperformed the LK Optical Flow method in terms of forecasting skill. For the 1- and 3-year time series lengths, the DL models trained with the Hea-P dataset exhibited reduced ETSs and sharply increased FAR scores. Meanwhile, the DL models trained with the Hea-P datasets with the 5- and 8-year time series lengths also sharply increased in FAR scores. Overall, use of the Hea-P datasets improved the POD and resulted in sharply increased FAR values for the DL models because the DL models tended to generate more high-value precipitation nowcasting owing to training samples with high-value precipitation in the Hea-P dataset.

The ETSs of hourly graded precipitation for different thresholds are listed in Table 4. The DL models trained with the Hea-P dataset improved markedly. The PredRNN model trained with the training set with the 8-year time series length obtained notable improvement for thresholds exceeding 2.5 mm, i.e., the ETS of the PredRNN model for thresholds of 2.5, 8, and 15 mm for a 1-h lead time increased sharply by 0.1141, 0.0702, and 0.0648, respectively. The average ETS for thresholds below 2.5 mm for the 2-h lead time also increased by more than 0.05. Importantly, the ETS for the threshold of 8.0 mm increased by an average of 0.0676 for a 0–2-h lead time and by 0.0925 for a 1-h lead time. Additionally, the ETS of the ConvLSTM model for the threshold of 2.5 mm increased by 0.0801, while that for below the threshold of 1.0 mm decreased (note: the ETSs of the ConvLSTM model are not listed in the table).

**Table 4.** ETSs of hourly graded precipitation for the PredRNN model trained with two types of datasets for different thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm).


Comprehensive evaluation revealed the DL models with the performance in terms of ETS. The PredRNN model trained with the Hea-P dataset with the 8-year time series length obtained the best performance of all models in terms of ETS for hourly graded precipitation (Table 5). The PredRNN model outperformed the LK Optical Flow method for all thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm), especially thresholds exceeding 2.5 mm. Although the ConvLSTM model was improved by training with the Hea-P dataset, it still performed a lower ETS than that of the other two models for all thresholds.

**Table 5.** ETSs of hourly graded precipitation for the LK Optical Flow method and the DL models trained with the Hea-P dataset for different thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm).


The spatial distributions of ETSs for hourly graded precipitation derived using the LK Optical Flow method and the DL models trained with the Hea-P dataset with the 8-year time series length for the 1-h lead time are shown in Figure 5. For thresholds of 0.1, 1, 2.5, and 8 mm, the higher the threshold, the lower the ETS of hourly graded precipitation; however, the ETS for the threshold of 8 mm exhibited a sharp decrease. Areas with high ETSs for thresholds of 0.1 and 1 mm mainly occurred in eastern, southern, and southwestern parts of Guizhou, whereas high ETSs for thresholds of 2.5, 8, and 15 mm mainly occurred in north-central and southern parts of Guizhou. For thresholds of 8 and 15 mm, the ETS of the ConvLSTM model was close to zero owing to many false predictions in areas of high-value precipitation and high values of POFD and FAR. This can result from cumulative error magnified by iterative calculation. For the threshold of 15 mm, the ETSs of the LK Optical Flow method and the PredRNN model showed improvement mainly in the southern and north-central parts of Guizhou. Additionally, the ETSs of the PredRNN

model with different thresholds (i.e., 0.1, 1, 2.5, 8 and 15 mm) indicate that the PredRNN model outperformed both the ConvLSTM model and the LK Optical Flow method.

**Figure 5.** Spatial distributions of ETSs for hourly graded precipitation for the 1-h lead time, derived from the LK Optical Flow method and the DL models trained with the Hea-P dataset with the 8-year time series length.

Figure 6 presents the spatial distributions of ETSs for the 2-h lead time. The ETSs of both the LK Optical Flow method and the PredRNN model for different thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm) decreased sharply for the 2-h lead time with the same pattern of spatial distribution as that shown for the 1-h lead time. The ETS of the PredRNN model for thresholds of 0.1 and 15 mm maintained values of 0.3–0.5 and 0.2–0.3, respectively, which were superior to the values of the other two models.

**Figure 6.** Spatial distributions of ETSs for hourly graded precipitation for the 2-h lead time, derived from the LK Optical Flow method and the DL models trained with the Hea-P dataset with the 8-year time series length.

In summary, the DL models trained with the Hea-P dataset improved notably in terms of ETS for 6-min and hourly graded precipitation, especially for thresholds of 2.5, 8, and 15 mm. We determined that the PredRNN model trained with the Hea-P dataset with an 8-year time series length performed best, outperforming the LK Optical Flow method for all thresholds. Furthermore, it performed well for thresholds of 2.5, 8, and 15 mm in the north-central and southern parts of Guizhou.

### *3.4. Structure Evaluation on a Rainstorm Case*

As a case for evaluation, we selected a severe rainstorm event that occurred on 6 June 2019, which produced a two-hour period (01:30–03:24) of intense precipitation [59]. We used the MODE approach to further assess the performance of both the LK Optical Flow method and the PredRNN model trained with the Hea-P dataset with the 8-year time series length in relation to this rainstorm case. The MODE consisted of an evaluation of the

6-min precipitation for the threshold of 0.01 mm and an evaluation of the hourly graded precipitation for thresholds of 0.1, 1, and 2.5 mm.

The attribute values of the 6-min precipitation objects for the observations and nowcasts, calculated for a threshold of 0.01 mm, included the area, axis angle, aspect ratio, zonal centroid, and meridional centroid, from which we obtained the total similarity of the precipitation objects. All attribute values and the total similarity of the 6-min precipitation objects for both the LK Optical Flow method and the PredRNN model are shown in Figure 7. The attribute values of the observations have a wide range and show rapid fluctuation during the rainstorm, resulting in both the shape and the position of the area of precipitation varying wildly and rapidly with time. Meanwhile, the zonal centroid and the meridional centroid indicate that the precipitation system first evolved toward the southwest and then moved rapidly toward the northeast. With the characteristic of linear variation, the attribute values of the LK Optical Flow method change smoothly with time, whereas the attribute values of the PredRNN model fluctuate widely and are much closer to those of the observations.

**Figure 7.** Attribute values of 6-min precipitation objects for the LK Optical Flow method and the PredRNN model trained with the Hea-P dataset with the 8-year time series length: (**a**) area, (**b**) axis angle, (**c**) aspect ratio, (**d**) zonal centroid, (**e**) meridional centroid, and (**f**) total similarity.

It is evident from Figure 8 that the precipitation nowcasting of the LK Optical Flow method for 6–12-min lead times is highly similar to the observations and shows only minimal linear movement in position and shape. Nevertheless, with such characteristics, the deviation between the observations and the nowcasting of the LK Optical Flow method increases from 18 to 30 min. Conversely, the nowcasting generated using the PredRNN model for 18 to 30 min is closer to the observations because of the ability of the model to capture the nonlinear changes in the observations, demonstrating especially good performance for thresholds of 8 and 15 mm.

We calculated attribute values of hourly graded precipitation for thresholds of 0.1, 1, and 2.5 mm (Table 6). The deviations between the observations and the nowcasting generated using the LK Optical Flow method are more pronounced for a 2-h lead time owing to the axial angle, aspect ratio, zonal centroid, and meridional centroid. Meanwhile, with an increase in the precipitation threshold, the PredRNN model trained with the Hea-P dataset with the 8-year time series length gradually increased in total similarity, substantially outperforming the LK Optical Flow method. For the threshold of 2.5 mm, the

PredRNN model maintains a value of total similarity of >0.9, whereas the value for the LK Optical Flow method rapidly drops below 0.75.

**Figure 8.** Distributions of observed and predicted precipitation for 6–30-min lead times for the evaluated rainstorm case.

**Table 6.** Hourly graded precipitation MODE for different thresholds (i.e., 0.1, 1, and 2.5 mm) for the LK Optical Flow method and the PredRNN model trained with the Hea-P dataset with the 8-year time series length.


The results of the above quantitative evaluation are well reflected in the spatial distributions of hourly precipitation shown in Figure 9. It is evident that the PredRNN model performs better than the LK Optical Flow method for thresholds exceeding 2.5 mm, i.e., it can well capture the nonlinear movement and evolution of the precipitation for a 2-h lead time. The LK Optical Flow method shows poor performance owing to the minimal changes in the shape and position of the area of precipitation with time, resulting in false predictions exceeding the threshold of 8.0 mm in northern Guizhou.

**Figure 9.** Distributions of observed and predicted hourly precipitation for 6–30 min lead times for the evaluated rainstorm case.

The characteristics of tiny movement and linear evolution mean that the precipitation nowcasting of the LK Optical Flow method is unable to capture nonlinear evolution that does not follow the rule of the latest optical flow vector in the observations. The errors also increased over time owing to the inability of the LK Optical Flow method to adapt to the rapid and extreme evolution of precipitation. However, the PredRNN model trained with the Hea-P dataset can overcome the deficiencies of the traditional LK Optical Flow method, capturing complex nonlinear evolution and establishing a reasonable mapping network during the drastic evolution of precipitation. Thus, the nowcasting of the PredRNN model is closer to the observations, especially for thresholds of 2.5, 8, and 15 mm.

#### **4. Conclusions and Discussion**

In this study, we formed two high-quality datasets (the Pre-P dataset and the Hea-P dataset) based on AWS precipitation observation data. Precipitation nowcasting at 6-min intervals for the lead times of 0–2 h was generated for Guizhou using the traditional LK Optical Flow method and the ConvLSTM and PredRNN DL models trained with the Pre-P and Hea-P datasets. Evaluations based on the RMSE, ETS, POD, and FAR were used to assess the performance of the different models in generating precipitation nowcasting. A rainstorm case was evaluated using the MODE approach to further examine the performance of the LK Optical Flow method and the PredRNN model trained with the Hea-P dataset with the 8-year time series length. The results obtained and the conclusions derived were as follows.

Reducing the number of missing values in the AWS observation data improved the quality of the DL training sample data. Data normalization also improved the training efficiency of the DL models. Both processes helped improve the forecasting skill of the DL models. The greater the length of the time series of the training dataset, the better the forecasting skill of the DL model. Increasing the time series length and the number of samples in the training data improved the POD and ETS of the DL models and enhanced the stability of the nowcasting over time. The PredRNN model was most improved for hourly graded precipitation, especially for the 1-h lead time.

Training with the Hea-P dataset further improved the forecasting skill of the DL models and sharply increased the ETSs for thresholds of 2.5, 8, and 15 mm, especially for the 1-h lead time. The improvement in the PredRNN model was greater than that in the ConvLSTM model. The PredRNN model trained with the Hea-P dataset with the 8-year time series length outperformed the traditional LK Optical Flow method for all thresholds (i.e., 0.1, 1, 2.5, 8, and 15 mm) and obtained the best performance in terms of ETS in comparison with the other DL models examined in this study. It also performed high ETSs for thresholds exceeding 2.5 mm in the north-central and southern parts of Guizhou. Additionally, the DL models tended to predict high-value precipitation owing to the high-value precipitation training samples in the Hea-P dataset, which is why most DL models showed a sharp increase in FAR values. The ETS of the ConvLSTM model for the threshold of 2.5 mm increased, but the ETS for the thresholds of 0.1 and 1 mm decreased owing to a sharp increase in the FAR. Owing to the magnified cumulative error through iterative calculation, the ConLSTM model performed many false predictions in high-value precipitation, and the ETS of the ConvLSTM model for thresholds of 8 and 15 mm was close to zero. Unlike radar data, observation instruments and measurement standards for precipitation observation data are unified in the industry. Generally, the spatial resolution of AWS data is higher than that of weather radar. Consequently, the DL models considered in this study could have reasonable transferability to other regions. Thus, the same approach could be used to construct highly accurate precipitation nowcasting DL models based on high-quality observation data.

The rainstorm case considered for evaluation revealed the characteristics of minimal movement and linear evolution in the traditional LK Optical Flow method. The generated precipitation nowcasting exhibited increasing error over time because the latest optical flow vector could not capture subsequent nonlinear evolution. Conversely, the PredRNN model trained with the Hea-P dataset could overcome the deficiencies of the traditional LK Optical Flow method and could capture the complex nonlinear evolution. Thus, the generated precipitation nowcasting was much closer to the observations. Specifically, the PredRNN model outperformed the traditional LK Optical Flow method for the threshold of 2.5 mm in the evaluated rainstorm case.

Currently, DL models using multisource observation data can overcome certain physical constraints and generate radar echoes that are physically more reasonable and of reference significance. For example, Li et al. [26] used a DL model and multisource data to produce radar echoes with physical characteristics that were improved in comparison with those derived using single-source observation data. In future work, DL models with multisource observations that include parameters such as temperature, pressure, and wind speed should be investigated. The inclusion of such dynamic and thermodynamic meteorological information will further improve the forecasting skill of DL models. Additionally, a comparative discussion with radar echo and more reasonable evaluation metrics will be added, such as radially averaged power spectral density [26,60,61].

**Author Contributions:** D.K. and X.Z. (Xiefei Zhi). contributed to the conception and design of the study. C.Y. and Y.T. contributed to the analysis. X.Z. (Xiaotuan Zeng), G.L., Y.J., D.K. and Y.W. organized the database. All authors contributed to writing and revising the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The study was jointly supported by the Natural Science Foundation of Hebei Province, China (Grant No.: D2021304003) and The Meteorological Science, Technology Open Research Fund of Guizhou Meteorological Bureau (Grant No.: KF201704), the Provincial and Municipal Joint Fund

Project of Guizhou Province Meteorological Bureau "Research on precipitation forecasting of short-term in Southwest Guizhou", the Provincial and Municipal Joint Fund Project of Guizhou Province Meteorological Bureau "Research on AI-based intelligent weather monitoring and dispatching technology".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The observation data used in this study were provided by the Guizhou Meteorological Information Center and are not public in nature. The dataset was newly established in this study.

**Acknowledgments:** The authors thank the Guizhou Meteorological Information Center for its cooperation and data support.

**Conflicts of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as potential conflict of interest.

### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
