Comparison of River Basin Water Level Forecasting Methods: Sequential Neural Networks and Multiple-Input Functional Neural Networks

To precisely forecast downstream water levels in catchment areas during typhoons, the deep learning artificial neural networks were employed to establish two water level forecasting models using sequential neural networks (SNNs) and multiple-input functional neural networks (MIFNNs). SNNs, which have a typical neural network structure, are network models constructed using sequential methods. To develop a network model capable of flexibly consolidating data, MIFNNs are employed for processing data from multiple sources or with multiple dimensions. Specifically, when images (e.g., radar reflectivity images) are used as input attributes, feature extraction is required to provide effective feature maps for model training. Therefore, convolutional layers and pooling layers were adopted to extract features. Long short-term memory (LSTM) layers adopted during model training enabled memory cell units to automatically determine the memory length, providing more useful information. The Hsintien River basin in northern Taiwan was selected as the research area and collected relevant data from 2011 to 2019. The input attributes comprised one-dimensional data (e.g., water levels at river stations, rain rates at rain gauges, and reservoir release) and two-dimensional data (i.e., radar reflectivity mosaics). Typhoons Saola, Soudelor, Dujuan, and Megi were selected, and the water levels 1 to 6 h after the typhoons struck were forecasted. The results indicated that compared with linear regressions (REG), SNN using dense layers (SNN-Dense), and SNN using LSTM layers (SNN-LSTM) models, superior forecasting results were achieved for the MIFNN model. Thus, the MIFNN model, as the optimal model for water level forecasting, was identified.


