INVERSÃO DE ATRIBUTOS SÍSMICOS PARA CARACTERIZAÇÃO … · Figura 19 Faixas de ângulos consideradas para a inversão no modelo =5%... 51 Figura 20 Faixas de ângulos consideradas

PETRE CANDELERO GRISCENCO

INVERSÃO DE ATRIBUTOS SÍSMICOS PARA CARACTERIZAÇÃO EMONITORAMENTO DE RESERVATÓRIOS

São Paulo

2015

Versão Corrigida. O original encontra-se disponível na Unidade.



Dissertação apresentada ao Instituto de Astro-nomia, Geofísica e Ciências Atmosféricas daUniversidade de São Paulo para obtenção dotítulo de Mestre em Ciências

São Paulo

2015



Dissertação apresentada ao Instituto de Astro-nomia, Geofísica e Ciências Atmosféricas daUniversidade de São Paulo para obtenção dotítulo de Mestre em Ciências

Área de concentração:Geofísica Aplicada

Orientadora:Liliana Alcazar Diogo

São Paulo

2015

AGRADECIMENTOS

Agradeço primeiramente a Deus. Agradeço também ao meu pai Paulo e a minha mãe Edna,

por todas as oportunidades e por todo o apoio, imensuráveis ao longo dos meus estudos e de

toda minha vida, a eles dedico essa dissertação. Aos meus irmãos, Yuri e Erick, pela amizade,

exemplos e conselhos, mas principalmente pela oportunidade de compartilharmos uma vida.

À Milena por tudo o que ela é e será. Aos amigos de graduação e pós-graduação, Gabriel,

Thaenan, João, Paulo, André, Bruno, Lucas e Henrique, pelo companheirismo, amizade e

conhecimento compartilhado nas várias discussões, úteis ou inúteis.

À Professora Liliana Alcazar Diogo, pela oportunidade, orientação e discussões. Aos co-

legas de corredor, aos professores, em especial ao Professor Renato, e a todo o pessoal do

Departamento por providenciarem um ótimo ambiente de estudo e trabalho.

E finalmente à Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES)

pela bolsa concedida.

RESUMO

Este trabalho tem como objetivo implementar uma inversão sísmica quantitativa e não linear

para investigar a relação entre as propriedades petrofísicas de um reservatório de hidrocarbo-

netos com os seus atributos sísmicos. Utilizando dados sísmicos sintéticos, foi invertido o

coeficiente de reflexão da onda P (RPP ) no topo de um modelo de reservatório com o obje-

tivo de recuperar os parâmetros de velocidade da onda P, da onda S e a densidade na camada

reservatório considerada (VP , VS e ⇢ respectivamente), para tal foi empregada a equação exata

de Zoeppritz. Também agregou-se informação da inversão do tempo de trânsito duplo da re-

flexão normal (t0) na base do reservatório e da inversão da distância crítica (xc) no topo do

reservatório objetivando-se uma melhora na estimativa de VP e consequente melhora na esti-

mativa conjunta dos 3 parâmetros (VP , VS, ⇢). Foram realizados testes quanto à composição

da Função Objetivo envolvendo mais do que um atributo sísmico e da influência da inclusão

de diferentes faixas de ângulos dos coeficientes de reflexão na inversão. Além disso, foram

desenvolvidas análises referentes à performance do algoritmo de inversão, elaborado a partir de

uma implementação do Algoritmo Genético, e foi avaliado como erros na estimativa de alguns

dos parâmetros da camada selante assumidos como conhecidos influencia no desempenho da

inversão. Por fim, utilizando a teoria de Gassmann, foi verificada a viabilidade da recuperação

da variação da saturação de óleo no reservatório em diferentes cenários de porosidade com o

intuito de monitorar a produção do reservatório.

A incorporação do atributo da distância crítica junto com a informação dos coeficientes de

reflexão mostrou resultados promissores, mas são necessários novos testes para recomendar o

seu uso em dados reais. A faixa de ângulos abrangendo a região do ângulo crítico ofereceu

maior sensibilidade e melhores resultados na inversão. Os reservatórios mais porosos são menos

suscetíveis a existência de erros nos parâmetros da camada selante. A recuperação da satu-

ração de óleo só foi possível na situação ideal; pequenas variações nos parâmetros invertidos

foram amplificadas fortemente na determinação da saturação, o que reforça a importância da

realização de estudos numéricos controlados para a inversão de parâmetros físicos da rocha,

com o intuito de avaliar quais fatores exercem maior influência na imprecisão nos resultados.

Palavras-chave: Sísmica 4D. Inversão Sísmica. Petrofísica. Geofísica Aplicada. Sísmica

Time-Lapse.

ABSTRACT

This work aims to implement a quantitative and nonlinear seismic inversion to investigate

the relationship between the petrophysical properties of a hydrocarbon reservoir and its seismic

attributes. Using synthetic seismic data, the P wave reflection coefficient (RPP ) on the top

of a reservoir model was inverted using the exact Zoeppritz equation. The parameters to be

estimated were the P-wave and S-wave velocities and the density of the reservoir layer (VP , VSe ⇢ respectively). It was also added information in inversion of the two-way-traveltime reflection

(t0) on the bottom of the reservoir and of the critical distance (xc) on the top of the reservoir

in order to improve the estimation of VP and consequently improve the joint estimation of the

3 parameters (VP , VS, ⇢). Some tests were conducted about the composition of the Objective

Function for more than one seismic attribute, and about the influence of using different ranges

of incident angles in reflection coefficient inversion. It was also analyzed the performance

of the inversion algorithm, which was based on an implementation of the genetic algorithm.

In addition, it was evaluated how errors in some cap rock estimated parameters, which were

assumed to be known, influences the inversion results. Finally, using the Gassmann theory,

the feasibility of recovering the oil saturation variation in the reservoir for different porosity

scenarios was verified in order to monitor reservoir’s production.

The incorporation of the critical distance attribute along with the information of the re-

flection coefficients showed promising results, but it needs further tests to recommend its use

on real data. The range of angles including the critical angle region showed higher sensitivity

and better results in inversion. The more porous reservoirs was less susceptible to the errors

in the cap rock parameters. The recovery of oil saturation was only possible in the ideal

situation; small variations in the inverted parameters were strongly amplified in determining

the saturation, which reinforces the importance of carrying out controlled numerical studies

for the inversion of the physical parameters of the rock, in order to assess which factors have

more influence in the inaccuracy of the results.

Keywords: 4D Seismic. Seismic Inversion. Rock Physics. Applied Geophysics. Time Lapse

Seismic

LISTA DE ILUSTRAÇÕES

Figura 1 Meio poroso saturado e seus diferentes constituintes. . . . . . . . . . . 16

Figura 2 Módulos de bulk efetivos de uma mistura entre 2 fases utilizando as

médias de Voigt e de Reuss. (Adaptado de (MAVKO; MUKERJI;

DVORKIN, 2009)) . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Figura 3 Partição da energia do campo de onda na interface entre os meios. . . 23

Figura 4 Tempos de chegada para um refletor horizontal. . . . . . . . . . . . . 25

Figura 5 Modelo geométrico do reservatório a ser monitorado. . . . . . . . . . 26

Figura 6 Variação da distância crítica em função da velocidade na camada re-

servatório. (Adaptado de (RODRÍGUEZ, 2013)) . . . . . . . . . . 27

Figura 7 Esquematização do problema direto. . . . . . . . . . . . . . . . . . . 28

Figura 8 Esquematização do problema inverso. . . . . . . . . . . . . . . . . . . 28

Figura 9 Esquematização do uso da FOBJ . . . . . . . . . . . . . . . . . . . . . 29

Figura 10 Esquema do modelo de reservatório. . . . . . . . . . . . . . . . . . . 33

Figura 11 Velocidade da onda P em função da saturação para diferentes �’s. . . 34

Figura 12 Velocidade da onda S em função da saturação para diferentes �’s. . . 35

Figura 13 Densidade em função da saturação para diferentes �’s. . . . . . . . . . 35

Figura 14 Dados gerados - Coeficientes de Reflexão (RPP ) em função do ângulo

de incidência (✓), distância crítica de refração (xc) e tempo normal

de reflexão (t0) - para os parâmetros de controle (Tabela 3) nos

dois cenários de porosidade. . . . . . . . . . . . . . . . . . . . . . 37

Figura 15 Análises dos resultados da inversão - � = 5% (ver texto). . . . . . . . 41

Figura 16 Análises dos resultados da inversão - � = 20% (ver texto). . . . . . . 42

Figura 17 Dados da Figura 14 com adição de ruído. . . . . . . . . . . . . . . . . 43

Figura 18 Soluções da inversão para os dados sintéticos com ruído. . . . . . . . . 45

Figura 19 Faixas de ângulos consideradas para a inversão no modelo � = 5%. . . 51

Figura 20 Faixas de ângulos consideradas para a inversão no modelo � = 20%. . 52

Figura 21 Velocidade da onda P em função da saturação em detalhe. . . . . . . 54

LISTA DE TABELAS

Tabela 1 Valores do modelo de reservatório: parâmetros elásticos (módulos de

bulk e de cisalhamento e densidades) e profundidades do topo (ht)

e da base (hb). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Tabela 2 Valores do folhelho selante: velocidades das ondas P e S e densidade. . 33

Tabela 3 Parâmetros de controle para os dois cenários de porosidade. . . . . . . 38

Tabela 4 Cenários de execução do algoritmo de inversão (I) e resultados dos

critérios adotados para a análise de performance da inversão. (ver

texto) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Tabela 5 Resultados das análises de performance para as inversões a partir de

dados ruidosos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

Tabela 6 Análises da composição da FOBJ a partir de dados ideais: erro mínimo

encontrado durante a execução do Algoritmo de Inversão (I). . . . 48

Tabela 7 Contribuição de cada componente do erro na composição total de ✏

para cada FOBJ a partir de dados ruidosos. . . . . . . . . . . . . . 49


para cada FOBJ a partir de dados ruidosos apenas no atributo de

amplitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50


para cada RPP na inversão das amplitudes em diferentes faixas de

ângulos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Tabela 10 Influência, na recuperação dos parâmetros do reservatório, de erros na

estimativa prévia de cada parâmetro da camada selante. . . . . . 53

Tabela 11 Variações de saturação no reservatório. . . . . . . . . . . . . . . . . . 56

