A Nonlinear Grid Transformation Method for Extrapolating and Predicting the Convective Echo of Weather Radar

Abstract

A nonlinear grid transformation (NGT) method is proposed for weather radar convective echo extrapolation prediction. The change in continuous echo images is regarded as a nonlinear transformation process of the grid. This process can be reproduced by defining and solving a 2 × 6 transformation matrix, and this approach can be applied to image prediction. In ideal experiments with numerical and path changes of the target, NGT produces a prediction result closer to the target than does a conventional optical flow (OF) method. In the presence of convection lines in real cases, NGT is superior to OF: the critical success index (CSI) for 40 dBZ of the echo prediction at 60 min is approximately 0.2 higher. This is due to the better estimation of the movement of the whole cloud system in the NGT results since it reflects the continuous change in the historical images. For the case with a mesoscale convective complex, the NGT results are better than the OF results, and a deep learning result is cited from a previous study for the same case for 20 and 30 dBZ. However, the result is the opposite for 40 dBZ, where the deep learning method may produce an overestimation of the stronger echo.

1. Introduction

A weather radar is an important means to monitor and warn of severe convective weather [1,2,3]. Although data assimilation of weather radar helps to improve numerical forecasting [4,5], extrapolation prediction based on weather radar echoes remains the most direct method for nowcasting a convective event, especially for mesoscale convective weather [6,7] with rapid development and dramatic changes. By means of the current observed radar data, the movement and variation of the echo are obtained, and the position and intensity of the echo in the next few minutes or longer period are estimated.

The methods of extracting radar echo moving features from horizontal two-dimensional radar images can be divided into two categories according to the target. One focuses on identifying and tracking a whole cloud, such as the Thunderstorm Identification, Tracking, Analysis, and Nowcasting algorithm (TITAN) [8]. The cloud is determined by the maximum circumscribed ellipse under a given reflectivity threshold to obtain its trajectory and variation characteristics. Another category of methods focuses on obtaining the vector field of echo movement under a uniform grid. The Tracking Radar Echo by Correlation algorithm (TREC) [9], a common algorithm in the past, determines the echo moving vector on the specified grid point based on the cross-correlation between the subunits around adjacent time points. Both TITAN and TREC contribute to a classical auto-nowcast system [3]. The Optical Flow (OF) method can also produce a vector field based on images around adjacent time points. By assuming the target gray level is constant and introducing the global smoothing condition (known as Horn-Schunk OF, [10]) or the local consistency condition (known as Lucas-Kanade OF, [11]), moving vectors can be solved based on the local change and spatial gradient of the images. A recent study shows that different nowcasting schemes based on OF can both produce better results than TREC to some extent [12]. In addition, an improved prediction method combines OF with an assimilation system and considers water content and principal component analysis [13].

However, some challenges remain in the feature extraction and extrapolation of convective clouds in the above methods. The alternating development of convective cells often occurs in convective systems. The moving direction of the whole cloud system typically depends on the combination of individual cell movement and the direction of new cell emergence, which may not be correctly reproduced by a method based on a few images at adjacent times, such as TREC and OF. Although TITAN can provide the echo movement trajectory based on longer historical echo image samples, each convective cell in a cloud system may not only move in different directions but also merge, resulting in complex topological relationships that are difficult to extrapolate. In addition, the shape and intensity of convective echoes can vary locally and rapidly, which makes it difficult to capture classical local image features such as Scale-Invariant Feature Transform (SIFT) [14] and Harris corner detectors [15]. Therefore, other image recognition and analysis algorithms based on those features cannot be effectively applied to the extrapolation of weather radar echoes. In general, there is a lack of algorithms that consider multiple historical times and do not rely on identifying individual targets.

There are also some learning methods that have been applied to radar echo and precipitation nowcasting [16,17,18,19,20,21,22,23] in recent years. These methods allow multiple observed image inputs and predict the image in the future, which appears to address the limitations mentioned above. However, since the model parameters and modeling process are implicit and there is a requirement for a large number of samples for training, the reliability of these methods requires further verification in application.

In this paper, a new method is proposed for the extrapolation and prediction of weather radar echo images. The change in continuous echo images is regarded as a nonlinear transformation process of the grid by a transformation matrix. By defining and solving this matrix, the extrapolation and prediction of the echo image can be realized. The process of establishing the new method is described in Section 2. The echo prediction experiments are shown in Section 3, including a set of ideal cases and three real cases. A summary and discussion are provided in Section 4.

2. Nonlinear Grid Transformation Method

2.1. Matrix for Grid Transformation

The concept of the transformation matrix is first introduced in this section, including the conventional transformation matrix and the extended transformation matrix used in this paper. The connection between the transformation matrix and the echo image prediction is introduced in Section 2.2 and Section 2.3.