Introduction
Taiwan is located in the Northwest Pacific, where typhoons frequently strike, and most of its rivers are short and steep. During typhoons, short-term heavy rains and hyperconcentrated streamflows cause river surges in downstream areas, which increases the possibility of flood loss in areas along rivers. The Hsintien River basin ( Figure 1) in northern Taiwan is the main source of service water to the Taipei Metropolitan Area (with a population of approximately four million); hence, water source management and conservation within the basin are essential. Typhoon-related disasters in the Hsintien River basin, such as Typhoon Soudelor in 2015, prompted us to develop highly accurate water level forecasting models to inform flood prevention. In research on flood-stage estimation and prediction, physically-based models often yield accurate results, but developing such a forecasting system requires extensive computational time and various hydrogeomorphological data [1]. Examples of such models are MIKE Flood [2], MODFLOW [3], HEC-RAS [4,5], and CCCMMOC [6,7]. Another common method is data-driven forecasting models; Rainfall-runoff and river-stage forecasting models are established using machine learning (ML). The advantages of these methods are that they are less subject to the constraints of physical mechanisms and that they effectively identify and learn correlated patterns between input data sets and the corresponding target values [8,9]. Such ML-based methods can be used to construct various algorithm models, such as support vector regression [10], multilayer feedforward perceptron [11], self-organizing map [12], k-nearest neighbor [13], and random forest [14], on the basis of rainfall data at precipitation stations and water level or water flow data at hydrographic stations.
In general, artificial neural network (ANN)-based models can process all types of sequential data. Time series data, such as streamflow and river stage, are a special type of sequence data ordered by timestamps [15,16]. Among ANN-based models, sophisticated ML-based time series models, including recurrent neural networks (RNNs), have been developed to solve the problem of time series data extending the memory length on the time axis. Accordingly, on the basis of RNN, models capable of automatically determining the memory length (e.g., LSTM neural networks) have been established. Researchers have employed RNN-based models to solve flood-stage prediction problems [17][18][19].
In most studies on flood-stage estimation and prediction, ML-based models are created using ground observation data, namely one-dimensional (1-D) data such as rainfall, water flow, and water level. According to Gires et al. [20] and Ochoa-Rodriguez et al. [21], rain gauges provide relatively accurate point rainfall estimates near the ground surface; nonetheless, they cannot effectively capture the spatial variability of rainfall, which significantly affects hydrological systems and thus runoff modeling. With the advancement and popularization of remote-sensing equipment, remote-sensing data have become cheaper and more accessible. Therefore, radar-sensing data can add value to hydrological data. Remote-sensing data (e.g., ground radar reflectivity data, which are commonly applied) can be employed to estimate ground rainfall conditions. Radar reflectivity is the signal reflected by the electromagnetic wave actively emitted by radar through water vapor particles in the atmosphere [22][23][24]. Radar reflectivity images exhibit various colors that display the intensity of the echo signals received by radars, and these images provide information on water vapor in space. Researchers over the past few decades have adopted radar reflectivity data for rainfall estimation [25][26][27][28][29][30][31]. Radar reflectivity data are essential remote-sensing data that facilitate the prediction of In research on flood-stage estimation and prediction, physically-based models often yield accurate results, but developing such a forecasting system requires extensive computational time and various hydrogeomorphological data [1]. Examples of such models are MIKE Flood [2], MODFLOW [3], HEC-RAS [4,5], and CCCMMOC [6,7]. Another common method is data-driven forecasting models; Rainfall-runoff and river-stage forecasting models are established using machine learning (ML). The advantages of these methods are that they are less subject to the constraints of physical mechanisms and that they effectively identify and learn correlated patterns between input data sets and the corresponding target values [8,9]. Such ML-based methods can be used to construct various algorithm models, such as support vector regression [10], multilayer feedforward perceptron [11], self-organizing map [12], k-nearest neighbor [13], and random forest [14], on the basis of rainfall data at precipitation stations and water level or water flow data at hydrographic stations.
In general, artificial neural network (ANN)-based models can process all types of sequential data. Time series data, such as streamflow and river stage, are a special type of sequence data ordered by timestamps [15,16]. Among ANN-based models, sophisticated ML-based time series models, including recurrent neural networks (RNNs), have been developed to solve the problem of time series data extending the memory length on the time axis. Accordingly, on the basis of RNN, models capable of automatically determining the memory length (e.g., LSTM neural networks) have been established. Researchers have employed RNN-based models to solve flood-stage prediction problems [17][18][19].
In most studies on flood-stage estimation and prediction, ML-based models are created using ground observation data, namely one-dimensional (1-D) data such as rainfall, water flow, and water level. According to Gires et al. [20] and Ochoa-Rodriguez et al. [21], rain gauges provide relatively accurate point rainfall estimates near the ground surface; nonetheless, they cannot effectively capture the spatial variability of rainfall, which significantly affects hydrological systems and thus runoff modeling. With the advancement and popularization of remote-sensing equipment, remote-sensing data have become cheaper and more accessible. Therefore, radar-sensing data can add value to hydrological data. Remote-sensing data (e.g., ground radar reflectivity data, which are commonly applied) can be employed to estimate ground rainfall conditions. Radar reflectivity is the signal reflected by the electromagnetic wave actively emitted by radar through water vapor particles in the atmosphere [22][23][24]. Radar reflectivity images exhibit various colors that display the intensity of the echo signals received by radars, and these images provide information on water vapor in space. Researchers over the past few decades have adopted radar reflectivity data for rainfall estimation [25][26][27][28][29][30][31]. Radar reflectivity data are essential remote-sensing data that facilitate the prediction of regional ground rainfall as well as the rainfall-runoff process. Borga [32] described the promising potential of radar rainfall observations because their high spatial and temporal resolution and extensive areal coverage provide detailed information on precipitation events. Accordingly, the application of two-dimensional (2-D) radar images compensates for the insufficient 1-D spatial rainfall data collected from land-based observation stations.
ML-based neural network architecture has been utilized for image recognition, object detection, and computer vision. Networks with such architecture, including convolutional neural networks (CNNs) [33] and deep belief networks [34], are called deep learning (DL) models because they require additional hidden layers for the processing of images. These DL-based neural network models have emerged as powerful methods for learning feature representations automatically from image data [35,36]. The most commonly employed CNNs are composed mainly of multiple convolutional and pooling layers, which extract and strengthen image features to further identify objects. CNN algorithms have also been successfully applied to hydrological problems. For example, Wang et al. [37] and Kimura et al. [38] have adopted CNNs and input 1-D water level data for water level forecasting. Van et al. [39] developed a CNN with 1-D data for rainfall-runoff modeling. Kabir et al. [40] developed a CNN model for the prediction of fluvial flood inundation, and the CNN model was trained with outputs from a hydraulic model to predict water depths. In the aforementioned findings, 1-D time-series data retrieved from land-based observation stations were initially collated into image data sets during preprocessing and were subsequently analyzed using CNN models. In addition, Wang et al. [41] developed a CNN for flood susceptibility mapping by using various 2-D images (such as altitude, curvature, distance to rivers, and rainfall).
The above-mentioned studies have used either 2-D images or 1-D attribute data alone to construct DL-based models. To the best of my knowledge, in the literature on water level prediction for typhoons, using both 2-D (e.g., radar reflectivity images) and 1-D (e.g., attributes obtained from land-based observation stations) data simultaneously have seldom been used as input variables for DL-based neural networks.
The development of a DL-based model using data from multiple sources and in multiple dimensions is thus valuable. Accordingly, this study used DL-based neural network models for accurately forecasting the downstream water level in the catchment area during typhoons and selected the Hsintien River basin as the research area. In general, multiple-input functional neural networks (MIFNNs) developed in this study offer the following contributions: (1) The increasing accessibility of ground observation and remote-sensing data enabled us to flexibly consolidate data from multiple sources in two dimensions (1-D and 2-D) for training DL models. The ground rainfall, water level, and reservoir release for 1-D data and radar reflectivity images for 2-D data were adopted.
(2) Feature extractions are required during the construction of MIFNN models when image files (e.g., radar reflectivity images) serve as input attributes; this provides appropriate feature maps for model training. Hence, the convolutional and pooling layers for extracting image features were employed.
(3) Because of time delays in water levels, this study applied LSTM layers, which enabled memory block units to automatically determine the memory length during the training of ANN-based models to guarantee the most effective memory length for time series data on the time axis.

Region and Material
The Hsintien River is located in northern Taiwan (Figure 1), with a drainage area of approximately 916 km 2 and a total length of 84.6 km [42]. The Hsintien River catchment area contains two main tributaries, Nanshi Creek and Peishi Creek. The source of the Peishi Creek is at an altitude of approximately 600 m, and the Feitsui Reservoir is located within its basin. The Feitsui Reservoir, the second-largest reservoir in Taiwan with a total storage capacity of approximately 406 million m 3 , Remote Sens. 2020, 12, 4172 4 of 24 was established in 1981 to ensure the stability of the tap water supply to the Taipei Metropolitan Area [43]. The source of Nanshi Creek is at an altitude of approximately 1300 m. It passes along the Lan-Sheng Bridge and the Shang-Gui-Shan Bridge water level stations and then the Chu-Chih and Hsiu-Lung stations after merging with the Peishi Creek. Table 1 lists the elevation information regarding these water stations. Located in the densely populated Taipei Metropolitan Area, the Hsiu-Lung station monitors the river flood stage to prevent overflows and floods. Therefore, the water level measured at the Hsiu-Lung station as the research objective was selected.

