1 Introduction↩︎

Driving safely towards a goal destination is a key functionality in automated vehicle (AV) design. However, assessing safety is an intricate matter, as even a motion plan that is collision free may not feel safe for the passengers of the AV. Furthermore, determining whether a (future) motion plan is collision free requires precise knowledge of the current and future states of obstacles and other road participants. Such knowledge, however, is generally only partial in nature as measurement and estimation introduces uncertainties. To account for uncertainties, risk functions are commonly used to assess the safety of a motion plan, where risk is understood as stochastic quantity describing potential harm given limited knowledge, i.e., uncertainty, about current and future driving states.
Risk functions based on this rough definition are typically constituted of two aspects: one related to the uncertainty on current and future driving states and the other to the severity of an interaction, e.g., a collision. Such risk measures have been shown to correlate with perceived risk by humans [1]–[3], where [1] has also shown that constraining risk within a motion planning algorithm in a static environment (i.e., without other dynamic actors) led to human-comparable vehicle behavior. Furthermore, minimizing or constraining risk within a motion planner is established to generate safe AV behavior (see, e.g., [4]–[8]).
A key challenge is computationally efficient risk estimation, which is generally desired in real-time motion planning applications. As analytic methods are scarce, some works (see, e.g., [2], [8]–[11]) resort to modelling risk by empirically parameterizing model equations, thus not guaranteeing that the uncertainty is correctly characterized in a probabilistic sense, and generalizability beyond the scope of data used for tuning the model is unclear. Other works choose to over-approximate the risk by representing each actor by a circular shape (see, e.g., [4], [5]), which simplifies the calculations significantly; however, it clearly undermines the accuracy of estimation. Furthermore, in [5], [12], [13], the risk is calculated with Monte Carlo sampling, which itself is causing a significant computational burden. An important observation is that if the severity term is constant, one is left to compute the probability of collision. Recently, [14] proposes a new method based on a multi-circular shape approximation of each actor to estimate the collision probability. That method is particularly tailored to optimization-based motion planning as it is guaranteed to over-approximate the collision probability value while being computationally efficient and supporting the solver in generating smooth, reproducible trajectories. From a risk estimation perspective, modeling a vehicle by individual areas, e.g., multiple overlapping circles, is a promising approach not only for uncertainty estimation but also for the calculation of the severity, since the collision severity is known to depend on the location of impact for each actor, i.e., the collision constellation [15]–[17].
This paper builds upon the multi-circular approach for collision probability estimation of [14] to develop a novel and accurate methodology for risk estimation. This approach, in which all road actors are represented by multiple overlapping circles, allows to assign different severity functions for each pair of colliding circles. Here, a particular circle may be related to the front, middle, rear, etc., of a vehicle, and, hence, e.g., a front-to-side collision can be represent by assigning individual severity functions for such a circle-to-circle pair. To combine many individual severity functions, we introduce an averaging method to combine individual severity functions to accurately represent different collision constellations by a single severity value. Our approach retains the over-approximate guaranties in the uncertainty aspect and also inherits the advantageous properties for optimization-based motion planning from [14], by using the same algorithmic structure and multi-circular approach. Finally, we present an exemplary implementation of the proposed risk estimation approach for Gaussian uncertainties and a kinetic energy model to assign collision severity. The resulting algorithm is computing sufficiently fast to support real-time motion planning applications and is provided in open-source.
The remainder of this article is structured as follows. Section 2 introduces preliminary notation and states the problem of risk estimation. Following, Section 3 derives the risk based on a multi-cirucular shape approximation that over-approximates the uncertainty aspects and includes a model to combine individual severity functions to account for different collision constellations. Section 4 presents an implementation for Gaussian uncertainties with a kinematic-energy severity model, for which we then demonstrate the risk estimation in Section 5. Lastly, we conclude this article in Section 6.

2 Problem Statement↩︎