The general form of the transformation matrix used in this paper is a matrix M_p×p, as shown in Equation (1):

$$M_{p\times p}=\left[\begin{array}{c}{m}_{1,1} m_{1,2} \dots m_{1,p-1} m_{1,p}\\ m_{2,1} m_{2,2} \dots m_{2,p-1} m_{2,p}\\ \dots \\ m_{p-1,1} m_{p-1,2}\dots m_{p-1,p-1} m_{p-1,p}\\ 0 0 \dots 0 1\end{array}\right]$$

(1)

where the first p − 1 rows are independent elements, the value in row p and column p is 1, and the other values in row p are 0. When p is 3, a two-dimensional point set [X₁, Y₁] can be transformed to [X₂, Y₂] as follows:

$${\left[X_2,Y_2,1\right]}^T=M_{3\times 3}{\left[X_1,Y_1,1\right]}^T$$

(2)

where the “1” at the third location in the square bracket provides a nonhomogeneous term when performing linear transformation of the point coordinates. Equation (2) is known as “affine transformation” [24], which is a common algorithm of image geometric transformation. With this matrix, the original point set can be translated, rotated, scaled, and deformed (sheared, stretched, etc.). For the convenience of observation, [X₁, Y₁] is set to an orthogonal grid, and Figure 1a shows a comprehensive example, which includes changes such as moving to the right, magnification, counterclockwise rotation, and stretching the upper right corner. This transformation is repeatable. After N-1 times left multiplication by M_3×3, the transformation result [X_N, Y_N] is obtained in the same way:

$${\left[X_{\mathrm{N}},Y_{\mathrm{N}},1\right]}^T={\left(M_{3\times 3}\right)}^{\mathrm{N}-1}{\left[X_1,Y_1,1\right]}^T$$

(3)

Related principles and examples of such common 2D image transformation are seen in [24]. Note that only the first two rows of M_3×3 can actually be taken if only one transformation is needed. In addition, the inverse matrix of M_3×3 represents the opposite transformation. Usually, when at least three groups of X and Y before and after the transformation are known, the six parameters in M_3×3 can be solved to achieve a similar transformation for other point sets, or the expected transformation can be achieved by directly specifying the 6 parameters.

However, this conventional affine transformation has a limited ability to characterize image deformation. Note that in Figure 1a, although the shape has changed to some extent, the opposite edges are still parallel, so it is obviously not suitable for characterizing the changes in weather echo images. The transformation could be expanded by adding more parameters. It is easy to extend M in the form of Equation (1). When expanding the formula with coordinates multiplied by M, the newly added variables must be known or knowable and cannot be completely linearly related to the point sequence of the original X₁ and Y₁. A nonlinear XY term is proposed to be added, and the transformation matrix M is 4 × 4 as follows:

$${\left[X_2,Y_2,X_2{Y}_2,1\right]}^T=M_{4\times 4}{\left[X_1,Y_1,X_1{Y}_1,1\right]}^T$$

(4)

An example of the transformation of Equation (4) is given in Figure 1b, which shows that it can already characterize the case where the contrast is not parallel. However, every edge of the orthogonal grid after transformation is still a straight line. To further improve the deformation ability of the transformation, the nonlinear terms X² and Y² are added, and a 6 × 6 transformation matrix is used:

$${\left[X_2,Y_2,X_2{Y}_2,X_2^2,Y_2^2,1\right]}^T=M_{6\times 6}{\left[X_1,Y_1,X_1{Y}_1,X_1^2,Y_1^2,1\right]}^T$$

(5)

Now, the transformation can characterize not only the translation of the grid points but also the flexible bending deformation of the original orthogonal grid (Figure 1c), which has the potential to characterize the movement and shape change of the weather echo. Therefore, Equation (5) is used as the form of nonlinear grid transformation discussed in the following sections.

2.2. Estimation of the Nonlinear Grid Transformation Matrix

As mentioned in Section 1, it is difficult to extract the target feature points that do not change with time in the weather echo image, so a pair of known and corresponding [X₁, Y₁] and [X₂, Y₂] is hard to find. Thus, a major challenge is how to solve Equation (5) to obtain the transformation matrix M_6×6 when given two or more echo images.

Assuming the original image values (reflectivity values) at two adjacent times are Q₁ and Q₂, [X₁, Y₁] is the original orthogonal grid of the image, and [X₂, Y₂] is the point set after transformation by M_6×6. Since [X₂, Y₂] probably does not coincide with the original orthogonal grid, let Q_2* be the estimated value on [X₂, Y₂] of Q₂. Here, Figure 1c can be taken as a simple schematic diagram, where the blue grids are [X₁, Y₁] and the red grids are [X₂, Y₂]. The pixels of Q₁ and Q₂ are 2D fields corresponding to the blue grids, while the pixels of Q_2* correspond to the red grids. The goal of the expected transformation is to make Q₁ look like Q₂, where the difference between Q₁ and Q_2* is as small as possible after moving the image values of Q₁ from [X₁, Y₁] to [X₂, Y₂]. The above objective can be expressed as minimizing Equation (6):

