TOMOGRAFIA DE ATENUACÃO ATRAVÉS DA DECOMPOSICÃO EM … · pada em duas classes: tomogra a cinemática, que utiliza como dado de entrada os tempos de trânsito entre cada par fonte

UNIVERSIDADE FEDERAL DA BAHIA

INSTITUTO DE GEOCIÊNCIAS

CURSO DE GRADUAÇÃO EM GEOFÍSICA

GEO213 � TRABALHO DE GRADUAÇÃO

TOMOGRAFIA DE ATENUACÃO

ATRAVÉS DA DECOMPOSICÃO EM

VALORES SINGULARES

CAIO MANOEL LIRA DA COSTA FONTES

SALVADOR � BAHIA

FEVEREIRO � 2014

Tomogra�a de Atenuação através da Decomposição em Valores Singulares

por

Caio Manoel Lira da Costa Fontes

GEO213 � TRABALHO DE GRADUAÇÃO

Departamento de Geofísica

do

Instituto de Geociências

da

Universidade Federal da Bahia

Comissão Examinadora

Dr. Amin Bassrei - Orientador

Dr. Carlos da Silva Vilar

Dr. Thierry Jacques Lemaire

Data da aprovação: 14/02/2014

O importante não é

vencer todos os dias,

mas lutar sempre.

RESUMO

Os problemas inversos da Geofísica geralmente são formulados como um sistema de

equações lineares. Tais problemas inversos são, na maioria dos casos, mal-postos, pois sua

solução pode não ser única, gerando ambiguidade, ou até mesmo a solução pode não existir.

Um método de imageamento originário da Medicina, bastante conhecido em outras

áreas, inclusive em Geofísica de Reservatórios é a tomogra�a. A tomogra�a pode ser agru-

pada em duas classes: tomogra�a cinemática, que utiliza como dado de entrada os tempos

de trânsito entre cada par fonte receptor, e tomogra�a dinâmica, que utiliza a forma de onda

registrada em cada receptor como dado de entrada.

Neste Trabalho de Graduação é utilizada a tomogra�a de atenuação e para o processo

inverso propriamente dito é utilizada a técnica de decomposição em valores singulares ou

SVD do inglês singular value decomposition.

No contexto deste Trabalho a tomogra�a de atenuação se encaixa como um problema

não-linear tal qual a maioria das aplicações de tomogra�a. Foi utilizada a abordagem linear

no processo inverso, onde a forma de onda é necessária para se calcular a distribuição do fator

de atenuação α, de modo que a tomogra�a de atenuação pode ser classi�cada no contexto

dinâmico. Na modelagem direta calcula-se a amplitude da onda que irá decair durante a

propagação, e para tanto é necessário o conhecimento da atenuação do meio.

Foram realizadas simulações em dois conjuntos de dados sintéticos. Com o propósito de

avaliar a robustez do método foi inserido ruído aleatório aos valores de atenuação observados

nos receptores. Também foi realizado um estudo para diferentes níveis de truncamento de

valores singulares na montagem da matriz inversa generalizada, uma vez os pequenos valores

singulares comprometem a qualidade da solução. As simulações realizadas com resultados

satisfatórios permitiram validar a metodologia em questão. O controle dos experimentos

numéricos foi através do cálculo do erro entre parâmetro de modelo verdadeiro e o parâmetro

de modelo estimado.

iii

ABSTRACT

Geophysical inverse problems are usually formulated as a system of linear equations.

Such inverse problems are ill-posed in most cases, since the solution can be non-unique,

creating ambiguity, or the solution may even not exist.

A well-known method now in other areas, including Reservoir Geophysics, tomography

was originally developed in Medicine. Tomography can be grouped into two classes: kine-

matic tomography, which uses as input data the transit times between each source receiver

pair, and dynamic tomography, which uses the waveform recorded in each receiving as input

data.

In this Graduation Work attenuation tomography is studied and the inverse process

itself is performed SVD or singular value decomposition.

In the context of this work attenuation tomography is de�ned as a nonlinear problem,

like most tomography applications. The linear approach in the inverse process was used,

where the waveform is required to calculate the distribution of attenuation factor α, so that

attenuation tomography can be classi�ed in the dynamic context. In the forward model-

ing one calculates the wave amplitude which will decay during propagation, so that the

knowledge of the medium attenuation is necessary.

There were performed simulations in two sets of synthetic data. In order to evaluate

the robustness of the method random noise was added to the attenuation data observed in

the receptors. Also it was made a study for di�erent levels of truncation of singular values

when building the generalized inverse matrix, since small singular values may compromise

the quality of the solution. The performed simulations with satisfactory results allowed

the validation of the proposed methodology. The control of the numerical experiments was

made by calculating the error between the true model parameter and the estimated model

parameter.

iv

ÍNDICE

RESUMO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii

ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv

ÍNDICE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

ÍNDICE DE TABELAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

ÍNDICE DE FIGURAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

CAPÍTULO 1 Teoria da Inversão . . . . . . . . . . . . . . . . . . . . . . . 3

1.1 Formulação do Problema Inverso . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2 Sistemas Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2.1 Estudo da Solução dos Sistemas Lineares . . . . . . . . . . . . . . . . 5

1.3 Problemas Inversos e as Questões de Condicionamento . . . . . . . . . . . . 5

1.3.1 Existência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.3.2 Unicidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.3.3 Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.4 Solução de Sistemas Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.4.1 Método dos Mínimos Quadrados . . . . . . . . . . . . . . . . . . . . . 8

1.4.2 Método dos Mínimos Quadrados Amortecidos . . . . . . . . . . . . . 9

1.4.3 Decomposição de Valores Singulares . . . . . . . . . . . . . . . . . . . 9

CAPÍTULO 2 Geotomogra�a . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.1 Tomogra�a Sísmica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.1.1 Tomogra�a de Tempos de Trânsito . . . . . . . . . . . . . . . . . . . 13

2.1.2 Tomogra�a de Atenuação . . . . . . . . . . . . . . . . . . . . . . . . 15

2.2 Modelagem Tomográ�ca Direta . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.2.1 Traçado de Raios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.2.2 Discretização do Meio . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.2.3 Ligação entre fonte e receptor . . . . . . . . . . . . . . . . . . . . . . 23

2.2.4 Modelagem Direta da Tomogra�a de Atenuação . . . . . . . . . . . . 25

2.3 Modelagem Tomográ�ca Inversa . . . . . . . . . . . . . . . . . . . . . . . . . 27

v

CAPÍTULO 3 Simulações e Resultados . . . . . . . . . . . . . . . . . . . . 28

3.0.1 Critérios de Seleção de Valores Singulares . . . . . . . . . . . . . . . 30

3.1 Modelo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.1.1 β = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.1.2 β = 10−3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.1.3 β = 10−1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.2 Modelo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.2.1 β = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.2.2 β = 10−3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.2.3 β = 10−1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

CAPÍTULO 4 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Agradecimentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Referências Bibliográ�cas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

vi

ÍNDICE DE TABELAS

3.1 Tabela de erro RMS dos parâmetros do dado e do modelo para diferentes

niveis de corte, sem a presença de ruído, no 'Modelo 1' . . . . . . . . . . . . 38


niveis de corte, com o nível de ruído β = 10−3, no 'Modelo 1' . . . . . . . . . 42




niveis de corte, sem a presença de ruído, no 'Modelo 2' . . . . . . . . . . . . 51





vii

ÍNDICE DE FIGURAS

2.1 Campo discretizado e interpolação bilinear dos índices de refração, onde cada

cruzamento da retícula corresponde a um índice de refração das coordenadas

respectivas e P (∆x,∆z) representa o ponto interpolado (Terra, 2007) . . . . 18

2.2 Interpretação aproximada da in�uência da curvatura, d2r/dl2, do ângulo ins-

tantâneo, α, entre o vetor tangente, dr/dl, e o eixo horizontal, x, e do novo

passo, r(l + ∆l). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.3 Exemplo de um modelo discretizado . . . . . . . . . . . . . . . . . . . . . . . 23

2.4 Esquema do método do canhoneio. . . . . . . . . . . . . . . . . . . . . . . . 24

2.5 Esquema do método do ligação. . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.1 Fluxograma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.2 Vetor de dados, em (a) sem in�uência de ruído e em (b) sob in�uência de

ruído da ordem de β = 10−1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.3 'Modelo 1' verdadeiro, sendo a escala em cores os valores do coe�ciente de

atenuação em m−1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.4 Grá�co (a) da amplitude e (b) da derivada do valor singular do 'Modelo 1',

em escala monolog. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.5 'Modelo 1' recuperado sem ruído onde o valor do σcorte = 103. A escala em

cores representa o valor do coe�ciente de atenuação (α) em m−1. Em (a)

a faixa dos α varia de acordo com o modelo recuperado e em (b) varia em

comparação com o modelo verdadeiro. . . . . . . . . . . . . . . . . . . . . . 35





3.7 'Modelo 1' recuperado sem ruído onde o valor do σcorte = 10−3. A escala

em cores representa o valor do coe�ciente de atenuação (α) em m−1. Em (a)







viii





3.10 'Modelo 1' recuperado, com ruído de amplitude β = 10−3, onde o valor do

σcorte = 103. A escala em cores representa o valor do coe�ciente de atenuação

(α) em m−1. Em (a) a faixa de α varia de acordo com o modelo recuperado

e em (b) varia em comparação com o modelo verdadeiro. . . . . . . . . . . . 39






σcorte = 10−3. A escala em cores representa o valor do coe�ciente de atenuação























3.18 'Modelo 2' verdadeiro, sendo a escala em cores os valores do coe�ciente de

atenuação em m−1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

ix

3.19 Grá�co (a) da amplitude e (b) da derivada do valor singular do 'Modelo 2',

em escala monolog. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47



a faixa de α varia de acordo com o modelo recuperado e em (b) varia em


































x



(α) em m−1. Em (a) a faixa dos α varia de acordo com o modelo recuperado










xi

INTRODUÇÃO

A tomogra�a é um processo de imageamento de grande utilidade em diversas áreas da

ciência que relaciona os dados obtidos em campo e os parâmetros físicos da subsuperfície.

Ela é uma técnica consagrada na medicina e desde a década de 1980 vem sendo utilizada na

Geofísica, incorporada na sísmica como método de inversão de dados e conhecida pelo seu

alto poder de resolução, quando comparada com a sísmica convencional. A principal idéia da

tomogra�a é a determinação de características desconhecidas de um meio através do estudo

de sua resposta quando o mesmo é perturbado por algum tipo de energia com intensidade

conhecida. Nos útimos anos a geotomogra�a tem sido bastante utilizada na construção

das imagens em subsuperfícies, em particular na caracterização e no monitoramento de

reservatórios, isso se deve pelo fato da mesma tratar-se de uma ferramenta que não apresenta

restrições quanto a complexidade das estruturas geológicas.

A tomogra�a de atenuação, utilizada neste trabalho, emprega como dados de entrada os

