“ In-Process Type ” Dynamic Muskingum Model Parameter Estimation Method

This paper discusses the Muskingum model as a novel parameter estimation method. Sixty representative floods over the past four decades serve as research objects; a linear Muskingum model and Pigeon-inspired optimization (PIO) algorithm are used to obtain the parameters of each flood. The proposed “in-process type” dynamic parameter estimation (IP-DPE) method is used to establish the characteristic attributes set of 50 floods. The characteristic attributes set refers to a set of parameters that could describe the shape, magnitude, and duration of the flood before flood peak; they are the input, whereas parameters K and x of each flood are the output to establish a Neural Network model. Then we input flood characteristic attributes to obtain flood parameters when estimating flood parameters practically. Ten floods were used to test the parameter estimation and flood routing efficacy. The results show that the IP-DPE method can quickly identify parameters and facilitate accurate river flood forecasting.


Introduction
The Muskingum model has been commonly utilized for river flood forecasting across various parts of the world since it was first proposed in the 1930s [1,2].For the model to work properly, its parameters must be accurately estimated.Several methods for such parameter estimation have been proposed in recent years, including trial-and-error [3], the least-square method [4][5][6], and nonlinear programming [7,8].In practice, however, these approaches come with high computational complexity, poor universality, and susceptibility to local optima [9].Advancements in computer technology have allowed for innovative parameter estimation optimization algorithms such as the Genetic Algorithm [10,11], Particle Swarm Optimization Algorithm [12,13], Chaotic Particle Swarm Optimization Algorithm [14], and Ant Colony Algorithm [15].
Both traditional approaches and optimized algorithms have assisted in developing a useful Muskingum model parameter estimation method [2].However, there still are important problems to be resolved: (1) At present, most methods are developed with the sole focus on improving the simulation accuracy of one flood [16][17][18], such as in [16], where the average relative errors obtained from four algorithms were 1.16%, 1.11%, 1.13%, and 1.10%, respectively, and the improvement was inherently limited.Also, although the obtained flood parameters could forecast the selected

Introduction to Muskingum Model Theories
The flood routing of a river course is meant to deduce the flood process of the downstream section according to the process of the upstream section.The Muskingum model is an effective channel flood routing model that has been widely used.This method is mainly used to carry out flood routing of a river course by establishing a river channel storage equation and solving the water balance equation simultaneously.The basic equation for flood routing in the Muskingum model is: in which W is the river course storage volume, t is the time, I and Q are the inflow and outflow of the river reach, respectively, x is the weighting factor, and K is a constant of the Muskingum model for the river reach.The differential Equation (1) can be discretized as follows: Water 2017, 9, 849 where Q(i) and Q(i) are the calculated outflow and the observed outflow, respectively, at the calculation time i, I(i) is the inflow at the calculation time i, n is the total amount of time, c 0 , c 1 and c 2 are flow routing coefficients, and c i = f (K, x, ∆t), i = 0, 1, 2 satisfy the following: The flood flow routing equation of the Muskingum model consists of Equations ( 2) and (3).As discussed above, parameter estimation is one of the most important aspects of correctly applying the Muskingum model.It involves the estimation of parameters K and x (or parameters c 0 , c 1 , c 2 ).

Establishing the Muskingum Parameter Optimization Model
According to Equations ( 2) and (3), the parameter optimization problem of the Muskingum model can be transformed as follows: In the Equation ( 5), the objective function minF refers to minimizing the relative error sum between calculated outflow and observed outflow.In the constraints, the range of c 0 , c 1 and c 2 = 1 − c 0 − c 1 are all from 0 to 1 [9,16,23], and, considering that if the value of c 2 is less than 0, it will not satisfy the constraints c 2 ∈ [0, 1], we added the penalty function h(g 3 (c 0 , c 1 )) to ensure that c 2 meet the constraints of g 3 .
The optimization function is: where h(g 3 (c 0 , c 1 )) is the penalty term.When the constraint g 3 is satisfied, the value is 0; otherwise, the value is 10 6 .Equation ( 6) is the objective function of the Muskingum parameter optimization model based on the PIO algorithm.

Single-Flood Muskingum Model Parameter Estimation Based on Pigeon-Inspired Optimization Algorithm
Reliable input and output samples are required for dynamic Muskingum model parameter estimation.The output samples are the parameter estimation values of each flood.We used the Pigeon-inspired optimization (PIO) algorithm to estimate the each flood parameter.

Pigeon-Inspired Optimization (PIO) Algorithm
Pigeons are naturally excellent navigators.Inspired by pigeons' unique homing characteristics, Duan et al. [24] first proposed the PIO intelligent swarm optimization algorithm.Guilford et al. [25] found that pigeons use different navigation tools at different stages of their journey, relying on the Earth's magnetic field in the early stages and later relying on landmarks.Whiten et al. [26] proved that the sun is also an important factor that influences pigeons' navigation ability.Pigeons initially sense magnetic fields then draw a mental map accordingly, towards which they constantly adjust their flight direction.As they approach their destination, they use familiar landmarks as their primary navigation tool (or follow the others who are familiar with the landmarks).To simulate these two stages of homing behavior, the PIO algorithm uses the Map and Compass operator and the Landmark operator.While solving the optimization problem, the Map and Compass operator is first used to conduct calculations until the maximum set number of iterations, then the currently obtained positions of all individuals in this stage are passed to the Landmark operator to continue the calculation in the second stage.A detailed introduction to the PIO algorithm can be found in our references, (e.g., [27]).
The PIO algorithm solution procedure for the Muskingum model parameter estimation problem is as follows.
Step 1 Generate the initial population, in which the position of each individual i is X i = {x i1 , x i2 , . . . ,x iD }, which means the values of c 0 , c 1 , c 2 , the initial velocity and the position X i satisfy the upper and lower bound (i.e., range = [0, 1]) of the problem to be solved.That is to say, the value range of c 0 , c 1 , c 2 is from 0 to 1.The maximum iteration of the Map and Compass operator is MC max , the maximum iteration of the Landmark operator is L max , the population size is N, the dimension of the problem is D, and the compass factor is R.

