ABSTRACT

Land seismic data is contaminated commonly with coherent and high amplitude back-scattered noise generated from roughness in the surface topography; these events make generally the key information on deeper layers unclear. In Colombian foothill areas with rough topography and high lateral velocity variation, this noise has been difficult to attenuate. Conventional methods aimed at this purpose usually yield unsatisfactory results. We present a strategy based in prediction and subtraction of the unwanted waves. Assuming knowledge of the source wavelet and the shallow velocity model we use a finite-element solution of the acoustic wave equation to model the back-scattered noise; this modeled response is then subtracted from the prestack data, resulting in a noticeable attenuation of noise in field seismograms.

The method was applied to prestack real data from colombian foothills in order to observe the enhancement in seismic records, planning in a close future to show results on stacked data.

Keywords: dispersion, seismic, topography, foothills, seismogram.

RESUMEN

Comúnmente, los datos sísmicos adquiridos en tierra están contaminados por ruido coherente y dispersivo de gran amplitud generado a partir de las rugosidades de la topografía, que opaca la información proveniente de las capas más profundas. En las zonas de piedemonte en Colombia, con topografía agreste y alta variación lateral de velocidad, este ruido ha sido difícil de atenuar. Los métodos convencionales para atacar este tipo de problemas usualmente dan resultados poco satisfactorios. Con el propósito de contribuir a la solución del problema, presentamos una estrategia basada en la predicción y remoción de las ondas no deseadas. Conocidos la ondícula de la fuente y el modelo de velocidad de las capas someras, se simula el comportamiento de la onda acústica mediante un algoritmo en elementos finitos para reproducir la respuesta de estas capas. Posteriormente, sustraemos tal respuesta de los datos preapilados, obteniendo una atenuación sustancial del ruido en los registros de campo.

El método se aplicó a datos preapilados del piedemonte colombiano a fin de observar una mejoría de la imagen en los registros sísmicos, planeando en un futuro cercano mostrar resultados en secciones apiladas.

RESUMEN

O método foi aplicado a dados pré-empilhados do pé de monte colombiano a fim de observar uma melhoria da imagem nos registros sísmicos, planejando em um futuro próximo mostrar resultados em sessões empilhadas.

INTRODUCTION

As a result from the findings of important oil reservoirs, exploration of hydrocarbons in the Colombian foothills has increased during the last decades. These areas are characterized by difficult conditions of the terrain, like rugged and steep topography, complex lithology and lateral velocity variation in the shallow layers of the subsoil. The seismic images are usually of lesser quality when compared with those from neighboring zones with milder conditions. The presence of coherent and high amplitude backscattered noise in the field data, as generated by the roughness of the terrain and multiple reflections occurring between the land-air interface and the low speed layer, causes poor-quality images. Due to the fact that the contrasts in velocity in these interfaces are usually high, the acoustic impedances thereof are high, thus generating a coherent noise of great amplitude and behaving like a guided wave trapped in the low speed layer. On account of its high amplitude and extensive duration, this coherent noise darkens - totally or partially- the signal coming from the deepest horizons of interest.

In the conventional processing, attempts are made to eliminate or attenuate these events of noise through the application of methods such as in-frequency or multi-channel filtering (Yilmaz, 1991). Attempts are usually made to eliminate topography-related effects, through the application of static corrections to take the source and detectors to a level or datum (Yilmaz, 1991). This procedure is successful when the wave propagates vertically; that is, when the elevation is small and detectors are close to the source. Another approach used to remove the topography effect is to continue the wave field, up to a certain level, through a procedure called "wave equation datuming" (Bevc, 1996); or, to migrate the shot records from topography before stacking (Yang et al., 1999). In foothill areas, results from the application of these methods are usually unsatisfactory. This is attributed to the aforementioned conditions of extreme variability in the medium, when some conditions that are required or expected are not met.

In order to overcome the back-scattered noise and to follow-up on some recent papers (Ernst et al., 2002; Guan et al., 2000; Fu et al., 1999), an alternative strategy is hereby introduced for the elimination or attenuation of noise events. Such events are estimated via modeling, by means of a solution for finite elements of the acoustic wave equation and further elimination by data subtraction. The algorithms of finite elements are particularly suited to the presentation of physical models with irregular geometry and therefore, are an adequate alternative for our purpose.

Because of the deterministic nature of the proposed method, this is based on two premises: first, knowledge of the in-depth velocity model for shallow layers; supposedly, this model consists of a weathered low-speed layer upon a semi infinite underlying layer. The second premise is the appropriate estimation of the source wavelet in each seismic record. We also suppose that the coherent noise of interest is generated only by wave propagation in shallow layers, so it can be separated from other events of interest present in the data. Furthermore, we believe that we don't need to include in-depth information on the model's deeper layers, in order to predict and further extract noise events.

A trial was carried out for the developed method, using seismic data from a Colombian foothill region, obtaining significant noise elimination levels at field record levels, and significant increases in reflector-related information.

THEORETICAL BASIS

