MODELAGEM SÍSMICA EM MEIOS COMPLEXOS

Eldues Oliveira Martins

TESE SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS

PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE

FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS

NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE EM CIÊNCIAS EM

ENGENHARIA CIVIL.

Aprovada por:

Prof. Luiz Landau, D.Sc.

Dr. Djalma Manoel Soares Filho, D.Sc.

Prof. José Luis Drummond Alves, D.Sc.

Prof. Liliana Alcazar Diogo, D.Sc.

Prof. Wilson Mouzer Figueiró, D.Sc.

RIO DE JANEIRO, RJ - BRASIL

MAIO - 2003

ii

MARTINS, ELDUES OLIVEIRA

Modelagem Sísmica em Meios

Complexos [Rio de Janeiro] 2003.

XVI, 132p. 29,7 cm (COPPE/UFRJ,

M.Sc., Engenharia Civil, 2003)

Tese – Universidade Federal do Rio

de Janeiro, COPPE

1. Geofísica;

2. Modelagem Sísmica;

3. Anisotropia;

I. COPPE/UFRJ II. Título (série).

iii

“À minha querida e paciente esposa Luciana, que me deu meus muito amados

filhos Pedro, Mariana e João, pelo seu amor, suporte, atenção e compreensão

que foram fundamentais para a conclusão desta tese; e enfim,

a todos que amo e admiro.”

iv

AGRADECIMENTOS

Agradeço ao coordenador e orientador Prof. Luiz Landau e a todo corpo técnico-

administrativo do Laboratório de Métodos Computacionais em Engenharia – LAMCE pela

oportunidade da realização deste projeto.

Agradeço ao meu orientador Dr. Djalma M. Soares Filho pelo seu incansável esforço

no desenvolvimento da pesquisa científica, pela sua brilhante dedicação na orientação desta

tese e valorosa contribuição na minha formação profissional.

Agradeço ao geolólogo e amigo Ricardo Bedregal pelas proveitosas discussões sobre a

geologia do petróleo, pelos incentivos e apoio incondicional durante a realização desta tese.

Agradeço ao Prof. José Luis Drummond Alves e ao colega Dênis Araujo Figueiras de

Souza pelas contribuições na otimização visualização e desenvolvimento dos algoritimos em

Fortran.

Agradeço aos sempre amigos Josias J. da Silva, Alexandre Lopes, Ricardo Alencar,

Luis Fernando, Magda Almada, Mônica Caruso e a todos os outros que com certeza, sem eles,

a caminhada seria muito mais difícil.

Agradeço ao geofísico da PETROBRAS/CENPES/PDEXP/GEOF João Cláudio

Conceição da Silva pela colaboração para a realização deste trabalho.

v

Resumo da Tese apresentada a COPPE/UFRJ como parte dos requisitos necessários para a

obtenção do grau de Mestre em Ciências (M. Sc.)

MODELAGEM SÍSMICA EM MEIOS COMPLEXOS

Eldues Oliveira Martins

Maio/2003

Orientadores: Djalma Manoel Soares Filho

Luiz Landau

Programa: Engenharia Civil.

Foram desenvolvidas e implementadas duas metodologias para simulações numéricas de levantamentos sísmicos em meios complexos através da solução do sistema de equações que regem o comportamento ondulatório em meios transversalmente isotrópicos bidimensionais, usando o método das diferenças finitas. Para meios terrestres, foi generalizado o algoritmo proposto por Zahradnik e Priolo [1994], no qual os componentes do campo de deslocamento das partículas são propagados e os parâmetros elásticos são introduzidos por integrações ao longo das linhas definidas pela malha, suposta regular. Utilizam-se aproximações de segunda ordem para as derivadas parciais na geração do operador de propagação. Para abordagens marinhas, modificou-se o algoritmo apresentado em Levander [1986], no qual os campos de velocidade e tensão, assim como os parâmetros elásticos, são definidos em malhas intercaladas regulares. Para este método, utilizam-se aproximações de segunda ordem para as derivadas temporais e quarta ordem para as espaciais, na geração dos operadores que realizam a propagação dos campos de velocidades e tensão. Os métodos não estão restritos a qualquer tipo de geometria de aquisição (fontes e receptores podem estar posicionados em qualquer lugar no modelo), assim como não há qualquer limitação no que diz respeito ao modelo geológico, desde que as condições que garantem a estabilidade e ausência de dispersões numéricas sejam satisfeitas. No que diz respeito às aplicações, inicialmente, foi estudado o comportamento ondulatório em função do grau de anisotropia TIV. Especificamente, verificou-se em instantâneos, obtidos em modelagens para meios homogêneos, a influência dos parâmetros de Thomsen [1986] na propagação de ondas compressionais e cisalhantes. Em seguida, foram obtidos sismogramas e instantâneos em simulações sísmicas em dois modelos com altos contrastes de velocidades, típicos de bacias brasileiras. Os resultados foram analisados e comparados com os obtidos em simulações elásticas isotrópicas.

vi

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

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

SEISMIC MODELING ON COMPLEX MEDIA

Eldues Oliveira Martins

Maio/2003

Advisors: Djalma Manoel Soares Filho

Luiz Landau

Departament: Civil Engineering.

Two finite differences type methods are introduced for seismic modeling on transversely isotropic media. The first one is designed for seismic acquisition simulations on on-shore geologic models and consists of a modification of the Zahradník and Priolo (1994) proposal to accommodate seismic modeling on anisotropic media. In this method, the components of the displacement vector are discretized in a regular grid and the elastic tensor components are introduced through vertical and horizontal integrations along the grid lines. Second order approximations are used for spatial and temporal derivatives and Vacuum Formalism, as it is posed in Zahradník, Moczo and Heron (1993), is used as boundary condition along the upper interface. This technique is not restricted to flat observation surfaces, but it cannot be used on off-shore models, since it gives rise to instabilities where shear modulus is small. The second method is based on Levander (1986) proposal technique and is indicated to seismic simulations on off-shore geologic models. As it was done with the Zahradník and Priolo (1986) method, we also adapt this second one to cope with transversely isotropic media. In this case, the components of the velocity field and the elastic parameters are defined in a staggered grid and the conditions for preventing spurious numerical dispersions do not depend on shear modulus value. This fact makes off-shore seismic simulations feasible using this approach. In relation to partial derivatives approximations, the temporal derivatives are approximated in second order whereas the spatial ones are fourth-order approximated. Both methods are not restricted to any particular acquisition pattern (sources and receivers may be located anywhere) and may be applied to any geologic case (as far as the conditions for preventing spurious dispersions are satisfied). To illustrate, we show snapshots and seismograms obtained on several simulations using the proposed methods and indicate some of the characteristics of the wave field propagation on realist transversely isotropic media.

Índice

Página de Assinaturas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i

Ficha Catalográfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii

Dedicatória . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii

Agradecimentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv

Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi

Índice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

Índice de Figuras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x

Índice de Tabelas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi

1 Introdução 1

1.1 Elasticidade do Meio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2 Anisotropia do Meio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.3 Heterogeneidade do Meio . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Modelagem Sísmica para Meios Elásticos TIV: Caso Terrestre 10

2.1 Desenvolvimento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.1.1 Discretização das Equações da Elastodinâmica para meios TIV . . . 12

vii

2.1.2 Condições de Contorno para as Bordas Laterais e Inferior . . . . . . 34

2.2 Termo Fonte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.3 Definição da malha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

2.3.1 Condições para Não Ocorrência de Dispersão Numérica . . . . . . . 37

2.3.2 Condições de Estabilidade Numérica . . . . . . . . . . . . . . . . . 38

2.4 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

2.5 Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

2.5.1 Primeira Aplicação: Propagação de Onda em Meios Homogêneos

Anisotrópicos nos Modelos Propostos em Thomsen (1986) . . . . . 42

2.5.2 Segunda Aplicação: Modelo Formado por Duas Camadas Separadas

por Interface Curva . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

2.5.3 Terceira Aplicação: Modelo com Altos Contrastes de Impedância . 58

3 Modelagem Sísmica para Meios Elásticos TIV: Caso Marinho 66

3.1 Desenvolvimento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.1.1 Discretização pelo Método de Diferenças Finitas . . . . . . . . . . . 69

3.2 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

3.3 Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

3.3.1 Primeira Aplicação: Modelagem em Meios Acústicos com operado-

res Elásticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

3.3.2 Segunda Aplicação: Modelagem Sísmica em Modelo Típico de Mar-

gens Passivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

viii

Conclusão 92

Referências 96

Apêndices 101

A Tensor de Elasticidade (Cijkl) Expresso no Sistema de Coordenadas Prin-

cipal 102

A.1 Componentes da tensão e deformação . . . . . . . . . . . . . . . . . . . . . 103

B Discretização por Diferenças Finitas 113

B.1 Malhas Intercaladas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

ix

Lista de Figuras

1.1 Direção de polarização e propagação definida para um caso anisotrópico de

meio com finas camadas paralelas. . . . . . . . . . . . . . . . . . . . . . . . 6

1.2 Modelos litológicos equivalentes para diferentes tipos de simetrias. . . . . . 7

2.1 Linhas da malha onde são calculadas as médias geométricas dos parâmetros. 16

2.2 Modelo da malha do método PS2 de diferenças finitas para aproximação

das derivadas em relação a uma única variável para pontos da malha na

superfície. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3 Modelo da malha do método PS2 de diferenças finitas para aproximação

das derivadas em relação a uma única variável para pontos da malha no

interior do modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4 Modelo da malha do método PS2 de diferenças finitas para aproximação das

derivadas em relação às variáveis x e z para pontos da malha na superfície

SE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.5 Modelo da malha do método PS2 de diferenças finitas para aproximação das

derivadas em relação às variáveis x e z para pontos da malha na superfície

SW. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

x

2.6 Modelo da malha do método PS2 de diferenças finitas para aproximação das

derivadas em relação às variáveis x e z para pontos da malha na superfície

NE e NW. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.7 Modelo da malha do método PS2 de diferenças finitas para aproximação

das derivadas em relação às variáveis x e z para pontos da malha no interior

WN, WS, EN e ES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.8 Gráfico da função fonte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

