Application of wavelet theory to analyze seismic signs in the space time-frequency 

Zaida Chirinos Fuentes¹, Juan Infante², Montoya Dan El² y Marcos Figueroa² 

¹ Procesamiento de Señales, Ingeomin, Caracas, Venezuela.

² Universidad Central de Venezuela (UCV). Caracas, Venezuela. Telf. 0212-5970826. zaidachirinos@yahoo.com

Abstract 

A seismic signal can be represented like the sum of geologic signal and noise. The noise causes problems in the geophysical data processing because it limits the identification of the geologic characteristics. In this sense, it is important to reduce the associated noise to the measurement trying to preserve the useful information. In Geophysical, the deconvolución is the technique commonly used in the seismic data processing; this technique is based on the theory of filter of Wiener-Levinson and is used to attenuate or to eliminate manifolds or reverberations that appear like noise in the signal. In this work a technique is presented that increases the performance of the Wiener filter in the dominion of the frequency by means of the application of a mixed scheme based on the technique of thresholding of the coefficients of the Transformed of wavelet and the Wiener filtrate in the dominate of wavelet. In order to evaluate the performance of the proposal, the results are compared in the dominion of the frequency with the classic Wiener filtrate. 

Key words: Wavelet transforms, data processing, seismic data, Wiener Filter.

Aplicación de la teoría de ondículas para analizar señales sísmicas en el espacio tiempo-frecuencia 

Resumen 

Una señal sísmica puede ser representada como la suma de señal geológica y ruido. El ruido causa problemas en el procesamiento de datos geofísicos porque limita la identificación de los rasgos geológicos. En este sentido, es importante reducir el ruido asociado a la medición tratando de preservar la información útil. En Geofísica, la deconvolución es la técnica comúnmente usada en el procesamiento de datos sísmicos, esta técnica está fundamentada en la teoría de filtro de Wiener-Levinson y es usada para atenuar o eliminar múltiples o reverberaciones que se presentan como ruido en la señal. En este trabajo es presentada una técnica que incrementa el desempeño del filtro Wiener en el dominio de la frecuencia mediante la aplicación de un esquema mixto basado en la técnica de umbralización de los coeficientes de la Transformada de Ondícula y el filtrado Wiener en el domino de la ondícula. Para evaluar el desempeño de la propuesta, se comparan los resultados con el filtrado Wiener clásico en el dominio de la frecuencia. 

Palabras clave:  Transformada de ondícula, procesamiento de datos, señal sísmica, filtro Wiener.

Recibido el 20 de Noviembre de 2007 En forma revisada el 09 de Febrero  de  2009 

Introducción 

Cuando se miden señales sísmicas, éstas se presentan perturbadas o contaminadas por lo que suelen representarse a través de la suma de señales geológicas puras y ruido. El análisis de una señal sísmica con ruido puede entorpecer la identificación de los rasgos geológicos. En este sentido, es importante reducir el ruido asociado a la medición tratando de preservar la información útil. El ruido se define, según Yilmaz [1], como aquella interferencia que pueda distorsionar la relación de proporcionalidad, de aquellos eventos, que son reflexiones primarias del subsuelo. En Geofísica, la deconvolución es usada para atenuar o eliminar múltiples o reverberaciones. Con ello se logra mejorar la forma de la ondícula, disminuyendo su anchura y haciéndola más resolutiva. La deconvolución, idealmente, comprime las componentes de la ondícula eliminando los múltiples y dejando sólo la reflectividad de la tierra en la traza sísmica. Según Foufola [2], la energía de los múltiples o reverberaciones es ruido de tipo coherente y es fuente de confusión en la interpretación sísmica. La técnica de deconvolución comúnmente usada en el procesamiento de datos sísmicos se basa en el filtro de Wiener-Levinson. Este filtro, en teoría, provee el mejor método lineal para remover ruido gaussiano estacionario sumado a un proceso lineal.