Consider an arbitrary traffic scene with two actors: an automated vehicle, called the ego, and another road user (e.g., another vehicle or a cyclist), called the object, see Figure 1. For clarity of exposition, we restrict ourselves to only one object; however, our method also applies to any number of objects. From the perspective of the ego, the object is characterized by its configuration \(\mathbf{y} := (\mathbf{x}, \theta)\in \mathcal{C} := \mathbb{R}^3\), which is composed of the Cartesian position of the object’s geometrical center with respect to the ego’s geometrical center, \(\mathbf{x} := (x, y) \in\mathbb{R}^2\), and the object’s heading angle measured with respect to the ego’s longitudinal axis, \(\theta \in \mathbb{R}\) (\(\mathcal{C}\) is called the configuration space). Note that we define \(\theta\) on \(\mathbb{R}\), as for later calculations it is beneficial to represent angles as multiples of \(2\pi\). Figure 1 depicts all the described coordinates. In the sequel, a configuration at the discrete-time instant \(k\) is denoted by \(\mathbf{y}_{k}\).
We define risk as the expected severity of a risk-inducing event \(\mathcal{E}\) (e.g., collisions or violating traffic regulation), as in [5]. In this paper, we only consider collisions and, hence, our objective is to estimate collision risk. To estimate the expected severity of a collision (i.e., collision risk), one needs to derive the conditions on the relative position and heading angle of the vehicles that, when satisfied, imply the vehicles are colliding and, if so, what the severity of that collision is.
Let \(\mathbf{z}_k \in \mathbb{R}^{n_\mathbf{z}}\) denote the variables involved in determining the occurrence and severity of a collision event \(\mathcal{E}\) at time \(k\), these are, e.g., positions, orientations, and velocities of both actors. In general \(\mathbf{z}_k\) is a vector of random variables with associated probability density function (PDF) \(p_{\mathbf{z}}(\cdot ;\mathbf{w}_k)\). Here, we assume that the random variables are independent and that \(p_{\mathbf{z}}\) has the same functional form for all \(k\), but time-varying parameters denoted by \(\mathbf{w}_k\) (e.g., time varying mean and variances), which characterizes the time dependence of the PDF. Let the event set \(\mathcal{B}_{\mathcal{E}} \subset \mathbb{R}^{n_\mathbf{z}}\) be such that if \(\mathbf{z}_k\in \mathcal{B}_{\mathcal{E}}\), then the event \(\mathcal{E}\), i.e., a collision, has occurred. Note that in order to determine whether a collision has occurred, only information about the configuration \(\mathbf{y}\) is required, as \(\mathbf{y}\) suffices to check if both actors’ footprints intersect. However, to estimate the severity, additional information about the kinematic variables of both actors, e.g., the velocities \(v_e, v_o\), see Figure 1, may also be needed.
Based on accidentology research, see e.g., [15]–[17], we assume that the severity of a collision depends on the collision constellation, i.e., the specific sections of each vehicle which are colliding. Therefore, consider that the ego vehicle consists of \(N_e\) sections (e.g., \(N_e = 3\): front, center, rear) and the object is constituted of \(N_o\) sections. To each specific section pair (e.g., ego front section to object rear section), we allocate an individual severity function that assigns the severity for that section pair, based on the actors’ velocities \(v_e, v_o\) and their relative directionality \(\theta\), see the collisions among different vehicular sections in Figure 1. Hence, to the \(j\)-th ego section, \(j\in \{1,2, ...N_e\}\), colliding with the \(l\)-th object section, \(l\in \{1,2,...,N_o\}\), we assign the severity with \(s_{j,l}\) such that \((\theta, v_e, v_o) \mapsto s_{j,l}(\theta, v_e, v_o) \in \mathbb{R}_{\geq0}\). Moreover, we collect all individual severity functions in a set \(\mathcal{S} = \{ s_{j,l} \mid \forall j \in \{1, 2, \dots, N_e\}, \forall l\in \{1, 2, \dots, N_o\}\}\). With the collection of all severity functions in \(\mathcal{S}\), we define a function \(f_s\) that combines all \(N_eN_o\) individual severity functions into a total severity function. In practice, this function may, e.g., represent the sum of all individual severities or the maximum severity. To retain generality, we define the total severity as \[\label{eq:total95severity} s = f_s(\mathcal{S}), \text{ such that } s \geq 0.\tag{1}\] As a result, an element \(\mathbf{z} \in \mathcal{B}_{\mathcal{E}}\) contains \(x, y, \theta, v_e, v_o\), of which \(x, y, \theta\) are needed for determination of collision probability and \(\theta, v_e, v_o\) are used to assess the severity. Note that although \(\theta\) is needed for both, its range for which the event \(\mathcal{E}\) occurs is restricted by the set \(\mathcal{A}_{rec} \in \mathcal{C}\) of configurations that lead to a collision among two vehicles’ rectangular footprints. Finally, we assume that all kinematic variables associated with the ego vehicle are known without any uncertainty, thus, \(v_e\) will be treated as a parameter to our risk estimation. Hence, \(\mathcal{B}_\mathcal{E} := \mathcal{A}_{rec} \times \mathcal{V}_o\), where \(\mathcal{V}_o\) is a subset of \(\mathbb{R}\) and represents the possible collision velocities of the object. Therefore, the risk \(R_k\) associated with \(\mathcal{E}\) at time \(k\) is defined as \[\label{eq:risk95general951} R^{rec}(v_e,\mathbf{w}_k) := \int_{\mathcal{B}_{\mathcal{E}}}s(\cmbtheta, v_e, \cmbv_{o})p_{\mathbf{z}}(\cmbbz;\mathbf{w}_k)\mathrm{d}\cmbbz.\tag{2}\] Note that we distinguish integration variables with the computer modern bright font, hence, we use \(\cmbtheta, \cmbv_{o}, \cmbbz\) instead of \(\theta,v_o, \mathbf{z}\) in (3 ), as we integrate over all random variables in \(\mathbf{z}\). The objective of this paper is to estimate \(R_k\) in (2 ) by over-approximating both actors’ shapes with multiple overlapping circles, giving a larger set of possible collision configurations, i.e., \(\mathcal{A}_{cir} \supset \mathcal{A}_{rec}\). The over-approximation is necessary, as the multi-circular approach allows for a more efficient computation than rectangular approaches, and under-approximation is generally undesired in safety-critical applications. Here, the multi-circular over-approximation naturally yields the individual severity functions \(s_{j,l}\) by considering individual ego and object circle pairs, e.g., a front object circle colliding with a rear ego circle may represent a rear end collision. The individual severity functions \(s_{j,l}\) are then combined by \(f_s\) in (1 ). Thus, the domain of integration for the multi-circular shape representation is given by \(\tilde{\mathcal{B}}_\mathcal{E} := \mathcal{A}_{cir} \times \mathcal{V}_o\). Hence, our objective is to estimate \[\label{eq:risk95general} R^{cir}(v_e; \mathbf{w}_k) := \int_{\tilde{\mathcal{B}}_{\mathcal{E}}}s(\cmbtheta, v_e, \cmbv_{o})p_{\mathbf{z}}(\cmbbz;\mathbf{w}_k)\mathrm{d}{\cmbbz}.\tag{3}\]

3 Methodology for Risk Estimation↩︎