$$J=\sum {\left(Q_1-Q_{2*}\right)}^2$$

(6)

where the summation is for all grids of the image. However, Q_2* cannot be obtained by interpolation since the transformation and [X₂, Y₂] are unknown in advance. Thus, an estimate of Q_2* is proposed using a Taylor first-order expansion as follows:

$$Q_{2*}\approx Q_2+\frac{\partial Q_2}{\partial x}\left(X_2-X_1\right)+\frac{\partial Q_2}{\partial y}\left(Y_2-Y_1\right)$$

(7)

The approximation of Equation (7) is clearly rough and does not always hold, so some additional conditions are needed. First, the grid spacing should be larger than the displacement of the object in the images. Considering that the jet velocity in convective weather is often larger than 20 m s⁻¹, as is the outflow formed by convective precipitation, clouds may move 4.8~7.2 km or more with the wind within the interval of an operational weather radar, whose volume scan often takes 4~6 min. Therefore, the grid spacing selected for echo extrapolation should be greater than this order of magnitude, for example, at least 10 km × 10 km. Then, the echo image should be smoothed to some extent, such as via small-scale two-dimensional Gaussian filtering, to smooth the local jagged texture that may exist in the image so that the first-order partial derivative of the space in Equation (7) is representative.

After substituting Equation (7) into Equation (6), the problem of minimizing (6) can be transformed into solving a system of equations on the original orthogonal grid, where each equation is:

$$\Delta Q_{2-1}+\frac{\partial Q_2}{\partial x}\left(X_2-X_1\right)+\frac{\partial Q_2}{\partial y}\left(Y_2-Y_1\right)=0 \left(\Delta Q_{2-1}\ne 0\right)$$

(8)

where ΔQ₂₋₁ is the difference of Q₂ minus Q₁ on each grid. The expressions of X₂ and Y₂ can be obtained by expanding Equation (5) as follows:

$$X_2=m_{1,1}{X}_1+m_{1,2}{Y}_1+m_{1,3}{X}_1{Y}_1+m_{1,4}{X}_1^2+m_{1,5}{Y}_1^2+m_{1,6}$$

(9)

$$Y_2=m_{2,1}{X}_1+m_{2,2}{Y}_1+m_{2,3}{X}_1{Y}_1+m_{2,4}{X}_1^2+m_{2,5}{Y}_1^2+m_{2,6}$$

(10)

After substituting Equations (9) and (10) into Equation (8), the nonhomogeneous linear equations of the parameters to be solved can be obtained:

$$B M^{*}=\frac{\partial Q_2}{\partial x}{X}_1+\frac{\partial Q_2}{\partial y}{Y}_1-\Delta Q_{2-1}$$

(11)

where the known quantity B and the parameter M* to be solved are shown in Equations (12) and (13) as follows:

$$B=\left[\frac{\partial Q_2}{\partial x}{X}_1, \frac{\partial Q_2}{\partial x}{Y}_1, \frac{\partial Q_2}{\partial x}{X}_1{Y}_1, \frac{\partial Q_2}{\partial x}{X}_1^2, \frac{\partial Q_2}{\partial x}{Y}_1^2, \frac{\partial Q_2}{\partial x},\frac{\partial Q_2}{\partial y}{X}_1, \frac{\partial Q_2}{\partial y}{Y}_1, \frac{\partial Q_2}{\partial y}{X}_1{Y}_1, \frac{\partial Q_2}{\partial y}{X}_1^2, \frac{\partial Q_2}{\partial y}{Y}_1^2, \frac{\partial Q_2}{\partial y}\right]$$

(12)

$$M^{*}={\left[m_{1,1}, m_{1,2}, m_{1,3}, m_{1,4}, m_{1,5}, m_{1,6}, m_{2,1}, m_{2,2}, m_{2,3}, m_{2,4}, m_{2,5}, m_{2,6}\right]}^T$$

(13)

As long as more than 6 grids are satisfied, the system of equations consisting of Equation (11) can be solved by Least Squares Estimation (LSE). The obtained M* contains the elements of the first two rows of M_6×6, and [X₂, Y₂] can then be obtained by Equation (14):

