Acoustic localisation in a two-dimensional granular medium

We focus on localizing the source of acoustic emissions within a compressed two-dimensional granular ensemble of photoelastic disks, having as main information the arrival times of the acoustic signal to 6 sensors located in the boundaries of the system. By estimating, thanks to the photoelasticity of the grains, the wave speed at every point of the structure, we are able to compute the arrival times from every point of the system to the sensors. A comparison between the arrival time differences between every set of computed values to those from the actual measurements allows finding the source of the acoustic emissions.


Introduction
Besides its practical and industrial relevance, granular materials play an important role as simplified models of complex phenomena. Jamming transition [1][2][3], Selforganized criticality [4,5] and earthquake dynamics [6][7][8][9] are some significant examples. In the case of real seismicity, acoustics is the main source of information, and with three noncollinear seismometers it is -in principlepossible to localise the epicentre of an earthquake. Building earthquake models having acoustics as a main measured variable is therefore essential into establishing relevant links between the model and the real phenomenon. However, in the case of models based on granular materials, localising the source of an acoustic emission is a challenging task.
The reason for that is simple: the fronts of the sound waves are not spheres, but they get affected by the heterogeneous character of the structure. Indeed, loads are transferred through contact mechanisms between neighbouring particles creating a network of force chains [10], which is responsible for a huge degree of heterogeneity inside the granular structure [11]. Both the speed and amplitude of sound waves depend on the local value of the force network [12,13]. Consequently, waves propagate faster within the strong links of the network, and get attenuated [14,15] within the weak links.
In our quest of using a sheared granular system as an analogue to a tectonic fault [9], we focus this work on localising the source of acoustic emissions within a compressed two-dimensional granular system, having as main information the arrival time of the acoustic signal to six sensors located in the boundaries of the system. By estimating the wave speed at every point of the structure, we are able to compute the arrivals times from every point of the system to the sensors. A comparison between the arrival time differences between every set of computed val-

Experimental setup
The experimental setup consists of a one layer of roughly 2200 disks sandwiched between two vertical and parallel acrylic plates 4.4 mm apart, and confined into a rectangular cell of initial size of 335 × 295 mm 2 (W × H) (see figure 1). The disks have been 3D-printed in Durus-white material -which is photoelastic-using an Objet30 printer. They are a bi-disperse mixture of 6.4 mm and 7 mm diameter (in equal proportion), and have 4 mm thickness. The bottom and lateral walls are fixed, while the upper one is mobile allowing the compression of the granular medium Figure 1. Schematic of the experimental setup. Grains are arranged in a 2D layer held by a shape made of insulating foam. Constant pressure is applied on the top and acoustic signals are recorded by 6 peripheral sensors. Artificial waves are created by a shaker, hitting one grain with a metallic pin, driven with a square signal. with a dead weight of 30 kg. All walls are covered with insulating foam. This guarantees that the acoustics propagates exclusively within the granular structure, as the mechanical coupling between the grains and the acrylic plates is rather weak. Six acoustic sensors (VP-1.5 from CTS Valpey Corp.), in contact with the granular structure, are placed in holes in the lateral and top walls. One controlled acoustic emission is generated thanks to a thin metallic rod connected to a permanent magnet shaker (LDS V201 from Brüel Kjaer) that hit one grain at the bottom of the granular ensemble. The shaker is driven with a squared signal. The acoustic emission is captured by the six sensors and recorded at a frequency of 2 MHz using a NI-USB-6366 card.

