1. Introduction
In a 2D seismic survey all sources and receivers are usually located on a straight line. Each shot is simultaneously recorded by an entire group of equally spaced receivers. The receivers are located either on one side of the source only, or on both sides of it. Each receiver records a trace that constitutes a time series of ground motion or ground acceleration amplitudes. The final imaging result is a cross section of the subsurface structure, which extends vertically below the acquisition line. Since data is recorded shot-by-shot, a dataset initially consists of a sequence of so-called shot gathers. However, it is common practice in seismic imaging to reorganize the traces into so-called common midpoint (CMP) gathers, each composed of records from different experiments which have source and receiver centered around the same midpoint, but a variable source–receiver distance, called offset.
In the very basic scenario of having a planar measurement surface and a parallel planar interface separating two homogeneous subsurface layers, we would record the same trace for every midpoint location. The reflection traveltime would depend only on the offset. If we represent all measured traces in a midpoint, offset, traveltime space, all reflection amplitudes would lie on a plane with hyperbolic curvature in offset direction and zero curvature in midpoint direction. The shortest reflection traveltimes would be found in the smallest-offset plane. If we had a reflector dip relative to the measurement surface, we would still find the same traveltime surface but tilted with respect to the midpoint axis by an angle proportional to the reflector dip. If we go one step further and assume a circularly shaped reflector segment with a finite curvature and a dip relative to the surface, we obtain a traveltime surface tilted with respect to the midpoint axis by an angle corresponding to the reflector dip, and with a curvature in midpoint direction proportional to the reflector curvature. For an infinite curvature of the reflector, i.e., for a diffraction point, we would obtain a bell-shaped traveltime surface, which is hyperbolic both in offset and midpoint direction.
Let us now assume that our records not only contain the reflected signal energy of the source pulse but also noise energy from ambient sources, like nearby street traffic or construction activities. These uncorrelated noises can be reduced by stacking amplitudes of traces with (a) different offsets and (b) different reflection points, located on the same reflector within a certain neighborhood, which translates in practice into a range of midpoints. To facilitate structural interpretation, the stack amplitude is placed at the extrapolated traveltime that would be registered by a notional experiment with coincident source and receiver positions. The corresponding Zero-Offset (ZO) ray hits the reflector normally. It ascends on the same trajectory towards the receiver on which it descended from the source down to the reflector. If the projected midpoint stacking aperture on the reflector, which is an area around the Normal-Incidence Point (NIP) of the ZO ray, does not exceed the size of the first interface Fresnel zone of the reflected wavefield, stacking can be expected to increase the Signal-to-Noise (S/N) ratio without harming resolution [
1]. Obviously, only a smaller part of the amplitudes stemming from a laterally extensive continuous reflector can be stacked into a single ZO traveltime. This is different for a diffraction where in principle the entire traveltime surface can be collapsed to a point.
In real applications, the structure of the subsurface is unknown. Therefore, scanning for reflection and diffraction events, so-called velocity analysis, must precede each stacking process. Velocity analysis and stacking need to be conducted for all traveltimes before a simulated ZO trace is obtained, and for all midpoints before the final ZO section is completed. While stacking creates a simulated ZO section of high S/N ratio, velocity analysis delivers geometrical information regarding the reflecting interface and the wave propagation velocities within the strata set above it, blended into one or more parameters. The simplistic two-layer scenarios discussed previously in this section illustrated that for 2D seismic reflection events are surfaces in a 3D multi-coverage dataspace. This makes clear that both velocity analysis and stacking should be performed along these surfaces in midpoint and offset direction.
In the pre-digital age of seismic prospecting, it was impracticable, if not unthinkable, to search for spatial reflection events in 3D pre-stack data volumes. After the introduction of multi-coverage data, Common-Midpoint (CMP) stacking [
2] was considered a breakthrough and maintained large popularity for a long time. The method limits both stack and velocity analysis to single CMP gathers, where one parameter, the Normal-Moveout (NMO) velocity, is enough to parameterize the moveout hyperbola along which the recorded amplitudes are shifted to ZO traveltime and summed up. The difference between the traveltime of a real finite-offset ray and the traveltime of the notional ZO ray is called normal-moveout.
For a planar surface and a parallel reflector separating two homogeneous layers, the midpoint of the ZO trace coincides with the lateral location of the reflection point, while the reflection traveltime multiplied by the NMO velocity is equal to the length of the ray path or twice the depth of the reflector. The increase in traveltime with offset is exactly hyperbolic and is described by the well-known NMO formula [
3], that reads
with
h being the half-offset, i.e., half of the offset,
t0 being the ZO traveltime and
vNMO being the NMO velocity. If we correct the traveltimes respectively, we obtain a flat NMO corrected event that we can stack in offset direction. Using an extra- and interpolated continuous velocity model, this process is applied for all traveltimes of the stacked ZO trace to be constructed.
Figure 1 shows the different stages of NMO velocity analysis and CMP stacking, using an example from ground penetrating radar. In
Figure 1a, a CMP gather with six distinct reflection events is displayed, overlain by those hyperbolas that provide the closest fit to these reflections. The latter were determined by coherence analysis, whereas semblance [
4] was used to measure the coherence of the data along various candidate hyperbolas. During this process, the semblance spectrum depicted in
Figure 1b, was obtained, where subsequently the best fitting velocities were picked for the six reflections and extrapolated and interpolated to a continuous NMO velocity model.
Figure 1c, shows the CMP gather after NMO correction. Stacking this gather in offset direction and placing the stack results at the traveltimes which the velocity model predicts for zero-offset produces the stacked ZO trace depicted in
Figure 1d. To better cope with dipping reflectors, Dip-Moveout Correction (DMO) (see [
5] and references given therein) was later added to CMP stacking.
NMO correction with subsequent CMP stacking was often reviewed, updated, and extended (see, e.g., [
7,
8,
9,
10]). The Common Reflection Surface (CRS) stack method (see [
11,
12,
13] and references given therein) constitutes one of the latest and probably the most successful generalizations of this technique. The CRS method extends CMP stacking towards the off-CMP dimension by adding a relative midpoint coordinate to the traveltime approximation. It provides a second-order approximation of the two-way traveltime of rays, reflected on a common reflection surface (3D seismic) or common reflection segment (2D seismic), respectively. The latter is depicted in
Figure 2a as the red area around the normal incidence point R of the central ZO ray (straight blue line), which is defined by its emergence location
x0 and traveltime
t0. The same situation, but with inverted time axis, can be found in
Figure 2b, where a collection of neighboring CMP gathers from near-surface seismic data is depicted side-by-side, in order to illustrate the continuation of reflection events in both midpoint and offset direction.
For 2D acquisition on a plane surface without topography overlying an inhomogeneous but isotropic medium, according to [
15], the hyperbolic CRS stacking operator in midpoint displacement, ∆
xm, and half-offset,
h, coordinates can be written as
with
. The near-surface velocity in the vicinity of
,
, is assumed to be constant and a priori known. It plays the role of a constant of proportionality, allowing to parametrize a hyperbolic traveltime surface by three independent kinematic properties of the wavefield measured in
: α, the emergence angle of the central ZO ray and
and
, the wavefront radii of two notional eigenwaves known as normal incident point (NIP) wave and normal (N) wave, as illustrated in
Figure 3.
These so-called CRS parameters represent the generalized counterpart to the NMO velocity which can be expressed as
Looking closer at how CRS stacking is actually implemented, we observe a basic difference to the conventional NMO/DMO correction plus CMP stack. Following the paradigm of data-driven imaging, an automated search for the optimum stacking parameter triple
) is performed for every sample of the output section instead of stacking with NMO velocities extrapolated and interpolated from a limited number of visually detectable reflection events, as portrayed in
Figure 1. In other words, for every ZO sample
, those values for
,
, and
are searched for which parameterize the CRS traveltime surface, given by Equation (2), in such a way that it fits best to a measured reflection in the data. If
does not lie on a reflection event, the attribute triple with the maximum coherence will nevertheless achieve only a relatively low coherence value and stacking will not result in a constructive summation of amplitudes.
A positive side-effect of treating every sample of the reflection event independently is that the so-called NMO-stretch (see
Figure 1c, event four), which deteriorates the resulting stack section, is avoided [
17]. This is because, unlike the conventional NMO velocity model shown in
Figure 1b, the NMO velocity calculated from the CRS stacking parameters is generally neither constant for the entire duration of a reflection event nor does it increase due to velocity extra- or interpolation. Instead, it generally decreases with time from the top to the bottom of the event [
16]. Why do we need NMO velocities that decrease over the wavelet to avoid the NMO stretch? Our hypothesis is that this results from the slightly unphysical way traveltimes are measured, according to which the last sample of the wavelet was not only recorded later but also has a larger traveltime than the first sample and thus has apparently traveled with a lower speed. Consequently, the NMO velocity must decrease towards the tail end of reflection events, as the apparent traveltime increases while the depth of the reflecting interface remains constant. This can also be seen from Equation (3), assuming a homogeneous layer with velocity
over a horizontal interface. In this case, the ZO ray emerges vertically, i.e., α = 0. The NMO velocity is equal to
if, and only if,
is the exact time between the firing of the shot and the arrival of the head of the wavelet at the surface; only then
holds. For the entire tail of the wavelet, which arrives later and corresponds to a larger
, the NMO-velocity will be smaller than
.
Fitting a surface in a three-dimensional data space can hardly be done by an interpreter and decisively not for every single image point of the ZO section. That is why such an entirely data-driven procedure has not been considered before the advent of data digitization and powerful computers. These technical breakthroughs were aided by scientific breakthroughs in the field of dynamic ray theory (see, [
18] and references given therein) which finally allowed the substitution of the generic stacking parameters
A,
B,
C, resulting from a second-order Taylor expansion in offset and midpoint, with physically interpretable quantities. The importance of this cannot be overstated since, for physically interpretable quantities, search ranges can be defined in a rational way and correlated properties, like the size of the projected Fresnel zone, the Common Reflection Point (CRP) trajectory, or the time-migrated image point location can be calculated (see, [
19]). Having three instead of one physically meaningful stacking parameter also enabled the NIP wave tomography [
20], the extension for rough surface topography [
21], CRS-stack-based residual static corrections [
22], pre-stack seismic data enhancement [
23], and finally the CRP time migration without velocity model, presented in [
24], which might be considered the capstone on top of the numerous applications that emerged from the CRS method. This paper is limited to the 2D case. For an overview of the 3D case, we like to refer the reader to [
18,
25,
26] and references given therein.
3. Shallow P-Wave Data Example
In the following we will conduct a detailed comparison between the three search strategies presented in the previous section. For this purpose, we reprocessed seismic line 2 of a P-wave data set acquired, processed and interpreted by [
46]. CRS stack results for the less noisy line 1, obtained using the pragmatic search strategy, are discussed in [
47]. The acquisition, which took place in the South-East of Sardinia, the second largest of Italy’s islands, can be categorized as a hydrogeological survey. In Sardinia, water is a precious resource, frequently in short supply during the long, hot, and dry summer, creating distribution conflicts between agriculture and tourism, as well as environmental problems, such as the salinization of coastal aquifers. The two lines, having a length of 1.1 km each, cross the coastal plain of the Flumendosa River, on the southeastern coast of Sardinia, Italy, roughly parallel to the coastline of the Tyrrhenian Sea. The survey was aiming both at imaging the Paleozoic bedrock topography, and at obtaining detailed structural and stratigraphic information on the sequence of largely fluvial sediments extending from the surface down to bedrock. The study aimed to gain a better understanding of the geological and hydrogeological factors influencing a highly productive aquifer located in thick Quaternary deposits and subject to extensive saltwater intrusion. For a detailed description of the geological and hydrogeological features of the target area we like to refer to [
46].
The 89 shot-gathers were recorded using a standard common-midpoint (CMP) roll-along technique in an end-on configuration with 48 active 50-Hz geophones. At every shot point, 0.25-kg of explosives were fired. The sources were placed at a depth of 2 m, which lay in general below the water table. Geophone spacing of 5 m and source spacing of 10 m provided twelvefold CMP coverage with a CMP spacing of 2.5 m. The recorded traces had a sample interval of 0.5 ms and a record length of 1024 ms. Expecting a maximum reflector depth of 200–300 m a maximum source receiver offset of 245 m was chosen, which was right for line 2, where the maximum depth of the basement turned out to be ca. 200 m.
Three representative shot gathers, for the sake of simplicity called records in the following, are shown in
Figure 5. In all records, the airwave (event
e) is clearly visible. An up-dip reflector at the beginning of the line is presumably related to event
f, with clear reversed moveout from record 3 and visible up to record 11. On the other end of the line where record 72 was selected as representative sample record, event
g can be interpreted as refracted arrivals, stemming from a very shallow refractor, corresponding to the interface between the bedrock and the overlying sediments.
Line 2 terminates on the southern bank of the Flumendosa River, 140 m from a Paleozoic outcrop. We chose line 2, because, of the two lines, it has the lower S/N ratio and is therefore better suited to highlight the benefits of a spatial velocity analysis. The annotated stack section produced by [
46], and depicted in
Figure 6, shows two dominant reflections. The first, event α, is the deepest coherent reflection in the section and was interpreted as the top of the Paleozoic bedrock, in accordance with the outcrop geology beyond the end of the line. It has an apparent dip to the north-northeast from about 290–165 ms two-way-traveltime (twtt) between CMP locations 10 and 60, before it gently shallows with a dip angle of about 10°, reaching 125 ms twtt at CMP 160. Even though less continuous and with lower amplitudes, event α can further be identified, from CMP 160 to the end of the line, with its change in reflection character probably caused by a change in reflectivity of the overlying sedimentary sequence. After CMP 160 event α abruptly starts descending until CMP 200, forming a trough with flat bottom between CDP 200–250, before it gently up-dips from 160 ms to 130 ms twtt at CMP 300. From CMP 300 to the end of the line, it continues to shallow, reaching 66 ms twtt, or a depth of approximately 52 m. The outcrop geology, and the velocity model depicted in Figure 13b, obtained by [
46] from combining reflection-derived interval velocities and those derived from first-arrival information using refraction tomography, support the interpretation of event α as the top of the bedrock. The second of the two dominant reflection events is called event γ and is interpreted as the boundary between Holocene and Pleistocene alluvia. It starts at the beginning of the line and extends at least until CMP 340, where event α pinches out. It is nearly horizontal, with onset at 50–70 ms twtt. Event γ is overlain by seismic unit A that extends to the surface and consists of Holocene alluvial deposits composed of conglomerates, sands, and clays, according to samples obtained by drilling. Not all events in unit A can be interpreted as genuine reflections, because this very shallow part of the recordings is contaminated by various types of source-generated, high-amplitude coherent noise such as surface and refracted waves. Finally, at the left side of the section, unit B can be found between event γ and event α. Ref. [
46] describes this unit as being composed of Pleistocene alluvial deposits.
For creating the series of CRS stacking results shown in
Figure 7,
Figure 8,
Figure 9 and
Figure 10, we used the following processing parameters. We chose a laterally constant near-surface velocity of 1600 m/s, consistent with the P-wave velocity of the water saturated sediments close to the ground surface. The total width of the coherence band was set to 33 samples, i.e., 0.0165 s. The offset aperture was 30 m for traveltimes up to 0.03 s and then increased linearly until the full offset range of 245 m was reached at 0.1 s. Roughly, this means that offsets do not become larger than two-times the reflector depth. The midpoint aperture increased linearly from the smallest to the largest traveltimes starting with 5 m halfwidth, i.e., 5 CDPs, and reaching finally 15 m halfwidth, i.e., 13 CDPs, for the maximum traveltime of 0.5 s. This means for the target area above the bedrock, i.e., 0 to 0.25 s, we have a midpoint aperture of 5 to 9 CDPs.
The CRS results presented in the following section are created in a fully data-driven way without the need of setting more than the above processing parameter values. We begin our comparison of the stack results with the automatic CMP stack result depicted in
Figure 7. This section is the basis for the two subsequent search steps of the cascaded search strategy. The result of the complete three-step strategy is depicted in
Figure 8.
Figure 9 shows the stack result obtained from CRS stack using cascaded search followed by a local three-parameter optimization.
Figure 10 shows the stack result obtained from CRS stack using cascaded search plus three iterations of event-consistent smoothing, each followed by a local three-parameter optimization. Now, coming to the simultaneous multi-parameter search strategies using spatial operators, we see in
Figure 11 the result of the CRS stack using the hybrid diffraction/reflection search algorithm. Finally, all these results can be compared to the result of the CRS stack using simultaneous three-parameter search displayed in
Figure 12.
From
Figure 7,
Figure 8,
Figure 9 and
Figure 10, it is obvious how the stack quality slowly improves with every step of the cascaded imaging strategy, including event-consistent smoothing. The comparison between the stack section displayed in
Figure 10 and those displayed in
Figure 11 and
Figure 12 shows clearly the superior reflection imaging capability of both global simultaneous search strategies. Not only reflections, but also diffractions, such as those at the left border (house wall) and in the center of the stack section (shallow object at the right end of the line), are more clearly imaged after velocity analysis with spatial operators.
Regarding the velocity analysis, the authors of [
46] state: ”As is often the case, detailed velocity analysis proved to be the most important and time-consuming processing step. For both Lines 1 and 2, the initial velocity models were developed through integrated analysis of constant velocity gathers, constant velocity stacks, and semblance plots. Hand-picked direct arrivals supplied additional information on near-surface velocities.” This laborious and often quite subjective process must be kept in mind, when confronting the output shown in
Figure 13 to the data-driven CRS results depicted in
Figure 7,
Figure 8,
Figure 9,
Figure 10,
Figure 11,
Figure 12,
Figure 14,
Figure 15,
Figure 16 and
Figure 17. In
Figure 14, the NMO velocities obtained from the cascaded search followed by local optimization are displayed, while in
Figure 15 we display the NMO velocities obtained from the cascaded search plus three iterations of event-consistent smoothing, each followed by a local optimization. In
Figure 16, the NMO velocities from the hybrid diffraction/reflection search are shown and in
Figure 17 the NMO velocities from the full simultaneous three-parameter search. In all the CRS stack examples, the NMO velocities were calculated according to Equation (3) from the estimated CRS parameters.
Comparing the NMO velocities in
Figure 14 and
Figure 15 and those in
Figure 16 and
Figure 17 demonstrates the superior accuracy of the global simultaneous search strategies in estimating NMO velocities and the underlying stacking parameters α and R
NIP. In addition, diffractions, such as those at the left border (house wall) and in the center of the stack section (shallow object at the right end of the line), are much more clearly imaged.
4. Discussion
If we compare
Figure 7 and
Figure 13b, it is obvious that the automatic CMP stack used as the first step of the pragmatic search strategy is far less successful than the conventional CMP stack methodology. For the making of this judgment, we assume a minor influence of the DMO procedure that was applied in
Figure 13b, but not in
Figure 14. In general, one can conclude from
Figure 7 that the automatic CMP stack was only successful in imaging event α, and mainly on the left half of the section. Interestingly, the relatively strong reflection imaged above event α, around CMP 140, may be connected to event γ, which is not present in the original CMP stack results. The general picture gets a little bit clearer in
Figure 8, where the CRS operator was used for stacking based on the parameter values estimated by the cascaded search strategy. The improvement of the S/N ratio is notable, which helps to uncover some parts of event α on the right side of the section. Particularly there, on the right side above event α, traces of event γ emerge from the noise. The enhancement obtained from a subsequent local three-parameter optimization is clearly visible in
Figure 9. Reflections become slightly more pronounced, and the traces of shallow events show up above event α. Going on to
Figure 10, the improved lateral event-continuity resulting from the three iterations of event-consistent smoothing of the stacking parameters followed by local three-parameter optimization is notable. In
Figure 11, where the CRS stack result using the hybrid diffraction/reflection search algorithm is displayed, the picture changes relatively drastically. Event γ becomes quite clearly observable over the whole section, except for the area between CMP 100 and 130. In addition, event α can be traced now through the whole section. Even the gap around CMP 240 is covered, but is still not as well defined as the rest of event α. Finally, the stack section displayed in
Figure 12, resulting from the full simultaneous three-parameter search, is the best image we obtained for this data, but outperforms the stack result displayed in
Figure 11, only in tiny details. In both
Figure 11 and
Figure 12, an ultra-shallow horizontal event, which might belong to seismic unit A, can be seen throughout the whole section except in the area between CMPs from 140 to 160. This gap in the continuity of the event is probably due to the bad quality of seismic signals, which not even the simultaneous three-parameter CRS stack was able to improve. However, as mentioned in the original publication of [
46], these very early traveltimes are not very reliable since the usable stacking aperture in offset direction is extremly narrow. Furthermore, comparing
Figure 12 (or even
Figure 11) and
Figure 13a, it is obvious that, relative to the NMO/DMO stack, the CRS stack is better able to image the internal layering of unit B, at least beyond CMP 160. In the first part of the section the internal structures of unit B are not as evident, probably due to a lack of data resolution and/or low data quality.
Coming to the NMO velocities, we will start with
Figure 14, which displays NMO velocities calculated from CRS parameters estimated using the cascaded search strategy. Here, we have first to state that, in this section, the deep blue color on the upper right edge is used as no-value flag, while the deep brownish color in the rest of the section can be seen as the no-useful-value flag, since the coherence analysis tends to favor very high NMO velocities, i.e., a flat traveltime assumption, in the absence of detectable events. Having said this, we can state that, in
Figure 14, well defined velocities are mainly available for event α and predominantly on the left side of the section. The picture could even lead to the wrong impression that there would be a steeper reflector below the event α pinching out at CDP 280. This continues to be true for
Figure 15, where the NMO velocities after cascaded search plus three iterations of event-consistent smoothing, each followed by a local three-parameter optimization, are displayed. A certain success of the strategy is clearly observable, as could be expected from the improved stack result. Nevertheless, there are still too many samples for which no proper value could be determined despite the elaborate optimization strategy.
Figure 16 and
Figure 17 shows a totally different picture. Searching directly with a spatial operator that spans over several CMP gathers creates a continuity of the NMO velocity values, which is consistent with the integral nature of this parameter that makes drastic variations from one sample to the other impossible. The dark brownish samples at the borders of the section still correspond to no useful values.