Tabela 12 Resultados das recuperações da variação de saturação (�S1), onde

os índices "i" e "f" referem-se, respectivamente, aos estados ini-

cial e final; e ✏S corresponde ao erro (dado pela equação 30) nos

resultados da inversão. . . . . . . . . . . . . . . . . . . . . . . . . 56





LISTA DE ABREVIATURAS E SIGLAS

GA Algoritmo Genético (Genetic Algorithm)

I Algoritmo de Inversão

M-R Média de Reuss

M-V Média de Voigt

M-VRH Média de Voigt-Reuss-Hill

LISTA DE SÍMBOLOS

VP Velocidade da onda compressional

VS Velocidade da onda cisalhante⇢ Densidade

K Módulo de bulkµ Módulo de cisalhamento� Porosidade

S Saturação de óleo na rocha reservatório

Z Impedância acústica

EI Impedância elástica

RPP Coeficiente de reflexão da onda PPt0 Tempo normal de reflexãoxc Distância críticaht Profundidade até o topo do reservatório

hb Profundidade até a base do reservatório

SUMÁRIO

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

2 FíSICA DE ROCHAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.1 MEIOS POROSOS E PARÂMETROS ELÁSTICOS . . . . . . . . . . . . . . . . 16

2.1.1 Média de Voigt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.1.2 Média de Reuss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.1.3 Média de Voigt-Reuss-Hill . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.2 EQUAÇÕES DE GASSMANN . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.2.1 Módulos elásticos da rocha seca . . . . . . . . . . . . . . . . . . . . . . . 202.2.2 Módulo de bulk do fluido saturante . . . . . . . . . . . . . . . . . . . . . 212.2.3 Densidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3 ATRIBUTOS SíSMICOS . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1 AMPLITUDE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2 TEMPOS DE CHEGADA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.3 DISTÂNCIA CRíTICA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4 O PROBLEMA DE INVERSÃO EM GEOFíSICA . . . . . . . . 284.1 PROBLEMA DIRETO E PROBLEMA INVERSO . . . . . . . . . . . . . . . . . . 28

4.2 FUNÇÃO OBJETIVO E O PROBLEMA DE OTIMIZAÇÃO . . . . . . . . . . . . 29

4.3 ALGORITMO GENÉTICO E ESTRATÉGIA DE INVERSÃO . . . . . . . . . . . . 30

5 MODELO E GERAÇÃO DE DADOS . . . . . . . . . . . . . . . . . 325.1 MODELO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

5.2 GERAÇÃO DE DADOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

6 ANÁLISES E RESULTADOS . . . . . . . . . . . . . . . . . . . . . . . 366.1 PERFORMANCE DO ALGORITMO DE INVERSÃO . . . . . . . . . . . . . . . . 36

6.2 ANÁLISE DA COMPOSIÇÃO DAS FUNÇÕES OBJETIVO . . . . . . . . . . . . 46

6.3 ANÁLISES PARA FAIXAS DISTINTAS DE ÂNGULOS . . . . . . . . . . . . . . . 51

6.4 ANÁLISE DA INFLUÊNCIA DE ERROS NOS PARÂMETROS ASSUMIDOS COMO

CONHECIDOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

6.5 SENSIBILIDADE DO ALGORITMO À VARIAÇÕES DE SATURAÇÃO . . . . . . 54

7 CONCLUSÕES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

Referências . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

12

1 INTRODUÇÃO

Tradicionalmente, na indústria de óleo e gás, a Geofísica sempre desempenhou papel funda-

mental na exploração de novos campos de hidrocarbonetos através do uso do método sís-

mico. Nas últimas décadas, entretanto, vem crescendo o uso do método sísmico na parte

pós-exploratória, no desenvolvimento dos campos, principalmente na caracterização e no mo-

nitoramento de campos já em fase de produção através da chamada "Sísmica 4D" ou "Sísmica

time-lapse".

O método utilizado pela Sísmica 4D consiste em se efetuar aquisições sísmicas (2D ou

3D) em diferentes momentos da vida do reservatório e comparar as aquisições buscando

pelas diferenças entre elas. Deve-se, entretanto, buscar apenas as diferenças provenientes

de variações petrofísicas dentro do reservatório e não diferenças decorrentes do fato de que as

aquisições foram tomadas em momentos diferentes com condições diferentes. Para vencer esse

problema são utilizadas diversas técnicas para aumentar a repetibilidade entre as aquisições

sísmicas e critérios para avaliar numericamente a repetibilidade dos dados.

As primeiras inversões 4D apresentavam caráter qualitativo, eram utilizadas como uma

ferramenta de auxílio no desenvolvimento dos campos. Nos últimos anos cada vez mais

avança-se no sentido de tornar a sísmica 4D mais quantitativa e ela tem papel cada vez mais

central no desenvolvimento dos campos. (GROCHAU et al., 2013)

Diversos casos de sucesso do uso de sísmica 4D em dados reais já foram observados (MAC-

BETH; FLORICICH; SOLDO, 2006; THORE; HUBANS, 2012; VEIRE; BORGOS; LANDRØ,

2006), entretanto a grande maioria desses casos de sucesso ocorre em reservatórios areníticos

menos rígidos e mais porosos, poucos são os casos de sucesso em reservatórios mais rígidos e

menos porosos.

A maioria dos trabalhos que utilizam a informação das amplitudes das reflexões sísmicas

para estudos 4D, assim como nas análises AVO convencionais, empregam aproximações para as

relações que descrevem a partição de energia do campo de onda nas interfaces em subsuperfície

(AKI; RICHARDS, 1980; CASTAGNA; BACKUS, 1993). Essas aproximações são válidas para

ângulos de incidência abaixo do ângulo crítico, região em que o coeficiente de reflexão apresenta

menor sensibilidade para as variações da impedância das camadas. A equação de Zoeppritz

para o coeficiente de reflexão permitiria incorporar a informação da onda refletida em grandes

ângulos de incidência, principalmente ao redor do ângulo crítico, aonde a sensibilidade às

variações de impedância do meio é maior (LANDRØ et al., 2004; DICEZARE; DIOGO,

2013).

Ao se estudar reservatórios mais rígidos e menos porosos, é usual que a impedância acústica

do reservatório seja maior do que a da rocha selante. Nesses casos, existe o ângulo crítico

13

de incidência na interface, e além das variações das amplitudes refletidas devido às mudanças

relativas da produção do reservatório, também pode-se procurar extrair informação sobre as

mudanças de velocidade da onda P na rocha reservatório, a partir da variação da distância

crítica de refração (afastamento fonte-receptor associado ao ângulo crítico de incidência na

interface). Em geral, a onda refletida apresenta amplitude máxima quando incide na interface

com ângulo crítico. Desta forma é possível identificar a posição da reflexão crítica.

O primeiro trabalho a sugerir o monitoramento das variações da distância crítica para averi-

guação das mudanças de velocidade no reservatório foi o de Landrø et al. (2004). Os autores

destacam o fato de que pequenas variações da velocidade geram variações bem maiores no

valor da distância crítica.

Para tornar a Sísmica 4D mais quantitativa é imperativo estabelecer a relação teórica entre

os atributos sísmicos e as características petrofísicas do reservatório (condições in-situ do

reservatório) como, por exemplo: saturação de óleo, pressão, temperatura, etc. Neste trabalho,

será dada atenção apenas às variações de saturação de fluidos durante a produção. O método

mais efetivo e mais utilizado para calcular o efeito da substituição de fluidos nas propriedades

sísmicas baseia-se na teoria desenvolvida por Gassmann (1955).

Neste trabalho foi investigado o problema de inversão não linear das amplitudes refletidas

no topo do reservatório com o objetivo de recuperar os parâmetros de velocidade da onda

P, da onda S e a densidade na camada reservatório (VP , VS e ⇢ respectivamente). Também

agregou-se informação da inversão do tempo de trânsito duplo da reflexão normal (t0) na

base do reservatório e da inversão da distância crítica (xc) no topo do reservatório com o

objetivo de uma possível melhora na estimativa de VP e consequente melhora na estimativa

conjunta dos 3 parâmetros (VP , VS, ⇢). Ao introduzir-se a inversão do atributo distância

crítica, procura-se obter benefício ao considerar reservatórios mais rígidos e utilizar a maior

sensibilidade à pequenas variações na velocidade da onda P dentro do reservatório, como

explorado por Landrø et al. (2004).

O algoritmo de inversão utilizado foi desenvolvido a partir de uma implementação do Algo-

ritmo Genético, que utiliza conceitos da biologia evolutiva para a resolução do problema de

otimização obtido ao efetuar-se a inversão dos atributos sísmicos.

Foi considerada uma situação simplificada do modelo de reservatório para o desenvolvimento

de análises da viabilidade e sensibilidade do processo de inversão visando a extração dos

parâmetros da rocha reservatório (VP , VS, ⇢) com o intuito final de obter-se a saturação de

óleo S na rocha reservatório.

A partir do modelo de reservatório considerado, foram estudados diversos cenários de confi-

guração do algoritmo e de formulação do problema de otimização para o estudo da performance

da inversão sísmica dos parâmetros do reservatório. Ao longo do trabalho foram abordados

14

como os seguintes fatores influenciam na performance e no resultado da inversão:

• Diferentes configurações iniciais de parâmetros de controle do algoritmo de inversão;

• Diferentes composições da Função Objetivo utilizando os 3 atributos sísmicos invertidos;

• Influência da introdução de diferentes faixas de ângulos (uma com ângulos corresponden-tes a afastamentos curtos, outra correspondente até afastamentos intermediários e outra

correspondente até afastamentos longos) dos valores de amplitudes considerados para a

inversão.

A presença de ruído e a influência de erros nos parâmetros assumidos como conhecidos também

foram avaliados no resultado final da recuperação dos parâmetros da inversão.

Ao final do trabalho, abordou-se o conceito de monitoramento do reservatório utilizando

a seguinte lógica: foram considerados dois momentos distintos de sua vida produtiva, para

cada um deles foi efetuada uma inversão sísmica recuperando as propriedades do reservatório

(VP , VS, ⇢). Utilizando a equação de Gassmann e alguns conceitos de Física das Rochas, foi

desenvolvida uma conexão entre as propriedades invertidas e a saturação de óleo no reservatório

considerado em cada um desses momentos. A viabilidade desse metódo foi analisada com o