In order to determine the risk at a certain time instant \(k\), the integral in (3 ) needs to be evaluated. In general, analytic solutions are not available, hence, the risk is either over-approximated [7] or estimated with random sampling, i.e., Monte Carlo methods [5], [12]. In the sequel, we will construct the risk for an averaged sum of severities in (1 ), which is averaged in cases a collision among multiple circles occurs. Since the risk is estimated in the same way for all times \(k\), we omit time indexing in the sequel for notational clarity. As we assume independence of the random variables in \(\mathbf{z}\), we can decompose the PDF into its components as \(p_{\mathbf{z}}(\mathbf{z};\mathbf{w}) = p_{\mathbf{y}}(\mathbf{y}; \mathbf{w}_{\mathbf{y}})p_{v_o}(v_o;\mathbf{w}_{v_o}) = p_{\mathbf{x}}(x, y, \mathbf{w}_{\mathbf{x}})p_{\theta}(\mathtt{\theta}; \mathbf{w}_{\theta})p_{v_o}(v_o;\mathbf{w}_{v_o})\).

3.1 From Collision Probability to Risk of Collision↩︎

In fact, computing the risk in (3 ) is closely related to calculating the probability of collision between a vehicle and an object. Considering the relative object configuration \(\mathbf{y}\) to be uncertain and \(\mathcal{A}_{rec}\) to be the set of all relative object configurations leading to a collision with the ego, the probability of collision is given by \[\label{eq:general95POC} \mathbb{P}\{\mathbf{y} \in\mathcal{A}_{rec}\} = \int\displaylimits_{ \mathcal{A}_{rec}} p_{\mathbf{y}}(\cmbby; \mathbf{w}_{\mathbf{y}}) \text{d} \cmbby.\tag{4}\] In [14], (4 ) is over-approximated by enclosing both the vehicle’s and object’s rectangular footprints in multiple overlapping circles. This approach allows for fast computation of the collision probability, while guaranteeing over-approximation. We will extend that approach from collision probability estimation to risk estimation.

Suppose that the vehicle and object shapes are (tightly) covered by \(N_{e}\) and \(N_{o}\) overlapping circles with radii \(r_e, r_o\), respectively (see Figure 2 (a)). The circles of each actor are placed equidistantly apart by distances \(d_{e}\) and \(d_{o}\) between the circle centers for the ego and object, respectively. Then, the furthest distance between both actors’ geometric centers, for which a collision of the multi-circular shape approximations can occur, is given by \[\label{eq:radial95bound} \overline{\rho} = r_e + r_o + \frac{d_{o}}{2}\big(N_{o} - 1\big) + \frac{d_{e}}{2}\big(N_{e} -1\big),\tag{5}\] see [14]. We use a polar coordinate frame because that allows us to restrict the maximum collision distance (5 ) solely in the radial coordinate, hence, we define a change of coordinates as \[\label{eq:polar95transform} \begin{align} & CT: \mathbb{R} \times [0, 2\pi) \rightarrow \mathbb{R}^2,\\ & (\rho, \phi) \mapsto (x, y) = (\rho \cos{\phi}, \rho \sin{\phi}). \end{align}\tag{6}\] Then, for all angles \(\phi \in [0, 2\pi)\), collisions are possible for \(\rho \in [0, \overline{\rho}]\). We denote \(\mathbf{y}_{p} := (\rho, \phi, \theta)\) a polar representation of the configuration, for which it holds that \(\mathbf{y}_{p}= (\sqrt{x^2 + y^2}, \text{atan2}(y, x), \theta)\), where ‘\(\text{atan2}\)’ designates the two-argument arctangent function. Recall that \(p_{\mathbf{y}}(x, y, \theta; \mathbf{w}_{\mathbf{y}}) = p_{\mathbf{x}}(x, y; \mathbf{w}_{\mathbf{x}})p_{\theta}(\theta; \mathbf{w}_\theta)\), where we apply the same change of coordinates to \(p_{\mathbf{x}}\), giving \(p_{\phi, \rho}\), and, hence, \(p_{\mathbf{y}}(x, y, \theta; \mathbf{w}_{\mathbf{y}}) = p_{\phi, \rho}(\phi, \rho; \mathbf{w}_{\phi, \rho})p_{\theta}(\theta; \mathbf{w}_\theta)\) with the parameterization \(\mathbf{w}_{\phi, \rho}\) for the position.
We denote by \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho) = [\underline{\theta}_{j,l}(\phi, \rho), \overline{\theta}_{j,l}(\phi, \rho)]\) an intersection angle interval, i.e., an interval for the heading angle \(\theta\) for which the object circle \(l\) intersects with the ego circle \(j\) for a specific position \((\phi, \rho) \in [0, 2\pi] \times [0, \overline{\rho}]\). For the detailed derivation of expressions for \(\mathbb{I}_{j,l}^{\theta}(\phi, \rho)\), we refer to [14]. With all \(N_e N_o\) intersection angle intervals, for a specific \((\phi, \rho) \in [0, 2\pi) \times [0, \overline{\rho}]\), the region of integration of the PDF for the object’s heading angle \(\theta\) is given by \[\label{eq:collision95set} \breve{\mathbb{I}}^{\theta}(\phi, \rho) := \bigcup_{j = 1}^{N_e}\bigcup_{l = 1}^{N_o} \mathbb{I}^{\theta}_{j,l}(\phi, \rho) = \bigcup_{i=1}^{N_I}\breve{\mathbb{I}}^{\theta}_{i}(\phi, \rho).\tag{7}\] Note that, in (7 ), several intervals intervals \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\) may overlap, which is illustrated in Example 1 below. The right-most expression in (7 ) therefore denotes that there might be up to \(N_I\) disjoint intervals \(\breve{\mathbb{I}}^{\theta}_{i}(\phi, \rho)\). Thus, the set of all polar object configurations leading to a collision with the ego vehicle is given by \[\label{eq:Acir} \mathcal{A}_{cir, p} = \{(\phi, \rho, \theta) \in \mathcal{C}_p\! \mid \! \phi \in [0, 2\pi), \rho \in [0, \overline{\rho}], \theta \in \breve{\mathbb{I}}^{\theta}(\phi, \rho) \},\tag{8}\] where \(\mathcal{C}_p\) is the polar transformation of the configuration space, i.e., \(\mathcal{C}_p := CT^{-1}(\mathbb{R}^2) \times \mathbb{R}\). Note that the same transformation can also be applied to \(\mathcal{A}_{rec}\), giving \(\mathcal{A}_{rec,p}\). With these considerations, (4 ) can be over-approximated, i.e., \(\mathcal{A}_{cir, p} \supset \mathcal{A}_{rec,p}\), as \[\label{eq:poc95overapproximation} \begin{align} &\mathbb{P}\{\mathbf{y}_{p} \in \mathcal{A}_{cir, p}\} = \\ &\int_0^{2\pi}\int_0^{\overline{\rho}} p_{\phi, \rho}(\cmbphi, \cmbrho; \mathbf{w}_{\phi, \rho} ) \left[ \sum_{i=1}^{N_{I}} \int_{\breve{\mathbb{I}}^{\theta}_{i}(\cmbphi, \cmbrho)} p_{\theta}(\cmbtheta; \mathbf{w}_{\theta})\text{d}\cmbtheta \right]\text{d}\cmbrho \text{d}\cmbphi. \end{align}\tag{9}\]

