Learning traffic as videos: A spatio-temporal VAE approach to periodic traffic raster data imputation

2.1 Prediction-based method

The prediction-based method tends to build a connection between historical and future missing value. Auto-regressive integrated moving average (ARIMA), one of the most popular prediction models for time series, is utilized in [7] for the imputation of traffic counts for highway agencies. [10] presents an approach to selecting and estimating ARIMA models for traffic volume imputation. However, the shortcoming of this method is obvious that the information following the missing data are not employed during the imputation, making it less effective.

2.2 Interpolation-based method

Interpolation-based approach fills the missing value with the neighboring value. As for spatio-temporal data, neighboring value can either be the temporal-liking value or spatial-liking value. [8] proposes a data-driven imputation method for sections of road based on their spatial and temporal correlation using a modified K-Nearest Neighbor method. [9] compares Local Least Squares (LLS) with K-Nearest when dealing with missing traffic data imputation. Both the research mentioned above use interpolation-based approaches. Nevertheless, this kind of method is determinative and cannot extract stochastic variations in spatio-temporal data.

2.3 Decomposition-based method

The decomposition-based approach intends to decompose the original tensor into multiple lower-dimension tensors, and then fills in the missing value on the lower-dimension tensors. Liu J, Musialski P, Wonka P, et al first come up with High Accuracy Low Rank Tensor Completion (HaLRTC), which extends the matrix case to the tensor case by proposing the definition of the trace norm for tensors for the first time and then building a working algorithm [11, 12] estimates traffic speed for ITS by combining HaLRTC and k-means++. Though the accuracy of this kind of approach is high, the computation complexity is sometimes unbearable.

2.4 Statistic-based method

The statistic-based approach learns the statistical rule of the data directly, and undertakes the imputation task based on the rule. The most widely-used method is Principal Component Analysis (PCA), and many improvements have been conducted on it [13] employs PPCA (Probabilistic Principal Component Analysis) to conduct probabilistic interpretation on PCA when dealing with traffic flow value imputation, while [14] adds Bayesian network to PCA to solve the same problem, which is known as BPCA (Bayesian Principal Component Analysis) [15] uses KPPCA (Kernel Probabilistic Principal Component Analysis) as imputing methods without any temporal or spatial dependency, and replaces linear kernel with non-linear kernel in PCA. However, this kind of approach focuses too much on constructing probability distribution model, thus loses detailed information varying in different situations.

2.5 Deep learning method

In recent years, the deep Learning methods have made noticeable achievement from all aspects, and various DL approaches have been improved to deal with data imputation issue [16] utilizes BRITS, an approach based on recurrent neural networks for imputing missing values in time series data without any specific assumption. Though the bidirectional temporal dependency within time series is considered during data imputation, little attention is given to spatial dependency among different time series [17] first transforms raw data into spatial-temporal images and then implements data imputation task with the help of convolution neural network (CNN). In [18], a deep auto-encoder for missing spatio-temporal data imputation is put forward, and the auto-encoder is designed as a combination of CNN and Bi-LSTM (Bi-Long Short Term Memory). Moreover, [19] deploys CNN based GAN (Generative Adversarial Network) to reconstruct missing multivariate time series. All these deep learning methods are limited to 2D-convolution when dealing with spatial feature extraction, and long-short term temporal dependencies are not considered.

3.1 Missing data type

According to [20], missing data can be summarized into three categories: Missing Completely at Random (MCAR), Missing at Random (MAR) and Missing Not at Random (MNAR).

If the missing is irrelevant to observed or unobserved data, that is MCAR, which may happen in any region at any time randomly. MAR refers to the missing that can be described using observed data. For example, in-flow data missing in one region may be related to the out-flow data missing in adjacent regions. When the missing depends on unobserved attribute or on the missing attribute itself [21], we call it MNAR. MNAR occurs when devices are in long-term failure. The study mainly focuses on the former two kinds of missing, for it is of little significance to populate the missing value in MNAR.

3.2 Traffic raster data

As [20] depicts, an urban traffic network can be regarded as a raster map G with longitude and latitude containing H rows and W columns. Each grid in G means one region. For each region ( h , w ) , C different channels are measured at each timestamp. In this paper, C is set to 2 for in-flow and out-flow of one region. Therefore, if the map G is considered as a whole, the measurements over a period of time, which involves T intervals, can be considered as a tensor x ∈ ℝ C * T * H * W .

3.3 Problem statement

