Improvement in Ridge Coefficient Optimization Criterion for Ridge Estimation-Based Dynamic System Response Curve Method in Flood Forecasting

Abstract

The ridge estimation-based dynamic system response curve (DSRC-R) method, which is an improvement of the dynamic system response curve (DSRC) method via the ridge estimation method, has illustrated its good robustness. However, the optimization criterion for the ridge coefficient in the DSRC-R method still needs further study. In view of this, a new optimization criterion called the balance and random degree criterion considering the sum of squares of flow errors (BSR) is proposed in this paper according to the properties of model-simulated residuals. In this criterion, two indexes, namely, the random degree of simulated residuals and the balance degree of simulated residuals, are introduced to describe the independence and the zero mean property of simulated residuals, respectively. Therefore, the BSR criterion is constructed by combining the sum of squares of flow errors with the two indexes. The BSR criterion, L-curve criterion and the minimum sum of squares of flow errors (MSSFE) criterion are tested on both synthetic cases and real-data cases. The results show that the BSR criterion is better than the L-curve criterion in minimizing the sum of squares of flow residuals and increasing the ridge coefficient optimization speed. Moreover, the BSR criterion has an advantage over the MSSFE criterion in making the estimated rainfall error more stable.

1. Introduction

Flood forecasting, an important non-structural measure, plays an important role in regional flood control, flood warning, risk decision making, etc. [1,2,3]. The hydrological model simplifies and conceptualizes the flood process with a set of equations aiming at obtaining the outlet flow. However, not all problems can be solved with such a model as the flood forecasting accuracy is often hampered and influenced by many error factors existing in the hydrologic system, including the errors in the model inputs, the errors in the model initial condition, the errors in the model simplification and the errors in the model parameters. Therefore, many scholars have devoted themselves to the research of error correction methods. For example, the autoregressive (AR) model estimates the flow error existing in a certain forecasting period by using the correlation of error series, and it was later developed into improved methods such as the recursive autoregressive model and the forgetting factor recursive autoregressive model [4,5,6]; Kalman filtering (KF) technology is widely used to update hydrological element time series in flood forecasting, and many improved types have been gradually formed, including the extended Kalman filter (EKF) [7] and the ensemble Kalman filter (EnKF) [8]. Data assimilation technology has also shown satisfying results in improving the prediction accuracy of models, including dynamic identifiability analysis (DYNIA) [9], and the Bayesian recursive estimation technique (BaRE) [10]. Robust theory [11] and the comprehensive correction method [12] are also used in flood forecasting error correction. However, error correction technology still needs to be improved; for example, the autoregressive model assumes that the error series has a linear correlation, but its performance near the flood peak is often not satisfying [12]. Additionally, some studies have also shown that the EnKF does not perform very well when the error structure is far away from a Gaussian distribution [13]. Scholars are still trying to solve the above problems.

A new error correction method called the dynamic system response curve (DSRC) method has been proposed by Bao et al. [14]. This method constructs a feedback model which is conceptualized on updating the hydrologic element series by tracing back to the source of the error. With the help of the first-order Taylor linearization to approximate the hydrologic model, error correction is achieved by solving the corresponding equations using the least square method. The DSRC method was initially applied to correct single hydrological elements, including runoff [15], rainfall [16] and model state variables [17]. Then, it was used to correct several hydrological elements comprehensively [18]. However, some studies [17,18,19,20] found that the correction results are not always stable, reflected in the excessive correction of hydrological element series and none-smooth simulated flow hydrographs. To solve the above problems, the DSRC-R method was developed by Si et al. [19] from the point of the regularization known as ridge estimation, and this method has improved the stability of correction results to some degree. Nevertheless, the selection criterion of the ridge coefficient still needs further study and improvement. Previous studies [17,18,19] often chose the ridge coefficient based on the L-curve criterion [21]. However, it has been found in practice that two aspects still need attention. One is that the L-curve seems insufficient to reflect the properties of model-simulated residual errors including the independence and the zero mean property of simulated residuals, which hinders the performance of DSRC-R in some cases; the other is that the application of the L-curve criterion takes a long time, which is not conducive to the real-time performance of flood forecasting. The L-curve criterion involves derivative calculation, and the difference method makes the operation efficiency lower.

Therefore, in this paper, we analyze the L-curve criterion and introduce the concepts of the random degree of simulated residuals and the balance degree of simulated residuals to describe the properties of the model-simulated residuals, and then a new criterion called the balance and random degree criterion considering the sum of squares of flow errors (BSR) is proposed. The new criterion takes the independence and the zero mean property of simulated residuals into account, which is conducive to obtaining a ridge coefficient in line with the statistical characteristics of residuals. The new criterion does not involve derivative calculation; thus, it can greatly shorten the search time of the ridge coefficient in optimization, improve the operational efficiency and enhance the real-time flood forecasting performance.

2. Methodology

2.1. DSRC Method

The main idea behind the DSRC method is that it firstly retrieves the rainfall errors from the outlet flow errors, then updates the rainfall series and finally reruns the model with the updated rainfall series. In this method, given a hydrological model $Q=Q(P)$ that generates outlet flow Q as a function of rainfall P, the variation process of the outlet flow change by the unit perturbation in rainfall is called the system response curve. Based on this, multi-time system response curves form the system response matrix $S$, and this matrix sets up the relation between rainfall errors $\Delta P$ and flow errors $\Delta Q$; then, the estimation of rainfall errors ${\widehat{\Delta P}}_{LS}$ can be computed via the least square method. In this round, the input rainfall series is updated with ${\widehat{\Delta P}}_{LS}$, and then the model is rerun with the updated rainfall series to correct the forecasting results. In this study, the DSRC method was combined with the Xinanjiang (XAJ) model, a hydrological model constructed by Professor Zhao Renjun of Hohai University which is widely used in flood forecasting in humid areas of China [22]. Additionally, here, we mainly talk about the calculation process of the DSRC method shown in Figure 1 rather than the XAJ model, as this paper will introduce it in Section 3. Additionally, the theoretical derivation of the DSRC method can be gained from [16].

According to [16], the rainfall error estimation ${\widehat{\Delta P}}_{LS}$ is expressed as the following Equation (1):