Typhoons and Raw Data
The water level data after 2011 collected at the water level stations located in the Hsintien River basin are comprehensive; hence, this study retrieved the hydrological data related to historical typhoon events from 2011 to 2019. As demonstrated in Table 2, 16 typhoons caused major rainfall in the research area, which substantially increased the downstream water level in the catchment area. This study collected the water level records of typhoons from the Hsiu-Lung, Chu-Chih, Shang-Gui-Shan Bridge, and Lan-Sheng Bridge stations using 1 h as the time interval, and the statistical properties of the collected data are displayed in Table 3. However, the rainfall data in the catchment area of Nanshi Creek were separately collected because it is an unregulated river. The data from two rain gauges in the catchment area, namely the Da-Tong-Shan gauge (at an elevation of 500 m) and the Fu-Shan gauge (at an elevation of 916 m), were obtained. This study also examined reservoir-regulated flows discharged from the Feisui Reservoir, which affected Peishi Creek. The rainfall data and the reservoir release values are listed in Table 3.
The water level stations and rain gauges within the research area are currently within the jurisdiction of the 10th River Management Office of the Water Resources Agency, and the Feitsui Reservoir is governed by the Taipei Feitsui Reservoir Administration of Taipei City Government.

Radar Imagery
The 2952 radar mosaics with resolutions of 1024 × 1024 pixels were collected. The hourly radar reflectivity images were obtained from the Central Weather Bureau of Taiwan. As depicted in Table 2

Radar Imagery
The 2952 radar mosaics with resolutions of 1024 × 1024 pixels were collected. The hourly radar reflectivity images were obtained from the Central Weather Bureau of Taiwan. As depicted in Table  2

Typhoon Paths
The influence of typhoon paths on the rainfall and water level in the research area were examined. As depicted in Figure 3a, the path of the eye of Typhoon Saola, one of the aforementioned four typhoons, moves northward after touching the east coast of Taiwan and then moves along the northeast coast of Taiwan. In the radar reflectivity image in Figure 2a, the typhoon circulation covered the location of the research area. Because of the relatively long time for which it hovered over northern Taiwan, total rainfall values of 666 and 878 m were measured at the two precipitation stations (i.e., the Da-Tong-Shan and Fu-Shan stations), which were the highest rainfall values of the 16 typhoons. Moreover, the paths of typhoons Soudelor, Dujuan, and Megi were similar to those presented in Figure 3b-d. They continued to travel to the west after landing on the east coast of Taiwan. When typhoon circulations are lifted by mountain terrains, heavy rains are likely to occur in windward areas, resulting in exceptionally high rainfall in the catchment area.

Typhoon Paths
The influence of typhoon paths on the rainfall and water level in the research area were examined. As depicted in Figure 3a, the path of the eye of Typhoon Saola, one of the aforementioned four typhoons, moves northward after touching the east coast of Taiwan and then moves along the northeast coast of Taiwan. In the radar reflectivity image in Figure 2a, the typhoon circulation covered the location of the research area. Because of the relatively long time for which it hovered over northern Taiwan, total rainfall values of 666 and 878 m were measured at the two precipitation stations (i.e., the Da-Tong-Shan and Fu-Shan stations), which were the highest rainfall values of the 16 typhoons. Moreover, the paths of typhoons Soudelor, Dujuan, and Megi were similar to those presented in Figure 3b-d. They continued to travel to the west after landing on the east coast of Taiwan. When typhoon circulations are lifted by mountain terrains, heavy rains are likely to occur in windward areas, resulting in exceptionally high rainfall in the catchment area.

Model Development
The SNN and MIFNN forecasting models for water level estimation during typhoons were developed (the frameworks are displayed in Figure 4).

Model Development
The SNN and MIFNN forecasting models for water level estimation during typhoons were developed (the frameworks are displayed in Figure 4).

SNNs
SNNs, a type of network model composed of sequential multilayers, have a typical neural network structure. As depicted in Figure 4a, the input layer is available for the input of 1-D array data (i.e., the attribute data of water level, rainfall, and reservoir discharge). Multiple dense layers, or LSTM layers, can be applied as the hidden layers (or core layers) in the middle [44]. After the dropout layer, the simulation results are produced from the output layer. The dropout layer is to randomly drop the weights of a certain proportion of nodes between the hidden layers to prevent overfitting.