The finite-element method is a widely accepted numerical procedure used to solve differential equations of physics-mathematics; the method uses an integral formulation to generate an algebraic system of linked equations, one for each element (Segerlind, 1984). According the Galerkin method (Zienkiewicz, 1992) in each finite element, the wave equation is stated as Helmholtz equation:

Where

[M] :  represents the so called mass matrix

[K] :  the rigidity matrix

{ f } :  is the force or source vector

{u} :  the wave field in the element (Segerlind, 1984)

Equation 1 is solved through Euler's implicit scheme (Langtangen, 1999), obtaining in each model element, the wave field at the time step l:

In (2) Δt is the sampling interval, {u}^l, {u}^l-1 and {u}^l-2 are the wave-field solutions in the current (l), previous (l-1), and prior to previous (l-2) time intervals. The recursive scheme stated in Equation 2 is consistent for values l > 0, but for l = 0 it is not; therefore it is necessary to take into account the initial conditions. Therefore in addition to Equation 2, the discrete form of the initial condition ¶{u} / ¶t = 0 has to be applied, thus implying that.

An algorithm that intend to simulate the wave behavior by finite elements must set initial wave field values for every model element, according Equation 4, and further apply expression (2) each time interval to update the wave field value for every model element.

With the purpose of minimizing reflections on the model borders, a damped zone was defined (Sarma et al., 1998). The acoustic wave equation in this zone contains an additional term of attenuation in Equation 1, i.e.:

Where [C] is the dampening matrix. The solution in the wave equation in the attenuated zone is given by:

The determination of the dampening matrix it is hard to get by mathematics; for this reason, it is approximated by a linear combination of the rigidity and mass matrixes [C] = a [M] + b [K] (Zienkiewicz, 1992). In this paper the a = 11,67; b = 0,00065 values were used, which correspond to the optimum values previously estimated (Sarma et al., 1998).

The set of Equations 2, 3, and 5 represents a recursive procedure. This algorithm was implemented in C++ under the Diffpack version 1.4 tool, on a Linux platform, according to the guidelines set forth by Langtangen (Langtangen, 1996). The code was partially modified to model a wave transmission through a non homogeneous media, so as to allow the use of an external wavelet as source when simulating real seismic experiments. To shorten the program's execution times and to limit memory requirements, the software stores the wave field values in each time interval, only in those nodes where geophones are located. To guarantee the solution's stability found for each discrete model element, it is required that a relation between the Δx, Δy grid dimensions for all the model's finite elements, the sampling interval Δt and the medium velocity field V(x,y) is duly satisfied. This condition, called the Courant condition, establishes the following:

MODEL ESTIMATE

The model's estimation is carried out via first breaks' tomography. The outcome is made up of the stratum geometry and the velocity field. A seismic line was taken from the Colombian foothills in the Department of Casanare. From such a line, a 1 km segment was selected which featured 24 shots and 66 reception points for each. Five out of these shots were selected. Based on this data, a preliminary model was estimated, taking into consideration depth and uphole time. This model may be updated through a refraction static calculation procedure in the ProMAX software. Finally, the model was fine-tuned through the application of refraction tomography; this last stage demanded Hampson & Russell's GLI3D software.

According to the line's topographical data, the obtained model has a rugged surface and a shallow layer over a semi infinite medium (Figure 1). The first velocity is associated to the weathered layer while the second one is related to refraction velocity.

The model was divided into finite elements, each one the size of Δx = Δy = 7 m , resulting in a total of 70 350 nodes and 69 774 elements. An average velocity of 1500 m/s was assigned to the shallow layer while a velocity of 2000 m/s was appointed to the second stratum. These values together with the sampling interval (0,002 s) satisfy Courant's conditions, as stated below:

WAVELET EXTRACTION

The response simulation of the shallow layers must have a source that disturbs the medium; to do so, the form of the wavelet must be known. The real shots were fired with dynamite, so the source wavelet features a minimum phase and therefore, possesses a unique spectrum. The seismograms are usually noise contaminated (groundroll, air wave, etc); therefore, the source spectrum is going to be different from the extracted wavelet and, in consequence, the estimated form differs from the source's own form. To attenuate these events without affecting the data phase, a zero-phase bandpass filter was applied. The process to extract the wavelet implies averaging the spectrum from all the traces within a window to generate its spectrum. The result is transferred to the time domain and the wavelet is presented in trace form. Therefore, the spectrum of a minimum-phase wavelet can be obtained, since the phase spectrum and amplitude (in logarithmic scale) are pairs of a Hilbert transform (Kou, Yuan et al., 2000).

The procedure was applied with a noise factor of 0,1% and an operator length of 120 ms, obtaining the form displayed in Figure 2 for seismogram No. 179.

APPLICATION

The program designed in C++ under Diffpack simulates the propagation phenomena, calculates the wave field in every model's element every 2 ms, and stores only the wave field from the detecting nodes in a file.