La Transformada de Ondícula ha sido usada exitosamente para suprimir el ruido de señales sísmicas. Así lo reflejan los trabajos de Charneley [3] y Zhou [4]. En estos trabajos el ruido es atenuado suavizando los coeficientes de detalle de la Transformada de Ondícula Discreta de la señal mediante umbralización. Esta técnica ha reportado mejores resultados que el filtrado pasabanda en el dominio de la frecuencia. Siyuan [5] implementó la Transformada de Ondícula en el filtrado del ruido de señales sísmicas a través del esquema realzado (Lifting Scheme Transform) conocido como el análisis ondícula de segunda generación, utilizando predicciones espaciales lineales y no lineales. Yan [6] aplicó un filtrado mediante umbralización combinado con la correlación en el dominio tiempo-frecuencia para remover el ruido obteniendo resultados satisfactorios. Por otra parte se ha aplicado la teoría del filtrado Wiener a través de la Transformada de Ondícula. Los trabajos de Ghael [7] y Hyeokho [8] apuntan en esta dirección. Estos se basan en la atenuación de los coeficientes de aproximación y de detalle de una señal ruidosa a través de un filtro Wiener caracterizado por las densidades espectrales de potencia de la señal deseada y del ruido en el dominio de la Transformada de Ondícula. 

El trabajo que a continuación se presenta tiene como objetivo reducir el ruido presente en una señal sísmica mediante filtrado en el dominio de la Transformada de Ondícula. En esta investigación se propone combinar la técnica de umbralización de los coeficientes de la Transformada de Ondícula y el filtrado Wiener en el dominio de esta transformada. Adicionalmente se comparan los resultados del esquema de filtrado propuesto con el clásico filtrado Wiener en el dominio de la frecuencia. 

Filtrado Wiener 

Dado un proceso estacionario x(t)=s(t)+r(t) donde s(t) es la señal original no distorsionada por perturbaciones, r(t) es el proceso aleatorio que define el ruido estadísticamente y x(t) la señal contaminada con ruido, se puede estimar la señal s(t) usando un filtro Wiener. Si s(t) y r(t) son señales estadísticamente estacionarias, con media cero y no están correlacionadas, el filtro Wiener óptimo está definido en el dominio de la frecuencia por la siguiente ecuación [9]:

Tomando en cuenta que la relación señal a ruido RSR es el cociente entre la densidad espectral de la señal deseada y la densidad espectral del ruido se puede expresar el filtro Wiener en función de la relación RSR obteniéndose:

En la práctica se puede estimar y conocida la densidad espectral de potencia de la señal x(t), Px(w). Suponiendo que las señales s(t) y r(t) no están correlacionadas se puede demostrar que [9]: 

De la ecuación (3), si se tiene caracterizado el ruido estadísticamente, se puede estimar la densidad espectral de potencia mediante la diferencia entre  y . 

Transformada de Ondícula 

Si una señal f(t) es analógica con energía finita, entonces la Transformada de Fourier, , representa el espectro de la señal y se define como:

La Transformada de Fourier no es útil cuando la señal carece de regularidad estadística determinada o cuando tiene variaciones de frecuencia locales, ya que no contempla variaciones de la frecuencia con el tiempo [10]. Para solventar este problema se aplica la Transformada de Fourier de Tiempo Corto (STFT) que consiste en dividir la señal en pequeñas porciones contiguas, de tal manera que se pueda suponer en cada segmento la señal estacionaria y así calcular la Transformada de Fourier de cada porción. La Transformada de Fourier de Tiempo Corto se define [11]: 

El análisis a través de la STFT implica trabajar con una resolución fija en el espacio tiempo frecuencia. Esto puede representar un problema para señales que requieren buena resolución en ambos dominios ya que, según el Teorema de Incertidumbre de Heisenberg, al aumentar la resolución en tiempo disminuye en frecuencia y viceversa. Una manera de solventar esta limitante es a través del análisis multiresolución provisto por la Transformada de Ondícula. 

La Transformada de Ondícula se ha convertido en la última década en una herramienta común en aplicaciones que involucran análisis de señales. Dada una función variante con el tiempo f(t) la transformada de ondícula puede ser vista como el cálculo de los coeficientes del producto interno de la señal con las funciones basadas en familias de ondículas. Por definición las familias de funciones ondículas se originan al escalar y desplazar en el tiempo una ondícula madre o prototipo ψ(t). En este sentido la transformada de ondícula continua se define como [11]: 

La familia de ondículas es:

Donde ψ(t) es la función ondícula que tiene características de función pasabanda. El factor 1/a^1/2 tiene como finalidad preservar o mantener la energía de la señal. Si se eligen para el escalamiento y el desplazamiento en tiempo los valores a=a₀^m y b=na₀^mb₀ con m y n valores enteros y se aproxima la integral a través de una sumatoria se obtiene la fórmula de la Transformada de Ondícula Discreta.

Discretizando el tiempo t=kTd, con Td igual al período de discretización o muestreo, y asignando los valores a₀=2 y b₀=1 se obtiene, además de una computación eficiente, una representación compacta de la Transformada de Ondícula Discreta.