$${\widehat{\Delta P}}_{LS}={\left(S^{T}S\right)}^{-1}{S}^T\Delta Q$$

(1)

where $\Delta Q$ is the model-simulated deviation series of the outlet flow, $\Delta Q=Q_O-Q_C$; $Q_O={[Q_{O,1},Q_{O,2},\cdots Q_{O,M}]}^T$ is the observed flow series; and $Q_C={[Q_{C,1},Q_{C,2},\cdots Q_{C,\mathrm{M}}]}^T$ is the simulated flow series computed from the observed rainfall series.

In Equation (1), $S$ is the system response matrix defined as

$$S=\left[\begin{array}{cccc}\frac{\partial Q_1(P)}{\partial p_1}& \frac{\partial Q_1(P)}{\partial p_2}& \cdots & \frac{\partial Q_1(P)}{\partial p_N}\\ \frac{\partial Q_2(P)}{\partial p_1}& \frac{\partial Q_2(P)}{\partial p_2}& \cdots & \frac{\partial Q_2(P)}{\partial p_N}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial Q_M(P)}{\partial p_1}& \frac{\partial Q_M(P)}{\partial p_2}& \cdots & \frac{\partial Q_M(P)}{\partial p_N}\end{array}\right]$$

(2)

where $p_1,p_2\cdots p_n$ are the initial rainfall values; and the indices M and N represent the lengths of observed flow and rainfall, respectively (M $\ge$ N). $\frac{\partial Q_i(P)}{\partial p_j}$ represents the influence of the j-th rainfall on the i-th outlet flow. When i < j, it is obvious that $\frac{\partial Q_i(P)}{\partial p_j}=0$ because rainfall does not affect the outlet flow that occurs before it. $\frac{\partial Q_i(P)}{\partial p_j}$ is generally obtained by the difference, which is $\frac{\partial Q_i(P)}{\partial p_j}=\frac{Q(p_1,\dots ,p_j+\Delta p,\dots ,p_N)-Q(p_1,\dots ,p_j,\dots ,p_N)}{\Delta p}$.

Accordingly, the updated rainfall series $P^{’}$ can be expressed as Equation (3):

$$P^{’}=P_O+\widehat{\Delta P}$$

(3)

where $P^{’}$ is the updated rainfall series; $P_O$ is the initial rainfall series; and $\widehat{\Delta P}$ is the rainfall error estimation series. In the DSRC method, $\widehat{\Delta P}$ is replaced by ${\widehat{\Delta P}}_{LS}$ mentioned above.

In order to improve the forecasting accuracy, the updated rainfall series $P^{’}$ is introduced into the hydrological model for recalculation, and then the updated model-simulated flow series $Q_C^{’}$ is obtained by Equation (4):

$$Q_C^{’}=Q(P^{’})$$

(4)

where $Q_C^{’}={[Q_{C,1}^{’},Q_{C,2}^{’},\cdots ,Q_{C,N}^{’}]}^T$ is the updated model-simulated flow series.

2.2. DSRC-R Method

Correction results from the DSRC method are sometimes unstable. Relevant studies [19,20] have pointed out that the DSRC method is prone to be ill-conditioned, which generates unstable results when this method is applied to small basins, or when the length of flow information is short. Therefore, Si et al. [19] combined the DSRC method with the ridge estimation method and proposed the DSRC-R method, which is more robust than the DSRC method. In previous studies, the ridge coefficient $\beta$ was often selected via the L-curve criterion [17,18,19]; however, the correction results were not always stable. Thus, the criterion for obtaining the appropriate ridge coefficient $\beta$ in the DSRC-R method needs further study. Here, we directly provide the formula of rainfall error estimation ${\widehat{\Delta P}}_{RE}$ as the following Equation (5). For more details about the derivation process of DSRC-R, please refer to [19].

$${\widehat{\Delta P}}_{RE}={\left(S^{T}S+\beta I\right)}^{-1}{S}^T\Delta Q$$

(5)

where $\beta$ is the ridge coefficient; $I$ is the identity matrix; and ${\widehat{\Delta P}}_{RE}$ is the rainfall error estimation series of the DSRC-R method.

${\widehat{\Delta P}}_{RE}$ can be introduced into Equation (3) to update the rainfall series, and then the model can be rerun with $P^{’}$ to correct the forecasted flow. The flow chart of the DSRC-R method is shown in Figure 2.

2.3. L-Curve Criterion

The ridge coefficient $\beta$ is significantly important in the DSRC-R method. In some previous studies [17,18,19], the L-curve criterion was used to determine the value of $\beta$. According to [23], the main idea behind the L-curve criterion in selecting the appropriate ridge estimation coefficient $\beta$ in the DSRC-R method can be summarized as: the balance of $\lg \Vert Q(P_O+{\widehat{\Delta P}}_{RE})-Q_O\Vert$ and $\log \Vert {\widehat{\Delta P}}_{RE}\Vert$, and the appropriate ridge coefficient lies at the corner of the curve, usually with the largest curvature. That is, the technique to find the appropriate ridge coefficient can be expressed as Equation (6).

$$\underset{\beta }{\max }\left(\frac{\left|f^{’}{g}^{’’}-f^{’’}{g}^{’}\right|}{{\left[{(f^{’})}^2+{(g^{’})}^2\right]}^{3/2}}\right)$$

(6)

where $f(\beta )=\log {\Vert {\widehat{\Delta P}}_{RE}\Vert }_2^2$; $g(\beta )=\lg {\Vert Q(P_O+{\widehat{\Delta P}}_{RE})-Q_O\Vert }_2^2$; and ${\Vert •\Vert }_2^2$ is the modular square.

It can be proved (see Appendix A and Appendix B) that $f(\beta )$ and $g(\beta )$ can be expressed as Equations (7) and (8).

$$f(\beta )=\log [{\displaystyle \sum _{i=1}^n{(\frac{k_i}{{\lambda }_i+\beta })}^2}]$$

(7)

$$g(\beta )=\lg {\Vert Q(P_O+{\displaystyle \sum _{i=1}^n\frac{k_i}{{\lambda }_i+\beta }}{v}_i)-Q_O\Vert }_2^2$$

(8)

