Uma metodologia para extração de ângulos de reflexão em …objdig.ufrj.br/60/teses/coppe_m/WilsonSouzaDuarte.pdf · 2011. 9. 26. · UMA METODOLOGIA PARA EXTRAC˘AO DE~ ANGULOS

UMA METODOLOGIA PARA EXTRACAO DE ANGULOS DE REFLEXAO

EM PROFUNDIDADE UTILIZANDO MATRIZES DE TEMPO DE TRANSITO

Wilson Souza Duarte

Dissertacao de Mestrado apresentada ao

Programa de Pos-graduacao em Engenharia

Civil, COPPE, da Universidade Federal do

Rio de Janeiro, como parte dos requisitos

necessarios a obtencao do tıtulo de Mestre

em Engenharia Civil.

Orientadores: Webe Joao Mansur

Cleberson Dors

Rio de Janeiro

Fevereiro de 2011



Wilson Souza Duarte

DISSERTACAO SUBMETIDA AO CORPO DOCENTE DO INSTITUTO

ALBERTO LUIZ COIMBRA DE POS-GRADUACAO E PESQUISA DE

ENGENHARIA (COPPE) DA UNIVERSIDADE FEDERAL DO RIO DE

JANEIRO COMO PARTE DOS REQUISITOS NECESSARIOS PARA A

OBTENCAO DO GRAU DE MESTRE EM CIENCIAS EM ENGENHARIA

CIVIL.

Examinada por:

Prof. Webe Joao Mansur, Ph.D.

Dr. Cleberson Dors, D.Sc.

Prof. Jose Antonio Fontes Santiago, D.Sc.

Prof. Carlos Andres Reyna Vera-Tudela, D.Sc.

RIO DE JANEIRO, RJ – BRASIL

FEVEREIRO DE 2011

Duarte, Wilson Souza

Uma metodologia para extracao de angulos de reflexao

em profundidade utilizando matrizes de tempo de

transito/Wilson Souza Duarte. – Rio de Janeiro:

UFRJ/COPPE, 2011.

XX, 140 p.: il.; 29, 7cm.


Cleberson Dors

Dissertacao (mestrado) – UFRJ/COPPE/Programa de

Engenharia Civil, 2011.

Referencias Bibliograficas: p. 128 – 136.

1. Modelagem Sısmica. 2. Extracao de Angulos em

Sısmica. 3. Angulo de Reflexao. 4. Analises de AVA. I.

Mansur, Webe Joao et al. II. Universidade Federal do Rio

de Janeiro, COPPE, Programa de Engenharia Civil. III.

Tıtulo.

iii

Dedico aquela que muito se

esforcou e se esforca para que eu

chegue cada vez mais longe.

Obrigado D. Valdinete Duarte

iv

Agradecimentos

Agradeco a Deus, por permitir chegar ao fim de mais uma jornada.

Agradeco ao meus orientadores Prof. Webe Joao Mansur e, sobretudo, ao Dr.

Cleberson Dors, o qual confiou a mim o tema desenvolvido neste trabalho.

Agradeco em especial ao professor e amigo Valdomiro Neves Lima. Obrigado

professor! Sem a sua ajuda seria muito difıcil concluir os meus estudos.

Agradeco aos meus professores de graduacao que muito contribuiram para a

minha formacao com seus ensinamentos e conselhos: Rosane Ferreira, Angel Ramon

Sanchez Delagado e Carlos Andres.

Aos amigos do LAMEC, manifesto meus sinceros agradecimentos pela companhia

e discussoes. Ressalto que foi um grande prazer esta ao lado de todos voces durante

estes dois anos. Especialmente a Elias da Conceicao, Gilmar Texeira, Raphael Cor-

reia, Edivaldo Junior, Leandro Di Bartolo, Viviane Ferreira, Israel Nunes, Wilian

Jeronimo, Franciane Peters, Pablo Oyarzun, Cid Monteiro, Edmundo Costa e Ro-

drigo Dias.

Obrigado Ivone, por tudo.

Agradeco ao professor da pos-graduacao Josias Silva, o qual forneceu os elementos

basicos essenciais para que eu trabalhase com o tema desenvolvido, relacioando a

modelagem sısmica, ataves da displina Introducao ao Metodo Sısmico.

Agradeco aos colegas de alojamento Raimundo Freire (Neto), Bruno de Souza,

Gustavo Pereira, Abraao Viana, Hugo Machado e Bruno Jucelino pelo companhei-

rismo e amizade ao longo dos meus estudos na UFRRJ.

Meu muito obrigado a minha famılia, que mesmo de longe me motiva a dar

ousados passos na vida.

Minha sincera gratidao e meu muito obrigado a minha namorada Fabıola Gon-

calves, pela compreensao e apoio nos momentos difıceis durante toda esta jornada.

Obrigado CAPES, pelo apoio financeiro!

v

Resumo da Dissertacao apresentada a COPPE/UFRJ como parte dos requisitos

necessarios para a obtencao do grau de Mestre em Ciencias (M.Sc.)



Wilson Souza Duarte

Fevereiro/2011


Cleberson Dors

Programa: Engenharia Civil

Na interpretacao de dados sısmicos, existe um crescente interesse nas chama-

das analises de AVA (Amplitude versus Angulo) que relacionam a amplitude das

ondas refletidas em subsuperfıcie com os angulos de reflexao correspondentes, pois

conhecendo-se estas relacoes pode-se estimar melhor os parametros petrofısicos do

meio. Dentro deste contexto, o presente trabalho visa apresentar uma metodologia

para obter informacoes de angulos de reflexao utilizando o macro-modelo de veloci-

dades e a matriz de tempo de transito associada ao campo de ondas deste modelo.

O procedimento para alcancar tal fim e via aplicacao do operador gradiente tanto

na matriz de tempo de transito quanto no macro-modelo de velocidades. A meto-

dologia proposta e aplicada para meios acusticos e isotropicos, sendo as matrizes de

tempo de transito obtidas via solucao da equacao iconal e pela modelagem baseada

na solucao por diferencas finitas da equacao completa da onda. Como resultados,

apresentam-se neste trabalho curvas de angulos associadas a incidencia ou passagem

da onda direta em cada ponto do refletor em subsuperfıcie. Avaliam-se diferen-

tes esquemas de calculo das derivadas numericas, necessarias para a aplicacao da

metodologia desenvolvida.

vi

Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the

requirements for the degree of Master of Science (M.Sc.)

A METHODOLOGY FOR EXTRACTING REFLECTION ANGLES IN DEPTH

BY USING TRAVELTIME MATRICES

Wilson Souza Duarte

February/2011

Advisors: Webe Joao Mansur

Cleberson Dors

Department: Civil Engineering

On seismic data analysis there is an increasing interest on the so called AVA

analysis (Amplitude versus Angle) that relate the amplitude of reflected waves on

subsurface with the corresponding reflection angles because by knowing these re-

lations one can better estimate the petrophysic parameters of the media. Within

this context, this research aims to present a methodology to obtain reflection angle

informations by using the macro-model of velocities and the traveltime matrix asso-

ciated to the wavefield of this model. The procedure to achieve this goal consists on

the application of the gradient operator on the traveltime matrix as well as on the

macro-model of velocities to obtain the angle in depth. The proposed methodology

is applied to acoustic and isotropic media and the traveltime matrices are obtained

by solving the Eikonal equation and by the solution of the complete wave equa-

tion with finite difference operators. As results, in this research are presented angle

curves associated to the incidence or direct wave path in each point of the reflector

on subsurface. Different numerical schemes are explored to evaluate the derivatives,

since they are necessary to apply the proposed methodology.

vii

Sumario

Lista de Figuras xi

Lista de Tabelas xv

Lista de Sımbolos xvi

Lista de Abreviaturas xx

1 Introducao 1

1.1 Consideracoes Preliminares . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Revisao Bibliografica . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.3 Objetivo e Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.4 Estrutura da Dissertacao . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Geracao de Matrizes de Tempo de Transito Via Metodos Numeri-

cos 11

2.1 Tempo de Transito Via Propagacao do Campo de Ondas . . . . . . . 12

2.1.1 Equacao da Onda Acustica . . . . . . . . . . . . . . . . . . . . 12

2.1.1.1 Aproximacoes Discretas para Derivadas Espaciais e

Temporais . . . . . . . . . . . . . . . . . . . . . . . . 13

2.1.1.2 Equacao da Onda Discretizada por Diferencas Finitas 14

2.1.1.3 Dispersao Numerica e Estabilidade . . . . . . . . . . 15

2.1.1.4 Termo Fonte . . . . . . . . . . . . . . . . . . . . . . 16

2.1.1.5 Tipos de Condicoes de Contorno . . . . . . . . . . . 17

2.1.2 Equacao nao Reflexiva da Onda . . . . . . . . . . . . . . . . . 20

2.1.3 Obtencao do Tempo de Transito . . . . . . . . . . . . . . . . . 20

2.2 Tempo de Transito Via Solucao da Equacao Iconal . . . . . . . . . . 22

2.2.1 Equacao Iconal . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.2.2 Obtencao do Tempo de Transito . . . . . . . . . . . . . . . . . 23

2.2.2.1 Fast Marching Method . . . . . . . . . . . . . . . . . 23

2.2.2.2 Multistencils Fast Marching Method 2D . . . . . . . 27

2.3 Suavizacao do Campo de Vagarosidade . . . . . . . . . . . . . . . . . 30

viii

3 Obtencao de Angulos de Incidencia Utilizando MTT 32

3.1 Conceitos Preliminares . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.1.1 Frente de Onda . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.1.2 Operador Gradiente . . . . . . . . . . . . . . . . . . . . . . . 34

3.1.3 Reflexao e Refracao das Ondas Acusticas . . . . . . . . . . . . 35

3.2 Metodologia Proposta para Obtencao de Angulo de Incidencia . . . . 37

3.2.1 Angulo da MTT . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.2.2 Angulo do Modelo de Velocidades . . . . . . . . . . . . . . . . 37

3.2.3 Composicao do Angulo de Incidencia/Reflexao . . . . . . . . . 38

3.3 Metodos Numericos para Calculo de Gradientes . . . . . . . . . . . . 39

3.3.1 Gradiente Via Operadores de Diferencas Finitas . . . . . . . . 40

3.3.2 Gradiente Via Funcoes de Base Radial . . . . . . . . . . . . . 40

3.3.3 Gradiente Via Spline Cubica Local . . . . . . . . . . . . . . . 41

3.3.4 Gradiente Via Propagacao do Campo de Ondas . . . . . . . . 45

4 Solucoes de Referencia 48

4.1 Conceitos Preliminares . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4.1.1 Princıpio de Fermat e Lei de Snell . . . . . . . . . . . . . . . . 48

4.1.2 Abordagens Tradicionais do Tracado de Raios . . . . . . . . . 51

4.1.2.1 Shooting Method . . . . . . . . . . . . . . . . . . . . 51

4.1.2.2 Bending Method . . . . . . . . . . . . . . . . . . . . 52

4.2 Geracao de Solucoes de Referencia . . . . . . . . . . . . . . . . . . . . 53

5 Exemplos 57

5.1 Analise dos Metodos de Extracao de Gradiente a Partir de Dados

Semi-Analıticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5.1.1 Extracao do Angulo do Modelo de Velocidades . . . . . . . . . 58

5.1.2 Extracao do Angulo da Matriz de Tempo de Transito . . . . . 71

5.1.3 Composicao do Angulo de Incidencia/Reflexao . . . . . . . . . 75

5.2 Extracao de Angulos a Partir de Dados Numericos . . . . . . . . . . . 78

5.2.1 Extracao do Angulo a Partir da MTT Obtida Pela Equacao

da Onda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

5.2.1.1 Influencia do Tamanho de Incremento de Tempo . . 78

5.2.1.2 Influencia da Frequencia de Corte . . . . . . . . . . . 88

5.2.1.3 Influencia do Contraste de Velocidades . . . . . . . . 99

5.2.1.4 Influencia do Grau de Suavizacao do Modelo . . . . . 102

5.2.2 Extracao do Angulo a Partir da matriz de tempo de transito

Obtida Pela Equacao Iconal . . . . . . . . . . . . . . . . . . . 109

5.2.3 Aplicacao para um Modelo com Varias Camadas . . . . . . . 116

5.2.3.1 Angulo de Incidencia Oriundo da Propagacao de Ondas117

ix

5.2.3.2 Angulo de Incidencia Oriundo da Solucao Iconal . . . 121

6 Conclusoes e Trabalhos Futuros 124

Referencias Bibliograficas 128

A Operadores de Diferencas Finitas 137

A.1 Expansao em Serie de Taylor . . . . . . . . . . . . . . . . . . . . . . . 137

A.2 Operadores para a Primeira Derivada . . . . . . . . . . . . . . . . . . 138

A.3 Operadores para a Segunda Derivada . . . . . . . . . . . . . . . . . . 139

x

Lista de Figuras

2.1 Representacao grafica do operador de diferencas finitas com quarta

ordem no espaco e segunda no tempo. . . . . . . . . . . . . . . . . . . 15

2.2 Funcao fonte dada pela segunda derivada da Gaussiana (CUNHA,

1997) para fcorte = 30Hz. . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3 Camada de amortecimento do lado direito da malha ilustrando a re-

ducao do campo de pressao no ponto (i, j). . . . . . . . . . . . . . . . 19

2.4 Ilustracao esquematica do funcionamento do FMM. . . . . . . . . . . 25

2.5 Ilustracao dos pontos envolvidos nos estencis S 1 e S 2 para um ponto

central x. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.6 Ilustracao do efeito produzido pelo processo de suavizacao com nps =

10 pontos por direcao no macro-modelo de velocidades. . . . . . . . . 31

3.1 Frentes de onda definidas pelas isocronas associadas a amplitude ma-

xima correspondente a propagacao do campo de ondas. . . . . . . . . 33

3.2 Orientacao da forma como e obtido o angulo θ. . . . . . . . . . . . . 34

3.3 Raio incidente gerando um raio refletido e transmitido no meio acustico. 35

3.4 Ilustracao do surgimento das head waves. . . . . . . . . . . . . . . . . 36

3.5 Raio de propagacao dado pela direcao do vetor gradiente e isocronas

representando as frentes de onda. . . . . . . . . . . . . . . . . . . . . 37

3.6 Vetor gradiente e suas componentes em um certo ponto do refletor. . 38

3.7 Exemplo de um raio de propagacao incidindo sobre um refletor incli-

nado onde as componentes do vetor gradiente da matriz de tempo de

transito sao positivos. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.8 RBFs tradicionais. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.9 Malha e celula computacional mostrando os pontos que determinam

o polinomio bicubico. . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.1 Ilustracao geometrica do princıpio de Fermat. . . . . . . . . . . . . . 50

4.2 Ilustracao do shooting method. . . . . . . . . . . . . . . . . . . . . . . 52

4.3 Ilustracao do bending method. . . . . . . . . . . . . . . . . . . . . . . 53

4.4 Ilustracao de modelos discretizados para obtencao de dados. . . . . . 55

xi

5.1 Ilustracao de modelos sem e com suavizacao. . . . . . . . . . . . . . . 59

5.2 Resultados na aproximacao do angulo de mergulho de 10◦ com o mo-

delo de velocidades sem suavizacao. . . . . . . . . . . . . . . . . . . . 60


delo de velocidades suavizado com 5 pontos. . . . . . . . . . . . . . . 61


delo de velocidades suavizado com 10 pontos. . . . . . . . . . . . . . 62



5.6 Modelo de velocidades sem suavizacao constituıdo por um refletor

com angulo de mergulho de Φ = 20◦. . . . . . . . . . . . . . . . . . . 64


delo de velocidades suavizado com 5 pontos. . . . . . . . . . . . . . . 65



5.9 Modelo de velocidades submetido ao processo de suavizacao. . . . . . 67

5.10 Resultados na aproximacao do angulo de mergulho para o modelo de

velocidades da figura 5.9(b). . . . . . . . . . . . . . . . . . . . . . . . 68


velocidades da figura 5.9(c). . . . . . . . . . . . . . . . . . . . . . . . 69


velocidades da figura 5.9(d). . . . . . . . . . . . . . . . . . . . . . . . 70

5.13 Matrizes de tempo de transito semi-analıticas obtidas para dois dife-

rentes modelos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.14 Resultados na aproximacao do angulo α para o modelo de velocidades

