Study on Forecasting Break-Up Date of River Ice in Heilongjiang Province Based on LSTM and CEEMDAN

Abstract

In spring, rivers at middle and high latitudes in the Northern Hemisphere are prone to ice jams, which threaten the safety of hydraulic structures in rivers. Heilongjiang Province is located on the highest latitude in China, starting at 43°26′ N and reaching 53°33′ N. Rivers in Heilongjiang Province freeze in winter and break up in spring. Forecasting the break-up date of river ice accurately can provide an important reference for the command, dispatch, and decision-making of ice flood preventing and shipping. Based on the observed break-up date series of river ice from seven representative hydrological stations in Heilongjiang Province, the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) was used to decompose the observed break-up date series of river ice into several subsequences, and the long-short term memory neural network (LSTM) was used to forecast the subsequences decomposed by CEEDMAN. Then, the forecast results of each subsequence were summed to obtain the forecasting value for the break-up date of river ice proceeded by CEEMDAN-LSTM. Compared with the LSTM, the forecast accuracy of CEEMDAN-LSTM for the break-up date of river ice had been significantly improved, with the mean absolute error reduced from 0.80–6.40 to 0.75–3.40, the qualification rate increased from 60–100% to 80–100%, the root-mean-square difference reduced from 1.37–5.97 to 0.95–1.69, the correlation coefficient increased from 0.51–0.97 to 0.97–0.98, and the Taylor skill score increased from 0.87–0.99 to 0.99. CEEMDAN-LSTM performed well in forecasting the break-up date of river ice in the Heilongjiang Province, which can provide important information for command, dispatch, and decision-making of ice flood control.

1. Introduction

River ice plays a key role in the erosion of river banks and beds in cold regions and the regulation of river flow regimes [1,2]. In middle and high-latitude rivers in the Northern Hemisphere, the break-up of river ice is an important seasonal hydrological event [3]. The break-up of river ice has important ecological and socioeconomic impacts [4], and destroys fish spawning grounds [5], leads to an increase in the concentration of suspended sediment [6,7], pounds hydraulic structures, and threatens property safety [8,9,10]. Forecasting the break-up date of river ice is one of the important contents of the ice regime forecast, and provides an important scientific basis for command, dispatch, and decision-making in ice flood control.

Research regarding the forecasting of the break-up date of river ice began in the 1970s [11]. Numerical simulation methods and mathematical–statistical methods are used to forecast the break-up date of river ice [12]. One-dimensional numerical models (such as RIVJAM, ICEJAM, RIVER1D, ICEPRO, MIKE11, and RIVICE, etc.) [13,14,15,16,17,18,19], two-dimensional numerical models (such as CRISSP2D/DYNARICE, etc.) based on the principles of ice dynamics and hydrodynamics considering water level and flow data [20,21], and other numerical simulation methods were evaluated for simulating the breakup of river ice. However, due to the lack of some physical parameter values or boundary conditions required in the numerical simulation method, the application of numerical simulation methods for forecasting the break-up date of river ice in some rivers was limited [22]. For making up for the shortcoming of numerical simulation methods, mathematical–statistical methods were used to attempt to figure out the break-up date of river ice by the correlation between the break-up date of river ice and the meteorological or hydrological elements, or by the autocorrelation of the break-up date series of river ice. For example, the correlation between the thermal, dynamic, and river factors and break-up date was described by a fuzzy optimization neural network [23], the correlation between the flow velocity of upstream and downstream, and temperature and break-up date described by a neural network [24], the correlation between the autumn weather process and break-up date [12] was used to forecast the break-up date of rivers. However, the meteorological and hydrological elements affecting the break-up date of rivers are different at different hydrological stations, which need to be analyzed and recognized separately. Therefore, it is complex and laborious to use the correlation between the break-up date of river ice and the meteorological or hydrological elements in forecasting the break-up date of river ice for regions with numerous hydrological stations. In addition, it is difficult to measure or guarantee the measurement accuracy for some hydrological elements, such as water level, flow velocity, and so on, influenced by the river ice during the freezing period. Therefore, methods using the autocorrelation of the break-up date series, namely time-series forecast methods, have been found to be simple and effective ways to forecast the break-up date of river ice. Among the time-series forecast methods, artificial neural networks (ANN) with structural advantages perform well [25]. In ANNs, the long short-term memory network (LSTM), which can be used in time-series forecasting, overcomes the problem of traditional recurrent neural networks (RNNs) with exploding and vanishing gradient problems, by adding cell states [26]. Furthermore, the break-up date series of river ice contain nonlinear and nonstationary information that increases the difficulty of forecasting. Signal decomposition techniques, producing cleaner series as model inputs, can improve forecasting model performance in comparison to traditional approaches [27,28,29,30]. Compared with other signal decomposition techniques, complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) resolves the mode-mixing problem of empirical mode decomposition (EMD) and achieves negligible reconstruction errors [31], and is widely used in hydrology [32,33,34,35].

Therefore, to forecast the break-up date of river ice (the break-up date, for short, in the following) in the Heilongjiang Province, CEEMDAN is used to decompose the original break-up date series into subsequences (some intrinsic mode functions (IMFs) and one residual series) according to different fluctuations and frequencies, and LSTM is used to forecast the subsequences decomposed by CEEMDAN, then the forecast results of each subsequence are summed to obtain the forecast break-up date, and the method is named as CEEMDAN-LSTM. The remainder of the paper is organized as follows. Section 2 describes the details of the LSTM, CEEMDAN-LSTM, study area, and data. Section 3 presents the forecasting results of the break-up date obtained from LSTM and CEEMDAN-LSTM. Section 4 discusses the reasons for the difference in the accuracy of LSTM and CEEMDAN-LSTM, and the reasons for the poor forecast accuracy of the break-up date at the two special stations of Baoqing and Qiqihar stations. Conclusions are drawn in Section 5.