$${\left[X_2,Y_2\right]}^T=\left[\begin{array}{c}{m}_{1,1}, m_{1,2}, m_{1,3}, m_{1,4}, m_{1,5}, m_{1,6}\\ m_{2,1}, m_{2,2}, m_{2,3}, m_{2,4}, m_{2,5}, m_{2,6}\end{array}\right]{\left[X_1,Y_1,X_1{Y}_1,X_1^2,Y_1^2,1\right]}^T$$

(14)

For a single transformation, the 2 × 6 transformation matrix in Equation (14) is sufficient. For the remaining 3~5 rows in M_6×6, each row can be obtained by LSE after substituting [X₂, Y₂] into Equation (5). Thus, a method for solving the nonlinear transformation matrix is obtained when two images are given. This approach is also applicable to multiple images. For example, for the case of three images, Equation (11) on each grid using the first and second images is first listed in the system of equations, then that of the second and third images are appended to the system of equations, and finally, they are solved together via LSE. This process enables the final transformation matrix to reflect the continuous movement and deformation characteristics of multiple images.

There are four more points to note:

(1)

Quadratic coordinates and distance derivatives are involved in the operation of Equation (11), which makes it possible for the value of the coefficient to span many orders of magnitude. Therefore, all the variables involved in the actual calculation should be at least double-precision floating-point numbers.

(2)

An underestimation of image change will occur if points of ΔQ₂₋₁ = 0 are considered when preparing Equation (11) because when only points of ΔQ₂₋₁ = 0 are used, the final obtained M_6×6 will be a unit matrix, that is, the image is unchanged.

(3)

Although the time complexity for solving Equation (11) appears to increase quadratically with the number of grids, the size of the original image does not become a problem because the method is not based on detailed texture and rather requires a certain degree of coarse grid spacing. A current personal computer can complete the calculation in 1–10 s when the historical sequence length is approximately 10, and the number of grid points is less than 100 × 100.

(4)

After obtaining the transformation matrix and before extrapolating using Equation (14), the image and coordinates are replaceable. For example, the filtered coarse-resolution image can be replaced by the original fine-resolution image. However, the unit of coordinates cannot be replaced because the transformation matrix is based on the coordinates used. For example, the coordinates cannot be changed from rectangular coordinate distance to latitude and longitude once the transformation matrix is obtained or from the latitude and longitude of a region to another different region.

2.3. Practical Steps for Extrapolation Prediction

An ideal extrapolation prediction scheme is to first find M_6×6 and obtain the transformation matrix after multiplication according to the time N that needs to be extrapolated, which is similar to Equation (3). Then, the inverse matrix is produced, and the discrete point set of reverse extrapolation is obtained and used for backward interpolation. Finally, the predicted image is obtained directly by linear interpolation. However, during the experiment for constructing the new method in this paper, some anomalies caused by nonlinear effects when N was large (such as greater than 10) were found. For example, the image does not further deform or move as expected but rolls back. Therefore, a more reliable extrapolation method is proposed.

First, only M_2×6 is obtained using Equation (11). Then, [X₂, Y₂] are obtained using Equation (14), and [X₁, Y₁] are subtracted from [X₂, Y₂] to find a vector field u and v to represent the one-time grid transformation. According to the time N to be extrapolated, the original grid points are moved in the opposite direction of u and v in N steps, and u and v on the new discrete points are interpolated in each step, which is known as backward interpolation. The values of the final backward positions are obtained by linear interpolation and are deemed to appear on the orthogonal grid point after N steps; thus, the predicted image is obtained. Although this process is somewhat different from a direct transformation similar to Equation (3), it can reflect the grid movement and deformation characteristics generated by the solved transformation matrix, and there is no abnormality caused by nonlinear effects.

The above introduces the image extrapolation method proposed in this paper, called the Nonlinear Grid Transformation (NGT) method. The main steps of echo extrapolation prediction are supplemented and summarized as follows:

(1)

Prepare a set of time-continuous historical echo images prior to the start time of the prediction. Perform resolution reduction and two-dimensional Gaussian filtering if necessary.

(2)

For each pair of adjacent images, each Equation (11) from a grid in which ΔQ₂₋₁ ≠ 0 is selected continuously to form a system of linear equations.

(3)

Solve the system of linear equations in the last step to obtain the transformation matrix M_2×6. Obtain [X₂, Y₂] according to Equation (14), and subtract the origin orthogonal grid [X₁, Y₁] to obtain grid transformation vector fields u and v. Here, the u and v are generated one-time.

(4)

If a resolution reduction is applied in the first step, u and v are linearly interpolated to the original resolution.

(5)

Based on the original image at the start time of the prediction, u and v are used to conduct backward interpolation, and the extrapolation image at N steps is obtained.

3. Weather Radar Echo Prediction Experiments

3.1. Experimental Setting

3.1.1. Case Setting and Statistics

