Pedro Mário Cruz e Silva Visualização Volumétrica de ...webserver2.tecgraf.puc-rio.br/~mgattass/teses/2004TesePedroMario.pdf · Todos os direitos reservados. É proibida a reprodução

Pedro Mário Cruz e Silva

Visualização Volumétrica de Horizontes em Dados Sísmicos 3D

Tese de Doutorado

Tese apresentada como requisito parcial para obtenção do título de Doutor pelo Programa de Pós-Graduação em Informática da PUC-Rio.

Orientadores: Marcelo Gattass Paulo Cezar Carvalho

Rio de Janeiro, junho de 2004

Todos os direitos reservados. É proibida a reprodução total ou parcial do trabalho sem autorização da universidade, do autor e do orientador.

Pedro Mário Cruz e Silva Pedro Mário Cruz e Silva graduou-se bacharel em matemática na UFPE (Universidade Federal de Pernambuco) em 1995. Concluiu o mestrado em matemática também na UFPE em 1998. Participa de projetos de pesquisa no Tecgraf/PUC-Rio (Laboratório de computação gráfica) desde 2000. Atualmente participa do desenvolvimento de um sistema de visualização de dados sísmicos realizado em parceria com a gerência de Tecnologia em Geofísica da Petrobras.

Ficha Catalográfica

Silva, Pedro Mário Cruz

Descrição.

Inf. Técnicas

Natureza.

Incluí referências bibliográficas.

Visualização volumétrica; Dados sísmicos 3D; Rastreamento de Horizontes; Orientação local; Atributos Sísmicos Instantâneos.

À minha família.

Agradecimentos

Primeiramente a Deus pelo caminho que me foi reservado até aqui.

Aos meus pais Mário e Maria da Saúde pela educação recebida, pelo

carinho, pelos cuidados e pelo amparo de sempre. Ao meu irmão Mário Sérgio,

tios, primos e toda minha família pelo ótimo ambiente no qual cresci e vivo até

hoje.

À minha futura esposa Patrícia agradecimentos mais que especiais pelo

incentivo, carinho, compreensão, apoio e companheirismo.

Ao meu orientador Professor Marcelo Gattass por toda confiança

depositada em mim, toda paciência, todo suporte acadêmico e principalmente

pela oportunidade de participar deste trabalho de pesquisa. Ao meu co-

orientador Professor Paulo Cezar pelas sugestões e comentários sempre

precisos e relevantes.

A todos os professores e colegas do Departamento de Informática da

PUC-Rio e do IMPA pelo ótimo ambiente acadêmico.

A todos os amigos do laboratório Tecgraf, principalmente à Carolina

Alfaro pelas revisões feitas no texto da tese.

A todos os amigos da Petrobras, em especial ao Marcos Machado pela

parceria no desenvolvimento desta pesquisa.

A todos os amigos e amigas que de alguma forma contribuíram para a

realização deste trabalho.

Ao CNPq e à PUC-Rio pelos auxílios concedidos sem os quais este

trabalho não poderia ter sido realizado.

Resumo Silva, Pedro. Visualização Volumétrica de Horizontes em Dados Sísmicos 3D. Rio de Janeiro, 2004. 104p. Tese de Doutorado - Departamento de Informática, Pontifícia Universidade Católica do Rio de Janeiro.

Neste trabalho apresentamos os aspectos da visualização volumétrica de

horizontes em dados sísmicos 3D. Consideramos as abordagens de visualização

volumétrica direta e indireta. Na abordagem direta investigamos o problema da

seleção de horizontes usando a função de transferência. Apresentamos a

técnica de opacidade 2D que busca aumentar a habilidade de seleção dos

horizontes para visualização. Comparamos a utilização dos atributos de fase

instantânea, fase ajustada e fase desenrolada como segunda dimensão,

enquanto a primeira é a amplitude sísmica. Ainda na abordagem direta,

mostramos que o gradiente da amplitude sísmica não aproxima bem os vetores

normais nos horizontes sísmicos. Sugerimos o gradiente da fase instantânea

como solução para este problema. Na abordagem de visualização volumétrica

indireta introduzimos uma modelagem de otimização para o problema de

rastreamento de horizontes. Sugerimos um método heurístico baseado em uma

estratégia gulosa para encontrar soluções que são boas aproximações para os

horizontes mesmo na presença de estruturas geológicas complexas.

Palavras-chave Visualização volumétrica; Dados sísmicos 3D; Rastreamento de

Horizontes; Orientação local; Iluminação de Horizontes; Atributos Sísmicos

Instantâneos.

Abstract Silva, Pedro. Volume Visualization of Horizons in 3-D Seismic Data. Rio de Janeiro, 2004. 104p. Tese de Doutorado - Departamento de Informática, Pontifícia Universidade Católica do Rio de Janeiro.

This work presents the aspects of volume visualization of seismic horizons

in 3-D seismic data. We consider both the direct and indirect approach of volume

visualization. In the direct approach we investigate the problem of selecting

horizons using the transfer function. We present the 2-D opacity technique which

tries to increase ability to select horizons for visualization. We compare the use of

instantaneous phase, adjusted phase and unwrapped phase as the second

dimension, while seismic amplitude is the first dimension. Also in the direct

approach, we show that the seismic amplitude gradient is not a good

approximation for the normal vectors in seismic horizons. We suggest the

gradient of instantaneous phase as a solution to this problem. In the indirect

volume visualization approach we introduce a new optimization model for the

seismic horizon tracking problem. We present a heuristic method based on a

greedy strategy for finding solutions that are good approximations for the horizon

of interest, even for complex geological structures.

Keywords Volume Rendering; 3-D Seismic data; Seismic Horizon Tracking; Local

Orientation; Seismic Horizon Illumination; Instantaneous Seismic Attributes.

Sumário

1 Introdução 14 1.1. Objetivos 14 1.2. Estrutura do Documento da Tese 15

2 Conceitos Básicos 17 2.1. Dados Sísmicos 3D 17 2.1.1. Aquisição 18 2.1.2. Processamento 19 2.1.3. Interpretação 21 2.2. Visualização Volumétrica direta com texturas 3D 21 2.2.1. Amostragem da textura 23 2.2.2. Composição 23 2.2.3. Freqüência de amostragem 24 2.2.4. Gerenciando a memória de textura 24 2.2.5. Iluminação 25 2.2.6. Função de Transferência 28 2.2.7. Algoritmo 29 2.2.8. Funções de transferência multidimensionais 33

3 Opacidade 2D 39 3.1. Motivação 39 3.2. Trabalhos Relacionados 42 3.3. Soluções Propostas 43 3.3.1. Amplitude Sísmica X Fase Instantânea 44 3.3.2. Amplitude Sísmica X Fase Ajustada 49 3.3.3. Amplitude Sísmica X Fase Desenrolada 55 3.4. Conclusões Parciais 62

4 Iluminação de Dados Sísmicos 63 4.1. Motivação 63 4.2. Trabalhos Relacionados 70 4.3. Solução Proposta 70

4.4. Resultados 72 4.5. Conclusões Parciais 77

5 Rastreamento de Horizontes Sísmicos 79 5.1. Trabalhos Relacionados 80 5.2. Modelagem de Otimização 80 5.2.1. Volumétrico Sísmico Discreto 80 5.2.2. Conjunto de traços do dado 81 5.2.3. Horizonte Sísmico Discreto 82 5.2.4. Conjunto de Sementes 82 5.2.5. Função de Similaridade 83 5.2.6. Função de Semente Associada 84 5.2.7. Problema de Otimização 84 5.3. Método Proposto 85 5.3.1. Inicialização 86 5.3.2. Aprimoramento da Solução Corrente 86 5.3.3. Convergência 88 5.3.4. Análise de Complexidade 88 5.4. Resultados 90 5.5. Conclusões Parciais 95

6 Conclusões e Trabalhos Futuros 96 6.1. Conclusões 96 6.2. Trabalhos Futuros 97

7 Bibliografia 98

Apêndice A: Transformada de Hilbert 102

Lista de Símbolos

Capítulo 3

( )X t Traço sísmico. Função de amplitude sísmica medida na direção

vertical (tempo).

( )Y t Transformada de Hilbert de um traço sísmico.

( )h t Núcleo da transformada de Hilbert.

∗ Sinal de convolução.

i Identidade complexa. 1i = − .

( )Z t Traço analítico. ( ) ( ) ( )Z t X t iY t= + ∈ .

( )A t Amplitude instantânea (envelope).

( )tθ Fase instantânea.

( )tθ ′ Fase ajustada.