da figura 5.1(a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.15 Resultados na aproximacao do angulo α para o modelo de velocidades

da figura 5.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

5.16 Resultados na aproximacao do angulo γ de incidencia apos compo-

sicao do angulo α e β sobre o refletor do modelo de velocidades da

figura 5.1(a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

5.17 Resultados na aproximacao do angulo γ de incidencia apos compo-

sicao do angulo α e β sobre o refletor do modelo de velocidades da

figura 5.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5.18 Comparacao entre matrizes de tempo de transito obtidas com dife-

rentes equacoes e criterios. . . . . . . . . . . . . . . . . . . . . . . . . 81

5.19 Resultados na aproximacao do angulo α (em graus) quando a matriz

de tempo de transito e obtida com ∆t = 10−3 s. . . . . . . . . . . . . . 82

5.20 Resultados na aproximacao do angulo α em regioes acima do refletor. 83

xii

5.21 Resultados na aproximacao do angulo α em regioes abaixo do refletor. 84


de tempo de transito e obtida com ∆t = 25 × 10−6 s. . . . . . . . . . . 85

5.23 Resultados na aproximacao do angulo α em regioes acima do refletor. 86

5.24 Resultados na aproximacao do angulo α em regioes abaixo do refletor. 87

5.25 Modelo de velocidades para analise da frequencia de corte, fcorte. . . . 88


de tempo de transito e obtida com fcorte = 30Hz. . . . . . . . . . . . . 90

5.27 Resultados na aproximacao do angulo α em regioes acima do refletor,

utilizando fcorte = 30Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . 91

5.28 Resultados na aproximacao do angulo α em regioes abaixo do refletor,



de tempo de transito e obtida com fcorte = 60Hz. . . . . . . . . . . . . 93






de tempo de transito e obtida com fcorte = 120Hz. . . . . . . . . . . . 96


utilizando fcorte = 120Hz. . . . . . . . . . . . . . . . . . . . . . . . . . 97


utilizando fcorte = 120Hz. . . . . . . . . . . . . . . . . . . . . . . . . . 98

5.35 Resultados na aproximacao do angulo α para um contraste baixo de

velocidades entre as camadas. . . . . . . . . . . . . . . . . . . . . . . 100

5.36 Resultados na aproximacao do angulo α para um contraste alto de

velocidades entre as camadas. . . . . . . . . . . . . . . . . . . . . . . 101

5.37 Resultados na aproximacao do angulo α, para a MTT obtida pela

modelagem sısmica a partir do modelo de velocidades suavizado com

5 pontos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103


para a MTT obtida pela modelagem sısmica a partir do modelo de

velocidades suavizado com 5 pontos. . . . . . . . . . . . . . . . . . . 104




xiii

5.40 Resultados na aproximacao do angulo α, para a MTT obtida pela

modelagem sısmica a partir do modelo de velocidades suavizado com

10 pontos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106







5.43 Resultados na aproximacao do angulo α, para a MTT obtida pelo

MSFM a partir do modelo de velocidades sem suavizacao. . . . . . . 110

5.44 Resultados na aproximacao do angulo α em regioes acima e abaixo

do refletor, para a MTT obtida pelo MSFM a partir do modelo de

velocidades sem suavizacao. . . . . . . . . . . . . . . . . . . . . . . . 111


MSFM a partir do modelo de velocidades suavizado com 3 pontos. . . 112





MSFM a partir do modelo de velocidades suavizado com 5 pontos. . . 114




5.49 Modelo de velocidades utilizado para obtencao das curvas de incidencia.116

5.50 Modelo de velocidades suavizado com 10 pontos, utilizado na obten-

cao da MTT pela modelagem sısmica. . . . . . . . . . . . . . . . . . . 118

5.51 Resultados na aproximacao do angulo γ sobre os refletores 1 e 2. . . . 119


5.53 Modelo de velocidades suavizado com 5 pontos, utilizado na obtencao

da MTT pelo MSFM. . . . . . . . . . . . . . . . . . . . . . . . . . . . 121



xiv

Lista de Tabelas

2.1 Valores da funcao g = g(h) dependente do tipo de orientacao e es-

quema numerico utilizado. . . . . . . . . . . . . . . . . . . . . . . . . 29

2.2 Coeficientes de condicao upwind para ambos os estencis S 1 e S 2. . . . 30

3.1 RBFs tradicionais. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.1 Parametros utilizados na analise de um modelo contendo um refletor

plano-inclinado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

5.2 Parametros utilizados nas modelagens do modelo da figura 5.1(a). . . 79

5.3 Parametros utilizados na modelagem do modelo da figura 5.25. . . . . 88

5.4 Parametros utilizados na modelagem sısmica do modelo de velocida-

des da figura 5.50. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

xv

Lista de Sımbolos

A Amplitude da serie assintotica de ordem zero, p. 22

Ak Amplitudes como parte da serie assintotica, p. 22

C Campo de velocidades, p. 12

Cmax Maior velocidade do modelo, p. 16

Cmin Menor velocidade do modelo, p. 16

D+i, j Operador de diferencas finitas progressivas, p. 23

D−i, j Operador de diferencas finitas atrasadas, p. 23

F Polinomio de terceiro grau, p. 44

G Funcao utilizada na definicao do operador gradiente, p. 34

Gx Derivada na direcao x, p. 34

Gz Derivada na direcao z, p. 34

Itotal Numero total de pontos na horizontal, p. 13

Jtotal Numero total de pontos na vertical, p. 13

MV Media movel, p. 30

Ntotal Numero total de passos temporais, p. 13

P Campo de pressao no domınio espaco-tempo, p. 12

Preg Matriz de amplitude registrada, p. 20

Rr Coeficiente de reflexao, p. 36

Rt Coeficiente de transmissao, p. 36

S 1 Estencil com orientacao do sistema de coordenadas naturais,

p. 27

xvi

S 2 Estencil com orientacao nas direcoes diagonais, p. 27

U Campo de pressao no domınio espaco-frequencia, p. 22

∆t Incremento temporal, p. 13

∆x Espacamento da malha na direcao x, p. 13

∆z Espacamento da malha na direcao z, p. 13

Φ Angulo de mergulho, p. 55

Θ Angulo entre vetores direcionais, p. 28

α Angulo entre o raio incidente e a reta vertical, p. 37

N Numero de pontos na definicao de sub-malha, p. 42

β Angulo entre a reta normal e vertical, p. 38

δ Denota aplicacao da fonte pontual, p. 13

` Termos de polinomio de terceiro grau, p. 44

η Normal externa ao contorno, p. 18

γ Angulo de incidencia, p. 35

γc Angulo crıtico, p. 36

κ Razao entre velocidades, p. 36

λ Comprimento de onda, p. 88

µ Numero empırico que define a quantidade mınima de pontos

da malha por comprimento de onda, p. 16

ω Frequencia angular, p. 22

µ Numero que define a quantidade de amostras temporais que

serao necessarias para que a frente de onda percorra uma dis-

tancia equivalente ao espacamento da malha, p. 16

τ Derivada direcional, p. 27

g Coeficientes de condicao upwind no MSFM, p. 30

−→r Vetores unitarios, p. 27

∂ Operador de derivada parcial, p. 12

xvii

φ Angulo de transmissao, p. 35

ψ Funcao de base radial, p. 40

ρ Inverso do produto da volatilidade e distancia, p. 45

σ Pesos da funcao de interpolacao, p. 40

τ Tempo de transito, p. 22

x Ponto no domınio, p. 12

x f Ponto de aplicacao da fonte, p. 13

θ Angulo entre vetor gradiente e reta vertical, p. 34

u Funcao de interpolacao, p. 40

υ Distancia de determinado ponto ao contorno, p. 20

ε Fator de volatilidade, p. 45

$ Soma dos quadrados da distancia, p. 45

% Parametro de raio, p. 35

ζ Parametro de forma da funcao RBF, p. 41

a Indice de somatorio, p. 42

b Indice de somatorio, p. 42

f Termo fonte, p. 13

fcorte Frequencia de corte, p. 16

fc Frequencia central da fonte, p. 17

fmu Fator multiplicativo, p. 20

f at Fator de amortecimento, p. 20

g Tipo de orientacao e esquema numerico utilizado no MSFM,

p. 29

h Espacamento uniforme da malha, p. 14

i f Posicao da fonte na direcao horizontal da malha, p. 14

j f Posicao da fonte na direcao vertical da malha, p. 14

xviii

k Indice de somatorio, p. 30

m Inclinacao de segmento de reta, p. 44

n Indice de discretizacao temporal, p. 13

namort Numero de pontos da camada de amortecimento, p. 20

nps Numero de pontos de suavizacao, p. 30

q Parametro de forma para expoente da funcao RBF, p. 41

ri j Distancia entre os pontos xi e x j, p. 40

t Tempo, p. 12

td Translacao temporal da fonte no tempo, p. 17

t f Perıodo da funcao Gaussiana, p. 17

w Peso da funcao bicubica, p. 42

x Coordenada horizontal, p. 12

x f Coordenada horizontal da fonte, p. 13

y′

Derivada de polinomio de terceiro grau, p. 44

z Coordenada vertical, p. 12

z f Coordenada vertical da fonte, p. 13

xix

Lista de Abreviaturas

ART Asymptotic Ray Theory (Teoria Assintotica do Raio), p. 4

AVA Amplitude versus Angle (Amplitude versus Angulo), p. 2

AVO Amplitude versus Offset (Amplitude versus Afastamento), p. 3

CMP Common Mid Point (Ponto Medio Comum), p. 8

FIM Fast Iterative Method (Metodo Rapido Iterativo), p. 125

FMM Fast Marching Method (Metodo Rapido de Marcha), p. 6

FSM Fast Sweeping Method (Metodo Rapido de Varredura), p. 7

GMM Group Marching Method (Grupo de Metodo de Marcha), p. 6

HAFMM Higher Accuracy Fast Marching Method (Metodo Rapido de

Marcha de Alta Precisao), p. 6

MEF Metodo dos Elementos Finitos, p. 40

MLSA Moving Least-Squares Approximation (Aproximacao por Mıni-

mos Quadrados Moveis), p. 126

MSFM Multistencils Fast Marching Method (Metodo Rapido de Mar-

cha com Multi-Estencis), p. 7

MTT Matriz de Tempo de Transito, p. 4

RBF Radial Basis Functions (Funcoes de Base Radial), p. 39

RNAs Redes Neurais Artificiais, p. 40

RTM Reverse Time Migration (Migracao Reversa no Tempo), p. 1

xx

Capıtulo 1

Introducao

1.1 Consideracoes Preliminares

A aplicacao da tecnica de imageamento sısmico, tambem conhecida como migra-

cao sısmica, durante a fase de processamento e uma das etapas mais importantes

no processo de exploracao de petroleo, pois permite delimitar/verificar o posiciona-

mento das diversas interfaces (refletores) que formam a geologia da subsuperfıcie.

O seu desenvolvimento nos ultimos anos se tornou ainda mais crucial para produzir

imagens mais acuradas, uma vez que os refletores de interesse estao localizados em

altas profundidades e imersos em meios com grande grau de complexidade, como

aqueles envolvendo flancos salinos, onde pode existir expressivas variacoes laterais

de velocidade e mergulhos bastante acentuados.

Gracas aos esforcos de diversos pesquisadores, atualmente, existem varios meto-

dos de migracao sısmica que podem ser aplicados pela industria petrolıfera para a

geracao de imagens sısmicas bi ou tridimensionais, a depender do volume de dados

e da qualidade que se deseja obter. Dentre os diversos metodos de migracao existen-

tes, pode-se citar como o mais indicado para imagear meios altamente complexos o

metodo denominado por Migracao Reversa no Tempo (RTM - Reverse Time Migra-

tion) introduzido por MCMECHAN (1983), WHITMORE (1983) e BAYSAL et al.

(1983) por produzir imagens sısmicas precisas em relacao aos demais, como a migra-

cao Kirchhoff (SCHNEIDER, 1978; WIGGINS, 1984), por exemplo. A habilidade

da RTM em produzir imagens de boa qualidade e devida ao fato desta ser baseada

na equacao completa da onda, nao admitindo simplificacoes na sua formulacao, o

que permite contemplar os diversos eventos de outras ondas que surgem no decorrer

do campo de propagacao, como multiplas reflexoes, refracoes, difracoes, ondas eva-

nescentes, etc. (ZHANG e SUN, 2008), o que e tratado simplificadamente ou mesmo

negligenciado na maioria das demais metodologias.

Os processos de migracao agem, basicamente, em dois aspectos distintos do dado

1

sısmico: tempo de transito (traveltime), que trata da parte cinematica, e amplitude

que trata da parte dinamica. Para o tempo de transito, que carrega informacoes

sobre o posicionamento das estruturas geologicas do meio e de suas velocidades, a

migracao atua de forma a reposicionar as reflexoes tanto em suas coordenadas de

superfıcie quanto nas coordenadas de tempo ou nas coordenadas de profundidade

relativas a origem das reflexoes. Para a amplitude, que traz informacoes sobre a lito-

logia e fluidos na subsuperfıcie, e mais especificamente informacoes sobre o contraste

das propriedades das camadas que definem a interface de reflexao, a migracao age

de modo a corrigir efeitos de propagacao da onda e fornecer valores de amplitude

que, relativamente, representam o conjunto das propriedades petrofısicas do meio

(SANTIAGO, 2004).

No inıcio, o imageamento visava obter somente a posicao dos refletores, sem se

preocupar com a amplitude associada a reflexao gerada por cada refletor. Fato este

ocorrido devido aos grandes volumes de oleo encontrados naquela epoca nao exigi-

rem grande processamento/interpretacao para serem localizados. Entretanto, com

o aumento da complexidade e profundidade das geologias envolvidas foi necessario

levar em conta a importancia das amplitudes associadas as reflexoes, exatamente

por serem estas capazes de fornecer informacoes a respeito da subsuperfıcie. A via-

bilizacao de tal estudo foi possıvel gracas ao aumento dos recursos computacionais,

que possibilitaram realizar as analises sısmicas com tecnicas mais apuradas, como o

uso da equacao completa da onda.

Um dos estudos mais notaveis sobre as amplitudes e aquele relacionado a cons-

trucao de curvas de refletividade, como Amplitude versus Afastamento (AVO - Am-

plitude versus Offset) e Amplitudes versus Angulo (AVA - Amplitude versus Angle),

ja que as amplitudes estao diretamente ligadas as propriedades litologicas do meio

(VEEKEN, 2007). Em funcao destes estudos, passou a ser crucial preservar as am-

plitudes nas imagens migradas. Caracterizando que o objetivo e obter nao so o

posicionamento correto dos refletores em profundidade como tambem a construcao

das curvas de refletividade, para posterior processo de inversao. Tendo em vista isso,

a determinacao das amplitudes com maior acuracia tem se tornado fundamental para

estudos de AVO e AVA (GUITTON et al., 2007) e o desenvolvimento de tecnicas

de migracao que as fornecam corretamente apos a migracao, em meios complexos,

continua sendo um desafio (DENG e MCMECHAN, 2007).

A geracao de imagens, em profundidade, a partir de dados sısmicos utilizando

a RTM fornece, em geral, o valor do coeficiente de reflexao, a menos das perdas

por transmissao, que sao funcao dos parametros geofısicos (densidade e velocidade

das ondas compressionais e cisalhantes) e do angulo de incidencia (Lei de Snell). O

processo de empilhamento tradicional pos migracao nao leva em conta a influencia

do angulo sob o qual uma imagem foi gerada, o que pode acarretar na perda de

2

precisao na visualizacao dos refletores que estao localizados em altas profundida-

des devido aos poucos registros da energia refletida por estes causados pelas perdas

por transmissao em comparacao aqueles em baixa profundidade. Alem de que ha

o surgimento do efeito “sorriso” de migracao decorrente da reflexao total apos um

certo angulo crıtico de incidencia (devido as head waves1). Logo e extremamente

importante delimitar a regiao, para cada geologia complexa, onde ainda e possıvel

obter imagem com razoavel padrao de qualidade em funcao do angulo. Natural-

mente, tambem e desejavel que o coeficiente de reflexao estimado apos a migracao

seja em funcao apenas dos parametros geofısicos e nao mais dependentes dos angulos

de incidencia/reflexao, pois e atraves dos parametros que e possıvel caracterizar os

meios geologicos nos quais estao inseridos os refletores.

Desta forma, obter informacoes de angulo de reflexao corretas, permite melhorar

a qualidade das imagens obtidas pela migracao possibilitando avaliar as amplitudes

registradas nas imagens de maneira mais coerente.

Alem disto, conhecer o angulo de reflexao corretamente, possibilita tracar as

curvas de refletividade de AVA, realizar filtragem de ruıdos, fazer o empilhamento

de imagens no domınio do angulo, etc. auxiliando inclusive no chamado estudo de

anomalias de amplitude. O estudo das anomalias de amplitudes tem sido bastante

usado na industria petrolıfera como uma ferramenta auxiliar tanto na determinacao

de novas reservas de petroleo e gas quanto na caracterizacao daquelas que estao em

processo de producao. Tradicionalmente, esse estudo e feito atraves da analise de

AVO, ja que, a princıpio, nao e possıvel calcular o angulo correto em cada ponto

de reflexao. Essa analise seria mais coerente se fosse feita em relacao ao angulo de

incidencia e nao ao afastamento fonte-receptor, uma vez que a amplitude refletida e

funcao do angulo de incidencia e dos parametros do meio em cada lado da interface

(PAUL, 2004; RESNICK et al., 1987).

Em geral, as analises de AVA sao feitas aplicando-se uma transformada offset

para angulo em dados AVO. Rotineiramente, a mesma e feita supondo-se que o

modelo e constituıdo de camadas planas horizontais e utilizando-se as velocidades

intervalares oriundas do campo de velocidades (determinado atraves da analise de

velocidades) e a Lei de Snell para obter os angulos em profundidade. Assim, a supo-

sicao do modelo de camadas planas horizontais, aplicada em areas onde o mergulho

das camadas geologicas e acentuado, leva a obtencao de valores de angulos de refle-

xao maiores do que os reais, gerando interpretacoes nao confiaveis da presenca ou

ausencia de hidrocarbonetos nos reservatorios (SANTIAGO, 2004). Embora para

meios suaves, os quais nao apresentam variacoes laterais de impedancia, usar AVO

ou AVA nao faca muita diferenca, pois de certa forma offset e angulo de incidencia

estao relacionados, em meios que apresentam altas complexidades estruturais acima

1Ondas reflatadas sob angulo crıtico (SHERIFF, 2002).

3

ou no entorno dos objetivos, AVO e AVA nao possuem uma relacao simples (SILVA,

2009), e portanto nao conduzem a resultados similares.

O tracado de raios (ray tracing) constitui uma das ferramentas principais dentre

as utilizadas na obtencao dos angulos necessarios para fazer a transformacao de

offset para angulo (RESNICK et al., 1987; VEEKEN e RAUCH-DAVIES, 2006;

VEEKEN, 2007). Pelo tracamento de raios, o angulo e obtido naturalmente, uma

vez que este constitui um dos principais parametros do sistema de equacoes. Assim,

o angulo e determinado ao longo do caminho do raio em profundidade.

A abordagem pelo tracado de raios, baseada na teoria assintotica do raio (ART

- Asymptotic Ray Theory), e possibilitada gracas a divisao da equacao da onda em

duas partes: cinematica e dinamica (CERVENY, 1985). Sendo de enorme interesse,

a parte cinematica, nao so por esta permitir obter os angulos de incidencia, mas

tambem pela geracao dos tempos de transito. Embora apresente a vantagem de

possuir um baixo custo computacional, o tracado de raios produz algumas ambigui-

dades em estruturas complexas quanto a existencia de regioes nas quais os raios se

interceptam (caminhos multiplos - multipath) e/ou mesmo naquelas onde estes nao

chegam (zonas de sombras - shadow zones). Salienta-se que os problemas relacio-

nados aos caminhos multiplos sao possıveis de serem resolvidos aplicando a tecnica

denominada por feixes gaussianos (Gaussian beams) (POPOV, 2002).

O calculo dos tempos de transito e de grande importancia sobre varios aspectos

no processamento sısmico (PESTANA e PIMENTEL, 1996; ZHANG et al., 2009)

podendo ser obtido de varias maneiras, como por exemplo, pela solucao direta da

equacao iconal (eikonal) (SETHIAN e POPOVICI, 1999; VIDALE, 1988) ou a solu-

cao dos sistemas de equacoes derivados desta (BLEISTEIN, 1984; CERVENY, 1987;

CHANG e MCMECHAN, 1986), pela extrapolacao da equacao completa (two-way

wave equation) (BULCAO, 2004; LOEWENTHAL e HU, 1991) ou unidirecional

(one-way wave equation) (CLAERBOUT, 1971, 1985) da onda, etc.

A geracao do tempo de transito esta relacionada a ideia de identificar um certo

evento, tal como a chegada da primeira frente de onda em um dado ponto do macro-

modelo de velocidades. Esta ideia e antiga, sendo que o uso do tempo de transito

associado a onda primaria como condicao de imagem constitui a base dos meto-

dos de migracao, ja que e utilizada para ligar um evento registrado no sismograma

(reflexao) com um ponto em profundidade que o gerou (refletor). A obtencao do

tempo de transito utilizando a extrapolacao da equacao completa da onda baseia-

se em registrar o momento em que a frente de onda incide pela primeira vez em

cada ponto do macro-modelo, assim como as amplitudes associadas, de acordo com

um criterio pre-estabelecido. Os tempos de transito sao armazenados em uma es-

trutura de dados denominada por matriz de tempo de transito (MTT), a qual e

geralmente referida como condicao de imagem (quando se esta realizando a etapa

4

de migracao). Ressalta-se que, mesmo nao sendo o objetivo, o procedimento des-

crito anteriormente tem como principal dificuldade a determinacao do coeficiente de

reflexao/transmissao em funcao do angulo.

A determinacao correta do angulo possibilita a recuperacao das perdas por trans-

missao, ocorridas durante a extrapolacao do campo de ondas oriundas da fonte (pro-

pagacao) e dos receptores (depropagacao) (DENG e MCMECHAN, 2007). Uma vez

que ha uma grande dificuldade de obtencao dos angulos de incidencia/reflexao com

precisao, para modelagens utilizando a equacao da onda, o presente trabalho pro-

poe uma nova metodologia para obte-los baseado na matriz de tempo de transito e

no macro-modelo de velocidades. A metodologia proposta considera o meio como

sendo isotropico com modelagem direta baseada na solucao, por diferencas finitas, da

equacao da onda acustica sem simplificacoes na sua formulacao original. Destaca-se

que a metodologia proposta, pode ser facilmente estendida para o caso da equacao

elastica da onda.

A seguir, faz-se uma breve revisao bibliografica sobre os principais metodos de

obtencao de tempo de transito que estao diretamente relacionados com as matrizes

de tempo de transito utilizadas neste trabalho, bem como as metodologias existentes

para a obtencao dos angulos de incidencia/reflexao.

1.2 Revisao Bibliografica

Uma vez que a metodologia desenvolvida neste trabalho e baseada em matrizes

de tempo de transito, cabe revisar as formas de obte-las. Faz-se-ao duas descricoes

diferentes: a primeira sera aquela onde os tempos de transitos sao obtidos direta-

mente da equacao iconal via esquemas de diferencas finitas e outra e atraves de um

criterio numerico, embutido na modelagem sısmica, capaz de registrar os tempos

durante a extrapolacao do campo de ondas.

Em se tratando da obtencao de tempos de transito atraves da solucao da equacao

iconal, RESHEF e KOSLOFF (1986) foram os primeiros a propor calcular os tempos

de transitos diretamente sobre uma malha regular atraves de uma extrapolacao

plana. Desta forma, usando uma aproximacao local de diferencas finitas de primeira

ordem para a equacao iconal, os autores obtem os tempos para todo o macro-modelo

de velocidades. Porem, devido a forma de extrapolacao utilizada neste esquema, o

mesmo apresenta resultados inadequados, mesmo em modelos homogeneos.

Visando resolver esta questao, VIDALE (1988, 1990) desenvolveu, no caso 2D,

uma extrapolacao retangular e, no caso 3D, em forma de cubo, baseado tambem em

diferencas finitas, o qual simula de forma mais adequada o avanco da propagacao

de frentes de ondas circulares e esfericas, respectivamente. Neste esquema, os tem-

pos de transito obtidos eliminam a necessidade do processo de interpolacao usado

5

comumente nos metodos anteriores, como no tracado de raios (CERVENY, 1987),

por exemplo. Caso se deseje fazer o tracamento de raios, este e feito utilizando o

gradiente do tempo de transito. Embora o esquema do autor funcione bem em uma

variedade de modelos, este nao satisfaz o princıpio da causalidade quando aplicado

em modelos contendo fortes contrastes de velocidades.

VAN TRIER e SYMES (1991) propuseram resolver a equacao iconal atraves da

extrapolacao do gradiente do tempo de transito aplicando o metodo upwind finite

difference (ENGQUIST e OSHER, 1980), ja que as mudancas nas componentes do

gradiente do tempo de transito sao descritas pela lei de conservacao hiperbolica

da mecanica dos fluidos. Para melhorar os problemas de precisao, ZHANG (1991)

desenvolveu esta metodologia fazendo a extrapolacao do tempo de transito em co-

ordenadas polares. Uma extensao do esquema de VAN TRIER e SYMES (1991)

foi feita para o caso 3D por POPOVICI (1991) que, mais tarde, foi modificado por

SCHNEIDER (1995) para reduzir problemas de estabilidade usando malhas adap-

tativas.

Inspirado no metodo de Vidale, PODVIN e LECOMTE (1991) sugeriram um

esquema estavel para a obtencao de tempos de transito baseado no princıpio de

Huygens e na aproximacao de frente de ondas planas. Nesta metodologia e possıvel

obter os tempos ate em modelos contendo fortes contrastes de velocidades.

SCHNEIDER et al. (1992) propuseram uma metodologia para a obtencao dos

tempos de transito que utiliza a tecnica de propagacao dinamica. O metodo e

baseado no princıpio de Fermat e usa tecnicas de calculo simples juntamente com

um esquema de mapeamento sistematico para calcular os tempos em uma malha

regular uniforme. O mapeamento e feito atraves de duas abordagens diferentes;

uma chamada de mapeamento de forca bruta (brute force), a qual obtem os tempos

extrapolando coluna por coluna (ou por linha), e a outra e feita de forma mais

natural, atraves da expansao na forma retangular. O metodo consegue obter os

tempos nos casos em que o campo de velocidade apresente grandes contrastes na

distribuicao de velocidades.

TSITSIKLIS (1995), e mais tarde SETHIAN (1996, 1999b) propuseram uma

nova metodologia de solucao para a equacao iconal, popularmente conhecida por

Fast Marching Method (FMM). A abordagem de Tsitsiklis esta voltada para a ob-

tencao de solucoes aproximadas de problemas com trajetorias de tempos mınimos,

enquanto a de Sethian esta direcionada a solucoes de problemas que incorporam fren-

tes evolutivas, isto e, problemas de propagacao de ondas (VAZ, 2004). Utilizando

a discretizacao upwind de GODUNOV (1959), a ideia do FMM e que a solucao em

cada ponto dependa apenas do ponto vizinho que possui o menor tempo, isto e, a

extrapolacao e feita seguindo o princıpio da causalidade (Huygens). Destaca-se que,

a precisao do FMM e dependente do tamanho do espacamento da malha e da ordem

6

do esquema de diferencas finitas utilizado (SETHIAN e POPOVICI, 1999).

Uma desvantagem do FMM e que ele possui uma complexidade computacional

muito elevada, da ordem de O(Nlog(N)), onde N e o numero de pontos da malha,

ja que, a cada passo do esquema, obtem-se a solucao para um unico ponto. Para

tornar o metodo mais rapido, KIM (2001) propos o Group Marching Method (GMM)

onde a solucao e obtida para um grupo de pontos em um unico passo de tempo, os

quais representam o avanco da frente de ondas. Dessa forma, o autor reduziu a

complexidade computacional para O(N).Para melhorar a precisao do FMM, SETHIAN (1999a) sugeriu o Higher Accu-

racy Fast Marching Method (HAFMM) no qual utilizam-se operadores de diferencas

finitas de segunda ordem. Ainda assim, o metodo continuou apresentando graves

erros no inıcio do processo de extrapolacao devido a grande curvatura das frentes

de onda. Para sanar este problema, ALKHALIFAH e FOMEL (2001) propuseram

aplicar o FMM em coordenadas esfericas (3D) e polares (2D). A implementacao em

coordenadas polares melhorou a precisao da solucao sem, praticamente, nenhuma

perda da eficiencia computacional.

Mais recentemente, HASSOUNA e FARAG (2007) propuseram uma versao me-

lhorada do FMM que permite trabalhar com espacamentos diferentes nas duas di-

recoes x e z (∆x = ∆z)), chamada de Multistencils Fast Marching Method (MSFM)

objetivando melhorar a precisao da solucao nas direcoes diagonais da malha. A solu-

cao da equacao iconal e obtida em cada ponto da malha usando dois estencis, sendo

escolhida como solucao aquela que melhor satisfaz a condicao upwind. Os autores

estenderam esta abordagem para o caso 3D onde seis estencis sao usados para cobrir

26 pontos vizinhos, enquanto que no 2D, dois estencis sao suficientes para cobrir 8

pontos.

Motivado pelos trabalhos de BOUE e DUPUIS (1999) e ZHAO et al. (2000),

ZHAO (2004) apresentou um algoritmo iterativo, chamado de Fast Sweeping Method

(FSM) para resolver a equacao iconal sobre uma malha retangular. O metodo usa

a discretizacao upwind de Godunov e Gauss-Seidel para resolver o sistema discreti-

zado, tendo uma complexidade otima de O(N).Agora, da-se-a uma breve revisao bibliografica dos principais trabalhos relacio-

nados a forma de obtencao da matriz de tempo de transito utilizando um criterio

numerico, incorporado na modelagem sısmica.

LOEWENTHAL e HU (1991) desenvolveram dois procedimentos para obter a

matriz de tempo de transito baseada nos seguintes criterios: a matriz de tempo e

formada com base na amplitude maxima da grandeza associada ao campo (pressao

no caso acustico) de ondas e a outra e baseada na amplitude registrada na pri-

meira quebra2 (criterio do tempo mınimo), mas desconsiderando as ondas frontais

2Registro do primeiro evento de pico maximo do sinal sısmico.

7

(head waves). Os autores mostraram que o criterio da amplitude maxima fornece os

melhores resultados.

Inspirado nos trabalhos de LOEWENTHAL e HU (1991), BULCAO (2004)

propos um novo criterio para a extracao do tempo de transito baseado na amplitude

maxima nas proximidades da primeira quebra, isto e, no first break3 do sinal prove-

niente da fonte. Com o emprego deste novo criterio a matriz de tempo de transito

apresenta um numero bem inferior de descontinuidades para as regioes distantes do

ponto de aplicacao da fonte sısmica.

Na sequencia, segue uma rapida descricao dos principais trabalhos relacionados

a extracao de angulos de incidencia/reflexao em profundidade.

RESNICK et al. (1987) foram os primeiros a discutir sobre o procedimento sim-

plificado de calcular os angulos atraves de uma relacao trigonometrica obtida a partir

de um triangulo retangulo fazendo a suposicao de que todas as camadas sao plano

horizontais. Ele mostrou como obter o angulo quando o refletor estava em posicao

inclinada, contanto que o angulo de mergulho do refletor seja conhecido a priori.

BLEISTEIN (1987), TYGEL et al. (1993) e BLEISTEIN et al. (1999) mostraram

como obter o angulo de incidencia e de mergulho do refletor sem conhecimento ne-

nhum acerca da trajetoria do raio incidente ou das propriedades do refletor, usando

para tanto a Migracao para Afastamento Nulo (MZO - Zero Offset Migration) ba-

seada na migracao Kirchhoff e na tecnica de empilhamento ponderado ao longo de

curvas de difracoes para controle das amplitudes. Como demonstrado por VAS-

QUEZ et al. (2003) e TYGEL et al. (1999), a metodologia e aplicada fazendo um

reagrupamento dos dados obtidos pelo experimento sısmico para o domınio CMP

(Common Mid Point) resultando em graficos de AVO. Na sequencia e aplicada uma

segunda MZO, com pesos distintos da primeira, obtendo-se os angulos de incidencia

e consequentemente os graficos de AVA.

O primeiro trabalho a utilizar o tempo de transito na obtencao dos angulos de

incidencia foi proposto por DENG e MCMECHAN (2007), os quais desenvolveram

uma forma de obter os angulos atraves de um processo de interpolacao na matriz

de tempo de transito. Porem os autores nao exploraram adequadamente as tecnicas

e aproximacoes utilizadas, bem como nao abordaram as questoes relacionadas ao

tratamento de descontinuidades destas matrizes, justificando que se trata de um

problema local para geometrias complexas.

Ressalta-se que, as referencias correspondentes aos metodos de obtencao de gra-

dientes utilizados para a aplicacao da metodologia desenvolvida, serao abordadas ao

longo deste trabalho.

3Primeiro evento considerado como sinal proveniente da fonte sısmica (SHERIFF, 2002).

8

1.3 Objetivo e Metodologia

O objetivo deste trabalho e apresentar uma metodologia desenvolvida para extra-

cao de angulos de incidencia/reflexao em profundidade, fazendo-se uma composicao

de angulos obtidos a partir da aplicacao do operador gradiente no macro-modelo de

velocidades e na matriz de tempo de transito. Utilizam-se duas abordagens diferen-

tes para a geracao das matrizes de tempo, a saber: modelagem sısmica da equacao

completa da onda acustica, por diferencas finitas, e solucao direta da equacao iconal,

mediante o Multistencils Fast Marching Method.

Relativamente a obtencao dos gradientes necessarios para a aplicacao da meto-

dologia proposta, sao adotados os operadores de diferencas finitas de quarta ordem,

a interpolacao atraves das funcoes de base radial e a spline cubica local de Akima.

Ressalta-se que foi adicionalmente desenvolvido um esquema numerico para a gera-

cao do gradiente diretamente da frente de ondas propagada, correspondente aquele

obtido da matriz de tempo de transito, aplicando-se para tanto os operadores de

diferencas finitas de quarta ordem ao campo de pressoes, durante a propagacao do

sinal sısmico emitido pela fonte.

1.4 Estrutura da Dissertacao

Este trabalho esta estruturado como segue:

No capıtulo 2, sao apresentados os metodos para obtencao de tempo de transito

pela equacao completa da onda, atraves de um criterio numerico embutido na

modelagem sısmica, e pela solucao da equacao iconal, via Multistencils Fast

Marching Method. Ao final do capıtulo, e descrito o procedimento de suavi-

zacao do macro-modelo de velocidades, necessario para melhores estimativas

dos gradientes ao longo dos refletores que constituem o modelo.

No capıtulo 3, apresentam-se a metodologia desenvolvida para obter os angulos de

incidencia/reflexao em profundidade e os metodos utilizados para calculo dos

gradientes. Inicialmente, e feita uma descricao sobre frentes de onda e como

ocorre o fenomeno de particao de energia, quando esta incide sobre um refletor

contendo diferentes propriedades acusticas. No desfecho do capıtulo, e mos-

trado uma outra abordagem, desenvolvida adicionalmente neste trabalho, para

geracao das curvas de angulos, correspondentes aquelas oriundas da matriz de

tempo de transito, aplicando o operador de diferencas finitas de quarta ordem

diretamente sobre o campo de pressoes durante a propagacao.

No capıtulo 4, e exposto o procedimento empregado para a geracao de solucoes de

referencia (semi-analıticas), que sao utilizadas a tıtulo de comparacao com as

9

solucoes obtidas pelos metodos de gradientes. Sendo o procedimento empre-

gado baseado no tracamento cinematico de raios, e feita uma breve descricao

de duas metodologias tradicionais da teoria de raios, o shooting method e o

bending method. Ainda, com o objetivo de descrever estas duas abordagens,

e presentado, inicialmente, o princıpio de Fermat e a Lei de Snell, necessarios

para compreender os demais conceitos desenvolvidos no capıtulo.

No capıtulo 5, sao apresentados os resultados e respectivas analises da metodolo-

gia desenvolvida neste trabalho para os casos de macro-modelos velocidades

que possuem solucoes de referencia, visando avaliar a precisao dos resulta-

dos numericos obtidos. Destaca-se que e avaliada a influencia dos principais

parametros da modelagem sısmica e do modelo de velocidades que podem

contribuir significativamente na melhor estimativa das curvas de angulos de

incidencia/reflexao em profundidade.

No capıtulo 6, sao apresentados as conclusoes e alguns comentarios, bem como as

sugestoes de trabalhos futuros, envolvendo a metodologia desenvolvida neste

trabalho.

Ao final deste trabalho, no apendice A, e descrito a forma de obtencao dos opera-

dores de diferencas finitas, que sao utilizados ao longo deste trabalho, atraves

da serie de Taylor.

10

Capıtulo 2

Geracao de Matrizes de Tempo de

Transito Via Metodos Numericos

O calculo dos tempos de transito e de grande importancia sobre varios aspectos

no processamento sısmico. Por exemplo, os metodos de migracao e modelagem

Kirchhoff requerem o calculo da funcao de Green, que por sua vez depende do calculo

do tempo de percurso (tempo de transito) entre as posicoes de registro na superfıcie

e os pontos em profundidade no modelo de velocidade (PESTANA e PIMENTEL,

1996).

Da mesma forma, na RTM que pode ser entendida como um problema de con-

dicao de contorno associado a imposicao de uma condicao de imagem, o tempo de

transito e o responsavel por permitir a geracao da imagem em profundidade. Neste

caso, a matriz de tempo de transito contem simplesmente o tempo gasto pelo campo

de ondas para percorrer a distancia entre o ponto de aplicacao da fonte e cada ponto

da malha em profundidade (tempo da onda direta).

Observando-se um evento (frente de onda) de maneira isolada, a solucao da equa-

cao da onda que governa a propagacao do evento possui duas diferentes propriedades

basicas, uma propriedade cinematica e outra dinamica. A informacao cinematica,

que esta associada ao tempo de propagacao, e utilizada no imageamento para li-

gar um evento registrado no sismograma (reflexao) com um ponto em profundidade

que o gerou (refletor), sendo descrita pelas trajetorias de propagacao (raios) e pelos

tempos de transito. A parte dinamica trata da distribuicao espacial da energia pro-

pagada, ou seja, traz informacoes sobre amplitudes, que estao relacionadas com os

coeficientes de reflexao dos refletores.

Em tecnicas como o tracado de raios (BLEISTEIN, 1984; CERVENY, 1987;

POPOV, 2002), que e baseado na solucao da equacao da onda por ondas planas,

as parcelas cinematica e dinamica da solucao podem ser obtidas separadamente

atraves da solucao da equacao iconal e de transporte, respectivamente. Assim, para

obter a matriz de tempo de transito nestes casos, basta resolver a equacao iconal do

11

problema relacionado, como sera visto na secao 2.2. Ja quando se adotam metodos

de discretizacao da equacao da onda (como o MDF ou MEF) associadas a tecnicas

discretas de marcha no tempo, a matriz de tempo de transito e geralmente obtida

aplicando-se as chamadas condicoes de imagem, que nada mais sao do que criterios

baseados na amplitude do campo de ondas propagante, como sera visto na secao

2.1.3.

No presente trabalho, alem da matriz de tempo de transito possibilitar a cons-

trucao de imagens, motivo pelo qual a mesma e associada a chamada condicao de

imagem, ela se constitui como um dado essencial para a geracao de informacoes de

angulo de incidencia em profundidade para cada ponto do macro-modelo de veloci-

dades. Portanto, como ficara mais evidente no decorrer deste capıtulo, a precisao

na obtencao do angulo dependera diretamente da qualidade da matriz de tempo de

transito associada, uma vez que na metodologia proposta o angulo de incidencia

vira, em parte, do gradiente numerico da matriz de tempo de transito. Neste traba-

lho, duas tecnicas distintas serao utilizadas para se obter as matrizes de tempo de

transito, a saber: propagacao de ondas via o metodo das diferencas finitas e solucao

da equacao iconal.

Neste capıtulo serao apresentados os metodos para obtencao de tempo de tran-

sito pela equacao completa da onda, atraves de um criterio numerico embutido na

modelagem sısmica, e pela solucao da equacao iconal, via MSFM. Ao final, e descrito

o procedimento de suavizacao do modelo de macro-velocidades.

2.1 Tempo de Transito Via Propagacao do

Campo de Ondas

2.1.1 Equacao da Onda Acustica

A equacao da onda acustica utilizada neste trabalho e aquela ja amplamente

difundida e utilizada para a modelagem e migracao sısmica, dada pela equacao (2.1)

para o caso 2D (BERKHOUT, 1987; BULLEN e BRUCE, 1985).

∂2P(x, t)∂x2 +

∂2P(x, t)∂z2 −

1C2(x)

∂2P(x, t)∂t2 = 0, (2.1)

onde P(x, t) representa o campo de pressao (pressao acustica) no ponto x = (x, z)no instante de tempo t do meio descrito pelo campo discreto de velocidades de

propagacao C(x).Para o problema geofısico existe ainda a necessidade de se prescrever o termo

12

fonte, denotado por f (t)δ(x − x f )δ(z − z f ) resultando na equacao (2.2).

∂2P(x, t)∂x2 +

∂2P(x, t)∂z2 −

1C2(x)

∂2P(x, t)∂t2 = f (t)δ(x − x f )δ(z − z f ), (2.2)

onde (x f , z f ) = x f representa o ponto de aplicacao da fonte sısmica e δ o delta de

Dirac. Para resolver a equacao (2.2) numericamente, deve-se discretiza-la no espaco

e no tempo, sendo as coordenadas discretas:

x −→ i∆x,

z −→ j∆z,

t −→ n∆t.

Dessa forma:

P(x, t) −→ Pni, j,

f (t) −→ f n.

Onde:

i = 1, 2, ..., Itotal,

j = 1, 2, ..., Jtotal,

n = 1, 2, ...,Ntotal,

e ∆x, ∆z e ∆t sao os respectivos espacamentos da malha na direcao x, z e temporal.

Isto permitira escrever as demais derivadas existentes na equacao (2.2) de forma

discreta atraves dos operadores de diferencas finitas, conforme sera apresentado na

secao 2.1.1.1. Destaca-se que no caso da geofısica, as condicoes iniciais sao dadas

por:

P0i, j = 0, (2.3)(

∂P∂t

)0

i, j= 0 (2.4)

e as condicoes de contorno sao dadas conforme descrito na secao 2.1.1.5.

2.1.1.1 Aproximacoes Discretas para Derivadas Espaciais e Temporais

A obtencao dos operadores de diferencas finitas pode ser feita tanto por expansao

em serie de Taylor como por interpolacao polinomial (FORTUNA, 2000), sendo que

a tecnica mais utilizada nos processos de modelagem sısmica e aquela baseada na

13

serie de Taylor por ser considerada uma ferramenta matematica basica na definicao

de aproximacoes para as derivadas.

Entao, utilizando quarta ordem de precisao para as aproximacoes das derivadas

espaciais e segunda ordem para a derivada temporal de P, chega-se a (conforme

apendice A):

∂2P(x, t)∂x2 = (Pxx)n

i, j =1

12(∆x)2

[−Pn

i−2, j + 16(Pn

i−1, j + Pni+1, j

)− 30Pn

i, j − Pni+2, j

],(2.5)

∂2P(x, t)∂z2 = (Pzz)n

i, j =1

12(∆z)2

[−Pn

i, j−2 + 16(Pn

i, j−1 + Pni, j+1

)− 30Pn

i, j − Pni, j+2

], (2.6)

∂2P(x, t)∂t2 = (Ptt)n

i, j =1

(∆t)2

(Pn−1

i, j − 2Pni, j + Pn+1

i, j

). (2.7)

A justificativa para a adocao destes operadores discretos vem dos inumeros tra-

balhos encontrados na literatura (ALFORD et al., 1974; BULCAO, 2004; CUNHA,

1997; FARIA, 1986; MUFTI, 1990; SILVA, 2009) indicando uma boa relacao entre

a precisao dos resultados e o custo computacional envolvido, comparativamente a

outras aproximacoes.

2.1.1.2 Equacao da Onda Discretizada por Diferencas Finitas

Fazendo a substituicao das equacoes (2.5), (2.6) e (2.7) em (2.2), explicitando-se

o campo de pressoes no tempo n + 1 em funcao dos tempos anteriores e fazendo-

se ∆x = ∆z = h, obtem-se a equacao da onda na versao discretizada, referente ao

equilıbrio em um ponto (i, j) da malha computacional, dada por:

Pn+1i, j = −

112

(Ci, j∆t

h

)2 [Pn

i−2, j + Pni+2, j + Pn

i, j−2 + Pni, j+2

−16(Pn

i−1, j + Pni+1, j + Pn

i, j−1 + Pni, j+1

)+ 60Pn

i, j

]+2Pn

i, j − Pn−1i, j + f nδi,i f δ j, j f . (2.8)

A equacao (2.8) implica que o (n + 1)-esimo passo de tempo depende somente

dos dois tempos anteriores (n)-esimo e (n − 1)-esimo, conforme indicado na figura

2.1, ou seja, esta equacao representa um esquema de marcha recursivo explicito

no tempo. Desta forma, possui ∆t crıtico e, consequentemente, requerendo um

criterio de estabilidade. Contrariamente aos esquemas implıcitos de marcha, que

possibilitam ∆t teoricamente ilimitados. Optou-se por utilizar o esquema explıcito,

uma vez que este e frequentemente adotado para problemas de propagacao de ondas

14

com relativo baixo custo computacional, comparativamente aos esquemas implıcitos

que geram um acoplamento entre equacoes onde o custo computacional e espaco para

armazenamento matricial podem ser muito maiores do que aqueles para esquemas

explıcitos devido a necessidade de solucao do sistema de equacoes envolvido (DORS,

2007).

Figura 2.1: Representacao grafica do operador de diferencas finitas com quartaordem no espaco e segunda no tempo.

2.1.1.3 Dispersao Numerica e Estabilidade

Conforme FORTUNA (2000), a dispersao numerica ocorre quando o metodo

numerico faz com que ondas, de diferentes comprimentos, propaguem-se com velo-

cidades distintas, conduzindo a presenca de oscilacoes na solucao computacional.

Ainda, de acordo com FORTUNA (2000), um metodo numerico estavel e aquele

no qual quaisquer perturbacoes ou erros na solucao nao sao amplificados sem limite.

Os exemplos mais diretos para essas perturbacoes e erros sao condicoes de fronteira

ou iniciais aproximadas de forma incorreta, e acumulo de erros de arredondamento

cometidos durante os calculos. O primeiro caso e evitado com uma correta discreti-

zacao das condicoes auxiliares; o segundo nao pode ser evitado, devendo, portanto,

ser controlado. Os erros de arredondamento nao devem crescer sem limites, a ponto

de influir desastrosamente na solucao numerica. O acumulo desses erros pode ser evi-

tado seguido-se os criterios de estabilidade dos metodos numericos, isto e, condicoes

que garantam que o metodo numerico seja estavel.

Para se evitar a dispersao numerica, deve-se limitar o espacamento entre os

pontos da malha de tal forma que este seja pequeno comparado com o comprimento

15

de onda gerado pela fonte. Para o caso da equacao (2.8), isso e conseguido atraves

da relacao (2.9) (BULCAO, 2004; MUFTI, 1990):

h ≤Cmin

µ fcorte, (2.9)

onde h e o espacamento da malha uniforme, Cmin e a menor velocidade do modelo,

fcorte e a frequencia de corte da fonte ou frequencia maxima utilizada e µ e um

numero empırico que determina quantos pontos da malha serao empregados para

representar o menor comprimento de onda.

Quanto ao intervalo de tempo ∆t, em geral metodos explıcitos apresentam um

limite superior no valor do mesmo que deve ser utilizado para garantir a estabilidade

do metodo. Este e, normalmente, expresso em funcao dos valores de ∆x, ∆z e da

velocidade de propagacao C. Assim, a condicao de estabilidade e dada pela expressao

(2.10) (BULCAO, 2004; FARIA, 1986; MUFTI, 1990):

∆t ≤h

µCmax, (2.10)

onde Cmax e a maior velocidade do modelo e µ e um numero que define quantos inter-

valos de tempo serao necessarios para que a frente de onda percorra uma distancia

equivalente ao espacamento entre os pontos da malha. Para o operador dado pela

equacao (2.8), os valores dos parametros µ e µ sao 5 e 4, respectivamente (BULCAO,

2004).

Destaca-se que, o criterio dado por (2.10) deve ser respeitado para que nao ocorra

instabilidade numerica no processo de marcha no tempo.

2.1.1.4 Termo Fonte

Na secao 2.1.1 foi incorporado um termo fonte a equacao da onda na posicao

x f , cujo objetivo e iniciar a propagacao no tempo e no espaco. Com o proposito de

simular uma fonte sısmica real, este termo fonte precisa ter algumas caracterısticas,

como ser preferencialmente limitado, tanto no domınio do tempo, para simular uma

fonte sısmica do tipo explosiva, quanto no domınio da frequencia, para que se tenha

um controle sobre a frequencia maxima ao qual o modelo numerico esta sujeito

(BULCAO, 2004; COSTA, 2006). Como pode ser observado na expressao (2.9), a

frequencia influencia diretamente no grau de refinamento da discretizacao da malha

computacional.

O termo fonte utilizado neste trabalho para a geracao do sinal sısmico e dado

16

pela derivada segunda da Gaussiana (CUNHA, 1997), como segue:

f (t) =[1 − 2π (π fctd)2

]e−π(π fctd)2

, (2.11)

onde td = n∆t − t f representa uma translacao temporal da fonte no tempo, e sendo

t f o perıodo da funcao Gaussiana, dado por:

t f =2√π

fc(2.12)

e fc a frequencia central da fonte, dada por:

fc =fcorte

3√π. (2.13)

A figura 2.2 ilustra um exemplo da funcao fonte no domınio do tempo para

fcorte = 30Hz.

Figura 2.2: Funcao fonte dada pela segunda derivada da Gaussiana (CUNHA, 1997)para fcorte = 30Hz.

2.1.1.5 Tipos de Condicoes de Contorno

As denominadas condicoes de contorno (ou de fronteira), sao definidas nas bor-

das das regioes que delimitam a malha computacional na qual se deseja resolver a

equacao (2.2).

A condicao de contorno mais comum em levantamentos sısmicos e a chamada

condicao de contorno essencial ou de Dirichlet, e corresponde a prescrever a pressao

Pni, j = 0, (2.14)

17

no respectivo contorno. Esta condicao e utilizada para dados marıtimos simulando

as condicoes de reflexao na superfıcie do mar ou em dados terrestres (CUNHA, 1997),

simulando a superfıcie do solo.

Outra condicao de contorno consiste em especificar o valor da derivada, ou da

taxa de modificacao, de P na direcao normal ao contorno, sendo esta chamada de

condicao de contorno natural ou de Neumann, dada por:(∂P∂η

)n

i, j= 0, (2.15)

onde η corresponde a normal externa ao referido contorno.

Assim, fica caracterizado que e possıvel a ocorrencia no mesmo problema de

condicoes de contorno onde P e prescrito em uma parte do contorno e a sua derivada

normal especificada em outra parte, isto e, o problema pode apresentar condicoes

de fronteira mistas.

Outra condicao de contorno bastante utilizada em modelagem geofısica, devido

a introducao de contornos artificiais, e aquela denominada por condicao de contorno

nao-reflexiva ou condicao de contorno absorsora (REYNOLDS, 1978). Esta e uti-

lizada quando se deseja simular problemas com domınios infinitos ou semi-infinitos

visando garantir um correto truncamento do modelo numerico sem o surgimento de

reflexoes artificiais nos respectivos bordos. Assim, pode-se dizer que as condicoes

de contorno nao-reflexivas visam absorver a energia da frente de ondas quando esta

incide no contorno.

Cabe destacar que, as reflexoes indesejadas geradas pelos contornos artificiais

sao muito problematicas para o resultado da modelagem, pois interagem com as

reflexoes provenientes dos eventos em subsuperfıcie, afetando a qualidade do dado

registrado (SILVA, 2008).

O modo mais simples de se derivar uma condicao de contorno nao-reflexiva e

atraves da fatoracao da equacao da onda unidimensional, tal que

∂2P∂x2 −

1C2

∂2P∂t2 =

(∂

∂x+

1C∂

∂t

) (∂

∂x−

1C∂

∂t

)[P] = 0. (2.16)

O primeiro termo do segundo membro da equacao (2.16) e relativo a ondas viajando

no sentido negativo do eixo x e o segundo termo, no sentido positivo. Assim, para

eliminar as reflexoes no contorno esquerdo e direito basta empregar, respectivamente

as equacoes (2.17) e (2.18) abaixo:(∂

∂x+

1C∂

∂t

)[P] = 0 (2.17)

18

(∂

∂x−

1C∂

∂t

)[P] = 0. (2.18)

Discretizando as equacoes (2.17) e (2.18) usando os operadores de diferencas

finitas de primeira ordem no espaco e no tempo, indicados no apendice A, tem-se

que:

Pn+11, j = Pn

1, j +∆tC1, j

∆x

(Pn

2, j − Pn1, j

), (2.19)

Pn+1Itotal, j = Pn

1, j −∆tCItotal, j

∆x

(Pn

Itotal, j − PnItotal−1, j

). (2.20)

Vale ressaltar que as equacoes (2.19) e (2.20) sao formuladas via suposicao de

que a direcao de propagacao da onda e normal ao contorno. Logo, para aplicar

adequadamente a condicao nao-reflexiva para problemas 2D e 3D, seria necessario

decompor o campo de ondas em sua parte normal e tangencial, de modo a absorver

somente a parcela normal do mesmo.

Outra forma de se reduzir as reflexoes espurias geradas pelo truncamento dos

contornos, e aplicar uma camada de amortecimento numerico para reduzir a inten-

sidade da onda antes desta incidir sobre o contorno (CERJAN et al., 1985).

O procedimento e feito aplicando a expressao (2.21) a cada ponto da malha per-

tencente a camada de amortecimento (figura 2.3) durante a propagacao do campo

de ondas. A reducao da amplitude no ponto (i, j) e feita gracas aos fatores multipli-

cativos gerados pela expressao (2.22).

Figura 2.3: Camada de amortecimento do lado direito da malha ilustrando a reducaodo campo de pressao no ponto (i, j).

Pni, j = fmu(υ)Pn

i, j, (2.21)

sendo

fmu(υ) = e−[ f at(namort−υ)]2

, (2.22)

19

onde fmu e o fator multiplicativo para atenuar o campo de pressao, υ e um ındice

que denota a distancia de determinado ponto em relacao ao contorno do modelo,

f at e o fator de amortecimento referente a intensidade da reducao de amplitude

da pressao que se deseja e namort e o numero de pontos que constitui a camada de

amortecimento.

2.1.2 Equacao nao Reflexiva da Onda

Uma alternativa a equacao da onda acustica tradicional (2.1) que visa melhorar

a qualidade da matriz de tempo de transito e a chamada equacao da onda acustica

nao-reflexiva (two way nonreflecting wave equation) (BAYSAL et al., 1984) que pode

ser expressa via equacao (2.23) assumindo impedancia constante ao longo do meio.

C(x)∂

∂x

(C(x)

∂P(x, t)∂x

)+ C(x)

∂

∂z

(C(x)

∂P(x, t)∂z

)−∂2P(x, t)∂t2 = 0. (2.23)

Em meios homogeneos a equacao (2.23) se torna identica a equacao acustica

tradicional, (2.1), ja que a primeira e obtida via modificacao da segunda. Para

meios heterogeneos, BAYSAL et al. (1984) mostrou que a equacao nao-reflexiva

tem a vantagem de fornecer o coeficiente de reflexao igual a zero para incidencia

normal e menor aquele gerado pela equacao reflexiva para incidencias diferentes de

zero. Portanto a equacao acustica nao-reflexiva possibilita atenuar as reflexoes nas

proximidades do refletor, o que pode, a princıpio, melhorar as matrizes de tempo de

transito obtidas, diminuindo as descontinuidades nestas.

2.1.3 Obtencao do Tempo de Transito

Utilizou-se neste trabalho o criterio da amplitude maxima desenvolvida por LO-

EWENTHAL e HU (1991) e o criterio da amplitude maxima nas proximidades da

primeira quebra (first break) proposto por BULCAO (2004). Destaca-se que uma

grande vantagem deste ultimo criterio, e que as matrizes de tempo de transito apre-

sentam um comportamento mais suave (contınuo), para regioes distantes do ponto

de aplicacao da fonte.

Segue abaixo o Algoritmo 01 referente ao criterio da amplitude maxima. Basica-

mente, neste algoritmo, a partir do sinal sısmico gerado pela fonte, para cada ponto

do modelo (i, j), a matriz de tempo de transito (MTT ), registrara o tempo de chegada

da frente de onda de maxima amplitude naquele ponto. Juntamente com o calculo

da matriz MTT , tambem e feito o registro da maxima amplitude naquele ponto pela

matriz de amplitude, Preg, uma vez que esta faz parte da condicao necessaria para

20

registrar os tempos de transito (linha 1).

Algoritmo 01

1 se (|P(i, j, n)| > |Preg(i, j)|) entao

2 Preg(i, j)← |P(i, j, n)|3 MTT (i, j)← n ∗ ∆t

4 fim se

fim

onde:

P(i, j, n) e o campo de ondas incidente em cada ponto da malha no intervalo

de tempo n;

Preg(i, j) e a matriz que registra a amplitude em cada ponto.

O segundo criterio leva em consideracao o intervalo de tempo associado ao com-

primento de onda da fonte sısmica, dada pela expressao (2.12). Para detectar a

amplitude maxima nas proximidades da primeira quebra, segue o Algoritmo 02,

abaixo.

Algoritmo 02

1 cond1 ←[(t − MTT (i, j)) ≤ 1, 5 ∗ t f

]2 cond2 ←

[|P(i, j, n)| > |Preg(i, j)|

]3 cond3 ←

[Preg(i, j) = 0

]4 se (cond2 e (cond1 ou cond3)) entao

5 Preg(i, j)← |P(i, j, n)|6 MTT (i, j)← n ∗ ∆t

7 fim se

fim

onde cond1, cond2 e cond3 sao variaveis logicas e, em ambos os casos, t = n∆t.

A cada passo de tempo n, quando o campo em cada ponto (i, j) e recalculado,

os algoritmos analisam se o campo nesse novo passo satisfaz as condicoes impostas.

Caso sejam, ele grava o tempo t, substituindo o valor de tempo de transito para este

ponto ate encontrar um valor do campo de onda que volte a satisfazer as condicoes

substituindo pelo valor registrado anteriormente como sendo o campo de maxima

amplitude ou de amplitude maxima nas proximidades da primeira quebra da onda

direta.

21

2.2 Tempo de Transito Via Solucao da Equacao

Iconal

2.2.1 Equacao Iconal

Aplicando a transformada de Fourier em ambos os lados da equacao (2.1), obtem-

se a equacao da onda acustica 2D no domınio da frequencia (equacao de Helmholtz)

dada por (BLEISTEIN, 1984; CERVENY, 1987; POPOV, 2002):

∂2U(x, ω)∂x2 +

∂2U(x, ω)∂z2 −

ω2

C2(x)U(x, ω) = 0, (2.24)

onde U e o campo de pressao no domınio espaco-frequencia e ω representa a frequen-

cia angular.

A equacao da onda acustica no domınio da frequencia (2.24) admite uma solucao

baseada em uma serie assintotica (JEFFREYS e JEFFREYS, 1950), como segue:

U(x, ω) = eiωτ(x)∞∑

k=0

Ak(x)(−iω)k , (2.25)

onde Ak(x) e τ(x) sao funcoes que fornecem amplitude e o tempo de transito, respec-

tivamente. A serie descrita em (2.25) e conhecida como serie de raios (ray series).

Considerando a aproximacao de ordem zero da serie (2.25) (altas frequencias -

ω −→ ∞) e fazendo A0(x) = A(x), tem-se:

U(x, ω) = A(x)eiωτ(x). (2.26)

Substituindo (2.26) na equacao de Helmholtz (2.24), obtem-se a equacao iconal(∂τ(x)∂x

)2

+

(∂τ(x)∂z

)2

=1

C2(x)(2.27)

e a equacao de transporte

2(∂τ(x)∂x

∂A(x)∂x

+∂τ(x)∂z

∂A(x)∂z

)+ A(x)

(∂2τ(x)∂2x

+∂2τ(x)∂2z

)= 0. (2.28)

A partir da EDP, nao linear e nao-homogenea (2.27) e, geralmente, obtido um

sistema nao linear de equacoes diferenciais ordinarias (EDO) de primeira ordem

atraves do metodo das caracterısticas (ZAUDERER, 1989). Este sistema e conhecido

como sistema de tracado de raios ou equacoes cinematicas do raio.

22

2.2.2 Obtencao do Tempo de Transito

Diversas abordagens de solucao da equacao iconal (2.27) tem sido propostas na

literatura, isto se deve ao fato da mesma surgir em varias areas de aplicacao, como

computacao grafica (BRUSS, 1982), processamento de imagens (SETHIAN, 1999a),

geociencias (RAWLINSON e SAMBRIDGE, 2005), remocao de ruıdos (MALLADI

e SETHIAN, 1996), etc.

Como ha um grande numero de estrategias numericas diferentes para resolver

a equacao iconal, existem varias formas de classificar as tecnicas necessarias para

resolve-la. Mas, basicamente, pode-se identificar tres abordagens: a primeira e

aquela ja tradicionalmente conhecida que resolve a iconal utilizando o metodo das ca-

racterısticas (BLEISTEIN, 1984; CERVENY, 1987; ZAUDERER, 1989), a segunda

e aproximando a iconal via o metodo das diferencas finitas (SCHNEIDER et al.,

1992; VIDALE, 1988, 1990) e a terceira e aquela baseada na expansao de frente de

ondas usando esquemas de diferencas finitas avancados (upwind) (SETHIAN, 1996;

SETHIAN e POPOVICI, 1999). Neste ultimo caso esta enquadrado o metodo de-

nominado por Multistencil Fast Marching Method (MSFM) desenvolvido por HAS-

SOUNA e FARAG (2007) que e uma melhoria do Fast Marching Method (FMM)

proposto inicialmente por SETHIAN (1996, 1999b).

Como neste trabalho a obtencao do tempo de transito utilizando a iconal sera

feita atraves do MSFM e este e baseado no FMM, faz-se na proxima secao uma

descricao de ambas as metodologias.

2.2.2.1 Fast Marching Method

O FMM tem como principais vantagens sobre os demais metodos o fato de ser

altamente eficiente e incondicionalmente estavel, uma vez que segue o princıpio da

causalidade de maneira sequencial, isto e, a solucao e atualizada ponto a ponto em

ordem crescente, como explicado abaixo. Como este metodo ainda tem algumas

deficiencias, modificacoes objetivando melhorar sua precisao (RAWLINSON e LIN,

2003) e aumentar sua eficiencia (YATZIV et al., 2006) tem sido propostas.

Seguindo os trabalhos de HASSOUNA e FARAG (2007); SETHIAN (1996,

1999a) o lado esquerdo da equacao (2.27) e escrito como:

max(D−xi, j τ,−D+x

i, j τ, 0)2 + max(D−zi, jτ,−D+z

i, jτ, 0)2 =1

C2i, j

, (2.29)

onde D−i, j e D+i, j representam os operadores de diferencas finitas atrasadas, (A.6), e

progressivas, (A.7), de primeira ordem aplicados no correspondente ponto (i, j) da

malha. Como D−i, j e D+i, j sao aproximados, normalmente, utilizando operadores de

23

primeira ordem, a equacao (2.29) pode ser reescrita como:

max(τi, j − τ1

∆x, 0

)2+ max

(τi, j − τ2

∆z, 0

)2=

1C2

i, j

, (2.30)

sendo

τ1 = min(τi−1, j, τi+1, j

), (2.31)

τ2 = min(τi, j−1, τi, j+1

). (2.32)

A solucao da equacao (2.30) e dada como segue:

1. τi, j > max(τ1, τ2), entao τi, j e a solucao de maximo da equacao quadratica

(τi, j − τ1

∆x

)2+

(τi, j − τ2

∆z

)2=

1C2

i, j

;

2. τ2 > τi, j > τ1, entao τi, j = τ1 + ∆xCi, j

;

3. τ1 > τi, j > τ2, entao τi, j = τ2 + ∆zCi, j

.

A ideia do FMM e gerar frentes de ondas de tal forma que o princıpio de Huygens

seja satisfeito em cada ponto da malha, isto e, as frentes de onda se propaguem

baseadas na relacao de causalidade, o qual implica em dizer que o tempo de transito

τi, j em qualquer ponto depende apenas dos pontos vizinhos que tem o menor valor

de tempo registrado.

Para compreender como a frente de onda se propaga ou, para ser mais claro,

como os tempos de transitos sao obtidos, considere a figura 2.4(a), onde a frente

de onda sera propagada a partir da origem (valor conhecido). A origem esta sendo

representada por um ponto na cor preta, enquanto que as posicoes na cor clara

representam os pontos cujo valores do tempo de transito sao desconhecidos. Ao

resolver a equacao (2.30) para os pontos vizinhos a origem, como mostrado na figura

2.4(b), obtem-se uma possıvel solucao para os mesmos nesta regiao, sendo estes

representados por um quadrado na cor cinza. Estes pontos sao inseridos numa zona

fina, conhecida como narrow band. Agora, a frente avanca sobre um dos quadrados,

como segue: seleciona-se o ponto que possui o menor tempo de transito, τ, do narrow

band, suponha que seja o ponto A conforme figura 2.4(c). Daı, aplica-se novamente

a equacao (2.30) para encontrar os possıveis valores dos pontos vizinhos a A, figura

2.4(d). E assim, segue esse mesmo processo para os demais pontos, como ilustram

as figuras 2.4(e) e 2.4(f) para o caso do ponto B.

Como durante a propagacao da frente de onda sao obtidos tempos de transito

com valor temporario (pontos de cor cinza), estes sao armazenados em uma estrutura

24

(a) Busca dos pontos mais proximo a ori-gem.

(b) Calculo do tempo de transito dospontos vizinhos.

(c) Localiza o ponto A com menor tempode transito.

(d) A partir do ponto fixo A e calculadoo tempo de transito dos pontos vizinhos.

(e) Localiza o ponto B com menor tempode transito.

(f) A partir do ponto fixo B e calculadoo tempo dos vizinhos.

Figura 2.4: Ilustracao esquematica do funcionamento do FMM.

25

de dados do tipo arvore binaria, chamada heap. A cada passo do algoritmo um ponto

com menor valor de tempo e removido e outros sao adicionados.

O FMM honra o princıpio de Huygens, ja que o algoritmo so deve avancar par-

tindo do ponto de menor tempo ao longo da frente de onda, garantindo que cada

ponto da frente se comporte como uma nova fonte (fonte secundaria). Esta estra-

tegia assegura que no FMM as frentes de onda que se propagam no meio sejam

causais, estabelecendo um esquema estavel incondicionalmente (VAZ, 2004).

Uma outra abordagem do FMM e utilizar operadores de diferencas finitas de

segunda ordem para D−i, j e D+i, j, isto e, (A.8) e (A.9), respectivamente . Neste caso o

metodo e chamado Higher Accuracy Fast Marching Method (HAFMM) (SETHIAN,

1999a). Neste caso, a equacao (2.30) e reescrita como:

max[3(τi, j − τ1)

2∆x, 0

]2

+ max[3(τi, j − τ2)

2∆z, 0

]2

=1

C2i, j

, (2.33)

onde

τ1 = min(4τi−1, j − τi−2, j

3,

4τi+1, j − τi+2, j

3

), (2.34)

τ2 = min(4τi, j−1 − τi, j−2

3,

4τi, j+1 − τi, j+2

3

). (2.35)

Cabe destacar que, a primeira aproximacao do FMM, baseada nos esquemas de

diferencas finitas de primeira ordem pode resultar em grandes erros nos tempos de

transito, principalmente nos dois seguintes itens:

1. quando ha uma grande curvatura das frentes de onda, o que ocorre geralmente

em torno da fonte (ALKHALIFAH e FOMEL, 2001);

2. quando as frentes da onda se propagam em direcao diagonal em relacao a

orientacao da malha.

De modo geral, isto e melhorado com a utilizacao do HAFMM, o qual permite

uma reducao da ordem do erro no calculo das derivadas envolvidas, mas implicando

em um inevitavel aumento do custo computacional.

Para minimizar o erro no primeiro item, ALKHALIFAH e FOMEL (2001) pro-

puseram o FMM em coordenadas esfericas (3D) e polares (2D). Segundo os autores,

a implementacao em coordenadas polares melhorou a precisao da solucao, pratica-

mente, sem nenhuma perda da eficiencia computacional. Quanto ao segundo item,

HASSOUNA e FARAG (2007) propuseram o MSFM, que obtem a solucao pela ico-

nal em cada ponto da malha usando dois estencis e entao seleciona a solucao que

satisfaz a condicao de upwind, como explicado abaixo.

26

2.2.2.2 Multistencils Fast Marching Method 2D

De acordo com HASSOUNA e FARAG (2007), dois procedimentos sao possıveis

para minimizar o erro cometido na direcao diagonal. O primeiro consiste em usar

apenas um estencil centrado em x e alinhado com o sistema de coordenadas naturais.

Este sistema de coordenadas e, entao, rotacionado de tal forma a interceptar os

pontos das diagonais. Assim, a equacao iconal e resolvida ao longo do sistema xy

e ao longo do sistema xy. Logo, duas solucoes sao encontradas para o tempo de

transito em x, sendo que aquela que satisfaz a relacao de causalidade do FMM e

selecionada.

Esse procedimento esta limitado a malhas uniformes, onde o espacamento na

direcao x e igual aquele na direcao z (∆x = ∆z).

O segundo procedimento e mais geral, pois usa dois estencis centrados em x

e inclui todos os pontos vizinhos das direcoes horizontais e verticais mais aqueles

das diagonais. Dessa forma, os gradientes nas diagonais sao aproximados usando

derivadas direcionais.

Para descrever o procedimento, considere o estencil S 1, que cobre os pontos na

direcao do sistema de coordenadas naturais, e S 2, que cobre os pontos nas direcoes

diagonais, os quais interceptam os pontos da malha p1, p2, q1 e q2, conforme ilustrado

na figura 2.5. Seja os vetores unitarios −→r 1 = [r11, r12]T e −→r 2 = [r21, r22]T ao longo dos

segmentos de reta p2 p1 e q2q1, respectivamente, e τ1 e τ2 as derivadas direcionais ao

longo de −→r 1 e −→r 2, respectivamente, que podem ser escritas como:

(a) Estencil S 1. (b) Estencil S 2.

Figura 2.5: Ilustracao dos pontos envolvidos nos estencis S 1 e S 2 para um pontocentral x.

τ1 = −→r 1 · ∇τ(x) = r11τx + r12τz, (2.36)

τ2 = −→r 2 · ∇τ(x) = r21τx + r22τz, (2.37)

27

esta duas expressoes podem ser reescrita em forma matricial, entao: τ1

τ2

=

r11 r12

r21 r22

τx

τz

. (2.38)

Assim, tem-se que:

τ = R∇τ(x) → ∇τ(x) = R−1τ,

aplicando a operacao de transposicao em ambos lados na ultima expressao:

[∇τ(x)]T = (R−1τ)T → [∇τ(x)]T = τT R−T .

Uma vez que

|∇τ(x)|2 = [∇τ(x)]T∇τ(x),

entao

|∇τ(x)|2 = τT R−T R−1τ = τT (RRT )−1τ =1

C2(x). (2.39)

Se Θ e o angulo entre os vetores direcionais, entao RRT e dado como: ||−→r 1||2 −→r 1·

−→r 2−→r 1·−→r 2 ||

−→r 2||2

=

1 cos(Θ)cos(Θ) 1

(2.40)

e

(RRT )−1 =−1

sen2(Θ)

−1 cos(Θ)cos(Θ) −1

. (2.41)

Substituindo (2.41) em (2.39), chega-se a uma expressao fechada que relaciona

as derivadas direcionais do tempo de transito ao longo de um estencil arbitrario que

e obtida como segue:

τ21 − 2τ1τ2cos(Θ) + τ2

2 =sen2(Θ)

C2i, j

, (2.42)

A aproximacao de primeira ordem para a derivada direcional τv que obedece a

28

solucao da equacao (2.29) e expressa como:

τv = max(τi, j − τv

||xi, j − xv||, 0

), v = 1, 2, (2.43)

onde τi, j, τ1 e τ2 estao dados em (2.31) e (2.32) e xv e um ponto da malha conhecido

em que τv e mınimo.

Por outro lado, a aproximacao de segunda ordem que tambem obedece (2.29) e

expressa como sendo:

τv = max(

3(τi, j − τv)2‖xi, j − xv‖

, 0), v = 1, 2, (2.44)

onde τ1 e τ2 estao dados em (2.34) e (2.35) e assim como mencionado anteriormente,

xv e um ponto da malha conhecido em que τv e mınimo.

De acordo com HASSOUNA e FARAG (2007), para facilitar a compreensao,

toma-se uma malha uniformemente distribuıda, isto e, ∆x = ∆z = h que implica

Θ = π2 . Neste caso, a solucao na direcao do sistema de coordenadas naturais pode

ser obtida utilizando as aproximacoes de primeira ordem, dadas em (2.30), ou de

segunda, dadas em (2.33). Por outro lado, para as direcoes diagonais, tem-se as

expressoes (2.45) e (2.46) para os esquemas de primeira e segunda ordem, respecti-

vamente.

max(τi, j − τ1√

2h, 0

)2

+ max(τi, j − τ2√

2h, 0

)2

=1

C2i, j

, (2.45)

max[

3(τi, j − τ1)

2√

∆x2 + ∆z2, 0

]2

+ max[

3(τi, j − τ2)

2√

∆x2 + ∆z2, 0

]2

=1

C2i, j

. (2.46)

Para ambos os estencis, se τi, j > max (τ1, τ2), entao as expressoes (2.30), (2.33),

(2.45) e (2.46) podem ser reduzidas a uma equacao do segundo grau na forma:

g(h)[(τi, j − τ1

)2+

(τi, j − τ2

)2]

=1

C2i, j

, (2.47)

onde g esta dado na tabela 2.1.

Tabela 2.1: Valores da funcao g = g(h) dependente do tipo de orientacao e esquemanumerico utilizado.

Estencil Primeira ordem Segunda ordemS 1 g = 1/h2 g = 9/4h2

S 2 g = 1/2h2 g = 9/8h2

A condicao upwind para o caso mais geral (para detalhes sobre a deducao desta

29

condicao, ver HASSOUNA e FARAG (2007)), e dada por:

|τ1 − τ2| ≤g(∆x,∆z)sen(Θ)

Ci, j. (2.48)

Neste caso, o valor de g = g(∆x,∆z) para ambos os estencis e apresentado na

tabela 2.2.

Tabela 2.2: Coeficientes de condicao upwind para ambos os estencis S 1 e S 2.

Estencil Primeira ordem Segunda ordemS 1 g = min(∆x,∆z) g = 2 min(∆x,∆z)S 2 g =

√∆x2 + ∆z2 g = 2

√∆x2 + ∆z2

2.3 Suavizacao do Campo de Vagarosidade

E altamente desejavel que a matriz de tempo de transito seja o mais suave pos-

sıvel, isto e, apresente um reduzido numero de descontinuidades para melhorar o

calculo dos gradientes associados. Sabe-se que, estas descontinuidades sao causadas

durante a propagacao do campo de ondas pelas diversas reflexoes e reverberacoes

provenientes principalmente na presenca de altos contrastes de impedancia entre as

interfaces que podem interagir construtivamente nos valores das amplitudes (BUL-

CAO, 2004). Como observado nos algoritmos mostrados na secao 2.1.3, um efeito

construtivo nas amplitudes significa que a matriz de tempo, MTT , podera regis-

trar tempos que nao o da onda direta, isto e, tempos maiores do que deveria ser

registrado.

Dentro deste contexto, um procedimento que pode ser empregado para minimizar

as reflexoes decorrentes da propagacao do campo de ondas, visando melhorar os

resultados associados a obtencao da matriz de tempo de transito atraves da reducao

das descontinuidades presentes na mesma, e a suavizacao do modelo de velocidades,

para diminuir os contrastes de impedancia acustica ao longo do modelo.

A suavizacao e feita utilizando o campo de vagarosidade (inverso da velocidade),

conforme mostra a expressao (2.49), de modo a garantir que o tempo de transito

necessario para se alcancar um determinado ponto da malha nao seja alterado sig-

nificativamente (LOEWENTHAL et al., 1987).

MVi, j =1

2 ∗ nps + 1

nps∑k=−nps

1Ci+k, j

, (2.49)

onde MVi, j e a media movel atribuıda a cada no (i, j) da malha considerando nps

pontos em uma dada direcao i e/ou direcao j, sendo Ci, j a velocidade em cada no

30

(i, j).A fim de mostrar o efeito do processo de suavizacao sobre o modelo de velocida-

des, o conjunto de figuras 2.6 ilustram a reducao do contraste de impedancia entre

as camadas quando aplicado uma suavizacao com nps = 10 pontos tanto na direcao

i quanto na direcao j. Ressalta-se que, embora o processo de suavizacao seja crucial

para a reducao do contrate de impedancia, o valor de nps adotado nao deve ser alto

a tal ponto que influencie significativamente na direcao de propagacao da onda.

(a) Modelo sem suavizacao. (b) Perfil de velocidades do modelo ao lado.

(c) Modelo com suavizacao. (d) Perfil de velocidades do modelo ao lado.

Figura 2.6: Ilustracao do efeito produzido pelo processo de suavizacao com nps = 10pontos por direcao no macro-modelo de velocidades.

31

Capıtulo 3

Obtencao de Angulos de

Incidencia Utilizando MTT

Durante anos a tecnica de AVO tem sido utilizada com ou sem sucesso, normal-

mente como um caracterizador da litologia, mas em muitos casos como um indicador

da existencia de hidrocarbonetos (VASQUEZ, 1999). A caracterizacao da subsuper-

fıcie atraves dos estudos de AVO pode ser crucial na fase do processamento sısmico

para evidenciar melhor a existencia de hidrocarbonetos. Sendo que, o seu grande

exito e devido ao fato de ser aplicado em meios com baixa ou sem nenhuma varia-

cao lateral de impedancia onde usar AVO ou AVA nao faz muita diferenca (SILVA,

2009), daı o porque da amplitude refletida ser medida em funcao do afastamento

e nao do angulo de incidencia/reflexao. No entanto, a analise de AVO e originaria

dos estudos de AVA (YILMAZ, 2001), como pode ser observado nas equacoes de Zo-

eppritz (Zoeppritz (1919) apud CASTAGNA e BACKUS (1993)) onde a amplitude

refletida e funcao do angulo e dos parametros do meio (PAUL, 2004).

Cada vez mais os refletores de interesse estao localizados em altas profundidades

e imersos em meios com grande grau de complexidade, como aqueles envolvendo

flancos salinos, expressivas variacoes laterais de velocidade e mergulhos bastante

acentuados, requerendo que o estudo das curvas de refletividades sejam feitas atraves

da tecnica de AVA, uma vez que a caracterıstica mais importante do processo de

reflexao e sua dependencia com o angulo, ou seja, a quantidade de energia que e

refletida em uma interface depende do angulo de incidencia do campo de ondas com

relacao a normal do refletor (SILVA, 2009).

A partir deste contexto, desenvolve-se neste trabalho uma metodologia de extra-

cao de angulos de incidencia/reflexao que consiste na utilizacao da matriz de tempo

de transito (MTT) e do macro-modelo de velocidades. Basicamente, a ideia e fazer

a aplicacao do operador gradiente na MTT para extrair o angulo entre o campo

incidente e a componente vertical do mesmo e, de igual modo, repete-se o mesmo

procedimento para o macro-modelo de velocidades, onde e obtido o angulo formado

32

pela reta normal ao refletor e a componente vertical do gradiente. A partir daı,

faz-se a composicao destes dois angulos de tal modo a obter o angulo incidente em

cada ponto ao longo da superfıcie refletora.

Neste capıtulo, serao apresentados a metodologia desenvolvida para obter os

angulos de incidencia e os metodos utilizados para calcular as derivadas numerica-

mente. Mas, primeiramente, faz-se uma descricao sobre frente de onda e como ocorre

o fenomeno de particao de energia quando esta incide sobre um refletor contendo

diferentes propriedades acusticas. No final do capıtulo e mostrado uma outra abor-

dagem desenvolvida neste trabalho para obter os angulos, correspondentes aqueles

oriundos da matriz de tempo de transito, aplicando o operador de diferencas finitas

de quarta ordem durante a propagacao do campo de ondas.

3.1 Conceitos Preliminares

3.1.1 Frente de Onda

A frente de onda e dada pelas curvas de nıvel da variavel tempo, isto e, t =

τ(x) = constante. Matematicamente, as curvas de nıvel descrevem curvas no plano,

e superfıcies no espaco, ao longo das quais uma funcao assume valores constantes.

As curvas de nıvel sao bastante utilizadas em mapas meteorologicos para fazer a

representacao de pontos de mesma temperatura, chamadas isotermas ou da mesma

pressao, chamadas isobaricas. Outro tipo comum de curvas de nıvel sao aquelas

utilizadas com frequencia para representar mapas topograficos, onde as curvas de

nıvel representam colecoes de pontos de mesma altura (LARSON et al., 1998). Daı,

fica caracterizado que as frentes de onda descritas pelo tempo sao curvas de nıvel

no espaco bidimensional e superfıcies no tridimensional (POPOV, 2002).

(a) Isocronas correspondentes a propa-gacao ao lado.

(b) Propagacao de campo de ondasem diferentes instantes de tempo.

Figura 3.1: Frentes de onda definidas pelas isocronas associadas a amplitude maximacorrespondente a propagacao do campo de ondas.

Sendo assim, as curvas de nıvel que representam pontos que possuem o mesmo

33

tempo sao chamadas isocronas. Conforme ilustrado na figura 3.1, durante o processo

de propagacao, para o caso deste trabalho, as isocronas estarao sempre associadas

aos registros da amplitude maxima ou da amplitude maxima nas proximidades da

primeira quebra do sinal da fonte sısmica que passa por cada ponto da malha.

3.1.2 Operador Gradiente

O operador gradiente e o elemento chave da metodologia proposta neste trabalho,

pois e atraves da aplicacao deste operador sobre o macro modelo de velocidades e so-

bre a matriz de tempo de transito que serao obtidos os angulos auxiliares necessarios

para compor o angulo de incidencia/reflexao em profundidade.

A aplicacao do operador gradiente sobre uma funcao G(x, z) no ponto (x, z), de-

notado por−→∇G(x, z), e definido como:

−→∇G(x, z) = (Gx,Gz) =

(∂G(x, z)∂x

,∂G(x, z)∂z

), (3.1)

e sua magnitude∣∣∣−→∇G(x, z)

∣∣∣ dada por:

∣∣∣−→∇G(x, z)∣∣∣ =

√G2

x + G2z . (3.2)

Figura 3.2: Orientacao da forma como e obtido o angulo θ.

Utilizando a equacao (3.1), pode-se entao calcular o angulo associado a direcao

do gradiente com a componente horizontal ou vertical. Conforme ilustrado na figura

(3.2), o angulo sera obtido em relacao a componente vertical pela seguinte expressao:

θ(x, z) = tan−1(Gx

Gz

). (3.3)

34

Ressalta-se que o calculo numerico das derivadas na expressao (3.1) esta direta-

mente relacionado com a precisao do angulo obtido, sendo extremamente desejavel,

portanto, utilizar aproximacoes adequadas nestes calculos.

3.1.3 Reflexao e Refracao das Ondas Acusticas

Considere a figura 3.3 onde um raio de propagacao1 de uma onda compressional

incide obliquamente sobre um refletor plano horizontal. Pelo princıpio de Huygens,

havera o fenomeno de particao de energia (amplitudes) devido a diferenca de pro-

priedades acusticas (velocidade entre os meios - C1 e C2) levando ao surgimento de

um raio de onda refletido e um transmitido. No meio acustico os raios incidente, re-

fletido e transmitido sao sempre perpendiculares as suas respectivas frentes de onda

(CERVENY, 1987; POPOV, 2002).

Associado a cada raio de propagacao, ha naturalmente o surgimento dos angulos

definidos entre os raios de propagacao e a normal ao refletor, sendo que o angulo de

incidencia e igual ao angulo de reflexao.

Figura 3.3: Raio incidente gerando um raio refletido e transmitido no meio acustico.

Como sera mostrado na secao 4.1.1, os angulos γ e φ podem ser relacionados

atraves da Lei de Snell a qual estabelece que os senos dos angulos de incidencia e

transmissao sao diretamente proporcionais as velocidades da onda nos respectivos

meios, ou seja:

C2sen(γ) = C1sen(φ), (3.4)

ou, reescrevendo na forma tradicional:

sen(γ)C1

=sen(φ)

C2= %, (3.5)

1Linha imaginaria na direcao de propagacao da onda.

35

onde % e conhecido como parametro do raio.

A partir da Lei de Snell e possıvel descrever o surgimento das head waves que

ocorrem a partir do angulo crıtico quando toda a energia incidente passar a ser

refletida. As head waves surgem para equilibrar o campo de pressoes em angulos

maiores que o crıtico, onde as frentes de onda abaixo e acima do refletor se separam

causando uma descontinuidade, conforme ilustra figura 3.4. Isto ocorre so se C1 < C2

e, naturalmente, para φ = 90 ◦. Nesta situacao a equacao (3.5) e escrita na forma

da expressao (3.6), onde γc e o angulo crıtico.

C1

C2= sen(γc). (3.6)

Figura 3.4: Ilustracao do surgimento das head waves.

Uma forma de medir o quanto de energia e refletido e transmitido quando o raio

de propagacao incide sobre um certo ponto do refletor e atraves dos coeficientes de

reflexao e transmissao. O coeficiente de reflexao e dado por:

Rr(x, z, γ) =κcos(γ) −

[1 − κ2sen2(γ)

]1/2

κcos(γ) +[1 − κ2sen2(γ)

]1/2 , (3.7)

e o de transmissao

Rt(x, z, γ) = 1 + Rr(x, z, γ) =2κcos(γ)

κcos(γ) +[1 − κ2sen2(γ)

]1/2 , (3.8)

onde κ = C2C1

.

Na teoria dos raios, os coeficientes de reflexao e transmissao governam as ampli-

tudes do raio quando este se reflete ou se transmite atraves de uma interface suave.

36

Isto ocorre porque na vizinhanca da interface, esta se comporta como plano tangente

e a frente de onda do raio incidente e aproximada como uma frente de onda plana

(VASQUEZ, 1999).

3.2 Metodologia Proposta para Obtencao de An-

gulo de Incidencia

A metodologia proposta baseia-se em obter angulos auxiliares a partir da MTT e

do modelo de velocidades para compor o angulo final de incidencia em profundidade,

conforme explicado a seguir.

3.2.1 Angulo da MTT

Novamente, considere um vetor de propagacao ou raio de propagacao para uma

dada isocrona. Uma vez que a direcao de propagacao e perpendicular as isocronas,

o vetor gradiente tambem cumpre o mesmo papel, ou seja, fornece a direcao de

propagacao. Sendo assim, pode-se obter o angulo entre a direcao do raio de pro-

pagacao e a componente vertical do vetor gradiente (eixo z do modelo) usando as

componentes do gradiente obtido a partir das isocronas, conforme figura a 3.5. Este

e o primeiro angulo auxiliar necessario para obter o angulo de incidencia/reflexao, o

qual e chamado aqui de angulo α.

Figura 3.5: Raio de propagacao dado pela direcao do vetor gradiente e isocronasrepresentando as frentes de onda.

3.2.2 Angulo do Modelo de Velocidades

Em se tratando do modelo, considere um refletor plano em situacao inclinada,

conforme figura 3.6. Da mesma forma que descrito anteriormente, pode-se obter o

37

angulo entre a direcao do vetor gradiente e a vertical quando o operador gradiente e

aplicado sobre o macro-modelo de velocidades, ja que este vetor e paralelo a direcao

normal do refletor. Portanto, a aplicacao do operador gradiente fornece o segundo

angulo auxiliar dado como sendo entre a reta normal ao refletor e a componente

vertical do vetor gradiente, o qual e chamado de angulo β.

Figura 3.6: Vetor gradiente e suas componentes em um certo ponto do refletor.

3.2.3 Composicao do Angulo de Incidencia/Reflexao

Para fazer a composicao do angulo de incidencia, vai-se sintetizar as duas con-

sideracoes feitas nas duas secoes anteriores. Para tanto, considere a figura 3.7 que

ilustra os angulos α, sendo formado entre o raio incidente e a reta vertical, e β, o an-

gulo entre a reta normal ao refletor e a reta vertical. Ambos os angulos sao oriundos

da aplicacao das equacoes (3.1), (3.2) e (3.3), sendo que, como ja visto, o primeiro

vem da matriz de tempo de transito e o segundo do macro-modelo de velocidades.

Como resultado final, faz-se a aplicacao da equacao (3.9) para obter o angulo γ de

incidencia/reflexao em cada ponto da malha.

γ =∣∣∣α − β∣∣∣ (3.9)

Os demais casos em que sao consideradas as diversas formas que um raio pode

incidir sobre o refletor, assim como as posicoes que este se encontra em subsuperfı-

cie (inclinado, vertical e horizontal) nao foram mostrados por se tratar de simples

situacoes semelhantes ao caso mostrado aqui.

38

(a) Raio de propagacao incidindo pelo ladodireito em relacao a normal.

(b) Raio de propagacao incidindo pelo ladoesquerdo em relacao a normal.

Figura 3.7: Exemplo de um raio de propagacao incidindo sobre um refletor inclinadoonde as componentes do vetor gradiente da matriz de tempo de transito sao positivos.

3.3 Metodos Numericos para Calculo de Gradi-

entes

A principal dificuldade da aplicacao da metodologia desenvolvida e, numerica-

mente, a obtencao das derivadas que compoem o vetor gradiente. Isto porque tanto

a matriz de tempo de transito quanto o macro-modelo de velocidades possuırem

fortes descontinuidades e a operacao de derivacao leva em conta a influencia local

das propriedades sendo bastante sensıvel a pequenas alteracoes.

Para diminuir as fortes descontinuidades e aplicar os procedimentos de obtencao

de derivadas explicados nas secoes seguintes, faz-se uso do procedimento de suavi-

zacao mostrado na secao 2.3, pois este reduz as reflexoes advindas da propagacao

do campo de ondas durante a obtencao da matriz de tempo de transito. Sendo as-

sim, para a matriz de tempo, menos descontinuidades significa menos contraste de

impedancia acustica por parte do modelo que a gerou, assim como reducao da forte

descontinuidade deste, como mostrado na figura 2.6.

Serao abordados, tres procedimentos de obtencao das derivadas, a saber: os

tradicionais operadores de diferencas finitas (FORTUNA, 2000; STRIKWERDA,

2004); as funcoes de base radial (RBFs - Radial Basis Functions) (BUHMANN, 2003;

LIU, 2002) e splines cubicas locais conhecidas como polinomio de Akima (AKIMA,

1974a,b, 1996), os quais serao mostrados nas proximas secoes.

Uma outra abordagem desenvolvida neste trabalho foi a extracao do gradiente,

atraves do operador de diferencas finitas de quarta ordem, durante a propagacao do

campo de ondas, isto e, diretamente do campo de pressao, como sera visto na secao

3.3.4.

39

3.3.1 Gradiente Via Operadores de Diferencas Finitas

No apendice A, encontra-se o desenvolvimento de como sao obtidos os principais

operadores de diferencas finitas utilizados ao longo deste trabalho, assim como o

operador de derivada primeira de quarta ordem utilizado na obtencao das derivadas

numericas.

3.3.2 Gradiente Via Funcoes de Base Radial

As funcoes de base radial sao muito conhecidas no contexto das redes neurais

artificiais (RNAs) (BISHOP, 1996) e mais recentemente vem sendo usadas nas apli-

cacoes que necessitam resolver EDPs, impondo-se como uma boa alternativa ao

metodo dos elementos finitos (MEF), por exemplo, ja que nao ha a necessidade do

processo de geracao de malha (LIU, 2002).

Embora o metodo das RBFs para aproximacao de funcoes seja hoje uma das

tecnicas mais usadas para aproximar dados dispersos (BUHMANN, 2003), o seu

uso neste trabalho se deu devido a facilidade de implementacao e facil controle dos

pontos de influencia que corroboram para encontrar o valor de um certo ponto em

questao, como sera visto a seguir.

Para facilitar a exposicao da obtencao das derivadas utilizando as RBFs, con-

sidere um conjunto de dados arbitrarios e distintos a serem interpolados, isto e,

{(x1, u1), (x2, u2),...,(xi, ui),...,(xN , uN)}, onde ui e o valor da funcao u amostrada no

ponto xi ∈ R2 e N corresponde ao numero total de pontos da malha. Deseja-se

encontrar uma funcao de interpolacao ui = u(xi) = u(xi) da forma:

ui =

N∑j=1

σ jψi j (3.10)

e, posteriormente, suas derivadas parciais:

∂ui

∂x=

N∑j=1

σ j∂ψi j

∂x, (3.11)

∂ui

∂z=

N∑j=1

σ j∂ψi j

∂z, (3.12)

sendo σ j os pesos da funcao de interpolacao a serem determinados, ψi j = ψ(ri j) uma

funcao de base radial, ou seja, uma funcao radial com respeito a distancia euclidiana

em torno do ponto xi (centro), onde ri j e dado por:

ri j =

√(xi − x j)2 + (zi − z j)2. (3.13)

40

Pertencente a uma classe especial de funcoes, a RBF ψ = ψ(r) esta associada a

uma funcao unidimensional cujos valores, em geral, crescem ou decrescem monoto-

nicamente em relacao a distancia r.

Para possibilitar a determinacao dos pesos σ j em (3.10), MICCHELLI (1986)

mostrou a existencia de uma gama de funcoes adequadas para interpolacao que ge-

ram um sistema linear de equacoes consistente e determinado. A tabela 3.1 mostra

as tradicionais funcoes utilizadas como RBFs (BUHMANN, 2003; LIU, 2002; SCHA-

BACK, 2000). Algumas delas possuem parametros livres chamados de parametros

de forma os quais controlam a forma da RBF, como pode ser visto na figura 3.8.

Estes parametros, portanto, precisam ser ajustados para se conseguir uma melhor

performace (LIU, 2002), como e o caso do parametro ζ que controla o raio de in-

fluencia de cada funcao e que pode alterar significativamente a qualidade da solucao.

Tabela 3.1: RBFs tradicionais.Nome RBF Parametro de forma

Multiquadrica ψ(r) = (r2 + ζ)q ζ, qGaussiana ψ(r) = e−ζr ζ

Logarıtmica ψ(r) = rζlog r ζ

Quando q = ±1/2 a funcao multiquadrica e classicamente conhecida como mul-

tiquadrica de HARDY (1971). Destaca-se que, neste trabalho foi utilizada a funcao

multiquadrica tomando o caso particular quando q = −1/2, a qual e chamada de

multiquadrica inversa. Dessa forma, para uma malha com dimensao total de N pon-

tos, e gerada uma matriz de interpolacao densa da ordem de N2 que apresenta um

custo computacional proibitivo para os problemas a serem analisados. Entao, para

contornar este problema, ao inves de fazer uma interpolacao global com um custo

computacional elevado, adotou-se por usar sub-regioes com tamanho N = 19x19pontos e obter as derivadas parciais apenas no ponto central destas.

Embora existam alguns procedimentos sugeridos na escolha do parametro de

forma ζ a ser usado na funcao multiquadrica (FASSHAUER, 2002; FRANKE, 1982;

HARDY, 1971; KANSA, 1990; RIPPA, 1999), adotou-se um esquema semelhante ao

leave-one-out cross validation cuja ideia e selecionar o parametro atraves da mini-

mizacao de uma funcao custo que coleciona erros para uma sequencia de previsoes

parciais dos dados (LIMA, 2009). Da serie de testes realizados, observou-se que

ζ = 10−4 gera os melhores resultados.

3.3.3 Gradiente Via Spline Cubica Local

A fim de evitar resolver um sistema linear de ordem elevada, diminuir a com-

plexidade e o tempo computacional para fazer a interpolacao de um conjunto de

41

(a) Multiquadrica (q = 1/2). (b) Multiquadrica inversa (q = −1/2).

(c) Gaussiana. (d) Logarıtmica.

Figura 3.8: RBFs tradicionais.

dados, AKIMA (1970) apresentou um novo conceito de interpolacao, inicialmente,

unidimensional, mais tarde estendido para o caso bidimensional aplicado em malhas

retangulares (AKIMA, 1974a,b) e triangulares (AKIMA, 1978a,b).

O grande diferencial do metodo de interpolacao desenvolvido pelo autor e que,

ao inves de ser aplicado no conjunto de dados como um todo, e aplicado localmente

utilizando um polinomio de terceiro grau. Essa caracterıstica o torna muito diferente

da spline cubica tradicional a qual e aplicada em todo conjunto de dados gerando

um sistema linear de ordem elevada. Por ser baseado em um polinomio cubico, a

interpolacao usando o polinomio de Akima costuma ser chamada de spline de Akima

constituindo-se numa classe de spline.

Com o proposito de fazer a explicacao do procedimento de extracao das derivadas

pelo polinomio de Akima, considere o mesmo conjunto de dados definidos na secao

3.3.2 a serem interpolados. O metodo e baseado em uma funcao contınua por partes

composta por um conjunto de polinomios bicubicos ui j da forma (AKIMA, 1974a,b):

ui j =

3∑a=0

3∑b=0

wab(x − xi)a(z − z j)b, (3.14)

42

onde wab sao os pesos do polinomio.

Cada polinomio bicubico e aplicado em uma celula retangular, (x, z) ∈

[xi, xi+1]x[z j, z j+1], sendo determinado atraves dos valores dos dados, ui, j, e pelos

valores das derivadas parciais (3.15), (3.16) e (3.17) nos quatro cantos da celula, isto

e, nos pontos (i, j), (i, j + 1), (i + 1, j + 1) e (i + 1, j) mostrados na figura 3.9.

∂ui j

∂x=

3∑a=0

3∑b=0

awab(x − xi)a−1(z − z j)b, (3.15)

∂ui j

∂z=

3∑a=0

3∑b=0

bwab(x − xi)a(z − z j)b−1, (3.16)

∂2ui j

∂x∂z=

3∑a=0

3∑b=0

abwab(x − xi)a−1(z − z j)b−1. (3.17)

Figura 3.9: Malha e celula computacional mostrando os pontos que determinam opolinomio bicubico.

Como evidenciado, a interpolacao pelo metodo Akima requer que sejam determi-

nadas as derivadas parciais para, posteriormente, poder resolver um sistema linear

localmente. A forma como sao calculadas as derivadas parciais e que torna este

metodo unico. Como neste trabalho se utiliza apenas as derivadas parciais (3.15)

e (3.16), faz-se abaixo a descricao de como as mesmas sao calculadas. Esta descri-

cao sera feita seguindo os trabalhos de AKIMA (1991) para o caso bidimensional

e AKIMA (1996) para o caso unidimensional, nos quais Akima apresenta um novo

procedimento para obter as derivadas.

Como se trata apenas das derivadas nas direcoes do sistema de coordenadas

cartesianas, toma-se, como exemplo, para descrever o procedimento apenas a direcao

do eixo x, como consta em AKIMA (1991). Sendo assim, a equacao (3.14) e reescrita

como:

u =

3∑a=0

wa(x − xi)a. (3.18)

43

Usando os valores dos dados correspondentes ao pontos i e i + 1, Akima deter-

minou que os coeficientes do polinomio cubico sao dados por:

w0 = ui

w1 = u′

i,

w2 =3mi − 2u

′

i − u′

i+1

xi+1 − xi,

w3 =u′

i + u′

i+1 − 2mi

xi+1 − xi,

onde mi e a inclinacao do segmento de reta que liga os pontos i e i + 1, ou seja:

mi =ui+1 − ui

xi+1 − xi.

A derivada, u′

i, e encontrada atraves da expressao:

u′

i =

∑4k=1 ρky

′

k∑4k=1 ρk

. (3.19)

onde y′

k e a derivada calculada em i por um polinomio de terceiro grau passando

pelo ponto i e outros tres pontos vizinhos (figura 3.9):

y′

1 = F(i, i − 3, i − 2, i − 1), (3.20)

y′

2 = F(i, i − 2, i − 1, i + 1), (3.21)

y′

3 = F(i, i − 1, i + 1, i + 2), (3.22)

y′

4 = F(i, i + 1, i + 2, i + 3), (3.23)

sendo F um polinomio de grau tres aplicado nos respectivos pontos acima utilizado

na determinacao de y′

k.

Para os pontos i, i + 1, i + 2 e i + 3, por exemplo, o polinomio e dado por:

F(i, i + 1, i + 2, i + 3) =`1 + `2 + `3

4, (3.24)

onde

44

`1 = (ui+1 − ui)(xi+2 − xi)2(xi+3 − xi)2xi+3xi+2,

`2 = (ui+2 − ui)(xi+3 − xi)2(xi+1 − xi)2(xi+1 − xi+3),

`3 = (ui+3 − ui)(xi+1 − xi)2(xi+2 − xi)2(xi+2 − xi+1),

4 = (xi+1 − xi)(xi+2 − xi)(xi+3 − xi)(xi+2 − xi+1)(xi+3 − xi+2)(xi+3 − xi+1).

Os pesos, ρk, sao inversamente proporcionais ao produto do que Akima chamou

de fator de volatilidade (volatility) e fator distancia:

ρk =1

εk$k. (3.25)

O fator distancia, $k, e a soma dos quadrados da distancia do ponto i aos outros

tres pontos vizinhos e o fator de volatilidade, εk, e a soma dos quadrados dos desvios

obtidos via regressao linear nesses mesmos pontos (o mesmo conjunto de pontos

mostrados nas equacoes (3.20) a (3.23)).

Neste trabalho, utilizou-se a sub-rotina RGPD3P (AKIMA, 1996) adaptada para

obter apenas as duas derivadas parciais, sendo que a mesma se encontra disponıvel

no repositorio Netlib2.

3.3.4 Gradiente Via Propagacao do Campo de Ondas

Neste trabalho, desenvolveu-se um procedimento alternativo para obtencao do

gradiente e, consequentemente, dos angulos auxiliares correspondentes aqueles oriun-

dos da matriz de tempo de transito. A ideia deste esquema e utilizar os mesmos

algoritmos usados para obter os tempos de transito (ver secao 2.1.3), isto e, usar as

condicoes que possibilitam registrar os tempos como condicao auxiliar para obter as

derivadas parciais do campo de pressao durante a propagacao do campo de ondas

em cada ponto da malha.

Partindo do Algoritmo 01, por exemplo, tem-se o Algoritmo 03, onde observa-se

que o procedimento requer que seja determinado o ponto de maxima amplitude,

Preg(i, j). Tao logo registrada tal amplitude em um certo ponto (i, j), da-se inıcio ao

processo do calculo das derivadas do campo de pressao, DPXatual e DPZatual, aplicado

sobre a parte restante da onda que passara naquele ponto. Isto e possibilitado gracas

a variavel T janela(i, j) que, ao inves de armazenar o tempo da amplitude maxima,

armazena o tempo que falta para todo restante do campo passar atraves do ponto

em questao, ou seja, a partir do campo de pressao maximo registrado, abre-se uma

janela de tempo para cada ponto da malha.

2http://netlib.sandia.gov/toms/index.html

45

Algoritmo 03

1 se (|P(i, j, n)| > |Preg(i, j)|) entao

2 Preg(i, j)← |P(i, j, n)|3 T janela(i, j)← n ∗ ∆t + t f

4 fim se

5 se (T janela(i, j) > n ∗ ∆t) entao

6 DPXatual ←∂P(i, j)∂x

7 DPZatual ←∂P(i, j)∂z

8 GRADatual ← (DPXatual)2 + (DPZatual)2

9 GRADreg ← (DPXreg(i, j))2 + (DPZreg(i, j))2

10 se (GRADatual > GRADreg) entao

11 DPXreg ← DPXatual

12 DPZreg ← DPZatual

13 fim se

14 fim se

fim

onde:

T janela(i, j) e uma matriz que contem os valores da janela de tempo criada

para cada ponto a partir do instante de tempo em que e registrada a amplitude

maior do que a anterior ate alcancar a maxima armazenada, Preg(i, j);

DPXatual e DPZatual sao os valores das derivadas parciais nas direcoes x e z do

campo de pressao no instante de tempo atual obtidos atraves dos operado-

res de diferencas finitas de quarta ordem, (A.13), assim como DPXreg(i, j) e

DPZreg(i, j) que contem as derivadas ja registradas;

GRADatual e GRADreg armazenam a soma dos quadrados das derivadas parciais,

respectivamente, no instante de tempo atual e naquele ja registrado.

Dentro da janela de tempo criada, varios valores para as derivadas parciais sao

calculados, selecionando-se as derivadas cuja soma dos quadrados seja maior do que

aquela ja registrada anteriormente (linha 10). Apos o campo de ondas percorrer todo

o modelo de velocidades, calcula-se os angulos auxiliares correspondentes aqueles

obtidos pela matriz de tempo de transito, conforme descrito na secao 3.2.

De forma opcional, no procedimento descrito anteriormente, tambem pode-se

obter a matriz de tempo, MTT . Nos exemplos mostrados no capıtulo 5 isso sera

realizado para efeito de comparacao.

Nota-se que obter o gradiente do campo de pressao no decorrer da propagacao do

campo de ondas com a metodologia descrita implica em um aumento consideravel

46

de tempo de processamento, assim como em um aumento da quantidade de memoria

utilizada. Estes sao fatores que, computacionalmente, podem tornar este procedi-

mento menos atraente, haja vista que a dimensao dos macro-modelos de velocidades

utilizados nas modelagens sısmicas terem dimensoes, geralmente, elevadas.

47

Capıtulo 4

Solucoes de Referencia

Com o intuito de comparar os resultados obtidos com a metodologia proposta

neste trabalho, implementou-se o tracamento cinematico de raios para a geracao dos

angulos de incidencia e das matrizes de tempo de transito de forma semi-analıtica

baseado na ideia do tracado de raios entre dois pontos (two-point ray tracing), mas

resolvido pelo shooting method (tracado de raios tradicional - ray tracing). Este tipo

de problema e considerado como tracado de raios de valor de contorno (boundary-

value ray tracing), cuja solucao e obtida pelo bending method, mas, devido as di-

ficuldades de implementacao e custo computacional, e amplamente resolvido pelo

shooting method (CERVENY, 1987).

Neste capıtulo, faz-se uma breve descricao tanto do shooting method quanto do

bending method para dar uma ideia destas duas abordagens, sendo que mais detalhes

destes dois metodos podem ser encontrados em CERVENY (1987) ou CERVENY

(2001). Com o objetivo de descrever as duas abordagens, apresenta-se inicialmente

o princıpio de Fermat e a Lei de Snell, necessarios para compreender os demais

conceitos desenvolvidos no capıtulo.

4.1 Conceitos Preliminares

4.1.1 Princıpio de Fermat e Lei de Snell

Para descrever inicialmente a teoria de raios, um princıpio bastante utilizado e

aquele formulado por Pierre Fermat conhecido como princıpio de Fermat (tempo

mınimo) (CERVENY, 1987), cuja ideia foi desenvolvida inicialmente para o caso

de raios opticos, mas que e igualmente aplicada para o caso dos raios de ondas

sısmicas. Este princıpio e baseado no fato de que para um conjunto de caminhos

percorridos pelos raios emitidos de um ponto a outro, aquele de caminho mais curto

nao corresponde a chegada de menor tempo entre os dois pontos.

Antigamente, pensava-se que o caminho percorrido pela luz indo de um ponto a

48

outro atraves de uma superfıcie refletora era o de caminho mais curto possıvel. Mais

tarde, descobriu-se que um feixe de luz atravessando uma interface nao toma uma

linha reta ou um caminho mais curto entre um ponto no meio incidente e um no

meio de transmissao. Fermat, propos seu princıpio de tempo mınimo valido tanto

para reflexao quanto refracao onde, basicamente, o caminho real entre dois pontos

tomados por um feixe de luz e aquele que “corre” em menor tempo.

Uma consequencia direta do princıpio de Fermat e a Lei de Snell, ja que o tracado

de raios corresponde aos tempos de transito estacionarios. Neste caso, entende-se a

estacionalidade como mınimo tempo, embora existam importantes situacoes fısicas

para as quais o tempo de transito e maximo (MARGRAVE, 2003). Uma vez que

se esta tratando apenas o caso acustico, a Lei de Snell e uma relacao que rege os

angulos de incidencia e transmissao da frente de onda compressional ao atravessar

uma interface separada por dois meios com propriedades acusticas distintas.

Como exemplo da aplicacao do princıpio de Fermat para o caso da refracao,

considere a figura 4.1(a) onde se deseja encontrar qual curva fornece o tempo de

transito para um campo de onda elementar (raio) que vai do ponto A ao ponto B.

Em outras palavras, deseja-se encontrar o tempo mınimo t com relacao a variavel L.

Sendo assim (HECHT, 2002):

t =AOC1

+OBC2

, (4.1)

isto e,

t =(H2 + L2)1/2

C1+

[H

2+

(L − L

)2]1/2

C2. (4.2)

Minimizando a equacao (4.2) ( dtdL = 0) e observando na figura 4.1(a) que

sen(γ) =L

(H2 + L2)1/2 ,

sen(φ) =L − L[

H2

+(L − L

)2]1/2 ,

tem-se:

sen(γ)C1

=sen(φ)

C2, (4.3)

a qual corresponde a Lei de Snell.

Desta forma, a Lei de Snell relaciona os angulos de incidencia, γ, e transmissao,

φ, da frente de onda de modo que os senos dos angulos de incidencia e transmissao

49

(a) Princıpio de Fermat aplicado a reflacao.

(b) Raio se propagando atraves de um meio estratificado.

Figura 4.1: Ilustracao geometrica do princıpio de Fermat.

sao diretamente proporcionais as velocidades das ondas nos respectivos meios (C1 e

C2). Portanto se um feixe de luz ou um campo de onda viaja do ponto A ao ponto

B no menor tempo possıvel, deve satisfazer a Lei de Snell.

No caso em que se tem um material composto de k camadas estratificadas, cada

uma com diferentes velocidades, como na figura 4.1(b), o tempo de transito de A

ate B sera dado por:

t =

k∑k=1

dk

Ck, (4.4)

onde dk e Ck denotam a extensao e a velocidade, respectivamente, associadas com a

k-esima camada.

50

Claramente que, para uma distribuicao contınua da velocidade, C(x), o tempo

de transito de um raio sısmico ao longo de um possıvel caminho que liga o ponto A

ao ponto B e dado pelo funcional de Fermat:

τ(x) =

∫ B

Adt =

∫ B

A

1C(x)

dl, (4.5)

onde dl denota o comprimento infinitesimal da curva ao longo do caminho.

Um exemplo de aplicacao, e o caso da tomografia sısmica interpocos onde, apro-

veitando a existencia de pocos ja perfurados, coloca-se a fonte em um poco e os

receptores em outro. Emitindo um campo de ondas sısmicas a partir da fonte,

registra-se os tempos de chegada das mesmas.

4.1.2 Abordagens Tradicionais do Tracado de Raios

Na teoria de raios, duas abordagens principais sao tradicionalmente aplicadas

(PIMENTEL, 1994; ZHANG et al., 2009): a primeira e mais usada e conhecida como

shooting method e a segunda como bending method. No caso da primeira abordagem,

a direcao e o ponto de partida do raio sao conhecidos e, assim, constituem um sistema

completo de condicoes iniciais para a determinacao da trajetoria do raio (initial-value

ray tracing). Ja no segundo caso, que tambem e conhecido como tracado de raios

de dois pontos (two-point ray tracing ou boundary-value ray tracing), a direcao do

raio nao e conhecida, exceto o ponto de partida (fonte) e o de chegada (receptor).

Uma forma de encontrar a trajetoria que satisfaca o sistema de tracado de raios e,

neste caso, determinar uma trajetoria inicial a qual e submetida a um processo de

pertubacao ou minimizacao ate encontrar o melhor caminho.

4.1.2.1 Shooting Method

Classicamente o shooting method padrao e conhecido por metodo de tracado de

raios (ray tracing) sendo baseado nas equacoes (2.27) e (2.28) advindas da ART.

Como ja dito anteriormente, a parte cinematica e responsavel por gerar os tempos

de transito, a frente de onda e as trajetorias do raio enquanto que a parte dinamica

e responsavel pelas amplitudes, a qual caracteriza a energia da onda. No entanto,

e possıvel tratar a parte cinematica sem fazer uso da teoria assintotica. De forma

muito mais simples, pode-se utilizar o princıpio de Fermat e a Lei de Snell, os quais

foram usados ao longo de muito tempo na sismologia para a geracao de resultados

bem apreciaveis (CERVENY, 2001). Para um conjunto de angulos de partida, o tra-

cado de raios gera um conjunto de caminhos onde o tempo de transito e os angulos

de incidencia sao obtidos ao longo das trajetorias (figura 4.2). Como dificuldades,

este metodo apresenta problemas quando aplicado em modelos contendo geometrias

51

complexas nas quais podem existir multiplos caminhos (multipath) e a geracao de

zonas de sombra (shadow zones). Como vantagem a relativa facilidade de imple-

mentacao e baixo custo computacional comparativamente a propagacao do campo

de ondas da equacao completa da onda (via MDF, por exemplo).

Figura 4.2: Ilustracao do shooting method.

4.1.2.2 Bending Method

Em relacao ao bending method, de acordo com CERVENY (2001), a ideia e tracar

um raio inicial que posteriormente deve ser perturbado iterativamente ate alcancar

uma solucao que satisfaca o sistema de equacoes do raio ou que minimize o tempo de

transito (figura 4.3). O raio inicial e apenas uma referencia auxiliar para ligar os dois

pontos e, alem disso, nao precisa satisfazer o princıpio de Fermat e nem mesmo a Lei

de Snell. Se o raio inicial arbitrado estiver longe do raio que efetivamente deveria

ser a solucao, o metodo diverge. Esta trajetoria inicial e estimada atraves de um

algoritmo independente chamado de estimador de raio (ray estimator). O bending

method pode ser aplicado como um pos-processador (postprocessor), corrigindo a

trajetoria inicial.

Por um lado, o bending method e geralmente mais rapido do que o shooting

method (CERVENY, 1987), mas tendo como desvantagem o fato de que os tempos

de transito obtidos nao podem ser garantidos como mınimos globais (ZHANG et al.,

2009). Alem disso, o bending method exige estimadores de raios adequados, caso

contrario, os raios e os tempos de transito podem resultar imprecisos ou mesmo

incorretos.

52

Figura 4.3: Ilustracao do bending method.

4.2 Geracao de Solucoes de Referencia

Para gerar as solucoes de referencia, implementou-se um procedimento iterativo

convergente para a obtencao dos angulos de incidencia e das matrizes de tempo de

transito, bem como das respectivas derivadas parciais, baseado no tracado de raios

de dois pontos com procedimento de solucao idealizado no shooting method. Assim,

sao conhecidos o ponto e a direcao inicial de partida do raio, alem disso, tambem

e conhecido o ponto em profundidade que o raio deve atingir. O esquema foi im-

plementado para macro-modelos de velocidades contendo varias camadas paralelas

horizontais e inclinadas.

Para descrever o procedimento, considere as figuras 4.4(a) e 4.4(b), as quais

apresentam um modelo homogeneo e outro constituıdo por tres camadas paralelas

horizontais homogeneas, respectivamente. A obtencao dos angulos, dos tempos de

transito e dos gradientes para o primeiro modelo e feita de forma imediata. Como

pode ser visto neste modelo, os dados analıticos para um certo ponto (xi, zi) sao

dados por:

γ = cos −1(|zi − z f |

di

), (4.6)

τi =di

Ci, (4.7)

∂τi

∂x=

xi − x f

di, (4.8)

53

∂τi

∂z=

zi − z f

di, (4.9)

onde γ e o angulo em relacao a vertical (angulo incidente quando o ponto (xi, zi)esta sobre a interface), τi e o tempo de transito, ∂τi

∂x e ∂τi∂z sao as derivadas parciais

que compoem o vetor gradiente ∇τ, e di e a distancia entre o ponto fonte (x f , z f ) e o

ponto em questao.

Para o caso do segundo modelo onde ha mais camadas e se deseja que um certo

raio atinja um determinado ponto da malha, o procedimento e iterativo. Para fazer

a descricao, suponha que o objetivo e tracar um raio que va da fonte ate o ponto

fixo (x3, z3) localizado na terceira camada. Basicamente, a ideia e:

1. lancar um raio inicial com um certo angulo em relacao a vertical (angulo

incidente), γ1;

2. identificar a posicao que este raio atinge a primeira interface, (x1, z1), atraves

da expressao (4.10), ja que, a priori, se sabe a profundidade da interface,

x1 = x f − |z1 − z f | tan(γ1); (4.10)

3. utilizar a Lei de Snell para conhecer o angulo de partida, φ1, que neste caso

corresponde ao angulo de incidencia na segunda interface, γ2, e, repetir o passo

2, identificando a posicao da incidencia via expressao

x2 = x1 − |z2 − z1| tan(γ2); (4.11)

4. novamente, e aplicado a Lei de Snell e verificado se o raio atingiu o ponto

desejado, (x3, z3), via expressao

x3 = x2 − |z3 − z2| tan(γ3). (4.12)

As expressoes (4.7) a (4.9) sao reescritas para o caso tratado, como segue:

τ3 =d1

C1+

d2

C2+

d3

C3, (4.13)

∂τ3

∂x=

x3 − x2

d3, (4.14)

∂τ3

∂z=

z3 − z2

d3. (4.15)

54

(a) Modelo homogeneo.

(b) Modelo de camadas paralelas horizontais.

(c) Modelo de camadas inclinadas.

Figura 4.4: Ilustracao de modelos discretizados para obtencao de dados.

55

Como, em geral, o raio nao atinge o ponto desejado na primeira tentativa, o

processo precisa ser repetido via incremento ou decremento do angulo de partida

na primeira camada, ∆γ. Alem do fato de o raio nao atingir o ponto em questao,

isto e, nao satisfazer o criterio de tolerancia estabelecido, o processo tambem pode

ser repetido caso o angulo de incidencia em uma dada interface ultrapasse o angulo

crıtico.

Alem do acrescimo ou decrescimo do angulo de partida, a garantia de convergen-

cia e feita atraves do tamanho do incremento, ∆γ, que, dependendo da necessidade,

assume valores cada vez menores para que seja possıvel alcancar o ponto desejado.

Para efeito de reducao do custo computacional, o procedimento e efeito apenas do

lado esquerdo ou direito em relacao a posicao da fonte, dependendo o quao proximo

esta se encontra do bordo, e a solucao nos demais pontos sao encontrados atraves

de simetria.

No caso de modelos contendo camadas paralelas inclinadas com um certo angulo

de mergulho, Φ, como mostrado na figura 4.4(c), os dados sao obtidos via rotacao

inversa em torno da fonte. Por exemplo, dado um ponto (x′i , z′i) no modelo atual, e

possıvel obter o ponto (xi, zi) como se o modelo estivesse na configuracao horizontal

atraves da expressao (4.16) e aplicando-se o procedimento descrito anteriormente. xi

zi

=

cos(Φ) sen(Φ)−sen(Φ) cos(Φ)

+

x′i − x f

z′i − z f

+

x f

z f

. (4.16)

56

Capıtulo 5

Exemplos

Neste capıtulo, serao apresentados os resultados e respectivas analises da metodo-

logia desenvolvida neste trabalho para os casos de macro-modelos de velocidades que

possuem solucao de referencia, visando avaliar a precisao dos resultados numericos

obtidos. Desta forma, sera possıvel avaliar a influencia dos principais parametros

da modelagem sısmica e do modelo de velocidades na estimativa dos angulos de

incidencia/reflexao em profundidade.

Inicialmente, serao testadas as tecnicas para a obtencao dos angulos de refle-

xao utilizando para tal solucoes analıticas para o angulo de mergulho/inclinacao do

macro-modelo de velocidades e solucoes semi-analıticas para a matriz de tempo de

transito correspondente, geradas pela aplicacao do procedimento descrito na secao

4.2. O proposito e avaliar as estrategias desenvolvidas para fazer a composicao dos

angulos α (MTT) e β (modelo) para a geracao dos angulos finais de reflexao, γ, assim

como expor a dificuldade de obtencao das derivadas numericas atraves dos metodos

de gradientes descritos na secao 3.3, relembrando que as mesmas sao necessarias

para a aplicacao da metodologia desenvolvida. Posteriormente, serao desenvolvi-

dos e apresentados exemplos visando avaliar as situacoes onde a matriz de tempo

de transito e obtida via propagacao do campo de ondas e pela solucao da equacao

iconal. Por fim, sera apresentado um exemplo, para um modelo de 5 camadas plana-

inclinadas, no qual sao obtidas as curvas de angulos de reflexao pelas metodologias

propostas.

5.1 Analise dos Metodos de Extracao de Gradi-

ente a Partir de Dados Semi-Analıticos

Pretende-se mostrar nesta secao a capacidade dos metodos de extracao de gradi-

entes quando estes sao aplicados sobre dados com solucao conhecida, visando avaliar

a precisao intrınseca dos mesmos. Para este proposito, serao fornecidos, como dados

57

de entrada, o modelo de velocidades, onde e conhecido o angulo formado pela reta

normal ao refletor e a vertical; e a matriz de tempo de transito semi-analıtica, da

qual e conhecido o angulo formado pela direcao do campo incidente em relacao a

vertical. Ressalta-se que a importancia destes testes, da-se, uma vez que os metodos

de gradientes serao utilizados em dados numericos, nos exemplos que seguem nas

proximas secoes, tornando-se necessario avaliar a qualidade dos resultados obtidos

por estes, quando aplicados em dados analıticos e semi-analıticos, para que se possa

ter uma ideia da suscetibilidade a erros intrınseca de cada um.

Observa-se que, para facilitar a descricao nas figuras que ilustrarao os resultados

nas proximas secoes, os metodos de obtencao de gradientes descritos na secao 3.3,

serao abreviados como segue: (i) operadores de diferencas finitas de quarta ordem,

ODF 4; (ii) interpolacao atraves das funcoes de base radial, Int RBF, e (iii) spline

cubica local de Akima, SCL Akima.

5.1.1 Extracao do Angulo do Modelo de Velocidades

Como ja dito anteriormente, com a aplicacao do operador gradiente sobre o

macro-modelo de velocidades extrai-se o angulo auxiliar, β, necessario para compor

o angulo de incidencia (subsecao 3.2.2). Em funcao disto, ha a necessidade de se

adequar o modelo para obter este angulo. A adequacao do modelo e feita aplicando

o procedimento de suavizacao mostrado na secao 2.3, com o proposito de analisar a

sua influencia nas metodologias de obtencao de gradientes.

Para exibir os resultados das analises realizadas nesta subsecao, inicialmente, sera

considerado um modelo de velocidades constituıdo por duas camadas homogeneas

separadas por um refletor em posicao inclinada. O objetivo e mostrar a habilidade

dos metodos de obtencao de gradientes para estimar os angulos de mergulho adota-

dos, de Φ = 10◦ e Φ = 20◦, quando diferentes numeros de pontos de suavizacao sao

utilizados, mantendo-se o espacamento da malha constante. Destaca-se que, para

os casos em que o refletor encontra-se em posicao plano-inclinada, estimar o angulo

formado pelo vetor gradiente e sua componente vertical e equivalente a estimar o

angulo de mergulho, isto e, implica em Φ = β ao longo de todo o refletor. Ressalta-

se tambem que este angulo de mergulho possui solucao, uma vez que sao utilizadas

funcoes especıficas para a geracao do refletor, como retas e senoides.

A tabela 5.1 mostra os parametros numericos utilizados nesta primeira analise.

As figuras 5.1(a) e 5.1(b) mostram o modelo de velocidades sem e com suavi-

zacao, respectivamente. Destaca-se que, o processo de suavizacao do modelo de

velocidades visa regularizar a curva que forma o refletor fazendo com que os “pa-

tamares” locais presentes na formacao da superfıcie refletora (figura 5.1(a)) sejam

praticamente eliminados, conforme observado na figura 5.1(b).

58

Tabela 5.1: Parametros utilizados na analise de um modelo contendo um refletorplano-inclinado.

Parametros Valor Significado/UnidadeItotal 101 Numero de pontos na direcao xJtotal 101 Numero de pontos na direcao z

h 10 Espacamento da malha (m) (h = ∆x = ∆z)nps 5 Numero de pontos de suavizacao

Φ 10 Angulo de mergulho (graus)

(a) Modelo de velocidades sem suavizacao constituıdo por um refletor com angulode mergulho de Φ = 10◦.

(b) Modelo de velocidades suavizado com 5 pontos. A reta tracejada indica aposicao analıtica do refletor/interface.

Figura 5.1: Ilustracao de modelos sem e com suavizacao.

59

Observando-se a figura 5.1(b) pode-se notar que o processo de suavizacao trans-

forma a linha de interface em uma camada contınua em torno da reta colocada ao

longo do refletor, faixa esta claramente proporcional a quantidade de pontos uti-

lizados na suavizacao. Comparando-se os graficos da figura 5.2 com aqueles das

figuras 5.3 a 5.5, nota-se que a nao suavizacao do modelo pode causar erros locais

elevados no calculo da direcao da normal do refletor. Portanto, para os demais

exemplos, o processo de suavizacao do modelo de velocidades sera um procedimento

indispensavel na obtencao das melhores aproximacoes numericas para as derivadas

(gradientes).

(a) Angulo de mergulho estimado, Φ.

(b) Erro absoluto na estimativa do angulo de mergulho.

Figura 5.2: Resultados na aproximacao do angulo de mergulho de 10◦ com o modelode velocidades sem suavizacao.

Avaliando agora a influencia do grau de suavizacao na qualidade dos resultados,

observe inicialmente a figura 5.3, a qual apresenta o angulo de mergulho analıtico

(linha reta) e aqueles estimados pelos tres metodos de gradientes ao longo do refletor

60

da figura 5.1(b). Observa-se que a curva de angulo obtida para cada metodo exibe

um comportamento caracterıstico, isto e, com oscilacoes periodicas. No caso da SCL

Akima e do ODF 4, as oscilacoes ocorrem em torno da solucao analıtica, sendo que

no primeiro, o erro maximo encontra-se em pouco mais de um grau, enquanto que

no segundo, mesmo apresentando comportamento similares ao anterior, apresenta

pontos com erro de 3, 5◦, como pode ser observado na figura 5.3(b). Por outro lado,

a Int RBF mostra um resultado mais regular e consistente, apresentando um erro

maximo em torno de um grau. Vale ressaltar que a RBF utilizada e a multiquadrica

inversa (subsecao 3.3.2) com parametro de forma ζ = 10−4, o qual, neste exemplo e

nos demais, fornece os melhores resultados.



Figura 5.3: Resultados na aproximacao do angulo de mergulho de 10◦ com o modelode velocidades suavizado com 5 pontos.

O que os resultados apresentam em comum e o fato dos tres metodos, mesmo com

a suavizacao de 5 pontos, ainda conseguirem detectar os“patamares” locais presentes

61

ao longo do refletor (observado pela periodicidade das oscilacoes), significando que

e necessario mais pontos na suavizacao do modelo para estreitamento das oscilacoes

e, consequentemente, melhora da precisao dos resultados.

Desta forma, aplicando-se agora uma suavizacao com nps = 10 pontos, a qual

regulariza ainda mais o contraste de velocidades entre as camadas, Observa-se atra-

ves da figura 5.4, que, de modo geral, houve uma reducao nos erros cometidos pelos

tres metodos na estimativa do angulo de mergulho, a menos da SCL Akima que em

algumas regioes amplificou os erros. Tambem e importante notar que, a Int RBF

apresenta o melhor resultado, com um erro de menos de 0, 5◦ enquanto que o ODF

4 e a SCL Akima exibem erros, ainda, de forma bastante oscilatoria em torno da

solucao analıtica.




62

Para uma ultima analise, utiliza-se o valor limite de nps = 20 pontos para a

suavizacao (resultados na figura 5.5). Observando-se que valores maiores do que

este poderiam causar interferencia entre os refletores, conforme descrito a frente.




Observando entao os resultados da figura 5.5, nota-se que o erro, como um todo,

foi reduzido quase pela metade em relacao ao caso anterior, com a Int RBF apre-

sentando um resultado ainda mais preciso. Embora esta reducao seja bastante sig-

nificativa, adotar um modelo com elevado numero de pontos de suavizacao, pode

ocasionar a perda de precisao na estimativa dos angulos ao longo do refletor, justifi-

cada, em modelos complexos ou mesmo simples, pela possıvel interferencia entre as

superfıcies refletoras, que contem, por exemplo, elevado grau de curvatura, a qual

pode ser eliminada pela suavizacao. Portanto, deve-se utilizar um numero de pontos

63

adequados de tal forma que o angulo de incidencia obtido no final da aplicacao da

metodologia desenvolvida nao seja gravemente alterado, conforme sera exemplificado

no ultimo exemplo desta subsecao.

Assim, para os tres casos avaliados anteriormente, ficou claro que o aumento

do numero de pontos utilizados na suavizacao do modelo melhora a estimativa dos

angulos de mergulho. A seguir, apresenta-se os resultados para o segundo angulo de

mergulho (Φ = 20◦) ilustrado na figura 5.6.

O objetivo agora e mostrar que a medida que a inclinacao do refletor aumenta,

se aproximando do angulo de 45◦, menos pontos na suavizacao do modelo sao neces-

sarios para estimar o angulo com mais precisao. Como sera verificado, isto e valido

para os tres metodos de gradientes.

Figura 5.6: Modelo de velocidades sem suavizacao constituıdo por um refletor comangulo de mergulho de Φ = 20◦.

A figura 5.7 apresenta os resultados obtidos considerando o modelo suavizado

com nps = 5 pontos. Comparando este resultado com aquele mostrado na figura

5.3, ve-se que os erros produzidos para cada metodo, de modo geral, sao menores.

Isto e justificado pela reducao dos “patamares” locais que formam o refletor, o que

leva a necessidade de menos pontos de suavizacao para estimar os angulos.

Apenas para confirmar o resultado anterior, fez-se uma suavizacao com nps = 10pontos. Na figura 5.8 e possıvel observar os resultados obtidos, os quais podem ser

comparados com aqueles da figura 5.4, de forma a verificar que ha uma reducao dos

erros estimados pelos tres metodos. Um detalhe que e observado na comparacao

das referidas figuras e que a SCL Akima apresenta mais oscilacoes e erros do que

ODF 4, o que caracteriza necessidade de mais pontos na suavizacao do modelo para

64




se conseguir um melhor resultado para este metodo. De modo geral, como sera

visto na subsecao 5.1.2 e na secao 5.2, a SCL Akima consegue obter muito bons

resultados quando o dado de entrada fornecido e a matriz de tempo de transito

semi-analıtica (subsecao 5.1.2), mas, por outro lado, mostra-se bastante sensıvel aos

dados numericos, acarretando em erros consideravelmente maiores (secao 5.2).

Ja no caso da Int RBF, como ficou claro nos resultados apresentados, esta

mostrou-se muito mais adequada na estimativa do angulo de mergulho (Φ = β),

por gerar curvas de angulo praticamente sem oscilacoes.

65




Para finalizar as analises apresentadas, neste exemplo, referentes a obtencao

de angulos de mergulho dos refletores, apresentam-se os resultados para o modelo

de velocidades da figura 5.9(a). Tal modelo visa demostrar a influencia de uma

suavizacao excessiva na precisao do calculo do angulo de mergulho.

Para tanto, o modelo da figura 5.9(a), constituıdo por dois refletores e tres ca-

madas homogeneas, foi submetido as suavizacoes de 5, 10 e 20 pontos. Ressalta-se

que, a geracao da superfıcie curva foi feita por uma senoide, o que implica que os

angulos analıticos sejam dados facilmente pela expressao:

Φ(x) = tan−1(cos(x)).

66

(a) Modelo de velocidades sem suavizacao.

(b) Modelo de velocidades suavizado com nps = 5pontos.

(c) Modelo de velocidades suavizado com nps = 10pontos.

(d) Modelo de velocidades suavizado com nps =

20 pontos.

Figura 5.9: Modelo de velocidades submetido ao processo de suavizacao.

67

A figura 5.10 apresenta os resultados numericos para o modelo da figura 5.9(b).

Como era de se esperar, uma suavizacao com 5 pontos faz com que a SCL Akima e o

ODF 4 apresentem oscilacoes ao longo da curva, enquanto que a Int RBF apresente

uma curva mais regular. Nota-se que na regiao do lado esquerdo da figura 5.10(a)

ha um afastamento maior entre as curvas numericas e teoricas, em relacao ao lado

direito, devido a interferencia ocorrida entre os refletores, ocasionada pelo processo

de suavizacao.



Figura 5.10: Resultados na aproximacao do angulo de mergulho para o modelo develocidades da figura 5.9(b).

68

Adotando-se agora uma suavizacao com 10 pontos (modelo da figura 5.9(c)),

chega-se a resultados menos oscilantes e mais precisos na parte direita do grafico,

como pode ser visto na figura 5.11. Entretanto, em relacao a solucao analıtica, e

comprovada a deterioracao dos resultados apresentados pelas curvas de angulos na

parte esquerda da figura, para os tres metodos de gradientes. Como ja observado,

este fato e ainda mais evidenciado quando e feita uma suavizacao com 20 pontos

(modelo da figura 5.9(d)), conforme apresentado na figura 5.12. Ressalta-se, nesta

figura, que os tres metodos apresentam solucoes quase praticamente identicas, in-

clusive na deteccao da interferencia ocorrida entre os refletores.



Figura 5.11: Resultados na aproximacao do angulo de mergulho para o modelo develocidades da figura 5.9(c).

69

Portanto, conclui-se que, quando o modelo de velocidades e constituıdo por refle-

tores contendo um elevado grau de curvatura, ou proximos entre si, o procedimento

de suavizacao do modelo de velocidades pode ocasionar a deterioracao dos resulta-

dos, e assim implicar na perda de precisao na estimativa dos angulos ao longo do

refletor da respectiva regiao. Sendo assim, deve-se limitar o numero maximo de pon-

tos utilizados na suavizacao, de tal forma que nao haja interferencias ou perda de

curvaturas. Dos exemplos analisados, tira-se como numero adequado a faixa dada

por 0 < nps 6 10.



Figura 5.12: Resultados na aproximacao do angulo de mergulho para o modelo develocidades da figura 5.9(d).

70

5.1.2 Extracao do Angulo da Matriz de Tempo de Transito

Nesta subsecao, o objetivo e aplicar os metodos propostos para a extracao de

gradientes, a fim de obter o angulo α a partir da matriz de tempo de transito. Tal

angulo, como ja mencionado, e utilizado para compor o angulo de reflexao em profun-

didade, conforme descrito na subsecao 3.2.1. Para tanto, serao utilizadas as matrizes

de tempo de transito semi-analıticas, obtidas mediante o procedimento apresentado

na secao 4.2, visando avaliar os referidos metodos na ausencia de ruıdos numericos,

gerados geralmente quando utilizam-se modelagens ou a solucao via equacao iconal

para obter a matriz de tempo de transito.

Assim, os exemplos utilizam as matrizes de tempo de transito mostradas na figura

5.13, correspondentes aos dois primeiros modelos de velocidades considerados na

subsecao 5.1.1, ou seja, aqueles das figuras 5.1(a) e 5.6, porem sem aplicar suavizacao.

A figura 5.14 exibe o angulo α estimado e o erro absoluto cometido por cada metodo

ao longo do refletor da figura 5.1(a), que possui angulo de mergulho de 10◦, com a

fonte localizada a 500 metros de distancia e 50 metros de profundidade em relacao

a origem do sistema de eixos mostrado no respectivo modelo.

Observa-se na figura 5.14(a) que a SCL Akima exibe um resultado bastante

apreciavel ja que os erros, presentes na curva de angulo, sao bem reduzidos, ao

contrario do ODF 4, que mostra oscilacoes e amplificacao dos erros proximo aos

dois extremos da curva. Ja no caso da Int RBF, embora a curva de angulo seja

mais suave em relacao ao ODF 4, a mesma mostra tambem um aumento dos erros

quando se aproxima dos extremos. Para se ter uma ideia geral de quao boa e a

aproximacao de cada metodo em todo modelo, a figura 5.14(c) apresenta as curvas

de nıvel referentes aos tres metodos, juntamente com a solucao semi-analıtica. E

possıvel verificar que os tres metodos conseguem estimar o angulo α em cada ponto

do modelo de forma muito satisfatoria, a menos dos pontos localizados na vizinhanca

e sobre o refletor.

Agora, considerando o caso em que o refletor tem uma inclinacao de 20◦ (modelo

da figura 5.6), e possıvel ver na figura 5.15 que, em relacao ao caso anterior, o

metodo ODF 4 apresenta muitas oscilacoes quando os angulos sao avaliados sobre o

refletor, o que se constitui em uma caracterıstica geral dos operadores de diferencas

finitas quando aplicados nas regioes onde ha fortes descontinuidades. Por outro

lado, a Int RBF mantem um certo padrao quanto a suavidade da curva de angulos,

embora a curva apresente um pequeno deslocamento em relacao a solucao semi-

analıtica. Quanto a SCL Akima, os resultados ainda sao bastantes satisfatorios para

a situacao em questao.

Assim, de modo geral, ve-se que nos dois modelos avaliados (Φ = 10◦ e 20◦)anteriormente, os tres metodos utilizados conseguiram obter os angulos com muita

71

(a) Matriz de tempo de transito semi-analıtica correspondente ao modelo dafigura 5.1(a).

(b) Matriz de tempo de transito semi-analıtica correspondente ao modelo dafigura 5.6.

Figura 5.13: Matrizes de tempo de transito semi-analıticas obtidas para dois dife-rentes modelos.

precisao nas regioes homogeneas onde nao ha problemas de descontinuidades na

matriz de tempo de transito (figuras 5.14(c) e 5.15(c)). Porem, quando aplicados

em regioes localizadas sobre ou nas proximidades dos refletores, que sao justamente

as regioes de maior interesse, o erro cometido e bastante significativo, a excecao da

SCL Akima que mostra bons resultados tanto sobre, quanto nas proximidades do

refletor.

72

(a) Angulo α estimado ao longo do refletor.

(b) Erro absoluto na estimativa do angulo α.

(c) Matriz de angulos representada por curvas de nıveis.

Figura 5.14: Resultados na aproximacao do angulo α para o modelo de velocidadesda figura 5.1(a).

73



(c) Matriz de angulos representada por curvas de nıveis.

Figura 5.15: Resultados na aproximacao do angulo α para o modelo de velocidadesda figura 5.6.

74

5.1.3 Composicao do Angulo de Incidencia/Reflexao

Para complementar os resultados das subsecoes 5.1.1 e 5.1.2, apresentam-se aqui

os graficos referentes aos angulos de reflexao em profundidade, obtidos pela com-

posicao dos angulos α (MTT) e β (modelo), utilizando-se o procedimento descrito

na subsecao 3.2.3. Vale ressaltar mais uma vez que, o objetivo desta subsecao, bem

como das duas anteriores, e aplicar o procedimento, ou melhor, os passos necessarios

para obter o angulo de incidencia em profundidade a partir de dados de entrada li-

vres de ruıdos, com o intuito de demonstrar a grande dificuldade existente no calculo

de derivadas de funcoes descontınuas utilizando processos numericos.

Neste contexto, os resultados da composicao entre os angulos, oriundos do modelo

de velocidades e da matriz de tempo de transito, sao realizados fazendo-se uma

correspondencia entre os resultados gerados por cada metodo de gradiente. Entao,

combinando as curvas de angulos procedentes do modelo com angulo de mergulho de

10◦, obtidas utilizando suavizacao com 10 pontos (figura 5.5), com aquelas oriundas

da matriz de tempo de transito da figura 5.14(a), obtem-se as curvas de angulo de

incidencia, indicadas pela figura 5.16. Da mesma forma, para o caso das figuras

5.8(a) com 5.15(a) resulta na figura 5.17 onde, neste caso, o angulo de mergulho e

de 20◦ obtido empregando suavizacao com 10 pontos.

Observando os resultados supracitados, evidencia-se que, apos o procedimento

ser aplicado, as curvas de angulos de reflexao geradas pela SCL Akima sao bas-

tante proximas da solucao de referencia, sendo este o melhor resultado entre todos

os metodos. Embora as curvas produzidas pela Int RBF nao sejam tao precisas

quanto aquelas geradas pela SCL Akima, e a que se destaca por apresentar curvas

mais suaves. Enquanto que aquelas concebidas pelo ODF 4 apresentam oscilacoes e

aumento da magnitude do erro de forma mais acentuada, a medida que se aproxima

dos extremos.

Ressalta-se que, o fato das curvas de angulos obtidas pela SCL Akima, atraves da

matriz de tempo de transito semi-analıtica gerada a partir do modelo nao suavizado,

serem as melhores concebidas entre todos os metodos, nao sera evidenciado adiante.

Quando este metodo sera submetido a matrizes de tempo de transito oriundas da

modelagem sısmica, uma vez que estes introduzem ruıdos numericos na matriz de

tempo.

Relativamente a reducao da magnitude do erro, bem como das possıveis oscilacoes

que possam ser produzidas pelos metodos de gradientes, na proxima secao, serao

avaliadas algumas estrategias inerentes ao procedimento de obtencao da matriz de

tempo de transito visando produzir curvas de angulos mais suaves.

75

(a) Angulo de incidencia, γ.

(b) Erro absoluto.

Figura 5.16: Resultados na aproximacao do angulo γ de incidencia apos composicaodo angulo α e β sobre o refletor do modelo de velocidades da figura 5.1(a).

76

(a) Angulo de incidencia, γ.

(b) Erro absoluto.

Figura 5.17: Resultados na aproximacao do angulo γ de incidencia apos composicaodo angulo α e β sobre o refletor do modelo de velocidades da figura 5.6.

77

5.2 Extracao de Angulos a Partir de Dados Nu-

mericos

O objetivo nesta secao e avaliar a precisao da tecnica proposta para obter os an-

gulos do campo incidente considerando os dados de entrada numericos. Isto significa

que a matriz de tempo de transito sera obtida atraves da propagacao do campo de

ondas (secao 2.1), mediante a aplicacao de um dos criterios descritos na subsecao

2.1.3, e tambem pela solucao da equacao iconal (secao 2.2), utilizando o MSFM

exposto na subsecao 2.2.2.2.

5.2.1 Extracao do Angulo a Partir da MTT Obtida Pela

Equacao da Onda

Igualmente ao objetivo da subsecao 5.1.2, aqui serao mostrados os resultados da

obtencao do angulo do campo incidente com a reta vertical, α, a partir da matriz de

tempo de transito oriunda da propagacao do campo de ondas via equacao da onda.

Uma vez que a qualidade de tal matriz depende dos parametros de modelagem,

torna-se importante avaliar os resultados de angulo em funcao de cada um destes.

Desta forma, os parametros avaliados sao:

1. tamanho do incremento temporal, ∆t;

2. valor da frequencia de corte, fcorte, da fonte sısmica;

3. magnitude do contraste de velocidades do modelo;

4. grau de suavizacao do modelo para obtencao da MTT.

5.2.1.1 Influencia do Tamanho de Incremento de Tempo

Para determinar a influencia do tamanho de incremento de tempo, inicialmente,

e avaliada a diferenca entre as matrizes de tempo de transito obtidas pelas equacoes

acustica tradicional e acustica nao-reflexiva. A fim de fazer esta avaliacao, serao

aplicados os dois criterios vistos na subsecao 2.1.3 (Algoritmos 01 e 02), com o

objetivo de mostrar que os resultados obtidos pelos mesmos sao quase praticamente

indiferentes em se tratando de modelos como o da figura 5.1(a) e que, portanto, para

os exemplos analisados adiante, podera ser adotada apenas uma das duas equacoes

e um dos dois algoritmos citados. Ressalta-se que, tal resultado pode nao ser valido

para modelos muito complexos e que, deste modo, aplica-se somente para os modelos

analisados neste trabalho. A tabela 5.2 mostra os parametros numericos utilizados

nas modelagens.

78

Tabela 5.2: Parametros utilizados nas modelagens do modelo da figura 5.1(a).



Φ 10 Angulo de mergulho (graus)x f 500 Posicao da fonte na direcao horizontal (m)z f 50 Posicao da fonte na direcao vertical (m)∆t 10−3 Incremento temporal (s)

ttotal 0,8 Tempo total de analise (s)fcorte 30 Frequencia de corte (Hz)

A figura 5.18(a) mostra uma comparacao entre as matrizes de tempo de transito,

atraves de isocronas, oriundas da solucao de referencia e da propagacao do campo de

ondas mediante a equacao acustica tradicional. Destaca-se que, foram utilizados os

criterios da amplitude maxima e amplitude maxima nas proximidades da primeira

quebra (Algoritmos 01 e 02), os quais sao abreviados como AMPM e AMPMQ.

Similarmente, a figura 5.18(b) exibe as isocronas obtidas mediante a equacao acustica

nao-reflexiva com os mesmo criterios.

Comparativamente a solucao de referencia observa-se, nas figuras 5.18(a) e

5.18(b) que as isocronas obtidas atraves da propagacao do campo de ondas para

cada uma das equacoes, mantem a mesma conformacao ao longo do modelo de velo-

cidades, o que indica que a direcao de propagacao da onda e praticamente a mesma

da solucao de referencia. Nas mesmas figuras citadas, evidencia-se que o resultado

da equacao nao-reflexiva esta sobre aquele da equacao tradicional.

Ressalta-se que, a solucao do campo de ondas pela equacao nao-reflexiva apre-

senta um custo computacional relativamente maior quando comparado com aquele

da equacao acustica tradicional, uma vez que a primeira contem um termo a mais

em relacao a segunda. Em funcao disto e tambem do fato que os modelos de ve-

locidades aqui utilizados nao contem grandes complexidades geometricas, adotou-se

a equacao tradicional nos testes realizados neste trabalho. A mesma justificativa

e utilizada em relacao a adocao do criterio da AMPM em detrimento a AMPMQ.

Sendo assim, seguem os resultados para a equacao acustica tradicional empregando

o criterio da amplitude maxima.

As figuras 5.19, 5.20 e 5.21 mostram o resultado obtido para os parametros descri-

tos na tabela 5.2, destacando-se que, especificamente quanto ao resultado da figura

5.19, e importante compara-lo com aquele da figura 5.14, que mostra os angulos ob-

tidos a partir da matriz de tempo de transito semi-analıtica. Da observacao destas

duas figuras, fica bastante perceptıvel que a estimativa dos angulos feita pela Int

79

RBF mostra uma curva de angulo relativamente suave e bastante semelhante entre

elas. Em contrapartida, observa-se um erro bem mais acentuado produzido pelos

gradientes numericos calculados pelos outros dois metodos, sendo a SCL Akima bas-

tante sensıvel as mudancas ocorridas na matriz de tempo de transito devido ao uso

da modelagem sısmica, mostrando uma curva de angulo bem mais dispersa. Tam-

bem nota-se que a curva gerada por ODF 4 apresenta muitas oscilacoes, chegando a

mostrar o mesmo comportamento que a curva obtida pela SCL Akima em algumas

regioes ao longo do refletor.

Na figura 5.19, alem de serem mostrados os angulos obtidos pelos metodos de

gradientes, tambem exibe-se a curva de angulo determinada atraves da aplicacao do

ODF 4 durante a propagacao do campo de ondas, isto e, obtida pelo gradiente da

frente de ondas, conforme procedimento descrito na subsecao 3.3.4. Confrontando-se

esta curva com aquelas geradas por ODF 4 e a SCL Akima, nota-se que a mesma

apresenta um resultado ligeiramente melhor, ja que produziu curvas mais suaves.

No geral, ve-se que a matriz de tempo de transito obtida pela modelagem levou

os metodos de gradientes a produzirem uma matriz de angulos com muitos ruıdos ao

longo do modelo de velocidades, como pode ser observado na figura 5.19(c). Ainda

assim, cabe destacar que o gradiente numerico produzido pela Int RBF mostrou um

melhor resultado na regiao acima e sobre o refletor, como e indicado pela figura 5.20

que mostra os angulos estimados para a segunda e primeira linhas acima do refletor.

Por outro lado, nas proximidades abaixo do refletor (figura 5.21), isto e, para os

angulos de refracao, as curvas sao melhor computadas quando obtidas diretamente

pela frente de onda, ja que esta exibe uma curva sem oscilacoes e mais proxima

da curva de referencia. Ressalta-se que os erros mostrados nas figuras 5.20 e 5.21

ocorrem somente nas regioes proximas ao refletor.

Assim, para reduzir o espalhamento e as oscilacoes nas curvas de angulos,

realizam-se testes reduzindo o incremento ∆t, ate encontrar-se um valor adequado

para o mesmo, que foi de ∆t = 25×10−6, o qual e chamado de ∆t otimo. Por questoes

de simplicidade e concisao entao, visando somente demonstrar a respectiva melhora

nas curvas de angulos, repetiu-se a modelagem feita anteriormente utilizando este

∆t. As figuras 5.22, 5.23 e 5.24 exibem os resultados obtidos.

Comparando os resultados entre as respectivas figuras provenientes da primeira

modelagem com aquelas desta segunda, para as quatro curvas de angulos, na regiao

sobre e acima do refletor, ve-se que a Int RBF e uma importante metodologia na de-

terminacao das derivadas numericas, uma vez que esta mostra-se quase praticamente

independe do incremento de tempo, possibilitando a geracao de curvas mais suaves

e proximas da solucao de referencia. Ja nos demais metodos, a obtencao de cur-

vas suaves e sem oscilacoes depende do valor de incremento adotado na modelagem

sısmica, o que fica evidenciado nos resultados apresentados.

80

Vale destacar, na comparacao entre os resultados produzidos atraves das duas

modelagens, que para os resultados obtidos utilizando o ∆t otimo, em regioes afasta-

das do refletor, foi melhorado de forma bastante significativa por todos os metodos,

como pode ser observado ao confrontar as figuras 5.19(c) e 5.22(c).

(a) Isocronas obtidas pelos criterios AMPM e AMPMQ via equacao acusticatradicional.

(b) Isocronas obtidas pelos criterios AMPM e AMPMQ via equacao acusticanao-reflexiva.

Figura 5.18: Comparacao entre matrizes de tempo de transito obtidas com diferentesequacoes e criterios.

81



(c) Matriz de angulos representado por curvas de nıveis.

Figura 5.19: Resultados na aproximacao do angulo α (em graus) quando a matrizde tempo de transito e obtida com ∆t = 10−3 s.

82

(a) Angulo α estimado para a segunda linha acima do refletor.


(c) Angulo α estimado para a primeira linha acima do refletor.

(d) Erro absoluto na estimativa do angulo α.

Figura 5.20: Resultados na aproximacao do angulo α em regioes acima do refletor.

83

(a) Angulo α estimado para a primeira linha abaixo do refletor.


(c) Angulo α estimado para a segunda linha abaixo do refletor.


Figura 5.21: Resultados na aproximacao do angulo α em regioes abaixo do refletor.

84




Figura 5.22: Resultados na aproximacao do angulo α (em graus) quando a matrizde tempo de transito e obtida com ∆t = 25 × 10−6 s.

85





Figura 5.23: Resultados na aproximacao do angulo α em regioes acima do refletor.

86





Figura 5.24: Resultados na aproximacao do angulo α em regioes abaixo do refletor.

87

5.2.1.2 Influencia da Frequencia de Corte

Um outro parametro da modelagem que influencia nas curvas de angulo e a

frequencia de corte utilizada para controlar o comprimento de onda emitido pela

fonte sısmica (subsecao 2.1.1.4). Para avaliar esta influencia, realizou-se tres mode-

lagens com os respectivos valores de frequencia de 30, 60 e 120Hz. As modelagens

foram feitas no modelo da figura 5.25, o qual contem o mesmo numero de pontos

utilizados no modelo da figura 5.1, mas com o espacamento da malha reduzido para

4m devido ao criterio de dispersao numerica, conforme descrito na subsecao 2.1.1.3.

A tabela 5.3 descreve os parametros utilizados.

Tabela 5.3: Parametros utilizados na modelagem do modelo da figura 5.25.



Φ 10 Angulo de mergulho (graus)x f 200 Posicao da fonte na direcao horizontal (m)z f 20 Posicao da fonte na direcao vertical (m)∆t 25 × 10−6 Incremento temporal (s)

ttotal 0,6 Tempo total de analise (s)fcorte 30, 60 e 120 Variacao na frequencia de corte (Hz)

Figura 5.25: Modelo de velocidades para analise da frequencia de corte, fcorte.

O primeiro resultado exibido, corresponde a frequencia de corte de 30Hz nas

figuras 5.26, 5.27 e 5.28. Este resultado e confrontado com aqueles das figuras 5.29,

88

5.30 e 5.31, referentes a frequencia de 60Hz, e aqueles das figuras 5.32, 5.33 e 5.34

relativas a frequencia de 120Hz. As observacoes podem ser agrupadas como segue:

(i) quanto a precisao das curvas obtidas sobre o refletor:

E possıvel observar nas figuras 5.26, 5.29 e 5.32 que, o aumento da frequencia

de corte, fez com que a Int RBF produzisse curvas de angulos com maior

precisao, pois e notado uma reducao na magnitude do erro. No entanto, tal

resultado nao e evidenciado para as curvas geradas por ODF 4, SCL Akima

e para as curvas obtidas atraves da frente de onda, sendo que, nos resultados

apresentados por estes, a magnitude do erro foi ampliada, apos o acrescimo

da frequencia, como pode ser visto nas respectivas figuras que indicam o erro

absoluto. Destaca-se ainda, o surgimento de oscilacoes nas curvas geradas por

estes tres metodos para o ultimo valor de frequencia adotado.

(ii) quanto a precisao das curvas obtidas na regiao acima do refletor:

Contemplando as figuras 5.27, 5.30 e 5.33, respectivamente, ve-se que todos

os metodos originaram curvas suficientemente satisfatorias, apresentando uma

magnitude no erro cada vez menor com o aumento da frequencia. Salienta-se

os bons resultados gerados pela Int RBF na frequencia de corte de 120Hz.

(iii) quanto a precisao das curvas obtidas na regiao abaixo do refletor:

Consultando as respectivas figuras 5.28, 5.31 e 5.34, verifica-se que a magnitude

do erro produzido pelos metodos ODF 4, SCL Akima e pela frente de ondas e

progressivamente reduzida com o acrescimento da frequencia. Por outro lado,

isto nao e comprovado nas curvas originadas a partir da Int RBF.

Por fim, destaca-se que nas figuras 5.26(c), 5.29(c) e 5.32(c), a medida que o valor

da frequencia aumenta, ocorre a contracao e amplificacao das oscilacoes presentes

nas curvas de nıvel nas proximidades do refletor, assinalando que o erro esta mais

concentrado naquela regiao. Tal evento, e razoavel de se esperar, uma vez que o

acrescimo da frequencia implica em reducao do comprimento de onda, fato este

evidenciado pela relacao λ = C/ fcorte. Portanto, pode-se intuir que um aumento

ainda maior na frequencia, poderia contrair ainda mais as oscilacoes. Entretanto,

este procedimento provoca uma diminuicao no valor de incremento de tempo e no

espacamento da malha, e consequentemente um acrescimo no custo computacional

de forma consideravel inviabilizando a estrategia.

Assim, como de modo geral, nao se observa um ganho efetivo de precisao nos

resultados nas vizinhancas proximas ou sobre o refletor com o aumento da frequen-

cia, conclui-se que a melhor opcao e adotar a frequencia de 30Hz para as analises

numericas, uma vez que este representara um menor custo computacional (∆t, ∆x e

∆z maiores) correspondente, inclusive, ao valor comumente utilizado em geofısica.

89




Figura 5.26: Resultados na aproximacao do angulo α (em graus) quando a matrizde tempo de transito e obtida com fcorte = 30Hz.

90





Figura 5.27: Resultados na aproximacao do angulo α em regioes acima do refletor,utilizando fcorte = 30Hz.

91





Figura 5.28: Resultados na aproximacao do angulo α em regioes abaixo do refletor,utilizando fcorte = 30Hz.

92





93






94






95



(c) Matriz de angulos representado por curvas de nıvel.


96






97






98

5.2.1.3 Influencia do Contraste de Velocidades

Para avaliar a influencia da magnitude do contraste de velocidades existente entre

as camadas, realizou-se duas modelagens, como ja dito, com a equacao acustica

tradicional via aplicacao do criterio AMPM, no mesmo modelo da figura 5.1(a),

utilizando-se o ∆t otimo e a frequencia de corte de 30Hz obtidos nas subsecoes

anteriores. A unica mudanca feita no modelo foi alterar a velocidade da camada

inferior de 2000m/s para 1550m/s e, posteriormente, para 3000m/s, criando-se um

contraste baixo e alto, respectivamente.

Destaca-se que, como ja foi observado nos resultados anteriores, que a obtencao

do angulo α nas regioes homogeneas do modelo de velocidades e feita com bastante

precisao, entao, e natural que para um contraste baixo entre as camadas os angulos

sobre o refletor sejam melhor estimados, como pode ser visto na figura 5.35. Por

outro lado, quando a matriz de tempo e obtida a partir de um contraste mais elevado,

a determinacao dos angulos sobre e nas proximidades do refletor apresenta resultados

menos precisos. Neste caso a matriz de tempo possui uma forte descontinuidade

gerada pela grande diferenca de velocidades entre as camadas, o que acarreta em

forte descontinuidade na matriz de angulo, sobre e nas proximidades do refletor,

como pode ser visto na figura 5.36.

Em funcao destes resultados, ve-se que a geracao das curvas de angulos sao

bastante difıceis de serem avaliadas, sobre e nas vizinhancas do refletor, quando

a matriz de tempo de transito e obtida de um modelo de velocidades dotado de

alto contrastes entre as camadas. Destaca-se que isto se trata de um problema

numerico inerente a todos os metodos na obtencao das derivadas que compoem o

vetor gradiente.

Como consequencia, conclui-se que para modelos contendo altos contrastes, deve-

se utilizar o procedimento de suavizacao, para que nao se produzam grandes des-

continuidades na matriz de tempo de transito, e assim melhorem as estimativas

das curvas de angulos nas regioes de maior interesse. Na proxima subsecao e ana-

lisada esta questao do numero de pontos que deve ser utilizado neste processo de

suavizacao.

99




Figura 5.35: Resultados na aproximacao do angulo α para um contraste baixo develocidades entre as camadas.

100




Figura 5.36: Resultados na aproximacao do angulo α para um contraste alto develocidades entre as camadas.

101

5.2.1.4 Influencia do Grau de Suavizacao do Modelo

Em funcao dos resultados verificados na subsecao 5.2.1.3, envolvendo um modelo

simples de velocidades, pode-se observar que, para os angulos serem melhor esti-

mados, e necessario que haja um gradiente suave entre as camadas que delimitam

a superfıcie refletora, evidenciado pelos resultados do modelo com baixo contraste

entre as camadas. Entao, nesta subsecao e avaliada a influencia da quantidade de

pontos usados na suavizacao do modelo para obter a matriz de tempo de transito

durante a modelagem. Para tanto, foram feitas duas modelagens de forma ana-

loga ao caso da subsecao anterior (5.2.1.3), mas utilizando o modelo de velocidades

suavizado com 5 e 10 pontos.

As figuras 5.37, 5.38 e 5.39 mostram as estimativas nos angulos obtidos atraves

da aplicacao dos metodos de gradientes, sobre a matriz de tempo de transito obtida

pela propagacao do campo de ondas no modelo de velocidades da figura 5.1(b),

suavizado com 5 pontos. Comparando estes resultados com aqueles das figuras 5.22,

5.23 e 5.24, nota-se uma significativa reducao na ordem do erro das curvas geradas

pelos metodos, exceto a Int RBF que, em um primeiro momento, mostrou-se menos

sensıvel a suavizacao em relacao aos demais metodos. Destaca-se que, nas regioes

abaixo do refletor, houve um aumento da magnitude do erro nas curvas geradas

por ODF 4, SCL Akima e aquela obtida pela frente de onda. Tal fato, tambem e

evidenciado para o caso das curvas obtidas quando se utiliza a suavizacao com 10

pontos, inclusive para aquelas geradas pela Int RBF, como observado na figura 5.42.

Comparativamente aos resultados obtidos com a suavizacao de 5 pontos, nas

figuras 5.40 e 5.41, percebe-se que, de modo geral, houve uma reducao no erro

produzido por todos os metodos, sobretudo na regiao sobre o refletor quando utiliza-

se 10 pontos. Tambem e notado que, usando 10 pontos na suavizacao, as curvas

de angulos, a princıpio, tendem para um certo limite onde elas se equiparam, ou

seja, as mesmas se tornam praticamente iguais para os tres metodos de extracao de

gradientes. Isto caracteriza que ha um limite para o numero de pontos utilizados

na suavizacao, a partir do qual os resultados produzidos pelos tres metodos de

gradientes nao mostram praticamente nenhuma distincao. Este fato tambem sera

observado na subsecao 5.2.2 onde a matriz de tempo de transito sera obtida pela

solucao da equacao iconal.

Assim, como o objetivo e determinar os melhores parametros que propiciem

curvas de angulos mais suaves e mais proximas da solucao de referencia sobre o

refletor, passa-se a adotar 10 pontos para a suavizacao do modelo, em se tratando

da obtencao da matriz de tempo de transito originada a partir da modelagem sısmica.

Tal valor, sera utilizado na geracao das curvas de angulos de incidencia que serao

mostradas na subsecao 5.2.3.1.

102




Figura 5.37: Resultados na aproximacao do angulo α, para a MTT obtida pelamodelagem sısmica a partir do modelo de velocidades suavizado com 5 pontos.

103





Figura 5.38: Resultados na aproximacao do angulo α em regioes acima do refletor,para a MTT obtida pela modelagem sısmica a partir do modelo de velocidadessuavizado com 5 pontos.

104





Figura 5.39: Resultados na aproximacao do angulo α em regioes abaixo do refletor,para a MTT obtida pela modelagem sısmica a partir do modelo de velocidadessuavizado com 5 pontos.

105




Figura 5.40: Resultados na aproximacao do angulo α, para a MTT obtida pelamodelagem sısmica a partir do modelo de velocidades suavizado com 10 pontos.

106





Figura 5.41: Resultados na aproximacao do angulo α em regioes acima do refletor,para a MTT obtida pela modelagem sısmica a partir do modelo de velocidadessuavizado com 10 pontos.

107





Figura 5.42: Resultados na aproximacao do angulo α em regioes abaixo do refletor,para a MTT obtida pela modelagem sısmica a partir do modelo de velocidadessuavizado com 10 pontos.

108

5.2.2 Extracao do Angulo a Partir da matriz de tempo de

transito Obtida Pela Equacao Iconal

Nesta subsecao, o objetivo e mostrar os resultados obtidos pela aplicacao dos

metodos de gradientes sobre a matriz de tempo de transito oriunda da solucao da

equacao iconal pelo MSFM, para a obtencao de angulos de reflexao em profundidade.

A ideia aqui e apresentar os resultados para o modelo de velocidades da figura 5.1

considerando a mesma posicao da fonte daquela dada na tabela 5.2.

Os primeiros resultados obtidos sao aqueles mostrados nas figuras 5.43 e 5.44,

considerando a matriz de tempo de transito gerada pelo MSFM a partir do modelo de

velocidades sem suavizacao. Nestes resultados, nota-se que os angulos acima e sobre

o refletor sao melhor estimados em relacao aqueles obtidos pela modelagem sısmica

com o modelo de velocidades sem suavizacao (observar figuras 5.22, 5.23 e 5.24).

Ve-se tambem que a curva de angulo que mais se aproxima da solucao de referencia

e aquela gerada pela SCL Akima, a menos dos saltos ocorridos exatamente onde

esta o degrau presente na formacao do refletor. Como ja observado nos resultados

anteriores, tal situacao tambem acontece nas curvas geradas pelo ODF 4.

Tambem e possıvel notar que os angulos, para uma linha acima do refletor (figura

5.44(a)), sao muito bem estimados pelos tres metodos de gradientes, com amplo

destaque para as curvas geradas pela SCL Akima e o ODF 4. Entretanto, quando se

trata das curvas geradas para os angulos situados uma linha abaixo do refletor (figura

5.44(c)), as curvas apresentam uma oscilacao em torno da solucao de referencia,

exceto para a curva gerada pela Int RBF. Entao, como ja foi visto na subsecao

5.2.1.4, para melhorar a solucao e tornar as curvas de angulos mais suaves e sem

oscilacao pode-se similarmente aplicar uma suavizacao no modelo antes de aplicar o

MSFM. Assim, as figuras 5.45 e 5.46 exibem os resultados no caso em que o modelo

de velocidades foi suavizado com 3 pontos. Delas observa-se um fato, ja comentado

na subsecao 5.2.1.4, que refere-se a tendencia apresentada pelas curvas de angulos a

se aproximarem a tal ponto que praticamente nao ha distincao entre os resultados

gerados pelos tres metodos de gradientes. Tal fato tambem e evidenciado nas figuras

5.47 e 5.48, que exibem a caracterizacao do mesmo acontecimento para o caso de

uma suavizacao com 5 pontos.

Assim, como o valor de suavizacao de 5 pontos possibilitou a geracao de curvas de

angulos mais suaves, ele sera utilizado na geracao das curvas de angulos de reflexao

que serao apresentadas na subsecao 5.2.3.2.

109



(c) Matriz de angulos representados por curvas de nıveis.

Figura 5.43: Resultados na aproximacao do angulo α, para a MTT obtida peloMSFM a partir do modelo de velocidades sem suavizacao.

110

(a) Angulo α estimado para a primeira linha acima do refletor.


(c) Angulo α estimado para a primeira linha abaixo do refletor.


Figura 5.44: Resultados na aproximacao do angulo α em regioes acima e abaixodo refletor, para a MTT obtida pelo MSFM a partir do modelo de velocidades semsuavizacao.

111




Figura 5.45: Resultados na aproximacao do angulo α, para a MTT obtida peloMSFM a partir do modelo de velocidades suavizado com 3 pontos.

112





Figura 5.46: Resultados na aproximacao do angulo α em regioes acima e abaixo dorefletor, para a MTT obtida pelo MSFM a partir do modelo de velocidades suavizadocom 3 pontos.

113




Figura 5.47: Resultados na aproximacao do angulo α, para a MTT obtida peloMSFM a partir do modelo de velocidades suavizado com 5 pontos.

114





Figura 5.48: Resultados na aproximacao do angulo α em regioes acima e abaixo dorefletor, para a MTT obtida pelo MSFM a partir do modelo de velocidades suavizadocom 5 pontos.

115

5.2.3 Aplicacao para um Modelo com Varias Camadas

Neste ultimo exemplo, objetiva-se apresentar diretamente as curvas de angulos

de incidencia/reflexao para as matrizes de tempo de transito obtidas, tanto pela

modelagem sısmica quanto pela solucao da equacao iconal, a partir de um modelo

de velocidades contendo varias camadas. Especificamente, o proposito e avaliar as

curvas sobre os refletores a medida que estes estao localizados em maior profundi-

dade, isto e, localizam-se mais distantes do ponto de aplicacao da fonte. Para tanto,

serao utilizados os valores dos parametros que mais influenciaram significativamente

na melhor estimativa dos angulos oriundos da matriz de tempo de transito, como foi

visto nas subsecoes 5.2.1 e 5.2.2. Em suma, serao mostradas as curvas de angulos

de reflexao apos a composicao dos angulos α e β, igualmente ao que foi feito na

subsecao 5.1.3.

As curvas de angulos de incidencia/reflexao serao apresentadas sobre os 4 refle-

tores presentes no modelo da figura 5.49, os quais estao dispostos em forma planos

inclinados com angulo de mergulho de 10◦ e separados por 5 camadas homogeneas.

Figura 5.49: Modelo de velocidades utilizado para obtencao das curvas de incidencia.

116

5.2.3.1 Angulo de Incidencia Oriundo da Propagacao de Ondas

Para a obtencao das curvas de angulos de incidencia utilizando a matriz de

tempo de transito gerada atraves da propagacao do campo de ondas, a tabela 5.4

sintetiza os parametros empregados. Na figura 5.50, e possıvel observar o modelo,

correspondente aquele da figura 5.49, suavizado com o valor de nps = 10 pontos.

Ressalta-se que os angulos (β) oriundos do modelo foram obtidos a partir daquele

mostrado na figura 5.50.

As figuras 5.51 e 5.52 exibem as curvas de angulos obtidas ao longo dos refletores,

sendo que, na primeira figura, sao mostradas as curvas para os refletores 1 e 2,

enquanto na segunda, para os refletores 3 e 4.

A partir dos resultados apresentados, ve-se que as curvas de incidencia geradas

por todos os metodos de gradientes sao muito proximas da solucao de referencia,

destacando-se que as mesmas apresentam-se bastante suaves ao longo de cada refle-

tor. Alem disso, cabe salientar que aquelas determinadas pela Int RBF sao levemente

mais regulares. Como caracterıstica geral, evidencia-se que, a medida que a profun-

didade aumenta, a magnitude do erro de todas as curvas diminui. Tal caracterıstica

e devida ao fato de, quanto mais distante da fonte estiver o refletor, mais plana e a

configuracao da frente de ondas nas proximidades deste, o que implica na geracao

de curvas de tempo (isocronas) mais planas e suaves naquela regiao.

As curvas finais, apresentadas nas figuras supracitadas, sao os melhores resulta-

dos obtidos atraves da metodologia desenvolvida, a qual depende diretamente dos

metodos estimadores de gradientes que, por sua vez, dependem da qualidade dos

dados de entrada fornecidos.

117

Tabela 5.4: Parametros utilizados na modelagem sısmica do modelo de velocidadesda figura 5.50.



Φ 10 Angulo de mergulho (graus)x f 500 Posicao da fonte na direcao horizontal (m)z f 20 Posicao da fonte na direcao vertical (m)∆t 25 × 10−6 Incremento temporal (s)

ttotal 1,2 Tempo total de analise (s)fcorte 30 Variacao na frequencia de corte (Hz)

Figura 5.50: Modelo de velocidades suavizado com 10 pontos, utilizado na obtencaoda MTT pela modelagem sısmica.

118

(a) Angulo γ estimado sobre o refletor 1.

(b) Erro absoluto na estimativa do angulo γ.

(c) Angulo γ estimado sobre o refletor 2.

(d) Erro absoluto na estimativa do angulo γ.

Figura 5.51: Resultados na aproximacao do angulo γ sobre os refletores 1 e 2.

119






120

5.2.3.2 Angulo de Incidencia Oriundo da Solucao Iconal

Para a matriz de tempo de transito originada a partir da solucao da equacao

iconal, a figura 5.53 apresenta o mesmo modelo de velocidades visto anteriormente,

mas suavizado com 5 pontos, no qual o MSFM foi aplicado. Observa-se que os

angulos originados do modelo foram obtidos da mesma forma descrita na subsecao

anterior.

As figuras 5.54 e 5.55 indicam os resultados obtidos, evidenciando-se aqui a

mesma conclusao feita na subsecao 5.2.3.1, exceto pelo fato de as curvas apresenta-

rem suaves oscilacoes nas proximidades da solucao de referencia.

Destaca-se, quanto a suavidade e precisao das curvas de angulos, que a obtencao

das matrizes de tempo de transito pelo MSFM e mais adequada, em relacao aquelas

obtidas pela modelagem sısmica, a qual requer um modelo suavizado com numero

de pontos superior ao primeiro, assim como apresenta dependencia de parametros

inerentes a modelagem que influem diretamente na qualidade da matriz de tempo.

Figura 5.53: Modelo de velocidades suavizado com 5 pontos, utilizado na obtencaoda MTT pelo MSFM.

121






122






123

Capıtulo 6

Conclusoes e Trabalhos Futuros

Neste trabalho, foi apresentada uma metodologia desenvolvida para extracao de

angulos de incidencia/reflexao em profundidade, a partir dos angulos obtidos atraves

da aplicacao do operador gradiente no macro-modelo de velocidades e na matriz

de tempo de transito. Utilizaram-se duas abordagens diferentes para a geracao das

matrizes de tempo: modelagem sısmica da equacao completa da onda, por diferencas

finitas, e solucao direta da equacao iconal, atraves do Multistencils Fast Marching

Method.

Relativamente a obtencao dos gradientes necessarios para a aplicacao da me-

todologia proposta, foram adotados os operadores de diferencas finitas de quarta

ordem, a interpolacao atraves das funcoes de base radial e a spline cubica local de

Akima. Ressalta-se que foi adicionalmente desenvolvido um esquema numerico para

a geracao do gradiente diretamente da frente de ondas propagada, correspondente

aquele obtido da matriz de tempo de transito, aplicando-se para tanto os operadores

de diferencas finitas de quarta ordem ao campo de pressoes, durante a propagacao

do sinal sısmico emitido pela fonte.

A partir dos resultados evidenciados ao longo deste trabalho, viu-se que o pro-

cedimento de suavizacao do modelo de velocidades e uma importante ferramenta

que deve ser utilizada para obtencao das derivadas numericas. Quando tal procedi-

mento e empregado, todos os metodos de gradientes geram curvas de angulos, tanto

do modelo quanto da matriz de tempo de transito, mais suaves e quase sem osci-

lacoes. Entretanto, deve-se utilizar um numero de pontos adequados, de tal forma

que nao ocorra interferencia entre os diferentes refletores que formam o modelo,

assim como perda de curvaturas, nos casos em que os refletores sao formados por

superfıcies muito curvas.

As funcoes de base radial mostram-se como uma notavel metodologia na aquisi-

cao dos gradientes, sobretudo naqueles obtidos do macro-modelo de velocidades, pois

sao capazes de gerar curvas com pouquıssimas oscilacoes e mais contınuas, compara-

tivamente aos demais metodos. Destaca-se que, para se conseguir tal suavidade nas

124

curvas produzidas pelos outros metodos, faz-se necessario aplicar o procedimento de

suavizacao ao modelo.

A aplicacao da spline cubica local de Akima nas matrizes de tempo de transito

obtidas pelo Multistencils Fast Marching Method a partir do modelo de velocidades

sem suavizacao, consegue estimar as curvas de angulo com muita precisao, o que

nao ocorre com as demais metodologias. Igual observacao e evidenciada quando da

adocao da matriz de tempo de transito semi-analıtica, obtida pelo procedimento im-

plementado neste trabalho. Porem, tal evento nao e observado quando e empregado

a matriz de tempo oriunda da modelagem sısmica, devido a introducao de ruıdos

numericos na mesma.

No que tange aos parametros da modelagem sısmica e do macro-modelo de velo-

cidades na estimativa de curvas de angulo mais precisas e com maior suavidade, as

seguintes conclusoes e consideracoes sao pertinentes:

• influencia do incremento temporal: o valor do incremento de tempo, ∆t, ado-

tado nas modelagens, para a obtencao das matrizes de tempo de transito, e

um parametro muito importante para a geracao de curvas de angulo, espe-

cialmente, para os metodos baseados nos operadores de diferencas finitas e

para a spline de Akima. Entretanto, no caso das funcoes de base radial, estas

possibilitam a geracao de curvas mais suaves e sem oscilacoes quase indepen-

dentemente do incremento temporal adotado.

• influencia da frequencia de corte da fonte sısmica: embora o aumento da

frequencia produza contracao das oscilacoes presentes nas curvas de nıvel da

matriz de angulo nas proximidades do refletor, de modo geral, nao foi obser-

vado um ganho efetivo de precisao nos resultados gerados pelas metodologias

de gradientes, a excecao daqueles obtidos pelas funcoes de base radial, para

as regioes sobre e acima do refletor. Destaca-se que o acrescimo da frequencia

de corte implica em reducao do comprimento de onda. Portanto, pode-se in-

tuir que um aumento ainda maior na frequencia, poderia contrair ainda mais

as oscilacoes. Entretanto, tal procedimento provoca uma diminuicao no valor

de incremento de tempo (∆t) e no espacamento da malha (∆x e ∆z), e con-

sequentemente um acrescimo no custo computacional de forma consideravel,

inviabilizando a aplicacao da metodologia desenvolvida.

• influencia do contraste de velocidades: a diferenca de velocidades entre as

camadas e um fator muito importante na determinacao dos gradientes, nas

proximidades e sobre o refletor, pois esta diferenca representa, na verdade, o

grau de descontinuidade existente naquela regiao. Salienta-se que para um

contraste elevado a geracao das curvas de angulos torna-se mais complicada,

125

apresentando erros mais acentuados. Isto se trata de um problema numerico

inerente a todos os metodos aqui estudados para a obtencao das derivadas que

compoem o vetor gradiente.

• influencia do grau de suavizacao do modelo: diretamente relacionado ao con-

traste de velocidades entre as camadas, esta o numero de pontos utilizados

no processo de suavizacao, tanto para obter os gradientes do modelo quanto

da matriz de tempo de transito. Assumindo um contraste de velocidades de

500m/s (1500m/s para 2000m/s), foi observado que usando 10 e 5 pontos para

a obtencao da matriz de tempo de transito pela modelagem sısmica e pelo

Multistencils Fast Marching Method, respectivamente, as curvas de angulos se

tornaram suaves e praticamente iguais. Isto ocorre para todos os metodos de

extracao de gradientes abordados neste trabalho. Portanto, ha um limite para

o numero de pontos utilizados na suavizacao, a partir do qual os resultados

produzidos pelos metodos de gradiente sao praticamente iguais.

Enfatiza-se, quanto a suavidade e precisao das curvas de angulos, que a obtencao

das matrizes de tempo de transito pelo Multistencils Fast Marching Method e mais

adequada, em relacao aquelas obtidas pela modelagem sısmica, a qual requer um

modelo suavizado com numero de pontos superior ao primeiro, assim como apre-

senta dependencia de parametros inerentes a modelagem que influem diretamente

na qualidade da matriz de tempo.

Destaca-se ainda que, apesar de nao terem sido apresentados resultados de custo

computacional para a obtencao das matrizes de tempo de transito, foi observado que

o Multistencils Fast Marching Method, como ja e conhecido na literatura, apresenta

custo muito inferior as solucoes baseadas na propagacao via equacao da onda.

Em suma, destaca-se que, como caracterıstica da propagacao de ondas, a qual

e refletida na matriz de tempo de transito, o fato de que a precisao das curvas de

angulos se tornar cada vez maior a medida que ocorre o aumento da profundidade,

devido a configuracao mais plana das isocronas.

Sugere-se como trabalhos futuros duas abordagens, a saber:

(i) quanto a melhoria da precisao dos resultados obtidos:

1. Utilizacao de novas metodologias de obtencao de matrizes de tempo de

transito, como o Fast Iterative Method (FIM) (JEONG e WHITAKER,

2008).

2. Avaliacao de outras metodologias de obtencao dos parametros de forma

para as funcoes de base radial (FASSHAUER, 2002; KANSA, 1990;

RIPPA, 1999), assim como outras funcoes base (LIU, 2002; SCHABACK,

2000).

126

3. Utilizacao do metodo dos mınimos quadrados moveis aproximados

(MLSA - Moving Least-Squares Approximation) (LANCASTER e SAL-

KAUSKAS, 1981) para avaliacao dos gradientes.

4. Utilizacao de Galerkin descontınuo (JOHNSON, 1995) para obter os gra-

dientes em regioes com altos contrastes de impedancia.

5. Utilizacao de transformadas para avaliacao dos gradientes;

(ii) quanto a aplicacao da metodologia desenvolvida:

1. Avaliacao da metodologia em modelos mais complexos, utilizando tracado

de raios como solucao de referencia.

2. Desenvolvimento de um procedimento de “mute” automatico, utilizando-

se os angulos em profundidade para eliminacao de ruıdos nas imagens

produzidas pela RTM.

3. Construcao de curvas de AVA.

4. Extensao da metodologia desenvolvida para o meio elastico.

5. Construcao de secoes de impedancia de reflexao.

6. Construcao de secoes de angulo comum.

7. Migracao no domınio do angulo.

127

Referencias Bibliograficas

AKIMA, H., 1970, “A new method of interpolation and smooth curve fitting based

on local procedures”, Journal of the ACM, v. 17, n. 4, pp. 589–602.

AKIMA, H., 1974a, “A method of bivariate interpolation and smooth surface fitting

based on local procedures”, Communications of the ACM, v. 17, n. 1,

pp. 18–20.

AKIMA, H., 1974b, “Algorithm 474, Bivariate Interpolation and Smooth Surface

Fitting Based on Local Procedures”, Communications of the ACM, v. 17,

n. 1, pp. 26–31.

AKIMA, H., 1978a, “A method of bivariate interpolation and smooth surface fitting

for irregularly Distributed Data Points”, ACM Transactions on Mathema-

tical Software, v. 4, n. 2, pp. 148–159.

AKIMA, H., 1978b, “Algorithm 525. Bivariate interpolation and smooth surface

fitting for irregularly Distributed Data Points”, ACM Transactions on

Mathematical Software, v. 4, n. 2, pp. 160–164.

AKIMA, H., 1991, “A method of univariate interpolation that has the accuracy of a

third-degree polynomial”, ACM Transactions on Mathematical Software,

v. 17, n. 3, pp. 341–366.

AKIMA, H., 1996, “Algorithm 760: Rectangular-Grid-Data Surface Fitting that

Has the Accuracy of a Bicubic Polynomial”, ACM Transactions on Mathe-

matical Software, v. 22, n. 3, pp. 357–361.

ALFORD, R. M., KELLY, K. R., BOORE, D. M., 1974, “Accuracy of finite diffe-

rence modeling of the acoustic wave equation”, Geophysics, v. 39, n. 6,

pp. 834–842.

ALKHALIFAH, T., FOMEL, S., 2001, “Implementing the fast marching eikonal

solver: spherical versus cartesian coordinates”, Geophysical, v. 49, pp. 165–

178.

128

BAYSAL, E., KOSLOFF, D., SHERWOOD, J. W. C., 1983, “Reverse Time Mi-

gration”, Geophysics, v. 48, pp. 1514–1525.

BAYSAL, E., KOSLOFF, D., SHERWOOD, J. W. C., 1984, “A two-way nonre-

flecting wave equation”, Geophysics, v. 49, pp. 132–14.

BERKHOUT, A. J., 1987, Applied Seismic Wave Theory. New York, Elsevier

Publishing Company.

BISHOP, C. M., 1996, Neural Networks for Pattern Recognition. Oxford, Oxford

University Press.

BLEISTEIN, N., 1984, Mathematical Methods for Wave Phenomena. New York,

Academic Press Inc.

BLEISTEIN, N., 1987, “On the imaging of reflectors in the earth”, Geophysical

Prospecting, v. 52, pp. 931–942.

BLEISTEIN, N., COHEN, J., JARAMILLO, H., 1999, “True-amplitude transfor-

mation to zero offset of data from curved reflectors”, Geophysics, v. 64,

pp. 112–129.

BOUE, M., DUPUIS, P., 1999, “Markov chain approximations for deterministic

control problems with affine dynamics and quadratic cost in the control”,

SIAM Journal on Scientific Computing, v. 36, n. 3, pp. 667–695.

BRUSS, A., 1982, “The eikonal equation: some results applicable to computer

vision”, Journal of Math. Phy., v. 23, n. 5, pp. 890–896.

BUHMANN, M. D., 2003, Radial Basis Functions: Theory and Implementations

(Cambridge Monographs on Applied and Computational Mathematics).

New York, Cambridge University Press.

BULCAO, A., 2004, Modelagem e Migracao Reversa no Tempo Empregando Opera-

dores Elasticos e Acusticos. Tese de Doutorado, UFRJ, Rio de Janeiro-RJ,

Outubro.

BULLEN, K. E., BRUCE, A. B., 1985, An Introduction to the Theory of Seismo-

logy. 4 ed. Ney York, Cambridge University Press.

CASTAGNA, J. P., BACKUS, M. M., 1993, Offset-dependent reflectivity - theory

and practice of AVO analysis, v. 8. Society of Exploration Geophysicists.

CERJAN, C., KOSLOFF, D., RESHEF, M., 1985, “A Nonreflecting Boundary

Condition for Discrete Acoustic and Elastic Wave Equations”, Geophysics,

v. 50, pp. 705–708.

129

CERVENY, V., 1985, “Ray Synthetic Seismograms for Complex Two- and Three-

Dimensional Structures”, J. Geophysics, v. 58, pp. 2–26.

CERVENY, V., 1987, Ray methods for three-dimensional seismic modeling.

Trondheim, Petroleum Industry Course.

CERVENY, V., 2001, Seismic ray theory. New York, Cambridge University Press.

CHANG, W., MCMECHAN, G., 1986, “Reverse-Time Migration of Offset Verti-

cal Seismic Profiling Data Using the Excitation Time Imaging Condion”,

Geophysics, v. 51, pp. 67–84.

CLAERBOUT, J., 1971,“Toward a unified theory of reflector imaging”, Geophysics,

v. 36, pp. 469–481.

CLAERBOUT, J., 1985, Imaging the earth’s interior. Inc, Blackwell Scientific

Publications.

COSTA, J. L., 2006, Migracao reversa no tempo de super-sismograma visando obje-

tivos localizados nos flancos e abaixo de domos salinos. Tese de Mestrado,

UFRJ, Rio de Janeiro-RJ, Dezembro.

CUNHA, P. E. M., 1997, Estrategias eficientes para migracao reversa no tempo

pre-empilhamento 3-D em profundidade pelo metodo das diferencas finitas.

Tese de Mestrado, UFBA, Salvador-Bahia, Dezembro.

DENG, F., MCMECHAN, G., 2007, “True Amplitude prestack depth migration”,

Geophysics, v. 72, pp. 155–166.

DORS, C., 2007, Propagacao de Ondas Elasticas Utilizando Funcoes de Green

Numericas Locais em Modelos Discretizados por Elementos Finitos. Tese

de Doutorado, UFRJ, Rio de Janeiro-RJ, Abril.

ENGQUIST, B., OSHER, S., 1980, “Stable and entropy satisfying approximati-

ons for transonic flow calculations”, Mathematics of Computation, v. 34,

pp. 45–75.

FARIA, E. L., 1986, Migracao antes do empilhamento utilizando propagacao reversa

no tempo. Tese de Mestrado, UFBA, Salvador-Bahia, Abril.

FASSHAUER, G. E., 2002, “Newton iteration with multiquadrics for the solution

of non-linear pdes”, Computers and Mathematics with Applications, v. 43,

n. 3-5, pp. 423–438.

130

FORTUNA, A. O., 2000, Tecnicas Computacionais para Dinamica dos Fluidos:

Conceitos Basicos e Aplicacoes. Sao Paulo, Editora Edusp.

FRANKE, R., 1982, “Scattered data interpolation tests of some methods”, Mathe-

matics of Computation, v. 38, n. 157, pp. 181–200.

GODUNOV, S., 1959, “Finite difference method for numerical computation of dis-

continuous solutions of the equation of fluid dynamics”, Matematicheskii

Sbornik, v. 47, n. 3, pp. 271–306.

GUITTON, A., VALENCIANO, A., BEVC, D., et al., 2007, “Smoothing image

condition for shot-profile migration”, Geophysics, v. 72, n. 3, pp. 149–154.

HARDY, R. L., 1971, “Multiquadric equations of topography and other irregular

surfaces”, Journal of Geophysical Research, v. 176, pp. 1905–1915.

HASSOUNA, M. S., FARAG, A. A., 2007, “Multistencils Fast Marching Methods:

A Highly Accurate Solution to the Eikonal Equation on Cartesian Do-

mains”, IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MA-

CHINE INTELLIGENCE, v. 29, n. 9, pp. 1–11.

HECHT, E., 2002, Optics. 4 ed. San Francisco, Addison Wesley.

JEFFREYS, H., JEFFREYS, B., 1950, Methods of Mathematical Physics (Cam-

bridge Mathematical Library). 2 ed. Cambridge, Cambridge University

Press.

JEONG, W. W., WHITAKER, R., 2008, “A Fast Iterative Method for Eikonal

Equations”, Society for Industrial and Applied Mathematics, v. 30, n. 5,

pp. 2512–2534.

JOHNSON, C., 1995, Numerical Solution of PDEs by the Finite Element Method.

Cambridge University Press.

KANSA, E. J., 1990, “Multiquadrics- A scattered data approximation scheme with

applications to computational fluid-dynamics: I. Surface approximations

and partial derivative estimates”, Computers and Mathematics with Ap-

plications, v. 19, n. 6-8, pp. 127–145.

KIM, S., 2001, “An O(N) level set method for eikonal equation”, SIAM Journal on

Scientific Computing, v. 22, n. 6, pp. 2193–2001.

LANCASTER, P., SALKAUSKAS, K., 1981, “Surfaces generated by moving least

squares methods”, Mathematical of Computational, v. 37, n. 155, pp. 141–

158.

131

LARSON, R. E., HOSTETLER, R. P., EDWARDS, B. H., 1998, Calculo com

Geometria Analıtica, v. 2. 5 ed. Rio de Janeiro, LTC.

LIMA, C. M. M., 2009, Um estudo sobre o metodo mınimos quadrados moveis por

aproximacoes iteradas. Tese de Mestrado, PUC-Rio, Rio de Janeiro-RJ,

Setembro.

LIU, G. R., 2002, Mesh Free Methods: Moving Beyond the Finite Element Method.

New York, CRC Press.

LOEWENTHAL, D., HU, L., 1991, “Two Methods for computing the imaging

condition for common-shot prestak migration”, Geophysics, v. 56, pp. 380–

399.

LOEWENTHAL, D., STOFFA, P. L., FARIA, E. L., 1987, “Suppressing the

Unwanted Refletions of the Full Wave Equation”, Geophysics, pp. 1007–

1012.

MALLADI, R., SETHIAN, J., 1996, “A unified approach to noise removal, image

enhancement, and shape recovery”, IEEE Trans. on Image Processing,

v. 5, n. 11, pp. 1554–1568.

MARGRAVE, G. F., 2003, Numerical Methods of Exploration Seismology with al-

gorithms in Matlab. Calgary-Canada, The University of Calgary.

MCMECHAN, G. A., 1983, “Migration by extrapolation by time dependent boun-

dary values”, Geophysical Prospecting, v. 31, pp. 413–420.

MICCHELLI, C. A., 1986, “Interpolation of Scatter Data: Distance Matrices and

Conditionally Positive Definite Functions”, Constructive Approximation,

v. 2, pp. 11–22.

MUFTI, I. R., 1990, “Large-scale three-dimensional seismic models and their inter-

pretive significance”, Geophysics, v. 55, pp. 1166–1182.

PAUL, M. K., 2004, “Rock properties and amplitude versus angle”, The Leading

Edge, v. 23, n. 10, pp. 986–988.

PESTANA, M., PIMENTEL, A., 1996, “Calculo dos tempos de transito por dife-

rencas finitas: uma comparacao”, Revista Brasileira de Geofısica, v. 14,

pp. 115–130.

PIMENTEL, A. L., 1994, Tracado de raios usando o princıpio da reciplocidade.

Tese de Mestrado, UFBA, Salvador-Bahia, Agosto.

132

PODVIN, P., LECOMTE, I., 1991, “Finite-difference computation of traveltimes

in very contrasted velocity models: A massively parallel approach and its

associated tools”, Geophys. J. Int., v. 105, pp. 271–284.

POPOV, M. M., 2002, Ray theory and gaussian beam method for geophysicists.

SALVADOR-BAHIA, EDUFBA.

POPOVICI, A. M., 1991, Finite-difference traveltime maps. Relatorio Tecnico 70,

Stanford Exploration.

RAWLINSON, N., SAMBRIDGE, M., 2005, “The fast marching method: an effec-

tive tool for tomographic imaging and tracking multiple phases in complex

layered media”, Exploration Geophysics, v. 36, pp. 341–350.

RAWLINSON, P. E., LIN, Q., 2003, “A Modified Fast Marching Method”, Proc.

Scandinavian Conf. Image Analysis, pp. 1154–1161.

RESHEF, P., KOSLOFF, D., 1986, “Migration of Common Shot Gathers”, Ge-

ophysics, v. 51, pp. 324–331.

RESNICK, J., PATRICK, N., LARNER, K., 1987, “Amplitude versus offset analy-

sis in the presence of dip”, 57th Ann. Internat. Mtg. - SEG Expanded

Abstracts, pp. 617–620.

REYNOLDS, A. C., 1978, “Boundary conditions for the numerical solution of wave

propagation problems”, Geophysics, v. 43, pp. 1099–1110.

RIPPA, S., 1999, “An Algorithm for selecting a good value for the parameter c in

radial basis function interpolation”, Adv. Comput. Math, v. 11, pp. 193–

210.

SANTIAGO, T. M. G., 2004, Migracao Kirchhoff 2,5D em Tempo no Domınio

de Angulo Comum e em Amplitude Verdadeira. Tese de Mestrado, UNI-

CAMP, Campinas-Sao Paulo, Dezembro.

SCHABACK, R., 2000, “Remarks on Meshless Local Construction of Surfaces”,

Proceedings of the 9th IMA Conference on the Mathematics of Surfaces,

pp. 34–58.

SCHNEIDER, W. A., 1995, “Robust and efficient upwind finite difference travel-

time calculations in 3d”, Geophysics, v. 60, pp. 1108–1117.

SCHNEIDER, W. A. A., RANZINGER, K. A., BALCH, A. H., et al., 1992, “A

dynamic programming approach to first traveltime computation in media

with arbitrarily distributed velocities”, Geophysics, v. 57, pp. 39–50.

133

SCHNEIDER, W. S., 1978, “Integral formulation for migration in two and three

dimensions”, Geophysics, v. 43, n. 1, pp. 49–76.

SETHIAN, J. A., 1996, “A fast marching level set method for monotonically advan-

cing fronts”, Proceedings of National Academy of Sciences, v. 93, pp. 1591–

1595.

SETHIAN, J. A., 1999a, Level Sets Methods and Fast Marching Methods. 2 ed.

New York, Cambridge University Press.

SETHIAN, J. A., 1999b, “Fast marching methods”, SIAM Review, v. 41, n. 2,

pp. 199–235.

SETHIAN, J. A., POPOVICI, A. M., 1999,“Three-dimensional Traveltime Compu-

tation Using the Fast Marching Method”, Geophysics, v. 22, pp. 516–523.

SHERIFF, R. E., 2002, Encyclopedic Dictionary of Applied Geophysics (Geophysi-

cal References, V. 13). USA, Society of Exploration.

SILVA, J. J., 2009, Migracao Reversa no Tempo na Determinacao das Amplitudes

de Reflexao em Funcao do Angulo de Incidencia. Tese de Doutorado,

UFRJ, Rio de Janeiro-RJ, Marco.

SILVA, M. W. X., 2008, Estudo da Variacao de Parametros de Aquisicao de Da-

dos Sısmicos Associado ao Imageamento de Falhas Utilizando Migracao

Reversa no Tempo. Tese de Mestrado, UFRJ, Rio de Janeiro-RJ, Marco.

STRIKWERDA, J., 2004, Finite Difference Schemes and Partial Differential Equa-

tions. 2 ed. Philadelphia, SIAM: Society for Industrial and Applied Mathe-

matics.

TSITSIKLIS, J., 1995, “Efficient Algorithms for Globally Optimal Trajectories”,

IEEE Trans. Automatic Control, v. 40, n. 9, pp. 1528–1538.

TYGEL, M., SCHLEICHER, J., HUBRAL, P., et al., 1993, “Multiple weights in

diffraction stack migration”, Geophysics, v. 58, n. 12, pp. 1820–1830.

TYGEL, M., SANTOS, L. T., JCHLEICHER, J., et al., 1999, “Kirchhoff imaging

as a tool for AVO/AVA analysis”, The Leading Edge, v. 18, pp. 940–945.

VAN TRIER, J., SYMES, W., 1991, “Upwind finite-difference calculation of tra-

veltimes”, Geophysics, v. 56, n. 6, pp. 812–821.

VASQUEZ, A. C. R., 1999, Recuperacao de Atributos Sısmicos Utilizando a Migra-

cao para Afastamento Nulo. Tese de Mestrado, UNICAMP, Campinas-Sao

Paulo, Dezembro.

134

VASQUEZ, A. C. R., OLIVEIRA, A. S., TYGEL, M., 2003, “Recuperacao de Atri-

butos Sısmicos Utilizando Migracao para Afastamento Nulo”, Brasilian

Journal of Geophysics, v. 20, pp. 59–65.

VAZ, J. R. P., 2004, Solucao da Equacao Iconal por Diferencas Finitas em 3D em

Meios Anisotropicos. Tese de Mestrado, UFPA, Belem-Para, Outubro.

VEEKEN, P., RAUCH-DAVIES, M., 2006, “AVO attribute analysis and seismic

reservoir characterization”, First Break, v. 24, pp. 41–52.

VEEKEN, P. C. H., 2007, Seismic Stratigraphy, Basin Analysis and Reservoir

Characterisation, v. 37, Handbook of Geophysical Exploration: Seismic

Exploration. Amsterdam, Elsevier Science.

VIDALE, J., 1988, “Finite-difference calculation of traveltimes”, Bull. Seis. Soc.

Am., v. 78, n. 6, pp. 2062–2076.

VIDALE, J., 1990,“Finite-difference calculation of traveltimes in three dimensions”,

Geophysics, v. 55, n. 5, pp. 521–526.

WHITMORE, N. D., 1983, “Iterative depth migration by backward time propaga-

tion”, 53rd Ann. Internat. Mtg. SEG - Expanded Abstracts, pp. 225–244.

WIGGINS, J. W., 1984,“Kirchhoff integral extrapolation an migration of nonplanar

data”, Geophysics, v. 49, pp. 1239–1248.

YATZIV, L., BARTESAGHI, A., SAPIRO, G., 2006, “O(N) Implementation of the

Fast Marching Algorithm”, J. Computational Physics, v. 212, pp. 393–399.

YILMAZ, O., 2001, Seismic Data Analysis: Processing, Inversion and Interpreta-

tion of Seismic Data. Tulsa, Society of Exploration Geophysicits.

ZAUDERER, R., 1989, Partial Differential Equations of Applied Mathematics.

John Wiley and Sons.

ZHANG, L., 1991, “Finite difference calculation of Green’s function”, SEP, v. 72,

pp. 153–170.

ZHANG, M., FU, L., LI, X. L. X., 2009, “2D efficient ray tracing with a modified

shortest path method”, Exploration Geophysics, v. 40, pp. 301–307.

ZHANG, Y., SUN, J., 2008, “Practical issues of reverse time migration: true-

amplitude gathers, noise removal and harmonic-source encoding”, 70th

EAGE Conference & Exhibition.

135

ZHAO, H., 2004, “A Fast Sweeping Method for Eikonal Equations”, Mathematics

of Computation, v. 74, n. 250, pp. 603–627.

ZHAO, H. K., OSHER, S., MERRIMAN, B., et al., 2000, “Implicit and non-

parametric shape reconstruction from unorganized points using variational

level set method”, Computer Vision and Image Understanding, v. 80, n. 3,

pp. 295–319.

136

Apendice A

Operadores de Diferencas Finitas

Neste apendice, sera descrito a forma de obtencao dos operadores de diferencas

finitas, utilizados ao longo deste trabalho, atraves da serie de Taylor.

O metodo das diferencas finitas pode ser utilizado para resolver problemas de

valor de contorno ou valor inicial, envolvendo equacoes diferenciais ordinarias ou

parciais (FORTUNA, 2000; STRIKWERDA, 2004). A tecnica consiste em substituir

cada derivada ou diferencial das equacoes diferenciais por aproximacoes de diferencas

finitas ou acrescimos finitos das variaveis.

A.1 Expansao em Serie de Taylor

Como os operadores de diferencas finitas tem como base a expansao em serie de

Taylor de uma funcao contınua, supoem-se que tal funcao seja ϕ, definida em um

intervalo [a, b] de interesse e que possua derivadas ate ordem ~ contınuas no intervalo.

Assim, o Teorema de Taylor nos permite escrever, para todo ponto x ∈ [a, b]:

ϕ(x ±h∆x) =

∞∑~=0

(±h ∆x)~ϕ~(x)~!

. (A.1)

Considerando h = 1, tem-se:

ϕ(x−∆x) = ϕ(x)−∆xdϕ(x)

dx+

(∆x)2

2!d2ϕ(x)

dx2 −(∆x)3

3!d3ϕ(x)

dx3 +(∆x)4

4!d4ϕ(x)

dx4 −O(∆x)~, (A.2)

ϕ(x+∆x) = ϕ(x)+∆xdϕ(x)

dx+

(∆x)2

2!d2ϕ(x)

dx2 +(∆x)3

3!d3ϕ(x)

dx3 +(∆x)4

4!d4ϕ(x)

dx4 +O(∆x)~, (A.3)

137

por outro lado, considerando h = 2:

ϕ(x−2∆x) = ϕ(x)−2∆xdϕ(x)

dx+

(2∆x)2

2!d2ϕ(x)

dx2 −(2∆x)3

3!d3ϕ(x)

dx3 +(2∆x)4

4!d4ϕ(x)

dx4 −O(∆x)~,(A.4)

ϕ(x+2∆x) = ϕ(x)+2∆xdϕ(x)

dx+

(2∆x)2

2!d2ϕ(x)

dx2 +(2∆x)3

3!d3ϕ(x)

dx3 +(2∆x)4

4!d4ϕ(x)

dx4 +O(∆x)~,(A.5)

onde O(∆x)~ representa a ordem do erro local de truncamento.

A.2 Operadores para a Primeira Derivada

Explicitando-se a primeira derivada de ϕ e omitindo o termo relacionado a ordem

do erro em (A.2) e (A.3), obtem-se o operador de primeira ordem dado por:

dϕ(x)dx

=ϕ(x) − ϕ(x − ∆x)

∆x(A.6)

e

dϕ(x)dx

=ϕ(x + ∆x) − ϕ(x)

∆x. (A.7)

sendo que (A.6) representa o operador de diferencas atrasadas (backward differences)

e (A.7) de diferencas progressivas (forward differences), ambos chamados de formula

de dois pontos.

Para obter os operadores de diferencas atrasadas e progressivas de segunda or-

dem, faz-se: (A.4)-4(A.2) e (A.5)-4(A.3), respectivamente. Apos explicitar o termo

relacionado a derivada primeira e omitir a ordem do erro, tem-se:

dϕ(x)dx

=3ϕ(x) − 4ϕ(x − ∆x) + ϕ(x − 2∆x)

∆x(A.8)

e

dϕ(x)dx

= −3ϕ(x) − 4ϕ(x + ∆x) + ϕ(x + 2∆x)

2∆x. (A.9)

Uma outra possibilidade para o operador de primeira derivada de segunda ordem

e combinar as expressoes (A.2) e (A.3) de forma a eliminar a segunda derivada de

ϕ. Para isso basta subtrair estas duas expressoes e explicitar a primeira derivada

omitindo o termo relacionado a ordem do erro. Daı, tem-se o classicamente chamado

138

operador de diferencas centrais (central differences):

dϕ(x)dx

=ϕ(x + ∆x) − ϕ(x − ∆x)

2∆x. (A.10)

Para obter um operador de quarta ordem para a derivada primeira, expande-se os

termos da serie nas expressoes (A.2) a (A.5) ate a derivada quarta para que a ordem

do erro seja O(∆x)5. Entao, faz-se, respectivamente, (A.2)-(A.3) e (A.4)-(A.5):

ϕ(x − ∆x) + ϕ(x + ∆x) = −2(∆x)dϕ(x)

dx−

2(∆x)3

6d3ϕ(x)

dx3 − O(∆x)5, (A.11)

ϕ(x − 2∆x) − ϕ(x + 2∆x) = −4(∆x)dϕ(x)

dx−

2(∆x)3

6d3ϕ(x)

dx3 − O(∆x)5. (A.12)

Daı, o passo final e fazer (A.12)-8(A.11) para eliminar o termo relacionado a terceira

derivada. E assim, omitindo a ordem do erro e explicitando a derivada primeira,

chega-se ao operador de quarta ordem:

dϕ(x)dx

=1

12(∆x)[−ϕ(x + 2∆x) + 8 (ϕ(x + ∆x) − ϕ(x − ∆x)) + ϕ(x − 2∆x)

]. (A.13)

A.3 Operadores para a Segunda Derivada

Utilizando as expansoes (A.2) e (A.3), pode-se combina-las para que a primeira

derivada de ϕ seja eliminada e assim obter um operador de segunda ordem para a

segunda derivada. Esse procedimento e feito somando-se (A.2) com (A.3) e, posteri-

ormente, explicitando a segunda derivada sem o termo relacionado a ordem do erro,

tem-se:

d2ϕ(x)dx2 =

ϕ(x + ∆x) − 2ϕ(x) + ϕ(x − ∆x)(∆x)2 . (A.14)

Por outro lado, se o objetivo e obter um operador de quarta ordem para a segunda

derivada, expande-se os termos da serie nas expressoes (A.2) a (A.5) ate a derivada

quinta, a fim de que a ordem do erro seja O(∆x)6. Daı, faz-se, respectivamente,

(A.2)+(A.3) e (A.4)+(A.5):

ϕ(x + ∆x) + ϕ(x − ∆x) = 2ϕ(x) + (∆x)2 d2ϕ(x)dx2 +

(∆x)4

12d4ϕ(x)

dx4 + O(∆x)6, (A.15)

139

ϕ(x + 2∆x) + ϕ(x − 2∆x) = 2ϕ(x) + 4(∆x)2 d2ϕ(x)dx2 +

16(∆x)4

12d4ϕ(x)

dx4 + O(∆x)6. (A.16)

Sendo assim, o ultimo procedimento e fazer 16(A.15) - (A.16) para eliminar o termo

de derivada quarta. Entao, suprimindo a ordem do erro e isolando a derivada se-

gunda, tem-se o operador de quarta ordem:

d2ϕ(x)dx2 =

112(∆x)2

[−ϕ(x + 2∆x) + 16 (ϕ(x + ∆x) + ϕ(x − ∆x)) − 30ϕ(x) − ϕ(x − 2∆x)

].

(A.17)

140

Documents

Uma metodologia para extração de ângulos de reflexão em …objdig.ufrj.br/60/teses/coppe_m/WilsonSouzaDuarte.pdf · 2011. 9. 26. · UMA METODOLOGIA PARA EXTRAC˘AO DE~ ANGULOS