2.9 Diagrama representando a 1ametodologia para modelagem em meios elás-

ticos TIV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2.10 Instantâneos do componente vertical no tempo de 3 ms. a) alumínio-lucita.

b) folhelho anisotrópico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

2.11 Instantâneos do componente vertical no tempo de 3 ms. a) apatita. b)

biotita. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

2.12 Instantâneos do componente vertical no tempo de 3 ms. a) calcita. b)

arenito calcarenoso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

2.13 Instantâneos do componente vertical no tempo de 3 ms. a) folhelho anisotrópico-

ss. b) folhelho argiloso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

2.14 Instantâneos do componente vertical no tempo de 3 ms. a) calcáreo síltico.

b) arenito-folhelho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

2.15 Instantâneos do componente vertical no tempo de 3 ms. a) arenito. b)

cristal de quartzo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

2.16 Instantâneos do componente vertical no tempo de 3 ms. a) folhelho com

óleo. b) folhelho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

xi

2.17 Instantâneos do componente vertical no tempo de 3 ms. a) folhelho anisotrópico-

ls. b) calcáreo-folhelho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

2.18 Instantâneos do componente vertical no tempo de 3 ms. a) siltito laminado.

b) arenito imaturo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

2.19 Instantâneos do componente vertical no tempo de 3 ms. a) gelo b) material

de gipso intemperizado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

2.20 Instantâneos do componente vertical no tempo de 3 ms. a) gas-arenito-

água-arenito. b) cristal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

2.21 Instantâneos do componente vertical no tempo de 3 ms. a) cinza vulcânica. 49

2.22 Modelo de velocidade abordado na segunda aplicação. . . . . . . . . . . . . 51

2.23 Instantâneos do componente vertical do campo de onda, considerando fonte

explosiva na posição central domodelo homogêneo, VP=3000m/s e VS=1732

m/s. a) Isotrópico. b) anisotrópico ε = δ = 0, 1. c) anisotrópico ε = δ =

0, 2. d) anisotrópico ε = δ = 0, 3. . . . . . . . . . . . . . . . . . . . . . . . 52

2.24 Instantâneos do componente horizontal do campo de onda, considerando

fonte explosiva na posição central do modelo homogêneo, VP=3000 m/s e

VS=1732 m/s. a) Isotrópico. b) anisotrópico ε = δ = 0, 1. c) anisotrópico

ε = δ = 0, 2. d) anisotrópico ε = δ = 0, 3. . . . . . . . . . . . . . . . . . . . 53

2.25 Sismogramas relativos ao componente vertical do campo de ondas. . . . . . 54

2.26 Sismogramas relativos ao componente horizontal do campo de ondas. . . . 55

2.27 Sismogramas relativos ao componente vertical do campo de ondas, sem a

presença das ondas superficiais. . . . . . . . . . . . . . . . . . . . . . . . . 56

xii

2.28 Sismogramas relativos ao componente horizontal do campo de ondas, sem

a presença das ondas superficiais. . . . . . . . . . . . . . . . . . . . . . . . 57

2.29 Modelo de velocidades da bacia sedimentar terrestre estudada. . . . . . . . 59

2.30 Distribuição dos parâmetros de Thomsen no modelo para o caso terrestre. . 59

2.31 Instantâneos relativos ao componente vertical para os modelos isotrópico

(à esquerda) e anisotrópico nos tempos de propagação 0,5; 0,75; 1,00 e 1,25 s. 61

2.32 Instantâneos relativos ao componente horizontal para os modelos isotrópico

(à esquerda) e anisotrópico nos tempos de propagação 0,5; 0,75; 1,00 e 1,25 s. 62

2.33 Sismogramas relativos ao componente vertical para os modelos isotrópico (à

esquerda) e anisotrópico (ε = δ = 0, 2, na camada que apresenta Vp = 2500

m/s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

2.34 Sismogramas relativos ao componente horizontal para os modelos isotrópico

(à esquerda) e anisotrópico (² = δ = 0, 2, na camada que apresenta Vp =

2500 m/s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

2.35 Ondas superficiais, especificamente as ondas diretas P e S, e as ondas

Rayleigh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

3.1 Instantâneos relativos ao componente vertical, ilustrando os problemas de

dispersão numérica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.2 Diagrama de velocidade. Malha intercalada de diferenças finitas e diagrama

espacial para as atualizações das velocidades. . . . . . . . . . . . . . . . . . 70

3.3 Diagrama de tensões. Malha intercalada de diferenças finitas e diagrama

espacial para as atualizações de tensões. . . . . . . . . . . . . . . . . . . . 72

xiii

3.4 Diagrama representando a 2ametodologia para modelagem em meios elás-

ticos TIV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

3.5 Instantâneos relativos aos componentes horizontal e vertical do campo de

velocidade (derivada do campo de deslocamento), U e V, aos componentes

do tensor de esforços τxx, τxz e τ zz, respectivamente, considerando ummeio

acústico (µ = 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

3.6 Modelo de velocidades da bacia sedimentar marítima estudada com os va-

lores de velocidade da onda P em cada camada. . . . . . . . . . . . . . . . 82

3.7 Modelo de densidade no modelo para o caso marinho. . . . . . . . . . . . . 82

3.8 Distribuição dos parâmetros de Thomsen no modelo para o caso marinho. . 83

3.9 Instantâneos relativos ao componente vertical do campo de velocidades

obtidos na simulação sísmica no modelo para o caso marinho. . . . . . . . 84

3.10 Instantâneos relativos ao componente horizontal do campo de velocidades

obtidos na simulação sísmica no modelo para o caso marinho. . . . . . . . 85

3.11 Instantâneos relativos ao componente τxx do campo de tensões obtidos na

simulação sísmica no modelo para caso marinho. . . . . . . . . . . . . . . . 86

3.12 Instantâneos relativos ao componente τxz do campo de tensões obtidos na

simulação sísmica no modelo para o caso marinho. . . . . . . . . . . . . . . 87

3.13 Instantâneos relativos ao componente τ zz do campo de tensões obtidos na

simulação sísmica no modelo para o caso marinho. . . . . . . . . . . . . . . 88

3.14 Sismogramas relativos ao componente vertical do campo de velocidades

obtidos na simulação sísmica no modelo para o caso marinho. . . . . . . . 89

xiv

3.15 Sismogramas relativos ao componente horizontal do campo de velocidades

obtidos na simulação sísmica no modelo para o caso marinho. . . . . . . . 89

3.16 Sismogramas relativos ao componente τxx do campo de tensões obtidos na

simulação sísmica no modelo para o caso marinho. . . . . . . . . . . . . . . 90

3.17 Sismogramas relativos ao componente τxz do campo de tensões obtidos na

simulação sísmica no modelo para o caso marinho. . . . . . . . . . . . . . . 90

3.18 Sismogramas relativos ao componente τ zz do campo de tensões obtidos na

simulação sísmica no modelo para o modelo marinho. . . . . . . . . . . . . 91

xv

Lista de Tabelas

2.1 Parâmetros da malha usada nas modelagens na primeira aplicação para o

caso terrestre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

2.2 Parâmetros de Thomsen usados na primeira aplicação para o caso terrestre. 43

2.3 Parâmetros da malha usada nas modelagens na segunda aplicação para o

caso terrestre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

2.4 Parâmetros da malha usada nas modelagens na terceira aplicação para o

caso terrestre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.1 Parâmetros da malha usada nas modelagens na primeira aplicação para o

caso marinho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

3.2 Parâmetros da malha usada nas modelagens na segunda aplicação para o

caso marinho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

xvi

Capítulo 1

Introdução

Amodelagem sísmica tem sido uma das ferramentas mais importantes à dis-

posição dos profissionais envolvidos nos processos de exploração e explotação de hidrocar-

bonetos através de métodos sísmicos. De fato, de posse de um conhecimento satisfatório

da resposta sísmica de cada componente do sistema petrolífero, intérpretes (geólogos,

geofísicos e engenheiros) podem trabalhar com mais segurança; através de simulações

numéricas de aquisições sísmicas, os geofísicos responsáveis pela definição dos parâmetros

de levantamento podem otimizar os custos envolvidos com o mesmo; e os profissionais

engajados no desenvolvimento de novas tecnologias de processamento e interpretação po-

dem inferir a efetividade de cada técnica com aplicações em dados simulados, explorando

o conhecimento exato do modelo geológico que os geraram. Com relação às ferramentas

atualmente utilizadas no processamento de dados de campo, a modelagem encontra-se no

cerne dos algoritmos. Por exemplo, na migração pré-empilhamento em profundidade, a

solução da equação da onda é necessária nas etapas de depropagação do campo de onda

registrado e na geração das matrizes de tempo de trânsito.

1

O contínuo e acelerado desenvolvimento dos métodos de aquisição, proces-

samento e interpretação sísmica têm permitido um aumento da recuperação e produção

das reservas petrolíferas. Na última década a exploração, explotação e caracterização de

reservas em estruturas e estratigrafias complexas e em condições especiais como fratu-

ramentos, por exemplo, tem sido intensificado, requisitando assim a necessidade de se

desenvolver novas tecnologias capazes de se reconstruir modelos adequados de subsuper-

fície e determinando propriedades específicas, como a orientação e densidade de fraturas,

porosidade, permeabilidade, presença e saturação de fluidos, etc [1].

Uma propriedade fundamental de subsuperfície é a anisotropia, que atual-

mente na indústria raramente é levada em consideração em função do dispêndio computa-

cional, que aumenta os custos de processamento, na inversão dos parâmetros elásticos e

migração. Se entende que a variação das propriedades físicas de acordo com a orientação

na qual se realiza a medida para a modelagem geológica é de grande importância, como

por exemplo, no fenômeno de birrefrigência das ondas de cisalhamento, no qual a onda

sísmica se divide em duas frentes de onda que se propagam com velocidades diferentes (o

movimento das partículas na onda mais rápida é paralela à direção de fraturamento do

meio).

Os dados de onda convertida implicam um novo paradigma de processa-

mento, pois as reflexões sobre eventos em subsuperfície ocorrem de forma totalmente

assimétrica com respeito às reflexões de um só modo de propagação. Estas reflexões

dependem em grande parte do quociente de velocidade das ondas P e S, assim como

da profundidade do objetivo, no qual aumenta a complexidade do processamento para

eventos em subsuperfície.

2

O começo da sismologia usando diferenças finitas para resolver problemas

com propagação de onda tem aproximadamente 50 anos. Geologias complexas, prove-

nientes de dobramentos e falhamentos decorrentes de variações de velocidades laterais

fortes e relevos de topografias severas são os principais responsáveis por dificultar o ima-

geamento sísmico em subsuperfície fidedignamente (Soares Filho et al) [13].

Problemas modelando geofísica de exploração quase sempre dividem-se em

estruturas complexas ou estratigrafias complexas caracterizados por interfaces irregulares

separando unidades litológicas maiores. Então para ser efetivo, um algoritmo de mode-

lagem numérica deve automaticamente aproximar as condições de contorno através destas

interfaces. Cada um algoritmo heterogêneo de diferenças finitas soluciona um sistema de

equações de segunda ordem em uma malha singular para malha normal, aplicado usual-

mente em modelagem terrestre, que foi introduzido por Alford (1974) [14] e estendido por

Kelly et al (1976) [15] ou quarta ordem em malha intercalada, aplicado usualmente em

modelagem marítma, que foi introduzido por Levander (1990) [4].

A técnica de diferenças finitas particularmente é atrativa para geometrias em

subsuperfícies estruturalmente complexas porque as condições de contorno dos contatos

litológicos são estimadas implicitamente pelo operador de diferenças finitas. O método de

diferenças finitas aproxima a equação da onda pela diferencial de segunda ordem, estas são

resolvidas iterativamente em cima de uma malha de discretização espacial em esquema

de pequenas extrapolações.

A utilização do método de diferenças finitas para criar sismogramas sin-

téticos não é novo. Cada vez mais surgem trabalhos de sismologia usando esta técnica

e discutindo em detalhes aspectos particulares. O programa para modelagem terrestre

3

em malha normal, desenvolvido neste trabalho, utiliza aproximação de segunda ordem

para ambos os intervalos, espacial e temporal, enquanto que o desenvolvido para a mode-

lagem marítima usa aproximação de quarta ordem no espaço e segunda ordem no tempo

em malha intercalada. Em ambos os programas foram utilizadas condições de absorção

de energia. O método de diferenças finitas nos permite obter fotografias (“snapshots”)

do campo de ondas em qualquer tempo da propagação, (“common shot point”) com sis-

mogramas em formato de dados padrões e ele tem a capacidade de recorrer a múltiplos

sismogramas sintéticos em um ponto comum em profundidade (“common depth point”)

simulando a aquisição dos dados.

Foi utilizado para aplicação de filtros e visualização de snapshots e sismo-

gramas o pacote de processamento sísmico “Seismic Unix”, desenvolvido pelo “Center for

Wave Phenomena - Colorado School of Mines”.

Este pacote é livremente distribuído, juntamente com todos os códigos

fontes.

Simogramas sintéticos de modelos complexos são úteis na exploração, servin-

do para ilustrar como a técnica pode ajudar o intérprete. Como exemplo podemos citar

a definição dos vários parâmetros na implementação da modelagem por diferenças finitas,

que envolvem dispersão da malha, reflexões artificiais vindo da borda do modelo, definição

dos incrementos das amostras espaciais e temporais do modelo dentre outros,Kelly et al

(1976) [15].

4

1.1 Elasticidade do Meio

As perturbações que se propagam em meio elástico denominam-se ondas

elásticas. Este meio é constituído de qualquer material que tenda a preservar seu compri-

mento, forma e volume contra as forças externas.

Os materiais que possuem forças restauradoras que tendem a retornar o

material à sua condição original após a remoção das forças externas são meios elásticos.

A força restauradora é característica do material e tem origem nas forças de ligação entre

seus átomos ou moléculas individuais.

Refere-se como meio não-dispersivo aquele em que a forma da onda não

se altera à medida que a onda se propaga e sua velocidade é constante desde que sejam

fixadas as características de elasticidade e a densidade do meio, pois a velocidade depende

em geral, destes parâmetros. A onda sonora no ar é um exemplo de onda que não sofre

dispersão.

1.2 Anisotropia do Meio

A anisotropia é um termo geral que se refere às mudanças nas propriedades

físicas resultantes das mudanças na direção ou polarização (Figura 1.1). Anisotropia na

velocidade sísmica é uma característica comum em muitos materiais da Terra, mas, em

geral, não é levado em consideração no processamento de dados sísmicos.

Virtualmente, todos os atuais algoritmos de processamento em uso assumem

que a Terra é isotrópica, porém, na verdade, a Terra é anisotrópica em quase todas as

suas características.

5

Direção da linha

Direção da linha

Figura 1.1: Direção de polarização e propagação definida para um caso anisotrópico demeio com finas camadas paralelas. As polarizações são definidas em termos da linhade direção indicada e são consideradas ambas as direções de propagação horizontal evertical. Na Figura 1.1 (a) é considerado um acamamento horizontal, enquanto na Figura1.1 (b) o acamamento é vertical. As Figuras mostram o desmembramento das diferentesvelocidades com diferentes direções de propagação e polarização. Figura adaptada deTatham and McCormack (1991) [16].

6

a)

c)

b)