intuito de monitorar o reservatório.

Os conceitos teóricos necessários ao entendimento do problema estudado, são abordados ao

longo dos capítulos 2, 3 e 4 da dissertação. No capítulo 2 são abordados conceitos de Física

de Rochas e apresentada a equação de Gassmann; já no capítulo 3 apresenta-se com mais

detalhe os atributos sísmicos a serem invertidos; no quarto capítulo o problema da inversão

sísmica em Geofísica é apresentado.

Nos capítulos 5 e 6 são apresentados o modelo de reservatório considerado, os resultados

obtidos e as análises efetuadas.

15

2 FÍSICA DE ROCHAS

Ao se estudar mudanças dinâmicas no reservatório (saturação, pressão, temperatura) decor-

rentes de sua produção utilizando dados sísmicos, é de fundamental importância estabelecer

uma firme conexão entre os dados sísmicos e as propriedades petrofísicas do reservatório. Essa

conexão pode ser estabelecida através da teoria da Física de Rochas.

Para utilizar a teoria de Física de Rochas de forma mais quantitativa é fundamental que

também sejam estabelecidas algumas relações empíricas calibradas in situ com dados de poço

e de produção especificamente para cada campo.

A influência da substituição de fluidos no comportamento dos módulos elásticos de um

meio poroso é comumente investigada utilizando a equação de Gassmann (GASSMANN,

1951). A teoria de Gassmann é bem estabelecida e muitos autores (BATZLE; WANG, 1992;

BERRYMAN, 1999; WANG, 2001; RUSSELL et al., 2003; HAN; BATZLE, 2004) têm discutido

as formulações e limitações da substituição de fluidos de Gassmann.

Na substituição de fluido relacionam-se as propriedades elásticas da rocha com seu espaço

poroso preenchido por determinado fluido (por exemplo, óleo) com as novas propriedades

elásticas que a rocha teria caso seus poros estivessem preenchidos por um fluido diferente (por

exemplo, uma mistura água-óleo). O objetivo é relacionar a velocidade extraída da sísmica em

uma determinada condição petrofísica do reservatório (por exemplo: pressão, temperatura,

porosidade, tipo mineral) com as variações de seus módulos elásticos.

A teoria da elasticidade relaciona os módulos elásticos de um meio homogêneo e isotrópico

com as velocidades de propagação das ondas de corpo P (primárias, ou compressionais, ou

longitudinais) e S (secundárias, ou cisalhantes, ou transversais) através de:

VP =

sK + (4/3)µ

⇢(1)

VS =

rµ

⇢(2)

onde K, µ e ⇢ são respectivamente o módulo de bulk, o módulo de cisalhamento e a densidade

do meio em que se propaga a onda.

Uma rápida análise nas equações acima mostra que ao se obter as velocidades e as densidades

de informações extraídas da sísmica como, por exemplo, de uma inversão sísmica, podemos

encontrar o valor dos módulos elásticos e então estabelecer uma conexão com parâmetros

petrofísicos da rocha reservatório através das equações de Gassmann.

16

2.1 MEIOS POROSOS E PARÂMETROS ELÁSTICOS

As equações 1 e 2 são derivadas pelas teorias de propagação de ondas e da elasticidade,

relacionando tensões e deformações dos meios. Considere-se uma rocha porosa, como um

arenito, saturada por um fluido. Devido a diferenças nas propriedades elásticas de cada um

dos constituintes do meio, uma onda propagando-se através dele faz com que o fluido, os

grãos e o esqueleto rochoso sejam submetidos a diferentes graus de deformação.

Rocha Saturada

Mineral

Poros/Fluidos

Esqueleto Rochoso(poros drenados)

Figura 1 - Meio poroso saturado e seus diferentes constituintes.

Nessas condições é interessante estimar um módulo elástico efetivo desse meio, considerando

de alguma maneira os módulos elásticos de todos os seus constituintes. Diversos esforços

tem sido empregados nesse tipo de análise e várias teorias de meio efetivo atualmente são

consideradas para contornar esse problema (MAVKO; MUKERJI; DVORKIN, 2009). A maioria

dos modelos propostos na literatura empenham-se em descrever teoricamente o módulo elástico

efetivo considerando:

1. A fração do volume de cada uma das fases constituintes,

2. O módulo elástico de cada uma das fases constituintes,

3. A disposição geométrica de como as fases constituintes do meio se relacionam entre si.

Entretanto, na prática, a situação de incluir as informações geométricas a respeito da consti-

tuição do meio no modelo teórico, dificilmente é solucionada de maneira adequada. A maioria

das tentativas é direcionada no sentido de simplificação e aproximação dessa geometria e ao

fazê-lo os modelos perdem sua generalidade e previsibilidade tornando-se específicos para as

geometrias incorporadas.

Uma abordagem muito interessante e útil é considerar apenas a fração volumétrica de

cada uma das fases e seus respectivos módulos elásticos. Na literatura, de maneira geral,

consideram-se três tipos de médias: Voigt, Reuss e Voigt-Reuss-Hill. Cada uma dessas médias

17

delineia limites para os módulos elásticos efetivos do meio e associa os módulos elásticos de

cada uma das fases constituintes de maneira diferente, considerando-se a Lei de Hooke da

teoria da elasticidade:

⌧ = M✏ (3)

aonde ⌧ é o esforço aplicado a um meio, ✏ é a deformação decorrente do esforço e M é o

módulo elástico do meio que pode ser tanto o módulo de bulk quanto o de cisalhamento.

2.1.1 Média de Voigt

A Média de Voigt (M-V) considera que o meio é submetido a um esforço ⌧ qualquer de

maneira que todos os seus constituintes sejam submetidos a mesma deformação ✏. Nesse

cenário, na Lei de Hooke considera-se o i-ésimo elemento que compõe o material contribuindo

para o esforço total com (Mi✏) ponderado pela sua fração fi do meio:

⌧ = f1(M1✏) + f2(M2✏) + . . .+ fn(Mn✏)

⌧ = (f1M1)✏+ (f2M2)✏+ . . .+ (fnMn)✏

⌧ =nX

i=1

(fiMi)✏ = MV ✏

portanto o meio composto por n elementos tem seu módulo elástico efetivo dado pela média

de Voigt MV :

MV =nX

i=1

(fiMi) (4)

A média de Voigt conduz aos maiores valores teóricos possíveis do módulo elástico efetivo

do meio e, portanto, é conhecida também como limite superior de Voigt. Entretanto, a

média de Voigt dificilmente representa sistemas físicos reais pois eles nunca podem ser tão

incompressíveis quanto ela prevê (exceto quando existe apenas uma fase no sistema).

2.1.2 Média de Reuss

Na Média de Reuss (M-R) considera-se que todos os elementos são submetidos ao mesmo

esforço, ao invés da mesma deformação como na média de Voigt. Novamente considerando-

se a lei de Hooke, tem-se que o i-ésimo elemento constituinte do material contribui para

a deformação total do meio com uma deformação⇣

⌧Mi

⌘poderada pela fração fi de cada

elemento que constitui o meio. De modo que:

✏ = f1

✓⌧

M1

◆+ f2

✓⌧

M2

◆+ . . .+ fn

✓⌧

Mn

◆

18

✏ =

✓f1M1

◆⌧ +

✓f2M2

◆⌧ + . . .+

✓fnMn

◆⌧

✏ =nX

i=1

✓fiMi

◆⌧ =

⌧

MR

Assim o meio composto por n elementos tem seu módulo elástico efetivo dado pela média

de Reuss MR:

MR =

nX

i=1

✓fiMi

◆!�1(5)

A média de Reuss conduz aos menores valores teóricos do módulo elástico efetivo do meio,

sendo assim o limite inferior. Ela descreve exatamente um sistema no qual todas as fases

são submetidas ao mesmo esforço como por exemplo, grãos sólidos estão em suspensão num

fluido, uma mistura de gases ou de fluidos.

2.1.3 Média de Voigt-Reuss-Hill

A Média de Voigt-Reuss-Hill (M-VRH) é definida como uma simples média aritmética das

médias de Voigt e de Reuss:

MV RH =1

2(MV +MR)

MV RH =1

2

0

@nX

i=1

(fiMi) +

nX

i=1

✓fiMi

◆!�11

A (6)

A média de Voigt-Reuss-Hill é usada quando faz-se necessário estimar o módulo elástico

efetivo do meio ao invés de se obter uma faixa de valores, sendo um resultado estritamente

heurístico. Na Figura 2 compara-se a aplicação das médias de Voigt e Reuss para a obtenção

de um módulo de bulk efetivo Keff em um sistema de duas fases, uma com módulo de bulk

K1 e outra com módulo de bulk K2; a média de Voigt-Reuss-Hill, obviamente, encontraria-se

entre as duas curvas apresentadas.

19

Figura 2 - Módulos de bulk efetivos de uma mistura entre 2 fases utilizando as médias de

Voigt e de Reuss. (Adaptado de (MAVKO; MUKERJI; DVORKIN, 2009))

2.2 EQUAÇÕES DE GASSMANN

As equações de Gassmann relacionam o módulo de bulk de uma rocha, com as propriedades

da sua estrutura (matriz mineral e porosidade) e da dos fluidos contidos em seus poros. O

módulo de bulk de uma rocha saturada é dado segundo a teoria de Gassmann (GASSMANN,

1951; HAN; BATZLE, 2004) por:

Ksat = Kd +K0(1�Kd/K0)2

1� ��Kd/K0 + �K0/Kfl(7)

onde Ksat, Kd, K0, Kfl são, respectivamente, o módulo de bulk da rocha saturada, do

esqueleto rochoso ou da rocha seca (esqueleto drenado de qualquer fluido nos poros), da

matriz mineral, dos fluidos saturantes e � é a porosidade (em fração).

As equações de Gassmann, são válidas quando aplicadas em um meio elástico homogêneo,

isotrópico e monominerálico no qual não há interação química entre o fluido e o esqueleto

rochoso.

Também é importante destacar que o trabalho de Gassmann trata de processos quasi-

estáticos, processos que podem ser tratados como uma infinidade de transições diferenciais

entre estados de equilíbrio do sistema. Dessa forma, também existe a necessidade de equilíbrio

da pressão entre os poros ao longo de todo o esqueleto rochoso durante a passagem do campo

de onda; tal fato implica que quanto menor a frequência do campo de onda, melhor essa

