NGC6388 (*) | $8.3\times 10^{4}$ | $6.0\times 10^{-13}$ | NGC6397 | $2.5\times 10^{5}$ | $4.1\times 10^{-12}$ | NGC6401 (*) | $1.6\times 10^{5}$ | $6.5\times 10^{-13}$
NGC6402 (*) | $1.5\times 10^{5}$ | $9.7\times 10^{-13}$ | NGC6426 | $8.6\times 10^{5}$ | $8.6\times 10^{-13}$ | NGC6440 (*) | $4.6\times 10^{4}$ | $1.1\times 10^{-13}$
NGC6441 (*) | $5.7\times 10^{4}$ | $2.6\times 10^{-14}$ | NGC6453 | $1.9\times 10^{5}$ | $5.2\times 10^{-11}$ | NGC6496 | $6.3\times 10^{5}$ | $2.7\times 10^{-13}$
NGC6517 (*) | $8.0\times 10^{4}$ | $3.8\times 10^{-13}$ | NGC6522 (*) | $1.2\times 10^{5}$ | $5.9\times 10^{-13}$ | NGC6528 | $1.9\times 10^{5}$ | $5.7\times 10^{-11}$
NGC6535 | $4.8\times 10^{5}$ | $2.9\times 10^{-12}$ | NGC6539 | $2.6\times 10^{5}$ | $2.5\times 10^{-12}$ | NGC6540 | $6.8\times 10^{5}$ | $2.0\times 10^{-13}$
NGC6541 (*) | $1.7\times 10^{5}$ | $6.1\times 10^{-13}$ | NGC6544 (*) | $1.0\times 10^{5}$ | $6.4\times 10^{-13}$ | NGC6553 | $1.9\times 10^{5}$ | $7.8\times 10^{-11}$
NGC6558 (*) | $1.4\times 10^{5}$ | $7.0\times 10^{-13}$ | NGC6569 (*) | $1.6\times 10^{5}$ | $1.5\times 10^{-13}$ | NGC6584 | $4.0\times 10^{5}$ | $4.7\times 10^{-13}$
NGC6624 (*) | $1.8\times 10^{5}$ | $4.3\times 10^{-13}$ | NGC6626 (*) | $6.4\times 10^{4}$ | $2.5\times 10^{-12}$ | NGC6637 | $1.8\times 10^{5}$ | $5.4\times 10^{-11}$
NGC6638 (*) | $9.7\times 10^{4}$ | $7.8\times 10^{-13}$ | NGC6642 (*) | $1.0\times 10^{5}$ | $1.3\times 10^{-12}$ | NGC6652 (*) | $1.3\times 10^{5}$ | $3.3\times 10^{-13}$
NGC6656 (*) | $1.8\times 10^{5}$ | $1.5\times 10^{-13}$ | NGC6681 (*) | $1.5\times 10^{5}$ | $4.3\times 10^{-13}$ | NGC6712 | $1.9\times 10^{5}$ | $4.5\times 10^{-11}$
NGC6715 (*) | $9.1\times 10^{4}$ | $1.4\times 10^{-12}$ | NGC6717 | $4.7\times 10^{5}$ | $7.4\times 10^{-13}$ | NGC6723 | $2.8\times 10^{5}$ | $2.1\times 10^{-11}$
NGC6749 | $4.2\times 10^{5}$ | $7.0\times 10^{-13}$ | NGC6752 | $2.4\times 10^{5}$ | $2.9\times 10^{-11}$ | NGC6760 | $2.4\times 10^{5}$ | $6.7\times 10^{-12}$
NGC6779 | $2.3\times 10^{5}$ | $1.3\times 10^{-11}$ | NGC6809 | $4.3\times 10^{5}$ | $1.7\times 10^{-12}$ | NGC6838 | $6.7\times 10^{5}$ | $7.0\times 10^{-14}$
NGC6864 (*) | $8.6\times 10^{4}$ | $9.1\times 10^{-13}$ | NGC6934 | $3.2\times 10^{5}$ | $6.4\times 10^{-12}$ | NGC6981 | $5.7\times 10^{5}$ | $7.0\times 10^{-13}$
NGC7006 | $5.1\times 10^{5}$ | $9.1\times 10^{-13}$ | NGC7078 (*) | $1.1\times 10^{5}$ | $4.1\times 10^{-13}$ | NGC7089 (*) | $1.4\times 10^{5}$ | $1.3\times 10^{-13}$
NGC7099 | $3.0\times 10^{5}$ | $5.2\times 10^{-14}$ | NGC7492 | $1.9\times 10^{6}$ | $6.0\times 10^{-13}$ | Pal1 | $2.1\times 10^{6}$ | $3.0\times 10^{-13}$
Pal10 | $4.0\times 10^{5}$ | $4.6\times 10^{-13}$ | Pal11 | $2.0\times 10^{6}$ | $3.5\times 10^{-13}$ | Pal12 | $4.4\times 10^{6}$ | $1.1\times 10^{-12}$
Pal13 | $1.3\times 10^{7}$ | $3.0\times 10^{-13}$ | Pal14 | $1.7\times 10^{7}$ | $7.4\times 10^{-13}$ | Pal15 | $6.3\times 10^{6}$ | $9.2\times 10^{-14}$
Pal2 | $4.9\times 10^{5}$ | $4.1\times 10^{-13}$ | Pal3 | $1.1\times 10^{7}$ | $8.1\times 10^{-13}$ | Pal4 | $8.8\times 10^{6}$ | $6.6\times 10^{-13}$
Pal5 | $1.5\times 10^{7}$ | $6.2\times 10^{-13}$ | Pal6 (*) | $1.6\times 10^{5}$ | $1.0\times 10^{-12}$ | Pal8 | $5.6\times 10^{5}$ | $1.1\times 10^{-12}$
Pyxis | $7.2\times 10^{6}$ | $6.0\times 10^{-13}$ | Rup106 | $2.2\times 10^{6}$ | $9.4\times 10^{-13}$ | SagittariusII | $1.9\times 10^{7}$ | $7.3\times 10^{-13}$
Ter1 (*) | $8.3\times 10^{4}$ | $7.1\times 10^{-13}$ | Ter10 | $1.8\times 10^{5}$ | $7.3\times 10^{-11}$ | Ter12 | $2.1\times 10^{5}$ | $3.1\times 10^{-11}$
Ter2 | $2.4\times 10^{5}$ | $1.2\times 10^{-11}$ | Ter3 | $9.8\times 10^{5}$ | $1.2\times 10^{-12}$ | Ter4 | $3.4\times 10^{5}$ | $5.9\times 10^{-12}$
Ter5 (*) | $7.7\times 10^{4}$ | $3.5\times 10^{-13}$ | Ter6 (*) | $4.9\times 10^{4}$ | $1.5\times 10^{-12}$ | Ter7 | $3.2\times 10^{6}$ | $7.4\times 10^{-13}$
Ter8 | $4.1\times 10^{6}$ | $9.4\times 10^{-13}$ | Ter9 (*) | $7.7\times 10^{4}$ | $4.7\times 10^{-13}$ | Ton2 | $3.8\times 10^{5}$ | $1.1\times 10^{-13}$
UKS1 | $2.8\times 10^{5}$ | $2.0\times 10^{-12}$ | VVV-CL001 (*) | $1.4\times 10^{5}$ | $1.1\times 10^{-12}$ | Whiting1 | $1.4\times 10^{7}$ | $1.6\times 10^{-13}$

To choose the optimal time-step $\Delta t$ for the simulations, we quantified
the energy conservation of the orbit integration, in the case where the 159
clusters evolve in isolation, that is in the case where each test-particle in
a cluster feels the gravitational attraction of the cluster itself, but not
that of the Galaxy. Since in the case of an isolated cluster, the
gravitational potential is time-independent, the total energy $E_{i}$ of each
particle in the system (sum of the particle kinetic and potential energy) must
be conserved. Any departure from energy conservation is thus a check of the
quality of the integration and in particular of the choice of the adopted
time-step.

