ABSTRACT

A flash bang is a non-lethal explosive device that delivers intensely loud bangs and bright lights to suppress potentially dangerous targets. It is usually used in crowd control, hostage rescue and numerous other missions. We construct a model for assessing quantitatively the risk of hearing loss injury caused by multiple flash bangs. The model provides a computational framework for incorporating the effects of the key factors defining the situation and for testing various sub-models for these factors. The proposed model includes 1) uncertainty in the burst point of flash bang mortar, 2) randomness in the dispersion of multiple submunitions after the flash bang mortar burst, 3) decay of acoustic impulse from a single submunition to an individual subject along the ground surface, 4) the effective combined sound exposure level on an individual subject caused by multiple submunitions at various distances from the subject, and 5) randomness in the spatial distribution of subjects in the crowd. With the mathematical model formulated, we seek to characterize the overall effect of flash bang mortar in the form of an effective injury area. We carry out simulations to study the effects of uncertainty and randomness on the risk of hearing loss injury of the crowd. The proposed framework serves as a starting point for a comprehensive assessment of hearing loss injury risk, taking into consideration all realistic and relevant features of flash bang mortar. It also provides a platform for testing and updating component models.

Keywords:

Risk of Significant Hearing Loss, Mathematical Framework for Assessing Injury Risk, Effective Injury Area, Decay of Acoustic Impulse along Ground Surface, Dose-Response Relation, Fluctuations in Actual Injury Numbers

1. Introduction

In recent years, armed gangs, militias, and terrorist cells have become more prevalent in modern asymmetric conflict and irregular warfare. In the 1950s noncombatants accounted for about one-half of all US military operation casualties and the rate rose to about 80% in the 1980s [1]. Military forces today must be able to execute missions across a large range of military operations. This spans from stability operations, disaster response and humanitarian assistance to full-scale armed combat. Non-lethal weapons can allow for tailored responses to targets and situations across this continuum and can provide commanders the flexibility with escalation-of-force options to minimize civilian casualties and collateral property damage [2]. Further, according to the US Department of Defense Non-Lethal Weapons program, while non-lethal weapons traditionally have supported operations such as peacekeeping and humanitarian assistance, there is a growing appreciation for these weapons, devices and munitions in irregular warfare operations such as counterinsurgency, counterterrorism, stability operations, and counter-piracy [3].

Non-lethal weapons have been successfully employed in the engagements with potential threats, including counterinsurgency operations, peacekeeping operations, humanitarian efforts, crowd and riot control, and crisis management [6] [7]. For example, in 1995 US forces in Somalia successfully utilized non-lethal weapons to preclude injury to civilians in support of humanitarian operations. In 2000, a US military police unit used non-lethal weapons to disperse a violent rocking-throwing and stick-wielding crowd and to provide protection for the peacekeeping personnels during international peacekeeping operations in Kosovo.

Non-lethal weapons are explicitly designed and primarily employed to incapacitate targeted personnel or material immediately, while minimizing fatalities, permanent injury to personnel, and undesired damage to property in the target environment. A key characteristic of non-lethal weapons is that they are intended to have reversible effects on personnel and material [4]. Generally speaking, conventional lethal weapons, such as explosive-filled warheads, damage or kill their targets through blast, penetration and fragmentation [5]. In contrast, non-lethal weapons employ means other than catastrophic physical destruction to interrupt the opponent’s normal functions. In order to use non-lethal weapons judiciously and effectively, it is important to be able to assess the risks and damages associated with applying various non-lethal weapons. From a mathematical point of view, predicting the outcome of an area non-lethal weapon used against a crowd can be challenging since in both achieving the desired effect and causing undesired injury, human effects play an important role and need to be taken into consideration [8].

One widely used type of non-lethal weapon is the flash bang munition. Flash bang devices are designed to deny access into and out of an area, move individuals through an area, or cause suppressive effects. These devices deliver a bright flash and loud bang used for crowd control and room clearing missions. The primary potential risk of injury from flash bang devices comes from impulse-noise-induced permanent hearing loss. Recently hearing loss injury associated with multiple acoustic impulses was examined in a study where an empirical logistic dose-response model was developed [9]. This empirical injury model was theoretically interpreted in our previous studies, from the point of view of immunity [10] and biovariability [11] [13] , respectively.

In this study we aim at building a mathematical framework that takes into consideration the key factors describing the situation of applying flash bang mortars on a crowd. These factors include 1) uncertainty in the burst point of flash bang mortar, 2) randomness in the dispersion of multiple submunitions after the burst of flash bang mortar, 3) decay of acoustic impulse from a single submunition to an individual subject along the ground surface, 4) the effective combined sound exposure level (SELA) on an individual subject caused by multiple submunitions at various distances from the subject, and 5) randomness in the spatial distribution of subjects in the crowd. With the framework developed, we study the risk and the fluctuations in the occurrences of significant hearing loss injury caused by multiple submunitions on a crowd.

2. Background and Formulation

A flash bang munition is a non-lethal explosive device that emits a dazzling flash of light and a thunderous noise impulse to temporarily disorient the senses of affected individual subjects. In particular, the bright flash can induce temporary flash blindness that lasts seconds whereas the loud blast is at 170 decibels or more and can cause temporary large threshold shifts in hearing [14]. The US military is developing a long range non-lethal mortar round that delivers a flash bang payload. When the mortar round gets close to the targeted area, it bursts to release multiple submunitions that are dispersed over an elliptical area. After falling to the ground, the submunitions are ignited to generate optical and acoustic impulses. In this study, we focus on the hearing loss due to the acoustic impulses from the submunitions. Below we first describe the model components for constructing a comprehensive modeling framework.

1) The ground surface and the spatial distribution of subjects in the crowd.

2) The burst of flash bang mortar and the spatial dispersion of submunitions.

3) The decay of acoustic impulse vs distance along ground surface.

4) Effective combined SELA caused by multiple submunitions.

5) Logistic dose-response relation predicting injury probability from SELA.

We then assemble these components into a computational framework for assessing the risk and fluctuations in the occurrences of hearing loss injury of a crowd, caused by multiple submunitions dispersed in the air over the crowd.

2.1. The Ground Surface and the Crowd

We establish the coordinate system as follows: the x-axis is the range direction from the mortar launch position to the target (center of the crowd), the y-axis the deflection direction, and the z-axis the vertical direction. We consider a ground surface, with possible variation in height given by function $z=H\left(x,y\right)$ . A crowd on the ground surface is distributed in the $\left(x\mathrm{,}y\right)$ dimensions according to a given probability density, such as, a normal distribution, or a uniform distribution inside a bounded region. As an example, we consider a crowd uniformly distributed in a circle. We put the origin of the coordinate system at the center of crowd. Let

1) $n_c$ : number of subjects in the crowd.

2) $d_c$ : diameter of the circle formed by the crowd.

3) $h_{\text{ear}}$ : height of ears from the ground surface, for a random subject in the crowd.

The uniform distribution of $n_c$ subjects inside the circle of diameter $d_c$ is sampled as

