Embed Size (px)
Citation preview
Apostila de Tomografia Sısmica
Professor Luiz Alberto Santos1
1Universidade Federal Fluminense Email: [email protected]
1
INDICE
1- Objetivos
2- Introducao
3- Velocidades sısmicas
4- Metodos de estimativa de velocidades do meio geologico
5- Modelagem sısmica por tracado de raios
6- Fundamentos de tomografia sısmica
7- Formas de discretizacao do meio fısico para tomografia
8- A Obtencao da matriz de sensibilidade (L)
9- Inversao tomografica com fidelidade geologica
10- Acesso aos exercıcios praticos via WEB.
11- Referencias Bibliograficas
2
1- OBJETIVOS E MOTIVACAO
Primeiramente proporcionar entendimento dos princıpios e tecnicas envolvidos na inversao
tomografica de velocidades sısmicas. Em seguida, resolver o problema tomografico esti-
mando velocidades crıveis para o meio geologico.
3
2- INTRODUCAO
Tomografia significa descricao (grafia) das partes (tomo (do grego)). Essa descricao pode se
referir a qualquer propriedade: velocidade P (Vp), velocidade S (Vs), densidade ou massa
especıfica (ρ), susceptibilidade magntica, concentracao de um elemento quımico (uranio
XU, potassio XK, torio XTh, ouro XAu, cobre XCu, ferro XFe etc..), idade, porosidade (φ),
permeabilidade (k), litologia ou tipo de rocha (ıgnea, metamorfica ou sedimentar), mergulho
(Dip) e Azimute de acamamento sedimentar (S0), densidade de fraturas etc.
A propriedade a ser descrita neste trabalho e a velocidade sısmica da onda P ou com-
pressional. Como se vera, os metodos matematicos empregados no processo de inversao
tomografica sao amplos e genericos. Aplicam-se a qualquer area da ciencia. Ao longo de
todo texto, contudo, o metodo tomografico e particularizado para o problema da estimativa
de velocidade P no meio geologico. Por conta destas especificidades, ha um capıtulo ded-
icado a definicao das diferentes velocidades sısmicas P e o seu comportamento em funcao
do ambiente fısico.
Como sera visto mais a frente, o problema inverso nao linear necessita de uma ferramenta
de resolucao do problema direto (ao longo de todo o texto os termos simulacao, modelagem
e problema direto estarao sendo empregados como sinonimos). O problema a ser resolvido
na modelagem consiste simplesmente em calcular o tempo que uma onda P se propaga de
um ponto a outro. A ferramenta a ser empregada e o tracado de raios. De forma generica
o problema direto pode ser definido com a igualdade d = Lm. Onde d constitui os dados,
L um operador e m o modelo. Particularizando para o tema em escopo, d e o tempo de
transito, m e a velocidade P e L o operador que representa o tracado de raios.
O problema inverso pode tambem ser definido genericamente pela igualdade m = L−1d,
4
ou seja, a estimativa do modelo a partir dos dados. Particularizando, estimar a velocidade a
partir dos temps de transito. Somente na metade seguinte desta obra, foca-se no problema
tomografico inverso para a onda P. Primeiramente com a fundamentacao matematica, as
formas de representacao do campo de velocidades, as formas correspondentes de calculo
da matriz de sensibilidade e sugestoes para estimativa de velocidades crıveis para o mdeio
geologico.
5
3- VELOCIDADES SISMICAS
Neste topico serao definidos os diferentes termos relacionando velocidades sısmicas pre-
sentes no dia a dia da industria do petroleo. Serao aqui abordadas seis velocidades cujas
definicoes sao encontradas em livros texto de geofısica, entre eles Sheriff e Geldart (1999):
velocidade de fase; velocidade de grupo; velocidade media; velocidade RMS; velocidade de
empilhamento ou NMO; velocidade intervalar.
Ao se admitir que um meio e dispersivo, remete-se ao fato de que a velocidade depende da
frequencia. Entao, cada um dos harmonicos, ou fases, que compoe um pulso sısmico viaja
com uma velocidade distinta, denominada velocidade de fase. A envoltoria dos diversos
harmonicos viaja com uma velocidade que e chamada de velocidade de grupo vg. Neste
trabalho assume-se que o meio geologico nao e dispersivo, logo, todas as fases e a envoltoria
(ou grupo) possuem a mesma velocidade.
A velocidade media (vm) e um conceito emprestado da mecanica e e definida como a
razao entre o distancia (modulo da posicao final menos posicao inicial) e o tempo decorrido
para um pulso viajar do ponto inicial ao final. Nao importa a trajetoria.
A velocidade RMS (vrms) e alcancada pela media quadratica (Sheriff and Geldart, 1999),
definida pela expressao:
vrms,n =
√
∑ni=1 v
2i ti
∑ni=1 ti
, (1)
onde vi e a velocidade intervalar ao longo do intervalo ti.
A velocidade de empilhamento ou NMO (Normal Move Out) (vnmo) e obtida durante
o processamento sısmico e e aquela capaz de horizontalizar eventos hiperbolicos nas secoes
6
de afastamento medio comum (Common Mid Point Gather - CMP). Em geral a velocidade
NMO se aproxima da velocidade RMS, sendo esta ultima pouco menor que a primeira. Para
um modelo de camadas plano-paralelas, vnmo > vrms > vm.
VARIAVEIS QUE CONTROLAM AS VELOCIDADES SISMICAS
Este trabalho dedica-se a estimativa de velocidades sısmicas a partir de inversao to-
mografica. Assim, e desejavel que sejam conhecidas a natureza e as caracterısticas gerais
do modelo a fim de se reduzir a multiplicidade de solucoes no Problema Inverso. O conhec-
imento a priori do modelo m norteia a melhor forma de representa-lo, seja na insercao de
vınculos (formas de restricao do espaco de modelos inseridas nas equacoes de inversao) nas
equacoes de inversao, seja na regularizacao (processo voltado para estabilizacao da inversao
por meio da imposicao de vınculos e restricoes que induzem a solucao).
Rocha e qualquer massa agregada firme e coerente de minerais que constitui parte de um
planeta (Skinner e Porter, 1987). Mineral constitui qualquer substancia com composicao
quımica e estrutura cristalina definida, formada por processos naturais.
Alem dos minerais constituintes, as rochas contem poros, que nada mais sao que o
espaco vazio entre os graos minerais (ou polimineralicos). Dada uma amostra de rocha, sua
porosidade e definida como a razao entre o volume de vazios e o seu volume total (Mavko
et al, 2009). No meio geologico as rochas sedimentares possuem seus poros preenchidos por
fluidos (e.g. ar, agua (doce ou salgada), gas, condensado e/ou oleo).
Os fatores intrınsecos que ditam as propriedades fısicas das rochas sao a composicao
mineral, o contato entre os graos, a forma dos graos, a porosidade e os fluidos interstici-
ais (composicao e pressao de poros). Externamente, temperatura e pressao confinante ou
litostatica influenciam fortemente a velocidade sısmica nos litotipos (Mavko et al., 2009 e
7
Thomas, 2000). A primeira tende a reduzir a velocidade e os modulos de rigidez a medida
que a temperatura aumenta. A pressao confinante crescente age no sentido de aumentar a
velocidade pela reducao de porosidade. A Figura 1 exibe a curva de variacao de velocidade
de onda P com a pressao confinante. Os dados extraıdos de Carmichael (1982) referem-se
a um arenito saturado com agua e sob pressao de poro igual a zero.
Figure 1: Curva de variacao da velocidade P em funcao da pressao de um arenito
extraıda de Carmichael (1982).
8
Em funcao das variaveis citadas no paragrafo anterior, uma mesma litologia experimenta
larga faixa de variacao de velocidades P. Em Mavko et al. (2009) e Carmichael (1982) ha
extensas listas contendo as velocidades de onda P para diversos litotipos sob diferentes
condicoes de temperatura, pressao confinante, tipos de fluidos e saturacao.
As rochas sedimentares formam-se a partir de sedimentos e estes sao depositados sob a
forma de camadas. Uma camada constitui um intervalo de rocha limitado superiormente
e inferiormente por superfıcies (interfaces), estas normalmente representando um intervalo
de nao deposicao ou erosao (Park, 1997). Sheriff (2002) define camada como um intervalo
cujas propriedades diferem do que esta acima e do que esta abaixo. A definicao de Park
(1997) tem carater genetico, informa a causa da distincao entre as camadas, enquanto a
definicao de Sheriff (2002) expressa a consequencia, que e a propria diferenca de propriedade
entre estratos. As ondas sısmicas sao sensıveis aos contrastes de propriedades mecanicas
denunciados pelas reflexoes destas ao atingirem as interfaces. A atitude (relacao espacial
de uma feicao com a horizontal, como o mergulho de uma camada, Sheriff (2002)) original
das rochas sedimentares, ou acamamento sedimentar, e predominantemente horizontal a
subhorizontal.
Em funcao das condicoes de fluxo do meio no qual os sedimentos estavam imersos e
das demais condicoes ambientais, uma mesma camada pode apresentar variacoes laterais de
propriedades ou facies (caracterısticas que distinguem uma rocha de outra adjacente, seja
na litologia, estruturas sedimentares, granulometria, etc. (Sheriff, 2002)).
Ainda nao explorado nestas anotacoes, por ora, resume-se o efeito das diversas variaveis
que controlam a velocidade P nas rochas na figura a seguir.
9
Figure 2: Efeito da pressao, temperatura e litologia na velocidade.
10
Ao longo da evolucao geologica, as rochas sedimentares podem experimentar episodios
de deformacao de diversas origens (tectonica ou nao) e intrusoes ıgneas, halinas ou mesmo
de argila, que alteram a geometria e atitude dos estratos gerando dobras e falhas. A Figura
3 resume as estruturas basicas encontradas em bacias sedimentares. Qualquer estrutura
mais complexa pode ser obtida com a combinacao de um ou mais tipos encontrados na
Figura 3.
11
Figure 3: Relacao de estruturas basicas e propriedades (cores) observadas em rochas
sedimentares. (A) Acamamento sedimentar plano-paralelo; (B) Acamamento sedimen-
tar plano-paralelo com variacao lateral de faceis; (C) Dobras; (D) Falhamento normal;
(E) Falhamento reverso; (F) Intrusao discordante.
12
4- METODOS DE ESTIMATIVA DE VELOCIDADES NO MEIO
GEOLOGICO
4.1- Introducao
4.2- Sismos Naturais
4.3- Experimento sısmico
4.4- Exercıcios
INTRODUCAO
Existe uma mirıade de meios para estimar a velocidade no meio geologico.
ANALISE DE VELOCIDADE DE EMPILHAMENTO + DIX
Em dados sısmicos de reflexao, seja uma famılia CMP, representada fisicamente pela
Figura 4. Para um modelo com camadas plano paralelas horizontais e sem variacao lateral
de velocidade, a equacao que relaciona o tempo do evento com a posicao x e dada pela
equacao 2.
t(x) =
√
t20 +x2
v2nmo
, . (2)
Sob as condicoes elencadas no paragrafo anterior, um exemplo de CMP para um modelo
de quatro camadas, e apresentado na Figura 5. Depois de aplicar a velocidade NMO correta
para cada evento, equacao 2, obtem-se o CMP corrigido, Figura 6.
13
Figure 4: Representacao pictorica de um CMP.
14
Figure 5: Famılia CMP antes da correcao NMO.
15
Figure 6: Famılia CMP depois da correcao NMO.
16
Considerando-se a velocidade vnmo muito proxima da vrms, as velocidades intervalares
podem ser obtidas da vnmo, bastando para isso, substituir vrms por vnmo na equacao 1, Dix
(1955).
TOMOGRAFIA DE REFLEXAO
Em Pollet et al. (1983) sao reportadas as chamadas anomalias de velocidade de empil-
hamento. Fazendo referencia aos levantamentos marıtimos, Bishop et al.(1985) menciona
que as citadas anomalias se fazem presentes quando ha variacoes de velocidade com com-
primento de onda inferior a dimensao dos cabos de aquisicao marıtima (streamers). Em
funcao de maiores complexidades no campo de velocidade, nasceu a busca por tecnicas
tomograficas voltadas para determinacao desse campo.
Destaca-se a publicacao pioneira de Bishop et al. (1985) versando sobre tomografia de
reflexao. Neste trabalho o dado observado compreende os tempos de transito das reflexoes
mapeadas em domınios nao empilhados. A modelagem e efetuada por tracado de raios e
a diferenca entre os tempos calculados e os tempos observados constitui a funcao objetivo,
alvo da minimizacao. O algoritmo empregado para inversao e o de Gauss-Newton, um
metodo de gradiente, no qual a cada iteracao o resıduo e reduzido a medida em que o
campo de velocidade e atualizado juntamente com a superfıcie dos refletores.
Sword (1987), que baseou seu trabalho em Rieber (1936) e Riabinkin (1957), emprega
a informacao contida no parametro p (razao entre o seno do angulo de incidencia e a
velocidade intervalar) da fonte e dos receptores, alem do tempo de transito, como dados
para inversao tomografica. Esse trabalho contem as formulacoes que alicercam a tecnica
chamada estereotomografia, citada adiante nesta secao.
Woodward et al. (2008) fazem uma retrospectiva dos metodos tomograficos a partir da
17
decada de 1990. Os autores citam que neste perıodo ha avancos na resolucao dos modelos,
tais como a incorporacao dos efeitos anisotropicos devido a utilizacao de maiores lancos
e abertura azimutal e ainda, ao maior emprego de inversoes mais guiadas pelo dado que
inversoes governadas por interpretacoes e modelos. Os autores informam tambem que todos
esses avancos somente foram possıveis com o aumento do poder computacional, tanto em
processamento como em capacidade de armazenamento (memoria) de dados.
De fato, os anos 1990 experimentam grandes avancos nas tecnicas de inversao to-
mografica, motivados pelo crescente emprego da migracao em profundidade pre-empilhamento
(PSDM) e pelo aumento dos recursos computacionais. Observa-se nesse perıodo a seg-
mentacao da abordagem tomografica em duas vertentes principais: tomografia guiada pre-
dominantemente pelos dados, com poucas informacoes advindas de interpretacoes e; tomo-
grafia guiada por modelos, com fortes vınculos provenientes de interpretacoes e focada na
determinacao de gradientes de velocidade entre camadas. Esta ultima abordagem e aplicada
nas inversoes camada a camada (layer-stripping), na qual os parametros de velocidade sao
determinados separadamente, da camada mais rasa para a mais profunda. Observou-se que
os resıduos acumulados nas camadas mais rasas prejudicam a determinacao de parametros
nas camadas mais profundas (Soares Filho, 1994).
Billete e Lambare (1998) alem de publicar pioneiramente a estereotomografia colige
um rico historico sobre os metodos tomograficos, utilizado para confeccao deste capıtulo.
A tecnica estereotomografica adicionou mais parametros ao universo de dados. Assim,
constituem dados passıveis de minimizacao: o tempo de transito, posicao da fonte, posicao
dos receptores, parametro p (razao entre o seno do angulo de incidencia e a velocidade) da
fonte, parametro p do receptor. A modelagem a cada passo e realizada por meio do tracado
de raios paraxiais e os parametros a serem invertidos sao a posicao dos refletores e o campo
18
de velocidades.
ANALISE DE VELOCIDADE DE POR MIGRACAO
Em funcao dos avancos nos recursos computacionais a partir da decada de 1990, parale-
lamente aos avancos da tomografia registram-se numerosos trabalhos voltados para analise
de velocidade para migracao (migration velocity analysis - MVA). Na sısmica de reflexao,
um mesmo evento (ponto) pode ser amostrado por diferentes famılias de tiro. Depois da
migracao, caso o campo de velocidades esteja correto, deve se situar na mesma posicao. A
diferenca na posicao dos eventos pode ser mensurada comparando-se secoes de afastamento
comum migradas (common offset gathers - COG), Trier (1990) e Deregowski (1990), ou
implicitamente avaliando-se o alinhamento de CRPs (common reflection point), Al-Yahya
(1989), Symes e Carazzone (1991) e Jin e Madariaga (1994) e constitui a funcao objetivo
a ser minimizada. Seja qual for o domınio, almeja-se, implıcita ou explicitamente, alin-
har eventos nos paineis CRP. Stork (1992), resume o processo de estimativa do campo de
velocidades em domınios pos-migrados, similares a MVA.
TOMOGRAFIA EM POCOS
Adicionam-se os trabalhos relacionados a tomografia em dispositivos de aquisicao poco
a poco e Vertical Seismic Profiling (VSP), entre eles McMechan (1983), Soares Filho et
al. (1997a, 1997b) empregados para determinacao do campo de velocidades e geracao de
imagens com maior resolucao que a sısmica convencional.
TOMOGRAFIA SOBRE DIFRACOES E OPERADORES FOCAIS
Berkhout (1997a, 1997b) lancou o conceito de operadores de focalizacao (Common Focus
Point Operator), CFPO, que representam a funcao de Green de uma fonte virtual impulsiva
deflagrada em profundidade e registrada na superfıcie de aquisicao. Thorbecke (1997) e
19
Bolte (2003) demonstram que os operadores focais sao passıveis de obtencao a partir de
dados sısmicos de reflexao. Assim, por volta da segunda metade da decada de 1990 seria
possıvel empregar, nao os tradicionais tempos de transito das reflexoes, mas os tempos
simples de funcoes de Green ascendentes - dos pontos sobre uma superfıcie refletora ate a
superfıcie.
Nos anos 1990 destacam-se os trabalhos de Hegge e Fokkema (1996), Hegge (1997),
Hegge et al. (1998), Hegge e Bolte (1999) e Kabir (1997) versando sobre inversao tomografica
a partir de operadores focais. As teses de doutorado de Kabir (1997), Hegge (2000) e Cox
(2004) empregam a inversao tomografica sobre CFPO. Embora os princıpios envolvidos
ja fossem utilizados na sismologia para estimativa dos focos de terremotos e o campo de
velocidades crustal, os trabalhos de Kabir (1997), Hegge (2000) e Cox (2004) empregam a
inversao tomografica sobre CFPO (funcoes de Green) no ambito da sısmica para petroleo.
Neste processo os dados de entrada D sao os tempos de transito dos operadores focais e
na inversao almeja-se determinar a posicao dos difratores e o campo de velocidades, que
compoem o espaco de modelos m. O presente trabalho desenvolve-se tendo como ponto de
partida as teses citadas.
INVERSAO DO CAMPO DE ONDAS COMPLETO
Na decada de 2010, o emprego da tecnica Full Wavefield Inversion (FWI), Inversao
do campo de ondas completo, cresceu rapidamente no meio academico e, sem seguida, na
industria do petroleo. Embora as bases teoricas da tecnica remonte ao inıcio da decada de
1980 (Tarantola, 1984), somente nos anos 2010, a ciencia computacional ofereceria poder
de processamento que permitisse o uso industrial do FWI. O FWI e um metodo de inversao
iterativo onde os dados sao os sismogramas registrados, e o modelo a ser atualizado, o
20
campo de velocidade. O modelo de velocidade e atualizado a cada iteracao a partir do
resıduo depropagado obtido com a diferenca entre sismogramas observado e modelado. Esta
definicao, resume de forma simploria o FWI. Maiores detalhes sobre o metodo podem ser
obtidos em etc. etc. etc. - BIBLIOGRAFIAS FWI !!!!
21
5- MODELAGEM SISMICA POR TRACADO DE RAIOS
Seja um modelo de velocidade v(x, y, z) e dois pontos A e B separados espacialmente.
Deseja-se computar o tempo τ de um raio que parta de A e chegue a B.
Considerando que a trajetoria, seja γ, tenha um comprimento s, o tempo de percurso
sera dado por
τ =
∫
s
ds
v(x, y, z), (3)
onde ds e o segmento de raio ao longo da trajetoria s.
Ainda nao se conhece a trajetoria γ percorrida por este raio, mas de acordo com o
Princıpio de Fermat (1601 a 1665), o tempo de A para B deve ser mınimo. Em outras
palavras, o sinal se propaga do ponto A ao ponto B ao longo da curva que fornece um valor
estacionario para a integral descrita pela equacao 3.
Considerando r como o vetor posicao de cada ponto da curva γ por uma variavel σ
tem-se:
r(σ) = (x(σ), y(σ), z(σ)), (4)
ds pode ser calculado por:
ds =√
(drdr), (5)
ou ainda
ds =
√
(dr
dσ
dr
dσ)dσ, (6)
Substituindo
22
ds =
√
(dx
dσ)2 + (
dy
dσ)2 + (
dz
dσ)2dσ, (7)
Substituindo em 3
τ =
∫
s
√
(dxdσ
)2 + ( dydσ
)2 + ( dzdσ
)2dσ
v(x, y, z), (8)
Considerando
L(x, y, z, x2, y2, z2) =
√
(dxdσ
)2 + ( dydσ
)2 + ( dzdσ
)2
v(x, y, z). (9)
O tempo τ e um funcional dado por
τ =
∫
sLdσ. (10)
Busca-se o extremo do funcional τ , que, portanto, obedece a condicao
δτ = δ
∫
sLdσ = 0. (11)
Para atender 11 as equacoes de Euler-Lagrange devem ser satisfeitas:
d
dσ(∂L
∂x) − ∂L
∂x= 0; (12)
d
dσ(∂L
∂y) − ∂L
∂y= 0; (13)
d
dσ(∂L
∂z) − ∂L
∂z= 0. (14)
Considerando:
px =∂L
∂x; (15)
py =∂L
∂y; (16)
pz =∂L
∂ze, (17)
23
observando as devidas substituicoes chega-se a
dx
ds= pxv(x, y, z), (18)
dy
ds= pyv(x, y, z), (19)
dz
ds= pzv(x, y, z), (20)
e
dpx
ds=
∂
∂x(
1
v(x, y, z)), (21)
dpy
ds=
∂
∂y(
1
v(x, y, z)), (22)
dpz
ds=
∂
∂z(
1
v(x, y, z)), (23)
para o parametro de raio igual ao espaco s.
Considerando-se o tempo τ como parametro de raio chegamos ao seguinte sistema.
dx
dτ= pxv
2(x, y, z), (24)
dy
dτ= pyv
2(x, y, z), (25)
dz
dτ= pzv
2(x, y, z), (26)
e
dpx
dτ=∂ − ln(v(x, y, z))
∂x, (27)
dpy
dτ=∂ − ln(v(x, y, z))
∂y, (28)
dpz
dτ=∂ − ln(v(x, y, z))
∂z. (29)
24
6- FUNDAMENTOS DE TOMOGRAFIA SISMICA
6.1- Introducao
6.2- Fundamentacao matematica empregada
6.3- Exercıcios
Conceitos basicos:
Filosofia D = Lm
Classificacao de problemas inversos
- Lineares
- Nao Lineares
Formas de resolucao do problema inverso
Deducao
Regularizacao
Conceitos basicos
*Funcional: Define-se como funcional, toda funcao cujo domınio e um espaco vetorial e a
imagem e um escalar.
* Norma de uma matriz L: A norma de L mede o maior valor que um vetor x (autovetor
ou nao) e amplificado quando se tem um produto L.x. A norma de L e igual a raiz quadrada
do maior autovalor de
norma =LTL
||L|| (30)
25
* Autovalor e Autovetor
Seja uma matriz quadrada L. Os autovalores desta matriz sao os valores de λ que atendem
a igualdade
L − λI = 0. (31)
Se a dimensao de L e n, a resolucao do determinante da equacao 31 nos fornece n λs.
Os autovetores de L sao calculados resolvendo-se o sistema
[L − λI].x = 0, (32)
para x.
* Numero de condicao de uma matriz: Razao entre o maior autovalor e o menor auto-
valor.
* Base: A base de um espao de vetores e um conjunto de vetores que devem ser linear-
mente independentes.
* Dimensao de um espaco de vetores: Qualquer base para o espaco V contem um mesmo
nmero de vetores, o qual e denominado dimensao de V. A dimenso diz respeito aos graus
de liberdade do espaco.
* QUANTO A SOLUCAO DO PROBLEMA INVERSO:
- Existencia: Dado um sistema d=L.m, L deve ser inversıvel.
- Robustez: A robustez diz respeito ao nıvel de sensibilidade da solucao face a presenca
26
de um pequeno nmero de erros grandes (ruıdos ou outliers) nos dados. Quanto maior a
robustez menos sensıvel e a solucao com relacao aos outliers.
- Unicidade: Uma solucao e dita unica se ao mudarmos de um modelo m1 para m2, os dados
tambm mudam de d1 para d2.
- Estabilidade: A estabilidade indica como pequenos erros nos dados propagam-se para o
modelo. A instabilidade pode induzir a nao unicidade.
* CLASSIFICACAO DO PROBLEMA INVERSO - Problema bem posto: Problema
d=L.m que possui solucao unica e estavel.
- Problema mal posto: Problema d=L.m tal que nao possui pelo menos uma das carac-
terısticas: unicidade, estabilidade e existencia.
* Regularizacao: Tecnicas para transformar problemas mal postos em problemas bem
postos.
* Resıduo: Diferenca entre o dado real e o dado calculado.
* Erro: Diferenca entre o modelo real e o modelo calculado.
Se o problema direto e dado por
d = Lm, (33)
O problema inverso e dado por
m = L−1d, (34)
27
Regressao linear
Seja o resıduo dado pela equacao a seguir.
r = dobs − dcalc. (35)
Considerando a equacao 33, a equacao 35, nos leva a
r = dobs − (Lm). (36)
Por ora, apresenta-se a solucao de m considerando o modulo do resıduo, norma 1, ou o
quadrado do resıduo, norma 2.
Funcional com norma 1 do resıduo
O emprego da norma 1 de r nos leva a
||d − Lm||1 =N
∑
i=1
|di − (Lm)i)|. (37)
Trata-se de uma solucao menos sensıvel aos outliers (dados altamente inconsistentes).
A solucao deste funcional e complexa. Uma alternativa pratica e o emprego de IRLS
(Iteratively Reweighted Least Square). Ela resolve uma sequencia de mınimos quadrados
ponderados cujas solucoes convergem para solucao de mınimos para norma 1.
Se
f(m) = ||r||, (38)
f(m) =N
∑
i=1
|ri|. (39)
28
A derivada de f em relacao e dada por
∂f(m)
∂mk
=N
∑
i=1
∂|ri|∂mk
, (40)
que leva a
∂f(m)
∂mk
=N
∑
i=1
Li,ksgn(ri). (41)
Considerando
sgn(ri) =ri
|ri|, (42)
a equacao 41 torna-se
∂f(m)
∂mk
=N
∑
i=1
Li,k
1
|ri|ri. (43)
Assim, o gradiente de f e dado por:
∇f(m) = LTRr = LTR(d − Lm). (44)
R e uma matriz diagonal de pesos dada por:
Ri, i =1
|ri|. (45)
Continuando, a busca pelo mınimo nos leva a
LTR(d − Lm) = 0, (46)
de onde se obtem
LTRLm = LTRd. (47)
Como R depende de m, a equacao 47 constitui um sistema nao linear que nao pode
ser resolvido diretamente. A seguirm apresenta-se a solucao da equacao 47 empregando a
29
tecnica IRLS. Trata-se de um algoritmo iterativo simples para calcular os pesos R e m.
Passo 1: m0 = mL2 (chute inicial)
Passo 2: r0 = d− Lm0
Passo 3: LTRLm = LTRd
Passo 4: mk+1
Passo 5: Teste ||mk+1−mk||21+||mk+1||2
< τ . Se sim siga para (6), se nao retorne para (2),
Passo 6: mfim
Funcional com norma 2 do resıduo
O emprego da norma 2 de r nos leva a
||d − Lm||2 =
√
√
√
√
N∑
i=1
di − (Lm)i)2. (48)
A solucao por mınimos quadrados minimiza o funcional dado pela norma 2 do resıduo:
m = [LTL]−1LTd. (49)
Trata-se de uma boa solucao, pois e estatisticamente mais provavel se os erros embutidos
nos dados tem distribuicao normal.
Metodos de resolucao do problema inverso nao linear
Abordam-se neste topico os principais metodos de resolucao de Problemas Inversos
fracamente nao lineares por meio de informacoes advindas ou relacionadas ao gradiente da
funcao objetivo.
30
Metodo da maxima declividade
Sendo o gradiente da funcao objetivo, ∇Θ = ∂Θ∂m
, busca-se o mınimo iterativamente pela
equacao
mk+1 = mk − a∇Θ, (50)
onde a e um escalar positivo. Ou seja, o modelo no passo k+1 e obtido a partir do modelo no
ponto k, mk na direcao contraria ao gradiente da funcao objetivo ∇Θ. A figura 7 ilustra a
funcao objetivo e o caminho das iteracoes, entre o modelo inicial m0 e a solucao m, segundo
o metodo da maxima declividade. Este constitui o metodo de gradiente mais simples, exibe
baixa taxa de convergencia (Flechter, 1980), contudo ela e assegurada (Aster, 2005), desde
que o modulo do gradiente seja ajustado segundo uma constante ad hoc.
31
Figure 7: Funcao objetivo com duas variaveis m1 e m2 ilustrando as iteracoes pelo metodo
de maxima declividade ate a convergencia para solucao em m (modificado de Shewchuck,
1994).
32
Metodo de Newton
A expansao em serie de Taylor (truncada no termo de ordem 2) da funcao objetivo nas
vizinhancas de mk fornece:
Θ(mk+1) = Θ(mk) + ∇Θ(mk)[mk+1 −mk] +1
2∇2Θ(mk)[mk+1 −mk]
2. (51)
Considere ∆m = mk+1 −mk. Almeja-se o ponto de mınimo de 51 em relacao a ∆m.
Suprimindo-se os passos intermediarios obtem-se a resolucao do modelo mk+1 pelo metodo
de Newton (Tarantola, 2002) por:
mk+1 = mk − [∇2Θ(mk)]−1∇Θ(mk) (52)
A aplicacao do metodo de Newton para a resolucao de Problema Inverso demanda o
calculo do Hessiano ∇2Θ ou ∂2Θ∂m2 , a derivada segunda da funcao objetivo em relacao ao
modelo m. Em alguns problemas da geofısica este calculo e complexo. Adicionalmente nao
ha garantia de convergencia (Aster et al., 2005). Para problemas passıveis de linearizacao
o metodo de Gauss-Newton e mais factıvel.
Metodo de Gauss-Newton
A deducao do metodo em escopo pode ser encontrada em diversos livros textos, entre eles,
Flechter (1980), Menke (1984), Nocedal e Wright (2000) e Aster et al. (2005), contudo ela e
apresentada a seguir completa e comentada para compreensao do leitor e futuras consultas
do autor.
Avancando na equacao 35 obtem-se:
Θ =1
N
N∑
i=1
[
dobsi − dcalc
i
]2(53)
33
Expandindo o termo dcalci da equacao acima em serie de Taylor nas proximidades de um
modelo m0,
Dcalci = Dcalc
i (m0) +M∑
j=1
∂di(m0)
∂mj
∆mj, (54)
onde ∆mj e uma perturbacao no modelo nas proximidades de m0. Derivando a equacao
53 em relacao ao modelo,
∂Θ
∂mk
=2
N
N∑
i=1
(
dobsi − dcalc
i
) ∂
∂mk
(
dobsi − dicalc
)
(55)
No ponto de mınimo ∂Θ∂mk
= 0. Adicionalmente, os dados observados, dobsi , nao dependem
do modelo m. Assim a equacao 55 torna-se
0 =N
∑
i=1
(
dobsi − dcalc
i
) ∂dcalci
∂mk
(56)
Substituindo a equacao 54 em 56 obtem-se:
N∑
i=1
dobsi − dcalc
i (m0) −M∑
j=1
∂di(m0)
∂mj∆mj
∂dcalci
∂mk
= 0 (57)
Considerando dobsi - dcalc
i (m0) = ∆di:
N∑
i=1
∆di,k −M∑
j=1
∂dcalci
∂mj∆mj
∂dcalci
∂mk
= 0. (58)
Rearranjando os termos,
N∑
i=1
∆di,k∂dcalc
i
∂mk
−N
∑
i=1
M∑
j=1
∂dcalci
∂mj
∆mj∂dcalc
i
∂mk
= 0. (59)
A equacao 59 pode ser escrita na forma matricial:
LT ∆d− LTL∆m = 0, (60)
34
ou,
LTL∆m = LT ∆d. (61)
Aplicando (LTL)−1 em ambos os lados:
∆m = (LTL)−1LT ∆d. (62)
A equacao 62 e a forma basica mais empregada nos Problemas Inversos considerando-se
a funcao objetivo com norma 2. Ela assemelha-se a solucao da equacao 52 aproximando
o Hessiano pelo termo LTL (Nocedal e Wright, 2000). Contudo, ha que se relembrar
duas premissas embutidas em sua solucao: (1) A modelagem e o operador L que regem o
fenomeno fısico devem ser passıveis de linearizacao nas vizinhancas de um modelo m0; e
(2) O modelo m0 deve estar proximo da solucao.
A implicacao direta da primeira premissa e que, nas aproximacoes extremamente nao
lineares o metodo falha. Embora a direcao apontada pelos cossenos diretores de ∆m nas
vizinhancas de m0 sejam corretas, seu modulo pode superar as condicoes de linearizacao
embutidas na expansao em Taylor. Neste caso, se a magnitude de ∆m e maior que este
limite, deve-se reduzir o seu modulo utilizando de um parametro ad-hoc, seja ǫ. Este
parametro e positivo, menor que 1 (Tarantola, 2002) e pode variar em cada iteracao. Assim,
a equacao 62 torna-se
∆m = ǫ(LTL)−1LT ∆d, (63)
aumentando as chances de convergencia. Um valor desnecessariamente baixo de ǫ pode
tornar a convergencia extremamente lenta.
35
Metodo de Levenberg-Marquardt
A aproximacao do Hessiano por LTL contribui para as incertezas na magnitude de
∆m no metodo de Gauss-Newton. No metodo de Levenberg-Marquardt evita-se o emprego
da constante ǫ, uma especie de busca linear (line search), empregando-se uma estrategia
de regiao de confianca (Nocedal e Wright, 2000). Esta regiao de confianca e obtida pela
da adicao de uma matriz diagonal, fruto do produto de uma constante ψ com a matriz
identidade I, a matriz LTL.
∆m = (LTL + ψI)−1LT ∆d, (64)
Uma vantagem adicional dessa estrategia e o fato de ela contornar possıveis deficiencias
no rank (numero de linhas ou colunas linearmente independentes) da matriz L, uma das
fraquezas do metodo de Gauss-Newton (Nocedal e Wright, 2000). Em outras palavras,
o metodo em escopo forca que a matriz LTL seja positiva definida garantindo sua nao
singularidade.
Pseudo-Inversa de Moore-Penrose
Para resolver a equacao 34, e necessario inverter a matriz L. Contudo, caso esta matriz
nao seja inversıvel, em lugar de se calcular m ou ∆m pela equacao 34, L e decomposta em
valores singulares (Singular value decomposition (SVD)) (Aster et al. 2005).
A matriz com valores singulares, Σ, e uma matriz diagonal com valores positivos e
decrescentes de acordo com o aumento de i (Σi=1 > Σi=2 > Σi=3 > ..). Ela pode possuir
alguns elementos muito baixos ou nulos. Se somente os p primeiros elementos nao nulos de
Σ forem considerados, esta matriz pode ser particionada (Aster et al., 2005).
36
A equacao resultante e dada por
L = WpΣpVTp . (65)
O produto WpΣpVTp e passıvel de inversao de acordo com Moore (1920) e Penrose
(1950). Assim a inversa generalizada de L ou pseudo-inversa de Moore e Penrose (Aster et
al., 2005) e calculada por
L† = VpΣ−1p WT
p . (66)
No lugar da equacao 34 emprega-se
m = L†d, (67)
desde que o posto de L seja igual ao numero de parametros do modelo m e inferior ao
numero de registros ou dados n.
Em problemas fracamente nao lineares a mesma relacao e valida para perturbacoes no
modelo e no dado,
∆m = L†∆d, (68)
, e a solucao para o modelo m e alcancada iterativamente,
mk+1 = mk + ∆m, (69)
ate que dado modelado e dado observado se ajustem.
Metodo dos gradientes conjugados
O metodo em escopo faz parte da famılia de metodos denominada direcoes conjugadas
(dois vetores di e dj sao conjugados ou L-ortogonais, se diLdj = 0, onde L e uma matriz
compatıvel com o produto (Shewchuck, 1994)). Tratam-se de metodos aplicaveis a matrizes
37
positiva definidas (para qualquer vetor x diferente de zero tem-se xTLx > 0) e que sao mais
efetivos no tratamento de matrizes esparsas (Shewchuck, 1994).
O ponto de partida dos metodos de direcoes conjugadas nao difere do metodo de gradi-
entes. Contudo, enquanto neste ultimo observam-se varios passos alternados seguindo uma
mesma direcao ate que se atinja a solucao (Figura 7), os metodos das direcoes conjugadas
empregam buscas lineares (line search) em direcoes especıficas, e, ao longo destas, busca-se
o mınimo. Este processo reduz o numero de iteracoes.
No metodo dos gradientes conjugados, descrito no pseudo-algoritmo a seguir, algoritmo
74, a direcao de busca e construıda atraves da ortogonalizacao dos resıduos segundo a matriz
L (L-ortogonalizacao) (Shewchuck, 1994).
∆m0 = r0 = Dobs − Lm0
loop while (r > tol)
αk =rTk rk
∆mTkL∆mk
(70)
mk+1 = mk + αk∆mk (71)
rk+1 = rk − αkL∆mk (72)
βk+1 =rTk+1rk+1
rTk rk
(73)
∆mk+1 = rk+1 + βk+1∆mk (74)
end loop
Na primeira iteracao assume-se que a atualizacao do modelo (∆m0) e igual ao resıduo
(r0). A constante α calculada a cada iteracao k pondera o modulo da atualizacao por
38
(αk∆mk). O calculo βk culmina o processo de ortogonalizacao apos o qual calcula-se o
valor efetivo da atualizacao do modelo ∆mk+1. O processo iterativo continua ate que o
resıduo rk esteja abaixo de uma tolerancia tol.
Regularizacao
Grande parte dos problemas inversos sao mal postos. Neste topico sao apresentados metodos
que tendem a melhorar as condicoes de resolucao, sobretudo para reduzir sua instabilidade,
entre eles a regularizacao (processo no qual sao impostos vınculos e restricoes que induzem
a solucao (Aster et al., 2005)). De acordo com a bibliografia consultada (Flechter, 1980;
Nocedal e Wright, 2000; Aster et al., 2005) os processos voltados para reducao da instabil-
idade e inducao a solucao podem atuar:
- na parametrizacao do modelo;
- na concepcao da funcao objetivo;
- na manipulacao da matriz de sensibilidade (tecnicas SVD);
Descrevem-se a seguir as diferentes formas voltadas para geracao de problemas inversos
bem postos.
Parametrizacao do modelo
Aster et al. (2005) citam que pela selecao cuidadosa da dimensao da malha (grid) e
possıvel produzir um sistema de equacoes bem condicionado. Os autores nomeiam esse
procedimento de regularizacao por discretizacao.
Afirma-se que as parametrizacoes, por camadas constitui forma de parametrizacao util
39
para descricao de modelos de velocidade em meios sedimentares. Neste caso o vazio de
informacoes gerado pela discretizacao por celulas e contornado. Adicionalmente, elas em-
pregam um menor numero de parametros pois simulam o comportamento esperado das
rochas em ambientes sedimentares.
Regularizacao atuando na concepcao da funcao objetivo
A funcao objetivo mais simples possıvel e alcancada atraves da equacao 53. A insercao de
termos adicionais na funcao objetivo e o instrumento para aplicacao de vınculos e restricoes.
Se a funcao objetivo Θ for adicionada uma funcao de restricao ou vınculos ψ1c1(m), onde
ψ1 e um escalar, obtem-se uma nova funcao
L(m,ψ1) = Θ(m) + ψ1c1(m) (75)
Adicionando mais funcoes ci ponderadas por ψi obtem-se:
L(m,ψi) = Θ(m) +n
∑
i=1
ψici(m). (76)
Demonstra-sebf que nos pontos de solucao m∗ existem escalares ψ∗i tais que ∇mL(m∗, ψ∗
i ) =
0 (Nocedal e Wright, 2000). Assim, nos pontos estacionarios de L(m∗, ψ∗i ) observa-se:
∇mΘ(m∗) = −n
∑
i=1
ψ∗i ∇mci(m
∗). (77)
Os escalares ψ∗i sao os parametros de regularizacao.
Sobre as funcoes ψici citadas de forma generica no paragrafo anterior e que devem ser
construıdas as restricoes e vınculos que relacionam o processo matematico de inversao com
a realidade geologica. Uma das aplicacoes e a imposicao de uma relacao entre o modelo mk
em uma iteracao k com um modelo de referencia mref :
c1(m) = ‖Wm(m− mref )‖2 (78)
40
onde Wm e uma matriz MxM (onde M e o numero de parametros do modelo) que pondera
o vetor diferenca mk −mref . Na tomografia cinematica para obtencao do campo de veloci-
dade, o modelo de referencia mref pode ser obtido com o conhecimento do comportamento
mais frequente das rochas:
- o meio geologico sedimentar apresenta maior conteudo de altas frequencias espaciais na
direcao vertical que nas direcoes horizontais;
- devido a compactacao, um mesmo litotipo em uma profundidade Z tera menor porosidade
em uma profundidade Z+dZ - desde que dZ > 0. Observadas as condicoes de compactacao
em meio drenado, a velocidade do meio equivalente em Z + dZ sera maior que em Z.
Assim, o modelo de referencia toma a forma de uma funcao do tipo mref = mref (Z)
em situacoes onde o conhecimento de subsuperfıcie e limitado. Supondo que o modelo de
referencia e o modelo inicial, a solucao iterativa sera dada por
∆m = (WdLTL+ ψWm)−1[LTWd∆D + ψWm(m−mref )] (79)
ondeWd e um operador para ponderar os dados. Quando ha mais de uma funcao de restricao
a equacao 79 generaliza-se para:
∆m = (WdLTL+
p∑
i=1
ψiWmi)−1[LTWd∆D +
p∑
i=1
ψiWmi(m−mref )] (80)
Na pratica a regularizacao em epıgrafe adiciona valores aos termos da matriz LTL
(equacao 79).
E oportuno destacar o efeito que a regularizacao causa no processo de inversao. Se o
valor de ψ for muito maior que os termos da matriz LTL, a soma (LTL+ψI) tende a ψI e
a equacao 53 se transforma em
∆m = (ψWm)−1LT ∆D. (81)
41
O termo LT ∆D e o gradiente do resıduo (∆D). Depreende-se que o sistema tende ao
metodo de maior declividade (Aster et al., 2005) ponderado por ψ−1.
De forma contraria, se o valor ψ for muito menor que os termos da matriz LTL, nao se
observa efeito de regularizacao e o sistema tende para o metodo de Gauss-Newton, equacao
62.
E desejavel que a adicao de valores a funcao objetivo, a torne cada vez mais concava.
Um funcional que possua calhas e pontos de cela desestabiliza o processo de inversao alem de
proporcionar multiplas solucoes. Assim a adicao de uma funcao que, literalmente, deforme o
funcional, melhora a estabilidade e reduz a ambiguidade. Este fato e observado no exemplo a
seguir com um modelo representado por dois parametros, m1 e m2, e um funcional calculado
por
Θ = m21. (82)
A Figura 8 exibe o funcinal Θ da equacao 82. Observa-se que, neste exemplo, Θ dispoe-
se segundo uma calha vertical com uma serie de mınimos ao longo do eixo m2. Com este
funcional a solucao sempre depende do modelo inicial (m1o ,m2o). Para um modelo inicial
(m1o ,m2o) =(2.0,2.0) a solucao e (0.0,2.0) e; para um modelo inicial (m1o ,m2o) =(-2.0,-2.0)
a solucao e (0.0,-2.0) (Figura 8). Neste caso ha solucoes ambıguas.
42
Fu
nc
tio
na
l v
alu
e
m1
m2
Figure 8: Funcao objetivo com duas variaveis m1 e m2 dada pela equacao
82. Ilustra-se a convergencia para o mınimo a partir de dois modelos ini-
ciais (m1o ,m2o) = (2.0,2.0) e (m1o ,m2o) = (-2.0,-2.0) gerando as respectivas
solucoes (0.0,2.0) e (0.0,-2.0). Neste caso ha solucoes ambıguas.
43
Se na equacao 82 adicionar-se um termo quadratico m22, o novo funcional Θ′ sera
Θ′ = m21 +m2
2. (83)
Sua forma esta representada na Figura 9. A solucao e unica, qualquer que seja o modelo
inicial (m1o ,m2o). Partindo dos modelos iniciais (m1o ,m2o) =(2.0,2.0) ou (m1o ,m2o) =(-2.0,-
2.0) a solucao econtra-se sempre em (m1o ,m2o) =(0.0,0.0). Nao se verifica ambiguidade da
solucao.
Fu
nc
tio
na
l v
alu
e
m1
m2
Figure 9: Funcao objetivo com duas variaveis m1 e m2 dada pela equacao
82. Ilustra-se a convergencia para o mınimo a partir de dois modelos iniciais
(m1o ,m2o) = (2.0,2.0) e (m1o ,m2o) = (-2.0,-2.0) gerando a mesma solucao
(0.0,0.0). Neste caso nao ha ambiguidade.
44
Neste topico foram apresentadas as formas genericas de regularizacao atraves da con-
cepcao da funcao objetivo.
Regularizacao atuando sobre a matriz de sensibilidade - Tikhonov de or-
dem zero
A pseudo-inversa calculada com a equacao 68 nao possui projecao sobre o espaco nulo de
L (N(L)). Contudo alguns valores singulares podem ser tao pequenos que, com a inversao
de Σ (equacao 66), estes podem gerar solucoes anomalamente altas. Na pior situacao eles
podem ser amplificadores de ruıdos (Aster el al. (2005)).
Segundo a condicao de Picard (Aster et al., 2005), quando o produto entre a matriz UT
(obtido pela decomposicao SV de L, secao ) e o vetor de dados D ((UT D)i) decai mais
rapidamente que Σi ha estabilidade da solucao. A plotagem do produto (UT D)i contra Σi
constitui uma poderosa ferramenta de analise para estabilidade com a decomposicao de val-
ores singulares. Caso a condicao de Picard nao seja atendida, a solucao consiste em eliminar
os valores singulares extremamente baixos ao preco de um menor ajuste do modelo. Esse
procedimento e denominado regularizacao por truncamento dos valores singulares (TSVD).
Com a regularizacao de Tikhonov (Tikhonov e Arsenin, 1977) busca-se uma solucao que,
alem de fornecer o menor resıduo, tambem apresente a menor norma. Na sua deducao, que
pode ser conduzida a partir dos multiplicadores de Lagrange (Aster et al., 2005), elevam-se
os valores singulares na busca por maior estabilidade. Na regularizacao de ordem zero, cada
valor singular Σi deve ser multiplicado por um fator de filtro (Aster et al., 2005) fi tal que
fi =Σ2
i
Σ2i + ψ2
i
, (84)
onde ψ e o parametro de regularizacao.
45
Aplicacao de restricoes na matriz de sensibilidade e controle da norma da perturbacao no
modelo
Em problemas inversos fracamente nao lineares, as perturbacoes no modelo relacionam-
se linearmente com as perturbacoes nos dados 68 (Sen, 2006). A atualizacao no modelo e
calculada pela equacao 69 e a solucao final e obtida quando o resıduo se torna menor que
um valor pre-determinado de tolerancia.
Ainda para problemas fracamente nao lineares, a matriz de sensibilidade L deve ser cal-
culada em cada iteracao assim como sua pseudo-inversa L†. Quando o espaco nulo de dados
[N(LT )] e nao trivial e o espaco nulo de modelos [N(L)] e trivial, a solucao por equacoes
normais [(LTL)−1] e obtida pela pseudo-inversa de Moore-Penrose sao exatamente iguais.
Apesar da solucao ser unica, os dados nao se ajustam perfeitamente. Quando [N(LT )] e
[N(L)] sao nao-triviais, a solucao utilizando a pseudo-inversa e um mınimo quadratico e
tambem possui menor norma. Nesse caso a solucao nao e unica e os dados nao se ajustam
perfeitamente (Aster et al., 2005).
O Problema Inverso estudado nesta tese pertence a uma das situacoes citadas no ultimo
paragrafo. O processo iterativo descrito pelas equacoes 68 e 69 possui as caracterısticas
do metodo de Gauss-Newton. A ausencia de calculo explıcito do Hessiano pode causar
inacuracia na estimativa da perturbacao do vetor modelo ∆m. Consequentemente, uma
norma incorreta de ∆m pode causar nao convergencia (Figura 10a). Neste exemplo, em
um modelo bidimensional m, com norma inadequada, nao se observa convergencia (Figuras
10a e 10b).
46
0.0
0.2
0.4
0.6
0.8
0.01
0.21
11019876543210
rebmun noitaretI
Re
sid
ua
l
Re
sid
ua
l
(a) (b)
0.0
0.2
0.4
0.6
0.8
0.01
0.21
11019876543210
rebmun noitaretI
Re
sid
ua
l
Re
sid
ua
l
(c) (d)
Figure 10: (a) Funcao objetivo com duas variaveis m1 e m2. A linha negra tracejada
representa o caminho da solucao ieterativa. A convergencia nao e observada. (b) Resıduo
em cada iteracao observado em (a). (c) Funcao objetivo com duas variaveis m1 e m2. A
linha negra tracejada representa o caminho da solucao iterativa. A convergencia para o
mınimo e suave e direta. (d) Resıduo em cada iteracao observado em (c).
47
Portanto, um fator de estabilizacao Ω e necessario para reduzir a norma do vetor ∆m.
Assim, a equacao 68 se torna
∆m =1
ΩL†∆D. (85)
Com o mesmo problema da Figura 10a, a norma do vetor ∆m e reduzida por um fator Ω
e a convergencia e assegurada (Figuras 10c e 10d).
O fator Ω demanda calibracao na primeira iteracao ao modo de uma busca linear (linear
search), mas na iteracoes seguintes Ω e ajustado. Se o resıduo em um passo subsequente
k+ 1 e maior que no passo k, um fator de reducao f e aplicado sobre Ω−1. Diferentemente,
se em uma iteracao posterior o resıduo e menor que na iteracao anterior, Ω−1 e aumentado
segunto este fator f . O melhor valor de f empregado nos Problemas Inversos encontrados
neste trabalho e igual a 5%.
Adicionalmente, restricoes devem ser impostas na matriz de sensibilidade para incor-
porar informacoes sobre o modelo. Santos et al. (2010) apresentam a ponderacao de cada
coluna da matriz L multiplicando-a por uma constante b variando entre 0 e 1 e encerrada
na matriz diagonal B. Parametros conhecidos e, portanto restritos, tem as correspondentes
colunas da matriz L multiplicados por b proximas ou iguais a zero. Parametros desconheci-
dos e nao restritos tem os correspondentes termos da matriz de sensibilidade multiplicados
por b igual a 1. Assim, em cada iteracao, a perturbacao do modelo e calculada por
∆m =1
Ω(LB)†∆D. (86)
Regularizacao na pratica
Retornemos ao problema linear direto d = Lm. Se a solucao e dada por
m = [LTL]−1LTd, (87)
48
o regularizador R e aplicado da seguinte forma:
m = [LTL + λRTR]−1LTd, (88)
onde R precisa ter o mesmo numero de colunas de L e qualquer numero de linhas.
Podemos considerar uma matriz L.
L L L L ... L L L L L
L L L L ... L L L L L
L L L L ... L L L L L
L L L L ... L L L L L
R R R R ... R R R R R
R R R R ... R R R R R
R R R R ... R R R R R
R R R R ... R R R R R
e um vetor d como
d
0
onde 0 possui o mesmo numero de linhas de R. Com estas condicoes existe uma solucao
por mınimos quadrados de d=L m.
49
7- FORMAS DE DISCRETIZACAO DO MEIO FISICO PARA
TOMOGRAFIA
7.1- Introducao
7.2- Polıgonos
7.3- Splines
7.4- Funcoes harmonicas
7.5- Funcoes de base radial
7.6- Por camadas
7.7- Exercıcios
Este topico dedica-se as formas de representacao do campo de propriedades m. Este
tema, por si, constitui ampla area de pesquisa e nesta secao apenas as principais formas
basicas de representacao 2D sao revisitadas.
Discretizacao por blocos ou celulas em malha regular
A discretizacao por blocos e o mais intuitivo e simples dos metodos de discretizacao.
Basicamente ela consiste em dividir o espaco de modelos em uma malha, preferencial-
mente regular, composta por diminutos polıgonos encerrando propriedades homogeneas, ou
variaveis, contınuas. A Tabela 2 exibe a representacao pictorica de um modelo por uma
malha regular de retangulos.
Este tipo de discretizacao e o mais adequado para modelos extremamente complexos,
sempre que a dimensao das celulas for menor que os menores detalhes do modelo. Sob estas
condicoes, qualquer combinacao de estruturas presentes na Figura 3 pode ser representada
50
por celuas. Bishop et al. (1985) e Fei e McMechan (2006), entre outros, usam a discretizacao
por celulas para estimativa da velocidade.
C1 C2 C3 C4 C5
C6 C7 C8 C9 C10
C11 C12 C13 C14 C15
C16 C17 C18 C19 C20
C21 C22 C23 C24 C25
Table 1: Representacao pictorica de um modelo discretizado por 25 retangulos.
Discretizacao por blocos ou celulas em malha nao regular
Malhas nao regulares constituem uma forma de discretizacao bastante atraente. Os
polıgonos nao possuem dimensoes constantes, assim, e possıvel empregar um numero re-
duzido de elementos nas regioes com pouca complexidade estrutural, e elevar o seu numero
concomitantemente a reducao de suas dimensoes nas areas mais complexas e/ou de eleva-
dos gradientes de velocidade. Rotinas de modelagem ou inversao implementadas segundo o
metodo dos elementos finitos permitem o emprego de malhas nao regulares (Cook, 1995).
Cox (2004) empregou triangulos de Delaunay para discretizar o campo de velocidades no
processo de inversao tomografica sobre operadores focais.
51
Discretizacao por polinomios
Nesta forma de discretizacao, sem termos cruzados, o campo de velocidades e represen-
tado pela funcao polinomial
v(x, z) = v0 +n
∑
i=1
Aixi +
m∑
j=1
Bjzj , (89)
onde v(x, z) e o campo de velocidade e, Ai e Bj os respectivos coeficientes aplicados as
coordenadas x e z e v0 o termo independente.
Equacoes de primeiro grau, v(z) = c0 +B1z , aliam a funcao compactacao (progressiva
perda de porosidade com o aumento da pressao litostatica ou profundidade) com a funcao
velocidade. Em outras palavras, para um mesmo litotipo, quanto menor a porosidade, maior
sera a velocidade de ondas compressionais. Uma serie de leis empıricas relaciona velocidade
de ondas compressionais com a porosidade. Para maiores detalhes sugere-se a leitura de
Mavko et al. (2009). Um exemplo de aplicacao do emprego de polinomio para descricao
do campo de velocidades encontra-se em Siddiqui et al. (2004). Os autores constroem
um modelo de velocidades consistente empregando equacoes de primeiro grau para a pilha
sedimentar clastica do Golfo do Mexico.
Discretizacao por funcao exponencial ou logarıtmica
Embora pouco comum, unidades sedimentares nao deformadas podem ser representadas
localmente por funcoes exponenciais ou logarıtmicas para grandes pacotes sedimentares
homogeneos. Assim uma grande sequencia arenıtica pode ser representada por uma funcao
exponencial
v(x, z) = v0 +B1e(qz), (90)
52
se ajustando ao inverso da funcao compactacao, onde q, v0 e B1 sao coeficientes que depen-
dem da rocha em estudo.
Algumas litologias, compostas predominantemente por graos nas fracoes silte e/ou
argila, perdem a porosidade muito rapidamente (com elevado gradiente) nas primeiras cen-
tenas de metros e apresentam gradientes mais suaves em maiores profundidades (Allen e
Allen, 2005). Uma vez que a velocidade compressional relaciona-se com o inverso da funcao
compactacao, e possıvel e factıvel representar espessos pacotes de folhelhos, argilitos ou
siltitos argilosos com funcoes logarıtmicas do tipo.
v(x, z) = v0 +B1ln(z), (91)
onde v0 e B1 sao coeficientes que dependerao da litologia em estudo. Nao se constatou o
emprego destas funcoes (exponencial e logarıtmica) para representacao do modelo de veloci-
dade em tomografia nas referencias consultadas. Todavia, a Figura 1 que ilustra a variacao
de velocidade de um arenito sob diferentes pressoes confinantes se ajusta perfeitamente a
uma funcao logarıtmica.
Discretizacao por funcao de base radial gaussiana
As funcoes de base radial gaussiana (Buhmann, 2003) i.e.
v(x, z) =n
∑
i=1
Aie− 1
2
[
(x−xi)2
σ2xi
+(z−zi)
2
σ2zi
]
, (92)
onde Ai sao parametros relacionados ao peso ou apice da gaussiana, σxie σzi
relacionam-se
a abertura da superfıcie nas direcoes x e z e, finalmente xi e zi constituem as coordenadas
do centro da gaussiana. Se os parametros σxie σzi
forem constantes para ambas as direcoes
e para qualquer i, os unicos coeficientes necessarios para descrever um modelo resumem-se
a Ai.
53
Discretizacao por funcoes harmonicas
De forma semelhante ao topico anterior, o modelo m pode ser representado por funcoes
harmonicas como a serie de Fourier:
v(x, z) =m
∑
jz=0
n∑
jx=0
Ajx,jzei2π(kjxx+kjzz), (93)
onde Ajx,jz sao coeficientes, i =√−1, kjx e kjz sao as frequencias espaciais respectivamente
nas direcoes x e z.
A soma de cossenos e senos (equacao 94) sem a parcela imaginaria da equacao 93,
tambem pode exibir bons resultados para representacao de modelos contınuos, a depender
do numero de parametros empregados.
v(x, z) =n
∑
i=0
Ai [cos(kx,ix) + sen(kx,ix)] +m
∑
j=0
Aj [cos(kz,jz) + sen(kz,jz)], (94)
onde i e j sao ındices variando de 0 a n e 0 a m respectivamente, Ai o coeficiente relativo a
direcao x, Aj o coeficiente na direcao z, kx,i e o numero de onda na direcao x, kz,j o numero
de onda na direcao z.
Santos e Figueiro (2006) descrevem modelos de velocidade atraves de polinomios trigonometricos,
um somatorio de funcoes harmonicas com coeficientes calculados de forma propria e sem
envolvimento de parcelas imaginarias. Os autores conseguem representar feicoes complexas
reduzindo em mais de 6 vezes (de 494 para 80) o numero de parametros necessarios para
discretizar o mesmo modelo atraves de celulas.
54
Discretizacao por splines cubicas (1D e 2D)
Observe-se primeiramente o caso 1D. Seja uma funcao C unidimensional, discreta,
definida no domınio x = [xa, xb] e representada por v=[v1, v2, v3, v4, v5, .... vn] = [vj], a
qual se deseja aproximar por funcoes de interpolacao cubica fj(x) (spline cubica) em cada
subdomınio [xj , xj+1] e com afastamento hx entre os nos do domınio. Sendo v o campo
de velocidade amostrado pontualmente, a funcao velocidade em cada subdomınio [xj , xj+1]
sera dada por:
fj(x) = Aj0 +Aj1(x− xj) +Aj2(x− xj)2 +Aj3(x− xj)
3. (95)
De acordo com as caracterısticas de spline em continuidade da funcao e suas derivadas
primeira e segunda (Kreyszig, 1999), advem:
Aj0 = fj(xj) = vj , (96)
Aj1 = f′
j(xj) = bj, (97)
Aj2 = f′′
j (xj) =3
h2x
(vj+1 − vj) −1
hx(bj+1 + 2bj), (98)
Aj3 = f′′′
j (xj) =2
h3x
(vj − vj+1) −1
h2x
(bj+1 + bj), (99)
onde f′
j, f′′
j e f′′′
j correspondem respectivamente as derivadas primeira, segunda e terceira
de fj em relacao a x. Os parametros bj sao calculados a partir da resolucao do sistema
resumido a seguir:
bj−1 + 4bj + bj+1 =3
hx(vj+1 − vj−1), (100)
ou, de forma expandida:
A extensao para o caso 2D e alcancada por meio de uma nova funcao obtida com o
produto de duas funcoes 1D, uma na direcao x e outra na direcao z concorrendo em um
55
mesmo no (x, z). Assim, para um domınio definido por x = [xa, xb] e z = [za, zb] com nos
em x e z definidos em v = [vi,j], o calculo do valor da funcao em cada domınio e realizado
por:
fi,j(x, z) = Ai0,j0 +Ai1,j1(x−xj)(z− zi)+Ai2,j2(x−xj)2(z− zi)
2 +Ai3,j3(x−xj)3(z− zi)3,
(101)
onde i e j variam de zero a n-1 e m-1 respectivamente, os coeficientes bj e Ai,j sao obtidos
analogamente ao caso 1D.
Discretizacao camada a camada
Este tipo de representacao demanda a existencia de camadas individualizaveis e o con-
hecimento dos limites destas. Hegge e Fokkema (1996) empregam a inversao tomografica
global sobre operadores focais em um modelo com 5 camadas homogeneas.
Qualquer uma das funcoes apresentadas anteriormente pode ser empregada isoladamente
ou em conjunto para compor o modelo de velocidade de uma camada. Kabir (1997) aplicou
a inversao tomografica sobre operadores focais e representou cada camada utilizando funcoes
lineares capazes de comportar variacoes laterais de velocidade em cada estrato. A funcao
utilizada por Kabir (1997) e definida como
vi(x, z) = v0,i + [zc,i(x) − z0]Bi. (102)
Esta funcao emprega dois parametros por camada: v0,i, uma velocidade media e o gradiente
Bi. Na equacao 102 zc,i(x) constitui, para cada x, a profundidade do centro da camada i;
z0 a profundidade de referencia, podendo ser o nıvel do mar.
56
8- A OBTENCAO DA MATRIZ DE SENSIBILIDADE (L)
57
9- INVERSAO TOMOGRAFICA COM FIDELIDADE GEOLOGICA
9.1- Estrategias para estimativa de velocidades geologicamente crıveis
- Escolha de um modelo inicial adequado e proximo da solucao
- Regularizacao
- Discretizacao adequada (analise da matriz de sensibilidade)
58
Vel
oci
ty (
m/s
)
b)a)
IV
V
Sea floor Sea floor
HZ A
HZ B
BSR
HZ A
HZ B
BSR
SW NE SW NE
994995 997
⊗⊗ ⊗
I
II
III
IV
I
II
III
V
IV
V
Figure 11: Exemplo de inversao em Blake Ridge (Santos et al., 2012).
No exemplo a seguir o modelo foi discretizado com splines 2D, o modelo inicial constituiu
um v(f(x)) onde f(x) e a topografia do fundo do mar e foram impostos vınculos na matriz
de sensibilidade. Vınculos que, na pratica, reduziram a quase zero, o efeito de atualizacao
do modelo na camada de agua.
59
No exemplo a seguir foi empregada a atitude das camadas para inserir vınculos que
privilegiassem os gradientes de velocidade com correspondencia nos dados sısmicos.
60
Figure 12: Exemplo de inversao tomografica guiada por estruturas (Chen et al., 2014).
61
Figure 13: Exemplo de inversao com integracao geologica (Chen et al., 2017).Imagem e
velocidade sem os corpos de folhelho (a) comparada com corpos geobodies de folhelho
inseridos no modelo (b).
Neste exemplo, tomografia e outros metodos, inclusive FWI, foram empregados de forma
integrada para obtencao de um campo de velocidade veraz geologicamente.
62
Este exemplo de inversao inseriu uma serie de vınculos que privilegiaram as estruturas
geologicas.
63
Figure 14: Exemplo de inversao tomografica levando em conta o mergulho das camadas
(Costa et al., 2008).
64
Na sequencia das Figuras 15 a 18 temos um exemplo sintetico onde se tem o modelo
sobre o qual foram obtidos os tempos de transito por modelagem 15. Esses dados foram
submetidos a diferentes metodos de inversao. Inversao tomografica camada a camada em
que cada uma e descrita por spline 1D. O resultado desta inversao esta na Figura 16.
A Figura 17 e resultado da inversao camada a camada na qual cada uma foi representada
por uma funcao v(f(x)), onde f(x) e a geometria do horizonte de topo da camada.
Uma vez conhecidas as geometrias de cada horizonte, a partir da inversao que gerou a
Figura 16, foi relizadas a inversao global representando o modelo de velocidade por spline
2D. O resultado encontra-se na Figura 18.
65
Figure 15: Modelo original (Santos, 2012).
66
Figure 16: Resultado de inversao tomografica camada a camada com discrertizacao por
spline 1D (Santos, 2012).
67
Figure 17: Resultado de inversao tomografica camada a camada empregando discretizacao
LASMOD topo (Santos, 2012).
68
Figure 18: Exemplo de inversao tomografica global com discretizacao de todo o modelo por
spline 2D (Santos, 2012).
69
10- ACESSO AOS EXERCICIOS PRATICOS VIA WEB
EXERCICIO 1
Definicoes e conceitos
a) Autovalor e autovetor - exiba as formulas para seus calculos.
b) Calcule os autovalores e autovetores das matrizes a seguir:
b.1)
4 1
6 3
b.2)
1 0 0
1 2 0
2 3 3
..
c) Vetores ortonormais, o que sao (formula)?
d) Para uma dada matriz quadrada, defina numero de condicao.
e) Quais os 3 principais quesitos que devem ser atendidos simultaneamente para que um
problema inverso seja bem posto ?
f) O que e regularizacao ?
g) Aponte as diferencas entre problema inverso linear e nao linear.
70
EXERCICIO 2
Para um problema de inversao tomografica:
a) Escreva um algoritmo para inversao linear.
b) Escreva um algoritmo para inversao nao linear.
c) Aponte as diferencas entre as inversoes linear e nao-linear.
EXERCICIO 3
Com base no sismograma da Figura 19, eixo vertical em tempo (segundos) e eixo horizontal
representando o afastamento fonte receptor (em metros):
a) Estime a velocidade NMO para cada um dos eventos hiperbolicos.
b) Calcule as velocidades intervalares.
c) Desenhe esquematicamente o modelo em profundidade. Justifique.
71
Figure 19: CMP com onda direta e 2 eventos hiperbolicos.
72
EXERCICIO 4
v1 v2
v3 v4
Table 2: Modelo simples com quatro celulas com dimensoes s=100m
O script script gmi1 (em alusao as iniciais de Gera (modelo), Modela e Inverte, primeira
versao) tem por objetivo criar um modelo simples de velocidade, modelar os tempos de
transito e, finalmente, estimar o campo de velocidade por inversao. Trata-se de um script
auto-explicativo, no qual sao empregadas duas rotinas (gcv1 e produtomatrizvetor - de au-
toria de Luiz Santos), composto pelas quatro partes descritas abaixo.
A) A primeira parte, intitula-se (”1-CRIANDO O MODELO”). Nela sao definidas 4
variaveis, v1, v2, v3 e v4, correspondentes as 4 celulas do modelo (vide script). Note
que entre as linhas 9 a 12 do script sao definidas as vagarosidades que sao simplesmente o
inverso da velocidade definida entre as linhas 4 e 7. Na linha 25, e definida a variavel s que
encerra a dimensao das celulas que sao quadradas.
B) Na segunda parte (”2- GERANDO A MATRIZ L (TRANSPOSTA)”) devemos gerar a
matriz L, ja transposta. Ela depende da geometria de aquisicao. Como temos 4 incognitas,
sao necessarios 4 raios. Os raios lancados foram os seguintes:
raio 1: atravessa as celulas 1 e 2 horizontalmente fornecendo o tempo T1;
raio 2: atravessa a celula 2 horizontalmente fornecendo o tempo T2;
raio 3: atravessa as celulas 3 e 4 horizontalmente fornecendo o tempo T3;
raio 4: atravessa a celula e 4 horizontalmente fornecendo o tempo T4;
73
C) Na terceira parte (”3- CALCULANDO TEMPOS DE TRANSITO”) o tempo de transito
e calculado a partir do produto entre a matriz L, gerada no passo anterior, e o vetor de
vagarosidades gerado no passo A.
D) Finalmente, na quarta parte (”4- INVERTENDO - RECUPERANDO MODELO DE
VELOCIDADE”) realiza-se o processo inverso. A partir da matriz L, e dos tempos de
transito, estima-se o campo de velocidade. Monta-se o sistema L.T e obtem-se a solucao
que e o campo de velocidade.
Ao executar o script, digitando no terminal ./script gmi1, todos os 4 passos sao realiza-
dos. O modelo gerado no passo 1 foi:
v1=1000.0
v2=2000.0
v3=1500.0
v4=2500.0
e a inversao estimou:
v1=9.9999993896e+02
v2=2.0000001221e+03
v3=1.4999246826e+03
v4=2.5000004883e+03
Este bom resultado atesta que o dispositivo de aquisicao foi adequado para recuperar o
modelo de velocidade a partir da inversao.
Perguntas:
1) Avalie as seguintes situacoes para o mesmo modelo de velocidade:
74
a) Dispositivo de aquisicao com dois raios:
raio 1: atravessa as celulas 1 e 2 horizontalmente e fornecendo o tempo T1;
raio 2: atravessa as celulas 3 e 4 horizontalmente e fornecendo o tempo T2;
Que resultados voce obteve ? Por que ?
b) Dispositivo de aquisicao com quatro raios:
raio 1: atravessa as celulas 1 e 2 horizontalmente e fornecendo o tempo T1;
raio 2: atravessa as celulas 1 e 2 obliquamente com angulo de 32,8 graus com a horizontal
e fornecendo o tempo T2;
raio 3: atravessa as celulas 3 e 4 horizontalmente e fornecendo o tempo T3;
raio 4: atravessa as celulas 3 e 4 obliquamente com angulo de 32,8 graus com a horizontal
e fornecendo o tempo T4;
Que resultados voce obteve ? Por que ?
Comente.
c) Dispositivo de aquisicao com quatro raios:
raio 1: atravessa as celulas 1 e 2 horizontalmente e fornecendo o tempo T1;
raio 2: atravessa as celulas 1 e 4 obliquamente com angulo de 45 graus com a horizontal,
partindo do canto superior esquerdo da celula 1, e fornecendo o tempo T2;
raio 3: atravessa as celulas 3 e 4 horizontalmente e fornecendo o tempo T3;
raio 4: atravessa as celulas 3 e 2 obliquamente com angulo de 45 graus com a horizontal,
partindo do canto inferior esquerdo da celula 3, e fornecendo o tempo T4;
Que resultados voce obteve ? Por que ?
Comente.
d) Em sua opiniao, qual seria o melhor dispositivo de aquisicao para o modelo proposto?
75
EXERCICIO 5: Inversao linear
Com o auxılio das rotinas em inversao linear4.tar (explicada em sala):
a) Avalie a influencia do espacamento entre fontes e entre receptores na solucao do problema
inverso. Inclua em sua analise observacoes sobre a matriz AtA e seus respectivos numeros de
condicao para cada caso. O calculo do numero de condicao pode ser alcancado no diretorio
de analise, com o script svd.
EXERCICIO 6
A ausencia de solucao e observada pela baixa densidade de raios no exemplo anterior. Com
o auxılio do script lusolver avalie e discuta a regularizacao de ordem zero.
EXERCICIO 7: Inversao nao linear para eventos com profundidade conhecida
Com o auxılio das rotinas em inversao n linear5.tar (explicada em sala) estime o modelo a
partir dos tempos de transito. Avalie o resultado da inversao testando:
a) o modelo inicial
b) o tamanho do passo
c) a convergencia
d) a tolerancia
EXERCICIO 8: Inverso nao linear para eventos com profundidade desconhecida
Com o auxılio das rotinas em inversao n linear6.tar (explicada em sala) estime o modelo a
partir dos tempos de transito. Avalie o resultado da inversao testando:
a) o modelo inicial
76
b) o tamanho do passo
c) a convergencia
d) a tolerancia
77
11- REFERENCIAS
78