For the isolated simulations, we thus proceeded as follows. We modeled each
cluster as a set of $N$=100 000 test-particles, subject to the cluster
potential, only. To model this latter, for each cluster we adopted the same
Plummer sphere distribution, with same characteristic radius $r_{c}$ and total
masse $M_{GC}$, as those adopted for the simulations described in Sect. 2 (see
in particular Sect. 2.2). A first set of 159 isolated simulations was run
adopting a $\Delta t=10^{5}$ yr and a total number of steps $N_{steps}=50000$
for all clusters, for a total simulated time interval of 5 Gyr. We then
quantified the energy conservation by calculating the median error per time-
step, $m_{err}$, of $(\Delta E/E)_{i}=|(E_{fin,i}-E_{ini,i})/E_{ini,i}|$,
where $E_{ini,i}$ and $E_{fin,i}$ are, respectively, the initial (time $t=0$)
and final (time t=5 Gyr) energy of each particle in the system. While for most
of the simulated clusters (109 over 159), this choice of the time-step was
sufficient to guarantee an excellent energy conservation (with $m_{err}$
typically of the order of $10^{-11}-10^{-12}$), for 50 clusters the
corresponding $m_{err}$ values were found to be above $10^{-10}$. This was the
case for all clusters with crossing times
$t_{cross}=\sqrt{{r_{c}}^{3}/GM_{GC}}$ lower than $2\times 10^{5}$ yr, that is
about twice the time-step. For these clusters, we hence reduced the $\Delta t$
of a factor 10, rerunning the simulations with $\Delta t=10^{4}$ yr and
$N_{steps}=500000$. With such a choice, the corresponding energy conservation
turned out to be excellent (below $10^{-10}$ per step). In Table 4, we
summarize the result of this study, reporting the cluster name, crossing time,
and median error, $m_{err}$, in energy conservation obtained for all isolated
cluster simulations. Clusters for which a $\Delta t=10^{4}$ yr, and a
corresponding number $N_{steps}$ have been used, are indicated in the Table
with an asterisk. The values of $\Delta t$ adopted for the isolated
simulations, and the associated number of time steps, $N_{steps}$, are also
those used to run the simulations of the same clusters orbiting in the
gravitational field of the Milky Way.

##  Appendix C Extra-tidal features generated by all the simulated clusters

In this section, we report all the extra-tidal features as predicted by our
models. For each cluster, we show (from Figs. 14 to 33) the probability
density of finding associated extra-tidal features in the sky, by calculating
the 2D histogram of the escaped particles. Each particle per cluster is
present, that is all 100,000 particles originating each of the 50 Monte-Carlo
realizations plus the case with the best values (see Sect. 2.1). Thus, each
map is a 500x500 histogram that bins 5.1$\times 10^{6}$ particles. The
retrieved cumulative 2D histogram is then normalized to its maximum value and
plotted in the following figures, in logarithmic scale.

Figure 14: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 15: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 16: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 17: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 18: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 19: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 20: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 21: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 22: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 23: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 24: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 25: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 26: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 27: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 28: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 29: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 30: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 31: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 32: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

Figure 33: Projected density distribution in the $(\ell,b)$ plane of a subset
of simulated globular clusters, as indicated at the top of each panel. In each
panel, the red circle indicates the current position of the cluster. The
densities have been normalized to their maximum value.

##  Appendix D Globular clusters classification

Fig. 34 (top panel) shows the distribution of the arctangent of the
$z_{max}/R_{max}$ ratio, with $z_{max}$ and $R_{max}$ being respectively the
maximum height above or below the Galactic plane and their maximum in-plane
distance from the Galactic center, reached in the past 5 Gyr of orbital
evolution in the Galactic potential adopted in this paper (see Sect. 2.3). As
already noticed for field stars (see Haywood et al. 2018), also the GCs
distribution shows a dip at about $10^{\circ}$, which separates clusters with
flattened orbits (arctan($z_{max}/R_{max}$) $\leq 10^{\circ}$) from thicker
ones. We thus define a first set of clusters (the _disk GCs_) as that
containing all globular clusters with arctan($z_{max}/R_{max}$) $\leq
10^{\circ}$. This first set contains 21 clusters. Of the remaining 138, we
distinguish between a _inner GCs_ sample, and a _outer GCs_ sample, on the
basis of the maximum 3D distance ($r_{max}$) that the cluster reaches from the
Galactic center. Inner GCs are those with $r_{max}\leq 10$ kpc and outer GCs
are those with $r_{max}>8.34$ kpc, which is the value of the distance of the
Sun to the Galactic Center used in this experiment. Such a value allows to
discriminate between two classes of tidal debris, for inner clusters are
necessarily restricted in latitude and longitude whereas outer cluster can
fill the sky.

Finally the third and bottom panels of Fig. 34 show the distribution of these
three defined groups in the $R_{max}-z_{max}$ and $E-L_{z}$ planes. We note
that, since disk clusters are uniquely defined on the basis of the ratio
between the maximum vertical and in-plane orbital excursion, and not on the
circularity of their orbits (as done, for example, by Massari et al. 2019),
some of our disk GCs have elongated (i.e. radial) orbits ($L_{z}\sim 0$) or
even retrograde ones ($L_{z}>0$). Our definition of disk clusters is purely
related to a morphological criterium: disk GCs are those whose orbits are
confined close to the Galactic plane, independently on their eccentricity.