Hyperbolic detection
Our main goal is localising an acoustic emission spontaneously generated within the medium, due to a reorganisation event. In the controlled case of injecting an acoustic signal into the structure that we are analysing in this article, we do not consider the emission time of an event, but we only measure the arrival times to each sensor. Let us suppose that we have access to the arrival times (t 1 , t 2 ) at two different locations (S 1 , S 2 ), in a continuous isotropic medium of wave velocity c. Then, the possible emission location P are located on an hyperbole, as the following relationship is followed: In our case, we have access to the sign of the time difference as we know which signal came first. This reduces the full hyperbole to only one of its branches, eliminating absolute values in the equation. Adding a third sensor S 3 , the number of different sensor pairs increases to three. Three hyperbole branches can then be computed, theoretically pinpointing the source accurately. Practically speaking, precision increases with every additional sensor due to uncertainty. Indeed, at 2 MHz and with a sound velocity approximated to 800 m/s, a 16-points duration in the acoustic signal translates to a propagation distance of ∼1 disk diameter for the wave. As our medium cannot be modeled as either homogeneous, isotropic or continuous, the locus of source point is not a true hyperbolic branch but rather a pseudo-hyperbolic branch. However, the two will be used interchangeably from now on.

Propagation model
We chose a simple ballistic propagation model with two key ingredients. The first one is the topology of the pile: sound propagation can only occur within disks. The second one is heterogeneous sound velocity, caused by the heterogeneous force chains [9]. The pile is modeled as a network with nodes being points of contacts between grains and edges being paths between the contacts of a given grain, as represented in figure 2. The choice of contacts instead of grain centers as nodes is motivated by two reasons: (i) the origin of an acoustic emission is always a node because sounds are produced by grains sliding and hitting each others, and (ii) the edges between such nodes only go through a single grain, making them easier to model. Grains and their contacts are determined through the Hough circle detection method [16]. Each edge must then be attributed its own propagation velocity. Its estimation is based on the relations described in [13] between applied force F and sound speed c: by using the same 3D-printed grains within a disordered array and low force values, c ∼ √ F. Force values are roughly estimated using the photoelasticity of the grains. Only one light fringe is developed during the deformation of the grains, which is a consequence of the high Young modulus of the material Y ∼ 100 MPa, compared to classical experiments using photoelastic disks with Y = 4 MPa [17]. As the deformations are small, and therefore only one light fringe is developed, the gradient square (G 2 ) method delivers a good approximation of the force values for each disk [17]. The distribution of G 2 is then mapped to a speed distribution by the formula: V = α G 2 + β, with α and β calibrated in the setup. Finally, by dividing an edge length by its associated velocity, we can finally obtain the wave flight time along that edge. To know the total propagation time of the wave's ballistic part from a given node to a given sensor, we compute the fastest path between the two points in this flight time network with the Dijkstra algorithm [18] and sum the flight times along all edges in the path. For a given potential source point, the difference in arrival time between two sensors follows with N and M being the number of nodes of the fastest paths for each of the two sensors. This makes ∆t comp ("computed") dependent on up to a hundred or more measurements depending on how far the considered source is, which is an important source of uncertainty.

Single branch likelihood
To evaluate whether a point A k in the network can be the source of a measured event according to a pair of sensors (i, j), we want to know how likely it is that point sits on the pseudo-hyperbolic branch stemmed from that sensor pair. We can model each ∆t value (measured or computed) as a random variable drawn from a normal distribution N( ∆t, σ 2 ), centered on the true value ∆t with variance equal to estimated uncertainty σ. Defining the difference E ≡ ∆t meas −∆t comp , it follows that E can be modelled with N( E, σ tot ) with E the true value and σ 2 tot = σ 2 comp + σ 2 meas . What we then want to know is, how likely it is that E is close to 0 -or more formally, that | E| < for a given small and positive : The above value depends on which pair of sensors i, j and which node k are being considered and can be written as P(i, j, k). We drop from the notation as it only plays a role similar to contrast ratio, but does not change ordering.
To makeP a proper probability, we can normalize by its sum over all nodes: Figure 3 represents this probability field overlaid on a photo of the experiment for sensor pair (3,5). As expected, the maximums are similar to a hyperbolic branch, with one "root" between the two sensors. The branch widens as it gets further away from the sensors, due to the distancedependent uncertainty discussed in section 3.2. Its interesting to note that due to the discrete nature of the network and the real-valued measurements, edges have slightly different and almost certainly unique flight times. With perfectly accurate values, a single pair of sensors could provide all the information needed to find the source as the probability maximum over every grain is almost certainly uniquely defined. practically speaking however, uncertainty forces us to use and combine information from all sensors.