To examine the proposed NGT method, a set of ideal cases and three real cases are used. In the ideal experiment, a simple two-dimensional normal distribution image with an extreme value region is used to simulate the echo of a precipitation cloud. Some simple changes, such as moving or numerically changing, are added to the images to test whether NGT can achieve the expected results. The three experiments using real radar observations are based on different radars and convective weather events. The original input echo images have a range of 600 km × 600 km and a resolution of 1 km × 1 km. The data form of these images is Constant Altitude Plan Position Indicator (CAPPI) or composite reflectivity, whose details are provided in Section 3.3, Section 3.4 and Section 3.5.

For real cases, three commonly used skill scores, namely, the critical success index (CSI, also known as the threat score), the probability of detection (POD), and the false alarm ratio (FAR), are used to examine the prediction results [25]. These skill scores are calculated as follows:

$$\mathrm{CSI}=\frac{{\mathrm{N}}_{\mathrm{A}}}{{\mathrm{N}}_{\mathrm{A}}+{\mathrm{N}}_{\mathrm{B}}+{\mathrm{N}}_{\mathrm{C}}}$$

(15)

$$\mathrm{POD}=\frac{{\mathrm{N}}_{\mathrm{A}}}{{\mathrm{N}}_{\mathrm{A}}+{\mathrm{N}}_{\mathrm{C}}}$$

(16)

$$\mathrm{FAR}=\frac{{\mathrm{N}}_{\mathrm{B}}}{{\mathrm{N}}_{\mathrm{A}}+{\mathrm{N}}_{\mathrm{B}}}$$

(17)

where N_A is the number of grid points with successful predictions (the hit number) for which both the observation and prediction are greater than a given threshold; N_B is the number of grid points with false alarms for which the observation is less than the given threshold, and the prediction is greater than the given threshold; and N_C is the number of grid points with unsuccessful predictions, for which the condition is opposite to N_B. Referring to previous studies [18,19], the thresholds of radar reflectivity (Z) on the echo image are set to 20, 30, and 40 dBZ. A perfect prediction result has a CSI and POD close to 1 and a FAR close to 0. CSI is the overall statistic, while POD and FAR help to further explain, for example, whether a lower CSI is due mainly to a low hit rate or a high false alarm rate.

When comparing the CSI of different methods, a static reference value is additionally calculated, where the CSI of this static reference is obtained without processing the image at the starting time of the prediction. The purpose of doing so is to test the effectiveness of the CSI in the prediction results. Specifically, for clouds with a large range or inconspicuous movement, higher CSI may be obtained even without any echo extrapolation. Therefore, an effective prediction result should have a CSI larger than the static reference at the same time; otherwise, it should be regarded as a quantitative result with a negative effect.

3.1.2. The Method for Comparison

A method that can also output the image moving vector field is to be selected for comparison to examine whether the image prediction is improved by the newly proposed NGT method. Note that TREC and OF can both fulfill that. However, OF has lower computational complexity than TREC at the same resolution, and some recent studies indicate that OF makes better predictions [12]. Therefore, a simple and feasible OF scheme is adopted. A global smoothing constraint condition of Horn-Schunk OF [10,26] is implemented, and the minimizing function is Equation (18):

$$J=\underset{\mathsf{\Omega }}{\overset{}{{\displaystyle \int }}}\left[{\left(I_{x}u+I_{y}v+I_t\right)}^2+{\alpha }^2\left({\left|\nabla u\right|}^2+{\left|\nabla v\right|}^2\right)\right]d\mathsf{\Omega }$$

(18)

where Ω is the area domain of grid points, I_t, I_x, and I_y are the temporal and spatial partial derivatives that synthesize two adjacent times, and α² is a smoothing parameter. The iterative scheme to solve the optical flow vector fields u and v is as follows:

$$u^{n+1}={\overline{u}}^n-I_x\frac{I_x{\overline{u}}^n+I_y{\overline{v}}^n+I_t}{{\alpha }^2+I_x^2+I_y^2}$$

(19)

$$v^{n+1}={\overline{v}}^n-I_y\frac{I_x{\overline{u}}^n+I_y{\overline{v}}^n+I_t}{{\alpha }^2+I_x^2+I_y^2}$$

(20)

where $\overline{u}$ and $\overline{v}$ are the regional average values of u and v, respectively, and the superscripts n and n + 1 represent the order of iterations. Here, α² is set to the global average of I_x² + I_y². The number of iterations is 1000 to ensure the solution instead of setting a breakdown condition. More details about the calculation scheme of I_t, I_x, I_y, $\overline{u}$, $\overline{v}$ and the boundary are provided in [26].

