In physical acoustic laboratories, wave propagation experiments often suffer from unwanted reflections at the boundaries of the experimental setup. We propose using multidimensional deconvolution (MDD) to post-process recorded experimental data such that the scattering imprint related to the domain boundary is completely removed and only the Green's functions associated with a scattering object of interest are obtained. The application of the MDD method requires in/out wavefield separation of data recorded along a closed surface surrounding the object of interest, and we propose a decomposition method to separate such data for arbitrary curved surfaces. The MDD results consist of the Green's functions between any pair of points on the closed recording surface, fully sampling the scattered field. We apply the MDD algorithm to post-process laboratory data acquired in a two-dimensional acoustic waveguide to characterize the wavefield scattering related to a rigid steel block while removing the scattering imprint of the domain boundary. The experimental results are validated with synthetic simulations, corroborating that MDD is an effective and general method to obtain the experimentally desired Green's functions for arbitrary inhomogeneous scatterers.

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,1 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):

p̂(xir,ω)=S(p̂(xr,ω)Ĝiv,q(xr,xir,ω)v̂i(xr,ω)Ĝp,q(xr,xir,ω))nidS(xr),
(1)

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 ni) in the physical experiment, while v̂i(xr) is particle velocity recorded at S(xr) in the xi 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.

FIG. 1.

(Color online) Geometries for physical and desired experiments. All the denoted wave quantities are in the time domain (without ̂). (a) A physical experiment with a source (red star at xs) generating wave energy. The blue dots (xr) represent receivers around the surface S (dashed blue line) with normal n. The blue triangle at xir denotes a receiver inside the experimental volume V. The solid black circle denotes the rigid boundary surrounding the experimental volume, while the black arrows denote wave propagation between the source and receivers. (b) The desired experiment associated with the Green's functions in Eq. (1) of which source and receiver positions are at xir and xr, respectively. The short notation Gin represents the in-going Green's functions Giv,q,in and Gp,q,in that are zero in this desired medium (i.e., homogeneous outside V). The dashed black circle represents a transparent boundary compared to the rigid boundary in the physical experiment. (c) The wavefields recorded around the receiver surface S are decomposed into in-going and out-going components, denoted by the black arrows.

FIG. 1.

(Color online) Geometries for physical and desired experiments. All the denoted wave quantities are in the time domain (without ̂). (a) A physical experiment with a source (red star at xs) generating wave energy. The blue dots (xr) represent receivers around the surface S (dashed blue line) with normal n. The blue triangle at xir denotes a receiver inside the experimental volume V. The solid black circle denotes the rigid boundary surrounding the experimental volume, while the black arrows denote wave propagation between the source and receivers. (b) The desired experiment associated with the Green's functions in Eq. (1) of which source and receiver positions are at xir and xr, respectively. The short notation Gin represents the in-going Green's functions Giv,q,in and Gp,q,in that are zero in this desired medium (i.e., homogeneous outside V). The dashed black circle represents a transparent boundary compared to the rigid boundary in the physical experiment. (c) The wavefields recorded around the receiver surface S are decomposed into in-going and out-going components, denoted by the black arrows.

Close modal

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

p̂(xir)=S[(p̂in(xr)+p̂out(xr))Ĝiv,q(v̂iin(xr)+v̂iout(xr))Ĝp,q]nidS(xr).
(2)

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:

p̂(xir)=S(p̂in(xr)Ĝiv,qv̂iin(xr)Ĝp,q)nidS(xr).
(3)

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:

p(xir,t)=2Spin(xr,t)*Giv,q(xr,t|xir,0)nidS(xr),
(4)

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,fni 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

p(xir,t)=2Spin(xr,t)*G¯d(xir,t|xr,0)dS(xr),
(5)

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 ptotalpin 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.