2. Materials and Methods

2.1. Study Area and Data

Heilongjiang Province is the northernmost and highest-latitude province in China. Heilongjiang Province, spanning 14 longitudes from east to west and 10 latitudes from north to south, geographically extends from 121°11′ E to 135°05′ E and 43°26′ N to 53°33′ N. Heilongjiang Province covers a total area of 473,000 km². There are three mountain ranges (Greater Khingan Mountains, Lesser Khingan Mountains, and the northern remnants of Changbai Mountain named Wanda Mountains) and two plains (Songnen plain and Sanjiang plain). The geographical location and topography of Heilongjiang Province are shown in Figure 1. Heilongjiang Province is located in cold and middle temperate zones with an average annual temperature of 3–5 °C. Heilongjiang province has a long and cold winter with an average temperature of −17.1 °C, and a windy and dry spring with an average temperature of 4.8 °C. The rivers in Heilongjiang province freeze up at the end of October and break up in April of the next year.

The break-up date series of Mohe hydrological station (MH), Yichun hydrological station (YC), Mudanjiang hydrological station (MDJ), Qiqihar hydrological station (QQHR), Harbin hydrological station (HRB), Fujin hydrological station (FJ), and Baoqing hydrological station (BQ) were selected to study the break-up date forecast. Among these hydrological stations, MH, YC, and MDJ are located in Greater Khingan Mountains, Lesser Khingan Mountains, and Wanda Mountains, respectively. QQHR and HRB are located in the Songnen plain. FJ and BQ are located in the Sanjiang plain. Except for MH, which is located in the cold temperate zone, the other hydrological stations are located in a middle temperate zone. The length, time span, mean, standard deviation (SD), and range of the break-up date series of each hydrological station are shown in

Table 1. Characteristic of break-up date series at each station.

Station	HRB	BQ	FJ	MH	MDJ	QQHR	YC
	(1)	(2)	(3)	(4)	(5)	(6)	(7)
Time span	1951–2019	1961–2019	1960–2019	1958–2019	1960–2019	1984–2019	1973–2019
Series length (year)	69	59	69	62	60	36	47
Mean	Day 98	Day 97	Day 107	Day 119	Day 101	Day 98	Day 107
Standard deviation (day)	5.14	6.84	5.21	5.07	5.44	6.82	5.20
Range (day)	25	45	23	25	25	34	19

Note: The break-up date in this paper is the day when the river surface is mainly covered by open water and the area of open water exceeds 80%. The time span of the break-up date series is the beginning year of the series to the end year of the series, and the unit of the time span of the break-up date series is the “year”. The length of the break-up date series is the number of samples that make up the series, and there is one sample a year for the break-up of rivers in Heilongjiang province. The unit of the length of the break-up date series is the “year”. In this paper, the break-up date of the river was converted to a numeric value, which is shown as the days after 1 January (1 January is the first day). For example, day 5 in this paper means the break-up date of the river with 5 January. After converting the break-up dates into corresponding numeric values, the break-up date series is converted into numeric value series, and the mean, standard deviation (SD), and range of the break-up date series are the mean value, the SD value, and the range value of the numerical value series, respectively. The units of mean, SD, and range of the break-up date series are the same, all of which are the “day”.

. As shown in Table 1, the break-up date series of all hydrological stations ended in 2019. Except for QQHR (length of 36) and YC (length of 47), the length of the break-up date series of other stations was about 60 years. Except for FJ, MH, and YC, breaking up in mid-late April (day 107–day 119), the break-up date of other hydrological stations is in early April (day 97 to day 101). Apart from BQ (standard deviation of 6.84, range of 45) and QQHR (standard deviation of 6.82, range of 34), the standard deviation of the break-up date series of other hydrological stations was less than 6, and the range (the difference between the earliest break-up date and the latest break-up date in the respective series) was less than 25 days.

2.2. CEEMDAN

In this paper, CEEMDAN, a signal decomposition technique, is used to decompose the observed break-up date series into subsequences according to different fluctuations and frequencies. The steps of CEEMDAN are as follows:

Step 1: Add adaptive Gaussian white noise n_j(t) to the break-up date series x(t), and obtain the signal x_j(t) used to decompose:

x_j(t) = x(t) + p_in_j(t)(i = 1, 2, ⋯, M; j = 1, 2, ⋯, N)

(1)

where i represents the number of the IMF in Step 2; j represents the realization number, that is, the number of times white noise was added to x(t) or r_k(t), the value of j in the hundreds will lead to a very good result [36], is set to 500 in this paper; p_i is the standard deviation in the noise, which controls the signal-to-noise ratio of n_j(t) to x(t), the value of p_i is set to 0.2 times the standard deviation of the break-up date series, which is recommended by Colominas, Schlotthauer, and Torres [37].

Step 2: Use the EMD [38] to decompose the signal obtained from Step 1 x_j(t) to obtain the IMF component $IM{F}_i^j \left(t\right)$. $IM{F}_i^j \left(t\right)$ is the ith IMF component of x_j(t) after EMD decomposition.

Step 3: Repeat Step 1 and Step 2 N times, and add different Gaussian white noise n_j(t) each time. The maximum iteration, which is not limited [39], is set to 5000 by experience to ensure that all extracted IMFs are valid before the EMD stops.

Step 4: Set the mean values of the IMF components after N time decompositions as the first IMF. The formula used to calculate the first IMF component is as follows:

$$\overline{IM{F}_1\left(t\right)}=\frac{1}{N}{\sum }_{j=1}^{N}IM{F}_1^j \left(t\right)$$

(2)

The residual signal is as follows:

$$r_1\left(t\right)=x\left(t\right)-\overline{IM{F}_1\left(t\right)}$$

(3)

Then, e₁(t) is defined as a process of the first IMF component $\overline{IM{F}_1\left(t\right)}$ by EMD and the sequence r₁(t) + p₂e₁(n_j(t)) to obtain the second IMF component as follows:

$$\overline{IM{F}_2\left(t\right)}=\frac{1}{N}{\sum }_{j=1}^N{e}_2\left(r_1\left(t\right)+p_2{e}_1\left(n_j\left(t\right)\right)\right)$$

(4)

The residual signal is as follows:

$$r_2\left(t\right)=r_1\left(t\right)-\overline{IM{F}_2\left(t\right)}$$

(5)

Similar to the above steps, the kth residual signal can be expressed as follows:

$$r_k\left(t\right)=r_{k-1}\left(t\right)-\overline{IM{F}_k\left(t\right)}$$

(6)

The (k + 1)th IMF component is as follows:

$$\overline{IM{F}_{k+1}\left(t\right)}=\frac{1}{N}{\sum }_{j=1}^N{e}_{k+1}\left(r_k\left(t\right)+p_{k+1}{e}_k\left(n_j\left(t\right)\right)\right)$$

(7)

Finally, the above steps are repeated until the residual signal $r_k\left(t\right)$ becomes a constant function or monotonic function. Supposing that there are M IMF components, the original sequence can be expressed as follows:

$$x\left(t\right)={\sum }_{i=1}^M\overline{IM{F}_i\left(t\right)}+r\left(t\right)$$

(8)

where r(t) represents the final residual signal that has become a constant function or monotonic function.

2.3. LSTM Network

The LSTM network is a special kind of RNN which reduces the decreased learning ability of traditional RNNs for long-term series [40]. There are three gates in the structure of LSTM designed to overcome the weakness in traditional RNNs to learn long-term dependencies.

The first gate is the forget gate, which controls the element of the cell state vector that will be forgotten. The formula used to calculate the output vector of the forget gate is shown in Formula (9)

$$f_t=\sigma \left(W_f{x}_t{}^{*}+U_f{h}_{t-1}+b_f\right)$$

(9)

where $x_t{}^{*}$ represents the normalized break-up date series (the normalization method is shown in Formula (10)); $f_t$ represents the output vector of the forget gate, whose values are in the range of (0, 1); $\sigma (\cdot )$ represents the logistic sigmoid function; h_t represents the hidden state; $W_f$, $U_f$ and $b_f$ represent the learnable parameters for the forget gate, $W_f$ and $U_f$ represent two adjustable weight matrices, and $b_f$ represents a bias vector.

$$x_t{}^{*}=\frac{x_t-\mu }{\sigma }$$

(10)

where $x_t$ is the sample value of the break-up date series, $\mu$ is the mean of the break-up date series, $\sigma$ is the standard deviation of the break-up date series.

The second gate is the input gate. This gate defines which information (and to what degree) is used to update the cell state at the current time, as shown in Formula (11):

$$i_t=\sigma \left(W_i{x}_t{}^{*}+U_i{h}_t{}_{-1}+b_i\right)$$

(11)

where $i_t$ represents a vector, and its values range from 0 to 1, $W_i and U_i$ represent two adjustable weight matrices, and $b_i$ represents a bias vector for the input gate.

$$\widehat{C_t}=\tan \mathrm{h}\left(W_c·\left[h_{t-1},x_t{}^{*}\right]+b_c\right)$$

(12)

$$C_t=f_t{C}_{t-1}+i_t\widehat{C_t}$$

(13)

where $C_t$ and $\widehat{C_t}$ are the cell state at the current moment and the previous moment, respectively. $\widehat{C_t}$ is the candidate status of the inputs. tanh is a hyperbolic tangent activation function.

The last gate is the output gate, which controls the information of the cell state $C_t$ that flows into a new hidden state. The output of the output gate is calculated by Formula (14):

$$O_t=\sigma \left(W_o{x}_t{}^{*}+U_o{h}_{t-1}+b_o\right)$$

(14)

$$h_t=O_t·\tan \mathrm{h}\left(C_t\right)$$

(15)

where $O_t$ represents a vector, and its values range from 0 to 1, $W_{o }and U_o$ represent two adjustable weight matrices, and $b_o$ represents a bias vector for the output gate.

In LSTM, the number of hidden layers is set to 18. The loss function is a mean square error function. More details of the LSTM can be found in Gers et al. (2000) [41]. In the original LSTM algorithm, a custom-designed approximate gradient calculation is used to update the weights after each timestep [42]. However, the backpropagation through time (BPTT) method is used to reverse-calculate the mean square error between the output value of each LSTM neuron and the observed value [43]. According to the mean square error, the gradient of each weight is calculated, and the gradient descent algorithm is applied to update the weight.

In addition to the model structure and weight updating, the window length is a part of LSTM. The window length is different for different model inputs and outputs [44]. For time-series forecast, the autocorrelation coefficient can be used to select the window length of LSTM [25]. The autocorrelation coefficient describes the degree of correlation between two periods of the random signal. The greater the absolute value of the autocorrelation coefficient, the stronger the correlation between the random signals over the two periods. The order (m) corresponding to the maximum value of the autocorrelation coefficient (maximum autocorrelation order) of the break-up date series represents the strongest correlation between the break-up date series x(t) and the break-up date series x(t − m) [45]. The value of the maximum autocorrelation order is set as the window length of the LSTM. The maximum autocorrelation order of the break-up date series at each hydrological station is shown in

