# GMM Clustering Based on WOA Optimization and Space-Time Coupled Urban Rail Traffic Flow Prediction by CEEMD-SE-BiGRU-AM

Authors: Qiong Jiang
Journal: Mobile Information Systems
 (2022)

---

## Abstract

Accurately predicting the short-term passenger flow of urban subways is very important for urban subway stations to formulate passenger flow organization and evacuation plans effectively and rationally plan passenger travel routes. This paper establishes a novel framework to predict the hourly inbound and outbound passenger flow of subway stations based on WOA-GMM station classifiers, CEEMD-SE noise reduction, and BiGRU optimized by attention. Firstly, this paper classifies subway stations using the improved Gaussian mixture model (GMM) with Whale Optimization Algorithm (WOA) to realize the feature extraction of different types of subway stations. Secondly, this paper uses the Complementary Ensemble Empirical Mode Decomposition (CEEMD) to decompose the noise reduction of each station’s hourly inbound and outbound passenger flow. It combines the empirical modal components by calculating the sample entropy (SE), which makes the time series stable and reduces the time cost of forecasting. Finally, a Bidirectional Gated Recurrent Unit (BiGRU) model improved by attention mechanism (AM) is established for each station’s inbound and outbound passenger flow. The prediction model established in this paper is verified by subway passenger flow data in Shanghai, China, within 24 days. Finally, it is concluded that the model can predict the passenger flow of subway stations. Compared with the traditional Backpropagation Neural Network (BP), the Long Short-Term Memory (LSTM), and the normal BiGRU model, the model proposed in this paper has an average reduction of 65.90%, 64.54%, and 49.06% in Mean Absolute Percentage Error (MAPE) in the prediction of the hourly inbound and outbound passenger flow of each type of station, respectively.


---

## Body

## 1. Introduction

The ability to accurately forecast the short-term passenger flow of city subways is critical to inhabitants’ quality of life. With the rapid advancement of global urbanization, many urban rail transit systems have achieved explosive growth [1]. The number of newly opened transportation systems has surged during the past decade. In 2019 alone, the total number of passengers transported by subway and light rail globally reached 70.794 billion. By the end of 2020, 538 cities in 77 nations and regions have inaugurated urban rail transportation, with a total operating mileage of over 33,346 km and over 34,220 stations [2]. The rapid expansion of the rail transit network has also caused passenger flow organization and train scheduling problems. Excessive passenger flow on the subway and excessive inbound and outbound passenger flow of the station may cause safety problems such as stampede incidents and are not conducive to passengers reaching their destination according to the target time. However, due to the different external characteristics such as spatial layout and commercial effect of each station, excessive attention to these characteristics will complicate the prediction of passenger flow. Therefore, making the best use of the time and space features of subway stations and accurately anticipating the passenger flow of rail transit stations is a pressing issue that must be addressed.


## 2. Literature Review

Considering that the geographical location of different stations and the surrounding economic development conditions are different, the passenger flow between adjacent stations greatly correlates. Therefore, it is necessary to perform cluster analysis on stations from time and space perspectives. Zhou et al. [2] thoroughly considered each station’s topology structure and passenger flow information, which used complex network evaluation and the k-means method to cluster the stations. It proposed a method to balance and control passenger flow. However, in practical applications, the k-means algorithm [3], Affinity Propagation (AP) algorithm [4], Spectral Clustering (SC) algorithm [5], and other methods are more suitable for processing spherical data, and the processing effect for data with irregular shapes is not satisfactory. In contrast, the GMM algorithm assumes that a combination of multiple Gaussian distributions forms the data points, and there is no restriction on the shape of the data. Multiple Gaussian distribution functions can be used to fit the shape of the data points to the greatest extent. However, since the GMM algorithm is based on the Expectation-Maximum (EM) algorithm, it is simple to attain a local optimum and it is extremely sensitive to the initial value. For such a problem, Guo et al. [7] presented utilizing the Particle Swarm Optimization (PSO) optimized GMM as a correction approach for the parameters of the electromagnetic transient model. It also verifies the effectiveness of optimizing the initial parameters of the GMM algorithm through the heuristic algorithm. WOA [8] is a new variant of the heuristic optimization algorithm that outperforms the PSO and Grey Wolf Optimizer (GWO) algorithm [9, 10] in terms of performance. Abbas et al. [11] used the WOA to extract essential features in the breast cancer dataset, which greatly improved the prediction accuracy. Gadekallu et al. [7] used the combination of Principal Component Analysis (PCA) and WOA to accurately extract the important features of tomato disease data and used a deep neural network to identify the disease types. Among them, using the WOA improves the model’s accuracy by 4% on the test set. Yan et al. [13] used WOA to improve the extreme gradient boosting (XGB) model and found that the optimized model predicted daily reference evapotranspiration considerably better than the basic FAO-56 PM model. To overcome the problems, the parameters of GMM based on the EM algorithm attain the local optimum and rely too much on the selection of the initial value. The WOA is used to optimize the initial parameters of GMM in this study. At the same time, it also solves the problem of low fitting accuracy and robustness of the k-means algorithm.In addition, changes in the external environment will have a particular impact on passenger flow data. The data noise created in this manner is unpredictable, and too much random fluctuation will diminish the predictability of the results. Therefore, it is necessary to decompose and denoise the original data. In recent years, researchers have developed various processing methods to identify the characteristics of nonlinear sequences to maximize the universality of the prediction results and reduce the impact of external environmental factors on the results. These include wavelet transform (WT) [14], Singular Spectrum Analysis (SSA) [15], Empirical Mode Decomposition (EMD) [16], Ensemble Empirical Mode Decomposition (EEMD) [17], etc. Some researchers have also proved the higher accuracy of the hybrid prediction model with decomposition and noise reduction processing. Zhu et al. [18] used the Autoregressive Integrated Moving Average model (ARIMA) to predict the passenger flow and combined WT to decompose and denoise the data, which achieved higher accuracy than the ARIMA model alone. However, the wavelet transform needs to specify the wavelet basis function artificially, and the decomposition function lacks directionality [19]. In contrast, the empirical mode decomposition method is more adaptive. Zhao et al. [20] used the empirical mode decomposition method to convert the time series into a series of eigenmode functions and residuals. They used a long and short-term memory neural network to estimate subway passenger flow, and the results were more accurate than previous predictions. However, EMD is prone to modal aliasing, and CEEMD [21] solves this problem by introducing complementary noise. Therefore, this paper uses CEEMD with high computational efficiency and high reconstruction accuracy of the original signal to denoise the original data.Parametric and nonparametric models are the two types of traditional short-term passenger flow prediction models. Therefore, the traditional short-term passenger flow prediction models are mainly divided into parametric and nonparametric models. One of the representative parametric models is ARIMA [22]. ARIMA considers the periodic characteristics of time, and the prediction method is direct and does not require too much preprocessing. However, ARIMA ignores the influence of randomness on the overall forecasting results, and the forecasting accuracy of time series with strong randomness is not high. Compared with parametric models, nonparametric models are more flexible and can deal with more complex situations, including Support Vector Machine (SVM) [23], LSTM [24], and Gated Recurrent Unit (GRU) [25]. The WT-SVM model proposed by Sun et al. [26] complemented the advantages of WT and SVM to predict the subway passenger flow. Therefore, SVM still has a big limitation in kernel selection. Yang et al. [27] proposed the Wave-LSTM model, which used the method of controlling variables to determine the model parameters and verified the effectiveness and prediction accuracy of the hybrid model by comparing it with traditional models. In contrast, Deep Neural Networks (DNNs) [28] have a greater advantage in sequence learning, so the improved LSTM model captures the long-term time dependence of sequences mainly through gated recursive units. [29], used LSTM and GRU to forecast the short-term passenger flow of urban rail transit and found that 5 minutes was the best temporal granularity for both models to predict short-term passenger flow. Under this time granularity, the overall performance of GRU is better than that of LSTM. However, it is often difficult for recurrent neural network models to capture local features accurately. In this regard, Wu and Tan [30] used a combination of CNN and LSTM to predict traffic flow, but this method would increase the difficulty of training and be less efficient for long-term sequence prediction. The AM [31] method extracts valuable features globally. Reference [32] improved the CNN model mainly used for local feature perception by using the features with a strong global perception of AM and effectively extracting global features. Therefore, this paper adopts the BiGRU model improved by AM to make the final passenger flow prediction, which greatly improves information utilization and accuracy.


## 3. Paper Contribution

This paper establishes the WOA-GMM-CEEMD-SE-BiGRU-AM model to predict the short-term passenger flow of urban rail transit and compares it with other models to verify its accuracy and effectiveness. The contributions of this paper are as follows:(1)
Since different stations have different spatiotemporal characteristics, this paper employs GMM to group stations with comparable characteristics. However, the EM algorithm is the core part of the GMM, and its clustering effect largely depends on the initial value of the EM algorithm. Therefore, the initial parameters of the GMM are optimized by WOA. The stations are finally classified according to their spatiotemporal characteristics.(2)
The passenger flow data of rail transit fluctuates greatly, and the regularity is poor. Therefore, this paper decomposes and denoises the passenger flow data of each type of station through CEEMD and uses a more regular component as the target variable of prediction to increase the prediction accuracy. However, considering the time cost of prediction, similar variables are combined by calculating SE.(3)
In order to predict the passenger flow of rail transit, this paper firstly establishes the BiGRU, but to extract the local features of the time series more accurately, AM is used to improve the BiGRU. Finally, compared with other models, it can be verified that the model established in this paper has high accuracy.The rest of the chapters are arranged as follows: Chapter 4 introduces the model, 4.1 introduces the processing flow of the station clustering algorithm, 4.2 introduces the CEEMD-SE method for decomposing noise reduction, and 4.3 introduces the building process of the BiGRU-AM model. The third chapter carries on the calculation of the example and the validation of the model, and the fourth chapter concludes this paper 4.


## 4. Model Introduction

### 4.1. Clustering Model Based on WOA-GMM

#### 4.1.1. Whale Optimization Algorithm