coe�cientes de atenuação medidos entre fontes e receptores, e a saída da inversão é o campo

de atenuações do meio bi-dimensional discretizado.

Na inversão desses dados geofísicos recorremos a um processo de inversão matricial

que foi realizado utilizando a técnica da decomposição em valores singulares, ou SVD (do

inglês, singular value decomposition), sendo a matriz a ser invertida denominada de matriz

tomográ�ca. Contudo, os problemas inversos são mal-postos que é o equivalente a dizer que

a solução do problema não existe, não é única e/ou não depende continuamente dos dados

de entrada. O problema inverso se deve ao fato desta matriz ser mal-condicionada, pois as

soluções obtidas são sensíveis aos dados de entrada, ou seja, uma pequena perturbação nos

dados não implica necessariamente numa pequena perturbação do modelo, sendo assim o

problema é considerado mal-posto. Os problemas geofísicos não são lineares, ou seja, algum

parâmetro necessário para se realizar a inversão do problema é função de algum parâmetro

que se deseja determinar, logo, a relação entre os dados de entrada e os parâmetros a serem

estimados é não-linear. Torna-se, então, imprescindível linearizar o problema com o propósito

de tirar essa condição de não linearidade.

Para contornar este problema, utilizamos uma quantidade ótima de valores singulares

que elimina os valores de baixa amplitude a �m de se estimar um modelo o mais próximo

possível da realidade, pois valores singulares pequenos comprometem a solução do sistema.

O presente Trabalho de Graduação está organizado da seguinte forma:

1

2

No Capítulo 1 são apresentados os problemas inversos referentes à teoria da inverão.

No Capítulo 2 são apresentados os fundamentos teóricos da tomogra�a de atenuação,

é abordada a fundamentação da teoria do raio, descrevendo os tipos de traçado de raios

utilizados, bem como a descrição dos arranjos de fonte-receptor utilizados

No Capítulo 3 serão apresentados e discutidos os resultados e simulações desse Trabalho.

No Capítulo 4 são apresentadas as principais conclusões obtidas com o Trabalho.

CAPÍTULO 1

Teoria da Inversão

A inversão é a técnica empregada que estima os parâmetros de um determinado modelo

utilizando os dados medidos com o dado de entrada.

A técnica da modelagem direta calcula os dados decorrentes de um determinado fenô-

meno modelado, conhecendo-se os parâmetros do modelo, ou seja, parte-se dos parâmetros

para a obtenção dos dados.

O problema inverso é relativamente mais complicado, uma vez que, em situações reais

podemos obter in�nitos modelos que se adequam aos mesmos dados. O objetivo do problema

inverso é determinar qualquer entrada ou o sistema que causa as medidas da saída.

Uma vez obtido os parâmetros do modelo (estimado), utilizamos os parâmetros de dados

para analisar as propriedades deste modelo e o que ele preserva do modelo real, tal como

erro e ruídos.

1.1 Formulação do Problema Inverso

Na analise de dados geofísico o ponto inicial é descreve-los, onde uma forma prática de

representar esses valores seria através de um vetor. O problema inverso pode ser representado

como:

d = Gm, (1.1)

sendo

d = [d1, d2, d3, d4, ..., dM ]T , (1.2)

e

m = [m1,m2,m3,m4, ...,mN ]T , (1.3)

onde d representa um vetor dos dados observados, m representa um vetor dos parâmetros

do modelo e G é uma matriz M×N de coe�cientes que relaciona os M dados observados aos

3

4

N parâmetros do modelo. G pode ser representada como:

G =

g11 g12 g13 g14 · · · g1N

g21 g22 g23 g24 · · · g2N...

......

.... . .

...

gM1 gM2 gM3 gM4 · · · gMN

. (1.4)

Na maior parte dos problemas geofísicos, o vetorm de parâmetros do modelo é estimado

a partir de dados observados e a matriz G é uma aproximação de um operador g que pode

ser não linear:

dcalc = g(mest), (1.5)

tal que dcalc é um vetor de dados calculados com parâmetros estimados. Supondo a matriz

G conhecida, pode-se resolver o sistema utilizando-se uma matriz inversa, ou seja, se:

d = Gm, (1.6)

então

m = G−1d. (1.7)

Deste modo, pode-se estimar os parâmetros de um modelo utilizando-se dados obser-

vados, realizando-se assim uma inversão de dados. No entanto, a matriz G só admite a

inversa G−1 se for quadrada e não-singular, ou seja, com posto completo. Como este caso é

muito raro e especí�co em problemas geofísicos reais, tornam-se necessários procedimentos

de resolução para matrizes não-quadradas ou com posto incompleto. Para o caso em que o

sistema tem mais equações do que incógnitas, utiliza-se o método dos mínimos quadrados,

obtendo-se a solução que minimiza o quadrado do erro. O método que será utilizado para a

resolução de problemas inversos é o da decomposição por valores singulares.

1.2 Sistemas Lineares

Na formulção do problema inverso, visto anteriormente, a equação explicita linear (1.6), na

notação matricial, pode ser reescrita como um sistema de equações lineares com M equações

e N incógnitas: g11m1 + g12m2 + · · · + g1NmN = d1

g21m1 + g22m2 + · · · + g2NmN = d2... +

... +. . . +

... =...

gM1m1 + gM2m2 + · · · + gMNmN = dM

(1.8)

Cujos elementos gij são os coe�cientes do sistema acima, podendo assumir valores complexos.

Resolver o sistema signi�ca encontrar os valores das incógnitas que satisfazem simultanea-

mente todas as equações, ou seja, um vetor representado por:

5

mest = [m1,m2, ...,mN ]T , (1.9)

que satisfaça o sistema acima. Tal vetor é chamado de solução do sistema linear.

1.2.1 Estudo da Solução dos Sistemas Lineares

A classi�cação do problema linear está baseado no fornecimento su�ciente de informação

para a determinação dos parâmetros do modelo, ou incógnitas do sistema. Sendo G a matriz

M ×N , o problema será:

• Subdeterminado - quando não prover informação su�ciente para determinar os pa-

râmetros do modelo, ou seja, os problemas indeterminados ocorrem quando existem

mais incógnitas do que dados, isto é M < N . Neste caso, existem várias soluções que

satisfazem o sistema;

• Determinado - quando existe informação su�ciente e exata, isto é M = N , temos,

então, única solução;

• Sobredeterminado - quando tem mais dados que incógnitas, isto é M > N .

No nosso caso particular, em se tratando de tomogra�a, é comum mais equações do que

incógnitas, o que nos leva a um sistema do tipo sobredeterminado, e estas as vezes ainda

serem incompatíveis devido a erros de arredondamentos inerentes a medição, tornando ainda

mais comum a aplicação de métodos como o Método dos Mínimos Quadrados ou de SVD.

1.3 Problemas Inversos e as Questões de Condicionamento

O problema inverso é considerado bem-posto se satisfaz as condições de existência, unici-

dade e estabilidade; e considerado mal-posto se alguma destas não seja satisfeita, ou seja, o

problema é considerado mal-posto se sua solução não existe, não é única e/ou não depende

continuamente dos dados de entrada.

1.3.1 Existência

O problema deve ter uma solução.

Essa condição pode ser violada por equações inconsistentes entre si, podendo ser aten-

dida através de uma reformulação do problema tal como Mínimos Quadrados ou SVD.

6

1.3.2 Unicidade

Deve existir apenas uma solução para o problema.

Essa condição é mais crítica, podendo também ser atendida utilizando uma reformulação

do problema, tipicamente incluindo requisitos adicionais ao problema, tal como buscar uma

solução de norma mínima. Se estes requesitos forem escolhidos adequadamente, obtém-se a

unicidade da solução.

1.3.3 Estabilidade

A solução deve variar continuamente com os dados.

Essa condição é mais difícil de ser atendida pela modi�cação de problemas originalmente

mal-postos, porque a violação da mesma implica no fato de que pequenas perturbações nos

dados podem produzir grandes perturbações nas soluções obtidas.

Para atender esta condição, faz-se necessário reformular o problema de modo a se obter

um novo problema que seja menos sensível às perturbações nos dados. Esta reformulação é

denominada condicionamento, estabilização ou regularização do problema e pode ser obtida,

por exemplo, através de requesitos adicionais de suavidade da solução. A regularização é

de grande relevância em problemas inversos, uma vez que os mesmos são frequentemente

mal-postos, requerendo regularização para que soluções realistas sejam obtidas.

Número de Condição de uma Matriz

Partindo da equação (1.6), ao tentar resolver problemas desse tipo, podemos ter problemas de

condicionamento e de instabilidade numérica. Os problemas de estabilidade numérica estão

relacionadas com o algoritmo que usamos para resolver o sistema. No entanto para problemas

mal condicionados, o sistema será sempre numericamente estável, então é interessante o fato

de identi�car quais sistemas poderão nos trazer problemas de condicionamento, ou seja, se

pegarmos a expressão anterior onde G e d são respectivamente uma matriz e um vetor, e m

o vetor solução do sistema, então aplicando o operador G−1 a esquerda e a direita, temos:

G−1d = G−1Gm, (1.10)

o que resultará em:

m = G−1d. (1.11)

Considerando que m+ δm representa a solução de um sistema perturbado, temos:

7

d + δd = G(m + δm). (1.12)

A partir de d = Gm e δm = G−1d, pode-se deduzir pela Desigualdade de Shwartz que:

||d|| ≤ ||G||||m||, (1.13)

e

||δm|| ≤ ||G−1||||δd||. (1.14)

Correlacionando as duas expressões, obtém-se que:

||δm|||m||

≤ ||G||||G−1|| ||δd|||d||

. (1.15)

O produto das normas matriciais da equação anterior trata da de�nição do número de condi-

ção. Além disso, temos que ||δd||/||d|| e ||δm||/||m|| tratam do erro relativo entre os dados

e os modelos respectivamente, logo podemos reescrever a equação anterior como:

Emodelo ≤ NC Edados, (1.16)

onde Emodelo e Edados são os erros relativos entre os parâmetros do modelo e dos dados,

respectivamente. Portanto o erro entre os parâmetros do modelo será controlado pelo valor

do NC considerado para o sistema e da qualidade dos dados de entrada.

Uma forma de se contornar este problema é encontrando um novo sistema em que as

soluções sejam menos sensíveis às perturbações nos dados, ou seja, um sistema em que o

valor de NC seja o menor possível.

1.4 Solução de Sistemas Lineares

A tomogra�a sísmica, frequentemente recai na solução de grandes sistemas lineares. Estes

sistemas podem exceder a memória do computador e requerem um elevado custo computa-

cional para serem resolvidos. Métodos para resolver estes problemas, que serão discutidos a

seguir, são de grande importância prática para a inversão tomográ�ca.

8

1.4.1 Método dos Mínimos Quadrados

O método clássico dos mínimos quadrados é aplicado normalmente a sistemas cujo número

de equações é maior do que o número de incognitas, sistemas portanto sobredeterminados.

Este método permite a obtenção da solução cujo somatório do quadrado dos erros é mínimo.