Under the MCAR assumption that data are missing completely at random [21], it seems that incomplete traffic raster data that are i.i.d. can be described as a set 𝑿 = { x i } i = 1 n ∈ ℝ C ∗ T ∗ H ∗ W . Correspondingly, another set 𝑴 = { m i } i = 1 n ∈ M n is defined as a mask set for 𝑿 , where M = { 0 , 1 } C ∗ T ∗ H ∗ W when m i ⁢ [ c , t , h , w ] = 0, whereas x i ⁢ [ c , t , h , w ] is missing, m i ⁢ [ c , t , h , w ] = 1.

For tensor-based traffic raster data imputation that this paper discusses on, an incomplete tensor x in the space 𝑿 is divided into an observed part x 𝑜𝑏𝑠 and a missing part x 𝑚𝑖𝑠 . The objective of this paper is to populate the missing value in x 𝑚𝑖𝑠 through x 𝑜𝑏𝑠 .

1) Pre-processing. It is known that traffic raster data show strong hourly, daily, weekly, and yearly periodicity when regarded as a whole [40]. In order to reflect this multi-periodic feature (mentioned in Section 1), the original traffic raster data are segmented into two kinds of samples 𝑿 l , 𝑿 s with different sizes (i.e. hour and day), corresponding to short-term periodicity and long-term periodicity. Then two kinds of samples are sent into two streams (as shown in Fig. 1) to learn short-term periodicity and long-term periodicity features respectively. However, this model is just an illustrative example and can be generalized to more streams as needed. The segment size is recommended to follow the natural time period. For example, supposing that the sampling frequency of sensor is q times per hour, then the segment size for short-term periodicity can be assumed as q and that for long-term periodicity is 24 * q correspondingly. In other words, the former samples represent hourly data, and the latter represent daily data. More details about pre-processing of datasets used in this paper will be introduced in Section 5.1.

2) Encoder. After Pre-processing, both two streams of traffic raster data transfer to the Encoder part and flow along four ECBs (elaborated in Section B) in sequence to generate the deterministic hidden representations denoted as h . Then h is flattened before entering the fully-connected layers to parameterize the mean μ and log deviation 𝑙𝑜𝑔 ⁢ σ 2 for the original latent variable z 0 , which will be developed by the reparametrization trick [21] as follows, where σ is sampled from 𝒩 ⁢ ( 0 , 1 ) .

z 0 = μ + σ ⊗ ε(1)

3) Sylvester NFs. However, the output of the encoder z 0 obeys the probabilistic posterior q 0 ( z 0 | x ) ∼ 𝒩 ( μ , σ 2 I ) , which may limit the expression of vanilla VAE [21]. Therefore, Berg et al. [6] proposed the Sylvester NFs which could enhance the scalability of the distribution q ψ ( z | x ) . A series of convertible mappings f K ∘ ⁢ … ∘ ⁢ f k ∘ ⁢ … ∘ ⁢ f 1 ⁢ ( z 0 ) are utilized to translate z 0 to z k , where f k ( z k - 1 ) = z k - 1 + 𝑄𝑅 tan ( R ∼ Q T z k - 1 + b ) . R and R ′′ are upper triangular parameter matrices, whereas Q and b are column orthogonal matrix and parameter vector respectively. All the parameters can be calculated by a dependent dense layer. The log posterior log q K ( z K | x ) of final iteration can be expressed in the Eq. (2), where det ⁢ ( ∂ ⁡ f k ⁢ ( z k - 1 ) ∂ ⁡ z k - 1 ) is the Jacobian determinant of the k-transformation f k . It is assumed that q φ ( z | x ) := q K ( z K | x ) and then z K is fed into the decoder part.

log q K ( z K | x ) = log q 0 ( z K | x ) - ∑ k = 1 K log | det ( ∂ ⁡ f k ⁢ ( z k - 1 ) ∂ ⁡ z k - 1 ) |(2)

4) Decoder. The decoder deploys single dense layer and reshaping operation to transform z K back to a tensor, and then utilizes three TCBs for up-sampling (elaborated in Section B). Then, the reconstructed sample can be obtained through the last 3D transposed convolution and tanh function. In the end, reconstructed samples X l ^ and X s ^ derived from the two streams are fused together to obtain the final reconstructed traffic raster data X ^ by the Eq. (3), where W l and W s are parameter matrices that need to be trained through the whole network.

X ^ = W l ⊙ X ^ l + W s ⊙ X s ^(3)

1) Extraction-and-Calibration Block. As the major component in the encoder part, Extraction-and-Calibration Block (ECB), mainly involves two sub-modules: dilated P3D gated convolution and multi-horizon attention mechanism. Assume that a feature map F l - 1 ∈ ℝ C ∗ T ∗ H ∗ W is the input of l-th ECB, by applying dilated P3D gated convolution 𝑮 , ECB initially extracts spatial-temporal information in F l - 1 and pays more attention to the valid ones to extract important features. Next, it further calibrates the feature maps by utilizing multi-horizon attention mechanism that combines channel attention M c , temporal attention M t with cross-spatial attention M c ⁢ s . Therefore, it is named as Extraction-and-Calibration Block as shown in Fig. 2. The overall procedure of ECB is explained in the following equations:

