Time Delay Estimation of Traffic Congestion Propagation Due to Accidents Based on Statistical Causality

Although our proposed method is the first to estimate time delays for congestion propagation in road traffic networks, time delay estimation (TDE) is not a new problem. In digital signal processing, TDE refers to the task of ascertaining the differences in arrival times between signals received at sensor array. The most widespread approach for TDE is cross-correlation [9, 10]. Supposing that signals are received from two sensors, the delay between the sensors can be estimated using the time lag that maximizes the cross-correlation between filtered versions of the received signal.

Limited attempts have been made to perform cross-correlation analyses with traffic speed data. Conventional TDE methods based on cross-correlation have been applied to a real road vehicle pass-by measurement to enable traffic monitoring using passive acoustic sensors [11]. Similarly, [12] conducted a cross-correlation analysis to prove the existence of a significant relationship between the current value of speed at a specific station, as well as past speed values at upstream and downstream stations in a freeway traffic network. We refer to this method as time-lagged cross-correlation (TLCC). This approach is not applicable in the presence of nonstationarity.

To quantify the TLCC level between two nonstationary time series at different scales, [13] proposed a time-lagged detrended cross-correlation analysis approach. This method, referred to as DCCA, divides an entire time series into overlapping boxes to handle nonstationarity [14].

The method proposed in the present study was compared with TLCC and DCCA for evaluation purposes.

Various traffic causal analysis methods have been developed to identify causal relationships among congested roads and detect congestion propagation patterns. The authors of [15] forecasted future traffic flow by ranking input variables to identify a subset of the Bayesian network as the set of cause nodes using the Pearson correlation coefficient. A two-step mining architecture has been proposed to capture the origin points of road traffic congestion [16]. TE was used to reveal the delay propagation network among multiple airports with a time series of airport delays [17]. However, the aforementioned approaches primarily focus on revealing causal relationships and do not provide information regarding estimated time delays during congestion propagation.

Suppose that $\{X_{t}\}_{t\geq 1}$ is a stationary Markov chain with a finite state space $S=\{s_{1},\ldots,s_{n}\}$, where $n\in\mathbb{N}$. Let $\mathbf{P}=(p_{ij})\in\mathbb{R}^{n\times n}$ be the transition probability matrix of the chain and the stationary distribution by $\boldsymbol{\pi}=(\pi_{1},\dots,\pi_{n})$. Thus, for any $1\leq i,j\leq n$, $p_{ij}=P(X_{t+1}=s_{j}|X_{t}=s_{i})$ and $\pi_{i}=P(X_{t}=s_{i})$. Given a time series $\{X_{1},\ldots,X_{L}\}$ of size $L$ from a stationary Markov chain, $\pi_{i}$ and $p_{ij}$ can be estimated as

The bootstrap observations $\left\{X_{1}^{*},\ldots,X_{L}^{*}\right\}$ can now be generated using the estimated transition matrix and marginal distribution in Eq. (3.1) [18].

1.      1)

   Generate a random variable $X_{1}^{*}$ from the discrete distribution on $\{1,\ldots,n\}$ that assigns mass $\hat{\pi}_{i}$ to $s_{i}$, $1\leq i\leq n$.

2.      2)

   Generate a random variable $X_{t+1}^{*}$ from the discrete distribution on $\{1,\ldots,n\}$ that assigns mass $\hat{p}_{ij}$ to $j$, $1\leq j\leq n$, where $s_{i}$ is the value of $X_{t}^{*}$.

3.      3)

   Repeat Step 2) until a simulated time series $\{X_{1}^{*},\ldots,X_{L}^{*}\}$ has been obtained.

TE is a measurement of directed information flow [19] based on the concept of Shannon entropy [20] in the field of information theory. For a discrete random variable $I$ with probability distribution $p(i)$, the Shannon entropy represents the average number of bits required to optimally encode independent draws, calculated as follows:

Eq. (3.2) can be easily extended to the concept of conditional entropy using two discrete random variables $I$ and $J$:

This equation can be used to measure information flow between two discrete random variables.

Consider two discrete random variables, $I$ and $J$, with marginal probability distributions $p(i)$ and $p(j)$, and joint probability $p(i,j)$. Suppose both variables represent stationary Markov processes of orders $k$ and $l$, respectively. For the order $k$ Markov process $I$, Eq. (3.2) can be extended to