d)

Isotropia TransversalVertical

Isotropia TransversalHorizontal

Ortorrombica Monoclinica

Figura 1.2: Modelos litológicos equivalentes para diferentes tipos de simetrias.

Vale a pena ressaltar que as considerações de isotropia são feitas já muito

cedo no desenvolvimento das definições da equação que rege a propagação do campo

de onda. Se uma anisotropia média for considerada, o problema se torna muito mais

complexo. Para se manter o problema controlado, algumas considerações sobre a simetria

e geometria da anisotropia são comumente feitas.

Neste trabalho foi usada isotropia transversalmente vertical (TIV), como

mostra a Figura 1.2 (a), nestes experimentos os materiais anisotrópicos assumiram isotropia

transversa (sem variação azimutal na velocidade) onde as ondas S são descritas pela po-

larização como ondas SV e SH. As ondas SV e SH podem ter diferentes velocidades.

7

1.3 Heterogeneidade do Meio

Na formulação para meios heterogêneos as propriedades elásticas podem ser

especificadas para cada ponto da malha de diferenças finitas, e as condições de contorno

são satisfeitas implicitamente, Kelly et al (1976) [15].

Omeio elástico pode ser considerado como um conjunto de regiões litológicas

homogêneas, cada uma caracterizada por valores constantes de densidade e parâmetros

elásticos.

A propagação do campo de onda em cada região pode ser apresentada por

uma representação de diferenças finitas apropriada para as equações elásticas correspon-

dentes. Por este método, todas as condições de contorno através das interfaces de se-

paração de diferentes regiões podem ser satisfeitas explicitamente. As aproximações são

incorporadas às condições de contorno implicitamente, construindo representação de di-

ferenças finitas para equações elásticas gerais que são válidas para regiões heterogêneas.

Estas aproximações tornam possível associar diferentes densidades e valores

de parâmetros elásticos a todo ponto da malha, trazendo flexibilidade na variação de

geometrias de superfícies complexas.

A dissertação aqui apresentada tem o objetivo de introduzir dois métodos

para modelagens sísmicas em meios transversalmente isotrópicos, compatíveis com as es-

pecificidades encontradas nos modelos geológicos abordados nos processos de exploração e

explotação de hidrocarbonetos. Para alcançar este intento foi preciso resolver diretamente

a equação diferencial que rege o comportamento ondulatório por estratégias que não en-

volvessem soluções particulares, tais como, as utilizadas nos métodos da refletividade e

8

traçamento de Raios. Especificamente, apresentamos soluções baseadas no método de

diferenças finitas, no qual as derivadas parciais espaciais e temporais foram discretizadas

com aproximações de segunda e quarta ordens. Para o caso de abordagens terrestres,

generalizamos o algoritmo PS2 de Zahradník e Priolo (1994) [2] e para o caso marinho,

implementamos o método apresentado em Levander (1990) [4], com modificações para

que o mesmo pudesse tratar casos anisotrópicos. Em ambos os métodos, desde que sejam

satisfeitas as condições necessárias para garantir a não existência de dispersões numéricas

e estabilidade, não existe qualquer limitação com respeito à geometria das camadas que

compõem o modelo. Assim como não há qualquer restrição à variação dos parâmetros

elásticos no interior da mesma.

Esta dissertação apresenta três capítulos. O primeiro capítulo diz respeito

a esta introdução; o segundo apresenta o método que generaliza o de Zahradník e Priolo

(1994) [2], indicado para modelagens terrestres; e o terceiro apresenta a modificação do

método de Levander (1990) [4], visando modelagens marinhas. No final do segundo e

terceiro capítulos, são apresentados exemplos de aplicação dos métodos propostos. Após

o terceiro capítulo, apresentamos as conclusões, comentários e propostas para trabalhos

futuros nesta área. Concluindo: apresentamos dois apêndices nos quais exibimos o de-

senvolvimento matemático empregado para expressar o tensor de elasticidade em meios

transversalmente isotrópicos em termos do sistema de coordenadas que melhor explora

esta simetria (simetria principal), e explicitamos a discretização em malhas intercaladas.

9

Capítulo 2

Modelagem Sísmica para Meios

Elásticos TIV: Caso Terrestre

Neste capítulo apresentamos uma generalização do algoritmo PS2 intro-

duzido por Zahradník e Priolo (1994) [2] visando modelagens sísmicas 2-D por diferenças

finitas em meios transversalmente isotrópicos, com eixo de simetria vertical. O algoritmo

é indicado para simulações de aquisições sísmicas em modelos terrestres e pode ser aplica-

do nos casos de superfícies de observação irregulares. Quanto à complexidade estrutural e

estratigráfica dos modelos que podem ser abordados por essa técnica, não existe qualquer

restrição com relação ao número de camadas, que podem ser separadas por interfaces

arbitrariamente irregulares, assim como não há qualquer limitação quanto à variação dos

