ABSTRACT

RawX-STAR sub-bottom profiler data are stored on tape in a native format. From this format the envelope can be easily obtained. The envelope of the data, which contains amplitude information only, is often used for visualisation. However, it is not suitable for further advanced digital signal processing. Fortunately, the native format can also be converted into a time series, containing amplitude as well as phase information. The data quality can be improved considerably in terms of resolution -both in a temporal and a spatial sense- and signal to noise ratio by applying signal processing techniques that are well known in the field of seismic exploration.

Examples of such techniques are amplitude recovery, which corrects for various amplitude-loss mechanisms, deconvolution, which enhances the temporal resolution, various types of noise reduction, and seismic migration, which increases the spatial resolution and puts dipping reflectors on their correct spatial locations. Once these processes have been applied, a more accurate qualitative and quantitative interpretation can be carried out, which is illustrated here with profiles recorded in various Dutch waterways.

INTRODUCTION

Today, the most advanced algorithms for processing echo-acoustic data can be found in the field of seismic oil exploration. Obviously the huge costs associated with drilling wells have pushed the development of acoustic imaging techniques to their limits. Such advanced techniques cannot be carried out in real-time. Many of these techniques are based on multi-trace processing, i.e. to compute the result for a specific location, data from neighbouring locations is required. Furthermore, skilled personnel are necessary to make proper choices with respect to the parameters that guide the signal processing.

The algorithms that are in use in the seismic industry can be applied to data acquired with sub-bottom profilers without significant problems. A necessary condition is that the profiler data is available as a full waveform time series. Systems that store a sampled envelope signal don’t satisfy this condition. Although the X-STAR sub-bottom profiler stores the raw data in a native format from which the envelope can be easily obtained, it turns out that after some processing the time series can be made available as well. The X-STAR data can therefore be further processed with seismic processing techniques. In this paper the improvements of the X-STAR results that can be realised by applying algorithms from the seismic exploration industry will be discussed and illustrated with several examples.

X-STAR SUB-BOTTOM PROFILER

The EdgeTech X-STAR full spectrum digital sub-bottom profiler is a wideband FM high resolution profiler, generating cross sections of the seabed and collecting normal incidence reflection data over many frequency ranges. It transmits an FM pulse that is linearly swept over the frequency range, the so called ‘chirp’. The acoustic signal subsequently received at the hydrophone is matched and filtered with the outgoing pulse for optimal resolution. It uses proprietary special amplitude and phase weighting functions for the transmitted pulse and a pulse compression filter that maximises signal to noise ratios.

Because the FM chirp is generated by a digital-to-analogue converter with a high dynamic range and a transmitter with linear components, the energy, amplitude and phase characteristics of the acoustic chirp can be precisely controlled to produce high repeatability and signal definition.

There are several towfish available, operating at various frequency ranges between 0.5 and 24kHz. The chirp length is between 20 and 40msec. The A/D converter operates at a sample frequency between 16 and 48kHz and uses a 16 bit sigma-delta technique. 

DATA STORAGE

The measured data is stored on tape (or disk) in the so-called SEG-Y file format, a digital tape standard known from the exploration geophysics industry. This standard was developed in 1974 by the Society of Exploration Geophysicists, so that seismic data and the like (e.g. acoustic profiler data or ground penetrating radar data) can be exchanged easily between different computers. In this SEG-Y format, the data measured at every spatial location is stored as a separate trace. Each trace is preceded by a 240 byte long trace header, in which important attributes of that trace are stored. Examples of such attributes are the spatial coordinates, date and time of the measurements, frequency range of the source, sampling interval, etc.

Every SEG-Y file is preceded by two tape headers of 3200 and 400 bytes in length, in which more general information about the dataset at hand is stored. The first header contains user supplied comment in ASCII format, the second header contains general header information, such as the sample format code, in binary form. In the case of X-STAR data, and as sample format code three states, the samples themselves are stored in a 2 byte IBM integer format. Per trace there are 3976 samples, of which only the first 3200 are live. Each trace (including its header) has a length of  (3976 x 2 + 240 = 8192 bytes).

