Identifying multiply scattered wavepaths in strongly scattering and dispersive media

: The ability to extract information from scattered waves is usually limited to singly scattered energy even if multiple scattering might occur in the medium. As a result, the information in arrival times of higher-order scattered events is underexplored. This information is extracted using ﬁngerprinting theory. This theory has never previously been applied successfully to real measurements, particularly when the medium is dispersive. The theory is used to estimate the arrival times and scattering paths of multiply scattered waves in a thin sheet using an automated scheme in a dispersive medium by applying an additional dispersion compensation method. Estimated times and paths are compared with predictions based on a sequence of straight ray paths for each scattering event given the known scatterer locations. Additionally, numerical modelling is performed to verify the interpretations of the compensated data. Since the source also acts as a scatterer in these experiments, initially, the predictions and the numerical results did not conform to the experimental observations. By reformulating the theory and the processing scheme and adding a source scatterer in the modelling, it is shown that predictions of all observed scattering events are possible with both prediction methods, verifying that the methods are both effective and practically achievable.


I. INTRODUCTION
Scattering is one of the universal physical phenomena that can occur during wave propagation. It is caused by the presence of heterogeneities referred to as scatterers, which in an elastic medium might be due to density and/or velocity contrasts. Scattering plays an important role in fields ranging from the study of atoms (Gilhaus et al., 1988), nondestructive testing applications (Darmon et al., 2009), biomedical imaging (Nguyen et al., 2011), geophysical exploration (Gibson and Levander, 1988), and even space exploration (Gordon, 1958), and many of the discoveries in physics involve scattering experiments (Godbole, 2011).
In geophysics, for example, ground-penetrating radar (GPR) uses diffracted scattered waves to detect underground pipes or other structures (Gibson and Levander, 1988). In seismic methods, multiply scattered energy is usually considered undesirable since it might mask reflections from a subsurface target of interest (Pasasa et al., 1998). In seismology, scattered energy is used to study inhomogeneities in the lithosphere (Wu and Aki, 1985). Thus, the ability to extract information from scattered energy is beneficial. Most applications of scattering theory assume only single-scattering interactions to simplify the applicable theory (e.g., the Born or Rytov approximation), whereas multiple scattering always occurs in real media (e.g., Foldy, 1945). Several works have shown that neglecting multiple scattering in recorded data may lead to errors in interpretation, whereas adding information about multiple scattering can lead to improvements in the results (Gao et al., 1983;Bordier et al., 1991).
A method for predicting the arrival times and scattering paths of higher-order multiple diffractions has recently been proposed and was demonstrated on synthetic datasets (Meles and Curtis, 2014). The method takes advantage of the fact that waves from each diffracting scatterer in a medium will have a unique set of kinematics or moveout, which is referred to as its "fingerprint." The theory to find these and use them to identify or predict travel times and paths of multiply scattered events is what we refer to here as "fingerprinting theory" after Meles and Curtis (2014). To facilitate the application of the theory to real data, L€ oer et al. (2015) proposed an automated scheme that relies on several typical seismic data processing techniques such as semblance analysis and stacking. They also performed laboratory measurements to test the scheme using steel rods as scatterers embedded in a nondispersive homogeneous medium of polyvinyl alcohol gel. Due to wave attenuation and related dispersion observed in the experimental data, the scheme only detected the first-order and a part of the second-order scattering correctly, and both are essential for estimating the higher-order scattering.
In this paper, we show that the method can be used to predict or infer real scattering paths and arrival times of multiply scattered waves in a laboratory. By augmenting the method, we achieve this in dispersive media. In what follows, we first introduce the experimental setup and then the wave types and methods used to address their dispersion. We next explain the fingerprinting theory and show the results of applying it to real data. Finally we discuss the implications of this work before concluding.

