INVERSAO TOMOGR˜ AFICA SEQUENCIAL PARA O CAMPO DE

INVERSAO TOMOGRAFICA SEQUENCIAL PARA O CAMPO DE VELOCIDADES

SISMICAS BASEADA EM DIFRACOES E CRITERIOS GEOLOGICOS

Luiz Alberto Santos

Tese de Doutorado 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 Doutor em Engenharia

Civil.

Orientadores: Webe Joao Mansur

George Albert McMechan

Rio de Janeiro

Junho de 2012



Luiz Alberto Santos

TESE 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 DOUTOR EM

CIENCIAS EM ENGENHARIA CIVIL.

Examinada por:

Prof. Webe Joao Mansur, Ph.D.

Prof. George Albert McMechan, Ph.D.

Prof. Otto Correa Rotunno Filho, Ph.D.

Profa. Valeria Cristina Ferreira Barbosa, D.Sc.

Prof. Jesse Carvalho Costa, D.Sc.

Dr. Djalma Manoel Soares Filho, D.Sc.

Dr. Gerson Luis da Silva Ritter, D.Sc.

Dr. Andre Bulcao, D.Sc.

RIO DE JANEIRO, RJ – BRASIL

JUNHO DE 2012

Santos, Luiz Alberto

Inversao tomografica sequencial para o campo de

velocidades sısmicas baseada em difracoes e criterios

geologicos/Luiz Alberto Santos. – Rio de Janeiro:

UFRJ/COPPE, 2012.

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



Tese (doutorado) – UFRJ/COPPE/Programa de Engenharia

Civil, 2012.

Referencias Bibliograficas: p. 119 – 128.

1. Problema Inverso. 2. Operadores focais (CFP). 3.

Tomografia. 4. Difracoes. I. Mansur, Webe Joao et al.. II.

Universidade Federal do Rio de Janeiro, COPPE, Programa de

Engenharia Civil. III. Tıtulo.

iii

Aos meus queridos Luizes: Carlos,

onde tudo comecou; Henrique e

Eduardo, mola mestra e energia

viva do gerundio e do futuro. A

minha amada Claudia, parceira na

edificacao de nossas vidas. A irma

Patrıcia pela visao rebelde e

vanguardista, inspiracao do porvir.

A minha mae Ercine Terezinha,

somente por tudo.

iv

Agradecimentos

Os sabios dizem que, na vida, ha tres fases do aprendizado: aprendiz, praticante e

mestre. Na primeira o indivıduo e um aluno e somente le e estuda. Na segunda, ele aplica

o que aprendeu e expande seus conhecimentos com a pratica. Na terceira, o indivıduo

passa aos outros aquilo que sabe. Em mim, prevalece o aprendiz que nao avancaria sem a

valiosa ajuda dos diversos colegas, mestres e familiares. Por isso gostaria de agradecer as

pessoas a seguir.

Primeiramente a minha mae, Ercine Terezinha Santos, minha primeira, e eterna mes-

tra. Nao consigo pensar nela sem me emocionar. Obrigado, minha mae.

A minha esposa, Claudia e aos nossos filhos, os Luizes - Henrique e Eduardo - pela

alegria e motivacao e, acima de tudo, pelos seus sacrifıcios, criados pelas ausencias deste

marido e pai em funcao do trabalho e dos estudos. Meu muito obrigado a voces, e tambem

minhas desculpas.

Ao sempre amigo, irmao de caminhada, geologo, poeta e empresario Jason Carneiro,

pela amizade, pelas correcoes gramaticais e dicas no texto da tese. Valeu, Bro !!

Ao professor e orientador Webe Joao Mansur pela orientacao, pelos conselhos, por

proporcionar um ambiente prolıfico para criacao e por transmitir tranquilidade nos mo-

mentos mais duros da pesquisa.

Ao professor e orientador George A. McMechan pela orientacao e bons exemplos.

Constatei com ele a maxima aristotelica que o exemplo, de fato, e a principal ferramenta

de educacao e ensino.

A Carlos Eduardo Theodoro, colega no Centro de Pesquisas da Petrobras (CENPES),

pelas discussoes, orientacoes, bons exemplos e transmissao de suas ricas experiencias no

universo da geofısica.

Ao colega de pesquisa na Petrobras, Djalma Manoel Soares Filho, grande mestre e

guru que me iniciou no mundo da modelagem sısmica e na tomografia. O gerundio e

v

nosso tempo preferido. Valeu, Mestre.

Ao tambem colega de pesquisa no CENPES, Andre Bulcao, pelas instrucoes sobre

migracao RTM, ondulatoria, programacao e amizade. Saiba que meu primeiro codigo de

RTM acustico saiu depois de um bate papo regado a acaraje la em Salvador. Obrigado,

brother.

Nao poderia deixar de citar o colega Paulo Eduardo Miranda Cunha pelo auxılio nas

deducoes mais complicadas, pelos bate-papos abordando filosofia, religiao e cultura e

ainda pelas primeiras aulas sobre Seismic Unix. O primeiro comando desta util ferra-

menta, a2b (conversao de arquivos ascii para binarios), me foi ensinada por ti. Valeu,

meu irmao.

Ao colega Ricardo Braganca pelas aulas de scripts e dicas do ambiente UNIX e Linux.

Ainda no universo Linux, agradeco ao colega Ricardo F. Chartuni Cabral da Cruz pelas

instalacoes que muito me ajudaram na conclusao da tese.

Ao colega e professor Eduardo Filpo responsavel pelos meus primeiros passos na se-

ara de migracao Kirchhoff ainda em 2003. Ainda tenho as anotacoes do curso introdutorio

de geofısica da Petrobras.

Aos colegas e professores Andre Romanelli e Carlos Cunha Filho pelas sugestoes e

discussoes.

Ao geofısico e gerente do CENPES/PDEXP/GEOF ate dezembro de 2010, Dr. Edu-

ardo Lopes de Faria por proporcionar as condicoes para pesquisa, pelo apoio na liberacao

para o doutorado e, por que nao, pelos questionamentos. Eles me fizeram amadurecer

muitas das ideias.

Ao atual gerente do CENPES/PDGEO/GEOF, Rui Pinheiro Silva, pelas otimas

condicoes proporcionadas na conclusao deste trabalho.

Ao colega de doutorado Rodrigo Dias, pelas discussoes, trocas bibliograficas e auxılio

nas deducoes mais complicadas.

Ao prestativo Raul Flores e a Franciane Peters, pelas dicas de programacao no

domınio da frequencia e pela ajuda nos preparativos da apresentacao da defesa.

Ao professor Eduardo Gomes Dutra do Carmo, pelas fundamentais aulas de problema

inverso e pelos bate-papos sobre religiao e filantropia. Saiba que meu interesse por esses

temas menos materiais renasceram com nossas conversas. Receba meu muito obrigado,

Dudu.

vi

Ao Dr. Cleberson Dors, pesquisador da COPPE, pelas aulas de corredor, muitıssimo

uteis ao longo do desenvolvimento do trabalho.

A Ivone, nossa administradora do LAMEC (Laboratorio de Mecanica Computacio-

nal), sem a qual nada andaria. Muito obrigado, Ivone.

Aos colegas do Departamento de Geociencias da Universidade do Texas em Dallas,

em especial a Hu Jin, Ernesto Oropeza, Roberto Falcon, Sylvia Pacheco, Xinfa Zhu,

pelas discussoes tecnicas e a Sharon Edwards, sempre facilitando a vida nos assuntos

administrativos.

Enfim, estendo os meus agradecimentos a instituicao Petrobras, representada pelos

seus diretores, alta gerencia e demais colegas de trabalho. Um corpo tecnico de altıssimo

nıvel que constitui uma fonte de inspiracao para mim.

Aos colegas da pos-graduacao na Engenharia Civil da COPPE-UFRJ e aos professores

pelo ambiente agradavel e de cortesia que impera em todas as salas e corredores.

A todos voces, citados explıcita ou implicitamente, meus sinceros agradecimentos. E,

por favor, sintam-se isentos de qualquer incorrecao presente neste trabalho. Sao arestas

ainda nao desbastadas neste humilde aprendiz.

vii

Resumo da Tese apresentada a COPPE/UFRJ como parte dos requisitos necessarios para

a obtencao do grau de Doutor em Ciencias (D.Sc.)



Luiz Alberto Santos

Junho/2012



Programa: Engenharia Civil

Esta tese avanca no emprego dos operadores focais (Common Focus Point Opera-

tor (CFPO)) para inversao tomografica de velocidades e estruturas de subsuperfıcie. Os

CFPOs constituem um neologismo para a funcao de Green de uma fonte deflagrada em

profundidade e registrada na superfıcie. O processo de inversao divide-se em duas partes:

a obtencao dos CFPOs e a tomografia sobre estes. A estimativa dos CFPOs de forma

convencional e revisitada e dois outros esquemas sao propostos nesta tese: por varredura

e atraves de difracoes. Utilizando-se somente operadores focais de difracao (Diffraction

Common Focus Operator (DCFPO)) demonstra-se que e possıvel estimar o campo de ve-

locidades a partir de dados sısmicos monocanal. Devido as caracterısticas geometricas do

dispositivo de inversao dos CFPOs ou DCFPOs, verifica-se que as formas mais estaveis

matematicamente para parametrizacao do modelo fornecem bons resultados estruturais,

mas o campo de velocidades e pouco realıstico. A consistencia do modelo invertido e

quantificada atraves do Indice de Incoerencia Geologica (IIG), proposto nesta tese, e que

relaciona campo de velocidades e estruturas. Estudos realizados com as inversoes glo-

bal, onde todo o modelo e invertido simultaneamente, e camada a camada, demonstram

que a inversao sequencial em pelo menos 2 passos proporciona modelos geologicamente

mais realısticos, com baixo IIG. Nomeou-se este processo como tomografia sequencial.

Aplicacoes em exemplos sinteticos e reais comprovam a viabilidade das teses propostas.

viii

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

for the degree of Doctor of Science (D.Sc.)

SEQUENTIAL TOMOGRAPHIC INVERSION FOR SEISMIC VELOCITY FIELD

BASED ON DIFFRACTIONS AND GEOLOGIC CRITERIA

Luiz Alberto Santos

June/2012

Advisors: Webe Joao Mansur


Department: Civil Engineering

This thesis has advances in the use of focusing operators (Common Focus Point Oper-

ator (CFPO)) for tomographic inversion for velocity and subsurface structures estimation.

The CFPO is a neologism for a Green function of a point source at depth and recorded

at the surface. The inversion process is divided in two steps: the CFPO estimation and

the tomography over CFPO. The conventional CFPO estimation is revisited and two other

schemes for CFPO estimation are proposed: by scanning and using diffractions. By us-

ing only diffractions CFPO (DCFPO), it is shown that the velocity field can be obtained

with single channel seismic data. Because of the geometric characteristics of the CFPO

or DCFPO device, the most mathematically stable model parameterization delivers good

structural results, but the velocity field is not geologically realistic. The consistency of the

inverted model is quantified through the Geological Incoherence Index (GII), that relates

velocity field and structures. By studying global inversion, in which the model is entirely

inverted simultaneously, and layer-stripping one, it is shown that the sequencial inversion

in at least 2 steps delivers more geologically realistic models, with low GII. This inver-

sion process is called sequential tomogaphy. Synthetic and field examples comprove the

proposed theses.

ix

Sumario

Lista de Figuras xiv

1 Introducao Geral 1

1.1 Estrutura da tese e suas contribuicoes cientıficas . . . . . . . . . . . . . . 2

1.2 Fomulacao do problema e objetivos . . . . . . . . . . . . . . . . . . . . 3

1.3 Historico de pesquisas voltados para inversao de dados sısmicos . . . . . 5

1.3.1 Inversao cinematica de dados sısmicos . . . . . . . . . . . . . . . 5

2 O modelo (m) 11

2.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.2 Comportamento das rochas . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3 Coerencia geologica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.4 Discretizacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.4.1 Discretizacao por blocos ou celulas em malha regular . . . . . . . 19

2.4.2 Discretizacao por blocos ou celulas em malha nao regular . . . . 19

2.4.3 Discretizacao por polinomios . . . . . . . . . . . . . . . . . . . 20

2.4.4 Discretizacao por funcao exponencial ou logarıtmica . . . . . . . 20

2.4.5 Discretizacao por funcao de base radial gaussiana . . . . . . . . . 21

2.4.6 Discretizacao por funcoes harmonicas . . . . . . . . . . . . . . . 21

2.4.7 Discretizacao por splines cubicas (1D e 2D) . . . . . . . . . . . . 22

2.4.8 Discretizacao camada a camada . . . . . . . . . . . . . . . . . . 23

2.5 Formas de calculo do erro . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.6 Comparacao entre as formas de discretizacao . . . . . . . . . . . . . . . 24

2.7 Discussao e propostas para discretizacao geologicamente coerente . . . . 25

2.8 Conclusoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

x

3 Modelagem (D = Lm) 34

3.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.2 Resolucao da equacao da onda por metodos integrais . . . . . . . . . . . 34

3.3 Resolucao da equacao da onda por metodos diferenciais . . . . . . . . . . 37

3.4 Aproximacao da equacao da onda por tracado de raios . . . . . . . . . . 38

3.5 Modelos cinematicamente equivalentes . . . . . . . . . . . . . . . . . . 39

3.6 Comentarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

4 Inversao (m=L−1 D) 42

4.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4.2 Metodos de resolucao do problema inverso . . . . . . . . . . . . . . . . . 44

4.3 Resolucao do problema inverso pelos metodos de gradientes . . . . . . . 45

4.3.1 Metodo da maxima declividade . . . . . . . . . . . . . . . . . . 45

4.3.2 Metodo de Newton . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.3.3 Metodo de Gauss-Newton . . . . . . . . . . . . . . . . . . . . . 47

4.3.4 Metodo de Levenberg-Marquardt . . . . . . . . . . . . . . . . . 49

4.3.5 Pseudo-Inversa de Moore-Penrose . . . . . . . . . . . . . . . . . 49

4.3.6 Metodo dos gradientes conjugados . . . . . . . . . . . . . . . . . 51

4.4 Influencia do dispositivo de aquisicao . . . . . . . . . . . . . . . . . . . 52

4.5 Transformando um problema mal posto em bem posto . . . . . . . . . . 56

4.5.1 Parametrizacao do modelo . . . . . . . . . . . . . . . . . . . . . 56

4.5.2 Regularizacao atuando na concepcao da funcao objetivo . . . . . 57

4.5.3 Regularizacao atuando sobre a matriz de sensibilidade - Tikhonov

de ordem zero . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.5.3.1 Aplicacao de restricoes na matriz de sensibilidade e

controle da norma da perturbacao no modelo . . . . . . 61

5 Obtendo operadores focais 65

5.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.2 Estimativa de operadores focais pela abordagem convencional . . . . . . 65

5.3 Estimativa de operadores focais por varredura . . . . . . . . . . . . . . . 69

5.4 Estimativa de operadores focais a partir de difracoes . . . . . . . . . . . . 75

5.5 Fatores a serem considerados na tomografia sobre DCFPO . . . . . . . . 76

xi

5.5.1 O efeito do afastamento na estimativa da velocidade . . . . . . . 76

5.5.2 Eventos que devem ser empregados na tomografia DCFPO . . . . 79

6 Tomografia sobre operadores focais 86

6.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

6.2 Derivada em relacao a posicao do ponto focal . . . . . . . . . . . . . . . 87

6.3 Derivada em relacao ao campo de velocidades . . . . . . . . . . . . . . . 87

6.3.1 Derivada em relacao a v sob discretizacao por blocos . . . . . . . 88

6.3.2 Derivada em relacao a v sob discretizacao por polinomios . . . . 88

6.3.3 Derivada em relacao a v sob discretizacao por funcoes exponen-

cial e logarıtmica . . . . . . . . . . . . . . . . . . . . . . . . . . 89

6.3.4 Derivada em relacao a v sob discretizacao por funcao de base radial 89

6.3.5 Derivada em relacao a v sob discretizacao por funcoes harmonicas 90

6.3.6 Derivada em relacao a v sob discretizacao por splines cubicas (1D

e 2D) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

6.3.7 Derivada em relacao a v sob discretizacao camada a camada . . . 91

6.4 Consideracoes sobre a abertura dos operadores focais . . . . . . . . . . . 92

6.5 Consideracoes sobre a abordagem do problema inverso - inversao global

ou camada a camada . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

6.6 Procedimento para estimativa do melhor MCE com baixo IIG: tomografia

sequencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

6.7 Aplicacoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

6.7.1 Inversao tomografica baseada em operadores CFP e DCFP - da-

dos sinteticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

6.7.2 Inversao tomografica baseada em operadores DCFP - dados

sısmicos monocanal . . . . . . . . . . . . . . . . . . . . . . . . 108

6.8 Discussao sobre as aplicacoes . . . . . . . . . . . . . . . . . . . . . . . 110

7 Discussao e conclusoes 116

7.1 Discussao final . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

7.2 Conclusoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

Referencias Bibliograficas 119

xii

A Experimento sismo 1 129

B Experimento discretizacao 134

B.1 Funcao polinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

B.2 Funcao gaussiana 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

B.3 Funcoes harmonicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

B.4 Funcao spline cubica 2D . . . . . . . . . . . . . . . . . . . . . . . . . . 138

B.5 Resumo dos resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

C Obtendo operadores CFP pela abordagem convencional 142

D Posicao x do tempo mınimo e velocidade aparente a partir de difracoes regis-

tradas em secoes de afastamento nao nulo 148

xiii

Lista de Figuras

1.1 Descricao do processo de analise de velocidade para migracao (MVA),

extraıdo de Stork (1992). . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1 Curva de variacao da velocidade P em funcao da pressao de um arenito

extraıda de Carmichael (1982). . . . . . . . . . . . . . . . . . . . . . . . 13

2.2 Relacao de estruturas basicas e propriedades (cores) observadas em ro-

chas sedimentares. (A) Acamamento sedimentar plano-paralelo; (B)

Acamamento sedimentar plano-paralelo com variacao lateral de facies;

(C) Dobras; (D) Falhamento normal; (E) Falhamento reverso; (F) In-

trusao discordante. As estruturas mais complexas sao obtidas atraves da

combinacao de uma ou mais situacoes observadas acima. . . . . . . . . . 14

2.3 (a) Registro de velocidade de uma pilha de rochas em uma bacia sedimen-

tar brasileira. (b) Versao em baixa frequencia espacial da curva (a). . . . . 16

2.4 (a) Camada com velocidade variando verticalmente limitada no topo e na

base por superfıcies horizontais (linhas tracejadas negra e branca respec-

tivamente no topo e na base). (b) θv calculado com a equacao 2.2 para (a)

gerando o angulo de 90o para todo o modelo, indicando o gradiente verti-

cal de velocidade. (c) θs (veja o texto para maiores detalhes) para o mo-

delo em (a). A normal aponta para a vertical (90o). (d) Camada com ve-

locidade variando lateralmente limitada no topo e na base por superfıcies

horizontais (linhas tracejadas brancas).(e) θv calculado com equacao 2.2

para (d) gerando o angulo de 0o para todo o modelo, o que indica gradi-

ente lateral de velocidade. (f) θs (veja texto para maiores detalhes) para o

modelo em (d). A normal aponta para direcao vertical (90o). . . . . . . . 18

xiv

2.5 Campo de velocidade obtido com a equacao 2.18 de Kabir, para v0=2000

e B=-1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.6 Campo de velocidade obtido com a equacao 2.22 para v0=2000 e B=-1.0. 29



3.1 Representacao de um campo de ondas propagando-se para cima e re-

gistrado na superfıcie S 0. Deseja-se estimar o campo de ondas em A

conhecendo-se o campo registrado na superfıcie S 0. . . . . . . . . . . . 36

3.2 Representacao de um modelo de duas celulas com velocidades v1 e v2.

A estrela representa a fonte de um raio, que viaja desta ate o receptor, o

triangulo. A dimensao da malha quadrada e DX. . . . . . . . . . . . . . . 40

4.1 Funcao objetivo com duas variaveis m1 e m2 ilustrando as iteracoes pelo

metodo de maxima declividade ate a convergencia para solucao em m

(modificado de Shewchuck, 1994). . . . . . . . . . . . . . . . . . . . . . 46

4.2 (a) Modelo de velocidade discretizado por celulas quadradas; (b)

Aquisicao sısmica hipotetica avaliando cada celula com um par fonte-

receptor, representados respectivamente pela estrela estrela vermelha e

o triangulo preto. As flechas pontilhadas indicam o caminho percorrido

pelo raio da fonte para o receptor. . . . . . . . . . . . . . . . . . . . . . 53

4.3 Diagrama kx vs kz para um dispositivo de aquisicao do tipo modelo explo-

sivo. Todos os numeros de onda sao amostrados igualmente. . . . . . . . 54

4.4 Diagrama kx vs kz para um dispositivo de aquisicao do tipo split-

spread. Observam-se heterogeneidades e vazios de amostragem de alguns

numeros de onda. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.5 Funcao objetivo com duas variaveis m1 e m2 dada pela equacao 4.42.

Ilustra-se a convergencia para o mınimo a partir de dois modelos inici-

ais (m1o ,m2o) = (2.0,2.0) e (m1o ,m2o) = (-2.0,-2.0) gerando as respectivas

solucoes (0.0,2.0) e (0.0,-2.0). Neste caso ha solucoes ambıguas. . . . . . 59

xv

4.6 Funcao objetivo com duas variaveis m1 e m2 dada pela equacao 4.42.

Ilustra-se a convergencia para o mınimo a partir de dois modelos iniciais

(m1o ,m2o) = (2.0,2.0) e (m1o ,m2o) = (-2.0,-2.0) gerando a mesma solucao

(0.0,0.0). Neste caso nao ha ambiguidade. . . . . . . . . . . . . . . . . . 60

4.7 (a) Funcao objetivo com duas variaveis m1 e m2. A linha negra trace-

jada representa o caminho da solucao ieterativa. A convergencia nao e

observada. (b) Resıduo em cada iteracao observado em (a). (c) Funcao

objetivo com duas variaveis m1 e m2. A linha negra tracejada representa

o caminho da solucao iterativa. A convergencia para o mınimo e suave e

direta. (d) Resıduo em cada iteracao observado em (c). . . . . . . . . . . 63

5.1 Representacao fısica de uma famılia CFP. Ela simula uma fonte na

posicao A com receptores espalhados na superfıcie. A famılia CFP regis-

tra a onda direta (linha pontilhada) e a onda refletida na interface abaixo

de A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.2 Fluxograma para estimativa do operador focal para um ponto imagem. O

processo tem inıcio com a definicao de um numero inicial m de modelos

mi sobre os quais sao calculados os respectivos operadores gi. As funcoes,

ui e Φi sao calculadas de forma a compor a matriz Ψ. O modelo final e

aquele apontado pelo maximo do modulo da funcao Ψ no lag 0. . . . . . 71

5.3 Secao de afastamento nulo exibindo dois pontos imagens em x = 12500

m. Ponto 1 no fundo do mar e ponto 2 abaixo do fundo do mar. . . . . . . 72

5.4 Painel Ψ(mi, t) para o ponto 1. O eixo horizontal corresponde a veloci-

dade variando entre 1300 e 1800 m/s. O eixo vertical representa o lag

da correlacao cruzada. A amplitude maxima no lag zero ocorre com a

velocidade de 1500 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.5 Painel Ψ(mi, t) para o ponto 2. O eixo horizontal corresponde a veloci-

dade limitada entre 1400 e 1800 m/s. O eixo vertical representa o lag

da correlacao cruzada. A amplitude maxima no lag zero ocorre com a

velocidade de 1520 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . 74

5.6 Par fonte (S ) Receptor (R) amostrando um difrator (A) na profundidade

h em um meio nao homogeneo. . . . . . . . . . . . . . . . . . . . . . . . 76

xvi

5.7 (a) Modelo de um refletor explosivo senoidal com amplitude de 25 m e

comprimento de onda igual 250 m em um meio com velocidade constante

de 2000 m/s. (b) Sismograma correspondente ao modelo em (a). A linha

branca tracejada mapeia uma difracao real no pico correspondente. Picos

e cavados da curva senoidal comportam-se como difractores; (c) Modelo

com uma serie de difratores, representados por cırculos brancos, localiza-

dos nos picos e cavados da mesma curva senoidal em (a); (d) Sismograma

correspondente ao modelo em (c). A linha branca tracejada mapeia uma

difracao real. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.8 (a) Modelo de um refletor explosivo senoidal com amplitude de 25 m e

comprimento de onda igual 500 m em um meio com velocidade constante

de 2000 m/s. (b) Sismograma correspondente ao modelo em (a). A linha

branca tracejada mapeia uma difracao real no pico correspondente. Pi-

cos e cavados da curva senoidal nao se comportam como difractores; (c)

Modelo com uma serie de difratores, representados por cırculos negros,

localizados nos picos e cavados da mesma curva senoidal em (a); (d) Sis-

mograma correspondente ao modelo em (c). A linha branca tracejada

mapeia uma difracao real. . . . . . . . . . . . . . . . . . . . . . . . . . 82

5.9 Cubo de dados sısmicos; os eixos SG e RG representam famılias de tiro

comum e receptores comuns respectivamente. O eixo vertical e o tempo,

crescente do topo para base. O plano COG representa a famılia de afasta-

mento comum para afastamento nulo. As curvas GS G, GRG e GCOG repre-

sentam o registro cinematico de uma difracao nas famılias de tiro comum,

receptor comum e afastamento (nulo) comum, respectivamente. . . . . . . 84

6.1 Experimento empregando o Projection Slice Theorem para estimar o

efeito do uso de diferentes aberturas: (a) 20000m; (b) 2000 m; (c) 1000

m; e (d) 200 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

6.2 Modelo original (modelo diapiro)). Os losangos em cada interface repre-

sentam os pontos focais. . . . . . . . . . . . . . . . . . . . . . . . . . . 96

6.3 Modelo inicial para inversao global. Os losangos em cada interface repre-

sentam os pontos focais no modelo a priori. . . . . . . . . . . . . . . . . 97

xvii

6.4 Modelo final para inversao global apos 24 iteracoes. Os losangos em cada

interface representam os pontos focais estimados com a inversao. . . . . . 97

6.5 Modelo final para inversao camada a camada. . . . . . . . . . . . . . . . 98

6.6 Modelo final estimado com a abordagem camada a camada utilizando-

se a equacao 2.22, topo-concordante, para representar cada camada. O

modelo a priori empregado foi o resultado da inversao camada a camada

presente na Figura 6.5. Os losangos em cada interface representam os

pontos focais estimados com a inversao. . . . . . . . . . . . . . . . . . . 100

6.7 Modelo final estimado com a inversao global empregando-se como mo-

delo a priori o resultado da inversao camada a camada presente na Figura

6.5. Os losangos em cada interface representam os pontos focais estima-

dos com a inversao. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

6.8 Overlay da secao de afastamento nulo com operadores CFP convencio-

nais (linhas brancas tracejadas) e operadores CFP obtidos por difracoes

(linhas brancas contınuas). O sismograma esta com tempo simples para

comportar os CFPOs. Os artefatos de modelagem nas bordas (as quatro

curvas hiperbolicas cujos apices estao nos tempos de 0.65 s, 1.05 s, 1.45

s and 1.80 s) a esquerda e a direita nao sao consideradas na analise. . . . . 104

6.9 (a) Resıduo a cada iteracao do experimento 1 para os dados sinteticos. (b)

Resıduo a cada iteracao do experimento 2 para os dados sinteticos. . . . 106

6.10 (a) Overlay do modelo real e dos CFPS calculados para o experimento

1. (b) Modelo de velocidade estimado a partir da inversao dos operado-

res CFP convencionais. (c) Overlay do modelo de velocidade real e dos

CFPs calculados no experimento 2, exibindo melhor resolucao estrutural

que o experimento 1. (d) Modelo de velocidade calculado pela inversao

dos CFPOs convencionais e DCFPOs conjuntamente. Todas as figuras

apresentam a mesma escala de cinza para velocidade. . . . . . . . . . . . 107

6.11 Linha sısmica 09 e pocos em Blake Ridge. . . . . . . . . . . . . . . . . 108

6.12 Linha sısmica 09 em tempo simples e as curvas de difracao mapeadas