The next step is to extend the probability of collision toward a risk measure by also considering the collision severity, see (3 ). Note that the circles of one vehicle can overlap (see the circular placement per vehicle in Figure 2). Hence, there is one or more areas (depending on the respective number of circles \(N_e\) and \(N_o\)) with two severity functions, essentially just creating another collision area which the function \(f_s\) must account for. Combining all severity functions \(s_{j,l}\) within (1 ) and using the total severity function and the characterization for the probability of collision (9 ), the risk in (3 ) can be estimated (in fact, over-approximated) with \[\label{eq:ego95risk95full} \begin{align} &R^{cir}(v_e; \mathbf{w})= \int_0^{2\pi}\!\int_0^{\overline{\rho}}\!p_{\phi, \rho}(\cmbphi, \cmbrho; \mathbf{w}_{\phi, \rho})\Bigg[\sum_{i=1}^{N_{I}}\\ & \int_{\breve{\mathbb{I}}^{\theta}_{i}(\cmbphi, \cmbrho)} \int_{\mathcal{V}_o}\!p_{\theta}(\cmbtheta; \mathbf{w}_{\theta})p_{v_o}(\cmbv_o; \mathbf{w}_{v_o})s(\cmbtheta, v_e, \cmbv_o)\text{d}\cmbv_o \text{d} \cmbtheta\Bigg]\text{d}\cmbrho \text{d}\cmbphi, \end{align}\tag{10}\] where \(\mathbf{w} = (\mathbf{w}_{\phi, \rho}, \mathbf{w}_{\theta}, \mathbf{w}_{v_o})\).

3.2 Combining Severity Functions↩︎

In (10 ), the total severity function \(s\) is integrated over \(\theta\) for \(N_I\) disjoint intersection angle intervals \(\breve{\mathbb{I}}^{\theta}_{i}(\phi, \rho)\). However, our individual severity functions \(s_{j,l}\) are defined on a circle-to-circle basis. Due to the problem of overlap (see Example 1), integrating \(\theta\) over all \(N_eN_o\) intervals \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\) would result in a much larger over-approximation as the areas within an overlap are accounted for multiple times. In (9 ), this is resolved by taking the union over all intersection angle intervals \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\) in (7 ), since the information about which circle-to-circle pair is colliding is not needed for the estimation of only the probability of collision. However, such information is now of importance for the severity \(f_s\), since otherwise one would not be able to differentiate between collision constellations. Hence, we require a method to combine the individual severity functions \(s_{j,l}\) into a total severity value \(s\), that can be used in (10 ), such that the domain of \(\theta\) remains disjoint but the information about the collision constellation is not discarded.
We approach deriving a total severity value by using a sum of the individual severities, which is averaged whenever multiple intersection angle intervals \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\) intersect. Therefore, we integrate the individual severity function \(s_{j,l}\) over the object’s velocity \(v_o\), as this variable is not considered in collision determination, giving the conditional expected individual severity with respect to the object’s velocity: \[\label{eq:expected95severity} \mathbb{E}[s_{j,l}(\theta, v_e, v_o)\mid v_o] = \int_{\mathcal{V}_o}p_{v_o}(\cmbv_o; \mathbf{w}_{v_o})s_{j,l}(\theta, v_e, \cmbv_o)\text{d}\cmbv_o.\tag{11}\] Using \(\mathbb{E}[s_{j,l}(\theta, v_e, v_o) \mid v_o]\), we construct the sum of average severities by following two rules. First, if a subinterval of a particular \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\) is disjoint from all other subintervals of \(\mathbb{I}^{\theta}_{u,v}(\phi, \rho)\), with \((j,l)\neq (u,v)\), then it’s severity is \(s_{j,l}\). Second, if an intersection among all \(N_eN_o\) intervals \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\) is found, the severity for the intersecting part is the average of all severity functions in that intersection. Thus, the severity will be given by a sum of disjoint and intersecting interval parts. We illustrate our approach to construct the sum of average severities with an example.