II. EXPERIMENTAL SETUP
In order to apply fingerprinting theory, a laboratory study was performed at Eidgen€ ossische Technische Hochschule (ETH) Z€ urich to examine how experimental design parameters (e.g., diffractor sizes, transducer types, and wavelet types) affect the results and, hence, how to choose the right settings for the main experiment. Several key instruments are used, such as a robotic arm equipped with three laser Doppler vibrometers (LDVs) and a signal generator, a signal amplifier, and an aluminum plate with a thickness of 1.45 mm and size of 2 Â 2 m in width and length acts as a host medium for the propagating waves.
Particle velocities (displacements) on the medium surface are measured by the LDV heads from Polytec 1 (model PSV-500-3D, Germany), which is attached to a robotic arm manufactured by KUKA 2 (Switzerland), as shown in Fig. 1. The robot can move along multiple axes and so can point the scan heads toward any desired scan area. In the following, we use the term "receiver" to refer to the location at which the LDV pointed and measured the displacements. The LDV heads produce three laser beams that point to a single location on the medium's surface and by using Doppler's principle can determine the velocity of the points in that medium in the three directions of the laser beams. Particle displacements are obtained by integrating the measured velocities in time. A key advantage of using three lasers is that from the velocity measured along the three beam directions, the three components of the displacement (x, y, and z) can be determined. Therefore, it is comparable with having a three-component geophone in the typical seismic refraction or reflection measurements but without any of the usual coupling issues because the LDV provides a noncontact measurement. The LDV used in the laboratory has a sensitivity of the velocity measurements up to 0.1 lm/ s with a sampling frequency of up to 25 MHz. Because of its accuracy and ability to measure without direct contact, LDVs are used extensively in many fields, ranging from structural health monitoring (Staszewski et al., 2012) and dynamic testing of microstructures (Ngoi et al., 2000) to biomedical applications that involve testing on the human body (Casaccia et al., 2018;Chan et al., 2013;Tabatabai et al., 2013). To introduce waves in the medium, we use a transducer that only excites vibrations in the z-component (the out of plane direction) and has a surface size of 5 Â 5 mm and length of 3.5 mm with a travel range 1 mm (630%) from its central point and resonance frequency at around 2 kHz (PhysikInstrumente, 2018). The transducer converts voltages to displacements, which are amplified by an amplifier. After testing several combinations of voltages, wavelet types, and central frequencies of the wavelets, and assessing how these affect the wavelet of the propagating waves, we settled on a combination of 100 mV for the voltage and a source signal consisting of one period of a sine wave as the wavelet with a central frequency of 1.5 kHz for all experiments, including for the final measurements when applying fingerprinting theory. To create scattering, we use magnetic cubes of size 0.5 cm 3 . To measure the source signature, we put the source transducer on the back side of the plate, whereas for the final experiment, we move the transducer to the front of the plate to make it easier to move the source. The schematic of the laboratory setup used to run these experiments is shown in Fig. 1 (left).