SIGNAL RECONSTRUCTION

The X-STAR is capable of storing data in processed or raw form. Usually the data is kept in raw form. In this case, the data is not stored as a time series but as the real and imaginary parts of an analytic signal. From this analytic signal, the envelope can be constructed easily by taking the square root of the sum of the squares of the real and imaginary part. The envelope is normally used for display on screen and on camera plots. In this case, from 3200 live analytic samples, 1600 envelope samples are obtained. The sampling interval is then twice the original sampling interval of the A/D converter. Envelope samples are always positive and contain no phase information.

It is however possible to reconstruct the original time signal from the analytic data by means of the Hilbert transform. In this case both amplitude and phase information is obtained and the sampling interval remains the same as the original sampling interval of the A/D converter. Software to accomplish this has been written by the first author as a module to be used with the Seismic Unix package (freely available from the Center of Wave Phenomena of the Colorado School of Mines). 

Fig. 1: analytic samples (left) and reconstructed time trace and envelope (right) (twt – two way travel time)

Figure 1 shows an example of the different data and signal types discussed above. It is a very small part (only 45 time samples) of one measured trace. The real and imaginary parts of the original analytic trace can be seen left. The sampling interval is 62ms. From these samples it is possible to calculate the envelope trace with the same sampling interval, as seen on the right. It is also possible to construct the time signal with a sampling interval of 31ms, the same as the interval of the A/D converter

Fig. 2a: profile with envelope traces

Fig. 2b: profile with reconstructed time traces

Figure 2 illustrates a larger part of the same data as used in Figure 1. It is part of a profile measured in the Ketelmeer in the Netherlands and shows the envelope data (2a) and the reconstructed time signals (2b). It is very clear that the reconstructed time section shows much more detail than the envelope section.

DATA PROCESSING

The advantage of envelope data is that the reflection strength can easily be made visible on a computer screen or on a camera plot. EdgeTech has incorporated a number of interactive processing and interpretation tools in their system so that real-time data analysis can be performed. The disadvantage of envelope data is that phase information is absent and therefore further signal processing cannot be performed. Also the temporal/depth resolution is relatively low. 

An important advantage of reconstructed time signals is that they inherently have a higher temporal resolution. Not only because phase information is now present, but also because the time signals have the same sample interval as the A/D converter, whereas the envelope data has a sample interval twice that. However, the most important advantage is that it is now possible to enhance the quality of the data by means of various digital signal processing techniques as used in the field of seismic data processing.  Use can be made of commercially available powerful software packages, specially designed to process, analyse and interpret this kind of geophysical data. This way, significant improvement of the interpretability can be expected.

Examples of such processing steps are:

AMPLITUDE RECOVERY

Many processing techniques rely heavily on the extraction of statistical information from stationary time series. Roughly similar amplitude levels throughout the data are therefore required. Also for interpretation purposes, the data should be brought into a displayable dynamic range. Due to various amplitude-loss mechanisms, such as spherical divergence, absorption, scattering and transmission losses, the amplitudes in the measured data tend to strongly decrease with recording time. There are various amplitude decay compensation or equalisation techniques available to restore this. They are either data independent, such as a constant gain function, or they are data dependent, such as automatic gain control, whereby scaling is applied dependent on the energy in a range of overlapping time windows. In the latter case however, lateral amplitude variations which can contain valuable information, may get lost.

Fig. 3: profile with reconstructed time traces, unprocessed

Figure 3, shows part of a profile recorded in the Westerschelde river in the Netherlands. The amplitudes are not laterally balanced, giving it a stripy appearance. After amplitude balancing and amplitude decay compensation, as depicted in Figure 4, the section reveals more detail in depth.

DECONVOLUTION

