Closed-aperture unbounded acoustics experimentation using multidimensional deconvolution

Acoustic wave propagation in air or water can be studied physically through laboratory experiments such as in air-filled cavities or water tanks (e.g., Cassereau and Fink, 1992; Fink, 1992). However, waves reflected from the boundaries of the experimental setup often cause an adverse effect, particularly for experiments that aim to study only the scattering within the experimental domain. The undesired boundary reflections can interfere with, or mask, the wavefields related to the interior scatterers. Ideally, the acoustic waves studied in such laboratory experiments should propagate in a boundless domain of infinite extent, which is impossible to realize experimentally (Larose et al., 2010; McDonald et al., 1983; Mo et al., 2015).

One common solution to the problem of boundary reflections in laboratory experiments is to deploy absorbing materials around the experimental setup such that an anechoic chamber is formed, absorbing the outward propagating acoustic waves (e.g., Beranek and Sleeper, 1946). The materials acting as absorbers can be pyramid-shaped wedges of fiberglass (Munjal, 2002; Trinh et al., 2019), for example, or self-similar acoustic meta-materials (Shao et al., 2019; Zhang et al., 2020). However, deploying such passive absorbers tends to be impractical for experiments carried out in the low-to-mid frequency range (e.g., 1–10 kHz), since a thick (order of meters) and costly layer of the absorber is required to enclose the experiments. Also, in practice, low-frequency outward propagating waves are not perfectly absorbed using this approach. Another solution is to build an active anechoic chamber where acoustic sources are densely deployed around the boundary of the experimental domain to cancel the outward propagating waves (Beyene and Burdisso, 1997; Guicking and Lorenz, 1984; Habault et al., 2017; Smith et al., 1999). These active sources, sometimes installed together with passive absorbers, can effectively render the boundary transparent for low-frequency experiments (<1 kHz) but are not efficient for high-frequency experiments. A recently developed technology called immersive wave experimentation can achieve broadband absorption of outward propagating waves (e.g., 1–20 kHz) (Becker et al., 2018, 2020; Vasmel et al., 2013). Simultaneously with the physical experiments, waves measured inside the experimental domain are extrapolated outward to the domain boundary ahead of time such that sources densely deployed on the boundary can cancel the outward propagating waves based on the prediction.

Instead of using the above-mentioned methods that attempt to change the boundary condition of the experimental setup during laboratory experiments, the problem caused by boundary reflections can also be solved by removing the scattering component related to the domain boundary from the recorded data afterwards (i.e., in post-processing). Conventionally, the arrivals or waveforms involving boundary reflections in the recorded data can be identified and removed using a time windowing approach (e.g., Blum et al., 2011). However, the determination of a suitable time window highly depends on separating the boundary-related waveforms from the waveforms related to the interior scattering in terms of the arrival times. Hence, the applicability of this technique is often restricted to high-frequency (e.g., MHz range) experiments. The post-processing of data recorded in the laboratory such that the boundary-related scattering is removed can also be achieved by means of multidimensional deconvolution (MDD), also called Rayleigh–Betti deconvolution, widely used in exploration seismology (Amundsen, 2001; Amundsen et al., 2001; Wapenaar et al., 2011). The seismic data recorded in the subsurface (e.g., in a horizontal well), or on the ocean bottom, can be post-processed such that the scattered waveforms related to the overlying medium, or a reflecting top surface,¹ are removed from the recorded data. One key feature of the MDD method, as used to remove the overlying-medium- or surface-related multiples, is that neither the medium below the subsurface receiver array nor the source wavelet is required to be known (Amundsen, 2001). Similarly, in a laboratory experimental setup, the material properties of the scattering medium (e.g., a rock sample) inside a closed-surface receiver array do not need to be known. Rakotonarivo et al. (2013) propose an MDD method that can be applied to laboratory experiments such that the structural impedance [i.e., the in vacuo response due to excitation forces, see Zhou et al. (2019)] of a physical object can be estimated from a random noise field. Williams et al. (2017) apply this method to real laboratory experiments to compute the structural impedance of an elastic object placed in a non-anechoic room with active sources randomly distributed outside the object. Following the same principles, Sternini et al. (2019) construct a laboratory facility to determine the scattering properties of an arbitrary object from the measurement of responses at the surface of the object, which is placed in a noise field. The scattering imprint caused by the boundary of the experimental domain is removed from data recorded at the surface of the object, and one obtains the structural impedance matrix as the MDD result. These laboratory experiments are similar to the seismic MDD applications discussed by Poletto et al. (2014) and Weemstra et al. (2017b) where a reflecting boundary condition or virtual reflector is applied at the level of the subsurface receiver array replacing all the scatterers above, including the top free surface.

