1 Introduction

Nowcasting usually refers to the 0–6-h forecasting of thunderstorms and accompanying catastrophic convective weather. At present, operational 0–1-h nowcasting is mainly based on thunderstorm detection and automatic extrapolation forecast technology using Doppler radar observations (Yu et al., 2012). Objective extrapolation of meteorological radar data for 30 min, 60 min, and longer predictions of thunderstorm location is often more effective than subjective extrapolation (Yu et al., 2012). In the 1950s, Ligda (1953) pioneered in the study of radar echo extrapolation. Since then, automatic radar echo extrapolation technology has been developed continuously.

Major extrapolation methods include the single-centroid tracking method, cross-correlation method, and optical flow method. In the single-centroid tracking method (Crane, 1979; Dixon and Wiener, 1993; Johnson et al., 1998), a storm is identified, tracked, and extrapolated as a three-dimensional monomer for prediction. In the cross-correlation method, correlations of radar echoes at consecutive times are computed to establish the best fitting relation of radar echoes (Rinehart and Garvey, 1978; Wilson et al, 1998; Chen et al., 2007) and obtain the echo motion vectors, which are used for extrapolation to determine radar echo locations and patterns. The single-centroid tracking method is mainly used for precipitation forecasting by means of extrapolation. The cross-correlation method can be employed to predict convective precipitation, and it is also able to track stratiform precipitation. Thereby, this method has been widely applied in operational weather forecasting (Rasmussen et al., 2001; Mueller et al., 2003). However, the failure incidence of extrapolation via the cross-correlation method significantly increases for locally generated echoes and echoes that change rapidly in intensity and pattern (Han et al., 2008). In recent years, application of the optical flow method for nowcasting has gradually increased and its performance is regarded as satisfactory (Han et al., 2008; Cheung and Yeung, 2012; Cao et al., 2015).

An evaluation of precipitation forecasting in South China (Cao et al., 2015) suggested that the optical flow method was better than the cross-correlation method for forecasting locally generated echoes and thunderstorm echoes with rapid changes in intensity and pattern. For tropical precipitation systems, however, forecasts produced by the cross-correlation method were better than those derived from the optical flow method. This implies that, to a certain extent, the optical flow method remains limited. The likelihood that we can apply high spatiotemporal resolution numerical products for operational nowcasting in the near future is still not high (Zheng et al., 2015). Thereby, nowcasting will rely heavily on the extrapolation of radar data (Wilson et al., 2010), at least before 2020. However, how to combine various nowcasting methods to achieve optimum forecasting is a challenging issue at present. To overcome the shortcomings of nowcasting and make up for the insufficiencies of nowcasting techniques for 0–1-h forecasts, a particle filter blending algorithm is developed in the present study on the basis of the optical flow method to obtain the optimum estimation of the echo motion. This estimate is then used for nowcasting to improve the ability of 0–1-h forecasting.

2 Data and method 2.1 Data

Radar reflectivity mosaics at an altitude of 2.5 km from the 12 S-band radars in Guangdong Province, China (specifically, at Guangzhou, Shenzhen, Shaoguan, Zhuhai, Qingyuan, Yangjiang, Heyuan, Shanwei, Shantou, Meizhou, Zhanjiang, and Zhaoqing), are used in the present study. For the convenience of computation, radar data are interpolated from the polar coordinate system to the three-dimensional Cartesian coordinate system by using the nearest neighbor method in the horizontal direction and the linear interpolation method in the vertical direction (Xiao and Liu, 2006).

2.2 Radar data quality control

The obtained radar base data contain real information on precipitation echoes and false information from various external noises. Such noises can significantly influence the track extrapolation of echoes. In the process of radar echo tracking, the particle filter algorithm requires high quality echo data as the input. Thereby, in addition to the above quality control, a much stricter quality control of echo noises is necessary to obtain more realistic smooth radar echoes and echo motion vector fields.

2.2.1 Quality control of radar echo noise—bilateral filtering

The purpose of filtering is to remove noises in radar echoes but retain as much real information as possible. The mean filter is one of the simplest filters, which performs a local averaging operation to replace each image pixel with the average of all the values in the local neighborhood. The weakness of this method is obvious. The median filter (Zhao et al., 2007) replaces each pixel value with the median of pixel values in the local neighborhood. Pixels in each window are sorted into ascending order and the pixel value in the middle is selected as the new value for a particular pixel. The median filter can effectively remove large variances caused by singularity at certain points, and is thus an important filter used in the quality control of radar base data. The disadvantage of the median filter is that it weakens the peak while reducing the noise, and is therefore unfavorable for tracking severe weather processes. The Gaussian filter (Si et al., 2014) is a convolution of each pixel of the input image with the convolution kernel, and the sum of the final results is taken as the pixel value of the output image. Spatially, radar reflectivity changes slowly and continuously, and thus the pixels in the local neighborhood will not change significantly within 6 min. However, noises are not spatially associated, often differing from one another in the local neighborhood. The Gaussian filter makes use of this feature to remove the noises but retain the real signals. The effect of the Gaussian filter is better than that of the median filter. However, strong echoes might be located at the edges of the echo belt for some severe weather processes, such as a squall line. In this situation, the Gaussian filter might smooth the echo edges, which can seriously affect the reliability of the estimated echo motion. The bilateral filter (Overton and Weymouth, 1979) is a weighted Gaussian filter, which retains the information at the echo edges while removing the noises. Thereby, it can overcome the disadvantages of Gaussian filter. Compared to the Gaussian filter, an extra variable, a Gaussian variance ${\sigma _{\rm s}}$ is added to the bilateral filter. The bilateral filter algorithm can be expressed as (Zhou et al., 2014)