where $i_{t-1}^{(k)}=(i_{t-1},\ldots,i_{t-k})$. Analogously, the information flow from $J$ to $I$ is measured by quantifying the deviation from the generalized Markov property $p(i_{t}|i_{t-1}^{(k)})=p(i_{t}|i_{t-1}^{(k)},j_{t-u}^{(l)})$ for an arbitrary source–target lag $u$, as follows:

Eq. (3.4) preserves the computational interpretation of TE as an information transfer, which is the only relevant option in keeping with Wiener’s principle of causality [21]. The transfer entropy is known to be biased in small samples [22]. To correct any bias, [22] proposed the effective transfer entropy ($ETE$):

where $T_{J_{\text{shuffled}}\rightarrow I}^{(k,l)}$ indicates the transfer entropy using a shuffled version of time series $J$. The shuffling process randomly draws values from the original time series and realigns them to generate a new time series. Thus, shuffling eliminates any time series dependencies of $J$, as well as statistical dependencies between $J$ and $I$. Note that $T_{J_{\text{shuffled}}\rightarrow I}^{(k,l)}$ converges to zero as the sample size increases, and any nonzero value of $T_{J_{\text{shuffled}}\rightarrow I}^{(k,l)}(t,u)$ is a result of the small sample effect. To ensure estimation consistency, shuffling is repeated, and the average of the shuffled transfer entropy estimates across all iterations serves as an estimator for the small sample bias, which is subsequently subtracted from the Shannon or Rényi transfer entropy estimate to correct any bias.

Suppose that congestion information is transferred from a source road to a destination road with a time delay of $u$. The objective of time delay estimation is to estimate $u$ given previously observed traffic speed data on the two roads, denoted as $\{X_{t}\}^{L}_{t=1}$ and $\{Y_{t}\}^{L}_{t=1}$, respectively. Therefore, the time delay estimation task entails formulating a function $f(\cdot)$ that computes the source–target lag $u$, $\left[\{X_{t}\}_{t=1}^{L},\{Y_{t}\}_{t=1}^{L}\right]\xrightarrow[]{f(\cdot)}u.$

The proposed time delay estimation algorithm comprises three steps. First, bootstrapping for each time series is performed using preprocessing methods. Next, the transfer entropy computation provides the estimated time delay lag $u^{*}$. Finally, the distribution of time delay lag is estimated to determine the existence of a statistical causal relationship.

Consider a time series of congested traffic speed data, $\{X_{t}\}^{L}_{t=1}$, which has the properties of scale dependence, nonlinearity, and nonstationarity. To effectively identify the causal relationship within such a time series, appropriate preprocessing methods are essential.

To this end, we first decompose the time series into a trend and its residual, as follows:

where $\mathcal{T}_{t}$ and $R_{t}$ are the trend and residual components, respectively. The trend component is a moving average of order $m$, representing the mean forefront value at time $t$. The purpose of the trend component is to smooth the time series for estimating the underlying trend. After extracting the underlying trend from $\{X_{t}\}^{L}_{t=1}$, the residual time series $\{R_{t}\}^{L}_{t=1}$ is assumed to be a stationary Markov process. The assumption of Markovian property in traffic speed is not a novel notion. Many prior traffic speed prediction and modeling studies have been conducted under this assumption [23, 24, 25, 26, 27]. Based on $\{R_{t}\}^{L}_{t=1}$, we can generate the bootstrap residuals $\{R^{*(b)}_{t}\}^{L}_{t=1}$ as explained in Section 3.1. Subsequently, we can easily obtain a bootstrap time series $\{X^{*(b)}_{t}\}^{L}_{t=1}$ by $\mathcal{T}_{t}+R^{*(b)}_{t}$ for $t=1,\ldots,L$.

Nonlinear normalization with a sliding window is then applied to the obtained time series to address the scale-dependency, nonlinearity, and nonstationarity of traffic speed data. To ensure the data are scale-independent and close to linear [28], the standard normal cumulative distribution function $\Phi$ is applied. A sliding window technique has similarly been employed to handle a nonstationary time series in [29]. Let $\mathbf{X}_{t,w}=\left\{X^{*(b)}_{k}\right\}_{k=t-w+1}^{t}$ be the forefront sequence of $X^{*(b)}_{t}$ with a sliding window size of $w$, and $F_{25,t}$, $F_{50,t}$, and $F_{75,t}$ be the 25th, 50th, and 75th percentiles of $\mathbf{X}_{t,w}$, respectively. Note that these percentiles depend on the location of the sliding window. Then, a normalized time series $\left\{\tilde{X}^{*(b)}_{t}\right\}^{L}_{t=1}$ can be obtained by