Step 2
Set the best position of the individuals X p = X i , in which X p means the optimal value of c 0 , c 1 , c 2 , and iterations k = 1, then calculate the fitness value of all individuals of the population (i.e., the solution of Equation ( 6)) and obtain the optimum value X best , which is the optimal value of the optimization function.

Step 3
Start the iteration procedure of the Map and Compass operator.Update the position X i and velocity V i of each individual in each iteration procedure, and note that each X i should satisfy the range.Run this iteration until the iterations k = MC max , then update X p and X best .

Step 4
Take the population obtained in Step 3 as the input of the iteration procedure of the Landmark operator.In each iteration, the individual population is first ranked by the fitness value; the first half is selected as the reference population, then the center position of the better population X center is calculated as the landmark.Next, update the position of each individual until the iterations reach the maximum MC max , then stop the iteration procedure and update the best position X p and best fitness X best .
Step 5 X best is the required global optimum.

Single-Flood Muskingum Model Parameter Estimation
The river channel from Chenggouwan to Linqing is located at the South Grand Canal of the Haihe River Basin.It is about 83.8 km long and contains no branch import.We selected a flood event which occurred in this river channel in 1961 as the flood routing example and applied the PIO algorithm, Least-Square (L-S) method [28], Accelerated Genetic Algorithm (AGA) [29], and Ant Colony Optimization (ACO) algorithm [17] to solve a single-flood Muskingum model parameter estimation problem.We then evaluated the estimation performance and flood routing accuracy of all four algorithms by comparison.Our PC infrastructure includes an Intel Core i7-3537U CPU with a clock speed of 2.50 GHz, a Windows 8.1 operating system, and the Java programming language.The PIO algorithm parameter set included population size N = 100, dimension D = 3, compass factor R = 0.5, maximum iterations of the Map and Compass operator MC max = 500, maximum iterations of the Landmark operator L max = 200, and upper and lower bounds of the search space range = (0, 1).The results are provided below.
The absolute error of the simulated results obtained from the four algorithms is shown in Figure 1.

Absolute error was calculated as follows:
Absolute error = |simulated value − observed value| (7) The relative error of the simulated results obtained from the four algorithms is shown in Figure 2, as calculated by: relative error = simulated value − observed value observed value The simulation results of four algorithms are shown in Figure 3.As shown in Figure 3, the PIO algorithm performs well in terms of solving the parameter estimation problem of the Muskingum model.The parameter estimation results of the four algorithms, their average absolute error, and the average relative error among their simulation results are shown in Table 1.
The simulation results of four algorithms are shown in Figure 3.As shown in Figure 3, the PIO algorithm performs well in terms of solving the parameter estimation problem of the Muskingum model.The parameter estimation results of the four algorithms, their average absolute error, and the average relative error among their simulation results are shown in Table 1.The simulation results of four algorithms are shown in Figure 3.As shown in Figure 3, the PIO algorithm performs well in terms of solving the parameter estimation problem of the Muskingum model.The parameter estimation results of the four algorithms, their average absolute error, and the average relative error among their simulation results are shown in Table 1.The simulation results of four algorithms are shown in Figure 3.As shown in Figure 3, the PIO algorithm performs well in terms of solving the parameter estimation problem of the Muskingum model.The parameter estimation results of the four algorithms, their average absolute error, and the average relative error among their simulation results are shown in Table 1.The PIO algorithm results have the lowest average absolute error and average relative error, i.e., PIO possesses better solution accuracy and better optimization performance than the three other algorithms we tested.Also, its computational time was 1.625 s.Thus, the PIO algorithm is an excellent approach to solving the Muskingum model parameter estimation problem.The artificial neural network is a simulation of human physiological structures and mechanisms which depart from symbolic reasoning or logical thinking [30,31].It plays a role in mitigating imperfect associative memory, defective feature pattern recognition, automatic rule-learning, and other functions.It plays an especially crucial role in solving complex problems [32][33][34] which are difficult to model, such as the computer-aided molecular design problem [32], spatial distribution of soil heavy metals problem [33], and modulation transfer function estimation problem [34].
The Back-Propagation (BP) Neural Network is one of the most commonly used neural network models.The research group led by Rumelhart and McClelland [35] analyzed the error back-propagation algorithm of multilayer feed-forward networks, which possess a nonlinear continuous transfer function, in the book Parallel Distributed Processing.The basic idea of the BP algorithm is that the learning process consists of two processes: forward propagation of signal and backward propagation of error.In the process of forward propagation, the input samples are imported from the input layer and then transmitted to the output layer after being processed in the hidden layer.If the errors between actual outputs of the output layer and the expected outputs do not meet the requirements, they then transfer to the error backward propagation.The error back-propagation is used to transfer the output errors to the input layer through the hidden layer, and then apportion the error to all the units of each layer in order to obtain the error signal of each unit, which is the basis for correcting each unit.The process of adjusting the weight of each layer in signal forward-propagation and error back-propagation goes on again and again, and the procedure of weight adjustment is the so-called learning and training process of the network.This process will not stop until the error of the network output is reduced to an acceptable level or to the preset number of studies.
Its advantages include: (1) a simple structure and strong operability; (2) nonlinear mapping from input to output, which realizes any complex nonlinear mapping; and (3) strong self-learning and summarization abilities [36].We selected a three-layer BP network for dynamic Muskingum model parameter estimation.The three layer BP network topology is shown in Figure 4.
The network contains an input layer, a hidden layer, and an output layer, and the learning procedure of the BP-ANN model is shown as below.
Step 1 Network initialization.The number of node in the input layer is n, in the hidden layer is p, and in the output layer the node number is m.The connection weights between the input layer and the hidden layer, the hidden layer and the output layer is W j i and V k j , respectively.The threshold of hidden layer is θ j (i = 1, 2, . . ., p) and the threshold of output layer is Step 2 Calculate the output h j of hidden layer.The calculation formula is: where f is the excitation function of hidden layer, and X i is the input of node i.