Figure 34:  _From the top left to the bottom right:_ (_First panel_) Distribution of the arctangent of the $z_{max}/R_{max}$ ratio for all simulated GCs. The values are expressed in degrees. The vertical dashed line at $10^{\circ}$ separates disk clusters (arctan($z_{max}/R_{max}$) $\leq 10^{\circ}$) from the rest of the population of GCs; (_Second panel_) Distribution of the maximum 3D distance, $r_{max}$, from the Galactic center, reached by the GCs orbits in the last 5 Gyr. The main plot shows this distribution for $r_{max}\leq 20$ kpc, where as the inset shows the whole distribution, which extends at $r_{max}>100$ kpc. In both panels the vertical dashed line, at the solar radius r⊙, separates the group of inner GCs from the group of outer GCs. Note that the clusters in one of these two groups which also satisfy the criterion to be disk clusters are classified as disk GCs, and are not in the inner/outer GCs groups. (_Third panel_) Distribution of disk GCs (magenta points), inner GCs (turquoise points), and outer GCs (dark turquoise points) in the $R_{max}-z_{max}$ plane. The main panel shows the distribution of the GCs having $R_{max}\leq 20$ kpc, the inset the distribution of the whole GC sample. (_Fourth panel_) Distribution of the disk GCs (magenta points), inner GCs (turquoise points), and outer GCs (dark turquoise points) in the $E-L_{z}$ plane. The dashed grey lines correspond, for any given energy E, to the angular momentum of the corresponding circular orbit. Prograde orbits correspond to negative $L_{z}$ values, retrograde orbits to positive $L_{z}$.  Table 5: Classification of the 159 globular clusters studied in this paper in disk clusters (”D”), inner clusters (”I”) and outer clusters (”O”). The values of $r_{max}$ and the angle in degrees of the $arctan(z_{max}/R_{max})$ are also given. Cluster | $r_{max}$ | angle | class | Cluster | $r_{max}$ | angle | class | Cluster | $r_{max}$ | angle | class
---|---|---|---|---|---|---|---|---|---|---|---
2MASS-GC01 | 5.45 | 0.85 | D | 2MASS-GC02 | 7.38 | 17.15 | I | AM1 | 123.19 | 36.56 | O
AM4 | 26.97 | 43.53 | O | Arp2 | 42.44 | 44.54 | O | BH140 | 10.52 | 5.54 | D
BH261 | 3.83 | 21.02 | I | Crater | 147.37 | 44.46 | O | Djor1 | 11.74 | 5.77 | D
Djor2 | 2.39 | 16.74 | I | E3 | 12.91 | 24.07 | O | ESO280-SC06 | 13.82 | 33.04 | O
ESO452-SC11 | 2.53 | 41.47 | I | Eridanus | 109.07 | 44.83 | O | FSR1716 | 5.46 | 17.07 | I
FSR1735 | 4.43 | 17.10 | I | FSR1758 | 14.35 | 23.13 | O | HP1 | 2.70 | 53.56 | I
IC1257 | 19.69 | 19.93 | O | IC1276 | 7.98 | 5.71 | D | IC4499 | 27.17 | 42.77 | O
Laevens3 | 70.74 | 41.35 | O | Liller1 | 1.20 | 5.84 | D | Lynga7 | 4.69 | 17.18 | I
NGC104 | 7.71 | 24.70 | I | NGC1261 | 20.99 | 36.68 | O | NGC1851 | 19.94 | 29.46 | O
NGC1904 | 19.57 | 22.81 | O | NGC2298 | 16.69 | 26.93 | O | NGC2419 | 96.76 | 38.13 | O
NGC2808 | 14.92 | 14.57 | O | NGC288 | 12.50 | 39.00 | O | NGC3201 | 25.53 | 22.62 | O
NGC362 | 12.27 | 36.03 | O | NGC4147 | 24.90 | 45.18 | O | NGC4372 | 7.36 | 16.30 | I
NGC4590 | 27.99 | 32.55 | O | NGC4833 | 8.13 | 9.07 | D | NGC5024 | 22.32 | 44.24 | O
NGC5053 | 18.08 | 44.39 | O | NGC5139 | 7.14 | 21.89 | I | NGC5272 | 16.02 | 40.27 | O
NGC5286 | 13.32 | 21.11 | O | NGC5466 | 41.16 | 43.51 | O | NGC5634 | 22.18 | 43.13 | O
NGC5694 | 51.29 | 35.50 | O | NGC5824 | 32.44 | 40.00 | O | NGC5897 | 9.18 | 41.61 | O
NGC5904 | 24.76 | 42.83 | O | NGC5927 | 5.83 | 7.91 | D | NGC5946 | 5.31 | 23.80 | I
NGC5986 | 5.00 | 29.36 | I | NGC6093 | 4.04 | 49.12 | I | NGC6101 | 32.24 | 30.56 | O
NGC6121 | 6.84 | 4.48 | D | NGC6139 | 3.72 | 34.95 | I | NGC6144 | 4.41 | 43.66 | I
NGC6171 | 3.95 | 33.40 | I | NGC6205 | 8.96 | 39.68 | O | NGC6218 | 4.95 | 31.01 | I
NGC6229 | 30.26 | 37.44 | O | NGC6235 | 8.37 | 34.55 | O | NGC6254 | 4.98 | 29.32 | I
NGC6256 | 2.68 | 14.83 | I | NGC6266 | 2.57 | 21.99 | I | NGC6273 | 5.56 | 37.78 | I
NGC6284 | 6.51 | 36.21 | I | NGC6287 | 6.50 | 30.33 | I | NGC6293 | 3.40 | 37.21 | I
NGC6304 | 3.38 | 14.75 | I | NGC6316 | 3.80 | 22.98 | I | NGC6325 | 2.57 | 32.34 | I
NGC6333 | 9.07 | 28.89 | O | NGC6341 | 10.90 | 40.72 | O | NGC6342 | 2.48 | 38.09 | I
NGC6352 | 4.53 | 9.72 | D | NGC6355 | 3.55 | 29.38 | I | NGC6356 | 8.83 | 28.77 | O
NGC6362 | 5.41 | 33.16 | I | NGC6366 | 6.04 | 17.80 | I | NGC6380 | 2.35 | 16.74 | I
NGC6388 | 3.91 | 19.73 | I | NGC6397 | 6.61 | 27.48 | I | NGC6401 | 3.70 | 21.02 | I
NGC6402 | 3.99 | 32.64 | I | NGC6426 | 16.84 | 22.68 | O | NGC6440 | 1.53 | 25.31 | I
NGC6441 | 4.67 | 16.48 | I | NGC6453 | 2.71 | 36.68 | I | NGC6496 | 5.71 | 26.65 | I
NGC6517 | 3.31 | 21.34 | I | NGC6522 | 1.97 | 42.54 | I | NGC6528 | 2.89 | 16.38 | I
NGC6535 | 4.92 | 16.72 | I | NGC6539 | 3.64 | 39.67 | I | NGC6540 | 2.54 | 12.61 | I
NGC6541 | 4.78 | 29.39 | I | NGC6544 | 5.93 | 18.30 | I | NGC6553 | 4.36 | 4.14 | D
NGC6558 | 2.75 | 21.35 | I | NGC6569 | 2.85 | 26.50 | I | NGC6584 | 20.31 | 35.07 | O
NGC6624 | 1.61 | 64.40 | I | NGC6626 | 3.22 | 20.88 | I | NGC6637 | 2.11 | 56.81 | I
NGC6638 | 2.34 | 31.98 | I | NGC6642 | 2.20 | 26.61 | I | NGC6652 | 3.15 | 49.58 | I
NGC6656 | 10.87 | 21.01 | O | NGC6681 | 6.33 | 41.70 | I | NGC6712 | 5.15 | 28.49 | I
NGC6715 | 38.71 | 45.11 | O | NGC6717 | 2.48 | 33.11 | I | NGC6723 | 4.26 | 47.80 | I
NGC6749 | 5.05 | 3.24 | D | NGC6752 | 5.72 | 20.72 | I | NGC6760 | 5.95 | 6.19 | D
NGC6779 | 13.46 | 30.72 | O | NGC6809 | 6.50 | 36.15 | I | NGC6838 | 7.34 | 5.82 | D
NGC6864 | 16.06 | 33.89 | O | NGC6934 | 37.41 | 20.95 | O | NGC6981 | 21.54 | 35.97 | O
NGC7006 | 47.26 | 32.50 | O | NGC7078 | 10.86 | 25.72 | O | NGC7089 | 18.77 | 34.91 | O
NGC7099 | 8.76 | 40.99 | O | NGC7492 | 25.78 | 45.99 | O | Pal1 | 18.77 | 13.96 | O
Pal10 | 12.06 | 4.56 | D | Pal11 | 8.69 | 23.83 | O | Pal12 | 41.33 | 42.80 | O
Pal13 | 49.09 | 42.54 | O | Pal14 | 88.95 | 35.97 | O | Pal15 | 46.62 | 44.52 | O
Pal2 | 38.20 | 8.94 | D | Pal3 | 104.49 | 43.88 | O | Pal4 | 105.89 | 43.58 | O
Pal5 | 17.54 | 42.94 | O | Pal6 | 3.53 | 27.77 | I | Pal8 | 3.99 | 20.14 | I
Pyxis | 73.65 | 47.40 | O | Rup106 | 32.08 | 31.98 | O | SagittariusII | 75.32 | 41.80 | O
Ter1 | 2.92 | 2.13 | D | Ter10 | 6.55 | 32.40 | I | Ter12 | 4.23 | 17.31 | I
Ter2 | 1.49 | 18.30 | I | Ter3 | 3.81 | 26.97 | I | Ter4 | 1.11 | 28.72 | I
Ter5 | 2.04 | 7.10 | D | Ter6 | 1.93 | 8.77 | D | Ter7 | 41.39 | 45.10 | O
Ter8 | 47.16 | 45.55 | O | Ter9 | 2.92 | 7.94 | D | Ton2 | 4.39 | 24.81 | I
UKS1 | 7.95 | 4.51 | D | VVV-CL001 | 5.04 | 7.09 | D | Whiting1 | 52.30 | 43.09 | O


[1]https://doi.org/#1

.tiffpng.pngconvert #1 .tiff

A deep learning approach for lower back-pain risk prediction during manual
lifting

Kristian Snyder1*, Brennan Thomas1, Ming-Lun Lu2, Rashmi Jha1, Menekse S.
Barim2‡, Marie Hayden2‡, Dwight Werren2‡

1 Department of Electrical Engineering and Computing Systems, University of
Cincinnati, Cincinnati, Ohio, United States of America

2 National Institute for Occupational Safety and Health, Cincinnati, Ohio,
United States of America

These authors contributed equally to this work.

‡These authors also contributed equally to this work.

¤Current Address: Department of Engineering and Computer Science, University
of Cincinnati, Cincinnati, Ohio, United States of America