III. LAMB WAVES AND DISPERSION COMPENSATION
The experiments use so-called Lamb waves, which are the typical guided waves used to investigate thin structures. However, Lamb waves are dispersive, meaning that different frequencies propagate with different velocities. Fundamentally, there are two modes, which are called the S0 (symmetric) mode and A0 (anti-symmetric) mode, and the dispersion of both modes can be found analytically for known medium properties. Using the open-source Waveform Revealer software 3 (Shen and Giurgiutiu, 2014) with input 1.45 mm for the thickness of the plate and density 2700 kg/m 3 , Poisson ratio 0.33, and Young's modulus 69 GPa for a typical aluminum plate, we obtain the analytical dispersion curves for both modes shown in Fig. 2 (black and red curve). We also performed a dispersion experiment to obtain the experimental dispersion curves using the multianalysis of surface waves (MASW) method, which allows us to pick the frequency-velocity relation for the dispersive waves (Xia et al., 1999). More specifically, 15 receivers were arranged in a linear fashion with a spacing of 1 cm and the first receiver coinciding with the source position (which was located on the back side of the plate). The comparison of both the analytical and experimentally measured dispersion curves using the MASW analysis are shown in Fig. 2, superimposed with the data from the dispersion experiment transformed to the frequency-wavenumber (x À k) domain. The MASW analysis reveals that only the A0 mode is excited and propagates in the plate as the experimental curve extracted from the transformed data matches the A0 analytical curve. The spectrum of the transformed data confirms this further as it overlays both of the A0 mode curves.
Fingerprinting theory and the automated scattered wavepath analysis scheme used in this study were previously tested and applied to nondispersive datasets where scattering events had relatively similar waveforms (amplitude could vary due to energy loss and geometrical spreading), whereas our experimental data are highly dispersive, which causes difficulties for the identification of arriving scattered waves (often referred to as events) due to highly deformed waveforms. We therefore perform dispersion compensation prior to identifying multiply scattered events.
To start, we acquire a data set involving a single source, receiver, and scatterer. The receiver is placed in the middle between the source and the scatterer, separated equally by 7.5 cm to ensure a (temporal) separation in the data between the direct wave and the reflection event from the nearest boundary of the plate, which was located approximately 45 cm away from the scatterer. Since the relation between the frequency and wavenumber kðxÞ is known from the dispersion curves (analytical or experimental), we can convert our dispersive data in the frequency domain G d ðxÞ to the wavenumber domain G d ðkðxÞÞ with the relation (1) The dispersion curve is now represented by the nonlinear relation of the wavenumber with respect to frequency. If it was linear, the dispersion would be removed, therefore, in this method, the nonlinearity is approximated by a linear relation and used to resample the data as if the relationship was linear (Cai et al., 2017). Following Wang et al. (2014), the nonlinear wavenumber relation can be expanded by a Taylor series on the basis of its central frequency (x c ). Taking the first two terms of the series, we end up with the approximation for a linear dispersive wavenumber (k l ) of the waveform formulated as where c g is the group velocity. The resulting time domain trace, g l ðtÞ, corresponding to this approximation can be written as where G l ðxÞ ¼ Gðk l ðxÞÞ: If we retain only the second term in the Taylor series, we produce another approximation for a linear nondispersive wavenumber (k n ), with the resulting time domain trace, where G n ðxÞ ¼ Gðk n ðxÞÞ: Taking the analytical curve for the A0 mode, the central frequency of our selected wavelet, and evaluating Eqs. (2) FIG. 2. Experimental data transformed to frequency-wavenumber domain using the multianalysis of surface waves (MASW) method (Xia et al., 1999) and overlain with analytical dispersion and experimental curves.
and (5), we might expect intuitively that both relations are linear as plotted in Fig. 3 (blue and red curve). However, when converting those relations into frequencies and velocities, only k n gives constant group and phase velocities for all frequencies (see the red line and cross in Fig. 4), which results in a truly nondispersive waveform g n ðtÞ-only if those velocities are constant will dispersion be compensated. In contrast, k l results in waveforms g l ðtÞ, in which all events propagate with constant group velocity but their phases change inside their envelopes due to the nonconstant phase velocity as a function of frequency (see the dashed and solid blue lines in Fig. 4). This leads us to choose the second approximation, k n , for the dispersion compensation in the remainder of this paper.
Applying the algorithm to the data from our linear test setup, we first transform the data into the frequency domain [G d ðxÞ] as shown in Fig. 5 (top left) and use Eq. (1), together with the analytical dispersion curve, to obtain the wavenumber spectrum [G d ðkÞ] shown in Fig. 5 (top right). Using the k n approximation, we then obtain the resampled wavenumber spectrum [G d ðk n ðxÞÞ] and its corresponding frequency spectrum [G n ðxÞ] plotted in the bottom right and bottom left of Fig. 5, respectively. Finally, applying the inverse Fourier transform to resampled frequency spectrum [see Eq. (6)], we obtain the nondispersive signal as shown in Fig. 6 (bottom panel). Comparing the resulting trace to the original signal in the top panel of Fig. 6, it is clear that the proposed algorithm is able to compress the dispersion, making the identification of each event easier. For example, it is much easy to distinguish energy corresponding to the scattered waves and energy corresponding to the reflection from the nearest boundary. Indeed, not only are those events compressed, but all events reflected from other boundaries and their interaction with the scatterer are nicely recorded, compressed, and separated (e.g., the events after 400 ls).

