Deriving a joint risk estimate from dynamic data collected at motorcycle rides

.


Introduction
As a fun and recreational activity, motorcycle usage has been steadily increasing. Given the motivation for these kinds of rides, they are unlikely to be replaced by any kind of automated vehicle, as controlling the bike itself is part of the enjoyment (see e.g. Will et al., 2019;Broughton and Stradling, 2005). Yet, when it comes to safety, motorcycle riders are much more vulnerable to suffer serious injuries from accidents (see EC, 2020). For cars, the advent of advanced assistant systems has affected the accident rate positively (see Golias et al., 2002). The same level of improvement has not been demonstrated for motorcycle rides in the past (see Savino et al., 2020), leading to a need for developing further systems to improve motorcyclist's safety and ultimately prevent fatalities among motorcycle riders in general (see for instance . Several physical injury prevention measures are being considered, such as helmet usage (see Siebert and Lin, 2020;Tabary et al., 2021;Meng et al., 2020;Patalak et al., 2020), vehicle design (see Xiao et al., 2020) or back protection (see Afquir et al., 2020). Yet, determining suitable interventions to the system motorcycle-rider, or even recognizing a critical situation for the rider in the first place is an ongoing challenge. Several studies aim to remedy these issues (see Ryder et al., 2016;Cuenca et al., 2018). In particular, using previous accident data will play a key role in trying to find recognizable features of potential accident risk spots (see examples in Li et al., 2021;Bíl et al., 2013;Yalcin and Duzgun, 2015), alongside with traffic conflict detection and rider self-reports. Yet, while these studies focus on a mix of driving, accident and infrastructure data (see for instance Afghari et al., 2020), the complex relation between individual motorcycle riding behaviour and the occurrence of accident hot spots makes it very difficult to obtain generalizable risk estimates and predict potential future accident hot spots.
We propose a novel method for the determination of motorcycle riding risks from recorded motorcycle dynamics data, obtained by rides of a set of motorcycle riders, based on a combination of several well established methods known from machine learning. Our approach derives individual risk thresholds and weights from the measured dynamic data of all riders. In continuation of this approach it is possible to use the "dangerous spot" classifications of each ride and rider and build an overlay/local sum of the various estimates. Combining this overlay with another bound for joint danger estimates, it is possible to derive a joint risk estimate across riders and individual rides to determine overall risk estimates associated to a specific road section, rather than to a specific rider or ride.
The necessary specific dynamics data can be measured with the help of a motorcycle equipped with special systems that allow the collection of high-quality GPS-signals and vehicle-specific dynamic data over time. We will discuss however, how this condition might be alleviated in the future, by using more readily available systems.
The proposed method of detecting potentially dangerous street sections (hot spots) could be used to improve motorcycle rider safety in several ways: Firstly, to provide riders with a map of upcoming known and predicted dangerous spots during a ride, to help them focus their attention in advance on an upcoming potential danger. Secondly, the method of fitted individual risk indicators could be used to determine individual danger or close-to-danger estimates online and provide assistance to the rider, provided that suitable non-distracting assistant systems can be devised. Moreover, dangerous spots jointly located by several riders can be used to investigate and evaluate these hot spots in regard of the reasons for their classification as dangerous. By providing that information to road maintainers they can possibly remove the causes and can make the road infrastructure safer for motorcycle riders.

Data acquisition
The presented approach to risk spot detection requires data acquisition of dynamic signals during motorcycle rides. Such signals may include position data acquired by a GPS-receiver (Geopositioning systems) and dynamic data measured by IMUs (Inertial Measurement Unit, see Waegli et al., 2008) attached to the motorcycle. Built-in IMUs, which are already available on premium motorcycles, can be accessed by a Controller Area Network (CAN) bus, see International Organization for Standardization, 2015. The collected time series measured by an IMU are typically the accelerations in x-, y-, and z-direction and the angular velocities (rates) about the roll-(x-), the pitch-(y-), and the yaw-(z-) axis (see Tseng et al., 2007). The resolution of the obtained time series will of course depend on the sampling rate of the data acquisition by the measuring system.
In this work, data from 5 different drivers on an instrumented motorcycle will be used to illustrate a possible realisation of this procedure. Details on the vehicle can be found in Ecker and Saleh (2016) and Schwieger et al. (2018). Measurements took place on 6 different sections of various regular country roads (all of them popular motorcyclist routes) in (Lower) Austria. Not all riders drove each section, rather at least 3 different riders drove at least 5 times on every track. The length and name of each track can be found in Table 1.
Anonymized accident data (for the years 2012 to 2015), comprising location, accident type and driving direction (if known), was retrieved from the national accident records for the 6 relevant road section and will be represented as a set AC with elements ac ∈ AC. All quantities pertaining to a particular accident location will be indexed by ac.