Partindo-se de um modelo linear:

d = Gm, (1.17)

o erro entre o dado medido d e o dado calculado Gm, pode ser escrito como:

e = d−Gm. (1.18)

Calculando o somatório do quadrado dos erros, tem-se:

S = eTe, (1.19)

que pode ser expressa como

S = (d−Gm)T (d−Gm), (1.20)

S = (dT −mTGT )(d−Gm), (1.21)

ou

S = dTd− dTGm−mTGTd + mTGTGm. (1.22)

Derivando S em relação a m e igualando a zero, com o intuito de se encontrar um erro

mínimo, temos:∂S

∂m= 0− dTG−GTd + 2GTGm = 0. (1.23)

Como,

dTG = GTd, (1.24)

então a solução para o problema dos mínimos quadrados é dado por um sistema de equações

denominadas normais, que é representada por:

GTGm = GTd. (1.25)

Como se trata de uma matriz G não singular, multiplicamos os dois lados por (GTG)−1,

logo:

m = (GTG)−1GTd, (1.26)

se a matriz G é não singular. Porem, podem ocorrer casos em que a matriz GTG é quadrada,

mas é singular, não podendo então ser obtida a inversa (GTG)−1. Nestes casos, utilizamos o

método dos mínimos quadrados amortecidos.

9

1.4.2 Método dos Mínimos Quadrados Amortecidos

Uma alternativa para a solução de sistemas subdeterminados é a aplicaç

Para os casos nos quais a matriz GTG é singular, minimiza-se a seguinte função:

S = eTe + ε2mTm. (1.27)

A solução obtida dessa situação tem o quadrado de erro mínimo, embora a funçõ passe

a ter o termo adicional ε2mTm, além do termo eTe utilizado na expressão do método dos

mínimos quadrados convencional. Essa solução terá o erro mínimo se o fator escalar ε2

multiplicado pelo quadrado do modelo, não afete signi�cativamente a função, assim pode-se

adotar a solução obtida.

Minimizando S, obtemos a seguinte expressão:

GTd = (GTG+ ε2I)m. (1.28)

Pré-multiplicando ambos os lados por (GTG+ ε2I)−1 , obtemos:

m = (GTG+ ε2I)−1GTd, (1.29)

uma vez que ε pode ser tomado de tal modo que o termo (GTG + ε2I) da equação seja

não-singular, logo, inversível. Deste modo, pode-se resolver o sistema linear mesmo quando

GTG é singular.

A escolha do fator ε2 deve manter um compromisso entre um valor pequeno que não

afete de forma a comprometer a solução e um valor su�cientemente grande para que o sistema

possa ser resolvido.

1.4.3 Decomposição de Valores Singulares

A decomposição de valores singulares ou SVD, do inglês singular value decomposition, é uma

técnica empregada para a simpli�cação de matriz associada a uma transformação linear. As

matrizes que representam problemas geofísicos são, via de regra, retangulares, não existindo

portanto uma matriz inversa. Nessas condições pode-se determinar a pseudo-inversa de uma

determinada matriz G utilizando-se a SVD.

Seja G uma matriz real M × N . A matriz G+, de dimensões (N × M) será a sua

pseudo-inversa com as seguintes propriedades:

• GG+G = G,

10

• G+GG+ = G+,

• (GG+)T = GG+,

• (G+G)T = G+G.

Se essas propriedades forem satisfeitas, a pseudo-inversa G+ será única. Para calcular

a inversa a partir da matriz original G pode-se utilizar a SVD.

Supondo uma matriz retangular GM×N de posto k, sua SVD é da seguinte forma:

G = UΣV T , (1.30)

tal que: UM×M é a matriz que contém os autovetores ortonormalizados de GGT ; ΣM×N

é a matriz que contém as raízes quadradas dos autovalores de GTG, denominados valores

singulares, sendo estes colocados em ordem decrescente, ou seja, σ1 ≥ σ2 ≥ ... ≥ σk > 0;

VN×N é a matriz que contém os autovetores ortonormalizados de GTG

Então a pseudo-inversa, ou inversa generalizada, a inversa natural, é uma matriz N×M :

G+ = V Σ+UT , (1.31)

onde Σ+, de dimensão N ×M , é a matriz que contém os recíprocos dos valores singulares

não nulos de G da seguinte forma:

Σ+ =

[E 0

0 0

], (1.32)

onde E é a matriz diagonal k× k cujo i-ésimo elemento diagonal é eii = σ−1i para 1≤ i ≤ k:

E =

σ−11 0 · · · 0

0 σ−12 · · · 0...

.... . .

...

0 0 · · · σ−1k

. (1.33)

Em se tratando de matrizes quadradas, de posto completo, a inversa clássica e a pseudo-

inversa são idênticas, sendo a pseudo-inversa uma generalização da inversa clássica, podendo

ser aplicada na inversão de dados em substituiçõ a inversa clássica. A pseudo-inversa G+ em

substituição direta a inversa clássica G−1, diferindo da notação adotada em Menke (1984) que

de�ne G+ = (GTG)−1GT para o caso de mínimos quadrados. Embora esta última notação

seja mais concisa, a substituiçõ explícita da inversa pela pseudo-inversa adotada a seguir

torna mais clara a similaridade entre expressões análogas.

11

Assim, partindo-se das expressões:

m = G−1d, (1.34)

e

m = (GTG)−1GTd, (1.35)

obtém-se as expressões análogas, utilizando-se pseudo-inversa, respectivas:

m = GTd, (1.36)

e

m = (GTG)+GTd. (1.37)

CAPÍTULO 2

Geotomogra�a

A tomogra�a é uma técnica de reconstrução de imagens através do mapeamento de uma

propriedade física, a partir da soma dos valores dessas propriedades em determinadas direções

(projeções). As projeções são as quantidades mensuráveis que são funções de propriedades

físicas de um objeto. Assim a reconstrução tomográ�ca é um tipo especial de problema

inverso que permite estimar uma função.

A tomogra�a clássica requer uma distribuição das fontes e dos receptores em torno do

objeto a ser imageado. Em técnicas de imageamento, tais como a tomogra�a médica de

raios-X, a fonte e o receptor giram em torno do objeto e a maioria das medidas de raios-X

e reconstrução de sistemas usam a energia transmitida (e não re�etida). Assume-se que

estes raios-X se propagam em feixes retos e possuem uma geometria dos raios relativamente

simples. Contrariamente, a tomogra�a sísmica, por utilizar geometrias de aquisição com

fontes e receptores dentro do poço, possui uma abertura angular limitada do objeto em

análise, e as ondas sísmicas, por facilmente se re�etirem, refratarem e difratarem, são muito

menos comportadas em vários sentidos, principalmente por possuirem uma geometria dos

raios mais esparsa e complexa.

2.1 Tomogra�a Sísmica

Em tomogra�a sísmica os dados utilizados são usualmente medidos dos tempos de trânsito,

amplitude da onda ou alguma estimativa da dispersão associada a um caminho que liga uma

fonte a um receptor. A tomogra�a sísmica pode fazer uso das ondas re�etidas, transmitidas,

refratadas e difratadas. A escolha depende da geometria de aquisição, da qualidade dos

dados e dos objetivos a serem atingidos.

Um raio é um �uxo de energia indo através do meio, medido entre a fonte e o receptor.

Para meios homogêneos e isotrópicos o caminho do raio é uma reta. Em geral, este caminho

depende do meio e portanto o raio é curvo.

A geometria de aquisição admite várias con�gurações possíveis. Se as fontes e os recep-

tores estão na superfície, e se considera apenas as ondas re�etidas, chama-se de tomogra�a

12

13

sísmica de re�exão. Se as fontes estão na superfície e os receptores num poço, ou ambos

conjuntos instalados em poços adjacentes, pode-se estudar o tempo de trânsito da onda

transmitida, e neste caso chama-se tomogra�a sísmica de transmissão. Para utilizar mais

informações e amenizar a presença de heterogeneidade, pode-se usar as ondas difratadas,

este tipo é chamado de tomogra�a sísmica de difração.

2.1.1 Tomogra�a de Tempos de Trânsito

A tomogra�a de tempos de trânsito, ou tomogra�a de raio, utiliza o tempo de trânsito en-

tre as fontes e os receptores como vetor, d, de dados observados na inversão; na forma da

onda são desconsiderados os parâmetros amplitude×tempo. A matriz G, utilizada n inver-

são da tomogra�a de tempos de trânsito, descreve a geometria dos raios em um problema

tomográ�co. Os parâmetros do modelo a serem estimados compõem o vetor m, e correspon-

dem a vagarosidade (inverso da velocidade). O tempo de trânsito é a integral de linha da

vagarosidade ao longo do raio:

t =

∫r

s(x, z)dl, (2.1)

onde t é o tempo de trânsito; r é o raio ao longo do qual é realizada a integração; dl é o

elemento de raio e s(x, z) é a vagarosidade do meio no ponto(x, z), onde x é a coordenada

horizontal e z é a coordenada vertical.

Como o caminho do raio depende da vagarosidade, a equação do tempo de trânsito é

não-linear. Para simpli�car a resolução do problema, reduz-se o mesmo a um modelo linear.

Podemos escrever o tempo t em função de uma variável g:

tr = g[s(x, z)], (2.2)

onde:

g[s(x, z)] =

∫r

s(x, z)dl. (2.3)

Os tempos de trânsito t, de M raios, podem ser representados pelo vetor t. Por seu

turno a distribuição de vagarosidade s(x, z) pode ser representada pelo vetor s.

Expandindo-a em série de Taylor em torno de uma vagarosidade de referência s0:

t = t0 +∂g

∂s

∣∣∣∣s=s

0

(s− s0) + · · · . (2.4)

Truncando a série na primeira derivada, obtém-se:

14

t− t0 =∂g

∂s

∣∣∣∣s=s

0

(s− s0), (2.5)

ou seja,

∆t = G∆s, (2.6)

onde o vetor ∆t corresponde as diferenças entre os tempos de trânsito calculados e os tempos

de trânsito observados para o modelo; o vetor ∆s corresponde as diferenças entre as vagaro-

sidades do modelo inicial e as vagarosidades estimadas; a matriz G contém os elementos gijque correspondem as distâncias que o i-ésimo raio percorre o j-ésimo bloco.

As equações que descrevem a tomogra�a de tempos de trânsito, aproximada na forma

de uma equação matricial, podem então ser formuladas como:

t1

t2...

tM

=

g11 g12 . . . g1N

g21 g22 . . . g2N...

.... . .

...

gM1 gM2 . . . gMN

s1

s2...

sN

. (2.7)

Podendo ser reescrita na forma de equação de um modelo linear como:

t = Gs, (2.8)

ou, utilizando a notação que vem sendo utilizada neste presente trabalho:

d = Gm, (2.9)

sendo o vetor m equivalente a vagarosidade, que é o recíproco da velocidade, isto é:

mi = si =1

vi, (2.10)

onde t é o vetor de tempos de trânsito, que correspondem ao vetor de dados observados

