# A Conjunction Method of Wavelet Transform-Particle Swarm Optimization-Support Vector Machine for Streamflow Forecasting

Authors: Fanping Zhang; Huichao Dai; Deshan Tang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Streamflow forecasting has an important role in water resource management and reservoir operation. Support vector machine (SVM) is an appropriate and suitable method for streamflow prediction due to its best versatility, robustness, and effectiveness. In this study, a wavelet transform particle swarm optimization support vector machine (WT-PSO-SVM) model is proposed and applied for streamflow time series prediction. Firstly, the streamflow time series were decomposed into various details (D
s) and an approximation (A
3) at three resolution levels (21-22-23) using Daubechies (db3) discrete wavelet. Correlation coefficients between each D subtime series and original monthly streamflow time series are calculated. D
s components with high correlation coefficients (D
3) are added to the approximation (A
3) as the input values of the SVM model. Secondly, the PSO is employed to select the optimal parameters, C, ε, and σ, of the SVM model. Finally, the WT-PSO-SVM models are trained and tested by the monthly streamflow time series of Tangnaihai Station located in Yellow River upper stream from January 1956 to December 2008. The test results indicate that the WT-PSO-SVM approach provide a superior alternative to the single SVM model for forecasting monthly streamflow in situations without formulating models for internal structure of the watershed.


---

## Body

## 1. Introduction

The accuracy of streamflow forecasting is a key factor for reservoir operation and water resource management. However, streamflow is one of the most complex and difficult elements of the hydrological cycle due to the complexity of the atmospheric process. The elements affecting streamflow forecasting precision include catchment, geomorphologic and climate characteristics, and so forth [1]. The process of streamflow is extremely complex due to the influence of these variables and their combinations. Therefore, there are many forecasting techniques that have been proposed for streamflow forecasting [2–4].Among them, the most popular and widely known statistical method used in time series forecasting is autoregressive integrated moving average (ARIMA) model due to its superiority of forecasting capabilities and richness of information on time-related changes [5]. Several studies have shown that ARIMA can be trusted as a reliable model in water resources time series analysis [6]. For example, Lee and Tong [7] proposed a hybrid model for nonlinear time series forecasting by combining ARIMA and genetic programming and demonstrated the effectiveness of the proposed forecasting model. But the ARIMA models are a class of linear model and thus only suitable for capturing linear features of data time series [8]. In recent years, gray model, artificial neural network (ANN), and support vector machine (SVM) have been frequently used to predict the nonlinear time series and achieved good results [9–11]. For instance, Kişi [12] used three different ANN techniques, namely, feed forward neural networks, generalized regression neural networks and, radial basis ANN in one-month-ahead streamflow forecasting. However, there are some disadvantages of ANN due to its network structure, which is hard to determine and usually established using a trial-and-error approach [13].Support vector machines (SVM) were suggested by Vapnik [14] as one of the soft computational techniques and are widely used for classification and regression based on statistical learning theory (SLT). The basic idea of SVM for regression is to introduce a kernel function, map the input data into a high-dimensional feature space by a nonlinear mapping, and then perform linear regression in the feature space [5]. Currently, SVM were frequently applied in a number of different fields, such as fault diagnosis [16], pattern recognition [17], and classification [18]. In the hydrology context, SVM has been successfully applied to forecast the flood stage [19–21], to predict future water levels in Lake Erie [22], and to forecast discharges [23, 24]. Previous studies have indicated that SVM is an effective method for streamflow forecasting [5, 23–25].More recently, the conjunction model of wavelet and SVM has drawn increasing interest and has displayed advantages over a single SVM model in terms of prediction accuracy. Wavelet analysis (WA) is an advanced method proposed by Morlet et al. [26] in signal processing and has attracted much attention due to its ability to reveal simultaneously both spectral and temporal information within one signal [27]. The application of WA in the areas of hydrology and water resource research mainly includes these aspects: identification of hydrologic series deterministic components such as trend, periods, and change points [28–31]; wavelet denoising in hydrologic series [30, 31]; and hydrologic series simulation and prediction based on wavelet [27, 32, 33]. Wavelet analysis can be used to decompose an observed time series (such as streamflow time series) into various components so that the new time series can be used as inputs for SVM models [34].SVM implements the principle of structure risk minimization in place of experiential risk minimization, which makes it have excellent generalization ability in the situation of small sample. However, the practicability of SVM is affected by the difficulty of selecting appropriate SVM parameters [35]. At present, the most common parameters selection method for SVM is the cross validation method but it is time-consuming [36]. Recently, some intelligent algorithms have been applied for parameters selecting. Compared with cross validation, genetic algorithm (GA) is less time-consuming and can obtain the optimal solution well, but the operation of genetic algorithm is difficult with the steps of choosing, crossover, and mutation for different optimal problems [37]. As a new global optimizing algorithm, particle swarm optimization (PSO), proposed by Kennedy and Eberhart in 1995, is based on swarm intelligent by generating a random decision variable set called “particles” [35]. PSO is a versatile algorithm and can be used to solve different optimizing problems. In recent years, because of the best global searching ability and the simple implementing procedure, PSO has been successfully applied for function optimization [38], data mining [39], and other engineering optimization problems [15, 40] and achieved good results. Therefore, the PSO can be applied to optimize the parameters of SVM model for streamflow forecasting in this paper.This paper is organized as follows. Section2 introduces the principle theory of wavelet analysis, parameter selection method of SVM models based on PSO, and SVM regression forecasting model. The study area and streamflow time series analysis are introduced in Section 3. The forecasting results of the conjunction model with the real streamflow time series data sets from Tangnaihai hydrology station in China are analyzed in Section 4. Finally, the conclusion is presented in Section 5.


## 2. Methodology

### 2.1. Support Vector Machine (SVM)

