Retrieving reflection responses by crosscorrelating transmission responses from deterministic transient sources: Application to ultrasonic data

The cross correlation of observed wave fields (acoustic or elastic) at two points retrieves the Green’s function at one of those points as if there were an impulsive source at the other. This was proven by researchers in different fields of science using different underlying theories.^1–7 Some of the theories require that the recorded wave fields are diffuse. These diffuse wave fields can be produced by uncorrelated noise sources and∕or by multiple scattering inside or at the boundaries of a medium. Until now, most of the applications to measured data made use of diffuse fields. By crosscorrelating recorded coda waves and ambient noise, researches could retrieve surface waves,^8,9 diving $P$-waves,¹⁰ and reflected body waves.¹¹ It was shown that for good retrieval results one requires recordings several orders longer than the normal recordings between the points of interest would be. This long recording time is required to reduce the amount of correlation between the different noise sources as much as possible or, in case of multiple scattering, to approximate the equipartitioned state as well as possible.

Other theories do not require diffuse wave fields, but instead make use of multitude of independent deterministic transient source responses. One of the advantages of using this method is that the recordings to be cross correlated do not need to be disproportionally long: The recording length is dictated by the expected arrival time of the latest event one wants to retrieve. The number of the recorded responses from transient sources multiplied by the length of each recording will still be much shorter than the required recording time in case of diffuse wavefields. To retrieve the Green’s function between two receivers in an open system, these receivers should be completely surrounded by sources.^12,13 When this is not the case, spurious arrivals may appear.¹⁴ To avoid this, a waveguide can be utilized.^15–17 Another solution is to retrieve parts of the Green’s function from specific source-receiver geometries. One such possibility is to retrieve the reflection response at a free surface from the cross correlation of transmission responses from sources illuminating the receivers from one side only. This can be very useful, for example, in seismology, where one can crosscorrelate transmission responses from subsurface sources (earthquakes) measured at the surface along receiver arrays to retrieve the earth’s reflection response as if from sources at the surface. Having reflection responses would allow the application of advanced imaging techniques to obtain images of the subsurface with a higher resolution than those obtained by tomographic inversion of the transmission data. In the following, we test the method of retrieval of reflection responses through cross correlation of transmission data from deterministic transient sources. We apply the method to ultrasonic data measured on an inhomogeneous granite sample.

Let one have three-dimensional (3D) inhomogeneous, lossless, and source-free medium bounded by a free surface. In this medium, one chooses a domain $D$ enclosed by a boundary $∂D=∂D0∪∂D1$ with an outward-pointing normal vector. The boundary $∂D0$ coincides with (a part of) the free surface, while $∂D1$ is an arbitrary-shaped surface inside the medium. If there would be an impulsive force source $f$ acting in the $xq$-direction (⁠$q=1$⁠, 2, or 3) at $xB=(xB,1,xB,2,xB,3)$⁠, then at $xA=(xA,1,xA,2,xA,3)$ one could record the impulse response, i.e., Green’s function of the medium $Gp,qv,f(xA,xB,t)$⁠, representing the particle velocity $v$ observed in the $xp$-direction. If $xA$ and $xB$ are chosen in $D$ or on $∂D0$⁠, and the medium along and outside $∂D1$ is taken homogeneous with propagation velocity $cP$ for $P$-waves and $cS$ for $S$-waves and density $ρ$⁠, then it can be shown that¹³ 

As the Green’s function is a causal function of time, it does not interfere with its time-reversed version and one can simply take the causal part in the left-hand side of Eq. (1) to obtain the Green’s function. In the right-hand side, the Green’s functions denote observed particle velocities in the $xp$- and $xq$-directions at $xA$ and $xB$⁠, respectively, due to impulsive sources at $x=(x1,x2,x3)$ along $∂D1$⁠; the superscripts $P$ and $S$ denote sources of $P$-waves and sources of $S$-waves with different polarizations along the three axes; these polarizations are denoted by the subscript $k$ and a summation is carried out over this subscript according to Einstein’s summation convention for repeated subscripts. The asterisk denotes convolution, but as one convolves a signal with a time-reversed signal, one performs correlation. Hence, Eq. (1) shows that to retrieve the Green’s function between the observation points $xA$ and $xB$⁠, one needs at these points separate recordings from $P$- and $S$-wave sources. In actual implementations, this would require decomposition of the recorded wave fields into wave fields from $P$- and $S$-wave sources. When this is not done, spurious events will appear in the retrieved results.