*<EMAIL_ADDRESS>
## Abstract

Occupationally-induced back pain is a leading cause of reduced productivity in
industry. Detecting when a worker is lifting incorrectly and at increased risk
of back injury presents significant possible benefits. These include increased
quality of life for the worker due to lower rates of back injury and fewer
workers’ compensation claims and missed time for the employer. However,
recognizing lifting risk provides a challenge due to typically small datasets
and subtle underlying features in accelerometer and gyroscope data. A novel
method to classify a lifting dataset using a 2D convolutional neural network
(CNN) and no manual feature extraction is proposed in this paper; the dataset
consisted of 10 subjects lifting at various relative distances from the body
with 720 total trials. The proposed deep CNN displayed greater accuracy
(90.6%) compared to an alternative CNN and multilayer perceptron (MLP). A deep
CNN could be adapted to classify many other activities that traditionally pose
greater challenges in industrial environments due to their size and
complexity.

## Introduction

Back pain, especially when occupationally-induced, is an extremely common
ailment. 23.2% of the world’s population is estimated to be affected in any
given month [1]. It is considered the leading cause of job-related disability
and missed work days, resulting in massive losses of productivity [2, 3]. Back
pain is, in the United States, the largest contributor to total workers’
compensation costs as of 2015, representing over 20% of costs and 13.7 billion
USD annually [4].

The revised National Institute for Occupational Safety and Health (NIOSH)
lifting equation [5] (RNLE) is currently considered the leading method of
measuring back pain risk involved in both single and repeated lifting of
objects in the workplace [6, 7]. Given the mass, relative source distance,
relative destination distance (to the person) of an object, frequency of the
lifting tasks, and the work-rest pattern, the RNLE determines the relative
level of risk to the lifter. Taking measurements for using the RNLE in the
field presents a challenge because the analyst needs to interrupt work
activity for measuring several lifting variables to recognize the
characteristics of lifting tasks. In a world of activity recognition and
wearable sensors becoming prominent [8, 9, 10, 11], there is a distinct lack
of research into automatically detecting lift risk, especially as it relates
to back pain.

Real-time recognition of unsafe lifting would provide immediate feedback to
the users, something currently impractical with traditional methods such as a
specialist constantly monitoring the workers. Therefore, the ability to, in
real-time, detect unsafe lifting behavior in an industrial setting would
provide significant benefits. Closing the feedback loop allows for the wearer
to be quickly alerted to any risk and help prevent further stress. A response
time in seconds presents a significant improvement over typical pain feedback,
which does not always present itself quickly and may be alongside debilitating
injury [12]. Additionally, an automated approach is far more scalable and can
be inexpensively rolled out to an entire set of workers for minimal cost
compared to treatment and loss of productivity due to back injury.

Classifying a person lifting objects of various distance from themselves
presents an inherent challenge because even to an observer, the different
movements are quite similar. Additionally, sourcing data for specialized
activities has its own challenges, given a dearth of existing datasets and the
increased expense of independently collecting data. Consequently, datasets are
typically small; the NIOSH lifting dataset consists of 720 total trials [13].
The challenge, then, is to develop a versatile model that can distinguish
between very similar activities and operates on datasets possibly magnitudes
smaller than for similar problems [14].

Significant advancements in deep learning approaches have been achieved in
recent years, broadly categorized into video-based and sensor-based
classification [15]. Video-based approaches have focused mostly on
surveillance and gait recognition, although there is recent research that has
successfully labeled live video data with a high degree of accuracy without
assistance from manual body part labeling or tracking devices [15, 16]. While
impressive, use of video classification in many workplaces requires a wide
array of surveillance cameras placed to view all workers, which can be
expensive and require complicated installation. In other areas, such as most
construction and industrial sites, video is near-impossible due to the
frequently temporary and constantly fluctuating workplaces. Sensors placed on
the body do not require the person to be in a specific place and are
relatively inexpensive and less privacy-intrusive than constant video capture
[15], making them more practical for this use case.

Leading sensor-based approaches typically classify more distinct activities,
such as standing, walking, and running with relative success [17, 18, 14].
However, classification of visually similar activities (such as walking up and
down stairs) displays notably lower performance [17]. Additionally, other
attempts to classify accelerometer data focused on $n>500$ samples for each
class and commonly in excess of thousands of samples [10, 17, 14, 19, 20]. The
existence of other datasets with multiple accelerometers suggests a place for
a model that can work on very detailed activities generating much more data
with dedicated sensors than a smartwatch or smartphone [10].

In this paper, a model utilized to classify the above lifting activities is
presented as a generalizable solution to the problem of classifying small
datasets of very similar activity classes. The model uses 2-dimensional
convolutional layers along with average pooling to classify the activities.
Accelerometer and gyroscope data is first preprocessed with a Butterworth
filter and manipulated into a 3D matrix resembling an image before being
trained on. The structure of the model is examined by comparing it to
variations in pooling, regularization, and complexity to display the
theoretical underpinnings, increasing its adaptability to other problems. The
network’s various hyperparameters are also examined to justify their specific
values. Finally, accuracy and other statistics of the proposed model are
compared to traditional activity recognition models and their published
results.

Expanding the list of activities to a wider subset of possible movement
unlocks a multitude of possible benefits. In collaboration with the National
Institute for Occupational Safety and Health (NIOSH), a deep convolutional
neural network (CNN) model to classify relative risk level of lifting objects
is developed. Such a model would assist in preventing serious, chronic back
injury in workers and significantly improve their quality of life.

## Materials and methods

### Data collection

For purposes of analysis, data utilized for model development was sourced from
a previous study by researchers at NIOSH to examine body motion for two-handed
lifting tasks similar to those in the workplace [13]. Lifts were performed
using the American Conference of Governmental Industrial Hygienists (ACGIH)
Threshold Limit Values (TLV) for lifting, which defines 12 zones relative to
the body in the sagittal plane shown in Fig 1. Save for Zones 1-3, 4, 7, and
10, all lifts began in the midpoint of these zones. In the aforementioned
exceptions, initial points in the zone were altered to provide realistic
motion for the subjects’ ranges of motion. The object lifted consisted of a
$36\text{cm}\times 12$ cm wire grid weighing 0.45 kg with two handles,
simulating a crate or box. To prevent injury, the weight was kept small. The
order of lifting from each zone was randomized to prevent bias from ordering
the lifts. A total of 720 trials with 6 trials for each subject in all zones
were available as input data.

Fig 1: ACGIH lifting zone system depicting the relative areas collected for
analysis. (Source: NIOSH)

Five male and five female subjects (mean and SD: $170\pm 7.4$ cm for height
and $85.7\pm 20.2$ kg for weight) participated in the lifting process. All
participants had six inertial measurement unit (IMU) sensors (Kinetic Inc.)
attached to their bodies on the upper back (T12), each wrist, the dominant
upper arm, waist, and the dominant thigh during all lifting. Each sensor
consisted of a tri-axial gyroscope and accelerometer sampling at a rate of 25
Hz. All sensor data was calibrated and synchronized prior to data collection
[13]. An example of this data is shown in Fig 2.

Fig 2: Plots of the accelerometer and gyroscope sensors for Subject 1’s first
lift in Zone 1 (high risk) from [13]. Not all data collected is shown; the two
vertical black lines show the beginning and end of the actual lift.

### Feature extraction

In training the model, 540 of the 720 total trials were utilized for training
while the remaining 25% (180 trials) formed the test set. Training and testing
subsets were sampled randomly without replacement. Features were the
preprocessed lift, with zero-padding applied to lifts that did not reach the
full time. All training utilized all of the 6 sensor areas, each containing an
accelerometer and gyroscope, resulting in 18 total measurements for each point
in time. The twelve ACGIH lifting zones were mapped to three risk levels: low,
medium, and high-risk. The mapping is based on a slight modification to the
Los Alamos National Laboratory recommendations to simplify the ACGIH zones;
the zones mapped are listed in Table 1.