d; G é a matriz núcleo; s é o vetor vagarosidade (inverso da velocidade), que representa os

parâmetros do modelos a serem determinados na tomogra�a de tempos de trânsito.

O caminho a ser percorrido pelo raio que sai da fonte e vai em direção ao receptor é

dado pelos elementos da matriz G:

G =

g11 g12 . . . g1N

g21 g22 . . . g2N...

.... . .

...

gM1 gM2 . . . gMN

, (2.11)

onde M é o número total de raios, ou seja, o produto do número de fontes pelo número de

receptores; N representa o número total de blocos na qual a área de estudo foi dividida.

15

Portanto, cada linha de G corresponde ao caminho de um raio. A matriz G é esparsa

pois um dado raio intercepta apenas uma pequena parte do modelo.

2.1.2 Tomogra�a de Atenuação

Historicamente, as grandes descobertas de acumulações de hidrocarbonetos no Brasil, em

especial na Bacia de Campos, foram realizadas com base na análise de mapas de amplitudes

sísmica. Esses mapas reproduzem contrastes de impedância elástica entre camadas, como

também o nível de atenuação que as ondas elásticas sofrem os se propagar nas camadas ro-

chosas. Outra técnica de análise do dado sísmico, a técnica da amplitude versus afastamento

ou AVO (amplitude versus o�set), também é baseada no estudo da energia do sinal sísmico.

No Brasil, pouco se estudou a respeito dos coe�cientes de atenuação das rochas que ocorrem

nas nossas bacias sedimentares.

As rochas são materiais perfeitamente elásticos. À medida que se propagam pelas

rochas, as ondas têm suas amplitudes reduzidas pelo processo de atenuação e suas formas

alteradas, sendo suas velocidades dependentes da frequência devido ao processo de dispersão.

Estes fenômenos são interligados e seus efeitos combinados são chamados de absorção sísmica,

atribuída ao comportamento parcialmente elástico das rochas durante a propagação das

ondas sísmicas. Um efeito da absorção no domínio do tempo é a mudança da forma no sinal,

enquanto no domínio da frequência é a deformação dos espectros de amplitude e de fase.

Informações como litologia, estado físico, grau e tipo de saturação das rochas podem

ser estimadas a partir de suas propriedades de atenuação, tornando a atenuação importante

para muitas aplicações na indústria do petróleo.

Estudos petrofísicos nessa área ainda são escassos. Entretanto, Moreira (2009) em seu

trabalho fez a madição de diferentes litogias em testemunhos de sondagens. Com isso, ob-

teve o coe�ciente de atenuação para diferentes litogias. Quan et. al (1997) �zeram uma

abordagem multifrequencia sobre a tomogra�a de tempos de trânsito e tomogra�a de atenu-

ação poço a poço. Plessix (2006) também fez um estudo onde foi analisado a tomogra�a de

tempos de trânsito e a tomogra�a de atenuação.

Na tomogra�a de atenuação, a propriedade física a ser estudada é a atenuação. Comu-

mente, mede-se a atenuação pelo coe�ciente de atenuação α, que é a constante de decaimento

exponencial da amplitude de uma onda plana em um meio homogêneo, ou pelo fator de qua-

lidade Q (ou seu inverso 1/Q, o fator de dissipação), ou ainda pelo decremento logarítmico

δ. Tais grandezes relacionam-se segundo:

1

Q=αv

πf=δ

π, (2.12)

onde v é a velocidade (m/s) e f é a frequência (Hz). O coe�ciente de atenuação tem como

16

unidade de medida o inverso da distância (m−1) e o fator de qualidade é uma grandeza

adimensional.

Podemos escrever a amplitude de uma onda plana, em função do espaço e do tempo,

propagando-se em um meio elástico homogêneo, de coe�ciente de atenuação α, no caso

unidimensional, como:

A(x, t) = A(x)ei(kx−ωt), (2.13)

onde k é o número de onda (k = 2π/λ, onde λ é o comprimento de onda), ω é a frequência

angular (ω = 2πf , onde f é a frequência linear).

Nos casos em que a atenuação da onda não é considerada, temos que A(x) = A0. Para

os demais casos, tem-se que:

A(x) = A0e−αx, (2.14)

onde A0 é a amplitude inicial da onda e α é o coe�ciente de atenuação do meio.

Uma forma de obtermos o coe�ciente de atenuação é através de dois registros de am-

plitudes A(x1) e A(x2), em dois pontos consecutivos, x1 e x2 (sendo que x1 < x2). Com isso,

temos:

α =ln[A(x1)A(x2)

]x1 − x2

. (2.15)

A unidade de α é nepers por unidade de comprimento ou, apenas, o inverso do comprimento

(m−1).

Na tomogra�a de atenuação, a propriedade de atenuação da Terra é relatada através

da amplitude (ak) da onda. A tomogra�a é a inversão desta integral de linha obtendo a

atenuação:

ln

(A0

Ak

)=

∫Rk

αdl, (2.16)

onde dl é o elemento do caminho, Rk é o caminho do raio e k=1,...,M . Que seria uma

equação na qual podemos notar a relação entre ela e a equação (2.15).

2.2 Modelagem Tomográ�ca Direta

Uma maneira simples de parametrização de um meio seria dividir a região de estudo em

pequenas células, denominados pixels no caso bidimensional, e atribuir valores constantes a

17

propriedade física em estudo em cada célula.

A propagação de sinais sísmicos em ummeio pode ser representada de diversas maneiras,

tais como modelos analíticos, diferenças �nitas e traçcados de raios. A modelagem analí-

tica é bastante e�ciente computacionalmente, mas restringe-se a casos relativamente simples,

tornando-a limitada. Através das diferenças �nitas podemos modelar meios complexos, en-

tretanto o seu custo computacional é mais elevado. A modelagem utilizando traçados de

raios permite modelar meios com exatidão e e�ciência computacional.

Utilizou-se neste trabalho a modelagem 2-D e a teoria do traçado de raios como mode-

lagem da propagação da energia através do meio estudado.

2.2.1 Traçado de Raios

Um raio é o caminho percorrido pela energia que vai de uma fonte a um receptor. Para meios

homogêneos, o ângulo de transmissão do raio, ao atravessar as interfaces dos meio, pode se

considerado constante, concedendo assim uma aproximação por raios retos. Entretanto, o

meio geológico que é encontrado, raramente é simples, logo a idéia dos raios retos não se

aplica, já que a geometria dos raios é mais complexa. Então obtêm-se o comprimeto do arco

do raio percorrido percorrido em cada célula do meio, dando uma aproximação de propagação

por raios curvos. O traçado de raios, considerando a curvatura do raio, geralmente é utilizada

para métodos iterativos, onde a matriz é atualizada a cada iteração.

O traçado de raios curvos pode ser analisado pelo princípio de Fermat, que diz que

a energia se propaga ao longo de caminhos que torna o tempo de trânsito mínimo, sendo

este caminho denomidado como raio. Matematicamente, o tempo de trânsito de um raio,

viajando de um ponto p1 até um ponto p2, é dado por:

t =

∫ p2

p1

n(x, z)

cdl, (2.17)

onde t é o tempo de trânsito, n(x, z) = c/v(x, z) é a distribuição dos índices de refração

no meio bidimensional, sendo c a velocidade do som em um meio de referência e v(x, z) a

velocidade do som no ponto (x, z) e dl o diferencial do comprimento do arco ao longo do

raio.

Pelo princípio de Fermat, o caminho será aquele para o qual a integral acima, da equação

(2.17), assume um valor estacionário, logo:

I =

∫ p2

p1

n(x, z)dl, (2.18)

onde I é o caminho do raio; P1 e P2 são os dois pontos conectados por um determinado

raio, sendo usualmente esses pontos as posições respectivas da fonte e do receptor; n(x, z) =

18

c/v(x, z) é o índice de refração bi-dimensional correspondente à posição (x, z) do meio para

uma dada velocidade de referência c e uma determinada velocidade real v(x, z) obtida a

partir do campo de velocidades; e dl é o diferencial do comprimento de arco ao longo do raio.

O campo de índices de refração, obtido a partir do campo de velocidades do meio, é

discretizado em uma malha, a partir do qual interpola-se o índice de refração para obter

uma estimativa do mesmo em um ponto atravessado pelo raio, como mostra a Figura (2.2.1)

abaixo.

Figura 2.1: Campo discretizado e interpolação bilinear dos índices de refração, onde

cada cruzamento da retícula corresponde a um índice de refração das

coordenadas respectivas e P (∆x,∆z) representa o ponto interpolado

(Terra, 2007)

Para um dado ponto nas coordenadas (∆x,∆z), onde ∆x e ∆z são normalizados entre

0 e 1 na região de interesse, o indíce de refração interpolado é dado por:

P (∆x,∆z) = A(1−∆x)(1−∆z) +B∆x(1−∆z) + c(1−∆x)∆z +D∆x∆z. (2.19)

Deste modo, obtém-se uma estimativa doíndice de refração no ponto atual utilizando-se

os índices de refração conhecidos nos pontos A, B, C e D das células vizinhas, permitindo

19

calcular a trajetória do raio entre a fonte e receptor passando por diversos pontos em uma

mesma célula do campo de velocidades discretizado.

A equação de Euler é uma condição necessária para a existência de um valor extremo

da integral variacional para o cálculo do comprimento acústico, atendendo ao princípio de

Fermat. Utilizando-se esta equação pode-se obter uma equação diferencial para uma família

de raios em um meio homogêneo:

d

dl

(ndrdl

)= ∇n, (2.20)

onde n(x, z) é o campo de índices de refração; r é o vetor posição do raio; dr/dl é o vetor

unitário tangente ao raio em (x, z); dl é o elemento de comprimento da trajetória de raio;

∇n é o gradiente do campo de índices de refração.

Esta equação diferencial é denominada como equação do raio. Para uma certa vizi-

nhança regular, ou seja,em que o índice de refração varia suavemente, sua solução representa

para os raios de menor comprimento acústico. A especi�cação dos raios desejados é obtida

através da imposição de condições de contorno (Andersen e Kak, 1982), tais como a posição

da fonte e do receptor bem como o ângulo de partida do raio.

Desenvolvendo a equação do raio, obtém-se:

dn

dl

drdl

+ nd2rdl2

= ∇n. (2.21)

Fazendo a seguinte substituição:

dn

dl=dn

dr· drdl

= ∇n · drdl, (2.22)

então: (∇n · dr

dl

)drdl

+ nd2rdl2

= ∇n. (2.23)

Expandindo o vetor posição em série de Taylor no ponto (l+∆l), onde ∆l é o incremento

do raio, e truncando a série na segunda derivada, obtém-se:

r(l + ∆l) ≈ r(l) +drdl

∆l +1

2

d2rdl2

∆l2. (2.24)

Isolando o vetor curvatura, d2r/dl2, da equação (2.23) e substituindo na equação (2.24),

obtém-se a seguinte aproximação:

r(l + ∆l) ≈ r(l) +drdl

∆l +1

2n

[∇n−

(∇n.dr

dl

)drdl

]∆l2. (2.25)

20

Esta expressão permite estimar a nova posição do raio a partir do conhecimento de

sua posição atual r(l), do vetor tangente dr/dl e do gradiente da distribução dos índices de

refração, ∇n, naquela posição, constituindo assim numa expressão para a implementação

dos algoritmos do traçado de raio. A partir desta expressão, pode ser derivada uma outra

expressão utilizando o ângulo instantâneo (α) entre a tangente do raio e o eixo x, como pode

ser visualizado na Figura (2.2).

Figura 2.2: Interpretação aproximada da in�uência da curvatura, d2r/dl2, do ângulo

instantâneo, α, entre o vetor tangente, dr/dl, e o eixo horizontal, x, e

do novo passo, r(l + ∆l).

Considerando o vetor unitário na direção de propagação do raio para a k-ésima iteração

como:

drdl

= (cosαk) i+ (senαk) j, (2.26)

onde αk é o ângulo entre a tangente do raio e o eixo x para a iteração k do algoritmo de

traçado de raios.

Para um dado ângulo instantâneo αk do raio, se obtém a seguinte aproximação entre

as iterações k − 1 e k:

21

senαk =zk − zk−1

∆l, (2.27)

e

cosαk =xk − xk−1

∆l. (2.28)

O gradiente do índice de refração é de�nido por:

∇n =

(∂n

∂x

)i +

(∂n

∂z

)j. (2.29)

O produto interno da equação (2.25) e que será denominado de d, pode ser reescrito

em função das equações (2.26) e (2.29), como:

d = ∇n.drdl

= nxcosα + nzsenα (2.30)

onde:

nx =∂n

∂x, (2.31)

e

nz =∂n

∂z. (2.32)

As derivadas direcionais do meio discretizado são aproximadas utilizando-se diferenças

�nitas:

nx(i, j) =n(i+ 1, j)− n(i− 1, j)

2∆x, (2.33)

e

nz(i, j) =n(i, j + 1)− n(i, j − 1)

2∆z. (2.34)

O ponto seguinte do raio é estimado pelas equações:

xk+1 = xk + cosαk∆l +1

2nk(nk,x − dkcosαk)∆l2, (2.35)

e

zk+1 = zk + senαk∆l +1

2nk(nk,z − dksenαk)∆l2, (2.36)

onde para um dado ponto (xk, zk) pertencente ao raio, na k-ésima iteração: nk é o índice de

refração; nk,x é a derivada do índice de refração na direção x; nk,z é a derivada do índice de

refração na direção z; ∆l é constante e de�ne o deslocamento na trajetória do raio em cada

iteração do algoritmo; e dk é a derivada direcional.

Considera-se a equação:

22

nk =c

vk= csk, (2.37)

onde sk é a vagarosidade correspondente à velocidade vk. Substituindo esta equação nas

equações (2.35) e (2.36), obtém-se:

xk+1 = xk + cosαk∆l +1

2sk(sk,x − dkcosαk)∆l2, (2.38)

e

zk+1 = zk + senαk∆l +1

2sk(sk,z − dksenαk)∆l2. (2.39)

Onde dk foi de�nido como:

dk = sk,xcosαk + sk,zsenαk. (2.40)

Portanto, de forma iterativa, partindo-se de um determinado ponto inicial (x0, z0) cor-

respondente à posição da fonte e de um ângulo inicial (α0), obtém-se sucessivamente os

pontos seguintes do raio a partir das expressões que estimam o ponto seguinte do raio. Vale

salientar que o ângulo αk não precisa ser calculado explicitamente, uma vez que apenas os

valores do seno e do cosseno do ângulo αk são utilizados, podendo ser calculado diretamente,

conforme as equações (2.27) e (2.28).

Em muitas aplicações práticas, a discretização em malha da distribuição de índices de

refração e os erros relacionados ao cálculo das derivadas parciais ∂n/∂x e ∂n/∂z impõem

algumas limitações para este método iterativo de traçado de raios (Andersen e Kak, 1982).

Então, como a posição do ponto seguinte é calculada a partir do ponto atual e do índice

de refração corrente e suas respectivas derivadas parciais, erros causados pela discretização

ou transições bruscas de velocidade podem se tornar cumulativos. Para minimizar este

problema, deve-se adotar uma malha com resolução su�ciente para que o meio seja amostrado

adequadamente, resultando em transições mais suaves de velocidade. Também podem ser

utilizadas suavização do campo de velocidades e interpolação bilinear dos índices de refração

e de suas derivadas parciais respectivas, segundo Andersen e Kak (1982), resultando em uma

exatidão maior, pois o raio pode percorrer diversos pontos diferentes em uma mesma célula.

Depois de �xadas as posições das fontes e dos receptores, resta ainda predizer o ângulo

inicial do raio para que uma determinada fonte seja atingida pelo mesmo.

2.2.2 Discretização do Meio

Para uma abordagem tomográ�ca necessitamos discretizar o meio, ou seja, dividi-lo para

que este possa ser tratado computacionalmente. No caso 2D, o procedimento consiste em

dividir o meio em blocos quadrados de dimensões δx e δz, como mostrado na Figura (2.3).

23

A partir desta divisão, podemos determinar com a equação do raio, a distância que

cada raio percorre em cada bloco e, desta forma, montar a matriz tomográ�ca, ou seja,

dividiremos o meio em células.

Figura 2.3: Exemplo de um modelo discretizado

Apesar, da frente de onda percorrer todo o meio, os raios não o farão, visto que este

é apenas um ente matemático que indica o caminho que a onda percorre. Logo, torna-se

necessário uma boa disposição de fontes e receptores e uma boa discretização do meio, a �m

de que não ocorra uma má iluminação do meio.

2.2.3 Ligação entre fonte e receptor

O maior problema do traçado de raios é o de encontrar as coordenadas em que o receptor se

encontra através do raio que partiu da fonte. Considerando um meio homogêneo, a ligação

entre fonte e receptor é tida como um raio reto, tornando em um problema rápido e fácil

de ser resolvido. Porém, considerando um meio heterogêneo deve-se considerar a curvatura

agregada ao raio. Nesta situação, o problema condiz em predizer o ângulo de lançamento do

raio que parte da fonte para que ele alcance um determinado receptor.

O processo de ligação entre a fonte e o receptor é denominado de ray linking. Esse

problema pode ser resolvidos por alguns métodos, como exemplo temos:

24

Método de canhoneio (shooting method)

Dados um ângulo inicial e outro ângulo �nal, sendo um valor máximo (θmax) e o outro o

valor mínimo (θmin), traçam-se os raios compreendidos entre essa região, limitada por esses

dois extremos, com um incremento angular ∆θ, tal que o ângulo de emissão não saia desse

região delimitada. O ângulo de lançamento a ser adotado, será aquele para o qual o raio

esteja mais próximo do receptor desejado. Sendo:

θi = θmin + i∆θ, para i = 0, 1, 2, ..., n, (2.41)

sendo n um número inteiro de valor máximo em que θi = θmax. A Figura (2.4) faz uma

representaçãode como o método funciona.

Figura 2.4: Esquema do método do canhoneio.

Percebe-se que o custo computacional desse método é elevado a depender do incremento

∆θ tomado.

Método de encurvamento (bending method)

Neste método são �xados os dois extremos do raio (a fonte e o receptor), há uma juste da

curvatura do raio até minimizar o tempo de trânsito total ao longo do trajeto.

25

Método de ligação (linking method)

Esse método de conexão é baseado no traçado de três raios a partir de um raio inicial, com

um ângulo inicial de lançamento θ2 e um incremento angular ∆θ, tal que os dois rais traçados

posterior a esse primeiro terão os ângulos de lançamentos dados por:

θ1 = θ2 −∆θ (2.42)

e

θ3 = θ1 + ∆θ, (2.43)

sendo o ângulo de partida (θ2) formado entre a linha reta que liga a fonte e o receptor e a

horizontal.

Depois de traçados esses três raios, podemos aplicar o re�namento de Newton-Raphson,

e então é encontrada uma estimativa para o próximo ângulo de lançamento (θ′2), do raio 2,

que estará mais próximo ao receptor, sendo dado por:

θ′

2 = θ2 ±2d2rd13

∆θ, (2.44)

onde d2r é a distância entre as coordenadas �nais do raio 2 com o ângulo inicial θ2 e as

coordenadas do receptor; d13 é a distância entre as coordenadas �nais do raio 1 e 3; o sinal ±depende da direçãodo crescimento angular e das coordenadas do raio 2 e do receptor (o sinal

escolhido para a equação (2.44) será positivo caso encaminhe o novo ângulo, θ′2, no sentido

do receptor e negativo caso contrário), como mostrado na Figura (2.5).

2.2.4 Modelagem Direta da Tomogra�a de Atenuação

A modelagem direta pode ser descrita, a partir da equação (2.16), matematicamente, na

forma matricial, como:

dobs = Gαver, (2.45)

onde dobs é o vetor de dados observados na qual tem os valorer do logaritmo neperiano da

razão entre as amplitudes inicial e �nal que chega no receptor; G é a matriz tomográ�ca das

distâncias, a qual foi calculada pelo traçado de raio, e αver ó vetor que representa o modelo

verdadeiro do campo de atenuações do meio geológico em estudo.

Uma outra forma de representar a modelagem direta é partindo da equação (2.14),

tendo a equação na seguinte forma:

Aj = A0e−

∑Ni=1αigji , (2.46)

26

Figura 2.5: Esquema do método do ligação.

onde Aj é a amplitude no j-ésimo receptor; A0 é a amplitude inicial; αi é o coe�ciente de

atenuação no i-ésimo bloco e gji é o elemento da matriz G distância percorrida entre um

ponto e outro.

Que pode ser reescrita na seguinte forma:

AjA0

= e−∑N

i=1αigji . (2.47)

Sendo simpli�cada por:

ln

(A0

Aj

)=

N∑i=1

αigji. (2.48)

Na modelagem direta os dados de entrada do modelo serão o vetor (αver) que contém

os valores do coe�ciente de atenuação em cada bloco e a matriz G das distâncias que os

raios percorrem entre fonte-receptor; o dado observado (dobs) será o logaritmo neperiano da

razão entre o amplitude inicial (A0), que é aquela lida na fonte, e a amplitude registrada na

j-ésima fonte.

Sendo que o vetor de dados será composto pelos valores de cada dado lido, estes que

são representados matematicamente como:

dj = ln

(A0

Aj

). (2.49)

27

2.3 Modelagem Tomográ�ca Inversa

A tomogra�a inversa ou, simplesmente, inversão consiste em estimar, a partir dos dados me-

didos (dobs), no caso desse trabalho a amplitude inicial e a amplitude registrada no receptor,

o valor do campo de atenuações do meio bi-dimensional. Esse modelo 2-D é rasterizado em

um vetor 1-D que consiste em varrer linha por linha os elementos do modelo 2-D e armazená-

