# Random Modeling of Daily Rainfall and Runoff Using a Seasonal Model and Wavelet Denoising

Authors: Chien-ming Chou
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Instead of Fourier smoothing, this study applied wavelet denoising to acquire the smooth seasonal mean and corresponding perturbation term from daily rainfall and runoff data in traditional seasonal models, which use seasonal means for hydrological time series forecasting. The denoised rainfall and runoff time series data were regarded as the smooth seasonal mean. The probability distribution of the percentage coefficients can be obtained from calibrated daily rainfall and runoff data. For validated daily rainfall and runoff data, percentage coefficients were randomly generated according to the probability distribution and the law of linear proportion. Multiplying the generated percentage coefficient by the smooth seasonal mean resulted in the corresponding perturbation term. Random modeling of daily rainfall and runoff can be obtained by adding the perturbation term to the smooth seasonal mean. To verify the accuracy of the proposed method, daily rainfall and runoff data for the Wu-Tu watershed were analyzed. The analytical results demonstrate that wavelet denoising enhances the precision of daily rainfall and runoff modeling of the seasonal model. In addition, the wavelet denoising technique proposed in this study can obtain the smooth seasonal mean of rainfall and runoff processes and is suitable for modeling actual daily rainfall and runoff processes.


---

## Body

## 1. Introduction

The modeling of rainfall and runoff processes is extremely important in the field of hydrology. Accurate and immediate hydrological forecasting (e.g., watershed immediate flood forecasting) is essential because it can be used as a reference for disaster decision-making. Long-term hydrological forecasting (e.g., annual runoff) can provide a critical reference for water resource planning. Seasonal models that use seasonal means to simulate and forecast hydrological time series are typically employed to model and forecast daily rainfall and runoff. However, estimated seasonal means generally produce oscillation because of noise when data are not collected for a long duration. The Fourier smoothing method is employed in the traditional seasonal model. Stationary hydrological time series data can be denoised by applying Fourier smoothing. However, because rainfall and runoff are nonstationary processes, the effectiveness of denoising using the Fourier smoothing method is limited. Usually, when Fourier smoothing is applied, data distortion is generated, and the true meaning of the data is lost [1].Other conventional denoising schemes include the Wiener filter and the Kalman filter, which are suitable for linear systems. However, the Kalman filter requires the establishment of state functions. In practice, hydrological systems are nonlinear and constructing appropriate state space models is difficult. Therefore, the Wiener filter and the Kalman filter have limited applications in hydrological modeling. Although the extended Kalman filter (EKF) can be employed for nonlinear dynamics, it becomes unstable when systemic nonlinearity is significant [2, 3]. Thus, the ensemble Kalman filter (EnKF), which is a Monte Carlo-based Kalman filter, was developed. The EnKF has recently gained popularity in the field of hydrology. This is partially because the increased computing power available can finally facilitate ensemble simulations and partially because the filter is easily implemented. However, this method does require the development of techniques for producing ensemble model simulations [4]. In principle, the EnKF is optimal only for Gaussian error statistics. Because EnKF only propagates the first two moments of error statistics, its effectiveness for significantly nonlinear uncertainty evolutions is limited [5].Wavelet techniques are effective for denoising; their numerous applications include image research [6]. Based on the Mallat algorithm of multiresolution analysis (MRA) of wavelet transformation, the original hydrological time series data can be decomposed into the approximated signals with low frequency and detailed signals with high frequency. Usually, useful signals appear as low frequency or approximate signals whereas noise usually appears as high frequency or detailed signals. For different scales of wavelet transformation, useful signals and noises have opposite properties [7]. The wavelet transformation coefficient of useful signals increases while that of noises decreases when the scale increases [7]. Thus, after several wavelet transformations, the wavelet transformation coefficient of noises has been eliminated or very small and the coefficient of useful signals remains [8]. The concept of wavelet denoising is to decompose noisy data using wavelet transform, threshold the detailed signals with high frequency, and finally obtain denoised data by wavelet reconstruction. The wavelet denoising algorithm can satisfy demands for various denoising procedures. Wavelet denoising is superior to conventional methods and has been applied recently in the field of hydrology. Lim and Lye [9] applied wavelet-based denoising to correct a series of high temporal resolution data for a streamflow after the data had been corrupted by tidal data. Their study confirmed that there is a tidal influence at the fluvial flow gauging site. Their method also demonstrated the potential use of wavelet analysis for solving similar problems. Liu et al. [10] applied wavelet analysis to decompose runoff time series data into approximate data and detailed data. Before reconstructing the wavelet method, some thresholds were added to various details so as to reduce noise in the original runoff data. The convergence capability, learning precision, and network generalization are greatly improved in a back-propagation (BP) network model with denoised runoff data.Wang and Fei [11] applied wavelet analysis to obtain the yearly periodic components in series data for hydrologic runoff. After eliminating periodic components, the remaining series data were denoised by wavelet to obtain the dependent stochastic components. An autoregression model was then constructed using the dependent stochastic series data. Their modeling results showed that, compared to traditional stochastic methods, simulation by wavelet method obtained hydrological series parameters that were closer to those of the measured series. Wang et al. [12] indicated that, during wavelet analysis of hydrologic series, denoising methods should be used to eliminate the effects of noise. The affected range of the hydrologic series data should be discarded before analysis, and the anomalous series data should be used to highlight the actual undulation of the hydrologic series. Cui et al. [13] applied the gray topological prediction method based on wavelet denoising in forecasting precipitation. Their computational results showed that their model was simpler and more accurate than the basic gray topological prediction model and therefore provided an important tool for forecasting precipitation and for preventing and mitigating disasters.Wang et al. [14] developed a new wavelet transform method for the synthetic generation of daily streamflow sequences. The advantage of their method is that the generated sequences can capture the dependence structure and statistical properties present in the data. Thus, wavelet denoising can potentially be combined with synthetic data generation. Chou [15] developed a novel framework for considering wavelet denoising in linear perturbation models (LPMs) and simple linear models (SLMs). Denoised rainfall and runoff time series data were applied to a SLM and used as the smooth seasonal mean employed by a LPM. Noise was employed as the perturbation term in the LPM. The denoised runoff and estimated runoff noise were summed to estimate the overall runoff in the LPM. The analytical results demonstrated that wavelet denoising enhances the rainfall-runoff modeling precision of the LPM.Instead of Fourier smoothing, this study applied wavelet denoising to acquire the smooth seasonal average and corresponding perturbation term from rainfall and runoff data in the seasonal model. The remainder of this paper is organized as follows. First, the wavelet denoising is described. The structures of the seasonal model are then defined. A case study of a small Taiwan watershed is presented to demonstrate the effectiveness of the proposed method. Finally, analytical results are discussed and conclusions are provided.


## 2. Materials and Methods

### 2.1. Wavelet Denoising

In practice, useful signals usually appear as low frequency or approximate signals whereas noise usually appears as high frequency or detailed signals. Wavelet transforms can separate low and high frequency signals. Signals or time series data can then be decomposed according to their specific characteristics. Wavelet reconstruction can derive denoised signals from processed decomposed signals and an approximate signal [16].