(linhas brancas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

xviii

6.13 Evolucao da inversao desde o modelo a priori ate o modelo final referente

aos dados sısmicos de Blake Ridge. A escala de velocidade em m/s, a

direita, e a mesma para todos modelos. (a) Modelo a priori; (b) quarta

iteracao; (c) oitava iteracao; (d) decima segunda iteracao; (e) decima sexta

iteracao; (f) vigesima iteracao; (g) vigesima quarta iteracao; (h) vigesima

oitava iteracao. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

6.14 Resıduos da inversao tomografica dos dados monocanal de Blake Ridge. . 112

6.15 (a) Secao sısmica (linha 09) convertida para profundidade com as velo-

cidades calculadas durante a inversao. Abaixo do fundo do mar ha tres

horizontes interpretados, com linhas pretas, HZA, HZB e BSR (o bot-

tom simulating reflector (BSR) sobre a zona de gas). Os cırculos brancos

constituem as posicoes dos difratores determinadas durante a inversao dos

DCFPOs. (b) Secao composta contendo o resultado da inversao sobre os

DCFPOs da linha 09 e as velocidades projetadas dos pocos 994, 995 e

997 da esquerda (SW) para a direita (NE). Todas as interpretacoes estao

repetidas em ambas as Figuras (a) e (b). A legenda do campo de velo-

cidade e valida somente para (b). A zona hachurada na base representa

a regiao nao invertida por ausencia de dados. No texto encontram-se as

observacoes sobre os algarismos I, II, III, IV e V. . . . . . . . . . . . . . 113

A.1 Campo de velocidade sobre o qual foi realizado o experimento sismo 1.

Um campo de velocidade que obedece a equacao A.1. As linhas pretas

representam os raios emanados da fonte na posicao (x=250, z=500). . . . 130

A.2 Tempos de transito registrado ao longo da superfıcie para uma fonte de-

tonada em (x=250, z=500) no campo de velocidade da Figura A.1. . . . . 131

A.3 Espaco de resıduos com parametros a e b respectivamente na horizontal e

vertical. As curvas de nıvel indicam o valor do resıduo. . . . . . . . . . . 132

B.1 Campo de de velocidade real. . . . . . . . . . . . . . . . . . . . . . . . . 135

B.2 Campo de de velocidade do modelo 1 representado por polinomio de or-

dem 7 nas dimensoes x e z. . . . . . . . . . . . . . . . . . . . . . . . . . 136

B.3 Campo de de velocidade do modelo 1 representado por polinomio de or-

dem 8 nas dimensoes x e z. . . . . . . . . . . . . . . . . . . . . . . . . . 136

xix

B.4 Campo de de velocidade do modelo 1 representado por funcao gaussiana

em 2D (vide texto para detalhes de parametrizacao). . . . . . . . . . . . . 137

B.5 Campo de de velocidade do modelo 1 representado por serie de Fourier

restrita aos 20 primeiros harmonicos por direcao. . . . . . . . . . . . . . 139

B.6 Campo de de velocidade do modelo 1 representado por soma de senos e

cossenos restrita a 200 coeficientes por direcao. . . . . . . . . . . . . . . 139

B.7 Campo de de velocidade do modelo 1 representado por spline cubica 2D

restrita a 200 coeficientes em cada direcao. . . . . . . . . . . . . . . . . . 140

C.1 Secao de afastamento nulo destacando o ponto I (cırculo vermelho) sobre

o qual se deseja estimar o operador focal (CFP). . . . . . . . . . . . . . . 143

C.2 Painel de tiro comum (shot gather) em tempo simples onde se estima o

operador a inicial g0 destacado pela linha tracejada verde. . . . . . . . . . 144

C.3 Famılia CFP para o ponto imagem I. A linha vermelha representa o ope-

rador g0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

C.4 Painel DTS na primeira iteracao. A linha tracejada amarela representa

o evento relacionado ao ponto I. Observa-se que ele nao se encontra

alinhado ao longo do lag zero. . . . . . . . . . . . . . . . . . . . . . . . 146

C.5 Painel DTS na terceira iteracao. Observa-se que o evento em estudo esta

alinhado ao longo do lag zero. . . . . . . . . . . . . . . . . . . . . . . . 147

xx

Capıtulo 1

Introducao Geral

“Obter propriedades fısicas em regioes inacessıveis” e o objetivo dos estudos

geofısicos de acordo com muitos livros textos (Sheriff, 2002). Se a lei fısica L e as pro-

priedades fısicas m sao definidas, e possıvel calcular um campo d (o problema direto). As

ferramentas geofısicas medem campos d e a tarefa do geofısico e estimar as propriedades

m uma vez que a lei fısica L seja conhecida (o problema inverso).

Durante a historia da geofısica, uma serie de metodos foram propostos para estima-

tiva da velocidade e a sua distribuicao nas rochas. Esses metodos estao bem documenta-

dos em revistas especializadas, principalmente a Geophysics, desde 1936, e Geophysical

Prospecting, a partir de 1953. Em cada um dos metodos propostos, a lei L que relaci-

ona os parametros de rocha m e o campo d deve ser conhecida, e essa relacao e definida

genericamente por

d = Lm (1.1)

ou

m = L−1d. (1.2)

Enquanto d e calculado, uma vez que os parametros m sao conhecidos, pela equacao

1.1, usa-se a equacao 1.2, para o calculo dos parametros m a partir das medidas d.

O campo de estudos voltado para solucao da equacao 1.1 e denominado Modelagem,

tambem referido neste trabalho como Problema Direto ou ainda Modelagem Direta. O

campo de estudo voltado para solucao da equacao 1.2 e denominado Problema Inverso,

tambem referido neste trabalho como Inversao.

A resolucao da equacao 1.2 e mais complexa do que a da equacao 1.1. Se os

parametros m e a lei L sao conhecidos, o vetor d e calculado a partir de solucoes

1

analıticas ou numericas. Contudo, se o campo d e a lei L sao conhecidos e nenhuma

informacao a priori sobre m e dada, e possıvel provar que existe uma infinidade de veto-

res m aproximam-se da resposta d. Este fenomeno representa a nao unicidade observada

nos problemas inversos.

1.1 Estrutura da tese e suas contribuicoes cientıficas

A estrutura desta tese e inspirada nas equacoes 1.1 e 1.2. Assim, ha um segmento

dedicado ao modelo m, concentrado no Capıtulo 2. Nele, o comportamento das rochas em

um ambiente sedimentar e apresentado; o conceito de incoerencia geologica e proposto

assim como a equacao para sua quantificacao; diferentes formas para representacao do

campo de velocidade sao apresentadas e analisadas. Ao termino deste capıtulo ha secoes

devotadas a discussao e conclusoes.

O Capıtulo 3 e dedicado ao problema direto, ou modelagem (equacao 1.1). Tres

metodos para solucao da equacao da onda, e que sao diretamente empregados nesta tese,

sao revisitados: o metodo integral, metodo das diferencas finitas e o tracado de raios. O

capıtulo e concluıdo com uma secao de comentarios.

Os Capıtulos 4, 5 e 6 sao dedicados a equacao 1.2. No Capıtulo 4, os conceitos

basicos acerca da teoria de inversao e os metodos de inversao com enfase nos metodos de

gradiente sao revisitados. Na sequencia sao apresentadas formas para transformar proble-

mas mal postos em bem postos englobando a parametrizacao do modelo e os metodos de

regularizacao. No final do capıtulo, apresenta-se a equacao utilizada em todas as rotinas

de inversao tomografica nao linear desenvolvidas nesta tese.

No Capıtulo 5 ha tres metodos para calculo de operadores focais (Berkhout, 1997a):

o classico, por varredura, e por difracoes. Os dois ultimos sao propostas originais desta

tese.

O capıtulo 6 apresenta em detalhe a construcao da matriz de sensibilidade para a to-

mografia baseada em operadores focais envolvendo as derivadas parciais do modelo. Adi-

cionalmente, sao estudadas as limitacoes enfrentadas pela inversao tomografica devida a

abertura dos dados de entrada. Na sequencia, apresenta-se uma proposta de inversao de-

nominada tomografia sequencial, capaz de fornecer modelos verazes do ponto de vista

geologico. O capıtulo e concluıdo com aplicacoes sobre dados sinteticos e reais empre-

2

gando difracoes para estimativa do campo de velocidades.

O capıtulo 7 apresenta as discussoes finais e conclusoes. Em seguida listam-se

as referencias bibliograficas consultadas durante a pesquisa desta tese. Finalmente,

concentrou-se nos Apendices os experimentos e as deducoes.

As constribuicoes ineditas e originais angariadas com o desenvolvimento deste traba-

lho estam distribuıdas ao longo de todo o texto. Nas linhas a seguir estas contribuicoes

sao destacadas. No Capıtulo 3 sao propostas as parametrizacoes topo-concordante, base-

concordante e concordante com o centro da camada. Neste capıtulo tambem se propoe o

Indice de Incoerencia Geologica (IIG), como se ve em capıtulos posteriores, uma ferra-

menta analıtica passıvel de ser empregada na inversao.

O conceito de Modelo Cinematicamente Equivalente (MCE), uma especificacao da

ambiguidade em problemas inversos que lidam somente com o tempo de transito como

dado de entrada, e apresentada no Capıtulo 3.

No Capıtulo 4 exibe-se a equacao empregada para inversao tomografica (equacao

4.46) que, de acordo com a bibliografia consultada, e aqui aplicada de forma inedita em

problemas inversos nao lineares.

No Capıtulo 5 apresentam-se duas formas inovadoras para o calculo de operadores fo-

cais, sao elas: por varredura e atraves das difracoes. A consequencia natural da estimativa

de operadores focais a partir de difracoes e a possibilidade de se obter o campo de veloci-

dades em dados sısmicos monocanal, ate agora nada explorado seja academicamente, seja

pela industria do petroleo. Aplicacoes em dados sintetico e real monocanal no capıtulo

6 atestam a viabilidade da tecnica proposta. Ainda no Capıtulo 6, propoe-se a inversao

tomografica sequencial, uma tecnica original desenvolvida no ambito desta tese.

1.2 Fomulacao do problema e objetivos

A aquisicao de dados sısmicos voltada para investigacao do subsolo remonta o fi-

nal do seculo XIX e tinha como objetivo principal o delineamento estrutural das rochas

em subsuperfıcie. Nos anos 1920 a sısmica de refracao foi largamente aplicada para o

mapeamento de domos salinos nos EUA com finalidades prospectivas para hidrocarbone-

tos. Na decada seguinte, gradativamente, a sısmica de reflexao superou a refracao no que

tange a exploracao petrolıfera (Sheriff e Geldart, 1999).

3

Basicamente, a sısmica de reflexao com finalidades prospectivas para petroleo e do-

tada de um dispositivo composto por uma fonte mecanica com conteudo proprio de

frequencia (variando de unidades a uma centena de Hz), e uma serie de receptores dis-

tribuıdos na superfıcie. A onda emitida pela fonte viajara atraves do subsolo e a cada

contraste de propriedades elasticas ocorrera particao do sinal. A parte refletida do sinal

ascendera e sera captada pelos receptores. A outra parcela, transmitida, continuara no

caminho descendente. Novamente, a cada contraste de propriedades ocorra particao e os

receptores resgatarao na superfıcie o sinal ascendente.

Ha tres informacoes independentes contidas nos registros sısmicos (sismograma):

1) posicao: sao conhecidas as posicoes da fonte e dos receptores a cada registro;

2) tempo: cada evento, reflexao, difracao, onda direta etc. registrado nos receptores

ocorrera em um tempo. Ao longo deste texto diz-se que o tempo de um evento carrega

sua informacao cinematica; e

3) amplitude: os eventos registrados a cada tempo ocorrerao com uma determi-

nada amplitude ou magnitude do sinal. A magnitude e polaridade do sinal encerram

informacoes dos contrastes de propriedades entre as rochas. Ao longo deste texto as

informacoes de amplitude sao ditas dinamicas.

A modelagem (Lm) adequada para simular a propagacao de ondas sısmicas e repre-

sentada pela equacao 1.3 ou equacao vetorial da onda (Aki e Richards, 1980).

ρ∂2P∂t2= f + (λ + 2μ)∇(∇.P) − μ∇ × (∇ × P), (1.3)

onde P = (ux, uy, uz); ux, uy e uz sao os deslocamentos respectivamente nas direcoes x, y e

z; f sao as forcas de corpo; ∇ = ( ∂∂x ,

∂∂y ,

∂∂z ); ρ a massa especıfica e; λ e μ os parametros de

Lame.

Esta tese se concentra na resolucao do problema cinematico a partir dos tempos

de transito das funcoes de Green extraıdas dos proprios dados sısmicos ao modo do

que foi desenvolvido em Kabir (1997), Hegge (2000) e Cox (2004). Nesse contexto,

posicao e tempo sao suficientes para a inversao do campo de velocidades compressio-

nais. Nesta tese, contudo, ha uma preocupacao maior com a obtencao de um campo de

velocidades mais veraz no sentido geologico. Essa preocupacao nasce da necessidade

de se obter um modelo de velocidade que, alem de ser util para o processo de migracao

sısmica, constitua um insumo util para estimativa de litologias, esteja de acordo com as

informacoes geologicas disponıveis, e, finalmente, seja um bom modelo a priori para

4

inversoes dinamicas (eg. inversao total do campo de ondas).

1.3 Historico de pesquisas voltados para inversao de da-

dos sısmicos

A busca pela base teorica que fundamentou o estudo do Problema Inverso remete o

leitor aos grandes matematicos do seculo XIX. No tocante ao ajuste de curvas (funcoes)

citam-se indistintamente os criterios de mınimo valor absoluto de Laplace em (1799) e

mınimos quadrados de Legendre (1801) e Gauss (1809) (Tarantola, 2005). O objetivo

consiste essencialmente em se reduzir a diferenca entre dado observado e dado calculado.

A busca por mınimos (ou maximos) de funcoes e objeto de estudo da Otimizacao. Muitas

das ferramentas desenvolvidas nesse ramo da matematica se prestam para resolucao de

Problemas Inversos.

Desde ja se empregam os termos resıduo e erro respectivamente como a diferenca

entre dados observado (ou real) e calculado, e a diferenca entre modelos real e estimado,

(Aster et al. 2005). Em problemas geofısicos reais o erro nao e passıvel de calculo. Ele

so pode ser obtido em problemas sinteticos teoricos.

1.3.1 Inversao cinematica de dados sısmicos

Esta revisao nao deve comecar sem antes se definir os diferentes tipos de veloci-

dades relacionadas aos dados sısmicos e citadas ao longo deste trabalho. Citam-se aqui

seis velocidades cujas definicoes sao encontradas em livros texto de geofısica, entre eles

Sheriff e Geldart (1999): velocidade de fase; velocidade de grupo; velocidade media;

velocidade RMS; velocidade de empilhamento ou NMO; velocidade intervalar.

Ao se admitir que um meio e dispersivo, remete-se ao fato de que a velocidade de-

pende da frequencia. Entao, cada um dos harmonicos, ou fases, que compoe um pulso

sısmico viaja com uma velocidade distinta, denominada velocidade de fase. A envoltoria

dos diversos harmonicos viaja com uma velocidade que e chamada de velocidade de grupo

vg. Neste trabalho assume-se que o meio geologico nao e dispersivo, logo, todas as fases

e a envoltoria (ou grupo) possuem a mesma velocidade.

A velocidade media (vm) e um conceito emprestado da mecanica e e definida como a

5

razao entre o distancia (modulo da posicao final menos posicao inicial) e o tempo decor-

rido para um pulso viajar do ponto inicial ao final. Nao importa a trajetoria.

A velocidade RMS (vrms) e alcancada pela media quadratica (Sheriff and Geldart,

1999), definida pela expressao:

vrms,n =

√∑ni=1 v2

i ti∑ni=1 ti

, (1.4)

onde vi e a velocidade intervalar ao longo do intervalo ti.

A velocidade de empilhamento ou NMO (Normal Move Out) (vnmo) e obtida durante o

processamento sısmico e e aquela capaz de horizontalizar eventos hiperbolicos nas secoes

de afastamento medio comum (Common Mid Point Gather - CMP). Em geral a velocidade

NMO se aproxima da velocidade RMS, sendo esta ultima pouco menor que a primeira.

Para um modelo de camadas plano-paralelas, vnmo > vrms > vm.

Considerando-se camadas sedimentares com velocidades constantes, a velocidade de

cada camada e chamada de intervalar (vi). Segundo essa mesma premissa de camadas

com velocidades constantes, e supondo ser a velocidade vnmo muito proxima da vrms, Dix

(1955) demonstrou que as velocidades intervalares podem ser obtidas da vvnmo, bastando

para isso, substituir vrms por vnmo na equacao 1.4.

A complexidade do meio geologico nao encerra uma unica velocidade em uma ca-

mada. Sabe-se que esta propriedade pode variar bastante ao longo de um estrato. Assim,

ao longo deste texto o termo velocidade destituıdo de qualquer outro adjetivo diz respeito

a propriedade intrınseca da rocha.

Embora a investigacao das velocidades no subsolo remonte ao ano de 1927 com a

primeira perfilagem em um poco no Kansas (EUA) (Sheriff e Geldart, 1999), a extracao

do campo de velocidades das rochas a partir dos dados sısmicos evoluiu em desarmonia

com os estudos petrofısicos. Observa-se que depois da contribuicao de Dix (1955) houve

poucos avancos para determinacao do campo de velocidades ate meados da decada de

1980. Ha tres relevantes fatores que explicam este fato: na decada de 1960 comecou a

entrar em uso o registro digital de dados sısmicos (Sheriff e Geldart, 1999); nessa mesma

decada, viabilizou-se o emprego da migracao em tempo (Bednar, 2005) e; o fato de as

areas prospectadas nao apresentarem grandes complexidades estruturais.

A migracao em tempo envolveu esforcos na otimizacao dos recursos computacionais

e algoritmos para obtencao de boas imagens. Os principais insumos para a migracao

6

em tempo eram (sao) os dados sısmicos empilhados e o campo de velocidades NMO,

vnmo. Simultaneamente, as areas prospectadas ao redor do mundo, em geral, exibiam

poucas complexidades estruturais, e respondiam muito bem para o citado tipo de fluxo

de processamento. Grande parte dos destacados trabalhos versando sobre migracao em

tempo foram publicados nas decadas de 1950, 1960 e principalmente 1970. Citam-se

Hagedoorn (1954), Musgrave (1961), Schneider (1971), Schneider (1978) e Stolt (1978).

Na decada de 1970 tambem sao publicadas as bases teoricas voltadas para migracao

em profundidade (Claerbout 1970 e 1971) que tem como insumo os dados pre-empilhados

e o campo de velocidades em profundidade. Em 1976 e realizada a primeira migracao

em profundidade baseada no tracado de raios (Sheriff e Geldart, 1999). O campo de

velocidades obtido por meio da aplicacao da formula de Dix sobre a vnmo e a conversao

da velocidade intervalar em tempo (vi em tempo) para vi em profundidade era suficiente

para areas com baixa complexidade estrutural.

Em Pollet et al. (1983) e Bishop et al. (1985) sao reportadas as chamadas anomalias de

velocidade de empilhamento. Fazendo referencia aos levantamentos marıtimos, Bishop

et al.(1985) menciona que as citadas anomalias se fazem presentes quando ha variacoes

de velocidade com comprimento de onda inferior a dimensao dos cabos de aquisicao

marıtima (streamers). Em funcao de maiores complexidades no campo de velocidade,

nasceu a busca por tecnicas tomograficas voltadas para determinacao desse campo.

Destaca-se a publicacao pioneira de Bishop et al. (1985) versando sobre tomografia de

reflexao. Neste trabalho o dado observado compreende os tempos de transito das reflexoes

mapeadas em domınios nao empilhados. A modelagem e efetuada por tracado de raios e a

diferenca entre os tempos calculados e os tempos observados constitui a funcao objetivo,

alvo da minimizacao. O algoritmo empregado para inversao e o de Gauss-Newton, um

metodo de gradiente, no qual a cada iteracao o resıduo e reduzido a medida em que o

campo de velocidade e atualizado juntamente com a superfıcie dos refletores.

Sword (1987), que baseou seu trabalho em Rieber (1936) e Riabinkin (1957), emprega

a informacao contida no parametro p (razao entre o seno do angulo de incidencia e a

velocidade intervalar) da fonte e dos receptores, alem do tempo de transito, como dados

para inversao tomografica. Esse trabalho contem as formulacoes que alicercam a tecnica

chamada estereotomografia, citada adiante nesta secao.

Woodward et al. (2008) fazem uma retrospectiva dos metodos tomograficos a partir da

7

decada de 1990. Os autores citam que neste perıodo ha avancos na resolucao dos modelos,

tais como a incorporacao dos efeitos anisotropicos devido a utilizacao de maiores lancos

e abertura azimutal e ainda, ao maior emprego de inversoes mais guiadas pelo dado que

inversoes governadas por interpretacoes e modelos. Os autores informam tambem que

todos esses avancos somente foram possıveis com o aumento do poder computacional,

tanto em processamento como em capacidade de armazenamento (memoria) de dados.

De fato, os anos 1990 experimentam grandes avancos nas tecnicas de inversao

tomografica, motivados pelo crescente emprego da migracao em profundidade pre-

empilhamento (PSDM) e pelo aumento dos recursos computacionais. Observa-se nesse

perıodo a segmentacao da abordagem tomografica em duas vertentes principais: tomo-

grafia guiada predominantemente pelos dados, com poucas informacoes advindas de

interpretacoes e; tomografia guiada por modelos, com fortes vınculos provenientes de

interpretacoes e focada na determinacao de gradientes de velocidade entre camadas. Esta

ultima abordagem e aplicada nas inversoes camada a camada (layer-stripping), na qual os

parametros de velocidade sao determinados separadamente, da camada mais rasa para a

mais profunda. Observou-se que os resıduos acumulados nas camadas mais rasas preju-

dicam a determinacao de parametros nas camadas mais profundas (Soares Filho, 1994).

Billete e Lambare (1998) alem de publicar pioneiramente a estereotomografia co-

lige um rico historico sobre os metodos tomograficos, utilizado para confeccao deste

capıtulo. A tecnica estereotomografica adicionou mais parametros ao universo de da-

dos. Assim, constituem dados passıveis de minimizacao: o tempo de transito, posicao da

fonte, posicao dos receptores, parametro p (razao entre o seno do angulo de incidencia e

a velocidade) da fonte, parametro p do receptor. A modelagem a cada passo e realizada

por meio do tracado de raios paraxiais e os parametros a serem invertidos sao a posicao

dos refletores e o campo de velocidades.

Em funcao dos avancos nos recursos computacionais a partir da decada de 1990, pa-

ralelamente aos avancos da tomografia registram-se numerosos trabalhos voltados para

analise de velocidade para migracao (migration velocity analysis - MVA). Na sısmica de

reflexao, um mesmo evento (ponto) pode ser amostrado por diferentes famılias de tiro.

Depois da migracao, caso o campo de velocidades esteja correto, deve se situar na mesma

posicao. A diferenca na posicao dos eventos pode ser mensurada comparando-se secoes

de afastamento comum migradas (common offset gathers - COG), Trier (1990) e Dere-

8

gowski (1990), ou implicitamente avaliando-se o alinhamento de CRPs (common reflec-

tion point), Al-Yahya (1989), Symes e Carazzone (1991) e Jin e Madariaga (1994) e cons-

titui a funcao objetivo a ser minimizada. Seja qual for o domınio, almeja-se, implıcita ou

explicitamente, alinhar eventos nos paineis CRP. A Figura 1.1, extraıda de Stork (1992),

resume o processo de estimativa do campo de velocidades em domınios pos-migrados,

similares a MVA.

Figura 1.1: Descricao do processo de analise de velocidade para migracao (MVA), ex-

traıdo de Stork (1992).

9

Adicionam-se os trabalhos relacionados a tomografia em dispositivos de aquisicao

poco a poco e Vertical Seismic Profiling (VSP), entre eles McMechan (1983), Soares

Filho et al. (1997a, 1997b) empregados para determinacao do campo de velocidades e

geracao de imagens com maior resolucao que a sısmica convencional.

Como citado anteriormente, entre os anos 1990 e 2008, os avancos experimentados

na determinacao do campo de velocidade concentraram-se no aumento da resolucao dos

modelos, incorporacao multi-parametrica (anisotropia, absorcao etc.) e maior emprego de

inversoes guiadas pelo dado em detrimento de inversoes guiadas pelo modelo (Woodward

et al., 2008). Contudo, e digna de destaque uma mudanca substancial nos dados de en-

trada voltados para inversao tomografica. Berkhout (1997a, 1997b) lancou o conceito de

operadores de focalizacao (Common Focus Point Operator), CFPO, que representam a

funcao de Green de uma fonte virtual impulsiva deflagrada em profundidade e registrada

na superfıcie de aquisicao. Thorbecke (1997) e Bolte (2003) demonstram que os opera-

dores focais sao passıveis de obtencao a partir de dados sısmicos de reflexao. Assim, por

volta da segunda metade da decada de 1990 seria possıvel empregar, nao os tradicionais

tempos de transito das reflexoes, mas os tempos simples de funcoes de Green ascendentes

- dos pontos sobre uma superfıcie refletora ate a superfıcie.

Nos anos 1990 destacam-se os trabalhos de Hegge e Fokkema (1996), Hegge (1997),

Hegge et al. (1998), Hegge e Bolte (1999) e Kabir (1997) versando sobre inversao to-

mografica a partir de operadores focais. As teses de doutorado de Kabir (1997), Hegge

(2000) e Cox (2004) empregam a inversao tomografica sobre CFPO. Embora os princıpios

envolvidos ja fossem utilizados na sismologia para estimativa dos focos de terremotos e o

campo de velocidades crustal, os trabalhos de Kabir (1997), Hegge (2000) e Cox (2004)

empregam a inversao tomografica sobre CFPO (funcoes de Green) no ambito da sısmica

para petroleo. Neste processo os dados de entrada D sao os tempos de transito dos ope-

radores focais e na inversao almeja-se determinar a posicao dos difratores e o campo de

velocidades, que compoem o espaco de modelos m. O presente trabalho desenvolve-se

tendo como ponto de partida as teses citadas.

10

Capıtulo 2

O modelo (m)

2.1 Introducao

Neste capıtulo o comportamento das rochas em ambientes sedimentares e revisi-

tado. Ele embasa as formas de discretizacao mais adequadas para ambientes sedimenta-

res.

2.2 Comportamento das rochas

Diz-se que um problema apresenta de solucao unica quando ha alteracao de um

modelo m1 para outro modelo m2, os dados tambem mudam de d1 para d2 de tal forma

que d1 � d2 (Sen, 2006). A nao unicidade de solucao, contudo, e uma indesejavel ca-

racterıstica encontrada na resolucao de problemas inversos. O experimento apresentado

no Apendice A descreve claramente a nao unicidade de solucao em exemplo bastante

simples ja envolvendo o conceito de operador focal (CFPO) e ponto focal. Nesse experi-

mento uma parcela do espaco de resıduos foi construıda (Figura A.3) e nela se observam,

ao longo da diagonal, uma serie de mınimos que satisfazem os tempos de transito. Ou

seja, todos os modelos, coordenadas no espaco de resıduos, que se situam nos pontos de

mınimo sao solucoes possıveis pois satisfazem os dados observados.

Assim, e desejavel que sejam conhecidas a natureza e as caracterısticas gerais do

modelo a fim de se reduzir a multiplicidade de solucoes no Problema Inverso. O conheci-

mento a priori do modelo m norteia a melhor forma de representa-lo, seja na insercao de

vınculos (formas de restricao do espaco de modelos inseridas nas equacoes de inversao)

11

nas equacoes de inversao, seja na regularizacao (processo voltado para estabilizacao da

inversao por meio da imposicao de vınculos e restricoes que induzem a solucao).

Rocha e qualquer massa agregada firme e coerente de minerais que constitui parte

de um planeta (Skinner e Porter, 1987). Mineral constitui qualquer substancia com

composicao quımica e estrutura cristalina definida, formada por processos naturais.

Alem dos minerais constituintes, as rochas contem poros, que nada mais sao que o

espaco vazio entre os graos minerais (ou polimineralicos). Dada uma amostra de rocha,

sua porosidade e definida como a razao entre o volume de vazios e o seu volume total

(Mavko et al, 2009). No meio geologico as rochas sedimentares possuem seus poros

preenchidos por fluidos (e.g. ar, agua (doce ou salgada), gas, condensado e/ou oleo).

Os fatores intrınsecos que ditam as propriedades fısicas das rochas sao a composicao

mineral, o contato entre os graos, a forma dos graos, a porosidade e os fluidos intersticiais

(composicao e pressao de poros). Externamente, temperatura e pressao confinante ou

litostatica influenciam fortemente a velocidade sısmica nos litotipos (Mavko et al., 2009 e

Thomas, 2000). A primeira tende a reduzir a velocidade e os modulos de rigidez a medida

que a temperatura aumenta. A pressao confinante crescente age no sentido de aumentar

a velocidade pela reducao de porosidade. A Figura 2.1 exibe a curva de variacao de

velocidade de onda P com a pressao confinante. Os dados extraıdos de Carmichael (1982)

referem-se a um arenito saturado com agua e sob pressao de poro igual a zero.

Em funcao das variaveis citadas no paragrafo anterior, uma mesma litologia experi-

menta larga faixa de variacao de velocidades P. Em Mavko et al. (2009) e Carmichael

(1982) ha extensas listas contendo as velocidades de onda P para diversos litotipos sob

diferentes condicoes de temperatura, pressao confinante, tipos de fluidos e saturacao.

As rochas sedimentares formam-se a partir de sedimentos e estes sao depositados

sob a forma de camadas. Uma camada constitui um intervalo de rocha limitado supe-

riormente e inferiormente por superfıcies (interfaces), estas normalmente representando

um intervalo de nao deposicao ou erosao (Park, 1997). Sheriff (2002) define camada

como um intervalo cujas propriedades diferem do que esta acima e do que esta abaixo.

A definicao de Park (1997) tem carater genetico, informa a causa da distincao entre as

camadas, enquanto a definicao de Sheriff (2002) expressa a consequencia, que e a propria

diferenca de propriedade entre estratos. As ondas sısmicas sao sensıveis aos contrastes

de propriedades mecanicas denunciados pelas reflexoes destas ao atingirem as interfaces.

12

��

��

��

�

��

��

��

��

� ��

��

��

Figura 2.1: Curva de variacao da velocidade P em funcao da pressao de um arenito

extraıda de Carmichael (1982).

A atitude (relacao espacial de uma feicao com a horizontal, como o mergulho de uma

camada, Sheriff (2002)) original das rochas sedimentares, ou acamamento sedimentar, e

predominantemente horizontal a subhorizontal.

Em funcao das condicoes de fluxo do meio no qual os sedimentos estavam imersos e

das demais condicoes ambientais, uma mesma camada pode apresentar variacoes laterais

de propriedades ou facies (caracterısticas que distinguem uma rocha de outra adjacente,

seja na litologia, estruturas sedimentares, granulometria, etc. (Sheriff, 2002)).

Ao longo da evolucao geologica, as rochas sedimentares podem experimentar

episodios de deformacao de diversas origens (tectonica ou nao) e intrusoes ıgneas, ha-

linas ou mesmo de argila, que alteram a geometria e atitude dos estratos gerando dobras

e falhas. A Figura 2.2 resume as estruturas basicas encontradas em bacias sedimentares.

Qualquer estrutura mais complexa pode ser obtida com a combinacao de um ou mais tipos

encontrados na Figura 2.2.

13

(A)

Acamamento plano-paralelo

(B)

Acamamento plano-paralelo com

variação lateral de fácies

(C)

Dobras

(D)

Falha normal

(E)

Falha reversa

(F)

Intrusão discordante

Figura 2.2: Relacao de estruturas basicas e propriedades (cores) observadas em rochas

sedimentares. (A) Acamamento sedimentar plano-paralelo; (B) Acamamento sedimen-

tar plano-paralelo com variacao lateral de faceis; (C) Dobras; (D) Falhamento normal;

(E) Falhamento reverso; (F) Intrusao discordante.

14

2.3 Coerencia geologica

Rochas, em ambientes sedimentares, sao principalmente distribuıdas em camadas,

deformadas ou nao e, devido a compactacao, espera-se que a massa especıfica de uma de-

terminada litologia aumente com a profundidade. Admitindo que, para um meio poroso,

a velocidade e proporcional a massa especıfica (Gardner, 1974), a velocidade tambem

aumenta com a profundidade para uma rocha de composicao homogenea.

Na Figura 2.3a, expoe-se o registro bruto de velocidade extraıdo de um poco que

atravessou espessa coluna de rochas em uma bacia sedimentar brasileira. A Figura 2.3b

exibe a versao em baixa frequencia da curva mostrada na Figura 2.3a, e representa a

tendencia de compactacao desta localidade. Para uma pilha sedimentar indeformada ou

levemente deformada, a curva de velocidade em baixa frequencia espacial aumenta do

topo para a base em uma mesma camada constituıda por litotipo homogeneo (Figura

2.3b).

Devido as condicoes de sedimentacao e diagenese, o gradiente de velocidade tende a

ser paralelo ao vetor aceleracao da gravidade, ortogonal a uma superfıcie equipotencial.

Para condicoes de leito plano, o vetor de maxima variacao (aumento) de velocidade e

ortogonal a base da camada. Para uma sequencia sedimentar sobre a qual nao ha qualquer

episodio de erosao ou deformacao, o vetor gradiente de velocidade tende a ser paralelo a

normal do topo ou da base da camada.

Uma vez estabelecidas as premissas que relacionam camadas com o campo de velo-

cidade, e possvel quantificar o ındice de incoerencia geologica (IIG), um ındice original

desta tese, que deve ser computado da seguinte forma.

Primeiramente o gradiente do campo de velocidade (v) deve ser calculado.

∇v =(∂v∂x,∂v∂z

). (2.1)

Entao, o angulo θv que o vetor gradiente faz com a horizontal e obtido por

θv = arctan

⎛⎜⎜⎜⎜⎜⎝ ∂v∂z∂v∂x

⎞⎟⎟⎟⎟⎟⎠ . (2.2)

Em um terceiro passo, o vetor normal de cada ponto da interface de topo ou base de

uma camada e calculado. Em seguida, o angulo que este vetor faz com a horizontal

e extrapolado para a camada correspondente. A escolha do topo ou base da camada

deve obedecer a um criterio geologico que privilegie a interface que encerre as condicoes

15

0

1000

2000

3000

4000

5000

6000

900 1400 1900 2400 2900 3400

Profundidade (m)

Ve

loc

ida

dd

e (

m/s

)

0

1000

2000

3000

4000

5000

6000

900 1400 1900 2400 2900 3400

Profundidade (m)

Ve

loc

ida

dd

e (

m/s

)

(a)

(b)

Figura 2.3: (a) Registro de velocidade de uma pilha de rochas em uma bacia sedimentar

brasileira. (b) Versao em baixa frequencia espacial da curva (a).

16

reinantes durante a deposicao. Depois que a extrapolacao e feita para todo o modelo, ha

um campo θs que possui o mesmo significado de θv.

De acordo com as premissas previamente citadas, θs deveria ser igual a θv. O ındice

IG e computado por

IIG =∑n

i=1 abs(θvi − θsi)

n, (2.3)

onde n e o numero de celulas do modelo.

Computando θvi e θsi variando entre 0 e 180o, a situacao mais incoerente ocorre quando

θvi e θsi sao ortogonais entre si em cada celula. Neste caso, o IIG atinge o valor maximo

de 90o. Se, diferentemente, θvi e θsi sao paralelos em todo o modelo, IIG tem seu valor

mınimo em zero. Assim, o IIG varia entre 0 e 90 o para menor e maior incoerencia res-

pectivamente. A analise comparando grupos de modelos deve ser realizada sob o mesmo

numero de celulas.

Na Figura 2.4a ha um modelo representado por uma camada limitada no topo e na

base por superfıcies horizontais e o campo de velocidade exibe um gradiente vertical. θv

e θs sao comparados e sao apresentados nas Figuras 2.4b e 2.4c. O IIG, calculado com

a equacao 2.3, e 0o. Para o modelo na Figura 2.4d, os limites da camada sao os mesmos

observados em 2.4a, mas o gradiente de velocidade e lateral. Assim, para a situacao

ilustrada na Figura 2.4d, os campos θv e θs sao aqueles presentes respectivamente nas

Figuras 2.4e e 2.4f, o que causa um IIG de 90o. O modelo da Figura 2.4a e geologicamente

mais coerente que aquele observado em 2.4d.

17

Figura 2.4: (a) Camada com velocidade variando verticalmente limitada no topo e na base

por superfıcies horizontais (linhas tracejadas negra e branca respectivamente no topo e na

base). (b) θv calculado com a equacao 2.2 para (a) gerando o angulo de 90o para todo o

modelo, indicando o gradiente vertical de velocidade. (c) θs (vide o texto para maiores

detalhes) para o modelo em (a). A normal aponta para a vertical (90o). (d) Camada com

velocidade variando lateralmente, limitada no topo e na base por superfıcies horizontais

(linhas tracejadas brancas). (e) θv calculado com equacao 2.2 para (d) gerando o angulo

de 0o para todo o modelo, o que indica gradiente lateral de velocidade. (f) θs (veja texto

para maiores detalhes) para o modelo em (d). A normal aponta para direcao vertical (90o).

18

2.4 Discretizacao