The OF adopts a similar method as the NGT, where the vector fields u and v for extrapolation are obtained first, and then backward interpolation is used on the image at the predicted time. For real cases, the solving of u and v is performed under an average grid of 10 km × 10 km. An additional two-dimensional Gaussian filter with a standard deviation of 1 grid and a radius of 10 grid is applied to NGT. The prediction is under the resolution of the original 1 km × 1 km grid. Images of the prediction start time and the previous time are used as historical samples in OF, while NGT uses historical samples one hour ahead of the prediction start time. The impact of historical sample selection on OF is briefly discussed in Section 3.2, and the impact of historical sample length on NGT is briefly discussed at the end of Section 3.3.

3.2. Ideal Cases

As shown in Figure 2a, both NGT and OF can make basically accurate predictions when only a simple translation exists. This also shows that the OF extrapolation scheme used for comparison in this paper is reliable for simple prediction.

When the image value has a time increment (Figure 2b), although the results of the two methods are both biased, the result of NGT is closer to the target object, while OF has almost more than twice the deviation, possibly due to the time increment, which does not conform to the basic assumption of the OF method itself.

When the moving path of the target object has an arc with rotation characteristics (Figure 2c), the result of NGT is closer to the target, while the result of OF seems to be only a linear extrapolation of the historical sample trajectory, resulting in an obvious deviation. In addition, it can be expected that the change in the historical sample selection will produce different OF results, where the predicted red circle will appear at the extension lines of the central positions of the two historical samples, which is still farther from the target. That is, the OF method is weak at extrapolating this path turning, regardless of how historical samples are selected. Note that the use of more than two historical samples for OF in [12] is expected also not helpful under this condition since it does not change the algorithm feature. In summary, compared with the conventional OF method, NGT can perceive more image changes from historical samples and can obtain more accurate extrapolation prediction results when the target value and path change.

The deviation of the NGT should also be noted. It is seen from Section 2.2 that only the spatial derivatives and coordinates of the image are involved in NGT, while the simple target used in Figure 2 has similar spatial derivatives in all directions, which might bring confusion and cause deviation. However, since this defect does not exist in real cases, it is no longer discussed below.

3.3. Real Case 1: A Convective Line with Shape Evolution

A convective weather event with obvious movement and shape change in radar echo images is selected first. A squall line event occurred in the U.S. on 20 May 2011, which was analyzed in a previous study [27]. The stage is selected in which a convective line is detected by a NEXRAD S-band radar (the station code is KTLX, located at 97.28°W, 35.33°N) when this convective line with shape-shifting moves to the radar. This radar performs a volume scan consisting of 15 elevations in 255 s. The CAPPI at a height of 3 km is obtained by triple linear interpolation from the volume scan data, which can not only reflect the main movement of the echo but also automatically shield most of the clutter in the boundary layer.

As shown in Figure 3, the convection line and its strong echo front are first in an inverted L-shape. As the cloud system moves northeastward, the convection line gradually evolves into a bow shape. Due to the north side of the cloud body moving out of the radar detection range, the original inverted L-shaped east–west echo band also gradually disappears. The visual moving directions of cells, new cells emerging and the whole cloud system are marked in Figure 3a. The horizontal wind U and V from ERA5 reanalysis data (Figure 4) are cited as a physical wind field for comparison with NGT and OF to analyze their features. The direction of new cell emergence is toward the southeast side, which is the windward direction of the relatively low level (e.g., 950 hPa), while the direction of the cell movement at a height of 3 km is along the north and slightly eastward, basically the same as that of the wind field at a similar height (from 850 hPa to 700 hPa). The combination of the above two directions eventually makes the whole cloud system move northeast.

Taking a 60 min prediction as an example, the data in Figure 3b are taken as the starting time, and Figure 3c is taken as the prediction target where the convection line moves and deforms. The results of NGT are basically consistent with the observed front edge of the convective line (Figure 5a). The grid transformation vectors mainly point to the northeast, which is consistent with the overall movement direction of the cloud system. In addition, these vectors point to the west at the northernmost point, which is also consistent with the image change characteristics of the weakening and disappearance of the east–west echo band on the north side in Figure 3. On the other hand, the OF results obviously underestimate the movement of the convection line (Figure 5b). The vector value of OF on the southeast side is small, which may be related to the interference of small-scale scattered clouds in the clear sky region at the front. More importantly, although the vector field of OF is similar to the wind field at the front of the convective line, pointing north and similar to the direction of monomer movement, this is not the direction of the overall cloud movement, which may be the main reason for the underestimation of movement.