$\{\begin{array}{l}{x}^{\left(k\right)}=\frac{1}{2}{d}_c\sqrt{u_k^{\left(1\right)}}\cos \left(2\text{π}{u}_k^{\left(2\right)}\right),\hfill \\ y^{\left(k\right)}=\frac{1}{2}{d}_c\sqrt{u_k^{\left(1\right)}}\sin \left(2\text{π}{u}_k^{\left(2\right)}\right),\hfill \end{array}k=1,2,\cdots ,n_c$ (1)

where $\left\{u_k^{\left(1\right)}\right\}$ and $\left\{u_k^{\left(2\right)}\right\}$ are independent samples of the uniform distribution in $\left[\mathrm{0,1}\right]$ .

For each subject, the point of focus for assessing hearing loss injury is the ears. The z-coordinate of ears of each subject is specified by

$z^{\left(k\right)}=H\left(x^{\left(k\right)}\mathrm{,}{y}^{\left(k\right)}\right)+h_{\text{ear}}^{\left(k\right)}$ (2)

where $h_{\text{ear}}$ is the relative height of ears from the ground surface. In general, $h_{\text{ear}}$ is a random variable described by a probability distribution. This can accommodate the situation where some subjects in the crowd are standing while some others are sitting or squatting. The distribution of $h_{\text{ear}}$ can also be used to model the heterogeneity in height among subjects when they are all standing. For example, we can set $h_{\text{ear}}=1.7\text{m}$ to model the idealized situation where all subjects are standing and all have the standard height. On the other hand, to model a crowd in which 60% of subjects are sitting while the rest 40% are standing, we can use

$h_{\text{ear}}=\{\begin{array}{ll}1.7\text{m}\mathrm{,}\hfill & \text{with}\text{probability}0.4\hfill \\ 1.25\text{m}\mathrm{,}\hfill & \text{with}\text{probability}0.6\hfill \end{array}$

2.2. Burst of Mortar Round and Dispersion of Submunitions

Next we describe the spatial distribution of the multiple submunitions released when a mortar round bursts. After being released and dispersed, the submunitions controlled by time delay fuses, by design, will fall to the ground surface before being ignited to discharge optical and acoustic energy. In the current study, we focus on this idealized situation. In subsequent studies, we will include the possible failure of time delay fuses, which leads to uncertainty in the altitude of submunition ignition. Here, we only need to consider the $\left(x\mathrm{,}y\right)$ coordinates of a submunition at its ignition, which is the outcome of the random dispersion; its z coordinate at ignition is calculated from $\left(x\mathrm{,}y\right)$ coordinates based on the ground surface: $z=H\left(x,y\right)$ . The dispersion of submunitions in the $\left(x\mathrm{,}y\right)$ dimensions is characterized by an elliptical area aligned with the range and the deflection directions. Let

1) $X_d=\left(x_d\mathrm{,}{y}_d\mathrm{,}{z}_d\right)$ : the actual burst point of main munition releasing sub-munitions.

2) $\sigma$ : standard deviation of Gaussian aiming error in each dimension of $X_d$ .

3) $n_s$ : number of submunitions released from a mortar round.

4) $d_{\text{rng}}$ : range dimension of the elliptical dispersion area of submunitions.

5) $d_{\text{defl}}$ : deflection dimension of the elliptical dispersion area of submunitions.

The actual burst point $X_d$ is a random variable. Relative to the intended burst point ( $X_d^{\left(\text{intended}\right)}$ ), the actual burst point is modeled as a normal distribution.

$X_d=X_d^{\left(\text{intended}\right)}+N\left(\left(\mathrm{0,0,0}\right)\mathrm{,}{\sigma }^2{I}_3\right)$

The elliptical area of dispersion is specified by its major and minor axes, and is concisely represented as $d_{\text{rng}}\times d_{\text{defl}}$ . We model the dispersion of submunitions as concentrated in a fat elliptical annulus that is contained in the elliptical area centered at the actual burst point $X_d$ . The spatial dispersion relative to $X_d$ has the distribution density:

$\rho \left(x\mathrm{,}y\right)\propto exp\left(-\frac{2}{{\delta }^2}{\left[\sqrt{{\left(\frac{2x}{d_{\text{rng}}}\right)}^2+{\left(\frac{2y}{d_{\text{defl}}}\right)}^2}-\left(1-\frac{\delta }{2}\right)\right]}^2\right)$ (3)

where $\delta$ is the relative breath of the annulus, which is the ratio of breadth to the outer radius of the annulus. The relative breadth $\delta$ measures how fat the annulus is. Figure 1 illustrates the dispersion of submunition as given in (3) with parameters $d_{\text{rng}}=25\text{m}$ , $d_{\text{defl}}=15\text{m}$ and $\delta =0.8$ . The dashed white line depicts the $d_{\text{rng}}\times d_{\text{defl}}$ elliptical area that characterizes the dispersion.

In Monte Carlo simulations, we need to draw independent samples for the dispersion of submunitions. This issue is addressed in Appendix where we derive an exact and efficient algorithm for drawing independent samples from distribution (3).

Another reasonable yet simpler approach for modeling the dispersion of submunitions is to treat them as uniformly distributed in the elliptical area. That is, after dispersion, the positions of submunitions relative to $X_d$ are independent samples of the uniform distribution over the $d_{\text{rng}}\times d_{\text{defl}}$ elliptical area. In that case, samples of submunitions dispersion are computationally generated by