IV. IDENTIFYING MULTIPLY SCATTERED WAVES
After finalizing the experimental setup and the dispersion compensation algorithm, we then perform our main measurements with the setup illustrated in Fig. 7. It consists of two scatterers that are separated by 15 cm, 4 sources (spacing 3 cm), and 49 receivers (spacing 0.03 cm). To increase the signal-to-noise ratio for each gather, the measurement is repeated 100 times and stacked (i.e., summed). Stacking a higher number of repeat experiments does not significantly improve the quality of the acquired signals (Seim, 2018).
Methods to identity multiply scattered waves from point diffractors in a common-source gather (CSG-data recorded on all receivers from a single or common energy source point) and common-receiver gather (CRG-data recorded on a single receiver from all energy source points) are explained in detail in Meles and Curtis (2014). By examining moveouts (kinematics) of waves, they identify patterns which correspond to either the first or the last scatterer of the wave scattering paths, the so-called fingerprints of each scatterer. The concept is that the first-order (primary) scattering of each individual scatterer has a unique moveout; energy from a multiple scattering path which has the same last scatterer (scatterer 1) will have the same moveout in a CSG as the primary but will be delayed in time. This is depicted in the CSG panel in Fig. 8, where the primary's moveout (light purple) is similar to the secondary or secondorder scattering (light dashed purple), but the latter is delayed due to the differences in the early path of the secondary as shown by the dark purple arrows before it arrives at the same last scatterer (note that "primary" and "secondary" are also commonly used as nouns to refer to an individual arriving wave of each type). This identifies the last scatterer of the path. The same shape observed in a FIG. 3. Comparison between the analytical (black), linear dispersive (k l ; blue), and linear nondispersive (k n ; red) curves of wavenumber with respect to frequency. Note that the linear dispersion approximation does not start from zero due to the constant term in the Taylor expansion of kðxÞ. This causes nonlinearity in the phase velocity.
FIG. 4. The result of calculating the group and phase velocities given the curves in Fig. 3. Only the nondispersive (k n ) wavenumber gives constant velocity for both phase and group velocities for all frequencies.
CRG identifies the first scatterer as depicted in the CRG panel where the primary (light gray) has the same shape as the secondary (dashed darker green), which shared the first scatterer (scatterer 1). Then, by combining information from both CSGs and CRGs, the travel time of higher-order scattering can be calculated as depicted in the lower section of Fig. 8, where the travel times of two secondaries from CSGs and CRGs are added; the travel time of a primary from the CSG is then subtracted from the result. This calculation must be done in a common trace that shared the source (S) and receiver (R).
Instead of identifying and extracting the multiply scattered energy as in Curtis (2014), L€ oer et al. (2015) introduced an automated scheme to identify fingerprints. They shift the observed fingerprints (moveouts) of primaries down the time axis of the gathers, and stack the energy (calculate the semblance) along that moveout. The arrival times of multiply scattered arrivals that share the  same fingerprint stand out as times with high semblance. Results of this process are illustrated in the semblance panel of Fig. 8, where a high amplitude appears when it detects a similar shape to the primary-in this case, the secondary which then gives the estimated primary-to-secondary delay time (DT). This scheme yields superior results on synthetic data.
When applying both the theory and scheme to our data, we face three main challenges. First, both the theory and the scheme are explained and tested using nondispersive data.
As we will show below, this challenge can be addressed by the above dispersion compensation algorithm. Second, our data only consist of CSGs instead of both CSGs and CRGs. Third, each of our gathers contains the signature of an additional scatterer that is not included in the theory, which, it turns out, is due to the physical source also acting as a scatterer. Here, we show how to address those challenges by adapting the theory so as still to allow the use of an automated scheme.