parâmetros no interior das camadas, desde que sejam satisfeitas as condições que garan-

tem a inexistência de dispersões numéricas e instabilidades. O método proposto é capaz

de abordar modelos geológicos realistas, sendo, portanto, uma valiosa ferramenta a ser

usada por geólogos e geofísicos envolvidos nos processos exploratórios e explotatórios de

10

hidrocarbonetos.

Na formulação do método apresentado neste capítulo, temos:

(1) as equações que regem o campo de deslocamento das partículas são

discretizadas, numa malha regular, com aproximações de segunda ordem nas derivadas

parciais;

(2) os parâmetros são introduzidos por integrações ao longo das linhas da

malha;

(3) o termo fonte é dado pela segunda derivada da função gaussiana, como

apresentado em Cunha (1997) [6];

(4) as condições de contorno não reflexivas são as propostas por Emerman

e Stephen (1990) [10] combinadas com as bordas de absorção sugeridas por Cerjan et al

(1985) [7] e

(5) aplica-se a condição de contorno conhecida por Formalismo do Vácuo,

apresentada em Zahradník, Moczo e Hron (1993) [3], na superfície de observação.

O programa Fortran desenvolvido tem como dados de entrada a densidade ρ,

os parâmetros anisotrópicos ε e δ de Thomsen (1986) [5], e as velocidades das ondas com-

pressionais VP e cisalhantes VS na direção vertical, ao longo do modelo. Os componentes

do tensor elástico Cijkl são computados a partir destes parâmetros e são integrados ao

longo das linhas da malha. Os dados de saída são instantâneos dos componentes do cam-

po de deslocamento e sismogramas, especificamente os componentes vertical e horizontal

do campo de deslocamentos em todos os tempos registrados na superfície de observação.

A generalização para modelagens 3-D pode ser realizada explorando-se as facilidades de

paralelização do código.

11

Com relação às aplicações:

• Com o intuito de estudar o efeito da anisotropia na propagação do campo de ondas,

foram realizadas modelagens numéricas para simulação da propagação de onda em 23

meios dos presentes no artigo de Thomsen (1986). Como resultado destas simulações

serão apresentados instantâneos da componente vertical da velocidade para todos

os modelos abordados.

• Com o objetivo de estimar os efeitos decorrentes da anisotropia nas reflexões sísmicas

foram realizadas simulações em ummodelo composto por duas camadas homogêneas

sendo a primeira anisotrópica, separadas por uma interface curva (“senoidal”).

• Para mostrar a capacidade do método proposto emmodelagens emmeios complexos,

foram obtidos instantâneos e sismogramas para os componentes vertical e horizontal

do campo de ondas num modelo com alta complexidade geológica.

2.1 Desenvolvimento

2.1.1 Discretização das Equações da Elastodinâmica para meios

TIV

Considere o sistema de equações que rege o comportamento ondulatório em

meios elásticos heterogêneos transversalmente isotrópicos (com eixo de simetria vertical),

12

bidimensional, isto é:

ρ∂2u

∂t2=

∂

∂x

µC11

∂u

∂x

¶+

∂

∂x

µC13

∂v

∂z

¶+

+∂

∂z

µC44

∂u

∂z

¶+

∂

∂z

µC44

∂v

∂x

¶, (2.1)

ρ∂2v

∂t2=

∂

∂x

µC44

∂v

∂x

¶+

∂

∂z

µC11

∂u

∂x

¶+

+∂

∂z

µC33

∂v

∂z

¶+

∂

∂x

µC44

∂u

∂z

¶. (2.2)

Onde u e v denotam, respectivamente, os componentes horizontal e vertical

do campo de deslocamento e C11, C13, C44 e C33 são os componentes do tensor de elasti-

cidade, definidos nas Equações A.19, A.20, A.18 e A.17 apresentadas no apêndice, e ρ é a

densidade. Temos que: n−1, n e n+1 representam os índices referentes, respectivamente,

ao passo de tempo anterior, atual e posterior de valores de u e v de um dado ponto da

malha, e h representa a variação espacial.

Introduzindo os índices citados anteriormente i, j e n, teremos:

u(x, y, t) = uni,j = ui,j , (2.3)

e

v(x, y, t) = vni,j = vi,j , (2.4)

13

com

n = 0, 1, 2, 3, ... . (2.5)

Denotando o lado direito da Equação 2.1 por $1(x, z), ou seja escrevendo

ρ∂2u

∂t2= $1(x, z), (2.6)

e aproximando em 2aordem a derivada parcial temporal temos:

un+1i,j − 2uni,j + un−1i,jDT 2

∼= 1ρi,j$n1(i, j) , (2.7)

ou ainda,

un+1i,j =DT 2

ρi,j$n1(i, j) + 2u

ni,j − un−1i,j , (2.8)

onde,

i = 1, ..., Nx;

j = 1, ..., Nz;

n = 1, ..., Nt; e

DT é o intervalo de tempo usado na discretização a ser definido na Equação

2.76.

14

Analogamente, podemos escrever:

vn+1i,j =DT 2

ρi,j$2(i, j) + 2v

ni,j − vn−1i,j , (2.9)

onde $2(i, j) representa as derivadas espaciais da Equação 2.2.

Na discretização dos termos $1(i, j) e $2(i, j), consideramos dois tipos de

derivadas parciais espaciais, denominadas simples e mista. As derivadas simples envolvem

uma única variável espacial, isto é, são do tipo,

∂

∂x

µa∂f

∂x

¶,

∂

∂z

µa∂f

∂z

¶.

As derivadas mistas são do tipo,

∂

∂z

µa∂f

∂x

¶,

∂

∂x

µa∂f

∂z

¶.

Onde a e f são funções de x e z, no caso específico um dos parâmetros elásticos e um dos

campos de deslocamento, respectivamente.

Assim sendo serão discretizados esses dois tipos de derivadas. Inicialmente

considere a do tipo simples,

∂

∂x

µa∂f

∂x

¶. (2.10)

15

E W

N

S

i,j+1

i,j-1

i-1,j i+1,j

Figura 2.1: Linhas da malha onde são calculadas as médias geométricas dos parâmetros.

A discretização de

∂

∂z

µa∂f

∂z

¶, (2.11)

é, naturalmente, similar.

A derivada simples pode ser reescrita como

∂

∂xg,

onde

g = a∂f

∂x. (2.12)

16

(i-1/2,j) (i+1/2,j)

(i,j)

h

h

(i-1/2,j) (i+1/2,j)

(i,j)

h

h

Figura 2.2: Modelo da malha do método PS2 de diferenças finitas para aproximação dasderivadas em relação a uma única variável para pontos da malha na superfície.

Assim podemos aproximá-la por,

∂

∂xg ∼= g(i+

12, j)− g(i− 1

2, j)

h, (2.13)

onde h é a distância horizontal e vertical entre dois pontos consecutivos na malha regular,

Figura 2.2, sendo:

g(i± 12, j) ∼= a∂f

∂x|(i± 1

2,j) . (2.14)

Os valores de g, Equação 2.12, nos pontos (i± 12, j) serão estimados a partir das seguintes

aproximações:

g(i+1

2, j) ∼=

R (i+1,j)(i,j)

∂f

∂xdxR (i+1,j)

(i,j)

1

adx

=f(i+ 1, j)− f(i, j)R (i+1,j)

(i,j)

1

adx

, (2.15)

17

(i-1,j) (i+1,j)(i,j)

(i,j-1)

(i,j+1)

N

S

W E

(i-1,j) (i+1,j)(i,j)

(i,j-1)

(i,j+1)

N

S

W E

Figura 2.3: Modelo da malha do método PS2 de diferenças finitas para aproximação dasderivadas em relação a uma única variável para pontos da malha no interior do modelo.

ou

g(i+1

2, j) ∼= aE

h(fi+1,j − fi,j) , (2.16)

onde aE é definido por

aE =hR (i+1,j)

(i,j)

1

adx=

hRE

1

adx, (2.17)

na qual o símbolo E na integral denota integração ao longo do segmento horizontal que

vai de (i, j) a (i+ 1, j), Figura 2.3. Analogamente, o valor de g em (i− 12, j) é estimado

18

por

g(i− 12, j) ∼= aW

h(fi,j − fi−1,j) , (2.18)

no qual

aW =hR (i,j)

(i−1,j)1

adx=

hRW

1

adx, (2.19)

e W indica que a integração de1

ano segmento horizontal que vai de (i − 1, j) a (i, j),

Figura 2.3.

Assim sendo, a aproximação usada no cálculo das derivadas parciais simples

em relação a variável x será:

∂

∂x

µa∂f

∂x

¶∼= 1h2{aE(fi+1,j − fi,j)− aW (fi,j − fi−1,j)}. (2.20)

De forma análoga, pode-se encontrar a fórmula que aproximará as derivadas simples em

relação a variável z, isto é,

∂

∂z

µa∂f

∂z

¶∼= 1h2{aS(fi,j+1 − fi,j)− aN(fi,j − fi,j−1)}, (2.21)

onde aN e aS são dados por:

aS =hR (i,j+1)

(i,j)

1

adx=

hRS

1

adx, (2.22)

19

e

aN =hR (i,j)

(i,j−1)1

adx=

hRN

1

adx, (2.23)

e S e N denotam integrações ao longo dos segmentos de reta verticais de (i, j) a (i, j+1)

e (i, j − 1) a (i, j), respectivamente, Figura 2.3.

Com relação às derivadas mistas,

∂

∂z

µa∂f

∂x

¶,

∂

∂x

µa∂f

∂z

¶,

usaremos duas aproximações. A primeira, que será referida por forma curta, será usada

no cálculo do campo de deslocamento no interior do modelo. A segunda, desenvolvida

por forma cheia, será usada na superfície (j=1), em conjunto com o chamado Formalismo

do Vácuo, como sugerido em Zahradník e Priolo (1994).

Inicialmente, será deduzido a forma curta da aproximação usada para,

∂

∂z

µa∂f

∂x

¶.

Assim, mais uma vez usando

g ∼= a∂f∂x

20

e

∂g

∂z∼= g(i, j + 1/2)− g(i, j − 1/2)

h, (2.24)

onde