A través del análisis multiresolución se pueden calcular los coeficientes que definen a la DWT a baja frecuencia mediante el filtrado de la señal con un filtro pasabajo h(n) definido a partir de la función escalamiento Ф(t) asociada a la ondícula ψ(t). Los coeficientes correspondientes a la alta frecuencia se obtienen procesando la señal con un filtro pasabanda g(n) definido a partir de la función ondícula ψ(t). Este proceso se indica en las siguientes ecuaciones [12]: 

En las ecuaciones anteriores la DWT correspondiente a las bajas frecuencias se expresa a través de los coeficientes a_j+1 y la DWT correspondiente a las altas frecuencia se expresa a través de los coeficientes d_j+1. El índice j se asocia con el nivel de descomposición ondícula y corresponde al escalamiento 2^j y a la aproximación de la señal a la resolución 2^-j. 

Dado que la representación de una señal a través de los coeficientes DWT es ortogonal se puede recuperar la señal original a partir de los coeficientes calculados. Se puede demostrar que la recuperación de la señal puede obtenerse a partir de la siguiente ecuación [12]:

Filtrado en el dominio de la transformada de ondícula 

El método propuesto para suprimir el ruido de una señal sísmica combina la técnica de umbralización de los coeficientes de la Transformada de Ondícula y el filtrado Wiener en el dominio de esta transformada. 

Filtrado mediante umbralización de los coeficientes ondícula 

La técnica de umbralización se basa en la atenuación de los coeficientes de detalle de la Transformada de Ondícula según una regla o ley. En general para una señal unidimensional en tiempo discreto las altas frecuencias afectan los coeficientes de detalle de los primeros niveles mientras que, las bajas frecuencias, afectan los coeficientes de detalle de los niveles más altos y los coeficientes de las aproximaciones. Si se analiza una señal de ruido blanco gaussiano en el dominio de la Transformada de Ondícula se encontrará que los coeficientes de detalles en niveles sucesivos decrecen en magnitud al incrementarse el nivel de análisis. 

Para la supresión del ruido presente en una señal mediante umbralización se deben seguir los siguientes etapas [13]: 

• Descomposición: En esta etapa se calcula la descomposición ondícula de la señal ruidosa hasta un nivel N con una ondícula madre específica. 

• Umbralización de los coeficientes de detalle: En esta etapa se modifican los coeficientes de detalle de la descomposición ondícula desde el nivel 1 hasta el N aplicando alguna regla de umbralización. 

• Reconstrucción: En esta etapa se debe reconstruir la señal mediante la aplicación de la Transformada de Ondícula Inversa a partir de los coeficientes de aproximación originales y los coeficientes de de detalle modificados en el paso anterior. 

Existen dos tipos de umbralización: La dura, en donde se asigna el valor cero a aquellos coeficientes cuyo valor absoluto este por debajo del umbral, y la suave, en donde se descartan aquellos coeficientes cuyo valor absoluto esta por debajo del umbral comprimiéndose, según una regla, los coeficientes restantes [13]. Las reglas de umbralización que se usan comúnmente en la supresión del ruido son: La función de pérdidas cuadráticas (QLF) que selecciona el umbral de manera adaptativa usando el principio de estimación de riesgo de Stein, la función , la función SURE que es una variante heurística de la función de pérdidas cuadráticas y La función umbral MINIMAX [13]. 

Filtrado Wiener en el dominio de la transformada de ondícula 

El ruido presente en una señal puede modelarse aditivo. Si x(i) representa la iésima muestra de la señal ruidosa esta puede ser escrita como x(i)=s(i)+r(i) en donde s(i) y r(i) representan la iésima muestra de la señal original y el ruido respectivamente. Si aplicamos la Transformada de Ondícula Discreta (DWT) a las señales x y s denotando DWT_x(i), DWT_s(i) los iésimos coeficientes de la DWT de las señales respectivas y suponemos que el ruido es blanco gaussiano con desviación estándar s, el filtro Wiener en el dominio de la Transformada de Ondícula puede expresarse como [7]: 

Al igual que el filtro Wiener en el dominio de la frecuencia, el filtro Wiener en el dominio ondícula necesita estimaciones de DWT_s²(i) y s². Si se tiene caracterizado el ruido estadísticamente mediante la desviación estándar s y se supone que las señales aleatorios s y r no están correlacionadas se puede estimar DWT_s²(i) mediante la ecuación [7]: 

Metodología 