${F\!_{\rm B}}\!\left(I \right) \!=\! \frac{{\sum\limits_{i\!,\, j = - w}^w \!\!\!\!{{G_{{\sigma _{\rm s}}}}\!\left({x, y;{x_i}, {x_j}} \right)\!{G_{{\sigma _{\rm r}}}}\!\left({x, y;{x_i}, {x_j}} \right)\!I\!\left({{x_i}, {x_j}} \right)} }}{{\sum\limits_{i\!,\,j = - w}^w {{G_{{\sigma _{\rm s}}}}\left({x, y;{x_i}, {x_j}} \right){G_{{\sigma _{\rm r}}}}\left({x, y;{x_i}, {x_j}} \right)} }}, $

(1)

where

${G\!_{{\sigma _{\rm s}}}}\left({x, y;{x_i}, {x_j}} \right) = \exp \left\{ { - \frac{{{{\left({x - {x_i}} \right)}^2} + {{\left({y - {y_j}} \right)}^2}}}{{2\sigma _{\rm s}^2}}} \right\},$

(2)

and

${G\!_{{\sigma _{\rm r}}}}\left({x, y;{x_i}, {x_j}} \right) = \exp \left\{ { - \frac{{I\left({x, y} \right) - I{{\left({{x_i}, {x_j}} \right)}}}}{{2\sigma _{\rm r}^2}}} \right\}.$

(3)

In Eq. (1), I and ${F_{\rm B}}\!\!\left(I \right)$ are the input image and the image after filtering, respectively. The Gaussian kernel function ${G_{{\sigma _{\rm s}}}}\left({x, y;{x_i}, {x_j}} \right)$ represents the spatial similarity of interior points in rectangles centered at $\left({x, y} \right)$ with a radius of w; ${\sigma _{\rm s}}$ is the variance parameter, which is set to 3; ${G_{{\sigma _{\rm r}}}}\left({x, y;{x_i}, {x_j}} \right)$ represents the pixel similarity of interior points in rectangles centered at $\left({x, y} \right)$ with a radius of w; and ${\sigma _r}$ is set to 0.1.

The above algorithm suggests that the bilateral filter is based on the weighting of the luminance difference between the surrounding pixels and the central pixel. The intensity value at each pixel in an image is replaced by a weighted average of intensity values from nearby pixels. Similar pixels have higher weights and dissimilar pixels have smaller weights. The bilateral filter can retain excellent signals while greatly reducing noises and fully preserving the echo edges.

2.2.2 Effects of bilateral filtering

Two different filters were applied to radar echoes on 13 April 2016. Before filtering, the edges of the strong echoes, the convective echoes within the strong echoes, and the echoes of some stratiform clouds, were scattered and fragmentary with serrated edges (Fig. 1a). After median filtering and bilateral filtering (Figs. 1b, c), the echo pattern, the strong echo center, and the echoes of stratiform clouds in the rear, became smoother and more orderly. Compared to those before filtering, the echoes at the front edge were weakened by the median filter but well preserved by the bilateral filter. The median filter improved the quality of the echo but weakened the echo edge, especially the echo edge of the squall line. As a result, wind vectors at the echo edges could only be roughly estimated in the computation of motion vectors of radar echoes, which would affect the estimation quality of echo motions. In contrast, the bilateral filter improved the echo quality and preserved the echo edges well. Thereby, the wind vectors at the echo edges could be obtained directly with relatively reliable quality in the computation of echo motions. The above analysis indicates that the effects of the bilateral filter are better than those of the median filter. After bilateral filtering, the echo patterns remained realistic and became smoother with a clear echo edge. The noises were effectively removed.

2.3 Introduction to the particle filter blending algorithm

The particle filter method was first proposed by Carpenter et al. (1997). It obtains recursive Bayesian (Bayes filter) estimation mainly by non-parametric Monte Carlo simulation. Bayesian estimation is the theoretical basis of particle filtering (Li et al., 2015). The core idea of the particle filter is to represent the posterior probability density distribution by using a random sample set (particle set) with corresponding weights (Arulampalam et al., 2002). Specifically, the probability density distribution is approximated first by a sample set of weighted particles, and the realistic distribution is then obtained by updating the weight of each particle in the sample based on measurements. Finally, importance sampling is conducted to achieve a uniform weight distribution (Lyu et al., 2013). The requirements of the particle filter on the system are relatively low, and previous studies (Wang G. et al., 2013) have shown that the accuracy of particle filter blending can approximately reach the optimal motion estimation. In recent years, the particle filter algorithm has become a hot research topic. It is widely applied in automatic control, navigation, and motion tracking (Wang et al., 2015). The particle filter algorithm is simple and easy to implement, which makes it suitable for trajectory blending of most nonlinear moving objects. It can also be used for the blending of echo radar motion vectors. Motion vectors of radar echoes are required for extrapolated forecasts of radar echoes. Different to that using a single set of motion vectors for weather forecasting, two sets of motion vectors are obtained from the same radar echoes with two different estimation methods. These two sets of motion vectors are then blended by using the particle filter algorithm and extrapolated for prediction. Details of the particle filter algorithm are described as follows.

In the particle filter algorithm, echo vectors can be written as (Wang et al., 2015)

