Virtual acoustics in inhomogeneous media with single-sided access

A virtual acoustic source inside a medium can be created by emitting a time-reversed point-source response from the enclosing boundary into the medium. However, in many practical situations, like in non-destructive testing, near-field acoustic holography or geophysical holography, the medium can be accessed from one side only. In those cases the time-reversal approach is not exact, and it breaks down in inhomogeneous media with strong impedance contrasts. Here, we demonstrate the experimental design and use of complex focusing functions to create virtual acoustic sources and virtual receivers inside an inhomogeneous medium with single-sided access. The retrieved virtual acoustic responses between those sources and receivers mimic the complex propagation and multiple scattering paths of waves that would be ignited by physical sources and recorded by physical receivers inside the medium. The possibility to predict complex virtual acoustic responses between any two points inside an inhomogeneous medium, without needing a detailed model of the medium, has large potential for holographic imaging and monitoring of objects with single-sided access, ranging from photoacoustic medical imaging to the monitoring of induced-earthquake waves all the way from the source to the earth's surface.

2 Kees Wapenaar, Joeri Brackenhoff, Jan Thorbecke, Joost van der Neut, Evert Slob, Eric Verschuur these interfaces. The top panel shows the time-reversal of the response V (x, s, t) to a point source at s in the third layer of the medium, as a function of receiver position x = (x, z) along the boundary and time t. V stands for the normal component of the particle velocity. Only the response at the upper boundary is shown, but the response is available along the entire enclosing boundary S. The time-reversed response V (x, s, −t) is fed to sources (the red dots) at the original positions of the receivers, which emit the wave field back into the medium. The other panels in column (a) show "snapshots" (i.e., wave fields frozen at constant time) of the wave field propagating through the medium. For negative time (· · · −t 2 , −t 1 · · · ), the field follows the same paths as the original field, but in opposite direction. Then, at t = 0, the field focuses at the position s of the original source. Because there is no sink to absorb the focused field, the wave field continues its propagation, away from the focal point. Hence, the focal point acts as a virtual source. The snapshots for positive time (· · · +t 1 , +t 2 · · · ) show the response to this virtual source. The virtual source is omni-directional and radiates a perfect replica of the original field into the inhomogeneous medium. Mathematically, time-reversal acoustics is formulated as follows (Derode et al. 2003): G(r, s, t) + G(r, s, −t) = 2 S G(r, x, t) "propagator" * V (x, s, −t) "secondary sources"