The strategy was tried by applying the algorithm on the selected seismograms from a foothill 2D line. The line was fired with dynamite with a group interval of 15 m and a minimum offset of 7,5 m. Each shot was recorded with 600 channels. Although the acquisition time was 5 s, for this trial only 1 s was considered. The objective consisted in observing qualitatively the response attenuation of the shallow layer in various seismograms. The line's seismic processing was done via Landmark's ProMAX software. The wavelet extracted from the field record NO. 179 (Figure 2), was introduced for the execution of the program running on the model.

Figure 3 shows seismogram No. 179 with AGC after applying a bandpass filter to attenuate ground roll and high frequency noise.

A great content of dispersive coherent noise is observed.

The response obtained from the simulation record. Results are shown in Figure 4.

Several events are observed in the shallow part and the influence of the direct wave is evident, although it is easy to eliminate through the stacking process. As for the other seismograms, the process outcomes where consistent with that obtained for seismogram No. 179.

A noticeable difference is observed between original data in Figure 3, as affected by the strong presence of the dispersive coherent noise, and data in Figure 4, following the application of the proposed procedure. The dispersive noise eliminated from data can be seen in Figure 3 as a series of events with linear trends, occasionally called "reverberations". These replicate the first arrival tendencies throughout the whole seismogram. The absence of such events in data included in Figure 4 helps to appreciate other events of greater interest which were previously overshadowed by the presence of a spread out noise.

CONCLUSIONS

• Multiple dispersions and reflections caused by a rugged topography deteriorate the quality of seismic images. To remove this effect from those seismograms achieved under such conditions, a method based on the solution of the wave equation via finite elements is presented.
• Results from this assay featuring synthetic data support the validity of the response prediction strategy for the shallow layer and its further extraction from the overall seismogram.
• The method proved to be effective in reaching the objective, and robust in solving the problem, thus allowing for the clarification of the deepest reflectors.
• Through the application of the procedure to pre-stack seismic data, an improvement of the image could be observed in certain sectors of the seismogram. The multiple reflections and dispersion due to the topography disappear, highlighting the deep reflections.
• By using the wave theory to describe this dispersion phenomenon, the method can be applied to 3D data, for which the application is also designed; in such a case, it is necessary to create a 3D model in finite elements.
• The strategy being deployed and presented in this document is based on two suppositions: first, to define the model's geometry and its interval velocity field; and, secondly, to extract the wavelet from the field data source.

The authors wish to express their gratitude to the Universidad Nacional de Colombia for its support to the development of this project, which is co-sponsored by the Instituto Colombiano para el avance de la ciencia y la tecnología "Francisco José de Caldas" - Colciencias  through contract RC-271-2000; and to ICP - Ecopetrol, for its involvement in this project.

Also, the authors extend their words of gratitude to Petrobras, for providing field data for the trials; to Halliburton LatinAmerica S.A., for lending the ProMAX software; and, to Infopetrol Ltda., for its logistics and conceptual support throughout the research.

BIBLIOGRAPHY

Bevc, D., 1996." Flooding the topography: wave equation datuming of land data with rugged acquisition topography". Geophysics, 61 (5): 1558-1569 .         [ Links ]

Ernst, F., Herman, G. C. and Ditzel, A., 2002."Removal of scattered guided waves from seismic data: geophysics". SEG, 67 (4): 1240-1248 .         [ Links ]

Fu, L. Y.,Guan, H. and Wu, R. S., 1999. "Removing rugged-topography scattering effects in surface seismic data". 69th Ann. Internat. Mtg: SEG, 453-456.         [ Links ]

Guan, H., Wu, R. S. and Fu, L. Y., 2000. "Removing scattering effects of rugged topography using finite-difference method". 70th Ann. Internat. Mtg: SEG, 2189-2192.         [ Links ]

Kou-Yuan, H. and Shen-Pyng, W., 2000. "Neural networks for seismic wavelet extraction and clustering". 70^th Ann. Internat. Mtg: SEG, 741-744 .         [ Links ]

Langtangen, H. P., 1999. "Computational partial differential equations: numerical methods and diffpack programming". Springer Verlag.         [ Links ]

Langtangen, H. P., 1996. "Efficient element solution of the linear wave equation in diffpack". The Diffpack version 1.4 Report Series, SINTEF, University of Oslo .         [ Links ]

Sarma, G. S., Mallick, K. and GadhinGlajkar, V. R., 1998. "Nonreflecting boundary condition in finite-element formulation for an elastic wave equation". Geophysics, 63 (3): 1006-1016 .         [ Links ]

Segerlind, L. J., 1984. "Applied finite element analysis". John Wiley & sons.         [ Links ]

Yang, K., Wang, H. and Ma, Z., 1999. "Wave equation datuming from irregular surface using finite difference scheme", SEG Technical Program with Biographies, 69th Annual Meeting, Houston, Texas, 1465-1568.         [ Links ]

Yilmaz, O., 1991. "Seismic data processing". SEG, Investigations in Geophysics, (2) .         [ Links ]

Zienkiewicz, O., 1992. "El método de los elementos finitos". Vol. 2, McGraw Hill.         [ Links ]