1. MDD for a laboratory configuration: Full-aperture sampling of the scattering 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 xir 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 Srec (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 Srec, 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 Srec (“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

G¯dscat=G¯dhetG¯dhom,
(6)

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 Srec 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 Srec.

FIG. 2.

(Color online) (a) Schematic of the acquisition geometry used for MDD experiments. The emitting surface Semt is composed of evenly spaced sources (red stars) located at the rigid boundary (solid black circle) surrounding a 2D waveguide. Two circular layers of densely spaced receivers (blue dots) make up the recording surface Srec. (b) 2D waveguide used for laboratory experiments. The scatterer under study is a rigid steel block placed inside Srec.

FIG. 2.

(Color online) (a) Schematic of the acquisition geometry used for MDD experiments. The emitting surface Semt is composed of evenly spaced sources (red stars) located at the rigid boundary (solid black circle) surrounding a 2D waveguide. Two circular layers of densely spaced receivers (blue dots) make up the recording surface Srec. (b) 2D waveguide used for laboratory experiments. The scatterer under study is a rigid steel block placed inside Srec.

Close modal

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 Srec 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.

FIG. 3.

(Color online) (a) Selected group of receivers (red dots) in the inner layer of the recording surface Sr1 (blue dots). A local coordinate x,z is constructed to include these receivers, and the green dot denotes the central receiver. (b) Zoom-in of the local coordinate in the graph (a) with a different orientation. The blue circles denote receivers in the outer layer of the recording surface Sr2.

FIG. 3.

(Color online) (a) Selected group of receivers (red dots) in the inner layer of the recording surface Sr1 (blue dots). A local coordinate x,z is constructed to include these receivers, and the green dot denotes the central receiver. (b) Zoom-in of the local coordinate in the graph (a) with a different orientation. The blue circles denote receivers in the outer layer of the recording surface Sr2.

Close modal

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 kn is the regularly sampled horizontal wavenumber (still with n1,,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,

ptot(f)=[p̂(x1r1,z1r1,f),,p̂(xNr1,zNr1,f),p̂(x1r2,z1r2,f),,p̂(xNr2,zNr2,f)]T,
(7a)
ssep(f)=[s̃in(f,k1),,s̃in(f,kN),s̃out(f,k1),,s̃out(f,kN)]T.
(7b)

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

L(f)=[Lr1out(f)Lr1in(f)Lr2out(f)Lr2in(f)]
(8)

with the submatrices

Lr1out(f)=[ei2π(kmxnr1kz(km)|znr1|)]n,m=1,,N,
(9a)
Lr1in(f)=[ei2π(kmxnr1+kz(km)|znr1|)]n,m=1,,N,
(9b)
Lr2out(f)=[ei2π(kmxnr2kz(km)|znr2|)]n,m=1,,N,
(9c)
Lr2in(f)=[ei2π(kmxnr2+kz(km)|znr2|)]n,m=1,,N,
(9d)

and kz(km)=(f/c)2km2 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,

ptot=Lssep,
(10)

and its block form for each frequency

[pr1(f)pr2(f)]=[Lr1out(f)Lr1in(f)Lr2out(f)Lr2in(f)][sout(f)sin(f)].
(11)

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),

ssep(LHL+ϵI)1LHptot,
(12)

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):

pr1in(f)=Lr1insin.
(13)

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 Srec contain small artefacts (i.e., presence of unphysical events) toward the beginning and end of time traces, which are suppressed using tapering.

FIG. 4.

Flow chart describing the wavefield decomposition method.

FIG. 4.

Flow chart describing the wavefield decomposition method.

Close modal

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.

FIG. 5.

(Color online) Acquisition geometry of the synthetic experiment validating the MDC relation. A pressure source (red star) is located at the rigid boundary (solid black circle). Two circular layers of receivers (blue dots) make up the recording surface Srec, and an extra circular array of receivers (magenta dots for Sirec) is placed inside Srec.

FIG. 5.

(Color online) Acquisition geometry of the synthetic experiment validating the MDC relation. A pressure source (red star) is located at the rigid boundary (solid black circle). Two circular layers of receivers (blue dots) make up the recording surface Srec, and an extra circular array of receivers (magenta dots for Sirec) is placed inside Srec.

Close modal
TABLE I.

Specifications for the synthetic and physical experiments.

ParameterDefinitionValue
c0 Acoustic velocity (air) 347 m/s 
ρ0 Air density 1 kg/m3 
Remt Radius of 2D waveguide and emitting surface Semt 0.6621 m 
Rrec1 Radius of recording surface (inner layer) Sr1 0.3529 m 
Rrec2 Radius of recording surface (outer layer) Sr2 0.3729 m 
Rirec Radius of interior receiver array Sirec 0.15 m 
Nemt Number of sources 52 
Nrec Number of receivers at the recording surface (two layers) Srec 228 
Nirec Number of receivers at the interior receiver array Sirec 114 
fp Peak frequency of Ricker wavelet 2 kHz 
tmax Time length of physical/synthetic experiments 37.5 ms 
dt Time step 6.375 ×105s 
ParameterDefinitionValue
c0 Acoustic velocity (air) 347 m/s 
ρ0 Air density 1 kg/m3 
Remt Radius of 2D waveguide and emitting surface Semt 0.6621 m 
Rrec1 Radius of recording surface (inner layer) Sr1 0.3529 m 
Rrec2 Radius of recording surface (outer layer) Sr2 0.3729 m 
Rirec Radius of interior receiver array Sirec 0.15 m 
Nemt Number of sources 52 
Nrec Number of receivers at the recording surface (two layers) Srec 228 
Nirec Number of receivers at the interior receiver array Sirec 114 
fp Peak frequency of Ricker wavelet 2 kHz 
tmax Time length of physical/synthetic experiments 37.5 ms 
dt Time step 6.375 ×105s 

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 fp = 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 (fp = 2000 Hz) used in the simulations [Fig. 6(a)]. Hence, the MDC equation can be recast as

q(t)*p(xir,t)=2Spin(xr,t)*(G¯d(xir,t,xr,0)*q(t))dS(xr).
(14)
FIG. 6.

(Color online) (a) Dipolar Green's function convolved with a Ricker wavelet q(t). The upper panel shows the model used to simulate the Green's function with an absorbing boundary (dashed black circle). The cyan dot denotes the source location at the recording surface Srec, while receivers are at the interior recording surface Sirec (circular green array). (b) The separated in-going pressure wavefield. The red star denotes the source location at the rigid boundary (solid black circle), and the receivers are on Srec (blue circle). (c) The computed right-hand side of Eq. (14). (d) The computed left-hand side of Eq. (14).

FIG. 6.

(Color online) (a) Dipolar Green's function convolved with a Ricker wavelet q(t). The upper panel shows the model used to simulate the Green's function with an absorbing boundary (dashed black circle). The cyan dot denotes the source location at the recording surface Srec, while receivers are at the interior recording surface Sirec (circular green array). (b) The separated in-going pressure wavefield. The red star denotes the source location at the rigid boundary (solid black circle), and the receivers are on Srec (blue circle). (c) The computed right-hand side of Eq. (14). (d) The computed left-hand side of Eq. (14).

Close modal

The wavefield separation algorithm proposed in Sec. II B was applied to the simulated wavefields ptot(xr,t) recorded at the recording surface Srec 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.5c0da1 (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 Srec 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 fp = 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.

FIG. 7.

(Color online) Laboratory MDD results. (a) A homogeneous experimental domain with a source (red star) located at the rigid boundary (solid black circle) and receivers on Srec (blue circle) with azimuth ϕ. (b) and (c) The data recorded for the scenario (a) and their separated in-going components. (d)–(f) A heterogeneous experimental domain with the steel block (black rectangle), the recorded data, and separated in-going components. (g) and (h) A desired waveguide experiment including the scatterer with a radiation boundary condition (dashed black line) and the corresponding MDD result. The cyan star denotes a dipolar source. (j) and (k) A desired homogeneous waveguide experiment and the MDD result. (i) The scattering Green's functions obtained by subtracting the MDD result in the graph (k) from that in the graph (e). (l) The synthetic reference solution.

FIG. 7.

(Color online) Laboratory MDD results. (a) A homogeneous experimental domain with a source (red star) located at the rigid boundary (solid black circle) and receivers on Srec (blue circle) with azimuth ϕ. (b) and (c) The data recorded for the scenario (a) and their separated in-going components. (d)–(f) A heterogeneous experimental domain with the steel block (black rectangle), the recorded data, and separated in-going components. (g) and (h) A desired waveguide experiment including the scatterer with a radiation boundary condition (dashed black line) and the corresponding MDD result. The cyan star denotes a dipolar source. (j) and (k) A desired homogeneous waveguide experiment and the MDD result. (i) The scattering Green's functions obtained by subtracting the MDD result in the graph (k) from that in the graph (e). (l) The synthetic reference solution.

Close modal

The data recorded on the surface Srec 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., ptotalpin 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||Axb||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=106, 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 fp = 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 Srec. 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 Srec [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 Srec with full aperture.

FIG. 8.

(Color online) Zoom-in of laboratory MDD results. (e) The comparison between the picked traces (dashed blue and red lines) in the graphs (b) and (d). Otherwise, key as in Fig. 7.

FIG. 8.

(Color online) Zoom-in of laboratory MDD results. (e) The comparison between the picked traces (dashed blue and red lines) in the graphs (b) and (d). Otherwise, key as in Fig. 7.

Close modal

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 (fp = 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.

FIG. 9.

(Color online) Synthetic MDD results. Key as in Fig. 7.

FIG. 9.

(Color online) Synthetic MDD results. Key as in Fig. 7.

Close modal
FIG. 10.

(Color online) Zoom-in of synthetic MDD results. Key as in Fig. 9.

FIG. 10.

(Color online) Zoom-in of synthetic MDD results. Key as in Fig. 9.

Close modal

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 Srec. 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 Srec from all directions. The receivers placed around Srec 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 Srec. 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 Srec placed inside the experimental domain to the active domain boundary Semt. 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 Srec, 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):

p(x,t)t+ρ(x)c2(x)·v(x,t)=ρ(x)c2(x)q(x,t)
(A1)

and

ρ(x)v(x,t)t+p(x,t)=f(x,t),
(A2)

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):

jωκ(x,ω)p̂(x,ω)+iv̂i(x,ω)=q̂(x,ω)
(A3)

and

jωρ(x)v̂i(x,ω)+ip̂(x,ω)=f̂i(x,ω),
(A4)

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

κ(x)=ρ1(x)c2(x)

is the compressibility. The subscript i for v̂i(x) and i denotes the scalar value and its spatial gradient in the xi 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

i(p̂Av̂iBv̂iAp̂B)

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

S(p̂Av̂iBv̂iAp̂B)nidS(xr)=V[p̂Aq̂Bv̂iAf̂iBq̂Ap̂B+f̂iAv̂iB]dV(x)+jωV[(κAκB)p̂Ap̂B(ρAρB)v̂iAv̂iB]dV(x).
(A5)
FIG. 11.

(Color online) Configurations of the acoustic states A and B. The wave equations for states A and B are expressed in a short form with the spatial and angular frequency variables omitted compared to Eqs. (A3) and (A4). State-A: The red dot denotes a source located at xa outside the volume V, which is enclosed by the surface S with outward-pointing normal vector n. The blue dot xr is at the surface S. State-B: The source (red dot) is located at xb inside V.

FIG. 11.

(Color online) Configurations of the acoustic states A and B. The wave equations for states A and B are expressed in a short form with the spatial and angular frequency variables omitted compared to Eqs. (A3) and (A4). State-A: The red dot denotes a source located at xa outside the volume V, which is enclosed by the surface S with outward-pointing normal vector n. The blue dot xr is at the surface S. State-B: The source (red dot) is located at xb inside V.

Close modal

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 xi direction). Applying the definitions in Table II to Eq. (A5), we obtain

p̂A(xb,ω)=S(p̂A(xr,ω)Ĝiv,q(xr,xb,ω)v̂iA(xr,ω)Ĝp,q(xr,xb,ω))nidS(xr),
(A6)

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).