Este topico dedica-se as formas de representacao do campo de propriedades m.

Este tema, por si, constitui ampla area de pesquisa e nesta secao apenas as principais

formas basicas de representacao 2D sao revisitadas.

2.4.1 Discretizacao por blocos ou celulas em malha regular

A discretizacao por blocos e o mais intuitivo e simples dos metodos de

discretizacao. Basicamente ela consiste em dividir o espaco de modelos em uma malha,

preferencialmente regular, composta por diminutos polıgonos encerrando propriedades

homogeneas, ou variaveis, contınuas. A Tabela 2.1 exibe a representacao pictorica de um

modelo por uma malha regular de retangulos.

Este tipo de discretizacao e o mais adequado para modelos extremamente complexos,

sempre que a dimensao das celulas for menor que os menores detalhes do modelo. Sob

estas condicoes, qualquer combinacao de estruturas presentes na Figura 2.2 pode ser re-

presentada por celuas. Bishop et al. (1985) e Fei e McMechan (2006), entre outros, usam

a discretizacao por celulas para estimativa da velocidade.

C1 C2 C3 C4 C5

C6 C7 C8 C9 C10

C11 C12 C13 C14 C15

C16 C17 C18 C19 C20

C21 C22 C23 C24 C25

Tabela 2.1: Representacao pictorica de um modelo discretizado por 25 retangulos.

2.4.2 Discretizacao por blocos ou celulas em malha nao regular

Malhas nao regulares constituem uma forma de discretizacao bastante atraente. Os

polıgonos nao possuem dimensoes constantes, assim, e possıvel empregar um numero

reduzido de elementos nas regioes com pouca complexidade estrutural, e elevar o seu

numero concomitantemente a reducao de suas dimensoes nas areas mais complexas e/ou

de elevados gradientes de velocidade. Rotinas de modelagem ou inversao implementadas

19

segundo o metodo dos elementos finitos permitem o emprego de malhas nao regulares

(Cook, 1995). Cox (2004) empregou triangulos de Delaunay para discretizar o campo de

velocidades no processo de inversao tomografica sobre operadores focais.

2.4.3 Discretizacao por polinomios

Nesta forma de discretizacao, sem termos cruzados, o campo de velocidades e

representado pela funcao polinomial

v(x, z) = v0 +

n∑i=1

Aixi +

m∑j=1

Bjz j, (2.4)

onde v(x, z) e o campo de velocidade e, Ai e Bj os respectivos coeficientes aplicados as

coordenadas x e z e v0 o termo independente.

Equacoes de primeiro grau, v(z) = c0+B1z , aliam a funcao compactacao (progressiva

perda de porosidade com o aumento da pressao litostatica ou profundidade) com a funcao

velocidade. Em outras palavras, para um mesmo litotipo, quanto menor a porosidade,

maior sera a velocidade de ondas compressionais. Uma serie de leis empıricas relaciona

velocidade de ondas compressionais com a porosidade. Para maiores detalhes sugere-se a

leitura de Mavko et al. (2009). Um exemplo de aplicacao do emprego de polinomio para

descricao do campo de velocidades encontra-se em Siddiqui et al. (2004). Os autores

constroem um modelo de velocidades consistente empregando equacoes de primeiro grau

para a pilha sedimentar clastica do Golfo do Mexico.

2.4.4 Discretizacao por funcao exponencial ou logarıtmica

Embora pouco comum, unidades sedimentares nao deformadas podem ser repre-

sentadas localmente por funcoes exponenciais ou logarıtmicas para grandes pacotes sedi-

mentares homogeneos. Assim uma grande sequencia arenıtica pode ser representada por

uma funcao exponencial

v(x, z) = v0 + B1e(qz), (2.5)

se ajustando ao inverso da funcao compactacao, onde q, v0 e B1 sao coeficientes que

dependem da rocha em estudo.

Algumas litologias, compostas predominantemente por graos nas fracoes silte e/ou

argila, perdem a porosidade muito rapidamente (com elevado gradiente) nas primeiras

20

centenas de metros e apresentam gradientes mais suaves em maiores profundidades (Allen

e Allen, 2005). Uma vez que a velocidade compressional relaciona-se com o inverso

da funcao compactacao, e possıvel e factıvel representar espessos pacotes de folhelhos,

argilitos ou siltitos argilosos com funcoes logarıtmicas do tipo.

v(x, z) = v0 + B1ln(z), (2.6)

onde v0 e B1 sao coeficientes que dependerao da litologia em estudo. Nao se constatou

o emprego destas funcoes (exponencial e logarıtmica) para representacao do modelo de

velocidade em tomografia nas referencias consultadas. Todavia, a Figura 2.1 que ilus-

tra a variacao de velocidade de um arenito sob diferentes pressoes confinantes se ajusta

perfeitamente a uma funcao logarıtmica.

2.4.5 Discretizacao por funcao de base radial gaussiana

As funcoes de base radial gaussiana (Buhmann, 2003) i.e.

v(x, z) =

n∑i=1

Aie− 1

2

[(x−xi)

2

σ2xi+

(z−zi)2

σ2zi

], (2.7)

onde Ai sao parametros relacionados ao peso ou apice da gaussiana, σxi e σzi relacionam-

se a abertura da superfıcie nas direcoes x e z e, finalmente xi e zi constituem as coorde-

nadas do centro da gaussiana. Se os parametros σxi e σzi forem constantes para ambas as

direcoes e para qualquer i, os unicos coeficientes necessarios para descrever um modelo

resumem-se a Ai.

2.4.6 Discretizacao por funcoes harmonicas

De forma semelhante ao topico anterior, o modelo m pode ser representado por

funcoes harmonicas como a serie de Fourier:

v(x, z) =

m∑jz=0

n∑jx=0

Ajx, jzei2π(k jx x+k jzz), (2.8)

onde Ajx, jz sao coeficientes, i =√−1, k jx e k jz sao as frequencias espaciais respectivamente

nas direcoes x e z.

A soma de cossenos e senos (equacao 2.9) sem a parcela imaginaria da equacao 2.8,

tambem pode exibir bons resultados para representacao de modelos contınuos, a depender

21

do numero de parametros empregados (anexo B).

v(x, z) =

n∑i=0

Ai[cos(kx,ix) + sen(kx,ix)

]+

m∑j=0

Aj

[cos(kz, jz) + sen(kz, jz)

], (2.9)

onde i e j sao ındices variando de 0 a n e 0 a m respectivamente, Ai o coeficiente relativo

a direcao x, A j o coeficiente na direcao z, kx,i e o numero de onda na direcao x, kz, j o

numero de onda na direcao z.

Santos e Figueiro (2006) descrevem modelos de velocidade atraves de polinomios

trigonometricos, um somatorio de funcoes harmonicas com coeficientes calculados de

forma propria e sem envolvimento de parcelas imaginarias. Os autores conseguem repre-

sentar feicoes complexas reduzindo em mais de 6 vezes (de 494 para 80) o numero de

parametros necessarios para discretizar o mesmo modelo atraves de celulas.

2.4.7 Discretizacao por splines cubicas (1D e 2D)

Observe-se primeiramente o caso 1D. Seja uma funcao C unidimensional, discreta,

definida no domınio x = [xa, xb] e representada por v=[v1, v2, v3, v4, v5, .... vn] = [v j], a

qual se deseja aproximar por funcoes de interpolacao cubica f j(x) (spline cubica) em cada

subdomınio [x j, x j+1] e com afastamento hx entre os nos do domınio. Sendo v o campo de

velocidade amostrado pontualmente, a funcao velocidade em cada subdomınio [x j, x j+1]

sera dada por:

f j(x) = Aj0 + Aj1(x − x j) + Aj2(x − x j)2 + Aj3(x − x j)

3. (2.10)

De acordo com as caracterısticas de spline em continuidade da funcao e suas derivadas

primeira e segunda (Kreyszig, 1999), advem:

Aj0 = f j(x j) = v j, (2.11)

Aj1 = f′j (x j) = bj, (2.12)

Aj2 = f′′j (x j) =

3

h2x(v j+1 − v j) − 1

hx(bj+1 + 2bj), (2.13)

Aj3 = f′′′j (x j) =

2

h3x(v j − v j+1) − 1

h2x(bj+1 + bj), (2.14)

22

onde f′j , f

′′j e f

′′′j correspondem respectivamente as derivadas primeira, segunda e terceira

de f j em relacao a x. Os parametros bj sao calculados a partir da resolucao do sistema

resumido a seguir:

bj−1 + 4bj + bj+1 =3

hx(v j+1 − v j−1), (2.15)

ou, de forma expandida: ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

b0+ 4b1+ b2 =3hx

(v2 − v0)

b1+ 4b2+ b3 =3hx

(v3 − v1)

b2+ 4b3+ b4 =3hx

(v4 − v2)

... ... ... ...

bn−2+ 4bn−1+ bn =3hx

(vn − vn−2)

(2.16)

A extensao para o caso 2D e alcancada por meio de uma nova funcao obtida com o

produto de duas funcoes 1D, uma na direcao x e outra na direcao z concorrendo em um

mesmo no (x, z). Assim, para um domınio definido por x = [xa, xb] e z = [za, zb] com nos

em x e z definidos em v = [vi, j], o calculo do valor da funcao em cada domınio e realizado

por:

fi, j(x, z) = Ai0, j0+Ai1, j1(x− x j)(z−zi)+Ai2, j2(x− x j)2(z−zi)

2+Ai3, j3(x− x j)3(z−zi)

3, (2.17)

onde i e j variam de zero a n-1 e m-1 respectivamente, os coeficientes bj e Ai, j sao obtidos

analogamente ao caso 1D.

2.4.8 Discretizacao camada a camada

Este tipo de representacao demanda a existencia de camadas individualizaveis e

o conhecimento dos limites destas. Hegge e Fokkema (1996) empregam a inversao to-

mografica global sobre operadores focais em um modelo com 5 camadas homogeneas.

Qualquer uma das funcoes apresentadas anteriormente pode ser empregada isolada-

mente ou em conjunto para compor o modelo de velocidade de uma camada. Kabir (1997)

aplicou a inversao tomografica sobre operadores focais e representou cada camada utili-

zando funcoes lineares capazes de comportar variacoes laterais de velocidade em cada

estrato. A funcao utilizada por Kabir (1997) e definida como

vi(x, z) = v0,i + [zc,i(x) − z0]Bi. (2.18)

23

Esta funcao emprega dois parametros por camada: v0,i, uma velocidade media e o gra-

diente Bi. Na equacao 2.18 zc,i(x) constitui, para cada x, a profundidade do centro da

camada i; z0 a profundidade de referencia, podendo ser o nıvel do mar.

2.5 Formas de calculo do erro

Nesta tese emprega-se o termo erro para se referir a um desajuste entre o modelo

real e o modelo calculado. Listam-se a seguir os tipos de erros empregados ao longo do

texto.

O erro quadratico medio, empregado no topico a seguir, e calculado pela soma dos

quadrados da diferenca entre modelo real vobs e modelo calculado vcalc. Para um modelo

de velocidade com nx e nz celulas respectivamente nas direcoes x e z o erro quadratico

medio (quadrado da norma L2) e calculado por

ε =

(1

nxnz

) nz∑i=1

nx∑j=1

(vobsi, j − vcalci, j

)2. (2.19)

O erro total, empregado nas aplicacoes apresentadas no capıtulo 4, e calculado por

εT =

nz∑i=1

nx∑j=1

∣∣∣vobsi, j − vcalci, j

∣∣∣, (2.20)

e o erro medio por ponto ou celula constitui o proprio erro total dividido pelo numero de

celulas (norma L1)

εm =

(1

nxnz

) nz∑i=1

nx∑j=1

∣∣∣vobsi, j − vcalci, j

∣∣∣. (2.21)

2.6 Comparacao entre as formas de discretizacao

No Apendice B, um campo de velocidades NMO, vnmo, de uma bacia sedimentar

marıtima brasileira (Modelo 1) contem 150 celulas na horizontal e 150 celulas na vertical

(Figura B.1). Este modelo foi representado por funcoes polinomiais, de base radial gaus-

siana, harmonicas e por splines cubicas em duas dimensoes, todas descritas nos topicos

anteriores, e comparadas com os dados originais qualitativamente e quantitativamente

atraves da funcao erro na equacao 2.19, onde vobsi, j e a velocidade real do modelo 1 e

vcalci, j o campo de velocidade calculado.

Nao ha duvidas de que as funcoes por domınio, as splines, apresentam melhores re-

sultados com os menores erros quadraticos (Tabela B.1 no Apendice B). Essa assertiva

24

pode ser generalizada para demais tipos de modelos que apresentem continuidade dos

gradientes e das derivadas segunda.

A segunda melhor opcao de discretizacao, em termos de erro quadratico, embora tenha

disponibilizado um modelo bastante suavizado, foi a funcao gaussiana. A terceira melhor

forma de representacao foi a obtida pela soma de senos e cossenos.

Em relacao ao erro quadratico, a quarta melhor forma de representacao foi a baseada

em polinomios de ordem 7 (7 termos em cada direcao x e z, mais o termo indepen-

dente, totalizando 15 parametros). Contudo, o reduzido numero permitido e necessario

de parametros os tornam atrativos para calculo de modelos em baixa frequencia espacial.

A elevacao do numero de parametros associada a ordem dos polinomios nao significa me-

lhor representacao do modelo, pois os coeficientes se tornam muito pequenos e instaveis,

proporcionando o fenomeno de Runge (Berrut e Trefethen, 2003). Os ensaios apontaram

que ha um numero otimo de ordem polinomial que deve ser pesquisado previamente -

neste caso, a ordem 7. A adicao de um grau na representacao polinomial elevou o erro

relativo em uma ordem de grandeza (Tabela B.1). Essa observacao se deve ao fenomeno

de Runge, responsavel por oscilacoes nos limites de um intervalo quando se emprega uma

interpolacao polinomial (Berrut e Trefethen, 2004).

A representacao do Modelo 1 por serie de Fourier truncada, de forma equivalente, a 20

termos em cada direcao x e z, apresentou as oscilacoes inerentes as funcoes harmonicas.

Esse fato elevou o erro quadratico deste tipo de funcao, tornando-o inadequado para ca-

racterizar convenientemente o Modelo 1.

2.7 Discussao e propostas para discretizacao geologica-

mente coerente

Nesta secao almeja-se unir as caracterısticas do Modelo m, listadas na secao 2.2,

de acordo com os aspectos estruturais e a distribuicao de velocidades, com as formas de

representacao apresentadas na secao 2.4. A natureza pode gerar infinitas possibilidades de

estruturas e distribuicao de velocidades. Assim, quando nao ha restricoes quanto aos graus

de liberdade (numero de parametros) para a caracterizacao do modelo de velocidade,

a melhor forma de descreve-lo e atraves de celulas com funcoes contınuas por partes.

Contudo, para inversao numerica e desejavel reduzir ao maximo o numero de parametros

25

sem prejudicar a representacao do modelo.

O emprego de funcoes contınuas de base e util para representacao de modelos uma

vez que e possıvel extrair facilmente derivadas e integrais analıticas (Apendice B).

Os resultados discutidos em detalhe no Apendice B, apontam que, dentre as funcoes

testadas, splines cubicas 2D apresentam melhores resultados para modelos de velocidade

contınuos e suaves, caracterısticas observadas em um campo de velocidades NMO. Con-

tudo, como mencionado na secao 2.2, em uma bacia sedimentar, a distribuicao das rochas

ocorre principalmente sob a forma de camadas. Camadas sao entidades cujas proprie-

dades diferem do que esta acima e abaixo (Sheriff, 2002). Para o metodo sısmico, uma

camada e uma entidade com propriedades mecanicas - velocidades de ondas P e S alem da

massa especıfica em meios isotropicos elasticos. Se o contraste entre meios adjacentes e

abrupto, a interface que separa estas camadas constitui uma descontinuidade matematica.

Em outras palavras, a derivada ao longo da direcao normal desta interface e indeterminada

ou infinita. Nenhuma das funcoes apresentadas pode representar uma descontinuidade.

Entao, a combinacao de funcoes matematicas com a representacao por camadas e a

melhor maneira de reduzir o numero de parametros honrando a principal caracterstica do

meio geologico sedimentar (a distribuicao por estratos). Kabir (1997) e Hegge (2000)

obtem otimos resultados quando combinam funcoes matematicas simples por camadas

aplicando a inversao tomografica baseada em operadores focais (CFPO). Kabir (1997)

utiliza apenas dois parametros por camada para representar extensos domınios levando

em conta variacoes laterais de velocidade. Contudo, seus resultados nao sao consistentes

sob o ponto de vista geologico. Os testes a seguir ilustram a afirmacao citada no perıodo

anterior utilizando um modelo hipotetico. O campo de velocidade de uma camada e calcu-

lado segundo a equacao de Kabir (1997) (2.18), empregando a profundidade de referencia

z0=0, o gradiente B = -1.0 e a velocidade v0=2000.0 m/s. Com estes dados obtem-se o

modelo da Figura 2.5.

26

0

50

100

150

200

250

300

350

Z(m

)

0 100 200 300 400

X(m)

1500

1600

1700

1800

(m/s

)

Figura 2.5: Campo de velocidade obtido com a equacao 2.18 de Kabir, para v0=2000

e B=-1.0.

27

De fato ha variacao de velocidade na camada, mas esta se da somente com contraste

lateral, proporcionando baixa consistencia geologica.

Nos proximos paragrafos desta secao, sao apresentadas tres formas de representacao

do modelo, que constituem variacoes da equacao 2.18 de Kabir. Nessas representacoes,

a profundidade de referencia e substituıda por uma interface que nao e necessariamente

horizontal. Esta pode ser o topo da camada, o centro da camada ou a base da camada.

Adicionalmente, em cada uma das representacoes ha uma velocidade de referencia que

coincide com o topo, com o centro ou a base do estrato.

Na primeira proposta, denominada modelo topo-concordante, empregam-se as mes-

mas condicoes de gradiente, B=-1.0, e v0=2000.0 m/s, mas considerando o topo da ca-

mada como interface de referencia, via a equacao

vi(x, z) = v0,i +[ztopo,i(x) − z

]Bi, (2.22)

onde vi(x,z) e a velocidade em qualquer ponto no interior da camada, v0,i e a velocidade

de referencia , z e a profundidade de um ponto no interior da camada e ztopo,i(x) e a

profundidade do topo da camada. Com a equacao 2.22 calcula-se o modelo presente na

Figura 2.6.

28

0

50

100

150

200

250

300

350

Z(m

)

0 100 200 300 400

X(m)

1600

1800

2000

2200

(m/s

)

Figura 2.6: Campo de velocidade obtido com a equacao 2.22 para v0=2000 e B=-1.0.

29

0

50

100

150

200

250

300

350

Z(m

)

0 100 200 300 400

X(m)

1600

1800

2000

(m/s

)


A representacao com a equacao 2.22 tem consistencia geologica e e adequado para a

unidade subjacente constituir uma intrusao.

Empregando os mesmos parametros B=-1.0 e v0=2000.0 m/s, mas considerando o

centro da camada como o horizonte de referencia, com a equacao 2.23,

vi(x, z) = v0,i +[zc,i(x) − z

]Bi, (2.23)

onde vi(x, z) e a velocidade para cada ponto, v0,i a velocidade de referencia, zc,i(x) a pro-

fundidade do centro da camada e z a profundidade de um ponto dentro da camada, o

modelo na Figura 2.7 e calculado.

30

0

50

100

150

200

250

300

350

Z(m

)

0 100 200 300 400

X(m)

1600

1800

2000

(m/s

)


A representacao com a equacao 2.23 honra a velocidade de referencia ao longo do

centro da camada. Adicionalmente ela melhor distribui o campo de velocidade no interior

da camada. Quando topo e base da camada possuem mergulhos diferentes, o padrao

do campo de velocidade obtido com a equacao 2.23 exibe gradientes que representam a

media dos mergulhos de topo e base da camada.

Finalmente, utilizando o base da camada como interface de referencia, aqui nomeados

modelo base-concordante, e empregando B=-1.0 e v0=2000.0 m/s na equacao

vi(x, z) = v0,i +[zbase,i(x) − z

]Bi, (2.24)

onde vi(x,z) e a velocidade para um ponto no interior da camada, v0,i a velocidade de

referencia, zbase,i(x) e a profundidade da base da camada e z a profundidade de um ponto

no interior da camada, calcula-se o modelo presente na Figura 2.8.

31

Dentre as quatro formas de representacao de modelos, a ultima e aquela que melhor

representa camadas sedimentares deformadas originalmente depositadas em leito plano.

O campo de velocidade e paralelo a atitude da base da camada e pode nao ser paralelo ao

topo. Idealmente espera-se que, durante a deposicao dos sedimentos, o acamamento seja

paralelo ao leito ou embasamento. Mesmo que haja episodios de deformacao com dobra-

mentos, situacao sugerida pela Figura 2.8, as propriedades tendem a seguir a geometria

da base.

A representacao pelas equacoes 2.18, 2.22, 2.23 e 2.24 indicam que e possıvel inse-

rir informacoes geologicas e reduzir o numero de parametros. Todas as representacoes

listadas neste topico sao capazes de descrever o campo de velocidades P de qualquer das

estruturas A, B, C e D presentes na Figura 2.2.

Falhas reversas e intrusoes, respectivamente estruturas E e F da Figura 2.2, encerram

complexidades difıceis de representar por camadas. Nesses casos, e necessario tratar essas

feicoes separadamente; empregando-se polıgonos com formas complexas.

2.8 Conclusoes

O experimento de discretizacao por funcoes contınuas apontou que a representacao

por spline cubica 2D e a melhor uma vez que apresentou o menor resıduo quadratico

comparadas as funcoes polinomial, base radial gaussiana, soma de funcoes harmonicas e

serie de Fourier truncada. Contudo, o reduzido numero de coeficientes necessarios para

representacao polinomial indica que essa famılia de funcoes torna-se interessante para

reduzir o numero de parametros (espaco de modelos) no processo de inversao.

A combinacao de funcoes lineares ajustando adequadamente a profundidade de re-

ferencia nas equacoes 2.22, 2.23 and 2.24 ilustram que e possıvel construir modelos con-

sistentes sob o ponto de vista geologico. Assim, a equacao 2.24 e util para descrever

camadas com gradiente de compactacao e variacao lateral deformadas devido a uma in-

trusao; a equacao 2.23 distribui o campo de velocidades crescrente do topo para a base

de uma camada tambem encerrando variacoes laterais; a equac ao 2.24 e adequada para

representar o gradiente de compactacao e mudancas laterais de velocidade em camadas

cujo topo e uma discordancia.

32

A combinacao de funcoes polinomiais de baixa ordem em cada camada constitui

poderosa ferramenta para representacao de modelos geologicamente complexos unindo

informacoes estruturais e distribuicao de velocidades.

33

Capıtulo 3

Modelagem (D = Lm)

3.1 Introducao

Esta parte da tese dedica-se a solucao do Problema Direto, equacao 1.1 voltado para

propagacao de ondas sısmicas em meios acusticos com densidade constante. Apresenta-

se, portanto, formas de resolver a equacao 3.1, um caso particular da equacao 1.1.

∇2u − 1

v2

∂2u∂t2= F, (3.1)

onde v e a velocidade de propagacao do meio, F a fonte, t o tempo, e u o campo de pressao

no domınio tempo-espaco.

Os metodos de modelagem apresentados abaixo constituem o nucleo (kernel) do pro-

cesso de inversao descrito nas proximas secoes. Assim, o principal objetivo deste capıtulo

e prover informacoes basicas que sao citadas nas secoes posteriores. Todos os metodos

para solucao da equacao 3.1 (integral, diferencas finitas e tracado de raios) sao descritos

de forma resumida com referencias classicas e recentes para que o leitor possa se apro-

fundar nesses temas. Adicionalmente, neste capıtulo se propoe o conceito de Modelo

Cinematicamente Equivalente (MCE).

3.2 Resolucao da equacao da onda por metodos integrais

Esta secao e dedicada a descricao da equacao integral voltada para simulacao do

campo de ondas acustico. Esse metodo tem maior aplicacao no ambito da sısmica para

extrapolacao reversa do campo de ondas. Ou seja, conhecido o campo de ondas na su-

34

perfıcie, calcula-se o campo de ondas no interior da Terra; ou ainda, calcula-se sua res-

posta, alterando-se virtualmente a posicao da fonte. As deducoes referentes aos metodos

integrais, resolucao do campo de ondas a partir de equacoes de integracao, baseiam-se

no Teorema de Gauss ou Princıpio da Divergencia (Mansur, 1983). O Teorema de Gauss

relaciona um campo vetorial P dentro de um volume V, limitado por uma superfıcie S,

cuja normal e n, com o fluxo P · n que atravessa S,∫V∇.PdV =

∫S

P · n dS . (3.2)

Considerando o campo vetorial P como

P = U∇G −G∇U, (3.3)

onde U e G sao dois campos escalares. Substituindo a equacao 3.3 na equacao 3.2 obtem-

se a Segunda Identidade de Green:∫V

(U∇2G −G∇2U

)dV =

∫S

(U∇G −G∇U) · n dS . (3.4)

A equacao 3.4 tem aplicacao generalizada. Em problemas acusticos U e G sao governados

pela equacao de Helmholtz,

∇2U +ω2

v2(x)U = 0, (3.5)

onde ω e a frequencia angular, x o vetor posicao e v(x) o campo de velocidade.

As deducoes a seguir estao desenvolvidas para meios com massa especıfica cons-

tante, o que constitui uma simplicacao nao considerada em Berkhout e Wapenaar (1989)

e Wapenaar et al. (1989). Apesar do rigor matematico empregado nos citados trabalhos, a

aproximacao acustica em meios geologicos nao honra as amplitudes observadas em dados

reais. Assim, as deducoes a seguir buscam honrar somente as informacoes cinematicas

trazidas pelos operadores CFP. Entao, a aproximacao para um meio com massa especıfica

constante e suficiente.

Considere um campo de ondas propagando-se para cima e sendo registrado em uma

superfıcie nao reflexiva S 0 (Figura 3.1). O meio abaixo de S 0 e caracterizado por um

campo de velocidade [v(x)] suave e sua massa especıfica e constante. O contorno S 1

encontra-se no infinito, de tal forma que qualquer campo que o atinja nao tem influencia

sobre A (Figura 3.1), em acordo com a condicao de radiacao de Sommerfeld (Mansur,

1983; Wapenaar et al., 1989). Com estas premissas, tem-se a seguinte funcao de Green

35

∇2G +ω2

v(x)2G = �[δ(x − xA)], (3.6)

onde o operador � corresponde a transformada direta de Fourier, xA significa o vetor

posicao do ponto A, local da fonte virtual abaixo de S 0.

Figura 3.1: Representacao de um campo de ondas propagando-se para cima e regis-

trado na superfıcie S 0. Deseja-se estimar o campo de ondas em A conhecendo-se o

campo registrado na superfıcie S 0.

Substituindo as equacoes 3.5 e 3.6 na equacao 3.4 tem-se:∫V

U�[δ(x − xA)] − U[ω2

v(x)2G] −G[− ω

2

v(x)2U] dV =

∫S 0

(U∇G −G∇U) · n dS (3.7)

e simplificando o lado esquerdo,∫V

U�[δ(x − xA)] dV =∫

S 0

(U∇G −G∇U) · n dS . (3.8)

Empregando as propriedades da funcao delta de Dirac na posicao xA, o lado esquerdo se

torna:

U(xA) =

∫S 0

(U∇G −G∇U) · n dS , (3.9)

onde U(xA) e o campo escalar de ondas em A.

Para propagacao inversa (Wapenaar et al., 1989) faz-se necessario utilizar o conjugado

G∗, de G. Assim, a equacao 3.9 torna-se

U(xA) =

∫S 0

(U∇G∗ −G∗∇U) · n dS . (3.10)

36

Como a superfıcie S 0 e plana, sua normal e paralela ao eixo z

U(xA) =

∫ ∞

−∞

∫ ∞

−∞U∂G∗

∂z−G∗

∂U∂z

dxdy. (3.11)

Seguindo Wapenaar et al. (1989), a equacao 3.10 torna-se

U(xA) = −2

∫ ∞

−∞

∫ ∞

−∞U∂G∗

∂zdxdy. (3.12)

Considerando G=�(g) e U=�(u), como u e g sao respectivamente o campo escalar de

ondas e a funcao de Green no domınio do tempo, e tomando a equacao 3.12 no domınio

do tempo,

u(xA) = −2

∫ ∞

−∞

∫ ∞

−∞u ⊗t∂g∂z

dxdy. (3.13)

3.3 Resolucao da equacao da onda por metodos diferen-

ciais

Na secao anterior observou-se que, conhecendo-se a funcao de Green de um meio

e o campo de ondas ao longo de uma superfıcie infinita S0, e possıvel calcular o campo

de ondas subjacente a ela utilizando de um metodo integral. Pelo Metodo das Diferencas

Finitas (MDF), e possıvel calcular o campo de ondas u, equacao 3.1, conhecendo-se as

condicoes inciais e de contorno.

A modelagem pelo MDF demanda a discretizacao do campo de velocidades por meio

de celulas, regulares ou nao, descritas na secao 2.4.1. Na formulacao mais simples

possıvel empregando-se operadores de segunda ordem no tempo e no espaco (trunca-

mento no termo de ordem 1 na serie de Taylor) em meios tridimensionais, a equacao 3.1

e resolvida iterativamente da seguinte forma:

ux,y,z,t+1 = v2x,y,zΔt2[(ux+1,y,z,t − 2ux,y,z,t + ux−1,y,z,t

Δx2

)

+

(ux,y+1,z,t − 2ux,y,z,t + ux,y−1,z,t

Δy2

)

+

(ux,y,z+1,t − 2ux,y,z,t + ux,y,z−1,t

Δz2

)

] + 2ux,y,z,t − ux,y,z,t−1.

(3.14)

37

Foge ao escopo desta tese avancar na seara de modelagem. A equacao apresentada

acima e empregada nas rotinas de simulacao que geraram alguns dos dados sinteticos

utilizados neste trabalho.

O aprofundamento do tema modelagem sısmica por MDF demanda a consulta de

bibliografias basicas como Kelly et al. (1976) que, embora seja voltada para ondas

elasticas, ilustra didaticamente o processo de implementacao para meios homogeneos e

heterogeneos. Virieux (1986) e Levander (1988), ambos tambem dedicados a simulacao

em meios elasticos, descrevem a implementacao da modelagem com o emprego de malha

intercalada.

Estudos de estabilidade, emprego das condicoes de contorno assim como formas de

aplicacao e tipos de fontes sao abordados em Cunha (1997) que, lidando com o tema da

migracao reversa no tempo, estudou os citados temas detalhadamente. Bulcao (2004), nas

secoes e apendices dedicados a modelagem, apresenta estudos de estabilidade, emprego

de operadores de ordem superior, alem de detalhes de implementacao para modelagem

acustica e elastica em meios heterogeneos em duas e tres dimensoes. Carcione (2007)

aborda a teoria e formas de implementacao de meios complexos envolvendo, alem de

meio acustico, o elastico, o poroelastico, o visco-elastico e anisotropias. Destaca-se,

adicionalmente Kelly e Marfurt (1987) e, mais recentemente, o suplemento dedicado a

modelagem introduzido por Robertsson et al. (2007).

3.4 Aproximacao da equacao da onda por tracado de

raios

O metodo de tracado de raios permite calcular o campo de ondas a partir da solucao

assintotica de alta frequencia da equacao 3.1. No ambito da tomografia estudada neste

trabalho, interessam somente as informacoes cinematicas trazidas pelos raios: tempo de

transito, posicao e direcao do raio. As deducoes referentes ao tracado cinematico de

raios sao encontradas em Botelho (1986), Soares Filho (1994), Cerveny (2001) e Cunha

(2005). As equacoes diferenciais para modelagem direta do tracado de raios a partir de

um campo de velocidade conhecido, e estabelecidas as condicoes iniciais de posicao e

38

cossenos diretores sao:

dxdl= vp, (3.15)

dpdl= ∇

(1

v

), (3.16)

dτdl=

(1

v

), (3.17)

onde x e o vetor posicao (x = (x, y, z) ), dl e o comprimento de raio infinitesimal, v e

a velocidade na posicao (x, y, z) (v = v(x, y, z)), p e o parametro de raio contendo os

cossenos diretores (p = (px, py, pz)) e τ e o tempo de transito.

Durante o desenvolvimento desta tese, uma rotina de modelagem baseada em tracado

de raios foi implementada. Ela emprega a formulacao de Runge-Kutta (Cunha, 2005)

e constitui a ferramenta de modelagem, ou kernel, das rotinas de inversao descritas nas

proximas secoes.

3.5 Modelos cinematicamente equivalentes

Na secao anterior foram apresentadas formas de resolver a equacao da onda.

Mesmo que os metodos integrais e diferenciais apresentem estimativas corretas do tempo