Table 2. The maximum autocorrelation order of the break-up date series at each hydrological station.

Model	Station	HRB	BQ	FJ	MH	MDJ	QQHR	YC
LSTM	Maximum autocorrelation order	14	3	1	16	3	6	2

Note: The value of the maximum autocorrelation order is set as the window length of the LSTM.

. As shown in Table 2, the window length of the LSTM of HRB, BQ, FJ, MH, MDJ, QQHR, and YC are 14, 3, 1, 16, 3, 6, and 2, respectively. The last five years’ (2015–2019) break-up dates from each hydrological station are used as the forecast data of the LSTM, and the remaining break-up dates are used as training data. The forecast break-up date obtained by the LSTM is the calculation result of the LSTM after inverse normalization.

2.4. Fundamental Frameworks of CEEMDAN-LSTM

The coupling of CEEMDAN and LSTM, named CEEMDAN-LSTM, is used to forecast the break-up date of each station. The basic idea of CEEMDAN-LSTM is the normalized IMFs and residual series are obtained by CEEMDAN as input data for the LSTM. Same as the LSTM, the last five years’ break-up dates from each hydrological station are used as the forecast data of the CEEMDAN-LSTM, and the remaining break-up dates are used as training data. The value of the maximum autocorrelation order of each IMF and residual series is the window length of the LSTM. The forecast break-up date obtained by the CEEMDAN-LSTM is the sum of the calculation results of the LSTM after inverse normalization of all IMFs and residual. The specific steps of the CEEMDAN-LSTM are as follows:

(1)

Obtain the break-up date series and decompose it into subsequences (multiple IMFs and a residual series) by CEEMDAN.

(2)

Divide the IMFs and residual series into training data and forecast data, and normalize them.

(3)

Calculate the maximum autocorrelation order of each IMF and residual series (shown in

Table 3. The maximum autocorrelation order of each IMF and residual series of the break-up date series of each hydrological station.

Model	Station	HRB	BQ	FJ	MH	MDJ	QQHR	YC
CEEMDAN-LSTM	IMF1	1	2	2	1	2	1	2
	IMF2	3	3	3	3	3	3	3
	IMF3	1	1	1	4	1	1	1
	IMF4	1	1	1	1	1	1	1
	IMF5	1	1	1	1	1	-	1
	IMF6	-	1	-	-	-	-	-
	Residual	1	1	1	1	1	1	1

Note: “-” means no corresponding IMF.

), and build the LSTM to forecast the value of IMFs and residual series.

(4)

Denormalize and evaluate the forecast result.

2.5. Performance Evaluation

The Taylor diagram was adopted to assess and compare the performance of LSTM and CEEMDAN-LSTM in the break-up date forecast. The Taylor diagram is a polar coordinate diagram composed of standard deviation (SD), the root-mean-square difference (RMSD), and the correlation coefficient (R) of the forecast break-up date series ($\widehat{x_t}$) and the observed break-up date series (x_t) [46]. The Taylor diagram is constructed by the relationship between the four statistical quantities above shown in Formula (16). The SD, RMSD, and R can be calculated by Formula (17)–(20).

RMSD²= SD_o²+ SD_f² − 2SD_oSD_fR

(16)

$$S{D}_o=\sqrt{\frac{1}{N}{\sum }_{t=1}^N{\left(x_t-\overline{x_t}\right)}^2}$$

(17)

$$S{D}_f=\sqrt{\frac{1}{N}{\sum }_{t=1}^N{\left(\widehat{x_t}-\overline{\widehat{x_t}}\right)}^2}$$

(18)

$$RMSD=\sqrt{\frac{1}{N}{\sum }_{t=1}^N{\left[\left(\widehat{x_t}-\stackrel{=}{\widehat{x_t}}\right)-\left(x_t-\overline{x_t}\right)\right]}^2}$$

(19)

$$R=\frac{\frac{1}{N}{\sum }_{t=1}^N\left(\widehat{x_t}-\stackrel{=}{\widehat{x_t}}\right)\left(x_t-\overline{x_t}\right)}{S{D}_{o}S{D}_f}$$

(20)

where $S{D}_o$ is the standard deviation of observed break-up date series ($x_t$), $S{D}_f$ is the standard deviation of forecast break-up date series ($\widehat{x_t}$), N is the series length.

For quantitatively comparing the model performance, the Taylor skill score was introduced to evaluate the pros and cons of model capabilities. The calculation formula of the Taylor skill score is shown in Formula (21).

$$S=\frac{4{\left(1+R\right)}^4}{{\left(\frac{S{D}_f}{S{D}_o}+\frac{S{D}_o}{S{D}_f}\right)}^2{\left(1+R_{max}\right)}^4}$$

(21)

where S is the Taylor skill score. The larger the S, the better the model performs. R_max is the maximum value of R between the forecast break-up date and the observed break-up date.

In addition to the Taylor diagram and the Taylor skill score, the mean average error (MAE) and qualification ratio (QR) are selected to evaluate the model performance; the calculation formula of the MAE and QR are shown in Formula (22) and Formula (23), respectively.

$$MAE=\frac{1}{N}{\sum }_{t=1}^N\left|{\widehat{x}}_t-x_t\right|$$

(22)

When MAE = 0, all forecast break-up dates are equal to the observed break-up dates. The closer the MAE is to 0, the better the model performance.