Although the prediction effect decreases with time, there are differences between the two methods (Figure 6). The CSI of NGT is lower than that of OF only in the first 1–3 times of prediction and is better than that of OF in other prediction periods. For the higher 40 dBZ echo, NGT has the most obvious advantage over OF, and the CSI can be 0.2 higher than that of OF in the 40~60 min prediction period. The CSI of OF for 30 and 40 dBZ is lower than the static references within 60 min, which means that the prediction accuracy loses its value, while the CSI of NGT decreases more smoothly, indicating that it has the capability to perform echo prediction over a longer prediction period. In addition, due to the overall deviation of the position of the convective line predicted by OF, its prediction skill scores are all lower than those of NGT (

Table 1. Skill scores of the two methods for 60 min prediction. The target time is 10:05 UTC on 20 May 2011.

Threshold (dBZ)	POD		FAR		CSI
	NGT	OF	NGT	OF	NGT	OF
20	0.84	0.75	0.23	0.26	0.68	0.60
30	0.41	0.29	0.38	0.44	0.33	0.23
40	0.42	0.08	0.61	0.90	0.25	0.04

).

Figure 7 briefly shows the effect of historical sample length on NGT. The CSI is lower when only two historical samples (less than 5 min in this case) are taken. With an increase in the length of the historical sample, the CSI rises in oscillation first and then rises slowly when the sample is longer than 20~30 min, indicating that a longer sample length leads to a better NGT prediction effect. However, since the duration of convective weather is usually only a few hours or less, it is not practical to use historical samples longer than 60 min. Therefore, only 60 min of the historical sample for NGT is used and discussed in most parts of this paper.

3.4. Real Case 2: A Convective Event with Convective Line Formation

To examine whether the NGT is always better than OF in different stages of a convective event, it is necessary to select a case that has been detected by radar for a more complete life cycle. The selected case is a convective event that occurred in Shandong Province, China, on 17 May 2020. The radar data are from a CINRAD-SA type S-band radar (the station number is Z9532, located at 120.23°E, 35.99°N). This radar performs a volume scan consisting of 9 elevations within 342 s. Other information is provided in [28]. The CAPPI at 3 km height is used for analysis, similar to the last case in Section 3.3.

As shown in Figure 8, scattered convective clouds appear at the northwest side of the radar first, develop into a severe convective cell group, and then form a convective line, passing over the radar origin and moving southeastward. In the following time, the stratiform cloud area at the back of the convection line gradually expands, and the strong echo at the front gradually weakens, which is no longer displayed together here.

Figure 9 shows the 60 min prediction skill scores of almost the entire weather process. Clearly, NGT is not always superior to OF in echo prediction for 20 and 30 dBZ. However, the CSI of NGT for 40 dBZ is higher than that of OF at the target time of 14~16 UTC, which is the stage where the convective line exists.

A one-hour prediction is taken as an example for analysis similar to Section 3.3, where the data in Figure 8c are taken as the starting time, and Figure 8d is taken as the prediction target when the convection line persists. Note that 57 min is used here to represent a one-hour period since there is no data point exactly at the next 60 min. The NGT results show that the predicted convective lines are basically consistent with the observations (Figure 10a). Note that there is a large grid transformation vector over the normal wind speed in the clear sky area on the southeast side of the convective line pointing to the southeast, which indicates that the image is greatly stretched southeastward in the prediction. However, in the OF results (Figure 10b), an underestimation of the movement of the convective line remains. This underestimation is similar to that in Section 3.3 and is thus also due to the deviation of the moving direction of the whole cloud system and cells. In this 60 min prediction example, the CSI of NGT is better than that of OF (Figure 11), and the advantage is most obvious for 40 dBZ. Regarding the other skill scores, only the FAR of NGT is inferior to that of OF for 20 and 30 dBZ (

Table 2. Skill scores of the two methods for 60 min prediction. The target time is 15:22 UTC on 17 May 2020.

Threshold (dBZ)	POD		FAR		CSI
	NGT	OF	NGT	OF	NGT	OF
20	0.61	0.47	0.15	0.08	0.55	0.45
30	0.47	0.34	0.27	0.25	0.40	0.31
40	0.33	0.08	0.55	0.83	0.24	0.06

). This may be due to the increase in the false alarm rate of NGT caused by a displacement greater than the expected echo at X = 0~100 km and Y = −200~−100 km in Figure 10a, but this does not affect the overall advantage of NGT.

3.5. Real Case 3: A Mesoscale Convective Complex and Comparison with Previous Deep Learning Results