Step 3
Calculate the output Y k of output layer.It is calculated by Step 4 Update the weights through the equations as follows: where η is the learning rate and η > 0, D(t) = −∂J/∂W k j (t), D (t) = −∂J/∂V k j (t), and β is momentum factor and 0 ≤ β < 1.
Step 5 Threshold update.Update the value of θ j and µ k according to the error between the network output Y k and the desired output y k .The formulas to calculate θ j and µ k are: Step 6 Determine whether the iteration ends.If it doesn't end, go back to step 2.
Water 2017, 9, 849 7 of 19 where f is the excitation function of hidden layer, and i X is the input of node i.
Step 3 Calculate the output k Y of output layer.It is calculated by Step 4 Update the weights through the equations as follows: where  is the learning rate and , and  is momentum factor and 0 1    .
Step 5 Threshold update.Update the value of j  and k  according to the error between the network output k Y and the desired output k y .The formulas to calculate j  and k  are: Step 6 Determine whether the iteration ends.If it doesn't end, go back to step 2.

Basic Theory of "In-Process Type" Dynamic Parameter Estimation
As discussed above, the Muskingum model has continually increased in popularity since its proposal in 1938; the most important function of the model is the selection of parameters K and x.Hydrologists such as Xu et al. [16,18,37] have taken a single set of parameters as the Muskingum model parameters in one river channel, which can be determined by the estimation of one of the floods or the mean estimated value of multiple floods.Modern optimization algorithms, as opposed to traditional methods like trial-and error, have contributed substantially to the development of today's single-flood parameter estimation models.However, a river with only a single set of model parameters does not readily satisfy the application requirements [9].Thus, Liang [38] proposed that the parameter K should be determined by the flow classification method-that is to say, the peak discharge with varying magnitudes should adopt different K values while the parameter x is still determined by the estimated mean.Although this remedies the parameter estimation models.However, a river with only a single set of model parameters does not readily satisfy the application requirements [9].Thus, Liang [38] proposed that the parameter K should be determined by the flow classification method-that is to say, the peak discharge with varying magnitudes should adopt different K values while the parameter x is still determined by the estimated mean.Although this remedies the application defect of a single parameter set to some extent, the problem remains fundamentally unresolved due to the large nonlinearity between the peak discharge and the values of K and x.We suggest that the flood model parameters should be determined by the specific flood condition, and that real-time dynamic parameter estimation of the Muskingum model for river course flood should also be realized.
Most of the commonly used Muskingum model parameter estimation methods are "posterior type" [9], where flood parameters are estimated via simulation of floods which have already occurred.Of course, the accuracy of this type of flood simulation is generally high, but the actual forecast may have insufficient accuracy even when the parameter accuracy of the posterior estimation is sufficient.Each flood has its own characteristics, which causes variations in the model parameters (or sometimes even markedly diverge them).The "in-process type" parameter estimation method better accounts for this; it estimates parameters while a flood is occurring, then assigns a set of unique parameter sets to the flood.
As to the real-time flood forecasting, it corrects the results and parameters in the process of model prediction.The research of real-time flood forecasting was carried out earlier, and Japanese scholar M. Hino [39,40] proposed that the method of Kalman Filter could be applied to hydrological forecasting, and he developed it in his later papers in 1973.And Todini [41], a professor of University of Bologna in Italy, established the mutually interactive state-parameter (MISP) algorithm in 1978, becoming the first-generation real-time flood forecasting model.Kleanidis and Bras [42] improved the Sacramento model in 1980, realizing real-time prediction through the application of the Kalman filtering model.In 1987, Georgakakos [43] presented a stochastic and dynamic meteorological and hydrological forecast model.
In 1995, Song [44] first introduced and improved the method of recursive least-squares parameter identification raised by T.R. Fortesuce et al., based on variable forgetting factors, which provided a powerful tool for parameter identification of real-time flood forecasting model.Then, in 1995, Zhong et al. [45] tried to introduce artificial neural network technology into real-time flood forecasting.In 1996, Guo et al. [46] developed an integrated automation system of reservoir operation for the Tuo Xi hydropower station in Hunan.From then on, a decision support system for real-time flood forecasting, flood control, and power generation dispatching had been developed, and the automation of forecast dispatching had also been realized.In 2011, Luo et al. [9] thought that floods from the same river should have different parameters, and presented the Muskingum model with adaptive parameter using the time series similarity search method.However, this article also has some shortcomings; only after the flood from upstream has finished can the parameters be obtained, thus the forecast period of the flood is shortened.From this point of view, this paper puts forward the research idea.
In this study, based on the previous single-flood parameter estimation method in Section 2.2.2, we first estimated the flood parameter and extracted the characteristic attributes of several floods, then took these attributes as inputs and flood parameters as the outputs to train the BP-ANN model.After that, if we could obtain the attributes of the happening flood, then the flood parameters could also be obtained from the BP-ANN model.Thus, the Muskingum model parameters could be dynamically estimated in real-time.The proposed method was implemented as follows.
Step 1 Select multiple floods and use the PIO algorithm to estimate their respective parameters.
Step 2 Extract the characteristic attributes of each flood according to the specific river conditions.Step 3 Take the extracted characteristic attributes of each flood as the input, and parameters K and x of the Muskingum model as outputs to train the Neural Network model; use the PIO algorithm to estimate the parameters.
Water 2017, 9, 849 9 of 19 Step 4 Extract the characteristic attributes of the flood to be forecasted and plug them into the Neural Network to calculate parameters K and x of the Muskingum model, then use the Muskingum model to perform the flood routing computation.
The Neural Network model is not invariable once it has been established.New data should be updated into the database and the model continually trained to improve its prediction accuracy.