los sequencialmente no vetor α, ou seja, o último elemento de cada linha do modelo 2-D

passa a ser adjacente ao primeiro elemento do linha seguinte.

A solução de αver é obtida através da inversão de matrizes utilizando-se a técnica

já descrita anteriormente no item 1.4.3, a chamada SVD. A equação (1.37) descreve essa

inversão, que no caso da atenuação descrevemos como:

αest = (GTG)+GTdobs, (2.50)

onde αest será o vetor estimado que conterá os valores das atenuações; G é a matriz que

corresponde ao caminho percorrido por um raio; dobs é o vetor de dados observados, como

já visto anteriormente; e (GTG)+ é a pseudo-inversa da matriz GTG.

CAPÍTULO 3

Simulações e Resultados

Neste cap�'tulo serão apresentados e analisados as simulações obtidas com a tomogra�a

de atenuação. Os estudos foram feitos partindo de duas simulações. No primeiro modelo,

para testar e validar o método, foi feito em um modelo simples, onde existe um meio com

atenuação praticamente nula e apenas uma camada com atenuação considerável. Um segundo

modelo foi feito baseado em um dos modelos de Oliveira (2010), onde o modelo é faz uma

simulação de um meio geológico que representa um anticlinal.

Inicialmente o modelo verdadeiro é criado com informações a priori que são os valores

dos coe�cientes de atenuação do meio e a sua respectiva posição. O vetor de dados calcu-

lados (dobs) são obtidos a partir da equação (2.49). Partindo da dimensão do modelo que é

trabalhado, é montada a matriz G, que é representada pela equação (2.11). Para ser feita

a modelagem inversa é necessário obtermos a inversão da matriz G, que neste trabalho foi

realizada pela decomposição de valores singulares, mostrada no item (1.4.3) deste trabalho.

Com isso, pela equação (2.50) é obtido o modelo estimado de atenuação do meio.

A forma de resolução do problema foi esquematizada na �gura (3.1) em que é apresen-

tado esquematicamente como foi realizado esse trabalho.

Inserção de Ruído

Para testarmos a solução de uma maneira um pouco mais real, inserimos um ruído gaussiano

no vetor de dados. Esse vetor de ruído aleatório pode ser representado como:

d∗j = dj + βrjdj, (3.1)

onde d∗j é o vetor de dados contaminado com ruído; dj é o vetor de dados sem ruído; β é o

fator de amplitude do ruído e rj é a sequência de números aleatórios que variam se encontram

no conjunto de valores de (-0,5 ; 0,5).

A exemplo, a �gura (3.2) mostra o vetor de dados sem in�uência de ruído (a), e sob a

in�uência de ruído (b).

28

29

Inicializa o vetor atenuação 𝛂ver = [α1, α2, … , αN]

T

a partir de informações a priori.

Calcula a matriz G, utilizando traçado de raios.

Calcula o vetor dado observado:

dj = ln(A0

Aj),

e coloca em forma de vetor:

𝐝𝐜𝐚𝐥𝐜 = [d1, d2, … , dM]T.

Realiza a inversão da matriz G pela decomposição em valores singulares, achando a matriz pseudo-

inversa da matriz G, que é a matriz 𝐆+.

Seleciona a quantidade de valores singulares convenientes.

Calcula o modelo de atenuação estimado: 𝛂est = 𝐆+𝐝𝐜𝐚𝐥𝐜

A imagem do tomograma é satisfatória?

S

N

Calcula o erro entre 𝛂ver e 𝛂est

O erro é satisfatório?

N

S

FIM

Figura 3.1: Fluxograma

Os critérios de seleção dos valores singulares são discutidas na seção seguinte desse

capítulo.

30

(a) (b)

Figura 3.2: Vetor de dados, em (a) sem in�uência de ruído e em (b) sob in�uência

de ruído da ordem de β = 10−1

3.0.1 Critérios de Seleção de Valores Singulares

O critério de seleção de corte utilizado nesse trabalho foram baseados em alguns dos critérios

utilizados por Silva (2009). Serão utilizados apenas o critério de amplitude do valor singular

e a derivada do valor singular.

Para ser feita uma comparação e segurança na escolha do ponto de corte ótimo, utili-

zamos também o erro RMS entre os modelos verdadeiro e estimado. Entretanto, este não

pode ser considerado como um critério, já que, em um caso real, o modelo verdadeiro não é

conhecido, logo não é possível determinar tal função.

No capítulo 1 deste trabalho houve uma abordagem dos problemas inversos mal-postos

e foi visto que a estabilidade do problema diminui com o aumento de uma grandeza escalar

conhecida como número de condição (NC), esta é de�nida como a razão entre o maior e o

menor valor singular utilizados na composição da matriz referente ao sistema linear. Logo,

valores singulares muito pequenos comprometerão a condição de estabilidade do problema.

Então, devemos eliminá-los para se obter um sistema estável. Contudo, se tirarmos todos

os valores singulares pequenos, a ponto de encontrar um número de condição muito próximo

da unidade, com a �nalidade de aumentar a estabilidade do problema, não conseguiremos

encotrar uma matriz pseudo-inversa muito boa, pois ao diminuirmos a quantidade de valores

singulares, diminuimos, também, a quantidade de autoimagens quer irão compor a matriz

inversa generalizada e, desta forma, reduzimos a informação contida no sistema. Então

torna-se necessário encontrar uma quantidade ótima de valores singulares que nos forneça

uma boa pseudo-inversa e que também amenize o mal-condicionamento do problema, ou

seja, encontrar um ponto intermediário que forneça a melhor solução factível.

31

Amplitude do valor singular

Nesse método, o interesse é gerar uma curva davariação da amplitude do valor singular em

função do índice destes valores, com o propósito de se analisar o decaimento dos valores

singulares e encontrar um ponto de corte ótimo que, quando aplicado à inversão, possa

fornecer uma solução satisfatória.

Nesta curva, a variação da amplitude dos valores singulares se dá em uma escala mo-

nologarítmica, para que desta forma possamos encontrar uma região anômala que poderá

sugerir o ponto de corte ótimo. Esta regi¿o pode ser determinada por, relativamente, uma

alta variação na amplitude dos valores singulares, sugerindo assim uma grande mudança no

comportamento destes.

Derivada do valor singular

Este é um critério que consiste em determinar a derivada de todos os valores singulares e,

posteriormente, se gerar uma curva desta em função do índice do valor singular. A derivada

para cada ponto i é:

f′(xi) = lim

xi→xi+1

f(xi+1) − f(xi)xi+1 − xi

, (3.2)

onde f(xi) representa a amplitude do valor singular e xi o índice do valor singular, sendo,

portanto, valores inteiros. Além disso, a variação do índice do valor singular é da ordem da

unidade, logo, a equação anterior pode ser reescrita como:

f′(xi) = f(xi+1) − f(xi). (3.3)

Contudo, a amplitude dos valores singulares decresce com a quantidade de valores.

Logo, podemos perceber que os valores da derivada serão negativos. No trabalho foi utilizado

o módulo dos valores da derivada do valor singular, como o propósito de se trabalhar com

valores positivos.

O objetivo de se utilizar a derivada do valor singular como critério de seleção é o

mesmo da amplitude, de se encontrar um ponto onde a curva tem uma mudança brusca

de direção indicando que os valores singulares alteram seu comportamento. Logo, o ponto

dein�exão poderá mostrar o número z'otimo, pois indica que os valores singulares aumentam

consideravelmete a instabilidade.

É importante ressaltar que a utilização dos critérios da amplitude e derivada do valor

singular, são in�uenciados apenas pelo valor singular e ao com a quantidade de valores

utilizadas na inversão, ou seja, um aumento no erro entre os parâmetros de dados pode

modi�car a quantidade de valores singulares que deve ser utilizada para se obter um modelo

32

compatível e variações no erro não podem ser vista nestes critérios, tornando seu uso um

pouco mais limitado. Essa mudança na quantidade de valores singulares a ser utilizada pode

ser explicada a partir da relação entre os erros dos parâmetros de dados e de modelo e o

valor do número de condição, conforme pode ser visto na equação (1.16).

Erro RMS dos parâmetros de dados

A partir dos dados calculados com o modelo estimado e dos dados observados, podemos

determinar o erro RMS relativo, escrita em porcentagem, entre os dados, para a tomogra�a

de atenuação escrita como:

Ed =

√∑Mi=1

(dobsi − dcalci

)2√∑Mi=1

(dobsi)2 × 100 (3.4)

ondeM é o número de raios que percorre o meio, dobsi e dcalci são os vetores de dados observados

e calculados, respectivamente, que o raio i leva para percorrer a distância fonte-receptor.

Erro RMS dos parâmetros do modelo

A partir dos resultados da inversão tomográ�ca, é necessário ponderar o erro entre os parâ-

metros de modelo verdadeiros e os parâmetros do modelo estimado (no caso de um modelo

sintético). Contudo, como dito anteriormente, se este fosse considerado como um critério

de seleção não teria grande signi�cado, pois o modelo verdadeiro não é conhecido em um

caso real. Portanto, esta função é utilizada apenas para �ns de validação das metodologias

empregadas.

O erro RMS relativo, escrito em porcentagem, entre os modelos, pode ser escrita como:

Em =

√∑Ni=1 (mver

i −mesti )2√∑N

i=1 (mveri )2

× 100, (3.5)

onde N é o número de blocos da malha, mveri e mest

i são as propriedades físicas, que neste

trabalho representa o coe�ciente de atenuacdao (α), verdadeira e estimada, respectivamente,

do bloco i.

33

3.1 Modelo 1

O primeiro modelo trata-se de um modelo relativamente simples, onde o meio tem atenuação

muito baixa, porém ele contém uma camada com uma atenuação de 0.25m−1 que pode ser

um possível reservatório alvo. Esse modelo é utilizado para testar a validade do método de

inversão. O meio do modelo, representado na �gura (3.3) foi discretizado em uma malha 20

× 20, ou seja, ele foi dividido em 400 células. A geometria de aquisição dos dados é do tipo

poço a poço, sendo colocado 20 fontes em um poç e 20 receptores em outro, totalizando um

total de 400 raios descritos. Os elementos da matriz G representam a distância em que o

j-ésimo raio percorre o i-ésimo bloco, sendo obtida pela técnica de traçado de raios (descrita

no item 2.2.1), sendo de ordem 400 × 400, consequentemente a sua matriz pseudo-inversa

G+ também será da mesma ordem.

Figura 3.3: 'Modelo 1' verdadeiro, sendo a escala em cores os valores do coe�ciente

de atenuação em m−1.

Depois de realizado o processo da inversão, o grá�co da amplitude do valor singular

e da derivada do valor singular foram obtidas com o intuito de ter uma noção dos valores

singulares e, a partir disso, fazer uma seleção dos valores singulares que entrarão na solução

do problema.

Para analisar a qualidade da inversão, foi feita a inversão do modelo verdadeiro com

diferentes valores de corte de valores singulares (σcorte), que variaram de 10−13 até 105 com

