(PDF) Near-surface diffractor detection at archaeological sites based on an interferometric workflow

GEOPHYSICS, VOL. 86, NO. 3 (MAY-JUNE 2021); P. WB1–WB11, 12 FIGS., 1 TABLE. 10.1190/GEO2020-0142.1 Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms Near-surface diffractor detection at archaeological sites based on an interferometric workflow Jianhuan Liu1, Deyan Draganov1, Ranajit Ghose1, and Quentin Bourgeois2 ABSTRACT revealed by a combination of seismic interferometry and nonsta- tionary adaptive subtraction, and then further enhanced through Detecting small-size objects is a primary challenge at crosscoherence-based supervirtual interferometry. A diffraction archaeological sites due to the high degree of heterogeneity image is then computed by a spatial summation of the revealed present in the near surface. Although high-resolution reflection diffractions. We use the phase-weighted stack to enhance the DOI:10.1190/geo2020-0142.1 seismic imaging often delivers the target resolution of the sub- coherent summation of weak diffraction signals. Using synthetic surface in different near-surface settings, the standard process- data, we show that our scheme is robust in locating diffractors ing for obtaining an image of the subsurface is not suitable from data dominated by strong Love waves. We test our method to map local diffractors. This happens because shallow seismic- on field data acquired at an archaeological site. The resulting dis- reflection data are often dominated by strong surface waves that tribution of shallow diffractors agrees with the location of anoma- might cover weaker diffractions and because traditional common- lous objects identified in the V S model obtained by elastic SH/Love midpoint moveout corrections are only optimal for reflection full-waveform inversion using the same field data. The anomalous events. We propose an approach for imaging subsurface objects objects correspond to the position of a suspected burial, also de- using masked diffractions. These masked diffractions are first tected in an independent magnetic survey and corings. INTRODUCTION events is usually weak (Klem-Musatov, 1994) and difficult to detect in field seismic data due to the usually dominant presence Archaeology is the study of past human cultures through the of other coherent events such as surface waves (SWs) and specular analysis and interpretation of artifacts and material remains (Smith, reflections. 2014). The material remains, which possess certain physical pro- Several methods have been developed to detect various near- perties (e.g., elastic impedance) that contrast with the subsurface surface buried features using diffracted waves. Landa and Keydar background medium, can be detected by means of noninvasive (1998) develop an algorithm to construct a so-called diffraction- near-surface geophysical surveys and more specifically with the point-section (D-section) by concentrating the diffracted signals seismic methods. For detecting small localized objects, the usual spread from diffractor points. In this D-section, high-amplitude normal-moveout stacking is generally not useful because traditional anomalies indicate the existence of local heterogeneities. Shtivel- common-midpoint moveout corrections are only optimal for reflec- man and Keydar (2005) present a multipath summation approach tion events. Although common-offset gathers can potentially show to stack diffracted signals along all possible diffraction trajectories such objects more reliably (Ghose et al., 1998), these gathers focus for detection of shallow inhomogeneities. Kaslilar et al. (2013) on the primary reflected waves and mostly ignore the diffracted propose an approach inspired by seismic interferometry (SI) to es- waves. The conspicuous appearance of diffracted waves can be timate the location of a near-surface tunnel by traveltime inversion used to locate an object buried in the heterogeneous subsurface of crosscorrelated diffracted arrivals. The crosscorrelation pro- at an archaeological site. However, the amplitude of these diffracted cedure eliminates the common raypath between the source and Manuscript received by the Editor 5 March 2020; revised manuscript received 10 September 2020; published ahead of production 19 January 2021; published online 11 March 2021. 1 Delft University of Technology, Department of Geoscience and Engineering, Stevinweg 1, Delft 2628 CN, The Netherlands. E-mail:

(corresponding author);

;

. 2 Leiden University, Faculty of Archaeology, Einsteinweg 2, Leiden 2333 CC, The Netherlands. E-mail:

. © 2021 Society of Exploration Geophysicists. All rights reserved. WB1 WB2 Liu et al. the subsurface diffractor. This makes the traveltime dependent only acronym used in this paper in Table 1. Next, we present the theory on properties between the receivers and the diffractors. of each of these steps. Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms In this paper, we introduce a workflow for imaging subsurface objects using masked diffractions, that is, weak diffractions covered SI + AS for SW suppression by other signals. This workflow consists of three main steps. We first reveal masked diffractions by suppression of the dominant SI refers to a process of retrieving the seismic response between Love waves through a combination of SI and nonstationary adap- two receivers by crosscorrelating and integrating the wavefields re- tive subtraction (AS) (Dong et al., 2006; Halliday et al., 2007; corded at these receivers from a boundary of sources (Wapenaar Konstantaki et al., 2015; Liu et al., 2018). We then enhance the re- et al., 2008; Schuster, 2009). When the sources are located along vealed weak diffraction signal through crosscoherence-based super- the earth’s surface, the response retrieved by SI would be dominated virtual interferometry (SVI) (Dai et al., 2011; Nakata et al., 2011; by SWs (Dong et al., 2006; Halliday et al., 2007). The dominance of An and Hu, 2016; Place et al., 2019). Finally, a diffraction image, SWs is observed in passive seismology (Shapiro et al., 2005), ex- which can be used to interpret the locations of subsurface diffrac- ploration seismology (Dong et al., 2006; Balestrini et al., 2020), and tors, is generated by a spatial summation (Shtivelman and Keydar, near-surface seismology (Konstantaki et al., 2015). Depending on 2005; Landa et al., 2006; Shtivelman et al., 2009) of enhanced the source type (passive or active), the retrieval steps for SI are dif- diffractions. ferent. In the context of this paper, we focus only on active-source In the following, we first describe the theory of each of the above- interferometry. mentioned steps. We then illustrate the feasibility of our approach in The steps of interreceiver SW retrieval using active sources are detecting near-surface objects using synthetic seismic data domi- illustrated in Figure 2a. In Figure 2a, ↝ represents SWs propagat- nated by Love waves. Finally, we test our workflow on field seismic ing along the earth’s surface from an active source to the receivers. data acquired at an archaeological site in the Veluwe, the Nether- By crosscorrelating the trace recorded at receiver A with that at B, lands, with the aim of mapping the distribution of shallow diffrac- we obtain a virtual trace at B as if caused by a virtual source located tors of archaeological importance. at A. The traveltime of this virtual trace is denoted as τAB (red ↝ in Figure 2a), and it is independent of the source position as long as the source falls inside the stationary-phase region (Snieder, 2004). For DOI:10.1190/geo2020-0142.1 METHODOLOGY this reason, the virtual trace at B from each of the sources can be In this section, we present a workflow (Figure 1) for near-surface stacked constructively to give the interreceiver estimate of the SWs diffractor detection from data dominated by strong SWs (in our propagating from A to B. Mathematically, it can be formulated in case, Love waves), which consists of three main steps. These steps the frequency domain as (Wapenaar and Fokkema, 2006; Halliday include SWs suppression by SI + AS, diffraction enhancement et al., 2007) through SVI, and diffraction imaging. We list the full names of each X N CX B X A ¼ u ðXA ; X i ÞuðX B ; X i Þ: (1) i¼1 In equation 1, CXB XA represents interreceiver SWs propagating from A to B retrieved by SI. The seismic data generated by a source at Xi and recorded by receivers at X A and X B are denoted as uðXA ; Xi Þ and uðXB ; Xi Þ, respectively. The superscript indicates complex conjugation, and N is the number of active sources avail- able for stacking. If we ignore the presence of noise, in the fre- quency domain the seismic wavefield also can be represented by the multiplication of a source wavelet and a Green’s function: uðX; Xi Þ ¼ WðX i ÞGðX; X i Þ; (2) Table 1. A list of acronyms used in this paper. Acronyms Full names AS Adaptive subtraction f-k Frequency wavenumber FWI Full-waveform inversion PWS Phase-weighted stack Figure 1. Flowchart for raw-data processing for diffraction imag- ing. To reveal weak diffractions, a combination of SI and AS is first SI Seismic interferometry used to suppress the dominating SWs. A crosscoherence-based SVI SVI Supervirtual interferometry is then applied to further enhance the diffractions. Near-surface diffractor detection WB3 where WðXi Þ is the source wavelet generated by the source located In this paper, we use crosscoherence instead of crosscorrelation at Xi and GðX; Xi Þ denotes the Green’s function between Xi and X. (similar to equation 1) to retrieve the virtual diffraction HXB XA Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms Substituting equation 2 into equation 1, we have (Nakata et al., 2011; Place et al., 2019): X N X N u ðXA ; X i ÞuðX B ; X i Þ CXB XA ¼ W 2 ðX i Þ G ðX A ; X i ÞGðX B ; Xi Þ: (3) H XB XA ¼ i¼1 i¼1 ju ðX A ; X i ÞjjuðxB ; Xi Þj þ η X N G ðX A ; X i ÞGðXB ; X i Þ To accurately retrieve the amplitude of the seismic response using ¼ ; (4) i¼1 jG ðX A ; X i ÞjjGðxB ; X i Þj þ η equation 3, according to the theory, the sources should effectively surround the receivers and illuminate them homogeneously, the source boundary should be a sphere with sufficiently large radius, where η is a small number for stabilizing the division and j · j in- and the medium should be lossless (Draganov et al., 2006; Wape- dicates that the amplitude spectrum is used. If we consider the effect naar and Fokkema, 2006). For conventional 2D near-surface seis- of the source wavelet, which is denoted by WðX i Þ, a virtual diffrac- mic surveys, this requirement cannot be fulfilled because active tion obtained by crosscorrelation will contain an amplitude factor sources are placed only at the surface and usually along the line proportional to W 2 ðX i Þ (see equation 3). The corresponding super- connecting the receivers (this limitation is also true for 2D seismic virtual diffraction will then include an amplitude factor proportional surveys at the exploration scale); thus, the retrieved SWs will be to W 3 ðX i Þ, due to the additional convolution involved in the second characterized by amplitude errors. step. This W 3 ðX i Þ factor may lead to wavelet distortion of the en- To account for the amplitude difference between the retrieved hanced diffractions. However, SI by crosscoherence (equation 4) SWs and the dominant SWs from the physical, active source, a eliminates (theoretically) the contribution of the source wavelet matching filter (see Appendix A) is estimated via the regularized WðXi Þ and makes the virtual diffraction contain only the medium nonstationary regression technique proposed by Fomel (2009). This response. The resulting supervirtual diffraction will then contain is done by using shaping regularization (Fomel, 2007) to constrain an amplitude factor proportional to WðX i Þ, which avoids wavelet DOI:10.1190/geo2020-0142.1 the continuity and smoothness of the filter coefficients. The re- distortion and thus improves the precision of the reconstructed trieved SWs are then convolved with this estimated matching filter diffractions. and are subsequently subtracted from the field data. Diffraction imaging SVI for diffraction enhancement The diffractions that are revealed by SI + AS and enhanced by After the suppression of the dominant SWs, the hidden diffrac- SVI can now be used for detecting near-surface diffractors. We de- tion events may become identifiable. To further enhance the diffrac- sign a diffraction algorithm (based on studies from Landa and Key- tion energy, we use a crosscoherence-based SVI (Dai et al., 2011). dar [1998] and Shtivelman and Keydar [2005]) to coherently focus The principle of SVI consists of two steps and is illustrated in the diffraction energy back to its original position; all other events in Figure 2b and 2c. First, the traces recorded at receivers A and B are crosscorrelated. The a) common raypath from the source to a subsurface object (the black straight line) is subtracted; thus, the traveltime of the obtained virtual diffraction can be denoted as τOB þ τOA . This traveltime is the same for all active-source positions of a b) survey, so stacking the virtual diffraction at B from N sources will enhance pffiffiffiffi its signal-to-noise ratio (S/N) by a factor of N assuming uncorre- lated noise. Next, SI by convolution is applied (Slob et al., 2007) and the stacked virtual diffrac- tion is convolved with an actual trace originally c) recorded at a receiver position A from the source at X to produce a supervirtual trace at B. This supervirtual trace represents the seismic response from a subsurface diffractor, recorded by receiver B from the source positioned at X. The supervir- tual trace is kinematically equivalent (with the Figure 2. (a) The step for retrieving dominant SWs between two receivers by SI. (b) and traveltime equal to τXO þ τOB ) for all receiver (c) The steps for creating supervirtual diffractions with an increased S/N. Note that the positions A located between the source X and SW component also can be retrieved/enhanced by the SVI step. Hence, it is beneficial to receiver B. Thus, the retrieved supervirtual traces apply SVI for diffraction L enhancement only after suppression of the dominant SWs. The at N receiver positions A could be stacked to ob- symbols ⊗, ⊙, and denote the crosscorrelation-, crosscoherence-, and crossconvo- lution-based operators, respectively. The symbol ↝ represents SWs propagating along tain a final tracepwith ffiffiffiffi the S/N increased again by the surface from an active source to the receivers. The symbol → denotes the propa- as much as Oð N Þ. gation of diffracted waves from an active source to the receivers. WB4 Liu et al. the shot-gather domain will be suppressed by the used directional and destructive interference of the amplitudes contributing along summation. Such a procedure will make the true diffractions appear each diffraction trajectory. Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms as high-amplitude anomalies in the resulting sections, which can be The practical implementation steps of this diffraction-stacking used to indicate the locations of the diffractors. method include the following. For a specific common-source Consider a seismic wave emitted from a source X0 (Figure 3). gather, we first assume that the diffractor is located directly under When it encounters a subsurface diffractor along its propagation the first receiver and the depth of the diffractor is known as R. This path, it would generate a secondary wave that spreads from this will define a specific diffraction traveltime curve according to equa- point in all directions. In Figure 3, X i is the image point and Ω tion 5. Next, we obtain a diffraction-moveout corrected gather by denotes the diffracted wavefront element surrounding this image applying the diffraction-moveout correction to this common-source point. We assume that the velocity variations are small and the gather. Such a traveltime correction is repeated for all other combi- propagation distance between source and receiver is relatively short, nations of source-receiver pairs. We then resort these moveout-cor- so that the diffracted wavefronts can be approximated as an arc of a rected common-source gathers into the common-receiver domain circle with a radius of R. This radius R has a specific physical mean- and stack traces within each common-receiver gather to produce ing, which is the depth of the diffractor below the image point. one trace per receiver position. To enhance coherent summation Thus, the kinematic response (T sr ) of the diffractor recorded by of weak diffraction signals, a phase-weighted stack (PWS) (see Ap- receivers can be written as (Landa and Keydar, 1998) pendix B) method is used. The stacked traces from each common- receiver gather are then assembled into a diffraction image (Landa pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and Keydar, 1998; Liu et al., 2019). For better visualization of high- ðX 0 − X i Þ2 þ R2 þ ðX k − X i Þ2 þ R2 2R T sr ¼ − þ 2td ; amplitude anomalies, we produce the final diffraction image using a V0 V0 coherency measure calculated from each diffraction image resulting (5) from different R. The coherency measure that we choose is the un- normalized crosscorrelation sum (Neidell and Taner, 1971), due to its high sensitivity to coherent weak signals. This coherency func- where V 0 is the average velocity of the medium above the diffractor tion (C) can be expressed as and td is the vertical time above the diffractor. The terms X 0 , Xi , and XX 2 DOI:10.1190/geo2020-0142.1 X k denote the lateral positions of the source (s), image point, and M receiver (r), respectively. C¼ f i;tðiÞ ; (6) If we know the exact value of R (the depth of the diffractor), a t i¼1 diffraction image can be obtained by stacking the seismic energy along the diffraction surface defined by equation 5. If the value where f i;tðiÞ denotes the amplitude value of the ith trace at two-way of R is unknown, it can be estimated by maximizing the semblance traveltime tðiÞ and M is the number of traces. The outer summation function (Taner and Koehler, 1969) calculated from seismic records is performed over the two-way zero-offset time samples t within a within a time window along the traveltime surface defined by equa- time gate. The length of the gate should be approximately equal to tion 5 (similar to the traditional velocity analysis). The alternative the main wavelength of the seismic signal. approach for diffraction imaging, not requiring the specification of the radius R, is the multipath summation (Shtivelman and Keydar, 2005; Shtivelman et al., 2009). This is done by stacking seismic SYNTHETIC TEST energy (with unit weights) along diffraction trajectories defined To verify the effectiveness of our method for near-surface diffrac- by equation 5, calculated for various values of R within a specified tor detection, we first test it on data from 2D synthetic modeling. range. The resulting diffraction image would be close to the one The model shown in Figure 4 is a shear wave (S-wave) velocity (V S ) produced by stacking with the correct radius, due to the constructive model. This model consists of two layers. The first layer has a lower velocity (100 m/s) and thickness of 0.5 m. Below this layer, we place a high-velocity (150 m/s) half-space. Two circular diffractors are also embedded in the model at a depth of 5 m. These diffractors have a certain impedance contrast with respect to the background Location (m) 15 20 25 0 250 S-wave velocity (m/s) 2 Depth (m) 200 4 6 150 8 10 100 Figure 3. Schematic view of diffracted waves. The terms X 0 , Xi , and Xk denote the lateral position of the source (the red star), image Figure 4. The S-wave velocity (V S ) model used to generate the syn- point (the blue triangle), and receiver (the black triangle), respec- thetic common-source gathers. Two shallow diffractors with a given tively. The term Ω is the diffracted wavefront element surrounding impedance contrast with respect to the background media are in- the image point (the blue triangle at the lateral position Xi ). cluded in the model. Near-surface diffractor detection WB5 medium. The size of these diffractors is 1.0 m. The receiver array, equation 3, such as elastic media and active sources distributed only which is located at the surface, consists of 40 geophones aligned in at the surface, causes the amplitude of the estimated Love waves in Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms the horizontal direction from 10 to 29.5 m at an interval of 0.5 m. Figure 5b to be erroneous compared to the amplitude of the Love During data generation, the receiver array is kept fixed whereas the waves in the original shot gather in Figure 5a. source, also deployed at the surface, advances at a step of 1 m. The To account for the amplitude, phase, and frequency differences first source is located 5 m to the left of the first geophone, whereas between Figure 5a and 5b, a matching filter is estimated via the the last source is placed 4.5 m to the right of the last geophone. With regularized nonstationary regression technique proposed by Fomel this acquisition geometry, 30 common-source gathers are computed (2007). Figure 5d shows the mean coefficients of the matching filter using an elastic 2D SH finite-difference algorithm (Bohlen, 2002; determined by minimizing the difference between Figure 5a and 5b Dokter et al., 2017). We discuss the reason for choosing SH-wave in the least-squares sense. From Figure 5d, we can see that the filter processing in the “Field example” section. At the top boundary, a coefficients vary in time and space, which agrees well with the vari- free-surface boundary condition is realized by the image technique ability of the original shot gather (Figure 5a). We then convolve the (Robertsson, 1996) for accurate SH-wave mod- eling. For the other boundaries, convolutional perfectly matched layers absorbing boundary a) Receiver number b) Receiver number 10 20 30 40 10 20 30 40 conditions are used (Komatitsch and Martin, 0 0 2007). The source signature is a band-limited spike with the central frequency at 40 Hz. 0.05 0.05 Revealing weak diffractions Time (s) 0.10 0.10 Figure 5a shows an example of the synthetic shot gathers for the source located at horizontal position 13 m, that is, at the seventh receiver. 0.15 0.15 This shot gather (Figure 5a) is dominated by dis- DOI:10.1190/geo2020-0142.1 persive Love waves because of the presence of 0.20 0.20 the high-velocity half-space in the velocity model. The amplitudes of the diffracted waves c) 10 20 30 40 d) 10 20 30 40 0 0 in Figure 5a, which represent the seismic re- sponse from the diffractors embedded in the 0.006 model in Figure 4, are so weak that it is difficult 0.05 0.05 to identify them directly. 0.004 Time (s) Before we introduce our SI + AS scheme for 0.10 0.10 0.002 the suppression of Love waves, we first apply 0 conventional frequency-wavenumber (f-k) filter- –0.002 ing to eliminate the dominant Love waves, while 0.15 0.15 –0.004 preserving the weak diffraction signals. Figure 6a displays the f-k spectrum of the synthetic gather 0.20 0.20 from Figure 5a. We can see that this spectrum is dominated by two clusters of Love-wave energy Figure 5. Steps for revealing weak diffractions dominated by strong Love waves. (a) A (the black arrows). We design a fan filter to reject synthetic SH shot gather computed for the model shown in Figure 4, (b) retrieved Love waves from the gather as shown in Figure 5a, (c) result after convolution of the data in the Love-wave energy in the f-k domain (Fig- (b) with the matching filter (d), and (d) mean filter coefficients estimated by minimizing ure 6b). This is followed by inverse mapping the difference between (a) and (b), using the regularized nonstationary regression back to the shot-gather domain. As shown in method. the resulting filtered shot gather (Figure 7a), weak diffraction signals begin to be identifiable, due to the signifi- cant suppression of the Love waves. However, we can still observe a) Wavenumber (1/m) b) Wavenumber (1/m) some remaining Love waves (the red box in Figure 7a). This is be- –1.0 –0.5 0 0.5 –1.0 –0.5 0 0.5 0 0 cause Love waves and diffractions are not well separated in the f-k domain (Figure 6a) due to their similar apparent velocities. Thus, it 20 1.2 20 0.12 is difficult to design an efficient fan filter that can completely reject 1.0 40 0.10 Frequency (Hz) 40 Amplitude the Love-wave energy while preserving the diffraction events. 0.8 0.08 60 60 We then make use of SI to compute a virtual common-source 0.6 0.06 80 0.4 80 0.04 gather for the receiver positioned at 13 m (this receiver becomes 100 0.2 100 0.02 the virtual source), following the scheme presented in Figure 2a. 0 0 As shown in Figure 5b, the main kinematic characteristics of the 120 120 Love waves in Figure 5a are retrieved well. However, due to the source term W 2 ðX i Þ involved in the SI procedure (equation 3), Figure 6. (a) The f-k spectrum of the synthetic SH shot gather from the wavelet of the Love waves in Figure 5b is slightly broader than Figure 5a, and (b) the f-k spectrum of the same shot gather, but after that in Figure 5a. Further, the interferometric approximation used in rejecting Love-wave energy indicated by the black arrows in (a). WB6 Liu et al. estimated matching filter (Figure 5d) with Figure 5b to compensate Diffraction imaging for the amplitude, phase, and frequency distortions in Figure 5b We now apply our diffraction-focusing approach to the data Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms mentioned previously, which leads to the result presented in Fig- shown in Figure 7, that is, after the suppression of the Love waves ure 5c. Comparing Figure 5a with Figure 5c, we can now see that by different approaches, which gives the corresponding diffraction the dominant Love waves in Figure 5c and 5a match very well. images (Figure 8). In this figure, the horizontal axes show the lateral Next, we subtract Figure 5c from Figure 5a, which gives the result location (in m) as shown in Figure 4, whereas the vertical axes de- shown in Figure 7b. As shown in Figure 7b, two diffraction events note the approximate depth (in m) that we derive using an average with negative moveouts can now be easily recognized. This is due to velocity (V 0 ¼ 150 m/s) from the model in Figure 4. We also use the significant suppression of Love waves in Figure 5a through the this average velocity (V 0 ¼ 150 m/s) in equation 5 to describe the SI + AS procedure (Figure 2a). moveouts of possible diffractions. The range of R (depth of poten- During the procedure of AS, the dominant Love waves dictate the tial diffractors) used for diffraction stacking is from 0.5 to 8 m with parameters of the matching filter; hence, these waves will be most a step of 0.5 m. Such a range covers the area of interest for near- effectively suppressed. The weak diffraction signals also might be surface diffractor detection. To mitigate the near-field effects as affected, but to a lesser extent. As shown in Figure 7b, certain parts much as possible, we mute early arrivals before producing the final of diffractions that overlap with the Love wave are also regarded as diffraction images. Love-wave energy by this algorithm (Figure 2a) and suppressed. To Figure 8a–8c represents the diffraction images from the data after recover this lost diffraction energy and further enhance the ampli- Love-wave elimination by the f-k filtering, SI + AS, and tude of the diffraction events, we then apply SVI to the data as in SI + AS + SVI, respectively. Every second common-source gather Figure 7b, obtaining results as shown in Figure 7c. For comparison, is used. We apply PWS to enhance the coherent summation of weak we also show a reference shot gather (Figure 7d) containing only the signals. From Figure 8a–8c, we can clearly identify two prominent seismic response from the diffractors. This shot gather is obtained anomalies (the red color). We interpret the maximum amplitudes of by taking the difference of synthetic data modeled with and without these anomalies as the centers of the detected diffractors, whose the diffractors. This process removes any arrivals other than the dif- lateral locations agree well with those of the objects embedded fractions. Comparing Figure 7c with Figure 7b, we can see that the in the synthetic model (Figure 4). There are some errors in the esti- diffractions in Figure 7c show more complete moveouts (see the DOI:10.1190/geo2020-0142.1 mated depths of the detected diffractors (approximately 6 m). This blue boxes in Figure 7b and 7c), and their apices are recovered well. can be explained by the fact that a constant veloc- ity (V 0 ¼ 150 m/s) may not be accurate enough to describe the travel path of diffracted waves and a) b) hence may cause errors in the estimated depths. Note that the shapes of the anomalies in Fig- ure 8a–8c do not necessarily indicate the actual shapes of the true objects in Figure 4. Figure 8d– 8f shows diffraction images as in Figure 8a–8c, respectively, but using all modeled common- source gathers. We show this result to investigate the effect of the number of common-source gath- ers on the resolution of the diffraction images. We can see that having half as many sources did not affect the resolution of the obtained dif- fraction images. This is very encouraging for field applications, where due to operational rea- c) d) 10 20 30 40 sons sparser source points would be available. In 0 the following field data example, we have sources every 2 m. Figure 8g–8i shows diffrac- tion images as in Figure 8d–8f, respectively, but 0.05 using linear stack for weak diffraction summa- tions. We can see that the resolution of the anomalies in Figure 8d–8f is higher than in Fig- 0.10 ure 8g–8i. This is because PWS is more efficient for incoherent-noise reduction than linear stack, thus reducing the amount of incoherent noise 0.15 present in the resulting diffraction image. 0.20 FIELD EXAMPLE Figure 7. (a) Result after Love-wave suppression in the data in Figure 5a by f-k filtering, (b) result after adaptive subtraction of Figure 5b from Figure 5a, (c) result after applying Site overview SVI to the result in (b) for diffraction enhancement, and (d) reference shot gather show- ing only diffracted arrivals obtained from subtraction of modeled data using models with In 2017, we acquired seismic data at the Epe- and without diffractors. Niersen barrow alignment, an archaeological site Near-surface diffractor detection WB7 located in the Veluwe, the Netherlands. With our survey, we wanted Near-surface diffractor detection to obtain more information about such monuments (burial mounds) Figure 9a shows a typical raw SH-wave common-source gather Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms using noninvasive geophysical methods, including seismic imaging, with the aim to minimize or even completely eliminate the need for when using the sledgehammer source, in this case at lateral position excavation. To do that, we investigated one burial of 0 m. Figure 9b shows an example shot gather when using the mound from the Epe-Niersen barrow alignment, Location (m) Location (m) Location (m) a) 15 20 25 b) 15 20 25 c) 15 20 25 which is known as Mound 4749 (Bourgeois, 0 0 0 Approx. depth (m) 2012). 2 2 2 Amplitude 4 4 0.02 4 2 200 6 6 6 Seismic data acquisition 8 8 8 0 0 0 We carried out a seismic survey over the top of x10 –6 Mound 4749. We used two kinds of active d) 15 20 25 e) 15 20 25 f) 15 20 25 0 0 0 sources — a sledgehammer and a high-fre- Approx. depth (m) 2 2 2 Amplitude quency S-wave vibrator (Ghose et al., 1996; 4 2 4 0.02 4 200 Brouwer et al., 1997) — to excite seismic waves 6 6 6 8 8 8 recorded by horizontal, 10 Hz single-component 0 0 0 geophones. The horizontal geophones were ori- x10 –6 ented in the crossline direction, whereas the ac- g) 15 20 25 h) 15 20 25 i) 15 20 25 0 0 0 Approx. depth (m) tive sources (hammer and vibrator) were used in 2 2 0.10 2 Amplitude the SH mode, that is, also oriented in the cross- 4 5 4 4 500 6 6 0.05 6 line direction. Under such an acquisition system, 8 8 8 we can generate and record SH-waves. 0 0 0 We used SH-waves because they have several x10 –6 important advantages over compressional waves Figure 8. Comparison of results from diffraction imaging by different approaches. (a– DOI:10.1190/geo2020-0142.1 (P-waves) and/or SV-waves. The first advantage c) Diffraction imaging of the data after Love-wave suppression by f-k filtering, SI + AS, is that they can offer higher resolution of subsur- and SI + AS + SVI, respectively. Every second common-source gather is used. PWS is face structures than P/SV-waves, given the same applied to enhance the coherent summation of weak diffraction signals. (d–f) Similar to frequency content. This is due to the relatively (a–c), but using all of the modeled common-source gathers. (g–i) Similar to (d–f), but the low propagation velocity of SH-waves in soft ordinary linear stack is used to stack the weak diffraction signals. soils and the resulting short wavelength. The sec- ond benefit of SH-waves is that they are directly linked to the small-strain rigidity and are quite sensitive to subtle changes in the subsoil mechanical properties (Ghose, 2003; Ghose and Goudswaard, 2004; Ghose et al., 2013). An- other advantage is that, when SH-waves encoun- ter a diffractive object in the subsurface, the diffracted wavefield will mainly consist of SH- wave diffractions. For P/SV-waves, however, such discontinuity will cause a complex dif- fracted wavefield, which includes P-P, P-SV, SV-P, and SV-SV diffractions. Of all these dif- fraction modes, the SH-wave diffractions have the largest amplitude and most coherent phase characteristic along the traveltime hyperbola, making it advantageous for diffraction imaging (Lellouch and Reshef, 2017; Peterie et al., 2020). Our receiver array consists of 120 geo- phones planted with a 0.25 m interval. During data acquisition, we kept this receiver array fixed and moved only the source at an interval of 2 m. The first source position was at 4 m to the left of the first geophone, whereas the last source posi- tion was at 4 m to the right of the last geophone. Figure 9. (a) A typical preprocessed SH-wave shot gather from the field data acquired At each source position, four recordings were ac- using a sledgehammer source. Due to irregular surface condition and vegetation cover quired and stacked to yield one common-source rendering geophone coupling nonuniform over the array, the field data were rather noisy. The shot gather is dominated by strong Love waves. (c) Result after Love-wave sup- gather. This was done to reduce the source-inco- pression by SI + AS; (e) result after diffraction enhancement by SVI; (b), (d), and (f) are herent noise and increase the S/N of the re- the same as (a), (c), and (e), respectively, but using the S-wave vibrator as the source. corded data. The first 101 traces are displayed. WB8 Liu et al. vibrator as the source at the same source position, after crosscorre- suppress other coherent signals, we then apply SVI to the data lation of the raw vibrograms with the estimated groundforce as shown in Figure 9c and 9d. The results are shown in Figure 9e Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms (Ghose, 2002). The main preprocessing steps that we apply to these and 9f, respectively. As in the synthetic example (Figure 6), Love two types of seismic data are trace editing, top muting, geometrical waves and masked diffraction signals in the preprocessed data (Fig- spreading correction, static correction, trace normalization, and ure 9a and 9b) map to the similar regions in the f-k domain. In such band-pass filtering (5–80 Hz). Due to the soft-soil condition in cases, it is difficult to design an efficient filter to suppress the Love the near surface of our study area, we can see that the raw seismo- waves and preserve the weak diffractions at the same time. Because grams (Figure 9a and 9b) are dominated by distinct, dispersive Love of this fact, we choose not to apply conventional f-k filtering to sup- waves. The strong presence of these Love waves makes it difficult press Love waves in the field data. to identify the diffraction events, if they exist, directly from the raw Figure 10 illustrates the four final diffraction images obtained data. To reveal possible weak diffracted energy, these dominant from the hammer and vibrator data sets (after SI + AS, SI + AS + SVI) Love waves have to be suppressed whereas the diffraction events using our diffraction-focusing approach described in the “Method- should also be preserved. ology” section. The average velocity that we use for diffraction To attenuate the Love waves, we first retrieve the dominant Love- stacking is 90 m/s. We base this velocity value on the results from wave energy from the preprocessed data by SI in a data-driven way, an iterative full-waveform inversion (FWI) algorithm (Tarantola, as explained previously, and then adaptively subtract the retrieved 1984; Virieux and Operto, 2009). Note that the diffraction-focusing result from the raw data (Figure 9a and 9b); that is, we apply SI + AS approach does not require detailed subsurface information, but it which leads to the results shown in Figure 9c and 9d, respectively. assumes a homogeneous value characteristic of the velocity close Comparing the latter with the respective gathers illustrated in Fig- to the surface (V 0 in equation 5). This assumption is only valid ure 9a and 9b, some meaningful diffraction events caused by hetero- at sites with gradual velocity variation. Figure 11 shows the inverted geneities in the subsurface (e.g., the red circles in Figure 9c and 9d) V S model by elastic SH/Love FWI using the hammer data. From the can be identified. To further enhance these diffracted events and velocity model (Figure 11) obtained by FWI, we think that here the near surface (<1.0 m) is quite homogeneous lat- erally and can be, on average, described by a sin- DOI:10.1190/geo2020-0142.1 gle velocity (V 0 = 90 m/s). We stack the diffraction events in the shot domain over diffrac- tor depths from 0.5 to 5.0 m (R in equation 5), which mainly covers the depth range of our in- terest in the shallow subsurface. In Figure 10, the horizontal axes show lateral locations of the receiver array that we deployed over the Mound 4749, whereas the vertical axes indicate the approximate depth (in m) converted from time using an average velocity of 90 m/s. Figure 10a and 10c displays diffraction images obtained from sledgehammer and S-wave vibra- tor data after Love-wave suppression by SI + AS. Figure 10b and 10d represents diffraction images from sledgehammer and S-wave vibrator data after Love-wave suppression by SI + AS and dif- Figure 10. Diffraction imaging result obtained from the field seismic data following our fraction enhancement by SVI. Two clusters of proposed workflow 1: (a) sledgehammer data after Love-wave suppression by SI + AS high-amplitude anomalies (indicated by ellipses) and (b) sledgehammer data after Love-wave suppression by SI + AS and diffraction can be identified at similar positions (approxi- enhancement by SVI; (c) and (d) are the same as (a) and (b), respectively, but using S-wave vibrator data. High-amplitude anomalies (the solid and dashed ellipses) indicate mately 15 and 22 m horizontal distance) in these the potential locations of subsurface diffractor-like objects. four diffractions images. From the V S model (Figure 11) obtained by FWI, we can see that areas at similar positions also show high-velocity (the blue ellipses) contrasts with the background medium. These facts give us more confidence to interpret these distinct anomalies (ellipses in Figure 10) as potential buried objects of archaeological importance. In an earlier separate field work, a magnetic survey was con- ducted at this site (Lambers et al., 2017). This survey detected mag- netic anomalies in the range of 5 ∼ 8 nT, which is nearly the strongest value among their measurements made within the mound. The configuration of the anomalies and their position within the burial mound suggest that these are probably traces of a burial Figure 11. The 2D subsurface V S model obtained by elastic FWI of underneath the mound. Additional corings also confirm the place the sledgehammer data. of a pit at this position — most likely a grave. Earlier excavations Near-surface diffractor detection WB9 of burial mounds in the direct vicinity of this mound have found pretation regarding the presence of these diffractors and their evidence for large stones that were incorporated within the structure approximate lateral locations. The estimated depths of the detected Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms of the graves (Bourgeois et al., 2009). The high-amplitude anoma- diffractors have inherent uncertainty and are connected to the aver- lies detected in our diffraction images (Figure 10) might correspond age velocity V 0 in equation 5. Assuming a constant velocity may to such stones. not be enough to describe the travel path of the diffracted waves, During a seismic survey, when seismic waves encounter such especially at sites with strong lateral or vertical velocity variations. stones, the waves would be diffracted due to the strong impedance This may contribute to errors in the estimated depths. The estimated contrast of the stones with the background medium. In the shot depth information of the detected diffractors could be verified by gathers after SI + AS and SI + AS + SVI, we can identify diffraction other geophysical methods (such as core profiling) once their hori- events (the red circles in Figure 9), with apices located at approx- zontal locations are determined. imately 15 m horizontal distance. However, the diffraction event at After Love-wave suppression and diffraction enhancement, dif- approximately 22 m horizontal distance is not recognizable in the fraction events (the red ellipses in Figure 9) with clear apices could shot gathers. A possible explanation is that the S/N of this diffrac- be identified. These diffraction events are representative of subsur- tion is low. However, after coherent summation of this weak diffrac- face heterogeneities, and they are useful to map the distribution of tion in the following diffraction-stacking procedure, its amplitude buried objects at an archaeological site. However, as mentioned pre- becomes strong enough and can be easily identified in the final dif- viously, a diffraction image can only indicate the approximate lo- fraction images (the dashed ellipses in Figure 10). Comparing the cation of the subsurface objects. To obtain a more complete picture diffraction images from data after SI + AS and SI + AS + SVI, we regarding the subsurface objects of archaeological importance, the can see that the high-amplitude anomalies in the SI + AS + SVI use of 3D seismic imaging and near-surface FWI could be important diffraction image (the ellipses in Figure 10b and 10d) seem to options, which is the direction of our future research. be more easily identifiable than those in the SI + AS diffraction image (the ellipses in Figure 10a and 10c). CONCLUSION DISCUSSION We present a workflow for imaging shallow subsurface objects DOI:10.1190/geo2020-0142.1 When a propagating seismic wave encounters a subsurface object represented by weak diffractions hidden behind dominant SWs (in or a velocity perturbation of size comparable to the wavelength, it our case, Love waves). Our workflow includes three main steps. will be diffracted. For a 2D seismic survey, the imaging of the target The masked diffractions are first revealed after the suppression objects is reliable only when the seismic data are acquired along a of the dominant Love waves. This is done by retrieving the Love line above such targets. If this is not the case, the imaging of dif- waves in a data-driven way by SI and then adaptively subtracting fractors will be negatively affected (spatially smeared or estimated them from the raw data. Second, we enhance the revealed weak at an incorrect location). diffraction signals further through crosscoherence-based SVI. In the modeled example, we saw that f-k filtering could damage Third, we produce a diffraction section by a spatial summation the desired diffraction arrivals. In our case, it indeed resulted in such of the revealed diffraction energy, where no specific velocity for a damage to the right sides of the diffractions (Figure 7a). Still, the the subsurface is needed. We introduce PWS to enhance the coher- energy at and around the apices of the diffractions was preserved, ent summation of weak diffraction signals. Using synthetic data, we and the diffraction stacking gave good results, possibly even better illustrate that our workflow is robust in detecting and imaging weak than the results from our proposed methodology. Figure 7b (and the diffractions. We apply our workflow on seismic data sets acquired filter in Figure 5d) shows that SI + AS had damaged the diffraction in an archaeological site using two different active seismic sources. apices with the result of less strong diffraction stacking. The appli- Our results show two prominent diffraction objects in the subsur- cation of SVI partly compensated for that; thus, the diffraction face at our test location. Our workflow has the potential to be used stacking produced better results. Note that this is also a result of to map the spatial locations of shallow heterogeneities in near-sur- the modeling because the diffractions are relatively strong and thus face seismic surveys, when no detailed subsurface velocity is the SI step retrieves not only the dominant Love waves, but also the available. apices of the diffractions. In the case of the field data, the diffraction events are not clear at all. In such cases, the damaging effect of the f- k filter on diffraction events, which map into the same region in the ACKNOWLEDGMENTS f-k domain as the Love-wave energy, would lead to a lower diffrac- tion-stacking image. Contrary to this, our proposed methodology, The first author is sponsored by the China Scholarship Council making use of a data-driven Love-wave suppression and diffrac- (file no. 201604910851). We thank D. Ngan-Tillard for the help tion-event enhancement through SI + AS + SVI, would not result during the field data acquisition. The authors would like to thank in damaging but in enhancement, and thus would produce better the associate editor A. Lellouch and the four anonymous reviewers diffraction-stacking images. for their constructive comments that have significantly improved The two data sets corresponding to two different seismic sources, the manuscript. sledgehammer and S-wave vibrator, used in this study were inde- pendently acquired and processed. The horizontal locations of dif- fractors are well-constrained by surface seismic methods (Peterie DATA AND MATERIALS AVAILABILITY et al., 2020). The very close lateral locations of the diffractors de- rived from the diffraction focusing applied on SI + AS and Data associated with this research are available and can be ob- SI + AS + SVI processed data sets (Figure 10) validate our inter- tained by contacting the corresponding author. WB10 Liu et al. APPENDIX A a) b) NONSTATIONARY MATCHING FILTER Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms ESTIMATION Matched filtering is a method to measure the time-varying differences between two signals (Fomel, 2009). In this paper, we use it for the AS of SWs. Consider two time series — prediction pðtÞ and observation bðtÞ. A stationary matching filter fðγÞ can be obtained by stretching pðtÞ to different scales to match pðtÞ and dðtÞ. This can be formalized as the following least-squares inversion problem: X Figure B-1. Illustrations of the summation of exp½iϕ1 ðtÞ and min j pðγtÞ fðγÞ − dðtÞj2 ; (A-1) exp½iϕ2 ðtÞ in the complex plane. The red arrows denote these γ two vectors, whereas the blue vector is the addition of these two vectors. (a) When two signals have the same instantaneous phase (ϕ1 ðtÞ ¼ ϕ2 ðtÞ), the amplitude of the phase sum (the blue ar- where denotes the convolution operator and γ is a stretching var- row) will be two and the corresponding value of cðtÞ will be iable. For a particular value of the stretching variable γ, equation A- one. (b) When two signals have significantly different instantaneous 1 gives a single measurement fðγÞ for all of the time coordinates. To phases, the amplitude of the phase sum (the blue arrow) and the corresponding value of cðtÞ will be approximately zero. Modified handle the inherent nonstationary seismic data, a time-varying from Schimmel and Paulssen (1997). matching filter fðγ; tÞ is needed. As equation A-1 shows, this fðγ; tÞ can be obtained by minimizing the following objective func- qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tion: AðtÞ ¼ s2 ðtÞ þ ℋ2 ðsðtÞÞ; (B-3) X min j pðγtÞ fðγ; tÞ − dðtÞj2 : (A-2) DOI:10.1190/geo2020-0142.1 γ and ℋðsðtÞÞ Equation A-2 is an ill-posed problem because it contains more ΦðtÞ ¼ arctan : (B-4) sðtÞ unknown variables than constraints. One remedy is to add addi- tional constraints, that is, regularization, to constrain the variability of the filter coefficients. With shaping regularization (Fomel, 2007), Schimmel and Paulssen (1997) define the following phase stack equation A-2 can be solved as cðtÞ, where no amplitudes of complex traces are explicitly involved: fðγ; tÞ ¼ ½λ2 I þ SðPT P − λ2 IÞ−1 SPT d; (A-3) 1 XN where λ is a scaling coefficient, which is defined as λ ¼ jpðtÞj2 . The cðtÞ ¼ j exp½iϕj ðtÞj; (B-5) N j¼1 term P is the data matrix composed of pðγtÞ. The term S denotes the shaping operator, which is chosen as a triangle smoothing operator in this paper. where N is the number of traces involved in stacking and j is the index of each trace. The amplitude of the phase stack cðtÞ varies between 0 and 1, as schematically illustrated by Figure B-1. As APPENDIX B shown in Figure B-1a, the amplitude of the phase stack cðtÞ equals 1, if the instantaneous phase of the two traces is exactly the same THE PWS METHOD (coherent) at time t. However, when the instantaneous phase of the two traces varies quite significantly, cðtÞ will be approximately 0 PWS is a technology first proposed by Schimmel and Paulssen (Figure B-1b). Thus, we can use the phase stack to weight the sum- (1997) to detect weak but coherent arrivals. The basic idea under- mation of all the traces — gðtÞ is expressed as lying PWS is to suppress components of stacked signals, which do not share the same instantaneous phase. Following the notations from Schimmel and Paulssen (1997), a complex trace SðtÞ is con- ½cðtÞν X N gðtÞ ¼ s ðtÞ: (B-6) structed from a seismic trace sðtÞ and its Hilbert transform ℋ½sðtÞ: N j¼1 j SðtÞ ¼ sðtÞ þ iℋðsðtÞÞ; (B-1) which also can be expressed by amplitude AðtÞ and instantaneous In equation B-6, ½cðtÞν is used to enhance coherent signals, that phase ΦðtÞ: is, signals with similar phase. The exponent ν controls the transition between more coherent and less incoherent signal summations. In SðtÞ ¼ AðtÞ exp½iΦðtÞ; (B-2) our paper, we use ν ¼ 2, as suggested by Schimmel and Paulssen (1997). When ν equals 0, equation B-6 reduces to the conventional where linear stack. Near-surface diffractor detection WB11 REFERENCES Lambers, L., J. W. E. Fassbinder, K. Lambers, and Q. Bourgeois, 2017, The iron-age burial mounds of Epe-Niersen, the Netherlands: Results from magnetometry in the range of ± 1.0 nT: 12th International Conference An, S., and T. Hu, 2016, Suppression of seismic surface waves based on Downloaded 03/16/21 to 145.94.64.230. Redistribution subject to SEG license or copyright; see Terms of Use at https://library.seg.org/page/policies/terms of Archaeological Prospection, 132–134. adaptive weighted super-virtual interferometry: Science China Earth Sci- Landa, E., S. Fomel, and T. Moser, 2006, Path-integral seismic imaging: ences, 59, 2179–2188, doi: 10.1007/s11430-015-5525-7. Geophysical Prospecting, 54, 491–503, doi: 10.1111/j.1365-2478.2006 Balestrini, F., D. Draganov, A. Malehmir, P. Marsden, and R. Ghose, 2020, .00552.x. Improved target illumination at Ludvika Mines of Sweden through seis- Landa, E., and S. Keydar, 1998, Seismic monitoring of diffraction images for mic-interferometric surface-wave suppression: Geophysical Prospecting, detection of local heterogeneities: Geophysics, 63, 1093–1100, doi: 10 68, 200–213, doi: 10.1111/1365-2478.12890. .1190/1.1444387. Bohlen, T., 2002, Parallel 3-D viscoelastic finite difference seismic model- Lellouch, A., and M. Reshef, 2017, Shallow diffraction imaging in an SH- ling: Computers & Geosciences, 28, 887–899, doi: 10.1016/S0098-3004 wave crosshole configurationcrosshole SH diffraction imaging: Geophys- (02)00006-7. ics, 82, no. 1, S9–S18, doi: 10.1190/geo2016-0154.1. Bourgeois, Q., 2012, Monuments on the horizon: The formation of the bar- Liu, J., Q. Bourgeois, R. Ghose, and D. Draganov, 2019, Detection row landscape throughout the 3rd and 2nd millennium BC: Sidestone of near-surface heterogeneities at archaeological sites using seismic Press. diffractions: First Break, 37, 93–97, doi: 10.3997/2214-4609.201902 Bourgeois, Q., L. Amkreutz, and R. Panhuysen, 2009, The Niersen Beaker 407. burial: A renewed study of a century-old excavation: Journal of Archae- Liu, J., D. Draganov, and R. Ghose, 2018, Seismic interferometry facilitating ology in the Low Countries, 1, 83–105. the imaging of shallow shear-wave reflections hidden beneath surface Brouwer, J., R. Ghose, K. Helbig, and V. Nijhof, 1997, The improvement of waves: Near Surface Geophysics, 16, 372–382, doi: 10.3997/1873- geotechnical subsurface models through the application of S-wave reflec- 0604.2018013. tion seismic exploration: 3rd Annual Meeting, EEGS, Extended Ab- Nakata, N., R. Snieder, T. Tsuji, K. Larner, and T. Matsuoka, 2011, Shear stracts, 103–106. wave imaging from traffic noise using seismic interferometry by cross- Dai, W., T. Fei, Y. Luo, and G. T. Schuster, 2011, Super-virtual interfero- coherence: Geophysics, 76, no. 6, SA97–SA106, doi: 10.1190/ metric diffractions as guide stars: 81st Annual International Meeting, geo2010-0188.1. SEG, Expanded Abstracts, 3819–3823, doi: 10.1190/1.3628002. Neidell, N. S., and M. T. Taner, 1971, Semblance and other coherency mea- Dokter, E., D. Köhn, D. Wilken, D. De Nil, and W. Rabbel, 2017, Full wave- sures for multichannel data: Geophysics, 36, 482–497, doi: 10.1190/1 form inversion of SH- and Love-wave data in near-surface prospecting: .1440186. Geophysical Prospecting, 65, 216–236, doi: 10.1111/1365-2478.12549. Peterie, S. L., R. D. Miller, J. Ivanov, and S. D. Sloan, 2020, Shallow tunnel Dong, S., R. He, and G. T. Schuster, 2006, Interferometric prediction and detection using SH-wave diffraction imaging: Geophysics, 85, no. 2, least squares subtraction of surface waves: 76th Annual International EN29–EN37, doi: 10.1190/geo2018-0731.1. Meeting, SEG, Expanded Abstracts, 2783–2786, doi: 10.1190/1 Place, J., D. Draganov, A. Malehmir, C. Juhlin, and C. Wijns, 2019, Cross- .2370102. coherence-based interferometry for the retrieval of first arrivals and sub- Draganov, D., K. Wapenaar, and J. Thorbecke, 2006, Seismic interferom- sequent tomographic imaging of differential weathering: Geophysics, 84, etry: Reconstructing the earth’s reflection response: Geophysics, 71, DOI:10.1190/geo2020-0142.1 no. 4, Q37–Q48, doi: 10.1190/geo2018-0405.1. no. 4, SI61–SI70, doi: 10.1190/1.2209947. Robertsson, J. O., 1996, A numerical free-surface condition for elastic/vis- Fomel, S., 2007, Shaping regularization in geophysical-estimation prob- coelastic finite-difference modeling in the presence of topography: Geo- lems: Geophysics, 72, no. 2, R29–R36, doi: 10.1190/1.2433716. physics, 61, 1921–1934, doi: 10.1190/1.1444107. Fomel, S., 2009, Adaptive multiple subtraction using regularized nonstation- Schimmel, M., and H. Paulssen, 1997, Noise reduction and detection of ary regression: Geophysics, 74, no. 1, V25–V33, doi: 10.1190/1.3043447. weak, coherent signals through phase-weighted stacks: Geophysical Jour- Ghose, R., 2002, High-frequency shear wave reflections from shallow sub- nal International, 130, 497–505, doi: 10.1111/j.1365-246X.1997.tb05664 soil layers using a vibrator source: Sweep cross-correlation versus decon- .x. volution with ground-force derivative: 72nd Annual International Schuster, G. T., 2009, Seismic interferometry: Cambridge University Press. Meeting, SEG, Expanded Abstracts, 1408–1411, doi: 10.1190/1 Shapiro, N. M., M. Campillo, L. Stehly, and M. H. Ritzwoller, 2005, High- .1816924. resolution surface-wave tomography from ambient seismic noise: Sci- Ghose, R., 2003, High-frequency shear-wave reflections to monitor lateral ence, 307, 1615–1618, doi: 10.1126/science.1108339. variations in soil, supplementing downhole geotechnical tests, in J. Sa- Shtivelman, V., and S. Keydar, 2005, Imaging shallow subsurface inhomo- veur, ed., Reclaiming the underground space: Swets & Zeitlinger, geneities by 3D multipath diffraction summation: First Break, 23, 39–42, 827–833. doi: 10.3997/1365-2397.2005001. Ghose, R., J. Brouwer, and V. Nijhof, 1996, A portable S-wave vibrator for Shtivelman, V., S. Keydar, and M. Mikenberg, 2009, Imaging near- high resolution imaging of the shallow subsurface: 58th Annual surface inhomogeneities using weighted multipath summation: International Conference and Exhibition, EAGE, Extended Abstracts, Near Surface Geophysics, 7, 171–177, doi: 10.3997/1873-0604.200 M037, doi: 10.3997/2214-4609.201408721. 9005. Ghose, R., J. Carvalho, and A. Loureiro, 2013, Signature of fault zone de- Slob, E., D. Draganov, and K. Wapenaar, 2007, Interferometric electromag- formation in near-surface soil visible in shear wave seismic reflections: netic Green’s functions representations using propagation invariants: Geophysical Research Letters, 40, 1074–1078, doi: 10.1002/grl.50241. Geophysical Journal International, 169, 60–80, doi: 10.1111/j.1365- Ghose, R., and J. Goudswaard, 2004, Integrating S-wave seismic-reflection 246X.2006.03296.x. data and cone penetration test data using a multiangle multiscale ap- Smith, C., 2014, Encyclopedia of global archaeology: Springer Reference. proach: Geophysics, 69, 440–459, doi: 10.1190/1.1707064. Snieder, R., 2004, Extracting the Green’s function from the correlation of Ghose, R., V. Nijhof, J. Brouwer, Y. Matsubara, Y. Kaida, and T. Takahashi, coda waves: A derivation based on stationary phase: Physical Review E, 1998, Shallow to very shallow, high-resolution reflection seismic using a 69, 046610, doi: 10.1103/PhysRevE.69.046610. portable vibrator system: Geophysics, 63, 1295–1309, doi: 10.1190/1 Taner, M. T., and F. Koehler, 1969, Velocity spectra-digital computer der- .1444431. ivation applications of velocity functions: Geophysics, 34, 859–881, doi: Halliday, D. F., A. Curtis, J. O. Robertsson, and D.-J. van Manen, 2007, 10.1190/1.1440058. Interferometric surface-wave isolation and removal: Geophysics, 72, Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic no. 5, A69–A73, doi: 10.1190/1.2761967. approximation: Geophysics, 49, 1259–1266, doi: 10.1190/1.1441754. Kaslilar, A., U. Harmankaya, K. Wapenaar, and D. Draganov, 2013, Esti- Virieux, J., and S. Operto, 2009, An overview of full-waveform inversion in mating the location of a tunnel using correlation and inversion of Rayleigh exploration geophysics: Geophysics, 74, no. 6, WCC1–WCC26, doi: 10 wave scattering: Geophysical Research Letters, 40, 6084–6088, doi: 10 .1190/1.3238367. .1002/2013GL058462. Wapenaar, K., D. Draganov, and J. O. Robertsson, 2008, Seismic interfer- Klem-Musatov, K., 1994, Theory of seismic diffractions: SEG. ometry: History and present status: SEG. Komatitsch, D., and R. Martin, 2007, An unsplit convolutional perfectly Wapenaar, K., and J. Fokkema, 2006, Green’s function representations for matched layer improved at grazing incidence for the seismic wave equa- seismic interferometry: Geophysics, 71, no. 4, SI33–SI46, doi: 10.1190/1 tion: Geophysics, 72, no. 5, SM155–SM167, doi: 10.1190/1.2757586. .2213955. Konstantaki, L., D. Draganov, R. Ghose, and T. Heimovaara, 2015, Seismic interferometry as a tool for improved imaging of the heterogeneities in the body of a landfill: Journal of Applied Geophysics, 122, 28–39, doi: 10 .1016/j.jappgeo.2015.08.008. Biographies and photographs of the authors are not available.