#### 2.1.1. Wavelet Transforms and Daubechies Wavelet Coefficients

The discrete wavelet transform of a vector is the outcome of a linear transformation that generates a new vector of dimension equal to that of the primeval vector. This transformation, also called decomposition, can also be performed efficiently by using the Mallet MRA algorithm [17]. Two discrete filters, decomposition low-pass filter  H and decomposition high-pass filter  G, are needed to compute discrete wavelet transforms. This work applied the DAUB4 filters proposed by Daubechies. The DAUB4 has the following four coefficients: c
0, c
1, c
2, and c
3. A decomposed matrix,   F, can be applied to decompose the hydrological data vector,   X, as follows [18]:
(1)
F
X

=

[
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
⋮
⋮
⋱
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
c
2
c
3
c
0
c
1
c
1
-
c
0
c
3
-
c
2
]
[
x
1
x
2
⋮
x
n
]
.Here, a blank space represents a zero value. Equation (1) gives the one resolution level decomposition. The term “resolution level” refers to the available decomposition level for the data. Let the filter {
c
0
,
c
1
,
c
2
,
c
3
} be smoothing filter  H. The odd elements of the  F
X output are obtained by convolving the  H with the X, which gives an approximation of the original data,   S. The filter {
c
3
,
-
c
2
,
c
1
,
-
c
0
}, denoted as  G, cannot be a smoothing filter due to its negative values. The even elements of the  F
X output are obtained by convolving the  G with the  X, which obtains  W, the detailed signal for the original data [18]. For a three-resolution-level decomposition, the original data (S
0) is decomposed into S
1 and W
1, and S
1 is then decomposed into S
2 and W
2. Finally, S
2 is decomposed into S
3 and W
3 (Figure 1(a)). In Figure 1(a), ↓ denotes the downsampling process, which retains alternate samples. The  H and  G filters associated with decomposition form the “analysis filter bank” [19].Filter banks structure for a three-level discrete wavelet transform. (a) Analysis bank and (b) synthesis bank [19].
(a)
(b)To apply the features of the approximate and detailed signals, data with lengthN should be restructured from the approximation with length N
/
2 and from the detailed signal with length N
/
2. The transform matrix should be orthogonal, meaning that the inverse matrix of the decomposed matrix should be equal to the transposed matrix of the decomposed matrix. The resulting reconstruction matrix   F
′ is obtained as follows [18]:
(2)
[
c
0
c
3
…
c
2
c
1
c
1
-
c
2
…
c
3
-
c
0
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
-
c
2
⋱
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
c
2
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
-
c
2
]
.The result obtained by carrying out 1-level wavelet decomposition and then carrying out reconstruction is equal to original data; that is,F
′
(
F
X
)
=
X. Three resolution levels are reconstructed as shown in Figure 1(b) where ↑ denotes the upsampling process and zeros are inserted between adjacent samples. The H
′ and   G
′ filters that effect reconstruction form the “synthesis filter bank” [19].The coefficient values in the filter{
c
0
,
c
1
,
c
2
,
c
3
} estimated by Daubechies are as follows [18]:
(3)
c
0
=
1
+
3
4
2
,
c
1
=
3
+
3
4
2
,
c
2
=
3
-
3
4
2
,
c
3
=
1
-
3
4
2
.


#### 2.1.2. Wavelet Denoising Procedure