um incremento de 101, em diferentes níveis de ruído variando a sua amplitude (β). Com isso

foi feita a tabela (3.1) que contém os valores de σcorte, de amplitude β do ruído e o erro RMS

dos parâmetros de dados e o erro dos parâmetros do modelo.

34

(a)

(b)

Figura 3.4: Grá�co (a) da amplitude e (b) da derivada do valor singular do 'Modelo

1', em escala monolog.

Foram recuperados alguns modelos em diferentes níveis de ruído (β), variando o valor

do corte de valores singulares (σcorte).

35

3.1.1 β = 0

Os modelos recuperados, sem a adição de ruído no vetor de dados, foram:

(a) (b)

Figura 3.5: 'Modelo 1' recuperado sem ruído onde o valor do σcorte = 103. A escala

em cores representa o valor do coe�ciente de atenuação (α) em m−1.

Em (a) a faixa dos α varia de acordo com o modelo recuperado e em

(b) varia em comparação com o modelo verdadeiro.

(a) (b)





36

(a) (b)

Figura 3.7: 'Modelo 1' recuperado sem ruído onde o valor do σcorte = 10−3. A escala




(a) (b)

Figura 3.8: 'Modelo 1' recuperado sem ruído onde o valor do σcorte = 10−6. A escala




37

(a) (b)

Figura 3.9: 'Modelo 1' recuperado sem ruído onde o valor do σcorte = 10−10. A

escala em cores representa o valor do coe�ciente de atenuação (α) em

m−1. Em (a) a faixa dos α varia de acordo com o modelo recuperado e

em (b) varia em comparação com o modelo verdadeiro.

Como bem visto nas �guras (3.5), (3.6), (3.7), (3.8) e (3.9), e con�rmado pela tabela

(3.1), para os modelos invertidos sem ruído, podemos perceber que os modelos recuperados

foram satisfatórios visualmente, mesmo com os erros dos parâmetros do modelo nos valores

de 43,70%, 34,29%, 28,04%, 28,03% e 18,33%, respectivamente. O valor de corte que teve

o erro mínimo foi o σcorte = 10−10, onde entraram na solução 357 valores singulares, como

comprovação disso, a imagem (3.9), de todos os modelos recuperados, visualmente é o que

mais se aproxima do modelo verdadeiro, �gura (3.3).

A partir do valor de σcorte = 10−12, onde foram incluídos 372 valores singulares, a

recuperação do modelo já não foi mais satisfatórias e o erro dos parâmetros do modelo au-

mentaram signi�cativamente, e isso pode ser notado na �gura (3.4(a)) onde a amplitude do

valor singular após esse valor atrapalharam a solução. A solução ideal para esse método

proposto, admitiria 355 valores singulares na solução. Entretanto, conseguiria recuperar mo-

delos de modo satisfatório, através da análise dos erros dos parâmetros do modelo, qualquer

modelo que admitisse entre σcorte = 104 e σcorte = 10−11, ou seja, admitisse na solução de

100 a 362 valores singulares.

38

Tabela 3.1: Tabela de erro RMS dos parâmetros do dado e do modelo para diferentes

niveis de corte, sem a presença de ruído, no 'Modelo 1'

3.1.2 β = 10−3

Agora, recuperando os modelos com a adição de ruído de nível baixo (com amplitude β =

10−3) no vetor de dados, foram obtidos os seguintes modelos:

Como bem visto nas �guras (3.10), (3.11), (3.12) e (3.13), e con�rmado pela tabela (3.2),

para os modelos invertidos com amplitude de ruí da ordem de 10−3, podemos perceber que

os modelos recuperados foram satisfatórios visualmente, mesmo com os erros dos parâmetros

do modelo nos valores de 43,70%, 34,30%, 28,52% e 42,49%, respectivamente. Entretanto,

na �gura (3.14) já foi possível notar que não houve uma recuperação satisfatória do modelo,

con�rmado pelo erro dos parâmetros do modelo que foi de 333,34%. O valor de corte que

teve o erro mínimo foi o σcorte = 10−3, onde entraram na solução 351 valores singulares,

como comprovação disso, a imagem (3.13), de todos os modelos recuperados, visualmente é

o que mais se aproxima do modelo verdadeiro, �gura (3.3).

39

(a) (b)

Figura 3.10: 'Modelo 1' recuperado, com ruído de amplitude β = 10−3, onde o valor

do σcorte = 103. A escala em cores representa o valor do coe�ciente

de atenuação (α) em m−1. Em (a) a faixa de α varia de acordo com

o modelo recuperado e em (b) varia em comparação com o modelo

verdadeiro.

(a) (b)





verdadeiro.

A presença do ruído, mesmo que em baixo nível, in�uência na recuperação do modelo

e no aumento do erro dos parâmetros do modelo. A partir do valor de σcorte = 10−7, onde

foram incluídos 354 valores singulares, a recuperação do modelo já não foi mais satisfatórias

e o erro dos parâmetros do modelo aumentaram signi�cativamente. A solução ideal para

40

(a) (b)


do σcorte = 10−3. A escala em cores representa o valor do coe�ciente



verdadeiro.

(a) (b)





verdadeiro.

41

(a) (b)





verdadeiro.

esse método proposto, admitiria 354 valores singulares na solução. Entretanto, é possivel

recuperar modelos de modo satisfatório, através da análise dos erros dos parâmetros do

modelo, qualquer modelo que admitisse entre σcorte = 104 e σcorte = 10−6, ou seja, admitisse

na solução de 100 a 352 valores singulares.

42


niveis de corte, com o nível de ruído β = 10−3, no 'Modelo 1'

3.1.3 β = 10−1

Recuperando os modelos com a adição de ruído de nível médio a alto (com amplitude β =

10−1) no vetor de dados, foram obtidos os seguintes modelos:

Como bem visto nas �guras (3.15) e (3.16), e con�rmado pela tabela (3.3), para os

modelos invertidos com amplitude de ruí da ordem de 10−1, podemos perceber que os modelos

recuperados foram satisfatórios visualmente, mesmo com os erros dos parâmetros do modelo

nos valores de 43,87% e 71,28%, respectivamente. Entretanto, na �gura (3.17) já foi possível

notar que não houve uma recuperação satisfatória do modelo, con�rmado pelo erro dos

parâmetros do modelo que foi de 521,03%. O valor de corte que teve o erro mínimo foi o

σcorte = 102, onde entraram na solução 259 valores singulares.

A presença do ruído, mesmo que de médio a alto nível, in�uência na recuperação do

modelo e no aumento do erro dos parâmetros do modelo. A partir do valor de σcorte =

43

(a) (b)





verdadeiro.

(a) (b)





verdadeiro.

10−1, onde foram incluídos 344 valores singulares, a recuperação do modelo já não foi mais

satisfatórias e o erro dos parâmetros do modelo aumentaram signi�cativamente. A solução

ideal para esse método proposto, admitiria 259 valores singulares na solução. Entretanto, é

possivel recuperar modelos de modo satisfatório, através da análise dos erros dos parâmetros

44

(a) (b)





verdadeiro.

do modelo, qualquer modelo que admitisse entre σcorte = 104 e σcorte = 100, ou seja, admitisse


45



3.2 Modelo 2

Esse modelo foi adaptado de Oliveira (2010), representado na �gura (3.18) é estruturalmente

descrito por um anticlinal assimétrico de origem tectônica que pode ter sofrido mecanismos

de deformação anti e pós-deposicional, não obstante o objetivo da inversão é recuperar bem

as estruturas sem realizar interpretação geológica. Tais situações geolólicas são de rele para

exploração de petróleo por gerar trapas, mais precisamente trapas estruturais dominadas por

dobras. As trapas são arranjos geométricos de rochas que permitem signi�cante acumulação

de hidrocarbonetos em subsuperfície. A trapa muitas vezes inclui a rocha reservatório que

armazena o hidrocarboneto e a rocha selante que impede a migração para fora do reservatório

(Biddle e Wielchowsky, 1994). Neste modelo, a nossa rocha reservatório está representada

pela camada de arenito poroso e permeável, apresentando possível acumulação de hidrocar-

boneto e a rocha selante representada pelo folhelho impermeável. Os valores dos coe�cientes

46

de atenuação foram: 2, 4 × 10−3m−1 para as camadas de folhelho; 3, 0 × 10−3m−1 para ca-

madas de arenito; 3, 3× 10−3m−1 para camadas de arenito reservatório; 2, 4× 10−3m−1 para

camada de folhelho com caráter selante; 1, 5×10−3m−1 para camadas de arenito reservatório

que contém hidrocarboneto nos seu espaço poroso e 2, 6×10−3m−1 para a camada de folhelho

gerador.

O modelo é limitado lateralmente por poços com 30 fontes num poço e 30 receptores

no outro poco, de modo a se ter 900 raios ou 900 equaçoes, consequentemente o problema

é dito sobredeterminado. A discretização se deu com 800 blocos quadrados com coe�cientes

de atenuação constante em cada bloco o quais possuem formato quadrado de dimensão 20

m, sendo 20 blocos na horizontal e 40 blocos na vertical constituindo um modelo com 800 m

na vertical e 400 m na horizontal.

Figura 3.18: 'Modelo 2' verdadeiro, sendo a escala em cores os valores do coe�ciente

de atenuação em m−1.

Depois de realizado o processo da inversão, o grá�co da amplitude do valor singular e

da derivada do valor singular foram obtidas, assim como no Modelo 1, com o intuito de ter

uma noção dos valores singulares e, a partir disso, fazer uma seleção dos valores singulares

que entrarão na solução do problema.

Assim como no Modelo 1, para analisar a qualidade da inversão, foi feita a inversão do

modelo verdadeiro com diferentes valores de corte de seleção de valores singulares (σcorte), que

variaram de 10−15 até 105 com um incremento de 101, em diferentes níveis de ruído variando

a sua amplitude (β). Com isso foi feita a tabela (3.2) que contém os valores de σcorte, de

amplitude β do ruído e o erro RMS dos parâmetros de dados e o erro dos parâmetros do

modelo.

Foram recuperados alguns modelos em diferentes níveis de ruído (β), variando o valor

do corte de valores singulares (σcorte).

47

(a)

(b)

Figura 3.19: Grá�co (a) da amplitude e (b) da derivada do valor singular do 'Modelo

2', em escala monolog.

3.2.1 β = 0

Os modelos recuperados, sem a adição de ruído no vetor de dados, foram:

Como bem visto nas �guras (3.20), (3.21), (3.22), (3.23), (3.24) e (3.25), e con�rmado

pela tabela (3.4), para os modelos invertidos sem ruído, podemos perceber que os modelos

recuperados foram satisfatórios visualmente, com os erros dos parâmetros do modelo nos

valores de 11,57%, 7,84%, 5,00%, 3,12%, 1,95% e 2,16%, respectivamente. O valor de corte

que teve o erro mínimo foi o σcorte = 10−10, onde entraram na solução 749 valores singulares.