A. Practical aspects and improvements
As an initial step, we apply the above dispersion compensation method for all recorded data. Then, following L€ oer et al. (2015), the first step of the automated scheme is to isolate the primaries. This is done via a sequential process of cross correlating two common source gathers (A and B, say) to find the delay between the two that yields maximum similarity; this corresponds to all primary events from the same scatterer. For optimal results, the direct waves are muted so that they do not interfere with the cross-correlation process. The initial cross correlation is defined by Aðm; nÞB 1 ðm þ i; nÞ; where R  Illustration of CSG and CRG domains for the secondary and primary arrivals, and the resulting moveout in the CRG and CSG. The solid dark purple and green arrows denote the scattering paths that cause the delay of the primaries. The red squares illustrate sources and the one with annotation S illustrates the specific source used to select a common trace as for the receiver (the cyan triangle) with the annotation R. The semblance panel illustrates the semblance analysis for the primary from the CSG panel with the delay between primary and secondary denoted by DT. The bottom section shows how travel times of secondaries from the CSG and CRG settings and of the primary from the CSG settings are used to estimate travel times of a higher-order scattered wave (tertiary) in a common trace expressed by the path from a shared source (S) to receiver (R).
index in a CRG. Taking the maximum value of R 1 ðiÞ corresponds to time-shift i ¼ i 1 , which is used to shift gather B. The absolute value of the shifted gather is used to perform element-wise multiplication with gather A, yielding crossgather C 1 , Cross-gather C 1 is expected to have maximum values for the primary travel times, whereas other elements should be close to zero. However, if the data are noisy, additional steps need to be taken to improve the quality of the estimated primary by taking another gather (B 2 ) and then cross correlating it with the cross-gather C 1 . The general formula for additional cross correlations and element-wise multiplications is then and the cross-gather is calculated by where for j ¼ 0 we use gather A as C 0 . The travel time curve that corresponds to the primary in the cross-gather can then be picked by taking the time index of the maximum values in the cross-gather, Equation (12) is the final step required to obtain the travel time curve which corresponds to the primary of the first moveout. The moveout of this primary is then the fingerprint of a scatterer, which will also guide us to find other similar moveouts that correspond to higher-order scattering events. To isolate the second moveout, the first primary contained in the initial/main gather C 0 is muted prior to another initial cross correlation [Eq. (8)], followed by the same sequence of steps up to Eq. (12). The process can be applied iteratively to identify the rest of the primaries.
In our analysis, we add one more step after evaluating Eq. (12): we perform polynomial regression to the picked travel time curve since simply taking the maximum value of the final cross-gather might not correspond to the primary, especially in the case of low signal-to-noise ratio data or if there are primaries which cross kinematically as we have in our data. Furthermore, other criteria could be included in this step so that outliers are excluded. The next step is to calculate semblances in order to identify secondaries. A secondary represents a second-order scattering event. These must share similar first or last scatterers as particular primaries, but since they take a longer path, their fingerprint will be similar to the primaries but delayed in time. Semblance is calculated by taking the ratio of energies around particular kinematic moveouts (Yilmaz, 2001). In our case, we apply semblance methods to find energy with the kinematics of each primary in turn; the earliest wave in a CRG that shares the particular kinematic fingerprint as primary on the CSG must be the secondary between the two fingerprinted scatterers as illustrated in the CSG and CRG settings in Fig. 8.
In the original automated scheme, after primaries and secondaries in both CSG and CRG are determined, arrival times of higher-order scattering paths can be estimated by adding times of the secondaries from CSGs and CRGs and subtracting the travel time corresponding to one of the primaries. However, using similar concepts, we find that this also works for our data, which consist of CSGs only, as illustrated in Fig. 9. By adding times of secondary 21 (second-order scattering where the first scatterer is number 2 and the last scatterer is scatterer 1) with secondary 12 and subtracting times of primary 2, the travel time of a thirdorder scattering event denoted by 121 can be determined. This resolves the second challenge mentioned earlier.
For the third challenge, by taking the direct waves in place of a primary and using it to calculate the semblance in the above algorithm, we can obtain the secondaries corresponding to events scattered from the source. This is because the moveout of the direct waves is identical to all scattering paths where the source scatterer acts as the last scatterer as illustrated in Fig. 10.