condição é satisfeita. Para a faixa de frequências adquiridas na sísmica profunda a teoria de

Gassmann mostra ótimos resultados (MAVKO; MUKERJI, 1998; BERRYMAN, 1999).

Como ondas cisalhantes não interagem com os fluidos presentes nos poros e uma das con-

dições para a validade das equações de Gassmann é a não interação química entre o fluido e o

20

esqueleto rochoso, o módulo de cisalhamento é independente do fluido nos poros e é mantido

constante durante as substituições de fluido (GASSMANN, 1951; BERRYMAN, 1999).

µsat = µd (8)

aonde µsat é o módulo de cisalhamento da rocha saturada pelo fluido e µd é o módulo de

cisalhamento da rocha seca.

Observando a equação 7 pode-se notar que para o cálculo do módulo elástico da rocha

saturada, é necessário calcular/estimar parâmetros elásticos da rocha seca, da matriz mineral

e dos fluidos saturantes. As densidades, por exemplo, podem simplesmente ser calculadas

como uma média aritmética simples das densidades de cada um dos constituintes ponderada

pelo respectivo volume ocupado pelo constituinte. Já os módulos de bulk requerem o uso das

médias discutidas no capítulo anterior.

Nas seções seguintes, discutem-se fórmulas para calcular o módulo de bulk e densidades da

matriz mineral, do fluido saturante e do esqueleto rochoso (rocha seca).

2.2.1 Módulos elásticos da rocha seca

O módulo de bulk do esqueleto rochoso ou módulo de bulk da rocha seca Kd, é obtido atra-

vés de relações empíricas para cada tipo de rocha. Muitos estudos (HAN; NUR; MORGAN,

1986; HAN; BATZLE, 2004; MAVKO; MUKERJI, 1995; MURPHY; REISCHER; HSU, 1993)

abordam diferentes técnicas e métodos para a a obtenção dessas relações; a grande maioria

deles relacionando o módulo com a porosidade do reservatório.

Em nosso caso, utilizaremos uma relação empírica derivada de um modelo de Voigt mo-

dificado (incorporando o conceito de porosidade crítica) que se aplica muito bem a arenitos

proposto por Han e Batzle (HAN; BATZLE, 2004), dado por:

Kd = K0(1� �/�c) = K0(1� 2.5�) (9)

onde K0 é o módulo de bulk do mineral, � é a porosidade e �c é a porosidade crítica. Nesse

trabalho, será utilizado o valor de �c = 40%, que é um valor típico de porosidade crítica para

arenitos clássicos (NUR et al., 1998).

De maneira análoga, tem-se para o módulo de cisalhamento da rocha seca (HAN; BATZLE,

2004):

µd = µ0(1� �/�c) = µ0(1� 2.5�) (10)

aonde µ0 é o módulo de cisalhamento da matriz mineral, � é a porosidade e �c = 40% é a

porosidade crítica.

21

2.2.2 Módulo de bulk do fluido saturante

Já para o módulo de bulk do fluido, conforme discutido em 2.1.2, utilizaremos a média de

Reuss (equação 5) dos módulos de bulk da água e do óleo:

Kfl =

(1� S)KH2O

+S

Koil

��1(11)

onde KH2O é o módulo de bulk da água, Koil é o módulo de bulk do óleo e S é a saturação

de óleo na rocha reservatório que varia de 0 (100% água) a 1 (100% óleo). Note que a média

calculada na equação assume um fluido bifásico composto apenas de água e óleo. Se houver

uma fase gás, por exemplo, a média pode ser calculada do mesmo modo considerando-se as

frações de saturação de cada fluido.

2.2.3 Densidades

A densidade do fluido ⇢fl é calculada diretamente como uma média aritmética simples das

densidades de cada um das fases constituintes do fluido ponderadas pela respectiva fração de

cada um deles. Então:

⇢fl = S⇢oil + (1� S)⇢W (12)

onde S é a saturação de óleo no reservatório, ⇢oil é a densidade do óleo e ⇢W é a densidade

da água.

Admitindo-se um reservatório homogêneo e com os poros completamente saturados pelo

fluido, com base na densidade do fluido ⇢fl pode-se calcular a densidade média do reservatório

⇢ novamente como uma média aritmética simples ponderada pela fração que cada constituinte

ocupa no reservatório:

⇢ = (1� �)⇢0 + �⇢fl (13)

onde ⇢0 é a densidade da matriz mineral do reservatório e � é a porosidade do reservatório.

22

3 ATRIBUTOS SÍSMICOS

Atributo sísmico é qualquer propriedade extraída do dado sísmico, seja diretamente ou por

manipulação matemática, que isola ou enfatiza alguma característica específica dos dados,

favorecendo a interpretação geofísica e ou geológica dos dados.

Dentre os diversos tipos de atributos, alguns podem ser usados para indicação direta da pre-

sença de hidrocarbonetos, já outros podem ser usados para realçar determinada característica

geológica ou para indicar alguma tendência geral nos dados. Ou seja, dependendo do objetivo

e do tipo de análise devemos escolher atributos sísmicos diferentes.

Em um trabalho com o intuito de monitorar um reservatório a partir de dados sísmicos,

quais seriam os atributos mais adequados para serem utilizados? Deve-se procurar utilizar

atributos que relacionam-se com as propriedades petrofísicas monitoradas da maneira mais

direta e mais sensível possível. Neste trabalho, também leva-se em conta o fato de procurar

utilizar os atributos sísmicos mais convenientes para serem simulados numericamente, uma

vez que não serão utilizados dados reais e sim dados sintéticos.

Em seguida apresentam-se e discutem-se os atributos escolhidos para o desenvolvimento do

trabalho.

3.1 AMPLITUDE

Comumente, um atributo utilizado para extrair algum tipo de informação litológica das

camadas é a amplitude do sinal refletido. Essa amplitude está diretamente relacionada à

partição da energia das ondas sísmicas atravessando uma interface representada por uma

wavelet no traço sísmico.

A partição de energia ocorre devido ao contraste entre as impedâncias das camadas que

definem tal interface. As impedâncias, acústica da onda P (ZP ) e da onda S (ZS), de um

meio são definidas por:

ZP = ⇢VP (14)

ZS = ⇢VS (15)

aonde VP e VS são a velocidade das ondas P e S, respectivamente e ⇢ é a densidade do meio.

23

Figura 3 - Partição da energia do campo de onda na interface entre os meios.

Como pode-se observar na Figura 3, uma frente de onda incidente P (PPi) origina 4 frentes

refletidas/refratadas após atingir a interface: 2 de mesma fase P (refletida (PPr) e refratada

(PPt)) e outras duas convertidas para a fase S (refletida (PSr) e refratada (PSt)).

As equações que governam o fenômeno da partição da energia do campo de onda sísmico

na interface entre dois meios com diferentes propriedades elásticas são as equações de Zoep-

pritz (ZOEPPRITZ, 1919). A partir das condições de contorno na interface: continuidade

dos esforços e das deformações, ambos normais e de cisalhamento, Zoeppritz deduziu um

sistema de 4 equações relacionando as amplitudes das ondas incidente, refletida na interface

e transmitida entre os meios, a partir do qual é possível isolar as 4 incógnitas que definem os

coeficientes de reflexão e transmissão da partição de energia ilustrada na Figura 3.

As equações de Zoeppritz para os coeficientes de reflexão são complicadas para resoluções

estritamente algébricas e não fornecem intuitivamente noções de como se comportam as am-

plitudes em função do ângulo de incidência (CASTAGNA; BACKUS, 1993). De maneira geral,

para se extrair essas noções por trás das equações e simplificar seu tratamento algébrico, são

utilizadas linearizações e aproximações ao se trabalhar com as equações de Zoeppritz na prá-

tica. As aproximações mais conhecidas e utilizadas são as de Aki-Richards (AKI; RICHARDS,

1980) e Shuey (SHUEY, 1985). Essas aproximações, embora facilitem o tratamento e forne-

çam uma melhor percepção física do fenômeno, têm limites de validade e geralmente valem

para baixos contrastes de impedância acústica entre as camadas e para ângulos de incidên-

cia menores do que 30 graus (SHUEY, 1985) e abaixo do ângulo crítico, no caso em que a

velocidade da onda aumenta da camada superior para a inferior.

Entretanto, na região próxima ao ângulo crítico, a sensibilidade para variações da velocidade

da camada inferior é consideravelmente maior (LANDRØ et al., 2004; DICEZARE; DIOGO,

2013). E como procura-se um atributo que apresente forte sensibilidade às propriedades

monitoradas é interessante utilizar as equações exatas de Zoeppritz que contemplam toda a

24

faixa de ângulos de incidência disponíveis.

Optou-se por implementar computacionalmente a equação de Zoeppritz para o coeficiente

de reflexão (RPP ), escrita por Ikelle e Amundsen (IKELLE; AMUNDSEN, 2005), como:

RPP =c1d2 � c2d4d1d2 + d3d4

(16)

aonde

c1 = 2p2�µ(qP1 + qP2)� (⇢1qP2 � ⇢2qP1)

c2 = �p⇥2�µ

�qP1qS2 � p2

�+�⇢

⇤

d1 = 2p2�µ (qP1 � qP2) + (⇢1qP2 + ⇢2qP1)

d2 = 2p2�µ (qS1 � qS2) + (⇢1qS2 + ⇢2qS1)

d3 = p⇥2�µ

�qP1qS2 + p

2��⇢

⇤

d4 = p⇥2�µ

�qP2qS1 + p

2��⇢

⇤

e qP1 =qV �2P1 � p2, qS1 =

qV �2S1 � p2, qP2 =

qV �2P2 � p2 e qS2 =

qV �2S2 � p2 e

p = sin ✓iVP1 .

Os índices 1 e 2 referem-se respectivamente ao meio superior e ao meio inferior à interface

refletora sendo VP , VS, ⇢i e ✓i a velocidade da onda P, a velocidade da onda S, a densidade

do meio i e o ângulo de incidência da onda compressional na interface.

3.2 TEMPOS DE CHEGADA

Outro atributo sísmico utilizado na caracterização do reservatório é o tempo de chegada

duplo (two-way-traveltime - TWT), que é o intervalo de tempo gasto para a onda percor-

rer a distância entre a fonte e a interface refletora e voltar, consequentemente depende da

