Geophysical Prospecting, 1998, 46, 105–120
Very high resolution 3D marine seismic data
processing for geotechnical applications1
B. Marsset,2 T. Missiaen,2 Y.-H. De Roeck,2 M. Noble,3 W. Versteeg4
and J.P. Henriet4
Abstract
The processing of a small-scale, very high resolution (VHR) shallow marine 3D data
volume is described. The data were acquired over a small clay diapir, on the river
Schelde, in 1990. Using an array of 12 dual-channel microstreamers towed from a
catamaran, a network of 1 m × 1 m bins could be produced over an area of
50 m × 180 m (< 100 m penetration). Positioning was performed with an auto-
tracking laser ranging system, assuring an absolute accuracy of a few decimetres.
Preliminary processing steps included tidal correction and multiple removal. An
important step concerned the application of 3D prestack Kirchhoff depth migration.
Indeed this processing allows easy handling of the exact positions of both source and
receivers as the latter were not set out on a conventional regular grid due to navigation
difficulties. Because of the restricted data volume and the more-or-less stratified
medium, a 1D velocity model could be used. This allowed a considerable
simplification of the migration algorithm, based on summation. Traveltimes were
calculated only once, using a 2D time grid with 0.1 m intervals.
This migration method proved very efficient, greatly improving the seismic image,
and involved only limited CPU time on a small computer (Sparc 10 workstation). It
clearly demonstrates that advanced seismic processing can form a valuable and
economically feasible tool for VHR shallow subsurface 3D seismics, as long as the
velocity field is not too complex. This method should therefore no longer be restricted
to large computers and hydrocarbon exploration, but should also become a routine
for VHR 3D shallow seismic work.
Introduction
Very high resolution (VHR) marine seismic data are widely used both in the field of
engineering geophysics as well as in geological and environmental site investigations.
This technology has progressed quickly. Evolution towards cheaper and more
1
Received March 1996, revision accepted August 1997.
2
IFREMER, Centre de Brest, BP 70, 29280 Plouzane, France.
3
Ecole des Mines de Paris, 77305 Fontainebleau, France.
4
RCMG, Universiteit Gent, 9000 Gent, Belgium.
q 1998 European Association of Geoscientists & Engineers 105
106 B. Marsset et al.
powerful computers is making VHR multichannel recording more popular. Whereas
2D surveys are now in common use, most 3D acquisition is still carried out only in the
hydrocarbon exploration realm. It involves complex field and processing procedures,
making it very costly and time-consuming.
In recent years, the Renard Centre of Marine Geology (University of Gent,
Belgium) has scaled down the shallow marine 3D method to very high resolution,
thereby entering the world of small-scale geological structures. Using a relatively
simple acquisition system, a small 3D data volume (lateral dimensions 50 m × 180 m,
penetration < 100 m) was acquired over a clay diapir under the river Schelde. Time
slices obtained after stacking showed a resolution never achieved previously in 3D
practice (Henriet, Verschuren and Versteeg 1992).
Up to now, the processing sequence of shallow seismic data has often seemed to
stagnate around the 1D assumption of the common depth point. Indeed very few
people in this field go beyond NMO and stack. The sea-bed and its near-subsurface,
however, are, by nature, 2- and often 3-dimensional structures and should be treated
as such. This can be done by applying migration algorithms to the data. This imaging
technique, coming directly from the oil industry, has the well-established reputation
of being difficult to use and requiring large computer resources. Nowadays this is no
longer true, due not only to faster processors but also to new algorithms.
The aim of this paper is not to give a full detailed discussion of the 3D migration
technique. Several authors have written excellent articles and books covering this
topic, and in many cases have done so with considerable rigour. Our motivation for
writing this paper came from the belief that practical advice for the application of 3D
prestack migration for near-subsurface seismic data is in short supply. We prove here
that migration forms an essential processing step involving only moderate
supplementary costs, and it should therefore become a routine procedure for VHR
3D shallow seismic work.
Target
The target consisted of a small diapir under the river Schelde near Antwerp, close to
the site proposed for a metro tunnel. The clay beds are marked by a large number of
concretions, with diameters ranging between 0.5 m and 1 m and thicknesses between
0.2 m and 0.3 m. Both the structural deformation and the presence of concretion beds
motivated a detailed seismic and geotechnical survey of the tunnel site (Schittekat,
Henriet and Vandenberghe 1983). These detailed site investigations made the area a
well-documented test ground to evaluate the capabilities of VHR 3D seismic
reflection imaging in geotechnical applications. This survey was carried out in 1990,
in the framework of an EEC-supported programme for the stability evaluation of
potential offshore sites.
The diapir in the Rupelian (Oligocene—middle Tertiary) clay has an apparent
diameter of 60 m, and a vertical amplitude increasing from a few decimetres at 50 m
depth to a few metres at about 25 m depth. As the lowest horizons do not appear to be
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
VHR 3D marine seismic data processing 107
Figure 1. Schematic overview of the 3D seismic acquisition system ‘SEISCAT’.
significantly deformed below the diapir, the effect of velocity pull-up can probably be
ruled out. The origin of the diapir is thought most likely to be related to decom-
paction, the sediment load being reduced by erosion of the river bed (Verschuren
1992). The concretions in the clay beds, which have a diameter of 0.5–1 m and a
thickness of 0.2–0.3 m, locally stand out as diffraction clusters on analogue boomer
profiles recorded during former seismic surveys (Hemerijckx et al. 1983).
Data acquisition
In order to allow accurate imaging of the diapiric deformation and the larger clay
concretions, the 3D acquisition on the river Schelde was planned with a bin size of 1
metre (Henriet et al. 1992). Since the depth of the target is extremely shallow (about
25 m), this would give moderate oversampling of the first Fresnel zone of 4 m for an
average signal frequency of 1 kHz.
The receiver configuration consisted of an array of 12 dual-channel microstrea-
mers, towed from a source-bearing modified catamaran (Fig. 1). The streamer
spacing and the channel spacing within the streamers were both 1 m. Source and
streamers were towed near the surface and offsets ranged between 6 m and 9 m. The
catamaran was towed some 10 m behind the ship. During the survey a boomer and a
modified watergun source were used. About 20 000 shots were recorded, each shot
generating 24 seismic traces in a band 5.5 m wide along the tracks. Recording was
done with an EGG seismograph, using a sampling rate of 0.25 ms. The shooting rate
was 1 second.
In order to allow advanced processing, very accurate positioning was needed. This
was done with an auto-tracking laser ranging theodolite, set up on the river bank,
which continuously tracked a crown of reflector prisms mounted on top of the
catamaran. This location reference point proved very stable due to the absence of
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
108 B. Marsset et al.
Figure 2. Geometry of source and receiver positions over the area of interest.
wave action (river environment). The measurement of transmission delay was taken
into account when calculating the array positions. This full procedure resulted in a
horizontal and vertical positioning accuracy of a few decimetres.
The strong tidal currents on the river, however, caused serious problems for the
array stability. In order to keep the towed equipment well stretched behind the vessel,
the seismic tracks were sailed against the current, alternately forwards and backwards,
the latter simply by decreasing the speed. This procedure minimized relative
positioning errors of the array. The absolute displacements over the ground varied
considerably, resulting in a shot interval between 0.5 m and 2 m. Due to the currents
and the sometimes busy river traffic, it was difficult to steer along the preplanned
regularly spaced sail-lines. The latter were therefore adjusted in order to obtain a
coverage as complete as possible.
Preliminary processing
The study presented here focused on a 50 m × 180 m area covering the diapir. Only
the boomer data set (2000 shots, total data volume 70 Mb) will be discussed. The
source and receiver locations were calculated from the prism positions taking into
account the ship’s movement, including the skew due to the currents. Due to the
navigation problems mentioned above, the coverage with the boomer data was rather
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
VHR 3D marine seismic data processing 109
Figure 3. Top: example of one shot record after applying a band-pass filter and top mute. Below:
corresponding amplitude spectrum.
variable. Although the coverage over the area of interest was fairly good in most parts,
some areas (especially along the margin) were marked by relatively poor coverage
(Fig. 2).
Data editing
The boomer source has a frequency range from a few hundred Hz up to well over
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
110 B. Marsset et al.
Figure 4. Top: example of one CMP (after NMO correction) before tidal correction. Middle:
corresponding tidal amplitude. Below: same CMP after tidal correction.
2000 Hz. In view of the sampling rate, an antialias filter of 1400 Hz was applied
during the survey. This limitation of the sampled frequency band has undoubtedly
restricted the possible imaging resolution. It probably explains the absence on the 3D
data of the above-mentioned clay concretions, the latter being automatically filtered
out during digital recording. Since analogue recording can involve a much broader
frequency range, it allowed observation of the diffraction hyperbolae on the paper
records.
Band-pass filtering was carried out using a zero-phase Butterworth filter with half-
power frequency of 400 Hz, and 40 dB/decade roll off at the slope. As the first part of
each trace contained direct wave arrivals, a mute was applied before the river-bed
reflection. Figure 3 shows a shot record of 24 traces after band-pass filtering (no
NMO correction), and its corresponding average power spectrum. Due to the small
offsets involved, the variation in arrival times between the different receivers is barely
visible.
The sometimes strong variations in amplitude are due both to the lack of
repetitivity of the source (inherent in broad-band, very high-frequency, seismic
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
VHR 3D marine seismic data processing 111
Figure 5. Stack fold coverage over the area of interest.
sources) and to the shallow immersion of source and receivers which cannot easily be
monitored. Indeed, variations of a few decimetres in the towing depth will cause the
original signal to interfere with its sea-surface ghost leading to frequency-dependent
destructive summation. Therefore, no amplitude corrections were carried out on the
data set.
Tidal correction
Tidal action on the river Schelde was quite large, with amplitudes up to 5.5 m (Fig. 4),
whereas the dominant signal wavelength was less than 1 metre. Therefore the arrival
times were corrected to a common level by adding a timeshift based on the vertical
coordinate of the positioning prism crown (Fig. 4). These corrections were carried
out on the zero-offset data, after the application of NMO assuming a constant velocity
of 1500 m/s.
Multiple removal
The upper part of the clay diapir was marked by a relatively strong river-bed multiple.
Due to the different tidal corrections, this multiple will be largely suppressed in the
process of stacking. Nevertheless multiple removal was carried out through predictive
deconvolution, using a Wiener–Levinson algorithm, applied on a trace-by-trace basis.
The method applied proved to be quite efficient, resulting in a good suppression of
the river-bed multiple.
Stacking
Taking into account the irregular navigation, NMO and stack had to be carried out
using the exact locations of each source–receiver pair. Once the geometrical midpoint
position had been computed for each trace, a stack grid was set out in bins of
1 m × 1 m (Fig. 5). The resulting stack fold coverage appears rather heterogeneous
over the survey area. According to the theoretical acquisition geometry a four-fold
coverage should be obtained. Due to the irregular grid, however, the coverage was
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
112 B. Marsset et al.
Figure 6. Results after stack. The time section is taken from a central row of bins. The time slices
are separated by 0.5 ms.
high in some parts (up to 24 traces per bin), and consequently relatively poor in other
parts.
The lack of large-offset data did not prevent the calculation of velocities using
NMO-like analysis. Moreover, the small offsets and shallow depth of the target led to
extremely small and almost constant Dt over the streamer array (from 4 to 6 samples).
Therefore the NMO was carried out using velocity values extracted from geotechnical
investigations (1500 m/s for the water column and 1600 m/s for the entire sediment
column) (Heldens 1983). Such small variations in the P-wave velocity, as expected
from the shallow depth, have very little effect on the imaging process. The actual
stacking process was carried out per bin, resulting in a total of over 8000 stacks. The
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
VHR 3D marine seismic data processing 113
full sequence of NMO and tidal corrections, geometry computations and stacking
required 35 min on a Sparc 10 workstation.
The stacked data still displayed some high-frequency variations in the first-break
arrivals. In order to correct for these time delays, and having no bathymetric control,
the river-bed reflectors were smoothed by low-passing the first-break spatial
wavelengths that were larger than the diameter of the Fresnel zone estimated for
the average bathymetry, thus the traces were only shifted in time with a constant
value. Although this resulted in a clear improvement in the upper few metres, the data
quality decreased dramatically with depth. This suggests that the observed delays
were most likely related to slight errors in the positioning of the acquisition array.
The results after stacking are shown in Fig. 6. The time section is chosen from a
central row of bins, across the centre of the clay diapir. The time slices were taken
from the lower part of the diapir (indicated by the arrow), and are separated by
0.5 ms. Due to the above-mentioned lack of information concerning the amplitudes,
the time slices are simply represented in black and white, indicating the positive and
negative phases.
The quality of the stacked section is good in the lower part of the diapir, with a large
number of continuous and energetic reflectors visible. The middle and upper parts of
the diapir, however, still remain quite disturbed, marked by discontinuous and
incoherent reflectors. The time slices show the concentric reflector patterns growing
towards the base of the diapir. Time slices above those shown in Fig. 6 were less and less
coherent, which confirms the chaotic reflector pattern observed on the time section.
3D migration
Introduction
In order to obtain an accurate image of the clay diapir, and to improve the
interpretability, the application of 3D migration is essential. The main benefit of
migration concerns the removal of the propagation distortion from the reflectivity and
the repositioning of the reflectors in their correct positions. The resulting
improvement is most obvious for complex geological structures. For our specific
data set acquired over a more-or-less 1D stratified medium, we do not expect
dramatic geometrical changes nor quantitative information about the reflectivity.
A second benefit of migration concerns noise. The treatment of noise by migration
algorithms has sometimes been neglected at the expense of the treatment of the signal,
in spite of the fact that it can be of great importance. Migration reduces the amplitude
of background noise without affecting the signal amplitudes. This is clearly
demonstrated on our data set.
Migration processing is conventionally carried out after stacking. Post-stack
migration is less costly as it involves a reduced volume of data, and is generally quite
robust. In general, prestack migration gives more accurate results as it removes the
cosine dependence of the stacking velocity on dip, forces spatial alignment and allows
events of different dip at the same unmigrated time to stack correctly (Hatton,
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
114 B. Marsset et al.
Worthington and Makin 1986). However, this technique is less used in practice as it is
often difficult to use and involves much higher processing costs due to the larger data
volume. With increasingly powerful computers and new algorithms, the latter has
become less of a disadvantage.
The question remains whether the results justify the costs. Algorithmically, full-
waveform 3D depth migration before stack can be implemented, but it would incur an
excessive computational cost. Fortunately, the objectives of shallow subsurface and
geotechnical seismics are often not as exigent as those required by the oil industry. In
terms of migration this means that we are mainly interested in the kinematics of how
to reposition the reflectors of interest correctly. This allows a few approximations,
which, in turn, allow a reduction in CPU time and the use of smaller computers,
enabling us to circumvent the economic problem of migration.
Migration method
There are several approaches to migration. All have their own advantages and
disadvantages, the proper choice depending upon the type of data to be migrated. For
VHR small-scale 3D work, the main requirement is that it should be economical, and
honour the kinematics. Furthermore, it should be able to handle the complex
acquisition geometry (irregularly sampled spatial grid) easily. Considering these
parameters, the Kirchhoff summation approach was the logical choice.
The Kirchhoff method can be interpreted as summing along diffraction curves. We
look upon every point in the depth section as a possible diffractor. Since it is assumed that
we know the velocity above this point, we can find the apex of its hypothetical diffraction
hyperbola. Then we sum along this curve and place the value obtained at the given point.
When the medium does not have a constant velocity, the latter must be generalized.
Essentially this amounts to replacing the straight ray traveltime by a traveltime
accurately computed for a variable velocity medium. In this particular case, we used
the algorithm introduced by Podvin and Lecomte (1991) which computes the first-
arrival traveltimes by solving the eikonal equation by finite differences. This algorithm
was used for both post-stack and prestack migration processing.
Migration processing requires a velocity model. The latter can normally be
deduced from the stacking velocities obtained through conventional NMO
processing. As mentioned above, this was not possible here due to the very short
offsets which did not allow the picking of reliable stacking velocities. Instead the same
1D stratified velocity model was used as mentioned above for the NMO correction,
based on geotechnical measurements at the tunnel site (Heldens 1983). Although this
certainly represents an oversimplification of reality, the extremely shallow environ-
ment seems to justify this approximation.
Prestack migration
In order to prepare the data for prestack migration, an inverse NMO algorithm was
applied to the tidal-corrected zero-offset data. The overall prestack migration
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
VHR 3D marine seismic data processing 115
Figure 7. Geometry of the prestack migration (t ¼ traveltime).
operation involved the following steps: (1) computation of geometry parameters; (2)
computation of the first time arrivals, and (3) summation migration.
In order to take into account our particular navigation pattern, we have tailored the
Kirchhoff migration to meet our needs. Again, it should be emphasized how easily this
particular migration handles both irregular shot and line spacing.
Prior to the migration itself, a grid with 1 m × 1 m bins is created. For a given aperture,
the maximum offset (source-migrated trace or receiver-migrated trace) is calculated in
order to compute the traveltime grid only once for the 1D velocity model we intend to
use. This 2D time grid uses 0.1 m grid cells; the size of the cells was chosen with respect
to the average wavelength of the seismic signal (less than one tenth).
Figure 8. Post-migration stack fold coverage over the area of interest.
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
116 B. Marsset et al.
Figure 9. Results after post-stack migration. The depth section is taken from a central row of
bins. The depth slices are separated by 0.5 m.
The migration is then carried out on a shot-by-shot basis. For each trace, migrated
traces are calculated for the bins falling within the area centred around the
corresponding CMP and extending up to the given aperture. In the case of an
aperture of 108, this roughly corresponds to 25 bins. Each migrated depth sample is
evaluated as the sum of the 24 different trace values picked at the traveltimes
corresponding to the (source–migrated trace–receiver) triplet (Fig. 7).
Migrated traces issued from different shots and located within the same bin are
then summed together to form the final migrated volume. The resulting fold coverage
was again rather heterogeneous over the survey area, with a maximum of 34 traces per
bin (Fig. 8).
Preliminary tests were carried out using different apertures. The use of a large
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
VHR 3D marine seismic data processing 117
Figure 10. Results after prestack migration. The depth section is taken from a central row of
bins. The depth slices are separated by 0.5 m.
aperture (more than 208) dramatically increases the computing time while not adding
new information to the results. Therefore a small aperture (108) was used to produce
the results discussed here. For the given aperture, the total processing time for inverse
NMO and prestack migration on a Sparc 10 workstation was about 60 min, which
includes the geometry calculations.
Post-stack migration
In order to compare the different migration methods, we also applied a 3D
migration routine to the stacked data volume. The time grid used for the
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
118 B. Marsset et al.
Figure 11. Different processing steps applied to the VHR 3D shallow marine data set.
calculation of the first-arrival traveltimes was the same as mentioned above. We
therefore used an aperture equal to the coverage used for prestack migration.
The full sequence of stacking and post-stack migration required about 45 min on a
Sparc 10 workstation, including NMO correction and the necessary stack geometry
calculations.
Discussion
Figures 9 and 10 show the results after post-stack and prestack migration, in both
cases for an aperture of 108. The input data have been band-pass filtered and
deconvolved to remove the multiple. The depth sections were chosen from a central
row of bins, across the centre of the diapir. The depth slices, separated by 0.5 m, were
taken from the lower part of the diapir (see arrow).
From these figures it is evident that both migration operations clearly improved the
seismic image, even when based on a simple 1D stratified velocity model. In the latter
case, however, the data improvement is mainly due to noise reduction and, to a lesser
extent, the repositioning of the reflectors. The depth sections are marked by a large
number of clearly distinguishable and coherent reflectors, which were still weak and
distorted in the stack. The depth slices agree with these observations, as they are
clearly marked by an improved resolution, with the concentric reflector pattern being
better defined.
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
VHR 3D marine seismic data processing 119
When comparing the post-stack and prestack results, we can see that the latter has
better results. Indeed the continuity of the observed reflectors is better on the prestack
migrated section, especially towards the middle part of the diapir. This is also the case
for the corresponding depth slices, where the concentric reflector pattern is areally
more restricted.
The uppermost part of the diapir, however, still remains somewhat chaotic. This
seems to suggest that this part may have been subject to several gravitational collapses.
Indeed the diapir growth may have induced stresses in the less undercompacted upper
part of the clayey formation, causing overlying horizons to collapse into the initially
undercompacted zone and destroying its internal strucure (Verschuren 1992).
Figure 11 gives an overview of the different processing steps applied to the VHR 3D
shallow marine data set. The processing times involved were far from being excessive,
even when carried out on a modest Sparc 10 workstation. In fact the extra CPU time
for prestack migration was almost negligible compared to the post-stack sequence,
whereas the latter produced less good results.
Conclusions
Very high resolution 3D seismics is of the greatest interest for detailed small-scale
geological and geotechnical site investigations. The acquisition of VHR shallow 3D
data has been demonstrated to be feasible, in a modest and cost-effective way, without
requiring complex field procedures, and with bin sizes of the order of 1 metre. The
processing of this small-scale 3D data volume implied challenging aspects related to
the algorithmic treatment of the data, taking into account both the geotechnical goals
and the processing facilities.
Prestack depth migration processing formed an important step in the processing
sequence of the VHR 3D data volume. The algorithm used here was based on the
Kirchhoff summation approach, using a 2D traveltime grid based on a 1D velocity
model. The latter was possible due to the more-or-less stratified medium and the
restricted data volume involved. With these approximations the prestack migration
method applied proved very efficient, involving limited CPU time on a small
workstation (Sparc 10), and clearly improving the interpretability of the data. This
demonstrates that advanced seismic processing forms a valuable and economically
feasible step in VHR 3D near-subsurface seismics, and there should no longer be any
reason not to apply it.
Acknowledgements
The present study was partly carried out in the framework of an EEC MAST-II
Research Fellowship (Project CT/5024/94). The development of the 3D acquisi-
tion system was supported by the EEC Directorate for Energy (Project TH/06036/
87). The cooperation of the officers and crew of the RV Belgica is gratefully
acknowledged.
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120
120 B. Marsset et al.
References
Hatton L., Worthington M.H. and Makin J. 1986. Seismic Data Processing. Blackwell Scientific
Publications.
Heldens P. 1983. Een seismische studie van de klei van Boom en de klei van Ieper. PhD thesis,
University of Gent, Belgium.
Hemerijckx E., Carpentier R., De Schrijver P., Henriet J.P. and Heldens P. 1983. Preliminary
investigations of concretion horizons in a Tertiary clay layer in view of the construction of a
pre-metro tunnel under the river Scheldt at Antwerp (Belgium). Symposium on Engineering
Geology and Underground construction, Lisbon, Proceedings I, 53–73.
Henriet J.P., Verschuren M. and Versteeg W. 1992. Very high resolution 3D seismic reflection
imaging of small-scale structural deformation. First Break 10, 81–88.
Podvin P. and Lecomte I. 1991. Finite difference computation of traveltimes in very contrasted
velocity models: a massively parallel approach and its associated tools. Geophysics 105, 271–
284.
Schittekat J., Henriet J.P. and Vandenberghe N. 1983. Geology and geotechnique of the Scheldt
surge barrier – characteristics of an oversolidated clay. Proceedings of the 8th International
Harbour Congress, Antwerp (B), I.121–I.134.
Verschuren M. 1992. An integrated 3D approach to clay tectonic deformation. PhD thesis,
University of Gent, Belgium.
q 1998 European Association of Geoscientists & Engineers, Geophysical Prospecting, 46, 105–120