To test Eq. (1) we use a controlled laboratory experiment with an Oshima granite block with dimensions $300×300×100mm$⁠, in which there is a cylindrical hole of $150mm$ in length and $15mm$ in diameter filled with epoxy [see Fig. 1(a)]. The Oshima granite consists mainly of quartz, plagioclase, and biotite with grain sizes mostly from $1to5mm$⁠. The estimated average $P$-wave and $S$-wave propagation velocities for the block are 4500 and $2700m∕s$⁠, respectively. The estimated $P$-wave and $S$-wave velocities for the epoxy are 2300 and $1400m∕s$⁠, respectively. For simplicity, we consider a two-dimensional (2D) experiment, using observation and source line arrays at opposite sides of the block lying in a plane perpendicular to the long axis of the cylinder. This means that in the second integral in Eq. (1), $k$ takes on only the value of 2 meaning that we will need to record the response from $P$- and $SV$-sources. The observation array of 101 measurement points is placed along the free surface on the front wall of the granite sample; the distance between the observation points is $1mm$⁠. At these points we measure only the normal component of the particle velocity. This means that in Eq. (1) the subscripts $p$ and $q$ take on the value of 3 [note the coordinate system in Fig. 1(a)]. The measurements at the observation points are performed using a laser Doppler vibrometer (LDV). Details about the LDV can be found in Nishizawa et al.¹⁸ and in Sivaji et al.¹⁹ 

Compared to measurements in the field, in a controlled laboratory experiment one can use ultrasonic compressional and shear transducers to approximate separate $P$- and $S$-wave sources. The source array, consisting of 21 source points with $5mm$ distance between the individual points, is placed along the back wall of the granite block. The sources are fed with a single-cycle $250kHz$ sine wavelet. For this frequency, the wavelengths for $P$-waves and $S$-waves are 18 and $11mm$⁠, respectively. At this scale, the granite sample is effectively heterogeneous and both compressional and shear waves are diffracted from the grains. Due to the finite diameter $(5mm)$ of the source transducers, their radiation pattern will not be like the desired point-source pattern.^18,20 This means that the compressional transducer will produce some amount of $S$-wave energy, especially at larger angles. Still, for the chosen source-receiver geometry the main energy will be radiated in the form of $P$-waves. A similar effect holds for the shear transducers, which will produce some $P$-wave energy.

Using stationary phase analysis¹³ it can be shown that for the chosen 2D setup, sources placed along the walls perpendicular to the observation array would contribute strongly to the retrieval of surface waves, while sources placed along the wall opposite to the observation array (as in our experiment) will mainly contribute to the retrieval of reflection arrivals. We choose to concentrate on retrieving reflection arrivals as they can further be used to obtain a better image of the interior of the sample.