Capítulo 4

),,( yxtX

Campo escalar de amplitudes sísmicas.

),,( yxtY Campo escalar das transformadas de Hilbert dos traços sísmicos.

),,( yxtθ Campo escalar de fase instantânea.

)(tf Freqüência instantânea.

n Vetor normal.

),,( BGR Representação de cor RGB.

Capítulo 5

Conjunto dos números inteiros.

3 Espaço discreto { }( , , ) | , ,i j k i j k ∈ .

U Subconjunto de 3Z onde o dado sísmico está amostrado. 3[0, 1] [0, 1] [0, 1]U NT NX NY= − × − × − ⊂ Z

U~ Espaço discreto que representa o conjunto de traços do dado.

[0, 1] [0, 1]U NX NY= − × − .

( , , )X i j k Valor da amplitude sísmica amostrada em um ponto de U .

:X U → R

NT Número de amostras em cada traço do dado.

NX Número de crosslines do dado.

NY Número de inlines do dado.

S Conjunto de voxels de U escolhidos como sementes.

v Traço ( , )v x y U= ∈ que contém o voxel ( , , )v t x y U= ∈ .

σ Função de similaridade entre voxels. :U Uσ × → R .

H Conjunto de voxels contidos num horizonte sísmico. Pode ser

interpretado como um par ( , )H W τ= .

H + Conjunto de voxels candidatos a entrar na solução H .

s Função que associa uma semente a cada voxel de um horizonte.

:s H S→ .

Introdução 11

Lista de figuras

Figura 1 – Etapas da exploração baseada na sísmica de reflexão. 17 Figura 2 – Aquisição sísmica. Adaptada de (Gerhardt, 1998). 18 Figura 3 - Traço sísmico (esquerda), Linha sísmica (centro) e volume sísmico

(direita). 20 Figura 4 – Modelo de convolução. Adaptado de (Gerhardt, 1998). 20 Figura 5 – Modelo geológico. Adaptada de (Robinson e Treitel, 1980). 21 Figura 6 – Geometrias de amostragem. Adaptada de (Blythe e McReynolds,

2000). 23 Figura 7 – Volume de dados dividido em tijolos. 25 Figura 8 – Iluminação de uma superfície de nível de um volume. Adaptado de

(Hadwiger et al, 2002) 26 Figura 9 – Gradiente de uma superfície de nível. Adaptado de (Westermann e

Ertl, 1998). 27 Figura 10 – Voxels com densidade baixa (esquerda), voxels com densidade alta

(centro), voxels com densidade baixa e alta (direita). 28 Figura 11 – Mapeamento de dados da memória principal para a placa de vídeo.31 Figura 12 – Cálculo de iluminação feito por pixel durante a visualização. 33 Figura 13 – Função de Transferência Bidimensional. Adaptada de (Pfister et al,

2001) 34 Figura 14 – Realce de fronteiras. Adaptada de (Pfister et al, 2001) 35 Figura 15 – Mapeamento de dados. Função de transferência 2D. 36 Figura 16 – Função de transferência 2D. 38 Figura 17 - Superposição dos intervalos de amplitudes dos eventos sísmicos 40 Figura 18 - Dado sintético simulando uma estrutura de domo (anticlinal). 41 Figura 19 - Efeito da superposição de intervalos de amplitude na visualização

volumétrica de dados sísmicos. 42 Figura 20 - Amplitude sísmica e fase instantânea de um traço sísmico real 45 Figura 21 - Função de transferência 2D (esquerda) e visualização volumétrica do

dado do Domo (direita). 46 Figura 22 - Evento de interesse isolado graças a opacidade 2D. Função de

transferência (esquerda) e visualização volumétrica do evento (direita). 47 Figura 23 - Volume de dados real (esquerda). Eventos sísmicos reais

destacados (direita). 48

Introdução 12

Figura 24 - Função de opacidade 1D aplicada a um dado sísmico real (acima).

Função de opacidade 2D aplicada a um dado sísmico real (abaixo). 49 Figura 25 - Artefatos na visualização volumétrica da Fase igual a zero. 51 Figura 26 - Função de ajuste de fase. Expressão (esquerda) e gráfico (direita). 51 Figura 27 - Histogramas 2D de Amplitude X Fase Instantânea (esquerda) e

Amplitude X Fase Ajustada (direita). 52 Figura 28 – Comparação entre Fase Instantânea e Fase Ajustada. 53 Figura 29 – Visualização volumétrica dos eventos positivos e negativos

separadamente usando a fase ajustada. 54 Figura 30 – Fase desenrolada em 1D. 56 Figura 31 – Fase desenrolada em 2D. Dado sintético. 58 Figura 32 – Fase desenrolada em 2D. Dado real. 59 Figura 33 – Comparação da opacidade 2D com fase instantânea e fase

desenrolada. Dado sintético. 60 Figura 34 – Comparação da opacidade 2D com fase instantânea e fase

desenrolada. Dado real. 61 Figura 35 - Traço gerador do modelo domo com os eventos destacados. 64 Figura 36 – Visualização volumétrica do modelo plano. 65 Figura 37 – Visualização volumétrica direta com iluminação aplicada ao modelo

plano. 66 Figura 38 – Modelo domo modificado. Modelo sintético com variação de

amplitudes ao longo dos refletores. 68 Figura 39 – Iluminação aplicada a um horizonte com variação lateral de

amplitudes. 69 Figura 40 – Superfícies de nível do modelo sintético do domo modificado. 69 Figura 41 – Imagem do modelo plano obtida pelo método proposto (esquerda).

Imagem obtida pelo método tradicional (direita). 73 Figura 42 – Gráfico da amplitude sísmica, fase instantânea e freqüência

instantânea do traço gerador do modelo plano. 73 Figura 43 – Visualização do campo vetorial normal. Campo normal à superfície

de uma esfera (esquerda). Campo normal a um horizonte sísmico obtido pela

solução proposta (centro), campo obtido pelo método tradicional (direita). 75 Figura 44 – Visualização volumétrica do dado sintético com variação de

amplitudes ao longo dos horizontes: sem iluminação (equerda), iluminado

usando a solução proposta (centro) e iluminado da maneira tradicional (direita).76 Figura 45 – Visualização volumétrica de um dado sísmico 3D real: sem

Introdução 13

iluminação (esquerda), iluminação usando o método proposto (centro) e

iluminação usando a abordagem tradicional (direita). 76 Figura 46 – Visualização volumétrica. Adaptado de (Silva et al, 2003). 77 Figura 47 – Horizonte sísmico em um dado real. 90 Figura 48 – Evolução do algoritmo. 91 Figura 49 – Sementes e função de semente associada. 92 Figura 50 – Estrutura temporal do horizonte. 92 Figura 51 – Horizontes diferentes alinhados por uma falha geológica. 93 Figura 52 – Horizonte “A” rastreado usando uma ferramenta comercial. 94 Figura 53 – Horizonte “A” rastreado usando nosso método de rastreamento. 95

Introdução 14

1 Introdução

1.1. Objetivos

O objetivo principal deste trabalho é analisar a aplicação das técnicas de

visualização volumétrica a dados sísmicos 3D. A visualização volumétrica é um

processo de síntese de imagens que trabalha com modelos descritos

implicitamente. Estes modelos, chamados de dados volumétricos ou

simplesmente volumes, apresentam-se como amostras de um campo escalar

obtidas em um grid regular. Os dados obtidos por tomografia computadorizada e

ressonância magnética são exemplos de volumes e são utilizados principalmente

em medicina, dados volumétricos também são originados como resultado de

simulações de fenômenos físicos. Para os dados desta natureza os algoritmos

de visualização volumétrica conhecidos hoje em dia são capazes de produzir

imagens bastante informativas e com alta qualidade visual.

Os dados sísmicos volumétricos são modelos implícitos da subsuperficie

que integram as características estruturais e estratigráficas das subcamadas. A

industria de óleo e gás vem utilizando a visualização volumétrica em volumes

sísmicos na área de exploração, tanto para controle de qualidade quanto para

auxilio a interpretação, pois esta permite a obtenção de uma visão prévia das

feições geológicas contidas no dado de maneira relativamente simples e

interativa.

Apesar dos dados sísmicos volumétricos se adequarem ao modelo

esperado pelos algoritmos de visualização volumétrica, isto é, serem um

conjunto de amostras obtidas em cima do campo escalar de amplitude sísmica

em uma grade regular, vários problemas ocorrem ao se tentar aplicar

diretamente as técnicas de visualização volumétrica aos dados de natureza

sísmica.

Neste trabalho consideramos problemas que ocorrem nas etapas de

classificação e iluminação que integram a maioria dos algoritmos de visualização

volumétrica direta. A etapa de classificação consiste na aplicação da função de

transferência ao volume, é a função de transferência que determina qual porção

Introdução 15

do dado estará visível na imagem final e defina as propriedades óticas dessa

porção, como por exemplo, a cor e a opacidade. Para selecionar regiões de

interesse em dados sísmicos, que são os eventos sísmicos, sugerimos a

utilização de funções de transferência multidimensionais. Analisando a natureza

oscilatória do dado sísmico propomos funções de transferência bidimensionais,

definidas nas dimensões de amplitude (dimensão original) e fase instantânea

(atributo sísmico calculado a partir da amplitude).

Na etapa de iluminação o gradiente do volume freqüentemente é usado

como vetor normal à superfície de interesse para o cálculo de iluminação local. É

um dos resultados deste trabalho a verificação de que esta abordagem não se

aplica a dados sísmicos novamente devido a sua natureza oscilatória. Propomos

aqui a utilização do gradiente do atributo de fase instantânea como estimador da

orientação local dos horizontes sísmicos e validamos esta técnica em volumes

sísmicos sintéticos e reais.

O rastreamento de horizontes pode ser considerado, do ponto de vista de

visão computacional, como uma técnica de segmentação ou, do ponto de vista

de modelagem e visualização, como uma técnica de extração de superfícies. A

finalidade do rastreamento vai além da simples visualização do horizonte. Mas,

dentro do contexto deste trabalho de tese, consideramos o rastreamento de

horizontes como sendo uma técnica de visualização volumétrica direta adaptada

para dados sísmicos. Apresentamos neste trabalho uma abordagem de

otimização para o rastreamento de horizontes como uma opção de visualização

volumétrica indireta de horizontes sísmicos.

1.2. Estrutura do Documento da Tese

No capítulo 2 fazemos uma pequena revisão de alguns conceitos básicos

que são importantes para a compreensão do restante do documento. Este

capítulo está dividido basicamente em duas partes. Na primeira parte

apresentamos os dados sísmicos 3D, sua origem e características básicas.

Também definimos o que são horizontes sísmicos que são o principal objeto de

estudo deste trabalho. Esta parte é destinada principalmente aos pesquisadores

da área de computação gráfica que não estão familiarizados com os dados desta

natureza. Na segunda parte descrevemos com um certo nível de detalhe o

método de visualização volumétrica baseada em texturas 3D e a técnica de

funções de transferência multidimensionais. Julgamos esta parte importante pois

Introdução 16

nela está descrita a implementação do sistema utilizado para testar as idéias dos

capítulos 3 e 4. Esta parte é destinada principalmente aos pesquisadores da

área de geologia e geofísica que não estão familiarizados com esta técnica de

computação gráfica.

No capítulo 3 apresentamos a técnica de opacidade 2D. A técnica de

opacidade 2D é uma aplicação de funções de transferência multidimensionais à

dados sísmicos. Esta aplicação foi publicada por nós em um trabalho anterior.

Neste capítulo sugerimos o uso de atributos sísmicos instantâneos juntamente

com a amplitude sísmica na opacidade 2D. Fazemos uma comparação dos

resultados obtidos com o atributo de fase instantânea, fase ajustada e fase

desenrolada. O critério de comparação entre os atributos é a capacidade de se

isolar horizontes sísmicos usando o atributo na opacidade 2D.

No capítulo 4 analisamos a etapa de iluminação dos algoritmos de

visualização volumétrica direta, apresentamos os problemas que ocorrem nesta

etapa ao adotar o gradiente do campo escalar como vetor normal no modelo de

iluminação local. Descrevemos a utilização do gradiente do atributo de fase

instantânea como solução para este problema. Apresentamos os resultados

obtidos em dados sísmicos sintéticos e reais.

No capítulo 5 sugerimos uma modelagem de otimização para o problema

do rastreamento de horizontes sísmicos. Em seguida propomos um método

heurístico baseado em uma estratégia gulosa para encontrar soluções quase-

ótimas para o problema. Mostramos através de exemplos que as soluções

quase-ótimas obtidas são boas aproximações para os horizontes sísmicos.

No capítulo 6 faremos as considerações finais sobre o trabalho de

pesquisa desenvolvido e apresentamos a pesquisa em andamento que poderá

resultar em outros trabalhos no futuro.

O Apêndice A contém um resumo da teoria da Transformada de Hilbert

que é usada no cálculo dos atributos sísmicos instantâneos.

Conceitos Básicos 17

2 Conceitos Básicos

O caráter interdisciplinar deste trabalho torna necessária a apresentação

de conceitos básicos da área de exploração sísmica e de computação gráfica.

Na seção seguinte faremos uma breve descrição do processo de obtenção de

dados sísmicos 3D. Em seguida introduziremos os atributos sísmicos

instantâneos têm um papel importante no nosso estudo.

As discussões e resultados relativos à visualização volumétrica direta de

dados sísmicos contidos neste documento (capítulos 3 e 4) são independentes

de implementação, ou seja, são pertinentes as implementações de sistemas de

visualização volumétrica de dados sísmicos em geral. Porém, para validar

nossas idéias desenvolvemos um sistema de visualização volumétrica baseado

na utilização de texturas 3D. Escolhemos esta opção por esta permitir um bom

nível de interatividade com o modelo. Atualmente existe uma vasta literatura

sobre o uso de textura 3D na visualização volumétrica direta. Na seção 2.2.1

faremos um apanhado das técnicas usadas na nossa implementação.

2.1. Dados Sísmicos 3D

A sísmica de reflexão é um método indireto de exploração da subsuperfície

da terra. Este método vem sendo largamente utilizado pelo fato de ser capaz de

cobrir grandes áreas e ser extremamente mais econômico quando comparado

com um método direto, como por exemplo, a perfuração de poços. Segundo

(Robinson e Treitel, 1980) a exploração de hidrocarbonetos, óleo e gás, baseada

em sísmica pode ser dividida em três etapas: aquisição, processamento e

interpretação.

Aquisição Processamento Interpretação

Exploração (Óleo e Gás) baseada em Sísmica

Aquisição Processamento Interpretação

Exploração (Óleo e Gás) baseada em Sísmica Figura 1 – Etapas da exploração baseada na sísmica de reflexão.


Nas subseções seguintes faremos uma breve descrição de cada etapa.

Para uma visão mais aprofundada sugerimos consultar os trabalhos (Gerhardt,

1998), (Machado, 2000), (Robinson e Treitel, 1980) e (Yilmaz, 1987).

2.1.1. Aquisição

A aquisição é feita usando uma fonte para gerar ondas sísmicas que se

propagam abaixo da superfície da terra. Em aquisições terrestres é comum usar

explosões de dinamite como fonte, em aquisições marinhas são usados

normalmente dispositivos pneumáticos como canhões de ar. Quando a onda

sísmica alcança uma interface entre duas camadas de rocha com valores de

impedância acústica diferentes, parte da onda é refratada e continua viajando

para baixo; outra parte é refletida. A porção da energia refletida é proporcional a

diferença de impedância acústica entre os dois meios. A parte refletida da onda

retorna à superfície, é captada nos receptores (geofones em aquisições

terrestres ou hidrofones em marinhas) e gravada nos sismógrafos. O sismógrafo

armazena tanto o tempo de chegada da onda quanto a intensidade medida

neste momento. Após várias detonações variando a posição da fonte e dos

receptores, todo o dado armazenado é enviado para ser processado. A Figura 2

ilustra os processos de aquisição terrestre e marinha.

Figura 2 – Aquisição sísmica. Adaptada de (Gerhardt, 1998).


2.1.2. Processamento

Na etapa de processamento alguns erros inerentes ao fenômeno físico são

corrigidos no levantamento. Além disso, o dado é reorganizado para formar uma

grade tridimensional com uma amostra de amplitude sísmica em cada vértice da

grade (voxel). Duas das dimensões do dado são direções espaciais e estão

relacionadas com as posições das fontes e dos receptores. Uma das

transformações realizadas no dado durante o processamento faz com que as

posições da fonte e do receptor sejam a mesma. Também graças a esta

transformação podemos considerar a terceira dimensão do dado como sendo o

tempo e que a propagação da onda é feita apenas na direção vertical. Como

podemos considerar que a fonte e o receptor estão na mesma posição na

superfície, o tempo de cada amostra corresponde ao tempo que a onda leva

para viajar até uma interface mais o tempo da volta à superfície. Uma coluna de

amostras com mesmas coordenadas espaciais, variando apenas o tempo, é

chamada de traço sísmico.

A organização das amostras em um dado sísmico é mostrada na Figura 3.

Do lado esquerdo temos a função de amplitudes sísmicas do traço sísmico a

única dimensão é a temporal (1D). No centro temos uma seção vertical do dado,

esta seção vertical é formada por um conjunto de traços sísmicos e é chamada

de linha sísmica (2D); uma dimensão é espacial e a outra temporal. No caso do

dado sísmico 3D (volume sísmico), formado por várias linhas sísmicas, temos

duas direções temporais, que são chamadas de inline (direção das linhas

sísmicas) e crossline (direção perpendicular as linhas sísmicas), e uma direção

temporal.


AMPLITUDES

13

017

1925313743495561677379859197

103109115121

TEM

PO

AMPLITUDES

13

017

1925313743495561677379859197

103109115121

TEM

PO

AMPLITUDES

13

017

1925313743495561677379859197

103109115121

TEM

PO

Figura 3 - Traço sísmico (esquerda), Linha sísmica (centro) e volume sísmico (direita).

Um modelo matemático interessante que descreve bem o efeito do

processamento sísmico realizado sobre o dado é o modelo de convolução,

ilustrado na Figura 4. Neste modelo consideramos a função de amplitude

sísmica de cada traço do dado como sendo o resultado da convolução de um

impulso sísmico com uma função refletividade, a rigor a função refletividade é

uma distribuição de coeficientes de reflexão. Os coeficientes de reflexão são

proporcionais à diferença de impedância acústica entre camadas geológicas

adjacentes.

Cam

adas

Geo

lógi

cas

Coe

ficie

ntes

de R

efle

xão

∗

Am

plitu

des S

ísm

icas

Traço

ImpulsoSísmicoC

amad

asG

eoló

gica

s

Coe

ficie

ntes

de R

efle

xão

∗

Am

plitu

des S

ísm

icas

Traço

ImpulsoSísmico

Figura 4 – Modelo de convolução. Adaptado de (Gerhardt, 1998).


2.1.3. Interpretação

Na etapa de interpretação o intérprete, em geral um geólogo, analisa os

dados sísmicos e tenta criar um modelo que represente a geologia contida na

área do levantamento. A Figura 5 mostra um modelo geológico que poderia ser

resultante da interpretação de uma linha sísmica. A interpretação sísmica pode

ser classificada em dois tipos de acordo com o foco. A interpretação estrutural

basicamente tenta identificar as camadas geológicas, ou de forma equivalente as

interfaces entre as camadas, bem como as falhas geológicas que recortam as

camadas. Na interpretação estratigráfica o foco do trabalho está em entender a

maneira como as camadas foram se formando ao longo do tempo.

Figura 5 – Modelo geológico. Adaptada de (Robinson e Treitel, 1980).

(Sheriff, 1991) define um horizonte sísmico como sendo: a superfície

separando duas camadas diferentes de rocha, onde tal superfície (mesmo sem

ter sido identificada) está associada com uma reflexão que se estende por uma

grande área. Um horizonte sísmico se manifesta em um dado sísmico como uma

série de picos (ou vales) de amplitudes sísmica que aparecem de forma

consistente traço a traço. O mapeamento dos horizontes do dado é uma das

tarefas mais importantes da interpretação sísmica.

2.2. Visualização Volumétrica direta com texturas 3D

A visualização volumétrica é um processo de síntese de imagens onde os

dados de entrada são modelados matematicamente como um campo escalar

tridimensional. A versão discreta deste modelo implícito é um conjunto de

amostras tomadas em uma grade regular tridimensional. Cada amostra


armazena um número real correspondente ao valor do campo em um dado ponto

da grade. A grande maioria das aplicações quantiza as amostras para

armazená-las como inteiros de 8 bits. Como exemplos de modelos volumétricos

podemos citar os dados obtidos por tomografia computadorizada, ressonância

magnética, resultados de simulações físicas, equações matemáticas implícitas e

dados sísmicos.

Existem duas estratégias distintas na visualização volumétrica. Na

estratégia indireta visualizamos superfícies extraídas do modelo implícito. Na

abordagem direta utilizamos transparências para conseguir visualizar estruturas

tridimensionais complexas dentro do dado. Esta abordagem é indicada para

modelos onde a extração de superfícies é difícil.

Os algoritmos de visualização volumétrica direta produzem imagens

associando cor e opacidade a cada voxel. Esta associação é feita indiretamente,

isto é, uma tabela de cor e opacidade é criada para os valores armazenados nos

voxels. Por exemplo, em um algoritmo que usa uma representação de 8-bits

para os valores armazenados as tabelas de cor e opacidade possuem 256

entradas. A função que associa a cada valor do campo escalar uma cor e uma

opacidade é chamada de função de transferência.

Existem vários algoritmos de visualização volumétrica direta, entre eles

podemos citar ray casting (Levoy, 1988), shear-warp (Lacroute and Levoy, 1994)

e splatting (Westover, 1990) como alguns dos mais importantes.

O algoritmo de visualização volumétrica baseado em texturas é

equivalente ao algoritmo de ray casting e produz os mesmos resultados.

Enquanto no algoritmo de ray casting a imagem é construída lançando um raio

para cada pixel, na abordagem por textura vários raios são processados

simultaneamente, uma fatia 2D por vez. O dado é armazenado em textura 3D e

a fatia é utilizada para amostrar a textura. As fatias com textura aplicada são

combinadas para formar a imagem final. Como esta abordagem pode ser

acelerado por hardware, o algoritmo usando textura é mais rápido do que o ray

casting. Em sistemas onde o hardware gráfico não implementa texturas 3D uma

versão mais limitada pode ser implementada com texturas 2D. Porém, a

implementação com textura 3D é mais simples, e produz imagens com qualidade

melhor.


2.2.1. Amostragem da textura

A textura 3D em geral é amostrada usando planos perpendiculares a

direção de visualização. Porém, quando o observador está muito próximo ao

volume, ou mesmo dentro do volume, uma amostragem mais correta da textura

pode ser obtida utilizando conchas esféricas centradas na posição do

observador. A Figura 6 ilustra as geometria de amostragem formadas por fatias

planas e por conchas esféricas. A geometria de amostragem escolhida é

decomposta em polígonos recortados pelos limites do volume. A textura 3D

contendo o dado volumétrico é aplicada aos polígonos recortados. As fatias

texturizadas são compostas de trás para frente na direção do observador. A

composição é feita utilizando o operador over (Drebin et al, 1988).

Fatias planasFatias planas

(a)

ConchasConchas

(b)

Figura 6 – Geometrias de amostragem. Adaptada de (Blythe e McReynolds, 2000).

2.2.2. Composição

O operador over é a maneira mais comum de compor as amostras

transparentes na visualização volumétrica. Este operador também é utilizado em

outros algoritmos de visualização volumétrica direta como, por exemplo, o ray

casting. O operador over aproxima a integral de iluminação definida pelo

transporte da luz através do volume transparente. A opacidade (transparência)

de uma amostra está armazenada no valor do canal alfa do texel. Texels com

valores de alfa alto obstruem a visualização dos texels que estão atrás.

Segundo (Wittenbrink, 1998) para evitar artefatos de interpolação as cores

das amostras devem ser multiplicadas pela opacidade antes de serem


interpoladas. Portanto, o operador over pode ser implementado definindo a

função de blending do OpenGL da seguinte forma:

glBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);

2.2.3. Freqüência de amostragem

Segundo (Blythe e McReynolds, 2000) existem vários fatores que devem

ser levados em conta na hora de escolher o número de fatias (planos ou

conchas esféricas) usadas na amostragem da textura.

A performance que se deseja obter é um dos fatores mais importantes,

quanto maior é o número de fatias menor a performance. Em geral os sistemas

de visualização volumétrica baseados em textura operam em dois modos: um

modo ‘interativo’ onde são usadas poucas fatias, para garantir maior

performance; e um modo ‘detalhado’ onde são usadas mais fatias, para obter

uma imagem com maior qualidade.

Para evitar artefatos de re-amostragem é importante considerar os texels

como sendo cúbicos, ou seja, temos uma amostragem uniforme em todas

direções. Portanto, é interessante ter uma freqüência de amostragem entre as

fatias igual à freqüência de amostragem do dado. Isto pode ser conseguido

definindo o número de fatias igual a soma das dimensões do dado (em texels).

Aumentar exageradamente o número de fatias pode tornar difícil a

visualização da porção mais afastada do dado. O operador over é não-linear, por

isso é difícil compensar o aumento do número de fatias com uma atenuação da

opacidade das fatias.

Quando a visualização é feita com projeção em perspectiva a freqüência

de amostragem deve aumentar com a distância ao observador. Devido a

distorção da perspectiva a parte mais afastada do dado deve parecer mais

densa.

2.2.4. Gerenciando a memória de textura

O algoritmo de visualização volumétrica baseada em texturas 3D é capaz

de tratar volumes grandes que não cabem na memória de textura. Para isto o

volume é dividido em vários pequenos volumes, chamados de tijolos. A Figura 7


ilustra a subdivisão de um volume em tijolos. O tamanho dos tijolos é escolhido

para que cada um caiba na memória de textura.

Figura 7 – Volume de dados dividido em tijolos.

Durante a visualização cada tijolo é carregado na memória de textura,

amostrado por fatias e o resultado é composto, como descrito anteriormente.

Assim como as fatias, os tijolos são processados de trás para frente na direção

do observador.

Para evitar erros na re-amostragem do dado nas emendas dos tijolos, o

recorte dos polígonos e suas coordenadas de textura devem ser ajustadas para

não usarem os texels das faces dos tijolos. Além disso, cada tijolo é construído

de forma que os texels de suas faces sejam comuns, respectivamente, a cada

um dos seus 6 vizinhos. Com isso garantimos que os tijolos encaixem

corretamente na imagem resultante.

Podemos encontrar mais informações sobre paginação (paging) de

texturas em geral na seção 6.8 de (Blythe e McReynolds, 2000).

2.2.5. Iluminação

Os modelos de iluminação local aproximam pontualmente a interação entre

a luz e a superfície de um objeto. Em um modelo local, como o modelo de

Phong, a intensidade da luz refletida em um ponto é calculada em função da

orientação da superfície; da posição e direção da fonte de iluminação; e de

propriedades óticas do material do objeto. Nestes modelos a orientação local da

superfície é estimada pelo vetor normal em cada ponto. Portanto, para aplicar o

um modelo de iluminação local à visualização volumétrica é necessário estimar o


vetor normal em cada voxel. Para a maioria dos campos escalares o vetor

gradiente pode ser um substituto apropriado para o vetor normal da superfície,

uma vez que o gradiente é o próprio vetor normal da superfície de nível em cada

voxel. A Figura 8 mostra o modelo de Phong aplicado na visualização de uma

superfície de nível de um volume de dado médico. Na figura temos: (a) apenas a

componente ambiente; (b) componentes ambiente e difusa; (c) componentes

ambiente, difusa e especular.

(a) (b)

(c)

Figura 8 – Iluminação de uma superfície de nível de um volume. Adaptado de (Hadwiger

et al, 2002)

Existem diversas maneiras de calcular o gradiente em um dado

volumétrico. Em nossa implementação utilizamos o método diferenças centrais,

que é derivado dos primeiros termos da série de Taylor. Na visualização

volumétrica baseada em textura o campo gradiente é pré-computado e

armazenado como uma textura 3D. As três componentes do gradiente

normalizado são transformadas do intervalo [ ]1,1− para o intervalo [0,1] e

armazenadas como triplas RGB.

( , , )I x y z (1)

( 1, , ) ( 1, , )( , 1, ) ( , 1, )( , , 1) ( , , 1)

x

y

z

I I x y z I x y zI I I x y z I x y z

I I x y z I x y z

+ − − ∇ = = + − − + − −

(2)


x

y

z

nIn nI

n

∇ = = ∇

(3)

1

1 12

1

x

y

z

R nG nB n

+ = + +

(4)

Na Figura 9 temos uma visualização do campo vetorial do gradiente em

uma superfície de nível. As cores, definidas pelo mapeamento da eq. (4),

representam os vetores do campo gradiente do volume.

Figura 9 – Gradiente de uma superfície de nível. Adaptado de (Westermann e Ertl, 1998).

Em (Westermann e Ertl, 1998) é apresentada uma implementação para a

visualização de superfícies de nível, com cálculo de iluminação feito dentro da

placa gráfica. Ao desenhar os polígonos da geometria de amostragem o teste de