En este trabajo se presenta una técnica que incrementa desempeño del filtro Wiener en el dominio de la frecuencia mediante la aplicación de un esquema mixto basado en la técnica de umbralización y el filtrado Wiener en el dominio de la Transformada de Ondícula. La señal sísmica a ser analizada será modelada sintéticamente con la ondícula de Ricker y la perturbación o ruido agregado se modelará sintéticamente como un proceso aleatorio blanco gaussiano de media cero y desviación estándar σ, de esta manera la densidad espectral de potencia asociada al ruido será directamente σ². Para evaluar el desempeño del esquema de filtrado propuesto se compararán los resultados con el filtrado Wiener clásico en el dominio de la frecuencia. La señal sísmica se muestrea con una frecuencia de 2500 Hz y el ruido agregado se define según los valores de desviación estándar 0,025; 0,125; 0,225 y 0,325. Estos valores se definieron tomando en cuenta que la desviación estándar de la ondícula de Ricker es 0,2236. Para las pruebas se diseño un programa en la Plataforma Matlab que modela los procesos sísmicos y el filtrado. 

Para definir los parámetros del filtrado por umbralización de los coeficientes ondícula se fijó la descomposición al primer nivel, como lo sugiere Ghael [7] dado que el ruido se asocia a los primeros coeficientes de detalle. El tipo de umbralización se varió entre duro y suave y la regla de umbralización entre QLF, SQRLOG, SURE y MINIMAX. De esta manera evaluando el error respecto a la señal ideal, se escogerá aquella combinación de parámetros en donde se presente el mínimo error. De la misma manera se varió en el análisis la ondícula madre entre Daubechies y Symlet con momentos de desvanecimiento 1, 3, 6 y 9. Dado que estas ondículas son de soporte compacto su incorporación en el análisis permitirá evaluar resultados con el mínimo número de coeficientes con la variación de los momentos de desvanecimiento de las mismas. 

Una vez escogido los mejores parámetros del filtrado por umbralización se deben definir los parámetros para el filtrado Wiener en el dominio de la Transformada de Ondícula. En este sentido, por las razones expuestas en la selección de los parámetros de la técnica de umbralización, la descomposición de la señal se realizará con las mismas ondículas madres Daubechies y Symlet para el primer nivel de descomposición y para los mismos momentos de desvanecimiento 1, 3, 6 y 9. Una vez encontrados los mejores parámetros para el filtrado Wiener en el dominio ondícula se compararán los resultados obtenidos con los del filtrado Wiener en el dominio de la frecuencia. 

Análisis de Resultados 

Las pruebas para la determinación del mejor filtro de umbralización se realizaron variando los parámetros de ajuste: tipo de umbralización, regla de umbralización y ondícula madre. El mejor filtro será aquel que presente menor error entre la señal filtrada y la señal sísmica original. El ruido añadido tiene una desviación estándar igual a la desviación estándar de la ondícula de Ricker (σ= 0,2236). Para el cálculo del error se realizó la sumatoria del valor absoluto de la diferencia entre los puntos del muestreo de la respuesta del filtro y la señal original. Estos resultados se muestran en las Tablas 1 y 2. 

de los resultados obtenidos podemos comentar que el menor error si se usa la ondícula de Daubechies se obtuvo a través del tipo de umbralización suave con la regla MINIMAX para la ondícula madre Daubechies1. Para el caso de la ondícula de Symlet se obtuvo el menor error usando el mismo esquema de umbralización (tipo: suave, regla:MINIMAX) pero a través de la ondicula madre Symlet1. Estos errores fueron respectivamente 0,2813 y 0,2451.

Para el ajuste de los parámetros del filtro Wiener en el dominio ondícula se procesó la respuesta del filtrado por umbralización del esquema con mínimo error (tipo: suave, regla: MINIMAX, ondícula, ondícula madre: Db1, Db3, Db9, Sym1, Sym3, Sym9) variando la desviación estándar del ruido y la ondícula madre en el análisis. Los mínimos errores se encontraron para las ondículas madres Daubechies6 y Symlet6 respectivamente. En las Tablas 3 y 4 se presentan los errores obtenidos al procesar la señal con el filtro Wiener en el dominio de la frecuencia, el filtrado por umbralización ondícula y el filtrado a través del esquema combinado umbralización-Wiener ondícula. Para efectos de comparación, también se presenta el error inicial entre la señal ruidosa y la señal sísmica original.