F l - 1 ′ = G ⁢ ( F l - 1 )(4) F l = M c ⁢ ( F ′ l - 1 ) ⊙ M t ⁢ ( F ′ l - 1 ) ⊙ M c ⁢ s ⁢ ( F ′ l - 1 ) ⊙ F ′ l - 1(5)

where F l - 1 ′ ∈ ℝ C ′ * T ′ * H ′ * W ′ is an intermediate feature map, while F l ∈ ℝ C ′ * T ′ * H ′ * W ′ is the output of l-th ECB and ⊙ denotes element-wise product. M c , M t , and M c ⁢ s can be calculated in Eqs (6), (7), and (8) respectively. Furthermore, the dilated P3D gated convolution and multi-horizon attention mechanisms are elaborated in the following.

OPEN IN VIEWER

OPEN IN VIEWER

a) Dilated P3D gated convolution. In order to extract spatio-temporal features from data (like video), one natural way is to extend convolution kernels in CNN from 2D to 3D, like what [26] does to deal with video in-painting issue. To be more specific, in each layer of 3D CNN, a set of N maps of size D * H * W are convoluted by M sets of K D * K H * K W size filters, as shown in Fig. 3. However, traditional 3D CNN requires huge amount of computation and memory overheads [28], thus [27] comes up with a novel Pseudo-3D Residual Net (P3D ResNet) to solve this problem. The key point of P3D lies in the utilization of 2D convolutions to encode spatial information and 1D convolutional filters for temporal dimension. Moreover, [27] proposes three modes for P3D, as depicted in Fig. 3. MSST-VAE is mainly based on design (a), which considers stacked architecture by making temporal 1D filters (T) follow spatial 2D filters (S) in cascaded manner [27].

Different from [27], in dilated P3D gated convolution, two parts (S and T) of P3D CNN are optimized respectively to capture global-local-granularity spatial dependencies and long-short-term temporal dependencies. Besides, an additional gating mechanism is utilized as an attention map on the output features.

OPEN IN VIEWER

The workflow of dilated P3D gated convolution is displayed in Fig. 4, in which the upper part is for nonlinear spatio-temporal feature learning while the bottom part is used as a gating for soft weights calibration to the upper part by a sigmoid function. Finally, both parts are fused by the element-wise product. Dilated convolution was first proposed in [29] to increase the receptive filed in convolution operation.

It should be noticed that spatial 2D filters of P3D CNN in [27] are replaced with 2D-dilated convolution to extract global-local-granularity spatial features, and temporal 1D filters are substitute for 1D-dilated convolution for long-short-term temporal features. More details about these two parts will be discussed later.

According to [38], treating traffic raster data as images and simply applying plain CNN may not achieve the best performance, for correlations among ‘far away’ regions in the image are not taken into account. Inspired by [39], 2D-dilated spatial convolution (i.e., the blue block in Fig. 5) is utilized in dilated P3D gated convolution to catch local-granularity and global-granularity spatial features. Its structure is shown in Fig. 6. For each layer of ECB, traffic raster images with time interval t are convoluted with different kernel size. Different feature maps derived from various sizes of kernels are concatenated at channel dimension. Then, the concatenated feature maps will be fed into the next layer with the same set of kernel size but a smaller dilation factor, so on and so forth.

OPEN IN VIEWER

According to [37], time series is characterized by its numerical and continuous nature, so it is always seen as a whole instead of individual numerical fields. Traffic raster data can be viewed as 2D time series. Therefore, 1D-dilated temporal convolution (i.e. the yellow block in Fig. 5) is utilized in dilated P3D gated convolution to catch long-term and short-term temporal features. Its structure is similar to 2D-dilated temporal convolution, as shown in Fig. 5. Feature maps derived from 2D-dilated spatial convolution is convoluted by 1D kernel at temporal dimension. Kernel sizes are the same among most different layers, but differ in one certain layer. When the network goes deeper, the dilation factor gets smaller gradually.

b) Multi-horizon attention mechanism. Channel attention and spatial attention in Convolution Block Attention Module (CBAM) [22] are employed to image classification task. In this study, we extend CBAM from spatial attention to cross-spatial attention, and temporal attention is also added as shown in Fig. 6. In the multi-horizon attention mechanism, channel attention is utilized to learn the channel relations, such as inflows and outflows of traffic crowd, while temporal attention and cross-spatial attention are employed to enlarge or weaken the influence of features computed by dilated P3D gated convolution. Moreover, max-pooling and average-pooling are introduced to accelerate the computing as well as avoiding over-fitting.