alfa do OpenGL é usado para filtrar os voxels da superfície de nível desejada. A

matriz de cor do OpenGL é utilizada para armazenar a direção da fonte de

iluminação. Após a projeção dos polígonos os pixels da imagem são copiados do

frame buffer para ele mesmo. A cor, contendo os vetores normais, dos pixels

copiados é multiplicada pela matriz de cor resultando no cálculo da equação de

iluminação.

A arquitetura da placa de vídeo GeForce3 da NVIDIA implementa em

hardware os combinadores de registro (register combiners) além da textura 3D

(Spitzer, 2001). Os combinadores de registro permitem que se façam operações

sobre texturas. Uma implementação do cálculo da iluminação feita no hardware,

baseada nesta arquitetura é apresentada em (Rezk-Salama, 2000). A

implementação do nosso sistema é baseada nesta abordagem. As placas de


vídeo da ATI, a partir da RADEON 8500 implementam um recurso de hardware

chamado de framgment shaders que são similares aos combinadores de registro

da NVIDIA.

Atualmente as placas gráficas programáveis permitem uma implementação

simples do cálculo de iluminação. Esta implementação pode ser feita usando a

linguagem Cg (Fernando e Kilgard, 2003).

2.2.6. Função de Transferência

A função de transferência tem um papel fundamental na visualização

volumétrica direta. É através da função de transferência que podemos

determinar as porções visíveis do volume e definir a cor e a transparência com

que estas irão aparecer na imagem. A Figura 10 mostra a visualização

volumétrica de um dado obtido por tomografia computadorizada de uma peça

mecânica. Na imagem à esquerda foi aplicada uma função de transferência que

associa uma cor ciano e opacidade alta (igual a um) aos valores de baixa

densidade, tornando-os totalmente opacos, e associa opacidade baixa (zero) aos

valores de densidade alta tornando-os invisíveis. A imagem do centro foi obtida

com tornando os voxels de densidade baixa associando a estes: opacidade zero,

e opacidade igual a um para os valores de densidade alta que também

receberam a cor vermelha. Na imagem da direita os valores de densidade baixa

tornaram-se transparentes recebendo um valor de opacidade intermediário

(entre zero e um).

Figura 10 – Voxels com densidade baixa (esquerda), voxels com densidade alta (centro),

voxels com densidade baixa e alta (direita).

A função de transferência pode ser vista como uma tabela que associa cor

e transparência (RGBA) a cada valor do volume. A função de transferência

podes ser implementada utilizando uma das tabelas (look-up tables) do OpenGL.

A extensão SGI_texture_color_table consulta a tabela durante a aplicação


da textura, após o valor do fragmento ter sido interpolado. Usando a extensão

EXT_color_table podemos usar uma tabela para converter os valores do

volume para RGBA durante o carregamento da textura, porém a textura será

uma textura RGBA. Se a extensão EXT_paletted_texture estiver disponível

podemos armazenar o volume como uma textura de índices. Esta é uma opção

interessante, pois podemos ajustar a imagem rapidamente a cada modificação

na função de transferência. Atualmente, com as placas de vídeo programáveis

podemos implementar a função de transferência usando um Fragment Program.

2.2.7. Algoritmo

O quadro abaixo apresenta, em linhas gerais, os passos básicos de um

algoritmo de visualização volumétrica direta baseada em textura 3D (Blythe e

McReynolds, 2000). Aqui consideramos a função de transferência implementada

usando a extensão EXT_paletted_texture e que a iluminação é feita por

meio de combinadores de registro. Consideramos também que existe uma única

fonte de luz que é direcional e branca (monocromática). Apenas as componentes

ambiente e difusa do modelo de Phong serão levadas em conta no cálculo da

iluminação.

1. Carregar o dado volumétrico em uma textura 3D.

2. Carregar o campo de vetores normais em uma textura 3D.

3. Definir o número de fatias.

4. Habilitar o operador de composição (over).

5. Carregar a direção da iluminação na cor primária.

6. Configurar os combinadores de registro (iluminação).

7. Carregar a função de transferência.

8. Encontrar a posição do observador e a direção de visualização.

9. Computar os polígonos da geometria de amostragem. Usar geração

automática de coordenadas de textura.

10. Usar a matriz de transformação de texturas para orientar corretamente a

aplicação da textura 3D nas fatias.

11. Visualizar os polígono de trás para frente na direção do observador.


A Figura 11 ilustra como os dados que descrevem a visualização

volumétrica são mapeados para a placa de vídeo. O dado volumétrico é

carregado em uma textura 3D de índices, ou seja, o valor de cada amostra deve

ser mapeado para um inteiro de 8 bits. O campo de vetores normais deve ser

mapeado em uma textura 3D de valores RGB, cada componente, após ter sido

transformada usando a eq. (4), deve ser convertida para um inteiro de 8 bits. A

função de transferência é mapeada em uma tabela com 256 entradas onde cada

entrada contém um quádrupla RGBA, esta tabela é definida como palheta da

textura de índices usando funções da extensão EXT_paletted_texture. A

direção da iluminação também deve ser transformada usando a eq. (4) e em

seguida armazenada em RGB como cor primária. Desta forma os combinadores

de registro podem acessar a direção de iluminação para usá-la no cálculo da

componente difusa.


Dado volumétrico I TEX0 – Textura 3D de Índices

Campo vetorial de normais n TEX1 – Textura 3D (RGB)

R

ABGR

ABGR

ABG

Função de Transferência 1D Palette

COR_PRIMÁRIA

Fonte de luz direcional

Memória Principal Placa de Vídeo

Figura 11 – Mapeamento de dados da memória principal para a placa de vídeo.

O calcula da iluminação é feita por pixel no momento em que cada

polígono é desenhado. A Figura 12 apresenta um esboço do processamento que

é feito pela placa de vídeo durante a visualização. A geração automática de

coordenadas de textura determina para cada vértice do polígono uma tripla

( , , )i i is t r de coordenadas. Durante a rasterização, para cada pixel do polígono é

calculada uma tripla de coordenadas de textura ( , , )s t r interpolando-se as


coordenadas de textura dos vértices. As coordenadas de textura interpoladas

definem um texel em cada unidade de textura.

No caso da figura abaixo temos duas unidades de textura. O texel da

unidade de textura TEX0 determinado pelas coordenadas de textura ( , , )s t r é

um índice que por sua vez determina uma entrada RGBA na tabela da função de

transferência. Essa entrada RGBA contendo a cor e a opacidade da amostra é

então mandada para o combinador final. Na unidade de textura TEX1, o texel

determinado pelas coordenadas de textura ( , , )s t r contem uma tripla RGB que

armazena o vetor normal da amostra. Dentro do combinador de registros o vetor

normal e a cor primária são convertidos de valores RGB novamente para

coordenas e é calculado o produto interno entre eles. O resultado do produto

interno é armazenado em cada componente de outra tripla RGB que também

segue para o combinador final.

No combinador final, a cor vinda da função de transferência (palheta) é

multiplicada pelo resultado do produto interno vindo do combinador de registros.

Esta multiplicação é feita coordenada a coordenada, por isso o combinador de

registros repete o valor do produto interno nas três componentes da sua saída

RGB. O combinador final também pode ser usado para adicionar a componente

de cor ambiente do modelo de Phong. O combinador final não altera o valor da

opacidade da amostra. Finalmente, a saída RGBA do combinador final é

combinada com o RGBA que está no Frame Buffer pelo operador de composição

(over).


Figura 12 – Cálculo de iluminação feito por pixel durante a visualização.

2.2.8. Funções de transferência multidimensionais

As funções de transferências multidimensionais são usadas em dados

volumétricos onde temos mais de um valor por amostra. Estes dados podem ser

campos vetoriais propriamente ditos, uma reunião de diferentes medidas feitas

nas amostras ou uma reunião de diferentes grandezas calculadas nas amostras.

A dimensão da função de transferência é igual ao número de valores por

amostra. O uso de funções de transferência multidimensionais pode aumentar a

capacidade do usuário selecionar a porção visível do dado porém, a

complexidade da tarefa de especificar uma função de transferência de dimensão

alta pode tornar sua aplicação inviável.

O modelo de visualização volumétrica descrito em (Levoy, 1988) já incluía

uma função de transferência 2D. Neste modelo cada amostra possui dois

valores: o próprio valor do campo escalar como primeira dimensão, e o módulo

do gradiente como segunda dimensão. Esta abordagem usa uma função de

transferência fixa que serve para realçar as bordas entre regiões homogêneas

do dado. Extensões para esta abordagem são apresentadas em (Kindlmann e

Durkin, 1998; Kniss, 2002), nestas extensões o usuário pode manipular a função

de transferência. Na Figura 13 e na Figura 14 temos a aplicação da função de


transferência 2D a um dado de tomografia computadorizada de um dente. Na

Figura 13a vemos a distribuição das amostras no plano definido pelos valores

campo escalar e pelo módulo do gradiente do campo. Cada pixel da imagem

corresponde a um par de valores do dado. Com uma ferramenta similar a um

sistema de edição de imagens o usuário pode definir cores para os pares de

valores, como na Figura 13b. Nesta imagem cada região colorida representa

contém os voxels da fronteira entre duas regiões homogêneas do volume.

(a) (b)

Figura 13 – Função de Transferência Bidimensional. Adaptada de (Pfister et al, 2001)

Na Figura 14 temos o resultado da visualização volumétrica com função de

transferência 2D aplicada ao dado do dente. A função de cor utilizada foi a

função da Figura 13b. No canto superior esquerdo de cada imagem temos a

opacidade 2D utilizada. Os tons de cinza representam a opacidade da amostra,

sendo branco para mais opaco e preto para menos opaco. Podemos ver que a

função de transferência bidimensional torna possível o realce das fronteiras entre

as regiões homogêneas do campo escalar.

I∇

IValor do Campo Escalar

Mód

ulo

do G

radi

ente

I∇


Mód

ulo

do G

radi

ente

I∇


Mód

ulo

do G

radi

ente

I∇


Mód

ulo

do G

radi

ente


(a)

(b)

(c)

(d) Figura 14 – Realce de fronteiras. Adaptada de (Pfister et al, 2001)

A função de transferência bidimensional pode ser implementada utilizando

o recurso de textura dependente disponível nas placas de vídeo da NVIDIA

(Spitzer, 2001) ou usando a programabilidade das placas de vídeo mais

modernas. A Figura 15 mostra como é feito o mapeamento dos dados da

memória principal para a placa de vídeo na implementação da função de

transferência 2D. Comparando a Figura 15 com a Figura 11 vemos o que muda

da implementação da função de transferência 1D para a 2D. O dado volumétrico

e o módulo do gradiente do dado são mapeados para uma única textura 3D de

valores RGB. O dado volumétrico para o canal G e o módulo do gradiente para o

canal B. O campo de normais e a direção da iluminação são mapeados da

mesma forma. A função de transferência 2D é mapeada para uma textura 2D de

valores RGBA.


Dado volumétrico I

G

B

G

B

TEX0 – Textura 3D (RGB)

Módulo do gradiente I∇

Campo vetorial de normais n TEX1 – Textura 3D (RGB)

Função de Transferência 2D TEX2 – Textura 2D (RGBA)

COR_PRIMÁRIA

Fonte de luz direcional

Memória Principal Placa de Vídeo

Figura 15 – Mapeamento de dados. Função de transferência 2D.


Na nossa implementação da função de transferência 2D utilizamos 3

unidades de textura. A tabela da função de transferência 1D é substituída por

uma textura dependente. O recurso de textura dependente, disponível nas

placas de vídeo mais recentes da NVIDIA, permite que os canais de cor de um

texel de uma textura sejam utilizados como índices para a textura dependente.

Na nossa implementação a unidade de textura TEX2, contendo a função de

transferência, é dependente da unidade de textura TEX0, contendo o dado e o

módulo do gradiente. Portanto, os canais G e B da unidade TEX0 são usados

como índice na unidade dependente TEX2, ou seja, a cor e a opacidade da

amostra dependem tanto do valor da amostra quanto do módulo do gradiente na

mesma amostra.

Comparando a Figura 16 com a Figura 12 vemos que na implementação

da função de transferência 2D o cálculo da iluminação é o mesmo. Porém, a

maneira de determinar a cor e a opacidade da amostra é diferente. A unidade de

textura TEX0 não armazena mais uma textura de índices, agora a textura em

TEX0 é 3D de valores RGB. As coordenadas de textura ( , , )s t r definem o valor

RGB do texel correspondente. O canais G e B são convertidos em coordenadas

de textura ( , )s t que são usadas para determinar um texel na unidade de textura

2D dependente TEX2. O valor RGBA correspondendo ao texel da unidade de

tetura TEX2 representa a cor e a opacidade da amostra que segue para o

combinador fianal.


(s0,t0,r0) (s1,t1,r1)

(s2,t2,r2)

(s,t,r)

Combinador Final

RGBA

RGB

(s,t,r)

Combinador

COR_PRIMÁRIA

RGB

PRO

DU

TO

INTE

RN

O

Frame Buffer

Mod

elo

de

Phon

g

RGBA

RGB

TEX0 TEX1

G

TEX2

B

(s0,t0,r0) (s1,t1,r1)

(s2,t2,r2)

(s,t,r)

Combinador Final

RGBA

RGBRGB

(s,t,r)

Combinador

COR_PRIMÁRIA

RGB

PRO

DU

TO

INTE

RN

O

Frame Buffer

Mod

elo

de

Phon

g

RGBA

RGB

TEX0 TEX1

G

TEX2

B

TEX2

B

Figura 16 – Função de transferência 2D.

Opacidade 2D 39

3 Opacidade 2D

Neste capítulo apresentaremos o estudo feito sobre a aplicação de funções

de transferência multidimensionais à visualização volumétrica de horizontes em

dados sísmicos 3D. Numa abordagem de visualização volumétrica direta a

função de transferência tem um papel muito importante, pois é ela quem define

como as estruturas tridimensionais contidas no dado aparecerão na imagem

final. A técnica de Opacidade 2D sugere o uso de uma função de transferência

bidimensional, onde a primeira dimensão é sempre o dado de amplitude sísmica.

A cor depende do valor da amplitude sísmica. A opacidade depende da

amplitude sísmica e de um outro atributo sísmico. Apresentaremos neste

capítulo a comparação entre 3 atributos diferentes. Para isto escolhemos uma

estrutura para este capítulo um pouco diferente da tradicional. Na seção seguinte

apresentaremos a motivação da pesquisa. Em seguida citaremos alguns

trabalhos relacionados mais importantes. Dentro da seção sobre soluções

propostas introduziremos os atributos na ordem em que foram testados na nossa

pesquisa, tentando motivar cada um a partir da deficiência do anterior. Após a

descrição de cada atributo apresentaremos os resultados obtidos pela

Opacidade 2D com aquele atributo. Ao final do capítulo faremos conclusões em

cima dos resultados obtidos. A implementação do sistema usado para gerar os

resultados deste capítulo está descrita no capítulo sobre conceitos básicos.

3.1. Motivação

Ao tentar utilizar técnicas de visualização volumétrica direta em dados

sísmicos 3D podemos experimentar alguns problemas. Por exemplo, não existe

um modelo “a priori” para a função de transferência, como acontece para dados

médicos. No caso de dados médicos, de acordo com a aquisição utilizada,

existem tabelas que classificam os valores dos voxels, que em geral são valores

de densidade, em tipos de tecidos. Com isso, para dados médicos, o trabalho de

especificar uma função de transferência para visualização, seja ele de

tomografia computadorizada ou de ressonância magnética, pode ser bastante

Opacidade 2D 40

simplificado. Outro problema surge da própria natureza oscilatória dos dados

sísmicos. O caráter oscilatório tem como conseqüência a superposição dos

intervalos de amplitude dos eventos sísmicos. A Figura 17 ilustra o problema da

superposição dos intervalos de amplitudes de dois eventos positivos (picos) em

um traço sísmico.

127

0

-127

TEMPO

A

B

127

0

-127

TEMPO

A

B

Figura 17 - Superposição dos intervalos de amplitudes dos eventos sísmicos

O primeiro evento possui valores intermediários de amplitude, estes

valores estão compreendidos no intervalo “A”. O segundo evento possui valores

de amplitude mais alta, porém, o seu intervalo de valores se estende por todo o

conjunto de valores positivos de amplitude, desta forma o intervalo de amplitudes

do segundo evento se superpõe ao intervalo “A”. A região “B” destaca a

interseção entre os dois intervalos de amplitude, ou seja, todos os voxels do

segundo evento que possuem valores de amplitude dentro do intervalo de

amplitudes do primeiro evento. Uma função de transferência definida da maneira

usual para visualizar apenas o evento de amplitude intermediária associaria

valores de opacidade igual a um dentro do intervalo “A” e opacidade igual a zero

fora de “A”. A aplicação de tal função de transferência em um algoritmo de

visualização volumétrica daria como resultado uma imagem onde apareceriam

os voxels do evento de amplitude intermediária, porém nesta imagem também

estariam os voxels da região “B” da superposição dos dois eventos.

Para exemplificar o problema da superposição dos intervalos de amplitude

dos eventos utilizamos um dado sintético que simula uma estrutura de domo

(anticlinal). Este dado é o mesmo que foi utilizado para testes em (Gerhardt,

2001). Este dado foi obtido tomando-se um traço sísmico de um dado real e

replicando este traço lateralmente para gerar o dado 3D. Ao replicar o traço

lateralmente foram feitos deslocamentos na direção vertical para curvar os

Opacidade 2D 41

eventos em forma de domo. O dado apresenta um evento com amplitudes altas

(abaixo) e destaca dois eventos de interesse, evento “A” e evento “B”, com

amplitudes intermediárias. Chamaremos este dado de modelo domo. A Figura 18

mostra a visualização do modelo domo com todos os voxels opacos.

Evento “A”

Evento “B”

Evento com amplitudes altas

-127 0 127

Evento “A”

Evento “B”


-127 0 127-127 0 127 Figura 18 - Dado sintético simulando uma estrutura de domo (anticlinal).

A Figura 19 mostra a imagem obtida na visualização volumétrica do

modelo domo. Esta imagem foi gerada na tentativa de visualizar apenas o

evento “B” (com amplitudes intermediárias). Foi utilizada a função de

transferência que aparece no canto superior esquerdo, a função aplica

opacidade igual a um para os valores de amplitude do evento “B” e opacidade

igual a zero para os valores de amplitude fora deste intervalo. Como o intervalo

de amplitudes do evento vizinho (com amplitudes mais altas) se sobrepõe ao

intervalo do evento “B”, podemos ver na parte inferior da imagem os voxels

referentes a essa região de interseção.

Opacidade 2D 42

Evento “B”


Amplitudes

Opa

cida

de

Imagem da V.V.

Evento “B”


Amplitudes

Opa

cida

de

Amplitudes

Opa

cida

de