de transito (informacao cinematica) e amplitudes (informacao dinamica), o objetivo deste

trabalho e estimar o campo de velocidades P a partir dos tempos extraıdos de operadores

focais. Assim, o kernel da rotina de inversao tomografica adotada e uma sub-rotina de

tracado de raios.

O tempo de transito para um raio viajar de um ponto A para um ponto B e governado

pela equacao 3.17 ao longo da trajetoria deste raio. Ha diferentes campos de velocidade

que fornecem o mesmo tempo de transito causando, portanto, ambiguidades. Modelo

Cinematicamente Equivalente (MCE) e qualquer modelo contınuo de velocidade e estru-

turas cujos dados calculados se ajustam aos dados registrados em um experimento.

Mesmo um simples ensaio representado por duas celulas (Figura 3.2) ha infinitos

MCEs. O tempo de transito t de um raio horizontal viajando entre a fonte no lado es-

39

v1 v2

DX DX

Figura 3.2: Representacao de um modelo de duas celulas com velocidades v1 e v2. A es-

trela representa a fonte de um raio, que viaja desta ate o receptor, o triangulo. A dimensao

da malha quadrada e DX.

querdo (estrela) ate o receptor no lado direito (triangulo) e dado por:

t =DXv1

+DXv2

, (3.18)

onde DX e a dimensao da celula, v1 a velocidade da celula a esquerda e v2 a velocidade

da celula a direita. Qualquer outro par de velocidades vA e vB cuja soma dos inversos

satisfaca a condicao tDX e um MCE desse modelo hipotetico.

O experimento descrito no Apendice A apresenta um exemplo de MCE encerrando o

conceito de operador focal.

3.6 Comentarios

Todas as rotinas de modelagem sısmica lidando com metodos integrais (secao 3.2),

MDF (secao 3.3) e tracado de raios (secao 3.4) foram desenvolvidas e implementadas

durante o desenvolvimento desta tese. Contudo, nao ha contribuicoes alem do que ja fora

desenvolvido nos trabalhos citados nas secoes deste capıtulo.

Este capıtulo encerra a proposta e definicao de conceito de Modelo Cinematicamente

Equivalente. Ele especifica a ambiguidade, ou multiplicidade de modelos que satisfazem

40

um conjunto de dados, que ocorre durante a inversao de operadores focais.

41

Capıtulo 4

Inversao (m=L−1 D)

4.1 Introducao

Esta secao dedica-se a solucao da equacao 1.2 voltada para problemas cinematicos

acusticos da sısmica de reflexao. Antes de se visitar os metodos de resolucao, apresentam-

se a seguir informacoes e conceitos gerais relacionados ao Problema Inverso.

Hadamart (1902) elenca tres propriedades que um modelo matematico deve possuir

para representacao de um fenomeno fısico:

1) Existencia de solucao: deve-se averiguar se o sistema m = L−1d, equacao 1.2, tem

solucao. Em, outras palavras, se ha um m tal que d = Lm.

2) Unicidade: a nao unicidade implica em dizer que diferentes modelos mi sao capazes

de gerar um mesmo dado d.

3) Estabilidade: a estabilidade indica como pequenos erros nos dados se propagam no

modelo calculado. Quanto mais estavel e o problema menos sensıvel ele e a pequenos

erros nas condicoes iniciais.

Modelos matematicos que possuem as propriedades acima sao ditos bem postos.

Como citado anteriormente, a solucao da equacao 1.2 e mais complexa, pois Problemas

Inversos, normalmente, sao mal postos. Em geral os Problemas geofısicos apresentam

grande instabilidade, principalmente devido ao mal condicionamento (operadores L com

elevado numero de condicao). A estabilizacao pode ser alcancada pela regularizacao do

problema, processo no qual sao impostos vınculos que induzem a solucao.

O erro foi previamente definido na secao 1.3 como a diferenca entre modelos real e

42

calculado. Assim, calcula-se o vetor erro por

ε = mreal −mcalc, (4.1)

onde mreal e mcalc sao respectivamente os modelos real e calculado.

Para um universo de m amostras o erro total e calculado por

εt =

m∑i=1

∣∣∣mreali − mcalci

∣∣∣ , (4.2)

e o erro relativo por

εt =1

m

m∑i=1

∣∣∣mreali − mcalci

∣∣∣ . (4.3)

O vetor resıduo r constitui a diferenca entre o dado observado e o dado calculado.

r = dobs − dcalc. (4.4)

Substituindo a equacao 1.1 na equacao 4.4 obtem-se:

r = dobs − Lm. (4.5)

A resolucao do problema inverso reside na busca pelos extremos de um funcional que

meca a norma entre dados observados e dados calculados, a funcao objetivo (Wolfram,

1999). Uma funcao objetivo Θ simples e desprovida de qualquer informacao a priori

sobre o modelo pode ser obtida a partir da norma 2 do vetor resıduo (r):

Θ =[dobs − Lm

]T [dobs − Lm

].. (4.6)

Em funcao da relacao que existe entre o dado d e o modelo m classificam-se os pro-

blemas inversos em: lineares e nao lineares. Em problemas inversos lineares o dado d e

o modelo m relacionam-se linearmente por meio do operador L, tal como representado

pela equacao 1.1 (Aster et al., 2005; Menke, 1984). Em problemas nao lineares a relacao

expressa pela equacao 1.1 e valida, contudo, o operador L e funcao dos parametros m

(Aster et al., 2005).

Problemas fracamente nao lineares sao aqueles cuja perturbacao nos dados Δd

relaciona-se com a perturbacao no modelo por L pela equacao 4.7 (Sen, 2006).

Δd = LΔm. (4.7)

43

Nesta classe de problemas inversos (fracamente nao lineares), a funcao objetivo pode ser

aproximada por uma funcao quadratica na vizinhanca de um modelo mk (Menke, 1984),

levando a uma solucao mk+1 dentro da regiao de convergencia quadratica. Uma serie de

passos localmente lineares constitui uma aproximacao para a solucao fracamente nao-