where ${\lambda }_i(i=1,\cdots ,N)$ is the eigenvalue of the matrix $S^{T}S$; $v_i(i=1,\cdots ,N)$, orthogonal to each other, is the unit eigenvector corresponding to ${\lambda }_i(i=1,\cdots ,N)$; and $k_i(i=1,\cdots ,n)$ is a group of coefficients that enable $S^T\Delta Q$ to be linearly expressed by $v_i(i=1,\cdots ,N)$, that is, $S^T\Delta Q={\displaystyle \sum _{i=1}^n{k}_i}{v}_i$.

The L-curve criterion has a good effect on selecting the ridge coefficient in the DSRC-R method, but there are still some problems that are worthy of attention. First, the result is sometimes unsatisfactory. The reason is that the L-curve criterion seems to insufficiently reflect the properties of the model-simulated residuals, although it pays attention to the balance of $\lg \Vert Q(P_O+{\widehat{\Delta P}}_{RE})-Q_O\Vert$ and $\log \Vert {\widehat{\Delta P}}_{RE}\Vert$. Second, the optimization of the ridge coefficient consumes too much time. The L-curve criterion involves first-order and second-order derivatives, as shown in Equation (6). Additionally, the explicit expression of $g(\beta )$ cannot be obtained at present, meaning its derivative can only be obtained by a difference method; thus, it will take a long time and is not conducive to the real-time performance of flood forecasting.

2.4. New Optimization Criterion (BSR)

In previous studies [17,18,19], the L-curve criterion was generally adopted to find the suitable ridge estimation coefficient $\beta$ in the DSRC-R method. Nevertheless, the L-curve criterion has some shortcomings such as insufficient consideration of model-simulated residuals, imperfect utilization of information and huge consumption of time. Therefore, this study takes the properties of the model-simulated residuals into consideration and then explores a new optimization criterion which is more suitable for the DSRC-R method.

For any model, it is always expected that the simulated residual series satisfies the zero mean property and non-correlative statistical property; in other words, let the mean of the residual series and correlation coefficient be as small as possible. This shows us that the criterion for determining the ridge coefficient should consider Equations (9) and (10). We use Equation (9) to express the zero mean property of the residuals, and this indicator is called the balance degree of simulated residuals (BDSR). We use Equation (10) to express the correlation of the residual series, and its reciprocal is called the random degree of simulated residuals (RDSR), which is shown in Equation (11). Additionally, RDSR indicates the independence of the residuals, and we generally expect a smaller value of the BDSR indicator and a bigger value of the RDSR indicator.

$$BDSR=abs\left[{\displaystyle \sum _{i=1}^m\left(Q_{O,i}-Q_{C,i}^{’}\right)}\right]$$

(9)

$$r_e=abs[\frac{{\displaystyle \sum _{i=1}^{m-1}(\Delta Q_i^{’}-\overline{\Delta Q_X})(\Delta Q_{i+1}^{’}-\overline{\Delta Q_Y})}}{\sqrt{{\displaystyle \sum _{i=1}^{m-1}(\Delta Q_i^{’}-\overline{\Delta Q_X}}{)}^2{\displaystyle \sum _{i=1}^{m-1}(\Delta Q_{i+1}^{’}-\overline{\Delta Q_Y}}{)}^2}}]$$

(10)

$$RDSR=\frac{1}{r_e}$$

(11)

where abs represents the absolute value sign; $r_e$ is the absolute value of the correlation coefficient of adjacent residuals;$\Delta Q_i^{’}=Q_{O,i}-Q_{C,i}^{’}$; $\overline{\Delta Q_X}=\frac{1}{m-1}{\displaystyle \sum _{i=1}^{m-1}\Delta Q_i^{’}}$; and $\overline{\Delta Q_Y}=\frac{1}{m-1}{\displaystyle \sum _{i=1}^{m-1}\Delta Q_{i+1}^{’}}$.

Traditional methods often take the minimum sum of squares of flow errors (MSSFE) as the objective function for calibration parameters, but the effective information contained in this method is not sufficient to obtain the ridge coefficient. The reason is that the derivation process of the DSRC-R method utilizes the least square method, and when the system is linear, the least sum of squares of flow errors equals the least square method; thus, the value of $\beta$ should be zero in this circumstance. Although the DSRC-R method belongs to non-linear system inversion methods, the value of $\beta$ still has a decreasing tendency, and this is not conducive to the stability of the method. Therefore, we need to further excavate more useful information in the simulated errors.

In this paper, we take the independence and the zero mean property of simulated residuals into account, combine these two points with the sum of squares of flow errors (SSFE) and lastly explore a new criterion called the balance and random degree criterion considering the sum of squares of flow errors (BSR criterion). Generally, we hope to find a large value of RDSR which is more consistent with the property of the residuals and is conducive to avoiding system errors; moreover, we hope to find a small value of BDSR which can satisfy the zero mean property of residuals and can decrease the flood volume errors. Based on this, considering the SSFE indicator, we propose the BSR criterion, the mathematical form of which is provided in Equation (12). Overall, the new BSR criterion considers the traditional sum of squares of residuals; furthermore, it is possible to find the ridge coefficient which satisfies the properties of flow residuals.

$$\underset{\beta }{\min }\frac{(BDSR+1)SSFE}{RDSR}$$

(12)

where “+1” is used to avoid the value of BDSR being zero; and $SSFE={{\displaystyle \sum _{i=1}^m\left(Q_{O,i}-Q_{C,i}^{’}\right)}}^2$.

Since the BSR criterion pays more attention to the independence and the zero mean property of simulated residuals than the L-curve criterion, it is more likely to select a ridge coefficient $\beta$ which satisfies the properties of residuals, and thus a better performance can be achieved with the DSRC-R method; moreover, the BSR criterion does not involve derivation calculation, meaning it can improve operational efficiency and save much more time.

In real-time flood forecasting, in order to quickly obtain an appropriate $\beta$, we need to utilize an automatic optimization method. Additionally, in this paper, we adopt the particle swarm optimization algorithm, which was first proposed by Kennedy and Eberhart and constructed on the concept of mimicking the social behavior of birds [24,25,26,27]. This algorithm has been widely used in many types of optimization problems. For more details about particle swarm algorithms, please refer to [24].

2.5. The Entire Research Process