Improving the temporal resolution by reducing the effective length of the pulse (or wavelet) can be achieved by deconvolution techniques. In this process the wavelet is changed (generally shortened) and its amplitude spectrum is shaped (generally broadened) so that the frequency content increases and there is less ringing. This enhances the overall resolution capability of the wavelet. 

An effective implementation (in terms of robustness and ease of use) of signal shortening and spectral broadening is statistical deconvolution using Wiener prediction-error filters. In other words, a filter is derived which will predict the tail of the (ringing) wavelet and then subtract it from that wavelet, thus effectively shortening it. The position where the tail starts is the prediction distance or gap and can be retrieved from the autocorrelation of the signal. As a rule of thumb, the second zero crossing is chosen as the gaplength.

Fig. 4: profile after amplitude processing

Figure 4 shows the original profile after amplitude restoration. The sea bottom reflection indicates much ringing and also that the event at trace no. 550 at 17ms is not well resolved. By means of predictive deconvolution the vertical resolution is enhanced. From the autocorrelation, as

Fig. 5a: autocorrelation of profile before deconvolution

Fig. 5b: autocorrelation of profile after deconvolution

shown in Figure 5a, the prediction lag of 0.2ms was chosen and the filterlength was set at 1ms. After deconvolution, the autocorrelation is plotted again in Figure 5b and clearly shows much less ringing. This is also expressed in the amplitude spectrum before (Figure 6a) and after (Figure 6b) deconvolution. The broader and flatter spectrum after deconvolution indicates that a shorter wavelet is now present. The outcome of the deconvolution on the profile is seen in Figure 7. The sea bottom reflector is much sharper now and the event at trace no. 560 at 17ms is well resolved. 

Fig. 6a: amplitude spectrum of profile before deconvolution

Fig. 6b: amplitude spectrum of profile after deconvolution

Fig. 7: profile after deconvolution

NOISE REDUCTION

Several noise reduction schemes can be applied to the data in a 1D or 2D sense, aiming at reducing random and/or coherent noise and enhancing the primary reflection energy. Bandpass filtering is an example of a single channel process, which can be applied when outside the signal band, noise is so dominant compared with the decaying signal content that elimination of these frequencies will reduce the overall noise content whilst preserving wavelet resolution. But often signal and noise share a common bandwidth and more sophisticated methods, working in a 2D sense, are required to enhance the data quality. 

In the example shown below fxdecon, a particular type of deconvolution implementation, followed by median filtering was applied to enhance the spatial resolution whilst reducing the random noise. Fxdecon applied a Fourier transform to each trace of a dataset, applied a complex Wiener prediction filter along the horizontal direction for each frequency and then inversely transformed each resulting frequency trace back to the time domain. The signal was assumed to have similarity from one trace to the next, i.e. it was spatially predictable. Differences between the predicted waveform and the actual one can be classified as noise and removed.

Next, a 2D mean filter was applied. Per time sample, it sorted the sample values and averaged a range of values centred around the median. This way, outlier points (noise samples) were discarded and the remaining ones smoothened.

Fig. 8: profile before noise reduction

Fig. 9: profile after noise reduction

Figure 8 shows part of an original profile as recorded in the harbour of Amsterdam. After noise reduction had been applied according to the recipe described above, the result was plotted and Figure 9 was obtained. The reflector continuity had been improved and the noise level had been reduced considerably. To appreciate the quality improvement of processing, the original envelope data of the same part of the profile was plotted in Figure 10.

Figure 10: profile with original envelope data

MIGRATION

The sub-bottom profiler recordings, are used to determine the structure of the bottom from measurements obtained at the surface. The processing of these measurements is aimed at rendering them more comprehensible geologically. Migration, also called Synthetic Aperture Sonar or SAS, is one of these processes. It is an inverse operation involving rearrangement of (seismic) information elements so that reflections and diffractions are positioned at their true locations.  The need for this arises since variable velocities and dipping horizons cause these elements to be recorded at surface positions different from their subsurface positions. In other words, the purpose is to ‘migrate’ the recorded events to their correct spatial positions by backward projection or de-propagation based on wave theoretical considerations. In summary: migration spatially ‘sharpens’ the recorded events and puts them in depth correctly. In exploration seismology, migration is considered as the most powerful imaging algorithm.