As a newly proposed method, NGT must be tested in different types of convective weather, and it is best to compare it with other newer popular methods. We have identified a case that can approximately meet this goal. A convective weather process dominated by a Mesoscale Convective Complex (MCC) occurred in Hebei Province, China, on 21 June 2017. A 30 min prediction result was presented in previous deep learning studies by Liang et al. [18] and Han et al. [19], where 8 years of data were used for training and 30 min of historical data were inputted to generate predictions. They mainly used a U-Net model, which is a convolutional neural network and is constructed by stacking downsampling and upsampling convolutional modules yielding a unique U-shape network architecture. The table of skill scores was given in [18] for this prediction example. Therefore, it is possible to compare NGT with the previous deep learning method by collecting the same radar data and repeating only the 30-min predictions during the same prediction period instead of collecting a large number of training datasets to reproduce the deep learning result or directly citing the resulting images, which is not applicable. The radar echo image used in [18] is composite reflectivity, which is a combination of six operational radars. The station numbers are Z9010, Z9220, Z9311, Z9313, Z9314, and Z9335, where Z9313 and Z9314 are CINRAD-CB type, and others are CINRAD-SA type. The time interval of the image is 6 min.

The observation at the target time and the prediction results of NGT and OF are shown in Figure 12, while the results of the deep learning method under the same prediction conditions are seen in [18]. The differences between the observation and the prediction results from both NGT and OF seem no longer as obvious as shown in the above two sections. For the CSI (Figure 13), the results of NGT and OF are very close, where NGT is slightly higher for 20 and 30 dBZ. In addition, the CSIs of NGT and OF are only slightly higher than the static reference, indicating that the prediction of echo movement is limited. Considering that the OF method is sensitive to the instantaneous movement of the cloud, these results indicate that the movement itself of MCC is not as obvious as the squall lines referred to in the last sections. More importantly, the deep learning method appears to be a more advanced method with results better than those of NGT and OF for 40 dBZ; however, its CSI for 20 and 30 dBZ is not only worse than that of NGT and OF but also worse than the static reference, which indicates that it may not correctly reflect the movement of the echo.

As shown in the skill scores (

Table 3. Skill scores of two methods for 30 min prediction and the comparison with results from a previous deep learning study. The target time is 11:30 UTC on 21 June 2017.

Threshold (dBZ)	POD			FAR			CSI
	NGT	OF	Liang 2021 *	NGT	OF	Liang 2021 *	NGT	OF	Liang 2021 *
20	0.75	0.77	0.84	0.13	0.18	0.33	0.67	0.66	0.54
30	0.61	0.63	0.67	0.27	0.32	0.49	0.49	0.49	0.39
40	0.28	0.31	0.51	0.62	0.62	0.57	0.19	0.21	0.30

* The results of “Liang 2021” are cited and reorganized from a previous study [18] as a deep learning method for comparison based on the best scores derived from the U-Net model and TrajGRU model.

), although the POD of this deep learning method at 20 and 30 dBZ is the highest, the FAR is also higher than that of the other two methods. That is, both the hit rate and false rate are high, indicating an overestimation of the reflectivity of the echo images. This overestimation is seen in the image in Liang et al. (2021) and might be one of the reasons that the deep learning method obtains a better skill score for a larger echo threshold.

4. Conclusions

The NGT method is proposed for weather radar convective echo extrapolation prediction. A set of ideal experiments and three real echo extrapolation experiments were performed to demonstrate the performance of the new method. The results are compared with those from an optical flow (OF) method, and a brief comparison with a previous study based on a deep learning method is conducted in a real case. The main conclusions are as follows.

(1)

In the ideal experiments for simple targets, NGT can achieve better prediction results than OF when there are changes in the value and path of the image target.

(2)

For cases with a convective line, NGT is superior to OF, where the CSI for 40 dBZ of the echo prediction at 60 min is approximately 0.2 higher.

(3)

For the case of MCC, NGT is better than OF and a deep learning method for 20 and 30 dBZ, while the deep learning method produces the best skill score for 40 dBZ. However, the good performance of the deep learning method may be due to the overestimation of the stronger echo.

The conventional OF method is based on a few historical data samples, tending to present the instantaneous movement direction of convective cells. Therefore, although the results of OF are closer to the physical motion of a cloud, there is a deviation compared with the movement of the overall cloud system. On the other hand, NGT can be modeled based on historical 60 min data. Although the grid transformation vector of the NGT is generally different from the physical wind field, the results are closer to the continuous variation characteristics of the echo image. Therefore, the movement and variation pattern of the whole convection line can be predicted more approximately.

In general, NGT, as a new method, improves the prediction effect of at least one type of convective system. This improvement in nowcasting capability can help to respond to severe convective weather, such as making corresponding changes to outdoor activities and transportation plans and better arranging artificial hail suppression operations.

However, many studies must be conducted in the future. The sensitivity of the prediction effect to the parameters in the prediction method needs to be further studied. More importantly, the CSIs of NGT and OF for stronger echo are still low and are markedly less than those for weaker echoes, which may be caused by the spatial nonuniformity and the strength variation of convective echoes. Therefore, an extrapolation of the echo strength may be needed to improve the result. Moreover, since deep learning methods have the potential to better predict the increase in echo strength, more operational observation cases are needed for further study to take full advantage of different prediction methods.