In order to make the whole research process more clear, we created a flow chart, as shown in Figure 3. The figure shows the entire research process including the proposal of the BSR criterion, the research of a synthetic case and a real case and the comparison of three criteria (BSR, L-curve and MSSFE).

3. Case Study

This research projecting synthetic and real-data studies aimed at comparing the performance of the DSRC-R method under three criteria which include the BSR criterion, the L-curve criterion and the MSSFE criterion.

3.1. Model Description

The selected hydrological model in this research is the Xinanjiang (XAJ) model by Professor Zhao Renjun of Hohai University [22], which is one of the most widely used conceptual hydrological models in China. The XAJ model, including the inputs of observed precipitation as well as pan evaporation, and the outputs of forecasted flow as well as evaporation, can be used in different spatial and temporal scales and be divided into four layers: the first layer utilizing the three-layer evapotranspiration (TLE) model to realize basin evaporation; the second layer utilizing the saturated runoff production (SRP) model to realize the basin runoff production; the third layer utilizing the free water storage model to realize runoff separation; and the fourth layer utilizing the linear reservoir method and the Muskingum method to realize the basin flow concentration. When applying the XAJ model, firstly, divide the basin into several sub-basins and then compute the runoff and outlet flow in every sub-basin; lastly, gather the flow of each sub-basin at the outlet of the basin. For more details about the XAJ model, please refer to [22]. The structure of the XAJ model is shown in Figure 4, and the meaning of parameters in each layer is shown in

Table 1. The parameters of the XAJ model.

Layer	Function	Parameter	Meaning
First layer	Evaporation	K	Ratio of potential evapotranspiration to pan evaporation
		WUM	Areal mean tension water capacity of the upper layer
		WLM	Areal mean tension water capacity of the lower layer
		WDM	Areal mean tension water capacity of the deeper layer
		C	Coefficient of deep evapotranspiration
Second layer	Runoff production	IM	Ratio of impervious area
		WM	Areal mean tension water capacity
		B	Exponent of the tension water capacity distribution curve
Third layer	Runoff separation	SM	Areal mean free water capacity of the surface soil layer
		EX	Exponent of the free water capacity curve
		KI	Outflow coefficients of the free water storage to interflow
		KG	Outflow coefficients of the free water storage to groundwater
Fourth layer	Flow Concentration	CS	Recession constant of the surface water storage
		CI	Recession constant of the interflow storage
		CG	Recession constant of the groundwater storage
		KE	Storage time constant
		XE	Weight factor

.

3.2. Synthetic Case

The synthetic case was to design a typical artificial basin whose catchment area, station distribution, model parameters, hydrological features of the basin condition, “observed” precipitation, evaporation, “observed” outlet flow, error factors and any other information about the basin are all known in order to compare different schemes expediently. The synthetic case included 10,000 synthetic precipitations and corresponding floods in order to compare the performance of the DSRC-R method under three criteria (the BSR criterion, the L-curve criterion and the MSSFE criterion).

A major point of the synthetic case was to obtain the “observed” flow. Here, we referred to [16] and utilized Equation (13) to obtain it. By using different $P_O$ and $\Delta P$, we can obtain a different $Q_O$.

$$Q_O=Q(P_O+\Delta P)+e$$

(13)

where $Q_O$ is the “observed” flow series; $P_O$ is the initial precipitation series; $\Delta P$ is the given error series, and each value in the series $\Delta P$ does not exceed 30% of the corresponding value in $P_O$; and $e$ is Gaussian white noise which cannot exceed 5% of the initial value in this study.

This synthetic case assumed that the basin area is 1000 km² and there are 8 precipitation stations in the basin. The value of parameters from each layer is shown in the following

Table 2. The parameters of the XAJ model in the synthetic basin.

Parameter	K	WM	WUM	WLM	WDM	IM	B	C	SM
Value	1.1	150	20	80	50	0.01	0.3	0.16	10
Parameter	EX	KI	KG	CS	CI	CG	KE	XE
Value	1.5	0.35	0.35	0.78	0.865	0.995	1.50	0.380

.

3.2.1. Data

A major point of the synthetic case was to generate different initial series of precipitations $P_O$. In order to increase the diversity of $P_O$, we applied the following method. Firstly, we chose 55 typical areal precipitation processes from a real basin, then transformed the position of the rainfall peak in each precipitation and eventually formed 500 synthetic typical precipitation processes. When generating synthetic rainfall, we selected a synthetic typical precipitation, randomly adjusted each rainfall period ranging less than 30% of the typical rainfall and then obtained the proportion of each time interval of the rainfall series; then, we randomly generated the total rainfall and allocated it to each time interval according to the above proportion, thus forming one synthetic precipitation, that is, one initial series of precipitation $P_O$. Then, we introduced $P_O$ and the given $\Delta P$ into Equation (13) to obtain the “observed” flow series. In this case, we constructed 10,000 synthetic precipitations and corresponding floods, and the total rainfall of each flood ranged from 10 to 200 mm.

3.2.2. Statistical Indicators

In the synthetic case, what we consider most is the performance of the DSRC-R method under different criteria rather than the contrast of the results between the DSRC and DSRC-R methods, as was accomplished in [19]. The criteria include the BSR criterion, the L-curve criterion and the MSSFE criterion. The relevant statistical indicators include the relative error of flood peak (RPF), relative error of runoff depth (RRD), Nash–Sutcliffe effiency coefficent (NSE), time needed to update a flood (TU) and root mean square error (RMSE). RMSE can be utilized to evaluate the robustness of the DSRC-R method under different criteria, and this index was applied in [17]. The smaller the value of the RMSE indicator, the more robust the DSRC-R method will be.

The statistical indicators between synthetic cases and real-data cases are different. In synthetic cases, RMSE is one of the indicators; however, it is not covered in real-data cases as there is no way to obtain the “true” precipitation. The mathematical definitions of each statistical indicator are expressed as follows:

Relative error of flood peak (RPF):

$$\mathrm{RPF} =(Q_{OP}-Q_{CP})/Q_{OP}\times 100\%$$

(14)

Relative error of runoff depth (RRD):

$$\mathrm{RRD} =(R_O-R_C)/R_O\times 100\%$$

(15)

Nash–Sutcliffe effiency coefficent (NSE):