g(i, j ± 12) ∼= a∂f

∂x|(i,j±1

2), (2.25)

g(i, j + 1/2) ∼=R (i,j+1)(i,j)

∂f

∂xdzR (i,j+1)

(i,j)

1

adz

(2.26)

e

g(i, j − 1/2) 'R (i,j)(i,j−1)

∂f

∂xdzR (i,j)

(i,j−1)1

adz, (2.27)

ou ainda, usando as definições de aS e aN ,

g(i, j + 1/2) ' aSh

Z (i,j+1)(i,j)

∂f

∂xdz (2.28)

e

g(i, j − 1/2) ' aNh

Z (i,j+1)(i,j)

∂f

∂xdz, (2.29)

21

ou

g(i, j + 1/2) ' aSh

∂f

∂x|(i,j+1/2) h (2.30)

e

g(i, j − 1/2) ' aNh

∂f

∂x|(i,j−1/2) h. (2.31)

Assim,

g(i, j + 1/2) ' aS fi+1,j+1/2 − fi−1,j+1/22h

(2.32)

e

g(i, j − 1/2) ' aN fi+1,j−1/2 − fi−1,j−1/22h

. (2.33)

Então, usando as aproximações

fi±1,j+1/2 ' fi±1,j + fi±1,j+12

(2.34)

e

fi±1,j−1/2 ' fi±1,j + fi±1,j−12

, (2.35)

22

chega-se

g(i, j + 1/2) ' aS4h[fi+1,j+1 + fi+1,j − fi−1,j+1 − fi−1,j] , (2.36)

g(i, j − 1/2) ' aN4h[fi+1,j + fi+1,j−1 − fi−1,j − fi−1,j−1] , (2.37)

portanto

∂

∂z

µa∂f

∂x

¶' 14h2

{aS[fi+1,j+1 + fi+1,j − fi−1,j+1 − fi−1,j]

− aN [fi+1,j + fi+1,j−1 − fi−1,j − fi−1,j−1]} . (2.38)

Analogamente, a forma curta de aproximações de ∂∂x

¡a∂f∂z

¢é dada por:

∂

∂x

µa∂u

∂z

¶' 14h2

{aE[fi+1,j+1 + fi,j+1 − fi+1,j−1 − fi,j−1]

− aW [fi,j+1 + fi−1,j+1 − fi,j−1 − fi−1,j−1]} . (2.39)

A forma cheia de aproximação da derivada ∂∂z

¡a∂f∂x

¢é obtida de maneira similar, contudo

a aproximação inicial é um pouco diferente, ou seja,

∂g

∂z' 1

2h[g(i− 1/2, j + 1/2) + g(i+ 1/2, j + 1/2)]

−[g(i− 1/2, j − 1/2) + g(i+ 1/2, j − 1/2)] , (2.40)

onde

23

g(i+ 1/2, j + 1/2) = a∂f

∂x|(i+1/2,j+1/2) (2.41)

que é aproximada por

g(i+ 1/2, j + 1/2) 'R (i+1,j+1/2)(i,j+1/2)

∂f

∂xdxR (i+1,j+1/2)

(i,j+1/2)

1

adx

=aSEh

Z (i+1,j+1/2)(i,j+1/2)

∂f

∂xdx (2.42)

onde aSE é definido por

aSE ' hR (i+1,j+1/2)(i,j+1/2)

1

adx=

hRSE

1

adx

(2.43)

e SE indica que a integração de1

aé realizada no segmento de reta que vai de (i, j +1/2)

a (i+ 1, j + 1/2), como mostrado na Figura 2.4.

Assim sendo,

g(i+ 1/2, j + 1/2) ' aSEh[fi+1,j+1/2 − fi,j+1/2] , (2.44)

ou ainda

g(i+ 1/2, j + 1/2) ' aSE2h[fi+1,j+1 + fi+1,j − fi,j+1 − fi,j] . (2.45)

24

(i+1/2,j+1/2)

(i,j)

SE(i+1/2,j+1/2)

(i,j)

SE

Figura 2.4: Modelo da malha do método PS2 de diferenças finitas para aproximação dasderivadas em relação às variáveis x e z para pontos da malha na superfície SE.

De forma análoga, encontramos:

g(i− 1/2, j + 1/2) ' aSW2h

[fi,j+1 + fi,j − fi−1,j+1 − fi−1,j] , (2.46)

onde

aSW ' hR (i,j+1/2)(i−1,j+1/2)

1

adx=

hRSW

1

adx

(2.47)

em que SW indica que a integração é efetuada no segmento de reta que vai de (i−1, j+1/2)

a (i, j + 1/2) como mostrado na Figura 2.5.

g(i+ 1/2, j − 1/2) ' aEN2h[fi+1,j + fi+1,j−1 − fi,j − fi,j−1] , (2.48)

25

(i-1/2,j+1/2)

(i,j)

SW(i-1/2,j+1/2)

(i,j)

SW

Figura 2.5: Modelo da malha do método PS2 de diferenças finitas para aproximação dasderivadas em relação às variáveis x e z para pontos da malha na superfície SW.

e

g(i− 1/2, j − 1/2) ' aWN2h

[fi,j + fi,j−1 − fi−1,j − fi−1,j−1] , (2.49)

onde

aNE ' hR (i+1,j−1/2)(i,j−1/2)

1

adx=

hRNE

1

adx, (2.50)

e

aNW ' hR (i,j−1/2)(i−1,j−1/2)

1

adx=

hRNW

1

adx, (2.51)

onde NE e NW indicam que as integrações são realizadas nos segmentos apresentados

na Figura 2.6.

26

(i-1/2,j-1/2)

(i,j)

WN EN(i-1/2,j-1/2)

(i,j)

WN EN

Figura 2.6: Modelo da malha do método PS2 de diferenças finitas para aproximação dasderivadas em relação às variáveis x e z para pontos da malha na superfície NE e NW.

Com isso, a forma cheia de ∂∂z

¡a∂f∂x

¢escreve-se:

∂

∂z

µa∂f

∂x

¶' 1

4h2{aSE(fi+1,j+1 + fi+1,j − fi,j+1 − fi,j)

+aSW (fi,j+1 + fi,j − fi−1,j+1 − fi−1,j)

−aNE(fi+1,j + fi+1,j−1 − fi,j − fi,j−1)

−aNW (fi,j + fi,j−1 − fi−1,j − fi−1,j−1)}. (2.52)

Finalmente, de forma análoga pode-se escrever:

∂

∂x

µa∂f

∂z

¶' 1

4h2{aES(fi+1,j+1 + fi,j+1 − fi+1,j − fi,j)

+aEN(fi+1,j + fi,j − fi+1,j−1 − fi,j−1)

−aWS(fi,j+1 + fi−1,j+1 − fi,j − fi−1,j)

−aWN(fi,j + fi−1,j − fi,j−1 − fi−1,j−1)} (2.53)

27

(i,j)

WN EN

WS ES

(i,j)

WN EN

WS ES

Figura 2.7: Modelo da malha do método PS2 de diferenças finitas para aproximação dasderivadas em relação às variáveis x e z para pontos da malha no interior WN, WS, EN eES.

ES, EN , WS e WN denotam os segmentos de reta usados nas integrações que definem

os parâmetros aES, aEN , aWS e aWN , respectivamente, Figura 2.7.

aES ' hR (i+1/2,j+1)(i+1/2,j)

1

adz, aEN ' hR (i+1/2,j)

(i+1/2,j−1)1

adz

(2.54)

aWS ' hR (i−1/2,j+1)(i−1/2,j)

1

adz, aWN ' hR (i−1/2,j)

(i−1/2,j−1)1

adz

(2.55)

Agora podemos escrever ás Equações 2.8 e 2.9 para o cálculo das com-

ponemtes do campo de onda de deslocamento em todos os pontos da malha, especifica-

mente no interior (i = 2,...., Nx-1) e (j = 2,....,Nz-1) e na superfície (borda superior, ou

seja, i = 2, ....Nx-1 e j = 1).

Assim, tendo em conta que $n1(i, j) e $n2(i, j) representam a soma das

28

derivadas parciais espaciais nas Equações 2.8 e 2.9, respectivamente, ou seja:

$n1(i, j) =∂

∂x

µC11

∂u

∂x

¶(i, j) +

∂

∂x

µC13

∂v

∂z

¶(i, j)

+∂

∂z

µC44

∂v

∂x

¶(i, j) +

∂

∂z

µC44

∂u

∂z

¶(i, j) (2.56)

e

$n2(i, j) =∂

∂x

µC44

∂v

∂x

¶(i, j) +

∂

∂x

µC44

∂u

∂z

¶(i, j)

+∂

∂z

µC11

∂u

∂x

¶(i, j) +

∂

∂z

µC33

∂v

∂z

¶(i, j) , (2.57)

onde usando as aproximações das derivadas encontradas acima (usando a forma curta

para as derivadas mistas), temos:

∂

∂x

µC11

∂u

∂x

¶n(i, j) ' 1

h2{C11E(i, j)[ui+1,j − ui,j]

−C11w(i, j)[ui,j − ui−1,j]} (2.58)

∂

∂x

µC13

∂v

∂z

¶n(i, j) ' 1

