Recognition of Vehicles Entering Expressway Service Areas and Estimation of Dwell Time Using ETC Data

To scientifically and effectively evaluate the service capacity of expressway service areas (ESAs) and improve the management level of ESAs, we propose a method for the recognition of vehicles entering ESAs (VeESAs) and estimation of vehicle dwell times using electronic toll collection (ETC) data. First, the ETC data and their advantages are described in detail, and then the cleaning rules are designed according to the characteristics of the ETC data. Second, we established feature engineering according to the characteristics of VeESA and proposed the XGBoost-based VeESA recognition (VR-XGBoost) model. Studied the driving rules in depth, we constructed a kinematics-based vehicle dwell time estimation (K-VDTE) model. The field validation in Part A/B of Yangli ESA using real ETC transaction data demonstrates that the effectiveness of our proposal outperforms the current state-of-the-art. Specifically, in Part A and Part B, the recognition accuracies of VR-XGBoost are 95.9% and 97.4%, respectively, the mean absolute errors (MAEs) of dwell time are 52 and 14 s, respectively, and the root mean square errors (RMSEs) are 69 and 22 s, respectively. In addition, the confidence level of controlling the MAE of dwell time within 2 min is more than 97%. This work can effectively recognize the VeESA and accurately estimate the dwell time, which can provide a reference idea and theoretical basis for the service capacity evaluation and layout optimization of the ESA.


Introduction
By the end of 2021, the total mileage of expressways in China was approximately 169,100 km, ranking first in the world [1]. As an essential and critical core node of expressways, ESA is of great significance in regulating road traffic flow and relieve traffic pressure. However, the infrastructure of most ESAs built in the early years has been unable to meet the demand of the rising traffic volume, resulting in frequent queues and congestion [2,3]. Therefore, a scientific and reasonable evaluation of the service capacity of ESA and further quantitative suggestions for the reconstruction and extension of ESA have become urgent issues at present [4,5]. The pause rate and dwell time of vehicles are essential parameters in the operation and management of ESA. It is not only an important metric for operation evaluation but also a premise for layout optimization. Therefore, it is of great practical significance and application value to accurately estimate the pause rate and dwell time for the quantification of the reconstruction and extension of ESA.
Currently, the main methods for pause rate estimation are the elastic coefficient method [6,7] and feature engineering method [8][9][10][11][12][13][14][15][16][17][18]. The elasticity coefficient method, proposed by the Japanese expressway design standards, mainly uses ESA pause rate survey data and national economic data to establish the pause rate trend model to estimate the future pause rate. Although this method is simple in principle and convenient in 1.
We proposed a VR-XGBoost model for recognizing vehicles entering expressway service areas based on ETC data, which not only achieves an effective estimation of the pause rate but also accurately recognizes individual vehicles driving into ESA.

2.
Taking into full consideration the driving pattern of vehicles entering/exiting the ESA, we proposed a K-VDTE model for vehicle dwell time estimation. 3.
The validity of the proposed method is verified by using real ETC data, which can provide a more scientific and reasonable reference basis for ESA reconstruction and extension.
The remainder of this work is organized as follows: Section 2 reviews related work regarding ESA pause rate and vehicle dwell time estimation. The proposed method, including the framework, data preprocessing, feature engineering, the VR-XGBoost model, and the K-VDTE model, is described in Section 3. Section 4 shows the experimental results and analysis. Finally, the conclusion is presented in Section 5.

Pause Rate Estimation
In this section, an overview of pause rate estimation methods is presented. The elastic coefficient method (ECM) was proposed in early Japanese expressway design standards for calculating the pause rate of various types of VeESAs [6]. Drawing on relevant experience in Japan, Sun et al. [7] concluded that ECM was also applicable to the development pattern of the ESA pause rate in Guangdong Province, China, and used the ECM to estimate the average growth rate of the pause rate to achieve prediction.
Considering the close relationship between the pause rate and ESA spacing, Cui et al. [8] proposed a new method for determining the pause rate based on the continuous vehicle travel time. Through an in-depth analysis of the relationship between the pause rate and traffic flow parameters [9,10], Chen et al. [11] proposed a pause rate estimation method Entropy 2022, 24, 1208 3 of 21 based on a traditional linear regression model, which provided an important reference basis for the layout optimization and function design of ESA. In response to the low accuracy of pause rate prediction, a BP neural network-based ESA pause rate prediction model was constructed [12]. On the basis of previous work, Shen et al. [13] extracted multidimensional feature vectors from the data and constructed a tree-level BP neural network for pause rate prediction, which further improved the prediction accuracy. To further optimize the essential parameters of the wavelet neural network (WNN), some scholars introduced evolutionary algorithms, such as particle swarm optimization (PSO) [14] and genetic algorithm (GA) [15], to optimize the initial parameters of the WNN. The improved WNNbased pause rate prediction models were established, and the validity and reliability were verified on a real dataset. Under the premise of fully investigating the global optimal search capability of particles, Sun et al. [16] improved the topology of traditional PSO and fused it with the XGBoost algorithm to form a combined model for ESA traffic flow prediction. Experiments have demonstrated that the combined model has higher prediction accuracy and stronger generalization ability than a single model.
In the past few decades, deep learning methods [27,28], such as long short-term memory (LSTM) and convolutional neural networks (CNN), have achieved good performance in the field of transportation and are widely used in traffic flow prediction. Wang et al. [17] built a model based on LSTM for ESA instantaneous population analysis and prediction. The experimental results showed that it was able to accurately predict population mobility despite the relatively large population fluctuations. Zhao et al. [18] extracted spatiotemporal features using CNN, LSTM, and attention mechanism models and proposed a short-term traffic flow prediction model based on STL-OMS to achieve an accurate prediction of ESA traffic flow.