Since GMM is classified based on the EM algorithm, this model can easily achieve local convergence [33]. Therefore, this paper uses WOA to solve this problem. The whale optimization algorithm is a metaheuristic algorithm based on swarm intelligence. By simulating the predation behavior of humpback whales, the random learning strategy is used to achieve global shrinkage, and the convergence rate is fast, which is suitable for solving global optimization problems. The predation process of the humpback whale is as follows:Step1: Determine the prey location and select the hunting method.(1)
Shrinkage and encirclement mechanism:When the humpback whale surrounds the prey, it will choose to swim toward the humpback whale with the best position, identify the position of the prey, and surround it. Such hunting is realized by reducing the value ofa.(1)D′=CXpt−XtXt+1=Xpt−AD′,r<0.5A=2ar−a,C=2r,a=2−2t/tmax,A<1whereXt is the position of the t-th iteration of the humpback whale. Xpt is the optimal position currently searched; t is the current iteration number; D is the distance from the humpback whale to the prey. tmax is the maximum number of iterations; a is a value linearly decreasing from 2 to 0; r is a random vector in the range 0 to 1.(2)
Spiral spit mechanismCalculate the distance between itself and the optimal individual in the current population, and spit out bubbles while swimming in a spiral.(2)D″=Xpt−XtXt+1=D″eblcos2πl+Xpt,r≥0.5b is the logarithmic spiral shape constant; l is a random number in [−1, 1].Step 2: Search for prey.When humpback whales search for prey, they will randomly swim according to the position of the same species to achieve a global search.(3)D=CXrandt−XtXt+1=Xrandt−AD,A≥1.Xrandt is a random position.


#### 4.1.2. Gaussian Mixture Model

The Gaussian mixture model (GMM) [34] assumes that all data points come from different distributions, and by finding a mixture of multidimensional Gaussian probability distributions, it eventually generates a distribution that can satisfy any shape. Data points are judged to be from the same distribution in one group. Compared with the k-means algorithm to calculate the Euclidean distance between data points, the GMM algorithm can measure clusters’ distribution probability, which makes the fitting accuracy higher. In addition, GMM has lower time complexity than traditional clustering algorithms and is more in line with the central limit theorem under the condition of large samples.In the multidimensional case, the objective function of GMM is(4)L=logpX=∑i=1Nlog∑k=1Kpxi|μk,∑kpμk,∑kpxi|μk,∑k=12πn/2∑k1/2exp−12xi−μkT∑k−1xi−μk∼Nμk,∑k.pX is the product of the likelihood functions of all Gaussian distributions. xi is the i-th sample (dimension p); pxi|μk,∑k is the conditional distribution under the specified parameters, which obeys the expectation μk (dimension p), and the covariance matrix is ∑k (dimension p × p) which is usually distributed; L is the log-likelihood function; K is the number of clusters; N is the number of samples.The calculation steps of GMM are as follows:Step 1: Set the number of sample clusters to K, and the samples obey a Gaussian distribution:(5)pxi=∑i=1Kαipxi|μi,∑i,wherepxi is a multivariate Gaussian distribution; αi is the probability of the i-th mixture component.Step 2: Calculate the probability ofxj:(6)γij=αipxj|μi,Σi∑l=1kαlpxj|μl,Σl,whereγij represents the probability that the j-th sample belongs to the i-th class.Step 3: Compute the maximum likelihood estimate for each cluster to update the model parameters, including the expectation, variance, and probability of each data belonging to class i of Gaussian mixture distribution.(7)μi′=∑j=1mγijxj∑j=1mγij,Σi′=∑j=1mγijxj−μi′xj−μi′T∑j=1mγij,αi′=∑j=1mγijm.Step 4: Repeat steps 2 and 3 until the objective function converges.Step 5: Thej-th template is divided into the corresponding category according to the rule of taking the maximum value of γij.


#### 4.1.3. WOA-GMM

GMM is established based on the EM algorithm. The EM algorithm is an expectation-maximization algorithm, but it is easy to fall into a local optimum, and the result strongly depends on the initial parameters. Therefore, this paper uses the WOA to initialize the parameters of the EM algorithm and then establishes the GMM.Step 1: Determine the number of clusters.Step 2: Construct a Gaussian mixture model objective function.Step 3: WOA is used to solve the maximum likelihood estimation of the Gaussian mixture model constructed in Step 2.Step 4: The optimization parameters obtained in Step 3 are used as the initial parameters of the EM algorithm.Step 5: Iterate until optimal.The model flow is shown in Figure1. t1 and itermax1 are the current and maximum iterations of the WOA. t2 and itermax2 are the current number of iterations and the maximum number of iterations of the GMM algorithm.Figure 1
Flowchart of station clustering model.


### 4.2. Data Noise Reduction Model Based on CEEMD-SE

#### 4.2.1. Complementary Ensemble Empirical Mode Decomposition

EMD is a decomposition noise reduction method used to smooth and linearize data. However, due to the discontinuity of the Intrinsic Mode Function (IMF) components, the method will produce modal aliasing during the decomposition process, and the improved EEMD model solves this problem. At the same time, due to the white noise introduced by EEMD, its residual parts will affect the original data to a large extent [35]. As a result, this work uses the CEEMD model, including complimentary noise to improve the computational efficiency and signal reconstruction accuracy. The calculation process of the CEEMD model is as follows:Step 1: A set of positive and negative Gaussian white noise sequencesϵi+t and ϵi−t is added, which obtains a new passenger flow sequence based on the original passenger flow sequence xt.(8)xi+t=xt+ϵi+t,xi−t=xt+ϵi−t.xi+t and xi−t represent the result of adding positive and negative white noise to the original passenger flow sequence.Step 2: Perform empirical mode decomposition on the sequences in Step 1 separately to obtain n-1 IMF components and one residual component.(9)xi+t=∑j=1ncij+t,xi−t=∑j=1ncij−t,cij+t and cij−t are, respectively, the j-th modal component obtained by empirical modal decomposition after adding Gaussian white noise for the i-th time; when j = n, it is the residual component.Step 3: Replace the Gaussian white noise sequence and repeat Step 1 and Step 2 until N groups of eigenmode components are obtained.Step 4: Calculate the average value of N groups of IMF components in Step 3 to obtain the final decomposition result.(10)cit=12N∑j=1Ncij+t+cij−t.cit is the final i-th decomposition result.


#### 4.2.2. Sample Entropy

After CEEMD decomposes the time series data, multiple IMF components can be obtained. SE [36] can measure the complexity of different IMF components, which improves the Approximate Entropy (AE), reduces the dependence on the length of the time series in the calculation process, and reduces the impact of the loss of original data on the complexity calculation. The sample entropy calculation steps are as follows:Step1: Let a particular IMF component obtained by CEEMD decomposition beX=x1,x2,…,xN. Take all consecutive u-dimensional vectors Xi=xi,xi+1,…,xi+u−1,i=1,2,…,N−u+1 in X.Step 2: Define the distance betweenXi and Xj as(11)dXi,Xj=maxk∈0,u−1xi+k−xj+k.Step 3: Given a similarity tolerancer (r > 0), count the number of dXi,Xj<r for each vector Xi, denoted by Gi, and calculate its ratio to N-u expressed as Giur.(12)Giur=1N−uGi.Step 4: Calculate the mean of allGiur expressed as Gur.(13)Gur=1N−u+1∑i=1N−u+1Giur.Step 5: Increase theu dimension to u + 1 dimension, repeat Step1∼Step4, and get Gu+1r.Step 6: Calculate SE.(14)SampEnu,r=limN⟶∞−lnGu+1rGur.However, when SE is applied to actual data, N is always finite, thus giving an estimate of sample entropy as(15)SampEnu,r,N=−lnGu+1rGur.


### 4.3. Predictive Model (BiGRU-AM)

#### 4.3.1. Bidirectional Gated Recurrent Unit

GRU is a gated recurrent neural network model that can solve problems such as long-term memory and gradient in backpropagation and is suitable for processing time series data. Compared with LSTM, the GRU model replaces the input gate and forget gate functions with only the update gate. The formula of the hidden layer of the model is as follows:(16)zt=sigmoidwzht−1,xt+bz,rt=sigmoidwrht−1,xt+br,ht˜=tanhwhrtht−1,xt+bh,ht=1−ztht−1+ztht˜,where zt is the update gate. rt is the reset gate; ht˜ is the newly generated information after processing the historical information through the update gate. wz, wr, and wh are the weight parameters of the GRU network; bz, br, bh are the thresholds of the GRU network.The GRU structure diagram is shown in Figure2.Figure 2
GRU structure diagram.The hidden state of the GRU model layer only considers the impact of historical information on the current input, and the impact of future information is equally important, so this paper uses BiGRU to solve this problem. BiGRU consists of two unidirectional and opposite GRUs. The model integrates the information before and after the current state, and the final output is the weighted summation of the output results of the forward GRU and the backward GRU and is linearly superimposed with the threshold. The main calculation formula of BiGRU is(17)hft=GRUht−1,xt,hbt=GRUht+1,xt,ht=αhft+βhbt+γ,where hft and hbt are the forward and backward GRU output vectors for the current state, and β are the weight parameters of the forward and backward GRU; γ is the linear threshold. The BiGRU network structure is shown in Figure 3.Figure 3
BiGRU network structure diagram.


#### 4.3.2. Attention Mechanism

Although BiGRU improves the efficiency of information utilization to a certain extent, it still has the problem that local features cannot be accurately extracted when dealing with long-term sequence data. Therefore, when building the model, this paper adds an attention layer to the hidden layer to solve this problem. The attention layer further extracts the key information of the time series by calculating the weight of the output result in the BiGRU layer, thereby reducing the length of the output result and improving the robustness of the model:(18)vit=ustanhwhit+b,αit=evit∑i=1nevit,ct=∑i=1nαithit,vit is the unnormalized attention score. The higher the attention score, the better the match between the input vector and the target vector; us is the randomly initialized attention matrix; w is the weight coefficient; b is the bias coefficient; ct is the attention vector; hit is the value represented by the i-th output result at the t-th time.


## 4.1. Clustering Model Based on WOA-GMM

### 4.1.1. Whale Optimization Algorithm

Since GMM is classified based on the EM algorithm, this model can easily achieve local convergence [33]. Therefore, this paper uses WOA to solve this problem. The whale optimization algorithm is a metaheuristic algorithm based on swarm intelligence. By simulating the predation behavior of humpback whales, the random learning strategy is used to achieve global shrinkage, and the convergence rate is fast, which is suitable for solving global optimization problems. The predation process of the humpback whale is as follows:Step1: Determine the prey location and select the hunting method.(1)
Shrinkage and encirclement mechanism:When the humpback whale surrounds the prey, it will choose to swim toward the humpback whale with the best position, identify the position of the prey, and surround it. Such hunting is realized by reducing the value ofa.(1)D′=CXpt−XtXt+1=Xpt−AD′,r<0.5A=2ar−a,C=2r,a=2−2t/tmax,A<1whereXt is the position of the t-th iteration of the humpback whale. Xpt is the optimal position currently searched; t is the current iteration number; D is the distance from the humpback whale to the prey. tmax is the maximum number of iterations; a is a value linearly decreasing from 2 to 0; r is a random vector in the range 0 to 1.(2)
Spiral spit mechanismCalculate the distance between itself and the optimal individual in the current population, and spit out bubbles while swimming in a spiral.(2)D″=Xpt−XtXt+1=D″eblcos2πl+Xpt,r≥0.5b is the logarithmic spiral shape constant; l is a random number in [−1, 1].Step 2: Search for prey.When humpback whales search for prey, they will randomly swim according to the position of the same species to achieve a global search.(3)D=CXrandt−XtXt+1=Xrandt−AD,A≥1.Xrandt is a random position.