dx
(1) (see Appendices). The time-reversed field V (x, s, −t) is propagated through the medium by the Green's function G(r, x, t) from the sources at x on the boundary S to any receiver position r inside the medium (the asterisk denotes convolution). The integral is taken along all sources x on the closed boundary. Note that the right-hand side resembles Huygens' principle, which states that each point of an incident wave field acts as a secondary source, except that here the secondary sources on S consist of time-reversed measurements rather than an actual incident field. On the left-hand side, the time-reversed Green's function G(r, s, −t) represents the wave field at negative time that converges to the focal point s; the Green's function G(r, s, t) is the response at positive time to the virtual source at s.
In many practical situations, the medium can be accessed from one side only. Figure 1 (b) shows what happens when the time-reversed response is emitted into the medium by sources (red dots) at the upper boundary only. The field still focuses at t = 0, but in addition several ghost foci occur at t = 0. The field at positive time is a virtual-source response, contaminated by artefacts, caused by the ghost foci. Moreover, because the focal point is illuminated mainly from above, the virtual source is far from isotropic: it radiates mainly downward.   (Marchenko 1955;Lamb 1980;Chadan & Sabatier 1989) and its applications in 1D autofocusing (Rose 2001(Rose , 2002Broggini & Snieder 2012). The upper panel in Figure 1(c) shows a 2D focusing function F (x, s, t), for the same focal point s as in the time-reversal example. Note that the main event (indicated in blue) is the same as that in V (x, s, −t) in the upper panel in Figure 1(b), but the other events in F (x, s, t) come after the main event (instead of preceding it, like in V (x, s, −t)). The snapshots in Figure 1(c) show the propagation of this focusing function through the medium.
Mathematically, the emission of the focusing function F (x, s, t) into the medium by sources at x at the upper boundary S 0 is described by (see Appendices). The right-hand side resembles again Huygens' principle, this time with the focusing function defining secondary sources on S 0 only. The left-hand side represents the virtual-source response G(r, s, t), contaminated by artefacts that are anti-symmetric in time.
Because the anti-symmetric term vanishes at t = 0, the panel at t = 0 in Figure 1(c) shows a "clean" focus. Like in the time-reversal method, the focused field acts as a virtual source.
The snapshots at positive time show that this virtual source radiates mainly downward.
Next, we symmetrize both sides of equation (2), by adding the time-reversal. This suppresses the anti-symmetric artefacts: (see Appendices). Note that the left-hand side is identical to that in equation (1). However, unlike equation (1), the right-hand side of equation (3) contains an integral along the accessible boundary S 0 only. Symmetrization implies addition of the snapshots at negative times in Figure 1(c) to those at the corresponding positive times and vice versa, see Figure 1(d). Note that these snapshots are nearly identical to those obtained by emitting the time-reversed response into the medium from the entire enclosing boundary (Figure 1(a)). The remaining artefacts are caused by the finite source aperture and the fact that evanescent waves are neglected in equations (2) and (3) (see Appendices).
The snapshots in Figure 1 (for both methods) were obtained by numerically modelling the medium's response to fields emitted from (parts of) its boundary. These snapshots nicely visualise the propagation, scattering, focusing and defocusing of the fields inside the medium.
In practical situations these fields are not visible, unless receivers would be placed through- out the medium, which is of course not feasible. However, our focusing methodology can be extended to create not only virtual sources, but also virtual receivers anywhere inside the medium. As input we need the reflection response of the medium, measured with sources and receivers at the accessible boundary S 0 only (hence, no physical sources nor receivers are needed inside the medium). In the following we apply this to ultrasonic physical model data and seismic reflection data. The size of the model is 70×600×600 mm. The model is placed in a watertank and probed with ultrasound, emitted and received by piezo-electric transducers in the water. The acquisition is carried out along a horizontal diagonal line (indicated in Figure 2(a)), 12 mm above the upper boundary of the model and perpendicular to its main structures. A 2D cross-section of the model below the acquisition line is shown in Figure 2 The recorded response is deconvolved for the sweep signal, effectively compressing the source signal to a short zero-phase pulse with a central frequency of 1.1 MHz (Blacquière & Koek 1997). This experiment is repeated 106 times, with the source at the same position and the receiver moving along the acquisition line in steps of 1.25 mm. Next, the source is moved 1.25 mm along the line and again 106 traces are recorded. This whole process is carried out 301 times, leading to a recorded reflection response consisting of 301 × 106 = 31 906 traces. Figure   2(b) shows 51 of those traces, for 3 source positions and 17 receivers per source position. Before further processing, source-receiver reciprocity is applied, effectively doubling the number of traces, and the data are interpolated to a twice as dense spatial grid (source and receiver spacing 0.625 mm) to suppress spatial aliasing.
We denote the recorded reflection response by Green's function G(x , x, t), where x denotes the variable position of the source and x that of the receiver (actually the recorded response is the Green's function convolved with the compressed source pulse, but for the sake of simplicity we treat the recorded data as a Green's function). Consider the following variant of equation (see Appendices). This expression shows how the recorded data G(x , x, t), measured above the model, are transformed into G(r, x, t) and its time-reversal, being the response to a real source at x, observed by virtual receivers at r anywhere inside the medium. The focusing function F (x , r, t), required for this transformation, can be derived from the recorded data wave between x and r as an initial estimate of the focusing function F (x , r, t). This direct arrival, in turn, is based on an estimate of the propagation velocity of the medium (note that this does not require information about the layer interfaces). Then, by evaluating equation (4) we obtain G(r, x, t) for many virtual receiver positions r inside the medium. Next, using the retrieved virtual-receiver data G(r, x, t) in the right-hand side of equation (3), we obtain G(r, s, t) and its time-reversal, being the response to a virtual source at s, observed by virtual receivers at r. Figure 3 shows snapshots of the virtual acoustic response G(r, s, t)+G(r, s, −t), wave field through the medium, including scattering at the layer interfaces. Imperfections are explained by the finite aperture, the limited radiation angles of the piezo-electric transducers, the negligence of evanescent waves and the fact that we used a 2D method to retrieve this virtual wave field in a 3D medium.
The proposed methodology can be applied to reflection data at a wide range of scales.
Next we apply our methodology to vintage seismic reflection data, acquired in 1994 over the Vøring Basin by SAGA Petroleum A.S. (currently part of Statoil ASA). Figure 4 shows snapshots of G(r, s, t) + G(r, s, −t) obtained from these seismic data. Again, the evolution of the retrieved wave field clearly includes the primary and multiply scattered events, which have been obtained directly from the recorded reflection data. In the background these snapshots show an independently obtained seismic image of the interfaces between the geological layers, for visualisation only. Note the consistency between the position of these interfaces and the apparent origin of scattering in the snapshots.
The ability to create virtual sources and receivers inside a medium from single-sided Another exciting potential application is the investigation of induced seismicity. By acquiring high-resolution seismic reflection data in areas prone to induced seismicity, our virtual acoustic approach could forecast the wave field caused by possible future earthquakes and the associated ground motion. Moreover, when the same acquisition system is also used to passively record the response to actual induced earthquakes, our method could be used to create virtual seismometers in the subsurface around the actual earthquake and use these to retrieve accurate knowledge of the source mechanism of the earthquake, insight in the evolution of the geomechanical state of the subsurface (horizontal and vertical stress distribution, fault and fracture properties etc.), and deep understanding of the link between the earthquake and the observed ground motion.

A1 Definition of the homogeneous Green's function
Consider an inhomogeneous lossless acoustic medium with compressibility κ(x) and mass density ρ(x). Here x denotes the Cartesian coordinate vector x = (x 1 , x 2 , x 3 ); the x 3 axis is pointing downward. In this medium a space (x) and time (t) dependent source distribution q(x, t) is present, with q defined as the volume-injection rate density. The acoustic wave field, caused by this source distribution, is described in terms of the acoustic pressure p(x, t) and the particle velocity v i (x, t). These field quantities obey the equation of motion and the stressstrain relation, according to Here ∂ t and ∂ i stand for the temporal and spatial differential operators ∂/∂t and ∂/∂x i , respectively. The summation convention applies to repeated subscripts. When q is an impulsive source at x = s and t = 0, according to we find that the Green's function G(x, s, t) obeys the following wave equation Wave equation (A.5) is symmetric in time, except for the source on the right-hand side, which is anti-symmetric. Hence, the time-reversed Green's function G(x, s, −t) obeys the same wave equation, but with opposite sign for the source. By summing the wave equations for G(x, s, t) and G(x, s, −t), the sources on the right-hand sides cancel each other, hence, the function obeys the homogeneous equation Therefore G h (x, s, t), as defined in equation (A.6), is called the homogeneous Green's function.

A2 Reciprocity theorems
We define the temporal Fourier transform of a space-and time-dependent quantity p(x, t) as where ω is the angular frequency and j the imaginary unit. To keep the notation simple, we denote quantities in the time and frequency domain by the same symbol. In the frequency domain, equations (A.1) and (A.2) transform to We introduce two independent acoustic states, which will be distinguished by subscripts A and B. Rayleigh's reciprocity theorem is obtained by considering the quantity applying the product rule for differentiation, substituting equations (A.9) and (A.10) for both states, integrating the result over a spatial domain V enclosed by boundary S with outward pointing normal n i , and applying the theorem of Gauss (de Hoop 1988;Fokkema & van den Berg 1993). Assuming that in V the medium parameters κ(x) and ρ(x) in the two states are identical, this yields Rayleigh's reciprocity theorem of the convolution type We derive a second form of Rayleigh's reciprocity theorem for time-reversed wave fields. In the frequency domain, time-reversal is replaced by complex conjugation. When p is a solution of equations (A.9) and (A.10) with source distribution q (and real-valued medium parameters), then p * obeys the same equations with source distribution −q * (the superscript * denotes complex conjugation). Making these substitutions for state A in equation (A.11) we obtain Rayleigh's reciprocity theorem of the correlation type (Bojarski 1983)

A3 Representation of the homogeneous Green's function
We choose point sources in both states, according to q A (x, ω) = δ(x − s) and q B (x, ω) = δ(x − r), with s and r both in V. The fields in states A and B are thus expressed in terms of Green's functions, according to with G(x, s, ω) and G(x, r, ω) being the Fourier transforms of G(x, s, t) and G(x, r, t), respectively. Making these substitutions in equation (A.12) and using source-receiver reciprocity of the Green's functions gives (Porter 1970;Oristaglio 1989;Wapenaar 2004; van Manen et al.

2005)
where G h (r, s, ω) is the homogeneous Green's function in the frequency domain, defined as where denotes the real part. Equation (A.15) is an exact representation for the homogeneous Green's function G h (r, s, ω).
When S is sufficiently smooth and the medium outside S is homogeneous (with mass density ρ 0 and compressibility κ 0 ), the two terms under the integral in equation (A.15) are nearly identical (but with opposite signs), hence The main approximation is that evanescent waves are neglected at S (Wapenaar et al. 2011).
Using equations (A.9) and (A.13) this becomes where V (x, s, ω) = v i (x, s, ω)n i stands for the normal component of the particle velocity at x on S, due to a source at s. In the time domain this becomes where the inline asterisk ( * ) denotes temporal convolution. This is equation (1)  with a vertical axis through s and infinite radius. This cylindrical boundary exists between S 0 and S s and closes the boundary S.
The contribution of the boundary integral over S cyl in equations (A.11) and (A.12) vanishes (Wapenaar & Berkhout 1989). This implies that we can restrict the integration to the boundaries S 0 and S s . Note that n = (0, 0, −1) on S 0 and n = (0, 0, +1) on S s . On the boundaries S 0 and S s we decompose the fields in both states into downgoing and upgoing fields, according to where superscripts + and − stand for downgoing (i.e., propagating in the positive x 3 -direction) and upgoing (i.e., propagating in the negative x 3 -direction), respectively. Substituting these expressions into equations (A.11) and (A.12) and using the one-way wave equations for downgoing and upgoing waves at S 0 and S s , we obtain (Wapenaar & Berkhout 1989) Vs and (B.4) respectively. In the latter equation, evanescent waves at S 0 and S s are ignored. For state A we introduce a focusing function f 1 (x, s, ω). Here s denotes the focal point; it lies at S s which we will call the focal plane. The focusing function is defined in a sourcefree truncated medium, which is identical to the actual medium above the focal plane S s but homogeneous below this plane, see Figure S2. For x on the boundaries S 0 and S s (and in the homogeneous half-spaces above S 0 and below S s ), the focusing function is written as a superposition of downgoing and upgoing components, according to The downgoing focusing function f + 1 is incident to the inhomogeneous medium from the homogeneous half-space above S 0 , and the upgoing function f − 1 is its response. The focusing function propagates and scatters in the truncated actual medium in V s , focuses at s on S s , and continues as a downgoing wave field f + 1 in the homogeneous lossless half-space below S s . The focusing conditions at the focal plane S s are defined as Here x H and x H,A stand for the horizontal coordinates of x and s, respectively.
For state B we take again the Green's function G(x, r, ω), which, for x at S 0 and S s , is written as Here r can be chosen anywhere below the upper boundary S 0 . When r lies above s (and below S 0 , as illustrated in Figure S1), it is by definition situated in V s . When r lies below s, it is situated outside V s . Because the upper half-space above S 0 is homogeneous, we have [G + (x, r, ω)] x 3 =x 3,0 = 0. (B.9) Substituting p A (x, ω) = f 1 (x, s, ω), p ± A (x, ω) = f ± 1 (x, s, ω), q A (x, ω) = 0, p ± B (x, ω) = G ± (x, r, ω) and q B (x, ω) = δ( This is equation (3) in the main paper.
Note that the Green's function G(r, x, ω) on the right-hand side of equation (B.16) can be obtained from a similar representation. To see this, replace in the right-hand side of equation (B.16) S 0 by S 0 just above S 0 , replace x on S 0 by x on S 0 , r inside the medium by x on S 0 and s by r. This gives a representation for G h (x, r, ω). Using source-receiver reciprocity we finally get In the time domain this becomes

(B.19)
This is equation (4) in the main paper, where for simplicity the prime in S 0 is dropped.