4h2{C13E(i, j)[vi+1,j+1 − vi,j+1

−vi+1,j−1 − vi,j−1]

−C13w(i, j)[vi,j+1 + vi−1,j+1

−vi,j−1 − vi−1,j−1] (2.59)

29

∂

∂z

µC44

∂v

∂x

¶n(i, j) ' 1

4h2{C44S(i, j)[vi+1,j+1 + vi+1,j

−vi−1,j+1 − vi−1,j]

−C44N(i, j)[vi+1,j + vi+1,j−1

−vi−1,j − vi−1,j−1] (2.60)

∂

∂z

µC44

∂u

∂z

¶n(i, j) ' 1

h2{C44S(i, j)[ui,j+1 − ui,j]

−C44N(i, j)[ui,j − ui,j−1]} (2.61)

∂

∂x

µC44

∂v

∂x

¶n(i, j) ' 1

h2{C44E(i, j)[vi+1,j − vi,j]

−C44w(i, j)[vi,j − vi−1,j]} (2.62)

∂

∂x

µC44

∂u

∂z

¶n(i, j) ' 1

4h2{C44E(i, j)[ui+1,j+1 + ui,j+1

−ui+1,j−1 − ui,j−1]

−C44w(i, j)[ui,j+1 + ui−1,j+1

−ui,j−1 − ui−1,j−1] (2.63)

30

∂

∂z

µC11

∂u

∂x

¶n(i, j) ' 1

4h2{C11S(i, j)[ui+1,j+1 + ui+1,j

−ui−1,j+1 − ui−1,j]

−C11N(i, j)[ui+1,j + ui+1,j−1

−ui−1,j − ui−1,j] (2.64)

∂

∂z

µC33

∂v

∂z

¶n(i, j) ' 1

h2{C33S(i, j)[vi,j+1 − vi,j]

−C33N(i, j)[vi,j − vi,j−1]} . (2.65)

Desta forma, as Equações 2.1 e 2.2 podem ser escritas como:

un+1i,j =DT 2

ρ(i, j)h2{C11E(i, j)(ui+1,j − ui,j)− C11E(i− 1, j)(ui,j − ui−1,j)

+0, 25[C13E(i, j)(vi+1,j+1 − vi,j+1 − vi+1,j−1 − vi,j−1)

−C13E(i− 1, j)(vi,j+1 + vi−1,j+1 − vi,j−1 − vi−1,j−1)]

+0, 25[C44S(i, j)(vi+1,j+1 + vi+1,j − vi−1,j+1 − vi−1,j)

−C44S(i, j − 1)(vi+1,j + vi+1,j−1 − vi−1,j − vi−1,j−1)]

+C44S(i, j)(ui,j+1 − ui,j)− C44S(i, j − 1)(ui,j−ui,j−1)}

+2ui,j − un−1i,j (2.66)

31

e

vn+1i,j =DT 2

ρ(i, j)h2{C44E(i, j)(vi+1,j − vi,j)− C44E(i− 1, j)(vi,j − vi−1,j)

+0, 25[C44E(i, j)(ui+1,j+1 + ui,j+1 − ui+1,j−1 − ui,j−1)

−C44E(i− 1, j)(ui,j+1 + ui−1,j+1 − ui,j−1 − ui−1,j−1)]

+0, 25[C11S(i, j)(ui+1,j+1 + ui+1,j − ui−1,j+1 − ui−1,j)

−C11S(i, j − 1)(ui+1,j + ui+1,j−1 − ui−1,j − ui−1,j−1)]

+C33S(i, j)(vi,j+1 − vi,j)− C33S(i, j − 1)(vi,j−vi,j−1)}

+2vi,j − vn−1i,j , (2.67)

para i = 2,...., Nx-1; j = 2,....Nz-1 e k = 1,...., Nt.

Usando o fato que CIJW(i, j) = CIJE(i − 1, j) e CIJN(i, j) = CIJE(i − 1, j),

reduzimos pela metade o número de parâmetros envolvidos.

Para o cálculo dos componentes uni,j, νni,j em j = 1, utilizamos os operado-

res cheios nas aproximações das derivadas mistas, como sugerido em Zahradník e Priolo

(1994), além das hipóteses relativas ao chamado Formalismo do Vácuo (Zahradník, Moczo

e Hron (1993)), ou seja, para j = 0 assumimos:

uni,j = vni,j = CIJ (i, j) = ρ(i, j) = 0 . (2.68)

32

Assim, temos:

un+1i,1 =DT 2

ρ(i, 1)h2{C11E(i, 1)(ui+1,1 − ui,1)− C11E(i− 1, 1)(ui,1 − ui−1,1)

+0, 25[C13ES(i, 1)(vi+1,2 − vi,2 − vi+1,1 − vi,1)

−C13ES(i− 1, 1)(vi,2 + vi−1,2 − vi,1 − vi−1,1)]

+C44S(i, j)(ui,2 − ui,1)

+0, 25[C44SE(i− 1, j)(vi,2 + vi,1 − vi−1,2 − vi−1,1)

−C44SE(i, 1)(vi+1,2 + vi+1,1 − vi,2 − vi,1)]

+2ui,1 − un−1i,1 (2.69)

e

vn+1i,1 =DT 2

ρ(i, 1)h2{0, 25[C44ES(i, j)(ui+1,2 + ui,2 − ui+1,1 − ui,1)

−C44ES(i− 1, 1)(ui,2 + ui−1,2 − ui,1 − ui−1,1)]

−C44E(i, j)(vi+1,1 − vi,1)− C44E(i− 1, 1)(vi,1 − vi−1,1)

+0, 25[C13SE(i− 1, 1)(ui,2 + ui,1 − ui−1,2 − ui−1,1)

+C13SE(i, 1)(ui+1,2 + ui+1,1 − ui,2 − ui,1)]

+C33S(i, j)(vi,2 − vi,j)}

+2vi,j − vi,j , (2.70)

para i = 2,....,Nx-1, e k=1,....,Nt.

33

2.1.2 Condições de Contorno para as Bordas Laterais e Inferior

Serão aqui apresentados, os operadores envolvidos na estratégia usada para

atenuação das reflexões nas bordas do modelo nas modelagens elásticas, ou seja, as

equações introduzidas por Cerjan (1985), Emerman e Stephen (1983) e as equações não

reflexivas (equações da onda one-way).

As condições de bordas de absorção propostas por Cerjan (1985) [7] com os

operadores propostos por Emerman e Stephen (1990) [10] são implementadas para solu-

cionar os problemas que mais surgem na aplicação de métodos de soluções discretizadas

para cálculo de propagação de onda, a presença de reflexões e envoltórias vindas das bor-

das da malha numérica. Estes eventos indesejados prejudicam o sinal sísmico real que se

propaga na região modelada.

Uma das soluções para evitar o efeito das bordas é aumentar a malha numéri-

ca, assim as reflexões das bordas e as envoltórias têm o mesmo tempo envolvido na mo-

delagem. Obviamente esta solução aumenta consideravelmente o custo computacional.

As condições de contorno para bordas não reflexivas, introduzidas para o

método de diferenças finitas (Clayton and Enquist (1977) [11]; Reynolds (1978) [12]) são

baseadas na troca da equação de onda numa faixa da região das bordas por uma equação

unidirecional que não permite que a energia de propagação que vem da reflexão da borda

volte para a malha numérica.

34

2.2 Termo Fonte

As fontes de energia sísmica mais utilizadas são a dinamite e o vibrador em

levantamentos terrestres e canhões de ar comprimido em levantamentos marítimos. Cada

uma destas fontes emite um pulso característico conhecido como assinatura da fonte que

se propaga em todas as direções com objetivo principal de encontrar reservatórios de óleo

e gás pelas propriedades reflexivas das rochas no interior da Terra. Estes pulsos elásticos

ou detonações são de duração ou comprimento muito pequeno, da ordem de 200 ms, e se

refletem e se refratam em cada uma das camadas geológicas em profundidade, retornando

a superfície, no caso de sísmica de superfície, ou para o poço, no caso de sísmica interpoços,

com informações valiosas para pesquisa de petróleo.

Para iniciar o processo de propagação de ondas sísmicas, foi adicionado às

equações que regem o campo de deslocamento das partículas, um termo fonte, especifica-

mente em termo da densidade de força, como sugerido em Cunha (1997) [6]:

f (t) =£1− 2π (πfct)2

¤e−π(πfct)

2

, (2.71)

realizando a discretização, substituindo t→ (n− 1)∆t− tf , fica:

f ((n− 1)∆t) = £1− 2π (πfc((n− 1)∆t− tf))2¤ e−π(πfc((n−1)∆t−tf ))2 , (2.72)onde fc, frequência central está em função da frequência de corte fcorte:

fc =fcorte3√π. (2.73)

35

0,00 0,02 0,04 0,06 0,08 0,10 0,12-1,2

-1,0

-0,8

-0,6

-0,4

-0,2

0,0

0,2

0,4

0,6

Ampl

itude

Tempo (s)

Função Fonte (Fc=60 Hz)

Figura 2.8: Gráfico da função fonte: Amplitude x Tempo

Para frequência de corte de 60 Hz e 0,00015 s de intervalo de tempo temos

o gráfico da funcão fonte apresentado na Figura 2.8.

O comprimento do pulso (Nf) e o período do pulso (tf) no domínio do

tempo são:

Nf =4√π

Dtfmax, tf =

2√π

fcorte. (2.74)

A adequada simulação da fonte requer um tratamento especial para os

casos homogêneo e heterogêneo. Os receptores utilizados para registrar as reflexões destes

pulsos são basicamente de dois tipos: eletromagnético (geofones) para registro em terra,

e de pressão (hidrofones) para levantamentos na água.

Os hidrofones utilizam cristais piezoelétricos, que geram uma corrente elétri-

36

ca proporcional à variação de pressão produzida pela onda. Estes receptores, a exemplo

dos geofones, devem reproduzir o mais fielmente possível as vibrações mecânicas na forma

de oscilações elétricas.

Estas oscilações elétricas são transmitidas até ao sismógrafo, onde são digi-

talizadas, multiplexadas e registradas (ou retransmitidas via satélite para uma central de

computadores), após severo depuramento e amplificação eletrônica [25].

2.3 Definição da malha

2.3.1 Condições para Não Ocorrência de Dispersão Numérica

As dimensões da malha são de importância vital para o método das dife-

renças finitas. As figuras mostrando os detalhes do modelo, como as dimensões totais em

ambas as direções, as velocidades em cada camada e as velocidades máximas e mínimas

para cada modelo, serão apresentadas na seção aplicações.

Considera-se que a função velocidade ci,j é discretizada dentro de um valor

médio para cada quadrado da malha. Esta hipótese é válida desde que os espaçamentos

da malha sejam pequenos comparados com o comprimento de onda da propagação.

Uma relação entre a menor velocidade utilizada no modelo (cmin) e a fre-

qüência máxima (fmax), limita o máximo valor do espaçamento da malha de forma a

não se ter excessiva dispersão de energia [31], lembrando que neste modelo utiliza-se

37

h = ∆x = ∆z:

h =cmin

ω.fmax, (2.75)

onde ω representa o número máximo de amostras por comprimento de onda correspon-

dente à freqüência máxima. O valor ótimo encontrado de maneira empírica para este

número é 5.

2.3.2 Condições de Estabilidade Numérica

Um problema muito importante que deve ser considerado é a estabilidade

numérica. Para solucionar este problema foi desenvolvida uma relação para controle

dos valores dos intervalos do tempo de amostragem para evitar que o sistema se torne

numericamente instável, onde cmax é a maior velocidade adotada no modelo e µ é uma

constante definida da mesma forma que na dispersão da malha [31]. O melhor valor

encontrado para esta constante é 5.

Dt =h

µ.cmax, (2.76)

2.4 Metodologia

A primeira metodologia é definida nos ítens abaixo:

1- Definição da malha:

• Maior valor de velocidade das ondas compressionais e menor valor de velocidade das

38

ondas cisalhantes, no modelo fornecido pelo geofísico intérprete.

• Frequência máxima da fonte sísmica (fmax).

• Obtenção de h e Dt empregando as Equações 2.75 e 2.76.

2- Entrada de Dados:

Modelo sísmico:

Considerando i = 1,...,Nx e j = 1,...,Nz:

• Campo de velocidades das ondas compressionais: Vp(i,j);

• Campo de velocidades das ondas cisalhantes: Vs(i,j);

• Campo de densidade: ρ(i,j);

• Campos dos parâmetros de Thomsen: ε(i,j) e δ(i,j);

Obs: A malha deve ser pelo menos dez vezes mais densa do que a encontrada

com a condição de dispersão, para a realização das integrações ao longo das linhas da

malha, Equação 2.75.

Fonte sísmica:

• Frequência de corte;

• Tipo de fonte (explosiva ou direcional);

• Posição da fonte em relação à malha.

39

Geometria de Aquisição:

• Número e posição dos receptores em relação à malha.

Intervalo de amostragem temporal para saída:

• Tempo de registro;

• Definição da malha.

3- Integração dos parâmetros:

• Obtenção dos parâmetros elásticos C11(i,j), C13(i,j), C33(i,j) e C44(i,j);

• Obtenção dos parâmetros C11E(i,j), C11W (i,j), C11S(i,j), C11N(i,j), C13E(i,j), C13W (i,j),

C33S(i,j), C33N(i,j), C44S(i,j), C44W (i,j), C44N(i,j) e C44E(i,j) através de integrações

ao longo das linhas horizontais e verticais;

• Obtenção dos parâmetros C13ES(i,j), C44SE(i,j);

Estes parâmetros são definidos em relação à malha definida.

4- Execução do programa principal:

Como mostra a Figura 2.9 os dados de entradas, de integração dos parâmet-

ros servem como entrada para o programa principal.

5- Saídas:

• Sismogramas relativos aos componentes vertical e horizontal do campo de desloca-

mento das partículas;

• Instantâneos e filmes apresentando a propagação dos campos de deslocamento.

40

—

Figura 2.9: Diagrama representando a 1ametodologia para modelagem em meios elásticosTIV.

41

2.5 Aplicações

2.5.1 Primeira Aplicação: Propagação de Onda em Meios Ho-

mogêneos Anisotrópicos nosModelos Propostos emThom-

sen (1986)

O objetivo desta aplicação é exibir a assinatura sísmica das frentes de onda

qP e qSV emmeios transversalmente isotrópicos com eixo de simetria vertical encontrados

na natureza. Especificamente, escolhemos 23 modelos dos apresentados em Thomsen

(1986), listados na Tablela 2.2. Os parâmetros usados nas modelagens estão descritos na

Tabela 2.1. As Figuras 2.10 a 2.21 exibem os resultados.

Tabela 2.1: Parâmetros da malha usada nas modelagens na primeira aplicação.

h 5 distância entre pontos consecutivos nas direções x e z (malha regular)Dt(s) 0,0005 tempo de amostragem de diferenças finitas

fcorte(Hz) 60 frequência de corteNx 601 número de pontos da malha na direção xNz 601 número de pontos da malha na direção z

Ntotal 2000 número total de passos de tempoixf 301 índice relativo à posição da fonte explosiva na direção xjzf 301 índice relativo à posição da fonte explosiva na direção zNa 50 número de colunas ou linhas da borda não-reflexivajobs 3 índice relativo à posição da superfície de observação (receptores)

42

Tabela 2.2: Parâmetros de Thomsen usados na primeira aplicação (obtidos em Thomsen,1986).

Amostra Vp (m/s) Vs (m/s) ε δ γArenito 3368 1829 0,110 -0,127 0,255Folhelho 4529 2703 0,034 0,250 0,046

Arenito Imaturo 4476 2814 0,097 0,051 0,051Calcáreo Síltico 4972 2899 0,056 -0,041 0,067Folhelho argiloso 3928 2055 0,334 0,818 0,575Siltito Laminado 4449 2585 0,091 0,688 0,046Arenito Calcarenoso 5460 3219 0,000 -0,345 -0,007Folhelho com Óleo 4231 2539 0,200 0,000 0,145Cinza Vulcânica 4846 1856 0,200 -0,003 0,105

Cristal 4420 2091 1,120 -1,230 2,280Cristal de Quartzo 6096 4481 -0,096 0,169 -0,159

Calcita 5334 3353 0,369 0,127 0,169Biotita 4054 1341 1,222 -1,437 6,120Apatita 6340 4389 0,097 0,257 0,079Gelo 3627 1676 -0,038 -0,100 0,031

Composto Alumínio-Lucita 2868 1350 0,970 -0,890 1,300Arenito-Folhelho 3009 1654 0,013 -0,010 0,035

Folhelho Anisotrópico-SS 3009 1654 0,059 -0,042 0,163Calcáreo-Folhelho 3306 1819 0,134 -0,094 0,156

Folhelho Anisotrópico-LS 3306 1819 0,169 -0,123 0,271Folhelho Anisotrópico 2745 1508 0,103 -0,073 0,345

Gás-Arenito-Água-Arenito 1409 780 0,022 -0,002 0,004Material de Gipso intemperizado 1911 795 1,161 -1,075 2,781

43

0 200100 300 400 500 6000

100

200

300

400

500

600

a) b)

0 200100 300 400 500 6000

100

200

300

400

500

600

Figura 2.10: Instantâneos do componente vertical no tempo de 3 ms. a) alumínio-lucitaε = 0, 970, δ = −0, 890. b) folhelho anisotrópico ε = 0, 103, δ = −0, 073.

0 200100 300 400 500 6000

100

200

300

400

500

600

a) b)