linear (equacao 4.7. Devido ao carater nao linear desta classe de problemas inversos, o

operador L muda a cada iteracao.

Nesta tese, todos os problemas abordados sao considerados lineares ou fracamente

nao lineares.

4.2 Metodos de resolucao do problema inverso

Ha uma mirıade de metodos (Menke, 1984; Tarantola, 1984; Tarantola, 2002;

Aster, 2005; Tarantola, 2005; Sen, 2006) voltados para solucao do Problema Inverso

(equacao 1.2), que pode ser resumida a otimizacao, obtencao de mınimo ou maximo da

funcao objetivo. Por ora classificam-se as solucoes nos seguintes grupos basicos, aqui

modificada a partir de Tarantola (2002):

- Metodos baseados na busca sistematica, ou varredura, no espaco de modelos;

- Metodos de tentativa e erro;

- Metodos de busca aleatoria no espaco de solucoes (ex. Monte Carlo, Simulating

Annealing, algoritmo genetico);

- Metodos baseados em gradientes.

Os metodos baseados na busca sistematica, ou varredura, no espaco de modelos con-

sistem em se construir uma malha regular nesse espaco e computar o valor da funcao

objetivo em todos os pontos desta malha (Tarantola, 2002). Espera-se que aquele modelo

que possuir menor resıduo seja a solucao do problema inverso. O experimento presente

no Apendice A descreve a construcao da funcao objetivo pela varredura do espaco de

modelos. Esse metodo e indicado quando o numero de parametros e pequeno.

O metodo de tentativa e erro consiste em se alterar as variaveis do modelo em uma

base intuitiva e calcular os dados ate que a diferenca entre dados reais e calculados, o

resıduo, situe-se dentro de uma tolerancia aceitavel (Tarantola, 2002). Este metodo, de-

finitivamente, nao e passıvel de realizacao quando o espaco de modelos e muito grande.

44

Adicionalmente, ele depende da experiencia e conhecimento do profissional que esta rea-

lizando a inversao.

A busca aleatoria no espaco de modelos, ou metodos tipo Monte Carlo, sao aqueles

que empregam em algum estagio do processo inverso a escolha aleatoria de m para con-

secutivo calculo dos dados (Tarantola, 2002). Essa classe de metodos e particularmente

util para solucao de problemas altamente nao lineares. Uma segunda virtude dos metodos

tipo Monte Carlo e o fato de eles constituirem algoritmos de busca global, fato aumentam

as chances de convergencia para um mınimo global.

Agrupou-se nos metodos baseados em gradiente, aqueles que fazem uso do gradiente

da funcao objetivo (resıduo, funcao probabilidade ou outra medida que reflita o desajuste

entre dado observado e dado calculado) na busca do seu mınimo ou maximo global. Essa

classe de metodos baseia-se no metodo de gradiente de Cauchy, datado de 1847 (Kreyszig,

1999). Dentre os metodos de gradiente citam-se: Newton, Gauss-Newton e gradiente

conjugado.

Nesta tese, foram explorados os metodos baseados na busca sistematica no espaco

de modelos e os metodos baseados em gradientes. A proxima secao detalha os metodos

baseados em gradientes.

4.3 Resolucao do problema inverso pelos metodos de gra-

dientes

Abordam-se neste topico os principais metodos de resolucao de Problemas Inversos

fracamente nao lineares por meio de informacoes advindas ou relacionadas ao gradiente

da funcao objetivo.

4.3.1 Metodo da maxima declividade

Sendo o gradiente da funcao objetivo, ∇Θ = ∂Θ∂m , busca-se o mınimo iterativamente

pela equacao

mk+1 = mk − a∇Θ, (4.8)

onde a e um escalar positivo. Ou seja, o modelo no passo k + 1 e obtido a partir do

modelo no ponto k, mk na direcao contraria ao gradiente da funcao objetivo ∇Θ. A

45

Figura 4.1: Funcao objetivo com duas variaveis m1 e m2 ilustrando as iteracoes pelo

metodo de maxima declividade ate a convergencia para solucao em m (modificado de

Shewchuck, 1994).

figura 4.1 ilustra a funcao objetivo e o caminho das iteracoes, entre o modelo inicial m0

e a solucao m, segundo o metodo da maxima declividade. Este constitui o metodo de

gradiente mais simples, exibe baixa taxa de convergencia (Flechter, 1980), contudo ela

e assegurada (Aster, 2005), desde que o modulo do gradiente seja ajustado segundo uma

constante ad hoc.

4.3.2 Metodo de Newton

A expansao em serie de Taylor (truncada no termo de ordem 2) da funcao objetivo

nas vizinhancas de mk fornece:

Θ(mk+1) = Θ(mk) + ∇Θ(mk)[mk+1 − mk] +1

2∇2Θ(mk)[mk+1 − mk]

2. (4.9)

Considere Δm = mk+1 − mk. Almeja-se o ponto de mınimo de 4.9 em relacao a

Δm. Suprimindo-se os passos intermediarios obtem-se a resolucao do modelo mk+1 pelo

46

metodo de Newton (Tarantola, 2002) por:

mk+1 = mk − [∇2Θ(mk)]−1∇Θ(mk) (4.10)

A aplicacao do metodo de Newton para a resolucao de Problema Inverso demanda o

calculo do Hessiano ∇2Θ ou ∂2Θ∂m2 , a derivada segunda da funcao objetivo em relacao ao mo-

delo m. Em alguns problemas da geofısica este calculo e complexo. Adicionalmente nao

ha garantia de convergencia (Aster et al., 2005). Para problemas passıveis de linearizacao

o metodo de Gauss-Newton e mais factıvel.

4.3.3 Metodo de Gauss-Newton

A deducao do metodo em escopo pode ser encontrada em diversos livros textos,

entre eles, Flechter (1980), Menke (1984), Nocedal e Wright (2000) e Aster et al. (2005),

contudo ela e apresentada a seguir completa e comentada para compreensao do leitor e

futuras consultas do autor.

Avancando na equacao 4.4 obtem-se:

Θ =1

N

N∑i=1

[dobs

i − dcalci

]2(4.11)

Expandindo o termo dcalci da equacao acima em serie de Taylor nas proximidades de um

modelo m0,

Dcalci = Dcalc

i (m0) +

M∑j=1

∂di(m0)

∂mjΔmj, (4.12)

onde Δm j e uma perturbacao no modelo nas proximidades de m0. Derivando a equacao

4.11 em relacao ao modelo,

∂Θ

∂mk=

2

N

N∑i=1

(dobs

i − dcalci

) ∂∂mk

(dobs

i − dicalc

)(4.13)

No ponto de mınimo ∂Θ∂mk= 0. Adicionalmente, os dados observados, dobs

i , nao dependem

do modelo m. Assim a equacao 4.13 torna-se

0 =N∑

i=1

(dobs

i − dcalci

) ∂dcalci

∂mk(4.14)

Substituindo a equacao 4.12 em 4.14 obtem-se:

N∑i=1

⎛⎜⎜⎜⎜⎜⎜⎝dobsi − dcalc

i (m0) −M∑j=1

∂di(m0)

∂mjΔmj

⎞⎟⎟⎟⎟⎟⎟⎠ ∂dcalci

∂mk= 0 (4.15)

47

Considerando dobsi - dcalc

i (m0) = Δdi:

N∑i=1

⎛⎜⎜⎜⎜⎜⎜⎝Δdi,k −M∑j=1

∂dcalci

∂mjΔmj

⎞⎟⎟⎟⎟⎟⎟⎠ ∂dcalci

∂mk= 0. (4.16)

Rearranjando os termos,

N∑i=1

Δdi,k∂dcalc

i

∂mk−

N∑i=1

M∑j=1

∂dcalci

∂mjΔmj∂dcalc

i

∂mk= 0. (4.17)

A equacao 4.17 pode ser escrita na forma matricial se for considerado

L =

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

∂dcalc1

∂m1

∂dcalc1

∂m2...

∂dcalc1

∂mM

∂dcalc2

∂m1

∂dcalc2

∂m2...

∂dcalc2

∂mM

... ... ... ...

∂dcalcN∂m1

∂dcalcN∂m2

...∂dcalc

N∂mM

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠, (4.18)

e

LT =

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

∂dcalc1

∂m1

∂dcalc2

∂m1...

∂dcalcN∂m1

∂dcalc1

∂m2

∂dcalc2

∂m2...

∂dcalcN∂m2

... ... ... ...

∂dcalc1

∂mM

∂dcalc2

∂mM...

∂dcalcN∂mM

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠. (4.19)

Logo,

LTΔd − LT LΔm = 0, (4.20)

ou,

LT LΔm = LTΔd. (4.21)

Aplicando (LT L)−1 em ambos os lados:

Δm = (LT L)−1LTΔd. (4.22)

A equacao 4.22 e a forma basica mais empregada nos Problemas Inversos

considerando-se a funcao objetivo com norma 2. Ela assemelha-se a solucao da equacao

4.10 aproximando o Hessiano pelo termo LTL (Nocedal e Wright, 2000). Contudo, ha

que se relembrar duas premissas embutidas em sua solucao: (1) A modelagem e o opera-

dor L que regem o fenomeno fısico devem ser passıveis de linearizacao nas vizinhancas

de um modelo m0; e (2) O modelo m0 deve estar proximo da solucao.

A implicacao direta da primeira premissa e que, nas aproximacoes extremamente nao

lineares o metodo falha. Embora a direcao apontada pelos cossenos diretores de Δm nas

48

vizinhancas de m0 sejam corretas, seu modulo pode superar as condicoes de linearizacao

embutidas na expansao em Taylor. Neste caso, se a magnitude de Δm e maior que este

limite, deve-se reduzir o seu modulo utilizando de um parametro ad-hoc, seja ε. Este

parametro e positivo, menor que 1 (Tarantola, 2002) e pode variar em cada iteracao. As-

sim, a equacao 4.22 torna-se

Δm = ε(LT L)−1LTΔd, (4.23)

aumentando as chances de convergencia. Um valor desnecessariamente baixo de ε pode

tornar a convergencia extremamente lenta.

4.3.4 Metodo de Levenberg-Marquardt

A aproximacao do Hessiano por LT L contribui para as incertezas na magnitude

de Δm no metodo de Gauss-Newton. No metodo de Levenberg-Marquardt evita-se o

emprego da constante ε, uma especie de busca linear (line search), empregando-se uma

estrategia de regiao de confianca (Nocedal e Wright, 2000). Esta regiao de confianca e

obtida pela da adicao de uma matriz diagonal, fruto do produto de uma constante ψ com

a matriz identidade I, a matriz LT L.

Δm = (LT L + ψI)−1LTΔd, (4.24)

Uma vantagem adicional dessa estrategia e o fato de ela contornar possıveis de-

ficiencias no rank (numero de linhas ou colunas linearmente independentes) da matriz

L, uma das fraquezas do metodo de Gauss-Newton (Nocedal e Wright, 2000). Em outras

palavras, o metodo em escopo forca que a matriz LT L seja positiva definida garantindo

sua nao singularidade.

4.3.5 Pseudo-Inversa de Moore-Penrose

Para resolver a equacao 1.2 ou 4.7, e necessario inverter a matriz L. Contudo, caso

esta matriz nao seja inversıvel, em lugar de se calcular m ou Δm pela equacao 1.2, L

e decomposta em valores singulares (Singular value decomposition (SVD)) (Aster et al.

2005).

A matriz com valores singulares, Σ, e uma matriz diagonal com valores positivos e

decrescentes de acordo com o aumento de i (Σi=1 > Σi=2 > Σi=3 > ..). Ela pode possuir

49

alguns elementos muito baixos ou nulos. Se somente os p primeiros elementos nao nulos

de Σ forem considerados, esta matriz pode ser particionada da seguinte forma (Aster et

al., 2005): ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣Σp 0

0 0

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦ , (4.25)

onde Σp corresponde a uma matriz quadrada diagonal composta por apenas valores sin-

gulares positivos nao nulos. A matriz L pode ser representada por (Aster et al. 2005),

L =W

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣Σp 0

0 0

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦ VT . (4.26)

Se as matrizes U e V na equacao 4.26 forem expandidas, obtem-se:

L =

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

W1,1 W1,2 ... W1,p ... W1,m

W2,1 W2,2 ... W2,p ... W2,m

... ... ... ... ... ...

Wp,1 Wp,2 ... Wp,p ... Wp,m

... ... ... ... ... ...

Wm,1 Wm,2 ... Wm,p ... Wm,m

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣Σp 0

0 0

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

V1,1 V1,2 ... V1,p ... V1,n

V2,1 V2,2 ... V2,p ... V2,n

... ... ... ... ... ...

Vp,1 Vp,2 ... Vp,p ... Vp,n

... ... ... ... ... ...

Vn,1 Vn,2 ... Vn,p ... Vn,n

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

T

.

(4.27)

Observa-se que na multiplicacao matricial de W por Σ, todos os elementos de Wi, j tais

que j > p estarao multiplicados pelos valores nulos de Σ. Da mesma forma se sucedera

com a matriz V. Assim, empregando a notacao de Aster et al. (2005), reescreve-se o

produto matricial levando-se em conta somente os elementos de W, Σ e V que darao re-

sultados nao nulos. Neste caso, somente as p primeiras colunas de W, Wp, as p primeiras

colunas de V, Vp e somente os elementos nao nulos de Σ, Σp. Logo L pode ser calculada

por:

L =

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

W1,1 W1,2 ... W1,p

W2,1 W2,2 ... W2,p

... ... ... ...

Wp,1 Wp,2 ... Wp,p

... ... ... ...

Wm,1 Wm,2 ... Wm,p

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

[Σp

]

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

V1,1 V1,2 ... V1,p

V2,1 V2,2 ... V2,p

... ... ... ...

Vp,1 Vp,2 ... Vp,p

... ... ... ...

Vn,1 Vn,2 ... Vn,p

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

T

. (4.28)

50

A equacao acima escrita de forma compacta torna-se

L =WpΣpVTp . (4.29)

O produto WpΣpVTp e passıvel de inversao de acordo com Moore (1920) e Penrose

(1950). Assim a inversa generalizada de L ou pseudo-inversa de Moore e Penrose (Aster

et al., 2005) e calculada por

L† = VpΣ−1p WT

p . (4.30)

No lugar da equacao 1.2 emprega-se

m = L†d, (4.31)

desde que o posto de L seja igual ao numero de parametros do modelo m e inferior ao

numero de registros ou dados n.

Em problemas fracamente nao lineares a mesma relacao e valida para perturbacoes no

modelo e no dado,

Δm = L†Δd, (4.32)

, e a solucao para o modelo m e alcancada iterativamente,

mk+1 = mk + Δm, (4.33)

ate que dado modelado e dado observado se ajustem.

4.3.6 Metodo dos gradientes conjugados

O metodo em escopo faz parte da famılia de metodos denominada direcoes conju-

gadas (dois vetores di e dj sao conjugados ou L-ortogonais, se diLdj = 0, onde L e uma

matriz compatıvel com o produto (Shewchuck, 1994)). Tratam-se de metodos aplicaveis

a matrizes positiva definidas (para qualquer vetor x diferente de zero tem-se xT Lx > 0) e

que sao mais efetivos no tratamento de matrizes esparsas (Shewchuck, 1994).

O ponto de partida dos metodos de direcoes conjugadas nao difere do metodo de gra-

dientes. Contudo, enquanto neste ultimo observam-se varios passos alternados seguindo

uma mesma direcao ate que se atinja a solucao (Figura 4.1), os metodos das direcoes

conjugadas empregam buscas lineares (line search) em direcoes especıficas, e, ao longo

destas, busca-se o mınimo. Este processo reduz o numero de iteracoes.

51

No metodo dos gradientes conjugados, descrito no pseudo-algoritmo a seguir, algo-

ritmo 4.34, a direcao de busca e construıda atraves da ortogonalizacao dos resıduos se-

gundo a matriz L (L-ortogonalizacao) (Shewchuck, 1994).

Δm0 = r0 = Dobs − Lm0

loop while (r > tol)

αk =rT

k rk

ΔmTk LΔmk

mk+1 = mk + αkΔmk

rk+1 = rk − αkLΔmk

βk+1 =rT

k+1rk+1

rTk rk

Δmk+1 = rk+1 + βk+1Δmk (4.34)

end loop

Na primeira iteracao assume-se que a atualizacao do modelo (Δm0) e igual ao resıduo

(r0). A constante α calculada a cada iteracao k pondera o modulo da atualizacao por

(αkΔmk). O calculo βk culmina o processo de ortogonalizacao apos o qual calcula-se o

valor efetivo da atualizacao do modelo Δmk+1. O processo iterativo continua ate que o

resıduo rk esteja abaixo de uma tolerancia tol.

4.4 Influencia do dispositivo de aquisicao

Todos os metodos de gradiente apresentados na secao anterior demandam o calculo

do operador L, aqui referido como matriz de sensibilidade, que avalia a influencia de cada

parametro sobre os dados. No experimento hipotetico representado pela Figura 4.2, um

modelo e discretizado segundo quatro celulas quadradas cujos lados posssuem dimensao

l.

52

Figura 4.2: (a) Modelo de velocidade discretizado por celulas quadradas;

(b) Aquisicao sısmica hipotetica avaliando cada celula com um par fonte-

receptor, representados respectivamente pela estrela estrela vermelha e o

triangulo preto. As flechas pontilhadas indicam o caminho percorrido pelo

raio da fonte para o receptor.

Segundo experimento da Figura 4.2, o sistema a ser resolvido para calculo do modelo

de velocidade e dado por ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

m1

m2

m3

m4

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦=

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

l l 0 0

0 l 0 l

0 0 l l

l 0 0 l

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

−1 ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

D1

D2

D3

D4

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦(4.35)

onde os tempos de transito estao descritos no vetor de dados com termos Di, o modelo

e representado pelo vetor mi, que contem as vagarosidades (inverso da velocidade), e a

matriz de sensibilidade (lei que rege a fısica do problema) e quadrada e contem a distancia

l entre as faces de cada celula. Neste problema o campo de velocidade e completamente

determinado pois o dispositivo de aquisicao amostra todos os numeros de onda (1l e 1

2l )

necessarios para descrever o modelo nas direcoes vertical e horizontal.

Expandindo o numero de celulas do modelo da Figura 4.2, seja nx na direcao horizon-

tal e nz na direcao vertical, se todos os pontos deste modelo se comportam como fonte e

receptor em um experimento aqui denominado modelo explosivo e; se a cada detonacao

todos os demais pontos se comportam como receptores; entao seria possıvel recuperar

53

completamente o modelo de velocidade. Todos os comprimentos de onda seriam amos-

trados e a matriz de sensilidade seria inversıvel. Essa afirmacao pode ser avaliada pela

aplicacao do Teorema da Projecao de Fatias (Projection Slice Theorem), TPF, que postula

o seguinte: a transformada de Fourier do dado em uma superfıcie de observacao e igual

a uma fatia do modelo descrito no domınio do numero de onda (Lo and Inderwiessen,

1994), desde que os raios sejam retilıneos, ou seja, o modelo de fundo deve ser constante.

Atraves do TPF, representa-se em um diagrama o numero de onda na direccao x versus

numero de onda na direcao z (kx vs. kz), as frequencias espaciais de um modelo que

sao passıveis de recuperacao em um problema inverso 2D. O diagrama kx vs. kz para o

experimento do modelo explosivo esta representado na Figura 4.3. Observa-se que todos

os comprimentos de onda sao igualmente amostrados em todas as direcoes.

Figura 4.3: Diagrama kx vs kz para um dispositivo de aquisicao do tipo

modelo explosivo. Todos os numeros de onda sao amostrados igualmente.

A Figura 4.3 representa a situacao ideal para amostragem de um modelo. Contudo, o

dispositivo de aquisicao sısmica convencional, fontes e receptores na superfıcie, com ar-

ranjo split-spread (receptores a esquerda e a direita da fonte em um arranjo 2D) apresenta

54

deficiencia de amostragem (Figura 4.4). Observa-se que os numeros de onda kz baixos,

entre 0 e 0,25 para este exemplo, nao sao amostrados. Acamamentos sedimentares hori-

zontais espessos nao sao passıveis de inversao. As duas areas em branco delimitadas na

Figura 4.4 pelos semi-cırculos, portanto nao amostradas, demonstram o quao severo e o

corte de frequencias espaciais com este dispositivo de aquisicao split-spread.

-2.0

-1.5

-1.0

-0.5

z_

k

-2 -1 0 1k_x

0

50

100

150

200

ed

utilp

ma

Figura 4.4: Diagrama kx vs kz para um dispositivo de aquisicao do tipo split-

spread. Observam-se heterogeneidades e vazios de amostragem de alguns

numeros de onda.

Os vazios de informacao em alguns comprimentos de onda aumentam o espaco nulo

da matriz de sensibilidade, fonte de instabilidades e aumento da ambiguidade na solucao

do Problema Inverso. No proximo topico apresentam-se formas de contornar este pro-

blema. Pontua-se que existem outras causas de mal condicionamento da matriz de sen-

sibilidade, entre elas a baixa iluminacao provocada pelas estruturas em subsuperfıcie.

Contudo, afirma-se que a principal causa dos espacos nulos no operador L se deve ao

dispositivo de aquisicao e a limitacoes instrumentais (receptores e fontes) que nao sao

55

capazes de registrar todas as frequencias e os numeros de onda associados.

4.5 Transformando um problema mal posto em bem

posto

Grande parte dos problemas inversos sao mal postos. Neste topico sao apresenta-

dos metodos que tendem a melhorar as condicoes de resolucao, sobretudo para reduzir

sua instabilidade, entre eles a regularizacao (processo no qual sao impostos vınculos e

restricoes que induzem a solucao (Aster et al., 2005)). De acordo com a bibliografia

consultada (Flechter, 1980; Nocedal e Wright, 2000; Aster et al., 2005) os processos

voltados para reducao da instabilidade e inducao a solucao podem atuar:

- na parametrizacao do modelo;

- na concepcao da funcao objetivo;

- na manipulacao da matriz de sensibilidade (tecnicas SVD);

Descrevem-se a seguir as diferentes formas voltadas para geracao de problemas inver-

sos bem postos.

4.5.1 Parametrizacao do modelo

Aster et al. (2005) citam que pela selecao cuidadosa da dimensao da malha (grid)

e possıvel produzir um sistema de equacoes bem condicionado. Os autores nomeiam esse

procedimento de regularizacao por discretizacao. Como citado na secao 2.7, a melhor

forma de parametrizar um modelo e por meio de celulas pequenas o suficiente para des-

crever os seus menores detalhes. Contudo, pelas razoes apresentadas nas secoes 2.7 e 4.4,

almeja-se a reducao do numero de parametros (Backus e Gilbert, 1968).

Afirma-se que as parametrizacoes, por camadas, representadas pelas equacoes 2.22,

2.23 e 2.24 constituem formas de regularizacao uteis para descricao de modelos de ve-

locidade em meios sedimentares. Essas funcoes preenchem o vazio de informacoes ge-

rado pela discretizacao por celulas. Adicionalmente, elas empregam um menor numero

de parametros pois simulam o comportamento esperado das rochas em ambientes sedi-

mentares. As parametrizacoes pelas equacoes 2.22, 2.23 e 2.24 constituem contribuicoes

56

tecnico-cientıficas desta tese.

4.5.2 Regularizacao atuando na concepcao da funcao objetivo

A funcao objetivo mais simples possıvel e alcancada atraves da equacao 4.11.

A insercao de termos adicionais na funcao objetivo e o instrumento para aplicacao de

vınculos e restricoes. Se a funcao objetivo Θ for adicionada uma funcao de restricao ou

vınculos ψ1c1(m), onde ψ1 e um escalar, obtem-se uma nova funcao

L (m, ψ1) = Θ(m) + ψ1c1(m)

Adicionando mais funcoes ci ponderadas por ψi obtem-se:

L (m, ψi) = Θ(m) +

n∑i=1

ψici(m). (4.36)

Demonstra-se que nos pontos de solucao m∗ existem escalares ψ∗i tais que

∇mL (m∗, ψ∗i ) = 0 (Nocedal e Wright, 2000). Assim, nos pontos estacionarios de

L (m∗, ψ∗i ) observa-se:

∇mΘ(m∗) = −n∑

i=1

ψ∗i∇mci(m∗). (4.37)

Os escalares ψ∗i sao os parametros de regularizacao.

Sobre as funcoes ψici citadas de forma generica no paragrafo anterior e que devem ser

construıdas as restricoes e vınculos que relacionam o processo matematico de inversao

com a realidade geologica. Uma das aplicacoes e a imposicao de uma relacao entre o

modelo mk em uma iteracao k com um modelo de referencia mre f :

c1(m) =∥∥∥Wm(m −mre f )

∥∥∥2

(4.38)

onde Wm e uma matriz MxM (onde M e o numero de parametros do modelo) que pon-

dera o vetor diferenca mk −mre f . Na tomografia cinematica para obtencao do campo de

velocidade, o modelo de referencia mre f pode ser obtido com o conhecimento do com-

portamento mais frequente das rochas:

- o meio geologico sedimentar apresenta maior conteudo de altas frequencias espaciais na

direcao vertical que nas direcoes horizontais;

- devido a compactacao, um mesmo litotipo em uma profundidade Z tera menor poro-

sidade em uma profundidade Z + dZ - desde que dZ > 0. Observadas as condicoes de

57

compactacao em meio drenado, a velocidade do meio equivalente em Z + dZ sera maior

que em Z.

Assim, o modelo de referencia toma a forma de uma funcao do tipo mre f = mre f (Z)

em situacoes onde o conhecimento de subsuperfıcie e limitado. Supondo que o modelo

de referencia e o modelo inicial, a solucao iterativa sera dada por

Δm = (WdLT L + ψWm)−1[LT WdΔD + ψWm(m − mre f )] (4.39)

onde Wd e um operador para ponderar os dados. Quando ha mais de uma funcao de

restricao a equacao 4.39 generaliza-se para:

Δm = (WdLT L +p∑

i=1

ψiWmi)−1[LT WdΔD +

p∑i=1

ψiWmi(m − mre f )] (4.40)

Na pratica a regularizacao em epıgrafe adiciona valores aos termos da matriz LT L

(equacao 4.39).

E oportuno destacar o efeito que a regularizacao causa no processo de inversao. Se o

valor de ψ for muito maior que os termos da matriz LT L, a soma (LT L+ψI) tende a ψI e

a equacao 4.11 se transforma em

Δm = (ψWm)−1LTΔD. (4.41)

O termo LTΔD e o gradiente do resıduo (ΔD). Depreende-se que o sistema tende ao

metodo de maior declividade (Aster et al., 2005) ponderado por ψ−1.

De forma contraria, se o valor ψ for muito menor que os termos da matriz LT L, nao

se observa efeito de regularizacao e o sistema tende para o metodo de Gauss-Newton,

equacao 4.22.

E desejavel que a adicao de valores a funcao objetivo, a torne cada vez mais concava.

Um funcional que possua calhas e pontos de cela desestabiliza o processo de inversao

alem de proporcionar multiplas solucoes. Assim a adicao de uma funcao que, literal-

mente, deforme o funcional, melhora a estabilidade e reduz a ambiguidade. Este fato e

observado no exemplo a seguir com um modelo representado por dois parametros, m1 e

m2, e um funcional calculado por

Θ = m21. (4.42)

Este funcional assemelha-se aquele apresentado no apendice A onde se ilustra o efeito

dos modelos cinematicamente equivalentes. A Figura 4.5 exibe o funcinal Θ da equacao

58

4.42. Observa-se que, neste exemplo, Θ dispoe-se segundo uma calha vertical com uma

serie de mınimos ao longo do eixo m2. Com este funcional a solucao sempre depende

do modelo inicial (m1o ,m2o). Para um modelo inicial (m1o ,m2o) =(2.0,2.0) a solucao e

(0.0,2.0) e; para um modelo inicial (m1o ,m2o) =(-2.0,-2.0) a solucao e (0.0,-2.0) (Figura

4.5). Neste caso ha solucoes ambıguas.

Figura 4.5: Funcao objetivo com duas variaveis m1 e m2 dada pela equacao

4.42. Ilustra-se a convergencia para o mınimo a partir de dois modelos

iniciais (m1o ,m2o) = (2.0,2.0) e (m1o ,m2o) = (-2.0,-2.0) gerando as respectivas

solucoes (0.0,2.0) e (0.0,-2.0). Neste caso ha solucoes ambıguas.

59

Se na equacao 4.42 adicionar-se um termo quadratico m22, o novo funcional Θ′ sera

Θ′ = m21 + m2

2. (4.43)

Sua forma esta representada na Figura 4.6. A solucao e unica, qualquer que seja o modelo

inicial (m1o ,m2o). Partindo dos modelos iniciais (m1o ,m2o) =(2.0,2.0) ou (m1o ,m2o) =(-2.0,-

2.0) a solucao econtra-se sempre em (m1o ,m2o) =(0.0,0.0). Nao se verifica ambiguidade

da solucao.

Figura 4.6: Funcao objetivo com duas variaveis m1 e m2 dada pela equacao

4.42. Ilustra-se a convergencia para o mınimo a partir de dois modelos

iniciais (m1o ,m2o) = (2.0,2.0) e (m1o ,m2o) = (-2.0,-2.0) gerando a mesma

solucao (0.0,0.0). Neste caso nao ha ambiguidade.

60

Neste topico foram apresentadas as formas genericas de regularizacao atraves da

concepcao da funcao objetivo.

4.5.3 Regularizacao atuando sobre a matriz de sensibilidade - Tikho-

nov de ordem zero

A pseudo-inversa calculada com a equacao 4.32 nao possui projecao sobre o espaco

nulo de L (N(L)). Contudo alguns valores singulares podem ser tao pequenos que, com

a inversao de Σ (equacao 4.30), estes podem gerar solucoes anomalamente altas. Na pior

situacao eles podem ser amplificadores de ruıdos (Aster el al. (2005)).

Segundo a condicao de Picard (Aster et al., 2005), quando o produto entre a matriz UT

(obtido pela decomposicao SV de L, secao 4.3.5) e o vetor de dados D ((UT D)i) decai mais

rapidamente que Σi ha estabilidade da solucao. A plotagem do produto (UT D)i contra

Σi constitui uma poderosa ferramenta de analise para estabilidade com a decomposicao

de valores singulares. Caso a condicao de Picard nao seja atendida, a solucao consiste

em eliminar os valores singulares extremamente baixos ao preco de um menor ajuste do

modelo. Esse procedimento e denominado regularizacao por truncamento dos valores

singulares (TSVD).

Com a regularizacao de Tikhonov (Tikhonov e Arsenin, 1977) busca-se uma solucao

que, alem de fornecer o menor resıduo, tambem apresente a menor norma. Na sua

deducao, que pode ser conduzida a partir dos multiplicadores de Lagrange (Aster et al.,

2005), elevam-se os valores singulares na busca por maior estabilidade. Na regularizacao

de ordem zero, cada valor singular Σi deve ser multiplicado por um fator de filtro (Aster

et al., 2005) fi tal que

fi =Σ2

i

Σ2i + ψ

2i

, (4.44)

onde ψ e o parametro de regularizacao.

4.5.3.1 Aplicacao de restricoes na matriz de sensibilidade e controle da norma da

perturbacao no modelo

Em problemas inversos fracamente nao lineares, (secao 4.1) as perturbacoes no

modelo relacionam-se linearmente com as perturbacoes nos dados 4.32 (Sen, 2006). A

atualizacao no modelo e calculada pela equacao 4.33 e a solucao final e obtida quando o

61

resıduo se torna menor que um valor pre-determinado de tolerancia.

Ainda para problemas fracamente nao lineares, a matriz de sensibilidade L deve ser

calculada em cada iteracao assim como sua pseudo-inversa L†. Quando o espaco nulo

de dados [N(LT )] e nao trivial e o espaco nulo de modelos [N(L)] e trivial, a solucao

por equacoes normais [(LT L)−1] e obtida pela pseudo-inversa de Moore-Penrose sao exa-

tamente iguais. Apesar da solucao ser unica, os dados nao se ajustam perfeitamente.

Quando [N(LT )] e [N(L)] sao nao-triviais, a solucao utilizando a pseudo-inversa e um

mınimo quadratico e tambem possui menor norma. Nesse caso a solucao nao e unica e os

dados nao se ajustam perfeitamente (Aster et al., 2005).

O Problema Inverso estudado nesta tese pertence a uma das situacoes citadas no ultimo

paragrafo. O processo iterativo descrito pelas equacoes 4.32 e 4.33 possui as carac-

terısticas do metodo de Gauss-Newton. A ausencia de calculo explıcito do Hessiano pode

causar inacuracia na estimativa da perturbacao do vetor modelo Δm. Consequentemente,

uma norma incorreta de Δm pode causar nao convergencia (Figura 4.7a). Neste exemplo,

em um modelo bidimensional m, com norma inadequada, nao se observa convergencia

(Figuras 4.7a e 4.7b).

62

0.0

0.2

0.4

0.6

0.8

0.01

0.21

9876543210

rebmun noitaretI

Re

sid

ua

l

Re

sid

ua

l

(a) (b)

0.0

0.2

0.4

0.6

0.8

0.01

0.21

9876543210

rebmun noitaretI

Re

sid

ua

l

Re

sid

ua

l

(c) (d)

Figura 4.7: (a) Funcao objetivo com duas variaveis m1 e m2. A linha negra tracejada

representa o caminho da solucao ieterativa. A convergencia nao e observada. (b) Resıduo

em cada iteracao observado em (a). (c) Funcao objetivo com duas variaveis m1 e m2. A

linha negra tracejada representa o caminho da solucao iterativa. A convergencia para o

mınimo e suave e direta. (d) Resıduo em cada iteracao observado em (c).

63

Portanto, um fator de estabilizacao Ω e necessario para reduzir a norma do vetor Δm.

Assim, a equacao 4.32 se torna

Δm =1

ΩL†ΔD. (4.45)

Com o mesmo problema da Figura 4.7a, a norma do vetor Δm e reduzida por um fator Ω

e a convergencia e assegurada (Figuras 4.7c e 4.7d).

O fator Ω demanda calibracao na primeira iteracao ao modo de uma busca linear

(linear search), mas na iteracoes seguintes Ω e ajustado. Se o resıduo em um passo

subsequente k + 1 e maior que no passo k, um fator de reducao f e aplicado sobre Ω−1.

Diferentemente, se em uma iteracao posterior o resıduo e menor que na iteracao anterior,

Ω−1 e aumentado segunto este fator f . O melhor valor de f empregado nos Problemas

Inversos encontrados neste trabalho e igual a 5%.

Adicionalmente, restricoes devem ser impostas na matriz de sensibilidade para incor-

porar informacoes sobre o modelo. Santos et al. (2010) apresentam a ponderacao de cada

coluna da matriz L multiplicando-a por uma constante b variando entre 0 e 1 e encerrada

na matriz diagonal B. Parametros conhecidos e, portanto restritos, tem as correspon-

dentes colunas da matriz L multiplicados por b proximas ou iguais a zero. Parametros

desconhecidos e nao restritos tem os correspondentes termos da matriz de sensibilidade

multiplicados por b igual a 1. Assim, em cada iteracao, a perturbacao do modelo e calcu-

lada por

Δm =1

Ω(LB)†ΔD. (4.46)

64

Capıtulo 5

Obtendo operadores focais

5.1 Introducao

Nesta capıtulo apresentam-se metodos dedicados a estimativa dos operadores focais

(CFP). Na secao 5.2, o metodo de Berkhout e revisitado. Nas duas ultimas secoes, 5.3

e 5.4, sao apresentadas duas formas originais e inovadoras para calculo dos operadores

CFP. Na secao 5.5 sao elencados os fatores que devem ser considerados para emprego das

difracoes no ambito da inversao tomografica.

5.2 Estimativa de operadores focais pela abordagem con-

vencional

Berkhout (1997a,b) define operador de focalizacao (Common Focus Point Operator

(CFPO)) como a funcao de Green devida a uma fonte impulsiva deflagrada em profundi-

dade e registrada na superfıcie. Thorbecke (1997) e Bolte (2003) mostram que os CFPOs

podem ser diretamente extraıdos dos dados sısmicos de reflexao. Nesta secao descreve-se

como os operadores de focalizacao sao extraıdos dos dados sısmicos, a teoria envolvida

alem das aproximacoes e premissas adotadas.

O ponto de partida para estimativa dos CFPOs e a equacao 3.13 deduzida no topico

65

3.2 e reproduzida a seguir

u(xA) = −2

∫ ∞

−∞

∫ ∞

−∞u ⊗t∂g∂z

dxdy.

Esta equacao permite calcular o campo de ondas u(xA) no ponto A (Figura 3.1)

conhecendo-se o campo de ondas a qualquer tempo ao longo da superfıcie S 0.

Na equacao 3.13, o sımbolo ⊗t siginifica correlacao cruzada traco a traco no tempo;

u sao as famılias de tiro comum, e u(xscr, xrec, t) e um campo de ondas registrado nas

coordenadas xrec ao longo da superfıcie S 0 para uma fonte na coordenada xsrc, em um

tempo t. De acordo com Berkhout (1997a), a focalizacao das famılias de tiro u(xsrc, xrec, t)

gera famılias CFP. Tomando o ponto A como um difrator, ou um ponto sobre um refletor,

a famılia CFP deste ponto e dada pela equacao 3.13, onde u(xA) e a famılia CFP do ponto

A. Ela corresponde a uma fonte virtual deflagrada no ponto A e registrada na superfıcie

S 0 (Figura 5.1). Expandindo a equacao 3.13, obtem-se a famılia CFP:

u(xA, xsrc, t) = −2

∫ ∞

−∞

∫ ∞

−∞u(xsrc, xrec, t) ⊗t

∂g(xA, xrec, t)∂z

dxrecdy. (5.1)

Aplicando a focalizacao sobre as famılias CFP obtem-se paineis migrados pre-

empilhamento

u(xA, t)mig = −2

∫ ∞

−∞

∫ ∞

−∞u(xA, xsrc, t) ⊗t

∂g(xA, xsrc, t)∂z

dxsrcdy, (5.2)

onde g(xA, xsrc, t) e a funcao de Green ou CFPO obtido quando se assumem receptores

nas mesmas estacoes das fontes.

A estimativa de um operador CFP faz uso da equacao 3.13, como o kernel de um

Problema Inverso no qual o objetivo consiste em calcular g(xA, xsrc, t) uma vez que

u(xsrc, xrec, t) e conhecido. Na pratica, devido ao Teorema de Gauss, a partir das reflexoes

observadas nos dados sısmicos, apenas uma aproximacao da funcao de Green e calculada

(pois os receptores nao estao distribuıdos em uma superfıcie horizontal infinita).

A inversao iterativa tem inıcio com a selecao de um ponto imagem A sobre uma re-

flexao. Um CFPO a priori, gk=0, deve ser estimado para propagacao entre A e a superfıcie

S 0. Atraves deste operador inicial, a equacao 5.1 e empregada para a primeira estimativa

da famılia CFP uk=0 como

uk(xA, xsrc, t) = −2

∫ ∞

−∞

∫ ∞


∂gk(xA, xrec, t)∂z

dxrecdy, (5.3)

66

onde uk(xA, xsrc, t) e a famılia CFP calculada na iteracao k-esima com a funcao de Green

gk. A equacao 5.3 e igual a equacao 5.1, mas com o objetivo neste instante de estimar g

(a funcao de Green). Em um passo subsequente, realiza-se a correlacao cruzada traco a

traco do j-esimo operador focal com a famılia CFP na k-esima iteracao. De acordo com

o Princıpio do tempo de transito equivalente, “o tempo envolvido na propagacao reversa

do campo da fonte e a reposta do ponto focal na propagacao direta sao iguais” (Berkhout,

1997a). Se o operador inicial esta correto, a correlacao Φk(xA, xsrc, t) entre o operador

gk(xA, xsrc, t) e a famılia CFP uk(xA, xsrc, t), obtido por

Φk(xA, xsrc, t) = uk(xA, xsrc, t) ⊗t gk(xA, xsrc, t). (5.4)

ira exibir um evento horizontal no lag-zero. Berkhout (1997b) nomeia Φk(xA, xsrc, t)

(equation 5.4) como painel de deslocamento diferencial no tempo (Differential Time Shift

- DTS). Rastreando o evento correspondente ao ponto imagem A no painel DTS gerara

um vetor σk, na k-esima iteracao. Se for adotado

rk = L0 − σk, (5.5)

como medida de resıduo, representando a diferenca entre o vetor lag-zero (o evento L0) e

o evento σk no painel DTS (equacao 5.5) e o minimizando iterativamente, o operador real

g(xA, xsrc, t) sera achado.

O objetivo com este procedimento e determinar a resposta cinematica do operador

CFP. Considerando

rk(xsrc, t) = δ(x − xsrc, t − tr), (5.6)

e extraindo a metade de rk(xsrc, t) (metade do resıduo estimado no passo k),

rk, 12(xsrc, t) = δ(x − xsrc, t − tr

2). (5.7)

A minimizacao do resıduo e alcancada via a atualizacao de gk(xA, xsrc, t) pela

convolucao de gk(xA, xsrc, t) e rk, 12(xsrc, t)

gk+1(xA, xsrc,t) = gk(xA, xsrc,t) ∗ rk, 12(xsrc, t). (5.8)

O processo iterativo continua ate que rk(xsrc, t) se torne uma funcao

rk(xsrc, t) = δ(x − xsrc, t). (5.9)

Sob esta condicao, o princıpio do tempo equivalente e satisfeito.

67

Figura 5.1: Representacao fısica de uma famılia CFP. Ela simula uma fonte na posicao

A com receptores espalhados na superfıcie. A famılia CFP registra a onda direta (linha

pontilhada) e a onda refletida na interface abaixo de A.

68

Apresentou-se neste topico a forma de obtencao do operador CFP originalmente

proposta por Berkhout (1997a,b) que pode ser sintetizada pelo pseudo-algoritmo a seguir:

1- escolha de um ponto imagem;

2- estimativa de um operador g0;

3- geracao da famılia CFP uk(xA, xsrc, t) com equacao 5.3;

4- calculo do painel DTS Φk com a equacao 5.4;

5- rastreamento do evento σk no painel DTS;

6- calculo do resıduo rk com equacao 5.5;

7- atualizacao do operador gk+1 com a equacao 5.8;

8- repetir os passos de (3) a (7) ate obter a condicao da equacao 5.9.

No Apendice C ilustra-se a obtencao do operador CFP a partir dos passos descritos

acima.

5.3 Estimativa de operadores focais por varredura

A necessidade de picagem em varias iteracoes e obtencao de σi para diversos pontos

imagem de uma reflexao em um universo de varios refletores/reflexoes constitui tarefa

interativa com elevado consumo de tempo. Sob o escopo desta tese apresenta-se uma

forma de estimativa de operadores CFP recomendada para areas pouco deformadas ou

com mediana complexidade estrutural. Nesta abordagem interessa-se somente pelo tempo

de transito dos operadores focais, denominados doravante tambem por g.

A essa tecnica nomeia-se Metodo de Varredura de Operadores CFP (CFPO scanning

method) e ela constitui uma contribuicao original desta tese. Ela consiste em, primeira-

mente, escolher uma funcao que represente a parcela cinematica do operador focal g(x).

Esta funcao encerra o modelo de velocidade entre o difrator e a superfıcie atraves de

coeficientes Aj presentes em g(x). Em seguida devem ser definidos os limites mınimo e

maximo para cada um dos coeficientes Aj, assim como intervalos incrementais ΔAj. Cada

conjunto de coeficientes compoe um numero m limitado de modelos mi que definira uma

funcao gi(x).

Para cada modelo mi calcula-se a famılia CFP de um ponto, seja A, em subsuperfıcie

69

u(xA, xsrc, t) de acordo com a equacao 5.3 adequadamente reproduzida a seguir

ui(xA, xsrc, t) = −2

∫ ∞

−∞

∫ ∞


∂gi(xA, xrec, t)∂z

dxrecdy.

A correlacao entre as funcoes de Green, calculadas por gi(x), e a famılia CFP correspon-

dente gera o painel Φi

Φi(xsrc, t) = gi(xA, xsrc,t) ⊗t ui(xA, xsrc, t), (5.10)

que tambem pode ser representado por:

Φi(mi, xsrc, t) = gi(xA, xsrc,t) ⊗t ui(xA, xsrc, t). (5.11)

Em um passo adicional a funcao Φi(mi, xsrc, t) e integrada ao longo de xsrc e desta

nova funcao extraem-se os valores absolutos (operador abs).

Ψ(mi, t) = abs[

∫xsrc

Φ(mi, xsrc, t)dxsrc]. (5.12)

A busca pelo maximo (max) da funcao Ψ(mi, t) no tempo t = 0 (lag 0) definira o

melhor modelo mfinal. Isso ocorre porque, caso a varredura compreenda o modelo real, o

maximo da funcao sera igual ao modelo real ou estara muito proximo a este. Somente sob

essa condicao ocorrera o alinhamento em fase (do evento referente ao difrator no ponto

A) ao longo do lag 0 gerando, portanto, o maior valor no empilhamento (equacao 5.13).

mfinal = m{max[Ψ(mi, t = 0)]}, (5.13)

onde: max e o operador para extracao do maximo de funcoes.

O fluxograma da Figura 5.2 resume o processo para obtencao de CFPO por varredura.

Nas proximas linhas apresenta-se uma aplicacao da varredura de CFPOs.

A Figura 5.3 representa uma secao de afastamento nulo (zero-offset) de um modelo

sintetico realıstico em ambiente marinho. Ha dois pontos imagens nas posicoes x 12500

m, um esta no fundo do mar (ponto 1) e o outro encontra-se na mesma coordenada x, mas

em maior profundidade abaixo do fundo do mar (ponto 2). A funcao g(x) representa o

operador

gi = (t20 +

x2

v2i

)12 + αix, (5.14)

70

Definição do modelo

m1

< mi< m

n

Cálculo de

u1(x

A,x

src,t)

Cálculo de

u2(x

A,x

src,t)

Cálculo de

un(x

A,x

src,t)

Cálculo de

Ф1(m

1,x

src,t)

Cálculo de

Ф2(m

2,x

src,t)

Cálculo de

Фn(m

n,x

src,t)

Cálculo de

Ψ(mi,t)

Modelo final:

m = m{ max [Ψ(mi,t=0)] }

Definição do modelo

m1

< mi< m

n

Cálculo de

u1(x

A,x

src,t)

Cálculo de

u2(x

A,x

src,t)

Cálculo de

un(x

A,x

src,t)

Cálculo de

Ф1(m

1,x

src,t)

Cálculo de

Ф2(m

2,x

src,t)

Cálculo de

Фn(m

n,x

src,t)

Cálculo de

Ψ(mi,t)

Modelo final:

m = m{ max [Ψ(mi,t=0)] }

Figura 5.2: Fluxograma para estimativa do operador focal referente a um difrator (difrat).

O processo tem inıcio com a definicao de um numero inicial n de modelos mi sobre os

quais sao calculados os respectivos operadores Gi. As funcoes, ui e Φi sao calculadas de

forma a compor a matriz Ψ. O modelo final e aquele apontado pelo maximo do modulo

da funcao Ψ no lag 0.

71

3.0

4.0

5.0

8000 10000 12000 14000 16000 18000 20000

Zero-offset section

3.0

4.0

5.0

8000 10000 12000 14000 16000 18000 20000

Figura 5.3: Secao de afastamento nulo exibindo dois pontos imagens em x = 12500

m. Ponto 1 no fundo do mar e ponto 2 abaixo do fundo do mar.

72

onde gi e o tempo de transito do operador, t0 e o tempo correspondente ao ponto imagem,

x e a posicao na direcao horizontal, vi e uma constante com dimensao de velocidade, e

αi e o coeficiente angular. O termo αix executa uma inclinacao da hiperbole de forma a

representar variacoes laterais suaves do campo de velocidade.

Para o ponto 1, no modelo mi, vi foi limitado a 1300 < vi < 1800, com Δvi = 100 m/s;

α = 0 simulando a nao existencia de mudanca lateral de velocidade. Assim, neste caso, vi

representa a velocidade RMS. O fluxograma da Figura 5.2 e seguido e a matriz Ψ(mi, t)

(Figura 5.4) e calculada. A amplitude maxima no painel Ψ(mi, t) para o ponto ocorre sob

a velocidade de 1500 m/s, correspondendo a coluna d’agua.

Ψ

Figura 5.4: Painel Ψ(mi, t) para o ponto 1. O eixo horizontal corresponde a velocidade

variando entre 1300 e 1800 m/s. O eixo vertical representa o lag da correlacao cruzada.

A amplitude maxima no lag zero ocorre com a velocidade de 1500 m/s.

Para o ponto 2, o modelo mi tem a velocidade vi limitada a 1400 < vi < 1800, com

Δvi = 20m/s; e α = 0, tambem nao considerando variacao lateral de velocidade na pilha

73

ΨFigura 5.5: Painel Ψ(mi, t) para o ponto 2. O eixo horizontal corresponde a velocidade

limitada entre 1400 e 1800 m/s. O eixo vertical representa o lag da correlacao cruzada.

A amplitude maxima no lag zero ocorre com a velocidade de 1520 m/s.

sedimentar rasa. A matriz ψ(mi, t) e calculada e esta ilustrada na Figura 5.5. A maxima

amplitude no painel ψ(mi, t) ocorre sob a velocidade de 1520 m/s, que e compatıvel com

a velocidade RMS para aquele nıvel.

74

5.4 Estimativa de operadores focais a partir de difracoes

Levantamentos sısmicos em areas complexas estruturalmente geram volumes de

dados sısmicos ricos em difracoes cuja origem se deve a grande curvatura dos refletores e

descontinuidades das interfaces. Nestes ambientes geologicos, a estimativa tradicional de

operadores CFP, reproduzida a partir de Thorbecke (1997) e Bolte (2003) no topico 5.2, e

difıcil de aplicar. O rastreamento de um evento σi no painel DTS e comprometido se este

nao possuir continuidade lateral devido a interferencia de difracoes. Contudo as difracoes

podem ajudar na estimativa dos CFPOs.

Considere um par fonte-receptor com afastamento Δx em um meio homogeneo (Fi-

gura 5.6). O tempo de transito duplo, na trajetoria fonte (S), difrator (A) ate o receptor (R)

e obtido pela integral de linha ao longo da curva SA (lS A), com segmentos infinitesimais

dlsi,di f ract, mais a integral de linha ao longo da curva AR (lAR), com segmentos de linha

infinitesimais dldi f ract,rq;

t(xsi , xrq) =

∫lS A

dlsi,di f f ract

v(x, y, z)+

∫lAR

dldi f f ract,rq

v(x, y, z). (5.15)

Fazendo Δx tender a zero, o angulo agudo θ (Figura 5.6) tambem tende a zero, e pode

ser facilmente deduzido que∫lS A

dlsi,di f f ract

v(x, y, z)=

∫lAR

dldi f f ract,rq

v(x, y, z). (5.16)

A equacao 5.16 representa o princıpio da reciprocidade aplicado a dados de afastamento

nulo. Entao , para um difrator, o tempo duplo em uma secao de afastamento nulo e dado

por

t(xsi = xrq) = 2

∫lS A

dlsi,di f f ract

v(x, y, z), (5.17)

e a parte cinematica da funcao de Green e a metade do tempo duplo t

g(xdi f f ract, xsi) =t(xsi = xrq)

2. (5.18)

Atraves da extracao dos tempos das difracoes nas secoes de mınimo afastamento, os CF-

POs podem ser calculados pela equacao 5.18. Isto e possıvel porque o mınimo afasta-

mento e muito menor que a profundidade do difrator. Assim, a aproximacao de que o

caminho do difrator ate a fonte e igual ao do difrator ate o receptor e valida. Para levan-

tamentos marinhos em areas com complexas estruturas em profundidade, esta tecnica e

especialmente aplicavel.

75

Figura 5.6: Par fonte (S ) Receptor (R) amostrando um difrator (A) na profundidade h

em um meio nao homogeneo.

5.5 Fatores a serem considerados na tomografia sobre

DCFPO

5.5.1 O efeito do afastamento na estimativa da velocidade

A principal premissa para estimativa dos operadores focais a partir das difracoes

(DCFPO) e que o afastamento seja muito menor que a profundidade do difrator.

Considere-se um ponto difrator localizado em (x, h), onde x=x0, e h=profundidade, em

um meio isotropico com apenas variacao de velocidade na direcao vertical v(z). O tempo

de transito duplo t para um levantamento de afastamento nulo, onde o par fonte-receptor

localiza-se em (x,0), e dado por

t(x) =

√t20+

4(x − x0)2

v2, (5.19)

onde t0 e o tempo de transito para o receptor (e a fonte) localizados em x0, acima do

difrator, e v e a velocidade RMS em (x, h).

Se erradamente considerar-se um levantamento monocanal com afastamento offset

Δx como sendo de afastamento nulo, qual sera o erro na estimativa da velocidade? A

partir das deducoes constantes no Apendice D, o mınimo tempo de transito para esse

levantamento e dado por

tmin(x) =

√4h2 + Δx2

v2, (5.20)

76

e esse tempo ocorre no receptor em

xrecmin = x0 +Δx2. (5.21)

Uma velocidade aparente, vap e calculada a partir da equacao D.12 obtida sobre uma

secao de afastamento comum (Common Offset Gather (COG)) com afastamento nao nulo.

Esta claro que a velocidade aparente nao e constante, mas varia de acordo com a relacao

entre xrec e as coordenadas do difrator (x0, h). Assim, a curva de tempo de transito quando

o afastamento nao e nulo, nao e uma hiperbole, mesmo em meios homogeneos. Contudo,

essa curva pode ser aproximada por uma funcao hiperbolica e, a partir desta, estimar-se

uma velocidade, vest, atraves de regressao.

Considere-se um levantamento hipotetico onde x0=2500 m, h=1000 m e velocidade

v=1500 m/s, no qual 501 receptores sao alinhados ao longo da superfıcie, desde a estacao

em 0 m ate 5000 m, equiespacados a 10 m. A velocidade e estimada para diferentes afas-

tamentos assumindo-se que possuem afastamento nulo. A Tabela 5.1 exibe os resıduos,

os erros na estimativa da velocidade, e o erro nas coordenadas para cada afastamento que

variou de zero a 100 m a intervalos de 10 m.

Para cada afastamento, vest e estimado por regressao da equacao 5.19, mas

empregando-se xrecmin no lugar de x0, e tmin no lugar de t0. O resıduo do tempo de transito

r em cada afastamento e calculado por

r =n∑

i=1

1

n

∣∣∣ti − tmini

∣∣∣, (5.22)

onde n e o numero de receptores.

A posicao do difrator na direcao x e estimada pela equacao 5.21 e a profundidade pelo

mınimo tempo de transito duplo tmin em cada afastamento, multiplicado pela respectiva

velocidade,vest, divida por 2.

Os resıduos aumentam com o afastamento, mas todos sao inferiores a 1 ms (menos que

o intervalo de amostragem em levantamentos sısmicos). A velocidade estimada tambem

aumenta com o afastamento, mas mesmo com o maximo afastamento de 100 m, o erro e

de apenas 0.027%. O erro na profundidade com o afastamento maximo (100 m) e igual

a 0.15%. Finalmente, o erro na coordenada x aumenta com o afastamento atingindo o

maximo de 50 m. Assim os erros na profundidade e velocidade sao negligenciaveis no

experimento apresentado. O erro na coordenada x do difrator deve, contudo, ser corrigido

ou reduzido pelo emprego da equacao 5.21.

77

Afastamento (m) Resıduo (s) εv (m/s) εh (m) εx (m)

0 0.00 0.000 0.000 0.000

10 1.58 × 10−5 0.000 0.025 5.000

20 1.92 × 10−5 0.000 0.050 10.000

30 2.47 × 10−5 0.050 0.158 15.000

40 3.79 × 10−5 0.100 0.233 20.000

50 5.43 × 10−5 0.100 0.392 25.000

60 7.55 × 10−5 0.200 0.550 30.000

70 9.38 × 10−5 0.200 0.758 35.000

80 1.28 × 10−4 0.300 0.966 40.000

90 1.62 × 10−4 0.300 1.258 45.000

100 1.85 × 10−4 0.400 1.516 50.000

Tabela 5.1: Efeito do afastamento em um meio de velocidade constante para estimativa

da velocidade e coordenadas do difrator. A velocidade real e a coordenada do difrator

sao respectivamente 1500 m/s e (2500, 1000) m. A tabela exibe os resıduos, o erro na

velocidade (εv vide texto para maiores detalhes sobre a regressao), erros na profundidade

(εh, e erros na coordenada horizontal x (εx).

78

5.5.2 Eventos que devem ser empregados na tomografia DCFPO

A tomografia sobre DCFPO falhara quando os eventos mapeados na COG nao

forem difracoes reais (geradas por pontos). As difracoes devem ser realcadas uma vez

que sua energia e mais baixa que as reflexoes (Berkovitch et al., 2009). Neste caso, as re-

flexoes especulares devem ser filtradas nas secoes COG (Fomel, 2002). Esses metodos po-

dem produzir eventos que se assemelham a difracoes, mas nao mostram as caracterısticas

de difratores pontuais. Superfıcies curvas ou mesmo segmentos de planos geram curvas

hiperbolicas nas secoes sısmicas. No entanto, seus formatos nao sao compatıveis com o

campo de velocidade real e, portanto, nao encerram informacoes cinematicas passıveis

de emprego na tomografia. Adicionamente, uma feicao geologica, como um canal ou

dobra, pode produzir uma frente de onda cilındrica que, se amostrada obliquamente a

direcao de mergulho, tera todas as caracterısticas de curva de difracao, mas seu formato

sera incompatıvel com a velocidade reinante.

Considere-se, um levantamento hipotetico com dispositivo de aquisicao paralelo a

direcao de mergulho das feicoes geologicas. Se o raio de curvatura de uma interface

dobrada e muito menor que o comprimento de onda, ou se esta interface e abruptamente

terminada por falhas, pinchout ou discordancia etc., esta comportar-se-a como difratores

(Sheriff and Geldart, 1999).

A curvatura K(x) de uma funcao 2D f (x) (Wolfram, 1999) e dada por

K(x) =

⎧⎪⎪⎨⎪⎪⎩d2 fdx2

[1 + (d fdx )2]3/2

⎫⎪⎪⎬⎪⎪⎭ . (5.23)

O raio de curvatura rc(x) e calculado pelo modulo do inverso de K(x),

rc(x) =

∣∣∣∣∣ 1

K(x)

∣∣∣∣∣ . (5.24)

Para uma interface descrita por uma funcao senoidal

f (x) = A0 sin

(2πxλ0

), (5.25)

seu raio de curvatura em cada posicao x e

rc(x) =

∣∣∣∣∣∣∣∣∣∣∣

[1 + A2

0

(2πλ0

)2 (cos 2πx

λ0

)2] 3

2

−A0

(2πλ0

)2sin

(2πxλ0

)∣∣∣∣∣∣∣∣∣∣∣. (5.26)

Pelo Princıpio de Huygens, cada ponto de uma interface, depois que ela e estimulada,

age como uma fonte pontual individual. Se em uma superfıcie suave e contınua, com

79

pontos adjacentes suficientemente proximos, e houver a coalescencia das frentes de onda,

o efeito observado e o de uma reflexao contınua. Se a interface e uma senoide (equacao

5.25), e se os picos e cavados estao afastados o suficiente na direcao vertical, ambos se

comportarao como difratores isolados, nao ocasionando interferencia construtiva entre si.

E possıvel simular esses dados relacionando a amplitude A0 e o comprimento de onda λ0

da estrutura senoidal [funcao f (x)] com a equacao 5.26, por

λ0 = α A0. (5.27)

Uma interface senoidal com amplitude (A0) de 25 m e α=10 (ou λ0=250 m) tem um

raio de curvatura muito pequeno nos picos e cavados. Um experimemto de refletor ex-

plosivo e simulado sobre esta estrutura (Figura 5.7a) empregando-se uma wavelet Ricker

de 20 Hz, em um meio com velocidade de 2000 m/s. Ao longo da superfıcie, foram

distribuıdos receptores equiespacados a cada 1 m resultando no sismograma presente na

Figura 5.7b. Em uma segunda simulacao sısmica, empregando-se a mesma geometria do

experimento anterior, mas com pontos difratores agindo como fontes somente localizadas

nos picos e cavados (Figura 5.7c). O sismograma resultante (Figura 5.7d) exibe, portanto

difracoes reais.

Uma das difracoes (a linha tracejada na Figura 5.7d) foi destacada e foi superimposta

na Figura 5.7b; a parcela cinematica da difracao real se ajusta perfeitamente as difracoes

produzidas pelos picos e cavados da superfıcie senoidal (Figura 5.7b). Esta claro que

quando α=10, picos e cavados se comportam cinematicamente como difratores.

Em um outro teste, o experimento descrito no paragrafo anterior foi repetido,

mantendo-se a amplitude da interface em 25 m mas utilizando um comprimento de onda

de 500 m (Figuras 5.8a para o refletor senoidal e 5.8c para difratores), portanto com α=20.

A Figura 5.8b ilustra o efeito do refletor explosivo senoidal com comprimento de onda de

500 m, e a Figura 5.8d exibe o sismograma para o experimento com fontes pontuais nos

picos e cavados para o mesmo comprimento de onda. As difracoes (curva branca trace-

jada na Figura 5.8b e 5.8d) nao se ajustam a curva hiperbolica correspondente ao pico da

funcao seno na Figura 5.8b. Na verdade, a curva obtida e o resultado da soma das resposta

sısmica do apice e flancos. Devido a esse somatorio, o formato da curva e mais largo que

o daquela obtida a partir de difracoes pontuais (a curva branca tracejada na Figura 5.8d)

e, portanto, nao e consistente com a velocidade de 2000 m/s.

Uma serie de experimentos foram realizados com α variando entre 10 e 20 mantendo-

80

360

380

400

420

440

460

Dep

th (m

)

0 500 1000 1500Horizontal position (m)

360

380

400

420

440

460D

epth

(m)


0.12

0.14

0.16

0.18

0.20

0.22

0.24

0.26

0.28

Tim

e (s

)


0.12

0.14

0.16

0.18

0.20

0.22

0.24

0.26

0.28

Tim

e (s

)


a) c)