Agglomerated likelihood
We have seen how to derive a probability for node k from a pair of sensors i, j, but 15 different pairs are available. How to efficiently exploit all this information is not obvious.
As described in [19], several techniques exist: variable-wise, based on the whole distribution, linear, logarithmic, and so on. We chose a product based pooling with a formula of the form with the pooled probability P A being a product combination of each pair probability. Weights are path-dependent w = w i, j, k . The only condition on weights is that P A (k) must be defined -in particular, they do not have to sum to 1. The weights dependency on each pair (i, j) makes sense as not all signals received will have the same signalto-noise ratio. The less obvious choice of "expert weight", depending on the point k , is motivated by the uncertainty described in 3.2. When propagating to points further away from a sensor, errors accumulate and may counterbalance the gain from an otherwise "good" pair with accurate signals.
To favor more precise information, weights have been set to w i, j (k) = 1/σ tot . Another popular choice to combine probabilities is linear combination. This was decided against because product pooling generally gives unimodal distribution, more adapted to our case as we look for a single source. Product pooling also has the desirable property of eliminating any outcome which is strongly voted against (w i, j (k) ∼ 1 and P i, j (k) 1) by a constituent probability. The result of the method using all six sensors is displayed on figure 4, localizing the source at 2 grain diameter distance from its actual position. If the same excitation is reproduced, the acoustic signal is quite robust, and therefore we obtain exactly the same result. In other configurations and using different acoustic sources, the accuracy of the method remains under a few grain diameters.

Discussion
While the presented method has good localisation capacity, several caveats must be addressed. First, the influence of the parameter has been described as minimal -this is true for a given range, typically not more than two orders of magnitude away (× 0.01 ∼ 100) from typical σ values. For very large (respectively, very small) , P A (k) converges to 1 (respectively, 0) for all k and ordering can change in these extreme ranges. A possible way to eliminate the issue can be to set = σ i, j (k), but further computational analysis has to be made to check whether this choice could introduce unwanted bias. In addition, the influence of distance to the sensor (through uncertainty) could not be tested for very long range detection due to the limited dimensions of the setup.
One main limitation of the method is the difficulties in measuring with enough accuracy the ballistic arrival time of the acoustic signal to the sensors. As the acoustic signal is generated by the collision between a metallic rod and one disk, the increment of its amplitude is relatively slow. This makes it difficult to define a proper threshold either on the amplitude or on its time derivative to define the arrival of the signal.
We are currently tackling this issue with different approaches: first by modifying the source of the acoustic emission with the addition of an active grain. The first attempt of the device is shown in figure 3. It is capable of generating a controlled acoustic emission in the central part of the cell, far from its boundaries. A circular plastic box contains a hollow metallic cylinder that will move in the presence of a local magnetic field. A rapid movement provokes that the cylinder hits the walls of the plastic box, generating a sound. However, the quality of the signal is worse than the one with the shaker and the accuracy of the localisation is about five grain-diameters. Other more controlled sources of acoustic emission are currently under inspection, based on the fracture of a pencil lead inside a 3D-printed grain and triggered by the deformation of the grain. This can increase the characteristic frequency of the signal, which will improve the accuracy of the localisation. The second approach is simulating the time-reversal of the acoustic signals [20] that arrive to the sensors. This will require a new modelling of the system, accounting for wave propagation. This approach has two advantages: there would be no need to measure any timing as all signals are retro-propagated in sync; and all the information contained in the signals would be used, instead of a single point in the case of ballistic approach.

Conclusion
In a compressed 2D granular layer, sound propagates preferentially through force chains due to heterogeneities in wave velocity. This local velocity field can be estimated thanks to photoelasticity. Using this information, we can model the pile as a network in which ballistic wave-front propagates. Six acoustic sensors, located on the boundaries of the system, provide arrival timings of this wavefront. Comparing the pattern of timing differences to what the model would predict they would be for all possible source points, we can map the acoustic source probability and find the actual source down to a 2 grain diameters.
We acknowledge financial support from the AAP-iLM-2020.