Three–Dimensional Seismic Diffraction Imaging for Detecting Near-Surface Inhomogeneities

Abstract

One of the problems encountered in a variety of near-surface investigations is detecting and mapping localized inhomogeneities. Typical examples of such inhomogeneous sources are cavities, caves and tunnels. Different methods for detecting shallow subsurface sources utilizing seismic waves diffracted by these sources were proposed by many researchers in the last three decades. Most of these methods suggest that every subsurface point is a possible location of a point diffractor. Imaging of the diffractors is based on a spatial summation of the diffracted wavefield along diffraction time surfaces (defined by source-receiver geometry) in 2D or 3D space. The summation is performed with a fixed velocity value estimated from velocity analysis of the diffraction data. In this study, we present a path integral summation approach, where for every subsurface point the wavefield is stacked together along all possible diffraction time surfaces having a common apex at a given time. The result of the imaging is a 3D volume in which prominent diffraction anomalies appear at spatial locations close to the imaged sources. This path integral summation approach has been successfully tested on synthetic data and further applied at several sites with known subsurface sources.

1. Introduction

The problem of detecting localized near-surface inhomogeneities is encountered in a variety of applications, such as engineering site investigation, environmental studies, archeology and security problems. Typical examples of such sources of inhomogeneities are cavities, caves and tunnels and other natural and man-made local inhomogeneities.

In the past decades, many attempts have been made to solve the problem using various geophysical methods. The methods in use can be roughly subdivided into two groups. The first group includes low-resolution bulk techniques such as microgravity, geoelectric method and seismic surface wave method. A typical drawback of these methods is their low detectability, limiting them to relatively large sources, that is, such sources whose dimension-to-depth ratio exceeds 1/2 [1, 2, 3, 4]. The second group includes the high-resolution ground – penetrating radar (GPR) method. The main drawback of this method is its high sensitivity to electrical conductivity of near-surface material which makes it limited to very shallow depths usually not exceeding a number of meters. Furthermore, our experience shows that the general problem of all the existing techniques is a high degree of ambiguity of data interpretation and, as a result, a high percentage of errors. One of the methods that may prove very effective for this purpose is active seismic, which is widely employed by seismic exploration industry in its extensive and detailed studies of the subsurface. Adaptation of these methods for the purpose of searching for local underground heterogeneities is not an easy task. The primary aim of seismic exploration is to study geological sources of large lateral extent, such as boundaries of geologic formations, oil and gas deposits, etc. To this end, seismic exploration (“reflection seismic”) utilizes seismic waves emitted from man-made (active) seismic sources and reflected from the boundaries of interest back to the surface. The methodology, the data acquisition and processing techniques employed by the reflection seismic are largely based on the model of the subsurface consisting of a number of layers divided by smooth curvilinear boundaries with gradual changes of parameters within layers. Virtually, the same model of the subsurface is behind the refraction seismic method, another technology used primarily to study the shallow part of the subsurface (first tens of meters). Clearly, such methods are not directly applicable to a search for local subsurface inhomogeneities of small lateral extent (10–100 m). The amount of energy reflected by such small sources would be very small, and it would be indistinguishable from other geological features of the subsurface. In order to be effective, the search for local inhomogeneities has been employed to image the features of seismic wavefield, which is primarily associated with the presence of local heterogeneity in the subsurface. Such features are known as scattered (diffracted) waves. In a standard seismic prospecting, these waves are usually treated as noise. But for locating local sources these waves provide the most valuable source of information [5, 6, 7, 8, 9, 10, 11, 12].

2. Methods for locating subsurface heterogeneities

2.1. Wavefield

Let us look at the resulting wavefield obtained from a typical underground heterogeneity in a typical layered geological environment. The wavefield could be divided into a number of different kinds of waves such as primary and multiple reflected waves, refracted waves and scattered/diffracted waves. The reflected waves originate from the layered boundaries, the refracted waves from boundaries of high contrast while the scattered/diffracted waves are generated by local heterogeneities. The reflected and refracted wavefields would only slightly be affected by the presence of the local heterogeneities. However, since the scattered/diffracted waves are generated by local heterogeneities, they may be used for the imaging of those heterogeneities.

2.2. Diffraction imaging

Imaging of the local heterogeneities in 3D space can be performed by a spatial summation of scattered/diffracted waves along the travel-time surfaces of those waves [9].