The basic idea of SVM for regression is to introduce a kernel function, map the input data into a high-dimensional feature space by a nonlinear mapping, and then perform linear regression in the feature space [5]. Supposing that there is a training dataset D
=
{
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
…
,
(
x
n
,
y
n
)
}
∈
R
p
×
R, x is the input vector, y is the excepted output, n is the number of data, and p is the total number of data patterns. By nonlinear mapping function Ф, x is mapped into a feature space in which a linear estimate function is defined as
(1)
y
=
f
(
x
,
ω
)
=
〈
ω
,
Φ
(
x
)
〉
+
b
;
ω
,
x
∈
R
p
,
b
∈
R
,
where Ф
(
x
) represents the high-dimensional feature spaces, which is nonlinearly mapped from the input space x; ω and b are coefficients that have to be estimated from the input data. By introducing the slack variables ξ and ξ
* and following the regularization theory, parameters ω and b are estimated by minimizing the cost function
(2)
min
⁡
J
(
ω
,
ξ
)
=
1
2
∥
ω
∥
2
+
C
∑
i
=
1
n
(
ξ
i
+
ξ
i
*
)
subject to the constraints:
(3)
y
i
-
〈
ω
,
x
i
〉
-
b
≤
ε
+
ξ
i
〈
ω
,
x
i
〉
+
b
-
y
i
≤
ε
+
ξ
i
*
ξ
i
,
ξ
i
*
≥
0
.The first term(
1
/
2
)
∥
ω
∥
2 is weight vector norm; C is referred to as the regularized constraint determining the tradeoff between the empirical error and the regularized term; and ε is the insensitive loss function.By using Lagrange multiplier techniques, the minimization of (2) leads to the following dual-optimization problem:
(4)
max
⁡
W
(
α
i
,
α
i
*
)

=

-
1
2
∑
i
=
1
n
∑
j
=
1
n
[
(
α
i
-
α
i
*
)
(
α
j
-
α
j
*
)

K
(
x
i
,
x
j
)
]
11111111111111111
-
ε
∑
i
=
1
n
(
α
i
+
α
i
*
)
+
∑
i
=
1
n
y
i
(
α
i
-
α
i
*
)
s
.
t
.
∑
i
=
1
n
(
α
i
-
α
i
*
)
=
0
;
α
i
,
α
i
*
∈
[
0
,
C
]
,
where (
α
i
*
,
α
i
) are coefficients determined by training and K
(
x
i
,
x
j
) is the kernel function which can be expressed as inner product:
(5)
K
(
x
i
,
x
)
=
Φ
(
x
i
)
T
Φ
(
x
)
.The decision function takes the form(6)
f
(
x
)
=
∑
i
=
1
n
(
α
i
*
-
α
i
)
K
(
x
,
x
i
)
+
b
.The selection of an appropriate kernel function plays an important role in SVM regression since the kernel function defines the feature space. Gaussian radial basis function kernel has received significant attention from the machine learning community. Gaussian radial basis function (RBF) kernel is defined as(7)
K
(
x
,
x
i
)
=
exp
⁡
(
-
∥
x
-
x
i
∥
2
2
σ
2
)
.Hereσ is the kernel parameter.


### 2.2. Wavelet Analysis

In wavelet analysis, the signals are analyzed in both the time and the frequency domain by decomposing the original signals in different frequency bands using wavelet functions. The wavelet transform (WT) uses the scalable windowing technique for analyzing local variation in the time series [41]. WT provides useful decompositions of original time series, so that wavelet-transformed data improve the ability of a forecasting model by capturing useful information on various resolution levels [42]. The time series data are preprocessed using wavelet transformation techniques to obtain decomposed wavelet coefficients that are used as inputs in the forecasting models.The basic objective of WT is to achieve a complete timescale representation of localized and transient phenomena occurring at different timescales [4, 43]. The continuous wavelet transform is defined as the sum over all time of the signal multiplied by scale and shifted versions of wavelet function ψ:
(8)
W
(
a
,
b
)
=
1
a
∫
-
∞
+
∞
f
(
x
)
ψ
(
x
-
b
a
)
d
x
,
where a is a scale parameter; b is a position parameter; and ψ corresponds to the complex conjugate. The coefficient plots of the continuous wavelet transform are precisely the timescale view of the signal. However, calculating wavelet coefficients at every possible scale is time-consuming and generates large amount of information. Thus, the use of the continuous wavelet transform for forecasting is not practically possible.In hydrology, observed hydrologic series are often expressed as discrete series, so the discrete wavelet transform is usually employed to decompose a hydrologic series into a set of coefficients and subsignals under different scales, and then guide other time series analyses [44]. The DWT is defined as the following form:
(9)
f
(
t
)
=
∑
C
j
o
,
k
ϕ
j
o
,
k
(
t
)
+
∑
j
>
j
o
∑
ω
j
,
k
2
j
/
2
ψ
(
2
j
t
-
k
)
,
where j is the dilation or level index, k is the translation or scaling index, and ϕ
j
o
,
k is a scaling function of coarse scale coefficients. C
j
o
,
k, ω
j
,
k, is the scaling function of detailed (fine scale) coefficients and all functions of ψ
(
2
j
t
-
k
) are orthonormal.The original time series are decomposed into various details (D
s) and an approximation (A
s) at different resolution levels using DWT. The approximations are the high-scale, low frequency components of the signal and the details are the low-scale, high frequency components. Normally, the low frequency component of the signal (
A
) is the most important part which demonstrates the signal identity [45]. The choice of wavelet type is an important issue. The Daubechies wavelets are one of the widely used in wavelet family, which are written dbN, where db is the “surname” and N is the order of the wavelet [46]. Daubechies wavelets exhibit good tradeoff between parsimony and information richness [34], so in this study the Daubechies wavelets were employed as the mother wavelet to decompose the time series.


### 2.3. Parameters Selection of SVM Based on PSO

#### 2.3.1. The Principle of PSO