To verify the effectiveness of the nonlinear normalization method expressed in Eq. (4.2), we compared its performance with that of existing normalization methods [28, 29], including the min-max method $\left(\tilde{X}^{*(b)}_{t}=\frac{X^{*(b)}_{t}}{\max\mathbf{X}_{t,w}}\right)$ and the z-score method $\left(\tilde{X}^{*(b)}_{t}=\frac{X^{*(b)}_{t}-\mu(\mathbf{X}_{t,w})}{\sigma(\mathbf{X}_{t,w})}\right)$ with and without a sliding window technique (see Section 5).

Using Eq. (3.5), the time lag in a causal relationship $J\rightarrow I$ can be estimated by solving the following optimization problem:

In this study, we assume $k=\ell=1$.

To compute the lag-specific ETE in Eq. (4.3), we discretize continuous data using symbolic encoding. This discretization can be performed by partitioning the data into a finite number of bins. We denote the bounds specified for the $n$ bins by $q_{1},q_{2},\ldots,q_{n-1}$, where $q_{1}<q_{2}<\cdots<q_{n-1}$. For the normalized time series in Eq. (4.2), we obtain the encoded time series $\{J^{*(b)}_{t}\}^{L}_{t=1}$ by the following equation:

The choice of bins depends on the distribution of data. In the case where tail observations are of particular interest, binning is typically based on empirical quantiles, such that the left and right tail observations are allocated into separate bins. In this study, we implemented symbolic encoding with $n=3$ based on 5% and 95% empirical quantiles, thereby emphasizing speed extremes caused by dynamic speed changes and traffic accidents.

Consequently, we obtain $\{J^{*(b)}_{t}\}^{L}_{t=1}$ from $\{\tilde{X}^{*(b)}_{t}\}^{L}_{t=1}$, and $\{I^{*(b)}_{t}\}^{L}_{t=1}$ from $\{\tilde{Y}^{*(b)}_{t}\}^{L}_{t=1}$, respectively, for $b=1,\ldots,B$. Given $\{J^{*(b)}_{t}\}^{L}_{t=1}$ and $\{I^{*(b)}_{t}\}^{L}_{t=1}$, $b=1,\ldots,B$, we obtain bootstrap observations of the time lag, $u^{*(1)},\ldots,u^{*(B)}$ using Eq. (4.3).

Suppose that bootstrap observations of the time lag follow a distribution $\mathcal{G}$,

which is unknown in practice. Let $\mu$ and $\sigma^{2}$ denote the mean and variance of $\mathcal{G}$, respectively, which can be estimated by

Proposition 1 implies that 1) the bootstrap estimate $\hat{\mu}_{B}$ is an unbiased estimate of $\mu$, and 2) $\frac{1}{B}\hat{\sigma}_{B}^{2}$ quantifies the uncertainty of $\hat{\mu}_{B}$. That is, we can evaluate the uncertainty of the bootstrap estimate $\hat{\mu}_{B}$ using $\frac{1}{B}\hat{\sigma}_{B}^{2}$, which is practically useful because $\mu$ is unknown. This approach can be applied to hyperparameter tuning. In this study, we used a grid search to determine a set of hyperparameters (length of time series ($L$) and sliding window size ($w$)) that minimizes $\frac{1}{B}\hat{\sigma}_{B}^{2}$.

To determine whether the bootstrap estimate $\hat{\mu}_{B}$ is reliable, we compare $\hat{\sigma}_{B}^{2}$ with a predetermined threshold $\sigma_{T}^{2}$. That is, if $\hat{\sigma}_{B}^{2}>\sigma_{T}^{2}$, we can conclude that $\hat{\mu}_{B}$ is not reliable, and congestion is therefore not propagated on the corresponding road segment.

To determine an appropriate value of $\sigma_{T}^{2}$, we employ the concept of the tolerance interval ($TI$). The $TI$ is a statistical interval in which a specified proportion $\gamma$ of a population will fall with a certain level of confidence $(1-\alpha)$. By definition, a $(\gamma,1-\alpha)$-$TI$ of $\hat{\mu}_{B}$,