Ethics statement
The drivers gave informed consent to participate in this experiment and allow the collected data to be presented. All of them were experienced motorcycle riders and were monitored for safety during their rides.

Variables of interest and preprocessing
In the presented approach dynamic variables of interest were the Yaw-rate Υ, Roll-rate Р, Pitch-rate Ψ (physical units in use: deg/s, see Tseng et al., 2007), together with the time-derivative κ of the driven curvature. These variables describe the (signed) speed of change along the principal rotations of the vehicle: the Yaw-rate is the change of the yaw angle in time, i.e. the angular velocity about the z-axis (the axis pointing towards the ground, i.e. a left or right turning movement), the Roll-rate is the change of the roll angle in time, i.e. the angular velocity about the x-axis (the axis pointing towards the vehicle's front, i.e. a left or right leaning movement) and the Pitch-rate is the change of the pitch angle in time, i.e. the angular velocity about the y-axis (the axis pointing towards the right of the moving vehicle, i.e. a raising or lowering of the vehicle's front). For the calculation of κ the constant of 0.1 is added in the denominator to avoid numerical problems when a division by 0 would occur: In order to compare the variations in the driving behaviour of multiple rides on the same road section, a possible approach is to project the time series of the measured variables on a longitudinal scale, such that a per-distance value is calculated for the road section of interest. Thereby, the independent variable "time" of the measured series is changed to the new independent variable "distance" d. This mapping clearly leaves some degrees of freedom in the choice of how this projection is carried out. A standard method would be, to choose for each increment of the longitudinal unit (e.g. one Meter) on the new scale a reference point. Then one can assign to each unit increment the mean values of all measured parameters. This will be the approach in use and as denoted above, one Meter will be used as the unit increment.
Together, these variables form a vector for every given longitudinal increment in the time series and form a 4 dimensional vector from those. Once the measured dynamic parameters are averaged and represented as series of values on a longitudinal scale, data smoothing is beneficial to reduce the signal noise caused by measurement errors. A common method to achieve this goal is the computation of a rolling average over a small number r 1 of data points.
Eq. (3) calculates the simple moving average (SMA), which is the unweighted mean of the previous r 1 of data points. For a risk prediction system working in real-time on a vehicle, this would be the most simple filter to be applied. With off-line data processing, more sophisticated filtering methods may be used.
We will refer to the per distance smoothed data as Ξ(d) for a given position d on the track: (4)

Table 1
The 6 names of the tracks with the length of the respective tracks in meters.
It is convenient to split the time series in 2 series, one containing only positive and one for negative data values, respectively. This separation allows for different weighting factors for increasing, or decreasing parameter values and is helpful with roll angle and yaw angle data processing as well as with the direction of motion to the left or to the right. This data splitting process yields new vectors for positive signed data and for negative signed data Based on these movement trajectories, first differences of the variables can be calculated as approximate derivatives and used to characterize local rider motions more precisely: The remark above about weighting positive and negative valued time series separately, applies here too, hence giving Furthermore, with ID denoting the set of all drivers, in order to make reference to the driving dynamics of a particular driver id ∈ ID, we use the notation Ξ id (d) (and accordingly for the components and derived quantities). All quantities pertaining to a specific driver will be indexed by id.

Clustering and dynamic value
The smoothed values and derivatives of an individual driver can now be used to find clusters (see Scrucca et al., 2016) within the data. By applying e.g. a k-nearest neighbour algorithm (see Ling, 1972), the cluster centers can be used to define "central" or "standard" riding behaviour of a certain motorcycle rider. The resulting cluster centers are summarized in the set Clust id C with elements q id ∈ Clust id C . Cluster membership of individual elements of the time series can be stored and illustrated as in Fig. 1.
These cluster centers can then be compared to the same rider's values at a given accident site ac. A separation between the values at accident sites and the standard motions can be obtained by a multitude of separation methods. For instance, a linear model or a neural network (see Bishop, 2006) can be fit to the dynamics data by mapping the standard motions of a rider to 0 and the values at accident spots to 1. This will yield weights on the effects of the dynamics variables, and by using these weights, it is possible to associate a value to every set of dynamics values in the data set. The mapping derived by this extension of a separation of dangerous and non-dangerous spots, will be denoted by ℓ id , with its values named "dynamic value".