Imagem da V.V. Figura 19 - Efeito da superposição de intervalos de amplitude na visualização

volumétrica de dados sísmicos.

Com isto podemos observar que a aplicação direta dos métodos de

visualização volumétrica tradicional utilizada para dados sísmicos 3D, não é

suficiente para que o usuário seja capaz de controlar as regiões de interesse que

deveriam aparecer na imagem final. Dentro do contexto de visualização de

dados sísmicos, podemos dizer que a aplicação direta da visualização

volumétrica tradicional é incapaz de isolar os eventos de interesse. Nas seções

subseqüentes apresentaremos uma nova abordagem para a visualização

volumétrica de dados sísmicos. Esta abordagem é baseada na utilização de

funções de transferência multidimensionais. Mais especificamente, no emprego

de uma função de opacidade 2D para tentar melhorar a capacidade de se

visualizar apenas os eventos de interesse.

3.2. Trabalhos Relacionados

O uso de funções de opacidade 2D na visualização de campos escalares

3D foi introduzido em (Levoy, 1988). Nesse trabalho são apresentados dois tipos

de função de opacidade, ambas bidimensionais, e ambas utilizando o módulo do

gradiente como a segunda dimensão. O primeiro tipo de função de opacidade é

destinado à visualização das interfaces entre materiais, ou seja, aplica-se à

visualização das fronteiras entre regiões homogêneas do dado. O outro tipo de

Opacidade 2D 43

função de opacidade é destinado à visualização de superfícies de nível em

dados com variação suave. Em (Kindlmann & Durkin, 1998; Kindlmann, 1999) é

apresentado um método semi-automático de geração de funções de opacidade

2D, usando a derivada na direção do gradiente como segunda dimensão, e

geração de funções de opacidade 3D acrescentando a segunda derivada,

também na direção do gradiente, como terceira dimensão. Em (Gerhardt et al.,

2001) fazemos uma adaptação da geração semi-automática de funções de

opacidade bidimensional à visualização de dados sísmicos 3D. Nesse trabalho

utilizamos a derivada do dado sísmico na direção do tempo (vertical) como

segunda dimensão. Nossa implementação atual é baseada em (Kniss, 2001).

Esse trabalho propõe elementos de interface com o usuário para manipulação de

funções de transferência multidimensionais, propõe também o uso de placas

gráficas programáveis e textura 3D na implementação destas funções. Nos

trabalhos citados até aqui é criada uma função de transferência preliminar

baseada em histogramas bidimensionais do campo vetorial, e em uma etapa

posterior a função de transferência pode ser editada para tentar se alcançar à

imagem desejada.

3.3. Soluções Propostas

A idéia central é estender o domínio da função de transferência para duas

dimensões, ou seja, os valores de cor e opacidade dependem de dois valores e

não mais de apenas um como usual. Os dados volumétricos, por sua vez, devem

armazenar dois valores por cada voxel. De certa forma passamos do problema

da visualização de um campo escalar para a visualização de um campo vetorial

(com vetores bidimensionais), no caso particular onde existe uma relação

conhecida entre as dimensões. Os resultados obtidos em (Gerhardt et al., 2001),

em dados sintéticos, confirmaram que a utilização da opacidade 2D aumenta a

capacidade de seleção de regiões de interesse. Porém esta técnica não

apresentou resultados significativos para dados reais. Por isso partimos para

uma abordagem mais adaptada à realidade da sísmica. Mantendo a amplitude

sísmica como primeira dimensão da função de opacidade, resolvemos

experimentar atributos sísmicos instantâneos como segunda dimensão.

Opacidade 2D 44

3.3.1. Amplitude Sísmica X Fase Instantânea

Em (Silva, 2003) descrevemos um método totalmente interativo para a

especificação de funções de opacidade 2D, onde o usuário define a função de

opacidade com o auxílio de um histograma também bidimensional do dado. Este

método utiliza como segunda dimensão o atributo sísmico de fase instantânea. O

atributo de fase instantânea pode ser calculado a partir da amplitude sísmica

( )X t de cada traço sísmico como.

( )X t (5)

Primeiro calculamos a função ( )Y t que é a transformada de Hilbert dos

valores de amplitude do traço sísmico ( )X t . De maneira simplificada podemos

dizer que a transformada de Hilbert de uma função é o resultado da convolução

desta função com o núcleo de Hilbert (mais detalhes no Apêndice A).

( ) ( ) ( )Y t h t X t= ∗ (6)

Em seguida definimos o traço analítico (ou traço complexo) como sendo a

função complexa ( )Z t , com parte real igual a ( )X t e parte imaginária igual a

( )Y t .

( ) ( ) ( )Z t X t iY t= + (7)

Tomando a forma polar do traço analítico obtemos a função real do atributo

de envelope (ou amplitude instantânea) ( )A t e a função real de fase instantânea

( )tθ , como na eq.(8).

( )( ) ( ) i tZ t A t e θ= (8)

2 2( ) ( ) ( )A t X t Y t= + (9)

1 ( )( ) tan( )

Y ttX t

θ − =

(10)

Opacidade 2D 45

Portanto o atributo de fase instantânea pode ser calculado como o

argumento complexo do traço analítico. A Figura 20 mostra os gráficos das

funções de amplitude e fase instantânea de um traço sísmico de um dado real.

-4

-3

-2

-1

0

1

2

3

4

1 7 13 19 25 31 37 43 49 55 61 67 73 79 85 91 97 103 109 115 121

Amplitude Fase instantânea-4

-3

-2

-1

0

1

2

3

4

1 7 13 19 25 31 37 43 49 55 61 67 73 79 85 91 97 103 109 115 121

-4

-3

-2

-1

0

1

2

3

4

1 7 13 19 25 31 37 43 49 55 61 67 73 79 85 91 97 103 109 115 121

Amplitude Fase instantâneaAmplitude Fase instantânea

Figura 20 - Amplitude sísmica e fase instantânea de um traço sísmico real

Os valores do atributo sísmico de fase instantânea estão compreendidos

no intervalo [−π ,π ]. Na Figura 20 observamos que onde o gráfico da função de

fase instantânea cruza o eixo do tempo, isto é, assume o valor zero temos um

máximo local com valores positivos na função de amplitude sísmica. Como um

máximo local com valores positivos de amplitude dentro de um traço sísmico

representa um evento sísmico concluímos que o atributo de fase instantânea é

um bom indicador de eventos sísmicos. As linhas tracejadas no gráfico da

função de fase instantânea representam os pontos onde a função salta de π

para -π .

3.3.1.1. Resultados

A Figura 21 mostra (do lado esquerdo) uma função de transferência

bidimensional definida em cima do histograma, também bidimensional, de

Amplitudes versus Fase Instantânea. O histograma bidimensional foi obtido a

partir do dado do domo sintético, usado como primeira dimensão (eixo horizontal

da figura), e calculando o volume do atributo sísmico de fase instantânea, usado

como segunda dimensão (eixo vertical da figura). Ambos os volumes foram

quantizados para 8-bits e a ocorrência dos pares foi computada no histograma.

Opacidade 2D 46

As freqüências (ocorrências) são apresentadas na figura em tons de cinza, em

preto para os pares com maior ocorrência, branco para nenhuma ocorrência e os

tons de cinza para ocorrências intermediárias. Em cima do histograma está

representada a função de transferência bidimensional definida pelo usuário. A

opacidade 2D está definida com o valor 1 dentro do retângulo e zero fora do

mesmo. As cores da função de transferência dependem apenas do valor da

primeira dimensão do dado, ou seja, a cor depende apenas do valor de

amplitude do voxel, como mostra o esquema de cores situado abaixo do

histograma bidimensional.

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Figura 21 - Função de transferência 2D (esquerda) e visualização volumétrica do dado

do Domo (direita).

Do lado direito da Figura 21 temos o resultado da visualização volumétrica

com a função de opacidade 2D mostrada do lado esquerdo. Os voxels que

aparecem na visualização são aqueles cujos valores de amplitude e fase

instantânea estão dentro da região indicada pelo retângulo. Este retângulo

representa a região onde a amplitude é positiva e a fase instantânea é próxima

de zero. Como foi observado anteriormente é na vizinhança dos picos que o

valor da fase instantânea vai para zero. Portanto temos todos os eventos com

valores de amplitude positiva do dado identificados na imagem. Podemos agora

utilizar uma outra função de opacidade 2D que torne visíveis apenas os voxels

do evento sísmico de interesse “B” da Figura 18. Tal função de opacidade 2D

pode ser construída utilizando um retângulo menor que limite os valores de

amplitude aos valores encontrados no evento “B”.

A Figura 22 mostra (do lado esquerdo) a função de opacidade 2D usada

para isolar o evento sísmico “B” de amplitudes intermediárias mostrado na figura

Opacidade 2D 47

5. Do lado direito da figura 8 vemos a visualização volumétrica do evento “B”

sem a interferência do evento de amplitudes altas situado logo abaixo. Com isso

mostramos que a utilização da opacidade 2D combinando amplitude sísmica e

fase instantânea foi capaz de solucionar o problema causado pela superposição

de intervalos de amplitude dos eventos sísmicos, para o dado do domo sintético.

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Figura 22 - Evento de interesse isolado graças a opacidade 2D. Função de transferência

(esquerda) e visualização volumétrica do evento (direita).

No caso de um dado sísmico real a dificuldade de isolar eventos é ainda

maior. O dado do domo sintético utilizado anteriormente possui uma

particularidade que raramente pode ser encontrada em um dado sísmico real,

este dado não possui variação lateral das amplitudes das amostras que

compõem o centro do evento, ou seja, os valores de amplitude praticamente são

constantes nas amostras que estão próximas dos máximos (ou mínimos) locais.

A Figura 23 mostra um volume de dados real onde naturalmente há variação

lateral de amplitudes. Do lado esquerdo é mostrada a visualização volumétrica

do dado utilizando uma função de opacidade igual a um para todos os valores e

uma tabela de cores definida para realçar a variação de valores positivos. Do

lado direito estão destacados dois eventos vizinhos, em ambos os eventos a

variação lateral de amplitudes positivas é visível. No evento situado acima

podemos observar uma predominância de amplitudes positivas altas (amostras

em vermelho) enquanto no evento situado mais abaixo há uma predominância

de amplitudes positivas intermediárias (amostras em verde). Novamente aqui,

como na experiência realizada com o dado sintético, queremos visualizar

isoladamente o evento de amplitudes intermediárias.

Opacidade 2D 48

Figura 23 - Volume de dados real (esquerda). Eventos sísmicos reais destacados

(direita).

Na Figura 24 comparamos o resultado da visualização volumétrica

utilizando opacidade 1D aplicada a um dado real com o resultado da opacidade

2D aplicada ao mesmo dado real. Na porção superior da figura estão a função

de transferência 1D usada na visualização (centro) e o resultado da visualização

(direita). A função de opacidade 1D foi construída atribuindo opacidade igual a

um para os valores de amplitude intermediária que correspondem aos valores de

amplitude do evento de interesse, e zero para os demais valores. Assim como

para o dado sintético, utilizando apenas a opacidade 1D baseada nos valores de

amplitudes não é possível remover totalmente o evento de amplitudes mais

altas. Ainda estão presentes na imagem obtida com a opacidade 1D os voxels

de amplitudes intermediárias que estão na vizinhança do evento de amplitudes

predominantemente altas. Estes voxels obstruem a visualização do evento de

interesse. Na parte inferior da Figura 24 vemos a função de transferência 2D

utilizada (centro) e a imagem obtida na visualização volumétrica (direita). A

função de opacidade 2D foi definida com valores iguais para os valores de

amplitude e fase dentro do retângulo destacado e zero fora dele. Este retângulo

delimita a região do plano: Amplitude X Fase instantânea, onde, a fase

instantânea é próxima de zero e os valores de amplitude estão na faixa de

amplitudes intermediárias, que compreende os valores de amplitude do evento

de interesse. A imagem resultante da aplicação da opacidade 2D mostra que a

maior parte do evento de amplitudes altas foi removida. A parte que foi removida

corresponde exatamente aos voxels com valores mais altos de amplitude

sísmica. Como ambos os eventos possuem uma certa variação lateral de valores

Opacidade 2D 49

de amplitude, boa parte do evento de ainda fica visível e uma porção do evento

de interesse situada no centro foi removida como pode ser verificado na imagem

resultante. Com isso concluímos que a utilização da opacidade 2D em dados

sísmicos reais também se mostrou mais eficiente que a abordagem tradicional

com opacidade 1D.

Amplitude

Opa

cida

de

Opacidade1D

Amplitude

Opa

cida

de

Amplitude

Opa

cida

de

Opacidade1D

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

iasOpacidade

2D

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

iasOpacidade

2D

Figura 24 - Função de opacidade 1D aplicada a um dado sísmico real (acima). Função

de opacidade 2D aplicada a um dado sísmico real (abaixo).

3.3.2. Amplitude Sísmica X Fase Ajustada

Apesar de resolver, de certa forma, o problema da superposição dos

intervalos de amplitude dos eventos sísmicos, a opacidade 2D definida no

domínio Amplitude Sísmica X Fase Instantânea, como foi descrito, possui uma

deficiência na visualização de eventos de amplitudes negativas, ou seja, eventos

definidos por mínimos locais na curva de amplitudes sísmicas. Esta deficiência é

decorrente do fato da função do atributo de fase instantânea ser descontínua.

Como foi visto na Figura 20, a função de fase instantânea possui pontos de

descontinuidade que coincidem com os mínimos locais da função de amplitude

sísmica. Deste modo os eventos de amplitude sísmica negativas aparecem no

histograma 2D (de Amplitude X Fase) distribuídos em torno das retas: Fase = π

Opacidade 2D 50

e Fase = -π . Na verdade, Fase = π e Fase = -π , são uma única reta já que π

e –π representam o mesmo valore de fase instantânea. Portanto a topologia

correta para representar o histograma 2D de Amplitude X Fase é a topologia de

um cilindro obtido identificando-se as retas: Fase = π e Fase = -π .

O salto de π para –π dado pela fase instantânea acarreta, pelo menos,

dois problemas na visualização de dados sísmicos 3D utilizando opacidade 2D.

O primeiro problema, não tão grave, está na especificação da função de

opacidade feita pelo usuário no caso do evento de interesse ser um evento de

amplitudes negativas. Como os eventos de amplitude sísmica negativas

aparecem divididos em torno das duas retas o usuário é obrigado a definir uma

região desconexa para selecionar o evento, isto é, o usuário é obrigado a definir

duas áreas (retângulos) para selecionar um único evento, o que torna a tarefa

menos intuitiva. O segundo problema reside na interpolação do dado

volumétrico, o que afeta diretamente o cerne dos algoritmos de visualização de

dados volumétricos, inclusive a utilização da opacidade 2D.

A Figura 25 mostra o problema da descontinuidade da fase instantânea na

visualização volumétrica. Do lado esquerdo está a função de transferência 2D

usada na visualização. Esta função foi definida com valores de opacidade um

dentro do retângulo azul e zero fora dele. A escala de cores está associada com

a primeira dimensão que é a amplitude sísmica. A segunda dimensão utilizada é

a fase instantânea. Fica evidente a existência de artefatos na imagem resultante,

pois a superfície de nível de fase instantânea igual a zero não contem voxels

com valores negativos, ou seja, todos os voxels que aparecem em tons de cinza

na imagem são artefatos. Em torno dos eventos negativos existem voxels

adjacentes um com valor de fase igual a π e outro com valor igual a –π . Ao

utilizar a opacidade 2D para selecionar os voxels com fase instantânea próxima

de zero deveríamos estar visualizando apenas voxels com valor de amplitude

sísmica positivos, porém na visualização volumétrica direta a interpolação do

dado introduz artefatos (voxels com amplitude negativa) na vizinhança dos

mínimos locais onde a fase salta de –π para π .

Opacidade 2D 51

π/2 π

−π −π/2

π/2

−π/2

Fase (θ)

Fase Ajustada (θ’)

π/2 π

−π −π/2

π/2

−π/2

Fase (θ)

Fase Ajustada (θ’)

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Amplitude

Fase

Inst

antâ

nea

π

−π

Freq

üênc

ias

Figura 25 - Artefatos na visualização volumétrica da Fase igual a zero.

Para contornar o problema da descontinuidade da função de fase

instantânea sugerimos um ajuste que deve ser feito durante o cálculo.

Acrescentando este ajuste, na verdade, estamos calculando um novo atributo

sísmico o qual chamamos de atributo de fase ajustada. O atributo de fase

ajustada ou função de fase ajustada é a composição da função de fase com a

função de ajuste de fase. A função de ajuste de fase aplica uma reflexão em

torno das retas de fase = π /2 e fase = -π /2 como mostra a Figura 26.

Figura 26 - Função de ajuste de fase. Expressão (esquerda) e gráfico (direita).

Pelo próprio gráfico da função de ajuste de fase podemos observar que os

valores de fase instantânea no intervalo [-π /2, π /2] ficam inalterados. Este

intervalo diz respeito, justamente, aos valores de fase instantânea dos eventos

de amplitude sísmica positiva. Os valores de fase instantânea maiores que π /2

e menores que –π /2 é que são realmente ajustados, pois estes estão

relacionados aos eventos de amplitude sísmica negativa.

−<−−

>−

=

contráriocaso

se

se

,2

,2

,

'

θ

πθθπ

πθθπ

θ

Opacidade 2D 52

A Figura 27 mostra, do lado esquerdo, uma função de opacidade 2D que

atribui opacidade igual a um aos voxels com valor de fase instantânea próximo

de zero, independente do valor de amplitude. Do lado direito da figura vemos a

mesma função e transferência porem, definida em cima do histograma 2D de

Amplitude X Fase ajustada. Podemos observar na figura que os eventos de

amplitude sísmica negativas que apareciam, no histograma 2D de Amplitude X

Fase instantânea, em torno das retas fase = π e fase = -π agora aparecem, no

histograma 2D de Amplitude X Fase ajustada, em torno da reta fase igual a zero.

Portando, tanto os eventos positivos quanto os negativos ocorrem na vizinhança

de fase ajustada igual a zero. Com o uso da fase ajustada como segunda

dimensão, na visualização volumétrica com opacidade 2D, conseguimos

contornar o problema do aparecimento de artefatos.

Amplitude

Fase

Inst

antâ

nea

π

−πFr

eqüê

ncia

sAmplitude

Fase

Aju

stad

a

π/2

−π/2

Figura 27 - Histogramas 2D de Amplitude X Fase Instantânea (esquerda) e Amplitude X

Fase Ajustada (direita).

3.3.2.1. Resultados

Na Figura 28 comparamos os resultados obtidos usando fase instantânea

e fase ajustada na opacidade 2D. Na parte superior da figura temos, do lado

esquerdo, a função de transferência 2D definida em cima do histograma 2D de

amplitude sísmica e fase instantânea. Do lado direito temos a imagem da