Donoho [20] proposed a wavelet denoising method based on thresholds for acquiring correct denoised results. This method, which is now the most common method of wavelet denoising, is performed as follows [16].(i)
Choose an appropriate wavelet function and number of resolution levelM. The original one-dimensional time series S
0 is decomposed into an approximation at resolution level M and detailed signals at various resolution levels by using wavelet transform.
(ii)
Below certain thresholds, the absolute values of detailed signals,w
j
(
t
)
(
j
=
1
,
2
,
…
,
M
), are set to zero at each resolution level. The subscript j represents the jth resolution levels. The absolute values of detailed signals that exceed certain thresholds are treated as the difference between the values of detailed signals and thresholds as follows [16]:
(4)
w
^
j
(
t
)
=
{
sgn
⁡
(
w
j
(
t
)
)
(
|
w
j
(
t
)
|
-
T
)
|
w
j
(
t
)
|
>
T
0
|
w
j
(
t
)
|
≤
T
.
Equation (4) gives the threshold quantifications used to obtain the processed detailed signals at each resolution level during wavelet denoising. The approximation usually does not perform threshold quantifications.
(iii)
Wavelet reconstruction can derive the denoised time series data from the approximation at resolution levelM and processed detailed signals (w
^
j
(
t
)) at all resolution levels.


#### 2.1.3. Determining Thresholds

In this work, the constant thresholdT at each resolution level was determined as follows [16]:
(5)
T
=
σ
2
ln
⁡
n
j
,
where n
j is the length of detailed signals at each resolution level j and σ is the noise strength of detailed signals at each resolution level. Noise strength was obtained as follows [16]:
(6)
σ
=
1
0.6745
n
j
∑
t
=
1
n
j
|
w
j
(
t
)
|
,
where w
j
(
t
) are the detailed signals at resolution level j. The T varies according to the length of detailed signals and noise strength.


### 2.2. Random Rainfall and Runoff Modeling Using a Seasonal Model

#### 2.2.1. The Seasonal Model

The primary concept of a seasonal model is that the hydrological time series are divided into seasonal means and perturbation terms. The seasonal mean has the strongest effect on the behavior of the original time series. Seasonal means are essential for assessing the hydrological properties of a river basin, which should be a smooth process. Additionally, seasonal means can be applied to forecast the hydrological time series. However, estimated seasonal means generally produce oscillation because of the noise when data are not collected for a long duration. The seasonal mean should be smooth in the seasonal model.The seasonal model procedure used in this study is explained below [21].(i)
A considerable amount of useful information, such as the seasonal mean, is hidden in daily rainfall or runoff time series. The seasonal mean of the rainfall or runoff time seriesI
d can be obtained from observed data [21]:
(7)
I
d
=
1
n
(
I
d
,
1
+
I
d
,
2
+
⋯
+
I
d
,
i
+
⋯
+
I
d
,
n
)
,
where I
d is the seasonal mean of rainfall or runoff, I
d
,
i is the rainfall or runoff on day d and year i, and n  is the years of when the data were obtained.
(ii)
The estimated seasonal mean typically produces vibrations with noise when the time span for the data is inadequate. Thus, an appropriate mathematical method must be adopted to smooth seasonal meansI
d into I
~
d. The method most commonly used in seasonal models is the Fourier series presented below [21]:
(8)
I
~
d
=
I
-
d
+
∑
j
=
1
L
[
A
I
,
j
cos
⁡
⁡
(
2
π
j
d
365
)
+
B
I
,
j
sin
⁡
(
2
π
j
d
365
)
]
,
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
(
d
=
1,2
,
…
,
365
)
,
where I
d
- is the average, A
I
,
j and B
I
,
j are Fourier coefficients, and j is the order of the harmonic function; that is [21],
(9)
I
-
d
=
1
365
∑
d
=
1
365
[
I
d
]
(10)
A
I
,
j
=
2
365
∑
d
=
1
L
[
I
d
cos
⁡
(
2
π
j
d
365
)
]
(11)
B
I
,
j
=
2
365
∑
d
=
1
L
[
I
d
sin
(
2
π
j
d
365
)
]
.
When the harmonic function in (8) has limited terms, the smooth seasonal mean I
~
d can be obtained. The L is typically set to 4 or 5.
(iii)
This study assumes thatI
(
t
) represents the observed rainfall or runoff data. The perturbation term P
(
t
) can be calculated using [21]:
(12)
P
(
t
)
=
I
(
t
)
-
I
~
d
.


#### 2.2.2. Random Generation of Percentage Coefficients

The seasonal mean is an essential attribute of river basins. Under ideal conditions, actual daily rainfall or runoff should be extremely similar to the daily seasonal mean of the basin, with obvious seasonal variation rules. In practice, observed daily rainfall or runoff fluctuates against the seasonal mean because of disturbance from external factors. However, actual observed daily rainfall or runoff is typically similar to the seasonal mean. The relationship between the fluctuation size and the occurrence frequency is standard. The ratio of the fluctuation size (i.e., perturbation term) to the seasonal mean is defined as the percentage coefficientsk, as shown below [21]:
(13)
k
=
P
(
t
)
I
~
d
.The probability distribution of the percentage coefficientsk can be obtained from calibrated daily rainfall or runoff data. Generally, the range of percentage coefficients k is divided into several intervals. The cumulative probability in each interval can be calculated to determine the probability distribution of the percentage coefficients k. To validate the daily rainfall or runoff data, we generated a series of random numbers r in the interval [
0,1
], which is similar to the concept of cumulative probability. Based on each r value (i.e., cumulative probability) and the law of linear proportion, we can determine a new corresponding interval for each r value. After specifying the new interval, a random number p can be generated between the lower and upper limits of the specific interval, which represents the estimated percentage coefficients k
^. The estimated percentage coefficients k
^ multiplied by the smooth seasonal mean I
~
d equals the estimated perturbation term P
^
(
t
). The estimated rainfall or runoff I
^
(
t
) can be obtained from the sum of the estimated perturbation term P
^
(
t
) and the smooth seasonal mean I
~
d; that is, I
^
(
t
)
=
P
^
(
t
)
+
I
~
d.


### 2.3. The Proposed Wavelet-Based Method

Instead of employing a Fourier series, this study applied wavelet denoising to determine the smooth seasonal mean and perturbation term from rainfall and runoff data in the seasonal model. The procedures of the proposed wavelet-based method are shown below.(i)
First, rainfall and runoff time series data were decomposed into approximations and detailed signals at various resolution levels using the DAUB4 wavelet transform or (1) and (3).
(ii)
Threshold quantifications were performed at each resolution level using (4)–(6).
(iii)
The denoised rainfall and runoff time series data were obtained from the approximation at the final resolution level and the processed detailed signals by threshold quantifications at all rainfall and runoff resolution levels using wavelet reconstruction or (2) and (3).
(iv)
The denoised rainfall and runoff time series data were considered the smooth seasonal means employed in the seasonal model using (8)–(11).
(v)
The value of the original time series minus the denoised time series using (12) was considered the perturbation term P
(
t
) employed in the seasonal model.
(vi)
Based on calibrated rainfall and runoff data, the probability distribution of the percentage coefficients can be determined.
(vii)
After determining the probability distribution of the percentage coefficients, the estimated percentage coefficientsk
^ in the validated rainfall and runoff data can be obtained using (13) through random generation in a specific range of the probability distribution.
(viii)
The estimated percentage coefficientsk
^ multiplied by the smooth seasonal mean I
~
d using wavelet denoising equals the estimated perturbation term P
^
(
t
).
(ix)
The estimated rainfall or runoffI
^
(
t
) can be obtained from the sum of the estimated perturbation term P
^
(
t
) and the smooth seasonal mean I
~
d; that is, I
^
(
t
)
=
P
^
(
t
)
+
I
~
d.


### 2.4. Application and Analysis

#### 2.4.1. Study Basin

This study confirmed the efficacy of the proposed wavelet-based method for modeling daily rainfall and runoff using the Wu-Tu watershed in northern Taiwan as an example. The upstream area of the watershed is 203 km2 (Figure 2), and the mean annual precipitation in this watershed is 2500 mm. Because of the topography of this watershed, the runoff pathlines are short and steep, and the rainfall is nonuniform in both time and space. Large floods that develop rapidly in the middle to downstream reaches of this watershed can cause severe damage. After collecting daily rainfall and runoff data from 1966 to 1994, data from 1966 to 1980 were used to calibrate the proposed model. Data from 1981 to 1994 were then used to assess the performance of the proposed method. Observed daily rainfall data were obtained using the kriging method to average the values obtained at the Jui-Fang, Huo-Shao-Liao, and Wu-Tu weather stations.Figure 2
A map of the Wu-Tu watershed study area located near Taipei, Taiwan.


#### 2.4.2. Comparison of Model Performances

To quantitatively compare Fourier smoothing with wavelet denoising in a seasonal model, the validated results were evaluated based on the mean square error (MSE). The MSE is defined as(14)
MSE
=
∑
i
=
1
N
[
q
(
i
)
-
q
^
(
i
)
]
2
N
,
where q
^
(
i
) denotes the simulated discharge of the hydrograph for time period i
(
m
3
/
s
), q
(
i
) is the discharge of the observed hydrograph for time period i
(
m
3
/
s
), and N is number of data.


## 2.1. Wavelet Denoising

In practice, useful signals usually appear as low frequency or approximate signals whereas noise usually appears as high frequency or detailed signals. Wavelet transforms can separate low and high frequency signals. Signals or time series data can then be decomposed according to their specific characteristics. Wavelet reconstruction can derive denoised signals from processed decomposed signals and an approximate signal [16].


### 2.1.1. Wavelet Transforms and Daubechies Wavelet Coefficients

The discrete wavelet transform of a vector is the outcome of a linear transformation that generates a new vector of dimension equal to that of the primeval vector. This transformation, also called decomposition, can also be performed efficiently by using the Mallet MRA algorithm [17]. Two discrete filters, decomposition low-pass filter  H and decomposition high-pass filter  G, are needed to compute discrete wavelet transforms. This work applied the DAUB4 filters proposed by Daubechies. The DAUB4 has the following four coefficients: c
0, c
1, c
2, and c
3. A decomposed matrix,   F, can be applied to decompose the hydrological data vector,   X, as follows [18]:
(1)
F
X

=

[
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
⋮
⋮
⋱
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
c
2
c
3
c
0
c
1
c
1
-
c
0
c
3
-
c
2
]
[
x
1
x
2
⋮
x
n
]
.Here, a blank space represents a zero value. Equation (1) gives the one resolution level decomposition. The term “resolution level” refers to the available decomposition level for the data. Let the filter {
c
0
,
c
1
,
c
2
,
c
3
} be smoothing filter  H. The odd elements of the  F
X output are obtained by convolving the  H with the X, which gives an approximation of the original data,   S. The filter {
c
3
,
-
c
2
,
c
1
,
-
c
0
}, denoted as  G, cannot be a smoothing filter due to its negative values. The even elements of the  F
X output are obtained by convolving the  G with the  X, which obtains  W, the detailed signal for the original data [18]. For a three-resolution-level decomposition, the original data (S
0) is decomposed into S
1 and W
1, and S
1 is then decomposed into S
2 and W
2. Finally, S
2 is decomposed into S
3 and W
3 (Figure 1(a)). In Figure 1(a), ↓ denotes the downsampling process, which retains alternate samples. The  H and  G filters associated with decomposition form the “analysis filter bank” [19].Filter banks structure for a three-level discrete wavelet transform. (a) Analysis bank and (b) synthesis bank [19].
(a)
(b)To apply the features of the approximate and detailed signals, data with lengthN should be restructured from the approximation with length N
/
2 and from the detailed signal with length N
/
2. The transform matrix should be orthogonal, meaning that the inverse matrix of the decomposed matrix should be equal to the transposed matrix of the decomposed matrix. The resulting reconstruction matrix   F
′ is obtained as follows [18]:
(2)
[
c
0
c
3
…
c
2
c
1
c
1
-
c
2
…
c
3
-
c
0
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
-
c
2
⋱
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
c
2
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
-
c
2
]
.The result obtained by carrying out 1-level wavelet decomposition and then carrying out reconstruction is equal to original data; that is,F
′
(
F
X
)
=
X. Three resolution levels are reconstructed as shown in Figure 1(b) where ↑ denotes the upsampling process and zeros are inserted between adjacent samples. The H
′ and   G
′ filters that effect reconstruction form the “synthesis filter bank” [19].The coefficient values in the filter{
c
0
,
c
1
,
c
2
,
c
3
} estimated by Daubechies are as follows [18]:
(3)
c
0
=
1
+
3
4
2
,
c
1
=
3
+
3
4
2
,
c
2
=
3
-
3
4
2
,
c
3
=
1
-
3
4
2
.


### 2.1.2. Wavelet Denoising Procedure

Donoho [20] proposed a wavelet denoising method based on thresholds for acquiring correct denoised results. This method, which is now the most common method of wavelet denoising, is performed as follows [16].(i)
Choose an appropriate wavelet function and number of resolution levelM. The original one-dimensional time series S
0 is decomposed into an approximation at resolution level M and detailed signals at various resolution levels by using wavelet transform.
(ii)
Below certain thresholds, the absolute values of detailed signals,w
j
(
t
)
(
j
=
1
,
2
,
…
,
M
), are set to zero at each resolution level. The subscript j represents the jth resolution levels. The absolute values of detailed signals that exceed certain thresholds are treated as the difference between the values of detailed signals and thresholds as follows [16]:
(4)
w
^
j
(
t
)
=
{
sgn
⁡
(
w
j
(
t
)
)
(
|
w
j
(
t
)
|
-
T
)
|
w
j
(
t
)
|
>
T
0
|
w
j
(
t
)
|
≤
T
.
Equation (4) gives the threshold quantifications used to obtain the processed detailed signals at each resolution level during wavelet denoising. The approximation usually does not perform threshold quantifications.
(iii)
Wavelet reconstruction can derive the denoised time series data from the approximation at resolution levelM and processed detailed signals (w
^
j
(
t
)) at all resolution levels.


### 2.1.3. Determining Thresholds

In this work, the constant thresholdT at each resolution level was determined as follows [16]:
(5)
T
=
σ
2
ln
⁡
n
j
,
where n
j is the length of detailed signals at each resolution level j and σ is the noise strength of detailed signals at each resolution level. Noise strength was obtained as follows [16]:
(6)
σ
=
1
0.6745
n
j
∑
t
=
1
n
j
|
w
j
(
t
)
|
,
where w
j
(
t
) are the detailed signals at resolution level j. The T varies according to the length of detailed signals and noise strength.


## 2.1.1. Wavelet Transforms and Daubechies Wavelet Coefficients

The discrete wavelet transform of a vector is the outcome of a linear transformation that generates a new vector of dimension equal to that of the primeval vector. This transformation, also called decomposition, can also be performed efficiently by using the Mallet MRA algorithm [17]. Two discrete filters, decomposition low-pass filter  H and decomposition high-pass filter  G, are needed to compute discrete wavelet transforms. This work applied the DAUB4 filters proposed by Daubechies. The DAUB4 has the following four coefficients: c
0, c
1, c
2, and c
3. A decomposed matrix,   F, can be applied to decompose the hydrological data vector,   X, as follows [18]:
(1)
F
X

=

[
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
⋮
⋮
⋱
c
0
c
1
c
2
c
3
c
3
-
c
2
c
1
-
c
0
c
2
c
3
c
0
c
1
c
1
-
c
0
c
3
-
c
2
]
[
x
1
x
2
⋮
x
n
]
.Here, a blank space represents a zero value. Equation (1) gives the one resolution level decomposition. The term “resolution level” refers to the available decomposition level for the data. Let the filter {
c
0
,
c
1
,
c
2
,
c
3
} be smoothing filter  H. The odd elements of the  F
X output are obtained by convolving the  H with the X, which gives an approximation of the original data,   S. The filter {
c
3
,
-
c
2
,
c
1
,
-
c
0
}, denoted as  G, cannot be a smoothing filter due to its negative values. The even elements of the  F
X output are obtained by convolving the  G with the  X, which obtains  W, the detailed signal for the original data [18]. For a three-resolution-level decomposition, the original data (S
0) is decomposed into S
1 and W
1, and S
1 is then decomposed into S
2 and W
2. Finally, S
2 is decomposed into S
3 and W
3 (Figure 1(a)). In Figure 1(a), ↓ denotes the downsampling process, which retains alternate samples. The  H and  G filters associated with decomposition form the “analysis filter bank” [19].Filter banks structure for a three-level discrete wavelet transform. (a) Analysis bank and (b) synthesis bank [19].
(a)
(b)To apply the features of the approximate and detailed signals, data with lengthN should be restructured from the approximation with length N
/
2 and from the detailed signal with length N
/
2. The transform matrix should be orthogonal, meaning that the inverse matrix of the decomposed matrix should be equal to the transposed matrix of the decomposed matrix. The resulting reconstruction matrix   F
′ is obtained as follows [18]:
(2)
[
c
0
c
3
…
c
2
c
1
c
1
-
c
2
…
c
3
-
c
0
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
-
c
2
⋱
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
c
2
c
2
c
1
c
0
c
3
c
3
-
c
0
c
1
-
c
2
]
.The result obtained by carrying out 1-level wavelet decomposition and then carrying out reconstruction is equal to original data; that is,F
′
(
F
X
)
=
X. Three resolution levels are reconstructed as shown in Figure 1(b) where ↑ denotes the upsampling process and zeros are inserted between adjacent samples. The H
′ and   G
′ filters that effect reconstruction form the “synthesis filter bank” [19].The coefficient values in the filter{
c
0
,
c
1
,
c
2
,
c
3
} estimated by Daubechies are as follows [18]:
(3)
c
0
=
1
+
3
4
2
,
c
1
=
3
+
3
4
2
,
c
2
=
3
-
3
4
2
,
c
3
=
1
-
3
4
2
.


## 2.1.2. Wavelet Denoising Procedure

Donoho [20] proposed a wavelet denoising method based on thresholds for acquiring correct denoised results. This method, which is now the most common method of wavelet denoising, is performed as follows [16].(i)
Choose an appropriate wavelet function and number of resolution levelM. The original one-dimensional time series S
0 is decomposed into an approximation at resolution level M and detailed signals at various resolution levels by using wavelet transform.
(ii)
Below certain thresholds, the absolute values of detailed signals,w
j
(
t
)
(
j
=
1
,
2
,
…
,
M
), are set to zero at each resolution level. The subscript j represents the jth resolution levels. The absolute values of detailed signals that exceed certain thresholds are treated as the difference between the values of detailed signals and thresholds as follows [16]:
(4)
w
^
j
(
t
)
=
{
sgn
⁡
(
w
j
(
t
)
)
(
|
w
j
(
t
)
|
-
T
)
|
w
j
(
t
)
|
>
T
0
|
w
j
(
t
)
|
≤
T
.
Equation (4) gives the threshold quantifications used to obtain the processed detailed signals at each resolution level during wavelet denoising. The approximation usually does not perform threshold quantifications.
(iii)
Wavelet reconstruction can derive the denoised time series data from the approximation at resolution levelM and processed detailed signals (w
^
j
(
t
)) at all resolution levels.


## 2.1.3. Determining Thresholds

In this work, the constant thresholdT at each resolution level was determined as follows [16]:
(5)
T
=
σ
2
ln
⁡
n
j
,
where n
j is the length of detailed signals at each resolution level j and σ is the noise strength of detailed signals at each resolution level. Noise strength was obtained as follows [16]:
(6)
σ
=
1
0.6745
n
j
∑
t
=
1
n
j
|
w
j
(
t
)
|
,
where w
j
(
t
) are the detailed signals at resolution level j. The T varies according to the length of detailed signals and noise strength.


## 2.2. Random Rainfall and Runoff Modeling Using a Seasonal Model

### 2.2.1. The Seasonal Model

The primary concept of a seasonal model is that the hydrological time series are divided into seasonal means and perturbation terms. The seasonal mean has the strongest effect on the behavior of the original time series. Seasonal means are essential for assessing the hydrological properties of a river basin, which should be a smooth process. Additionally, seasonal means can be applied to forecast the hydrological time series. However, estimated seasonal means generally produce oscillation because of the noise when data are not collected for a long duration. The seasonal mean should be smooth in the seasonal model.The seasonal model procedure used in this study is explained below [21].(i)
A considerable amount of useful information, such as the seasonal mean, is hidden in daily rainfall or runoff time series. The seasonal mean of the rainfall or runoff time seriesI
d can be obtained from observed data [21]:
(7)
I
d
=
1
n
(
I
d
,
1
+
I
d
,
2
+
⋯
+
I
d
,
i
+
⋯
+
I
d
,
n
)
,
where I
d is the seasonal mean of rainfall or runoff, I
d
,
i is the rainfall or runoff on day d and year i, and n  is the years of when the data were obtained.
(ii)
The estimated seasonal mean typically produces vibrations with noise when the time span for the data is inadequate. Thus, an appropriate mathematical method must be adopted to smooth seasonal meansI
d into I
~
d. The method most commonly used in seasonal models is the Fourier series presented below [21]:
(8)
I
~
d
=
I
-
d
+
∑
j
=
1
L
[
A
I
,
j
cos
⁡
⁡
(
2
π
j
d
365
)
+
B
I
,
j
sin
⁡
(
2
π
j
d
365
)
]
,
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
(
d
=
1,2
,
…
,
365
)
,
where I
d
- is the average, A
I
,
j and B
I
,
j are Fourier coefficients, and j is the order of the harmonic function; that is [21],
(9)
I
-
d
=
1
365
∑
d
=
1
365
[
I
d
]
(10)
A
I
,
j
=
2
365
∑
d
=
1
L
[
I
d
cos
⁡
(
2
π
j
d
365
)
]
(11)
B
I
,
j
=
2
365
∑
d
=
1
L
[
I
d
sin
(
2
π
j
d
365
)
]
.
When the harmonic function in (8) has limited terms, the smooth seasonal mean I
~
d can be obtained. The L is typically set to 4 or 5.
(iii)
This study assumes thatI
(
t
) represents the observed rainfall or runoff data. The perturbation term P
(
t
) can be calculated using [21]:
(12)
P
(
t
)
=
I
(
t
)
-
I
~
d
.


### 2.2.2. Random Generation of Percentage Coefficients

The seasonal mean is an essential attribute of river basins. Under ideal conditions, actual daily rainfall or runoff should be extremely similar to the daily seasonal mean of the basin, with obvious seasonal variation rules. In practice, observed daily rainfall or runoff fluctuates against the seasonal mean because of disturbance from external factors. However, actual observed daily rainfall or runoff is typically similar to the seasonal mean. The relationship between the fluctuation size and the occurrence frequency is standard. The ratio of the fluctuation size (i.e., perturbation term) to the seasonal mean is defined as the percentage coefficientsk, as shown below [21]:
(13)
k
=
P
(
t
)
I
~
d
.The probability distribution of the percentage coefficientsk can be obtained from calibrated daily rainfall or runoff data. Generally, the range of percentage coefficients k is divided into several intervals. The cumulative probability in each interval can be calculated to determine the probability distribution of the percentage coefficients k. To validate the daily rainfall or runoff data, we generated a series of random numbers r in the interval [
0,1
], which is similar to the concept of cumulative probability. Based on each r value (i.e., cumulative probability) and the law of linear proportion, we can determine a new corresponding interval for each r value. After specifying the new interval, a random number p can be generated between the lower and upper limits of the specific interval, which represents the estimated percentage coefficients k
^. The estimated percentage coefficients k
^ multiplied by the smooth seasonal mean I
~
d equals the estimated perturbation term P
^
(
t
). The estimated rainfall or runoff I
^
(
t
) can be obtained from the sum of the estimated perturbation term P
^
(
t
) and the smooth seasonal mean I
~
d; that is, I
^
(
t
)
=
P
^
(
t
)
+
I
~
d.


## 2.2.1. The Seasonal Model

The primary concept of a seasonal model is that the hydrological time series are divided into seasonal means and perturbation terms. The seasonal mean has the strongest effect on the behavior of the original time series. Seasonal means are essential for assessing the hydrological properties of a river basin, which should be a smooth process. Additionally, seasonal means can be applied to forecast the hydrological time series. However, estimated seasonal means generally produce oscillation because of the noise when data are not collected for a long duration. The seasonal mean should be smooth in the seasonal model.The seasonal model procedure used in this study is explained below [21].(i)
A considerable amount of useful information, such as the seasonal mean, is hidden in daily rainfall or runoff time series. The seasonal mean of the rainfall or runoff time seriesI
d can be obtained from observed data [21]:
(7)
I
d
=
1
n
(
I
d
,
1
+
I
d
,
2
+
⋯
+
I
d
,
i
+
⋯
+
I
d
,
n
)
,
where I
d is the seasonal mean of rainfall or runoff, I
d
,
i is the rainfall or runoff on day d and year i, and n  is the years of when the data were obtained.
(ii)
The estimated seasonal mean typically produces vibrations with noise when the time span for the data is inadequate. Thus, an appropriate mathematical method must be adopted to smooth seasonal meansI
d into I
~
d. The method most commonly used in seasonal models is the Fourier series presented below [21]:
(8)
I
~
d
=
I
-
d
+
∑
j
=
1
L
[
A
I
,
j
cos
⁡
⁡
(
2
π
j
d
365
)
+
B
I
,
j
sin
⁡
(
2
π
j
d
365
)
]
,
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
(
d
=
1,2
,
…
,
365
)
,
where I
d
- is the average, A
I
,
j and B
I
,
j are Fourier coefficients, and j is the order of the harmonic function; that is [21],
(9)
I
-
d
=
1
365
∑
d
=
1
365
[
I
d
]
(10)
A
I
,
j
=
2
365
∑
d
=
1
L
[
I
d
cos
⁡
(
2
π
j
d
365
)
]
(11)
B
I
,
j
=
2
365
∑
d
=
1
L
[
I
d
sin
(
2
π
j
d
365
)
]
.
When the harmonic function in (8) has limited terms, the smooth seasonal mean I
~
d can be obtained. The L is typically set to 4 or 5.
(iii)
This study assumes thatI
(
t
) represents the observed rainfall or runoff data. The perturbation term P
(
t
) can be calculated using [21]:
(12)
P
(
t
)
=
I
(
t
)
-
I
~
d
.


## 2.2.2. Random Generation of Percentage Coefficients

The seasonal mean is an essential attribute of river basins. Under ideal conditions, actual daily rainfall or runoff should be extremely similar to the daily seasonal mean of the basin, with obvious seasonal variation rules. In practice, observed daily rainfall or runoff fluctuates against the seasonal mean because of disturbance from external factors. However, actual observed daily rainfall or runoff is typically similar to the seasonal mean. The relationship between the fluctuation size and the occurrence frequency is standard. The ratio of the fluctuation size (i.e., perturbation term) to the seasonal mean is defined as the percentage coefficientsk, as shown below [21]:
(13)
k
=
P
(
t
)
I
~
d
.The probability distribution of the percentage coefficientsk can be obtained from calibrated daily rainfall or runoff data. Generally, the range of percentage coefficients k is divided into several intervals. The cumulative probability in each interval can be calculated to determine the probability distribution of the percentage coefficients k. To validate the daily rainfall or runoff data, we generated a series of random numbers r in the interval [
0,1
], which is similar to the concept of cumulative probability. Based on each r value (i.e., cumulative probability) and the law of linear proportion, we can determine a new corresponding interval for each r value. After specifying the new interval, a random number p can be generated between the lower and upper limits of the specific interval, which represents the estimated percentage coefficients k
^. The estimated percentage coefficients k
^ multiplied by the smooth seasonal mean I
~
d equals the estimated perturbation term P
^
(
t
). The estimated rainfall or runoff I
^
(
t
) can be obtained from the sum of the estimated perturbation term P
^
(
t
) and the smooth seasonal mean I
~
d; that is, I
^
(
t
)
=
P
^
(
t
)
+
I
~
d.


## 2.3. The Proposed Wavelet-Based Method

Instead of employing a Fourier series, this study applied wavelet denoising to determine the smooth seasonal mean and perturbation term from rainfall and runoff data in the seasonal model. The procedures of the proposed wavelet-based method are shown below.(i)
First, rainfall and runoff time series data were decomposed into approximations and detailed signals at various resolution levels using the DAUB4 wavelet transform or (1) and (3).
(ii)
Threshold quantifications were performed at each resolution level using (4)–(6).
(iii)
The denoised rainfall and runoff time series data were obtained from the approximation at the final resolution level and the processed detailed signals by threshold quantifications at all rainfall and runoff resolution levels using wavelet reconstruction or (2) and (3).
(iv)
The denoised rainfall and runoff time series data were considered the smooth seasonal means employed in the seasonal model using (8)–(11).
(v)
The value of the original time series minus the denoised time series using (12) was considered the perturbation term P
(
t
) employed in the seasonal model.
(vi)
Based on calibrated rainfall and runoff data, the probability distribution of the percentage coefficients can be determined.
(vii)
After determining the probability distribution of the percentage coefficients, the estimated percentage coefficientsk
^ in the validated rainfall and runoff data can be obtained using (13) through random generation in a specific range of the probability distribution.
(viii)
The estimated percentage coefficientsk
^ multiplied by the smooth seasonal mean I
~
d using wavelet denoising equals the estimated perturbation term P
^
(
t
).
(ix)
The estimated rainfall or runoffI
^
(
t
) can be obtained from the sum of the estimated perturbation term P
^
(
t
) and the smooth seasonal mean I
~
d; that is, I
^
(
t
)
=
P
^
(
t
)
+
I
~
d.


## 2.4. Application and Analysis

### 2.4.1. Study Basin

This study confirmed the efficacy of the proposed wavelet-based method for modeling daily rainfall and runoff using the Wu-Tu watershed in northern Taiwan as an example. The upstream area of the watershed is 203 km2 (Figure 2), and the mean annual precipitation in this watershed is 2500 mm. Because of the topography of this watershed, the runoff pathlines are short and steep, and the rainfall is nonuniform in both time and space. Large floods that develop rapidly in the middle to downstream reaches of this watershed can cause severe damage. After collecting daily rainfall and runoff data from 1966 to 1994, data from 1966 to 1980 were used to calibrate the proposed model. Data from 1981 to 1994 were then used to assess the performance of the proposed method. Observed daily rainfall data were obtained using the kriging method to average the values obtained at the Jui-Fang, Huo-Shao-Liao, and Wu-Tu weather stations.Figure 2
A map of the Wu-Tu watershed study area located near Taipei, Taiwan.


### 2.4.2. Comparison of Model Performances

To quantitatively compare Fourier smoothing with wavelet denoising in a seasonal model, the validated results were evaluated based on the mean square error (MSE). The MSE is defined as(14)
MSE
=
∑
i
=
1
N
[
q
(
i
)
-
q
^
(
i
)
]
2
N
,
where q
^
(
i
) denotes the simulated discharge of the hydrograph for time period i
(
m
3
/
s
), q
(
i
) is the discharge of the observed hydrograph for time period i
(
m
3
/
s
), and N is number of data.


## 2.4.1. Study Basin

This study confirmed the efficacy of the proposed wavelet-based method for modeling daily rainfall and runoff using the Wu-Tu watershed in northern Taiwan as an example. The upstream area of the watershed is 203 km2 (Figure 2), and the mean annual precipitation in this watershed is 2500 mm. Because of the topography of this watershed, the runoff pathlines are short and steep, and the rainfall is nonuniform in both time and space. Large floods that develop rapidly in the middle to downstream reaches of this watershed can cause severe damage. After collecting daily rainfall and runoff data from 1966 to 1994, data from 1966 to 1980 were used to calibrate the proposed model. Data from 1981 to 1994 were then used to assess the performance of the proposed method. Observed daily rainfall data were obtained using the kriging method to average the values obtained at the Jui-Fang, Huo-Shao-Liao, and Wu-Tu weather stations.Figure 2
A map of the Wu-Tu watershed study area located near Taipei, Taiwan.


## 2.4.2. Comparison of Model Performances

To quantitatively compare Fourier smoothing with wavelet denoising in a seasonal model, the validated results were evaluated based on the mean square error (MSE). The MSE is defined as(14)
MSE
=
∑
i
=
1
N
[
q
(
i
)
-
q
^
(
i
)
]
2
N
,
where q
^
(
i
) denotes the simulated discharge of the hydrograph for time period i
(
m
3
/
s
), q
(
i
) is the discharge of the observed hydrograph for time period i
(
m
3
/
s
), and N is number of data.


## 3. Results and Discussion

In the conventional seasonal model, the smooth seasonal mean and corresponding perturbation terms are obtained using the Fourier series from rainfall and runoff data. The smooth seasonal means obtained from calibrated and validated rainfall and runoff data using the Fourier series are shown in Figures3 and 4, respectively. This study applied the DAUB4 wavelet to decompose the observed rainfall and runoff time series data into approximate signal data and detailed signal data. The line plots denoting the wavelet in Figures 3 and 4 show that the smooth seasonal mean was obtained using wavelet denoising for the three-resolution-level decomposition of rainfall and runoff data. Loss of the true signal, which typically occurs when smoothing using the Fourier series, does not occur during wavelet denoising, as shown in Figures 3 and 4.(a) The smooth seasonal means obtained using the Fourier series and wavelet denoising for a three-resolution-level decomposition from calibrated rainfall data. (b) The smooth seasonal means obtained using the Fourier series and wavelet denoising for a three-resolution-level decomposition from calibrated runoff data.
(a)
(b)(a) The smooth seasonal means obtained using the Fourier series and wavelet denoising for a three-resolution-level decomposition from validated rainfall data. (b) The smooth seasonal means obtained using the Fourier series and wavelet denoising for a three-resolution-level decomposition from validated runoff data.
(a)
(b)Based on the calibrated data from 1966 to 1980, the minimum and maximum percentage coefficients obtained from rainfall and runoff data using Fourier smoothing were shown in Table1(a). The range of k was divided into 6 intervals. The minimum and maximum percentage coefficients obtained from rainfall and runoff data using wavelet denoising for various resolution levels were shown in Table 1(b). Decompositions involving more than 4 resolution levels are not performed because of the short length of decomposed data and because of the substantial boundary effects [15]. For example, the lengths of decomposed data for 1, 2, 3, and 4 resolution levels are 184, 92, 46, and 23 days, respectively. To compare Fourier smoothing with wavelet denoising in the seasonal model, the range of k was also divided into 6 intervals. Table 2 shows the probability distribution of the percentage coefficients obtained from rainfall and runoff data using Fourier smoothing and wavelet denoising, respectively.Table 1
(a) The minimum and maximum percentage coefficients obtained from rainfall and runoff data using Fourier smoothing. (b) The minimum and maximum percentage coefficients obtained from rainfall and runoff datausing wavelet denoising.
(a)
Types
Minimum
Maximum
Rainfall
−0.982
3.000
Runoff
−0.885
4.853
(b)
Types
Resolution level M
Minimum
Maximum
Rainfall
1
−0.879
1.757
2
−0.960
1.588
3
−0.966
2.919
4
−0.982
2.085
Runoff
1
−0.640
1.313
2
−0.731
1.389
3
−0.744
1.708
4
−0.861
2.056Table 2
The probability distribution of the percentage coefficients obtained from calibrated daily rainfall and runoff data using Fourier smoothing and wavelet denoising.
Rainfall data using Fourier smoothing
M
Range ofk
-
1.0
<
k
≦
-
0.5
-
0.5
<
k
≦
-
0.3
-
0.3
<
k
≦
0.0
0.0
<
k
≦
0.5
0.5
<
k
≦
1.0
1.0
<
k
≦
3.0
0
Probability
0.139726
0.158904
0.252055
0.306849
0.104110
0.038356
Cumulative probability
0.139726
0.298630
0.550685
0.857534
0.961644
1
Runoff data using Fourier smoothing
M
Range ofk
-
1.0
<
k
≦
-
0.5
-
0.5
<
k
≦
-
0.3
-
0.3
<
k
≦
0.0
0.0
<
k
≦
0.5
0.5
<
k
≦
1.0
1.0
<
k
≦
5.0
0
Probability
0.098630
0.147945
0.323288
0.315069
0.068493
0.046575
Cumulative probability
0.098630
0.246575
0.569863
0.884932
0.953425
1
Rainfall data using wavelet denoising
M
Range ofk
-
1.0
<
k
≦
-
0.5
-
0.5
<
k
≦
-
0.3
-
0.3
<
k
≦
0.0
0.0
<
k
≦
0.5
0.5
<
k
≦
1.0
1.0
<
k
≦
3.0
1
Probability
0.035616
0.063014
0.380822
0.490411
0.024658
0.005479
Cumulative probability
0.035616
0.098630
0.479452
0.969863
0.994521
1
2
Probability
0.071233
0.090411
0.347945
0.405479
0.068494
0.016438
Cumulative probability
0.071233
0.161644
0.509589
0.915068
0.983562
1
3
Probability
0.104110
0.131506
0.298630
0.364384
0.087671
0.013699
Cumulative probability
0.104110
0.235616
0.534246
0.898630
0.986301
1
4
Probability
0.109589
0.134247
0.279452
0.356164
0.093151
0.027397
Cumulative probability
0.109589
0.243836
0.523288
0.879452
0.972603
1
Runoff data using wavelet denoising
M
Range ofk
-
1.0
<
k
≦
-
0.5
-
0.5
<
k
≦
-
0.3
-
0.3
<
k
≦
0.0
0.0
<
k
≦
0.5
0.5
<
k
≦
1.0
1.0
<
k
≦
2.1
1
Probability
0.002740
0.032877
0.468493
0.476712
0.016438
0.002740
Cumulative probability
0.002740
0.035617
0.504110
0.980822
0.997260
1
2
Probability
0.019178
0.046575
0.468494
0.435616
0.021918
0.008219
Cumulative probability
0.019178
0.065753
0.534247
0.969863
0.991781
1
3
Probability
0.027397
0.079452
0.449315
0.402740
0.030137
0.010959
Cumulative probability
0.027397
0.106849
0.556164
0.958904
0.989041
1
4
Probability
0.049315
0.112329
0.383561
0.378082
0.063014
0.013699
Cumulative probability
0.049315
0.161644
0.545205
0.923287
0.986301
1Then, to quantitatively compare Fourier smoothing with wavelet denoising in the seasonal model, the number of times of random generationr was calculated as 10. In addition, random generation p was conducted 3 times for each random generation r. Figure 5 displays a representative validation result for the estimated rainfall and runoff obtained using Fourier smoothing and wavelet denoising (M
=
3) in the seasonal model. Regarding the MSE criterion, the average value of the validated results obtained from rainfall and runoff data using wavelet denoising for a three-resolution-level decomposition (78.91 and 18465) outperforms the validated results obtained using Fourier smoothing (112.50 and 62925). The MSE of the validated results obtained using wavelet denoising also outperforms that obtained using Fourier smoothing for each random generation.(a) The validation results for the estimated rainfall obtained using Fourier smoothing and wavelet denoising (M
=
3) in the seasonal model (the random generation for r and p is number 1 and number 2, resp.). (b) The validation results for the estimated runoff obtained using Fourier smoothing and wavelet denoising (M
=
3) in the seasonal model (the random generation for r and p is number 1 and number 2, resp.).
(a)
(b)In this study, the number of resolution levels appropriate for wavelet decomposition was set to three [15]. To confirm the efficacy of using three resolution levels for decomposition, the validation results are compared among 1, 2, 3, and 4 resolution levels. Figures 6 and 7 show various smoothing effects using 1, 2, 3, and 4 resolution levels for wavelet decomposition for smooth seasonal means obtained from calibrated and validated rainfall and runoff data. The averages of MSE of the rainfall validation results obtained using 1, 2, 3, and 4 resolution levels in the seasonal model are 38.25, 64.38, 78.91, and 90.78, respectively. The averages of MSE of the runoff validation results obtained using 1, 2, 3, and 4 resolution levels in the seasonal model are 8688, 13798, 18465, and 22969, respectively. Based on the MSE criterion, the validation results obtained using 1, 2, 3, and 4 resolution levels are better in the proposed wavelet denoising method than in the conventional seasonal model. Generally, the smoothing of the approximation correlates with the number of resolution levels. However, computing errors arise during decomposition, and they increase with the number of resolution levels. The error means the difference between the estimation and the observations when the decomposed series are used to the modeling or prediction of hydrological time series. The number of resolution levels should be determined when applying a wavelet transform to the hydrological time series modeling or rainfall-runoff processes. The smoothing of the detailed signals and the residual improve as the number of resolution levels increases. However, greater resolution level takes more models for modeling and leads to the error propagation. The error increases with the number of resolution levels and affects the accuracy of hydrological time series modeling. Therefore, the number of resolution levels used in the wavelet decomposition may not be too many or too few. In this study, the best result can be obtained using only one resolution level. However, the smoothing effect is insufficient for denoising, as Figures 6 and 7 show. To assess smoothing performance and accuracy based on the MSE criterion, the number of resolution levels in the proposed wavelet denoising method was set to three. Thus, the level chosen is a data dependent tradeoff between the desired amount of smoothing and the performance of the proposed wavelet denoising method in comparison with the conventional seasonal model.(a) Various smoothing effects using 1, 2, 3, and 4 resolution levels for wavelet decomposition for smooth seasonal means obtained from calibrated rainfall data. (b) Various smoothing effects using 1, 2, 3, and 4 resolution levels for wavelet decomposition for smooth seasonal means obtained from calibrated runoff data.
(a)
(b)(a) Various smoothing effects using 1, 2, 3, and 4 resolution levels for wavelet decomposition for smooth seasonal means obtained from validated rainfall data. (b) Various smoothing effects using 1, 2, 3, and 4 resolution levels for wavelet decomposition for smooth seasonal means obtained from validated runoff data.
(a)
(b)The validation results show that smoothing the seasonal mean using wavelet transforms instead of the Fourier series in the seasonal model enhances the daily rainfall and runoff modeling accuracy. Compared to the conventional seasonal model, the proposed wavelet denoising method models daily rainfall and runoff processes more accurately because it calculates the smooth seasonal means and corresponding perturbation terms using wavelet transforms to separate the approximated and detailed signals.The uncertainty of any model mainly depends on input data, the parameter value, and the model structure. Assessment of data uncertainty and its effect on the evaluation should be emphasized more and more during each assessment and management [22, 23]. Otherwise, it may still be possible to get the “right answer for the wrong reasons” [24]. The wavelet denoising applied in the study actually reduces the uncertainty due to the input data.In general, observed rainfall and runoff data can be applied to the modeling or forecasting of rainfall or runoff time series. The main problem is how to select appropriate temporal scales for time series forecasting. Chou [25] provided a solution to determine the appropriate temporal scales in rainfall or runoff time series forecasting. Chou [25] applied sample entropy (SampEn) to rainfall and runoff time series to investigate the complexity of different temporal scales. Rainfall and runoff time series with intervals of 1, 10, 30, 90, and 365 days for the Wu-Tu upstream watershed were used. Thereafter, SampEn was computed for the five rainfall and runoff time series. Temporal scales with lower complexity (i.e., higher predictability) are used to provide a reference for choosing the appropriate temporal scales in the analysis and forecasting of rainfall and runoff time series. His results show that daily and annual data with low complexity and high predictability are recommended for use in rainfall and runoff time series forecasting. In this study, it is suitable to select daily rainfall and runoff data for analysis.Since a wavelet transform is a convolution operation, the duration of data may be insufficient in wavelet decomposition. If the duration of data is insufficient, wavelet decomposition will be adversely affected by the boundary effect. The 365-day duration of the study data in this work resulted in the boundary effect in decomposed and restructured processes of wavelet transforms. Thus, the duration of data was extended to 368 days via symmetrical extension [18].The seasonal means of rainfall and runoff processes, which are essential for assessing the hydrological properties of a river basin, should be smooth processes. Estimated seasonal means of rainfall and runoff processes typically produce oscillation because of the noise when data are not collected for a long duration. The novel wavelet denoising technique proposed in this study can obtain the smooth seasonal mean of rainfall and runoff processes and is suitable for modeling actual daily rainfall and runoff processes.Because the rainfall and runoff modeling in the seasonal model are obtained directly from rainfall and runoff data, and considering that the physical sense is not contained in parameters, there are some limitations in the applications of the seasonal model. However, the developments for the seasonal model are potential; its computation method is convenient and intuitive, its structure is simple, and its modeling scheme is suitable for the calculations of annual or monthly time series [21]. In addition, the test of the statistical distribution for rainfall or runoff time series can be avoided in the seasonal model because it can be applied to time series with any statistical distribution.


## 4. Conclusions

This study employed wavelet denoising instead of the Fourier series to obtain smoothed rainfall and runoff data. The data were considered the smooth seasonal mean in the seasonal model. Based on the MSE criterion, the validation results of the proposed wavelet denoising method are superior to those obtained using Fourier smoothing in the seasonal models. Multiresolution wavelet analysis effectively separates detailed data from approximate data. Noise is effectively removed from detailed data using a constant threshold, which enhances the modeling accuracy of daily rainfall and runoff processes.Instead of Fourier smoothing, this study applied wavelet denoising to acquire the smooth seasonal mean and corresponding perturbation term from daily rainfall and runoff data in traditional seasonal models, which use seasonal means for hydrological time series forecasting. With different decomposed resolution levels, the influence of different denoising effects on the results of the modeling of rainfall and runoff time series can be understood. The analytical results demonstrate that wavelet denoising not only obtains sufficient smoothing effect for denoising but also enhances the precision of daily rainfall and runoff modeling of the seasonal model, thereby being suitable for modeling actual daily rainfall and runoff processes.



---

*Source: 917365-2014-06-24.xml*