$$\mathrm{NSE}=1-\frac{{\displaystyle \sum _{i=1}^N{(Q_{C,i}-Q_{O,i})}^2}}{{\displaystyle \sum _{i=1}^N{(Q_{O,i}-{\overline{Q}}_O)}^2}}$$

(16)

Root mean square error (RMSE):

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{(p_{T,i}-p_i^{’})}^2}}$$

(17)

where $Q_{CP}$ is the forecasted value of the flood peak; $Q_{OP}$ is the observed value of the flood peak; $R_C$ is the forecasted depth of runoff; $R_O$ is the observed depth of runoff; ${\overline{Q}}_O$ is the average value of flow; $Q_{C,i}$ and $Q_{O,i}$ are the forecasted flow and observed flow in the i-th time interval; and $p_{T,i}$ and $p_i^{’}$ are the original precipitation and the updated precipitation in the i-th time interval. In the synthetic case, $p_{T,i}=p_{O,i}+\Delta p_i$, and $p_i^{’}=p_{O,i}+{\widehat{\Delta p}}_i$. The larger the NSE (NSE $\le 1$), the higher the forecasting accuracy, the smaller the RMSE and the more robust the DSRC-R method will be.

3.2.3. Computational Process of DSRC Method and DSRC-R Method

The mechanism of the DSRC method is that it firstly utilizes the error information of the outlet flow to invert and estimate the rainfall error, then updates the original rainfall series and lastly reruns the model with the updated rainfall series to correct the forecasting result. The specific steps of the DSRC method are as follows:

• Add the additional precipitation in the i-th time interval $\Delta p_i$ to the precipitation in the i-th time interval $p_{O,i}$ while keeping the precipitation in the j-th time interval ($j\ne i$)$p_{O,j}$ unchanged; then, obtain the new precipitation series $P_O+\Delta p_i$.
• Introduce the original precipitation series $P_O$ and new precipitation series $P_O+\Delta p_i$ into the model and obtain the series $Q(P_O)$ and $Q(P_O+\Delta p_i)$, respectively. Then, $S_i$ is obtained by the equation $S_i=\left[Q(P_O+\Delta p_i)-Q(P_O)\right]/\Delta p_i$, where $S_i$ is the dynamic system response curve of the i-th rainfall, that is, the i-th column of matrix S.
• Cycle Steps 1 and 2 n times and obtain the precipitation dynamic system response matrix $S$.
• Add the estimated precipitation error series ${\widehat{\Delta P}}_{LS}$ to the original precipitation series $P_O$ and obtain the updated precipitation series $P^{’}$.
• Introduce the updated precipitation series $P^{’}$ into the model in order to obtain the updated forecasted flow $Q_C^{’}$.

According to relevant research [19,20], when the flow data are insufficient, the DSRC method will tend to be unstable, characterized by wide fluctuations in the precipitation error estimated series and potentially oscillation, which will influence the flow correction effect. Therefore, a more robust method, DSRC-R, is proposed through combination with the ridge estimation method. Although the DSRC-R method ensures the stability of the error estimate, the ridge coefficient varies with floods. Thus, the ridge coefficient selection should be in accordance with the instant data. The steps of ridge coefficient optimization are as follows:

• Initialize the ridge coefficient $\beta$;
• Obtain the precipitation error estimation series ${\widehat{\Delta P}}_{RE}={\left(S^{T}S+\beta I\right)}^{-1}{S}^T\Delta Q$. Add ${\widehat{\Delta P}}_{RE}$ to $P_O$ in order to obtain $P^{’}$;
• Rerun the model with $P^{’}$ and obtain the updated flow process $Q_C^{’}$;
• Judge whether the results meet the criteria (the criteria adopted in this essay include the BSR criterion, the L-curve criterion and the MSSFE criterion). If yes, turn to Step 6; if no, go back to Step 5;
• Adjust the ridge coefficient $\beta$ according to the optimization algorithm (this essay applied the particle swarm optimization algorithm) and then turn to Step 2;
• Finish the optimization process and acquire the optimal ridge coefficient $\beta$.

In order to evaluate the correction effect of the three criteria, the synthetic case not only includes some flow indicators (RPF, RRD, NSE) but also includes some precipitation indicators (RMSE). Additionally, RMSE is proposed to quantitatively describe the robustness of the method and the oscillation phenomena. The smaller the RMSE, the more robust the precipitation error estimation will be. RMSE has been applied in previous studies and has been successful. However, in the actual case, due to the inability to acquire the “true” value of rainfall, the performance of the three criteria cannot be evaluated by the RMSE indicator.

3.2.4. Results and Discussion

The performance of the DSRC-R method on 10,000 synthetic floods under three criteria is shown in

Table 3. The results of the synthetic case.

Items 1	RPF	RRD	NSE	RDSR	BDSR	TU	$\mathit{\beta }$	RMSE
Before correction	9.38	5.12	0.985	——	——	——	——	——
L-curve	2.19	1.30	0.998	5.57	80.8	12.21	985.11	0.941
MSSFE	1.97	1.27	0.999	9.31	74.8	3.99	64.26	1.040
BSR	1.90	1.01	0.999	12.88	66.9	4.11	821.35	0.759

¹ The values of indicators (RPF, RRD, NSE, RDSR, BDSR, TU, $\beta$, RMSE) in the table are the average values of 10,000 synthetic floods.

, and one typical synthetic flood is shown in Figure 5. As can be seen in Table 3, although the three criteria (BSR, L-curve and MSSFE) have a certain effect, some differences still exist. In terms of the indicators RPF and RRD, the BSR criterion has the best performance (1.90% and 1.01%). The MSSFE criterion takes second place in this regard (1.97% and 1.27%), and the L-curve criterion has the worst performance (2.19% and 1.30%). In terms of the indicator NSE, the BSR criterion and MSSFE criterion have the same result, with a value of 0.999, which outnumbers that of the L-curve criterion, with 0.001. In terms of operational efficiency, the L-curve criterion consumes much more time, where the average TU of a flood is 12.21 s. The value of the ridge coefficient $\beta$ under the three criteria has a big difference, and $\beta$ tends to be smaller (average value of 64.26) when it applies the MSSFE criterion. In terms of the indicator RMSE, the value under the BSR criterion is 0.759, which is significantly less than the value (1.040) under the MSSFE criterion, indicating that the BSR criterion is more conducive to improving the robustness of the DSRC-R method than the MSSFE criterion.