V. RESULTS
The measured waveforms are first processed with the dispersion compensation algorithm (using the linear FIG. 9. Illustration of the method used to estimate the travel time of the first tertiary (121): add the travel times of secondaries 21 and 12, and then subtract the travel time of primary 2 from the result. nondispersive wavenumber approximation) and bandpassed with a Butterworth filter of 1.5-4 kHz. Applying the automated scheme to the compensated data from source number 2, we are able to isolate the cross-gathers and pick the travel time curves that represent the primaries as depicted in the left and middle panels of Fig. 11 with the semblances, including the one for the case with a source scatterer as shown in Fig. 12. Taking the shifted time of the first amplitude maximum on each semblance (excluding the initial peak) as the additional time needed for its secondary to arrive, we then plot the primaries together with the secondary in the right panel of Fig. 11. Note that primary 1 (dashed purple) has the same shape as secondary 1 (dashed red) and similarly for primary 2 (dashed green) and secondary 2 (dashed light blue), which confirms the theory. Additionally, the placement of the scatterers results in a crossing pattern of moveouts.
To validate our dispersion compensation results and assess the result of the automated scheme, we performed numerical modelling using the Foldy method implemented in the code of Galetti et al. (2013). This is a specific acoustic wavefield modelling code based on the multiple scattering theory formulated by Foldy (1945) and calculated in the frequency domain under the assumption of idealized isotropic point scatterers. Using the same geometry as for the experiment above, together with the estimated velocity, we model the scattering using a Ricker wavelet with a central frequency of 2.2 kHz. We choose this high frequency since this is the peak of the frequency spectrum that results from resampling the wavenumber from the linear nondispersive algorithm (see the top left and bottom left panels of Fig. 5). Note that we model the scattering using a different wavelet from the real data. Below we show, on the individual traces, that our observed and modelled traces have different events due to differences in source signature: our data have a broader source signature, which results in more interference in events.
We then compare the original data, the result of the dispersion compensation, and the results of the numerical modelling using two scatterers and three scatterers (the source also acting as one of the scatterers) in Fig. 13, superimposed with the primaries and secondaries estimated using the automated scheme, which show a very good match with the modelled gathers. The original gather is very dispersive, which makes events difficult to identify, but applying the compensation method makes events nicely identifiable FIG. 10. Illustration of the method used to estimate the travel time of a tertiary event that involves a source scatterer (1s2): add the travel times of the secondary that involves scatterer 1 and the source scatterer s denoted (1 s) with primary 2, denoted by (2), and then subtract the travel time of the direct waves (dw).
FIG. 11. The result of our automated scheme. (Left) The travel time curve for primary 1, (middle) the travel time curve for primay 2, and (right) the result of using those primaries to search for secondaries using semblance (Fig. 12).
before being subject to the fingerprinting theory. Looking at the gathers obtained from numerical modelling, the one obtained using three scatterers (including the source) is closer to the compensated gather than that for the case of two scatterers especially for the secondaries that involve the source scatterer (at 1 s and 2 s), which are marked by the red arrows. The compensated and modelled gathers are adaptively gained with time, which reveals artifacts that are almost horizontal (starting from about 200 ls). Those artifacts are the cause of the high amplitudes in the direct wave semblances shown in the right panel of Fig. 12.
Given the primaries and secondaries, we then estimate the travel time for higher-order scattered waves, which we limit to fourth-order scattering (by using CSG only) in this research. Plotting the result of the estimation on the semblances (all with zero time moved temporarily to the travel time of their respective primaries), we can see that only primaries (at time 0) and secondaries have considerably higher amplitude as seen in Fig. 12. When each semblance of primaries 1 and 2 (left and middle panels) has only one obvious secondary, the direct wave's semblance (right panel) will give rise to two clear peaks of secondaries due to secondorder scattering that involves the source with each of the other scatterers (see Fig. 10), which are observed as 1 s and 2 s (red arrows in Fig. 13). The high amplitudes after 100 ls are due to artifacts from the compensation process that have FIG. 12. The semblances obtained from primary 1 (left), primary 2 (middle), and the direct wave (right). Each semblance is also marked with the times of secondaries, together with the higher-order scattering path arrival time that comes from the same last scatterer (time axis is zeroed at the arrival times of each respective primary). See the main text for details of the notation used in the key.
FIG. 13. (Left to right) The original gather, original gather after dispersion compensation, modelled gather using three scatterers, and the modelled gather using two scatterers. The red arrows mark the secondary events involving the source scatterer. horizontal apparent moveout almost identical to that of the direct waves. For clarity, to read the legend in Fig. 12, the first letters P, S, T, and Q are primary, secondary, tertiary, and quaternary without interaction with the source scatterer, respectively, whereas the remaining number indicates the last scatterer visited by the scattering path before it gets recorded at the receiver. For example, T1 means that it is a tertiary or third-order scattered wave with the scattering path from the source to scatterer 1, then 2, then back to 1, before finally being recorded by the receivers (121). On the other hand, the rest (i.e., without letters P, S, T, and Q) denote scattering events that involve the source scatterer. For example, 1s2s means that it is a fourth-order scattering event where the scattering path is from the source to scatterer 1, back to the source (acting as a scatterer), then to scatterer 2, before finally hitting the source again as the last scatterer and being recorded by the receivers.
Additionally, since the coordinates of all components (sources, scatterers, and receivers) as well as the estimated constant group velocity are known after applying the dispersion compensation algorithm, and given that the central frequency of our wavelet is 1.5 kHz means that the waves travel with constant phase and group velocities of approximately 2310 m/s (see Fig. 4), we are able to estimate all arrival times by assuming straight ray paths for all scattering paths (ray theory).
To obtain more detailed comparisons, we plot a single trace taken from receiver five from gathers shown in Figs. 13 and 14 from both the observed/recorded compensated and modelled gathers, for the modelled case using two scatterers (top) and three scatterers (bottom), superimposed with the arrival time predictions (from both straight ray and fingerprinting theory).
The top panel shows that the modelled trace/seismogram (red curve) matches with several unique events in the compensated trace (black), together with the estimations from the fingerprinting theory (black, purple, and blue upside-down triangles) and ray theory (red upside-down triangles), where all estimates and the modelled trace only take into account two scatterers present in the medium (ignoring the source scatterer). However, as also indicated in Fig. 13, the modelled trace and the scattering path estimates do not account for the full information especially in relation to all of the scattering paths that involve the source acting as a scatterer.
Adding the source scatterer to the numerical modelling (see bottom panel), the trace fits the compensated trace better compared to the result using two scatterers. Travel time estimates that account for the source scatterer (light blue and green triangles for estimates using fingerprinting and red for straight ray estimates) match the events in both traces better compared to those in the top panel. Note that all of the observed and modelled gathers have been adaptively gained such that the scattering events can be visually inspected and compared to the modelling; this causes the amplitude to be exaggerated at later times. Furthermore, since the effective wavelets (source signature) of our experiment (sine wavelets) are considerably broader than the modelled wavelet, we find more interference in our observed trace, which gives less time separation compared with both modelled traces. However, the general trends of the FIG. 14. A comparison of (top) a trace taken from the compensated data and numerical modelling using two scatterers and (bottom) a trace taken from compensated data and numerical modelling using three scatterers. All traces are taken from receiver number 5 of the gathers displayed in Fig. 13 and superimposed with the travel time predictions (assuming straight rays and using fingerprinting theory). The gain is adaptively increased at later times, amplifying later-arriving energy. Note that for the modelled traces, having three scatterers in our modelling results in more events compared to using only two scatterers. amplitude and event timing seem to be well matched between the observed and modelled traces especially when compared with the results of both travel time predictions and, particularly, up to 200 ls when we still identify some of the fourth-order scattering. Beyond that, we find it difficult to distinguish the remaining predicted events in both the observed and modelled traces.