profundidade e da velocidade média dos meios acima do refletor.

25

Figura 4 - Tempos de chegada para um refletor horizontal.

Na Figura 4 observa-se o esquema geométrico da trajetória dos raios sísmicos e a corres-

pondente curva dos tempos de chegada em função do afastamento x da fonte, dados pela

expressão:

t2(x) =

✓2z

V

◆2+⇣ xV

⌘2= t20 +

⇣ xV

⌘2(17)

aonde z é a profundidade do refletor, V é a velocidade do meio acima da interface refletora

e t0 é o tempo de chegada normal, aquele em que fonte e receptor encontram-se na mesma

posição.

Note que a velocidade é considerada constante no meio acima do refletor. Se o meio acima

do refletor for constituído por diversas camadas, pode-se utilizar a equação da hipérbole (17)

como uma aproximação para os tempos de chegada das reflexões, substituindo a velocidade

da expressão (17) pela velocidade quadrática média (DIX, 1955; SHERIFF; GELDART, 1995)

VRMS dada por:

VRMS =

sPni=1 v

2i ⌧iPn

i=1 ⌧i(18)

aonde n é o número de camadas, vi é a velocidade da i-ésima camada e ⌧i é o tempo simples

ou tempo duplo de trânsito dentro da i-ésima camada.

26

Camada Reservatório

Figura 5 - Modelo geométrico do reservatório a ser monitorado.

Na Figura 5 apresenta-se o esquema geométrico do reservatório a ser monitorado. Como o

intuito é monitorar a velocidade da camada reservatório V2, considerando a velocidade média

V1 acima do topo do reservatório como já conhecida, tem-se o tempo normal de reflexão na

base do reservatório t0 dado por:

t0 =2htV1

+2(hb � ht)

V2(19)

3.3 DISTÂNCIA CRÍTICA

Para as situações geológicas em que a velocidade no reservatório for maior do que na rocha

selante, pode-se lançar mão de uma nova abordagem apresentada por Landrø (LANDRØ et

al., 2004; LANDRØ, 2006) que utiliza um atributo sísmico incomum para a recuperação de

informação da velocidade na camada reservatório: a distância crítica.

A distância crítica é o afastamento fonte-receptor que corresponde ao ângulo crítico de

incidência, aquele a partir do qual não há mais energia transmitida para a camada inferior.

Para um modelo geométrico, como o da Figura 5 a distância crítica é dada por:

xc =2htr⇣V2V1

⌘2� 1

(20)

Na Figura 6 é ilustrado como o ângulo crítico, e consequentemente, a distância crítica,

variam devido a uma mudança na velocidade do reservatório V2.

27

Figura 6 - Variação da distância crítica em função da velocidade na camada reservatório.

(Adaptado de (RODRÍGUEZ, 2013))

Para o caso de várias camadas acima do topo do reservatório, Landro et al. (2004) sugere

uma aproximação para a expressão (20), substituindo V1 por VRMS.

A onda refletida, em geral, apresenta amplitude máxima quando incide na interface com

ângulo crítico. Desta forma é possível identificar a posição da reflexão crítica. Mas, ao invés do

valor da amplitude, será o afastamento, o atributo utilizado para estimar o valor da velocidade

do reservatório.

A distância crítica pode trazer informação altamente sensível à variação de velocidade da

camada reservatório, uma pequena variação na velocidade da camada reservatório V2 conduz

a uma grande variação na distância crítica xc.

28

4 O PROBLEMA DE INVERSÃO EM GEOFÍSICA

Um problema central em geofísica é o de efetuar-se inferências sobre sistemas físicos a partir

de dados observados que, supostamente, sejam causados pelo sistema físico em questão. Esse

problema, amplamente abordado na área de geofísica, é o problema de inversão.

4.1 PROBLEMA DIRETO E PROBLEMA INVERSO

O sistema físico sobre o qual deseja-se obter informações é desconhecido e o objetivo da

inversão é extrair tais informações dos dados. Entretanto, em geral, existe algum tipo de infor-

mação adicional (geralmente proporcionada pela geologia) que fornece uma série de hipóteses

sobre o sistema.

A partir dessas hipóteses, pode-se construir um modelo do sistema físico estudado e um

modelo matemático que descreva como esse modelo gera dados que simulem os dados observa-

dos. O modelo do sistema físico é descrito por parâmetros, e o modelo matemático é descrito

por equações que relacionem os parâmetros aos dados. Portanto, definir as relações funcionais

entre o modelo parametrizado (construído a partir das hipóteses) com dados que simulem os

dados observados através de um modelo matemático é o que define-se por problema direto.

hipóteses modelo dados

Figura 7 - Esquematização do problema direto.

Utilizando os mesmo conceitos, o problema inverso consiste em utilizar os dados observados

e as hipóteses para inferir o modelo paramétrico que origina dados mais semelhantes aos dados

efetivamente observados. Nota-se, entretanto, que o caminho inverso é bem mais complicado.

Na natureza, as relações matemáticas que descrevem como os dados relacionam-se ao modelo

são, na maioria das vezes, não-lineares. Isso, aliado ao fato do modelo ser multiparamétrico,

faz com que existam vários modelos que geram dados semelhantes aos dados observados.

hipóteses

modelo dados

Figura 8 - Esquematização do problema inverso.

29

Mais formalmente pode-se definir o problema direto como:

�!d0 = f (�!m) (21)

aonde�!d0 = (d01, d

02, . . . , d

0n) é um conjunto de n dados que simulam os dados observados,

�!m = (p1, p2, . . . , pk) é um conjunto de k parâmetros que descrevam o modelo e f é umfuncional que relaciona o conjunto de dados ao conjunto de parâmetros.

Já o problema inverso seria o de desenvolver um método para encontrar um modelo�!m0 que

produza dados�!d0 = f

⇣�!m0⌘

mais próximos possíveis dos dados efetivamente observados �!d .Note que ao contrário do que pode parecer em um primeiro momento, o problema inverso,

em geral, não é o de encontrar um funcional f�1 que forneça o modelo �!m a partir dos dadosobservados, pois como mencionado anteriormente, em sistemas físicos presentes na natureza,

as relações matemáticas que descrevem o funcional f são não-lineares, o que, por vezes, torna

impossível definir univocamente e/ou de forma estável o funcional f�1.

4.2 FUNÇÃO OBJETIVO E O PROBLEMA DE OTIMIZAÇÃO

O problema de inversão conforme apresentado na seção anterior pode ser considerado um

problema de otimização. Problemas de otimização são aqueles aonde faz-se necessário en-

contrar estados ótimos, ou seja, a melhor solução possível dentre todas as soluções viáveis.

Um problema de otimização é, por exemplo, encontrar o caminho mais curto entre dois pon-

tos; ou ainda encontrar a melhor geometria para um veículo oferecer a menor resistência ao

deslocamento.

Conforme as definições da seção anterior, o problema de otimização, no caso da inversão,

é o de encontrar um modelo�!m0 que gere dados

�!d0 = f

⇣�!m0⌘

mais próximos possíveis dos

dados �!d efetivamente observados.

dadoscalculadosmodelo

dadoscomparação

Figura 9 - Esquematização do uso da FOBJ .

Para comparar os dados observados �!d e os dados calculados�!d0 é necessário estabelecer um

critério de semelhança entre os dois conjuntos de dados, em problemas inversos, esse critério é

30

estabelecido através de uma função objetivo (FOBJ). De maneira geral, pode-se utilizar uma

função de mínimos quadrados para definir a função objetivo:

FOBJ =

"nX

i=1

⇣di � d

0

i

⌘2#1/2

(22)

aonde di é o i-ésimo dado do conjunto de n dados e d0i é o i-ésimo dado calculado utilizando

o modelo�!m0.

Resolver o problema inverso é, portanto, encontrar o modelo (�!m 0)min que minimiza a funçãoobjetivo, ou, em outras palavras, maximiza a semelhança entre os conjuntos de dados. Ao

encontrar o modelo (�!m 0)min, temos o conjunto de parâmetros que melhor descreve o sistemafísico estudado.

4.3 ALGORITMO GENÉTICO E ESTRATÉGIA DE INVERSÃO

Deve-se, adequar a cada problema estudado uma estratégia de inversão adequada para

minimizar a função objetivo. Neste trabalho, os dados invertidos serão os atributos sísmicos

apresentados no Capítulo 3, mais especificamente os atributos definidos teoricamente pelas

equações 16, 19 e 20. A inversão sísmica a ser efetuada deve valer-se da construção de funções

objetivos adequadas para encontrar o modelo que gera os dados que mais se aproximam dos

dados observados.

Para resolver o problema de otimização das funções objetivo, ou seja, para encontrar o

valor do mínimo dessas funções será utilizado o Algoritmo Genético (Genetic Algorithm)

(GA). O uso do Algoritmo Genético na inversão de problemas inversos em geofísica é rela-

tivamente recente e vem apresentando bons resultados em comparação com outros métodos

tradicionalmente empregados para encontrar os mínimos da FOBJ (BOSCHETTI; DENTITH;

LIST, 1996; LOUIS; CHEN; PULLAMMANAPPALLIL, 1999; SAMBRIDGE; DRIJKONINGEN,

1992).

O uso do algoritmo genético é extremamente adequado à resolução de problemas de oti-

mização não lineares e de minimização global; ao contrário de técnicas de minimização local

(que podem fornecer soluções presas em um mínimo local), o algoritmo genético evita o uso

de qualquer tipo de informação a respeito da curvatura da função objetivo. Isso significa que

derivadas de funções não precisam ser avaliadas e, portanto, pode-se utilizar qualquer tipo de

função objetivo obtendo bons resultados da mesma maneira.

A maioria dos métodos interativos de busca trabalham com apenas um modelo e, através

de perturbações, vão melhorando esse modelo a cada passo. A grande vantagem do algoritmo

genético é trabalhar ao mesmo tempo com uma população de modelos e utilizar conceitos da

31

biologia evolutiva (como heredietariedade, mutação, recombinação e, principalmente , seleção

natural) para "evoluir" a população de modelos a cada "geração" em busca dos "indivíduos

ótimos".

A seguir apresenta-se um resumo de como o algoritmo genético funciona:

• O algoritmo cria aleatoriamente uma população inicial de modelos/indivíduos.

• O algoritmo cria uma sequência de novas populações. A cada passo, a geração atual deindivíduos/modelos é utilizada para criar a próxima seguindo os seguintes passos:

1. Cada membro da população recebe um valor de uma função fitness (quanto menor o

valor de fitness, melhor é o modelo/indivíduo)

2. Membros da população são escolhidos, chamados de "pais", baseados no seu valor de

função fitness

3. Os melhores membros da população, chamados de "elite", são passados automaticamente

para a próxima geração

4. Os "pais" produzem a "prole" através de mudanças aleatórias em um único "pai" (mu-

tação) ou através de uma combinação entre entre as entradas dos vetores associados a

um par de modelos/indivíduos (recombinação).

5. A população atual é então substituída pela elite e pela prole.

• O Algoritmo é encerrado quando algum critério de parada é atingido, em geral quando amudança entre os valores de fitness dos "pais" e da "prole" não varia mais do que um

valor pré-estabelecido.

O algoritmo Genético utilizado neste trabalho foi o fornecido no pacote de otimização do

Matlab1. Todos os demais programas, para os testes, simulações e modelagens, realizados

neste trabalho foram efetuados utilizando esse software.

1Matlab R2014a.

32

5 MODELO E GERAÇÃO DE DADOS

5.1 MODELO

O modelo de reservatório a ser considerado obedece as premissas estabelecidas pelas relações

de Gassmann da Física de Rochas apresentadas no Capítulo 2, tendo em vista que elas serão

utilizadas para a obtenção de diversos parâmetros do modelo para o estudo do reservatório.

Essas premissas são sumarizadas a seguir:

• O reservatório é poroso, isotrópico, elástico, monominerálico e homogêneo;

• Os poros são conectados e estão em equilíbrio de pressão;

• O meio é um sistema fechado sem troca de fluidos através de suas bordas;

• Não há interação química entre os fluidos e o esqueleto rochoso do reservatório.

Os valores dos parâmetros do modelo do reservatório, apresentados na Tabela 1, foram base-

ados em dados de um arenito eocênico da Bacia de Campos (ROMANELLI, 2010), conside-

rando seus poros saturados por apenas duas fases: água e óleo. Considerou-se uma camada

reservatório composta por esses arenitos relativamente limpos com matriz mineral de quartzo,

portanto o índice "0" dos módulos da tabela referem-se ao módulo de bulk e densidade do

quartzo. Já os índices "oil" e "H2O" referem-se as duas fases respectivamente: óleo e a

água.

Parâmetro Valor

K0 36.4GPa

KH2O 2.95GPa

Koil 1.05GPa

µ0 40GPa

⇢0 2.65g/cm3

⇢H2O 1.04g/cm3

⇢oil 0.715g/cm3

ht 0.5km

hb 0.6km

Tabela 1 - Valores do modelo de reservatório: parâmetros elásticos (módulos de bulk e de

cisalhamento e densidades) e profundidades do topo (ht) e da base (hb).

33

A camada selante foi admitida como sendo um folhelho com os valores de velocidades e

densidade apresentados na Tabela 2 (valores citados em (CASTAGNA; SWAN, 1997)). Na

Figura 10 observa-se uma esquematização do modelo de reservatório estudado.

Parâmetro Valor

VP 2.9km/s

VS 1.33km/s

⇢ 2.29g/cm3

Tabela 2 - Valores do folhelho selante: velocidades das ondas P e S e densidade.

ArenitoReservatório

FolhelhoSelante

Figura 10 - Esquema do modelo de reservatório.

5.2 GERAÇÃO DE DADOS

Utilizando os parâmetros do modelo definidos nas Tabelas 1 e 2 e a teoria apresentada nos

Capítulos 2 e 3, pode-se modelar os dados sintéticos que serão invertidos.

Inicialmente foram utilizadas as equações 7, 9 e 11 para o cálculo do módulo de bulk da

rocha saturada Ksat e as equações 12 e 13 para o cálculo da densidade média no reservatório

saturado, ambas como funções da saturação de óleo S e da porosidade � do reservatório.

Com os módulos elásticos e com a densidade como função de S e �, utilizaram-se as

equações 1 e 2 para o cálculo das velocidades. Tem-se dessa maneira uma modelagem das

velocidades e da densidade do reservatório em função de sua saturação de óleo e de sua

porosidade:

VP = f(S,�) (23)

VS = f(S,�) (24)

34

⇢ = f(S,�) (25)

Nas Figuras 11, 12 e 13 pode-se observar o comportamento dos parâmetros modelados

(equações 23, 24 e 25) com a variação da saturação de óleo S para alguns valores de porosidade

do reservatório.

A partir de 23, 24 e 25 calculou-se os valores dos atributos sísmicos modelados como função

da saturação S e da pororosidade �, empregando-se respectivamente, as equações 16, 19 e

20, obtendo-se:

RPP = f(S,�) (26)

t0 = f(S,�) (27)

xc = f(S,�) (28)

RPP , t0 e xc são portanto os dados sintéticos que podem ser gerados para qualquer com-

binação de S e �.

Foram ainda incluidos ruídos no dados, representados pelas equações 26, 27 e 28, para

simular uma situação de atributos observados da sísmica adquirida. Para a inclusão de ruídos

utilizou-se a função rand() do Matlab para construir um dado ruidoso yr a partir de um dado

y utilizando a seguinte expressão:

yr = y + ⌘ = y + (f ⇤ y) ⇤ (rand()� 0.5) ⇤ 2) (29)

onde f é a fração que mede quanto os dados ruidosos devem se dispersar ao redor do dado

y e rand() é uma função aleatória que gera números de 0 a 1.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 13.2

3.4

3.6

3.8

4

4.2

4.4

4.6

4.8

5

5.2

Saturação

Velo

cidade d

a o

nda P

Velocidade da onda P para diferentes saturações

5%10%20%30%

Figura 11 - Velocidade da onda P em função da saturação para diferentes �’s.

35

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 11.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6

Saturação

Velo

cidade d

a o

nda S

Velocidade da onda S para diferentes saturações

5%10%20%30%

Figura 12 - Velocidade da onda S em função da saturação para diferentes �’s.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 12

2.1

2.2

2.3

2.4

2.5

2.6

2.7

Saturação

Densi

dade d

a c

am

ada r

ese

rvató

rio

Densidades para diferentes saturações

5%10%20%30%

Figura 13 - Densidade em função da saturação para diferentes �’s.

36

6 ANÁLISES E RESULTADOS

Inicialmente serão desenvolvidas análises referentes à performance do algoritmo de inversão

baseado no algoritmo genético; em seguida, com o intuito de avaliar a sensibilidade na recupe-

ração dos parâmetros de inversão, serão abordadas análises referentes à composição da Função

Objetivo e à faixa de afastamentos (ou ângulos de incidência) considerada. Presença de ruído

e influência de erros nos parâmetros assumidos como conhecidos também foram avaliados. E,

por fim, será testado o método desenvolvido para a recuperação da variação da saturação de

óleo no reservatório em dois cenários de porosidade com o intuito de verificar sua eficiência e

viabilidade para o monitoramento do reservatório.

6.1 PERFORMANCE DO ALGORITMO DE INVERSÃO

Devido à natureza estocástica do algoritmo genético, os valores de VP , VS e ⇢ retornados

pelo algoritmo como ponto de mínimo da Função Objetivo não serão os mesmos para cada

execução do algoritmo. Entretanto, espera-se uma precisão mínima apresentada pelo algoritmo

de otimização; para isso foi desenvolvido um Algoritmo de Inversão (I) que executa n vezes

o algoritmo genético para uma dada configuração de população e do número máximo de

gerações iteradas pelo GA.

Conforme os conceitos discutidos em 4.3, a população é o conjunto de modelos (indivíduos)

disponíveis em cada iteração (geração) e o limite de gerações é um dos critérios de parada das

iterações dentro de uma execução do GA. A população inicial de modelos é gerada de maneira

aleatória dentro de um determinado espaço inicial de busca. Esse espaço inicial de busca dos

parâmetros foi calibrado de acordo com o esperado dos modelos apresentados nas Figuras 11,

12 e 13. É importante apontar que mesmo em uma situação real, esse tipo de informação para

calibrar o espaço de busca está disponível, principalmente no caso do 4D aonde testemunhos,

informações de poços e modelos de velocidade estão disponíveis.

Para a análise da performance do algoritmo de inversão, foram analisados 5 cenários de

configurações iniciais do Algoritmo de Inversão para dois valores de porosidade do reservatório:

5% e 20%.

As análises realizadas consistiram em efetuar a inversão de dados gerados pela modelagem

direta a partir de valores de p(C) =⇣V (C)P , V

(C)S , ⇢

(C)⌘

previamente estabelecidos, denominados

de parâmetros controle, e em seguida, compará-los com os parâmetros recuperados p(R) =⇣V (R)P , V

(R)S , ⇢

(R)⌘

por essa inversão. Para a porosidade � = 5% temos os parâmetros de

controle na Tabela 3a e os dados gerados na Figura 14a; para a porosidade de � = 20%

os parâmetros de controle estão na Tabela 3b e os dados na Figura 14b. Note, que em um

37

primeiro momento, não foram incluídos ruídos nos dados para verificar a eficiência intrínseca

do algoritmo em um caso ideal.

0 10 20 30 40 50 60 70 80 90−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

θ

RP

P

Dados Controle

xc=1.1568(km)

t0=0.3970(s)

(a) � = 5%.

0 10 20 30 40 50 60 70 80 90−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

θ

RP

P

Dados Controle

xc=1.7866

t0=0.4050(s)

(b) � = 20%.

Figura 14 - Dados gerados - Coeficientes de Reflexão (RPP ) em função do ângulo de incidência

(✓), distância crítica de refração (xc) e tempo normal de reflexão (t0) - para os parâmetros de

controle (Tabela 3) nos dois cenários de porosidade.

38

p(C) valor

VP 3.8333

VS 1.2497

⇢ 2.5614

(a) � = 5%.

p(C) valor

VP 3.3234

VS 1.3201

⇢ 2.2955

(b) � = 20%.

Tabela 3 - Parâmetros de controle para os dois cenários de porosidade.

As inversões dos dados sintéticos foram então efetuadas utilizando o algoritmo de inver-

são I com 100 execuções do algoritmo genético n = 100 para cada um dos 5 cenários de

configurações iniciais do Algoritmo de Inversão avaliados.