$$QR=\frac{n_{abs\left(AE\right)\le 7}}{N}$$

(23)

where $n_{abs\left(AE\right)\le 7}$ is the number of the absolute error less than or equal to 7 days (according to the standard for hydrological information and hydrological forecasting, for ice regime forecasting with a foreseen period greater than 15 days, the maximum allowable error is 7 days), N is the length of series, AE = $x_t-{\widehat{x}}_t$.

3. Result

3.1. Results of Decomposing the Observed Break-Up Date Series Using CEEMDAN

The observed break-up date series of each hydrological station (except for the last five years’ data) was taken as the input of CEEMDAN; the decomposition result of CEEMDAN is shown in Figure 2.

As shown in Figure 2, the observed break-up date series is decomposed into several IMFs and a residual series. Combining Figure 2 and Table 1, it can be seen that the observed break-up date series with a larger standard deviation, compared with other similar break-up date series in length, obtains relatively more decomposed subsequences. Except for the break-up date series of QQHR and YC stations, the lengths of the break-up date series of the other stations are similar (Table 1). The standard deviation of the break-up date series of BQ station is the largest (Table 1), and the break-up date series of BQ station are decomposed into 7 components (6 IMFs and 1 residual series shown in Figure 2), while the break-up date series of the other hydrological stations are decomposed into 6 components (Figure 2).

As shown in Figure 2, the subsequences (IMFs and residual series) have different but simple intrinsic oscillation modes. Each IMF has a different frequency and fluctuation degree. With the decomposition processing, the frequency and fluctuation degree decreased. IMF1 has the highest fluctuation degree and frequency. The residual series has the lowest fluctuation degree and frequency. The residual series is close to linear and varied slightly around the long-term average, whose values are much larger than the values of IMFs. The residual series represents the trend in the original break-up date series [47].

3.2. Result of LSTM Applied to Forecast the Break-Up Date

The last five years’ break-up dates from each hydrological station were used as the forecast data of the LSTM, and the remaining break-up dates were used as input data. The absolute errors (AE) for the seven stations are shown in Figure 3. The mean absolute error (MAE) and qualification ratio (QR) and the range of AE are shown in

Table 4. The absolute error range, MAE, and QR of LSTM at each station.

Period	Performance Evaluation Index	HRB	BQ	FJ	MH	MDJ	QQHR	YC
Training	Range of AE	[−8, 9]	[−9, 9]	[−6, 5]	[−13, 9]	[−5, 4]	[−11, 12]	[−8, 9]
	MAE	2.53	2.61	1.69	2.48	0.80	3.97	2.71
	QR(%)	95.31	96.30	100.00	91.23	100.00	75.19	95.24
Forecast	Range of AE	[−6, 0]	[−10, 8]	[0, 4]	[−2, 4]	[1, 2]	[−13, −1]	[0, 3]
	MAE	1.80	6.40	2.00	2.67	2.17	6.33	1.67
	QR(%)	100.00	60.00	100.00	100.00	100.00	80.00	100.00

. The Taylor diagram is shown in Figure 4. The standard deviation (SD), the root-mean-square difference (RMSD), and the correlation coefficient (R) used to draw the Taylor diagram are shown in

Table 5. The S, RMSD, R, and SD of LSTM at each station.

Station	HRB	BQ	FJ	MH	MDJ	QQHR	YC
S	0.93	0.95	0.99	0.95	0.99	0.87	0.99
RMSD	3.58	3.96	2.27	3.87	1.37	5.97	3.10
R	0.72	0.82	0.91	0.66	0.97	0.51	0.83
SDo	5.14	6.84	5.21	5.07	5.44	6.82	5.20
SDf	3.95	5.50	5.23	4.05	5.65	4.64	5.27

.

As shown in Figure 3, among the absolute errors of break-up date forecasting for the seven representative hydrological stations, only the absolute errors of FJ and MDJ stations fall within the range of the maximum allowable error. As shown in Table 4, the absolute errors of all stations range from −13 to 12. The absolute errors in the training period range from −13 to 12, and the absolute errors in the forecast period range from −13 to 8, which is smaller than the absolute error range in the training period. In the training period, the absolute error range of QQHR is the largest, from −11 to 12, and the absolute error range of MDJ is the smallest, from −5 to 4. In the forecast period, the absolute error range of BQ is the largest, from −10 to 8, which is slightly larger than the absolute error range of QQHR from −13 to −1, and the absolute error range of MDJ is the smallest, from 1 to 2.

As shown in Table 4, the MAE values of all stations obtained by LSTM range from 0.80 to 6.40, and the MAE values of QQHR and BQ are larger than the others. In the training period, the MAE value of QQHR with 3.97 is the largest, which is slightly larger than the MAE value of BQ with 2.61. The MAE value of MDJ with 0.80 is the smallest in all stations. In the forecast period, the MAE value of BQ with 6.40 is the largest, which is slightly larger than the MAE value of QQHR with 6.33. The MAE value of YC with 1.67 is the smallest in all stations.

According to the maximum allowable error stipulated by the standard for hydrological information and hydrological forecasting, for an ice regime forecasting with a foreseen period greater than 15 days, the maximum allowable error is 7 days. Namely, a forecasting with the absolute error fewer than 7 days is considered as a qualified forecasting. The qualification ratio (QR) is listed in Table 4. As shown in Table 4, the QR values for all stations are above 60%. Furthermore, in the training period, except for the QR value of QQHR with 75.19%, the QR values are above 90%. In the forecast period, except for the QR value of QQHR with 80% and the QR value of BQ with 60%, the QR values are 100%.