As it is depicted in Figure 5a–c, the estimated value of the rainfall error under the BSR criterion is the closest to the “true” value, the total difference is 0.8 mm (1.2 mm errors under the L-curve criterion and 0.9 mm errors under the MSSFE criterion), the correlation coefficient reaches 0.962 (0.875 under the L-curve criterion and 0.896 under the MSSFE criterion) and the spots are evenly distributed on both sides of the 1:1 line. All of the above contribute to the DSRC-R method achieving the best performance under the BSR criterion. As it is shown in Figure 5d, the indicators RPT and RRD have values of 0.9% and 0.6%, respectively, under the BSR criterion. RPT has a value of 2.7% and RRD a value of 1.4% under the MSSFE criterion. Lastly, RPT has a value of 4.5% and RRD a value of 3.8% under the L-curve criterion. In a typical flood, the optimal value of $\beta$ under the MSSFE criterion is 2.47, which is significantly less than 651.33 under the BSR criterion. This result shows that the instability of the method is not sufficiently alleviated. Therefore, the points in Figure 5a are scattered. The L-curve criterion does not fully consider the properties of the simulated residuals, meaning the flow correction result in Figure 5d is not satisfactory.

Compared with the L-curve criterion, the BSR criterion improves the performance of the DSRC-R method. This is because the BSR criterion takes more account of the properties of model-simulated residuals, including the mutual independence between residuals and the zero mean property of residuals. As is shown in Table 3, the average BDSR under the BSR criterion is 66.9 m³/s, which is significantly less than 80.8 m³/s under the L-curve criterion, indicating that the former can better reflect the zero mean property of the simulated residuals; moreover, the average RDSR under the BSR criterion is 12.88, which is significantly greater than 5.57 under the L-curve criterion, indicating that the former is more conducive to the mutual independence of simulated residuals. The above analysis shows that the BSR criterion is more conducive to the optimization of the ridge coefficient, in line with the properties of simulated residuals, than the L-curve criterion. Therefore, the BSR criterion improves the performance of the DSRC-R method. In terms of operation efficiency, the BSR criterion takes less time than the L-curve criterion because the former does not involve derivative calculation, while the latter involves derivative difference calculation.

Compared with the BSR criterion, the MSSFE criterion tends to make the value of $\beta$ smaller (average value of 64.26), which is not conducive to improving the robustness of the DSRC-R method. The average RSME under the MSSFE criterion is 1.040, which is significantly greater than 0.759 under the BSR criterion, indicating that the BSR criterion makes the DSRC-R method more stable. This is because the MSSFE method is equivalent to the least square method in linear systems, meaning the value of $\beta$ should be 0 if the DSRC-R method is applied in a linear system. Although the XAJ model is a non-linear system, the value of $\beta$ will still tend to be small, which is not conducive to improving ill-conditioned problems. However, the BSR criterion introduces BDSR and RDSR, and this makes more use of the effective information contained in the simulated errors, which is conducive to avoiding a value of $\beta$ that is too small, making the DSRC-R method more stable.

To sum up, the BSR criterion has greater advantages among the three criteria. It pays more attention to extracting effective information from simulated errors of the outlet flow (in fact, any type of error will eventually be reflected here). It takes more account of the mutual independence and zero mean property of residuals, which is conducive to selecting a more reasonable value of $\beta$ and improving the performance of the DSRC-R method.

3.3. Real Case

The research basin in this study is Tankeng basin, with a total area of 3330 km² and 15 precipitation stations, which is located at the tributary of the Ou River in Zhejiang Province. Tankeng basin is in the subtropical monsoon climate zone, and it enjoys a temperate climate with well-marked seasons and plenty of rainfall and sunshine. Runoff from Tankeng basin is mainly supplied by precipitation. The annual average precipitation of Tankeng basin is between 1500 and 2100 mm, and multi-annual average evaporation is 969.9 mm. The precipitation spatial distribution is uneven, as is the annual precipitation temporal distribution; based on this, the entire year can be divided into three parts with the first part “spring rain” ranging from March to April, the second part “plum rain” ranging from May to June and the third part “thunderstorm” ranging from July to September. In Tankeng basin, multi-annual average flow is 120 m³/s, and the maximal runoff happens in June, with the proportion of the entire annual runoff reaching up to 19.7%. More details about Tankeng basin and its rainfall station distribution are illustrated in Figure 6.

3.3.1. Data

There are two time scales for observed data including the day scale and the hour scale. Evaporation data are day-scale information, while precipitation and flow data are both day-scale information and hour-scale information. This research collected observed data such as evaporation, observed flow and precipitation from 1980 to 2005. The rainfall data were from 15 rainfall stations, evaporation data came from one evaporation station and observed flow data came from the outlet flow of Tankeng basin.

Tankeng basin has been used in previous studies using the XAJ model [20,28], and thus no calibration was required here. The parameters of the XAJ model at Tankeng basin are shown in

Table 4. Parameters of the XAJ model for Tankeng basin.

Parameter	K	WM	WUM	WLM	WDM	IM	B	C	SM	EX
Value	1.296	150	20	80	50	0.01	0.3	0.16	10	1.5
Parameter	KI	KG	CS	CI	CG	KE	XE
Value	0.35	0.35	0.65	0.865	0.95	1.466	0.380

.

3.3.2. Results and Discussion

The statistical indicators of the real case are different from those of the synthetic case. In the real case, the “true” value of rainfall cannot be obtained, meaning RMSE cannot be applied. Here, the indicators include RPF, RRD, NSE and TU.

In this study, 31 historical floods in the basin were selected to compare the results of three criteria including the L-curve criterion, MSSFE criterion and BSR criterion. In order to show the performance of the DSRC-R method under each criterion, the statistical indicators of each flood are listed in

Table A1. The performance of the DSRC-R method on each flood in Tankeng basin.

Flood Number	Before Updating			L-Curve						MSSFE						BSR
	RPF	RRD	NSE	RPF	RRD	NSE	BDSR	RDSR	TU	RPF	RRD	NSE	TU	BDSR	RDSR	RPF	RRD	NSE	TU	BDSR	RDSR