OPEN IN VIEWER

For channel attention, we first perform max-pooling and average-pooling operations to generate two channel-wise variables v l - 1 max ∈ ℝ C ∗ 1 ∗ 1 ∗ 1 and v l - 1 𝑎𝑣𝑔 ∈ ℝ C ∗ 1 ∗ 1 ∗ 1 from feature map F l - 1 ′ . Then both two variables are fed into two com mon dense layers independently, one for channel reduction from C to C / r ( r is the reduction ratio), the other for channel recovery back to C . Finally, two variables are fused and calculated by a sigmoid function to get the channel attention map M c ∈ ℝ C ∗ 1 ∗ 1 ∗ 1 . The computation formulas are as follows:

M c = 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 ⁢ ( 𝑀𝐿𝑃 ⁢ ( 𝑀𝑎𝑥𝑃𝑜𝑜𝑙𝑖𝑛𝑔 ⁢ ( F l - 1 ′ ) ) + 𝑀𝐿𝑃 ⁢ ( 𝐴𝑣𝑔𝑃𝑜𝑜𝑙𝑖𝑛𝑔 ⁢ ( F l - 1 ′ ) ) ) = 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 ⁢ ( 𝑀𝐿𝑃 ⁢ ( v l - 1 max ) + 𝑀𝐿𝑃 ⁢ ( v l - 1 𝑎𝑣𝑔 ) )(6)

For temporal attention, the architecture of which is the same with that of channel attention, and the temporal attention map M t ∈ ℝ 1 ∗ T ∗ 1 ∗ 1 can be calculated correspondingly as follows:

M t = 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 ⁢ ( 𝑀𝐿𝑃 ⁢ ( 𝑀𝑎𝑥𝑃𝑜𝑜𝑙𝑖𝑛𝑔 ⁢ ( F l - 1 ′ ) ) + 𝑀𝐿𝑃 ⁢ ( 𝐴𝑣𝑔𝑃𝑜𝑜𝑙𝑖𝑛𝑔 ⁢ ( F l - 1 ′ ) ) ) = 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 ⁢ ( 𝑀𝐿𝑃 ⁢ ( v l - 1 max ) + 𝑀𝐿𝑃 ⁢ ( v l - 1 𝑎𝑣𝑔 ) )(7)

For cross-spatial attention, max-pooling and average-pooling are operated on feature map F l - 1 ′ to develop two 2D maps u l - 1 max ∈ ℝ 1 ∗ 1 ∗ H ∗ W and u l - 1 𝑎𝑣𝑔 ∈ ℝ 1 ∗ 1 ∗ H ∗ W . Next, both two are concatenated in channel dimension and sent to cross convolution part followed by a sigmoid function to create a cross-spatial attention map M c ⁢ s ∈ ℝ 1 ∗ 1 ∗ H ∗ W .

M c ⁢ s = 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 ⁢ ( 𝐶𝑟𝑜𝑠𝑠𝐶𝑜𝑛𝑣 ⁢ ( [ 𝑀𝑎𝑥𝑃𝑜𝑜𝑙𝑖𝑛𝑔 ⁢ ( F l - 1 ′ ) , 𝐴𝑣𝑔𝑃𝑜𝑜𝑙𝑖𝑛𝑔 ⁢ ( F l - 1 ′ ) ] ) ) = 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 ⁢ ( 𝐶𝑟𝑜𝑠𝑠𝐶𝑜𝑛𝑣 ⁢ ( [ u l - 1 max ; u l - 1 𝑎𝑣𝑔 ] ) )(8)

c) Cross-Convolution. Cross-Convolution mainly consists of two operations: Lat-Conv and Lon-Conv, just as Fig. 7 depicts. The former is utilized for learning spatial attention of latitude dimension and the latter is for spatial attention of longitude dimension, which is typical for traffic raster data. A 𝑙𝑎𝑡 and A 𝑙𝑜𝑛 can be calculated by the following:

F l - 1 ′′ = α ⋅ 𝐿𝑎𝑡𝐶𝑜𝑛𝑣 ⁢ ( F l - 1 ′ ) + F l - 1 ′(9) 𝐶𝑟𝑜𝑠𝑠𝐶𝑜𝑛 ⁢ ( F l - 1 ′ ) = β ⋅ 𝐿𝑜𝑛𝐶𝑜𝑛𝑣 ⁢ ( F l - 1 ′′ ) + F l - 1 ′′(10)