visualização volumétrica com opacidade 2D usando esta função de

transferência. Nesta imagem podemos observar os artefatos (voxels em azul)

criados pela descontinuidade da fase instantânea.

Opacidade 2D 53

Figura 28 – Comparação entre Fase Instantânea e Fase Ajustada.

Na parte inferior da Figura 28, do lado esquerdo, vemos a função de

transferência definida em cima do histograma 2D de amplitude sísmica X fase

ajustada. Esta função define opacidade um apenas para os voxels com valor de

fase ajustada próximo de zero independente do valor de amplitude sísmica. Do

lado direito temos a imagem da visualização volumétrica com opacidade 2D do

dado. Podemos observar que esta imagem não apresenta artefatos oriundos da

interpolação já que a fase ajustada é contínua.

O fato da função de fase ajustada colapsar os eventos positivos e

negativos em torno do zero, ou seja, tanto as amostras dos eventos positivos

quanto as amostras dos eventos negativos possuem fase ajustada próxima de

zero, não cria um novo problema para a visualização volumétrica com opacidade

2D, pois ainda podemos separá-los pelo valor da amplitude sísmica como mostra

a Figura 29.

Amplitude

Fase

Aju

stad

aFa

se In

stan

tâne

a

π

−πAmplitude

Amplitude

Fase

Aju

stad

a

Amplitude

Fase

Aju

stad

aFa

se In

stan

tâne

a

π

−πAmplitude

Fase

Inst

antâ

nea

π

−πAmplitude

Opacidade 2D 54

Amplitude

Fase

Aju

stad

a

Amplitude

Fase

Aju

stad

a

Amplitude

Fase

Aju

stad

a

Amplitude

Fase

Aju

stad

a

Amplitude

Fase

Aju

stad

a

Amplitude

Fase

Aju

stad

a

Figura 29 – Visualização volumétrica dos eventos positivos e negativos separadamente

usando a fase ajustada.

Na parte superior da Figura 29 vemos, do lado esquerdo, a função de

transferência que seleciona os voxels com fase ajustada próxima de zero e

amplitude sísmica positiva, que representam os eventos de amplitude positiva.

Do lado direito temos a visualização volumétrica com opacidade 2D usando

amplitude sísmica e fase ajustada. Na parte inferior da figura temos, do lado

esquerdo, a função de transferência associa opacidade igual a um aos voxels

com fase ajustada próxima de zero e amplitude sísmica negativa, este conjunto

equivale aos eventos de amplitude negativa. Do lado direito temos a visualização

volumétrica com opacidade 2D usando esta função de transferência. Na imagem

estão apenas os eventos de amplitude negativa do modelo.

Com isso mostramos que a fase ajustada é um atributo mais apropriado

para ser usado em conjunto com a amplitude sísmica na opacidade 2D, uma vez

que não temos problema na visualização de eventos de amplitude negativa.

Opacidade 2D 55

3.3.3. Amplitude Sísmica X Fase Desenrolada

Para dados sísmicos reais a estratégia de visualização volumétrica com

opacidade 2D, tanto utilizando fase instantânea quanto fase ajustada, mostrou-

se eficiente para separar os horizontes sísmicos do restante do dado. Porém,

quando necessitamos de isolar um evento específico dos demais isso nem

sempre é possível. Dentro dessas duas abordagens só conseguimos

individualizar um evento no caso em que seus voxels possuem valores de

amplitude dentro de um intervalo disjunto do intervalo de amplitudes dos voxels

dos demais eventos.

Na tentativa de melhorar a capacidade de individualizar um único horizonte

sísmico na visualização volumétrica fizemos um estudo em (Gerhardt et al, 2004)

da aplicação do atributo de fase desenrolada1 como segunda dimensão na

opacidade 2D. A abordagem usando a fase desenrolada é uma extensão natural

da abordagem usando fase instantânea sugerida em (Silva et al, 2003). Em

(Stark, 2003) o desenrolamento da fase é apresentado como um método para

obter um volume de tempos para as amostras, neste volume de tempos

relativos, segundo o autor, os horizontes se apresentam como superfícies de

valor constante. A fase instantânea obtida a partir do traço complexo é medida

módulo 2π , como foi definido na eq. (10). A fase instantânea medida módulo

2π também é conhecida como VP (valor principal). A fase desenrolada φ pode

assumir qualquer valor real e pode ser obtida adicionando-se múltiplos de 2π à

fase instantânea. Na eq. (11) temos a formulação 1D do problema de desenrolar

a fase instantânea.

( ) ( ) 2 ( ) ( )t t n t n tφ θ π= + ∴ ∈ (11)

O valor ( )n t é chamado de índice de rotação ou número de voltas. Os

algoritmos de desenrolar fase trabalham naturalmente na versão discreta do

problema procurando por uma função de índices de rotação que tornem a fase

desenrolada contínua ou, pelo menos minimizem as suas descontinuidades. O

algoritmo de Itoh (Ghiglia & Pritt, 1998) é um método de desenrolar fase em 1D

1 Usamos a expressão fase desenrolada como tradução para a expressão em

inglês Unwrapped Phase.

Opacidade 2D 56

que o faz integrando a derivada da fase instantânea. AFigura 30 mostra o

resultado do algoritmo de Itoh aplicado a um traço sísmico real.

Figura 30 – Fase desenrolada em 1D.

Opacidade 2D 57

Os principais problemas que podem comprometer os algoritmos de

desenrolamento de fase são: a presença de ruído no sinal de fase instantânea e

o aliasing causado pela amostragem pobre do sinal. Podemos encontrar na

literatura várias versões de métodos de desenrolamento de fase em 2D (Ghiglia

& Pritt, 1998; Spagnolini, 1993). A partir da dimensão dois o problema de

desenrolar fase torna-se um problema variacional e sua solução depende da

integração do gradiente do campo escalar de fase instantânea. Sob certas

condições essa integral pode ser calculada de forma independente do caminho

(Ghiglia & Pritt, 1998). Um dos métodos de desenrolamento de fase em 2D, que

pode ser facilmente estendido para um espaço N-dimensional, é o método

baseado em mapa de qualidade. Este algoritmo define uma função que atribui

um valor de qualidade (número real) a cada amostra do campo escalar de fase

instantânea, mapa de qualidade. Em seguida encontramos a amostra com maior

valor de qualidade. A partir daí é utilizada uma estratégia do tipo de crescimento

de região. A cada passo encontramos a amostra da vizinhança da região atual,

com maior valor de qualidade. Esta amostra é acrescentada à região. Em

(Ghiglia & Pritt, 1998) são comparados vários mapas de qualidade 2D.

Com o atributo de fase desenrolada calculado para um dado sísmico

podemos utilizar a estratégia de opacidade 2D e individualizar um único evento

na visualização volumétrica. Isso se deve ao fato de que, assim como na fase

instantânea, um horizonte é uma superfície de nível do atributo de fase

desenrolada. Além disso, no atributo de fase instantânea todos os horizontes

formam uma única superfície de nível, a saber ( , , ) 0t x yθ = , como foi visto,

porém no atributo de fase desenrolada cada horizonte é uma superfície de nível

diferente. Na seção seguinte mostraremos os resultados obtidos utilizando a

estratégia de opacidade 2D com amplitude sísmica como primeira dimensão e

fase desenrolada como segunda dimensão.

3.3.3.1. Resultados

A Figura 31 ilustra a fase desenrolada bidimensional. Da esquerda para a

direita temos: uma seção transversal (inline) do dado do domo sintético, com os

valores de amplitude sísmica desta seção; os valores do atributo de fase

instantânea da seção; os valores do atributo de fase desenrolada da seção,

esses valores estão expressos em radianos; por fim, os valores dos índices de

rotação da seção. Esta imagem foi obtida aplicando o algoritmo de

Opacidade 2D 58

desenrolamento da fase em 3D ao modelo do domo. Na imagem da fase

instantânea podemos observar que as descontinuidades (saltos de π para π− )

coincidem com os eventos de amplitudes negativas do dado. Também podemos

observar que a fase desenrolada é contínua na seção. Assim como em 1D a

função de índice de rotação divide o domínio em regiões onde o índice de

rotação é constante, as fronteiras entre estas regiões também coincidem com as

descontinuidades (saltos) da função de fase instantânea. Alguns algoritmos de

desenrolamento de fase operam identificando estas fronteiras, porém em dados

reais essa identificação quase sempre é impraticável.

-0.958 0 0.958-0.958 0 0.958 -π +π0-π +π0 1.58− 32.251.58− 32.25 0 1 2 3 4 50 1 2 3 4 5 Figura 31 – Fase desenrolada em 2D. Dado sintético.

Na Figura 32 vemos o desenrolamento da fase aplicado a um dado real.

Esta figura também foi obtida da aplicação do algoritmo de desenrolamento de

fase em 3D ao dado. Na figura vemos uma seção transversal (crossline) do dado

real 3D. Na parte superior esquerda temos a imagem dos valores de amplitude

da seção. Nesta imagem podemos observar que a estrutura geológica da região

é relativamente complexa, contendo falhas geológicas e regiões com bastante

ruído. Na parte superior direita temos a imagem dos valores de fase instantânea

da seção. Nesta imagem podemos observar que as descontinuidades da fase

instantânea não formam um padrão bem comportado a ponto de podermos

utilizá-las para delimitar regiões onde o índice de rotação é constante. Na parte

inferior esquerda temos a imagem dos valores da fase desenrolada da seção.

Nesta imagem podemos observar que as curvas de nível da imagem de

acompanham, até certo ponto, o comportamento da geologia. Porém

examinando com mais cuidado, podemos constatar que estas curvas de nível

não mantêm uma relação um para um com os horizontes sísmicos como é

desejável.

Opacidade 2D 59

0-1769.34 1769.340-1769.34 1769.34 0−π +π0−π +π

-16.02 324.99-16.02 324.99 3− 523− 52 Figura 32 – Fase desenrolada em 2D. Dado real.

Na parte inferior direita da Figura 32 temos a imagem dos valores dos

índices de rotação da seção (crossline) do dado real 3D. Pela imagem e pelo

intervalo de valores desse dado podemos observar que existem por volta de 55

índices de rotação diferentes. A eq. (11) nos mostra que existe uma variação

grande dos valores de fase desenrolada entre as amostras com um mesmo

índice de rotação. Portanto, a quantização do atributo de fase desenrolada para

8-bits pode acarretar um problema de aliasing e prejudicar a visualização.

Na Figura 33 comparamos a opacidade 2D usando amplitude sísmica

como primeira dimensão e fase instantânea como segunda, com a opacidade 2D

usando amplitude sísmica como primeira dimensão e a fase desenrolada como

segunda. Na imagem superior esquerda temos a função de transferência

Opacidade 2D 60

definida sobre o histograma 2D de amplitude sísmica e fase instantânea. Na

parte superior direita vemos a visualização volumétrica do dado sintético usando

esta função de transferência.

+π-π

Am

plitu

des

Fase Instantânea

Freq

üênc

ias

-0.96

0.96

+π-π

Am

plitu

des

Fase Instantânea

Freq

üênc

ias

-0.96

0.96

Am

plitu

des

Fase Desenrolada

Freq

üênc

ias

-0.96

0.96

-1.58 32.25

Am

plitu

des

Fase Desenrolada

Freq

üênc

ias

-0.96

0.96

-1.58 32.25

Figura 33 – Comparação da opacidade 2D com fase instantânea e fase desenrolada.

Dado sintético.

Na parte inferior da Figura 33 temos, do lado esquerdo, a função de

transferência definida em cima do histograma 2D de amplitude sísmica e fase

desenrolada do dado sintético do domo, e do lado direito temos a visualização

volumétrica do dado sintético do domo utilizando esta função de transferência.

Enquanto no histograma 2D de amplitude sísmica e fase instantânea as

amostras dos horizontes aparecem aglomerados, i.e., as amostras referentes a

todos os horizontes de amplitudes positivas estão em torno da reta 0θ = , no

histograma 2D de amplitude sísmica e fase desenrolada as amostras estão em

uma curva que se assemelha ao próprio traço sísmico.

Opacidade 2D 61

Na Figura 34 comparamos os resultados obtidos aplicando a opacidade 2D

a um dado real, usando a fase instantânea e usando a fase desenrolada como

segunda dimensão. Na parte superior temos: do lado esquerdo, uma função de

transferência definida sobre o histograma 2D de amplitude sísmica e fase

instantânea; do lado direito, a visualização volumétrica usando esta função de

transferência. Na imagem do histograma 2D de amplitude sísmica e fase

instantânea vemos que as amostras dos horizontes positivos aglomeram-se em

torno do segmento de reta fase = 0 e amplitude sísmica > 0. Por isso não

conseguimos distinguí-los entre si. Na parte inferior da Figura 34 temos: do

esquerdo, uma função de transferência definida sobre o histograma 2D de

amplitude sísmica e fase desenrolada; do lado direito, a visualização volumétrica

usando esta função de transferência. Com a opacidade 2D utilizando a fase

desenrolada como segunda dimensão conseguimos individualizar um horizonte.

-π

Am

plitu

des

Fase Instantânea

Freq

üênc

ias

+π-1769

1769

-π

Am

plitu

des

Fase Instantânea

Freq

üênc

ias

+π-1769

1769

Am

plitu

des

Fase Desenrolada

Freq

üênc

ias

-16.02 325.01-1769

1769

Am

plitu

des

Fase Desenrolada

Freq

üênc

ias

-16.02 325.01-1769

1769

Figura 34 – Comparação da opacidade 2D com fase instantânea e fase desenrolada.

Dado real.

Opacidade 2D 62

3.4. Conclusões Parciais

Neste capítulo concluímos que o atributo sísmico de fase instantânea é

uma boa opção para ser usada como segunda dimensão na abordagem de

opacidade 2D, devido a sua já conhecida característica de ser um bom indicador

de continuidade para horizontes sísmicos. Com a opacidade 2D definida sobre a

amplitude sísmica e a fase instantânea conseguimos isolar os horizontes

sísmicos positivos do restante do dado.

Porém, devido ao fato da fase instantânea ser um atributo sísmico

descontínuo, i.e., possuir saltos de π para π− , e esses saltos estarem

localizados exatamente nos horizontes sísmicos negativos, concluímos que a

fase instantânea não serve para ser usada como segunda dimensão na

opacidade 2D, no caso da visualização de horizontes sísmicos negativos.

O atributo de fase ajustada sugerido nas seções anteriores é uma solução

para contornar o problema da descontinuidade da fase instantânea na

visualização de horizontes negativos. Com a opacidade 2D definida sobre a

amplitude sísmica e a fase ajustada conseguimos isolar os horizontes positivos e

também os horizontes negativos do restante do dado. Além disso, conseguimos

separa os horizontes positivos dos negativos.

Quando dois horizontes possuem intervalos de amplitude que se

intersectam nem usando a fase instantânea nem a fase ajustada na opacidade

2D conseguimos distinguir estes eventos. Para tal feito o atributo de fase

desenrolada se apresenta como uma opção mais apropriada. Mais do que isso,

para podermos visualizar isoladamente cada horizonte de um dado sísmico 3D

basta ter o atributo de fase desenrolada corretamente calculado.

Iluminação de Dados Sísmicos 63

4 Iluminação de Dados Sísmicos

Neste capítulo enunciaremos alguns dos problemas encontrados na

visualização volumétrica direta de dados sísmicos, mais especificamente na

etapa de iluminação, e apresentaremos uma solução para estimar a orientação

local dos horizontes em cada voxel do dado. Na seção seguinte apresentaremos

a motivação para a pesquisa feita. Na seção 4.2 destacaremos alguns trabalhos

relacionados. Na seção 4.3 apresentaremos a nossa solução. Na seção 4.4

apresentaremos os resultados obtidos usando nossa abordagem para dados

sintéticos e reais. Ao final do capítulo faremos considerações sobre os

resultados obtidos. A implementação do sistema usado para gerar os resultados

deste capítulo está descrita no capítulo sobre conceitos básicos.

4.1. Motivação

O efeito da iluminação em um ambiente traz para o observador

informações muito importantes sobre a geometria e a topologia dos objetos

contidos neste ambiente. Por isso, na grande maioria dos métodos de síntese de

imagens conhecidos na computação gráfica, a etapa de iluminação tem um

papel importante no que diz respeito ao realismo e na sensação de

tridimensionalidade das imagens resultantes. No caso da visualização de dados

sísmicos 3D não é diferente, a iluminação pode acrescentar bastante informação

visual revelando detalhes estruturais, geometria e topologia aparente dos

eventos sísmicos. Porém, como já foi observado em trabalhos anteriores

(Gerhardt, 1998), a aplicação direta das técnicas de visualização volumétrica

direta, comumente utilizadas em dados de origem médicas, em dados sísmicos

resulta em imagens onde os eventos são iluminados erradamente.

Para ilustrar o problema da iluminação de dados sísmicos 3D

consideraremos inicialmente um modelo sintético simples de camadas planas

horizontais. Esse modelo simples chamado de modelo plano foi obtido a partir

de um único traço de um dado real. Para gerar o modelo plano tomamos este


traço sísmico real, e replicamos este traço lateralmente sem deslocá-lo na

vertical.