31850823	10.1	−14.41	0.934	6.42	−9.26	0.974	468.0	1.32	25.2	6.25	−8.31	0.975	8.1	959.6	1.30	6.41	−7.38	0.985	15.3	6.2	1.33
31860520	2.77	−11.18	0.743	2.56	9.71	0.972	2958.1	1.59	42.5	2.52	8.89	0.973	21.1	2718.9	1.83	1.38	−0.02	0.954	18.7	2461.2	3.25
31870910	−11.7	−7.38	0.952	8.93	3.38	0.963	1256.9	1.06	29.4	−1.72	−4.77	0.986	6.4	1418.6	1.09	0.74	−3.35	0.984	9.9	40.6	2.35
31880520	−16.44	−10.91	0.762	6.93	7.49	0.888	2070.8	1.10	43.3	0.95	−6.17	0.903	11.1	236.0	1.31	0.65	−6.57	0.904	10.0	32.1	1.48
31880922	−5.48	−13.33	0.755	−1.19	6.99	0.976	2314.3	1.13	28.6	−0.66	5.06	0.979	30.6	2530.1	1.16	1.15	0.03	0.979	35.3	149.5	1.56
31890527	−17.28	−3.54	0.805	3.65	3.97	0.971	673.3	1.05	46.6	2.29	0.56	0.971	9.2	11.1	1.02	2.25	0.54	0.975	14.7	5.5	1.33
31890618	−15.11	−14.14	0.753	6.26	12.78	0.966	1210.2	1.08	24.9	−4.18	−3.05	0.977	7.8	2479.9	1.05	1.63	−0.02	0.973	8.8	2.78	1.13
31890721	−13.78	−17.3	0.866	1.87	7.78	0.982	1387.1	1.24	32.8	1.46	6.31	0.984	10.0	1623.5	1.26	0.76	−0.12	0.986	13.3	49.0	1.35
31920830	−13.34	9.85	0.959	−1.57	6.98	0.992	1847.3	1.40	39.3	−1.82	6.88	0.992	15.5	1304.8	1.61	−2.94	6.5	0.992	17.0	1161.8	1.73
31940821	13.96	−15.18	0.935	4.64	1.61	0.991	2249.6	1.29	17.6	4.76	1.31	0.991	9.5	1934.0	1.33	5.76	−0.82	0.99	8.5	12.5	1.37
31960801	−17.42	−1.81	0.937	−2.26	0.22	0.971	4241.8	1.99	24.3	−2.31	0.2	0.971	5.4	4231.6	1.98	−1.9	0.4	0.971	9.7	4344.9	2.06
31970703	−4.69	−0.37	0.781	1.32	14.73	0.983	3240.7	1.15	29.4	2.12	5.89	0.983	8.0	3179.3	1.16	1.2	5.64	0.983	13.5	3068.3	1.17
31980618	−5.91	4.28	0.96	1.43	15.92	0.997	1120.8	2.18	26.4	1.43	2.41	0.997	10.9	1127.6	2.22	1.32	2.23	0.997	7.9	1012.9	3.78
31990525	−18.28	−8.19	0.917	10.67	7.84	0.947	400.8	1.06	31.7	−0.45	0.48	0.983	8.3	135.7	1.08	−0.84	0.22	0.983	7.6	2.2	1.09
31990711	−2.94	−8.89	0.867	2.23	3.46	0.967	1430.2	1.13	24.6	2.27	−0.54	0.967	9.9	1133.0	1.14	2.71	−2.64	0.99	6.4	8.3	1.14
31000609	−20.05	−7.3	0.946	−8.96	0.14	0.964	671.1	1.78	32.3	−10.56	−1.68	0.968	12.1	663.2	1.80	−9.03	0.03	0.971	10.3	623.6	1.80
31000823	8.75	68.53	0.907	3.14	7.29	0.991	1621.2	1.06	26.2	2.81	−1.37	0.992	15.6	990.1	1.63	2.81	−0.87	0.991	8.6	967.0	1.66
31030624	−11.61	58.79	0.866	−1.83	6.08	0.981	659.4	1.03	24.9	−1.07	−1.97	0.982	8.0	522.1	1.02	0.7	−0.81	0.986	5.2	32.6	1.89
31040812	19.29	36.59	0.84	20.29	5.93	0.896	1214.4	1.52	15.2	11.58	1.11	0.93	4.0	1203.5	1.49	7.68	−1.46	0.941	6.5	1205.4	1.53
31850604	−13.11	−17.32	0.754	6.49	3.3	0.934	573.5	1.14	11.7	5.14	2.3	0.935	5.3	443.1	1.15	3.87	1.31	0.935	6.8	4.4	1.28
31850626	−10.91	−3.17	0.821	−2.03	11.97	0.896	59.6	2.01	16.3	−2.93	−1.49	0.906	10.5	388.3	1.08	−2.61	−0.7	0.905	9.7	22.9	2.12
31860330	−27.8	−18.3	0.68	3.71	0.96	0.968	28.7	1.15	38.7	4.23	1.23	0.969	36.9	112.8	1.14	3.47	0.86	0.968	33.0	8.7	1.15
31860428	41.36	19.54	0.02	10.92	11.8	0.825	1488.3	1.09	19.7	9.04	10.89	0.828	6.8	517.8	1.32	9.02	10.87	0.849	5.9	515.1	1.32
31870411	−2.08	−6.7	0.724	13.18	5.89	0.875	2131.9	1.29	14.5	10.36	4.31	0.879	5.4	1975.7	1.22	6.34	0.99	0.879	14.1	1806.3	1.37
31870527	−23.87	−9.43	0.738	−10.41	3.49	0.847	298.6	1.19	16.3	−10.2	3.62	0.847	9.7	225.3	1.20	−15.5	0.5	0.836	16.8	225.1	1.20
31880327	−19.21	−6.61	0.869	−10.11	−4.66	0.901	404.8	1.02	24.8	−10.87	−4.94	0.901	7.6	168.0	1.06	−10.08	−4.65	0.903	8.5	8.36	1.06
31880524	8.66	−2.92	0.429	1.59	0.09	0.926	50.4	1.03	19.9	1.61	0.1	0.926	10.3	23.3	1.05	−1.94	−0.83	0.931	6.4	4.1	1.21
31880613	−29.34	−10.44	0.193	1.07	0.91	0.595	353.5	1.23	16.1	−0.91	−0.15	0.597	9.0	431.7	1.20	0.91	0.82	0.595	5.3	4.8	1.24
31880627	0.97	47.39	0.412	2.51	4.76	0.967	876.5	1.17	18.4	2.84	5.15	0.967	5.3	411.9	2.04	−1.09	−0.07	0.971	17.0	480.5	2.70
31890413	14.65	10.1	0.125	2.18	8.35	0.931	617.2	1.03	20.8	0.54	8.11	0.932	3.8	516.5	1.07	0.27	8.07	0.932	8.6	492.5	2.16
31980513	28.8	57.29	0.098	7.8	19.32	0.889	1438.5	1.03	13.8	2.63	17.14	0.892	8.1	1312.3	1.65	1.25	16.6	0.893	16.2	1332.4	2.58
Average 1	14.51	16.81	0.721	5.42	6.68	0.933	1269.6	1.28	25.7	3.95	4.08	0.938	10.7	1126.7	1.34	3.49	2.77	0.940	12.1	648.1	1.69

