Neural Network-Based Hybrid Forecasting Models for Time-Varying Passenger Flow of Intercity High-Speed Railways

Abstract

Time-varying passenger flow is the input data in the optimization design of intercity high-speed railway transportation products, and it plays an important role. Therefore, it is necessary to predict the origin-destination (O-D) passenger flow at different times of the day in combination with the stable time-varying characteristics. In this paper, three neural network-based hybrid forecasting models are designed and compared, named Variational Mode Decomposition-Multilayer Perceptron (VMD-MLP), Variational Mode Decomposition-Gated Recurrent Unit Neural Network (VMD-GRU), and Variational Mode Decomposition-Bidirectional Long Short-Term Memory Neural Network (VMD-Bi-LSTM). First, the time-varying characteristics of passenger travel demand under different time granularities are analyzed and extracted by the VMD method. Second, three neural network prediction models are constructed to predict the passenger flow sequence after VMD decomposition and reconstruction. Experimental analysis is performed on the Guangzhou Zhuhai intercity high-speed railway in China, and the passenger flow at different time periods of the day under different time granularities is predicted. The following results were found: (i) The number of hidden neurons and the number of iterations of the hybrid forecasting model have a great impact on the prediction accuracy. The error of the VMD-MLP model fluctuates less and it performs more smoothly than both the VMD-GRU model and the VMD-Bi-LSTM model. (ii) The VMD-MLP, VMD-GRU, and VMD-Bi-LSTM models can basically reduce the MAPE error to less than 10%. With the increase of time granularity, RMSE and MAE errors tend to gradually increase, while the MAPE error tends to gradually decrease. (iii) For passenger flow under a smaller time granularity, the prediction accuracy of the VMD-MLP model is higher, while for passenger flow under a larger time granularity, the prediction accuracy of the VMD-GRU and VMD-Bi-LSTM models is higher. (iv) The proposed neural network-based hybrid models outperform the existing models and the hybrid models perform better than the single models.

1. Introduction

The rapid development of China’s intercity high-speed railway not only alleviates the pressure of urban traffic, but also realizes the high-speed connection between adjacent cities. The high-speed and high-density transportation organization of the intercity high-speed railway gives passengers more choices of departure time. The passenger travel demand for the intercity high-speed railway at different departure periods in one day shows obvious fluctuation characteristics, which have a certain stability in a short period of time. Such a demand can be called the time-varying passenger flow. Therefore, this stable fluctuation law of passenger flow can be used to predict the time-varying passenger flow of the intercity high-speed railway in the short term and provide decision support for the optimization design of intercity high-speed railway transportation products.

In recent years, the time-varying passenger flow of the intercity high-speed railway has been used as input data in the optimization of timetables, line planning, and differential pricing by various researchers [1,2,3,4,5]. The prediction accuracy of passenger demand with different departure times in one day directly impacts the optimization solutions for high-speed railway operation management. Therefore, the prediction of the time-varying passenger flow of the intercity high-speed railway becomes an important role.

Traditional railway passenger flow forecasting is mainly aimed at the daily, monthly, and yearly passenger flow, which lacks research on the time-varying passenger flow of the intercity high-speed railway [6,7,8]. Although, some short-term methods were developed to predict passenger flow at the time granularities of 5 min, 15 min, and half an hour [9,10], but they were mainly focused on the prediction of subway passenger flow. There are quite a few differences between the subway and high-speed railways.

Considering the above research gaps, the objective of this paper is to develop methods for the intercity high-speed railway to predict passenger flow under different time granularities in one day. Prediction methods have been developed for decades and are usually divided into parametric and non-parametric approaches. Traditional parametric methods [11,12,13] are simple and fast but less accurate for nonlinear real-time data. Non-parametric methods have been widely researched in recent years, among which, neural networks are hot topics as they are more suitable for passenger flow prediction with non-linear and more complex data [14,15]. Considering that single prediction methods always have some defects, some hybrid methods are proposed to further improve the prediction accuracy, especially the neural network-based hybrid forecasting method [15,16].

Based on the above analysis, the research ideas of this paper are as follows: passenger travel demand is analyzed based on historical passenger flow, with an operating time of one day divided into time periods with different time granularities of 1 h, 3 h, 5 h, and 16 h, and the origin-destination (O-D) passenger flow time series under different time granularities are obtained. Then, the non-stationary passenger flow time series is decomposed into a stationary passenger flow time series by using the Variational Mode Decomposition (VMD), to extract the time fluctuation characteristics of passenger flow. Finally, three neural network-based hybrid prediction models are proposed, named VMD-MLP, VMD-GRU, and VMD-Bi-LSTM, to predict the reconstructed passenger flow time series and obtain the O-D passenger travel demand at different times of the day in the future.

The novelty of the research is the proposal of three neural network-based hybrid prediction models for the intercity high-speed railway. The designed methods not only use the stable fluctuation law of passenger flow but also predict the O-D passenger flow under different time granularities in one day, which have not been sufficiently studied for the intercity high-speed railway in the current literature.

The main contributions of this paper are as follows: (1) the fluctuation characteristics of the O-D passenger flow of the intercity high-speed railway under different time granularity in one day are analyzed, and (2) based on the time fluctuation characteristics of passenger flow, three neural network-based hybrid prediction models: VMD-MLP, VMD-GRU, and VMD-Bi-LSTM, are designed to predict the O-D passenger flow of the intercity high-speed railway under different time granularities and different periods in one day.

This paper is arranged as follows: Section 2 summarizes the literature, Section 3 analyzes the characteristics of passenger travel demand and proposes a method of passenger flow time series decomposition; in Section 4, three neural network-based hybrid prediction models are designed, Section 5 presents the experimental analysis, and Section 6 is the conclusion.

2. Related Literature Review

The research methods of railway passenger flow forecasting can be divided into three categories: parametric, non-parametric, and hybrid forecasting models. Parametric forecasting models mainly include the Time Series Model, Historical Average Model, and the Autoregressive Comprehensive Moving Average Model [13,17]. The parametric models have been applied in passenger flow forecasting before, but they are sensitive to the linear relationship between variables and cannot grasp the nonlinear relationship of time series, for which they have certain limitations in the application range of passenger flow forecasting. Milenkovi, Libor, and Melicha [18] used the Autoregressive Comprehensive Moving Average method to predict the passenger flow of railway stations in the short term at the granularity of a month.

Non-parametric forecasting models mainly include the Support Vector Machine, Neural Network, Decision Tree, etc. Non-parametric forecasting models are flexible, can fit a large number of different function forms, and have more advantages in terms of their training ability and performance than parametric prediction models. Zhang, Cheng, and Gao [19] proposed a prediction method based on multi-layer LSTM for short-term prediction of railway station passenger flow with a time granularity of 0.5 h. This method combined multiple traffic data and feature selection based on Spearman correlation and had a high prediction accuracy. Jing and Yin [20] used a BP neural network to predict the short-term passenger flow of intercity high-speed railway stations at the granularity of one hour. They also adopted an improved step-updating method to prevent the error function from falling into a local minimum and oscillation during the weight-updating process, which showed that this model had good performance. Peng, Bai, and Wu [21] predicted the short-term passenger flow of railway stations with the time granularity of one hour, and proposed an improved LSTM model, which had higher accuracy than the traditional BPNN module and the LSTM model without improvement.