d)b)

Figura 5.7: (a) Modelo de um refletor explosivo senoidal com amplitude de 25 m e

comprimento de onda igual 250 m em um meio com velocidade constante de 2000

m/s. (b) Sismograma correspondente ao modelo em (a). A linha branca tracejada ma-

peia uma difracao real no pico correspondente. Picos e cavados da curva senoidal

comportam-se como difractores; (c) Modelo com uma serie de difratores, representa-

dos por cırculos brancos, localizados nos picos e cavados da mesma curva senoidal em

(a); (d) Sismograma correspondente ao modelo em (c). A linha branca tracejada mapeia

uma difracao real.

81

0.12

0.14

0.16

0.18

0.20

0.22

0.24

0.26

0.28

Tim

e (s

)


360

380

400

420

440

460

Dep

th (m

)


360

380

400

420

440

460

Dep

th (m

)


0.12

0.14

0.16

0.18

0.20

0.22

0.24

0.26

0.28

Tim

e (s

)


a) c)

d)b)

Figura 5.8: (a) Modelo de um refletor explosivo senoidal com amplitude de 25 m e

comprimento de onda igual 500 m em um meio com velocidade constante de 2000 m/s.

(b) Sismograma correspondente ao modelo em (a). A linha branca tracejada mapeia

uma difracao real no pico correspondente. Picos e cavados da curva senoidal nao se

comportam como difractores; (c) Modelo com uma serie de difratores, representados

por cırculos negros, localizados nos picos e cavados da mesma curva senoidal em (a);

(d) Sismograma correspondente ao modelo em (c). A linha branca tracejada mapeia

uma difracao real.

82

se a frequencia dominante da fonte sob o mesmo modelo de velocidade (Tabela 5.2).

Conclui-se que, para estruturas senoidais com α maior ou igual a 16 (α crıtico, αc =16),

o comportamento de difracao pontual nao e observado. O comportamento de difracoes

pontuais ocorre se o raio de curvatura e menor que 1/16 do comprimento de onda domi-

nante da fonte. Analogamente ao criterio de resolucao de Rayleigh (Sheriff and Geldart,

1999), comprimentos de onda dominantes mais altos proporcionam αc mais baixos. As-

sim, espera-se que varias estruturas, tais como pinchouts, falhas, dobras kink, e canais

(vales incisos) possam gerar sinais como difracoes pontuais que sao uteis na tomografia

sobre CFPO.

Modelo α Comp. de onda da fonte(m) Comportamento de difrator

1 10 100 Sim

2 12 100 Sim

3 14 100 Sim

4 16 100 Nao

5 18 100 Nao

6 20 100 Nao

Tabela 5.2: Experimentos com refletores explosivos em interfaces senoidais para di-

ferentes α, empregando-se as mesmas fontes e com comprimento de onda de 100 m.

Comportam-se como difratores pontuais os picos e cavados nos modelos que apresentam

α menor que 16.

A analise descrita acima e importante para propositos de modelagem e e util na analise

de dados reais quando ha alguma sorte de conhecimento acerca das estruturas em subsu-

perfıcie. Caso nao haja nenhuma informacao de subsuperfıcie, uma das formas de se

averiguar se uma curva provem ou nao de difratores pontuais e por meio de sua per-

sistencia em forma e posicao do apice em diferentes famılias de tiro (Figura 5.9). A

parcela cinematica de uma difracao real exibira a mesma forma em diferentes famılias de

tiro (SGs) e receptores (RGs).

Matematicamente, qualquer registro da parcela cinematica de uma difracao localizada

em xA num painel de tiro comum (gS G) e obtido pela convolucao de uma funcao delta

δ(t − tsrc,A), que corresponde ao tempo de transito entre a fonte e o difrator, e o tempo de

83

Figura 5.9: Cubo de dados sısmicos; os eixos SG e RG representam famılias de tiro

comum e receptores comuns respectivamente. O eixo vertical e o tempo, crescente do topo

para base. O plano COG representa a famılia de afastamento comum para afastamento

nulo. As curvas GS G, GRG e GCOG representam o registro cinematico de uma difracao nas

famılias de tiro comum, receptor comum e afastamento (nulo) comum, respectivamente.

transito da funcao de Green entre o difrator e o receptor gA(xA, xrec,t).

gS G(xsrc, xrec, t) = δ(t − tsrc,A) ∗t gA(xA, xrec, t), (5.28)

onde o sımbolo ∗t representa a convolucao traco a traco no domınio temporal.

Landa et al. (1987) detectam difracoes utilizando correlacao e tecnicas estatısticas em

dados sısmicos em famılias de ponto medio comum (CMP), afastamento comum (COG)

e tiro comum (CSG). De forma semelhante, ao executar a correlacao cruzada entre a

resposta cinematica do difrator A na famılia de tiro 1 com a mesma resposta na famılia de

tiro 2 obtem-se

CS G(xsrc1,src2, xrec, t) = gS G(xsrc1, xrec, t)⊗t gS G(xsrc2, xrec, t) = δ(t−(tsrc1,A−tsrc2,A)). (5.29)

onde o sımbolo ⊗t significa correlacao cruzada traco a traco no domınio temporal.

84

O painel CS G(xsrc1,src2, xrec, t) exibira um evento (uma linha) horizontal no tempo

tsrc1,A-tsrc2,A quando o evento em analise for uma difracao real. Esse tempo (tsrc1,A-tsrc2,A)

constitui a diferenca de tempo entre fonte 1 e difrator A e, fonte 2 e este difrator A (equa-

tion 5.29).

Finalmente, as difracoes mapeadas devem ser comparadas com os operadores CFPs

adjacentes obtidos de forma convencional (Berkhout, 1997a,b) na mesma famılia de afas-

tamento comum. A forma de um operador DCFP deve ser semelhante a de um CFPO vi-

zinho. Mudancas abruptas observadas em operadores adjacentes representam acentuado

contraste de velocidade. Se as informacoes a priori nao corroborarem essa caracterıstica,

operadores anomalos, com raio de curvatura muito maiores que CFPOs adjacentes, devem

ser descartados.

85

Capıtulo 6

Tomografia sobre operadores focais

6.1 Introducao

Nesta secao descreve-se o metodo de estimativa do campo de velocidades e estruturas

de subsuperfıcie a partir dos operadores focais (CFPO). A parcela cinematica dos opera-

dores focais (CFPO) constitui os dados de entrada, dobsi na equacao 4.11, e esta representa

a funcao objetivo do problema contendo n registros. dcalci na mesma equacao 4.11 repre-

senta os dados calculados. Na tomografia sobre operadores focais eles representam os

tempos de transito simples entre cada ponto focal e os receptores na superfıcie. Nas roti-

nas desenvolvidas nesta tese, os tempos de transito sao calculados por rotinas de tracado

de raios cinematicos com resolucao pelo metodo de Runge-Kutta.

Empregando os mesmos passos e premissas descritos no topico 4.3.3 para solucao do

Problema Inverso pelo metodo de Gauss-Newton, a cada iteracao resolve-se a atualizacao

no modelo Δm atraves da equacao 4.22. A menos das restricoes e regularizacoes, a

resolucao pelo metodo dos mınimos quadrados empregando Gauss-Newton e a forma

basica de resolucao encontrada em Kabir (1997), Hegge (2000), e Cox (2004).

Nesta tese, o calculo das atualizacoes no modelo Δm empregou a pseudo-inversa de

Moore-Penrose utilizando a equacao 4.46, e na qual a matriz de sensibilidade e decom-

posta em valores singulares. Embora envolva um produto matricial adicional, esse metodo

permite calcular as atualizacoes no modelo reduzindo-se o numero de condicao da ma-

triz a ser invertida. Enquanto se inverte (LT L) nas abordagens tradicionais (Kabir, 1997;

86

Hegge, 2000; Cox, 2004), com a pseudo-inversa inverte-se a matriz L truncada no seu

p-esimo valor singular.

Em ambas as abordagens faz-se necessario calcular a matriz de sensibilidade L, des-

crita no topico a seguir.

6.2 Derivada em relacao a posicao do ponto focal

A matriz de sensibilidade contem as derivadas dos dados modelados em relacao

aos parametros do modelo. Entao, cada termo do vetor de dados dcalci depende do campo

de velocidade e da posicao dos pontos focais.

De acordo com Cox (2004), para um raio partindo de um ponto focal fazendo um

angulo θ com a vertical, a derivada do tempo de transito (parte cinematica da funcao de

Green), g, em relacao a coordenada x, e

dgdx=

senθv, (6.1)

e a derivada em relacao a vertical, eixo z, sob as mesmas condicoes e igual a:

dgdz=

cosθv. (6.2)

Em ambas as equacoes, v corresponde a velocidade na posicao do ponto focal.

A derivada do tempo do CFPO em relacao a velocidade v depende da forma como o

meio geologico e discretizado. A seguir, descrevem-se cada uma das formas de derivacao

em funcao das formas de discretizacao visitadas na secao 2.4.

6.3 Derivada em relacao ao campo de velocidades

O tempo de transito modelado por tracado de raios e calculado pela integral da

equacao 3.17 em relacao ao comprimento de raio infinitesimal dl. Assim, o tempo de

transito g e dado por:

g =∫

l

dlv

(6.3)

onde v e a velocidade.

87

Derivando-se a equacao 6.3 em relacao a velocidade obtem-se:

∂g∂v= −

∫l

dlv2

(6.4)

Se a velocidade for uma funcao dependente de um parametro, diga-se f , a derivada

do tempo em relacao a esse parametro sera dada por

∂g∂ f= −

∫l

dvd f

dlv2

(6.5)

que constitui a aplicaccao da regra da cadeia (Kreyszig, 1999). O que se apresenta a

seguir e a aplicacao dessa regra para as funcoes de velocidade visitadas na secao 2.4.

6.3.1 Derivada em relacao a v sob discretizacao por blocos

Supondo que o modelo tenha sido discretizado por nx blocos na horizontal e nz

blocos na vertical, o espaco de modelos de velocidade e composto por nx x nz dimensoes.

Cada bloco possui uma velocidade constante, seja v j, e dimensoes regulares dx e dz,

respectivamente na horizontal e vertical. Cada raio tracado entre o difrator e a superfıcie

atravessara diversos blocos, e a intersecao do raio com estes possui comprimento variavel

Δlk. O tempo total de cada raio e calculado de forma discreta pela soma das razoes Δlkvk

ao

longo de todo o raio com comprimento l:

g =n∑

k=1

Δlk

vk. (6.6)

Empregando-se um reduzido numero de parametros, e factıvel calcular as derivadas

perturbando-se a velocidade em cada uma das celulas. Assim computa-se a modelagem

e calcula-se a diferenca de tempos dos dados antes e depois da perturbacao e divide-se

o resultado pela perturbacao. Este procedimento e proibitivo em modelos grandes com

elevado numero de incognitas.

6.3.2 Derivada em relacao a v sob discretizacao por polinomios

Nesta forma de discretizacao, sem termos cruzados, o campo de velocidades e

representado pela equacao 2.4. Subsituindo a equacao 2.4 na equacao 6.3 e derivando o

resultado obtem-se as derivadas em relacao aos coeficientes do polinomio:

∂g∂v0

= −∫

l

[dl

(∑n

i=1 Aixi +∑m

j=1 Bjz j)2

], (6.7)

88

∂g∂Ai= −

∫l

[xi dl

(∑n

i=1 Aixi +∑m

j=1 Bjz j)2

], (6.8)

∂g∂Bj= −

∫l

[z j dl

(∑n

i=1 Aixi +∑m

j=1 Bjz j)2

]. (6.9)

6.3.3 Derivada em relacao a v sob discretizacao por funcoes exponen-

cial e logarıtmica

Substituindo a equacao 2.5 na equacao 6.3 e tomando a derivada em relacao a v0 e

B1 obtem-se:

∂g∂v0

= −∫

l

[dl

(v0 + B1e(qz))2

], (6.10)

∂g∂B1

= −∫

l

[e(qz) dl

(v0 + B1e(qz))2

], (6.11)

∂g∂q= −B1

∫l

[ze(qz) dl

(v0 + B1e(qz))2

]. (6.12)

Analogamente, a funcao logarıtmica 2.6 substituıda na equacao 6.3 e derivada em

relacao aos parametros torna-se:

∂g∂v0

= −∫

l

[dl

(v0 + B1ln(z))2

], (6.13)

∂g∂B1

= −∫

l

[ln(z)

dl(v0 + B1ln(z))2

], (6.14)

6.3.4 Derivada em relacao a v sob discretizacao por funcao de base

radial

Substituindo a equacao 2.7 na equacao 6.3, assumindo-se que os parametros σxi e

σzi sao constantes, as derivadas em relacao aos coeficientes Ai serao dadas por:

∂g∂Ai= −

∫l

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣(

(x−xi)2

σ2xi+

(z−zi)2

σ2zi

)∑n

i=1 Aie−0.5

((x−xi)2

σ2xi+

(z−zi)2

σ2zi

)⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ dl. (6.15)

89

6.3.5 Derivada em relacao a v sob discretizacao por funcoes

harmonicas

Se o campo de velocidade e representado pela equacao 2.8, a derivada em relacao

ao campo de aos parametros sera dada por

∂g∂Aj= −

∫l

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣ ei2π(k jx x+k jzz)(∑nj=1 Ajei2π(k jx x+k jzz)

)2

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦ dl. (6.16)

Analogamente, a derivada em relacao aos parametros de velocidade representados

pela equacao 2.9 e dada por:

∂g∂Ai= −

∫l

⎧⎪⎪⎪⎨⎪⎪⎪⎩ cos(kixx) + sen(kixx)[∑ni=1 Ai(cos(kixx) + sen(kixx)) +

∑mj=1 Aj(cos(k jzz) + sen(k jzz))

]2

⎫⎪⎪⎪⎬⎪⎪⎪⎭ dl,

(6.17)

e

∂g∂Aj= −

∫l

⎧⎪⎪⎪⎨⎪⎪⎪⎩cos(k jzz) + sen(k jzz)[∑n

i=1 Ai(cos(kixx) + sen(kixx)) +∑m

j=1 Aj(cos(k jzz) + sen(k jzz))]2

⎫⎪⎪⎪⎬⎪⎪⎪⎭ dl.

(6.18)

6.3.6 Derivada em relacao a v sob discretizacao por splines cubicas

(1D e 2D)

Veja-se primeiramente o caso 1D. No topico 2.4.7 foi descrita a representacao

do campo de velocidade por splines cubicas. Observou-se entao que a equacao 2.10

representa a funcao velocidade em um determinado subdomınio (x j, x j+1). De acordo com

as deducoes apresentadas, para a interpolacao por splines devem ser calculados quatro

coeficientes por subdomınio, Aj0, Aj1, Aj2 e Aj3, atraves das equacoes 2.11, 2.12, 2.13 e

2.14, respectivamente.

Analisando as quatro equacoes, 2.11, 2.12, 2.13 e 2.14, observa-se que o coeficiente

Aj0 e a propria velocidade v j em um determinado no do subdomınio. Adicionalmente,

nota-se que os coeficientes Aj1 e Aj2 sao dependentes das velocidades nos nos c j e c j+1,

e o coeficiente Aj1 e calculado de forma independente, pois e funcao somente de k j. Em

funcao dessas observacoes advoga-se que o coeficiente Aj0 e suficiente para descrever

90

as splines por subdomınio. Em outras palavras, a construcao da matriz de sensibilidade

em um modelo descrito por spline cubica 1D, e feito por meio da derivada do tempo de

transito em relacao a velocidade Aj0 no no da spline.

Assim, deduz-se que a derivada do tempo de transito em relacao a velocidade no

subdomınio j e calculada por:

∂g∂v j= −

∫l

{1

v2

[1 − (x − x j)

2 3

h2x+ (x − x j)

3 2

h3x

]}dl. (6.19)

Analogmanente, com as devidas substituicoes em funcao do acrescimo de outra di-

mensao, deduz-se a derivada do tempo em relacao a velocidade do no para o caso 2D por:

∂g∂vi, j

= −∫

l

{1

v2

[1 − (x − x j)

2(z − zi)2 3

h2xh2

z+ (x − x j)

3(z − zi)3 2

h3xh3

z

]}dl, (6.20)

onde hx corresponde ao espacamento entre os nos na direcao x e hz corresponde ao

espacamento entre os nos na direcao z.

6.3.7 Derivada em relacao a v sob discretizacao camada a camada

Este tipo de representacao demanda a existencia de camadas individualizaveis e

o conhecimento dos seus limites. Qualquer uma das funcoes apresentadas pode ser em-

pregada para caracterizar as derivadas em relacao aos parametros de velocidade em uma

camada.

Considerando que uma camada i e discretizada pelas equacoes 2.22, 2.23 ou 2.24,

somente dois parametros sao necessarios para descreve-las: uma velocidade de referencia

v0,i e um gradiente Bi. Como citado na secao 2.7, as equacoes 2.22, 2.23 e 2.24 cons-

tituem variacoes da equacao 2.18 (Kabir, 1997), onde a profundidade de referencia e

substituıda por uma interface de referencia nao necessariamente horizontal. As interfaces

de referencia das equacoes 2.22, 2.23 e 2.24 sao respectivamente o topo, o centro e base

da camada.

Considerando-se que a interface de referencia e igual ao topo do estrato, as derivadas

em relacao a v0,i e Bi sao dados pelas equacoes a seguir:

∂v∂v0,i

= 1, (6.21)

∂v∂Bi= [ztopo,i(x) − z], (6.22)

91

onde ztopo,i(x) e a profundidade do topo da camada para cada posicao x.

De forma semelhante, considerando a profundidade de referencia como o centro da

camada - discretizacao utilizando a equacao 2.23 - as derivadas parciais em relacao a v0,i

e Bi sao dadas por:

∂v∂v0,i

= 1, (6.23)

∂v∂Bi= [zc,i(x) − z], (6.24)

onde zc,i(x) e a profundidade do centro da camada para cada posicao x.

Finalmente, as derivadas em relacao a v0,i e Bi para camadas descritas pela equacao

2.24, onde a interface de referencia e a base do estrato, sao:

∂v∂v0,i

= 1, (6.25)

∂v∂Bi= [zbase,i(x) − z]. (6.26)

onde zbase,i(x) e a profundidade da base da camada para cada posicao x.

Assim, a derivada do tempo de transito em relacao aos parametros v0,i e Bi e alcancada

pelas equacoes∂g∂v0,i

= −∫

l

{1

v2

}dl, (6.27)

∂g∂v0,i

= −∫

l

{1

v2

[zre f ,i(x) − z

]}dl, (6.28)

onde zre f ,i(x) e a profundidade de referencia podendo ser ztopo,i(x), zc,i(x) ou zbase,i(x).

6.4 Consideracoes sobre a abertura dos operadores fo-

cais

A influencia da abertura dos operadores no resultado da inversao pode ser avaliada

qualitativamente a luz do Teorema da Projecao de Fatias (Projection slice theorem) resu-

mido na secao 4.4. Por este teorema descreve-se o experimento a seguir: varios pontos

difratores foram posicionados na profundidade de 1000 m e, ao longo da superfıcie, re-

ceptores foram distribuıdos com aberturas de 20000, 2000, 1000 e 200 m. O ponto sujeito

92

a analise encontra-se na posicao (0; 500) que corresponde exatamente ao centro deste

modelo finito. A Figura 6.1 resume a analise e indica quais os tipos e a orientacao das

estruturas amostradas com o respectivo dispositivo e, portanto, passıveis de estimativa no

Problema Inverso de acordo com a densidade de pontos no diagrama kx × kz.

Figura 6.1: Experimento empregando o Projection Slice Theorem para estimar o efeito

do uso de diferentes aberturas: (a) 20000m; (b) 2000 m; (c) 1000 m; e (d) 200 m.

Como o eixo horizontal kx representa as estruturas verticais e kz as horizontais, sempre

que se empregar uma abertura infinita (aproximada, neste caso pela abertura de 20000

m), todos os comprimentos de onda na direcao vertical sao amostrados. Entao, e possıvel

identificar estruturas verticais cujos comprimentos de onda vao de 0 a infinito. Mas esse

93

tipo de geometria de aquisicao nao permite a amostragem e identificacao de estruturas

horizontais ou subhorizontais. Como consequencia, o acamamento sedimentar original,

de qualquer espessura, nao e amostrado. Apenas estruturas verticais podem ser resolvidas

em problemas inversos sob o uso deste tipo de geometria.

Reduzindo o afastamento para 2000 m, a cobertura no diagrama kx × kz e reduzida.

Consequentemente, os comprimentos de onda horizontais reduzem-se, indicando que ape-

nas estruturas verticais com longos comprimentos de onda sao amostradas e, portanto, sao

passıveis de inversao.

Empregando-se um afastamento ainda menor, com 1000 m, a cobertura no espaco

kx × kz reduz-se. O diagrama kx × kz indica uma restricao mais severa na amostragem de

estruturas verticais que as demais aberturas anteriores. Entre as aberturas 20000 e 1000 m

ha uma progressiva perda de resolucao horizontal (e tambem para as estruturas oblıquas).

Sob a abertua de 1000 m, apenas estruturas verticais e subverticais com grandes compri-

mentos de onda sao amostradas e, portanto sujeitas a recuperacao atraves de inversao.

Para uma abertua de apenas 100 m, mantendo-se todos os demais parametros, ape-

nas estruturas verticais com comprimentos de onda muito grandes podem ser registradas.

Nesse caso extremo, as equacoes utilizadas para resolucao do Problema Inverso tendem a

ser linearmente dependentes, uma vez que, no caso da tomografia com kernel de tracado

de raios, estes tendem a ser paralelos, ocasionando nao unicidade de solucao. Uma con-

sequencia imediata, no caso da tomografia sobre CFPO, seria a dualidade profundidade-

velocidade. O espaco de resıduos assemelha-se ao que se observa na Figura A.3.

A analise qualitativa apresentada, apesar da limitacao oferecida pelo uso de raios re-

tilıneos, e util para estimativa de quais estruturas sao passıveis de inversao por meio da

tomografia sobre CFPO. Do exposto nos paragrafos anteriores, conclui-se que quanto

maior a abertura dos operadores CFP, mais comprimentos de onda sao amostrados,

fato que favorece o processo de inversao reduzindo progressivamente a necessidade de

regularizacoes. Adicionalmente, e completamente impossıvel amostrar camadas horizon-

tais ou subhorizontais entre os pontos focais e a superfıcie. Os operadores CFP, indivi-

dualmente, favorecem a estimatva de modelos caracterizados por gradientes laterais de

velocidade, fisicamente representados por estruturas verticais a subverticais. O conjunto

de operadores CFP distribuıdos no modelo e que permite a estimativa da velocidade entre

camadas subhorizontais e horizontais.

94

6.5 Consideracoes sobre a abordagem do problema in-

verso - inversao global ou camada a camada

A partir dos dados de tempo de transito dos operadores CFP ha duas abordagens

para resolucao da equacao 1.2 iterativamente. Uma delas, referida como inversao global

em Soares Filho (1994), considera todo o conjunto de dados D e computa a inversao para

estimativa das propriedades de todo o modelo m. Essa abordagem e seguida em Cox

(2004).

Uma segunda abordagem consiste em dividir o modelo em subdomınios e perfazer

a inversao somente com o conjunto de dados relativos aquele subdomınio. No meio se-

dimentar, o subdomınio e caracterizado por camadas e assim essa inversao e nomeada

camada a camada (layer-stripping). Neste sentido, os operadores CFPs que definem a

base da primeira camada sao empregados para estimar o campo de velocidades desta e a

geometria de sua base. Em seguida, os operadores CFPs que definem a base da segunda

camada sofrem processo identico aos da primeira. Este procedimento e repetido ate que

se atinja a base da camada mais profunda. A tecnica de inversao camada a camada e

ilustrada em Kabir (1997). Em Hegge (2000), o modelo e discretizado por camadas mas

a inversao e realizada globalmente.

Nesta secao, as inversoes global e camada a camada sao comparadas. Os dados de

entrada foram gerados sobre o modelo diapiro ilustrado na Figura 6.2. Esse modelo exibe

complexidades estruturais, mormente nos flancos do diapiro cujos mergulhos sao de 60 o.

Os operadores CFPs foram obtidos a partir da modelagem por tracado de raios de pontos

focais distribuıdos ao longo das cinco superfıcies que limitam as camadas do modelo

diapiro.

Para a inversao global, empregou-se a discretizacao por splines cubicas 2D, aquela

apontada como sendo a melhor dentre funcoes testadas na secao 2.4. Cinco nos sao uti-

lizados na direcao horizontal e cinco na vertical, totalizando 25 nos na malha regular de

spline. Para a inversao camada a camada, o modelo foi dividido em cinco camadas, e,

em cada estrato empregou-se uma discretizacao intracamada por spline 1D com cinco nos

igualmente espacados. Pela discretizacao com o campo de velocidade variando lateral-

mente, cada camada apresenta a melhor forma de discretizacao de acordo com os estudos

PST constante nas secoes 4.4 e 6.4.

95

Figura 6.2: Modelo original (modelo diapiro). Os losangos em cada interface repre-

sentam os pontos focais.

Para a inversao global, o modelo inicial e homogeneo com todas as interfaces horizon-

tais, tal como ilustrado na Figura 6.3. Nao foram empregadas restricoes na equacao de in-

versao 4.46, assim foram aplicados b = 1. O fator de escala adequado foi Ω = 0.01. Com

os parametros citados, depois de 24 iteracoes o sistema convergiu para um resıduo igual

a 0.00601 s e o modelo resultante encontra-se na Figura 6.4. O erro total (equacao 2.20)

observado com este modelo foi de 9.0808e+08 km/s e o Indice de Incoerencia Geologica

(IIG) e igual a 22.00, enquanto o modelo original apresenta IIG igual a 6.45.

Para inversao camada a camada foi empregado um campo de velocidade homogeneo

de 1500 m/s para a primeira camada e a sua base foi estimada, a priori, como estando

a 400 m de profundidade. Depois da obtencao da velocidade em cada camada, a media

destas e utilizada para estimar a velocidade a priori da camada imediatamente abaixo. A

equacao 4.46 foi empregada utilizando-se Ω variando de 0.08 a 0.1 para cada camada.

Nao foram empregadas restricoes. Foram necessarias 5, 6, 10, 6 e 17 iteracoes para as ca-

madas 1, 2, 3, 4 e 5 respectivamente para convergencia que resultou no modelo na Figura

96

Figura 6.3: Modelo inicial para inversao global. Os losangos em cada interface repre-

sentam os pontos focais no modelo a priori.

Figura 6.4: Modelo final para inversao global apos 24 iteracoes. Os losangos em cada

interface representam os pontos focais estimados com a inversao.

97

6.5. Os resıduos relativos para cada camada 1, 2, 3, 4 e 5 sao 0.00119491, 0.0010786,

0.00258387, 0.00216280 e 0.00398514 s respectivamente. O erro total (equacao 2.20) e

igual a 4.4853e+10 m/s e o IIG e 66.70.

Figura 6.5: Modelo final para inversao camada a camada. Os losangos em cada inter-

face representam os pontos focais estimados com a inversao.

A abordagem de inversao global proporcionou erros e resıduos altos. O IIG indica alta

coerencia, mas este fato se deve a um modelo caracterizado por interfaces subhorizontais e

gradientes de velocidade verticais. Essa relacao entre estruturas e gradiente de velocidade

se aproxima do acamamento original que denota alta coerencia geologica. Contudo, a alta

coerencia ocorrendo com altos erros e resıduos tornam o resultado da inversao inaceitavel.

Uma das premissas citadas na secao 4.3.3 como necessaria para solucao do Problema

Inverso pelo metodo dos gradientes nao e honrada, uma vez que o modelo a priori esta

muito longe da solucao. Esta solucao esta longe da regiao da funcao objetivo que pode

ser aproximada por uma funcao quadratica.

A inversao camada a camada forneceu erro e resıduo menores que aqueles obtidos

com a inversao global. O resıduo tendeu a crescer da camada mais rasa para a camada

mais profunda, caracterıstica esperada para a inversao camada a camada. O alto ındice de

incoerencia geologica indica que, a despeito do baixo resıduo, o modelo resultante nao e

realıstico. Trata-se de um bom modelo cinematicamente equivalente (MCE) e, portanto,

98

pode ser empregado para migracao de dados sısmicos, mas nenhuma outra informacao

dinamica pode ser extraıda dele. O seu emprego como modelo inicial em um processo de

inversao total do campo de ondas pode causar convergencia para mınimos locais.

6.6 Procedimento para estimativa do melhor MCE com

baixo IIG: tomografia sequencial

Sob condicoes iniciais similares, a abordagem camada a camada discutida no topico

anterior forneceu melhores resultados estruturais que a inversao global. A primeira e mais

estavel, uma vez que o espaco de dados e dedicado para cada camada e essa discretizacao

mostra-se mais adequada ao tipo de dispositivo de aquisicao. Para melhorar os resul-

tados da inversao, reduzindo-se o IIG e o resıduo, propoe-se um fluxo de trabalho que

envolve a inversao tomografica sequencial (ITS). Pela ITS, a abordagem de inversao ca-

mada a camada com parametrizacao das camadas por splines 1D e realizada no primeiro

estagio deste fluxo. Em seguida, os resultados estruturais e de velocidade obtidos no

primeiro estagio sao empregados em um segundo estagio de inversao utilizando-se nova

parametrizacao.

Nos dois exemplos a seguir, as inversoes tomograficas global e camada da camada sao

empregadas nesse segundo estagio, e os modelos de entrada para ambas e o resultado da

inversao camada a camada descrito na secao anterior (Figura 6.5). Para o exemplo A, a

equacao 2.22 e utilizada para representar o campo de velocidade na abordagem camada

a camada. No exemplo B, utiliza-se a mesma parametrizacao por spline 2D descrita na

secao anterior.

No exemplo A, a inversao convergiu para o modelo da Figura 6.6 depois de 22

iteracoes apresentando resıduos de 0.00110168, 0.000997521, 0.00287459, 0.00216722 e

0.00187692 s para as camadas 1, 2, 3, 4 e 5 respectivamente. O erro total (equacao 2.20) e

igual a 3.05589e+10 m/s e o IIG e 24.22. Durante as iteracoes as principais modificacoes

no modelo ocorreram no campo de velocidade.