$\{\begin{array}{l}{x}_j=\frac{1}{2}{d}_{\text{rng}}\sqrt{u_j^{\left(1\right)}}cos\left(2\text{π}{u}_j^{\left(2\right)}\right)\mathrm{,}\hfill \\ y_j=\frac{1}{2}{d}_{\text{defl}}\sqrt{u_j^{\left(1\right)}}sin\left(2\text{π}{u}_j^{\left(2\right)}\right)\mathrm{,}\hfill \end{array}j=\mathrm{1,2,}\cdots \mathrm{,}{n}_s$ (4)

where $\left\{u_j^{\left(1\right)}\right\}$ and $\left\{u_j^{\left(2\right)}\right\}$ are independent samples of the uniform distribution in $\left[\mathrm{0,1}\right]$ .

2.3. Decay of Acoustic Impulse vs Distance along Ground Surface

The A-weighted sound exposure level (SELA) is an effective single metric for predicting the injury risk [15]. For that reason, SELA is selected as the dose in the dose-response relation for predicting injury probability. At a subject’s ears, the SELA value caused by a single submunition varies with the distance between the two. We consider two models. Model A is an empirical model from [16] , which does not explicitly count for the energy dissipation in acoustic wave propagation. We propose Model B, a revised model that includes both the power law decay of energy per area due to the expansion of the spherical wave and the exponential decay due to energy dissipation in wave propagation. We fit Model B to the experimental data from [16] to determine the model parameters.

Let r denote the distance between the sound source and the target. To distinguish models A and B, the SELA at the target is denoted, respectively, as $F_A\left(r\right)$ and $F_B\left(r\right)$ .

Model A from [16] :

$\begin{array}{l}\text{Model}\text{A}:F_A\left(r\right)=-11.59\ln \left(r\right)+150.36\\ =-26.69{\log }_{10}\left(r\right)+150.36\end{array}$ (5)

Model B proposed in this study:

From a submunition viewed as a point source, the emitted acoustic impulse propagates out as a spherical wave. At distance r from the point source, the area of spherical wave front is proportional to r². In the absence of energy dissipation, the energy per area at distance r is inversely proportional to r². The sound energy dissipation does occur during wave propagation, partly due to viscous dissipation in the fluid medium (the air). When both the sound source and the receiving target are on the ground, however, the interaction between the sound wave and the ground surface may cause a sound energy loss significantly larger than that attributed to the viscous dissipation in open air. We model the sound energy loss phenomenologically as an exponential decay with the distance traveled for all relevant frequencies in human hearing loss injury. A more detailed model for acoustic wave propagation would have a frequency dependent energy dissipation [12]. Combining the spherical wave expansion and the energy dissipation, we model the sound energy per area at distance r as

$\left(\text{Energy}\text{per}\text{area}\right)=E_0{r}^{-2}exp\left(-\kappa r\right)$

It follows that the SELA in dBA decays with the distance in the form of

$\begin{array}{c}\left(\text{SELA}\text{in}\text{dBA}\right)=10{\log }_{10}\left(\text{Energy}\text{per}\text{area}\right)+\text{const}\\ =-20{\log }_{10}\left(r\right)-c_{1}r+c_0\end{array}$ (6)

To determine the coefficients $c_0$ and $c_1$ , we fit the proposed model (6) to the data measured in [16]. The result is Model B displayed below.

$\text{Model}\text{B}\mathrm{<mo>\; :\; </mo>}{F}_B\left(r\right)=-20lo{g}_{10}\left(r\right)-0.5r+148.5$ (7)

Figure 2 plots SELA vs distance, respectively, for the experimental measurements from [16] , for Model A from [16] , and for Model B proposed in this study. The comparison in Figure 2 is for distance between 0.5m and 10m, the distance range of the experimental measurements in [16]. Over this short distance, the difference between Model A and Model B is insignificant. Both models match the data points fairly well. Over a longer distance, however, the linear decrease component of Model B in SELA (corresponding to exponential decay in energy) will set the two models apart. Figure 3 compares the two models over the distance range up to 80m. Clearly, over a longer distance, Model B predicts a SELA value significantly smaller than the one from Model A.

The two models above are both empirical models, each obtained by fitting a given function form to experimental measurements. But they are based on different model forms. Model A assumes that the acoustic energy per area decays as a power of distance without specifying the exponent of the power law. Mathematically, in Model A, the SELA in dBA decays linearly with respect to log (distance). The coefficients in the linear function are determined by fitting to short range data. The estimated parameters in Model A suggest that energy per area at distance r is phenomenologically proportional to $r^{-2.67}$ .

$\text{Model}\text{A}\mathrm{<mo>\; :\; </mo>}\left(\text{Energy}\text{per}\text{area}\right)=E_0{r}^{-2.67}$

In contrast, in Model B, the decrease caused by spherical wave expansion is modeled as being proportional to $r^{-2}$ based on the physical model of wave expansion. To count for the observed decrease faster than predicted by $r^{-2}$ , Model B contains a separate exponentially decaying factor attributed to the viscous dissipation and to the interaction with ground surface. To test models A and B in real experiments, more measurements are needed on the energy decay of an acoustic impulse propagating along the ground surface over distances beyond 10 meters.

2.4. Effective Combined SELA Caused by Multiple Submunitions

After being dispersed in the air over the crowd and falling to the ground, the $n_s$ submunitions are ignited individually to generate acoustic and optical impulses. The ignitions of $n_s$ submunitions are not completely synchronized but they occur within a short period of time (typically within 3 seconds [16] ). Let

1) $X^{\left(k\right)}$ : position of subject k

2) $Y_j$ : position of submunition j

3) $S_j^{\left(k\right)}$ : SELA on subject k caused by submunition j

4) $S_{\text{comb}}^{\left(k\right)}$ : effective combined SELA on subject k from all $n_s$ submunitions.

$S_j^{\left(k\right)}$ is determined by the distance between the submunition and the subject.

$S_j^{\left(k\right)}=F\left(\left|X^{\left(k\right)}-Y_j\right|\right)$ (8)

where function $F(\cdot )$ is either Model A or Model B described above. The combined SELA from all $n_s$ submunitions, $S_{\text{comb}}^{\left(k\right)}$ , is calculated using the dose combination rule [10].

$S_{\text{max}}^{\left(k\right)}=max\left\{S_j^{\left(k\right)}\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}{n}_s\right\}$

$S_{\text{comb}}^{\left(k\right)}=S_{\text{max}}^{\left(k\right)}+\frac{\lambda }{ln10}ln\left(\sum _{j=1}^{n_s}exp\left[\frac{ln10}{\lambda }\left(S_j^{\left(k\right)}-S_{\text{max}}^{\left(k\right)}\right)\right]\right)$ (9)

where, for combining 25 or fewer impulses, coefficient $\lambda$ takes the value $\lambda =3.44$ [9] , which describes the situation of flash bang mortars.

In Monte Carlo simulations, once the positions $\left\{Y_j,j=1,2,\cdots ,n_s\right\}$ of $n_s$ submunitions after dispersion are generated, the effective combined SELA from all submunitions on subject k located at $X^{\left(k\right)}$ is calculated using Equation (9). The combined SELA, $S_{\text{comb}}^{\left(k\right)}$ , is then used in the dose-response relation to predict the injury probability.

2.5. Logistic Dose-Response Relation

We review the injury model from [9]. Here we make it into a short sub-section for readers’ convenience and for facilitating the flow of presentation.

In [9] , an empirical dose-response relation was developed based on extensive data from hearing loss experiments using a chinchilla model [17] [18] [19]. The model on chinchillas was then calibrated with data of rifle noise on human and scaled up for describing human hearing loss. In the injury model, the dose is the effective combined SELA $S_{\text{comb}}$ , and the response is, $P\left(S_{\text{comb}}\right)$ , the injury probability of a subject exposed to dose $S_{\text{comb}}$ .

The logistic dose-response relation established in [9] has the form :

$P\left(S_{\text{comb}}\right)=\frac{1}{1+exp\left(-\alpha \left(S_{\text{comb}}-{\text{ID}}_{50}\right)\right)}$ (10)

where ID₅₀ is the median injury dose and coefficient $\alpha$ controls the steepness of function. A hearing loss injury is characterized as a permanent threshold shift (PTS) above a given cut-off, for example, PTS ³ 30 dB. Table 1 lists the median

injury doses ID₅₀ for PTS cut-offs from 10 dB to 60 dB. These values of ID₅₀ are based on fittings of the logistic model to chinchilla data presented in [9]. The shape parameter $\alpha$ remains approximately unchanged for all PTS cut-offs: $\alpha =0.1/\text{dB}$ [9].

Let $P^{\left(k\right)}$ = injury probability of subject k caused by all submunitions, for PTS of a prescribed cut-off, say PTS ³ 30 dB. $P^{\left(k\right)}$ is calculated from $S_{\text{comb}}^{\left(k\right)}$ , the effective combined SELA of all submunitions, using the logistic dose-response relation:

$P^{\left(k\right)}=\frac{1}{1+exp\left(-\alpha \left(S_{\text{comb}}^{\left(k\right)}-{\text{ID}}_{50}\right)\right)}$ (11)

2.6. A Computational Framework

We assemble the modeling components set up in previous sub-sections into a computational framework for assessing the risk of hearing loss injury of a crowd, caused by multiple submunitions dispersed in the air over the crowd.

1) Generate, $\left\{X^{\left(k\right)}\right\}$ , the positions of $n_c$ subjects in the crowd

$X^{\left(k\right)}=\left(x^{\left(k\right)}\mathrm{,}{y}^{\left(k\right)}\mathrm{,}H\left(x^{\left(k\right)}\mathrm{,}{y}^{\left(k\right)}\right)+h_{\text{ear}}^{\left(k\right)}\right)\mathrm{,}k=\mathrm{1,2,}\cdots \mathrm{,}{n}_s$ (12)