Vehicle Dwell Time Estimation
In this section, an overview of dwell time estimation methods is presented. King et al. [19] conducted an early field survey at nine locations in the United States. The results showed that the average vehicle dwell time in rest areas was 11.4 min, with a standard deviation of 12.87 min, a minimum dwell time of 1 min and a maximum dwell time of 3 h and 31 min. Recently, the Japanese Institute of Expressway General Technology noted through actual statistics that the average dwell time of small vehicles in most ESAs exceeded 25 min [20], while the dwell time for families with elderly and children was extended by an average of 10~20 min in ESAs [21]. Furthermore, analysis of dwell time by vehicle type showed that heavy vehicles had the longest average dwell time, significantly longer than other vehicle types [22]. Analysis of dwell time by seasonal differences showed that all categories of vehicles had longer dwell times in summer than in any other season [23]. Analysis of dwell time by diurnal differences showed that the average dwell time was significantly longer at night than during the day [24].
In addition, Hirai et al. [29] estimated the total dwell time in the service area for the whole trip by mining the ETC trip data using the average travel speed method. The correlation analysis of the dwell time distribution characteristics and rest behavior [30] was expected to construct the next rest behavior model. At the same time, the driver's rest behavior was used to characterize the distribution of vehicle travel time [31] to further construct driving behavior characteristics [32]. A method for calculating the number of stranded vehicles across time was proposed through statistical analysis of vehicle dwell time and rest behavior characteristics [33], and then a mathematical model for ESA scale design was proposed [34], which was used to optimize the ESA layout [35,36].

Framework
In this section, we present the framework of this study, as shown in Figure 1. First, we perform data preprocessing, including extraction of required data, ETC trajectory construction, data cleaning, data fusion and forming of structured data. Second, we Entropy 2022, 24, 1208 4 of 21 consider features such as speed features, spatiotemporal features and external features to construct feature engineering, thus building an XGBoost-based VeESA recognition model. On this basis, a kinematics-based vehicle dwell time estimation model is proposed. This study not only enables the effective recognition of VeESA but also further estimates their dwell time in the service area.

Framework
In this section, we present the framework of this study, as shown in Figure 1. First, we perform data preprocessing, including extraction of required data, ETC trajectory construction, data cleaning, data fusion and forming of structured data. Second, we consider features such as speed features, spatiotemporal features and external features to construct feature engineering, thus building an XGBoost-based VeESA recognition model. On this basis, a kinematics-based vehicle dwell time estimation model is proposed. This study not only enables the effective recognition of VeESA but also further estimates their dwell time in the service area.

Data Overview
The experimental datasets in this work contain the ETC dataset and ESA dataset. The ETC data were collected by more than 1000 ETC gantries deployed in the whole road network of the Fujian Provincial Expressway. Specifically, as the world's largest IoV system, the ETC system uses radio frequency identification (RFID) technology to enable mobile vehicles equipped with an onboard unit (OBU) to communicate with roadside units (RSU) for data collection [37]. The collection period was from 3 to 10 September 2020. We obtained a total of 42,964,489 ETC data, including vehicle ID (after desensitization), transaction time, gantry ID, vehicle type, etc., as shown in Table 1. According to the classification of vehicle types and tolls of China's expressway, vehicles can be divided into 4 categories of buses, 6 categories of trucks and 6 categories of special operating vehicles. The total number of vehicles is approximately 1.72 million in the dataset. Specifically, each transaction data contains all field information.