We perform the following measurements. A source transducer is placed and ignited at a certain source point and the resulting wave fields are measured at a point along the observation array. The electronic circuit of the LDV causes incoherent high-frequency noise to appear in the recorded wave forms. To eliminate this noise, we repeat the ignition of the source transducer at the same position 2000 times and at the observation point the 2000 recorded wave forms are summed to obtain a final wave form recording at that observation point. By placing a compressional transducer and recording the wave forms at all observation points in the way described earlier, we obtain a $P$-source transmission panel. Similarly, placing a shear transducer and repeating the measurements, we obtain a $S$-source transmission panel. This procedure is repeated for each of the 21 source points. Fig. 1(b) shows an example $P$-source transmission panel (left) with a clear $P$-wave direct arrival with an apex at about $0.024ms$ and an example $S$-source transmission panel (right) with a clear $P$-wave direct arrival and a clear $S$-wave direct arrival with an apex at around $0.039ms$⁠. The transmission panels are $0.164ms$ long aiming at recording the first few multiples of the direct $P$- and $S$-wave arrivals, i.e., the arrivals that have been reflected by both the front and the back wall in Fig. 1(a). Even though the later-arriving conversions and multiples are not readily interpretable due to scattering, parts of the first and the second multiple arrivals of the direct $P$-wave can be observed in the left panel in Fig. 1(b) (with estimated apices around 0.07 and $0.115ms$⁠, respectively) and parts of the first multiple of the direct $S$-wave can be seen in the left-hand panel in Fig. 1(b) (with an estimated apex around $0.115ms$⁠). On the panels we also see inclined linear arrivals caused by the presence of the left and right walls in Fig. 1(a). After cross correlation and integration [Eq. (1)], such events would contribute to retrieval of similar inclined events, which interfere with retrieved reflections. For this reason, we apply a frequency–wave number $(f–k)$ filter to eliminate the inclined events. In the $P$-wave response, we also see nearly horizontal arrivals around 0.08 and $0.125ms$ that arise from the presence of the to wall in Fig. 1(a) and representing 3D (out-of-plane) events. These arrivals cannot be filtered out as their behavior in the $f–k$ domain is very similar to the direct transmission arrivals and their free-surface multiples. After cross correlation and integration, such out-of-plane events will cause nonphysical events to appear in the retrieved results.

By placing a source transducer at a position of an observation point and recording the resulting wave fields along the observation array, we obtain a reflection common-source panel. Such a panel will be used to verify possible retrieved reflections.

To implement Eq. (1), we do the following. We choose an observation point at a position where we would like to have a retrieved source, say $xB$⁠. We take one $P$-source transmission panel [left-hand panel in Fig. 1(b)] and extract the trace $G3v,P(xB,x,t)$ at the position of $xB$ (the source is at $x$⁠). Then we correlate the trace with the whole panel, representing $G3v,P(xA,x,t)$ for a variable $xA$ and fixed $x$⁠, and the result is deconvolved with a wavelet extracted around zero time from the autocorrelation trace. The result, representing the first integrant in Eq. (1), is a correlation $P$-source panel. We repeat this for all $P$-source transmission panels, i.e., for all $x$⁠, by always correlating with the trace at the chosen position $xB$⁠. Then, according to the first integral in Eq. (1), we sum the individual correlation $P$-source panels (i.e., integration along the source position $x$⁠). The previous procedure is repeated for all the $S$-source transmission panels as well. This is followed by summation of the outputs from the two integrals and extraction of the causal part – the result is a reflection common-source panel for a retrieved source position at $xB$⁠.

After a reflection common-source panel is retrieved, it should be compared to a directly observed reflection common-source panel. The left-hand panels in Figs. 2(a) and 2(b) show the observed reflection common-source panel for a compressional transducer (a source is represented by a star) at 25 and $50mm$⁠, respectively, after $f–k$ filtering of the surface waves and deconvolution. The removal of the surface waves is not a trivial task as the inhomogeneities close to the surface cause dispersion. Due to the physical presence of the source transducer the near offsets could not be recorded. The middle pictures in Figs. 2(a) and 2(b) show the retrieved result using Eq. (1) for retrieved source positions at 25 and $50mm$⁠, respectively. The presence of inhomogeneities in the Oshima granite makes both pictures look cluttered and makes the interpretation and comparison of reflection events difficult. To be able to interpret the main events and compare the panels, we performed 2D elastic numerical modeling with a finite-difference scheme using a homogeneous background velocity model with a circle in it representing the cross section of the epoxy cylinder. The dimensions and velocities of the model were the same as the dimensions and the estimated average velocities of the granite block. The right-hand pictures in Figs. 2(a) and 2(b) show the modeling results (without the surface waves) for the same source positions. Note that these modeling results should only be used for interpretation of expected reflection arrivals. The modeling results are not representative for the clutter due to the inhomogeneous background of the Oshima granite block. This means that a reflection arrival, visible in a modeled common-source panel, might be (partly) obscured by diffractions in the observed and the retrieved common-source panels.