Study River Channel
The Yellow River is located in Northern China.It has a total length of 5464 km and its basin covers 752,443 km 2 , making it the fifth longest river in the world and the second longest river in China.The annual precipitation in most parts of the Yellow River Basin is about 200 to 650 mm, and the annual precipitation in the lower reach is more than 650 mm.The Sunkou and Gaocun Hydrometric Stations, which are the national key hydrometric stations in the lower reaches of the Yellow River, are located in Shandong Province.There are about 197 km between the two hydrometric stations, in which there are no tributary imports or exports.We selected this river channel to validate the proposed dynamic parameter estimation method of the Muskingum model (Figure 5).The Yellow River is located in Northern China.It has a total length of 5464 km and its basin covers 752,443 km 2 , making it the fifth longest river in the world and the second longest river in China.The annual precipitation in most parts of the Yellow River Basin is about 200 to 650 mm, and the annual precipitation in the lower reach is more than 650 mm.The Sunkou and Gaocun Hydrometric Stations, which are the national key hydrometric stations in the lower reaches of the Yellow River, are located in Shandong Province.There are about 197 kilometers between the two hydrometric stations, in which there are no tributary imports or exports.We selected this river channel to validate the proposed dynamic parameter estimation method of the Muskingum model (Figure 5).

Extracting Flood Characteristic Attributes
There are several characteristic attributes that influence flood routing processes, including antecedent rainfall, soil moisture content, water level, flow velocity, and peak discharge.We analyzed the 60 floods which occurred over the past 40 years according to the Gaocun and Sunkou Hydrometric Stations, and we found that it takes about 6-65 h for the flood routing from Gaocun to Sunkou; a single flood in this portion of the Yellow River lasts for about 5-36 days.After a comprehensive analysis, we extracted the following characteristic attributes of the floods in this river reach via the proposed "in-process type" estimation method, which can be shown in Figure 6:

Extracting Flood Characteristic Attributes
There are several characteristic attributes that influence flood routing processes, including antecedent rainfall, soil moisture content, water level, flow velocity, and peak discharge.We analyzed the 60 floods which occurred over the past 40 years according to the Gaocun and Sunkou Hydrometric Stations, and we found that it takes about 6-65 h for the flood routing from Gaocun to Sunkou; a single flood in this portion of the Yellow River lasts for about 5-36 days.After a comprehensive analysis, we extracted the following characteristic attributes of the floods in this river reach via the proposed "in-process type" estimation method, which can be shown in Figure 6: (1) The initial water level (IWL), which corresponds to the rising point (the time at t 1 ) of the flood, reflects the rising conditions of the flood and the previous water volume of the river.The larger the previous water volume, the more prone the area of interest is to flooding.The water level also reflects flow velocity information and thus exerts some influence on flood routing.(2) The peak discharge (PD), which reflects the magnitude of the flood and has significant influence on flood routing.(3) The flood peak emergence time (FPET), i.e., the time at which the flood appears.FPET is usually an instantaneous value, so it cannot exactly reflect the characteristics of each flood along the entire route.Here, we consider the point at which the flood peak appears as the FPET; that is to say, the FPET = 1, 2, . . ., T, in which T is the duration time of the flood.
(4) The total flood volume on the first two days after a flood appears (FVFT) reflects the rising conditions of the flood.( 5) The flood volume on the day before the flood peak appears (FVBFP), which reflects the width of the flood peak.(6) The average discharge before the flood peak (ADBFP), i.e., the mean discharge prior to the flood peak, which reflects the ratio between the total flood volume and time.