Data Overview
The experimental datasets in this work contain the ETC dataset and ESA dataset. The ETC data were collected by more than 1000 ETC gantries deployed in the whole road network of the Fujian Provincial Expressway. Specifically, as the world's largest IoV system, the ETC system uses radio frequency identification (RFID) technology to enable mobile vehicles equipped with an onboard unit (OBU) to communicate with roadside units (RSU) for data collection [37]. The collection period was from 3 to 10 September 2020. We obtained a total of 42,964,489 ETC data, including vehicle ID (after desensitization), transaction time, gantry ID, vehicle type, etc., as shown in Table 1. According to the classification of vehicle types and tolls of China's expressway, vehicles can be divided into 4 categories of buses, 6 categories of trucks and 6 categories of special operating vehicles. The total number of vehicles is approximately 1.72 million in the dataset. Specifically, each transaction data contains all field information.
The ESA data were collected by the cameras at the entrance and exit of Yangli ESA Part A/B. Specifically, the camera uses the technology of license plate recognition to obtain information about the vehicles entering the service area [38]. The collection period is consistent with the ETC data. We obtained more than 30,000 data points, including vehicle ID (after desensitization), capture time, service area ID and entrance/exit information, as shown in Table 2. The total number of vehicles is approximately 18,000 in the dataset. It is worth noting that this dataset is only used for experimental validation to evaluate the recognition effect and the estimation accuracy of vehicle dwell time. EnTime entrance time 2020/9/5 00:00:00 6 GantryID gantry ID G000335001000120020 7 TradeTime transaction time 2020/9/5 01:00:00 8 Workday workday 0 In this work, only the data of Yangli ESA and two ETC gantries before and after it are used, whose deployment locations are shown in Figure 2. To facilitate the following explanation, we have made relevant definitions as follows. EnTime entrance time 2020/9/5 00:00:00 6 GantryID gantry ID G000335001000120020 7 TradeTime transaction time 2020/9/5 01:00:00 8 Workday workday 0 The ESA data were collected by the cameras at the entrance and exit of Yangli ESA Part A/B. Specifically, the camera uses the technology of license plate recognition to obtain information about the vehicles entering the service area [38]. The collection period is consistent with the ETC data. We obtained more than 30,000 data points, including vehicle ID (after desensitization), capture time, service area ID and entrance/exit information, as shown in Table 2. The total number of vehicles is approximately 18,000 in the dataset. It is worth noting that this dataset is only used for experimental validation to evaluate the recognition effect and the estimation accuracy of vehicle dwell time. CapTime capture time 2020/9/5 00:00:00 In this work, only the data of Yangli ESA and two ETC gantries before and after it are used, whose deployment locations are shown in Figure 2. To facilitate the following explanation, we have made relevant definitions as follows.   Expressway Section QD [39]: Each ETC gantry and the entrance/exit of an expressway toll station is collectively called a node G, and two adjacent nodes constitute an expressway section, referred to as QD: where G 1 and G 2 are the start and end points of QD.
Taking road upline as an example, it can be seen from Figure 2 that G 1 and G 2 constitute Section 1 (QD 1 ), G 2 and G 3 constitute Section 2 (QD 2 ), where the ESA is located, and G 3 and G 4 constitute Section 3 (QD 3 ). It can be found from the partial enlarged detail that the gantries always appear in pairs, which are distributed along the upline and downline of Entropy 2022, 24, 1208 6 of 21 the road, such as G 2 and G 2 . Therefore, the discrete ETC data need to be processed into vehicle trajectories and fused with the ESA data, as detailed in Section 3.2.2.

Data Preprocessing
The prerequisite for effective data mining is to ensure data quality. However, there is a large amount of "dirty" data in ETC data, which is caused by various objective factors such as equipment failure, wireless signal crosstalk and bad weather in the process of ETC data collection, transmission and storage, which seriously affects the potential value of ETC data mining. There are 3 main problems in the "dirty" data as follows: (1) Data Redundancy Generally, it is generated by repeated uploading of data in the transmission process or repeated copying in the storage process. This tends to cause an increase in the data scale and serious interference with data mining. In addition, the continuous communication between vehicle OBUs and ETC antennas due to traffic congestion and anchor failure within the antenna coverage area is also a cause of data redundancy. In general, it is sufficient to keep only one of the instances of data and delete the rest directly.
(2) Data Missing Due to equipment failure, bad weather and other reasons, the vehicle OBU does not communicate or communicates unsuccessfully with the ETC antenna, which results in missing data. At the same time, there is also the possibility of missing data due to network packet loss during data transmission.

