The applicability of real-time flood forecasting correction techniques coupled with the Muskingum method

Improving flood forecasting performance is critical for flood management. Real-time flood forecasting correction techniques (e.g., proportional correction (PC) and Kalman filter (KF)) coupled with the Muskingum method improve the forecasting performance but have limitations (e.g., short lead times and inadequate performance, respectively). Here, particle filter (PF) and combination forecasting (CF) are coupled with the Muskingum method and then applied to 10 flood events along the Shaxi River, China. Two indexes (overall consistency and permissible range) are selected to compare the performances of PC, KF, PF and CF for 3 h lead time. The changes in overall consistency for different lead times (1–6 h) are used to evaluate the applicability of PC, KF, PF and CF. The main conclusions are as follows: (1) for 3 h lead time, the two indexes indicate that the PF performance is optimal, followed in order by KF and PC; CF performance is close to PF and better than KF. (2) The performance of PC decreases faster than that of KF and PF with increases in the lead time. PC and PF are applicable for short (1–2 h) and long lead times (3–6 h), respectively. CF is applicable for 1–6 h lead times; however, it has no advantage over PC and PF for short and long lead times, respectively, which may be due to insufficient training and increase in cumulative errors. doi: 10.2166/nh.2019.128 s://iwaponline.com/hr/article-pdf/doi/10.2166/nh.2019.128/639092/nh2019128.pdf Ruixiang Yang Baodeng Hou (corresponding author) Weihua Xiao Xuelei Zhang Baoqi Li State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing, China E-mail: houbaodeng@163.com Ruixiang Yang Chuan Liang Haiying Yu State Key Laboratory of Hydraulics and Mountain River Engineering, College of Water Resource and Hydropower, Sichuan University, Chengdu, China Ruixiang Yang Pearl River Comprehensive Technology and Network Information Centre of Pearl River Water Resources Commission, Ministry of Water Resources, Guangzhou, China Haiying Yu Engineering College, Sichuan Normal University, Chengdu, China