¹ The average values of RPF and RRD actually represent the average of the absolute values of RPF and RRD, respectively.

(see Appendix C). In order to more clearly compare the results of the BSR criterion and the other two criteria, scatter diagrams of each indicator of 31 floods were constructed, as shown in Figure 7. As it is illustrated in Table A1, in terms of the flood peak, the average RPF is 5.42% under the L-curve criterion, 3.95% under the MSSFE criterion and 3.49% under the BSR criterion; therefore, it is apparent that the BSR criterion has a better performance. In terms of the runoff depth, the average RRD is 2.77% under the BSR criterion, which is 3.91% and 1.31% lower, respectively, than that under the L-curve and MSSFE criteria. In terms of NSE, the BSR criterion has the maximal NSE of 0.940, while the NSE under the L-curve criterion is 0.933 and 0.938 under the MSSFE criterion. In terms of operational efficiency, the average time consumed under the BSR criterion is 12.7 s, which is 25.7 s lower than that under the L-curve criterion, and this illustrates that the operation efficiency is dramatically improved when the BSR criterion is adopted.

Figure 7a–c show the comparison of the statistical indicators of 31 floods between the BSR criterion and L-curve criterion. Figure 7a,b are inclined to the lower side of the 1:1 line, and Figure 7c is inclined to the upper side of the 1:1 line, indicating that the BSR criterion is more effective than the L-curve criterion for most floods. In terms of the flood peak and runoff depth, the BSR criterion is more effective in 24 floods and 29 floods, respectively; additionally, the NSE of 23 floods is greater under the BSR criterion. Figure 7d–f show the comparison of flood indicators between the BSR criterion and MSSFE criterion. Similar to Figure 7a–c, under the BSR criterion, the forecasting accuracy of 21 flood peaks and 25 flood runoff depths is higher, and NSE of 18 floods is larger, indicating that the BSR criterion has more advantages over the MSSFE criterion for most floods.

Compared with the L-curve criterion, the BSR criterion improves the performance of the DSRC-R method. This is because the BSR criterion introduces BDSR and RDSR, thus more fully considering the properties of the model-simulated residuals including the zero mean property and the mutual independence. The results in Table A1 show that the average BDSR corresponding to the BSR criterion is 648.1 m³/s, which is far less than 1269.9 m³/s of the L-curve criterion, indicating that the BSR criterion is more conducive to the zero mean property of the simulated residuals. At the same time, the average RDSR corresponding to the BSR criterion is 1.69, which is greater than 1.28 of the L-curve criterion, indicating that the BSR criterion tends to meet the mutual independence of the simulated residuals. Therefore, under the BSR criterion, the DSRC-R method has a better performance. In the real case, the mean value of $\beta$ under the MSSFE criterion is 327.67, which is far less than 2048.97 under the BSR criterion. This is because the MSSFE criterion only considers the sum of squares of errors, which tends to make $\beta$ small, while the BSR criterion focuses on the nature of the simulated residuals and extracts more effective information from the outlet flow error information, which is conducive to avoiding a small value of $\beta$. Therefore, the application of the BSR criterion in the real case is better than that of the MSSFE criterion.

4. Conclusions and Prospect

A previous study [19] proposed the DSRC-R method and verified that it has stronger robustness. However, the selection criterion of the ridge coefficient, usually the L-curve criterion [17,18,19], has only received limited attention. This essay constructed the BSR criterion based on the properties of model-simulated residuals, utilizing the indicator RDSR to quantitatively describe the independence of the residual series while utilizing the indicator BDSR to quantitatively describe the zero mean property of residuals. Additionally, we then contrasted the performance of the DSRC-R method under three different criteria (BSR criterion, L-curve criterion and MSSFE criterion) through synthetic and real-data studies.

From the results, we found that among the three criteria, the BSR criterion is more suitable for the DSRC-R method. Compared with the L-curve criterion, the BSR criterion improves the performance of the DSRC-R method. This is because the BSR criterion introduces RDSR and BDSR, which quantitatively describe the mutual independence of model-simulated residuals and zero mean property of model-simulated residuals, respectively. Moreover, the BSR criterion saves more time than the L-curve criterion because the BSR criterion does not involve derivative calculation. In addition, compared with the traditional MSSFE criterion, the BSR criterion is more conducive to enhancing the robustness of the DSRC-R method. The MSSFE criterion tends to make the ridge coefficient $\beta$ smaller, which is unfavorable to the performance of the DSRC-R method. Meanwhile, the BSR criterion is conducive to avoiding a small ridge coefficient by extracting more effective information contained in the simulated errors, and this makes the DSRC-R method more robust and improves its performance.

Further research is needed. The BSR criterion proposed in this paper improves the performance of the DSRC-R method, and this seems to benefit from the rational use of outlet flow information. In recent years, data assimilation technologies combined with radar information and remote sensing information have been continuously emerging and have received a significant amount of attention. However, making full use of outlet flow information to update hydrological elements deserves more attention. This is because hydrological elements such as rainfall, evaporation and soil moisture will eventually be reflected in the outlet flow. Therefore, much more effort needs to be made in this regard.