The Taylor diagram provides a statistical summary of the matching degree between the forecasting break-up date series and the observed break-up date series in terms of correlation coefficient (R), standard deviation, and root-mean-square difference (RMSD) in a polar plot. In a Taylor diagram, the cosine of the azimuth angle represents the R between the forecasted break-up date series and the observed break-up date series, the distance from the origin (radius) represents the standard deviation of the forecasted (observed) break-up date series (SD_f or SD_o), and the distance from the reference point (azimuth angle = 0 and radius = SD_o), being in full agreement with observed break-up date series, represents the RMSD. In a Taylor diagram, the smaller the azimuth angle, the closer the radius to SD_o, and the closer the distance to the reference point, the better the LSTM performance.

As shown in Figure 4, the azimuth angles of the points representing the forecast break-up date series of FJ, MDJ, and YC are smaller, the radiuses of points representing the forecast break-up date series of FJ, MDJ, and YC are closer to their corresponding SD_o, and the points representing the forecast break-up date series of FJ, MDJ, and YC are closer to the corresponding reference points. The Taylor skill score (S) values of FJ, MDJ, and YC are all 0.99 (Table 5). The azimuth angle of the points representing the forecast break-up date series of QQHR is the largest, the radius of points representing the forecast break-up date series of QQHR is the farthest from the corresponding SD_o, and the points representing the forecast break-up date series of QQHR are the farthest from the corresponding reference points. The S value of QQHR is 0.87 (Table 5).

As shown in Table 5, the RMSD values range from 1.37 to 5.97. The RMSD value of QQHR with 5.97 is the largest, and is larger than the RMSD value of BQ with 3.96. The RMSD value of MDJ with 1.37 is the smallest. The values of R range from 0.51 to 0.97. The value of R of QQHR with 0.51 is the smallest. The value of R of MDJ with 0.97 is the largest. The SD_f values of FJ, YC, and MDJ are 5.23, 5.27, and 5.65, respectively, which are close to the SD_o values of FJ, YC, and MDJ with 5.21, 5.20, and 5.44, respectively. The difference between the SD_f and the SD_o of FJ, YC, and MDJ are 0.02, 0.07, and 0.21, respectively. The SD_f of QQRH is 4.64, which is quite different from the SD_o of QQHR with 6.82. The difference between the SD_f and the SD_o of QQHR is 2.18.

3.3. Result of CEEMDAN-LSTM Applied to Forecast the Break-Up Date

CEEMDAN-LSTM was used to forecast the break-up date of each station. The last five years’ samples of each IMF and residual series were used as the forecast data for a station, and the remaining data were used for training. The absolute errors (AE) of the simulation results of seven stations are shown in Figure 5. The mean absolute error (MAE) and qualification ratio (QR) and the range of AE are shown in

Table 6. The absolute error range, MAE, and QR of CEEMDAN-LSTM at each station.

Period	Performance Evaluation Index	HRB	BQ	FJ	MH	MDJ	QQHR	YC
Training	Range of AE	[−3, 3]	[−4, 4]	[−3, 4]	[−2, 3]	[−4, 3]	[−3, 3]	[−3, 4]
	MAE	0.75	1.97	0.79	0.77	0.87	1.35	1.01
	QR(%)	100.00	100.00	100.00	100.00	100.00	100.00	100.00
Forecast	Range of AE	[−1, 1]	[−6, −1]	[−2, 1]	[−2, 2]	[−1, 1]	[−4, 2]	[−2, 2]
	MAE	0.80	3.40	0.83	1.17	0.83	1.83	1.00
	QR(%)	100.00	100.00	100.00	100.00	100.00	100.00	100.00

. The Taylor diagram is shown in Figure 6. The standard deviation (SD), the root-mean-square difference (RMSD), and the correlation coefficient (R) used to draw the Taylor diagram are shown in

Table 7. The S, RMSD, R, and SD of CEEMDAN-LSTM at each station.

Station	HRB	BQ	FJ	MH	MDJ	QQHR	YC
S	0.99	0.99	0.99	0.99	0.99	0.99	0.99
RMSD	0.95	1.55	1.05	1.02	1.10	1.69	1.27
R	0.98	0.97 *	0.98	0.98	0.98	0.97 **	0.97 ***
SDo	5.14	6.84	5.21	5.07	5.44	6.82	5.20
SDf	5.23	6.45	5.31	5.15	5.42	6.73	5.51

Note: * represents the minimum R value (0.970 of BQ), ** represents the second minimum R value (0.971 of QQHR), *** represents the third minimum R value (0.974 of YC).

.

As shown in Figure 5, the absolute errors of all hydrological stations fall within the range of maximum allowable error. It can be seen that, compared with Figure 3, the absolute errors of CEEMDAN-LSTM of each station shown in Figure 5 are obviously reduced. As shown in Table 6, the absolute errors of all stations range from −6 to 4, the absolute errors in the training period range from −4 to 4, and the absolute errors in the forecast period range from −6 to 2. The absolute error range of each station in the training period is the same as the absolute error range of each station in the forecast period. In the training period, the absolute error range of BQ is the largest, from −4 to 4, and the absolute error range of MH is the smallest, from −2 to 3. In the forecast period, the absolute error range of QQHR is the largest, from −4 to 2, which is slightly larger than the absolute error range of BQ, from −6 to −1. The absolute error range of MDJ, is the smallest from −1 to 1. Compared with the absolute error range of LSTM (Table 4), the absolute error range of CEEMDAN-LSTM is reduced by 15 days: 17 days in the training period, and 13 days in the forecast period. Among all seven stations, the absolute error range of QQHR in the training period is reduced the most by 17 days, the absolute error range of BQ in the forecast period is reduced the most by 13 days, and followed by QQHR with 6 days. Similar to the LSTM, the absolute error range of CEEMDAN-LSTM in the training period is larger than the absolute error range in the forecast period. The absolute error range of BQ in the training period is the largest, and the absolute error range of QQHR in the forecast period is the largest, followed by BQ.