(3) Data Abnormality
With the influence of wireless signal crosstalk and other factors, the vehicle OBU of the vehicle traveling on the road upline communicates successfully with the ETC antenna deployed on the road downline, and the dataset generates records that do not comply with expressway driving rules.
However, due to the highly discrete characteristic of ETC data, it is difficult to achieve effective judgment of data abnormalities due to isolated data points. Therefore, it is necessary to rely on the trajectory semantic context formed by the topology of the expressway ETC gantry network to accurately detect and repair the above situation. For this purpose, we further define it as follows: ETC Trajectory eTr: The sequence of ETC gantry nodes formed by a vehicle passing through a continuous expressway Section QD 1 , QD 2 , . . . , QD n−1 is called an ETC trajectory eTr: eTr = tr 1 , tr 2 , . . . , tr n (2) where tr 1 and tr n are the start and end points of the trajectory, respectively. tr i is the transaction data when the vehicle travels through the ETC gantry, which contains information such as gantry ID tr i.N , transaction timestamp tr i.T , vehicle ID tr i.P , vehicle type tr i.C , entrance gross axle weight tr i.EW , entrance ID tr i.EID , entrance timestamp tr i.ET and workday tr i.H (consistent with Table 1). n indicates the total number of nodes that the vehicle passes through. The ETC data cleaning algorithm (Algorithm 1) includes the construction of the vehicle trajectory, data cleaning and data repair. First, the ETC data are grouped by vehicle ID tr i.P , entrance ID tr i.EID , and entrance timestamp tr i.ET . Second, we eliminate duplicate data after sorting by transaction timestamp tr i.T for each set of data. Third, we obtain two adjacent data in each set of data and judge the correctness by its topological information, which mainly includes the removal of redundant data generated by the opposite gantries and the repair of missing data. It is worth noting that the topology dataset includes two subsets: TP and TP , which is a collection of topologies (e.g., G 1 , G 2 ). Specifically, TP = { G 1 , G 2 , G 2 , G 3 , G 3 , G 4 , . . .} denotes normal topology data and TP = G 1 , G 2 , G 2 , G 3 , G 3 , G 4 , . . . denotes opposite topology data. The topologies in TP and TP always appear in pairs, such as G 1 , G 2 and G 1 , G 2 . Finally, the vehicle Entropy 2022, 24, 1208 7 of 21 trajectories that meet the requirements are added to the trajectory dataset eTRAJ. The specific algorithm is shown as follows: Algorithm 1: ETC data cleaning algorithm Input: ETC data eData, Topology data TP, Opposite topology data TP Output: ETC trajectory dataset eTRAJ 1: G_eData = eData.Groupby([tr p ,tr EID ,tr ET ]); # Grouping 2: For eTr j ∈ G_eData do: # Traversal operation for each set of data 3: eTr j ← eTr j .sorted(by = tr T ) # Sorted by transaction time 4: eTr j ← eTr j .drop_duplicates() # Data deduplication 5: While (i=1, i < len(eTr j )): IF tp ∈ TP: 8: i+ = 1; 9: continue; 10: Else IF tp ∈ TP : 11: tp ← tr i.N , tr i+2.N 12: IF : tp ∈ TP 13: delete tr i+1 # Delete opposite gantry transaction data 14: i+= 2; 15: Else : 16: IF tp ∈ TP&& tp ∈ TP: 18: The ETC driving trajectory through data cleaning also needs to be fused and matched with the service area traffic data as the label data for subsequent experiments. Therefore, we designed algorithm for fusion of ETC trajectory and ESA data (Algorithm 2). As seen from Section 3.2.1, the ESA is located in QD 2 . Therefore, only the ETC driving trajectory data and service area data vehicle data must be obtained, and at the same time, the service area entrance and exit capture time in the 2nd and 3rd gantry transaction time periods can match the VeESA to the corresponding ETC driving trajectory. The remaining unmatched driving trajectories are not driven into the service area trajectories.
Notably, the gantry system and the service area entrance/exit camera system appear to be clocked out of sync. Therefore, the time difference delta is set. We make the transaction time of G 2 ∆t hours ahead and the transaction time of G 3 ∆t hours behind, i.e., tr j 2.T − ∆t and tr j 3.T + ∆t. By expanding the time range, we ensure that VeESA is fully matched. After the experiments, the time difference in this work is set to 1 h, i.e., ∆t = 1 h. The specific algorithm is shown as follows: Through data cleaning and data fusion, a total of approximately 44,000 and 39,000 trajectories were obtained in Yangli Part A and Part B, respectively. The final data samples are shown in Table 3. In these trajectories, the total ETC trajectories of entering Part A and Part B are approximately 7800 and 6700, respectively, and the pause rates of both Parts A and B are approximately 17%. It is worth noting that due to equipment failure and other reasons, there is a missing situation of service area entrance/exit capture data in the experimental dataset. However, this problem does not affect the experiments on the recognition of VeESA in this work. In other words, only one valid capture of data needs to exist in the ESA entrance/exit data to complete the tagging work. Subsequent vehicle dwell time estimation experiments will be conducted by selecting the trajectories where both entrance and exit capture data exist.

Feature Vector Modeling
There are numerous factors that affect the pause rate and dwell time of ESA, which have highly nonlinear characteristics. Therefore, we summarize the previous research results [15,18] and construct feature vectors from 3 dimensions, such as speed features, spatiotemporal features, and external factors. The details are as follows: (1) Speed Features The speed features are the key features for the recognition of VeESA. When a vehicle enters the ESA, the average speed of ESA section (QD 2 ) will be significantly lower than QD 1 and QD 2 . Meanwhile, it will also be lower than the overall average speed of other vehicles of the same type in this section. Therefore, we construct the speed feature vector as follows: where v 1 ∼ v 3 represent the driving state of the individual vehicle during the whole trip. Among them, v 1 = d 1 /(tr 2.T − tr 1.T ) indicates the average speed of the vehicle in QD 1 , v 2 = d 2 /(tr 3.T − tr 2.T ) represents the average speed of the vehicle in QD 2 , v 3 = d 3 /(tr 4.T − tr 3.T ) represents the average speed of the vehicle in QD 3 , represent the total mileage of QD 1 , QD 2 and QD 3 , respectively, and v 4 = 1 n ∑ n j=1 v j 2 represents the overall average speed of vehicles of the same type, except for the vehicle in QD 3 . v 4 mainly avoids the disturbance caused by the reduction in v 2 in certain time periods due to special conditions or traffic congestion. (2) Spatiotemporal Features In general, the longer a vehicle spends on the expressway, the more demands on the ESA for drivers and passengers. Therefore, we construct the actual cumulative travel time of the vehicle from the entrance of the toll station to the ESA as one of the spatiotemporal features. At the same time, people's needs for ESA are also different during different times of the day and on non-workdays. For example, the pause rate of ESA is generally higher at meal times, after midnight and on non-workdays. Therefore, we construct the spatio-temporal features vector as shown below.
where γ 1 represents the actual cumulative travel time of the vehicle from the entrance of the toll station to the ESA, γ 2 represents the time period feature, which divides the whole day into 24 time periods by hour, whose value range is 0~23, and γ 3 is a variable for the workday, and its value is 0 (workday) or 1 (non-workday). (

3) External Features
Vehicle type is also an important feature in road traffic. Different types of vehicles have different demands on the ESA. At the same time, the difference in passenger/freight volume will also have some influence on the pause rate of ESA. For example, the more passengers a bus carries, the more stops it needs for rest, dining, etc. Fully loaded large trucks often require services such as breaks and water refills. Therefore, the feature vector is constructed as follows: where θ 1 represents vehicle type. From the data source of Section 3.2.1, the vehicle types are divided into 16 categories, θ 2 represents passenger/freight volume, and θ 3 represents the traffic flow of the same time slice. Feature vector modeling is completed by constructing all feature values into vector form.