PSO, deriving from the research for the movement of organisms in a bird flocking or fish schooling, performs searches using a population (called swarm) of individuals (called particles) that are updated from iteration to iteration [47, 48]. An equation (velocity update) controls the swarm in moving around the search space seeking the optimum state. In each iteration, the algorithm saves the local optimum and compares it with the global (best yet) optimum values. Definitely the criteria for being chosen as an optimum state depend on the fitness of the objective function. Candidate solutions (decision variables) of any particle calculate and remember its own fitness. The position of any particle accelerated towards the global best position by using (10) and (11) [49]. In any search step t, the ith particle is used to update its candidate solution’s current position x
i
,
j
(
t
) by using local best p
i
,
j
(
t
) and best p
g
,
j
(
t
) position achieved yet. Consider the following:
(10)
v
i
,
j
(
t
+
1
)
=
ω
v
i
,
j
(
t
)
+
c
1
r
1
[
p
i
,
j
-
x
i
,
j
(
t
)
]
+
c
2
r
2
[
p
g
,
j
-
x
i
,
j
(
t
)
]
(11)
x
i
,
j
(
t
+
1
)
=
x
i
,
j
(
t
)
+
v
i
,
j
(
t
+
1
)
,
j
=
1,2
,
…
,
d
,
where v
i
,
j is velocity measures for particles; ω is inertial weight controlling velocity direction; c
1 and c
2 are acceleration coefficients; r
1 and r
2 are random numbers uniformly distributed between [
0
,
1
]. x
i
,
j is the position of any particle.


#### 2.3.2. Parameters Selection of SVM Based on PSO

In the SVM regression model, three parameters, namely,C, ε, and σ, should be identified before forecasting. Therefore, PSO algorithm is used for optimizing the SVM parameters. The process of optimizing the SVM parameters with PSO is presented in Figure 1 and the steps are described as follows [47].Figure 1
The process of optimizing the SVM parameters with PSO.Step  1. Initialization: consider randomly initial particles and velocities of a population (every particle contains three variables, namely, C, ε, and σ).Step  2. Fitness evaluation: the fitness function of PSO is shown as follows:
(12)
fitness
=
η
MAPE
=
1
m
∑
i
=
1
m
|
y
i
-
y
^
i
|
y
i
,
where η
MAPE is the mean absolute percentage error; y
i is the actual value and y
^
i is the predicted value; and m is the number of subsets. The solution with a smaller η
MAPE has a smaller fitness value.Step  3. Update global and personal best according to fitness evaluation results.Step  4. Calculation of velocity: particle flies toward a new position by calculating the velocity of position change. Velocity of each particle is calculated by (10).Step  5. Update position value: each particle moves to its next position according to (11).Step  6. Termination: repeat the same procedures from Step2 to Step5 until stopping conditions are satisfied.


### 2.4. Model Evaluation

It is essential to evaluate the performance of the models by employing appropriate methods. In this study, the performance of the models is evaluated by the indexes of the correlation coefficients(
R
), root mean squared error (RMSE), mean absolute error (MAE), and mean absolute relative error (MARE). These indexes are respectively defined as follows.Correlation coefficients(
R
):
(13)
R
=
∑
i
=
1
n
(
y
i
-
y
-
i
)
(
y
^
i
-
y
^
-
i
)
∑
i
=
1
n
(
y
i
-
y
-
i
)
2
∑
i
=
1
n
(
y
^
i
-
y
^
-
i
)
2
.Root mean squared error (RMSE):(14)
RMSE
=
1
n
∑
i
=
1
n
(
y
i
-
y
^
i
)
2
.Mean absolute error (MAE):(15)
MAE
=
1
n
∑
i
=
1
n
|
y
i
-
y
^
i
|
.Mean absolute relative error (MARE):(16)
MARE
=
1
n
∑
i
=
1
n
|
y
i
-
y
^
i
y
i
|
;
where y
i stands for the observed data and y
^
i stands for the forecasting data. n is number of the data.


## 2.1. Support Vector Machine (SVM)

The basic idea of SVM for regression is to introduce a kernel function, map the input data into a high-dimensional feature space by a nonlinear mapping, and then perform linear regression in the feature space [5]. Supposing that there is a training dataset D
=
{
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
…
,
(
x
n
,
y
n
)
}
∈
R
p
×
R, x is the input vector, y is the excepted output, n is the number of data, and p is the total number of data patterns. By nonlinear mapping function Ф, x is mapped into a feature space in which a linear estimate function is defined as
(1)
y
=
f
(
x
,
ω
)
=
〈
ω
,
Φ
(
x
)
〉
+
b
;
ω
,
x
∈
R
p
,
b
∈
R
,
where Ф
(
x
) represents the high-dimensional feature spaces, which is nonlinearly mapped from the input space x; ω and b are coefficients that have to be estimated from the input data. By introducing the slack variables ξ and ξ
* and following the regularization theory, parameters ω and b are estimated by minimizing the cost function
(2)
min
⁡
J
(
ω
,
ξ
)
=
1
2
∥
ω
∥
2
+
C
∑
i
=
1
n
(
ξ
i
+
ξ
i
*
)
subject to the constraints:
(3)
y
i
-
〈
ω
,
x
i
〉
-
b
≤
ε
+
ξ
i
〈
ω
,
x
i
〉
+
b
-
y
i
≤
ε
+
ξ
i
*
ξ
i
,
ξ
i
*
≥
0
.The first term(
1
/
2
)
∥
ω
∥
2 is weight vector norm; C is referred to as the regularized constraint determining the tradeoff between the empirical error and the regularized term; and ε is the insensitive loss function.By using Lagrange multiplier techniques, the minimization of (2) leads to the following dual-optimization problem:
(4)
max
⁡
W
(
α
i
,
α
i
*
)

=