where α and β are scalar value, LatConv and LonConv are 2D convolution with kernel size of 1 * W and 1 * H respectively. Eventually, the feature map F ⁢ l of next layer can be obtained by element-wise product of 𝑀𝑐 , 𝑀𝑡 , 𝑀𝑐𝑠 and F l - 1 ′ , just as explained in Eq. (5) above.

OPEN IN VIEWER

2) Transposed Convolution Block. Transposed Convolution Block (TCB) is the main component for up-sampling in Decoder part of MSST-VAE. Compared with ECB, the architecture of TCB is relatively simpler, for it only consists of P3D gated transposed convolution, without multi-horizon attention mechanism. The workflow of P3D gated transposed convolution is similar with Fig. 5, and the convolution operation is replaced with the transposed convolution operation.

The application of MSST-VAE in engineering practice can be divided into two phrases: offline model training and online imputation. The former is for training the model to converge, and the latter is to utilize the trained model to impute missing value in real time. Therefore, both share the same network architecture and the models as well as parameters of the latter are derived from the former.

1) Offline Model Training. Since the true value of missing data are unavailable, the objective of MSST-VAE is to minimize the negative log likelihood of the observed data in offline model training. In Fig. 1, each stream is trained separately and then the fusion parameter in Decoder is fine-tuned as expressed in Eq. (3). For the i-th sample x i in each stream, the loss function is designed as in the following, where ( θ , ϕ ) are parameterized by the Encoder and Decoder respectively. In Eq. (11), the first item is the reconstruction error of the observed part x i 𝑜𝑏𝑠 , while the second item is from SNFs. The last item is the Kullback-Leibler distance [30] between two densities to avoid over-fitting by making conditional distribution of z 0 similar to the prior distribution of z K .

F i ( θ , ϕ ) = - ( log p ( x i 𝑜𝑏𝑠 , z k ; θ ) - log p ( z k | x i 𝑜𝑏𝑠 ; θ ) ) = - log p ⁢ ( x i 𝑜𝑏𝑠 , z k ; θ ) q k ( z k | x i 𝑜𝑏𝑠 ; ϕ ) + log p ( z k | x i 𝑜𝑏𝑠 ; θ ) q k ( z k | x i 𝑜𝑏𝑠 ; ϕ ) = 𝔼 q k [ 𝑙𝑜𝑔𝑝 ( x i 𝑜𝑏𝑠 , z k ; θ ) ] + 𝔼 q k [ 𝑙𝑜𝑔𝑞 k ( z k | x i 𝑜𝑏𝑠 ; ϕ ) ] - 𝔼 q k [ 𝑙𝑜𝑔 q k ( z k | x i 𝑜𝑏𝑠 ; ϕ ) p ( z k | x i 𝑜𝑏𝑠 ; θ ) ] ⩽ - 𝔼 q k [ 𝑙𝑜𝑔𝑝 ( x i 𝑜𝑏𝑠 , z k ; θ ) ] + 𝔼 q k [ 𝑙𝑜𝑔𝑞 k ( z k | x i 𝑜𝑏𝑠 ; ϕ ) ] = - 𝔼 q 0 [ 𝑙𝑜𝑔𝑝 ( x i 𝑜𝑏𝑠 , z k ; θ ) - ∑ n = 1 N log | 𝑑𝑒𝑡 ( ∂ ⁡ f n ⁢ ( z n - 1 ) ∂ ⁡ z n - 1 ) | + D 𝐾𝐿 [ q 0 ( z 0 | x i 𝑜𝑏𝑠 ; ϕ ) | | p ( z k ; θ ) ] ](11)

2) Online Imputation. Once the model is well-trained in offline model training, incomplete traffic raster sample x (e.g., Fig. 8) can be filled online. First of all, all the missing values x 𝑚𝑖𝑠 are initialized by zero. The imputation effectiveness is maintained because of the nonlinear transformation of neural network [31]. Then, a robust distribution of the latent variable is obtained through the encoder and SNFs, as introduced in Section 4. Lastly, according to the conditional distribution p θ ( x 𝑚𝑖𝑠 | x 𝑜𝑏𝑠 ) , missing value x 𝑚𝑖𝑠 can be imputed by sampling from a Markov chain [32] in decoder part.

p ( x 𝑚𝑖𝑠 | x 𝑜𝑏𝑠 ) = ∫ p θ ( x 𝑚𝑖𝑠 | Z K ) q ϕ ( Z K | x 𝑜𝑏𝑠 ) d Z K(12)

5.1 Datasets and experiment settings

1) Datasets. Experiments in this paper are conducted on three real-world traffic flow datasets: TaxiBJ [2], TaxiNYC [33] and BikeNYC [2], as shown in Table 1. TaxiBJ and TaxiNYC are GPS tracking data of taxis in Beijing and New York, BikeNYC contains trip data from New York bike rental system.