Extracting Flood Characteristic Attributes
There are several characteristic attributes that influence flood routing processes, including antecedent rainfall, soil moisture content, water level, flow velocity, and peak discharge.We analyzed the 60 floods which occurred over the past 40 years according to the Gaocun and Sunkou Hydrometric Stations, and we found that it takes about 6-65 h for the flood routing from Gaocun to Sunkou; a single flood in this portion of the Yellow River lasts for about 5-36 days.After a comprehensive analysis, we extracted the following characteristic attributes of the floods in this river reach via the proposed "in-process type" estimation method, which can be shown in Figure 6:   Based on the first 50 floods in our dataset, we used the PIO algorithm to estimate the Muskingum model parameters.We used the same PIO parameters as those defined in Section 2.2.2 and calculation period dt = 12 h, then as to each flood, we ran the algorithm five times to obtain five results of the parameter to make sure that the results are believable.We took the optimum of the objective function as the parameter estimation result.The last ten floods were used to validate the constructed dynamic parameter estimation model, while PIO was used to estimate parameters K and x.We can obtain the correlation between the characteristic attributes and K and x, as shown in Figure 7, according to the pre-extracted characteristic attributes of each flood.
Water 2017, 9, 849 10 of 19 (1) The initial water level (IWL), which corresponds to the rising point (the time at 1 t ) of the flood, reflects the rising conditions of the flood and the previous water volume of the river.The larger the previous water volume, the more prone the area of interest is to flooding.The water level also reflects flow velocity information and thus exerts some influence on flood routing.(2) The peak discharge (PD), which reflects the magnitude of the flood and has significant influence on flood routing.(3) The flood peak emergence time (FPET), i.e., the time at which the flood appears.FPET is usually an instantaneous value, so it cannot exactly reflect the characteristics of each flood along the entire route.Here, we consider the point at which the flood peak appears as the FPET; that is to say, , in which T is the duration time of the flood.
(4) The total flood volume on the first two days after a flood appears (FVFT) reflects the rising conditions of the flood.( 5) The flood volume on the day before the flood peak appears (FVBFP), which reflects the width of the flood peak.( 6) The average discharge before the flood peak (ADBFP), i.e., the mean discharge prior to the flood peak, which reflects the ratio between the total flood volume and time.
Based on the first 50 floods in our dataset, we used the PIO algorithm to estimate the Muskingum model parameters.We used the same PIO parameters as those defined in Section 2.2.2 and calculation period , then as to each flood, we ran the algorithm five times to obtain five results of the parameter to make sure that the results are believable.We took the optimum of the objective function as the parameter estimation result.The last ten floods were used to validate the constructed dynamic parameter estimation model, while PIO was used to estimate parameters K and x.We can obtain the correlation between the characteristic attributes and K and x, as shown in Figure 7, according to the pre-extracted characteristic attributes of each flood.Figure 7a,c,e,g,i,k indicate that the value of parameter K varies from 5 to 65.K decreases as the value of each characteristic attribute increases; it has a large variation range.These results suggest that error is inevitable regardless of whether the mean K value or the K value in accordance with the peak discharge magnitude is utilized.Figure 7a,c,e,g,i,k indicate that the value of parameter K varies from 5 to 65.K decreases as the value of each characteristic attribute increases; it has a large variation range.These results suggest that error is inevitable regardless of whether the mean K value or the K value in accordance with the peak discharge magnitude is utilized.
As shown in Figure 7b,d,f,h,j,l, the value of parameter x fluctuates between 0 and 0.5; in other words, it fluctuates considerably as each characteristic attribute increases in value.We found no obvious change trend in x.To this effect, substituting the mean value of x into the Muskingum model will cause an error.

The Flow Chart of the IP-DPE Method
To clearly explain the process of the proposed IP-DPE method, we provide a flow chart as Figure 8.
As shown in Figure 7b,d,f,h,j,l, the value of parameter x fluctuates between 0 and 0.5; in other words, it fluctuates considerably as each characteristic attribute increases in value.We found no obvious change trend in x.To this effect, substituting the mean value of x into the Muskingum model will cause an error.

The Flow Chart of the IP-DPE Method
To clearly explain the process of the proposed IP-DPE method, we provide a flow chart as Figure 8.