INTRODUCTION
Flood forecasting is a non-engineering measure for flood control development that resulted from the attempt to address the occurrence of flood disasters (Ryder ). Traditional forecasting (TF), based on a physically hydrological model using historical data, reflects only the general regularities of the forecasting area (Kan et al. ). TF errors may be generated and accumulated gradually for a variety of reasons, mainly including model structure errors, model parameter errors and the errors caused by water projects. These errors often increase over time and may eventually exceed the permissible range of precision specified by the standard or used in practice (Liu et al. ). With the development of automatic hydrological monitoring and information transmission technology, real-time flood forecasting is commonly achieved by proportional correction (PC) and Kalman filter (KF) to correct TF using real-time hydrological data (Calvo & Savi ; Liu et al. ).
These two correction techniques have numerous advantages that improve the forecasting performance to some extent; for example, the PC approach boasts a low computational cost, and KF has an unbiased minimum variance (Ocio et al. ). However, some limitations, such as the short lead time of PC and the inadequate correction performance of KF, have been observed in practical applications, and the direct correction of TF values using PC may cause unexpected and unstable deviations (Perumal & Sahoo ).
Particle filter (PF), which has undergone considerable development in recent years, can be applied to any nonlinear and non-Gaussian system represented by a statespace model (Weerts & El Serafy ). Compared with the above correction techniques, PF updates only forecasting value weights and then maintains the actual forecasting values, which avoids the situations in which forecasting values exceed the physical range during the updating process (Peter & Luc ). Therefore, the researches on PF are valuable in hydrology. Since Moradkhani et al. () applied PF to the evaluation of parameters and states in hydrological models, scholars have performed several studies on the same topic. Dechant & Moradkhani () and Bi et al. ()  Combination forecasting (CF) integrates the results of different models to obtain more stable and reliable results (Cloke & Pappenberger ). Many studies have reported that a single model can produce large forecasting uncertainty; accordingly, a multi-model approach has been shown to produce better results (Zsoter et al. 

Study area and flood characteristics
The studied river reach extends from the upstream boundary at the Yongan hydrological station to the downstream boundary at the Shaxian hydrological station, which is located in the middle and lower reaches of the Shaxi River, China (Figure 1). The average annual temperature and precipitation of the basin are 19.9 C and 1,706.5 mm, respectively, and surface water is derived mainly from precipitation. The study area is 2,619 km 2 , and the channel length is 78 km between the two stations. The section widths of the two stations are 186 and 242 m, respectively.
The average flood travel time is 6 h, and there are five runof-river hydropower stations with incomplete daily regulations. Ten typical flood events between the two stations, Based on the analysis of a total of 15 peaks from the 10 studied flood events between the two stations, 6 single-peak flood events, 3 double-peak flood events and 1 significant double-peak flood event with a hidden single-peak are observed. In addition, the relative size of a peak in the same double-peak flood event differs between the two stations.
The peak times for five peaks observed at the Yongan hydrological station and almost all of the rise times occur later than the peak times and rise times observed at the Shaxian hydrological station. After the rainfall forecasting is implemented, the hydropower stations are pre-discharged before the flood to increase their storage capacity, which results in an early flood rise time observed at the Shaxian hydrological station. Additionally, the peak generated by the pre-discharge process occurs earlier than that of a natural flood, causing the peak time at the Shaxian hydrological station to occur earlier than that at the Yongan hydrological station. However, long-duration flood events cannot be completely controlled because of the incomplete daily regulation capacities at the five run-of-river hydropower stations. Consequently, the flood events observed at the Shaxian hydrological station still largely retain natural flood characteristics.
The flood forecasting error distribution during the calibration period (before 2002) using TF is analysed with the prior probability density function of PF. The mean error value is À156.20, and the standard deviation is 384.46. The error distribution is considered Gaussian, but it is rejected at the 0.05 significance level. Therefore, this distribution is considered to be noncompliant with the Gaussian distribution, which is similar to the results of Zhang et al. ().
Considering the t location-scale distribution (representing the translation and expansion of the standard t distribution) as the error distribution, the probability density function is calculated by the following equation: where u represents the location parameter; σ represents the scale parameter and v represents the degree of freedom. The maximum likelihood estimation results of the three parameters are u ¼ À 130, σ ¼ 287 and v ¼ 4, which are accepted at the 0.05 significance level. Therefore, the error distribution is considered to fit the t location-scale distribution. Figure 2 shows the error frequency histograms and the corresponding Gaussian and t location-scale fits.
Muskingum method and real-time correction techniques A flood routing model using the Muskingum method is applied to TF in the study area (Zhao ). The basic principle of the Muskingum method is to divide a river reach into n unit reaches. After n À 1 iterations, the discharge of unit reach n is obtained. The Muskingum method is calculated as follows: where i represents the unit reach i, i ¼ 1, …, n À 1; t represents the time; Δt represents the time interval; c represents the forecasting value; Q represents the discharge (e.g., Q iþ1 c,tþΔt represents the forecasting discharge at time t þ Δt for unit reach i þ 1); and C 1 , C 2 and C 3 represent the model coefficients, which are estimated as follows: where x represents the weighting factor of the discharge; and k represents the storage constant, reflecting the flood travel time for the storage volume of a river course in a steady flow state. The TF discharge is obtained by adding the local inflow, which is derived from the related rainfall forecasting and rainfall-runoff model in practice, to the discharge of unit reach n.
PC is a real-time correction technology; this technique is used for practical forecasting in the study area. The main assumption is that the ratio of the observed value to the forecasting value is relatively stable over a short period of time.
The ratio at the current time is equal to that at a lead time before the current time. Multiplying the forecasting value by the ratio achieves the correction value, and the correction value for unit reach n is as follows: where u represents the correction value; o represents the observed value and t 0 represents the lead time.
KF is a modern analytical technology for dynamic systems that mainly includes multivariable controls, optimal controls, and estimation and adaptive controls. (5) estimate the filtering error covariance; (6) calculate the filtered value and (7) repeat steps (2)-(6) until the steps have been operated for all moments. Steps (3)-(5) are expressed in the following equations: where P represents the error covariance matrix; Φ represents the state transition matrix; Γ represents the state noise distribution matrix; O represents the state noise covariance matrix; K represents the gain matrix; H represents the observation matrix; R represents the observation noise covariance matrix; and I represents the identity matrix. The correction value for unit reach n is as follows: PF is an algorithm that uses the Monte Carlo algorithm to implement Bayesian estimation theory, which can be applied to any type of state-space model (Han & Li ).
The main idea is to select a set of weighted and random sample particles from the state space to achieve an approximation of the state probability density distribution and then replace the integral operation with the sample mean to obtain the minimum variance estimate of the state. The steps are as follows (van Leeuwen ): (1) particle initialization. This process considers the number of particles, which is set to N, and the weight is equal to 1/N. The perturbation term is added for this step.
(2) Particle weight updating. This process involves assuming that N particles can approximate the importance probability density function of the state values that is commonly expressed in practice with a prior probability density function. After obtaining the observed values, the weight of each particle is recalculated according to the difference between the particle and the observed value, where particles with small differences are assigned large weights. Then, the particle weights are normalized.
(3) Resampling (i.e., copying the particles according to the weights). This process is commonly performed when the number of effective particles is less than the set threshold. The greater the weight is, the larger the number of copies. After resampling, the weight of each particle is equally set to 1/N. (4) Calculate the filtered value. The state transfer equation of the system is used to forecast the state of each particle to obtain a forecasting value. Based on the forecasting value, the weight of the particle is updated using the same methods in step (2). (5) Steps (2)-(4) should be repeated until the steps have been performed for all moments. The correction value for unit reach n is as follows: where w represents the normalized weight of each particle and j represents a particle.
The multiple linear regression is used to construct where a m represents the weight of correction technique m and M represents the number of correction techniques employed.

Performance indexes
In this paper, the flood forecasting performance is evaluated with two indexes: the overall consistency and the permiss- The NSE is calculated by the following equation: where T represents the total time interval and frequencies not only increase the computational costs but also reduce the particle diversity (Zuo ). Accordingly, one flood event is selected here to analyse the computing time of PF. The number of particles is set to 100, 1,000 and 10,000, and the required computation times are 0.9, 2.8 and 19.2 s, respectively. However, the change in the NSE is negligible, and the largest number of effective particles is greater than two-thirds of the total. Therefore, the number of particles is set to 1,000, and the resampling threshold is two-third of the number of particles.

RESULTS AND DISCUSSION
Correction performance for 3 h lead time  Table 1.
For the NSE index of the overall consistency, the NSE of PC is larger than that of TF (except for the second flood For the peak discharge index of the permissible range, the qualified numbers of flood events using TF, PC, KF and PF are 12, 10, 13 and 14, respectively; the numbers of peaks with the absolute error percentages below 5% are 7, 4, 9 and 7, respectively; the numbers of peaks with percentages between 5 and 10% are 2, 3, 4 and 6, respectively. Therefore, the performances of KF and PF are similar in terms of the peak discharge index; KF performs slightly better, whereas PC performs even worse than TF. For the peak time index of the permissible range, the qualified numbers of flood events using TF, PC, KF and PF are 8, 12, 13 and 14, respectively, and the numbers of flood events with the minimum absolute error values are 4 (2 parallel), 5 (3 parallel), 6 (3 parallel) and 9 (6 parallel), respectively.
Therefore, in terms of the peak time index, the performance of PF is optimal; none of the remaining correction techniques shows obvious advantages, but they perform better than TF.
For the process discharge index of the permissible range, the numbers of flood events with qualified rates in the 75-90% range are 1, 2, 3 and 3 for TF, PC, KF and PF, respectively, and the numbers of flood events with qualified rates in the 60-75% range are 2, 6, 4 and 3, respectively, while the numbers of flood events with the highest qualified rates are 1 (1 parallel), 4 (1 parallel), 3 (1 parallel) and 4 (1 parallel), respectively. Therefore, for the process discharge index, the performances of PC and PF are better than that of KF. In conclusion, for the permissible range index, the order of performance is the same as that for the overall consistency.
The discharge processes and performances of TF, PC, KF and PF for the 10 flood events are shown in Figure 3.  for the ninth to tenth flood events is shown in Table 2.  however, the absolute error is only 1 h. Both flood events account for more than 70% of the process discharge index of the permissible range. In total, although there is no direct comparison between CF and KF in Table 2, the performance of CF is clearly better than KF.

Correction performance for different lead times
For different lead times (1-6 h), the overall consistency index values of the three correction techniques and CF are shown in Table 3. However, comparing these three techniques individually is difficult because there are 10 flood events and 6 lead times. Therefore, the following three indexes are generalized: the slope of the average decrease in NSE, the number   is multiplied by the error at time t and added to the forecasting value at time t þ t 0 to obtain the correction value at time t þ t 0 . For PF, the particle weight corresponding to the forecasting value is calculated according to the error at time t.
The error probability is obtained from the importance probability density function (i.e., the TF error distribution in the section 'Study area and flood characteristics'). In addition, the error probability is multiplied by the initial weight of 1/N and normalized to obtain the updated weight; and the correction value at time t þ t 0 is obtained by multiplying the forecasting value at time t þ t 0 by the updated weight and summing the results. Therefore, the correction performances of KF and PF are less affected by the lead times than that of PC.
One of the reasons why CF does not perform best for The main conclusions are as follows: 1. For 3 h lead time, the performance of PF in the overall consistency index, including the NSE value and the precision grade, is optimal, followed by those of KF and PC. For the permissible range index, the peak discharge performances using KF and PF are similar and those using PF are slightly superior; however, PC performs the worst among these techniques. The peak time performance using PF is better than that using the other two techniques, which exhibit similar performances.
The process discharge performance of PC is slightly better than that of PF and much better than that of KF.
Furthermore, the weighted average CF including the three correction techniques is constructed. The CF weights decrease in the order of PF, PC and KF, and its performance is similar to PF. Generally, the comprehensive performances of all the considered techniques (from adequate to poor) are in the order of PF, CF, KF and PC.
2. When the lead time increases from 1 to 6 h, the decline in the performance of PC is much faster than those of KF and PF, indicating that the influence of the lead time on PC is much larger than those on the other two techniques. The performance of PC is optimal for 1 and 2 h lead times, while the performance of PF is optimal for 3-6 h lead times. For 1 and 2 h lead times, the NSE values of CF are higher than those of PC, KF and PF; and the weights of PC are the largest and far greater than those of KF and PF. For 3-6 h lead times, the NSE values of CF are also the top 1 (but parallel); PC, KF, PF and CF are applicable for 1-3, 1-5, 1-6 and 1-6 h lead times, respectively. For short lead times (1 and 2 h), PC is applicable; its performance is the best except for CF, and the computation requires minimal time for accurate forecasting. For long lead times (3-6 h), PF is applicable because of its highest performance. In contrast, CF, with its long computation time, should not be selected for short and long lead times. PC and KF lead times are increased by 3 and 1 h, respectively, by using PF. Furthermore, PF coupled with the Muskingum method can provide high-precision forecasting results for longer lead times to support flood control strategies.
The importance of flood peak discharge and peak time forecasting is self-evident, but their performance is not satisfactory. Therefore, improving the performance of real-time correction techniques for peak discharge and peak time forecasting should be the focus of future studies. Comparisons among extended KF, Ensemble KF and PF applied to nonlinear systems should be performed, and the application of PF in channels with longer flood travel times should also be further studied. In addition, the impacts of the methods on the uncertainty estimates and the uncertainties arising from the quality and amount of observed data should be considered. A greater number of flood events are needed to study the performance of CF and multimethod CF (i.e., different indexes should be evaluated by a combination of results from different forecasting methods), which represent two methods of forecasting using different techniques.