$\left\{ {\begin{aligned}& {{x_t} = f\left( {{x_{t - 1}}} \right) + {m_{t - 1}}}\\& {{y_t} = h\left( {{x_t}} \right) + {n_t}}\end{aligned}} \right..$

(4)

Here, $f\left(x \right)$ and $h\left(x \right)$ are state equations and observation equations, respectively; and x_t, y_t, m_t, and n_t represent the state of the system, the observation, the process noise, and the observation noise, respectively. In this study, $f\left(x \right)$ is the motion vector obtained by the Lucas–Kanade optical flow (Cao et al., 2015). The basic method of the Lucas–Kanade optical flow for echo motion vector computation is that a small movement to $\left({x + \Delta x, y + \Delta y} \right)$ at $\left({t + \Delta t} \right)$ is assumed, the gray-scale intensity value remains unchanged during the time interval $\Delta t$ , and then the optical flow constraint equation can be expressed by

${I_x}u + {I_y}v + {I_t} = 0, $

(5)

where $u \!=\! \displaystyle\frac{{{\rm d}x}}{{{\rm d}t}}\!, v \!=\! \displaystyle\frac{{{\rm d}y}}{{{\rm d}t}}, {I_{\! x}} \!=\! \displaystyle\frac{{\partial I}}{{\partial x}}, {I_{\! y}}\! =\! \displaystyle\frac{{\partial I}}{{\partial y}}, {\rm and}\,{I_{\! t}} \!=\! \displaystyle\frac{{\partial I}}{{\partial t}}$ .

The Lucas–Kanade local constraint (Lucas and Kanade, 1981) is applied in the optical flow computation, which is computed as in Cao et al. (2015).

Moreover, $h\left(x \right)$ is based on the Harris corner detection algorithm (Gao et al., 2012) to track echoes and obtain motion vectors. In the study of motion tracking and estimation, corners provide rich information, and they can be detected from images at consecutive times. The definition of the Harris corner is based on the matrix of the second order derivative of image gray-scale intensity. A new second-order derivative image is obtained by deriving all the pixels of the original image. The Hessian matrix is written as

$H\left(P \right) = \left({\begin{array}{*{20}{c}}{\frac{{{\partial ^2}I}}{{\partial {x^2}}}}&{\frac{{{\partial ^2}I}}{{\partial x\partial y}}}\\{\frac{{{\partial ^2}I}}{{\partial x\partial y}}}&{\frac{{{\partial ^2}I}}{{\partial {y^2}}}}\end{array}} \right).$

(6)

The Harris corner is determined by using the autocorrelation matrix of the second derivative image within the neighborhood of each pixel. It can be expressed as

$M\left({x, y} \right) = \left[ {\begin{array}{*{20}{c}}{\sum\limits_{ - k \leqslant i, j \leqslant k} {{W_{i, j}}I_x^2\left({x + i, y + j} \right)} }&{\sum\limits_{ - k \leqslant i, j \leqslant k} {{W_{i, j}}{I_x}\left({x + i, y + j} \right){I_y}\left({x + i, y + j} \right)} }\\{\sum\limits_{ - k \leqslant i, j \leqslant k} {{W_{i, j}}{I_x}\left({x + i, y + j} \right){I_y}\left({x + i, y + j} \right)} }&{\sum\limits_{ - k \leqslant i, j \leqslant k} {{W_{i, j}}I_y^2\left({x + i, y + j} \right)} }\end{array}} \right].$

(7)

Here, ${W_{i, j}}$ is the weighting factor and $I\left({x + i, y + j} \right)$ is the second derivative image pixel. The Harris corner corresponds to where the autocorrelation matrix has the two largest eigenvalues. Since the second derivative is irresponsive to the uniform gradient, this method can take the strong echo center or rotated echoes as the Harris corner, while no corner can be detected in spatiotemporally uniform echoes. This indicates that the above method is appropriate for tracking strong echoes and rotated echoes, but other supplementary methods are necessary for tracking echoes of uniform motion. More details can be found in Chen et al. (2003).

After the observational value y_t is obtained, the conditional probability density function of the state variable x_t is estimated. For the radar echo tracking, assume that the state of the echo x_t at the current time is only relevant to that at the previous time ${x_{t - 1}}$ , and that the observational value y_t at the time t is only relevant to the state variable x, that is, observations are independent of each other.

The core idea of the particle filter is to obtain the posterior probability density distribution using a random sample set (particle set) with corresponding weights. The recursive Bayesian filtering theory provides a strict solution to this problem. The optimization algorithm is introduced to optimize the parameters of the particle filter, and the entire parameter space is efficiently searched to obtain the optimal solution. Repetitive iterations are conducted for the updated parameters to obtain the estimation of the optimal $\Omega \left({{x_t}\left| {{y_t}} \right.} \right)$ . Compared to other methods, the particle filter can obtain more optimal echo motion vectors. A series of spatially random samples are used in the particle filter to approximate the a posteriori probability density function, which is represented by a sample set of weighted particles ${X_k} = \left\{ {x_k^i, \varpi _k^i} \right\}_{i = 1}^N$ written as (Gong, 2010)

$\Omega \left({{x_t}\left| {{y_t}} \right.} \right) \approx {\Omega _N}\left({{x_t}\left| {{y_t}} \right.} \right) = \sum\limits_{i = 1}^N {\varpi _k^i\delta \left({{x_k} - x_k^i} \right)} .$

(8)

Here, $\varpi _k^i$ is the normalized weight that satisfies $\sum\limits_{i = 1}^N {\varpi _k^i \!=\!\! 1} $ , and N is the number of particles.

There are three major steps to approximate the a posteriori probability density function using the particle filter algorithm: the particles are sampled first from the prior distribution; the weights of the particles are then calculated; and finally, the samples are re-sampled to avoid particle degeneration. The procedure is as follows:

First, we need to initialize k = 0. Samples $x_0^i$ are collected based on the a priori density function $\Omega \left({{x_0}} \right)$ , and the weight of each particle is set to be $ \omega _0^i = 1/N, $ $ i = 1, \ldots, N$ . Next, when k ≥ 1, iteration is conducted following the steps below:

(1) Particle importance sampling. After particles $x_k^i \sim q\left({{x_k}\left| {x_{k - 1}^i} \right.} \right)$ are sampled, their corresponding weights are computed according to $\omega _k^i = \omega _{k - 1}^ip\left({{z_k}\left| {x_k^i} \right.} \right)$ and normalized by $\varpi _k^i = \omega _k^i/\sum\limits_{i = 1}^N {\omega _k^i} $ , i = 1, …, N. One key issue in the particle filter algorithm is the selection and re-acquisition of the particle importance distribution function. In the present study, a priori density is chosen as the importance distribution function, in which $q\left({ \cdot \left| \cdot \right.} \right)$ represents the probability density function that is easier to acquire, while $p\left({ \cdot \left| \cdot \right.} \right)$ represents the probability density function that is harder to acquire.

(2) Re-sampling. If the number of effective samples is smaller than thethreshold value ${N_{\rm {th}}}$ (where N = $1/{\sum\limits_{i = 1}^N {\left({\varpi _K^I} \right)} ^2}$ and ${N_{\rm {th}}}$ is the pre-specified threshold value, which is set to 3N/4 in this study), the particles $\left({x_k^i} \right)_{i = 1}^N$ are resampled to produce a new set of particle samples $\left({X_k^i} \right)_{i = 1}^N$ that satisfies the condition $P\left({X_k^i = x_k^i} \right) = \omega _k^i$ ( $P\left(\cdot \right)$ is the probability). Weight values are redefined as $\omega _k^i = 1/N, $ i = 1, ..., N. Otherwise, step (2) can be skipped and we can move from step (1) to (3) directly.

(3) Update of the state by

${x_k} = \sum\limits_{i = 1}^N {X_k^i\omega _k^i} .$

(9)

Recursive Bayesian filtering is used to optimize the two echo motion vector fields to seek the optimal solution of the posterior probability density function and obtain the best track of echo motions.

A realistic and smooth motion vector field is critical for a successful forecast by extrapolation. The core idea of cross-correlation is to determine the echo motion vector field based on the best fitting relation established between echoes of consecutive times. The optical flow method determines the echo motion vector field by the echo trajectory under certain constraints. The particle filter method blends the two types of echo motion vector field to approximate an optimal estimation of the echo motion and obtain a more realistic and more accurate motion vector field. The above three methods were employed to obtain the echo motion vectors of a squall line process at 0630 BT (Beijing Time) 13 April 2016, and the results are shown in Fig. 2. In the cross-correlation method and optical flow method, the median filter was applied for radar data quality control.

During the 1-h period from 0630 to 0730 BT, the squall line roughly moved along the northwest–southeast direction (figure omitted). The moving track predicted by the particle filter blending method was generally consistent with the observation. The observation indicated that the moving directions at various parts of the squall line were different. The southern part of the squall line moved southeastwards, the northern part moved more eastwards, and the southeastern part moved between the southeast and east. The motion vector field of the squall line determined by the particle filter blending method (Fig. 2a) suggested that the southern part of the squall line was under the control of northwesterly winds, the northern part was under the control of westerly winds, and the southeastern part was under the control of northwesterly winds. These results were basically consistent with the observed moving direction of the squall line. The particle filter blending method was also able to accurately describe the vector field. The results of the cross-correlation method indicated that westerly winds prevailed over the squall line (Fig. 2b), which implied an eastward movement of the squall line. The results of the optical flow method suggested that westerly winds prevailed over the western part of the squall line, southwesterly winds prevailed over the northern part (Fig. 2c), and the squall line was predicted to move northeastwards within an hour. Because of quality control, the motion vectors at the strong echo center or areas of high rotation would be smoothed by the surrounding motion vectors when tracking the echo motions via the cross-correlation and optical flow methods, which would lead to biases in the forecast of the echo motion path. For this reason, it is clear that the particle filter blending method can provide more realistic and accurate forecasts of echo motion vectors compared to the other two methods. In the particle filter blending method, the optimum estimation of the echo motion is obtained by combining various motion vectors.

2.4 Semi-Lagrangian extrapolation method

Extrapolated nowcasting of echo motions can be conducted after the motion vectors are obtained through the particle filter blending. Radar echo motions are usually not linear, and thus the semi-Lagrangian extrapolation method is often applied for predicting the echo motion, since this method can overcome the weakness in linear extrapolation that does not consider the rotation of the echo motion. Experimental results in previous studies indicate that the semi-Lagrangian method performs well in terms of both computation stability and accuracy, and thus has been widely applied in numerical weather forecasting and climate models (Tian et al., 2004; Zhang et al., 2015). The semi-Lagrangian advection can be expressed as (Wang G. L. et al., 2013)

$\overline F \left({{t_0} + \tau, x} \right) = F\left({{t_0}, x - \alpha } \right), $

(10)

where $\overline F $ and F are the forecast and observation of radar echoes, respectively; t₀ is the start time of the forecast; $\tau $ is the valid lead time of the forecast; and $\alpha $ is the displacement. In the semi-Lagrangian method discussed here, several extrapolations (N) are conducted over the valid time of the forecast ( $\tau $ ). The time step $\Delta t$ is set to 6 min, and forward extrapolation is utilized. The displacement $\alpha $ at each time step is obtained by iteration, as shown below (Wang G. L. et al., 2013):

$\alpha = \Delta tu\left({t_0, x - \frac{\alpha }{2}} \right), $

(11)

where $u\left({t_0, x - \displaystyle \frac{\alpha }{2}} \right)$ is the speed at grid $x - \displaystyle \frac{\alpha }{2}$ . The initial value of $\alpha $ is 0, and the total displacement is the sum of the displacements of the N time steps. Details of echo extrapolation can be found in Reich (2007).

The forecast map of radar echoes can then be obtained based on the semi-Lagrangian extrapolation of the echo motion vectors from the particle filter blending.

3 Case studies of extrapolated forecasts

At present, the cross-correlation method and optical flow method are used operationally in the nowcasting system of the platform of nowcasting decision supporting (PONDS) in the Meteorological Bureau of Shenzhen Municipality for nowcasting of severe convective weather in Guangdong. Like the cross-correlation method and optical flow method, the particle filter blending method is used for 60-min convective forecasts at 6-min intervals based on radar mosaics of 12 radars in Guangdong. After receiving all the radar data in Guangdong Province, the system starts to conduct quality control, produce radar mosaic data, calculate motion vectors, perform particle filter fusion blending and echo motion extrapolation, and so on. Finally, the forecast is presented in the forecast system. All the above processes from receiving the radar data to presenting the forecast results can be completed in 2.5 min. Therefore, the particle filter blending method can meet the requirements of nowcasting. Comparative experiments using the particle filter blending method, the cross-correlation method, and the optical flow method for nowcasting will be discussed next.

3.1 Squall line on 13 April 2016

From early to late morning of 13 April 2016, a powerful squall line moved from northwest to southeast and affected Guangdong. The 12 radars deployed in Guangdong recorded the complete evolution of this squall line generated by an upper-level trough and cold air mass. A 60-min extrapolated forecast of the squall line at 6-min intervals was produced by using the particle filter blending, cross-correlation, and optical flow methods. Radar echoes at 0630 BT were used as the initial field for the forecast. The results from the three different methods were compared.

Comparison of the above results indicates that the 60-min forecasts of the echo pattern, location, and area of strong echoes by the particle filter blending method (Fig. 3c) were consistent with observations at 0730 BT (Fig. 3b). The bow echoes in the southeastern section, the linear convective echoes in the northern section, and the strong echoes in the southern section of the squall line were well predicted. The location of the predicted echoes in the southern section of the squall line basically agreed with the observation. The area of strong echoes was overestimated by the cross-correlation method (Fig. 3d), and the echoes in the southern section of the squall line were located slightly to the north of observations. The area of strong echoes was underestimated by the optical flow method (Fig. 3e), and the predicted echoes in the southern and northern sections of the squall line were weaker than observations. The predicted echoes in the southern section of the squall line shifted more northwards than observed and the predicted moving speed of the squall line produced by the particle filter blending method was slower than observed. This was attributed to the slower movement of echo motion vectors predicted by the particle filter blending algorithm. Therefore, it is necessary to adjust the blending parameters.

By comparing Fig. 3b with Fig. 3c, it is clear that the particle filter blending method realistically predicted the characteristic location of the maximum reflectivity of the squall line. Strong echoes were largely located in the northern, southeastern, and southern sections of the squall line. However, biases also existed in the echo intensity forecast in various sections of the squall line. The forecast for the southern part of the squall line was close to observations, with the echo intensity within the range of 35–50 dBZ. The maximum forecast echo intensity in the southeastern section of the squall line was 53 dBZ, which was similar to the observation. The length of the strong echo belt was up to 150 km, which was long enough to cover the entire coastal region along both sides of the Pearl River estuary. Observed strong echoes were concentrated over the coastal region of Shenzhen, to the east of the Pearl River estuary. The length of the strong echo belt was relatively short in the observation while it was overestimated in the forecast. For the northern section of the squall line, the forecast echo intensity was stronger than observed by about 10 dBZ. Overall, the squall line intensity was overestimated by the cross-correlation method. The echo intensity predicted by the optical method was similar to the observation in the northern section of the squall line, but it was overestimated in the southeastern section and underestimated in the southern section. All three methods were able to successfully reproduce the echoes in the rear of the squall line corresponding to the trailing stratiform clouds. However, it was hard to predict the intensification, weakening, and disappearance of echoes since the mechanisms for echo intensity changes were not considered in the extrapolation.

Comparison of the 30-min forecast and observation of the squall line (figure omitted) suggested that the squall line location, echo pattern, intensity, and moving speed predicted by the particle filter blending method were all consistent with observations. Apparently, the extrapolated forecast produced by the particle filter method could provide valuable implications for 30-min forecasts of well-organized convective systems such as squall lines.

Compared to the cross-correlation and optical flow methods, the particle filter blending method performed better in the prediction of the echo pattern, intensity, and moving direction associated with the squall line process. It improved the extrapolated nowcasting of this squall line. This result indicates that it is practical to apply the particle filter method for extrapolated nowcasting.

3.2 Southwesterly-monsoon precipitation on 21 May 2016

On 21 May 2016, a heavy rainstorm and localized extremely heavy rainfall occurred in the coastal region of Guangdong due to the influence of intensified southwesterly monsoon. Rainfall was concentrated during the early to late morning of 21 May. The monsoon flow converged in the coastal region of South China, leading to heavy rainstorms and localized extremely heavy rainfall over large coastal areas from Yangjiang to Jieyang County. The maximum 24-h accumulative precipitation of 463 mm was recorded in Shenzhen, and the maximum hourly precipitation was 112.5 mm. This was a typical rainfall process solely induced by tropical system with no distinct cold air activities involved (Yu et al., 2012).

The radar echoes at 0200 BT were selected as the initial field for the extrapolated forecast. Comparing the radar mosaic at 0230 BT (Fig. 4b) for this rainfall process and the extrapolated forecast (Fig. 4c) produced by the particle filter blending method, the λ-shaped echo pattern over the ocean to the west of the Pearl River estuary, the two areas of strong echo nearby the Pearl River estuary, and the linear strong echoes in southeastern Shenzhen to the east of the Pearl River estuary, were all well predicted by the particle filter blending method. The predicted echo location, shape, and intensity were consistent with observations, except that the predicted echo shape in the eastern part of the squall line was thinner than observed. As shown in Fig. 4d, the convective echo area predicted by the cross-correlation method was larger than observed, and the echo intensity was slightly overestimated. The echo shape predicted by the optical flow method was similar to observations, but the echo intensity was underestimated (Fig. 4e).

After 0230 BT, the echoes intensified distinctly. Large biases of the echo location and shape were found in the 60-min extrapolated forecast by all three methods. Apparently, it was difficult for the extrapolation system to predict the evolution of the echoes. The above analysis indicates that the performance of the particle filter blending method was better than that of the cross-correlation and optical flow method for the extrapolated forecast of this southwesterly-monsoon precipitation event.

3.3 Local convective precipitation on 14 April 2016

In the morning of 14 April 2016, a local convective weather event developed in the Pearl River delta. The echo first occurred at around 0830 BT to the west of the Pearl River estuary in Jiangmen, and then slowly moved eastwards and gradually intensified. Two new echoes formed in the rear of the previous echo during its eastward movement, leading to slow-moving echoes. During this process, the thunderstorm moved slowly with heavy hourly precipitation. Comparative experiments of 60-min extrapolated forecasts of the echo motion were initialized at 1054 BT.

Figure 5 displays the 60-min extrapolated forecasts of this local heavy precipitation event using the three methods. Comparison of the echoes predicted by the particle filter blending method (Fig. 5c) with observations (Fig. 5b) indicated that the arc echo pattern in the northern section of the echoes was well predicted, while the northeast–southwest extending echo in the southern section was slightly overestimated. The predicted location, shape, and intensity were consistent with observations (Fig. 5b). However, certain biases could be found in the 60-min forecasts of echo location and shape by both the cross-correlation (Fig. 5d) and optical flow methods (Fig. 5e), and the echo intensity was also underestimated by the two methods.

Experiments of 30-min extrapolated forecasts (figure omitted) suggested that the particle filter method, cross-correlation method, and optical flow method yielded similar results. The predicted echo location and intensity were roughly consistent with observations. The above discussion indicates that it is practical to apply the particle filter blending method for extrapolated forecasts of locally generated echoes.

3.4 Precipitation of Typhoon Nida on 2 August 2016

Typhoon Nida, the fourth typhoon of the annual typhoon season in 2016, landed in Shenzhen at around 0400 BT 2 August. It then moved northwestwards across Shenzhen, bringing winds and precipitation weather to Guangdong. Figure 6 presents the 30-min extrapolated forecast of this typhoon process. The initial time of the forecast was 0430 BT.

Figure 6b shows radar echo observations at 0500 BT 2 August. The eye of the typhoon was located inside Dongguan. The eye area was a long clear strip and extended from southeast to northwest. Outside the eye area were large spiral areas of rainfall. Areas of high dBZ values were located to the east, south, and northwest of the typhoon center. Strong echoes to the east of the typhoon center covered the northeast–east of the first quadrant, strong echoes to the south covered the third and fourth quadrants, and the strong echoes to the northwest occupied the northwest–west of the second quadrant. Echoes were relatively weak in the northeast–north of the first quadrant and northwest–north of the fourth quadrant. Together they constituted the typhoon spiral rainbelt. The echo intensity within the spiral rainbelt was higher than 35 dBZ and the maximum intensity was greater than 40 dBZ. Outside the spiral rainbelt, about 400 km to the southwest of the typhoon eye, a northwest–southeast extending echo belt appeared with a maximum intensity greater than 40 dBZ, which corresponded to cumulus convection that was mainly related to the outer cloud system of Typhoon Nida.

The 30-min extrapolated forecast initialized at 0430 BT (Fig. 6a) predicted clearly the typhoon eye (Fig. 6c) that was also located inside Dongguan, and had a size consistent to that observed. The spiral rainbelt in the eyewall region agreed well with the radar-observed spiral rainbelt. Similar to the observations, the predicted strong convection in the eyewall region was also located to the south, east, and northwest of the typhoon eye. Areas of high echoes in the prediction were close to observations, with echo intensities greater than 35 dBZ. The strong northwest–southeast extending echo belt to the southwest of the typhoon, 400 km away from the eye, was also predicted. Compared to observations, the echo intensity predicted by the particle filter algorithm was roughly similar but the area of echoes was slightly overestimated. Strong echoes in the spiral rainbelt were also predicted by the cross-correlation method, but the echo intensity was overestimated (Fig. 6d). The 30-min prediction produced by the optical flow method was also satisfactory (Fig. 6e), although the prediction of the long northeast– southwest extending echo strip to the southeast of the typhoon eye and the echoes in the southwestern periphery of the typhoon was not as good as that by the particle filter algorithm. The above analysis suggests that the spiral rainbelt of the typhoon and the rainbelt in the periphery of the typhoon in the 30-min extrapolated forecast produced by the particle filter algorithm were all consistent with observations.

An experimental 1-h extrapolated forecast of Typhoon Nida was also conducted (figure omitted). The results indicated that the 1-h forecast of the spiral rainbelt and rainbelt in the periphery of the typhoon were roughly consistent with observations, as was the predicted echo intensity. Thereby, it is practical to apply the particle filter method to retrieve and track the spiral rainbelt of a tropical storm.

Comparative experiments for extrapolated forecasts of the above four weather events were conducted to assess the nowcasting effects of the three methods, that is, the particle filter blending method, the cross-correlation method, and the optical flow method. Analysis of these experimental results indicates that the three methods could predict the 30- and 60-min echo locations, shapes, and intensity of the thunderstorms in the four strong convective weather events. In general, the extrapolated forecast using the particle filter blending method was better than that using the cross-correlation or the optical flow method. Since the particle filter blending method can yield better forecasts than the other two methods, it is applicable to the nowcasting of echo motions, which has considerable positive implications for operational weather forecasting.

4 Verification of experimental results

In the previous section, comparative experiments of extrapolated forecasts of four types of convective weather events were conducted by using the particle filter blending, cross-correlation, and optical flow methods. The particle filter blending method performed better than the other two methods in the extrapolated forecast of the powerful squall line steered by westerly winds (the case on 13 April 2016). The echo distribution, shape, intensity, and moving direction were all well predicted. For the extrapolated forecast of the southwesterly monsoon precipitation on 21 May 2016, both the particle filter and optical flow methods performed better than the cross-correlation method. For the local convective precipitation event on 14 April 2016, the 1-h extrapolated forecast produced by the particle filter method was better than that produced by the cross-correlation or the optical flow method. The extrapolated forecast of Typhoon Nida on 2 August 2016 indicated that the particle filter blending method was applicable to the forecasting of a typhoon spiral rainbelt. These individual case studies together suggest that the particle filter method could be applied in the extrapolated nowcasting of echoes.

To quantitatively compare forecasts produced by the three methods, forecast scores were computed according to the forecasts produced by the three methods for the same weather process over the same geographic region and with the same method of data processing. Specifically, the probability of detection (POD), false alarm rate (FAR), and critical success index (CSI) were computed to evaluate the performance of the forecasts. The predicted echoes and the observations at corresponding times were projected to 1 km × 1 km grids and compared. If both the prediction and observation suggested that an echo existed in a grid, the forecast in this grid was regarded as successful; if there was no echo from the prediction but an echo did exist in the observation, the forecast in this grid was regarded as missing; and if an echo was predicted but the echo did not exist in the observation, the forecast in this grid was regarded as a false alarm. A similar approach was applied to the verification of the echo intensity. The POD, FAR, and CSI were calculated as follows:

${\rm{P \!\,\! O \,\!\! D}} = \frac{{{n_s}}}{{{n_s} + {n_m}}},\quad\quad$

(13)

${\rm{FAR}} = \frac{{{n_f}}}{{{n_s} + {n_m}}},\quad\quad$

(14)

${\rm{CSI}} = \frac{{{n_f}}}{{{n_s} + {n_m} + {n_f}}},$

(15)

where n_s, n_f, and n_m are the number of points for a successful forecast, a false alarm, and a miss.

Fourteen cases were selected for the experiments, including westerly-wind precipitation, southwesterly-monsoon precipitation, locally formed precipitation, easterly-wind precipitation, and so on. These cases are listed in Table 1. The 30- and 60-min extrapolated forecasts of the cases were conducted by using the particle filter blending, cross-correlation, and optical flow methods.

Figure 7 presents the averages of POD, FAR, and CSI for the forecasts of westerly-wind precipitation, southwesterly-monsoon precipitation, locally induced precipitation, and easterly-wind precipitation produced by the three methods. The averages of the scores of the five types of precipitation shown in Fig. 7a indicate that the POD of the 30-min forecast produced by the particle filter blending method was 77%, while that produced by the cross-correlation and optical flow methods was 69% and 70%, respectively. The particle filter blending method improved the POD by about 5%, and it also yielded better CSI and FAR scores than the other two methods. The overall performance of the particle filter method was better than those of the cross-correlation and optical flow methods.

For the forecast of precipitation induced by westerly winds (Fig. 7b), the POD of the particle filter method was 82%, which was 5% higher than that of the cross-correlation and optical flow methods; in addition, the CSI was 12% higher and the FAR was lower than that of the other two methods. The POD of the 60-min forecast for the particle filter method was up to 78%. Apparently, the 30- and 60-min forecasts produced by the particle filter method were both better than those produced by the cross-correlation and optical flow methods. For the southwesterly-monsoon precipitation (Fig. 7c), the PODs of the 30- and 60-min forecasts produced by the particle filter method were 75% and 73%, which were 2% higher than those produced by the other two methods. This result suggests that the performances of the three methods were similar. For the forecast of easterly-wind precipitation (Fig. 7d), the POD of the particle filter blending method was higher than that of the other two methods by 4%–5%. For the locally formed slow-moving echo system (Fig. 7e), the 30- and 60-min forecasts produced by the particle filter blending method were better than those produced by the cross-correlation and optical flow methods. The POD for locally formed fast-moving convective precipitation (Fig. 7f) was 66% for the particle filter method and 65% for the other two methods, which was lower than the POD for the forecasts of other types of precipitation by 9%–12%. This result implies that the extrapolation algorithm still has some limitations.

Verification of forecasts for echoes larger than 40 dBZ (Fig. 8) indicates that the PODs of the 30- and 60-min forecasts produced by the particle filter blending method were 48% and 42%, which were higher than the PODs of 42% and 38% produced by the cross-correlation method and the PODs of 43% and 38% produced by the optical flow method. The performance of the particle filter method was better than that of the cross-correlation and optical flow methods. For all three methods, the FAR was high and CSI was low.

In addition, the verification results indicate that, for all three methods, the POD and CSI reduced but the FAR increased with a longer forecast time. This indicates that the performance of the extrapolated forecast became worse with an increase in the forecast time.

5 Conclusions and discussion

This paper introduces the basic principles of the particle filter blending algorithm. A bilateral filter was applied for radar data quality control. Realistic echo motions were approximated by using the particle filter blending method. Semi-Lagrangian extrapolation was applied for the nowcasting of echo motions. The particle filter blending algorithm combines two types of echo motion vectors to obtain the optimum estimation of the echo motion, which is then used to achieve better nowcasting results. Comparative experiments were conducted for five types of precipitation process by using the particle filter blending method, the cross-correlation method, and the optical flow method, separately. Quantitative verification and comparison of the forecasts using the three methods yield the following conclusions.

(1) The bilateral filter method could be applied for radar data quality control and it could effectively remove echo noises. More realistic and smoother radar echoes could be obtained with use of the bilateral filter. Characteristic echo edges were well preserved after the bilateral filtering of the echoes on 13 April 2016. The echo edges were clear, without distinct distortion, suggesting that the filter yielded an ideal output and realistic radar echoes were obtained.

(2) A critical issue for extrapolated forecasting is to obtain realistic and smooth echo motion vectors. In this study, recursive Bayesian filtering was utilized in the particle filter blending algorithm to obtain the optimal track of echo motions. A case study of a squall line process on 13 April 2016 showed that, with the application of the bilateral filter for radar data quality control, the particle filter blending method could yield more realistic and more accurate echo motions.

(3) Extrapolated forecasting was conducted with the semi-Lagrangian method to extrapolate the echo motion vectors. Experimental forecasts of four cases that occurred in the spring of 2016 in Guangdong indicated that the particle filter blending method could yield satisfactory 30- and 60-min forecasts of radar echoes for all four cases. The predicted location, shape, and intensity of echoes were better than those by the cross-correlation and optical flow methods in most cases. Thereby, the particle filter blending method is applicable for the nowcasting of echoes, which has considerable positive implications for operational forecasting.

(4) Evaluation of forecasts of westerly-wind precipitation, southwesterly-monsoon precipitation, easterly-wind precipitation, and locally formed precipitation, produced by the particle filter blending method, cross-correlation method, and optical flow method, indicated that the particle filter method performed better overall than the other two methods. The POD was improved by more than 5%. For precipitation induced by westerly and easterly winds, and locally formed slow-moving precipitation, the forecast produced by the particle filter method was also better than that produced by the other two methods. For the southwesterly-monsoon precipitation, the three methods performed similarly. For locally formed fast-moving echoes, the three methods all failed to successfully forecast the echo motion. The applicability of the three methods decreased with an increase in the forecast lead time.

At present, radar echo extrapolation forecasting is solely based on the single factor of radar reflectivity, without any consideration of storm dynamic and thermodynamic impacts on the radar echo intensity (Zhang et al., 2014). Thereby, extrapolation echo nowcasting has no capability to predict the formation and disappearance of echoes. For thunderstorms that have a life span smaller than the forecast lead time, it makes no sense to conduct extrapolation forecasting (Yu et al., 2012). Extrapolation forecasting only works for those storms with a longer life span, for example, super single-cell storms, powerful squall lines that last for several hours, strong frontal precipitation, and typhoon precipitation. However, predictions of storm generation, intensity change, and disappearance still depend on subjective forecasting.

In the extrapolation forecasting of radar echoes, it is always assumed that the echo motion is steered by environmental flows within the valid forecast time period (1 h) and that the environmental flow varies slowly without any abrupt changes (Zeng et al., 2010). However, environmental flows are often weak during the life cycle of locally formed echoes, which can easily experience abrupt changes during their movement because the steering flow is weak. Thereby, the extrapolation forecasting of such echoes is often unsatisfactory. In addition, the echo movement is affected by both echo advection and propagation. If perturbations caused by echo propagation exist, for example, newly generated single-cell storms in the nearby area, the echo motion will be influenced by echo propagation, which can affect the extrapolated forecast.