satisfies

where $k_{\gamma,\alpha}$ is the tolerance factor [31]. Then, $\sigma_{T}^{2}$ can be determined by $k_{\gamma,\alpha}\sqrt{\frac{1}{B}\sigma_{T}^{2}}=1$ (min) to yield a $\pm 1$ minute $TI$. With $B=100$, $\gamma=0.9$, and $\alpha=0.01$, the present study uses $\sigma_{T}^{2}=\frac{100}{k^{2}_{0.9,0.01}}=5.05^{2}$.

Consider the propagation of traffic congestion caused by accidents. Here, we define the propagation path by a sequence of incoming roads in the direction opposite to the traffic flow, where $k$th element in the propagation path is denoted as Hop($k$-1). Let $\hat{\mu}_{B,k-1}$ and $\hat{\sigma}^{2}_{B,k-1}$ denote the bootstrap estimate and sample variance at Hop($k$-1), respectively. Hop$0$ corresponds to the road where the accident occurred. We consider Hop($k$-1) to be statistically significant if $\hat{\mu}_{B,k-1}$ and $\hat{\sigma}^{2}_{B,k-1}$ satisfy the following conditions: 1) $\hat{\sigma}^{2}_{B,k-1}<\sigma^{2}_{T}$ and 2) $\hat{\mu}_{B,k-2}<\hat{\mu}_{B,k-1}$. The second condition states that the time delay between Hop$0$ and Hop($k$-1) must exceed that between Hop$0$ and Hop($k$-2).

The accident occurred on September 8, 2020, at 06:44 AM. The blue star in Figure 4 denotes the exact location of the accident. Case 1 has one propagation path $[A,B,C,D]$. The black, red, blue, and green line segments in the figures indicate Hop$0$, Hop$1$, Hop$2$, and Hop$3$, respectively. Time delays were estimated using average speed data recorded at one-minute intervals from the previous hour to the subsequent two hours based on the time when the accident was reported.

Figures 5(a) and 5(b) show the results of time delay estimation for the propagation path $[A,B,C,D]$. It is apparent from the significantly smaller values of $\hat{\sigma}_{B}^{2}$ that the time series with nonlinear normalization produced a more consistent estimate of time delays that increased with each hop. We, therefore, conclude that the congestion effect of the accident propagated along the path to Hop1, Hop2, and Hop3 at 3.60, 7.30, and 19.97 min after the accident, respectively.

Unlike the proposed method, both TLLC and DCCA failed to identify the congestion propagation effect of the accident. Furthermore, the DCCA method produced inconsistent time delay estimates for varying values of $n$. Note that both the TLLC and DCCA methods do not provide uncertainty quantification of the time delay estimates.

The accident occurred on September 4, 2020 at 10:16 PM, and affected five propagation paths: $[A,B,C,D]$, $[A,E,F,G]$, $[A,H,I,J]$, $[A,H,K,M]$, and $[A,H,K,L]$. Figure 6 depicts the exact location of the accident, along with the five propagation paths. For each path, the previous hour and the subsequent two hours were considered for time delay estimation. Figure 7 and Table 6 present the results of time delay estimation. In Path 1, no specific causal relationship could be found, as shown in Figure 7(a). This finding is supported by the corresponding high values of $\hat{\sigma}_{B}^{2}$ ($>\sigma^{2}_{T}=5.05^{2}$) in Table 6. Similarly, we can conclude that the congestion effect of the accident at road $A$ propagated along the second hop of Paths 2 and 3, and the third hop of Paths 4 and 5. Moreover, the values of $\hat{\sigma}_{B}^{2}$ in Table 6 indicate that the congestion effect propagated along Path 3 up to Hop2, as depicted in Figure 7(b), and along Paths 4 and 5 up to Hop3, as depicted in Figure 7(c). For Paths 4 and 5, the time delay estimates are (8.23, 15.65, 22.06) and (8.22, 15.55, 20.75), respectively.

As in Case 1, both the TLCC and DCCA methods failed to estimate consistent time delays. The DCCA method with 30 overlapping boxes ($n=30$) obtained comparable results with the proposed method only for Path 3, as seen in Table 6. From the results of both cases, it can be concluded that the proposed method identifies causal relationships and estimates the time lag more accurately than the conventional TDE methods.