To calculate the sum of average severities, three operations are necessary. First, one must determine all intersections among all intervals \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\). Here, we distinguish the order of intersection index \(\tau\), where \(\tau\) indicates the number of intervals within the resulting disjoint interval, hence, \(0 <\tau \leq N_eN_o\) (e.g, \(\tau = 1\) for \(\tilde{\mathbb{I}}_5^{\theta}(\phi, \rho)\), or \(\tau = 3\) for \(\tilde{\mathbb{I}}_1^{\theta}(\phi, \rho)\) in Figure 3, last row). Second, one must ensure that each intersection is not included in another intersection. And, third, the average severity of each intersection must be calculated. To implement this procedure, we design an algorithm that starts identifying the intersection of the highest order \(\tau\) first and then identifies the remaining intersections in descending order, as lower-order intersections may already have been accounted for in a higher-order intersection. The algorithm is iterative. Every iteration ends when a new intersection or disjoint sub-inteval is found. Now, let \(m\) denote the current algorithm iteration. The three operations described above are performed as follows: first, we construct a set of tuples of indices \(\mathcal{L}_m\) preserving the information of which circle-to-circle pairs are intersecting. For example, \(\mathcal{L}_m = \{(j,l), (u,v)\}\) contains the indices of \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\) and \(\mathbb{I}^{\theta}_{u,v}(\phi, \rho)\). With \(\mathcal{L}_m\), we can identify the intersections of \(\mathbb{I}^{\theta}_{j,l}(\phi, \rho)\) as follows: \[\label{eq:intersections} \mathbb{I}_m^{\theta}(\phi, \rho) = \bigcap_{(j,l) \in \mathcal{L}_m }\mathbb{I}^\theta_{j,l}(\phi, \rho).\tag{12}\] Given \(\mathbb{I}_m^{\theta}(\phi, \rho)\), we execute the second operation, i.e., ensuring disjointedness of the intersection with any other intersections, \[\label{eq:disjoint95remainder} \tilde{\mathbb{I}}^\theta_m(\phi, \rho) = \mathbb{I}_m^{\theta}(\phi, \rho)\setminus\bigcup_{u \in \{1, 2, \dots, m -1\}}\tilde{\mathbb{I}}^\theta_{u}(\phi, \rho).\tag{13}\] For the third step, we average the severities by the cardinality of \(\mathcal{L}_m\), i.e., \(|\mathcal{L}_m| = \tau\), giving the average severity for that iteration \(m\), which, when considering the expected value of that, reads \[\label{eq:weighted95sum} \begin{align} \mathbb{E}[s_{m}(\theta, v_e, v_o) \mid v_o] = \frac{1}{|\mathcal{L}_m|} \sum_{(j,l) \in \mathcal{L}_m}\mathbb{E}[s_{j,l}(\theta, v_e, v_o)\mid v_o], \end{align}\tag{14}\] where \(\mathbb{E}[s_{j,l}(\theta, v_e, v_o)\mid v_o]\) is given by (11 ). With (11 ) and (14 ), we can express the inner integrals in (10 ) with respect to \(v_o, \theta\) as \[\label{eq:inner95integral} \begin{align} \sum_{m=1}^{N_{m}} \int_{\tilde{\mathbb{I}}^{\theta}_{m}(\cmbphi, \cmbrho)} \!p_{\theta}(\cmbtheta; \mathbf{w}_{\theta})\mathbb{E}&[s_m(\cmbtheta, v_e, v_o) \mid v_o] \text{d} \cmbtheta \\ &=: \mathbb{E}[s(\theta, v_e, v_o) \mid \theta, v_o, \phi, \rho \;]. \end{align}\tag{15}\]

We integrated operations (12 ) - (15 ) in Algorithm 4, which iterates over all intervals in descending order of intersections \(\tau\). Note that we abbreviated the expected values in Algorithm 4 and Table 1 for conciseness. Further, the function ‘permIndicies(\(\tau, m\))’ in line 6 of Algorithm 4 permutes the circle-to-circle indices \(j, l\), returns the set \(\mathcal{L}_m\), and assigns \(itr\) to check whether all possible permutations within the intersection order \(\tau\) have been calculated. In Table 1, we recursively apply Algorithm 4 to Example 2. Finally, the severity term in (10 ) is estimated with (11 ) - (14 ) and the result is re-introduced to (10 ) with (15 ), which gives the risk as \[\label{eq:ego95risk} \begin{align} &R^{cir}(v_e; \mathbf{w})= \int_0^{2\pi}\!\int_0^{\overline{\rho}}\!p_{\phi, \rho}(\cmbphi, \cmbrho; \mathbf{w}_{\phi, \rho})\\ & \Bigg[\sum_{m=1}^{N_{m}}\int_{\tilde{\mathbb{I}}^{\theta}_{m}(\cmbphi, \cmbrho)} \!p_{\theta}(\cmbtheta; \mathbf{w}_{\theta})\mathbb{E}[s_m(\cmbtheta, v_e, v_o) \mid v_o] \text{d} \cmbtheta\Bigg]\text{d}\cmbrho \text{d}\cmbphi. \end{align}\tag{16}\]

4 Risk Estimation for Gaussian Uncertainties and a Kinematic Severity Model↩︎

This section evaluates the risk from (16 ) for Gaussian uncertainties and a kinematic collision severity model, as the previous sections neither details the type of probability density functions nor the severity functions \(s_{j,l}\). We use Gaussians as those are often used (see, e.g., [11], [18]–[22]) to model kinematic uncertainties in motion planning, hence, we model \(x,y, v_o\) to be Gaussian random variables and assign a wrapped Gaussian distribution to \(\theta\). As a severity function, we will use a kinematic energy model that can accommodate different collision constellations, as the impact velocities and masses form a suitable predictor of collision severity (see, e.g., [23], [24]).

4.1 Gaussian Uncertainty↩︎