Results
We optima [47,48], we used the PIO algorithm to realize parameter estimation in the Neural Network model, using the same parameters referred to above.The structure diagram of the resulting BP-Neural Network is shown in Figure 8.
The first 50 floods were substituted into the Neural Network to train the model, and the last 10 floods were used to verify the model.We also ran the traditional Muskingum model parameter optimization methods, mean value (MV) method, and flow classification (FC) method for comparison against the proposed "in-process type" dynamic parameter estimation (IP-DPE) method.According to the Standard for Hydrological Information and Hydrological Forecasting (SL 250-2000), we used the following three indexes to evaluate the results.
(1) The mean absolute error (MAE) of the flood (m 3 /s): where T is the duration of the flood, Y oi is the observed flow of the flood at time i, and Y ci is the calculated flow of the flood at time i. (2) The average relative error (ARE) of the flood (%): (3) The peak flow relative error (PFRE) of the flood (%): where P o represents the observed peak flow of the flood and P c is the calculated peak flow of the flood.The evaluation results are shown in Table 2.The simulated performance of the three methods is shown in Figure 9.
Water 2017, 9, 849 13 of 19 comprises the "6-5-2 type" BP-Neural Network model.There are 42 parameters in the Neural Network model, as established.Because the traditional Newton Iteration and Gradient methods easily fall into local optima [47,48], we used the PIO algorithm to realize parameter estimation in the Neural Network model, using the same parameters referred to above.The structure diagram of the resulting BP-Neural Network is shown in Figure 8.The first 50 floods were substituted into the Neural Network to train the model, and the last 10 floods were used to verify the model.We also ran the traditional Muskingum model parameter optimization methods, mean value (MV) method, and flow classification (FC) method for comparison against the proposed "in-process type" dynamic parameter estimation (IP-DPE) method.According to the Standard for Hydrological Information and Hydrological Forecasting (SL 250-2000), we used the following three indexes to evaluate the results.
(1) The mean absolute error (MAE) of the flood (m 3 /s): where T is the duration of the flood, oi Y is the observed flow of the flood at time i, and ci Y is the calculated flow of the flood at time i.
(2) The average relative error (ARE) of the flood (%): (3) The peak flow relative error (PFRE) of the flood (%): where o P represents the observed peak flow of the flood and c P is the calculated peak flow of the flood.The evaluation results are shown in Table 2.The simulated performance of the three methods is shown in Figure 9.As shown in Table 2 and Figure 9, the flood forecasting results obtained from the three methods are all technically reasonable but do show marked differences in forecasting accuracy.The MAE value can well reflect the actual prediction error without the positive and negative values cancellation.Nine of the optimal MAE values of floods were obtained from the IP-DPE method, and As shown in Table 2 and Figure 9, the flood forecasting results obtained from the three methods are all technically reasonable but do show marked differences in forecasting accuracy.The MAE value can well reflect the actual prediction error without the positive and negative values cancellation.Nine of the optimal MAE values of floods were obtained from the IP-DPE method, and one optimum was obtained from the FC method.The MV method had the worst estimation results of six floods, and the simulated MAE of flood No. 090618 was almost three times that of IP-DPE method.The MAE range of the MV method was 16.177-216.32,that of the FC method was 16.294-160.83,and that of the IP-DPE method was 15.765-153.97.The IP-DPE method had the least fluctuation range and the error is generally low; in other words, IP-DPE had smaller fluctuations and was more stable than the other methods, and the MV method had the worst performance.
The ARE value can reflect the overall situation of the prediction error; nine optimal ARE values were obtained from the IP-DPE method and one optimal from the FC method.Seven of the MV method forecast results were the worst, and the simulated ARE of flood No. 091108 was nearly two times that of the IP-DPE method.The three floods estimated by the MV method fluctuated by more 10% in ARE, and one flood estimated by the FC method exceeded 10%.The ten floods estimated by the IP-DPE method fluctuated less than 9% in ARE values, again indicating that the proposed method has more stable performance.
Eight optimal PFRE values were obtained from the IP-DPE method.The estimation error of MV and FC methods for flood No. 100904 both exceeded 15%, indicating that their peak flow forecasting performance is poor.The IP-DPE method also has the optimal fluctuation range among them.
A comparison among the prediction results of the three methods is shown in Figure 10."Gaocun obs" refers to the observed values of Gaocun Hydrometric Station, while "Sunkou obs" refers to those of the Sunkou Hydrometric Station.
Water 2017, 9, 849 14 of 19 one optimum was obtained from the FC method.The MV method had the worst estimation results of six floods, and the simulated MAE of flood No. 090618 was almost three times that of IP-DPE method.The MAE range of the MV method was 16.177-216.32,that of the FC method was 16.294-160.83,and that of the IP-DPE method was 15.765-153.97.The IP-DPE method had the least fluctuation range and the error is generally low; in other words, IP-DPE had smaller fluctuations and was more stable than the other methods, and the MV method had the worst performance.
The ARE value can reflect the overall situation of the prediction error; nine optimal ARE values were obtained from the IP-DPE method and one optimal from the FC method.Seven of the MV method forecast results were the worst, and the simulated ARE of flood No. 091108 was nearly two times that of the IP-DPE method.The three floods estimated by the MV method fluctuated by more than 10% in ARE, and one flood estimated by the FC method exceeded 10%.The ten floods estimated by the IP-DPE method fluctuated less than 9% in ARE values, again indicating that the proposed method has more stable performance.
Eight optimal PFRE values were obtained from the IP-DPE method.The estimation error of MV and FC methods for flood No. 100904 both exceeded 15%, indicating that their peak flow forecasting performance is poor.The IP-DPE method also has the optimal fluctuation range among them.
A comparison among the prediction results of the three methods is shown in Figure 10."Gaocun obs" refers to the observed values of Gaocun Hydrometric Station, while "Sunkou obs" refers to those of the Sunkou Hydrometric Station.Our results indicate that the IP-DPE method had the best performance and the MV method had the poorest performance of the three methods we tested.
From Figure 11, we can find that flood has moderate attenuation from the upstream section to the downstream section, especially the decreases of flood peak and when the flood peak becomes chunky.In this regard, Zhang et al. [49] believe that if there is no bed scouring and siltation in the river channel, and the shapes of cross-sections are unchanged along the channel, the flow is assumed as clear water.Under the circumstance of constant roughness, the flood peak discharge attenuates and peak pattern grows chunky along the river because of the influence of the channel storage and side wall resistance.In the case that roughness decreases with time and the reduction rate of the upper reaches is greater than that of the lower reaches, the discharge increases abnormally along the river channel.In addition, the greater the reduction rate, the more obvious the peak increase.In this paper, we use the Muskingum model to carry out our research from the hydrological point of view.These effects have all been included in the model parameters, and the floods evolved by these parameters also verify this phenomenon.Our results indicate that the IP-DPE method had the best performance and the MV method had the poorest performance of the three methods we tested.
From Figure 11, we can find that the flood has moderate attenuation from the upstream section to the downstream section, especially the decreases of flood peak and when the flood peak becomes chunky.In this regard, Zhang et al. [49] believe that if there is no bed scouring and siltation in the river channel, and the shapes of cross-sections are unchanged along the channel, the flow is assumed as clear water.Under the circumstance of constant roughness, the flood peak discharge attenuates and peak pattern grows chunky along the river because of the influence of the channel storage and side wall resistance.In the case that roughness decreases with time and the reduction rate of the upper reaches is greater than that of the lower reaches, the discharge increases abnormally along the river channel.In addition, the greater the reduction rate, the more obvious the peak increase.In this paper, we use the Muskingum model to carry out our research from the hydrological point of view.These effects have all been included in the model parameters, and the floods evolved by these parameters also verify this phenomenon.Our results indicate that the IP-DPE method had the best performance and the MV method had the poorest performance of the three methods we tested.
From Figure 11, we can find that the flood has moderate attenuation from the upstream section to the downstream section, especially the decreases of flood peak and when the flood peak becomes chunky.In this regard, Zhang et al. [49] believe that if there is no bed scouring and siltation in the river channel, and the shapes of cross-sections are unchanged along the channel, the flow is assumed as clear water.Under the circumstance of constant roughness, the flood peak discharge attenuates and peak pattern grows chunky along the river because of the influence of the channel storage and side wall resistance.In the case that roughness decreases with time and the reduction rate of the upper reaches is greater than that of the lower reaches, the discharge increases abnormally along the river channel.In addition, the greater the reduction rate, the more obvious the peak increase.In this paper, we use the Muskingum model to carry out our research from the hydrological point of view.These effects have all been included in the model parameters, and the floods evolved by these parameters also verify this phenomenon.