-
1
2
∑
i
=
1
n
∑
j
=
1
n
[
(
α
i
-
α
i
*
)
(
α
j
-
α
j
*
)

K
(
x
i
,
x
j
)
]
11111111111111111
-
ε
∑
i
=
1
n
(
α
i
+
α
i
*
)
+
∑
i
=
1
n
y
i
(
α
i
-
α
i
*
)
s
.
t
.
∑
i
=
1
n
(
α
i
-
α
i
*
)
=
0
;
α
i
,
α
i
*
∈
[
0
,
C
]
,
where (
α
i
*
,
α
i
) are coefficients determined by training and K
(
x
i
,
x
j
) is the kernel function which can be expressed as inner product:
(5)
K
(
x
i
,
x
)
=
Φ
(
x
i
)
T
Φ
(
x
)
.The decision function takes the form(6)
f
(
x
)
=
∑
i
=
1
n
(
α
i
*
-
α
i
)
K
(
x
,
x
i
)
+
b
.The selection of an appropriate kernel function plays an important role in SVM regression since the kernel function defines the feature space. Gaussian radial basis function kernel has received significant attention from the machine learning community. Gaussian radial basis function (RBF) kernel is defined as(7)
K
(
x
,
x
i
)
=
exp
⁡
(
-
∥
x
-
x
i
∥
2
2
σ
2
)
.Hereσ is the kernel parameter.


## 2.2. Wavelet Analysis

In wavelet analysis, the signals are analyzed in both the time and the frequency domain by decomposing the original signals in different frequency bands using wavelet functions. The wavelet transform (WT) uses the scalable windowing technique for analyzing local variation in the time series [41]. WT provides useful decompositions of original time series, so that wavelet-transformed data improve the ability of a forecasting model by capturing useful information on various resolution levels [42]. The time series data are preprocessed using wavelet transformation techniques to obtain decomposed wavelet coefficients that are used as inputs in the forecasting models.The basic objective of WT is to achieve a complete timescale representation of localized and transient phenomena occurring at different timescales [4, 43]. The continuous wavelet transform is defined as the sum over all time of the signal multiplied by scale and shifted versions of wavelet function ψ:
(8)
W
(
a
,
b
)
=
1
a
∫
-
∞
+
∞
f
(
x
)
ψ
(
x
-
b
a
)
d
x
,
where a is a scale parameter; b is a position parameter; and ψ corresponds to the complex conjugate. The coefficient plots of the continuous wavelet transform are precisely the timescale view of the signal. However, calculating wavelet coefficients at every possible scale is time-consuming and generates large amount of information. Thus, the use of the continuous wavelet transform for forecasting is not practically possible.In hydrology, observed hydrologic series are often expressed as discrete series, so the discrete wavelet transform is usually employed to decompose a hydrologic series into a set of coefficients and subsignals under different scales, and then guide other time series analyses [44]. The DWT is defined as the following form:
(9)
f
(
t
)
=
∑
C
j
o
,
k
ϕ
j
o
,
k
(
t
)
+
∑
j
>
j
o
∑
ω
j
,
k
2
j
/
2
ψ
(
2
j
t
-
k
)
,
where j is the dilation or level index, k is the translation or scaling index, and ϕ
j
o
,
k is a scaling function of coarse scale coefficients. C
j
o
,
k, ω
j
,
k, is the scaling function of detailed (fine scale) coefficients and all functions of ψ
(
2
j
t
-
k
) are orthonormal.The original time series are decomposed into various details (D
s) and an approximation (A
s) at different resolution levels using DWT. The approximations are the high-scale, low frequency components of the signal and the details are the low-scale, high frequency components. Normally, the low frequency component of the signal (
A
) is the most important part which demonstrates the signal identity [45]. The choice of wavelet type is an important issue. The Daubechies wavelets are one of the widely used in wavelet family, which are written dbN, where db is the “surname” and N is the order of the wavelet [46]. Daubechies wavelets exhibit good tradeoff between parsimony and information richness [34], so in this study the Daubechies wavelets were employed as the mother wavelet to decompose the time series.


## 2.3. Parameters Selection of SVM Based on PSO

### 2.3.1. The Principle of PSO

PSO, deriving from the research for the movement of organisms in a bird flocking or fish schooling, performs searches using a population (called swarm) of individuals (called particles) that are updated from iteration to iteration [47, 48]. An equation (velocity update) controls the swarm in moving around the search space seeking the optimum state. In each iteration, the algorithm saves the local optimum and compares it with the global (best yet) optimum values. Definitely the criteria for being chosen as an optimum state depend on the fitness of the objective function. Candidate solutions (decision variables) of any particle calculate and remember its own fitness. The position of any particle accelerated towards the global best position by using (10) and (11) [49]. In any search step t, the ith particle is used to update its candidate solution’s current position x
i
,
j
(
t
) by using local best p
i
,
j
(
t
) and best p
g
,
j
(
t
) position achieved yet. Consider the following:
(10)
v
i
,
j
(
t
+
1
)
=
ω
v
i
,
j
(
t
)
+
c
1
r
1
[
p
i
,
j
-
x
i
,
j
(
t
)
]
+
c
2
r
2
[
p
g
,
j
-
x
i
,
j
(
t
)
]
(11)
x
i
,
j
(
t
+
1
)
=
x
i
,
j
(
t
)
+
v
i
,
j
(
t
+
1
)
,
j
=
1,2
,
…
,
d
,
where v
i
,
j is velocity measures for particles; ω is inertial weight controlling velocity direction; c
1 and c
2 are acceleration coefficients; r
1 and r
2 are random numbers uniformly distributed between [
0
,
1
]. x
i
,
j is the position of any particle.


### 2.3.2. Parameters Selection of SVM Based on PSO

In the SVM regression model, three parameters, namely,C, ε, and σ, should be identified before forecasting. Therefore, PSO algorithm is used for optimizing the SVM parameters. The process of optimizing the SVM parameters with PSO is presented in Figure 1 and the steps are described as follows [47].Figure 1
The process of optimizing the SVM parameters with PSO.Step  1. Initialization: consider randomly initial particles and velocities of a population (every particle contains three variables, namely, C, ε, and σ).Step  2. Fitness evaluation: the fitness function of PSO is shown as follows:
(12)
fitness
=
η
MAPE
=
1
m
∑
i
=
1
m
|
y
i
-
y
^
i
|
y
i
,
where η
MAPE is the mean absolute percentage error; y
i is the actual value and y
^
i is the predicted value; and m is the number of subsets. The solution with a smaller η
MAPE has a smaller fitness value.Step  3. Update global and personal best according to fitness evaluation results.Step  4. Calculation of velocity: particle flies toward a new position by calculating the velocity of position change. Velocity of each particle is calculated by (10).Step  5. Update position value: each particle moves to its next position according to (11).Step  6. Termination: repeat the same procedures from Step2 to Step5 until stopping conditions are satisfied.