Nas Tabelas 4a e 4b encontram-se sumarizados os cenários (colunas 1 a 3) e os critérios

adotados para a análise (colunas 4 a 6).

Para verificar a performance do algoritmo I em cada um dos cenários, o critério adotado para

mensurar a eficiência foi o de comparar os parâmetros recuperados com os parâmetros controle

em cada cenário. Existem vários outros critérios para quantificar a performance do algoritmo;

aqui não foi avaliado o melhor deles e sim foi escolhido um critério para verificar a perfor-

mance relativa entre os casos estudados. Neste trabalho, a comparação foi feita utilizando-se

a distância vetorial " (denominada de erro) no espaço de parâmetros entre os parâmetros

recuperados p(R) =⇣V (R)P , V

(R)S , ⇢

(R)⌘

e os parâmetros controle p(C) =⇣V (C)P , V

(C)S , ⇢

(C)⌘:

" = p(R) � p(C) =r⇣

V (R)P � V(C)P

⌘2+⇣V (R)S � V

(C)S

⌘2+ (⇢(R) � ⇢(C))2 (30)

Foram computados para cada cenário, como variavam, em função de n:

• o erro entre a média dos parâmetros recuperados até o passo n (⇣p(R)n

⌘

médio) e os

parâmetros controle (p(C)), chamado de "médio;

• o erro entre os parâmetros recuperados no passo n (⇣p(R)n

⌘

passo) e os parâmetros controle

(p(C)), chamado de "passo;