No exemplo B, depois de 14 iteracoes, a inversao convergiu para um resıduo de

0.001712 s, um erro total (equacao 2.20) de 2.45407e+10 m/s e um IIG de 21.55. Du-

rante as iteracoes, o modelo estrutural sofreu poucas alteracoes, em oposicao ao campo

de velocidade, que foi alterado severamente (Figura 6.7).

99

Figura 6.6: Modelo final estimado com a abordagem camada a camada utilizando-se

a equacao 2.22, topo-concordante, para representar cada camada. O modelo a priori

empregado foi o resultado da inversao camada a camada presente na Figura 6.5. Os

losangos em cada interface representam os pontos focais estimados com a inversao.

Figura 6.7: Modelo final estimado com a inversao global empregando-se como modelo

a priori o resultado da inversao camada a camada presente na Figura 6.5. Os losangos

em cada interface representam os pontos focais estimados com a inversao.

100

Esse procedimento de inversao sequencial utiliza as melhores propriedades de cada

tecnica de inversao e parametrizacao no estagio mais adequado de conhecimento do mo-

delo. No primeiro estagio desse fluxo, quando nao ha qualquer conhecimento acerca

do modelo, o espaco de dados e reduzido (em concordancia com as camadas) e assim ,

estima-se um modelo cinematicamente equivalente. Neste estagio a parametrizacao das

camadas por splines 1D deve ser utilizada. Esta parametrizacao e o melhor meio de repre-

sentar o modelo de acordo com a geometria dos dados de entrada, operadores CFP (vide

secao 6.4). Assim o processo de inversao e mais estavel, nao necessita de regularizacoes,

e disponibiliza resultados acurados do ponto de vista estrutural (Figura 6.5).

Em seguida, adotando parametrizacao compatıvel com o comportamento geologico

dominante, efetua-se novo estagio de inversao. Na abordagem por camadas, foi empre-

gada a equacao 2.22 para descrever cada uma das camadas, que consiste em uma forma

de estabilizacao via parametrizacao do modelo (topico 4.5.1).

Na abordagem global tambem se observam os benefıcios da inversao sequencial. Em

funcao da estimativa do aracabouco estrutural e do comportamento cinematico, o modelo

de entrada tende a estar proximo da solucao e, por isso, as chances de convergencia para

um mınimo local sao reduzidas.

Ambos os metodos, camada a camada e global, disponibilizaram resultados com-

patıveis quando usados no segundo estagio da inversao uma vez que resıduos, erros e

Indice de Incoerencia Geologica foram reduzidos em comparacao com o primeiro estagio

da inversao. Contudo, como o IIG do modelo real e 6.45, observa-se pequena vantagem

da abordagem global para os experimentos realizados, como mostra a Tabela 6.1.

101

Modelo Erro (m/s) Resıduo (s) Convergencia IIG Ω

modelo diapiro (original) 0.00 1.00 × 10−7 1 6.45 1.00

inversao global 9.08081 × 10+11 0.00601815 22 22.00 0.10

inversao camada a camada 4.48534 × 10+10 ly 1: 0.00119491 5 66.70 0.10

discretizado por spline 1D ly 2: 0.00107860 6 0.10

ly 3: 0.00258387 10 0.10

ly 4: 0.00216280 6 0.10

ly 5: 0.00398514 17 0.08

inversao camada a camada 3.05589 × 10+10 ly 1: 0.00110168 2 24.22 0.10

sequencial discretizado no ly 2: 0.00099752 3 0.05

segundo passo por modelo ly 3: 0.00287459 5 0.01

topo-concordante ly 4: 0.00216722 10 0.10

ly 5: 0.00187692 2 0.01

inversao sequencial 2.45471 × 10+10 0.00171291 14 21.55 0.10

discretizado no segundo

passo por spline 2D global

Tabela 6.1: Comparacao do erro total, resıduo, convergencia, IIG e Omega para cada

procedimento de inversao. Na coluna dedicada a resıduo, ly significa camada. O IIG foi

calculado entre as amostras verticais 200 e 1700 dos modelos. Vide texto para detalhes

adicionais.

102

6.7 Aplicacoes

Neste topico sao apresentadas aplicacoes da tomografia baseada em operadores

CFP obrtidos com difracoes.

6.7.1 Inversao tomografica baseada em operadores CFP e DCFP -

dados sinteticos

O modelo empregado para os dois experimentos a seguir representa uma regiao

governada por tectonica compressiva. Feicoes tectonicas rasas e detachment basal sao

observados ao longo de 12 km de secao. A secao de afastamento nulo desta regiao (Figura

6.8) exibe feicoes suavemente deformadas nas extremidades esquerda e direita, e uma area

de deformacao mais intensa na regiao central.

Em ambos os experimentos foi empregada a inversao camada a camada pois os ope-

radores estao presentes ao longo das interfaces, o que permite este procedimento. Do

topo para a base do modelo, o campo de velocidade e a base da camada sao calculados

para cada camada ate que a interface mais profunda com dados de operadores CFPs dis-

ponıveis (neste exemplo, a base da camada 4) e alcancada. Nenhuma restricao ou vınculo

e aplicado, entao os termos da diagonal da matriz B sao todos iguais a 1,0 (equacao 4.46).

Os passos das iteracoes foram ajustados adequadamente com Ω igual a 2,0 na equacao

4.46.

Devido aos altos contrastes de velocidade, e porque o kernel da inversao e baseado em

tracado de raios, uma suavizacao de 31 x 31 amostras e aplicada em ambos os experimen-

tos. Um benefıcio secundario dessa suavizacao e o fato de ela evitar o emprego direto da

lei de Snell nas interfaces (Hegge, 2000).

Para o experimento 1, os espaco de dados contem os tempos de transito (as linhas

brancas tracejadas na Figura 6.8); no experimento 2, os tempos de transito dos operado-

res DCFP sao adicionados ao espaco de dados (linhas brancas contınuas na Figura 6.8).

O espaco de modelos de ambos os experimentos consiste de camadas com velocidades

constantes e coordenadas dos pontos focais (x,z). O modelo a priori contem CFPs posi-

cionados horizontalmente ao longo da base das camadas nas profundidades de 1000 m,

1600 m, 2450 m e 3300 m. Os modelos a priori de velocidades da primeira (topo) para a

quarta (mais profunda) camada foram 1500 m/s, 1500 m/s, 1800 m/s e 2400 m/s.

103

Figura 6.8: Overlay da secao de afastamento nulo com operadores CFP convencio-

nais (linhas brancas tracejadas) e operadores CFP obtidos por difracoes (linhas brancas

contınuas). O sismograma esta com tempo simples para comportar os CFPOs. Os ar-

tefatos de modelagem nas bordas (as quatro curvas hiperbolicas cujos apices estao nos

tempos de 0.65 s, 1.05 s, 1.45 s and 1.80 s) a esquerda e a direita nao sao consideradas na

analise.

104

O resıduo medio aceitavel para o tempo de transito em cada camada foi ajustado para

0.6 ms; os resıduos medios sao calculados a partir dos CFPOs observados e calculados,

equacao 4.6, e divididos pelo numero de registros. No experimento 1, nao mais que 4

iteracoes sao necessarias para encontrar a convergencia em 0.6 ms (Figura 6.9a), esti-

mando, portanto, a velocidade e a posicao dos CFPs. A Figura 6.10a exibe um overlay do

modelo de velocidade real e as interfaces (estruturas) estimadas. As velocidades e estru-

turas sao apresentadas conjuntamente na Figura 6.10b. O erro medio da velocidade por

ponto, calculado por

e =1

M

∣∣∣vobs − vcalc∣∣∣ , (6.29)

e 0.038 m/s. M e a dimensao do modelo, vobs e o modelo real e vcalc o modelo de veloci-

dade estimado.

No experimento 2, os processos de inversao, parametrizacao e tolerancia para

o resıduo empregadas sao os mesmos empregados no experimento 1, contudo as

informacoes carreadas pelas difracoes (Figura 6.8) sao tambem empregadas, correspon-

dendo, portanto a um espaco de dados mais completo. Um importante passo consiste em

mesclar os operadores convencionais CFP e os oriundos de difracoes, ordenando-os se-

gundo sua distribuicao, neste caso da esquerda para a direita. Para tanto, a posicao inicial

do difrator no eixo x e estimada a partir da coordenada x do apice do operador. Esse

procedimento permite uma melhor estimativa e interpolacao do horizonte falhado durante

a inversao.

O comportamento da convergencia do experimento 2 esta registrado na Figura 6.9b.

O emprego dos DCFPOs aumenta a resolucao estrutural (Figura 6.10c) e disponibiliza um

modelo final mais proximo do real (Figura 6.10d). O erro calculado com a equacao 6.29

e 0.035 m/s.

105

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

11019876543210

rebmun noitaretI

Re

sid

ua

l (m

s)

4 reyal3 reyal2 reyal

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.01

019876543210

rebmun noitaretI

Resid

ual

(ms)

4 reyal3 reyal2 reyal

Figura 6.9: (a) Resıduo a cada iteracao do experimento 1 para os dados sinteticos. (b)

Resıduo a cada iteracao do experimento 2 para os dados sinteticos.

106

a) c)

b) d)

1500 2000 2500

Velocity (m/s)

Figura 6.10: (a) Overlay do modelo real e dos CFPS calculados para o experimento 1.

(b) Modelo de velocidade estimado a partir da inversao dos operadores CFP convencio-

nais. (c) Overlay do modelo de velocidade real e dos CFPs calculados no experimento

2, exibindo melhor resolucao estrutural que o experimento 1. (d) Modelo de velocidade

calculado pela inversao dos CFPOs convencionais e DCFPOs conjuntamente. Todas as

figuras apresentam a mesma escala de cinza para velocidade.

107

Figura 6.11: Linha sısmica 09 e pocos em Blake Ridge.

6.7.2 Inversao tomografica baseada em operadores DCFP - dados

sısmicos monocanal

O segundo exemplo de aplicacao da tomografia sobre Operadores Focais de Difracao

(DCFPO) e realizado em um segmento de 60 km de um perfil sısmico monocanal com

direcao SW-NE e levantado em Blake Ridge (costa leste dos Estados Unidos da America)

pelo United States Geological Survey (USGS) (Taylor et al., 1999) (Figura 6.11). Esse

perfil foi escolhido por ser rico em difracoes comparado aos demais perfis presentes na

regiao.

O perfil foi adquirido com afastamento fonte-receptor (offset) constante de 10 m, e

amostragem temporal de 2 ms (Taylor et al., 1999). Para o estudo de velocidades, foi

empregado o perfil nao migrado. A linha teve seus 4136 tracos equiespacados. Com base

no espacamento fonte-receptor e nos ensaios constantes na Tabela 5.1, considerou-se o

afastamento igual a zero.

O processamento dos dados incluiu filtragem preditiva para determinacao da wavelet,

108

Figura 6.12: Linha sısmica 09 em tempo simples e as curvas de difracao mapeadas (linhas

brancas).

filtragem passa-banda (8-110 Hz) e correcao para o espalhamento geometrico.

Para realcar os eventos nao especulares, as reflexoes rasas coerentes foram atenuadas

do perfil processado empregando-se a tecnica apresentada por Cao e McMechan (2010)

para eliminacao de multiplas, aplicada, contudo, para remocao de reflexoes especulares.

Utilizando as secoes processadas com e sem reflexoes, todas as difracoes pontuais foram

identificadas. Durante o mapeamento alguns eventos foram omitidos devido aos seus

elevados raios de curvatura. A abertura horizontal dos DCFPOs na superfıcie variou

entre 1800 e 2300 m. Essas aberturas sao largas o suficiente para permitir que muitos

dos operadores DCFP apresentem redundancia com seus pares adjacentes e, portanto,

reduzam as ambiguidades velocidade-profundidade (Figura 6.12).

O modelo e discretizado por splines cubica 2D, e a inversao global foi empregada.

Esta abordagem global se mostra mais adequada, pois as difracoes encontram-se espalha-

das e nao definem interfaces ou superfıcies como na aplicacao sintetica 1. No modelo a

109

priori, a velocidade da agua de 1500 m/s (Klitgord e Schneider, 1994) e empregada para

a parcela do modelo acima do fundo do mar. Um polinomino de ordem 5

v(x,Δz) =

5∑i=0

Ci(Δz)i, (6.30)

referenciado ao fundo do mar v(x,Δz), variando de 1550 a 1750 m/s (Carmichael, 1982),

e aplicado para a pilha sedimentar simulando os efeitos combinados da compactacao e da

zona de hidratos. Os coeficientes Ci foram calculados por regressao aplicada a media das

velocidades extraıdas dos pocos 994, 995 and 997 (Paull et al., 2000; Lu and McMechan,

2002) que distam de 9 a 11 km da linha 09. Os coeficientes C0, C1, C2, C3, C4 e C5

sao respectivamente 1516.19, 1.64838e-01, -1.24340e-02, 2.15813e-03, -5.59170e-05 e

3.92534e-07; Δz e a profundidade em metros a partir do fundo do mar para a posicao x.

Trata-se de um modelo topo-concordante (vide secao 2.7) empregando um polinomio de

ordem superior.

Durante a inversao, a coluna d’agua foi restrita pela aplicacao de zeros nos elementos

correspondentes na matriz B na equacao 4.46. O passo foi ajustado com Ω−1 igual a 0.2.

A Figura 6.13 apresenta a evolucao dos modelos desde a primeira iteracao ate o modelo

final. Apos 28 iteracoes, a inversao convergiu para um resıduo medio de 0.298 ms (Figura

6.14) e disponibilizando o modelo de velocidade das Figuras 6.13h e 6.15b.

6.8 Discussao sobre as aplicacoes

Para ambos os experimentos com dados sinteticos, nao mais que quatro iteracoes fo-

ram necessarias para se atingir a convergencia em cada camada. O erro observado no

campo de velocidades no experimento 1 e superior ao observado no experimento 2 de-

vido a ausencia de informacao em 1, exatamente onde ha mudancas abruptas no modelo

(Figuras 6.8 e 6.10a).

No experimento 1, as estruturas estimadas por inversao localizadas na base da camada

2 sao muito suaves (Figure 6.10a). O experimento 2 estima melhor a posicao dos pon-

tos focais (CFPs) de forma que a posicao do refletor esta bastante proxima ao modelo

real (Figure 6.10c). Isto ocorre porque os operadores CFPs obtidos convencionalmente e

pelas difracoes adicionam informacao independente ao sistema de inversao. Em ambos

110

Fig

ura

6.1

3:

Evolu

cao

da

inver

sao

des

de

om

odel

oa

prio

riat

eo

model

ofi

nal

refe

rente

aos

dad

os

de

Bla

ke

Rid

ge.

Aes

cala

de

vel

oci

dad

eem

m/s

ea

mes

ma

par

ato

dos

os

model

os.

(a)

Model

oa

prio

ri;

(b)

quar

tait

erac

ao;

(c)

oit

ava

iter

acao

;(d

)d

ecim

ase

gunda

iter

acao

;(e

)dec

ima

sexta

iter

acao

;(f

)vig

esim

ait

erac

ao;

(g)

vig

esim

aquar

tait

erac

ao;

(h)

vig

esim

aoit

ava

iter

acao

.

111

0.25

0.30

0.35

0.40

0.45

0.50

0 5 10 15 20 25

Iteration number

Resid

ual

(ms)

Figura 6.14: Resıduos da inversao tomografica dos dados monocanal de Blake Ridge.

os experimentos observa-se que a velocidade calculada e inferior ao modelo real. Este

fato e explicado pela abordagem de inversao camada a camada. Este tipo de inversao

tende a propagar o erro encontrado nas camadas rasas para estimativa das camadas mais

profundas.

Devido a inacuracia inerente ao mapeamento das reflexoes e difracoes e a baixa

frequencia dominante dos registros sısmicos (25 Hz), todos os COGs com afastamentos

entre 0 e 60 m geraram curvas de difracao similares. A equacao 5.21 deve ser empre-

gada para minimizar os erros na estimativa das coordenadas horizontais quando forem

empregados afastamentos superiores a 60 m.

Os resultados da inversao tomografica sobre DCFPO com os dados de campo de Blake

Ridge foram comparados com os resultados da inversao tomografica de Tinivella e Lodolo

(2000). Neste trabalho os autores empregam duas linhas sısmicas multicanal localizadas

a leste da linha 09. A secao sısmica na Figura 6.12 foi convertida para profundidade e tres

reflexoes foram interpretadas dentro da pilha sedimentar: HZA, HZB e BSR (o bottom

simulating reflector) (Figuras 6.15a e 6.15b). HZA e BSR correspondem aos horizontes

2 e 3 de Tinivella e Lodolo (2000).

Ao longo do flanco da estrutura, o campo de velocidades e paralelo a geometria dos

112

Velo

city

(m/s

)

b)a)

IVIV

V

994 995 997⊗ ⊗ ⊗

I

IIII

IIIIII

IVIVI

IIII

IIIIII

V

IVIV

V

Figura 6.15: (a) Secao sısmica (linha 09) convertida para profundidade com as velocidades

calculadas durante a inversao. Abaixo do fundo do mar ha tres horizontes interpretados,

com linhas pretas, HZA, HZB e BSR (o bottom simulating reflector (BSR) sobre a zona

de gas). Os cırculos brancos constituem as posicoes dos difratores determinadas durante

a inversao dos DCFPOs. (b) Secao composta contendo o resultado da inversao sobre os

DCFPOs da linha 09 e as velocidades projetadas dos pocos 994, 995 e 997 da esquerda

(SW) para a direita (NE). Todas as interpretacoes estao repetidas em ambas as Figuras (a)

e (b). A legenda do campo de velocidade e valida somente para (b). A zona hachurada

na base representa a regiao nao invertida por ausencia de dados. No texto encontram-se as

observacoes sobre os algarismos I, II, III, IV e V.

113

refletores entre o fundo do mar e o horizonte HZA, contudo observam-se gradientes la-

terais (em I na Figura 6.15b). Abaixo da posicao apical do fundo do mar, HZA cruza

velocidades mais baixas que as dominantes naquela posicao (em II na Figura 6.15b) que

se correlacionam com uma visıvel perturbacao das reflexoes (Figura 6.15a) e com algum

deslocamento no fundo do mar (vide um pequeno monte no apice da topografia subma-

rina na Figura 6.15a). As velocidades estimadas entre o fundo do mar e o horizonte HZA

foram um pouco mais baixas do que aquelas encontradas em Tinivella e Lodolo (2000).

Abaixo do horizonte HZB, as velocidades variam entre 1650 to 1750 m/s (Figura

6.15b). Guerrin et al. (1999) interpreta essas velocidades como causadas pela presenca de

hidratos de metano cimentando os espacos intergranulares, e assim elevando os modulos

elasticos dos sedimentos hospedeiros. Abaixo do HZB, entre as posicoes 0 e 20 km

existe uma zona com altas velocidades (III na Figura 6.15b) que aparece na secao sısmica

convertida em profundidade (Figura 6.15a) como uma sismofacie homogenea delimitada

por nıveis concordantes, com reflexoes nao especulares no topo e na base. Uma zona

de alta velocidade (em IV na Figura 6.15b) e observada abaixo do apice da topografia

do fundo do mar em profundidades situadas entre HZB e BSR. A proximidade do apice

promove uma maior disponibilidade de metano e, portanto, favorece uma concentracao

maior de hidratos na zona de estabilidade da hidratos de gas (Gas hydrate Stability Zone

- GHSZ), aumentando a velocidade.

Nao foi possıvel identificar a zona de baixa velocidade associada ao apice do BSR,

indicado pela alta amplitude negativa (em V nas Figuras 6.15a e b), e associada a um inter-

valo de gas livre no poco 997. Tinivella e Lodolo (2000), empregando uma discretizacao

camada a camada, subestimam a espessura da camada com gas livre em um terco me-

nos que aquela estimada por meio de dados de perfis sısmicos verticais (Vertical Seis-

mic Profile (VSP)). Os DCFPO disponıveis e a discretizacao por splines nao apresentam

resolucao suficiente para identificar esse intervalo. De acordo com Guerrin et al. (1999)

o BSR marca o topo da ocorrencia de gas livre, nao a base da GHSZ. Entao, abaixo do

BSR gas livre deve coexistir com hidratos. Assim, a velocidade estimada e uma media

das altas velocidades dos sedimentos portadores de hidratos com as baixas velocidades

dos sedimentos contendo gas livre.

A inversao realizada por Tinivella e Lodolo (2000) esta limitada ao horizonte 4, o mais

profundo mapeado pelos autores, e situado logo abaixo do BSR. Empregando as difracoes

114

disponıveis foi possıvel investigar uma regiao mais profunda que aquela estudada em Ti-

nivella e Lodolo (2000). Nesse sentido, a tomografia DCFPO e a tomografia por reflexao

em dados multicanal sao complementares.

115

Capıtulo 7

Discussao e conclusoes

7.1 Discussao final

Nesta tese e proposto o Indice de Coerencia Geologica (IIG) que mede a con-

sistencia entre as interfaces das camadas e o gradiente de velocidade. O IIG e util para am-

bientes sedimentares e pode ser incorporado como uma ferramenta de restricao e analise

dos resultados de inversao.

Dentre as funcoes polinomial, gaussiana, harmonica, exponencial, logarıtimca e spline

2D, esta ultima e a mais indicada para representacao de modelos contınuos 2D, de acordo

com os modelos analisados.

O conceito de Modelos Cinematicamente Equivalentes (MCE) e proposto para referir-

se a ambiguidades de modelos obtidos em inversoes tomograficas cujos dados sao tempos

de transito.

Os principais metodos de inversao baseados em gradientes da funcao objetivo foram

descritos.

A influencia da abertura do dispositivo de aquisicao sob o contexto PST (Projection

slice theorem) e estudado para avaliar a necessidade de regularizacoes.

Na parte dedicada a regularizacao, mostrou-se que esta se aplica na parametrizacao

do modelo, na concepcao da funcao objetivo, ou diretamente na matriz de sensibilidade.

O esquema de inversao empregado nesta tese e baseado na formulacao de Moore-Penrose

e as restricoes sao aplicadas diretamente na matriz de sensiblidade.

O metodo convencional de estimativa de CFPO foi revisitado, e duas outras formas

de calculo sao apresentadas. Uma baseia-se na varredura de CFPOs que se aplica para

116

regioes pouco deformadas. Avaliou-se sua aplicacao em um modelo sintetico realıstico.

O outro metodo para estimativa de CFPOs baseia-se na analise de difracoes, que e util

para areas complexas.

Metodos e equacoes para construcao da matriz de sensibilidade envolvendo estruturas

e diferentes representacoes de campos de velocidade sao apresentados: por celulas ou

grid, por funcoes (polinomial, exponencial, logarıtimica, gaussiana, harmonicas, splines

cubicas 1D e 2D), e por camadas.

O efeito da abertura dos operadores CFP e estudado. Empregando o PST, ilustra-se

que a ambiguidade e reduzida a medida que a abertura dos CFPOs aumenta.

A inversao tomografica sequencial e uma proposta original deste trabalho. Ela con-

siste em usar uma parametrizacao estavel (spline 1D por camada) na primeira etapa de

inversao e empregar os resultados desta como modelo a priori para uma segunda etapa de

inversao na qual se usa uma representacao geologicamente mais realista. Um experimento

(secao 6.5) aplicou essa tecnica e duas parametrizacoes foram testadas na segunda etapa

de inversao; camada a camada com parametrizacao concordante com o topo, e global

com discretizacao por spline cubica 2D. Ambas reduziram o resıduo e o IIG em relacao a

primeira etapa, mas a sequencia envolvendo parametrizacao por spline na segunda etapa

exibiu resultados melhores.

Duas aplicacoes de inversao tomografica abordando DCFPO sao apresentadas. A

primeira envolve dados sinteticos e exibe as vantagens adquiridas quando DCFPOs sao

incorporados aos demais dados convencionais em um problema inverso. Os resultados

apontaram que a incorporacao de DCFPO aos dados de CFPO melhoram a resolucao

estrutural. Em uma segunda aplicacao, somente empregando DCFPOs sobre dados re-

ais monocanal na regiao de Blake Ridge (offshore dos Estados Unidos da America), foi

possıvel extrair o campo de velocidade da pilha sedimentar. Os resultados concordaram

com estudos previamente realizados.

7.2 Conclusoes

As conclusoes alcancadas neste trabalho sao agrupadas a seguir de acordo com a

organizacao do texto: modelo, os dados para inversao (CFPO), equacoes de inversao e

procedimento de inversao.

117

Do modelo:

1) IIG e um ındice passıvel de incorporacao na restricao (e vınculos), analise e

avaliacao dos resultados de inversao.

2) Dentre as funcoes analisadas para representacao de modelos, splines 2D e a melhor

para representacao global de modelos.

Da inversao:

3) CFPOs com grande aberturas reduzem as ambiguidades no processo de inversao.

4) CFPO pode ser obtido automaticamente, sempre que as estruturas forem suaves,

pela tecnica de varredura de CFPs.

5) CFPOs de difracoes (DCFPO) sao uteis apenas onde ha complexidades estruturais

capazes de produzir difracoes nos registros sısmicos. Metodos convencionais de obtencao

de CFPOs sao preferenciais quando nao ha complexidades.

6) Duas propostas sao apresentadas para mapear corretamente difracoes reais, a saber:

por meio da consistencia da forma das difracoes em diferentes registros de tiro ou receptor

comum (equacao 5.29), e pela comparacao das curvas de difracoes com tempo simples e

CFPOs convencionais adjacentes.

7) Para levantamentos sısmicos monocanal, a estrategia de se utilizar curvas de

difracao com menores raios de curvatura e indicada para se evitar curvas com altos α

(Tabela 5.2).

Das equacoes de inversao:

8) A equacao 4.46 de fato tem aplicacao. Ela e indicada para se evitar altos numeros

de condicao no sistema a ser resolvido.

9) O emprego de umΩ−1 adaptativo, equacao 4.46, com ajuste automatico por f regula

os passos durante o processo iterativo. O uso de um f igual a 5% e indicado uma vez que

ele exibe ampla aplicacao sob diferentes volumes de dados e complexiedade de modelos.

Do procedimento de inversao:

10) A inversao tomografica sequencial constitui uma robusta forma de se estimar mo-

delos geologicamente crıveis. A inversao camada a camada com parametrizacao estavel

(campo de velocidade variando apenas lateralmente por camada) seguida por outra etapa

parametrizada de forma mais realıstica. As parametrizacoes camada a camada concor-

dante com o topo, e global com splines 2D foram testadas na segunda etapa de inversao.

As duas se mostraram equivalentes, mas a segunda apresenta pequenas vantagens.

118

Referencias Bibliograficas

[AKI and RICHARDS, 1980] AKI, K., and P. RICHARDS, 1980, Quantitative seismo-

logy theory and methods: W.H. Freeman and Company, 1.

[AL-YAHYA, 1989] AL-YAHYA, K., 1989, Velocity analysis by iterative profile migra-

tion: Geophysics, 54, 718–729.

[ALLEN and ALLEN, 2005] ALLEN, P. A., and J. R. ALLEN, 2005, Basin analysis:

principles and applications: Blackwell, 1.

[ASTER et al., 2005] ASTER, R. C., B. BORCHERS, and C. H. THURBER, 2005, Pa-

rameter estimation and inverse problems: Elsevier Academic Press.

[BACKUS and GILBERT, 1967] BACKUS, G., and F. GILBERT, 1967, Numerical ap-

plication of a formalism for geophysical inverse problems: Geophys. J. Roy.

Astr. Soc., 13, 247–276.

[BEDNAR, 2005] BEDNAR, B., 2005, A brief history of migration: Geophysics, 70,

3MJ–20MJ.

[BERKHOUT, 1997a] BERKHOUT, A. J., 1997a, Pushing the limits of seismic imaging,

part i: Prestack migration in terms of double dynamic focusing: Geophysics,

62, 937–954.

[BERKHOUT, 1997b] ——–, 1997b, Pushing the limits of seismic imaging, part ii: Inte-

gration of prestack migration, velocity estimation, and avo analysis: Geophy-

sics, 62, 954–969.

[BERKHOUT, 1997c] ——–, 1997c, Pushing the limits of seismic imaging, part ii: Inte-

gration of prestack migration, velocity estimation, and avo analysis: Geophy-

sics, 62, 954–969.

119

[BERKHOUT and WAPENAAR, 1989] BERKHOUT, A. J., and C. P. A. WAPENAAR,

1989, One-way versions of kirchhoff integrals: Geophysics, 54, 460–467.

[BERKOVITCH et al., 2009] BERKOVITCH, A., I. BELFER, Y. HASSIN, and L. E.,

2009, Diffraction imaging by multifocusing: Geophysics, 74, 75–81.

[BERRUT and TREFETHEN, 2004] BERRUT, J., and L. TREFETHEN, 2004, Barycen-

tric lagrange interpolation: Society for Industrial and Applied Mathematics,

46, 501–517.

[BILLETTE and LAMBARE, 1998] BILLETTE, F., and G. LAMBARE, 1998, Velocity

macro-model estimation from seismic reflection data by stereotomography:

Geophysical Journal International, 135, 671–690.

[BISHOP et al., 1985] BISHOP, T., R. T. BUBE, R. CUTLER, R. T. LANGAN, P.

LOVE, J. RESNICK, R. SHUEY, D. SPINDLER, and H. WYLD, 1985, To-

mographic determination of velocity and depth in laterally varying media: Ge-

ophysics, 50, 903–923.

[BOLTE, 2003] BOLTE, J., 2003, Estimation of focusing operators using the common

focal point method: PhD thesis, Delft University of Technology.

[BOTELHO, 1986] BOTELHO, M., 1986, Modelamento sısmico na bacia do reconcavo

usando a tecnica de tracamento dos raios: PhD thesis, Universidade Federal

da Bahia.

[BUHMANN, 2003] BUHMANN, M. D., 2003, Radial basis functions: Theory and im-

plementations: Cambridge University Press.

[BULCAO, 2004] BULCAO, A., 2004, Modelagem e migracao reversa no tempo empre-

gando operadores acusticos e elasticos: PhD thesis, Universidade Federal do

Rio de Janeiro.

[CAO and MCMECHAN, 2010] CAO, J., and G. A. MCMECHAN, 2010, Multiple pre-

diction and subtraction from apparent slowness relations in 2d synthetic and

field ocean-bottom cable data: Geophysics, 75, 89–99.

[CARCIONE et al., 2008] CARCIONE, J., G. HERMAN, and A. Ten KROODE, 2008,

Seismic modeling: Geophysics, 67, 1304–1325.

120

[CARCIONE, 2007] CARCIONE, J. M., 2007, Wavefields in real media. wave propa-

gation in anisotropic, anelastic, porous and eletromagnetic media: Elsevier,

38.

[CARMICHAEL, 1982] CARMICHAEL, R. S., 1982, Handbook of physical properties

of rocks v.2: CRC Press, 2.

[CERVENY, 2001] CERVENY, V., 2001, Seismic ray theory: Cambridge University

Press.

[CLAERBOUT, 1970] CLAERBOUT, J. F., 1970, Coarse grid calculations of waves in

inhomogeneous media with application to delineation of complicated seismic

structure: Geophysics, 35, 407–418.

[CLAERBOUT, 1971] ——–, 1971, Toward a unified theory of reflector mapping: Ge-

ophysics, 36, 467–481.

[COOK, 1995] COOK, R. D., 1995, Finite element modeling for stress analysis: Wiley,

1.

[COX, 2004] COX, B., 2004, Tomographic inversion of focusing operators: PhD thesis,

Delft University of Technology.

[CUNHA, 1997] CUNHA, P. E. M., 1997, Estrategias eficientes para migracao reversa

no tempo pre-empilhamento 3d em profundidade pelo metodo das diferencas

finitas: PhD thesis, Universidade Federal da Bahia.

[CUNHA, 2005] ——–, 2005, Imageamento sısmico por propagacao de ondas no limite

de altas e baixas frequencias: PhD thesis, Universidade Federal do Rio de

Janeiro.

[DEREGOWSKI, 1990] DEREGOWSKI, S. M., 1990, Common-offset migrations and

velocity analysis: First Break, 8, 225–234.

[DIX, 1955] DIX, C., 1955, Seismic velocities from surface measurements: Geophysics,

20, 68–86.

[FEI and MCMECHAN, 2006] FEI, W., and G. A. MCMECHAN, 2006, Crp-based seis-

mic migration velocity analysis: Geophysics, 71, U21–U28.

[FLECHTER, 1980] FLECHTER, R., 1980, Practical methdods of optimization: Wiley.

121

[FOMEL, 2002] FOMEL, S., 2002, Applications of plane-wave destruction filters: Ge-

ophysics, 67, 1946–1960.

[GARDNER et al., 1974] GARDNER, G. H. F., L. W. GARDNER, and A. GREGORY,

1974, Formation velocity and density - the diagnostic basis for stratigraphic

traps: Geophysics, 39, 770–780.