Discussion
Theoretically, for the same river channel, we can suppose that the values K and x are constant in the linear Muskingum model.In practice, however, the factors that impact the flood forecasting accuracy may contain the uncertainty of data; the rainfall conditions before the flood, runoff in the river channel, flood magnitude, and discharge hydrograph vary under so many influence factors that using one of parameters in the same river channel is no longer applicable.Although the factors that cause the forecasting errors may be known, it is quite difficult to analyze their exact influence or the exact extent of said influence.We use the Neural Network black-box model to synthesize the characteristic attributes of the flood to output the necessary parameters in such a way that these influence factors are minimized.
We can see from Figure 11, in the process of validation, that there is a small part of floods presenting forecasts of peak timing.After the analysis, we found that the inaccurate results mainly appear in the prediction of multi-peak floods, and there is a tiny fluctuation near some parts of the flood peaks of the multi-peak flood.However, in the course of flood routing, the next period data was obtained by superimposing the data of the previous period with the data of this period.Therefore, it is possible that the peak timing is not completely consistent after the

Discussion
Theoretically, for the same river channel, we can suppose that the values K and x are constant in the linear Muskingum model.In practice, however, the factors that impact the flood forecasting accuracy may contain the uncertainty of data; the rainfall conditions before the flood, runoff in the river channel, flood magnitude, and discharge hydrograph vary under so many influence factors that using one group of parameters in the same river channel is no longer applicable.Although the factors that cause the forecasting errors may be known, it is quite difficult to analyze their exact influence or the exact extent of said influence.We use the Neural Network black-box model to synthesize the characteristic attributes of the flood to output the necessary parameters in such a way that these influence factors are minimized.
We can see from Figure 11, in the process of validation, that there is a small part of floods presenting inaccurate forecasts of peak timing.After the analysis, we found that the inaccurate results mainly appear in the prediction of multi-peak floods, and there is a tiny fluctuation near some parts of the flood peaks of the multi-peak flood.However, in the course of flood routing, the next period data was obtained by superimposing the data of the previous period with the data of this period.Therefore, it is possible that the peak timing is not completely consistent after the superposition of the current and the previous period.Although there are errors between the actual values and the predicted values, the errors are very small and will almost not affect the flood forecast.
Though the results obtained through the proposed method are more accurate than through the traditional method, IP-DPE is not perfect.For example, the model can only accurately reflect floods with the specific characteristic attributes that are substituted into it.The same attributes may not be available for a model of a different river channel.Although the current simulation accuracy is higher than the other methods, some of the flood characteristics are made available very close the flood peak, and there may also be altogether better attributes than those we used in this study.We plan to continue researching IP-DPE using other river channels as case studies and with different characteristic attributes for parameter estimation in order to continue improving the model's efficiency and efficacy.

Conclusions
In researching the more commonly used Muskingum model parameter estimation methods, we observed a tendency to set parameters according to mean values or flow magnitude, in which the practicability of the model and the applicability of the parameter are not fully considered.We assert that because each flood has its own particular parameters, that an "in-process type" dynamic parameter estimation method based on the PIO algorithm and BP-Artificial Neural Network model is the better choice for building a highly accurate Muskingum model.
After a brief literature review, we established the proposed Muskingum model parameter optimization model via the PIO algorithm.We then used data of occurring in the river channel from Chenggouwan to Linqing of Haihe River Basin in the year 1961 as an example and ran the PIO algorithm, L-S method, AGA method, and ACO algorithm to estimate Muskingum model parameters to validate the applicability and excellent performance of the PIO algorithm by comparison among them.
We also proposed the "in-process type" dynamic parameter estimation method for the Muskingum model based on the BP-Artificial Neural Network.According to the single-flood parameter estimation results, we extracted six characteristic attributes of each flood: IWL, PD, FPET, FVFT, FVBFP, and ADBFP, which were taken as input to obtain Muskingum model parameters K and x, which we used to establish the BP-ANN model.
We next took a channel of the Yellow River from Sunkou Hydrometric Station to Gaocun Hydrometric Station as a case study.We selected 60 floods which occurred over the past 40 years for analysis and applied the PIO algorithm to estimate the Muskingum model parameters of each flood.These results were used to obtain the correlation between the characteristic attributes and parameters K and x.We found that K decreases nonlinearly as each characteristic attribute increases, while x has a larger fluctuation range with no obvious rule.Therefore, if we select the mean parameter values, or take the parameter values in accordance with peak discharge magnitude as the input of the river flood forecasting model, it will cause an error.
Finally, we established IP-DPE for the Muskingum model based on the "6-5-2 type" BP-Artificial Neural Network.We used the first 50 floods from the dataset for model training and the last ten for validation.We also ran the MV and FC parameter estimation methods for performance comparison.We took the average absolute error, average relative error, and peak flow relative error as evaluation criteria to evaluate the simulation performance of the three methods and give out the flood simulation results of the last ten floods.The results showed that the proposed IP-DPE method had the most satisfying simulation results and best performance.