Although a relative heading angle \(\theta\) is defined on \(\mathbb{R}\), to any observer it appears periodic on the interval \([0, 2\pi)\). To account for that, we use a wrapped normal distribution to assign a Gaussian-like probability density function to the relative heading angle. With \(\mu_{\theta}, \sigma^2_{\theta}\) being the mean and variance of \(\theta\) composed in \(\mathbf{w}_{\theta} = (\mu_\theta, \sigma_\theta)\), let the wrapped-Gaussian for the relative heading angle be given by \[\label{eq:inf95wrapped95gauss} p_{\theta}(\theta; \mathbf{w}_{\theta})\!=\!\frac{1}{\sqrt{2\pi} \sigma_{\theta}}\!\sum_{\beta = -\infty}^{\infty}\!\exp{\left[ -\frac{(\theta\!+\!2\pi \beta\! -\!\mu_{\theta})^2}{2\sigma_{\theta}^2} \right]}.\tag{17}\] As shown in [25], the infinite sum in (17 ) can be tightly approximated by truncating it within \(\beta \in \mathbb{Z}_{-N_\beta}^{N_\beta}\), when \(N_\beta \geq 3\), leaving \(2N_\beta + 1\) summands⁴.
Given the polar coordinate transformation (6 ), we retrieve the bivariate-Gaussian for the positional coordinates \(\mathbf{x}\) as a function of the polar coordinates \(\phi\) and \(\rho\) as follows: \[\label{eq:polar95gauss} \begin{align} & p_{\phi, \rho}(\phi, \rho; \mathbf{w}_{\phi, \rho}) = \frac{\rho}{2 \pi \sigma_{x}\sigma_{y}} \text{exp} \Bigg[ \frac{(\rho \cos{(\phi)} - \mu_{x})^2}{-2\sigma_{x}^2} \\ & -\frac{(\rho \sin{(\phi)} - \mu_{y})^2}{2\sigma_{y}^2} \Bigg], \end{align}\tag{18}\] where \(\mu_{x}, \mu_{y}, \sigma^2_{x}, \sigma^2_{y}\) are the means and variances of \(x, y\), composed in \(\mathbf{w}_{\phi, \rho} = ((\mu_x, \sigma_x),(\mu_y, \sigma_y))^T\) and we note that the integral over (18 ) in (16 ) requires numerical integration (see [26]).
Lastly, the uncertainty in the velocity is given by a univariate Gaussian with mean \(\mu_{v_o}\) and variance \(\sigma_{v_o}^2\), where the cumulative distribution function with \(\mathcal{V}_o = [\underline{v}_o, \overline{v}_o]\) reads \[\label{eq:analytic95univariant95gauss} \begin{align} &F_v(\mathcal{V}_o; \mathbf{w}_{v_o}) := \int_{\mathcal{V}_o}p_{v_o}(\cmbv_o; \mathbf{w}_{v_o})\text{d}\cmbv_o\\ &=\frac{1}{2}\left[\text{erf}\left(\frac{\overline{v}_o- \mu_{v_o}}{\sigma_{v_o} \sqrt{2}}\right)-\text{erf}\left(\frac{\underline{v}_o- \mu_{v_o}}{\sigma_{v_o} \sqrt{2}}\right)\right], \end{align}\tag{19}\] where \(\text{erf}(a)\) denotes the error function of \(a\). Note that (19 ) will be a useful expression to derive the conditional expected severity next.

4.2 Severity Function↩︎

Whether an analytical expression for (16 ) can be found depends on the form of the severity function \(s_{j,l}\) and the uncertainty characterization in terms of the PDFs. The directionality of each actor’s velocity is, however, challenging to consider when using Gaussian uncertainty, as an analytical solution to (15 ) may not be obtainable. The proposed multi-circular model, however, can accommodate different severity functions for different collision constellations regardless of the heading angle. We, therefore, propose a kinetic energy model which weighs collisions differently, depending on the specific circle-to-circle collision. Hence, we have \[\label{eq:severity} s_{j,l}(v_e, v_o) = \frac{w_{j,l}m_em_o}{2(m_e + m_o)}\cdot\begin{cases} (v^2_e + v^2_o) \text{ head-on coll.,}\\ v_e^2 \text{ ego-to-obj. side coll.,}\\ v_o^2 \text{ obj.-to-ego side coll.,}\\ (v^2_e - v^2_o) \text{ ego rear-end coll.,}\\ (v^2_o - v^2_e) \text{ obj. rear-end coll.,} \end{cases}\tag{20}\] with \(m_e, m_o\) being the ego and object mass, respectively. The cases in (20 ) are assigned based on the specific circle-to-circle collision pairs⁵. In addition, \(w_{j,l}\) denotes a positive scalar weighting factor that accounts for the severity of the \(j\)-th ego circle colliding with the \(l\)-th object circle. With the severity function and the integral of the univariate Gaussian in the velocity (19 ), we can derive the conditional expected value of (20 ), with respect to the object’s velocity \(v_o\), by using two integrals, \(I_e\) and \(I_o\), respectively. The first integral \(I_e\) represents the integral of the univariate Gaussian \(p_{v_o}\) multiplied with the constant \(v_e^2\), see (19 ). The second integral, \(I_o\), is a product of \(v_o^2\) and \(p_{v_o}\), where the solution is known, e.g., see [27]. Hence, the conditional expected circle-to-circle severity is given by: \[\label{eq:analytic95severity} \begin{align} &\mathbb{E}[s_{j,l}(v_e, v_o) \mid v_o] = c_1\cdot \begin{cases}(I_e + I_o)\text{ head-on coll.,}\\ I_o \text{ ego-to-obj. side coll.},\\ I_e \text{ obj.-to-ego side coll.}, \\ (I_e - I_o) \text{ ego rear-end coll.,}\\ (I_o - I_e) \text{ obj rear-end coll.,} \end{cases}\\ &\text{where: } I_e = \frac{v_e^2}{2}\left[\text{erf}\left( \frac{\overline{v}_o - \mu_{v_o}}{\sigma_{v_o} \sqrt{2}}\right) - \text{erf}\left( \frac{\underline{v}_o - \mu_{v_o}}{\sigma_{v_o} \sqrt{2}}\right) \right],\\ & I_o = \frac{1}{\sigma_{v_o} c_3^3}\Bigg[ \frac{c_2^2 + 1}{2} \left( \text{erf}\left(\frac{c_2 + c_3\overline{v}_o}{\sqrt{2}} \right) - \text{erf}\left(\frac{c_2 + c_3\underline{v}_o}{\sqrt{2}} \right)\right) \\ & \quad + \frac{1}{\sqrt{2\pi}}\Bigg((c_2 - c_3\overline{v}_o) \exp{\left(\frac{(c_2+c_3\overline{v}_o)^2}{-2} \right)}\\ & \quad - (c_2 - c_3\underline{v}_o)\exp{\left(\frac{(c_2+c_3\underline{v}_o)^2}{-2} \right)}\Bigg)\Bigg],\\ & c_1 = \frac{w_{j,l}m_em_o}{2(m_e + m_o)}, \quad c_2 = \frac{-\mu_{v_o}}{\sigma_{v_o}}, \quad c_3 = \frac{1}{\sigma_{v_o}}. \end{align}\tag{21}\] Note that we substitute (21 ) (albeit without the argument of \(\theta\)) in (14 ) to retrieve \(\mathbb{E}[s(v_e, v_o)\mid v_o]\) with Algorithm 4.