As shown in Table 6, the MAE values of all stations obtained by CEEMDAN-LSTM range from 0.75 to 3.40, and the MAE values of QQHR and BQ are larger than the others. In the training period, the MAE value of BQ with 1.97 is the largest, which is slightly larger than the MAE value of QQHR with 1.35. In the forecast period, the MAE value of BQ with 3.40 is the largest, which is slightly larger than the MAE value of QQHR with 1.83. Compared with LSTM, the range of MAE value obtained by CEEMDAN is reduced by 2.95. The MAE values of QQHR and BQ in the training period decrease by 2.62 and 0.64, respectively. The MAE values of QQHR and BQ in the forecast period decrease by 4.5 and 3.0, respectively.

In the training period and the forecast period, the QR values of all stations obtained by CEEMDAN-LSTM were 100%. Compared with LSTM, the QR value of QQHR obtained by CEEMDAN-LSTM in the training period increased by 25%, and the QR values of QQHR and BQ obtained by CEEMDAN-LSTM in the forecast period increased by 20% and 40%, respectively.

As shown in Figure 6, compared with the results of LSTM shown in Figure 4, the azimuth angles of the points representing the forecast break-up date series are smaller, the radiuses of the points representing the forecast break-up date series are closer to their corresponding SD_o values, and the points representing the forecast break-up date series are closer to the corresponding reference points. As shown in Table 7, the Taylor skill score of each station obtained by CEEMDAN-LSTM was improved to 0.99. The Taylor skill score of QQHR increased the most, by 0.12. The RMSD values obtained by CEEMDAN-LSTM range from 0.95 to 1.69. The RMSD value of QQHR is the largest with 1.69, followed by BQ with the RMSD value of 1.55. The R values obtained by CEEMAN-LSTM range from 0.97 to 0.98. The R values of QQHR and BQ are smaller with 0.97. Compared with Table 5, the RMSD values obtained by CEEMDAN-LSTM decrease by 0.27–4.28. The RMSD value of QQHR decreases the most with a decrease of 4.28. The R values obtained by CEEMDAN-LSTM improved by 0.01–0.46. The R value of QQHR improved the most with an increase of 0.46. Based on the analysis for Figure 6 and Table 7 above, the LSTM-CEEMDAN performs better than LSTM.

4. Discussion

4.1. Reasons for the Improvement of the Forecasting Accuracy of Break-Up Date by CEEMDAN-LSTM

Compared with the forecasting accuracy of LSTM shown in Figure 3 and Table 4, the forecasting accuracy of CEEMDAN-LSTM is obviously improved (Figure 5 and Table 6), which indicates that a decomposition series by CEEMDAN plays an important role in the improvement of the forecasting accuracy of break-up dates.

As shown in Figure 2, during the training period, with the decomposition of the series, the fluctuation degree of the subsequence gradually decreases. Generally speaking, the weaker the fluctuation degree, the weaker the randomness of the series, the easier it is to forecast, and the higher the forecast accuracy. As shown in Figure 7 and

Table 8. The MAE of IMFs and residual series at each station.

Station	IMF1	IMF2	IMF3	IMF4	IMF5	IMF6	Residual
HRB	0.293	0.155	0.126	0.095	0.042	-	0.037
BQ	0.741	0.511	0.288	0.265	0.073	0.035	0.053
FJ	0.342	0.213	0.105	0.089	0.004	-	0.043
MH	0.362	0.233	0.134	0.100	0.071	-	0.054
MDJ	0.368	0.157	0.185	0.082	0.029	-	0.043
QQHR	0.569	0.298	0.158	0.170		-	0.163
YC	0.455	0.222	0.163	0.083	0.086	-	0.005

, IMF1 and IMF2 have a high fluctuation degree and strong randomness. The MAE values of IMF1 and IMF2 of each station are large, with the order of magnitude of 10⁻¹. The IMF5, IMF6, and residual series have a low fluctuation degree and weak randomness. The MAE values of IMF5, IMF6, and the residual series are small with the order of magnitude about 10⁻² to 10⁻³. The MAE value with an order of magnitude about 10⁻² to 10⁻³ has negligible influence on the forecasting accuracy for break-up date. The Taylor diagram, which comprehensively shows the forecast accuracy of each subsequence of each station, is shown in Figure A1 in the Appendix A.

In addition, the reason for the improvement of the forecast accuracy of the break-up date by CEEMDAN-LSTM may be related to the sample values of the IMFs and residual series. Generally speaking, in the case of two series with a similar fluctuation degree of the same length, the larger the sample values, the greater the absolute error. Comparing the fluctuation degree and length shown in Figure 2 and the MAE values shown in Table 4 and Table 8, it is found that the fluctuation degree and length of the IMF1, IMF2, and the observed break-up date series are similar, but the sample values of the IMF1 and IMF2 are about 1/10 of the observed break-up date series, and the MAE values of IMF1 and IMF2 are less than the MAE value of the observed break-up date series. In addition, although the sample values of residual series are similar to the observed break-up date series (Figure 2), the fluctuation degree of residual series is low, the randomness of residual series is weak, and the MAE value of the residual series is small.