Modeling of Recognition of VeESA
XGBoost is an integrated learning method based on the boosting algorithm, whose learner usually chooses the decision tree model [40], as shown in Figure 3. The model learns the residuals of the true values and the predicted values of the decision tree by iteratively generating new decision trees. Eventually, the results of all trees are accumulated as the final result to obtain better classification accuracy, i.e., the weak classifiers are combined into a stronger classifier. Therefore, we introduce XGBoost to build a VeESA recognition model. Ω( ) represents the L1 regularizer, which is used to prevent the model from overfitting.
The formula for the regularizer is Equation (9): where represents the regularization penalty coefficient, which takes values in the range of [0, 1].
presents the number of leaves of the k-th tree and represents the leaf weight of the k-th tree. The XGBoost algorithm adopts an additive stepwise integration strategy in the training process. Tree-1 is optimized first, followed by Tree-2 until Tree-K has been optimized.
We improve the prediction accuracy by adding an incremental function to optimize the objective function during the iterative process, which is calculated as in Equation (14): where represents the constant term and ( ) denotes the predicted value in the − 1st iteration on sample . We abstracted a 10-dimensional feature vector from the raw ETC data with known label information to form the sample dataset. We set the dataset as S = {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n )}, where x i = (v 1 , v 2 , v 3 , v 4 , γ 1 , γ 2 , γ 3 , θ 1 , θ 2 , θ 3 ) T (i = 1, 2, . . . , N) represents the feature vector of the i-th sample. y i = 0/1(i = 1, 2, . . . , N) represents the classification label value corresponding to x i . We assume that VR-XGBoost integrates K decision trees, and the prediction result is shown in Equation (6) where K represents the number of trees, f k (x i ) represents the predicted value of the kth decision tree on sample x i , and F represents the integrated classifier composed of all decision trees. The objective function of XGBoost consists of the loss function and the regularization item, as shown in Equation (7): Entropy 2022, 24, 1208 11 of 21 where loss represents the logistic regression loss function used for classification.
loss y i ,ŷ i = y i ln(1 + e −ŷ i ) + (1 − y i )ln(1 + eˆy i ) (8) Ω( f k ) represents the L1 regularizer, which is used to prevent the model from overfitting. The formula for the regularizer is Equation (9): where α represents the regularization penalty coefficient, which takes values in the range of [0, 1]. T k presents the number of leaves of the k-th tree and w k represents the leaf weight of the k-th tree. The XGBoost algorithm adopts an additive stepwise integration strategy in the training process. Tree-1 is optimized first, followed by Tree-2 until Tree-K has been optimized.
We improve the prediction accuracy by adding an incremental function f k to optimize the objective function during the iterative process, which is calculated as in Equation (14): where c represents the constant term andŷ i (k−1) denotes the predicted value in the k − 1st iteration on sample x i .
Next, we expand the second-order Taylor formula and discard the constant term to speed up the solution and reduce the running time, which is calculated as Equation (15): where I j = {i|q(x i ) = j} denotes the sample set of leaf j, and g i and h i are the first derivative and the second derivative of the loss function, respectively. The objective function is transformed into a quadratic Obj (k) minimization problem on w j . Then, we obtain the optimal prediction of each leaf node and the minimum value of the objective function, that is, the optimal value: where G j = ∑ i∈I j g i , H j = ∑ i∈I j h i.