4.3 Risk Estimation↩︎

While (20 ) does not depend on \(\theta\) (which may not be case in general, see (15 )), it is still part of the uncertainty in the configuration. However, by assumption, \(\theta\) is independent of the positional coordinates \(\phi, \rho\). Hence, we can solve the integral of \(\theta\) in (16 ) irrespective of the severity term and the the positional coordinates \(\phi, \rho\). Therefore, we analytically integrate (17 ) over \(\bigcup_{m = 1}^{N_m}\tilde{\mathbb{I}}_{m}^{\theta}(\phi, \rho)\), given the approximation \(-N_\beta \leq \beta \leq N_\beta\) in the sum in (17 ). Hence, for the conditional expected average severity (15 ) we have \[\label{eq:analytic95warpped95gauss} \begin{align} & \mathbb{E}[s(v_e, v_o) \mid \theta, v_o, \phi, \rho \;] = \\ &\sum_{m = 1}^{N_m}\Bigg\{\sum_{\beta=-N_\beta}^{N_\beta}\Bigg[\text{erf}\left( \frac{\overline{\theta}_{m}(\phi, \rho)\!-\!\mu_{\theta}\!+\!2\pi \beta}{\sigma_{\theta} \sqrt{2}}\right)\\ & - \! \text{erf}\left( \frac{\underline{\theta}_{m}(\phi, \rho)\!-\!\mu_{\theta}\!+\!2\pi \beta}{\sigma_{\theta} \sqrt{2}}\right)\!\Bigg]\mathbb{E}[s_m(v_e, v_o) \mid v_o]\Bigg\}. \end{align}\tag{22}\] With (22 ), the risk from (16 ) is given by \[\label{eq:ego95risk95short} \begin{align} R^{cir}(v_e; \mathbf{w}) = \int_0^{2\pi}\int_0^{\overline{\rho}}&p_{\phi, \rho}(\cmbphi, \cmbrho; \mathbf{w}_{\phi, \rho})\\ & \mathbb{E}[s(v_e, v_o) \mid \theta, v_o, \cmbphi, \cmbrho \;] \text{d}\cmbrho \text{d}\cmbphi, \end{align}\tag{23}\] where we solve the remaining two-dimensional integral numerically⁶ using the expressions in (22 ) and (18 )⁷.

5 Numerical Case Studies↩︎

In this section, we present a numerical case study to demonstrate the proposed methodology for risk estimation in different collision scenarios and assess its computational performance.

5.1 Collision Scenarios↩︎

A key contribution of this paper is that the proposed risk function features a severity model that takes different impact constellations into account by associating them with different severity values. Therefore, we design five exemplary test cases that demonstrate the efficacy of the risk function in frequently occurring real-world accidents, i.e., head-on, rear-end, and side collisions with varying impact points. Throughout all scenarios, the motion of each vehicle is (assumed to be) unaffected by the motion of the other, allowing both vehicles to travel through each other. Note that the initial configurations, velocities, and uncertainties are provided in 3, where the configurations are in a global coordinate frame. Therefore, we distinguish the ego’s and object’s variables by \(e\) and \(o\) subscripts, respectively. All velocities are kept constant for each test case. Cases (I) and (II) examine the risk function in a head-on and a rear-end collision scenario, respectively, where the object is traveling along the ego’s longitudinal axis, see 5. Cases (III), (IV), and (V) consider side collisions, where the object collides with the ego at a \(90\,\text{deg}\) angle at three different impact points, i.e., front, center, and rear, respectively (see 6). Case (III) focuses on a collision with the front of the ego, whilst Case (IV) describes a collision with the ego’s driver-sided door. Lastly, Case (V) considers an impact towards the back of the ego. For each test case, we keep the uncertainty, i.e., \(\sigma_x, \sigma_y, \sigma_\theta, \sigma_v\), constant. The motion of each vehicle is modeled by constant velocity model, and each vehicle’s shape is approximated by three circles. The three-circle approximation of each actor’s shape yields nine severity weighting factors \(w_{j,l}\), that we collect in a matrix \(W\). We assign the weights such that front-to-side collisions are weighted the strongest, followed by frontal, and, lastly, rear collisions. All simulation parameters are provided in Table 4.

5.2 Computational Performance↩︎