Toque et al. [22] proposed the Gated Recurrent Unit Neural Network Model (GRU) to predict the station passenger flow of various rail transit modes in the Paris Business District in the short term, with a time granularity of 15 min. Compared with the Random Forest Model and LSTM Model, it had greater accuracy. Jérémy, Gérald, and Stéphane [23] proposed the Bayesian model to predict the passenger flow of subway stations in the short term with a time granularity of 10 min. Marc, Peter, and Dieter [24] proposed the DLR four-step model to predict the O-D passenger flow of an airport in the long term, with a time granularity of one year. Ulrich and Bozana [25] proposed the GVAR model to predict the passenger flow of an airport in the short term, with a time granularity of one month. Chan, Yash, and Sameer [26] proposed the Neural Granger Causality model to also predict the passenger flow of an airport in the short term, with a time granularity of one month.

At present, most of the passenger flow forecasting studies use hybrid models to predict passenger flow. There are different combinations of hybrid models. One approach is to combine two or more parametric or non-parametric models, and the other is to combine the passenger flow sequence decomposition model with parametric or non-parametric models to give full play to the advantages of each model and achieve a higher prediction accuracy.

(1) Passenger flow time series decomposition model mixed with parameter forecasting models or non-parameter forecasting models

Jiang, Zhang, and Chen [27] used a hybrid model, EEMD-GSVM, which decomposed the time series into several intrinsic mode functions using EEMD, and then used the particle swarm-calibrated GSVM to predict the daily O-D passenger flow of an intercity high-speed railway in the short term. Zhao and Mi [28] proposed a hybrid model, SSA-WPDCNN-SVR, to predict the daily O-D passenger flow of an intercity high-speed railway in the short term. First, the SSA method was used to decompose the passenger flow time series into a main series and several sub-series, then the Convolutional Neural Network Model was used to predict the main series, while the SVR model was used to predict the subseries.

Sun et al. [29] used the hybrid models VMD-SARIMA-MLP and VMD-SARIMA-LSTM. The SARIMA model was used to predict the periodic components, the LSTM was used to learn and predict the deterministic components, and the MLP network was used to predict the volatile components to predict the short-term daily passenger flow of the subway station. Xin et al. [30] proposed a wave-LSTM hybrid model to predict the short-term passenger flow of subway stations with 15 min as the time granularity. Wei and Chen [31] proposed a hybrid model, EMD-BPN, to predict the short-term passenger flow of subway stations, with 15 min as the time granularity. Jin et al. [32] used a hybrid model, VMD-ARMA/KELM-KELM, including Variational Mode Decomposition (VMD), Autoregressive Moving Average Model (ARMA), and Kernel Limit Learning Machine (KELM), which, respectively, predicted the short-term passenger flow of airlines with a monthly time granularity.

It is noted that Empirical Mode Decomposition (EMD) and Variational Mode Decomposition (VMD) are usually used as decomposition techniques to improve the prediction accuracy, which were proposed by Huang et al. in 1998 [33] and by Dragomiretskiy and Zosso in 2014 [34], respectively. Some researchers [35,36,37] also performed comparison studies and found that the performances of VMD were better than EMD, which was partly because EMD had some defects, such as mode mixing and an endpoint effect, while VMD could ensure decomposition optimality in the sense of minimum sum of modes’ bandwidths.

(2) Mixed parameter forecasting models or non-parameter forecasting models

Wen et al. [38] used a hybrid forecasting model, SARIMA-TrAdaboost, which decomposed the time series into linear and non-linear time series and predicted the short-term daily passenger flow of Beijing–Shanghai high-speed railway stations by embedding a random forest model. Glisovic, Milenkovi, and Bojovi [39] proposed two hybrid forecasting models, SARIMA-GA-ANN and SARIMA-ANN, to predict the short-term passenger flow of railway stations with a monthly time granularity. Zhu and Zhou [8] used the panel vector autoregression (PVAR) and neural network (NN) hybrid PVAR-NN prediction methods to predict yearly passenger flow in the railway system.

Jérémy, Stéphane, and Gérald [40] used a hybrid forecasting model, Bayesian and Gaussian methods, to predict the short-term passenger flow of subway stations, with 2 min as the time granularity. Ma, Guo, and Ma [41] proposed GCN-Bi-LSTM to predict the short-term passenger flow of subway stations with 5 min as the time granularity. Zhang, Chen, and Shen [42] proposed a hybrid forecasting model, CB-LSTM, to predict the short-term passenger flow of subway stations, with 1 to 30 min as the time granularity. Chen et al. [15] developed an EMD-LSTM hybrid forecasting model for predicting short-term metro-inbound passenger flow, with 15 min as the time granularity.

Robert, Yingqi, and Suzilah [43] used the hybrid forecasting model ADL-TVP-VAR to predict the passenger flow of an airport in the medium and long term with a time granularity of one year. Kim and Shin [44] used the hybrid forecasting model NAF-k-fold to predict the passenger flow of an airport in the short term with a time granularity of one month. Rodrigo [45] used the hybrid model MLEM to train multiple neural networks to predict the passenger flow of an airport in the long term with a time granularity of one month. Wai et al. [46] used the hybrid model ARIMA-SRIMA to predict the passenger flow of an airport in the short and long term with a time granularity of one month.

In conclusion, the prediction accuracy of the hybrid model is higher than that of the single model, and hence hybrid prediction methods are widely used in recent years. Railway passenger flow forecasting mainly forecasts short-term daily or monthly passenger flow, and only a few scholars have studied hourly passenger flow forecasting of railway stations. However, there is a lack of research on O-D passenger flow forecasting of different departure periods under different time granularities in one day. Therefore, this paper uses the time distribution characteristics of intercity high-speed railway passenger flow to design neural network-based hybrid prediction models, that is, the Variational Mode Decomposition-Multilayer Perceptron (VMD-MLP) model, the Variational Mode Decomposition-Gated Recurrent Unit Neural Network (VMD-GRU) model, and the Variational Mode Decomposition-Bidirectional Long Short-Term Memory Neural Network (VMD-Bi-LSTM) model, to predict the O-D passenger flow of an intercity high-speed railway at different times of one day under different time granularities.

3. Data Analysis and Decomposition

3.1. Time-Varying Passenger Flow

The passenger flow along the intercity high-speed railway is huge, the train departure frequency is high, and the passenger travel demand shows obvious time-varying characteristics, which means that the passenger flow at different departure periods in one day has obvious fluctuations and reflects certain periodic characteristics. According to the historical ticket sale data of the Guangzhou–Zhuhai intercity high-speed railway in China, the time-varying characteristics of passenger flow were analyzed. The passenger flow with time-varying characteristics can be called time-varying passenger flow.

Statistical analysis of the passenger flow in the dimensions of the departure period and departure week can obtain the distributions of O-D passenger flow under different time granularities every day in a week. The passenger demand before 7:00 and after 23:00 is small, so only the passenger flow between 7:00 and 23:00 was studied. For convenience, the time granularity is denoted as $\Delta t$, and $\Delta t=1,3,5,16$. This paper adopted time granularity $\Delta t$, dividing the operation time [7:00, 23:00] into several time periods, [7, 7 + $\Delta t$] as in period 1, [7 + $\Delta t$, 7 + $2\Delta t$] as in period 2, and so on. The passenger flow distributions from Guangzhou South station to Zhongshan North station under different time granularities are shown in Figure 1.

As can be seen from Figure 1:

• The fluctuation characteristics of passenger flow every day in a week have a certain periodicity.
• The fluctuation characteristics of the passenger flow at different departure times every day have a certain tendency.
• The fluctuation characteristics of passenger flow under different time granularities are obviously different. The smaller the time granularity is, the more detailed and more complex the time-varying characteristics of the passenger flow will be.

The actual passenger flow in the $m$th period of day $n$ is recorded as $x_n^m$, and the predicted passenger flow is recorded as ${\widehat{x}}_n^m, n=1,2\dots \dots N, m=1,2,3,\dots$. According to the time-varying characteristics of passenger flow in Figure 1, for an O-D pair, the passenger flows in the same period of the first 14 days are used to predict the passenger flow in the same period of the 15th day, which means that $x_{n-14}^m, x_{n-13}^m, \dots , x_{n-1}^m$ are used to predict ${\widehat{x}}_n^m$, as shown in Figure 2.

3.2. Variational Mode Decomposition Model

It can be seen from Figure 1 that the original passenger flow time series is non-stationary, which has an impact on the accuracy of the passenger flow prediction. To better predict the passenger flow and improve the prediction performance of the model, it is necessary to process the data first. In this paper, Variable Mode Decomposition (VMD) is used to process the original passenger flow time series, and $x_n^m$ is decomposed into multiple relatively stable sub-sequences with different frequency scales, $x_{n,1}^m,x_{n,2}^m,\dots ,x_{n,k}^m$, which is the Intrinsic Mode Function (IMF). The core idea of VMD is to construct and solve variational problems. The decomposition process is as follows [34].

• Each sub-sequence, $x_{n,k}^m$, is processed by Hilbert transform:

  $$\left(\delta t+\frac{j}{\pi t}\right)\ast x_{n,k}^m\left(t\right)$$

  (1)

  where $\delta t$ represents the Dirac distribution function, which is the pulse function, and $\ast$ represents the convolution operation.
• The transformed $x_{n,k}^m$ is multiplied by the exponential hybrid demodulation center frequency $e^{-j{\omega }_{k}t}$, and then the spectrum is transformed to a baseband:

  $$\left[\left(\delta t+\frac{j}{\pi t}\right)\ast x_{n,k}^m\left(t\right)\right]e^{-j{\omega }_{k}t}$$

  (2)
• The bandwidth of each $x_{n,k}^m$ can be obtained by Gaussian smooth estimation of the demodulated signal, which is the $L^2$ norm of the gradient, and the constrained variational mode is expressed as:

  $$\begin{array}{l}\underset{\left\{x_{n,k}^m\right\},\left\{{\omega }_k\right\}}{min}{\displaystyle {\displaystyle \sum }_k}||{\partial }_t\left[\left(\delta t+\frac{j}{\pi t}\right)\ast x_{n,k}^m\left(t\right)\right]e^{-j{\omega }_{k}t}|{|}_2^2\\ s.t. {\displaystyle {\displaystyle \sum }_k}{x}_{n,k}^m=x_n^m\end{array}$$

  (3)

  where ${\omega }_k$ is the center frequency corresponding to each mode, $k$ is the number of mode components, $x_n^m$ is the original passenger flow time series, and ${\partial }_t$ means finding the partial derivative.
• To solve the optimal solution of the constrained variational problem in step 3, the problem can be converted into unconstrained variational mode problems. A quadratic penalty factor, $\alpha$, and Lagrangian operator, $\lambda \left(t\right)$, are introduced, and an augmented Lagrangian expression in the following form is constructed:

  $$\begin{array}{r}L\left(\left\{x_{n,k}^m\right\},\left\{{\omega }_k\right\}, \lambda \right)=\alpha {\displaystyle {\displaystyle \sum }_k}||{\partial }_t[\left(\delta t+\frac{j}{\pi t}\right)\\ \ast x_{n,k}^m\left(t\right)]e^{-j{\omega }_{k}t}|{|}_2^2+\\ ||x_n^m\left(t\right)-{\displaystyle {\displaystyle \sum }_k}{x}_{n,k}^m\left(t\right)|{|}_2^2+\langle \lambda \left(t\right),x_n^m\left(t\right)-{\displaystyle {\displaystyle \sum }_k}{x}_{n,k}^m\left(t\right)\rangle \end{array}$$

  (4)
• To solve Equation (4), the saddle point of the above formula can be obtained by the alternating direction multiplier method in VMD, and by updating $x_{n,k}^m$, ${\omega }_k$, and $\lambda$ to find the optimal solution of the constrained variational mode, the updated method is as follows:

  $${\widehat{x}}_{n,k}^m{}^{c+1}\left(\omega \right)=\frac{{\widehat{x}}_n^m\left(\omega \right)-{\sum }_{i\ne k}{\widehat{x}}_{n,k}^m{}^c\left(\omega \right)+\frac{{\mathsf{\lambda }}^c\left(\omega \right)}{2}}{1+2\alpha {(\omega -{\omega }_k^c)}^2}$$

  (5)

  $${\omega }_k^{c+1}=\frac{{\int }_0^{\infty }\omega |{\widehat{x}}_{n,k}^m{}^{c+1}\left(\omega \right){|}^{2}d\omega }{{\int }_0^{\infty }|{\widehat{x}}_{n,k}^m{}^{c+1}\left(\omega \right){|}^{2}d\omega }$$

  (6)

  $${\hat{\mathsf{\lambda }}}^{c+1}\left(\omega \right)={\hat{\mathsf{\lambda }}}^c\left(\omega \right)+\tau [{\widehat{x}}_n^m\left(\omega \right)-{\displaystyle \sum }_k{\widehat{x}}_{n,k}^m{}^{c+1}\left(\omega \right)]$$

  (7)

  where ${\widehat{x}}_n^m\left(\omega \right),{\widehat{x}}_{n,k}^m\left(\omega \right),\mathrm{and} \widehat{\lambda }{\left(\omega \right)}^{}$ are the Fourier transforms of $x_n^m\left(t\right),x_{n,k}^m\left(t\right), \mathrm{and} \lambda \left(t\right)$, respectively, $\omega$ is the random frequency, $c$ is the number of iterations, and $\tau$ is the updated parameter of the Lagrange multiplier.
• Repeating step 5, the judgment condition for stopping the loop iteration is:

  $${\displaystyle \sum }_k\frac{\left|\right|{\widehat{x}}_{n,k}^m{}^{c+1}-{\widehat{x}}_{n,k}^m{}^c|{|}_2^2}{\left|\right|{\widehat{x}}_{n,k}^m{}^c|{|}_2^2}<\epsilon$$

  (8)