OPEN IN VIEWER

Missing values are manually generated with different missing rates at random. In this paper, traffic raster data within each day and hour is regarded as a sample for long-periodicity stream and short-periodicity stream respectively, considering the time interval of three datasets. Taking TaxiBJ as an example, if x t - 1 and x t are collected from the same hour, then one short-periodicity sample can be formulated as 𝑿 s = { x t - 1 , x t } , and the corresponding long-periodicity sample can be formulated as 𝑿 l = { x t - 25 , x t - 24 , … , x t - 1 , x t } . During data reconstruction phrase, Min-Max normalization is utilized to regulate data to the range of [ - 1 , 1 ] , because tanh (whose interval is also from - 1 to 1) is applied as the activation function of MSST-VAE in the final layer. All the datasets are divided training sets and test sets with the 9:1 scale. In case of over-fitting, 90% of the data in training sets are chosen for model training, and the rest is for validation. Model with the best performance on the validation set is saved, and the final results in Fig. 9 are the validation-set-saved model conducted on test sets.

2) Model Configuration. The hyper-parameters of MSST-VAE is set empirically as follows: Window size in long-periodicity stream is set to be 1 day, and 1 hour for short-periodicity stream. The dimension of each latent variable z k is fixed to 128. The length of Sylvester NFs is 5. The quantity of ECB and TCB is set to 4 and 3 respectively. The number of mid channels (output of 1D-dilated temporal convolution) for each ECB is 8, 44, 88, 176 in ascending order. The number of out channels (output of 2D-dilated spatial convolution) for each ECB is 16, 32, 64, 128 which is also in ascending order. In dilated P-3D convolution of ECB, the kernel sizes of each 2D-dilated spatial convolution are set to be 1 * 2 * 2, 1 * 3 *3, 1 * 6 * 6 and 1 * 7 * 7, and the stride is 1. Similarly, the kernel sizes of each 1D-dilated temporal convolution are set to be 2 * 1 * 1, 3 * 1 * 1, 6 * 1 * 1 and 7 * 1 * 1, and the stride is 1. The number of dilation factor of both two kinds of convolutions in each ECB is 8, 4, 2, 1 as in descending order. The number of out channels for each TCB is 64, 32, 16 and that of last 3D transposed convolution lay is 2. Kernel size for 3D transposed convolution in both TCB and last layer is set to 2, with padding of 1 in spatial and channel dimension, 0 in temporal dimension.

Besides, all the hyper-parameters and thresholds are set by our experience or at random initially, and hyper-parameter-searching is utilized to test all the assemblies to get the best performance hyper-parameter. The same operation is performed with other baseline approaches. All the hyper-parameters involved in some baseline methods of comparison experiments is shown as follows, and the hyper-parameters of the rest deep learning baseline methods are according to the configurations in the original papers.

The reduction ratio r of channel attention and temporal attention mechanism is set to be 2 and 8 respectively. The initial learning rate of Adam optimizer that employed in gradient descent is 10 - 3 . What’s more, the total number of epochs in short-term stream for TaxiBJ, TaxiNYC and BikeNYC are 700, 500 and 2000 separately. The total number of epochs in long-term stream for TaxiBJ, TaxiNYC and BikeNYC are 1200, 1000 and 3000 respectively. The total number of epochs in fusion model for TaxiBJ, TaxiNYC and BikeNYC are 100, 50 and 500 severally.

3) Evaluation Metrics. In order to compare our approach with other baseline methods, the matching of truth-ground data and imputed data are measured by Refined Normalized Mean Square Error (NMSE), Refined Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) expressed in Eqs (13), (14), and (15). NMSE is used to calculate absolute error of single-timestamp traffic raster data for each region in detail, while RMSE and MAE is for single-timestamp traffic raster data for all regions as a whole. Therefore, the former one can represents individual similarity and the latter two targets can be utilized to verify the global similarity. Moreover, RMSE is relatively more sensitive to ‘abnormal’ data when compared with MAE. So, high RMSE means textremely ‘bad’ padding value exists.

𝑅𝑀𝑆𝐸 = 1 T ⁢ ∑ t = 1 T ( 1 - m ) ⊗ ( x t - x t ^ ) 2(13) 𝑁𝑀𝑆𝐸 = 1 T ⁢ ∑ t = 1 T ∑ j = 1 J ∑ i = 1 I ( 1 - m ) ⊗ ( x i , j - x i , j ^ ) 2 ∑ j = 1 J ∑ i = 1 I ( 1 - m ) ⊗ ( x i , j ) 2(14) 𝑀𝐴𝐸 = 1 T ⁢ ∑ t = 1 T | ( 1 - m ) ⊗ ( x t - x t ^ ) |(15)