### 4.1.2. Gaussian Mixture Model

The Gaussian mixture model (GMM) [34] assumes that all data points come from different distributions, and by finding a mixture of multidimensional Gaussian probability distributions, it eventually generates a distribution that can satisfy any shape. Data points are judged to be from the same distribution in one group. Compared with the k-means algorithm to calculate the Euclidean distance between data points, the GMM algorithm can measure clusters’ distribution probability, which makes the fitting accuracy higher. In addition, GMM has lower time complexity than traditional clustering algorithms and is more in line with the central limit theorem under the condition of large samples.In the multidimensional case, the objective function of GMM is(4)L=logpX=∑i=1Nlog∑k=1Kpxi|μk,∑kpμk,∑kpxi|μk,∑k=12πn/2∑k1/2exp−12xi−μkT∑k−1xi−μk∼Nμk,∑k.pX is the product of the likelihood functions of all Gaussian distributions. xi is the i-th sample (dimension p); pxi|μk,∑k is the conditional distribution under the specified parameters, which obeys the expectation μk (dimension p), and the covariance matrix is ∑k (dimension p × p) which is usually distributed; L is the log-likelihood function; K is the number of clusters; N is the number of samples.The calculation steps of GMM are as follows:Step 1: Set the number of sample clusters to K, and the samples obey a Gaussian distribution:(5)pxi=∑i=1Kαipxi|μi,∑i,wherepxi is a multivariate Gaussian distribution; αi is the probability of the i-th mixture component.Step 2: Calculate the probability ofxj:(6)γij=αipxj|μi,Σi∑l=1kαlpxj|μl,Σl,whereγij represents the probability that the j-th sample belongs to the i-th class.Step 3: Compute the maximum likelihood estimate for each cluster to update the model parameters, including the expectation, variance, and probability of each data belonging to class i of Gaussian mixture distribution.(7)μi′=∑j=1mγijxj∑j=1mγij,Σi′=∑j=1mγijxj−μi′xj−μi′T∑j=1mγij,αi′=∑j=1mγijm.Step 4: Repeat steps 2 and 3 until the objective function converges.Step 5: Thej-th template is divided into the corresponding category according to the rule of taking the maximum value of γij.


### 4.1.3. WOA-GMM

GMM is established based on the EM algorithm. The EM algorithm is an expectation-maximization algorithm, but it is easy to fall into a local optimum, and the result strongly depends on the initial parameters. Therefore, this paper uses the WOA to initialize the parameters of the EM algorithm and then establishes the GMM.Step 1: Determine the number of clusters.Step 2: Construct a Gaussian mixture model objective function.Step 3: WOA is used to solve the maximum likelihood estimation of the Gaussian mixture model constructed in Step 2.Step 4: The optimization parameters obtained in Step 3 are used as the initial parameters of the EM algorithm.Step 5: Iterate until optimal.The model flow is shown in Figure1. t1 and itermax1 are the current and maximum iterations of the WOA. t2 and itermax2 are the current number of iterations and the maximum number of iterations of the GMM algorithm.Figure 1
Flowchart of station clustering model.


## 4.1.1. Whale Optimization Algorithm

Since GMM is classified based on the EM algorithm, this model can easily achieve local convergence [33]. Therefore, this paper uses WOA to solve this problem. The whale optimization algorithm is a metaheuristic algorithm based on swarm intelligence. By simulating the predation behavior of humpback whales, the random learning strategy is used to achieve global shrinkage, and the convergence rate is fast, which is suitable for solving global optimization problems. The predation process of the humpback whale is as follows:Step1: Determine the prey location and select the hunting method.(1)
Shrinkage and encirclement mechanism:When the humpback whale surrounds the prey, it will choose to swim toward the humpback whale with the best position, identify the position of the prey, and surround it. Such hunting is realized by reducing the value ofa.(1)D′=CXpt−XtXt+1=Xpt−AD′,r<0.5A=2ar−a,C=2r,a=2−2t/tmax,A<1whereXt is the position of the t-th iteration of the humpback whale. Xpt is the optimal position currently searched; t is the current iteration number; D is the distance from the humpback whale to the prey. tmax is the maximum number of iterations; a is a value linearly decreasing from 2 to 0; r is a random vector in the range 0 to 1.(2)
Spiral spit mechanismCalculate the distance between itself and the optimal individual in the current population, and spit out bubbles while swimming in a spiral.(2)D″=Xpt−XtXt+1=D″eblcos2πl+Xpt,r≥0.5b is the logarithmic spiral shape constant; l is a random number in [−1, 1].Step 2: Search for prey.When humpback whales search for prey, they will randomly swim according to the position of the same species to achieve a global search.(3)D=CXrandt−XtXt+1=Xrandt−AD,A≥1.Xrandt is a random position.


## 4.1.2. Gaussian Mixture Model

The Gaussian mixture model (GMM) [34] assumes that all data points come from different distributions, and by finding a mixture of multidimensional Gaussian probability distributions, it eventually generates a distribution that can satisfy any shape. Data points are judged to be from the same distribution in one group. Compared with the k-means algorithm to calculate the Euclidean distance between data points, the GMM algorithm can measure clusters’ distribution probability, which makes the fitting accuracy higher. In addition, GMM has lower time complexity than traditional clustering algorithms and is more in line with the central limit theorem under the condition of large samples.In the multidimensional case, the objective function of GMM is(4)L=logpX=∑i=1Nlog∑k=1Kpxi|μk,∑kpμk,∑kpxi|μk,∑k=12πn/2∑k1/2exp−12xi−μkT∑k−1xi−μk∼Nμk,∑k.pX is the product of the likelihood functions of all Gaussian distributions. xi is the i-th sample (dimension p); pxi|μk,∑k is the conditional distribution under the specified parameters, which obeys the expectation μk (dimension p), and the covariance matrix is ∑k (dimension p × p) which is usually distributed; L is the log-likelihood function; K is the number of clusters; N is the number of samples.The calculation steps of GMM are as follows:Step 1: Set the number of sample clusters to K, and the samples obey a Gaussian distribution:(5)pxi=∑i=1Kαipxi|μi,∑i,wherepxi is a multivariate Gaussian distribution; αi is the probability of the i-th mixture component.Step 2: Calculate the probability ofxj:(6)γij=αipxj|μi,Σi∑l=1kαlpxj|μl,Σl,whereγij represents the probability that the j-th sample belongs to the i-th class.Step 3: Compute the maximum likelihood estimate for each cluster to update the model parameters, including the expectation, variance, and probability of each data belonging to class i of Gaussian mixture distribution.(7)μi′=∑j=1mγijxj∑j=1mγij,Σi′=∑j=1mγijxj−μi′xj−μi′T∑j=1mγij,αi′=∑j=1mγijm.Step 4: Repeat steps 2 and 3 until the objective function converges.Step 5: Thej-th template is divided into the corresponding category according to the rule of taking the maximum value of γij.


## 4.1.3. WOA-GMM

GMM is established based on the EM algorithm. The EM algorithm is an expectation-maximization algorithm, but it is easy to fall into a local optimum, and the result strongly depends on the initial parameters. Therefore, this paper uses the WOA to initialize the parameters of the EM algorithm and then establishes the GMM.Step 1: Determine the number of clusters.Step 2: Construct a Gaussian mixture model objective function.Step 3: WOA is used to solve the maximum likelihood estimation of the Gaussian mixture model constructed in Step 2.Step 4: The optimization parameters obtained in Step 3 are used as the initial parameters of the EM algorithm.Step 5: Iterate until optimal.The model flow is shown in Figure1. t1 and itermax1 are the current and maximum iterations of the WOA. t2 and itermax2 are the current number of iterations and the maximum number of iterations of the GMM algorithm.Figure 1
Flowchart of station clustering model.


## 4.2. Data Noise Reduction Model Based on CEEMD-SE

### 4.2.1. Complementary Ensemble Empirical Mode Decomposition

EMD is a decomposition noise reduction method used to smooth and linearize data. However, due to the discontinuity of the Intrinsic Mode Function (IMF) components, the method will produce modal aliasing during the decomposition process, and the improved EEMD model solves this problem. At the same time, due to the white noise introduced by EEMD, its residual parts will affect the original data to a large extent [35]. As a result, this work uses the CEEMD model, including complimentary noise to improve the computational efficiency and signal reconstruction accuracy. The calculation process of the CEEMD model is as follows:Step 1: A set of positive and negative Gaussian white noise sequencesϵi+t and ϵi−t is added, which obtains a new passenger flow sequence based on the original passenger flow sequence xt.(8)xi+t=xt+ϵi+t,xi−t=xt+ϵi−t.xi+t and xi−t represent the result of adding positive and negative white noise to the original passenger flow sequence.Step 2: Perform empirical mode decomposition on the sequences in Step 1 separately to obtain n-1 IMF components and one residual component.(9)xi+t=∑j=1ncij+t,xi−t=∑j=1ncij−t,cij+t and cij−t are, respectively, the j-th modal component obtained by empirical modal decomposition after adding Gaussian white noise for the i-th time; when j = n, it is the residual component.Step 3: Replace the Gaussian white noise sequence and repeat Step 1 and Step 2 until N groups of eigenmode components are obtained.Step 4: Calculate the average value of N groups of IMF components in Step 3 to obtain the final decomposition result.(10)cit=12N∑j=1Ncij+t+cij−t.cit is the final i-th decomposition result.


### 4.2.2. Sample Entropy

After CEEMD decomposes the time series data, multiple IMF components can be obtained. SE [36] can measure the complexity of different IMF components, which improves the Approximate Entropy (AE), reduces the dependence on the length of the time series in the calculation process, and reduces the impact of the loss of original data on the complexity calculation. The sample entropy calculation steps are as follows:Step1: Let a particular IMF component obtained by CEEMD decomposition beX=x1,x2,…,xN. Take all consecutive u-dimensional vectors Xi=xi,xi+1,…,xi+u−1,i=1,2,…,N−u+1 in X.Step 2: Define the distance betweenXi and Xj as(11)dXi,Xj=maxk∈0,u−1xi+k−xj+k.Step 3: Given a similarity tolerancer (r > 0), count the number of dXi,Xj<r for each vector Xi, denoted by Gi, and calculate its ratio to N-u expressed as Giur.(12)Giur=1N−uGi.Step 4: Calculate the mean of allGiur expressed as Gur.(13)Gur=1N−u+1∑i=1N−u+1Giur.Step 5: Increase theu dimension to u + 1 dimension, repeat Step1∼Step4, and get Gu+1r.Step 6: Calculate SE.(14)SampEnu,r=limN⟶∞−lnGu+1rGur.However, when SE is applied to actual data, N is always finite, thus giving an estimate of sample entropy as(15)SampEnu,r,N=−lnGu+1rGur.