Seja )(0 tX a função que representa o traço gerador do modelo plano, i.e.,

a função que representa a variação da amplitude sísmica ao longo do traço. A

Figura 35 mostra o gráfico da função )(0 tX , nela estão destacados os eventos

(pontos de máximo e de mínimo) do traço.

-0,8-0,6-0,4-0,2

00,20,40,60,8

11,2

1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 50

- Traço gerador do modelo plano

Am

plitu

de S

ísm

ica

T

)(0 tX

-0,8-0,6-0,4-0,2

00,20,40,60,8

11,2

-0,8-0,6-0,4-0,2

00,20,40,60,8

11,2

1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 501 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 50

- Traço gerador do modelo plano

Am

plitu

de S

ísm

ica

T

)(0 tX

Figura 35 - Traço gerador do modelo domo com os eventos destacados.

A partir da função )(0 tX podemos definir o campo escalar de amplitudes

sísmicas do modelo plano como na eq. (12) abaixo.

0( , , ) ( )X t x y X t= (12)

Este dado sintético representa um dos modelos geológicos mais simples

que é um conjunto de camadas horizontais, planas e sem variação lateral

amplitudes. A Figura 36 mostra a visualização volumétrica do dado do modelo

plano.


-0.95 0 0.95-0.95 0 0.95 Figura 36 – Visualização volumétrica do modelo plano.

Os métodos de visualização volumétrica direta em geral utilizam modelos

locais para calcular a iluminação em cada voxel do volume. Os modelos locais

de iluminação, por sua vez, calculam a intensidade da iluminação em um ponto

de uma superfície como função do vetor normal da superfície naquele ponto. O

modelo de iluminação local mais utilizado é o modelo de Phong (Phong, 1975).

No caso da visualização de um dado volumétrico, o gradiente do campo escalar

é utilizado para aproximar o vetor normal de cada voxel.

Tomando o gradiente do campo escalar de amplitudes sísmicas do modelo

plano obtemos o campo vetorial descrito na eq. (13):

0 ( )( , , ) , , ,0,0X tX X XX t x yt x y t

∂∂ ∂ ∂ ∇ = = ∂ ∂ ∂ ∂ (13)

Da eq.(13) podemos observar que nos máximos e mínimos locais da

função geradora )(0 tX a derivada 0 ( )X t t∂ ∂ é zero e conseqüentemente o

gradiente é nulo. Portando, exatamente nos eventos sísmicos o gradiente é nulo

e o vetor normal não está bem definido. Por isso, as imagens obtidas por

visualização volumétrica direta do modelo plano são iluminadas incorretamente

como pode ser visto na Figura 37. Nesta figura temos imagens obtidas variando

o ângulo de incidência entre a direção da iluminação e a normal dos horizontes


(planos). Ângulo de incidência quase perpendicular aos horizontes (esquerda).

Ângulo de incidência intermediário (centro). Ângulo de incidência quase paralelo

aos horizontes (direita).

Figura 37 – Visualização volumétrica direta com iluminação aplicada ao modelo plano.

A Figura 37 foi obtida utilizando um algoritmo de visualização volumétrica

direta baseado no modelo de iluminação local de Phong. A equação de

iluminação utilizada leva em conta apenas as componentes ambiente e difusa,

por isso podemos observar a variação da tonalização das imagens da esquerda

para a direita. Apenas uma fonte de luz direcional e perpendicular ao plano da

imagem foi considerada no cálculo. Da eq. (13) observamos que os vetores

normais definidos como o gradiente, ou são nulos, ou apontam na direção

vertical (tempo), para cima ou para baixo de acordo com o sinal da função

)(0 tX . Desta forma a imagem à esquerda está mais iluminada (mais clara) pelo

fato dos horizontes estarem quase perpendiculares à direção de iluminação e a

imagem à direita está menos iluminada (mais escura) pelo fato dos horizontes

estarem quase paralelos à direção de iluminação. Na imagem ao centro

observamos linhas escuras coincidindo com a localização dos eventos

(horizontes) do dado. Estes artefatos indesejáveis se devem ao fato mencionado

anteriormente do gradiente ser identicamente nulo nos horizontes do modelo,

pois nestes pontos temos 0)(0 =∂

∂t

tX. Portanto, mesmo para um modelo bem

simplificado de dado sísmico 3D como o modelo plano, verificamos que a

utilização do gradiente para estimar o vetor normal no modelo de iluminação

local nos leva a resultados incorretos.


Ao utilizarmos o gradiente como estimativa para o vetor normal a superfície

em cada voxel cometemos um erro ainda mais grave que o erro apresentado

anteriormente no exemplo do modelo plano. Este erro tem origem na própria

natureza da informação que se deseja extrair do dado sísmico. Os eventos

sísmicos não são necessariamente superfícies de nível do campo escalar de

amplitude sísmica. Esta observação surgiu durante a nossa pesquisa na

utilização da Opacidade 2D. Por isso não podemos garantir que o gradiente

deste campo estará apontando corretamente na direção perpendicular a

superfície dos refletores. Para ilustrar este problema utilizaremos um outro

modelo, derivado do modelo domo, que foi obtido induzindo uma perturbação

radial nas amplitudes do modelo domo.

Seja 1( , , )X t x y o campo escalar de amplitudes sísmicas do modelo domo,

definimos o modelo domo modificado 2 ( , , )X t x y como equação abaixo:

2 21 2( , , ) ( , , ) cos( )X t x y X t x y x yα ω = + (14)

As constantes α e ω são escolhidas de forma a manter a posição espacial

e a geometria dos refletores do dado, porém introduzindo uma variação lateral

nas amplitudes ao longo destes. A Figura 38 mostra, do lado esquerdo, uma

seção transversal (inline) do modelo sintético do domo modificado. A seta indica

um horizonte onde pode ser observada a variação lateral de amplitudes. Do lado

direito da figura vemos a imagem do horizonte indicado obtida por visualização

volumétrica com opacidade 2D. É fácil ver que os horizontes contidos no volume

do domo modificado não são superfícies de nível do campo escalar de amplitude

sísmica.


Figura 38 – Modelo domo modificado. Modelo sintético com variação de amplitudes ao

longo dos refletores.

A Figura 39 mostra o resultado obtido com a visualização volumétrica do

dado do domo modificado utilizando o gradiente para estimar os vetores normais

nos voxels. Do lado esquerdo temos a imagem obtida sem iluminação, e do lado

direito com iluminação usando o gradiente da amplitude para estimar o vetor

normal nos voxels, como é usual. Podemos observar que o efeito de tonalização

causado pela iluminação não está coerente com a geometria da superfície que

forma o horizonte. O padrão gerado pela iluminação sugere um relevo

inexistente na superfície, uma falsa impressão de rugosidade. O efeito

indesejável observado nesse exemplo pode ser comparado à técnica de bump

mapping (Blinn, 1978) onde, é causada uma perturbação nos vetores normais

durante o cálculo da iluminação, por meio de texturas, para dar à superfície uma

aparência áspera, rugosa, ondulada etc. De fato a geometria do horizonte em

questão não possui as ondulações sugeridas pela iluminação, porém os vetores

normais utilizados na iluminação acompanham a geometria das superfícies de

nível do modelo.


Figura 39 – Iluminação aplicada a um horizonte com variação lateral de amplitudes.

A Figura 40 mostra a geometria de duas superfícies de nível que estão

situadas na mesma região do horizonte da Figura 39. A diferença entre a

geometria do horizonte e das superfícies de nível é evidente. Esta é uma razão

ainda mais forte para não utilizarmos o gradiente para estimar os vetores

normais.

2 1( , , )X t x y c= 2 2( , , )X t x y c=

1 2c c>2 1( , , )X t x y c= 2 2( , , )X t x y c=

1 2c c>

Figura 40 – Superfícies de nível do modelo sintético do domo modificado.

Os problemas observados nestes exemplos justificam o fato dos sistemas

de visualização de dados sísmicos que fornecem a opção de visualização

volumétrica não utilizarem iluminação.



Em (Barnes, 2003) é apresentada a técnica de relevo sombreado (shaded

reliefs). Esta técnica produz um atributo sísmico volumétrico que é obtido usando

a equação do modelo de iluminação de Phong. Este atributo é calculado

iluminando todos os voxels do dado, tratando cada voxel como um ponto em

uma superfície cuja orientação local é dada pelo mergulho e pelo azimute.

A abordagem do relevo sombreado é uma técnica comumente utilizada

para realce em mapas e difere da nossa principalmente em 2 pontos. Primeiro

na própria visualização, o autor desse trabalho sugere que o atributo calculado

seja visualizado como uma textura aplicada a fatias horizontais (time-slices) do

dado sísmico ou aplicada a horizontes já rastreados. Portanto, podemos

classificar o relevo sombreado como um método de visualização volumétrica

indireta. A própria maneira como a orientação local dos horizontes é calculada

nos voxel também é diferente. O atributo de iluminação é calculado a partir dos

atributos de mergulho e azimute locais, estes atributos por sua vez são

calculados a partir dos números de onda locais de cada voxel. O autor não

sugere a aplicação o relevo sombreado na visualização volumétrica direta.

4.3. Solução Proposta

Como foi visto anteriormente os eventos sísmicos não são

necessariamente superfícies de nível do campo escalar de amplitude sísmica.

Por isso o gradiente deste campo não é uma boa aproximação dos vetores

normais das superfícies definidas pelos eventos. Para fazer uso de algoritmos de

visualização volumétrica com cálculo de iluminação baseado em modelos locais,

precisamos calcular um campo vetorial, definido no mesmo domínio do dado,

que seja perpendicular aos eventos sísmicos. Em (Silva, 2003) sugerimos uma

maneira de calcular tal campo vetorial para iluminação com bons resultados

práticos tanto para os dados sintéticos quanto para dados reais. Essa nova

maneira de estimar vetores normais em dados sísmicos volumétricos surgiu da

observação, feita nos trabalhos sobre opacidade 2D, de que o atributo de fase

instantânea se mantém constante ao longo de um evento sísmico. Daí nasceu a

idéia de utilizar o gradiente do campo escalar de fase instantânea para

aproximar os vetores normais aos eventos sísmicos.


Calculamos o campo gradiente do atributo de fase instantânea da seguinte

forma:

( , , )X t x y (15)

( , , )Y t x y (16)

Na eq. (15) X é o campo escalar da amplitude sísmica, e na eq. (16) Y é

o campo das transformadas de Hilbert dos traços de X, ou seja:

( , , ) ( , , ) ( )Y t X t h t= ∗i i i i (17)

Cada traça do campo Y pode ser obtido pela convolução do traço

equivalente em X com o núcleo da transformada de Hilbert (Apêndice A), como

na eq. (17). Definimos então o campo escalar de fase instantânea como sendo:

( , , )( , , ) arctan( , , )

Y t x yt x yX t x y

θ

=

(18)

Utilizamos o campo escalar de fase instantânea com valores contidos nos

intervalos [-π , π ]. Devido ao fato do campo escalar de fase instantânea ter

descontinuidades, isto é, saltos de –π para π ao longo dos traços, utilizamos a

expansão da derivada do arco-tangente para definir o campo gradiente da fase

instantânea.

( , , ) arctan , arctan , arctanY Y Yt x yt X x X y X

θ ∂ ∂ ∂ ∇ = ∂ ∂ ∂

(19)

Substituindo a expressão da derivada do arco-tangente,

2

1arctan( )1

d sds s

=+

(20)

Por exemplo, substituindo a eq. (20) na primeira componente do gradiente

da eq. (19) e usando a regra da cadeia temos:


2 2

( , , ) 1arctan

( , , )1

t tXY X YY t x yt X t x y XY

X

− ∂ = ∂ +

(21)

Simplificando,

2 2

( , , )arctan( , , )

t tY t x y XY X Yt X t x y X Y

∂ −= ∂ +

(22)

Finalmente chegamos à expressão do campo gradiente da fase

instantânea abaixo.

( )2 2

1( , , ) , ,t t x x y yt x y XY X Y XY X Y XY X YX Y

θ∇ = − − −+

(23)

4.4. Resultados

Nesta seção apresentaremos os resultados obtidos utilizando nossa

proposta que utiliza o campo vetorial definido na eq. (23) para estimar o vetor

normal às superfícies formadas pelos eventos sísmicos.

Na Figura 41 temos do lado esquerdo o resultado obtido aplicando o

método proposto ao dado sintético do modelo plano, e do lado direito temos a

imagem obtida com a abordagem tradicional. Podemos observar na Figura 41

que o método proposto evita o aparecimento dos artefatos (linhas pretas) na

região onde estão localizados os horizontes. A seta indica artefatos de

visualização presentes na imagem gerada com nossa abordagem, porém

podemos observar que estes artefatos não estão localizados na região de

nenhum horizonte do dado.


-0.95 0 0.95-0.95 0 0.95 Figura 41 – Imagem do modelo plano obtida pelo método proposto (esquerda). Imagem

obtida pelo método tradicional (direita).

A solução proposta corrige o problema da iluminação no caso do modelo

plano porque, por definição, a derivada da fase instantânea é a freqüência

instantânea, e a freqüência instantânea é, naturalmente diferente de zero nos

eventos sísmicos (extremos) do traço gerador do modelo plano, como podemos

verificar na Figura 42. Porém a freqüência instantânea pode ir a zero, ou até

mesmo ser negativa, devido a uma amostragem pobre (aliasing) ou devido a

presença de ruído.

-4

-3

-2

-1

0

1

2

3

4

1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49

Amplitude Sísmica Fase Instantânea Freqüência Instantânea-4

-3

-2

-1

0

1

2

3

4

-4

-3

-2

-1

0

1

2

3

4

1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 491 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49

Amplitude Sísmica Fase Instantânea Freqüência Instantânea Figura 42 – Gráfico da amplitude sísmica, fase instantânea e freqüência instantânea do

traço gerador do modelo plano.


Seguindo a solução proposta definimos o campo escalar de fase

instantânea como sendo:

0( , , ) ( )t x y tθ θ= (24)

Onde )(0 tθ é a fase instantânea (unidimensional) calculada no traço

gerador do modelo plano. Tomando o gradiente do campo escalar de fase

instantânea do modelo plano temos:

0( , , ) ( ),0,0t x y tt

θ θ∂ ∇ = ∂ (25)

0 0( ) ( )f t ttθ∂

=∂

(26)

Portanto pela eq. (26) o gradiente do campo escalar de fase instantânea

está bem definido mesmo na posição onde estão localizados os eventos

sísmicos do modelo plano, graças ao fato de que a derivada da fase instantânea

é a freqüência instantânea )(0 tf .

Examinaremos agora os resultados obtidos com a solução proposta

aplicada ao dado sintético do modelo do domo modificado. Na Figura 43 temos,

do lado esquerdo, a visualização do campo vetorial normal à superfície de uma

esfera mapeado em cores, ou seja, as componentes dos vetores normais foram

convertidas em cores e aplicadas como uma textura em cima da superfície,

segundo as equações abaixo.

( , , )n x y z= (27)

( 1) / 2 (Canal vermelho)R x= + (28)

( 1) / 2 (Canal verde)G y= + (29)

( 1) / 2 (Canal azul)B z= + (30)

Considerando que o campo vetorial esteja normalizado temos que as

componentes dos vetores normais estão no intervalo [-1, 1], conseqüentemente,

após a transformação descrita nas equações (28) a (30), as intensidades dos

canais de cor estão no intervalo [0, 1]. No centro da Figura 43 temos a


visualização do campo vetorial que é o gradiente do campo escalar de fase

instantânea (solução proposta) do dado do modelo do domo modificado,

enquanto do lado direito temos o campo vetorial do gradiente da amplitude

sísmica do mesmo dado. O tom azulado significa que os vetores estão próximos

da direção T (direção vertical) apontando para cima. Observamos que a

geometria da superfície definida pelo horizonte do modelo assemelha-se ao topo

da esfera, por isso vemos que o campo vetorial definido pelo gradiente da fase

instantânea (centro) está coerente com a geometria do horizonte enquanto o

campo gradiente da amplitude sísmica (direita) não reflete a geometria do

horizonte e apresenta um padrão de interferência que é causado nitidamente

pela variação lateral de amplitudes do modelo.

Figura 43 – Visualização do campo vetorial normal. Campo normal à superfície de uma

esfera (esquerda). Campo normal a um horizonte sísmico obtido pela solução proposta

(centro), campo obtido pelo método tradicional (direita).

Na Figura 44 vemos o resultado do método proposto aplicado ao dado

sintético do modelo do domo modificado (imagem ao centro), comparado com a

visualização volumétrica sem iluminação do dado (imagem à direita) e com a

visualização volumétrica com iluminação com cálculo baseado no gradiente da

amplitude sísmica. Podemos observar que a iluminação obtida utilizando o

gradiente da fase instantânea está coerente com a geometria do horizonte. Na

Figura 44 a direção da fonte de iluminação é a direção perpendicular ao plano da

imagem e o modelo local de iluminação está considerando apenas a

componente difusa.


Figura 44 – Visualização volumétrica do dado sintético com variação de amplitudes ao

longo dos horizontes: sem iluminação (equerda), iluminado usando a solução proposta

(centro) e iluminado da maneira tradicional (direita).

Aplicando o método de iluminação proposto a dados sísmicos 3D reais

obtivemos ótimos resultados. Para ressaltar a importância da iluminação

definimos na Figura 45 uma função de transferência com tabela de cor

monocromática e uma função de opacidade 2D que destaca um horizonte

sísmico. Comparando a imagem (a) obtida sem iluminação com a imagem (b)

obtida pelo método proposto, podemos observar que o efeito da iluminação

revela detalhes estruturais do horizonte destacado. Na imagem (c) obtida com

iluminação calculada usando o gradiente da amplitude sísmica para estimar os

vetores normais, podemos observar o mesmo padrão de interferência que foram

observados no dado sintético devido à variação de amplitudes ao longo do

horizonte.

(a) (b) (c)(a) (b) (c) Figura 45 – Visualização volumétrica de um dado sísmico 3D real: sem iluminação

(esquerda), iluminação usando o método proposto (centro) e iluminação usando a

abordagem tradicional (direita).

Na Figura 46 utilizamos para teste um dado real de aquisição marinha. Nas

quatro imagens da figura a cor representa o valor da amplitude sísmica. Na


imagem (a) todos os voxels estão com opacidade igual a um, a parte branca

superior corresponde à lâmina d’água do dado. Na imagem (b) temos a

visualização volumétrica sem iluminação, o que faz com que detalhes da

topologia do fundo do mar sejam perdidos.

(a) (c) (d)(b)

1280-127

(a) (c) (d)(b)

1280-127

(b)

1280-127 Figura 46 – Visualização volumétrica. Adaptado de (Silva et al, 2003).

A Figura 46c foi obtida utilizando o método proposto, ou seja, o gradiente

da fase instantânea foi utilizado para aproximar os vetores normais nos voxels

para o cálculo da iluminação, enquanto na imagem (d) foi utilizado o gradiente da

amplitude sísmica para estimar o vetor normal (abordagem tradicional).


Neste capítulo concluímos que a abordagem tradicional para o cálculo da

iluminação na visualização volumétrica direta, quando aplicada a dados

sísmicos, nos leva a resultados incorretos. O uso do gradiente do dado

volumétrico de amplitude sísmica para estimar o vetor normal (orientação local)

não pode ser aplicado na visualização de horizontes sísmicos. Apresentamos os

problemas da abordagem tradicional aplicando-a a dois modelos sintéticos de

dados sísmicos e examinando os resultados. No caso do modelo mais simples

de camadas planas e horizontais, verificamos que a abordagem tradicional falha

porque o gradiente vai a zero exatamente onde estão localizados os horizontes

sísmicos. No caso do modelo com forma de domo (anticlinal) e com variação

lateral de amplitudes verificamos que a abordagem tradicional falha porque o

gradiente acompanha as superfícies de nível da amplitude sísmica e estas não

coincidem com os horizontes sísmicos.

Também concluímos que a nossa abordagem, usando o campo vetorial do

gradiente da fase instantânea, calculado segundo a eq. (23), é uma maneira


mais apropriada de estimar o vetor normal ao longo dos horizontes sísmicos do

dado. Verificamos isto aplicando a nossa abordagem e a abordagem tradicional

aos dados sintéticos e reais e comparando os resultados. Esta comparação foi

feita baseada em critérios de qualidade visual das imagens resultantes e na

presença (ou ausência) de artefatos de visualização.

Rastreamento de Horizontes Sísmicos 79

5 Rastreamento de Horizontes Sísmicos

Uma técnica importante de interpretação sísmica é o rastreamento de

horizontes. O rastreamento de um horizonte em um dado sísmico 3D consiste

em identificar em quais traços o horizonte esta presente e, nestes traços,

identificar quais amostras pertencem ao horizonte.

As dificuldades encontradas no processo de rastreamento de horizontes

são principalmente a baixa razão-sinal-ruído em algumas regiões do dado, a

proximidade entre vários horizontes distintos, a própria complexidade da

geologia e a presença de falhas geológicas. Grande parte dos algoritmos tende a

confundir horizontes diferentes como se fossem um único durante o

rastreamento. O problema de rastreamento de horizontes torna-se ainda mais

difícil devido ao fato do horizonte sísmico não ser uma superfície de nível dentro

do volume de amplitudes sísmicas, por isso os métodos tradicionais de extração

de superfícies como o Marching Cubes (Lorensen e Cline, 1987) não podem ser

aplicados diretamente.

O rastreamento de horizontes pode ser considerado, do ponto de vista de

visão computacional, como uma técnica de segmentação ou, do ponto de vista

de modelagem e visualização, como uma técnica de extração de superfícies. A

finalidade do rastreamento vai além da simples visualização do horizonte. Mas,

dentro do contexto deste trabalho de tese, consideramos o rastreamento de

horizontes como sendo uma técnica de visualização volumétrica direta adaptada

para dados sísmicos.

Neste capítulo apresentamos um novo método de extração de horizontes

sísmicos. Este método é baseado em uma abordagem de otimização e foi

implementado e testado em dados sintéticos e reais. Na seção 5.1 citamos

alguns trabalhos relacionados com a nossa pesquisa. Na seção 5.2

apresentamos um modelo matemático de otimização para o problema de

rastreamento. Na seção seguinte propomos um método heurístico que produz

soluções quase-ótimas para o problema de otimização. Na seção 5.4

apresentamos os resultados obtidos pela aplicação do nosso método em um

dado real. Na seção 5.5 fazemos considerações sobre os resultados obtidos.



Existem diversas abordagens diferentes para resolver o problema de

rastreamento de horizontes sísmicos de maneira automática ou semi-automática.

(Woodham et al., 1995) sugere uma abordagem probabilística que aplica filtros

de associação de dados probabilísticos (probabilistic data association – PDA) ao

problema. (Alberts et al., 2000) usa redes neurais artificiais para rastrear vários

horizontes simultaneamente. Em (Dorn, 1998) os métodos de rastreamento

automático são classificados em duas categorias: métodos baseados em

feições, que procuram por máximos (mínimos ou cruzamentos de zero) de

amplitude; e métodos baseados em correlação, que usam a medida de

correlação cruzada para guiar o rastreamento de um traço sísmico para o outro.

A nossa abordagem está mais próxima dos métodos baseados em correlação.

5.2. Modelagem de Otimização

Nesta seção apresentaremos uma modelagem matemática onde o

problema de rastreamento de horizontes pode ser considerado como um

problema de otimização. Nesta abordagem um horizonte sísmico discreto pode

ser aproximado pela solução deste problema de otimização. Nas subseções

seguintes introduziremos algumas definições e notações, que serão usadas na

modelagem do problema.

5.2.1. Volumétrico Sísmico Discreto

Para modelar o rastreamento de horizontes sísmicos como um problema

de otimização, consideramos como um dos dados de entrada um conjunto de

amostras (voxels) de amplitudes sísmicas, em valores reais, tomadas numa

grade regular tridimensional U . De acordo com a notação padrão usada em

geofísica representamos o dado de amplitude sísmica por uma função X

definida sobre o espaço U .

:X U → R (31)

3[1, ] [1, ] [1, ]U NT NX NY= × × ⊂ Z (32)


Os valores NT , NX e NY são, respectivamente, o número de amostras

em cada traço, o número de crosslines do dado e o número de inlines do dado.

Dado um voxel ( , , )i j kv v v v U= ∈ os inteiros [1, ]iv NT∈ , [1, ]jv NX∈ e

[1, ]kv NY∈ são os índices na direção vertical (tempo), direção crossline e

direção inline do voxel, respectivamente.

5.2.2. Conjunto de traços do dado

Definimos U como sendo o conjunto de todos os traços sísmicos contidos

no dado volumétrico em questão. Definir o problema de otimização utilizando

este conjunto, e não o conjunto U servirá para garantir a condição de que um

horizonte sísmico intercepta um traço apenas uma vez. Voltaremos à discussão

desta condição mais adiante. Cada traço do dado está completamente

determinado pelo seu par de coordenadas espaciais, por isso definimos:

[1, ] [1, ]U NX NY= × (33)

É importante observar que existe uma relação natural entre as amostras

(voxels) do conjunto U e os traços do conjunto U , i.e., cada amostra está

contida em um único traço. A partir desta relação natural entre estes conjuntos,

que pode ser considerada como uma projeção de U em U , definimos a

notação abaixo.

:

( , , ) ( , )v U Uv i j k v j k

→= → =

(34)

Portanto, dada uma amostra v contida em U usamos a notação v para

representar o traço em U contendo v .


5.2.3. Horizonte Sísmico Discreto

Definimos um horizonte sísmico discreto H como sendo simplesmente um

subconjunto de voxels com a condição de que este conjunto possui apenas um

único voxel por traço. Esta condição pode ser formalizada como abaixo:

, ej j k k i i

H Uu v H v u v u v u

⊂∈ ∴ = = ⇒ =

(35)

Na definição acima temos que um horizonte sísmico discreto H é um

conjunto de voxels (subconjunto de U ) satisfazendo a seguinte condição: se u

e v são dois voxels de H , tais que j jv u= e k kv u= , ou seja, u e v pertencem

ao mesmo traço, temos que i iv u= (u e v possuem o mesmo índice na direção

do tempo), ou melhor, u e v são o mesmo voxel. Definido-o desta forma

garantimos que o horizonte sísmico discreto corta um traço do volume uma única

vez.

A condição de que um horizonte sísmico discreto corta um traço do volume

uma única vez não possui uma versão equivalente para horizontes reais. Em

várias configurações geológicas reais esta condição não seria satisfeita, porém

para a modelagem que estamos introduzindo impor esta condição sobre os

horizontes sísmicos esta condição é essencial, é ela que impõe a maior restrição

para que uma solução do problema de otimização contenha um único horizonte

sísmico. Portanto, como estamos interessados em encontras apenas boas

aproximações dos horizontes sísmicos reais achamos que é mais vantajoso

considerar tal condição.

5.2.4. Conjunto de Sementes

Um outro dado de entrada para o problema é um conjunto de voxels,

chamado de conjunto de sementes, escolhido pelo usuário. Na modelagem do

problema de otimização restringiremos o conjunto de soluções viáveis para o

conjunto de horizontes sísmicos discretos que contenham o conjunto de

sementes.

1{ ,..., }MS s s U= ⊂ (36)


Consideramos que o conjunto de sementes é escolhido de forma

consistente, ou seja, não pode ser vazio, não pode conter duas sementes em um

mesmo traço e deve ser escolhido sobre um mesmo horizonte sísmico.

5.2.5. Função de Similaridade

A função de similaridade associa um número real a cada par de voxels do

dado volumétrico, este número tenta aproximar a probabilidade de dois voxels

estarem no mesmo horizonte sísmico, ou seja, quanto maior for o valor da

função de similaridade para dois voxels maior é a chance destes estarem no

mesmo evento sísmico.

:U Uσ × → R (37)

Em nosso método utilizamos o coeficiente de correlação como função de

similaridade entre os voxels. Dados dois voxels ( , , )i j ku u u u= e ( , , )i j kv v v v=

definimos:

22 2

( , , ) ( , , ) ( , , ) ( , , )( , )

( , , ) ( , , ) ( , , ) ( , , )

w w w

i j k i j k i j k i j kd w d w d w

w w w w

i j k i j k i j k i j kd w d w d w d w

X u d u u X v d v v X u d u u X v d v vu v

X u d u u X u d u u X v d v v X v d v v

σ =− =− =−

=− =− =− =−

+ + − + +

=

+ − + + − +

∑ ∑ ∑

∑ ∑ ∑ ∑2

(38)

Na eq. (38) acima w define duas janelas de correlação com largura 2 1w +

em torno dos voxels u e v . Esta janela está restrita aos traços u e v contendo,

respectivamente, u e v . Esta medida de similaridade é comumente usada para

registro de imagens em sistemas de informação geográficas (Brown, 1992). O

valor da função objetivo está compreendido no intervalo [0,1] , quando o valor

está próximo de 1 os voxels são similares e quando o valor está próximo de zero

não são similares.


5.2.6. Função de Semente Associada

A função de semente associada, como o próprio nome diz, associa um

elemento do conjunto de sementes a cada voxel de um horizonte sísmico

discreto.

:s H S→ (39)

A função de semente associada tem um papel importante na definição da

função objetivo do problema de otimização. Tentamos encontrar o horizonte

sísmico discreto onde a similaridade entre cada voxel e sua respectiva semente

associada é máxima. Escolhemos considerar a similaridade entre o voxel e sua

semente associada, e não a similaridade entre o voxel e seus vizinhos, para

diminuir a propagação de erro, já que a informação das sementes é a informação

segura da qual dispomos.

5.2.7. Problema de Otimização

Usando as definições e notações das seções anteriores modelamos o

problema de rastreamento de horizontes sísmicos como o problema de

otimização descrito abaixo:

• Determinar um horizonte sísmico discreto H e uma função de semente

associada s definida em H que

• Maximiza ( , ( ))v H

v s vσ∈∑

• E satisfaz:

o (Condição das sementes)

O horizonte contém todas as sementes.

S H⊂ .

o (Condição de caminho)

Para todo voxel no horizonte existe um caminho de voxels do

horizonte liga esse voxel a sua semente associada. A variação da


coordenada tempo (coordenada i ) entre os voxels do caminho é

limitada.

1, { ,..., } tal quenv H v v H∀ ∈ ∃ ⊂

1 1, ( ) e ( ) ( )n r i r iv v v s v v v δ+= = − ≤ .

o (Condição das amplitudes sísmicas)

As amplitudes sísmicas dos voxels do horizonte estão contidas

num certo intervalo de amplitudes predefinido.

[ ]min max( , , ) , ( , , ) ,i j k i j kv v v v H X v v v X X∀ = ∈ ∈ .

A condição de caminho e a condição das amplitudes ajudam a evitar que o

horizonte solução do problema contenha voxels que não sejam do horizonte de

interesse. A condição de caminho garante que para cada voxel v de uma

solução viável H existe um caminho 1{ ,..., }nv v de voxels de H ligando o voxel

v a sua semente associada ( )s v . O número real δ restringe a variação das

coordenadas na direção do tempo ( )r iv dos voxels rv do caminho. Com a

condição de caminho uma solução viável é composta de componentes conexas

por caminho. Cada componente conexa corresponde a imagem inversa de uma

semente pela função de semente associada 1( )uH s u−= .

Com a condição das amplitudes sísmicas restringimos os valores das

amplitudes dos voxels de uma solução viável. Por exemplo, se queremos

rastrear um horizonte com amplitudes sísmicas positivas podemos definir

min 0X = e maxX = +∞ . Por isso, também é importante que o intervalo de

amplitudes seja definido de forma coerente com a escolha das sementes, ou

seja, se escolhemos sementes com valores de amplitude positiva o intervalo de

amplitudes deve ser definido de acordo.

5.3. Método Proposto

Apesar de não termos conseguido uma prova formal, acreditamos

fortemente que o problema definido na seção anterior é NP-completo. Porém,

ainda existe a possibilidade de se encontrar uma solução exata para o problema

com complexidade polinomial. Nesta seção apresentamos uma solução

heurística baseada em uma estratégia gulosa para encontrar soluções quase-


ótimas para o problema. A idéia básica do nosso método é a cada passo

expandir a solução viável corrente adicionando novos voxels. Para isto

mantemos dois conjuntos de voxels, um conjunto H contendo os voxels da

solução viável corrente e um outro conjunto H + contendo os voxels que são

candidatos a entrar na solução corrente.

Se considerarmos o conjunto H de traços contendo os voxels do

horizonte H este método assemelha-se a uma estratégia do tipo crescimento de

região (region growing), definido sobre uma imagem U . Porém, como H é uma

superfície mergulhada em um volume, necessitamos de uma definição mais

precisa para o conceito de vizinhança para definir os conjunto de voxels

candidatos a entrar na região.

5.3.1. Inicialização

A solução viável inicial deve ser definida de modo a satisfazer a condição

das sementes. Por isto o horizonte inicial é igual ao conjunto de sementes.

H S= (40)

Definimos a função de semente associada como sendo a identidade do

conjunto H .

( ) ,r r rs s s s H= ∈ (41)

Como estamos assumindo que o conjunto de sementes (voxels) é

escolhido de forma consistente, a solução inicial definida desta forma é

trivialmente viável.

5.3.2. Aprimoramento da Solução Corrente

Para estender a solução viável corrente definimos um conjunto H + de

voxels na vizinhança de H . Estes voxels são candidatos a entrar no horizonte

da solução corrente.


Um voxel v pertence ao conjunto H + se o traço v contendo v não é um

traço de H , existe um voxel u H∈ tal que u e v são traços vizinhos e a

variação entre as coordenadas na direção do tempo de u e v está dentro do

limite. Mais precisamente:

[ ]min max

( ) and ( ) ,

i i

v Hv H u H v N u v u

X v X Xδ+

∈

∈ ⇔ ∃ ∈ ∴ ∈ − ≤ ∈

(42)

Na definição acima ( )N u representa a vizinhança do traço u dentro do

conjunto de traços. Esta vizinhança pode ser definida 4-conectada ou 8-

conectada de forma similar às vizinhanças de pixels em uma imagem.

Pela própria construção do conjunto H + podemos dizer que: para cada

voxel v H +∈ existe pelo menos um voxel u H∈ tal que ( )v N u∈ os traços de

u e v são vizinhos e i iv u δ− ≤ a variação da coordenada na direção do tempo

está dentro do limite. Dentre todos as opções escolhemos o voxel u cujo o valor

( , ( ))v s uσ é máximo. Então definimos para cada v H +∈ :

( , ( ))v v s uσ σ= (43)

Em cada passo do método escolhemos o voxel candidato v H +∈ onde

0vσ > e vσ é máximo. Se não existe um voxel candidato com 0vσ > então o

algoritmo é encerrado. Caso contrário, o voxel candidato v escolhido entra no

horizonte H , o conjunto de voxels candidatos H + é atualizado e definimos a

semente associada de v como sendo:

( ) ( )s v s u= (44)

Usando indução sobre o processo descrito acima para selecionar um novo

voxel v , podemos mostrar facilmente que a solução obtida adicionando um novo

voxel escolhido desta forma a uma solução corrente viável ainda é uma solução

viável.


5.3.3. Convergência

O processo para estender uma solução corrente viável, descrito acima, nos

fornece uma seqüência de soluções viáveis onde o valor da função objetivo é

estritamente crescente. Dada a natureza discreta e finita do problema de

otimização em questão, o valor da função objetivo está naturalmente limitado.

Logo, o método converge.

5.3.4. Análise de Complexidade

Apresentamos abaixo um pseudo-código para auxiliar na análise de

complexidade do método proposto.

/* Inicialização */

1 H S← O(1)

2 Para cada rs H∈ O(1)

3 ( )r rs s s= O(1)

4 Atualiza(H + , rs ) O(log(NXxNY))

/* Seleção do Candidato */

5 { }( , , ) arg max |v vv u v Hσ σ +← ∈ O(1)

/* Estendendo a solução corrente */

6 Enquanto ( v∃ e 0vσ > ) O(1)

7 { }H H v← ∪ O(1)

8 ( ) ( )s v s u= O(1)

9 Atualiza(H + , rs ) O(log(NXxNY))

10 { }( , , ) arg max |v vv u v Hσ σ +← ∈ O(1)

Código 1 – Principal.

Na nossa análise consideramos o tamanho do dado igual a NT NX NY× ×

e que foram definidas M sementes. O laço da linha 2 é repedido M vezes (uma


para cada semente). Como o número de sementes é insignificante comparado

com o número de traços NX NY× podemos considerá-lo como uma constante.

Portanto a complexidade da etapa de inicialização (linhas 2 a 4) é igual a

complexidade da função Atualiza. A complexidade desta função será

analisada mais adiante.

A linha 5 pode ser implementada para ser executada em tempo constante

se utilizarmos uma estrutura de dados do tipo heap para armazenar as triplas

( , , )vv u σ . Estas triplas, definidas como em (43), representam os voxels

candidatos que são os elementos do conjunto H + . O heap representando o

conjunto de voxels candidatos está em ordem decrescente dos valores vσ .

O laço da linha 6 é executado a cada vez que um novo voxel é extraído do

heap H + e entra na solução corrente. O número total de voxels que podem

entrar na solução é no máximo NX NY× pois só podemos ter um único voxel por

traço do dado. Portanto, o número de repetições do laço da linha 6 é O(NXxNY).

Sempre que um novo voxel v entra na solução corrente o conjunto de

voxels candidatos precisa ser atualizado. Isto é feito pela função Atualiza.

Nesta função varremos todos os voxels vizinhos de v e incluímos no heap H +

todos aqueles que se enquadram na definição (42). A função Atualiza também

é responsável por gerenciar elementos repetidos no conjunto H + .

Atualiza(H + ,u)

1 Para cada v t.q. ( ( )v N u∈ e i iu v δ− ≤ ) O(1)

2 Se ( v H∉ e min max( )X X v X≤ ≤ ) O(1)

3 ( , ( ))v v s uσ σ← O(1)

4 {( , , )}vH H v u σ+ +← ∪ O(log(NXxNY))

Código 2 – Função Atualiza.

A função Atualiza tem tempo de execução (log( ))O NX NY× pois

consiste basicamente em uma inserção no heap. Logo, o método como um todo

tem complexidade ( log( ))O N N⋅ onde N NX NY= × é a dimensão espacial

(número de traços) do dado visto em mapa.


5.4. Resultados

O método apresentado na seção anterior foi implementado na linguagem

de programação C++ e usa as primitivas básicas do OpenGL para visualizar os

resultados.

O algoritmo foi implementado usando uma estrutura de dados heap para

manipular de forma eficiente o conjunto de voxels candidatos. Com esta

otimização no código o método leva poucos segundos para encontrar um

horizonte em um dado sísmico volumétricos de tamanho 255 256 256× × rodando

em um Pentium 4 com configuração comum.

Para avaliar a qualidade das soluções geradas pelo método proposto

apresentaremos os resultados obtidos com um dado sísmico público que é parte

de um levantamento real. O tamanho do dado é 255 256 256× × amostras.

Figura 47 mostra duas vistas de um dado sísmico volumétrico. Em cada

imagem podemos ver 3 seções correspondendo as faces de um cubo. As cores

representam os valores de amplitude sísmica dos voxels e são: azul para valores

negativos, branco para valores próximos de zero e vermelho para valores

positivos. As curvas em verde, apontadas pelas setas, marcam o horizonte

sísmico de interesse no dado.

- 0 +- 0 + Figura 47 – Horizonte sísmico em um dado real.


As imagens (a), (b), (c) e (d) da Figura 48 ilustram a evolução do nosso

algoritmo durante o rastreamento da Figura 47. Os parâmetros utilizados para

algoritmo de otimização foram: (i) uma função de similaridade com janela de

correlação de 15 voxels de largura ( 7w = ), (ii) o intervalo [ ]0,+∞ para restringir

os valores de amplitude sísmica, e (iii) o limite máximo para a variação na

coordenada do tempo igual a 1 ( 1δ = ).

Utilizamos um conjunto com 3 sementes para rastrear este horizonte.

Comparando a Figura 48b (50%) com a Figura 48c (75%) podemos observar que

o método encontra um boa solução para tratar a falha geológica presente no

dado. Ao invés de tentar atravessar a falha, o método segue um caminho

contornando-a. Isto se deve ao fato dos voxel na vizinhança da falha terem baixa

correlação com qualquer das sementes, portanto os voxels na vizinhança da

falha só entrarão no horizonte solução após todos os outros voxels com

correlação maior serem incluídos. Esta habilidade do método de conseguir evitar

cruzar uma falha movendo-se ao redor dela por um caminho mais apropriado, se

existe um, é uma característica muito útil.

(a) 25% (b) 50%

(c) 75% (d) 100% Figura 48 – Evolução do algoritmo.

A Figura 49 mostra as três sementes usadas no rastreamento do

horizonte. Cada semente está apontada por uma seta. Os voxels do horizonte

estão coloridos de acordo com sua semente associada. Temos três cores

(vermelho, verde e azul), cada uma correspondendo a uma das três sementes. É


interessante observar que as regiões definidas pelas sementes azul e verde

estão separadas por uma falhas menor situada na porção superior da imagem,

esta falha divide o horizonte em duas partes. Por sua vez, a falha maior que está

bem visível na porção direita Figura 48d, não divide o dado e ambos os lados

estão relacionados com a mesma semente, a verde.

Figura 49 – Sementes e função de semente associada.

Na Figura 50 podemos examinar a estrutura temporal do horizonte

rastreado. Nesta figura as cores representam a coordenada na direção do tempo

dos voxels do horizonte.

Figura 50 – Estrutura temporal do horizonte.


Uma falha geológica pode alinhar horizontes sísmicos distintos tornando o

rastreamento uma tarefa extremamente difícil. A Figura 51 mostra duas seções

verticais do dado (a) inline e (b) crossline. Em ambas temos a falha maior (linha

tracejada preta) e três horizontes “A”, “B” e “C”, de cima para baixo, interpretados

manualmente. O horizonte sísmico “A” é o mesmo mostrado nas figuras

anteriores. Como podemos ver, os horizontes “A” e “B” foram alinhados pela

falha, e o mesmo ocorre com os horizontes “B” e “C”.

BC

AB

C

A

(a) inline

AB

C

AB

C

(b) crossline

Figura 51 – Horizontes diferentes alinhados por uma falha geológica.


Rastreamos o horizonte sísmico “A” usando uma ferramenta comercial de

interpretação sísmica e comparamos o resultado com o horizonte rastreado com

o nosso método. Utilizamos apenas uma semente neste teste. Na Figura 52

temos quatro vistas diferentes do horizonte sísmico “A” que foi rastreado usando

um ferramenta comercial de rastreamento. A curva em verde marca a posição

correta do horizonte, e as setas apontam os locais onde o rastreamento está

incorreto. Podemos ver claramente nas imagens (a), (c) e (d) que o rastreamento

saltou do horizonte “A” para o horizonte “B” que está logo abaixo.

(a) (b)

(c) (d)

Figura 52 – Horizonte “A” rastreado usando uma ferramenta comercial.

O horizonte que foi obtido usando nosso método de rastreamento está

mostrado na Figura 53. Tentamos reproduzir nesta figura as mesmas vistas da

figura anterior. Novamente as curvas verdes marcam a posição correta do

horizonte. Nas áreas apontadas pelas setas podemos verificar que nosso

método rastreou o horizonte “A” corretamente, como esperávamos.


(a) (b)

(c) (d)

Figura 53 – Horizonte “A” rastreado usando nosso método de rastreamento.


Neste capítulo mostramos que o rastreamento de horizontes sísmicos

pode ser modelado como um problema de otimização. Apresentamos também

um método heurístico baseado em uma estratégia gulosa para encontrar

soluções quase-ótimas para o problema de otimização.

Os resultados obtidos em dados sísmicos reais nos mostram que as

soluções quase-ótimas fornecidas pelo nosso método são boas aproximações do

horizonte sísmico de interesse. Testamos o nosso método em diversos dados,

porém estes resultados não puderam ser incluídos neste documento devido a

restrições de segurança da informação da Petrobras.

Conclusões e Trabalhos Futuros 96

6 Conclusões e Trabalhos Futuros

6.1. Conclusões

Concluímos que a opacidade 2D aumenta a habilidade de selecionar os

horizontes sísmicos na visualização volumétrica direta. Utilizamos a amplitude

sísmica como a primeira dimensão na opacidade 2D. Utilizando o atributo de

fase instantânea como segunda dimensão conseguimos isolar os horizontes com

amplitude positivas do restante do dado. Observamos que os horizontes

sísmicos são superfícies de nível do atributo de fase instantânea. Devido às

descontinuidades do atributo de fase instantânea este se mostrou ineficiente

para isolar os horizontes sísmicos negativos. O atributo de fase ajusta idealizado

por nós é uma solução possível para este problema. Com o atributo de fase

ajustada conseguimos isolar tanto os horizontes positivos quanto os negativos

do restante do dado, além disso, conseguimos separar os horizontes positivos

dos negativos. Porém, para individualizar um único horizonte vimos que o

atributo sísmico ideal é o atributo ideal é o atributo de fase desenrolada.

Concluímos que o gradiente da amplitude sísmica não aproxima

corretamente os vetores normais aos horizontes sísmicos. Utilizamos modelos

sintéticos simples para mostrar que este problema se deve ao fato de que um

horizonte sísmico não é uma superfície de nível do volume de amplitudes

sísmicas. Por isso quando aplicamos a visualização volumétrica direta a volumes

de amplitudes sísmicas usando a abordagem tradicional para o cálculo de

iluminação obtemos resultados incorretos. Para resolver este problema

sugerimos a utilização do gradiente do atributo de fase instantânea para

aproximar os vetores normais dos horizontes. Este abordagem se mostrou

visivelmente mais eficiente. Consideramos este resultado como a maior

contribuição da nossa pesquisa até aqui.

Concluímos que o problema do rastreamento de horizontes sísmicos pode

ser modelado matematicamente como um problema de otimização. Acreditamos

que o problema definido é NP-Completo. Para este problema de otimização

apresentamos um método heurístico baseado em uma estratégia gulosa. O

Conclusões e Trabalhos Futuros 97

método proposto por nós fornece soluções quase-ótimas que são boas

aproximações do horizonte de interesse. Mostramos que a complexidade do

algoritmo de rastreamento proposto é ( log( ))O N N⋅ onde N é a dimensão

espacial do dado em questão, ou seja, N NX NY= × é o número de traços

sísmicos do dado.

6.2. Trabalhos Futuros

Sugerimos como trabalho futuro a utilização da opacidade 2D combinada

com outros atributos sísmicos, como por exemplo um atributo de coerência que

seja um indicador de descontinuidades para visualização de falhas.

Sugerimos fortemente a investigação de métodos mais robustos e

eficientes para o desenrolamento da fase instantânea em 3D. Acreditamos que o

atributo de fase desenrolada é de extrema utilidade tanto para visualização

volumétrica direta quanto para a indireta.

Sugerimos uma comparação entre o gradiente da fase instantânea e

outras abordagens para aproximar o vetor normal nos horizontes. Por exemplo,

comparar o gradiente da fase instantânea com os números de onda da

transformada de Fourier 3D.

Sugerimos a utilização da informação de orientação obtida com o gradiente

da fase instantânea para implementar filtros direcionais para dados sísmicos.

Estes filtros podem ser de suavização para eliminação de ruído, ou mesmo filtros

de detecção de bordos para tentar mapear falhas geológicas automaticamente.

Sugerimos a investigação da complexidade do problema de otimização

para rastreamento de horizontes. Encontrar uma prova formal de que o problema

é NP-Completo ou elaborar um método eficiente (polinomial) que encontre

soluções ótimas para o problema.

Sugerimos a extensão da abordagem de otimização para trabalhar com

múltiplos atributos.

Bibliografia 98

7 Bibliografia

ALBERTS, P.; WARNER, M.; LISTER, D. Artificial Neural Networks for Simultaneous Multihorizon Tracking Across Discontinuities, 70th Annual

Meeting: Society of Exploration Geophysicists – SEG, p. 651—653, 2000.

Proceedings.

BAHORICH, M. S.; FARMER, S. L. 3-D seismic discontinuity for faults and

stratigraphic features: The coherence cube. The Leading Edge, v.12, n.10, p.

1053—1058. October 1995.

BARNES, A. E. Theory of 2-D complex seismic trace analysis. Geophysics,

v.61, n.1, p. 264—272, January/February 1996.

BARNES, A. E. The complex seismic trace made simple. The Leading Edge,

v.17, n.4, p. 473—478, April 1998.

BLINN, J. Simulation of Wrinkled Surfaces. In: International Conference on

Computer Graphics and Interactive Techniques (SIGGRAPH’78), 5th, p. 286—

292, August 1978. Proceedings.

BLYTHE, D; McREYNOLDS, T. Advanced Graphics Programming Techniques Using OpenGL, In: SIGGRAPH 2000, Course Notes, No. 32

(2000).

BROWN, L. S. A Survey Of Image Registration Techniques, ACM Computing Surveys, v.24, n.4, p. 325—376, 1992.

DORN, G. A. Modern 3-D Seismic Interpretation, The Leading Edge, v.17, n.9,

p. 1262—1272. 1998.

DREBIN, R. A.; CARPENTER, L.; HANRAHAN, P. Volume rendering, In:

SIGGRAPH ’88, 15th, v.22, p. 65—74, August 1988. Proceedings.

Bibliografia 99

FERNANDO, R.; KILGARD, M. J. The Cg Tutorial: The Definitive Guide to Programmable Real-Time Graphics, Addison-Wesley Pub Co, February 26,

2003.

GERHARDT, A. Aspectos da Visualização Volumétrica de Dados Sísmicos.

1998. Dissertação de Mestrado, Pontifícia Universidade Católica do Rio de

Janeiro, Rio de Janeiro - RJ, 1998.

GERHARDT, A. et al. Two-Dimensional Opacity Functions for Improved Volume Rendering of Seismic Data. In: Congresso da Sociedade Brasileira de

Geofísica, VII, 2001 Salvador, Proceedings.

GHIGLIA, D. C.; PRITT, M. D. Two-Dimensional Phase Unwrapping. Wiley,

NY, 1998.

HADWIGER, M.; KNISS, J. M.; ENGEL, k.; REZK-SALAMA, C. High-Quality Volume Graphics on Consumer PC Hardware, In: SIGGRAPH 2002, Course

Notes, No. 42 (2002).

KINDLMANN, G.; DURKIN, J. W. Semi-Automatic Generation of Transfer Functions for Direct Volume Rendering. In: IEEE Symposium On Volume

Visualization, 1998 North Carolina, p. 79—86. Proceedings.

KINDLMANN, G. Semi-Automatic Generation of Transfer Functions for Direct Volume Rendering. Masters of Science Thesis - Program of Computer

Graphics, Cornell University, Ithaca, NY, 1999.

KNISS, J.; KINDLMANN, G.; HANSEN, C. Interactive Volume Rendering Using Multi-Dimensional Transfer Functions and Direct Manipulation Widgets. In: IEEE Visualization 2001 Conference (VIS 2001), 12th, 2001 San

Diego. Proceedings.

KNISS, J.; KINDLMANN, G.; HANSEN, C. Multi-Dimensional Transfer Functions

for Interactive Volume Rendering. IEEE Transactions on Visualization and Computer Graphics, v.8, n.3, p. 270—285, July/September 2002.

Bibliografia 100

LACROUTE, P.; LEVOY, M. Fast volume rendering using a shear-warp factorization of the viewing transformation. In: SIGGRAPH ’94, pp. 451–458, 1994. Proceedings.

LEVOY, M. Display of Surfaces from Volume Data. IEEE Computer Graphics & Applications, v.8, n.3, p. 29—37, May/June 1988.

LORENSEN, W. E.; CLINE, H. E. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. In: International Conference on Computer

Graphics and Interactive Techniques (SIGGRAPH'87), 14th, v.21, n.4, August

1987, p. 163—169. Proceedings.

MACHADO, M. Segmentação de Dados Sísmicos Via Hyperstack para Visualização, Dissertação de Mestrado, Pontifícia Universidade Católica do Rio

de Janeiro, Rio de Janeiro - RJ, 2000.

MAX, N. Optical models for direct volume rendering. IEEE Transactions on Visualization and Computer Graphics, v.1, n.2, p. 99—108, June 1995.

MEIßNER, M.;HOFFMANN, U.; STRAßER, W. Enabling Classification and Shading for 3D-texture Based Volume Rendering Using OpenGL and Extensions. In: IEEE Visualization, 1999. Proceedings.

PFISTER, H. et al. The Transfer Function Bake-Off. IEEE Computer Graphics and Applications, v.21, n.3, p. 16—22, May/June 2001.

PHONG, B. T., Illumination for Computer Generated Pictures, Communications of the ACM, v.18, n.6, p. 311—317, June 1975.

REZK-SALAMA, C.; ENGEL, K.; BAUER, M.; GREINER, G.; ERTL, T.

Interactive Volume Rendering on Standard PC Graphics Hardware Using Multi-Textures and Multi-Stage Rasterization. In: SIGGRAPH/Eurographics

Workshop on Graphics Hardware, 2000. Proceedings.

ROBINSON, E. A.; TREITEL, S. Geophysical Signal Analysis, Prentice-Hall,

1980.

Bibliografia 101

SHERIFF, R. E. Encyclopedic Dictionary of Exploration Geophysics, 3. ed.,

Society of Exploration Geophysicists – SEG, Tulsa, OK, 1991.

SILVA, P. M.; MACHADO, M.; GATTASS, M. 3D Seismic Volume Visualization, In: Congresso da Sociedade Brasileira de Geofísica, VIII, 2003

Rio de Janeiro, Proceedings.

SPITZER, J. Programmable Texture Blending. In: Game Developers

Conference. 2001. http://www.nvidia.com/developer.

WESTERMANN, R.; ERTL, T. Efficiently using graphics hardware in volume rendering applications. In: SIGGRAPH ’98, pp. 169–178, 1998. Proceedings.

WESTOVER, L. Footprint evaluation for volume rendering, In: International

Conference on Computer Graphics and Interactive Techniques, 17th, 1990,

Dallas, TX, Proceedings.

WITTENBRINK, C. M.; MALZBENDER, T.; GOSS, M. E. Opacity-weighted color interpolation, for volume sampling, In: IEEE symposium on Volume

visualization, 1998, p. 135—142, Proceedings.

WOODHAM, C. A.; SANDHAM, W. A.; DURRANI, T. S. 3-D Seismic Tracking

with Probabilistic Data Association, Geophysics, v.60, n.4, p. 1088—1094. 1995.

YILMAZ, O. Seismic Data Processing, Society of Exploration Geophysicists –

SEG, Tulsa, OK, 1987.


Apêndice A: Transformada de Hilbert

A transformada de Hilbert foi derivada originalmente por Hilbert no contexto

da teoria de funções analíticas.

A transformada de Hilbert de um sinal real unidimensional ( )u t é definida

pela integral

1 ( )( ) u sv t VP ds

t sπ

+∞

−∞

=−∫ (45)

A transformada de Hilbert inversa é dada pela integral

1 ( )( ) v su t VP ds

t sπ

+∞

−∞

−=

−∫ (46)

As integrais definidas na eq. (45) e na eq. (46) são impróprias, pois seus

integrandos não estão definidos para s t= . Por isso as integrais estão definidas

usando o valor principal de Cauchy (denotado por VP ). Tomar o valor principal

de Cauchy significa que as equações devem ser calculadas usando limites como

abaixo:

0,

1 ( ) ( )( ) limt A

AA t

u s u sv t ds dst s t s

ε

εεπ

−

→ →∞− +

= + − −

∫ ∫ (47)

As definições da transformada de Hilbert e sua inversa também podem ser

escritas usando a notação de convolução.

( )( ) ( ) ( )f g t f s g t s ds+∞

−∞

∗ = −∫ (48)


1( ) ( )v t u ttπ

= ∗ (49)

1( ) ( )u t v ttπ

−= ∗ (50)

Ao contrário de outras transformadas, como a transformada de Fourier, a

transformada de Hilbert não muda o domínio da função. Por exemplo, uma

função de uma variável de tempo t é transformada para outra função de mesma

variável de tempo t , enquanto a transformada de Fourier muda uma função com

variável de tempo t para uma função com variável de freqüência f .

A transformada de Hilbert está intimamente relacionada com a

transformada de Fourier definida abaixo:

{ } 2( ) F ( ) ( ) ( ) i ftU f u t f u t e dtπ+∞

−

−∞

= = ∫ (51)

A transformada de Fourier do núcleo 1( )h ttπ

= da transformada de Hilbert

é:

2

( ) sinal( )i fteH f dt i ft

π

π

+∞ −

−∞

= = −∫ (52)

onde

1, 0

sinal( ) 0, 01, 0

ff f

f

+ >= =− <

(53)

Aplicando o teorema que diz que a transformada de Fourier de uma

convolução é o produto das transformadas à eq. (49) temos:

{ } { } { } { }F ( ) F ( ) ( ) ( ) F ( ) F ( )v f u t h t f u f h f= ∗ = ⋅ (54)

Aplicando a transformada de Fourier inversa temos:


{ }1( ) F sinal( ) ( ) ( )v t i f U f t−= − ⋅ (55)

A equação acima representa uma maneira elegante e precisa para o

cálculo da transformada de Hilbert de um sinal real unidimensional discreto.

Usando a transformada de Fourier direta e inversa, este cálculo pode ser

implementado de forma eficiente usando alguma implementação de FFT. Na

nossa implementação utilizamos a equação acima juntamente com a biblioteca

FFTW para implementar a transformada de Hilbert.

Documents

Pedro Mário Cruz e Silva Visualização Volumétrica de ...webserver2.tecgraf.puc-rio.br/~mgattass/teses/2004TesePedroMario.pdf · Todos os direitos reservados. É proibida a reprodução