## 2.3.1. The Principle of PSO

PSO, deriving from the research for the movement of organisms in a bird flocking or fish schooling, performs searches using a population (called swarm) of individuals (called particles) that are updated from iteration to iteration [47, 48]. An equation (velocity update) controls the swarm in moving around the search space seeking the optimum state. In each iteration, the algorithm saves the local optimum and compares it with the global (best yet) optimum values. Definitely the criteria for being chosen as an optimum state depend on the fitness of the objective function. Candidate solutions (decision variables) of any particle calculate and remember its own fitness. The position of any particle accelerated towards the global best position by using (10) and (11) [49]. In any search step t, the ith particle is used to update its candidate solution’s current position x
i
,
j
(
t
) by using local best p
i
,
j
(
t
) and best p
g
,
j
(
t
) position achieved yet. Consider the following:
(10)
v
i
,
j
(
t
+
1
)
=
ω
v
i
,
j
(
t
)
+
c
1
r
1
[
p
i
,
j
-
x
i
,
j
(
t
)
]
+
c
2
r
2
[
p
g
,
j
-
x
i
,
j
(
t
)
]
(11)
x
i
,
j
(
t
+
1
)
=
x
i
,
j
(
t
)
+
v
i
,
j
(
t
+
1
)
,
j
=
1,2
,
…
,
d
,
where v
i
,
j is velocity measures for particles; ω is inertial weight controlling velocity direction; c
1 and c
2 are acceleration coefficients; r
1 and r
2 are random numbers uniformly distributed between [
0
,
1
]. x
i
,
j is the position of any particle.


## 2.3.2. Parameters Selection of SVM Based on PSO

In the SVM regression model, three parameters, namely,C, ε, and σ, should be identified before forecasting. Therefore, PSO algorithm is used for optimizing the SVM parameters. The process of optimizing the SVM parameters with PSO is presented in Figure 1 and the steps are described as follows [47].Figure 1
The process of optimizing the SVM parameters with PSO.Step  1. Initialization: consider randomly initial particles and velocities of a population (every particle contains three variables, namely, C, ε, and σ).Step  2. Fitness evaluation: the fitness function of PSO is shown as follows:
(12)
fitness
=
η
MAPE
=
1
m
∑
i
=
1
m
|
y
i
-
y
^
i
|
y
i
,
where η
MAPE is the mean absolute percentage error; y
i is the actual value and y
^
i is the predicted value; and m is the number of subsets. The solution with a smaller η
MAPE has a smaller fitness value.Step  3. Update global and personal best according to fitness evaluation results.Step  4. Calculation of velocity: particle flies toward a new position by calculating the velocity of position change. Velocity of each particle is calculated by (10).Step  5. Update position value: each particle moves to its next position according to (11).Step  6. Termination: repeat the same procedures from Step2 to Step5 until stopping conditions are satisfied.


## 2.4. Model Evaluation

It is essential to evaluate the performance of the models by employing appropriate methods. In this study, the performance of the models is evaluated by the indexes of the correlation coefficients(
R
), root mean squared error (RMSE), mean absolute error (MAE), and mean absolute relative error (MARE). These indexes are respectively defined as follows.Correlation coefficients(
R
):
(13)
R
=
∑
i
=
1
n
(
y
i
-
y
-
i
)
(
y
^
i
-
y
^
-
i
)
∑
i
=
1
n
(
y
i
-
y
-
i
)
2
∑
i
=
1
n
(
y
^
i
-
y
^
-
i
)
2
.Root mean squared error (RMSE):(14)
RMSE
=
1
n
∑
i
=
1
n
(
y
i
-
y
^
i
)
2
.Mean absolute error (MAE):(15)
MAE
=
1
n
∑
i
=
1
n
|
y
i
-
y
^
i
|
.Mean absolute relative error (MARE):(16)
MARE
=
1
n
∑
i
=
1
n
|
y
i
-
y
^
i
y
i
|
;
where y
i stands for the observed data and y
^
i stands for the forecasting data. n is number of the data.


## 3. Study Area and Data Analysis

In this study, we examined the data obtained from the monthly streamflow of the Tangnaihai Hydrological station located at the upper catchment of the Yellow River in Qinghai Province of China. Tangnaihai Hydrological station is the upstream hydrological station of Longyang Gorge Reservoir which is the largest regulating reservoir in the upper catchment of the Yellow River, so there are few human impacts disordering streamflow regular. Location of Tangnaihai Hydrological station is shown in Figure2. The Yellow River catchment covers an area of 95,000 km2.Figure 2
Location map of Yellow River and Tangnaihai station.The monthly streamflow time series of Tangnaihai Hydrological station, consisting of 636 monthly records (January 1956 to December 2008), are used in this study. The dataset was split up into two parts: training and testing, where the first dataset consisting of 536 monthly records (January 1956 to August 2000) was used for training, while the final dataset contains 100 monthly records (September 2000 to December 2008). Training data were used exclusively for model development and testing data were used to measure the performance of the model on untrained data. The testing set was also used to evaluate the forecasting ability of the model and to compare the proposed model with others.


## 4. Results Analysis

### 4.1. Wavelet Decomposition of Streamflow Time Series