Kinematics-Based Dwell Time Estimation
The vehicle dwell time is also an essential parameter in the operation and management of ESA. Therefore, after the recognition of VeESA, we need to further estimate the dwell time in the ESA. From the location of the service area between Gantry 2 and Gantry 3, we know that the total travel time of the section consists of the actual travel time of the vehicle in the section and the dwell time. Therefore, the dwell time ∆t s can be obtained as follows: where ∆t QD2 = tr 3.T − tr 2.T , ∆t r represents the actual travel time, which is an unknown parameter. Therefore, the vehicle dwell time estimation is converted into the actual vehicle travel time estimation. Since the traffic conditions of the expressway are relatively smooth, the expressway can approximate the free-flow state in noncongested and nonemergency conditions. Vehicles usually travel smoothly on the highway, so the average speed of QD 1 and QD 3 can be used as the speed of QD 2 , and thus the actual travel time of QD 2 can be estimated: By substituting Equation (19) into Equation (20), we can obtain the following: Although the average speed method is simple and straightforward, it does not take into account the kinematics of the VeESA during the entrance/exit ramp. In general, VeESA goes through a total of five kinematic stages, including smooth driving upstream, decelerating into the ESA, dwelling in the service area, accelerating out of the ESA and smooth driving downstream, as shown in Figure 4.  Therefore, we construct a kinematics-based model for estimating the dwell time, where the actual travel time ∆ is redefined as follows: where ∆ ~∆ correspond to the time spent in each of the above five stages.
Stage 1: smooth driving upstream According to the principle of inertia, the driving state of this stage can be considered as the continuation of the previous section ( ). Therefore, we approximate Stage 1 as uniform motion. We take the average travel speed of as the travel speed of Stage 1, and we can obtain the time spent in Stage 1: where ∆ is the distance from Gantry 2 to the diversion point of the entrance ramp and Therefore, we construct a kinematics-based model for estimating the dwell time, where the actual travel time ∆t r is redefined as follows: where ∆t 1 ∼ ∆t 5 correspond to the time spent in each of the above five stages. Stage 1: smooth driving upstream According to the principle of inertia, the driving state of this stage can be considered as the continuation of the previous section (QD 1 ). Therefore, we approximate Stage 1 as Entropy 2022, 24, 1208 13 of 21 uniform motion. We take the average travel speed of QD 1 as the travel speed of Stage 1, and we can obtain the time spent in Stage 1: where ∆s 1 is the distance from Gantry 2 to the diversion point of the entrance ramp and v 1 is the average travel speed of QD 1 .
Stage 2: decelerating into the ESA To a certain extent, the ramps in ESA are similar to the ramps at the entrance and exit of the expressway toll station. However, the ramps at the entrances and exits of toll stations are usually designed with large curvature, while the ramps at service areas are generally of small curvature or even similar to straight lines. This makes the vehicle smoother when driving in/out of the service area. Therefore, we approximate Stage 2, i.e., the deceleration driving process into the service area entrance ramp, as uniform deceleration linear motion. From stage 1, the initial velocity of uniformly decelerating linear motion is v 1 . Let the velocity at the moment ∆t 2 be v ∆t 2 , the displacement be ∆s 2 and the acceleration be a − , which gives the following. v We combine Equations (23) and (24) to obtain the following.
where ∆s 2 is the distance from the diversion point of the entrance ramp to the service area.
Stage 4: accelerating out of the ESA Similarly, we approximate Stage 4, i.e., the service area exit ramp acceleration process, as uniformly accelerated rectilinear motion until the driving speed reaches a steady state. From Stage 4, it can be seen that the vehicle reaches a smooth state after moment ∆t 4 , whose driving speed is v 3 . Meanwhile, we assume the initial velocity v 30 and acceleration a + of uniformly accelerated linear motion. From Equation (23), we can obtain the time spent in Stage 4: In general [41], a + = 0.8 ∼ 1.2 m·s −2 . Stage 5: smooth driving downstream This stage is similar to the smooth driving upstream. Therefore, we approximate Stage 5 as uniform motion. The driving state of the next section (QD 3 ) can be considered a continuation of Stage 5. We take the average travel speed of QD 3 as the travel speed of Stage 5, and we can obtain the time spent in Stage 5 as follows: where v 3 is the average travel speed in the back section of the service area, and ∆s 5 is the distance from the smooth point in Stage 4 to Gantry 3, which is expressed as follows: We substitute Equations (22) and (25)- (27) into Equation (21) to obtain the following.
After finishing, we obtain the mathematical model for the estimation of vehicle dwell time based on kinematics.
It can be generally considered that the velocity v ∆t 2 in uniformly decelerating linear motion and the initial velocity v 30 in uniformly accelerating linear motion are both zero, which can be simplified as follows:

Experiments
The experimental platform is a Centos Linux release 7.9.2009 (Core) operating system based on an Intel(R) Core (TM) i9-10900K CPU @ 3.70 GHz and 64 GB RAM, and all experiments were implemented on the open-source web application Jupyter Notebook using Python version 3.8.8.