0 200100 300 400 500 6000

100

200

300

400

500

600

Figura 2.11: Instantâneos do componente vertical no tempo de 3 ms. a) apatita ε = 0, 097,δ = 0, 257. b) biotita ε = 1, 222, δ = −1, 437.

44

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.12: Instantâneos do componente vertical no tempo de 3 ms. a) calcita ε = 0, 369,δ = 0, 127. b) arenito calcarenoso ε = 0, 000, δ = −0, 345.

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.13: Instantâneos do componente vertical no tempo de 3 ms. a) folhelhoanisotrópico-ss ε = 0, 059, δ = −0, 042. b) folhelho argiloso ε = 0, 334, δ = 0, 818.

45

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.14: Instantâneos do componente vertical no tempo de 3 ms. a) calcáreo sílticoε = 0, 056, δ = −0, 041. b) arenito-folhelho ε = 0, 013, δ = −0, 010.

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.15: Instantâneos do componente vertical no tempo de 3 ms. a) arenito ε = 0, 110,δ = −0, 127. b) cristal de quartzo ε = −0, 096, δ = 0, 169.

46

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.16: Instantâneos do componente vertical no tempo de 3 ms. a) folhelho com óleoε = 0, 200, δ = 0, 000. b) folhelho ε = 0, 034, δ = 0, 250.

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.17: Instantâneos do componente vertical no tempo de 3 ms. a) folhelhoanisotrópico-ls ε = 0, 169, δ = −0, 123. b) calcáreo-folhelho ε = 0, 134, δ = −0, 094

47

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.18: Instantâneos do componente vertical no tempo de 3 ms. a) siltito laminadoε = 0, 091, δ = 0, 688. b) arenito imaturo ε = 0, 097, δ = 0, 051.

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.19: Instantâneos do componente vertical no tempo de 3 ms. a) gelo ε = −0, 038,δ = −0, 100. b) material de gipso intemperizado ε = 1, 161, δ = −1, 075.

48

0

100

200

300

400

500

600

a) b)

0 00 200100 300 400 500 6000

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.20: Instantâneos do componente vertical no tempo de 3 ms. a) gas-arenito-água-arenito ε = 0, 022, δ = −0, 002. b) cristal ε = 1, 120, δ = −1, 230.

0

100

200

300

400

500

600

0 00 200100 300 400 500 600

Figura 2.21: Instantâneos do componente vertical no tempo de 3 ms. a) cinza vulcânicaε = 0, 200, δ = −0, 003.

49

2.5.2 Segunda Aplicação: Modelo Formado por Duas Camadas

Separadas por Interface Curva

Nesta aplicação exibimos instantâneos e sismogramas encontrados em si-

mulações sísmicas realizadas nos modelos definidos na Figura 2.22, visando estimar os

efeitos nas reflexões da interface que separa os dois meios decorrentes da anisotropia. Os

parâmetros usados nas modelagens estão descritos na Tabela 2.3.

As Figuras 2.23 e 2.24 mostram instantâneos relativos a propagação do cam-

po de onda em meio homogêneo com Vp=3000 m/s e VS=1732 m/s. Neles observamos um

aumento gradual da elipsidade no que diz respeito à propagação das ondas compressionais

(quase compressionais, rigorosamente), com o aumento dos valores de ε e δ. Observamos

também um resquício de onda S originado por não usarmos uma fonte exploxiva exata.

Ao contrário da frente de ondas compressionais, as ondas cisalhantes exibem um caráter

mais isotrópico, independente dos valores de ε e δ assumidos nesta aplicação.

Tabela 2.3: Parâmetros da malha usada nas modelagens na segunda aplicação.