The WT-PSO-SVM model structure is shown in Figure3. For the SVM model inputs, the original time series are decomposed into subseries with an approximation (A
s) with low frequency and details (
D
1
,
D
2
,
…
,
D
s
) with high frequency by Daubechies DWT algorithm.Figure 3
The WA-LSSVM-PSO model structure.The optimal decomposition level of the streamflow time series in wavelet analysis plays an important role in preserving the information and reducing the distortion of the datasets [4]. The number of decomposition levels controls the streamflow approximation in the data. The general rule for the appropriate decomposition levels is that the largest levels should be shorter than the size of the testing data [50].In this case, the largest scales were chosen as three for the Tangnaihai station streamflow time series. Therefore, the flow data sets are decomposed into various details (D
s) and an approximation (
A
3
) at three resolution levels (21-22-23) using db3 DWT shown in Figure 4. The new decomposed subseries present variations of the original time series on different periods. MATLAB codes were developed using its library functions to perform wavelet decomposition of the time series data. The correlation coefficients between each D subtime series and original monthly streamflow time series are given in Table 1 for the Tangnaihai station. In the table, the D
t
-
1 and Q
t denote the D subtime series at time t
-
1 and measured streamflow at time t, respectively. These correlation values provide information for the determination of effective wavelet components on streamflow. It can be seen from Table 1 that the D
3 has the highest correlation (
R
=
0.259
) among D
s. The average correlation between A and Q
t is 0.363. According to the correlation analysis between D
s and the original current streamflow data, the effective component (D
3) is selected. Then, the new series obtained by adding the effective D
3 and approximation component are used as an input combination to the SVM model.Table 1
The correlation coefficients between subtime series and original streamflow data.
Discrete wavelet components
Correlations
D
t
-
1 and Q
t
D
t
-
2 and Q
t
D
t
-
3 and Q
t
Average
D
1
0.103
−0.003
0.105
0.070
D
2
−0.081
0.084
0.328
0.164
D
3
−0.023
0.295
0.459
0.259
A
0.211
0.390
0.488
0.363Figure 4
Decomposed wavelet subtime series components (Ds) of runoff data of Tangnaihai station.


### 4.2. Parameters Selection of SVM Based on PSO

In this study, RBF is employed as kernel function of SVM forecasting model, so three parameters, namely, balance parameterC, insensitive parameter ε, and kernel function, parameter σ should be selected. Some researchers have shown that different kernel functions have little impact on performance, but kernel function parameter σ is a key factor affecting performance of SVM. Among three parameters, σ precisely defines structure of highly dimensional space, so it controls complexity of ultimate solution; C determines complexity of model and punishment level of fitting deviation; ε indicates forecasting model’s expectation on estimating functions’ error of sample data, and the larger ε, the less support vector number and more sparse solution expression. But large ε can also reduce accuracy of SVM forecasting model.For monthly streamflow time seriesX
=
{
x
1
,
x
2
,
…
,
x
N
}, the flow data at time i
+
p is predicted based on the previous flow data. The general expression is shown as follows:
(17)
x
i
+
p
=
f
(
x
i
+
1
,
x
i
+
2
,
…
,
x
i
+
p
-
1
)
,
where f is a nonlinear function indicating relationship of monthly streamflow time series; x
i is streamflow data at time i,  i
=
1,2
,
…
,
N; and  p is the forecasting step (month), which is set as 3 in this paper. N
-
p monthly streamflow time series data sets are used for training and testing SVM forecasting model. The particle swarm optimization is employed to optimize the best parameters set (
σ
,
C
,
ε
) of SVM model [51].For the Tangnaihai station, three input combinations based on preceding monthly streamflows are evaluated to estimate current streamflow value. The input combinations evaluated in the study are as follows: (i)Q
t
-
1, Q
t
-
2, and Q
t
-
3; (ii) Q
t
-
1
(
A
), Q
t
-
2
(
A
), and Q
t
-
3
(
A
); (iii) Q
t
-
1
(
A
+
D
3
), Q
t
-
2
(
A
+
D
3
), and Q
t
-
3
(
A
+
D
3
). In all cases, the output is the discharge Q
t for the current month.In the training stage, firstly the parametersσ, C, and ε of SVM model are optimized by PSO, the validation error is measured by (12), and the adjusted parameters with minimum validation error are selected as the most appropriate parameters which are provided in Table 2. Then, the optimal parameters are utilized to train SVM and WSVM models.Table 2
Optimal parameters of SVM and WSVM models.
Model
Input variables
Output
Parameters
σ
C
ε
(i) SVM
(i)Q
t
-
1, Q
t
-
2, and Q
t
-
3
Q
t
0.90156
2.3149
0.01
(ii) WSVM1
(ii)Q
t
-
1
(
A
), Q
t
-
2
(
A
), and Q
t
-
3
(
A
)
0.01
46.8765
0.01
(iii) WSVM2
(iii)Q
t
-
1
(
A
+
D
3
), Q
t
-
2
(
A
+
D
3
), and Q
t
-
3
(
A
+
D
3
)
28.4876
6.0929
0.01Before the training process begins, data normalization is often performed. Runoff time series data was normalized in the range[
0
,
1
] by the following equation:
(18)
y
i
=
x
i
-
x
min
⁡
x
max
⁡
-
x
min
⁡
,
where y
i represents the normalized data, while x
i is the actual observation value and x
max
⁡, x
min
⁡, respectively, represent the maximum and minimum value among the actual observation values.


### 4.3. Streamflow Forecasting Based on WT-PSO-SVM