The summation can be implemented using one of the two approaches. In the first step (conventional approach), target waves are stacked along the time surfaces defined by some “optimal” velocity values. These values are usually estimated by a time-consuming procedure of velocity analysis involving visual examination and interpretation of intermediate seismic images. Recently, an alternative method, a more formal, path integral summation approach, has been proposed [7, 10, 11, 12]. Path integral summation is performed by stacking the target waves along all possible time surfaces having a common apex at the given point. This approach does not require any explicit information on velocities since the required multipath summation is performed for all possible velocity values within a wide specified range. Application of the path integral summation to synthetic and real data shows that it can provide reliable images of subsurface objects close to their true spatial locations [9]. However, the spatial resolution of the images obtained by the path integral summation is usually inferior to those obtained using the fixed velocity.

In the subsequent sections, we give a brief description of the main principles of the method and illustrate its application by two synthetic examples. Then, we present four case studies demonstrating how the method can be implemented for detecting various subsurface objects located in different geological environments.

2.3. Theoretical background

The method suggests that every subsurface point is a possible location of a point diffractor. Imaging of the diffractors is based on a spatial summation of the diffracted wavefield along diffraction time surfaces (defined by source-receiver geometry) in 3D space. The result of the imaging is a 3D volume in which prominent diffraction anomalies appear at spatial locations close to the imaged objects. The proposed method proved to be a robust high-resolution technique capable of providing reliable results under various geological conditions in a wide range of depths. The method seems to be especially suitable for detecting relatively small objects.

The proposed method is based on a spatial summation of diffracted waves in 3D space. The principles of the method can be elucidated using a simple scheme represented in Figure 1. Assume that a source of seismic energy is located at the point (X_S, Y_S, Z_S) and that there is a diffractor (i.e., a void) at the point (X_D, Y_D, Z_D) within a homogeneous medium with a velocity V. When seismic waves generated by the source reach the diffractor, it starts to operate as a secondary source generating diffracted waves. At a certain time the diffracted wave arrives at a receiver located at the point (X_R, Y_R, Z_R). For the given source-receiver pair and a fixed diffractor position, the arrival time can be represented as a sum of two terms: source-to-diffractor time and diffractor-to-receiver time. In a homogeneous medium, this sum is totally defined by spatial coordinates of the source, receiver and diffractor and the velocity in the medium. Assuming that the source and receiver are located at the surface Z = 0 (i.e., Z_S = Z_R = 0), the arrival time of the diffracted wave can be defined as follows [9]:

$$T_{\mathit{SR}}=\sqrt{{\left({\xi }_{\mathit{SD}}/V\right)}^2+{T_D}^2}+\sqrt{{\left({\xi }_{\mathit{DR}}/V\right)}^2+{T_D}^2},$$E1

where ${\xi }_{\mathit{SD}}=\sqrt{{\left(X_D-X_S\right)}^2+{\left(Y_D-Y_S\right)}^2}$ is the horizontal source-to-diffractor distance; ${\xi }_{\mathit{DR}}=\sqrt{{\left(X_R-X_D\right)}^2+{\left(Y_R-Y_D\right)}^2}$ is the horizontal diffractor-to-receiver distance; $T_D=Z_D/V$ is vertical time above the diffractor; and V is the velocity in the medium. The velocity defines the shape of the diffraction time surface, whereas time T_D defines its apex. For a spatial array of receivers, equation (1) defines a diffraction time surface with its apex at time 2T_D. If the velocity in the medium is known, the diffractor can be imaged by summing seismic energy arriving at different receivers along the time surface defined by equation (1). The principles of the imaging can be formulated as follows:

• A targeted source (a void) acts as a secondary source of waves (diffracted waves).

• Every point in the subsurface is considered as a possible diffractor.

• Seismic data acquisition is performed by a 3D survey involving a large amount of source and receiver points.

• The imaging is based on a spatial summation of diffracted waves arriving at all the receivers from all the sources.

• The summation is performed using Eq. (1) with a fixed velocity value estimated from velocity analysis of the diffraction data.

• Alternatively, it can be done using a path integral summation approach, where for every subsurface point, the wavefield is stacked together along all possible diffraction time surfaces having a common apex at a given time [10].