h 5 distância entre pontos consecutivos nas direções x e z (malha regular)Dt(s) 0,0005 tempo de amostragem de diferenças finitas

fcorte(Hz) 60 frequência de corteNx 400 número de pontos da malha na direção xNz 400 número de pontos da malha na direção z

Ntotal 2000 número total de passos de tempoixf 200 índice relativo à posição da fonte explosiva na direção xjzf 3 índice relativo à posição da fonte explosiva na direção zNa 50 número de colunas ou linhas da borda não-reflexivajobs 2 índice relativo à posição da superfície de observação (receptores)

As Figuras 2.27 e 2.28 são o resultado da subtração das ondas superficiais

50

0

Prof

undi

dade

(m)

0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900

Superfície (m)

200

400

600

800

1000

1200

1400

1600

1800

Figura 2.22: Representa os modelos sísmicos abordados na segunda aplicação. A primeiracamada apresenta VP=3000 m/s e VS= 1732 m/s; a segunda possui VP=3800 m/s eVS=2196 m/s. Nos quatro casos abordados a segunda camada é isotrópica. No primeirocaso a camada superior é isotrópica e nos demais transversalmente isotrópica com ε = δiguais a 0,1, 0,2 e 0,3, respectivamente. Em todos os casos a densidade foi consideradainvariável ao longo do modelo (2,4 g/cm3).

51

Prof

undi

dade

(m)

Superfície (m)0

0500 1000 1500 2000

500

1000

1500

2000

Componente V

Prof

undi

dade

(m)

Superfície (m)0

0500 1000 1500 2000

500

1000

1500

2000

Prof

undi

dade

(m)

Superfície (m)500 1000 1500 2000

1000

1500

2000

Prof

undi

dade

(m)

Superfície (m)500 1000 1500 2000

1000

1500

2000

00

500

00

500

Figura 2.23: Instantâneos do componente vertical do campo de onda, considerando fonteexplosiva na posição central do modelo homogêneo, VP=3000 m/s e VS=1732 m/s. a)Caso isotrópico. b) Caso anisotrópico, com ε = δ = 0, 1. c) Caso anisotrópico, comε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3. Observe o crescimento da elipsidadecom o aumento do valor dos parâmetros de Thomsen, no caso das ondas compressionais. Afrente de ondas cisalhantes (mais lenta e de menor amplitude) é decorrente da aproximaçãousada na introdução da fonte explosiva. Para os valores de anisotropia escolhidos as ondascisalhantes apresentam um caráter mais isotrópico em relação às compressionais.

52

Prof

undi

dade

(m)

Superfície (m)0

0500 1000 1500 2000

500

1000

1500

2000

Componente U

Prof

undi

dade

(m)

Superfície (m)0

0500 1000 1500 2000

500

1000

1500

2000

Prof

undi

dade

(m)

Superfície (m)500 1000 1500 2000

1000

1500

2000

Prof

undi

dade

(m)

Superfície (m)500 1000 1500 2000

1000

1500

2000

00

500

00

500

Figura 2.24: Instantâneos do componente horizontal do campo de onda, considerandofonte explosiva na posição central do modelo homogêneo, VP=3000 m/s e VS=1732 m/s.a) Caso isotrópico. b) Caso anisotrópico, com ε = δ = 0, 1. c) Caso anisotrópico, comε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3. Observe o crescimento da elipsidadecom o aumento do valor dos parâmetros de Thomsen, no caso das ondas compressionais. Afrente de ondas cisalhantes (mais lenta e de menor amplitude) é decorrente da aproximaçãousada na introdução da fonte explosiva. Para os valores de anisotropia escolhidos as ondascisalhantes apresentam um caráter mais isotrópico em relação às compressionais.

53

Tem

po d

e re

gist

ro (s

)Superfície (m)

00

500 1000 1500 2000

0.2

0.3

Componente V

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)0

0500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)0

0500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

00

Figura 2.25: Sismogramas relativos ao componente vertical do campo de ondas. Em todosos casos observa-se a predominância das ondas superficiais (ondas P e S diretas e ondasRayleigh) em relação as ondas refletidas na interface que delimita as duas camadas. a)Caso isotrópico. b) Caso anisotrópico, com ε = δ = 0, 1. c) Caso anisotrópico, comε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3.

54

Tem

po d

e re

gist

ro (s

)Superfície (m)

00

500 1000 1500 2000

0.2

0.3

Componente U

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)0

0500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)0

0500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

00

Figura 2.26: Sismogramas relativos ao componente horizontal do campo de ondas. Emtodos os casos observa-se a predominância das ondas superficiais (ondas P e S diretas eondas Rayleigh) em relação as ondas refletidas na interface que delimita as duas camadas.a) Caso isotrópico. b) Caso anisotrópico, com ε = δ = 0, 1. c) Caso anisotrópico, comε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3.

55

Tem

po d

e re

gist

ro (s

)Superfície (m)

00

500 1000 1500 2000

0.2

0.3

Componente V

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)0

0500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)0

0500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

00

Figura 2.27: Sismogramas relativos ao componente vertical do campo de ondas, sem apresença das ondas superficiais. a) Caso isotrópico. b) Caso anisotrópico, com ε = δ =0, 1. c) Caso anisotrópico, com ε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3.Na parte (a) indicamos os quatro modos de onda refletida presente nesta modelagem.Observamos variações consideráveis nos tempos de trânsito das ondas PP nos offsetslongos. As ondas SP também apresentam variações, porém inferiores às apresentadas nasondas PP. As ondas PS e SS praticamente não apresentam variações.

56

Tem

po d

e re

gist

ro (s

)Superfície (m)

00

500 1000 1500 2000

0.2

0.3

Componente U

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)0

0500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)0

0500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

Tem

po d

e re

gist

ro (s

)

Superfície (m)500 1000 1500 2000

0.2

0.3

0.1

0.4

0.5

0.6

0.7

0.8

00

Figura 2.28: Sismogramas relativos ao componente horizontal do campo de ondas, sem apresença das ondas superficiais. a) Caso isotrópico. b) Caso anisotrópico, com ε = δ =0, 1. c) Caso anisotrópico, com ε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3.Na parte (a) indicamos os quatro modos de onda refletida presente nesta modelagem.Observamos variações consideráveis nos tempos de trânsito das ondas PP nos offsetslongos. As ondas SP também apresentam variações, porém inferiores às apresentadas nasondas PP. As ondas PS e SS praticamente não apresentam variações.

57

dos sismogramas completos exibidos nas Figuras 2.25 e 2.26. Na parte (a) das figuras

indicamos os modos de onda relativos às reflexões. Observamos, como esperado, variações

consideráveis no tempo de trânsito das ondas PP nos offsets maiores. As ondas SP

também são modificadas com o valor crescente de ε e δ. As ondas PS e SS praticamente

não são alteradas com a variação dos parâmetros de Thomsen usados nesta aplicação.

2.5.3 Terceira Aplicação: Modelo com Altos Contrastes de Im-

pedância

Nesta aplicação mostramos os instantâneos e sismogramas obtidos na simu-

lação sísmica realizada num modelo que apresenta alto grau de complexidade estrutural,

descrito nas Figuras 2.29 e 2.30. O objetivo é mostrar a capacidade do método de abor-

dar modelos realistas, além de exibir os efeitos decorrentes da anisotropia. Os parâmetros

usados nas modelagens estão dispostos na Tabela 2.4 e os resultados são apresentados nas

Figuras 2.31 até 2.35.

Nas figuras relativas aos instantâneos, Figuras 2.31 e 2.32, observamos um

espalhamento lateral de energia superior na propagação na camada anisotrópica, também

constatado nos maiores offsets nas Figuras 2.33 e 2.34 as partes (a) e (b) destas figuras e-

xibem os simogramas completos relativos aos componentes vertical e horizontal do campo

de ondas, respectivamente; enquanto que as partes (c) e (d) mostram os mesmos sismo-

gramas subtraídos das superficiais (ondas P e S diretas e ondas Rayleigh). Nestas últimas

observamos o maior espalhamento lateral de energia no caso anisotrópico, como indicado

nos instantâneos. Tal fato tem o potencial de degradar o imageamento sísmico quando

58

Prof

undi

dade

(m)

Superfície (m)

00

250

500

750

1000

1250

1500

1750

2000

2250

2500

2750

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

Vp= 1700 m/s

Vp= 2500 m/s

Vp= 3200 m/s

Vp= 3700 m/s

Vp= 4700 m/s

Vp= 4200 m/s

Vp= 5200 m/s

Vp= 5000 m/s

Vp= 4000 m/s

Vp= 5100 m/sVp= 6345 m/s

Figura 2.29: Modelo de velocidades da bacia sedimentar terrestre estudada com os valoresde velocidade da onda P em cada camada. Os valores de velocidade da onda S podemser obtidos dividindo nos valores de VP por

√3. O valor de densidade é considerado

constante ao longo do modelo, igual a 2,4 g/cm3.

Vp=2500 m/sEpsilon = delta = 0,2

Prof

undi

dade

(m)

Superfície (m)

00

250

500

750

1000

1250

1500

1750

2000

2250

2500

2750

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

Figura 2.30: Distribuição dos parâmetros de Thomsen no modelo para o caso terrestre.Com exceção da camada que corresponde à velocidade das ondas compressionais igual a2500 m/s, na qual ε = δ = 0, 2, temos valores nulos para estes parâmetros.

59

o caráter anisotrópico da propagação do campo de onda não for levado em consideração

na estratégia de processamento dos dados, como mencionado em Rosa Filho (2002) [29].

A Figura 2.35 mostra em detalhe a onda tipo Rayleigh encontrada com a aplicação do

Formalismo do Vácuo na superfície que delimita o modelo na parte superior. Observe a

dispersão características na parte inferior da figura. Com relação ao tempo de processa-

mento gasto na simulação de um tiro no modelo em questão, utilizando um computador

pessoal com processador Pentium IV com velocidade de 1.8 GHz e 1.0 Gb de memória

RAM, foram necessários 146 minutos de processamento.

Finalmente, embora tenhamos considerado geometrias de aquisição com