where $\left\{\left(x^{\left(k\right)}\mathrm{,}{y}^{\left(k\right)}\right)\right\}$ are independent samples from distribution (1).

2) Generate, $X_d$ , a realized burst point of the mortar round

$X_d=X_d^{\left(\text{intended}\right)}+N\left(\left(\mathrm{0,0,0}\right)\mathrm{,}{\sigma }^2{I}_3\right)$ (13)

which is the superposition of the intended burst point and the aiming error.

3) Generate, $\left\{Y_j\right\}$ , the positions of $n_s$ submunitions after being dispersed in the mortar burst and falling to the ground surface

$Y_j=\left(x_d+x_j\mathrm{,}{y}_d+y_j\mathrm{,}H\left(x_d+x_j\mathrm{,}{y}_d+y_j\right)\right)$ (14)

where $\left\{\left(x_j\mathrm{,}{y}_j\right)\right\}$ are independent samples from distribution (3) or (4). $\left\{Y_j\right\}$ contain both the uncertainty in the actual burst point and the uncertainty in the dispersion of submunitions after mortar burst.

4) Calculate, $S_j^{\left(k\right)}$ , the SELA on subject k caused by submunition j, using the formula $S_j^{\left(k\right)}=F\left(\left|X^{\left(k\right)}-Y_j\right|\right)$ with Model A or Model B.

5) Calculate, $S_{\text{comb}}^{\left(k\right)}$ , the effective combined SELA on subject k from all submunitions, by applying the dose combination rule (9) on $\left\{S_j^{\left(k\right)}\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}{n}_s\right\}$ .

$\left\{S_j^{\left(k\right)}\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}{n}_s\right\}\stackrel{\text{Dose com bination rule}}{\to }{S}_{\text{comb}}^{(k)}$

6) Calculate, $P^{\left(k\right)}$ , the injury probability of subject k caused by all sub-munitions, using the dose-response relation (10) with $S_{\text{comb}}^{\left(k\right)}$ as the dose:

$S_{\text{comb}}^{\left(k\right)}\stackrel{\text{Dose-response relation}}{\to }{P}^{(k)}$

7) Let $I^{\left(k\right)}$ be the indicator function that subject k in the crowd is injured. For each $k=1,2,\cdots ,n_c$ , generate a sample of $I^{\left(k\right)}$ with injury probability $P^{\left(k\right)}$ .

$I^{\left(k\right)}=\{\begin{array}{ll}1,\hfill & u^{\left(k\right)}\le P^{\left(k\right)}\hfill \\ 0,\hfill & u^{\left(k\right)}>P^{\left(k\right)}\hfill \end{array}k=1,2,\cdots ,n_c$ (15)

where $\left\{u^{\left(k\right)}\right\}$ are independent samples of the uniform distribution in $\left[\mathrm{0,1}\right]$ .

8) In each Monte Carlo iteration, the set $\left\{I^{\left(k\right)},k=1,2,\cdots ,n_c\right\}$ contains the individual injury outcomes for the $n_c$ subjects in the crowd. For each set of parameter values, a large number of independent sample sets are generated for studying the statistical properties that characterize the crowd injury.

2.7. Setup and Parameters

Throughout this study, we use the setup and parameters below.

1) Ground surface: $H\left(x\mathrm{,}y\right)\equiv 0$ .

2) Center of the crowd in the $\left(x\mathrm{,}y\right)$ dimensions: $\left(\mathrm{0,0}\right)$ (by the design of the coordinate system).

3) Number of subjects in the crowd: $n_c=\text{variable}$ .

4) Diameter of the circle formed by the crowd: $d_c=\text{variable}$ .

5) Height of ears relative to ground: $h_{\text{ear}}=1.7\text{m}$ .

6) Intended burst point of the mortar round: $X_d=\left(\mathrm{0,0,}\ast \right)$ (in the air over the center of crowd).

7) Standard deviation of aiming error: $\sigma =\text{variable}$ .

8) Number of submunitions in a mortar round: $n_s=20$ .

9) Elliptical dispersion area of submunitions: $d_{\text{rng}}\times d_{\text{defl}}=25\text{m}\times 15\text{m}$ .

10) Acoustic strength of each submunition, described by SELA (dBA) vs distance (m):

Model A: $F_A\left(r\right)=-26.69lo{g}_{10}\left(r\right)+150.36$ , from [16].

Model B: $F_B\left(r\right)=-20lo{g}_{10}\left(r\right)-0.5r+148.5$ , (proposed model with energy dissipation, fitted to data in [16] ).

3. Injury Area Characterizing the Flash Bang Mortar’s Potential in Causing Injury

With the framework established, we carry out simulations to study the injury causing effect of multiple submunitions dispersed from a flash bang mortar. The objective of the simulations in this section is to study the risk of significant injury (RSI) as a function of location $\left(x\mathrm{,}y\right)$ . In terms of risk function $\text{RSI}\left(x\mathrm{,}y\right)$ , we compare the two models of SELA vs distance, introduced in previous section. Based on $\text{RSI}\left(x\mathrm{,}y\right)$ , we seek to find as a single metric to quantify the flash bang mortar’s potential in causing injury. We study the effective injury area as a candidate for the metric.

3.1. Risk of Significant Injury (RSI) as a Function of Location

We consider the injury caused by a flash bang mortar with multiple submunitions. We first introduce the risk of significant injury (RSI) at $\left(x\mathrm{,}y\right)$ , which is defined as the conditional probability of injury for a subject located at $\left(x\mathrm{,}y\right)$ given a realized dispersion of submunitions, averaged over all realizations of dispersion.

$\text{RSI}\left(x\mathrm{,}y\right)=E_{\text{dispersion}}\left[P\left(x\mathrm{,}y|\text{dispersion}\right)\right]$

We study the effect of submunition dispersion with respect to the mortar burst position and the effect of SELA vs distance decay. For that purpose, we shall set the aiming error to zero ( $\sigma =0$ ). Function $\text{RSI}\left(x\mathrm{,}y\right)$ is the injury risk for a subject fixed at $\left(x\mathrm{,}y\right)$ , which is unaffected by the number of subjects in the crowd ( $n_c$ ) or the radius of crowd distribution ( $d_c$ ). As a result, parameters $n_c$ and $d_c$ are irrelevant in the calculation of $\text{RSI}\left(x\mathrm{,}y\right)$ . All other parameters are as described in the previous section.

Figure 4 plots contour lines of risk function $\text{RSI}\left(x\mathrm{,}y\right)$ calculated based on, respectively, 1) Model A of SELA vs distance (left panel), and 2) Model B of SELA vs distance (right panel). The label values shown in Figure 4 are in percentage. For example, the curve labeled “0.2” near the outskirt of the right panel, represents the contour line of $\text{RSI}\left(x\mathrm{,}y\right)=0.2\%$ . Near the burst point (the center), the RSI values for the two models are comparable. As $\left(x\mathrm{,}y\right)$ gets away from the burst point, the RSI values decrease to zero. However, the pace of decrease is very different for the two models. With the exponentially decaying factor in the acoustic energy of Model B, the RSI based on Model B decreases rapidly, much faster than that of Model A. Near the burst point of the main

munition, the contour lines of $\text{RSI}\left(x\mathrm{,}y\right)$ are ellipses with longer diameters aligned with the range direction. This is attributed to the larger dispersion of submunitions along the range direction in comparison with that along the deflection direction. Away from the burst point, the contour lines of $\text{RSI}\left(x\mathrm{,}y\right)$ become more circular.