SNNs
SNNs, a type of network model composed of sequential multilayers, have a typical neural network structure. As depicted in Figure 4a, the input layer is available for the input of 1-D array data (i.e., the attribute data of water level, rainfall, and reservoir discharge). Multiple dense layers, or LSTM layers, can be applied as the hidden layers (or core layers) in the middle [44]. After the dropout layer, the simulation results are produced from the output layer. The dropout layer is to randomly drop the weights of a certain proportion of nodes between the hidden layers to prevent overfitting.
Core layers are generally applied as dense layers. A dense layer is a regular layer of neurons in a neural network. Each neuron receives input from all neurons in the previous layer, thus ensuring a dense connection. Dense layers are also called fully connected layers because the neuron nodes are connected to the neurons in the previous layer. Dense layer implements the operation: output = activation {dot product of (input, kernel) + bias}, where activation is the element-wise activation function passed as the activation argument, the kernel is a weight matrix created by the layer, and the bias is a bias vector created by the layer [45,46]. Dense layers are not appropriate for all types of problems, particularly time sequence problems.
RNNs are suitable for sequential information. In traditional neural networks such as multiple layer perceptrons, all inputs and outputs are assumed to be independent. RNNs are termed recurrent because they perform the same task for every element of a sequence, with the output depending on Core layers are generally applied as dense layers. A dense layer is a regular layer of neurons in a neural network. Each neuron receives input from all neurons in the previous layer, thus ensuring a dense connection. Dense layers are also called fully connected layers because the neuron nodes are connected to the neurons in the previous layer. Dense layer implements the operation: output = activation {dot product of (input, kernel) + bias}, where activation is the element-wise activation function passed as the activation argument, the kernel is a weight matrix created by the layer, and the bias is a bias vector created by the layer [45,46]. Dense layers are not appropriate for all types of problems, particularly time sequence problems.
RNNs are suitable for sequential information. In traditional neural networks such as multiple layer perceptrons, all inputs and outputs are assumed to be independent. RNNs are termed recurrent because they perform the same task for every element of a sequence, with the output depending on the previous computations [47,48]. However, RNNs are difficult to train and are detrimentally affected by vanishing or exploding gradients; therefore, they cannot solve the long-term dependency problem [49]. LSTM neural networks are a type of RNN-based model. An LSTM network has LSTM cell blocks instead of standard hidden layers [50]. The components of these cells are known as the input gate; the forget gate, and the output gate. The structural details of LSTM networks were elaborated by [51,52]. In LSTM networks, a gate mechanism is introduced to prevent backpropagated errors from vanishing or exploding, and this approach has subsequently been proved to be more effective than the use of conventional RNNs [53].

MIFNNs
As illustrated in Figure 4b, MIFNNs could be used for modeling with 1-D hydrographic data and 2-D radar imagery. Accordingly, radar reflectivity images were employed to enhance the rainfall data captured in the upstream catchment area because the application of only rainfall data retrieved at rain gauges in SNN models was insufficient to comprehensively describe the rainfall signals within the catchment area.
The input data had multiple sources and two dimensions; hence, this study developed a model capable of merging and integrating data into single inputs. As depicted in Figure 4b, the two input layers separately served as entrances for the 1-D and 2-D data. First, 2-D images were processed. The input layer received the original bitmaps and extracted features through image-processing techniques to provide more effective feature maps for model training. As revealed in Figure 4b, in the first stage (2-D feature extraction, which is marked by red dotted lines), a network structure comprising convolutional and pooling layers are generally adopted for feature extraction. For convolutional layers, input images are implemented through convolution by using various kernels [54]. Convolution is implemented in two steps: sliding and calculating dot products; that is, using filters that slide onto input images and continuing to calculate matrices to obtain dot products. The convoluted images are called feature maps [55]. However, features generated in convolutional layers are often excessively accurate; hence, the sensitivity of convolutional layer-generated features to edges decreases after passing through pooling layers, which reduces the number of subsequent calculations required. Detailed equation derivation can be found in the studies of [56,57].
The second stage (marked by blue dotted lines) in the figure is to infuse 1-D and 2-D data. First, flattened layers are used to flatten both 1-D and 2-D data; specifically, dimension transmissions are processed by flattened layers. Subsequently, concatenate layers (also called shared layers) are employed to convert data into 1-D arrays. Eventually, sequential modeling (marked by green-dotted lines) is conducted to verify model parameters. The structure in this part is identical to that of the SNN model.

Model Constructions
The presented models were implemented with the open-source scikit-learn and Keras libraries in Python 3.7 (Python Software Foundation, Wilmington, DE, USA) [58,59]. The MIFNN and SNN models were constructed with a functional application programming interface (API) and a Keras sequential API, respectively. Because the proposed DL-based models use two-dimensional matrices to calculate images, the advanced computer processing units of graphics cards are required. The computation environment comprised GeForce RTX 2080 Ti graphics card in the study.

Data Splitting
Before model construction, data sets were divided into training, validation, and testing data set. Four typhoons (i.e., typhoons Saola, Soudelor, Dujuan, and Megi) were selected as the testing data set that significantly influenced the water level at the Hsiu-Lung station. As presented in Table 2, the maximum flood stages of the remaining 12 typhoons not in the testing data set were relatively low, which may compromise the accuracy of the estimated maximum flood-stage values if used for model training.
Several new cross-validation procedures were introduced to overcome the defects of the old ones [60]. One of the methods, the leave-p-out method (where p denotes the cardinality of the testing data set), could resolve the problem I encountered during data splitting. The leave-p-out method involves the use of p observations as the testing data set and the use of the remaining observations as the training-validation data set. This is repeated in all ways to cut the original sample on a testing data set of p observations and a training-validation data set. To increase the data in the training-validation data set, I let p equal 1. That is, one of the four biggest typhoons as the test data and the remaining 15 typhoons were used as the training and validation data sets. The remaining 15 typhoons were then used to train the models through 10-fold cross-validation. Therefore, four models were built.