## 4.2.1. Complementary Ensemble Empirical Mode Decomposition

EMD is a decomposition noise reduction method used to smooth and linearize data. However, due to the discontinuity of the Intrinsic Mode Function (IMF) components, the method will produce modal aliasing during the decomposition process, and the improved EEMD model solves this problem. At the same time, due to the white noise introduced by EEMD, its residual parts will affect the original data to a large extent [35]. As a result, this work uses the CEEMD model, including complimentary noise to improve the computational efficiency and signal reconstruction accuracy. The calculation process of the CEEMD model is as follows:Step 1: A set of positive and negative Gaussian white noise sequencesϵi+t and ϵi−t is added, which obtains a new passenger flow sequence based on the original passenger flow sequence xt.(8)xi+t=xt+ϵi+t,xi−t=xt+ϵi−t.xi+t and xi−t represent the result of adding positive and negative white noise to the original passenger flow sequence.Step 2: Perform empirical mode decomposition on the sequences in Step 1 separately to obtain n-1 IMF components and one residual component.(9)xi+t=∑j=1ncij+t,xi−t=∑j=1ncij−t,cij+t and cij−t are, respectively, the j-th modal component obtained by empirical modal decomposition after adding Gaussian white noise for the i-th time; when j = n, it is the residual component.Step 3: Replace the Gaussian white noise sequence and repeat Step 1 and Step 2 until N groups of eigenmode components are obtained.Step 4: Calculate the average value of N groups of IMF components in Step 3 to obtain the final decomposition result.(10)cit=12N∑j=1Ncij+t+cij−t.cit is the final i-th decomposition result.


## 4.2.2. Sample Entropy

After CEEMD decomposes the time series data, multiple IMF components can be obtained. SE [36] can measure the complexity of different IMF components, which improves the Approximate Entropy (AE), reduces the dependence on the length of the time series in the calculation process, and reduces the impact of the loss of original data on the complexity calculation. The sample entropy calculation steps are as follows:Step1: Let a particular IMF component obtained by CEEMD decomposition beX=x1,x2,…,xN. Take all consecutive u-dimensional vectors Xi=xi,xi+1,…,xi+u−1,i=1,2,…,N−u+1 in X.Step 2: Define the distance betweenXi and Xj as(11)dXi,Xj=maxk∈0,u−1xi+k−xj+k.Step 3: Given a similarity tolerancer (r > 0), count the number of dXi,Xj<r for each vector Xi, denoted by Gi, and calculate its ratio to N-u expressed as Giur.(12)Giur=1N−uGi.Step 4: Calculate the mean of allGiur expressed as Gur.(13)Gur=1N−u+1∑i=1N−u+1Giur.Step 5: Increase theu dimension to u + 1 dimension, repeat Step1∼Step4, and get Gu+1r.Step 6: Calculate SE.(14)SampEnu,r=limN⟶∞−lnGu+1rGur.However, when SE is applied to actual data, N is always finite, thus giving an estimate of sample entropy as(15)SampEnu,r,N=−lnGu+1rGur.


## 4.3. Predictive Model (BiGRU-AM)

### 4.3.1. Bidirectional Gated Recurrent Unit

GRU is a gated recurrent neural network model that can solve problems such as long-term memory and gradient in backpropagation and is suitable for processing time series data. Compared with LSTM, the GRU model replaces the input gate and forget gate functions with only the update gate. The formula of the hidden layer of the model is as follows:(16)zt=sigmoidwzht−1,xt+bz,rt=sigmoidwrht−1,xt+br,ht˜=tanhwhrtht−1,xt+bh,ht=1−ztht−1+ztht˜,where zt is the update gate. rt is the reset gate; ht˜ is the newly generated information after processing the historical information through the update gate. wz, wr, and wh are the weight parameters of the GRU network; bz, br, bh are the thresholds of the GRU network.The GRU structure diagram is shown in Figure2.Figure 2
GRU structure diagram.The hidden state of the GRU model layer only considers the impact of historical information on the current input, and the impact of future information is equally important, so this paper uses BiGRU to solve this problem. BiGRU consists of two unidirectional and opposite GRUs. The model integrates the information before and after the current state, and the final output is the weighted summation of the output results of the forward GRU and the backward GRU and is linearly superimposed with the threshold. The main calculation formula of BiGRU is(17)hft=GRUht−1,xt,hbt=GRUht+1,xt,ht=αhft+βhbt+γ,where hft and hbt are the forward and backward GRU output vectors for the current state, and β are the weight parameters of the forward and backward GRU; γ is the linear threshold. The BiGRU network structure is shown in Figure 3.Figure 3
BiGRU network structure diagram.


### 4.3.2. Attention Mechanism

Although BiGRU improves the efficiency of information utilization to a certain extent, it still has the problem that local features cannot be accurately extracted when dealing with long-term sequence data. Therefore, when building the model, this paper adds an attention layer to the hidden layer to solve this problem. The attention layer further extracts the key information of the time series by calculating the weight of the output result in the BiGRU layer, thereby reducing the length of the output result and improving the robustness of the model:(18)vit=ustanhwhit+b,αit=evit∑i=1nevit,ct=∑i=1nαithit,vit is the unnormalized attention score. The higher the attention score, the better the match between the input vector and the target vector; us is the randomly initialized attention matrix; w is the weight coefficient; b is the bias coefficient; ct is the attention vector; hit is the value represented by the i-th output result at the t-th time.


## 4.3.1. Bidirectional Gated Recurrent Unit

GRU is a gated recurrent neural network model that can solve problems such as long-term memory and gradient in backpropagation and is suitable for processing time series data. Compared with LSTM, the GRU model replaces the input gate and forget gate functions with only the update gate. The formula of the hidden layer of the model is as follows:(16)zt=sigmoidwzht−1,xt+bz,rt=sigmoidwrht−1,xt+br,ht˜=tanhwhrtht−1,xt+bh,ht=1−ztht−1+ztht˜,where zt is the update gate. rt is the reset gate; ht˜ is the newly generated information after processing the historical information through the update gate. wz, wr, and wh are the weight parameters of the GRU network; bz, br, bh are the thresholds of the GRU network.The GRU structure diagram is shown in Figure2.Figure 2
GRU structure diagram.The hidden state of the GRU model layer only considers the impact of historical information on the current input, and the impact of future information is equally important, so this paper uses BiGRU to solve this problem. BiGRU consists of two unidirectional and opposite GRUs. The model integrates the information before and after the current state, and the final output is the weighted summation of the output results of the forward GRU and the backward GRU and is linearly superimposed with the threshold. The main calculation formula of BiGRU is(17)hft=GRUht−1,xt,hbt=GRUht+1,xt,ht=αhft+βhbt+γ,where hft and hbt are the forward and backward GRU output vectors for the current state, and β are the weight parameters of the forward and backward GRU; γ is the linear threshold. The BiGRU network structure is shown in Figure 3.Figure 3
BiGRU network structure diagram.


## 4.3.2. Attention Mechanism

Although BiGRU improves the efficiency of information utilization to a certain extent, it still has the problem that local features cannot be accurately extracted when dealing with long-term sequence data. Therefore, when building the model, this paper adds an attention layer to the hidden layer to solve this problem. The attention layer further extracts the key information of the time series by calculating the weight of the output result in the BiGRU layer, thereby reducing the length of the output result and improving the robustness of the model:(18)vit=ustanhwhit+b,αit=evit∑i=1nevit,ct=∑i=1nαithit,vit is the unnormalized attention score. The higher the attention score, the better the match between the input vector and the target vector; us is the randomly initialized attention matrix; w is the weight coefficient; b is the bias coefficient; ct is the attention vector; hit is the value represented by the i-th output result at the t-th time.


## 5. Data Analysis

### 5.1. Data Sources and Descriptions

This paper selects Automatic Fare Collection (AFC) data of Shanghai Metro Lines 1, 11, 2, 8, 9, 12, and 16 from April 1, 2015, to April 24, 2015, to test the established subway station passenger flow prediction model. The dataset contains timestamps, station names, and passenger arrival and departure records.


### 5.2. Station Clustering

Considering the large differences in passenger flow data in different time periods, the original data needs to be standardized first. Second, the optimal number of clusters needs to be determined. Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) are always used for trading off the complexity of the clustering model and the goodness of fitting the data:(19)AIC=−2lnL+2k,(20)BIC=−2lnL+lnNK.L is the log-likelihood function in equation (4). K is the number of clusters; N is the sample size.The second term in (19) and (20) is the penalty term. It can be seen that the penalty term of BIC considers the sample size, and the overall value is larger than that of AIC. Therefore, BIC can more effectively balance the model accuracy and complexity when the number of samples is too large, preventing the model from becoming too complex. Since 162 stations need to be clustered and due to the large sample size, only BIC is used to determine the optimal number of classifications. The smaller the BIC, the better the model classification effect. Set the number of clusters to a range from 1 to 20, and take the number of clusters corresponding to the minimum BIC value. From Figure 4, it can be determined that the optimal number of clusters is 5.Figure 4
Bayesian Information Criterion to determine the optimal number of clusters.Secondly, use WOA to find GMM parameters for the aggregated passenger flow data of 5 types of stations: the number of humpback whales is 80, and the maximum number of iterations is 500. Finally, the optimal initial weights of each Gaussian distribution are 0.09, 0.35, 0.31, 0.09, and 0.16.To illustrate the advantages and efficiency of the WOA, we still use AIC and BIC as evaluation metrics to compare the GMM and WOA-GMM. As shown in Table1, the model optimized by the WOA is smaller in both AIC and BIC indicators, which proves its excellent improvement effect.Table 1
WOA-GMM vs. GMM.
ModelAICBICWOA-GMM−16764.0204−788.7970GMM−16562.1785−586.9550Figures5 and 6, respectively, represent the 24-hour cumulative passenger flow changes in 5 clusters within 24 days. The typical characteristics of these clusters of stations are shown in Table 2.Figure 5
The inbound volume of 5 clusters.Figure 6
The outbound volume of 5 clusters.Table 2
Typical features of stations.
Number of stationsPassenger flow characteristicsRepresentative stationsSpatial characteristics of the stationCluster115The number of such stations is small, but the passenger flow is large. The morning and evening peaks of passenger flow are obvious, and the inbound volume in the evening peak is higher than in the morning peak. The outbound volume in the morning peak reaches an obvious peak, and the inbound volume is still large after 20 : 00.Qibao Station, Shanghai Railway Station, Zhongshan Park Station, People’s Square Station, etc.Most of these stations are subway transfer stations, railway stations, airport transportation hubs, or located in the city center, which belong to the city’s core business district. After the evening peak, the passenger flow is still relatively large, which further proves that the commercial development around these stations is relatively good, and more people work overtime at night and participate in entertainment activities.Cluster 257The number of such sites is large, and the traffic is low. Compared with other categories, the morning and evening peaks of cluster 2 are not very prominent, and the overall passenger flow is relatively low.Wild Animal Park Station, Donglu Road Station, Zhongchun Road Station, Lingang Avenue Station, etc.Such stations are mostly located in underdeveloped areas, where buildings are mostly nonresidential, nonoffice, or in scenic spots.Cluster 350The passenger flow of such stations has apparent peaks in the morning and evening. The number of inbound volumes in the morning peak is the highest among several types of stations, and the number of outbound stations in the evening peak is also larger.Sanlin Station, Shanghai Indoor Stadium Station, Century Park Station, Beixinjing Station, etc.The areas where such stations are located are primarily residential, and some are located in the outer suburbs, with many houses and relatively low house prices. Therefore, many citizens go to work during the morning rush hour and leave work during the evening rush hour.Cluster 414The overall passenger flow of such stations is low, and the inbound volume in the morning peak is the most obvious.Shanghainan Railway Station, Jiuting Station, Songjiang University Town Station, Shendu Highway Station, etc.This station is also surrounded by residential areas, similar to the third cluster, but the overall passenger flow is low.Cluster 526The morning peak departure volume and evening peak arrival volume of such stations are extremely obvious.Dongchang Road Station, China Art Museum Station, Hechuan Road Station, Shangcheng Road Station, etc.There are many R&D and technology buildings around such stations and rich industrial parks.