To further compare the decay of RSI between Model A and Model B, in Figure 5, we plot RSIs of both models along the range (x-direction) over longer distances. The left panel plots RSI in linear scale showing its behavior near the burst point and its decay over intermediate distance. It is already evident that the RSI of Model B decreases significantly faster than that of Model A. The right panel plots RSI in logarithmic scale, clearly demonstrating the exponential decay of RSI for Model B. The dashed line in the right panel is a fitting of $c{x}^{-0.8686}{\text{e}}^{-0.05x}$ to the RSI values of Model B. The fitting function form is based on an asymptotic analysis in the simplified case of a single submunition. The asymptotic analysis will be given in next subsection when we discuss effective injury areas. In contrast to the exponential decay of RSI for Model B, the RSI of Model A decreases much slower. The discrepancy in RSI between the two models is more pronounced at larger distance. This discrepancy is the mechanism behind an observation in next subsection that for a large crowd distributed with a constant density over a large area, Model A yields a diverging total number of injuries while in Model B, the total number of injuries remains finite, regardless of how large the area is.

3.2. Effective Injury Area of a Flash Bang Mortar

We use an effective injury area relative to the burst point as a single metric to characterize the injury causing effect of the flash bang mortar. The injury area is defined based on the risk function $\text{RSI}\left(x\mathrm{,}y\right)$ . In the simulations for studying injury areas, we shall set the aiming error to zero ( $\sigma =0$ ) since the injury area is relative to the burst point. Also, as pointed out in the previous subsection, parameters $n_c$ (number of subjects) and $d_c$ (radius of crowd distribution) are

irrelevant in the calculation of $\text{RSI}\left(x\mathrm{,}y\right)$ . Below we explore two candidates for the effective injury area.

1) Injury Risk Contour Area

Injury risk contour area of $q\mathrm{<mi>\; \%\; </mi>}$ risk, denoted by ${\text{Reg}}_{\text{Contour}}\left(q\mathrm{<mi>\; \%\; </mi>}\right)$ , is defined as the region enclosed by the contour of $\text{RSI}\left(x\mathrm{,}y\right)=q\mathrm{<mi>\; \%\; </mi>}$

${\text{Reg}}_{\text{Contour}}\left(q\mathrm{<mi>\; \%\; </mi>}\right)\equiv \left\{\left(x\mathrm{,}y\right)\mathrm{<mo>\; |\; </mo>}\text{RSI}\left(x\mathrm{,}y\right)\ge q\mathrm{<mi>\; \%\; </mi>}\right\}$ (16)

${\text{Reg}}_{\text{Contour}}\left(q\mathrm{<mi>\; \%\; </mi>}\right)$ may be empty if the threshold (q%) is specified above the maximum value of risk function $\text{RSI}\left(x\mathrm{,}y\right)$ .

2) Cookie-Cutter Type Effective Injury Area

Cookie-cutter type effective injury area of q% risk (abbreviated to cookie-cutter injury area) is denoted by ${\text{Reg}}_{\text{Cookie-cutter}}\left(q\mathrm{<mi>\; \%\; </mi>}\right)$ and is defined based on the total number of injuries. We consider a large crowd uniformly distributed with a constant density ${\rho }_0$ over a large region well beyond the effect of the mortar round. Mathematically, we treat it as an infinite crowd uniformly distributed with density ${\rho }_0$ in the infinite two-dimensional space ${ℝ}^2$ . The average total number of subjects injured has the expression

$\text{number}\text{of}\text{injuries}={\int }_{{ℝ}^2}{\rho }_0\text{RSI}\left(x\mathrm{,}y\right)\text{d}x\text{d}y$

Let $\text{Cir}\left(0\mathrm{,}{R}_B\right)$ denote the circle centered at $0$ with radius $R_B$ . We consider the situation of a hypothetical cookie-cutter style injury risk function ${\text{RSI}}_{\text{Cookie-cutter}}\left(x\mathrm{,}y\right)$ defined over circle $\text{Cir}\left(0\mathrm{,}{R}_B\right)$ as follows:

${\text{RSI}}_{\text{Cookie-cutter}}\left(x\mathrm{,}y\right)=\{\begin{array}{ll}q\mathrm{<mi>\; \%\; </mi>,}\hfill & \left(x\mathrm{,}y\right)\in \text{Cir}\left(0\mathrm{,}{R}_B\right)\hfill \\ \mathrm{0,}\hfill & \text{otherwise}\hfill \end{array}$ (17)

We set radius $R_B$ such that the number of injury for the hypothetical injury risk ${\text{RSI}}_{\text{Cookie-cutter}}\left(x\mathrm{,}y\right)$ equals to that for the actual injury risk $\text{RSI}\left(x\mathrm{,}y\right)$ .

${\int }_{{ℝ}^2}{\rho }_0\text{RSI}\left(x\mathrm{,}y\right)\text{d}x\text{d}y={\int }_{{ℝ}^2}{\rho }_0{\text{RSI}}_{\text{Cookie-cutter}}\left(x\mathrm{,}y\right)\text{d}x\text{d}y={\rho }_0\text{π}{R}_B^2\times q\mathrm{<mi>\; \%\; </mi>}$

Circle $\text{Cir}\left(0\mathrm{,}{R}_B\right)$ is named the cookie-cutter injury area of $q\mathrm{<mi>\; \%\; </mi>}$ risk:

$\begin{array}{l}{\text{Reg}}_{\text{Cookie-cutter}}\left(q\mathrm{<mi>\; \%\; </mi>}\right)\equiv \text{Cir}\left(0\mathrm{,}{R}_B\right)\mathrm{,}\\ \text{where}{R}_B=\sqrt{\frac{1}{\text{π}\times q\mathrm{<mi>\; \%\; </mi>}}{\int }_{{ℝ}^2}\text{RSI}\left(x\mathrm{,}y\right)\text{d}x\text{d}y}\end{array}$ (18)

In summary, cookie-cutter injury area of q% risk is the circle that would hypothetically yield the same total number of injuries if the injury risk were artificially set to q% inside the circle and zero outside.

While injury risk contour area may be empty if the injury risk threshold (q%) is set too high, cookie-cutter injury area is always non-empty. However, the integral in the definition of cookie-cutter injury area (18) diverges when SELA vs distance is governed by Model A. This can be seen by studying the simple case of a single submunition. Let $P\left(x\mathrm{,}y\right)$ be the injury probability for a subject at $\left(x\mathrm{,}y\right)$ caused by a single submunition at $\left(\mathrm{0,0}\right)$ . The divergence of integral is determined by the behavior of $P\left(x\mathrm{,0}\right)$ for large x. The distance from the submunition to the subject’s ears is

$r={\left(x^2+h_{\text{ear}}^2\right)}^{1/2}~x$

With Model A of SELA vs distance, the SELA at $\left(x\mathrm{,0}\right)$ is $\left(-11.59ln\left(r\right)+150.36\right)$ , which along with dose-response relation (10), yields the injury probability at $\left(x\mathrm{,0}\right)$ for Model A:

$\begin{array}{c}{P}_A\left(x\mathrm{,0}\right)=\frac{1}{1+exp\left(-\alpha \left(-11.59ln\left(r\right)+150.36-{\text{ID}}_{50}\right)\right)}\\ ~\frac{1}{1+c{x}^{11.59\alpha }}~x^{-11.59\alpha }\end{array}$ (19)

where the subscript in $P_A\left(x\mathrm{,}y\right)$ refers to Model A. For a single submunition at $\left(\mathrm{0,0}\right)$ , injury probability is axisymmetric. Using the symmetry, we write the 2-D integral in (18) over ${ℝ}^2$ as a 1-D integral

${\int }_{{ℝ}^2}{P}_A\left(x\mathrm{,}y\right)\text{d}x\text{d}y>{\int }_M^{+\infty }2\text{π}x{P}_A\left(x\mathrm{,0}\right)\text{d}x~{\int }_M^{+\infty }{x}^{\left(1-11.59\alpha \right)}\text{d}x$

The dose-response relation developed in [9] has parameter value $\alpha =0.1$ , which implies

${\int }_M^{+\infty }{x}^{\left(1-11.59\alpha \right)}\text{d}x={\int }_M^{+\infty }{x}^{-0.159}\text{d}x=+\infty$

Therefore, for a large crowd distributed over a large area with a fixed constant density, when SELA vs distance is described by Model A, the total number of injuries diverges to infinity. In contrast, Model B of SELA vs distance yields a well defined total number of injuries, depending only on the crowd distribution density, regardless of how large the crowd size is. This nice mathematical property of Model B is again demonstrated with an asymptotic analysis on the injury probability at large distance from the submunition.

With Model B of SELA vs distance, the SELA at $\left(x\mathrm{,0}\right)$ caused by a single submunition at $\left(\mathrm{0,0}\right)$ has the expression $\left(-20lo{g}_{10}\left(r\right)-0.5r+148.5\right)$ . Based on the SELA, the dose-response relation (10) yields the injury probability at $\left(x\mathrm{,0}\right)$ for Model B:

$\begin{array}{c}{P}_B\left(x\mathrm{,0}\right)=\frac{1}{1+exp\left(-\alpha \left(-20lo{g}_{10}\left(r\right)-0.5r+148.5-{\text{ID}}_{50}\right)\right)}\\ ~\frac{1}{1+c\cdot x^{8.686\alpha }\cdot {\text{e}}^{0.5\alpha x}}~x^{-8.686\alpha }\cdot {\text{e}}^{-0.5\alpha x}\end{array}$ (20)

For Model B, the 2-D integral in (18) over ${ℝ}^2$ is always convergent due to the presence of exponentially decaying factor ${\text{e}}^{-0.5\alpha x}$ ,

$\begin{array}{c}{\int }_{{ℝ}^2}{P}_B\left(x\mathrm{,}y\right)\text{d}x\text{d}y={\int }_{\text{Cir}\left(\mathrm{0,}M\right)}{P}_B\left(x\mathrm{,}y\right)\text{d}x\text{d}y+{\int }_M^{+\infty }2\text{π}x{P}_B\left(x\mathrm{,0}\right)\text{d}x\\ ~C+{\int }_M^{+\infty }{x}^{\left(1-8.686\alpha \right)}\cdot {\text{e}}^{-0.5\alpha x}\text{d}x<+\infty \end{array}$

In conclusion, cookie-cutter injury area is finite only for Model B. On the other hand, injury risk contour area is always defined and finite for both models but it may be empty if the injury risk threshold is set too high. Figure 6 shows injury risk contour area and cookie-cutter injury area, respectively, of 2%, 5% and 10% risk. All five injury areas displayed in Figure 6 are based on Model B

for SELA vs distance. Injury risk contour area of 10% risk is empty since the maximum of risk function $\text{RSI}\left(x\mathrm{,}y\right)$ is only 8.33%, below the 10% threshold specified.

4. Results and Discussion

We study the hearing loss injury risk of a crowd caused by multiple submunitions dispersed from a flash bang mortar round. We carry out simulations to calculate the fraction of injured caused by a flash-bang mortar of $n_s=20$ submunitions on a crowd of $n_c$ subjects uniformly distributed in a circle of diameter $d_c$ . We will examine both the average fraction of injured (RSI) and Monte Carlo samples of the actual injury fraction based on individual realizations of flash bang mortar burst location, submunitions dispersion and crowd subjects distribution. The problem setup and the parameters used in simulations are described in Section 2.

4.1. Average Injury Fraction as a Function of Crowd Diameter, PTS Cut-Off, and Magnitude of Aiming Error

Table 2 displays the average injury fraction as a function of two variables: 1) PTS cut-off, and 2) diameter of the crowd distribution (d_c). Each value of average injury fraction is calculated based on 100000 Monte Carlo iterations. Table 2 is for the case of no aiming error ( $\sigma =0$ ). That is, the burst of the flash bang mortar is alway at the crowd center $\left(\mathrm{0,0}\right)$ . The dispersion of submunitions after the mortar burst and the distribution of crowd subjects, however, are still random in Monte Carlo simulations.

Table 2 shows that when the crowd diameter is in the range of $d_c\le 25\text{m}$ , an increase in the crowd distribution area does not yield a significant decrease in the injury fraction. The dispersion of submunitions covers an elliptical area of $15\text{m}\times 25\text{m}$ . As long as the crowd is within the acoustic range of submunitions from this elliptical area, the average injury fraction is approximately a constant, independent of how the crowd is distributed: the average injury fraction is 7% - 8% for $\text{PTS}\ge 30\text{dB}$ while it is 4% - 4.5% for $\text{PTS}\ge 40\text{dB}$ . As the crowd diameter increases, the fixed number of subjects in the crowd are spread out over a larger circle. When the crowd circle is much larger than the elliptical dispersion area of submunitions, a significant portion of crowd is outside the most effective reach of submunitions, and as a result, the overall injury fraction decreases significantly. For a crowd uniformly distributed in a circle of diameter 60 m, the average injury fraction drops to 2.9% and 1.6%, respectively, for $\text{PTS}\ge 30\text{dB}$ and $\text{PTS}\ge 40\text{dB}$ .

After examining injury fractions for various PTS cut-offs in the situation of zero aiming error, we focus on the injury of $\text{PTS}\ge 30\text{dB}$ and instead study the effects of aiming error and other variables. Table 3 lists the average injury fraction as a function of 1) standard deviation of aiming error (σ), and 2) diameter of the crowd distribution (d_c) for hearing loss injury of $\text{PTS}\ge 30\text{dB}$ . As in the case of Table 2, each value in Table 3 is based on 100,000 Monte Carlo iterations. From the results shown in Table 3, we see that an aiming error of standard deviation 5m does not seem to change the average injury fraction appreciably. An aiming error of larger magnitude reduces the average injury fraction monotonically. At a fixed crowd distribution diameter, the larger the

magnitude of aiming error, the larger the decrease in average injury fraction. On the other hand, the effect of an aiming error of fixed magnitude, measured as the relative reduction in the average injury fraction, decreases as the crowd distribution diameter increases. For example, for a crowd uniformly distributed in a circle of diameter 10m, an aiming error of standard deviation 20 m, reduces the average injury fraction from 7.98% to 2.54%, a 68% reduction. For a crowd of diameter 100m, however, the same aiming error reduces the average injury fraction from 1.26% to 1.10%, a mere 13% reduction. This behavior is reasonable and is consistent with intuition. When the crowd distribution area is large, the aiming error simply moves the elliptical area of submunitions dispersion from the center of crowd to another part still within the crowd, resulting in only small change in the average injury fraction. In contrast, when the crowd distribution area is small, the aiming error moves the elliptical area of submunitions dispersion completely off the crowd, leading to substantial reduction in the average injury fraction.

4.2. Fluctuations in the Actual Number of Injured