Construction of Feature Vector
We constructed the feature vector dataset for the training of VR-XGBoost by using 10 statistical features, as shown in Table 4. In the feature vector dataset, each vector contains 10 dimensions of attributes and its classification label l, where l = 0 represents non-VeESA and l = 1 represents VeESA. It is worth noting that the cumulative travel time γ 1 is not directly available in the ETC data. We replaced it with the cumulative travel time from the entrance of the toll station to the front gantry of ESA, i.e., γ 1 is the cumulative travel time from the entrance of the entrance to G 2 . The correlation heatmap is further inscribed for correlation analysis of the feature vectors, as shown in Figure 5. In the figure, blue indicates a positive correlation between vectors, and red indicates a negative correlation between feature vectors. At the same time, when the color is more prominent, the correlation between vectors is stronger. The speed features were positively correlated with the traffic flow and negatively correlated with the cumulative travel time γ 1 , vehicle type θ 1 and entrance gross axle weight θ 2 . Specifically, there is a strong positive correlation among v 1 , v 3 and v 4 , which are both weakly positively correlated with v 2 . The two features of γ 2 and γ 3 have a very low correlation with other features. Through heatmap analysis, we can clearly understand the correlation between feature vectors.
The correlation heatmap is further inscribed for correlation analysis of the feature vectors, as shown in Figure 5. In the figure, blue indicates a positive correlation between vectors, and red indicates a negative correlation between feature vectors. At the same time, when the color is more prominent, the correlation between vectors is stronger. The speed features were positively correlated with the traffic flow and negatively correlated with the cumulative travel time , vehicle type and entrance gross axle weight . Specifically, there is a strong positive correlation among , and , which are both weakly positively correlated with . The two features of and have a very low correlation with other features. Through heatmap analysis, we can clearly understand the correlation between feature vectors.

Parameters Selection
The XGBoost classification algorithm has numerous parameters, including the following three aspects.
(2) Booster Parameters: the number of decision trees (n_estimators), learning rate (learn_rate), maximum depth of the tree (max_depth), minimum weight in leaf nodes

Parameters Selection
The XGBoost classification algorithm has numerous parameters, including the following three aspects.
(1) General Parameters: booster, silent, nthread, etc. The general parameters and learning task parameters are set directly according to the model needs, while the booster parameters should be parameter-seeking by the tuning method. At present, the tuning method is mainly the grid search method, which is combined with the K-fold cross-validation method to achieve the optimal parameters [42]. In this work, we also used this method for tuning the parameters and set the cross-validation parameter K = 5. The search range, step size and optimal values of parameters for each parameter are shown in Table 5.

Comparative Analysis of Classification Models
To comprehensively evaluate the effectiveness of the VR-XGBoost model, this work introduced evaluation metrics such as accuracy, precision, recall, and F1-score, as shown below.
We compared and analyzed this experimental model with commonly used machine learning models (e.g., RF, GBDT, KNN), and the experimental results are shown in Table 6.
The experimental results showed that VR-XGBoost, RF and GBDT all obtained good recognition results with accuracy above 95%, while DT performed the worst due to its tendency to overfit. In particular, VR-XGBoost achieved the best results in evaluation metrics. Specifically, in Part B, the accuracy of VR-XGBoost was as high as 97.4%. This result showed the significant superiority of the VR-XGBoost model for the recognition of VeESA. Next, the feature contributions are further analyzed, as shown in Figure 6. As a whole, the feature contribution ranking from largest to smallest is speed features, external features, and spatiotemporal features. In particular, the contribution of the speed feature in Part A and Part B, both of which exceed 65%, is much higher than that of the spatiotemporal feature and external features. Specifically, the feature contribution of v 2 in speed features is more than 50%, indicating that the feature is the most important. In contrast, the contribution rates of features, such as the actual cumulative travel time γ 1 , the time period feature γ 2 , the passenger/freight volume θ 2 , and the traffic flow θ 3 , etc., are all less than 5%. These features seem less important. Through quantitative analysis of contribution rate, we can clearly know the importance of each feature.

K-VDTE Evaluation
We sliced in 5-min increments to count the dwell time, and the distribution is shown in Figure 7. It can be seen that the dwell times in Part A and Part B both exhibit a long-tailed distribution, which indicates that most vehicles stay in the ESA only temporarily and briefly. Specifically, the number of vehicles with a dwell time of 5~10 min is the greatest, and more than 90% of the vehicles have a dwell time of less than 1 h in the ESA. Furthermore, the average dwell time in Part A and Part B was approximately 30 min, with a standard deviation of approximately 70 min, a minimum dwell time of less than 30 s and a maximum dwell time of more than 12 h. Through the statistical analysis, we can clearly understand the general situation of the vehicle dwell time in each ESA.
Next, the feature contributions are further analyzed, as shown in Figure 6. As a whole, the feature contribution ranking from largest to smallest is speed features, external features, and spatiotemporal features. In particular, the contribution of the speed feature in Part A and Part B, both of which exceed 65%, is much higher than that of the spatiotemporal feature and external features. Specifically, the feature contribution of in speed features is more than 50%, indicating that the feature is the most important. In contrast, the contribution rates of features, such as the actual cumulative travel time , the time period feature , the passenger/freight volume , and the traffic flow , etc., are all less than 5%. These features seem less important. Through quantitative analysis of contribution rate, we can clearly know the importance of each feature.