• The result is a 3D image of the subsurface (a 3D volume) where the diffractor (the void) appears as a distinct anomaly.

  Figure 2a is a schematic illustration of a 3D volume containing a diffraction anomaly obtained by the imaging. Note that the vertical axis is two-way time rather than depth, since the imaging is performed in the time domain. Subsequently, the results of the imaging will be presented in two forms. The first one is obtained by vertical slicing of the 3D volume parallel to the X axis, as illustrated by Figure 2a. In this manner, we obtain a set of vertical diffraction sections computed for a number of Y values and represented in the X-T plane (we call them D-sections). The second form is obtained as follows (Figure 2b): for every output trace a(t) resulting from the spatial summation of the corresponding input traces, we compute its energy by summing squared values of the trace amplitudes and put the resulting sum at the corresponding trace location (X_D,Y_D) at the surface. In such a manner we obtain a map of trace energy (we will call it D-map) on which the diffraction anomaly is represented by its horizontal projection onto the X-Y plane, as illustrated by Figure 2b.

  As mentioned above, the imaging can be performed using either a fixed velocity or path integral approach. The path integral imaging was described and applied in our previous work [9]. In this work we applied the fixed velocity imaging (which usually provides a better resolution).

Figure 1.

Schematic illustration of basic diffraction geometry.

Figure 2.

(a) Schematic illustration of a 3D volume and D-sections in X-T plane and (b) schematic illustration of D-map in X-Y plane.

3. Synthetic examples

3.1. Point and linear diffractors

The synthetic examples presented below illustrate application of the method for imaging two kinds of objects: a point diffractor (which can be considered as approximation of a small quasi-isometric source, such as a small karstic cavity) and a linear diffractor (roughly approximating an elongated source, such as a tube or a tunnel).

3.1.1. Point diffractor

The model in the first example includes a point diffractor located at a point with the coordinates X = 40 m, Y = 10 m and Z = 10 m. The medium is assumed to be homogeneous with the velocity of 1000 m/s. The modeling was performed using a simulated acquisition system represented in Figure 3. The system included 2 receiver lines with 48 receivers on each line and 5 source lines with 48 source points on each line. The source and receiver spacing along the lines (i.e., in the X direction) was 2.5 m. The distance between the source lines was 5 m and between the receiver lines 20 m, so that the area covered by the system was about 120 m × 20 m. Using this system, synthetic seismograms of diffracted waves were computed for each source position, given the total of 240 shot records with 48 traces in each record. Then, a spatially uncorrelated random noise with the average s/n ratio of 0.03 was added to the seismograms. Figure 4 represents an example of one of the records. Due to the low s/n ratio, no traces of diffracted events were identified on the record.

Figure 3.

Source-receiver geometry used for a point diffractor modeling (plan view).

Figure 4.

Point diffractor model: an example of a synthetic record.

Figure 5a and b represents the results of the point diffractor imaging with the velocity of 1000 m/s. Figure 5a shows 5 D-sections computed along the lines defined by the Y values of 0. 5, 10, 15, and 20 m. The sections display a well-defined anomaly whose maximum is located at X = 40 m on the section with Y = 10 m. The amplitude of the anomaly is rapidly attenuated away from its maximum in both the X and Y directions. This can be seen even more clearly on the D-map (Figure 5b): the figure shows that the diffraction anomaly is very well localized on the X-Y plane around the imaged point diffractor.

Figure 5.

(a) Point diffractor model: D-sections and (b) point diffractor model: D-map.

3.1.2. Linear diffractor

The model in the second synthetic example contains a linear diffractor (which can be considered as a set of densely spaced point diffractors). The diffractor crosses the X axis at the distance of 40 m from the origin at the angle of 30° (Figure 6). The model parameters and the acquisition system were the same as in the previous example.

Figure 6.

Source-receiver geometry used for linear diffractor imaging (plan view).

The D-sections obtained for this model with the velocity of 1000 m/s (Figure 7a) display a well-defined anomaly whose location on the sections varies along the Y axis in accordance with the inclination of the imaged linear diffractor. This is further illustrated by the D-map (Figure 7b). The map shows a prolongated quasi-linear anomaly whose location and inclination on the X-Y plane correspond well with the imaged linear diffractor.

Figure 7.

(a) Linear diffractor model: D-sections and (b) linear diffractor model: D-map.

4. Case studies

4.1. Dead Sea-Israel