Table 1:  Mapping of ACGIH lifting zones to relative risk levels. Risk Level | Zones
---|---
Low | 4, 5
Medium | 6, 7, 8, 9
High | 1, 2, 3, 10, 11, 12

Finally, in each case the data was scaled with normalization to the bounds
[-1, 1] based on the training data, with the same scaling applied to the
testing data. This is in direct interest of increasing performance of the
neural network, which trains best on data normalized to these bounds. This
resulted in 720 total 27,000-dimension class-labelled vectors for training and
testing.

A maximum time window of 30 seconds was selected to train the model. Lifts
that did not reach this length of time were zero-padded to reach the full
time. Most trials did not reach this time period; the majority ended between
10 and 15 seconds. The long time slice length was selected to investigate
dependence on alignment of the starting and ending times, ideally mitigating
or eliminating any significant dependence. To reduce overall noise and drift,
a Butterworth filter with order 2, lower bound of 2 Hz and upper bound of 12
Hz was applied to each dimension (X, Y, Z) of the gyroscope and accelerometer
in each sensor, resulting in 36 total measurements for each point in time. A
Butterworth filter is a signal processing filter that has maximally flat
frequency response in the passband, preserving the original signal better than
other filters [21]. Each dimension of each of the 12 sensors (an acceleromter
and gyroscope each on the side, left wrist, right wrist, back, upper arm, and
thigh) has 750 data points in each trial, resulting in 27,000 total features
as shown in Eq 1.

$\centering\begin{split}&\quad\text{6 IMU devices}\\\ \times&\quad\text{2
sensors per IMU (gyro and accelerometer)}\\\ \times&\quad\text{3 dimensions
each}\\\ \times&\quad\text{30 sec max per trial}\\\ \times&\quad\frac{\text{25
frames}}{\text{sec}}\\\ =&\quad\text{27,000
features}\end{split}\@add@centering$ (1)

Parameters and the filter itself were chosen based on previous research [22],
[23], [8].

To prepare the data for ingestion to the model, each 27,000-feature vector was
reshaped to form a $95\times 95\times 3$ matrix. Specific reshaping is
performed by stacking each sensor to form a $36\times 750\times 3$ matrix,
placing each point of time into successive columns. The matrix is then line-
wrapped to form the final image. See Fig 3 for a visualization of the process.
Fig 4 displays an example of the final image after standardization.

Fig 3: Process of ingesting data into the model for training. All figures not
to scale. The resultant matrix could theoretically be any size; a square was
selected for highest compatibility with existing CNN research. Fig 4: Example
of an input image to the network. The image shown has had a Butterworth
bandpass filter of order 2 and bounds 2 and 12 Hz applied to it in addition to
a standardizing scaler. The grey block at the bottom represents padding to the
model that makes all inputs the same size.

This method produces a data format friendly to CNN models while preserving
time locality of the data as much as possible. Convolutions, therefore, will
be made more often between features occurring at similar times to help correct
for an imprecise window mislabeling the start and end of the lift. Finally,
the data is standardized for each sensor by scaling the data to a mean of zero
and standard deviation of one.

### Overall statistics

During testing and training of the proposed model, several class-specific
statistics were collected to help measure its performance. In addition to the
statistics below, two overall statistics were also collected: RK (primarily
utilized in hyperparameter tuning, see 3.1. Hyperparameter tuning) and overall
accuracy.

Precision (shown in Eq 2) is defined as the proportion of instances belonging
to the class (true positive or TP) over all instances (both TP and false
positive or FP) classified as that class.

$\displaystyle Precision=\frac{TP}{TP+FP}$ (2)

Recall (shown in Eq 3) is defined as the proportion of a elements in a
particular class classified as that class over all elements belonging to that
particular class (including TP and false negative or FN).

$\displaystyle Recall=\frac{TP}{TP+FN}$ (3)

F-measure (shown in Eq 4) is the harmonic mean of precision and recall, used
as an alternative to raw accuracy.

$\displaystyle\textit{F-measure}=\frac{2\times Recall\times
Precision}{Recall+Precision}$ (4)

## Model design

The proposed model is based on the Visual Geometry Group Network (VGGNet,
developed at the University of Oxford), a high-performing CNN model that is
notable for its high depth and use of additional layers and small kernel size
instead of fewer layers with a larger kernel size [24]. This helps to reduce
the number of parameters of the network as well, especially important for
small datasets which most models with many parameters struggle to converge on
[25].

Specifically, the model is based on variation B of VGGNet (VGGNet B), with max
pooling layers separated by groups of two convolutional layers with increasing
filter count. Most notably, the filter count in each group of convolutional
layers is smaller than for VGGNet B, ranging from 32 to 128 filters instead of
64 to 512. Additionally, the $2\times 2$ max pooling layers are replaced with
$2\times 2$ average pooling layers. Max pooling is employed in most CNN models
to increase contrast and preserve the most important information of an image
while decreasing dimensionality of the input; however, contrast for the images
is already quite high, with important features throughout the sample [26].
Average pooling retains more information from layer to layer because it
incorporates all source pixels in the output compared to max pooling rejecting
all but one of the pixels. When adaping VGGNet B to train on accelerometer
data, this most significantly improved performance on the NIOSH lifting
dataset. Table 2 contains a detailed description of the layers.

Table 2: Detailed specification of the layers involved in the proposed model.
All 2D convolution layers contain a ReLU activation layer. Model configuration
---
Layer type | Parameters
Input (95x95x3 matrix)
2D convolution | 32 filters, 3x3 kernel
Average pooling | 2x2 cell
Dropout | 25%
2D convolution | 64 filters, 3x3 kernel
2D convolution | 64 filters, 3x3 kernel
Average pooling | 2x2 cell
Dropout | 25%
2D convolution | 128 filters, 3x3 kernel
2D convolution | 128 filters, 3x3 kernel
Average pooling | 2x2 cell
Dropout | 25%
Flatten
Fully connected | 1024 units
Batch normalization
Dropout | 25%
Softmax output

In addition to the proposed model, a separate model (CNN+LSTM) was developed
as an alternative approach, utilizing a network of 1-dimensional convolutional
layers and long short-term memory (LSTM) layers, based on DeepConvLSTM by
Ordóñez et al [27]. This approach produces a far less complex network and
treats the dataset as a time series instead of a 3D matrix thanks to the LSTM
layers. LSTM layers utilize the current state of the network with their memory
units, building direct relationships between the currently analyzed data and
previously analyzed data. This makes it especially effective on temporally
organized data, such as accelerometer and gyroscope signals [28].

The CNN+LSTM model also uses a slightly data format. Instead of the 97x97x3
matrix, the 12x750x3 matrix is utilized. Features are extracted by sliding a
12x1x3 window (essentially a column of the matrix) over the trial, reducing
the dimensionality of the data from 3D to 2D and making it compatible with the
LSTM layers. See [27] for additional details on the structure of LSTM and 1D
convolutional layers.

Three other models were developed strictly for comparison: a simpler CNN model
with lower depth, a max pooling VGGNet B variation, and a multi-layer
perceptron (MLP) model. The simpler CNN model does not utilize any form of
pooling or a dense layer, consisting only of convolutional layers and a
softmax output layer. The max pooling VGGNet B model is identical to the
proposed model in all ways save the usage of max pooling layers instead of
average pooling.

Models were all trained in the same manner, utilizing ADAM [29] as a gradient
descent estimator, categorical cross-entropy as a loss function and, in place
of a specific number of epochs, utilizing early stopping to halt training. The
early stopping module monitored loss with a min delta of 0 and patience of 10
epochs. After loss failed to decrease, the best-performing (according to loss,
not accuracy) weights of the last 10 epochs were selected for testing.

To develop the models, Keras v2.2.4 was used as an interface to TensorFlow
v1.15.0. Data manipulation and standardization was performed with Scikit-learn
v0.21.3.

### Hyperparameter tuning

A difficult portion when developing machine learning models is tuning the
various hyperparameters, which are parameters for the model that are
statically set before training begins. These parameters can affect training
and testing results as significantly as alterations to the model structure and
so deserve their own discussion.