The use of migration is constrained however, by the accuracy of the navigation data. The position of every trace must be known with cm precision. Also the depth of the towfish must be known precisely and should be corrected for wave height, swell and tidal influences. Despite recent progress in positioning accuracy, in present day sub-bottom profiling operations this accuracy is usually not met and migration therefore cannot successfully be applied.

INTERPRETATION

After the data quality has been improved by processing, other tools and packages from the geophysical industry can be used to interpret the data. For example the water bottom can be picked automatically, reflector strength can be determined, deeper horizons may be tracked and various other attributes can be determined. After the integrity of the data, and the information derived from it has been secured with QC tools it is stored in a database for easy retrieval. The final products of the interpretation may be delivered in the form of maps or exported to spreadsheet format.

Figure 11: shale layer depth map

As an example, in Figure 11 a depth map of a shale layer beneath the water bottom is plotted from the data recorded in the harbour of Amsterdam. The reflection of this layer can be seen in Figures 8 and 9 at around 17ms. Semi-automatic picks of this layer (where present) were taken at some 79,000 traces and the picked times, together with the coordinates were stored in a database, from which the map was plotted. The area is approx. 4.3 by 5.6km.

CONCLUSIONS

It has been shown that the quality of the original (envelope) measurements of the X-STAR sub-bottom profiler can be further enhanced. Therefore the raw data samples, stored in an analytic form, should be transformed to real time samples which contain amplitude and phase information. Then seismic data processing techniques, readily available in industry standard software packages, can be applied to improve the data quality. Examples of such processing steps have been illustrated with real data from profiles recorded in various Dutch waterways. Subsequently, detailed interpretation of the data can be carried out, again using specialised software which ensures information quality and exportability. 

The fact that it is possible to reconstruct the original time signal from the recorded data gives the X-STAR a distinct advantage over systems whereby only the envelope is sampled and stored.

SUGGESTIONS FOR FURTHER READING 

The following references are suggested to further enhance the reader’s knowledge on the principles of geophysical data processing:

•    Ozdogan Yilmaz, Seismic Data Processing, 1987, The Society of Exploration Geophysicists, isbn 0-931830-40-0

•    Proceedings of the IEEE, Vol. 72, No 10, October 1984. Special issue on seismic signal processing.

•    Sheriff, R.E., Encyclopedic Dictionary of Exploration Geophysics, 1991, third edition, The Society of Exploration Geophysicists, isbn 1‑56080-018-6

Furthermore the following websites are recommended to be visited:

• www.edgetech.com (EdgeTech)

• www.cwp.mines.edu/cwpcodes (Colorado School of Mines)

• www.seg.org (Society of Exploration Geophysicists)

ACKNOWLEDGEMENTS

The authors wish to thank Rijkswaterstaat for providing the datasets used in this article. We especially thank Dr. N.A. Kinneging of the Survey Department of Rijkswaterstaat for his support and advice on the project. 

Remco Romijn

received an MSc in geophysics from Delft University of Technology in 1983. He then worked as a seismic data processor for an oilfield service company. In 1994 he joined the geophysics/petrophysics group of the faculty of applied earth sciences of the TU Delft. Since 2000 he workes at TNO-TPD in the field of ultrasonic transducers and systems

Gerrit Blacquière

received a PhD (1989) in geophysics from Delft University of Technology. He spent several years in the seismic industry before joining TNO TPD. Since 1999 he has been a part-time professor of acoustic instrumentation in Delft. He has been involved in many projects ranging from seismic processing to hydrographic surveying.