### 5.3. Decomposition and Noise Reduction of Station Passenger Flow Sequence

Next, this paper takes the inbound passenger flow of the first cluster of the station as an example to illustrate the process of data decomposition and noise reduction.Using CEEMD, the empirical mode decomposition process of adding complementary white noise is carried out step by step for the daily hourly inbound and outbound passenger flow data of five clusters of subway stations. Figure7 represents the decomposition of the corresponding data of the first cluster of subway station inbound passenger flow into 10 IMF components, the last of which represents the residual.Figure 7
The first cluster of station inbound passenger flow IMF.As shown in Figure7, many of the components have similar complexities. IMF1∼IMF5 have high complexity, IM6∼IMF7 have lower complexity and a particular trend, and IMF8∼R have obvious trends. However, conclusions based on observation alone are not accurate, so this paper uses sample entropy (SE) to determine the complexity of these components.Figure8 represents the changing trends of inbound passenger flow of the first cluster of the station, respectively. Next, to reduce the training time and improve the model prediction efficiency, this paper combines the IMF components of the inbound passenger flow of various types of stations into the combined components in Table 3 according to the IMF components with similar complexity as shown in Figure 8.Figure 8
SE variation trend of inbound IMF for the first cluster of the station.Table 3
Grouping of IMF components of inbound passenger flow at the first cluster of stations.
ComponentsIMFFrequency and randomnessComponent 1IMF1, IMF2High-frequency componentComponent 2IMF3, IMF4, IMF5High-frequency componentComponent 3IMF6, IMF7Low-frequency componentComponent 4IMF8, IMF9, RTrend componentThe trend of the combined components over time is shown in Figure9. In this paper, the recombination components Comp1 and Comp2 with high frequency and randomness are regarded as high-frequency components; the recombination components Comp3 with low frequency and certain randomness are regarded as low-frequency components; the recombination components Comp4 with the lowest frequency are regarded as trend components. The assignment results are shown in Table 3.Figure 9
The trend of the combined components over time through CEEMD-SE.


### 5.4. Passenger Flow Forecasting Process

This paper optimizes BiGRU based on the attention mechanism, implements it inPython, and uses the TensorFlow framework to build the model. The experimental hardware device is configured as an 11th Gen Intel(R) Core(TM) i7-1165G7 with 16G memory. According to the prediction module in the model process, the model parameters are first initialized. The number of nodes in the input layer is set to three according to the number of auxiliary prediction variables (date, hour, subway station type), and the number of hidden layers of the BiGRU-AM network is selected to be two. Adam is used to continuously adjust the hyperparameters of the BiGRU-AM model in the iterative process until the error value is small.In the process of inbound passenger flow prediction for the first station cluster, the Root Mean Square Error (RMSE) is used to evaluate the prediction accuracy, and its expression is shown in equation (21). The smaller the RMSE, the higher the prediction accuracy. Figure 10 shows the effect of setting the number of neurons in the two layers in BiGRU on the prediction accuracy of each combined component. It can be seen that when the number of neurons in the two layers is 13 and 11, respectively, the overall predicted RMSE reaches the minimum. The same method is used for the number of neurons in the BiGRU layer of other clusters of inbound and outbound passenger flow prediction models. Next, choose the optimal learning rate with the optimal number of neurons already determined. To simplify the model, a unified learning rate is determined separately in the prediction process of inbound passenger flow and outbound passenger flow. Taking the inbound passenger flow as an example, Figure 11 shows the variation of the inbound passenger flow prediction accuracy with the learning rate for each station cluster. It can be seen that the RMSE is the smallest when the learning rate is 0.006. Similarly, it can be found that the optimal learning rate for inbound passenger flow prediction should also be set to 0.006.Figure 10
Hidden units determination.Figure 11
Learning rate determination.This paper sets the step size to 20 and the batch size to 110. According to experience, the number of neurons in the attention layer is set to 64. In this process, the optimal number of iterations for inbound and outbound passenger flow and the number of neurons in the two GRU layers in BiGRU are shown in Tables4 and 5, respectively. Among them, the number of iterations is the average number of iterative predictions for each IMF combination, and the ratio of the training set and test set is 9 : 1.Table 4
Inbound passenger flow prediction hyperparameter settings.
Cluster 1Cluster 2Cluster 3Cluster 4Cluster 5Average number of iterations499421437322323BiGRU hidden units 11316221618BiGRU hidden units 21114211417Table 5
Outbound passenger flow prediction hyperparameter settings.
Cluster 1Cluster 2Cluster 3Cluster 4Cluster 5Average number of iterations278251486522503BiGRU hidden units 11316222018BiGRU hidden units 21114211917In order to verify the validity of the BiGRU-AM model in the prediction of subway passenger flow, experiments are carried out based on the above data. The prediction results of the combined components prediction results are shown in Figure12. Figure 12 shows that the model’s prediction results for Comp1∼Comp3 are close to the true value, while the prediction for Comp4 is slightly different from the true value. It further illustrates the advantages of decomposing the passenger flow series into more trending components. In later chapters, we will further evaluate the model with more accurate metrics.Figure 12
Prediction results of inbound traffic by the combined components of the first cluster of the station.


### 5.5. Predictive Model Comparative Analysis

In order to further illustrate the effectiveness of the overall model, after classifying the stations, BPNN, LSTM, and BiGRU are first selected for comparison. Secondly, the BiGRU model and the BiGRU-AM model are selected to compare with the CEEMD-BiGRU-AM model finally constructed in this paper to prove the accuracy of the basic model and the effect of model improvement.There are many error evaluation indicators in subway passenger flow forecasting. This paper adopts the following four commonly used error evaluation indicators: MAPE, RMSE, Mean Absolute Error (MAE), and Coefficient of Determination (R2). The expression is as follows:(21)MAPE=1n∑i=1ny^i−yiyi×100%,RMSE=1n∑i=1ny^i−yi2,MAE=1n∑i=1ny^i−yi,R2=1−∑i=1ny^i−yi2∑i=1ny¯−yi2,where n is the number of test set data; yi and y^i are the actual and predicted values at a time, respectively; y¯ is the sample mean.Among them, MAPE represents the difference between the actual value and the predicted value as a percentage of the predicted value itself, RMSE represents the average difference between the predicted value and the actual value, and MAE represents the average absolute difference between the predicted value and the actual value. The smaller these three indicators are, the better the model prediction effect is.R2 represents the proportion of the total variation of predicted passenger flow that can be explained by independent variables such as station location, entry, and exit time. The larger the R2, the better the prediction effect of the model.Figure13 compares the forecasting effects of different models for the first cluster of subway station inbound passenger flow in a short period of time. It can be seen from the figure that, compared with other models, the CEEMD-BIGRE-AM model has the highest fitting degree to the passenger flow sequence and is closer to the actual value. In comparison, other models are more different from the actual value. The basic hyperparameters of the comparison models in this paper are the same as those of the main model, and the settings mentioned in Section 5.4 are the same.Figure 13
Prediction effect of different models on the inbound volume of the first cluster of the station.Next, we will use the four model evaluation indicators MAPE, RMSE, |R2, and MAE to comprehensively analyze the prediction accuracy of different models. The prediction errors of different models for inbound and outbound passenger flow are shown in Tables 6 and 7.Table 6
Prediction results of inbound passenger flow by different models.
ModelMAPERMSEMAER2Cluster1BP0.5802519.7456329.70530.9685LSTM0.9802486.2718294.51460.9725BiGRU0.5616347.8101224.68940.9859BiGRU-AM0.2508278.0797181.26150.9910CEEMD-SE-BiGRU-AM0.1846106.761752.92220.9945Cluster 2BP1.489490.005664.19100.9367LSTM1.370277.620058.88330.9529BiGRU1.360979.483258.31510.9507BiGRU-AM0.867376.734254.21090.9540CEEMD-SE-BiGRU-AM0.530351.974535.69990.9789Cluster 3BP0.5435123.659789.28090.9763LSTM0.371974.509254.43060.9914BiGRU0.336974.496953.70100.9914BiGRU-AM0.444572.353550.88860.9919CEEMD-SE-BiGRU-AM0.138822.130814.69560.9992Cluster 4BP1.377207.373141.9360.9840LSTM1.0932151.2548100.17050.9915BiGRU0.7577133.33796.89340.9934BiGRU-AM0.8194116.705784.0850.9949CEEMD-SE-BiGRU-AM0.3983113.67157.77370.9952Cluster 5BP1.213245.1781174.34730.9744LSTM0.9732161.1152110.44290.9890BiGRU1.2524152.2723103.41690.9901BiGRU-AM1.7635147.83298.47950.9907CEEMD-SE-BiGRU-AM0.9436138.670889.93030.9918Table 7
Prediction results of inbound passenger flow by different models.
ModelMAPERMSEMAER2Cluster1BP0.5411529.8081354.00680.9624LSTM0.8076580.7955338.06790.9548BiGRU0.2570393.9422252.44330.9841BiGRU-AM0.6007442.1145207.97810.9843CEEMD-SE-BiGRU-AM0.1971330.0983217.47310.9854Cluster 2BP1.256984.762257.38320.9258LSTM1.201787.934663.81130.9202BiGRU0.742653.909740.05800.9700BiGRU-AM0.511749.360236.03840.9749CEEMD-SE-BiGRU-AM0.411448.001735.43400.9762Cluster 3BP0.4730143.0040103.49390.9525LSTM0.442592.735769.72110.9800BiGRU0.298285.015961.71470.9832BiGRU-AM0.207575.810754.93530.9867CEEMD-SE-BiGRU-AM0.063120.171014.30870.9991Cluster 4BP0.8992145.4933104.82150.9880LSTM1.2938132.956699.95090.9890BiGRU0.8559130.311387.76980.9904BiGRU-AM0.6187118.615484.53650.9922CEEMD-SE-BiGRU-AM0.256368.543841.44840.9974Cluster 5BP0.8703220.3739151.09480.9844LSTM0.6068140.215989.99980.9936BiGRU0.8587132.488388.12450.9943BiGRU-AM0.5882136.900488.40330.9940CEEMD-SE-BiGRU-AM0.164453.192833.09570.9991It can be seen from the comparison that the CEEMD-BiGRU-AM model outperforms the other five models on the four error evaluation indicators of MAPE, RMSE, MAE, andR2. For predicting each type of passenger flow, the model established in this paper has reached the minimum in MAPE, RMSE, and MAE indicators and the maximum in R2. Regardless of whether the overall passenger flow of the station is large, the morning and evening peak passenger flow is large, or the passenger flow is small.(1)
BP neural network easily falls into local optimum, resulting in problems such as gradient disappearance. Although LSTM solves this problem, its hidden layer only considers the impact of historical information on the current input, and GRU as a simplified version also has this problem. BiGRU considers the equally important future information by combining two GRUs in opposite directions. Tables6 and 7 show that the standard BiGRU model outperforms the BP and LSTM models in most prediction situations and metrics. It verifies the superiority of using the BiGRU model as the base model.(2)
The model can extract local features more accurately after using the attention mechanism to improve BiGRU. The BiGRU-AM model has been far superior to the BiGRU, LSTM, and BP neural network models for each station’s inbound and outbound passenger flow prediction cluster. For the prediction of the fourth type of station departures, the BiGRU-AM model reduced 31.19%, 52.18%, and 27.71% in MAPE indicators compared with BP neural network, LSTM, and BiGRU, respectively.(3)
The data becomes more stable after using the CEEMD-SE model to decompose and reconstruct the passenger flow time series data, and the prediction accuracy is further improved. In terms of the RMSE index, the station’s entry and exit prediction accuracy are reduced by more than 40% on average compared to BP and LSTM and increased by more than 0.8% on average.