VI. DISCUSSION
We have shown that fingerprinting theory is applicable to predict scattering paths and arrival times of higher-order scattering in practical laboratory experiments. This inspires future possible applications, such as better localization of scatterers, because waves leaving each scatterer on a wide distribution of azimuths may be analysed. Further, the information from multiple scattering may be useful to generate new algorithms for nonlinear seismic migration or imaging (Halliday and Curtis, 2010) that usually rely on single scattering theory (Smitha et al., 2016). The dispersion compensation method may also be used to process data from dispersive energy, such as surface waves, in new ways given estimates of their dispersion curves (Cao et al., 2020), for example, to separate surface wave modes in order to measure their phase velocities (Zhang et al., 2020). To add complexity to the experiments, we might use materials that simulate layered structures to include reflections from layer interfaces rather than only from the boundaries. All of these avenues are open questions for future research to explore.

VII. CONCLUSION
We have shown that it is possible to perform scattering experiments in strongly dispersive media, in this case, using a thin aluminum plate and to apply fingerprinting theory to identify multiply scattered wavepaths. Prior to fingerprinting, dispersion compensation is necessary such that the theory can be applied, requiring analytical or experimental dispersion curves. Using the automated scheme, we isolated all primaries and predict secondaries and higher-order scattering up to quaternaries. This prediction is compared with the results of calculating arrival times from straight ray paths which were found to match. We added a case when a source also acts as one of the scatterers to the theory and algorithms. The result when including this additional scatterer enables us to estimate more scattering paths and yields results that better match the numerical modelling for the data, which adds significant confidence in our application of fingerprinting theory. We expect both dispersion compensation, as well as fingerprinting, to act as significant enablers for future wave experiments.