Em seguida, foram tabelados os valores: do erro médio para o último passo ("médio)final =

("médio)n=100; o menor valor dentre todos os n valores do erro de cada passo ("passo)min; e o

tempo de execução do algoritmo para recuperar os parâmetros em cada cenário.

39

Cená

rioTa

man

hoda

Popu

laçã

oLi

mite

deG

eraç

ões

("méd

io) f

inal

("passo) m

inTe

mpo

deex

ecuç

ão(m

in)

150

300

8.11

28e-

061.

7653

e-06

20

215

030

06.

4516

e-05

2.30

11e-

0637

350

509.

5949

e-04

1.39

03e-

048

415

050

3.33

59e-

042.

6431

e-05

15

525

030

01.

1128

e-05

8.76

53e-

0752

(a) � = 5%.

Cená

rioTa

man

hoda

Popu

laçã

oLi

mite

deG

eraç

ões

("méd

io) f

inal

("passo) m

inTe

mpo

deex

ecuç

ão(m

in)

150

300

4.95

45e-

062.

8070

e-06

28

215

030

05.

3395

e-05

1.77

02e-

0544

350

509.

4492

e-04

4.10

96e-

049

415

050

7.77

40e-

043.

0079

e-04

14

525

030

08.

8269

e-06

1.73

99e-

0758

(b) � = 20%.

Tabela 4 - Cenários de execução do algoritmo de inversão (I) e resultados dos critérios adotados

para a análise de performance da inversão. (ver texto)

40

Para cada um dos cenários, foram efetuadas análises dos resultados conforme gráficos apre-

sentados nas Figuras 15 e 16: à esquerda, os dois gráficos relacionam-se ao comportamento

de "médio (superior) e "passo (inferior) em função do número de passos (n) do algoritmo I;

e os dois gráficos à direita ilustram o espaço dos 3 parâmetros da inversão (VP ,VS,⇢); no

superior estão representados os n = 100 valores recuperados para cada iteração do algoritmo

I (⇣p(R)n

⌘

passo) em escala de cores segundo seu respectivo valor de função objetivo, nele tam-

bém está representado por uma estrela preta o ponto correspondente aos parâmetros controle

(p(C)); e no inferior representa-se a média de todas as n = 100 iterações de I por um ponto

vermelho e o conjunto de parâmetros recuperados no passo correspondente ao menor erro

(("passo)min) por um ponto verde.

Observando-se, como exemplo os resultados obtidos para o cenário 1 nas Figuras 15 e

16, pode-se sumarizar o comportamento do algoritmo de inversão (que se repete para os 10

cenários analisados) como se segue:

• Nos gráficos da variação de "médio em função do passo (iteração de I) pode-se observarum corportamento instável para valores de n até aproximadamente n s 50; tal fato édevido à aleatoriedade no comportamento de cada execução do GA dentro do algotimo

de inversão I. Quanto maior o valor de n, maior é o espaço amostral aonde realiza-se a

média; portanto para n s 50 já tem-se um espaço amostral suficientemente grande paraacomodar as variações de cada iteração de I;

• Nos gráficos da variação de "passo para cada passo (iteração de I) pode-se observar aaleatoriedade comentada acima na distribuição dos valores do erro para cada passo n;

• Nas distribuições das recuperações nos espaços de parâmetros, nota-se que quanto maior ovalor da função objetivo mais distante está o ponto representando a recuperação (iteração)

dos parâmetros controle (p(C)); também nota-se que a distância entre o ponto controle

(estrela preta) e o ponto correspondente a média de todas as iterações (ponto vermelho

no gráfico inferior) é sempre maior que a distância entre o mesmo ponto controle (estrela

preta) e o ponto correspondente à recuperação com menor erro (ponto verde no gráfico

inferior); além disso, o ponto verde corresponde ao menor valor de FOBJ .

Pode-se inferir, então, que o ponto com menor valor de função objetivo é aquele que mais

se aproxima do ponto de controle. Observando-se também os parâmetros de configuração

utilizados no GA, os tempos de execução, e os erros de cada um dos 10 cenários nas Tabelas

4a e 4b, nota-se que quanto maior a população inicial, melhor é o resultado da inversão,

entretanto, o tempo de execução do algoritmo também aumenta consideravelmente.

41

010

20

30

40

50

60

70

80

90

100

2

2.53

3.54

4.55

5.56

6.5

x 10

−5

n

εmédio

ε e

m funçã

o n

o n

úm

ero

de p

ass

os

010

20

30

40

50

60

70

80

90

100

0123x

10

−4

n

εpasso

ε p

ara

cada p

ass

o

3.8

333

3.8

333

3.8

334

3.8

334

3.8

334

1.2

496

1.2

497

1.2

498

1.2

499

2.5

613

2.5

614

2.5

614

2.5

615

2.5

615

2.5

616

VP

Esp

aço

de P

arâ

metr

os

VS

ρ

min

max

p(C

)

itera

ções

3.8

333

3.8

333

3.8

334

3.8

334

3.8

334

1.2

496

1.2

497

1.2

498

1.2

499

1.2

52.5

612

2.5

613

2.5

614

2.5

615

2.5

616

VP

Esp

aço

de P

arâ

metr

os

VS

ρ

p(C

)

média

menor

εpass

o

Figura 15 - Análises dos resultados da inversão - � = 5% (ver texto).

42

010

20

30

40

50

60

70

80

90

100

2

2.53

3.54

4.55

x 10

−5

n

εmédio

ε e

m funçã

o n

o n

úm

ero

de p

ass

os

010

20

30

40

50

60

70

80

90

100

0123x

10

−4

n

εpasso

ε p

ara

cada p

ass

o

3.3

233

3.3

234

3.3

234

3.3

234

3.3

234

1.3

2

1.3

201

1.3

202

1.3

203

2.2

955

2.2

956

2.2

957

2.2

958

VP

Esp

aço

de P

arâ

metr

os

VS

ρ

min

max

p(C

)

itera

ções

3.3

233

3.3

233

3.3

234

3.3

234

3.3

234

1.3

2

1.3

201

1.3

202

1.3

203

1.3

204

2.2

955

2.2

956

2.2

957

2.2

958

VP

Esp

aço

de P

arâ

metr

os

VS

ρ

p(C

)

média

menor

εpass

o

Figura 16 - Análises dos resultados da inversão - � = 20% (ver texto).

43

Efeito de ruído nos dados

Os testes simulados foram executados a partir dos dados sintéticos sem inclusão alguma de

ruído, conforme pode-se observar nas Figuras 14a e 14b. Utilizando os mesmos parâmetros

controle nas Tabelas 3a e 3b foram gerados dados com a inclusão de ruídos corforme a equação

29, com f = 0.07, ou seja, 7% de ruído. Os novos dados estão nas Figuras 17a e 17b.

0 10 20 30 40 50 60 70 80 90−1.5

−1

−0.5

0

0.5

1

θ

RP

P

Dados Controle

dado sem ruídodado com ruído

xc=1.1067(km)

t0=0.3897(s)

(a) � = 5%.

0 10 20 30 40 50 60 70 80 90−1.5

−1

−0.5

0

0.5

1

θ

RP

P

Dados Controle

dado sem ruídodado com ruído

xc=1.8794(km)

t0=0.3965(s)

(b) � = 20%.

Figura 17 - Dados da Figura 14 com adição de ruído.

A partir dos novos dados sintéticos com ruído, foram repetidos os testes para o melhor

e o pior cenário de acordo com o critério do menor erro encontrado (("passo)min) em cada

um dos modelos de porosidade. Nos dois modelos o melhor cenário foi o 5 e o pior o 3

conforme as Tabelas 4a e 4b. Observa-se nos resultados apresentados na Tabela 5 que os

ruídos introduzidos nos dados acarretam em erros na solução da inversão maiores do que a

44

precisão do algoritmo de inversão (valores apresentados na Tabela 4). A precisão do algoritmo

foi definida como sendo o erro devido apenas à performance do próprio algoritmo, que também

será chamado de erro intrínseco do algoritmo.

Cenário ("médio)final ("passo)min Tempo de execução (min)

5 - Melhor 0.003924801320626 0.003925734302610 42

3 - Pior 0.004032656445837 0.004042420208598 9

(a) Melhor e pior cenário para � = 5%.

Cenário ("médio)final ("passo)min Tempo de execução (min)

5 - Melhor 0.002380838331268 0.002388961885877 35

3 - Pior 0.002362438330779 0.002371212699500 7

(b) Melhor e pior cenário para � = 20%.

Tabela 5 - Resultados das análises de performance para as inversões a partir de dados ruidosos.

Os testes realizados para inversão dos dados com ruído apresentam diferença nas soluções

de cada passo do algoritmo I, com variação na casa decimal da precisão do algoritmo, o que

pode ser verificado comparando-se os valores apresentados nas Tabelas 4 e 5. Essa observação

corrobora a definição acima do critério de erro intrínseco do algoritmo.

A utilização do algoritmo de inversão diversas vezes para os dados com ruído geram dife-

rentes valores de erro médio e erro mínimo, mas sempre na mesma ordem de grandeza dos

valores apresentados na Tabela 5, exceto em uma das execuções para o modelo de porosidade

20% com o pior cenário (3), que chegou a apresentar erros na segunda casa decimal, com

variação das soluções na quarta casa decimal.

Na Figura 18 observa-se o espaço de parâmetros de inversão com as soluções obtidas

empregando-se o cenário 3 para o modelo de porosidade � = 5%. Nota-se que o conjunto de

soluções encontra-se agrupado em torno de um ponto que está deslocado do ponto correto

da ordem de 10�3. Embora o ruído adicionado aos dados tenha incluído uma tendência

nas soluções da inversão, o conjunto das soluções está dentro da precisão aceitável para os

parâmetros buscados.

45

4.2124.2125

4.2134.2135

4.2144.2145

4.215

2.61

2.612

2.614

2.616

2.618

2.622.255

2.26

2.265

2.27

2.275

2.28

VP

Espaço de Parâmetros

VS

ρ

min

max

p(C)

p(R)

Figura 18 - Soluções da inversão para os dados sintéticos com ruído.

Critério Final

Sumarizando as observações dos testes e análises da performance do algoritmo:

1. Utilizando valores de n s 50 (número de passos do algoritmo I) já tem-se um espaçoamostral grande o suficiente para o comportamento do erro médio "médio estabilizar-se e

para buscar-se um "passo mínimo com uma precisão melhor do que a do ("médio)final;

2. Utilizando-se o pior caso de performance do algoritmo com o critério de erro ("passo)min,

obtem-se um erro da ordem de 10�5 para dados sem ruído;

3. Ruídos introduzidos nos dados conduzem a erros da ordem de 10�2 a 10�3 ;

4. O ponto de parâmetros recuperados pelo algoritmo (p(R)) correspondente à iteração de

I que apresentou o menor ("passo)min é o que apresenta menor valor de função objetivo;

Considerando os resultados acima a configuração, o critério escolhido a partir de agora será:

• O critério para a escolha do ponto recuperado p(R) aceito será o que apresenta menorvalor de FOBJ ; e o erro chamado de ✏ a partir de agora, será o relativo a este ponto;

• O cenário de configuração utilizado para calibrar o algoritmo I será o cenário 4 comn s 50 (tendo seu valor variado caso necessário) pois apresenta a precisão necessáriapara o estudo de dados ruidosos e menor tempo computacional para sua execução.

46

6.2 ANÁLISE DA COMPOSIÇÃO DAS FUNÇÕES OBJETIVO

Cada atributo sísmico definido no Capítulo 3, conduz a uma respectiva função objetivo:

uma para as amplitudes, uma para o tempo de chegada e uma para a distância crítica. Para

construír-se cada uma das funções objetivo utiliza-se a definição dada pela equação 22.

Utilizando-se a equação 16 para calcular a amplitude A de uma onda refletida com o ângulo

k, tem-se que a função objetivo das amplitudes FAOBJ é:

FAOBJ =

vuut✓finalX

k=✓0

(A(k)obs � A(k)calc)2 (31)

onde A(k)obs é o dado de amplitude observado no ângulo k e A(k)calc é a amplitude

calculada para os modelos testados pelo algoritmo de inversão nesse mesmo ângulo. A soma

ocorre com os valores de k variando do ângulo inicial de aquisição dos dados (✓0) até o ângulo

final de aquisição (✓final).

Utilizando a equação 19 para o cálculo do tempo de chegada normal t0, tem-se a função

objetivo dos tempos de chegada normal, para um sismograma CMP dada simplesmente por:

F t0OBJ =��tobs0 � tcalc0

�� (32)

onde tobs0 é o dado observado e tcalc0 é o respectivo valor iterado pelo algoritmo.

Analogamente para a componente referente à distância crítica (equação 20) tem-se:

F xcOBJ =��xobsc � xcalcc

�� (33)

onde xobsc é o dado observado da distância crítica e xcalcc é o respectivo valor calculado para

os modelos durante a inversão.

Serão testadas agora diferentes formas de associação e combinação das funções objetivos

com o intuito de analisar a influência dos respectivos atributos no processo de inversão. As

funções objetivo 31, 32 e 33 serão combinadas de 4 maneiras diferentes:

• Somatória simples

O primeiro modo de construção é uma função objetivo construída como sendo a soma das

funções objetivos individuais provenientes dos 3 atributos sísmicos:

F 1OBJ = FAOBJ + F

t0OBJ + F

xcOBJ

F 1OBJ =

vuut✓finalX

k=✓0

(A(k)obs � A(k)calc)2 +��tobs0 � tcalc0

��+��xobsc � xcalcc

�� (34)

47

• Mínimos Quadrados do conjunto total de atributos

Para a construção da função objetivo desta segunda maneira, os dados de distância crítica e

de tempo de reflexão normal foram tratados como sendo apenas mais uma informação dentro

de um vetor de dados já contendo os valores das amplitudes. Em outras palavras, foi criado

um vetor de dados d = (A(✓), t0, xc) que foi minimizado a partir da diferença dos mínimos

quadrados:

F 2OBJ = FA,t0,xcOBJ =

q(dobs � dcalc)2

F 2OBJ = FA,t0,xcOBJ =

vuut✓finalX

k=✓0

(A(k)obs � A(k)calc)2 +�tobs0 � tcalc0

�2+ (xobsc � xcalcc )

2 (35)

• Soma Mínimos Quadrados das Amplitudes e do conjunto parcial do tempo de chegada edistância crítica

Nesta terceira forma de construir a função objetivo utilizou-se a soma da função objetivo das

amplitudes (equação 31) com a função objetivo correspondente aos mínimos quadrados do

conjunto parcial do tempo de chegada com a distância crítica:

F 3OBJ = FAOBJ + F

t0,xcOBJ

F 3OBJ =

vuut✓finalX

k=✓0

(A(k)obs � A(k)calc)2 +q�

tobs0 � tcalc0�2

+ (xobsc � xcalcc )2 (36)

• Inversão somente com as Amplitudes

Foi também executado um exemplo de inversão com a função objetivo somente com a con-

tribuição das amplitudes.

F 4OBJ = FAobj

F 4OBJ =

vuut✓finalX

k=✓0

(A(k)obs � A(k)calc)2 (37)

As inversões foram executadas novamente utilizando os mesmos dados e parâmetros con-

trole das Tabelas 3a e 3b e das Figuras 14a e 14b para os casos de � = 5% e � = 20%

respectivamente.

48

FOBJ ✏

F 1obj 4.4367e-06

F 2obj 6.3564e-06

F 3obj 4.8745e-06

F 4obj 1.9545e-06

(a) � = 5%.

FOBJ ✏

F 1obj 6.7337e-06

F 2obj 2.8270e-06

F 3obj 1.9280e-05

F 4obj 8.8707e-06

(b) � = 20%.

Tabela 6 - Análises da composição da FOBJ a partir de dados ideais: erro mínimo encontrado

durante a execução do Algoritmo de Inversão (I).

Pode-se observar nos resultados da Tabela 6 que a precisão no caso ideal (sem ruído nos

dados) de todas as FOBJ é da ordem de 10�6 a 10�5 conforme já havia se observado na seção

de análises de performance. na qual a FOBJ utilizada nos testes foi a F 3OBJ .

Efeito de ruído nos dados

As inversões realizadas até aqui para as FOBJ utilizaram dados sem ruído. Os erros das

soluções obtidas da inversão utilizando-se os dados ruidosos das Figuras 17a e 17b para a

execução de I são apresentados na Tabela 7.

Note que o erro ✏, utilizado como critério de escolha das melhores FOBJ , é composto pelas

componentes de erro na recuperação de VP , VS e ⇢, conforme a equação 30. É interessante

quantificar a contribuição de ✏VP , ✏VS e ✏⇢ na composição do erro total ✏. Pode-se explicitar

as componentes que contribuem para a composição do erro reescrevendo a equação 30 como:

✏ = p(R) � p(C) =r⇣

V (R)P � V(C)P

⌘2+⇣V (R)S � V

(C)S

⌘2+ (⇢(R) � ⇢(C))2

✏ = p(R) � p(C) =q(✏VP )

2 + (✏VS)2 + (✏⇢)

2 (38)

Na Tabela 7, também estão quantificadas as contribuições de cada componente para o erro

total.

Observando-se os resultados presentes nas Tabelas 7a e 7b pode-se notar claramente uma

equivalência entre o par de funções objetivo F 1OBJ e F 3OBJ e entre o par F 2OBJ e F 4OBJ .

Levando-se em conta as equações 34, 35, 36 e 37 conclui-se que essa equivalênca ocorre

devido ao maior peso dado aos atributos xc e t0 nas F 1OBJ e F 3OBJ e o menor peso dado a eles

nas F 2OBJ e F 4OBJ .

O par que conduz aos menores erros, ao contrário do esperado, é o par F 2OBJ e F 4OBJ ,

aquele que dá menos importância ao xc e ao t0 . O esperado seria que, ao introduzir-se

49

FOBJ ✏ ✏VP ✏VS ✏⇢

F 1obj 0.0082 2.7157e-04 0.0035 0.0074

F 2obj 0.0081 3.2434e-04 0.0028 0.0076

F 3obj 0.0084 3.0104e-04 0.0038 0.0075

F 4obj 0.0078 3.0159e-05 0.0034 0.0070

(a) � = 5%.

FOBJ ✏ ✏VP ✏VS ✏⇢

F 1obj 0.0038 9.2302e-04 0.0020 0.0031

F 2obj 0.0036 4.4834e-04 0.0023 0.0027

F 3obj 0.0039 1.0687e-04 0.0025 0.0030

F 4obj 0.0037 9.8623e-04 0.0019 0.0030

(b) � = 20%.

Tabela 7 - Contribuição de cada componente do erro na composição total de ✏ para cada

FOBJ a partir de dados ruidosos.

atributos mais sensíveis na função objetivo, principalmente o

Documents

INVERSÃO DE ATRIBUTOS SÍSMICOS PARA CARACTERIZAÇÃO … · Figura 19 Faixas de ângulos consideradas para a inversão no modelo =5%... 51 Figura 20 Faixas de ângulos consideradas