As shown in Table2, the optimal parameters for SVM (input model (i)) and WSVM (input models (ii) and (iii)) models are σ
1
=
0.90156, C
1
=
2.3149, and ε
1
=
0.01; σ
2
=
0.01, C
2
=
46.8765, and ε
2
=
0.01; and σ
3
=
28.4876, C
3
=
6.0929, and ε
3
=
0.01. The optimal parameters are used to examine the accuracy of the SVM and WSVM forecasting models with the testing data sets. Table 3 shows the performance results obtained in the training and testing periods of the SVM and WSVM models for Tangnaihai station.Table 3
Forecasting performance indexes of SVM and WSVM models in Tangnaihai station.
Model
Training
Testing
R
RMSE (m3/s)
MAE (m3/s)
MARE (%)
R
RMSE (m3/s)
MAE (m3/s)
MARE (%)
(i) SVM
0.853
287.315
196.524
14.33
0.768
317.035
247.480
17.81
(ii) WSVM1
0.548
306.554
217.883
18.24
0.613
331.420
256.368
23.49
(iii) WSVM2
0.808
201.784
141.730
7.46
0.806
243.268
173.20
11.52Results obtained from the three models for 3-month-advance flow forecasting at Tangnaihai station are presented in Table3. Performance of the three models was compared by evaluating indexes of R, RMSE, MAE, and MARE. It is observed that the forecasting accuracy of the WSVM2 model was much better than that of corresponding SVM and WSVM1 models. Comparing the forecasting results of SVM and WSVM1, the R, RMSE, MAE, and MARE in the testing period were 0.768–0.613, 317.035–331.420 (m3/s), 247.480–256.368 (m3/s), and 17.81–23.49 (%), respectively. It is observed that SVM is more superior than WSVM1 because some useful details (D
s) of original streamflow series were eliminated in the model WSVM1. Comparing the forecasting results of SVM and WSVM2, the values of R, RMSE, MAE, and MARE in the testing period were 0.768–0.806, 317.035–243.268, 247.480–173.20, and 17.81–11.52, respectively. It is obvious from Table 3 and Figure 5 that the WSVM2 performs better than the SVM model. These results indicated that D
3 was an effective component for the runoff series and the D
1 and D
2 were the noise that should be eliminated before streamflow forecasting. Wavelet transform is a necessary process of data preprocessing for improving predicting accuracy.Figure 5
Predicted and observed streamflow in the testing period by SVM, WSVM1, and WSVM2.


## 4.1. Wavelet Decomposition of Streamflow Time Series

The WT-PSO-SVM model structure is shown in Figure3. For the SVM model inputs, the original time series are decomposed into subseries with an approximation (A
s) with low frequency and details (
D
1
,
D
2
,
…
,
D
s
) with high frequency by Daubechies DWT algorithm.Figure 3
The WA-LSSVM-PSO model structure.The optimal decomposition level of the streamflow time series in wavelet analysis plays an important role in preserving the information and reducing the distortion of the datasets [4]. The number of decomposition levels controls the streamflow approximation in the data. The general rule for the appropriate decomposition levels is that the largest levels should be shorter than the size of the testing data [50].In this case, the largest scales were chosen as three for the Tangnaihai station streamflow time series. Therefore, the flow data sets are decomposed into various details (D
s) and an approximation (
A
3
) at three resolution levels (21-22-23) using db3 DWT shown in Figure 4. The new decomposed subseries present variations of the original time series on different periods. MATLAB codes were developed using its library functions to perform wavelet decomposition of the time series data. The correlation coefficients between each D subtime series and original monthly streamflow time series are given in Table 1 for the Tangnaihai station. In the table, the D
t
-
1 and Q
t denote the D subtime series at time t
-
1 and measured streamflow at time t, respectively. These correlation values provide information for the determination of effective wavelet components on streamflow. It can be seen from Table 1 that the D
3 has the highest correlation (
R
=
0.259
) among D
s. The average correlation between A and Q
t is 0.363. According to the correlation analysis between D
s and the original current streamflow data, the effective component (D
3) is selected. Then, the new series obtained by adding the effective D
3 and approximation component are used as an input combination to the SVM model.Table 1
The correlation coefficients between subtime series and original streamflow data.
Discrete wavelet components
Correlations
D
t
-
1 and Q
t
D
t
-
2 and Q
t
D
t
-
3 and Q
t
Average
D
1
0.103
−0.003
0.105
0.070
D
2
−0.081
0.084
0.328
0.164
D
3
−0.023
0.295
0.459
0.259
A
0.211
0.390
0.488
0.363Figure 4
Decomposed wavelet subtime series components (Ds) of runoff data of Tangnaihai station.


## 4.2. Parameters Selection of SVM Based on PSO

In this study, RBF is employed as kernel function of SVM forecasting model, so three parameters, namely, balance parameterC, insensitive parameter ε, and kernel function, parameter σ should be selected. Some researchers have shown that different kernel functions have little impact on performance, but kernel function parameter σ is a key factor affecting performance of SVM. Among three parameters, σ precisely defines structure of highly dimensional space, so it controls complexity of ultimate solution; C determines complexity of model and punishment level of fitting deviation; ε indicates forecasting model’s expectation on estimating functions’ error of sample data, and the larger ε, the less support vector number and more sparse solution expression. But large ε can also reduce accuracy of SVM forecasting model.For monthly streamflow time seriesX
=
{
x
1
,
x
2
,
…
,
x
N
}, the flow data at time i
+
p is predicted based on the previous flow data. The general expression is shown as follows:
(17)
x
i
+
p
=
f
(
x
i
+
1
,
x
i
+
2
,
…
,
x
i
+
p
-
1
)
,
where f is a nonlinear function indicating relationship of monthly streamflow time series; x
i is streamflow data at time i,  i
=
1,2
,
…
,
N; and  p is the forecasting step (month), which is set as 3 in this paper. N
-
p monthly streamflow time series data sets are used for training and testing SVM forecasting model. The particle swarm optimization is employed to optimize the best parameters set (
σ
,
C
,
ε
) of SVM model [51].For the Tangnaihai station, three input combinations based on preceding monthly streamflows are evaluated to estimate current streamflow value. The input combinations evaluated in the study are as follows: (i)Q
t
-
1, Q
t
-
2, and Q
t
-
3; (ii) Q
t
-
1
(
A
), Q
t
-
2
(
A
), and Q
t
-
3
(
A
); (iii) Q
t
-
1
(
A
+
D
3
), Q
t
-
2
(
A
+
D
3
), and Q
t
-
3
(
A
+
D
3
). In all cases, the output is the discharge Q
t for the current month.In the training stage, firstly the parametersσ, C, and ε of SVM model are optimized by PSO, the validation error is measured by (12), and the adjusted parameters with minimum validation error are selected as the most appropriate parameters which are provided in Table 2. Then, the optimal parameters are utilized to train SVM and WSVM models.Table 2
Optimal parameters of SVM and WSVM models.
Model
Input variables
Output
Parameters
σ
C
ε
(i) SVM
(i)Q
t
-
1, Q
t
-
2, and Q
t
-
3
Q
t
0.90156
2.3149
0.01
(ii) WSVM1
(ii)Q
t
-
1
(
A
), Q
t
-
2
(
A
), and Q
t
-
3
(
A
)
0.01
46.8765
0.01
(iii) WSVM2
(iii)Q
t
-
1
(
A
+
D
3
), Q
t
-
2
(
A
+
D
3
), and Q
t
-
3
(
A
+
D
3
)
28.4876
6.0929
0.01Before the training process begins, data normalization is often performed. Runoff time series data was normalized in the range[
0
,
1
] by the following equation:
(18)
y
i
=
x
i
-
x
min
⁡
x
max
⁡
-
x
min
⁡
,
where y
i represents the normalized data, while x
i is the actual observation value and x
max
⁡, x
min
⁡, respectively, represent the maximum and minimum value among the actual observation values.