K-VDTE Evaluation
We sliced in 5-min increments to count the dwell time, and the distribution is shown in Figure 7. It can be seen that the dwell times in Part A and Part B both exhibit a longtailed distribution, which indicates that most vehicles stay in the ESA only temporarily and briefly. Specifically, the number of vehicles with a dwell time of 5~10 min is the greatest, and more than 90% of the vehicles have a dwell time of less than 1 h in the ESA. Furthermore, the average dwell time in Part A and Part B was approximately 30 min, with a standard deviation of approximately 70 min, a minimum dwell time of less than 30 s and a maximum dwell time of more than 12 h. Through the statistical analysis, we can clearly understand the general situation of the vehicle dwell time in each ESA. To evaluate the effectiveness of the K-VDTE model, the estimation errors are quantified using the evaluation metrics of root mean square error (RMSE), mean absolute error (MAE), and R coefficient, as shown below: To evaluate the effectiveness of the K-VDTE model, the estimation errors are quantified using the evaluation metrics of root mean square error (RMSE), mean absolute error (MAE), and R coefficient, as shown below: whereŷ i denotes the estimated dwell time obtained using the model, y i denotes the true dwell time, y i is the average dwell time, and n denotes the amount of data. We compared the proposed K-VDTE model with the traditional averaging speed method and commonly used machine learning models, as shown in Table 7. The experimental results show that the proposed K-VDTE method performs the best, while the machine learning model performs the worst. Specifically, taking Part B as an example, the MAE of the K-VDTE model was only 14 s, which was not only more than one times better than the average speed method but also at least four times better than the machine learning models. This demonstrated the higher accuracy of our method. Moreover, comparing the RMSE of each model, the proposed K-VDTE model improved at least one order of magnitude over the machine learning models, which indicated that the proposed method was more robust. Moreover, the integrated learning models, such as XGBoost, RF and GBDT, among machine learning models, perform better on evaluation metrics, while the single models, such as Lasso and KNN, obtain very poor results on all evaluation metrics. This result indicates that single models may not be suitable for vehicle dwell time estimation.
To investigate the estimated errors in depth, we performed a statistical analysis of the MAE of the dwell times and carved out the distribution of the cumulative probabilities, as shown in Figure 8. It can be seen that the distribution curves of the cumulative probabilities in Part A and Part B all exhibited a rapid increase with the increase in the dwell time estimation error until they stabilized after 2 min. Specifically, P{MAE ≤ 120 s} > 95% indicates that the probability of keeping the MAE within 2 min is more than 95%. Specifically, taking Part B as an example, the probability of controlling the MAE within 1 min and 2 min are P B {MAE ≤ 60 s} > 97% and P B {MAE ≤ 120 s} > 99.8%, respectively. The results further validate that the K-VDTE model has strong robustness.
ties, as shown in Figure 8. It can be seen that the distribution curves of the cumulative probabilities in Part A and Part B all exhibited a rapid increase with the increase in the dwell time estimation error until they stabilized after 2 min. Specifically, { ≤ 120 s} > 95% indicates that the probability of keeping the MAE within 2 min is more than 95%. Specifically, taking Part B as an example, the probability of controlling the MAE within 1 min and 2 min are { ≤ 60 s} > 97% and { ≤ 120 s} > 99.8%, respectively. The results further validate that the K-VDTE model has strong robustness.

Conclusions
In this work, we proposed a method for the recognition of vehicles entering expressway service areas and the estimation of dwell time based on ETC data. This method provides reference ideas for scientific and reasonable evaluation of the service capacity of the ESA, which can also provide decision support for the optimization of the layout when reconstructing and extending the ESA. The specific conclusions are as follows:

Conclusions
In this work, we proposed a method for the recognition of vehicles entering expressway service areas and the estimation of dwell time based on ETC data. This method provides reference ideas for scientific and reasonable evaluation of the service capacity of the ESA, which can also provide decision support for the optimization of the layout when reconstructing and extending the ESA. The specific conclusions are as follows: (1) Experiments were conducted using real ETC data with a user penetration rate of over 80%. It not only solves the issue of insufficient data volume but also solves the geographical differences existing in different service areas in vehicle dwell time estimation. It can provide a more scientific and reasonable reference basis for the evaluation of the service capacity of ESA.
(2) Considering multidimensional information such as speed features, spatiotemporal features and external features, we constructed a VR-XGBoost model. This model can achieve not only the estimation of the overall pause rate of ESA but also the accurate recognition of vehicles entering the service area. (3) After an in-depth study of the driving pattern of vehicles in the process of driving in/out of the ESA, we proposed a K-VDTE to realize vehicle dwell time estimation. The estimation accuracy of vehicle dwell time can be further improved by considering vehicle kinematics.
However, the present method also has certain limitations, whose expressway traffic state must approximate free-flow conditions. In the future, we will further explore the vehicle driving characteristics and laws under nonfree flow conditions to form a more scientific and reliable evaluation system. Author Contributions: Conceptualization, F.Z. and Q.C.; methodology, Q.C.; software, Q.C., Z.Z. and N.L.; validation, D.Y., Q.C. and F.Z.; formal analysis, Q.C. and F.Z.; investigation, F.G. and Q.C.; resources, F.Z.; data curation, Q.C. and D.Y.; writing-original draft preparation, Q.C., Z.Z. and N.L; writing-review and editing, F.Z., Q.C., D.Y., Z.Z. and N.L; visualization, Z.Z. and N.L.; supervision, Q.C.; project administration, F.G. and Q.C.; funding acquisition, F.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
Restrictions apply to the availability of these data. Data were obtained from Fujian Expressway Information Technology Co., Ltd. and are available from the authors with the permission of Fujian Expressway Information Technology Co., Ltd.

Conflicts of Interest:
The authors declare no conflict of interest.