As the risk function is based on a method for collision probability estimation, we compare the computational performance of the risk function to the algorithm for collision probability estimation of [14]. To estimate the risk we use the implementation of [14], however, we substitute Algorithm 3 "Sort Intervals ‘SortInts’" in [14] with our Algorithm 4, to calculate risk instead. Note that both resulting algorithms consist of an initialization step and an estimation step. The initialization step only has to be performed once for a given geometry and number of circles of both vehicles, whereas the estimation step is typically carried out at high frequency at runtime (e.g., in a motion planning context). The initialization step calculates the disjoint intersection angle intervals \(\breve{\mathbb{I}}^{\theta}(\phi, \rho)\) and, for the risk estimation algorithm, the expected conditional severity per interval \(\mathbb{E}[s_m(\cmbtheta, v_e, v_o) \mid v_o]\). Once initialized, the estimation step is only left with computing (22 ) and (23 ). We calculate the average initialization time over 100 initializations and the average estimation time over \(10^4\) estimations. We measure both times for a two, three, and four circle vehicle approximation. Furthermore, we randomly vary both vehicle configurations, velocities, and uncertainties between each function call, so as not to calculate the same risk and collision probability for each measurement. As a computational platform we use a notebook using an Intel Core i7-12800HX processor clocked at up to \(4.80\,\text{GHz}\) with \(32\,\text{GB}\) of memory.

5.3 Results↩︎

5.3.0.1 Risk estimation for varying collision constellations

The results for all test cases are given in 5 and 6. Here, we let the vehicles traverse through the respective Cases (I) - (V) and estimate the risk, after a single initialization for each test case, at a frequency of 1 ms. The graphs are aligned such that at 0 s both vehicles touch for the first time for the mean motion of the object, i.e., the collision begins. Regarding Cases (I) and (II), we report maximum risk values of \(3.3 \cdot 10^5\) in the head-on case and \(3.1 \cdot 10^5\) in the rear-end case. Note that both vehicles are traveling through each other, and, hence, at collision impact, i.e., at 0 s, the rear-end case reports a risk of \(1.5 \cdot 10^5\) and the head-on case is estimating a considerably higher risk at \(2.0 \cdot 10^5\).
The vertical markings over the \(x\)-axis in Figures 5 and 6 give the respective collision endings for the mean object motion, i.e., at that time, both vehicles stop colliding. The head-on collision reports the highest maximum risk shortly after the start of the collision, as (21 ) measures the head-on case with the absolute kinematic energy, and a collision between a front and a center circle is weighted the strongest (which starts to occur shortly after the collision begins). The graph of the head-on cases rises later as the rear-end case due to the higher relative velocity, i.e., the distance between both vehicles reduces in less time. Note that in the rear-end case, a head-on collision is also occurring when the ego travels through the object (hence, both front circles intersect at a later point); however, here both vehicles are also colliding with their other, lower-weighted sections, where the total severity value is then averaged. This combination explains the delayed, slightly lower peak in the rear-end case. Whether such behavior is desirable depends on the application of the risk estimation, e.g., one may restrict risk estimation to some brief interval until the start of a collision.
We report a similar effect for Cases (III)-(V), in which the different criticalities of potential collisions become visible. Associated with the highest risk is Case (IV) at \(1.21\cdot 10^5\), followed by Case (III) with \(1.07\cdot 10^5\) and, finally, Case (V) at \(0.69\cdot 10^5\). The reduction of about \(12\,\%\) from side-to-front collisions and another \(64\,\%\) from front-to-rear impacts showcases the design of the severity matrix \(W\), that weighs front-to-side collisions as the most severe.

5.3.0.2 Computational performance

2 reports the average initialization and estimation times for risk and collision probability calculations, as well as their difference. For the initialization time, we find an increase by a factor of about \(4.4\) for the circle added to the ego and the object. Concerning risk and collision probability estimation times, a speed-up of \(60\,\text{\%}\) is found for reducing the number of circles used for vehicle approximation per vehicle from three to two circles. Meanwhile, the differences between risk and collision probability calculation times are negligibly low at maximum \(8.13\,\mu\text{s}\). Note that for many applications, e.g., motion planning, the initialization times are inconsequential, as collision probability or risk estimation algorithms are initialized before run-time. Furthermore, it is possible to prepare multiple risk functions, which are then only called from memory instead of being initialized whenever needed. Given the fact that both risk and collision probability estimation are calculated in little time, both algorithms can be suitable for use in real-time motion planning.

6 Conclusions and Future Work↩︎

In this paper, we propose a risk estimation methodology for automated vehicles, in which risk reflects both the probability and severity of a collision event. The risk estimation approach is shown have good computational performance and is hence suitable to be deployed in, e.g., motion planning algorithms or for general safety assessment. The proposed method allows us to assign different collision severities for different collision constellations. To achieve this, we extend the method for collision probability estimation of [14] to risk estimation, inheriting us computational efficiency and estimation accuracy as we utilize the same algorithmic structure. Furthermore, we provide an exemplary derivation for risk estimation for the scenario of one vehicle facing another vehicle with Gaussian uncertainties in the position, heading angle, and velocity, along with a kinematic energy model that considers the collision severity differently between different collision constellations. Given the severity model and Gaussian uncertainties, we can reduce the computational effort of risk estimation significantly by reducing the four-dimensional integral to a two-dimensional integral, allowing us to leverage established efficient numerical integration methods. We demonstrate in a representative example that our risk estimation approach indeed allows us to assign different collision severities for different collision constellations. In addition, we show that our algorithm is computationally efficiently enough for real-time motion planning applications. Lastly, we provide the programming code for the risk function as open source.
For future work, we envision to investigate other models for the uncertainties, which we assumed to be Gaussian for exemplary purposes and due to their frequent occurrence. Similarly, other and higher-fidelity severity models may be employed to represent collision severity more accurately. Here, we note that especially the resulting combination of probability density functions and severity functions is interesting, as the computational efficiency depends on the functional forms of both, i.e., whether arising integrals can be analytically solved. As our results demonstrate the efficacy with respect to computational efficiency and severity estimation under different collision constellations, we clearly recommend the use for motion planning applications.

References↩︎