Selection of Lag Times
In this study, the water levels at the Hsiu-Lung station were influenced by the natural processes of rainfall-runoff and control flow from upstream reservoirs. Accordingly, for the proposed water level prediction models, the hourly water levels predicted were assumed to be a function of the respective related attributes and could be expressed as: where L H,t+i represents the forecast water levels at the Hsiu-Lung station for the i-h ahead forecast (i ranges from 1 to 6 h in this research); P F and P D denote the rain rates at Da-Tong-Shan and Fu-Shan rain gauges, respectively; R F is the release from Feitsui Reservoir; L H , L C , L S , and L L are the water levels at Hsiu-Lung, Chu-Chih, Shang-Gui-Shan Bridge, and Lan-Sheng Bridge stations, respectively; and t and ∆t a -∆t g are the indices of the time period and lag-time length, respectively. Note that in Taiwan, the most common release strategy for reservoir flood control operations involves predefined operation rules, in which releases are expressed as a function of reservoir variables (water levels) and hydrological inputs (reservoir inflows) [61,62]. Moreover, numerous studies presented in-depth reviews of reservoir operations and modeling on this topic [63][64][65][66][67]. Because this study focused on developing the downstream water levels prediction models, for simplicity, the records of reservoir outflow were employed as an attribute when the modeling process.
Acceptable lag times (∆t a -∆t g ) for the water level forecasts were first identified. Figure 5 illustrates the correlation coefficient (R) between water level with lead times (1-to 6-h ahead forecast) at Hsiu-Lung station and various attributes with lag times. The results indicated that a lag time of 1 h yielded the highest R values for the four water level stations and the Feitsui Reservoir because of the small catchment. For the two rain gauges, lag times of 4, 3, 2, 1, 0, and 0 h, respectively, 1 to 6 h ahead yielded the highest R values. Note that lag time = 0 refers to the current time.

Training SNNs
Section 3.1 details the use of dense and LSTM layers as hidden layers of SNN models. For the modeling in this section, this study adopted two types of hidden layers to establish models and named them the SNN-Dense and SNN-LSTM models. The calculation parameters of the SNN-Dense

Training SNNs
Section 3.1 details the use of dense and LSTM layers as hidden layers of SNN models. For the modeling in this section, this study adopted two types of hidden layers to establish models and named them the SNN-Dense and SNN-LSTM models. The calculation parameters of the SNN-Dense model were set as follows: the loss function of dense layers = mean squared errors; activation function = sigmoid; the dropout rate of the dropout layer = 0.2. The parameters of the SNN-LSTM model were set as follows: the loss function of LSTM layers = mean squared error, cell activation function = rectified linear units; and gate activation function = sigmoid. The settings of the output layer of the SNN-Dense and SNN-LSTM models were loss function = mean squared error and activation function = sigmoid.
The optimizer method used in SNN-based models was the adaptive moment estimation (Adam) optimization algorithm. Adam [68] optimizes the momentum and learning rate. The merit of the Adam optimizer is its ability to compute individual adaptive learning rates for various parameters from estimates of the first and second moments of the gradients [69].
Hyperparameters were calibrated for both SNN-based models, comprising the number of hidden layers and neuron nodes in a hidden layer. In Figure 6, the testing typhoon, Saola, was used as an example. The number of hidden layers was calibrated from 1 to 4. The optimal lengths of the hidden layers were mostly 3 and 2 for SNN-Dense and SNN-LSTM models, respectively (Figure 6a,b). Subsequently, the number of nodes in a hidden layer, adjusting them to values between 10 and 100, were calibrated. In Figure 6c,d, the optimal numbers of nodes were within the ranges of 20-30 and 30-60 for the SNN-Dense and SNN-LSTM models, respectively.

Training MIFNNs
In Section 3.2, 2-D radar reflectivity images were incorporated into the MIFNN model. As

Training MIFNNs
In Section 3.2, 2-D radar reflectivity images were incorporated into the MIFNN model. As

Training MIFNNs
In Section 3.2, 2-D radar reflectivity images were incorporated into the MIFNN model. As described in Section 2.1, because the Peishi Creek is regulated by the Feitsui Reservoir, I cropped the corresponding position of the image according to the longitudinal and latitudinal range of the Nanshi Creek catchment. The longitude of the Nanshi Creek catchment is 121.43° E-121.60° E, and its latitude is 24.68° N-24.89° N. Its corresponding area on radar images could be represented by approximately 28 × 28 pixels. Figure 7 presents the cropped radar reflectivity images of the Nanshi Creek catchment for 16 consecutive hours on 8 August 8, 2015, when Typhoon Soudelor made landfall in Taiwan on 8 August 2015, from 0000 LST to 1500 LST. The figure illustrates how the intensity of the radar reflectivity signals dBZ changes with the interaction between typhoon circulation and the land. The parameter settings of the MIFNN model are as follows: In the first stage (2-D feature extraction), the input layer receives a radar reflectivity of 28 × 28 pixels as input; this study then used a sequence of three convolutional and pooling layers as feature extractors. For the 1-to 6-h predictions, this study employed the difference input image on the basis of the closest correlation of lag times. According to Section 4, lag times of 4, 3, 2, 1, 0, and 0 h, respectively at 1 to 6 h ahead, yielded the greatest R values. For the convolutional layer setting, the kernel size = 3 × 3, and stride The parameter settings of the MIFNN model are as follows: In the first stage (2-D feature extraction), the input layer receives a radar reflectivity of 28 × 28 pixels as input; this study then used a sequence of three convolutional and pooling layers as feature extractors. For the 1-to 6-h predictions, this study employed the difference input image on the basis of the closest correlation of lag times. According to Section 4, lag times of 4, 3, 2, 1, 0, and 0 h, respectively at 1 to 6 h ahead, yielded the greatest R values. For the convolutional layer setting, the kernel size = 3 × 3, and stride (the number of grids to move when sliding the window) = 1. For the pooling layer, a max-pooling method was employed for the subsampling process.
In the second stage (infusing of 1-D and 2-D data), the 2-D feature maps from these feature extractors with 1-D data were flattened into vectors, concatenated into one long vector, and passed on to a fully connected layer for the following modeling. In the third stage (sequential modeling), LSTM layers were applied because during the modeling of SNN-based models; this study discovered that the errors for LSTM layers were lower than those for dense layers. The LSTM layer environmental setups are the same as in the previous section.
Hyperparameters were calibrated for MIFNN models, comprising the number of hidden layers and of neuron nodes in a hidden layer. In Figure 8, the testing typhoon, Saola, was used as an example. The optimal lengths of the hidden layers were 3, and the optimal number of nodes in a hidden layer was in the range of 120-160 nodes.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 25 (the number of grids to move when sliding the window) = 1. For the pooling layer, a max-pooling method was employed for the subsampling process.
In the second stage (infusing of 1-D and 2-D data), the 2-D feature maps from these feature extractors with 1-D data were flattened into vectors, concatenated into one long vector, and passed on to a fully connected layer for the following modeling. In the third stage (sequential modeling), LSTM layers were applied because during the modeling of SNN-based models; this study discovered that the errors for LSTM layers were lower than those for dense layers. The LSTM layer environmental setups are the same as in the previous section.
Hyperparameters were calibrated for MIFNN models, comprising the number of hidden layers and of neuron nodes in a hidden layer. In Figure 8, the testing typhoon, Saola, was used as an example. The optimal lengths of the hidden layers were 3, and the optimal number of nodes in a hidden layer was in the range of 120-160 nodes.