We study the actual numbers of injured among individual Monte Carlo realizations, for a crowd of $n_c$ subjects uniformly distributed in a circle of diameter $d_c=20\text{m}$ . Again, we focus on the hearing loss injury of $\text{PTS}\ge 30\text{dB}$ .

Figure 7 shows the actual numbers of injured out of the $n_c$ subjects,

respectively, for a crowd of $n_c=100$ subjects (left panel) and for a crowd of $n_c=10$ subjects (right panel). For each crowd, histograms are compared among 3 aiming errors for the mortar burst position: 1) aiming error of standard deviation $\sigma =0$ (no aiming error), 2) $\sigma =5\text{m}$ , and 3) $\sigma =10\text{m}$ . In Monte Carlo simulations, each histogram is computed based on 500,000 iterations. As the number of subjects in the crowd ( $n_c$ ) increases, the histogram converges toward a normal distribution. At $n_c=10$ (right panel), the histogram is still very far from a normal distribution while for $n_c=100$ (left panel) it is getting close to a normal distribution.

The main effect of aiming error is to shift the distribution of actual injury fraction toward the lower end, and to reduce the overall average injury fraction. Table 4 compares the average injury fractions (RSI) of the two crowds ( $n_c=10$ , $n_c=100$ ) for aiming errors of various magnitudes. The average injury fractions calculated in Monte Carlo simulations are virtually the same for the two crowds, confirming the theoretical prediction that the average injury fraction is independent of the number of subjects. In contrast, the shape of the distribution of the actual injury number varies significantly from $n_c=10$ (right panel) to $n_c=100$ (left panel) as shown in Figure 7. Next we explore the possibility of predicting the behavior of actual number of injured based on the average injury fraction (RSI) and the number of subjects ( $n_c$ ).

4.3. Predictions Based on the Average Injury Fraction and the Binomial Distribution Model

Given a particular realized distribution of submunitions $\left\{Y_j\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}{n}_s\right\}$ , the injury occurrence for subject k is independent of the injury occurrence for other subjects and independent of tag k since all subjects are independently and identically distributed. Given $\left\{Y_j\right\}$ , the injury occurrence for a subject is a Bernoulli distribution with the conditional average injury probability $p\mathrm{<mo>\; |\; </mo>}\left\{Y_j\right\}$ ; the actual number of injured out of $n_c$ subjects is a random variable of binomial distribution with parameters $n_c$ and $p\mathrm{<mo>\; |\; </mo>}\left\{Y_j\right\}$ . If the randomness in the realized submunition distribution $\left\{Y_j\right\}$ does not significantly change the conditional injury probability $p\mathrm{<mo>\; |\; </mo>}\left\{Y_j\right\}$ among individual realizations, then it is reasonable to expect that the actual number of injured out of $n_c$ subjects for a random realization of $\left\{Y_j\right\}$ also has a binomial distribution with parameters $n_c$ and RSI where RSI is the overall average injury fraction

$\text{RSI}=E_{\left\{Y_j\right\}}\left(p\mathrm{<mo>\; |\; </mo>}\left\{Y_j\right\}\right)$

Let $N_{inj}$ be the number of injured out of $n_c$ subjects. The binomial distribution predicts the probability distribution of $N_{inj}$ as

$P\left[N_{inj}=k\right]=C\left(n_c\mathrm{,}k\right){\left(\text{RSI}\right)}^k{\left(1-\text{RSI}\right)}^{n_c-k}$

For a crowd of $n_c=100$ subjects, Figure 8 compares the Monte Carlo simulation results and theoretical predictions using binomial distribution in four cases of different aiming errors. Based on the results shown in Figure 8, the binomial distribution approximation is accurate when the aiming error is small. As the aiming error increases, the binomial distribution approximation deviates from the true distribution of actual injury number.

Figure 9 does the same comparison for a crowd of only $n_c=10$ subjects, instead of 100 subjects. The results of Figure 9 suggest that for smaller number of subjects, the binomial distribution approximation is relatively more accurate. The mathematical theory behind this seemingly unintuitive valid approximation for small number of subjects will be investigated in a subsequent study, in which we will also explore alternative approximations when the binomial distribution approximation breaks down.

When the binomial distribution approximation is valid, we can use it to predict probabilities of events characterizing the actual number of injured. For example, we consider the case of no aiming error ( $\sigma =0$ ) for the mortar burst position, and calculate the risk of having an actual injury number above 12 out of $n_c=100$ subjects uniformly distributed in a circle of diameter 20 m. In this case, the injury fraction averaged over all realizations is $\text{RSI}=\mathrm{7.74\; <mi>\; \%\; </mi>}$ (Table 4). We study the probability of the actual injury number exceeding a prescribed threshold $k_{\text{TLV}}$ among individual realizations.

Table 5 compares the values of this probability calculated in Monte Carlo simulations and those predicted based on binomial distribution with parameters $n_c$ and RSI. The predicted probability has the expression

$P^{\left(\text{pred}\right)}\left[N_{inj}>k_{\text{TLV}}\right]=\sum _{k=k_{\text{TLV}}+1}^{n_c}C\left(n_c\mathrm{,}k\right){\left(\text{RSI}\right)}^k{\left(1-\text{RSI}\right)}^{n_c-k}$ (21)

The results displayed in Table 5 indicate that for small or no aiming error, the

binomial distribution approximation is valid and yields reasonably accurate predictions for the true probability: $P^{\left(\text{pred}\right)}\left[N_{inj}>k_{\text{TLV}}\right]\approx P^{\left(\text{true}\right)}\left[N_{inj}>k_{\text{TLV}}\right]$ .

5. Concluding Remarks

Flash bangs are one of the commonly used anti-personnel non-lethal weapons with dual civil-military applications. In this paper, we have developed a mathematical model for computing the risk of hearing loss injury caused by multiple flash bang submunitions on a crowd. Our model includes the effects of 1) aiming error in the burst point of flash bang mortar, 2) uncertainty in the dispersion of multiple submunitions after the burst of flash bang mortar, 3) propagation of acoustic impulse from a single submunition to an individual subject along the ground surface, 4) effective combined sound exposure level on an individual subject caused by multiple submunitions ignited at various distances from the subject, and 5) randomness in the spatial distribution of subjects of the crowd. Based on the mathematical model, we explored two effective injury areas as two candidates for a single metric characterizing the overall injury causing potential of flash bang mortar. We conducted numerical simulations to study the dependence of the average injury fraction on i) magnitude of aiming error, ii) diameter of crowd distribution, and iii) PTS cut-off in defining injury. We examined the random actual injury fraction among individual realizations of mortar burst position, submunitions dispersion, and crowd subjects distribution. In the case of small or no aiming error, we found that the behavior of actual injury fraction is well predicted using a binomial distribution. This observation gives us a simple way of characterizing the actual injury fraction among individual realizations: when the binomial distribution approximation is valid, the behavior of random actual injury fraction is completely described by the number of subjects in the crowd and the overall average injury fraction. The proposed mathematical framework serves as a starting point for a comprehensive assessment of hearing loss injury risk, taking into consideration all realistic and relevant features of flash bang mortar. It provides a platform for testing and updating component models for various aspects in flash bang’s injury causing process.

Disclaimer

The authors thank the Joint Non-Lethal Weapons Directorate of US Department of Defense for supporting this work. The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the US Government.

Cite this paper

Wang, H., Burgei, W.A. and Zhou, H. (2018) Risk of Hearing Loss Injury Caused by Multiple Flash Bangs on a Crowd. American Journal of Operations Research, 8, 239-265. https://doi.org/10.4236/ajor.2018.84014