In developing the model, three hyperparameters were focused on due to their
great effect on accuracy: L2 regularization importance ($\lambda$), dropout
percentage, and learning rate ($\alpha$). Searches were not performed for
global minima due to the extreme effort involved and were tested with both
single-value variation and multiple-value variation to evaluate any effects
dependent on multiple hyperparameters.

$\lambda$ configures the importance of the L2 regularizer, which is applied
solely to the final softmax dense layer in both the activity and output
portions. Regularization in general is the practice of adding a loss function
to the complexity of the layer and incentivizing training of a sparser and
less complex model. This helps to reduce the prevalence of overfitting to the
training data [30]. Importance is the weight given to this regularizer over
the basic loss function used in training; a higher value signifies higher
importance. Most values for $\lambda$ are between 0.1 and $10^{-10}.$

Dropout percentage determines the portion of the weights at each layer that,
for each epoch, are temporarily removed from the model and inserted after
training the epoch. Dropout helps to further reduce overfitting by forcing the
network to learn the same features multiple times as the portions of the model
that previously trained are randomly removed [31]. Increased dropout can
further prevent overfitting but also introduce instability when training,
justifying tuning this parameter.

$\alpha$ configures the degree to which the optimizer tweaks the weights at
each layer to descend the gradient. A higher value results in typically faster
descent but increased instability. Overshooting the target, increasing the
loss, and having to backtrack is more common. Lower $\alpha$ can help to
prevent this, but typically requires a far longer training time and may lead
to the optimizer converging poorly.

When comparing these models performance, a single value was desired to prevent
subjectivity. However, traditional measures like F-measure and accuracy can
fail to sufficiently capture performance, especially for multiclass and
imbalanced datasets as present here. To correct for this and utilize a more
descriptive statistic, RK correlation was utilized [32]. RK is a
generalization of Pearson and Matthews correlation for multiclass problems. It
has been shown to be superior to Cohen’s kappa and standard metrics like
f-measure (the latter specifically due to its inclusion of true negatives) as
a single-metric comparison tool [33].

Eq 5 displays a discretized version of RK that allows for its usage on a
single confusion matrix [32]. $C_{lk}$ represents elements of the matrix.

$\displaystyle
R_{K}=\frac{\sum\limits_{klm}C_{kk}C_{lm}-C_{kl}C_{mk}}{\sqrt{\sum\limits_{k}\sum\limits_{l}C_{kl}\sum\limits_{\begin{subarray}{c}l^{\prime}\\\
k^{\prime}\neq
k\end{subarray}}C_{k^{\prime}l^{\prime}}}\sqrt{\sum\limits_{k}\sum\limits_{l}C_{lk}\sum\limits_{\begin{subarray}{c}l^{\prime}\\\
k^{\prime}\neq k\end{subarray}}C_{l^{\prime}k^{\prime}}}}$ (5)

## Results

720 total trials were used for training and testing. Table 3 lists the number
of trials for each class label.

Table 3:  Number of trials for each class of lift. Class | Count
---|---
Low Risk | 120
Medium Risk | 240
High Risk | 360

### Hyperparameter results

Fig 5 displays the results of training the proposed models with different
values for the stated hyperparameters. Fig 5 displays how RK can report
different results due to its consideration of true negatives compared to
accuracy, making it a more well-rounded single metric for comparing various
models.

Fig 5: Comparison of performance for various hyperparameters set on the
proposed model. A: comparison of RK statistics. B: RK statistic compared with
accuracy. Accuracy ranges from 0 to 1 while RK ranges from -1 to 1. Values
from -1 to 0 are not shown due to no results in that range.

However, when selecting the resulting hyperparameters, a qualitative selection
of the second-best performing parameters was made. This is due to the behavior
of the optimizer at such a high learning rate; it behaved in a more unstable
pattern and terminated at a relatively high categorical cross-entropy compared
to a slightly lower $\alpha$. Fig 6 displays the behavior involved. Therefore,
the parameters utilized for the proposed model and variant with max pooling
are $\lambda=10^{-5},\alpha=10^{-3},\text{dropout percentage}=25\%$.

Fig 6: Loss gradient descent for $\alpha$ of 0.01 and 0.001. Each point
represents the training categorical cross-entropy at the completion of each
epoch.

### Classification results

Although various model structures were used, all relied on the same basic
layer unit: the 2D convolutional layer. Models varied in complexity and depth,
with some containing LSTM layers in an attempt to formally capture time-based
data.

Table 4 compares precision, recall, and f-measure for the proposed model, a
typical simplified CNN model, and a fully-connected multilayer perceptron
network.

Table 4: Summary of classification results for proposed model and alternatives. Risk Level | Proposed Model | CNN+LSTM | VGGNet B (Max Pool) | CNN | MLP
---|---|---|---|---|---
Low F-measure | 0.776 | 0.8 | 0.746 | 0.444 | 0.367
Medium F-measure | 0.862 | 0.783 | 0.692 | 0.635 | 0.500
High F-measure | 0.964 | 0.923 | 0.815 | 0.765 | 0.667
Accuracy | 0.906 | 0.861 | 0.766 | 0.667 | 0.567
RK | 0.839 | 0.77 | 0.65 | 0.53 | 0.28

The modified VGGNet B model performed the best out of all models on medium and
high-risk lifts and only slightly underperformed the max pooling variant in
low-risk lifts. It also significantly outperformed the simplified CNN and
multilayer perceptron in f-measure for all classes and in overall accuracy.
Fig 7 displays the distribution of predictions for the proposed model on the
testing data. The results for the proposed model are displayed in more detail
in Fig 8.

Fig 7: Swarm plot of the testing results by the proposed model. The x-axis
represents the true labeling and the y-axis the model output. The y-axis is
divided into three zones that define the resulting class for each value,
labeled in the top-left of each box. Fig 8: Heatmap plot of the testing
results by the proposed model. Each row has been normalized so that each class
has the same color scale.

## Discussion

Small datasets are notoriously difficult to train on for machine learning
models and that is borne out here, with typical CNN and MLP models barely
outperforming randomness for some guesses and delivering no more than a 2/3
accuracy. The proposed VGGNet B variant significantly outperformed the other
two models and max pooling variation, both in RK (0.862) and accuracy (90.6%).
Notably, it also displayed the best f-measure for every single class as well,
although it improved the least in low-risk classification.

As RK is typically only considered satisfactory above an 0.7 correlation, the
proposed model and CNN+LSTM models were also the only methods that displayed
acceptable performance.

Notably, the data trained on (see Table 3) is an example of an imbalanced
dataset, with significantly varying numbers of samples in each risk level.
This typically poses a challenge for machine learning models that
automatically become biased toward the larger class [25]. This is obvious from
examining the f-measures for all models, which increase with the number of
trials for the class.

While low-risk lifts were classified more poorly than other classes, the
proposed model improves upon all alternatives and helps to avoid issues such
as over and under-sampling, which can lead to worse performance when testing
[25].

### Average pooling

Average pooling, other than parameter changes, is the major departure from
VGGNet B for the model. It resulted in a 5.3 basis point increase in accuracy
with major benefits to both medium and high-risk accuracy. It is theorized
that, in this case, average pooling outperforms max pooling because it passes
more information to the next layer by using all 4 cells in the pool instead of
selecting the highest value. This could increase generalizability of the model
to testing data. However, max pooling is typically the selected model for
successful CNN-based classifiers [24].

These models typically are trained on largely unprocessed image data that may
have low contrast and significant redundancy in a given region, reducing the
information loss of max pooling. Max pooling picking a blue pixel from the sky
in an image is still representative. Average pooling is similar to various
image downsampling methods, albeit simplified in its attempts to preserve
information. As accelerometer data, especially at low frequencies (the dataset
was recorded at 25 Hz), is vulnerable to loss of information, preserving this
information may have resulted in the increase in accuracy. Use of max pooling
would simulate a very rough downsampling of the data, which could clip many of
the important fluctuations in rotation and acceleration from the sensor data.
Average pooling smoothes out this downsampling and, while it still eliminates
information, performs it less severely than max pooling.

### Saliency mapping