Where m is the mask tensor of observations defined in Section 3, ⊗ denotes Hadmard product, T denotes number of timestamp, I and J denotes the number of regions.

According to [13, 44, 45], to test the effectiveness and robustness of MSST-VAE as well as other comparing approaches in all kinds of situations, different missing ratios ranging from 10% to 50% with 10% interval are considered. Therefore, even when some sensors are in failure, data imputation can ensure the integrity of data and avoid system failure of ITS before the related sensors to be repaired.

In this section, a comparison is made among MSST-VAE and ten baseline methods: ARIMA [7], LLS [9], PPCA [13], HaLRTC [34], SDAE [35], MLP-VAE [5], CombCN [36], 3DConvGAN [46], DeepSTN+ [47], and STVAE [41], which can be classified into two groups. ARIMA, LLS, PPCA, and HaLRTC are traditional methods corresponding to the four sub-categories in Section 2. SDAE, CombCN, MLP-VAE, 3DConvGAN, DeepSTN+ and STVAE are deep learning-based methods. Moreover, MLP-VAE and STVAE are also VAE-based model for data imputation. Figures 9 and 10 show the comparison results of our approach with other baseline methods on three datasets at different missing rates.

OPEN IN VIEWER

ARIMA [7] directly predicts the missing value based on data from past timestamp at the same location. However, spatial dependency among different locations is not considered. LLS [9] imputes the missing value with the most “similar” value on the same raster map, whereas temporal dependency among data from different timestamps is not taken into consideration. Apparently, taking only temporal or spatial dependency into account is far from satisfactory in three measurements according to the results in Figs 9 and 10. It’s also worth noticing that the prediction ability of ARIMA degrades badly on account of the missing value according to the bad performance of MAE.

OPEN IN VIEWER

PPCA [13] tends to build a Gaussian latent variable model and then fill in the missing value according to the model. Yet, the relationship between latent and observed variables is linear, making it incapable of discovering complex spatio-temporal features. That is why PPCA performs worse at NMSE than MAE and RMSE in Figs 9 and 10, since global imputation ability on spatial-temporal dimension is limited. On the contrary, as a non-linear model, MSST-VAE outperforms PPCA on all datasets obviously.

HaLRTC [34] is a typical decomposition-based method mentioned in Section 2. It performs well at low missing rates on TaxiNYC and BikeNYC, for spatial dimension and temporal dimension are both considered. Nevertheless, it is not robust when missing rates rise, and even shows over-fitting when the missing rate reaches 90%. In contrast, regularization terms are added to the loss function Eq. (11) of MSST-VAE to prevent over-fitting. The results in Figs 9 and 10 also confirm that MSST-VAE performs relatively more stable under different missing rates. Moreover, as a decomposition-based method, the computation complexity of HaLRTC is higher than MSST-VAE (discussed in Section 3.2), which can also be verified by the results in Table 3.

Denosing Autoencoder (SDAE) [35] is a stacking deep learning method based on reconstruction architecture. Although encoder-decoder model is utilized in SDAE as MSST-VAE does, the hidden feature of SDAE being extracted is determinative. Contrast with VAE used in our model, which shows more generative ability, SDAE behave even worse than some traditional methods at RMSE, for there may exists some extremely bad imputed value. Besides, SDAE requires the completeness of data when training model, whereas MSST-VAE can be trained with incomplete traffic raster data, which is more suitable for practical engineering applications.

CombCN [36] is originally designed for video inpainting tasks, which applies a 3D convolution subnetwork for temporal consistency and a 2D convolution subnetwork for spatial modeling. In MSST-VAE, further improvements have been made, which is dilated P3D convolution. Compared with CombCN, global-local-granularity spatial features and long-short-term temporal features, which are the typical character of traffic raster data in ITS, can be extracted through MSST-VAE efficiently. The results of three measurements in Figs 9 and 10 also illustrate the superiority of dilated P3D gated convolution, comparing with traditional 2D and 3D convolution in CombCN.

MLP-VAE [5] is a deep Bayesian model that could learn stochastic mappings between observed data and latent variables. Due to the stochastic representations of VAE, the overall performance of MLP-VAE is second only to MSST-VAE, demonstrating the advantage of VAE with generalization ability in data imputation tasks. However, it is still not as good as MSST-VAE for two reasons: On the one hand, hidden variable in MLP-VAE obeys the Gaussian distribution, making the model less generalized, when compared with Sylvester NFs used in MSST-VAE. On the other hand, multi-layer perceptions are less satisfactory than dilated P3 convolution when extracting spatio-temporal features.