After VMD decomposition, the original passenger flow time series, $x_n^m$, is decomposed into $k$ IMF components, which are $x_{n,1}^m,x_{n,2}^m,\dots ,x_{n,k}^m$, and the $k$ IMF components decomposed in each time period are reconstructed separately. There are many ways of reconstruction, such as component addition, component weighted average RE addition, or high-correlation and low-correlation components are added, respectively. This paper adopts the method of component addition, which is ${x^{\prime }}_n^m=x_{n,1}^m+x_{n,2}^m+\dots +x_{n,k}^m$, the reconstructed new passenger flow time series is ${x^{\prime }}_n^m$, and then the sequence ${x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m$ is used as the input of the forecasting model.

4. Neural Network-Based Hybrid Forecasting Models

The hybrid model can integrate the advantages of multiple models and improve the accuracy of passenger flow forecasting. Based on this, this paper combined the VMD method with a single passenger flow forecasting method, and designed hybrid models to predict the time-varying passenger flow of an intercity high-speed railway, which are the VMD-MLP model, the VMD-GRU model, and the VMD-Bi-LSTM model.

4.1. Hybrid Model Design

The hybrid forecasting models mainly include two steps, as follows.

Step 1: Data preprocessing. First, the data were cleaned and screened for repeated values, abnormal values, and missing values, and then the appropriate super-parameters were selected based on the decomposition principle of VMD. The original passenger flow time series $x_n^m$ was decomposed into $k$ IMF eigenmode functions to weaken the noise interference on the forecasting models and improve the stability of the input time series, and it was finally reconstructed into a new passenger flow time series ${x^{\prime }}_n^m$.

Step 2: Passenger flow forecasting. The reconstructed time series ${x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m$ were substituted into the neural network forecasting model for prediction. First, the training time step, prediction time step, number of training samples, and division into a training set and test set of the forecasting models were determined. Second, the neural network forecasting model was constructed based on the deep learning framework Keras, and appropriate network hyperparameters were selected, such as hidden neurons and iteration times. The output value is ${\widehat{x}}_n^m$. Finally, the predicted value was compared with the real value to calculate the prediction error. There are three commonly used prediction errors, namely Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE), as shown in Formulas (9)–(11), where $s$ is the number of samples in the dataset and $M$ is the number of time periods.

$$RMSE=\frac{1}{M}{\displaystyle \sum }_{m=1}^M\sqrt{\frac{{\sum }_{n=1}^s|x_n^m-{\widehat{x}}_n^m{|}^2}{s}}$$

(9)

$$MAE=\frac{1}{M}{\displaystyle \sum }_{m=1}^M\frac{1}{s}{\displaystyle \sum }_{n=1}^s\left|x_n^m-{\widehat{x}}_n^m\right|$$

(10)

$$MAPE=\frac{1}{M}{\displaystyle \sum }_{m=1}^M\frac{1}{s}{\displaystyle \sum }_{n=1}^s|(\frac{x_n^m-{\widehat{x}}_n^m}{x_n^m})|$$

(11)

4.2. Introduction of Neural Network Forecasting Model

4.2.1. Multi-Layer Perceptron Model

Multi-Layer Perceptron (MLP) is a specific artificial neural network. To overcome the limitations of the single-layer perceptron, it was proposed to add a hidden layer between the input layer and the output layer, so as to turn a single-layer perceptron into a multi-layer perceptron, and the input signal propagates forward through the network layer-by-layer [47]. The MLP model structure in this paper is shown in Figure 3.

The input vector ${x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m$ propagates through the neurons in the input layer of the MLP to all neurons in the hidden layer, and the input and output of the $h$-th hidden neuron are $y_h$ and $Y_h$, respectively. The output of the output layer is the predicted value ${\widehat{x}}_n^m$.

$$y_h={\displaystyle \sum }_{i=1}^{14}{W}_{ih}{x^{\prime }}_{n-i}^m$$

(12)

$$Y_h=f\left(y_h\right)$$

(13)

$${\widehat{x}}_n^m=f\left({\displaystyle \sum }_{h=1}^H{Y}_h{W}_h\right)$$

(14)

where $i$ is the $i$-th neuron of the input layer, $h$ is the $h$-th neuron of the hidden layer, $W_{ih}$ represents the weight between the $i$-th neuron of the input layer and the $h$-th neuron of the hidden layer, and $W_h$ represents the weight between the $h$-th neuron of the hidden layer and the output layer.

4.2.2. Gated Recurrent Unit Neural Network Model

The Gated Recurrent Unit Neural Network (GRU) is an improved version of the Recurrent Neural Network. The GRU model has two gates, that is, a reset gate and an update gate. The reset gate controls the length of forgotten information, and the update gate can forget as well as select memory [47], as shown in Figure 4. The GRU model is composed of a simpler gate mechanism, which makes the calculation speed of learning and prediction faster.

In Figure 4, $t$ is the iteration time, $H_{t-1}$ is the hidden state at the previous time, $H_t$ is the current time output, $\sigma$ and $tanh$ are activation functions, $r_t$ is the reset gate, $z_t$ is the update gate, and ${\tilde{H}}_t$ is the hidden layer information, as follows:

$$\sigma \left(t\right)=\frac{1}{1+e^{-t}}$$

(15)

$$tanh\left(t\right)=\frac{e^t-e^{-t}}{e^t+e^{-t}}$$

(16)

The output of the function $\sigma \left(t\right)$ does not consider the information learned at the previous time. The $tanh\left(t\right)$ function compresses the previously learned information to stabilize the value. The calculation process of the GRU model includes the following key steps:

• Calculation of reset gate, $r_t$. ${x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m$ represent the inputs of the current cell. The reset gate is used to control the amount of information that needs to be forgotten about the previous cell. It will read $H_{t-1}$ and ${x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m$, and the process can be expressed as:

  $$r_t=\sigma \left(W_r\left[H_{t-1},{x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m\right]+b_r\right)$$

  (17)

  where $W_r$ is the weight matrix of the reset gate and $b_r$ is the deviation.
• Calculation of update gate, ${\tilde{H}}_t$. The update gate determines which information is to be discarded from the previous cell and which new information is to be added to the current cell in the GRU model, to reduce the risk of gradient disappearance. This process can be expressed as:

  $$z_t=\sigma \left(W_z\left[H_{t-1},{x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m\right]+b_z\right)$$

  (18)

  where $W_z$ is the weight matrix of the reset gate and $b_z$ is the deviation.
• Calculation of hidden layer information, ${\tilde{H}}_t$. The information of the cell passing through the reset gate, $r_t$, and inputs ${x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m$ together are processed by the tanh function, and the output is the hidden layer information. The process can be expressed as:

  $${\tilde{H}}_t=tanh\left(W_{\tilde{H}}\left[r_t{H}_{t-1},{x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m\right]+b_H\right)$$

  (19)

  where $W_{\tilde{H}}$ is the weight matrix and $b_H$ is the deviation.
• Calculation of output, $H_t$: First, the hidden layer state at the previous time $H_{t-1}$ multiplies the treated $z_t$, and then $z_t$ multiplies ${\tilde{H}}_t;$ finally, the two are added to obtain $H_t$. Specifically, this can be expressed as:

  $$H_t=\left(1-z_t\right)H_{t-1}+z_t{\tilde{H}}_t$$

  (20)

The above four steps are the loop body of the GRU model. The iteration times are usually determined in experiments. The output $H_t$ obtained after completing the given iteration times is the last predicted value, ${\widehat{x}}_n^m$.

4.2.3. Bi-Directional Long Short-Term Memory Model

The Bi-Directional Long Short-Term Memory Model (Bi-LSTM) is a model designed based on the LSTM model. This model can obtain the data feature information in both directions of the hidden layer in the calculation process, which helps to improve the prediction accuracy. The structure of the Bi-LSTM model is shown in Figure 5, which contains six weight matrices, named $W_1–W_6$. The forward layer performs forward calculation from time 1 to time t to obtain and save data at each time, while the backward layer reverses the calculation from time t to time 1 to obtain and save the data at each time. $h_t$ is the final output value.

The operation formulas of the Bi-LSTM model are as follows:

$$s_t=f\left(W_1\left[{x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m\right]+W_2{s}_{t-1}\right)$$

(21)

$${s_t}^{\prime }=f\left({W_3\left[{x^{\prime }}_{n-14}^m,{x^{\prime }}_{n-13}^m,\dots \dots {x^{\prime }}_{n-1}^m\right]+W_5{s}_{t-1}}^{\prime }\right)$$

(22)

$$h_t=g\left({W_4{s}_t+W_6{s}_t}^{\prime }\right)$$

(23)

5. Experimental Analysis

Based on the historical ticket sale data of the Guangzhou–Zhuhai intercity high-speed railway in China from June 2013 to June 2016, five typical O-Ds were selected for experimental analysis: Guangzhou South railway station to Zhongshan North railway station (D1), Guangzhou South railway station to Zhuhai railway station (D2), Guangzhou South railway station to Xiaolan railway station (D3), Guangzhou South railway station to Zhuhai North railway station (D4), and Zhongshan North railway station to Zhuhai North railway station (D5). Data during national legal holidays were deleted and a total of 1047 days of historical ticket sale data remained. The first 80% (838 days) of the dataset was selected as the training set, and the last 20% (209 days) was selected as the test set.

5.1. Decomposition Results of VMD Model

The VMD method was used to decompose the non-stationary original passenger flow time series $\left\{x_n^m\right\}$ into stable sub-sequences $\left\{x_{n,1}^m\right\},\left\{x_{n,2}^m\right\},\dots \dots \left\{x_{n,k}^m\right\}$ for the five selected O-Ds. When the number of modes, $k$, is small, some important information in the original passenger flow time series $\left\{x_n^m\right\}$ will be filtered, affecting the accuracy of the subsequent prediction. However, when the number of modes, $k$, is large, it will lead to repeated modes or extra noise, and excessive decomposition will lead to overlapping signal components. After many experiments, it was determined that the number of modes, $k$, of the VMD model was 10, and the output sub-sequences were $\left\{x_{n,1}^m\right\},\left\{x_{n,2}^m\right\},\dots \dots \left\{x_{n,10}^m\right\};$ that is, $IMF1,IMF2,$…$IMF10$. The bandwidth, $a$, was 7000, the noise tolerance, $tau$, was 0, the non-constant value of the synthetic signal was 0, and the control error constant, $tol$, was $e^{-7}$. The VMD model decomposed the four passenger flow time series under different time granularities for each O-D. The frequency of IMF components is listed from low to high. The decomposition results of the passenger flow time series $\left\{x_1^1,x_2^1,\dots \dots x_{1047}^1\right\}$ from Guangzhou South to Zhuhai in time period 1 are shown in Figure 6.

The low-frequency component IMF1 mainly reflects the overall change trend of the original passenger flow sequence. Several peaks in IMF2 were caused by fluctuations around holidays (excluded). The middle-frequency components from IMF3 to IMF6 exhibited certain periodic characteristics. The high-frequency components from IMF7 to IMF10 mainly reflect other fluctuation characteristics of the original passenger flow sequence, without obvious periodicity.

5.2. Prediction Results of Hybrid Model

To accurately predict the passenger flow time series, for datasets of different O-Ds, the ranges of iteration times, optimizer, hidden neurons, and hidden layers were adjusted appropriately during training, to obtain better prediction results. For the GRU and Bi-LSTM models, the parameters were set as follows: the time step was 14, the prediction step was 1, the number of iterations was in the interval [10, 100], the learning rate was 0.01, the loss function was RMSE, the optimizer was the Adam optimizer, and the number of hidden layers was 1. For the MLP model, the numbers of neurons in the input layer and the output layer were 14 and 1, respectively. The other super-parameter settings were the same as for GRU and Bi-LSTM. There are many methods to determine the number of hidden neurons in neural networks, such as the empirical method, repeated experiment method, growth method, and genetic algorithm [47]. In this paper, the empirical method was used, the number of neurons in the hidden layer was $2^i$, and $i$ was a positive integer.

Based on the above parameter settings, the hybrid models VMD-MLP, VMD-GRU, and VMD-Bi-LSTM were used to predict the time-varying passenger flow of the Guangzhou–Zhuhai intercity high-speed railway. The comparison between the prediction value (the prediction value corresponding to the hybrid model with the smallest MAPE error was selected) and the actual value from Guangzhou South station to Zhongshan North station under different time granularities every day in a week is shown in Figure 7.

In Figure 7, the distribution characteristics of the passenger flow prediction values in each period are consistent with the actual values overall, which shows that the hybrid models can better fit the time-varying characteristics of passenger travel demand and are suitable for predicting the time-varying passenger flow of the intercity high-speed railway.

Table 1. Prediction error of 1 h (RMSE/MAE/MAPE).

Model	NHN	NI	D1	D2	D3	D4	D5
		10	18.7/14.6/0.10	67.1/53.7/0.21	22.3/17.6/0.13	5.5/4.3/0.05	2.8/2.2/0.02
	64	50	21.3/16.5/0.09	60.3/48.2/0.25	19.5/15.7/0.12	5.6/4.4/0.05	1.9/1.6/0.01
		100	17.2/13.6/0.11	47.0/35.8/0.16	22.3/17.8/0.14	5.0/3.9/0.04	3.6/2.9/0.04
		10	20.2/16.0/0.11	52.1/41.0/0.19	14.4/11.2/0.08	5.3/4.2/0.05	8.5/8.2/0.08
VMD	128	50	22.6/17.8/0.09	54.7/43.2/0.22	14.8/11.7/0.09	5.0/4.0/0.05	1.9/1.5/0.01
-MLP		100	19.4/15.6/0.13	45.4/34.6/0.16	20.0/16.0/0.12	5.0/4.0/0.05	3.6/3.2/0.04
		10	19.4/15.5/0.12	54.6/43.0/0.22	16.0/12.8/0.10	5.6/4.3/0.05	4.8/3.9/0.05
	256	50	17.8/13.6/0.09	48.0/37.4/0.18	15.5/12.1/0.09	5.5/4.4/0.06	2.5/2.1/0.02
		100	17.5/13.7/0.11	55.8/43.9/0.21	32.6/28.2/0.22	5.3/4.3/0.05	4.9/4.0/0.06
		10	75.7/72.6/0.64	141.3/132.6/0.74	52.0/49.6/0.39	6.3/4.9/0.07	6.0/6.0/0.08
	64	50	32.5/26.3/0.17	66.4/53.3/0.28	104.3/102.6/0.81	13.1/12.0/0.15	3.9/3.9/0.05
		100	21.5/17.3/0.13	73.0/59.0/0.31	39.3/36.1/0.28	12.3/10.7/0.13	2.3/2.0/0.02
		10	25.6/20.4/0.14	146.2/131.3/0.68	28.3/24.1/0.19	26.3/25.6/0.32	2.4/2.1/0.03
VMD	128	50	31.4/27.4/0.27	69.6/55.0/0.27	49.0/40.3/0.32	5.9/4.7/0.06	6.3/6.2/0.08
-GRU		100	25.1/19.9/0.15	77.6/59.7/0.31	36.9/27.3/0.21	15.0/12.9/0.15	1.6/1.3/0.01
		10	19.2/14.8/0.10	113.2/97.4/0.53	22.6/19.6/0.15	20.1/19.4/0.25	1.5/1.3/0.01
	256	50	24.2/19.3/0.13	140.2/130.7/0.62	22.1/18.0/0.14	9.7/8.4/0.11	1.7/1.4/0.02
		100	24.1/18.6/0.10	73.8/59.0/0.30	29.8/25.3/0.19	8.9/7.2/0.09	6.0/5.4/0.07
		10	19.1/14.7/0.11	51.1/39.4/0.20	13.4/10.7/0.08	6.2/4.8/0.06	1.9/1.5/0.02
	64	50	22.0/16.3/0.12	57.2/44.3/0.22	7.4/3.9/0.03	6.1/4.7/0.06	1.1/0.9/0.01
VMD		100	24.4/16.3/0.13	69.5/53.6/0.31	13.3/4.2/0.03	6.1/4.7/0.06	4.8/3.0/0.04
-Bi-L		10	19.7/14.9/0.11	59.2/44.9/0.29	11.9/9.0/0.07	5.9/4.6/0.06	7.4/6.6/0.07
STM	128	50	22.3/17.1/0.14	64.9/49.6/0.22	10.2/7.3/0.05	6.3/5.0/0.06	2.2/1.8/0.02
		100	24.0/18.8/0.15	77.6/53.9/0.23	4.6/3.0/0.02	6.6/5.1/0.06	5.8/3.1/0.04
		10	19.8/15.1/0.11	56.1/43.2/0.21	14.7/10.7/0.08	6.8/5.4/0.07	5.5/2.8/0.04
	256	50	22.4/17.3/0.15	67.1/52.5/0.24	9.5/5.8/0.04	5.9/4.6/0.06	3.9/2.9/0.04
		100	24.9/19.7/0.18	62.1/48.4/0.23	13.5/8.8/0.06	6.2/4.8/0.06	9.9/5.8/0.08

Note: The minimum prediction error of each model for each O-D pair is underlined.

,

Table 2. Prediction error of 3 h (RMSE/MAE/MAPE).

Model	NHN	NI	D1	D2	D3	D4	D5
		10	21.1/15.8/0.02	170.6/130.3/0.11	31.5/24.0/0.03	14.8/4.6/0.09	2.0/1.5/0.03
	64	50	20.7/15.4/0.02	151.1/113.2/0.11	36.8/29.3/0.04	14.7/4.2/0.08	1.5/1.1/0.01
		100	21.3/16.9/0.02	153.5/115.8/0.11	34.0/26.6/0.03	14.7/4.2/0.08	1.8/1.5/0.02
		10	20.6/16.2/0.02	183.3/144.0/0.12	31.2/23.6/0.03	14.9/4.5/0.09	2.0/1.6/0.05
VMD	128	50	20.5/15.7/0.02	148.2/111.4/0.11	37.1/29.7/0.04	14.7/4.2/0.08	1.5/1.1/0.01
-MLP		100	20.9/16.5/0.02	152.6/113.9/0.10	33.5/26.2/0.03	14.7/4.2/0.08	1.8/1.4/0.02
		10	22.7/18.5/0.02	170.8/130.7/0.11	34.8/27.4/0.03	14.8/4.8/0.09	2.2/1.8/0.10
	256	50	20.9/16.4/0.02	149.8/113.0/0.11	37.1/29.7/0.04	14.7/4.2/0.09	1.5/1.1/0.01
		100	22.8/16.9/0.02	153.5/114.8/0.11	34.1/26.8/0.03	14.7/4.2/0.08	1.5/1.1/0.01
		10	30.8/25.9/0.03	507.0/476.3/0.42	50.5/40.1/0.05	15.1/6.8/0.13	5.0/4.7/0.04
	64	50	26.2/21.0/0.03	216.0/174.9/0.15	36.4/28.1/0.03	15.2/7.5/0.10	4.2/3.9/0.03
		100	37.5/29.1/0.04	321.8/260.8/0.23	35.8/28.2/0.04	14.8/5.0/0.10	2.2/1.8/0.02
		10	27.4/24.6/0.03	162.0/122.8/0.14	55.4/47.9/0.06	15.1/6.5/0.13	6.3/5.9/0.06
VMD	128	50	30.0/15.7/0.02	270.6/215.6/0.18	92.7/86.2/0.12	15.7/8.6/0.17	6.9/6.6/0.07
-GRU		100	66.2/53.4/0.08	516.3/442.7/0.49	40.6/32.2/0.04	7.1/4.5/0.40	4.9/4.5/0.03
		10	23.5/18.2/0.02	175.4/134.8/0.14	51.3/41.4/0.05	16.4/9.4/0.19	2.2/1.8/0.01
	256	50	30.6/25.1/0.03	235.7/179.8/0.15	47.4/37.8/0.05	15.3/7.8/0.15	2.6/2.1/0.01
		100	39.6/33.1/0.04	237.3/189.6/0.20	43.7/34.6/0.04	14.8/5.0/0.10	3.5/3.5/0.02
		10	24.5/18.8/0.02	165.7/126.8/0.11	30.9/23.8/0.03	14.7/4.2/0.08	2.1/1.7/0.01
	64	50	23.3/18.4/0.02	216.9/156.9/0.16	37.9/30.0/0.04	14.7/4.2/0.08	1.4/1.0/0.01
VMD		100	23.8/18.7/0.02	212.6/155.2/0.16	35.3/28.3/0.04	14.7/4.2/0.08	1.4/1.0/0.01
-Bi-L		10	25.3/20.1/0.02	176.8/137.6/0.12	30.3/23.5/0.03	14.7/4.2/0.08	1.7/1.3/0.01
STM	128	50	23.4/18.2/0.02	195.0/146.5/0.15	36.0/27.6/0.03	14.7/4.2/0.08	1.4/1.1/0.01
		100	23.9/18.5/0.02	205.8/150.0/0.15	35.9/28.0/0.04	14.7/4.2/0.08	1.3/1.0/0.01
		10	26.8/21.0/0.03	176.5/134.9/0.12	30.5/23.5/0.03	14.7/4.2/0.08	2.0/1.6/0.01
	256	50	22.4/17.5/0.02	222.0/162.4/0.15	36.9/28.2/0.04	14.7/4.2/0.08	1.3/1.0/0.01
		100	22.6/17.2/0.02	218.3/154.4/0.14	37.8/28.1/0.04	14.7/4.2/0.08	1.3/0.9/0.01

Note: The minimum prediction error of each model for each O-D pair is underlined.

,

Table 3. Prediction error of 5 h (RMSE/MAE/MAPE).

Model	NHN	NI	D1	D2	D3	D4	D5
		10	38.0/30.9/0.02	242.1/187.6/0.09	35.2/29.2/0.02	20.2/14.0/0.04	1.4/1.1/0.04
	64	50	26.8/20.6/0.01	205.7/150.2/0.07	26.0/20.1/0.01	21.2/14.0/0.04	1.2/0.8/0.04
		100	29.3/22.6/0.01	255.3/200.7/0.09	27.0/21.2/0.01	21.6/14.0/0.04	1.0/0.8/0.03
		10	42.3/35.0/0.02	251.8/197.6/0.09	35.8/29.8/0.02	20.2/13.7/0.04	1.3/1.0/0.04
VMD	128	50	26.8/20.6/0.01	204.7/148.8/0.07	27.1/21.2/0.01	21.2/13.6/0.04	1.0/0.8/0.04
-MLP		100	29.4/22.7/0.01	264.3/209.8/0.10	27.2/21.4/0.01	21.6/14.2/0.04	1.1/0.8/0.03
		10	40.5/33.4/0.02	262.7/209.5/0.10	34.1/28.0/0.02	20.6/14.1/0.04	1.5/1.2/0.06
	256	50	27.0/20.8/0.01	205.5/151.2/0.07	28.4/22.5/0.01	21.3/21.3/0.04	1.3/0.7/0.04
		100	34.1/27.2/0.02	260.2/205.1/0.10	31.9/25.2/0.01	21.4/13.8/0.04	1.1/0.8/0.03
		10	91.5/83.3/0.06	527.6/462.7/0.21	99.6/86.3/0.07	23.4/17.6/0.07	5.4/5.1/0.12
	64	50	100.0/93.5/0.07	240.9/184.2/0.08	37.5/31.5/0.02	22.1/17.5/0.07	3.5/2.1/0.11
		100	48.5/39.0/0.02	246.0/184.8/0.08	83.6/71.8/0.05	24.9/19.8/0.09	4.8/4.5/0.14
		10	84.9/77.9/0.05	817.4/778.7/0.03	36.5/30.1/0.02	34.0/28.5/0.12	10.0/9.8/0.17
VMD	128	50	54.1/44.2/0.03	271.8/204.0/0.09	38.0/28.4/0.02	44.5/41.0/0.20	2.0/1.3/0.06
-GRU		100	128/118./0.09	503.7/423.2/0.21	77.3/67.7/0.05	20.9/14.2/0.05	4.7/4.3/0.11
		10	69.4/56.1/0.04	654.9/612.2/0.28	125.6/121.6/0.09	22.7/15.8/0.05	10.0/ 9.0/0.15
	256	50	101.0/88.8/0.06	282.1/209.9/0.10	44.4/37.3/0.03	37.2/32.7/0.14	7.1/6.6/0.13
		100	48.1/39.0/0.02	452.2/364.0/0.17	77.5/69.6/0.05	36.8/32.6/0.15	2.4/2.0/0.10
		10	31.8/24.5/0.01	226.6/166.5/0.08	31.9/26.0/0.02	18.6/12.5/0.04	1.3/1.0/0.05
	64	50	40.9/30.7/0.02	330.8/225.3/0.10	28.6/22.4/0.01	18.6/12.9/0.04	1.2/0.9/0.06
		100	40.6/30.3/0.02	337.1/210.0/0.09	29.4/22.9/0.01	20.6/13.4/0.04	1.3/1.0/0.05
		10	32.0/24.0/0.01	231.0/169.2/0.10	31.9/26.0/0.02	18.7/12.8/0.04	1.2/1.0/0.05
VMD	128	50	37.4/28.2/0.02	348.7/227.7/0.10	31.6/25.4/0.02	18.5/12.2/0.04	1.0/0.8/0.05
-BiLST		100	39.7/29.1/0.02	322.4/222.1/0.10	32.3/24.8/0.01	23.2/13.6/0.04	1.1/0.8/0.05
M		10	31.7/24.3/0.01	238.1/179.0/0.08	32.3/26.2/0.02	18.8/12.7/0.04	1.2/0.9/0.05
	256	50	36.7/27.3/0.02	304.8/217.5/0.10	29.6/23.3/0.01	18.2/11.8/0.04	1.0/0.8/0.05
		100	36.7/27.8/0.02	307.5/215.0/0.10	28.7/22.6/0.01	19.4/12.4/0.04	1.1/0.9/0.05

Note: The minimum prediction error of each model for each O-D pair is underlined.

and

Table 4. Prediction error of 16 h (RMSE/MAE/MAPE).

	NHN	NI	D1	D2	D3	D4	D5
		10	310.0/242.0/0.06	711.4/532.4/0.07	295.6/217.4/0.01	45.6/31.9/0.01	13.0/9.7/0.10
	64	50	310.0/242.0/0.06	713.1/534.9/0.07	298.0/219.7/0.01	45.6/31.8/0.01	13.0/9.7/0.10
		100	310.0/242.0/0.06	713.1/534.8/0.07	298.0/219.7/0.01	45.6/31.8/0.01	13.0/9.7/0.10
		10	310.0/242.0/0.06	710.4/531.1/0.07	295.9/217.7/0.01	45.6/31.9/0.01	13.0/9.6/0.10
VMD	128	50	310.0/242.0/0.06	712.1/534.5/0.07	297.2/219.0/0.01	45.6/31.9/0.01	13.0/9.7/0.10
-MLP		100	310.0/242.0/0.06	712.1/534.4/0.07	297.2/219.0/0.01	45.6/31.9/0.01	13.0/9.7/0.10
		10	310.0/242.0/0.06	709.1/530.0/0.07	295.4/217.2/0.01	45.6/31.9/0.01	13.0/9.7/0.10
	256	50	310.0/242.0/0.06	710.5/534.0/0.07	296.2/218.1/0.01	45.6/31.9/0.01	13.0/9.7/0.09
		100	310.0/242.0/0.06	710.5/534.0/0.07	296.2/218.1/0.01	45.6/31.9/0.01	13.0/9.7/0.09
		10	608.0/388.0/0.10	885.4/717.9/0.09	1005.4/898.8/0.02	62.5/38.7/0.02	14.5/11.1/0.15
	64	50	2405.9/2308.0/0.58	3327.0/3015.7/0.41	875.0/738.8/0.01	433.3/361.1/0.15	147.5/122.1/0.34
		100	31131.9/31056.3/7.89	11591.1/10296.1/1.40	2827.3/2474.7/0.06	1072.3/1041.3/0.48	244.0/238.9/0.44
		10	866.5/770.3/0.19	3555.2/3478.1/0.46	1354.5/1289.5/0.15	749.7/743.0/0.35	73.9/71.3/0.21
VMD	128	50	13535.1/13396.5/3.44	4694.1/4311.4/0.57	5953.5/5933.7/0.15	404.7/384.3/0.18	281.7/274.6/0.68
-GRU		100	5907.2/5718.6/1.43	15964.9/15562.6/2.10	10365.2/10166.7/0.26	287.1/203.6/0.10	115.8/107.5/0.31
		10	1199.1/825.0/0.21	4870.4/4813.4/0.63	3032.8/2032.1/0.08	73.6/53.3/0.23	18.4/14.6/0.17
	256	50	991.0/790.4./ 0.21	18295.5/18103.8/2.41	1747.7/1431.1/0.03	601.8/503.4/0.25	74.4/69.2/0.23
		100	82387.7/82249.7/21.0	49438.6/49291.2/6.58	8795.4/8444.0/0.22	2914.0/2871.3/0.13	988.2/987.1/10.35
		10	246.8/182.0/0.04	605.5/437.3/0.05	278.9/193.3/0.04	44.1/30.5/0.01	11.8/8.5/0.11
	64	50	970.9/513.9/0.13	10400.8/4135.7/0.54	2638.5/703.7/0.17	67.9/43.9/0.02	16.0/12.2/0.09
		100	335.5/193.0/0.04	7510.5/6636.9/0.88	270.6/183.3/0.04	259.1/174.8/0.14	11.5/8.4/0.07
		10	245.2/183.8/0.04	756.8/603.2/0.07	292.2/196.8/0.04	43.1/29.3/0.01	12.8/9.6/0.13
VMD	128	50	246.3/186.6/0.04	712.7/531.1/0.07	11781.0/8718.8/2.25	215.6/179.8/0.08	8.9/5.8/0.06
-Bi-L		100	2337.8/1346.6./ 0.33	6339.7/5467.7/0.72	310.8/209.2/0.05	340.6/290.4/0.14	193.2/179.2/0.17
STM		10	246.0/180.5/0.04	620.0/447.7/0.05	300.1/197.1/0.05	44.8/29.3/0.01	10.8/8.7/0.11
	256	50	287.6/211.6/0.05	8333.9/8180.7/1.10	1923.3/1650.0/0.40	210.3/176.9/0.08	31.5/22.9/0.14
		100	1399.3/1290.9/0.32	8333.9/8180.7/1.10	302.1/197.1/0.05	1421.4/1349.0/0.65	473.5/467.7/0.50

Note: The minimum prediction error of each model for each O-D pair is underlined.

show the prediction errors of the three hybrid models corresponding to different parameters at four time granularities for each O-D pair, where the minimum prediction error of each model for each O-D pair is underlined, NHN represents the number of hidden neurons, and NI represents the number of iterations. For convenience, when MAPE reached the minimum, the final tuning of hyperparameters of the three forecasting models under different time granularities were as shown in

Table 5. The final tuning of hyperparameters (NHN/NI) of the three forecasting models under different time granularities.

Model	$\mathit{\Delta }\mathit{t}$	D1	D2	D3	D4	D5
	1	64/50	64/100	128/10	64/100	64/50
VMD	3	64/10	128/100	64/10	64/100	64/50
-MLP	5	64/50	64/50	256/50	64/10	128/100
	16	64/10	64/10	64/10	64/10	64/10
	1	256/10	128/50	256/50	128/50	128/100
VMD	3	128/50	128/10	64/50	64/50	256/10
-GRU	5	256/100	128/10	64/50	128/100	128/50
	16	64/10	64/10	64/50	64/10	64/10
	1	64/10	64/10	128/100	64/50	64/50
VMD	3	64/10	64/10	64/10	64/10	64/10
-Bi-	5	64/10	64/10	64/50	64/10	64/10
LSTM	16	64/10	64/10	64/10	64/10	128/50

. Figure 8 shows a comparison of the minimum values of the RMSE, MAE, and MAPE prediction errors of the three hybrid models for each O-D pair under different time granularities.

The following findings were obtained based on the above comparison results from Table 1, Table 2, Table 3, Table 4 and Table 5 and Figure 8.

• When the number of hidden neurons was constant, the error tended to decrease with the increase of the number of iterations, and when the number of iterations was constant, the error also tended to decrease with the increase of the number of hidden neurons. However, when the two increased to a certain amount, the error increased, which indicates that it is not the case that a greater number of hidden neurons and number of iterations is better. Using too few neurons or too few iterations in the hidden layer led to underfitting. However, using too many hidden neurons or too many iterations led to overfitting and increased the training time.
• The minimum MAPE error of the three hybrid models can be controlled within 10% generally, which indicates that the prediction accuracy of the three hybrid models is high. The MAPE prediction error of O-D D2 under time granularity $\Delta t=1$ was slightly higher. This can be attributed to the complexity of the time-varying characteristics of the passenger flow of O-D D2 under this time granularity, which led to the reduced applicability of the model.
• With the increase of time granularity, RMSE and MAE errors tended to gradually increase, while MAPE error tended to gradually decrease. The reasons are as follows: (i) The larger the time granularity was, the larger the passenger flow in each time period was, according to Equations (9) and (10), and then the larger the corresponding RMSE and MAE errors were. (ii) The smaller the time granularity was, the more fully the time-varying characteristics of passenger flow were reflected, and some irregular fluctuation characteristics easily occurred, which was not convenient for the extraction of hybrid prediction models, resulting in the decrease of the prediction accuracy, that is, the increase of the MAPE prediction error.
• Under different parameters, the errors of the VMD-GRU and VMD-Bi-LSTM models fluctuated greatly, and the errors of the VMD-MLP model changed relatively smoothly. For passenger flow under a smaller time granularity, the prediction error of the VMD-MLP model was smaller, while for the passenger flow under a larger time granularity, the prediction errors of the VMD-GRU and VMD-Bi-LSTM models were smaller. This shows that the applicability of the models was different under different time granularities.

To measure the prediction performance, seven baselines (VMD-LSTM, VMD-ARIMA, MLP, GRU, Bi-LSTM, LSTM, and ARIMA) were adopted to compare with the proposed models, and the results for O-D D1 and $\Delta t=1$ are shown in

Table 6. Comparison of the proposed models and baselines.

Model	RMSE	MAE	MAPE
VMD-MLP	17.2	13.6	0.09
VMD-GRU	19.2	14.8	0.10
VMD-Bi-LSTM	19.1	14.7	0.11
VMD-LSTM	20.1	15.6	0.12
VMD-ARIMA	21.2	15.4	0.15
MLP	25.5	19.1	0.17
GRU	26.9	20.2	0.18
Bi-LSTM	34	21.1	0.18
LSTM	31.0	24.2	0.22
ARIMA	33.2	22.9	0.20

. From Table 6, the prediction accuracy of the VMD-MLP, VMD-GRU, and VMD-Bi-LSTM models was obviously higher than the baselines and the hybrid models performed better than the single models. Hence, the proposed neural network-based hybrid models outperformed the existing models for predicting the time-varying passenger flow of an intercity high-speed railway.

6. Conclusions

In this paper, three neural network-based hybrid models: VMD-MLP, VMD-GRU, and VMD-Bi-LSTM, were designed to predict the time-varying passenger flow of an intercity high-speed railway. According to the historical ticket sale data of the Guangzhou–Zhuhai intercity high-speed railway in China, the time-varying characteristics of passenger travel demand were analyzed. The non-stationary passenger flow time series was decomposed into several stationary passenger flow time series by the VMD method, and three neural network forecasting models, MLP, GRU, and Bi-LSTM, were then used to predict the decomposed time series to obtain the passenger travel demand of O-D pairs at different time periods of the day under different time granularities. After experimental analysis, the following conclusions were drawn.

• The number of hidden neurons and the number of iterations of the neural network had a great impact on the prediction error of the hybrid models. Within a certain value range, with the increase of the number of hidden neurons and the number of iterations, the error tended to decrease, but when it increased to a certain extent, the error tended to increase. It is necessary to calibrate the parameters in combination with the change of prediction errors.
• The prediction accuracies of the three hybrid models were high, and the MAPE error could generally be controlled within 10%. With the increase of time granularity, RMSE and MAE errors tended to gradually increase, while the MAPE error tended to gradually decrease.
• In the optimization design of the model parameters, the errors of the VMD-MLP model fluctuated less, and it performed more smoothly than the VMD-GRU and VMD-Bi-LSTM models. The VMD-MLP model had a higher accuracy in passenger flow forecasting under a smaller time granularity, and the VMD-GRU and VMD-Bi-LSTM models had a higher accuracy in passenger flow forecasting under a larger time granularity.
• The proposed neural network-based hybrid models outperformed the existing models for predicting the time-varying passenger flow of an intercity high-speed railway. The hybrid models performed better than the single models.

In this study, only historical time series data were used to predict passenger flow, and in a real-world setting, multi-source data can be found to improve the prediction accuracy by fusing spatial and temporal features of passenger flow. Besides, only neural network-based hybrid models were proposed. There are also some other prediction models, such as the deep learning approach, which are more flexible in architecture. In the future, we plan to design a deep learning approach by considering multi-source data, such as weather, spatial structure, and temporal features, to predict the time-varying passenger flow of the high-speed railway or a large-scale metro system.