De los resultados obtenidos en la Tabla 3 se puede decir que el mejor esquema combinado se obtiene con la ondícula Db1 en la umbralización y la ondícula Db6 en el filtrado Wiener. Si se usa Db1 en la umbralización se reduce el error entre 6% y 37% a medida que la potencia del ruido disminuye. Por otra parte, cuando se usa esta combinación (umbralización Db1, Wiener Db6) la reducción del error inicial varía del 25 al 60% con la disminución de la potencia del ruido. En contraste, la reducción del error usando el filtro Wiener en el dominio de la frecuencia varía de 0,3% al 28% con misma disminución de la potencia del ruido. 

De la Tabla 4 se puede concluir que el mejor esquema combinado se obtiene usando la ondícula de Symlet1 en la umbralización y la ondícula Sym 6 en el filtrado Wiener. Si se usa Sym1 en la umbralización se reduce el error entre 1% y 27% a medida que la potencia del ruido disminuye. Cuando se usa esta combinación la reducción del error inicial varía del 6 al 54% con la disminución de la potencia del ruido. En contraste, la reducción del error usando el filtro Wiener en el dominio de la frecuencia varía de 0,3% al 28% con misma disminución de la potencia del ruido.

Conclusiones

Del análisis de resultados se puede concluir que: 

• La técnica de filtrado por umbralización en el dominio ondícula mejora la recuperación de la señal Ricker, contaminada con ruido blanco gaussiano aditivo de desviación estándar 0,025; 0,125; 0,225 y 0.325, respecto al filtrado Wiener en el dominio de la frecuencia. 

• La técnica de filtrado mediante un esquema combinado de umbralización y filtrado Wiener en el dominio ondícula mejora la recuperación de la señal Ricker, contaminada con ruido blanco gaussiano aditivo de desviación estándar 0,025; 0,125; 0,225 y 0,325, respecto al filtrado Wiener en el dominio de la frecuencia.

• El mejor filtrado usando el esquema combinado de umbralización y filtrado con la ondícula de Daubechies se obtuvo para las ondículas madres Daubechies1 y Daubechies6 respectivamente. 

• El mejor filtrado usando el esquema combinado de umbralización y filtrado con la ondícula Symlet se obtuvo para las ondículas madres Symlet1 y Symlet6 respectivamente. 

Referencias Bibliográficas 

1.  Yilmaz, O. “Seismic Data Processing. Investigation of Geophysics”, vol 2, 1997.          [ Links ]

2. Foufoula, Georgiou and Praveen, K. “Wavelet Analysis for Geophysical applications”. Reviews of Geophysics, 35 (4): (November, 1997).          [ Links ]

3.  Charneley A. y Alexander N.: “An approach to seismic correction which includes wavelet de-noising”. Proceedings of the sixth Conference on Computational Structures Technology, (2002), 107-108.          [ Links ]

4. Zhou Y y Garossino P. “High-energy noise attenuation of seismic data in the wavelet-transform domain”. Integrated Computer-Aided Engineering –IOS Press., 12 (1): (2005). 57-67.          [ Links ]

5. Siyuan C y Xiangpeng Ch. “The second-generation wavelet transform and its application in denoising of seismic data”. Applied Geophysics, Chinese Geophysical Society, 2 (2): (2005). 70-74.          [ Links ]

6. Yan F y Chunguin Z. “A new seismic data de-noising method based on wavelet”. Proceedings of the International Conference on Computer Science and Information Technology, ICCSIT (2008), 92-96.          [ Links ]

7.  Ghael S., Sayeed A. y Baraniuk R.: “Improved Wavelet Denoising via Empirical Wiener Filtering”. Proceedings of SPIE (1997), 17-27.          [ Links ]

8.  Hyeokho Ch. y Baraniuk R.: “Analysis of Wavelet-Domain Wiener Filters”. Proceedings of the IEEE-SP International Symposium on Time-frequency and Time-scale Analysis (1998), 6-9.          [ Links ]

9.  Wiener, N.: “Selected papers of Norbert Wiener”, published by Society for Industrial and Applied Mathematics and The Technology, Cambridge Massachusetts. 1964.          [ Links ]

10.  Bracewell R. “The Fourier Transform and its applications”, McGraw Hill Kogakusha, LTD., second edition, 1978.          [ Links ]

11.  Chui Ch. “An introduction to wavelets, Wavelets Analysis and its applications”, Academic Press, Texas. 1998          [ Links ]

12.  Mallat S.: “Wavelet Tour of Signal Processing”, Academic Pres, San Diego, 1998.          [ Links ]

13.  Misiti M., Misiti Y., Oppenheim G., Poggi J.: ”Wavelet Toolbox: For Use with Matlab”, The MathWorks, Inc., Mass, 1996.        [ Links ]