DeepSTN+ [47] is based on ST-ResNet to fulfil traffic data prediction task, in which different temporal periodicity streams are employed in the model (e.g day of week, hour of day), which is similar to the two VAE streams in MSST-VAE shown in Fig. 1. However, though the network of DeepSTN+ is deep, the operation within it is mainly based on plain 2D-Convolution, which is far from enough to learn the complex spatio-temporal features in traffic raster data. That’s why MSST-VAE outperforms it in all indexes according to Figs 9 and 10.

3DConvGAN [46] is an generative model based on GAN to deal with traffic imputation task, and 3D-Convolution is utilized to extract spatio-temporal correlations to improve the recovery accuracy. On the one hand, compared with VAE, GAN in 3DConGAN is hard to converge during training phrase since the input of generator in 3D-Convolution is randomly sampled from Gaussian distribution, which can be verified by the results in Table 3 (epoch of 3DConvGan can be sometimes 10 times more than ST-VAE in offline training phrase). On the other hand, dilated P3D convolution in MSST-VAE can learn short-long-term temporal features and global-local-granularity spatial features, while plain 3D-Convolution in 3DConvGAN can only gain simple spatio-temporal features. Also, 3D-Convolution is more time-consuming and resource-consuming than P3D-Convolution [27].

STVAE [41] is our previous work, compared with which, periodic features are considered by deploying two kinds of streams in MSST-VAE. What’s more, not only global-local-granularity spatial dependencies and long-short-term temporal dependencies are more emphasized but also multi-horizon attention mechanism is introduced in ECB, making MSST-VAE outperforms STVAE on three measurements according to Figs 9 and 10.

From the discussion above, we can draw a conclusion that MSST-VAE fulfills a more competitive imputation performance on traffic raster data than all the baseline models including the state-of-the-art VAE-based model ST-VAE, especially when the missing rate is high. Moreover, MSST-VAE shows steady performance on three indexes: Lowest NMSE and MAE demonstrating good local and global imputation ability, and lowest RMSE showing little extremely bad imputed value. What’s more, MSST-VAE displays strong robustness even at different missing rates.

In order to validate the effectiveness of major techniques in MSST-VAE, ablation experiments are conducted in this section. We reconfigure MSST-VAE to create four variants described as follows: 1) single-stream MSST-VAE. 2) MSST-VAE w/o SNFs. 3) MSST-VAE w/o dilated P3D convolution. 4) MSST-VAE w/o multi-horizon attention mechanism. Table 4 presents the NMSE and RMSE of MSST-VAE and four variants on three datasets at different missing rates.

In order to validate the effectiveness of major techniques in MSST-VAE, ablation experiments are conducted in this section. We reconfigure MSST-VAE to create four variants described as follows: 1) single-stream MSST-VAE. 2) MSST-VAE w/o SNFs. 3) MSST-VAE w/o dilated P3D convolution. 4) MSST-VAE w/o multi-horizon attention mechanism. Table 4 presents the NMSE and RMSE of MSST-VAE and four variants on three datasets at different missing rates.

Single-stream MSST-VAE. One of the input streams ( X L in Fig. 1) is removed from MSST-VAE, and only the short-periodicity data X s is fed into the network. It shows the second worst performance in comparison with the rest four approaches, demonstrating the importance of learning periodic feature for traffic raster data.

MSST-VAE w/o SNFs. SNFs in removed from MSST-VAE, and the hidden variables of VAE obey Gaussian distribution. Contrasting with MSST-VAE, when the generalization for hidden variable in VAE is weakened, the performance of imputation surely degrades according to Table 4.

MSST-VAE w/o dilated P3D convolution. Dilated P3D convolution in ECB and TCB is replaced with plain 3D convolution to extract spatio-temporal features. It behaves the worst of the four, because plain 3D extract global-local granularity spatial feature and long-short-term temporal feature less effectively and efficiently compared than dilated P3D convolution. In other words, the importance of spatio-temporal dependency for traffic raster data imputation task should be noted.

MSST-VAE w/o multi-horizon attention. Multi-horizon attention mechanism, which is used for better representation of spatio-temporal features, is removed from ECB in this part of experiment. Though the effectiveness is not as obvious as that in P3D, it is acceptable under such circumstance.

OPEN IN VIEWER

OPEN IN VIEWER

In order to evaluate the effectiveness of various parameterizations of MSST-VAE, hyper-parameter sensitive study is conducted on public datasets with average missing rates. Two major hyper-parameter is selected as the object of the study: number of ECB n and length of NFs l , which determinate the basic architecture of MSST-VAE. The results are as shown in Figs 11 and 12. It can be observed from the results that too many ECB and too deep NFs may not improve the performance of the MSST-VAE, therefore 4 is set for n and 5 is set for l .