[GUERRIN et al., 1999] GUERRIN, G., D. GOLDBERG, and A. MELTSER, 1999,

Characterization of in situ elastic properties of gas hydrate-bearing sediments

on the blake ridge: Journal of Geophysical Research, 104, 17781–17795.

[HADAMARD, 1902] HADAMARD, J., 1902, Sur les problemes aux derives partielles

et leur signification physique: Princeton University Bulletin, 13, 49–52.

[HAGEDOORN, 1954] HAGEDOORN, J., 1954, A process of seismic reflection inter-

pretation: Geophysical Prospecting, 2, 86–127.

[HEGGE, 1997] HEGGE, J. F., 1997, Velocity model estimation by inversion of focusing

operators: Presented at the 67th Annual International Meeting, Soc. Expl. Ge-

ophys.

[HEGGE, 2000] ——–, 2000, Seismic macromodel estimation by inversion of focusing

operators: PhD thesis, Delft University of Technology.

[HEGGE and BOLTE, 1999] HEGGE, J. F., and V. D. F. J. T. D. A. BOLTE, J.F.B., 1999,

Cfp operator estimation and inversion demonstrated on a field data set - part ii:

velocity estimation: Presented at the 69th Annual International Meeting, Soc.

Expl. Geophys.

[HEGGE et al., 1998] HEGGE, J. F., J. FOKKEMA, and A. DUIJNDAM, 1998, Using

focusing operators to estimate velocity model: Presented at the 68th Annual

International Meeting, Soc. Expl. Geophys.

[HEGGE and FOKKEMA, 1996] HEGGE, J. F., and J. T. FOKKEMA, 1996, Macro mo-

del updating by global inversion of cfp operators: Presented at the 66th Annual

International Meeting, Soc. Expl. Geophys.

[JIN and MADARIAGA, 1994] JIN, S., and R. MADARIAGA, 1994, Nonlinear velocity

inversion by a two-step monte carlo method: Geophysics, 59, 577–590.

122

[KABIR, 1997] KABIR, M., 1997, Velocity estimation of the complex subsurface using

the common focus point technology: PhD thesis, Delft University of Techno-

logy.

[KABIR and VERSCHURR, 2000] KABIR, M., and D. VERSCHURR, 2000, A cons-

trained parametric inversion for velocity analysis based on cfp technology:

Geophysics, 65, 1210–1222.

[KELLY and MARFURT, 1987] KELLY, K., and K. MARFURT, 1987, Numerical mo-

deling of seismic wave propagation: SEG, Geophysics Reprint Series.

[KELLY et al., 1976] KELLY, K., R. W. WARD, S. TREITEL, and R. M. ALFORD,

1976, Synthetic seismograms: a finite-difference approach: Geophysics, 41,

2–27.

[KHAIDUKOV et al., 2004] KHAIDUKOV, V., E. LANDA, and T. J. MOSER, 2004,

Diffraction imaging by focusing-defocusing: An outlook on seismic superre-

solution: Geophysics, 69, 1478–1490.

[KLITGORD and SCHNEIDER, 1994] KLITGORD, K., and C. SCHNEIDER, 1994,

Geophysical database of the east coast of the united states northern atlantic

margin: velocity analysis, in U.S. Geological Survey Open-File Report: 94–

192.

[KREYSZIG, 1999] KREYSZIG, E., 1999, Advanced engineering mathematics: Wiley.

[LANDA et al., 1987] LANDA, E., V. SHTIVELMAN, and B. GELCHINSKY, 1987, A

method for detection of diffracted waves on common-offset sections: Geophy-

sical Prospecting, 35, 359–373.

[LEVANDER, 1988] LEVANDER, A. R., 1988, Fourth-order finite-difference p-w seis-

mograms: Geophysics, 53, 1425–1436.

[LIU et al., 2006] LIU, J., Y. LUO, and Y. YAO, 2006, Tomographic inversion based

on cfp technology: Presented at the 76th Annual International Meeting, Soc.

Expl. Geophys.

[LO and INDERWIESEN, 1994] LO, T., and P. INDERWIESEN, 1994, Fundamentals

of seismic tomography: Society of Exploration Geophysicist.

123

[LU and MCMECHAN, 2002] LU, S., and G. MCMECHAN, 2002, Estimation of gas

hydrate and free gas saturation, concentration and distribution from seismic

data: Geophysics, 67, 582–593.

[MANSUR, 1983] MANSUR, W., 1983, A time-stepping technique to solve wave propa-

gation problems using the boundary element method: PhD thesis, University

of Southhampton.

[MARQUARDT, 1963] MARQUARDT, D., 1963, An algorithm for least-squares esti-

mation of nonlinear parameters: J. Soc. Ind. Appl. Math., 2, 601–612.

[MAVKO et al., 2009] MAVKO, G., T. MUKERJI, and J. DVORKIN, 2009, The rock

physics handbook. tools for seismic analysis of porous media: Cambridge

University Press.

[MCKAY and ABMA, 1992] MCKAY, S., and R. ABMA, 1992, Imaging and velocity

estimation with depth-focusing analysis: Geophysics, 57, 1608–1622.

[MCMECHAN, 1983] MCMECHAN, G. A., 1983, Seismic tomography in boreholes:

Geophysical Journal of the Royal Astronomical Society, 74, 601–612.

[MENKE, 1984] MENKE, W., 1984, Geophysical data analysis: discrete inverse theory:

Academic Press.

[MOORE, 1920] MOORE, E., 1920, On the reciprocal of the general algebraic matrix:

Bulletin of the American Mathematical Society, 26, 394–395.

[MOSER and HOWARD, 2008] MOSER, T., and C. HOWARD, 2008, Diffraction ima-

ging in depth: Geophysical Prospecting, 56, 627–641.

[MUSGRAVE, 1961] MUSGRAVE, A. W., 1961, Wave-front charts and three dimensi-

onal migrations: Geophysics, 26, 738–753.

[NOCEDAL and WRIGHT, 2000] NOCEDAL, J., and S. J. WRIGHT, 2000, Numerical

optimization: Springer.

[PARK, 1997] PARK, R. G., 1997, Foundations of structural geology: Chapman and

Hall.

[PAULL et al., 2000] PAULL, C., R. MATSUMOTO, P. J. WALLACE, and W. P. Dillon,

2000, Proceedings of the ocean drilling program, scientific results, 167, 1–

318.

124

[PENROSE, 1955] PENROSE, R., 1955, A generalized inverse for matrices: Procee-

dings of the Cambridge Philosophical Society, 51, 406–413.

[POLLET, 1974] POLLET, A., 1974, Simple velocity modelling and the continuous ve-

locity section: Presented at the 44th Annual International Meeting, Soc. Expl.

Geophys.

[RIEBER, 1936] RIEBER, F., 1936, A new reflection system with controlled directional

sensitivity: Geophysics, 1, 97–106.

[ROBERTSSON et al., 2007] ROBERTSSON, J. O. A., B. BEDNAR, J. BLANCH, and

D. J. KOSTOV, c. MANEN, 2007, Introduction to the supplement on seis-

mic modeling with applications to acquisition, processing, and interpretation:

Geophysics, 72, 1–4.

[SANTOS et al., 2010] SANTOS, L. A., D. SOARES FILHO, A. BULCAO, C. THEO-

DORO, W. MANSUR, and G. MCMECHAN, 2010, Tomographic inversion

over diffraction cfp operators in structurally complex areas: Presented at the

80th Annual International Meeting, Soc. Expl. Geophys.

[SANTOS and FIGUEIRO, 2006] SANTOS, R. H. M., and W. M. FIGUEIRO, 2006,

Modelagem acustica bidimensional usando diferentes parametrizacoes de

campos de velocidades: Revista Brasileira de Geofısica, 24, 103–115.

[SCHNEIDER, 1971] SCHNEIDER, W. A., 1971, Developments in seismic data proces-

sing and analysis (1968-1970): Geophysics, 36, 1043–1073.

[SCHNEIDER, 1978] ——–, 1978, Integral formulation for migration in two and three

dimensions: Geophysics, 43, 49–76.

[SEN, 2006] SEN, M. K., 2006, Seismic inversion: Society of Petroleum Engineers.

[SHERIFF, 2002] SHERIFF, R. E., 2002, Encyclopedic dictionary of exploration ge-

ophysics: Society of Exploration Geophysicist, 1.

[SHERIFF and GELDART, 1999] SHERIFF, R. E., and L. P. GELDART, 1999, Explo-

ration seismology: Cambridge University Press.

[SHEWCHUK, 1994] SHEWCHUK, J. R., 1994, An introduction to the

conjugate gradient method without the agonizing pain: http :

//www.cs.cmu.edu/ jrs/ jrspapers.html.

125

[SIDDIQUI, 2004] SIDDIQUI, K., 2004, Velocity model building methodology and

psdm in deep water gulf of mexico: A case history: Presented at the 74th

Annual International Meeting, Soc. Expl. Geophys.

[SKINNER and PORTER, 1987] SKINNER, B. J., and S. C. PORTER, 1987, Physical

geology: Wiley.

[SOARES FILHO, 1994] SOARES FILHO, D. M., 1994, Inversao dos tempos das ondas

refletidas: combinacoes de tomografia, migracao e tecnicas tipo dix: PhD

thesis, Universidade Federal da Bahia.

[SOARES FILHO and ECKHARDT, 1997a] SOARES FILHO, D. M., and B. R. S.

ECKHARDT, W., 1997a, Seismic crosswell tomography: A gauss-newton

type method based on b-splines for velocity field parametrization. 1st part:

Fundamentals: Presented at the 5o Congresso Internacional da SBGf, Socie-

dade Brasileira de Geofısica.

[SOARES FILHO and ECKHARDT, 1997b] ——–, 1997b, Seismic crosswell tomo-

graphy: A gauss-newton type method based on b-splines for velocity field

parametrization. 2st part: Application: Presented at the 5o Congresso Interna-

cional da SBGf, Sociedade Brasileira de Geofısica.

[STOLT, 1978] STOLT, R. H., 1978, Migration by fourier transform: Geophysics, 43,

23–48.

[STORK, 1992] STORK, C., 1992, Reflection tomography in the postmigrated domain:

Geophysics, 57, 680–692.

[SWORD, 1987] SWORD, C. H., 1987, Tomographic determination of interval veloci-

ties from reflection seismic data: the method of controlled directional recep-

tion: PhD thesis, Stanford University.

[SYMES and CARAZZONE, 1991] SYMES, W. W., and J. J. CARAZZONE, 1991, Ve-

locity inversion by differential semblance optimization: Geophysics, 56, 654–

663.

[TARANTOLA, 1984] TARANTOLA, A., 1984, Inversion of seismic reflection data in

the acoustic approximation: Geophysics, 49, 1259–1266.

126

[TARANTOLA, 2002] ——–, 2002, Inverse problem theory: methods for data fitting

and model parameter estimation: Elsevier.

[TARANTOLA, 2005] ——–, 2005, Inverse problem theory and methods for model pa-

rameter estimation: SIAM.

[TAYLOR et al., 1999] TAYLOR, M., W. DILLON, C. ANTON, and W. DANFORTH,

1999, Seismic-reflection surveys of the blake ridge, r/v cape hatteras 1992 and

1995; data acquisition, navigation and processing, in U.S. Geological Survey

Open-File Report 99-372.

[THOMAS, 2000] THOMAS, J. E., 2000, Velocidades sısmicas: Petrobras.

[THORBECKE, 1997] THORBECKE, J. W., 1997, Common focus point technology:

PhD thesis, Delft University of Technology.

[TIKHONOV and ARSENIN, 1977] TIKHONOV, A., and V. ARSENIN, 1977, Soluti-

ons of ill-posed problems: V.H. Winstons and Sons.

[TINIVELLA and LODELO, 2000] TINIVELLA, U., and E. LODELO, 2000, The blake

ridge bottom-simulating reflector transect: tomographic velocity field and the-

oretical model to estimate methane hydrate quantities, in Proceedings of the

ocean drilling program, scientific results,, 273–318.

[TRIER, 1990] TRIER, J. A., 1990, Tomographic determination of structural velocities

from depth-migrated seismic data: PhD thesis, Stanford University.

[VIRIEUX, 1986] VIRIEUX, J., 1986, P-sv wave propagation in heterogenous media:

Velocity-stress finite-difference method: Geophysics, 51, 889–901.

[WAPENAAR et al., 1989] WAPENAAR, C. P. A., G. L. Peels, V. Budejicky, and A. J.

Berkhout, 1989, Inverse extrapolation of primary seismic waves: Geophysics,

54, 853–863.

[WOLFRAM, 1999] WOLFRAM, S., 1999, Mathematica book: Cambridge Academic

Press.

[WOODWARD et al., 2008] WOODWARD, M. J., D. NICHOLS, O. ZDRAVEVA, P.

WHITFIELD, and T. Johns, 2008, A decade of tomography: Geophysics, 73,

VE5–VE11.

127

[ZHANG and MCMECHAN, 2011] ZHANG, Q., and G. A. MCMECHAN, 2011, Di-

rect vector-field method to obtain angle-domain common-image gathers

from isotropic acoustic and elastic reverse-time migration: Geophysics, 76,

WB136–WB149.

128

Apendice A

Experimento sismo 1

O objetivo deste estudo e obter o modelo de velocidades a partir de dados de tempos de

transito por meio da varredura do espaco de resıduos (diferenca entre os dados observados

e os dados reais). Aquele modelo que apresentar o menor resıduo constituira a solucao

almejada.

Para este experimento foi criado um modelo simples cujo campo de velocidade obe-

dece a equacao A.1.

v(x, z) = 1500.0 + 5x, (A.1)

ou seja, trata-se de um modelo de velocidade variante lateralmente, que obedece a equacao

de uma reta (v = a + bx) cujos coeficientes linear e angular sao respectivamente iguais a

a = 1500 e b = 5 (Figura A.1). Portanto, somente dois parametros definem o espaco de

modelos m.

129

Figura A.1: Campo de velocidade sobre o qual foi realizado o experimento sismo 1. Um

campo de velocidade que obedece a equacao A.1. As linhas pretas representam os raios

emanados da fonte na posicao (x=250, z=500).

130

Figura A.2: Tempos de transito registrado ao longo da superfıcie para uma fonte detonada

em (x=250, z=500) no campo de velocidade da Figura A.1.

Uma fonte pontual na posicao (x=250, z=500) e detonada e ao longo da superfıcie

sao distribuıdos receptores em diferentes posicoes X. Nestes receptores sao registrados os

tempos de chegada da onda direta como mostra a Figura A.2.

131

0.0003

0.0003

0.0003

0.0006

0.0006

0.0009

0.0009

0.0012

0.0015

0.0018

0.0018

0.0021

0.0021

0.0027

0.0027

0.0033

0.0039

0.00450.006

0

2

4

6

8

10

b (

m/s

/m)

1000 1200 1400 1600 1800

a (m/s)

Figura A.3: Espaco de resıduos com parametros a e b respectivamente na horizontal e

vertical. As curvas de nıvel indicam o valor do resıduo.

O experimento pelo metodo de varredura consistiu em construir diversos modelos

com parametros a e b variando respectivamente de 1000 a 1700, e de 0 a 10, calculando

os resıduos sobre os tempos de transito de cada um. O produto deste ensaio esta na figura

A.3.

132

Observa-se que a coordenada (a,b)=(1500,5) na figura A.3 constitui um mınimo do

espaco de resıduos, como esperado, uma vez que estes dois parametros correspondem ao

modelo real. Todavia, e possıvel observar uma serie de mınimos que satisfazem os tempos

de transito registrados na superfıcie. Na verdade, neste simples ensaio, ha uma infinidade

de solucoes, ou modelos, que satisfazem os dados da Figura A.2. Todos os modelos que

seguem a diagonal com resıduo proximo a zero da Figura A.3 sao cinematicamente equi-

valentes. Essa infinidade de solucoes e propria da ambiguidade inerente aos problemas

inversos. Somente a insercao de restricoes e vınculos permite a reducao da ambiguidade.

133

Apendice B

Experimento discretizacao

Neste ensaio o campo de velocidades NMO de uma bacia sedimentar marıtima bra-

sileira composto por 150 celulas na horizontal e 150 celulas na vertical (Figura B.1) foi

representado por algumas das funcoes apresentadas no capıtulo 2.4. As funcoes expo-

nencial e logarıtmica nao se mostraram adequadas para representacao do modelo de ve-

locidade escolhido devido a sua complexidade. Uma combinacao dessas funcoes com as

demais pode disponibilizar melhores resultados, mas estes nao foram testados no ambito

desse trabalho. A representacao por malha nao regular nao foi desenvolvida nesta tese.

Assim, o experimento de discretizacao apresentado a seguir restringiu-se as funcoes

polinomial, de base radial gaussiana, harmonicas e por splines cubicas em duas dimensoes

para descrever o modelo 1 mostrado na Figura B.1. Para as funcoes polinomial, gaussiana

e harmonicas foi empregada uma rotina que utiliza o metodo dos mınimos quadrados para

o calculo dos coeficientes adequados de cada funcao. Para essas funcoes, o experimento

constitui um pequeno ensaio de inversao linear de ajuste de superfıcies, o campo de ve-

locidades NMO. Para a funcao spline foi utilizada uma rotina propria desenvolvida nesta

tese construıda com as equacoes 2.10 a 2.17. Em todos os experimentos foi calculado o

erro quadratico, e, utilizando a equacao 2.19.

A representacao por celulas expoe todos os detalhes do campo de velocidades (Fi-

gura B.1). Contudo, encerra elevado numero de parametros - um total de 22500 celulas

(nx=150 e nz=150).

134

Figura B.1: Campo de de velocidade real - modelo 1.

B.1 Funcao polinomial

Em todas as formas de representacao almejou-se utilizar o mesmo numero de

parametros igual a 400, contudo com os polinomios nao foi possıvel empregar mais que

15 parametros (n=m=7 na equacao 2.4). Ao elevar a ordem dos polinomios, a equacao

2.4 torna-se instavel para inversao e consecutiva obtencao dos coeficientes da equacao.

Observa-se, portanto, que existe um limite associado a ordem do polinomio que permite

melhor representar um modelo. Ultrapassando-se este limite a representacao se afasta do

modelo original. Esta instabilidade presente nas interpolacoes polinomiais de ordem su-

perior e denominada fenomeno de Runge (Berrut e Trefethen, 2003). A Figura B.2 ilustra

o resultado alcancado com um polinomio de ordem 7 para cada dimensao x e z alem de

um coeficiente independente. Utilizando-se polinomios de ordem 8 para cada dimensao

obtem-se resultado nao veraz (Figura B.3).

135

Figura B.2: Campo de de velocidade do modelo 1 representado por polinomio de ordem

7 nas dimensoes x e z.

Figura B.3: Campo de de velocidade do modelo 1 representado por polinomio de ordem

8 nas dimensoes x e z.

136

O erro quadratico obtido com polinomio de ordem 7, como esperado, foi inferior ao

calculado com o polinomio de ordem 8. O primeiro apresentou erro igual a 2,701E+04 e

o segundo 1,574E+05.

B.2 Funcao gaussiana 2D

A representacao do modelo 1 por de uma funcao gaussiana em duas dimensoes,

equacao 2.7, empregando σ igual a 4,87 (2σ2xi= 2σ2

xi= 95), e utilizando 400 parametros

distribuıdos regularmente ao longo da malha, originalmente contendo 22500 celulas, le-

vou ao modelo da figura B.4. Observa-se que o resultado alcancado constitui uma versao

bastante suavizada do modelo original, exibindo, contudo, as grandes tendencias estrutu-

rais de distribuicao do campo de velocidades.

Figura B.4: Campo de de velocidade do modelo 1 representado por funcao gaussiana em

2D (vide texto para detalhes de parametrizacao).

Durante os ensaios com este tipo de representacao observou-se que o parametro σ

deve ser testado, caso contrario os resultados alcancados com funcao gaussiana serao

pouco realısticos. Uma vez que σ esta relacionado a abertura da gaussiana, valores muito

baixos favorecem a descricao de estruturas pequenas enquanto valores mais altos favore-

137

cem a representacao de estruturas maiores. Para este modelo1, empregando 400 amostras

regularmente espacadas, os valores de 2σ2xi

proximos a 2/3 da dimensao do modelo, fo-

ram os que melhor responderam, de forma a apontar tendencia geral, portanto em baixa

freqencia espacial, do campo de velocidade.

O erro quadratico obtido com a funcao gaussiana foi igual a 2,150E+04, inferior aos

observados com as funcoes polinomiais.

B.3 Funcoes harmonicas

A filragem de altas frequencias nas direcoes x e z limitando-as aos 20 primeiros

harmonicos (suavizacao de 25 amostras 0-0-10-35) forneceu o resultado constante na

Figura B.5. Claramente a filtragem, compatıvel com o emprego de 20 harmonicos em

cada direcao (perfazendo um total de 400 coeficientes, n=m=20 na equacao 2.8), des-

tacou o indesejavel carater oscilatorio. Embora, visualmente, o resultado obtido com a

serie de Fourier restrita aos 20 primeiros harmonicos - ainda que se tenha tido cautela

para se evitar efeito de Gibbs - guarde os principais aspectos estruturais do modelo 1, o

erro quadratico obtido com este tipo de representacao foi igual a 6,351E+04, superior aos

anteriores.

Utilizando o simples somatorio de funcoes seno e cosseno de acordo com a equacao

2.9, limitado a 200 coeficientes por direcao, obtemos o resultado contido na Figura B.6

O erro quadratico calculado para este tipo de representacao foi igual a 2,587E+04,

inferior ao de todas as equacoes anteriores.

B.4 Funcao spline cubica 2D

A representacao do modelo 1 por splines empregou 400 parametros. Os pontos

de amostragem, regularmente distribuıdos, representam os domınios das funcoes spline

(Kreyszig, 1999). O resultado da representacao por splines esta ilustrado na figura B.7.

Dentre todas as funcoes testadas, o uso de splines cubicas 2D apresentou o menor erro

quadratico: 5,465E+03.

138

Figura B.5: Campo de de velocidade do modelo 1 representado por serie de Fourier

restrita aos 20 primeiros harmonicos por direcao.

Figura B.6: Campo de de velocidade do modelo 1 representado por soma de senos e

cossenos restrita a 200 coeficientes por direcao.

139

Figura B.7: Campo de de velocidade do modelo 1 representado por spline cubica 2D

restrita a 200 coeficientes em cada direcao.

B.5 Resumo dos resultados

A tabela B.1 resume os erros quadraticos de cada uma das funcoes testadas sob as

mesmas condicoes possıveis. Nao ha duvidas que as funcoes por domınio, tal como as

splines no experimento em escopo, apresentam melhores resultados, em funcao do menor

erro quadratico. Essa assertiva pode ser generalizada para demais tipos de modelos que

apresentem continuidade dos gradientes e das derivadas segunda.

A segunda melhor opcao de discretizacao, em termos de erro quadratico, embora te-

nha disponibilizado um modelo bastante suavizado, foi a funcao gaussiana. Destaca-se,

contudo, que houve uma pesquisa previa para determinacao do melhor parametro σ (igual

a 4,87).

A terceira melhor forma de representacao foi atraves da soma de senos e cossenos.

Em relacao ao erro quadratico, a quarta melhor forma de representacao foi atraves de po-

linomios de ordem 7. Contudo, o reduzido numero de parametros permitido e necessario

para os polinomios os torna atrativos para calculo de modelos em baixa freqencia espacial.

Os ensaios apontaram que ha um numero otimo de ordem polinomial (neste caso igual a

140

Funcao Erro quadratico Numero de parametros

polinomial de ordem 8 1,574E+05 17

polinomial de ordem 7 2,702E+04 15

gaussiana 2,150E+04 400

serie de Fourier 6,352E+04 equivalente a 400

spline cubica 2D 5,465E+03 400

Tabela B.1: Tabela comparativa dos erros quadraticos de cada funcao na representacao de

modelo 1.

ordem 7) que deve ser pesquisada previamente. A adicao de um grau na representacao

polinomial elevou o erro relativo em uma ordem de grandeza (tabela B.1).

Finalmente a representacao do modelo 1 por serie de Fourier truncada de forma equi-

valente a 20 termos apresentou as oscilacoes inerentes as funcoes harmonicas. Este fato

elevou o erro quadratico deste tipo de funcao tornando-o inadequado para caracterizar o

modelo 1.

141

Apendice C

Obtendo operadores CFP pela

abordagem convencional

Neste Apendice os passos listados no final da secao 5.2, listados abaixo, sao seguidos

para estimativa do operador CFP em um dado sintetico de acordo com a metodologia

apresentada por Berkhout (1997a, 1997b).

1- escolha de um ponto imagem;

2- estimativa de um operador g0;

3- geracao da famılia CFP uk(xA, xsrc, t) com equacao 5.3;

4- calculo do painel DTS Φk com a equacao 5.4;

5- rastreamento do evento σk no painel DTS;

6- calculo do resıduo rk com equacao 5.5;

7- atualizacao do operador gk+1 com a equacao 5.8;

8- repetir os passos de (3) a (7) ate obter a condicao da equacao 5.9.

142

1) Escolha do ponto imagem. A Figura C.1 apresenta uma secao de afastamento nulo

onde se mostra o ponto I em que se deseja estimar o operador CFP, destacado pelo cırculo

vermelho.

3.0

4.0

5.0

8000 10000 12000 14000 16000 18000 20000

Figura C.1: Secao de afastamento nulo destacando o ponto I (cırculo vermelho) sobre o

qual se deseja estimar o operador focal (CFP).

143

2) Estimativa de um operador g0. Diversas sao as formas de obtencao de operadores focais

a priori. A mais direta e, naturalmente obtida durante o processamento de dados sısmicos,

e atraves da velocidade de empilhamento ou NMO. Aqui apresenta-se um metodo alter-

nativo que consiste em mapear o evento em analise no painel de tiro comum (shot gather)

em tempo simples (Figura C.2). Na Figura C.2, o evento mapeado, e representado pela

linha tracejada verde, e a primeira estimativa do operador focal g0. Essa abordagem de

mapeamento no shot gather com tempo simples tem a vantagem de fornecer um oper-

dor CFP a priori que engloba complexidades difıceis de serem reproduzidas por funcoes

matematicas hiperbolicas. Uma desvantagem, e a necessidade de mapeamento, fato que

diminui a automacao.

Figura C.2: Painel de tiro comum (shot gather) em tempo simples onde se estima o

operador a inicial g0 destacado pela linha tracejada verde.

144

3) Geracao da famılia CFP uk(xA, xsrc, t) com a equacao 5.3. Empregando a equacao 5.3

obtem-se a famılia CFP, ou cfp gather que esta ilustrado na Figura C.3.

Figura C.3: Famılia CFP para o ponto imagem I. A linha vermelha representa o operador

g0.

145

4) Calculo do painel DTS (Φk) com a equacao 5.4. A equacao 5.4 realiza a correlacao

traco a traco entre o operador g0 e o cfp gather na Figura C.3. Esta operacao gera o painel

DTS constante na Figura C.4.

5) Rastreamento do evento σk no painel DTS. O evento relacionado ao ponto I e mapeado

e constitui a variavel σk destacada pela linha tracejada amarela.

6) Calculo do resıduo rk com a equacao 5.5. O resıduo e obtido com a diferenca entre σk

(curva amarela tracejada na Figura C.4) e o lag zero para cada traco do painel na Figura

C.4.

Figura C.4: Painel DTS na primeira iteracao. A linha tracejada amarela representa o

evento relacionado ao ponto I. Observa-se que ele nao se encontra alinhado ao longo do

lag zero.

146

7) Atualizacao do operador gk+1 com a equacao 5.8. Com o novo operador, os passos de

(3) a (6) foram seguidos para obtencao do DTS.

Em duas outras iteracoes o evento relacionado ao ponto I horizontalizou-se ao longo do

lagzero, como pode ser observado na Figura C.5. Nesta situacao obteve-se o operador

correto.

Figura C.5: Painel DTS na terceira iteracao. Observa-se que o evento em estudo esta

alinhado ao longo do lag zero.

147

Apendice D

Posicao x do tempo mınimo e velocidade

aparente a partir de difracoes

registradas em secoes de afastamento

nao nulo

Nesta secao, as equacoes de Normal Move Out (NMO) (Sheriff e Geldart, 1999)

adaptadas para difracoes, e Landa (1987), para calculo da posicao x do tempo mınimo

duplo de uma difracao em uma secao de afastamento comum (Common Offset Gather

- COG), com afastamento maior que zero. Os resultados constituem a base para as

deducoes da velocidade aparente vap que ilustra o comportamento nao hiperbolico das

curvas de transito de difracoes para secoes de afastamento nao nulo, mesmo para meios

homgeneos (Landa, 1987).

Considere-se um levantamento sısmico hipotetico com varios pares de fontes e recep-

tores com afastamento ΔxS R ao longo de uma superfıcie horizontal. Tiros e receptores sao

posicionados em xs e xr respectivamente. O modelo e composto por um unico difrator

localizado na coordenada (x0, h) e imerso em um meio homogeneo com velocidade cons-

tante v. O tempo de transito t para um raio que parte da fonte em (S) passa pelo difrator

(A) e chega ao receptor (R) (Figura 5.6) e

t =1

v

(√Δx2

s + h2 +

√Δx2

r + h2

), (D.1)

onde Δxs = x0 − xs e Δxr = xr − x0. xs, xr e x0 sao as coordenadas x da fonte, receptor e

148

difrator respectivamente.

Como ΔxS R = Δxs + Δxr, a equacao D.1 pode ser escrita por

t =1

v

( √(x0 − xs)2 + h2 +

√(xs + ΔxS R − x0)2 + h2

). (D.2)

A derivada da equacao D.2 relativa a posicao da fonte, xs, e

∂t∂xs=

1

vxs + ΔxS R − x0√

(xs + ΔxS R − x0)2 + h2− 1

vx0 − xs√

(x0 − xs})2 + h2. (D.3)

O tempo de transito mınimo ocorre quando a equacao D.3 iguala-se a zero. Entao, a

relacao entre x0, ΔxS R, xs, xr, Δxs e Δxr e

Δxs√Δx2

s + h2=

Δxr√Δx2

r + h2. (D.4)

Para satisfazer a igualdade na equacao D.4, Δxs deve ser igual a Δxr, o que significa

Δxs = Δxr =ΔxS R

2. O tempo duplo mınimo de uma difracao em um meio homogeneo

ocorre quando x0 e o ponto medio entre a a posicao da fonte e a posicao do receptor. Por

conseguinte, o tempo mınimo e

τ =

√4h2 + Δx2

S R

v(D.5)

que ocorre em

xm = x0 +ΔxS R

2. (D.6)

Se erroneamente o afastamento nao nulo ΔxS R for considerado como afastamento

zero, a velocidade aparente (vap) sera calculada por uma suposta distancia dupla entre

fonte (ou receptor) e difrator dividida por um tempo observado t.

vap = 2

√h2

w + (x − xm)2

t, (D.7)

onde hw e a suposta profundidade do difrator. Com alguma manipulacao algebrica,

t2 = 4h2

w

v2ap+ 4

(x − xm)2

v2ap

. (D.8)

O primeiro termo a direita da equacao D.8 e o tempo de transito mınimo observado τ,

descrito na equacao D.5. Entao, a velocidade aparente e:

v2ap = 4

(x − xm)2

t2 − τ2. (D.9)

149

Substituindo as equacoes D.6 e D.5 na equacao D.9,

v2ap = 4

(x0 +ΔxS R

2− x)2

t2 − 4h2+Δx2S R

v2

. (D.10)

O tempo de transito observado para este modelo e dado pela equacao D.1. Substi-

tuindo a equacao D.1 na equacao D.10:

v2ap = 4v2

(x0 +ΔxS R

2− x)2( √

Δx2s + h2 +

√Δx2

r + h2)2 − 4h2 + Δx2

S R

. (D.11)

Empregando as definicoes na equacao D.1 e considerando Δxs, Δxr, x0, xr e xs, a

equacao D.11 se torna

v2ap = 4v2

(x0 +ΔxS R

2− x)2( √

(x0 − xr + ΔxS R)2 + h2 +√

(x0 − xr)2 + h2)2 − 4h2 + Δx2

S R

. (D.12)

A partir da equacao D.12 ilustra-se que a a velocidade aparente vap nao e constante,

mas varia de acordo com a posicao do receptor e do afastamento ΔxS R. Se a diferenca

x0 − xr e muito maior que o afastamento ΔxS R (se os receptores estao muito afastados

dos difratores), vap tende a v. Igualmente, se ΔxS R tende a zero, vap tende a velocidade

real v. Todas as deducoes apresentadas para vap sao validas para meios isotropicos com

variacoes verticais da velocidade. Contudo, v deve ser substituıdo pela velocidade RMS.

150

Documents

INVERSAO TOMOGR˜ AFICA SEQUENCIAL PARA O CAMPO DE