Figure 1 .
Figure 1.Absolute error results of four algorithms.

Figure 2 .
Figure 2. Relative error results of four algorithms.

Figure 1 .
Figure 1.Absolute error results of four algorithms.

Figure 1 .
Figure 1.Absolute error results of four algorithms.

Figure 2 .
Figure 2. Relative error results of four algorithms.

Figure 2 .
Figure 2. Relative error results of four algorithms.

Figure 1 .
Figure 1.Absolute error results of four algorithms.

Figure 2 .
Figure 2. Relative error results of four algorithms.

2. 3 .
"In-Process Type" Dynamic Parameter Estimation of the Muskingum Model Based on BP-Neural Network 2.3.1.Basic Principles of BP-Neural Network

Figure 5 .
Figure 5.The study river channel.

Figure 5 .
Figure 5.The study river channel.

Figure 5 .
Figure 5.The study river channel.

Figure 7 .
Figure 7. Correlation between each flood characteristic attribute and Muskingum model parameters.(a) Correlation between IWL and parameter K; (b) Correlation between IWL and parameter x; (c) Correlation between PD and parameter K; (d) Correlation between PD and parameter x; (e) Correlation between FPET and parameter K; (f) Correlation between FPET and parameter x; (g) Correlation between FVFT and parameter K; (h) Correlation between FVFT and parameter x; (i) Correlation between FVBFP and parameter K; (j) Correlation between FVBFP and parameter x; (k) Correlation between ADBFP and parameter K; (l) Correlation between ADBFP and parameter x.

Figure 7 .
Figure 7. Correlation between each flood characteristic attribute and Muskingum model parameters.(a) Correlation between IWL and parameter K; (b) Correlation between IWL and parameter x; (c) Correlation between PD and parameter K; (d) Correlation between PD and parameter x; (e) Correlation between FPET and parameter K; (f) Correlation between FPET and parameter x; (g) Correlation between FVFT and parameter K; (h) Correlation between FVFT and parameter x; (i) Correlation between FVBFP and parameter K; (j) Correlation between FVBFP and parameter x; (k) Correlation between ADBFP and parameter K; (l) Correlation between ADBFP and parameter x.

Figure 8 .
Figure 8.The flow chart of the "In-Process Type" Dynamic Parameter Estimation (IP-DPE) method.

Figure 8 .
Figure 8.The flow chart of the "In-Process Type" Dynamic Parameter Estimation (IP-DPE) method.
established the dynamic parameter estimation method of the Muskingum model based on the BP-Neural Network, as discussed above.The characteristic attributes of each flood were input to the model.The input layer of the Neural Network has six nodes, the hidden layer has five nodes, and the outputs of the Neural Network are the Muskingum model parameters K and x-this comprises the "6-5-2 type" BP-Neural Network model.There are 42 parameters in the Neural Network model, as established.Because the traditional Newton Iteration and Gradient methods easily fall into local Water 2017, 9, 849 13 of 19

Figure 10 .
Figure 10.Performance comparison between the three methods.(a) Mean absolute error comparison between three methods; (b) Average relative error comparison between three methods; (c) Peak flow relative error comparison between three methods.

Figure 10 .
Figure 10.Performance comparison between the three methods.(a) Mean absolute error comparison between three methods; (b) Average relative error comparison between three methods; (c) Peak flow relative error comparison between three methods.

Figure 10 .
Figure 10.Performance comparison between the three methods.(a) Mean absolute error comparison between three methods; (b) Average relative error comparison between three methods; (c) Peak flow relative error comparison between three methods.

Figure 11 .
Figure 11.Prediction results of three methods.(a) Prediction results of No. 071008 flood; (b) Prediction results of No. 080620 flood; (c) Prediction results of No. 090618 flood; (d) Prediction results of No. 090727 flood; (e) Prediction results of No. 090920 flood; (f) Prediction results of No. 091108 flood; (g) Prediction results of No. 100619 flood; (h) Prediction results of No. 100726 flood; (i) Prediction results of No. 100812 flood; (j) Prediction results of No. 100904 flood.

Figure 11 .
Figure 11.Prediction results of three methods.(a) Prediction results of No. 071008 flood; (b) Prediction results of No. 080620 flood; (c) Prediction results of No. 090618 flood; (d) Prediction results of No. 090727 flood; (e) Prediction results of No. 090920 flood; (f) Prediction results of No. 091108 flood; (g) Prediction results of No. 100619 flood; (h) Prediction results of No. 100726 flood; (i) Prediction results of No. 100812 flood; (j) Prediction results of No. 100904 flood.

Table 1 .
Parameter estimation results and performance comparison among four algorithms.

Table 2 .
Parameter estimation results obtained from mean value (MV) method, flow classification (FC) method, and IP-DPE method.

Table 2 .
Parameter estimation results obtained from mean value (MV) method, flow classification (FC) method, and IP-DPE method.