Saliency mapping is the process of determining the input features that the
model recognizes as the most significant to the output class. First defined by
Simonyan et. al. in [34], saliency is a multistep process. Given a final score
matrix $S_{c}(I)$ for a class $c$ selected by the model, the linear score can
be represented in Eq 6 as

$\displaystyle S_{c}(I)=w_{c}^{T}I+b_{c}$ (6)

where $I$ is the pixels of the image, $w_{c}$ is the weights for the class,
and $b_{c}$ is the bias of the class. As the model is non-linear in nature due
to the activation functions and so could not be easily computed, the
approximation is computed instead in Eq 8.

$\displaystyle S_{c}(I)\approx w^{T}I+b$ (7)

$w^{T}$ is the derivative of the score matrix $S_{c}$ to the image matrix $I$
at the point $I_{0}$ (the image itself), shown in Eq 8 [34].

$\displaystyle w=\frac{\delta S_{c}}{\delta I}\bigg{\rvert}_{I_{0}}$ (8)

This approximation determines how much each pixel, if changed, would affect
the class score and assigns a value to them.

Here, the primary purpose was to determine whether the model was truly
recognizing portions of the input as the lift or simply tweaking the weights
to fit on noise, which is a possibility. However, this is not the behavior
displayed in the saliency plots obtained from the proposed model.

The low and medium saliency plots (Fig 9 and Fig 9) are nearly identical but
this is not an indicator that they cannot be differentiated. Instead, the
model considers the same parts of the lifting motion as important for each
class. For high-risk lifts, shown in Fig 9, the model is drawing on two very
specific regions, with very low weighting for the rest of the image. We
theorize that it may have determined two specific points of the lift: when the
object is lifted and when the lifter accelerates back to a neutral, upright
position. These regions are also present on the other two saliency plots but
are not as clearly delineated. This difference in quality is the clearest
indicator that further improvements are possible beyond simply examining
classification accuracy. However, this would most likely require alterations
to the structure of the model and not simple hyperparameter tuning.

Fig 9: Saliency plots for final softmax layer of network. A: low-risk
saliency. B: medium-risk saliency. C: high-risk saliency. Bright green/yellow
represents the highest weighting; dark purple represents the lowest weighting.

Producing these plots also provides a separate benefit: examining what sensors
contribute most effectively to the results. While more sensor data can assist
in classification, many of the sensors have difficult or impractical
placement, such as those on the thigh and upper back. On the high-risk
classification, the two hot areas center around the wrist and back sensors.
Side, upper arm, and thigh sensors contributed to the classification but, as
shown in Fig 9, are not as bright and so are candidates for possible removal
in future research.

A significant advantage of the CNN+LSTM model is, due to the data input shape,
saliency analysis is far clearer here. Fig 10 displays the saliency for high-
risk lifts as well, with the various sensors labeled. Here, it is far clearer
that the back and wrist sensors are the most significant. While all sensors
contribute to the prediction (in at least one dimension), removing the thigh,
upper arm, and side sensors may still leave enough information to sufficiently
classify lifting risk level.

Fig 10: Saliency plot for high-risk lift trials obtained from the CNN+LSTM
model. Scale ranges from deep blue as the lowest significance and deep red as
the highest significance. The x-axis is frames of the input data and y-axis is
the sensor data, where A/G is accelerometer or gyroscope and x, y, z are the
dimensions for the sensor.

Additionally, we see two general peaks here as well, suggesting that lifting
behavior may have an initial acceleration and final acceleration as general
features. As both models appear to be examining the same region of the data,
it is increasingly likely that the data is sufficiently separable and contains
true features instead of simple noise.

### Feature extraction

Many other examples of accelerometer classification employ manual feature
extraction with various measurements including means, variances, zero-crossing
rates, and various other statistics [35].

This method only requires a signal filter and standardization before reshaping
and feeding to the model. This could be performed very quickly and benefits
from significant previous efforts in signal processing to allow streamed and
real-time transformations. While many more features are utilized, handcrafted
features typically require significant domain experience and tuning for high
performance. The preprocessing proposed could easily be applied to many
different datasets with minimal alterations.

### Possible improvements and future research

90.6% accuracy, while a significant improvement over other models, still
provides opportunity for further advancement. It is doubtful that
hyperparameter tweaking could significantly increase accuracy, as Fig 5
displays only small increases once 85% accuracy was reached. Model alterations
are most likely necessary to reach overall 95% testing accuracy. However,
high-risk classification, considered most important to ensure worker safety,
is excellent, with 96.9% accuracy.

One possible region of interest is configuring the number of convolutional
layers between the average pooling layers. Deeper into the network, additional
convolutional layers may provide high-level feature extraction. Nearer to the
input layer, more layers typically increase recognition of granular details.
However, as stated in the introduction, a high number of parameters for a
small dataset such as this can result in the optimizer failing to converge.
Therefore, simply adding layers may fail to significantly increase accuracy
and require alterations to the general network structure as well.

The presence of a fully-connected layer at the end of the network is also a
point of interest. This layer especially provides a significant number of
parameters and so is a target for optimization. Altering the number of units
and configuring a regularizer on that layer may assist, especially on
improving low-risk lifts due to optimizing for a less complex model.

For real-world use, minimizing the number of sensors will significantly
advance the practicality, reducing cost and eliminating the awkward placement
of several sensors. Ideally, one sensor would provide sufficient data to
classify lifting. However, the two wrist sensors and a side sensor are simple
enough to attach that they may be an acceptable alternative. The strong
emphasis, unfortunately, by both models on the back sensor may indicate that
this sensor is necessary. Redundancies may nonetheless exist in the sensor
data and may provide sufficient separability for a similar model.

Finally, based on the saliency maps generated, the window size for the lift
could be significantly reduced with, in all likelihood, minimal reduction in
accuracy. This would increase the speed of activity recognition in real-world
use simply due to the reduced input size. A reduction in input size also
decreases the number of parameters in the network by eliminating features,
possibly allowing for further depth in the network. However, overzealous input
size reduction could cut off significant parts of lifts and so needs to be
performed carefully.

## Conclusion

Classifying accelerometer data is traditionally difficult due to the
requirement of most machine learning models requiring large datasets.
Therefore, much of the existing research focuses on typical activity and
exercise classification that can draw on existing datasets or be compiled from
many subjects. Specialized activities with lesser impact, then, have been
neglected due to the difficulty involved in compiling enough data for
traditional models.

The proposed model was able to quickly and accurately (90.6% accuracy, 0.839
RK) classify a small accelerometer dataset provided by NIOSH with minimal
feature extraction and significantly greater performance than other models
tested. Specifically, the usage of a CNN that would normally classify images
along with the alteration to use average pooling over max pooling provided the
greatest benefit. Hyperparameter tuning was also shown to have significant
effects on the performance, but in lower magnitudes. It is very likely that a
similar model could be trained on other small and/or unbalanced datasets to
make their classification feasible where other models have failed. The
proposed technique provides excellent monitoring of risks involved in the
various types of manual lifts in an industrial setup to provide timely
interventions.

## Acknowledgments

This study was funded, in part, by the National Institute for Occupational
Safety and Health (NIOSH) and the Research Experiences for Undergraduates
Program of the National Science Foundation (REU NSF) REU supplement for NSF
award no.: ECCS 1556294. Disclaimer: the findings and conclusions in this
report are those of the authors and do not necessarily represent the official
position of NIOSH, Centers for Disease Control and Prevention (NIOSH/CDC) or
NSF. Mention of any company or product does not constitute endorsement by
NIOSH/CDC or NSF.