Individual risk estimate and individual threshold
The target value of ℓ id at accident spots would be 1, yet it is clear that an exact fit will not be possible, nor reasonable, given that the same dangerous spot is traversed several times with different motions. We obtain an initial separation value (termed "individual risk threshold") between risky and non-risky dynamics ℓ id b via the average ℓ id value at accident spots i.e.
with |AC| denoting the number of accident spots in AC (counted once for every drive-by in this situation). Using this individual cut-off, we could classify the dynamics vector of a position d on the track as being in the "risky" domain iff

Fig. 1.
Illustration of cluster membership of points along an example ride for a single ride, displayed in the dimension of Yaw-rate Υ(d). These example clusters correspond to a strongly negative Yaw-rate (cluster 1, Pink), a sign-change between spiking Yaw-rates (cluster 2, Green), a strongly positive Yawrate (cluster 3, Red), a moderately negative Yaw-rate (cluster 4, yellow) and a moderately sized evolution of the Yaw-rate around zero (cluster 5, Black). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) This defines a binary map on the track being considered, with the value 1 being assumed at "risky" spots and the value 0 being assumed everywhere else. We go one step further and use the boundary value ℓ id b to define a non-binary map l id of an individual local accident risk: This map was defined such that the relative values of how close the risk estimate was to the bounding ℓ id b ) could be included in the risk estimate in a non-binary way, to avoid the case of a "close to threshold" drive contributing the same 0 amount as a "far from threshold" value.
The fundamental data components Υ(d), Υ + (d), Υ − (d) as well as l id and l id b are illustrated in Fig. 2.

Joint estimate of accident risk spot locations
Following the determination of the individual risk estimates for all rides and riders, all given risk estimates for a given position on the track can be summed up and divided by the number of times a valid measurement was obtained on that position. The values of this overlay of individual limit transgressions need to be compared to a final "group level threshold b", determining whether the observed individual risk estimates were sufficient to qualify a position on the track as potentially dangerous. This threshold can be determined by using the trade-off between accident sites correctly found as dangerous spots, which can be encoded as a probability Θ and the amount of the track being classified as dangerous, encoded as a probability A , with the aim of keeping the latter low but the former high. An example output of this approach can be seen in Fig. 3. Values of the overlay that lie above this threshold are then classified as dangerous and a map of dangerous spots on the track has been obtained, as is shown in Fig. 4.