A partir do valor de σcorte = 10−12, onde foram incluídos 782 valores singulares, a

recuperação do modelo já não foi mais satisfatórias e o erro dos parâmetros do modelo au-

mentaram signi�cativamente, e isso pode ser notado na �gura (3.19(a)) onde a amplitude

do valor singular após esse valor atrapalharam a solução. A solução ideal para esse método

48

(a) (b)



Em (a) a faixa de α varia de acordo com o modelo recuperado e em


(a) (b)



Em (a) a faixa de α varia de acordo com o modelo recuperado e em


proposto, admitiria 749 valores singulares na solução. Entretanto, conseguiria recuperar mo-

delos de modo satisfatório, através da análise dos erros dos parâmetros do modelo, qualquer

modelo que admitisse entre σcorte = 104 e σcorte = 10−11, ou seja, admitisse na solução de

49

(a) (b)



m−1. Em (a) a faixa de α varia de acordo com o modelo recuperado

e em (b) varia em comparação com o modelo verdadeiro.

(a) (b)



m−1. Em (a) a faixa dos α varia de acordo com o modelo recuperado


255 a 750 valores singulares.

50

(a) (b)





(a) (b)





51


niveis de corte, sem a presença de ruído, no 'Modelo 2'

3.2.2 β = 10−3

Os modelos recuperados, com a adição de ruído (com amplitude β = 10−3) no vetor de

dados, foram:

Como bem visto nas �guras (3.26), (3.27) e (3.28), e con�rmado pela tabela (3.5), para

os modelos invertidos com amplitude de ruí da ordem de 10−3, podemos perceber que os

modelos recuperados foram satisfatórios visualmente, mesmo com os erros dos parâmetros

do modelo nos valores de 11,57%, 7,93% e 14,77%, respectivamente. Entretanto, na �gura

(3.29) já foi possível notar que não houve uma recuperação satisfatória do modelo, con�rmado

pelo erro dos parâmetros do modelo que foi de 178,82%. O valor de corte que teve o erro

mínimo foi o σcorte = 10−1, onde entraram na solução 707 valores singulares.

52

(a) (b)





verdadeiro.

(a) (b)





verdadeiro.

53

(a) (b)





verdadeiro.

(a) (b)



de atenuação (α) em m−1. Em (a) a faixa dos α varia de acordo com


verdadeiro.

54



A presença do ruído, mesmo que em baixo nível, in�uência na recuperação do modelo

e no aumento do erro dos parâmetros do modelo. A partir do valor de σcorte = 10−4, onde

foram incluídos 735 valores singulares, a recuperação do modelo já não foi mais satisfatórias

e o erro dos parâmetros do modelo aumentaram signi�cativamente. A solução ideal para

esse método proposto, admitiria 707 valores singulares na solução. Entretanto, é possivel

recuperar modelos de modo satisfatório, através da análise dos erros dos parâmetros do

modelo, qualquer modelo que admitisse entre σcorte = 104 e σcorte = 10−3, ou seja, admitisse


55

3.2.3 β = 10−1

Os modelos recuperados, com a adição de ruído (com amplitude β = 10−3) no vetor de

dados, foram:

(a) (b)





verdadeiro.

Como bem visto na �gura (3.30), e con�rmado pela tabela (3.6), para o modelo inver-

tido com amplitude de ruí da ordem de 10−1, podemos perceber que o modelo recuperado

foi satisfatório visualmente, com os erros dos parâmetros do modelo no valor de 13,22%.

Entretanto, na �gura (3.31) já foi possível notar que não houve uma recuperação satisfatória

do modelo, con�rmado pelo erro dos parâmetros do modelo que foi de 118,71%. O valor

de corte que teve o erro mínimo foi o σcorte = 103, onde entraram na solução 475 valores

singulares, como comprovação disso, a imagem (3.30), visualmente, se aproxima do modelo

verdadeiro, �gura (3.18).

A presença do ruído, mesmo que de médio a alto nível, in�uência na recuperação do

modelo e no aumento do erro dos parâmetros do modelo. A partir do valor de σcorte =

10−1, onde foram incluídos 707 valores singulares, a recuperação do modelo já não foi mais

satisfatórias e o erro dos parâmetros do modelo aumentaram signi�cativamente. A solução

ideal para esse método proposto, admitiria 475 valores singulares na solução. Entretanto, é

possivel recuperar modelos de modo satisfatório, através da análise dos erros dos parâmetros

do modelo, qualquer modelo que admitisse entre σcorte = 104 e σcorte = 101, ou seja, admitisse


56

(a) (b)





verdadeiro.

57



CAPÍTULO 4

Conclusões

A tomogra�a de atenuação é um problema inverso considerado mal posto. Para ameni-

zar este mal-condicionamento e obter uma solução satisfatória foi utilizado uma metodologia

em que se variava o número de corte dos valores singulares, assim, geralmente,cortando os

pequenos valores singulares.

As simulações numéricas efetuadas para os dois modelos, com diferentes níveis de ruiído

, realizadas com a �nalidade de comparação, permitiram validar as metodologias que foram

empregadas neste trabalho, que se baseiam em uma extração visual (grá�co da amplitude e

grá�co da derivada dos valores singulares) e numérica (calculando o erro RMS dos parâmetros

do modelo) do ponto de corte dos valores singulares com o intuito de minimizar, na medida

do possível, o erro do modelo. Os critérios se mostraram coerentes e tem com uma relação

entre eles.

O método de corte de valores singulares foi satisfatório, mas ele não é extremamente

preciso, pois ele seleciona uma faixa de valores para ser feito o corte, e não o valor exato

do número de valores singulares que devem entrar na solução do sistema. Dessa maneira,

é sugerido em trabalhos futuros que seja feito uma seleção do número de valores singulares

com o intuito de ter uma técnica mais precisa na inversão pelo método da decomposição em

valores singulares na tomogra�a de atenuação.

Os resultados das simulações numéricas para o 'Modelo 1' puderam validar o método

para situações com dados ruidosos, assim como livre de ruídos. O 'Modelo 2' foi a tentativa

de simular um meio geológico um pouco mais complexo. Em ambos os casos, foram obtidos

resultados satisfatórios, entretanto quanto mais ruído era adicionado, menor era o número

de valores singulares que entravam na solução, com isso, comprometendo a qualidade da

recuperação do modelo estimado. Isso �cou claro nas imagens em que a recuperação não era

coerente com o modelo verdadeiro para o mesmo valor de corte de valores singulares para

diferentes tipos de ruído.

A função de erro entre os parâmetros do modelo não pode ser considerada como critério

devido àfalta de conhecimento do campo verdadeiro de fatores de atenuação. Contudo,

mostrou ser importante para que pudéssemos validar a metodologia empregada.

58

59

A inversão por meio da decomposição em valores singulares foi validada para a tomo-

gra�a de atenuação de modo satisfatório.

Agradecimentos

Gostaria de agradecer aos meus pais, Julio Cesar Rocha Fontes e Ivete Helena Lira da

Costa Fontes, por terem me dado todo o apoio em minhas escolhas e o suporte necessário

em vários aspectos.

Aos meus irmãos, Julio Filho e Paloma Helena, pelo companheirismo e cumplicidade.

A toda minha família.

A minha namorada e parceira, Michelle, por todo seu apoio não só durante o curso,

mas ao longo desse caminho que começamos a trilhar juntos. A minha sogra, e grande amiga

Virginia, pelo suporte e todos os conselhos.

A Amin pela oportunidade, atenção, conhecimentos transmitidos e por toda a disponi-

bilidade durante esses anos.

Ao professores Thierry e Vilar por terem aceitado prontamente o convite para participar

da banca deste trabalho.

Aos meus amigos de curso, Rodrigo (Routo), Paulo Augusto (Curió), Leonardo (Mo-

cita), Daniel (Cavanha) e Wilker (Duz¿o) que percorreram esse árduo caminho ao longo

desses 4 anos, por inumeras noites perdidas de estudos e resenhas. E a Rafael Manenti, por

todo o apoio na parte computacional, inumeras caronas, dicas e estudos.

A todos os professores que contribuiram nessa minha graduação.

Aos funcionários, que sempre nos propiciaram um bom ambiente de estudos e trabalho.

A Rede CT-PETRO de Geofísica de Exploração - Fase 3 - Imageamento Sísmico

sob Quebra de Plataforma Continental - Projeto Complementar pela bolsa (de 11/2011

a 07/2012).

A FAPESB/PIBIC/UFBA pela bolsa de iniciação cientí�ca (de 08/2012 a 02/2014) no

projeto de Imageamento Tomográ�co aplicado ao Armazenamento Geológico de CO2.

60

Referências Bibliográ�cas

Andersen, A. H. e Kak, A. C. (1982) Digital ray tracing in two-dimensional refractive

�elds, Journal of Acoustical Society of America, 72(5):1593-1606.

Biddle, K. T. e Wielchowsky, C. C. (1994) Hydrocarbon Traps, In: The Petroleum

System - From Source to Trap, vol. M60, cap. 13, Exxon Exploration Company, Houston,

Texas, U.S.A.

Levenberg, K. (1944) A method for the solution of the certain non-linear problems in

last squares, Quartely of Applied Mathematics, 2, 164-168.

Marquardt, D.W. (1963) An algorithm for least-squares of nonlinear parameters, Jornal

of Society for Industrial and Applied Mathematics, 11, 431-441.

Menke, W. (1984) Geophysical Data Analysis: Discrete Inverse Theory, Academic Press,

Orlando, Florida, 260.

Moreira, Z. M. (2009) Caracterização do Coe�ciente de Atenuação de Ondas Elásticas

em Rochas Sedimentares. Dissertação de Mestrado. Universidade Federal do Rio de Janeiro,

Rio de Janeiro, RJ.

Oliveira, N. P. (2010). Seleção do Parâmetro de Regularização Através da Validação

Cruzada Generalizada: Aplicações em Geotomogra�a Acústica e Eletromagnética. Trabalho

de Graduação. Universidade Federal da Bahia, Salvador, BA.

Plessix, R. (2006) Estimation of velocity and attenuation coe�cient maps from crosswell

seismic data. Geophysics, 71: S235-S240.

Quan, Y. e Harris, J. M. (1997) Seismic attenuation tomography using the frequency

shift method. Geophysics, 62: 895-905.

Silva, C. J. M. G. (2009) Critérios de Seleção de Valores Singulares em Tomogra�a

Acústica de Tempos de Trânsito. Trabalho de Graduação. Universidade Federal da Bahia,

Salvador, BA.

Terra, F. A. (2007) Aplicação da Curva L em Problemas Inversos: Metodologias de

Extração do Parâmetro Ótimo de Regularização. Trabalho de Graduação. Universidade

Federal da Bahia, Salvador, BA.

61

Documents

TOMOGRAFIA DE ATENUACÃO ATRAVÉS DA DECOMPOSICÃO EM … · pada em duas classes: tomogra a cinemática, que utiliza como dado de entrada os tempos de trânsito entre cada par fonte