Above all, IMFs and residual series obtained by CEEMDAN have smaller sample values, weaker fluctuation degree, and stronger regularity than observed series, which leads to a reduction of learning difficulty for LSTM and an improvement in forecasting accuracy. Therefore, the forecasting accuracy of CEEMDAN-LSTM is higher than that of LSTM.

4.2. Reasons for the Lower Forecasting Accuracy of BQ and QQHR by LSTM

As shown in Table 4 in Section 3.2, the forecasting accuracy of the break-up date at BQ and QQHR by LSTM is not as good as others in the forecast period. The standard deviation of BQ is the largest, followed by that of QQHR, and the length of the break-up date series of QQHR is the shortest, as shown in Table 1 in Section 2.1. It therefore seems that the standard deviation and the length of the observed break-up date series are the reasons for the lower forecast accuracy of BQ and QQHR.

The standard deviation, as is known, describes the fluctuation degree of the series. The larger the standard deviation is, the higher the fluctuation degree of a series, and the stronger the randomness of a series. Generally speaking, the higher the fluctuation degree and the stronger the randomness of the series, the more difficult it is to obtain a good forecast Compared to the standard deviation of each station shown in Table 1, the standard deviation of BQ and QQHR is larger than the standard deviation of other stations, which indicates that the observed break-up date series of BQ and QQHR had the larger fluctuation degree and stronger randomness. Therefore, the forecasting accuracy of BQ and QQHR is lower.

Not only the larger standard deviation but also the short series length is a reason for the lower forecasting accuracy of QQHR (Table 1). Compared with other stations, the series length of QQHR is shorter, and the training samples of QQHR are fewer. The size of the training sample plays an important role in the forecasting accuracy of artificial neural networks. The shorter the series length, the fewer features of the observed break-up date series learned by LSTM, and the weaker the forecasting ability of the LSTM. Therefore, the forecasting accuracy of QQHR is low. It is worth noting that, comparing the forecasting accuracy obtained by LSTM and CEEMDAN-LSTM, the forecasting accuracy of QQHR has a greater improvement, which indicates that CEEMDAN can reduce the influence of few samples on the forecasting accuracy of LSTM to some extent.

4.3. The Performance of LSTM in Forecasting the Black Swan Events

As shown in Figure 3 in Section 3.2 and Figure 5 in Section 3.3, the absolute errors of forecasting results for BQ and QQHR in 2019 are significantly larger than that of other years, because the black swan event, which is accidental, greatly influential, very difficult to forecast, and predictable after occurring, occurred in 2019 at BQ and QQHR. Namely, the break-up dates of BQ and QQHR in 2019 were much earlier than other years for all the observed break-up date series. Moreover, the LSTM, introducing a gate to determine whether the input data is important enough to be remembered, is a neural network that can remember long short-term information in history and forecast important events with very long intervals and delays in series. For the black swan events, the LSTM performs poorly because the extreme characteristics of the black swan events are not learned from the observed samples. The subsequences, decomposed from the observed break-up date series without the black swan event by CEEMDAN, cannot describe the extreme characteristics of the black swan event. Therefore, the LSTM and CEEMDAN-LSTM perform poorly in the forecast of the break-up date of BQ and QQHR in 2019. Above all, using the LSTM and CEEMDAN-LSTM to forecast the break-up date, the extreme samples in history should be included in the training data to provide the information of extreme events for LSTM and improve the forecasting ability for extreme events.

5. Conclusions

LSTM and CEEMDAN-LSTM methods were used to study the forecasting of the break-up date of seven representative hydrological stations in Heilongjiang Province based on topographic and climatic features. The following conclusions are obtained.

(1)

CEEMDAN decomposed the observed break-up date series into subsequences according to different fluctuations and frequencies. The observed break-up date series with a larger standard deviation, compared with a similar break-up date series in length, had relatively more decomposed subsequences. With the decomposition processing, the frequency and fluctuation degree of subsequence decreased, and the sample values of subsequence increased. The residual series had the lowest fluctuation degree and frequency, which was close to linear and varied slightly around the long-term average.

(2)

The subsequence decomposed by CEEMDAN with a lower fluctuation degree or smaller sample values compared with the observed series for LSTM obtained a higher forecasting accuracy. The IMF1 and IMF2 had smaller values, the MAE values for forecasting results of IMF1 and IMF2 were small with the order of magnitude of 10⁻¹. The IMF5, IMF6, and residual series had lower fluctuation degrees, and the MAE values of forecasting results of IMF5, IMF6, and residual series were small with the order of magnitude of 10⁻² and 10⁻³.

(3)

Among the performance evaluation of the LSTM for all seven stations, the absolute error ranged from −13 to 12, the MAE values ranged from 0.80 to 6.40, the QR values were above 60%, the RMSD values ranged from 1.37 to 5.97, the R values ranged from 0.51 to 0.97, and the S values ranged from 0.87 to 0.99.

(4)

The forecasting accuracy was obviously improved by LSTM coupled with CEEMDAN. CEEMDAN-LSTM performed better than LSTM. In the performance evaluation of the CEEMDAN-LSTM for all seven stations, the absolute error ranged from −6 to 4, the MAE values ranged from 0.75 to 3.40, the QR values improved to 100%, the RMSD values ranged from 0.95 to 1.69, the R values ranged from 0.97 to 0.98, and the S values improved to 0.99.

(5)

CEEMDAN can reduce the influence of the few samples on the forecasting accuracy of LSTM. The forecasting accuracy by LSTM was obviously improved after decomposing the observed break-up date series of QQHR with a short length by CEEMDAN. The MAE value of forecasting results for QQHR decreased from 6.33 to 1.83, the QR value was improved from 80% to 100%, and the S value was improved from 0.87 to 0.99.