Result
We used a procedure of initially linear separation (see Bishop, 2006;MacQueen, 1967), similar to the classic linear perceptron (see Rosenblatt, 1958), to obtain an initial ℓ id per driver, followed by a non-linear optimisation to get the best possible trade-off between minimizing the area classified as dangerous and maximising the proportion of accident spots classified correctly as dangerous via our model. The implementation was carried out in the language and environment for statistical computing, named R (R Core Team, 2018; Zeileis and Grothendieck,  The linear models (one for each id ∈ ID) to determine the coefficients yielding 17 coefficients (per driver). These coefficients were used to define ℓ id as functionals on the driving dynamics data Ξ id (d) for every position d: Given the ℓ id , we can derive the l id , as defined in Eq. (11). If we sum up the l id (Ξ id (d)) of all recorded rides on the track and divide by the number N id (d) of valid measurements per driver d id 1 , d id 2 , …, d id N id (d) on each position d of the track, we get a weight function of risk estimates on elements of the track. This weight function was smoothed by a rolling average of r 2 distance units again, to provide a regular transition between high risk and low risk domains. The weight function together with an additional boundary b yields a final map of "risk spots" d r on the track of interest. Every risk spot d r is assigned a "surrounding area" A r of r 1 +r 2 distance units before and after.
If |A| denotes the number of elements of a set A and the set of individual positions on a track is defined as D track , we can obtain a fraction of distance units A track ∈ [0, 1] included in "risk spot areas". Letting ⋃ denote the union of different sets, the fraction of risky area distance units of the whole track (all distance units) is For multiple tracks one can proceed in the same way and define A jointly over all tracks.
Finally, we explore which accident spots ac where included in any of the areas A r and calculate the fraction Θ ∈ [0, 1] of accident spots hit by the model. To that end, we use the set intersection ∩ of accident spots AC with the risky area output by the model: Analogously, Θ track can be defined using only risk spots and accidents for a given track.
After the above set-up steps have been taken care of, the final step in getting the model estimates consists of an optimisation procedure, with the goal of maximizing Θ while minimizing A i.e. obtain precise estimates which nevertheless include most known accident spots in their warnings. We choose a target function C to encode this trade-off in the form: where ‖.‖ 2 denotes the standard ℓ 2 -norm (root mean square). The input variables for the estimation were all ℓ id b and the 17 regression coefficients per rider and the value b for the joint estimate. The target was to minimize C. We used a smoothing range of r 1 = 5 distance units and r 2 = 5 distance units, after precursory experiments with several ranges from 5 to 40 distance units (in steps of 5). We used an implementation of the above model in R and carried out the optimisation with the "optim" procedure included in R, at a maximum of 30000 iterations. We investigated the issue of whether estimates derived on a subset of tracks could be used to infer dangerous spots on tracks not previously used to fit the model. To this end we used any combination of 5 tracks to fit the model as above and then carry out predictions on all 6 tracks. The resulting quality of fit (with A and Θ expressed in percent)  Fig. 4. Example track (in turquoise) with the resulting map of risk areas (light red for risk and dark red for high risk) derived from the overlay of risk estimates shown just before, together with a small lead up and follow up domain for each risk spot. Known accident spots are indicated by black hexagons. Areas of high risk are highlighted additionally by red stars, at the last meter of the high risk area. Driving direction is indicated by a black arrow, going from the upper left to the lower right on this map. on 6 tracks in two directions each can be seen in Table 2, left table. We then used the outcomes of the 6 "leave-one-out" fits as additional starting values for the fit of the model on all tracks.
For comparison, the values obtained when leaving out e.g. the Adamstal track can be found in Table 2 (middle) and they are close to the original 6 track fit, though the values of the 6 track fit produce a smaller risky area on most tracks, while finding a similar proportion of accidents, making it perform slightly better overall.
Conversely, leaving out the Exelberg track leads to very poor performance on this track (see Table 2, right table). This implies, that in the current 6 track data set, not all fits can be generalized to the other tracks.

Interpretation of result and outlook
We have presented a methodology for predicting dangerous spots for motorcycle rides from the vehicle dynamics data of several test riders on 6 given tracks. This was achieved by the comparison of vehicle dynamics data at known accident spots and representatives of individual clusters of the dynamics values on the rest of the tracks. The methodology relies on the use of smoothing, clustering, regression and classification/separation, with a follow-up step of using sum statistics. To show transferability the methodology was trained on 5 of the 6 tracks and the quality of the derived criteria was demonstrated empirically. Criteria for "dangerous spot"-like dynamics were derived individually and "dangerous" classifications were summed up into an overlay between all rides on a given track, to achieve a measure of riding risk independent of a given individual rider. The fit to the given accident data was much improved over earlier versions of the method (see Hula et al., 2019), with the model now typically finding 90% of the accident spots, while covering only between 6% and 38% (Kalte Kuchl track) of the tracks with a "potentially risky" qualifier.
With regards to transferability of the criteria from one track to another, it appears that some tracks can be more readily removed from the training data than others (see the case of leaving out the Exelbergtrack in the results section). This suggests that a set of exemplary tracks should be derived in more extensive data studies in the future and that the study of which track features determine the degree of transferability would be a potential further research question.
The presented methodology provides an significantly optimised extension of the method presented in Hula et al. (2019) and Schwieger et al. (2020), in that the coefficients obtained in the initial separation approach were optimised again, based on the criterion in Eq. (18). Also, contrary to the cited earlier work, risky spots were no longer defined as maxima of the overlay of individual risk estimates, but rather all values above a certain threshold were defined as risky.
The risk thresholds, both individual and the group level risk threshold, were obtained from the optimisation of our target function in Eq. (18). In this idiosyncratic form they are not a priori associated to risk considerations from naturalistic driving studies in the literature (see Singh and Kathuria, 2021;Khattak et al., 2021;Arvin and Khattak, 2020;Eenink et al., 2014). Such an association could however be very beneficial to the interpretation (see "Limitations" below) and further development of the presented results, hence it would be an important venue for further research.
The resultant individual regression functions or weights and "risk thresholds" provide a means of creating a profile of a given rider and could be a first step to actively providing assistance to the rider even online during a motorcycle trip and augment other real time systems (for instance conflict detection such as in Formosa et al., 2020). As noted, there is some difficulty as to what measures of assistance may be feasible. Warnings to the rider may prove distracting but the same could be said about any noticeable change in driving dynamics, leaving this an open issue.
The obtained map of dangerous spots can be used to warn motorcycle riders of upcoming risk spots. As this warning can be readily given to drivers in advance, compared to a potentially surprising warning based on current driving data, this use-case appears to have fewer issues in its practical implementation compared to an online estimate. It could help to ensure that riders focus their attention on the upcoming risk spot and thus improve safety significantly.
Furthermore, the map of predicted dangerous spots gives an indication as to where a road safety inspection to assess safety risks through expert judgement can be recommended, thus allowing for a more targeted investigation of road infrastructure, in particular with respect to popular motorcycling tracks.
It would be desirable to allow large data sets of motorcycle rides to be investigated, rather than use an experimental set up with few riders on a specially equipped bike. This will require the methodology to work with less sophisticated measurement systems, which will be addressed in future work. The methodology as presented has been granted a patent (see Hula, 2020) for the purpose of the detection of critical road sections.

Limitations
The approach is limited by the quality of the data from accident reporting (in our case location and accident type, also see Hauer, 2020 for a perspective on crash cause investigation) as well as, as of now, a quite limited data sample (6 tracks, 5 drivers). The first issue needs careful auditing of the data, possibly including an in-depth reconstruction of accident causes. The latter issue could be resolved by a large-scale data collection, in particular the use of on-board systems (see Chaczko et al., 2017) like IMUs and in many cases GPS. That way it might be possible to apply the given methodology to larger data sets and create individual profiles of a large set of drivers. Of course, it may be necessary to create overlays not on a per-distance basis but GPS location proximity. This is in principle well possible, although further work will have to investigate this.
While transferability on the 6 different tracks has been investigated by "leave one out" model fits and subsequent evaluation on the whole data set, transferability on larger and more diverse sets of tracks and riders will require ongoing testing and adaptation and possibly new methods find well generalizable subsets (see another methodology in Tang et al., 2020 in the context of other risk models and Feng et al., 2020 for a perspective on international transferability).
We note that our 5 drivers being experienced motorcyclists likely affects the collected data in that we would expect less experienced riders to be more variable/less secure in the driving patterns that we observe. This may have made it easier to derive a joint estimate of risky areas and might be more challenging in a broad data collection based on commercial on-board systems, as suggested above. However, in this case the larger empirical basis should stabilise the estimates despite the additional variability and if that were to prove insufficient, the individual driver models could be refined to deal with a larger number of "representative" movements or a non-linear separation model could be employed (see below).
Furthermore, the classification of a spot on the track as "dangerous" will not immediately yield an interpretable reason for the spot being found. Turning the contributions of different dynamics variables and their approximate derivatives into an interpretation of why a given spot is dangerous will require further work. Also note, that we have not yet included the effects of speed on accident severity and/or accident risk, which has been discussed in recent studies (see Murphy and Morris, 2020), but rather treated all previous accident spots as equal regardless of accident severity. Additionally, accidents might have to be split into groups according to weather conditions and daytime (see Robbins and Fotios, 2020) at the time of occurence.
The optimisation problem used in this work to fit the model is quite high-dimensional and the procedure risks being stuck in local minima. We suggest defining a larger set of starting values to avoid this issue.
The separation model in use is based on an underlying linear separation, but this could potentially be replaced by nonlinear separation approaches (see Cortes and nd Vapnik, 1995;Tipping, 2001), although it would have to be investigated carefully, whether this added complexity and computational cost would improve estimates enough to warrant such a change.
Finally, continued usage on a changing road infrastructure will require a mechanism for repeatedly adjusting the model estimates to infrastructure changes, such as discounting further back estimates compared to more recent ones or running explicit new calibrations.

Declaration of Competing Interest
A patent has been obtained on the usage of the described method in the detection of critical road sections (see Hula, 2020). The motorcycle used in this study has been provided by KTM. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.