### 5.6. Statistical Testing of Predictive Models

Considering that the above four evaluation indicators may have generalization and randomness, this paper further uses the statistical test method to verify the model’s overall prediction effect. In order to verify the overall statistical significance of the combined model proposed in this paper, the Friedman test, Nemenyi test [37], and paired t-test [38] were used to verify the differences between different models.The Friedman test can determine whether different algorithms perform significantly on multiple datasets without the need for normality assumptions. In this section, this paper uses the inbound and outbound traffic of each type of passenger flow as a whole to consider 10 sets of data sets and studies the differences between different models on the four measurement indicators MAPE, RMSE, MAE, andR2, under the condition of a significance level α=0.05. The null hypothesis is that there is no significant difference in the performance of the five algorithms, and the alternative hypothesis is that the performance of the five models is significantly different. It is necessary to sort each algorithm and its corresponding evaluation index corresponding to each data, which set it as 1, 2, 3, 4, and 5 from the 1st to the 5th.Due to the small sample size, it cannot be calculated directly using the chi-square statistic, but the F statistic is constructed from it:(22)τχ2=k−1k⋅12Nk2−1∑i=1kri−k+122τF=N−1τχ2Nk−1−τχ2,k is the number of algorithms to be compared, which is 5 here; N is the number of datasets, 10 here. ri is the average ordinal value of the i-th algorithm. τχ2 is the chi-square statistic. τF is the F statistic and obeys the F distribution with k-1 and (k-1) (N-1) degrees of freedom.The F-test here has a critical value ofF0.054,36=2.634 and a significance level of 0.05. Through calculation, the p-values (probability that the test statistic is greater than the critical value) corresponding to different metrics are shown in Table 8.Table 8
p-value of the Friedman test.
MAPEMAERMSER2p-value1.745e − 051.256e − 072.394e − 071.977e − 07It can be found that, for all indicators, thep-values of all indicators are much less than 0.01. Therefore, the null hypothesis is rejected, and the prediction effects among the five models are considered to be significantly different.On this basis, it is necessary to continue to compare whether there are significant differences between the two models. In this regard, the Shapiro-Wilk test [39] is firstly performed on the evaluation index corresponding to each algorithm to judge whether it conforms to the normality assumption. The null hypothesis for this test is that the data follow a normal distribution. So if the p-value is less than 0.05, the null hypothesis is rejected, indicating that the data are not from a normal distribution. Its test statistic is(23)W=∑i=1naixi2∑i=1nxi−x¯2,xi is the i-th smallest sample value in the data; x¯ is the average value of the sample; xi is the i-th sample value; ai is a constant that meets certain conditions, which will not be repeated here.Next, using the Shapiro-Wilk test, it can be seen that, among the four indicators, only the five models corresponding to the RMSE after logarithmic transformation do not reject the assumption of normal distribution. Thep-values of the tests are shown in Table 9.Table 9
p-value of the Shapiro-Wilk test.
BPLSTMBiGRUBiGRU-AMCEEMD-SE-BiGRU-AMMAPE0.16650.55670.40620.02190.0478MAE0.02430.00080.00730.02930.0007RMSE0.00900.00060.00530.00500.0037R20.23900.02240.00600.00200.0370ln(RMSE)0.32500.07580.28610.50420.7861As can be seen in Table9, most of the models have p-values less than 0.05 for MAPE, MAE, and metrics. Therefore, we performed a further test using Nemenyi’s test, which does not require normality to be assumed. The null hypothesis is that there is no significant difference between the two models, and the alternative hypothesis is that there is a significant difference between the two models. The statistic is(24)CD=qαkk+16N,qα is the critical value of α of the Tukey distribution, where α=0.05 is taken and the value is 2.728. k is the number of algorithms to be compared, which is 5. N is the number of datasets, 10 here.The critical value is 1.9290 by calculation, and the average sequence value of the five models is shown in Table10. The difference in mean ordinal values between the two pairs is shown in Figure 14.Table 10
The average sequence value of five models.
BPLSTMBiGRUBiGRU-AMCEEMD-SE-BiGRU-AMMAPE4.43.92.92.81MAE4.94.12.92.01.1R24.84.053.052.11Figure 14
Nemenyi test of MAPE, MAE, andR2 (CSBA replaces the CEEMD-SE-BiGRU-AM model in the figure), (a) MAPE, (b) MAE, (c) R2.
(a)(b)(c)In the case of a small amount of data, the Nemenyi test rejected the null hypothesis that the two models are inconsistent. When the average ordinal difference between the two models is only greater than the critical value of 1.929, the null hypothesis is rejected, and there is a significant difference between the two contrasting models. It can be seen from Figure14 that the combined model established in this paper is significantly different from other models, and the differences between the BiGRU-AM model and BP and LSTM are also significant, while other models cannot directly reject the null hypothesis.Then, for the ln(RMSE) that meets the normality assumption, a pairedt-test with the same null hypothesis is performed between the two models, and its statistic is(25)t=d¯s/N,d¯ is the mean of the sample differences; N is the sample size of 10; s is the standard deviation of the sample differences.The statistic (25) follows a t-distribution with 9 degrees of freedom. It is obtained that the significance level is α=0.05 the p-value of the paired t-test between the models, as shown in Table 11.Table 11
p-value of the paired t-test.
BPLSTMBiGRUBiGRU-AMCEEMD-SE-BiGRU-AMBP10.0092<0.0001<0.00010.0004LSTM0.009210.02560.01040.0007BiGRU<0.00010.025610.07360.0036BiGRU-AM<0.00010.01040.073610.0050CEEMD-SE-BiGRU-AM0.00040.00070.00360.00501At a significance levelα=0.05, the above test can prove the following:(1)
The model after comprehensively using CEEMD and AM to improve BiGRU is significantly different from other comparable models.(2)
BiGRU already has a significant difference compared with BP and LSTM and has a significant difference with BiGRU-AM under the significance levelα=0.05.(3)
Combining the results of Tables9 and 10, it can be found that the model established in this paper is significantly better than other comparative models. AM and CEEMD have a noticeable improvement effect on the BiGRU model.


## 5.1. Data Sources and Descriptions

This paper selects Automatic Fare Collection (AFC) data of Shanghai Metro Lines 1, 11, 2, 8, 9, 12, and 16 from April 1, 2015, to April 24, 2015, to test the established subway station passenger flow prediction model. The dataset contains timestamps, station names, and passenger arrival and departure records.


## 5.2. Station Clustering