## 4.3. Streamflow Forecasting Based on WT-PSO-SVM

As shown in Table2, the optimal parameters for SVM (input model (i)) and WSVM (input models (ii) and (iii)) models are σ
1
=
0.90156, C
1
=
2.3149, and ε
1
=
0.01; σ
2
=
0.01, C
2
=
46.8765, and ε
2
=
0.01; and σ
3
=
28.4876, C
3
=
6.0929, and ε
3
=
0.01. The optimal parameters are used to examine the accuracy of the SVM and WSVM forecasting models with the testing data sets. Table 3 shows the performance results obtained in the training and testing periods of the SVM and WSVM models for Tangnaihai station.Table 3
Forecasting performance indexes of SVM and WSVM models in Tangnaihai station.
Model
Training
Testing
R
RMSE (m3/s)
MAE (m3/s)
MARE (%)
R
RMSE (m3/s)
MAE (m3/s)
MARE (%)
(i) SVM
0.853
287.315
196.524
14.33
0.768
317.035
247.480
17.81
(ii) WSVM1
0.548
306.554
217.883
18.24
0.613
331.420
256.368
23.49
(iii) WSVM2
0.808
201.784
141.730
7.46
0.806
243.268
173.20
11.52Results obtained from the three models for 3-month-advance flow forecasting at Tangnaihai station are presented in Table3. Performance of the three models was compared by evaluating indexes of R, RMSE, MAE, and MARE. It is observed that the forecasting accuracy of the WSVM2 model was much better than that of corresponding SVM and WSVM1 models. Comparing the forecasting results of SVM and WSVM1, the R, RMSE, MAE, and MARE in the testing period were 0.768–0.613, 317.035–331.420 (m3/s), 247.480–256.368 (m3/s), and 17.81–23.49 (%), respectively. It is observed that SVM is more superior than WSVM1 because some useful details (D
s) of original streamflow series were eliminated in the model WSVM1. Comparing the forecasting results of SVM and WSVM2, the values of R, RMSE, MAE, and MARE in the testing period were 0.768–0.806, 317.035–243.268, 247.480–173.20, and 17.81–11.52, respectively. It is obvious from Table 3 and Figure 5 that the WSVM2 performs better than the SVM model. These results indicated that D
3 was an effective component for the runoff series and the D
1 and D
2 were the noise that should be eliminated before streamflow forecasting. Wavelet transform is a necessary process of data preprocessing for improving predicting accuracy.Figure 5
Predicted and observed streamflow in the testing period by SVM, WSVM1, and WSVM2.


## 5. Conclusion

This study developed a WT-PSO-SVM hybrid model to forecast monthly streamflow. The WT-PSO-SVM model was obtained by combining three methods, discrete wavelet transform-particle swarm optimization, and support vector machine regression. The combined model integrated the advantages of best versatility, robustness and effectiveness of SVM, the best global searching ability and the simple implementing procedure of PSO for parameter selection, and the ability of WT to reveal simultaneously both spectral and temporal information within one signal. This hybrid approach was successfully applied to simulate streamflow time series of Tangnaihai Hydrology station in the Yellow River.The streamflow time series were decomposed into various details (D
s) and an approximation (
A
3
) at three resolution levels (21-22-23) by using db3 DWT of the wavelet function of Daubechies 3 (db3). The correlation coefficients between each D subtime series and original monthly streamflow time series were calculated. D
s
(
D
3
) components with high correlation coefficients (
R
=
0.259
) were added to the approximation (
A
) as the input values of SVM model. The input combinations evaluated in the study are as follows: (i) Q
t
-
1, Q
t
-
2, and Q
t
-
3; (ii) Q
t
-
1
(
A
), Q
t
-
2
(
A
), and Q
t
-
3
(
A
); (iii) Q
t
-
1
(
A
+
D
3
), Q
t
-
2
(
A
+
D
3
), and Q
t
-
3
(
A
+
D
3
). The PSO was employed to select the optimal parameters, C, ε, and σ, of the three input models which were used to test the accuracy of the SVM model. Three different input combinations of SVM predicting results indicated that the discrete wavelet transform can significantly increase the accuracy of the SVM model in forecasting monthly streamflow. In addition, particle swarm optimization can determine suitable parameters to forecast streamflow as well. Predicting accuracy was evaluated by indexes of R, RMSE, MAE, and MARE. At the Tangnahai station, the best predictions belong to WSVM2 model. WSVM2 model increased the prediction R by 0.038 and 0.193 with respect to the SVM and WSVM1 models and reduced MARE by 6.29% and 11.97%, respectively, in the testing period. These results indicated that D
3 was an effective component for the runoff series and the D
1and D
2 were the noise that should be eliminated before streamflow forecasting. Wavelet transform is a necessary process of data preprocessing for improving predicting accuracy. The test results indicated that PSO-WT-SVM approach provides a superior alternative to the single SVM model for forecasting monthly streamflow in situations, without formulating models for the internal structure of the watershed.



---

*Source: 910196-2014-05-22.xml*