TABLE II.

Definitions of states A and B.

State A: Physical experimentState B: Desired experiment
Wavefields p̂A(x,ω),viA(x,ω) p̂B=Gp,q(x,xb,ω),viB=Giv,q(x,xb,ω) 
Sources q̂A(x,ω)=δ(xxa)â(ω) q̂B(x,ω)=δ(xxb) 
 f̂iA(x,ω)=0 f̂iB(x,ω)=0 
Medium κA(x)=κB(x),ρA(x)=ρB(x) for xV 
State A: Physical experimentState B: Desired experiment
Wavefields p̂A(x,ω),viA(x,ω) p̂B=Gp,q(x,xb,ω),viB=Giv,q(x,xb,ω) 
Sources q̂A(x,ω)=δ(xxa)â(ω) q̂B(x,ω)=δ(xxb) 
 f̂iA(x,ω)=0 f̂iB(x,ω)=0 
Medium κA(x)=κB(x),ρA(x)=ρB(x) for xV 
1

In other words, the traction free surface on land or the water-air interface (pressure-release surface) in a marine environment.

1.
Afanasiev
,
M.
,
Boehm
,
C.
,
van Driel
,
M.
,
Krischer
,
L.
,
Rietmann
,
M.
,
May
,
D. A.
,
Knepley
,
M. G.
, and
Fichtner
,
A.
(
2019
). “
Modular and flexible spectral-element waveform modelling in two and three dimensions
,”
Geophys. J. Int.
216
(
3
),
1675
1692
.
2.
Amundsen
,
L.
(
1993
). “
Wavenumber-based filtering of marine point-source data
,”
Geophysics
58
(
9
),
1335
1348
.
3.
Amundsen
,
L.
(
2001
). “
Elimination of free-surface related multiples without need of the source wavelet
,”
Geophysics
66
(
1
),
327
341
.
4.
Amundsen
,
L.
,
Ikelle
,
L. T.
, and
Berg
,
L. E.
(
2001
). “
Multidimensional signature deconvolution and free-surface multiple elimination of marine multicomponent ocean-bottom seismic data
,”
Geophysics
66
(
5
),
1594
1604
.
5.
Amundsen
,
L.
, and
Robertsson
,
J. O. A.
(
2014
). “
Wave equation processing using finite-difference propagators, part 1: Wavefield dissection and imaging of marine multicomponent seismic data
,”
Geophysics
79
(
6
),
T287
T300
.
6.
Arfken
,
G.
(
1985
).
Mathematical Methods for Physicists
, 3rd ed. (
Academic Press
,
Orlando, FL
).
7.
Aster
,
R.
,
Borchers
,
B.
, and
Thurber
,
C.
(
2013
).
Parameter Estimation and Inverse Problems
(
Elsevier Science
,
Amsterdam
).
8.
Becker
,
T. S.
,
Börsing
,
N.
,
Haag
,
T.
,
Bärlocher
,
C.
,
Donahue
,
C. M.
,
Curtis
,
A.
,
Robertsson
,
J. O. A.
, and
van Manen
,
D.-J.
(
2020
). “
Real-time immersion of physical experiments in virtual wave-physics domains
,”
Phys. Rev. Applied
13
,
064061
.
9.
Becker
,
T. S.
,
van Manen
,
D.-J.
,
Donahue
,
C. M.
,
Bärlocher
,
C.
,
Börsing
,
N.
,
Broggini
,
F.
,
Haag
,
T.
,
Robertsson
,
J. O. A.
,
Schmidt
,
D. R.
,
Greenhalgh
,
S. A.
, and
Blum
,
T. E.
(
2018
). “
Immersive wave propagation experimentation: Physical implementation and one-dimensional acoustic results
,”
Phys. Rev. X
8
(
3
),
031011
.
10.
Bellezza
,
C.
, and
Poletto
,
F.
(
2014
). “
Multidimensional deconvolution and processing of seismic-interferometry Arctic data
,”
Geophysics
79
(
3
),
WA25
WA38
.
11.
Beranek
,
L. L.
, and
Sleeper
,
H. P.
(
1946
). “
The design and construction of anechoic sound chambers
,”
J. Acoust. Soc. Am.
18
(
1
),
140
150
.
12.
Beyene
,
S.
, and
Burdisso
,
R. A.
(
1997
). “
A new hybrid passive-active noise absorption system
,”
J. Acoust. Soc. Am.
101
(
3
),
1512
1515
.
13.
Blum
,
T. E.
,
van Wijk
,
K.
,
Snieder
,
R.
, and
Willis
,
M. E.
(
2011
). “
Laser excitation of a fracture source for elastic waves
,”
Phys. Rev. Lett.
107
,
275501
.
14.
Cassereau
,
D.
, and
Fink
,
M.
(
1992
). “
Time-reversal of ultrasonic fields. III. theory of the closed time-reversal cavity
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Cont.
39
(
5
),
579
592
.
15.
Day
,
A.
,
Klüver
,
T.
,
Slöner
,
W.
,
Tabti
,
H.
, and
Carlson
,
D.
(
2013
). “
Wavefield-separation methods for dual-sensor towed-streamer data
,”
Geophysics
78
(
2
),
WA55
WA70
.
16.
de Hoop
,
A. T.
(
1988
). “
Time-domain reciprocity theorems for acoustic wave fields in fluids with relaxation
,”
J. Acoust. Soc. Am.
84
(
5
),
1877
1882
.
17.
de Hoop
,
A.
(
1995
).
Handbook of Radiation and Scattering of Waves: Acoustic Waves in Fluids, Elastic Waves in Solids, Electromagnetic Waves
(
Academic
,
San Diego
).
18.
Duijndam
,
A. J. W.
,
Schonewille
,
M. A.
, and
Hindriks
,
C. O. H.
(
1999
). “
Reconstruction of band-limited signals, irregularly sampled along one spatial direction
,”
Geophysics
64
(
2
),
524
538
.
19.
Ferber
,
R.
, and
van Manen
,
D.-J.
(
2017
). “
Noise transfer in variable-depth streamer deghosting
,”
Geophys. Prospect.
65
(
4
),
903
912
.
20.
Fink
,
M.
(
1992
). “
Time reversal of ultrasonic fields. I. Basic principles
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Cont.
39
(
5
),
555
566
.
21.
Fokkema
,
J.
, and
van den Berg
,
P.
(
2013
).
Seismic Applications of Acoustic Reciprocity
(
Elsevier Science
,
Amsterdam
).
22.
Guicking
,
D.
, and
Lorenz
,
E.
(
1984
). “
An active sound absorber with porous plate
,”
J. Vib. Acoust. Stress. Reliab. Des.
106
(
3
),
389
392
.
23.
Habault
,
D.
,
Friot
,
E.
,
Herzog
,
P.
, and
Pinhede
,
C.
(
2017
). “
Active control in an anechoic room: Theory and first simulations
,”
Acta Acust. United Acust.
103
(
3
),
369
378
.
24.
Holvik
,
E.
, and
Amundsen
,
L.
(
2005
). “
Elimination of the overburden response from multicomponent source and receiver seismic data, with source designature and decomposition into PP-, PS-, SP-, and SS-wave responses
,”
Geophysics
70
(
2
),
S43
S59
.
25.
Hulsebos
,
E. M.
(
2004
). “
Auralization using wave field synthesis
,” Ph.D. thesis,
Technische Universiteit Delft
,
Delft, Netherlands
.
26.
Hunziker
,
J. W.
,
Fan
,
Y.
,
Slob
,
E. C.
,
Wapenaar
,
K.
, and
Snieder
,
R.
(
2010
). “
Solving spatial sampling problems in 2D-CSEM interferometry using elongated sources
,”
Proceedings of the 72nd EAGE Conference and Exhibition Incorporating SPE EUROPEC 2010
, June 14–17, Barcelona, Spain.
27.
Idesman
,
A.
, and
Pham
,
D.
(
2014
). “
Accurate finite element modeling of acoustic waves
,”
Comput. Phys. Commun.
185
(
7
),
2034
2045
.
28.
Kosloff
,
R.
, and
Kosloff
,
D.
(
1986
). “
Absorbing boundaries for wave propagation problems
,”
J. Comput. Phys.
63
(
2
),
363
376
.
29.
Larose
,
E.
,
Planes
,
T.
,
Rossetto
,
V.
, and
Margerin
,
L.
(
2010
). “
Locating a small change in a multiple scattering environment
,”
Appl. Phys. Lett.
96
(
20
),
204101
.
30.
McDonald
,
J.
,
Gardner
,
G.
, and
Hilterman
,
F.
(
1983
).
Seismic Studies in Physical Modeling
(
Springer
,
Amsterdam)
.
31.
Minato
,
S.
,
Matsuoka
,
T.
,
Tsuji
,
T.
,
Draganov
,
D.
,
Hunziker
,
J.
, and
Wapenaar
,
K.
(
2011
). “
Seismic interferometry using multidimensional deconvolution and crosscorrelation for crosswell seismic reflection data without borehole sources
,”
Geophysics
76
(
1
),
SA19
SA34
.
32.
Mo
,
Y.
,
Greenhalgh
,
S. A.
,
Robertsson
,
J. O.
, and
Karaman
,
H.
(
2015
). “
The development and testing of a 2D laboratory seismic modelling system for heterogeneous structure investigations
,”
J. Appl. Geophys.
116
,
224
235
.
33.
Munjal
,
M. L.
(
2002
).
IUTAM Symposium on Designing for Quietness
(
Kluwer Academic Publishers
,
Dordrecht, Netherlands
).
34.
Paige
,
C. C.
, and
Saunders
,
M. A.
(
1982
). “
LSQR: An algorithm for sparse linear equations and sparse least squares
,”
ACM Trans. Math. Softw.
8
(
1
),
43
71
.
35.
Poletto
,
F.
,
Bellezza
,
C.
, and
Farina
,
B.
(
2014
). “
Virtual reflector multidimensional deconvolution: Inversion issues for convolutive-type interferometry
,”
Geophys. J. Int.
196
(
2
),
1018
1024
.
36.
Rakotonarivo
,
S. T.
,
Kuperman
,
W. A.
, and
Williams
,
E. G.
(
2013
). “
Prediction of a body's structural impedance and scattering properties using correlation of random noise
,”
J. Acoust. Soc. Am.
134
(
6
),
4401
4411
.
37.
Ramírez
,
A. C.
, and
Weglein
,
A. B.
(
2009
). “
Green's theorem as a comprehensive framework for data reconstruction, regularization, wavefield separation, seismic interferometry, and wavelet estimation: A tutorial
,”
Geophysics
74
(
6
),
W35
W62
.
38.
Ravasi
,
M.
,
Meles
,
G.
,
Curtis
,
A.
,
Rawlinson
,
Z.
, and
Yikuo
,
L.
(
2015
). “
Seismic interferometry by multidimensional deconvolution without wavefield separation
,”
Geophys. J. Int.
202
(
1
),
1
16
.
39.
Ravasi
,
M.
, and
Vasconcelos
,
I.
(
2020
). “
PyLops—A linear-operator Python library for scalable algebra and optimization
,”
SoftwareX
11
,
100361
.
40.
Redwood
,
M.
(
1960
).
Mechanical Waveguides: The Propagation of Acoustic and Ultrasonic Waves in Fluids and Solids with Boundaries
(
Pergamon
,
New York
).
41.
Riyanti
,
C. D.
,
van Borselen
,
R. G.
,
van den Berg
,
P. M.
, and
Fokkema
,
J. T.
(
2008
). “
Pressure wave-field deghosting for non-horizontal streamers
,” in
SEG Technical Program Expanded Abstracts 2008
(
Society of Exploration Geophysicists
,
Tulsa, OK
), pp.
2652
2656
.
42.
Robertsson
,
J. O.
,
van Manen
,
D.-J.
,
Schmelzbach
,
C.
,
Van Renterghem
,
C.
, and
Amundsen
,
L.
(
2015
). “
Finite-difference modelling of wavefield constituents
,”
Geophys. J. Int.
203
(
2
),
1334
1342
.
43.
Ruigrok
,
E.
,
Campman
,
X.
,
Draganov
,
D.
, and
Wapenaar
,
K.
(
2010
). “
High-resolution lithospheric imaging with seismic interferometry
,”
Geophys. J. Int.
183
(
1
),
339
357
.
44.
Shao
,
C.
,
Long
,
H.
,
Cheng
,
Y.
, and
Liu
,
X.
(
2019
). “
Low-frequency perfect sound absorption achieved by a modulus-near-zero metamaterial
,”
Sci. Rep.
9
(
1
),
13482
.
45.
Smith
,
J. P.
,
Johnson
,
B. D.
, and
Burdisso
,
R. A.
(
1999
). “
A broadband passive-active sound absorption system
,”
J. Acoust. Soc. Am.
106
(
5
),
2646
2652
.
46.
Sternini
,
S.
,
Sarkar
,
J.
,
Rakotonarivo
,
S.
,
Bottero
,
A.
,
Williams
,
E. G.
,
Tippmann
,
J. D.
, and
Kuperman
,
W. A.
(
2019
). “
Determining structural acoustic properties from noise-based holographic measurements
,”
J. Acoust. Soc. Am.
146
(
4
),
2943
2943
.
47.
Thomsen
,
H. R.
,
van Manen
,
D.-J.
, and
Robertsson
,
J.
(
2018
). “
Exact wavefield separation on an elastic free surface with sharp corners
,” in
SEG Technical Program Expanded Abstracts 2018
(
Society of Exploration Geophysicists
,
Tulsa, OK
), pp.
5017
5021
.
48.
Thomson
,
C. J.
(
2012
). “
Research Note: Internal/external seismic source wavefield separation and cancellation
,”
Geophys. Prospect.
60
(
3
),
581
587
.
49.
Trinh
,
V. H.
,
Langlois
,
V.
,
Guilleminot
,
J.
,
Perrot
,
C.
,
Khidas
,
Y.
, and
Pitois
,
O.
(
2019
). “
Tuning membrane content of sound absorbing cellular foams: Fabrication, experimental evidence and multiscale numerical simulations
,”
Materials Design
162
,
345
361
.
50.
van Dalen
,
K. N.
,
Wapenaar
,
K.
, and
Halliday
,
D. F.
(
2014
). “
Surface wave retrieval in layered media using seismic interferometry by multidimensional deconvolution
,”
Geophys. J. Int.
196
(
1
),
230
242
.
51.
van der Neut
,
J.
, and
Herrmann
,
F. J.
(
2013
). “
Interferometric redatuming by sparse inversion
,”
Geophys. J. Int.
192
(
2
),
666
670
.
52.
Vasmel
,
M.
,
Robertsson
,
J. O. A.
,
van Manen
,
D.-J.
, and
Curtis
,
A.
(
2013
). “
Immersive experimentation in a wave propagation laboratory
,”
J. Acoust. Soc. Am.
134
(
6
),
EL492
EL498
.
53.
Wang
,
Y.
(
2015
). “
Frequencies of the ricker wavelet
,”
Geophysics
80
(
2
),
A31
A37
.
54.
Wapenaar
,
C.
, and
Berkhout
,
A.
(
1989
). “
Advances in exploration geophysics
,” in
Elastic Wave Field Extrapolation: Redatuming of Single- and Multi-component Seismic Data
(
Elsevier
,
Amsterdam
).
55.
Wapenaar
,
K.
, and
Fokkema
,
J.
(
2006
). “
Green's function representations for seismic interferometry
,”
Geophysics
71
(
4
),
SI33
SI46
.
56.
Wapenaar
,
K.
,
van der Neut
,
J.
, and
Ruigrok
,
E.
(
2008
). “
Passive seismic interferometry by multidimensional deconvolution
,”
Geophysics
73
(
6
),
A51
A56
.
57.
Wapenaar
,
K.
,
van der Neut
,
J.
,
Ruigrok
,
E.
,
Draganov
,
D.
,
Hunziker
,
J.
,
Slob
,
E.
,
Thorbecke
,
J.
, and
Snieder
,
R.
(
2011
). “
Seismic interferometry by crosscorrelation and by multidimensional deconvolution: A systematic comparison
,”
Geophys. J. Int.
185
(
3
),
1335
1364
.
58.
Weemstra
,
C.
,
Draganov
,
D.
,
Ruigrok
,
E. N.
,
Hunziker
,
J.
,
Gomez
,
M.
, and
Wapenaar
,
K.
(
2017a
). “
Application of seismic interferometry by multidimensional deconvolution to ambient seismic noise recorded in Malargüe, Argentina
,”
Geophys. J. Int.
208
(
2
),
693
714
.
59.
Weemstra
,
C.
,
Wapenaar
,
K.
, and
van Dalen
,
K. N.
(
2017b
). “
Reflecting boundary conditions for interferometry by multidimensional deconvolution
,”
J. Acoust. Soc. Am.
142
(
4
),
2242
2257
.
60.
Williams
,
E. G.
,
Tippmann
,
J. D.
,
Rakotonarivo
,
S. T.
,
Waters
,
Z. J.
,
Roux
,
P.
, and
Kuperman
,
W. A.
(
2017
). “
Experimental estimation of in vacuo structural admittance using random sources in a non-anechoic room
,”
J. Acoust. Soc. Am.
142
(
1
),
103
109
.
61.
Zhang
,
X.
,
Qu
,
Z.
, and
Wang
,
H.
(
2020
). “
Engineering acoustic metamaterials for sound absorption: From uniform to gradient structures
,”
iScience
23
(
5
),
101110
.
62.
Zhou
,
F.
,
Wang
,
B.
,
Fan
,
J.
, and
Peng
,
Z.
(
2019
). “
Theoretical and numerical studies on in vacuo structural admittance of an infinite, coated cylindrical shell
,”
Acoust. Phys.
65
(
1
),
14
22
.