Considering the large differences in passenger flow data in different time periods, the original data needs to be standardized first. Second, the optimal number of clusters needs to be determined. Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) are always used for trading off the complexity of the clustering model and the goodness of fitting the data:(19)AIC=−2lnL+2k,(20)BIC=−2lnL+lnNK.L is the log-likelihood function in equation (4). K is the number of clusters; N is the sample size.The second term in (19) and (20) is the penalty term. It can be seen that the penalty term of BIC considers the sample size, and the overall value is larger than that of AIC. Therefore, BIC can more effectively balance the model accuracy and complexity when the number of samples is too large, preventing the model from becoming too complex. Since 162 stations need to be clustered and due to the large sample size, only BIC is used to determine the optimal number of classifications. The smaller the BIC, the better the model classification effect. Set the number of clusters to a range from 1 to 20, and take the number of clusters corresponding to the minimum BIC value. From Figure 4, it can be determined that the optimal number of clusters is 5.Figure 4
Bayesian Information Criterion to determine the optimal number of clusters.Secondly, use WOA to find GMM parameters for the aggregated passenger flow data of 5 types of stations: the number of humpback whales is 80, and the maximum number of iterations is 500. Finally, the optimal initial weights of each Gaussian distribution are 0.09, 0.35, 0.31, 0.09, and 0.16.To illustrate the advantages and efficiency of the WOA, we still use AIC and BIC as evaluation metrics to compare the GMM and WOA-GMM. As shown in Table1, the model optimized by the WOA is smaller in both AIC and BIC indicators, which proves its excellent improvement effect.Table 1
WOA-GMM vs. GMM.
ModelAICBICWOA-GMM−16764.0204−788.7970GMM−16562.1785−586.9550Figures5 and 6, respectively, represent the 24-hour cumulative passenger flow changes in 5 clusters within 24 days. The typical characteristics of these clusters of stations are shown in Table 2.Figure 5
The inbound volume of 5 clusters.Figure 6
The outbound volume of 5 clusters.Table 2
Typical features of stations.
Number of stationsPassenger flow characteristicsRepresentative stationsSpatial characteristics of the stationCluster115The number of such stations is small, but the passenger flow is large. The morning and evening peaks of passenger flow are obvious, and the inbound volume in the evening peak is higher than in the morning peak. The outbound volume in the morning peak reaches an obvious peak, and the inbound volume is still large after 20 : 00.Qibao Station, Shanghai Railway Station, Zhongshan Park Station, People’s Square Station, etc.Most of these stations are subway transfer stations, railway stations, airport transportation hubs, or located in the city center, which belong to the city’s core business district. After the evening peak, the passenger flow is still relatively large, which further proves that the commercial development around these stations is relatively good, and more people work overtime at night and participate in entertainment activities.Cluster 257The number of such sites is large, and the traffic is low. Compared with other categories, the morning and evening peaks of cluster 2 are not very prominent, and the overall passenger flow is relatively low.Wild Animal Park Station, Donglu Road Station, Zhongchun Road Station, Lingang Avenue Station, etc.Such stations are mostly located in underdeveloped areas, where buildings are mostly nonresidential, nonoffice, or in scenic spots.Cluster 350The passenger flow of such stations has apparent peaks in the morning and evening. The number of inbound volumes in the morning peak is the highest among several types of stations, and the number of outbound stations in the evening peak is also larger.Sanlin Station, Shanghai Indoor Stadium Station, Century Park Station, Beixinjing Station, etc.The areas where such stations are located are primarily residential, and some are located in the outer suburbs, with many houses and relatively low house prices. Therefore, many citizens go to work during the morning rush hour and leave work during the evening rush hour.Cluster 414The overall passenger flow of such stations is low, and the inbound volume in the morning peak is the most obvious.Shanghainan Railway Station, Jiuting Station, Songjiang University Town Station, Shendu Highway Station, etc.This station is also surrounded by residential areas, similar to the third cluster, but the overall passenger flow is low.Cluster 526The morning peak departure volume and evening peak arrival volume of such stations are extremely obvious.Dongchang Road Station, China Art Museum Station, Hechuan Road Station, Shangcheng Road Station, etc.There are many R&D and technology buildings around such stations and rich industrial parks.


## 5.3. Decomposition and Noise Reduction of Station Passenger Flow Sequence

Next, this paper takes the inbound passenger flow of the first cluster of the station as an example to illustrate the process of data decomposition and noise reduction.Using CEEMD, the empirical mode decomposition process of adding complementary white noise is carried out step by step for the daily hourly inbound and outbound passenger flow data of five clusters of subway stations. Figure7 represents the decomposition of the corresponding data of the first cluster of subway station inbound passenger flow into 10 IMF components, the last of which represents the residual.Figure 7
The first cluster of station inbound passenger flow IMF.As shown in Figure7, many of the components have similar complexities. IMF1∼IMF5 have high complexity, IM6∼IMF7 have lower complexity and a particular trend, and IMF8∼R have obvious trends. However, conclusions based on observation alone are not accurate, so this paper uses sample entropy (SE) to determine the complexity of these components.Figure8 represents the changing trends of inbound passenger flow of the first cluster of the station, respectively. Next, to reduce the training time and improve the model prediction efficiency, this paper combines the IMF components of the inbound passenger flow of various types of stations into the combined components in Table 3 according to the IMF components with similar complexity as shown in Figure 8.Figure 8
SE variation trend of inbound IMF for the first cluster of the station.Table 3
Grouping of IMF components of inbound passenger flow at the first cluster of stations.
ComponentsIMFFrequency and randomnessComponent 1IMF1, IMF2High-frequency componentComponent 2IMF3, IMF4, IMF5High-frequency componentComponent 3IMF6, IMF7Low-frequency componentComponent 4IMF8, IMF9, RTrend componentThe trend of the combined components over time is shown in Figure9. In this paper, the recombination components Comp1 and Comp2 with high frequency and randomness are regarded as high-frequency components; the recombination components Comp3 with low frequency and certain randomness are regarded as low-frequency components; the recombination components Comp4 with the lowest frequency are regarded as trend components. The assignment results are shown in Table 3.Figure 9
The trend of the combined components over time through CEEMD-SE.


## 5.4. Passenger Flow Forecasting Process

This paper optimizes BiGRU based on the attention mechanism, implements it inPython, and uses the TensorFlow framework to build the model. The experimental hardware device is configured as an 11th Gen Intel(R) Core(TM) i7-1165G7 with 16G memory. According to the prediction module in the model process, the model parameters are first initialized. The number of nodes in the input layer is set to three according to the number of auxiliary prediction variables (date, hour, subway station type), and the number of hidden layers of the BiGRU-AM network is selected to be two. Adam is used to continuously adjust the hyperparameters of the BiGRU-AM model in the iterative process until the error value is small.In the process of inbound passenger flow prediction for the first station cluster, the Root Mean Square Error (RMSE) is used to evaluate the prediction accuracy, and its expression is shown in equation (21). The smaller the RMSE, the higher the prediction accuracy. Figure 10 shows the effect of setting the number of neurons in the two layers in BiGRU on the prediction accuracy of each combined component. It can be seen that when the number of neurons in the two layers is 13 and 11, respectively, the overall predicted RMSE reaches the minimum. The same method is used for the number of neurons in the BiGRU layer of other clusters of inbound and outbound passenger flow prediction models. Next, choose the optimal learning rate with the optimal number of neurons already determined. To simplify the model, a unified learning rate is determined separately in the prediction process of inbound passenger flow and outbound passenger flow. Taking the inbound passenger flow as an example, Figure 11 shows the variation of the inbound passenger flow prediction accuracy with the learning rate for each station cluster. It can be seen that the RMSE is the smallest when the learning rate is 0.006. Similarly, it can be found that the optimal learning rate for inbound passenger flow prediction should also be set to 0.006.Figure 10
Hidden units determination.Figure 11
Learning rate determination.This paper sets the step size to 20 and the batch size to 110. According to experience, the number of neurons in the attention layer is set to 64. In this process, the optimal number of iterations for inbound and outbound passenger flow and the number of neurons in the two GRU layers in BiGRU are shown in Tables4 and 5, respectively. Among them, the number of iterations is the average number of iterative predictions for each IMF combination, and the ratio of the training set and test set is 9 : 1.Table 4
Inbound passenger flow prediction hyperparameter settings.
Cluster 1Cluster 2Cluster 3Cluster 4Cluster 5Average number of iterations499421437322323BiGRU hidden units 11316221618BiGRU hidden units 21114211417Table 5
Outbound passenger flow prediction hyperparameter settings.
Cluster 1Cluster 2Cluster 3Cluster 4Cluster 5Average number of iterations278251486522503BiGRU hidden units 11316222018BiGRU hidden units 21114211917In order to verify the validity of the BiGRU-AM model in the prediction of subway passenger flow, experiments are carried out based on the above data. The prediction results of the combined components prediction results are shown in Figure12. Figure 12 shows that the model’s prediction results for Comp1∼Comp3 are close to the true value, while the prediction for Comp4 is slightly different from the true value. It further illustrates the advantages of decomposing the passenger flow series into more trending components. In later chapters, we will further evaluate the model with more accurate metrics.Figure 12
Prediction results of inbound traffic by the combined components of the first cluster of the station.


## 5.5. Predictive Model Comparative Analysis

In order to further illustrate the effectiveness of the overall model, after classifying the stations, BPNN, LSTM, and BiGRU are first selected for comparison. Secondly, the BiGRU model and the BiGRU-AM model are selected to compare with the CEEMD-BiGRU-AM model finally constructed in this paper to prove the accuracy of the basic model and the effect of model improvement.There are many error evaluation indicators in subway passenger flow forecasting. This paper adopts the following four commonly used error evaluation indicators: MAPE, RMSE, Mean Absolute Error (MAE), and Coefficient of Determination (R2). The expression is as follows:(21)MAPE=1n∑i=1ny^i−yiyi×100%,RMSE=1n∑i=1ny^i−yi2,MAE=1n∑i=1ny^i−yi,R2=1−∑i=1ny^i−yi2∑i=1ny¯−yi2,where n is the number of test set data; yi and y^i are the actual and predicted values at a time, respectively; y¯ is the sample mean.Among them, MAPE represents the difference between the actual value and the predicted value as a percentage of the predicted value itself, RMSE represents the average difference between the predicted value and the actual value, and MAE represents the average absolute difference between the predicted value and the actual value. The smaller these three indicators are, the better the model prediction effect is.R2 represents the proportion of the total variation of predicted passenger flow that can be explained by independent variables such as station location, entry, and exit time. The larger the R2, the better the prediction effect of the model.Figure13 compares the forecasting effects of different models for the first cluster of subway station inbound passenger flow in a short period of time. It can be seen from the figure that, compared with other models, the CEEMD-BIGRE-AM model has the highest fitting degree to the passenger flow sequence and is closer to the actual value. In comparison, other models are more different from the actual value. The basic hyperparameters of the comparison models in this paper are the same as those of the main model, and the settings mentioned in Section 5.4 are the same.Figure 13
Prediction effect of different models on the inbound volume of the first cluster of the station.Next, we will use the four model evaluation indicators MAPE, RMSE, |R2, and MAE to comprehensively analyze the prediction accuracy of different models. The prediction errors of different models for inbound and outbound passenger flow are shown in Tables 6 and 7.Table 6
Prediction results of inbound passenger flow by different models.
ModelMAPERMSEMAER2Cluster1BP0.5802519.7456329.70530.9685LSTM0.9802486.2718294.51460.9725BiGRU0.5616347.8101224.68940.9859BiGRU-AM0.2508278.0797181.26150.9910CEEMD-SE-BiGRU-AM0.1846106.761752.92220.9945Cluster 2BP1.489490.005664.19100.9367LSTM1.370277.620058.88330.9529BiGRU1.360979.483258.31510.9507BiGRU-AM0.867376.734254.21090.9540CEEMD-SE-BiGRU-AM0.530351.974535.69990.9789Cluster 3BP0.5435123.659789.28090.9763LSTM0.371974.509254.43060.9914BiGRU0.336974.496953.70100.9914BiGRU-AM0.444572.353550.88860.9919CEEMD-SE-BiGRU-AM0.138822.130814.69560.9992Cluster 4BP1.377207.373141.9360.9840LSTM1.0932151.2548100.17050.9915BiGRU0.7577133.33796.89340.9934BiGRU-AM0.8194116.705784.0850.9949CEEMD-SE-BiGRU-AM0.3983113.67157.77370.9952Cluster 5BP1.213245.1781174.34730.9744LSTM0.9732161.1152110.44290.9890BiGRU1.2524152.2723103.41690.9901BiGRU-AM1.7635147.83298.47950.9907CEEMD-SE-BiGRU-AM0.9436138.670889.93030.9918Table 7
Prediction results of inbound passenger flow by different models.
ModelMAPERMSEMAER2Cluster1BP0.5411529.8081354.00680.9624LSTM0.8076580.7955338.06790.9548BiGRU0.2570393.9422252.44330.9841BiGRU-AM0.6007442.1145207.97810.9843CEEMD-SE-BiGRU-AM0.1971330.0983217.47310.9854Cluster 2BP1.256984.762257.38320.9258LSTM1.201787.934663.81130.9202BiGRU0.742653.909740.05800.9700BiGRU-AM0.511749.360236.03840.9749CEEMD-SE-BiGRU-AM0.411448.001735.43400.9762Cluster 3BP0.4730143.0040103.49390.9525LSTM0.442592.735769.72110.9800BiGRU0.298285.015961.71470.9832BiGRU-AM0.207575.810754.93530.9867CEEMD-SE-BiGRU-AM0.063120.171014.30870.9991Cluster 4BP0.8992145.4933104.82150.9880LSTM1.2938132.956699.95090.9890BiGRU0.8559130.311387.76980.9904BiGRU-AM0.6187118.615484.53650.9922CEEMD-SE-BiGRU-AM0.256368.543841.44840.9974Cluster 5BP0.8703220.3739151.09480.9844LSTM0.6068140.215989.99980.9936BiGRU0.8587132.488388.12450.9943BiGRU-AM0.5882136.900488.40330.9940CEEMD-SE-BiGRU-AM0.164453.192833.09570.9991It can be seen from the comparison that the CEEMD-BiGRU-AM model outperforms the other five models on the four error evaluation indicators of MAPE, RMSE, MAE, andR2. For predicting each type of passenger flow, the model established in this paper has reached the minimum in MAPE, RMSE, and MAE indicators and the maximum in R2. Regardless of whether the overall passenger flow of the station is large, the morning and evening peak passenger flow is large, or the passenger flow is small.(1)
BP neural network easily falls into local optimum, resulting in problems such as gradient disappearance. Although LSTM solves this problem, its hidden layer only considers the impact of historical information on the current input, and GRU as a simplified version also has this problem. BiGRU considers the equally important future information by combining two GRUs in opposite directions. Tables6 and 7 show that the standard BiGRU model outperforms the BP and LSTM models in most prediction situations and metrics. It verifies the superiority of using the BiGRU model as the base model.(2)
The model can extract local features more accurately after using the attention mechanism to improve BiGRU. The BiGRU-AM model has been far superior to the BiGRU, LSTM, and BP neural network models for each station’s inbound and outbound passenger flow prediction cluster. For the prediction of the fourth type of station departures, the BiGRU-AM model reduced 31.19%, 52.18%, and 27.71% in MAPE indicators compared with BP neural network, LSTM, and BiGRU, respectively.(3)
The data becomes more stable after using the CEEMD-SE model to decompose and reconstruct the passenger flow time series data, and the prediction accuracy is further improved. In terms of the RMSE index, the station’s entry and exit prediction accuracy are reduced by more than 40% on average compared to BP and LSTM and increased by more than 0.8% on average.