The first case is from a sinkhole site located in the Dead Sea area of Israel. A well drilled at the site revealed a water-filled cavity located in a salt layer within the depth of 22–28 m. The cavity was a target of the diffraction survey carried out at the site. The acquisition system used in the survey was similar to that described in the earlier synthetic examples. Imaging was performed with the velocity of 1200 m/s estimated from the velocity analysis of the data. The resulting D-sections (Figure 8a) show a prominent anomaly located in the area close to the well location. The depth of the anomaly was estimated to be about 25 m which is quite close to the depth of the cavity as revealed in the well. The D-map (Figure 8b) displays an isometric anomaly implying that the cavity might have a quasi-spherical shape.

Figure 8.

(a) Case 1: D-sections and (b) case 1: D-map.

4.2. Limestone quarry

The second case is from a limestone quarry where a number of karstic cavities were revealed by drilling shallow boreholes. The depth of the cavities varied between 8 and 23 m. At several sites in the quarry, diffraction surveys were carried out with the aim of detecting the cavities from the surface. We present an example from one of the surveys. The acquisition system used in this case was half as long as the system specified earlier (see Figure 3), that is, it included only 24 receivers and 24 source points per line, so that the area covered by the survey was about 60 m × 20 m. The well located near the point (X = 35 m, Y = 10 m) revealed a cavity at the depth of 8 m. Imaging was performed with the velocity of 1500 m/s. The D-sections (Figure 9a) show a well-defined anomaly with its maximum located very close to the well location. The amplitude of the anomaly decreases rapidly away from the maximum. The depth of the anomaly was estimated to be about 7 m. Thus, the estimated spatial location of the target object is in good agreement with the location of the cavity revealed by drilling. This conclusion is further confirmed by the D-map (Figure 9b) that shows an isometric anomaly located in close vicinity around the well.

Figure 9.

(a) Case 2: D-sections and (b) case 2: D-map.

4.3. Chalky hills

The third case refers to an ancient subsurface aqueduct which has been recently discovered in Chalky hills, some 23 km east of the Mediterranean shore of Israel (Figure 10a). The aqueduct was part of a water supply system operating during the Roman and Byzantine periods. In the excavated part of the aqueduct, it is about 1.8 m high and 0.8 m wide and its depth from the surface is from 8 to 10 m. A diffraction survey was carried out at the site located above an extension of the aqueduct beyond its excavated part, at the distance of about 200 m from the nearest excavated shaft. The acquisition system was similar to that specified in Figure 3. The aqueduct was supposed to cross the system at a distance of 80–90 m from the beginning of the lines. Imaging was performed with the velocity of 1200 m/s. The D-sections (Figure 10b) show a prominent anomaly located in the area which is close to the estimated extension of the aqueduct. The D-map (Figure 10c) displays an elongated anomaly whose inclination on the X-Y plane is in good agreement with the general direction of the aqueduct.

Figure 10.

(a) Ancient subsurface aqueduct, (b) case 3: D-sections and (c) case 3: D-map.

4.4. Modern operating aqueduct

The fourth case is from a site located above a modern operating aqueduct which supplies water to a large part of Israel. The diameter of the aqueduct is 3 m. At the investigated site it is located at the depth of 27 m within the limestone. In this case our target object was quite a small one, since the ratio of its dimension (diameter) to the depth to its center is about 1/10, making it quite challenging.

The diffraction survey carried out at the site used the acquisition system including 2 receiver lines with 24 receivers on each line and 4 source lines with 24 source points on each line. The source spacing was 5 m in both the X and Y directions; the receiver spacing was 5 m in the X direction and 15 m in Y direction. The survey crossed the aqueduct at approximately the right angle at a distance of about 40 m from the beginning of the lines. Imaging was performed with the velocity of 1100 m/s. The D-sections obtained for the data (Figure 11a) show a distinct anomaly at X = 45 m. The estimated depth of the anomaly is about 24 m. Thus, the estimated spatial location of the target object corresponds fairly well with the actual location of the aqueduct. This is confirmed by the D-map (Figure 11b) which displays a quasi-linear anomaly whose location on the X-Y plane is close to the aqueduct location.

Figure 11.

(a) Case 4: D-sections and (b) case 4: D-map.

5. Summary

The main points of our work can be summarized as follows: a method for detecting and mapping localized near-surface inhomogeneites was presented.

The method is based on 3D subsurface imaging using spatial summation of diffracted waves along their time surfaces.

The method was successfully tested on synthetic data with a very low s/n ratio.

Application of the method to real data demonstrates that it can provide a reliable detection of relatively small subsurface sources located under different geographical conditions at a wide range of depths.

Acknowledgments

The authors are grateful to the Geophysical Institute of Israel for granting permission to publish the datasets.