Simulations
A simulation analysis of each of the selected typhoons was implemented. Four testing typhoons were employed to identify the most suitable water level prediction models in this section. Subsequently, a conventional regression-based (REG) water level forecasting model was established to serve as a benchmark model.

Simulations
A simulation analysis of each of the selected typhoons was implemented. Four testing typhoons were employed to identify the most suitable water level prediction models in this section. Subsequently, a conventional regression-based (REG) water level forecasting model was established to serve as a benchmark model.

Testing Typhoons and Predictions
Water level changes at the Hsiu-Lung station mainly result from the upstream rainfall and the water released by the reservoir. Hence, a rainfall and reservoir release time-series image was developed. Figure 9 (typhoons Saola and Soudelor) and Figure 10 (typhoons Dujuan and Megi) depict the hyetograph at the Da-Tong-Shan gauge and the Fu-Shan gauge, the reservoir released at the Feitsui Reservoir, and the water levels at the Hsiu-Lung station during each typhoon. Figures 9 and 10 show that the peak water level of the Hsiu-Lung station is contributed primarily by upstream rainfall (because the peak water level at the Hsiu-Lung station occurred after maximum rainfall) and secondarily by reservoir release (because the maximum water discharge times of the four tested typhoons all occurred after the water level peaked at the Hsiu-Lung station). This study further confirmed that, for the four testing typhoons, the difference between the maximum precipitation time and the peak water level time at the Hsiu-Lung station (i.e., the lag time) was minimal (between 3 and 4 h). Furthermore, this study discovered that, for the four testing typhoons, the maximum precipitation values estimated at the Da-Tong-Shan gauge and the Fu-Shan gauge differed greatly. The results were as follows: the maximum rainfall values at the Da-Tong-Shan gauge and the Fu-Shan gauge were (50, Figures 11 and 12, when lead time = 1 h, no major deviations were visible between the predicted results and the observed values. However, when the forecasting horizons were longer, the differences between the models' overall forecast performance, actual peak water levels, and arrival times increased.

Evaluation
To evaluate the deviations between the predicted and observed values, the root-mean-squared error (RMSE) deviation results for lead times of 1-6 h were analyzed. The RMSE is calculated as follows: where n is the total number of observations, P t is the predicted water levels at time t, and O t is the observed water levels at time t.
For forecasting horizons of 1-6 h, the RMSE evaluation values demonstrated that the ordering of the sizes of the models' deviations ( Figure 13) was as follows, from largest to smallest: the REG model, the SNN-Dense and SNN-LSTM models, and then the MIFNN model. Figure 14 illustrates the scatter diagrams of the observed and predicted values for all tested typhoons. The figure reveals that, when the lead time was larger, the MIFNN model's predicted values were closer to its observed values compared with those for the REG, SNN-Dense, and SNN-LSTM models. In terms of correlation, the ordering of the values for the models' coefficients of determination (R 2 ) in the figure was as follows, from largest to smallest: the MIFNN model, the SNN-LSTM and SNN-Dense models, and then the REG model. Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 25 Figure 11 (typhoons Saola and Soudelor) and Figure 12 (typhoons Dujuan and Megi) depict the water level time series graphs of observations and the predictions of the proposed models. The thick gray solid lines (observation) represent the observed water levels of the Hsin-Lung station, and the remaining four lines are the forecasting results of the REG, SNN-Dense, SNN-LSTM, and MIFNN models. As seen in Figures 11 and 12, when lead time = 1 h, no major deviations were visible between the predicted results and the observed values. However, when the forecasting horizons were longer, the differences between the models' overall forecast performance, actual peak water levels, and arrival times increased. Figure 11. Simulation of water levels for typhoons Saola and Soudelor at lead times of (a) 1 h, (b) 3 h, and (c) 6 h.

Evaluation
To evaluate the deviations between the predicted and observed values, the root-mean-squared error (RMSE) deviation results for lead times of 1-6 h were analyzed. The RMSE is calculated as follows: Figure 11. Simulation of water levels for typhoons Saola and Soudelor at lead times of (a) 1 h, (b) 3 h, and (c) 6 h.
where n is the total number of observations, is the predicted water levels at time t, and is the observed water levels at time t. For forecasting horizons of 1-6 h, the RMSE evaluation values demonstrated that the ordering of the sizes of the models' deviations ( Figure 13) was as follows, from largest to smallest: the REG model, the SNN-Dense and SNN-LSTM models, and then the MIFNN model. Figure 14 illustrates the scatter diagrams of the observed and predicted values for all tested typhoons. The figure reveals that, when the lead time was larger, the MIFNN model's predicted values were closer to its observed values compared with those for the REG, SNN-Dense, and SNN-LSTM models. In terms of correlation, the ordering of the values for the models' coefficients of determination (R 2 ) in the figure was as follows, from largest to smallest: the MIFNN model, the SNN-LSTM and SNN-Dense models, and then the REG model.

Overall, Performance
To estimate the model deviation for all four testing typhoons, the mean absolute error (MAE), relative MAE (rMAE), RMSE, and relative RMSE (rRMSE) were adopted as indicators. Moreover, the efficiency coefficient (CE) indicating the goodness-of-fit measures and R were used. The relevant equations are as follows: was as follows, from largest to smallest: the MIFNN model, the SNN-LSTM and SNN-Dense models, and then the REG model.

Overall, Performance
To estimate the model deviation for all four testing typhoons, the mean absolute error (MAE), relative MAE (rMAE), RMSE, and relative RMSE (rRMSE) were adopted as indicators. Moreover, the efficiency coefficient (CE) indicating the goodness-of-fit measures and R were used. The relevant equations are as follows:

Overall, Performance
To estimate the model deviation for all four testing typhoons, the mean absolute error (MAE), relative MAE (rMAE), RMSE, and relative RMSE (rRMSE) were adopted as indicators. Moreover, the efficiency coefficient (CE) indicating the goodness-of-fit measures and R were used. The relevant equations are as follows: where P t is the average of all predicted water levels, and O t is the average of all observed water levels.
In Figure 15, model performance is compared within a lead time of 1-6 h. Moreover, to examine peak water levels and the deviations of their arrival times, the error percentage of peak water levels (EWp) and time error of peak water levels (TWp) were defined:  Table 4 presents the mean measured values of the deviation indicators within a lead time of 1-6 h. Table 4 indicates that among the four models, the MIFNN model most effectively reduced forecasting errors. The improvement rate of the REG model by using its CE results was defined: After calculation, the SNN-Dense, SNN-LSTM, and MIFNN models exhibited respective improvements of 7.49%, 10.99%, and 12.87% compared with the REG model. Moreover, to examine peak water levels and the deviations of their arrival times, the error percentage of peak water levels (EW p ) and time error of peak water levels (TW p ) were defined: where P k is the predicted peak water level for typhoon event k; O k is the observed peak water level for typhoon event k; m is the total number of typhoon events; T P k is the arrival time of the predicted peak water level for typhoon event k; and T O k is the arrival time of the observed peak water level for typhoon event k. Figure 16 depicts the results of EW p and TW p . (1) According to the EW p values of the models, the peak water level deviation rate of the REG model was within the range of 3.87%-13.19%, that of the SNN-Dense model was within 5.46%-11.21%, that of the SNN-LSTM model was within 4.74%-8.47%, and that of the MIFNN model was within 3.66-7.97%. (2) According to the models' TW p values, the arrival times of peak water levels of the REG model were within 1-6 h for a lead time of 1-6 h, those of the SNN-Dense were within 0.75-3.5 h, those of the SNN-LSTM model were within 0.5-2.75 h, and those of the MIFNN model were within 0.25-2.5 h. Table 4 further indicates that, on average, the MIFNN model more accurately forecasts peak water levels and arrival times than REG, SNN-Dense, and SNN-LSTM models.

Discussion
As described, the results obtained from the MIFNN model indicated considerable changes in the spatial distribution of rainfall in the catchment area examined. However, the limited ground observation data failed to comprehensively address the spatial changes in rainfall. This study thus employed radar reflectivity images to effectively describe the spatial changes in water vapor and to provide additional hydrophilic data for the catchment area. For REG, SNN-Dense, SNN-LSTM, and MIFNN models, the RMSE values for our water level forecasts were 0.653, 0.480, 0.383, and 0.317 m, respectively. The results indicate that the MIFNN had superior forecasting accuracy to the REG, SNN-Dense, and SNN-LSTM models without the incorporation of radar reflectivity images.

Discussion
As described, the results obtained from the MIFNN model indicated considerable changes in the spatial distribution of rainfall in the catchment area examined. However, the limited ground observation data failed to comprehensively address the spatial changes in rainfall. This study thus employed radar reflectivity images to effectively describe the spatial changes in water vapor and to provide additional hydrophilic data for the catchment area. For REG, SNN-Dense, SNN-LSTM, and MIFNN models, the RMSE values for our water level forecasts were 0.653, 0.480, 0.383, and 0.317 m, respectively. The results indicate that the MIFNN had superior forecasting accuracy to the REG, SNN-Dense, and SNN-LSTM models without the incorporation of radar reflectivity images.
In addition, the application of LSTM layers as the network structure of the MIFNN model enhanced the forecasting ability of the algorithm. The LSTM layer is a neural network containing memory cell blocks. Each memory cell block comprises memory cell units that retain state across time-steps as well as three types of specialized gate units that learn to protect, utilize, or destroy this state as appropriate [44,70]. Memory blocks have a gate for determining whether the input data are valuable enough to be remembered and output. Such blocks can memorize information from various time lengths to automatically determine the amount of information that is useful [71]. Therefore, in typhoon water level forecasting through time series data analysis, determining for how long preceding memories should be retained is crucial. By contrast, because the REG and SNN-Dense models lacked memory functions, limited data were available for use. The experimental results also verified that the two aforementioned models yielded less accurate forecasting performance compared with the SNN-LSTM and MIFNN models, which had memory functions.
Limitations and suggestions are as follows: in the present study, this study used the trial-anderror method to calibrate the model hyperparameters. According to Maier et al. [72], as the number of model parameters increases, it is the difficulty in finding the parameter value or combination of parameter values that results in the smallest model error. Therefore, future studies are likely to focus In addition, the application of LSTM layers as the network structure of the MIFNN model enhanced the forecasting ability of the algorithm. The LSTM layer is a neural network containing memory cell blocks. Each memory cell block comprises memory cell units that retain state across time-steps as well as three types of specialized gate units that learn to protect, utilize, or destroy this state as appropriate [44,70]. Memory blocks have a gate for determining whether the input data are valuable enough to be remembered and output. Such blocks can memorize information from various time lengths to automatically determine the amount of information that is useful [71]. Therefore, in typhoon water level forecasting through time series data analysis, determining for how long preceding memories should be retained is crucial. By contrast, because the REG and SNN-Dense models lacked memory functions, limited data were available for use. The experimental results also verified that the two aforementioned models yielded less accurate forecasting performance compared with the SNN-LSTM and MIFNN models, which had memory functions.
Limitations and suggestions are as follows: in the present study, this study used the trial-and-error method to calibrate the model hyperparameters. According to Maier et al. [72], as the number of model parameters increases, it is the difficulty in finding the parameter value or combination of parameter values that results in the smallest model error. Therefore, future studies are likely to focus on conducting a suitable optimization algorithm, such as the genetic algorithm, to increase the ability to find global optima, although this is generally at the expense of computational efficiency.
In the modeling process, this study used data that were collected from 2-D images of weather radars. Radar instruments can easily monitor and provide information on the spatial distribution of rainfall over larger areas [73]. However, the radar rainfall estimates face many acquisition limitations, such as the radar coverage is limited due to terrain complexity and the shading problem [74]. In addition, there are huge areas where only satellite observations available, so any global solution for precipitation prediction needs to incorporate satellite data [75]. Hence, the use of data from satellite measurements or data fusion of radar and satellite measurements is suggested for constructing DL-based models in the future.
Furthermore, in Taiwan, other existing forecast methods were using physical-based models such as hydrodynamics-hydrology models (e.g., SOBEK model [76]) and conceptual models (e.g., PWRI tank model [77]). Wei [78] demonstrated the multilayer perceptron neural networks achieved superior prediction accuracy compared with the conceptual PWRI rainfall-runoff model for flood forecasting in the Tsengwen Reservoir watershed in Southern Taiwan. The inputs of multilayer perceptron neural networks comprised the observed rainfalls, reservoir inflows, and ground weather data. This was similar to the model case of SNN-Dense in the paper. This study may suggest that applying the proposed models using deep learning combined with radar reflectivity images in other regions or compared the physical-based models with proposed models in the Hsintien River catchment in future research.

Conclusions
To precisely predict the downstream water level in catchment areas during typhoons, the water level-forecasting models were developed by using DL-based ANN models. In addition, due to the increasing accessibility of ground observation and remote-sensing data, methods are required for integrating data from multiple sources and in two dimensions (1-D and 2-D data) for DL-model training.
Two types of models: SNN and MIFNN models, were designed. The SNN, a typical neural network structure, was a network model constructed using sequential methods. The MIFNN model, however, could process multisource or multidimensional data, meaning that a network model was required to flexibly consolidate data. Specifically, when image files (e.g., radar reflectivity images) were employed as input attributes, feature extraction was necessary to provide effective feature maps for model training. Therefore, this study adopted convolutional and pooling layers to extract features and applied LSTM layers during model training to enable memory cell units to automatically determine the length of memory and to provide more useful preliminary information during training.
The Hsintien River in northern Taiwan was selected as the research area, and data concerning the changes in rainfall and water levels caused by 16 typhoons in the area from 2011 to 2019 were collected. Furthermore, the following input attributes were applied: 1-D data (i.e., water levels at river station, rain rate at rain gauges, and reservoir release of Feitsui Reservoir) and 2-D data (i.e., radar reflectivity mosaics obtained from the central weather bureau). I established models to forecast the water level at the Hsiu-Lung station and examined the results for Typhoon Saola in 2012, Typhoon Soudelor in 2015, Typhoon Dujuan in 2015, and Typhoon Megi in 2015.
This study forecasted the results for lead times of 1-6 h and reached the following conclusions. First, in terms of the predictions of water level for a typhoon, on average, the RMSE values of the REG, SNN-Dense, SNN-LSTM, and the MIFNN models were, respectively 0.653, 0.480, 0.383, and 0.317 m. This indicated that the MIFNN model had an apparently superior forecasting performance to the REG, SNN-Dense, and SNN-LSTM models. Compared with the REG model, the SNN-Dense, SNN-LSTM, and MIFNN models achieved improvements in the efficiency of 7.49%, 10.99%, and 12.87%, respectively. In addition, the peak water levels and their arrival times estimated by the REG, SNN-Dense, SNN-LSTM, and MIFNN models were 7.832 m and 3.333 h, 10.979 m and 1.792 h, 7.686 m and 1.417 h, and 6.150 m and 1.167 h, respectively. In brief, this study successfully developed a MIFNN model able to simultaneously process data from multiple sources and in multiple dimensions, provided a new technique for establishing water level forecasting models, and further enhanced these models' accuracy.