## 5.6. Statistical Testing of Predictive Models

Considering that the above four evaluation indicators may have generalization and randomness, this paper further uses the statistical test method to verify the model’s overall prediction effect. In order to verify the overall statistical significance of the combined model proposed in this paper, the Friedman test, Nemenyi test [37], and paired t-test [38] were used to verify the differences between different models.The Friedman test can determine whether different algorithms perform significantly on multiple datasets without the need for normality assumptions. In this section, this paper uses the inbound and outbound traffic of each type of passenger flow as a whole to consider 10 sets of data sets and studies the differences between different models on the four measurement indicators MAPE, RMSE, MAE, andR2, under the condition of a significance level α=0.05. The null hypothesis is that there is no significant difference in the performance of the five algorithms, and the alternative hypothesis is that the performance of the five models is significantly different. It is necessary to sort each algorithm and its corresponding evaluation index corresponding to each data, which set it as 1, 2, 3, 4, and 5 from the 1st to the 5th.Due to the small sample size, it cannot be calculated directly using the chi-square statistic, but the F statistic is constructed from it:(22)τχ2=k−1k⋅12Nk2−1∑i=1kri−k+122τF=N−1τχ2Nk−1−τχ2,k is the number of algorithms to be compared, which is 5 here; N is the number of datasets, 10 here. ri is the average ordinal value of the i-th algorithm. τχ2 is the chi-square statistic. τF is the F statistic and obeys the F distribution with k-1 and (k-1) (N-1) degrees of freedom.The F-test here has a critical value ofF0.054,36=2.634 and a significance level of 0.05. Through calculation, the p-values (probability that the test statistic is greater than the critical value) corresponding to different metrics are shown in Table 8.Table 8
p-value of the Friedman test.
MAPEMAERMSER2p-value1.745e − 051.256e − 072.394e − 071.977e − 07It can be found that, for all indicators, thep-values of all indicators are much less than 0.01. Therefore, the null hypothesis is rejected, and the prediction effects among the five models are considered to be significantly different.On this basis, it is necessary to continue to compare whether there are significant differences between the two models. In this regard, the Shapiro-Wilk test [39] is firstly performed on the evaluation index corresponding to each algorithm to judge whether it conforms to the normality assumption. The null hypothesis for this test is that the data follow a normal distribution. So if the p-value is less than 0.05, the null hypothesis is rejected, indicating that the data are not from a normal distribution. Its test statistic is(23)W=∑i=1naixi2∑i=1nxi−x¯2,xi is the i-th smallest sample value in the data; x¯ is the average value of the sample; xi is the i-th sample value; ai is a constant that meets certain conditions, which will not be repeated here.Next, using the Shapiro-Wilk test, it can be seen that, among the four indicators, only the five models corresponding to the RMSE after logarithmic transformation do not reject the assumption of normal distribution. Thep-values of the tests are shown in Table 9.Table 9
p-value of the Shapiro-Wilk test.
BPLSTMBiGRUBiGRU-AMCEEMD-SE-BiGRU-AMMAPE0.16650.55670.40620.02190.0478MAE0.02430.00080.00730.02930.0007RMSE0.00900.00060.00530.00500.0037R20.23900.02240.00600.00200.0370ln(RMSE)0.32500.07580.28610.50420.7861As can be seen in Table9, most of the models have p-values less than 0.05 for MAPE, MAE, and metrics. Therefore, we performed a further test using Nemenyi’s test, which does not require normality to be assumed. The null hypothesis is that there is no significant difference between the two models, and the alternative hypothesis is that there is a significant difference between the two models. The statistic is(24)CD=qαkk+16N,qα is the critical value of α of the Tukey distribution, where α=0.05 is taken and the value is 2.728. k is the number of algorithms to be compared, which is 5. N is the number of datasets, 10 here.The critical value is 1.9290 by calculation, and the average sequence value of the five models is shown in Table10. The difference in mean ordinal values between the two pairs is shown in Figure 14.Table 10
The average sequence value of five models.
BPLSTMBiGRUBiGRU-AMCEEMD-SE-BiGRU-AMMAPE4.43.92.92.81MAE4.94.12.92.01.1R24.84.053.052.11Figure 14
Nemenyi test of MAPE, MAE, andR2 (CSBA replaces the CEEMD-SE-BiGRU-AM model in the figure), (a) MAPE, (b) MAE, (c) R2.
(a)(b)(c)In the case of a small amount of data, the Nemenyi test rejected the null hypothesis that the two models are inconsistent. When the average ordinal difference between the two models is only greater than the critical value of 1.929, the null hypothesis is rejected, and there is a significant difference between the two contrasting models. It can be seen from Figure14 that the combined model established in this paper is significantly different from other models, and the differences between the BiGRU-AM model and BP and LSTM are also significant, while other models cannot directly reject the null hypothesis.Then, for the ln(RMSE) that meets the normality assumption, a pairedt-test with the same null hypothesis is performed between the two models, and its statistic is(25)t=d¯s/N,d¯ is the mean of the sample differences; N is the sample size of 10; s is the standard deviation of the sample differences.The statistic (25) follows a t-distribution with 9 degrees of freedom. It is obtained that the significance level is α=0.05 the p-value of the paired t-test between the models, as shown in Table 11.Table 11
p-value of the paired t-test.
BPLSTMBiGRUBiGRU-AMCEEMD-SE-BiGRU-AMBP10.0092<0.0001<0.00010.0004LSTM0.009210.02560.01040.0007BiGRU<0.00010.025610.07360.0036BiGRU-AM<0.00010.01040.073610.0050CEEMD-SE-BiGRU-AM0.00040.00070.00360.00501At a significance levelα=0.05, the above test can prove the following:(1)
The model after comprehensively using CEEMD and AM to improve BiGRU is significantly different from other comparable models.(2)
BiGRU already has a significant difference compared with BP and LSTM and has a significant difference with BiGRU-AM under the significance levelα=0.05.(3)
Combining the results of Tables9 and 10, it can be found that the model established in this paper is significantly better than other comparative models. AM and CEEMD have a noticeable improvement effect on the BiGRU model.


## 6. Conclusion

This paper proposes a clustering algorithm based on the GMM and WOA. On this basis, CEEMD decomposes the passenger flow sequence and then uses the AM-optimized BiGRU model to predict the short-term passenger flow. This paper uses the AFC data of some lines of the Shanghai Metro from April 1, 2015, to April 24, 2015, as a calculation example. This paper draws the following conclusions:(1)
Using WOA to optimize the parameters of the GMM, the stations are divided into 5 categories according to their spatiotemporal characteristics to save computational costs.(2)
Use the CEEMD method combined with SE to stabilize the time series and decompose and denoise the time series of inbound and outbound passenger flow in this method. The actual measurement example shows that, compared with the BiGRU-AM model without decomposition and noise reduction, the method improves the model’s measurement of MAPE by an average of 49.92%.(3)
Using BiGRU as the prediction model’s main body and the attention mechanism to capture local features of long-term sequences, the model reduces the MAPE indicators by 31.19%, 52.18%, and 27.71%, respectively, compared to BP, LSTM, and BiGRU.This paper predicts the hourly inbound and outbound passenger flow of subway stations in different geographic locations.


### 6.1. Future Work

Some further research work needs to be considered for the current research framework. For example, a model is established based on the structure of the subway station network to form a more complex spatiotemporal system to improve the efficiency of this theoretical framework and broaden its application field, but it may increase the time cost. Of course, further verification of this theoretical approach is needed.


## 6.1. Future Work

Some further research work needs to be considered for the current research framework. For example, a model is established based on the structure of the subway station network to form a more complex spatiotemporal system to improve the efficiency of this theoretical framework and broaden its application field, but it may increase the time cost. Of course, further verification of this theoretical approach is needed.



---

*Source: 7846630-2022-08-21.xml*