## References

  * 1. Hoy D, Bain C, Williams G, March L, Brooks P, Blyth F, et al.  A systematic review of the global prevalence of low back pain.  Arthritis and Rheumatism. 2012;64(6):2028–2037.  {10.1002/art.34347}.
  * 2. Bernhard BP.  A critical review of epidemiological evidence for work-related musculoskeletal disorders of the neck, upper extremity and low back.  Musculoskeletal disorders and workplace factors. 1997;.
  * 3. Council NR, Others.  Musculoskeletal disorders and the workplace: low back and upper extremities.  National Academies Press; 2001.
  * 4. Liberty Mutual Research Institute for Safety.  Liberty Mutual Workplace Safety Index.  Liberty Mutual Insurance; 2018.  Available from: https://business.libertymutualgroup.com/business-insurance/Documents/Services/WorkplaceSafetyIndex.pdf.
  * 5. Waters T, Putz-Anderson V, Garg A.  Applications Manual for the Revised NIOSH Lifting Equation.  DHHS (NIOSH). 1994;94-110:1–164.
  * 6. Waters TR, Lu ML, Piacitelli LA, Werren D, Deddens JA.  Efficacy of the revised NIOSH lifting equation to predict risk of low back pain due to manual lifting: Expanded cross-sectional analysis.  Journal of Occupational and Environmental Medicine. 2011;53(9):1061–1067.  {10.1097/JOM.0b013e31822cfe5e}.
  * 7. Dempsey PG, McGorry RW, Maynard WS.  A survey of tools and methods used by certified professional ergonomists.  Applied Ergonomics. 2005;36(4 SPEC. ISS.):489–503.  {10.1016/j.apergo.2005.01.007}.
  * 8. Music J, Stancic I, Zanchi V.  Is it possible to detect mobile phone user’s attention based on accelerometer measurment of gait pattern?  Proceedings - International Symposium on Computers and Communications. 2013; p. 522–527.  {10.1109/ISCC.2013.6754999}.
  * 9. Liu Y, Nie L, Liu L, Rosenblum DS.  From action to activity: Sensor-based activity recognition.  Neurocomputing. 2016;181:108–115.  {10.1016/j.neucom.2015.08.096}.
  * 10. Bao L, Intille SS.  Activity recognition from user-annotated acceleration data.  Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics). 2004;3001:1–17.  {10.1007/978-3-540-24646-6_1}.
  * 11. Preece SJ, Goulermas JY, Kenney LPJ, Howard D, Meijer K, Crompton R.  Activity identification using body-mounted sensors - A review of classification techniques.  Physiological Measurement. 2009;30(4).  {10.1088/0967-3334/30/4/R01}.
  * 12. Mayer JM, Mooney V, Matheson LN, Erasala GN, Verna JL, Udermann BE, et al.  Continuous Low-Level Heat Wrap Therapy for the Prevention and Early Phase Treatment of Delayed-Onset Muscle Soreness of the Low Back: A Randomized Controlled Trial.  Archives of Physical Medicine and Rehabilitation. 2006;87(10):1310–1317.  {10.1016/j.apmr.2006.07.259}.
  * 13. Barim MS, Lu ML, Feng S, Hughes G, Hayden M, Werren D.  Accuracy of An Algorithm Using Motion Data Of Five Wearable IMU Sensors For Estimating Lifting Duration And Lifting Risk Factors.  Proceedings of the Human Factors and Ergonomics Society Annual Meeting. 2019;63(1):1105–1111.  {10.1177/1071181319631367}.
  * 14. Weiss GM, Lockhart JW, Pulickal TT, McHugh PT, Ronan IH, Timko JL.  Actitracker: A smartphone-based activity recognition system for improving health and well-being.  Proceedings - 3rd IEEE International Conference on Data Science and Advanced Analytics, DSAA 2016. 2016; p. 682–688.  {10.1109/DSAA.2016.89}.
  * 15. Cook D, Feuz KD, Krishnan NC.  Transfer learning for activity recognition: a survey.  Knowledge and Information Systems. 2013;36(3):537–556.  {10.1007/s10115-013-0665-3}.
  * 16. Jalal A, Kim YH, Kim YJ, Kamal S, Kim D.  Robust human activity recognition from depth video using spatiotemporal multi-fused features.  Pattern Recognition. 2017;61:295–308.  {10.1016/j.patcog.2016.08.003}.
  * 17. Kuspa K, Pratkanis T.  Classification of Mobile Device Accelerometer Data for Unique Activity Identification; 2013.  Available from: http://cs229.stanford.edu/proj2013/PratkanisKuspa-ClassificationOfMobileDeviceAccelerometerDataforUniqueActivityIdentification.pdf.
  * 18. Kwapisz JR, Weiss GM, Moore SA.  Activity recognition using cell phone accelerometers.  ACM SIGKDD Explorations Newsletter. 2011;12(2):74.  {10.1145/1964897.1964918}.
  * 19. Maurer U, Smailagic A, Siewiorek DP, Deisher M.  Activity Recognition and Monitoring Using Multiple Sensors on Different Body Positions.  In: International Workshop on Wearable and Implantable Body Sensor Networks (BSN’06). IEEE; 2006. p. 113–116.  Available from: http://ieeexplore.ieee.org/document/1612909/.
  * 20. Hammerla NY, Halloran S, Plötz T.  Deep, convolutional, and recurrent models for human activity recognition using wearables.  IJCAI International Joint Conference on Artificial Intelligence. 2016;2016-Janua:1533–1540.
  * 21. Erer KS.  Adaptive usage of the Butterworth digital filter.  Journal of Biomechanics. 2007;40(13):2934–2943.  {10.1016/j.jbiomech.2007.02.019}.
  * 22. Kwon S, Lee J, Chung GS, Park KS.  Validation of heart rate extraction through an iPhone accelerometer.  In: 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE; 2011. p. 5260–5263.
  * 23. Mayagoitia RE, Nene AV, Veltink PH.  Accelerometer and rate gyroscope measurement of kinematics: An inexpensive alternative to optical motion analysis systems.  Journal of Biomechanics. 2002;35(4):537–542.  {10.1016/S0021-9290(01)00231-7}.
  * 24. Simonyan K, Zisserman A.  Very deep convolutional networks for large-scale image recognition.  3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings. 2015; p. 1–14.
  * 25. Kotsiantis S, Kanellopoulos D, Pintelas P.  Handling imbalanced datasets : A review.  Science. 2006;30(1):25–36.  {10.1007/978-0-387-09823-4_45}.
  * 26. Nagi J, Ducatelle F, Di Caro GA, Cireşan D, Meier U, Giusti A, et al.  Max-pooling convolutional neural networks for vision-based hand gesture recognition.  2011 IEEE International Conference on Signal and Image Processing Applications, ICSIPA 2011. 2011; p. 342–347.  {10.1109/ICSIPA.2011.6144164}.
  * 27. Ordóñez FJ, Roggen D.  Deep convolutional and LSTM recurrent neural networks for multimodal wearable activity recognition.  Sensors (Switzerland). 2016;16(1).  {10.3390/s16010115}.
  * 28. Sak H, Senior A, Beaufays F.  Long Short-Term Memory Based Recurrent Neural Network Architectures for Large Vocabulary Speech Recognition.  INTERSPEECH-2014. 2014; p. 338–342.
  * 29. Kingma DP, Ba JL.  Adam: A method for stochastic optimization.  3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings. 2015; p. 1–15.
  * 30. Cortes C, Mohri M, Rostamizadeh A.  L2 Regularization for Learning Kernels.  Optimization. 2012;.
  * 31. Hinton GE, Srivastava N, Krizhevsky A, Sutskever I, Salakhutdinov RR.  Improving neural networks by preventing co-adaptation of feature detectors. 2012; p. 1–18.
  * 32. Gorodkin J.  Comparing two K-category assignments by a K-category correlation coefficient.  Computational Biology and Chemistry. 2004;28(5-6):367–374.  {10.1016/j.compbiolchem.2004.09.006}.
  * 33. Delgado R, Tibau XA.  Why Cohen’s Kappa should be avoided as performance measure in classification.  PLoS ONE. 2019;14(9):1–26.  {10.1371/journal.pone.0222916}.
  * 34. Simonyan K, Vedaldi A, Zisserman A.  Deep inside convolutional networks: Visualising image classification models and saliency maps.  2nd International Conference on Learning Representations, ICLR 2014 - Workshop Track Proceedings. 2014; p. 1–8.
  * 35. Arif M, Kattan A.  Physical activities monitoring using wearable acceleration sensors attached to the body.  PLoS ONE. 2015;10(7):1–16.  {10.1371/journal.pone.0130851}.