References

Appendix: Drawing Independent Samples from Distribution (3)

Our goal is to design an algorithm for drawing independent samples from distribution $\rho \left(x\mathrm{,}y\right)$ given in (3). We construct the algorithm, step by step, starting with two simple distributions. [Step 1:]

Step 1: The standard normal distribution:

${\rho }_1\left(s\right)=\frac{1}{\sqrt{\text{2π}}}exp\left(\frac{-s^2}{2}\right)$ (22)

We draw samples of ${\rho }_1\left(s\right)$ directly using the built-in function “randn” in Matlab.

Step 2: A distribution of $s>0$ :

${\rho }_2\left(s\right)=\{\begin{array}{ll}s\cdot \exp \left(\frac{-s^2}{2}\right),\hfill & s>0\hfill \\ 0,\hfill & s<0\hfill \end{array}$ (23)

We notice that random variable $\frac{1}{2}{s}^2$ has the exponential distribution

$\text{Prob}\left(\frac{1}{2}{s}^2>\eta \right)=\exp \left(-\eta \right),\text{for}\eta >0$

It follows that random variable $exp\left(-\frac{1}{2}{s}^2\right)$ is uniformly distributed in $\left[\mathrm{0,1}\right]$ .

$\text{Prob}\left(\exp \left(-\frac{1}{2}{s}^2\right)<q\right)=q,\text{for}q\in \left[0,1\right]$

Thus, independent samples of distribution ${\rho }_2\left(s\right)$ in (23) are generated by

$s_j=\sqrt{-2\log (uj)}$

where $\left\{u_j\right\}$ are independent samples of the uniform distribution in $\left[\mathrm{0,1}\right]$ .

Step 3: A mixture of distributions ${\rho }_1\left(s\right)$ and ${\rho }_2\left(s\right)$ :

${\rho }_3\left(s\right)=\frac{1}{1+c_0}\left(c_0{\rho }_1\left(s\right)+{\rho }_2\left(s\right)\right)$ (24)

Independent samples of distribution ${\rho }_3\left(s\right)$ are generated by mixing samples from ${\rho }_1\left(s\right)$ and ${\rho }_2\left(s\right)$ , respectively with probabilities $\frac{c_0}{1+c_0}$ and $\frac{1}{1+c_0}$

$s_j=\{\begin{array}{ll}{s}_j^{\left(1\right)}\text{from}{\rho }_1\left(s\right),\hfill & \text{if}{u}_j<\frac{c_0}{1+c_0}\hfill \\ s_j^{\left(2\right)}\text{from}{\rho }_2\left(s\right),\hfill & \text{otherwise}\hfill \end{array}$ (25)

where $\left\{s_j^{\left(1\right)}\right\}$ are independent samples from density ${\rho }_1\left(s\right)$ , $\left\{s_j^{\left(2\right)}\right\}$ samples from ${\rho }_2\left(s\right)$ , and $\left\{u_j\right\}$ samples from the uniform distribution in $\left[\mathrm{0,1}\right]$ .

Step 4: Scaling and shifting of ${\rho }_3\left(s\right)$ :

We scale s from distribution ${\rho }_3\left(s\right)$ by $a>0$ snd shift the result by $b>0$ to form a new variable $\xi =a\cdot s+b$ . Random variable $\xi$ has the distribution

$\begin{array}{c}{\rho }_4\left(\xi \right)=\frac{1}{a}{\rho }_3\left(\frac{\xi -b}{a}\right)\propto \left(c_0{\rho }_1\left(\frac{\xi -b}{a}\right)+{\rho }_2\left(\frac{\xi -b}{a}\right)\right)\\ \propto \left(\frac{c_0}{\sqrt{\text{2π}}}+\max \left\{0,\frac{\xi -b}{a}\right\}\right)\cdot \exp \left(\frac{-{\left(\xi -b\right)}^2}{2a^2}\right)\\ \propto \left(\frac{c_{0}a}{\sqrt{2\text{π}}}-b+\max \left\{b,\xi \right\}\right)\cdot \exp \left(\frac{-{\left(\xi -b\right)}^2}{2a^2}\right)\end{array}$ (26)

We select $c_0$ such that $\frac{c_{0}a}{\sqrt{2\text{π}}}-b=0$ . The distribution of random variable $\xi$ becomes

${\rho }_4\left(\xi \right)\propto max\left\{b\mathrm{,}\xi \right\}\cdot exp\left(\frac{-{\left(\xi -b\right)}^2}{2a^2}\right)\equiv f_4\left(\xi \right)$ (27)

Independent samples of distribution ${\rho }_4\left(\xi \right)\propto f_4\left(\xi \right)$ are generated by scaling and shifting samples from distribution ${\rho }_3\left(r\right)$ .

Step 5: Parametric form of distribution (3):

Let $\left(x\mathrm{,}y\right)$ be the vector random variable with density (3). We write $\left(x\mathrm{,}y\right)$ as

$\begin{array}{l}x=\frac{d_{\text{rng}}}{2}r\cos \left(\theta \right)\\ y=\frac{d_{\text{defl}}}{2}r\sin \left(\theta \right)\end{array}$ (28)

Scalar random variable $\theta$ is uniformly distributed in $\left[\mathrm{0,2}\text{π}\right]$ since density $\rho \left(x\mathrm{,}y\right)$ in (3) does not vary with $\theta$ . The distribution of scalar random variable r has the form

${\rho }_5\left(r\right)\propto r\cdot \exp \left(-\frac{{\left(r-b\right)}^2}{2a^2}\right)\equiv f_5\left(r\right),\text{for}r>0$ (29)

where constants a and b are related to parameter $\delta$ in (3) as

$a\equiv \frac{\delta }{2}\mathrm{,}b\equiv \left(1-\frac{\delta }{2}\right)$ (30)

Note that function $f_5(\cdot )$ in density (29) is dominated by $f_4(\cdot )$ in density (27).

$f_5\left(\xi \right)\le f_4(\xi )$

This observation enables us to draw from distribution (29) using the rejection sampling. Let $\left\{{\xi }_j\right\}$ be independent samples from distribution (27).

${\xi }_j\text{is}\{\begin{array}{ll}\text{accepted}\text{as}\text{a}\text{sample}\text{of}r,\hfill & \text{if}{u}_j\le \frac{f_5\left({\xi }_j\right)}{f_4\left({\xi }_j\right)}\hfill \\ \text{rejected},\hfill & \text{otherwise}\hfill \end{array}$ (31)

where $\left\{u_j\right\}$ are independent samples of the uniform distribution in $\left[\mathrm{0,1}\right]$ . The remaining subset of $\left\{{\xi }_j\right\}$ after the rejection selection (31) gives us a set of independent samples for distribution (29).

Step 6: Drawing from distribution (3):

We use the parametric representation of $\left(x\mathrm{,}y\right)$ given in (28). Independent samples of distribution $\left(x\mathrm{,}y\right)$ are generated by drawing samples of r from distribution (29) using the rejection sampling algorithm described above, and drawing samples of $\theta$ from the uniform distribution in $\left[\mathrm{0,2}\text{π}\right]$ .

$\begin{array}{l}{x}_j=\frac{d_{\text{rng}}}{2}{r}_{j}cos\left({\theta }_j\right)\\ y_j=\frac{d_{\text{defl}}}{2}{r}_{j}sin\left({\theta }_j\right)\end{array}$ (32)