In this paper, we focus on post-processing of a dataset acquired in an air-filled waveguide using an MDD method such that the scattering effects of the rigid boundary surrounding the experimental setup are completely removed while only the Green's functions associated with the scattering object inside the experimental volume are obtained in the absence of an exterior reflecting boundary. The acquisition geometry involves a receiver array, completely enclosing an experimental volume, which enables full-aperture sampling of the scattering related to the medium inside the receiver array. Also, this MDD method requires wavefield decomposition such that the in-going wavefield can be derived from the total wavefield recorded at the receiver array. Hence, we propose a wavefield separation method that enables the wavefields recorded around general, curved, closed receiver arrays to be separated into in- and out-going components. The general methodology proposed in this paper based on MDD thus has the potential to replace the current paradigm of absorption- and active cancellation-based removal of boundary reflections. The new paradigm involves proper sampling of boundary reflections combined with MDD to enable fully enclosing reconstruction of scattering wavefields for true radiation boundary conditions.

In Sec. II, we present the MDD theory and propose the wavefield separation method. In Sec. III, we demonstrate the accuracy of the wavefield decomposition method with a numerical example using the forward, multidimensional convolution (MDC) relationship derived in Sec. II. We then demonstrate the ability of the MDD method to remove the scattering related to the boundary of a physical two-dimensional (2D) waveguide and further confirm the results (i.e., obtained Green's functions) with numerical simulations. In Sec. IV, we discuss the general applicability of MDD in acoustic wave propagation experiments and compare the MDD method to immersive wave experimentation. Section V summarizes our conclusions.

We consider a physical experiment bounded by a rigid boundary and a corresponding desired experiment without the rigid boundary, as shown in Figs. 1(a) and 1(b). The two experiments are connected via the convolution-type integral of the two-way representation theorem (Ravasi et al., 2015):

where all the quantities are expressed in the frequency domain (denoted by the symbol $̂$⁠), and ω is the angular frequency. In Eq. (1), the wave quantities $p̂(xir)$ and $p̂(xr)$ are pressures recorded at the interior receiver (⁠$xir$⁠) and the receiver surface $S(xr)$ (with normal vector component n_i) in the physical experiment, while $v̂i(xr)$ is particle velocity recorded at $S(xr)$ in the x_i direction. The Green's functions are both associated with the desired medium in which $Giv,q(xr,xir,ω)$ represents the ith component of the particle velocity recording at $xr$ due to an impulsive point source of volume injection rate density at $xir$ (superscript q for monopole), and $Gp,q(xr,xir,ω)$ represents the pressure recording due to an impulsive point source of volume injection rate density. The detailed derivation of Eq. (1) is given in the  Appendix.

The medium properties inside the volume V are the same in the physical and desired experiments, as shown in Fig. 1. Outside V, the medium for the desired experiment is homogeneous [i.e., Fig. 1(b)], namely such that the surface S effectively satisfies a radiation boundary condition so that waves propagating out of S are never scattered back into V. Since we are interested in an expression for a single type of Green's function, we apply the chosen medium properties and simplify the representation integral in Eq. (1) as follows. First, we decompose the wavefields associated with the physical experiment on the right-hand side of Eq. (1) into in- and out-going components: $p̂=p̂in+p̂out$ and $v̂i=v̂iin+v̂iout$⁠, as shown in Fig. 1(c). Hence, Eq. (1) is recast as

Note that the Green's functions in Eq. (2) can also be decomposed into in- and out-going components (i.e., $G=Gin+ Gout$⁠), but due to the radiation boundary condition, the in-going component, $Gin$⁠, is always zero, and hence $G=Gout$⁠. The two purely out-going terms $p̂outĜiv,q$ and $v̂ioutĜp,q$ cancel each other, since their stationary points yield exactly opposite contributions to the surface integral (Ravasi et al., 2015; Wapenaar and Berkhout, 1989; Wapenaar and Fokkema, 2006). Hence, Eq. (2) is simplified:

It can be shown, as in Ravasi et al. (2015) and Wapenaar et al. (2011), that the two terms $p̂in(xr)Ĝiv,q$ and $−v̂iin(xr)Ĝp,q$ in Eq. (3) yield the same contribution to the surface integral. Here, we choose to keep the term $p̂in(xr)Ĝiv,q$ and transform the simplified equation back to the time domain:

where the symbol $*$ refers to convolution in time t. In Eq. (4), the Green's function $Giv,q(xr,t| xir,0) ni$ denotes the normal particle velocity at the surface $S(xr)$ for an impulsive volume source injection at $xir$⁠. We further apply source-receiver reciprocity to Eq. (4), i.e., $Giv,q(xr,t| xir,0) ni=−G,ip,f(xir,t| xr,0) ni$ (Fokkema and van den Berg, 2013), where $[.],i,f$ denotes the impulse response due to a body force source, f (i.e., a dipole), in the $xi$ direction. For convenience, the Green's function $G,ip,f ni$ can be replaced by $G¯d$ for a normally oriented body force source and a pressure wavefield (i.e., dipolar Green's function). Hence, we finally obtain

which constitutes a MDC relationship between the wavefields in the physical and desired experiments.

Note that the dipolar Green's functions on the right-hand side of Eq. (5) are defined in the desired, unbounded experiment. Hence, the Green's functions are to be treated as unknowns that we aim to solve for using knowledge of wavefield quantities such as $p(xir,t)$ recorded in a physical experiment and $pin(xr,t)$ separated from the recorded data. However, Eq. (5) is a Fredholm integral of the first kind (Aster et al., 2013), the inversion of which is highly ill-posed if wavefields p and $pin$ are considered for only a single source outside the receiver surface (Wapenaar et al., 2008). Hence, one needs to excite many sources outside S, record and separate the data for each of these sources individually and compose a system of equations by considering multiple versions of Eq. (5) simultaneously. This system of equations, denoted as the $ptotal−pin$ scheme in the following, can then be inverted using MDD (Amundsen, 2001; van der Neut and Herrmann, 2013; Wapenaar et al., 2008).

The sought-after Green's functions in Eq. (5) correspond to the impulsive responses of a desired medium that only includes the physical scatterers of the real medium inside the surface S; the exterior scatterers, such as the rigid boundary of the experimental domain, do not exist in the desired medium. Hence, MDD can be used to remove undesired scattering effects caused by the physical boundary in laboratory experiments. The recorded data are (post-)processed such that only waveforms related to the interior scatterers of interest are kept. Furthermore, note that the source that generates the waves in the physical experiment remains undefined in the theory leading up to Eq. (5). Therefore, the sources outside the receiver surface S do not need to be densely or regularly spaced and can be put sparsely, without consideration of the Nyquist sampling criterion (Weemstra et al., 2017b). In fact, they can even be a noise wavefield (e.g., Sternini et al., 2019). Thus, the source spectra or temporal signatures do not need to be known in the physical experiment, and, in fact, any source characteristics, including radiation patterns, will be removed by MDD. Finally, because MDD is based on a convolutional representation, as opposed to a correlational representation, MDD works for dissipative media and for different types of wave physics inside V (e.g., elastic and viscoelastic) (Wapenaar et al., 2011).

We discuss the practical aspects of solving a system of equations constructed from multiple equations of the form as in Eq. (5) later (in Sec. III B). First, we discuss how MDD can be applied to a more typical laboratory configuration for a complete, full-aperture sampling of the scattering related to the medium inside V.

Equation (5), giving the physical wavefield in the interior in terms of a MDC of the in-going part of the physical wavefield on the surface S and the dipolar Green's functions for the same medium, but with radiation boundary conditions, is valid for any point x_ir inside S. Having this ability, through MDD, to remove the boundary-related scattering events from recorded data at receivers inside S seems very useful. However, the requirement of interior receivers, i.e., within the recording surface, is impractical as such interior receivers impede placing large solid objects there. We therefore propose to apply Eq. (5) to the configuration shown in Fig. 2(a) in the limit where the receivers, ${xir}$⁠, approaching the surface from the interior, lie on S^rec (i.e., on the inner blue circle). However, note that an integrable singularity exists for the dipolar Green's function on the right-hand side of Eq. (5), and pole artefacts may be problematic in the inverted Green's functions at and around the location where the source and receiver positions are co-located (i.e., when $xr=xir$⁠). Also, the Green's functions $G¯d$ in Eq. (5) contain the direct arrivals between all locations on the recording surface S^rec, which are often not of interest when studying the scattering properties of the interior object in the laboratory experiments (e.g., when they bypass the scatterers completely). The pole artefacts and direct arrivals can be removed from the MDD result obtained from a laboratory experiment with the scattering object placed inside the recording surface S^rec (“heterogeneous” case) by carrying out a second experiment with a homogeneous experimental domain (“homogeneous” case) and subtracting the MDD result from that in the heterogeneous case. Thus, we carry out two laboratory MDD experiments with and without the interior scattering object for two sets of Green's functions, $G¯dhet$ and $G¯dhom$⁠, and we finally obtain

where $G¯dscat$ only contains the waveforms that are related to the scattering caused by the interior object of interest. The two-layer recording surface S^rec shown in Fig. 2(a) allows for wavefield decomposition, which will be discussed later, and has nothing to do with the singularity of the Green's functions in Eq. (5), which is considered when the two blue circles in Fig. 2(a) coincide at the inner receiver circle, i.e., assuming that the wavefield decomposition can be done for data recorded along a one-layer receiver array. Note that sources as shown in Fig. 2(a), in the MDD theory, are not necessarily located at the rigid boundary surrounding the experimental domain and can be placed anywhere between the rigid boundary and the recording surface S^rec.

The geometry shown in Fig. 2(a) is similar to the acquisition geometry for ocean-bottom data processed using MDD in marine seismic applications (Amundsen, 2001; Wapenaar et al., 2011). However, compared to those applications, with almost exclusively one-sided acquisition geometries (Amundsen et al., 2001; Bellezza and Poletto, 2014; Minato et al., 2011; van Dalen et al., 2014; Weemstra et al., 2017a), a physical experiment in the laboratory such as the one we propose in Fig. 2(a) has the advantage of a fully enclosing acquisition geometry. The sources located (sparsely) outside the recording surface (in our case on the rigid boundary) can illuminate the scatterers inside S^rec from all directions, which alleviates the problem of rank deficiency in the inversion of the MDC equations, particularly compared to experiments with a one-sided acquisition geometry. Nevertheless, the inverse problem is still ill-posed and needs to be regularized (Rakotonarivo et al., 2013). Rank deficiency is also affected by the number of sources and the frequency bandwidths (Ruigrok et al., 2010; Hunziker et al., 2010), but a further discussion is beyond the scope of this paper. Note that when performing MDD for the configuration shown in Fig. 2(a), since Green's functions between every pair of points on the recording surface are obtained, not only are all the boundary-related scattering events removed from the data, the Green's functions enable full-aperture reconstruction of the scattering wavefield related to the object(s) under investigation.

The MDD theory, given in Sec. II A, requires wavefield separation of data recorded along closed, circular receiver arrays as in the experimental setup shown in Fig. 2(a). However, as pointed out in Ravasi et al. (2015), wavefield separation along a closed-aperture receiver array is still a challenge in MDD applications. The commonly used plane-wave decomposition method for data recorded on a planar receiver array (e.g., Amundsen, 1993; Day et al., 2013; Fokkema and van den Berg, 2013; Ramírez and Weglein, 2009; Wapenaar and Berkhout, 1989) does not work for data recorded along a curved surface. Wavefield separation can be achieved using the technique of cylindrical harmonics (Hulsebos, 2004), which is restricted to circular receiver arrays. Also, wavefield injection in numerical modeling can be used to separate wavefields along arbitrary curved surfaces (e.g., Amundsen and Robertsson, 2014; Robertsson et al., 2015). However, for data recorded in a closed aperture, wavefield injection cannot decompose the data into in- and out-going components, which are often obtained as up- and down-going wavefields for data recorded in experiments with a one-sided acquisition geometry (see Thomsen et al., 2018; Thomson, 2012).

In this paper, we propose a new wavefield separation method for curved, closed surfaces, inspired by a seismic separation method for waves recorded on curved but open surfaces developed in Riyanti et al. (2008) and Ferber and van Manen (2017). The method for open surfaces involves decomposing the data recorded in a marine seismic survey on a non-planar receiver array. The method relies on setting up a linear relation between the full (up-going + down-going) pressure data (in the space-frequency domain) recorded at an array with an arbitrary depth profile and the desired, up-going pressure data (in the wavenumber-frequency domain) that would have been recorded on a flat, reference surface. The basic idea behind this method is that it is easier to formulate the (inverse) transform from a regularly sampled (plane-wave) spectrum in the wavenumber-frequency domain to the irregularly sampled data in the space-frequency domain than vice versa. Therefore, a least-squares formulation is intuitive: the linear relation is inverted to estimate the regularly sampled up-going wave spectrum on the flat reference surface. Finally, from the regularly sampled spectrum, the up-going wavefield can be obtained at any desired regular or irregular location in space.

Consider a group of neighboring receivers on the inner layer of the recording surface $Sr1$⁠, as shown in Fig. 3(a), which are indexed clockwise as $1,2,…,N$⁠. For this group of receivers, we construct a local coordinate system $x′,z′$ with the x$′$ axis tangential to the inner circle $Sr1$ and z$′$ axis pointing in the normal (radial) direction. The origin in the local coordinate system coincides with the central receiver of this receiver group. Figure 3(b) shows the selected receivers in this local coordinate system, as well as the corresponding receivers in the outer layer of the recording surface $Sr2$ (i.e., with the same angular range). In the local coordinate system $x′,z′$ as shown in Fig. 3, the locations of the selected receivers at the inner and outer circles are denoted as $(xnr1,znr1)$ and $(xnr2,znr2)$⁠, and we treat the x axis or $z′=0$ as the construction surface of the separated wavefields.

In the pressure wavefield separation method of Ferber and van Manen (2017), the construction surface is located at the sea surface, and the linear operator also contains a model for wave reflection at the sea surface, relating the up- and down-going wavefields below the sea surface, halving the number of unknowns. For wavefield separation in the laboratory configurations such as shown in Fig. 2(a), there is no simple, general model for reflection at the boundaries of the experimental domain, and thus no such relation exists between the in- and out-going waves. Therefore, either pressure data at two offset recording surfaces are needed to invert for the individual in- and out-going wavefields, or both the pressure and normal particle velocity are needed on a single surface. We process a subset of data corresponding to the receivers at $(xnr1,znr1)$ and $(xnr2,znr2)$ in the local coordinate system to separate the wavefields at the construction surface $z′=0$⁠. In total, the number of recorded traces in this subset thus is 2N. Because only the central receiver located in the inner circle $Sr1$ is also part of the wavefields inverted at the construction surface, we only keep the needed separated (in-going) wavefield at the central receiver. The same process is then repeated for all the receivers at $Sr1$ as the “central receiver.”

The pressure data recorded at the selected receivers in the local coordinate system $x′,z′$⁠, as shown in Fig. 3(b), are denoted as $p(xnr1,znr1,t)$ and $p(xnr2,znr2,t)$ (⁠$n=1,2,…,N$⁠). We first transform the time-domain pressure data into the frequency (f) domain, giving $p̂(xnr1,znr1,f)$ and $p̂(xnr1,znr1,f)$ (⁠$n=1,2,…,N$⁠) with $f=ω/2π$⁠. In this paper, we define the Fourier transform of a time-dependent wavefield quantity, e.g., $p(x,t)$⁠, as $p̂(x,ω)=∫−∞+∞ exp (−jωt)p(x,t) dt$⁠, where j is the imaginary unit. We further define the desired separated in- and out-going components in the wavenumber-frequency domain as $s̃in(f,kn)$ and $s̃out(f,kn)$⁠, where k_n is the regularly sampled horizontal wavenumber (still with $n∈1,…,N$⁠). Note that the wavenumber spacing is chosen in accordance with the spatial Nyquist sampling criterion and that only wavenumbers corresponding to propagating waves are considered (i.e., $kn<f/c$ with c as acoustic wave speed), implying an increasing number of unknowns (⁠$≤N$⁠) with frequency. For each frequency, we express the total pressure and the in- and out-going components at the construction surface as two column vectors,

Then, for each frequency, we construct a Hermitian matrix, i.e., the Gram matrix $L(f)$⁠, as

with the submatrices

and $kz(km)=(f/c)2−km2$ is the vertical wavenumber (Ferber and van Manen, 2017). The Gram matrix L relates the out- and in-going components at the construction surface, i.e., $ssep=[sout,sin]T$⁠, to the total pressure waves, i.e., $ptot=[pr1,pr2]T$⁠, through a linear system,

and its block form for each frequency

As mentioned previously, the total pressure $ptot$ is in the space-frequency domain with data points $(xnr1,2,znr1,2)$ irregularly sampled along $x′$ and $z′$ axes in the local coordinate system, while the separated in- and out-going components $ssep$ are in the wavenumber-frequency domain with regularly sampled wavenumbers.

Equation (10) can be solved for $ssep$ using the method of least squares (Aster et al., 2013),

where the superscript H denotes conjugate transpose, I is an identity matrix, and ϵ is a small value compared to the largest eigenvalue of the square matrix $LHL$⁠. The regularization term $ϵI$ helps alleviate the problem of rank deficiency when inverting the matrix $LHL$⁠. The obtained in- and out-going components $ssep=[sout,sin]T$ for each frequency are still in the regularly sampled wavenumber domain, and for the desired separated wavefields (i.e., $pr1in$⁠) that are sampled irregularly along $x′$ and $z′$ axes in the space domain, we can reuse the individual Gram submatrix to do the conversion (Duijndam et al., 1999):

The separated wavefield $pr1in$ is taken at the construction surface that is tangent to the inner receiver array $Sr1$ and finally transformed back to the time domain. Finally, the recovered wavefields are band limited to the bandwidth of the source wavefield using a bandpass filter.

Figure 4 summarizes the steps in the wavefield decomposition algorithm used to process both synthetic and laboratory data. In practice, the separated wavefields along the recording surface S^rec contain small artefacts (i.e., presence of unphysical events) toward the beginning and end of time traces, which are suppressed using tapering.

Since wavefield separation is key to the proposed change of boundary conditions using MDD, we first demonstrate the performance of the wavefield decomposition method using the MDC relation given in Eq. (5). In a numerical simulation, we consider the configuration shown in Fig. 5. Both sides of Eq. (5) can be modeled and computed such that their comparison can demonstrate the correctness of the wavefield separation. This synthetic example was simulated with spectral-element modeling using the software salvus (Afanasiev et al., 2019), and the parameters used are summarized in Table I.

To generate $p(xir,t)$ and $pin(xr,t)$ in Eq. (5) in numerical simulations, we put a pressure source at the rigid boundary. The source signature corresponds to a so-called Ricker wavelet (i.e., the second derivative of a Gaussian), q(t), with peak frequency f_p = 2000 Hz (Wang, 2015). The simulations of the reference Green's functions in Eq. (5) involve the replacement of the rigid boundary with an absorbing boundary condition (Kosloff and Kosloff, 1986) in the synthetic configuration shown in Fig. 2(a). Since the Green's function in two dimensions has an infinitely long tail, we convolve the simulated Green's functions for visualization purposes with the same Ricker wavelet (f_p = 2000 Hz) used in the simulations [Fig. 6(a)]. Hence, the MDC equation can be recast as

The wavefield separation algorithm proposed in Sec. II B was applied to the simulated wavefields $ptot(xr,t)$ recorded at the recording surface S^rec to obtain the in-going component $pin(xr,t)$⁠. The size of the selected receiver group was N = 13, and the regularization factor ϵ used in Eq. (12) is equal to 0.001 multiplied with the largest eigenvalue of the square matrix $LHL$ per frequency. The values of N and ϵ are chosen by trial and error using the evaluation of the MDC relation [Eq. (14)] as an effective tool to assess the quality of the wavefield separation. Also, a bandpass filter (i.e., a Butterworth filter) between 200 and 5500 Hz was applied to the separated wavefield.

Figure 6 compares the computed left-hand side of Eq. (14) to the right-hand side, and their good agreement demonstrates the performance of the wavefield separation method for closed, curved surfaces.

Figure 2(b) shows a 2D waveguide used for laboratory MDD experiments, corresponding to the schematic shown in Fig. 2(a). The key properties of the physical experiments are summarized in Table I. The air-filled waveguide consists of two parallel polymethyl methacrylate (PMMA) plates spaced $da=$ 2.5 cm apart. Waves propagating in the waveguide between the two plates can be considered 2D as long as the frequency content of the waves is below the cutoff frequency of the fundamental mode, $fc=0.5c0da−1$ (Redwood, 1960). With sound speed $c0=347$ m/s, the cutoff frequency for the waveguide shown in Fig. 2(b) is $fc=6.9$ kHz. The rigid circular outer boundary of the waveguide consists of synthetic resin and is placed between the two plates. The loudspeakers (i.e., sources) are installed on the circular boundary, while the pressure-sensitive microphones (i.e., receivers) are mounted flush with the inside of the PMMA plates. Note that the recording surface S^rec comprises two circular arrays, as required for wavefield separation, each consisting of 114 electret condenser microphones. To demonstrate our approach, a physical object of interest, i.e., a steel block with dimensions 30 cm × 6 cm × 2.5 cm, is placed at the centre of the waveguide. This steel block can be regarded as a perfectly rigid scatterer due to the large contrast in acoustic impedance (ρc) between air and steel.

In the MDD experiment, each source on the circular boundary of the waveguide is excited sequentially. The source signature corresponds to a Ricker wavelet with peak frequency f_p = 2 kHz. For each source, the experiments are repeated 30 times and averaged to improve the signal-to-noise ratio of recorded data. Figures 7(a) and 7(d) illustrate the two experiments (for one source location) without and with the steel scatterer (i.e., the homogeneous and heterogeneous cases). Figures 7(b) and 7(e) show the corresponding recorded data, which contain a large number of boundary reflections. Due to the strong dominance of the direct wave and the multiple reflection from the rigid boundary, the imprint of the physical scatterer is only barely visible when inspecting the two datasets [Figs. 7(b) and 7(e)]. The recorded data appear somewhat patchy due to some systematic, low-frequency noise introduced by the power supply units (PSUs). As multiple receivers are connected to the same PSUs, all these receivers show the same systematic noise. Note that the actual source signatures for the sources outside the recording surface [see Fig. 2(b)] are slightly different from one another due to variability in manufacturing and mounting of the sources. However, the MDD method works very well also in the presence of such unknown source characteristics.

The data recorded on the surface S^rec are separated to isolate the in-going part for both the experiments without and with the steel scatterer [Figs. 7(c) and 7(f)]. In the wavefield separation, both the regularization factor ϵ used in Eq. (12) and the bandpass filter for the separated wavefields follow the choices in Sec. III A.

The system of equations assembled from Eq. (5) [i.e., $ptotal−pin$ scheme] can be solved as an inverse problem in the Fourier domain (Amundsen, 2001; Holvik and Amundsen, 2005; Poletto et al., 2014; Ravasi et al., 2015; Wapenaar et al., 2008; Weemstra et al., 2017b) or in other domains such as the curvelet domain (van der Neut and Herrmann, 2013). In this paper, we solve the system of equations using an iterative damped least-squares method (Paige and Saunders, 1982), which carries out forward MDC in the frequency domain to minimize a regularized $ℓ2$ normal cost function (e.g., $minx||Ax−b||2+d2||x||2$⁠) evaluated in the time domain (i.e., PyLops MDD solver, see Ravasi and Vasconcelos, 2020). The damping parameter was set to $d=10−6$⁠, and 20 iterations were used to invert for the Green's functions.

Figures 7(h) and 7(k) show examples of the inverted Green's functions, which are convolved with a Ricker wavelet with peak frequency f_p = 2 kHz for visualization purposes. These Green's functions correspond to the configurations shown in Figs. 7(g) and 7(j), where the rigid boundary of the waveguide is removed, and the dipolar sources of the Green's functions are located at the recording surface S^rec. Note that the inverted Green's functions contain undesired direct arrivals as well as pole artefacts caused by the integrable singularity that exists for co-located source and receiver positions of the Green functions. Hence, we carry out the subtraction as in Eq. (6) to obtain the scattering Green's functions that are only related to the steel object, as shown in Fig. 7(i). This result is compared to a reference MDD solution (with amplitude normalization) obtained from a synthetic experiment [Fig. 7(l)]. Note that in the laboratory, we cannot obtain a true reference solution, since this would require removing the rigid boundary and carrying out the physical experiment in a largely extended waveguide with a radius about 6 m. Instead, we will have a synthetic reference solution described later.

Comparing Figs. 7(h) and 7(e) shows that almost all the boundary reflections are removed. We further focus on the waveforms in the obtained Green's functions. Figure 8 shows a zoom-in of the laboratory MDD results. In the heterogeneous MDD result including the steel scatterer, the direct arrival and pole artefacts mask the imprint of the scatterer [Fig. 8(a)]. The MDD result obtained in a homogeneous experimental domain [Fig. 8(c)] additionally shows a non-acoustic imprint, which is likely caused by the systematic acquisition noise present in the acquired data [see Fig. 7(b)]. The power-supply-related noise causes systematic energy on neighboring traces that is not related to acoustic arrivals. This energy interacts with physical acoustic arrivals through the MDD process and maps to coherent-looking arrivals in the inverted Green's functions. [i.e., non-acoustic imprint in Figs. 8(b) and 8(c)]. This noise is in the same frequency range as the experiments and hence cannot be suppressed by bandpass filtering. As the noise does not linearly generate artefacts in the MDD process, its imprints in the inverted Green's functions do not cancel out during the subtraction [i.e., Eq. (6)], and some undesired components exist in the obtained scattering Green's function [Fig. 8(b)]. Figure 8(e) compares the waveforms of the experimental Green's functions and the synthetic solution for a single receiver. Their match demonstrates that the desired Green's functions related to the steel scatterer without the scattering effects caused by the boundary of the waveguide can be recovered by the MDD method. Note that the inverted Green's functions shown in Figs. 7 and 8 only correspond to one single dipolar source located on the recording surface S^rec [e.g., Fig. 7(g)]. The complete MDD result (not shown here) corresponds to an entire set of Green's functions that completely sample the scattering related only to the steel block placed inside S^rec with full aperture.

Synthetic experiments were carried out using time-domain finite-element modeling with the Acoustics Module of comsol Multiphysics^® (Idesman and Pham, 2014). The model used is shown in Fig. 2(a), and the parameter values are given in Table I for the same configuration as used in the laboratory MDD experiments. In the synthetic MDD experiments, the physical scatterer (i.e., the steel block) used in the laboratory experiments is simulated with a perfectly rigid surface, and we deploy monopolar pressure sources at the exterior rigid boundary of the simulated 2D waveguide with a Ricker wavelet (f_p = 2 kHz) as the source signature for all the sources. The MDD scheme, wavefield separation, and other data processing steps such as bandpass filtering follow those in the laboratory experiments except that the number of iterations in the MDD solver is increased to 200.

Figure 9 shows the data, separated in-going wavefields, and MDD results for the synthetic waveguide experiments. The simulated data shown in Figs. 9(b) and 9(e) do not contain the systematic noise present in the laboratory experiments, and hence the MDD solver can run for a large number of iterations (i.e., 200) without risk of over-fitting the noise in the data. Figure 10 shows a zoom-in of the synthetic MDD results. The reference solution (i.e., simulated Green's functions), as shown in Fig. 9(l) and also in Fig. 8(d), contains a vertical-line-like modeling artefact that is present at early times before the arrival related to the interior steel scatterer. This modeling artefact is caused by the way that the simulation software comsol attempts to handle co-located sources and receivers. However, the modeling artefact is only significant at a single trace in each inverted Green's function and does not influence our results. Although the difference is small, Fig. 10(e) shows that the MDD-obtained waveform does not exactly match the reference solution [also in Fig. 8(e)]. Solving MDD equations is inherently an ill-posed inverse problem with rank deficiency (Ravasi et al., 2015; Wapenaar et al., 2011), and the inverted waveforms will be slightly distorted due to inversion errors.

The data obtained in the laboratory show a reduction in amplitude with time [e.g., Fig. 7(b)], which is not present in the synthetic experiments [e.g., Fig. 9(b)]. The energy loss of the wavefields propagating in the laboratory 2D waveguide is primarily caused by wave energy leaking through the imperfect rigid boundaries of the waveguide.

The MDD method can effectively transform recorded data contaminated by the imprint of the domain boundary to the Green's functions that are only related to the interior scatterer(s). Thus, the scattering properties of the interior medium can be studied without interference caused by the boundary of the experimental setup (e.g., a 2D waveguide). One critical requirement of this MDD application is that the sources that generate wave energy for experiments should be located outside the recording surface S^rec. The sources can be of any characteristics, and their distribution can be sparse; namely, they do not need to be distributed according to the spatial Nyquist sampling criterion. However, the source distribution needs to illuminate the medium inside S^rec from all directions. The receivers placed around S^rec are densely distributed such that the spatial Nyquist criterion is satisfied. Hence in MDD experiments, the number of sources used can be smaller than the number of receivers on S^rec. When such a number of required channels for sources and receivers cannot be satisfied, an alternative is to exploit acoustic source-receiver reciprocity to acquire the data. A movable source positioned on a recording surface in the waveguide together with a small number of receivers replacing the sources can be used. Theoretically, if the receiver is also movable, a single source-receiver pair is sufficient for recording all required data to carry out MDD. Note that exploiting source-receiver reciprocity in the laboratory requires the use of physical monopolar sources, since the radiation patterns of sources and receivers will be exchanged and the MDD method requires omnidirectional receivers.

The MDD method used to remove boundary reflections can be compared to the recently developed immersive wave experimentation application where active sources are deployed around the boundary of the experimental domain such that the emitted wavefield cancels the outward propagating waves (Becker et al., 2018, 2020). The emitted wavefield is determined by means of extrapolating the waves measured at a recording surface S^rec placed inside the experimental domain to the active domain boundary S^emt. This extrapolation has to happen in real time during a physical experiment. Hence, immersive experimentation, as an online, real-time method, is highly complicated with a high-performance control and computing system required for the wavefield extrapolation. MDD, on the other hand, is significantly simpler to operate, since it is carried out after the physical experiment. However, the MDD method and immersive wave experiments are different in several aspects. In particular, immersive experimentation allows for nonlinear, time-variant media placed inside S^rec, while the MDD method only works with linear, time-invariant media due to the linearity for the definition of the (inverted) Green's functions and the repeated measurements for each source. In addition to canceling outward propagating waves, immersive experimentation enables replacing boundary reflections with interactions from an arbitrary numerical environment, thus immersing the physical experiment in the numerical environment in real time.

A MDD method was proposed and applied to 2D acoustic wave propagation experiments such that the scattering imprint related to the experimental domain boundary is completely removed. The MDD approach enables the full-aperture reconstruction of the wavefields related to scatterers placed inside a laboratory in the absence of unwanted boundary reflections that usually contaminate such datasets. In our experiments, the laboratory data were recorded along a circular recording surface. A wavefield decomposition method for in/out separation on a general curved surface was formulated and used prior to applying the MDD method. The quality of the wavefield separation was demonstrated using a (forward) MDC relation.

The MDD scheme used in this paper involves in-going pressure data multidimensionally deconvolved from the total pressure data. For the laboratory experiment including a rigid scatterer, the MDD result contains pole artefacts and direct arrivals, both of which were successfully removed by subtracting the MDD result for a homogeneous experimental domain. By considering the limit where the evaluation point in the MDC integral approaches the recording surface, the scattering Green's functions for any pair of source and receiver on the recording surface can be obtained. Thus, the proposed MDD method allows for both removal of any exterior scattering imprints and fully sampling the scattering wavefield related to the enclosed volume.

We would like to thank Professor Kees Wapenaar and an anonymous reviewer for their positive and constructive comments. We thank the Mondaic team for their support in using the salvus software package. We thank Professor Lasse Amundsen and Professor Andrew Curtis for the discussion about MDD work in exploration geophysics. We thank Dr. Marc Serra-Garcia for the discussion about comsol. We would like to thank Thomas Haag and Christoph Bärlocher for setting up the acquisition system in the laboratory. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 694407).

The derivation of the two-way representation theorem as given in Eq. (1) starts from the acoustic wave equations. In an acoustic medium, propagating waves can be described through pressures $p(x,t)$ and particle velocities $v(x,t)$ governed by the two following first-order partial differential equations (Fokkema and van den Berg, 2013):

and

where $ρ(x)$ and $c(x)$ are medium density and wave velocity, respectively. The source term $q(x,t)$ on the right-hand side of Eq. (A1) corresponds to a distribution of volume injection sources (or monopoles), while the source term $f(x,t)$ on the right-hand side of Eq. (A2) corresponds to a distribution of body force sources (or dipoles). The spatial variable x represents the Cartesian coordinate.

The acoustic wave equations expressed as Eqs. (A1) and (A2) in the space-time $(x,t)$ domain can be transformed into the space-frequency $(x,ω)$ domain (Wapenaar and Fokkema, 2006):

and

where the symbol $̂$ in $p̂(x)$ and $v̂(x)$ denotes the equivalent counterparts in the frequency domain, and

is the compressibility. The subscript i for $v̂i(x)$ and $∂i$ denotes the scalar value and its spatial gradient in the x_i direction. Einstein's summation convention applies to repeated subscripts in this appendix, and the derivation of the following equations is done in the space-frequency domain for the convenience of expressing convolution in the time domain as multiplication in the frequency domain.

Wavefield propagation can be defined in two independent states, A and B, with their corresponding configurations shown in Fig. 11. The two states A and B are linked through the interaction quantity

composed of cross-convolution terms. One can substitute the acoustic quantities appearing in the equations for states A and B, given in Fig. 11, into this interaction quantity. The result is further integrated over the volume V enclosed by the boundary S. The volume integral is converted to a surface integral using the divergence theorem of Gauss (Arfken, 1985), which gives

Equation (A5) is called Rayleigh's reciprocity theorem of the convolution type (de Hoop, 1988, 1995; Fokkema and van den Berg, 2013). Table II gives the definitions of source fields and medium properties of states A and B. We identify state A as a physical experiment, including a physical source, and state B as a desired experiment with an impulsive, point volume injection source. The medium properties inside the integration volume V are the same in states A and B, and hence the second volume integral on the right-hand side of Eq. (A5) will vanish. Since an impulsive source exists in state B, the pressure wavefield in state B becomes Green's function $Gp,q(x,xb,ω)$ with the pressure recording at x and an impulsive point source of volume injection rate density at $xb$ (superscript q for monopole), following the notations in de Hoop (1995). The particle velocity wavefield in state B becomes the Green's function $Giv,q(x,xb,ω)$ (subscript i for recording at x_i direction). Applying the definitions in Table II to Eq. (A5), we obtain

where the quantities associated with state A are measured in a physical experiment, as shown in Fig. 1(a), with a physical source located at $xs=xa$ and a receiver inside the surface S at $xir=xb$⁠. Note that the derivation of Eq. (A6) does not require the definition of the source type in state A (physical experiment) as long as the source location $xa$ is outside the volume V. The Green's functions $Ĝiv,q$ and $Ĝp,q$ in Eq. (A6) are associated with state B as a desired experiment. Note that in state B, the medium property outside the volume V has not been specified, and the property can and will be different from that in state A, e.g., without and with an exterior rigid boundary (Wapenaar et al., 2011).