In the numerically modeled results in Fig. 2, the most prominent events are the top-of-cylinder reflection arrivals with apices around $0.018ms$⁠, the bottom-of-cylinder reflection arrivals with apexes around $0.035ms$⁠, and the reflection arrivals from the granite’s backwall, which arrive around $0.045ms$⁠. The backwall reflection event is followed by a train of conversions. Due to the missing offsets and the remnants of the surface waves, it is hard to interpret the top-of-cylinder reflection in the directly observed reflection common-source panels. Parts of this reflection might be present for the source position at $50mm$⁠. Another factor making the interpretation of the top-of-cylinder reflection difficult is the interference with linear events in the directly observed data coming from the left and right walls of the granite block and remaining after the $f–k$ filtering. Such events are easier to remove in the transmission panels and after cross correlation and summation their effect on the retrieved results is much smaller. On the retrieved reflection common-source panels, the apices of the top-of-cylinder reflection are readily visible and the wings of the reflection can be traced from $0mm$ until the apex in Fig. 2(a) and from the apex until about $85mm$ in Fig. 2(b). The backwall reflection in the retrieved common-source panels is at least as good as the reflection in the observed common-source panels. In Fig. 2(b) the retrieved backwall reflection is even more continuous than in the directly observed reflection panel. The bottom-of-cylinder reflections are not interpretable in both the directly observed and the retrieved common-source panels. In the retrieved reflection common-source pictures, we can see linear nonphysical events, like the ones at around $0.05ms$⁠, which are a result from the presence of the out-of-plane events in the transmission data.

Having retrieved reflection common-source panels at the position of every receiver, we can apply advanced imaging techniques to look into the interior of the granite sample. The left-hand panel in Fig. 3 shows the imaging result from the application of prestack depth migration to the retrieved reflection events using a homogeneous model with a velocity of $4500m∕s$⁠, i.e., the estimated average propagation velocity of $P$-waves through the granite sample. The image is quite cluttered, as it can be expected for an image of an inhomogeneous sample, and because of this it is difficult to describe it in terms of imaged grains of the granite sample, imaged epoxy cylinder, imaged back wall of the sample, and artifacts due to the out-of-plane events. To verify if we could image grains in the sample, we would need directly observed reflection measurements from a source at each of the receiver positions, which we do not have. To facilitate the interpretation of the left-hand image in Fig. 3, in the right-hand panel in Fig. 3 we show the obtained image from the migration of the numerically modeled results. This image can be used to identify only the main reflectors in the granite block—the epoxy cylinder and the backwall of the granite, but is not representative for imaged inhomogeneities (grains) inside the granite block. Comparing the two images in Fig. 3, we see that by using only a few sources at the backwall of the block, we have successfully imaged the top of the cylinder and the sample’s backwall.

We showed how one can retrieve the Green’s function between two receivers from the cross correlation of responses from deterministic transient sources. We showed results from the application of this method to ultrasonic data from an inhomogeneous granite block with an epoxy cylinder in it. We successfully retrieved reflection arrivals from the cross correlation of transmission measurements from separate compressional and shear sources. The retrieved reflection results allowed us to use an advanced imaging technique. We applied imaging using a homogeneous velocity model and successfully identified the main reflectors in the granite sample, i.e., the top of the epoxy cylinder and the back wall of the granite sample.

This research is supported by The Netherlands Research Centre for Integrated Solid Earth Sciences ISES; the Technology Foundation STW, Applied Science Division of The Netherlands Organization for Scientific Research NWO and the technology program of the Ministry of Economic Affairs (Grant No. DTN.4915). The measurement of the data was partly sponsored by NWO through a travel grant. The authors would like to thank Dr. J. Spetzler for helping with the measurements and Dr. X. Campman for the very helpful discussion.